source: LMDZ6/branches/contrails/libf/phylmd/cv_routines.f90 @ 5421

Last change on this file since 5421 was 5346, checked in by fhourdin, 3 weeks ago

Debut de replaysation de la convection profonde.

Regroupement de cvparam, cv3param et cvthermo (récemment
passés de statut de .h à module, dans un unique module
lmdz_cv_ini.f90
  • 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: 50.5 KB
Line 
1
2! $Id: cv_routines.f90 5346 2024-11-28 19:41:47Z fhourdin $
3
4SUBROUTINE cv_param(nd)
5  USE lmdz_cv_ini, ONLY : alpha,betad,coeffr,coeffs,cu,damp,delta,dtmax,elcrit,entp,minorig,nl,nlm,nlp,noff,omtrain,omtsnow,sigd,sigs,tlcrit
6
7  IMPLICIT NONE
8
9  ! ------------------------------------------------------------
10  ! Set parameters for convectL
11  ! (includes microphysical parameters and parameters that
12  ! control the rate of approach to quasi-equilibrium)
13  ! ------------------------------------------------------------
14
15  ! *** ELCRIT IS THE AUTOCONVERSION THERSHOLD WATER CONTENT (gm/gm) ***
16  ! ***  TLCRIT IS CRITICAL TEMPERATURE BELOW WHICH THE AUTO-        ***
17  ! ***       CONVERSION THRESHOLD IS ASSUMED TO BE ZERO             ***
18  ! ***     (THE AUTOCONVERSION THRESHOLD VARIES LINEARLY            ***
19  ! ***               BETWEEN 0 C AND TLCRIT)                        ***
20  ! ***   ENTP IS THE COEFFICIENT OF MIXING IN THE ENTRAINMENT       ***
21  ! ***                       FORMULATION                            ***
22  ! ***  SIGD IS THE FRACTIONAL AREA COVERED BY UNSATURATED DNDRAFT  ***
23  ! ***  SIGS IS THE FRACTION OF PRECIPITATION FALLING OUTSIDE       ***
24  ! ***                        OF CLOUD                              ***
25  ! ***        OMTRAIN IS THE ASSUMED FALL SPEED (P/s) OF RAIN       ***
26  ! ***     OMTSNOW IS THE ASSUMED FALL SPEED (P/s) OF SNOW          ***
27  ! ***  COEFFR IS A COEFFICIENT GOVERNING THE RATE OF EVAPORATION   ***
28  ! ***                          OF RAIN                             ***
29  ! ***  COEFFS IS A COEFFICIENT GOVERNING THE RATE OF EVAPORATION   ***
30  ! ***                          OF SNOW                             ***
31  ! ***     CU IS THE COEFFICIENT GOVERNING CONVECTIVE MOMENTUM      ***
32  ! ***                         TRANSPORT                            ***
33  ! ***    DTMAX IS THE MAXIMUM NEGATIVE TEMPERATURE PERTURBATION    ***
34  ! ***        A LIFTED PARCEL IS ALLOWED TO HAVE BELOW ITS LFC      ***
35  ! ***    ALPHA AND DAMP ARE PARAMETERS THAT CONTROL THE RATE OF    ***
36  ! ***                 APPROACH TO QUASI-EQUILIBRIUM                ***
37  ! ***   (THEIR STANDARD VALUES ARE  0.20 AND 0.1, RESPECTIVELY)    ***
38  ! ***                   (DAMP MUST BE LESS THAN 1)                 ***
39
40  INTEGER nd
41  CHARACTER (LEN=20) :: modname = 'cv_routines'
42  CHARACTER (LEN=80) :: abort_message
43
44  ! noff: integer limit for convection (nd-noff)
45  ! minorig: First level of convection
46
47  noff = 2
48  minorig = 2
49
50  nl = nd - noff
51  nlp = nl + 1
52  nlm = nl - 1
53
54  elcrit = 0.0011
55  tlcrit = -55.0
56  entp = 1.5
57  sigs = 0.12
58  sigd = 0.05
59  omtrain = 50.0
60  omtsnow = 5.5
61  coeffr = 1.0
62  coeffs = 0.8
63  dtmax = 0.9
64
65  cu = 0.70
66
67  betad = 10.0
68
69  damp = 0.1
70  alpha = 0.2
71
72  delta = 0.01 ! cld
73
74  RETURN
75END SUBROUTINE cv_param
76
77SUBROUTINE cv_prelim(len, nd, ndp1, t, q, p, ph, lv, cpn, tv, gz, h, hm)
78  USE lmdz_cv_ini, ONLY : cl,clmcpv,cpd,cpv,epsim1,hrd,lv0,nlp,t0
79
80  IMPLICIT NONE
81
82  ! =====================================================================
83  ! --- CALCULATE ARRAYS OF GEOPOTENTIAL, HEAT CAPACITY & STATIC ENERGY
84  ! =====================================================================
85
86  ! inputs:
87  INTEGER len, nd, ndp1
88  REAL t(len, nd), q(len, nd), p(len, nd), ph(len, ndp1)
89
90  ! outputs:
91  REAL lv(len, nd), cpn(len, nd), tv(len, nd)
92  REAL gz(len, nd), h(len, nd), hm(len, nd)
93
94  ! local variables:
95  INTEGER k, i
96  REAL cpx(len, nd)
97
98  DO k = 1, nlp
99    DO i = 1, len
100      lv(i, k) = lv0 - clmcpv*(t(i,k)-t0)
101      cpn(i, k) = cpd*(1.0-q(i,k)) + cpv*q(i, k)
102      cpx(i, k) = cpd*(1.0-q(i,k)) + cl*q(i, k)
103      tv(i, k) = t(i, k)*(1.0+q(i,k)*epsim1)
104    END DO
105  END DO
106
107  ! gz = phi at the full levels (same as p).
108
109  DO i = 1, len
110    gz(i, 1) = 0.0
111  END DO
112  DO k = 2, nlp
113    DO i = 1, len
114      gz(i, k) = gz(i, k-1) + hrd*(tv(i,k-1)+tv(i,k))*(p(i,k-1)-p(i,k))/ph(i, &
115        k)
116    END DO
117  END DO
118
119  ! h  = phi + cpT (dry static energy).
120  ! hm = phi + cp(T-Tbase)+Lq
121
122  DO k = 1, nlp
123    DO i = 1, len
124      h(i, k) = gz(i, k) + cpn(i, k)*t(i, k)
125      hm(i, k) = gz(i, k) + cpx(i, k)*(t(i,k)-t(i,1)) + lv(i, k)*q(i, k)
126    END DO
127  END DO
128
129  RETURN
130END SUBROUTINE cv_prelim
131
132SUBROUTINE cv_feed(len, nd, t, q, qs, p, hm, gz, nk, icb, icbmax, iflag, tnk, &
133    qnk, gznk, plcl)
134   USE lmdz_cv_ini, ONLY : minorig,nl,nlm,nlp
135
136    IMPLICIT NONE
137
138  ! ================================================================
139  ! Purpose: CONVECTIVE FEED
140  ! ================================================================
141
142
143  ! inputs:
144  INTEGER len, nd
145  REAL t(len, nd), q(len, nd), qs(len, nd), p(len, nd)
146  REAL hm(len, nd), gz(len, nd)
147
148  ! outputs:
149  INTEGER iflag(len), nk(len), icb(len), icbmax
150  REAL tnk(len), qnk(len), gznk(len), plcl(len)
151
152  ! local variables:
153  INTEGER i, k
154  INTEGER ihmin(len)
155  REAL work(len)
156  REAL pnk(len), qsnk(len), rh(len), chi(len)
157
158  ! -------------------------------------------------------------------
159  ! --- Find level of minimum moist static energy
160  ! --- If level of minimum moist static energy coincides with
161  ! --- or is lower than minimum allowable parcel origin level,
162  ! --- set iflag to 6.
163  ! -------------------------------------------------------------------
164
165  DO i = 1, len
166    work(i) = 1.0E12
167    ihmin(i) = nl
168  END DO
169  DO k = 2, nlp
170    DO i = 1, len
171      IF ((hm(i,k)<work(i)) .AND. (hm(i,k)<hm(i,k-1))) THEN
172        work(i) = hm(i, k)
173        ihmin(i) = k
174      END IF
175    END DO
176  END DO
177  DO i = 1, len
178    ihmin(i) = min(ihmin(i), nlm)
179    IF (ihmin(i)<=minorig) THEN
180      iflag(i) = 6
181    END IF
182  END DO
183
184  ! -------------------------------------------------------------------
185  ! --- Find that model level below the level of minimum moist static
186  ! --- energy that has the maximum value of moist static energy
187  ! -------------------------------------------------------------------
188
189  DO i = 1, len
190    work(i) = hm(i, minorig)
191    nk(i) = minorig
192  END DO
193  DO k = minorig + 1, nl
194    DO i = 1, len
195      IF ((hm(i,k)>work(i)) .AND. (k<=ihmin(i))) THEN
196        work(i) = hm(i, k)
197        nk(i) = k
198      END IF
199    END DO
200  END DO
201  ! -------------------------------------------------------------------
202  ! --- Check whether parcel level temperature and specific humidity
203  ! --- are reasonable
204  ! -------------------------------------------------------------------
205  DO i = 1, len
206    IF (((t(i,nk(i))<250.0) .OR. (q(i,nk(i))<=0.0) .OR. (p(i,ihmin(i))< &
207      400.0)) .AND. (iflag(i)==0)) iflag(i) = 7
208  END DO
209  ! -------------------------------------------------------------------
210  ! --- Calculate lifted condensation level of air at parcel origin level
211  ! --- (Within 0.2% of formula of Bolton, MON. WEA. REV.,1980)
212  ! -------------------------------------------------------------------
213  DO i = 1, len
214    tnk(i) = t(i, nk(i))
215    qnk(i) = q(i, nk(i))
216    gznk(i) = gz(i, nk(i))
217    pnk(i) = p(i, nk(i))
218    qsnk(i) = qs(i, nk(i))
219
220    rh(i) = qnk(i)/qsnk(i)
221    rh(i) = min(1.0, rh(i))
222    chi(i) = tnk(i)/(1669.0-122.0*rh(i)-tnk(i))
223    plcl(i) = pnk(i)*(rh(i)**chi(i))
224    IF (((plcl(i)<200.0) .OR. (plcl(i)>=2000.0)) .AND. (iflag(i)==0)) iflag(i &
225      ) = 8
226  END DO
227  ! -------------------------------------------------------------------
228  ! --- Calculate first level above lcl (=icb)
229  ! -------------------------------------------------------------------
230  DO i = 1, len
231    icb(i) = nlm
232  END DO
233
234  DO k = minorig, nl
235    DO i = 1, len
236      IF ((k>=(nk(i)+1)) .AND. (p(i,k)<plcl(i))) icb(i) = min(icb(i), k)
237    END DO
238  END DO
239
240  DO i = 1, len
241    IF ((icb(i)>=nlm) .AND. (iflag(i)==0)) iflag(i) = 9
242  END DO
243
244  ! Compute icbmax.
245
246  icbmax = 2
247  DO i = 1, len
248    icbmax = max(icbmax, icb(i))
249  END DO
250
251  RETURN
252END SUBROUTINE cv_feed
253
254SUBROUTINE cv_undilute1(len, nd, t, q, qs, gz, p, nk, icb, icbmax, tp, tvp, &
255    clw)
256  USE lmdz_cv_ini, ONLY : cl,clmcpv,cpd,cpv,eps,epsi,lv0,minorig,rrv,t0
257
258  IMPLICIT NONE
259
260
261  ! inputs:
262  INTEGER len, nd
263  INTEGER nk(len), icb(len), icbmax
264  REAL t(len, nd), q(len, nd), qs(len, nd), gz(len, nd)
265  REAL p(len, nd)
266
267  ! outputs:
268  REAL tp(len, nd), tvp(len, nd), clw(len, nd)
269
270  ! local variables:
271  INTEGER i, k
272  REAL tg, qg, alv, s, ahg, tc, denom, es, rg
273  REAL ah0(len), cpp(len)
274  REAL tnk(len), qnk(len), gznk(len), ticb(len), gzicb(len)
275
276  ! -------------------------------------------------------------------
277  ! --- Calculates the lifted parcel virtual temperature at nk,
278  ! --- the actual temperature, and the adiabatic
279  ! --- liquid water content. The procedure is to solve the equation.
280  ! cp*tp+L*qp+phi=cp*tnk+L*qnk+gznk.
281  ! -------------------------------------------------------------------
282
283  DO i = 1, len
284    tnk(i) = t(i, nk(i))
285    qnk(i) = q(i, nk(i))
286    gznk(i) = gz(i, nk(i))
287    ticb(i) = t(i, icb(i))
288    gzicb(i) = gz(i, icb(i))
289  END DO
290
291  ! ***  Calculate certain parcel quantities, including static energy   ***
292
293  DO i = 1, len
294    ah0(i) = (cpd*(1.-qnk(i))+cl*qnk(i))*tnk(i) + qnk(i)*(lv0-clmcpv*(tnk(i)- &
295      273.15)) + gznk(i)
296    cpp(i) = cpd*(1.-qnk(i)) + qnk(i)*cpv
297  END DO
298
299  ! ***   Calculate lifted parcel quantities below cloud base   ***
300
301  DO k = minorig, icbmax - 1
302    DO i = 1, len
303      tp(i, k) = tnk(i) - (gz(i,k)-gznk(i))/cpp(i)
304      tvp(i, k) = tp(i, k)*(1.+qnk(i)*epsi)
305    END DO
306  END DO
307
308  ! ***  Find lifted parcel quantities above cloud base    ***
309
310  DO i = 1, len
311    tg = ticb(i)
312    qg = qs(i, icb(i))
313    alv = lv0 - clmcpv*(ticb(i)-t0)
314
315    ! First iteration.
316
317    s = cpd + alv*alv*qg/(rrv*ticb(i)*ticb(i))
318    s = 1./s
319    ahg = cpd*tg + (cl-cpd)*qnk(i)*ticb(i) + alv*qg + gzicb(i)
320    tg = tg + s*(ah0(i)-ahg)
321    tg = max(tg, 35.0)
322    tc = tg - t0
323    denom = 243.5 + tc
324    IF (tc>=0.0) THEN
325      es = 6.112*exp(17.67*tc/denom)
326    ELSE
327      es = exp(23.33086-6111.72784/tg+0.15215*log(tg))
328    END IF
329    qg = eps*es/(p(i,icb(i))-es*(1.-eps))
330
331    ! Second iteration.
332
333    s = cpd + alv*alv*qg/(rrv*ticb(i)*ticb(i))
334    s = 1./s
335    ahg = cpd*tg + (cl-cpd)*qnk(i)*ticb(i) + alv*qg + gzicb(i)
336    tg = tg + s*(ah0(i)-ahg)
337    tg = max(tg, 35.0)
338    tc = tg - t0
339    denom = 243.5 + tc
340    IF (tc>=0.0) THEN
341      es = 6.112*exp(17.67*tc/denom)
342    ELSE
343      es = exp(23.33086-6111.72784/tg+0.15215*log(tg))
344    END IF
345    qg = eps*es/(p(i,icb(i))-es*(1.-eps))
346
347    alv = lv0 - clmcpv*(ticb(i)-273.15)
348    tp(i, icb(i)) = (ah0(i)-(cl-cpd)*qnk(i)*ticb(i)-gz(i,icb(i))-alv*qg)/cpd
349    clw(i, icb(i)) = qnk(i) - qg
350    clw(i, icb(i)) = max(0.0, clw(i,icb(i)))
351    rg = qg/(1.-qnk(i))
352    tvp(i, icb(i)) = tp(i, icb(i))*(1.+rg*epsi)
353  END DO
354
355  DO k = minorig, icbmax
356    DO i = 1, len
357      tvp(i, k) = tvp(i, k) - tp(i, k)*qnk(i)
358    END DO
359  END DO
360
361  RETURN
362END SUBROUTINE cv_undilute1
363
364SUBROUTINE cv_trigger(len, nd, icb, cbmf, tv, tvp, iflag)
365   USE lmdz_cv_ini, ONLY : dtmax
366    IMPLICIT NONE
367
368  ! -------------------------------------------------------------------
369  ! --- Test for instability.
370  ! --- If there was no convection at last time step and parcel
371  ! --- is stable at icb, then set iflag to 4.
372  ! -------------------------------------------------------------------
373
374
375  ! inputs:
376  INTEGER len, nd, icb(len)
377  REAL cbmf(len), tv(len, nd), tvp(len, nd)
378
379  ! outputs:
380  INTEGER iflag(len) ! also an input
381
382  ! local variables:
383  INTEGER i
384
385
386  DO i = 1, len
387    IF ((cbmf(i)==0.0) .AND. (iflag(i)==0) .AND. (tvp(i, &
388      icb(i))<=(tv(i,icb(i))-dtmax))) iflag(i) = 4
389  END DO
390
391  RETURN
392END SUBROUTINE cv_trigger
393
394SUBROUTINE cv_compress(len, nloc, ncum, nd, iflag1, nk1, icb1, cbmf1, plcl1, &
395    tnk1, qnk1, gznk1, t1, q1, qs1, u1, v1, gz1, h1, lv1, cpn1, p1, ph1, tv1, &
396    tp1, tvp1, clw1, iflag, nk, icb, cbmf, plcl, tnk, qnk, gznk, t, q, qs, u, &
397    v, gz, h, lv, cpn, p, ph, tv, tp, tvp, clw, dph)
398   USE lmdz_cv_ini, ONLY : nl
399    USE print_control_mod, ONLY: lunout
400  IMPLICIT NONE
401
402
403  ! inputs:
404  INTEGER len, ncum, nd, nloc
405  INTEGER iflag1(len), nk1(len), icb1(len)
406  REAL cbmf1(len), plcl1(len), tnk1(len), qnk1(len), gznk1(len)
407  REAL t1(len, nd), q1(len, nd), qs1(len, nd), u1(len, nd), v1(len, nd)
408  REAL gz1(len, nd), h1(len, nd), lv1(len, nd), cpn1(len, nd)
409  REAL p1(len, nd), ph1(len, nd+1), tv1(len, nd), tp1(len, nd)
410  REAL tvp1(len, nd), clw1(len, nd)
411
412  ! outputs:
413  INTEGER iflag(nloc), nk(nloc), icb(nloc)
414  REAL cbmf(nloc), plcl(nloc), tnk(nloc), qnk(nloc), gznk(nloc)
415  REAL t(nloc, nd), q(nloc, nd), qs(nloc, nd), u(nloc, nd), v(nloc, nd)
416  REAL gz(nloc, nd), h(nloc, nd), lv(nloc, nd), cpn(nloc, nd)
417  REAL p(nloc, nd), ph(nloc, nd+1), tv(nloc, nd), tp(nloc, nd)
418  REAL tvp(nloc, nd), clw(nloc, nd)
419  REAL dph(nloc, nd)
420
421  ! local variables:
422  INTEGER i, k, nn
423  CHARACTER (LEN=20) :: modname = 'cv_compress'
424  CHARACTER (LEN=80) :: abort_message
425
426
427  DO k = 1, nl + 1
428    nn = 0
429    DO i = 1, len
430      IF (iflag1(i)==0) THEN
431        nn = nn + 1
432        t(nn, k) = t1(i, k)
433        q(nn, k) = q1(i, k)
434        qs(nn, k) = qs1(i, k)
435        u(nn, k) = u1(i, k)
436        v(nn, k) = v1(i, k)
437        gz(nn, k) = gz1(i, k)
438        h(nn, k) = h1(i, k)
439        lv(nn, k) = lv1(i, k)
440        cpn(nn, k) = cpn1(i, k)
441        p(nn, k) = p1(i, k)
442        ph(nn, k) = ph1(i, k)
443        tv(nn, k) = tv1(i, k)
444        tp(nn, k) = tp1(i, k)
445        tvp(nn, k) = tvp1(i, k)
446        clw(nn, k) = clw1(i, k)
447      END IF
448    END DO
449  END DO
450
451  IF (nn/=ncum) THEN
452    WRITE (lunout, *) 'strange! nn not equal to ncum: ', nn, ncum
453    abort_message = ''
454    CALL abort_physic(modname, abort_message, 1)
455  END IF
456
457  nn = 0
458  DO i = 1, len
459    IF (iflag1(i)==0) THEN
460      nn = nn + 1
461      cbmf(nn) = cbmf1(i)
462      plcl(nn) = plcl1(i)
463      tnk(nn) = tnk1(i)
464      qnk(nn) = qnk1(i)
465      gznk(nn) = gznk1(i)
466      nk(nn) = nk1(i)
467      icb(nn) = icb1(i)
468      iflag(nn) = iflag1(i)
469    END IF
470  END DO
471
472  DO k = 1, nl
473    DO i = 1, ncum
474      dph(i, k) = ph(i, k) - ph(i, k+1)
475    END DO
476  END DO
477
478  RETURN
479END SUBROUTINE cv_compress
480
481SUBROUTINE cv_undilute2(nloc, ncum, nd, icb, nk, tnk, qnk, gznk, t, q, qs, &
482    gz, p, dph, h, tv, lv, inb, inb1, tp, tvp, clw, hp, ep, sigp, frac)
483  USE lmdz_cv_ini, ONLY : nl,cl,clmcpv,cpd,cpv,elcrit,eps,epsi,lv0,minorig,nlp,rrv,sigs,t0,tlcrit
484
485  IMPLICIT NONE
486
487  ! ---------------------------------------------------------------------
488  ! Purpose:
489  ! FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
490  ! &
491  ! COMPUTE THE PRECIPITATION EFFICIENCIES AND THE
492  ! FRACTION OF PRECIPITATION FALLING OUTSIDE OF CLOUD
493  ! &
494  ! FIND THE LEVEL OF NEUTRAL BUOYANCY
495  ! ---------------------------------------------------------------------
496
497  ! inputs:
498  INTEGER ncum, nd, nloc
499  INTEGER icb(nloc), nk(nloc)
500  REAL t(nloc, nd), q(nloc, nd), qs(nloc, nd), gz(nloc, nd)
501  REAL p(nloc, nd), dph(nloc, nd)
502  REAL tnk(nloc), qnk(nloc), gznk(nloc)
503  REAL lv(nloc, nd), tv(nloc, nd), h(nloc, nd)
504
505  ! outputs:
506  INTEGER inb(nloc), inb1(nloc)
507  REAL tp(nloc, nd), tvp(nloc, nd), clw(nloc, nd)
508  REAL ep(nloc, nd), sigp(nloc, nd), hp(nloc, nd)
509  REAL frac(nloc)
510
511  ! local variables:
512  INTEGER i, k
513  REAL tg, qg, ahg, alv, s, tc, es, denom, rg, tca, elacrit
514  REAL by, defrac
515  REAL ah0(nloc), cape(nloc), capem(nloc), byp(nloc)
516  LOGICAL lcape(nloc)
517
518  ! =====================================================================
519  ! --- SOME INITIALIZATIONS
520  ! =====================================================================
521
522  DO k = 1, nl
523    DO i = 1, ncum
524      ep(i, k) = 0.0
525      sigp(i, k) = sigs
526    END DO
527  END DO
528
529  ! =====================================================================
530  ! --- FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
531  ! =====================================================================
532
533  ! ---       The procedure is to solve the equation.
534  ! cp*tp+L*qp+phi=cp*tnk+L*qnk+gznk.
535
536  ! ***  Calculate certain parcel quantities, including static energy   ***
537
538
539  DO i = 1, ncum
540    ah0(i) = (cpd*(1.-qnk(i))+cl*qnk(i))*tnk(i) + qnk(i)*(lv0-clmcpv*(tnk(i)- &
541      t0)) + gznk(i)
542  END DO
543
544
545  ! ***  Find lifted parcel quantities above cloud base    ***
546
547
548  DO k = minorig + 1, nl
549    DO i = 1, ncum
550      IF (k>=(icb(i)+1)) THEN
551        tg = t(i, k)
552        qg = qs(i, k)
553        alv = lv0 - clmcpv*(t(i,k)-t0)
554
555        ! First iteration.
556
557        s = cpd + alv*alv*qg/(rrv*t(i,k)*t(i,k))
558        s = 1./s
559        ahg = cpd*tg + (cl-cpd)*qnk(i)*t(i, k) + alv*qg + gz(i, k)
560        tg = tg + s*(ah0(i)-ahg)
561        tg = max(tg, 35.0)
562        tc = tg - t0
563        denom = 243.5 + tc
564        IF (tc>=0.0) THEN
565          es = 6.112*exp(17.67*tc/denom)
566        ELSE
567          es = exp(23.33086-6111.72784/tg+0.15215*log(tg))
568        END IF
569        qg = eps*es/(p(i,k)-es*(1.-eps))
570
571        ! Second iteration.
572
573        s = cpd + alv*alv*qg/(rrv*t(i,k)*t(i,k))
574        s = 1./s
575        ahg = cpd*tg + (cl-cpd)*qnk(i)*t(i, k) + alv*qg + gz(i, k)
576        tg = tg + s*(ah0(i)-ahg)
577        tg = max(tg, 35.0)
578        tc = tg - t0
579        denom = 243.5 + tc
580        IF (tc>=0.0) THEN
581          es = 6.112*exp(17.67*tc/denom)
582        ELSE
583          es = exp(23.33086-6111.72784/tg+0.15215*log(tg))
584        END IF
585        qg = eps*es/(p(i,k)-es*(1.-eps))
586
587        alv = lv0 - clmcpv*(t(i,k)-t0)
588        ! print*,'cpd dans convect2 ',cpd
589        ! print*,'tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd'
590        ! print*,tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd
591        tp(i, k) = (ah0(i)-(cl-cpd)*qnk(i)*t(i,k)-gz(i,k)-alv*qg)/cpd
592        ! if (.not.cpd.gt.1000.) then
593        ! print*,'CPD=',cpd
594        ! stop
595        ! endif
596        clw(i, k) = qnk(i) - qg
597        clw(i, k) = max(0.0, clw(i,k))
598        rg = qg/(1.-qnk(i))
599        tvp(i, k) = tp(i, k)*(1.+rg*epsi)
600      END IF
601    END DO
602  END DO
603
604  ! =====================================================================
605  ! --- SET THE PRECIPITATION EFFICIENCIES AND THE FRACTION OF
606  ! --- PRECIPITATION FALLING OUTSIDE OF CLOUD
607  ! --- THESE MAY BE FUNCTIONS OF TP(I), P(I) AND CLW(I)
608  ! =====================================================================
609
610  DO k = minorig + 1, nl
611    DO i = 1, ncum
612      IF (k>=(nk(i)+1)) THEN
613        tca = tp(i, k) - t0
614        IF (tca>=0.0) THEN
615          elacrit = elcrit
616        ELSE
617          elacrit = elcrit*(1.0-tca/tlcrit)
618        END IF
619        elacrit = max(elacrit, 0.0)
620        ep(i, k) = 1.0 - elacrit/max(clw(i,k), 1.0E-8)
621        ep(i, k) = max(ep(i,k), 0.0)
622        ep(i, k) = min(ep(i,k), 1.0)
623        sigp(i, k) = sigs
624      END IF
625    END DO
626  END DO
627
628  ! =====================================================================
629  ! --- CALCULATE VIRTUAL TEMPERATURE AND LIFTED PARCEL
630  ! --- VIRTUAL TEMPERATURE
631  ! =====================================================================
632
633  DO k = minorig + 1, nl
634    DO i = 1, ncum
635      IF (k>=(icb(i)+1)) THEN
636        tvp(i, k) = tvp(i, k)*(1.0-qnk(i)+ep(i,k)*clw(i,k))
637        ! print*,'i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)'
638        ! print*, i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)
639      END IF
640    END DO
641  END DO
642  DO i = 1, ncum
643    tvp(i, nlp) = tvp(i, nl) - (gz(i,nlp)-gz(i,nl))/cpd
644  END DO
645
646  ! =====================================================================
647  ! --- FIND THE FIRST MODEL LEVEL (INB1) ABOVE THE PARCEL'S
648  ! --- HIGHEST LEVEL OF NEUTRAL BUOYANCY
649  ! --- AND THE HIGHEST LEVEL OF POSITIVE CAPE (INB)
650  ! =====================================================================
651
652  DO i = 1, ncum
653    cape(i) = 0.0
654    capem(i) = 0.0
655    inb(i) = icb(i) + 1
656    inb1(i) = inb(i)
657  END DO
658
659  ! Originial Code
660
661  ! do 530 k=minorig+1,nl-1
662  ! do 520 i=1,ncum
663  ! if(k.ge.(icb(i)+1))then
664  ! by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
665  ! byp=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
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  ! endif
672  ! endif
673  ! 520    continue
674  ! 530  continue
675  ! do 540 i=1,ncum
676  ! byp=(tvp(i,nl)-tv(i,nl))*dph(i,nl)/p(i,nl)
677  ! cape(i)=capem(i)+byp
678  ! defrac=capem(i)-cape(i)
679  ! defrac=max(defrac,0.001)
680  ! frac(i)=-cape(i)/defrac
681  ! frac(i)=min(frac(i),1.0)
682  ! frac(i)=max(frac(i),0.0)
683  ! 540   continue
684
685  ! K Emanuel fix
686
687  ! call zilch(byp,ncum)
688  ! do 530 k=minorig+1,nl-1
689  ! do 520 i=1,ncum
690  ! if(k.ge.(icb(i)+1))then
691  ! by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
692  ! cape(i)=cape(i)+by
693  ! if(by.ge.0.0)inb1(i)=k+1
694  ! if(cape(i).gt.0.0)then
695  ! inb(i)=k+1
696  ! capem(i)=cape(i)
697  ! byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
698  ! endif
699  ! endif
700  ! 520    continue
701  ! 530  continue
702  ! do 540 i=1,ncum
703  ! inb(i)=max(inb(i),inb1(i))
704  ! cape(i)=capem(i)+byp(i)
705  ! defrac=capem(i)-cape(i)
706  ! defrac=max(defrac,0.001)
707  ! frac(i)=-cape(i)/defrac
708  ! frac(i)=min(frac(i),1.0)
709  ! frac(i)=max(frac(i),0.0)
710  ! 540   continue
711
712  ! J Teixeira fix
713
714  CALL zilch(byp, ncum)
715  DO i = 1, ncum
716    lcape(i) = .TRUE.
717  END DO
718  DO k = minorig + 1, nl - 1
719    DO i = 1, ncum
720      IF (cape(i)<0.0) lcape(i) = .FALSE.
721      IF ((k>=(icb(i)+1)) .AND. lcape(i)) THEN
722        by = (tvp(i,k)-tv(i,k))*dph(i, k)/p(i, k)
723        byp(i) = (tvp(i,k+1)-tv(i,k+1))*dph(i, k+1)/p(i, k+1)
724        cape(i) = cape(i) + by
725        IF (by>=0.0) inb1(i) = k + 1
726        IF (cape(i)>0.0) THEN
727          inb(i) = k + 1
728          capem(i) = cape(i)
729        END IF
730      END IF
731    END DO
732  END DO
733  DO i = 1, ncum
734    cape(i) = capem(i) + byp(i)
735    defrac = capem(i) - cape(i)
736    defrac = max(defrac, 0.001)
737    frac(i) = -cape(i)/defrac
738    frac(i) = min(frac(i), 1.0)
739    frac(i) = max(frac(i), 0.0)
740  END DO
741
742  ! =====================================================================
743  ! ---   CALCULATE LIQUID WATER STATIC ENERGY OF LIFTED PARCEL
744  ! =====================================================================
745
746  ! initialization:
747  DO i = 1, ncum*nlp
748    hp(i, 1) = h(i, 1)
749  END DO
750
751  DO k = minorig + 1, nl
752    DO i = 1, ncum
753      IF ((k>=icb(i)) .AND. (k<=inb(i))) THEN
754        hp(i, k) = h(i, nk(i)) + (lv(i,k)+(cpd-cpv)*t(i,k))*ep(i, k)*clw(i, k &
755          )
756      END IF
757    END DO
758  END DO
759
760  RETURN
761END SUBROUTINE cv_undilute2
762
763SUBROUTINE cv_closure(nloc, ncum, nd, nk, icb, tv, tvp, p, ph, dph, plcl, &
764    cpn, iflag, cbmf)
765  USE lmdz_cv_ini, ONLY : alpha,damp,dtmax,minorig,rrd
766
767  IMPLICIT NONE
768
769  ! inputs:
770  INTEGER ncum, nd, nloc
771  INTEGER nk(nloc), icb(nloc)
772  REAL tv(nloc, nd), tvp(nloc, nd), p(nloc, nd), dph(nloc, nd)
773  REAL ph(nloc, nd+1) ! caution nd instead ndp1 to be consistent...
774  REAL plcl(nloc), cpn(nloc, nd)
775
776  ! outputs:
777  INTEGER iflag(nloc)
778  REAL cbmf(nloc) ! also an input
779
780  ! local variables:
781  INTEGER i, k, icbmax
782  REAL dtpbl(nloc), dtmin(nloc), tvpplcl(nloc), tvaplcl(nloc)
783  REAL work(nloc)
784
785  ! -------------------------------------------------------------------
786  ! Compute icbmax.
787  ! -------------------------------------------------------------------
788
789  icbmax = 2
790  DO i = 1, ncum
791    icbmax = max(icbmax, icb(i))
792  END DO
793
794  ! =====================================================================
795  ! ---  CALCULATE CLOUD BASE MASS FLUX
796  ! =====================================================================
797
798  ! tvpplcl = parcel temperature lifted adiabatically from level
799  ! icb-1 to the LCL.
800  ! tvaplcl = virtual temperature at the LCL.
801
802  DO i = 1, ncum
803    dtpbl(i) = 0.0
804    tvpplcl(i) = tvp(i, icb(i)-1) - rrd*tvp(i, icb(i)-1)*(p(i,icb(i)-1)-plcl( &
805      i))/(cpn(i,icb(i)-1)*p(i,icb(i)-1))
806    tvaplcl(i) = tv(i, icb(i)) + (tvp(i,icb(i))-tvp(i,icb(i)+1))*(plcl(i)-p(i &
807      ,icb(i)))/(p(i,icb(i))-p(i,icb(i)+1))
808  END DO
809
810  ! -------------------------------------------------------------------
811  ! --- Interpolate difference between lifted parcel and
812  ! --- environmental temperatures to lifted condensation level
813  ! -------------------------------------------------------------------
814
815  ! dtpbl = average of tvp-tv in the PBL (k=nk to icb-1).
816
817  DO k = minorig, icbmax
818    DO i = 1, ncum
819      IF ((k>=nk(i)) .AND. (k<=(icb(i)-1))) THEN
820        dtpbl(i) = dtpbl(i) + (tvp(i,k)-tv(i,k))*dph(i, k)
821      END IF
822    END DO
823  END DO
824  DO i = 1, ncum
825    dtpbl(i) = dtpbl(i)/(ph(i,nk(i))-ph(i,icb(i)))
826    dtmin(i) = tvpplcl(i) - tvaplcl(i) + dtmax + dtpbl(i)
827  END DO
828
829  ! -------------------------------------------------------------------
830  ! --- Adjust cloud base mass flux
831  ! -------------------------------------------------------------------
832
833  DO i = 1, ncum
834    work(i) = cbmf(i)
835    cbmf(i) = max(0.0, (1.0-damp)*cbmf(i)+0.1*alpha*dtmin(i))
836    IF ((work(i)==0.0) .AND. (cbmf(i)==0.0)) THEN
837      iflag(i) = 3
838    END IF
839  END DO
840
841  RETURN
842END SUBROUTINE cv_closure
843
844SUBROUTINE cv_mixing(nloc, ncum, nd, icb, nk, inb, inb1, ph, t, q, qs, u, v, &
845    h, lv, qnk, hp, tv, tvp, ep, clw, cbmf, m, ment, qent, uent, vent, nent, &
846    sij, elij)
847  USE lmdz_cv_ini, ONLY : cpd,cpv,entp,minorig,nl,nlp,rrv
848
849  IMPLICIT NONE
850
851
852  ! inputs:
853  INTEGER ncum, nd, nloc
854  INTEGER icb(nloc), inb(nloc), inb1(nloc), nk(nloc)
855  REAL cbmf(nloc), qnk(nloc)
856  REAL ph(nloc, nd+1)
857  REAL t(nloc, nd), q(nloc, nd), qs(nloc, nd), lv(nloc, nd)
858  REAL u(nloc, nd), v(nloc, nd), h(nloc, nd), hp(nloc, nd)
859  REAL tv(nloc, nd), tvp(nloc, nd), ep(nloc, nd), clw(nloc, nd)
860
861  ! outputs:
862  INTEGER nent(nloc, nd)
863  REAL m(nloc, nd), ment(nloc, nd, nd), qent(nloc, nd, nd)
864  REAL uent(nloc, nd, nd), vent(nloc, nd, nd)
865  REAL sij(nloc, nd, nd), elij(nloc, nd, nd)
866
867  ! local variables:
868  INTEGER i, j, k, ij
869  INTEGER num1, num2
870  REAL dbo, qti, bf2, anum, denom, dei, altem, cwat, stemp
871  REAL alt, qp1, smid, sjmin, sjmax, delp, delm
872  REAL work(nloc), asij(nloc), smin(nloc), scrit(nloc)
873  REAL bsum(nloc, nd)
874  LOGICAL lwork(nloc)
875
876  ! =====================================================================
877  ! --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS
878  ! =====================================================================
879
880  DO i = 1, ncum*nlp
881    nent(i, 1) = 0
882    m(i, 1) = 0.0
883  END DO
884
885  DO k = 1, nlp
886    DO j = 1, nlp
887      DO i = 1, ncum
888        qent(i, k, j) = q(i, j)
889        uent(i, k, j) = u(i, j)
890        vent(i, k, j) = v(i, j)
891        elij(i, k, j) = 0.0
892        ment(i, k, j) = 0.0
893        sij(i, k, j) = 0.0
894      END DO
895    END DO
896  END DO
897
898  ! -------------------------------------------------------------------
899  ! --- Calculate rates of mixing,  m(i)
900  ! -------------------------------------------------------------------
901
902  CALL zilch(work, ncum)
903
904  DO j = minorig + 1, nl
905    DO i = 1, ncum
906      IF ((j>=(icb(i)+1)) .AND. (j<=inb(i))) THEN
907        k = min(j, inb1(i))
908        dbo = abs(tv(i,k+1)-tvp(i,k+1)-tv(i,k-1)+tvp(i,k-1)) + &
909          entp*0.04*(ph(i,k)-ph(i,k+1))
910        work(i) = work(i) + dbo
911        m(i, j) = cbmf(i)*dbo
912      END IF
913    END DO
914  END DO
915  DO k = minorig + 1, nl
916    DO i = 1, ncum
917      IF ((k>=(icb(i)+1)) .AND. (k<=inb(i))) THEN
918        m(i, k) = m(i, k)/work(i)
919      END IF
920    END DO
921  END DO
922
923
924  ! =====================================================================
925  ! --- CALCULATE ENTRAINED AIR MASS FLUX (ment), TOTAL WATER MIXING
926  ! --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
927  ! --- FRACTION (sij)
928  ! =====================================================================
929
930
931  DO i = minorig + 1, nl
932    DO j = minorig + 1, nl
933      DO ij = 1, ncum
934        IF ((i>=(icb(ij)+1)) .AND. (j>=icb(ij)) .AND. (i<=inb(ij)) .AND. (j<= &
935            inb(ij))) THEN
936          qti = qnk(ij) - ep(ij, i)*clw(ij, i)
937          bf2 = 1. + lv(ij, j)*lv(ij, j)*qs(ij, j)/(rrv*t(ij,j)*t(ij,j)*cpd)
938          anum = h(ij, j) - hp(ij, i) + (cpv-cpd)*t(ij, j)*(qti-q(ij,j))
939          denom = h(ij, i) - hp(ij, i) + (cpd-cpv)*(q(ij,i)-qti)*t(ij, j)
940          dei = denom
941          IF (abs(dei)<0.01) dei = 0.01
942          sij(ij, i, j) = anum/dei
943          sij(ij, i, i) = 1.0
944          altem = sij(ij, i, j)*q(ij, i) + (1.-sij(ij,i,j))*qti - qs(ij, j)
945          altem = altem/bf2
946          cwat = clw(ij, j)*(1.-ep(ij,j))
947          stemp = sij(ij, i, j)
948          IF ((stemp<0.0 .OR. stemp>1.0 .OR. altem>cwat) .AND. j>i) THEN
949            anum = anum - lv(ij, j)*(qti-qs(ij,j)-cwat*bf2)
950            denom = denom + lv(ij, j)*(q(ij,i)-qti)
951            IF (abs(denom)<0.01) denom = 0.01
952            sij(ij, i, j) = anum/denom
953            altem = sij(ij, i, j)*q(ij, i) + (1.-sij(ij,i,j))*qti - qs(ij, j)
954            altem = altem - (bf2-1.)*cwat
955          END IF
956          IF (sij(ij,i,j)>0.0 .AND. sij(ij,i,j)<0.9) THEN
957            qent(ij, i, j) = sij(ij, i, j)*q(ij, i) + (1.-sij(ij,i,j))*qti
958            uent(ij, i, j) = sij(ij, i, j)*u(ij, i) + &
959              (1.-sij(ij,i,j))*u(ij, nk(ij))
960            vent(ij, i, j) = sij(ij, i, j)*v(ij, i) + &
961              (1.-sij(ij,i,j))*v(ij, nk(ij))
962            elij(ij, i, j) = altem
963            elij(ij, i, j) = max(0.0, elij(ij,i,j))
964            ment(ij, i, j) = m(ij, i)/(1.-sij(ij,i,j))
965            nent(ij, i) = nent(ij, i) + 1
966          END IF
967          sij(ij, i, j) = max(0.0, sij(ij,i,j))
968          sij(ij, i, j) = min(1.0, sij(ij,i,j))
969        END IF
970      END DO
971    END DO
972
973    ! ***   If no air can entrain at level i assume that updraft detrains
974    ! ***
975    ! ***   at that level and calculate detrained air flux and properties
976    ! ***
977
978    DO ij = 1, ncum
979      IF ((i>=(icb(ij)+1)) .AND. (i<=inb(ij)) .AND. (nent(ij,i)==0)) THEN
980        ment(ij, i, i) = m(ij, i)
981        qent(ij, i, i) = q(ij, nk(ij)) - ep(ij, i)*clw(ij, i)
982        uent(ij, i, i) = u(ij, nk(ij))
983        vent(ij, i, i) = v(ij, nk(ij))
984        elij(ij, i, i) = clw(ij, i)
985        sij(ij, i, i) = 1.0
986      END IF
987    END DO
988  END DO
989
990  DO i = 1, ncum
991    sij(i, inb(i), inb(i)) = 1.0
992  END DO
993
994  ! =====================================================================
995  ! ---  NORMALIZE ENTRAINED AIR MASS FLUXES
996  ! ---  TO REPRESENT EQUAL PROBABILITIES OF MIXING
997  ! =====================================================================
998
999  CALL zilch(bsum, ncum*nlp)
1000  DO ij = 1, ncum
1001    lwork(ij) = .FALSE.
1002  END DO
1003  DO i = minorig + 1, nl
1004
1005    num1 = 0
1006    DO ij = 1, ncum
1007      IF ((i>=icb(ij)+1) .AND. (i<=inb(ij))) num1 = num1 + 1
1008    END DO
1009    IF (num1<=0) GO TO 789
1010
1011    DO ij = 1, ncum
1012      IF ((i>=icb(ij)+1) .AND. (i<=inb(ij))) THEN
1013        lwork(ij) = (nent(ij,i)/=0)
1014        qp1 = q(ij, nk(ij)) - ep(ij, i)*clw(ij, i)
1015        anum = h(ij, i) - hp(ij, i) - lv(ij, i)*(qp1-qs(ij,i))
1016        denom = h(ij, i) - hp(ij, i) + lv(ij, i)*(q(ij,i)-qp1)
1017        IF (abs(denom)<0.01) denom = 0.01
1018        scrit(ij) = anum/denom
1019        alt = qp1 - qs(ij, i) + scrit(ij)*(q(ij,i)-qp1)
1020        IF (scrit(ij)<0.0 .OR. alt<0.0) scrit(ij) = 1.0
1021        asij(ij) = 0.0
1022        smin(ij) = 1.0
1023      END IF
1024    END DO
1025    DO j = minorig, nl
1026
1027      num2 = 0
1028      DO ij = 1, ncum
1029        IF ((i>=icb(ij)+1) .AND. (i<=inb(ij)) .AND. (j>=icb( &
1030          ij)) .AND. (j<=inb(ij)) .AND. lwork(ij)) num2 = num2 + 1
1031      END DO
1032      IF (num2<=0) GO TO 783
1033
1034      DO ij = 1, ncum
1035        IF ((i>=icb(ij)+1) .AND. (i<=inb(ij)) .AND. (j>=icb( &
1036            ij)) .AND. (j<=inb(ij)) .AND. lwork(ij)) THEN
1037          IF (sij(ij,i,j)>0.0 .AND. sij(ij,i,j)<0.9) THEN
1038            IF (j>i) THEN
1039              smid = min(sij(ij,i,j), scrit(ij))
1040              sjmax = smid
1041              sjmin = smid
1042              IF (smid<smin(ij) .AND. sij(ij,i,j+1)<smid) THEN
1043                smin(ij) = smid
1044                sjmax = min(sij(ij,i,j+1), sij(ij,i,j), scrit(ij))
1045                sjmin = max(sij(ij,i,j-1), sij(ij,i,j))
1046                sjmin = min(sjmin, scrit(ij))
1047              END IF
1048            ELSE
1049              sjmax = max(sij(ij,i,j+1), scrit(ij))
1050              smid = max(sij(ij,i,j), scrit(ij))
1051              sjmin = 0.0
1052              IF (j>1) sjmin = sij(ij, i, j-1)
1053              sjmin = max(sjmin, scrit(ij))
1054            END IF
1055            delp = abs(sjmax-smid)
1056            delm = abs(sjmin-smid)
1057            asij(ij) = asij(ij) + (delp+delm)*(ph(ij,j)-ph(ij,j+1))
1058            ment(ij, i, j) = ment(ij, i, j)*(delp+delm)*(ph(ij,j)-ph(ij,j+1))
1059          END IF
1060        END IF
1061      END DO
1062783 END DO
1063    DO ij = 1, ncum
1064      IF ((i>=icb(ij)+1) .AND. (i<=inb(ij)) .AND. lwork(ij)) THEN
1065        asij(ij) = max(1.0E-21, asij(ij))
1066        asij(ij) = 1.0/asij(ij)
1067        bsum(ij, i) = 0.0
1068      END IF
1069    END DO
1070    DO j = minorig, nl + 1
1071      DO ij = 1, ncum
1072        IF ((i>=icb(ij)+1) .AND. (i<=inb(ij)) .AND. (j>=icb( &
1073            ij)) .AND. (j<=inb(ij)) .AND. lwork(ij)) THEN
1074          ment(ij, i, j) = ment(ij, i, j)*asij(ij)
1075          bsum(ij, i) = bsum(ij, i) + ment(ij, i, j)
1076        END IF
1077      END DO
1078    END DO
1079    DO ij = 1, ncum
1080      IF ((i>=icb(ij)+1) .AND. (i<=inb(ij)) .AND. (bsum(ij, &
1081          i)<1.0E-18) .AND. lwork(ij)) THEN
1082        nent(ij, i) = 0
1083        ment(ij, i, i) = m(ij, i)
1084        qent(ij, i, i) = q(ij, nk(ij)) - ep(ij, i)*clw(ij, i)
1085        uent(ij, i, i) = u(ij, nk(ij))
1086        vent(ij, i, i) = v(ij, nk(ij))
1087        elij(ij, i, i) = clw(ij, i)
1088        sij(ij, i, i) = 1.0
1089      END IF
1090    END DO
1091789 END DO
1092
1093  RETURN
1094END SUBROUTINE cv_mixing
1095
1096SUBROUTINE cv_unsat(nloc, ncum, nd, inb, t, q, qs, gz, u, v, p, ph, h, lv, &
1097    ep, sigp, clw, m, ment, elij, iflag, mp, qp, up, vp, wt, water, evap)
1098  USE lmdz_cv_ini, ONLY : cl,coeffr,coeffs,cpd,g,ginv,nl,omtrain,omtsnow,sigd
1099
1100  IMPLICIT NONE
1101
1102  ! inputs:
1103  INTEGER ncum, nd, nloc
1104  INTEGER inb(nloc)
1105  REAL t(nloc, nd), q(nloc, nd), qs(nloc, nd)
1106  REAL gz(nloc, nd), u(nloc, nd), v(nloc, nd)
1107  REAL p(nloc, nd), ph(nloc, nd+1), h(nloc, nd)
1108  REAL lv(nloc, nd), ep(nloc, nd), sigp(nloc, nd), clw(nloc, nd)
1109  REAL m(nloc, nd), ment(nloc, nd, nd), elij(nloc, nd, nd)
1110
1111  ! outputs:
1112  INTEGER iflag(nloc) ! also an input
1113  REAL mp(nloc, nd), qp(nloc, nd), up(nloc, nd), vp(nloc, nd)
1114  REAL water(nloc, nd), evap(nloc, nd), wt(nloc, nd)
1115
1116  ! local variables:
1117  INTEGER i, j, k, ij, num1
1118  INTEGER jtt(nloc)
1119  REAL awat, coeff, qsm, afac, sigt, b6, c6, revap
1120  REAL dhdp, fac, qstm, rat
1121  REAL wdtrain(nloc)
1122  LOGICAL lwork(nloc)
1123
1124  ! =====================================================================
1125  ! --- PRECIPITATING DOWNDRAFT CALCULATION
1126  ! =====================================================================
1127
1128  ! Initializations:
1129
1130  DO i = 1, ncum
1131    DO k = 1, nl + 1
1132      wt(i, k) = omtsnow
1133      mp(i, k) = 0.0
1134      evap(i, k) = 0.0
1135      water(i, k) = 0.0
1136    END DO
1137  END DO
1138
1139  DO i = 1, ncum
1140    qp(i, 1) = q(i, 1)
1141    up(i, 1) = u(i, 1)
1142    vp(i, 1) = v(i, 1)
1143  END DO
1144
1145  DO k = 2, nl + 1
1146    DO i = 1, ncum
1147      qp(i, k) = q(i, k-1)
1148      up(i, k) = u(i, k-1)
1149      vp(i, k) = v(i, k-1)
1150    END DO
1151  END DO
1152
1153
1154  ! ***  Check whether ep(inb)=0, if so, skip precipitating    ***
1155  ! ***             downdraft calculation                      ***
1156
1157
1158  ! ***  Integrate liquid water equation to find condensed water   ***
1159  ! ***                and condensed water flux                    ***
1160
1161
1162  DO i = 1, ncum
1163    jtt(i) = 2
1164    IF (ep(i,inb(i))<=0.0001) iflag(i) = 2
1165    IF (iflag(i)==0) THEN
1166      lwork(i) = .TRUE.
1167    ELSE
1168      lwork(i) = .FALSE.
1169    END IF
1170  END DO
1171
1172  ! ***                    Begin downdraft loop                    ***
1173
1174
1175  CALL zilch(wdtrain, ncum)
1176  DO i = nl + 1, 1, -1
1177
1178    num1 = 0
1179    DO ij = 1, ncum
1180      IF ((i<=inb(ij)) .AND. lwork(ij)) num1 = num1 + 1
1181    END DO
1182    IF (num1<=0) GO TO 899
1183
1184
1185    ! ***        Calculate detrained precipitation             ***
1186
1187    DO ij = 1, ncum
1188      IF ((i<=inb(ij)) .AND. (lwork(ij))) THEN
1189        wdtrain(ij) = g*ep(ij, i)*m(ij, i)*clw(ij, i)
1190      END IF
1191    END DO
1192
1193    IF (i>1) THEN
1194      DO j = 1, i - 1
1195        DO ij = 1, ncum
1196          IF ((i<=inb(ij)) .AND. (lwork(ij))) THEN
1197            awat = elij(ij, j, i) - (1.-ep(ij,i))*clw(ij, i)
1198            awat = max(0.0, awat)
1199            wdtrain(ij) = wdtrain(ij) + g*awat*ment(ij, j, i)
1200          END IF
1201        END DO
1202      END DO
1203    END IF
1204
1205    ! ***    Find rain water and evaporation using provisional   ***
1206    ! ***              estimates of qp(i)and qp(i-1)             ***
1207
1208
1209    ! ***  Value of terminal velocity and coeffecient of evaporation for snow
1210    ! ***
1211
1212    DO ij = 1, ncum
1213      IF ((i<=inb(ij)) .AND. (lwork(ij))) THEN
1214        coeff = coeffs
1215        wt(ij, i) = omtsnow
1216
1217        ! ***  Value of terminal velocity and coeffecient of evaporation for
1218        ! rain   ***
1219
1220        IF (t(ij,i)>273.0) THEN
1221          coeff = coeffr
1222          wt(ij, i) = omtrain
1223        END IF
1224        qsm = 0.5*(q(ij,i)+qp(ij,i+1))
1225        afac = coeff*ph(ij, i)*(qs(ij,i)-qsm)/(1.0E4+2.0E3*ph(ij,i)*qs(ij,i))
1226        afac = max(afac, 0.0)
1227        sigt = sigp(ij, i)
1228        sigt = max(0.0, sigt)
1229        sigt = min(1.0, sigt)
1230        b6 = 100.*(ph(ij,i)-ph(ij,i+1))*sigt*afac/wt(ij, i)
1231        c6 = (water(ij,i+1)*wt(ij,i+1)+wdtrain(ij)/sigd)/wt(ij, i)
1232        revap = 0.5*(-b6+sqrt(b6*b6+4.*c6))
1233        evap(ij, i) = sigt*afac*revap
1234        water(ij, i) = revap*revap
1235
1236        ! ***  Calculate precipitating downdraft mass flux under     ***
1237        ! ***              hydrostatic approximation                 ***
1238
1239        IF (i>1) THEN
1240          dhdp = (h(ij,i)-h(ij,i-1))/(p(ij,i-1)-p(ij,i))
1241          dhdp = max(dhdp, 10.0)
1242          mp(ij, i) = 100.*ginv*lv(ij, i)*sigd*evap(ij, i)/dhdp
1243          mp(ij, i) = max(mp(ij,i), 0.0)
1244
1245          ! ***   Add small amount of inertia to downdraft              ***
1246
1247          fac = 20.0/(ph(ij,i-1)-ph(ij,i))
1248          mp(ij, i) = (fac*mp(ij,i+1)+mp(ij,i))/(1.+fac)
1249
1250          ! ***      Force mp to decrease linearly to zero
1251          ! ***
1252          ! ***      between about 950 mb and the surface
1253          ! ***
1254
1255          IF (p(ij,i)>(0.949*p(ij,1))) THEN
1256            jtt(ij) = max(jtt(ij), i)
1257            mp(ij, i) = mp(ij, jtt(ij))*(p(ij,1)-p(ij,i))/ &
1258              (p(ij,1)-p(ij,jtt(ij)))
1259          END IF
1260        END IF
1261
1262        ! ***       Find mixing ratio of precipitating downdraft     ***
1263
1264        IF (i/=inb(ij)) THEN
1265          IF (i==1) THEN
1266            qstm = qs(ij, 1)
1267          ELSE
1268            qstm = qs(ij, i-1)
1269          END IF
1270          IF (mp(ij,i)>mp(ij,i+1)) THEN
1271            rat = mp(ij, i+1)/mp(ij, i)
1272            qp(ij, i) = qp(ij, i+1)*rat + q(ij, i)*(1.0-rat) + &
1273              100.*ginv*sigd*(ph(ij,i)-ph(ij,i+1))*(evap(ij,i)/mp(ij,i))
1274            up(ij, i) = up(ij, i+1)*rat + u(ij, i)*(1.-rat)
1275            vp(ij, i) = vp(ij, i+1)*rat + v(ij, i)*(1.-rat)
1276          ELSE
1277            IF (mp(ij,i+1)>0.0) THEN
1278              qp(ij, i) = (gz(ij,i+1)-gz(ij,i)+qp(ij,i+1)*(lv(ij,i+1)+t(ij, &
1279                i+1)*(cl-cpd))+cpd*(t(ij,i+1)-t(ij, &
1280                i)))/(lv(ij,i)+t(ij,i)*(cl-cpd))
1281              up(ij, i) = up(ij, i+1)
1282              vp(ij, i) = vp(ij, i+1)
1283            END IF
1284          END IF
1285          qp(ij, i) = min(qp(ij,i), qstm)
1286          qp(ij, i) = max(qp(ij,i), 0.0)
1287        END IF
1288      END IF
1289    END DO
1290899 END DO
1291
1292  RETURN
1293END SUBROUTINE cv_unsat
1294
1295SUBROUTINE cv_yield(nloc, ncum, nd, nk, icb, inb, delt, t, q, u, v, gz, p, &
1296    ph, h, hp, lv, cpn, ep, clw, frac, m, mp, qp, up, vp, wt, water, evap, &
1297    ment, qent, uent, vent, nent, elij, tv, tvp, iflag, wd, qprime, tprime, &
1298    precip, cbmf, ft, fq, fu, fv, ma, qcondc)
1299  USE lmdz_cv_ini, ONLY : g,lv0,nl,rrd,sigd,betad,cl,cpd,cpv,cu,delta
1300
1301  IMPLICIT NONE
1302
1303
1304  ! inputs
1305  INTEGER ncum, nd, nloc
1306  INTEGER nk(nloc), icb(nloc), inb(nloc)
1307  INTEGER nent(nloc, nd)
1308  REAL delt
1309  REAL t(nloc, nd), q(nloc, nd), u(nloc, nd), v(nloc, nd)
1310  REAL gz(nloc, nd)
1311  REAL p(nloc, nd), ph(nloc, nd+1), h(nloc, nd)
1312  REAL hp(nloc, nd), lv(nloc, nd)
1313  REAL cpn(nloc, nd), ep(nloc, nd), clw(nloc, nd), frac(nloc)
1314  REAL m(nloc, nd), mp(nloc, nd), qp(nloc, nd)
1315  REAL up(nloc, nd), vp(nloc, nd)
1316  REAL wt(nloc, nd), water(nloc, nd), evap(nloc, nd)
1317  REAL ment(nloc, nd, nd), qent(nloc, nd, nd), elij(nloc, nd, nd)
1318  REAL uent(nloc, nd, nd), vent(nloc, nd, nd)
1319  REAL tv(nloc, nd), tvp(nloc, nd)
1320
1321  ! outputs
1322  INTEGER iflag(nloc) ! also an input
1323  REAL cbmf(nloc) ! also an input
1324  REAL wd(nloc), tprime(nloc), qprime(nloc)
1325  REAL precip(nloc)
1326  REAL ft(nloc, nd), fq(nloc, nd), fu(nloc, nd), fv(nloc, nd)
1327  REAL ma(nloc, nd)
1328  REAL qcondc(nloc, nd)
1329
1330  ! local variables
1331  INTEGER i, j, ij, k, num1
1332  REAL dpinv, cpinv, awat, fqold, ftold, fuold, fvold, delti
1333  REAL work(nloc), am(nloc), amp1(nloc), ad(nloc)
1334  REAL ents(nloc), uav(nloc), vav(nloc), lvcp(nloc, nd)
1335  REAL qcond(nloc, nd), nqcond(nloc, nd), wa(nloc, nd) ! cld
1336  REAL siga(nloc, nd), ax(nloc, nd), mac(nloc, nd) ! cld
1337
1338
1339  ! -- initializations:
1340
1341  delti = 1.0/delt
1342
1343  DO i = 1, ncum
1344    precip(i) = 0.0
1345    wd(i) = 0.0
1346    tprime(i) = 0.0
1347    qprime(i) = 0.0
1348    DO k = 1, nl + 1
1349      ft(i, k) = 0.0
1350      fu(i, k) = 0.0
1351      fv(i, k) = 0.0
1352      fq(i, k) = 0.0
1353      lvcp(i, k) = lv(i, k)/cpn(i, k)
1354      qcondc(i, k) = 0.0 ! cld
1355      qcond(i, k) = 0.0 ! cld
1356      nqcond(i, k) = 0.0 ! cld
1357    END DO
1358  END DO
1359
1360
1361  ! ***  Calculate surface precipitation in mm/day     ***
1362
1363  DO i = 1, ncum
1364    IF (iflag(i)<=1) THEN
1365      ! c            precip(i)=precip(i)+wt(i,1)*sigd*water(i,1)*3600.*24000.
1366      ! c     &                /(rowl*g)
1367      ! c            precip(i)=precip(i)*delt/86400.
1368      precip(i) = wt(i, 1)*sigd*water(i, 1)*86400/g
1369    END IF
1370  END DO
1371
1372
1373  ! ***  Calculate downdraft velocity scale and surface temperature and  ***
1374  ! ***                    water vapor fluctuations                      ***
1375
1376  DO i = 1, ncum
1377    wd(i) = betad*abs(mp(i,icb(i)))*0.01*rrd*t(i, icb(i))/(sigd*p(i,icb(i)))
1378    qprime(i) = 0.5*(qp(i,1)-q(i,1))
1379    tprime(i) = lv0*qprime(i)/cpd
1380  END DO
1381
1382  ! ***  Calculate tendencies of lowest level potential temperature  ***
1383  ! ***                      and mixing ratio                        ***
1384
1385  DO i = 1, ncum
1386    work(i) = 0.01/(ph(i,1)-ph(i,2))
1387    am(i) = 0.0
1388  END DO
1389  DO k = 2, nl
1390    DO i = 1, ncum
1391      IF ((nk(i)==1) .AND. (k<=inb(i)) .AND. (nk(i)==1)) THEN
1392        am(i) = am(i) + m(i, k)
1393      END IF
1394    END DO
1395  END DO
1396  DO i = 1, ncum
1397    IF ((g*work(i)*am(i))>=delti) iflag(i) = 1
1398    ft(i, 1) = ft(i, 1) + g*work(i)*am(i)*(t(i,2)-t(i,1)+(gz(i,2)-gz(i, &
1399      1))/cpn(i,1))
1400    ft(i, 1) = ft(i, 1) - lvcp(i, 1)*sigd*evap(i, 1)
1401    ft(i, 1) = ft(i, 1) + sigd*wt(i, 2)*(cl-cpd)*water(i, 2)*(t(i,2)-t(i,1))* &
1402      work(i)/cpn(i, 1)
1403    fq(i, 1) = fq(i, 1) + g*mp(i, 2)*(qp(i,2)-q(i,1))*work(i) + &
1404      sigd*evap(i, 1)
1405    fq(i, 1) = fq(i, 1) + g*am(i)*(q(i,2)-q(i,1))*work(i)
1406    fu(i, 1) = fu(i, 1) + g*work(i)*(mp(i,2)*(up(i,2)-u(i,1))+am(i)*(u(i, &
1407      2)-u(i,1)))
1408    fv(i, 1) = fv(i, 1) + g*work(i)*(mp(i,2)*(vp(i,2)-v(i,1))+am(i)*(v(i, &
1409      2)-v(i,1)))
1410  END DO
1411  DO j = 2, nl
1412    DO i = 1, ncum
1413      IF (j<=inb(i)) THEN
1414        fq(i, 1) = fq(i, 1) + g*work(i)*ment(i, j, 1)*(qent(i,j,1)-q(i,1))
1415        fu(i, 1) = fu(i, 1) + g*work(i)*ment(i, j, 1)*(uent(i,j,1)-u(i,1))
1416        fv(i, 1) = fv(i, 1) + g*work(i)*ment(i, j, 1)*(vent(i,j,1)-v(i,1))
1417      END IF
1418    END DO
1419  END DO
1420
1421  ! ***  Calculate tendencies of potential temperature and mixing ratio  ***
1422  ! ***               at levels above the lowest level                   ***
1423
1424  ! ***  First find the net saturated updraft and downdraft mass fluxes  ***
1425  ! ***                      through each level                          ***
1426
1427  DO i = 2, nl + 1
1428
1429    num1 = 0
1430    DO ij = 1, ncum
1431      IF (i<=inb(ij)) num1 = num1 + 1
1432    END DO
1433    IF (num1<=0) GO TO 1500
1434
1435    CALL zilch(amp1, ncum)
1436    CALL zilch(ad, ncum)
1437
1438    DO k = i + 1, nl + 1
1439      DO ij = 1, ncum
1440        IF ((i>=nk(ij)) .AND. (i<=inb(ij)) .AND. (k<=(inb(ij)+1))) THEN
1441          amp1(ij) = amp1(ij) + m(ij, k)
1442        END IF
1443      END DO
1444    END DO
1445
1446    DO k = 1, i
1447      DO j = i + 1, nl + 1
1448        DO ij = 1, ncum
1449          IF ((j<=(inb(ij)+1)) .AND. (i<=inb(ij))) THEN
1450            amp1(ij) = amp1(ij) + ment(ij, k, j)
1451          END IF
1452        END DO
1453      END DO
1454    END DO
1455    DO k = 1, i - 1
1456      DO j = i, nl + 1
1457        DO ij = 1, ncum
1458          IF ((i<=inb(ij)) .AND. (j<=inb(ij))) THEN
1459            ad(ij) = ad(ij) + ment(ij, j, k)
1460          END IF
1461        END DO
1462      END DO
1463    END DO
1464
1465    DO ij = 1, ncum
1466      IF (i<=inb(ij)) THEN
1467        dpinv = 0.01/(ph(ij,i)-ph(ij,i+1))
1468        cpinv = 1.0/cpn(ij, i)
1469
1470        ft(ij, i) = ft(ij, i) + g*dpinv*(amp1(ij)*(t(ij,i+1)-t(ij, &
1471          i)+(gz(ij,i+1)-gz(ij,i))*cpinv)-ad(ij)*(t(ij,i)-t(ij, &
1472          i-1)+(gz(ij,i)-gz(ij,i-1))*cpinv)) - sigd*lvcp(ij, i)*evap(ij, i)
1473        ft(ij, i) = ft(ij, i) + g*dpinv*ment(ij, i, i)*(hp(ij,i)-h(ij,i)+t(ij &
1474          ,i)*(cpv-cpd)*(q(ij,i)-qent(ij,i,i)))*cpinv
1475        ft(ij, i) = ft(ij, i) + sigd*wt(ij, i+1)*(cl-cpd)*water(ij, i+1)*(t( &
1476          ij,i+1)-t(ij,i))*dpinv*cpinv
1477        fq(ij, i) = fq(ij, i) + g*dpinv*(amp1(ij)*(q(ij,i+1)-q(ij, &
1478          i))-ad(ij)*(q(ij,i)-q(ij,i-1)))
1479        fu(ij, i) = fu(ij, i) + g*dpinv*(amp1(ij)*(u(ij,i+1)-u(ij, &
1480          i))-ad(ij)*(u(ij,i)-u(ij,i-1)))
1481        fv(ij, i) = fv(ij, i) + g*dpinv*(amp1(ij)*(v(ij,i+1)-v(ij, &
1482          i))-ad(ij)*(v(ij,i)-v(ij,i-1)))
1483      END IF
1484    END DO
1485    DO k = 1, i - 1
1486      DO ij = 1, ncum
1487        IF (i<=inb(ij)) THEN
1488          awat = elij(ij, k, i) - (1.-ep(ij,i))*clw(ij, i)
1489          awat = max(awat, 0.0)
1490          fq(ij, i) = fq(ij, i) + g*dpinv*ment(ij, k, i)*(qent(ij,k,i)-awat-q &
1491            (ij,i))
1492          fu(ij, i) = fu(ij, i) + g*dpinv*ment(ij, k, i)*(uent(ij,k,i)-u(ij,i &
1493            ))
1494          fv(ij, i) = fv(ij, i) + g*dpinv*ment(ij, k, i)*(vent(ij,k,i)-v(ij,i &
1495            ))
1496          ! (saturated updrafts resulting from mixing)               ! cld
1497          qcond(ij, i) = qcond(ij, i) + (elij(ij,k,i)-awat) ! cld
1498          nqcond(ij, i) = nqcond(ij, i) + 1. ! cld
1499        END IF
1500      END DO
1501    END DO
1502    DO k = i, nl + 1
1503      DO ij = 1, ncum
1504        IF ((i<=inb(ij)) .AND. (k<=inb(ij))) THEN
1505          fq(ij, i) = fq(ij, i) + g*dpinv*ment(ij, k, i)*(qent(ij,k,i)-q(ij,i &
1506            ))
1507          fu(ij, i) = fu(ij, i) + g*dpinv*ment(ij, k, i)*(uent(ij,k,i)-u(ij,i &
1508            ))
1509          fv(ij, i) = fv(ij, i) + g*dpinv*ment(ij, k, i)*(vent(ij,k,i)-v(ij,i &
1510            ))
1511        END IF
1512      END DO
1513    END DO
1514    DO ij = 1, ncum
1515      IF (i<=inb(ij)) THEN
1516        fq(ij, i) = fq(ij, i) + sigd*evap(ij, i) + g*(mp(ij,i+1)*(qp(ij, &
1517          i+1)-q(ij,i))-mp(ij,i)*(qp(ij,i)-q(ij,i-1)))*dpinv
1518        fu(ij, i) = fu(ij, i) + g*(mp(ij,i+1)*(up(ij,i+1)-u(ij, &
1519          i))-mp(ij,i)*(up(ij,i)-u(ij,i-1)))*dpinv
1520        fv(ij, i) = fv(ij, i) + g*(mp(ij,i+1)*(vp(ij,i+1)-v(ij, &
1521          i))-mp(ij,i)*(vp(ij,i)-v(ij,i-1)))*dpinv
1522        ! (saturated downdrafts resulting from mixing)               ! cld
1523        DO k = i + 1, inb(ij) ! cld
1524          qcond(ij, i) = qcond(ij, i) + elij(ij, k, i) ! cld
1525          nqcond(ij, i) = nqcond(ij, i) + 1. ! cld
1526        END DO ! cld
1527        ! (particular case: no detraining level is found)            ! cld
1528        IF (nent(ij,i)==0) THEN ! cld
1529          qcond(ij, i) = qcond(ij, i) + (1.-ep(ij,i))*clw(ij, i) ! cld
1530          nqcond(ij, i) = nqcond(ij, i) + 1. ! cld
1531        END IF ! cld
1532        IF (nqcond(ij,i)/=0.) THEN ! cld
1533          qcond(ij, i) = qcond(ij, i)/nqcond(ij, i) ! cld
1534        END IF ! cld
1535      END IF
1536    END DO
15371500 END DO
1538
1539  ! *** Adjust tendencies at top of convection layer to reflect  ***
1540  ! ***       actual position of the level zero cape             ***
1541
1542  DO ij = 1, ncum
1543    fqold = fq(ij, inb(ij))
1544    fq(ij, inb(ij)) = fq(ij, inb(ij))*(1.-frac(ij))
1545    fq(ij, inb(ij)-1) = fq(ij, inb(ij)-1) + frac(ij)*fqold*((ph(ij, &
1546      inb(ij))-ph(ij,inb(ij)+1))/(ph(ij,inb(ij)-1)-ph(ij, &
1547      inb(ij))))*lv(ij, inb(ij))/lv(ij, inb(ij)-1)
1548    ftold = ft(ij, inb(ij))
1549    ft(ij, inb(ij)) = ft(ij, inb(ij))*(1.-frac(ij))
1550    ft(ij, inb(ij)-1) = ft(ij, inb(ij)-1) + frac(ij)*ftold*((ph(ij, &
1551      inb(ij))-ph(ij,inb(ij)+1))/(ph(ij,inb(ij)-1)-ph(ij, &
1552      inb(ij))))*cpn(ij, inb(ij))/cpn(ij, inb(ij)-1)
1553    fuold = fu(ij, inb(ij))
1554    fu(ij, inb(ij)) = fu(ij, inb(ij))*(1.-frac(ij))
1555    fu(ij, inb(ij)-1) = fu(ij, inb(ij)-1) + frac(ij)*fuold*((ph(ij, &
1556      inb(ij))-ph(ij,inb(ij)+1))/(ph(ij,inb(ij)-1)-ph(ij,inb(ij))))
1557    fvold = fv(ij, inb(ij))
1558    fv(ij, inb(ij)) = fv(ij, inb(ij))*(1.-frac(ij))
1559    fv(ij, inb(ij)-1) = fv(ij, inb(ij)-1) + frac(ij)*fvold*((ph(ij, &
1560      inb(ij))-ph(ij,inb(ij)+1))/(ph(ij,inb(ij)-1)-ph(ij,inb(ij))))
1561  END DO
1562
1563  ! ***   Very slightly adjust tendencies to force exact   ***
1564  ! ***     enthalpy, momentum and tracer conservation     ***
1565
1566  DO ij = 1, ncum
1567    ents(ij) = 0.0
1568    uav(ij) = 0.0
1569    vav(ij) = 0.0
1570    DO i = 1, inb(ij)
1571      ents(ij) = ents(ij) + (cpn(ij,i)*ft(ij,i)+lv(ij,i)*fq(ij,i))*(ph(ij,i)- &
1572        ph(ij,i+1))
1573      uav(ij) = uav(ij) + fu(ij, i)*(ph(ij,i)-ph(ij,i+1))
1574      vav(ij) = vav(ij) + fv(ij, i)*(ph(ij,i)-ph(ij,i+1))
1575    END DO
1576  END DO
1577  DO ij = 1, ncum
1578    ents(ij) = ents(ij)/(ph(ij,1)-ph(ij,inb(ij)+1))
1579    uav(ij) = uav(ij)/(ph(ij,1)-ph(ij,inb(ij)+1))
1580    vav(ij) = vav(ij)/(ph(ij,1)-ph(ij,inb(ij)+1))
1581  END DO
1582  DO ij = 1, ncum
1583    DO i = 1, inb(ij)
1584      ft(ij, i) = ft(ij, i) - ents(ij)/cpn(ij, i)
1585      fu(ij, i) = (1.-cu)*(fu(ij,i)-uav(ij))
1586      fv(ij, i) = (1.-cu)*(fv(ij,i)-vav(ij))
1587    END DO
1588  END DO
1589
1590  DO k = 1, nl + 1
1591    DO i = 1, ncum
1592      IF ((q(i,k)+delt*fq(i,k))<0.0) iflag(i) = 10
1593    END DO
1594  END DO
1595
1596
1597  DO i = 1, ncum
1598    IF (iflag(i)>2) THEN
1599      precip(i) = 0.0
1600      cbmf(i) = 0.0
1601    END IF
1602  END DO
1603  DO k = 1, nl
1604    DO i = 1, ncum
1605      IF (iflag(i)>2) THEN
1606        ft(i, k) = 0.0
1607        fq(i, k) = 0.0
1608        fu(i, k) = 0.0
1609        fv(i, k) = 0.0
1610        qcondc(i, k) = 0.0 ! cld
1611      END IF
1612    END DO
1613  END DO
1614
1615  DO k = 1, nl + 1
1616    DO i = 1, ncum
1617      ma(i, k) = 0.
1618    END DO
1619  END DO
1620  DO k = nl, 1, -1
1621    DO i = 1, ncum
1622      ma(i, k) = ma(i, k+1) + m(i, k)
1623    END DO
1624  END DO
1625
1626
1627  ! *** diagnose the in-cloud mixing ratio   ***            ! cld
1628  ! ***           of condensed water         ***            ! cld
1629  ! ! cld
1630  DO ij = 1, ncum ! cld
1631    DO i = 1, nd ! cld
1632      mac(ij, i) = 0.0 ! cld
1633      wa(ij, i) = 0.0 ! cld
1634      siga(ij, i) = 0.0 ! cld
1635    END DO ! cld
1636    DO i = nk(ij), inb(ij) ! cld
1637      DO k = i + 1, inb(ij) + 1 ! cld
1638        mac(ij, i) = mac(ij, i) + m(ij, k) ! cld
1639      END DO ! cld
1640    END DO ! cld
1641    DO i = icb(ij), inb(ij) - 1 ! cld
1642      ax(ij, i) = 0. ! cld
1643      DO j = icb(ij), i ! cld
1644        ax(ij, i) = ax(ij, i) + rrd*(tvp(ij,j)-tv(ij,j)) & ! cld
1645          *(ph(ij,j)-ph(ij,j+1))/p(ij, j) ! cld
1646      END DO ! cld
1647      IF (ax(ij,i)>0.0) THEN ! cld
1648        wa(ij, i) = sqrt(2.*ax(ij,i)) ! cld
1649      END IF ! cld
1650    END DO ! cld
1651    DO i = 1, nl ! cld
1652      IF (wa(ij,i)>0.0) &          ! cld
1653        siga(ij, i) = mac(ij, i)/wa(ij, i) & ! cld
1654        *rrd*tvp(ij, i)/p(ij, i)/100./delta ! cld
1655      siga(ij, i) = min(siga(ij,i), 1.0) ! cld
1656      qcondc(ij, i) = siga(ij, i)*clw(ij, i)*(1.-ep(ij,i)) & ! cld
1657        +(1.-siga(ij,i))*qcond(ij, i) ! cld
1658    END DO ! cld
1659  END DO ! cld
1660
1661  RETURN
1662END SUBROUTINE cv_yield
1663
1664SUBROUTINE cv_uncompress(nloc, len, ncum, nd, idcum, iflag, precip, cbmf, ft, &
1665    fq, fu, fv, ma, qcondc, iflag1, precip1, cbmf1, ft1, fq1, fu1, fv1, ma1, &
1666    qcondc1)
1667   USE lmdz_cv_ini, ONLY : nl
1668    IMPLICIT NONE
1669
1670
1671  ! inputs:
1672  INTEGER len, ncum, nd, nloc
1673  INTEGER idcum(nloc)
1674  INTEGER iflag(nloc)
1675  REAL precip(nloc), cbmf(nloc)
1676  REAL ft(nloc, nd), fq(nloc, nd), fu(nloc, nd), fv(nloc, nd)
1677  REAL ma(nloc, nd)
1678  REAL qcondc(nloc, nd) !cld
1679
1680  ! outputs:
1681  INTEGER iflag1(len)
1682  REAL precip1(len), cbmf1(len)
1683  REAL ft1(len, nd), fq1(len, nd), fu1(len, nd), fv1(len, nd)
1684  REAL ma1(len, nd)
1685  REAL qcondc1(len, nd) !cld
1686
1687  ! local variables:
1688  INTEGER i, k
1689
1690  DO i = 1, ncum
1691    precip1(idcum(i)) = precip(i)
1692    cbmf1(idcum(i)) = cbmf(i)
1693    iflag1(idcum(i)) = iflag(i)
1694  END DO
1695
1696  DO k = 1, nl
1697    DO i = 1, ncum
1698      ft1(idcum(i), k) = ft(i, k)
1699      fq1(idcum(i), k) = fq(i, k)
1700      fu1(idcum(i), k) = fu(i, k)
1701      fv1(idcum(i), k) = fv(i, k)
1702      ma1(idcum(i), k) = ma(i, k)
1703      qcondc1(idcum(i), k) = qcondc(i, k)
1704    END DO
1705  END DO
1706
1707  RETURN
1708END SUBROUTINE cv_uncompress
1709
Note: See TracBrowser for help on using the repository browser.