/[lmdze]/trunk/Sources/phylmd/Orography/gwprofil.f
ViewVC logotype

Contents of /trunk/Sources/phylmd/Orography/gwprofil.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 158 - (show annotations)
Tue Jul 21 14:44:45 2015 UTC (8 years, 9 months ago) by guez
File size: 5242 byte(s)
Subroutine sugwd sets variables of module yoegwd. Better to put it
into module yoegwd.

Variables of module yoegwd other than NKTOPG, NSTRA can be symbolic
constants. sugwd now only sets NKTOPG, NSTRA. Simplified the
computation of NKTOPG, NSTRA by making the local variable zpm1r an
array instead of a scalar and calling ifirstloc.

1 module gwprofil_m
2
3 IMPLICIT NONE
4
5 contains
6
7 SUBROUTINE gwprofil(nlon, nlev, ktest, kkcrith, kcrit, paphm1, prho, pstab, &
8 pvph, pri, ptau, pdmod, psig, pvar)
9
10 ! 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 ! linearly with height from the top of the blocked layer to 3 *
14 ! varor (kknu), to simulate lee waves or nonlinear gravity wave
15 ! breaking. Above it is constant, except when the wave encounters
16 ! a critical level (kcrit) or when it breaks.
17
18 ! Reference. See ECMWF research department documentation of the
19 ! "I.F.S."
20
21 ! Modifications. Passage of the new gwdrag TO I.F.S. (F. LOTT,
22 ! 22/11/93)
23
24 USE dimphy, ONLY : klev, klon
25 USE yoegwd, ONLY : gkdrag, grahilo, grcrit, gssec, gtsec, nstra
26
27 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 REAL, INTENT (IN) :: psig(nlon)
35 REAL, INTENT (IN) :: pvar(nlon)
36
37 ! Local:
38 INTEGER ilevh, jl, jk
39 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
44 !-----------------------------------------------------------------------
45
46 ! 1. INITIALIZATION
47
48 ! COMPUTATIONAL CONSTANTS.
49
50 ilevh = klev/3
51
52 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
59 DO jk = klev, 2, -1
60 ! 4.1 CONSTANT WAVE STRESS UNTIL TOP OF THE
61 ! BLOCKING LAYER.
62 DO jl = 1, klon
63 IF (ktest(jl)==1) THEN
64 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 END IF
70 end DO
71
72 ! 4.2 WAVE DISPLACEMENT AT NEXT LEVEL.
73 DO jl = 1, klon
74 IF (ktest(jl)==1) THEN
75 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 END IF
81 end DO
82
83 ! 4.3 WAVE RICHARDSON NUMBER, NEW WAVE DISPLACEMENT
84 ! AND STRESS: BREAKING EVALUATION AND CRITICAL
85 ! LEVEL
86 DO jl = 1, klon
87 IF (ktest(jl)==1) THEN
88 IF (jk<kkcrith(jl)) THEN
89 IF ((ptau(jl, jk+1)<gtsec) .OR. (jk<=kcrit(jl))) THEN
90 ptau(jl, jk) = 0.0
91 ELSE
92 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+zalfa*zsqr)**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 END IF
106 END IF
107 END IF
108 end DO
109 end DO
110
111 ! REORGANISATION OF THE STRESS PROFILE AT LOW LEVEL
112
113 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
120 DO jk = 1, klev
121 DO jl = 1, klon
122 IF (ktest(jl)==1) THEN
123 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 END IF
130 end DO
131
132 ! REORGANISATION IN THE STRATOSPHERE
133 DO jl = 1, klon
134 IF (ktest(jl)==1) THEN
135 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 END IF
141 end DO
142
143 ! REORGANISATION IN THE TROPOSPHERE
144 DO jl = 1, klon
145 IF (ktest(jl)==1) THEN
146 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 END IF
153 end DO
154 end DO
155
156 END SUBROUTINE gwprofil
157
158 end module gwprofil_m

  ViewVC Help
Powered by ViewVC 1.1.21