diff --git a/src/common/micro/compute_frac_ice.func.h b/src/common/micro/compute_frac_ice.func.h
index f63cfe5548499d148d97b967ea89afe934f7857c..8c6d4e617d519e2277d3a7defe3b11c95513cafc 100644
--- a/src/common/micro/compute_frac_ice.func.h
+++ b/src/common/micro/compute_frac_ice.func.h
@@ -30,13 +30,13 @@ CHARACTER(LEN=1), INTENT(IN) :: HFRAC_ICE ! scheme to use
TYPE(NEB_t), INTENT(IN) :: NEB
REAL, INTENT(IN) :: PT ! temperature
REAL, INTENT(INOUT) :: PFRAC_ICE ! Ice fraction (1 for ice only, 0 for liquid only)
-INTEGER, INTENT(OUT) :: KERR ! Error code in return
+INTEGER, OPTIONAL, INTENT(OUT) :: KERR ! Error code in return
!
!------------------------------------------------------------------------
! 1. Compute FRAC_ICE
!
-KERR=0
+IF (PRESENT(KERR)) KERR=0
SELECT CASE(HFRAC_ICE)
CASE ('T') !using Temperature
PFRAC_ICE = MAX( 0., MIN(1., (( NEB%XTMAXMIX - PT ) / ( NEB%XTMAXMIX - NEB%XTMINMIX )) ) ) ! freezing interval
@@ -48,7 +48,7 @@ SELECT CASE(HFRAC_ICE)
! (almost) nothing to do
PFRAC_ICE = MAX( 0., MIN(1., PFRAC_ICE ) )
CASE DEFAULT
- KERR=1
+ IF (PRESENT(KERR)) KERR=1
END SELECT
END SUBROUTINE COMPUTE_FRAC_ICE
diff --git a/src/common/turb/mode_compute_entr_detr.F90 b/src/common/turb/mode_compute_entr_detr.F90
index e0a5add791f07b7c2c43a1269580ccb8d1fedb81..9315bb751087c4f0eb781c4ea1bb600dee3c2ab4 100644
--- a/src/common/turb/mode_compute_entr_detr.F90
+++ b/src/common/turb/mode_compute_entr_detr.F90
@@ -79,8 +79,6 @@ USE MODD_CST, ONLY: CST_t
USE MODD_NEB, ONLY: NEB_t
USE MODD_PARAM_MFSHALL_n, ONLY: PARAM_MFSHALL_t
!
-USE MODE_TH_R_FROM_THL_RT_1D, ONLY: TH_R_FROM_THL_RT_1D
-!
USE PARKIND1, ONLY : JPRB
USE YOMHOOK , ONLY : LHOOK, DR_HOOK
@@ -172,6 +170,7 @@ REAL, DIMENSION(D%NIT) :: ZDZ_STOP,& ! Exact Height of the LCL above f
ZTHV_PLUS_HALF ! Thv at flux point(kk+kkl)
REAL :: ZDZ ! Delta Z used in computations
INTEGER :: JI, JLOOP
+REAL, DIMENSION(D%NIT, 16) :: ZBUF
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!----------------------------------------------------------------------------------
@@ -286,10 +285,11 @@ ENDDO
!but we are dealing with updraft and not mixture
ZRCMIX(:)=PRC_UP(:)
ZRIMIX(:)=PRI_UP(:)
- CALL TH_R_FROM_THL_RT_1D(CST,NEB,D%NIT,HFRAC_ICE,ZFRAC_ICE,&
+ CALL TH_R_FROM_THL_RT(CST,NEB,D%NIT,HFRAC_ICE,ZFRAC_ICE,&
PPRE_PLUS_HALF,PTHL_UP,PRT_UP,&
ZTHMIX,ZRVMIX,ZRCMIX,ZRIMIX,&
- ZRSATW, ZRSATI,OOCEAN=.FALSE.)
+ ZRSATW, ZRSATI,OOCEAN=.FALSE.,&
+ PBUF=ZBUF)
ZTHV_UP_F2(:) = ZTHMIX(:)*(1.+ZRVORD*ZRVMIX(:))/(1.+PRT_UP(:))
! Integral buoyancy for cloudy part
@@ -363,19 +363,21 @@ DO JLOOP=1,SIZE(OTEST)
ZMIXRT(JLOOP) = 0.1
ENDIF
ENDDO
- CALL TH_R_FROM_THL_RT_1D(CST,NEB,D%NIT,HFRAC_ICE,ZFRAC_ICE,&
+ CALL TH_R_FROM_THL_RT(CST,NEB,D%NIT,HFRAC_ICE,ZFRAC_ICE,&
ZPRE,ZMIXTHL,ZMIXRT,&
ZTHMIX,ZRVMIX,PRC_MIX,PRI_MIX,&
- ZRSATW, ZRSATI,OOCEAN=.FALSE.)
+ ZRSATW, ZRSATI,OOCEAN=.FALSE.,&
+ PBUF=ZBUF)
ZTHVMIX(:) = ZTHMIX(:)*(1.+ZRVORD*ZRVMIX(:))/(1.+ZMIXRT(:))
! Compute cons then non cons. var. of mixture at the flux level KK+KKL with initial ZKIC
ZMIXTHL(:) = ZKIC_INIT * 0.5*(PTHLM(:,KK)+PTHLM(:,KK+KKL))+(1. - ZKIC_INIT)*PTHL_UP(:)
ZMIXRT(:) = ZKIC_INIT * 0.5*(PRTM(:,KK)+PRTM(:,KK+KKL))+(1. - ZKIC_INIT)*PRT_UP(:)
- CALL TH_R_FROM_THL_RT_1D(CST,NEB,D%NIT,HFRAC_ICE,ZFRAC_ICE,&
+ CALL TH_R_FROM_THL_RT(CST,NEB,D%NIT,HFRAC_ICE,ZFRAC_ICE,&
PPRE_PLUS_HALF,ZMIXTHL,ZMIXRT,&
ZTHMIX,ZRVMIX,PRC_MIX,PRI_MIX,&
- ZRSATW, ZRSATI,OOCEAN=.FALSE.)
+ ZRSATW, ZRSATI,OOCEAN=.FALSE.,&
+ PBUF=ZBUF)
ZTHVMIX_F2(:) = ZTHMIX(:)*(1.+ZRVORD*ZRVMIX(:))/(1.+ZMIXRT(:))
!Computation of mean ZKIC over the cloudy part
@@ -445,5 +447,8 @@ DO JLOOP=1,SIZE(OTEST)
ENDDO
IF (LHOOK) CALL DR_HOOK('COMPUTE_ENTR_DETR',1,ZHOOK_HANDLE)
+CONTAINS
+INCLUDE "th_r_from_thl_rt.func.h"
+INCLUDE "compute_frac_ice.func.h"
END SUBROUTINE COMPUTE_ENTR_DETR
END MODULE MODE_COMPUTE_ENTR_DETR
diff --git a/src/common/turb/mode_compute_mf_cloud_bigaus.F90 b/src/common/turb/mode_compute_mf_cloud_bigaus.F90
index bfe4f3a91e5014bfe13f43b342ad8209e04f6c62..da0e4f92bb261bd7e98c8221284b574255498d60 100644
--- a/src/common/turb/mode_compute_mf_cloud_bigaus.F90
+++ b/src/common/turb/mode_compute_mf_cloud_bigaus.F90
@@ -121,7 +121,7 @@ CALL MZF_MF(D, PFRAC_ICE_UP(:,:), ZFRAC_ICE_UP_M(:,:))
!computation of omega star up
ZOMEGA_UP_M(:)=0.
-DO JK=D%NKB,D%NKE,D%NKL
+DO JK=D%NKB,D%NKE-D%NKL,D%NKL
!Vertical integration over the entire column but only buoyant points are used
!ZOMEGA_UP_M(:)=ZOMEGA_UP_M(:) + &
! ZEMF_M(:,JK) * &
diff --git a/src/common/turb/mode_compute_updraft.F90 b/src/common/turb/mode_compute_updraft.F90
index 76c8df61cea7da482e8713da09a24d16b0b3337f..f5e2e5a74ce253e2e13d4f510ccf07bf2a2f58ac 100644
--- a/src/common/turb/mode_compute_updraft.F90
+++ b/src/common/turb/mode_compute_updraft.F90
@@ -74,7 +74,6 @@ USE MODD_TURB_n, ONLY: TURB_t
USE MODD_CTURB, ONLY: CSTURB_t
!
USE MODE_COMPUTE_ENTR_DETR, ONLY: COMPUTE_ENTR_DETR
-USE MODE_TH_R_FROM_THL_RT_1D, ONLY: TH_R_FROM_THL_RT_1D
USE MODI_SHUMAN_MF, ONLY: MZM_MF, MZF_MF, GZ_M_W_MF
USE MODE_COMPUTE_BL89_ML, ONLY: COMPUTE_BL89_ML
@@ -194,6 +193,7 @@ REAL, DIMENSION(D%NIT) :: ZSURF
REAL, DIMENSION(D%NIT,D%NKT) :: ZSHEAR,ZDUDZ,ZDVDZ ! vertical wind shear
!
REAL, DIMENSION(D%NIT,D%NKT) :: ZWK
+REAL, DIMENSION(D%NIT,16) :: ZBUF
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
@@ -298,9 +298,10 @@ IF (OENTR_DETR) THEN
! (all or nothing ajustement)
PRC_UP(:,D%NKB)=0.
PRI_UP(:,D%NKB)=0.
- CALL TH_R_FROM_THL_RT_1D(CST, NEB, D%NIT, HFRAC_ICE,PFRAC_ICE_UP(:,D%NKB),ZPRES_F(:,D%NKB), &
+ CALL TH_R_FROM_THL_RT(CST, NEB, D%NIT, HFRAC_ICE,PFRAC_ICE_UP(:,D%NKB),ZPRES_F(:,D%NKB), &
PTHL_UP(:,D%NKB),PRT_UP(:,D%NKB),ZTH_UP(:,D%NKB), &
- PRV_UP(:,D%NKB),PRC_UP(:,D%NKB),PRI_UP(:,D%NKB),ZRSATW(:),ZRSATI(:),OOCEAN=.FALSE.)
+ PRV_UP(:,D%NKB),PRC_UP(:,D%NKB),PRI_UP(:,D%NKB),ZRSATW(:),ZRSATI(:), OOCEAN=.FALSE., &
+ PBUF=ZBUF(:,:))
! compute updraft thevav and buoyancy term at KKB level
PTHV_UP(:,D%NKB) = ZTH_UP(:,D%NKB)*((1+ZRVORD*PRV_UP(:,D%NKB))/(1+PRT_UP(:,D%NKB)))
@@ -492,9 +493,10 @@ DO JK=D%NKB,D%NKE-D%NKL,D%NKL
! Compute non cons. var. at level JK+KKL
ZRC_UP(:)=PRC_UP(:,JK) ! guess = level just below
ZRI_UP(:)=PRI_UP(:,JK) ! guess = level just below
- CALL TH_R_FROM_THL_RT_1D(CST, NEB, D%NIT, HFRAC_ICE,PFRAC_ICE_UP(:,JK+D%NKL),ZPRES_F(:,JK+D%NKL), &
+ CALL TH_R_FROM_THL_RT(CST, NEB, D%NIT, HFRAC_ICE,PFRAC_ICE_UP(:,JK+D%NKL),ZPRES_F(:,JK+D%NKL), &
PTHL_UP(:,JK+D%NKL),PRT_UP(:,JK+D%NKL),ZTH_UP(:,JK+D%NKL), &
- ZRV_UP(:),ZRC_UP(:),ZRI_UP(:),ZRSATW(:),ZRSATI(:), OOCEAN=.FALSE.)
+ ZRV_UP(:),ZRC_UP(:),ZRI_UP(:),ZRSATW(:),ZRSATI(:), OOCEAN=.FALSE., &
+ PBUF=ZBUF(:,:))
WHERE(GTEST)
PRC_UP(:,JK+D%NKL)=ZRC_UP(:)
PRV_UP(:,JK+D%NKL)=ZRV_UP(:)
@@ -588,5 +590,8 @@ IF(OENTR_DETR) THEN
ENDIF
IF (LHOOK) CALL DR_HOOK('COMPUTE_UPDRAFT',1,ZHOOK_HANDLE)
+CONTAINS
+INCLUDE "th_r_from_thl_rt.func.h"
+INCLUDE "compute_frac_ice.func.h"
END SUBROUTINE COMPUTE_UPDRAFT
END MODULE MODE_COMPUTE_UPDRAFT
diff --git a/src/common/turb/mode_compute_updraft_raha.F90 b/src/common/turb/mode_compute_updraft_raha.F90
index 3d35ad43ac14ccbb7a2db3975e8d75cd22da4544..f6818168e1722e9383e6aa1469a6e2c25f2f6a08 100644
--- a/src/common/turb/mode_compute_updraft_raha.F90
+++ b/src/common/turb/mode_compute_updraft_raha.F90
@@ -63,7 +63,6 @@ USE MODD_CST, ONLY: CST_t
USE MODD_NEB, ONLY: NEB_t
USE MODD_PARAM_MFSHALL_n, ONLY: PARAM_MFSHALL_t
!
-USE MODE_TH_R_FROM_THL_RT_1D, ONLY: TH_R_FROM_THL_RT_1D
USE MODI_SHUMAN_MF, ONLY: MZM_MF
!
USE PARKIND1, ONLY : JPRB
@@ -187,6 +186,7 @@ REAL, DIMENSION(D%NIT) :: ZA,ZB,ZQTM,ZQT_UP
REAL :: ZDEPTH_MAX1, ZDEPTH_MAX2 ! control auto-extinction process
REAL :: ZTMAX,ZRMAX, ZEPS ! control value
+REAL, DIMENSION(D%NIT,16) :: ZBUF
REAL(KIND=JPRB) :: ZHOOK_HANDLE
IF (LHOOK) CALL DR_HOOK('COMPUTE_UPDRAF_RAHA',0,ZHOOK_HANDLE)
@@ -293,9 +293,10 @@ GTEST = (ZW_UP2(:,D%NKB) > ZEPS)
PRC_UP(:,D%NKB)=0.
PRI_UP(:,D%NKB)=0.
-CALL TH_R_FROM_THL_RT_1D(CST,NEB, D%NIT, HFRAC_ICE,PFRAC_ICE_UP(:,D%NKB),ZPRES_F(:,D%NKB), &
+CALL TH_R_FROM_THL_RT(CST,NEB, D%NIT, HFRAC_ICE,PFRAC_ICE_UP(:,D%NKB),ZPRES_F(:,D%NKB), &
PTHL_UP(:,D%NKB),PRT_UP(:,D%NKB),ZTH_UP(:,D%NKB), &
- PRV_UP(:,D%NKB),PRC_UP(:,D%NKB),PRI_UP(:,D%NKB),ZRSATW(:),ZRSATI(:),OOCEAN=.FALSE.)
+ PRV_UP(:,D%NKB),PRC_UP(:,D%NKB),PRI_UP(:,D%NKB),ZRSATW(:),ZRSATI(:),OOCEAN=.FALSE.,&
+ PBUF=ZBUF)
! compute updraft thevav and buoyancy term at KKB level
PTHV_UP(:,D%NKB) = ZTH_UP(:,D%NKB)*((1+ZRVORD*PRV_UP(:,D%NKB))/(1+PRT_UP(:,D%NKB)))
@@ -472,9 +473,10 @@ DO JK=D%NKB,D%NKE-D%NKL,D%NKL
ZRC_UP(:)=PRC_UP(:,JK) ! guess = level just below
ZRI_UP(:)=PRI_UP(:,JK) ! guess = level just below
ZRV_UP(:)=PRV_UP(:,JK)
- CALL TH_R_FROM_THL_RT_1D(CST,NEB, D%NIT, HFRAC_ICE,PFRAC_ICE_UP(:,JK+D%NKL),ZPRES_F(:,JK+D%NKL), &
+ CALL TH_R_FROM_THL_RT(CST,NEB, D%NIT, HFRAC_ICE,PFRAC_ICE_UP(:,JK+D%NKL),ZPRES_F(:,JK+D%NKL), &
PTHL_UP(:,JK+D%NKL),PRT_UP(:,JK+D%NKL),ZTH_UP(:,JK+D%NKL), &
- ZRV_UP(:),ZRC_UP(:),ZRI_UP(:),ZRSATW(:),ZRSATI(:),OOCEAN=.FALSE.)
+ ZRV_UP(:),ZRC_UP(:),ZRI_UP(:),ZRSATW(:),ZRSATI(:),OOCEAN=.FALSE.,&
+ PBUF=ZBUF)
WHERE(GTEST)
ZT_UP(:) = ZTH_UP(:,JK+D%NKL)*PEXNM(:,JK+D%NKL)
ZCP(:) = CST%XCPD + CST%XCL * ZRC_UP(:)
@@ -585,6 +587,10 @@ WHERE (GWORK2)
ENDWHERE
IF (LHOOK) CALL DR_HOOK('COMPUTE_UPDRAF_RAHA',1,ZHOOK_HANDLE)
-
+!
+CONTAINS
+INCLUDE "th_r_from_thl_rt.func.h"
+INCLUDE "compute_frac_ice.func.h"
+!
END SUBROUTINE COMPUTE_UPDRAFT_RAHA
END MODULE MODE_COMPUTE_UPDRAFT_RAHA
diff --git a/src/common/turb/mode_compute_updraft_rhcj10.F90 b/src/common/turb/mode_compute_updraft_rhcj10.F90
index 57b31c12da1e44ce6770fb825cbaccddf045f90d..b5bfe8cf2f6965333331f8e898efee07a5c25f52 100644
--- a/src/common/turb/mode_compute_updraft_rhcj10.F90
+++ b/src/common/turb/mode_compute_updraft_rhcj10.F90
@@ -68,7 +68,6 @@ USE MODD_PARAM_MFSHALL_n, ONLY: PARAM_MFSHALL_t
USE MODD_TURB_n, ONLY: TURB_t
USE MODD_CTURB, ONLY: CSTURB_t
!
-USE MODE_TH_R_FROM_THL_RT_1D, ONLY: TH_R_FROM_THL_RT_1D
USE MODI_SHUMAN_MF, ONLY: MZF_MF, MZM_MF, GZ_M_W_MF
USE MODE_COMPUTE_BL89_ML, ONLY: COMPUTE_BL89_ML
@@ -192,6 +191,7 @@ REAL :: ZTMAX,ZRMAX, ZEPS ! control value
REAL, DIMENSION(D%NIT,D%NKT) :: ZSHEAR,ZDUDZ,ZDVDZ ! vertical wind shear
!
REAL, DIMENSION(D%NIT,D%NKT) :: ZWK
+REAL, DIMENSION(D%NIT,16) :: ZBUF
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
IF (LHOOK) CALL DR_HOOK('COMPUTE_UPDRAFT_RHCJ10',0,ZHOOK_HANDLE)
@@ -310,9 +310,10 @@ ZW_UP2(:,D%NKB) = MAX(0.0001,(2./3.)*ZTKEM_F(:,D%NKB))
PRC_UP(:,D%NKB)=0.
PRI_UP(:,D%NKB)=0.
-CALL TH_R_FROM_THL_RT_1D(CST,NEB,D%NIT,HFRAC_ICE,PFRAC_ICE_UP(:,D%NKB),ZPRES_F(:,D%NKB), &
+CALL TH_R_FROM_THL_RT(CST,NEB,D%NIT,HFRAC_ICE,PFRAC_ICE_UP(:,D%NKB),ZPRES_F(:,D%NKB), &
PTHL_UP(:,D%NKB),PRT_UP(:,D%NKB),ZTH_UP(:,D%NKB), &
- PRV_UP(:,D%NKB),PRC_UP(:,D%NKB),PRI_UP(:,D%NKB),ZRSATW(:),ZRSATI(:),OOCEAN=.FALSE.)
+ PRV_UP(:,D%NKB),PRC_UP(:,D%NKB),PRI_UP(:,D%NKB),ZRSATW(:),ZRSATI(:),OOCEAN=.FALSE.,&
+ PBUF=ZBUF)
DO JI=1,D%NIT
! compute updraft thevav and buoyancy term at KKB level
@@ -418,9 +419,10 @@ DO JK=D%NKB,D%NKE-D%NKL,D%NKL
ZRC_UP(:) =PRC_UP(:,JK) ! guess
ZRI_UP(:) =PRI_UP(:,JK) ! guess
ZRV_UP(:) =PRV_UP(:,JK)
- CALL TH_R_FROM_THL_RT_1D(CST,NEB, D%NIT, HFRAC_ICE,PFRAC_ICE_UP(:,JK),&
+ CALL TH_R_FROM_THL_RT(CST,NEB, D%NIT, HFRAC_ICE,PFRAC_ICE_UP(:,JK),&
PPABSM(:,JK),PTHL_UP(:,JK),PRT_UP(:,JK),&
- ZTH_UP(:,JK),ZRV_UP,ZRC_UP,ZRI_UP,ZRSATW(:),ZRSATI(:),OOCEAN=.FALSE.)
+ ZTH_UP(:,JK),ZRV_UP,ZRC_UP,ZRI_UP,ZRSATW(:),ZRSATI(:),OOCEAN=.FALSE.,&
+ PBUF=ZBUF)
DO JI=1,D%NIT
IF (GTEST(JI)) THEN
@@ -518,9 +520,10 @@ DO JK=D%NKB,D%NKE-D%NKL,D%NKL
ZRC_UP(:)=PRC_UP(:,JK) ! guess = level just below
ZRI_UP(:)=PRI_UP(:,JK) ! guess = level just below
ZRV_UP(:)=PRV_UP(:,JK)
- CALL TH_R_FROM_THL_RT_1D(CST,NEB, D%NIT, HFRAC_ICE,PFRAC_ICE_UP(:,JK+D%NKL),ZPRES_F(:,JK+D%NKL), &
+ CALL TH_R_FROM_THL_RT(CST,NEB, D%NIT, HFRAC_ICE,PFRAC_ICE_UP(:,JK+D%NKL),ZPRES_F(:,JK+D%NKL), &
PTHL_UP(:,JK+D%NKL),PRT_UP(:,JK+D%NKL),ZTH_UP(:,JK+D%NKL), &
- ZRV_UP(:),ZRC_UP(:),ZRI_UP(:),ZRSATW(:),ZRSATI(:),OOCEAN=.FALSE.)
+ ZRV_UP(:),ZRC_UP(:),ZRI_UP(:),ZRSATW(:),ZRSATI(:),OOCEAN=.FALSE.,&
+ PBUF=ZBUF)
DO JI=1,D%NIT
IF(GTEST(JI)) THEN
@@ -614,6 +617,10 @@ DO JK=1, D%NKT
ENDDO
IF (LHOOK) CALL DR_HOOK('COMPUTE_UPDRAFT_RHCJ10',1,ZHOOK_HANDLE)
-
+!
+CONTAINS
+INCLUDE "th_r_from_thl_rt.func.h"
+INCLUDE "compute_frac_ice.func.h"
+!
END SUBROUTINE COMPUTE_UPDRAFT_RHCJ10
END MODULE MODE_COMPUTE_UPDRAFT_RHCJ10
diff --git a/src/common/turb/mode_th_r_from_thl_rt_1d.F90 b/src/common/turb/mode_th_r_from_thl_rt_1d.F90
deleted file mode 100644
index 1c974d8abab215de0b53bc4e82a2a465f781cc9f..0000000000000000000000000000000000000000
--- a/src/common/turb/mode_th_r_from_thl_rt_1d.F90
+++ /dev/null
@@ -1,211 +0,0 @@
-!MNH_LIC Copyright 2006-2022 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 MODE_TH_R_FROM_THL_RT_1D
-IMPLICIT NONE
-CONTAINS
- SUBROUTINE TH_R_FROM_THL_RT_1D(CST, NEB, KT, HFRAC_ICE,PFRAC_ICE,PP, &
- PTHL, PRT, PTH, PRV, PRL, PRI, &
- PRSATW, PRSATI, PRR, PRS, PRG, PRH,OOCEAN)
-! #################################################################
-!
-!
-!!**** *TH_R_FROM_THL_RT_1D* - computes the non-conservative variables
-!! from conservative variables
-!!
-!! PURPOSE
-!! -------
-!!
-!!** METHOD
-!! ------
-!!
-!!
-!! EXTERNAL
-!! --------
-!!
-!! IMPLICIT ARGUMENTS
-!! ------------------
-!!
-!!
-!! REFERENCE
-!! ---------
-!!
-!! AUTHOR
-!! ------
-!! Julien PERGAUD * Meteo-France *
-!!
-!! MODIFICATIONS
-!! -------------
-!! Original 13/03/06
-!! S. Riette April 2011 : ice added, allow ZRLTEMP to be negative
-!! we use dQsat/dT to help convergence
-!! use of optional PRR, PRS, PRG, PRH
-!! S. Riette Nov 2016: support for HFRAC_ICE='S'
-!!
-!! --------------------------------------------------------------------------
-!
-!* 0. DECLARATIONS
-! ------------
-!
-USE PARKIND1, ONLY : JPRB
-USE YOMHOOK , ONLY : LHOOK, DR_HOOK
-USE MODD_CST, ONLY : CST_t
-USE MODD_NEB, ONLY : NEB_t
-USE MODE_THERMO
-!
-IMPLICIT NONE
-!
-!
-!* 0.1 declarations of arguments
-!
-TYPE(CST_t), INTENT(IN) :: CST
-TYPE(NEB_t), INTENT(IN) :: NEB
-INTEGER, INTENT(IN) :: KT
-CHARACTER(LEN=1), INTENT(IN) :: HFRAC_ICE
-LOGICAL, INTENT(IN) :: OOCEAN ! switch for Ocean model version
-REAL, DIMENSION(KT), INTENT(INOUT) :: PFRAC_ICE
-REAL, DIMENSION(KT), INTENT(IN) :: PP ! Pressure
-REAL, DIMENSION(KT), INTENT(IN) :: PTHL ! thetal to transform into th
-REAL, DIMENSION(KT),INTENT(IN) :: PRT ! Total mixing ratios to transform into rv,rc and ri
-REAL, DIMENSION(KT),OPTIONAL,INTENT(IN) :: PRR, PRS, PRG, PRH
-REAL, DIMENSION(KT), INTENT(OUT):: PTH ! th
-REAL, DIMENSION(KT), INTENT(OUT):: PRV ! vapor mixing ratio
-REAL, DIMENSION(KT), INTENT(INOUT):: PRL ! vapor mixing ratio
-REAL, DIMENSION(KT), INTENT(INOUT):: PRI ! vapor mixing ratio
-REAL, DIMENSION(KT), INTENT(OUT) :: PRSATW ! estimated mixing ration at saturation over water
-REAL, DIMENSION(KT), INTENT(OUT) :: PRSATI ! estimated mixing ration at saturation over ice
-!
-!-------------------------------------------------------------------------------
-!
-! 0.2 declaration of local variables
-INTEGER :: II ! Loop control
-INTEGER :: JITER ! number of iterations
-INTEGER :: J
-REAL, DIMENSION(KT) :: ZEXN
-REAL, DIMENSION(KT) :: ZRVSAT,ZCPH,ZRLTEMP,ZCPH2
-REAL, DIMENSION(KT) :: ZT,ZLVOCPEXN,ZLSOCPEXN
-REAL, DIMENSION(KT) :: ZDRSATODT,ZDRSATODTW,ZDRSATODTI
-REAL, DIMENSION(KT) :: ZFOESW, ZFOESI
-REAL, DIMENSION(KT) :: ZLOGT, Z99PP, Z1PRT
-REAL :: ZVAR1, ZVAR2, ZTPOW2, ZDELT
-INTEGER, DIMENSION(KT) :: IERR
-
-REAL(KIND=JPRB) :: ZHOOK_HANDLE
-!----------------------------------------------------------------------------
-!
-!* 1 Initialisation
-! --------------
-!
-!
-IF (LHOOK) CALL DR_HOOK('TH_R_FROM_THL_RT_1D',0,ZHOOK_HANDLE)
-!
-!Number of iterations
-JITER=2
-!
-!Computation of ZCPH2 depending on dummy arguments received
-ZCPH2(:)=0
-IF(PRESENT(PRR)) ZCPH2(:)=ZCPH2(:) + CST%XCL*PRR(:)
-IF(PRESENT(PRS)) ZCPH2(:)=ZCPH2(:) + CST%XCI*PRS(:)
-IF(PRESENT(PRG)) ZCPH2(:)=ZCPH2(:) + CST%XCI*PRG(:)
-IF(PRESENT(PRH)) ZCPH2(:)=ZCPH2(:) + CST%XCI*PRH(:)
-!
-!Computation of an approximate state thanks to PRL and PRI guess
-ZEXN(:)=(PP(:)/CST%XP00) ** CST%RDSCPD
-
-DO J=1,KT
- Z99PP(J)=0.99*PP(J)
- PRV(J)=PRT(J)-PRL(J)-PRI(J)
- ZCPH(J)=CST%XCPD+ CST%XCPV * PRV(J)+ CST%XCL * PRL(J) + CST%XCI * PRI(J) + ZCPH2(J)
- ZVAR2=ZCPH(J)*ZEXN(J)
- ZDELT=(PTHL(J)*ZEXN(J))-CST%XTT
- ZLVOCPEXN(J) = (CST%XLVTT + (CST%XCPV-CST%XCL) * ZDELT) /ZVAR2
- ZLSOCPEXN(J) = (CST%XLSTT + (CST%XCPV-CST%XCI) * ZDELT) /ZVAR2
- PTH(J)=PTHL(J)+ZLVOCPEXN(J)*PRL(J)+ZLSOCPEXN(J)*PRI(J)
- Z1PRT(J)=1+PRT(J)
-ENDDO
-!
-!
-! 2 Iteration
-! ---------
-
-DO II=1,JITER
- IF (OOCEAN) THEN
- ZT=PTH
- ELSE
- ZT(:)=PTH(:)*ZEXN(:)
- END IF
- !Computation of liquid/ice fractions
- PFRAC_ICE(:) = 0.
- DO J=1, KT
- IF(PRL(J)+PRI(J) > 1.E-20) THEN
- PFRAC_ICE(J) = PRI(J) / (PRL(J)+PRI(J))
- ENDIF
- ENDDO
- CALL COMPUTE_FRAC_ICE(HFRAC_ICE,NEB,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
-
- ! Log does not vectorize on all compilers:
- ZLOGT(:)=LOG(ZT(:))
-
- DO J=1,KT
- ZFOESW(J) = MIN(EXP( CST%XALPW - CST%XBETAW/ZT(J) - CST%XGAMW*ZLOGT(J) ), Z99PP(J))
- ZFOESI(J) = MIN(EXP( CST%XALPI - CST%XBETAI/ZT(J) - CST%XGAMI*ZLOGT(J) ), Z99PP(J))
- PRSATW(J) = CST%XRD/CST%XRV*ZFOESW(J)/PP(J) / (1.+(CST%XRD/CST%XRV-1.)*ZFOESW(J)/PP(J))
- PRSATI(J) = CST%XRD/CST%XRV*ZFOESI(J)/PP(J) / (1.+(CST%XRD/CST%XRV-1.)*ZFOESI(J)/PP(J))
- ZTPOW2=ZT(J)**2
- ZDRSATODTW(J) = PRSATW(J) / (1.+(CST%XRD/CST%XRV-1.)*ZFOESW(J)/PP(J) ) &
- * (CST%XBETAW/ZTPOW2 - CST%XGAMW/ZT(J))*Z1PRT(J)
- ZDRSATODTI(J) = PRSATI(J) / (1.+(CST%XRD/CST%XRV-1.)*ZFOESI(J)/PP(J) ) &
- * (CST%XBETAI/ZTPOW2 - CST%XGAMI/ZT(J))*Z1PRT(J)
- !PRSATW(J) = QSAT(ZT(J),PP(J)) !qsatw
- !PRSATI(J) = QSATI(ZT(J),PP(J)) !qsati
- !ZDRSATODTW(J) = DQSAT(ZT(J),PP(J),PRSATW(J))*Z1PRT(J)
- !ZDRSATODTI(J) = DQSATI(ZT(J),PP(J),PRSATI(J))*Z1PRT(J)
- PRSATW(J) = PRSATW(J)*Z1PRT(J)
- PRSATI(J) = PRSATI(J)*Z1PRT(J)
- ZRVSAT(J) = PRSATW(J)*(1-PFRAC_ICE(J)) + PRSATI(J)*PFRAC_ICE(J)
- ZDRSATODT(J) = (ZDRSATODTW(J)*(1-PFRAC_ICE(J))+ &
- & ZDRSATODTI(J)*PFRAC_ICE(J))
-
- !Computation of new PRL, PRI and PRV
- !Correction term applied to (PRV(J)-ZRVSAT(J)) is computed assuming that
- !ZLVOCPEXN, ZLSOCPEXN and ZCPH don't vary to much with T. It takes into account
- !the variation (estimated linear) of Qsat with T
- ZRLTEMP(J)=(PRV(J)-ZRVSAT(J))/ &
- &(1 + ZDRSATODT(J)*ZEXN(J)* &
- & (ZLVOCPEXN(J)*(1-PFRAC_ICE(J))+ZLSOCPEXN(J)*PFRAC_ICE(J)))
- ZRLTEMP(J)=MIN(MAX(-PRL(J)-PRI(J), ZRLTEMP(J)),PRV(J))
- PRV(J)=PRV(J)-ZRLTEMP(J)
- PRL(J)=PRL(J)+PRI(J)+ZRLTEMP(J)
- PRI(J)=PFRAC_ICE(J) * (PRL(J))
- PRL(J)=(1-PFRAC_ICE(J)) * (PRT(J) - PRV(J))
-
- !Computation of Cph (as defined in Meso-NH doc, equation 2.2, to be used with mixing ratios)
- ZCPH(J)=CST%XCPD+ CST%XCPV * PRV(J)+ CST%XCL * PRL(J) + CST%XCI * PRI(J) + ZCPH2(J)
-
- !Computation of L/Cph/EXN, then new PTH
- ZVAR2=ZCPH(J)*ZEXN(J)
- ZLVOCPEXN(J) = (CST%XLVTT + (CST%XCPV-CST%XCL) * (ZT(J)-CST%XTT)) /ZVAR2
- ZLSOCPEXN(J) = (CST%XLSTT + (CST%XCPV-CST%XCI) * (ZT(J)-CST%XTT)) /ZVAR2
- PTH(J)=PTHL(J)+ZLVOCPEXN(J)*PRL(J)+ZLSOCPEXN(J)*PRI(J)
-
- !Computation of estimated mixing ration at saturation
- !To compute the adjustement a first order development was used
- ZVAR1=PTH(J)*ZEXN(J)-ZT(J)
- PRSATW(J)=PRSATW(J) + ZDRSATODTW(J)*ZVAR1
- PRSATI(J)=PRSATI(J) + ZDRSATODTI(J)*ZVAR1
- ENDDO
-ENDDO
-
-IF (LHOOK) CALL DR_HOOK('TH_R_FROM_THL_RT_1D',1,ZHOOK_HANDLE)
-
-!
-CONTAINS
-INCLUDE "compute_frac_ice.func.h"
-!
-END SUBROUTINE TH_R_FROM_THL_RT_1D
-END MODULE MODE_TH_R_FROM_THL_RT_1D
diff --git a/src/common/turb/mode_th_r_from_thl_rt_2d.F90 b/src/common/turb/mode_th_r_from_thl_rt_2d.F90
deleted file mode 100644
index c4911018a8fba842c010932302d1dd557b3c46e9..0000000000000000000000000000000000000000
--- a/src/common/turb/mode_th_r_from_thl_rt_2d.F90
+++ /dev/null
@@ -1,114 +0,0 @@
-!MNH_LIC Copyright 2006-2022 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 MODE_TH_R_FROM_THL_RT_2D
-IMPLICIT NONE
-CONTAINS
- SUBROUTINE TH_R_FROM_THL_RT_2D(HFRAC_ICE,PFRAC_ICE,PP, &
- PTHL, PRT, PTH, PRV, PRL, PRI, &
- PRSATW, PRSATI, PRR, PRS, PRG, PRH,OOCEAN)
-! #################################################################
-!
-!
-!!**** *TH_R_FROM_THL_RT_2D* - computes the non-conservative variables
-!! from conservative variables
-!!
-!! PURPOSE
-!! -------
-!!
-!!** METHOD
-!! ------
-!!
-!!
-!! EXTERNAL
-!! --------
-!!
-!! IMPLICIT ARGUMENTS
-!! ------------------
-!!
-!!
-!! REFERENCE
-!! ---------
-!!
-!! AUTHOR
-!! ------
-!! Julien PERGAUD * Meteo-France *
-!!
-!! MODIFICATIONS
-!! -------------
-!! Original 13/03/06
-!! Sébastien Riette April 2011: code moved in th_r_from_thl_rt_1D
-!!
-!! --------------------------------------------------------------------------
-!
-!* 0. DECLARATIONS
-! ------------
-!
-!
-USE MODE_TH_R_FROM_THL_RT_3D, ONLY: TH_R_FROM_THL_RT_3D
-USE PARKIND1, ONLY : JPRB
-USE YOMHOOK , ONLY : LHOOK, DR_HOOK
-USE MODD_CST, ONLY: CST
-USE MODD_NEB, ONLY: NEB
-
-IMPLICIT NONE
-!
-!
-!* 0.1 declarations of arguments
-!
-CHARACTER(LEN=1), INTENT(IN) :: HFRAC_ICE
-REAL, DIMENSION(:,:), INTENT(INOUT) :: PFRAC_ICE
-REAL, DIMENSION(:,:), INTENT(IN) :: PP ! Pressure
-REAL, DIMENSION(:,:), INTENT(IN) :: PTHL ! Liquid pot. temp.
-REAL, DIMENSION(:,:), INTENT(IN) :: PRT ! Total mixing ratios
-REAL, DIMENSION(:,:),OPTIONAL,INTENT(IN) :: PRR, PRS, PRG, PRH
-REAL, DIMENSION(:,:), INTENT(OUT):: PTH ! Potential temp.
-REAL, DIMENSION(:,:), INTENT(OUT):: PRV ! vapor mixing ratio
-REAL, DIMENSION(:,:), INTENT(INOUT):: PRL ! cloud mixing ratio
-REAL, DIMENSION(:,:), INTENT(INOUT):: PRI ! ice mixing ratio
-REAL, DIMENSION(:,:), INTENT(OUT) :: PRSATW ! estimated mixing ration at saturation over water
-REAL, DIMENSION(:,:), INTENT(OUT) :: PRSATI ! estimated mixing ration at saturation over ice
-LOGICAL, INTENT(IN) :: OOCEAN ! switch OCEAN version
-
-!
-!-------------------------------------------------------------------------------
-!
-! 0.2 declaration of local variables
-!
-!----------------------------------------------------------------------------
-!
-REAL, DIMENSION(SIZE(PP,1),SIZE(PP,2)) :: ZRR, ZRS, ZRG, ZRH
-REAL(KIND=JPRB) :: ZHOOK_HANDLE
-INTEGER :: JK
-!----------------------------------------------------------------------------
-!
-!* 1 Initialisation
-! --------------
-!
-IF (LHOOK) CALL DR_HOOK('TH_R_FROM_THL_RT_2D',0,ZHOOK_HANDLE)
-ZRR(:,:)=0.
-ZRS(:,:)=0.
-ZRG(:,:)=0.
-ZRH(:,:)=0.
-IF(PRESENT(PRR)) ZRR(:,:)=PRR(:,:)
-IF(PRESENT(PRS)) ZRS(:,:)=PRS(:,:)
-IF(PRESENT(PRG)) ZRG(:,:)=PRG(:,:)
-IF(PRESENT(PRH)) ZRH(:,:)=PRH(:,:)
-!
-!
-! 2 Call of 1d version
-! ------------------
-DO JK=1, SIZE(PTHL,2)
- CALL TH_R_FROM_THL_RT_1D(CST, NEB, SIZE(PTHL,2), HFRAC_ICE,PFRAC_ICE(:,JK),PP(:,JK), &
- PTHL(:,JK), PRT(:,JK), PTH(:,JK), &
- PRV(:,JK), PRL(:,JK), PRI(:,JK), &
- PRSATW(:,JK), PRSATI(:,JK), &
- ZRR(:,JK), ZRS(:,JK), ZRG(:,JK), ZRH(:,JK),OOCEAN)
-ENDDO
-
-IF (LHOOK) CALL DR_HOOK('TH_R_FROM_THL_RT_2D',1,ZHOOK_HANDLE)
-
-END SUBROUTINE TH_R_FROM_THL_RT_2D
-END MODULE MODE_TH_R_FROM_THL_RT_2D
-
diff --git a/src/common/turb/mode_th_r_from_thl_rt_3d.F90 b/src/common/turb/mode_th_r_from_thl_rt_3d.F90
deleted file mode 100644
index 2c1733ce9dda37de77a20a615c06feb43f40999b..0000000000000000000000000000000000000000
--- a/src/common/turb/mode_th_r_from_thl_rt_3d.F90
+++ /dev/null
@@ -1,110 +0,0 @@
-!MNH_LIC Copyright 2006-2022 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 MODE_TH_R_FROM_THL_RT_3D
-IMPLICIT NONE
-CONTAINS
- SUBROUTINE TH_R_FROM_THL_RT_3D(HFRAC_ICE,PFRAC_ICE,PP, &
- PTHL, PRT, PTH, PRV, PRL, PRI, &
- PRSATW, PRSATI, PRR, PRS, PRG, PRH,OOCEAN)
-! #################################################################
-!
-!
-!!**** *TH_R_FROM_THL_RT_3D* - computes the non-conservative variables
-!! from conservative variables
-!!
-!! PURPOSE
-!! -------
-!!
-!!** METHOD
-!! ------
-!!
-!!
-!! EXTERNAL
-!! --------
-!!
-!! IMPLICIT ARGUMENTS
-!! ------------------
-!!
-!!
-!! REFERENCE
-!! ---------
-!!
-!! AUTHOR
-!! ------
-!! Julien PERGAUD * Meteo-France *
-!!
-!! MODIFICATIONS
-!! -------------
-!! Original 15/03/2011
-!! S. Riette April 2011 : code moved in th_r_from_thl_rt_1d
-!!
-!! --------------------------------------------------------------------------
-!
-!* 0. DECLARATIONS
-! ------------
-!
-USE MODE_TH_R_FROM_THL_RT_1D, ONLY: TH_R_FROM_THL_RT_1D
-USE PARKIND1, ONLY : JPRB
-USE YOMHOOK , ONLY : LHOOK, DR_HOOK
-USE MODD_CST, ONLY : CST
-USE MODD_NEB, ONLY : NEB
-!
-IMPLICIT NONE
-!
-!
-!* 0.1 declarations of arguments
-!
-CHARACTER(LEN=1), INTENT(IN) :: HFRAC_ICE
-REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PFRAC_ICE
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PP ! Pressure
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHL ! thetal to transform into th
-REAL, DIMENSION(:,:,:),INTENT(IN) :: PRT ! Total mixing ratios to transform into rv,rc and ri
-REAL, DIMENSION(:,:,:),OPTIONAL,INTENT(IN) :: PRR, PRS, PRG, PRH
-REAL, DIMENSION(:,:,:), INTENT(OUT):: PTH ! th
-REAL, DIMENSION(:,:,:), INTENT(OUT):: PRV ! vapor mixing ratio
-REAL, DIMENSION(:,:,:), INTENT(INOUT):: PRL ! vapor mixing ratio
-REAL, DIMENSION(:,:,:), INTENT(INOUT):: PRI ! vapor mixing ratio
-REAL, DIMENSION(:,:,:), INTENT(OUT) :: PRSATW ! estimated mixing ration at saturation over water
-REAL, DIMENSION(:,:,:), INTENT(OUT) :: PRSATI ! estimated mixing ration at saturation over ice
-LOGICAL, INTENT(IN) :: OOCEAN ! switch OCEAN version
-!
-!-------------------------------------------------------------------------------
-!
-! 0.2 declaration of local variables
-REAL, DIMENSION(SIZE(PTHL,1),SIZE(PTHL,2),SIZE(PTHL,3)) :: ZRR, ZRS, ZRG, ZRH
-REAL(KIND=JPRB) :: ZHOOK_HANDLE
-INTEGER :: JJ, JK
-!----------------------------------------------------------------------------
-!
-!* 1 Initialisation
-! --------------
-!
-IF (LHOOK) CALL DR_HOOK('TH_R_FROM_THL_RT_3D',0,ZHOOK_HANDLE)
-ZRR(:,:,:)=0.
-ZRS(:,:,:)=0.
-ZRG(:,:,:)=0.
-ZRH(:,:,:)=0.
-IF(PRESENT(PRR)) ZRR(:,:,:)=PRR(:,:,:)
-IF(PRESENT(PRS)) ZRS(:,:,:)=PRS(:,:,:)
-IF(PRESENT(PRG)) ZRG(:,:,:)=PRG(:,:,:)
-IF(PRESENT(PRH)) ZRH(:,:,:)=PRH(:,:,:)
-!
-!
-! 2 Call of 1d version
-! ------------------
-DO JK=1, SIZE(PTHL,3)
- DO JJ=1, SIZE(PTHL,2)
- CALL TH_R_FROM_THL_RT_1D(CST, NEB, SIZE(PTHL,2), HFRAC_ICE,PFRAC_ICE(:,JJ,JK),PP(:,JJ,JK), &
- PTHL(:,JJ,JK), PRT(:,JJ,JK), PTH(:,JJ,JK), &
- PRV(:,JJ,JK), PRL(:,JJ,JK), PRI(:,JJ,JK), &
- PRSATW(:,JJ,JK), PRSATI(:,JJ,JK), &
- ZRR(:,JJ,JK), ZRS(:,JJ,JK), ZRG(:,JJ,JK), ZRH(:,JJ,JK),OOCEAN)
- ENDDO
-ENDDO
-
-IF (LHOOK) CALL DR_HOOK('TH_R_FROM_THL_RT_3D',1,ZHOOK_HANDLE)
-
-END SUBROUTINE TH_R_FROM_THL_RT_3D
-END MODULE MODE_TH_R_FROM_THL_RT_3D
diff --git a/src/common/turb/th_r_from_thl_rt.func.h b/src/common/turb/th_r_from_thl_rt.func.h
new file mode 100644
index 0000000000000000000000000000000000000000..62cbcd5c79fd7af3ac78ffaf24ed9c13783a754c
--- /dev/null
+++ b/src/common/turb/th_r_from_thl_rt.func.h
@@ -0,0 +1,197 @@
+!MNH_LIC Copyright 2006-2022 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.
+ SUBROUTINE TH_R_FROM_THL_RT(CST, NEB, KT, HFRAC_ICE,PFRAC_ICE,PP, &
+ PTHL, PRT, PTH, PRV, PRL, PRI, &
+ PRSATW, PRSATI, PRR, PRS, PRG, PRH, OOCEAN,&
+ PBUF)
+! ******* TO BE INCLUDED IN THE *CONTAINS* OF A SUBROUTINE, IN ORDER TO EASE AUTOMATIC INLINING ******
+! => Don't use drHook !!!
+! "compute_frac_ice.func.h" must be included at the same time
+! #################################################################
+!
+!
+!!**** *TH_R_FROM_THL_RT* - computes the non-conservative variables
+!! from conservative variables
+!!
+!! PURPOSE
+!! -------
+!!
+!!** METHOD
+!! ------
+!!
+!!
+!! EXTERNAL
+!! --------
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!!
+!!
+!! REFERENCE
+!! ---------
+!!
+!! AUTHOR
+!! ------
+!! Julien PERGAUD * Meteo-France *
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 13/03/06
+!! S. Riette April 2011 : ice added, allow ZRLTEMP to be negative
+!! we use dQsat/dT to help convergence
+!! use of optional PRR, PRS, PRG, PRH
+!! S. Riette Nov 2016: support for HFRAC_ICE='S'
+!!
+!! --------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODD_CST, ONLY : CST_t
+USE MODD_NEB, ONLY : NEB_t
+!
+IMPLICIT NONE
+!
+!
+!* 0.1 declarations of arguments
+!
+TYPE(CST_t), INTENT(IN) :: CST
+TYPE(NEB_t), INTENT(IN) :: NEB
+INTEGER, INTENT(IN) :: KT
+CHARACTER(LEN=1), INTENT(IN) :: HFRAC_ICE
+LOGICAL, INTENT(IN) :: OOCEAN ! switch for Ocean model version
+REAL, DIMENSION(KT), INTENT(INOUT) :: PFRAC_ICE
+REAL, DIMENSION(KT), INTENT(IN) :: PP ! Pressure
+REAL, DIMENSION(KT), INTENT(IN) :: PTHL ! thetal to transform into th
+REAL, DIMENSION(KT),INTENT(IN) :: PRT ! Total mixing ratios to transform into rv,rc and ri
+REAL, DIMENSION(KT),OPTIONAL,INTENT(IN) :: PRR, PRS, PRG, PRH
+REAL, DIMENSION(KT), INTENT(OUT):: PTH ! th
+REAL, DIMENSION(KT), INTENT(OUT):: PRV ! vapor mixing ratio
+REAL, DIMENSION(KT), INTENT(INOUT):: PRL ! vapor mixing ratio
+REAL, DIMENSION(KT), INTENT(INOUT):: PRI ! vapor mixing ratio
+REAL, DIMENSION(KT), INTENT(OUT) :: PRSATW ! estimated mixing ration at saturation over water
+REAL, DIMENSION(KT), INTENT(OUT) :: PRSATI ! estimated mixing ration at saturation over ice
+REAL, DIMENSION(KT, 16), INTENT(OUT) :: PBUF ! buffer to replace automatic arrays
+!
+!-------------------------------------------------------------------------------
+!
+! 0.2 declaration of local variables
+INTEGER :: II ! Loop control
+INTEGER :: JITER ! number of iterations
+INTEGER :: J
+INTEGER, PARAMETER :: IEXN=1, IRVSAT=2, ICPH=3, IRLTEMP=4, ICPH2=5, IT=6, ILVOCPEXN=7, ILSOCPEXN=8, &
+ & IDRSATODT=9, IDRSATODTW=10, IDRSATODTI=11, IFOESW=12, IFOESI=13, &
+ & ILOGT=14, I99PP=15, I1PRT=16
+REAL :: ZVAR1, ZVAR2, ZTPOW2, ZDELT
+
+!----------------------------------------------------------------------------
+!
+!* 1 Initialisation
+! --------------
+!
+!
+!
+!Number of iterations
+JITER=2
+!
+!Computation of PBUF(:, ICPH2) depending on dummy arguments received
+PBUF(:, ICPH2)=0
+IF(PRESENT(PRR)) PBUF(:, ICPH2)=PBUF(:, ICPH2) + CST%XCL*PRR(:)
+IF(PRESENT(PRS)) PBUF(:, ICPH2)=PBUF(:, ICPH2) + CST%XCI*PRS(:)
+IF(PRESENT(PRG)) PBUF(:, ICPH2)=PBUF(:, ICPH2) + CST%XCI*PRG(:)
+IF(PRESENT(PRH)) PBUF(:, ICPH2)=PBUF(:, ICPH2) + CST%XCI*PRH(:)
+!
+!Computation of an approximate state thanks to PRL and PRI guess
+PBUF(:, IEXN)=(PP(:)/CST%XP00) ** CST%RDSCPD
+
+DO J=1,KT
+ PBUF(J, I99PP)=0.99*PP(J)
+ PRV(J)=PRT(J)-PRL(J)-PRI(J)
+ PBUF(J, ICPH)=CST%XCPD+ CST%XCPV * PRV(J)+ CST%XCL * PRL(J) + CST%XCI * PRI(J) + PBUF(J, ICPH2)
+ ZVAR2=PBUF(J, ICPH)*PBUF(J, IEXN)
+ ZDELT=(PTHL(J)*PBUF(J, IEXN))-CST%XTT
+ PBUF(J, ILVOCPEXN) = (CST%XLVTT + (CST%XCPV-CST%XCL) * ZDELT) /ZVAR2
+ PBUF(J, ILSOCPEXN) = (CST%XLSTT + (CST%XCPV-CST%XCI) * ZDELT) /ZVAR2
+ PTH(J)=PTHL(J)+PBUF(J, ILVOCPEXN)*PRL(J)+PBUF(J, ILSOCPEXN)*PRI(J)
+ PBUF(J, I1PRT)=1+PRT(J)
+ENDDO
+!
+!
+! 2 Iteration
+! ---------
+
+DO II=1,JITER
+ IF (OOCEAN) THEN
+ PBUF(:, IT)=PTH(:)
+ ELSE
+ PBUF(:, IT)=PTH(:)*PBUF(:, IEXN)
+ END IF
+ !Computation of liquid/ice fractions
+ PFRAC_ICE(:) = 0.
+ DO J=1, KT
+ IF(PRL(J)+PRI(J) > 1.E-20) THEN
+ PFRAC_ICE(J) = PRI(J) / (PRL(J)+PRI(J))
+ ENDIF
+ ENDDO
+ CALL COMPUTE_FRAC_ICE(HFRAC_ICE,NEB,PFRAC_ICE(:),PBUF(:, IT))
+
+ !Computation of Rvsat and dRsat/dT
+ !In this version QSAT, QSATI, DQSAT and DQASATI functions are not used
+ !due to performance issue
+
+ ! Log does not vectorize on all compilers:
+ PBUF(:, ILOGT)=LOG(PBUF(:, IT))
+
+ DO J=1,KT
+ PBUF(J, IFOESW) = MIN(EXP( CST%XALPW - CST%XBETAW/PBUF(J, IT) - CST%XGAMW*PBUF(J, ILOGT) ), PBUF(J, I99PP))
+ PBUF(J, IFOESI) = MIN(EXP( CST%XALPI - CST%XBETAI/PBUF(J, IT) - CST%XGAMI*PBUF(J, ILOGT) ), PBUF(J, I99PP))
+ PRSATW(J) = CST%XRD/CST%XRV*PBUF(J, IFOESW)/PP(J) / (1.+(CST%XRD/CST%XRV-1.)*PBUF(J, IFOESW)/PP(J))
+ PRSATI(J) = CST%XRD/CST%XRV*PBUF(J, IFOESI)/PP(J) / (1.+(CST%XRD/CST%XRV-1.)*PBUF(J, IFOESI)/PP(J))
+ ZTPOW2=PBUF(J, IT)**2
+ PBUF(J, IDRSATODTW) = PRSATW(J) / (1.+(CST%XRD/CST%XRV-1.)*PBUF(J, IFOESW)/PP(J) ) &
+ * (CST%XBETAW/ZTPOW2 - CST%XGAMW/PBUF(J, IT))*PBUF(J, I1PRT)
+ PBUF(J, IDRSATODTI) = PRSATI(J) / (1.+(CST%XRD/CST%XRV-1.)*PBUF(J, IFOESI)/PP(J) ) &
+ * (CST%XBETAI/ZTPOW2 - CST%XGAMI/PBUF(J, IT))*PBUF(J, I1PRT)
+ !PRSATW(J) = QSAT(PBUF(J, IT),PP(J)) !qsatw
+ !PRSATI(J) = QSATI(PBUF(J, IT),PP(J)) !qsati
+ !PBUF(J, IDRSATODTW) = DQSAT(PBUF(J, IT),PP(J),PRSATW(J))*PBUF(J, I1PRT)
+ !PBUF(J, IDRSATODTI) = DQSATI(PBUF(J, IT),PP(J),PRSATI(J))*PBUF(J, I1PRT)
+ PRSATW(J) = PRSATW(J)*PBUF(J, I1PRT)
+ PRSATI(J) = PRSATI(J)*PBUF(J, I1PRT)
+ PBUF(J, IRVSAT) = PRSATW(J)*(1-PFRAC_ICE(J)) + PRSATI(J)*PFRAC_ICE(J)
+ PBUF(J, IDRSATODT) = (PBUF(J, IDRSATODTW)*(1-PFRAC_ICE(J))+ &
+ & PBUF(J, IDRSATODTI)*PFRAC_ICE(J))
+
+ !Computation of new PRL, PRI and PRV
+ !Correction term applied to (PRV(J)-PBUF(J, IRVSAT)) is computed assuming that
+ !PBUF(J, ILVOCPEXN), PBUF(J, ILSOCPEXN) and PBUF(J, ICPH) don't vary to much with T. It takes into account
+ !the variation (estimated linear) of Qsat with T
+ PBUF(J, IRLTEMP)=(PRV(J)-PBUF(J, IRVSAT))/ &
+ &(1 + PBUF(J, IDRSATODT)*PBUF(J, IEXN)* &
+ & (PBUF(J, ILVOCPEXN)*(1-PFRAC_ICE(J))+PBUF(J, ILSOCPEXN)*PFRAC_ICE(J)))
+ PBUF(J, IRLTEMP)=MIN(MAX(-PRL(J)-PRI(J), PBUF(J, IRLTEMP)),PRV(J))
+ PRV(J)=PRV(J)-PBUF(J, IRLTEMP)
+ PRL(J)=PRL(J)+PRI(J)+PBUF(J, IRLTEMP)
+ PRI(J)=PFRAC_ICE(J) * (PRL(J))
+ PRL(J)=(1-PFRAC_ICE(J)) * (PRT(J) - PRV(J))
+
+ !Computation of Cph (as defined in Meso-NH doc, equation 2.2, to be used with mixing ratios)
+ PBUF(J, ICPH)=CST%XCPD+ CST%XCPV * PRV(J)+ CST%XCL * PRL(J) + CST%XCI * PRI(J) + PBUF(J, ICPH2)
+
+ !Computation of L/Cph/EXN, then new PTH
+ ZVAR2=PBUF(J, ICPH)*PBUF(J, IEXN)
+ PBUF(J, ILVOCPEXN) = (CST%XLVTT + (CST%XCPV-CST%XCL) * (PBUF(J, IT)-CST%XTT)) /ZVAR2
+ PBUF(J, ILSOCPEXN) = (CST%XLSTT + (CST%XCPV-CST%XCI) * (PBUF(J, IT)-CST%XTT)) /ZVAR2
+ PTH(J)=PTHL(J)+PBUF(J, ILVOCPEXN)*PRL(J)+PBUF(J, ILSOCPEXN)*PRI(J)
+
+ !Computation of estimated mixing ration at saturation
+ !To compute the adjustement a first order development was used
+ ZVAR1=PTH(J)*PBUF(J, IEXN)-PBUF(J, IT)
+ PRSATW(J)=PRSATW(J) + PBUF(J, IDRSATODTW)*ZVAR1
+ PRSATI(J)=PRSATI(J) + PBUF(J, IDRSATODTI)*ZVAR1
+ ENDDO
+ENDDO
+
+END SUBROUTINE TH_R_FROM_THL_RT
diff --git a/src/mesonh/ext/ice_adjust_bis.f90 b/src/mesonh/ext/ice_adjust_bis.f90
index c6f72a0868eee72ac6bd97ad98e4d3bd692529b2..44ab0c680b6d689ab050c53ddd39ec799bf0b100 100644
--- a/src/mesonh/ext/ice_adjust_bis.f90
+++ b/src/mesonh/ext/ice_adjust_bis.f90
@@ -65,10 +65,10 @@ END MODULE MODI_ICE_ADJUST_BIS
!* 0. DECLARATIONS
! ------------
!
-USE MODD_CST, ONLY : XCPD, XRD, XP00
+USE MODD_CST, ONLY : XCPD, XRD, XP00, CST
+USE MODD_NEB, ONLY : NEB
!
USE MODI_COMPUTE_FUNCTION_THERMO
-USE MODE_TH_R_FROM_THL_RT_3D
USE MODI_THLRT_FROM_THRVRCRI
!
USE MODE_ll
@@ -89,6 +89,7 @@ REAL, DIMENSION(SIZE(PTH,1),SIZE(PTH,2),SIZE(PTH,3)) :: ZTHL, ZRW, ZRV, ZRC, &
ZRI, ZWORK
REAL, DIMENSION(SIZE(PTH,1),SIZE(PTH,2),SIZE(PTH,3)) :: ZFRAC_ICE, ZRSATW, ZRSATI
REAL, DIMENSION(SIZE(PTH,1),SIZE(PTH,2),SIZE(PTH,3)) :: ZT, ZEXN, ZLVOCPEXN,ZLSOCPEXN
+REAL, DIMENSION(SIZE(PTH,1),SIZE(PTH,2),SIZE(PTH,3), 16) :: ZBUF
INTEGER :: IRR
CHARACTER(LEN=1) :: YFRAC_ICE
!
@@ -127,10 +128,11 @@ CALL COMPUTE_FUNCTION_THERMO( IRR, &
CALL THLRT_FROM_THRVRCRI( IRR, PTH, PR, ZLVOCPEXN, ZLSOCPEXN,&
ZTHL, ZRW )
!
-CALL TH_R_FROM_THL_RT_3D(YFRAC_ICE,ZFRAC_ICE(:,:,:),PP(:,:,:), &
+CALL TH_R_FROM_THL_RT(CST, NEB, SIZE(ZFRAC_ICE), YFRAC_ICE,ZFRAC_ICE(:,:,:),PP(:,:,:), &
ZTHL(:,:,:), ZRW(:,:,:), PTH(:,:,:), &
ZRV(:,:,:), ZRC(:,:,:), ZRI(:,:,:), &
- ZRSATW(:,:,:), ZRSATI(:,:,:),OOCEAN=.FALSE.)
+ ZRSATW(:,:,:), ZRSATI(:,:,:),OOCEAN=.FALSE.,&
+ PBUF=ZBUF)
CALL ADD3DFIELD_ll( TZFIELDS_ll, PTH, 'ICE_ADJUST_BIS::PTH')
IF (IRR>=1) THEN
CALL ADD3DFIELD_ll( TZFIELDS_ll, ZRV, 'ICE_ADJUST_BIS::ZRV' )
@@ -152,4 +154,7 @@ PR(:,:,:,2) = ZRC(:,:,:)
IF (IRR>=4) &
PR(:,:,:,4) = ZRI(:,:,:)
!
+CONTAINS
+INCLUDE "th_r_from_thl_rt.func.h"
+INCLUDE "compute_frac_ice.func.h"
END SUBROUTINE ICE_ADJUST_BIS
diff --git a/src/mesonh/ext/prep_ideal_case.f90 b/src/mesonh/ext/prep_ideal_case.f90
index a1c1ec6c6a4d7ab52ea4549f93ffc4a35a2deb5b..7a1b8787e20d1a285ca8977ac578783f49181674 100644
--- a/src/mesonh/ext/prep_ideal_case.f90
+++ b/src/mesonh/ext/prep_ideal_case.f90
@@ -416,7 +416,6 @@ USE MODD_SUB_MODEL_n
USE MODE_MNH_TIMING
USE MODN_CONFZ
!JUAN
-USE MODE_TH_R_FROM_THL_RT_3D
!
USE MODI_VERSION
USE MODI_INIT_PGD_SURF_ATM
@@ -431,6 +430,7 @@ USE MODI_SET_RELFRC
!
USE MODI_INI_CST
USE MODI_INI_NEB
+USE MODD_NEB, ONLY: NEB
USE MODI_WRITE_HGRID
USE MODD_MPIF
USE MODD_VAR_ll
@@ -556,6 +556,7 @@ REAL, DIMENSION(:), ALLOCATABLE :: ZXHAT_ll, ZYHAT_ll
REAL, DIMENSION(:,:,:), ALLOCATABLE ::ZTHL,ZT,ZRT,ZFRAC_ICE,&
ZEXN,ZLVOCPEXN,ZLSOCPEXN,ZCPH, &
ZRSATW, ZRSATI
+REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: ZBUF
! variables for adjustement
REAL :: ZDIST
!
@@ -1667,6 +1668,7 @@ IF (CIDEAL == 'RSOU') THEN
ALLOCATE(ZFRAC_ICE(NIU,NJU,NKU))
ALLOCATE(ZRSATW(NIU,NJU,NKU))
ALLOCATE(ZRSATI(NIU,NJU,NKU))
+ ALLOCATE(ZBUF(NIU,NJU,NKU,16))
ZRT=XRT(:,:,:,1)+XRT(:,:,:,2)+XRT(:,:,:,4)
IF (LOCEAN) THEN
ZEXN(:,:,:)= 1.
@@ -1682,8 +1684,9 @@ ELSE
ZLVOCPEXN = (XLVTT + (XCPV-XCL) * (ZT-XTT))/(ZCPH*ZEXN)
ZLSOCPEXN = (XLSTT + (XCPV-XCI) * (ZT-XTT))/(ZCPH*ZEXN)
ZTHL=XTHT-ZLVOCPEXN*XRT(:,:,:,2)-ZLSOCPEXN*XRT(:,:,:,4)
- CALL TH_R_FROM_THL_RT_3D('T',ZFRAC_ICE,XPABST,ZTHL,ZRT,XTHT,XRT(:,:,:,1), &
- XRT(:,:,:,2),XRT(:,:,:,4),ZRSATW, ZRSATI,OOCEAN=.FALSE.)
+ CALL TH_R_FROM_THL_RT(CST, NEB, SIZE(ZFRAC_ICE), 'T',ZFRAC_ICE,XPABST,ZTHL,ZRT,XTHT,XRT(:,:,:,1), &
+ XRT(:,:,:,2),XRT(:,:,:,4),ZRSATW, ZRSATI,OOCEAN=.FALSE.,&
+ PBUF=ZBUF)
END IF
DEALLOCATE(ZEXN)
DEALLOCATE(ZT)
@@ -1692,6 +1695,7 @@ END IF
DEALLOCATE(ZLSOCPEXN)
DEALLOCATE(ZTHL)
DEALLOCATE(ZRT)
+ DEALLOCATE(ZBUF)
! Coherence test
IF ((.NOT. LUSERI) ) THEN
IF (MAXVAL(XRT(:,:,:,4))/= 0) THEN
@@ -1930,4 +1934,8 @@ WRITE(NLUOUT,FMT=*) '****************************************************'
!
CALL FINALIZE_MNH()
!
+!
+CONTAINS
+INCLUDE "th_r_from_thl_rt.func.h"
+INCLUDE "compute_frac_ice.func.h"
END PROGRAM PREP_IDEAL_CASE
diff --git a/src/mesonh/ext/set_rsou.f90 b/src/mesonh/ext/set_rsou.f90
index 8943cfbcde9ada76f156fb26f2a44b1a176ad123..352af8a53b1efdaca9aa28ef28c9658ef2d27ef5 100644
--- a/src/mesonh/ext/set_rsou.f90
+++ b/src/mesonh/ext/set_rsou.f90
@@ -283,7 +283,6 @@ USE MODI_PRESS_HEIGHT
USE MODI_SET_MASS
USE MODI_SHUMAN
USE MODI_THETAVPU_THETAVPM
-USE MODE_TH_R_FROM_THL_RT_1D
USE MODI_VERT_COORD
!
USE NETCDF ! for reading the NR files
@@ -375,6 +374,7 @@ REAL, DIMENSION(:), ALLOCATABLE :: ZFRAC_ICE ! ice fraction
REAL, DIMENSION(:), ALLOCATABLE :: ZRSATW, ZRSATI
REAL :: ZDZSDH,ZDZ1SDH,ZDZ2SDH ! interpolation
! working arrays
+REAL, DIMENSION(:,:), ALLOCATABLE :: ZBUF
!
INTEGER :: JK,JKLEV,JKU,JKM,JKT,JJ,JI,JO,JLOOP ! Loop indexes
INTEGER :: IKU ! Upper bound in z direction
@@ -1574,6 +1574,7 @@ ALLOCATE(ZFRAC_ICE(IKU))
ALLOCATE(ZRSATW(IKU))
ALLOCATE(ZRSATI(IKU))
ALLOCATE(ZMRT(IKU))
+ALLOCATE(ZBUF(IKU,16))
ZMRT=ZMRM+ZMRCM+ZMRIM
ZTHVM=ZTHLM
!
@@ -1593,8 +1594,9 @@ ELSE
DO JLOOP=1,20 ! loop for pression
CALL COMPUTE_EXNER_FROM_GROUND(ZTHVM,ZZMASS_PROFILE(:),ZEXNSURF,ZEXNFLUX,ZEXNMASS)
ZPRESS(:)=XP00*(ZEXNMASS(:))**(XCPD/XRD)
- CALL TH_R_FROM_THL_RT_1D(CST,NEB,SIZE(ZPRESS,1),'T',ZFRAC_ICE,ZPRESS,ZTHLM,ZMRT,ZTHM,ZMRM,ZMRCM,ZMRIM, &
- ZRSATW, ZRSATI,OOCEAN=.FALSE.)
+ CALL TH_R_FROM_THL_RT(CST,NEB,SIZE(ZPRESS,1),'T',ZFRAC_ICE,ZPRESS,ZTHLM,ZMRT,ZTHM,ZMRM,ZMRCM,ZMRIM, &
+ ZRSATW, ZRSATI,OOCEAN=.FALSE.,&
+ PBUF=ZBUF)
ZTHVM(:)=ZTHM(:)*(1.+XRV/XRD*ZMRM(:))/(1.+(ZMRM(:)+ZMRIM(:)+ZMRCM(:)))
ENDDO
ENDIF
@@ -1606,6 +1608,7 @@ DEALLOCATE(ZFRAC_ICE)
DEALLOCATE(ZRSATW)
DEALLOCATE(ZRSATI)
DEALLOCATE(ZMRT)
+DEALLOCATE(ZBUF)
!-------------------------------------------------------------------------------
!
!* 4. COMPUTE FIELDS ON THE MODEL GRID (WITH OROGRAPHY)
@@ -1630,5 +1633,8 @@ CONTAINS
CALL PRINT_MSG( NVERB_ERROR, 'IO', 'SET_RSOU', 'error at ' // Trim( yloc) // ': ' // NF90_STRERROR( ISTATUS ) )
END IF
END SUBROUTINE check
-!
+ !
+ INCLUDE "th_r_from_thl_rt.func.h"
+ INCLUDE "compute_frac_ice.func.h"
+ !
END SUBROUTINE SET_RSOU