From 73ec9d536204688230716131907da28ff9352c48 Mon Sep 17 00:00:00 2001
From: Quentin Rodier <quentin.rodier@meteo.fr>
Date: Fri, 25 Nov 2022 17:43:57 +0100
Subject: [PATCH] Quentin 25/11/2022: merging from MesoNH dev 5.5 to 5.6 (first
part).
---
src/arome/ext/aro_turb_mnh.F90 | 6 +-
src/common/aux/modd_cst.F90 | 2 +-
src/common/aux/modd_nsv.F90 | 15 ++
src/common/micro/modd_param_ice.F90 | 5 +-
src/common/micro/modd_rain_ice_descr.F90 | 16 +-
src/common/micro/mode_ice4_fast_rg.F90 | 6 +
src/common/micro/mode_ice4_fast_rh.F90 | 6 +
src/common/micro/mode_ice4_fast_rs.F90 | 58 ++++-
src/common/micro/mode_ice4_rsrimcg_old.F90 | 4 +
.../micro/mode_ice4_sedimentation_split.F90 | 74 +++++--
.../micro/mode_ice4_sedimentation_stat.F90 | 3 +
src/common/micro/mode_ice4_slow.F90 | 14 ++
src/common/micro/rain_ice.F90 | 38 +++-
src/common/turb/modi_turb.F90 | 2 +
src/common/turb/turb.F90 | 24 +++
src/mesonh/aux/sources_neg_correct.f90 | 152 +++++++++-----
src/mesonh/ext/deallocate_model1.f90 | 145 +++++++++++++
src/mesonh/ext/ground_paramn.f90 | 198 +++++++++++++++++-
src/mesonh/ext/ini_radar.f90 | 1 +
src/mesonh/ext/lesn.f90 | 7 +
src/mesonh/ext/phys_paramn.f90 | 14 +-
src/mesonh/ext/prep_ideal_case.f90 | 8 +
22 files changed, 694 insertions(+), 104 deletions(-)
diff --git a/src/arome/ext/aro_turb_mnh.F90 b/src/arome/ext/aro_turb_mnh.F90
index 8a0bf7851..61c0f8b85 100644
--- a/src/arome/ext/aro_turb_mnh.F90
+++ b/src/arome/ext/aro_turb_mnh.F90
@@ -69,6 +69,7 @@
! ------------
!
USE MODD_CONF
+USE MODD_NSV, ONLY: NSV_LIMA_NR, NSV_LIMA_NS, NSV_LIMA_NG, NSV_LIMA_NH
USE MODD_CST, ONLY:CST
USE MODD_CTURB, ONLY:CSTURB
USE MODD_LES, ONLY:LLES_CALL
@@ -425,8 +426,9 @@ ENDDO
OCOMPUTE_SRC=SIZE(PSIGS, 3)/=0
CALL TURB (CST,CSTURB,TBUCONF,TURBN, YLDIMPHYEX,&
& IMI, KRR, KRRL, KRRI, HLBCX, HLBCY,&
- & ISPLIT,IMI, KSV, KSV_LGBEG, KSV_LGEND, &
- & HPROGRAM, O2D, ONOMIXLG, OFLAT, LLES_CALL,OCOUPLES,OBLOWSNOW,&
+ & ISPLIT,IMI, KSV, KSV_LGBEG, KSV_LGEND, HPROGRAM,&
+ & NSV_LIMA_NR, NSV_LIMA_NS, NSV_LIMA_NG, NSV_LIMA_NH, &
+ & O2D, ONOMIXLG, OFLAT, LLES_CALL,OCOUPLES,OBLOWSNOW,&
& OCOMPUTE_SRC, 1.0, &
& OOCEAN,ODEEPOC, .FALSE., &
& 'NONE',CMICRO, &
diff --git a/src/common/aux/modd_cst.F90 b/src/common/aux/modd_cst.F90
index 7e647e4e6..544c75464 100644
--- a/src/common/aux/modd_cst.F90
+++ b/src/common/aux/modd_cst.F90
@@ -97,7 +97,7 @@ REAL :: XD2=23.
!REAL :: XD1=0.35
!REAL :: XD2=23.
-REAL :: XRHOLI ! Volumic mass of liquid water
+REAL :: XRHOLI ! Volumic mass of ice
!
INTEGER :: NDAYSEC ! Number of seconds in a day
!
diff --git a/src/common/aux/modd_nsv.F90 b/src/common/aux/modd_nsv.F90
index 63cab9dbf..ca85a9433 100644
--- a/src/common/aux/modd_nsv.F90
+++ b/src/common/aux/modd_nsv.F90
@@ -31,6 +31,7 @@
!! V. Vionnet 07/17 add blowing snow
! P. Wautelet 10/03/2021: add CSVNAMES and CSVNAMES_A to store the name of all the scalar variables
! B. Vie 06/2021: add prognostic supersaturation for LIMA
+! A. Costes 12/2021: add Blaze fire model smoke
!
!-------------------------------------------------------------------------------
!
@@ -132,6 +133,9 @@ INTEGER,DIMENSION(JPMODELMAX)::NSV_LIMA_CCN_FREE_A = 0 ! First Free CCN conc.
INTEGER,DIMENSION(JPMODELMAX)::NSV_LIMA_CCN_ACTI_A = 0 ! First Acti. CNN conc.
INTEGER,DIMENSION(JPMODELMAX)::NSV_LIMA_SCAVMASS_A = 0 ! Scavenged mass variable
INTEGER,DIMENSION(JPMODELMAX)::NSV_LIMA_NI_A = 0 ! First Ni var.
+INTEGER,DIMENSION(JPMODELMAX)::NSV_LIMA_NS_A = 0 ! First Ns var.
+INTEGER,DIMENSION(JPMODELMAX)::NSV_LIMA_NG_A = 0 ! First Ng var.
+INTEGER,DIMENSION(JPMODELMAX)::NSV_LIMA_NH_A = 0 ! First Nh var.
INTEGER,DIMENSION(JPMODELMAX)::NSV_LIMA_IFN_FREE_A = 0 ! First Free IFN conc.
INTEGER,DIMENSION(JPMODELMAX)::NSV_LIMA_IFN_NUCL_A = 0 ! First Nucl. IFN conc.
INTEGER,DIMENSION(JPMODELMAX)::NSV_LIMA_IMM_NUCL_A = 0 ! First Nucl. IMM conc.
@@ -145,6 +149,10 @@ INTEGER,DIMENSION(JPMODELMAX)::NSV_FFEND_A = 0 ! NSV_FFBEG_A...NSV_FFEND_A
#endif
!
INTEGER,DIMENSION(JPMODELMAX)::NSV_CO2_A = 0 ! index for CO2
+! Blaze smoke indexes
+INTEGER,DIMENSION(JPMODELMAX)::NSV_FIRE_A = 0 ! number of Blaze smoke scalar variables
+INTEGER,DIMENSION(JPMODELMAX)::NSV_FIREBEG_A = 0 ! with indices in the range :
+INTEGER,DIMENSION(JPMODELMAX)::NSV_FIREEND_A = 0 ! NSV_FIREBEG_A...NSV_FIREEND_A
!
INTEGER,DIMENSION(JPMODELMAX)::NSV_SNW_A = 0 ! number of blowing snow scalar
INTEGER,DIMENSION(JPMODELMAX)::NSV_SNWBEG_A = 0 ! with indices in the range :
@@ -237,6 +245,9 @@ INTEGER :: NSV_LIMA_CCN_FREE !
INTEGER :: NSV_LIMA_CCN_ACTI !
INTEGER :: NSV_LIMA_SCAVMASS !
INTEGER :: NSV_LIMA_NI !
+INTEGER :: NSV_LIMA_NS !
+INTEGER :: NSV_LIMA_NG !
+INTEGER :: NSV_LIMA_NH !
INTEGER :: NSV_LIMA_IFN_FREE !
INTEGER :: NSV_LIMA_IFN_NUCL !
INTEGER :: NSV_LIMA_IMM_NUCL !
@@ -250,6 +261,10 @@ INTEGER :: NSV_FFEND = 0 ! NSV_FFBEG...NSV_FFEND
#endif
!
INTEGER :: NSV_CO2 = 0 ! index for CO2
+! Blaze smoke
+INTEGER :: NSV_FIRE = 0 ! number of Blaze smoke scalar variables
+INTEGER :: NSV_FIREBEG = 0 ! with indices in the range :
+INTEGER :: NSV_FIREEND = 0 ! NSV_FIREBEG...NSV_FIREEND
!
INTEGER :: NSV_SNW = 0 ! number of blowing snow scalar variables
INTEGER :: NSV_SNWBEG = 0 ! with indices in the range :
diff --git a/src/common/micro/modd_param_ice.F90 b/src/common/micro/modd_param_ice.F90
index f35cb367d..fbf4c51df 100644
--- a/src/common/micro/modd_param_ice.F90
+++ b/src/common/micro/modd_param_ice.F90
@@ -77,6 +77,7 @@ CHARACTER(LEN=1) :: CFRAC_ICE_SHALLOW_MF ! ice fraction for shallow_mf
LOGICAL :: LSEDIM_AFTER ! sedimentation done before (.FALSE.) or after (.TRUE.) microphysics
!
REAL :: XSPLIT_MAXCFL ! Maximum CFL number allowed for SPLIT scheme
+LOGICAL :: LSNOW_T ! Snow parameterization from Wurtz (2021)
END TYPE PARAM_ICE_t
!
TYPE(PARAM_ICE_t), SAVE, TARGET :: PARAM_ICE
@@ -94,7 +95,8 @@ LOGICAL, POINTER :: LWARM => NULL(), &
LCRFLIMIT => NULL(), &
LADJ_BEFORE => NULL(), &
LADJ_AFTER => NULL(), &
- LSEDIM_AFTER => NULL()
+ LSEDIM_AFTER => NULL(),&
+ LSNOW_T => NULL()
REAL, POINTER :: XVDEPOSC => NULL(), &
XFRACM90 => NULL(), &
@@ -132,6 +134,7 @@ SUBROUTINE PARAM_ICE_ASSOCIATE()
LADJ_BEFORE => PARAM_ICE%LADJ_BEFORE
LADJ_AFTER => PARAM_ICE%LADJ_AFTER
LSEDIM_AFTER => PARAM_ICE%LSEDIM_AFTER
+ LSNOW_T => PARAM_ICE%LSNOW_T
!
XVDEPOSC => PARAM_ICE%XVDEPOSC
XFRACM90 => PARAM_ICE%XFRACM90
diff --git a/src/common/micro/modd_rain_ice_descr.F90 b/src/common/micro/modd_rain_ice_descr.F90
index d46a9c1b8..a7b8113de 100644
--- a/src/common/micro/modd_rain_ice_descr.F90
+++ b/src/common/micro/modd_rain_ice_descr.F90
@@ -67,11 +67,13 @@ REAL :: XALPHAC,XNUC,XALPHAC2,XNUC2, XLBEXC ! Cloud droplet distribution p
REAL,DIMENSION(2) :: XLBC ! Cloud droplet distribution parameters
REAL :: XALPHAR,XNUR,XLBEXR,XLBR ! Raindrop distribution parameters
REAL :: XALPHAI,XNUI,XLBEXI,XLBI ! Cloud ice distribution parameters
-REAL :: XALPHAS,XNUS,XLBEXS,XLBS ! Snow/agg. distribution parameters
+REAL :: XALPHAS,XNUS,XLBEXS,XLBS,XNS ! Snow/agg. distribution parameters
REAL :: XALPHAG,XNUG,XLBEXG,XLBG ! Graupel distribution parameters
REAL :: XALPHAH,XNUH,XLBEXH,XLBH ! Hail distribution parameters
!
-REAL :: XLBDAR_MAX,XLBDAS_MAX,XLBDAG_MAX ! Max values allowed for the shape
+REAL :: XFVELOS ! factor for snow fall speed after Thompson (2008)
+REAL :: XTRANS_MP_GAMMAS ! coefficient to convert lambdas for gamma function
+REAL :: XLBDAR_MAX,XLBDAS_MIN,XLBDAS_MAX,XLBDAG_MAX ! Max values allowed for the shape
! parameters (rain,snow,graupeln)
!
REAL,DIMENSION(:),ALLOCATABLE :: XRTMIN ! Min values allowed for the mixing ratios
@@ -145,6 +147,7 @@ REAL, POINTER :: XCEXVT => NULL(), &
XLBI => NULL(), &
XALPHAS => NULL(), &
XNUS => NULL(), &
+ XNS => NULL(), &
XLBEXS => NULL(), &
XLBS => NULL(), &
XALPHAG => NULL(), &
@@ -160,7 +163,10 @@ REAL, POINTER :: XCEXVT => NULL(), &
XLBDAG_MAX => NULL(), &
XCONC_SEA => NULL(), &
XCONC_LAND => NULL(), &
- XCONC_URBAN => NULL()
+ XCONC_URBAN => NULL(), &
+ XFVELOS => NULL(), &
+ XTRANS_MP_GAMMAS => NULL(), &
+ XLBDAS_MIN => NULL()
!
CONTAINS
SUBROUTINE RAIN_ICE_DESCR_ASSOCIATE()
@@ -243,6 +249,10 @@ SUBROUTINE RAIN_ICE_DESCR_ASSOCIATE()
XCONC_SEA => RAIN_ICE_DESCR%XCONC_SEA
XCONC_LAND => RAIN_ICE_DESCR%XCONC_LAND
XCONC_URBAN => RAIN_ICE_DESCR%XCONC_URBAN
+ XNS => RAIN_ICE_DESCR%XNS
+ XFVELOS => RAIN_ICE_DESCR%XFVELOS
+ XTRANS_MP_GAMMAS => RAIN_ICE_DESCR%XTRANS_MP_GAMMAS
+ XLBDAS_MIN => RAIN_ICE_DESCR%XLBDAS_MIN
END SUBROUTINE
!
SUBROUTINE RAIN_ICE_DESCR_ALLOCATE(KRR)
diff --git a/src/common/micro/mode_ice4_fast_rg.F90 b/src/common/micro/mode_ice4_fast_rg.F90
index 0d634dacd..a3457b89e 100644
--- a/src/common/micro/mode_ice4_fast_rg.F90
+++ b/src/common/micro/mode_ice4_fast_rg.F90
@@ -30,6 +30,7 @@ SUBROUTINE ICE4_FAST_RG(CST, PARAMI, ICEP, ICED, KPROMA, KSIZE, LDSOFT, LDCOMPUT
! P. Wautelet 26/04/2019: replace non-standard FLOAT function by REAL function
! P. Wautelet 29/05/2019: remove PACK/UNPACK intrinsics (to get more performance and better OpenACC support)
!! R. El Khatib 24-Aug-2021 Optimizations
+! J. Wurtz 03/2022: New snow characteristics with LSNOW_T
!
!
!* 0. DECLARATIONS
@@ -226,8 +227,13 @@ IF(.NOT. LDSOFT) THEN
WHERE(GDRY(1:KSIZE))
PRG_TEND(1:KSIZE, IRSWETG)=ICEP%XFSDRYG*ZZW(1:KSIZE) & ! RSDRYG
/ ICEP%XCOLSG &
+#if defined(REPRO48) || defined(REPRO55)
*(PLBDAS(1:KSIZE)**(ICED%XCXS-ICED%XBS))*( PLBDAG(1:KSIZE)**ICED%XCXG ) &
*(PRHODREF(1:KSIZE)**(-ICED%XCEXVT-1.)) &
+#else
+ *(PRST(1:KSIZE))*( PLBDAG(1:KSIZE)**ICED%XCXG ) &
+ *(PRHODREF(1:KSIZE)**(-ICED%XCEXVT)) &
+#endif
*( ICEP%XLBSDRYG1/( PLBDAG(1:KSIZE)**2 ) + &
ICEP%XLBSDRYG2/( PLBDAG(1:KSIZE) * PLBDAS(1:KSIZE) ) + &
ICEP%XLBSDRYG3/( PLBDAS(1:KSIZE)**2))
diff --git a/src/common/micro/mode_ice4_fast_rh.F90 b/src/common/micro/mode_ice4_fast_rh.F90
index 19dfe1417..3d13263d7 100644
--- a/src/common/micro/mode_ice4_fast_rh.F90
+++ b/src/common/micro/mode_ice4_fast_rh.F90
@@ -28,6 +28,7 @@ SUBROUTINE ICE4_FAST_RH(CST, PARAMI, ICEP, ICED, KPROMA, KSIZE, LDSOFT, LDCOMPUT
! P. Wautelet 26/04/2019: replace non-standard FLOAT function by REAL function
! P. Wautelet 29/05/2019: remove PACK/UNPACK intrinsics (to get more performance and better OpenACC support)
!! R. El Khatib 24-Aug-2021 Optimizations
+! J. Wurtz 03/2022: New snow characteristics with LSNOW_T
!
!
!* 0. DECLARATIONS
@@ -188,8 +189,13 @@ IF(.NOT. LDSOFT) THEN
!$mnh_expand_where(JL=1:KSIZE)
WHERE(GWET(1:KSIZE))
PRH_TEND(1:KSIZE, IRSWETH)=ICEP%XFSWETH*ZZW(1:KSIZE) & ! RSWETH
+#if defined(REPRO48) || defined(REPRO55)
*( PLBDAS(1:KSIZE)**(ICED%XCXS-ICED%XBS) )*( PLBDAH(1:KSIZE)**ICED%XCXH ) &
*( PRHODREF(1:KSIZE)**(-ICED%XCEXVT-1.) ) &
+#else
+ *( PRST(1:KSIZE))*( PLBDAH(1:KSIZE)**ICED%XCXH ) &
+ *( PRHODREF(1:KSIZE)**(-ICED%XCEXVT) ) &
+#endif
*( ICEP%XLBSWETH1/( PLBDAH(1:KSIZE)**2 ) + &
ICEP%XLBSWETH2/( PLBDAH(1:KSIZE) * PLBDAS(1:KSIZE) ) + &
ICEP%XLBSWETH3/( PLBDAS(1:KSIZE)**2) )
diff --git a/src/common/micro/mode_ice4_fast_rs.F90 b/src/common/micro/mode_ice4_fast_rs.F90
index e04523abd..d1a1afb4b 100644
--- a/src/common/micro/mode_ice4_fast_rs.F90
+++ b/src/common/micro/mode_ice4_fast_rs.F90
@@ -30,6 +30,7 @@ SUBROUTINE ICE4_FAST_RS(CST, PARAMI, ICEP, ICED, KPROMA, KSIZE, LDSOFT, LDCOMPUT
! P. Wautelet 26/04/2019: replace non-standard FLOAT function by REAL function
! P. Wautelet 29/05/2019: remove PACK/UNPACK intrinsics (to get more performance and better OpenACC support)
!! R. El Khatib 24-Aug-2021 Optimizations
+! J. Wurtz 03/2022: New snow characteristics with LSNOW_T
!
!
!* 0. DECLARATIONS
@@ -110,8 +111,14 @@ DO JL=1, KSIZE
PRS_TEND(JL, IFREEZ1)=PKA(JL)*(CST%XTT-PT(JL)) + &
&(PDV(JL)*(CST%XLVTT+(CST%XCPV-CST%XCL)*(PT(JL)-CST%XTT)) &
&*(CST%XESTT-PRS_TEND(JL, IFREEZ1))/(CST%XRV*PT(JL)) )
+#if defined(REPRO48) || defined(REPRO55)
PRS_TEND(JL, IFREEZ1)=PRS_TEND(JL, IFREEZ1)* (ICEP%X0DEPS* PLBDAS(JL)**ICEP%XEX0DEPS + &
& ICEP%X1DEPS*PCJ(JL)*PLBDAS(JL)**ICEP%XEX1DEPS )/ &
+#else
+ PRS_TEND(JL, IFREEZ1)=PRS_TEND(JL, IFREEZ1)* PRST(JL) *(ICEP%X0DEPS* PLBDAS(JL)**ICEP%XEX0DEPS + &
+ & ICEP%X1DEPS*PCJ(JL)*PLBDAS(JL)**(ICEP%XBS+ICEP%XEX1DEPS )* &
+ (1+0.5*(ICEP%XFVELOS/PLBDAS(JL))**ICED%XALPHAS)**(-ICEP%XNUS+ICEP%XEX1DEPS/ICED%XALPHAS))/ &
+#endif
&(PRHODREF(JL)*(CST%XLMTT-CST%XCL*(CST%XTT-PT(JL))))
PRS_TEND(JL, IFREEZ2)=(PRHODREF(JL)*(CST%XLMTT+(CST%XCI-CST%XCL)*(CST%XTT-PT(JL))) ) / &
&(PRHODREF(JL)*(CST%XLMTT-CST%XCL*(CST%XTT-PT(JL))))
@@ -151,7 +158,7 @@ IF(.NOT. LDSOFT) THEN
! 5.1.1 select the PLBDAS
!
DO JJ = 1, IGRIM
- ZVEC1(JJ) = PLBDAS(I1(JJ))
+ ZVEC1(JJ) = (PLBDAS(I1(JJ))**ICED%XALPHAS + ICED%XFVELOS**ICED%XALPHAS)**(1./ICED%XALPHAS)
END DO
!
! 5.1.2 find the next lower indice for the PLBDAS in the geometrical
@@ -180,8 +187,14 @@ IF(.NOT. LDSOFT) THEN
!$mnh_expand_where(JL=1:KSIZE)
WHERE (GRIM(1:KSIZE))
PRS_TEND(1:KSIZE, IRCRIMSS) = ICEP%XCRIMSS * ZZW(1:KSIZE) * PRCT(1:KSIZE) & ! RCRIMSS
+#if defined(REPRO48) || defined(REPRO55)
* PLBDAS(1:KSIZE)**ICEP%XEXCRIMSS &
* PRHODREF(1:KSIZE)**(-ICED%XCEXVT)
+#else
+ * PRST(1:KSIZE)*(1+(ICEP%XFVELOS/PLBDAS(1:KSIZE))**ICED%XALPHAS)**(-ICEP%XNUS+ICEP%XEXCRIMSS/ICED%XALPHAS) &
+ * PRHODREF(1:KSIZE)**(-ICED%XCEXVT+1.) &
+ * (PLBDAS(1:KSIZE)) ** (ICEP%XEXCRIMSS+ICEP%XBS)
+#endif
END WHERE
!$mnh_end_expand_where(JL=1:KSIZE)
!
@@ -213,8 +226,14 @@ IF(.NOT. LDSOFT) THEN
!$mnh_expand_where(JL=1:KSIZE)
WHERE(GRIM(1:KSIZE))
PRS_TEND(1:KSIZE, IRCRIMS)=ICEP%XCRIMSG * PRCT(1:KSIZE) & ! RCRIMS
+#if defined(REPRO48) || defined(REPRO55)
* PLBDAS(1:KSIZE)**ICEP%XEXCRIMSG &
* PRHODREF(1:KSIZE)**(-ICED%XCEXVT)
+#else
+ * PRST(1:KSIZE)*(1+(ICED%XFVELOS/PLBDAS(1:KSIZE))**(ICED%XALPHAS))**(-ICED%XNUS+ICED%XEXCRIMSG/ICED%XALPHAS) &
+ * PRHODREF(1:KSIZE)**(-ICED%XCEXVT+1.) &
+ * PLBDAS(1:KSIZE)**(ICEP%XBS+ICEP%XEXCRIMSG)
+#endif
END WHERE
!$mnh_end_expand_where(JL=1:KSIZE)
@@ -223,10 +242,18 @@ IF(.NOT. LDSOFT) THEN
!$mnh_expand_where(JL=1:KSIZE)
WHERE(GRIM(1:KSIZE))
ZZW6(1:KSIZE) = PRS_TEND(1:KSIZE, IRCRIMS) - PRS_TEND(1:KSIZE, IRCRIMSS) ! RCRIMSG
+#if defined(REPRO48) || defined(REPRO55)
PRS_TEND(1:KSIZE, IRSRIMCG)=ICEP%XSRIMCG * PLBDAS(1:KSIZE)**ICEP%XEXSRIMCG*(1.0-ZZW(1:KSIZE))
+#else
+ PRS_TEND(1:KSIZE, IRSRIMCG)=ICEP%XSRIMCG * PRST(1:KSIZE)*PRHODREF(1:KSIZE)*PLBDAS(1:KSIZE)**(ICEP%XEXSRIMCG+ICEP%XBS)*(1.0-ZZW(1:KSIZE))
+#endif
PRS_TEND(1:KSIZE, IRSRIMCG)=ZZW6(1:KSIZE)*PRS_TEND(1:KSIZE, IRSRIMCG)/ &
MAX(1.E-20, &
+#if defined(REPRO48) || defined(REPRO55)
ICEP%XSRIMCG3*ICEP%XSRIMCG2*PLBDAS(1:KSIZE)**ICEP%XEXSRIMCG2*(1.-ZZW2(1:KSIZE)) - &
+#else
+ ICEP%XSRIMCG3*ICEP%XSRIMCG2*PRST(1:KSIZE)*PRHODREF(1:KSIZE)*PLBDAS(1:KSIZE)***ICEP%XEXSRIMCG2*(1.-ZZW2(1:KSIZE)) - &
+#endif
ICEP%XSRIMCG3*PRS_TEND(1:KSIZE, IRSRIMCG))
END WHERE
!$mnh_end_expand_where(JL=1:KSIZE)
@@ -321,7 +348,11 @@ IF(.NOT. LDSOFT) THEN
!$mnh_expand_where(JL=1:KSIZE)
WHERE(GACC(1:KSIZE))
ZZW6(1:KSIZE) = & !! coef of RRACCS
+#if defined(REPRO48) || defined(REPRO55)
ICEP%XFRACCSS*( PLBDAS(1:KSIZE)**ICED%XCXS )*( PRHODREF(1:KSIZE)**(-ICED%XCEXVT-1.) ) &
+#else
+ ICEP%XFRACCSS*( PRST(1:KSIZE)*PLBDAS(1:KSIZE)**ICED%XBS )*( PRHODREF(1:KSIZE)**(-ICED%XCEXVT) ) &
+#endif
*( ICEP%XLBRACCS1/((PLBDAS(1:KSIZE)**2) ) + &
ICEP%XLBRACCS2/( PLBDAS(1:KSIZE) * PLBDAR(1:KSIZE) ) + &
ICEP%XLBRACCS3/( (PLBDAR(1:KSIZE)**2)) )/PLBDAR(1:KSIZE)**4
@@ -371,7 +402,11 @@ IF(.NOT. LDSOFT) THEN
!$mnh_expand_where(JL=1:KSIZE)
WHERE(GACC(1:KSIZE))
PRS_TEND(1:KSIZE, IRSACCRG) = ICEP%XFSACCRG*ZZW(1:KSIZE)* & ! RSACCRG
+#if defined(REPRO48) || defined(REPRO55)
( PLBDAS(1:KSIZE)**(ICED%XCXS-ICED%XBS) )*( PRHODREF(1:KSIZE)**(-ICED%XCEXVT-1.) ) &
+#else
+ ( PRST(1:KSIZE))*( PRHODREF(1:KSIZE)**(-ICED%XCEXVT) ) &
+#endif
*( ICEP%XLBSACCR1/((PLBDAR(1:KSIZE)**2) ) + &
ICEP%XLBSACCR2/( PLBDAR(1:KSIZE) * PLBDAS(1:KSIZE) ) + &
ICEP%XLBSACCR3/( (PLBDAS(1:KSIZE)**2)) )/PLBDAR(1:KSIZE)
@@ -415,11 +450,22 @@ DO JL=1, KSIZE
!
! compute RSMLT
!
- PRSMLTG(JL) = ICEP%XFSCVMG*MAX(0., (-PRSMLTG(JL) * (ICEP%X0DEPS* PLBDAS(JL)**ICEP%XEX0DEPS + &
- ICEP%X1DEPS*PCJ(JL)*PLBDAS(JL)**ICEP%XEX1DEPS) &
- -(PRS_TEND(JL, IRCRIMS) + PRS_TEND(JL, IRRACCS)) * &
- (PRHODREF(JL)*CST%XCL*(CST%XTT-PT(JL))) &
- ) / (PRHODREF(JL)*CST%XLMTT))
+ PRSMLTG(JL) = ICEP%XFSCVMG*MAX(0., (-PRSMLTG(JL) * &
+#if defined(REPRO48) || defined(REPRO55)
+ (ICEP%X0DEPS* PLBDAS(JL)**ICEP%XEX0DEPS + &
+ ICEP%X1DEPS*PCJ(JL)*PLBDAS(JL)**ICEP%XEX1DEPS) &
+#else
+ PRST(1:KSIZE)*PRHODREF(1:KSIZE) * &
+ (ICEP%X0DEPS* PLBDAS(JL)**(ICEP%XABS+ICEP%XEX0DEPS) + &
+ ICEP%X1DEPS*PCJ(JL)*(1+0.5*(ICED%XFVELOS/PLBDAS(1:KSIZE))**ICED%XALPHAS)**(-ICED%XNUS+ICED%XEX1DEPS/ICED%XALPHAS)*PLBDAS(1:KSIZE)**(ICED%XBS+ICED%XEX1DEPS)) &
+#endif
+ -(PRS_TEND(JL, IRCRIMS) + PRS_TEND(JL, IRRACCS)) * &
+ (PRHODREF(JL)*CST%XCL*(CST%XTT-PT(JL))) &
+ ) / (PRHODREF(JL)*CST%XLMTT))
+ !
+ ! note that RSCVMG = RSMLT*XFSCVMG but no heat is exchanged (at the rate RSMLT)
+ ! because the graupeln produced by this process are still icy!!!
+ !
! When T < XTT, rc is collected by snow (riming) to produce snow and graupel
! When T > XTT, if riming was still enabled, rc would produce snow and graupel with snow becomming graupel (conversion/melting) and graupel becomming rain (melting)
! To insure consistency when crossing T=XTT, rc collected with T>XTT must be transformed in rain.
diff --git a/src/common/micro/mode_ice4_rsrimcg_old.F90 b/src/common/micro/mode_ice4_rsrimcg_old.F90
index 9ee32f9ae..865ed7400 100644
--- a/src/common/micro/mode_ice4_rsrimcg_old.F90
+++ b/src/common/micro/mode_ice4_rsrimcg_old.F90
@@ -119,7 +119,11 @@ IF(.NOT. LDSOFT) THEN
!$mnh_expand_where(JL=1:KSIZE)
WHERE(GRIM(1:KSIZE))
PRSRIMCG_MR(1:KSIZE) = ICEP%XSRIMCG * PLBDAS(1:KSIZE)**ICEP%XEXSRIMCG & ! RSRIMCG
+#if defined(REPRO48) || defined(REPRO55)
* (1.0 - ZZW(1:KSIZE) )/PRHODREF(1:KSIZE)
+#else
+ * (1.0 - ZZW(1:KSIZE) )*PRST(1:KSIZE)
+#endif
PRSRIMCG_MR(1:KSIZE)=MIN(PRST(1:KSIZE), PRSRIMCG_MR(1:KSIZE))
END WHERE
!$mnh_end_expand_where(JL=1:KSIZE)
diff --git a/src/common/micro/mode_ice4_sedimentation_split.F90 b/src/common/micro/mode_ice4_sedimentation_split.F90
index 4864c1472..d2f06a7ba 100644
--- a/src/common/micro/mode_ice4_sedimentation_split.F90
+++ b/src/common/micro/mode_ice4_sedimentation_split.F90
@@ -8,7 +8,7 @@ IMPLICIT NONE
CONTAINS
SUBROUTINE ICE4_SEDIMENTATION_SPLIT(D, CST, ICEP, ICED, PARAMI, &
&PTSTEP, KRR, OSEDIC, PDZZ, &
- &PRHODREF, PPABST, PTHT, PRHODJ, &
+ &PRHODREF, PPABST, PTHT, PT, PRHODJ, &
&PRCS, PRCT, PRRS, PRRT, PRIS, PRIT, PRSS, PRST, PRGS, PRGT,&
&PINPRC, PINPRR, PINPRI, PINPRS, PINPRG, &
&PSEA, PTOWN, &
@@ -27,6 +27,7 @@ SUBROUTINE ICE4_SEDIMENTATION_SPLIT(D, CST, ICEP, ICED, PARAMI, &
!! -------------
!!
! P. Wautelet 10/04/2019: replace ABORT and STOP calls by Print_msg
+! J. Wurtz 03/2022: New snow characteristics with LSNOW_T
!
!
!* 0. DECLARATIONS
@@ -60,6 +61,7 @@ REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PDZZ ! Layer t
REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PRHODREF! Reference density
REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PPABST ! absolute pressure at t
REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PTHT ! Theta at time t
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PT ! Temperature at time t
REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PRHODJ ! Dry density * Jacobian
REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(INOUT) :: PRCS ! Cloud water m.r. source
REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PRCT ! Cloud water m.r. at t
@@ -207,9 +209,8 @@ ENDDO
!* 2.1 for cloud
!
IF (GSEDIC) THEN
- CALL INTERNAL_SEDIM_SPLI(D, CST, ICEP, ICED, KRR, &
- &PARAMI%XSPLIT_MAXCFL, &
- &PRHODREF, ZW, PDZZ, PPABST, PTHT, PTSTEP, &
+ CALL INTERNAL_SEDIM_SPLI(D, CST, ICEP, ICED, PARAMI, KRR, &
+ &PRHODREF, ZW, PDZZ, PPABST, PTHT, PT, PTSTEP, &
&2, &
&ZRCT, PRCS, PINPRC, ZPRCS, &
&ZRAY, ZLBC, ZFSEDC, ZCONC3D, PFPR=PFPR)
@@ -217,36 +218,32 @@ ENDIF
!
!* 2.2 for rain
!
- CALL INTERNAL_SEDIM_SPLI(D, CST, ICEP, ICED, KRR, &
- &PARAMI%XSPLIT_MAXCFL, &
- &PRHODREF, ZW, PDZZ, PPABST, PTHT, PTSTEP, &
+ CALL INTERNAL_SEDIM_SPLI(D, CST, ICEP, ICED, PARAMI, KRR, &
+ &PRHODREF, ZW, PDZZ, PPABST, PTHT, PT, PTSTEP, &
&3, &
&ZRRT, PRRS, PINPRR, ZPRRS, &
&PFPR=PFPR)
!
!* 2.3 for pristine ice
!
- CALL INTERNAL_SEDIM_SPLI(D, CST, ICEP, ICED, KRR, &
- &PARAMI%XSPLIT_MAXCFL, &
- &PRHODREF, ZW, PDZZ, PPABST, PTHT, PTSTEP, &
+ CALL INTERNAL_SEDIM_SPLI(D, CST, ICEP, ICED, PARAMI, KRR, &
+ &PRHODREF, ZW, PDZZ, PPABST, PTHT, PT, PTSTEP, &
&4, &
&ZRIT, PRIS, PINPRI, ZPRIS, &
&PFPR=PFPR)
!
!* 2.4 for aggregates/snow
!
- CALL INTERNAL_SEDIM_SPLI(D, CST, ICEP, ICED, KRR, &
- &PARAMI%XSPLIT_MAXCFL, &
- &PRHODREF, ZW, PDZZ, PPABST, PTHT, PTSTEP, &
+ CALL INTERNAL_SEDIM_SPLI(D, CST, ICEP, ICED, PARAMI, KRR, &
+ &PRHODREF, ZW, PDZZ, PPABST, PTHT, PT, PTSTEP, &
&5, &
&ZRST, PRSS, PINPRS, ZPRSS, &
&PFPR=PFPR)
!
!* 2.5 for graupeln
!
- CALL INTERNAL_SEDIM_SPLI(D, CST, ICEP, ICED, KRR, &
- &PARAMI%XSPLIT_MAXCFL, &
- &PRHODREF, ZW, PDZZ, PPABST, PTHT, PTSTEP, &
+ CALL INTERNAL_SEDIM_SPLI(D, CST, ICEP, ICED, PARAMI, KRR, &
+ &PRHODREF, ZW, PDZZ, PPABST, PTHT, PT, PTSTEP, &
&6, &
&ZRGT, PRGS, PINPRG, ZPRGS, &
&PFPR=PFPR)
@@ -254,9 +251,8 @@ ENDIF
!* 2.6 for hail
!
IF (IRR==7) THEN
- CALL INTERNAL_SEDIM_SPLI(D, CST, ICEP, ICED, KRR, &
- &PARAMI%XSPLIT_MAXCFL, &
- &PRHODREF, ZW, PDZZ, PPABST, PTHT, PTSTEP, &
+ CALL INTERNAL_SEDIM_SPLI(D, CST, ICEP, ICED, PARAMI, KRR, &
+ &PRHODREF, ZW, PDZZ, PPABST, PTHT, PT, PTSTEP, &
&7, &
&ZRHT, PRHS, PINPRH, ZPRHS, &
&PFPR=PFPR)
@@ -270,8 +266,8 @@ CONTAINS
!-------------------------------------------------------------------------------
!
!
-SUBROUTINE INTERNAL_SEDIM_SPLI(D, CST, ICEP, ICED, KRR, &
- &PMAXCFL, PRHODREF, POORHODZ, PDZZ, PPABST, PTHT, PTSTEP, &
+SUBROUTINE INTERNAL_SEDIM_SPLI(D, CST, ICEP, ICED, PARAMI, KRR, &
+ &PRHODREF, POORHODZ, PDZZ, PPABST,PTHT,PT,PTSTEP, &
&KSPE, PRXT, PRXS, PINPRX, PPRXS, &
&PRAY, PLBC, PFSEDC, PCONC3D, PFPR)
!
@@ -281,6 +277,7 @@ SUBROUTINE INTERNAL_SEDIM_SPLI(D, CST, ICEP, ICED, KRR, &
USE MODD_CST, ONLY: CST_t
USE MODD_RAIN_ICE_DESCR, ONLY: RAIN_ICE_DESCR_t
USE MODD_RAIN_ICE_PARAM, ONLY: RAIN_ICE_PARAM_t
+USE MODD_PARAM_ICE, ONLY: PARAM_ICE_t
!
IMPLICIT NONE
!
@@ -290,13 +287,14 @@ TYPE(DIMPHYEX_t), INTENT(IN) :: D
TYPE(CST_t), INTENT(IN) :: CST
TYPE(RAIN_ICE_PARAM_t), INTENT(IN) :: ICEP
TYPE(RAIN_ICE_DESCR_t), INTENT(IN) :: ICED
+TYPE(PARAM_ICE_t), INTENT(IN) :: PARAMI
INTEGER, INTENT(IN) :: KRR
-REAL, INTENT(IN) :: PMAXCFL ! maximum CFL allowed
REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PRHODREF ! Reference density
REAL, DIMENSION(D%NIT,D%NJT,D%NKTB:D%NKTE), INTENT(IN) :: POORHODZ ! One Over (Rhodref times delta Z)
REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PDZZ ! layer thikness (m)
REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PPABST
REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PTHT
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PT
REAL, INTENT(IN) :: PTSTEP ! total timestep
INTEGER, INTENT(IN) :: KSPE ! 1 for rc, 2 for rr...
REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(INOUT) :: PRXT ! mr of specy X
@@ -315,6 +313,7 @@ INTEGER, DIMENSION(D%NIT*D%NJT*D%NKT) :: I1,I2,I3 ! Used to replace the COUNT
LOGICAL :: GPRESENT_PFPR
REAL :: ZINVTSTEP
REAL :: ZZWLBDC, ZRAY, ZZT, ZZWLBDA, ZZCC
+REAL :: ZLBDA
REAL :: ZFSED, ZEXSED
REAL :: ZMRCHANGE
REAL, DIMENSION(D%NIT, D%NJT) :: ZMAX_TSTEP ! Maximum CFL in column
@@ -397,15 +396,44 @@ DO WHILE (ANY(ZREMAINT>0.))
& ALOG(PRHODREF(JI,JJ,JK)*PRXT(JI,JJ,JK)) )**ICEP%XEXCSEDI
ENDIF
ENDDO
+#if defined(REPRO48) || defined(REPRO55)
+#else
+ ELSEIF(KSPE==5) THEN
+ ! ******* for snow
+ ZWSED(:,:,:) = 0.
+ DO JL=1, ISEDIM
+ JI=I1(JL)
+ JJ=I2(JL)
+ JK=I3(JL)
+ IF(PRXT(JI,JJ,JK)> ICED%XRTMIN(KSPE)) THEN
+ IF (PARAMI%LSNOW_T .AND. PT(JI,JJ,JK)>263.15) THEN
+ ZLBDA = MAX(MIN(ICED%XLBDAS_MAX, 10**(14.554-0.0423*PT(JI,JJ,JK))),ICED%XLBDAS_MIN)*ICED%XTRANS_MP_GAMMAS
+ ELSE IF (PARAMI%LSNOW_T) THEN
+ ZLBDA = MAX(MIN(ICED%XLBDAS_MAX, 10**(6.226 -0.0106*PT(JI,JJ,JK))),ICED%XLBDAS_MIN)*ICED%XTRANS_MP_GAMMAS
+ ELSE
+ ZLBDA=MAX(MIN(ICED%XLBDAS_MAX, ICED%XLBS * ( PRHODREF(JI,JJ,JK) * PRXT(JI,JJ,JK) )**ICED%XLBEXS),ICED%XLBDAS_MIN)
+ END IF
+ ZWSED(JI, JJ, JK) = ICEP%XFSEDS * &
+ & PRXT(JI,JJ,JK)* &
+ & PRHODREF(JI,JJ,JK)**(1-ICED%XCEXVT) * &
+ & (1 + (ICED%XFVELOS/ZLBDA)**ICED%XALPHAS)** (-ICED%XNUS+ICEP%XEXSEDS/ICED%XALPHAS) * &
+ & ZLBDA ** (ICED%XBS+ICEP%XEXSEDS)
+
+ ENDIF
+ ENDDO
+#endif
ELSE
! ******* for other species
SELECT CASE(KSPE)
CASE(3)
ZFSED=ICEP%XFSEDR
ZEXSED=ICEP%XEXSEDR
+#if defined(REPRO48) || defined(REPRO55)
CASE(5)
ZFSED=ICEP%XFSEDS
ZEXSED=ICEP%XEXSEDS
+#else
+#endif
CASE(6)
ZFSED=ICEP%XFSEDG
ZEXSED=ICEP%XEXSEDG
@@ -434,7 +462,7 @@ DO WHILE (ANY(ZREMAINT>0.))
JJ=I2(JL)
JK=I3(JL)
IF(PRXT(JI,JJ,JK)>ICED%XRTMIN(KSPE) .AND. ZWSED(JI, JJ, JK)>1.E-20) THEN
- ZMAX_TSTEP(JI, JJ) = MIN(ZMAX_TSTEP(JI, JJ), PMAXCFL * PRHODREF(JI, JJ, JK) * &
+ ZMAX_TSTEP(JI, JJ) = MIN(ZMAX_TSTEP(JI, JJ), PARAMI%XSPLIT_MAXCFL * PRHODREF(JI, JJ, JK) * &
& PRXT(JI, JJ, JK) * PDZZ(JI, JJ, JK) / ZWSED(JI, JJ, JK))
ENDIF
ENDDO
diff --git a/src/common/micro/mode_ice4_sedimentation_stat.F90 b/src/common/micro/mode_ice4_sedimentation_stat.F90
index 609a57e37..d4e52e8fd 100644
--- a/src/common/micro/mode_ice4_sedimentation_stat.F90
+++ b/src/common/micro/mode_ice4_sedimentation_stat.F90
@@ -9,6 +9,7 @@ CONTAINS
SUBROUTINE ICE4_SEDIMENTATION_STAT(D, CST, ICEP, ICED, &
&PTSTEP, KRR, OSEDIC, PDZZ, &
&PRHODREF, PPABST, PTHT, PRHODJ, &
+ &PLBDAS, &
&PRCS, PRCT, PRRS, PRRT, PRIS, PRIT, &
&PRSS, PRST, PRGS, PRGT,&
&PINPRC, PINPRR, PINPRI, PINPRS, PINPRG, &
@@ -32,6 +33,7 @@ SUBROUTINE ICE4_SEDIMENTATION_STAT(D, CST, ICEP, ICED, &
!! Ryad El Khatib 09-Oct-2019 Substantial re-write for optimization
!! (outerunrolling, vectorization, memory cache saving, unrolling)
! P. Wautelet 21/01/2021: initialize untouched part of PFPR
+! J. Wurtz 03/2022: New snow characteristics with LSNOW_T
!
!
!* 0. DECLARATIONS
@@ -61,6 +63,7 @@ REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PRHODREF! Referen
REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PPABST ! absolute pressure at t
REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PTHT ! Theta at time t
REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PRHODJ ! Dry density * Jacobian
+REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PLBDAS ! lambda parameter for snow
REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(INOUT) :: PRCS ! Cloud water m.r. source
REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PRCT ! Cloud water m.r. at t
REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(INOUT) :: PRRS ! Rain water m.r. source
diff --git a/src/common/micro/mode_ice4_slow.F90 b/src/common/micro/mode_ice4_slow.F90
index fc82cebc6..1b0c6525e 100644
--- a/src/common/micro/mode_ice4_slow.F90
+++ b/src/common/micro/mode_ice4_slow.F90
@@ -24,6 +24,7 @@ SUBROUTINE ICE4_SLOW(CST, ICEP, ICED, KPROMA, KSIZE, LDSOFT, LDCOMPUTE, PRHODREF
!! -------------
!!
!! R. El Khatib 24-Aug-2021 Optimizations
+! J. Wurtz 03/2022: New snow characteristics with LSNOW_T
!
!
!* 0. DECLARATIONS
@@ -112,8 +113,15 @@ ENDDO
DO JL=1, KSIZE
IF(PRVT(JL)>ICED%XRTMIN(1) .AND. PRST(JL)>ICED%XRTMIN(5) .AND. LDCOMPUTE(JL)) THEN
IF(.NOT. LDSOFT) THEN
+#if defined(REPRO48) || defined(REPRO55)
PRVDEPS(JL) = ( PSSI(JL)/(PRHODREF(JL)*PAI(JL)) ) * &
( ICEP%X0DEPS*PLBDAS(JL)**ICEP%XEX0DEPS + ICEP%X1DEPS*PCJ(JL)*PLBDAS(JL)**ICEP%XEX1DEPS )
+#else
+ PRVDEPS(JL) = ( PRST(JL)*(PSSI(JL)/PAI(JL)) ) * &
+ ( ICEP%X0DEPS*PLBDAS(JL)**(ICED%XBS+ICEP%XEX0DEPS) + ICEP%X1DEPS*PCJ(JL) * &
+ (1+0.5*(ICED%XFVELOS/PLBDAS(JL))**ICED%XALPHAS)**(-ICED%XNUS+ICED%XEX1DEPS/ICED%XALPHAS) &
+ *(PLBDAS(JL))**(ICED%XBS+ICED%XEX1DEPS) )
+#endif
ENDIF
ELSE
PRVDEPS(JL) = 0.
@@ -127,8 +135,14 @@ DO JL=1, KSIZE
IF(.NOT. LDSOFT) THEN
PRIAGGS(JL) = ICEP%XFIAGGS * EXP( ICEP%XCOLEXIS*(PT(JL)-CST%XTT) ) &
* PRIT(JL) &
+#if defined(REPRO48) || defined(REPRO55)
* PLBDAS(JL)**ICEP%XEXIAGGS &
* PRHODREF(JL)**(-ICED%XCEXVT)
+#else
+ * PRST(JL) * (1+(ICED%XFVELOS/PLBDAS(JL))**ICED%XALPHAS)**(-ICED%XNUS+ICEP%XEXIAGGS/ICED%XALPHAS) &
+ * PRHODREF(JL)**(-ICED%XCEXVT+1.) &
+ * ((PLBDAS(JL))**(ICED%XBS+ICEP%XEXIAGGS))
+#endif
ENDIF
ELSE
PRIAGGS(JL) = 0.
diff --git a/src/common/micro/rain_ice.F90 b/src/common/micro/rain_ice.F90
index 9b8828c60..01ef96e1e 100644
--- a/src/common/micro/rain_ice.F90
+++ b/src/common/micro/rain_ice.F90
@@ -173,6 +173,7 @@
! P. Wautelet 02/2020: use the new data structures and subroutines for budgets
! P. Wautelet 25/02/2020: bugfix: add missing budget: WETH_BU_RRG
!! R. El Khatib 24-Aug-2021 Optimizations
+! J. Wurtz 03/2022: New snow characteristics with LSNOW_T
!-----------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -304,6 +305,7 @@ REAL, DIMENSION(D%NIJT,D%NKT) :: ZRST ! Snow/aggregate m.r. source at t
REAL, DIMENSION(D%NIJT,D%NKT) :: ZRGT ! Graupel m.r. source at t
REAL, DIMENSION(D%NIJT,D%NKT) :: ZRHT ! Hail m.r. source at t
REAL, DIMENSION(D%NIJT,D%NKT) :: ZCITOUT ! Output value for CIT
+REAL, DIMENSION(D%NIJT,D%NKT) :: ZLBDAS ! Modif !lbda parameter snow
!Diagnostics
REAL, DIMENSION(D%NIJT) :: ZINPRI ! Pristine ice instant precip
@@ -468,6 +470,34 @@ DO JK = D%NKTB,D%NKTE
ENDDO
ENDDO
!
+!Compute lambda_snow parameter
+!ZT en KELVIN
+DO JK = D%NKTB,D%NKTE
+ DO JIJ = D%NIJB,D%NIJE
+ ZLBDAS(JIJ,JK)=1000.
+ END DO
+END DO
+DO JK = D%NKTB,D%NKTE
+ DO JIJ = D%NIJB,D%NIJE
+ IF (PARAMI%LSNOW_T) THEN
+ IF (PRST(JIJ,JK)>ICED%XRTMIN(5)) THEN
+ IF(ZT(JIJ,JK)>CST%XTT-10.0) THEN
+ ZLBDAS(JIJ,JK) = MAX(MIN(ICED%XLBDAS_MAX, 10**(14.554-0.0423*ZT(JIJ,JK))),ICED%XLBDAS_MIN)*ICED%XTRANS_MP_GAMMAS
+ ELSE
+ ZLBDAS(JIJ,JK) = MAX(MIN(ICED%XLBDAS_MAX, 10**(6.226-0.0106*ZT(JIJ,JK))),ICED%XLBDAS_MIN)*ICED%XTRANS_MP_GAMMAS
+ END IF
+ END IF
+#if defined(REPRO48) || defined(REPRO55)
+#else
+ ELSE
+ IF (PRST(JIJ,JK).GT.ICED%XRTMIN(5)) THEN
+ ZLBDAS(JIJ,JK) = MAX(MIN(ICED%XLBDAS_MAX,ICED%XLBS*(PRHODREF(JIJ,JK)*PRST(JIJ,JK))**ICED%XLBEXS),ICED%XLBDAS_MIN)
+ END IF
+#endif
+ END IF
+ END DO
+END DO
+!
!-------------------------------------------------------------------------------
!
!* 2. COMPUTE THE SEDIMENTATION (RS) SOURCE
@@ -529,7 +559,7 @@ IF(.NOT. PARAMI%LSEDIM_AFTER) THEN
IF(KRR==7) THEN
CALL ICE4_SEDIMENTATION_SPLIT(D, CST, ICEP, ICED, PARAMI, &
&PTSTEP, KRR, OSEDIC, PDZZ, &
- &PRHODREF, PPABST, PTHT, PRHODJ, &
+ &PRHODREF, PPABST, PTHT, ZT, PRHODJ, &
&PRCS, PRCT, PRRS, PRRT, PRIS, PRIT, PRSS, PRST, PRGS, PRGT,&
&PINPRC, PINPRR, ZINPRI, PINPRS, PINPRG, &
&PSEA=PSEA, PTOWN=PTOWN, &
@@ -537,7 +567,7 @@ IF(.NOT. PARAMI%LSEDIM_AFTER) THEN
ELSE
CALL ICE4_SEDIMENTATION_SPLIT(D, CST, ICEP, ICED, PARAMI, &
&PTSTEP, KRR, OSEDIC, PDZZ, &
- &PRHODREF, PPABST, PTHT, PRHODJ, &
+ &PRHODREF, PPABST, PTHT, ZT, PRHODJ, &
&PRCS, PRCT, PRRS, PRRT, PRIS, PRIT, PRSS, PRST, PRGS, PRGT,&
&PINPRC, PINPRR, ZINPRI, PINPRS, PINPRG, &
&PSEA=PSEA, PTOWN=PTOWN, &
@@ -1612,7 +1642,7 @@ IF(PARAMI%LSEDIM_AFTER) THEN
IF(KRR==7) THEN
CALL ICE4_SEDIMENTATION_SPLIT(D, CST, ICEP, ICED, PARAMI, &
&PTSTEP, KRR, OSEDIC, PDZZ, &
- &PRHODREF, PPABST, PTHT, PRHODJ, &
+ &PRHODREF, PPABST, PTHT, ZT, PRHODJ, &
&PRCS, PRCT, PRRS, PRRT, PRIS, PRIT, PRSS, PRST, PRGS, PRGT,&
&PINPRC, PINPRR, ZINPRI, PINPRS, PINPRG, &
&PSEA=PSEA, PTOWN=PTOWN, &
@@ -1620,7 +1650,7 @@ IF(PARAMI%LSEDIM_AFTER) THEN
ELSE
CALL ICE4_SEDIMENTATION_SPLIT(D, CST, ICEP, ICED, PARAMI, &
&PTSTEP, KRR, OSEDIC, PDZZ, &
- &PRHODREF, PPABST, PTHT, PRHODJ, &
+ &PRHODREF, PPABST, PTHT, ZT, PRHODJ, &
&PRCS, PRCT, PRRS, PRRT, PRIS, PRIT, PRSS, PRST, PRGS, PRGT,&
&PINPRC, PINPRR, ZINPRI, PINPRS, PINPRG, &
&PSEA=PSEA, PTOWN=PTOWN, &
diff --git a/src/common/turb/modi_turb.F90 b/src/common/turb/modi_turb.F90
index 134cc60c3..e2367a4b0 100644
--- a/src/common/turb/modi_turb.F90
+++ b/src/common/turb/modi_turb.F90
@@ -7,6 +7,7 @@ INTERFACE
SUBROUTINE TURB(CST,CSTURB,BUCONF,TURBN,D, &
& KMI,KRR,KRRL,KRRI,HLBCX,HLBCY, &
& KSPLIT,KMODEL_CL,KSV,KSV_LGBEG,KSV_LGEND,HPROGRAM, &
+ & KSV_LIMA_NR, KSV_LIMA_NS, KSV_LIMA_NG, KSV_LIMA_NH, &
& O2D,ONOMIXLG,OFLAT,OLES_CALL,OCOUPLES,OBLOWSNOW, &
& OCOMPUTE_SRC, PRSNOW, &
& OOCEAN,ODEEPOC,ODIAG_IN_RUN, &
@@ -48,6 +49,7 @@ INTEGER, INTENT(IN) :: KRR ! number of moist var.
INTEGER, INTENT(IN) :: KRRL ! number of liquid water var.
INTEGER, INTENT(IN) :: KRRI ! number of ice water var.
INTEGER, INTENT(IN) :: KSV, KSV_LGBEG, KSV_LGEND ! number of scalar variables
+INTEGER, INTENT(IN) :: KSV_LIMA_NR,KSV_LIMA_NS,KSV_LIMA_NG,KSV_LIMA_NH
CHARACTER(LEN=4),DIMENSION(2),INTENT(IN):: HLBCX, HLBCY ! X- and Y-direc LBC
INTEGER, INTENT(IN) :: KSPLIT ! number of time-splitting
INTEGER, INTENT(IN) :: KMODEL_CL ! model number for cloud mixing length
diff --git a/src/common/turb/turb.F90 b/src/common/turb/turb.F90
index 8c445ed9d..6d00c37fb 100644
--- a/src/common/turb/turb.F90
+++ b/src/common/turb/turb.F90
@@ -6,6 +6,7 @@
SUBROUTINE TURB(CST,CSTURB,BUCONF,TURBN,D, &
& KMI,KRR,KRRL,KRRI,HLBCX,HLBCY, &
& KSPLIT,KMODEL_CL,KSV,KSV_LGBEG,KSV_LGEND,HPROGRAM, &
+ & KSV_LIMA_NR, KSV_LIMA_NS, KSV_LIMA_NG, KSV_LIMA_NH, &
& O2D,ONOMIXLG,OFLAT,OLES_CALL,OCOUPLES,OBLOWSNOW, &
& OCOMPUTE_SRC, PRSNOW, &
& OOCEAN,ODEEPOC,ODIAG_IN_RUN, &
@@ -295,6 +296,7 @@ INTEGER, INTENT(IN) :: KRR ! number of moist var.
INTEGER, INTENT(IN) :: KRRL ! number of liquid water var.
INTEGER, INTENT(IN) :: KRRI ! number of ice water var.
INTEGER, INTENT(IN) :: KSV, KSV_LGBEG, KSV_LGEND ! number of scalar variables
+INTEGER, INTENT(IN) :: KSV_LIMA_NR,KSV_LIMA_NS,KSV_LIMA_NG,KSV_LIMA_NH
CHARACTER(LEN=4),DIMENSION(2),INTENT(IN):: HLBCX, HLBCY ! X- and Y-direc LBC
INTEGER, INTENT(IN) :: KSPLIT ! number of time-splitting
INTEGER, INTENT(IN) :: KMODEL_CL ! model number for cloud mixing length
@@ -463,6 +465,9 @@ REAL, DIMENSION(D%NIJT) :: ZTAU11M,ZTAU12M, &
! Virtual Potential Temp. used
! in the Deardorff mixing length computation
!
+!with LIMA, do not change rain, snow, graupel and hail concentrations (mixing ratio is not changed)
+REAL, DIMENSION(D%NIJT,D%NKT,KSV) :: ZRSVS
+!
REAL :: ZEXPL ! 1-TURBN%XIMPL deg of expl.
REAL :: ZRVORD ! RV/RD
REAL :: ZEPS ! XMV / XMD
@@ -520,6 +525,10 @@ IF (TURBN%CTURBLEN=='BL89' .OR. TURBN%CTURBLEN=='RM17' .OR. TURBN%CTURBLEN=='ADA
ZRM(IIJB:IIJE,1:D%NKT,:) = PRT(IIJB:IIJE,1:D%NKT,:)
END IF
!
+!Save LIMA scalar variables sources
+ZRSVS(IIJB:IIJE,1:D%NKT,1:KSV)=PRSVS(IIJB:IIJE,1:D%NKT,1:KSV)
+!
+!
!----------------------------------------------------------------------------
!
!* 2. COMPUTE CONSERVATIVE VARIABLES AND RELATED QUANTITIES
@@ -1000,6 +1009,13 @@ CALL TURB_VER(D, CST,CSTURB,TURBN,KRR, KRRL, KRRI, &
PSSTFL, PSSTFL_C, PSSRFL_C,PSSUFL_C,PSSVFL_C, &
PSSUFL,PSSVFL )
+IF (HCLOUD == 'LIMA') THEN
+ IF (KSV_LIMA_NR.GT.0) PRSVS(:,:,KSV_LIMA_NR) = ZRSVS(:,:,KSV_LIMA_NR)
+ IF (KSV_LIMA_NS.GT.0) PRSVS(:,:,KSV_LIMA_NS) = ZRSVS(:,:,KSV_LIMA_NS)
+ IF (KSV_LIMA_NG.GT.0) PRSVS(:,:,KSV_LIMA_NG) = ZRSVS(:,:,KSV_LIMA_NG)
+ IF (KSV_LIMA_NH.GT.0) PRSVS(:,:,KSV_LIMA_NH) = ZRSVS(:,:,KSV_LIMA_NH)
+END IF
+
IF( BUCONF%LBUDGET_U ) CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_U), 'VTURB', PRUS(:,:) )
IF( BUCONF%LBUDGET_V ) CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_V), 'VTURB', PRVS(:,:) )
IF( BUCONF%LBUDGET_W ) CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_W), 'VTURB', PRWS(:,:) )
@@ -1093,6 +1109,14 @@ IF( TURBN%CTURBDIM == '3DIM' ) THEN
ZTRH, &
PRUS,PRVS,PRWS,PRTHLS,PRRS,PRSVS )
#endif
+ !
+ IF (HCLOUD == 'LIMA') THEN
+ IF (KSV_LIMA_NR.GT.0) PRSVS(:,:,KSV_LIMA_NR) = ZRSVS(:,:,KSV_LIMA_NR)
+ IF (KSV_LIMA_NS.GT.0) PRSVS(:,:,KSV_LIMA_NS) = ZRSVS(:,:,KSV_LIMA_NS)
+ IF (KSV_LIMA_NG.GT.0) PRSVS(:,:,KSV_LIMA_NG) = ZRSVS(:,:,KSV_LIMA_NG)
+ IF (KSV_LIMA_NH.GT.0) PRSVS(:,:,KSV_LIMA_NH) = ZRSVS(:,:,KSV_LIMA_NH)
+ END IF
+ !
IF( BUCONF%LBUDGET_U ) CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_U), 'HTURB', PRUS(:,:) )
IF( BUCONF%LBUDGET_V ) CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_V), 'HTURB', PRVS(:,:) )
IF( BUCONF%LBUDGET_W ) CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_W), 'HTURB', PRWS(:,:) )
diff --git a/src/mesonh/aux/sources_neg_correct.f90 b/src/mesonh/aux/sources_neg_correct.f90
index 9d2cffd55..a1e832734 100644
--- a/src/mesonh/aux/sources_neg_correct.f90
+++ b/src/mesonh/aux/sources_neg_correct.f90
@@ -51,7 +51,8 @@ use modd_budget, only: lbudget_th, lbudget_rv, lbudget_rc, lbudget_rr, lbudg
NBUDGET_RS, NBUDGET_RG, NBUDGET_RH, NBUDGET_SV1, &
tbudgets
use modd_cst, only: xci, xcl, xcpd, xcpv, xlstt, xlvtt, xp00, xrd, xtt
-use modd_nsv, only: nsv_c2r2beg, nsv_c2r2end, nsv_lima_beg, nsv_lima_end, nsv_lima_nc, nsv_lima_nr, nsv_lima_ni
+use modd_nsv, only: nsv_c2r2beg, nsv_c2r2end, nsv_lima_beg, nsv_lima_end, nsv_lima_nc, nsv_lima_nr,&
+ nsv_lima_ni, nsv_lima_ns, nsv_lima_ng, nsv_lima_nh
use modd_param_lima, only: lcold_lima => lcold, lrain_lima => lrain, lspro_lima => lspro, lwarm_lima => lwarm, &
xctmin_lima => xctmin, xrtmin_lima => xrtmin
@@ -78,6 +79,7 @@ integer :: jrmax
integer :: jsv
integer :: isv_lima_end
real, dimension(:, :, :), allocatable :: zt, zexn, zlv, zls, zcph, zcor
+logical, dimension(:, :, :), allocatable :: zmask
if ( krr == 0 ) return
@@ -260,67 +262,109 @@ CLOUD: select case ( hcloud )
!
!
case( 'LIMA' )
+ allocate( zmask ( Size( prths, 1 ), Size( prths, 2 ), Size( prths, 3 ) ) )
! Correction where rc<0 or Nc<0
- if ( lwarm_lima ) then
- where ( prrs(:, :, :, 2) < xrtmin_lima(2) / ptstep .or. prsvs(:, :, :, nsv_lima_nc) < xctmin_lima(2) / ptstep )
- prrs(:, :, :, 1) = prrs(:, :, :, 1) + prrs(:, :, :, 2)
- prths(:, :, :) = prths(:, :, :) - prrs(:, :, :, 2) * zlv(:, :, :) / &
- ( zcph(:, :, :) * zexn(:, :, :) )
- prrs(:, :, :, 2) = 0.
- prsvs(:, :, :, nsv_lima_nc) = 0.
- end where
- where ( prrs(:, :, :, 1) < 0. .and. prrs(:, :, :, 2) > 0. )
- prrs(:, :, :, 1) = prrs(:, :, :, 1) + prrs(:, :, :, 2)
- prths(:, :, :) = prths(:, :, :) - prrs(:, :, :, 2) * zlv(:, :, :) / &
- ( zcph(:, :, :) * zexn(:, :, :) )
- prrs(:, :, :, 2) = 0.
- prsvs(:, :, :, nsv_lima_nc) = 0.
- end where
- end if
+ if ( krr.GE.2 ) then
+ zmask(:,:,:)=(prrs(:, :, :, 2) < xrtmin_lima(2) / ptstep)
+ if (nsv_lima_nc.gt.0) zmask(:,:,:)=(zmask(:,:,:) .or. prsvs(:, :, :, nsv_lima_nc) < xctmin_lima(2) / ptstep )
+ where ( zmask(:,:,:) )
+ prrs(:, :, :, 1) = prrs(:, :, :, 1) + prrs(:, :, :, 2)
+ prths(:, :, :) = prths(:, :, :) - prrs(:, :, :, 2) * zlv(:, :, :) / &
+ ( zcph(:, :, :) * zexn(:, :, :) )
+ prrs(:, :, :, 2) = 0.
+ end where
+ where ( prrs(:, :, :, 1) < 0. .and. prrs(:, :, :, 2) > 0. )
+ prrs(:, :, :, 1) = prrs(:, :, :, 1) + prrs(:, :, :, 2)
+ prths(:, :, :) = prths(:, :, :) - prrs(:, :, :, 2) * zlv(:, :, :) / &
+ ( zcph(:, :, :) * zexn(:, :, :) )
+ prrs(:, :, :, 2) = 0.
+ end where
+ if (nsv_lima_nc.gt.0) then
+ where (prrs(:, :, :, 2) == 0.) prsvs(:, :, :, nsv_lima_nc) = 0.
+ end if
+ end if
! Correction where rr<0 or Nr<0
- if ( lwarm_lima .and. lrain_lima ) then
- where ( prrs(:, :, :, 3) < xrtmin_lima(3) / ptstep .or. prsvs(:, :, :, nsv_lima_nr) < xctmin_lima(3) / ptstep )
- prrs(:, :, :, 1) = prrs(:, :, :, 1) + prrs(:, :, :, 3)
- prths(:, :, :) = prths(:, :, :) - prrs(:, :, :, 3) * zlv(:, :, :) / &
- ( zcph(:, :, :) * zexn(:, :, :) )
- prrs(:, :, :, 3) = 0.
- prsvs(:, :, :, nsv_lima_nr) = 0.
- end where
- end if
+ if ( krr.GE.3 .and. hbudname.ne.'NETUR' ) then
+ zmask(:,:,:)=(prrs(:, :, :, 3) < xrtmin_lima(3) / ptstep)
+ if (nsv_lima_nr.gt.0) zmask(:,:,:)=(zmask(:,:,:) .or. prsvs(:, :, :, nsv_lima_nr) < xctmin_lima(3) / ptstep )
+ where ( zmask(:,:,:) )
+ prrs(:, :, :, 1) = prrs(:, :, :, 1) + prrs(:, :, :, 3)
+ prths(:, :, :) = prths(:, :, :) - prrs(:, :, :, 3) * zlv(:, :, :) / &
+ ( zcph(:, :, :) * zexn(:, :, :) )
+ prrs(:, :, :, 3) = 0.
+ end where
+ if (nsv_lima_nr.gt.0) then
+ where (prrs(:, :, :, 3) == 0.) prsvs(:, :, :, nsv_lima_nr) = 0.
+ end if
+ end if
! Correction where ri<0 or Ni<0
- if ( lcold_lima ) then
- where ( prrs(:, :, :, 4) < xrtmin_lima(4) / ptstep .or. prsvs(:, :, :, nsv_lima_ni) < xctmin_lima(4) / ptstep )
- prrs(:, :, :, 1) = prrs(:, :, :, 1) + prrs(:, :, :, 4)
- prths(:, :, :) = prths(:, :, :) - prrs(:, :, :, 4) * zls(:, :, :) / &
- ( zcph(:, :, :) * zexn(:, :, :) )
- prrs(:, :, :, 4) = 0.
- prsvs(:, :, :, nsv_lima_ni) = 0.
- end where
- if ( hbudname /= 'NETUR' ) then
- do jr = 5, Size( prrs, 4 )
- where ( prrs(:, :, :, jr) < 0. )
- prrs(:, :, :, 1) = prrs(:, :, :, 1) + prrs(:, :, :, jr)
- prths(:, :, :) = prths(:, :, :) - prrs(:, :, :, jr) * zls(:, :, :) / &
- ( zcph(:, :, :) * zexn(:, :, :) )
- prrs(:, :, :, jr) = 0.
- end where
- end do
- end if
- if(krr > 3) then
+ if ( krr.GE.4 ) then
+ zmask(:,:,:)=(prrs(:, :, :, 4) < xrtmin_lima(4) / ptstep)
+ if (nsv_lima_ni.gt.0) zmask(:,:,:)=(zmask(:,:,:) .or. prsvs(:, :, :, nsv_lima_ni) < xctmin_lima(4) / ptstep)
+ where ( zmask(:,:,:) )
+ prrs(:, :, :, 1) = prrs(:, :, :, 1) + prrs(:, :, :, 4)
+ prths(:, :, :) = prths(:, :, :) - prrs(:, :, :, 4) * zls(:, :, :) / &
+ ( zcph(:, :, :) * zexn(:, :, :) )
+ prrs(:, :, :, 4) = 0.
+ end where
allocate( zcor( Size( prths, 1 ), Size( prths, 2 ), Size( prths, 3 ) ) )
where ( prrs(:, :, :, 1) < 0. .and. prrs(:, :, :, 4) > 0. )
- zcor(:, :, :) = Min( -prrs(:, :, :, 1), prrs(:, :, :, 4) )
- prrs(:, :, :, 1) = prrs(:, :, :, 1) + zcor(:, :, :)
- prths(:, :, :) = prths(:, :, :) - zcor(:, :, :) * zls(:, :, :) / &
- ( zcph(:, :, :) * zexn(:, :, :) )
- prrs(:, :, :, 4) = prrs(:, :, :, 4) - zcor(:, :, :)
+ zcor(:, :, :) = Min( -prrs(:, :, :, 1), prrs(:, :, :, 4) )
+ prrs(:, :, :, 1) = prrs(:, :, :, 1) + zcor(:, :, :)
+ prths(:, :, :) = prths(:, :, :) - zcor(:, :, :) * zls(:, :, :) / &
+ ( zcph(:, :, :) * zexn(:, :, :) )
+ prrs(:, :, :, 4) = prrs(:, :, :, 4) - zcor(:, :, :)
end where
deallocate( zcor )
- end if
- end if
-
+ if (nsv_lima_ni.gt.0) then
+ where (prrs(:, :, :, 4) == 0.) prsvs(:, :, :, nsv_lima_ni) = 0.
+ end if
+ end if
+! Snow
+ if ( krr.GE.5 .and. hbudname.ne.'NETUR' ) then
+ zmask(:,:,:)=(prrs(:, :, :, 5) < xrtmin_lima(5) / ptstep)
+ if (nsv_lima_ns.gt.0) zmask(:,:,:)=(zmask(:,:,:) .or. prsvs(:, :, :, nsv_lima_ns) < xctmin_lima(5) / ptstep )
+ where ( zmask(:,:,:) )
+ prrs(:, :, :, 1) = prrs(:, :, :, 1) + prrs(:, :, :, 5)
+ prths(:, :, :) = prths(:, :, :) - prrs(:, :, :, 5) * zls(:, :, :) / &
+ ( zcph(:, :, :) * zexn(:, :, :) )
+ prrs(:, :, :, 5) = 0.
+ end where
+ if (nsv_lima_ns.gt.0) then
+ where (prrs(:, :, :, 5) == 0.) prsvs(:, :, :, nsv_lima_ns) = 0.
+ end if
+ end if
+! Graupel
+ if ( krr.GE.6 .and. hbudname.ne.'NETUR' ) then
+ zmask(:,:,:)=(prrs(:, :, :, 6) < xrtmin_lima(6) / ptstep)
+ if (nsv_lima_ng.gt.0) zmask(:,:,:)=(zmask(:,:,:) .or. prsvs(:, :, :, nsv_lima_ng) < xctmin_lima(6) / ptstep )
+ where ( zmask(:,:,:) )
+ prrs(:, :, :, 1) = prrs(:, :, :, 1) + prrs(:, :, :, 6)
+ prths(:, :, :) = prths(:, :, :) - prrs(:, :, :, 6) * zls(:, :, :) / &
+ ( zcph(:, :, :) * zexn(:, :, :) )
+ prrs(:, :, :, 6) = 0.
+ end where
+ if (nsv_lima_ng.gt.0) then
+ where (prrs(:, :, :, 6) == 0.) prsvs(:, :, :, nsv_lima_ng) = 0.
+ end if
+ end if
+! Hail
+ if ( krr.GE.7 .and. hbudname.ne.'NETUR' ) then
+ zmask(:,:,:)=(prrs(:, :, :, 7) < xrtmin_lima(7) / ptstep)
+ if (nsv_lima_nh.gt.0) zmask(:,:,:)=(zmask(:,:,:) .or. prsvs(:, :, :, nsv_lima_nh) < xctmin_lima(7) / ptstep )
+ where ( zmask(:,:,:) )
+ prrs(:, :, :, 1) = prrs(:, :, :, 1) + prrs(:, :, :, 7)
+ prths(:, :, :) = prths(:, :, :) - prrs(:, :, :, 7) * zls(:, :, :) / &
+ ( zcph(:, :, :) * zexn(:, :, :) )
+ prrs(:, :, :, 7) = 0.
+ end where
+ if (nsv_lima_nh.gt.0) then
+ where (prrs(:, :, :, 7) == 0.) prsvs(:, :, :, nsv_lima_nh) = 0.
+ end if
+ end if
+!
prsvs(:, :, :, nsv_lima_beg : isv_lima_end) = Max( 0.0, prsvs(:, :, :, nsv_lima_beg : isv_lima_end) )
-
+ deallocate(zmask)
end select CLOUD
diff --git a/src/mesonh/ext/deallocate_model1.f90 b/src/mesonh/ext/deallocate_model1.f90
index 0c7bfc0a8..4a940c6d8 100644
--- a/src/mesonh/ext/deallocate_model1.f90
+++ b/src/mesonh/ext/deallocate_model1.f90
@@ -69,6 +69,7 @@ END MODULE MODI_DEALLOCATE_MODEL1
! C. Lac 02/2019: add rain fraction as an output field
! P. Wautelet 07/06/2019: bugfix: deallocate XLSRVM only if allocated
! S. Riette 04/2020: XHL* fields
+! A. Costes 12:2021: Blaze Fire model variables
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -198,6 +199,10 @@ IF ( ASSOCIATED(XCLDFR) .AND. KCALL==2 ) THEN
DEALLOCATE(XCLDFR)
END IF
!
+IF ( ASSOCIATED(XICEFR) .AND. KCALL==2 ) THEN
+ DEALLOCATE(XICEFR)
+END IF
+!
IF ( ASSOCIATED(XRAINFR) .AND. KCALL==2 ) THEN
DEALLOCATE(XRAINFR)
END IF
@@ -220,6 +225,146 @@ END IF
IF (ASSOCIATED(XDUMMY_GR_FIELDS) .AND. KCALL==3 ) THEN
DEALLOCATE(XDUMMY_GR_FIELDS)
END IF
+
+IF (ASSOCIATED(XLSPHI)) THEN
+ DEALLOCATE(XLSPHI)
+END IF
+
+IF (ASSOCIATED(XBMAP)) THEN
+ DEALLOCATE(XBMAP)
+END IF
+
+IF (ASSOCIATED(XFMRFA)) THEN
+ DEALLOCATE(XFMRFA)
+END IF
+
+IF (ASSOCIATED(XFMWF0)) THEN
+ DEALLOCATE(XFMWF0)
+END IF
+
+IF (ASSOCIATED(XFMR0)) THEN
+ DEALLOCATE(XFMR0)
+END IF
+
+IF (ASSOCIATED(XFMR00)) THEN
+ DEALLOCATE(XFMR00)
+END IF
+
+IF (ASSOCIATED(XFMIGNITION)) THEN
+ DEALLOCATE(XFMIGNITION)
+END IF
+
+IF (ASSOCIATED(XFMFUELTYPE)) THEN
+ DEALLOCATE(XFMFUELTYPE)
+END IF
+
+IF (ASSOCIATED(XFIRETAU)) THEN
+ DEALLOCATE(XFIRETAU)
+END IF
+
+IF (ASSOCIATED(XFLUXPARAMH)) THEN
+ DEALLOCATE(XFLUXPARAMH)
+END IF
+
+IF (ASSOCIATED(XFLUXPARAMW)) THEN
+ DEALLOCATE(XFLUXPARAMW)
+END IF
+
+IF (ASSOCIATED(XFIRERW)) THEN
+ DEALLOCATE(XFIRERW)
+END IF
+
+IF (ASSOCIATED(XFMASE)) THEN
+ DEALLOCATE(XFMASE)
+END IF
+
+IF (ASSOCIATED(XFMAWC)) THEN
+ DEALLOCATE(XFMAWC)
+END IF
+
+IF (ASSOCIATED(XFMWALKIG)) THEN
+ DEALLOCATE(XFMWALKIG)
+END IF
+
+IF (ASSOCIATED(XFMFLUXHDH)) THEN
+ DEALLOCATE(XFMFLUXHDH)
+END IF
+
+IF (ASSOCIATED(XFMFLUXHDW)) THEN
+ DEALLOCATE(XFMFLUXHDW)
+END IF
+
+IF (ASSOCIATED(XFMHWS)) THEN
+ DEALLOCATE(XFMHWS)
+END IF
+
+IF (ASSOCIATED(XFMWINDU)) THEN
+ DEALLOCATE(XFMWINDU)
+END IF
+
+IF (ASSOCIATED(XFMWINDV)) THEN
+ DEALLOCATE(XFMWINDV)
+END IF
+
+IF (ASSOCIATED(XFMWINDW)) THEN
+ DEALLOCATE(XFMWINDW)
+END IF
+
+IF (ASSOCIATED(XFMGRADOROX)) THEN
+ DEALLOCATE(XFMGRADOROX)
+END IF
+
+IF (ASSOCIATED(XFMGRADOROY)) THEN
+ DEALLOCATE(XFMGRADOROY)
+END IF
+
+IF (ASSOCIATED(XGRADLSPHIX)) THEN
+ DEALLOCATE(XGRADLSPHIX)
+END IF
+
+IF (ASSOCIATED(XGRADLSPHIY)) THEN
+ DEALLOCATE(XGRADLSPHIY)
+END IF
+
+IF (ASSOCIATED(XFIREWIND)) THEN
+ DEALLOCATE(XFIREWIND)
+END IF
+
+IF (ASSOCIATED(XLSPHI2D)) THEN
+ DEALLOCATE(XLSPHI2D)
+END IF
+
+IF (ASSOCIATED(XGRADLSPHIX2D)) THEN
+ DEALLOCATE(XGRADLSPHIX2D)
+END IF
+
+IF (ASSOCIATED(XGRADLSPHIY2D)) THEN
+ DEALLOCATE(XGRADLSPHIY2D)
+END IF
+
+IF (ASSOCIATED(XGRADMASKX)) THEN
+ DEALLOCATE(XGRADMASKX)
+END IF
+
+IF (ASSOCIATED(XGRADMASKY)) THEN
+ DEALLOCATE(XGRADMASKY)
+END IF
+
+IF (ASSOCIATED(XSURFRATIO2D)) THEN
+ DEALLOCATE(XSURFRATIO2D)
+END IF
+
+IF (ASSOCIATED(XLSDIFFUX2D)) THEN
+ DEALLOCATE(XLSDIFFUX2D)
+END IF
+
+IF (ASSOCIATED(XLSDIFFUY2D)) THEN
+ DEALLOCATE(XLSDIFFUY2D)
+END IF
+
+IF (ASSOCIATED(XFIRERW2D)) THEN
+ DEALLOCATE(XFIRERW2D)
+END IF
!
!* 3. Module MODD_GRID$n
!
diff --git a/src/mesonh/ext/ground_paramn.f90 b/src/mesonh/ext/ground_paramn.f90
index e27d56f95..84022d589 100644
--- a/src/mesonh/ext/ground_paramn.f90
+++ b/src/mesonh/ext/ground_paramn.f90
@@ -10,12 +10,14 @@ MODULE MODI_GROUND_PARAM_n
INTERFACE
!
SUBROUTINE GROUND_PARAM_n(D, PSFTH, PSFRV, PSFSV, PSFCO2, PSFU, PSFV, &
- PDIR_ALB, PSCA_ALB, PEMIS, PTSRAD )
+ PDIR_ALB, PSCA_ALB, PEMIS, PTSRAD, KTCOUNT, TPFILE )
!
USE MODD_DIMPHYEX, ONLY: DIMPHYEX_t
!* surface fluxes
! --------------
!
+USE MODD_IO, ONLY: TFILEDATA
+!
TYPE(DIMPHYEX_t), INTENT(IN) :: D
REAL, DIMENSION(:,:), INTENT(OUT) :: PSFTH ! surface flux of potential temperature (Km/s)
REAL, DIMENSION(:,:), INTENT(OUT) :: PSFRV ! surface flux of water vapor (m/s*kg/kg)
@@ -33,6 +35,8 @@ REAL, DIMENSION(:,:,:), INTENT(OUT) :: PSCA_ALB ! diffuse albedo for each spect
REAL, DIMENSION(:,:,:), INTENT(OUT) :: PEMIS ! surface emissivity (-)
REAL, DIMENSION(:,:), INTENT(OUT) :: PTSRAD ! surface radiative temperature (K)
!
+INTEGER, INTENT(IN) :: KTCOUNT ! temporal iteration count
+TYPE(TFILEDATA), INTENT(IN) :: TPFILE ! Synchronous output file
END SUBROUTINE GROUND_PARAM_n
!
END INTERFACE
@@ -41,7 +45,7 @@ END MODULE MODI_GROUND_PARAM_n
!
! ######################################################################
SUBROUTINE GROUND_PARAM_n(D, PSFTH, PSFRV, PSFSV, PSFCO2, PSFU, PSFV, &
- PDIR_ALB, PSCA_ALB, PEMIS, PTSRAD )
+ PDIR_ALB, PSCA_ALB, PEMIS, PTSRAD, KTCOUNT, TPFILE )
! #######################################################################
!
!
@@ -114,6 +118,7 @@ END MODULE MODI_GROUND_PARAM_n
!! (Bielli S.) 02/2019 Sea salt : significant sea wave height influences salt emission; 5 salt modes
! P. Wautelet 20/05/2019: add name argument to ADDnFIELD_ll + new ADD4DFIELD_ll subroutine
! P. Wautelet 09/02/2022: bugfix: add missing XCURRENT_LEI computation
+! A. Costes 12/2021: Blaze Fire model
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -135,7 +140,12 @@ USE MODD_DIMPHYEX, ONLY : DIMPHYEX_t
USE MODD_PARAMETERS, ONLY : JPVEXT, XUNDEF
USE MODD_DYN_n, ONLY : XTSTEP
USE MODD_CH_MNHC_n, ONLY : LUSECHEM
-USE MODD_FIELD_n, ONLY : XUT, XVT, XWT, XTHT, XRT, XPABST, XSVT, XTKET, XZWS
+USE MODD_CH_M9_n, ONLY : CNAMES
+USE MODD_FIELD_n, ONLY : XUT, XVT, XWT, XTHT, XRT, XPABST, XSVT, XTKET, XZWS,&
+XLSPHI, XBMAP, XFMR0, XFMRFA, XFMWF0, XFMR00, XFMIGNITION, XFMFUELTYPE,&
+XFIRETAU, XFLUXPARAMH, XFLUXPARAMW, XFIRERW, XFMASE, XFMAWC, XFMWALKIG,&
+XFMFLUXHDH, XFMFLUXHDW, XRTHS, XRRS, XFMHWS, XFMWINDU, XFMWINDV, XFMWINDW, XGRADLSPHIX, &
+XGRADLSPHIY, XFIREWIND, XFMGRADOROX, XFMGRADOROY
USE MODD_METRICS_n, ONLY : XDXX, XDYY, XDZZ
USE MODD_DIM_n, ONLY : NKMAX
USE MODD_GRID_n, ONLY : XLON, XZZ, XDIRCOSXW, XDIRCOSYW, XDIRCOSZW, &
@@ -188,6 +198,14 @@ USE MODD_TIME
!
USE MODD_PARAM_LIMA, ONLY : MSEDC=>LSEDC
!
+USE MODD_FIRE
+USE MODD_FIELD
+USE MODI_FIRE_MODEL
+USE MODD_CONF, ONLY : NVERB, NHALO
+USE MODE_MNH_TIMING, ONLY : SECOND_MNH2
+USE MODE_MSG
+USE MODD_IO, ONLY: TFILEDATA
+!
IMPLICIT NONE
!
!
@@ -214,6 +232,8 @@ REAL, DIMENSION(:,:,:), INTENT(OUT) :: PSCA_ALB ! diffuse albedo for each spect
REAL, DIMENSION(:,:,:), INTENT(OUT) :: PEMIS ! surface emissivity (-)
REAL, DIMENSION(:,:), INTENT(OUT) :: PTSRAD ! surface radiative temperature (K)
!
+INTEGER, INTENT(IN) :: KTCOUNT ! temporal iteration count
+TYPE(TFILEDATA), INTENT(IN) :: TPFILE ! Synchronous output file
!
!-------------------------------------------------------------------------------
!
@@ -359,6 +379,16 @@ CHARACTER(LEN=6), DIMENSION(:), ALLOCATABLE :: YSV_SURF ! name of the scalar var
REAL :: ZTIMEC
INTEGER :: ILUOUT ! logical unit
!
+! Fire model
+REAL, DIMENSION(2) :: ZFIRETIME1, ZFIRETIME2 ! CPU time for Blaze perf profiling
+REAL, DIMENSION(2) :: ZGRADTIME1, ZGRADTIME2 ! CPU time for Blaze perf profiling
+REAL, DIMENSION(2) :: ZPROPAGTIME1, ZPROPAGTIME2 ! CPU time for Blaze perf profiling
+REAL, DIMENSION(2) :: ZFLUXTIME1, ZFLUXTIME2 ! CPU time for Blaze perf profiling
+REAL, DIMENSION(2) :: ZROSWINDTIME1, ZROSWINDTIME2 ! CPU time for Blaze perf profiling
+REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: ZFIREFUELMAP ! Fuel map
+CHARACTER(LEN=7) :: YFUELMAPFILE ! Fuel Map file name
+TYPE(LIST_ll), POINTER :: TZFIELDFIRE_ll ! list of fields to exchange
+
!-------------------------------------------------------------------------------
!
!
@@ -668,6 +698,158 @@ WHERE (ZSFU(:,:)/=XUNDEF .AND. ZWIND(:,:)>0.)
PSFV(:,:) = - SQRT(ZSFU**2+ZSFV**2) * ZVA(:,:) / ZWIND(:,:) / XRHODREF(:,:,IKB)
END WHERE
!
+
+!* 2.1 Blaze Fire Model
+! ----------------
+!
+IF (LBLAZE) THEN
+ ! get start time
+ CALL SECOND_MNH2( ZFIRETIME1 )
+
+ !* 2.1.1 Local variables allocation
+ ! --------------------------
+ !
+
+ ! Parallel fuel
+ NULLIFY(TZFIELDFIRE_ll)
+ IF (KTCOUNT <= 1) THEN
+ ! fuelmap
+ SELECT CASE (CPROPAG_MODEL)
+ CASE('SANTONI2011')
+ !
+ ALLOCATE( ZFIREFUELMAP(SIZE(XLSPHI,1), SIZE(XLSPHI,2), SIZE(XLSPHI,3), 22) );
+ ! Parallel fuel
+ CALL ADD4DFIELD_ll( TZFIELDFIRE_ll, ZFIREFUELMAP(:,:,:,1::22), 'MODEL_n::ZFIREFUELMAP' )
+ ! Default value
+ ZFIREFUELMAP(:,:,:,:) = 0.
+ END SELECT
+
+ !* 2.1.2 Read fuel map file
+ ! ------------------
+ !
+ ! Fuel map file name
+ YFUELMAPFILE = 'FuelMap'
+ !
+ CALL FIRE_READFUEL( TPFILE, ZFIREFUELMAP, XFMIGNITION, XFMWALKIG )
+
+ !* 2.1.3 Ignition LS function with ignition map
+ ! --------------------------------------
+ !
+ SELECT CASE (CFIRE_CPL_MODE)
+ CASE('2WAYCPL', 'ATM2FIR')
+ ! force ignition
+ WHERE (XFMIGNITION <= TDTCUR%XTIME ) XLSPHI = 1.
+ ! walking ignition
+ CALL FIRE_LS_RECONSTRUCTION_FROM_BMAP( XLSPHI, XFMWALKIG, 0.)
+ !
+ !* 2.1.4 Update BMAP
+ ! -----------
+ !
+ WHERE (XLSPHI >= .5 .AND. XBMAP < 0) XBMAP = TDTCUR%XTIME
+ !
+ CASE('FIR2ATM')
+ CALL FIRE_READBMAP(TPFILE,XBMAP)
+
+ END SELECT
+ !
+ !* 2.1.5 Compute R0, A, Wf0, R00
+ ! -----------------------
+ !
+ SELECT CASE (CPROPAG_MODEL)
+ CASE('SANTONI2011')
+ CALL FIRE_NOWINDROS( ZFIREFUELMAP, XFMR0, XFMRFA, XFMWF0, XFMR00, XFMFUELTYPE, XFIRETAU, XFLUXPARAMH, &
+ XFLUXPARAMW, XFMASE, XFMAWC )
+ END SELECT
+ !
+ !* 2.1.6 Compute orographic gradient
+ ! ---------------------------
+ CALL FIRE_GRAD_OROGRAPHY( XZS, XFMGRADOROX, XFMGRADOROY )
+ !
+ !* 2.1.7 Test halo size
+ ! --------------
+ IF (NHALO < 2 .AND. NFIRE_WENO_ORDER == 3) THEN
+ WRITE(ILUOUT,'(A/A)') 'ERROR BLAZE-FIRE : WENO3 fire gradient calculation needs NHALO >= 2'
+ !callabortstop
+ CALL PRINT_MSG(NVERB_FATAL,'GEN','GROUND_PARAM_n','')
+ ELSEIF (NHALO < 3 .AND. NFIRE_WENO_ORDER == 5) THEN
+ WRITE(ILUOUT,'(A/A)') 'ERROR : WENO5 fire gradient calculation needs NHALO >= 3'
+ !callabortstop
+ CALL PRINT_MSG(NVERB_FATAL,'GEN','GROUND_PARAM_n','')
+ END IF
+ !
+ END IF
+ !
+ !* 2.1.6 Compute grad of level set function phi
+ ! --------------------------------------
+ !
+ SELECT CASE (CFIRE_CPL_MODE)
+ CASE('2WAYCPL', 'ATM2FIR')
+ ! get time 1
+ CALL SECOND_MNH2( ZGRADTIME1 )
+ CALL FIRE_GRADPHI( XLSPHI, XGRADLSPHIX, XGRADLSPHIY )
+
+ ! get time 2
+ CALL SECOND_MNH2( ZGRADTIME2 )
+ XGRADPERF = XGRADPERF + ZGRADTIME2 - ZGRADTIME1
+ !
+ !* 2.1.7 Get horizontal wind speed projected on LS gradient direction
+ ! ------------------------------------------------------------
+ !
+ CALL FIRE_GETWIND( XUT, XVT, XWT, XGRADLSPHIX, XGRADLSPHIY, XFIREWIND, KTCOUNT, XTSTEP, XFMGRADOROX, XFMGRADOROY )
+ !
+ !* 2.1.8 Compute ROS XFIRERW with wind
+ ! -----------------------------
+ !
+ !
+ SELECT CASE (CPROPAG_MODEL)
+ CASE('SANTONI2011')
+ CALL FIRE_RATEOFSPREAD( XFMFUELTYPE, XFMR0, XFMRFA, XFMWF0, XFMR00, XFIREWIND, XGRADLSPHIX, XGRADLSPHIY, &
+ XFMGRADOROX, XFMGRADOROY, XFIRERW )
+ END SELECT
+ CALL SECOND_MNH2( ZROSWINDTIME2 )
+ XROSWINDPERF = XROSWINDPERF + ZROSWINDTIME2 - ZGRADTIME2
+ !
+ !* 2.1.8 Integrate model on atm time step to propagate
+ ! ---------------------------------------------
+ !
+ SELECT CASE (CPROPAG_MODEL)
+ CASE('SANTONI2011')
+ CALL FIRE_PROPAGATE( XLSPHI, XBMAP, XFMIGNITION, XFMWALKIG, XGRADLSPHIX, XGRADLSPHIY, XTSTEP, XFIRERW )
+ END SELECT
+ CALL SECOND_MNH2( ZPROPAGTIME2 )
+ XPROPAGPERF = XPROPAGPERF + ZPROPAGTIME2 - ZROSWINDTIME2
+ !
+ CASE('FIR2ATM')
+ !
+ CALL SECOND_MNH2( ZPROPAGTIME1 )
+ CALL FIRE_LS_RECONSTRUCTION_FROM_BMAP( XLSPHI, XBMAP, XTSTEP )
+ CALL SECOND_MNH2( ZPROPAGTIME2 )
+ XPROPAGPERF = XPROPAGPERF + ZPROPAGTIME2 - ZPROPAGTIME1
+ XGRADPERF(:) = 0.
+ !
+ END SELECT
+ !
+ !* 2.1.8 Compute fluxes
+ ! --------------
+ !
+ SELECT CASE (CFIRE_CPL_MODE)
+ CASE('2WAYCPL','FIR2ATM')
+ CALL SECOND_MNH2( ZFLUXTIME1 )
+ ! 2 way coupling
+ CALL FIRE_HEATFLUXES( XLSPHI, XBMAP, XFIRETAU, XTSTEP, XFLUXPARAMH, XFLUXPARAMW, XFMFLUXHDH, XFMFLUXHDW, XFMASE, XFMAWC )
+ ! vertical distribution of fire heat fluxes
+ CALL FIRE_VERTICALFLUXDISTRIB( XFMFLUXHDH, XFMFLUXHDW, XRTHS, XRRS, ZSFTS, XEXNREF, XRHODJ, XRT, XRHODREF )
+ !
+ CALL SECOND_MNH2( ZFLUXTIME2 )
+ XFLUXPERF = XFLUXPERF + ZFLUXTIME2 - ZFLUXTIME1
+ CASE DEFAULT
+ XFLUXPERF(:) = 0.
+ END SELECT
+ ! get end time
+ CALL SECOND_MNH2( ZFIRETIME2 )
+ ! add to Blaze time
+ XFIREPERF = XFIREPERF + ZFIRETIME2 - ZFIRETIME1
+END IF
!* conversion from H (W/m2) to w'Theta'
!
PSFTH(:,:) = ZSFTH(:,:) / XCPD / XRHODREF(:,:,IKB)
@@ -691,7 +873,7 @@ END IF
IF (LUSECHEM) THEN
DO JSV=NSV_CHEMBEG,NSV_CHEMEND
PSFSV(:,:,JSV) = ZSFTS(:,:,JSV) * XMD / ( XAVOGADRO * XRHODREF(:,:,IKB))
- IF ((LCHEMDIAG).AND.(CPROGRAM == 'DIAG ')) XCHFLX(:,:,JSV) = PSFSV(:,:,JSV)
+ IF ((LCHEMDIAG).AND.(CPROGRAM == 'DIAG ')) XCHFLX(:,:,JSV-NSV_CHEMBEG+1) = PSFSV(:,:,JSV)
END DO
ELSE
PSFSV(:,:,NSV_CHEMBEG:NSV_CHEMEND) = 0.
@@ -756,7 +938,7 @@ IF (LDIAG_IN_RUN) THEN
XCURRENT_SFCO2(:,:) = ZSFCO2(:,:)
XCURRENT_DSTAOD(:,:)=0.0
XCURRENT_SLTAOD(:,:)=0.0
- IF (CRAD=='ECMW') THEN
+ IF (CRAD/='NONE') THEN
XCURRENT_LWD (:,:) = XFLALWD(:,:)
XCURRENT_SWD (:,:) = SUM(XDIRSRFSWD(:,:,:)+XSCAFLASWD(:,:,:),DIM=3)
XCURRENT_LWU (:,:) = XLWU(:,:,IKB)
@@ -797,6 +979,12 @@ IF (LDIAG_IN_RUN) THEN
CALL CLEANLIST_ll(TZFIELDSURF_ll)
END IF
!
+IF (LBLAZE) THEN
+ IF (KTCOUNT <= 1) THEN
+ DEALLOCATE(ZFIREFUELMAP)
+ END IF
+ CALL CLEANLIST_ll(TZFIELDFIRE_ll)
+END IF
!==================================================================================
!
CONTAINS
diff --git a/src/mesonh/ext/ini_radar.f90 b/src/mesonh/ext/ini_radar.f90
index f0f5e0307..dbc94a726 100644
--- a/src/mesonh/ext/ini_radar.f90
+++ b/src/mesonh/ext/ini_radar.f90
@@ -177,6 +177,7 @@ XLBR = ( XAR*XCCR*MOMG(XALPHAR,XNUR,XBR) )**(-XLBEXR)
XLBEXI = 1.0/(-XBI)
XLBI = ( XAI*MOMG(XALPHAI,XNUI,XBI) )**(-XLBEXI)
!
+XNS = 1.0/(XAS*MOMG(XALPHAS,XNUS,XBS))
XLBEXS = 1.0/(XCXS-XBS)
XLBS = ( XAS*XCCS*MOMG(XALPHAS,XNUS,XBS) )**(-XLBEXS)
!
diff --git a/src/mesonh/ext/lesn.f90 b/src/mesonh/ext/lesn.f90
index 95392487a..7983b1217 100644
--- a/src/mesonh/ext/lesn.f90
+++ b/src/mesonh/ext/lesn.f90
@@ -109,6 +109,7 @@ REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZEW
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZINDCLD !indice cloud si rc>0
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZINDCLD2 !indice cloud rc>1E-5
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZCLDFR_LES! CLDFR on LES vertical grid
+REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZICEFR_LES! ICEFR on LES vertical grid
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZRAINFR_LES! RAINFR on LES vertical grid
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZMASSF ! massflux=rho*w
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZREHU ! relative humidity
@@ -349,9 +350,11 @@ END IF
IF (LUSERI) THEN
ALLOCATE(ZRI_LES (IIU,IJU,NLES_K))
ALLOCATE(ZIWP_LES(IIU,IJU))
+ ALLOCATE(ZICEFR_LES(IIU,IJU,NLES_K))
ELSE
ALLOCATE(ZRI_LES (0,0,0))
ALLOCATE(ZIWP_LES(0,0))
+ ALLOCATE(ZICEFR_LES(0,0,0))
END IF
IF (LUSERS) THEN
ALLOCATE(ZRS_LES (IIU,IJU,NLES_K))
@@ -602,6 +605,7 @@ IF (LUSERI) THEN
END DO
CALL LES_MEAN_ll ( ZIWP_LES, LLES_CURRENT_CART_MASK(:,:,1), &
XLES_IWP(NLES_CURRENT_TCOUNT) )
+ CALL LES_VER_INT( XICEFR(:,:,:) ,ZICEFR_LES )
END IF
IF (LUSERS) THEN
IRR = IRR + 1
@@ -816,6 +820,8 @@ END IF
IF (LUSERI) &
CALL LES_MEAN_ll ( ZRI_LES, LLES_CURRENT_CART_MASK, &
XLES_MEAN_Ri(:,NLES_CURRENT_TCOUNT,1) )
+ CALL LES_MEAN_ll ( ZICEFR_LES, LLES_CURRENT_CART_MASK, &
+ XLES_MEAN_If(:,NLES_CURRENT_TCOUNT,1) )
!
IF (LUSERS) &
CALL LES_MEAN_ll ( ZRS_LES, LLES_CURRENT_CART_MASK, &
@@ -1050,6 +1056,7 @@ DEALLOCATE(ZINDCLD2 )
DEALLOCATE(ZINDCLD2D )
DEALLOCATE(ZINDCLD2D2)
DEALLOCATE(ZCLDFR_LES)
+DEALLOCATE(ZICEFR_LES)
DEALLOCATE(ZRAINFR_LES)
DEALLOCATE(ZMASSF )
DEALLOCATE(ZTEMP )
diff --git a/src/mesonh/ext/phys_paramn.f90 b/src/mesonh/ext/phys_paramn.f90
index 1870c55dc..f84db2c54 100644
--- a/src/mesonh/ext/phys_paramn.f90
+++ b/src/mesonh/ext/phys_paramn.f90
@@ -237,6 +237,7 @@ END MODULE MODI_PHYS_PARAM_n
! C. Lac 11/2019: correction in the drag formula and application to building in addition to tree
! F. Auguste 02/2021: add IBM
! JL Redelsperger 03/2021: add the SW flux penetration for Ocean model case
+! A. Costes 12/2021: add Blaze fire model
!!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -289,7 +290,8 @@ USE MODD_NESTING, ONLY : XWAY,NDAD, NDXRATIO_ALL, NDYRATIO_ALL
USE MODD_NSV, ONLY : NSV, NSV_LGBEG, NSV_LGEND, &
NSV_SLTBEG,NSV_SLTEND,NSV_SLT,&
NSV_AERBEG,NSV_AEREND, &
- NSV_DSTBEG,NSV_DSTEND, NSV_DST
+ NSV_DSTBEG,NSV_DSTEND, NSV_DST,&
+ NSV_LIMA_NR,NSV_LIMA_NS,NSV_LIMA_NG,NSV_LIMA_NH
USE MODD_OCEANH
USE MODD_OUT_n
USE MODD_PARAM_C2R2, ONLY : LSEDC
@@ -644,7 +646,7 @@ IF (CRAD /='NONE') THEN
!
IF (CRAD =='ECMW' .OR. CRAD =='ECRA') THEN
IF (GRAD .AND. NRR.LE.3 ) THEN
- IF( MAXVAL(XCLDFR(:,:,:)).LE. 1.E-10 .AND. OCLOUD_ONLY ) THEN
+ IF( MAX(MAXVAL(XCLDFR(:,:,:)),MAXVAL(XICEFR(:,:,:))).LE. 1.E-10 .AND. OCLOUD_ONLY ) THEN
GRAD = .FALSE. ! only the cloudy verticals would be
! refreshed but there is no clouds
END IF
@@ -759,7 +761,7 @@ CALL SUNPOS_n ( XZENITH, ZCOSZEN, ZSINZEN, ZAZIMSOL )
COPWLW, COPILW, XFUDG, &
NDLON, NFLEV, NRAD_DIAG, NFLUX, NRAD, NAER, NSWB_OLD, NSWB_MNH, NLWB_MNH, &
NSTATM, NRAD_COLNBR, ZCOSZEN, XSEA, XCORSOL, &
- XDIR_ALB, XSCA_ALB, XEMIS, XCLDFR, XCCO2, XTSRAD, XSTATM, XTHT, XRT, &
+ XDIR_ALB, XSCA_ALB, XEMIS, MAX(XCLDFR,XICEFR), XCCO2, XTSRAD, XSTATM, XTHT, XRT, &
XPABST, XOZON, XAER,XDST_WL, XAER_CLIM, XSVT, &
XDTHRAD, XFLALWD, XDIRFLASWD, XSCAFLASWD, XRHODREF, XZZ , &
XRADEFF, XSWU, XSWD, XLWU, XLWD, XDTHRADSW, XDTHRADLW )
@@ -1251,7 +1253,7 @@ IF (CSURF=='EXTE') THEN
DEALLOCATE( ZSAVE_DIRFLASWD,ZSAVE_SCAFLASWD,ZSAVE_DIRSRFSWD)
END IF
CALL GROUND_PARAM_n(YLDIMPHYEX,ZSFTH, ZSFRV, ZSFSV, ZSFCO2, ZSFU, ZSFV, &
- ZDIR_ALB, ZSCA_ALB, ZEMIS, ZTSRAD )
+ ZDIR_ALB, ZSCA_ALB, ZEMIS, ZTSRAD, KTCOUNT, TPFILE )
!
IF (LIBM) THEN
WHERE(XIBM_LS(:,:,IKB,1).GT.-XIBM_EPSI)
@@ -1550,7 +1552,9 @@ LHARAT = .FALSE.
!
CALL TURB( CST,CSTURB, TBUCONF, TURBN,YLDIMPHYEX,&
IMI, NRR, NRRL, NRRI, CLBCX, CLBCY, 1, NMODEL_CLOUD, &
- NSV, NSV_LGBEG, NSV_LGEND,CPROGRAM, L2D, LNOMIXLG,LFLAT, &
+ NSV, NSV_LGBEG, NSV_LGEND,CPROGRAM, &
+ NSV_LIMA_NR, NSV_LIMA_NS, NSV_LIMA_NG, NSV_LIMA_NH, &
+ L2D, LNOMIXLG,LFLAT, &
LLES_CALL, LCOUPLES, LBLOWSNOW, &
GCOMPUTE_SRC, XRSNOW, &
LOCEAN, LDEEPOC, LDIAG_IN_RUN, &
diff --git a/src/mesonh/ext/prep_ideal_case.f90 b/src/mesonh/ext/prep_ideal_case.f90
index 7a1b8787e..4c75d2886 100644
--- a/src/mesonh/ext/prep_ideal_case.f90
+++ b/src/mesonh/ext/prep_ideal_case.f90
@@ -409,6 +409,9 @@ USE MODI_PGD_SURF_ATM
USE MODI_ICE_ADJUST_BIS
USE MODI_WRITE_PGD_SURF_ATM_n
USE MODI_PREP_SURF_MNH
+USE MODI_INIT_SALT
+USE MODI_AER2LIMA
+USE MODD_PARAM_LIMA
!
!JUAN
USE MODE_SPLITTINGZ_ll
@@ -716,6 +719,8 @@ IF (GFOUND) READ(UNIT=NLUPRE,NML=NAM_IBM_LSF )
CALL INI_FIELD_LIST(1)
!
CALL INI_FIELD_SCALARS()
+! Sea salt
+CALL INIT_SALT
!
IF( LEN_TRIM(CPGD_FILE) /= 0 ) THEN
! open the PGD_FILE
@@ -1734,6 +1739,9 @@ CALL IO_File_close(TZEXPREFILE) ! Close the EXPRE file
!
IF ( LCH_INIT_FIELD ) CALL CH_INIT_FIELD_n(1, NLUOUT, NVERB)
!
+! Initialization LIMA variables by ORILAM
+IF (CCLOUD == 'LIMA' .AND. ((LORILAM).OR.(LDUST).OR.(LSALT))) &
+ CALL AER2LIMA(XSVT, XRHODREF, XRT(:,:,:,1), XPABST, XTHT, XZZ)
!-------------------------------------------------------------------------------
!
!* 7. INITIALIZE LEVELSET FOR IBM
--
GitLab