/[lmdze]/trunk/libf/phylmd/orografi.f
ViewVC logotype

Annotation of /trunk/libf/phylmd/orografi.f

Parent Directory Parent Directory | Revision Log Revision Log

Revision 10 - (hide annotations)
Fri Apr 18 14:45:53 2008 UTC (16 years, 1 month ago) by guez
File size: 48760 byte(s) 
Added NetCDF directory "/home/guez/include" in "g95.mk" and
"nag_tools.mk".

Added some "intent" attributes in "PVtheta", "advtrac", "caladvtrac",
"calfis", "diagedyn", "dissip", "vlspltqs", "aeropt", "ajsec",
"calltherm", "clmain", "cltrac", "cltracrn", "concvl", "conema3",
"conflx", "fisrtilp", "newmicro", "nuage", "diagcld1", "diagcld2",
"drag_noro", "lift_noro", "SUGWD", "physiq", "phytrac", "radlwsw", "thermcell".

Removed the case "ierr == 0" in "abort_gcm"; moved call to "histclo"
and messages for end of run from "abort_gcm" to "gcm"; replaced call
to "abort_gcm" in "leapfrog" by exit from outer loop.

In "calfis": removed argument "pp" and variable "unskap"; changed
"pksurcp" from scalar to rank 2; use "pressure_var"; rewrote
computation of "zplev", "zplay", "ztfi", "pcvgt" using "dyn_phy";
added computation of "pls".

Removed unused variable in "dynredem0".

In "exner_hyb": changed "dellta" from scalar to rank 1; replaced call
to "ssum" by call to "sum"; removed variables "xpn" and "xps";
replaced some loops by array expressions.

In "leapfrog": use "pressure_var"; deleted variables "p", "longcles".

Converted common blocks "YOECUMF", "YOEGWD" to modules.

Removed argument "pplay" in "cvltr", "diagetpq", "nflxtr".

Created module "raddimlw" from include file "raddimlw.h".

Corrected call to "new_unit" in "test_disvert".
1 guez 3 !
2     ! $Header: /home/cvsroot/LMDZ4/libf/phylmd/orografi.F,v 1.4 2005/12/01 11:27:29 fairhead Exp $
3     !
4     SUBROUTINE drag_noro (nlon,nlev,dtime,paprs,pplay,
5     e pmea,pstd, psig, pgam, pthe,ppic,pval,
6     e kgwd,kdx,ktest,
7     e t, u, v,
8     s pulow, pvlow, pustr, pvstr,
9     s d_t, d_u, d_v)
10     c
11     use dimens_m
12     use dimphy
13     use YOMCST
14     IMPLICIT none
15     c======================================================================
16     c Auteur(s): F.Lott (LMD/CNRS) date: 19950201
17     c Objet: Frottement de la montagne Interface
18     c======================================================================
19     c Arguments:
20     c dtime---input-R- pas d'integration (s)
21     c paprs---input-R-pression pour chaque inter-couche (en Pa)
22     c pplay---input-R-pression pour le mileu de chaque couche (en Pa)
23     c t-------input-R-temperature (K)
24     c u-------input-R-vitesse horizontale (m/s)
25     c v-------input-R-vitesse horizontale (m/s)
26     c
27     c d_t-----output-R-increment de la temperature
28     c d_u-----output-R-increment de la vitesse u
29     c d_v-----output-R-increment de la vitesse v
30     c======================================================================
31     c
32     c ARGUMENTS
33     c
34     INTEGER nlon,nlev
35     REAL dtime
36     REAL, intent(in):: paprs(klon,klev+1)
37 guez 10 REAL, intent(in):: pplay(klon,klev)
38 guez 3 REAL pmea(nlon),pstd(nlon),psig(nlon),pgam(nlon),pthe(nlon)
39     REAL ppic(nlon),pval(nlon)
40     REAL pulow(nlon),pvlow(nlon),pustr(nlon),pvstr(nlon)
41     REAL t(nlon,nlev), u(nlon,nlev), v(nlon,nlev)
42     REAL d_t(nlon,nlev), d_u(nlon,nlev), d_v(nlon,nlev)
43     c
44     INTEGER i, k, kgwd, kdx(nlon), ktest(nlon)
45     c
46     c Variables locales:
47     c
48     REAL zgeom(klon,klev)
49     REAL pdtdt(klon,klev), pdudt(klon,klev), pdvdt(klon,klev)
50     REAL pt(klon,klev), pu(klon,klev), pv(klon,klev)
51     REAL papmf(klon,klev),papmh(klon,klev+1)
52     c
53     c initialiser les variables de sortie (pour securite)
54     c
55     DO i = 1,klon
56     pulow(i) = 0.0
57     pvlow(i) = 0.0
58     pustr(i) = 0.0
59     pvstr(i) = 0.0
60     ENDDO
61     DO k = 1, klev
62     DO i = 1, klon
63     d_t(i,k) = 0.0
64     d_u(i,k) = 0.0
65     d_v(i,k) = 0.0
66     pdudt(i,k)=0.0
67     pdvdt(i,k)=0.0
68     pdtdt(i,k)=0.0
69     ENDDO
70     ENDDO
71     c
72     c preparer les variables d'entree (attention: l'ordre des niveaux
73     c verticaux augmente du haut vers le bas)
74     c
75     DO k = 1, klev
76     DO i = 1, klon
77     pt(i,k) = t(i,klev-k+1)
78     pu(i,k) = u(i,klev-k+1)
79     pv(i,k) = v(i,klev-k+1)
80     papmf(i,k) = pplay(i,klev-k+1)
81     ENDDO
82     ENDDO
83     DO k = 1, klev+1
84     DO i = 1, klon
85     papmh(i,k) = paprs(i,klev-k+2)
86     ENDDO
87     ENDDO
88     DO i = 1, klon
89     zgeom(i,klev) = RD * pt(i,klev)
90     . * LOG(papmh(i,klev+1)/papmf(i,klev))
91     ENDDO
92     DO k = klev-1, 1, -1
93     DO i = 1, klon
94     zgeom(i,k) = zgeom(i,k+1) + RD * (pt(i,k)+pt(i,k+1))/2.0
95     . * LOG(papmf(i,k+1)/papmf(i,k))
96     ENDDO
97     ENDDO
98     c
99     c appeler la routine principale
100     c
101     CALL orodrag(klon,klev,kgwd,kdx,ktest,
102     . dtime,
103     . papmh, papmf, zgeom,
104     . pt, pu, pv,
105     . pmea, pstd, psig, pgam, pthe, ppic,pval,
106     . pulow,pvlow,
107     . pdudt,pdvdt,pdtdt)
108     C
109     DO k = 1, klev
110     DO i = 1, klon
111     d_u(i,klev+1-k) = dtime*pdudt(i,k)
112     d_v(i,klev+1-k) = dtime*pdvdt(i,k)
113     d_t(i,klev+1-k) = dtime*pdtdt(i,k)
114     pustr(i) = pustr(i)
115     cIM BUG . +rg*pdudt(i,k)*(papmh(i,k+1)-papmh(i,k))
116     . +pdudt(i,k)*(papmh(i,k+1)-papmh(i,k))/RG
117     pvstr(i) = pvstr(i)
118     cIM BUG . +rg*pdvdt(i,k)*(papmh(i,k+1)-papmh(i,k))
119     . +pdvdt(i,k)*(papmh(i,k+1)-papmh(i,k))/RG
120     ENDDO
121     ENDDO
122     c
123     RETURN
124     END
125     SUBROUTINE orodrag( nlon,nlev
126     i , kgwd, kdx, ktest
127     r , ptsphy
128     r , paphm1,papm1,pgeom1,ptm1,pum1,pvm1
129     r , pmea, pstd, psig, pgamma, ptheta, ppic, pval
130     c outputs
131     r , pulow,pvlow
132     r , pvom,pvol,pte )
133    
134     use dimens_m
135     use dimphy
136     use YOMCST
137 guez 10 use yoegwd
138 guez 3 implicit none
139    
140     c
141     c
142     c**** *gwdrag* - does the gravity wave parametrization.
143     c
144     c purpose.
145     c --------
146     c
147     c this routine computes the physical tendencies of the
148     c prognostic variables u,v and t due to vertical transports by
149     c subgridscale orographically excited gravity waves
150     c
151     c** interface.
152     c ----------
153     c called from *callpar*.
154     c
155     c the routine takes its input from the long-term storage:
156     c u,v,t and p at t-1.
157     c
158     c explicit arguments :
159     c --------------------
160     c ==== inputs ===
161     c ==== outputs ===
162     c
163     c implicit arguments : none
164     c --------------------
165     c
166     c implicit logical (l)
167     c
168     c method.
169     c -------
170     c
171     c externals.
172     c ----------
173     integer ismin, ismax
174     external ismin, ismax
175     c
176     c reference.
177     c ----------
178     c
179     c author.
180     c -------
181     c m.miller + b.ritter e.c.m.w.f. 15/06/86.
182     c
183     c f.lott + m. miller e.c.m.w.f. 22/11/94
184     c-----------------------------------------------------------------------
185     c
186     c
187     c-----------------------------------------------------------------------
188     c
189     c* 0.1 arguments
190     c ---------
191     c
192     c
193     integer nlon, nlev, klevm1
194     integer kgwd, jl, ilevp1, jk, ji
195     real zdelp, ztemp, zforc, ztend
196     real rover, zb, zc, zconb, zabsv
197     real zzd1, ratio, zbet, zust,zvst, zdis
198     real pte(nlon,nlev),
199     * pvol(nlon,nlev),
200     * pvom(nlon,nlev),
201     * pulow(klon),
202     * pvlow(klon)
203     real pum1(nlon,nlev),
204     * pvm1(nlon,nlev),
205     * ptm1(nlon,nlev),
206     * pmea(nlon),pstd(nlon),psig(nlon),
207     * pgamma(nlon),ptheta(nlon),ppic(nlon),pval(nlon),
208     * pgeom1(nlon,nlev),
209     * papm1(nlon,nlev),
210     * paphm1(nlon,nlev+1)
211     c
212     integer kdx(nlon),ktest(nlon)
213     c-----------------------------------------------------------------------
214     c
215     c* 0.2 local arrays
216     c ------------
217     integer isect(klon),
218     * icrit(klon),
219     * ikcrith(klon),
220     * ikenvh(klon),
221     * iknu(klon),
222     * iknu2(klon),
223     * ikcrit(klon),
224     * ikhlim(klon)
225     c
226     real ztau(klon,klev+1),
227     $ ztauf(klon,klev+1),
228     * zstab(klon,klev+1),
229     * zvph(klon,klev+1),
230     * zrho(klon,klev+1),
231     * zri(klon,klev+1),
232     * zpsi(klon,klev+1),
233     * zzdep(klon,klev)
234     real zdudt(klon),
235     * zdvdt(klon),
236     * zdtdt(klon),
237     * zdedt(klon),
238     * zvidis(klon),
239     * znu(klon),
240     * zd1(klon),
241     * zd2(klon),
242     * zdmod(klon)
243     real ztmst, ptsphy, zrtmst
244     c
245     c------------------------------------------------------------------
246     c
247     c* 1. initialization
248     c --------------
249     c
250     100 continue
251     c
252     c ------------------------------------------------------------------
253     c
254     c* 1.1 computational constants
255     c -----------------------
256     c
257     110 continue
258     c
259     c ztmst=twodt
260     c if(nstep.eq.nstart) ztmst=0.5*twodt
261     klevm1=klev-1
262     ztmst=ptsphy
263     zrtmst=1./ztmst
264     c ------------------------------------------------------------------
265     c
266     120 continue
267     c
268     c ------------------------------------------------------------------
269     c
270     c* 1.3 check whether row contains point for printing
271     c ---------------------------------------------
272     c
273     130 continue
274     c
275     c ------------------------------------------------------------------
276     c
277     c* 2. precompute basic state variables.
278     c* ---------- ----- ----- ----------
279     c* define low level wind, project winds in plane of
280     c* low level wind, determine sector in which to take
281     c* the variance and set indicator for critical levels.
282     c
283     200 continue
284     c
285     c
286     c
287     call orosetup
288     * ( nlon, ktest
289     * , ikcrit, ikcrith, icrit, ikenvh,iknu,iknu2
290     * , paphm1, papm1 , pum1 , pvm1 , ptm1 , pgeom1, pstd
291     * , zrho , zri , zstab , ztau , zvph , zpsi, zzdep
292     * , pulow, pvlow
293     * , ptheta,pgamma,pmea,ppic,pval,znu ,zd1, zd2, zdmod )
294     c
295     c
296     c
297     c***********************************************************
298     c
299     c
300     c* 3. compute low level stresses using subcritical and
301     c* supercritical forms.computes anisotropy coefficient
302     c* as measure of orographic twodimensionality.
303     c
304     300 continue
305     c
306     call gwstress
307     * ( nlon , nlev
308     * , ktest , icrit, ikenvh, iknu
309     * , zrho , zstab, zvph , pstd, psig, pmea, ppic
310     * , ztau
311     * , pgeom1,zdmod)
312     c
313     c
314     c* 4. compute stress profile.
315     c* ------- ------ --------
316     c
317     400 continue
318     c
319     c
320     call gwprofil
321     * ( nlon , nlev
322     * , kgwd , kdx , ktest
323     * , ikcrith, icrit
324     * , paphm1, zrho , zstab , zvph
325     * , zri , ztau
326     * , zdmod , psig , pstd)
327     c
328     c
329     c* 5. compute tendencies.
330     c* -------------------
331     c
332     500 continue
333     c
334     c explicit solution at all levels for the gravity wave
335     c implicit solution for the blocked levels
336    
337     do 510 jl=1,klon
338     zvidis(jl)=0.0
339     zdudt(jl)=0.0
340     zdvdt(jl)=0.0
341     zdtdt(jl)=0.0
342     510 continue
343     c
344     ilevp1=klev+1
345     c
346     c
347     do 524 jk=1,klev
348     c
349     c
350     c do 523 jl=1,kgwd
351     c ji=kdx(jl)
352     c Modif vectorisation 02/04/2004
353     do 523 ji=1,klon
354     if(ktest(ji).eq.1) then
355    
356     zdelp=paphm1(ji,jk+1)-paphm1(ji,jk)
357     ztemp=-rg*(ztau(ji,jk+1)-ztau(ji,jk))/(zvph(ji,ilevp1)*zdelp)
358     zdudt(ji)=(pulow(ji)*zd1(ji)-pvlow(ji)*zd2(ji))*ztemp/zdmod(ji)
359     zdvdt(ji)=(pvlow(ji)*zd1(ji)+pulow(ji)*zd2(ji))*ztemp/zdmod(ji)
360     c
361     c controle des overshoots:
362     c
363     zforc=sqrt(zdudt(ji)**2+zdvdt(ji)**2)+1.E-12
364     ztend=sqrt(pum1(ji,jk)**2+pvm1(ji,jk)**2)/ztmst+1.E-12
365     rover=0.25
366     if(zforc.ge.rover*ztend)then
367     zdudt(ji)=rover*ztend/zforc*zdudt(ji)
368     zdvdt(ji)=rover*ztend/zforc*zdvdt(ji)
369     endif
370     c
371     c fin du controle des overshoots
372     c
373     if(jk.ge.ikenvh(ji)) then
374     zb=1.0-0.18*pgamma(ji)-0.04*pgamma(ji)**2
375     zc=0.48*pgamma(ji)+0.3*pgamma(ji)**2
376     zconb=2.*ztmst*gkwake*psig(ji)/(4.*pstd(ji))
377     zabsv=sqrt(pum1(ji,jk)**2+pvm1(ji,jk)**2)/2.
378     zzd1=zb*cos(zpsi(ji,jk))**2+zc*sin(zpsi(ji,jk))**2
379     ratio=(cos(zpsi(ji,jk))**2+pgamma(ji)*sin(zpsi(ji,jk))**2)/
380     * (pgamma(ji)*cos(zpsi(ji,jk))**2+sin(zpsi(ji,jk))**2)
381     zbet=max(0.,2.-1./ratio)*zconb*zzdep(ji,jk)*zzd1*zabsv
382     c
383     c simplement oppose au vent
384     c
385     zdudt(ji)=-pum1(ji,jk)/ztmst
386     zdvdt(ji)=-pvm1(ji,jk)/ztmst
387     c
388     c projection dans la direction de l'axe principal de l'orographie
389     cmod zdudt(ji)=-(pum1(ji,jk)*cos(ptheta(ji)*rpi/180.)
390     cmod * +pvm1(ji,jk)*sin(ptheta(ji)*rpi/180.))
391     cmod * *cos(ptheta(ji)*rpi/180.)/ztmst
392     cmod zdvdt(ji)=-(pum1(ji,jk)*cos(ptheta(ji)*rpi/180.)
393     cmod * +pvm1(ji,jk)*sin(ptheta(ji)*rpi/180.))
394     cmod * *sin(ptheta(ji)*rpi/180.)/ztmst
395     zdudt(ji)=zdudt(ji)*(zbet/(1.+zbet))
396     zdvdt(ji)=zdvdt(ji)*(zbet/(1.+zbet))
397     end if
398     pvom(ji,jk)=zdudt(ji)
399     pvol(ji,jk)=zdvdt(ji)
400     zust=pum1(ji,jk)+ztmst*zdudt(ji)
401     zvst=pvm1(ji,jk)+ztmst*zdvdt(ji)
402     zdis=0.5*(pum1(ji,jk)**2+pvm1(ji,jk)**2-zust**2-zvst**2)
403     zdedt(ji)=zdis/ztmst
404     zvidis(ji)=zvidis(ji)+zdis*zdelp
405     zdtdt(ji)=zdedt(ji)/rcpd
406     c pte(ji,jk)=zdtdt(ji)
407     c
408     c ENCORE UN TRUC POUR EVITER LES EXPLOSIONS
409     c
410     pte(ji,jk)=0.0
411    
412     endif
413     523 continue
414    
415     524 continue
416     c
417     c
418     return
419     end
420     SUBROUTINE orosetup
421     * ( nlon , ktest
422     * , kkcrit, kkcrith, kcrit
423     * , kkenvh, kknu , kknu2
424     * , paphm1, papm1 , pum1 , pvm1 , ptm1 , pgeom1, pstd
425     * , prho , pri , pstab , ptau , pvph ,ppsi, pzdep
426     * , pulow , pvlow
427     * , ptheta, pgamma, pmea, ppic, pval
428     * , pnu , pd1 , pd2 ,pdmod )
429     c
430     c**** *gwsetup*
431     c
432     c purpose.
433     c --------
434     c
435     c** interface.
436     c ----------
437     c from *orodrag*
438     c
439     c explicit arguments :
440     c --------------------
441     c ==== inputs ===
442     c ==== outputs ===
443     c
444     c implicit arguments : none
445     c --------------------
446     c
447     c method.
448     c -------
449     c
450     c
451     c externals.
452     c ----------
453     c
454     c
455     c reference.
456     c ----------
457     c
458     c see ecmwf research department documentation of the "i.f.s."
459     c
460     c author.
461     c -------
462     c
463     c modifications.
464     c --------------
465     c f.lott for the new-gwdrag scheme november 1993
466     c
467     c-----------------------------------------------------------------------
468     use dimens_m
469     use dimphy
470     use YOMCST
471 guez 10 use yoegwd
472 guez 3 implicit none
473     c
474    
475    
476     c-----------------------------------------------------------------------
477     c
478     c* 0.1 arguments
479     c ---------
480     c
481     integer nlon
482     integer jl, jk
483     real zdelp
484    
485     integer kkcrit(nlon),kkcrith(nlon),kcrit(nlon),
486     * ktest(nlon),kkenvh(nlon)
487    
488     c
489     real paphm1(nlon,klev+1),papm1(nlon,klev),pum1(nlon,klev),
490     * pvm1(nlon,klev),ptm1(nlon,klev),pgeom1(nlon,klev),
491     * prho(nlon,klev+1),pri(nlon,klev+1),pstab(nlon,klev+1),
492     * ptau(nlon,klev+1),pvph(nlon,klev+1),ppsi(nlon,klev+1),
493     * pzdep(nlon,klev)
494     real pulow(nlon),pvlow(nlon),ptheta(nlon),pgamma(nlon),pnu(nlon),
495     * pd1(nlon),pd2(nlon),pdmod(nlon)
496     real pstd(nlon),pmea(nlon),ppic(nlon),pval(nlon)
497     c
498     c-----------------------------------------------------------------------
499     c
500     c* 0.2 local arrays
501     c ------------
502     c
503     c
504     integer ilevm1, ilevm2, ilevh
505     real zcons1, zcons2,zcons3, zhgeo
506     real zu, zphi, zvt1,zvt2, zst, zvar, zdwind, zwind
507     real zstabm, zstabp, zrhom, zrhop, alpha
508     real zggeenv, zggeom1,zgvar
509     logical lo
510     logical ll1(klon,klev+1)
511     integer kknu(klon),kknu2(klon),kknub(klon),kknul(klon),
512     * kentp(klon),ncount(klon)
513     c
514     real zhcrit(klon,klev),zvpf(klon,klev),
515     * zdp(klon,klev)
516     real znorm(klon),zb(klon),zc(klon),
517     * zulow(klon),zvlow(klon),znup(klon),znum(klon)
518     c
519     c ------------------------------------------------------------------
520     c
521     c* 1. initialization
522     c --------------
523     c
524     c print *,' entree gwsetup'
525     100 continue
526     c
527     c ------------------------------------------------------------------
528     c
529     c* 1.1 computational constants
530     c -----------------------
531     c
532     110 continue
533     c
534     ilevm1=klev-1
535     ilevm2=klev-2
536     ilevh =klev/3
537     c
538     zcons1=1./rd
539     cold zcons2=g**2/cpd
540     zcons2=rg**2/rcpd
541     cold zcons3=1.5*api
542     zcons3=1.5*rpi
543     c
544     c
545     c ------------------------------------------------------------------
546     c
547     c* 2.
548     c --------------
549     c
550     200 continue
551     c
552     c ------------------------------------------------------------------
553     c
554     c* 2.1 define low level wind, project winds in plane of
555     c* low level wind, determine sector in which to take
556     c* the variance and set indicator for critical levels.
557     c
558     c
559     c
560     do 2001 jl=1,klon
561     kknu(jl) =klev
562     kknu2(jl) =klev
563     kknub(jl) =klev
564     kknul(jl) =klev
565     pgamma(jl) =max(pgamma(jl),gtsec)
566     ll1(jl,klev+1)=.false.
567     2001 continue
568     c
569     c Ajouter une initialisation (L. Li, le 23fev99):
570     c
571     do jk=klev,ilevh,-1
572     do jl=1,klon
573     ll1(jl,jk)= .FALSE.
574     ENDDO
575     ENDDO
576     c
577     c* define top of low level flow
578     c ----------------------------
579     do 2002 jk=klev,ilevh,-1
580     do 2003 jl=1,klon
581     lo=(paphm1(jl,jk)/paphm1(jl,klev+1)).ge.gsigcr
582     if(lo) then
583     kkcrit(jl)=jk
584     endif
585     zhcrit(jl,jk)=ppic(jl)
586     zhgeo=pgeom1(jl,jk)/rg
587     ll1(jl,jk)=(zhgeo.gt.zhcrit(jl,jk))
588     if(ll1(jl,jk).neqv.ll1(jl,jk+1)) then
589     kknu(jl)=jk
590     endif
591     if(.not.ll1(jl,ilevh))kknu(jl)=ilevh
592     2003 continue
593     2002 continue
594     do 2004 jk=klev,ilevh,-1
595     do 2005 jl=1,klon
596     zhcrit(jl,jk)=ppic(jl)-pval(jl)
597     zhgeo=pgeom1(jl,jk)/rg
598     ll1(jl,jk)=(zhgeo.gt.zhcrit(jl,jk))
599     if(ll1(jl,jk).neqv.ll1(jl,jk+1)) then
600     kknu2(jl)=jk
601     endif
602     if(.not.ll1(jl,ilevh))kknu2(jl)=ilevh
603     2005 continue
604     2004 continue
605     do 2006 jk=klev,ilevh,-1
606     do 2007 jl=1,klon
607     zhcrit(jl,jk)=amax1(ppic(jl)-pmea(jl),pmea(jl)-pval(jl))
608     zhgeo=pgeom1(jl,jk)/rg
609     ll1(jl,jk)=(zhgeo.gt.zhcrit(jl,jk))
610     if(ll1(jl,jk).neqv.ll1(jl,jk+1)) then
611     kknub(jl)=jk
612     endif
613     if(.not.ll1(jl,ilevh))kknub(jl)=ilevh
614     2007 continue
615     2006 continue
616     c
617     do 2010 jl=1,klon
618     kknu(jl)=min(kknu(jl),nktopg)
619     kknu2(jl)=min(kknu2(jl),nktopg)
620     kknub(jl)=min(kknub(jl),nktopg)
621     kknul(jl)=klev
622     2010 continue
623     c
624    
625     210 continue
626     c
627     c
628     cc* initialize various arrays
629     c
630     do 2107 jl=1,klon
631     prho(jl,klev+1) =0.0
632     pstab(jl,klev+1) =0.0
633     pstab(jl,1) =0.0
634     pri(jl,klev+1) =9999.0
635     ppsi(jl,klev+1) =0.0
636     pri(jl,1) =0.0
637     pvph(jl,1) =0.0
638     pulow(jl) =0.0
639     pvlow(jl) =0.0
640     zulow(jl) =0.0
641     zvlow(jl) =0.0
642     kkcrith(jl) =klev
643     kkenvh(jl) =klev
644     kentp(jl) =klev
645     kcrit(jl) =1
646     ncount(jl) =0
647     ll1(jl,klev+1) =.false.
648     2107 continue
649     c
650     c* define low-level flow
651     c ---------------------
652     c
653     do 223 jk=klev,2,-1
654     do 222 jl=1,klon
655     if(ktest(jl).eq.1) then
656     zdp(jl,jk)=papm1(jl,jk)-papm1(jl,jk-1)
657     prho(jl,jk)=2.*paphm1(jl,jk)*zcons1/(ptm1(jl,jk)+ptm1(jl,jk-1))
658     pstab(jl,jk)=2.*zcons2/(ptm1(jl,jk)+ptm1(jl,jk-1))*
659     * (1.-rcpd*prho(jl,jk)*(ptm1(jl,jk)-ptm1(jl,jk-1))/zdp(jl,jk))
660     pstab(jl,jk)=max(pstab(jl,jk),gssec)
661     endif
662     222 continue
663     223 continue
664     c
665     c********************************************************************
666     c
667     c* define blocked flow
668     c -------------------
669     do 2115 jk=klev,ilevh,-1
670     do 2116 jl=1,klon
671     if(jk.ge.kknub(jl).and.jk.le.kknul(jl)) then
672     pulow(jl)=pulow(jl)+pum1(jl,jk)*(paphm1(jl,jk+1)-paphm1(jl,jk))
673     pvlow(jl)=pvlow(jl)+pvm1(jl,jk)*(paphm1(jl,jk+1)-paphm1(jl,jk))
674     end if
675     2116 continue
676     2115 continue
677     do 2110 jl=1,klon
678     pulow(jl)=pulow(jl)/(paphm1(jl,kknul(jl)+1)-paphm1(jl,kknub(jl)))
679     pvlow(jl)=pvlow(jl)/(paphm1(jl,kknul(jl)+1)-paphm1(jl,kknub(jl)))
680     znorm(jl)=max(sqrt(pulow(jl)**2+pvlow(jl)**2),gvsec)
681     pvph(jl,klev+1)=znorm(jl)
682     2110 continue
683     c
684     c******* setup orography axes and define plane of profiles *******
685     c
686     do 2112 jl=1,klon
687     lo=(pulow(jl).lt.gvsec).and.(pulow(jl).ge.-gvsec)
688     if(lo) then
689     zu=pulow(jl)+2.*gvsec
690     else
691     zu=pulow(jl)
692     endif
693     zphi=atan(pvlow(jl)/zu)
694     ppsi(jl,klev+1)=ptheta(jl)*rpi/180.-zphi
695     zb(jl)=1.-0.18*pgamma(jl)-0.04*pgamma(jl)**2
696     zc(jl)=0.48*pgamma(jl)+0.3*pgamma(jl)**2
697     pd1(jl)=zb(jl)-(zb(jl)-zc(jl))*(sin(ppsi(jl,klev+1))**2)
698     pd2(jl)=(zb(jl)-zc(jl))*sin(ppsi(jl,klev+1))*cos(ppsi(jl,klev+1))
699     pdmod(jl)=sqrt(pd1(jl)**2+pd2(jl)**2)
700     2112 continue
701     c
702     c ************ define flow in plane of lowlevel stress *************
703     c
704     do 213 jk=1,klev
705     do 212 jl=1,klon
706     if(ktest(jl).eq.1) then
707     zvt1 =pulow(jl)*pum1(jl,jk)+pvlow(jl)*pvm1(jl,jk)
708     zvt2 =-pvlow(jl)*pum1(jl,jk)+pulow(jl)*pvm1(jl,jk)
709     zvpf(jl,jk)=(zvt1*pd1(jl)+zvt2*pd2(jl))/(znorm(jl)*pdmod(jl))
710     endif
711     ptau(jl,jk) =0.0
712     pzdep(jl,jk) =0.0
713     ppsi(jl,jk) =0.0
714     ll1(jl,jk) =.false.
715     212 continue
716     213 continue
717     do 215 jk=2,klev
718     do 214 jl=1,klon
719     if(ktest(jl).eq.1) then
720     zdp(jl,jk)=papm1(jl,jk)-papm1(jl,jk-1)
721     pvph(jl,jk)=((paphm1(jl,jk)-papm1(jl,jk-1))*zvpf(jl,jk)+
722     * (papm1(jl,jk)-paphm1(jl,jk))*zvpf(jl,jk-1))
723     * /zdp(jl,jk)
724     if(pvph(jl,jk).lt.gvsec) then
725     pvph(jl,jk)=gvsec
726     kcrit(jl)=jk
727     endif
728     endif
729     214 continue
730     215 continue
731     c
732     c
733     c* 2.2 brunt-vaisala frequency and density at half levels.
734     c
735     220 continue
736     c
737     do 2211 jk=ilevh,klev
738     do 221 jl=1,klon
739     if(ktest(jl).eq.1) then
740     if(jk.ge.(kknub(jl)+1).and.jk.le.kknul(jl)) then
741     zst=zcons2/ptm1(jl,jk)*(1.-rcpd*prho(jl,jk)*
742     * (ptm1(jl,jk)-ptm1(jl,jk-1))/zdp(jl,jk))
743     pstab(jl,klev+1)=pstab(jl,klev+1)+zst*zdp(jl,jk)
744     pstab(jl,klev+1)=max(pstab(jl,klev+1),gssec)
745     prho(jl,klev+1)=prho(jl,klev+1)+paphm1(jl,jk)*2.*zdp(jl,jk)
746     * *zcons1/(ptm1(jl,jk)+ptm1(jl,jk-1))
747     endif
748     endif
749     221 continue
750     2211 continue
751     c
752     do 2212 jl=1,klon
753     pstab(jl,klev+1)=pstab(jl,klev+1)/(papm1(jl,kknul(jl))
754     * -papm1(jl,kknub(jl)))
755     prho(jl,klev+1)=prho(jl,klev+1)/(papm1(jl,kknul(jl))
756     * -papm1(jl,kknub(jl)))
757     zvar=pstd(jl)
758     2212 continue
759     c
760     c* 2.3 mean flow richardson number.
761     c* and critical height for froude layer
762     c
763     230 continue
764     c
765     do 232 jk=2,klev
766     do 231 jl=1,klon
767     if(ktest(jl).eq.1) then
768     zdwind=max(abs(zvpf(jl,jk)-zvpf(jl,jk-1)),gvsec)
769     pri(jl,jk)=pstab(jl,jk)*(zdp(jl,jk)
770     * /(rg*prho(jl,jk)*zdwind))**2
771     pri(jl,jk)=max(pri(jl,jk),grcrit)
772     endif
773     231 continue
774     232 continue
775    
776     c
777     c
778     c* define top of 'envelope' layer
779     c ----------------------------
780    
781     do 233 jl=1,klon
782     pnu (jl)=0.0
783     znum(jl)=0.0
784     233 continue
785    
786     do 234 jk=2,klev-1
787     do 234 jl=1,klon
788    
789     if(ktest(jl).eq.1) then
790    
791     if (jk.ge.kknub(jl)) then
792    
793     znum(jl)=pnu(jl)
794     zwind=(pulow(jl)*pum1(jl,jk)+pvlow(jl)*pvm1(jl,jk))/
795     * max(sqrt(pulow(jl)**2+pvlow(jl)**2),gvsec)
796     zwind=max(sqrt(zwind**2),gvsec)
797     zdelp=paphm1(jl,jk+1)-paphm1(jl,jk)
798     zstabm=sqrt(max(pstab(jl,jk ),gssec))
799     zstabp=sqrt(max(pstab(jl,jk+1),gssec))
800     zrhom=prho(jl,jk )
801     zrhop=prho(jl,jk+1)
802     pnu(jl) = pnu(jl) + (zdelp/rg)*
803     * ((zstabp/zrhop+zstabm/zrhom)/2.)/zwind
804     if((znum(jl).le.gfrcrit).and.(pnu(jl).gt.gfrcrit)
805     * .and.(kkenvh(jl).eq.klev))
806     * kkenvh(jl)=jk
807    
808     endif
809    
810     endif
811    
812     234 continue
813    
814     c calculation of a dynamical mixing height for the breaking
815     c of gravity waves:
816    
817    
818     do 235 jl=1,klon
819     znup(jl)=0.0
820     znum(jl)=0.0
821     235 continue
822    
823     do 236 jk=klev-1,2,-1
824     do 236 jl=1,klon
825    
826     if(ktest(jl).eq.1) then
827    
828     znum(jl)=znup(jl)
829     zwind=(pulow(jl)*pum1(jl,jk)+pvlow(jl)*pvm1(jl,jk))/
830     * max(sqrt(pulow(jl)**2+pvlow(jl)**2),gvsec)
831     zwind=max(sqrt(zwind**2),gvsec)
832     zdelp=paphm1(jl,jk+1)-paphm1(jl,jk)
833     zstabm=sqrt(max(pstab(jl,jk ),gssec))
834     zstabp=sqrt(max(pstab(jl,jk+1),gssec))
835     zrhom=prho(jl,jk )
836     zrhop=prho(jl,jk+1)
837     znup(jl) = znup(jl) + (zdelp/rg)*
838     * ((zstabp/zrhop+zstabm/zrhom)/2.)/zwind
839     if((znum(jl).le.rpi/2.).and.(znup(jl).gt.rpi/2.)
840     * .and.(kkcrith(jl).eq.klev))
841     * kkcrith(jl)=jk
842    
843     endif
844    
845     236 continue
846    
847     do 237 jl=1,klon
848     kkcrith(jl)=min0(kkcrith(jl),kknu2(jl))
849     kkcrith(jl)=max0(kkcrith(jl),ilevh*2)
850     237 continue
851     c
852     c directional info for flow blocking *************************
853     c
854     do 251 jk=ilevh,klev
855     do 252 jl=1,klon
856     if(jk.ge.kkenvh(jl)) then
857     lo=(pum1(jl,jk).lt.gvsec).and.(pum1(jl,jk).ge.-gvsec)
858     if(lo) then
859     zu=pum1(jl,jk)+2.*gvsec
860     else
861     zu=pum1(jl,jk)
862     endif
863     zphi=atan(pvm1(jl,jk)/zu)
864     ppsi(jl,jk)=ptheta(jl)*rpi/180.-zphi
865     end if
866     252 continue
867     251 continue
868     c forms the vertical 'leakiness' **************************
869    
870     alpha=3.
871    
872     do 254 jk=ilevh,klev
873     do 253 jl=1,klon
874     if(jk.ge.kkenvh(jl)) then
875     zggeenv=amax1(1.,
876     * (pgeom1(jl,kkenvh(jl))+pgeom1(jl,kkenvh(jl)-1))/2.)
877     zggeom1=amax1(pgeom1(jl,jk),1.)
878     zgvar=amax1(pstd(jl)*rg,1.)
879     cmod pzdep(jl,jk)=sqrt((zggeenv-zggeom1)/(zggeom1+zgvar))
880     pzdep(jl,jk)=(pgeom1(jl,kkenvh(jl)-1)-pgeom1(jl, jk))/
881     * (pgeom1(jl,kkenvh(jl)-1)-pgeom1(jl,klev))
882     end if
883     253 continue
884     254 continue
885    
886     260 continue
887    
888     return
889     end
890     SUBROUTINE gwstress
891     * ( nlon , nlev
892     * , ktest, kcrit, kkenvh
893     * , kknu
894     * , prho , pstab , pvph , pstd, psig
895     * , pmea , ppic , ptau
896     * , pgeom1 , pdmod )
897     c
898     c**** *gwstress*
899     c
900     c purpose.
901     c --------
902     c
903     c** interface.
904     c ----------
905     c call *gwstress* from *gwdrag*
906     c
907     c explicit arguments :
908     c --------------------
909     c ==== inputs ===
910     c ==== outputs ===
911     c
912     c implicit arguments : none
913     c --------------------
914     c
915     c method.
916     c -------
917     c
918     c
919     c externals.
920     c ----------
921     c
922     c
923     c reference.
924     c ----------
925     c
926     c see ecmwf research department documentation of the "i.f.s."
927     c
928     c author.
929     c -------
930     c
931     c modifications.
932     c --------------
933     c f. lott put the new gwd on ifs 22/11/93
934     c
935     c-----------------------------------------------------------------------
936     use dimens_m
937     use dimphy
938     use YOMCST
939 guez 10 use yoegwd
940 guez 3 implicit none
941    
942     c-----------------------------------------------------------------------
943     c
944     c* 0.1 arguments
945     c ---------
946     c
947     integer nlon, nlev
948     integer kcrit(nlon),
949     * ktest(nlon),kkenvh(nlon),kknu(nlon)
950     c
951     real prho(nlon,nlev+1),pstab(nlon,nlev+1),ptau(nlon,nlev+1),
952     * pvph(nlon,nlev+1),
953     * pgeom1(nlon,nlev),pstd(nlon)
954     c
955     real psig(nlon)
956     real pmea(nlon),ppic(nlon)
957     real pdmod(nlon)
958     c
959     c-----------------------------------------------------------------------
960     c
961     c* 0.2 local arrays
962     c ------------
963     integer jl
964     real zblock, zvar, zeff
965     logical lo
966     c
967     c-----------------------------------------------------------------------
968     c
969     c* 0.3 functions
970     c ---------
971     c ------------------------------------------------------------------
972     c
973     c* 1. initialization
974     c --------------
975     c
976     100 continue
977     c
978     c* 3.1 gravity wave stress.
979     c
980     300 continue
981     c
982     c
983     do 301 jl=1,klon
984     if(ktest(jl).eq.1) then
985    
986     c effective mountain height above the blocked flow
987    
988     if(kkenvh(jl).eq.klev)then
989     zblock=0.0
990     else
991     zblock=(pgeom1(jl,kkenvh(jl))+pgeom1(jl,kkenvh(jl)+1))/2./rg
992     endif
993    
994     zvar=ppic(jl)-pmea(jl)
995     zeff=amax1(0.,zvar-zblock)
996    
997     ptau(jl,klev+1)=prho(jl,klev+1)*gkdrag*psig(jl)*zeff**2
998     * /4./pstd(jl)*pvph(jl,klev+1)*pdmod(jl)*sqrt(pstab(jl,klev+1))
999    
1000     c too small value of stress or low level flow include critical level
1001     c or low level flow: gravity wave stress nul.
1002    
1003     lo=(ptau(jl,klev+1).lt.gtsec).or.(kcrit(jl).ge.kknu(jl))
1004     * .or.(pvph(jl,klev+1).lt.gvcrit)
1005     c if(lo) ptau(jl,klev+1)=0.0
1006    
1007     else
1008    
1009     ptau(jl,klev+1)=0.0
1010    
1011     endif
1012    
1013     301 continue
1014     c
1015     return
1016     end
1017     SUBROUTINE GWPROFIL
1018     * ( NLON, NLEV
1019     * , kgwd, kdx , ktest
1020     * , KKCRITH, KCRIT
1021     * , PAPHM1, PRHO , PSTAB , PVPH , PRI , PTAU
1022     * , pdmod , psig , pvar)
1023    
1024     C**** *GWPROFIL*
1025     C
1026     C PURPOSE.
1027     C --------
1028     C
1029     C** INTERFACE.
1030     C ----------
1031     C FROM *GWDRAG*
1032     C
1033     C EXPLICIT ARGUMENTS :
1034     C --------------------
1035     C ==== INPUTS ===
1036     C ==== OUTPUTS ===
1037     C
1038     C IMPLICIT ARGUMENTS : NONE
1039     C --------------------
1040     C
1041     C METHOD:
1042     C -------
1043     C THE STRESS PROFILE FOR GRAVITY WAVES IS COMPUTED AS FOLLOWS:
1044     C IT IS CONSTANT (NO GWD) AT THE LEVELS BETWEEN THE GROUND
1045     C AND THE TOP OF THE BLOCKED LAYER (KKENVH).
1046     C IT DECREASES LINEARLY WITH HEIGHTS FROM THE TOP OF THE
1047     C BLOCKED LAYER TO 3*VAROR (kKNU), TO SIMULATES LEE WAVES OR
1048     C NONLINEAR GRAVITY WAVE BREAKING.
1049     C ABOVE IT IS CONSTANT, EXCEPT WHEN THE WAVE ENCOUNTERS A CRITICAL
1050     C LEVEL (KCRIT) OR WHEN IT BREAKS.
1051     C
1052     C
1053     C
1054     C EXTERNALS.
1055     C ----------
1056     C
1057     C
1058     C REFERENCE.
1059     C ----------
1060     C
1061     C SEE ECMWF RESEARCH DEPARTMENT DOCUMENTATION OF THE "I.F.S."
1062     C
1063     C AUTHOR.
1064     C -------
1065     C
1066     C MODIFICATIONS.
1067     C --------------
1068     C PASSAGE OF THE NEW GWDRAG TO I.F.S. (F. LOTT, 22/11/93)
1069     C-----------------------------------------------------------------------
1070     use dimens_m
1071     use dimphy
1072     use YOMCST
1073 guez 10 use yoegwd
1074 guez 3 implicit none
1075     C
1076    
1077     C
1078    
1079    
1080     C-----------------------------------------------------------------------
1081     C
1082     C* 0.1 ARGUMENTS
1083     C ---------
1084     C
1085     integer nlon,nlev
1086     INTEGER KKCRITH(NLON),KCRIT(NLON)
1087     * ,kdx(nlon) , ktest(nlon)
1088    
1089     C
1090     REAL PAPHM1(NLON,NLEV+1), PSTAB(NLON,NLEV+1),
1091     * PRHO (NLON,NLEV+1), PVPH (NLON,NLEV+1),
1092     * PRI (NLON,NLEV+1), PTAU(NLON,NLEV+1)
1093    
1094     REAL pdmod (NLON) , psig(NLON),
1095     * pvar(NLON)
1096    
1097     C-----------------------------------------------------------------------
1098     C
1099     C* 0.2 LOCAL ARRAYS
1100     C ------------
1101     C
1102     integer ilevh, ji, kgwd, jl, jk
1103     real zsqr, zalfa, zriw, zdel, zb, zalpha,zdz2n
1104     real zdelp, zdelpt
1105     REAL ZDZ2 (KLON,KLEV) , ZNORM(KLON) , zoro(KLON)
1106     REAL ZTAU (KLON,KLEV+1)
1107     C
1108     C-----------------------------------------------------------------------
1109     C
1110     C* 1. INITIALIZATION
1111     C --------------
1112     C
1113     c print *,' entree gwprofil'
1114     100 CONTINUE
1115     C
1116     C
1117     C* COMPUTATIONAL CONSTANTS.
1118     C ------------- ----------
1119     C
1120     ilevh=KLEV/3
1121     C
1122     c DO 400 ji=1,kgwd
1123     c jl=kdx(ji)
1124     c Modif vectorisation 02/04/2004
1125     DO 400 jl=1,klon
1126     if (ktest(jl).eq.1) then
1127     Zoro(JL)=Psig(JL)*Pdmod(JL)/4./max(pvar(jl),1.0)
1128     ZTAU(JL,KLEV+1)=PTAU(JL,KLEV+1)
1129     endif
1130     400 CONTINUE
1131    
1132     C
1133     DO 430 JK=KLEV,2,-1
1134     C
1135     C
1136     C* 4.1 CONSTANT WAVE STRESS UNTIL TOP OF THE
1137     C BLOCKING LAYER.
1138     410 CONTINUE
1139     C
1140     c DO 411 ji=1,kgwd
1141     c jl=kdx(ji)
1142     c Modif vectorisation 02/04/2004
1143     do 411 jl=1,klon
1144     if (ktest(jl).eq.1) then
1145     IF(JK.GT.KKCRITH(JL)) THEN
1146     PTAU(JL,JK)=ZTAU(JL,KLEV+1)
1147     C ENDIF
1148     C IF(JK.EQ.KKCRITH(JL)) THEN
1149     ELSE
1150     PTAU(JL,JK)=GRAHILO*ZTAU(JL,KLEV+1)
1151     ENDIF
1152     endif
1153     411 CONTINUE
1154     C
1155     C* 4.15 CONSTANT SHEAR STRESS UNTIL THE TOP OF THE
1156     C LOW LEVEL FLOW LAYER.
1157     415 CONTINUE
1158     C
1159     C
1160     C* 4.2 WAVE DISPLACEMENT AT NEXT LEVEL.
1161     C
1162     420 CONTINUE
1163     C
1164     c DO 421 ji=1,kgwd
1165     c jl=kdx(ji)
1166     c Modif vectorisation 02/04/2004
1167     do 421 jl=1,klon
1168     if(ktest(jl).eq.1) then
1169     IF(JK.LT.KKCRITH(JL)) THEN
1170     ZNORM(JL)=gkdrag*PRHO(JL,JK)*SQRT(PSTAB(JL,JK))*PVPH(JL,JK)
1171     * *zoro(jl)
1172     ZDZ2(JL,JK)=PTAU(JL,JK+1)/max(ZNORM(JL),gssec)
1173     ENDIF
1174     endif
1175     421 CONTINUE
1176     C
1177     C* 4.3 WAVE RICHARDSON NUMBER, NEW WAVE DISPLACEMENT
1178     C* AND STRESS: BREAKING EVALUATION AND CRITICAL
1179     C LEVEL
1180     C
1181    
1182     c DO 431 ji=1,kgwd
1183     c jl=Kdx(ji)
1184     c Modif vectorisation 02/04/2004
1185     do 431 jl=1,klon
1186     if(ktest(jl).eq.1) then
1187    
1188     IF(JK.LT.KKCRITH(JL)) THEN
1189     IF((PTAU(JL,JK+1).LT.GTSEC).OR.(JK.LE.KCRIT(JL))) THEN
1190     PTAU(JL,JK)=0.0
1191     ELSE
1192     ZSQR=SQRT(PRI(JL,JK))
1193     ZALFA=SQRT(PSTAB(JL,JK)*ZDZ2(JL,JK))/PVPH(JL,JK)
1194     ZRIW=PRI(JL,JK)*(1.-ZALFA)/(1+ZALFA*ZSQR)**2
1195     IF(ZRIW.LT.GRCRIT) THEN
1196     ZDEL=4./ZSQR/GRCRIT+1./GRCRIT**2+4./GRCRIT
1197     ZB=1./GRCRIT+2./ZSQR
1198     ZALPHA=0.5*(-ZB+SQRT(ZDEL))
1199     ZDZ2N=(PVPH(JL,JK)*ZALPHA)**2/PSTAB(JL,JK)
1200     PTAU(JL,JK)=ZNORM(JL)*ZDZ2N
1201     ELSE
1202     PTAU(JL,JK)=ZNORM(JL)*ZDZ2(JL,JK)
1203     ENDIF
1204     PTAU(JL,JK)=MIN(PTAU(JL,JK),PTAU(JL,JK+1))
1205     ENDIF
1206     ENDIF
1207     endif
1208     431 CONTINUE
1209    
1210     430 CONTINUE
1211     440 CONTINUE
1212    
1213     C REORGANISATION OF THE STRESS PROFILE AT LOW LEVEL
1214    
1215     c DO 530 ji=1,kgwd
1216     c jl=kdx(ji)
1217     c Modif vectorisation 02/04/2004
1218     do 530 jl=1,klon
1219     if(ktest(jl).eq.1) then
1220     ZTAU(JL,KKCRITH(JL))=PTAU(JL,KKCRITH(JL))
1221     ZTAU(JL,NSTRA)=PTAU(JL,NSTRA)
1222     endif
1223     530 CONTINUE
1224    
1225     DO 531 JK=1,KLEV
1226    
1227     c DO 532 ji=1,kgwd
1228     c jl=kdx(ji)
1229     c Modif vectorisation 02/04/2004
1230     do 532 jl=1,klon
1231     if(ktest(jl).eq.1) then
1232    
1233    
1234     IF(JK.GT.KKCRITH(JL))THEN
1235    
1236     ZDELP=PAPHM1(JL,JK)-PAPHM1(JL,KLEV+1 )
1237     ZDELPT=PAPHM1(JL,KKCRITH(JL))-PAPHM1(JL,KLEV+1 )
1238     PTAU(JL,JK)=ZTAU(JL,KLEV+1 ) +
1239     . (ZTAU(JL,KKCRITH(JL))-ZTAU(JL,KLEV+1 ) )*
1240     . ZDELP/ZDELPT
1241    
1242     ENDIF
1243    
1244     endif
1245     532 CONTINUE
1246    
1247     C REORGANISATION IN THE STRATOSPHERE
1248    
1249     c DO 533 ji=1,kgwd
1250     c jl=kdx(ji)
1251     c Modif vectorisation 02/04/2004
1252     do 533 jl=1,klon
1253     if(ktest(jl).eq.1) then
1254    
1255    
1256     IF(JK.LT.NSTRA)THEN
1257    
1258     ZDELP =PAPHM1(JL,NSTRA)
1259     ZDELPT=PAPHM1(JL,JK)
1260     PTAU(JL,JK)=ZTAU(JL,NSTRA)*ZDELPT/ZDELP
1261    
1262     ENDIF
1263    
1264     endif
1265     533 CONTINUE
1266    
1267     C REORGANISATION IN THE TROPOSPHERE
1268    
1269     c DO 534 ji=1,kgwd
1270     c jl=kdx(ji)
1271     c Modif vectorisation 02/04/2004
1272     do 534 jl=1,klon
1273     if(ktest(jl).eq.1) then
1274    
1275    
1276     IF(JK.LT.KKCRITH(JL).AND.JK.GT.NSTRA)THEN
1277    
1278     ZDELP=PAPHM1(JL,JK)-PAPHM1(JL,KKCRITH(JL))
1279     ZDELPT=PAPHM1(JL,NSTRA)-PAPHM1(JL,KKCRITH(JL))
1280     PTAU(JL,JK)=ZTAU(JL,KKCRITH(JL)) +
1281     * (ZTAU(JL,NSTRA)-ZTAU(JL,KKCRITH(JL)))*ZDELP
1282     . /ZDELPT
1283    
1284     ENDIF
1285     endif
1286     534 CONTINUE
1287    
1288    
1289     531 CONTINUE
1290    
1291    
1292     RETURN
1293     END
1294     SUBROUTINE lift_noro (nlon,nlev,dtime,paprs,pplay,
1295     e plat,pmea,pstd, ppic,
1296     e ktest,
1297     e t, u, v,
1298     s pulow, pvlow, pustr, pvstr,
1299     s d_t, d_u, d_v)
1300     c
1301     use dimens_m
1302     use dimphy
1303     use YOMCST
1304     IMPLICIT none
1305     c======================================================================
1306     c Auteur(s): F.Lott (LMD/CNRS) date: 19950201
1307     c Objet: Frottement de la montagne Interface
1308     c======================================================================
1309     c Arguments:
1310     c dtime---input-R- pas d'integration (s)
1311     c paprs---input-R-pression pour chaque inter-couche (en Pa)
1312     c pplay---input-R-pression pour le mileu de chaque couche (en Pa)
1313     c t-------input-R-temperature (K)
1314     c u-------input-R-vitesse horizontale (m/s)
1315     c v-------input-R-vitesse horizontale (m/s)
1316     c
1317     c d_t-----output-R-increment de la temperature
1318     c d_u-----output-R-increment de la vitesse u
1319     c d_v-----output-R-increment de la vitesse v
1320     c======================================================================
1321     c
1322     c ARGUMENTS
1323     c
1324     INTEGER nlon,nlev
1325     REAL dtime
1326     REAL, intent(in):: paprs(klon,klev+1)
1327 guez 10 REAL, intent(in):: pplay(klon,klev)
1328 guez 3 REAL, intent(in):: plat(nlon)
1329     real pmea(nlon)
1330     REAL pstd(nlon)
1331     REAL ppic(nlon)
1332     REAL pulow(nlon),pvlow(nlon),pustr(nlon),pvstr(nlon)
1333     REAL t(nlon,nlev), u(nlon,nlev), v(nlon,nlev)
1334     REAL d_t(nlon,nlev), d_u(nlon,nlev), d_v(nlon,nlev)
1335     c
1336     INTEGER i, k, ktest(nlon)
1337     c
1338     c Variables locales:
1339     c
1340     REAL zgeom(klon,klev)
1341     REAL pdtdt(klon,klev), pdudt(klon,klev), pdvdt(klon,klev)
1342     REAL pt(klon,klev), pu(klon,klev), pv(klon,klev)
1343     REAL papmf(klon,klev),papmh(klon,klev+1)
1344     c
1345     c initialiser les variables de sortie (pour securite)
1346     c
1347     DO i = 1,klon
1348     pulow(i) = 0.0
1349     pvlow(i) = 0.0
1350     pustr(i) = 0.0
1351     pvstr(i) = 0.0
1352     ENDDO
1353     DO k = 1, klev
1354     DO i = 1, klon
1355     d_t(i,k) = 0.0
1356     d_u(i,k) = 0.0
1357     d_v(i,k) = 0.0
1358     pdudt(i,k)=0.0
1359     pdvdt(i,k)=0.0
1360     pdtdt(i,k)=0.0
1361     ENDDO
1362     ENDDO
1363     c
1364     c preparer les variables d'entree (attention: l'ordre des niveaux
1365     c verticaux augmente du haut vers le bas)
1366     c
1367     DO k = 1, klev
1368     DO i = 1, klon
1369     pt(i,k) = t(i,klev-k+1)
1370     pu(i,k) = u(i,klev-k+1)
1371     pv(i,k) = v(i,klev-k+1)
1372     papmf(i,k) = pplay(i,klev-k+1)
1373     ENDDO
1374     ENDDO
1375     DO k = 1, klev+1
1376     DO i = 1, klon
1377     papmh(i,k) = paprs(i,klev-k+2)
1378     ENDDO
1379     ENDDO
1380     DO i = 1, klon
1381     zgeom(i,klev) = RD * pt(i,klev)
1382     . * LOG(papmh(i,klev+1)/papmf(i,klev))
1383     ENDDO
1384     DO k = klev-1, 1, -1
1385     DO i = 1, klon
1386     zgeom(i,k) = zgeom(i,k+1) + RD * (pt(i,k)+pt(i,k+1))/2.0
1387     . * LOG(papmf(i,k+1)/papmf(i,k))
1388     ENDDO
1389     ENDDO
1390     c
1391     c appeler la routine principale
1392     c
1393     CALL OROLIFT(klon,klev,ktest,
1394     . dtime,
1395     . papmh, zgeom,
1396     . pt, pu, pv,
1397     . plat,pmea, pstd, ppic,
1398     . pulow,pvlow,
1399     . pdudt,pdvdt,pdtdt)
1400     C
1401     DO k = 1, klev
1402     DO i = 1, klon
1403     d_u(i,klev+1-k) = dtime*pdudt(i,k)
1404     d_v(i,klev+1-k) = dtime*pdvdt(i,k)
1405     d_t(i,klev+1-k) = dtime*pdtdt(i,k)
1406     pustr(i) = pustr(i)
1407     cIM BUG . +RG*pdudt(i,k)*(papmh(i,k+1)-papmh(i,k))
1408     . +pdudt(i,k)*(papmh(i,k+1)-papmh(i,k))/RG
1409     pvstr(i) = pvstr(i)
1410     cIM BUG . +RG*pdvdt(i,k)*(papmh(i,k+1)-papmh(i,k))
1411     . +pdvdt(i,k)*(papmh(i,k+1)-papmh(i,k))/RG
1412     ENDDO
1413     ENDDO
1414     c
1415     RETURN
1416     END
1417     SUBROUTINE OROLIFT( NLON,NLEV
1418     I , KTEST
1419     R , PTSPHY
1420     R , PAPHM1,PGEOM1,PTM1,PUM1,PVM1
1421     R , PLAT
1422     R , PMEA, PVAROR, ppic
1423     C OUTPUTS
1424     R , PULOW,PVLOW
1425     R , PVOM,PVOL,PTE )
1426    
1427     C
1428     C**** *OROLIFT: SIMULATE THE GEOSTROPHIC LIFT.
1429     C
1430     C PURPOSE.
1431     C --------
1432     C
1433     C** INTERFACE.
1434     C ----------
1435     C CALLED FROM *lift_noro
1436     C ----------
1437     C
1438     C AUTHOR.
1439     C -------
1440     C F.LOTT LMD 22/11/95
1441     C
1442     use dimens_m
1443     use dimphy
1444     use YOMCST
1445 guez 10 use yoegwd
1446 guez 3 implicit none
1447     C
1448     C
1449     C-----------------------------------------------------------------------
1450     C
1451     C* 0.1 ARGUMENTS
1452     C ---------
1453     C
1454     C
1455     integer nlon, nlev
1456     REAL PTE(NLON,NLEV),
1457     * PVOL(NLON,NLEV),
1458     * PVOM(NLON,NLEV),
1459     * PULOW(NLON),
1460     * PVLOW(NLON)
1461     REAL PUM1(NLON,NLEV),
1462     * PVM1(NLON,NLEV),
1463     * PTM1(NLON,NLEV)
1464     real, intent(in):: PLAT(NLON)
1465     real PMEA(NLON),
1466     * PVAROR(NLON),
1467     * ppic(NLON),
1468     * PGEOM1(NLON,NLEV),
1469     * PAPHM1(NLON,NLEV+1)
1470     C
1471     INTEGER KTEST(NLON)
1472     real ptsphy
1473     C-----------------------------------------------------------------------
1474     C
1475     C* 0.2 LOCAL ARRAYS
1476     C ------------
1477     logical lifthigh
1478     integer klevm1, jl, ilevh, jk
1479     real zcons1, ztmst, zrtmst,zpi, zhgeo
1480     real zdelp, zslow, zsqua, zscav, zbet
1481     INTEGER
1482     * IKNUB(klon),
1483     * IKNUL(klon)
1484     LOGICAL LL1(KLON,KLEV+1)
1485     C
1486     REAL ZTAU(KLON,KLEV+1),
1487     * ZTAV(KLON,KLEV+1),
1488     * ZRHO(KLON,KLEV+1)
1489     REAL ZDUDT(KLON),
1490     * ZDVDT(KLON)
1491     REAL ZHCRIT(KLON,KLEV)
1492     C-----------------------------------------------------------------------
1493     C
1494     C* 1.1 INITIALIZATIONS
1495     C ---------------
1496    
1497     LIFTHIGH=.FALSE.
1498    
1499     IF(NLON.NE.KLON.OR.NLEV.NE.KLEV)STOP
1500     ZCONS1=1./RD
1501     KLEVM1=KLEV-1
1502     ZTMST=PTSPHY
1503     ZRTMST=1./ZTMST
1504     ZPI=ACOS(-1.)
1505     C
1506     DO 1001 JL=1,klon
1507     ZRHO(JL,KLEV+1) =0.0
1508     PULOW(JL) =0.0
1509     PVLOW(JL) =0.0
1510     iknub(JL) =klev
1511     iknul(JL) =klev
1512     ilevh=klev/3
1513     ll1(jl,klev+1)=.false.
1514     DO 1000 JK=1,KLEV
1515     PVOM(JL,JK)=0.0
1516     PVOL(JL,JK)=0.0
1517     PTE (JL,JK)=0.0
1518     1000 CONTINUE
1519     1001 CONTINUE
1520    
1521     C
1522     C* 2.1 DEFINE LOW LEVEL WIND, PROJECT WINDS IN PLANE OF
1523     C* LOW LEVEL WIND, DETERMINE SECTOR IN WHICH TO TAKE
1524     C* THE VARIANCE AND SET INDICATOR FOR CRITICAL LEVELS.
1525     C
1526     C
1527     C
1528     DO 2006 JK=KLEV,1,-1
1529     DO 2007 JL=1,klon
1530     IF(KTEST(JL).EQ.1) THEN
1531     ZHCRIT(JL,JK)=amax1(Ppic(JL)-pmea(JL),100.)
1532     ZHGEO=PGEOM1(JL,JK)/RG
1533     ll1(JL,JK)=(ZHGEO.GT.ZHCRIT(JL,JK))
1534     IF(ll1(JL,JK).neqv.ll1(JL,JK+1)) THEN
1535     iknub(JL)=JK
1536     ENDIF
1537     ENDIF
1538     2007 CONTINUE
1539     2006 CONTINUE
1540     C
1541     do 2010 jl=1,klon
1542     IF(KTEST(JL).EQ.1) THEN
1543     iknub(jl)=max(iknub(jl),klev/2)
1544     iknul(jl)=max(iknul(jl),2*klev/3)
1545     if(iknub(jl).gt.nktopg) iknub(jl)=nktopg
1546     if(iknub(jl).eq.nktopg) iknul(jl)=klev
1547     if(iknub(jl).eq.iknul(jl)) iknub(jl)=iknul(jl)-1
1548     ENDIF
1549     2010 continue
1550    
1551     C do 2011 jl=1,klon
1552     C IF(KTEST(JL).EQ.1) THEN
1553     C print *,' iknul= ',iknul(jl),' iknub=',iknub(jl)
1554     C ENDIF
1555     C2011 continue
1556    
1557     C PRINT *,' DANS OROLIFT: 2010'
1558    
1559     DO 223 JK=KLEV,2,-1
1560     DO 222 JL=1,klon
1561     ZRHO(JL,JK)=2.*PAPHM1(JL,JK)*ZCONS1/(PTM1(JL,JK)+PTM1(JL,JK-1))
1562     222 CONTINUE
1563     223 CONTINUE
1564     C PRINT *,' DANS OROLIFT: 223'
1565    
1566     C********************************************************************
1567     C
1568     C* DEFINE LOW LEVEL FLOW
1569     C -------------------
1570     DO 2115 JK=klev,1,-1
1571     DO 2116 JL=1,klon
1572     IF(KTEST(JL).EQ.1) THEN
1573     if(jk.ge.iknub(jl).and.jk.le.iknul(jl)) then
1574     pulow(JL)=pulow(JL)+PUM1(JL,JK)*(PAPHM1(JL,JK+1)-PAPHM1(JL,JK))
1575     pvlow(JL)=pvlow(JL)+PVM1(JL,JK)*(PAPHM1(JL,JK+1)-PAPHM1(JL,JK))
1576     zrho(JL,klev+1)=zrho(JL,klev+1)
1577     * +zrho(JL,JK)*(PAPHM1(JL,JK+1)-PAPHM1(JL,JK))
1578     end if
1579     ENDIF
1580     2116 CONTINUE
1581     2115 CONTINUE
1582     DO 2110 JL=1,klon
1583     IF(KTEST(JL).EQ.1) THEN
1584     pulow(JL)=pulow(JL)/(PAPHM1(JL,iknul(jl)+1)-PAPHM1(JL,iknub(jl)))
1585     pvlow(JL)=pvlow(JL)/(PAPHM1(JL,iknul(jl)+1)-PAPHM1(JL,iknub(jl)))
1586     zrho(JL,klev+1)=zrho(JL,klev+1)
1587     * /(PAPHM1(JL,iknul(jl)+1)-PAPHM1(JL,iknub(jl)))
1588     ENDIF
1589     2110 CONTINUE
1590    
1591    
1592     200 CONTINUE
1593    
1594     C***********************************************************
1595     C
1596     C* 3. COMPUTE MOUNTAIN LIFT
1597     C
1598     300 CONTINUE
1599     C
1600     DO 301 JL=1,klon
1601     IF(KTEST(JL).EQ.1) THEN
1602     ZTAU(JL,KLEV+1)= - GKLIFT*ZRHO(JL,KLEV+1)*2.*ROMEGA*
1603     C * (2*PVAROR(JL)+PMEA(JL))*
1604     * 2*PVAROR(JL)*
1605     * SIN(ZPI/180.*PLAT(JL))*PVLOW(JL)
1606     ZTAV(JL,KLEV+1)= GKLIFT*ZRHO(JL,KLEV+1)*2.*ROMEGA*
1607     C * (2*PVAROR(JL)+PMEA(JL))*
1608     * 2*PVAROR(JL)*
1609     * SIN(ZPI/180.*PLAT(JL))*PULOW(JL)
1610     ELSE
1611     ZTAU(JL,KLEV+1)=0.0
1612     ZTAV(JL,KLEV+1)=0.0
1613     ENDIF
1614     301 CONTINUE
1615    
1616     C
1617     C* 4. COMPUTE LIFT PROFILE
1618     C* --------------------
1619     C
1620    
1621     400 CONTINUE
1622    
1623     DO 401 JK=1,KLEV
1624     DO 401 JL=1,klon
1625     IF(KTEST(JL).EQ.1) THEN
1626     ZTAU(JL,JK)=ZTAU(JL,KLEV+1)*PAPHM1(JL,JK)/PAPHM1(JL,KLEV+1)
1627     ZTAV(JL,JK)=ZTAV(JL,KLEV+1)*PAPHM1(JL,JK)/PAPHM1(JL,KLEV+1)
1628     ELSE
1629     ZTAU(JL,JK)=0.0
1630     ZTAV(JL,JK)=0.0
1631     ENDIF
1632     401 CONTINUE
1633     C
1634     C
1635     C* 5. COMPUTE TENDENCIES.
1636     C* -------------------
1637     IF(LIFTHIGH)THEN
1638     C
1639     500 CONTINUE
1640     C PRINT *,' DANS OROLIFT: 500'
1641     C
1642     C EXPLICIT SOLUTION AT ALL LEVELS
1643     C
1644     DO 524 JK=1,klev
1645     DO 523 JL=1,KLON
1646     IF(KTEST(JL).EQ.1) THEN
1647     ZDELP=PAPHM1(JL,JK+1)-PAPHM1(JL,JK)
1648     ZDUDT(JL)=-RG*(ZTAU(JL,JK+1)-ZTAU(JL,JK))/ZDELP
1649     ZDVDT(JL)=-RG*(ZTAV(JL,JK+1)-ZTAV(JL,JK))/ZDELP
1650     ENDIF
1651     523 CONTINUE
1652     524 CONTINUE
1653     C
1654     C PROJECT PERPENDICULARLY TO U NOT TO DESTROY ENERGY
1655     C
1656     DO 530 JK=1,klev
1657     DO 530 JL=1,KLON
1658     IF(KTEST(JL).EQ.1) THEN
1659    
1660     ZSLOW=SQRT(PULOW(JL)**2+PVLOW(JL)**2)
1661     ZSQUA=AMAX1(SQRT(PUM1(JL,JK)**2+PVM1(JL,JK)**2),GVSEC)
1662     ZSCAV=-ZDUDT(JL)*PVM1(JL,JK)+ZDVDT(JL)*PUM1(JL,JK)
1663     IF(ZSQUA.GT.GVSEC)THEN
1664     PVOM(JL,JK)=-ZSCAV*PVM1(JL,JK)/ZSQUA**2
1665     PVOL(JL,JK)= ZSCAV*PUM1(JL,JK)/ZSQUA**2
1666     ELSE
1667     PVOM(JL,JK)=0.0
1668     PVOL(JL,JK)=0.0
1669     ENDIF
1670     ZSQUA=SQRT(PUM1(JL,JK)**2+PUM1(JL,JK)**2)
1671     IF(ZSQUA.LT.ZSLOW)THEN
1672     PVOM(JL,JK)=ZSQUA/ZSLOW*PVOM(JL,JK)
1673     PVOL(JL,JK)=ZSQUA/ZSLOW*PVOL(JL,JK)
1674     ENDIF
1675    
1676     ENDIF
1677     530 CONTINUE
1678     C
1679     C 6. LOW LEVEL LIFT, SEMI IMPLICIT:
1680     C ----------------------------------
1681    
1682     ELSE
1683    
1684     DO 601 JL=1,KLON
1685     IF(KTEST(JL).EQ.1) THEN
1686     DO JK=KLEV,IKNUB(JL),-1
1687     ZBET=GKLIFT*2.*ROMEGA*SIN(ZPI/180.*PLAT(JL))*ztmst*
1688     * (PGEOM1(JL,IKNUB(JL)-1)-PGEOM1(JL, JK))/
1689     * (PGEOM1(JL,IKNUB(JL)-1)-PGEOM1(JL,KLEV))
1690     ZDUDT(JL)=-PUM1(JL,JK)/ztmst/(1+ZBET**2)
1691     ZDVDT(JL)=-PVM1(JL,JK)/ztmst/(1+ZBET**2)
1692     PVOM(JL,JK)= ZBET**2*ZDUDT(JL) - ZBET *ZDVDT(JL)
1693     PVOL(JL,JK)= ZBET *ZDUDT(JL) + ZBET**2*ZDVDT(JL)
1694     ENDDO
1695     ENDIF
1696     601 CONTINUE
1697    
1698     ENDIF
1699    
1700     RETURN
1701     END
1702     SUBROUTINE SUGWD(NLON,NLEV,paprs,pplay)
1703     C
1704     C**** *SUGWD* INITIALIZE COMMON YOEGWD CONTROLLING GRAVITY WAVE DRAG
1705     C
1706     C PURPOSE.
1707     C --------
1708     C INITIALIZE YOEGWD, THE COMMON THAT CONTROLS THE
1709     C GRAVITY WAVE DRAG PARAMETRIZATION.
1710     C
1711     C** INTERFACE.
1712     C ----------
1713     C CALL *SUGWD* FROM *SUPHEC*
1714     C ----- ------
1715     C
1716     C EXPLICIT ARGUMENTS :
1717     C --------------------
1718     C PSIG : VERTICAL COORDINATE TABLE
1719     C NLEV : NUMBER OF MODEL LEVELS
1720     C
1721     C IMPLICIT ARGUMENTS :
1722     C --------------------
1723     C COMMON YOEGWD
1724     C
1725     C METHOD.
1726     C -------
1727     C SEE DOCUMENTATION
1728     C
1729     C EXTERNALS.
1730     C ----------
1731     C NONE
1732     C
1733     C REFERENCE.
1734     C ----------
1735     C ECMWF Research Department documentation of the IFS
1736     C
1737     C AUTHOR.
1738     C -------
1739     C MARTIN MILLER *ECMWF*
1740     C
1741     C MODIFICATIONS.
1742     C --------------
1743     C ORIGINAL : 90-01-01
1744     C ------------------------------------------------------------------
1745 guez 10 use yoegwd
1746 guez 3 implicit none
1747     C
1748     C -----------------------------------------------------------------
1749     C ----------------------------------------------------------------
1750     C
1751     integer nlon,nlev, jk
1752     REAL, intent(in):: paprs(nlon,nlev+1)
1753 guez 10 REAL, intent(in):: pplay(nlon,nlev)
1754 guez 3 real zpr,zstra,zsigt,zpm1r
1755     C
1756     C* 1. SET THE VALUES OF THE PARAMETERS
1757     C --------------------------------
1758     C
1759     100 CONTINUE
1760     C
1761     PRINT *,' DANS SUGWD NLEV=',NLEV
1762     GHMAX=10000.
1763     C
1764     ZPR=100000.
1765     ZSTRA=0.1
1766     ZSIGT=0.94
1767     cold ZPR=80000.
1768     cold ZSIGT=0.85
1769     C
1770     DO 110 JK=1,NLEV
1771     ZPM1R=pplay(nlon/2,jk)/paprs(nlon/2,1)
1772     IF(ZPM1R.GE.ZSIGT)THEN
1773     nktopg=JK
1774     ENDIF
1775     ZPM1R=pplay(nlon/2,jk)/paprs(nlon/2,1)
1776     IF(ZPM1R.GE.ZSTRA)THEN
1777     NSTRA=JK
1778     ENDIF
1779     110 CONTINUE
1780     c
1781     c inversion car dans orodrag on compte les niveaux a l'envers
1782     nktopg=nlev-nktopg+1
1783     nstra=nlev-nstra
1784     print *,' DANS SUGWD nktopg=', nktopg
1785     print *,' DANS SUGWD nstra=', nstra
1786     C
1787     GSIGCR=0.80
1788     C
1789     GKDRAG=0.2
1790     GRAHILO=1.
1791     GRCRIT=0.01
1792     GFRCRIT=1.0
1793     GKWAKE=0.50
1794     C
1795     GKLIFT=0.50
1796     GVCRIT =0.0
1797     C
1798     C
1799     C ----------------------------------------------------------------
1800     C
1801     C* 2. SET VALUES OF SECURITY PARAMETERS
1802     C ---------------------------------
1803     C
1804     200 CONTINUE
1805     C
1806     GVSEC=0.10
1807     GSSEC=1.E-12
1808     C
1809     GTSEC=1.E-07
1810     C
1811     C ----------------------------------------------------------------
1812     C
1813     RETURN
1814     END
  ViewVC Help
Powered by ViewVC 1.1.21