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

Last change on this file since 5705 was 5705, checked in by yann meurdesoif, 5 weeks ago

Convection GPU porting : suppress potential dependency between columns, may change results (cv_closure)

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: 51.5 KB
Line 
1
2! $Id: cv_routines.f90 5705 2025-06-13 19:47:39Z 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  !ym do not do that, independance between column
253  !ym icbmax = 2
254  !ym DO i = 1, len
255  !ym  icbmax = max(icbmax, icb(i))
256  !ym 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,nl
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  !ym bad dependance between column => icbmax computed in cv_feed
309!ym  DO k = minorig, icbmax - 1
310  DO k = minorig, nd
311    DO i = 1, len
312      IF (k <= MAX(2,icb(i))-1) THEN
313        tp(i, k) = tnk(i) - (gz(i,k)-gznk(i))/cpp(i)
314        tvp(i, k) = tp(i, k)*(1.+qnk(i)*epsi)
315      ENDIF
316    END DO
317  END DO
318
319  ! ***  Find lifted parcel quantities above cloud base    ***
320
321  DO i = 1, len
322    tg = ticb(i)
323    qg = qs(i, icb(i))
324    alv = lv0 - clmcpv*(ticb(i)-t0)
325
326    ! First iteration.
327
328    s = cpd + alv*alv*qg/(rrv*ticb(i)*ticb(i))
329    s = 1./s
330    ahg = cpd*tg + (cl-cpd)*qnk(i)*ticb(i) + alv*qg + gzicb(i)
331    tg = tg + s*(ah0(i)-ahg)
332    tg = max(tg, 35.0)
333    tc = tg - t0
334    denom = 243.5 + tc
335    IF (tc>=0.0) THEN
336      es = 6.112*exp(17.67*tc/denom)
337    ELSE
338      es = exp(23.33086-6111.72784/tg+0.15215*log(tg))
339    END IF
340    qg = eps*es/(p(i,icb(i))-es*(1.-eps))
341
342    ! Second iteration.
343
344    s = cpd + alv*alv*qg/(rrv*ticb(i)*ticb(i))
345    s = 1./s
346    ahg = cpd*tg + (cl-cpd)*qnk(i)*ticb(i) + alv*qg + gzicb(i)
347    tg = tg + s*(ah0(i)-ahg)
348    tg = max(tg, 35.0)
349    tc = tg - t0
350    denom = 243.5 + tc
351    IF (tc>=0.0) THEN
352      es = 6.112*exp(17.67*tc/denom)
353    ELSE
354      es = exp(23.33086-6111.72784/tg+0.15215*log(tg))
355    END IF
356    qg = eps*es/(p(i,icb(i))-es*(1.-eps))
357
358    alv = lv0 - clmcpv*(ticb(i)-273.15)
359    tp(i, icb(i)) = (ah0(i)-(cl-cpd)*qnk(i)*ticb(i)-gz(i,icb(i))-alv*qg)/cpd
360    clw(i, icb(i)) = qnk(i) - qg
361    clw(i, icb(i)) = max(0.0, clw(i,icb(i)))
362    rg = qg/(1.-qnk(i))
363    tvp(i, icb(i)) = tp(i, icb(i))*(1.+rg*epsi)
364  END DO
365 
366  !ym bad dependance between column => ibmax computed in cv_feed
367!ym  DO k = minorig, icbmax
368  DO k = minorig, nd
369    DO i = 1, len
370      IF (k <= MAX(2,icb(i))) THEN
371        tvp(i, k) = tvp(i, k) - tp(i, k)*qnk(i)
372      ENDIF
373    END DO
374  END DO
375
376  RETURN
377END SUBROUTINE cv_undilute1
378
379SUBROUTINE cv_trigger(len, nd, icb, cbmf, tv, tvp, iflag)
380   USE lmdz_cv_ini, ONLY : dtmax
381    IMPLICIT NONE
382
383  ! -------------------------------------------------------------------
384  ! --- Test for instability.
385  ! --- If there was no convection at last time step and parcel
386  ! --- is stable at icb, then set iflag to 4.
387  ! -------------------------------------------------------------------
388
389
390  ! inputs:
391  INTEGER len, nd, icb(len)
392  REAL cbmf(len), tv(len, nd), tvp(len, nd)
393
394  ! outputs:
395  INTEGER iflag(len) ! also an input
396
397  ! local variables:
398  INTEGER i
399
400
401  DO i = 1, len
402    IF ((cbmf(i)==0.0) .AND. (iflag(i)==0) .AND. (tvp(i, &
403      icb(i))<=(tv(i,icb(i))-dtmax))) iflag(i) = 4
404  END DO
405
406  RETURN
407END SUBROUTINE cv_trigger
408
409SUBROUTINE cv_compress(len, nloc, ncum, nd, iflag1, nk1, icb1, cbmf1, plcl1, &
410    tnk1, qnk1, gznk1, t1, q1, qs1, u1, v1, gz1, h1, lv1, cpn1, p1, ph1, tv1, &
411    tp1, tvp1, clw1, iflag, nk, icb, cbmf, plcl, tnk, qnk, gznk, t, q, qs, u, &
412    v, gz, h, lv, cpn, p, ph, tv, tp, tvp, clw, dph)
413   USE lmdz_cv_ini, ONLY : nl
414    USE print_control_mod, ONLY: lunout
415  IMPLICIT NONE
416
417
418  ! inputs:
419  INTEGER len, ncum, nd, nloc
420  INTEGER iflag1(len), nk1(len), icb1(len)
421  REAL cbmf1(len), plcl1(len), tnk1(len), qnk1(len), gznk1(len)
422  REAL t1(len, nd), q1(len, nd), qs1(len, nd), u1(len, nd), v1(len, nd)
423  REAL gz1(len, nd), h1(len, nd), lv1(len, nd), cpn1(len, nd)
424  REAL p1(len, nd), ph1(len, nd+1), tv1(len, nd), tp1(len, nd)
425  REAL tvp1(len, nd), clw1(len, nd)
426
427  ! outputs:
428  INTEGER iflag(nloc), nk(nloc), icb(nloc)
429  REAL cbmf(nloc), plcl(nloc), tnk(nloc), qnk(nloc), gznk(nloc)
430  REAL t(nloc, nd), q(nloc, nd), qs(nloc, nd), u(nloc, nd), v(nloc, nd)
431  REAL gz(nloc, nd), h(nloc, nd), lv(nloc, nd), cpn(nloc, nd)
432  REAL p(nloc, nd), ph(nloc, nd+1), tv(nloc, nd), tp(nloc, nd)
433  REAL tvp(nloc, nd), clw(nloc, nd)
434  REAL dph(nloc, nd)
435
436  ! local variables:
437  INTEGER i, k, nn
438  CHARACTER (LEN=20) :: modname = 'cv_compress'
439  CHARACTER (LEN=80) :: abort_message
440
441
442  DO k = 1, nl + 1
443    nn = 0
444    DO i = 1, len
445      IF (iflag1(i)==0) THEN
446        nn = nn + 1
447        t(nn, k) = t1(i, k)
448        q(nn, k) = q1(i, k)
449        qs(nn, k) = qs1(i, k)
450        u(nn, k) = u1(i, k)
451        v(nn, k) = v1(i, k)
452        gz(nn, k) = gz1(i, k)
453        h(nn, k) = h1(i, k)
454        lv(nn, k) = lv1(i, k)
455        cpn(nn, k) = cpn1(i, k)
456        p(nn, k) = p1(i, k)
457        ph(nn, k) = ph1(i, k)
458        tv(nn, k) = tv1(i, k)
459        tp(nn, k) = tp1(i, k)
460        tvp(nn, k) = tvp1(i, k)
461        clw(nn, k) = clw1(i, k)
462      END IF
463    END DO
464  END DO
465
466  IF (nn/=ncum) THEN
467    WRITE (lunout, *) 'strange! nn not equal to ncum: ', nn, ncum
468    abort_message = ''
469    CALL abort_physic(modname, abort_message, 1)
470  END IF
471
472  nn = 0
473  DO i = 1, len
474    IF (iflag1(i)==0) THEN
475      nn = nn + 1
476      cbmf(nn) = cbmf1(i)
477      plcl(nn) = plcl1(i)
478      tnk(nn) = tnk1(i)
479      qnk(nn) = qnk1(i)
480      gznk(nn) = gznk1(i)
481      nk(nn) = nk1(i)
482      icb(nn) = icb1(i)
483      iflag(nn) = iflag1(i)
484    END IF
485  END DO
486
487  DO k = 1, nl
488    DO i = 1, ncum
489      dph(i, k) = ph(i, k) - ph(i, k+1)
490    END DO
491  END DO
492
493  RETURN
494END SUBROUTINE cv_compress
495
496SUBROUTINE cv_undilute2(nloc, ncum, nd, icb, nk, tnk, qnk, gznk, t, q, qs, &
497    gz, p, dph, h, tv, lv, inb, inb1, tp, tvp, clw, hp, ep, sigp, frac)
498  USE lmdz_cv_ini, ONLY : nl,cl,clmcpv,cpd,cpv,elcrit,eps,epsi,lv0,minorig,nlp,rrv,sigs,t0,tlcrit
499
500  IMPLICIT NONE
501
502  ! ---------------------------------------------------------------------
503  ! Purpose:
504  ! FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
505  ! &
506  ! COMPUTE THE PRECIPITATION EFFICIENCIES AND THE
507  ! FRACTION OF PRECIPITATION FALLING OUTSIDE OF CLOUD
508  ! &
509  ! FIND THE LEVEL OF NEUTRAL BUOYANCY
510  ! ---------------------------------------------------------------------
511
512  ! inputs:
513  INTEGER ncum, nd, nloc
514  INTEGER icb(nloc), nk(nloc)
515  REAL t(nloc, nd), q(nloc, nd), qs(nloc, nd), gz(nloc, nd)
516  REAL p(nloc, nd), dph(nloc, nd)
517  REAL tnk(nloc), qnk(nloc), gznk(nloc)
518  REAL lv(nloc, nd), tv(nloc, nd), h(nloc, nd)
519
520  ! outputs:
521  INTEGER inb(nloc), inb1(nloc)
522  REAL tp(nloc, nd), tvp(nloc, nd), clw(nloc, nd)
523  REAL ep(nloc, nd), sigp(nloc, nd), hp(nloc, nd)
524  REAL frac(nloc)
525
526  ! local variables:
527  INTEGER i, k
528  REAL tg, qg, ahg, alv, s, tc, es, denom, rg, tca, elacrit
529  REAL by, defrac
530  REAL ah0(nloc), cape(nloc), capem(nloc), byp(nloc)
531  LOGICAL lcape(nloc)
532
533  ! =====================================================================
534  ! --- SOME INITIALIZATIONS
535  ! =====================================================================
536
537  DO k = 1, nl
538    DO i = 1, ncum
539      ep(i, k) = 0.0
540      sigp(i, k) = sigs
541    END DO
542  END DO
543
544  ! =====================================================================
545  ! --- FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
546  ! =====================================================================
547
548  ! ---       The procedure is to solve the equation.
549  ! cp*tp+L*qp+phi=cp*tnk+L*qnk+gznk.
550
551  ! ***  Calculate certain parcel quantities, including static energy   ***
552
553
554  DO i = 1, ncum
555    ah0(i) = (cpd*(1.-qnk(i))+cl*qnk(i))*tnk(i) + qnk(i)*(lv0-clmcpv*(tnk(i)- &
556      t0)) + gznk(i)
557  END DO
558
559
560  ! ***  Find lifted parcel quantities above cloud base    ***
561
562
563  DO k = minorig + 1, nl
564    DO i = 1, ncum
565      IF (k>=(icb(i)+1)) THEN
566        tg = t(i, k)
567        qg = qs(i, k)
568        alv = lv0 - clmcpv*(t(i,k)-t0)
569
570        ! First iteration.
571
572        s = cpd + alv*alv*qg/(rrv*t(i,k)*t(i,k))
573        s = 1./s
574        ahg = cpd*tg + (cl-cpd)*qnk(i)*t(i, k) + alv*qg + gz(i, k)
575        tg = tg + s*(ah0(i)-ahg)
576        tg = max(tg, 35.0)
577        tc = tg - t0
578        denom = 243.5 + tc
579        IF (tc>=0.0) THEN
580          es = 6.112*exp(17.67*tc/denom)
581        ELSE
582          es = exp(23.33086-6111.72784/tg+0.15215*log(tg))
583        END IF
584        qg = eps*es/(p(i,k)-es*(1.-eps))
585
586        ! Second iteration.
587
588        s = cpd + alv*alv*qg/(rrv*t(i,k)*t(i,k))
589        s = 1./s
590        ahg = cpd*tg + (cl-cpd)*qnk(i)*t(i, k) + alv*qg + gz(i, k)
591        tg = tg + s*(ah0(i)-ahg)
592        tg = max(tg, 35.0)
593        tc = tg - t0
594        denom = 243.5 + tc
595        IF (tc>=0.0) THEN
596          es = 6.112*exp(17.67*tc/denom)
597        ELSE
598          es = exp(23.33086-6111.72784/tg+0.15215*log(tg))
599        END IF
600        qg = eps*es/(p(i,k)-es*(1.-eps))
601
602        alv = lv0 - clmcpv*(t(i,k)-t0)
603        ! print*,'cpd dans convect2 ',cpd
604        ! print*,'tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd'
605        ! print*,tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd
606        tp(i, k) = (ah0(i)-(cl-cpd)*qnk(i)*t(i,k)-gz(i,k)-alv*qg)/cpd
607        ! if (.not.cpd.gt.1000.) then
608        ! print*,'CPD=',cpd
609        ! stop
610        ! endif
611        clw(i, k) = qnk(i) - qg
612        clw(i, k) = max(0.0, clw(i,k))
613        rg = qg/(1.-qnk(i))
614        tvp(i, k) = tp(i, k)*(1.+rg*epsi)
615      END IF
616    END DO
617  END DO
618
619  ! =====================================================================
620  ! --- SET THE PRECIPITATION EFFICIENCIES AND THE FRACTION OF
621  ! --- PRECIPITATION FALLING OUTSIDE OF CLOUD
622  ! --- THESE MAY BE FUNCTIONS OF TP(I), P(I) AND CLW(I)
623  ! =====================================================================
624
625  DO k = minorig + 1, nl
626    DO i = 1, ncum
627      IF (k>=(nk(i)+1)) THEN
628        tca = tp(i, k) - t0
629        IF (tca>=0.0) THEN
630          elacrit = elcrit
631        ELSE
632          elacrit = elcrit*(1.0-tca/tlcrit)
633        END IF
634        elacrit = max(elacrit, 0.0)
635        ep(i, k) = 1.0 - elacrit/max(clw(i,k), 1.0E-8)
636        ep(i, k) = max(ep(i,k), 0.0)
637        ep(i, k) = min(ep(i,k), 1.0)
638        sigp(i, k) = sigs
639      END IF
640    END DO
641  END DO
642
643  ! =====================================================================
644  ! --- CALCULATE VIRTUAL TEMPERATURE AND LIFTED PARCEL
645  ! --- VIRTUAL TEMPERATURE
646  ! =====================================================================
647
648  DO k = minorig + 1, nl
649    DO i = 1, ncum
650      IF (k>=(icb(i)+1)) THEN
651        tvp(i, k) = tvp(i, k)*(1.0-qnk(i)+ep(i,k)*clw(i,k))
652        ! print*,'i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)'
653        ! print*, i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)
654      END IF
655    END DO
656  END DO
657  DO i = 1, ncum
658    tvp(i, nlp) = tvp(i, nl) - (gz(i,nlp)-gz(i,nl))/cpd
659  END DO
660
661  ! =====================================================================
662  ! --- FIND THE FIRST MODEL LEVEL (INB1) ABOVE THE PARCEL'S
663  ! --- HIGHEST LEVEL OF NEUTRAL BUOYANCY
664  ! --- AND THE HIGHEST LEVEL OF POSITIVE CAPE (INB)
665  ! =====================================================================
666
667  DO i = 1, ncum
668    cape(i) = 0.0
669    capem(i) = 0.0
670    inb(i) = icb(i) + 1
671    inb1(i) = inb(i)
672  END DO
673
674  ! Originial Code
675
676  ! do 530 k=minorig+1,nl-1
677  ! do 520 i=1,ncum
678  ! if(k.ge.(icb(i)+1))then
679  ! by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
680  ! byp=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
681  ! cape(i)=cape(i)+by
682  ! if(by.ge.0.0)inb1(i)=k+1
683  ! if(cape(i).gt.0.0)then
684  ! inb(i)=k+1
685  ! capem(i)=cape(i)
686  ! endif
687  ! endif
688  ! 520    continue
689  ! 530  continue
690  ! do 540 i=1,ncum
691  ! byp=(tvp(i,nl)-tv(i,nl))*dph(i,nl)/p(i,nl)
692  ! cape(i)=capem(i)+byp
693  ! defrac=capem(i)-cape(i)
694  ! defrac=max(defrac,0.001)
695  ! frac(i)=-cape(i)/defrac
696  ! frac(i)=min(frac(i),1.0)
697  ! frac(i)=max(frac(i),0.0)
698  ! 540   continue
699
700  ! K Emanuel fix
701
702  ! call zilch(byp,ncum)
703  ! do 530 k=minorig+1,nl-1
704  ! do 520 i=1,ncum
705  ! if(k.ge.(icb(i)+1))then
706  ! by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k)
707  ! cape(i)=cape(i)+by
708  ! if(by.ge.0.0)inb1(i)=k+1
709  ! if(cape(i).gt.0.0)then
710  ! inb(i)=k+1
711  ! capem(i)=cape(i)
712  ! byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1)
713  ! endif
714  ! endif
715  ! 520    continue
716  ! 530  continue
717  ! do 540 i=1,ncum
718  ! inb(i)=max(inb(i),inb1(i))
719  ! cape(i)=capem(i)+byp(i)
720  ! defrac=capem(i)-cape(i)
721  ! defrac=max(defrac,0.001)
722  ! frac(i)=-cape(i)/defrac
723  ! frac(i)=min(frac(i),1.0)
724  ! frac(i)=max(frac(i),0.0)
725  ! 540   continue
726
727  ! J Teixeira fix
728
729  byp(1:ncum) = 0
730
731  DO i = 1, ncum
732    lcape(i) = .TRUE.
733  END DO
734  DO k = minorig + 1, nl - 1
735    DO i = 1, ncum
736      IF (cape(i)<0.0) lcape(i) = .FALSE.
737      IF ((k>=(icb(i)+1)) .AND. lcape(i)) THEN
738        by = (tvp(i,k)-tv(i,k))*dph(i, k)/p(i, k)
739        byp(i) = (tvp(i,k+1)-tv(i,k+1))*dph(i, k+1)/p(i, k+1)
740        cape(i) = cape(i) + by
741        IF (by>=0.0) inb1(i) = k + 1
742        IF (cape(i)>0.0) THEN
743          inb(i) = k + 1
744          capem(i) = cape(i)
745        END IF
746      END IF
747    END DO
748  END DO
749  DO i = 1, ncum
750    cape(i) = capem(i) + byp(i)
751    defrac = capem(i) - cape(i)
752    defrac = max(defrac, 0.001)
753    frac(i) = -cape(i)/defrac
754    frac(i) = min(frac(i), 1.0)
755    frac(i) = max(frac(i), 0.0)
756  END DO
757
758  ! =====================================================================
759  ! ---   CALCULATE LIQUID WATER STATIC ENERGY OF LIFTED PARCEL
760  ! =====================================================================
761
762  ! initialization:
763!ym very bad 
764!ym  DO i = 1, ncum*nlp
765!ym    hp(i, 1) = h(i, 1)
766!ym  END DO
767  DO k=1,nlp
768    DO i=1,ncum
769      hp(i, k) = h(i, k)
770    ENDDO
771  ENDDO
772
773  DO k = minorig + 1, nl
774    DO i = 1, ncum
775      IF ((k>=icb(i)) .AND. (k<=inb(i))) THEN
776        hp(i, k) = h(i, nk(i)) + (lv(i,k)+(cpd-cpv)*t(i,k))*ep(i, k)*clw(i, k &
777          )
778      END IF
779    END DO
780  END DO
781
782  RETURN
783END SUBROUTINE cv_undilute2
784
785SUBROUTINE cv_closure(nloc, ncum, nd, nk, icb, tv, tvp, p, ph, dph, plcl, &
786    cpn, iflag, cbmf)
787  USE lmdz_cv_ini, ONLY : alpha,damp,dtmax,minorig,rrd
788
789  IMPLICIT NONE
790
791  ! inputs:
792  INTEGER ncum, nd, nloc
793  INTEGER nk(nloc), icb(nloc)
794  REAL tv(nloc, nd), tvp(nloc, nd), p(nloc, nd), dph(nloc, nd)
795  REAL ph(nloc, nd+1) ! caution nd instead ndp1 to be consistent...
796  REAL plcl(nloc), cpn(nloc, nd)
797
798  ! outputs:
799  INTEGER iflag(nloc)
800  REAL cbmf(nloc) ! also an input
801
802  ! local variables:
803  INTEGER i, k, icbmax
804  REAL dtpbl(nloc), dtmin(nloc), tvpplcl(nloc), tvaplcl(nloc)
805  REAL work(nloc)
806
807  ! -------------------------------------------------------------------
808  ! Compute icbmax.
809  ! -------------------------------------------------------------------
810 
811!ym independance betwwen column
812!ym  icbmax = 2
813!ym  DO i = 1, ncum
814!ym    icbmax = max(icbmax, icb(i))
815!ym  END DO
816
817  ! =====================================================================
818  ! ---  CALCULATE CLOUD BASE MASS FLUX
819  ! =====================================================================
820
821  ! tvpplcl = parcel temperature lifted adiabatically from level
822  ! icb-1 to the LCL.
823  ! tvaplcl = virtual temperature at the LCL.
824
825  DO i = 1, ncum
826    dtpbl(i) = 0.0
827    tvpplcl(i) = tvp(i, icb(i)-1) - rrd*tvp(i, icb(i)-1)*(p(i,icb(i)-1)-plcl( &
828      i))/(cpn(i,icb(i)-1)*p(i,icb(i)-1))
829    tvaplcl(i) = tv(i, icb(i)) + (tvp(i,icb(i))-tvp(i,icb(i)+1))*(plcl(i)-p(i &
830      ,icb(i)))/(p(i,icb(i))-p(i,icb(i)+1))
831  END DO
832
833  ! -------------------------------------------------------------------
834  ! --- Interpolate difference between lifted parcel and
835  ! --- environmental temperatures to lifted condensation level
836  ! -------------------------------------------------------------------
837
838  ! dtpbl = average of tvp-tv in the PBL (k=nk to icb-1).
839!ym independance betwwen column
840!ym  DO k = minorig, icbmax
841  DO k = minorig, nd
842    DO i = 1, ncum
843      IF (k<=MAX(2,icb(i))) THEN
844        IF ((k>=nk(i)) .AND. (k<=(icb(i)-1))) THEN
845          dtpbl(i) = dtpbl(i) + (tvp(i,k)-tv(i,k))*dph(i, k)
846        END IF
847      ENDIF
848    END DO
849  END DO
850  DO i = 1, ncum
851    dtpbl(i) = dtpbl(i)/(ph(i,nk(i))-ph(i,icb(i)))
852    dtmin(i) = tvpplcl(i) - tvaplcl(i) + dtmax + dtpbl(i)
853  END DO
854
855  ! -------------------------------------------------------------------
856  ! --- Adjust cloud base mass flux
857  ! -------------------------------------------------------------------
858
859  DO i = 1, ncum
860    work(i) = cbmf(i)
861    cbmf(i) = max(0.0, (1.0-damp)*cbmf(i)+0.1*alpha*dtmin(i))
862    IF ((work(i)==0.0) .AND. (cbmf(i)==0.0)) THEN
863      iflag(i) = 3
864    END IF
865  END DO
866
867  RETURN
868END SUBROUTINE cv_closure
869
870SUBROUTINE cv_mixing(nloc, ncum, nd, icb, nk, inb, inb1, ph, t, q, qs, u, v, &
871    h, lv, qnk, hp, tv, tvp, ep, clw, cbmf, m, ment, qent, uent, vent, nent, &
872    sij, elij)
873  USE lmdz_cv_ini, ONLY : cpd,cpv,entp,minorig,nl,nlp,rrv
874
875  IMPLICIT NONE
876
877
878  ! inputs:
879  INTEGER ncum, nd, nloc
880  INTEGER icb(nloc), inb(nloc), inb1(nloc), nk(nloc)
881  REAL cbmf(nloc), qnk(nloc)
882  REAL ph(nloc, nd+1)
883  REAL t(nloc, nd), q(nloc, nd), qs(nloc, nd), lv(nloc, nd)
884  REAL u(nloc, nd), v(nloc, nd), h(nloc, nd), hp(nloc, nd)
885  REAL tv(nloc, nd), tvp(nloc, nd), ep(nloc, nd), clw(nloc, nd)
886
887  ! outputs:
888  INTEGER nent(nloc, nd)
889  REAL m(nloc, nd), ment(nloc, nd, nd), qent(nloc, nd, nd)
890  REAL uent(nloc, nd, nd), vent(nloc, nd, nd)
891  REAL sij(nloc, nd, nd), elij(nloc, nd, nd)
892
893  ! local variables:
894  INTEGER i, j, k, ij
895  INTEGER num1, num2
896  REAL dbo, qti, bf2, anum, denom, dei, altem, cwat, stemp
897  REAL alt, qp1, smid, sjmin, sjmax, delp, delm
898  REAL work(nloc), asij(nloc), smin(nloc), scrit(nloc)
899  REAL bsum(nloc, nd)
900  LOGICAL lwork(nloc)
901
902  ! =====================================================================
903  ! --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS
904  ! =====================================================================
905
906!ym very bad
907!ym  DO i = 1, ncum*nlp
908!ym    nent(i, 1) = 0
909!ym    m(i, 1) = 0.0
910!ym  END DO
911  DO k = 1, nlp
912    DO i = 1, ncum
913      nent(i, k) = 0
914      m(i, k) = 0.0
915    ENDDO
916  ENDDO
917
918  DO k = 1, nlp
919    DO j = 1, nlp
920      DO i = 1, ncum
921        qent(i, k, j) = q(i, j)
922        uent(i, k, j) = u(i, j)
923        vent(i, k, j) = v(i, j)
924        elij(i, k, j) = 0.0
925        ment(i, k, j) = 0.0
926        sij(i, k, j) = 0.0
927      END DO
928    END DO
929  END DO
930
931  ! -------------------------------------------------------------------
932  ! --- Calculate rates of mixing,  m(i)
933  ! -------------------------------------------------------------------
934
935  work(1:ncum) = 0.
936 
937  DO j = minorig + 1, nl
938    DO i = 1, ncum
939      IF ((j>=(icb(i)+1)) .AND. (j<=inb(i))) THEN
940        k = min(j, inb1(i))
941        dbo = abs(tv(i,k+1)-tvp(i,k+1)-tv(i,k-1)+tvp(i,k-1)) + &
942          entp*0.04*(ph(i,k)-ph(i,k+1))
943        work(i) = work(i) + dbo
944        m(i, j) = cbmf(i)*dbo
945      END IF
946    END DO
947  END DO
948  DO k = minorig + 1, nl
949    DO i = 1, ncum
950      IF ((k>=(icb(i)+1)) .AND. (k<=inb(i))) THEN
951        m(i, k) = m(i, k)/work(i)
952      END IF
953    END DO
954  END DO
955
956
957  ! =====================================================================
958  ! --- CALCULATE ENTRAINED AIR MASS FLUX (ment), TOTAL WATER MIXING
959  ! --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING
960  ! --- FRACTION (sij)
961  ! =====================================================================
962
963
964  DO i = minorig + 1, nl
965    DO j = minorig + 1, nl
966      DO ij = 1, ncum
967        IF ((i>=(icb(ij)+1)) .AND. (j>=icb(ij)) .AND. (i<=inb(ij)) .AND. (j<= &
968            inb(ij))) THEN
969          qti = qnk(ij) - ep(ij, i)*clw(ij, i)
970          bf2 = 1. + lv(ij, j)*lv(ij, j)*qs(ij, j)/(rrv*t(ij,j)*t(ij,j)*cpd)
971          anum = h(ij, j) - hp(ij, i) + (cpv-cpd)*t(ij, j)*(qti-q(ij,j))
972          denom = h(ij, i) - hp(ij, i) + (cpd-cpv)*(q(ij,i)-qti)*t(ij, j)
973          dei = denom
974          IF (abs(dei)<0.01) dei = 0.01
975          sij(ij, i, j) = anum/dei
976          sij(ij, i, i) = 1.0
977          altem = sij(ij, i, j)*q(ij, i) + (1.-sij(ij,i,j))*qti - qs(ij, j)
978          altem = altem/bf2
979          cwat = clw(ij, j)*(1.-ep(ij,j))
980          stemp = sij(ij, i, j)
981          IF ((stemp<0.0 .OR. stemp>1.0 .OR. altem>cwat) .AND. j>i) THEN
982            anum = anum - lv(ij, j)*(qti-qs(ij,j)-cwat*bf2)
983            denom = denom + lv(ij, j)*(q(ij,i)-qti)
984            IF (abs(denom)<0.01) denom = 0.01
985            sij(ij, i, j) = anum/denom
986            altem = sij(ij, i, j)*q(ij, i) + (1.-sij(ij,i,j))*qti - qs(ij, j)
987            altem = altem - (bf2-1.)*cwat
988          END IF
989          IF (sij(ij,i,j)>0.0 .AND. sij(ij,i,j)<0.9) THEN
990            qent(ij, i, j) = sij(ij, i, j)*q(ij, i) + (1.-sij(ij,i,j))*qti
991            uent(ij, i, j) = sij(ij, i, j)*u(ij, i) + &
992              (1.-sij(ij,i,j))*u(ij, nk(ij))
993            vent(ij, i, j) = sij(ij, i, j)*v(ij, i) + &
994              (1.-sij(ij,i,j))*v(ij, nk(ij))
995            elij(ij, i, j) = altem
996            elij(ij, i, j) = max(0.0, elij(ij,i,j))
997            ment(ij, i, j) = m(ij, i)/(1.-sij(ij,i,j))
998            nent(ij, i) = nent(ij, i) + 1
999          END IF
1000          sij(ij, i, j) = max(0.0, sij(ij,i,j))
1001          sij(ij, i, j) = min(1.0, sij(ij,i,j))
1002        END IF
1003      END DO
1004    END DO
1005
1006    ! ***   If no air can entrain at level i assume that updraft detrains
1007    ! ***
1008    ! ***   at that level and calculate detrained air flux and properties
1009    ! ***
1010
1011    DO ij = 1, ncum
1012      IF ((i>=(icb(ij)+1)) .AND. (i<=inb(ij)) .AND. (nent(ij,i)==0)) THEN
1013        ment(ij, i, i) = m(ij, i)
1014        qent(ij, i, i) = q(ij, nk(ij)) - ep(ij, i)*clw(ij, i)
1015        uent(ij, i, i) = u(ij, nk(ij))
1016        vent(ij, i, i) = v(ij, nk(ij))
1017        elij(ij, i, i) = clw(ij, i)
1018        sij(ij, i, i) = 1.0
1019      END IF
1020    END DO
1021  END DO
1022
1023  DO i = 1, ncum
1024    sij(i, inb(i), inb(i)) = 1.0
1025  END DO
1026
1027  ! =====================================================================
1028  ! ---  NORMALIZE ENTRAINED AIR MASS FLUXES
1029  ! ---  TO REPRESENT EQUAL PROBABILITIES OF MIXING
1030  ! =====================================================================
1031
1032  bsum(1:ncum,1:nlp) = 0.
1033  DO ij = 1, ncum
1034    lwork(ij) = .FALSE.
1035  END DO
1036  DO i = minorig + 1, nl
1037
1038    num1 = 0
1039    DO ij = 1, ncum
1040      IF ((i>=icb(ij)+1) .AND. (i<=inb(ij))) num1 = num1 + 1
1041    END DO
1042!ym    IF (num1<=0) GO TO 789
1043    IF (num1<=0) CYCLE
1044
1045    DO ij = 1, ncum
1046      IF ((i>=icb(ij)+1) .AND. (i<=inb(ij))) THEN
1047        lwork(ij) = (nent(ij,i)/=0)
1048        qp1 = q(ij, nk(ij)) - ep(ij, i)*clw(ij, i)
1049        anum = h(ij, i) - hp(ij, i) - lv(ij, i)*(qp1-qs(ij,i))
1050        denom = h(ij, i) - hp(ij, i) + lv(ij, i)*(q(ij,i)-qp1)
1051        IF (abs(denom)<0.01) denom = 0.01
1052        scrit(ij) = anum/denom
1053        alt = qp1 - qs(ij, i) + scrit(ij)*(q(ij,i)-qp1)
1054        IF (scrit(ij)<0.0 .OR. alt<0.0) scrit(ij) = 1.0
1055        asij(ij) = 0.0
1056        smin(ij) = 1.0
1057      END IF
1058    END DO
1059    DO j = minorig, nl
1060
1061      num2 = 0
1062      DO ij = 1, ncum
1063        IF ((i>=icb(ij)+1) .AND. (i<=inb(ij)) .AND. (j>=icb( &
1064          ij)) .AND. (j<=inb(ij)) .AND. lwork(ij)) num2 = num2 + 1
1065      END DO
1066!ym      IF (num2<=0) GO TO 783
1067      IF (num2<=0) CYCLE
1068
1069      DO ij = 1, ncum
1070        IF ((i>=icb(ij)+1) .AND. (i<=inb(ij)) .AND. (j>=icb( &
1071            ij)) .AND. (j<=inb(ij)) .AND. lwork(ij)) THEN
1072          IF (sij(ij,i,j)>0.0 .AND. sij(ij,i,j)<0.9) THEN
1073            IF (j>i) THEN
1074              smid = min(sij(ij,i,j), scrit(ij))
1075              sjmax = smid
1076              sjmin = smid
1077              IF (smid<smin(ij) .AND. sij(ij,i,j+1)<smid) THEN
1078                smin(ij) = smid
1079                sjmax = min(sij(ij,i,j+1), sij(ij,i,j), scrit(ij))
1080                sjmin = max(sij(ij,i,j-1), sij(ij,i,j))
1081                sjmin = min(sjmin, scrit(ij))
1082              END IF
1083            ELSE
1084              sjmax = max(sij(ij,i,j+1), scrit(ij))
1085              smid = max(sij(ij,i,j), scrit(ij))
1086              sjmin = 0.0
1087              IF (j>1) sjmin = sij(ij, i, j-1)
1088              sjmin = max(sjmin, scrit(ij))
1089            END IF
1090            delp = abs(sjmax-smid)
1091            delm = abs(sjmin-smid)
1092            asij(ij) = asij(ij) + (delp+delm)*(ph(ij,j)-ph(ij,j+1))
1093            ment(ij, i, j) = ment(ij, i, j)*(delp+delm)*(ph(ij,j)-ph(ij,j+1))
1094          END IF
1095        END IF
1096      END DO
1097783 END DO
1098    DO ij = 1, ncum
1099      IF ((i>=icb(ij)+1) .AND. (i<=inb(ij)) .AND. lwork(ij)) THEN
1100        asij(ij) = max(1.0E-21, asij(ij))
1101        asij(ij) = 1.0/asij(ij)
1102        bsum(ij, i) = 0.0
1103      END IF
1104    END DO
1105    DO j = minorig, nl + 1
1106      DO ij = 1, ncum
1107        IF ((i>=icb(ij)+1) .AND. (i<=inb(ij)) .AND. (j>=icb( &
1108            ij)) .AND. (j<=inb(ij)) .AND. lwork(ij)) THEN
1109          ment(ij, i, j) = ment(ij, i, j)*asij(ij)
1110          bsum(ij, i) = bsum(ij, i) + ment(ij, i, j)
1111        END IF
1112      END DO
1113    END DO
1114    DO ij = 1, ncum
1115      IF ((i>=icb(ij)+1) .AND. (i<=inb(ij)) .AND. (bsum(ij, &
1116          i)<1.0E-18) .AND. lwork(ij)) THEN
1117        nent(ij, i) = 0
1118        ment(ij, i, i) = m(ij, i)
1119        qent(ij, i, i) = q(ij, nk(ij)) - ep(ij, i)*clw(ij, i)
1120        uent(ij, i, i) = u(ij, nk(ij))
1121        vent(ij, i, i) = v(ij, nk(ij))
1122        elij(ij, i, i) = clw(ij, i)
1123        sij(ij, i, i) = 1.0
1124      END IF
1125    END DO
1126789 END DO
1127
1128  RETURN
1129END SUBROUTINE cv_mixing
1130
1131SUBROUTINE cv_unsat(nloc, ncum, nd, inb, t, q, qs, gz, u, v, p, ph, h, lv, &
1132    ep, sigp, clw, m, ment, elij, iflag, mp, qp, up, vp, wt, water, evap)
1133  USE lmdz_cv_ini, ONLY : cl,coeffr,coeffs,cpd,g,ginv,nl,omtrain,omtsnow,sigd
1134
1135  IMPLICIT NONE
1136
1137  ! inputs:
1138  INTEGER ncum, nd, nloc
1139  INTEGER inb(nloc)
1140  REAL t(nloc, nd), q(nloc, nd), qs(nloc, nd)
1141  REAL gz(nloc, nd), u(nloc, nd), v(nloc, nd)
1142  REAL p(nloc, nd), ph(nloc, nd+1), h(nloc, nd)
1143  REAL lv(nloc, nd), ep(nloc, nd), sigp(nloc, nd), clw(nloc, nd)
1144  REAL m(nloc, nd), ment(nloc, nd, nd), elij(nloc, nd, nd)
1145
1146  ! outputs:
1147  INTEGER iflag(nloc) ! also an input
1148  REAL mp(nloc, nd), qp(nloc, nd), up(nloc, nd), vp(nloc, nd)
1149  REAL water(nloc, nd), evap(nloc, nd), wt(nloc, nd)
1150
1151  ! local variables:
1152  INTEGER i, j, k, ij, num1
1153  INTEGER jtt(nloc)
1154  REAL awat, coeff, qsm, afac, sigt, b6, c6, revap
1155  REAL dhdp, fac, qstm, rat
1156  REAL wdtrain(nloc)
1157  LOGICAL lwork(nloc)
1158
1159  ! =====================================================================
1160  ! --- PRECIPITATING DOWNDRAFT CALCULATION
1161  ! =====================================================================
1162
1163  ! Initializations:
1164
1165  DO i = 1, ncum
1166    DO k = 1, nl + 1
1167      wt(i, k) = omtsnow
1168      mp(i, k) = 0.0
1169      evap(i, k) = 0.0
1170      water(i, k) = 0.0
1171    END DO
1172  END DO
1173
1174  DO i = 1, ncum
1175    qp(i, 1) = q(i, 1)
1176    up(i, 1) = u(i, 1)
1177    vp(i, 1) = v(i, 1)
1178  END DO
1179
1180  DO k = 2, nl + 1
1181    DO i = 1, ncum
1182      qp(i, k) = q(i, k-1)
1183      up(i, k) = u(i, k-1)
1184      vp(i, k) = v(i, k-1)
1185    END DO
1186  END DO
1187
1188
1189  ! ***  Check whether ep(inb)=0, if so, skip precipitating    ***
1190  ! ***             downdraft calculation                      ***
1191
1192
1193  ! ***  Integrate liquid water equation to find condensed water   ***
1194  ! ***                and condensed water flux                    ***
1195
1196
1197  DO i = 1, ncum
1198    jtt(i) = 2
1199    IF (ep(i,inb(i))<=0.0001) iflag(i) = 2
1200    IF (iflag(i)==0) THEN
1201      lwork(i) = .TRUE.
1202    ELSE
1203      lwork(i) = .FALSE.
1204    END IF
1205  END DO
1206
1207  ! ***                    Begin downdraft loop                    ***
1208
1209
1210  wdtrain(1:ncum) = 0.
1211  DO i = nl + 1, 1, -1
1212
1213    num1 = 0
1214    DO ij = 1, ncum
1215      IF ((i<=inb(ij)) .AND. lwork(ij)) num1 = num1 + 1
1216    END DO
1217!ym    IF (num1<=0) GO TO 899
1218    IF (num1<=0) CYCLE
1219
1220
1221    ! ***        Calculate detrained precipitation             ***
1222
1223    DO ij = 1, ncum
1224      IF ((i<=inb(ij)) .AND. (lwork(ij))) THEN
1225        wdtrain(ij) = g*ep(ij, i)*m(ij, i)*clw(ij, i)
1226      END IF
1227    END DO
1228
1229    IF (i>1) THEN
1230      DO j = 1, i - 1
1231        DO ij = 1, ncum
1232          IF ((i<=inb(ij)) .AND. (lwork(ij))) THEN
1233            awat = elij(ij, j, i) - (1.-ep(ij,i))*clw(ij, i)
1234            awat = max(0.0, awat)
1235            wdtrain(ij) = wdtrain(ij) + g*awat*ment(ij, j, i)
1236          END IF
1237        END DO
1238      END DO
1239    END IF
1240
1241    ! ***    Find rain water and evaporation using provisional   ***
1242    ! ***              estimates of qp(i)and qp(i-1)             ***
1243
1244
1245    ! ***  Value of terminal velocity and coeffecient of evaporation for snow
1246    ! ***
1247
1248    DO ij = 1, ncum
1249      IF ((i<=inb(ij)) .AND. (lwork(ij))) THEN
1250        coeff = coeffs
1251        wt(ij, i) = omtsnow
1252
1253        ! ***  Value of terminal velocity and coeffecient of evaporation for
1254        ! rain   ***
1255
1256        IF (t(ij,i)>273.0) THEN
1257          coeff = coeffr
1258          wt(ij, i) = omtrain
1259        END IF
1260        qsm = 0.5*(q(ij,i)+qp(ij,i+1))
1261        afac = coeff*ph(ij, i)*(qs(ij,i)-qsm)/(1.0E4+2.0E3*ph(ij,i)*qs(ij,i))
1262        afac = max(afac, 0.0)
1263        sigt = sigp(ij, i)
1264        sigt = max(0.0, sigt)
1265        sigt = min(1.0, sigt)
1266        b6 = 100.*(ph(ij,i)-ph(ij,i+1))*sigt*afac/wt(ij, i)
1267        c6 = (water(ij,i+1)*wt(ij,i+1)+wdtrain(ij)/sigd)/wt(ij, i)
1268        revap = 0.5*(-b6+sqrt(b6*b6+4.*c6))
1269        evap(ij, i) = sigt*afac*revap
1270        water(ij, i) = revap*revap
1271
1272        ! ***  Calculate precipitating downdraft mass flux under     ***
1273        ! ***              hydrostatic approximation                 ***
1274
1275        IF (i>1) THEN
1276          dhdp = (h(ij,i)-h(ij,i-1))/(p(ij,i-1)-p(ij,i))
1277          dhdp = max(dhdp, 10.0)
1278          mp(ij, i) = 100.*ginv*lv(ij, i)*sigd*evap(ij, i)/dhdp
1279          mp(ij, i) = max(mp(ij,i), 0.0)
1280
1281          ! ***   Add small amount of inertia to downdraft              ***
1282
1283          fac = 20.0/(ph(ij,i-1)-ph(ij,i))
1284          mp(ij, i) = (fac*mp(ij,i+1)+mp(ij,i))/(1.+fac)
1285
1286          ! ***      Force mp to decrease linearly to zero
1287          ! ***
1288          ! ***      between about 950 mb and the surface
1289          ! ***
1290
1291          IF (p(ij,i)>(0.949*p(ij,1))) THEN
1292            jtt(ij) = max(jtt(ij), i)
1293            mp(ij, i) = mp(ij, jtt(ij))*(p(ij,1)-p(ij,i))/ &
1294              (p(ij,1)-p(ij,jtt(ij)))
1295          END IF
1296        END IF
1297
1298        ! ***       Find mixing ratio of precipitating downdraft     ***
1299
1300        IF (i/=inb(ij)) THEN
1301          IF (i==1) THEN
1302            qstm = qs(ij, 1)
1303          ELSE
1304            qstm = qs(ij, i-1)
1305          END IF
1306          IF (mp(ij,i)>mp(ij,i+1)) THEN
1307            rat = mp(ij, i+1)/mp(ij, i)
1308            qp(ij, i) = qp(ij, i+1)*rat + q(ij, i)*(1.0-rat) + &
1309              100.*ginv*sigd*(ph(ij,i)-ph(ij,i+1))*(evap(ij,i)/mp(ij,i))
1310            up(ij, i) = up(ij, i+1)*rat + u(ij, i)*(1.-rat)
1311            vp(ij, i) = vp(ij, i+1)*rat + v(ij, i)*(1.-rat)
1312          ELSE
1313            IF (mp(ij,i+1)>0.0) THEN
1314              qp(ij, i) = (gz(ij,i+1)-gz(ij,i)+qp(ij,i+1)*(lv(ij,i+1)+t(ij, &
1315                i+1)*(cl-cpd))+cpd*(t(ij,i+1)-t(ij, &
1316                i)))/(lv(ij,i)+t(ij,i)*(cl-cpd))
1317              up(ij, i) = up(ij, i+1)
1318              vp(ij, i) = vp(ij, i+1)
1319            END IF
1320          END IF
1321          qp(ij, i) = min(qp(ij,i), qstm)
1322          qp(ij, i) = max(qp(ij,i), 0.0)
1323        END IF
1324      END IF
1325    END DO
1326899 END DO
1327
1328  RETURN
1329END SUBROUTINE cv_unsat
1330
1331SUBROUTINE cv_yield(nloc, ncum, nd, nk, icb, inb, delt, t, q, u, v, gz, p, &
1332    ph, h, hp, lv, cpn, ep, clw, frac, m, mp, qp, up, vp, wt, water, evap, &
1333    ment, qent, uent, vent, nent, elij, tv, tvp, iflag, wd, qprime, tprime, &
1334    precip, cbmf, ft, fq, fu, fv, ma, qcondc)
1335  USE lmdz_cv_ini, ONLY : g,lv0,nl,rrd,sigd,betad,cl,cpd,cpv,cu,delta
1336
1337  IMPLICIT NONE
1338
1339
1340  ! inputs
1341  INTEGER ncum, nd, nloc
1342  INTEGER nk(nloc), icb(nloc), inb(nloc)
1343  INTEGER nent(nloc, nd)
1344  REAL delt
1345  REAL t(nloc, nd), q(nloc, nd), u(nloc, nd), v(nloc, nd)
1346  REAL gz(nloc, nd)
1347  REAL p(nloc, nd), ph(nloc, nd+1), h(nloc, nd)
1348  REAL hp(nloc, nd), lv(nloc, nd)
1349  REAL cpn(nloc, nd), ep(nloc, nd), clw(nloc, nd), frac(nloc)
1350  REAL m(nloc, nd), mp(nloc, nd), qp(nloc, nd)
1351  REAL up(nloc, nd), vp(nloc, nd)
1352  REAL wt(nloc, nd), water(nloc, nd), evap(nloc, nd)
1353  REAL ment(nloc, nd, nd), qent(nloc, nd, nd), elij(nloc, nd, nd)
1354  REAL uent(nloc, nd, nd), vent(nloc, nd, nd)
1355  REAL tv(nloc, nd), tvp(nloc, nd)
1356
1357  ! outputs
1358  INTEGER iflag(nloc) ! also an input
1359  REAL cbmf(nloc) ! also an input
1360  REAL wd(nloc), tprime(nloc), qprime(nloc)
1361  REAL precip(nloc)
1362  REAL ft(nloc, nd), fq(nloc, nd), fu(nloc, nd), fv(nloc, nd)
1363  REAL ma(nloc, nd)
1364  REAL qcondc(nloc, nd)
1365
1366  ! local variables
1367  INTEGER i, j, ij, k, num1
1368  REAL dpinv, cpinv, awat, fqold, ftold, fuold, fvold, delti
1369  REAL work(nloc), am(nloc), amp1(nloc), ad(nloc)
1370  REAL ents(nloc), uav(nloc), vav(nloc), lvcp(nloc, nd)
1371  REAL qcond(nloc, nd), nqcond(nloc, nd), wa(nloc, nd) ! cld
1372  REAL siga(nloc, nd), ax(nloc, nd), mac(nloc, nd) ! cld
1373
1374
1375  ! -- initializations:
1376
1377  delti = 1.0/delt
1378
1379  DO i = 1, ncum
1380    precip(i) = 0.0
1381    wd(i) = 0.0
1382    tprime(i) = 0.0
1383    qprime(i) = 0.0
1384    DO k = 1, nl + 1
1385      ft(i, k) = 0.0
1386      fu(i, k) = 0.0
1387      fv(i, k) = 0.0
1388      fq(i, k) = 0.0
1389      lvcp(i, k) = lv(i, k)/cpn(i, k)
1390      qcondc(i, k) = 0.0 ! cld
1391      qcond(i, k) = 0.0 ! cld
1392      nqcond(i, k) = 0.0 ! cld
1393    END DO
1394  END DO
1395
1396
1397  ! ***  Calculate surface precipitation in mm/day     ***
1398
1399  DO i = 1, ncum
1400    IF (iflag(i)<=1) THEN
1401      ! c            precip(i)=precip(i)+wt(i,1)*sigd*water(i,1)*3600.*24000.
1402      ! c     &                /(rowl*g)
1403      ! c            precip(i)=precip(i)*delt/86400.
1404      precip(i) = wt(i, 1)*sigd*water(i, 1)*86400/g
1405    END IF
1406  END DO
1407
1408
1409  ! ***  Calculate downdraft velocity scale and surface temperature and  ***
1410  ! ***                    water vapor fluctuations                      ***
1411
1412  DO i = 1, ncum
1413    wd(i) = betad*abs(mp(i,icb(i)))*0.01*rrd*t(i, icb(i))/(sigd*p(i,icb(i)))
1414    qprime(i) = 0.5*(qp(i,1)-q(i,1))
1415    tprime(i) = lv0*qprime(i)/cpd
1416  END DO
1417
1418  ! ***  Calculate tendencies of lowest level potential temperature  ***
1419  ! ***                      and mixing ratio                        ***
1420
1421  DO i = 1, ncum
1422    work(i) = 0.01/(ph(i,1)-ph(i,2))
1423    am(i) = 0.0
1424  END DO
1425  DO k = 2, nl
1426    DO i = 1, ncum
1427      IF ((nk(i)==1) .AND. (k<=inb(i)) .AND. (nk(i)==1)) THEN
1428        am(i) = am(i) + m(i, k)
1429      END IF
1430    END DO
1431  END DO
1432  DO i = 1, ncum
1433    IF ((g*work(i)*am(i))>=delti) iflag(i) = 1
1434    ft(i, 1) = ft(i, 1) + g*work(i)*am(i)*(t(i,2)-t(i,1)+(gz(i,2)-gz(i, &
1435      1))/cpn(i,1))
1436    ft(i, 1) = ft(i, 1) - lvcp(i, 1)*sigd*evap(i, 1)
1437    ft(i, 1) = ft(i, 1) + sigd*wt(i, 2)*(cl-cpd)*water(i, 2)*(t(i,2)-t(i,1))* &
1438      work(i)/cpn(i, 1)
1439    fq(i, 1) = fq(i, 1) + g*mp(i, 2)*(qp(i,2)-q(i,1))*work(i) + &
1440      sigd*evap(i, 1)
1441    fq(i, 1) = fq(i, 1) + g*am(i)*(q(i,2)-q(i,1))*work(i)
1442    fu(i, 1) = fu(i, 1) + g*work(i)*(mp(i,2)*(up(i,2)-u(i,1))+am(i)*(u(i, &
1443      2)-u(i,1)))
1444    fv(i, 1) = fv(i, 1) + g*work(i)*(mp(i,2)*(vp(i,2)-v(i,1))+am(i)*(v(i, &
1445      2)-v(i,1)))
1446  END DO
1447  DO j = 2, nl
1448    DO i = 1, ncum
1449      IF (j<=inb(i)) THEN
1450        fq(i, 1) = fq(i, 1) + g*work(i)*ment(i, j, 1)*(qent(i,j,1)-q(i,1))
1451        fu(i, 1) = fu(i, 1) + g*work(i)*ment(i, j, 1)*(uent(i,j,1)-u(i,1))
1452        fv(i, 1) = fv(i, 1) + g*work(i)*ment(i, j, 1)*(vent(i,j,1)-v(i,1))
1453      END IF
1454    END DO
1455  END DO
1456
1457  ! ***  Calculate tendencies of potential temperature and mixing ratio  ***
1458  ! ***               at levels above the lowest level                   ***
1459
1460  ! ***  First find the net saturated updraft and downdraft mass fluxes  ***
1461  ! ***                      through each level                          ***
1462
1463  DO i = 2, nl + 1
1464
1465    num1 = 0
1466    DO ij = 1, ncum
1467      IF (i<=inb(ij)) num1 = num1 + 1
1468    END DO
1469!ym    IF (num1<=0) GO TO 1500
1470    IF (num1<=0) CYCLE
1471
1472    amp1(1:ncum)=0.
1473    ad(1:ncum)=0.
1474
1475    DO k = i + 1, nl + 1
1476      DO ij = 1, ncum
1477        IF ((i>=nk(ij)) .AND. (i<=inb(ij)) .AND. (k<=(inb(ij)+1))) THEN
1478          amp1(ij) = amp1(ij) + m(ij, k)
1479        END IF
1480      END DO
1481    END DO
1482
1483    DO k = 1, i
1484      DO j = i + 1, nl + 1
1485        DO ij = 1, ncum
1486          IF ((j<=(inb(ij)+1)) .AND. (i<=inb(ij))) THEN
1487            amp1(ij) = amp1(ij) + ment(ij, k, j)
1488          END IF
1489        END DO
1490      END DO
1491    END DO
1492    DO k = 1, i - 1
1493      DO j = i, nl + 1
1494        DO ij = 1, ncum
1495          IF ((i<=inb(ij)) .AND. (j<=inb(ij))) THEN
1496            ad(ij) = ad(ij) + ment(ij, j, k)
1497          END IF
1498        END DO
1499      END DO
1500    END DO
1501
1502    DO ij = 1, ncum
1503      IF (i<=inb(ij)) THEN
1504        dpinv = 0.01/(ph(ij,i)-ph(ij,i+1))
1505        cpinv = 1.0/cpn(ij, i)
1506
1507        ft(ij, i) = ft(ij, i) + g*dpinv*(amp1(ij)*(t(ij,i+1)-t(ij, &
1508          i)+(gz(ij,i+1)-gz(ij,i))*cpinv)-ad(ij)*(t(ij,i)-t(ij, &
1509          i-1)+(gz(ij,i)-gz(ij,i-1))*cpinv)) - sigd*lvcp(ij, i)*evap(ij, i)
1510        ft(ij, i) = ft(ij, i) + g*dpinv*ment(ij, i, i)*(hp(ij,i)-h(ij,i)+t(ij &
1511          ,i)*(cpv-cpd)*(q(ij,i)-qent(ij,i,i)))*cpinv
1512        ft(ij, i) = ft(ij, i) + sigd*wt(ij, i+1)*(cl-cpd)*water(ij, i+1)*(t( &
1513          ij,i+1)-t(ij,i))*dpinv*cpinv
1514        fq(ij, i) = fq(ij, i) + g*dpinv*(amp1(ij)*(q(ij,i+1)-q(ij, &
1515          i))-ad(ij)*(q(ij,i)-q(ij,i-1)))
1516        fu(ij, i) = fu(ij, i) + g*dpinv*(amp1(ij)*(u(ij,i+1)-u(ij, &
1517          i))-ad(ij)*(u(ij,i)-u(ij,i-1)))
1518        fv(ij, i) = fv(ij, i) + g*dpinv*(amp1(ij)*(v(ij,i+1)-v(ij, &
1519          i))-ad(ij)*(v(ij,i)-v(ij,i-1)))
1520      END IF
1521    END DO
1522    DO k = 1, i - 1
1523      DO ij = 1, ncum
1524        IF (i<=inb(ij)) THEN
1525          awat = elij(ij, k, i) - (1.-ep(ij,i))*clw(ij, i)
1526          awat = max(awat, 0.0)
1527          fq(ij, i) = fq(ij, i) + g*dpinv*ment(ij, k, i)*(qent(ij,k,i)-awat-q &
1528            (ij,i))
1529          fu(ij, i) = fu(ij, i) + g*dpinv*ment(ij, k, i)*(uent(ij,k,i)-u(ij,i &
1530            ))
1531          fv(ij, i) = fv(ij, i) + g*dpinv*ment(ij, k, i)*(vent(ij,k,i)-v(ij,i &
1532            ))
1533          ! (saturated updrafts resulting from mixing)               ! cld
1534          qcond(ij, i) = qcond(ij, i) + (elij(ij,k,i)-awat) ! cld
1535          nqcond(ij, i) = nqcond(ij, i) + 1. ! cld
1536        END IF
1537      END DO
1538    END DO
1539    DO k = i, nl + 1
1540      DO ij = 1, ncum
1541        IF ((i<=inb(ij)) .AND. (k<=inb(ij))) THEN
1542          fq(ij, i) = fq(ij, i) + g*dpinv*ment(ij, k, i)*(qent(ij,k,i)-q(ij,i &
1543            ))
1544          fu(ij, i) = fu(ij, i) + g*dpinv*ment(ij, k, i)*(uent(ij,k,i)-u(ij,i &
1545            ))
1546          fv(ij, i) = fv(ij, i) + g*dpinv*ment(ij, k, i)*(vent(ij,k,i)-v(ij,i &
1547            ))
1548        END IF
1549      END DO
1550    END DO
1551    DO ij = 1, ncum
1552      IF (i<=inb(ij)) THEN
1553        fq(ij, i) = fq(ij, i) + sigd*evap(ij, i) + g*(mp(ij,i+1)*(qp(ij, &
1554          i+1)-q(ij,i))-mp(ij,i)*(qp(ij,i)-q(ij,i-1)))*dpinv
1555        fu(ij, i) = fu(ij, i) + g*(mp(ij,i+1)*(up(ij,i+1)-u(ij, &
1556          i))-mp(ij,i)*(up(ij,i)-u(ij,i-1)))*dpinv
1557        fv(ij, i) = fv(ij, i) + g*(mp(ij,i+1)*(vp(ij,i+1)-v(ij, &
1558          i))-mp(ij,i)*(vp(ij,i)-v(ij,i-1)))*dpinv
1559        ! (saturated downdrafts resulting from mixing)               ! cld
1560        DO k = i + 1, inb(ij) ! cld
1561          qcond(ij, i) = qcond(ij, i) + elij(ij, k, i) ! cld
1562          nqcond(ij, i) = nqcond(ij, i) + 1. ! cld
1563        END DO ! cld
1564        ! (particular case: no detraining level is found)            ! cld
1565        IF (nent(ij,i)==0) THEN ! cld
1566          qcond(ij, i) = qcond(ij, i) + (1.-ep(ij,i))*clw(ij, i) ! cld
1567          nqcond(ij, i) = nqcond(ij, i) + 1. ! cld
1568        END IF ! cld
1569        IF (nqcond(ij,i)/=0.) THEN ! cld
1570          qcond(ij, i) = qcond(ij, i)/nqcond(ij, i) ! cld
1571        END IF ! cld
1572      END IF
1573    END DO
15741500 END DO
1575
1576  ! *** Adjust tendencies at top of convection layer to reflect  ***
1577  ! ***       actual position of the level zero cape             ***
1578
1579  DO ij = 1, ncum
1580    fqold = fq(ij, inb(ij))
1581    fq(ij, inb(ij)) = fq(ij, inb(ij))*(1.-frac(ij))
1582    fq(ij, inb(ij)-1) = fq(ij, inb(ij)-1) + frac(ij)*fqold*((ph(ij, &
1583      inb(ij))-ph(ij,inb(ij)+1))/(ph(ij,inb(ij)-1)-ph(ij, &
1584      inb(ij))))*lv(ij, inb(ij))/lv(ij, inb(ij)-1)
1585    ftold = ft(ij, inb(ij))
1586    ft(ij, inb(ij)) = ft(ij, inb(ij))*(1.-frac(ij))
1587    ft(ij, inb(ij)-1) = ft(ij, inb(ij)-1) + frac(ij)*ftold*((ph(ij, &
1588      inb(ij))-ph(ij,inb(ij)+1))/(ph(ij,inb(ij)-1)-ph(ij, &
1589      inb(ij))))*cpn(ij, inb(ij))/cpn(ij, inb(ij)-1)
1590    fuold = fu(ij, inb(ij))
1591    fu(ij, inb(ij)) = fu(ij, inb(ij))*(1.-frac(ij))
1592    fu(ij, inb(ij)-1) = fu(ij, inb(ij)-1) + frac(ij)*fuold*((ph(ij, &
1593      inb(ij))-ph(ij,inb(ij)+1))/(ph(ij,inb(ij)-1)-ph(ij,inb(ij))))
1594    fvold = fv(ij, inb(ij))
1595    fv(ij, inb(ij)) = fv(ij, inb(ij))*(1.-frac(ij))
1596    fv(ij, inb(ij)-1) = fv(ij, inb(ij)-1) + frac(ij)*fvold*((ph(ij, &
1597      inb(ij))-ph(ij,inb(ij)+1))/(ph(ij,inb(ij)-1)-ph(ij,inb(ij))))
1598  END DO
1599
1600  ! ***   Very slightly adjust tendencies to force exact   ***
1601  ! ***     enthalpy, momentum and tracer conservation     ***
1602
1603  DO ij = 1, ncum
1604    ents(ij) = 0.0
1605    uav(ij) = 0.0
1606    vav(ij) = 0.0
1607    DO i = 1, inb(ij)
1608      ents(ij) = ents(ij) + (cpn(ij,i)*ft(ij,i)+lv(ij,i)*fq(ij,i))*(ph(ij,i)- &
1609        ph(ij,i+1))
1610      uav(ij) = uav(ij) + fu(ij, i)*(ph(ij,i)-ph(ij,i+1))
1611      vav(ij) = vav(ij) + fv(ij, i)*(ph(ij,i)-ph(ij,i+1))
1612    END DO
1613  END DO
1614  DO ij = 1, ncum
1615    ents(ij) = ents(ij)/(ph(ij,1)-ph(ij,inb(ij)+1))
1616    uav(ij) = uav(ij)/(ph(ij,1)-ph(ij,inb(ij)+1))
1617    vav(ij) = vav(ij)/(ph(ij,1)-ph(ij,inb(ij)+1))
1618  END DO
1619  DO ij = 1, ncum
1620    DO i = 1, inb(ij)
1621      ft(ij, i) = ft(ij, i) - ents(ij)/cpn(ij, i)
1622      fu(ij, i) = (1.-cu)*(fu(ij,i)-uav(ij))
1623      fv(ij, i) = (1.-cu)*(fv(ij,i)-vav(ij))
1624    END DO
1625  END DO
1626
1627  DO k = 1, nl + 1
1628    DO i = 1, ncum
1629      IF ((q(i,k)+delt*fq(i,k))<0.0) iflag(i) = 10
1630    END DO
1631  END DO
1632
1633
1634  DO i = 1, ncum
1635    IF (iflag(i)>2) THEN
1636      precip(i) = 0.0
1637      cbmf(i) = 0.0
1638    END IF
1639  END DO
1640  DO k = 1, nl
1641    DO i = 1, ncum
1642      IF (iflag(i)>2) THEN
1643        ft(i, k) = 0.0
1644        fq(i, k) = 0.0
1645        fu(i, k) = 0.0
1646        fv(i, k) = 0.0
1647        qcondc(i, k) = 0.0 ! cld
1648      END IF
1649    END DO
1650  END DO
1651
1652  DO k = 1, nl + 1
1653    DO i = 1, ncum
1654      ma(i, k) = 0.
1655    END DO
1656  END DO
1657  DO k = nl, 1, -1
1658    DO i = 1, ncum
1659      ma(i, k) = ma(i, k+1) + m(i, k)
1660    END DO
1661  END DO
1662
1663
1664  ! *** diagnose the in-cloud mixing ratio   ***            ! cld
1665  ! ***           of condensed water         ***            ! cld
1666  ! ! cld
1667  DO ij = 1, ncum ! cld
1668    DO i = 1, nd ! cld
1669      mac(ij, i) = 0.0 ! cld
1670      wa(ij, i) = 0.0 ! cld
1671      siga(ij, i) = 0.0 ! cld
1672    END DO ! cld
1673    DO i = nk(ij), inb(ij) ! cld
1674      DO k = i + 1, inb(ij) + 1 ! cld
1675        mac(ij, i) = mac(ij, i) + m(ij, k) ! cld
1676      END DO ! cld
1677    END DO ! cld
1678    DO i = icb(ij), inb(ij) - 1 ! cld
1679      ax(ij, i) = 0. ! cld
1680      DO j = icb(ij), i ! cld
1681        ax(ij, i) = ax(ij, i) + rrd*(tvp(ij,j)-tv(ij,j)) & ! cld
1682          *(ph(ij,j)-ph(ij,j+1))/p(ij, j) ! cld
1683      END DO ! cld
1684      IF (ax(ij,i)>0.0) THEN ! cld
1685        wa(ij, i) = sqrt(2.*ax(ij,i)) ! cld
1686      END IF ! cld
1687    END DO ! cld
1688    DO i = 1, nl ! cld
1689      IF (wa(ij,i)>0.0) &          ! cld
1690        siga(ij, i) = mac(ij, i)/wa(ij, i) & ! cld
1691        *rrd*tvp(ij, i)/p(ij, i)/100./delta ! cld
1692      siga(ij, i) = min(siga(ij,i), 1.0) ! cld
1693      qcondc(ij, i) = siga(ij, i)*clw(ij, i)*(1.-ep(ij,i)) & ! cld
1694        +(1.-siga(ij,i))*qcond(ij, i) ! cld
1695    END DO ! cld
1696  END DO ! cld
1697
1698  RETURN
1699END SUBROUTINE cv_yield
1700
1701SUBROUTINE cv_uncompress(nloc, len, ncum, nd, idcum, iflag, precip, cbmf, ft, &
1702    fq, fu, fv, ma, qcondc, iflag1, precip1, cbmf1, ft1, fq1, fu1, fv1, ma1, &
1703    qcondc1)
1704   USE lmdz_cv_ini, ONLY : nl
1705    IMPLICIT NONE
1706
1707
1708  ! inputs:
1709  INTEGER len, ncum, nd, nloc
1710  INTEGER idcum(nloc)
1711  INTEGER iflag(nloc)
1712  REAL precip(nloc), cbmf(nloc)
1713  REAL ft(nloc, nd), fq(nloc, nd), fu(nloc, nd), fv(nloc, nd)
1714  REAL ma(nloc, nd)
1715  REAL qcondc(nloc, nd) !cld
1716
1717  ! outputs:
1718  INTEGER iflag1(len)
1719  REAL precip1(len), cbmf1(len)
1720  REAL ft1(len, nd), fq1(len, nd), fu1(len, nd), fv1(len, nd)
1721  REAL ma1(len, nd)
1722  REAL qcondc1(len, nd) !cld
1723
1724  ! local variables:
1725  INTEGER i, k
1726
1727  DO i = 1, ncum
1728    precip1(idcum(i)) = precip(i)
1729    cbmf1(idcum(i)) = cbmf(i)
1730    iflag1(idcum(i)) = iflag(i)
1731  END DO
1732
1733  DO k = 1, nl
1734    DO i = 1, ncum
1735      ft1(idcum(i), k) = ft(i, k)
1736      fq1(idcum(i), k) = fq(i, k)
1737      fu1(idcum(i), k) = fu(i, k)
1738      fv1(idcum(i), k) = fv(i, k)
1739      ma1(idcum(i), k) = ma(i, k)
1740      qcondc1(idcum(i), k) = qcondc(i, k)
1741    END DO
1742  END DO
1743
1744  RETURN
1745END SUBROUTINE cv_uncompress
1746
1747END MODULE cv_routines_mod
Note: See TracBrowser for help on using the repository browser.