From 3c0546fdb070c5cfe6d8d95dda698bc65ada6d65 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Beno=C3=AEt=20Vi=C3=A9?= <benoit.vie@meteo.fr>
Date: Wed, 30 Nov 2022 09:52:50 +0100
Subject: [PATCH] bugfix add checks for concentration (vs XCTMIN) in LIMA
 processes

---
 src/mesonh/micro/ini_lima_cold_mixed.f90         |  4 ++--
 src/mesonh/micro/ini_rain_ice.f90                |  6 +++---
 .../micro/lima_collisional_ice_breakup.f90       |  3 ++-
 .../micro/lima_conversion_melting_snow.f90       |  4 ++--
 src/mesonh/micro/lima_droplets_accretion.f90     |  5 +++--
 src/mesonh/micro/lima_droplets_hom_freezing.f90  |  2 +-
 src/mesonh/micro/lima_droplets_riming_snow.f90   |  9 +++++----
 src/mesonh/micro/lima_drops_self_collection.f90  |  4 ++--
 src/mesonh/micro/lima_graupel.f90                | 16 ++++++++--------
 src/mesonh/micro/lima_rain_accr_snow.f90         | 12 +++++++++---
 src/mesonh/micro/lima_rain_freezing.f90          |  5 +++--
 11 files changed, 40 insertions(+), 30 deletions(-)

diff --git a/src/mesonh/micro/ini_lima_cold_mixed.f90 b/src/mesonh/micro/ini_lima_cold_mixed.f90
index e8af15d58..2e3d956a4 100644
--- a/src/mesonh/micro/ini_lima_cold_mixed.f90
+++ b/src/mesonh/micro/ini_lima_cold_mixed.f90
@@ -343,8 +343,8 @@ IF (GFLAG) THEN
   WRITE(UNIT=ILUOUT0,FMT='(" XLBEXH =",E13.6," XLBH =",E13.6)') XLBEXH,XLBH
 END IF
 !
-XLBDAS_MAX = 500000. ! used only before transforming lambda for non MP PSD
-XLBDAS_MIN = 1000. *1.E-10
+XLBDAS_MAX = 1.E7 ! (eq to r~1E-7kg/kg) (for non MP PSD, use conversion XTRANS_MP_GAMMAS)
+XLBDAS_MIN = 1.   ! (eq to r~0.18kg/kg) (for non MP PSD, use conversion XTRANS_MP_GAMMAS)
 XLBDAG_MAX = 100000.0
 !
 ZCONC_MAX  = 1.E6 ! Maximal concentration for falling particules set to 1 per cc
diff --git a/src/mesonh/micro/ini_rain_ice.f90 b/src/mesonh/micro/ini_rain_ice.f90
index 38f8fed02..2f3b2b1c4 100644
--- a/src/mesonh/micro/ini_rain_ice.f90
+++ b/src/mesonh/micro/ini_rain_ice.f90
@@ -442,11 +442,11 @@ XLBH   = ( XAH*XCCH*MOMG(XALPHAH,XNUH,XBH) )**(-XLBEXH)
 XLBDAS_MAX = 100000.0
 !
 ZCONC_MAX  = 1.E6 ! Maximal concentration for falling particules set to 1 per cc
-#if defined(REPRO48) || defined(REPRO55)
 IF(XCCS>0. .AND. XCXS>0. )XLBDAS_MAX = ( ZCONC_MAX/XCCS )**(1./XCXS)
+#if defined(REPRO48) || defined(REPRO55)
 #else
-XLBDAS_MAX = 1.E6
-XLBDAS_MIN = 1000.
+IF (LSNOW_T) XLBDAS_MAX = 1.E6
+XLBDAS_MIN = 1.E-10
 #endif
 !
 XCONC_SEA=1E8 ! 100/cm3
diff --git a/src/mesonh/micro/lima_collisional_ice_breakup.f90 b/src/mesonh/micro/lima_collisional_ice_breakup.f90
index 3c22bf960..a6848d143 100644
--- a/src/mesonh/micro/lima_collisional_ice_breakup.f90
+++ b/src/mesonh/micro/lima_collisional_ice_breakup.f90
@@ -128,7 +128,8 @@ REAL                               :: ZFACT1_XNDEBRIS, ZFACT2_XNDEBRIS
 !
 !-------------------------------------------------------------------------------
 
-GCIBU(:) = LCIBU .AND. (PRST(:)>XRTMIN(5)) .AND. (PRGT(:)>XRTMIN(6)) .AND. LDCOMPUTE(:)
+GCIBU(:) = LCIBU        .AND. PRST(:)>XRTMIN(5) .AND. PRGT(:)>XRTMIN(6) .AND. &
+           LDCOMPUTE(:) .AND. PCST(:)>XCTMIN(5) .AND. PCGT(:)>XCTMIN(6)
 ICIBU    = COUNT( GCIBU(:) )
 !
 P_RI_CIBU(:)=0.
diff --git a/src/mesonh/micro/lima_conversion_melting_snow.f90 b/src/mesonh/micro/lima_conversion_melting_snow.f90
index 454df3eba..ef46c794f 100644
--- a/src/mesonh/micro/lima_conversion_melting_snow.f90
+++ b/src/mesonh/micro/lima_conversion_melting_snow.f90
@@ -64,7 +64,7 @@ END MODULE MODI_LIMA_CONVERSION_MELTING_SNOW
 !              ------------
 !
 USE MODD_CST,              ONLY : XTT, XMV, XMD, XLVTT, XCPV, XCL, XESTT, XRV
-USE MODD_PARAM_LIMA,       ONLY : XRTMIN, XNUS, XALPHAS
+USE MODD_PARAM_LIMA,       ONLY : XRTMIN, XCTMIN, XNUS, XALPHAS
 USE MODD_PARAM_LIMA_MIXED, ONLY : XFSCVMG
 USE MODD_PARAM_LIMA_COLD,  ONLY : X0DEPS, XEX0DEPS, X1DEPS, XEX1DEPS, XBS, XFVELOS
 !
@@ -104,7 +104,7 @@ P_RS_CMEL(:)=0.
 P_CS_CMEL(:)=0.
 !
 ZW(:) = 0.0
-WHERE( (PRST(:)>XRTMIN(5)) .AND. (PT(:)>XTT) .AND. LDCOMPUTE(:) )
+WHERE( PRST(:)>XRTMIN(5) .AND. PCST(:)>XCTMIN(5) .AND. PT(:)>XTT .AND. LDCOMPUTE(:) )
    ZW(:) = PRVT(:)*PPRES(:)/((XMV/XMD)+PRVT(:)) ! Vapor pressure
    ZW(:) = PKA(:)*(XTT-PT(:)) +                                 &
               ( PDV(:)*(XLVTT + ( XCPV - XCL ) * ( PT(:) - XTT )) &
diff --git a/src/mesonh/micro/lima_droplets_accretion.f90 b/src/mesonh/micro/lima_droplets_accretion.f90
index c7a4232fd..d183953cd 100644
--- a/src/mesonh/micro/lima_droplets_accretion.f90
+++ b/src/mesonh/micro/lima_droplets_accretion.f90
@@ -116,6 +116,7 @@ ZW4(:) = 0.0
 IF ( LKHKO ) THEN
 !
    GACCR(:) = PRRT(:)>XRTMIN(3) .AND. &
+              PCRT(:)>XCTMIN(3) .AND. &
               PRCT(:)>XRTMIN(2) .AND. &
               PCCT(:)>XCTMIN(2)
 !
@@ -130,8 +131,8 @@ IF ( LKHKO ) THEN
    END WHERE
 !
 ELSE IF (NMOM_C.EQ.1 .AND. NMOM_R.EQ.1) THEN
-   GACCR(:) = PRRT(:)>XRTMIN(3) .AND. &
-              PRCT(:)>XRTMIN(2)
+   GACCR(:) = PRRT(:)>XRTMIN(3) .AND. PCRT(:)>XCTMIN(3) .AND. &
+              PRCT(:)>XRTMIN(2) .AND. PCCT(:)>XCTMIN(2)
    WHERE ( GACCR(:) )
       P_RC_ACCR(:) = - XFCACCR * PRCT(:)    &
                    * PLBDR(:)**XEXCACCR     &
diff --git a/src/mesonh/micro/lima_droplets_hom_freezing.f90 b/src/mesonh/micro/lima_droplets_hom_freezing.f90
index 58e12b5f1..b33d7a350 100644
--- a/src/mesonh/micro/lima_droplets_hom_freezing.f90
+++ b/src/mesonh/micro/lima_droplets_hom_freezing.f90
@@ -98,7 +98,7 @@ P_TH_HONC(:) = 0.
 P_RC_HONC(:) = 0.
 P_CC_HONC(:) = 0.
 !
-WHERE ( (PT(:)<XTT-35.0) .AND. (PCCT(:)>XCTMIN(2)) .AND. (PRCT(:)>XRTMIN(2)) )
+WHERE ( PT(:)<XTT-35.0 .AND. PCCT(:)>XCTMIN(2) .AND. PRCT(:)>XRTMIN(2) )
    ZTCELSIUS(:) = PT(:)-XTT                                    ! T [°C]
    !
    ZZW(:) = 0.0
diff --git a/src/mesonh/micro/lima_droplets_riming_snow.f90 b/src/mesonh/micro/lima_droplets_riming_snow.f90
index b9d765365..cd4668238 100644
--- a/src/mesonh/micro/lima_droplets_riming_snow.f90
+++ b/src/mesonh/micro/lima_droplets_riming_snow.f90
@@ -76,7 +76,7 @@ END MODULE MODI_LIMA_DROPLETS_RIMING_SNOW
 !              ------------
 !
 USE MODD_CST,              ONLY : XTT
-USE MODD_PARAM_LIMA,       ONLY : XRTMIN, XCEXVT, XNUS, XALPHAS, LMURAKAMI
+USE MODD_PARAM_LIMA,       ONLY : XRTMIN, XCTMIN, XCEXVT, XNUS, XALPHAS, LMURAKAMI
 USE MODD_PARAM_LIMA_MIXED, ONLY : NGAMINC, XRIMINTP1, XRIMINTP2, XGAMINC_RIM1, XGAMINC_RIM2, XGAMINC_RIM4, &
                                   XCRIMSS, XEXCRIMSS, XSRIMCG, XEXSRIMCG, XSRIMCG2, XSRIMCG3, XEXSRIMCG2, &
                                   XHMLINTP1, XHMLINTP2, XGAMINC_HMC, XHM_FACTS, XHMTMIN, XHMTMAX
@@ -129,7 +129,8 @@ DO JI = 1, SIZE(PRCT)
 !*       Cloud droplet riming of the aggregates  
 !        --------------------------------------
 !
-   IF ( PRCT(JI)>XRTMIN(2) .AND. PRST(JI)>XRTMIN(5) .AND. PT(JI)<XTT .AND. LDCOMPUTE(JI) ) THEN
+   IF ( PRCT(JI)>XRTMIN(2) .AND. PRST(JI)>XRTMIN(5) .AND. PT(JI)<XTT .AND. &
+        PCCT(JI)>XCTMIN(2) .AND. PCST(JI)>XCTMIN(5) .AND. LDCOMPUTE(JI) ) THEN
 !
       ZVEC1(JI) = PLBDS(JI)
       ZVEC1W(JI)= ( XFVELOS**XALPHAS + PLBDS(JI)**XALPHAS ) ** (1./XALPHAS) ! modified equivalent lambda
@@ -200,8 +201,8 @@ DO JI = 1, SIZE(PRCT)
 !*       Hallett-Mossop ice production (HMS)  
 !        -----------------------------------
 !
-   IF ( PT(JI)<XHMTMAX .AND. PT(JI)>XHMTMIN .AND.                      &
-        PRST(JI)>XRTMIN(5) .AND. PRCT(JI)>XRTMIN(2) .AND. LDCOMPUTE(JI) ) THEN
+   IF ( PRST(JI)>XRTMIN(5) .AND. PRCT(JI)>XRTMIN(2) .AND. PT(JI)<XHMTMAX .AND. PT(JI)>XHMTMIN .AND. &
+        PCST(JI)>XCTMIN(5) .AND. PCCT(JI)>XCTMIN(2) .AND. LDCOMPUTE(JI) ) THEN
 !
       ZVEC1(JI) = PLBDC(JI)
       ZVEC2(JI) = MAX( 1.0001, MIN( REAL(NGAMINC)-0.0001,           &
diff --git a/src/mesonh/micro/lima_drops_self_collection.f90 b/src/mesonh/micro/lima_drops_self_collection.f90
index 042cde084..3f064dfcd 100644
--- a/src/mesonh/micro/lima_drops_self_collection.f90
+++ b/src/mesonh/micro/lima_drops_self_collection.f90
@@ -105,12 +105,12 @@ ZW1(:) = 0.0
 ZW2(:) = 0.0
 ZW3(:) = 0.0
 !
-WHERE (PCRT(:)>XCTMIN(3) .AND. (ZW4(:)>1.E-4) .AND. LDCOMPUTE(:))  ! analytical integration
+WHERE ( PCRT(:)>XCTMIN(3) .AND. ZW4(:)>1.E-4 .AND. LDCOMPUTE(:))  ! analytical integration
    ZW1(:) = XSCBU2 * PCRT(:)**2 / PLBDR3(:)                        ! D>100 10-6 m
    ZW3(:) = ZW1(:)*ZSCBU(:)
 END WHERE
 !
-WHERE (PCRT(:)>XCTMIN(3) .AND. (ZW4(:)<=1.E-4) .AND. LDCOMPUTE(:))
+WHERE ( PCRT(:)>XCTMIN(3) .AND. ZW4(:)<=1.E-4 .AND. LDCOMPUTE(:))
    ZW2(:) = XSCBU3 *(PCRT(:) / PLBDR3(:))**2                       ! D<100 10-6 m
    ZW3(:) = ZW2(:)
 END WHERE
diff --git a/src/mesonh/micro/lima_graupel.f90 b/src/mesonh/micro/lima_graupel.f90
index 00ad1646a..8c96d2e09 100644
--- a/src/mesonh/micro/lima_graupel.f90
+++ b/src/mesonh/micro/lima_graupel.f90
@@ -146,7 +146,7 @@ END MODULE MODI_LIMA_GRAUPEL
 !              ------------
 !
 USE MODD_CST,              ONLY : XTT, XMD, XMV, XRD, XRV, XLVTT, XLMTT, XESTT, XCL, XCI, XCPV
-USE MODD_PARAM_LIMA,       ONLY : XRTMIN, XCEXVT, LHAIL
+USE MODD_PARAM_LIMA,       ONLY : XRTMIN, XCTMIN, XCEXVT, LHAIL
 USE MODD_PARAM_LIMA_MIXED, ONLY : XCXG, XDG, X0DEPG, X1DEPG, NGAMINC,                             &
                                   XFCDRYG, XFIDRYG, XCOLIG, XCOLSG, XCOLEXIG, XCOLEXSG,           &
                                   XFSDRYG, XLBSDRYG1, XLBSDRYG2, XLBSDRYG3, XKER_SDRYG,           &
@@ -323,7 +323,7 @@ END WHERE
 !*           1.b Collection of rs in the dry mode
 !            ------------------------------------
 !
-GDRY(:) = PRST(:)>XRTMIN(5) .AND. PRGT(:)>XRTMIN(6) .AND. LDCOMPUTE(:)
+GDRY(:) = PRST(:)>XRTMIN(5) .AND. PCST(:)>XCTMIN(5) .AND. PRGT(:)>XRTMIN(6) .AND. PCGT(:)>XCTMIN(6) .AND. LDCOMPUTE(:)
 !
 WHERE( GDRY )
 !
@@ -389,7 +389,7 @@ END WHERE
 !*           1.c  Collection of rr in the dry mode
 !            -------------------------------------
 !
-GDRY(:) = PRRT(:)>XRTMIN(3) .AND. PRGT(:)>XRTMIN(6) .AND. LDCOMPUTE(:)
+GDRY(:) = PRRT(:)>XRTMIN(3) .AND. PCRT(:)>XCTMIN(3) .AND. PRGT(:)>XRTMIN(6) .AND. PCGT(:)>XCTMIN(6) .AND. LDCOMPUTE(:)
 !
 WHERE( GDRY )
 !
@@ -462,7 +462,7 @@ ZRDRYG(:) = ZZW1(:) + ZZW2(:) + ZZW3(:) + ZZW4(:)
 !            ------------------------------
 !
 ZZW(:) = 0.0
-WHERE( PRGT(:)>XRTMIN(6) .AND. LDCOMPUTE(:) )
+WHERE( PRGT(:)>XRTMIN(6) .AND. PCGT(:)>XCTMIN(6) .AND. LDCOMPUTE(:) )
    ZZW5(:) = ZZW2(:) / (XCOLIG*EXP(XCOLEXIG*(PT(:)-XTT)) ) ! RIWETG
    ZZW6(:) = ZZW3(:) / (XCOLSG*EXP(XCOLEXSG*(PT(:)-XTT)) ) ! RSWETG
    ZZW6N(:)= ZZW3N(:)/ (XCOLSG*EXP(XCOLEXSG*(PT(:)-XTT)) ) ! NSWETG
@@ -488,7 +488,7 @@ END WHERE
 ZZW(:) = 0.0
 NHAIL = 0.
 IF (LHAIL) NHAIL = 1. 
-WHERE( LDCOMPUTE(:) .AND. PRGT(:)>XRTMIN(6) .AND. PT(:)<XTT .AND. &
+WHERE( LDCOMPUTE(:) .AND. PRGT(:)>XRTMIN(6) .AND. PCGT(:)>XCTMIN(6) .AND. PT(:)<XTT .AND. &
        (ZRDRYG(:)-ZZW2(:)-ZZW3(:))>=(ZRWETG(:)-ZZW5(:)-ZZW6(:)) .AND. ZRWETG(:)-ZZW5(:)-ZZW6(:)>0.0 ) 
 !
 ! Mass of rain and cloud droplets frozen by graupel in wet mode : RCWETG + RRWETG = RWETG - RIWETG - RSWETG
@@ -517,7 +517,7 @@ END WHERE
 !            1.g Dry mode
 !            ------------
 !
-WHERE( LDCOMPUTE(:) .AND. PRGT(:)>XRTMIN(6) .AND. PT(:)<XTT .AND.                  &
+WHERE( LDCOMPUTE(:) .AND. PRGT(:)>XRTMIN(6) .AND. PCGT(:)>XCTMIN(6) .AND. PT(:)<XTT .AND.                  &
        (ZRDRYG(:)-ZZW2(:)-ZZW3(:))<(ZRWETG(:)-ZZW5(:)-ZZW6(:)) .AND. ZRDRYG(:)>0.0 )
    !
    P_RC_DRYG(:) = - ZZW1(:)
@@ -540,7 +540,7 @@ END WHERE
 ! BVIE test ZRDRYG<ZZW ?????????????????????????
 !GDRY(:) = (PT(:)<XHMTMAX) .AND. (PT(:)>XHMTMIN)    .AND. (ZRDRYG(:)<ZZW(:))&
 GDRY(:) = PT(:)<XHMTMAX .AND. PT(:)>XHMTMIN .AND. PRGT(:)>XRTMIN(6) .AND. PRCT(:)>XRTMIN(2) .AND. LDCOMPUTE(:) .AND. &
-          (ZRDRYG(:)-ZZW2(:)-ZZW3(:))<(ZRWETG(:)-ZZW5(:)-ZZW6(:))
+          PCGT(:)>XCTMIN(6) .AND. PCCT(:)>XCTMIN(2) .AND. (ZRDRYG(:)-ZZW2(:)-ZZW3(:))<(ZRWETG(:)-ZZW5(:)-ZZW6(:))
 
 ZZX(:)=9999.
 ZVEC1(:)=0.
@@ -571,7 +571,7 @@ END WHERE
 !        -------------------
 !
 ZZX(:) = 0.0
-WHERE( (PRGT(:)>XRTMIN(6)) .AND. (PT(:)>XTT) .AND. LDCOMPUTE(:) )
+WHERE( PRGT(:)>XRTMIN(6) .AND. PCGT(:)>XCTMIN(6) .AND. PT(:)>XTT .AND. LDCOMPUTE(:) )
    ZZX(:) = PRVT(:)*PPRES(:)/((XMV/XMD)+PRVT(:)) ! Vapor pressure
    ZZX(:) = PKA(:)*(XTT-PT(:)) +                                 &
               ( PDV(:)*(XLVTT + ( XCPV - XCL ) * ( PT(:) - XTT )) &
diff --git a/src/mesonh/micro/lima_rain_accr_snow.f90 b/src/mesonh/micro/lima_rain_accr_snow.f90
index 803f123b2..a63ac24a4 100644
--- a/src/mesonh/micro/lima_rain_accr_snow.f90
+++ b/src/mesonh/micro/lima_rain_accr_snow.f90
@@ -70,7 +70,7 @@ END MODULE MODI_LIMA_RAIN_ACCR_SNOW
 !              ------------
 !
 USE MODD_CST,              ONLY : XTT
-USE MODD_PARAM_LIMA,       ONLY : XRTMIN, XCEXVT
+USE MODD_PARAM_LIMA,       ONLY : XRTMIN, XCTMIN, XCEXVT
 USE MODD_PARAM_LIMA_WARM,  ONLY : XBR
 USE MODD_PARAM_LIMA_COLD,  ONLY : XBS, XTRANS_MP_GAMMAS
 USE MODD_PARAM_LIMA_MIXED, ONLY : NACCLBDAS, XACCINTP1S, XACCINTP2S,             &
@@ -134,6 +134,12 @@ ZZW3(:) = 0.
 ZZW4(:) = 0.
 ZZW5(:) = 0.
 !
+ZZWC1(:) = 0.
+ZZWC2(:) = 0.
+ZZWC3(:) = 0.
+ZZWC4(:) = 0.
+ZZWC5(:) = 0.
+!
 IVEC1(:) = 0
 IVEC2(:) = 0
 ZVEC1(:) = 0.
@@ -145,7 +151,8 @@ ZVEC3(:) = 0.
 !
 !
 GACC(:) = .False.
-GACC(:) = (PRRT(:)>XRTMIN(3)) .AND. (PRST(:)>XRTMIN(5)) .AND. (PT(:)<XTT) .AND. LDCOMPUTE(:)
+GACC(:) = (PRRT(:)>XRTMIN(3)) .AND. (PRST(:)>XRTMIN(5)) .AND. (PT(:)<XTT) .AND. LDCOMPUTE(:) .AND. &
+          (PCRT(:)>XCTMIN(3)) .AND. (PCST(:)>XCTMIN(5))
 !
 WHERE( GACC )
 !
@@ -303,7 +310,6 @@ WHERE( GACC )
 !
 END WHERE
 !
-!
 !-------------------------------------------------------------------------------
 !
 CONTAINS
diff --git a/src/mesonh/micro/lima_rain_freezing.f90 b/src/mesonh/micro/lima_rain_freezing.f90
index a5a9225bc..a6c9504a1 100644
--- a/src/mesonh/micro/lima_rain_freezing.f90
+++ b/src/mesonh/micro/lima_rain_freezing.f90
@@ -65,7 +65,7 @@ END MODULE MODI_LIMA_RAIN_FREEZING
 !              ------------
 !
 USE MODD_CST,              ONLY : XTT
-USE MODD_PARAM_LIMA,       ONLY : XRTMIN, XCEXVT
+USE MODD_PARAM_LIMA,       ONLY : XRTMIN, XCTMIN, XCEXVT
 USE MODD_PARAM_LIMA_MIXED, ONLY : XICFRR, XEXICFRR, XRCFRI, XEXRCFRI
 !
 IMPLICIT NONE
@@ -111,7 +111,8 @@ P_CI_CFRZ(:)=0.
 ZW1(:)=0.
 ZW2(:)=0.
 !
-WHERE( (PRIT(:)>XRTMIN(4)) .AND. (PRRT(:)>XRTMIN(3)) .AND. (PT(:)<XTT) .AND. LDCOMPUTE(:) )
+WHERE( PRIT(:)>XRTMIN(4) .AND. PRRT(:)>XRTMIN(3) .AND. PT(:)<XTT .AND. &
+       PCIT(:)>XCTMIN(4) .AND. PCRT(:)>XCTMIN(3) .AND. LDCOMPUTE(:) )
 !
    ZW1(:) = XICFRR * PRIT(:) * PCRT(:)                    & ! RICFRRG
                                      * PLBDR(:)**XEXICFRR         &
-- 
GitLab