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

Last change on this file since 5744 was 5712, checked in by yann meurdesoif, 3 weeks ago

Convection GPU porting : Compression of active convection point is now optional (default remain to true). For GPU runs, convection is not compressed and is computed on each column. The update is done only for column where convection is active

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