!copyright (C) 2001 MSC-RPN COMM %%%RPNKUO%%%
***S/P LIN_KUOSYM
*
#include "phy_macros_f.h"
SUBROUTINE LIN_KUOSYM ( CTT,CQT,ilab,CCF,DBDT,CRR,CSR, 2,5
+ TP,TM,QP,QM,GZM,PSP,
+ SIGMA, TAU, NI, NK )
#include "impnone.cdk"
*
C
INTEGER NI,NK
REAL CTT(NI,NK),CQT(NI,NK)
INTEGER ilab(NI,NK)
REAL CCF(NI,NK),DBDT(NI),CRR(NI),CSR(NI)
REAL TP(NI,NK),TM(NI,NK),QP(NI,NK),QM(NI,NK),GZM(NI,NK)
REAL PSP(NI),SIGMA(NI,NK)
REAL TAU
*
*Authors
* Claude Girard and Gerard Pellerin 1995
*
*Revision
* 001 B.Bilodeau (Jan 2001) - Automatic arrays
* 002 J.-F. Mahfouf (Sept 2002) - Reecriture pour codage TL/AD versions
*
*Object
* To calculate the convective tendencies of T and Q
* using a scheme with a "symmetric Kuo-type closure".
* Geleyn's method is used to obtain the cloud profiles.
*
*Arguments
*
* - Outputs -
* CTT convective temperature tendency
* CQT convective specific humidity tendency
* ilab flag array: an indication of convective activity
* CCF estimated cumulus cloud fraction
* DBDT estimated averaged cloud fraction growth rate
* CRR rate of liquid convective precipitation
* CSR rate of solid convective precipitation
*
* - Inputs -
* TP temperature at (t+dt)
* TM temperature at (t-dt)
* QP specific humidity at (t+dt)
* QM specific humidity at (t-dt)
* GZM geopotential
* PSP surface pressure at (t+dt)
* SIGMA sigma levels
* TAU effective timestep (2*dt)
* NI horizontal dimension
* NK vertical dimension
*
*Notes
* The routine is divided into 5 parts:
* 1)allocation and position for work space
* 2)preliminary computations
* 3)cloud ascent and flagging
* 4)total dry and moist enthalpy accession calculations
* 5)cloud heating and moistening (drying) calculations
*
**
LOGICAL LO
INTEGER IS,IKA,IKB,jk,jkm1,jl
REAL ZTVC
REAL ENTRM,TAUCU,CHLS,DELTA2
REAL ZA1,ZA3,ZA4,ZSIG
REAL ZK,ZDH,ZTDIF,ZQDIF
*
C
************************************************************************
* AUTOMATIC ARRAYS
************************************************************************
*
AUTOMATIC ( LO1 , LOGICAL , (NI ))
*
AUTOMATIC ( ZCP , REAL , (NI ))
AUTOMATIC ( ZLDCP0 , REAL , (NI ))
AUTOMATIC ( ZPR , REAL , (NI ))
AUTOMATIC ( ZT1 , REAL , (NI ))
AUTOMATIC ( ZQ1 , REAL , (NI ))
AUTOMATIC ( ZP1 , REAL , (NI ))
AUTOMATIC ( ZPP , REAL , (NI,NK))
AUTOMATIC ( ZDSG , REAL , (NI,NK))
AUTOMATIC ( ZDP , REAL , (NI,NK))
AUTOMATIC ( ZSDP , REAL , (NI,NK))
AUTOMATIC ( ZQAC , REAL , (NI,NK))
AUTOMATIC ( ZTAC , REAL , (NI,NK))
AUTOMATIC ( ZSTAC , REAL , (NI,NK))
AUTOMATIC ( ZHAC , REAL , (NI,NK))
AUTOMATIC ( ZSHAC , REAL , (NI,NK))
AUTOMATIC ( ZQSE , REAL , (NI,NK))
AUTOMATIC ( ZTC , REAL , (NI,NK))
AUTOMATIC ( ZQC , REAL , (NI,NK))
AUTOMATIC ( ZTE , REAL , (NI,NK))
AUTOMATIC ( ZQE , REAL , (NI,NK))
AUTOMATIC ( ZTVE , REAL , (NI,NK))
AUTOMATIC ( ZDQ , REAL , (NI,NK))
AUTOMATIC ( ZDT , REAL , (NI,NK))
AUTOMATIC ( ZSDH , REAL , (NI,NK))
AUTOMATIC ( ZESE , REAL , (NI,NK))
*
************************************************************************
C
C* PHYSICAL CONSTANTS.
C -------- ----------
C
#include "consphy.cdk"
C
ENTRM = 5.E-6
TAUCU = 3600.
DELTA2 = CPV/CPD - 1.
CHLS = CHLC + CHLF
ZA1 = 610.78
C
C ------------------------------------------------------------------
C
C* 1. ALLOCATION AND POSITION FOR WORK SPACE.
C ---------- --- -------- --- ---- ------
C
C
C
C***
C
C METHOD.
C -------
C
C IN (3) A NEARLY ADIABATIC ASCENT IS ATTEMPTED FOR A CLOUD
C PARCEL STARTING FROM THE LOWEST MODEL LAYER. THIS CLOUD ASCENT
C IS COMPUTED IN TERMS OF TEMPERATURE AND SPECIFIC HUMIDITY.
C ENTRAINMENT IS SIMULATED VIA AN ENTRAINMENT PARAMETER.
C THE LAYERS ARE FLAGGED ACCORDING TO THE FOLLOWING CODE:
C 0 = STABLE OR INACTIVE LAYER,
C 1 = PART OF THE WELL MIXED BOUNDARY LAYER OR DRY UNSTABLE LAYER,
C 2 = MOIST UNSTABLE OR ACTIVE OR CLOUD LAYER.
C THE 1-FLAGS ARE RESET TO 0-FLAGS FOR THE NEXT SECTION.
C IN (4) THE INTEGRATED MOIST AND DRY ENTHALPY ACCESSIONS
C FOR EACH CLOUD LAYER ARE STORED INTO ALL THE CORRESPONDING
C LAYERS IF THE FIRST IS POSITIVE WHILE THE SECOND IS NEGATIVE,
C OTHERWISE, THE 2-FLAGS ARE ALSO RESET TO 0-FLAGS.
C IN (5) THE ACTUAL MODIFICATIONS OF TEMPERATURE AND SPECIFIC
C HUMIDITY ARE COMPUTED. A CLOUD-COVER VALUE IS ESTIMATED BY
C COMPARING THE TIME AT WHICH THE ENVIRONMENT WOULD REACH
C EQUILIBRIUM WITH THE CLOUD TO A PRESCRIBED CLOUD LIFE-TIME.
C
C ------------------------------------------------------------------
C
C* 2. PRELIMINARY COMPUTATIONS.
C ----------- -------------
C
C* 2.0 INITIALISATION OF OUTPUT FIELDS
C
ilab(:,:) = 0
CTT(:,:) = 0.0
CQT(:,:) = 0.0
CCF(:,:) = 0.0
CRR(:) = 0.0
CSR(:) = 0.0
C
C* 2.1 ENVIRONMENTAL PROFILES AND PARAMETERS,
C* DRY AND MOIST ENTHALPY ACCESSIONS (divided by cp)
C* AND INITIALIZATIONS.
C
DO jl=1,NI
ZDSG(jl,1)=0.5*(SIGMA(jl,2)-SIGMA(jl,1))
END DO
C
DO jk=2,NK-1
DO jl=1,NI
ZDSG(jl,jk)=0.5*(SIGMA(jl,jk+1)-SIGMA(jl,jk-1))
END DO
END DO
C
DO jl=1,NI
ZDSG(jl,NK)=0.5*(1.-SIGMA(jl,NK-1))+0.5*(1.-SIGMA(jl,NK))
END DO
C
DO jl=1,NI
DBDT(jl) = 0.
IF (TP(jl,NK) < TRPL) THEN
ZLDCP0(jl) = CHLS / CPD
ELSE
ZLDCP0(jl) = CHLC / CPD
ENDIF
END DO
C
DO jk=1,NK
DO jl=1,NI
ZPP(jl,jk)=SIGMA(jl,jk)*PSP(jl)
ZDP(jl,jk)=ZDSG(jl,jk)*PSP(jl)
ZTE(jl,jk)=TP(jl,jk)
IF (ZTE(jl,jk) > TRPL) THEN
ZA3=17.269
ZA4=35.860
ELSE
ZA3=21.875
ZA4= 7.660
ENDIF
ZESE(jl,jk)=ZA1*EXP(ZA3*(ZTE(jl,jk)-TRPL)/(ZTE(jl,jk)-ZA4))
ZQSE(jl,jk)=EPS1*ZESE(jl,jk)/(ZPP(jl,jk)-EPS2*ZESE(jl,jk))
IF (ZQSE(jl,jk) > QM (jl,jk)) THEN
ZQE(jl,jk) = QM(jl,jk)
ELSE
ZQE(jl,jk) = ZQSE(jl,jk)
ENDIF
ZTVE(jl,jk) = ZTE(jl,jk)*(1.0 + DELTA*ZQE(jl,jk))
C
ZTAC(jl,jk)=(TP(jl,jk)-TM(jl,jk))*ZDP(jl,jk)/TAU
ZQAC(jl,jk)=(QP(jl,jk)-ZQE(jl,jk))*ZDP(jl,jk)/TAU
C
ZHAC(jl,jk)= ZTAC(jl,jk) + ZLDCP0(jl)*ZQAC(jl,jk)
END DO
END DO
C
C* 2.2 SPECIFY TC AND QC AT THE LOWEST LAYER TO START THE
C* CLOUD ASCENT. CHECK FOR POSITIVE ACCESSION
C* BETWEEN SURFACE AND CLOUD BASE.
C* ZQC=0 INDICATES STABLE CONDITIONS.
C
DO jl=1,NI
ZPR(jl) = 0.
ZTC(jl,NK)=ZTE(jl,NK)
ZQC(jl,NK)=0.
IF (ZHAC(jl,NK) > 0.) THEN
ZQC(jl,NK)=ZQE(jl,NK)
ilab(jl,NK) = 1
ENDIF
END DO
C
C ------------------------------------------------------------------
C
C* 3. CLOUD ASCENT AND FLAGGING.
C ----- ------ --- ---------
C
C* 3.1 CALCULATE TC AND QC AT UPPER LEVELS BY DRY ADIABATIC
C* LIFTING FOLLOWED BY LATENT HEAT RELEASE WHEN REQUIRED.
C* CONDENSATION CALCULATIONS ARE DONE WITH TWO ITERATIONS.
C***
DO jk=NK-1,1,-1
C***
DO jl=1,NI
ZCP(jl)=CPD*(1.+DELTA2*ZQC(jl,jk+1))
IF ( ZTC(jl,jk+1) > ZTE(jl,jk+1) ) THEN
ZTDIF = ZTC(jl,jk+1)-ZTE(jl,jk+1)
ELSE
ZTDIF = 0.0
ENDIF
ZTC(jl,jk)=ZTC(jl,jk+1)+(GZM(jl,jk+1)-GZM(jl,jk))*
* (1./ZCP(jl)+ENTRM*ZTDIF)
IF ( ZQC(jl,jk+1) > ZQE(jl,jk+1) ) THEN
ZQDIF = ZQC(jl,jk+1)-ZQE(jl,jk+1)
ELSE
ZQDIF = 0.0
ENDIF
ZQC(jl,jk)=ZQC(jl,jk+1)+(GZM(jl,jk+1)-GZM(jl,jk))*
* ENTRM*ZQDIF
ZTVC = ZTC(jl,jk)*(1.0 + DELTA*ZQC(jl,jk))
LO= ZTVC.GT.ZTVE(jl,jk) .AND. ZQC(jl,jk).NE.0.
IF (LO) ilab(jl,jk) = 1
END DO
C
DO jl=1,NI
ZP1(jl) = ZPP(jl,jk)
ZT1(jl) = ZTC(jl,jk)
ZQ1(jl) = ZQC(jl,jk)
END DO
C
CALL LIN_ADJTQ ( ZT1, ZQ1, ZP1, NI )
C
DO jl=1,NI
LO1(jl) = ZT1(jl) /= ZTC(jl,jk)
ZTC(jl,jk) = ZT1(jl)
ZQC(jl,jk) = ZQ1(jl)
END DO
C
DO jl=1,NI
ZTVC=ZTC(jl,jk)*(1.0+DELTA*ZQC(jl,jk))
LO= ZTVC.GT.ZTVE(jl,jk) .AND. LO1(jl)
IF (LO) ilab(jl,jk) = 2
LO1(jl)=ilab(jl,jk).EQ.0
IF (LO1(jl)) THEN
ZTC(jl,jk) = ZTE(jl,jk)
ZQC(jl,jk) = 0.0
ENDIF
END DO
C
C* 3.2 IF NOT AT THE TOP CHECK FOR NEW LIFTING LEVEL, I.E.
C* ENTHALPY ACCESSION IN A STABLE LAYER.
C***
IF (jk.NE.1) THEN
DO jl=1,NI
IF (LO1(jl).AND.(ZHAC(jl,jk) > 0.)) THEN
ZTC(jl,jk) = ZTE(jl,jk)
ZQC(jl,jk) = ZQE(jl,jk)
ENDIF
END DO
ENDIF
C***
END DO
C***
C* 3.3 ilab=0 UNLESS ilab=2
C* IKA INDICATES THE HIGHEST TOP OF A CLOUD
C* (TO AVOID UNNECESSARY COMPUTATIONS LATER).
C
IKA=NK+1
C
DO jk=1,NK
C
DO jl=1,NI
IF (ilab(jl,jk) == 1) ilab(jl,jk) = 0
END DO
C
IF (IKA == NK+1) THEN
IS=0
DO jl=1,NI
IS=IS+ilab(jl,jk)
END DO
IF (IS.NE.0) IKA=jk
ENDIF
C
END DO
C***
IF (IKA == NK+1) RETURN
C***
C ------------------------------------------------------------------
C
C* 4. TOTAL ENERGY ACCESSION
C ----- ------ ---------
C
C* 4.1 CALCULATE TOTAL ENTHALPY ACCESSIONS REQUIRING THAT
C* - TOTAL MOIST ENTHALPY ACCESSION BE > 0
C* - TOTAL DRY ENTHALPY ACCESSION BE < 0
C* IKB IS AN UPDATE OF IKA.
C
DO jl=1,NI
ZSHAC(jl,:) = 0.0
ZSTAC(jl,:) = 0.0
ZSDP(jl,:) = 0.0
END DO
C
DO jk=NK-1,IKA,-1
DO jl=1,NI
IF (ilab(jl,jk) == 2) THEN
ZSHAC(jl,jk) = ZSHAC(jl,jk+1)+ZHAC(jl,jk)
ZSTAC(jl,jk) = ZSTAC(jl,jk+1)+ZTAC(jl,jk)
ZSDP(jl,jk) = ZSDP(jl,jk+1)+ZDP(jl,jk)
ENDIF
END DO
END DO
C
IKB=NK+1
C
DO jk=IKA,NK-1
jkm1=max0(jk-1,1)
C
DO jl=1,NI
IF ((ilab(jl,jk) == 2).AND.(ilab(jl,jkm1) == 2)) THEN
ZSHAC(jl,jk) = ZSHAC(jl,jkm1)
ZSTAC(jl,jk) = ZSTAC(jl,jkm1)
ZSDP(jl,jk) = ZSDP(jl,jkm1)
ENDIF
LO = ZSHAC(jl,jk).GT.0. .and. ZSTAC(jl,jk).LT.0.
& and. ZSDP(jl,jk).GT.0.
IF (.not.LO) ilab(jl,jk) = 0
END DO
C
IF (IKB == NK+1) THEN
IS=0
DO jl=1,NI
IS=IS+ilab(jl,jk)
END DO
IF (IS.NE.0) IKB=jk
ENDIF
C
END DO
C***
IF (IKB == NK+1) RETURN
C***
C ------------------------------------------------------------------
C
C* 5. HEATING AND MOISTENING
C ----------------------
C
C* 5.1 COMPUTE THE TOTAL CLOUD-ENVIRONMENT ENTHALPY
C* DIFFERENCE IN CLOUD LAYERS.
C
DO jl=1,NI
ZSDH(jl,:)=0.
END DO
C
DO jk=NK-1,IKB,-1
DO jl=1,NI
ZTVC = ZTC(jl,jk)*(1.0 + DELTA*ZQC(jl,jk))
ZDQ(jl,jk) = (ZQSE(jl,jk)-ZQE(jl,jk))*ZDP(jl,jk)
ZDT(jl,jk) = (ZTVC-ZTVE(jl,jk))*ZDP(jl,jk)
ZDH = ZDT(jl,jk)+ZLDCP0(jl)*ZDQ(jl,jk)
IF (ilab(jl,jk) == 2) THEN
ZSDH(jl,jk) = ZSDH(jl,jk+1) + ZDH
ENDIF
END DO
END DO
C
DO jk=IKB+1,NK-1
DO jl=1,NI
IF ((ilab(jl,jk) == 2).AND.(ilab(jl,jk-1) == 2)) THEN
ZSDH(jl,jk)= ZSDH(jl,jk-1)
ENDIF
END DO
END DO
C
C* 5.2 COMPUTE CONVECTIVE HEATING AND MOISTENING.
C* ESTIMATE CONVECTIVE CLOUD FRACTION.
C
DO jk=IKB,NK-1
DO jl=1,NI
IF (ilab(jl,jk) == 0) THEN
ZQAC(jl,jk) = 0.0
ZTAC(jl,jk) = 0.0
ZSHAC(jl,jk) = 0.0
ENDIF
IF (ZSDH(jl,jk) <= 0.) THEN
ZSDH(jl,jk) = -1.0
ENDIF
C
ZK = ZSHAC(jl,jk)/ZSDH(jl,jk)
C
CQT(jl,jk) = (ZK*ZDQ(jl,jk)-ZQAC(jl,jk))/ZDP(jl,jk)
CTT(jl,jk) = (ZK*ZDT(jl,jk)-ZTAC(jl,jk))/ZDP(jl,jk)
C
ZPR(jl) = ZPR(jl) + CTT(jl,jk)/ZLDCP0(jl)*ZDP(jl,jk)
C
IF (DBDT(jl) < ZK) THEN
DBDT(jl) = ZK
ENDIF
C
END DO
END DO
C
C* STORE SURFACE PRECIPITATION
C
DO jl=1,NI
IF (ZTE(jl,NK) > TRPL) THEN
CRR(jl) = ZPR(jl) / GRAV
ELSE
CSR(jl) = ZPR(jl) / GRAV
ENDIF
END DO
C
DO jl=1,NI
ZPR(jl) = ZPR(jl) / ( GRAV * 1.E3 )
IF (ZPR(jl) < 1.0E-12) THEN
ZPR(jl) = 1.0E-12
ENDIF
ZPR(jl) = 2.5 + 0.125 * ALOG(ZPR(jl))
IF (ZPR(jl) < DBDT(jl) * TAU) THEN
ZPR(jl) = DBDT(jl) * TAU
ENDIF
IF (ZPR(jl) > 0.8) THEN
ZPR(jl) = 0.8
ENDIF
END DO
C
DO jk=IKB,NK-1
DO jl=1,NI
IF (ilab(jl,jk) == 2) THEN
CCF(jl,jk) = ZPR(jl)
ELSE
CCF(jl,jk) = 0.0
ENDIF
ZSIG=(SIGMA(jl,jk)/0.8)**2
IF (ZSIG > 1.0) ZSIG=1.0
CCF(jl,jk) = CCF(jl,jk) * ZSIG
END DO
END DO
C
RETURN
END SUBROUTINE LIN_KUOSYM