phylmd/Orography/orodrag.f

    SUBROUTINE orodrag(nlon,nlev,ktest,ptsphy,paphm1,papm1,pgeom1, &
        ptm1,pum1,pvm1,pmea,pstd,psig,pgamma,ptheta,ppic,pval &
        ,pulow,pvlow,pvom,pvol,pte)

      USE dimens_m
      USE dimphy
      USE suphec_m
      USE yoegwd
      use gwprofil_m, only: gwprofil
      IMPLICIT NONE


!**** *gwdrag* - does the gravity wave parametrization.

!     purpose.
!     --------

!          this routine computes the physical tendencies of the
!     prognostic variables u,v  and t due to  vertical transports by
!     subgridscale orographically excited gravity waves

!**   interface.
!     ----------
!          called from *callpar*.

!          the routine takes its input from the long-term storage:
!          u,v,t and p at t-1.

!        explicit arguments :
!        --------------------
!     ==== inputs ===
!     ==== outputs ===

!        implicit arguments :   none
!        --------------------

!      implicit logical (l)

!     method.
!     -------

!     externals.
!     ----------
      INTEGER ismin, ismax
      EXTERNAL ismin, ismax

!     reference.
!     ----------

!     author.
!     -------
!     m.miller + b.ritter   e.c.m.w.f.     15/06/86.

!     f.lott + m. miller    e.c.m.w.f.     22/11/94
!-----------------------------------------------------------------------


!-----------------------------------------------------------------------

!*       0.1   arguments
!              ---------


      INTEGER nlon, nlev, klevm1
      INTEGER jl, ilevp1, jk, ji
      REAL zdelp, ztemp, zforc, ztend
      REAL rover, zb, zc, zconb, zabsv
      REAL zzd1, ratio, zbet, zust, zvst, zdis
      REAL pte(nlon,nlev), pvol(nlon,nlev), pvom(nlon,nlev), pulow(klon), &
        pvlow(klon)
      REAL pum1(nlon,nlev), pvm1(nlon,nlev), ptm1(nlon,nlev), pmea(nlon)
      REAL, INTENT (IN) :: pstd(nlon)
      REAL, INTENT (IN) :: psig(nlon)
      REAL pgamma(nlon), ptheta(nlon), ppic(nlon), pval(nlon), &
        pgeom1(nlon,nlev), papm1(nlon,nlev), paphm1(nlon,nlev+1)

      INTEGER ktest(nlon)
!-----------------------------------------------------------------------

!*       0.2   local arrays
!              ------------
      INTEGER icrit(klon), ikcrith(klon), ikenvh(klon), &
        iknu(klon), iknu2(klon), ikcrit(klon)

      REAL ztau(klon,klev+1), zstab(klon,klev+1), &
        zvph(klon,klev+1), zrho(klon,klev+1), zri(klon,klev+1), &
        zpsi(klon,klev+1), zzdep(klon,klev)
      REAL zdudt(klon), zdvdt(klon), zdtdt(klon), zdedt(klon), zvidis(klon), &
        znu(klon), zd1(klon), zd2(klon), zdmod(klon)
      REAL ztmst, zrtmst
      REAL, INTENT (IN) :: ptsphy

!------------------------------------------------------------------

!*         1.    initialization
!                --------------

!*         1.1   computational constants
!                -----------------------

      klevm1 = klev - 1
      ztmst = ptsphy
      zrtmst = 1./ztmst
!     ------------------------------------------------------------------

!*         1.3   check whether row contains point for printing
!                ---------------------------------------------

!*         2.     precompute basic state variables.
!*                ---------- ----- ----- ----------
!*                define low level wind, project winds in plane of
!*                low level wind, determine sector in which to take
!*                the variance and set indicator for critical levels.


      CALL orosetup(nlon,ktest,ikcrit,ikcrith,icrit,ikenvh,iknu,iknu2,paphm1, &
        papm1,pum1,pvm1,ptm1,pgeom1,pstd,zrho,zri,zstab,ztau,zvph,zpsi,zzdep, &
        pulow,pvlow,ptheta,pgamma,pmea,ppic,pval,znu,zd1,zd2,zdmod)


!***********************************************************


!*         3.      compute low level stresses using subcritical and
!*                 supercritical forms.computes anisotropy coefficient
!*                 as measure of orographic twodimensionality.

      CALL gwstress(nlon,nlev,ktest,icrit,ikenvh,iknu,zrho,zstab,zvph,pstd, &
        psig,pmea,ppic,ztau,pgeom1,zdmod)


!*         4.      compute stress profile.
!*                 ------- ------ --------

      CALL gwprofil(nlon,nlev,ktest,ikcrith,icrit,paphm1,zrho,zstab, &
        zvph,zri,ztau,zdmod,psig,pstd)


!*         5.      compute tendencies.
!*                 -------------------

!  explicit solution at all levels for the gravity wave
!  implicit solution for the blocked levels

      DO 510 jl = 1, klon
        zvidis(jl) = 0.0
        zdudt(jl) = 0.0
        zdvdt(jl) = 0.0
        zdtdt(jl) = 0.0
510   CONTINUE

      ilevp1 = klev + 1


      DO 524 jk = 1, klev


!  Modif vectorisation 02/04/2004
        DO 523 ji = 1, klon
          IF (ktest(ji)==1) THEN

            zdelp = paphm1(ji,jk+1) - paphm1(ji,jk)
            ztemp = -rg*(ztau(ji,jk+1)-ztau(ji,jk))/(zvph(ji,ilevp1)*zdelp)
            zdudt(ji) = (pulow(ji)*zd1(ji)-pvlow(ji)*zd2(ji))*ztemp/zdmod(ji)
            zdvdt(ji) = (pvlow(ji)*zd1(ji)+pulow(ji)*zd2(ji))*ztemp/zdmod(ji)

! controle des overshoots:

            zforc = sqrt(zdudt(ji)**2+zdvdt(ji)**2) + 1.E-12
            ztend = sqrt(pum1(ji,jk)**2+pvm1(ji,jk)**2)/ztmst + 1.E-12
            rover = 0.25
            IF (zforc>=rover*ztend) THEN
              zdudt(ji) = rover*ztend/zforc*zdudt(ji)
              zdvdt(ji) = rover*ztend/zforc*zdvdt(ji)
            END IF

! fin du controle des overshoots

            IF (jk>=ikenvh(ji)) THEN
              zb = 1.0 - 0.18*pgamma(ji) - 0.04*pgamma(ji)**2
              zc = 0.48*pgamma(ji) + 0.3*pgamma(ji)**2
              zconb = 2.*ztmst*gkwake*psig(ji)/(4.*pstd(ji))
              zabsv = sqrt(pum1(ji,jk)**2+pvm1(ji,jk)**2)/2.
              zzd1 = zb*cos(zpsi(ji,jk))**2 + zc*sin(zpsi(ji,jk))**2
              ratio = (cos(zpsi(ji,jk))**2+pgamma(ji)*sin(zpsi(ji, &
                jk))**2)/(pgamma(ji)*cos(zpsi(ji,jk))**2+sin(zpsi(ji,jk))**2)
              zbet = max(0.,2.-1./ratio)*zconb*zzdep(ji,jk)*zzd1*zabsv

! simplement oppose au vent

              zdudt(ji) = -pum1(ji,jk)/ztmst
              zdvdt(ji) = -pvm1(ji,jk)/ztmst

!  projection dans la direction de l'axe principal de l'orographie
!mod     zdudt(ji)=-(pum1(ji,jk)*cos(ptheta(ji)*rpi/180.)
!mod *              +pvm1(ji,jk)*sin(ptheta(ji)*rpi/180.))
!mod *              *cos(ptheta(ji)*rpi/180.)/ztmst
!mod     zdvdt(ji)=-(pum1(ji,jk)*cos(ptheta(ji)*rpi/180.)
!mod *              +pvm1(ji,jk)*sin(ptheta(ji)*rpi/180.))
!mod *              *sin(ptheta(ji)*rpi/180.)/ztmst
              zdudt(ji) = zdudt(ji)*(zbet/(1.+zbet))
              zdvdt(ji) = zdvdt(ji)*(zbet/(1.+zbet))
            END IF
            pvom(ji,jk) = zdudt(ji)
            pvol(ji,jk) = zdvdt(ji)
            zust = pum1(ji,jk) + ztmst*zdudt(ji)
            zvst = pvm1(ji,jk) + ztmst*zdvdt(ji)
            zdis = 0.5*(pum1(ji,jk)**2+pvm1(ji,jk)**2-zust**2-zvst**2)
            zdedt(ji) = zdis/ztmst
            zvidis(ji) = zvidis(ji) + zdis*zdelp
            zdtdt(ji) = zdedt(ji)/rcpd
!     pte(ji,jk)=zdtdt(ji)

!  ENCORE UN TRUC POUR EVITER LES EXPLOSIONS

            pte(ji,jk) = 0.0

          END IF
523     CONTINUE

524   CONTINUE


      RETURN
    END
1	guez	150	SUBROUTINE orodrag(nlon,nlev,ktest,ptsphy,paphm1,papm1,pgeom1, &
2	guez	23	ptm1,pum1,pvm1,pmea,pstd,psig,pgamma,ptheta,ppic,pval &
3			,pulow,pvlow,pvom,pvol,pte)
4
5			USE dimens_m
6			USE dimphy
7	guez	38	USE suphec_m
8	guez	23	USE yoegwd
9	guez	54	use gwprofil_m, only: gwprofil
10	guez	23	IMPLICIT NONE
11
12
13
14			!**** gwdrag - does the gravity wave parametrization.
15
16			! purpose.
17			! --------
18
19			! this routine computes the physical tendencies of the
20			! prognostic variables u,v and t due to vertical transports by
21			! subgridscale orographically excited gravity waves
22
23			!** interface.
24			! ----------
25			! called from callpar.
26
27			! the routine takes its input from the long-term storage:
28			! u,v,t and p at t-1.
29
30			! explicit arguments :
31			! --------------------
32			! ==== inputs ===
33			! ==== outputs ===
34
35			! implicit arguments : none
36			! --------------------
37
38			! implicit logical (l)
39
40			! method.
41			! -------
42
43			! externals.
44			! ----------
45			INTEGER ismin, ismax
46			EXTERNAL ismin, ismax
47
48			! reference.
49			! ----------
50
51			! author.
52			! -------
53			! m.miller + b.ritter e.c.m.w.f. 15/06/86.
54
55			! f.lott + m. miller e.c.m.w.f. 22/11/94
56			!-----------------------------------------------------------------------
57
58
59			!-----------------------------------------------------------------------
60
61			!* 0.1 arguments
62			! ---------
63
64
65			INTEGER nlon, nlev, klevm1
66	guez	150	INTEGER jl, ilevp1, jk, ji
67	guez	23	REAL zdelp, ztemp, zforc, ztend
68			REAL rover, zb, zc, zconb, zabsv
69			REAL zzd1, ratio, zbet, zust, zvst, zdis
70			REAL pte(nlon,nlev), pvol(nlon,nlev), pvom(nlon,nlev), pulow(klon), &
71			pvlow(klon)
72			REAL pum1(nlon,nlev), pvm1(nlon,nlev), ptm1(nlon,nlev), pmea(nlon)
73			REAL, INTENT (IN) :: pstd(nlon)
74			REAL, INTENT (IN) :: psig(nlon)
75			REAL pgamma(nlon), ptheta(nlon), ppic(nlon), pval(nlon), &
76			pgeom1(nlon,nlev), papm1(nlon,nlev), paphm1(nlon,nlev+1)
77
78	guez	150	INTEGER ktest(nlon)
79	guez	23	!-----------------------------------------------------------------------
80
81			!* 0.2 local arrays
82			! ------------
83	guez	150	INTEGER icrit(klon), ikcrith(klon), ikenvh(klon), &
84			iknu(klon), iknu2(klon), ikcrit(klon)
85	guez	23
86	guez	150	REAL ztau(klon,klev+1), zstab(klon,klev+1), &
87	guez	23	zvph(klon,klev+1), zrho(klon,klev+1), zri(klon,klev+1), &
88			zpsi(klon,klev+1), zzdep(klon,klev)
89			REAL zdudt(klon), zdvdt(klon), zdtdt(klon), zdedt(klon), zvidis(klon), &
90			znu(klon), zd1(klon), zd2(klon), zdmod(klon)
91			REAL ztmst, zrtmst
92			REAL, INTENT (IN) :: ptsphy
93
94			!------------------------------------------------------------------
95
96			!* 1. initialization
97			! --------------
98
99			!* 1.1 computational constants
100			! -----------------------
101
102			klevm1 = klev - 1
103			ztmst = ptsphy
104			zrtmst = 1./ztmst
105			! ------------------------------------------------------------------
106
107			!* 1.3 check whether row contains point for printing
108			! ---------------------------------------------
109
110			!* 2. precompute basic state variables.
111			!* ---------- ----- ----- ----------
112			!* define low level wind, project winds in plane of
113			!* low level wind, determine sector in which to take
114			!* the variance and set indicator for critical levels.
115
116
117			CALL orosetup(nlon,ktest,ikcrit,ikcrith,icrit,ikenvh,iknu,iknu2,paphm1, &
118			papm1,pum1,pvm1,ptm1,pgeom1,pstd,zrho,zri,zstab,ztau,zvph,zpsi,zzdep, &
119			pulow,pvlow,ptheta,pgamma,pmea,ppic,pval,znu,zd1,zd2,zdmod)
120
121
122
123			!***********************************************************
124
125
126			!* 3. compute low level stresses using subcritical and
127			!* supercritical forms.computes anisotropy coefficient
128			!* as measure of orographic twodimensionality.
129
130			CALL gwstress(nlon,nlev,ktest,icrit,ikenvh,iknu,zrho,zstab,zvph,pstd, &
131			psig,pmea,ppic,ztau,pgeom1,zdmod)
132
133
134			!* 4. compute stress profile.
135			!* ------- ------ --------
136
137	guez	150	CALL gwprofil(nlon,nlev,ktest,ikcrith,icrit,paphm1,zrho,zstab, &
138	guez	23	zvph,zri,ztau,zdmod,psig,pstd)
139
140
141			!* 5. compute tendencies.
142			!* -------------------
143
144			! explicit solution at all levels for the gravity wave
145			! implicit solution for the blocked levels
146
147			DO 510 jl = 1, klon
148			zvidis(jl) = 0.0
149			zdudt(jl) = 0.0
150			zdvdt(jl) = 0.0
151			zdtdt(jl) = 0.0
152			510 CONTINUE
153
154			ilevp1 = klev + 1
155
156
157			DO 524 jk = 1, klev
158
159
160			! Modif vectorisation 02/04/2004
161			DO 523 ji = 1, klon
162			IF (ktest(ji)==1) THEN
163
164			zdelp = paphm1(ji,jk+1) - paphm1(ji,jk)
165			ztemp = -rg(ztau(ji,jk+1)-ztau(ji,jk))/(zvph(ji,ilevp1)zdelp)
166			zdudt(ji) = (pulow(ji)zd1(ji)-pvlow(ji)zd2(ji))*ztemp/zdmod(ji)
167			zdvdt(ji) = (pvlow(ji)zd1(ji)+pulow(ji)zd2(ji))*ztemp/zdmod(ji)
168
169			! controle des overshoots:
170
171			zforc = sqrt(zdudt(ji)2+zdvdt(ji)2) + 1.E-12
172			ztend = sqrt(pum1(ji,jk)2+pvm1(ji,jk)2)/ztmst + 1.E-12
173			rover = 0.25
174			IF (zforc>=rover*ztend) THEN
175			zdudt(ji) = roverztend/zforczdudt(ji)
176			zdvdt(ji) = roverztend/zforczdvdt(ji)
177			END IF
178
179			! fin du controle des overshoots
180
181			IF (jk>=ikenvh(ji)) THEN
182			zb = 1.0 - 0.18pgamma(ji) - 0.04pgamma(ji)**2
183			zc = 0.48pgamma(ji) + 0.3pgamma(ji)**2
184			zconb = 2.ztmstgkwakepsig(ji)/(4.pstd(ji))
185			zabsv = sqrt(pum1(ji,jk)2+pvm1(ji,jk)2)/2.
186			zzd1 = zbcos(zpsi(ji,jk))2 + zcsin(zpsi(ji,jk))**2
187			ratio = (cos(zpsi(ji,jk))*2+pgamma(ji)sin(zpsi(ji, &
188			jk))*2)/(pgamma(ji)cos(zpsi(ji,jk))2+sin(zpsi(ji,jk))2)
189			zbet = max(0.,2.-1./ratio)zconbzzdep(ji,jk)zzd1zabsv
190
191			! simplement oppose au vent
192
193			zdudt(ji) = -pum1(ji,jk)/ztmst
194			zdvdt(ji) = -pvm1(ji,jk)/ztmst
195
196			! projection dans la direction de l'axe principal de l'orographie
197			!mod zdudt(ji)=-(pum1(ji,jk)cos(ptheta(ji)rpi/180.)
198			!mod * +pvm1(ji,jk)sin(ptheta(ji)rpi/180.))
199			!mod * cos(ptheta(ji)rpi/180.)/ztmst
200			!mod zdvdt(ji)=-(pum1(ji,jk)cos(ptheta(ji)rpi/180.)
201			!mod * +pvm1(ji,jk)sin(ptheta(ji)rpi/180.))
202			!mod * sin(ptheta(ji)rpi/180.)/ztmst
203			zdudt(ji) = zdudt(ji)*(zbet/(1.+zbet))
204			zdvdt(ji) = zdvdt(ji)*(zbet/(1.+zbet))
205			END IF
206			pvom(ji,jk) = zdudt(ji)
207			pvol(ji,jk) = zdvdt(ji)
208			zust = pum1(ji,jk) + ztmst*zdudt(ji)
209			zvst = pvm1(ji,jk) + ztmst*zdvdt(ji)
210			zdis = 0.5(pum1(ji,jk)2+pvm1(ji,jk)2-zust2-zvst*2)
211			zdedt(ji) = zdis/ztmst
212			zvidis(ji) = zvidis(ji) + zdis*zdelp
213			zdtdt(ji) = zdedt(ji)/rcpd
214			! pte(ji,jk)=zdtdt(ji)
215
216			! ENCORE UN TRUC POUR EVITER LES EXPLOSIONS
217
218			pte(ji,jk) = 0.0
219
220			END IF
221			523 CONTINUE
222
223			524 CONTINUE
224
225
226			RETURN
227			END