diff --git a/docs/TODO b/docs/TODO
index 2d1ac4b2bc01a4adb164dde29b6dc0339678d331..65fbe237e019935d57a2c54592276d3514b79996 100644
--- a/docs/TODO
+++ b/docs/TODO
@@ -10,10 +10,10 @@ Dependencies:
- liste dans document Interfaces
- pour AROME placés, en attendant, dans phyex/externals
-Clé de compilation REPRO48 ajoutée pour permettre de reproduire le cycle 48, elle:
-- contourne des corrections de bug
-- modifie l'organisation de calculs
-Cette clé devra être supprimée
+Clé de compilation REPRO48 + REPRO55 ajoutées pour permettre de reproduire le cycle 48 MNH-5.5.0, elles:
+- contournent des corrections de bug
+- modifient l'organisation de calculs
+Ces clés devront être supprimées
Merge pb:
- ice4_nucleation_wrapper:
@@ -45,7 +45,7 @@ Pb identifiés à corriger plus tard:
=> Ã corriger
- seules les options oper ont été testées, il manque des test pour sedim_after, nmaxiter, xmrstep, xtstep, autoconv, rainfr
-Répertoire arome/ext contient les codes non PHYEX qu'il faut modifier dans le pack pour qu'il puisse être compilé.
+Répertoire arome/ext et mesonh/ext contient les codes non PHYEX qu'il faut modifier dans le pack pour qu'il puisse être compilé.
Ce répertoire devra être vidé à la fin du phasage, les modifications nécessaires ayadevront avoir été fournies par ailleurs
Budgets/DDH
diff --git a/src/common/micro/condensation.F90 b/src/common/micro/condensation.F90
index 2e8112fc2fa0a54d8fbca37bc3fca5b4a82eab8c..49cf4e8964b4071a533139a323c902e939ea1e0a 100644
--- a/src/common/micro/condensation.F90
+++ b/src/common/micro/condensation.F90
@@ -173,7 +173,7 @@ 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
+INTEGER, DIMENSION(KIB:KIE) :: IERR
!
!* 0.3 Definition of constants :
!
@@ -207,6 +207,12 @@ PCLDFR(:,:,:) = 0. ! Initialize values
PSIGRC(:,:,:) = 0. ! Initialize values
ZPRIFACT = 1. ! Initialize value
ZCLDUM=-1. ! Initialize value
+! Init of the HALO (should be on HALO points only)
+#ifdef REPRO55
+PRC_OUT = PRC_IN
+PRV_OUT = PRV_IN
+PRI_OUT = PRI_IN
+#endif
IF(OCND2)ZPRIFACT = 0.
!
!
@@ -325,7 +331,7 @@ DO JK=IKTB,IKTE
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
+ CALL COMPUTE_FRAC_ICE(HFRAC_ICE, ZFRAC(:), PT(KIB:KIE,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) )
diff --git a/src/common/micro/ice_adjust.F90 b/src/common/micro/ice_adjust.F90
index 355cb63ccc4fc7de4d12decc5433d2a6b9bc094a..5196f0b9414f431fb48224bda5d180d2d93aebce 100644
--- a/src/common/micro/ice_adjust.F90
+++ b/src/common/micro/ice_adjust.F90
@@ -11,10 +11,9 @@
PPABST, PZZ, &
PEXN, PCF_MF, PRC_MF, PRI_MF, &
PRV, PRC, PRVS, PRCS, PTH, PTHS, PSRCS, PCLDFR, &
- PRR, PRI, PRIS, PRS, PRG, PRH, &
+ PRR, PRI, PRIS, PRS, PRG, TBUDGETS, KBUDGETS, PRH,&
POUT_RV, POUT_RC, POUT_RI, POUT_TH, &
PHLC_HRC, PHLC_HCF, PHLI_HRI, PHLI_HCF, &
- TBUDGETS, KBUDGETS, &
PICE_CLD_WGT)
! #########################################################################
!
diff --git a/src/common/micro/mode_icecloud.F90 b/src/common/micro/mode_icecloud.F90
index 43b394d211255e7ef887cecec79246067a4ac324..774dd688734a27d00b2fa1caaebe32d83ad07011 100644
--- a/src/common/micro/mode_icecloud.F90
+++ b/src/common/micro/mode_icecloud.F90
@@ -11,7 +11,8 @@ SUBROUTINE ICECLOUD &
USE YOMHOOK , ONLY : LHOOK, DR_HOOK
USE MODD_CST,ONLY : XCPD,XCPV,XLVTT,XLSTT,XG,XRD,XEPSILO
USE MODE_TIWMX, ONLY: ESATW, ESATI
- USE MODI_TIWMX
+ USE MODE_QSATMX_TAB
+! USE MODI_TIWMX
IMPLICIT NONE
!-----------------------------------------------------------------------
!
@@ -77,7 +78,7 @@ REAL :: ZSIGMAX,ZSIGMAY,ZSIGMAZ,ZXDIST,ZYDIST,&
INTEGER :: JK
! External function
-REAL :: QSATMX_TAB
+!REAL :: QSATMX_TAB
REAL(KIND=JPRB) :: ZHOOK_HANDLE
IF (LHOOK) CALL DR_HOOK('ICECLOUD',0,ZHOOK_HANDLE)
diff --git a/src/common/micro/modi_condensation.F90 b/src/common/micro/modi_condensation.F90
index 7163e82088c944853c5f5f65a30a6f8eeac923fa..4f1a548de37a6b80b26b50ac2eb705fe0eafe9b2 100644
--- a/src/common/micro/modi_condensation.F90
+++ b/src/common/micro/modi_condensation.F90
@@ -12,7 +12,6 @@ INTERFACE
PLV, PLS, PCPH, &
PHLC_HRC, PHLC_HCF, PHLI_HRI, PHLI_HCF, PICE_CLD_WGT)
!
-USE MODD_SPP_TYPE
!
INTEGER, INTENT(IN) :: KIU ! horizontal dimension in x
INTEGER, INTENT(IN) :: KJU ! horizontal dimension in y
diff --git a/src/common/micro/modi_ice_adjust.F90 b/src/common/micro/modi_ice_adjust.F90
index e949147f50114a2a109eb8b12ec78096fad83ffc..147a971cde73d72b10db6551860a4ff42c1fa037 100644
--- a/src/common/micro/modi_ice_adjust.F90
+++ b/src/common/micro/modi_ice_adjust.F90
@@ -11,10 +11,9 @@ INTERFACE
PPABST, PZZ, &
PEXN, PCF_MF, PRC_MF, PRI_MF, &
PRV, PRC, PRVS, PRCS, PTH, PTHS, PSRCS, PCLDFR, &
- PRR, PRI, PRIS, PRS, PRG, PRH, &
+ PRR, PRI, PRIS, PRS, PRG, TBUDGETS, KBUDGETS, PRH,&
POUT_RV, POUT_RC, POUT_RI, POUT_TH, &
PHLC_HRC, PHLC_HCF, PHLI_HRI, PHLI_HCF, &
- TBUDGETS, KBUDGETS, &
PICE_CLD_WGT)
USE MODD_BUDGET, ONLY: TBUDGETDATA
IMPLICIT NONE
diff --git a/src/mesonh/micro/resolved_cloud.f90 b/src/mesonh/ext/resolved_cloud.f90
similarity index 100%
rename from src/mesonh/micro/resolved_cloud.f90
rename to src/mesonh/ext/resolved_cloud.f90
diff --git a/src/mesonh/micro/ice_adjust_elec.f90 b/src/mesonh/micro/ice_adjust_elec.f90
new file mode 100644
index 0000000000000000000000000000000000000000..a52ffd2b375acda4dd911a74577409144202de0c
--- /dev/null
+++ b/src/mesonh/micro/ice_adjust_elec.f90
@@ -0,0 +1,651 @@
+!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.
+!-----------------------------------------------------------------
+! ###########################
+ MODULE MODI_ICE_ADJUST_ELEC
+! ###########################
+!
+INTERFACE
+!
+ SUBROUTINE ICE_ADJUST_ELEC (KRR, KMI, HRAD, HTURBDIM, HSCONV, HMF_CLOUD, &
+ OSUBG_COND, OSIGMAS, PTSTEP,PSIGQSAT, &
+ PRHODJ, PEXNREF, PSIGS, PPABST, PZZ, &
+ PMFCONV, PCF_MF, PRC_MF, PRI_MF, &
+ PRVT, PRCT, PRVS, PRCS, PTHS, PSRCS, PCLDFR , &
+ PRRT, PRRS, PRIT, PRIS, PRST, PRSS, PRGT, PRGS, &
+ PQPIT, PQPIS, PQCT, PQCS, &
+ PQRT, PQRS, PQIT, PQIS, PQST, PQSS, PQGT, PQGS, &
+ PQNIT, PQNIS, PRHT, PRHS, PQHT, PQHS )
+!
+INTEGER, INTENT(IN) :: KRR ! Number of moist variables
+INTEGER, INTENT(IN) :: KMI ! Model index
+CHARACTER(len=4), INTENT(IN) :: HTURBDIM ! Dimensionality of the
+ ! turbulence scheme
+CHARACTER(LEN=4), INTENT(IN) :: HSCONV ! Shallow convection scheme
+CHARACTER(LEN=4), INTENT(IN) :: HMF_CLOUD! Type of statistical cloud
+CHARACTER(len=4), INTENT(IN) :: HRAD ! Radiation scheme name
+LOGICAL, INTENT(IN) :: OSUBG_COND ! Switch for Subgrid
+ ! Condensation
+LOGICAL :: OSIGMAS ! Switch for Sigma_s:
+ ! use values computed in CONDENSATION
+ ! or that from turbulence scheme
+REAL, INTENT(IN) :: PTSTEP ! Double Time step
+ ! (single if cold start)
+REAL, INTENT(IN) :: PSIGQSAT ! coeff applied to qsat variance contribution
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODJ ! Dry density * Jacobian
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PEXNREF ! Reference Exner function
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PSIGS ! Sigma_s at time t
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PMFCONV ! convective mass flux
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PPABST ! Absolute Pressure at t
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PZZ ! height of model layer
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PCF_MF! Convective Mass Flux Cloud fraction
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRC_MF! Convective Mass Flux liquid mixing ratio
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRI_MF! Convective Mass Flux solid mixing ratio
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRVT ! Water vapor m.r. at t
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRCT ! Cloud water m.r. at t
+REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PRVS ! Water vapor m.r. source
+REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PRCS ! Cloud water m.r. source
+REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PTHS ! Theta source
+REAL, DIMENSION(:,:,:), INTENT(OUT) :: PSRCS ! Second-order flux
+ ! s'rc'/2Sigma_s2 at time t+1
+ ! multiplied by Lambda_3
+REAL, DIMENSION(:,:,:), INTENT(OUT) :: PCLDFR ! Cloud fraction
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRRS ! Rain water m.r. at t+1
+REAL, DIMENSION(:,:,:), INTENT(INOUT):: PRIS ! Cloud ice m.r. at t+1
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRSS ! Aggregate m.r. at t+1
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRGS ! Graupel m.r. at t+1
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRRT ! Rain water m.r. at t
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRIT ! Cloud ice m.r. at t
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRST ! Aggregate m.r. at t
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRGT ! Graupel m.r. at t
+REAL, DIMENSION(:,:,:), OPTIONAL, INTENT(IN) :: PRHT ! Hail m.r. at t
+REAL, DIMENSION(:,:,:), OPTIONAL, INTENT(IN) :: PRHS ! Hail m.r. at t+1
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PQPIT ! positive ion m.r. at t
+REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PQPIS ! source
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PQNIT ! negative ion m.r. at t
+REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PQNIS ! source
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PQCT ! Cloud water m.r. at t
+REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PQCS ! Cloud water m.r. source
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PQRS ! Rain water m.r. at t+1
+REAL, DIMENSION(:,:,:), INTENT(INOUT):: PQIS ! Cloud ice m.r. at t+1
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PQSS ! Aggregate m.r. at t+1
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PQGS ! Graupel m.r. at t+1
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PQRT ! Rain water m.r. at t
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PQIT ! Cloud ice m.r. at t
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PQST ! Aggregate m.r. at t
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PQGT ! Graupel m.r. at t
+REAL, DIMENSION(:,:,:), OPTIONAL, INTENT(IN) :: PQHT ! Hail m.r. at t
+REAL, DIMENSION(:,:,:), OPTIONAL, INTENT(IN) :: PQHS ! Hail m.r. at t+1
+!
+END SUBROUTINE ICE_ADJUST_ELEC
+END INTERFACE
+END MODULE MODI_ICE_ADJUST_ELEC
+!
+! ########################################################################
+ SUBROUTINE ICE_ADJUST_ELEC (KRR, KMI, HRAD, HTURBDIM, HSCONV, &
+ HMF_CLOUD, OSUBG_COND, OSIGMAS, PTSTEP,PSIGQSAT,&
+ PRHODJ, PEXNREF, PSIGS, PPABST, PZZ, &
+ PMFCONV, PCF_MF, PRC_MF, PRI_MF, &
+ PRVT, PRCT, PRVS, PRCS, PTHS, PSRCS, PCLDFR , &
+ PRRT, PRRS, PRIT, PRIS, PRST, PRSS, PRGT, PRGS, &
+ PQPIT, PQPIS, PQCT, PQCS, &
+ PQRT, PQRS, PQIT, PQIS, PQST, PQSS, PQGT, PQGS, &
+ PQNIT, PQNIS, PRHT, PRHS, PQHT, PQHS )
+! ########################################################################
+!
+!!**** *ICE_ADJUST_ELEC* - compute the ajustment of water vapor in mixed-phase
+!! clouds
+!!
+!! PURPOSE
+!! -------
+!! The purpose of this routine is to compute the fast microphysical sources
+!! through a saturation ajustement procedure in case of mixed-phase clouds.
+!!
+!!
+!!** METHOD
+!! ------
+!! Langlois, Tellus, 1973 for the cloudless version.
+!! When cloud water is taken into account, refer to book 1 of the
+!! documentation.
+!!
+!!
+!! EXTERNAL
+!! --------
+!! None
+!!
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! Module MODD_CST
+!! XP00 ! Reference pressure
+!! XMD,XMV ! Molar mass of dry air and molar mass of vapor
+!! XRD,XRV ! Gaz constant for dry air, gaz constant for vapor
+!! XCPD,XCPV ! Cpd (dry air), Cpv (vapor)
+!! XCL ! Cl (liquid)
+!! XCI ! Ci (ice)
+!! XTT ! Triple point temperature
+!! XLVTT ! Vaporization heat constant
+!! XLSTT ! Sublimation heat constant
+!! XALPW,XBETAW,XGAMW ! Constants for saturation vapor over liquid
+!! ! pressure function
+!! XALPI,XBETAI,XGAMI ! Constants for saturation vapor over ice
+!! ! pressure function
+!! Module MODD_CONF
+!! CCONF
+!!
+!!
+!! REFERENCE
+!! ---------
+!! Book 1 and Book2 of documentation ( routine ICE_ADJUST )
+!! Langlois, Tellus, 1973
+!!
+!! AUTHOR
+!! ------
+!! J.-P. Pinty * Laboratoire d'Aerologie*
+!!
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 2002
+!! C. Barthe 19/11/09 update to version 4.8.1
+!! M. Chong Mar. 2010 Add small ions
+!! J.Escobar : 15/09/2015 : WENO5 & JPHEXT <> 1
+! P. Wautelet 05/2016-04/2018: new data structures and calls for I/O
+! P. Wautelet 03/2020: use the new data structures and subroutines for budgets
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+use modd_budget, only: lbudget_th, lbudget_rv, lbudget_rc, lbudget_ri, lbudget_sv, &
+ NBUDGET_TH, NBUDGET_RV, NBUDGET_RC, NBUDGET_RI, NBUDGET_SV1, &
+ tbudgets
+USE MODD_CONF
+USE MODD_CST
+USE MODD_ELEC_DESCR, ONLY : XRTMIN_ELEC, XQTMIN, XFC, XFI, XECHARGE
+USE MODD_NSV, ONLY : NSV_ELECBEG, NSV_ELECEND
+USE MODD_PARAMETERS
+USE MODD_RAIN_ICE_DESCR, ONLY : XRTMIN, XBI
+
+use mode_budget, only: Budget_store_init, Budget_store_end
+use mode_tools_ll, only: GET_INDICE_ll
+
+USE MODI_CONDENSATION
+USE MODI_GET_HALO
+!
+IMPLICIT NONE
+!
+!
+!* 0.1 Declarations of dummy arguments
+!
+INTEGER, INTENT(IN) :: KRR ! Number of moist variables
+INTEGER, INTENT(IN) :: KMI ! Model index
+CHARACTER(len=4), INTENT(IN) :: HTURBDIM ! Dimensionality of the
+ ! turbulence scheme
+CHARACTER(LEN=4), INTENT(IN) :: HSCONV ! Shallow convection scheme
+CHARACTER(LEN=4), INTENT(IN) :: HMF_CLOUD! Type of statistical cloud
+CHARACTER(len=4), INTENT(IN) :: HRAD ! Radiation scheme name
+LOGICAL, INTENT(IN) :: OSUBG_COND ! Switch for Subgrid
+ ! Condensation
+LOGICAL :: OSIGMAS ! Switch for Sigma_s:
+ ! use values computed in CONDENSATION
+ ! or that from turbulence scheme
+REAL, INTENT(IN) :: PTSTEP ! Double Time step
+ ! (single if cold start)
+REAL, INTENT(IN) :: PSIGQSAT ! coeff applied to qsat variance contribution
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODJ ! Dry density * Jacobian
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PEXNREF ! Reference Exner function
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PSIGS ! Sigma_s at time t
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PMFCONV ! convective mass flux
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PPABST ! Absolute Pressure at t
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PZZ ! height of model layer
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PCF_MF! Convective Mass Flux Cloud fraction
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRC_MF! Convective Mass Flux liquid mixing ratio
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRI_MF! Convective Mass Flux solid mixing ratio
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRVT ! Water vapor m.r. at t
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRCT ! Cloud water m.r. at t
+REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PRVS ! Water vapor m.r. source
+REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PRCS ! Cloud water m.r. source
+REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PTHS ! Theta source
+REAL, DIMENSION(:,:,:), INTENT(OUT) :: PSRCS ! Second-order flux
+ ! s'rc'/2Sigma_s2 at time t+1
+ ! multiplied by Lambda_3
+REAL, DIMENSION(:,:,:), INTENT(OUT) :: PCLDFR ! Cloud fraction
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRRS ! Rain water m.r. at t+1
+REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PRIS ! Cloud ice m.r. at t+1
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRSS ! Aggregate m.r. at t+1
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRGS ! Graupel m.r. at t+1
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRRT ! Rain water m.r. at t
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRIT ! Cloud ice m.r. at t
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRST ! Aggregate m.r. at t
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRGT ! Graupel m.r. at t
+REAL, DIMENSION(:,:,:), OPTIONAL, INTENT(IN) :: PRHT ! Hail m.r. at t
+REAL, DIMENSION(:,:,:), OPTIONAL, INTENT(IN) :: PRHS ! Hail m.r. at t+1
+!
+REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PQPIT ! positive ion m.r. at t
+REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PQPIS ! source
+REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PQNIT ! negative ion m.r. at t
+REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PQNIS ! source
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PQCT ! Cloud water m.r. at t
+REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PQCS ! Cloud water m.r. source
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PQRS ! Rain water m.r. at t+1
+REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PQIS ! Cloud ice m.r. at t+1
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PQSS ! Aggregate m.r. at t+1
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PQGS ! Graupel m.r. at t+1
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PQRT ! Rain water m.r. at t
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PQIT ! Cloud ice m.r. at t
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PQST ! Aggregate m.r. at t
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PQGT ! Graupel m.r. at t
+REAL, DIMENSION(:,:,:), OPTIONAL, INTENT(IN) :: PQHT ! Hail m.r. at t
+REAL, DIMENSION(:,:,:), OPTIONAL, INTENT(IN) :: PQHS ! Hail m.r. at t+1
+!
+!
+!* 0.2 Declarations of local variables :
+!
+REAL :: ZEPS ! Mv/Md
+REAL :: ZT00,ZT0 ! Min and max temperature for the mixed phase liquid and solid water
+ ! for the coeff CND of the barycentric mixing ratio
+REAL, DIMENSION(SIZE(PEXNREF,1),SIZE(PEXNREF,2),SIZE(PEXNREF,3)) &
+ :: ZEXNS,& ! guess of the Exner functon at t+1
+ ZT, & ! guess of the temperature at t+1
+ ZCPH, & ! guess of the CPh for the mixing
+ ZLV, & ! guess of the Lv at t+1
+ ZLS, & ! guess of the Ls at t+1
+ ZW1,ZW2,ZW3,ZW4,ZW5,ZW6,ZW7,& ! Work arrays for intermediate fields
+ ZW1_IN, ZW2_IN, ZW3_IN, &
+ ZCND ! CND=(T-T00)/(T0-T00) cf sc doc and TAO etal (89)
+REAL, DIMENSION(SIZE(PEXNREF,1),SIZE(PEXNREF,2),SIZE(PEXNREF,3)) &
+ :: ZWE1, &
+ ZWE2
+REAL, DIMENSION(SIZE(PEXNREF,1),SIZE(PEXNREF,2),SIZE(PEXNREF,3)) &
+ :: ZION_NUMBER, & !nearly Nb of elementary charge
+ ! in hydrometeor charge
+ ZADD ! ratio (0 or 1) of ZION_NUMBER
+ ! to add to positive
+ ! or negative ion number
+REAL, DIMENSION(SIZE(PEXNREF,1),SIZE(PEXNREF,2)) :: ZSIGQSAT2D
+ !
+INTEGER :: IIU,IJU,IKU! dimensions of dummy arrays
+INTEGER :: IIB,IJB ! Horz index values of the first inner mass points
+INTEGER :: IIE,IJE ! Horz index values of the last inner mass points
+INTEGER :: IKB ! K index value of the first inner mass point
+INTEGER :: IKE ! K index value of the last inner mass point
+INTEGER :: JITER,ITERMAX ! iterative loop for first order adjustment
+!
+LOGICAL :: LPRETREATMENT, LNEW_ADJUST
+!
+!-------------------------------------------------------------------------------
+!
+!* 1. PRELIMINARIES
+! -------------
+!
+IIU = SIZE(PEXNREF,1)
+IJU = SIZE(PEXNREF,2)
+IKU = SIZE(PEXNREF,3)
+CALL GET_INDICE_ll (IIB,IJB,IIE,IJE)
+IKB = 1 + JPVEXT
+IKE = IKU - JPVEXT
+!
+ZEPS = XMV / XMD
+!
+ITERMAX = 1
+!
+LPRETREATMENT=.TRUE. ! FALSE to retreive the previous MASDEV4_1 version
+LNEW_ADJUST =.TRUE. ! FALSE to retreive the previous MASDEV4_1 version
+ZT0 = XTT ! Usefull if LPRETREATMENT=T or LNEW_ADJUST=T
+ZT00 = XTT-40. ! Usefull if LPRETREATMENT=T or LNEW_ADJUST=T
+!
+!-------------------------------------------------------------------------------
+if ( lbudget_th ) call Budget_store_init( tbudgets(NBUDGET_TH), 'DEPI', pths(:, :, :) * prhodj(:, :, :) )
+if ( lbudget_rv ) call Budget_store_init( tbudgets(NBUDGET_RV), 'DEPI', prvs(:, :, :) * prhodj(:, :, :) )
+if ( lbudget_rc ) call Budget_store_init( tbudgets(NBUDGET_RC), 'DEPI', prcs(:, :, :) * prhodj(:, :, :) )
+if ( lbudget_ri ) call Budget_store_init( tbudgets(NBUDGET_RI), 'DEPI', pris(:, :, :) * prhodj(:, :, :) )
+if ( lbudget_sv ) then
+ call Budget_store_init( tbudgets(NBUDGET_SV1 - 1 + nsv_elecbeg ), 'DEPI', pqpis(:, :, :) * prhodj(:, :, :) )
+ call Budget_store_init( tbudgets(NBUDGET_SV1 - 1 + nsv_elecend ), 'DEPI', pqnis(:, :, :) * prhodj(:, :, :) )
+ call Budget_store_init( tbudgets(NBUDGET_SV1 - 1 + nsv_elecbeg + 1 ), 'DEPI', pqcs (:, :, :) * prhodj(:, :, :) )
+ call Budget_store_init( tbudgets(NBUDGET_SV1 - 1 + nsv_elecbeg + 3 ), 'DEPI', pqis (:, :, :) * prhodj(:, :, :) )
+end if
+!
+!* 2. COMPUTE QUANTITIES WITH THE GUESS OF THE FUTURE INSTANT
+! -------------------------------------------------------
+!
+!* 2.1 estimate the pressure at t+1
+!
+ZEXNS(:,:,:) = ( PPABST(:,:,:) / XP00)**(XRD/XCPD)
+!
+! beginning of the iterative loop
+!
+DO JITER = 1, ITERMAX
+!
+!* 2.2 compute the intermediate temperature at t+1, T*
+!
+ ZT(:,:,:) = (PTHS(:,:,:) * PTSTEP) * ZEXNS(:,:,:)
+!
+!* 2.3 compute the latent heat of vaporization Lv(T*) at t+1
+! and the latent heat of sublimation Ls(T*) at t+1
+!
+ ZLV(:,:,:) = XLVTT + (XCPV - XCL) * (ZT(:,:,:) - XTT)
+ ZLS(:,:,:) = XLSTT + (XCPV - XCI) * (ZT(:,:,:) - XTT)
+!
+!* 2.4 compute the specific heat for moist air (Cph) at t+1
+!
+ IF ( KRR == 7 ) THEN
+ ZCPH(:,:,:) = XCPD + XCPV *PTSTEP* PRVS(:,:,:) &
+ + XCL *PTSTEP* (PRCS(:,:,:) + PRRS(:,:,:)) &
+ + XCI *PTSTEP* (PRIS(:,:,:) + PRSS(:,:,:) + PRGS(:,:,:) &
+ + PRHS(:,:,:))
+ ELSE IF( KRR == 6 ) THEN
+ ZCPH(:,:,:) = XCPD + XCPV *PTSTEP* PRVS(:,:,:) &
+ + XCL *PTSTEP* (PRCS(:,:,:) + PRRS(:,:,:)) &
+ + XCI *PTSTEP* (PRIS(:,:,:) + PRSS(:,:,:) + PRGS(:,:,:))
+ ELSE IF( KRR == 5 ) THEN
+ ZCPH(:,:,:) = XCPD + XCPV *PTSTEP* PRVS(:,:,:) &
+ + XCL *PTSTEP* (PRCS(:,:,:) + PRRS(:,:,:)) &
+ + XCI *PTSTEP* (PRIS(:,:,:) + PRSS(:,:,:))
+ ELSE IF( KRR == 3 ) THEN
+ ZCPH(:,:,:) = XCPD + XCPV *PTSTEP* PRVS(:,:,:) &
+ + XCL *PTSTEP* (PRCS(:,:,:) + PRRS(:,:,:))
+ ELSE IF( KRR == 2 ) THEN
+ ZCPH(:,:,:) = XCPD + XCPV *PTSTEP* PRVS(:,:,:) &
+ + XCL *PTSTEP* PRCS(:,:,:)
+ END IF
+!
+!
+!* 3. FIRST ORDER SUBGRID CONDENSATION SCHEME
+! ---------------------------------------
+!
+ IF (OSUBG_COND) THEN
+!
+!* 3.1 compute condensate, cloud fraction
+!
+ ! ZW3=water vapor ZW1=rc (OUT) ZW2=ri (OUT) PSRC= s'rci'/Sigma_s^2
+ ! ZW3_IN/ZW2_IN/ZW1_IN (IN)
+ ZW3_IN = PRVS * PTSTEP; ZW1_IN = PRCS * PTSTEP; ZW2_IN = PRIS * PTSTEP
+ ZW3=ZW3_IN; ZW2=ZW2_IN; ZW1=ZW1_IN
+ ZSIGQSAT2D(:,:)=PSIGQSAT
+ ZW4 = 1. ! PRODREF is not used if HL variables are not present
+!
+ CALL CONDENSATION( IIU, IJU, IKU, IIB, IIE, IJB, IJE, IKB, IKE,1, 'T', 'CB02', 'CB', &
+ PPABST, PZZ, ZW4, ZT, ZW3_IN, ZW3, ZW1_IN, ZW1, ZW2_IN, ZW2, &
+ PRSS*PTSTEP, PRGS*PTSTEP, &
+ PSIGS, PMFCONV, PCLDFR, PSRCS, .TRUE., &
+ OSIGMAS, .FALSE., ZSIGQSAT2D, PLV=ZLV, PLS=ZLS, PCPH=ZCPH )
+!
+!* 3.2 compute the variation of mixing ratio
+!
+ ! Rc - Rc*
+ ZW1(:,:,:) = (ZW1(:,:,:) / PTSTEP) - PRCS(:,:,:) ! Pcon = ----------
+ ! 2 Delta t
+
+ ZW2(:,:,:) = (ZW2(:,:,:) / PTSTEP) - PRIS(:,:,:) ! idem ZW1 but for Ri
+
+ ELSE
+!
+!
+!* 4. SECOND ORDER ALL OR NOTHING CONDENSATION SCHEME
+! FOR MIXED-PHASE CLOUD
+! -----------------------------------------------
+!
+!
+!* 4.1 Eventually pretreatment
+!
+ IF (LPRETREATMENT) THEN
+!
+! compute the saturation vapor pressures at t+1
+!
+ CALL GET_HALO(ZT)
+ ZW1(:,:,:) = EXP(XALPW - XBETAW/ZT(:,:,:) - XGAMW*ALOG(ZT(:,:,:))) ! e_sw
+ ZW2(:,:,:) = EXP(XALPI - XBETAI/ZT(:,:,:) - XGAMI*ALOG(ZT(:,:,:))) ! e_si
+ ZW1(:,:,:) = MIN(PPABST(:,:,:)/2.,ZW1(:,:,:)) ! safety limitation
+ ZW2(:,:,:) = MIN(PPABST(:,:,:)/2.,ZW2(:,:,:)) ! safety limitation
+!
+! compute the saturation mixing ratios at t+1
+!
+ ZW3(:,:,:) = ZW1(:,:,:) * ZEPS / &
+ ( PPABST(:,:,:) - ZW1(:,:,:)) ! r_sw
+ ZW4(:,:,:) = ZW2(:,:,:) * ZEPS / &
+ ( PPABST(:,:,:) - ZW2(:,:,:)) ! r_si
+!
+ WHERE(PRVS(:,:,:)*PTSTEP .LT. ZW4(:,:,:) .AND. &
+ PRCS(:,:,:) .GT. 0. .AND. ZT(:,:,:) .LT. XTT)
+!
+! Subsaturation case with respect to rsi(T,P) (and case rv<0):
+! Evaporation of rc>0 (while enough) to decrease the lack of vapor
+!
+ ZW5 (:,:,:)= MIN( PRCS , ZW4(:,:,:)/PTSTEP - PRVS(:,:,:) ) ! RVCNDC
+ PRVS(:,:,:)= PRVS(:,:,:) + ZW5(:,:,:)
+ PRCS(:,:,:)= PRCS(:,:,:) - ZW5(:,:,:)
+ PTHS(:,:,:)= PTHS(:,:,:) - ZW5(:,:,:) * ZLV(:,:,:) /(ZCPH(:,:,:)*PEXNREF(:,:,:))
+!
+ END WHERE
+!
+ WHERE (PRVS(:,:,:)*PTSTEP .GT. ZW3(:,:,:))
+!
+! Supersaturation case with respect to rsw(T,P):
+! Condensation of the vapor that is left
+!
+ ZW5 (:,:,:)= PRVS(:,:,:) - ZW3(:,:,:)/PTSTEP
+ PRVS(:,:,:)= PRVS(:,:,:) - ZW5(:,:,:) ! RVCNDC
+ PRCS(:,:,:)= PRCS(:,:,:) + ZW5(:,:,:)
+ PTHS(:,:,:)= PTHS(:,:,:) + ZW5(:,:,:) * ZLV(:,:,:) /(ZCPH(:,:,:)*PEXNREF(:,:,:))
+!
+ END WHERE
+!
+ WHERE (PRCS(:,:,:) .GT. 0. .AND. ZT(:,:,:) .LT. ZT00)
+!
+! Treatment of rc>0 if T<T00:
+!
+ PRIS(:,:,:)= PRIS(:,:,:) + PRCS(:,:,:)
+ PTHS(:,:,:)= PTHS(:,:,:) + PRCS(:,:,:) * &
+ (ZLS(:,:,:) - ZLV(:,:,:)) / (ZCPH(:,:,:) * PEXNREF(:,:,:))
+ PRCS(:,:,:)= 0.
+!
+ END WHERE
+!
+!* 4.2 compute the intermediate temperature at t+1, T*
+!
+ ZT(:,:,:) = (PTHS(:,:,:) * PTSTEP) * ZEXNS(:,:,:)
+!
+ END IF !end PRETREATMENT
+!
+!* 4.3 compute the saturation vapor pressures at t+1
+!
+ ZW1(:,:,:) = EXP(XALPW - XBETAW / ZT(:,:,:) - XGAMW * ALOG(ZT(:,:,:))) ! e_sw
+ ZW2(:,:,:) = EXP(XALPI - XBETAI / ZT(:,:,:) - XGAMI * ALOG(ZT(:,:,:))) ! e_si
+ ZW1(:,:,:) = MIN(PPABST(:,:,:)/2.,ZW1(:,:,:)) ! safety limitation
+ ZW2(:,:,:) = MIN(PPABST(:,:,:)/2.,ZW2(:,:,:)) ! safety limitation
+!
+!* 4.4 compute the saturation mixing ratios at t+1
+!
+ ZW3(:,:,:) = ZW1(:,:,:) * ZEPS &
+ / ( PPABST(:,:,:) - ZW1(:,:,:) ) ! r_sw
+ ZW4(:,:,:) = ZW2(:,:,:) * ZEPS &
+ / ( PPABST(:,:,:) - ZW2(:,:,:) ) ! r_si
+!
+!* 4.5 compute the saturation mixing ratio derivatives (r'_vs)
+!
+ ZW1(:,:,:) = (( XBETAW/ZT(:,:,:) - XGAMW ) / ZT(:,:,:)) & ! r'_sw
+ * ZW3(:,:,:) * ( 1. + ZW3(:,:,:)/ZEPS )
+ ZW2(:,:,:) = (( XBETAI/ZT(:,:,:) - XGAMI ) / ZT(:,:,:)) & ! r'_si
+ * ZW4(:,:,:) * ( 1. + ZW4(:,:,:)/ZEPS )
+!
+ IF (LNEW_ADJUST) THEN
+ ZCND(:,:,:)= (ZT(:,:,:) - ZT00) / (ZT0 - ZT00) ! Like Tao et al 89
+ ZCND(:,:,:)= MAX ( MIN(ZCND(:,:,:),1.) , 0. )
+ ELSE
+ WHERE ((PRCS(:,:,:)+PRIS(:,:,:)) .GT. 1.0E-20) &
+ ZCND(:,:,:)= PRCS(:,:,:) / (PRCS(:,:,:) + PRIS(:,:,:)) ! Like the original version
+ END IF
+!
+!* 4.5 compute L_v CND + L_s DEP and F'(T)
+!
+ WHERE ((PRCS(:,:,:)+PRIS(:,:,:)) .GT. 1.0E-20)
+!
+ ZW5(:,:,:) = ZLS(:,:,:) + (ZLV(:,:,:) - ZLS(:,:,:)) * ZCND(:,:,:)
+ ZW6(:,:,:) = ZCPH(:,:,:) * (PRCS(:,:,:) + PRIS(:,:,:)) + &
+ ZW5(:,:,:) * (PRCS(:,:,:) * ZW1(:,:,:) &
+ + PRIS(:,:,:) * ZW2(:,:,:))
+!
+!* 4.6 compute Delta 2
+!
+ ZW7(:,:,:) = (ZW5(:,:,:) / (ZW6(:,:,:) * ZT(:,:,:))) * &
+ (PRCS(:,:,:) * ZW1(:,:,:) * &
+ ((-2. * XBETAW + XGAMW * ZT(:,:,:)) / (XBETAW - XGAMW * ZT(:,:,:)) + &
+ (XBETAW - XGAMW * ZT(:,:,:)) * (1.0 + 2.0 * ZW3(:,:,:) / ZEPS) / ZT(:,:,:)) + &
+ PRIS(:,:,:) * ZW2(:,:,:) * &
+ ((-2. * XBETAI + XGAMI * ZT(:,:,:)) / (XBETAI - XGAMI * ZT(:,:,:)) + &
+ (XBETAI - XGAMI * ZT(:,:,:)) * (1.0 + 2.0 * ZW4(:,:,:) / ZEPS) / ZT(:,:,:)))
+!
+!* 4.7 compute Delta 1
+!
+ ZW6(:,:,:) = ZW5(:,:,:) * (PRCS(:,:,:) * ZW3(:,:,:) + PRIS(:,:,:) * ZW4(:,:,:) - &
+ PRVS(:,:,:) * PTSTEP * (PRCS(:,:,:) + PRIS(:,:,:))) / &
+ ZW6(:,:,:)
+!
+!* 4.8 compute the sources
+!
+ ZW3(:,:,:) = (ZCPH(:,:,:) / ZW5(:,:,:)) * &
+ (-ZW6(:,:,:) * (1.0 + 0.5 * ZW6(:,:,:) * ZW7(:,:,:))) / PTSTEP
+ ZW1(:,:,:) = ZW3(:,:,:) * ZCND(:,:,:) ! RVCNDC
+ ZW2(:,:,:) = ZW3(:,:,:) * (1.0 - ZCND(:,:,:)) ! RVDEPI
+!
+ ELSEWHERE
+!
+!* 4.9 special case when both r_c and r_i are zero
+!
+ ZW6(:,:,:) = ZCPH(:,:,:) + ZLV(:,:,:) * ZW1(:,:,:) ! F'(T)
+ ZW7(:,:,:) = (ZLV(:,:,:) / (ZW6(:,:,:) * ZT(:,:,:))) * & ! Delta 2
+ (ZW1(:,:,:) * &
+ ((-2. * XBETAW + XGAMW * ZT(:,:,:)) / (XBETAW - XGAMW * ZT(:,:,:)) + &
+ (XBETAW - XGAMW * ZT(:,:,:)) * (1.0 + 2.0 * ZW3(:,:,:) / ZEPS) / ZT(:,:,:)))
+ ZW6(:,:,:) = ZLV(:,:,:) * (ZW3(:,:,:) - PRVS(:,:,:) * PTSTEP) / ZW6(:,:,:) ! Delta 1
+ ZW1(:,:,:) = (ZCPH(:,:,:) / ZLV(:,:,:)) * & ! RVCNDC
+ (-ZW6(:,:,:) * ( 1.0 + 0.5 * ZW6(:,:,:) * ZW7(:,:,:))) / PTSTEP
+ ZW2(:,:,:) = 0.0 ! RVDEPI
+!
+ END WHERE
+ END IF
+!
+!* 5. COMPUTE THE SOURCES AND STORES THE CLOUD FRACTION
+! -------------------------------------------------
+!
+!* 5.1 compute the sources
+!
+!* 5.1.1 microphysics
+!
+ WHERE (ZW1(:,:,:) < 0.0)
+ ZW1(:,:,:) = MAX (ZW1(:,:,:), -PRCS(:,:,:)) ! Evaporation rate
+ ELSEWHERE
+ ZW1(:,:,:) = MIN (ZW1(:,:,:), PRVS(:,:,:)) ! Condensation rate
+ END WHERE
+!
+ PRVS(:,:,:) = PRVS(:,:,:) - ZW1(:,:,:)
+ PRCS(:,:,:) = PRCS(:,:,:) + ZW1(:,:,:)
+ PTHS(:,:,:) = PTHS(:,:,:) + &
+ ZW1(:,:,:) * ZLV(:,:,:) / (ZCPH(:,:,:) * PEXNREF(:,:,:))
+!
+ WHERE (ZW2(:,:,:) < 0.0)
+ ZW2(:,:,:) = MAX (ZW2(:,:,:), -PRIS(:,:,:)) ! Sublimation rate
+ ELSEWHERE
+ ZW2(:,:,:) = MIN (ZW2(:,:,:), PRVS(:,:,:)) ! Deposition rate
+ END WHERE
+!
+ PRVS(:,:,:) = PRVS(:,:,:) - ZW2(:,:,:)
+ PRIS(:,:,:) = PRIS(:,:,:) + ZW2(:,:,:)
+ PTHS(:,:,:) = PTHS(:,:,:) + &
+ ZW2(:,:,:) * ZLS(:,:,:) / (ZCPH(:,:,:) * PEXNREF(:,:,:))
+!
+!* 5.1.2 electricity
+!
+ ZWE1(:,:,:) = 0.
+ ZWE2(:,:,:) = 0.
+!
+! the electrical process due to condensation is removed and
+! capture of ions by cloud droplets is done in ion_attach_elec routine
+!
+! evaporation
+ WHERE (ABS(PRCT(:,:,:)) > XRTMIN_ELEC(2) .AND. &
+ ABS(PQCT(:,:,:)) > XQTMIN(2) .AND. &
+ ZW1(:,:,:) < -XRTMIN(1))
+ ZWE1(:,:,:) = (XFC / 3.) * (PQCT(:,:,:) / PRCT(:,:,:)) * (-ZW1(:,:,:))
+ ZION_NUMBER(:,:,:) = ABS(ZWE1(:,:,:)) / XECHARGE
+ ZADD(:,:,:) = 0.5 + SIGN(0.5, ZWE1(:,:,:))
+ PQPIS(:,:,:) = PQPIS(:,:,:) + &
+ ZADD(:,:,:) *ZION_NUMBER(:,:,:)
+ PQNIS(:,:,:) = PQNIS(:,:,:) + &
+ (1.-ZADD(:,:,:)) *ZION_NUMBER(:,:,:)
+ PQCS(:,:,:) = PQCS(:,:,:) - ZWE1(:,:,:)
+ END WHERE
+!
+! the electrical process due to deposition is removed and
+! capture of ions by raindropsp is done in ion_attach_elec routine
+!
+! sublimation
+ WHERE (ABS(PRIT(:,:,:)) > XRTMIN_ELEC(4) .AND. &
+ ABS(PQIT(:,:,:)) > XQTMIN(4) .AND. &
+ ZW2(:,:,:) < -XRTMIN(1))
+ ZWE2(:,:,:) = (XFI / XBI) * (PQIT(:,:,:) / PRIT(:,:,:)) * (-ZW2(:,:,:))
+ ZION_NUMBER(:,:,:) = ABS(ZWE2(:,:,:)) / XECHARGE
+ ZADD(:,:,:) = 0.5 + SIGN(0.5, ZWE2(:,:,:))
+ PQPIS(:,:,:) = PQPIS(:,:,:) + &
+ ZADD(:,:,:) *ZION_NUMBER(:,:,:)
+ PQNIS(:,:,:) = PQNIS(:,:,:) + &
+ (1.-ZADD(:,:,:)) *ZION_NUMBER(:,:,:)
+ PQIS(:,:,:) = PQIS(:,:,:) - ZWE2(:,:,:)
+ END WHERE
+END DO ! end of the iterative loop
+!
+!
+!* 5.2 compute the cloud fraction PCLDFR
+!
+IF (.NOT. OSUBG_COND) THEN
+ WHERE (PRCS(:,:,:) + PRIS(:,:,:) > 1.E-12 / PTSTEP)
+ PCLDFR(:,:,:) = 1.
+ ELSEWHERE
+ PCLDFR(:,:,:) = 0.
+ ENDWHERE
+ IF (SIZE(PSRCS,3) /= 0) THEN
+ PSRCS(:,:,:) = PCLDFR(:,:,:)
+ END IF
+ELSE
+ IF (HSCONV == 'EDKF' .AND. HMF_CLOUD == 'DIRE') THEN
+ PCLDFR(:,:,:) = MIN(1.,PCLDFR(:,:,:)+PCF_MF(:,:,:))
+ PRCS(:,:,:) = PRCS(:,:,:) + PRC_MF(:,:,:) / PTSTEP
+ PRIS(:,:,:) = PRIS(:,:,:)+PRI_MF(:,:,:)/PTSTEP
+ PRVS(:,:,:) = PRVS(:,:,:)- ( PRC_MF(:,:,:) + PRI_MF(:,:,:)) /PTSTEP
+ PTHS(:,:,:) = PTHS(:,:,:) + ( PRC_MF(:,:,:) * ZLV(:,:,:) + &
+ PRI_MF(:,:,:) * ZLS(:,:,:) ) / ZCPH(:,:,:) / &
+ PEXNREF(:,:,:) / PTSTEP
+ END IF
+ENDIF
+!
+!
+!
+!* 6. STORE THE BUDGET TERMS
+! ----------------------
+!
+if ( lbudget_th ) call Budget_store_end( tbudgets(NBUDGET_TH), 'DEPI', pths(:, :, :) * prhodj(:, :, :) )
+if ( lbudget_rv ) call Budget_store_end( tbudgets(NBUDGET_RV), 'DEPI', prvs(:, :, :) * prhodj(:, :, :) )
+if ( lbudget_rc ) call Budget_store_end( tbudgets(NBUDGET_RC), 'DEPI', prcs(:, :, :) * prhodj(:, :, :) )
+if ( lbudget_ri ) call Budget_store_end( tbudgets(NBUDGET_RI), 'DEPI', pris(:, :, :) * prhodj(:, :, :) )
+if ( lbudget_sv ) then
+ call Budget_store_end( tbudgets(NBUDGET_SV1 - 1 + nsv_elecbeg ), 'DEPI', pqpis(:, :, :) * prhodj(:, :, :) )
+ call Budget_store_end( tbudgets(NBUDGET_SV1 - 1 + nsv_elecend ), 'DEPI', pqnis(:, :, :) * prhodj(:, :, :) )
+ call Budget_store_end( tbudgets(NBUDGET_SV1 - 1 + nsv_elecbeg + 1 ), 'DEPI', pqcs (:, :, :) * prhodj(:, :, :) )
+ call Budget_store_end( tbudgets(NBUDGET_SV1 - 1 + nsv_elecbeg + 3 ), 'DEPI', pqis (:, :, :) * prhodj(:, :, :) )
+end if
+!------------------------------------------------------------------------------
+!
+END SUBROUTINE ICE_ADJUST_ELEC
diff --git a/src/mesonh/micro/lima_adjust_split.f90 b/src/mesonh/micro/lima_adjust_split.f90
index edaeec82007c9f14ed51c77c8aa092ad3bbc1d7b..ab87f141f0ee8d1dc4cb3bae80bf3d81b4259f8c 100644
--- a/src/mesonh/micro/lima_adjust_split.f90
+++ b/src/mesonh/micro/lima_adjust_split.f90
@@ -269,11 +269,12 @@ REAL, DIMENSION(SIZE(PRHODJ,1),SIZE(PRHODJ,2),SIZE(PRHODJ,3)) &
ZLV, & ! guess of the Lv at t+1
ZLS, & ! guess of the Ls at t+1
ZMASK,&
- ZRV, ZRV2, &
- ZRC, ZRC2, &
- ZRI, &
+ ZRV, ZRV2,ZRV_IN, &
+ ZRC, ZRC2,ZRC_IN, &
+ ZRI, ZRI_IN, &
ZSIGS, &
ZW_MF
+REAL, DIMENSION(SIZE(PRHODJ,1),SIZE(PRHODJ,2)) :: ZSIGQSAT2D
LOGICAL, DIMENSION(SIZE(PRHODJ,1),SIZE(PRHODJ,2),SIZE(PRHODJ,3)) &
:: GMICRO ! Test where to compute cond/dep proc.
INTEGER :: IMICRO
@@ -500,15 +501,20 @@ DO JITER =1,ITERMAX
!
ZRV=PRVS*PTSTEP
ZRC=PRCS*PTSTEP
+ ZRV_IN=ZRV
+ ZRC_IN=ZRC
+ ZRI_IN=0.
ZRV2=PRVT
ZRC2=PRCT
ZRI=0.
ZSIGS=PSIGS
+ ZSIGQSAT2D(:,:)=PSIGQSAT
CALL CONDENSATION(IIU, IJU, IKU, IIB, IIE, IJB, IJE, IKB, IKE, 1, 'S', &
HCONDENS, HLAMBDA3, &
- PPABST, PZZ, PRHODREF, ZT, ZRV, ZRC, ZRI, PRSS*PTSTEP, PRGS*PTSTEP, &
- ZSIGS, PMFCONV, PCLDFR, PSRCS, .FALSE., OSIGMAS, &
- PSIGQSAT, PLV=ZLV, PLS=ZLS, PCPH=ZCPH )
+ PPABST, PZZ, PRHODREF, ZT, ZRV_IN, ZRV, ZRC_IN, ZRC, ZRI_IN, ZRI,&
+ PRSS*PTSTEP, PRGS*PTSTEP, &
+ ZSIGS, PMFCONV, PCLDFR, PSRCS, .FALSE., OSIGMAS, .FALSE., &
+ ZSIGQSAT2D, PLV=ZLV, PLS=ZLS, PCPH=ZCPH )
PCLDFR(:,:,:) = MIN(PCLDFR(:,:,:) + PCF_MF(:,:,:) , 1.)
ZRV(:,:,:) = ZRV(:,:,:) - MAX(MIN(PRC_MF(:,:,:), ZRV(:,:,:)),0.)
ZRC(:,:,:) = ZRC(:,:,:) + MAX(MIN(PRC_MF(:,:,:), ZRV(:,:,:)),0.)
diff --git a/src/mesonh/micro/radtr_satel.f90 b/src/mesonh/micro/radtr_satel.f90
new file mode 100644
index 0000000000000000000000000000000000000000..ccad5e716ca659cfd4a058ce37ea5b38b06ad531
--- /dev/null
+++ b/src/mesonh/micro/radtr_satel.f90
@@ -0,0 +1,738 @@
+!MNH_LIC Copyright 2000-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.
+!-----------------------------------------------------------------
+! #######################
+ MODULE MODI_RADTR_SATEL
+! #######################
+INTERFACE
+!
+ SUBROUTINE RADTR_SATEL(KYEARF, KMONTHF, KDAYF, PSECF, &
+ KDLON, KFLEV, KSTATM, KRAD_COLNBR, PEMIS, PCCO2, &
+ PTSRAD, PSTATM, PTHT, PRT, PPABST, PZZ, &
+ PSIGS, PMFCONV, PCLDFR, OUSERI, OSIGMAS, &
+ OSUBG_COND, ORAD_SUBG_COND, PIRBT, PWVBT, KGEO,PSIGQSAT )
+!
+INTEGER, INTENT(IN) :: KYEARF ! year of Final date
+INTEGER, INTENT(IN) :: KMONTHF ! month of Final date
+INTEGER, INTENT(IN) :: KDAYF ! day of Final date
+REAL, INTENT(IN) :: PSECF ! number of seconds since date at 00 UTC
+!
+INTEGER, INTENT(IN) :: KDLON !number of columns where the
+ !radiation calculations are performed
+INTEGER, INTENT(IN) :: KFLEV !number of vertical levels where the
+ !radiation calculations are performed
+INTEGER, INTENT(IN) :: KSTATM !index of the standard atmosphere level
+ !just above the model top
+INTEGER, INTENT(IN) :: KRAD_COLNBR !factor by which the memory is split
+!
+REAL, DIMENSION(:,:), INTENT(IN) :: PEMIS !Surface IR EMISsivity
+REAL, INTENT(IN) :: PCCO2 !CO2 content
+REAL, DIMENSION(:,:), INTENT(IN) :: PTSRAD !RADiative Surface Temperature
+REAL, DIMENSION(:,:), INTENT(IN) :: PSTATM !selected standard atmosphere
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHT !THeta at t
+REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PRT !moist variables at t
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PPABST !pressure at t
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PZZ !Model level heights
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PSIGS ! Sigma_s from turbulence scheme
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PMFCONV! convective mass flux (kg /s m^2)
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PCLDFR ! cloud fraction
+!
+LOGICAL, INTENT(IN) :: OUSERI ! logical switch to compute both
+ ! liquid and solid condensate (OUSERI=.TRUE.)
+ ! or only liquid condensate (OUSERI=.FALSE.)
+LOGICAL, INTENT(IN) :: OSIGMAS! use present global Sigma_s values
+ ! or that from turbulence scheme
+LOGICAL, INTENT(IN) :: OSUBG_COND ! Switch for Subgrid Condensation
+ ! (prognotic mode)
+LOGICAL, INTENT(IN) :: ORAD_SUBG_COND ! Switch for Subgrid Condensation
+ ! (diagnostic mode)
+!
+REAL, DIMENSION(:,:), INTENT(OUT):: PIRBT !IR Brightness Temp. (K)
+REAL, DIMENSION(:,:), INTENT(OUT):: PWVBT !WV Brightness Temp. (K)
+!
+INTEGER, INTENT(IN) :: KGEO !SATELLITE INDEX
+REAL, INTENT(IN) :: PSIGQSAT ! use an extra "qsat" variance contribution (OSIGMAS case)
+!
+END SUBROUTINE RADTR_SATEL
+END INTERFACE
+END MODULE MODI_RADTR_SATEL
+! #####################################################################
+ SUBROUTINE RADTR_SATEL(KYEARF, KMONTHF, KDAYF, PSECF, &
+ KDLON, KFLEV, KSTATM, KRAD_COLNBR, PEMIS, PCCO2, &
+ PTSRAD, PSTATM, PTHT, PRT, PPABST, PZZ, &
+ PSIGS, PMFCONV, PCLDFR, OUSERI, OSIGMAS, &
+ OSUBG_COND, ORAD_SUBG_COND, PIRBT, PWVBT, KGEO,PSIGQSAT)
+! #####################################################################
+!
+!!**** *RADTR_SATEL* -
+!!
+!! PURPOSE
+!! -------
+!!
+!!** METHOD
+!! ------
+!!
+!! EXTERNAL
+!! --------
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!!
+!! REFERENCE
+!! ---------
+!! Chaboureau, J.-P., J.-P. Cammas, P. Mascart, J.-P. Pinty, C. Claud, R. Roca,
+!! and J.-J. Morcrette, 2000: Evaluation of a cloud system life-cycle simulated
+!! by Meso-NH during FASTEX using METEOSAT radiances and TOVS-3I cloud retrievals.
+!! Q. J. R. Meteorol. Soc., 126, 1735-1750.
+!! Chaboureau, J.-P. and P. Bechtold, 2002: A simple cloud parameterization from
+!! cloud resolving model data: Theory and application. J. Atmos. Sci., 59, 2362-2372.
+!!
+!! AUTHOR
+!! ------
+!! J.-P. Chaboureau *L.A.*
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 29/03/00
+!! J.-P. Chaboureau 15/04/03 add call to the subgrid condensation scheme
+!! J.Escobar : 15/09/2015 : WENO5 & JPHEXT <> 1
+!! G.Delautier 04/2016 : BUG JPHEXT
+!! S. Riette 11/2016 : Condensation interface changed
+! P. Wautelet 26/04/2019: replace non-standard FLOAT function by REAL function
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODD_CST
+USE MODD_PARAMETERS
+USE MODD_GRID_n
+!
+USE MODD_RAD_TRANSF
+USE MODE_ll
+!
+USE MODI_INIT_NBMOD
+USE MODI_DETER_ANGLE
+USE MODI_MAKE_RADSAT
+!
+USE MODI_CONDENSATION
+!
+IMPLICIT NONE
+!
+!* 0.1 DECLARATIONS OF DUMMY ARGUMENTS :
+!
+INTEGER, INTENT(IN) :: KYEARF ! year of Final date
+INTEGER, INTENT(IN) :: KMONTHF ! month of Final date
+INTEGER, INTENT(IN) :: KDAYF ! day of Final date
+REAL, INTENT(IN) :: PSECF ! number of seconds since date at 00 UTC
+!
+INTEGER, INTENT(IN) :: KDLON !number of columns where the
+ ! radiation calculations are performed
+INTEGER, INTENT(IN) :: KFLEV !number of vertical levels where the
+ ! radiation calculations are performed
+INTEGER, INTENT(IN) :: KSTATM !index of the standard atmosphere level
+ !just above the model top
+INTEGER, INTENT(IN) :: KRAD_COLNBR !factor by which the memory is split
+!
+REAL, DIMENSION(:,:), INTENT(IN) :: PEMIS !Surface IR EMISsivity
+REAL, INTENT(IN) :: PCCO2 !CO2 content
+REAL, DIMENSION(:,:), INTENT(IN) :: PTSRAD !RADiative Surface Temperature
+REAL, DIMENSION(:,:), INTENT(IN) :: PSTATM !selected standard atmosphere
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHT !THeta at t
+REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PRT !moist variables at t
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PPABST !pressure at t
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PZZ !Model level heights
+!
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PSIGS ! Sigma_s from turbulence scheme
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PMFCONV! convective mass flux (kg /s m^2)
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PCLDFR ! cloud fraction
+!
+LOGICAL, INTENT(IN) :: OUSERI ! logical switch to compute both
+ ! liquid and solid condensate (OUSERI=.TRUE.)
+ ! or only liquid condensate (OUSERI=.FALSE.)
+LOGICAL, INTENT(IN) :: OSIGMAS! use present global Sigma_s values
+ ! or that from turbulence scheme
+LOGICAL, INTENT(IN) :: OSUBG_COND ! Switch for Subgrid Condensation
+ ! (prognotic mode)
+LOGICAL, INTENT(IN) :: ORAD_SUBG_COND ! Switch for Subgrid Condensation
+ ! (diagnostic mode)
+!
+REAL, DIMENSION(:,:), INTENT(OUT):: PIRBT !IR Brightness Temp. (K)
+REAL, DIMENSION(:,:), INTENT(OUT):: PWVBT !WV Brightness Temp. (K)
+!
+INTEGER, INTENT(IN) :: KGEO !SATELLITE INDEX
+REAL, INTENT(IN) :: PSIGQSAT ! use an extra "qsat" variance contribution (OSIGMAS case)
+!
+!* 0.2 DECLARATIONS OF LOCAL VARIABLES
+!
+LOGICAL :: GPTDEP, GPVOIGT
+!
+! reference state
+!from inprof
+INTEGER :: IGL, ICABS, ING1, IUABS, IINIS, IENDS, ICONF, ICLOUD, IOVLP
+INTEGER :: IH2O, ICO2, IO3, ICNT, IN2O, ICH4, ICO, IC11, IC12, ICFC
+!
+LOGICAL, DIMENSION(KDLON) :: GDOIT_2D ! .TRUE. for the larger scale
+LOGICAL, DIMENSION(KDLON,KFLEV) :: GDOIT ! .TRUE. for all the levels of the
+ ! larger scale columns
+!
+INTEGER :: JI,JJ,JK,JK1,JK2,JKRAD ! loop indexes
+!
+INTEGER :: IIB,IIE ! I index value of the first/last inner mass point
+INTEGER :: IJB,IJE ! J index value of the first/last inner mass point
+INTEGER :: IKB,IKE ! K index value of the first/last inner mass point
+INTEGER :: IIU ! array size for the first index
+INTEGER :: IJU ! array size for the second index
+INTEGER :: IKU ! array size for the third index
+INTEGER :: IIJ ! reformatted array index
+INTEGER :: IKSTAE ! level number of the STAndard atmosphere array
+INTEGER :: IKUP ! vertical level above which STAndard atmosphere data
+INTEGER :: IDOIT_COL ! number of larger scale columns
+INTEGER :: IDOIT ! number of levels corresponding of the larger scale
+ ! columns are filled in
+INTEGER :: IDIM ! effective number of columns for which the radiation
+ ! code is run
+INTEGER, DIMENSION(SIZE(PZZ,1),SIZE(PZZ,3)) :: IKKOZ ! indice array used to
+ ! vertically interpolate the ozone content on the model grid
+!
+REAL, DIMENSION(:,:), ALLOCATABLE :: ZTAVE ! mean-layer temperature
+REAL, DIMENSION(:,:), ALLOCATABLE :: ZQVAVE ! mean-layer specific humidity
+REAL, DIMENSION(:,:), ALLOCATABLE :: ZO3AVE ! mean-layer ozone content
+REAL, DIMENSION(:,:), ALLOCATABLE :: ZPRES_HL ! half-level pressure
+REAL, DIMENSION(:,:), ALLOCATABLE :: ZT_HL ! half-level temperature
+REAL, DIMENSION(:,:), ALLOCATABLE :: ZCLDLD ! Downward cloud emissivity
+REAL, DIMENSION(:,:), ALLOCATABLE :: ZCLDLU ! Upward cloud emissivity
+REAL, DIMENSION(:), ALLOCATABLE :: ZVIEW ! cosecant of viewing angle
+REAL, DIMENSION(:), ALLOCATABLE :: ZREMIS ! Reformatted PEMIS array
+REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZEXNT ! Exner function
+REAL, DIMENSION(SIZE(PSTATM,1)) :: ZSTAZZ,ZSTAOZ ! STAndard atmosphere height
+ ! and OZone content
+REAL :: ZOZ ! variable used to interpolate the ozone profile
+!
+REAL, DIMENSION(:), ALLOCATABLE :: ZDT0 ! surface discontinuity
+REAL, DIMENSION(:,:), ALLOCATABLE :: ZRADBT
+REAL, DIMENSION(:,:), ALLOCATABLE :: ZRADBC
+REAL, DIMENSION(:,:), ALLOCATABLE :: ZRADFT
+REAL, DIMENSION(:), ALLOCATABLE :: ZULAT
+REAL, DIMENSION(:), ALLOCATABLE :: ZULON
+!
+REAL, DIMENSION(:,:), ALLOCATABLE :: ZZRADFT
+!
+REAL, DIMENSION(:), ALLOCATABLE :: ZWORK1, ZWORK3
+!
+! split arrays used to split the memory required by the ECMWF_radiation
+! subroutine, the fields have the same meaning as their complete counterpart
+REAL, DIMENSION(:), ALLOCATABLE :: ZREMIS_SPLIT
+REAL, DIMENSION(:,:), ALLOCATABLE :: ZO3AVE_SPLIT
+REAL, DIMENSION(:,:), ALLOCATABLE :: ZT_HL_SPLIT
+REAL, DIMENSION(:,:), ALLOCATABLE :: ZPRES_HL_SPLIT
+REAL, DIMENSION(:,:), ALLOCATABLE :: ZQVAVE_SPLIT
+REAL, DIMENSION(:,:), ALLOCATABLE :: ZTAVE_SPLIT
+REAL, DIMENSION(:,:), ALLOCATABLE :: ZCLDLD_SPLIT
+REAL, DIMENSION(:,:), ALLOCATABLE :: ZCLDLU_SPLIT
+REAL, DIMENSION(:), ALLOCATABLE :: ZVIEW_SPLIT
+REAL, DIMENSION(:), ALLOCATABLE :: ZDT0_SPLIT
+REAL, DIMENSION(:,:), ALLOCATABLE :: ZRADBT_SPLIT
+REAL, DIMENSION(:,:), ALLOCATABLE :: ZRADBC_SPLIT
+!
+INTEGER :: JI_SPLIT ! loop on the split array
+INTEGER :: INUM_CALL ! number of CALL of the radiation scheme
+INTEGER :: IDIM_EFF ! effective number of air-columns to compute
+INTEGER :: IDIM_RESIDUE ! number of remaining air-columns to compute
+INTEGER :: IBEG, IEND ! auxiliary indices
+!
+! Other arrays for emissivity
+REAL :: ZFLWP, ZFIWP, ZANGCOR, ZRADLP, ZMULTS, ZTMP, ZKI
+!
+! Other arrays for condensation
+REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZTEMP ! Temperature
+REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZSIGRC ! s r_c / sig_s^2
+REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZNCLD ! grid scale cloud fraction
+REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZRC_IN, ZRC_OUT ! grid scale r_c mixing ratio (kg/kg)
+REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZRI_IN, ZRI_OUT ! grid scale r_i (kg/kg)
+REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZRV_IN, ZRV_OUT ! grid scale r_v (kg/kg)
+REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZRHO
+REAL, DIMENSION(SIZE(PPABST,1),SIZE(PPABST,2)) :: ZSIGQSAT2D
+!----------------------------------------------------------------------------
+!
+!* 1. INITIALIZATION OF CONSTANTS FOR TRANSFERT CODE
+! ----------------------------------------------
+!
+CALL INIT_NBMOD(KFLEV, IGL, ICABS, ING1, IUABS, IINIS, IENDS, &
+ IH2O, ICO2, IO3, ICNT, IN2O, ICH4, ICO, IC11, IC12, ICFC, &
+ ICONF, ICLOUD, IOVLP, GPVOIGT, GPTDEP)
+X1CO2 = PCCO2 / 44.0 * XMD
+!
+!----------------------------------------------------------------------------
+!
+!* 2. COMPUTE DIMENSIONS OF ARRAYS AND OTHER INDICES
+! ----------------------------------------------
+!
+IIU = SIZE(PTHT,1)
+IJU = SIZE(PTHT,2)
+IKU = SIZE(PTHT,3)
+CALL GET_INDICE_ll (IIB,IJB,IIE,IJE)
+IKB = 1 + JPVEXT
+IKE = IKU - JPVEXT
+!
+IKSTAE = SIZE(PSTATM,1)
+IKUP = IKE-JPVEXT+1
+!
+!----------------------------------------------------------------------------
+!
+!* 3. INITIALIZES THE MEAN-LAYER VARIABLES
+! ------------------------------------
+!
+ALLOCATE(ZEXNT(SIZE(PTHT,1),SIZE(PTHT,2),SIZE(PTHT,3)))
+ZEXNT(:,:,:)= ( PPABST(:,:,:)/XP00 ) ** (XRD/XCPD)
+!
+ALLOCATE(ZTAVE(KDLON,KFLEV))
+ALLOCATE(ZQVAVE(KDLON,KFLEV))
+!
+ZQVAVE(:,:) = 0.0
+!
+DO JK=IKB,IKE
+ JKRAD = JK-JPVEXT
+ DO JJ=IJB,IJE
+ DO JI=IIB,IIE
+ IIJ = 1 + (JI-IIB) + (IIE-IIB+1)*(JJ-IJB)
+ ZTAVE(IIJ,JKRAD) = PTHT(JI,JJ,JK)*ZEXNT(JI,JJ,JK)
+ END DO
+ END DO
+END DO
+!
+! Check if the humidity mixing ratio is available
+!
+IF( SIZE(PRT(:,:,:,:),4) >= 1 ) THEN
+ DO JK=IKB,IKE
+ JKRAD = JK-JPVEXT
+ DO JJ=IJB,IJE
+ DO JI=IIB,IIE
+ IIJ = 1 + (JI-IIB) + (IIE-IIB+1)*(JJ-IJB)
+ ZQVAVE(IIJ,JKRAD) = PRT(JI,JJ,JK,1)
+ END DO
+ END DO
+ END DO
+END IF
+!
+! Standard atmosphere extension
+!
+DO JK=IKUP,KFLEV
+ JK1 = (KSTATM-1)+(JK-IKUP)
+ JK2 = JK1+1
+ ZTAVE(:,JK) = 0.5*( PSTATM(JK1,3)+PSTATM(JK2,3) )
+ ZQVAVE(:,JK) = 0.5*( PSTATM(JK1,5)/PSTATM(JK1,4)+ &
+ PSTATM(JK2,5)/PSTATM(JK2,4) )
+END DO
+!
+!----------------------------------------------------------------------------
+!
+!* 4. INITIALIZES THE HALF-LEVEL VARIABLES
+! ------------------------------------
+!
+ALLOCATE(ZPRES_HL(KDLON,KFLEV+1))
+ALLOCATE(ZT_HL(KDLON,KFLEV+1))
+!
+DO JK=IKB,IKE+1
+ JKRAD = JK-JPVEXT
+ DO JJ=IJB,IJE
+ DO JI=IIB,IIE
+ IIJ = 1 + (JI-IIB) + (IIE-IIB+1)*(JJ-IJB)
+ ZPRES_HL(IIJ,JKRAD) = XP00 * &
+ (0.5*(ZEXNT(JI,JJ,JK)+ZEXNT(JI,JJ,JK-1)))**(XCPD/XRD)
+ END DO
+ END DO
+END DO
+!
+! Standard atmosphere extension
+! begining at ikup+1 level allows to use a model domain higher than 50km
+!
+DO JK=IKUP+1,KFLEV+1
+ JK1 = (KSTATM-1)+(JK-IKUP)
+ ZPRES_HL(:,JK) = PSTATM(JK1,2)*100.0
+END DO
+!
+! Surface temperature at the first level
+DO JJ=IJB,IJE
+ DO JI=IIB,IIE
+ IIJ = 1 + (JI-IIB) + (IIE-IIB+1)*(JJ-IJB)
+ ZT_HL(IIJ,1) = PTSRAD(JI,JJ)
+ END DO
+END DO
+!
+! Temperature at half levels
+ZT_HL(:,2:IKE-JPVEXT) = 0.5*(ZTAVE(:,1:IKE-JPVEXT-1)+ZTAVE(:,2:IKE-JPVEXT))
+!
+DO JJ=IJB,IJE
+ DO JI=IIB,IIE
+ IIJ = 1 + (JI-IIB) + (IIE-IIB+1)*(JJ-IJB)
+ ZT_HL(IIJ,IKE-JPVEXT+1) = 0.5*PTHT(JI,JJ,IKE )*ZEXNT(JI,JJ,IKE ) &
+ + 0.5*PTHT(JI,JJ,IKE+1)*ZEXNT(JI,JJ,IKE+1)
+ END DO
+END DO
+!
+! Standard atmosphere extension
+! begining at ikup+1 level allows to use a model domain higher than 50km
+!
+DO JK=IKUP+1,KFLEV+1
+ JK1 = (KSTATM-1)+(JK-IKUP)
+ ZT_HL(:,JK) = PSTATM(JK1,3)
+END DO
+!
+!----------------------------------------------------------------------------
+!
+!* 5. INITIALIZES THE OZONE PROFILES from the standard atmosphere
+! ------------------------------
+!
+ALLOCATE(ZO3AVE(KDLON,KFLEV))
+!
+ZSTAOZ(:) = PSTATM(:,6)/PSTATM(:,4)
+ZSTAZZ(:) = 1000.0*PSTATM(:,1)
+!
+DO JJ = IJB,IJE
+ DO JK2 = IKB,IKE
+ JKRAD = JK2-JPVEXT
+ IKKOZ(:,JK2) = IKB-1
+ DO JK1 = 1,IKSTAE
+ DO JI = IIB,IIE
+ IKKOZ(JI,JK2)=IKKOZ(JI,JK2) + NINT(0.5 + SIGN(0.5, &
+ -ZSTAZZ(JK1)+0.5*(PZZ(JI,JJ,JK2)+PZZ(JI,JJ,JK2+1)) ))
+ END DO
+ END DO
+ DO JI = IIB,IIE
+ ZOZ=(0.5*(PZZ(JI,JJ,JK2)+PZZ(JI,JJ,JK2+1))- ZSTAZZ(IKKOZ(JI,JK2))) &
+ /( ZSTAZZ(IKKOZ(JI,JK2)+1) - ZSTAZZ(IKKOZ(JI,JK2)))
+ IIJ = 1 + (JI-IIB) + (IIE-IIB+1)*(JJ-IJB)
+ ZO3AVE(IIJ,JKRAD) =( (1.- ZOZ) * ZSTAOZ(IKKOZ(JI,JK2)) &
+ + ZOZ * ZSTAOZ(IKKOZ(JI,JK2)+1))
+ END DO
+ END DO
+END DO
+!
+DO JK=IKUP,KFLEV
+ JK1 = (KSTATM)+(JK-IKUP)
+ ZO3AVE(:,JK) = ZSTAOZ(JK1)
+END DO
+!
+!----------------------------------------------------------------------------
+!
+!* 6. CALLS THE E.C.M.W.F. RADIATION CODE
+! -----------------------------------
+!
+!* 6.1 INITIALIZES 2D AND SURFACE FIELDS
+!
+ALLOCATE(ZREMIS(KDLON))
+DO JJ=IJB,IJE
+ DO JI=IIB,IIE
+ IIJ = 1 + (JI-IIB) + (IIE-IIB+1)*(JJ-IJB)
+ ZREMIS(IIJ) = PEMIS(JI,JJ)
+ END DO
+END DO
+!
+! initializes surface discontinuity field
+ALLOCATE(ZDT0(KDLON))
+DO JJ=IJB,IJE
+ DO JI=IIB,IIE
+ IIJ = 1 + (JI-IIB) + (IIE-IIB+1)*(JJ-IJB)
+ ZDT0(IIJ) = PTSRAD(JI,JJ) - PTHT(JI,JJ,1)*ZEXNT(JI,JJ,1)
+ END DO
+END DO
+!
+ALLOCATE(ZULAT(KDLON))
+ALLOCATE(ZULON(KDLON))
+DO JJ=IJB,IJE
+ DO JI=IIB,IIE
+ IIJ = 1 + (JI-IIB) + (IIE-IIB+1)*(JJ-IJB)
+ ZULON(IIJ) = XLON(JI,JJ)
+ ZULAT(IIJ) = XLAT(JI,JJ)
+ END DO
+END DO
+ALLOCATE(ZVIEW(KDLON))
+CALL DETER_ANGLE(KGEO, KDLON, ZULAT, ZULON, ZVIEW)
+DEALLOCATE(ZULAT)
+DEALLOCATE(ZULON)
+!
+!
+ALLOCATE(ZCLDLD(KDLON,KFLEV))
+ALLOCATE(ZCLDLU(KDLON,KFLEV))
+ZCLDLD = 0.
+ZCLDLU = 0.
+!
+IF( SIZE(PRT(:,:,:,:),4) >= 2 ) THEN
+ ALLOCATE(ZNCLD(IIU,IJU,IKU))
+ ALLOCATE(ZRC_IN(IIU,IJU,IKU))
+ ALLOCATE(ZRC_OUT(IIU,IJU,IKU))
+ ZRC_IN=PRT(:,:,:,2)
+ ALLOCATE(ZRI_IN(IIU,IJU,IKU))
+ ALLOCATE(ZRI_OUT(IIU,IJU,IKU))
+ ZRI_IN=0.
+ IF( OUSERI ) ZRI_IN=PRT(:,:,:,4)
+ IF ( .NOT. OSUBG_COND .AND. ORAD_SUBG_COND) THEN
+ PRINT*,' THE SUBGRID CONDENSATION SCHEME IN DIAGNOSTIC MODE IS ACTIVATED'
+ ALLOCATE(ZTEMP(IIU,IJU,IKU))
+ ZTEMP=PTHT*ZEXNT
+ ALLOCATE(ZSIGRC(IIU,IJU,IKU))
+ ALLOCATE(ZRV_IN(IIU,IJU,IKU))
+
+ ZRV_IN=PRT(:,:,:,1)
+ ALLOCATE(ZRHO(IIU,IJU,IKU))
+ ZRHO=1. !unused
+ ZSIGQSAT2D(:,:)=PSIGQSAT
+ CALL CONDENSATION( IIU, IJU, IKU, IIB, IIE, IJB, IJE, IKB, IKE, 1, 'T', 'CB02', 'CB',&
+ PPABST, PZZ, ZRHO, ZTEMP, ZRV_IN, ZRV_OUT, ZRC_IN, ZRC_OUT, ZRI_IN, ZRI_OUT, PRT(:,:,:,5), PRT(:,:,:,6), PSIGS,&
+ PMFCONV, ZNCLD, ZSIGRC, OUSERI, OSIGMAS, .FALSE., PSIGQSAT=ZSIGQSAT2D )
+ DEALLOCATE(ZTEMP,ZSIGRC)
+ DEALLOCATE(ZRV_OUT)
+ ELSE
+ ZNCLD=PCLDFR
+ END IF
+ DO JK=IKB,IKE-1
+ JKRAD = JK-JPVEXT
+ DO JJ=IJB,IJE
+ DO JI=IIB,IIE
+ IIJ = 1 + (JI-IIB) + (IIE-IIB+1)*(JJ-IJB)
+ IF ( ZVIEW(IIJ) /= XUNDEF .AND. &
+ (ZRC_OUT(JI,JJ,JK) > 0. .OR. ZRI_OUT(JI,JJ,JK) > 0. ) ) THEN
+ ZFLWP = ZRC_OUT(JI,JJ,JK) / XG /MAX(1.E-10,ZNCLD(JI,JJ,JK)) &
+ * (PPABST(JI,JJ,JK)-PPABST(JI,JJ,JK+1))
+ ZFIWP = ZRI_OUT(JI,JJ,JK) / XG /MAX(1.E-10,ZNCLD(JI,JJ,JK)) &
+ * (PPABST(JI,JJ,JK)-PPABST(JI,JJ,JK+1))
+ ZANGCOR = ZVIEW(IIJ) / 1.66
+ !!!Parametrization following Ou and Chou, 1995 (Atmos. Res.)
+ ZTMP = ZTAVE(IIJ,JKRAD)-XTT !ZTMP in Celsius degree
+ ZRADLP = 326.3+12.42*ZTMP+0.197*(ZTMP**2)+0.0012*(ZTMP**3)
+ ZRADLP = MIN(140., MAX(20., ZRADLP))
+!!! Parametrization following Ebert and Curry, 1992 (JGR-d)
+ ZKI = 0.3 + 1290. / ZRADLP
+ ZCLDLD(IIJ,JKRAD) = ZNCLD(JI,JJ,JK)*(1.-EXP &
+ ( -158.*ZFLWP *ZANGCOR-ZKI*ZFIWP*ZVIEW(IIJ)))
+ ZCLDLU(IIJ,JKRAD) = ZNCLD(JI,JJ,JK)*(1.-EXP &
+ ( -130.*ZFLWP *ZANGCOR-ZKI*ZFIWP*ZVIEW(IIJ)))
+ END IF
+ END DO
+ END DO
+ END DO
+ DEALLOCATE(ZNCLD,ZRC_OUT,ZRI_OUT)
+END IF
+!
+DEALLOCATE(ZEXNT)
+!
+GDOIT_2D(:) = .FALSE.
+!
+! Flags the columns for which the computations have to be performed
+!
+DO JJ=IJB,IJE
+ DO JI=IIB,IIE
+ IIJ = 1 + (JI-IIB) + (IIE-IIB+1)*(JJ-IJB)
+ IF (ZVIEW(IIJ) /= XUNDEF) GDOIT_2D(IIJ) = .TRUE.
+ END DO
+END DO
+IDOIT_COL = COUNT( GDOIT_2D(:) ) ! number of larger scale columns
+!
+GDOIT(:,:) = SPREAD( GDOIT_2D(:),DIM=2,NCOPIES=KFLEV )
+IDOIT = IDOIT_COL*KFLEV
+ALLOCATE(ZWORK1(IDOIT))
+!
+! temperature profiles
+ZWORK1(:) = PACK( ZTAVE(:,:),MASK=GDOIT(:,:) )
+DEALLOCATE(ZTAVE)
+ALLOCATE(ZTAVE(IDOIT_COL,KFLEV))
+ZTAVE(:,:) = RESHAPE( ZWORK1(:),(/IDOIT_COL,KFLEV/) )
+!
+! vapor mixing ratio profiles
+ZWORK1(:) = PACK( ZQVAVE(:,:),MASK=GDOIT(:,:) )
+DEALLOCATE(ZQVAVE)
+ALLOCATE(ZQVAVE(IDOIT_COL,KFLEV))
+ZQVAVE(:,:) = RESHAPE( ZWORK1(:),(/IDOIT_COL,KFLEV/) )
+!
+! cloud emissivities
+ZWORK1(:) = PACK( ZCLDLD(:,:),MASK=GDOIT(:,:) )
+DEALLOCATE(ZCLDLD)
+ALLOCATE(ZCLDLD(IDOIT_COL,KFLEV))
+ZCLDLD(:,:) = RESHAPE( ZWORK1(:),(/IDOIT_COL,KFLEV/) )
+!
+ZWORK1(:) = PACK( ZCLDLU(:,:),MASK=GDOIT(:,:) )
+DEALLOCATE(ZCLDLU)
+ALLOCATE(ZCLDLU(IDOIT_COL,KFLEV))
+ZCLDLU(:,:) = RESHAPE( ZWORK1(:),(/IDOIT_COL,KFLEV/) )
+!
+! ozone content profiles
+ZWORK1(:) = PACK( ZO3AVE(:,:),MASK=GDOIT(:,:) )
+DEALLOCATE(ZO3AVE)
+ALLOCATE(ZO3AVE(IDOIT_COL,KFLEV))
+ZO3AVE(:,:) = RESHAPE( ZWORK1(:),(/IDOIT_COL,KFLEV/) )
+!
+! half-level variables
+ZWORK1(:) = PACK( ZPRES_HL(:,1:KFLEV),MASK=GDOIT(:,:) )
+DEALLOCATE(ZPRES_HL)
+ALLOCATE(ZPRES_HL(IDOIT_COL,KFLEV+1))
+ZPRES_HL(:,1:KFLEV) = RESHAPE( ZWORK1(:),(/IDOIT_COL,KFLEV/) )
+ZPRES_HL(:,KFLEV+1) = PSTATM(IKSTAE,2)*100.0
+!
+ZWORK1(:) = PACK( ZT_HL(:,1:KFLEV),MASK=GDOIT(:,:) )
+DEALLOCATE(ZT_HL)
+ALLOCATE(ZT_HL(IDOIT_COL,KFLEV+1))
+ZT_HL(:,1:KFLEV) = RESHAPE( ZWORK1(:),(/IDOIT_COL,KFLEV/) )
+ZT_HL(:,KFLEV+1) = PSTATM(IKSTAE,3)
+!
+! surface fields
+ALLOCATE(ZWORK3(IDOIT_COL))
+ZWORK3(:) = PACK( ZVIEW(:),MASK=GDOIT_2D(:) )
+DEALLOCATE(ZVIEW)
+ALLOCATE(ZVIEW(IDOIT_COL))
+ZVIEW(:) = ZWORK3(:)
+!
+ZWORK3(:) = PACK( ZREMIS(:),MASK=GDOIT_2D(:) )
+DEALLOCATE(ZREMIS)
+ALLOCATE(ZREMIS(IDOIT_COL))
+ZREMIS(:) = ZWORK3(:)
+!
+ZWORK3(:) = PACK( ZDT0(:),MASK=GDOIT_2D(:) )
+DEALLOCATE(ZDT0)
+ALLOCATE(ZDT0(IDOIT_COL))
+ZDT0(:) = ZWORK3(:)
+!
+DEALLOCATE(ZWORK1)
+DEALLOCATE(ZWORK3)
+!
+! radiation fields
+ALLOCATE(ZRADBC(IDOIT_COL,JPWVINT))
+ALLOCATE(ZRADBT(IDOIT_COL,JPWVINT))
+!
+IDIM = IDOIT_COL
+PRINT *,'KGEO =',KGEO,' IDIM =',IDIM
+!
+!* 6.2 CALLS THE ECMWF_RADIATION ROUTINES
+!
+! ***********************************************************
+! *CAUTION: Routine nbmvec is written in FORTRAN 77*
+! ***********************************************************
+!
+! mixing ratio -> specific humidity conversion
+ZQVAVE(:,:) = ZQVAVE(:,:) / (1.+ZQVAVE(:,:))
+!
+IF( IDIM <= KRAD_COLNBR ) THEN
+ !
+ ! there is less than KRAD_COLNBR verticals to be considered therefore
+ ! no split of the arrays is performed
+ !
+ CALL NBMVEC( 1, IDIM, IDIM, KFLEV, IGL, ICABS, ING1, IUABS, &
+ IH2O, ICO2, IO3, ICNT, IN2O, ICH4, ICO, IC11, IC12, ICFC, &
+ IINIS, IENDS, ICONF, ICLOUD, IOVLP, GPVOIGT, GPTDEP, &
+ ZTAVE, ZQVAVE, ZO3AVE, ZPRES_HL, ZT_HL, &
+ ZVIEW, ZCLDLD, ZCLDLU, ZDT0, ZREMIS, ZRADBC, ZRADBT)
+ELSE
+ !
+ ! the splitting of the arrays will be performed
+ !
+ INUM_CALL = CEILING( REAL( IDIM ) / REAL( KRAD_COLNBR ) )
+ IDIM_RESIDUE = IDIM
+ DO JI_SPLIT = 1 , INUM_CALL
+ IDIM_EFF = MIN( IDIM_RESIDUE,KRAD_COLNBR )
+ !
+ IF( JI_SPLIT == 1 .OR. JI_SPLIT == INUM_CALL ) THEN
+ ALLOCATE( ZREMIS_SPLIT(IDIM_EFF))
+ ALLOCATE( ZO3AVE_SPLIT(IDIM_EFF,KFLEV))
+ ALLOCATE( ZT_HL_SPLIT(IDIM_EFF,KFLEV+1))
+ ALLOCATE( ZPRES_HL_SPLIT(IDIM_EFF,KFLEV+1))
+ ALLOCATE( ZQVAVE_SPLIT(IDIM_EFF,KFLEV))
+ ALLOCATE( ZTAVE_SPLIT(IDIM_EFF,KFLEV))
+ ALLOCATE( ZCLDLU_SPLIT(IDIM_EFF,KFLEV))
+ ALLOCATE( ZCLDLD_SPLIT(IDIM_EFF,KFLEV))
+ ALLOCATE( ZVIEW_SPLIT(IDIM_EFF))
+ ALLOCATE( ZDT0_SPLIT(IDIM_EFF))
+ ALLOCATE( ZRADBT_SPLIT(IDIM_EFF,JPWVINT))
+ ALLOCATE( ZRADBC_SPLIT(IDIM_EFF,JPWVINT))
+ END IF
+ !
+ ! fill the split arrays with their values
+ ! taken from the full arrays
+ !
+ IBEG = IDIM-IDIM_RESIDUE+1
+ IEND = IBEG+IDIM_EFF-1
+ ZREMIS_SPLIT(:) = ZREMIS( IBEG:IEND )
+ ZO3AVE_SPLIT(:,:) = ZO3AVE( IBEG:IEND ,:)
+ ZT_HL_SPLIT(:,:) = ZT_HL( IBEG:IEND ,:)
+ ZPRES_HL_SPLIT(:,:) = ZPRES_HL( IBEG:IEND ,:)
+ ZQVAVE_SPLIT(:,:) = ZQVAVE( IBEG:IEND ,:)
+ ZTAVE_SPLIT(:,:) = ZTAVE ( IBEG:IEND ,:)
+ ZCLDLU_SPLIT(:,:) = ZCLDLU ( IBEG:IEND ,:)
+ ZCLDLD_SPLIT(:,:) = ZCLDLD ( IBEG:IEND ,:)
+ ZVIEW_SPLIT(:) = ZVIEW ( IBEG:IEND )
+ ZDT0_SPLIT(:) = ZDT0 ( IBEG:IEND )
+ !
+ ! call ECMWF_radiation with the split arrays
+ !
+ CALL NBMVEC( 1, IDIM_EFF, IDIM_EFF, KFLEV, IGL, ICABS, ING1, IUABS,&
+ IH2O, ICO2, IO3, ICNT, IN2O, ICH4, ICO, IC11, IC12, ICFC, &
+ IINIS, IENDS, ICONF, ICLOUD, IOVLP, GPVOIGT, GPTDEP, &
+ ZTAVE_SPLIT, ZQVAVE_SPLIT, ZO3AVE_SPLIT, &
+ ZPRES_HL_SPLIT, ZT_HL_SPLIT, &
+ ZVIEW_SPLIT, ZCLDLD_SPLIT, ZCLDLU_SPLIT, ZDT0_SPLIT, &
+ ZREMIS_SPLIT, ZRADBC_SPLIT, ZRADBT_SPLIT)
+ !
+ ! fill the full output arrays with the split arrays
+ !
+ ZRADBT( IBEG:IEND ,:) = ZRADBT_SPLIT(:,:)
+ ZRADBC( IBEG:IEND ,:) = ZRADBC_SPLIT(:,:)
+ !
+ IDIM_RESIDUE = IDIM_RESIDUE - IDIM_EFF
+ !
+ ! desallocation of the split arrays
+ !
+ IF( JI_SPLIT >= INUM_CALL-1 ) THEN
+ DEALLOCATE(ZREMIS_SPLIT)
+ DEALLOCATE(ZO3AVE_SPLIT)
+ DEALLOCATE(ZT_HL_SPLIT)
+ DEALLOCATE(ZPRES_HL_SPLIT)
+ DEALLOCATE(ZQVAVE_SPLIT)
+ DEALLOCATE(ZTAVE_SPLIT)
+ DEALLOCATE(ZCLDLU_SPLIT)
+ DEALLOCATE(ZCLDLD_SPLIT)
+ DEALLOCATE(ZVIEW_SPLIT)
+ DEALLOCATE(ZDT0_SPLIT)
+ DEALLOCATE(ZRADBT_SPLIT)
+ DEALLOCATE(ZRADBC_SPLIT)
+ END IF
+ END DO
+END IF
+!
+DEALLOCATE(ZTAVE,ZQVAVE,ZO3AVE)
+DEALLOCATE(ZPRES_HL,ZT_HL)
+DEALLOCATE(ZREMIS)
+DEALLOCATE(ZDT0)
+DEALLOCATE(ZCLDLD,ZCLDLU)
+DEALLOCATE(ZVIEW)
+!
+ZRADBT = ZRADBT / XPI
+ALLOCATE(ZRADFT(IDIM,JPCAN))
+CALL MAKE_RADSAT(KYEARF, KMONTHF, KDAYF, PSECF, &
+ KGEO, IDIM, ZRADBT, ZRADFT)
+DEALLOCATE(ZRADBT)
+DEALLOCATE(ZRADBC)
+!
+ALLOCATE(ZWORK1(IDIM*JPCAN))
+ZWORK1(:) = PACK( ZRADFT(:,:),MASK=.TRUE. )
+ALLOCATE(ZZRADFT(KDLON,JPCAN))
+ZZRADFT(:,:) = UNPACK( ZWORK1(:),MASK=GDOIT(:,1:JPCAN),FIELD=XUNDEF )
+DEALLOCATE(ZRADFT)
+DEALLOCATE(ZWORK1)
+!
+PIRBT = XUNDEF
+PWVBT = XUNDEF
+DO JJ=IJB,IJE
+ DO JI=IIB,IIE
+ IIJ = 1 + (JI-IIB) + (IIE-IIB+1)*(JJ-IJB)
+ PIRBT(JI,JJ) = ZZRADFT(IIJ,1)
+ PWVBT(JI,JJ) = ZZRADFT(IIJ,2)
+ END DO
+END DO
+DEALLOCATE(ZZRADFT)
+!
+END SUBROUTINE RADTR_SATEL
diff --git a/src/mesonh/mode_qsatmx_tab.F90 b/src/mesonh/mode_qsatmx_tab.F90
new file mode 100644
index 0000000000000000000000000000000000000000..01d697b19bdeeb1cc011c4826e8c37da037cb944
--- /dev/null
+++ b/src/mesonh/mode_qsatmx_tab.F90
@@ -0,0 +1,27 @@
+MODULE MODE_QSATMX_TAB
+IMPLICIT NONE
+CONTAINS
+FUNCTION QSATMX_TAB(P,T,FICE)
+
+ USE PARKIND1, ONLY : JPRB
+ USE MODD_CST ,ONLY : XEPSILO
+ USE MODE_TIWMX, ONLY : ESATI,ESATW
+
+ IMPLICIT NONE
+
+ REAL :: QSATMX_TAB
+ REAL, INTENT(IN) :: P,T,FICE
+
+ REAL :: ZES
+
+ ZES = ESATI(T)*FICE + ESATW(T)*(1.-FICE)
+ IF(ZES >= P)THEN ! temp > boiling point, condensation not possible.
+ ! Then this function lacks physical meaning,
+ ! here set to one
+ QSATMX_TAB = 1.
+ ELSE
+ QSATMX_TAB = XEPSILO*ZES/(P-ZES) !r
+ ENDIF
+
+END FUNCTION QSATMX_TAB
+END MODULE MODE_QSATMX_TAB
diff --git a/src/mesonh/turb/th_r_from_thl_rt_1d.f90 b/src/mesonh/turb/th_r_from_thl_rt_1d.f90
index f89a211d69d0e5a815f2d8f31da4bbb0d39c55a0..dbb4d564a66e743470b2204ecd9d0d813b092cc2 100644
--- a/src/mesonh/turb/th_r_from_thl_rt_1d.f90
+++ b/src/mesonh/turb/th_r_from_thl_rt_1d.f90
@@ -104,6 +104,7 @@ REAL, DIMENSION(SIZE(PTHL,1)) :: ZRVSAT,ZCPH,ZRLTEMP,ZCPH2
REAL, DIMENSION(SIZE(PTHL,1)) :: ZT,ZLVOCPEXN,ZLSOCPEXN
REAL, DIMENSION(SIZE(PTHL,1)) :: ZDRSATODT,ZDRSATODTW,ZDRSATODTI
REAL, DIMENSION(SIZE(PTHL,1)) :: ZFOESW, ZFOESI
+INTEGER, DIMENSION(SIZE(PTHL,1)) :: IERR
!----------------------------------------------------------------------------
!
!* 1 Initialisation
@@ -147,8 +148,7 @@ DO II=1,JITER
WHERE(PRL(:)+PRI(:) > 1.E-20)
PFRAC_ICE(:) = PRI(:) / (PRL(:)+PRI(:))
ENDWHERE
- CALL COMPUTE_FRAC_ICE(SIZE(PFRAC_ICE,1),HFRAC_ICE,PFRAC_ICE(:),ZT(:))
-
+ CALL COMPUTE_FRAC_ICE(HFRAC_ICE,PFRAC_ICE(:),ZT(:), IERR(:))
!Computation of Rvsat and dRsat/dT
!In this version QSAT, QSATI, DQSAT and DQASATI functions are not used
!due to performance issue