GCC Code Coverage Report


Directory: ./
File: phys/convect3.f90
Date: 2022-01-11 19:19:34
Exec Total Coverage
Lines: 0 535 0.0%
Branches: 0 296 0.0%

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