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

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

Parent Directory Parent Directory | Revision Log Revision Log

Revision 12 - (hide annotations)
Mon Jul 21 16:05:07 2008 UTC (15 years, 11 months ago) by guez
File size: 48826 byte(s) 
-- Minor modification of input/output:

Created procedure "read_logic". Variables of module "logic" are read
by "read_logic" instead of "conf_gcm". Variable "offline" of module
"conf_gcm" is read from namelist instead of "*.def".

Deleted arguments "dtime", "co2_ppm_etat0", "solaire_etat0",
"tabcntr0" and local variables "radpas", "tab_cntrl" of
"phyetat0". "phyetat0" does not read "controle" in "startphy.nc" any
longer. "phyetat0" now reads global attribute "itau_phy" from
"startphy.nc". "phyredem" does not create variable "controle" in
"startphy.nc" any longer. "phyredem" now writes global attribute
"itau_phy" of "startphy.nc". Deleted argument "tabcntr0" of
"printflag". Removed diagnostic messages written by "printflag" for
comparison of the variable "controle" of "startphy.nc" and the
variables read from "*.def" or namelist input.

-- Removing unwanted functionality:

Removed variable "lunout" from module "iniprint", replaced everywhere
by standard output.

Removed case "ocean == 'couple'" in "clmain", "interfsurf_hq" and
"physiq". Removed procedure "interfoce_cpl".

-- Should not change anything at run time:

Automated creation of graphs in documentation. More documentation on
input files.

Converted Fortran files to free format: "phyredem.f90", "printflag.f90".

Split module "clesphy" into "clesphys" and "clesphys2".

Removed variables "conser", "leapf", "forward", "apphys", "apdiss" and
"statcl" from module "logic". Added arguments "conser" to "advect",
"leapf" to "integrd". Added local variables "forward", "leapf",
"apphys", "conser", "apdiss" in "leapfrog".

Added intent attributes.

Deleted arguments "dtime" of "phyredem", "pdtime" of "flxdtdq", "sh"
of "phytrac", "dt" of "yamada".

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