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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 82 - (show annotations)
Wed Mar 5 14:57:53 2014 UTC (10 years, 2 months ago) by guez
File size: 7565 byte(s)
Changed all ".f90" suffixes to ".f".
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 suphec_m
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