source: LMDZ4/trunk/libf/phylmd/orografi_strato.F @ 1176

Last change on this file since 1176 was 1001, checked in by Laurent Fairhead, 16 years ago
  • Modifs sur le parallelisme: masquage dans la physique
  • Inclusion strato
  • mise en coherence etat0
  • le mode offline fonctionne maintenant en parallele,
  • les fichiers de la dynamiques sont correctement sortis et peuvent etre reconstruit avec rebuild
  • la version parallele de la dynamique peut s'executer sans MPI (sur 1 proc)
  • L'OPENMP fonctionne maintenant sans la parallelisation MPI.

YM
LF
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 57.3 KB
RevLine 
[1001]1      SUBROUTINE drag_noro_strato (nlon,nlev,dtime,paprs,pplay,
2     e                   pmea,pstd, psig, pgam, pthe,ppic,pval,
3     e                   kgwd,kdx,ktest,
4     e                   t, u, v,
5     s                   pulow, pvlow, pustr, pvstr,
6     s                   d_t, d_u, d_v)
7c
8      USE dimphy
9      IMPLICIT none
10c======================================================================
11c Auteur(s): F.Lott (LMD/CNRS) date: 19950201
12c Object: Mountain drag interface. Made necessary because:
13C 1. in the LMD-GCM Layers are from bottom to top,
14C    contrary to most European GCM.
15c 2. the altitude above ground of each model layers
16c    needs to be known (variable zgeom)
17c======================================================================
18c Explicit Arguments:
19c ==================
20c nlon----input-I-Total number of horizontal points that get into physics
21c nlev----input-I-Number of vertical levels
22c dtime---input-R-Time-step (s)
23c paprs---input-R-Pressure in semi layers    (Pa)
24c pplay---input-R-Pressure model-layers      (Pa)
25c t-------input-R-temperature (K)
26c u-------input-R-Horizontal wind (m/s)
27c v-------input-R-Meridional wind (m/s)
28c pmea----input-R-Mean Orography (m)
29C pstd----input-R-SSO standard deviation (m)
30c psig----input-R-SSO slope
31c pgam----input-R-SSO Anisotropy
32c pthe----input-R-SSO Angle
33c ppic----input-R-SSO Peacks elevation (m)
34c pval----input-R-SSO Valleys elevation (m)
35c
36c kgwd- -input-I: Total nb of points where the orography schemes are active
37c ktest--input-I: Flags to indicate active points
38c kdx----input-I: Locate the physical location of an active point.
39
40c pulow, pvlow -output-R: Low-level wind
41c pustr, pvstr -output-R: Surface stress due to SSO drag      (Pa)
42c
43c d_t-----output-R: T increment           
44c d_u-----output-R: U increment             
45c d_v-----output-R: V increment             
46c
47c Implicit Arguments:
48c ===================
49c
50c iim--common-I: Number of longitude intervals
51c jjm--common-I: Number of latitude intervals
52c klon-common-I: Number of points seen by the physics
53c                (iim+1)*(jjm+1) for instance
54c klev-common-I: Number of vertical layers
55c======================================================================
56c Local Variables:
57c ================
58c
59c zgeom-----R: Altitude of layer above ground
60c pt, pu, pv --R: t u v from top to bottom
61c pdtdt, pdudt, pdvdt --R: t u v tendencies (from top to bottom)
62c papmf: pressure at model layer (from top to bottom)
63c papmh: pressure at model 1/2 layer (from top to bottom)
64c
65c======================================================================
66cym#include "dimensions.h"
67cym#include "dimphy.h"
68#include "YOMCST.h"
69#include "YOEGWD.h"
70c
71c  ARGUMENTS
72c
73      INTEGER nlon,nlev
74      REAL dtime
75      REAL paprs(nlon,nlev+1)
76      REAL pplay(nlon,nlev)
77      REAL pmea(nlon),pstd(nlon),psig(nlon),pgam(nlon),pthe(nlon)
78      REAL ppic(nlon),pval(nlon)
79      REAL pulow(nlon),pvlow(nlon),pustr(nlon),pvstr(nlon)
80      REAL t(nlon,nlev), u(nlon,nlev), v(nlon,nlev)
81      REAL d_t(nlon,nlev), d_u(nlon,nlev), d_v(nlon,nlev)
82c
83      INTEGER i, k, kgwd,  kdx(nlon), ktest(nlon)
84c
85c LOCAL VARIABLES:
86c
87      REAL zgeom(klon,klev)
88      REAL pdtdt(klon,klev), pdudt(klon,klev), pdvdt(klon,klev)
89      REAL pt(klon,klev), pu(klon,klev), pv(klon,klev)
90      REAL papmf(klon,klev),papmh(klon,klev+1)
91c
92c INITIALIZE OUTPUT VARIABLES
93c
94      DO i = 1,klon
95         pulow(i) = 0.0
96         pvlow(i) = 0.0
97         pustr(i) = 0.0
98         pvstr(i) = 0.0
99      ENDDO
100      DO k = 1, klev
101      DO i = 1, klon
102         d_t(i,k) = 0.0
103         d_u(i,k) = 0.0
104         d_v(i,k) = 0.0
105         pdudt(i,k)=0.0
106         pdvdt(i,k)=0.0
107         pdtdt(i,k)=0.0
108      ENDDO
109      ENDDO
110c
111c PREPARE INPUT VARIABLES FOR ORODRAG (i.e., ORDERED FROM TOP TO BOTTOM)
112C CALCULATE LAYERS HEIGHT ABOVE GROUND)
113c
114      DO k = 1, klev
115      DO i = 1, klon
116         pt(i,k) = t(i,klev-k+1)
117         pu(i,k) = u(i,klev-k+1)
118         pv(i,k) = v(i,klev-k+1)
119         papmf(i,k) = pplay(i,klev-k+1)
120      ENDDO
121      ENDDO
122      DO k = 1, klev+1
123      DO i = 1, klon
124         papmh(i,k) = paprs(i,klev-k+2)
125      ENDDO
126      ENDDO
127      DO i = 1, klon
128         zgeom(i,klev) = RD * pt(i,klev)
129     .                  * LOG(papmh(i,klev+1)/papmf(i,klev))
130      ENDDO
131      DO k = klev-1, 1, -1
132      DO i = 1, klon
133         zgeom(i,k) = zgeom(i,k+1) + RD * (pt(i,k)+pt(i,k+1))/2.0
134     .               * LOG(papmf(i,k+1)/papmf(i,k))
135      ENDDO
136      ENDDO
137c
138c CALL SSO DRAG ROUTINES       
139c
140      CALL orodrag_strato(klon,klev,kgwd,kdx,ktest,
141     .            dtime,
142     .            papmh, papmf, zgeom,
143     .            pt, pu, pv,
144     .            pmea, pstd, psig, pgam, pthe, ppic,pval,
145     .            pulow,pvlow,
146     .            pdudt,pdvdt,pdtdt)
147C
148C COMPUTE INCREMENTS AND STRESS FROM TENDENCIES
149
150      DO k = 1, klev
151      DO i = 1, klon
152         d_u(i,klev+1-k) = dtime*pdudt(i,k)
153         d_v(i,klev+1-k) = dtime*pdvdt(i,k)
154         d_t(i,klev+1-k) = dtime*pdtdt(i,k)
155         pustr(i)        = pustr(i)
156     .                    +pdudt(i,k)*(papmh(i,k+1)-papmh(i,k))/rg
157         pvstr(i)        = pvstr(i)
158     .                    +pdvdt(i,k)*(papmh(i,k+1)-papmh(i,k))/rg
159      ENDDO
160      ENDDO
161c
162      RETURN
163      END
164
165      SUBROUTINE orodrag_strato( nlon,nlev
166     i                 , kgwd,  kdx, ktest
167     r                 , ptsphy
168     r                 , paphm1,papm1,pgeom1,ptm1,pum1,pvm1
169     r                 , pmea, pstd, psig, pgam, pthe, ppic, pval
170c outputs
171     r                 , pulow,pvlow
172     r                 , pvom,pvol,pte )
173     
174      USE dimphy
175      IMPLICIT NONE
176c
177c
178c**** *orodrag* - does the SSO drag  parametrization.
179c
180c     purpose.
181c     --------
182c
183c     this routine computes the physical tendencies of the
184c     prognostic variables u,v  and t due to  vertical transports by
185c     subgridscale orographically excited gravity waves, and to
186c     low level blocked flow drag.
187c
188c**   interface.
189c     ----------
190c          called from *drag_noro*.
191c
192c          the routine takes its input from the long-term storage:
193c          u,v,t and p at t-1.
194c
195c        explicit arguments :
196c        --------------------
197c     ==== inputs ===
198c nlon----input-I-Total number of horizontal points that get into physics
199c nlev----input-I-Number of vertical levels
200c
201c kgwd- -input-I: Total nb of points where the orography schemes are active
202c ktest--input-I: Flags to indicate active points
203c kdx----input-I: Locate the physical location of an active point.
204c ptsphy--input-R-Time-step (s)
205c paphm1--input-R: pressure at model 1/2 layer
206c papm1---input-R: pressure at model layer
207c pgeom1--input-R: Altitude of layer above ground
208c ptm1, pum1, pvm1--R-: t, u and v
209c pmea----input-R-Mean Orography (m)
210C pstd----input-R-SSO standard deviation (m)
211c psig----input-R-SSO slope
212c pgam----input-R-SSO Anisotropy
213c pthe----input-R-SSO Angle
214c ppic----input-R-SSO Peacks elevation (m)
215c pval----input-R-SSO Valleys elevation (m)
216
217      integer nlon,nlev,kgwd
218      real ptsphy
219
220c     ==== outputs ===
221c pulow, pvlow -output-R: Low-level wind
222c
223c pte -----output-R: T tendency
224c pvom-----output-R: U tendency
225c pvol-----output-R: V tendency
226c
227c
228c Implicit Arguments:
229c ===================
230c
231c klon-common-I: Number of points seen by the physics
232c klev-common-I: Number of vertical layers
233c
234c     method.
235c     -------
236c
237c     externals.
238c     ----------
239      integer ismin, ismax
240      external ismin, ismax
241c
242c     reference.
243c     ----------
244c
245c     author.
246c     -------
247c     m.miller + b.ritter   e.c.m.w.f.     15/06/86.
248c
249c     f.lott + m. miller    e.c.m.w.f.     22/11/94
250c-----------------------------------------------------------------------
251c
252c
253cym#include "dimensions.h"
254cym#include "dimphy.h"
255#include "YOMCST.h"
256#include "YOEGWD.h"
257c-----------------------------------------------------------------------
258c
259c*       0.1   arguments
260c              ---------
261c
262c
263      real  pte(nlon,nlev),
264     *      pvol(nlon,nlev),
265     *      pvom(nlon,nlev),
266     *      pulow(nlon),
267     *      pvlow(nlon)
268      real  pum1(nlon,nlev),
269     *      pvm1(nlon,nlev),
270     *      ptm1(nlon,nlev),
271     *      pmea(nlon),pstd(nlon),psig(nlon),
272     *      pgam(nlon),pthe(nlon),ppic(nlon),pval(nlon),
273     *      pgeom1(nlon,nlev),
274     *      papm1(nlon,nlev),
275     *      paphm1(nlon,nlev+1)
276c
277      integer  kdx(nlon),ktest(nlon)
278c-----------------------------------------------------------------------
279c
280c*       0.2   local arrays
281c              ------------
282      integer  isect(klon),
283     *         icrit(klon),
284     *         ikcrith(klon),
285     *         ikenvh(klon),
286     *         iknu(klon),
287     *         iknu2(klon),
288     *         ikcrit(klon),
289     *         ikhlim(klon)
290c
291      real   ztau(klon,klev+1),
292     *       zstab(klon,klev+1),
293     *       zvph(klon,klev+1),
294     *       zrho(klon,klev+1),
295     *       zri(klon,klev+1),
296     *       zpsi(klon,klev+1),
297     *       zzdep(klon,klev)
298      real   zdudt(klon),
299     *       zdvdt(klon),
300     *       zdtdt(klon),
301     *       zdedt(klon),
302     *       zvidis(klon),
303     *       ztfr(klon),
304     *       znu(klon),
305     *       zd1(klon),
306     *       zd2(klon),
307     *       zdmod(klon)
308
309
310c local quantities:
311
312      integer jl,jk,ji
313      real ztmst,zdelp,ztemp,zforc,ztend,rover               
314      real zb,zc,zconb,zabsv,zzd1,ratio,zbet,zust,zvst,zdis
315   
316c
317c------------------------------------------------------------------
318c
319c*         1.    initialization
320c                --------------
321c
322c        print *,' in orodrag'
323 100  continue
324c
325c     ------------------------------------------------------------------
326c
327c*         1.1   computational constants
328c                -----------------------
329c
330 110  continue
331c
332c     ztmst=twodt
333c     if(nstep.eq.nstart) ztmst=0.5*twodt
334      ztmst=ptsphy
335c     ------------------------------------------------------------------
336c
337 120  continue
338c
339c     ------------------------------------------------------------------
340c
341c*         1.3   check whether row contains point for printing
342c                ---------------------------------------------
343c
344 130  continue
345c
346c     ------------------------------------------------------------------
347c
348c*         2.     precompute basic state variables.
349c*                ---------- ----- ----- ----------
350c*                define low level wind, project winds in plane of
351c*                low level wind, determine sector in which to take
352c*                the variance and set indicator for critical levels.
353c
354
355  200 continue
356c
357c
358c
359      call orosetup_strato
360     *     ( nlon, nlev , ktest
361     *     , ikcrit, ikcrith, icrit, isect, ikhlim, ikenvh,iknu,iknu2
362     *     , paphm1, papm1 , pum1   , pvm1 , ptm1 , pgeom1, pstd
363     *     , zrho  , zri   , zstab  , ztau , zvph , zpsi, zzdep
364     *     , pulow, pvlow
365     *     , pthe,pgam,pmea,ppic,pval,znu  ,zd1,  zd2,  zdmod )
366
367
368c
369c
370c
371c***********************************************************
372c
373c
374c*         3.      compute low level stresses using subcritical and
375c*                 supercritical forms.computes anisotropy coefficient
376c*                 as measure of orographic twodimensionality.
377c
378  300 continue
379c
380      call gwstress_strato
381     *    ( nlon  , nlev
382     *    , ikcrit, isect, ikhlim, ktest, ikcrith, icrit, ikenvh, iknu
383     *    , zrho  , zstab, zvph  , pstd,  psig, pmea, ppic, pval
384     *    , ztfr   , ztau
385     *    , pgeom1,pgam,zd1,zd2,zdmod,znu)
386
387c
388c
389c*         4.      compute stress profile including
390c                  trapped waves, wave breaking,
391c                  linear decay in stratosphere.
392c
393  400 continue
394c
395c
396
397      call gwprofil_strato
398     *       (  nlon , nlev
399     *       , kgwd   , kdx  , ktest
400     *       , ikcrit, ikcrith, icrit  , ikenvh, iknu
401     *       ,iknu2 , paphm1, zrho   , zstab , ztfr   , zvph
402     *       , zri   , ztau
403 
404     *       , zdmod , znu    , psig  , pgam , pstd , ppic , pval)
405
406c
407c*         5.      Compute tendencies from waves stress profile.
408c                  Compute low level blocked flow drag.
409c*                 --------------------------------------------
410c
411  500 continue
412
413     
414c
415c  explicit solution at all levels for the gravity wave
416c  implicit solution for the blocked levels
417
418      do 510 jl=kidia,kfdia
419      zvidis(jl)=0.0
420      zdudt(jl)=0.0
421      zdvdt(jl)=0.0
422      zdtdt(jl)=0.0
423  510 continue
424c
425
426      do 524 jk=1,klev
427c
428
429C  WAVE STRESS
430C-------------
431c
432c
433      do 523 ji=kidia,kfdia
434
435      if(ktest(ji).eq.1) then
436
437      zdelp=paphm1(ji,jk+1)-paphm1(ji,jk)
438      ztemp=-rg*(ztau(ji,jk+1)-ztau(ji,jk))/(zvph(ji,klev+1)*zdelp)
439
440      zdudt(ji)=(pulow(ji)*zd1(ji)-pvlow(ji)*zd2(ji))*ztemp/zdmod(ji)
441      zdvdt(ji)=(pvlow(ji)*zd1(ji)+pulow(ji)*zd2(ji))*ztemp/zdmod(ji)
442c
443c Control Overshoots
444c
445
446      if(jk.ge.nstra)then
447        rover=0.10
448        if(abs(zdudt(ji)).gt.rover*abs(pum1(ji,jk))/ztmst)
449     C    zdudt(ji)=rover*abs(pum1(ji,jk))/ztmst*
450     C              zdudt(ji)/(abs(zdudt(ji))+1.E-10)
451        if(abs(zdvdt(ji)).gt.rover*abs(pvm1(ji,jk))/ztmst)
452     C    zdvdt(ji)=rover*abs(pvm1(ji,jk))/ztmst*
453     C              zdvdt(ji)/(abs(zdvdt(ji))+1.E-10)
454      endif
455
456      rover=0.25
457      zforc=sqrt(zdudt(ji)**2+zdvdt(ji)**2)       
458      ztend=sqrt(pum1(ji,jk)**2+pvm1(ji,jk)**2)/ztmst                     
459
460      if(zforc.ge.rover*ztend)then
461        zdudt(ji)=rover*ztend/zforc*zdudt(ji)
462        zdvdt(ji)=rover*ztend/zforc*zdvdt(ji)
463      endif
464c
465c BLOCKED FLOW DRAG:
466C -----------------
467c
468      if(jk.gt.ikenvh(ji)) then
469         zb=1.0-0.18*pgam(ji)-0.04*pgam(ji)**2
470         zc=0.48*pgam(ji)+0.3*pgam(ji)**2
471         zconb=2.*ztmst*gkwake*psig(ji)/(4.*pstd(ji))
472         zabsv=sqrt(pum1(ji,jk)**2+pvm1(ji,jk)**2)/2.
473         zzd1=zb*cos(zpsi(ji,jk))**2+zc*sin(zpsi(ji,jk))**2
474         ratio=(cos(zpsi(ji,jk))**2+pgam(ji)*sin(zpsi(ji,jk))**2)/
475     *   (pgam(ji)*cos(zpsi(ji,jk))**2+sin(zpsi(ji,jk))**2)
476         zbet=max(0.,2.-1./ratio)*zconb*zzdep(ji,jk)*zzd1*zabsv
477c
478c OPPOSED TO THE WIND
479c
480         zdudt(ji)=-pum1(ji,jk)/ztmst
481         zdvdt(ji)=-pvm1(ji,jk)/ztmst
482c
483c PERPENDICULAR TO THE SSO MAIN AXIS:
484C                           
485cmod     zdudt(ji)=-(pum1(ji,jk)*cos(pthe(ji)*rpi/180.)
486cmod *              +pvm1(ji,jk)*sin(pthe(ji)*rpi/180.))
487cmod *              *cos(pthe(ji)*rpi/180.)/ztmst
488cmod     zdvdt(ji)=-(pum1(ji,jk)*cos(pthe(ji)*rpi/180.)
489cmod *              +pvm1(ji,jk)*sin(pthe(ji)*rpi/180.))
490cmod *              *sin(pthe(ji)*rpi/180.)/ztmst
491C
492         zdudt(ji)=zdudt(ji)*(zbet/(1.+zbet))
493         zdvdt(ji)=zdvdt(ji)*(zbet/(1.+zbet))
494      end if
495      pvom(ji,jk)=zdudt(ji)
496      pvol(ji,jk)=zdvdt(ji)
497      zust=pum1(ji,jk)+ztmst*zdudt(ji)
498      zvst=pvm1(ji,jk)+ztmst*zdvdt(ji)
499      zdis=0.5*(pum1(ji,jk)**2+pvm1(ji,jk)**2-zust**2-zvst**2)
500      zdedt(ji)=zdis/ztmst
501      zvidis(ji)=zvidis(ji)+zdis*zdelp
502      zdtdt(ji)=zdedt(ji)/rcpd
503c
504c  NO TENDENCIES ON TEMPERATURE .....
505c
506c  Instead of, pte(ji,jk)=zdtdt(ji), due to mechanical dissipation
507c
508      pte(ji,jk)=0.0
509
510      endif
511
512  523 continue
513  524 continue
514c
515c
516  501 continue
517
518      return
519      end
520      SUBROUTINE orosetup_strato
521     *         ( nlon   , nlev  , ktest
522     *         , kkcrit, kkcrith, kcrit, ksect , kkhlim
523     *         , kkenvh, kknu  , kknu2
524     *         , paphm1, papm1 , pum1   , pvm1 , ptm1  , pgeom1, pstd
525     *         , prho  , pri   , pstab  , ptau , pvph  ,ppsi, pzdep
526     *         , pulow , pvlow 
527     *         , ptheta, pgam, pmea, ppic, pval
528     *         , pnu  ,  pd1  ,  pd2  ,pdmod  )
529C
530c**** *gwsetup*
531c
532c     purpose.
533c     --------
534c     SET-UP THE ESSENTIAL PARAMETERS OF THE SSO DRAG SCHEME:
535C     DEPTH OF LOW WBLOCKED LAYER, LOW-LEVEL FLOW, BACKGROUND
536C     STRATIFICATION.....
537c
538c**   interface.
539c     ----------
540c          from *orodrag*
541c
542c        explicit arguments :
543c        --------------------
544c     ==== inputs ===
545c
546c nlon----input-I-Total number of horizontal points that get into physics
547c nlev----input-I-Number of vertical levels
548c ktest--input-I: Flags to indicate active points
549c
550c ptsphy--input-R-Time-step (s)
551c paphm1--input-R: pressure at model 1/2 layer
552c papm1---input-R: pressure at model layer
553c pgeom1--input-R: Altitude of layer above ground
554c ptm1, pum1, pvm1--R-: t, u and v
555c pmea----input-R-Mean Orography (m)
556C pstd----input-R-SSO standard deviation (m)
557c psig----input-R-SSO slope
558c pgam----input-R-SSO Anisotropy
559c pthe----input-R-SSO Angle
560c ppic----input-R-SSO Peacks elevation (m)
561c pval----input-R-SSO Valleys elevation (m)
562
563c     ==== outputs ===
564c pulow, pvlow -output-R: Low-level wind
565c kkcrit----I-: Security value for top of low level flow
566c kcrit-----I-: Critical level
567c ksect-----I-: Not used
568c kkhlim----I-: Not used
569c kkenvh----I-: Top of blocked flow layer
570c kknu------I-: Layer that sees mountain peacks
571c kknu2-----I-: Layer that sees mountain peacks above mountain mean
572c kknub-----I-: Layer that sees mountain mean above valleys
573c prho------R-: Density at 1/2 layers
574c pri-------R-: Background Richardson Number, Wind shear measured along GW stress
575c pstab-----R-: Brunt-Vaisala freq. at 1/2 layers
576c pvph------R-: Wind in  plan of GW stress, Half levels.
577c ppsi------R-: Angle between low level wind and SS0 main axis.
578c pd1-------R-| Compared the ratio of the stress
579c pd2-------R-| that is along the wind to that Normal to it.
580c               pdi define the plane of low level stress
581c               compared to the low level wind.
582c see p. 108 Lott & Miller (1997).                     
583c pdmod-----R-: Norme of pdi
584
585c     === local arrays ===
586c
587c zvpf------R-: Wind projected in the plan of the low-level stress.
588
589c     ==== outputs ===
590c
591c        implicit arguments :   none
592c        --------------------
593c
594c     method.
595c     -------
596c
597c
598c     externals.
599c     ----------
600c
601c
602c     reference.
603c     ----------
604c
605c        see ecmwf research department documentation of the "i.f.s."
606c
607c     author.
608c     -------
609c
610c     modifications.
611c     --------------
612c     f.lott  for the new-gwdrag scheme november 1993
613c
614c-----------------------------------------------------------------------
615      USE dimphy
616      implicit none
617c
618
619cym#include "dimensions.h"
620cym#include "dimphy.h"
621#include "YOMCST.h"
622#include "YOEGWD.h"
623
624c-----------------------------------------------------------------------
625c
626c*       0.1   arguments
627c              ---------
628c
629      integer nlon,nlev
630      integer kkcrit(nlon),kkcrith(nlon),kcrit(nlon),ksect(nlon),
631     *        kkhlim(nlon),ktest(nlon),kkenvh(nlon)
632
633c
634      real paphm1(nlon,klev+1),papm1(nlon,klev),pum1(nlon,klev),
635     *     pvm1(nlon,klev),ptm1(nlon,klev),pgeom1(nlon,klev),
636     *     prho(nlon,klev+1),pri(nlon,klev+1),pstab(nlon,klev+1),
637     *     ptau(nlon,klev+1),pvph(nlon,klev+1),ppsi(nlon,klev+1),
638     *     pzdep(nlon,klev)
639       real pulow(nlon),pvlow(nlon),ptheta(nlon),pgam(nlon),pnu(nlon),
640     *     pd1(nlon),pd2(nlon),pdmod(nlon)
641      real pstd(nlon),pmea(nlon),ppic(nlon),pval(nlon)
642c
643c-----------------------------------------------------------------------
644c
645c*       0.2   local arrays
646c              ------------
647c
648c
649      integer ilevh ,jl,jk
650      real zcons1,zcons2,zhgeo,zu,zphi
651      real zvt1,zvt2,zdwind,zwind,zdelp
652      real zstabm,zstabp,zrhom,zrhop
653      logical lo
654      logical ll1(klon,klev+1)
655      integer kknu(klon),kknu2(klon),kknub(klon),kknul(klon),
656     *        kentp(klon),ncount(klon) 
657c
658      real zhcrit(klon,klev),zvpf(klon,klev),
659     *     zdp(klon,klev)
660      real znorm(klon),zb(klon),zc(klon),
661     *      zulow(klon),zvlow(klon),znup(klon),znum(klon)
662c       
663c     ------------------------------------------------------------------
664c
665c*         1.    initialization
666c                --------------
667c
668c       PRINT *,' in orosetup'
669 100  continue
670c
671c     ------------------------------------------------------------------
672c
673c*         1.1   computational constants
674c                -----------------------
675c
676 110  continue
677c
678      ilevh =klev/3
679c
680      zcons1=1./rd
681      zcons2=rg**2/rcpd
682c
683c
684c     ------------------------------------------------------------------
685c
686c*         2.
687c                --------------
688c
689 200  continue
690c
691c     ------------------------------------------------------------------
692c
693c*         2.1     define low level wind, project winds in plane of
694c*                 low level wind, determine sector in which to take
695c*                 the variance and set indicator for critical levels.
696c
697c
698c
699      do 2001 jl=kidia,kfdia
700      kknu(jl)    =klev
701      kknu2(jl)   =klev
702      kknub(jl)   =klev
703      kknul(jl)   =klev
704      pgam(jl) =max(pgam(jl),gtsec)
705      ll1(jl,klev+1)=.false.
706 2001 continue
707c
708c Ajouter une initialisation (L. Li, le 23fev99):
709c
710      do jk=klev,ilevh,-1
711      do jl=kidia,kfdia
712      ll1(jl,jk)= .false.
713      ENDDO
714      ENDDO
715c
716c*      define top of low level flow
717c       ----------------------------
718      do 2002 jk=klev,ilevh,-1
719      do 2003 jl=kidia,kfdia
720      if(ktest(jl).eq.1) then
721      lo=(paphm1(jl,jk)/paphm1(jl,klev+1)).ge.gsigcr
722      if(lo) then
723        kkcrit(jl)=jk
724      endif
725      zhcrit(jl,jk)=ppic(jl)-pval(jl)           
726      zhgeo=pgeom1(jl,jk)/rg
727      ll1(jl,jk)=(zhgeo.gt.zhcrit(jl,jk))
728      if(ll1(jl,jk).neqv.ll1(jl,jk+1)) then
729        kknu(jl)=jk
730      endif
731      if(.not.ll1(jl,ilevh))kknu(jl)=ilevh
732      endif
733 2003 continue
734 2002 continue
735      do 2004 jk=klev,ilevh,-1
736      do 2005 jl=kidia,kfdia
737      if(ktest(jl).eq.1) then
738      zhcrit(jl,jk)=ppic(jl)-pmea(jl)
739      zhgeo=pgeom1(jl,jk)/rg
740      ll1(jl,jk)=(zhgeo.gt.zhcrit(jl,jk))
741      if(ll1(jl,jk).neqv.ll1(jl,jk+1)) then
742        kknu2(jl)=jk
743      endif
744      if(.not.ll1(jl,ilevh))kknu2(jl)=ilevh
745      endif
746 2005 continue
747 2004 continue
748      do 2006 jk=klev,ilevh,-1
749      do 2007 jl=kidia,kfdia
750      if(ktest(jl).eq.1) then
751      zhcrit(jl,jk)=amin1(ppic(jl)-pmea(jl),pmea(jl)-pval(jl))
752      zhgeo=pgeom1(jl,jk)/rg
753      ll1(jl,jk)=(zhgeo.gt.zhcrit(jl,jk))
754      if(ll1(jl,jk).neqv.ll1(jl,jk+1)) then
755        kknub(jl)=jk
756      endif
757      if(.not.ll1(jl,ilevh))kknub(jl)=ilevh
758      endif
759 2007 continue
760 2006 continue
761c
762      do 2010 jl=kidia,kfdia 
763      if(ktest(jl).eq.1) then
764      kknu(jl)=min(kknu(jl),nktopg)
765      kknu2(jl)=min(kknu2(jl),nktopg)
766      kknub(jl)=min(kknub(jl),nktopg)
767      kknul(jl)=klev
768      endif
769 2010 continue     
770c
771 210  continue
772c
773cc*     initialize various arrays
774c
775      do 2107 jl=kidia,kfdia
776      prho(jl,klev+1)  =0.0
777cym correction en attendant mieux
778      prho(jl,1)  =0.0     
779      pstab(jl,klev+1) =0.0
780      pstab(jl,1)      =0.0
781      pri(jl,klev+1)   =9999.0
782      ppsi(jl,klev+1)  =0.0
783      pri(jl,1)        =0.0
784      pvph(jl,1)       =0.0
785      pvph(jl,klev+1)  =0.0
786cym correction en attendant mieux
787cym      pvph(jl,klev)    =0.0
788      pulow(jl)        =0.0
789      pvlow(jl)        =0.0
790      zulow(jl)        =0.0
791      zvlow(jl)        =0.0
792      kkcrith(jl)      =klev
793      kkenvh(jl)       =klev
794      kentp(jl)        =klev
795      kcrit(jl)        =1
796      ncount(jl)       =0
797      ll1(jl,klev+1)   =.false.
798 2107 continue
799c
800c*     define flow density and stratification (rho and N2)
801c      at semi layers.
802c      -------------------------------------------------------
803c
804      do 223 jk=klev,2,-1
805      do 222 jl=kidia,kfdia
806      if(ktest(jl).eq.1) then
807        zdp(jl,jk)=papm1(jl,jk)-papm1(jl,jk-1)
808        prho(jl,jk)=2.*paphm1(jl,jk)*zcons1/(ptm1(jl,jk)+ptm1(jl,jk-1))
809        pstab(jl,jk)=2.*zcons2/(ptm1(jl,jk)+ptm1(jl,jk-1))*
810     *  (1.-rcpd*prho(jl,jk)*(ptm1(jl,jk)-ptm1(jl,jk-1))/zdp(jl,jk))
811        pstab(jl,jk)=max(pstab(jl,jk),gssec)
812      endif
813  222 continue
814  223 continue
815c
816c********************************************************************
817c
818c*     define Low level flow (between ground and peacks-valleys)
819c      ---------------------------------------------------------
820      do 2115 jk=klev,ilevh,-1
821      do 2116 jl=kidia,kfdia
822      if(ktest(jl).eq.1)  then
823      if(jk.ge.kknu2(jl).and.jk.le.kknul(jl)) then
824        pulow(jl)=pulow(jl)+pum1(jl,jk)*(paphm1(jl,jk+1)-paphm1(jl,jk))
825        pvlow(jl)=pvlow(jl)+pvm1(jl,jk)*(paphm1(jl,jk+1)-paphm1(jl,jk))
826        pstab(jl,klev+1)=pstab(jl,klev+1)
827     c                   +pstab(jl,jk)*(paphm1(jl,jk+1)-paphm1(jl,jk))
828        prho(jl,klev+1)=prho(jl,klev+1)
829     c                   +prho(jl,jk)*(paphm1(jl,jk+1)-paphm1(jl,jk))
830      end if
831      endif
832 2116 continue
833 2115 continue
834      do 2110 jl=kidia,kfdia
835      if(ktest(jl).eq.1)  then
836      pulow(jl)=pulow(jl)/(paphm1(jl,kknul(jl)+1)-paphm1(jl,kknu2(jl)))
837      pvlow(jl)=pvlow(jl)/(paphm1(jl,kknul(jl)+1)-paphm1(jl,kknu2(jl)))
838      znorm(jl)=max(sqrt(pulow(jl)**2+pvlow(jl)**2),gvsec)
839      pvph(jl,klev+1)=znorm(jl)
840      pstab(jl,klev+1)=pstab(jl,klev+1)
841     c                /(paphm1(jl,kknul(jl)+1)-paphm1(jl,kknu2(jl)))
842      prho(jl,klev+1)=prho(jl,klev+1)
843     c                /(paphm1(jl,kknul(jl)+1)-paphm1(jl,kknu2(jl)))
844      endif
845 2110 continue
846
847c
848c*******  setup orography orientation relative to the low level
849C       wind and define parameters of the Anisotropic wave stress.
850c
851      do 2112 jl=kidia,kfdia
852      if(ktest(jl).eq.1)  then
853        lo=(pulow(jl).lt.gvsec).and.(pulow(jl).ge.-gvsec)
854        if(lo) then
855          zu=pulow(jl)+2.*gvsec
856        else
857          zu=pulow(jl)
858        endif
859        zphi=atan(pvlow(jl)/zu)
860        ppsi(jl,klev+1)=ptheta(jl)*rpi/180.-zphi
861        zb(jl)=1.-0.18*pgam(jl)-0.04*pgam(jl)**2
862        zc(jl)=0.48*pgam(jl)+0.3*pgam(jl)**2
863        pd1(jl)=zb(jl)-(zb(jl)-zc(jl))*(sin(ppsi(jl,klev+1))**2)
864        pd2(jl)=(zb(jl)-zc(jl))*sin(ppsi(jl,klev+1))
865     *                         *cos(ppsi(jl,klev+1))
866        pdmod(jl)=sqrt(pd1(jl)**2+pd2(jl)**2)
867      endif
868 2112 continue
869c
870c  ************ projet flow in plane of lowlevel stress *************
871C  ************ Find critical levels...                 *************
872c
873      do 213 jk=1,klev
874      do 212 jl=kidia,kfdia
875      if(ktest(jl).eq.1)  then
876        zvt1       =pulow(jl)*pum1(jl,jk)+pvlow(jl)*pvm1(jl,jk)
877        zvt2       =-pvlow(jl)*pum1(jl,jk)+pulow(jl)*pvm1(jl,jk)
878        zvpf(jl,jk)=(zvt1*pd1(jl)+zvt2*pd2(jl))/(znorm(jl)*pdmod(jl))
879      endif
880      ptau(jl,jk)  =0.0
881      pzdep(jl,jk) =0.0
882      ppsi(jl,jk)  =0.0
883      ll1(jl,jk)   =.false.
884  212 continue
885  213 continue
886      do 215 jk=2,klev
887      do 214 jl=kidia,kfdia
888      if(ktest(jl).eq.1) then
889        zdp(jl,jk)=papm1(jl,jk)-papm1(jl,jk-1)
890        pvph(jl,jk)=((paphm1(jl,jk)-papm1(jl,jk-1))*zvpf(jl,jk)+
891     *            (papm1(jl,jk)-paphm1(jl,jk))*zvpf(jl,jk-1))
892     *            /zdp(jl,jk)
893        if(pvph(jl,jk).lt.gvsec) then
894          pvph(jl,jk)=gvsec
895          kcrit(jl)=jk
896        endif
897      endif
898  214 continue
899  215 continue
900c
901c*         2.3     mean flow richardson number.
902c
903  230 continue
904c
905      do 232 jk=2,klev
906      do 231 jl=kidia,kfdia
907      if(ktest(jl).eq.1) then
908        zdwind=max(abs(zvpf(jl,jk)-zvpf(jl,jk-1)),gvsec)
909        pri(jl,jk)=pstab(jl,jk)*(zdp(jl,jk)
910     *          /(rg*prho(jl,jk)*zdwind))**2
911        pri(jl,jk)=max(pri(jl,jk),grcrit)
912      endif
913  231 continue
914  232 continue
915 
916c
917c
918c*      define top of 'envelope' layer
919c       ----------------------------
920
921      do 233 jl=kidia,kfdia
922      pnu (jl)=0.0
923      znum(jl)=0.0
924 233  continue
925     
926      do 234 jk=2,klev-1
927      do 234 jl=kidia,kfdia
928     
929      if(ktest(jl).eq.1) then
930       
931      if (jk.ge.kknu2(jl)) then
932         
933            znum(jl)=pnu(jl)
934            zwind=(pulow(jl)*pum1(jl,jk)+pvlow(jl)*pvm1(jl,jk))/
935     *            max(sqrt(pulow(jl)**2+pvlow(jl)**2),gvsec)
936            zwind=max(sqrt(zwind**2),gvsec)
937            zdelp=paphm1(jl,jk+1)-paphm1(jl,jk)
938            zstabm=sqrt(max(pstab(jl,jk  ),gssec))
939            zstabp=sqrt(max(pstab(jl,jk+1),gssec))
940            zrhom=prho(jl,jk  )
941            zrhop=prho(jl,jk+1)
942            pnu(jl) = pnu(jl) + (zdelp/rg)*
943     *            ((zstabp/zrhop+zstabm/zrhom)/2.)/zwind     
944            if((znum(jl).le.gfrcrit).and.(pnu(jl).gt.gfrcrit)
945     *                          .and.(kkenvh(jl).eq.klev))
946     *      kkenvh(jl)=jk
947     
948      endif   
949
950      endif
951     
952 234  continue
953     
954c  calculation of a dynamical mixing height for when the waves
955C  BREAK AT LOW LEVEL: The drag will be repartited over
956C  a depths that depends on waves vertical wavelength,
957C  not just between two adjacent model layers.
958c  of gravity waves:
959
960      do 235 jl=kidia,kfdia
961      znup(jl)=0.0
962      znum(jl)=0.0
963 235  continue
964
965      do 236 jk=klev-1,2,-1
966      do 236 jl=kidia,kfdia
967     
968      if(ktest(jl).eq.1) then
969
970            znum(jl)=znup(jl)
971            zwind=(pulow(jl)*pum1(jl,jk)+pvlow(jl)*pvm1(jl,jk))/
972     *            max(sqrt(pulow(jl)**2+pvlow(jl)**2),gvsec)
973            zwind=max(sqrt(zwind**2),gvsec)
974            zdelp=paphm1(jl,jk+1)-paphm1(jl,jk)
975            zstabm=sqrt(max(pstab(jl,jk  ),gssec))
976            zstabp=sqrt(max(pstab(jl,jk+1),gssec))
977            zrhom=prho(jl,jk  )
978            zrhop=prho(jl,jk+1)
979            znup(jl) = znup(jl) + (zdelp/rg)*
980     *            ((zstabp/zrhop+zstabm/zrhom)/2.)/zwind     
981            if((znum(jl).le.rpi/4.).and.(znup(jl).gt.rpi/4.)
982     *                          .and.(kkcrith(jl).eq.klev))
983     *      kkcrith(jl)=jk
984     
985      endif
986     
987 236  continue
988 
989      do 237 jl=kidia,kfdia
990      if(ktest(jl).eq.1) then
991      kkcrith(jl)=max0(kkcrith(jl),ilevh*2)
992      kkcrith(jl)=max0(kkcrith(jl),kknu(jl))
993      if(kcrit(jl).ge.kkcrith(jl))kcrit(jl)=1
994      endif
995 237  continue         
996c
997c     directional info for flow blocking *************************
998c
999      do 251 jk=1,klev   
1000      do 252 jl=kidia,kfdia
1001      if(ktest(jl).eq.1) then
1002      lo=(pum1(jl,jk).lt.gvsec).and.(pum1(jl,jk).ge.-gvsec)
1003      if(lo) then
1004        zu=pum1(jl,jk)+2.*gvsec
1005      else
1006        zu=pum1(jl,jk)
1007      endif
1008       zphi=atan(pvm1(jl,jk)/zu)
1009       ppsi(jl,jk)=ptheta(jl)*rpi/180.-zphi
1010      endif
1011 252  continue
1012 251  continue
1013
1014c      forms the vertical 'leakiness' **************************
1015
1016      do 254  jk=ilevh,klev
1017      do 253  jl=kidia,kfdia
1018      if(ktest(jl).eq.1) then
1019      pzdep(jl,jk)=0
1020      if(jk.ge.kkenvh(jl).and.kkenvh(jl).ne.klev) then
1021        pzdep(jl,jk)=(pgeom1(jl,kkenvh(jl)  )-pgeom1(jl,  jk))/
1022     *               (pgeom1(jl,kkenvh(jl)  )-pgeom1(jl,klev))
1023      end if
1024      endif
1025 253  continue
1026 254  continue
1027
1028      return
1029      end
1030      SUBROUTINE gwstress_strato
1031     *         (  nlon  , nlev
1032     *         , kkcrit, ksect, kkhlim, ktest, kkcrith, kcrit, kkenvh
1033     *         , kknu
1034     *         , prho  , pstab , pvph  , pstd, psig
1035     *         , pmea , ppic , pval  , ptfr  , ptau 
1036     *         , pgeom1 , pgamma , pd1  , pd2   , pdmod , pnu )
1037c
1038c**** *gwstress*
1039c
1040c     purpose.
1041c     --------
1042c  Compute the surface stress due to Gravity Waves, according
1043c  to the Phillips (1979) theory of 3-D flow above
1044c  anisotropic elliptic ridges.
1045
1046C  The stress is reduced two account for cut-off flow over
1047C  hill.  The flow only see that part of the ridge located
1048c  above the blocked layer (see zeff).
1049c
1050c**   interface.
1051c     ----------
1052c     call *gwstress*  from *gwdrag*
1053c
1054c        explicit arguments :
1055c        --------------------
1056c     ==== inputs ===
1057c     ==== outputs ===
1058c
1059c        implicit arguments :   none
1060c        --------------------
1061c
1062c     method.
1063c     -------
1064c
1065c
1066c     externals.
1067c     ----------
1068c
1069c
1070c     reference.
1071c     ----------
1072c
1073c   LOTT and MILLER (1997)  &  LOTT (1999)
1074c
1075c     author.
1076c     -------
1077c
1078c     modifications.
1079c     --------------
1080c     f. lott put the new gwd on ifs      22/11/93
1081c
1082c-----------------------------------------------------------------------
1083      USE dimphy
1084      implicit none
1085
1086cym#include "dimensions.h"
1087cym#include "dimphy.h"
1088#include "YOMCST.h"
1089#include "YOEGWD.h"
1090
1091c-----------------------------------------------------------------------
1092c
1093c*       0.1   arguments
1094c              ---------
1095c
1096      integer nlon,nlev
1097      integer kkcrit(nlon),kkcrith(nlon),kcrit(nlon),ksect(nlon),
1098     *        kkhlim(nlon),ktest(nlon),kkenvh(nlon),kknu(nlon)
1099c
1100      real prho(nlon,nlev+1),pstab(nlon,nlev+1),ptau(nlon,nlev+1),
1101     *     pvph(nlon,nlev+1),ptfr(nlon),
1102     *     pgeom1(nlon,nlev),pstd(nlon)
1103c
1104      real pd1(nlon),pd2(nlon),pnu(nlon),psig(nlon),pgamma(nlon)
1105      real pmea(nlon),ppic(nlon),pval(nlon)
1106      real pdmod(nlon)
1107c
1108c-----------------------------------------------------------------------
1109c
1110c*       0.2   local arrays
1111c              ------------
1112c  zeff--real: effective height seen by the flow when there is blocking
1113
1114      integer jl
1115      real zeff 
1116c
1117c-----------------------------------------------------------------------
1118c
1119c*       0.3   functions
1120c              ---------
1121c     ------------------------------------------------------------------
1122c
1123c*         1.    initialization
1124c                --------------
1125c
1126c      PRINT *,' in gwstress'
1127 100  continue
1128c
1129c*         3.1     gravity wave stress.
1130c
1131  300 continue
1132c
1133c
1134      do 301 jl=kidia,kfdia
1135      if(ktest(jl).eq.1) then
1136     
1137c  effective mountain height above the blocked flow
1138 
1139         zeff=ppic(jl)-pval(jl)
1140         if(kkenvh(jl).lt.klev)then
1141         zeff=amin1(GFRCRIT*pvph(jl,klev+1)/sqrt(pstab(jl,klev+1))
1142     c              ,zeff)
1143         endif
1144
1145     
1146        ptau(jl,klev+1)=gkdrag*prho(jl,klev+1)
1147     *     *psig(jl)*pdmod(jl)/4./pstd(jl)
1148     *     *pvph(jl,klev+1)*sqrt(pstab(jl,klev+1))
1149     *     *zeff**2
1150
1151
1152c  too small value of stress or  low level flow include critical level
1153c  or low level flow:  gravity wave stress nul.
1154               
1155c       lo=(ptau(jl,klev+1).lt.gtsec).or.(kcrit(jl).ge.kknu(jl))
1156c    *      .or.(pvph(jl,klev+1).lt.gvcrit)
1157c       if(lo) ptau(jl,klev+1)=0.0
1158     
1159c      print *,jl,ptau(jl,klev+1)
1160
1161      else
1162     
1163          ptau(jl,klev+1)=0.0
1164         
1165      endif
1166
1167  301 continue
1168
1169c      write(21)(ptau(jl,klev+1),jl=kidia,kfdia)
1170 
1171      return
1172      end
1173
1174      subroutine gwprofil_strato
1175     *         ( nlon, nlev
1176     *         , kgwd ,kdx  , ktest
1177     *         , kkcrit, kkcrith, kcrit ,  kkenvh, kknu,kknu2
1178     *         , paphm1, prho   , pstab , ptfr , pvph , pri , ptau
1179     *         , pdmod   , pnu   , psig ,pgamma, pstd, ppic,pval)
1180
1181C**** *gwprofil*
1182C
1183C     purpose.
1184C     --------
1185C
1186C**   interface.
1187C     ----------
1188C          from *gwdrag*
1189C
1190C        explicit arguments :
1191C        --------------------
1192C     ==== inputs ===
1193C
1194C     ==== outputs ===
1195C
1196C        implicit arguments :   none
1197C        --------------------
1198C
1199C     method:
1200C     -------
1201C     the stress profile for gravity waves is computed as follows:
1202C     it decreases linearly with heights from the ground
1203C     to the low-level indicated by kkcrith,
1204C     to simulates lee waves or
1205C     low-level gravity wave breaking.
1206C     above it is constant, except when the waves encounter a critical
1207C     level (kcrit) or when they break.
1208C     The stress is also uniformly distributed above the level
1209C     nstra.                                         
1210C
1211      USE dimphy
1212      IMPLICIT NONE
1213
1214cym#include "dimensions.h"
1215cym#include "dimphy.h"
1216#include "YOMCST.h"
1217#include "YOEGWD.h"
1218
1219C-----------------------------------------------------------------------
1220C
1221C*       0.1   ARGUMENTS
1222C              ---------
1223C
1224      integer nlon,nlev,kgwd
1225      integer kkcrit(nlon),kkcrith(nlon),kcrit(nlon)
1226     *       ,kdx(nlon),ktest(nlon)
1227     *       ,kkenvh(nlon),kknu(nlon),kknu2(nlon)
1228C
1229      real paphm1(nlon,nlev+1), pstab(nlon,nlev+1),
1230     *     prho  (nlon,nlev+1), pvph (nlon,nlev+1),
1231     *     pri   (nlon,nlev+1), ptfr (nlon), ptau(nlon,nlev+1)
1232     
1233      real pdmod (nlon) , pnu (nlon) , psig(nlon),
1234     *     pgamma(nlon) , pstd(nlon) , ppic(nlon), pval(nlon)
1235     
1236C-----------------------------------------------------------------------
1237C
1238C*       0.2   local arrays
1239C              ------------
1240C
1241      integer jl,jk
1242      real zsqr,zalfa,zriw,zdel,zb,zalpha,zdz2n,zdelp,zdelpt
1243
1244      real zdz2 (klon,klev) , znorm(klon) , zoro(klon)
1245      real ztau (klon,klev+1)
1246C
1247C-----------------------------------------------------------------------
1248C
1249C*         1.    INITIALIZATION
1250C                --------------
1251C
1252C      print *,' entree gwprofil'
1253 100  CONTINUE
1254C
1255C
1256C*    COMPUTATIONAL CONSTANTS.
1257C     ------------- ----------
1258C
1259      do 400 jl=kidia,kfdia
1260      if(ktest(jl).eq.1)then
1261      zoro(jl)=psig(jl)*pdmod(jl)/4./pstd(jl)
1262      ztau(jl,klev+1)=ptau(jl,klev+1)
1263c     print *,jl,ptau(jl,klev+1)
1264      ztau(jl,kkcrith(jl))=grahilo*ptau(jl,klev+1)
1265      endif
1266  400 continue
1267 
1268C
1269      do 430 jk=klev+1,1,-1
1270C
1271C
1272C*         4.1    constant shear stress until top of the
1273C                 low-level breaking/trapped layer
1274  410 CONTINUE
1275C
1276      do 411 jl=kidia,kfdia
1277      if(ktest(jl).eq.1)then
1278           if(jk.gt.kkcrith(jl)) then
1279           zdelp=paphm1(jl,jk)-paphm1(jl,klev+1)
1280           zdelpt=paphm1(jl,kkcrith(jl))-paphm1(jl,klev+1)
1281           ptau(jl,jk)=ztau(jl,klev+1)+zdelp/zdelpt*
1282     c                 (ztau(jl,kkcrith(jl))-ztau(jl,klev+1))
1283           else                   
1284           ptau(jl,jk)=ztau(jl,kkcrith(jl))
1285           endif
1286       endif
1287 411  continue             
1288C
1289C*         4.15   constant shear stress until the top of the
1290C                 low level flow layer.
1291 415  continue
1292C       
1293C
1294C*         4.2    wave displacement at next level.
1295C
1296  420 continue
1297C
1298  430 continue
1299
1300C
1301C*         4.4    wave richardson number, new wave displacement
1302C*                and stress:  breaking evaluation and critical
1303C                 level
1304C
1305                         
1306      do 440 jk=klev,1,-1
1307
1308      do 441 jl=kidia,kfdia
1309      if(ktest(jl).eq.1)then
1310      znorm(jl)=prho(jl,jk)*sqrt(pstab(jl,jk))*pvph(jl,jk)
1311      zdz2(jl,jk)=ptau(jl,jk)/amax1(znorm(jl),gssec)/zoro(jl)
1312      endif
1313  441 continue
1314
1315      do 442 jl=kidia,kfdia
1316      if(ktest(jl).eq.1)then
1317          if(jk.lt.kkcrith(jl)) then
1318          if((ptau(jl,jk+1).lt.gtsec).or.(jk.le.kcrit(jl))) then
1319             ptau(jl,jk)=0.0
1320          else
1321               zsqr=sqrt(pri(jl,jk))
1322               zalfa=sqrt(pstab(jl,jk)*zdz2(jl,jk))/pvph(jl,jk)
1323               zriw=pri(jl,jk)*(1.-zalfa)/(1+zalfa*zsqr)**2
1324               if(zriw.lt.grcrit) then
1325C                 print *,' breaking!!!',ptau(jl,jk)
1326                  zdel=4./zsqr/grcrit+1./grcrit**2+4./grcrit
1327                  zb=1./grcrit+2./zsqr
1328                  zalpha=0.5*(-zb+sqrt(zdel))
1329                  zdz2n=(pvph(jl,jk)*zalpha)**2/pstab(jl,jk)
1330                  ptau(jl,jk)=znorm(jl)*zdz2n*zoro(jl)
1331               endif
1332               
1333               ptau(jl,jk)=amin1(ptau(jl,jk),ptau(jl,jk+1))
1334                 
1335          endif
1336          endif
1337      endif
1338  442 continue
1339  440 continue
1340
1341C  REORGANISATION OF THE STRESS PROFILE AT LOW LEVEL
1342
1343      do 530 jl=kidia,kfdia
1344      if(ktest(jl).eq.1)then
1345         ztau(jl,kkcrith(jl)-1)=ptau(jl,kkcrith(jl)-1)
1346         ztau(jl,nstra)=ptau(jl,nstra)
1347      endif
1348 530  continue     
1349
1350      do 531 jk=1,klev
1351     
1352      do 532 jl=kidia,kfdia
1353      if(ktest(jl).eq.1)then
1354               
1355         if(jk.gt.kkcrith(jl)-1)then
1356
1357          zdelp=paphm1(jl,jk)-paphm1(jl,klev+1    )
1358          zdelpt=paphm1(jl,kkcrith(jl)-1)-paphm1(jl,klev+1    )
1359          ptau(jl,jk)=ztau(jl,klev+1    ) +
1360     .                (ztau(jl,kkcrith(jl)-1)-ztau(jl,klev+1    ) )*
1361     .                zdelp/zdelpt
1362     
1363        endif
1364      endif
1365           
1366 532  continue   
1367 
1368C  REORGANISATION AT THE MODEL TOP....
1369
1370      do 533 jl=kidia,kfdia
1371      if(ktest(jl).eq.1)then
1372
1373         if(jk.lt.nstra)then
1374
1375          zdelp =paphm1(jl,nstra)
1376          zdelpt=paphm1(jl,jk)
1377          ptau(jl,jk)=ztau(jl,nstra)*zdelpt/zdelp
1378c         ptau(jl,jk)=ztau(jl,nstra)               
1379
1380        endif
1381
1382      endif
1383
1384 533  continue
1385
1386 
1387 531  continue       
1388
1389
1390 123   format(i4,1x,20(f6.3,1x))
1391
1392
1393      return
1394      end
1395      subroutine lift_noro_strato (nlon,nlev,dtime,paprs,pplay,     
1396     i                   plat,pmea,pstd, psig, pgam, pthe, ppic,pval,
1397     i                   kgwd,kdx,ktest,
1398     i                   t, u, v,
1399     o                   pulow, pvlow, pustr, pvstr,
1400     o                   d_t, d_u, d_v)
1401c
1402      USE dimphy
1403      implicit none
1404c======================================================================
1405c Auteur(s): F.Lott (LMD/CNRS) date: 19950201
1406c Object: Mountain lift interface (enhanced vortex stretching).
1407c         Made necessary because:
1408C 1. in the LMD-GCM Layers are from bottom to top,
1409C    contrary to most European GCM.
1410c 2. the altitude above ground of each model layers
1411c    needs to be known (variable zgeom)
1412c======================================================================
1413c Explicit Arguments:
1414c ==================
1415c nlon----input-I-Total number of horizontal points that get into physics
1416c nlev----input-I-Number of vertical levels
1417c dtime---input-R-Time-step (s)
1418c paprs---input-R-Pressure in semi layers    (Pa)
1419c pplay---input-R-Pressure model-layers      (Pa)
1420c t-------input-R-temperature (K)
1421c u-------input-R-Horizontal wind (m/s)
1422c v-------input-R-Meridional wind (m/s)
1423c pmea----input-R-Mean Orography (m)
1424C pstd----input-R-SSO standard deviation (m)
1425c psig----input-R-SSO slope
1426c pgam----input-R-SSO Anisotropy
1427c pthe----input-R-SSO Angle
1428c ppic----input-R-SSO Peacks elevation (m)
1429c pval----input-R-SSO Valleys elevation (m)
1430c
1431c kgwd- -input-I: Total nb of points where the orography schemes are active
1432c ktest--input-I: Flags to indicate active points
1433c kdx----input-I: Locate the physical location of an active point.
1434
1435c pulow, pvlow -output-R: Low-level wind
1436c pustr, pvstr -output-R: Surface stress due to SSO drag      (Pa)
1437c
1438c d_t-----output-R: T increment
1439c d_u-----output-R: U increment
1440c d_v-----output-R: V increment
1441c
1442c Implicit Arguments:
1443c ===================
1444c
1445c iim--common-I: Number of longitude intervals
1446c jjm--common-I: Number of latitude intervals
1447c klon-common-I: Number of points seen by the physics
1448c                (iim+1)*(jjm+1) for instance
1449c klev-common-I: Number of vertical layers
1450c======================================================================
1451c Local Variables:
1452c ================
1453c
1454c zgeom-----R: Altitude of layer above ground
1455c pt, pu, pv --R: t u v from top to bottom
1456c pdtdt, pdudt, pdvdt --R: t u v tendencies (from top to bottom)
1457c papmf: pressure at model layer (from top to bottom)
1458c papmh: pressure at model 1/2 layer (from top to bottom)
1459c
1460c======================================================================
1461
1462cym#include "dimensions.h"
1463cym#include "dimphy.h"
1464#include "YOMCST.h"
1465#include "YOEGWD.h"
1466c
1467c ARGUMENTS
1468c
1469      INTEGER nlon,nlev
1470      REAL dtime
1471      REAL paprs(klon,klev+1)
1472      REAL pplay(klon,klev)
1473      REAL plat(nlon),pmea(nlon)
1474      REAL pstd(nlon),psig(nlon),pgam(nlon),pthe(nlon)
1475      REAL ppic(nlon),pval(nlon)
1476      REAL pulow(nlon),pvlow(nlon),pustr(nlon),pvstr(nlon)
1477      REAL t(nlon,nlev), u(nlon,nlev), v(nlon,nlev)
1478      REAL d_t(nlon,nlev), d_u(nlon,nlev), d_v(nlon,nlev)
1479c
1480      INTEGER i, k, kgwd,  kdx(nlon), ktest(nlon)
1481c
1482c Variables locales:
1483c
1484      REAL zgeom(klon,klev)
1485      REAL pdtdt(klon,klev), pdudt(klon,klev), pdvdt(klon,klev)
1486      REAL pt(klon,klev), pu(klon,klev), pv(klon,klev)
1487      REAL papmf(klon,klev),papmh(klon,klev+1)
1488c
1489c initialiser les variables de sortie (pour securite)
1490c
1491
1492c     print *,'in lift_noro'
1493      DO i = 1,klon
1494         pulow(i) = 0.0
1495         pvlow(i) = 0.0
1496         pustr(i) = 0.0
1497         pvstr(i) = 0.0
1498      ENDDO
1499      DO k = 1, klev
1500      DO i = 1, klon
1501         d_t(i,k) = 0.0
1502         d_u(i,k) = 0.0
1503         d_v(i,k) = 0.0
1504         pdudt(i,k)=0.0
1505         pdvdt(i,k)=0.0
1506         pdtdt(i,k)=0.0
1507      ENDDO
1508      ENDDO
1509c
1510c preparer les variables d'entree (attention: l'ordre des niveaux
1511c verticaux augmente du haut vers le bas)
1512c
1513      DO k = 1, klev
1514      DO i = 1, klon
1515         pt(i,k) = t(i,klev-k+1)
1516         pu(i,k) = u(i,klev-k+1)
1517         pv(i,k) = v(i,klev-k+1)
1518         papmf(i,k) = pplay(i,klev-k+1)
1519      ENDDO
1520      ENDDO
1521      DO k = 1, klev+1
1522      DO i = 1, klon
1523         papmh(i,k) = paprs(i,klev-k+2)
1524      ENDDO
1525      ENDDO
1526      DO i = 1, klon
1527         zgeom(i,klev) = RD * pt(i,klev)
1528     .                  * LOG(papmh(i,klev+1)/papmf(i,klev))
1529      ENDDO
1530      DO k = klev-1, 1, -1
1531      DO i = 1, klon
1532         zgeom(i,k) = zgeom(i,k+1) + RD * (pt(i,k)+pt(i,k+1))/2.0
1533     .               * LOG(papmf(i,k+1)/papmf(i,k))
1534      ENDDO
1535      ENDDO
1536c
1537c appeler la routine principale
1538c
1539
1540      CALL OROLIFT_strato(klon,klev,kgwd,kdx,ktest,
1541     .            dtime,
1542     .            papmh, papmf, zgeom,
1543     .            pt, pu, pv,
1544     .            plat,pmea, pstd, psig, pgam, pthe, ppic,pval,
1545     .            pulow,pvlow,
1546     .            pdudt,pdvdt,pdtdt)
1547C
1548      DO k = 1, klev
1549      DO i = 1, klon
1550         d_u(i,klev+1-k) = dtime*pdudt(i,k)
1551         d_v(i,klev+1-k) = dtime*pdvdt(i,k)
1552         d_t(i,klev+1-k) = dtime*pdtdt(i,k)
1553         pustr(i)        = pustr(i)
1554     .                    +pdudt(i,k)*(papmh(i,k+1)-papmh(i,k))/rg
1555         pvstr(i)        = pvstr(i)
1556     .                    +pdvdt(i,k)*(papmh(i,k+1)-papmh(i,k))/rg
1557      ENDDO
1558      ENDDO
1559
1560c     print *,' out lift_noro'
1561c
1562      RETURN
1563      END
1564      subroutine orolift_strato( nlon,nlev
1565     I                 , kgwd, kdx, ktest
1566     R                 , ptsphy
1567     R                 , paphm1,papm1,pgeom1,ptm1,pum1,pvm1
1568     R                 , plat
1569     R                 , pmea, pstd, psig, pgam, pthe,ppic,pval
1570C OUTPUTS
1571     R                 , pulow,pvlow
1572     R                 , pvom,pvol,pte )
1573
1574C
1575C**** *OROLIFT: SIMULATE THE GEOSTROPHIC LIFT.
1576C
1577C     PURPOSE.
1578C     --------
1579C this routine computes the physical tendencies of the
1580C prognostic variables u,v  when enhanced vortex stretching
1581C is needed.
1582C
1583C**   INTERFACE.
1584C     ----------
1585C          CALLED FROM *lift_noro
1586c        explicit arguments :
1587c        --------------------
1588c     ==== inputs ===
1589c nlon----input-I-Total number of horizontal points that get into physics
1590c nlev----input-I-Number of vertical levels
1591c
1592c kgwd- -input-I: Total nb of points where the orography schemes are active
1593c ktest--input-I: Flags to indicate active points
1594c kdx----input-I: Locate the physical location of an active point.
1595c ptsphy--input-R-Time-step (s)
1596c paphm1--input-R: pressure at model 1/2 layer
1597c papm1---input-R: pressure at model layer
1598c pgeom1--input-R: Altitude of layer above ground
1599c ptm1, pum1, pvm1--R-: t, u and v
1600c pmea----input-R-Mean Orography (m)
1601C pstd----input-R-SSO standard deviation (m)
1602c psig----input-R-SSO slope
1603c pgam----input-R-SSO Anisotropy
1604c pthe----input-R-SSO Angle
1605c ppic----input-R-SSO Peacks elevation (m)
1606c pval----input-R-SSO Valleys elevation (m)
1607c plat----input-R-Latitude (degree)
1608c
1609c     ==== outputs ===
1610c pulow, pvlow -output-R: Low-level wind
1611c
1612c pte -----output-R: T tendency
1613c pvom-----output-R: U tendency
1614c pvol-----output-R: V tendency
1615c
1616c
1617c Implicit Arguments:
1618c ===================
1619c
1620c klon-common-I: Number of points seen by the physics
1621c klev-common-I: Number of vertical layers
1622c
1623
1624C     ----------
1625C
1626C     AUTHOR.
1627C     -------
1628C     F.LOTT  LMD 22/11/95
1629C
1630       USE dimphy
1631       implicit none
1632C
1633C
1634cym#include "dimensions.h"
1635cym#include "dimphy.h"
1636#include "YOMCST.h"
1637#include "YOEGWD.h"
1638C-----------------------------------------------------------------------
1639C
1640C*       0.1   ARGUMENTS
1641C              ---------
1642C
1643C
1644      integer nlon,nlev,kgwd
1645      real ptsphy
1646      real  pte(nlon,nlev),
1647     *      pvol(nlon,nlev),
1648     *      pvom(nlon,nlev),
1649     *      pulow(nlon),
1650     *      pvlow(nlon)
1651      real  pum1(nlon,nlev),
1652     *      pvm1(nlon,nlev),
1653     *      ptm1(nlon,nlev),
1654     *      plat(nlon),pmea(nlon),
1655     *      pstd(nlon),psig(nlon),pgam(nlon),
1656     *      pthe(nlon),ppic(nlon),pval(nlon),
1657     *      pgeom1(nlon,nlev),
1658     *      papm1(nlon,nlev),
1659     *      paphm1(nlon,nlev+1)
1660C
1661      INTEGER  KDX(NLON),KTEST(NLON)
1662C-----------------------------------------------------------------------
1663C
1664C*       0.2   local arrays
1665
1666      integer jl,ilevh,jk
1667      real zhgeo,zdelp,zslow,zsqua,zscav,zbet
1668C              ------------
1669      integer  iknub(klon),
1670     *         iknul(klon)
1671      logical ll1(klon,klev+1)
1672C
1673      real   ztau(klon,klev+1),
1674     *       ztav(klon,klev+1),
1675     *       zrho(klon,klev+1)
1676      real   zdudt(klon),
1677     *       zdvdt(klon)
1678      real zhcrit(klon,klev)
1679
1680      logical lifthigh
1681      real zcons1,ztmst
1682
1683C-----------------------------------------------------------------------
1684C
1685C*         1.1  initialisations
1686C               ---------------
1687
1688      lifthigh=.false.
1689
1690      if(nlon.ne.klon.or.nlev.ne.klev)stop
1691      zcons1=1./rd
1692      ztmst=ptsphy
1693C
1694      do 1001 jl=kidia,kfdia
1695      zrho(jl,klev+1)  =0.0
1696      pulow(jl)        =0.0
1697      pvlow(jl)        =0.0
1698      iknub(JL)   =klev
1699      iknul(JL)   =klev
1700      ilevh=klev/3
1701      ll1(jl,klev+1)=.false.
1702      do 1000 jk=1,klev
1703      pvom(jl,jk)=0.0
1704      pvol(jl,jk)=0.0
1705      pte (jl,jk)=0.0
1706 1000 continue
1707 1001 continue
1708
1709C
1710C*         2.1     DEFINE LOW LEVEL WIND, PROJECT WINDS IN PLANE OF
1711C*                 LOW LEVEL WIND, DETERMINE SECTOR IN WHICH TO TAKE
1712C*                 THE VARIANCE AND SET INDICATOR FOR CRITICAL LEVELS.
1713C
1714C
1715C
1716      do 2006 jk=klev,1,-1
1717      do 2007 jl=kidia,kfdia
1718      if(ktest(jl).eq.1) then
1719      zhcrit(jl,jk)=amax1(ppic(jl)-pval(jl),100.)
1720      zhgeo=pgeom1(jl,jk)/rg
1721      ll1(jl,jk)=(zhgeo.gt.zhcrit(jl,jk))
1722      if(ll1(jl,jk).neqv.ll1(jl,jk+1)) then
1723        iknub(jl)=jk
1724      endif
1725      endif
1726 2007 continue
1727 2006 continue
1728C
1729
1730      do 2010 jl=kidia,kfdia
1731      if(ktest(jl).eq.1) then
1732      iknub(jl)=max(iknub(jl),klev/2)
1733      iknul(jl)=max(iknul(jl),2*klev/3)
1734      if(iknub(jl).gt.nktopg) iknub(jl)=nktopg
1735      if(iknub(jl).eq.nktopg) iknul(jl)=klev
1736      if(iknub(jl).eq.iknul(jl)) iknub(jl)=iknul(jl)-1
1737      endif
1738 2010 continue
1739
1740      do 223 jk=klev,2,-1
1741      do 222 jl=kidia,kfdia
1742        zrho(jl,jk)=2.*paphm1(jl,jk)*zcons1/(ptm1(jl,jk)+ptm1(jl,jk-1))
1743  222 continue
1744  223 continue
1745c     print *,'  dans orolift: 223'
1746
1747C********************************************************************
1748C
1749c*     define low level flow
1750C      -------------------
1751      do 2115 jk=klev,1,-1
1752      do 2116 jl=kidia,kfdia
1753      if(ktest(jl).eq.1) THEN
1754      if(jk.ge.iknub(jl).and.jk.le.iknul(jl)) then
1755        pulow(JL)=pulow(JL)+pum1(jl,jk)*(paphm1(jl,jk+1)-paphm1(jl,jk))
1756        pvlow(JL)=pvlow(JL)+pvm1(jl,jk)*(paphm1(jl,jk+1)-paphm1(jl,jk))
1757        zrho(JL,klev+1)=zrho(JL,klev+1)
1758     *                 +zrho(JL,JK)*(paphm1(jl,jk+1)-paphm1(jl,jk))
1759      endif
1760      endif
1761 2116 continue
1762 2115 continue
1763      do 2110 jl=kidia,kfdia
1764      if(ktest(jl).eq.1) then
1765      pulow(JL)=pulow(JL)/(paphm1(jl,iknul(jl)+1)-paphm1(jl,iknub(jl)))
1766      pvlow(JL)=pvlow(JL)/(paphm1(jl,iknul(jl)+1)-paphm1(jl,iknub(jl)))
1767      zrho(JL,klev+1)=zrho(Jl,klev+1)
1768     *               /(paphm1(jl,iknul(jl)+1)-paphm1(jl,iknub(jl)))
1769      endif
1770 2110 continue
1771
1772
1773200   continue
1774
1775C***********************************************************
1776C
1777C*         3.      COMPUTE MOUNTAIN LIFT
1778C
1779  300 continue
1780C
1781      do 301 jl=kidia,kfdia
1782      if(ktest(jl).eq.1) then
1783       ztau(jl,klev+1)= - gklift*zrho(jl,klev+1)*2.*romega*
1784c    *                 (2*pstd(jl)+pmea(jl))*
1785     *                 2*pstd(jl)*
1786     *                 sin(rpi/180.*plat(jl))*pvlow(jl)
1787       ztav(jl,klev+1)=   gklift*zrho(jl,klev+1)*2.*romega*
1788c    *                 (2*pstd(jl)+pmea(jl))*
1789     *                 2*pstd(jl)*
1790     *                 sin(rpi/180.*plat(jl))*pulow(jl)
1791      else
1792       ztau(jl,klev+1)=0.0
1793       ztav(jl,klev+1)=0.0
1794      endif
1795301   continue
1796
1797C
1798C*         4.      COMPUTE LIFT PROFILE         
1799C*                 --------------------   
1800C
1801
1802  400 continue
1803
1804      do 401 jk=1,klev
1805      do 401 jl=kidia,kfdia
1806      if(ktest(jl).eq.1) then
1807      ztau(jl,jk)=ztau(jl,klev+1)*paphm1(jl,jk)/paphm1(jl,klev+1)
1808      ztav(jl,jk)=ztav(jl,klev+1)*paphm1(jl,jk)/paphm1(jl,klev+1)
1809      else
1810      ztau(jl,jk)=0.0
1811      ztav(jl,jk)=0.0
1812      endif
1813401   continue
1814C
1815C
1816C*         5.      COMPUTE TENDENCIES.
1817C*                 -------------------
1818      if(lifthigh)then
1819C
1820  500 continue
1821C
1822C  EXPLICIT SOLUTION AT ALL LEVELS
1823C
1824      do 524 jk=1,klev
1825      do 523 jl=kidia,kfdia
1826      if(ktest(jl).eq.1) then
1827      zdelp=paphm1(jl,jk+1)-paphm1(jl,jk)
1828      zdudt(jl)=-rg*(ztau(jl,jk+1)-ztau(jl,jk))/zdelp
1829      zdvdt(jl)=-rg*(ztav(jl,jk+1)-ztav(jl,jk))/zdelp
1830      endif 
1831  523 continue
1832  524 continue
1833C
1834C  PROJECT PERPENDICULARLY TO U NOT TO DESTROY ENERGY
1835C
1836      do 530 jk=1,klev
1837      do 530 jl=kidia,kfdia
1838      if(ktest(jl).eq.1) then
1839
1840        zslow=sqrt(pulow(jl)**2+pvlow(jl)**2)
1841        zsqua=amax1(sqrt(pum1(jl,jk)**2+pvm1(jl,jk)**2),gvsec)
1842        zscav=-zdudt(jl)*pvm1(jl,jk)+zdvdt(jl)*pum1(jl,jk)
1843        if(zsqua.gt.gvsec)then
1844          pvom(jl,jk)=-zscav*pvm1(jl,jk)/zsqua**2
1845          pvol(jl,jk)= zscav*pum1(jl,jk)/zsqua**2
1846        else
1847          pvom(jl,jk)=0.0
1848          pvol(jl,jk)=0.0     
1849        endif 
1850        zsqua=sqrt(pum1(jl,jk)**2+pum1(jl,jk)**2)               
1851        if(zsqua.lt.zslow)then
1852          pvom(jl,jk)=zsqua/zslow*pvom(jl,jk)
1853          pvol(jl,jk)=zsqua/zslow*pvol(jl,jk)
1854        endif
1855
1856      endif 
1857530   continue
1858C
1859C  6.  LOW LEVEL LIFT, SEMI IMPLICIT:
1860C  ----------------------------------
1861
1862      else
1863
1864        do 601 jl=kidia,kfdia
1865        if(ktest(jl).eq.1) then
1866          do jk=klev,iknub(jl),-1
1867          zbet=gklift*2.*romega*sin(rpi/180.*plat(jl))*ztmst*
1868     *        (pgeom1(jl,iknub(jl)-1)-pgeom1(jl,  jk))/
1869     *        (pgeom1(jl,iknub(jl)-1)-pgeom1(jl,klev))
1870          zdudt(jl)=-pum1(jl,jk)/ztmst/(1+zbet**2)
1871          zdvdt(jl)=-pvm1(jl,jk)/ztmst/(1+zbet**2)
1872          pvom(jl,jk)= zbet**2*zdudt(jl) - zbet   *zdvdt(jl)
1873          pvol(jl,jk)= zbet   *zdudt(jl) + zbet**2*zdvdt(jl)   
1874          enddo
1875        endif
1876 601    continue
1877
1878      endif
1879
1880c     print *,' out orolift'
1881
1882      return
1883      end
1884      SUBROUTINE SUGWD_strato(NLON,NLEV,paprs,pplay)
1885C     
1886C
1887C**** *SUGWD* INITIALIZE COMMON YOEGWD CONTROLLING GRAVITY WAVE DRAG
1888C
1889C     PURPOSE.
1890C     --------
1891C           INITIALIZE YOEGWD, THE COMMON THAT CONTROLS THE
1892C           GRAVITY WAVE DRAG PARAMETRIZATION.
1893C    VERY IMPORTANT:
1894C    ______________
1895C           THIS ROUTINE SET_UP THE "TUNABLE PARAMETERS" OF THE
1896C           VARIOUS SSO SCHEMES
1897C
1898C**   INTERFACE.
1899C     ----------
1900C        CALL *SUGWD* FROM *SUPHEC*
1901C              -----        ------
1902C
1903C        EXPLICIT ARGUMENTS :
1904C        --------------------
1905C        PAPRS,PPLAY : Pressure at semi and full model levels
1906C        NLEV        : number of model levels
1907c        NLON        : number of points treated in the physics
1908C
1909C        IMPLICIT ARGUMENTS :
1910C        --------------------
1911C        COMMON YOEGWD
1912C-GFRCRIT-R:  Critical Non-dimensional mountain Height
1913C             (HNC in (1),    LOTT 1999)
1914C-GKWAKE--R:  Bluff-body drag coefficient for low level wake
1915C             (Cd in (2),     LOTT 1999)
1916C-GRCRIT--R:  Critical Richardson Number
1917C             (Ric, End of first column p791 of LOTT 1999)
1918C-GKDRAG--R:  Gravity wave drag coefficient
1919C             (G in (3),      LOTT 1999)
1920C-GKLIFT--R:  Mountain Lift coefficient
1921C             (Cl in (4),     LOTT 1999)
1922C-GHMAX---R:  Not used
1923C-GRAHILO-R:  Set-up the trapped waves fraction
1924C             (Beta , End of first column,  LOTT 1999)
1925C
1926C-GSIGCR--R:  Security value for blocked flow depth
1927C-NKTOPG--I:  Security value for blocked flow level
1928C-nstra----I:  An estimate to qualify the upper levels of
1929C             the model where one wants to impose strees
1930C             profiles
1931C-GSSECC--R:  Security min value for low-level B-V frequency
1932C-GTSEC---R:  Security min value for anisotropy and GW stress.
1933C-GVSEC---R:  Security min value for ulow
1934C         
1935C
1936C     METHOD.
1937C     -------
1938C        SEE DOCUMENTATION
1939C
1940C     EXTERNALS.
1941C     ----------
1942C        NONE
1943C
1944C     REFERENCE.
1945C     ----------
1946C     Lott, 1999: Alleviation of stationary biases in a GCM through...
1947C                 Monthly Weather Review, 127, pp 788-801.
1948C
1949C     AUTHOR.
1950C     -------
1951C        FRANCOIS LOTT        *LMD*
1952C
1953C     MODIFICATIONS.
1954C     --------------
1955C        ORIGINAL : 90-01-01 (MARTIN MILLER, ECMWF)
1956C        LAST:  99-07-09     (FRANCOIS LOTT,LMD)
1957C     ------------------------------------------------------------------
1958      USE dimphy
1959      USE mod_phys_lmdz_para
1960      USE mod_grid_phy_lmdz
1961      IMPLICIT NONE
1962C
1963C     -----------------------------------------------------------------
1964#include "YOEGWD.h"
1965C      ----------------------------------------------------------------
1966C
1967C  ARGUMENTS
1968      integer nlon,nlev
1969      REAL paprs(nlon,nlev+1)
1970      REAL pplay(nlon,nlev)
1971C
1972      INTEGER JK
1973      REAL ZPR,ZTOP,ZSIGT,ZPM1R
1974      REAL :: pplay_glo(klon_glo,nlev)
1975      REAL :: paprs_glo(klon_glo,nlev+1)
1976
1977C
1978C*       1.    SET THE VALUES OF THE PARAMETERS
1979C              --------------------------------
1980C
1981 100  CONTINUE
1982C
1983      PRINT *,' DANS SUGWD NLEV=',NLEV
1984      GHMAX=10000.
1985C
1986      ZPR=100000.
1987      ZTOP=0.001
1988      ZSIGT=0.94
1989cold  ZPR=80000.
1990cold  ZSIGT=0.85
1991C
1992      CALL gather(pplay,pplay_glo)
1993      CALL bcast(pplay_glo)
1994      CALL gather(paprs,paprs_glo)
1995      CALL bcast(paprs_glo)
1996
1997      DO 110 JK=1,NLEV
1998      ZPM1R=pplay_glo(klon_glo/2,jk)/paprs_glo(klon_glo/2,1)
1999      IF(ZPM1R.GE.ZSIGT)THEN
2000         nktopg=JK
2001      ENDIF
2002      ZPM1R=pplay_glo(klon_glo/2,jk)/paprs_glo(klon_glo/2,1)
2003      IF(ZPM1R.GE.ZTOP)THEN
2004         nstra=JK
2005      ENDIF
2006  110 CONTINUE
2007c
2008c  inversion car dans orodrag on compte les niveaux a l'envers
2009      nktopg=nlev-nktopg+1
2010      nstra=nlev-nstra
2011      print *,' DANS SUGWD nktopg=', nktopg
2012      print *,' DANS SUGWD nstra=', nstra
2013C
2014      GSIGCR=0.80
2015C
2016      GKDRAG=0.1875
2017      GRAHILO=0.1   
2018      GRCRIT=1.00
2019      GFRCRIT=1.00
2020      GKWAKE=0.50
2021C
2022      GKLIFT=0.25
2023      GVCRIT =0.1
2024
2025      WRITE(UNIT=6,FMT='('' *** SSO essential constants ***'')')
2026      WRITE(UNIT=6,FMT='('' *** SPECIFIED IN SUGWD ***'')')
2027      WRITE(UNIT=6,FMT='('' Gravity wave ct '',E13.7,'' '')')GKDRAG
2028      WRITE(UNIT=6,FMT='('' Trapped/total wave dag '',E13.7,'' '')')
2029     S      GRAHILO
2030      WRITE(UNIT=6,FMT='('' Critical Richardson   = '',E13.7,'' '')')
2031     S                  GRCRIT
2032      WRITE(UNIT=6,FMT='('' Critical Froude'',e13.7)') GFRCRIT
2033      WRITE(UNIT=6,FMT='('' Low level Wake bluff cte'',e13.7)') GKWAKE
2034      WRITE(UNIT=6,FMT='('' Low level lift  cte'',e13.7)') GKLIFT
2035
2036C
2037C
2038C      ----------------------------------------------------------------
2039C
2040C*       2.    SET VALUES OF SECURITY PARAMETERS
2041C              ---------------------------------
2042C
2043 200  CONTINUE
2044C
2045      GVSEC=0.10
2046      GSSEC=0.0001
2047C
2048      GTSEC=0.00001
2049C
2050      RETURN
2051      END
2052
Note: See TracBrowser for help on using the repository browser.