1 | ! $Header$ |
---|
2 | |
---|
3 | SUBROUTINE convect3(dtime, epmax, ok_adj, t1, r1, rs, u, v, tra, p, ph, nd, & |
---|
4 | ndp1, nl, ntra, delt, iflag, ft, fr, fu, fv, ftra, precip, icb, inb, & |
---|
5 | upwd, dnwd, dnwd0, sig, w0, mike, mke, ma, ments, qents, tps, tls, sigij, & |
---|
6 | cape, tvp, pbase, buoybase, & ! ccc * |
---|
7 | ! DTVPDT1,DTVPDQ1,DPLCLDT,DPLCLDR) |
---|
8 | dtvpdt1, dtvpdq1, dplcldt, dplcldr, & ! sbl |
---|
9 | ft2, fr2, fu2, fv2, wd, qcond, qcondc) ! sbl |
---|
10 | |
---|
11 | ! *** THE PARAMETER NA SHOULD IN GENERAL EQUAL ND *** |
---|
12 | |
---|
13 | ! ################################################################# |
---|
14 | ! Fleur Introduction des traceurs dans convect3 le 6 juin 200 |
---|
15 | ! ################################################################# |
---|
16 | USE dimphy |
---|
17 | USE infotrac_phy, ONLY: nbtr |
---|
18 | USE lmdz_yomcst |
---|
19 | |
---|
20 | IMPLICIT NONE |
---|
21 | INTEGER na |
---|
22 | PARAMETER (na = 60) |
---|
23 | |
---|
24 | REAL deltac ! cld |
---|
25 | PARAMETER (deltac = 0.01) ! cld |
---|
26 | |
---|
27 | INTEGER nent(na) |
---|
28 | INTEGER nd, ndp1, nl, ntra, iflag, icb, inb |
---|
29 | REAL dtime, epmax, delt, precip, cape |
---|
30 | REAL dplcldt, dplcldr |
---|
31 | REAL t1(nd), r1(nd), rs(nd), u(nd), v(nd), tra(nd, ntra) |
---|
32 | REAL p(nd), ph(ndp1) |
---|
33 | REAL ft(nd), fr(nd), fu(nd), fv(nd), ftra(nd, ntra) |
---|
34 | REAL sig(nd), w0(nd) |
---|
35 | REAL uent(na, na), vent(na, na), traent(na, na, nbtr), tratm(na) |
---|
36 | REAL up(na), vp(na), trap(na, nbtr) |
---|
37 | REAL m(na), mp(na), ment(na, na), qent(na, na), elij(na, na) |
---|
38 | REAL sij(na, na), tvp(na), tv(na), water(na) |
---|
39 | REAL rp(na), ep(na), th(na), wt(na), evap(na), clw(na) |
---|
40 | REAL sigp(na), b(na), tp(na), cpn(na) |
---|
41 | REAL lv(na), lvcp(na), h(na), hp(na), gz(na) |
---|
42 | REAL t(na), rr(na) |
---|
43 | |
---|
44 | REAL ft2(nd), fr2(nd), fu2(nd), fv2(nd) ! added sbl |
---|
45 | REAL u1(nd), v1(nd) ! added sbl |
---|
46 | |
---|
47 | REAL buoy(na) ! Lifted parcel buoyancy |
---|
48 | REAL dtvpdt1(nd), dtvpdq1(nd) ! Derivatives of parcel virtual |
---|
49 | ! temperature wrt T1 and Q1 |
---|
50 | REAL clw_new(na), qi(na) |
---|
51 | |
---|
52 | REAL wd, betad ! for gust factor (sb) |
---|
53 | REAL qcondc(nd) ! interface cld param (sb) |
---|
54 | REAL qcond(nd), nqcond(na), wa(na), maa(na), siga(na), axc(na) ! cld |
---|
55 | |
---|
56 | LOGICAL ice_conv, ok_adj |
---|
57 | PARAMETER (ice_conv = .TRUE.) |
---|
58 | |
---|
59 | ! ccccccccccccccccccccccccccccccccccccccccccccc |
---|
60 | ! declaration des variables a sortir |
---|
61 | ! cccccccccccccccccccccccccccccccccccccccccccc |
---|
62 | REAL mke(nd) |
---|
63 | REAL mike(nd) |
---|
64 | REAL ma(nd) |
---|
65 | REAL tps(nd) !temperature dans les ascendances non diluees |
---|
66 | REAL tls(nd) !temperature potentielle |
---|
67 | REAL ments(nd, nd) |
---|
68 | REAL qents(nd, nd) |
---|
69 | REAL sigij(klev, klev) |
---|
70 | REAL pbase ! pressure at the cloud base level |
---|
71 | REAL buoybase ! buoyancy at the cloud base level |
---|
72 | ! ccccccccccccccccccccccccccccccccccccccccccccc |
---|
73 | |
---|
74 | REAL :: cpv, cl, cpvmcl, eps, alv0, rdcp, pbcrit, ptcrit, sigd, spfac |
---|
75 | REAL :: tau, beta, alpha, dtcrit, dtovsh, ahm, rm, um, vm, dphinv |
---|
76 | REAL :: a2, x, tvx, tvy, plcl, pden, dpbase, tvpbase, tvbase, tdif |
---|
77 | REAL :: ath1, ath, delti, deltap, dcape, dlnp, sigold, dtmin, fac, w |
---|
78 | REAL :: amu, rti, cpd, bf2, anum, denom, dei, altem, cwat, stemp, qp |
---|
79 | REAL :: scrit, alt, smax, asij, wgh, sjmax, sjmin, smid, delp, delm |
---|
80 | REAL :: asum, bsum, csum, wflux, tinv, wdtrain, awat, afac, afac1, afac2 |
---|
81 | REAL :: bfac, pr1, pr2, sigt, b6, c6, revap, tevap, delth, amfac, amp2 |
---|
82 | REAL :: xf, tf, af, bf, fac2, ur, sru, d, ampmax, dpinv, am, amde, cpinv |
---|
83 | REAL :: amp1, ad, rat, ax, bx, cx, dx, ex, dsum |
---|
84 | INTEGER :: nk, i, j, nopt, jn, k, im, jm, n |
---|
85 | |
---|
86 | REAL dnwd0(nd) ! precipitation driven unsaturated downdraft flux |
---|
87 | REAL dnwd(nd), dn1 ! in-cloud saturated downdraft mass flux |
---|
88 | REAL upwd(nd), up1 ! in-cloud saturated updraft mass flux |
---|
89 | |
---|
90 | ! *** ASSIGN VALUES OF THERMODYNAMIC CONSTANTS *** |
---|
91 | ! *** THESE SHOULD BE CONSISTENT WITH *** |
---|
92 | ! *** THOSE USED IN CALLING PROGRAM *** |
---|
93 | ! *** NOTE: THESE ARE ALSO SPECIFIED IN SUBROUTINE TLIFT *** |
---|
94 | |
---|
95 | ! sb CPD=1005.7 |
---|
96 | ! sb CPV=1870.0 |
---|
97 | ! sb CL=4190.0 |
---|
98 | ! sb CPVMCL=CL-CPV |
---|
99 | ! sb RV=461.5 |
---|
100 | ! sb RD=287.04 |
---|
101 | ! sb EPS=RD/RV |
---|
102 | ! sb ALV0=2.501E6 |
---|
103 | ! cccccccccccccccccccccc |
---|
104 | ! constantes coherentes avec le modele du Centre Europeen |
---|
105 | ! sb RD = 1000.0 * 1.380658E-23 * 6.0221367E+23 / 28.9644 |
---|
106 | ! sb RV = 1000.0 * 1.380658E-23 * 6.0221367E+23 / 18.0153 |
---|
107 | ! sb CPD = 3.5 * RD |
---|
108 | ! sb CPV = 4.0 * RV |
---|
109 | ! sb CL = 4218.0 |
---|
110 | ! sb CPVMCL=CL-CPV |
---|
111 | ! sb EPS=RD/RV |
---|
112 | ! sb ALV0=2.5008E+06 |
---|
113 | ! ccccccccccccccccccccc |
---|
114 | ! on utilise les constantes thermo du Centre Europeen: (SB) |
---|
115 | |
---|
116 | cpd = rcpd |
---|
117 | cpv = rcpv |
---|
118 | cl = rcw |
---|
119 | cpvmcl = cl - cpv |
---|
120 | eps = rd / rv |
---|
121 | alv0 = rlvtt |
---|
122 | |
---|
123 | nk = 1 ! origin level of the lifted parcel |
---|
124 | |
---|
125 | ! ccccccccccccccccccccc |
---|
126 | |
---|
127 | ! *** INITIALIZE OUTPUT ARRAYS AND PARAMETERS *** |
---|
128 | |
---|
129 | DO i = 1, nd |
---|
130 | ft(i) = 0.0 |
---|
131 | fr(i) = 0.0 |
---|
132 | fu(i) = 0.0 |
---|
133 | fv(i) = 0.0 |
---|
134 | |
---|
135 | ft2(i) = 0.0 |
---|
136 | fr2(i) = 0.0 |
---|
137 | fu2(i) = 0.0 |
---|
138 | fv2(i) = 0.0 |
---|
139 | |
---|
140 | DO j = 1, ntra |
---|
141 | ftra(i, j) = 0.0 |
---|
142 | END DO |
---|
143 | |
---|
144 | qcondc(i) = 0.0 ! cld |
---|
145 | qcond(i) = 0.0 ! cld |
---|
146 | nqcond(i) = 0.0 ! cld |
---|
147 | |
---|
148 | t(i) = t1(i) |
---|
149 | rr(i) = r1(i) |
---|
150 | u1(i) = u(i) ! added sbl |
---|
151 | v1(i) = v(i) ! added sbl |
---|
152 | END DO |
---|
153 | DO i = 1, nl |
---|
154 | rdcp = (rd * (1. - rr(i)) + rr(i) * rv) / (cpd * (1. - rr(i)) + rr(i) * cpv) |
---|
155 | th(i) = t(i) * (1000.0 / p(i))**rdcp |
---|
156 | END DO |
---|
157 | |
---|
158 | ! ************************************************************ |
---|
159 | ! * CALCUL DES TEMPERATURES POTENTIELLES A SORTIR |
---|
160 | ! ************************************************************ |
---|
161 | DO i = 1, nd |
---|
162 | rdcp = (rd * (1. - rr(i)) + rr(i) * rv) / (cpd * (1. - rr(i)) + rr(i) * cpv) |
---|
163 | |
---|
164 | tls(i) = t(i) * (1000.0 / p(i))**rdcp |
---|
165 | END DO |
---|
166 | |
---|
167 | |
---|
168 | |
---|
169 | |
---|
170 | ! *********************************************************** |
---|
171 | |
---|
172 | precip = 0.0 |
---|
173 | wd = 0.0 ! sb |
---|
174 | iflag = 1 |
---|
175 | |
---|
176 | ! *** SPECIFY PARAMETERS *** |
---|
177 | ! *** PBCRIT IS THE CRITICAL CLOUD DEPTH (MB) BENEATH WHICH THE *** |
---|
178 | ! *** PRECIPITATION EFFICIENCY IS ASSUMED TO BE ZERO *** |
---|
179 | ! *** PTCRIT IS THE CLOUD DEPTH (MB) ABOVE WHICH THE PRECIP. *** |
---|
180 | ! *** EFFICIENCY IS ASSUMED TO BE UNITY *** |
---|
181 | ! *** SIGD IS THE FRACTIONAL AREA COVERED BY UNSATURATED DNDRAFT *** |
---|
182 | ! *** SPFAC IS THE FRACTION OF PRECIPITATION FALLING OUTSIDE *** |
---|
183 | ! *** OF CLOUD *** |
---|
184 | ! *** ALPHA AND BETA ARE PARAMETERS THAT CONTROL THE RATE OF *** |
---|
185 | ! *** APPROACH TO QUASI-EQUILIBRIUM *** |
---|
186 | ! *** (THEIR STANDARD VALUES ARE 1.0 AND 0.96, RESPECTIVELY) *** |
---|
187 | ! *** (BETA MUST BE LESS THAN OR EQUAL TO 1) *** |
---|
188 | ! *** DTCRIT IS THE CRITICAL BUOYANCY (K) USED TO ADJUST THE *** |
---|
189 | ! *** APPROACH TO QUASI-EQUILIBRIUM *** |
---|
190 | ! *** IT MUST BE LESS THAN 0 *** |
---|
191 | |
---|
192 | pbcrit = 150.0 |
---|
193 | ptcrit = 500.0 |
---|
194 | sigd = 0.01 |
---|
195 | spfac = 0.15 |
---|
196 | ! sb: |
---|
197 | ! EPMAX=0.993 ! precip efficiency less than unity |
---|
198 | ! EPMAX=1. ! precip efficiency less than unity |
---|
199 | |
---|
200 | ! jyg |
---|
201 | ! CC BETA=0.96 |
---|
202 | ! Beta is now expressed as a function of the characteristic time |
---|
203 | ! of the convective process. |
---|
204 | ! CC Old value : TAU = 15000. !(for dtime = 600.s) |
---|
205 | ! CC Other value (inducing little change) :TAU = 8000. |
---|
206 | tau = 8000. |
---|
207 | beta = 1. - dtime / tau |
---|
208 | ! jyg |
---|
209 | ! CC ALPHA=1.0 |
---|
210 | alpha = 1.5E-3 * dtime / tau |
---|
211 | ! Increase alpha in order to compensate W decrease |
---|
212 | alpha = alpha * 1.5 |
---|
213 | |
---|
214 | ! jyg (voir CONVECT 3) |
---|
215 | ! CC DTCRIT=-0.2 |
---|
216 | dtcrit = -2. |
---|
217 | ! gf&jyg |
---|
218 | ! CC DT pour l'overshoot. |
---|
219 | dtovsh = -0.2 |
---|
220 | |
---|
221 | |
---|
222 | ! *** INCREMENT THE COUNTER *** |
---|
223 | |
---|
224 | sig(nd) = sig(nd) + 1.0 |
---|
225 | sig(nd) = amin1(sig(nd), 12.1) |
---|
226 | |
---|
227 | ! *** IF NOPT IS AN INTEGER OTHER THAN 0, CONVECT *** |
---|
228 | ! *** RETURNS ARRAYS T AND R THAT MAY HAVE BEEN *** |
---|
229 | ! *** ALTERED BY DRY ADIABATIC ADJUSTMENT; OTHERWISE *** |
---|
230 | ! *** THE RETURNED ARRAYS ARE UNALTERED. *** |
---|
231 | |
---|
232 | nopt = 0 |
---|
233 | ! NOPT=1 ! sbl |
---|
234 | |
---|
235 | ! *** PERFORM DRY ADIABATIC ADJUSTMENT *** |
---|
236 | |
---|
237 | ! *** DO NOT BYPASS THIS EVEN IF THE CALLING PROGRAM HAS A *** |
---|
238 | ! *** BOUNDARY LAYER SCHEME !!! *** |
---|
239 | |
---|
240 | IF (ok_adj) THEN ! added sbl |
---|
241 | |
---|
242 | DO i = nl - 1, 1, -1 |
---|
243 | jn = 0 |
---|
244 | DO j = i + 1, nl |
---|
245 | IF (th(j)<th(i)) jn = j |
---|
246 | END DO |
---|
247 | IF (jn==0) GO TO 30 |
---|
248 | ahm = 0.0 |
---|
249 | rm = 0.0 |
---|
250 | um = 0.0 |
---|
251 | vm = 0.0 |
---|
252 | DO k = 1, ntra |
---|
253 | tratm(k) = 0.0 |
---|
254 | END DO |
---|
255 | DO j = i, jn |
---|
256 | ahm = ahm + (cpd * (1. - rr(j)) + rr(j) * cpv) * t(j) * (ph(j) - ph(j + 1)) |
---|
257 | rm = rm + rr(j) * (ph(j) - ph(j + 1)) |
---|
258 | um = um + u(j) * (ph(j) - ph(j + 1)) |
---|
259 | vm = vm + v(j) * (ph(j) - ph(j + 1)) |
---|
260 | DO k = 1, ntra |
---|
261 | tratm(k) = tratm(k) + tra(j, k) * (ph(j) - ph(j + 1)) |
---|
262 | END DO |
---|
263 | END DO |
---|
264 | dphinv = 1. / (ph(i) - ph(jn + 1)) |
---|
265 | rm = rm * dphinv |
---|
266 | um = um * dphinv |
---|
267 | vm = vm * dphinv |
---|
268 | DO k = 1, ntra |
---|
269 | tratm(k) = tratm(k) * dphinv |
---|
270 | END DO |
---|
271 | a2 = 0.0 |
---|
272 | DO j = i, jn |
---|
273 | rr(j) = rm |
---|
274 | u(j) = um |
---|
275 | v(j) = vm |
---|
276 | DO k = 1, ntra |
---|
277 | tra(j, k) = tratm(k) |
---|
278 | END DO |
---|
279 | rdcp = (rd * (1. - rr(j)) + rr(j) * rv) / (cpd * (1. - rr(j)) + rr(j) * cpv) |
---|
280 | x = (0.001 * p(j))**rdcp |
---|
281 | t(j) = x |
---|
282 | a2 = a2 + (cpd * (1. - rr(j)) + rr(j) * cpv) * x * (ph(j) - ph(j + 1)) |
---|
283 | END DO |
---|
284 | DO j = i, jn |
---|
285 | th(j) = ahm / a2 |
---|
286 | t(j) = t(j) * th(j) |
---|
287 | END DO |
---|
288 | 30 END DO |
---|
289 | |
---|
290 | END IF ! added sbl |
---|
291 | |
---|
292 | ! *** RESET INPUT ARRAYS IF ok_adj 0 *** |
---|
293 | |
---|
294 | IF (ok_adj) THEN |
---|
295 | DO i = 1, nd |
---|
296 | |
---|
297 | ft2(i) = (t(i) - t1(i)) / delt ! sbl |
---|
298 | fr2(i) = (rr(i) - r1(i)) / delt ! sbl |
---|
299 | fu2(i) = (u(i) - u1(i)) / delt ! sbl |
---|
300 | fv2(i) = (v(i) - v1(i)) / delt ! sbl |
---|
301 | |
---|
302 | ! T1(I)=T(I) ! commente sbl |
---|
303 | ! R1(I)=RR(I) ! commente sbl |
---|
304 | END DO |
---|
305 | END IF |
---|
306 | |
---|
307 | ! *** CALCULATE ARRAYS OF GEOPOTENTIAL, HEAT CAPACITY AND STATIC ENERGY |
---|
308 | |
---|
309 | gz(1) = 0.0 |
---|
310 | cpn(1) = cpd * (1. - rr(1)) + rr(1) * cpv |
---|
311 | h(1) = t(1) * cpn(1) |
---|
312 | DO i = 2, nl |
---|
313 | tvx = t(i) * (1. + rr(i) / eps - rr(i)) |
---|
314 | tvy = t(i - 1) * (1. + rr(i - 1) / eps - rr(i - 1)) |
---|
315 | gz(i) = gz(i - 1) + 0.5 * rd * (tvx + tvy) * (p(i - 1) - p(i)) / ph(i) |
---|
316 | cpn(i) = cpd * (1. - rr(i)) + cpv * rr(i) |
---|
317 | h(i) = t(i) * cpn(i) + gz(i) |
---|
318 | END DO |
---|
319 | |
---|
320 | ! *** CALCULATE LIFTED CONDENSATION LEVEL OF AIR AT LOWEST MODEL LEVEL *** |
---|
321 | ! *** (WITHIN 0.2% OF FORMULA OF BOLTON, MON. WEA. REV.,1980) *** |
---|
322 | |
---|
323 | IF (t(1)<250.0 .OR. rr(1)<=0.0) THEN |
---|
324 | iflag = 0 |
---|
325 | ! sb3d PRINT*,'je suis passe par 366' |
---|
326 | RETURN |
---|
327 | END IF |
---|
328 | |
---|
329 | ! jyg1 Utilisation de la SUBROUTINE CLIFT |
---|
330 | ! C RH=RR(1)/RS(1) |
---|
331 | ! C CHI=T(1)/(1669.0-122.0*RH-T(1)) |
---|
332 | ! C PLCL=P(1)*(RH**CHI) |
---|
333 | CALL clift(p(1), t(1), rr(1), rs(1), plcl, dplcldt, dplcldr) |
---|
334 | ! jyg2 |
---|
335 | ! sb3d PRINT *,' em_plcl,p1,t1,r1,rs1,rh ' |
---|
336 | ! sb3d $ ,PLCL,P(1),T(1),RR(1),RS(1),RH |
---|
337 | |
---|
338 | IF (plcl<200.0 .OR. plcl>=2000.0) THEN |
---|
339 | iflag = 2 |
---|
340 | RETURN |
---|
341 | END IF |
---|
342 | ! jyg1 |
---|
343 | ! Essais de modification de ICB |
---|
344 | |
---|
345 | ! *** CALCULATE FIRST LEVEL ABOVE LCL (=ICB) *** |
---|
346 | |
---|
347 | ! C ICB=NL-1 |
---|
348 | ! C DO 50 I=2,NL-1 |
---|
349 | ! C IF(P(I).LT.PLCL)THEN |
---|
350 | ! C ICB=MIN(ICB,I) ! ICB sup ou egal a 2 |
---|
351 | ! C END IF |
---|
352 | ! C 50 CONTINUE |
---|
353 | ! C IF(ICB.EQ.(NL-1))THEN |
---|
354 | ! C IFLAG=3 |
---|
355 | ! C RETURN |
---|
356 | ! C END IF |
---|
357 | |
---|
358 | ! *** CALCULATE LAYER CONTAINING LCL (=ICB) *** |
---|
359 | |
---|
360 | icb = nl - 1 |
---|
361 | ! sb DO 50 I=2,NL-1 |
---|
362 | DO i = 3, nl - 1 ! modif sb pour que ICB soit sup/egal a 2 |
---|
363 | ! la modification consiste a comparer PLCL a PH et non a P: |
---|
364 | ! ICB est defini par : PH(ICB)<PLCL<PH(ICB-!) |
---|
365 | IF (ph(i)<plcl) THEN |
---|
366 | icb = min(icb, i) |
---|
367 | END IF |
---|
368 | END DO |
---|
369 | IF (icb==(nl - 1)) THEN |
---|
370 | iflag = 3 |
---|
371 | RETURN |
---|
372 | END IF |
---|
373 | icb = icb - 1 ! ICB sup ou egal a 2 |
---|
374 | ! jyg2 |
---|
375 | |
---|
376 | |
---|
377 | |
---|
378 | ! *** SUBROUTINE TLIFT CALCULATES PART OF THE LIFTED PARCEL VIRTUAL |
---|
379 | ! *** |
---|
380 | ! *** TEMPERATURE, THE ACTUAL TEMPERATURE AND THE ADIABATIC |
---|
381 | ! *** |
---|
382 | ! *** LIQUID WATER CONTENT |
---|
383 | ! *** |
---|
384 | |
---|
385 | |
---|
386 | ! jyg1 |
---|
387 | ! make sure that "Cloud base" seen by TLIFT is actually the |
---|
388 | ! fisrt level where adiabatic ascent is saturated |
---|
389 | IF (plcl>p(icb)) THEN |
---|
390 | ! sb CALL TLIFT(P,T,RR,RS,GZ,PLCL,ICB,TVP,TP,CLW,ND,NL) |
---|
391 | CALL tlift(p, t, rr, rs, gz, plcl, icb, nk, tvp, tp, clw, nd, nl, & |
---|
392 | dtvpdt1, dtvpdq1) |
---|
393 | ELSE |
---|
394 | ! sb CALL TLIFT(P,T,RR,RS,GZ,PLCL,ICB+1,TVP,TP,CLW,ND,NL) |
---|
395 | CALL tlift(p, t, rr, rs, gz, plcl, icb + 1, nk, tvp, tp, clw, nd, nl, & |
---|
396 | dtvpdt1, dtvpdq1) |
---|
397 | END IF |
---|
398 | ! jyg2 |
---|
399 | |
---|
400 | ! ***************************************************************************** |
---|
401 | ! *** SORTIE DE LA TEMPERATURE DE L ASCENDANCE NON DILUE |
---|
402 | ! ***************************************************************************** |
---|
403 | DO i = 1, nd |
---|
404 | tps(i) = tp(i) |
---|
405 | END DO |
---|
406 | |
---|
407 | |
---|
408 | ! ***************************************************************************** |
---|
409 | |
---|
410 | |
---|
411 | ! *** SET THE PRECIPITATION EFFICIENCIES AND THE FRACTION OF *** |
---|
412 | ! *** PRECIPITATION FALLING OUTSIDE OF CLOUD *** |
---|
413 | ! *** THESE MAY BE FUNCTIONS OF TP(I), P(I) AND CLW(I) *** |
---|
414 | |
---|
415 | DO i = 1, nl |
---|
416 | pden = ptcrit - pbcrit |
---|
417 | |
---|
418 | ! jyg |
---|
419 | ! cc EP(I)=(P(ICB)-P(I)-PBCRIT)/PDEN |
---|
420 | ! sb EP(I)=(PLCL-P(I)-PBCRIT)/PDEN |
---|
421 | ep(i) = (plcl - p(i) - pbcrit) / pden * epmax ! sb |
---|
422 | |
---|
423 | ep(i) = amax1(ep(i), 0.0) |
---|
424 | ! sb EP(I)=AMIN1(EP(I),1.0) |
---|
425 | ep(i) = amin1(ep(i), epmax) ! sb |
---|
426 | sigp(i) = spfac |
---|
427 | |
---|
428 | ! *** CALCULATE VIRTUAL TEMPERATURE AND LIFTED PARCEL *** |
---|
429 | ! *** VIRTUAL TEMPERATURE *** |
---|
430 | |
---|
431 | tv(i) = t(i) * (1. + rr(i) / eps - rr(i)) |
---|
432 | ! cd1 |
---|
433 | ! . Keep all liquid water in lifted parcel (-> adiabatic CAPE) |
---|
434 | |
---|
435 | ! cc TVP(I)=TVP(I)-TP(I)*(RR(1)-EP(I)*CLW(I)) |
---|
436 | ! !!! sb TVP(I)=TVP(I)-TP(I)*RR(1) ! calcule dans tlift |
---|
437 | ! cd2 |
---|
438 | |
---|
439 | ! *** Calculate first estimate of buoyancy |
---|
440 | |
---|
441 | buoy(i) = tvp(i) - tv(i) |
---|
442 | END DO |
---|
443 | |
---|
444 | ! *** Set Cloud Base Buoyancy at (Plcl+DPbase) level buoyancy |
---|
445 | |
---|
446 | dpbase = -40. !That is 400m above LCL |
---|
447 | pbase = plcl + dpbase |
---|
448 | tvpbase = tvp(icb) * (pbase - p(icb + 1)) / (p(icb) - p(icb + 1)) + & |
---|
449 | tvp(icb + 1) * (p(icb) - pbase) / (p(icb) - p(icb + 1)) |
---|
450 | tvbase = tv(icb) * (pbase - p(icb + 1)) / (p(icb) - p(icb + 1)) + & |
---|
451 | tv(icb + 1) * (p(icb) - pbase) / (p(icb) - p(icb + 1)) |
---|
452 | |
---|
453 | ! test sb: |
---|
454 | ! @ WRITE(*,*) '++++++++++++++++++++++++++++++++++++++++' |
---|
455 | ! @ WRITE(*,*) 'plcl,dpbas,tvpbas,tvbas,tvp(icb),tvp(icb1) |
---|
456 | ! @ : ,tv(icb),tv(icb1)' |
---|
457 | ! @ WRITE(*,*) plcl,dpbase,tvpbase,tvbase,tvp(icb) |
---|
458 | ! @ L ,tvp(icb+1),tv(icb),tv(icb+1) |
---|
459 | ! @ WRITE(*,*) '++++++++++++++++++++++++++++++++++++++++' |
---|
460 | ! fin test sb |
---|
461 | buoybase = tvpbase - tvbase |
---|
462 | |
---|
463 | ! C Set buoyancy = BUOYBASE for all levels below BASE. |
---|
464 | ! C For safety, set : BUOY(ICB) = BUOYBASE |
---|
465 | DO i = icb, nl |
---|
466 | IF (p(i)>=pbase) THEN |
---|
467 | buoy(i) = buoybase |
---|
468 | END IF |
---|
469 | END DO |
---|
470 | buoy(icb) = buoybase |
---|
471 | |
---|
472 | ! sb3d PRINT *,'buoybase,tvp_tv,tvpbase,tvbase,pbase,plcl' |
---|
473 | ! sb3d $, buoybase,tvp(icb)-tv(icb),tvpbase,tvbase,pbase,plcl |
---|
474 | ! sb3d PRINT *,'TVP ',(tvp(i),i=1,nl) |
---|
475 | ! sb3d PRINT *,'TV ',(tv(i),i=1,nl) |
---|
476 | ! sb3d PRINT *, 'P ',(p(i),i=1,nl) |
---|
477 | ! sb3d PRINT *,'ICB ',icb |
---|
478 | ! test sb: |
---|
479 | ! @ WRITE(*,*) '@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@' |
---|
480 | ! @ WRITE(*,*) 'icb,icbs,inb,buoybase:' |
---|
481 | ! @ WRITE(*,*) icb,icb+1,inb,buoybase |
---|
482 | ! @ WRITE(*,*) 'k,tvp,tv,tp,buoy,ep: ' |
---|
483 | ! @ do k=1,nl |
---|
484 | ! @ WRITE(*,*) k,tvp(k),tv(k),tp(k),buoy(k),ep(k) |
---|
485 | ! @ enddo |
---|
486 | ! @ WRITE(*,*) '@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@' |
---|
487 | ! fin test sb |
---|
488 | |
---|
489 | |
---|
490 | |
---|
491 | ! *** MAKE SURE THAT COLUMN IS DRY ADIABATIC BETWEEN THE SURFACE *** |
---|
492 | ! *** AND CLOUD BASE, AND THAT LIFTED AIR IS POSITIVELY BUOYANT *** |
---|
493 | ! *** AT CLOUD BASE *** |
---|
494 | ! *** IF NOT, RETURN TO CALLING PROGRAM AFTER RESETTING *** |
---|
495 | ! *** SIG(I) AND W0(I) *** |
---|
496 | |
---|
497 | ! jyg |
---|
498 | ! CC TDIF=TVP(ICB)-TV(ICB) |
---|
499 | tdif = buoy(icb) |
---|
500 | ath1 = th(1) |
---|
501 | ! jyg |
---|
502 | ! CC ATH=TH(ICB-1)-1.0 |
---|
503 | ath = th(icb - 1) - 5.0 |
---|
504 | ! ATH=0. ! ajout sb |
---|
505 | ! IF (ICB.GT.1) ATH=TH(ICB-1)-5.0 ! modif sb |
---|
506 | IF (tdif<dtcrit .OR. ath>ath1) THEN |
---|
507 | DO i = 1, nl |
---|
508 | sig(i) = beta * sig(i) - 2. * alpha * tdif * tdif |
---|
509 | sig(i) = amax1(sig(i), 0.0) |
---|
510 | w0(i) = beta * w0(i) |
---|
511 | END DO |
---|
512 | iflag = 0 |
---|
513 | RETURN |
---|
514 | END IF |
---|
515 | |
---|
516 | |
---|
517 | |
---|
518 | ! *** IF THIS POINT IS REACHED, MOIST CONVECTIVE ADJUSTMENT IS NECESSARY |
---|
519 | ! *** |
---|
520 | ! *** NOW INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS |
---|
521 | ! *** |
---|
522 | |
---|
523 | DO i = 1, nl |
---|
524 | hp(i) = h(i) |
---|
525 | wt(i) = 0.001 |
---|
526 | rp(i) = rr(i) |
---|
527 | up(i) = u(i) |
---|
528 | vp(i) = v(i) |
---|
529 | DO j = 1, ntra |
---|
530 | trap(i, j) = tra(i, j) |
---|
531 | END DO |
---|
532 | nent(i) = 0 |
---|
533 | water(i) = 0.0 |
---|
534 | evap(i) = 0.0 |
---|
535 | b(i) = 0.0 |
---|
536 | mp(i) = 0.0 |
---|
537 | m(i) = 0.0 |
---|
538 | lv(i) = alv0 - cpvmcl * (t(i) - 273.15) |
---|
539 | lvcp(i) = lv(i) / cpn(i) |
---|
540 | DO j = 1, nl |
---|
541 | qent(i, j) = rr(j) |
---|
542 | elij(i, j) = 0.0 |
---|
543 | ment(i, j) = 0.0 |
---|
544 | sij(i, j) = 0.0 |
---|
545 | uent(i, j) = u(j) |
---|
546 | vent(i, j) = v(j) |
---|
547 | DO k = 1, ntra |
---|
548 | traent(i, j, k) = tra(j, k) |
---|
549 | END DO |
---|
550 | END DO |
---|
551 | END DO |
---|
552 | |
---|
553 | delti = 1.0 / delt |
---|
554 | |
---|
555 | ! *** FIND THE FIRST MODEL LEVEL (INB) ABOVE THE PARCEL'S *** |
---|
556 | ! *** LEVEL OF NEUTRAL BUOYANCY *** |
---|
557 | |
---|
558 | inb = nl - 1 |
---|
559 | DO i = icb, nl - 1 |
---|
560 | ! jyg |
---|
561 | ! CC IF((TVP(I)-TV(I)).LT.DTCRIT)THEN |
---|
562 | IF (buoy(i)<dtovsh) THEN |
---|
563 | inb = min(inb, i) |
---|
564 | END IF |
---|
565 | END DO |
---|
566 | |
---|
567 | |
---|
568 | |
---|
569 | |
---|
570 | |
---|
571 | ! *** RESET SIG(I) AND W0(I) FOR I>INB AND I<ICB *** |
---|
572 | |
---|
573 | IF (inb<(nl - 1)) THEN |
---|
574 | DO i = inb + 1, nl - 1 |
---|
575 | ! jyg |
---|
576 | ! CC SIG(I)=BETA*SIG(I)-2.0E-4*ALPHA*(TV(INB)-TVP(INB))* |
---|
577 | ! CC 1 ABS(TV(INB)-TVP(INB)) |
---|
578 | sig(i) = beta * sig(i) + 2. * alpha * buoy(inb) * abs(buoy(inb)) |
---|
579 | sig(i) = amax1(sig(i), 0.0) |
---|
580 | w0(i) = beta * w0(i) |
---|
581 | END DO |
---|
582 | END IF |
---|
583 | DO i = 1, icb |
---|
584 | ! jyg |
---|
585 | ! CC SIG(I)=BETA*SIG(I)-2.0E-4*ALPHA*(TV(ICB)-TVP(ICB))* |
---|
586 | ! CC 1 (TV(ICB)-TVP(ICB)) |
---|
587 | sig(i) = beta * sig(i) - 2. * alpha * buoy(icb) * buoy(icb) |
---|
588 | sig(i) = amax1(sig(i), 0.0) |
---|
589 | w0(i) = beta * w0(i) |
---|
590 | END DO |
---|
591 | |
---|
592 | ! *** RESET FRACTIONAL AREAS OF UPDRAFTS AND W0 AT INITIAL TIME *** |
---|
593 | ! *** AND AFTER 10 TIME STEPS OF NO CONVECTION *** |
---|
594 | |
---|
595 | IF (sig(nd)<1.5 .OR. sig(nd)>12.0) THEN |
---|
596 | DO i = 1, nl - 1 |
---|
597 | sig(i) = 0.0 |
---|
598 | w0(i) = 0.0 |
---|
599 | END DO |
---|
600 | END IF |
---|
601 | |
---|
602 | ! *** CALCULATE LIQUID WATER STATIC ENERGY OF LIFTED PARCEL *** |
---|
603 | |
---|
604 | DO i = icb, inb |
---|
605 | hp(i) = h(1) + (lv(i) + (cpd - cpv) * t(i)) * ep(i) * clw(i) |
---|
606 | END DO |
---|
607 | |
---|
608 | ! *** CALCULATE CONVECTIVE AVAILABLE POTENTIAL ENERGY (CAPE), *** |
---|
609 | ! *** VERTICAL VELOCITY (W), FRACTIONAL AREA COVERED BY *** |
---|
610 | ! *** UNDILUTE UPDRAFT (SIG), AND UPDRAFT MASS FLUX (M) *** |
---|
611 | |
---|
612 | cape = 0.0 |
---|
613 | |
---|
614 | DO i = icb + 1, inb |
---|
615 | ! jyg1 |
---|
616 | ! CC CAPE=CAPE+RD*(TVP(I-1)-TV(I-1))*(PH(I-1)-PH(I))/P(I-1) |
---|
617 | ! CC DCAPE=RD*BUOY(I-1)*(PH(I-1)-PH(I))/P(I-1) |
---|
618 | ! CC DLNP=(PH(I-1)-PH(I))/P(I-1) |
---|
619 | ! The interval on which CAPE is computed starts at PBASE : |
---|
620 | deltap = min(pbase, ph(i - 1)) - min(pbase, ph(i)) |
---|
621 | cape = cape + rd * buoy(i - 1) * deltap / p(i - 1) |
---|
622 | dcape = rd * buoy(i - 1) * deltap / p(i - 1) |
---|
623 | dlnp = deltap / p(i - 1) |
---|
624 | ! jyg2 |
---|
625 | ! sb3d PRINT *,'buoy,dlnp,dcape,cape',buoy(i-1),dlnp,dcape,cape |
---|
626 | ! test sb: |
---|
627 | ! @ WRITE(*,*) '############################################' |
---|
628 | ! @ WRITE(*,*) 'cape,rrd,buoy,deltap,p,pbase,ph:' |
---|
629 | ! @ : ,cape,rd,buoy(i-1),deltap,p(i-1),pbase,ph(i) |
---|
630 | ! @ WRITE(*,*) '############################################' |
---|
631 | |
---|
632 | ! fin test sb |
---|
633 | cape = amax1(0.0, cape) |
---|
634 | |
---|
635 | sigold = sig(i) |
---|
636 | dtmin = 100.0 |
---|
637 | DO j = icb, i - 1 |
---|
638 | ! jyg |
---|
639 | ! CC DTMIN=AMIN1(DTMIN,(TVP(J)-TV(J))) |
---|
640 | dtmin = amin1(dtmin, buoy(j)) |
---|
641 | END DO |
---|
642 | ! sb3d PRINT *, 'DTMIN, BETA, ALPHA, SIG = ',DTMIN,BETA,ALPHA,SIG(I) |
---|
643 | sig(i) = beta * sig(i) + alpha * dtmin * abs(dtmin) |
---|
644 | sig(i) = amax1(sig(i), 0.0) |
---|
645 | sig(i) = amin1(sig(i), 0.01) |
---|
646 | fac = amin1(((dtcrit - dtmin) / dtcrit), 1.0) |
---|
647 | ! jyg |
---|
648 | ! C Essais de reduction de la vitesse |
---|
649 | ! C FAC = FAC*.5 |
---|
650 | |
---|
651 | w = (1. - beta) * fac * sqrt(cape) + beta * w0(i) |
---|
652 | amu = 0.5 * (sig(i) + sigold) * w |
---|
653 | m(i) = amu * 0.007 * p(i) * (ph(i) - ph(i + 1)) / tv(i) |
---|
654 | |
---|
655 | ! --------- test sb: |
---|
656 | ! WRITE(*,*) '############################################' |
---|
657 | ! WRITE(*,*) 'k,amu,buoy(k-1),deltap,w,beta,fac,cape,w0(k)' |
---|
658 | ! WRITE(*,*) i,amu,buoy(i-1),deltap |
---|
659 | ! : ,w,beta,fac,cape,w0(i) |
---|
660 | ! WRITE(*,*) '############################################' |
---|
661 | ! --------- |
---|
662 | |
---|
663 | w0(i) = w |
---|
664 | END DO |
---|
665 | w0(icb) = 0.5 * w0(icb + 1) |
---|
666 | m(icb) = 0.5 * m(icb + 1) * (ph(icb) - ph(icb + 1)) / (ph(icb + 1) - ph(icb + 2)) |
---|
667 | sig(icb) = sig(icb + 1) |
---|
668 | sig(icb - 1) = sig(icb) |
---|
669 | ! jyg1 |
---|
670 | ! sb3d PRINT *, 'Cloud base, c. top, CAPE',ICB,INB,cape |
---|
671 | ! sb3d PRINT *, 'SIG ',(sig(i),i=1,inb) |
---|
672 | ! sb3d PRINT *, 'W ',(w0(i),i=1,inb) |
---|
673 | ! sb3d PRINT *, 'M ',(m(i), i=1,inb) |
---|
674 | ! sb3d PRINT *, 'Dt1 ',(tvp(i)-tv(i),i=1,inb) |
---|
675 | ! sb3d PRINT *, 'Dt_vrai ',(buoy(i),i=1,inb) |
---|
676 | ! jyg2 |
---|
677 | |
---|
678 | ! *** CALCULATE ENTRAINED AIR MASS FLUX (MENT), TOTAL WATER MIXING *** |
---|
679 | ! *** RATIO (QENT), TOTAL CONDENSED WATER (ELIJ), AND MIXING *** |
---|
680 | ! *** FRACTION (SIJ) *** |
---|
681 | |
---|
682 | DO i = icb, inb |
---|
683 | rti = rr(1) - ep(i) * clw(i) |
---|
684 | DO j = icb - 1, inb |
---|
685 | bf2 = 1. + lv(j) * lv(j) * rs(j) / (rv * t(j) * t(j) * cpd) |
---|
686 | anum = h(j) - hp(i) + (cpv - cpd) * t(j) * (rti - rr(j)) |
---|
687 | denom = h(i) - hp(i) + (cpd - cpv) * (rr(i) - rti) * t(j) |
---|
688 | dei = denom |
---|
689 | IF (abs(dei)<0.01) dei = 0.01 |
---|
690 | sij(i, j) = anum / dei |
---|
691 | sij(i, i) = 1.0 |
---|
692 | altem = sij(i, j) * rr(i) + (1. - sij(i, j)) * rti - rs(j) |
---|
693 | altem = altem / bf2 |
---|
694 | cwat = clw(j) * (1. - ep(j)) |
---|
695 | stemp = sij(i, j) |
---|
696 | IF ((stemp<0.0 .OR. stemp>1.0 .OR. altem>cwat) .AND. j>i) THEN |
---|
697 | anum = anum - lv(j) * (rti - rs(j) - cwat * bf2) |
---|
698 | denom = denom + lv(j) * (rr(i) - rti) |
---|
699 | IF (abs(denom)<0.01) denom = 0.01 |
---|
700 | sij(i, j) = anum / denom |
---|
701 | altem = sij(i, j) * rr(i) + (1. - sij(i, j)) * rti - rs(j) |
---|
702 | altem = altem - (bf2 - 1.) * cwat |
---|
703 | END IF |
---|
704 | |
---|
705 | IF (sij(i, j)>0.0 .AND. sij(i, j)<0.95) THEN |
---|
706 | qent(i, j) = sij(i, j) * rr(i) + (1. - sij(i, j)) * rti |
---|
707 | uent(i, j) = sij(i, j) * u(i) + (1. - sij(i, j)) * u(nk) |
---|
708 | vent(i, j) = sij(i, j) * v(i) + (1. - sij(i, j)) * v(nk) |
---|
709 | DO k = 1, ntra |
---|
710 | traent(i, j, k) = sij(i, j) * tra(i, k) + (1. - sij(i, j)) * tra(nk, k) |
---|
711 | END DO |
---|
712 | elij(i, j) = altem |
---|
713 | elij(i, j) = amax1(0.0, elij(i, j)) |
---|
714 | ment(i, j) = m(i) / (1. - sij(i, j)) |
---|
715 | nent(i) = nent(i) + 1 |
---|
716 | END IF |
---|
717 | sij(i, j) = amax1(0.0, sij(i, j)) |
---|
718 | sij(i, j) = amin1(1.0, sij(i, j)) |
---|
719 | END DO |
---|
720 | |
---|
721 | ! *** IF NO AIR CAN ENTRAIN AT LEVEL I ASSUME THAT UPDRAFT DETRAINS |
---|
722 | ! *** |
---|
723 | ! *** AT THAT LEVEL AND CALCULATE DETRAINED AIR FLUX AND PROPERTIES |
---|
724 | ! *** |
---|
725 | |
---|
726 | IF (nent(i)==0) THEN |
---|
727 | ment(i, i) = m(i) |
---|
728 | qent(i, i) = rr(nk) - ep(i) * clw(i) |
---|
729 | uent(i, i) = u(nk) |
---|
730 | vent(i, i) = v(nk) |
---|
731 | DO j = 1, ntra |
---|
732 | traent(i, i, j) = tra(nk, j) |
---|
733 | END DO |
---|
734 | elij(i, i) = clw(i) |
---|
735 | sij(i, i) = 1.0 |
---|
736 | END IF |
---|
737 | |
---|
738 | DO j = icb - 1, inb |
---|
739 | sigij(i, j) = sij(i, j) |
---|
740 | END DO |
---|
741 | |
---|
742 | END DO |
---|
743 | |
---|
744 | ! *** NORMALIZE ENTRAINED AIR MASS FLUXES TO REPRESENT EQUAL *** |
---|
745 | ! *** PROBABILITIES OF MIXING *** |
---|
746 | |
---|
747 | DO i = icb, inb |
---|
748 | IF (nent(i)/=0) THEN |
---|
749 | qp = rr(1) - ep(i) * clw(i) |
---|
750 | anum = h(i) - hp(i) - lv(i) * (qp - rs(i)) + (cpv - cpd) * t(i) * (qp - rr(i)) |
---|
751 | denom = h(i) - hp(i) + lv(i) * (rr(i) - qp) + (cpd - cpv) * t(i) * (rr(i) - qp) |
---|
752 | IF (abs(denom)<0.01) denom = 0.01 |
---|
753 | scrit = anum / denom |
---|
754 | alt = qp - rs(i) + scrit * (rr(i) - qp) |
---|
755 | IF (scrit<=0.0 .OR. alt<=0.0) scrit = 1.0 |
---|
756 | smax = 0.0 |
---|
757 | asij = 0.0 |
---|
758 | DO j = inb, icb - 1, -1 |
---|
759 | IF (sij(i, j)>1.0E-16 .AND. sij(i, j)<0.95) THEN |
---|
760 | wgh = 1.0 |
---|
761 | IF (j>i) THEN |
---|
762 | sjmax = amax1(sij(i, j + 1), smax) |
---|
763 | sjmax = amin1(sjmax, scrit) |
---|
764 | smax = amax1(sij(i, j), smax) |
---|
765 | sjmin = amax1(sij(i, j - 1), smax) |
---|
766 | sjmin = amin1(sjmin, scrit) |
---|
767 | IF (sij(i, j)<(smax - 1.0E-16)) wgh = 0.0 |
---|
768 | smid = amin1(sij(i, j), scrit) |
---|
769 | ELSE |
---|
770 | sjmax = amax1(sij(i, j + 1), scrit) |
---|
771 | smid = amax1(sij(i, j), scrit) |
---|
772 | sjmin = 0.0 |
---|
773 | IF (j>1) sjmin = sij(i, j - 1) |
---|
774 | sjmin = amax1(sjmin, scrit) |
---|
775 | END IF |
---|
776 | delp = abs(sjmax - smid) |
---|
777 | delm = abs(sjmin - smid) |
---|
778 | asij = asij + wgh * (delp + delm) |
---|
779 | ment(i, j) = ment(i, j) * (delp + delm) * wgh |
---|
780 | END IF |
---|
781 | END DO |
---|
782 | asij = amax1(1.0E-16, asij) |
---|
783 | asij = 1.0 / asij |
---|
784 | DO j = icb - 1, inb |
---|
785 | ment(i, j) = ment(i, j) * asij |
---|
786 | END DO |
---|
787 | asum = 0.0 |
---|
788 | bsum = 0.0 |
---|
789 | DO j = icb - 1, inb |
---|
790 | asum = asum + ment(i, j) |
---|
791 | ment(i, j) = ment(i, j) * sig(j) |
---|
792 | bsum = bsum + ment(i, j) |
---|
793 | END DO |
---|
794 | bsum = amax1(bsum, 1.0E-16) |
---|
795 | bsum = 1.0 / bsum |
---|
796 | DO j = icb - 1, inb |
---|
797 | ment(i, j) = ment(i, j) * asum * bsum |
---|
798 | END DO |
---|
799 | csum = 0.0 |
---|
800 | DO j = icb - 1, inb |
---|
801 | csum = csum + ment(i, j) |
---|
802 | END DO |
---|
803 | |
---|
804 | IF (csum<m(i)) THEN |
---|
805 | nent(i) = 0 |
---|
806 | ment(i, i) = m(i) |
---|
807 | qent(i, i) = rr(1) - ep(i) * clw(i) |
---|
808 | uent(i, i) = u(nk) |
---|
809 | vent(i, i) = v(nk) |
---|
810 | DO j = 1, ntra |
---|
811 | traent(i, i, j) = tra(nk, j) |
---|
812 | END DO |
---|
813 | elij(i, i) = clw(i) |
---|
814 | sij(i, i) = 1.0 |
---|
815 | END IF |
---|
816 | END IF |
---|
817 | END DO |
---|
818 | |
---|
819 | |
---|
820 | |
---|
821 | ! ************************************************************** |
---|
822 | ! * CALCUL DES MENTS(I,J) ET DES QENTS(I,J) |
---|
823 | ! ************************************************************* |
---|
824 | |
---|
825 | DO im = 1, nd |
---|
826 | DO jm = 1, nd |
---|
827 | |
---|
828 | qents(im, jm) = qent(im, jm) |
---|
829 | ments(im, jm) = ment(im, jm) |
---|
830 | END DO |
---|
831 | END DO |
---|
832 | |
---|
833 | ! ********************************************************** |
---|
834 | ! --- test sb: |
---|
835 | ! @ WRITE(*,*) '^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^' |
---|
836 | ! @ WRITE(*,*) 'inb,m(inb),ment(inb,inb),sigij(inb,inb):' |
---|
837 | ! @ WRITE(*,*) inb,m(inb),ment(inb,inb),sigij(inb,inb) |
---|
838 | ! @ WRITE(*,*) '^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^' |
---|
839 | ! --- |
---|
840 | |
---|
841 | |
---|
842 | |
---|
843 | |
---|
844 | |
---|
845 | ! *** CHECK WHETHER EP(INB)=0, IF SO, SKIP PRECIPITATING *** |
---|
846 | ! *** DOWNDRAFT CALCULATION *** |
---|
847 | |
---|
848 | IF (ep(inb)<0.0001) GO TO 405 |
---|
849 | |
---|
850 | ! *** INTEGRATE LIQUID WATER EQUATION TO FIND CONDENSED WATER *** |
---|
851 | ! *** AND CONDENSED WATER FLUX *** |
---|
852 | |
---|
853 | wflux = 0.0 |
---|
854 | tinv = 1. / 3. |
---|
855 | |
---|
856 | ! *** BEGIN DOWNDRAFT LOOP *** |
---|
857 | |
---|
858 | DO i = inb, 1, -1 |
---|
859 | |
---|
860 | ! *** CALCULATE DETRAINED PRECIPITATION *** |
---|
861 | |
---|
862 | wdtrain = 10.0 * ep(i) * m(i) * clw(i) |
---|
863 | IF (i>1) THEN |
---|
864 | DO j = 1, i - 1 |
---|
865 | awat = elij(j, i) - (1. - ep(i)) * clw(i) |
---|
866 | awat = amax1(awat, 0.0) |
---|
867 | wdtrain = wdtrain + 10.0 * awat * ment(j, i) |
---|
868 | END DO |
---|
869 | END IF |
---|
870 | |
---|
871 | ! *** FIND RAIN WATER AND EVAPORATION USING PROVISIONAL *** |
---|
872 | ! *** ESTIMATES OF RP(I)AND RP(I-1) *** |
---|
873 | |
---|
874 | wt(i) = 45.0 |
---|
875 | IF (i<inb) THEN |
---|
876 | rp(i) = rp(i + 1) + (cpd * (t(i + 1) - t(i)) + gz(i + 1) - gz(i)) / lv(i) |
---|
877 | rp(i) = 0.5 * (rp(i) + rr(i)) |
---|
878 | END IF |
---|
879 | rp(i) = amax1(rp(i), 0.0) |
---|
880 | rp(i) = amin1(rp(i), rs(i)) |
---|
881 | rp(inb) = rr(inb) |
---|
882 | IF (i==1) THEN |
---|
883 | afac = p(1) * (rs(1) - rp(1)) / (1.0E4 + 2000.0 * p(1) * rs(1)) |
---|
884 | ELSE |
---|
885 | rp(i - 1) = rp(i) + (cpd * (t(i) - t(i - 1)) + gz(i) - gz(i - 1)) / lv(i) |
---|
886 | rp(i - 1) = 0.5 * (rp(i - 1) + rr(i - 1)) |
---|
887 | rp(i - 1) = amin1(rp(i - 1), rs(i - 1)) |
---|
888 | rp(i - 1) = amax1(rp(i - 1), 0.0) |
---|
889 | afac1 = p(i) * (rs(i) - rp(i)) / (1.0E4 + 2000.0 * p(i) * rs(i)) |
---|
890 | afac2 = p(i - 1) * (rs(i - 1) - rp(i - 1)) / (1.0E4 + 2000.0 * p(i - 1) * rs(i - 1)) |
---|
891 | afac = 0.5 * (afac1 + afac2) |
---|
892 | END IF |
---|
893 | IF (i==inb) afac = 0.0 |
---|
894 | afac = amax1(afac, 0.0) |
---|
895 | bfac = 1. / (sigd * wt(i)) |
---|
896 | |
---|
897 | ! jyg1 |
---|
898 | ! CC SIGT=1.0 |
---|
899 | ! CC IF(I.GE.ICB)SIGT=SIGP(I) |
---|
900 | ! Prise en compte de la variation progressive de SIGT dans |
---|
901 | ! les couches ICB et ICB-1: |
---|
902 | ! Pour PLCL<PH(I+1), PR1=0 & PR2=1 |
---|
903 | ! Pour PLCL>PH(I), PR1=1 & PR2=0 |
---|
904 | ! Pour PH(I+1)<PLCL<PH(I), PR1 est la proportion a cheval |
---|
905 | ! sur le nuage, et PR2 est la proportion sous la base du |
---|
906 | ! nuage. |
---|
907 | pr1 = (plcl - ph(i + 1)) / (ph(i) - ph(i + 1)) |
---|
908 | pr1 = max(0., min(1., pr1)) |
---|
909 | pr2 = (ph(i) - plcl) / (ph(i) - ph(i + 1)) |
---|
910 | pr2 = max(0., min(1., pr2)) |
---|
911 | sigt = sigp(i) * pr1 + pr2 |
---|
912 | ! sb3d PRINT *,'i,sigt,pr1,pr2', i,sigt,pr1,pr2 |
---|
913 | ! jyg2 |
---|
914 | |
---|
915 | b6 = bfac * 50. * sigd * (ph(i) - ph(i + 1)) * sigt * afac |
---|
916 | c6 = water(i + 1) + bfac * wdtrain - 50. * sigd * bfac * (ph(i) - ph(i + 1)) * evap(i + 1) |
---|
917 | IF (c6>0.0) THEN |
---|
918 | revap = 0.5 * (-b6 + sqrt(b6 * b6 + 4. * c6)) |
---|
919 | evap(i) = sigt * afac * revap |
---|
920 | water(i) = revap * revap |
---|
921 | ELSE |
---|
922 | evap(i) = -evap(i + 1) + 0.02 * (wdtrain + sigd * wt(i) * water(i + 1)) / (sigd * (ph(i & |
---|
923 | ) - ph(i + 1))) |
---|
924 | END IF |
---|
925 | |
---|
926 | |
---|
927 | |
---|
928 | ! *** CALCULATE PRECIPITATING DOWNDRAFT MASS FLUX UNDER *** |
---|
929 | ! *** HYDROSTATIC APPROXIMATION *** |
---|
930 | |
---|
931 | IF (i==1) GO TO 360 |
---|
932 | tevap = amax1(0.0, evap(i)) |
---|
933 | delth = amax1(0.001, (th(i) - th(i - 1))) |
---|
934 | mp(i) = 10. * lvcp(i) * sigd * tevap * (p(i - 1) - p(i)) / delth |
---|
935 | |
---|
936 | ! *** IF HYDROSTATIC ASSUMPTION FAILS, *** |
---|
937 | ! *** SOLVE CUBIC DIFFERENCE EQUATION FOR DOWNDRAFT THETA *** |
---|
938 | ! *** AND MASS FLUX FROM TWO SIMULTANEOUS DIFFERENTIAL EQNS *** |
---|
939 | |
---|
940 | amfac = sigd * sigd * 70.0 * ph(i) * (p(i - 1) - p(i)) * (th(i) - th(i - 1)) / (tv(i) * th(i)) |
---|
941 | amp2 = abs(mp(i + 1) * mp(i + 1) - mp(i) * mp(i)) |
---|
942 | IF (amp2>(0.1 * amfac)) THEN |
---|
943 | xf = 100.0 * sigd * sigd * sigd * (ph(i) - ph(i + 1)) |
---|
944 | tf = b(i) - 5.0 * (th(i) - th(i - 1)) * t(i) / (lvcp(i) * sigd * th(i)) |
---|
945 | af = xf * tf + mp(i + 1) * mp(i + 1) * tinv |
---|
946 | bf = 2. * (tinv * mp(i + 1))**3 + tinv * mp(i + 1) * xf * tf + & |
---|
947 | 50. * (p(i - 1) - p(i)) * xf * tevap |
---|
948 | fac2 = 1.0 |
---|
949 | IF (bf<0.0) fac2 = -1.0 |
---|
950 | bf = abs(bf) |
---|
951 | ur = 0.25 * bf * bf - af * af * af * tinv * tinv * tinv |
---|
952 | IF (ur>=0.0) THEN |
---|
953 | sru = sqrt(ur) |
---|
954 | fac = 1.0 |
---|
955 | IF ((0.5 * bf - sru)<0.0) fac = -1.0 |
---|
956 | mp(i) = mp(i + 1) * tinv + (0.5 * bf + sru)**tinv + & |
---|
957 | fac * (abs(0.5 * bf - sru))**tinv |
---|
958 | ELSE |
---|
959 | d = atan(2. * sqrt(-ur) / (bf + 1.0E-28)) |
---|
960 | IF (fac2<0.0) d = 3.14159 - d |
---|
961 | mp(i) = mp(i + 1) * tinv + 2. * sqrt(af * tinv) * cos(d * tinv) |
---|
962 | END IF |
---|
963 | mp(i) = amax1(0.0, mp(i)) |
---|
964 | b(i - 1) = b(i) + 100.0 * (p(i - 1) - p(i)) * tevap / (mp(i) + sigd * 0.1) - & |
---|
965 | 10.0 * (th(i) - th(i - 1)) * t(i) / (lvcp(i) * sigd * th(i)) |
---|
966 | b(i - 1) = amax1(b(i - 1), 0.0) |
---|
967 | END IF |
---|
968 | |
---|
969 | |
---|
970 | |
---|
971 | ! *** LIMIT MAGNITUDE OF MP(I) TO MEET CFL CONDITION *** |
---|
972 | |
---|
973 | ampmax = 2.0 * (ph(i) - ph(i + 1)) * delti |
---|
974 | amp2 = 2.0 * (ph(i - 1) - ph(i)) * delti |
---|
975 | ampmax = amin1(ampmax, amp2) |
---|
976 | mp(i) = amin1(mp(i), ampmax) |
---|
977 | |
---|
978 | ! *** FORCE MP TO DECREASE LINEARLY TO ZERO *** |
---|
979 | ! *** BETWEEN CLOUD BASE AND THE SURFACE *** |
---|
980 | |
---|
981 | IF (p(i)>p(icb)) THEN |
---|
982 | mp(i) = mp(icb) * (p(1) - p(i)) / (p(1) - p(icb)) |
---|
983 | END IF |
---|
984 | 360 CONTINUE |
---|
985 | |
---|
986 | ! *** FIND MIXING RATIO OF PRECIPITATING DOWNDRAFT *** |
---|
987 | |
---|
988 | IF (i==inb) GO TO 400 |
---|
989 | rp(i) = rr(i) |
---|
990 | IF (mp(i)>mp(i + 1)) THEN |
---|
991 | rp(i) = rp(i + 1) * mp(i + 1) + rr(i) * (mp(i) - mp(i + 1)) + & |
---|
992 | 5. * sigd * (ph(i) - ph(i + 1)) * (evap(i + 1) + evap(i)) |
---|
993 | rp(i) = rp(i) / mp(i) |
---|
994 | up(i) = up(i + 1) * mp(i + 1) + u(i) * (mp(i) - mp(i + 1)) |
---|
995 | up(i) = up(i) / mp(i) |
---|
996 | vp(i) = vp(i + 1) * mp(i + 1) + v(i) * (mp(i) - mp(i + 1)) |
---|
997 | vp(i) = vp(i) / mp(i) |
---|
998 | DO j = 1, ntra |
---|
999 | trap(i, j) = trap(i + 1, j) * mp(i + 1) + trap(i, j) * (mp(i) - mp(i + 1)) |
---|
1000 | trap(i, j) = trap(i, j) / mp(i) |
---|
1001 | END DO |
---|
1002 | ELSE |
---|
1003 | IF (mp(i + 1)>1.0E-16) THEN |
---|
1004 | rp(i) = rp(i + 1) + 5.0 * sigd * (ph(i) - ph(i + 1)) * (evap(i + 1) + evap(i)) / mp(i + 1 & |
---|
1005 | ) |
---|
1006 | up(i) = up(i + 1) |
---|
1007 | vp(i) = vp(i + 1) |
---|
1008 | DO j = 1, ntra |
---|
1009 | trap(i, j) = trap(i + 1, j) |
---|
1010 | END DO |
---|
1011 | END IF |
---|
1012 | END IF |
---|
1013 | rp(i) = amin1(rp(i), rs(i)) |
---|
1014 | rp(i) = amax1(rp(i), 0.0) |
---|
1015 | 400 END DO |
---|
1016 | |
---|
1017 | ! *** CALCULATE SURFACE PRECIPITATION IN MM/DAY *** |
---|
1018 | |
---|
1019 | precip = wt(1) * sigd * water(1) * 8640.0 |
---|
1020 | |
---|
1021 | ! sb *** Calculate downdraft velocity scale and surface temperature and |
---|
1022 | ! *** |
---|
1023 | ! sb *** water vapor fluctuations |
---|
1024 | ! *** |
---|
1025 | ! sb (inspire de convect 4.3) |
---|
1026 | |
---|
1027 | ! BETAD=10.0 |
---|
1028 | betad = 5.0 |
---|
1029 | wd = betad * abs(mp(icb)) * 0.01 * rd * t(icb) / (sigd * p(icb)) |
---|
1030 | |
---|
1031 | 405 CONTINUE |
---|
1032 | |
---|
1033 | ! *** CALCULATE TENDENCIES OF LOWEST LEVEL POTENTIAL TEMPERATURE *** |
---|
1034 | ! *** AND MIXING RATIO *** |
---|
1035 | |
---|
1036 | dpinv = 1.0 / (ph(1) - ph(2)) |
---|
1037 | am = 0.0 |
---|
1038 | DO k = 2, inb |
---|
1039 | am = am + m(k) |
---|
1040 | END DO |
---|
1041 | IF ((0.1 * dpinv * am)>=delti) iflag = 4 |
---|
1042 | ft(1) = 0.1 * dpinv * am * (t(2) - t(1) + (gz(2) - gz(1)) / cpn(1)) |
---|
1043 | ft(1) = ft(1) - 0.5 * lvcp(1) * sigd * (evap(1) + evap(2)) |
---|
1044 | ft(1) = ft(1) - 0.09 * sigd * mp(2) * t(1) * b(1) * dpinv |
---|
1045 | ft(1) = ft(1) + 0.01 * sigd * wt(1) * (cl - cpd) * water(2) * (t(2) - t(1)) * dpinv / cpn(1) |
---|
1046 | fr(1) = 0.1 * mp(2) * (rp(2) - rr(1)) * & ! correction bug conservation eau |
---|
1047 | ! 1 DPINV+SIGD*0.5*(EVAP(1)+EVAP(2)) |
---|
1048 | dpinv + sigd * 0.5 * (evap(1) + evap(2)) |
---|
1049 | ! IM cf. SBL |
---|
1050 | ! 1 DPINV+SIGD*EVAP(1) |
---|
1051 | fr(1) = fr(1) + 0.1 * am * (rr(2) - rr(1)) * dpinv |
---|
1052 | fu(1) = fu(1) + 0.1 * dpinv * (mp(2) * (up(2) - u(1)) + am * (u(2) - u(1))) |
---|
1053 | fv(1) = fv(1) + 0.1 * dpinv * (mp(2) * (vp(2) - v(1)) + am * (v(2) - v(1))) |
---|
1054 | DO j = 1, ntra |
---|
1055 | ftra(1, j) = ftra(1, j) + 0.1 * dpinv * (mp(2) * (trap(2, j) - tra(1, & |
---|
1056 | j)) + am * (tra(2, j) - tra(1, j))) |
---|
1057 | END DO |
---|
1058 | amde = 0.0 |
---|
1059 | DO j = 2, inb |
---|
1060 | fr(1) = fr(1) + 0.1 * dpinv * ment(j, 1) * (qent(j, 1) - rr(1)) |
---|
1061 | fu(1) = fu(1) + 0.1 * dpinv * ment(j, 1) * (uent(j, 1) - u(1)) |
---|
1062 | fv(1) = fv(1) + 0.1 * dpinv * ment(j, 1) * (vent(j, 1) - v(1)) |
---|
1063 | DO k = 1, ntra |
---|
1064 | ftra(1, k) = ftra(1, k) + 0.1 * dpinv * ment(j, 1) * (traent(j, 1, k) - tra(1, k)) |
---|
1065 | END DO |
---|
1066 | END DO |
---|
1067 | |
---|
1068 | ! *** CALCULATE TENDENCIES OF POTENTIAL TEMPERATURE AND MIXING RATIO *** |
---|
1069 | ! *** AT LEVELS ABOVE THE LOWEST LEVEL *** |
---|
1070 | |
---|
1071 | ! *** FIRST FIND THE NET SATURATED UPDRAFT AND DOWNDRAFT MASS FLUXES *** |
---|
1072 | ! *** THROUGH EACH LEVEL *** |
---|
1073 | |
---|
1074 | DO i = 2, inb |
---|
1075 | dpinv = 1.0 / (ph(i) - ph(i + 1)) |
---|
1076 | cpinv = 1.0 / cpn(i) |
---|
1077 | amp1 = 0.0 |
---|
1078 | DO k = i + 1, inb + 1 |
---|
1079 | amp1 = amp1 + m(k) |
---|
1080 | END DO |
---|
1081 | DO k = 1, i |
---|
1082 | DO j = i + 1, inb + 1 |
---|
1083 | amp1 = amp1 + ment(k, j) |
---|
1084 | END DO |
---|
1085 | END DO |
---|
1086 | IF ((0.1 * dpinv * amp1)>=delti) iflag = 4 |
---|
1087 | ad = 0.0 |
---|
1088 | DO k = 1, i - 1 |
---|
1089 | DO j = i, inb |
---|
1090 | ad = ad + ment(j, k) |
---|
1091 | END DO |
---|
1092 | END DO |
---|
1093 | ft(i) = 0.1 * dpinv * (amp1 * (t(i + 1) - t(i) + (gz(i + 1) - gz(i)) * cpinv) - ad * (t(i) - t(i - & |
---|
1094 | 1) + (gz(i) - gz(i - 1)) * cpinv)) - 0.5 * sigd * lvcp(i) * (evap(i) + evap(i + 1)) |
---|
1095 | rat = cpn(i - 1) * cpinv |
---|
1096 | ft(i) = ft(i) - 0.09 * sigd * (mp(i + 1) * t(i) * b(i) - mp(i) * t(i - 1) * rat * b(i - 1)) * & |
---|
1097 | dpinv |
---|
1098 | ft(i) = ft(i) + 0.1 * dpinv * ment(i, i) * (hp(i) - h(i) + t(i) * (cpv - cpd) * (rr(i) - & |
---|
1099 | qent(i, i))) * cpinv |
---|
1100 | ft(i) = ft(i) + 0.01 * sigd * wt(i) * (cl - cpd) * water(i + 1) * (t(i + 1) - t(i)) * dpinv * & |
---|
1101 | cpinv |
---|
1102 | fr(i) = 0.1 * dpinv * (amp1 * (rr(i + 1) - rr(i)) - ad * (rr(i) - rr(i - 1))) |
---|
1103 | fu(i) = fu(i) + 0.1 * dpinv * (amp1 * (u(i + 1) - u(i)) - ad * (u(i) - u(i - 1))) |
---|
1104 | fv(i) = fv(i) + 0.1 * dpinv * (amp1 * (v(i + 1) - v(i)) - ad * (v(i) - v(i - 1))) |
---|
1105 | DO k = 1, ntra |
---|
1106 | ftra(i, k) = ftra(i, k) + 0.1 * dpinv * (amp1 * (tra(i + 1, k) - tra(i, & |
---|
1107 | k)) - ad * (tra(i, k) - tra(i - 1, k))) |
---|
1108 | END DO |
---|
1109 | DO k = 1, i - 1 |
---|
1110 | awat = elij(k, i) - (1. - ep(i)) * clw(i) |
---|
1111 | awat = amax1(awat, 0.0) |
---|
1112 | fr(i) = fr(i) + 0.1 * dpinv * ment(k, i) * (qent(k, i) - awat - rr(i)) |
---|
1113 | fu(i) = fu(i) + 0.1 * dpinv * ment(k, i) * (uent(k, i) - u(i)) |
---|
1114 | fv(i) = fv(i) + 0.1 * dpinv * ment(k, i) * (vent(k, i) - v(i)) |
---|
1115 | ! (saturated updrafts resulting from mixing) ! cld |
---|
1116 | qcond(i) = qcond(i) + (elij(k, i) - awat) ! cld |
---|
1117 | nqcond(i) = nqcond(i) + 1. ! cld |
---|
1118 | DO j = 1, ntra |
---|
1119 | ftra(i, j) = ftra(i, j) + 0.1 * dpinv * ment(k, i) * (traent(k, i, j) - tra(i, j & |
---|
1120 | )) |
---|
1121 | END DO |
---|
1122 | END DO |
---|
1123 | DO k = i, inb |
---|
1124 | fr(i) = fr(i) + 0.1 * dpinv * ment(k, i) * (qent(k, i) - rr(i)) |
---|
1125 | fu(i) = fu(i) + 0.1 * dpinv * ment(k, i) * (uent(k, i) - u(i)) |
---|
1126 | fv(i) = fv(i) + 0.1 * dpinv * ment(k, i) * (vent(k, i) - v(i)) |
---|
1127 | DO j = 1, ntra |
---|
1128 | ftra(i, j) = ftra(i, j) + 0.1 * dpinv * ment(k, i) * (traent(k, i, j) - tra(i, j & |
---|
1129 | )) |
---|
1130 | END DO |
---|
1131 | END DO |
---|
1132 | fr(i) = fr(i) + 0.5 * sigd * (evap(i) + evap(i + 1)) + 0.1 * (mp(i + 1) * (rp(i + & |
---|
1133 | 1) - rr(i)) - mp(i) * (rp(i) - rr(i - 1))) * dpinv |
---|
1134 | fu(i) = fu(i) + 0.1 * (mp(i + 1) * (up(i + 1) - u(i)) - mp(i) * (up(i) - u(i - 1))) * dpinv |
---|
1135 | fv(i) = fv(i) + 0.1 * (mp(i + 1) * (vp(i + 1) - v(i)) - mp(i) * (vp(i) - v(i - 1))) * dpinv |
---|
1136 | DO j = 1, ntra |
---|
1137 | ftra(i, j) = ftra(i, j) + 0.1 * dpinv * (mp(i + 1) * (trap(i + 1, j) - tra(i, & |
---|
1138 | j)) - mp(i) * (trap(i, j) - trap(i - 1, j))) |
---|
1139 | END DO |
---|
1140 | ! (saturated downdrafts resulting from mixing) ! cld |
---|
1141 | DO k = i + 1, inb ! cld |
---|
1142 | qcond(i) = qcond(i) + elij(k, i) ! cld |
---|
1143 | nqcond(i) = nqcond(i) + 1. ! cld |
---|
1144 | END DO ! cld |
---|
1145 | ! (particular case: no detraining level is found) ! cld |
---|
1146 | IF (nent(i)==0) THEN ! cld |
---|
1147 | qcond(i) = qcond(i) + (1 - ep(i)) * clw(i) ! cld |
---|
1148 | nqcond(i) = nqcond(i) + 1. ! cld |
---|
1149 | END IF ! cld |
---|
1150 | IF (nqcond(i)/=0.) THEN ! cld |
---|
1151 | qcond(i) = qcond(i) / nqcond(i) ! cld |
---|
1152 | END IF ! cld |
---|
1153 | END DO |
---|
1154 | |
---|
1155 | |
---|
1156 | |
---|
1157 | |
---|
1158 | ! *** MOVE THE DETRAINMENT AT LEVEL INB DOWN TO LEVEL INB-1 *** |
---|
1159 | ! *** IN SUCH A WAY AS TO PRESERVE THE VERTICALLY *** |
---|
1160 | ! *** INTEGRATED ENTHALPY AND WATER TENDENCIES *** |
---|
1161 | |
---|
1162 | ! test sb: |
---|
1163 | ! @ WRITE(*,*) '--------------------------------------------' |
---|
1164 | ! @ WRITE(*,*) 'inb,ft,hp,h,t,rr,qent,ment,water,waterp,wt,mp,b' |
---|
1165 | ! @ WRITE(*,*) inb,ft(inb),hp(inb),h(inb) |
---|
1166 | ! @ : ,t(inb),rr(inb),qent(inb,inb) |
---|
1167 | ! @ : ,ment(inb,inb),water(inb) |
---|
1168 | ! @ : ,water(inb+1),wt(inb),mp(inb),b(inb) |
---|
1169 | ! @ WRITE(*,*) '--------------------------------------------' |
---|
1170 | ! fin test sb: |
---|
1171 | |
---|
1172 | ax = 0.1 * ment(inb, inb) * (hp(inb) - h(inb) + t(inb) * (cpv - cpd) * (rr(inb) - qent(inb, & |
---|
1173 | inb))) / (cpn(inb) * (ph(inb) - ph(inb + 1))) |
---|
1174 | ft(inb) = ft(inb) - ax |
---|
1175 | ft(inb - 1) = ft(inb - 1) + ax * cpn(inb) * (ph(inb) - ph(inb + 1)) / (cpn(inb - 1) * (ph(inb & |
---|
1176 | - 1) - ph(inb))) |
---|
1177 | bx = 0.1 * ment(inb, inb) * (qent(inb, inb) - rr(inb)) / (ph(inb) - ph(inb + 1)) |
---|
1178 | fr(inb) = fr(inb) - bx |
---|
1179 | fr(inb - 1) = fr(inb - 1) + bx * (ph(inb) - ph(inb + 1)) / (ph(inb - 1) - ph(inb)) |
---|
1180 | cx = 0.1 * ment(inb, inb) * (uent(inb, inb) - u(inb)) / (ph(inb) - ph(inb + 1)) |
---|
1181 | fu(inb) = fu(inb) - cx |
---|
1182 | fu(inb - 1) = fu(inb - 1) + cx * (ph(inb) - ph(inb + 1)) / (ph(inb - 1) - ph(inb)) |
---|
1183 | dx = 0.1 * ment(inb, inb) * (vent(inb, inb) - v(inb)) / (ph(inb) - ph(inb + 1)) |
---|
1184 | fv(inb) = fv(inb) - dx |
---|
1185 | fv(inb - 1) = fv(inb - 1) + dx * (ph(inb) - ph(inb + 1)) / (ph(inb - 1) - ph(inb)) |
---|
1186 | DO j = 1, ntra |
---|
1187 | ex = 0.1 * ment(inb, inb) * (traent(inb, inb, j) - tra(inb, j)) / & |
---|
1188 | (ph(inb) - ph(inb + 1)) |
---|
1189 | ftra(inb, j) = ftra(inb, j) - ex |
---|
1190 | ftra(inb - 1, j) = ftra(inb - 1, j) + ex * (ph(inb) - ph(inb + 1)) / (ph(inb - 1) - ph(& |
---|
1191 | inb)) |
---|
1192 | END DO |
---|
1193 | |
---|
1194 | ! *** HOMOGINIZE TENDENCIES BELOW CLOUD BASE *** |
---|
1195 | |
---|
1196 | asum = 0.0 |
---|
1197 | bsum = 0.0 |
---|
1198 | csum = 0.0 |
---|
1199 | dsum = 0.0 |
---|
1200 | DO i = 1, icb - 1 |
---|
1201 | asum = asum + ft(i) * (ph(i) - ph(i + 1)) |
---|
1202 | bsum = bsum + fr(i) * (lv(i) + (cl - cpd) * (t(i) - t(1))) * (ph(i) - ph(i + 1)) |
---|
1203 | csum = csum + (lv(i) + (cl - cpd) * (t(i) - t(1))) * (ph(i) - ph(i + 1)) |
---|
1204 | dsum = dsum + t(i) * (ph(i) - ph(i + 1)) / th(i) |
---|
1205 | END DO |
---|
1206 | DO i = 1, icb - 1 |
---|
1207 | ft(i) = asum * t(i) / (th(i) * dsum) |
---|
1208 | fr(i) = bsum / csum |
---|
1209 | END DO |
---|
1210 | |
---|
1211 | ! *** RESET COUNTER AND RETURN *** |
---|
1212 | |
---|
1213 | sig(nd) = 2.0 |
---|
1214 | |
---|
1215 | DO i = 1, nd |
---|
1216 | upwd(i) = 0.0 |
---|
1217 | dnwd(i) = 0.0 |
---|
1218 | ! sb dnwd0(i) = - mp(i) |
---|
1219 | END DO |
---|
1220 | |
---|
1221 | DO i = 1, nl |
---|
1222 | dnwd0(i) = -mp(i) |
---|
1223 | END DO |
---|
1224 | DO i = nl + 1, nd |
---|
1225 | dnwd0(i) = 0. |
---|
1226 | END DO |
---|
1227 | |
---|
1228 | DO i = icb, inb |
---|
1229 | upwd(i) = 0.0 |
---|
1230 | dnwd(i) = 0.0 |
---|
1231 | |
---|
1232 | DO k = i, inb |
---|
1233 | up1 = 0.0 |
---|
1234 | dn1 = 0.0 |
---|
1235 | DO n = 1, i - 1 |
---|
1236 | up1 = up1 + ment(n, k) |
---|
1237 | dn1 = dn1 - ment(k, n) |
---|
1238 | END DO |
---|
1239 | upwd(i) = upwd(i) + m(k) + up1 |
---|
1240 | dnwd(i) = dnwd(i) + dn1 |
---|
1241 | END DO |
---|
1242 | END DO |
---|
1243 | |
---|
1244 | ! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
---|
1245 | ! DETERMINATION DE LA VARIATION DE FLUX ASCENDANT ENTRE |
---|
1246 | ! DEUX NIVEAU NON DILUE Mike |
---|
1247 | ! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
---|
1248 | |
---|
1249 | |
---|
1250 | ! sb do i=1,ND |
---|
1251 | ! sb Mike(i)=M(i) |
---|
1252 | ! sb enddo |
---|
1253 | |
---|
1254 | DO i = 1, nl |
---|
1255 | mike(i) = m(i) |
---|
1256 | END DO |
---|
1257 | DO i = nl + 1, nd |
---|
1258 | mike(i) = 0. |
---|
1259 | END DO |
---|
1260 | |
---|
1261 | DO i = 1, nd |
---|
1262 | ma(i) = 0 |
---|
1263 | END DO |
---|
1264 | |
---|
1265 | ! sb do i=1,nd |
---|
1266 | ! sb do j=i,nd |
---|
1267 | ! sb Ma(i)=Ma(i)+M(j) |
---|
1268 | ! sb enddo |
---|
1269 | ! sb enddo |
---|
1270 | |
---|
1271 | DO i = 1, nl |
---|
1272 | DO j = i, nl |
---|
1273 | ma(i) = ma(i) + m(j) |
---|
1274 | END DO |
---|
1275 | END DO |
---|
1276 | |
---|
1277 | DO i = nl + 1, nd |
---|
1278 | ma(i) = 0. |
---|
1279 | END DO |
---|
1280 | |
---|
1281 | DO i = 1, icb - 1 |
---|
1282 | ma(i) = 0 |
---|
1283 | END DO |
---|
1284 | |
---|
1285 | |
---|
1286 | |
---|
1287 | ! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
---|
1288 | ! ICB REPRESENTE DE NIVEAU OU SE TROUVE LA |
---|
1289 | ! BASE DU NUAGE , ET INB LE TOP DU NUAGE |
---|
1290 | ! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
---|
1291 | |
---|
1292 | DO i = 1, nd |
---|
1293 | mke(i) = upwd(i) + dnwd(i) |
---|
1294 | END DO |
---|
1295 | |
---|
1296 | |
---|
1297 | ! *** Diagnose the in-cloud mixing ratio *** ! cld |
---|
1298 | ! *** of condensed water *** ! cld |
---|
1299 | ! cld |
---|
1300 | DO i = 1, nd ! cld |
---|
1301 | maa(i) = 0.0 ! cld |
---|
1302 | wa(i) = 0.0 ! cld |
---|
1303 | siga(i) = 0.0 ! cld |
---|
1304 | END DO ! cld |
---|
1305 | DO i = nk, inb ! cld |
---|
1306 | DO k = i + 1, inb + 1 ! cld |
---|
1307 | maa(i) = maa(i) + m(k) ! cld |
---|
1308 | END DO ! cld |
---|
1309 | END DO ! cld |
---|
1310 | DO i = icb, inb - 1 ! cld |
---|
1311 | axc(i) = 0. ! cld |
---|
1312 | DO j = icb, i ! cld |
---|
1313 | axc(i) = axc(i) + rd * (tvp(j) - tv(j)) * (ph(j) - ph(j + 1)) / p(j) ! cld |
---|
1314 | END DO ! cld |
---|
1315 | IF (axc(i)>0.0) THEN ! cld |
---|
1316 | wa(i) = sqrt(2. * axc(i)) ! cld |
---|
1317 | END IF ! cld |
---|
1318 | END DO ! cld |
---|
1319 | DO i = 1, nl ! cld |
---|
1320 | IF (wa(i)>0.0) & ! cld |
---|
1321 | siga(i) = maa(i) / wa(i) * rd * tvp(i) / p(i) / 100. / deltac ! cld |
---|
1322 | siga(i) = min(siga(i), 1.0) ! cld |
---|
1323 | qcondc(i) = siga(i) * clw(i) * (1. - ep(i)) & ! cld |
---|
1324 | + (1. - siga(i)) * qcond(i) ! cld |
---|
1325 | END DO ! cld |
---|
1326 | |
---|
1327 | |
---|
1328 | ! @$$cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
---|
1329 | ! @$$ CALL writeg1d(1,klev,ma,'ma ','ma ') |
---|
1330 | ! @$$ CALL writeg1d(1,klev,upwd,'upwd ','upwd ') |
---|
1331 | ! @$$ CALL writeg1d(1,klev,dnwd,'dnwd ','dnwd ') |
---|
1332 | ! @$$ CALL writeg1d(1,klev,dnwd0,'dnwd0 ','dnwd0 ') |
---|
1333 | ! @$$ CALL writeg1d(1,klev,tvp,'tvp ','tvp ') |
---|
1334 | ! @$$ CALL writeg1d(1,klev,tra(1:klev,3),'tra3 ','tra3 ') |
---|
1335 | ! @$$ CALL writeg1d(1,klev,tra(1:klev,4),'tra4 ','tra4 ') |
---|
1336 | ! @$$ CALL writeg1d(1,klev,tra(1:klev,5),'tra5 ','tra5 ') |
---|
1337 | ! @$$ CALL writeg1d(1,klev,tra(1:klev,6),'tra6 ','tra6 ') |
---|
1338 | ! @$$ CALL writeg1d(1,klev,tra(1:klev,7),'tra7 ','tra7 ') |
---|
1339 | ! @$$ CALL writeg1d(1,klev,tra(1:klev,8),'tra8 ','tra8 ') |
---|
1340 | ! @$$ CALL writeg1d(1,klev,tra(1:klev,9),'tra9 ','tra9 ') |
---|
1341 | ! @$$ CALL writeg1d(1,klev,tra(1:klev,10),'tra10','tra10') |
---|
1342 | ! @$$ CALL writeg1d(1,klev,tra(1:klev,11),'tra11','tra11') |
---|
1343 | ! @$$ CALL writeg1d(1,klev,tra(1:klev,12),'tra12','tra12') |
---|
1344 | ! @$$ CALL writeg1d(1,klev,tra(1:klev,13),'tra13','tra13') |
---|
1345 | ! @$$ CALL writeg1d(1,klev,tra(1:klev,14),'tra14','tra14') |
---|
1346 | ! @$$ CALL writeg1d(1,klev,tra(1:klev,15),'tra15','tra15') |
---|
1347 | ! @$$ CALL writeg1d(1,klev,tra(1:klev,16),'tra16','tra16') |
---|
1348 | ! @$$ CALL writeg1d(1,klev,tra(1:klev,17),'tra17','tra17') |
---|
1349 | ! @$$ CALL writeg1d(1,klev,tra(1:klev,18),'tra18','tra18') |
---|
1350 | ! @$$ CALL writeg1d(1,klev,tra(1:klev,19),'tra19','tra19') |
---|
1351 | ! @$$ CALL writeg1d(1,klev,tra(1:klev,20),'tra20','tra20 ') |
---|
1352 | ! @$$ CALL writeg1d(1,klev,trap(1:klev,1),'trp1','trp1') |
---|
1353 | ! @$$ CALL writeg1d(1,klev,trap(1:klev,2),'trp2','trp2') |
---|
1354 | ! @$$ CALL writeg1d(1,klev,trap(1:klev,3),'trp3','trp3') |
---|
1355 | ! @$$ CALL writeg1d(1,klev,trap(1:klev,4),'trp4','trp4') |
---|
1356 | ! @$$ CALL writeg1d(1,klev,trap(1:klev,5),'trp5','trp5') |
---|
1357 | ! @$$ CALL writeg1d(1,klev,trap(1:klev,10),'trp10','trp10') |
---|
1358 | ! @$$ CALL writeg1d(1,klev,trap(1:klev,12),'trp12','trp12') |
---|
1359 | ! @$$ CALL writeg1d(1,klev,trap(1:klev,15),'trp15','trp15') |
---|
1360 | ! @$$ CALL writeg1d(1,klev,trap(1:klev,20),'trp20','trp20') |
---|
1361 | ! @$$ CALL writeg1d(1,klev,ftra(1:klev,1),'ftr1 ','ftr1 ') |
---|
1362 | ! @$$ CALL writeg1d(1,klev,ftra(1:klev,2),'ftr2 ','ftr2 ') |
---|
1363 | ! @$$ CALL writeg1d(1,klev,ftra(1:klev,3),'ftr3 ','ftr3 ') |
---|
1364 | ! @$$ CALL writeg1d(1,klev,ftra(1:klev,4),'ftr4 ','ftr4 ') |
---|
1365 | ! @$$ CALL writeg1d(1,klev,ftra(1:klev,5),'ftr5 ','ftr5 ') |
---|
1366 | ! @$$ CALL writeg1d(1,klev,ftra(1:klev,6),'ftr6 ','ftr6 ') |
---|
1367 | ! @$$ CALL writeg1d(1,klev,ftra(1:klev,7),'ftr7 ','ftr7 ') |
---|
1368 | ! @$$ CALL writeg1d(1,klev,ftra(1:klev,8),'ftr8 ','ftr8 ') |
---|
1369 | ! @$$ CALL writeg1d(1,klev,ftra(1:klev,9),'ftr9 ','ftr9 ') |
---|
1370 | ! @$$ CALL writeg1d(1,klev,ftra(1:klev,10),'ftr10','ftr10') |
---|
1371 | ! @$$ CALL writeg1d(1,klev,ftra(1:klev,11),'ftr11','ftr11') |
---|
1372 | ! @$$ CALL writeg1d(1,klev,ftra(1:klev,12),'ftr12','ftr12') |
---|
1373 | ! @$$ CALL writeg1d(1,klev,ftra(1:klev,13),'ftr13','ftr13') |
---|
1374 | ! @$$ CALL writeg1d(1,klev,ftra(1:klev,14),'ftr14','ftr14') |
---|
1375 | ! @$$ CALL writeg1d(1,klev,ftra(1:klev,15),'ftr15','ftr15') |
---|
1376 | ! @$$ CALL writeg1d(1,klev,ftra(1:klev,16),'ftr16','ftr16') |
---|
1377 | ! @$$ CALL writeg1d(1,klev,ftra(1:klev,17),'ftr17','ftr17') |
---|
1378 | ! @$$ CALL writeg1d(1,klev,ftra(1:klev,18),'ftr18','ftr18') |
---|
1379 | ! @$$ CALL writeg1d(1,klev,ftra(1:klev,19),'ftr19','ftr19') |
---|
1380 | ! @$$ CALL writeg1d(1,klev,ftra(1:klev,20),'ftr20','ftr20 ') |
---|
1381 | ! @$$ CALL writeg1d(1,klev,mp,'mp ','mp ') |
---|
1382 | ! @$$ CALL writeg1d(1,klev,Mke,'Mke ','Mke ') |
---|
1383 | |
---|
1384 | |
---|
1385 | |
---|
1386 | ! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
---|
1387 | |
---|
1388 | END SUBROUTINE convect3 |
---|
1389 | ! --------------------------------------------------------------------------- |
---|