Newer
Older
!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
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
!!** 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 =(/ &
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
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
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
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)
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
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
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
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