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

Last change on this file since 3900 was 1403, checked in by Laurent Fairhead, 14 years ago

Merged LMDZ4V5.0-dev branch changes r1292:r1399 to trunk.

Validation:
Validation consisted in compiling the HEAD revision of the trunk,
LMDZ4V5.0-dev branch and the merged sources and running different
configurations on local and SX8 machines comparing results.

Local machine: bench configuration, 32x24x11, gfortran

  • IPSLCM5A configuration (comparison between trunk and merged sources):
    • numerical convergence on dynamical fields over 3 days
    • start files are equivalent (except for RN and PB fields)
    • daily history files equivalent
  • MH07 configuration, new physics package (comparison between LMDZ4V5.0-dev branch and merged sources):
    • numerical convergence on dynamical fields over 3 days
    • start files are equivalent (except for RN and PB fields)
    • daily history files equivalent

SX8 machine (brodie), 96x95x39 on 4 processors:

  • IPSLCM5A configuration:
    • start files are equivalent (except for RN and PB fields)
    • monthly history files equivalent
  • MH07 configuration:
    • start files are equivalent (except for RN and PB fields)
    • monthly history files equivalent

Changes to the makegcm and create_make_gcm scripts to take into account
main programs in F90 files

Fusion de la branche LMDZ4V5.0-dev (r1292:r1399) au tronc principal

Validation:
La validation a consisté à compiler la HEAD de le trunk et de la banche
LMDZ4V5.0-dev et les sources fusionnées et de faire tourner le modéle selon
différentes configurations en local et sur SX8 et de comparer les résultats

En local: 32x24x11, config bench/gfortran

  • pour une config IPSLCM5A (comparaison tronc/fusion):
    • convergence numérique sur les champs dynamiques après 3 jours
    • restart et restartphy égaux (à part sur RN et Pb)
    • fichiers histoire égaux
  • pour une config nlle physique (MH07) (comparaison LMDZ4v5.0-dev/fusion):
    • convergence numérique sur les champs dynamiques après 3 jours
    • restart et restartphy égaux
    • fichiers histoire équivalents

Sur brodie, 96x95x39 sur 4 proc:

  • pour une config IPSLCM5A:
    • restart et restartphy égaux (à part sur RN et PB)
    • pas de différence dans les fichiers histmth.nc
  • pour une config MH07
    • restart et restartphy égaux (à part sur RN et PB)
    • pas de différence dans les fichiers histmth.nc

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