source: LMDZ5/trunk/libf/phylmd/convect1.F90 @ 2250

Last change on this file since 2250 was 1992, checked in by lguez, 11 years ago

Converted to free source form files in libf/phylmd which were still in
fixed source form. The conversion was done using the polish mode of
the NAG Fortran Compiler.

In addition to converting to free source form, the processing of the
files also:

-- indented the code (including comments);

-- set Fortran keywords to uppercase, and set all other identifiers
to lower case;

-- added qualifiers to end statements (for example "end subroutine
conflx", instead of "end");

-- changed the terminating statements of all DO loops so that each
loop ends with an ENDDO statement (instead of a labeled continue).

-- replaced #include by include.

  • 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: 19.8 KB
Line 
1
2! $Header$
3
4SUBROUTINE convect1(len, nd, ndp1, noff, minorig, t, q, qs, u, v, p, ph, &
5    iflag, ft, fq, fu, fv, precip, cbmf, delt, ma)
6  ! .............................START PROLOGUE............................
7
8  ! SCCS IDENTIFICATION:  @(#)convect1.f        1.1 04/21/00
9  ! 19:40:52 /h/cm/library/nogaps4/src/sub/fcst/convect1.f_v
10
11  ! CONFIGURATION IDENTIFICATION:  None
12
13  ! MODULE NAME:  convect1
14
15  ! DESCRIPTION:
16
17  ! convect1     The Emanuel Cumulus Convection Scheme
18
19  ! CONTRACT NUMBER AND TITLE:  None
20
21  ! REFERENCES: Programmers  K. Emanuel (MIT), Timothy F. Hogan, M. Peng
22  ! (NRL)
23
24  ! CLASSIFICATION:  Unclassified
25
26  ! RESTRICTIONS: None
27
28  ! COMPILER DEPENDENCIES: FORTRAN 77, FORTRAN 90
29
30  ! COMPILE OPTIONS: Fortran 77: -Zu -Wf"-ei -o aggress"
31  ! Fortran 90: -O vector3,scalar3,task1,aggress,overindex  -ei -r 2
32
33  ! LIBRARIES OF RESIDENCE: /a/ops/lib/libfcst159.a
34
35  ! USAGE: call convect1(len,nd,noff,minorig,
36  ! &                   t,q,qs,u,v,
37  ! &                   p,ph,iflag,ft,
38  ! &                   fq,fu,fv,precip,cbmf,delt)
39
40  ! PARAMETERS:
41  ! Name            Type         Usage            Description
42  ! ----------      ----------     -------  ----------------------------
43
44  ! len           Integer        Input        first (i) dimension
45  ! nd            Integer        Input        vertical (k) dimension
46  ! ndp1          Integer        Input        nd + 1
47  ! noff          Integer        Input        integer limit for convection
48  ! (nd-noff)
49  ! minorig       Integer        Input        First level of convection
50  ! t             Real           Input        temperature
51  ! q             Real           Input        specific hum
52  ! qs            Real           Input        sat specific hum
53  ! u             Real           Input        u-wind
54  ! v             Real           Input        v-wind
55  ! p             Real           Input        full level pressure
56  ! ph            Real           Input        half level pressure
57  ! iflag         Integer        Output       iflag on latitude strip
58  ! ft            Real           Output       temp tend
59  ! fq            Real           Output       spec hum tend
60  ! fu            Real           Output       u-wind tend
61  ! fv            Real           Output       v-wind tend
62  ! cbmf          Real           In/Out       cumulus mass flux
63  ! delt          Real           Input        time step
64  ! iflag         Integer        Output       integer flag for Emanuel
65  ! conditions
66
67  ! COMMON BLOCKS:
68  ! Block      Name     Type    Usage              Notes
69  ! --------  --------   ----    ------   ------------------------
70
71  ! FILES: None
72
73  ! DATA BASES: None
74
75  ! NON-FILE INPUT/OUTPUT: None
76
77  ! ERROR CONDITIONS: None
78
79  ! ADDITIONAL COMMENTS: None
80
81  ! .................MAINTENANCE SECTION................................
82
83  ! MODULES CALLED:
84  ! Name           Description
85  ! convect2        Emanuel cumulus convection tendency calculations
86  ! -------     ----------------------
87  ! LOCAL VARIABLES AND
88  ! STRUCTURES:
89  ! Name     Type    Description
90  ! -------  ------  -----------
91  ! See Comments Below
92
93  ! i        Integer loop index
94  ! k        Integer loop index
95
96  ! METHOD:
97
98  ! See Emanuel, K. and M. Zivkovic-Rothman, 2000: Development and evaluation
99  ! of a
100  ! convective scheme for use in climate models.
101
102  ! FILES: None
103
104  ! INCLUDE FILES: None
105
106  ! MAKEFILE: /a/ops/met/nogaps/src/sub/fcst/fcst159lib.mak
107
108  ! ..............................END PROLOGUE.............................
109
110
111  USE dimphy
112  IMPLICIT NONE
113
114  ! ym#include "dimensions.h"
115  ! ym#include "dimphy.h"
116
117  INTEGER len
118  INTEGER nd
119  INTEGER ndp1
120  INTEGER noff
121  REAL t(len, nd)
122  REAL q(len, nd)
123  REAL qs(len, nd)
124  REAL u(len, nd)
125  REAL v(len, nd)
126  REAL p(len, nd)
127  REAL ph(len, ndp1)
128  INTEGER iflag(len)
129  REAL ft(len, nd)
130  REAL fq(len, nd)
131  REAL fu(len, nd)
132  REAL fv(len, nd)
133  REAL precip(len)
134  REAL cbmf(len)
135  REAL ma(len, nd)
136  INTEGER minorig
137  REAL delt, cpd, cpv, cl, rv, rd, lv0, g
138  REAL sigs, sigd, elcrit, tlcrit, omtsnow, dtmax, damp
139  REAL alpha, entp, coeffs, coeffr, omtrain, cu
140
141  ! -------------------------------------------------------------------
142  ! --- ARGUMENTS
143  ! -------------------------------------------------------------------
144  ! --- On input:
145
146  ! t:   Array of absolute temperature (K) of dimension ND, with first
147  ! index corresponding to lowest model level. Note that this array
148  ! will be altered by the subroutine if dry convective adjustment
149  ! occurs and if IPBL is not equal to 0.
150
151  ! q:   Array of specific humidity (gm/gm) of dimension ND, with first
152  ! index corresponding to lowest model level. Must be defined
153  ! at same grid levels as T. Note that this array will be altered
154  ! if dry convective adjustment occurs and if IPBL is not equal to 0.
155
156  ! qs:  Array of saturation specific humidity of dimension ND, with first
157  ! index corresponding to lowest model level. Must be defined
158  ! at same grid levels as T. Note that this array will be altered
159  ! if dry convective adjustment occurs and if IPBL is not equal to 0.
160
161  ! u:   Array of zonal wind velocity (m/s) of dimension ND, witth first
162  ! index corresponding with the lowest model level. Defined at
163  ! same levels as T. Note that this array will be altered if
164  ! dry convective adjustment occurs and if IPBL is not equal to 0.
165
166  ! v:   Same as u but for meridional velocity.
167
168  ! tra: Array of passive tracer mixing ratio, of dimensions (ND,NTRA),
169  ! where NTRA is the number of different tracers. If no
170  ! convective tracer transport is needed, define a dummy
171  ! input array of dimension (ND,1). Tracers are defined at
172  ! same vertical levels as T. Note that this array will be altered
173  ! if dry convective adjustment occurs and if IPBL is not equal to 0.
174
175  ! p:   Array of pressure (mb) of dimension ND, with first
176  ! index corresponding to lowest model level. Must be defined
177  ! at same grid levels as T.
178
179  ! ph:  Array of pressure (mb) of dimension ND+1, with first index
180  ! corresponding to lowest level. These pressures are defined at
181  ! levels intermediate between those of P, T, Q and QS. The first
182  ! value of PH should be greater than (i.e. at a lower level than)
183  ! the first value of the array P.
184
185  ! nl:  The maximum number of levels to which convection can penetrate, plus
186  ! 1.
187  ! NL MUST be less than or equal to ND-1.
188
189  ! delt: The model time step (sec) between calls to CONVECT
190
191  ! ----------------------------------------------------------------------------
192  ! ---   On Output:
193
194  ! iflag: An output integer whose value denotes the following:
195  ! VALUE   INTERPRETATION
196  ! -----   --------------
197  ! 0     Moist convection occurs.
198  ! 1     Moist convection occurs, but a CFL condition
199  ! on the subsidence warming is violated. This
200  ! does not cause the scheme to terminate.
201  ! 2     Moist convection, but no precip because ep(inb) lt 0.0001
202  ! 3     No moist convection because new cbmf is 0 and old cbmf is 0.
203  ! 4     No moist convection; atmosphere is not
204  ! unstable
205  ! 6     No moist convection because ihmin le minorig.
206  ! 7     No moist convection because unreasonable
207  ! parcel level temperature or specific humidity.
208  ! 8     No moist convection: lifted condensation
209  ! level is above the 200 mb level.
210  ! 9     No moist convection: cloud base is higher
211  ! then the level NL-1.
212
213  ! ft:   Array of temperature tendency (K/s) of dimension ND, defined at
214  ! same
215  ! grid levels as T, Q, QS and P.
216
217  ! fq:   Array of specific humidity tendencies ((gm/gm)/s) of dimension ND,
218  ! defined at same grid levels as T, Q, QS and P.
219
220  ! fu:   Array of forcing of zonal velocity (m/s^2) of dimension ND,
221  ! defined at same grid levels as T.
222
223  ! fv:   Same as FU, but for forcing of meridional velocity.
224
225  ! ftra: Array of forcing of tracer content, in tracer mixing ratio per
226  ! second, defined at same levels as T. Dimensioned (ND,NTRA).
227
228  ! precip: Scalar convective precipitation rate (mm/day).
229
230  ! wd:   A convective downdraft velocity scale. For use in surface
231  ! flux parameterizations. See convect.ps file for details.
232
233  ! tprime: A convective downdraft temperature perturbation scale (K).
234  ! For use in surface flux parameterizations. See convect.ps
235  ! file for details.
236
237  ! qprime: A convective downdraft specific humidity
238  ! perturbation scale (gm/gm).
239  ! For use in surface flux parameterizations. See convect.ps
240  ! file for details.
241
242  ! cbmf: The cloud base mass flux ((kg/m**2)/s). THIS SCALAR VALUE MUST
243  ! BE STORED BY THE CALLING PROGRAM AND RETURNED TO CONVECT AT
244  ! ITS NEXT CALL. That is, the value of CBMF must be "remembered"
245  ! by the calling program between calls to CONVECT.
246
247  ! det:   Array of detrainment mass flux of dimension ND.
248
249  ! -------------------------------------------------------------------
250
251  ! Local arrays
252
253  INTEGER nl
254  INTEGER nlp
255  INTEGER nlm
256  INTEGER i, k, n
257  REAL delti
258  REAL rowl
259  REAL clmcpv
260  REAL clmcpd
261  REAL cpdmcp
262  REAL cpvmcpd
263  REAL eps
264  REAL epsi
265  REAL epsim1
266  REAL ginv
267  REAL hrd
268  REAL prccon1
269  INTEGER icbmax
270  REAL lv(klon, klev)
271  REAL cpn(klon, klev)
272  REAL cpx(klon, klev)
273  REAL tv(klon, klev)
274  REAL gz(klon, klev)
275  REAL hm(klon, klev)
276  REAL h(klon, klev)
277  REAL work(klon)
278  INTEGER ihmin(klon)
279  INTEGER nk(klon)
280  REAL rh(klon)
281  REAL chi(klon)
282  REAL plcl(klon)
283  INTEGER icb(klon)
284  REAL tnk(klon)
285  REAL qnk(klon)
286  REAL gznk(klon)
287  REAL pnk(klon)
288  REAL qsnk(klon)
289  REAL ticb(klon)
290  REAL gzicb(klon)
291  REAL tp(klon, klev)
292  REAL tvp(klon, klev)
293  REAL clw(klon, klev)
294
295  REAL ah0(klon), cpp(klon)
296  REAL tg, qg, s, alv, tc, ahg, denom, es, rg
297
298  INTEGER ncum
299  INTEGER idcum(klon)
300
301  cpd = 1005.7
302  cpv = 1870.0
303  cl = 4190.0
304  rv = 461.5
305  rd = 287.04
306  lv0 = 2.501E6
307  g = 9.8
308
309  ! *** ELCRIT IS THE AUTOCONVERSION THERSHOLD WATER CONTENT (gm/gm) ***
310  ! ***  TLCRIT IS CRITICAL TEMPERATURE BELOW WHICH THE AUTO-        ***
311  ! ***       CONVERSION THRESHOLD IS ASSUMED TO BE ZERO             ***
312  ! ***     (THE AUTOCONVERSION THRESHOLD VARIES LINEARLY            ***
313  ! ***               BETWEEN 0 C AND TLCRIT)                        ***
314  ! ***   ENTP IS THE COEFFICIENT OF MIXING IN THE ENTRAINMENT       ***
315  ! ***                       FORMULATION                            ***
316  ! ***  SIGD IS THE FRACTIONAL AREA COVERED BY UNSATURATED DNDRAFT  ***
317  ! ***  SIGS IS THE FRACTION OF PRECIPITATION FALLING OUTSIDE       ***
318  ! ***                        OF CLOUD                              ***
319  ! ***        OMTRAIN IS THE ASSUMED FALL SPEED (P/s) OF RAIN       ***
320  ! ***     OMTSNOW IS THE ASSUMED FALL SPEED (P/s) OF SNOW          ***
321  ! ***  COEFFR IS A COEFFICIENT GOVERNING THE RATE OF EVAPORATION   ***
322  ! ***                          OF RAIN                             ***
323  ! ***  COEFFS IS A COEFFICIENT GOVERNING THE RATE OF EVAPORATION   ***
324  ! ***                          OF SNOW                             ***
325  ! ***     CU IS THE COEFFICIENT GOVERNING CONVECTIVE MOMENTUM      ***
326  ! ***                         TRANSPORT                            ***
327  ! ***    DTMAX IS THE MAXIMUM NEGATIVE TEMPERATURE PERTURBATION    ***
328  ! ***        A LIFTED PARCEL IS ALLOWED TO HAVE BELOW ITS LFC      ***
329  ! ***    ALPHA AND DAMP ARE PARAMETERS THAT CONTROL THE RATE OF    ***
330  ! ***                 APPROACH TO QUASI-EQUILIBRIUM                ***
331  ! ***   (THEIR STANDARD VALUES ARE  0.20 AND 0.1, RESPECTIVELY)    ***
332  ! ***                   (DAMP MUST BE LESS THAN 1)                 ***
333
334  sigs = 0.12
335  sigd = 0.05
336  elcrit = 0.0011
337  tlcrit = -55.0
338  omtsnow = 5.5
339  dtmax = 0.9
340  damp = 0.1
341  alpha = 0.2
342  entp = 1.5
343  coeffs = 0.8
344  coeffr = 1.0
345  omtrain = 50.0
346
347  cu = 0.70
348  damp = 0.1
349
350
351  ! Define nl, nlp, nlm, and delti
352
353  nl = nd - noff
354  nlp = nl + 1
355  nlm = nl - 1
356  delti = 1.0/delt
357
358  ! -------------------------------------------------------------------
359  ! --- SET CONSTANTS
360  ! -------------------------------------------------------------------
361
362  rowl = 1000.0
363  clmcpv = cl - cpv
364  clmcpd = cl - cpd
365  cpdmcp = cpd - cpv
366  cpvmcpd = cpv - cpd
367  eps = rd/rv
368  epsi = 1.0/eps
369  epsim1 = epsi - 1.0
370  ginv = 1.0/g
371  hrd = 0.5*rd
372  prccon1 = 86400.0*1000.0/(rowl*g)
373
374  ! dtmax is the maximum negative temperature perturbation.
375
376  ! =====================================================================
377  ! --- INITIALIZE OUTPUT ARRAYS AND PARAMETERS
378  ! =====================================================================
379
380  DO k = 1, nd
381    DO i = 1, len
382      ft(i, k) = 0.0
383      fq(i, k) = 0.0
384      fu(i, k) = 0.0
385      fv(i, k) = 0.0
386      tvp(i, k) = 0.0
387      tp(i, k) = 0.0
388      clw(i, k) = 0.0
389      gz(i, k) = 0.
390    END DO
391  END DO
392  DO i = 1, len
393    precip(i) = 0.0
394    iflag(i) = 0
395  END DO
396
397  ! =====================================================================
398  ! --- CALCULATE ARRAYS OF GEOPOTENTIAL, HEAT CAPACITY & STATIC ENERGY
399  ! =====================================================================
400  DO k = 1, nl + 1
401    DO i = 1, len
402      lv(i, k) = lv0 - clmcpv*(t(i,k)-273.15)
403      cpn(i, k) = cpd*(1.0-q(i,k)) + cpv*q(i, k)
404      cpx(i, k) = cpd*(1.0-q(i,k)) + cl*q(i, k)
405      tv(i, k) = t(i, k)*(1.0+q(i,k)*epsim1)
406    END DO
407  END DO
408
409  ! gz = phi at the full levels (same as p).
410
411  DO i = 1, len
412    gz(i, 1) = 0.0
413  END DO
414  DO k = 2, nlp
415    DO i = 1, len
416      gz(i, k) = gz(i, k-1) + hrd*(tv(i,k-1)+tv(i,k))*(p(i,k-1)-p(i,k))/ph(i, &
417        k)
418    END DO
419  END DO
420
421  ! h  = phi + cpT (dry static energy).
422  ! hm = phi + cp(T-Tbase)+Lq
423
424  DO k = 1, nlp
425    DO i = 1, len
426      h(i, k) = gz(i, k) + cpn(i, k)*t(i, k)
427      hm(i, k) = gz(i, k) + cpx(i, k)*(t(i,k)-t(i,1)) + lv(i, k)*q(i, k)
428    END DO
429  END DO
430
431  ! -------------------------------------------------------------------
432  ! --- Find level of minimum moist static energy
433  ! --- If level of minimum moist static energy coincides with
434  ! --- or is lower than minimum allowable parcel origin level,
435  ! --- set iflag to 6.
436  ! -------------------------------------------------------------------
437  DO i = 1, len
438    work(i) = 1.0E12
439    ihmin(i) = nl
440  END DO
441  DO k = 2, nlp
442    DO i = 1, len
443      IF ((hm(i,k)<work(i)) .AND. (hm(i,k)<hm(i,k-1))) THEN
444        work(i) = hm(i, k)
445        ihmin(i) = k
446      END IF
447    END DO
448  END DO
449  DO i = 1, len
450    ihmin(i) = min(ihmin(i), nlm)
451    IF (ihmin(i)<=minorig) THEN
452      iflag(i) = 6
453    END IF
454  END DO
455
456  ! -------------------------------------------------------------------
457  ! --- Find that model level below the level of minimum moist static
458  ! --- energy that has the maximum value of moist static energy
459  ! -------------------------------------------------------------------
460
461  DO i = 1, len
462    work(i) = hm(i, minorig)
463    nk(i) = minorig
464  END DO
465  DO k = minorig + 1, nl
466    DO i = 1, len
467      IF ((hm(i,k)>work(i)) .AND. (k<=ihmin(i))) THEN
468        work(i) = hm(i, k)
469        nk(i) = k
470      END IF
471    END DO
472  END DO
473  ! -------------------------------------------------------------------
474  ! --- Check whether parcel level temperature and specific humidity
475  ! --- are reasonable
476  ! -------------------------------------------------------------------
477  DO i = 1, len
478    IF (((t(i,nk(i))<250.0) .OR. (q(i,nk(i))<=0.0) .OR. (p(i,ihmin(i))< &
479      400.0)) .AND. (iflag(i)==0)) iflag(i) = 7
480  END DO
481  ! -------------------------------------------------------------------
482  ! --- Calculate lifted condensation level of air at parcel origin level
483  ! --- (Within 0.2% of formula of Bolton, MON. WEA. REV.,1980)
484  ! -------------------------------------------------------------------
485  DO i = 1, len
486    tnk(i) = t(i, nk(i))
487    qnk(i) = q(i, nk(i))
488    gznk(i) = gz(i, nk(i))
489    pnk(i) = p(i, nk(i))
490    qsnk(i) = qs(i, nk(i))
491
492    rh(i) = qnk(i)/qsnk(i)
493    rh(i) = min(1.0, rh(i))
494    chi(i) = tnk(i)/(1669.0-122.0*rh(i)-tnk(i))
495    plcl(i) = pnk(i)*(rh(i)**chi(i))
496    IF (((plcl(i)<200.0) .OR. (plcl(i)>=2000.0)) .AND. (iflag(i)==0)) iflag(i &
497      ) = 8
498  END DO
499  ! -------------------------------------------------------------------
500  ! --- Calculate first level above lcl (=icb)
501  ! -------------------------------------------------------------------
502  DO i = 1, len
503    icb(i) = nlm
504  END DO
505
506  DO k = minorig, nl
507    DO i = 1, len
508      IF ((k>=(nk(i)+1)) .AND. (p(i,k)<plcl(i))) icb(i) = min(icb(i), k)
509    END DO
510  END DO
511
512  DO i = 1, len
513    IF ((icb(i)>=nlm) .AND. (iflag(i)==0)) iflag(i) = 9
514  END DO
515
516  ! Compute icbmax.
517
518  icbmax = 2
519  DO i = 1, len
520    icbmax = max(icbmax, icb(i))
521  END DO
522
523  ! -------------------------------------------------------------------
524  ! --- Calculates the lifted parcel virtual temperature at nk,
525  ! --- the actual temperature, and the adiabatic
526  ! --- liquid water content. The procedure is to solve the equation.
527  ! cp*tp+L*qp+phi=cp*tnk+L*qnk+gznk.
528  ! -------------------------------------------------------------------
529
530  DO i = 1, len
531    tnk(i) = t(i, nk(i))
532    qnk(i) = q(i, nk(i))
533    gznk(i) = gz(i, nk(i))
534    ticb(i) = t(i, icb(i))
535    gzicb(i) = gz(i, icb(i))
536  END DO
537
538  ! ***  Calculate certain parcel quantities, including static energy   ***
539
540  DO i = 1, len
541    ah0(i) = (cpd*(1.-qnk(i))+cl*qnk(i))*tnk(i) + qnk(i)*(lv0-clmcpv*(tnk(i)- &
542      273.15)) + gznk(i)
543    cpp(i) = cpd*(1.-qnk(i)) + qnk(i)*cpv
544  END DO
545
546  ! ***   Calculate lifted parcel quantities below cloud base   ***
547
548  DO k = minorig, icbmax - 1
549    DO i = 1, len
550      tp(i, k) = tnk(i) - (gz(i,k)-gznk(i))/cpp(i)
551      tvp(i, k) = tp(i, k)*(1.+qnk(i)*epsi)
552    END DO
553  END DO
554
555  ! ***  Find lifted parcel quantities above cloud base    ***
556
557  DO i = 1, len
558    tg = ticb(i)
559    qg = qs(i, icb(i))
560    alv = lv0 - clmcpv*(ticb(i)-273.15)
561
562    ! First iteration.
563
564    s = cpd + alv*alv*qg/(rv*ticb(i)*ticb(i))
565    s = 1./s
566    ahg = cpd*tg + (cl-cpd)*qnk(i)*ticb(i) + alv*qg + gzicb(i)
567    tg = tg + s*(ah0(i)-ahg)
568    tg = max(tg, 35.0)
569    tc = tg - 273.15
570    denom = 243.5 + tc
571    IF (tc>=0.0) THEN
572      es = 6.112*exp(17.67*tc/denom)
573    ELSE
574      es = exp(23.33086-6111.72784/tg+0.15215*log(tg))
575    END IF
576    qg = eps*es/(p(i,icb(i))-es*(1.-eps))
577
578    ! Second iteration.
579
580    s = cpd + alv*alv*qg/(rv*ticb(i)*ticb(i))
581    s = 1./s
582    ahg = cpd*tg + (cl-cpd)*qnk(i)*ticb(i) + alv*qg + gzicb(i)
583    tg = tg + s*(ah0(i)-ahg)
584    tg = max(tg, 35.0)
585    tc = tg - 273.15
586    denom = 243.5 + tc
587    IF (tc>=0.0) THEN
588      es = 6.112*exp(17.67*tc/denom)
589    ELSE
590      es = exp(23.33086-6111.72784/tg+0.15215*log(tg))
591    END IF
592    qg = eps*es/(p(i,icb(i))-es*(1.-eps))
593
594    alv = lv0 - clmcpv*(ticb(i)-273.15)
595    tp(i, icb(i)) = (ah0(i)-(cl-cpd)*qnk(i)*ticb(i)-gz(i,icb(i))-alv*qg)/cpd
596    clw(i, icb(i)) = qnk(i) - qg
597    clw(i, icb(i)) = max(0.0, clw(i,icb(i)))
598    rg = qg/(1.-qnk(i))
599    tvp(i, icb(i)) = tp(i, icb(i))*(1.+rg*epsi)
600  END DO
601
602  DO k = minorig, icbmax
603    DO i = 1, len
604      tvp(i, k) = tvp(i, k) - tp(i, k)*qnk(i)
605    END DO
606  END DO
607
608  ! -------------------------------------------------------------------
609  ! --- Test for instability.
610  ! --- If there was no convection at last time step and parcel
611  ! --- is stable at icb, then set iflag to 4.
612  ! -------------------------------------------------------------------
613
614  DO i = 1, len
615    IF ((cbmf(i)==0.0) .AND. (iflag(i)==0) .AND. (tvp(i, &
616      icb(i))<=(tv(i,icb(i))-dtmax))) iflag(i) = 4
617  END DO
618
619  ! =====================================================================
620  ! --- IF THIS POINT IS REACHED, MOIST CONVECTIVE ADJUSTMENT IS NECESSARY
621  ! =====================================================================
622
623  ncum = 0
624  DO i = 1, len
625    IF (iflag(i)==0) THEN
626      ncum = ncum + 1
627      idcum(ncum) = i
628    END IF
629  END DO
630
631  ! Call convect2, which compresses the points and computes the heating,
632  ! moistening, velocity mixing, and precipiation.
633
634  ! print*,'cpd avant convect2 ',cpd
635  IF (ncum>0) THEN
636    CALL convect2(ncum, idcum, len, nd, ndp1, nl, minorig, nk, icb, t, q, qs, &
637      u, v, gz, tv, tp, tvp, clw, h, lv, cpn, p, ph, ft, fq, fu, fv, tnk, &
638      qnk, gznk, plcl, precip, cbmf, iflag, delt, cpd, cpv, cl, rv, rd, lv0, &
639      g, sigs, sigd, elcrit, tlcrit, omtsnow, dtmax, damp, alpha, entp, &
640      coeffs, coeffr, omtrain, cu, ma)
641  END IF
642
643  RETURN
644END SUBROUTINE convect1
Note: See TracBrowser for help on using the repository browser.