/[lmdze]/trunk/libf/phylmd/Conflx/flxasc.f90
ViewVC logotype

Contents of /trunk/libf/phylmd/Conflx/flxasc.f90

Parent Directory Parent Directory | Revision Log Revision Log


Revision 71 - (show annotations)
Mon Jul 8 18:12:18 2013 UTC (10 years, 9 months ago) by guez
File size: 10572 byte(s)
No reason to call inidissip in ce0l.

In inidissip, set random seed to 1 beacuse PGI compiler does not
accept all zeros.

dq was computed needlessly in caladvtrac. Arguments masse and dq of
calfis not used.

Replaced real*8 by double precision.

Pass arrays with inverted order of vertical levels to conflx instead
of creating local variables for this inside conflx.

1 module flxasc_m
2
3 IMPLICIT none
4
5 contains
6
7 SUBROUTINE flxasc(pdtime, ptenh, pqenh, pten, pqen, pqsen, pgeo, pgeoh, &
8 pap, paph, pqte, pvervel, ldland, ldcum, ktype, klab, ptu, pqu, plu, &
9 pmfu, pmfub, pentr, pmfus, pmfuq, pmful, plude, pdmfup, kcbot, kctop, &
10 kctop0, kcum, pen_u, pde_u)
11
12 ! This routine does the calculations for cloud ascents for cumulus
13 ! parameterization.
14
15 USE dimphy, ONLY: klev, klon
16 use flxadjtq_m, only: flxadjtq
17 USE suphec_m, ONLY: rcpd, rd, retv, rg, rtt
18 USE yoecumf, ONLY: cmfcmin, cmfctop, cprcon, entrmid, lmfmid
19
20 REAL, intent(in):: pdtime
21 REAL, intent(in):: ptenh(klon, klev)
22 REAL, intent(in):: pqenh(klon, klev)
23 REAL, intent(in):: pten(klon, klev)
24 REAL, intent(in):: pqen(klon, klev)
25 REAL, intent(in):: pqsen(klon, klev)
26 REAL, intent(in):: pgeo(klon, klev), pgeoh(klon, klev)
27 REAL pap(klon, klev), paph(klon, klev+1)
28 REAL pqte(klon, klev)
29 REAL pvervel(klon, klev) ! vitesse verticale en Pa/s
30 LOGICAL ldland(klon)
31 LOGICAL, intent(inout):: ldcum(klon)
32 INTEGER, intent(inout):: ktype(klon)
33 integer klab(klon, klev)
34 REAL ptu(klon, klev), pqu(klon, klev), plu(klon, klev)
35 REAL pmfu(klon, klev)
36 REAL, intent(inout):: pmfub(klon)
37 real pentr(klon)
38 real pmfus(klon, klev)
39 REAL pmfuq(klon, klev), pmful(klon, klev)
40 REAL plude(klon, klev)
41 REAL pdmfup(klon, klev)
42 integer kcbot(klon), kctop(klon)
43 INTEGER kctop0(klon)
44 integer, intent(out):: kcum
45 REAL pen_u(klon, klev), pde_u(klon, klev)
46
47 ! Local:
48
49 REAL zqold(klon)
50 REAL zdland(klon)
51 LOGICAL llflag(klon)
52 INTEGER k, i, is, icall
53 REAL ztglace, zdphi, zqeen, zseen, zscde, zqude
54 REAL zmfusk, zmfuqk, zmfulk, zbuo, zdnoprc, zprcon, zlnew
55
56 REAL zpbot(klon), zptop(klon), zrho(klon)
57 REAL zdprho, zentr, zpmid, zmftest, zmfmax
58 LOGICAL llo1, llo2
59
60 REAL zwmax(klon), zzzmb
61 INTEGER klwmin(klon) ! level of maximum vertical velocity
62
63 !----------------------------------------------------------------------
64
65 ztglace = RTT - 13.
66
67 ! Chercher le niveau où la vitesse verticale est maximale :
68
69 DO i = 1, klon
70 klwmin(i) = klev
71 zwmax(i) = 0.0
72 ENDDO
73
74 DO k = klev, 3, -1
75 DO i = 1, klon
76 IF (pvervel(i, k) < zwmax(i)) THEN
77 zwmax(i) = pvervel(i, k)
78 klwmin(i) = k
79 ENDIF
80 ENDDO
81 ENDDO
82
83 ! Set default values:
84
85 DO i = 1, klon
86 IF (.NOT. ldcum(i)) ktype(i)=0
87 ENDDO
88
89 DO k=1, klev
90 DO i = 1, klon
91 plu(i, k)=0.
92 pmfu(i, k)=0.
93 pmfus(i, k)=0.
94 pmfuq(i, k)=0.
95 pmful(i, k)=0.
96 plude(i, k)=0.
97 pdmfup(i, k)=0.
98 IF(.NOT. ldcum(i).OR.ktype(i) == 3) klab(i, k)=0
99 IF(.NOT. ldcum(i).AND.paph(i, k) < 4.E4) kctop0(i)=k
100 ENDDO
101 ENDDO
102
103 DO i = 1, klon
104 IF (ldland(i)) THEN
105 zdland(i)=3.0E4
106 zdphi=pgeoh(i, kctop0(i))-pgeoh(i, kcbot(i))
107 IF (ptu(i, kctop0(i)) >= ztglace) zdland(i)=zdphi
108 zdland(i)=MAX(3.0E4, zdland(i))
109 zdland(i)=MIN(5.0E4, zdland(i))
110 ENDIF
111 ENDDO
112
113 ! Initialiser les valeurs au niveau d'ascendance
114
115 DO i = 1, klon
116 kctop(i) = klev-1
117 IF (.NOT. ldcum(i)) THEN
118 kcbot(i) = klev-1
119 pmfub(i) = 0.
120 pqu(i, klev) = 0.
121 ENDIF
122 pmfu(i, klev) = pmfub(i)
123 pmfus(i, klev) = pmfub(i) * (RCPD * ptu(i, klev)+pgeoh(i, klev))
124 pmfuq(i, klev) = pmfub(i) * pqu(i, klev)
125 ENDDO
126
127 DO i = 1, klon
128 ldcum(i) = .FALSE.
129 ENDDO
130
131 ! Do ascent: subcloud layer (klab=1), clouds (klab=2) by doing
132 ! first dry-adiabatic ascent and then by adjusting t, q and l
133 ! accordingly in flxadjtq, then check for buoyancy and set flags
134 ! accordingly.
135
136 DO k = klev - 1, 3, -1
137 IF (LMFMID .AND. k < klev - 1 .AND. k > klev / 2) THEN
138 DO i = 1, klon
139 IF (.NOT. ldcum(i) .AND. klab(i, k + 1) == 0 .AND. &
140 pqen(i, k) > 0.9 * pqsen(i, k)) THEN
141 ptu(i, k+1) = pten(i, k) +(pgeo(i, k)-pgeoh(i, k+1))/RCPD
142 pqu(i, k+1) = pqen(i, k)
143 plu(i, k+1) = 0.0
144 zzzmb = MAX(CMFCMIN, -pvervel(i, k)/RG)
145 zmfmax = (paph(i, k)-paph(i, k-1))/(RG * pdtime)
146 pmfub(i) = MIN(zzzmb, zmfmax)
147 pmfu(i, k+1) = pmfub(i)
148 pmfus(i, k+1) = pmfub(i) * (RCPD * ptu(i, k+1)+pgeoh(i, k+1))
149 pmfuq(i, k+1) = pmfub(i) * pqu(i, k+1)
150 pmful(i, k+1) = 0.0
151 pdmfup(i, k+1) = 0.0
152 kcbot(i) = k
153 klab(i, k+1) = 1
154 ktype(i) = 3
155 pentr(i) = ENTRMID
156 ENDIF
157 ENDDO
158 ENDIF
159
160 is = 0
161 DO i = 1, klon
162 is = is + klab(i, k+1)
163 IF (klab(i, k+1) == 0) klab(i, k) = 0
164 llflag(i) = .FALSE.
165 IF (klab(i, k+1) > 0) llflag(i) = .TRUE.
166 ENDDO
167 IF (is == 0) cycle
168
169 ! Calculer le taux d'entraînement et de détraînement :
170
171 DO i = 1, klon
172 pen_u(i, k) = 0.0
173 pde_u(i, k) = 0.0
174 zrho(i)=paph(i, k+1)/(RD * ptenh(i, k+1))
175 zpbot(i)=paph(i, kcbot(i))
176 zptop(i)=paph(i, kctop0(i))
177 ENDDO
178
179 DO i = 1, klon
180 IF(ldcum(i)) THEN
181 zdprho=(paph(i, k+1)-paph(i, k))/(RG * zrho(i))
182 zentr=pentr(i) * pmfu(i, k+1) * zdprho
183 llo1=k < kcbot(i)
184 IF(llo1) pde_u(i, k)=zentr
185 zpmid=0.5 * (zpbot(i)+zptop(i))
186 llo2=llo1.AND.ktype(i) == 2.AND. &
187 (zpbot(i)-paph(i, k) < 0.2E5.OR. &
188 paph(i, k) > zpmid)
189 IF(llo2) pen_u(i, k)=zentr
190 llo2=llo1.AND.(ktype(i) == 1.OR.ktype(i) == 3).AND. &
191 (k >= MAX(klwmin(i), kctop0(i)+2).OR.pap(i, k) > zpmid)
192 IF(llo2) pen_u(i, k)=zentr
193 llo1=pen_u(i, k) > 0..AND.(ktype(i) == 1.OR.ktype(i) == 2)
194 IF(llo1) THEN
195 zentr=zentr * (1.+3. * (1.-MIN(1., (zpbot(i)-pap(i, k))/1.5E4)))
196 pen_u(i, k)=pen_u(i, k) * (1.+3. * (1.-MIN(1., &
197 (zpbot(i)-pap(i, k))/1.5E4)))
198 pde_u(i, k)=pde_u(i, k) * (1.+3. * (1.-MIN(1., &
199 (zpbot(i)-pap(i, k))/1.5E4)))
200 ENDIF
201 IF(llo2.AND.pqenh(i, k+1) > 1.E-5) &
202 pen_u(i, k)=zentr+MAX(pqte(i, k), 0.)/pqenh(i, k+1) * &
203 zrho(i) * zdprho
204 ENDIF
205 end DO
206
207 ! Do adiabatic ascent for entraining/detraining plume
208
209 DO i = 1, klon
210 IF (llflag(i)) THEN
211 IF (k < kcbot(i)) THEN
212 zmftest = pmfu(i, k+1)+pen_u(i, k)-pde_u(i, k)
213 zmfmax = MIN(zmftest, (paph(i, k)-paph(i, k-1))/(RG * pdtime))
214 pen_u(i, k)=MAX(pen_u(i, k)-MAX(0.0, zmftest-zmfmax), 0.0)
215 ENDIF
216 pde_u(i, k)=MIN(pde_u(i, k), 0.75 * pmfu(i, k+1))
217 ! calculer le flux de masse du niveau k a partir de celui du k+1
218 pmfu(i, k)=pmfu(i, k+1)+pen_u(i, k)-pde_u(i, k)
219 ! calculer les valeurs Su, Qu et l du niveau k dans le
220 ! panache montant
221 zqeen=pqenh(i, k+1) * pen_u(i, k)
222 zseen=(RCPD * ptenh(i, k+1)+pgeoh(i, k+1)) * pen_u(i, k)
223 zscde=(RCPD * ptu(i, k+1)+pgeoh(i, k+1)) * pde_u(i, k)
224 zqude=pqu(i, k+1) * pde_u(i, k)
225 plude(i, k)=plu(i, k+1) * pde_u(i, k)
226 zmfusk=pmfus(i, k+1)+zseen-zscde
227 zmfuqk=pmfuq(i, k+1)+zqeen-zqude
228 zmfulk=pmful(i, k+1) -plude(i, k)
229 plu(i, k)=zmfulk * (1./MAX(CMFCMIN, pmfu(i, k)))
230 pqu(i, k)=zmfuqk * (1./MAX(CMFCMIN, pmfu(i, k)))
231 ptu(i, k)=(zmfusk * (1./MAX(CMFCMIN, pmfu(i, k)))- &
232 pgeoh(i, k))/RCPD
233 ptu(i, k)=MAX(100., ptu(i, k))
234 ptu(i, k)=MIN(400., ptu(i, k))
235 zqold(i)=pqu(i, k)
236 ELSE
237 zqold(i)=0.0
238 ENDIF
239 end DO
240
241 ! Do corrections for moist ascent by adjusting t, q and l
242
243 icall = 1
244 CALL flxadjtq(paph(1, k), ptu(1, k), pqu(1, k), llflag, icall)
245
246 DO i = 1, klon
247 IF(llflag(i).AND.pqu(i, k).NE.zqold(i)) THEN
248 klab(i, k) = 2
249 plu(i, k) = plu(i, k)+zqold(i)-pqu(i, k)
250 zbuo = ptu(i, k) * (1.+RETV * pqu(i, k))- &
251 ptenh(i, k) * (1.+RETV * pqenh(i, k))
252 IF (klab(i, k+1) == 1) zbuo=zbuo+0.5
253 IF (zbuo > 0. .AND. pmfu(i, k) >= 0.1 * pmfub(i)) THEN
254 kctop(i) = k
255 ldcum(i) = .TRUE.
256 zdnoprc = 1.5E4
257 IF (ldland(i)) zdnoprc = zdland(i)
258 zprcon = CPRCON
259 IF ((zpbot(i)-paph(i, k)) < zdnoprc) zprcon = 0.0
260 zlnew=plu(i, k)/(1.+zprcon * (pgeoh(i, k)-pgeoh(i, k+1)))
261 pdmfup(i, k)=MAX(0., (plu(i, k)-zlnew) * pmfu(i, k))
262 plu(i, k)=zlnew
263 ELSE
264 klab(i, k)=0
265 pmfu(i, k)=0.
266 ENDIF
267 ENDIF
268 end DO
269 DO i = 1, klon
270 IF (llflag(i)) THEN
271 pmful(i, k)=plu(i, k) * pmfu(i, k)
272 pmfus(i, k)=(RCPD * ptu(i, k)+pgeoh(i, k)) * pmfu(i, k)
273 pmfuq(i, k)=pqu(i, k) * pmfu(i, k)
274 ENDIF
275 end DO
276
277 end DO
278
279 ! Determine convective fluxes above non-buoyancy level (note:
280 ! cloud variables like t, q and l are not affected by detrainment
281 ! and are already known from previous calculations above).
282
283 DO i = 1, klon
284 IF (kctop(i) == klev-1) ldcum(i) = .FALSE.
285 kcbot(i) = MAX(kcbot(i), kctop(i))
286 ENDDO
287
288 ldcum(1)=ldcum(1)
289
290 is = 0
291 DO i = 1, klon
292 if (ldcum(i)) is = is + 1
293 ENDDO
294 kcum = is
295 IF (is /= 0) then
296 DO i = 1, klon
297 IF (ldcum(i)) THEN
298 k=kctop(i)-1
299 pde_u(i, k)=(1.-CMFCTOP) * pmfu(i, k+1)
300 plude(i, k)=pde_u(i, k) * plu(i, k+1)
301 pmfu(i, k)=pmfu(i, k+1)-pde_u(i, k)
302 zlnew=plu(i, k)
303 pdmfup(i, k)=MAX(0., (plu(i, k)-zlnew) * pmfu(i, k))
304 plu(i, k)=zlnew
305 pmfus(i, k)=(RCPD * ptu(i, k)+pgeoh(i, k)) * pmfu(i, k)
306 pmfuq(i, k)=pqu(i, k) * pmfu(i, k)
307 pmful(i, k)=plu(i, k) * pmfu(i, k)
308 plude(i, k-1)=pmful(i, k)
309 ENDIF
310 end DO
311 end IF
312
313 END SUBROUTINE flxasc
314
315 end module flxasc_m

  ViewVC Help
Powered by ViewVC 1.1.21