source: LMDZ6/trunk/libf/phylmd/cv_routines.f90 @ 5700

Last change on this file since 5700 was 5700, checked in by yann meurdesoif, 4 weeks ago

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