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

Last change on this file since 2188 was 2110, checked in by lguez, 10 years ago

Abort if surface pressure becomes negative. The call to abort_gcm
already was in the sequential version but not in dyn3dpar nor in
dyn3dmem. In dyn3dmem, there is a variable checksum_all, which is
always true. The call to MPI_ALLREDUCE which should update
checksum_all is commented out. I do not know why (performance?).

Non-ASCII characters in comments are not always rendered properly and
they risk being lost. See revision [1740].

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