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

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

Parent Directory Parent Directory | Revision Log Revision Log

Revision 13 - (hide annotations)
Fri Jul 25 19:59:34 2008 UTC (15 years, 10 months ago) by guez
File size: 48890 byte(s) 
-- Minor change of behaviour:

"etat0" does not compute "rugsrel" nor "radpas". Deleted arguments
"radpas" and "rugsrel" of "phyredem". Deleted argument "rugsrel" of
"phyetat0". "startphy.nc" does not contain the variable "RUGSREL". In
"physiq", "rugoro" is set to 0 if not "ok_orodr". The whole program
"etat0_lim" does not use "clesphys2".

-- Minor modification of input/output:

Created subroutine "read_clesphys2". Variables of "clesphys2" are read
in "read_clesphys2" instead of "conf_gcm". "printflag" does not print
variables of "clesphys2".

-- Should not change any result at run time:

References to module "numer_rec" instead of individual modules of
"Numer_rec_Lionel".

Deleted argument "clesphy0" of "calfis", "physiq", "conf_gcm",
"leapfrog", "phyetat0". Deleted variable "clesphy0" in
"gcm". "phyetat0" does not modify variables of "clesphys2".

The program unit "gcm" does not modify "itau_phy".

Added some "intent" attributes.

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