!MNH_LIC Copyright 2002-2021 CNRS, Meteo-France and Universite Paul Sabatier
!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
!MNH_LIC for details. version 1.
!-----------------------------------------------------------------
! ######spl
SUBROUTINE CONDENSATION( KIU, KJU, KKU, KIB, KIE, KJB, KJE, KKB, KKE, KKL, &
HFRAC_ICE, HCONDENS, HLAMBDA3, &
PPABS, PZZ, PRHODREF, PT, PRV_IN, PRV_OUT, PRC_IN, PRC_OUT, PRI_IN, PRI_OUT, &
PRS, PRG, PSIGS, PMFCONV, PCLDFR, PSIGRC, OUSERI, &
OSIGMAS, OCND2, PSIGQSAT, &
PLV, PLS, PCPH, &
PHLC_HRC, PHLC_HCF, PHLI_HRI, PHLI_HCF, &
PICE_CLD_WGT)
! ################################################################################
!
!!
!! PURPOSE
!! -------
!!** Routine to diagnose cloud fraction, liquid and ice condensate mixing ratios
!! and s'rl'/sigs^2
!!
!!
!!** METHOD
!! ------
!! Based on the large-scale fields of temperature, water vapor, and possibly
!! liquid and solid condensate, the conserved quantities r_t and h_l are constructed
!! and then fractional cloudiness, liquid and solid condensate is diagnosed.
!!
!! The total variance is parameterized as the sum of stratiform/turbulent variance
!! and a convective variance.
!! The turbulent variance is parameterized as a function of first-order moments, and
!! the convective variance is modelled as a function of the convective mass flux
!! (units kg/s m^2) as provided by the mass flux convection scheme.
!!
!! Nota: if the host model does not use prognostic values for liquid and solid condensate
!! or does not provide a convective mass flux, put all these values to zero.
!!
!!
!! EXTERNAL
!! --------
!! INI_CST
!!
!! IMPLICIT ARGUMENTS
!! ------------------
!! Module MODD_CST : contains physical constants
!!
!! REFERENCE
!! ---------
!! Chaboureau J.P. and P. Bechtold (J. Atmos. Sci. 2002)
!!
!! AUTHOR
!! ------
!! P. BECHTOLD * Laboratoire d'Aerologie *
!!
!! MODIFICATIONS
!! -------------
!! Original: 31.1.2002
!! modified : 21.3.2002
!! S.Malardel : 05.2006 : Correction sur le calcul de la fonction de
!! Bougeault F2
!! W. de Rooy: 06-06-2010: Modification in the statistical cloud scheme
!! more specifically adding a variance term
!! following ideas of Lenderink & Siebesma 2002
!! and adding a height dependence
!! S. Riette, 18 May 2010 : PSIGQSAT is added
!! S. Riette, 11 Oct 2011 : MIN function in PDF for continuity
!! modification of minimum value for Rc+Ri to create cloud and minimum value for sigma
!! Use of guess point as a starting point instead of liquid point
!! Better computation of ZCPH and dRsat/dT
!! Set ZCOND to zero if PCLDFR==0
!! Safety limitation to .99*Pressure for saturation vapour pressure
!! 2012-02 Y. Seity, add possibility to run with reversed vertical levels
!! 2014-11 K.I Ivarsson add possibility to run with OCND2 option
!! 2016 S.Riette Change INQ1
!! 2016-11 S. Riette: use HFRAC_ICE, output adjusted state
!! 2020-12 U. Andrae : Introduce SPP for HARMONIE-AROME
!! R. El Khatib 24-Aug-2021 Optimizations
!! 2021-01: SPP computations moved in aro_adjust (AROME/HARMONIE)
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
! ------------
!
USE PARKIND1, ONLY : JPRB
USE YOMHOOK , ONLY : LHOOK, DR_HOOK
USE MODD_CST
USE MODD_PARAMETERS
USE MODD_RAIN_ICE_PARAM, ONLY : XCRIAUTC, XCRIAUTI, XACRIAUTI, XBCRIAUTI
USE MODE_TIWMX, ONLY : ESATW, ESATI
USE MODE_ICECLOUD, ONLY : ICECLOUD
!
IMPLICIT NONE
!
!* 0.1 Declarations of dummy arguments :
!
!
INTEGER, INTENT(IN) :: KIU ! horizontal dimension in x
INTEGER, INTENT(IN) :: KJU ! horizontal dimension in y
INTEGER, INTENT(IN) :: KKU ! vertical dimension
INTEGER, INTENT(IN) :: KIB ! value of the first point in x
INTEGER, INTENT(IN) :: KIE ! value of the last point in x
INTEGER, INTENT(IN) :: KJB ! value of the first point in y
INTEGER, INTENT(IN) :: KJE ! value of the last point in y
INTEGER, INTENT(IN) :: KKB ! value of the first point in z
INTEGER, INTENT(IN) :: KKE ! value of the last point in z
INTEGER, INTENT(IN) :: KKL ! +1 if grid goes from ground to atmosphere top, -1 otherwise
CHARACTER(LEN=1), INTENT(IN) :: HFRAC_ICE
CHARACTER(LEN=4), INTENT(IN) :: HCONDENS
CHARACTER(LEN=*), INTENT(IN) :: HLAMBDA3 ! formulation for lambda3 coeff
REAL, DIMENSION(KIU,KJU,KKU), INTENT(IN) :: PPABS ! pressure (Pa)
REAL, DIMENSION(KIU,KJU,KKU), INTENT(IN) :: PZZ ! height of model levels (m)
REAL, DIMENSION(KIU,KJU,KKU), INTENT(IN) :: PRHODREF
REAL, DIMENSION(KIU,KJU,KKU), INTENT(INOUT) :: PT ! grid scale T (K)
REAL, DIMENSION(KIU,KJU,KKU), INTENT(IN) :: PRV_IN ! grid scale water vapor mixing ratio (kg/kg) in input
REAL, DIMENSION(KIU,KJU,KKU), INTENT(OUT) :: PRV_OUT! grid scale water vapor mixing ratio (kg/kg) in output
REAL, DIMENSION(KIU,KJU,KKU), INTENT(IN) :: PRC_IN ! grid scale r_c mixing ratio (kg/kg) in input
REAL, DIMENSION(KIU,KJU,KKU), INTENT(OUT) :: PRC_OUT! grid scale r_c mixing ratio (kg/kg) in output
REAL, DIMENSION(KIU,KJU,KKU), INTENT(IN) :: PRI_IN ! grid scale r_i (kg/kg) in input
REAL, DIMENSION(KIU,KJU,KKU), INTENT(OUT) :: PRI_OUT! grid scale r_i (kg/kg) in output
LOGICAL, INTENT(IN) :: OUSERI ! logical switch to compute both
! liquid and solid condensate (OUSERI=.TRUE.)
! or only solid condensate (OUSERI=.FALSE.)
LOGICAL, INTENT(IN) :: OSIGMAS! use present global Sigma_s values
! or that from turbulence scheme
LOGICAL, INTENT(IN) :: OCND2 ! logical switch to sparate liquid and ice
! more rigid (DEFALT value : .FALSE.)
REAL, DIMENSION(KIU,KJU), INTENT(IN) :: PSIGQSAT ! use an extra "qsat" variance contribution (OSIGMAS case)
! multiplied by PSIGQSAT
REAL, DIMENSION(KIU,KJU,KKU), INTENT(IN) :: PRS ! grid scale mixing ration of snow (kg/kg)
REAL, DIMENSION(KIU,KJU,KKU), INTENT(IN) :: PRG ! grid scale mixing ration of graupel (kg/kg)
REAL, DIMENSION(KIU,KJU,KKU), INTENT(IN) :: PSIGS ! Sigma_s from turbulence scheme
REAL, DIMENSION(:,:,:), INTENT(IN) :: PMFCONV! convective mass flux (kg /s m^2)
REAL, DIMENSION(KIU,KJU,KKU), INTENT(OUT) :: PCLDFR ! cloud fraction
REAL, DIMENSION(KIU,KJU,KKU), INTENT(OUT) :: PSIGRC ! s r_c / sig_s^2
REAL, DIMENSION(KIU,KJU,KKU), OPTIONAL, INTENT(IN) :: PLV ! Latent heat L_v
REAL, DIMENSION(KIU,KJU,KKU), OPTIONAL, INTENT(IN) :: PLS ! Latent heat L_s
REAL, DIMENSION(KIU,KJU,KKU), OPTIONAL, INTENT(IN) :: PCPH ! Specific heat C_ph
REAL, DIMENSION(KIU,KJU,KKU), OPTIONAL, INTENT(OUT) :: PHLC_HRC
REAL, DIMENSION(KIU,KJU,KKU), OPTIONAL, INTENT(OUT) :: PHLC_HCF ! cloud fraction
REAL, DIMENSION(KIU,KJU,KKU), OPTIONAL, INTENT(OUT) :: PHLI_HRI
REAL, DIMENSION(KIU,KJU,KKU), OPTIONAL, INTENT(OUT) :: PHLI_HCF
REAL, DIMENSION(KIU,KJU), OPTIONAL, INTENT(IN) :: PICE_CLD_WGT
!
!
!* 0.2 Declarations of local variables :
!
INTEGER :: JI, JJ, JK, JKP, JKM, IKTB, IKTE ! loop index
REAL, DIMENSION(KIU,KJU,KKU) :: ZTLK, ZRT ! work arrays for T_l and total water mixing ratio
REAL, DIMENSION(KIU,KJU,KKU) :: ZL ! length scale
INTEGER, DIMENSION(KIU,KJU) :: ITPL ! top levels of troposphere
REAL, DIMENSION(KIU,KJU) :: ZTMIN ! minimum Temp. related to ITPL
!
REAL, DIMENSION(KIU,KJU,KKU) :: ZLV, ZLS, ZCPD
REAL :: ZGCOND, ZAUTC, ZAUTI, ZGAUV, ZGAUC, ZGAUI, ZGAUTC, ZGAUTI, ZCRIAUTI ! Used for Gaussian PDF integration
REAL :: ZLVS ! thermodynamics
REAL, DIMENSION(KIB:KIE) :: ZPV, ZPIV, ZQSL, ZQSI ! thermodynamics
REAL :: ZLL, DZZ, ZZZ ! used for length scales
REAL :: ZAH, ZDRW, ZDTL, ZSIG_CONV ! related to computation of Sig_s
REAL, DIMENSION(KIB:KIE) :: ZA, ZB, ZSBAR, ZSIGMA, ZQ1 ! related to computation of Sig_s
REAL, DIMENSION(KIB:KIE) :: ZCOND
REAL, DIMENSION(KIB:KIE) :: ZFRAC ! Ice fraction
INTEGER :: INQ1
REAL :: ZINC
! related to OCND2 noise check :
REAL :: ZRSP, ZRSW, ZRFRAC, ZRSDIF, ZRCOLD
! related to OCND2 ice cloud calulation :
REAL, DIMENSION(KIB:KIE) :: ESATW_T
REAL :: ZDUM1,ZDUM2,ZDUM3,ZDUM4,ZPRIFACT
REAL, DIMENSION(KIU,KJU,KKU) :: TCLD
REAL :: ZDZ(KIB:KIE), ZARDUM(KIE-KIB+1),ZCLDUM(KIE-KIB+1)
! end OCND2
REAL(KIND=JPRB) :: ZHOOK_HANDLE
INTEGER, DIMENSION(KIU) :: IERR
!
!* 0.3 Definition of constants :
!
!-------------------------------------------------------------------------------
!
REAL,PARAMETER :: ZL0 = 600. ! tropospheric length scale
REAL,PARAMETER :: ZCSIGMA = 0.2 ! constant in sigma_s parameterization
REAL,PARAMETER :: ZCSIG_CONV = 0.30E-2 ! scaling factor for ZSIG_CONV as function of mass flux
!
REAL, DIMENSION(-22:11),PARAMETER :: ZSRC_1D =(/ &
0. , 0. , 2.0094444E-04, 0.316670E-03, &
4.9965648E-04, 0.785956E-03 , 1.2341294E-03, 0.193327E-02, &
3.0190963E-03, 0.470144E-02 , 7.2950651E-03, 0.112759E-01, &
1.7350994E-02, 0.265640E-01 , 4.0427860E-02, 0.610997E-01, &
9.1578111E-02, 0.135888E+00 , 0.1991484 , 0.230756E+00, &
0.2850565 , 0.375050E+00 , 0.5000000 , 0.691489E+00, &
0.8413813 , 0.933222E+00 , 0.9772662 , 0.993797E+00, &
0.9986521 , 0.999768E+00 , 0.9999684 , 0.999997E+00, &
1.0000000 , 1.000000 /)
!
!-------------------------------------------------------------------------------
!
!
IF (LHOOK) CALL DR_HOOK('CONDENSATION',0,ZHOOK_HANDLE)
IKTB=1+JPVEXT
IKTE=KKU-JPVEXT
PCLDFR(:,:,:) = 0. ! Initialize values
PSIGRC(:,:,:) = 0. ! Initialize values
ZPRIFACT = 1. ! Initialize value
ZCLDUM=-1. ! Initialize value
IF(OCND2)ZPRIFACT = 0.
!
!
!-------------------------------------------------------------------------------
! store total water mixing ratio
DO JK=IKTB,IKTE
DO JJ=KJB,KJE
DO JI=KIB,KIE
ZRT(JI,JJ,JK) = PRV_IN(JI,JJ,JK) + PRC_IN(JI,JJ,JK) + PRI_IN(JI,JJ,JK)*ZPRIFACT
END DO
END DO
END DO
!-------------------------------------------------------------------------------
! Preliminary calculations
! latent heat of vaporisation/sublimation
IF(PRESENT(PLV) .AND. PRESENT(PLS)) THEN
ZLV(:,:,:)=PLV(:,:,:)
ZLS(:,:,:)=PLS(:,:,:)
ELSE
DO JK=IKTB,IKTE
DO JJ=KJB,KJE
DO JI=KIB,KIE
! latent heat of vaporisation/sublimation
ZLV(JI,JJ,JK) = XLVTT + ( XCPV - XCL ) * ( PT(JI,JJ,JK) - XTT )
ZLS(JI,JJ,JK) = XLSTT + ( XCPV - XCI ) * ( PT(JI,JJ,JK) - XTT )
ENDDO
ENDDO
ENDDO
ENDIF
IF(PRESENT(PCPH)) THEN
ZCPD(:,:,:)=PCPH(:,:,:)
ELSE
DO JK=IKTB,IKTE
DO JJ=KJB,KJE
DO JI=KIB,KIE
ZCPD(JI,JJ,JK) = XCPD + XCPV*PRV_IN(JI,JJ,JK) + XCL*PRC_IN(JI,JJ,JK) + XCI*PRI_IN(JI,JJ,JK) + &
XCI*(PRS(JI,JJ,JK) + PRG(JI,JJ,JK) )
ENDDO
ENDDO
ENDDO
ENDIF
! Preliminary calculations needed for computing the "turbulent part" of Sigma_s
IF ( .NOT. OSIGMAS ) THEN
DO JK=IKTB,IKTE
DO JJ=KJB,KJE
DO JI=KIB,KIE
! store temperature at saturation
ZTLK(JI,JJ,JK) = PT(JI,JJ,JK) - ZLV(JI,JJ,JK)*PRC_IN(JI,JJ,JK)/ZCPD(JI,JJ,JK) &
- ZLS(JI,JJ,JK)*PRI_IN(JI,JJ,JK)/ZCPD(JI,JJ,JK)*ZPRIFACT
END DO
END DO
END DO
! Determine tropopause/inversion height from minimum temperature
ITPL(:,:) = KIB+1
ZTMIN(:,:) = 400.
DO JK = IKTB+1,IKTE-1
DO JJ=KJB,KJE
DO JI=KIB,KIE
IF ( PT(JI,JJ,JK) < ZTMIN(JI,JJ) ) THEN
ZTMIN(JI,JJ) = PT(JI,JJ,JK)
ITPL(JI,JJ) = JK
ENDIF
END DO
END DO
END DO
! Set the mixing length scale
ZL(:,:,KKB) = 20.
DO JK = KKB+KKL,KKE,KKL
DO JJ=KJB,KJE
DO JI=KIB,KIE
! free troposphere
ZL(JI,JJ,JK) = ZL0
ZZZ = PZZ(JI,JJ,JK) - PZZ(JI,JJ,KKB)
JKP = ITPL(JI,JJ)
! approximate length for boundary-layer
IF ( ZL0 > ZZZ ) ZL(JI,JJ,JK) = ZZZ
! gradual decrease of length-scale near and above tropopause
IF ( ZZZ > 0.9*(PZZ(JI,JJ,JKP)-PZZ(JI,JJ,KKB)) ) &
ZL(JI,JJ,JK) = .6 * ZL(JI,JJ,JK-KKL)
END DO
END DO
END DO
END IF
!-------------------------------------------------------------------------------
!
DO JK=IKTB,IKTE
JKP=MAX(MIN(JK+KKL,IKTE),IKTB)
JKM=MAX(MIN(JK-KKL,IKTE),IKTB)
DO JJ=KJB,KJE
IF (OCND2) THEN
ZDZ(KIB:KIE) = PZZ(KIB:KIE,JJ,JKP) - PZZ(KIB:KIE,JJ,JKP+1)
CALL ICECLOUD(KIE-KIB+1,PPABS(KIB,JJ,JK),PZZ(KIB,JJ,JK),ZDZ(KIB), &
& PT(KIB,JJ,JK),PRV_IN(KIB,JJ,JK),1.,-1., &
& ZCLDUM,1.,TCLD(KIB,JJ,JK), &
& ZARDUM,ZARDUM,ZARDUM,ZARDUM)
! latent heats
! saturated water vapor mixing ratio over liquid water and ice
DO JI=KIB,KIE
ESATW_T(JI)=ESATW(PT(JI,JJ,JK))
ZPV(JI) = MIN(ESATW_T(JI), .99*PPABS(JI,JJ,JK))
ZPIV(JI) = MIN(ESATI(PT(JI,JJ,JK)), .99*PPABS(JI,JJ,JK))
END DO
ELSE
! latent heats
! saturated water vapor mixing ratio over liquid water and ice
DO JI=KIB,KIE
ZPV(JI) = MIN(EXP( XALPW - XBETAW / PT(JI,JJ,JK) - XGAMW * LOG( PT(JI,JJ,JK) ) ), .99*PPABS(JI,JJ,JK))
ZPIV(JI) = MIN(EXP( XALPI - XBETAI / PT(JI,JJ,JK) - XGAMI * LOG( PT(JI,JJ,JK) ) ), .99*PPABS(JI,JJ,JK))
END DO
ENDIF
!Ice fraction
ZFRAC(:) = 0.
IF (OUSERI .AND. .NOT.OCND2) THEN
DO JI=KIB,KIE
IF (PRC_IN(JI,JJ,JK)+PRI_IN(JI,JJ,JK) > 1.E-20) THEN
ZFRAC(JI) = PRI_IN(JI,JJ,JK) / (PRC_IN(JI,JJ,JK)+PRI_IN(JI,JJ,JK))
ENDIF
END DO
CALL COMPUTE_FRAC_ICE(HFRAC_ICE, ZFRAC(:), PT(:,JJ,JK), IERR) !error code IERR cannot be checked here to not break vectorization
ENDIF
DO JI=KIB,KIE
ZQSL(JI) = XRD / XRV * ZPV(JI) / ( PPABS(JI,JJ,JK) - ZPV(JI) )
ZQSI(JI) = XRD / XRV * ZPIV(JI) / ( PPABS(JI,JJ,JK) - ZPIV(JI) )
! interpolate between liquid and solid as function of temperature
ZQSL(JI) = (1. - ZFRAC(JI)) * ZQSL(JI) + ZFRAC(JI) * ZQSI(JI)
ZLVS = (1. - ZFRAC(JI)) * ZLV(JI,JJ,JK) + &
& ZFRAC(JI) * ZLS(JI,JJ,JK)
! coefficients a and b
ZAH = ZLVS * ZQSL(JI) / ( XRV * PT(JI,JJ,JK)**2 ) * (XRV * ZQSL(JI) / XRD + 1.)
ZA(JI) = 1. / ( 1. + ZLVS/ZCPD(JI,JJ,JK) * ZAH )
ZB(JI) = ZAH * ZA(JI)
ZSBAR(JI) = ZA(JI) * ( ZRT(JI,JJ,JK) - ZQSL(JI) + &
& ZAH * ZLVS * (PRC_IN(JI,JJ,JK)+PRI_IN(JI,JJ,JK)*ZPRIFACT) / ZCPD(JI,JJ,JK))
END DO
! switch to take either present computed value of SIGMAS
! or that of Meso-NH turbulence scheme
IF ( OSIGMAS ) THEN
DO JI=KIB,KIE
IF (PSIGQSAT(JI,JJ)/=0.) THEN
ZSIGMA(JI) = SQRT((2*PSIGS(JI,JJ,JK))**2 + (PSIGQSAT(JI,JJ)*ZQSL(JI)*ZA(JI))**2)
ELSE
ZSIGMA(JI) = 2*PSIGS(JI,JJ,JK)
END IF
END DO
ELSE
DO JI=KIB,KIE
! parameterize Sigma_s with first_order closure
DZZ = PZZ(JI,JJ,JKP) - PZZ(JI,JJ,JKM)
ZDRW = ZRT(JI,JJ,JKP) - ZRT(JI,JJ,JKM)
ZDTL = ZTLK(JI,JJ,JKP) - ZTLK(JI,JJ,JKM) + XG/ZCPD(JI,JJ,JK) * DZZ
ZLL = ZL(JI,JJ,JK)
! standard deviation due to convection
ZSIG_CONV =0.
IF( SIZE(PMFCONV) /= 0) ZSIG_CONV = ZCSIG_CONV * PMFCONV(JI,JJ,JK) / ZA(JI)
! zsigma should be of order 4.e-4 in lowest 5 km of atmosphere
ZSIGMA(JI) = SQRT( MAX( 1.E-25, ZCSIGMA * ZCSIGMA * ZLL*ZLL/(DZZ*DZZ)*(&
ZA(JI)*ZA(JI)*ZDRW*ZDRW - 2.*ZA(JI)*ZB(JI)*ZDRW*ZDTL + ZB(JI)*ZB(JI)*ZDTL*ZDTL) + &
ZSIG_CONV * ZSIG_CONV ) )
END DO
END IF
DO JI=KIB,KIE
ZSIGMA(JI)= MAX( 1.E-10, ZSIGMA(JI) )
! normalized saturation deficit
ZQ1(JI) = ZSBAR(JI)/ZSIGMA(JI)
END DO
IF(HCONDENS == 'GAUS') THEN
DO JI=KIB,KIE
! Gaussian Probability Density Function around ZQ1
! Computation of ZG and ZGAM(=erf(ZG))
ZGCOND = -ZQ1(JI)/SQRT(2.)
!Approximation of erf function for Gaussian distribution
ZGAUV = 1 - SIGN(1., ZGCOND) * SQRT(1-EXP(-4*ZGCOND**2/XPI))
!Computation Cloud Fraction
PCLDFR(JI,JJ,JK) = MAX( 0., MIN(1.,0.5*ZGAUV))
!Computation of condensate
ZCOND(JI) = (EXP(-ZGCOND**2)-ZGCOND*SQRT(XPI)*ZGAUV)*ZSIGMA(JI)/SQRT(2.*XPI)
ZCOND(JI) = MAX(ZCOND(JI), 0.)
PSIGRC(JI,JJ,JK) = PCLDFR(JI,JJ,JK)
END DO
!Computation warm/cold Cloud Fraction and content in high water content part
IF(PRESENT(PHLC_HCF) .AND. PRESENT(PHLC_HRC))THEN
DO JI=KIB,KIE
IF(1-ZFRAC(JI) > 1.E-20)THEN
ZAUTC = (ZSBAR(JI) - XCRIAUTC/(PRHODREF(JI,JJ,JK)*(1-ZFRAC(JI))))/ZSIGMA(JI)
ZGAUTC = -ZAUTC/SQRT(2.)
!Approximation of erf function for Gaussian distribution
ZGAUC = 1 - SIGN(1., ZGAUTC) * SQRT(1-EXP(-4*ZGAUTC**2/XPI))
PHLC_HCF(JI,JJ,JK) = MAX( 0., MIN(1.,0.5*ZGAUC))
PHLC_HRC(JI,JJ,JK) = (1-ZFRAC(JI))*(EXP(-ZGAUTC**2)-ZGAUTC*SQRT(XPI)*ZGAUC)*ZSIGMA(JI)/SQRT(2.*XPI)
PHLC_HRC(JI,JJ,JK) = PHLC_HRC(JI,JJ,JK) + XCRIAUTC/PRHODREF(JI,JJ,JK) * PHLC_HCF(JI,JJ,JK)
PHLC_HRC(JI,JJ,JK) = MAX(PHLC_HRC(JI,JJ,JK), 0.)
ELSE
PHLC_HCF(JI,JJ,JK)=0.
PHLC_HRC(JI,JJ,JK)=0.
ENDIF
END DO
ENDIF
IF(PRESENT(PHLI_HCF) .AND. PRESENT(PHLI_HRI))THEN
DO JI=KIB,KIE
IF(ZFRAC(JI) > 1.E-20)THEN
ZCRIAUTI=MIN(XCRIAUTI,10**(XACRIAUTI*(PT(JI,JJ,JK)-XTT)+XBCRIAUTI))
ZAUTI = (ZSBAR(JI) - ZCRIAUTI/ZFRAC(JI))/ZSIGMA(JI)
ZGAUTI = -ZAUTI/SQRT(2.)
!Approximation of erf function for Gaussian distribution
ZGAUI = 1 - SIGN(1., ZGAUTI) * SQRT(1-EXP(-4*ZGAUTI**2/XPI))
PHLI_HCF(JI,JJ,JK) = MAX( 0., MIN(1.,0.5*ZGAUI))
PHLI_HRI(JI,JJ,JK) = ZFRAC(JI)*(EXP(-ZGAUTI**2)-ZGAUTI*SQRT(XPI)*ZGAUI)*ZSIGMA(JI)/SQRT(2.*XPI)
PHLI_HRI(JI,JJ,JK) = PHLI_HRI(JI,JJ,JK) + ZCRIAUTI*PHLI_HCF(JI,JJ,JK)
PHLI_HRI(JI,JJ,JK) = MAX(PHLI_HRI(JI,JJ,JK), 0.)
ELSE
PHLI_HCF(JI,JJ,JK)=0.
PHLI_HRI(JI,JJ,JK)=0.
ENDIF
END DO
ENDIF
ELSEIF(HCONDENS == 'CB02')THEN
DO JI=KIB,KIE
!Total condensate
IF (ZQ1(JI) > 0. .AND. ZQ1(JI) <= 2) THEN
ZCOND(JI) = MIN(EXP(-1.)+.66*ZQ1(JI)+.086*ZQ1(JI)**2, 2.) ! We use the MIN function for continuity
ELSE IF (ZQ1(JI) > 2.) THEN
ZCOND(JI) = ZQ1(JI)
ELSE
ZCOND(JI) = EXP( 1.2*ZQ1(JI)-1. )
ENDIF
ZCOND(JI) = ZCOND(JI) * ZSIGMA(JI)
!Cloud fraction
IF (ZCOND(JI) < 1.E-12) THEN
PCLDFR(JI,JJ,JK) = 0.
ELSE
PCLDFR(JI,JJ,JK) = MAX( 0., MIN(1.,0.5+0.36*ATAN(1.55*ZQ1(JI))) )
ENDIF
IF (PCLDFR(JI,JJ,JK)==0.) THEN
ZCOND(JI)=0.
ENDIF
INQ1 = MIN( MAX(-22,FLOOR(MIN(100., MAX(-100., 2*ZQ1(JI)))) ), 10) !inner min/max prevents sigfpe when 2*zq1 does not fit into an int
ZINC = 2.*ZQ1(JI) - INQ1
PSIGRC(JI,JJ,JK) = MIN(1.,(1.-ZINC)*ZSRC_1D(INQ1)+ZINC*ZSRC_1D(INQ1+1))
END DO
IF(PRESENT(PHLC_HCF) .AND. PRESENT(PHLC_HRC))THEN
PHLC_HCF(:,JJ,JK)=0.
PHLC_HRC(:,JJ,JK)=0.
ENDIF
IF(PRESENT(PHLI_HCF) .AND. PRESENT(PHLI_HRI))THEN
PHLI_HCF(:,JJ,JK)=0.
PHLI_HRI(:,JJ,JK)=0.
ENDIF
END IF !HCONDENS
IF(.NOT. OCND2) THEN
DO JI=KIB,KIE
PRC_OUT(JI,JJ,JK) = (1.-ZFRAC(JI)) * ZCOND(JI) ! liquid condensate
PRI_OUT(JI,JJ,JK) = ZFRAC(JI) * ZCOND(JI) ! solid condensate
PT(JI,JJ,JK) = PT(JI,JJ,JK) + ((PRC_OUT(JI,JJ,JK)-PRC_IN(JI,JJ,JK))*ZLV(JI,JJ,JK) + &
&(PRI_OUT(JI,JJ,JK)-PRI_IN(JI,JJ,JK))*ZLS(JI,JJ,JK) ) &
& /ZCPD(JI,JJ,JK)
PRV_OUT(JI,JJ,JK) = ZRT(JI,JJ,JK) - PRC_OUT(JI,JJ,JK) - PRI_OUT(JI,JJ,JK)*ZPRIFACT
END DO
ELSE
DO JI=KIB,KIE
PRC_OUT(JI,JJ,JK) = (1.-ZFRAC(JI)) * ZCOND(JI) ! liquid condensate
!
! This check is mainly for noise reduction :
! -------------------------
IF(ABS(PRC_IN(JI,JJ,JK)-PRC_OUT(JI,JJ,JK))>1.0E-12 .AND. ESATW_T(JI) < PPABS(JI,JJ,JK)*0.5)THEN
ZRCOLD = PRC_OUT(JI,JJ,JK)
ZRFRAC = PRV_IN(JI,JJ,JK) - ZCOND(JI) + PRC_OUT(JI,JJ,JK)
IF( PRV_IN(JI,JJ,JK) < ZRSW )THEN ! sub - saturation over water:
! Avoid drying of cloudwater leading to supersaturation with
! respect to water
ZRSDIF= MIN(0.,ZRSP-ZRFRAC)
ELSE ! super - saturation over water:
! Avoid depostition of water leading to sub-saturation with
! respect to water
! ZRSDIF= MAX(0.,ZRSP-ZRFRAC)
ZRSDIF= MAX(0.,ZRSP*PCLDFR(JI,JJ,JK) - ZRFRAC)
ENDIF
PRC_OUT(JI,JJ,JK) = ZCOND(JI) - ZRSDIF
ELSE
ZRCOLD = PRC_IN(JI,JJ,JK)
ENDIF
! end check
! compute separate ice cloud:
ZDUM1 = MIN(1.0,20.* PRC_OUT(JI,JJ,JK)*SQRT(ZDZ(JI))/ZQSL(JI)) ! clould liquid water
! factor
ZDUM3 = MAX(0.,TCLD(JI,JJ,JK)-PCLDFR(JI,JJ,JK)) ! pure ice cloud part
IF (JK==IKTB) THEN
ZDUM4 = PRI_IN(JI,JJ,JK)
ELSE
ZDUM4 = PRI_IN(JI,JJ,JK) + PRS(JI,JJ,JK)*0.5 + PRG(JI,JJ,JK)*0.25
ENDIF
ZDUM4 = MAX(0.,MIN(1.,PICE_CLD_WGT(JI,JJ)*ZDUM4*SQRT(ZDZ(JI))/ZQSI(JI))) ! clould ice+solid
! precip. water factor
ZDUM2 = (0.8*PCLDFR(JI,JJ,JK)+0.2)*MIN(1.,ZDUM1 + ZDUM4*PCLDFR(JI,JJ,JK))
! water cloud, use 'statistical' cloud, but reduce it in case of low liquid content
PCLDFR(JI,JJ,JK) = MIN(1., ZDUM2 + (0.9*ZDUM3+0.1)*ZDUM4) ! Rad cloud
! Reduce ice cloud part in case of low ice water content
PRI_OUT(JI,JJ,JK) = PRI_IN(JI,JJ,JK)
PT(JI,JJ,JK) = PT(JI,JJ,JK) + ((PRC_OUT(JI,JJ,JK)-ZRCOLD)*ZLV(JI,JJ,JK) + &
&(PRI_OUT(JI,JJ,JK)-PRI_IN(JI,JJ,JK))*ZLS(JI,JJ,JK) ) &
& /ZCPD(JI,JJ,JK)
PRV_OUT(JI,JJ,JK) = ZRT(JI,JJ,JK) - PRC_OUT(JI,JJ,JK) - PRI_OUT(JI,JJ,JK)*ZPRIFACT
END DO
END IF ! End OCND2
IF(HLAMBDA3=='CB')THEN
DO JI=KIB,KIE
! s r_c/ sig_s^2
! PSIGRC(JI,JJ,JK) = PCLDFR(JI,JJ,JK) ! use simple Gaussian relation
!
! multiply PSRCS by the lambda3 coefficient
!
! PSIGRC(JI,JJ,JK) = 2.*PCLDFR(JI,JJ,JK) * MIN( 3. , MAX(1.,1.-ZQ1(JI)) )
! in the 3D case lambda_3 = 1.
PSIGRC(JI,JJ,JK) = PSIGRC(JI,JJ,JK)* MIN( 3. , MAX(1.,1.-ZQ1(JI)) )
END DO
END IF
END DO
END DO
!
IF (LHOOK) CALL DR_HOOK('CONDENSATION',1,ZHOOK_HANDLE)
!
CONTAINS
INCLUDE "compute_frac_ice.func.h"
!
END SUBROUTINE CONDENSATION