source: lmdz_wrf/branches/LMDZ_WRFmeas/WRFV3/lmdz/orografi.F90 @ 27

Last change on this file since 27 was 1, checked in by lfita, 10 years ago
  • -- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

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