GCC Code Coverage Report
Directory: ./ Exec Total Coverage
File: phylmd/convect3.F90 Lines: 0 535 0.0 %
Date: 2023-06-30 12:56:34 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
! ---------------------------------------------------------------------------