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

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

Parent Directory Parent Directory | Revision Log Revision Log

Revision 23 - (hide 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 guez 23 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