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	SUBROUTINE orodrag(nlon,nlev,ktest,ptsphy,paphm1,papm1,pgeom1, &
2	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	USE suphec_m
8	USE yoegwd
9	use gwprofil_m, only: gwprofil
10	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	INTEGER jl, ilevp1, jk, ji
67	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	INTEGER ktest(nlon)
79	!-----------------------------------------------------------------------
80
81	!* 0.2 local arrays
82	! ------------
83	INTEGER icrit(klon), ikcrith(klon), ikenvh(klon), &
84	iknu(klon), iknu2(klon), ikcrit(klon)
85
86	REAL ztau(klon,klev+1), zstab(klon,klev+1), &
87	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	CALL gwprofil(nlon,nlev,ktest,ikcrith,icrit,paphm1,zrho,zstab, &
138	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