source: LMDZ4/trunk/libf/phy_IPCC_AR4/conccm.F @ 932

Last change on this file since 932 was 868, checked in by Laurent Fairhead, 17 years ago

Preparation du remplacement de la physique utilisee pour l'exercice IPCC_AR4
par la version de la physique avec thermique. On garde le repertoire phylmd
pour un petit moment pour que les utilisateurs ne soient pas trop perdus ...
phy_IPCC_AR4 = phylmd
LF
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 30.8 KB
RevLine 
[868]1!
2! $Header$
3!
4      SUBROUTINE conccm (dtime,paprs,pplay,t,q,conv_q,
5     s                   d_t, d_q, rain, snow, kbascm, ktopcm)
6c
7      USE dimphy
8      IMPLICIT none
9c======================================================================
10c Auteur(s): Z.X. Li (LMD/CNRS) date: le 14 mars 1996
11c Objet: Schema simple (avec flux de masse) pour la convection
12c        (schema standard du modele NCAR CCM2)
13c======================================================================
14cym#include "dimensions.h"
15cym#include "dimphy.h"
16#include "YOMCST.h"
17#include "YOETHF.h"
18c
19c Entree:
20      REAL dtime              ! pas d'integration
21      REAL paprs(klon,klev+1) ! pression inter-couche (Pa)
22      REAL pplay(klon,klev)   ! pression au milieu de couche (Pa)
23      REAL t(klon,klev)       ! temperature (K)
24      REAL q(klon,klev)       ! humidite specifique (g/g)
25      REAL conv_q(klon,klev)  ! taux de convergence humidite (g/g/s)
26c Sortie:
27      REAL d_t(klon,klev)     ! incrementation temperature
28      REAL d_q(klon,klev)     ! incrementation vapeur
29      REAL rain(klon)         ! pluie (mm/s)
30      REAL snow(klon)         ! neige (mm/s)
31      INTEGER kbascm(klon)    ! niveau du bas de convection
32      INTEGER ktopcm(klon)    ! niveau du haut de convection
33c
34      REAL pt(klon,klev)
35      REAL pq(klon,klev)
36      REAL pres(klon,klev)
37      REAL dp(klon,klev)
38      REAL zgeom(klon,klev)
39      REAL cmfprs(klon)
40      REAL cmfprt(klon)
41      INTEGER ntop(klon)
42      INTEGER nbas(klon)
43      INTEGER i, k
44      REAL zlvdcp, zlsdcp, zdelta, zz, za, zb
45c
46      LOGICAL usekuo ! utiliser convection profonde (schema Kuo)
47      PARAMETER (usekuo=.TRUE.)
48c
49      REAL d_t_bis(klon,klev)
50      REAL d_q_bis(klon,klev)
51      REAL rain_bis(klon)
52      REAL snow_bis(klon)
53      INTEGER ibas_bis(klon)
54      INTEGER itop_bis(klon)
55      REAL d_ql_bis(klon,klev)
56      REAL rneb_bis(klon,klev)
57c
58c initialiser les variables de sortie (pour securite)
59      DO i = 1, klon
60         rain(i) = 0.0
61         snow(i) = 0.0
62         kbascm(i) = 0
63         ktopcm(i) = 0
64      ENDDO
65      DO k = 1, klev
66      DO i = 1, klon
67         d_t(i,k) = 0.0
68         d_q(i,k) = 0.0
69      ENDDO
70      ENDDO
71c
72c preparer les variables d'entree (attention: l'ordre des niveaux
73c verticaux augmente du haut vers le bas)
74      DO k = 1, klev
75      DO i = 1, klon
76         pt(i,k) = t(i,klev-k+1)
77         pq(i,k) = q(i,klev-k+1)
78         pres(i,k) = pplay(i,klev-k+1)
79         dp(i,k) = paprs(i,klev+1-k)-paprs(i,klev+1-k+1)
80      ENDDO
81      ENDDO
82      DO i = 1, klon
83         zgeom(i,klev) = RD * t(i,1) / (0.5*(paprs(i,1)+pplay(i,1)))
84     .                      * (paprs(i,1)-pplay(i,1))
85      ENDDO
86      DO k = 2, klev
87      DO i = 1, klon
88         zgeom(i,klev+1-k) = zgeom(i,klev+1-k+1)
89     .              + RD * 0.5*(t(i,k-1)+t(i,k)) / paprs(i,k)
90     .                   * (pplay(i,k-1)-pplay(i,k))
91      ENDDO
92      ENDDO
93c
94      CALL cmfmca(dtime, pres, dp, zgeom, pt, pq,
95     $                  cmfprt, cmfprs, ntop, nbas)
96C
97      DO k = 1, klev
98      DO i = 1, klon
99         d_q(i,klev+1-k) = pq(i,k) - q(i,klev+1-k)
100         d_t(i,klev+1-k) = pt(i,k) - t(i,klev+1-k)
101      ENDDO
102      ENDDO
103c
104      DO i = 1, klon
105         rain(i) = cmfprt(i) * rhoh2o
106         snow(i) = cmfprs(i) * rhoh2o
107         kbascm(i) = klev+1 - nbas(i)
108         ktopcm(i) = klev+1 - ntop(i)
109      ENDDO
110c
111      IF (usekuo) THEN
112      CALL conkuo(dtime, paprs, pplay, t, q, conv_q,
113     s            d_t_bis, d_q_bis, d_ql_bis, rneb_bis,
114     s            rain_bis, snow_bis, ibas_bis, itop_bis)
115      DO k = 1, klev
116      DO i = 1, klon
117         d_t(i,k) = d_t(i,k) + d_t_bis(i,k)
118         d_q(i,k) = d_q(i,k) + d_q_bis(i,k)
119      ENDDO
120      ENDDO
121      DO i = 1, klon
122         rain(i) = rain(i) + rain_bis(i)
123         snow(i) = snow(i) + snow_bis(i)
124         kbascm(i) = MIN(kbascm(i),ibas_bis(i))
125         ktopcm(i) = MAX(ktopcm(i),itop_bis(i))
126      ENDDO
127      DO k = 1, klev ! eau liquide convective est
128      DO i = 1, klon ! dispersee dans l'air
129         zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*q(i,k))
130         zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*q(i,k))
131         zdelta = MAX(0.,SIGN(1.,RTT-t(i,k)))
132         zz = d_ql_bis(i,k) ! re-evap. de l'eau liquide
133         zb = MAX(0.0,zz)
134         za = - MAX(0.0,zz) * (zlvdcp*(1.-zdelta)+zlsdcp*zdelta)
135         d_t(i,k) = d_t(i,k) + za
136         d_q(i,k) = d_q(i,k) + zb
137      ENDDO
138      ENDDO
139      ENDIF
140c
141      RETURN
142      END
143      SUBROUTINE cmfmca(deltat, p, dp, gz,
144     $                  tb, shb,
145     $                  cmfprt, cmfprs, cnt, cnb)
146      USE dimphy
147      IMPLICIT none
148C-----------------------------------------------------------------------
149C Moist convective mass flux procedure:
150C If stratification is unstable to nonentraining parcel ascent,
151C complete an adjustment making use of a simple cloud model
152C
153C Code generalized to allow specification of parcel ("updraft")
154C properties, as well as convective transport of an arbitrary
155C number of passive constituents (see cmrb array).
156C----------------------------Code History-------------------------------
157C Original version:  J. J. Hack, March 22, 1990
158C Standardized:      J. Rosinski, June 1992
159C Reviewed:          J. Hack, G. Taylor, August 1992
160c Adaptation au LMD: Z.X. Li, mars 1996 (reference: Hack 1994,
161c                    J. Geophys. Res. vol 99, D3, 5551-5568). J'ai
162c                    introduit les constantes et les fonctions thermo-
163c                    dynamiques du Centre Europeen. J'ai elimine le
164c                    re-indicage du code en esperant que cela pourra
165c                    simplifier la lecture et la comprehension.
166C-----------------------------------------------------------------------
167cym#include "dimensions.h"
168cym#include "dimphy.h"
169      INTEGER pcnst ! nombre de traceurs passifs
170      PARAMETER (pcnst=1)
171C------------------------------Arguments--------------------------------
172C Input arguments
173C
174      REAL deltat                 ! time step (seconds)
175      REAL p(klon,klev)           ! pressure
176      REAL dp(klon,klev)          ! delta-p
177      REAL gz(klon,klev)          ! geopotential (a partir du sol)
178c
179      REAL thtap(klon)            ! PBL perturbation theta
180      REAL shp(klon)              ! PBL perturbation specific humidity
181      REAL pblh(klon)             ! PBL height (provided by PBL routine)
182      REAL cmrp(klon,pcnst)       ! constituent perturbations in PBL
183c
184c Updated arguments:
185c
186      REAL tb(klon,klev)         ! temperature (t bar)
187      REAL shb(klon,klev)        ! specific humidity (sh bar)
188      REAL cmrb(klon,klev,pcnst) ! constituent mixing ratios (cmr bar)
189C
190C Output arguments
191C
192      REAL cmfdt(klon,klev)    ! dT/dt due to moist convection
193      REAL cmfdq(klon,klev)    ! dq/dt due to moist convection
194      REAL cmfmc(klon,klev )   ! moist convection cloud mass flux
195      REAL cmfdqr(klon,klev)   ! dq/dt due to convective rainout
196      REAL cmfsl(klon,klev )   ! convective lw static energy flux
197      REAL cmflq(klon,klev )   ! convective total water flux
198      REAL cmfprt(klon)        ! convective precipitation rate
199      REAL cmfprs(klon)        ! convective snowfall rate
200      REAL qc(klon,klev)       ! dq/dt due to rainout terms
201      INTEGER cnt(klon)        ! top level of convective activity   
202      INTEGER cnb(klon)        ! bottom level of convective activity
203C------------------------------Parameters-------------------------------
204      REAL c0         ! rain water autoconversion coefficient
205      PARAMETER (c0=1.0e-4)
206      REAL dzmin       ! minimum convective depth for precipitation
207      PARAMETER (dzmin=0.0)
208      REAL betamn      ! minimum overshoot parameter
209      PARAMETER (betamn=0.10)
210      REAL cmftau      ! characteristic adjustment time scale
211      PARAMETER (cmftau=3600.)
212      INTEGER limcnv   ! top interface level limit for convection
213      PARAMETER (limcnv=1)
214      REAL tpmax       ! maximum acceptable t perturbation (degrees C)
215      PARAMETER (tpmax=1.50)
216      REAL shpmax      ! maximum acceptable q perturbation (g/g)
217      PARAMETER (shpmax=1.50e-3)
218      REAL tiny        ! arbitrary small num used in transport estimates
219      PARAMETER (tiny=1.0e-36)
220      REAL eps         ! convergence criteria (machine dependent)
221      PARAMETER (eps=1.0e-13)
222      REAL tmelt       ! freezing point of water(req'd for rain vs snow)
223      PARAMETER (tmelt=273.15)
224      REAL ssfac ! supersaturation bound (detrained air)
225      PARAMETER (ssfac=1.001)
226C
227C---------------------------Local workspace-----------------------------
228      REAL gam(klon,klev)     ! L/cp (d(qsat)/dT)
229      REAL sb(klon,klev)      ! dry static energy (s bar)
230      REAL hb(klon,klev)      ! moist static energy (h bar)
231      REAL shbs(klon,klev)    ! sat. specific humidity (sh bar star)
232      REAL hbs(klon,klev)     ! sat. moist static energy (h bar star)
233      REAL shbh(klon,klev+1)  ! specific humidity on interfaces
234      REAL sbh(klon,klev+1)   ! s bar on interfaces
235      REAL hbh(klon,klev+1)   ! h bar on interfaces
236      REAL cmrh(klon,klev+1)  ! interface constituent mixing ratio
237      REAL prec(klon)         ! instantaneous total precipitation
238      REAL dzcld(klon)        ! depth of convective layer (m)
239      REAL beta(klon)         ! overshoot parameter (fraction)
240      REAL betamx             ! local maximum on overshoot
241      REAL eta(klon)          ! convective mass flux (kg/m^2 s)
242      REAL etagdt             ! eta*grav*deltat
243      REAL cldwtr(klon)       ! cloud water (mass)
244      REAL rnwtr(klon)        ! rain water  (mass)
245      REAL sc  (klon)         ! dry static energy   ("in-cloud")
246      REAL shc (klon)         ! specific humidity   ("in-cloud")
247      REAL hc  (klon)         ! moist static energy ("in-cloud")
248      REAL cmrc(klon)         ! constituent mix rat ("in-cloud")
249      REAL dq1(klon)          ! shb  convective change (lower lvl)
250      REAL dq2(klon)          ! shb  convective change (mid level)
251      REAL dq3(klon)          ! shb  convective change (upper lvl)
252      REAL ds1(klon)          ! sb   convective change (lower lvl)
253      REAL ds2(klon)          ! sb   convective change (mid level)
254      REAL ds3(klon)          ! sb   convective change (upper lvl)
255      REAL dcmr1(klon)        ! cmrb convective change (lower lvl)
256      REAL dcmr2(klon)        ! cmrb convective change (mid level)
257      REAL dcmr3(klon)        ! cmrb convective change (upper lvl)
258      REAL flotab(klon)       ! hc - hbs (mesure d'instabilite)
259      LOGICAL ldcum(klon)     ! .true. si la convection existe
260      LOGICAL etagt0          ! true if eta > 0.0
261      REAL dt                 ! current 2 delta-t (model time step)
262      REAL cats     ! modified characteristic adj. time
263      REAL rdt      ! 1./dt
264      REAL qprime   ! modified specific humidity pert.
265      REAL tprime   ! modified thermal perturbation
266      REAL pblhgt   ! bounded pbl height (max[pblh,1m])
267      REAL fac1     ! intermediate scratch variable
268      REAL shprme   ! intermediate specific humidity pert.
269      REAL qsattp   ! saturation mixing ratio for
270C                   !  thermally perturbed PBL parcels
271      REAL dz       ! local layer depth
272      REAL b1       ! bouyancy measure in detrainment lvl
273      REAL b2       ! bouyancy measure in condensation lvl
274      REAL g     ! bounded vertical gradient of hb
275      REAL tmass ! total mass available for convective exchange
276      REAL denom ! intermediate scratch variable
277      REAL qtest1! used in negative q test (middle lvl)
278      REAL qtest2! used in negative q test (lower lvl)
279      REAL fslkp ! flux lw static energy (bot interface)
280      REAL fslkm ! flux lw static energy (top interface)
281      REAL fqlkp ! flux total water (bottom interface)
282      REAL fqlkm ! flux total water (top interface)
283      REAL botflx! bottom constituent mixing ratio flux
284      REAL topflx! top constituent mixing ratio flux
285      REAL efac1 ! ratio cmrb to convectively induced change (bot lvl)
286      REAL efac2 ! ratio cmrb to convectively induced change (mid lvl)
287      REAL efac3 ! ratio cmrb to convectively induced change (top lvl)
288c
289      INTEGER i,k  ! indices horizontal et vertical
290      INTEGER km1  ! k-1 (index offset)
291      INTEGER kp1  ! k+1 (index offset)
292      INTEGER m    ! constituent index
293      INTEGER ktp  ! temporary index used to track top
294      INTEGER is   ! nombre de points a ajuster
295C
296      REAL tmp1, tmp2, tmp3, tmp4
297      REAL zx_t, zx_p, zx_q, zx_qs, zx_gam
298      REAL zcor, zdelta, zcvm5
299C
300      REAL qhalf, sh1, sh2, shbs1, shbs2
301#include "YOMCST.h"
302#include "YOETHF.h"
303#include "FCTTRE.h"
304      qhalf(sh1,sh2,shbs1,shbs2) = MIN(MAX(sh1,sh2),
305     $                            (shbs2*sh1 + shbs1*sh2)/(shbs1+shbs2))
306C
307C-----------------------------------------------------------------------
308c pas de traceur pour l'instant
309      DO m = 1, pcnst
310      DO k = 1, klev
311      DO i = 1, klon
312         cmrb(i,k,m) = 0.0
313      ENDDO
314      ENDDO
315      ENDDO
316c
317c Les perturbations de la couche limite sont zero pour l'instant
318c
319      DO m = 1, pcnst
320      DO i = 1, klon
321         cmrp(i,m) = 0.0
322      ENDDO
323      ENDDO
324      DO i = 1, klon
325         thtap(i) = 0.0
326         shp(i) = 0.0
327         pblh(i) = 1.0
328      ENDDO
329C
330C Ensure that characteristic adjustment time scale (cmftau) assumed
331C in estimate of eta isn't smaller than model time scale (deltat)
332C
333      dt   = deltat
334      cats = MAX(dt,cmftau)
335      rdt  = 1.0/dt
336C
337C Compute sb,hb,shbs,hbs
338C
339      DO k = 1, klev
340      DO i = 1, klon
341         zx_t = tb(i,k)
342         zx_p = p(i,k)
343         zx_q = shb(i,k)
344           zdelta=MAX(0.,SIGN(1.,RTT-zx_t))
345           zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
346           zcvm5 = zcvm5 / RCPD / (1.0+RVTMP2*zx_q)
347           zx_qs= r2es * FOEEW(zx_t,zdelta)/zx_p
348           zx_qs=MIN(0.5,zx_qs)
349           zcor=1./(1.-retv*zx_qs)
350           zx_qs=zx_qs*zcor
351           zx_gam = FOEDE(zx_t,zdelta,zcvm5,zx_qs,zcor)
352         shbs(i,k) = zx_qs
353         gam(i,k) = zx_gam
354      ENDDO
355      ENDDO
356C
357      DO k=limcnv,klev
358         DO i=1,klon
359            sb (i,k) = RCPD*tb(i,k) + gz(i,k)
360            hb (i,k) = sb(i,k) + RLVTT*shb(i,k)
361            hbs(i,k) = sb(i,k) + RLVTT*shbs(i,k)
362         ENDDO
363      ENDDO
364C
365C Compute sbh, shbh
366C
367      DO k=limcnv+1,klev
368         km1 = k - 1
369         DO i=1,klon
370            sbh (i,k) =0.5*(sb(i,km1) + sb(i,k))
371            shbh(i,k) =qhalf(shb(i,km1),shb(i,k),shbs(i,km1),shbs(i,k))
372            hbh (i,k) =sbh(i,k) + RLVTT*shbh(i,k)
373         ENDDO
374      ENDDO
375C
376C Specify properties at top of model (not used, but filling anyway)
377C
378      DO i=1,klon
379         sbh (i,limcnv) = sb(i,limcnv)
380         shbh(i,limcnv) = shb(i,limcnv)
381         hbh (i,limcnv) = hb(i,limcnv)
382      ENDDO
383C
384C Zero vertically independent control, tendency & diagnostic arrays
385C
386      DO i=1,klon
387         prec(i)  = 0.0
388         dzcld(i) = 0.0
389         cnb(i)   = 0
390         cnt(i)   = klev+1
391      ENDDO
392
393      DO k = 1, klev
394        DO i = 1,klon
395         cmfdt(i,k)  = 0.
396         cmfdq(i,k)  = 0.
397         cmfdqr(i,k) = 0.
398         cmfmc(i,k)  = 0.
399         cmfsl(i,k)  = 0.
400         cmflq(i,k)  = 0.
401        ENDDO
402      ENDDO
403C
404C Begin moist convective mass flux adjustment procedure.
405C Formalism ensures that negative cloud liquid water can never occur
406C
407      DO 70 k=klev-1,limcnv+1,-1
408         km1 = k - 1
409         kp1 = k + 1
410         DO 10 i=1,klon
411            eta   (i) = 0.0
412            beta  (i) = 0.0
413            ds1   (i) = 0.0
414            ds2   (i) = 0.0
415            ds3   (i) = 0.0
416            dq1   (i) = 0.0
417            dq2   (i) = 0.0
418            dq3   (i) = 0.0
419C
420C Specification of "cloud base" conditions
421C
422            qprime    = 0.0
423            tprime    = 0.0
424C
425C Assign tprime within the PBL to be proportional to the quantity
426C thtap (which will be bounded by tpmax), passed to this routine by
427C the PBL routine.  Don't allow perturbation to produce a dry
428C adiabatically unstable parcel.  Assign qprime within the PBL to be
429C an appropriately modified value of the quantity shp (which will be
430C bounded by shpmax) passed to this routine by the PBL routine.  The
431C quantity qprime should be less than the local saturation value
432C (qsattp=qsat[t+tprime,p]).  In both cases, thtap and shp are
433C linearly reduced toward zero as the PBL top is approached.
434C
435            pblhgt = MAX(pblh(i),1.0)
436            IF (gz(i,kp1)/RG.LE.pblhgt .AND. dzcld(i).EQ.0.0) THEN
437               fac1   = MAX(0.0,1.0-gz(i,kp1)/RG/pblhgt)
438               tprime = MIN(thtap(i),tpmax)*fac1
439               qsattp = shbs(i,kp1) + RCPD/RLVTT*gam(i,kp1)*tprime
440               shprme = MIN(MIN(shp(i),shpmax)*fac1,
441     $                        MAX(qsattp-shb(i,kp1),0.0))
442               qprime = MAX(qprime,shprme)
443            ELSE
444               tprime = 0.0
445               qprime = 0.0
446            ENDIF
447C
448C Specify "updraft" (in-cloud) thermodynamic properties
449C
450            sc (i)    = sb (i,kp1) + RCPD*tprime
451            shc(i)    = shb(i,kp1) + qprime
452            hc (i)    = sc (i    ) + RLVTT*shc(i)
453            flotab(i) = hc(i) - hbs(i,k)
454            dz        = dp(i,k)*RD*tb(i,k)/RG/p(i,k)
455            IF (flotab(i).gt.0.0) THEN
456               dzcld(i) = dzcld(i) + dz
457            ELSE
458               dzcld(i) = 0.0
459            ENDIF
460   10    CONTINUE
461C
462C Check on moist convective instability
463C
464         is = 0
465         DO i = 1, klon
466            IF (flotab(i).GT.0.0) THEN
467               ldcum(i) = .TRUE.
468               is = is + 1
469            ELSE
470               ldcum(i) = .FALSE.
471            ENDIF
472         ENDDO
473C
474         IF (is.EQ.0) THEN
475            DO i=1,klon
476               dzcld(i) = 0.0
477            ENDDO
478            GOTO 70
479         ENDIF
480C
481C Current level just below top level => no overshoot
482C
483         IF (k.le.limcnv+1) THEN
484            DO i=1,klon
485            IF (ldcum(i)) THEN
486               cldwtr(i) = sb(i,k)-sc(i)+flotab(i)/(1.0+gam(i,k))
487               cldwtr(i) = MAX(0.0,cldwtr(i))
488               beta(i)   = 0.0
489            ENDIF
490            ENDDO
491            GOTO 20
492         ENDIF
493C
494C First guess at overshoot parameter using crude buoyancy closure
495C 10% overshoot assumed as a minimum and 1-c0*dz maximum to start
496C If pre-existing supersaturation in detrainment layer, beta=0
497C cldwtr is temporarily equal to RLVTT*l (l=> liquid water)
498C
499         DO i=1,klon
500         IF (ldcum(i)) THEN
501            cldwtr(i) = sb(i,k)-sc(i)+flotab(i)/(1.0+gam(i,k))
502            cldwtr(i) = MAX(0.0,cldwtr(i))
503            betamx = 1.0 - c0*MAX(0.0,(dzcld(i)-dzmin))
504            b1        = (hc(i) - hbs(i,km1))*dp(i,km1)
505            b2        = (hc(i) - hbs(i,k  ))*dp(i,k  )
506            beta(i)   = MAX(betamn,MIN(betamx, 1.0+b1/b2))
507            IF (hbs(i,km1).le.hb(i,km1)) beta(i) = 0.0
508         ENDIF
509         ENDDO
510C
511C Bound maximum beta to ensure physically realistic solutions
512C
513C First check constrains beta so that eta remains positive
514C (assuming that eta is already positive for beta equal zero)
515c La premiere contrainte de beta est que le flux eta doit etre positif.
516C
517         DO i=1,klon
518         IF (ldcum(i)) THEN
519            tmp1 = (1.0+gam(i,k))*(sc(i)-sbh(i,kp1) + cldwtr(i))
520     $            - (hbh(i,kp1)-hc(i))*dp(i,k)/dp(i,kp1)
521            tmp2 = (1.0+gam(i,k))*(sc(i)-sbh(i,k))
522            IF ((beta(i)*tmp2-tmp1).GT.0.0) THEN
523               betamx = 0.99*(tmp1/tmp2)
524               beta(i) = MAX(0.0,MIN(betamx,beta(i)))
525            ENDIF
526C
527C Second check involves supersaturation of "detrainment layer"
528C small amount of supersaturation acceptable (by ssfac factor)
529c La 2e contrainte est que la convection ne doit pas sursaturer
530c la "detrainment layer", Neanmoins, une petite sursaturation
531c est acceptee (facteur ssfac).
532C
533            IF (hb(i,km1).lt.hbs(i,km1)) THEN
534               tmp1 = (1.0+gam(i,k))*(sc(i)-sbh(i,kp1) + cldwtr(i))
535     $               - (hbh(i,kp1)-hc(i))*dp(i,k)/dp(i,kp1)
536               tmp1 = tmp1/dp(i,k)
537               tmp2 = gam(i,km1)*(sbh(i,k)-sc(i) + cldwtr(i)) -
538     $                 hbh(i,k) + hc(i) - sc(i) + sbh(i,k)
539               tmp3 = (1.0+gam(i,k))*(sc(i)-sbh(i,k))/dp(i,k)
540               tmp4 = (dt/cats)*(hc(i)-hbs(i,k))*tmp2
541     $               / (dp(i,km1)*(hbs(i,km1)-hb(i,km1))) + tmp3
542               IF ((beta(i)*tmp4-tmp1).GT.0.0) THEN
543                  betamx = ssfac*(tmp1/tmp4)
544                  beta(i)   = MAX(0.0,MIN(betamx,beta(i)))
545               ENDIF
546            ELSE
547               beta(i) = 0.0
548            ENDIF
549C
550C Third check to avoid introducing 2 delta x thermodynamic
551C noise in the vertical ... constrain adjusted h (or theta e)
552C so that the adjustment doesn't contribute to "kinks" in h
553C
554            g = MIN(0.0,hb(i,k)-hb(i,km1))
555            tmp3 = (hb(i,k)-hb(i,km1)-g)*(cats/dt) / (hc(i)-hbs(i,k))
556            tmp1 = (1.0+gam(i,k))*(sc(i)-sbh(i,kp1) + cldwtr(i))
557     $            - (hbh(i,kp1)-hc(i))*dp(i,k)/dp(i,kp1)
558            tmp1 = tmp1/dp(i,k)
559            tmp1 = tmp3*tmp1 + (hc(i) - hbh(i,kp1))/dp(i,k)
560            tmp2 = tmp3*(1.0+gam(i,k))*(sc(i)-sbh(i,k))/dp(i,k)
561     $            + (hc(i)-hbh(i,k)-cldwtr(i))
562     $             *(1.0/dp(i,k)+1.0/dp(i,kp1))
563            IF ((beta(i)*tmp2-tmp1).GT.0.0) THEN
564               betamx = 0.0
565               IF (tmp2.NE.0.0) betamx = tmp1/tmp2
566               beta(i) = MAX(0.0,MIN(betamx,beta(i)))
567            ENDIF
568         ENDIF
569         ENDDO
570C
571C Calculate mass flux required for stabilization.
572C
573C Ensure that the convective mass flux, eta, is positive by
574C setting negative values of eta to zero..
575C Ensure that estimated mass flux cannot move more than the
576C minimum of total mass contained in either layer k or layer k+1.
577C Also test for other pathological cases that result in non-
578C physical states and adjust eta accordingly.
579C
580   20    CONTINUE
581         DO i=1,klon
582         IF (ldcum(i)) THEN
583            beta(i) = MAX(0.0,beta(i))
584            tmp1 = hc(i) - hbs(i,k)
585            tmp2 = ((1.0+gam(i,k))*(sc(i)-sbh(i,kp1)+cldwtr(i)) -
586     $               beta(i)*(1.0+gam(i,k))*(sc(i)-sbh(i,k)))/dp(i,k) -
587     $              (hbh(i,kp1)-hc(i))/dp(i,kp1)
588            eta(i) = tmp1/(tmp2*RG*cats)
589            tmass = MIN(dp(i,k),dp(i,kp1))/RG
590            IF (eta(i).GT.tmass*rdt .OR. eta(i).LE.0.0) eta(i) = 0.0
591C
592C Check on negative q in top layer (bound beta)
593C
594            IF(shc(i)-shbh(i,k).LT.0.0 .AND. beta(i)*eta(i).NE.0.0)THEN
595               denom = eta(i)*RG*dt*(shc(i) - shbh(i,k))/dp(i,km1)
596               beta(i) = MAX(0.0,MIN(-0.999*shb(i,km1)/denom,beta(i)))
597            ENDIF
598C
599C Check on negative q in middle layer (zero eta)
600C
601            qtest1 = shb(i,k) + eta(i)*RG*dt*((shc(i) - shbh(i,kp1)) -
602     $               (1.0 - beta(i))*cldwtr(i)/RLVTT -
603     $               beta(i)*(shc(i) - shbh(i,k)))/dp(i,k)
604            IF (qtest1.le.0.0) eta(i) = 0.0
605C
606C Check on negative q in lower layer (bound eta)
607C
608            fac1 = -(shbh(i,kp1) - shc(i))/dp(i,kp1)
609            qtest2 = shb(i,kp1) - eta(i)*RG*dt*fac1
610            IF (qtest2 .lt. 0.0) THEN
611               eta(i) = 0.99*shb(i,kp1)/(RG*dt*fac1)
612            ENDIF
613         ENDIF
614         ENDDO
615C
616C
617C Calculate cloud water, rain water, and thermodynamic changes
618C
619         DO 30 i=1,klon
620         IF (ldcum(i)) THEN
621            etagdt = eta(i)*RG*dt
622            cldwtr(i) = etagdt*cldwtr(i)/RLVTT/RG
623            rnwtr(i) = (1.0 - beta(i))*cldwtr(i)
624            ds1(i) = etagdt*(sbh(i,kp1) - sc(i))/dp(i,kp1)
625            dq1(i) = etagdt*(shbh(i,kp1) - shc(i))/dp(i,kp1)
626            ds2(i) = (etagdt*(sc(i) - sbh(i,kp1)) +
627     $                RLVTT*RG*cldwtr(i) - beta(i)*etagdt*
628     $                (sc(i) - sbh(i,k)))/dp(i,k)
629            dq2(i) = (etagdt*(shc(i) - shbh(i,kp1)) -
630     $                RG*rnwtr(i) - beta(i)*etagdt*
631     $                (shc(i) - shbh(i,k)))/dp(i,k)
632            ds3(i) = beta(i)*(etagdt*(sc(i) - sbh(i,k)) -
633     $               RLVTT*RG*cldwtr(i))/dp(i,km1)
634            dq3(i) = beta(i)*etagdt*(shc(i) - shbh(i,k))/dp(i,km1)
635C
636C Isolate convective fluxes for later diagnostics
637C
638            fslkp = eta(i)*(sc(i) - sbh(i,kp1))
639            fslkm = beta(i)*(eta(i)*(sc(i) - sbh(i,k)) -
640     $                       RLVTT*cldwtr(i)*rdt)
641            fqlkp = eta(i)*(shc(i) - shbh(i,kp1))
642            fqlkm = beta(i)*eta(i)*(shc(i) - shbh(i,k))
643C
644C
645C Update thermodynamic profile (update sb, hb, & hbs later)
646C
647            tb (i,kp1) = tb(i,kp1)  + ds1(i) / RCPD
648            tb (i,k  ) = tb(i,k  )  + ds2(i) / RCPD
649            tb (i,km1) = tb(i,km1)  + ds3(i) / RCPD
650            shb(i,kp1) = shb(i,kp1) + dq1(i)
651            shb(i,k  ) = shb(i,k  ) + dq2(i)
652            shb(i,km1) = shb(i,km1) + dq3(i)
653            prec(i)    = prec(i)    + rnwtr(i)/rhoh2o
654C
655C Update diagnostic information for final budget
656C Tracking temperature & specific humidity tendencies,
657C rainout term, convective mass flux, convective liquid
658C water static energy flux, and convective total water flux
659C
660            cmfdt (i,kp1) = cmfdt (i,kp1) + ds1(i)/RCPD*rdt
661            cmfdt (i,k  ) = cmfdt (i,k  ) + ds2(i)/RCPD*rdt
662            cmfdt (i,km1) = cmfdt (i,km1) + ds3(i)/RCPD*rdt
663            cmfdq (i,kp1) = cmfdq (i,kp1) + dq1(i)*rdt
664            cmfdq (i,k  ) = cmfdq (i,k  ) + dq2(i)*rdt
665            cmfdq (i,km1) = cmfdq (i,km1) + dq3(i)*rdt
666            cmfdqr(i,k  ) = cmfdqr(i,k  ) + (RG*rnwtr(i)/dp(i,k))*rdt
667            cmfmc (i,kp1) = cmfmc (i,kp1) + eta(i)
668            cmfmc (i,k  ) = cmfmc (i,k  ) + beta(i)*eta(i)
669            cmfsl (i,kp1) = cmfsl (i,kp1) + fslkp
670            cmfsl (i,k  ) = cmfsl (i,k  ) + fslkm
671            cmflq (i,kp1) = cmflq (i,kp1) + RLVTT*fqlkp
672            cmflq (i,k  ) = cmflq (i,k  ) + RLVTT*fqlkm
673            qc    (i,k  ) =                (RG*rnwtr(i)/dp(i,k))*rdt
674         ENDIF
675   30    CONTINUE
676C
677C Next, convectively modify passive constituents
678C
679         DO 50 m=1,pcnst
680         DO 40 i=1,klon
681         IF (ldcum(i)) THEN
682C
683C If any of the reported values of the constituent is negative in
684C the three adjacent levels, nothing will be done to the profile
685C
686            IF ((cmrb(i,kp1,m).LT.0.0) .OR.
687     $          (cmrb(i,k,m).LT.0.0) .OR.
688     $          (cmrb(i,km1,m).LT.0.0)) GOTO 40
689C
690C Specify constituent interface values (linear interpolation)
691C
692            cmrh(i,k  ) = 0.5*(cmrb(i,km1,m) + cmrb(i,k  ,m))
693            cmrh(i,kp1) = 0.5*(cmrb(i,k  ,m) + cmrb(i,kp1,m))
694C
695C Specify perturbation properties of constituents in PBL
696C
697            pblhgt = MAX(pblh(i),1.0)
698            IF (gz(i,kp1)/RG.LE.pblhgt .AND. dzcld(i).EQ.0.) THEN
699               fac1 = MAX(0.0,1.0-gz(i,kp1)/RG/pblhgt)
700               cmrc(i) = cmrb(i,kp1,m) + cmrp(i,m)*fac1
701            ELSE
702               cmrc(i) = cmrb(i,kp1,m)
703            ENDIF
704C
705C Determine fluxes, flux divergence => changes due to convection
706C Logic must be included to avoid producing negative values. A bit
707C messy since there are no a priori assumptions about profiles.
708C Tendency is modified (reduced) when pending disaster detected.
709C
710            etagdt = eta(i)*RG*dt
711            botflx   = etagdt*(cmrc(i) - cmrh(i,kp1))
712            topflx   = beta(i)*etagdt*(cmrc(i)-cmrh(i,k))
713            dcmr1(i) = -botflx/dp(i,kp1)
714            efac1    = 1.0
715            efac2    = 1.0
716            efac3    = 1.0
717C
718            IF (cmrb(i,kp1,m)+dcmr1(i) .LT. 0.0) THEN
719               efac1 = MAX(tiny,ABS(cmrb(i,kp1,m)/dcmr1(i)) - eps)
720            ENDIF
721C
722            IF (efac1.EQ.tiny .OR. efac1.GT.1.0) efac1 = 0.0
723            dcmr1(i) = -efac1*botflx/dp(i,kp1)
724            dcmr2(i) = (efac1*botflx - topflx)/dp(i,k)
725C 
726            IF (cmrb(i,k,m)+dcmr2(i) .LT. 0.0) THEN
727               efac2 = MAX(tiny,ABS(cmrb(i,k  ,m)/dcmr2(i)) - eps)
728            ENDIF
729C
730            IF (efac2.EQ.tiny .OR. efac2.GT.1.0) efac2 = 0.0
731            dcmr2(i) = (efac1*botflx - efac2*topflx)/dp(i,k)
732            dcmr3(i) = efac2*topflx/dp(i,km1)
733C
734            IF (cmrb(i,km1,m)+dcmr3(i) .LT. 0.0) THEN
735               efac3 = MAX(tiny,ABS(cmrb(i,km1,m)/dcmr3(i)) - eps)
736            ENDIF
737C
738            IF (efac3.EQ.tiny .OR. efac3.GT.1.0) efac3 = 0.0
739            efac3    = MIN(efac2,efac3)
740            dcmr2(i) = (efac1*botflx - efac3*topflx)/dp(i,k)
741            dcmr3(i) = efac3*topflx/dp(i,km1)
742C
743            cmrb(i,kp1,m) = cmrb(i,kp1,m) + dcmr1(i)
744            cmrb(i,k  ,m) = cmrb(i,k  ,m) + dcmr2(i)
745            cmrb(i,km1,m) = cmrb(i,km1,m) + dcmr3(i)
746         ENDIF
747   40    CONTINUE
748   50    CONTINUE              ! end of m=1,pcnst loop
749C
750         IF (k.EQ.limcnv+1) GOTO 60 ! on ne pourra plus glisser
751c
752c Dans la procedure de glissage ascendant, les variables thermo-
753c dynamiques des couches k et km1 servent au calcul des couches
754c superieures. Elles ont donc besoin d'une mise-a-jour.
755C
756         DO i = 1, klon
757         IF (ldcum(i)) THEN
758            zx_t = tb(i,k)
759            zx_p = p(i,k)
760            zx_q = shb(i,k)
761              zdelta=MAX(0.,SIGN(1.,RTT-zx_t))
762              zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
763              zcvm5 = zcvm5 / RCPD / (1.0+RVTMP2*zx_q)
764              zx_qs= r2es * FOEEW(zx_t,zdelta)/zx_p
765              zx_qs=MIN(0.5,zx_qs)
766              zcor=1./(1.-retv*zx_qs)
767              zx_qs=zx_qs*zcor
768              zx_gam = FOEDE(zx_t,zdelta,zcvm5,zx_qs,zcor)
769            shbs(i,k) = zx_qs
770            gam(i,k) = zx_gam
771c
772            zx_t = tb(i,km1)
773            zx_p = p(i,km1)
774            zx_q = shb(i,km1)
775              zdelta=MAX(0.,SIGN(1.,RTT-zx_t))
776              zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
777              zcvm5 = zcvm5 / RCPD / (1.0+RVTMP2*zx_q)
778              zx_qs= r2es * FOEEW(zx_t,zdelta)/zx_p
779              zx_qs=MIN(0.5,zx_qs)
780              zcor=1./(1.-retv*zx_qs)
781              zx_qs=zx_qs*zcor
782              zx_gam = FOEDE(zx_t,zdelta,zcvm5,zx_qs,zcor)
783            shbs(i,km1) = zx_qs
784            gam(i,km1) = zx_gam
785C
786            sb (i,k  ) = sb(i,k  ) + ds2(i)
787            sb (i,km1) = sb(i,km1) + ds3(i)
788            hb (i,k  ) = sb(i,k  ) + RLVTT*shb(i,k)
789            hb (i,km1) = sb(i,km1) + RLVTT*shb(i,km1)
790            hbs(i,k  ) = sb(i,k  ) + RLVTT*shbs(i,k  )
791            hbs(i,km1) = sb(i,km1) + RLVTT*shbs(i,km1)
792C
793            sbh (i,k) = 0.5*(sb(i,k) + sb(i,km1))
794            shbh(i,k) = qhalf(shb(i,km1),shb(i,k)
795     $                       ,shbs(i,km1),shbs(i,k))
796            hbh (i,k) = sbh(i,k) + RLVTT*shbh(i,k)
797            sbh (i,km1) = 0.5*(sb(i,km1) + sb(i,k-2))
798            shbh(i,km1) = qhalf(shb(i,k-2),shb(i,km1),
799     $                    shbs(i,k-2),shbs(i,km1))
800            hbh (i,km1) = sbh(i,km1) + RLVTT*shbh(i,km1)
801         ENDIF
802         ENDDO
803C
804C Ensure that dzcld is reset if convective mass flux zero
805C specify the current vertical extent of the convective activity
806C top of convective layer determined by size of overshoot param.
807C
808   60    CONTINUE
809         DO i=1,klon
810            etagt0 = eta(i).gt.0.0
811            IF (.not.etagt0) dzcld(i) = 0.0
812            IF (etagt0 .and. beta(i).gt.betamn) THEN
813               ktp = km1
814            ELSE
815               ktp = k
816            ENDIF
817            IF (etagt0) THEN
818               cnt(i) = MIN(cnt(i),ktp)
819               cnb(i) = MAX(cnb(i),k)
820            ENDIF
821         ENDDO
822   70 CONTINUE        ! end of k loop
823C
824C determine whether precipitation, prec, is frozen (snow) or not
825C
826      DO i=1,klon
827         IF (tb(i,klev).LT.tmelt .AND. tb(i,klev-1).lt.tmelt) THEN
828             cmfprs(i) = prec(i)*rdt
829         ELSE
830             cmfprt(i) = prec(i)*rdt
831         ENDIF
832      ENDDO
833C
834      RETURN  ! we're all done ... return to calling procedure
835      END
Note: See TracBrowser for help on using the repository browser.