/[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 158 - (show annotations)
Tue Jul 21 14:44:45 2015 UTC (8 years, 9 months ago) by guez
File size: 7443 byte(s)
Subroutine sugwd sets variables of module yoegwd. Better to put it
into module yoegwd.

Variables of module yoegwd other than NKTOPG, NSTRA can be symbolic
constants. sugwd now only sets NKTOPG, NSTRA. Simplified the
computation of NKTOPG, NSTRA by making the local variable zpm1r an
array instead of a scalar and calling ifirstloc.

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 !* 3. COMPUTE MOUNTAIN LIFT
160
161 DO 301 jl = 1, klon
162 IF (ktest(jl)==1) THEN
163 ztau(jl,klev+1) = -gklift*zrho(jl,klev+1)*2.*romega*2*pvaror(jl)*sin &
164 (zpi/180.*plat(jl))*pvlow(jl)
165 ztav(jl,klev+1) = gklift*zrho(jl,klev+1)*2.*romega*2*pvaror(jl)* &
166 sin(zpi/180.*plat(jl))*pulow(jl)
167 ELSE
168 ztau(jl,klev+1) = 0.0
169 ztav(jl,klev+1) = 0.0
170 END IF
171 301 CONTINUE
172
173
174 !* 4. COMPUTE LIFT PROFILE
175 !* --------------------
176
177
178 DO jk = 1, klev
179 DO jl = 1, klon
180 IF (ktest(jl)==1) THEN
181 ztau(jl,jk) = ztau(jl,klev+1)*paphm1(jl,jk)/paphm1(jl,klev+1)
182 ztav(jl,jk) = ztav(jl,klev+1)*paphm1(jl,jk)/paphm1(jl,klev+1)
183 ELSE
184 ztau(jl,jk) = 0.0
185 ztav(jl,jk) = 0.0
186 END IF
187 end DO
188 end DO
189
190
191 !* 5. COMPUTE TENDENCIES.
192 !* -------------------
193 IF (lifthigh) THEN
194
195 ! PRINT *,' DANS OROLIFT: 500'
196
197 ! EXPLICIT SOLUTION AT ALL LEVELS
198
199 DO 524 jk = 1, klev
200 DO 523 jl = 1, klon
201 IF (ktest(jl)==1) THEN
202 zdelp = paphm1(jl,jk+1) - paphm1(jl,jk)
203 zdudt(jl) = -rg*(ztau(jl,jk+1)-ztau(jl,jk))/zdelp
204 zdvdt(jl) = -rg*(ztav(jl,jk+1)-ztav(jl,jk))/zdelp
205 END IF
206 523 CONTINUE
207 524 CONTINUE
208
209 ! PROJECT PERPENDICULARLY TO U NOT TO DESTROY ENERGY
210
211 DO jk = 1, klev
212 DO jl = 1, klon
213 IF (ktest(jl)==1) THEN
214
215 zslow = sqrt(pulow(jl)**2+pvlow(jl)**2)
216 zsqua = amax1(sqrt(pum1(jl,jk)**2+pvm1(jl,jk)**2),gvsec)
217 zscav = -zdudt(jl)*pvm1(jl,jk) + zdvdt(jl)*pum1(jl,jk)
218 IF (zsqua>gvsec) THEN
219 pvom(jl,jk) = -zscav*pvm1(jl,jk)/zsqua**2
220 pvol(jl,jk) = zscav*pum1(jl,jk)/zsqua**2
221 ELSE
222 pvom(jl,jk) = 0.0
223 pvol(jl,jk) = 0.0
224 END IF
225 zsqua = sqrt(pum1(jl,jk)**2+pum1(jl,jk)**2)
226 IF (zsqua<zslow) THEN
227 pvom(jl,jk) = zsqua/zslow*pvom(jl,jk)
228 pvol(jl,jk) = zsqua/zslow*pvol(jl,jk)
229 END IF
230
231 END IF
232 end DO
233 end DO
234
235 ! 6. LOW LEVEL LIFT, SEMI IMPLICIT:
236 ! ----------------------------------
237
238 ELSE
239
240 DO 601 jl = 1, klon
241 IF (ktest(jl)==1) THEN
242 DO jk = klev, iknub(jl), -1
243 zbet = gklift*2.*romega*sin(zpi/180.*plat(jl))*ztmst* &
244 (pgeom1(jl,iknub(jl)-1)-pgeom1(jl,jk))/ &
245 (pgeom1(jl,iknub(jl)-1)-pgeom1(jl,klev))
246 zdudt(jl) = -pum1(jl,jk)/ztmst/(1+zbet**2)
247 zdvdt(jl) = -pvm1(jl,jk)/ztmst/(1+zbet**2)
248 pvom(jl,jk) = zbet**2*zdudt(jl) - zbet*zdvdt(jl)
249 pvol(jl,jk) = zbet*zdudt(jl) + zbet**2*zdvdt(jl)
250 END DO
251 END IF
252 601 CONTINUE
253
254 END IF
255
256 RETURN
257 END

  ViewVC Help
Powered by ViewVC 1.1.21