source: LMDZ5/trunk/libf/phylmd/convect2.F90 @ 3895

Last change on this file since 3895 was 2346, checked in by Ehouarn Millour, 9 years ago

Physics/dynamics separation:

  • remove all references to dimensions.h from physics. nbp_lon (==iim) , nbp_lat (==jjm+1) and nbp_lev (==llm) from mod_grid_phy_lmdz should be used instead.
  • added module regular_lonlat_mod in phy_common to store information about the global (lon-lat) grid cell boundaries and centers.

EM
  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 41.5 KB
RevLine 
[1992]1
[1403]2! $Id: convect2.F90 2346 2015-08-21 15:13:46Z emillour $
[524]3
[1992]4SUBROUTINE convect2(ncum, idcum, len, nd, ndp1, nl, minorig, nk1, icb1, t1, &
5    q1, qs1, u1, v1, gz1, tv1, tp1, tvp1, clw1, h1, lv1, cpn1, p1, ph1, ft1, &
6    fq1, fu1, fv1, tnk1, qnk1, gznk1, plcl1, precip1, cbmf1, iflag1, delt, &
7    cpd, cpv, cl, rv, rd, lv0, g, sigs, sigd, elcrit, tlcrit, omtsnow, dtmax, &
8    damp, alpha, entp, coeffs, coeffr, omtrain, cu, ma)
9  ! .............................START PROLOGUE............................
10
11  ! SCCS IDENTIFICATION:  @(#)convect2.f        1.2 05/18/00
12  ! 22:06:22 /h/cm/library/nogaps4/src/sub/fcst/convect2.f_v
13
14  ! CONFIGURATION IDENTIFICATION:  None
15
16  ! MODULE NAME:  convect2
17
18  ! DESCRIPTION:
19
20  ! convect1     The Emanuel Cumulus Convection Scheme - compute tendencies
21
22  ! CONTRACT NUMBER AND TITLE:  None
23
24  ! REFERENCES: Programmers  K. Emanuel (MIT), Timothy F. Hogan, M. Peng
25  ! (NRL)
26
27  ! CLASSIFICATION:  Unclassified
28
29  ! RESTRICTIONS: None
30
31  ! COMPILER DEPENDENCIES: FORTRAN 77, FORTRAN 90
32
33  ! COMPILE OPTIONS: Fortran 77: -Zu -Wf"-ei -o aggress"
34  ! Fortran 90: -O vector3,scalar3,task1,aggress,overindex  -ei -r 2
35
36  ! LIBRARIES OF RESIDENCE: /a/ops/lib/libfcst159.a
37
38  ! USAGE: call convect2(ncum,idcum,len,nd,nl,minorig,
39  ! &                 nk1,icb1,
40  ! &                 t1,q1,qs1,u1,v1,gz1,tv1,tp1,tvp1,clw1,h1,
41  ! &                 lv1,cpn1,p1,ph1,ft1,fq1,fu1,fv1,
42  ! &                 tnk1,qnk1,gznk1,plcl1,
43  ! &                 precip1,cbmf1,iflag1,
44  ! &                 delt,cpd,cpv,cl,rv,rd,lv0,g,
45  ! &                 sigs,sigd,elcrit,tlcrit,omtsnow,dtmax,damp,
46  ! &                 alpha,entp,coeffs,coeffr,omtrain,cu)
47
48  ! PARAMETERS:
49  ! Name            Type         Usage            Description
50  ! ----------      ----------     -------  ----------------------------
51
52  ! ncum          Integer        Input        number of cumulus points
53  ! idcum         Integer        Input        index of cumulus point
54  ! len           Integer        Input        first dimension
55  ! nd            Integer        Input        total vertical dimension
56  ! ndp1          Integer        Input        nd + 1
57  ! nl            Integer        Input        vertical dimension for
58  ! convection
59  ! minorig       Integer        Input        First level where convection is
60  ! allow to begin
61  ! nk1           Integer        Input        First level of convection
62  ! ncb1          Integer        Input        Level of free convection
63  ! t1            Real           Input        temperature
64  ! q1            Real           Input        specific hum
65  ! qs1           Real           Input        sat specific hum
66  ! u1            Real           Input        u-wind
67  ! v1            Real           Input        v-wind
68  ! gz1           Real           Inout        geop
69  ! tv1           Real           Input        virtual temp
70  ! tp1           Real           Input
71  ! clw1          Real           Inout        cloud liquid water
72  ! h1            Real           Inout
73  ! lv1           Real           Inout
74  ! cpn1          Real           Inout
75  ! p1            Real           Input        full level pressure
76  ! ph1           Real           Input        half level pressure
77  ! ft1           Real           Output       temp tend
78  ! fq1           Real           Output       spec hum tend
79  ! fu1           Real           Output       u-wind tend
80  ! fv1           Real           Output       v-wind tend
81  ! precip1       Real           Output       prec
82  ! cbmf1         Real           In/Out       cumulus mass flux
83  ! iflag1        Integer        Output       iflag on latitude strip
84  ! delt          Real           Input        time step
85  ! cpd           Integer        Input        See description below
86  ! cpv           Integer        Input        See description below
87  ! cl            Integer        Input        See description below
88  ! rv            Integer        Input        See description below
89  ! rd            Integer        Input        See description below
90  ! lv0           Integer        Input        See description below
91  ! g             Integer        Input        See description below
92  ! sigs          Integer        Input        See description below
93  ! sigd          Integer        Input        See description below
94  ! elcrit        Integer        Input        See description below
95  ! tlcrit        Integer        Input        See description below
96  ! omtsnow       Integer        Input        See description below
97  ! dtmax         Integer        Input        See description below
98  ! damp          Integer        Input        See description below
99  ! alpha         Integer        Input        See description below
100  ! ent           Integer        Input        See description below
101  ! coeffs        Integer        Input        See description below
102  ! coeffr        Integer        Input        See description below
103  ! omtrain       Integer        Input        See description below
104  ! cu            Integer        Input        See description below
105
106  ! COMMON BLOCKS:
107  ! Block      Name     Type    Usage              Notes
108  ! --------  --------   ----    ------   ------------------------
109
110  ! FILES: None
111
112  ! DATA BASES: None
113
114  ! NON-FILE INPUT/OUTPUT: None
115
116  ! ERROR CONDITIONS: None
117
118  ! ADDITIONAL COMMENTS: None
119
120  ! .................MAINTENANCE SECTION................................
121
122  ! MODULES CALLED:
123  ! Name           Description
124  ! zilch        Zero out an array
125  ! -------     ----------------------
126  ! LOCAL VARIABLES AND
127  ! STRUCTURES:
128  ! Name     Type    Description
129  ! -------  ------  -----------
130  ! See Comments Below
131
132  ! i        Integer loop index
133  ! k        Integer loop index
134
135  ! METHOD:
136
137  ! See Emanuel, K. and M. Zivkovic-Rothman, 2000: Development and evaluation
138  ! of a
139  ! convective scheme for use in climate models.
140
141  ! FILES: None
142
143  ! INCLUDE FILES: None
144
145  ! MAKEFILE: /a/ops/met/nogaps/src/sub/fcst/fcst159lib.mak
146
147  ! ..............................END PROLOGUE.............................
148
149
150  USE dimphy
151  IMPLICIT NONE
152
153  INTEGER kmax2, imax2, kmin2, imin2
154  REAL ftmax2, ftmin2
155  INTEGER kmax, imax, kmin, imin
156  REAL ftmax, ftmin
157
158  INTEGER ncum
[2110]159  INTEGER len
[1992]160  INTEGER idcum(len)
161  INTEGER nd
162  INTEGER ndp1
163  INTEGER nl
164  INTEGER minorig
165  INTEGER nk1(len)
166  INTEGER icb1(len)
167  REAL t1(len, nd)
168  REAL q1(len, nd)
169  REAL qs1(len, nd)
170  REAL u1(len, nd)
171  REAL v1(len, nd)
172  REAL gz1(len, nd)
173  REAL tv1(len, nd)
174  REAL tp1(len, nd)
175  REAL tvp1(len, nd)
176  REAL clw1(len, nd)
177  REAL h1(len, nd)
178  REAL lv1(len, nd)
179  REAL cpn1(len, nd)
180  REAL p1(len, nd)
181  REAL ph1(len, ndp1)
182  REAL ft1(len, nd)
183  REAL fq1(len, nd)
184  REAL fu1(len, nd)
185  REAL fv1(len, nd)
186  REAL tnk1(len)
187  REAL qnk1(len)
188  REAL gznk1(len)
189  REAL precip1(len)
190  REAL cbmf1(len)
191  REAL plcl1(len)
192  INTEGER iflag1(len)
193  REAL delt
194  REAL cpd
195  REAL cpv
196  REAL cl
197  REAL rv
198  REAL rd
199  REAL lv0
200  REAL g
201  REAL sigs ! SIGS IS THE FRACTION OF PRECIPITATION FALLING OUTSIDE
202  REAL sigd ! SIGD IS THE FRACTIONAL AREA COVERED BY UNSATURATED DNDRAFT
203  REAL elcrit ! ELCRIT IS THE AUTOCONVERSION THERSHOLD WATER CONTENT (gm/gm)
204  REAL tlcrit ! TLCRIT IS CRITICAL TEMPERATURE BELOW WHICH THE AUTO-
205  ! CONVERSION THRESHOLD IS ASSUMED TO BE ZERO
206  REAL omtsnow ! OMTSNOW IS THE ASSUMED FALL SPEED (P/s) OF SNOW
207  REAL dtmax ! DTMAX IS THE MAXIMUM NEGATIVE TEMPERATURE PERTURBATION
208  ! A LIFTED PARCEL IS ALLOWED TO HAVE BELOW ITS LFC.
209  REAL damp
210  REAL alpha
211  REAL entp ! ENTP IS THE COEFFICIENT OF MIXING IN THE ENTRAINMENT FORMULATION
212  REAL coeffs ! COEFFS IS A COEFFICIENT GOVERNING THE RATE OF EVAPORATION OF
213  ! SNOW
214  REAL coeffr ! COEFFR IS A COEFFICIENT GOVERNING THE RATE OF EVAPORATION OF
215  ! RAIN
216  REAL omtrain ! OMTRAIN IS THE ASSUMED FALL SPEED (P/s) OF RAIN
217  REAL cu ! CU IS THE COEFFICIENT GOVERNING CONVECTIVE MOMENTUM TRANSPORT
218
219  REAL ma(len, nd)
220
221
222  ! *** ELCRIT IS THE AUTOCONVERSION THERSHOLD WATER CONTENT (gm/gm) ***
223  ! ***  TLCRIT IS CRITICAL TEMPERATURE BELOW WHICH THE AUTO-        ***
224  ! ***       CONVERSION THRESHOLD IS ASSUMED TO BE ZERO             ***
225  ! ***     (THE AUTOCONVERSION THRESHOLD VARIES LINEARLY            ***
226  ! ***               BETWEEN 0 C AND TLCRIT)                        ***
227  ! ***   ENTP IS THE COEFFICIENT OF MIXING IN THE ENTRAINMENT       ***
228  ! ***                       FORMULATION                            ***
229  ! ***  SIGD IS THE FRACTIONAL AREA COVERED BY UNSATURATED DNDRAFT  ***
230  ! ***  SIGS IS THE FRACTION OF PRECIPITATION FALLING OUTSIDE       ***
231  ! ***                        OF CLOUD                              ***
232  ! ***        OMTRAIN IS THE ASSUMED FALL SPEED (P/s) OF RAIN       ***
233  ! ***     OMTSNOW IS THE ASSUMED FALL SPEED (P/s) OF SNOW          ***
234  ! ***  COEFFR IS A COEFFICIENT GOVERNING THE RATE OF EVAPORATION   ***
235  ! ***                          OF RAIN                             ***
236  ! ***  COEFFS IS A COEFFICIENT GOVERNING THE RATE OF EVAPORATION   ***
237  ! ***                          OF SNOW                             ***
238  ! ***     CU IS THE COEFFICIENT GOVERNING CONVECTIVE MOMENTUM      ***
239  ! ***                         TRANSPORT                            ***
240  ! ***    DTMAX IS THE MAXIMUM NEGATIVE TEMPERATURE PERTURBATION    ***
241  ! ***        A LIFTED PARCEL IS ALLOWED TO HAVE BELOW ITS LFC      ***
242  ! ***    ALPHA AND DAMP ARE PARAMETERS THAT CONTROL THE RATE OF    ***
243  ! ***                 APPROACH TO QUASI-EQUILIBRIUM                ***
244  ! ***   (THEIR STANDARD VALUES ARE  0.20 AND 0.1, RESPECTIVELY)    ***
245  ! ***                   (DAMP MUST BE LESS THAN 1)                 ***
246
247  ! Local arrays.
248
249  REAL work(ncum)
250  REAL t(ncum, klev)
251  REAL q(ncum, klev)
252  REAL qs(ncum, klev)
253  REAL u(ncum, klev)
254  REAL v(ncum, klev)
255  REAL gz(ncum, klev)
256  REAL h(ncum, klev)
257  REAL lv(ncum, klev)
258  REAL cpn(ncum, klev)
259  REAL p(ncum, klev)
260  REAL ph(ncum, klev)
261  REAL ft(ncum, klev)
262  REAL fq(ncum, klev)
263  REAL fu(ncum, klev)
264  REAL fv(ncum, klev)
265  REAL precip(ncum)
266  REAL cbmf(ncum)
267  REAL plcl(ncum)
268  REAL tnk(ncum)
269  REAL qnk(ncum)
270  REAL gznk(ncum)
271  REAL tv(ncum, klev)
272  REAL tp(ncum, klev)
273  REAL tvp(ncum, klev)
274  REAL clw(ncum, klev)
275  ! real det(ncum,klev)
276  REAL dph(ncum, klev)
277  ! real wd(ncum)
278  ! real tprime(ncum)
279  ! real qprime(ncum)
280  REAL ah0(ncum)
281  REAL ep(ncum, klev)
282  REAL sigp(ncum, klev)
283  INTEGER nent(ncum, klev)
284  REAL water(ncum, klev)
285  REAL evap(ncum, klev)
286  REAL mp(ncum, klev)
287  REAL m(ncum, klev)
288  REAL qti
289  REAL wt(ncum, klev)
290  REAL hp(ncum, klev)
291  REAL lvcp(ncum, klev)
292  REAL elij(ncum, klev, klev)
293  REAL ment(ncum, klev, klev)
294  REAL sij(ncum, klev, klev)
295  REAL qent(ncum, klev, klev)
296  REAL uent(ncum, klev, klev)
297  REAL vent(ncum, klev, klev)
298  REAL qp(ncum, klev)
299  REAL up(ncum, klev)
300  REAL vp(ncum, klev)
301  REAL cape(ncum)
302  REAL capem(ncum)
303  REAL frac(ncum)
304  REAL dtpbl(ncum)
305  REAL tvpplcl(ncum)
306  REAL tvaplcl(ncum)
307  REAL dtmin(ncum)
308  REAL w3d(ncum, klev)
309  REAL am(ncum)
310  REAL ents(ncum)
311  REAL uav(ncum)
312  REAL vav(ncum)
313
314  INTEGER iflag(ncum)
315  INTEGER nk(ncum)
316  INTEGER icb(ncum)
317  INTEGER inb(ncum)
318  INTEGER inb1(ncum)
319  INTEGER jtt(ncum)
320
321  INTEGER nn, i, k, n, icbmax, nlp, j
322  INTEGER ij
323  INTEGER nn2, nn3
324  REAL clmcpv
325  REAL clmcpd
326  REAL cpdmcp
327  REAL cpvmcpd
328  REAL eps
329  REAL epsi
330  REAL epsim1
331  REAL tg, qg, s, alv, tc, ahg, denom, es, rg, ginv, rowl
332  REAL delti
333  REAL tca, elacrit
334  REAL by, defrac
335  ! real byp
336  REAL byp(ncum)
337  LOGICAL lcape(ncum)
338  REAL dbo
339  REAL bf2, anum, dei, altem, cwat, stemp
340  REAL alt, qp1, smid, sjmax, sjmin
341  REAL delp, delm
342  REAL awat, coeff, afac, revap, dhdp, fac, qstm, rat
343  REAL qsm, sigt, b6, c6
344  REAL dpinv, cpinv
345  REAL fqold, ftold, fuold, fvold
346  REAL wdtrain(ncum), xxx
347  REAL bsum(ncum, klev)
348  REAL asij(ncum)
349  REAL smin(ncum)
350  REAL scrit(ncum)
351  ! real amp1,ad
352  REAL amp1(ncum), ad(ncum)
353  LOGICAL lwork(ncum)
354  INTEGER num1, num2
355
356  ! print*,'cpd en entree de convect2 ',cpd
357  nlp = nl + 1
358
359  rowl = 1000.0
360  ginv = 1.0/g
361  delti = 1.0/delt
362
363  ! Define some thermodynamic variables.
364
365  clmcpv = cl - cpv
366  clmcpd = cl - cpd
367  cpdmcp = cpd - cpv
368  cpvmcpd = cpv - cpd
369  eps = rd/rv
370  epsi = 1.0/eps
371  epsim1 = epsi - 1.0
372
373  ! Compress the fields.
374
375  DO k = 1, nl + 1
376    nn = 0
377    DO i = 1, len
378      IF (iflag1(i)==0) THEN
379        nn = nn + 1
380        t(nn, k) = t1(i, k)
381        q(nn, k) = q1(i, k)
382        qs(nn, k) = qs1(i, k)
383        u(nn, k) = u1(i, k)
384        v(nn, k) = v1(i, k)
385        gz(nn, k) = gz1(i, k)
386        h(nn, k) = h1(i, k)
387        lv(nn, k) = lv1(i, k)
388        cpn(nn, k) = cpn1(i, k)
389        p(nn, k) = p1(i, k)
390        ph(nn, k) = ph1(i, k)
391        tv(nn, k) = tv1(i, k)
392        tp(nn, k) = tp1(i, k)
393        tvp(nn, k) = tvp1(i, k)
394        clw(nn, k) = clw1(i, k)
395      END IF
396    END DO
397    ! print*,'100 ncum,nn',ncum,nn
398  END DO
399  nn = 0
400  DO i = 1, len
401    IF (iflag1(i)==0) THEN
402      nn = nn + 1
403      cbmf(nn) = cbmf1(i)
404      plcl(nn) = plcl1(i)
405      tnk(nn) = tnk1(i)
406      qnk(nn) = qnk1(i)
407      gznk(nn) = gznk1(i)
408      nk(nn) = nk1(i)
409      icb(nn) = icb1(i)
410      iflag(nn) = iflag1(i)
411    END IF
412  END DO
413  ! print*,'150 ncum,nn',ncum,nn
414
415  ! Initialize the tendencies, det, wd, tprime, qprime.
416
417  DO k = 1, nl
418    DO i = 1, ncum
419      ! det(i,k)=0.0
420      ft(i, k) = 0.0
421      fu(i, k) = 0.0
422      fv(i, k) = 0.0
423      fq(i, k) = 0.0
424      dph(i, k) = ph(i, k) - ph(i, k+1)
425      ep(i, k) = 0.0
426      sigp(i, k) = sigs
427    END DO
428  END DO
429  DO i = 1, ncum
430    ! wd(i)=0.0
431    ! tprime(i)=0.0
432    ! qprime(i)=0.0
433    precip(i) = 0.0
434    ft(i, nl+1) = 0.0
435    fu(i, nl+1) = 0.0
436    fv(i, nl+1) = 0.0
437    fq(i, nl+1) = 0.0
438  END DO
439
440  ! Compute icbmax.
441
442  icbmax = 2
443  DO i = 1, ncum
444    icbmax = max(icbmax, icb(i))
445  END DO
446
447
448  ! =====================================================================
449  ! --- FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
450  ! =====================================================================
451
452  ! ---       The procedure is to solve the equation.
453  ! cp*tp+L*qp+phi=cp*tnk+L*qnk+gznk.
454
455  ! ***  Calculate certain parcel quantities, including static energy   ***
456
457
458  DO i = 1, ncum
459    ah0(i) = (cpd*(1.-qnk(i))+cl*qnk(i))*tnk(i) + qnk(i)*(lv0-clmcpv*(tnk(i)- &
460      273.15)) + gznk(i)
461  END DO
462
463
464  ! ***  Find lifted parcel quantities above cloud base    ***
465
466
467  DO k = minorig + 1, nl
468    DO i = 1, ncum
469      IF (k>=(icb(i)+1)) THEN
470        tg = t(i, k)
471        qg = qs(i, k)
472        alv = lv0 - clmcpv*(t(i,k)-273.15)
473
474        ! First iteration.
475
476        s = cpd + alv*alv*qg/(rv*t(i,k)*t(i,k))
477        s = 1./s
478        ahg = cpd*tg + (cl-cpd)*qnk(i)*t(i, k) + alv*qg + gz(i, k)
479        tg = tg + s*(ah0(i)-ahg)
480        tg = max(tg, 35.0)
481        tc = tg - 273.15
482        denom = 243.5 + tc
483        IF (tc>=0.0) THEN
484          es = 6.112*exp(17.67*tc/denom)
485        ELSE
486          es = exp(23.33086-6111.72784/tg+0.15215*log(tg))
487        END IF
488        qg = eps*es/(p(i,k)-es*(1.-eps))
489
490        ! Second iteration.
491
492        s = cpd + alv*alv*qg/(rv*t(i,k)*t(i,k))
493        s = 1./s
494        ahg = cpd*tg + (cl-cpd)*qnk(i)*t(i, k) + alv*qg + gz(i, k)
495        tg = tg + s*(ah0(i)-ahg)
496        tg = max(tg, 35.0)
497        tc = tg - 273.15
498        denom = 243.5 + tc
499        IF (tc>=0.0) THEN
500          es = 6.112*exp(17.67*tc/denom)
501        ELSE
502          es = exp(23.33086-6111.72784/tg+0.15215*log(tg))
503        END IF
504        qg = eps*es/(p(i,k)-es*(1.-eps))
505
506        alv = lv0 - clmcpv*(t(i,k)-273.15)
507        ! print*,'cpd dans convect2 ',cpd
508        ! print*,'tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd'
509        ! print*,tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd
510        tp(i, k) = (ah0(i)-(cl-cpd)*qnk(i)*t(i,k)-gz(i,k)-alv*qg)/cpd
511        ! if (.not.cpd.gt.1000.) then
512        ! print*,'CPD=',cpd
513        ! stop
514        ! endif
515        clw(i, k) = qnk(i) - qg
516        clw(i, k) = max(0.0, clw(i,k))
517        rg = qg/(1.-qnk(i))
518        tvp(i, k) = tp(i, k)*(1.+rg*epsi)
519      END IF
520    END DO
521  END DO
522
523  ! =====================================================================
524  ! --- SET THE PRECIPITATION EFFICIENCIES AND THE FRACTION OF
525  ! --- PRECIPITATION FALLING OUTSIDE OF CLOUD
526  ! --- THESE MAY BE FUNCTIONS OF TP(I), P(I) AND CLW(I)
527  ! =====================================================================
528
529  DO k = minorig + 1, nl
530    DO i = 1, ncum
531      IF (k>=(nk(i)+1)) THEN
532        tca = tp(i, k) - 273.15
533        IF (tca>=0.0) THEN
534          elacrit = elcrit
535        ELSE
536          elacrit = elcrit*(1.0-tca/tlcrit)
537        END IF
538        elacrit = max(elacrit, 0.0)
539        ep(i, k) = 1.0 - elacrit/max(clw(i,k), 1.0E-8)
540        ep(i, k) = max(ep(i,k), 0.0)
541        ep(i, k) = min(ep(i,k), 1.0)
542        sigp(i, k) = sigs
543      END IF
544    END DO
545  END DO
546
547  ! =====================================================================
548  ! --- CALCULATE VIRTUAL TEMPERATURE AND LIFTED PARCEL
549  ! --- VIRTUAL TEMPERATURE
550  ! =====================================================================
551
552  DO k = minorig + 1, nl
553    DO i = 1, ncum
554      IF (k>=(icb(i)+1)) THEN
555        tvp(i, k) = tvp(i, k)*(1.0-qnk(i)+ep(i,k)*clw(i,k))
556        ! print*,'i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)'
557        ! print*, i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)
558      END IF
559    END DO
560  END DO
561  DO i = 1, ncum
562    tvp(i, nlp) = tvp(i, nl) - (gz(i,nlp)-gz(i,nl))/cpd
563  END DO
564
565
566  ! =====================================================================
567  ! --- NOW INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS
568  ! =====================================================================
569
570  DO i = 1, ncum*nlp
571    nent(i, 1) = 0
572    water(i, 1) = 0.0
573    evap(i, 1) = 0.0
574    mp(i, 1) = 0.0
575    m(i, 1) = 0.0
576    wt(i, 1) = omtsnow
577    hp(i, 1) = h(i, 1)
578    ! if(.not.cpn(i,1).gt.900.) then
579    ! print*,'i,lv(i,1),cpn(i,1)'
580    ! print*, i,lv(i,1),cpn(i,1)
581    ! k=(i-1)/ncum+1
582    ! print*,'i,k',mod(i,ncum),k,'  cpn',cpn(mod(i,ncum),k)
583    ! stop
584    ! endif
585    lvcp(i, 1) = lv(i, 1)/cpn(i, 1)
586  END DO
587
588  DO i = 1, ncum*nlp*nlp
589    elij(i, 1, 1) = 0.0
590    ment(i, 1, 1) = 0.0
591    sij(i, 1, 1) = 0.0
592  END DO
593
594  DO k = 1, nlp
595    DO j = 1, nlp
596      DO i = 1, ncum
597        qent(i, k, j) = q(i, j)
598        uent(i, k, j) = u(i, j)
599        vent(i, k, j) = v(i, j)
600      END DO
601    END DO
602  END DO
603
604  DO i = 1, ncum
605    qp(i, 1) = q(i, 1)
606    up(i, 1) = u(i, 1)
607    vp(i, 1) = v(i, 1)
608  END DO
609  DO k = 2, nlp
610    DO i = 1, ncum
611      qp(i, k) = q(i, k-1)
612      up(i, k) = u(i, k-1)
613      vp(i, k) = v(i, k-1)
614    END DO
615  END DO
616
617  ! =====================================================================
618  ! --- FIND THE FIRST MODEL LEVEL (INB1) ABOVE THE PARCEL'S
619  ! --- HIGHEST LEVEL OF NEUTRAL BUOYANCY
620  ! --- AND THE HIGHEST LEVEL OF POSITIVE CAPE (INB)
621  ! =====================================================================
622
623  DO i = 1, ncum
624    cape(i) = 0.0
625    capem(i) = 0.0
626    inb(i) = icb(i) + 1
627    inb1(i) = inb(i)
628  END DO
629
630  ! Originial Code
631
632  ! do 530 k=minorig+1,nl-1
633  ! do 520 i=1,ncum
634  ! if(k.ge.(icb(i)+1))then
635  ! by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
636  ! byp=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
637  ! cape(i)=cape(i)+by
638  ! if(by.ge.0.0)inb1(i)=k+1
639  ! if(cape(i).gt.0.0)then
640  ! inb(i)=k+1
641  ! capem(i)=cape(i)
642  ! endif
643  ! endif
644  ! 520    continue
645  ! 530  continue
646  ! do 540 i=1,ncum
647  ! byp=(tvp(i,nl)-tv(i,nl))*dph(i,nl)/p(i,nl)
648  ! cape(i)=capem(i)+byp
649  ! defrac=capem(i)-cape(i)
650  ! defrac=max(defrac,0.001)
651  ! frac(i)=-cape(i)/defrac
652  ! frac(i)=min(frac(i),1.0)
653  ! frac(i)=max(frac(i),0.0)
654  ! 540   continue
655
656  ! K Emanuel fix
657
658  ! call zilch(byp,ncum)
659  ! do 530 k=minorig+1,nl-1
660  ! do 520 i=1,ncum
661  ! if(k.ge.(icb(i)+1))then
662  ! by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
663  ! cape(i)=cape(i)+by
664  ! if(by.ge.0.0)inb1(i)=k+1
665  ! if(cape(i).gt.0.0)then
666  ! inb(i)=k+1
667  ! capem(i)=cape(i)
668  ! byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
669  ! endif
670  ! endif
671  ! 520    continue
672  ! 530  continue
673  ! do 540 i=1,ncum
674  ! inb(i)=max(inb(i),inb1(i))
675  ! cape(i)=capem(i)+byp(i)
676  ! defrac=capem(i)-cape(i)
677  ! defrac=max(defrac,0.001)
678  ! frac(i)=-cape(i)/defrac
679  ! frac(i)=min(frac(i),1.0)
680  ! frac(i)=max(frac(i),0.0)
681  ! 540   continue
682
683  ! J Teixeira fix
684
685  CALL zilch(byp, ncum)
686  DO i = 1, ncum
687    lcape(i) = .TRUE.
688  END DO
689  DO k = minorig + 1, nl - 1
690    DO i = 1, ncum
691      IF (cape(i)<0.0) lcape(i) = .FALSE.
692      IF ((k>=(icb(i)+1)) .AND. lcape(i)) THEN
693        by = (tvp(i,k)-tv(i,k))*dph(i, k)/p(i, k)
694        byp(i) = (tvp(i,k+1)-tv(i,k+1))*dph(i, k+1)/p(i, k+1)
695        cape(i) = cape(i) + by
696        IF (by>=0.0) inb1(i) = k + 1
697        IF (cape(i)>0.0) THEN
698          inb(i) = k + 1
699          capem(i) = cape(i)
700        END IF
701      END IF
702    END DO
703  END DO
704  DO i = 1, ncum
705    cape(i) = capem(i) + byp(i)
706    defrac = capem(i) - cape(i)
707    defrac = max(defrac, 0.001)
708    frac(i) = -cape(i)/defrac
709    frac(i) = min(frac(i), 1.0)
710    frac(i) = max(frac(i), 0.0)
711  END DO
712
713  ! =====================================================================
714  ! ---   CALCULATE LIQUID WATER STATIC ENERGY OF LIFTED PARCEL
715  ! =====================================================================
716
717  DO k = minorig + 1, nl
718    DO i = 1, ncum
719      IF ((k>=icb(i)) .AND. (k<=inb(i))) THEN
720        hp(i, k) = h(i, nk(i)) + (lv(i,k)+(cpd-cpv)*t(i,k))*ep(i, k)*clw(i, k &
721          )
722      END IF
723    END DO
724  END DO
725
726  ! =====================================================================
727  ! ---  CALCULATE CLOUD BASE MASS FLUX AND RATES OF MIXING, M(I),
728  ! --- AT EACH MODEL LEVEL
729  ! =====================================================================
730
731  ! tvpplcl = parcel temperature lifted adiabatically from level
732  ! icb-1 to the LCL.
733  ! tvaplcl = virtual temperature at the LCL.
734
735  DO i = 1, ncum
736    dtpbl(i) = 0.0
737    tvpplcl(i) = tvp(i, icb(i)-1) - rd*tvp(i, icb(i)-1)*(p(i,icb(i)-1)-plcl(i &
738      ))/(cpn(i,icb(i)-1)*p(i,icb(i)-1))
739    tvaplcl(i) = tv(i, icb(i)) + (tvp(i,icb(i))-tvp(i,icb(i)+1))*(plcl(i)-p(i &
740      ,icb(i)))/(p(i,icb(i))-p(i,icb(i)+1))
741  END DO
742
743  ! -------------------------------------------------------------------
744  ! --- Interpolate difference between lifted parcel and
745  ! --- environmental temperatures to lifted condensation level
746  ! -------------------------------------------------------------------
747
748  ! dtpbl = average of tvp-tv in the PBL (k=nk to icb-1).
749
750  DO k = minorig, icbmax
751    DO i = 1, ncum
752      IF ((k>=nk(i)) .AND. (k<=(icb(i)-1))) THEN
753        dtpbl(i) = dtpbl(i) + (tvp(i,k)-tv(i,k))*dph(i, k)
754      END IF
755    END DO
756  END DO
757  DO i = 1, ncum
758    dtpbl(i) = dtpbl(i)/(ph(i,nk(i))-ph(i,icb(i)))
759    dtmin(i) = tvpplcl(i) - tvaplcl(i) + dtmax + dtpbl(i)
760  END DO
761
762  ! -------------------------------------------------------------------
763  ! --- Adjust cloud base mass flux
764  ! -------------------------------------------------------------------
765
766  DO i = 1, ncum
767    work(i) = cbmf(i)
768    cbmf(i) = max(0.0, (1.0-damp)*cbmf(i)+0.1*alpha*dtmin(i))
769    IF ((work(i)==0.0) .AND. (cbmf(i)==0.0)) THEN
770      iflag(i) = 3
771    END IF
772  END DO
773
774  ! -------------------------------------------------------------------
775  ! --- Calculate rates of mixing,  m(i)
776  ! -------------------------------------------------------------------
777
778  CALL zilch(work, ncum)
779
780  DO j = minorig + 1, nl
781    DO i = 1, ncum
782      IF ((j>=(icb(i)+1)) .AND. (j<=inb(i))) THEN
783        k = min(j, inb1(i))
784        dbo = abs(tv(i,k+1)-tvp(i,k+1)-tv(i,k-1)+tvp(i,k-1)) + &
785          entp*0.04*(ph(i,k)-ph(i,k+1))
786        work(i) = work(i) + dbo
787        m(i, j) = cbmf(i)*dbo
788      END IF
789    END DO
790  END DO
791  DO k = minorig + 1, nl
792    DO i = 1, ncum
793      IF ((k>=(icb(i)+1)) .AND. (k<=inb(i))) THEN
794        m(i, k) = m(i, k)/work(i)
795      END IF
796    END DO
797  END DO
798
799
800  ! =====================================================================
801  ! --- CALCULATE ENTRAINED AIR MASS FLUX (ment), TOTAL WATER MIXING
802  ! --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
803  ! --- FRACTION (sij)
804  ! =====================================================================
805
806
807  DO i = minorig + 1, nl
808    DO j = minorig + 1, nl
809      DO ij = 1, ncum
810        IF ((i>=(icb(ij)+1)) .AND. (j>=icb(ij)) .AND. (i<=inb(ij)) .AND. (j<= &
811            inb(ij))) THEN
812          qti = qnk(ij) - ep(ij, i)*clw(ij, i)
813          bf2 = 1. + lv(ij, j)*lv(ij, j)*qs(ij, j)/(rv*t(ij,j)*t(ij,j)*cpd)
814          anum = h(ij, j) - hp(ij, i) + (cpv-cpd)*t(ij, j)*(qti-q(ij,j))
815          denom = h(ij, i) - hp(ij, i) + (cpd-cpv)*(q(ij,i)-qti)*t(ij, j)
816          dei = denom
817          IF (abs(dei)<0.01) dei = 0.01
818          sij(ij, i, j) = anum/dei
819          sij(ij, i, i) = 1.0
820          altem = sij(ij, i, j)*q(ij, i) + (1.-sij(ij,i,j))*qti - qs(ij, j)
821          altem = altem/bf2
822          cwat = clw(ij, j)*(1.-ep(ij,j))
823          stemp = sij(ij, i, j)
824          IF ((stemp<0.0 .OR. stemp>1.0 .OR. altem>cwat) .AND. j>i) THEN
825            anum = anum - lv(ij, j)*(qti-qs(ij,j)-cwat*bf2)
826            denom = denom + lv(ij, j)*(q(ij,i)-qti)
827            IF (abs(denom)<0.01) denom = 0.01
828            sij(ij, i, j) = anum/denom
829            altem = sij(ij, i, j)*q(ij, i) + (1.-sij(ij,i,j))*qti - qs(ij, j)
830            altem = altem - (bf2-1.)*cwat
831          END IF
832          IF (sij(ij,i,j)>0.0 .AND. sij(ij,i,j)<0.9) THEN
833            qent(ij, i, j) = sij(ij, i, j)*q(ij, i) + (1.-sij(ij,i,j))*qti
834            uent(ij, i, j) = sij(ij, i, j)*u(ij, i) + &
835              (1.-sij(ij,i,j))*u(ij, nk(ij))
836            vent(ij, i, j) = sij(ij, i, j)*v(ij, i) + &
837              (1.-sij(ij,i,j))*v(ij, nk(ij))
838            elij(ij, i, j) = altem
839            elij(ij, i, j) = max(0.0, elij(ij,i,j))
840            ment(ij, i, j) = m(ij, i)/(1.-sij(ij,i,j))
841            nent(ij, i) = nent(ij, i) + 1
842          END IF
843          sij(ij, i, j) = max(0.0, sij(ij,i,j))
844          sij(ij, i, j) = min(1.0, sij(ij,i,j))
845        END IF
846      END DO
847    END DO
848
849    ! ***   If no air can entrain at level i assume that updraft detrains
850    ! ***
851    ! ***   at that level and calculate detrained air flux and properties
852    ! ***
853
854    DO ij = 1, ncum
855      IF ((i>=(icb(ij)+1)) .AND. (i<=inb(ij)) .AND. (nent(ij,i)==0)) THEN
856        ment(ij, i, i) = m(ij, i)
857        qent(ij, i, i) = q(ij, nk(ij)) - ep(ij, i)*clw(ij, i)
858        uent(ij, i, i) = u(ij, nk(ij))
859        vent(ij, i, i) = v(ij, nk(ij))
860        elij(ij, i, i) = clw(ij, i)
861        sij(ij, i, i) = 1.0
862      END IF
863    END DO
864  END DO
865
866  DO i = 1, ncum
867    sij(i, inb(i), inb(i)) = 1.0
868  END DO
869
870  ! =====================================================================
871  ! ---  NORMALIZE ENTRAINED AIR MASS FLUXES
872  ! ---  TO REPRESENT EQUAL PROBABILITIES OF MIXING
873  ! =====================================================================
874
875
876  CALL zilch(bsum, ncum*nlp)
877  DO ij = 1, ncum
878    lwork(ij) = .FALSE.
879  END DO
880  DO i = minorig + 1, nl
881
882    num1 = 0
883    DO ij = 1, ncum
884      IF ((i>=icb(ij)+1) .AND. (i<=inb(ij))) num1 = num1 + 1
885    END DO
886    IF (num1<=0) GO TO 789
887
888    DO ij = 1, ncum
889      IF ((i>=icb(ij)+1) .AND. (i<=inb(ij))) THEN
890        lwork(ij) = (nent(ij,i)/=0)
891        qp1 = q(ij, nk(ij)) - ep(ij, i)*clw(ij, i)
892        anum = h(ij, i) - hp(ij, i) - lv(ij, i)*(qp1-qs(ij,i))
893        denom = h(ij, i) - hp(ij, i) + lv(ij, i)*(q(ij,i)-qp1)
894        IF (abs(denom)<0.01) denom = 0.01
895        scrit(ij) = anum/denom
896        alt = qp1 - qs(ij, i) + scrit(ij)*(q(ij,i)-qp1)
897        IF (scrit(ij)<0.0 .OR. alt<0.0) scrit(ij) = 1.0
898        asij(ij) = 0.0
899        smin(ij) = 1.0
900      END IF
901    END DO
902    DO j = minorig, nl
903
904      num2 = 0
905      DO ij = 1, ncum
906        IF ((i>=icb(ij)+1) .AND. (i<=inb(ij)) .AND. (j>=icb( &
907          ij)) .AND. (j<=inb(ij)) .AND. lwork(ij)) num2 = num2 + 1
908      END DO
909      IF (num2<=0) GO TO 783
910
911      DO ij = 1, ncum
912        IF ((i>=icb(ij)+1) .AND. (i<=inb(ij)) .AND. (j>=icb( &
913            ij)) .AND. (j<=inb(ij)) .AND. lwork(ij)) THEN
914          IF (sij(ij,i,j)>0.0 .AND. sij(ij,i,j)<0.9) THEN
915            IF (j>i) THEN
916              smid = min(sij(ij,i,j), scrit(ij))
917              sjmax = smid
918              sjmin = smid
919              IF (smid<smin(ij) .AND. sij(ij,i,j+1)<smid) THEN
920                smin(ij) = smid
921                sjmax = min(sij(ij,i,j+1), sij(ij,i,j), scrit(ij))
922                sjmin = max(sij(ij,i,j-1), sij(ij,i,j))
923                sjmin = min(sjmin, scrit(ij))
924              END IF
925            ELSE
926              sjmax = max(sij(ij,i,j+1), scrit(ij))
927              smid = max(sij(ij,i,j), scrit(ij))
928              sjmin = 0.0
929              IF (j>1) sjmin = sij(ij, i, j-1)
930              sjmin = max(sjmin, scrit(ij))
931            END IF
932            delp = abs(sjmax-smid)
933            delm = abs(sjmin-smid)
934            asij(ij) = asij(ij) + (delp+delm)*(ph(ij,j)-ph(ij,j+1))
935            ment(ij, i, j) = ment(ij, i, j)*(delp+delm)*(ph(ij,j)-ph(ij,j+1))
936          END IF
937        END IF
938      END DO
939783 END DO
940    DO ij = 1, ncum
941      IF ((i>=icb(ij)+1) .AND. (i<=inb(ij)) .AND. lwork(ij)) THEN
942        asij(ij) = max(1.0E-21, asij(ij))
943        asij(ij) = 1.0/asij(ij)
944        bsum(ij, i) = 0.0
945      END IF
946    END DO
947    DO j = minorig, nl + 1
948      DO ij = 1, ncum
949        IF ((i>=icb(ij)+1) .AND. (i<=inb(ij)) .AND. (j>=icb( &
950            ij)) .AND. (j<=inb(ij)) .AND. lwork(ij)) THEN
951          ment(ij, i, j) = ment(ij, i, j)*asij(ij)
952          bsum(ij, i) = bsum(ij, i) + ment(ij, i, j)
953        END IF
954      END DO
955    END DO
956    DO ij = 1, ncum
957      IF ((i>=icb(ij)+1) .AND. (i<=inb(ij)) .AND. (bsum(ij, &
958          i)<1.0E-18) .AND. lwork(ij)) THEN
959        nent(ij, i) = 0
960        ment(ij, i, i) = m(ij, i)
961        qent(ij, i, i) = q(ij, nk(ij)) - ep(ij, i)*clw(ij, i)
962        uent(ij, i, i) = u(ij, nk(ij))
963        vent(ij, i, i) = v(ij, nk(ij))
964        elij(ij, i, i) = clw(ij, i)
965        sij(ij, i, i) = 1.0
966      END IF
967    END DO
968789 END DO
969
970  ! =====================================================================
971  ! --- PRECIPITATING DOWNDRAFT CALCULATION
972  ! =====================================================================
973
974  ! ***  Check whether ep(inb)=0, if so, skip precipitating    ***
975  ! ***             downdraft calculation                      ***
976
977
978  ! ***  Integrate liquid water equation to find condensed water   ***
979  ! ***                and condensed water flux                    ***
980
981
982  DO i = 1, ncum
983    jtt(i) = 2
984    IF (ep(i,inb(i))<=0.0001) iflag(i) = 2
985    IF (iflag(i)==0) THEN
986      lwork(i) = .TRUE.
987    ELSE
988      lwork(i) = .FALSE.
989    END IF
990  END DO
991
992  ! ***                    Begin downdraft loop                    ***
993
994
995  CALL zilch(wdtrain, ncum)
996  DO i = nl + 1, 1, -1
997
998    num1 = 0
999    DO ij = 1, ncum
1000      IF ((i<=inb(ij)) .AND. lwork(ij)) num1 = num1 + 1
1001    END DO
1002    IF (num1<=0) GO TO 899
1003
1004
1005    ! ***        Calculate detrained precipitation             ***
1006
1007    DO ij = 1, ncum
1008      IF ((i<=inb(ij)) .AND. (lwork(ij))) THEN
1009        wdtrain(ij) = g*ep(ij, i)*m(ij, i)*clw(ij, i)
1010      END IF
1011    END DO
1012
1013    IF (i>1) THEN
1014      DO j = 1, i - 1
1015        DO ij = 1, ncum
1016          IF ((i<=inb(ij)) .AND. (lwork(ij))) THEN
1017            awat = elij(ij, j, i) - (1.-ep(ij,i))*clw(ij, i)
1018            awat = max(0.0, awat)
1019            wdtrain(ij) = wdtrain(ij) + g*awat*ment(ij, j, i)
1020          END IF
1021        END DO
1022      END DO
1023    END IF
1024
1025    ! ***    Find rain water and evaporation using provisional   ***
1026    ! ***              estimates of qp(i)and qp(i-1)             ***
1027
1028
1029    ! ***  Value of terminal velocity and coeffecient of evaporation for snow
1030    ! ***
1031
1032    DO ij = 1, ncum
1033      IF ((i<=inb(ij)) .AND. (lwork(ij))) THEN
1034        coeff = coeffs
1035        wt(ij, i) = omtsnow
1036
1037        ! ***  Value of terminal velocity and coeffecient of evaporation for
1038        ! rain   ***
1039
1040        IF (t(ij,i)>273.0) THEN
1041          coeff = coeffr
1042          wt(ij, i) = omtrain
1043        END IF
1044        qsm = 0.5*(q(ij,i)+qp(ij,i+1))
1045        afac = coeff*ph(ij, i)*(qs(ij,i)-qsm)/(1.0E4+2.0E3*ph(ij,i)*qs(ij,i))
1046        afac = max(afac, 0.0)
1047        sigt = sigp(ij, i)
1048        sigt = max(0.0, sigt)
1049        sigt = min(1.0, sigt)
1050        b6 = 100.*(ph(ij,i)-ph(ij,i+1))*sigt*afac/wt(ij, i)
1051        c6 = (water(ij,i+1)*wt(ij,i+1)+wdtrain(ij)/sigd)/wt(ij, i)
1052        revap = 0.5*(-b6+sqrt(b6*b6+4.*c6))
1053        evap(ij, i) = sigt*afac*revap
1054        water(ij, i) = revap*revap
1055
1056        ! ***  Calculate precipitating downdraft mass flux under     ***
1057        ! ***              hydrostatic approximation                 ***
1058
1059        IF (i>1) THEN
1060          dhdp = (h(ij,i)-h(ij,i-1))/(p(ij,i-1)-p(ij,i))
1061          dhdp = max(dhdp, 10.0)
1062          mp(ij, i) = 100.*ginv*lv(ij, i)*sigd*evap(ij, i)/dhdp
1063          mp(ij, i) = max(mp(ij,i), 0.0)
1064
1065          ! ***   Add small amount of inertia to downdraft              ***
1066
1067          fac = 20.0/(ph(ij,i-1)-ph(ij,i))
1068          mp(ij, i) = (fac*mp(ij,i+1)+mp(ij,i))/(1.+fac)
1069
1070          ! ***      Force mp to decrease linearly to zero
1071          ! ***
1072          ! ***      between about 950 mb and the surface
1073          ! ***
1074
1075          IF (p(ij,i)>(0.949*p(ij,1))) THEN
1076            jtt(ij) = max(jtt(ij), i)
1077            mp(ij, i) = mp(ij, jtt(ij))*(p(ij,1)-p(ij,i))/ &
1078              (p(ij,1)-p(ij,jtt(ij)))
1079          END IF
1080        END IF
1081
1082        ! ***       Find mixing ratio of precipitating downdraft     ***
1083
1084        IF (i/=inb(ij)) THEN
1085          IF (i==1) THEN
1086            qstm = qs(ij, 1)
1087          ELSE
1088            qstm = qs(ij, i-1)
1089          END IF
1090          IF (mp(ij,i)>mp(ij,i+1)) THEN
1091            rat = mp(ij, i+1)/mp(ij, i)
1092            qp(ij, i) = qp(ij, i+1)*rat + q(ij, i)*(1.0-rat) + &
1093              100.*ginv*sigd*(ph(ij,i)-ph(ij,i+1))*(evap(ij,i)/mp(ij,i))
1094            up(ij, i) = up(ij, i+1)*rat + u(ij, i)*(1.-rat)
1095            vp(ij, i) = vp(ij, i+1)*rat + v(ij, i)*(1.-rat)
1096          ELSE
1097            IF (mp(ij,i+1)>0.0) THEN
1098              qp(ij, i) = (gz(ij,i+1)-gz(ij,i)+qp(ij,i+1)*(lv(ij,i+1)+t(ij, &
1099                i+1)*(cl-cpd))+cpd*(t(ij,i+1)-t(ij, &
1100                i)))/(lv(ij,i)+t(ij,i)*(cl-cpd))
1101              up(ij, i) = up(ij, i+1)
1102              vp(ij, i) = vp(ij, i+1)
1103            END IF
1104          END IF
1105          qp(ij, i) = min(qp(ij,i), qstm)
1106          qp(ij, i) = max(qp(ij,i), 0.0)
1107        END IF
1108      END IF
1109    END DO
1110899 END DO
1111
1112  ! ***  Calculate surface precipitation in mm/day     ***
1113
1114  DO i = 1, ncum
1115    IF (iflag(i)<=1) THEN
1116      ! c            precip(i)=precip(i)+wt(i,1)*sigd*water(i,1)*3600.*24000.
1117      ! c     &                /(rowl*g)
1118      ! c            precip(i)=precip(i)*delt/86400.
1119      precip(i) = wt(i, 1)*sigd*water(i, 1)*86400/g
1120    END IF
1121  END DO
1122
1123
1124  ! ***  Calculate downdraft velocity scale and surface temperature and  ***
1125  ! ***                    water vapor fluctuations                      ***
1126
1127  ! wd=beta*abs(mp(icb))*0.01*rd*t(icb)/(sigd*p(icb))
1128  ! qprime=0.5*(qp(1)-q(1))
1129  ! tprime=lv0*qprime/cpd
1130
1131  ! ***  Calculate tendencies of lowest level potential temperature  ***
1132  ! ***                      and mixing ratio                        ***
1133
1134  DO i = 1, ncum
1135    work(i) = 0.01/(ph(i,1)-ph(i,2))
1136    am(i) = 0.0
1137  END DO
1138  DO k = 2, nl
1139    DO i = 1, ncum
1140      IF ((nk(i)==1) .AND. (k<=inb(i)) .AND. (nk(i)==1)) THEN
1141        am(i) = am(i) + m(i, k)
1142      END IF
1143    END DO
1144  END DO
1145  DO i = 1, ncum
1146    IF ((g*work(i)*am(i))>=delti) iflag(i) = 1
1147    ft(i, 1) = ft(i, 1) + g*work(i)*am(i)*(t(i,2)-t(i,1)+(gz(i,2)-gz(i, &
1148      1))/cpn(i,1))
1149    ft(i, 1) = ft(i, 1) - lvcp(i, 1)*sigd*evap(i, 1)
1150    ft(i, 1) = ft(i, 1) + sigd*wt(i, 2)*(cl-cpd)*water(i, 2)*(t(i,2)-t(i,1))* &
1151      work(i)/cpn(i, 1)
1152    fq(i, 1) = fq(i, 1) + g*mp(i, 2)*(qp(i,2)-q(i,1))*work(i) + &
1153      sigd*evap(i, 1)
1154    fq(i, 1) = fq(i, 1) + g*am(i)*(q(i,2)-q(i,1))*work(i)
1155    fu(i, 1) = fu(i, 1) + g*work(i)*(mp(i,2)*(up(i,2)-u(i,1))+am(i)*(u(i, &
1156      2)-u(i,1)))
1157    fv(i, 1) = fv(i, 1) + g*work(i)*(mp(i,2)*(vp(i,2)-v(i,1))+am(i)*(v(i, &
1158      2)-v(i,1)))
1159  END DO
1160  DO j = 2, nl
1161    DO i = 1, ncum
1162      IF (j<=inb(i)) THEN
1163        fq(i, 1) = fq(i, 1) + g*work(i)*ment(i, j, 1)*(qent(i,j,1)-q(i,1))
1164        fu(i, 1) = fu(i, 1) + g*work(i)*ment(i, j, 1)*(uent(i,j,1)-u(i,1))
1165        fv(i, 1) = fv(i, 1) + g*work(i)*ment(i, j, 1)*(vent(i,j,1)-v(i,1))
1166      END IF
1167    END DO
1168  END DO
1169
1170  ! ***  Calculate tendencies of potential temperature and mixing ratio  ***
1171  ! ***               at levels above the lowest level                   ***
1172
1173  ! ***  First find the net saturated updraft and downdraft mass fluxes  ***
1174  ! ***                      through each level                          ***
1175
1176  DO i = 2, nl + 1
1177
1178    num1 = 0
1179    DO ij = 1, ncum
1180      IF (i<=inb(ij)) num1 = num1 + 1
1181    END DO
1182    IF (num1<=0) GO TO 1500
1183
1184    CALL zilch(amp1, ncum)
1185    CALL zilch(ad, ncum)
1186
1187    DO k = i + 1, nl + 1
1188      DO ij = 1, ncum
1189        IF ((i>=nk(ij)) .AND. (i<=inb(ij)) .AND. (k<=(inb(ij)+1))) THEN
1190          amp1(ij) = amp1(ij) + m(ij, k)
1191        END IF
1192      END DO
1193    END DO
1194
1195    DO k = 1, i
1196      DO j = i + 1, nl + 1
1197        DO ij = 1, ncum
1198          IF ((j<=(inb(ij)+1)) .AND. (i<=inb(ij))) THEN
1199            amp1(ij) = amp1(ij) + ment(ij, k, j)
1200          END IF
1201        END DO
1202      END DO
1203    END DO
1204    DO k = 1, i - 1
1205      DO j = i, nl + 1
1206        DO ij = 1, ncum
1207          IF ((i<=inb(ij)) .AND. (j<=inb(ij))) THEN
1208            ad(ij) = ad(ij) + ment(ij, j, k)
1209          END IF
1210        END DO
1211      END DO
1212    END DO
1213
1214    DO ij = 1, ncum
1215      IF (i<=inb(ij)) THEN
1216        dpinv = 0.01/(ph(ij,i)-ph(ij,i+1))
1217        cpinv = 1.0/cpn(ij, i)
1218
1219        ft(ij, i) = ft(ij, i) + g*dpinv*(amp1(ij)*(t(ij,i+1)-t(ij, &
1220          i)+(gz(ij,i+1)-gz(ij,i))*cpinv)-ad(ij)*(t(ij,i)-t(ij, &
1221          i-1)+(gz(ij,i)-gz(ij,i-1))*cpinv)) - sigd*lvcp(ij, i)*evap(ij, i)
1222        ft(ij, i) = ft(ij, i) + g*dpinv*ment(ij, i, i)*(hp(ij,i)-h(ij,i)+t(ij &
1223          ,i)*(cpv-cpd)*(q(ij,i)-qent(ij,i,i)))*cpinv
1224        ft(ij, i) = ft(ij, i) + sigd*wt(ij, i+1)*(cl-cpd)*water(ij, i+1)*(t( &
1225          ij,i+1)-t(ij,i))*dpinv*cpinv
1226        fq(ij, i) = fq(ij, i) + g*dpinv*(amp1(ij)*(q(ij,i+1)-q(ij, &
1227          i))-ad(ij)*(q(ij,i)-q(ij,i-1)))
1228        fu(ij, i) = fu(ij, i) + g*dpinv*(amp1(ij)*(u(ij,i+1)-u(ij, &
1229          i))-ad(ij)*(u(ij,i)-u(ij,i-1)))
1230        fv(ij, i) = fv(ij, i) + g*dpinv*(amp1(ij)*(v(ij,i+1)-v(ij, &
1231          i))-ad(ij)*(v(ij,i)-v(ij,i-1)))
1232      END IF
1233    END DO
1234    DO k = 1, i - 1
1235      DO ij = 1, ncum
1236        IF (i<=inb(ij)) THEN
1237          awat = elij(ij, k, i) - (1.-ep(ij,i))*clw(ij, i)
1238          awat = max(awat, 0.0)
1239          fq(ij, i) = fq(ij, i) + g*dpinv*ment(ij, k, i)*(qent(ij,k,i)-awat-q &
1240            (ij,i))
1241          fu(ij, i) = fu(ij, i) + g*dpinv*ment(ij, k, i)*(uent(ij,k,i)-u(ij,i &
1242            ))
1243          fv(ij, i) = fv(ij, i) + g*dpinv*ment(ij, k, i)*(vent(ij,k,i)-v(ij,i &
1244            ))
1245        END IF
1246      END DO
1247    END DO
1248    DO k = i, nl + 1
1249      DO ij = 1, ncum
1250        IF ((i<=inb(ij)) .AND. (k<=inb(ij))) THEN
1251          fq(ij, i) = fq(ij, i) + g*dpinv*ment(ij, k, i)*(qent(ij,k,i)-q(ij,i &
1252            ))
1253          fu(ij, i) = fu(ij, i) + g*dpinv*ment(ij, k, i)*(uent(ij,k,i)-u(ij,i &
1254            ))
1255          fv(ij, i) = fv(ij, i) + g*dpinv*ment(ij, k, i)*(vent(ij,k,i)-v(ij,i &
1256            ))
1257        END IF
1258      END DO
1259    END DO
1260    DO ij = 1, ncum
1261      IF (i<=inb(ij)) THEN
1262        fq(ij, i) = fq(ij, i) + sigd*evap(ij, i) + g*(mp(ij,i+1)*(qp(ij, &
1263          i+1)-q(ij,i))-mp(ij,i)*(qp(ij,i)-q(ij,i-1)))*dpinv
1264        fu(ij, i) = fu(ij, i) + g*(mp(ij,i+1)*(up(ij,i+1)-u(ij, &
1265          i))-mp(ij,i)*(up(ij,i)-u(ij,i-1)))*dpinv
1266        fv(ij, i) = fv(ij, i) + g*(mp(ij,i+1)*(vp(ij,i+1)-v(ij, &
1267          i))-mp(ij,i)*(vp(ij,i)-v(ij,i-1)))*dpinv
1268      END IF
1269    END DO
12701500 END DO
1271
1272  ! *** Adjust tendencies at top of convection layer to reflect  ***
1273  ! ***       actual position of the level zero cape             ***
1274
1275  DO ij = 1, ncum
1276    fqold = fq(ij, inb(ij))
1277    fq(ij, inb(ij)) = fq(ij, inb(ij))*(1.-frac(ij))
1278    fq(ij, inb(ij)-1) = fq(ij, inb(ij)-1) + frac(ij)*fqold*((ph(ij, &
1279      inb(ij))-ph(ij,inb(ij)+1))/(ph(ij,inb(ij)-1)-ph(ij, &
1280      inb(ij))))*lv(ij, inb(ij))/lv(ij, inb(ij)-1)
1281    ftold = ft(ij, inb(ij))
1282    ft(ij, inb(ij)) = ft(ij, inb(ij))*(1.-frac(ij))
1283    ft(ij, inb(ij)-1) = ft(ij, inb(ij)-1) + frac(ij)*ftold*((ph(ij, &
1284      inb(ij))-ph(ij,inb(ij)+1))/(ph(ij,inb(ij)-1)-ph(ij, &
1285      inb(ij))))*cpn(ij, inb(ij))/cpn(ij, inb(ij)-1)
1286    fuold = fu(ij, inb(ij))
1287    fu(ij, inb(ij)) = fu(ij, inb(ij))*(1.-frac(ij))
1288    fu(ij, inb(ij)-1) = fu(ij, inb(ij)-1) + frac(ij)*fuold*((ph(ij, &
1289      inb(ij))-ph(ij,inb(ij)+1))/(ph(ij,inb(ij)-1)-ph(ij,inb(ij))))
1290    fvold = fv(ij, inb(ij))
1291    fv(ij, inb(ij)) = fv(ij, inb(ij))*(1.-frac(ij))
1292    fv(ij, inb(ij)-1) = fv(ij, inb(ij)-1) + frac(ij)*fvold*((ph(ij, &
1293      inb(ij))-ph(ij,inb(ij)+1))/(ph(ij,inb(ij)-1)-ph(ij,inb(ij))))
1294  END DO
1295
1296  ! ***   Very slightly adjust tendencies to force exact   ***
1297  ! ***     enthalpy, momentum and tracer conservation     ***
1298
1299  DO ij = 1, ncum
1300    ents(ij) = 0.0
1301    uav(ij) = 0.0
1302    vav(ij) = 0.0
1303    DO i = 1, inb(ij)
1304      ents(ij) = ents(ij) + (cpn(ij,i)*ft(ij,i)+lv(ij,i)*fq(ij,i))*(ph(ij,i)- &
1305        ph(ij,i+1))
1306      uav(ij) = uav(ij) + fu(ij, i)*(ph(ij,i)-ph(ij,i+1))
1307      vav(ij) = vav(ij) + fv(ij, i)*(ph(ij,i)-ph(ij,i+1))
1308    END DO
1309  END DO
1310  DO ij = 1, ncum
1311    ents(ij) = ents(ij)/(ph(ij,1)-ph(ij,inb(ij)+1))
1312    uav(ij) = uav(ij)/(ph(ij,1)-ph(ij,inb(ij)+1))
1313    vav(ij) = vav(ij)/(ph(ij,1)-ph(ij,inb(ij)+1))
1314  END DO
1315  DO ij = 1, ncum
1316    DO i = 1, inb(ij)
1317      ft(ij, i) = ft(ij, i) - ents(ij)/cpn(ij, i)
1318      fu(ij, i) = (1.-cu)*(fu(ij,i)-uav(ij))
1319      fv(ij, i) = (1.-cu)*(fv(ij,i)-vav(ij))
1320    END DO
1321  END DO
1322
1323  DO k = 1, nl + 1
1324    DO i = 1, ncum
1325      IF ((q(i,k)+delt*fq(i,k))<0.0) iflag(i) = 10
1326    END DO
1327  END DO
1328
1329
1330  DO i = 1, ncum
1331    IF (iflag(i)>2) THEN
1332      precip(i) = 0.0
1333      cbmf(i) = 0.0
1334    END IF
1335  END DO
1336  DO k = 1, nl
1337    DO i = 1, ncum
1338      IF (iflag(i)>2) THEN
1339        ft(i, k) = 0.0
1340        fq(i, k) = 0.0
1341        fu(i, k) = 0.0
1342        fv(i, k) = 0.0
1343      END IF
1344    END DO
1345  END DO
1346  DO i = 1, ncum
1347    precip1(idcum(i)) = precip(i)
1348    cbmf1(idcum(i)) = cbmf(i)
1349    iflag1(idcum(i)) = iflag(i)
1350  END DO
1351  DO k = 1, nl
1352    DO i = 1, ncum
1353      ft1(idcum(i), k) = ft(i, k)
1354      fq1(idcum(i), k) = fq(i, k)
1355      fu1(idcum(i), k) = fu(i, k)
1356      fv1(idcum(i), k) = fv(i, k)
1357    END DO
1358  END DO
1359
1360  DO k = 1, nd
1361    DO i = 1, len
1362      ma(i, k) = 0.
1363    END DO
1364  END DO
1365  DO k = nl, 1, -1
1366    DO i = 1, ncum
1367      ma(i, k) = ma(i, k+1) + m(i, k)
1368    END DO
1369  END DO
1370
1371  RETURN
1372END SUBROUTINE convect2
1373
Note: See TracBrowser for help on using the repository browser.