!copyright (C) 2001 MSC-RPN COMM %%%RPNKUO%%%
***S/P KUO2
*
#include "phy_macros_f.h"
SUBROUTINE KUO2 ( TE,QE,CRR,CSR,ILAB,CCK,OMEGAP,CLDW, 2,6
* TP1,TM1,QP1,QM1,GZM1,PSP1,PSM1,KBL,
+ SIGMA, TAU, N, NI, NK, DBGKUO, SATUCO)
#include "impnone.cdk"
*
LOGICAL SATUCO
C
INTEGER N,NK,NI
INTEGER ILAB(NI,NK)
LOGICAL DBGKUO
REAL TAU
REAL TP1(N,NK),TM1(N,NK),QP1(N,NK),QM1(N,NK),GZM1(N,NK)
REAL PSP1(N),PSM1(N),KBL(N),SIGMA(NI,NK)
REAL TE(NI,NK),QE(NI,NK),CRR(NI),CSR(NI)
REAL CCK(NI,NK),OMEGAP(N,NK),CLDW(N,NK)
*
*Author
* J.F.Geleyn E.C.M.W.F. 13/05/82.
*
*Revision
* 001 J. Mailhot R.P.N. 08/02/85. Adaptation and
* modification for EFR model.
* 002 J. Mailhot (Mar 1987) base of condensation
* 003 G.Pellerin (Oct87)Adaptation to revised code
* 004 J. Mailhot (Mar 1988) threshold of evaporation
* 005 G.Pellerin(August90)Adaptation to thermo functions
* 006 C. Girard(Nov 90)
* Remove calculations related to the boundary layer
* Introduction of the entrainment for the cloud profile
* Modification of the heating/moistening partition:NFRAC**3
* 007 N. Brunet (May91)
* New version of thermodynamic functions
* and file of constants
* 008 B. Bilodeau (August 1991)- Adaptation to UNIX
* 009 J. Mailhot (Dec 1992) - Bug correction to the
* Newton Method
* 010 C. Girard (Nov92) - Clean-up, and implicit
* calculation of evaporation and precipitation.
* Recycling of evaporated precipitation for
* modelling the increase in humidity and calculation
* of cloud fraction
* 011 A. Methot (Sept 93) -
* -Bug correction : L/Cp missing in precip. evap.
* -Evap. set to zero : logical variable EVAP skips unnecessary
* calculations when no evaporation of precip.
* -Change heating/moistening partition :
* BETA=((1-RH)/(1-RHC))**3;(1-RHC)**-3=4;;RHC~.37
* -SATUCO set to .FALSE. locally by the use of SATUC.
* Change the value of CHLS depending on SATUC.
* 012 G. Lemay (Oct 93) - Dynamic memory allocation with stkmemw
* 013 R. Benoit (Dec 93) - Restore the ILAB output for use by
* Sundqvist scheme. Also use icvmgt to handle integer ILAB.
* 014 A. Methot (Dec 93) - SATUCO is back as before revision 011
* but the value of CHLS is determined by SATUCO
* Also ZTSATCO determine whether CHLS or CHLC is used
* -criteria on vertical velocity (OMEGAP)at SIGMAK1 SIGMAK2
* 015 B. Bilodeau (Feb 94) - Cleanup - Change name from KUO to KUO2
* 018 R. Sarrazin (June 95) evap, beta, liquid water
* 016 B. Bilodeau (June 94) - New physics interface
* 017 G. Pellerin (Jan 94) - Omitted bugfix for ZTSATCO
* 018 R. Sarrazin (Summer 95) - Corrections; add diagnostic cloud water.
* 019 B. Bilodeau (Jan 2001) - Automatic arrays
* 020 M. Lepine (March 2003) - CVMG... Replacements
* 021 G. Pellerin (Mai 03) - Conversion IBM
* - calls to optimized routine MFOQST
*
*Object
* to calculate the tendencies of T and Q by moist convective
* adjustment according to KUO equations
*
*Arguments
*
* - Output -
* TE temperature tendency
* QE specific humidity tendency
* CRR rate of liquid precipitation
* CSR rate of solid precipitation (not available)
* ILAB label array from Kuo scheme
*
* - Output -
* CCK cloud fraction due to the Kuo scheme
*
* - Input -
* OMEGAP vertical velocity in pressure coordinates
*
* - Output -
* CLDW diagnostic cloud water
*
* - Input/Output -
* TP1 temperature at (T+DT)
* TM1 temperature at (T-DT)
* QP1 specific humidity at (T+DT)
* QM1 specific humidity at (T-DT)
*
* - Input -
* GZM1 geopotential (N,NK)
* PSP1 surface pressure at (T+DT)
* PSM1 surface pressure at (T-DT)
* KBL index of 1st level of the boundary layer
* SIGMA sigma levels
* TAU FACTDT * timestep (see common block OPTIONS)
* N dimension NI*NJ
* NI 1st horizontal dimension
* NK vertical dimension
* DBGKUO .TRUE. to debug KUO routine
* .FALSE. to not debug
* SATUCO .TRUE. if water/ice phase for saturation
* .FALSE. if water phase only for saturation
*
*Notes
* The process accounts for the moisture convergence for the
* formation of cumulus clouds, the humidification of the
* environment and the formation of precipitation. There is
* an exact balance between these 3 features. The process
* adapts the T values at (T+DT) and Q values at (T-DT) (the
* difference between T-DT and T+DT will give the moisture
* convergence) and the proportional intensity to available
* convergence. The routine is divided into 5 different
* steps:
* 1)allocation and position for work space
* 2)preliminary computations
* 3)cloud ascent and flagging
* 4)total moisture convergence and mean beta-parameter
* 5)moistening, condensation and evaporation of rain/snow
*
*
*REFERENCE
* ECMWF RESEARCH MANUAL (VOLUME 3)
* PHYSICAL PARAMETRIZATION CHAPITRE 4
**
*
real wc0, wmr, pc1, pc2
*
parameter (wc0=5.0E-4)
parameter (wmr=5.0E-4)
parameter (pc1=300.)
parameter (pc2=0.5)
*
LOGICAL LO,EVAP
INTEGER IERR
INTEGER NKP1,NKM1
INTEGER IS,IKS,IKS2
INTEGER JK,JL,MODP
INTEGER HEAPERR
REAL TEMP1, TEMP2, TEMP3, temp4, temp5, temp6
*
*
*
************************************************************************
* AUTOMATIC ARRAYS
************************************************************************
*
AUTOMATIC (LO1 , LOGICAL, (NI ))
*
AUTOMATIC (IQCD , INTEGER, (NI ))
AUTOMATIC (KM1 , INTEGER, (NI ))
AUTOMATIC (KM2 , INTEGER, (NI ))
AUTOMATIC (KP1 , INTEGER, (NI ))
AUTOMATIC (KP2 , INTEGER, (NI ))
*
AUTOMATIC (ZCDP , REAL , (NI ))
AUTOMATIC (ZCPD , REAL , (NI ))
AUTOMATIC (ZRFL , REAL , (NI ))
AUTOMATIC (ZSFL , REAL , (NI ))
AUTOMATIC (ZCUCOV , REAL , (NI ))
AUTOMATIC (ZRFLN , REAL , (NI ))
AUTOMATIC (ZSFLN , REAL , (NI ))
AUTOMATIC (SIGMAK1, REAL , (NI ))
AUTOMATIC (SIGMAK2, REAL , (NI ))
AUTOMATIC (KWW1 , REAL , (NI ))
AUTOMATIC (KWW2 , REAL , (NI ))
AUTOMATIC (ZQCD , REAL , (NI ))
AUTOMATIC (ZQSATC , REAL , (NI,NK))
AUTOMATIC (ZLDCPE , REAL , (NI,NK))
AUTOMATIC (ZPP1 , REAL , (NI,NK))
AUTOMATIC (ZDSG , REAL , (NI,NK))
AUTOMATIC (ZDPP1 , REAL , (NI,NK))
AUTOMATIC (ZQAC , REAL , (NI,NK))
AUTOMATIC (ZTP1 , REAL , (NI,NK))
AUTOMATIC (ZQP1 , REAL , (NI,NK))
AUTOMATIC (ZQSATE , REAL , (NI,NK))
AUTOMATIC (ZBETA , REAL , (NI,NK))
AUTOMATIC (ZDQTOT , REAL , (NI,NK))
AUTOMATIC (ZTC , REAL , (NI,NK))
AUTOMATIC (ZQC , REAL , (NI,NK))
AUTOMATIC (ZTVP1 , REAL , (NI,NK))
AUTOMATIC (ZDQLOC , REAL , (NI,NK))
AUTOMATIC (ZDTLOC , REAL , (NI,NK))
AUTOMATIC (ZDTTOT , REAL , (NI,NK))
*
************************************************************************
*
REAL ZEVAP, ZMELT, ZCCTIM, ZRCYCL, ZEPFLM, ZEPFLS
REAL ZEPCDP, ZEPCOV, ZTMST, ZCONS2, ZCONS3
REAL ZCONS5, ZDPM1, ZLDCP0, ZCOR
REAL ZDP, ZDQK, ZCUPRO, ZCPDL, ZDTK, ZSQFLN, ZNIMP
REAL ZRFLN0, ZLVDCP, ENTRM
* real cevap, ccctim, cmelt
*
C* PHYSICAL CONSTANTS.
C -------- ----------
C
#include "comphy.cdk"
*
REAL CHLS, DELTA2
#include "consphy.cdk"
#include "dintern.cdk"
#include "fintern.cdk"
*
C *NSTAB* REFERS TO THE MINIMUM NUMBER OF STABLE LAYERS BETWEEN THE
C TOP AND THE CLOUD BASE.
C IF *EVAP* IS TRUE: THE EVAPORATION OF PRECIPITATION IS ACTIVATED
C OTHERWISE THERE IS NO EVAPORATION
C *ZEVAP* IS A CONSTANT FOR THE EVAPORATION OF PRECIPITATION
C *ZCCTIM* IS THE CUMULUS LIFE-TIME VALUE TO COMPUTE ITS CLOUD COVER
C *ZRCYCL* IS THE EVAPORATED PRECIP RECYCLED FRACTION
C
NKP1=NK+1
NKM1=NK-1
*
EVAP=.TRUE.
IF ( EVAP ) THEN
ZEVAP = 6.0E-05
ELSE
ZEVAP = 0.0
ENDIF
ZMELT = CMELT
ZCCTIM = CCCTIM
ZRCYCL = 1.0
ENTRM = 5.E-6
C
C
C* SECURITY PARAMETERS.
C --------------------
C
C *ZEPFLM* IS A MINIMUM FLUX TO AVOID DIVIDING BY ZERO IN THE IC
C PROPORTION CALCULATIONS, *ZEPCDP* AVOIDS DIVIDING BY ZERO IN THE
C ABSENCE OF CLOUD AND *ZEPCOV* IS A MINIMUM CLOUD COVER.
C
ZEPFLM=1.E-24
ZEPFLS=SQRT(ZEPFLM)
ZEPCDP=1.E-12
ZEPCOV=1.E-12
C
C* COMPUTATIONAL CONSTANTS.
C ------------- ----------
C
ZTMST = TAU
DELTA2 = CPV/CPD - 1.
IF (SATUCO) THEN
CHLS = CHLC + CHLF
ELSE
CHLS = CHLC
ENDIF
C
ZCONS2 = 1./(ZTMST*GRAV)
ZCONS3 = ZCCTIM/ZTMST
ZCONS5 = 1./ZTMST
C
C
C ------------------------------------------------------------------
C
C* 1. ALLOCATION AND POSITION FOR WORK SPACE.
C ---------- --- -------- --- ---- ------
C
100 CONTINUE
C
C
*
*
C***
C
C METHOD.
C -------
C
C IN ORDER TO KEEP THE CODE LINEAR THE ROUTINE HAS BEEN DIVIDED
C INTO THREE PARTS THAT WE SHALL CALL HERE (A) (B) AND (C).
C IN (A) PRELIMINARY COMPUTATIONS ARE FIRST PERFORMED. THEN A
C NEARLY ADIABATIC ASCENT IS ATTEMPTED FOR A CLOUD PARCEL STARTING
C FROM THE LOWEST MODEL LAYER. THIS CLOUD ASCENT IS COMPUTED
C IN TERMS OF TEMPERATURE AND SPECIFIC HUMIDITY.
C ENTRAINMENT IS SIMULATED VIA THE ENTRAINMENT PARAMETER.
C THE LEVELS ARE FLAGGED ACCORDING TO THE FOLLOWING CODE:
C 0 = STABLE,
C 1 = PART OF THE WELL MIXED BOUNDARY LAYER OR DRY UNSTABLE SLAB,
C 2 = MOIST UNSTABLE SLAB (I.E. CLOUD LAYER) AND
C 3 = LIFTING LEVEL IF DIFFERENT FROM THE GROUND.
C IN (B) THE TOTAL MOISTURE CONVERGENCE FOR EACH SLAB OF
C NON-0-FLAGS IS STORED INTO ALL THE CORRESPONDING LAYERS IF IT IS
C POSITIVE AND THE SAME IS DONE FOR THE MEAN OF 1-RELATIVE HUMIDITY
C FOR EACH CORRESPONDING SLAB OF 2-FLAGS.
C IN (C) THE ACTUAL MODIFICATIONS OF TEMPERATURE AND SPECIFIC
C HUMIDITY ARE COMPUTED. FIRST THE ENVIRONMENTAL MOISTENING IS
C TAKEN PROPORTIONAL TO THE MOISTURE CONVERGENCE WEIGHTED BY "BETA"
C AND TO THE SATURATION DEFICIT OF SPECIFIC HUMIDITY. AT THE SAME
C TIME THE MOISTURE CONVERGENCE IS TAKEN AWAY FROM THE HUMIDITY FIELD.
C SECOND THE FORMATION OF PRECIPITATION IS TAKEN PROPORTIONAL TO THE
C MOISTURE CONVERGENCE WEIGHTED BY "1 - BETA" AND TO THE
C TEMPERATURE DIFFERENCE BETWEEN THE CLOUD AND THE ENVIRONMENT.
C "BETA" , THE MOISTENING PARAMETER, IS A FUNCTION OF MEAN REL. HUM.
C A CLOUD-COVER VALUE IS OBTAINED BY COMPARING THE TIME AT WHICH THE
C ENVIRONMENT WOULD REACH EQUILIBRIUM WITH THE CLOUD TO A PRESCRIBED
C LIFE-TIME VALUE FOR THE CLOUD ITSELF.
C
C
C ------------------------------------------------------------------
C
C* 2. PRELIMINARY COMPUTATIONS.
C ----------- -------------
C
200 CONTINUE
C
C* 2.1 NECESSARY SETTING FOR THE CASE THERE WOULD NOT BE
C* ANY CONVECTIVE POINT AROUND THE LATITUDE CIRCLE.
C
210 CONTINUE
C
C* 2.2 MOISTURE ACCESSION, T+1 T,Q VARIABLES AND SATURATION
C* MIXING RATIO PLUS INITIAL VALUE FOR TEST FLAG, ETC.
C
C The moisture accession is set to zero for all levels when the
C vertical velocity in pressure coordinates OMEGAP is positive
C (downward) at both sigma levels SIGMAK1 and SIGMAK2.
C
C Since SIGMAK1 and 2 do not necessarly coincide with sigma levels
C of the model, OMEGAP is linearly interpolated at levels SIGMAK1 and 2
C using weighting factors KWW1 and KWW2. KP1 or KP2 and KM1 or KM2 are
C the indices of the model's sigma levels from which the interpolation
C takes place.
*
C SIGMAK1 and SIGMAK2 are the sigma levels at which OMEGAP is tested
DO 10 JL = 1,NI
SIGMAK1(JL) = 0.9
SIGMAK2(JL) = 0.7
KM1 (JL) = NK
KM2 (JL) = NK
10 CONTINUE
*
* general case
DO 20 JK = 1,NK
*
DO 30 JL = 1,NI
*
IF (SIGMA(JL,JK) .LE. SIGMAK1(JL)) KM1(JL) = JK
IF (SIGMA(JL,JK) .LE. SIGMAK2(JL)) KM2(JL) = JK
*
KP1(JL) = KM1(JL) + 1
KP2(JL) = KM2(JL) + 1
*
KWW1(JL)=(SIGMAK1(JL)-SIGMA(JL,KM1(JL)))/
+ (SIGMA(JL,KP1(JL))-SIGMA(JL,KM1(JL)))
*
KWW2(JL)=(SIGMAK2(JL)-SIGMA(JL,KM2(JL)))/
+ (SIGMA(JL,KP2(JL))-SIGMA(JL,KM2(JL)))
*
30 CONTINUE
*
20 CONTINUE
*
* loop 40 covers special cases
DO 40 JL = 1,NI
*
IF (SIGMA(JL,1).GT.SIGMAK1(JL)) THEN
SIGMAK1(JL) = SIGMA(JL,1)
KM1(JL) = 1
KP1(JL) = 1
KWW1(JL) = 1.0
ENDIF
IF (SIGMA(JL,1).GT.SIGMAK2(JL)) THEN
SIGMAK2(JL) = SIGMA(JL,1)
KM2(JL) = 1
KP2(JL) = 1
KWW2(JL) = 1.0
ENDIF
IF ( KM1(JL) .EQ. NK ) THEN
KP1(JL) = KM1(JL)
KWW1(JL) = 1.0
ENDIF
IF ( KM2(JL) .EQ. NK ) THEN
KP2(JL) = KM2(JL)
KWW2(JL) = 1.0
ENDIF
*
40 CONTINUE
*
*
do 220 jl=1,ni
ZDSG(jl,1)=0.5*(SIGMA(jl,2)-SIGMA(jl,1))
DO 225 JK=2,NKM1
ZDSG(jl,JK)=0.5*(SIGMA(jl,JK+1)-SIGMA(jl,JK-1))
225 CONTINUE
ZDSG(jl,NK)=0.5*(1.-SIGMA(jl,NKM1))+0.5*(1.-SIGMA(jl,NK))
220 continue
C
DO 221 JK=1,NK
DO 221 JL=1,NI
ZPP1(JL,JK)=SIGMA(jl,JK)*PSP1(JL)
ZTP1(JL,JK)=TP1(JL,JK)
ZQP1(JL,JK)=QP1(JL,JK)
ZDPP1(JL,JK)=ZDSG(jl,JK)*PSP1(JL)
ZDPM1 =ZDSG(jl,JK)*PSM1(JL)
ZQAC(JL,JK)=ZQP1(JL,JK)*ZDPP1(JL,JK)-QM1(JL,JK)*ZDPM1
LO=(OMEGAP(JL,KP1(JL))*KWW1(JL) +
+ OMEGAP(JL,KM1(JL))*(1-KWW1(JL))).GT.0.
LO=LO .AND. (OMEGAP(JL,KP2(JL))*KWW2(JL)+
+ OMEGAP(JL,KM2(JL))*(1-KWW2(JL))).GT.0.
if (LO) ZQAC(JL,JK)=0.
ILAB(JL,JK)=0
ZBETA(JL,JK)=0.
ZDQTOT(JL,JK)=-1.
221 CONTINUE
C
MODP=3
IF(SATUCO)THEN
CALL MFOQST(ZQSATE,ZTP1,SIGMA,ZPP1,MODP,NI,NK,NI)
ELSE
CALL MFOQSA(ZQSATE,ZTP1,SIGMA,ZPP1,MODP,NI,NK,NI)
ENDIF
C
C* 2.3 COMPUTATIONS WITHIN THE BOUNDARY LAYER.
C
230 CONTINUE
C
C* 2.4 SPECIFY TC AND QC AT THE LOWEST LAYER TO START THE
C* CLOUD ASCENT. CHECK FOR POSITIVE MOISTURE ACCESSION
C* BETWEEN SURFACE AND CLOUD BASE.
C* ZQC=0 INDICATES STABLE CONDITIONS.
C
240 CONTINUE
DO 241 JL=1,NI
LO=ZQAC(JL,NK).GT.0.
ZTC(JL,NK)=ZTP1(JL,NK)
if (LO) then
ZQC(JL,NK) = QM1(JL,NK)
else
ZQC(JL,NK) = 0.
endif
IF (LO) ILAB(JL,NK) = 1
241 CONTINUE
C
C
C ------------------------------------------------------------------
C
C* 3. CLOUD ASCENT AND FLAGGING.
C ----- ------ --- ---------
C
300 CONTINUE
C
C* 3.1 CALCULATE TC AND QC AT NEXT LEVEL BY DRY ADIABATIC
C* LIFTING AND CONSIDERING LATENT HEAT RELEASE (THE
C* GEOPOTENTIAL DIFFERENCE IS CORRECTED TO TAKE INTO
C* ACCOUNT THE VIRTUAL TEMPERATURE DIFFERENCE BETWEEN
C* THE ASCENT AND THE ENVIRONMENT). THE CONDENSATION
C* CALCULATIONS ARE DONE WITH TWO ITERATIONS.
C
310 CONTINUE
C
DO 311 JL=1,NI
ZCPD(JL)=CPD*(1.+DELTA2*ZQC(JL,NK))
ZTVP1(JL,NK) = FOTVT(ZTP1(JL,NK),QM1(JL,NK))
311 CONTINUE
C***
DO 322 JK=NKM1,1,-1
C***
DO 312 JL=1,NI
ZTC(JL,JK)=ZTC(JL,JK+1)+(GZM1(JL,JK+1)-GZM1(JL,JK))*
* (1./ZCPD(JL)+ENTRM*MAX(0.,ZTC(JL,JK+1)-ZTP1(JL,JK+1)))
ZQC(JL,JK)=ZQC(JL,JK+1)+(GZM1(JL,JK+1)-GZM1(JL,JK))*
* ( ENTRM*MAX(0.,ZQC(JL,JK+1)-QM1(JL,JK+1)))
TEMP1 = ZTP1(JL,JK)
ZTVP1(JL,JK) = FOTVT(ZTP1(JL,JK),QM1(JL,JK))
LO=(FOTVT(ZTC(JL,JK),ZQC(JL,JK)).GT.ZTVP1(JL,JK)).AND
* .(ZQC(JL,JK).NE.0.)
IF (LO) ILAB(JL,JK) = 1
LO = ZTC(JL,JK).LT.TRPL .and. SATUCO
! ZLDCPE = CVMGT(CHLS,CHLC,LO)/ZCPD(JL)
if (LO) then
ZLDCPE(jl,jk) = CHLS/ZCPD(JL)
else
ZLDCPE(jl,jk) = CHLC/ZCPD(JL)
endif
312 CONTINUE
IF(SATUCO)THEN
CALL MFOQST(ZQSATC(1,jk),ZTC(1,jk),SIGMA,ZPP1(1,jk),MODP,NI,1,NI)
DO 313 JL=1,NI
* ZQSATC(jl,jk)=FOQST(ZTC(JL,JK),ZPP1(JL,JK))
ZCOR=ZLDCPE(jl,jk)*FODQS(ZQSATC(jl,jk),ZTC(JL,JK))
ZQCD(jl)=AMAX1(0.,(ZQC(JL,JK)-ZQSATC(jl,jk))/(1.+ZCOR))
313 CONTINUE
ELSE
CALL MFOQSA(ZQSATC(1,jk),ZTC(1,jk),SIGMA,ZPP1(1,jk),MODP,NI,1,NI)
DO 317 JL=1,NI
* ZQSATC(jl,jk)=FOQSA(ZTC(JL,JK),ZPP1(JL,JK))
ZCOR=ZLDCPE(jl,jk)*FODQA(ZQSATC(jl,jk),ZTC(JL,JK))
ZQCD(jl)=AMAX1(0.,(ZQC(JL,JK)-ZQSATC(jl,jk))/(1.+ZCOR))
317 CONTINUE
ENDIF
DO JL=1,NI
LO=ZQCD(jl).EQ.0.
IQCD(JL) = 1
IF (LO) IQCD(JL) = 0
ZQC(JL,JK)=ZQC(JL,JK)-ZQCD(jl)
ZTC(JL,JK)=ZTC(JL,JK)+ZQCD(jl)*ZLDCPE(jl,jk)
LO1(JL)=.FALSE.
LO = ZTC(JL,JK).LT.TRPL
if (LO) then
ZLDCPE(jl,jk) = CHLS/ZCPD(JL)
else
ZLDCPE(jl,jk) = CHLC/ZCPD(JL)
endif
ENDDO
IS=0
DO 314 JL=1,NI
IS=IS+IQCD(JL)
314 CONTINUE
IF (IS.NE.0) THEN
IF(SATUCO)THEN
CALL MFOQST(ZQSATC(1,jk),ZTC(1,jk),SIGMA,ZPP1(1,jk),MODP,NI,1,NI)
DO 315 JL=1,NI
* ZQSATC(jl,jk)=FOQST(ZTC(JL,JK),ZPP1(JL,JK))
ZCOR=ZLDCPE(jl,jk)*FODQS(ZQSATC(jl,jk),ZTC(JL,JK))
ZQCD(jl)=(ZQC(JL,JK)-ZQSATC(jl,jk))/(1.+ZCOR)
315 CONTINUE
ELSE
CALL MFOQSA(ZQSATC(1,jk),ZTC(1,jk),SIGMA,ZPP1(1,jk),MODP,NI,1,NI)
DO 318 JL=1,NI
* ZQSATC(jl,jk)=FOQSA(ZTC(JL,JK),ZPP1(JL,JK))
ZCOR=ZLDCPE(jl,jk)*FODQA(ZQSATC(jl,jk),ZTC(JL,JK))
ZQCD(jl)=(ZQC(JL,JK)-ZQSATC(jl,jk))/(1.+ZCOR)
318 CONTINUE
ENDIF
DO JL=1,NI
LO1(JL)=IQCD(JL).NE.0
if (.not. LO1(JL)) ZQCD(jl) = 0.
ZQC(JL,JK)=ZQC(JL,JK)-ZQCD(jl)
ZTC(JL,JK)=ZTC(JL,JK)+ZQCD(jl)*ZLDCPE(jl,jk)
ENDDO
ENDIF
DO 316 JL=1,NI
LO=(FOTVT(ZTC(JL,JK),ZQC(JL,JK)).GT.ZTVP1(JL,JK))
* .AND.LO1(JL)
IF (LO) ILAB(JL,JK) = 2
LO1(JL)=ILAB(JL,JK).EQ.0
if (LO1(JL)) ZTC(JL,JK) = ZTP1(JL,JK)
if (LO1(JL)) ZQC(JL,JK) = 0.
316 CONTINUE
C
C
C* 3.2 IF NOT AT THE TOP CHECK FOR NEW LIFTING LEVEL, I.E.
C* MOISTURE CONVERGENCE IN A STABLE LAYER.
C
C***
IF (JK.NE.1) THEN
C***
DO 321 JL=1,NI
LO=LO1(JL).AND.(ZQAC(JL,JK).GT.0.)
if (LO) ZTC(JL,JK) = ZTP1(JL,JK)
if (LO) ZQC(JL,JK) = QM1(JL,JK)
ZCPD(JL)=CPD*(1.+DELTA2*ZQC(JL,JK))
321 CONTINUE
C***
ENDIF
322 CONTINUE
C
C***
C
C* 3.3 ILAB=0 FOR DRY UNSTABLE LAYERS IF NO CLOUD IS ABOVE
C* ILAB=3 FOR LIFTING LEVEL IF LAYER ABOVE IS UNSTABLE.
C* IKS2 INDICATES THE HIGHEST TOP OF A CLOUD AROUND THE
C* LATITUDE CIRCLE (TO AVOID UNNECESSARY COMPUTATIONS
C* LATER).
C
330 CONTINUE
IKS2=NKP1
DO 331 JL=1,NI
LO=ILAB(JL,1).EQ.1
IF (LO) ILAB(JL,1) = 0
331 CONTINUE
IS=0
DO 332 JL=1,NI
IS=IS+ILAB(JL,1)
332 CONTINUE
IF (IS.NE.0) IKS2=1
DO 335 JK=2,NK
DO 333 JL=1,NI
LO=(ILAB(JL,JK).EQ.1).AND.(ILAB(JL,JK-1).EQ.0)
IF (LO) ILAB(JL,JK) = 0
333 CONTINUE
IF (IKS2.EQ.NKP1) THEN
IS=0
DO 334 JL=1,NI
IS=IS+ILAB(JL,JK)
334 CONTINUE
IF (IS.NE.0) IKS2=JK
ENDIF
335 CONTINUE
C***
IF (IKS2.EQ.NKP1) GO TO 600
C***
DO 337 JK=NK,2,-1
DO 336 JL=1,NI
LO=(ILAB(JL,JK).EQ.0).AND.(ILAB(JL,JK-1).NE.0)
IF (LO) ILAB(JL,JK) = 3
336 CONTINUE
337 CONTINUE
C
C
C ------------------------------------------------------------------
C
C* 4. TOTAL MOISTURE CONVERGENCE AND MEAN BETA-PARAMETER.
C ----- -------- ----------- --- ---- ---------------
C
400 CONTINUE
C
C* 4.1 CALCULATE TOTAL MOISTURE ACCESSION FOR UNSTABLE
C* LAYERS AND THE PARTITION PARAMETER BETA
C* AVERAGED OVER CLOUD LAYERS.
C
410 CONTINUE
DO 411 JL=1,NI
LO=ILAB(JL,NK).GT.0
if (.not. LO) ZQAC(JL,NK) = 0.
ZQSATE(JL,NK) = AMAX1(ZQSATE(JL,NK),QM1(JL,NK))
ZBETA(JL,NK)=0.
ZCDP(JL)=0.
411 CONTINUE
DO 413 JK=NKM1,IKS2,-1
DO 412 JL=1,NI
LO=ILAB(JL,JK).GT.0
if (.not. LO) ZQAC(JL,JK) = 0.
LO=LO.AND.(ILAB(JL,JK).NE.3)
if (LO) ZQAC(JL,JK) = ZQAC(JL,JK+1)+ZQAC(JL,JK)
ZQSATE(JL,JK) = AMAX1(ZQSATE(JL,JK),QM1(JL,JK))
ZBETA(JL,JK) = AMAX1(0.,QM1(JL,JK))/ZQSATE(JL,JK)
LO=ILAB(JL,JK).EQ.2
if (.not. LO) ZBETA(JL,JK) = 0.
if (LO) then
ZDP = ZDPP1(JL,JK)
else
ZDP = 0.
endif
LO=LO.AND.(ILAB(JL,JK+1).EQ.2)
if (LO) then
ZBETA(JL,JK) = (ZCDP(JL)*ZBETA(JL,JK+1)+ZDP*ZBETA(JL,JK))
% /AMAX1(ZCDP(JL)+ZDP,ZEPCDP)
endif
if (LO) then
ZCDP(JL) = ZCDP(JL)+ZDP
else
ZCDP(JL) = ZDP
endif
412 CONTINUE
413 CONTINUE
C
C
C* 4.2 REPLACE THE MOISTURE ACCESSION AT CLOUD LAYERS BY THE
C* TOTAL MOISTURE ACCESSION OF THE WHOLE UNSTABLE LAYER
C* AND DO THE SAME FOR THE BETA-PARAMETER MEAN VALUE.
C* UPDATE IKS2 DURING THE PROCESS.
C
420 CONTINUE
IKS2=NKP1
DO 421 JL=1,NI
LO=ZQAC(JL,1).GT.0.
IF (.NOT.LO) ILAB(JL,1) = 0
421 CONTINUE
IS=0
DO 422 JL=1,NI
IS=IS+ILAB(JL,1)
422 CONTINUE
IF (IS.NE.0) IKS2=1
DO 425 JK=2,NK
DO 423 JL=1,NI
LO=(ILAB(JL,JK).EQ.2).AND.(ILAB(JL,JK-1).EQ.2)
if (LO) ZQAC(JL,JK) = ZQAC(JL,JK-1)
temp1=(1-ZBETA(JL,JK))*(1-ZBETA(JL,JK))
temp1=temp1*(1-ZBETA(JL,JK))
ZBETA(JL,JK) = MIN( 1.0, 6.0*temp1 )
if (LO) ZBETA(JL,JK) = ZBETA(JL,JK-1)
LO=(ZQAC(JL,JK).LE.0.).AND.(ILAB(JL,JK).EQ.2)
IF (LO) ILAB(JL,JK) = 0
LO=(ILAB(JL,JK).NE.2).AND.(ILAB(JL,JK-1).EQ.0)
IF (LO) ILAB(JL,JK) = 0
423 CONTINUE
IF (IKS2.EQ.NKP1) THEN
IS=0
DO 424 JL=1,NI
IS=IS+ILAB(JL,JK)
424 CONTINUE
IF (IS.NE.0) IKS2=JK
ENDIF
425 CONTINUE
C
C
C***
IF (IKS2.EQ.NKP1) GO TO 600
C***
C
C ------------------------------------------------------------------
C
C* 5. MOISTENING, CONDENSATION AND EVAPORATION OF RAIN/SNOW.
C ----------- ------------ --- ----------- -- ----------
C
500 CONTINUE
C
C* 5.1 COMPUTE THE TOTAL MOISTURE DEFICIT IN THE CLOUD
C* LAYERS.
C
510 CONTINUE
DO 511 JL=1,NI
ZDQLOC(JL,NK)=0.
ZDQTOT(JL,NK)=0.
511 CONTINUE
DO 513 JK=NKM1,IKS2,-1
DO 512 JL=1,NI
ZDQLOC(JL,JK) = ZQSATE(JL,JK)-QM1(JL,JK)
ZDQK = ZDQLOC(JL,JK)*ZDPP1(JL,JK)
LO=ILAB(JL,JK).EQ.2
if (LO) then
ZDQTOT(JL,JK) = ZDQK
else
ZDQTOT(JL,JK) = 0.
endif
LO=LO.AND.(ILAB(JL,JK+1).EQ.2)
if (LO) ZDQTOT(JL,JK) = ZDQTOT(JL,JK+1)+ZDQK
512 CONTINUE
513 CONTINUE
IF (IKS2.EQ.1) THEN
DO 514 JL=1,NI
LO=ZDQTOT(JL,1).GT.0.
if (.not. LO) ZDQTOT(JL,1) = -1.
514 CONTINUE
ENDIF
IKS=MAX0(2,IKS2)
DO 516 JK=IKS,NK
DO 515 JL=1,NI
LO=(ILAB(JL,JK).EQ.2).AND.(ILAB(JL,JK-1).EQ.2)
if (LO) ZDQTOT(JL,JK) = ZDQTOT(JL,JK-1)
LO=ZDQTOT(JL,JK).GT.0.
if (.not. LO) ZDQTOT(JL,JK) = -1.
515 CONTINUE
516 CONTINUE
C
C
C* 5.2 REMOVE THE MOISTURE ACCESSION IN THE RELEVANT LAYERS
C* (BY RETURNING TO THE PREVIOUS TIMESTEP VALUES)
C* AND DO THE ENVIRONMENTAL MOISTENING.
C
520 CONTINUE
DO 523 JK=IKS2,NK
DO 522 JL=1,NI
LO=ILAB(JL,JK).GT.0
ZDPM1 = ZDSG(jl,JK)*PSM1(JL)
if (LO) ZQP1(JL,JK) = QM1(JL,JK)*ZDPM1/ZDPP1(JL,JK)
ZCUPRO = ZDQLOC(JL,JK)*ZBETA(JL,JK)*ZQAC(JL,JK)/ZDQTOT(JL,JK)
LO=ILAB(JL,JK).EQ.2
if (.not. LO) ZCUPRO = 0.
ZQP1(JL,JK) = ZQP1(JL,JK)+ZCUPRO
522 CONTINUE
523 CONTINUE
C
C
C* 5.3 COMPUTATION OF THE TOTAL VIRTUAL TEMPERATURE DEFICIT
C* IN CLOUD LAYERS.
C
530 CONTINUE
DO 531 JL=1,NI
ZTC(JL,NK)=ZTP1(JL,NK)
ZQC(JL,NK)=QM1(JL,NK)
ZDTLOC(JL,NK)=0.
ZDTTOT(JL,NK)=0.
531 CONTINUE
DO 533 JK=NKM1,IKS2,-1
DO 532 JL=1,NI
LO = ZTC(JL,JK).LT.TRPL
! ZCPDL=CPD*(1.+DELTA2*QM1(JL,JK))/(CVMGT(CHLS,CHLC,LO)*(1.+DELTA
! * *QM1(JL,JK)))
if (LO) then
ZCPDL=CPD*(1.+DELTA2*QM1(JL,JK))/(CHLS*(1.+DELTA*QM1(JL,JK)))
else
ZCPDL=CPD*(1.+DELTA2*QM1(JL,JK))/(CHLC*(1.+DELTA*QM1(JL,JK)))
endif
LO=ILAB(JL,JK).EQ.2
! ZTC(JL,JK) = CVMGT(ZTC(JL,JK),ZTP1(JL,JK),LO)
if (.not. LO) ZTC(JL,JK) = ZTP1(JL,JK)
! ZQC(JL,JK) = CVMGT(ZQC(JL,JK),QM1(JL,JK),LO)
if (.not. LO) ZQC(JL,JK) = QM1(JL,JK)
! ZQSATE(JL,JK) = CVMGT(ZQC(JL,JK),ZQSATE(JL,JK),LO)
if (LO) ZQSATE(JL,JK) = ZQC(JL,JK)
ZDTLOC(JL,JK)=(FOTVT(ZTC(JL,JK),ZQC(JL,JK))-ZTVP1(JL,JK))
* *ZCPDL
ZDTK = ZDTLOC(JL,JK)*ZDPP1(JL,JK)
! ZDTTOT(JL,JK) = CVMGT(ZDTK,0.,LO)
if (LO) then
ZDTTOT(JL,JK) = ZDTK
else
ZDTTOT(JL,JK) = 0.
endif
LO=LO.AND.(ILAB(JL,JK+1).EQ.2)
if (LO) ZDTTOT(JL,JK) = ZDTTOT(JL,JK+1)+ZDTK
532 CONTINUE
533 CONTINUE
IF (IKS2.EQ.1) THEN
DO 534 JL=1,NI
LO=ZDTTOT(JL,1).GT.0.
if (.not. LO) ZDTTOT(JL,1) = -1.
534 CONTINUE
ENDIF
DO 536 JK=IKS,NK
DO 535 JL=1,NI
LO=(ILAB(JL,JK).EQ.2).AND.(ILAB(JL,JK-1).EQ.2)
if (LO) ZDTTOT(JL,JK) = ZDTTOT(JL,JK-1)
LO=ZDTTOT(JL,JK).GT.0.
if (.not. LO) ZDTTOT(JL,JK) = -1.
535 CONTINUE
536 CONTINUE
C
C
C* 5.4 FORMATION OF PRECIPITATIONS.
C
540 CONTINUE
DO 541 JL=1,NI
ZRFL(JL)=0.
ZSFL(JL)=0.
ZCUCOV(JL)=ZEPCOV
541 CONTINUE
C
C***
C
DO 582 JK=IKS2,NK
C***
DO 542 JL=1,NI
LO=ILAB(JL,JK).EQ.2
if (.not. LO) ZQAC(JL,JK) = 0.
temp6 = ((1.-ZBETA(JL,JK))*ZQAC(JL,JK)/ZDTTOT(JL,JK))*ZCONS3
ZCUCOV(JL) = AMIN1(AMAX1(ZCUCOV(JL),temp6),0.5)
ZCUPRO = AMAX1(ZDTLOC(JL,JK)*((1.-ZBETA(JL,JK))*ZQAC(JL,JK)
* /ZDTTOT(JL,JK)),0.)
ZRFLN(JL)= ZRFL(JL)+(ZCUPRO*ZDPP1(JL,JK)*ZCONS2)
ZSFLN(JL)= ZSFL(JL)+(ZCUPRO*ZDPP1(JL,JK)*ZCONS2)
*
*
* diagnostic d'eau liquide
* ------------------------
*
if( ilab(jl,jk) .eq. 2 ) then
*
temp6 = amax1(1.E-12,amin1(temp6,0.5))
temp4 = ((zrfln(jl)-zrfl(jl)) * (grav/zdpp1(jl,jk))) /
+ (temp6*wc0*wmr)
temp1 = temp4
*
temp5 = temp1*temp1
temp2 = temp1*(1.0-exp(-temp5)) - temp4
temp3 = 1.0 + (2.0*temp5 - 1.0)*exp(-temp5)
temp1 = amax1(temp1 - ( temp2/(temp3+1.E-12) ), 0.0)
*
temp5 = temp1*temp1
temp2 = temp1*(1.0-exp(-temp5)) - temp4
temp3 = 1.0 + (2.0*temp5 - 1.0)*exp(-temp5)
temp1 = amax1(temp1 - ( temp2/(temp3+1.E-12) ), 0.0)
*
if(ztc(jl,jk) .lt. 268.) then
temp2 = wmr / (1.0+pc2*(268.-ztc(jl,jk))**0.5)
else
temp2 = wmr
endif
*
temp2 = temp2 / (1.0 + pc1*(0.5*(zrfln(jl)+zrfl(jl)))**0.5)
*
cldw(jl,jk) = temp1 * temp6 * temp2
*
endif
*
542 CONTINUE
C
C***
IF (JK.GT.1 .AND. EVAP ) THEN
C***
DO 561 JL=1,NI
C
C* 5.5 EVAPORATION OF PRECIPITATIONS.
C WITH RECYCLING.
C
ZSQFLN = SQRT( ZRFLN(JL)/ZCUCOV(JL) )
TEMP1 = ZQSATE(JL,JK)
TEMP2 = ZTP1(JL,JK)
ZCPD(JL)=CPD*(1.+DELTA2*ZQC(JL,JK))
LO = TEMP2.LT.TRPL .AND. ZTP1(JL,NK).LT.TRPL
if (LO) then
ZLDCP0 = CHLS/ZCPD(JL)
else
ZLDCP0 = CHLC/ZCPD(JL)
endif
ZNIMP = 1. + 2.*(1.+ ZLDCP0*FODQS(TEMP1,TEMP2))
* *ZEVAP*ZSQFLN/ZCONS2
C
ZRFLN0 = ZRFLN(JL)
* ZRFLN(JL) = ZCUCOV(JL)*(AMAX1(0.,ZSQFLN-ZEVAP*ZDPP1(JL,JK)
* * *AMAX1(0.,ZQSATE(JL,JK)-ZQC(JL,JK))/ZNIMP ))**2
temp2=AMAX1(0.,ZQSATE(JL,JK)-ZQC(JL,JK))
temp3=AMAX1(0.,ZSQFLN-ZEVAP*ZDPP1(JL,JK)*temp2/ZNIMP)
ZRFLN(JL) = ZCUCOV(JL)*temp3*temp3
C
ZQP1(JL,JK) = ZQP1(JL,JK)-(1.-ZRCYCL)*
+ (ZRFLN(JL)-ZRFLN0)/(ZDPP1(JL,JK)*ZCONS2)
C
C* 5.6 MELTING/FREEZING OF PRECIPITATIONS.
C NO MORE CONSIDERED.
C
561 CONTINUE
C***
ENDIF
C***
C
C* 5.7 ADD CONVECTIVE TENDENCIES FOR T AND Q.
C
570 CONTINUE
C
DO 571 JL=1,NI
LO = ZTC(JL,JK).LT.TRPL .AND. ZTP1(JL,NK).LT.TRPL
if (LO) then
ZLVDCP = CHLS/ZCPD(JL)
else
ZLVDCP = CHLC/ZCPD(JL)
endif
TE(JL,JK)=ZLVDCP*(ZRFLN(JL)-ZRFL(JL))
+ *(GRAV/ZDPP1(JL,JK))
QE(JL,JK)=(ZQP1(JL,JK)-QP1(JL,JK))*ZCONS5
571 CONTINUE
C
C* 5.8 SWAP OF FLUXES, END OF VERTICAL LOOP
C
DO 581 JL=1,NI
ZRFL(JL)=ZRFLN(JL)
ZSFL(JL)=ZSFLN(JL)
581 CONTINUE
C***
582 CONTINUE
C
C***
C* 5.9 EFFECTS OF RECYCLING,
C* CONVECTIVE CLOUD FRACTION AND
C* CONVECTIVE RAIN RATE.
C
590 CONTINUE
C
DO 591 JK=IKS2,NK
DO 591 JL=1,NI
TE(JL,JK) = TE(JL,JK)+MAX(TE(JL,JK),0.)*ZRCYCL
+ *(ZSFL(JL)-ZRFL(JL))/MAX(1.E-12,ZSFL(JL))
LO=ILAB(JL,JK).EQ.2
if (lo) then
cck(jl,jk) = zcucov(jl)
else
cck(jl,jk) = 0.
endif
591 CONTINUE
*
C
DO 592 JL=1,NI
CRR(JL)=ZRFL(JL)+(ZSFL(JL)-ZRFL(JL))*ZRCYCL
CSR(JL)=0.0
592 CONTINUE
C
C
C
C ------------------------------------------------------------------
C
C* 6. NECESSARY COMPUTATIONS IF SUBROUTINE IS BY-PASSED.
C --------- ------------ -- ---------- -- ----------
C
600 CONTINUE
C***
C
C ------------------------------------------------------------------
C
C* 7. RETURN WORKSPACE.
C ------ ----------
C
700 CONTINUE
C
*
RETURN
END