phylmd/Orography/gwprofil.f

module gwprofil_m

  IMPLICIT NONE

contains

  SUBROUTINE gwprofil(nlon, nlev, ktest, kkcrith, kcrit, paphm1, prho, pstab, &
       pvph, pri, ptau, pdmod, psig, pvar)

    ! Method. The stress profile for gravity waves is computed as
    ! follows: it is constant (no gwd) at the levels between the
    ! ground and the top of the blocked layer (kkenvh).  It decreases
    ! linearly with height from the top of the blocked layer to 3 *
    ! varor (kknu), to simulate lee waves or nonlinear gravity wave
    ! breaking. Above it is constant, except when the wave encounters
    ! a critical level (kcrit) or when it breaks.

    ! Reference. See ECMWF research department documentation of the
    ! "I.F.S."

    ! Modifications. Passage of the new gwdrag TO I.F.S. (F. LOTT,
    ! 22/11/93)

    USE dimphy, ONLY : klev, klon
    USE yoegwd, ONLY : gkdrag, grahilo, grcrit, gssec, gtsec, nstra

    INTEGER, intent(in):: nlon, nlev
    INTEGER, intent(in):: ktest(nlon), kkcrith(nlon), kcrit(nlon)
    REAL, intent(in):: paphm1(nlon, nlev+1), prho(nlon, nlev+1)
    REAL, intent(in):: pstab(nlon, nlev+1)
    real, intent(in):: pvph(nlon, nlev+1), pri(nlon, nlev+1)
    real ptau(nlon, nlev+1)
    REAL, intent(in):: pdmod(nlon)
    REAL, INTENT (IN) :: psig(nlon)
    REAL, INTENT (IN) :: pvar(nlon)

    ! Local:
    INTEGER ilevh, jl, jk
    REAL zsqr, zalfa, zriw, zdel, zb, zalpha, zdz2n
    REAL zdelp, zdelpt
    REAL zdz2(klon, klev), znorm(klon), zoro(klon)
    REAL ztau(klon, klev+1)

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

    ! 1. INITIALIZATION

    ! COMPUTATIONAL CONSTANTS.

    ilevh = klev/3

    DO jl = 1, klon
       IF (ktest(jl)==1) THEN
          zoro(jl) = psig(jl)*pdmod(jl)/4./max(pvar(jl), 1.0)
          ztau(jl, klev+1) = ptau(jl, klev+1)
       END IF
    end DO

    DO jk = klev, 2, -1
       ! 4.1 CONSTANT WAVE STRESS UNTIL TOP OF THE
       ! BLOCKING LAYER.
       DO jl = 1, klon
          IF (ktest(jl)==1) THEN
             IF (jk>kkcrith(jl)) THEN
                ptau(jl, jk) = ztau(jl, klev+1)
             ELSE
                ptau(jl, jk) = grahilo*ztau(jl, klev+1)
             END IF
          END IF
       end DO

       ! 4.2 WAVE DISPLACEMENT AT NEXT LEVEL.
       DO jl = 1, klon
          IF (ktest(jl)==1) THEN
             IF (jk<kkcrith(jl)) THEN
                znorm(jl) = gkdrag * prho(jl, jk) * sqrt(pstab(jl, jk)) &
                     * pvph(jl, jk)* zoro(jl)
                zdz2(jl, jk) = ptau(jl, jk+1)/max(znorm(jl), gssec)
             END IF
          END IF
       end DO

       ! 4.3 WAVE RICHARDSON NUMBER, NEW WAVE DISPLACEMENT
       ! AND STRESS: BREAKING EVALUATION AND CRITICAL
       ! LEVEL
       DO jl = 1, klon
          IF (ktest(jl)==1) THEN
             IF (jk<kkcrith(jl)) THEN
                IF ((ptau(jl, jk+1)<gtsec) .OR. (jk<=kcrit(jl))) THEN
                   ptau(jl, jk) = 0.0
                ELSE
                   zsqr = sqrt(pri(jl, jk))
                   zalfa = sqrt(pstab(jl, jk)*zdz2(jl, jk))/pvph(jl, jk)
                   zriw = pri(jl, jk)*(1.-zalfa)/(1+zalfa*zsqr)**2
                   IF (zriw<grcrit) THEN
                      zdel = 4./zsqr/grcrit + 1./grcrit**2 + 4./grcrit
                      zb = 1./grcrit + 2./zsqr
                      zalpha = 0.5*(-zb+sqrt(zdel))
                      zdz2n = (pvph(jl, jk)*zalpha)**2/pstab(jl, jk)
                      ptau(jl, jk) = znorm(jl)*zdz2n
                   ELSE
                      ptau(jl, jk) = znorm(jl)*zdz2(jl, jk)
                   END IF
                   ptau(jl, jk) = min(ptau(jl, jk), ptau(jl, jk+1))
                END IF
             END IF
          END IF
       end DO
    end DO

    ! REORGANISATION OF THE STRESS PROFILE AT LOW LEVEL

    DO jl = 1, klon
       IF (ktest(jl)==1) THEN
          ztau(jl, kkcrith(jl)) = ptau(jl, kkcrith(jl))
          ztau(jl, nstra) = ptau(jl, nstra)
       END IF
    end DO

    DO jk = 1, klev
       DO jl = 1, klon
          IF (ktest(jl)==1) THEN
             IF (jk>kkcrith(jl)) THEN
                zdelp = paphm1(jl, jk) - paphm1(jl, klev+1)
                zdelpt = paphm1(jl, kkcrith(jl)) - paphm1(jl, klev+1)
                ptau(jl, jk) = ztau(jl, klev+1) &
                     + (ztau(jl, kkcrith(jl)) - ztau(jl, klev+1))*zdelp/zdelpt
             END IF
          END IF
       end DO

       ! REORGANISATION IN THE STRATOSPHERE
       DO jl = 1, klon
          IF (ktest(jl)==1) THEN
             IF (jk < nstra) THEN
                zdelp = paphm1(jl, nstra)
                zdelpt = paphm1(jl, jk)
                ptau(jl, jk) = ztau(jl, nstra) * zdelpt / zdelp
             END IF
          END IF
       end DO

       ! REORGANISATION IN THE TROPOSPHERE
       DO jl = 1, klon
          IF (ktest(jl)==1) THEN
             IF (jk<kkcrith(jl) .AND. jk > nstra) THEN
                zdelp = paphm1(jl, jk) - paphm1(jl, kkcrith(jl))
                zdelpt = paphm1(jl, nstra) - paphm1(jl, kkcrith(jl))
                ptau(jl, jk) = ztau(jl, kkcrith(jl)) &
                     + (ztau(jl, nstra) - ztau(jl, kkcrith(jl)))*zdelp/zdelpt
             END IF
          END IF
       end DO
    end DO

  END SUBROUTINE gwprofil

end module gwprofil_m
1	guez	54	module gwprofil_m
2	guez	23
3	guez	54	IMPLICIT NONE
4	guez	23
5	guez	54	contains
6	guez	23
7	guez	158	SUBROUTINE gwprofil(nlon, nlev, ktest, kkcrith, kcrit, paphm1, prho, pstab, &
8			pvph, pri, ptau, pdmod, psig, pvar)
9	guez	23
10	guez	54	! Method. The stress profile for gravity waves is computed as
11			! follows: it is constant (no gwd) at the levels between the
12			! ground and the top of the blocked layer (kkenvh). It decreases
13	guez	158	! linearly with height from the top of the blocked layer to 3 *
14			! varor (kknu), to simulate lee waves or nonlinear gravity wave
15	guez	54	! breaking. Above it is constant, except when the wave encounters
16			! a critical level (kcrit) or when it breaks.
17	guez	23
18	guez	158	! Reference. See ECMWF research department documentation of the
19			! "I.F.S."
20	guez	23
21	guez	158	! Modifications. Passage of the new gwdrag TO I.F.S. (F. LOTT,
22			! 22/11/93)
23	guez	23
24	guez	54	USE dimphy, ONLY : klev, klon
25			USE yoegwd, ONLY : gkdrag, grahilo, grcrit, gssec, gtsec, nstra
26	guez	23
27	guez	158	INTEGER, intent(in):: nlon, nlev
28			INTEGER, intent(in):: ktest(nlon), kkcrith(nlon), kcrit(nlon)
29			REAL, intent(in):: paphm1(nlon, nlev+1), prho(nlon, nlev+1)
30			REAL, intent(in):: pstab(nlon, nlev+1)
31			real, intent(in):: pvph(nlon, nlev+1), pri(nlon, nlev+1)
32			real ptau(nlon, nlev+1)
33			REAL, intent(in):: pdmod(nlon)
34	guez	54	REAL, INTENT (IN) :: psig(nlon)
35			REAL, INTENT (IN) :: pvar(nlon)
36	guez	23
37	guez	158	! Local:
38	guez	150	INTEGER ilevh, jl, jk
39	guez	54	REAL zsqr, zalfa, zriw, zdel, zb, zalpha, zdz2n
40			REAL zdelp, zdelpt
41			REAL zdz2(klon, klev), znorm(klon), zoro(klon)
42			REAL ztau(klon, klev+1)
43	guez	23
44	guez	54	!-----------------------------------------------------------------------
45	guez	23
46	guez	54	! 1. INITIALIZATION
47	guez	23
48	guez	54	! COMPUTATIONAL CONSTANTS.
49	guez	23
50	guez	54	ilevh = klev/3
51	guez	23
52	guez	54	DO jl = 1, klon
53			IF (ktest(jl)==1) THEN
54			zoro(jl) = psig(jl)*pdmod(jl)/4./max(pvar(jl), 1.0)
55			ztau(jl, klev+1) = ptau(jl, klev+1)
56			END IF
57			end DO
58	guez	23
59	guez	54	DO jk = klev, 2, -1
60			! 4.1 CONSTANT WAVE STRESS UNTIL TOP OF THE
61			! BLOCKING LAYER.
62			DO jl = 1, klon
63	guez	23	IF (ktest(jl)==1) THEN
64	guez	54	IF (jk>kkcrith(jl)) THEN
65			ptau(jl, jk) = ztau(jl, klev+1)
66			ELSE
67			ptau(jl, jk) = grahilo*ztau(jl, klev+1)
68			END IF
69	guez	23	END IF
70	guez	54	end DO
71	guez	23
72	guez	54	! 4.2 WAVE DISPLACEMENT AT NEXT LEVEL.
73			DO jl = 1, klon
74	guez	23	IF (ktest(jl)==1) THEN
75	guez	54	IF (jk<kkcrith(jl)) THEN
76			znorm(jl) = gkdrag * prho(jl, jk) * sqrt(pstab(jl, jk)) &
77			* pvph(jl, jk)* zoro(jl)
78			zdz2(jl, jk) = ptau(jl, jk+1)/max(znorm(jl), gssec)
79			END IF
80	guez	23	END IF
81	guez	54	end DO
82	guez	23
83	guez	54	! 4.3 WAVE RICHARDSON NUMBER, NEW WAVE DISPLACEMENT
84			! AND STRESS: BREAKING EVALUATION AND CRITICAL
85			! LEVEL
86			DO jl = 1, klon
87	guez	23	IF (ktest(jl)==1) THEN
88	guez	54	IF (jk<kkcrith(jl)) THEN
89			IF ((ptau(jl, jk+1)<gtsec) .OR. (jk<=kcrit(jl))) THEN
90			ptau(jl, jk) = 0.0
91	guez	23	ELSE
92	guez	54	zsqr = sqrt(pri(jl, jk))
93			zalfa = sqrt(pstab(jl, jk)*zdz2(jl, jk))/pvph(jl, jk)
94			zriw = pri(jl, jk)(1.-zalfa)/(1+zalfazsqr)**2
95			IF (zriw<grcrit) THEN
96			zdel = 4./zsqr/grcrit + 1./grcrit**2 + 4./grcrit
97			zb = 1./grcrit + 2./zsqr
98			zalpha = 0.5*(-zb+sqrt(zdel))
99			zdz2n = (pvph(jl, jk)zalpha)*2/pstab(jl, jk)
100			ptau(jl, jk) = znorm(jl)*zdz2n
101			ELSE
102			ptau(jl, jk) = znorm(jl)*zdz2(jl, jk)
103			END IF
104			ptau(jl, jk) = min(ptau(jl, jk), ptau(jl, jk+1))
105	guez	23	END IF
106	guez	54	END IF
107	guez	23	END IF
108	guez	54	end DO
109			end DO
110	guez	23
111	guez	54	! REORGANISATION OF THE STRESS PROFILE AT LOW LEVEL
112	guez	23
113	guez	54	DO jl = 1, klon
114			IF (ktest(jl)==1) THEN
115			ztau(jl, kkcrith(jl)) = ptau(jl, kkcrith(jl))
116			ztau(jl, nstra) = ptau(jl, nstra)
117			END IF
118			end DO
119	guez	23
120	guez	54	DO jk = 1, klev
121			DO jl = 1, klon
122	guez	23	IF (ktest(jl)==1) THEN
123	guez	54	IF (jk>kkcrith(jl)) THEN
124			zdelp = paphm1(jl, jk) - paphm1(jl, klev+1)
125			zdelpt = paphm1(jl, kkcrith(jl)) - paphm1(jl, klev+1)
126			ptau(jl, jk) = ztau(jl, klev+1) &
127			+ (ztau(jl, kkcrith(jl)) - ztau(jl, klev+1))*zdelp/zdelpt
128			END IF
129	guez	23	END IF
130	guez	54	end DO
131	guez	23
132	guez	54	! REORGANISATION IN THE STRATOSPHERE
133			DO jl = 1, klon
134	guez	23	IF (ktest(jl)==1) THEN
135	guez	54	IF (jk < nstra) THEN
136			zdelp = paphm1(jl, nstra)
137			zdelpt = paphm1(jl, jk)
138			ptau(jl, jk) = ztau(jl, nstra) * zdelpt / zdelp
139			END IF
140	guez	23	END IF
141	guez	54	end DO
142	guez	23
143	guez	54	! REORGANISATION IN THE TROPOSPHERE
144			DO jl = 1, klon
145	guez	23	IF (ktest(jl)==1) THEN
146	guez	54	IF (jk<kkcrith(jl) .AND. jk > nstra) THEN
147			zdelp = paphm1(jl, jk) - paphm1(jl, kkcrith(jl))
148			zdelpt = paphm1(jl, nstra) - paphm1(jl, kkcrith(jl))
149			ptau(jl, jk) = ztau(jl, kkcrith(jl)) &
150			+ (ztau(jl, nstra) - ztau(jl, kkcrith(jl)))*zdelp/zdelpt
151			END IF
152	guez	23	END IF
153	guez	54	end DO
154			end DO
155	guez	23
156	guez	54	END SUBROUTINE gwprofil
157	guez	23
158	guez	54	end module gwprofil_m