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