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