condensation.f90

!MNH_LIC Copyright 2002-2023 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
      MODULE MODI_CONDENSATION
!     ########################
!
INTERFACE
!
       SUBROUTINE CONDENSATION( KIU, KJU, KKU, KIB, KIE, KJB, KJE, KKB, KKE, KKL,&
          HFRAC_ICE, HCONDENS, HLAMBDA3, &
          PPABS, PZZ, PRHODREF, PT, PRV, PRC, PRI, PRS, PRG, PSIGS, PMFCONV, PCLDFR, PSIGRC, OUSERI,&
          OSIGMAS, PSIGQSAT, PLV, PLS, PCPH, PHLC_HRC, PHLC_HCF, PHLI_HRI, PHLI_HCF)
!
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(INOUT) :: PRV    ! grid scale water vapor mixing ratio (kg/kg)
REAL, DIMENSION(KIU,KJU,KKU), INTENT(INOUT) :: PRC    ! grid scale r_c mixing ratio (kg/kg)
REAL, DIMENSION(KIU,KJU,KKU), INTENT(INOUT) :: PRI    ! grid scale r_i (kg/kg)
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
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
REAL, INTENT(IN)                            :: PSIGQSAT ! use an extra "qsat" variance contribution (OSIGMAS case)
REAL, DIMENSION(KIU,KJU,KKU), OPTIONAL, INTENT(IN)    :: PLV
REAL, DIMENSION(KIU,KJU,KKU), OPTIONAL, INTENT(IN)    :: PLS
REAL, DIMENSION(KIU,KJU,KKU), OPTIONAL, INTENT(IN)    :: PCPH
REAL, DIMENSION(KIU,KJU,KKU), OPTIONAL, INTENT(OUT)   :: PHLC_HRC   !cloud water content in precipitating part
REAL, DIMENSION(KIU,KJU,KKU), OPTIONAL, INTENT(OUT)   :: PHLC_HCF   !precipitating part
REAL, DIMENSION(KIU,KJU,KKU), OPTIONAL, INTENT(OUT)   :: PHLI_HRI   !
REAL, DIMENSION(KIU,KJU,KKU), OPTIONAL, INTENT(OUT)   :: PHLI_HCF   !

END SUBROUTINE CONDENSATION
!
END INTERFACE
!
END MODULE MODI_CONDENSATION
!     ######spl
    SUBROUTINE CONDENSATION( KIU, KJU, KKU, KIB, KIE, KJB, KJE, KKB, KKE, KKL,         &
       HFRAC_ICE, HCONDENS, HLAMBDA3,                                                  &
       PPABS, PZZ, PRHODREF, PT, PRV, PRC, PRI, PRS, PRG, PSIGS, PMFCONV, PCLDFR, PSIGRC, OUSERI,&
       OSIGMAS, PSIGQSAT, PLV, PLS, PCPH, PHLC_HRC, PHLC_HCF, PHLI_HRI, PHLI_HCF )
!   ################################################################################
!
!!
!!    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
!!      2015   C.Lac   Change min value of ZSIGMA to be in agreement with AROME
!!      2016   G.Delautier   Restore min value of ZSIGMA (instability)
!!      2016   S.Riette Change INQ1 
!!      2016-11 S. Riette: use HFRAC_ICE, output adjusted state
!-------------------------------------------------------------------------------
!
!*       0.    DECLARATIONS
!              ------------
!
USE MODD_CST
USE MODD_PARAMETERS
USE MODD_RAIN_ICE_PARAM, ONLY : XCRIAUTC, XCRIAUTI, XACRIAUTI, XBCRIAUTI
!
#ifdef MNH_OPENACC
USE MODE_MNH_ZWORK,      ONLY: MNH_MEM_GET, MNH_MEM_POSITION_PIN, MNH_MEM_RELEASE
#endif
USE MODE_MPPDB
use mode_msg
!
USE MODI_COMPUTE_FRAC_ICE
#if defined(MNH_BITREP) || defined(MNH_BITREP_OMP)
USE MODI_BITREP
#endif
!
#if defined(MNH_COMPILER_CCE) && defined(MNH_BITREP_OMP)
!$mnh_undef(LOOP)
!$mnh_undef(OPENACC)
#endif
!
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(INOUT) :: PRV    ! grid scale water vapor mixing ratio (kg/kg)
REAL, DIMENSION(KIU,KJU,KKU), INTENT(INOUT) :: PRC    ! grid scale r_c mixing ratio (kg/kg)
REAL, DIMENSION(KIU,KJU,KKU), INTENT(INOUT) :: PRI    ! grid scale r_i (kg/kg)
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
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
REAL, INTENT(IN)                            :: PSIGQSAT ! use an extra "qsat" variance contribution (OSIGMAS case)

!
!
!*       0.2   Declarations of local variables :
!
INTEGER  :: JI, JJ, JK, JKP, JKM, IKTB, IKTE    ! loop index
INTEGER,  DIMENSION(:,:), POINTER, CONTIGUOUS  :: JKPP
REAL,    DIMENSION(:,:,:), POINTER, CONTIGUOUS :: ZTLK, ZRT       ! work arrays for T_l and total water mixing ratio
REAL,    DIMENSION(:,:,:), POINTER, CONTIGUOUS :: ZL              ! length scale
REAL,    DIMENSION(:,:,:), POINTER, CONTIGUOUS :: ZFRAC           ! Ice fraction
REAL                                           :: ZCRIAUTI        !
REAL,    DIMENSION(:,:,:), POINTER, CONTIGUOUS :: ZCRIAUTIP
INTEGER, DIMENSION(:,:),   POINTER, CONTIGUOUS :: ITPL            ! top levels of troposphere
REAL,    DIMENSION(:,:),   POINTER, CONTIGUOUS :: ZTMIN           ! minimum Temp. related to ITPL
!
REAL,    DIMENSION(:,:,:), POINTER, CONTIGUOUS :: ZLV, ZLS, ZCPD, ZCONDP
REAL                                           :: ZCOND
REAL :: ZGCOND, ZSBAR, ZSBARC, ZQ1, ZAUTC, ZAUTI, ZGAUV, ZGAUC, ZGAUI, ZGAUTC, ZGAUTI   ! Used for integration in Gaussian Probability Density Function
REAL,    DIMENSION(:,:,:), POINTER, CONTIGUOUS :: ZSBARP,ZQ1P,ZGCONDP,ZGAUVP,ZAUTCP,ZGAUTCP,ZGAUCP,ZAUTIP,ZGAUTIP,ZGAUIP
REAL :: ZTEMP, ZPV, ZQSL, ZPIV, ZQSI, ZLVS ! thermodynamics
REAL,    DIMENSION(:,:,:), POINTER, CONTIGUOUS :: ZPVP,ZQSLP,ZPIVP,ZQSIP,ZLVSP
REAL :: ZLL, DZZ
REAL,    DIMENSION(:,:),   POINTER, CONTIGUOUS :: ZZZP
REAL,    DIMENSION(:,:,:), POINTER, CONTIGUOUS :: ZDZZP,ZLLP                           ! used for length scales
REAL :: ZAH, ZA, ZB, ZSIGMA, ZDRW, ZDTL, ZSIG_CONV ! related to computation of Sig_s
REAL,    DIMENSION(:,:,:), POINTER, CONTIGUOUS :: ZAHP,ZAP,ZBP,ZSIGMAP,ZDRWP,ZDTLP,ZSIG_CONVP
REAL :: ZRCOLD, ZRIOLD
REAL,    DIMENSION(:,:,:), POINTER, CONTIGUOUS :: ZRCOLDP,ZRIOLDP
INTEGER  :: INQ1
INTEGER,    DIMENSION(:,:,:), POINTER, CONTIGUOUS :: INQ1P
REAL :: ZINC
REAL,    DIMENSION(:,:,:), POINTER, CONTIGUOUS :: ZINCP
LOGICAL :: GPRESENT_PLV, GPRESENT_PLS, GPRESENT_PCPH
LOGICAL, DIMENSION(:,:,:), POINTER, CONTIGUOUS :: GWORK
CHARACTER(LEN=4) :: YLAMBDA3 !Necessary to workaround NVHPC bug (version 21.7 if OpenACC enabled)
LOGICAL :: GPRESENT_PHLC_HCF, GPRESENT_PHLC_HRC, GPRESENT_PHLI_HCF, GPRESENT_PHLI_HRI
LOGICAL :: GYLAMBDA3_CB
!
INTEGER,    DIMENSION(:), POINTER, CONTIGUOUS :: JKPK,JKMK
!
!*       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     /)
!
!-------------------------------------------------------------------------------
!
YLAMBDA3 = HLAMBDA3

IF( YLAMBDA3 /='CB' .AND. YLAMBDA3 /='NONE' ) &
  call Print_msg( NVERB_FATAL, 'GEN', 'CONDENSATION', 'invalid value for YLAMBDA3: ' // TRIM( YLAMBDA3 ) )

GPRESENT_PLV  = PRESENT(PLV)
GPRESENT_PLS  = PRESENT(PLS)
GPRESENT_PCPH = PRESENT(PCPH)
GPRESENT_PHLC_HCF = PRESENT(PHLC_HCF)
GPRESENT_PHLC_HRC = PRESENT(PHLC_HRC)
GPRESENT_PHLI_HCF = PRESENT(PHLI_HCF)
GPRESENT_PHLI_HRI = PRESENT(PHLI_HRI)
GYLAMBDA3_CB = (YLAMBDA3=='CB')
!
IF (MPPDB_INITIALIZED) THEN
  !Check all IN arrays
  CALL MPPDB_CHECK3D(PPABS,"CONDENSATION beg:PPABS",PRECISION)
  CALL MPPDB_CHECK3D(PZZ,"CONDENSATION beg:PZZ",PRECISION)
  CALL MPPDB_CHECK3D(PRS,"CONDENSATION beg:PRS",PRECISION)
  CALL MPPDB_CHECK3D(PRG,"CONDENSATION beg:PRG",PRECISION)
  CALL MPPDB_CHECK3D(PSIGS,"CONDENSATION beg:PSIGS",PRECISION)
  CALL MPPDB_CHECK3D(PMFCONV,"CONDENSATION beg:PMFCONV",PRECISION)
  IF (GPRESENT_PLV) CALL MPPDB_CHECK3D(PLV,"CONDENSATION beg:PLV",PRECISION)
  IF (GPRESENT_PLS) CALL MPPDB_CHECK3D(PLS,"CONDENSATION beg:PLS",PRECISION)
  IF (GPRESENT_PCPH) CALL MPPDB_CHECK3D(PCPH,"CONDENSATION beg:PCPH",PRECISION)
  !Check all INOUT arrays
  CALL MPPDB_CHECK3D(PT,"CONDENSATION beg:PT",PRECISION)
  CALL MPPDB_CHECK3D(PRV,"CONDENSATION beg:PRV",PRECISION)
  CALL MPPDB_CHECK3D(PRC,"CONDENSATION beg:PRC",PRECISION)
  CALL MPPDB_CHECK3D(PRI,"CONDENSATION beg:PRI",PRECISION)
END IF

#ifndef MNH_OPENACC
allocate( ztlk (kiu, kju, kku ) )
allocate( zrt  (kiu, kju, kku ) )
allocate( zl   (kiu, kju, kku ) )
allocate( zfrac(kiu, kju, kku ) )

allocate( itpl (kiu, kju ) )
allocate( ztmin(kiu, kju ) )

allocate( zlv  (kiu, kju, kku ) )
allocate( zls  (kiu, kju, kku ) )
allocate( zcpd (kiu, kju, kku ) )
allocate( gwork(kiu, kju, kku ) )

allocate( ZZZP(kiu, kju ) )
allocate( JKPP(kiu, kju ) )
allocate( ZPVP(kiu, kju, kku ) )
allocate( ZQSLP(kiu, kju, kku ) )
allocate( ZPIVP(kiu, kju, kku ) )
allocate( ZQSIP(kiu, kju, kku ) )
allocate( ZLVSP(kiu, kju, kku ) )
allocate( ZAHP(kiu, kju, kku ) )
allocate( ZAP(kiu, kju, kku ) )
allocate( ZBP(kiu, kju, kku ) )
allocate( ZSBARP(kiu, kju, kku ) )
allocate( ZSIGMAP(kiu, kju, kku ) )
allocate( ZDZZP(kiu, kju, kku ) )
allocate( ZDRWP(kiu, kju, kku ) )
allocate( ZDTLP(kiu, kju, kku ) )
allocate( ZLLP(kiu, kju, kku ) )
allocate( ZSIG_CONVP(kiu, kju, kku ) )
allocate( ZQ1P(kiu, kju, kku ) )
allocate( ZGCONDP(kiu, kju, kku ) )
allocate( ZGAUVP(kiu, kju, kku ) )
allocate( ZCONDP(kiu, kju, kku ) )
allocate( ZAUTCP(kiu, kju, kku ) )
allocate( ZGAUTCP(kiu, kju, kku ) )
allocate( ZGAUCP(kiu, kju, kku ) )
allocate( ZCRIAUTIP(kiu, kju, kku ) )
allocate( ZAUTIP(kiu, kju, kku ) )
allocate( ZGAUTIP(kiu, kju, kku ) )
allocate( ZGAUIP(kiu, kju, kku ) )
allocate( INQ1P(kiu, kju, kku ) )
allocate( ZINCP(kiu, kju, kku ) )
allocate( ZRCOLDP(kiu, kju, kku ) )
allocate( ZRIOLDP(kiu, kju, kku ) )
allocate( JKPK(kku) )
allocate( JKMK(kku) )

#else
!Pin positions in the pools of MNH memory
CALL MNH_MEM_POSITION_PIN( 'CONDENSATION' )

CALL MNH_MEM_GET( ztlk,  kiu, kju, kku )
CALL MNH_MEM_GET( zrt,   kiu, kju, kku )
CALL MNH_MEM_GET( zl,    kiu, kju, kku )
CALL MNH_MEM_GET( zfrac, kiu, kju, kku )

CALL MNH_MEM_GET( itpl,  kiu, kju )
CALL MNH_MEM_GET( ztmin, kiu, kju )

CALL MNH_MEM_GET( zlv,   kiu, kju, kku )
CALL MNH_MEM_GET( zls,   kiu, kju, kku )
CALL MNH_MEM_GET( zcpd,  kiu, kju, kku )
CALL MNH_MEM_GET( gwork, kiu, kju, kku )

CALL MNH_MEM_GET( ZZZP, kiu, kju )
CALL MNH_MEM_GET( JKPP, kiu, kju )
CALL MNH_MEM_GET( ZPVP, kiu, kju, kku )
CALL MNH_MEM_GET( ZQSLP, kiu, kju, kku )
CALL MNH_MEM_GET( ZPIVP, kiu, kju, kku )
CALL MNH_MEM_GET( ZQSIP, kiu, kju, kku )
CALL MNH_MEM_GET( ZLVSP, kiu, kju, kku )
CALL MNH_MEM_GET( ZAHP, kiu, kju, kku )
CALL MNH_MEM_GET( ZAP, kiu, kju, kku )
CALL MNH_MEM_GET( ZBP, kiu, kju, kku )
CALL MNH_MEM_GET( ZSBARP, kiu, kju, kku )
CALL MNH_MEM_GET( ZSIGMAP, kiu, kju, kku )
CALL MNH_MEM_GET( ZDZZP, kiu, kju, kku )
CALL MNH_MEM_GET( ZDRWP, kiu, kju, kku )
CALL MNH_MEM_GET( ZDTLP, kiu, kju, kku )
CALL MNH_MEM_GET( ZLLP, kiu, kju, kku )
CALL MNH_MEM_GET( ZSIG_CONVP, kiu, kju, kku )
CALL MNH_MEM_GET( ZQ1P, kiu, kju, kku )
CALL MNH_MEM_GET( ZGCONDP, kiu, kju, kku )
CALL MNH_MEM_GET( ZGAUVP, kiu, kju, kku )
CALL MNH_MEM_GET( ZCONDP, kiu, kju, kku )
CALL MNH_MEM_GET( ZAUTCP, kiu, kju, kku )
CALL MNH_MEM_GET( ZGAUTCP, kiu, kju, kku )
CALL MNH_MEM_GET( ZGAUCP, kiu, kju, kku )
CALL MNH_MEM_GET( ZCRIAUTIP, kiu, kju, kku )
CALL MNH_MEM_GET( ZAUTIP, kiu, kju, kku )
CALL MNH_MEM_GET( ZGAUTIP, kiu, kju, kku )
CALL MNH_MEM_GET( ZGAUIP, kiu, kju, kku )
CALL MNH_MEM_GET( INQ1P, kiu, kju, kku )
CALL MNH_MEM_GET( ZINCP, kiu, kju, kku )
CALL MNH_MEM_GET( ZRCOLDP, kiu, kju, kku )
CALL MNH_MEM_GET( ZRIOLDP, kiu, kju, kku )
CALL MNH_MEM_GET( JKPK, kku )
CALL MNH_MEM_GET( JKMK, kku )

!$acc data present( PPABS, PZZ, PT, PRV, PRC, PRI, PRS, PRG, PSIGS, PMFCONV, PCLDFR, PSIGRC, &
!$acc &             ztlk, zrt, zl, zfrac, itpl, ztmin, zlv, zls,zcpd, gwork,&
!$acc &             ZZZP,JKPP,ZPVP,ZQSLP,ZPIVP,ZQSIP,ZLVSP,ZAHP,ZAP,ZBP,ZSBARP,ZSIGMAP,&
!$acc &             ZDZZP,ZDRWP,ZDTLP,ZLLP,ZSIG_CONVP,ZQ1P,ZGCONDP,ZGAUVP,ZCONDP,ZAUTCP,&
!$acc &             ZGAUTCP,ZGAUCP,ZCRIAUTIP,ZAUTIP,ZGAUTIP,ZGAUIP,INQ1P,ZINCP,ZRCOLDP,ZRIOLDP,&
!$acc &             JKPK,JKMK )
#endif

IKTB=1+JPVEXT
IKTE=KKU-JPVEXT
!
!
!$acc kernels present_cr(ZRT)
PCLDFR(:,:,:) = 0. ! Initialize values
PSIGRC(:,:,:) = 0. ! Initialize values
!
!-------------------------------------------------------------------------------
! store total water mixing ratio
  ZRT(KIB:KIE,KJB:KJE,IKTB:IKTE)  = PRV(KIB:KIE,KJB:KJE,IKTB:IKTE) + PRC(KIB:KIE,KJB:KJE,IKTB:IKTE) + PRI(KIB:KIE,KJB:KJE,IKTB:IKTE)
!$acc end kernels
!-------------------------------------------------------------------------------
! Preliminary calculations
! latent heat of vaporisation/sublimation
IF(GPRESENT_PLV .AND. GPRESENT_PLS) THEN
!$acc data present( PLV, PLS )
!$acc kernels present_cr(ZLV,ZLS)
  ZLV(:,:,:)=PLV(:,:,:)
  ZLS(:,:,:)=PLS(:,:,:)
!$acc end kernels
!$acc end data
ELSE
!$acc kernels present_cr(ZLV, ZLS)
  ! latent heat of vaporisation/sublimation
  ZLV(KIB:KIE,KJB:KJE,IKTB:IKTE) = XLVTT + ( XCPV - XCL ) * ( PT(KIB:KIE,KJB:KJE,IKTB:IKTE) - XTT )
  ZLS(KIB:KIE,KJB:KJE,IKTB:IKTE) = XLSTT + ( XCPV - XCI ) * ( PT(KIB:KIE,KJB:KJE,IKTB:IKTE) - XTT )
!$acc end kernels
ENDIF
IF(GPRESENT_PCPH) THEN
!$acc data present( PCPH )
!$acc kernels present_cr(ZCPD)
  ZCPD(:,:,:)=PCPH(:,:,:)
!$acc end kernels
!$acc end data
ELSE
!$acc kernels present_cr(ZCPD)
     ZCPD(KIB:KIE,KJB:KJE,IKTB:IKTE) = XCPD + XCPV*PRV(KIB:KIE,KJB:KJE,IKTB:IKTE) &
              + XCL*PRC(KIB:KIE,KJB:KJE,IKTB:IKTE) + XCI*PRI(KIB:KIE,KJB:KJE,IKTB:IKTE) &
              + XCI*(PRS(KIB:KIE,KJB:KJE,IKTB:IKTE) + PRG(KIB:KIE,KJB:KJE,IKTB:IKTE) )
!$acc end kernels
ENDIF
!-------------------------------------------------------------------------------
!acc kernels
! Preliminary calculations needed for computing the "turbulent part" of Sigma_s
IF ( .NOT. OSIGMAS ) THEN
!$acc kernels present_cr(ZTLK,ITPL,ZTMIN,ZZZP)
  ! store temperature at saturation
  ZTLK(KIB:KIE,KJB:KJE,IKTB:IKTE) =  PT(KIB:KIE,KJB:KJE,IKTB:IKTE) &
              - ZLV(KIB:KIE,KJB:KJE,IKTB:IKTE)*PRC(KIB:KIE,KJB:KJE,IKTB:IKTE)/ZCPD(KIB:KIE,KJB:KJE,IKTB:IKTE) &
              - ZLS(KIB:KIE,KJB:KJE,IKTB:IKTE)*PRI(KIB:KIE,KJB:KJE,IKTB:IKTE)/ZCPD(KIB:KIE,KJB:KJE,IKTB:IKTE)
  ! Determine tropopause/inversion  height from minimum temperature
  ITPL(KIB:KIE,KJB:KJE)  = KKB+KKL
  ZTMIN(KIB:KIE,KJB:KJE) = 400.
  !$acc loop seq
  DO JK = IKTB+1,IKTE-1
     WHERE ( PT(KIB:KIE,KJB:KJE,JK) < ZTMIN(KIB:KIE,KJB:KJE) ) 
          ZTMIN(KIB:KIE,KJB:KJE) = PT(KIB:KIE,KJB:KJE,JK)
          ITPL(KIB:KIE,KJB:KJE) = JK
     END WHERE
  END DO
  ! Set the mixing length scale
  ZL(KIB:KIE,KJB:KJE,KKB) = 20.
  !$acc loop seq
  DO JK = KKB+KKL,KKE,KKL
        ! free troposphere
        ZL(KIB:KIE,KJB:KJE,JK) = ZL0
        ZZZP(KIB:KIE,KJB:KJE) =  PZZ(KIB:KIE,KJB:KJE,JK) -  PZZ(KIB:KIE,KJB:KJE,KKB)
        JKPP(KIB:KIE,KJB:KJE) = ITPL(KIB:KIE,KJB:KJE)
        ! approximate length for boundary-layer
        WHERE ( ZL0 > ZZZP(KIB:KIE,KJB:KJE))
           ZL(KIB:KIE,KJB:KJE,JK) = ZZZP(KIB:KIE,KJB:KJE)
        END WHERE
        ! gradual decrease of length-scale near and above tropopause
        DO CONCURRENT (JI=KIB:KIE,JJ=KJB:KJE)
        IF ( ZZZP(JI,JJ) > 0.9*(PZZ(JI,JJ,JKPP(JI,JJ))-PZZ(JI,JJ,KKB)) ) &
           ZL(JI,JJ,JK) = .6 * ZL(JI,JJ,JK-KKL)
        END DO
  END DO
!$acc end kernels
END IF
!-------------------------------------------------------------------------------
!
!
!$acc kernels
!Ice fraction
ZFRAC(:,:,:) = 0.
IF (OUSERI) THEN
  !$mnh_do_concurrent ( JI = KIB : KIE, JJ = KJB : KJE, JK = KKB : KKE )
    IF ( PRC(JI,JJ,JK) + PRI(JI,JJ,JK) > 1.E-20 ) ZFRAC(JI,JJ,JK) = PRI(JI,JJ,JK) / ( PRC(JI,JJ,JK) + PRI(JI,JJ,JK) )
  !$mnh_end_do()
ENDIF
!$acc end kernels
!
#ifndef MNH_OPENACC
IF (OUSERI) CALL COMPUTE_FRAC_ICE(HFRAC_ICE, ZFRAC, PT)
#else
IF (OUSERI) CALL COMPUTE_FRAC_ICE3D_DEVICE(HFRAC_ICE, ZFRAC, PT)
#endif
!
!$acc kernels
!$mnh_do_concurrent ( JK=IKTB:IKTE )
  JKPK(JK)=MAX(MIN(JK+KKL,IKTE),IKTB)
  JKMK(JK)=MAX(MIN(JK-KKL,IKTE),IKTB)
!$mnh_end_do()
! Bypass CCE/14++ compile bug with index in the good order !!!
#if defined(MNH_COMPILER_CCE) && defined(MNH_BITREP_OMP)
DO CONCURRENT(JI=KIB:KIE,JJ=KJB:KJE,JK=IKTB:IKTE)
#else
!$mnh_do_concurrent(JI=KIB:KIE,JJ=KJB:KJE,JK=IKTB:IKTE)
#endif
   ! latent heats
      ! saturated water vapor mixing ratio over liquid water
#if !defined(MNH_BITREP) && !defined(MNH_BITREP_OMP)
      ZPVP(JI,JJ,JK)    = MIN(EXP( XALPW - XBETAW / PT(JI,JJ,JK) - XGAMW * LOG( PT(JI,JJ,JK) ) ), .99*PPABS(JI,JJ,JK))
#else
      ZPVP(JI,JJ,JK)    = MIN(BR_EXP( XALPW - XBETAW / PT(JI,JJ,JK) - XGAMW * BR_LOG( PT(JI,JJ,JK) ) ), .99*PPABS(JI,JJ,JK))
#endif
      ZQSLP(JI,JJ,JK)   = XRD / XRV * ZPVP(JI,JJ,JK) / ( PPABS(JI,JJ,JK) - ZPVP(JI,JJ,JK) )

      ! saturated water vapor mixing ratio over ice
#if !defined(MNH_BITREP) && !defined(MNH_BITREP_OMP)
      ZPIVP(JI,JJ,JK)   = MIN(EXP( XALPI - XBETAI / PT(JI,JJ,JK) - XGAMI * LOG( PT(JI,JJ,JK) ) ), .99*PPABS(JI,JJ,JK))
#else
      ZPIVP(JI,JJ,JK)   = MIN(BR_EXP( XALPI - XBETAI / PT(JI,JJ,JK) - XGAMI * BR_LOG( PT(JI,JJ,JK) ) ), .99*PPABS(JI,JJ,JK))
#endif
      ZQSIP(JI,JJ,JK)   = XRD / XRV * ZPIVP(JI,JJ,JK) / ( PPABS(JI,JJ,JK) - ZPIVP(JI,JJ,JK) )

      ! interpolate between liquid and solid as function of temperature
      ZQSLP(JI,JJ,JK) = (1. - ZFRAC(JI,JJ,JK)) * ZQSLP(JI,JJ,JK) + ZFRAC(JI,JJ,JK) * ZQSIP(JI,JJ,JK)
      ZLVSP(JI,JJ,JK) = (1. - ZFRAC(JI,JJ,JK)) * ZLV(JI,JJ,JK) + &
             & ZFRAC(JI,JJ,JK)      * ZLS(JI,JJ,JK)

      ! coefficients a and b
#if !defined(MNH_BITREP) && !defined(MNH_BITREP_OMP)
      ZAHP(JI,JJ,JK)  = ZLVSP(JI,JJ,JK) * ZQSLP(JI,JJ,JK) / ( XRV * PT(JI,JJ,JK)**2 ) * (XRV * ZQSLP(JI,JJ,JK) / XRD + 1.)
#else
      ZAHP(JI,JJ,JK)  = ZLVSP(JI,JJ,JK) * ZQSLP(JI,JJ,JK) / ( XRV * BR_P2(PT(JI,JJ,JK)) ) * (XRV * ZQSLP(JI,JJ,JK) / XRD + 1.)
#endif
      ZAP(JI,JJ,JK)   = 1. / ( 1. + ZLVSP(JI,JJ,JK)/ZCPD(JI,JJ,JK) * ZAHP(JI,JJ,JK) )
      ZBP(JI,JJ,JK)   = ZAHP(JI,JJ,JK) * ZAP(JI,JJ,JK)

      ZSBARP(JI,JJ,JK) = ZAP(JI,JJ,JK) * ( ZRT(JI,JJ,JK) - ZQSLP(JI,JJ,JK) + &
                   & ZAHP(JI,JJ,JK) * ZLVSP(JI,JJ,JK) * (PRC(JI,JJ,JK)+PRI(JI,JJ,JK)) / ZCPD(JI,JJ,JK))

      ! switch to take either present computed value of SIGMAS
      ! or that of Meso-NH turbulence scheme
      IF ( OSIGMAS ) THEN
        IF (PSIGQSAT/=0.) THEN
#if !defined(MNH_BITREP) && !defined(MNH_BITREP_OMP)
          ZSIGMAP(JI,JJ,JK) = SQRT((2*PSIGS(JI,JJ,JK))**2 + (PSIGQSAT*ZQSLP(JI,JJ,JK)*ZAP(JI,JJ,JK))**2)
#else
          ZSIGMAP(JI,JJ,JK) = BR_POW(BR_P2(2*PSIGS(JI,JJ,JK)) + BR_P2(PSIGQSAT*ZQSLP(JI,JJ,JK)*ZAP(JI,JJ,JK)) , 0.5)
#endif
        ELSE
          ZSIGMAP(JI,JJ,JK) = 2*PSIGS(JI,JJ,JK)
        END IF
      ELSE
        ! parameterize Sigma_s with first_order closure
        ZDZZP(JI,JJ,JK)   =  PZZ(JI,JJ,JKPK(JK)) - PZZ(JI,JJ,JKMK(JK))
        ZDRWP(JI,JJ,JK)   =  ZRT(JI,JJ,JKPK(JK)) - ZRT(JI,JJ,JKMK(JK))
        ZDTLP(JI,JJ,JK)   =  ZTLK(JI,JJ,JKPK(JK)) - ZTLK(JI,JJ,JKMK(JK)) + XG/ZCPD(JI,JJ,JK) * ZDZZP(JI,JJ,JK)
        ZLLP(JI,JJ,JK) = ZL(JI,JJ,JK)
        ! standard deviation due to convection
        ZSIG_CONVP(JI,JJ,JK) =0.
        IF( SIZE(PMFCONV) /= 0) &
             ZSIG_CONVP(JI,JJ,JK) = ZCSIG_CONV * PMFCONV(JI,JJ,JK) / ZAP(JI,JJ,JK)
        ! zsigma should be of order 4.e-4 in lowest 5 km of atmosphere
        ZSIGMAP(JI,JJ,JK) =  SQRT( MAX( 1.E-25, ZCSIGMA * ZCSIGMA * ZLLP(JI,JJ,JK)*ZLLP(JI,JJ,JK) &
                            /(ZDZZP(JI,JJ,JK)*ZDZZP(JI,JJ,JK))*(&
                            ZAP(JI,JJ,JK)*ZAP(JI,JJ,JK)*ZDRWP(JI,JJ,JK)*ZDRWP(JI,JJ,JK) - 2.*ZAP(JI,JJ,JK)*ZBP(JI,JJ,JK) &
                            * ZDRWP(JI,JJ,JK)*ZDTLP(JI,JJ,JK) &
                            + ZBP(JI,JJ,JK)*ZBP(JI,JJ,JK)*ZDTLP(JI,JJ,JK)*ZDTLP(JI,JJ,JK)) + &
                            ZSIG_CONVP(JI,JJ,JK) * ZSIG_CONVP(JI,JJ,JK) ) )
      END IF
      ZSIGMAP(JI,JJ,JK)= MAX( 1.E-10, ZSIGMAP(JI,JJ,JK) )
!     ZSIGMAP(JI,JJ,JK)= MAX( 1.E-12, ZSIGMAP(JI,JJ,JK) )

      ! normalized saturation deficit
      ZQ1P(JI,JJ,JK)   = ZSBARP(JI,JJ,JK)/ZSIGMAP(JI,JJ,JK)

      IF(HCONDENS == 'GAUS')THEN
        ! Gaussian Probability Density Function around ZQ1
        ! Computation of ZG and ZGAM(=erf(ZG))
        ZGCONDP(JI,JJ,JK) = -ZQ1P(JI,JJ,JK)/SQRT(2.)

        !Approximation of erf function for Gaussian distribution
        ZGAUVP(JI,JJ,JK) = 1 - SIGN(1., ZGCONDP(JI,JJ,JK)) * SQRT(1-EXP(-4*ZGCONDP(JI,JJ,JK)**2/XPI))

        !Computation Cloud Fraction
        PCLDFR(JI,JJ,JK) = MAX( 0., MIN(1.,0.5*ZGAUVP(JI,JJ,JK)))

        !Computation of condensate
        ZCONDP(JI,JJ,JK) = (EXP(-ZGCONDP(JI,JJ,JK)**2)-ZGCONDP(JI,JJ,JK)*SQRT(XPI)*ZGAUVP(JI,JJ,JK))*ZSIGMAP(JI,JJ,JK)/SQRT(2.*XPI)
        ZCONDP(JI,JJ,JK) = MAX(ZCONDP(JI,JJ,JK), 0.)

        PSIGRC(JI,JJ,JK) = PCLDFR(JI,JJ,JK)

        !Computation warm/cold Cloud Fraction and content in high water content part
        IF(GPRESENT_PHLC_HCF .AND. GPRESENT_PHLC_HRC)THEN
          IF(1-ZFRAC(JI,JJ,JK) > 1.E-20)THEN
            ZAUTCP(JI,JJ,JK) = (ZSBARP(JI,JJ,JK) - XCRIAUTC/(PRHODREF(JI,JJ,JK)*(1-ZFRAC(JI,JJ,JK))))/ZSIGMAP(JI,JJ,JK)
            ZGAUTCP(JI,JJ,JK) = -ZAUTCP(JI,JJ,JK)/SQRT(2.)
            !Approximation of erf function for Gaussian distribution
            ZGAUCP(JI,JJ,JK) = 1 - SIGN(1., ZGAUTCP(JI,JJ,JK)) * SQRT(1-EXP(-4*ZGAUTCP(JI,JJ,JK)**2/XPI))
            PHLC_HCF(JI,JJ,JK) = MAX( 0., MIN(1.,0.5*ZGAUCP(JI,JJ,JK)))
            PHLC_HRC(JI,JJ,JK) = (1-ZFRAC(JI,JJ,JK))*(EXP(-ZGAUTCP(JI,JJ,JK)**2)-ZGAUTCP(JI,JJ,JK) &
                                 *SQRT(XPI)*ZGAUCP(JI,JJ,JK))*ZSIGMAP(JI,JJ,JK)/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
        ENDIF

        IF(GPRESENT_PHLI_HCF .AND. GPRESENT_PHLI_HRI)THEN
          IF(ZFRAC(JI,JJ,JK) > 1.E-20)THEN
            ZCRIAUTIP(JI,JJ,JK)=MIN(XCRIAUTI,10**(XACRIAUTI*(PT(JI,JJ,JK)-XTT)+XBCRIAUTI))
            ZAUTIP(JI,JJ,JK) = (ZSBARP(JI,JJ,JK) - ZCRIAUTIP(JI,JJ,JK)/ZFRAC(JI,JJ,JK))/ZSIGMAP(JI,JJ,JK)
            ZGAUTIP(JI,JJ,JK) = -ZAUTIP(JI,JJ,JK)/SQRT(2.)
            !Approximation of erf function for Gaussian distribution
            ZGAUIP(JI,JJ,JK) = 1 - SIGN(1., ZGAUTIP(JI,JJ,JK)) * SQRT(1-EXP(-4*ZGAUTIP(JI,JJ,JK)**2/XPI))
            PHLI_HCF(JI,JJ,JK) = MAX( 0., MIN(1.,0.5*ZGAUIP(JI,JJ,JK)))
            PHLI_HRI(JI,JJ,JK) = ZFRAC(JI,JJ,JK)*(EXP(-ZGAUTIP(JI,JJ,JK)**2)-ZGAUTIP(JI,JJ,JK) &
                                 * SQRT(XPI)*ZGAUIP(JI,JJ,JK))*ZSIGMAP(JI,JJ,JK)/SQRT(2.*XPI)
            PHLI_HRI(JI,JJ,JK) = PHLI_HRI(JI,JJ,JK) + ZCRIAUTIP(JI,JJ,JK)*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
        ENDIF

      ELSEIF(HCONDENS == 'CB02')THEN
        !Cloud fraction
#if !defined(MNH_BITREP) && !defined(MNH_BITREP_OMP)
        PCLDFR(JI,JJ,JK) = MAX( 0., MIN(1.,0.5+0.36*ATAN(1.55*ZQ1P(JI,JJ,JK))) )
#else
        PCLDFR(JI,JJ,JK) = MAX( 0., MIN(1.,0.5+0.36*BR_ATAN(1.55*ZQ1P(JI,JJ,JK))) )
#endif

        !Total condensate
        IF (ZQ1P(JI,JJ,JK) > 0. .AND. ZQ1P(JI,JJ,JK) <= 2 ) THEN
#if !defined(MNH_BITREP) && !defined(MNH_BITREP_OMP)
          ZCONDP(JI,JJ,JK) = MIN(EXP(-1.)+.66*ZQ1P(JI,JJ,JK)+.086*ZQ1P(JI,JJ,JK)**2, 2.) ! We use the MIN function for continuity
#else
          ZCONDP(JI,JJ,JK) = MIN(BR_EXP(-1.)+.66*ZQ1P(JI,JJ,JK)+.086*BR_P2(ZQ1P(JI,JJ,JK)), 2.) ! We use the MIN function for continuity
#endif
        ELSE IF (ZQ1P(JI,JJ,JK) > 2.) THEN
          ZCONDP(JI,JJ,JK) = ZQ1P(JI,JJ,JK)
        ELSE
#if !defined(MNH_BITREP) && !defined(MNH_BITREP_OMP)
          ZCONDP(JI,JJ,JK) = EXP( 1.2*ZQ1P(JI,JJ,JK)-1. )
#else
          ZCONDP(JI,JJ,JK) = BR_EXP( 1.2*ZQ1P(JI,JJ,JK)-1. )
#endif
        END IF
        ZCONDP(JI,JJ,JK) = ZCONDP(JI,JJ,JK) * ZSIGMAP(JI,JJ,JK)
        
        INQ1P(JI,JJ,JK) = MIN( MAX(-22,FLOOR(MIN(100., MAX(-100., 2*ZQ1P(JI,JJ,JK)))) ), 10)  !inner min/max prevents sigfpe when 2*zq1 does not fit into an int
        ZINCP(JI,JJ,JK) = 2.*ZQ1P(JI,JJ,JK) - INQ1P(JI,JJ,JK)

        PSIGRC(JI,JJ,JK) =  MIN(1.,(1.-ZINCP(JI,JJ,JK))*ZSRC_1D(INQ1P(JI,JJ,JK))+ZINCP(JI,JJ,JK)*ZSRC_1D(INQ1P(JI,JJ,JK)+1))

        IF(GPRESENT_PHLC_HCF .AND. GPRESENT_PHLC_HRC)THEN
          PHLC_HCF(JI,JJ,JK)=0.
          PHLC_HRC(JI,JJ,JK)=0.
        ENDIF
        IF(GPRESENT_PHLI_HCF .AND. GPRESENT_PHLI_HRI)THEN
          PHLI_HCF(JI,JJ,JK)=0.
          PHLI_HRI(JI,JJ,JK)=0.
        ENDIF
      ENDIF

      IF ( ZCONDP(JI,JJ,JK) < 1.E-12 ) THEN
        ZCONDP(JI,JJ,JK) = 0.
        PCLDFR(JI,JJ,JK) = 0.
      ENDIF
      IF (PCLDFR(JI,JJ,JK)==0.) THEN
        ZCONDP(JI,JJ,JK)=0.
      ENDIF

      ZRCOLDP(JI,JJ,JK)=PRC(JI,JJ,JK)
      ZRIOLDP(JI,JJ,JK)=PRI(JI,JJ,JK)

      PRC(JI,JJ,JK) = (1.-ZFRAC(JI,JJ,JK)) * ZCONDP(JI,JJ,JK) ! liquid condensate
      PRI(JI,JJ,JK) = ZFRAC(JI,JJ,JK) * ZCONDP(JI,JJ,JK)   ! solid condensate
      
      PT(JI,JJ,JK) = PT(JI,JJ,JK) + ((PRC(JI,JJ,JK)-ZRCOLDP(JI,JJ,JK))*ZLV(JI,JJ,JK) + &
                                    &(PRI(JI,JJ,JK)-ZRIOLDP(JI,JJ,JK))*ZLS(JI,JJ,JK)   ) &
                                  & /ZCPD(JI,JJ,JK)
      PRV(JI,JJ,JK) = ZRT(JI,JJ,JK) - PRC(JI,JJ,JK) - PRI(JI,JJ,JK)

! 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.-ZQ1P(JI,JJ,JK)) )
! in the 3D case lambda_3 = 1.

      IF(GYLAMBDA3_CB)THEN
        PSIGRC(JI,JJ,JK) = PSIGRC(JI,JJ,JK)* MIN( 3. , MAX(1.,1.-ZQ1P(JI,JJ,JK)) )
      ENDIF
#if defined(MNH_COMPILER_CCE) && defined(MNH_BITREP_OMP)
ENDDO ! CONCURRENT
#else
!$mnh_end_do() !CONCURRENT
#endif
!$acc end kernels

!$acc end data

#ifndef MNH_OPENACC
deallocate( ztlk, zrt, zl, zfrac, itpl, ztmin, zlv, zls,zcpd, gwork )
deallocate( ZZZP, JKPP, ZPVP, ZQSLP, ZPIVP, ZQSIP, ZLVSP, ZAHP, ZAP, ZBP, ZSBARP, ZSIGMAP )
deallocate( ZDZZP, ZDRWP, ZDTLP, ZLLP, ZSIG_CONVP, ZQ1P, ZGCONDP, ZGAUVP, ZCONDP, ZAUTCP )
deallocate( ZGAUTCP, ZGAUCP, ZCRIAUTIP, ZAUTIP, ZGAUTIP, ZGAUIP, INQ1P, ZINCP, ZRCOLDP, ZRIOLDP )
deallocate( JKPK,JKMK )
#else
!Release all memory allocated with MNH_MEM_GET calls since last call to MNH_MEM_POSITION_PIN
CALL MNH_MEM_RELEASE( 'CONDENSATION' )
#endif

IF (MPPDB_INITIALIZED) THEN
  !Check all INOUT arrays
  CALL MPPDB_CHECK3D(PT,"CONDENSATION end:PT",PRECISION)
  CALL MPPDB_CHECK3D(PRV,"CONDENSATION end:PRV",PRECISION)
  CALL MPPDB_CHECK3D(PRC,"CONDENSATION end:PRC",PRECISION)
  CALL MPPDB_CHECK3D(PRI,"CONDENSATION end:PRI",PRECISION)
  !Check all OUT arrays
  CALL MPPDB_CHECK3D(PCLDFR,"CONDENSATION end:PCLDFR",PRECISION)
  CALL MPPDB_CHECK3D(PSIGRC,"CONDENSATION end:PSIGRC",PRECISION)
END IF

END SUBROUTINE CONDENSATION