/[lmdze]/trunk/libf/phylmd/Orography/orosetup.f90
ViewVC logotype

Annotation of /trunk/libf/phylmd/Orography/orosetup.f90

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
File size: 11419 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 orosetup(nlon, ktest, kkcrit, kkcrith, kcrit, kkenvh, kknu, kknu2, &
2     paphm1, papm1, pum1, pvm1, ptm1, pgeom1, pstd, prho, pri, pstab, ptau, &
3     pvph, ppsi, pzdep, pulow, pvlow, ptheta, pgamma, pmea, ppic, pval, pnu, &
4     pd1, pd2, pdmod)
5    
6     ! *gwsetup*
7     ! interface from *orodrag*
8     ! see ecmwf research department documentation of the "i.f.s."
9     ! modifications f.lott for the new-gwdrag scheme november 1993
10    
11     USE dimens_m
12     USE dimphy
13     USE yomcst
14     USE yoegwd
15    
16     IMPLICIT NONE
17    
18     ! 0.1 arguments
19    
20     INTEGER nlon
21     INTEGER jl, jk
22     REAL zdelp
23    
24     INTEGER kkcrit(nlon), kkcrith(nlon), kcrit(nlon), ktest(nlon), &
25     kkenvh(nlon)
26    
27     REAL paphm1(nlon, klev+1), papm1(nlon, klev), pum1(nlon, klev), &
28     pvm1(nlon, klev), ptm1(nlon, klev), pgeom1(nlon, klev), &
29     prho(nlon, klev+1), pri(nlon, klev+1), pstab(nlon, klev+1), &
30     ptau(nlon, klev+1), pvph(nlon, klev+1), ppsi(nlon, klev+1), &
31     pzdep(nlon, klev)
32     REAL pulow(nlon), pvlow(nlon), ptheta(nlon), pgamma(nlon), pnu(nlon), &
33     pd1(nlon), pd2(nlon), pdmod(nlon)
34     REAL, INTENT (IN) :: pstd(nlon)
35     REAL pmea(nlon), ppic(nlon), pval(nlon)
36    
37     ! 0.2 local arrays
38    
39     INTEGER ilevm1, ilevm2, ilevh
40     REAL zcons1, zcons2, zcons3, zhgeo
41     REAL zu, zphi, zvt1, zvt2, zst, zvar, zdwind, zwind
42     REAL zstabm, zstabp, zrhom, zrhop
43     REAL zggeenv, zggeom1, zgvar
44     LOGICAL lo
45     LOGICAL ll1(klon, klev+1)
46     INTEGER kknu(klon), kknu2(klon), kknub(klon), kknul(klon), kentp(klon), &
47     ncount(klon)
48    
49     REAL zhcrit(klon, klev), zvpf(klon, klev), zdp(klon, klev)
50     REAL znorm(klon), zb(klon), zc(klon), zulow(klon), zvlow(klon), &
51     znup(klon), znum(klon)
52    
53     !------------------------------------------------------------------
54    
55     !!print *, "Call sequence information: orosetup"
56     ! 1. initialization
57     ! 1.1 computational constants
58    
59     ilevm1 = klev - 1
60     ilevm2 = klev - 2
61     ilevh = klev/3
62    
63     zcons1 = 1./rd
64     !old zcons2=g**2/cpd
65     zcons2 = rg**2/rcpd
66     !old zcons3=1.5*api
67     zcons3 = 1.5*rpi
68    
69     ! 2.
70    
71     ! 2.1 define low level wind, project winds in plane of
72     ! low level wind, determine sector in which to take
73     ! the variance and set indicator for critical levels.
74    
75     DO jl = 1, klon
76     kknu(jl) = klev
77     kknu2(jl) = klev
78     kknub(jl) = klev
79     kknul(jl) = klev
80     pgamma(jl) = max(pgamma(jl), gtsec)
81     ll1(jl, klev+1) = .FALSE.
82     end DO
83    
84     ! Ajouter une initialisation (L. Li, le 23fev99):
85    
86     DO jk = klev, ilevh, -1
87     DO jl = 1, klon
88     ll1(jl, jk) = .FALSE.
89     END DO
90     END DO
91    
92     ! define top of low level flow
93    
94     DO jk = klev, ilevh, -1
95     DO jl = 1, klon
96     lo = (paphm1(jl, jk)/paphm1(jl, klev+1)) >= gsigcr
97     IF (lo) THEN
98     kkcrit(jl) = jk
99     END IF
100     zhcrit(jl, jk) = ppic(jl)
101     zhgeo = pgeom1(jl, jk)/rg
102     ll1(jl, jk) = (zhgeo>zhcrit(jl, jk))
103     IF (ll1(jl, jk) .NEQV. ll1(jl, jk+1)) THEN
104     kknu(jl) = jk
105     END IF
106     IF ( .NOT. ll1(jl, ilevh)) kknu(jl) = ilevh
107     end DO
108     end DO
109     DO jk = klev, ilevh, -1
110     DO jl = 1, klon
111     zhcrit(jl, jk) = ppic(jl) - pval(jl)
112     zhgeo = pgeom1(jl, jk)/rg
113     ll1(jl, jk) = (zhgeo>zhcrit(jl, jk))
114     IF (ll1(jl, jk) .NEQV. ll1(jl, jk+1)) THEN
115     kknu2(jl) = jk
116     END IF
117     IF ( .NOT. ll1(jl, ilevh)) kknu2(jl) = ilevh
118     end DO
119     end DO
120     DO jk = klev, ilevh, -1
121     DO jl = 1, klon
122     zhcrit(jl, jk) = amax1(ppic(jl)-pmea(jl), pmea(jl)-pval(jl))
123     zhgeo = pgeom1(jl, jk)/rg
124     ll1(jl, jk) = (zhgeo>zhcrit(jl, jk))
125     IF (ll1(jl, jk) .NEQV. ll1(jl, jk+1)) THEN
126     kknub(jl) = jk
127     END IF
128     IF ( .NOT. ll1(jl, ilevh)) kknub(jl) = ilevh
129     end DO
130     end DO
131    
132     DO jl = 1, klon
133     kknu(jl) = min(kknu(jl), nktopg)
134     kknu2(jl) = min(kknu2(jl), nktopg)
135     kknub(jl) = min(kknub(jl), nktopg)
136     kknul(jl) = klev
137     end DO
138    
139     !c* initialize various arrays
140    
141     DO jl = 1, klon
142     prho(jl, klev+1) = 0.0
143     pstab(jl, klev+1) = 0.0
144     pstab(jl, 1) = 0.0
145     pri(jl, klev+1) = 9999.0
146     ppsi(jl, klev+1) = 0.0
147     pri(jl, 1) = 0.0
148     pvph(jl, 1) = 0.0
149     pulow(jl) = 0.0
150     pvlow(jl) = 0.0
151     zulow(jl) = 0.0
152     zvlow(jl) = 0.0
153     kkcrith(jl) = klev
154     kkenvh(jl) = klev
155     kentp(jl) = klev
156     kcrit(jl) = 1
157     ncount(jl) = 0
158     ll1(jl, klev+1) = .FALSE.
159     end DO
160    
161     ! define low-level flow
162    
163     DO jk = klev, 2, -1
164     DO jl = 1, klon
165     IF (ktest(jl)==1) THEN
166     zdp(jl, jk) = papm1(jl, jk) - papm1(jl, jk-1)
167     prho(jl, jk) = 2. * paphm1(jl, jk) * zcons1 &
168     / (ptm1(jl, jk) + ptm1(jl, jk-1))
169     pstab(jl, jk) = 2. * zcons2 / (ptm1(jl, jk) + ptm1(jl, jk-1)) &
170     * (1. - rcpd * prho(jl, jk) &
171     * (ptm1(jl, jk) - ptm1(jl, jk - 1)) / zdp(jl, jk))
172     pstab(jl, jk) = max(pstab(jl, jk), gssec)
173     END IF
174     end DO
175     end DO
176    
177     ! define blocked flow
178    
179     DO jk = klev, ilevh, -1
180     DO jl = 1, klon
181     IF (jk>=kknub(jl) .AND. jk<=kknul(jl)) THEN
182     pulow(jl) = pulow(jl) + pum1(jl, jk)*(paphm1(jl, jk+1)-paphm1(jl, jk) &
183     )
184     pvlow(jl) = pvlow(jl) + pvm1(jl, jk)*(paphm1(jl, jk+1)-paphm1(jl, jk) &
185     )
186     END IF
187     end DO
188     end DO
189     DO jl = 1, klon
190     pulow(jl) = pulow(jl)/(paphm1(jl, kknul(jl)+1)-paphm1(jl, kknub(jl)))
191     pvlow(jl) = pvlow(jl)/(paphm1(jl, kknul(jl)+1)-paphm1(jl, kknub(jl)))
192     znorm(jl) = max(sqrt(pulow(jl)**2+pvlow(jl)**2), gvsec)
193     pvph(jl, klev+1) = znorm(jl)
194     end DO
195    
196     ! setup orography axes and define plane of profiles
197    
198     DO jl = 1, klon
199     lo = (pulow(jl)<gvsec) .AND. (pulow(jl)>=-gvsec)
200     IF (lo) THEN
201     zu = pulow(jl) + 2.*gvsec
202     ELSE
203     zu = pulow(jl)
204     END IF
205     zphi = atan(pvlow(jl)/zu)
206     ppsi(jl, klev+1) = ptheta(jl)*rpi/180. - zphi
207     zb(jl) = 1. - 0.18*pgamma(jl) - 0.04*pgamma(jl)**2
208     zc(jl) = 0.48*pgamma(jl) + 0.3*pgamma(jl)**2
209     pd1(jl) = zb(jl) - (zb(jl)-zc(jl))*(sin(ppsi(jl, klev+1))**2)
210     pd2(jl) = (zb(jl)-zc(jl))*sin(ppsi(jl, klev+1))*cos(ppsi(jl, klev+1))
211     pdmod(jl) = sqrt(pd1(jl)**2+pd2(jl)**2)
212     end DO
213    
214     ! define flow in plane of lowlevel stress
215    
216     DO jk = 1, klev
217     DO jl = 1, klon
218     IF (ktest(jl)==1) THEN
219     zvt1 = pulow(jl)*pum1(jl, jk) + pvlow(jl)*pvm1(jl, jk)
220     zvt2 = -pvlow(jl)*pum1(jl, jk) + pulow(jl)*pvm1(jl, jk)
221     zvpf(jl, jk) = (zvt1*pd1(jl)+zvt2*pd2(jl))/(znorm(jl)*pdmod(jl))
222     END IF
223     ptau(jl, jk) = 0.0
224     pzdep(jl, jk) = 0.0
225     ppsi(jl, jk) = 0.0
226     ll1(jl, jk) = .FALSE.
227     end DO
228     end DO
229     DO jk = 2, klev
230     DO jl = 1, klon
231     IF (ktest(jl)==1) THEN
232     zdp(jl, jk) = papm1(jl, jk) - papm1(jl, jk-1)
233     pvph(jl, jk) = ((paphm1(jl, jk)-papm1(jl, jk-1))*zvpf(jl, jk)+(papm1( &
234     jl, jk)-paphm1(jl, jk))*zvpf(jl, jk-1))/zdp(jl, jk)
235     IF (pvph(jl, jk)<gvsec) THEN
236     pvph(jl, jk) = gvsec
237     kcrit(jl) = jk
238     END IF
239     END IF
240     end DO
241     end DO
242    
243     ! 2.2 brunt-vaisala frequency and density at half levels.
244    
245     DO jk = ilevh, klev
246     DO jl = 1, klon
247     IF (ktest(jl)==1) THEN
248     IF (jk>=(kknub(jl)+1) .AND. jk<=kknul(jl)) THEN
249     zst = zcons2/ptm1(jl, jk)*(1.-rcpd*prho(jl, jk)*(ptm1(jl, &
250     jk)-ptm1(jl, jk-1))/zdp(jl, jk))
251     pstab(jl, klev+1) = pstab(jl, klev+1) + zst*zdp(jl, jk)
252     pstab(jl, klev+1) = max(pstab(jl, klev+1), gssec)
253     prho(jl, klev+1) = prho(jl, klev+1) + paphm1(jl, jk)*2.*zdp(jl, jk)* &
254     zcons1/(ptm1(jl, jk)+ptm1(jl, jk-1))
255     END IF
256     END IF
257     end DO
258     end DO
259    
260     DO jl = 1, klon
261     pstab(jl, klev + 1) = pstab(jl, klev + 1) &
262     / (papm1(jl, kknul(jl)) - papm1(jl, kknub(jl)))
263     prho(jl, klev + 1) = prho(jl, klev + 1) &
264     / (papm1(jl, kknul(jl)) - papm1(jl, kknub(jl)))
265     zvar = pstd(jl)
266     end DO
267    
268     ! 2.3 mean flow richardson number.
269     ! and critical height for froude layer
270    
271     DO jk = 2, klev
272     DO jl = 1, klon
273     IF (ktest(jl)==1) THEN
274     zdwind = max(abs(zvpf(jl, jk)-zvpf(jl, jk-1)), gvsec)
275     pri(jl, jk) = pstab(jl, jk)*(zdp(jl, jk)/(rg*prho(jl, jk)*zdwind))**2
276     pri(jl, jk) = max(pri(jl, jk), grcrit)
277     END IF
278     end DO
279     end do
280    
281     ! define top of 'envelope' layer
282    
283     DO jl = 1, klon
284     pnu(jl) = 0.0
285     znum(jl) = 0.0
286     end DO
287    
288     DO jk = 2, klev - 1
289     DO jl = 1, klon
290    
291     IF (ktest(jl)==1) THEN
292    
293     IF (jk>=kknub(jl)) THEN
294    
295     znum(jl) = pnu(jl)
296     zwind = (pulow(jl)*pum1(jl, jk)+pvlow(jl)*pvm1(jl, jk))/ &
297     max(sqrt(pulow(jl)**2+pvlow(jl)**2), gvsec)
298     zwind = max(sqrt(zwind**2), gvsec)
299     zdelp = paphm1(jl, jk+1) - paphm1(jl, jk)
300     zstabm = sqrt(max(pstab(jl, jk), gssec))
301     zstabp = sqrt(max(pstab(jl, jk+1), gssec))
302     zrhom = prho(jl, jk)
303     zrhop = prho(jl, jk+1)
304     pnu(jl) = pnu(jl) + (zdelp/rg)*((zstabp/zrhop+zstabm/zrhom)/2.)/ &
305     zwind
306     IF ((znum(jl)<=gfrcrit) .AND. (pnu(jl)>gfrcrit) .AND. (kkenvh( &
307     jl)==klev)) kkenvh(jl) = jk
308    
309     END IF
310    
311     END IF
312    
313     end DO
314     end do
315    
316     ! calculation of a dynamical mixing height for the breaking
317     ! of gravity waves:
318    
319     DO jl = 1, klon
320     znup(jl) = 0.0
321     znum(jl) = 0.0
322     end DO
323    
324     DO jk = klev - 1, 2, -1
325     DO jl = 1, klon
326    
327     IF (ktest(jl)==1) THEN
328    
329     znum(jl) = znup(jl)
330     zwind = (pulow(jl)*pum1(jl, jk)+pvlow(jl)*pvm1(jl, jk))/ &
331     max(sqrt(pulow(jl)**2+pvlow(jl)**2), gvsec)
332     zwind = max(sqrt(zwind**2), gvsec)
333     zdelp = paphm1(jl, jk+1) - paphm1(jl, jk)
334     zstabm = sqrt(max(pstab(jl, jk), gssec))
335     zstabp = sqrt(max(pstab(jl, jk+1), gssec))
336     zrhom = prho(jl, jk)
337     zrhop = prho(jl, jk+1)
338     znup(jl) = znup(jl) + (zdelp/rg)*((zstabp/zrhop+zstabm/zrhom)/2.)/ &
339     zwind
340     IF ((znum(jl)<=rpi/2.) .AND. (znup(jl)>rpi/2.) .AND. (kkcrith( &
341     jl)==klev)) kkcrith(jl) = jk
342    
343     END IF
344    
345     end DO
346     end DO
347    
348     DO jl = 1, klon
349     kkcrith(jl) = min0(kkcrith(jl), kknu2(jl))
350     kkcrith(jl) = max0(kkcrith(jl), ilevh*2)
351     end DO
352    
353     ! directional info for flow blocking
354    
355     DO jk = ilevh, klev
356     DO jl = 1, klon
357     IF (jk>=kkenvh(jl)) THEN
358     lo = (pum1(jl, jk)<gvsec) .AND. (pum1(jl, jk)>=-gvsec)
359     IF (lo) THEN
360     zu = pum1(jl, jk) + 2.*gvsec
361     ELSE
362     zu = pum1(jl, jk)
363     END IF
364     zphi = atan(pvm1(jl, jk)/zu)
365     ppsi(jl, jk) = ptheta(jl)*rpi/180. - zphi
366     END IF
367     end DO
368     end DO
369     ! forms the vertical 'leakiness'
370    
371     DO jk = ilevh, klev
372     DO jl = 1, klon
373     IF (jk>=kkenvh(jl)) THEN
374     zggeenv = amax1(1., (pgeom1(jl, kkenvh(jl))+pgeom1(jl, &
375     kkenvh(jl)-1))/2.)
376     zggeom1 = amax1(pgeom1(jl, jk), 1.)
377     zgvar = amax1(pstd(jl)*rg, 1.)
378     !mod pzdep(jl, jk)=sqrt((zggeenv-zggeom1)/(zggeom1+zgvar))
379     pzdep(jl, jk) = (pgeom1(jl, kkenvh(jl)-1)-pgeom1(jl, jk))/ &
380     (pgeom1(jl, kkenvh(jl)-1)-pgeom1(jl, klev))
381     END IF
382     end DO
383     end DO
384    
385     END SUBROUTINE orosetup

  ViewVC Help
Powered by ViewVC 1.1.21