phylmd/Orography/orodrag.f

module orodrag_m

  IMPLICIT NONE

contains

  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 gwstress_m, only: gwstress
    USE suphec_m
    USE yoegwd
    use gwprofil_m, only: gwprofil
    use orosetup_m, only: orosetup

    !**** *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.

    ! 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
    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), zvidis(klon), &
         znu(klon), zd1(klon), zd2(klon), zdmod(klon)
    REAL ztmst
    REAL, INTENT (IN) :: ptsphy

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

    !* 1. initialization

    !* 1.1 computational constants

    ztmst = ptsphy

    !* 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, 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, ikenvh, 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 jl = 1, klon
       zvidis(jl) = 0.0
       zdudt(jl) = 0.0
       zdvdt(jl) = 0.0
    end DO

    ilevp1 = klev + 1

    DO jk = 1, klev

       ! Modif vectorisation 02/04/2004
       DO 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)
             zvidis(ji) = zvidis(ji) + zdis*zdelp

             ! ENCORE UN TRUC POUR EVITER LES EXPLOSIONS

             pte(ji, jk) = 0.0

          END IF
       end DO

    end DO

    RETURN
  END SUBROUTINE orodrag

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