/[lmdze]/trunk/Sources/phylmd/Orography/orolift.f
ViewVC logotype

Contents of /trunk/Sources/phylmd/Orography/orolift.f

Parent Directory Parent Directory | Revision Log Revision Log

Revision 23 - (show annotations)
Mon Dec 14 15:25:16 2009 UTC (14 years, 5 months ago) by guez
Original Path: trunk/libf/phylmd/Orography/orolift.f90
File size: 7563 byte(s) 
Split "orografi.f": one file for each procedure. Put the created files
in new directory "Orography".

Removed argument "vcov" of procedure "sortvarc". Removed arguments
"itau" and "time" of procedure "caldyn0". Removed arguments "itau",
"time" and "vcov" of procedure "sortvarc0".

Removed argument "time" of procedure "dynredem1". Removed NetCDF
variable "temps" in files "start.nc" and "restart.nc", because its
value is always 0.

Removed argument "nq" of procedures "iniadvtrac" and "leapfrog". The
number of "tracers read in "traceur.def" must now be equal to "nqmx",
or "nqmx" must equal 4 if there is no file "traceur.def". Replaced
variable "nq" by constant "nqmx" in "leapfrog".

NetCDF variable for ozone field in "coefoz.nc" must now be called
"tro3" instead of "r".

Fixed bug in "zenang".
1 SUBROUTINE orolift(nlon,nlev,ktest,ptsphy,paphm1,pgeom1,ptm1,pum1,pvm1, &
2 plat,pmea,pvaror,ppic &
3 ,pulow,pvlow,pvom,pvol,pte)
4
5
6 !**** *OROLIFT: SIMULATE THE GEOSTROPHIC LIFT.
7
8 ! PURPOSE.
9 ! --------
10
11 !** INTERFACE.
12 ! ----------
13 ! CALLED FROM *lift_noro
14 ! ----------
15
16 ! AUTHOR.
17 ! -------
18 ! F.LOTT LMD 22/11/95
19
20 USE dimens_m
21 USE dimphy
22 USE yomcst
23 USE yoegwd
24 IMPLICIT NONE
25
26
27 !-----------------------------------------------------------------------
28
29 !* 0.1 ARGUMENTS
30 ! ---------
31
32
33 INTEGER nlon, nlev
34 REAL pte(nlon,nlev), pvol(nlon,nlev), pvom(nlon,nlev), pulow(nlon), &
35 pvlow(nlon)
36 REAL pum1(nlon,nlev), pvm1(nlon,nlev), ptm1(nlon,nlev)
37 REAL, INTENT (IN) :: plat(nlon)
38 REAL pmea(nlon)
39 REAL, INTENT (IN) :: pvaror(nlon)
40 REAL ppic(nlon), pgeom1(nlon,nlev), paphm1(nlon,nlev+1)
41
42 INTEGER ktest(nlon)
43 REAL, INTENT (IN) :: ptsphy
44 !-----------------------------------------------------------------------
45
46 !* 0.2 LOCAL ARRAYS
47 ! ------------
48 LOGICAL lifthigh
49 INTEGER klevm1, jl, ilevh, jk
50 REAL zcons1, ztmst, zrtmst, zpi, zhgeo
51 REAL zdelp, zslow, zsqua, zscav, zbet
52 INTEGER iknub(klon), iknul(klon)
53 LOGICAL ll1(klon,klev+1)
54
55 REAL ztau(klon,klev+1), ztav(klon,klev+1), zrho(klon,klev+1)
56 REAL zdudt(klon), zdvdt(klon)
57 REAL zhcrit(klon,klev)
58 !-----------------------------------------------------------------------
59
60 !* 1.1 INITIALIZATIONS
61 ! ---------------
62
63 lifthigh = .FALSE.
64
65 IF (nlon/=klon .OR. nlev/=klev) STOP
66 zcons1 = 1./rd
67 klevm1 = klev - 1
68 ztmst = ptsphy
69 zrtmst = 1./ztmst
70 zpi = acos(-1.)
71
72 DO 1001 jl = 1, klon
73 zrho(jl,klev+1) = 0.0
74 pulow(jl) = 0.0
75 pvlow(jl) = 0.0
76 iknub(jl) = klev
77 iknul(jl) = klev
78 ilevh = klev/3
79 ll1(jl,klev+1) = .FALSE.
80 DO 1000 jk = 1, klev
81 pvom(jl,jk) = 0.0
82 pvol(jl,jk) = 0.0
83 pte(jl,jk) = 0.0
84 1000 CONTINUE
85 1001 CONTINUE
86
87
88 !* 2.1 DEFINE LOW LEVEL WIND, PROJECT WINDS IN PLANE OF
89 !* LOW LEVEL WIND, DETERMINE SECTOR IN WHICH TO TAKE
90 !* THE VARIANCE AND SET INDICATOR FOR CRITICAL LEVELS.
91
92
93
94 DO 2006 jk = klev, 1, -1
95 DO 2007 jl = 1, klon
96 IF (ktest(jl)==1) THEN
97 zhcrit(jl,jk) = amax1(ppic(jl)-pmea(jl),100.)
98 zhgeo = pgeom1(jl,jk)/rg
99 ll1(jl,jk) = (zhgeo>zhcrit(jl,jk))
100 IF (ll1(jl,jk) .NEQV. ll1(jl,jk+1)) THEN
101 iknub(jl) = jk
102 END IF
103 END IF
104 2007 CONTINUE
105 2006 CONTINUE
106
107 DO 2010 jl = 1, klon
108 IF (ktest(jl)==1) THEN
109 iknub(jl) = max(iknub(jl),klev/2)
110 iknul(jl) = max(iknul(jl),2*klev/3)
111 IF (iknub(jl)>nktopg) iknub(jl) = nktopg
112 IF (iknub(jl)==nktopg) iknul(jl) = klev
113 IF (iknub(jl)==iknul(jl)) iknub(jl) = iknul(jl) - 1
114 END IF
115 2010 CONTINUE
116
117 ! do 2011 jl=1,klon
118 ! IF(KTEST(JL).EQ.1) THEN
119 ! print *,' iknul= ',iknul(jl),' iknub=',iknub(jl)
120 ! ENDIF
121 !2011 continue
122
123 ! PRINT *,' DANS OROLIFT: 2010'
124
125 DO 223 jk = klev, 2, -1
126 DO 222 jl = 1, klon
127 zrho(jl,jk) = 2.*paphm1(jl,jk)*zcons1/(ptm1(jl,jk)+ptm1(jl,jk-1))
128 222 CONTINUE
129 223 CONTINUE
130 ! PRINT *,' DANS OROLIFT: 223'
131
132 !********************************************************************
133
134 !* DEFINE LOW LEVEL FLOW
135 ! -------------------
136 DO 2115 jk = klev, 1, -1
137 DO 2116 jl = 1, klon
138 IF (ktest(jl)==1) THEN
139 IF (jk>=iknub(jl) .AND. jk<=iknul(jl)) THEN
140 pulow(jl) = pulow(jl) + pum1(jl,jk)*(paphm1(jl,jk+1)-paphm1(jl, &
141 jk))
142 pvlow(jl) = pvlow(jl) + pvm1(jl,jk)*(paphm1(jl,jk+1)-paphm1(jl, &
143 jk))
144 zrho(jl,klev+1) = zrho(jl,klev+1) + zrho(jl,jk)*(paphm1(jl,jk+1) &
145 -paphm1(jl,jk))
146 END IF
147 END IF
148 2116 CONTINUE
149 2115 CONTINUE
150 DO 2110 jl = 1, klon
151 IF (ktest(jl)==1) THEN
152 pulow(jl) = pulow(jl)/(paphm1(jl,iknul(jl)+1)-paphm1(jl,iknub(jl)))
153 pvlow(jl) = pvlow(jl)/(paphm1(jl,iknul(jl)+1)-paphm1(jl,iknub(jl)))
154 zrho(jl,klev+1) = zrho(jl,klev+1)/(paphm1(jl,iknul(jl)+1)-paphm1(jl, &
155 iknub(jl)))
156 END IF
157 2110 CONTINUE
158
159
160 200 CONTINUE
161
162 !***********************************************************
163
164 !* 3. COMPUTE MOUNTAIN LIFT
165
166 300 CONTINUE
167
168 DO 301 jl = 1, klon
169 IF (ktest(jl)==1) THEN
170 ztau(jl,klev+1) = -gklift*zrho(jl,klev+1)*2.*romega*2*pvaror(jl)*sin &
171 (zpi/180.*plat(jl))*pvlow(jl)
172 ztav(jl,klev+1) = gklift*zrho(jl,klev+1)*2.*romega*2*pvaror(jl)* &
173 sin(zpi/180.*plat(jl))*pulow(jl)
174 ELSE
175 ztau(jl,klev+1) = 0.0
176 ztav(jl,klev+1) = 0.0
177 END IF
178 301 CONTINUE
179
180
181 !* 4. COMPUTE LIFT PROFILE
182 !* --------------------
183
184
185 400 CONTINUE
186
187 DO 401 jk = 1, klev
188 DO 401 jl = 1, klon
189 IF (ktest(jl)==1) THEN
190 ztau(jl,jk) = ztau(jl,klev+1)*paphm1(jl,jk)/paphm1(jl,klev+1)
191 ztav(jl,jk) = ztav(jl,klev+1)*paphm1(jl,jk)/paphm1(jl,klev+1)
192 ELSE
193 ztau(jl,jk) = 0.0
194 ztav(jl,jk) = 0.0
195 END IF
196 401 CONTINUE
197
198
199 !* 5. COMPUTE TENDENCIES.
200 !* -------------------
201 IF (lifthigh) THEN
202
203 500 CONTINUE
204 ! PRINT *,' DANS OROLIFT: 500'
205
206 ! EXPLICIT SOLUTION AT ALL LEVELS
207
208 DO 524 jk = 1, klev
209 DO 523 jl = 1, klon
210 IF (ktest(jl)==1) THEN
211 zdelp = paphm1(jl,jk+1) - paphm1(jl,jk)
212 zdudt(jl) = -rg*(ztau(jl,jk+1)-ztau(jl,jk))/zdelp
213 zdvdt(jl) = -rg*(ztav(jl,jk+1)-ztav(jl,jk))/zdelp
214 END IF
215 523 CONTINUE
216 524 CONTINUE
217
218 ! PROJECT PERPENDICULARLY TO U NOT TO DESTROY ENERGY
219
220 DO 530 jk = 1, klev
221 DO 530 jl = 1, klon
222 IF (ktest(jl)==1) THEN
223
224 zslow = sqrt(pulow(jl)**2+pvlow(jl)**2)
225 zsqua = amax1(sqrt(pum1(jl,jk)**2+pvm1(jl,jk)**2),gvsec)
226 zscav = -zdudt(jl)*pvm1(jl,jk) + zdvdt(jl)*pum1(jl,jk)
227 IF (zsqua>gvsec) THEN
228 pvom(jl,jk) = -zscav*pvm1(jl,jk)/zsqua**2
229 pvol(jl,jk) = zscav*pum1(jl,jk)/zsqua**2
230 ELSE
231 pvom(jl,jk) = 0.0
232 pvol(jl,jk) = 0.0
233 END IF
234 zsqua = sqrt(pum1(jl,jk)**2+pum1(jl,jk)**2)
235 IF (zsqua<zslow) THEN
236 pvom(jl,jk) = zsqua/zslow*pvom(jl,jk)
237 pvol(jl,jk) = zsqua/zslow*pvol(jl,jk)
238 END IF
239
240 END IF
241 530 CONTINUE
242
243 ! 6. LOW LEVEL LIFT, SEMI IMPLICIT:
244 ! ----------------------------------
245
246 ELSE
247
248 DO 601 jl = 1, klon
249 IF (ktest(jl)==1) THEN
250 DO jk = klev, iknub(jl), -1
251 zbet = gklift*2.*romega*sin(zpi/180.*plat(jl))*ztmst* &
252 (pgeom1(jl,iknub(jl)-1)-pgeom1(jl,jk))/ &
253 (pgeom1(jl,iknub(jl)-1)-pgeom1(jl,klev))
254 zdudt(jl) = -pum1(jl,jk)/ztmst/(1+zbet**2)
255 zdvdt(jl) = -pvm1(jl,jk)/ztmst/(1+zbet**2)
256 pvom(jl,jk) = zbet**2*zdudt(jl) - zbet*zdvdt(jl)
257 pvol(jl,jk) = zbet*zdudt(jl) + zbet**2*zdvdt(jl)
258 END DO
259 END IF
260 601 CONTINUE
261
262 END IF
263
264 RETURN
265 END
  ViewVC Help
Powered by ViewVC 1.1.21