From cf41ae7e57865a34283711bb9ad2bee28c355d6e Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Beno=C3=AEt=20Vi=C3=A9?= <benoit.vie@meteo.fr>
Date: Wed, 27 Apr 2022 15:20:48 +0200
Subject: [PATCH] Adapt CIBU to new snow characteristics Add CIBU and RDSF to
LIMA_SPLIT Move CIBU and RDSF from modd_param_lima_cold to
modd_param_lima_mixed
---
src/MNH/ini_lima_cold_mixed.f90 | 4 +-
src/MNH/lima.f90 | 28 ++
src/MNH/lima_collisional_ice_breakup.f90 | 388 ++++++++++++++++++
src/MNH/lima_mixed_fast_processes.f90 | 22 +-
src/MNH/lima_raindrop_shattering_freezing.f90 | 152 +++++++
src/MNH/lima_tendencies.f90 | 58 ++-
src/MNH/modd_param_lima_cold.f90 | 34 --
src/MNH/modd_param_lima_mixed.f90 | 42 +-
8 files changed, 679 insertions(+), 49 deletions(-)
create mode 100644 src/MNH/lima_collisional_ice_breakup.f90
create mode 100644 src/MNH/lima_raindrop_shattering_freezing.f90
diff --git a/src/MNH/ini_lima_cold_mixed.f90 b/src/MNH/ini_lima_cold_mixed.f90
index 23b533c92..213425be2 100644
--- a/src/MNH/ini_lima_cold_mixed.f90
+++ b/src/MNH/ini_lima_cold_mixed.f90
@@ -986,13 +986,13 @@ XCIBUINTP_G = XALPHAG / LOG(ZRATE_G)
XCIBUINTP1_G = 1.0 + XCIBUINTP_G * LOG(XDCGLIM_CIBU_MIN/(XGAMINC_BOUND_CIBU_GMIN)**(1.0/XALPHAG))
!
! For ZNI_CIBU
-XFACTOR_CIBU_NI = (XPI / 4.0) * XCCG * XCCS * (ZRHO00**XCEXVT)
+XFACTOR_CIBU_NI = XLBS * (XPI / 4.0) * XCCG * (ZRHO00**XCEXVT)
XMOMGG_CIBU_1 = MOMG(XALPHAG,XNUG,2.0+XDG)
XMOMGG_CIBU_2 = MOMG(XALPHAG,XNUG,2.0)
XMOMGS_CIBU_1 = MOMG(XALPHAS,XNUS,XDS)
!
! For ZRI_CIBU
-XFACTOR_CIBU_RI = XAS * (XPI / 4.0) * XCCG * XCCS * (ZRHO00**XCEXVT)
+XFACTOR_CIBU_RI = XLBS * XAS * (XPI / 4.0) * XCCG * (ZRHO00**XCEXVT)
XMOMGS_CIBU_2 = MOMG(XALPHAS,XNUS,XBS)
XMOMGS_CIBU_3 = MOMG(XALPHAS,XNUS,XBS+XDS)
!
diff --git a/src/MNH/lima.f90 b/src/MNH/lima.f90
index 02b7a5b9c..330a5a3f5 100644
--- a/src/MNH/lima.f90
+++ b/src/MNH/lima.f90
@@ -246,6 +246,8 @@ REAL, DIMENSION(:), ALLOCATABLE :: &
Z_TH_ACC, Z_RR_ACC, Z_CR_ACC, Z_RS_ACC, Z_RG_ACC, & ! rain accretion on aggregates (ACC) : rr, Nr, rs, rg, th
Z_RS_CMEL, & ! conversion-melting (CMEL) : rs, rg=-rs
Z_TH_CFRZ, Z_RR_CFRZ, Z_CR_CFRZ, Z_RI_CFRZ, Z_CI_CFRZ, & ! rain freezing (CFRZ) : rr, Nr, ri, Ni, rg=-rr-ri, th
+ Z_RI_CIBU, Z_CI_CIBU, & ! collisional ice break-up (CIBU) : ri, Ni, rs=-ri
+ Z_RI_RDSF, Z_CI_RDSF, & ! rain drops freezing shattering (RDSF) : ri, Ni, rg=-ri
Z_TH_WETG, Z_RC_WETG, Z_CC_WETG, Z_RR_WETG, Z_CR_WETG, & ! wet growth of graupel (WETG) : rc, NC, rr, Nr, ri, Ni, rs, rg, rh, th
Z_RI_WETG, Z_CI_WETG, Z_RS_WETG, Z_RG_WETG, Z_RH_WETG, & ! wet growth of graupel (WETG) : rc, NC, rr, Nr, ri, Ni, rs, rg, rh, th
Z_TH_DRYG, Z_RC_DRYG, Z_CC_DRYG, Z_RR_DRYG, Z_CR_DRYG, & ! dry growth of graupel (DRYG) : rc, Nc, rr, Nr, ri, Ni, rs, rg, th
@@ -297,6 +299,8 @@ REAL, DIMENSION(:,:,:), ALLOCATABLE :: &
ZTOT_TH_ACC, ZTOT_RR_ACC, ZTOT_CR_ACC, ZTOT_RS_ACC, ZTOT_RG_ACC, & ! rain accretion on aggregates (ACC)
ZTOT_RS_CMEL, & ! conversion-melting (CMEL)
ZTOT_TH_CFRZ, ZTOT_RR_CFRZ, ZTOT_CR_CFRZ, ZTOT_RI_CFRZ, ZTOT_CI_CFRZ, & ! rain freezing (CFRZ)
+ ZTOT_RI_CIBU, ZTOT_CI_CIBU, & ! collisional ice break-up (CIBU)
+ ZTOT_RI_RDSF, ZTOT_CI_RDSF, & ! rain drops freezing shattering (RDSF)
ZTOT_TH_WETG, ZTOT_RC_WETG, ZTOT_CC_WETG, ZTOT_RR_WETG, ZTOT_CR_WETG, & ! wet growth of graupel (WETG)
ZTOT_RI_WETG, ZTOT_CI_WETG, ZTOT_RS_WETG, ZTOT_RG_WETG, ZTOT_RH_WETG, & ! wet growth of graupel (WETG)
ZTOT_TH_DRYG, ZTOT_RC_DRYG, ZTOT_CC_DRYG, ZTOT_RR_DRYG, ZTOT_CR_DRYG, & ! dry growth of graupel (DRYG)
@@ -462,6 +466,10 @@ if ( lbu_enable ) then
allocate( ZTOT_CR_CFRZ (size( ptht, 1), size( ptht, 2), size( ptht, 3) ) ); ZTOT_CR_CFRZ(:,:,:) = 0.
allocate( ZTOT_RI_CFRZ (size( ptht, 1), size( ptht, 2), size( ptht, 3) ) ); ZTOT_RI_CFRZ(:,:,:) = 0.
allocate( ZTOT_CI_CFRZ (size( ptht, 1), size( ptht, 2), size( ptht, 3) ) ); ZTOT_CI_CFRZ(:,:,:) = 0.
+ allocate( ZTOT_RI_CIBU (size( ptht, 1), size( ptht, 2), size( ptht, 3) ) ); ZTOT_RI_CIBU(:,:,:) = 0.
+ allocate( ZTOT_CI_CIBU (size( ptht, 1), size( ptht, 2), size( ptht, 3) ) ); ZTOT_CI_CIBU(:,:,:) = 0.
+ allocate( ZTOT_RI_RDSF (size( ptht, 1), size( ptht, 2), size( ptht, 3) ) ); ZTOT_RI_RDSF(:,:,:) = 0.
+ allocate( ZTOT_CI_RDSF (size( ptht, 1), size( ptht, 2), size( ptht, 3) ) ); ZTOT_CI_RDSF(:,:,:) = 0.
allocate( ZTOT_TH_WETG (size( ptht, 1), size( ptht, 2), size( ptht, 3) ) ); ZTOT_TH_WETG(:,:,:) = 0.
allocate( ZTOT_RC_WETG (size( ptht, 1), size( ptht, 2), size( ptht, 3) ) ); ZTOT_RC_WETG(:,:,:) = 0.
allocate( ZTOT_CC_WETG (size( ptht, 1), size( ptht, 2), size( ptht, 3) ) ); ZTOT_CC_WETG(:,:,:) = 0.
@@ -1022,6 +1030,10 @@ DO WHILE(ANY(ZTIME(IIB:IIE,IJB:IJE,IKTB:IKTE)<PTSTEP))
ALLOCATE(Z_CR_CFRZ(IPACK)) ; Z_CR_CFRZ = 0.
ALLOCATE(Z_RI_CFRZ(IPACK)) ; Z_RI_CFRZ = 0.
ALLOCATE(Z_CI_CFRZ(IPACK)) ; Z_CI_CFRZ = 0.
+ ALLOCATE(Z_RI_CIBU(IPACK)) ; Z_RI_CIBU = 0.
+ ALLOCATE(Z_CI_CIBU(IPACK)) ; Z_CI_CIBU = 0.
+ ALLOCATE(Z_RI_RDSF(IPACK)) ; Z_RI_RDSF = 0.
+ ALLOCATE(Z_CI_RDSF(IPACK)) ; Z_CI_RDSF = 0.
ALLOCATE(Z_TH_WETG(IPACK)) ; Z_TH_WETG = 0.
ALLOCATE(Z_RC_WETG(IPACK)) ; Z_RC_WETG = 0.
ALLOCATE(Z_CC_WETG(IPACK)) ; Z_CC_WETG = 0.
@@ -1094,6 +1106,8 @@ DO WHILE(ANY(ZTIME(IIB:IIE,IJB:IJE,IKTB:IKTE)<PTSTEP))
Z_TH_ACC, Z_RR_ACC, Z_CR_ACC, Z_RS_ACC, Z_RG_ACC, &
Z_RS_CMEL, &
Z_TH_CFRZ, Z_RR_CFRZ, Z_CR_CFRZ, Z_RI_CFRZ, Z_CI_CFRZ, &
+ Z_RI_CIBU, Z_CI_CIBU, &
+ Z_RI_RDSF, Z_CI_RDSF, &
Z_TH_WETG, Z_RC_WETG, Z_CC_WETG, Z_RR_WETG, Z_CR_WETG, &
Z_RI_WETG, Z_CI_WETG, Z_RS_WETG, Z_RG_WETG, Z_RH_WETG, &
Z_TH_DRYG, Z_RC_DRYG, Z_CC_DRYG, Z_RR_DRYG, Z_CR_DRYG, &
@@ -1399,6 +1413,10 @@ DO WHILE(ANY(ZTIME(IIB:IIE,IJB:IJE,IKTB:IKTE)<PTSTEP))
ZTOT_CR_CFRZ(I1(II),I2(II),I3(II)) = ZTOT_CR_CFRZ(I1(II),I2(II),I3(II)) + Z_CR_CFRZ(II) * ZMAXTIME(II)
ZTOT_RI_CFRZ(I1(II),I2(II),I3(II)) = ZTOT_RI_CFRZ(I1(II),I2(II),I3(II)) + Z_RI_CFRZ(II) * ZMAXTIME(II)
ZTOT_CI_CFRZ(I1(II),I2(II),I3(II)) = ZTOT_CI_CFRZ(I1(II),I2(II),I3(II)) + Z_CI_CFRZ(II) * ZMAXTIME(II)
+ ZTOT_RI_CIBU(I1(II),I2(II),I3(II)) = ZTOT_RI_CIBU(I1(II),I2(II),I3(II)) + Z_RI_CIBU(II) * ZMAXTIME(II)
+ ZTOT_CI_CIBU(I1(II),I2(II),I3(II)) = ZTOT_CI_CIBU(I1(II),I2(II),I3(II)) + Z_CI_CIBU(II) * ZMAXTIME(II)
+ ZTOT_RI_RDSF(I1(II),I2(II),I3(II)) = ZTOT_RI_RDSF(I1(II),I2(II),I3(II)) + Z_RI_RDSF(II) * ZMAXTIME(II)
+ ZTOT_CI_RDSF(I1(II),I2(II),I3(II)) = ZTOT_CI_RDSF(I1(II),I2(II),I3(II)) + Z_CI_RDSF(II) * ZMAXTIME(II)
ZTOT_TH_WETG(I1(II),I2(II),I3(II)) = ZTOT_TH_WETG(I1(II),I2(II),I3(II)) + Z_TH_WETG(II) * ZMAXTIME(II)
ZTOT_RC_WETG(I1(II),I2(II),I3(II)) = ZTOT_RC_WETG(I1(II),I2(II),I3(II)) + Z_RC_WETG(II) * ZMAXTIME(II)
ZTOT_CC_WETG(I1(II),I2(II),I3(II)) = ZTOT_CC_WETG(I1(II),I2(II),I3(II)) + Z_CC_WETG(II) * ZMAXTIME(II)
@@ -1566,6 +1584,10 @@ DO WHILE(ANY(ZTIME(IIB:IIE,IJB:IJE,IKTB:IKTE)<PTSTEP))
DEALLOCATE(Z_CR_CFRZ)
DEALLOCATE(Z_RI_CFRZ)
DEALLOCATE(Z_CI_CFRZ)
+ DEALLOCATE(Z_RI_CIBU)
+ DEALLOCATE(Z_CI_CIBU)
+ DEALLOCATE(Z_RI_RDSF)
+ DEALLOCATE(Z_CI_RDSF)
DEALLOCATE(Z_TH_WETG)
DEALLOCATE(Z_RC_WETG)
DEALLOCATE(Z_CC_WETG)
@@ -1702,6 +1724,8 @@ if ( lbu_enable ) then
call Budget_store_add( tbudgets(NBUDGET_RI), 'HMS', ztot_ri_hms (:, :, :) * zrhodjontstep(:, :, :) )
call Budget_store_add( tbudgets(NBUDGET_RI), 'CFRZ', ztot_ri_cfrz (:, :, :) * zrhodjontstep(:, :, :) )
call Budget_store_add( tbudgets(NBUDGET_RI), 'DEPI', ztot_ri_depi (:, :, :) * zrhodjontstep(:, :, :) )
+ call Budget_store_add( tbudgets(NBUDGET_RI), 'CIBU', ztot_ri_cibu (:, :, :) * zrhodjontstep(:, :, :) )
+ call Budget_store_add( tbudgets(NBUDGET_RI), 'RDSF', ztot_ri_rdsf (:, :, :) * zrhodjontstep(:, :, :) )
call Budget_store_add( tbudgets(NBUDGET_RI), 'WETG', ztot_ri_wetg (:, :, :) * zrhodjontstep(:, :, :) )
call Budget_store_add( tbudgets(NBUDGET_RI), 'DRYG', ztot_ri_dryg (:, :, :) * zrhodjontstep(:, :, :) )
call Budget_store_add( tbudgets(NBUDGET_RI), 'HMG', ztot_ri_hmg (:, :, :) * zrhodjontstep(:, :, :) )
@@ -1717,6 +1741,7 @@ if ( lbu_enable ) then
call Budget_store_add( tbudgets(NBUDGET_RS), 'HMS', ztot_rs_hms (:, :, :) * zrhodjontstep(:, :, :) )
call Budget_store_add( tbudgets(NBUDGET_RS), 'ACC', ztot_rs_acc (:, :, :) * zrhodjontstep(:, :, :) )
call Budget_store_add( tbudgets(NBUDGET_RS), 'CMEL', ztot_rs_cmel(:, :, :) * zrhodjontstep(:, :, :) )
+ call Budget_store_add( tbudgets(NBUDGET_RS), 'CIBU', -ztot_ri_cibu(:, :, :) * zrhodjontstep(:, :, :) )
call Budget_store_add( tbudgets(NBUDGET_RS), 'WETG', ztot_rs_wetg(:, :, :) * zrhodjontstep(:, :, :) )
call Budget_store_add( tbudgets(NBUDGET_RS), 'DRYG', ztot_rs_dryg(:, :, :) * zrhodjontstep(:, :, :) )
end if
@@ -1729,6 +1754,7 @@ if ( lbu_enable ) then
call Budget_store_add( tbudgets(NBUDGET_RG), 'CMEL', -ztot_rs_cmel(:, :, :) * zrhodjontstep(:, :, :) )
call Budget_store_add( tbudgets(NBUDGET_RG), 'CFRZ', ( -ztot_rr_cfrz(:, :, :) - ztot_ri_cfrz(:, :, :) ) &
* zrhodjontstep(:, :, :) )
+ call Budget_store_add( tbudgets(NBUDGET_RG), 'RDSF', -ztot_ri_rdsf(:, :, :) * zrhodjontstep(:, :, :) )
call Budget_store_add( tbudgets(NBUDGET_RG), 'WETG', ztot_rg_wetg(:, :, :) * zrhodjontstep(:, :, :) )
call Budget_store_add( tbudgets(NBUDGET_RG), 'DRYG', ztot_rg_dryg(:, :, :) * zrhodjontstep(:, :, :) )
call Budget_store_add( tbudgets(NBUDGET_RG), 'HMG', ztot_rg_hmg (:, :, :) * zrhodjontstep(:, :, :) )
@@ -1782,6 +1808,8 @@ if ( lbu_enable ) then
call Budget_store_add( tbudgets(idx), 'IMLT', -ztot_cc_imlt (:, :, :) * zrhodjontstep(:, :, :) )
call Budget_store_add( tbudgets(idx), 'HMS', ztot_ci_hms (:, :, :) * zrhodjontstep(:, :, :) )
call Budget_store_add( tbudgets(idx), 'CFRZ', ztot_ci_cfrz (:, :, :) * zrhodjontstep(:, :, :) )
+ call Budget_store_add( tbudgets(idx), 'CIBU', ztot_ci_cibu (:, :, :) * zrhodjontstep(:, :, :) )
+ call Budget_store_add( tbudgets(idx), 'RDSF', ztot_ci_rdsf (:, :, :) * zrhodjontstep(:, :, :) )
call Budget_store_add( tbudgets(idx), 'WETG', ztot_ci_wetg (:, :, :) * zrhodjontstep(:, :, :) )
call Budget_store_add( tbudgets(idx), 'DRYG', ztot_ci_dryg (:, :, :) * zrhodjontstep(:, :, :) )
call Budget_store_add( tbudgets(idx), 'HMG', ztot_ci_hmg (:, :, :) * zrhodjontstep(:, :, :) )
diff --git a/src/MNH/lima_collisional_ice_breakup.f90 b/src/MNH/lima_collisional_ice_breakup.f90
new file mode 100644
index 000000000..ffcb69c7c
--- /dev/null
+++ b/src/MNH/lima_collisional_ice_breakup.f90
@@ -0,0 +1,388 @@
+!MNH_LIC Copyright 2018-2021 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
+!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
+!MNH_LIC for details. version 1.
+!-------------------------------------------------------------------------------
+! ########################################
+ MODULE MODI_LIMA_COLLISIONAL_ICE_BREAKUP
+! ########################################
+!
+INTERFACE
+ SUBROUTINE LIMA_COLLISIONAL_ICE_BREAKUP (LDCOMPUTE, &
+ PRHODREF, &
+ PRIT, PRST, PRGT, PCIT, &
+ PLBDS, PLBDG, &
+ P_RI_CIBU, P_CI_CIBU )
+!
+LOGICAL, DIMENSION(:),INTENT(IN) :: LDCOMPUTE
+!
+REAL, DIMENSION(:), INTENT(IN) :: PRHODREF
+!
+REAL, DIMENSION(:), INTENT(IN) :: PRIT
+REAL, DIMENSION(:), INTENT(IN) :: PRST
+REAL, DIMENSION(:), INTENT(IN) :: PRGT
+REAL, DIMENSION(:), INTENT(IN) :: PCIT
+REAL, DIMENSION(:), INTENT(IN) :: PLBDS
+REAL, DIMENSION(:), INTENT(IN) :: PLBDG
+!
+REAL, DIMENSION(:), INTENT(OUT) :: P_RI_CIBU
+REAL, DIMENSION(:), INTENT(OUT) :: P_CI_CIBU
+!
+END SUBROUTINE LIMA_COLLISIONAL_ICE_BREAKUP
+END INTERFACE
+END MODULE MODI_LIMA_COLLISIONAL_ICE_BREAKUP
+!
+! #######################################################################
+ SUBROUTINE LIMA_COLLISIONAL_ICE_BREAKUP (LDCOMPUTE, &
+ PRHODREF, &
+ PRIT, PRST, PRGT, PCIT, &
+ PLBDS, PLBDG, &
+ P_RI_CIBU, P_CI_CIBU )
+! #######################################################################
+!
+!! PURPOSE
+!! -------
+!! Compute the collisional ice break-up (secondary ice production process)
+!!
+!! AUTHOR
+!! ------
+!! J.-P. Pinty * Laboratoire d'Aerologie*
+!! T. Hoarau * Laboratoire d'Aerologie*
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 04/2022 duplicate from original for LIMA_SPLIT
+! B. Vié 04/2022 Adapt to the new snow characteristics
+!
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODD_PARAM_LIMA, ONLY : LCIBU, XRTMIN, XCTMIN, XCEXVT, XALPHAS, XNUS, XNDEBRIS_CIBU
+
+USE MODD_PARAM_LIMA_COLD, ONLY : XBS, XCS, XDS, XFVELOS, XMNU0
+USE MODD_PARAM_LIMA_MIXED, ONLY : XCG, XDG, XCXG, &
+ XCIBUINTP_S, XCIBUINTP1_S, XCIBUINTP2_S, &
+ XCIBUINTP_G, XCIBUINTP1_G, &
+ XFACTOR_CIBU_NI, XFACTOR_CIBU_RI, &
+ XMOMGG_CIBU_1, XMOMGG_CIBU_2, &
+ XMOMGS_CIBU_1, XMOMGS_CIBU_2, XMOMGS_CIBU_3, &
+ NGAMINC, XGAMINC_CIBU_S, XGAMINC_CIBU_G
+!
+IMPLICIT NONE
+!
+!* 0.1 Declarations of dummy arguments :
+!
+!
+LOGICAL, DIMENSION(:),INTENT(IN) :: LDCOMPUTE
+!
+REAL, DIMENSION(:), INTENT(IN) :: PRHODREF
+!
+REAL, DIMENSION(:), INTENT(IN) :: PRIT
+REAL, DIMENSION(:), INTENT(IN) :: PRST
+REAL, DIMENSION(:), INTENT(IN) :: PRGT
+REAL, DIMENSION(:), INTENT(IN) :: PCIT
+REAL, DIMENSION(:), INTENT(IN) :: PLBDS
+REAL, DIMENSION(:), INTENT(IN) :: PLBDG
+!
+REAL, DIMENSION(:), INTENT(OUT) :: P_RI_CIBU
+REAL, DIMENSION(:), INTENT(OUT) :: P_CI_CIBU
+!
+!
+!* 0.2 Declarations of local variables :
+!
+LOGICAL, DIMENSION(SIZE(PRST)) :: GCIBU ! Test where to compute collision process
+LOGICAL, SAVE :: GFIRSTCALL = .TRUE. ! control switch for the first call
+!
+INTEGER :: ICIBU
+INTEGER, DIMENSION(:), ALLOCATABLE :: IVEC2_S1,IVEC2_S2 ! Snow indice vector
+INTEGER, DIMENSION(:), ALLOCATABLE :: IVEC2_G ! Graupel indice vector
+INTEGER, PARAMETER :: I_SEED_PARAM = 26032012
+INTEGER, DIMENSION(:), ALLOCATABLE :: I_SEED
+INTEGER :: NI_SEED
+!
+REAL, DIMENSION(:), ALLOCATABLE :: ZVEC1_S, ZVEC1_S1, ZVEC1_S2, & ! Work vectors
+ ZVEC1_S3, ZVEC1_S4, &
+ ZVEC1_S11, ZVEC1_S12, & ! for snow
+ ZVEC1_S21, ZVEC1_S22, &
+ ZVEC1_S31, ZVEC1_S32, &
+ ZVEC1_S41, ZVEC1_S42, &
+ ZVEC2_S1, ZVEC2_S2
+REAL, DIMENSION(:), ALLOCATABLE :: ZVEC1_G, ZVEC1_G1, ZVEC1_G2, & ! Work vectors
+ ZVEC2_G ! for graupel
+REAL, DIMENSION(:), ALLOCATABLE :: ZFRAGMENTS, ZHARVEST
+REAL, DIMENSION(SIZE(PRST)) :: ZINTG_SNOW_1, & ! incomplete gamma function
+ ZINTG_SNOW_2, & ! for snow
+ ZINTG_SNOW_3, &
+ ZINTG_SNOW_4
+REAL, DIMENSION(SIZE(PRST)) :: ZINTG_GRAUPEL_1, & ! incomplete gamma
+ ZINTG_GRAUPEL_2 ! function for graupel
+REAL, DIMENSION(SIZE(PRST)) :: ZNI_CIBU, ZRI_CIBU ! CIBU rates
+REAL, DIMENSION(SIZE(PRST)) :: ZFRAG_CIBU
+REAL :: ZFACT1_XNDEBRIS, ZFACT2_XNDEBRIS
+!
+!-------------------------------------------------------------------------------
+
+GCIBU(:) = LCIBU .AND. (PRST(:)>XRTMIN(5)) .AND. (PRGT(:)>XRTMIN(6)) .AND. LDCOMPUTE(:)
+ICIBU = COUNT( GCIBU(:) )
+!
+P_RI_CIBU(:)=0.
+P_CI_CIBU(:)=0.
+!
+IF (ICIBU > 0) THEN
+!
+! 1.3.0 randomization of XNDEBRIS_CIBU values
+!
+ IF (GFIRSTCALL) THEN
+ CALL RANDOM_SEED(SIZE=NI_SEED) ! get size of seed
+ ALLOCATE(I_SEED(NI_SEED))
+ I_SEED(:) = I_SEED_PARAM !
+ CALL RANDOM_SEED(PUT=I_SEED)
+ GFIRSTCALL = .FALSE.
+ END IF
+!
+ ALLOCATE(ZFRAGMENTS(ICIBU))
+!
+ IF (XNDEBRIS_CIBU >= 0.0) THEN
+ ZFRAGMENTS(:) = XNDEBRIS_CIBU
+ ELSE
+!
+! Mantissa gives the mean value (randomization around 10**MANTISSA)
+! First digit after the comma provides the full range around 10**MANTISSA
+!
+ ALLOCATE(ZHARVEST(ICIBU))
+!
+ ZFACT1_XNDEBRIS = AINT(XNDEBRIS_CIBU)
+ ZFACT2_XNDEBRIS = ABS(ANINT(10.0*(XNDEBRIS_CIBU - ZFACT1_XNDEBRIS)))
+!
+ CALL RANDOM_NUMBER(ZHARVEST(:))
+!
+ ZFRAGMENTS(:) = 10.0**(ZFACT2_XNDEBRIS*ZHARVEST(:) + ZFACT1_XNDEBRIS)
+!
+ DEALLOCATE(ZHARVEST)
+!
+! ZFRAGMENTS is a random variable containing the number of fragments per collision
+! For XNDEBRIS_CIBU=-1.2345 => ZFRAGMENTS(:) = 10.0**(2.0*RANDOM_NUMBER(ZHARVEST(:)) - 1.0)
+! and ZFRAGMENTS=[0.1, 10.0] centered around 1.0
+!
+ END IF
+!
+!
+! 1.3.1 To compute the partial integration of snow gamma function
+!
+! 1.3.1.0 allocations
+!
+ ALLOCATE(ZVEC1_S(ICIBU))
+ ALLOCATE(ZVEC1_S1(ICIBU))
+ ALLOCATE(ZVEC1_S2(ICIBU))
+ ALLOCATE(ZVEC1_S3(ICIBU))
+ ALLOCATE(ZVEC1_S4(ICIBU))
+ ALLOCATE(ZVEC1_S11(ICIBU))
+ ALLOCATE(ZVEC1_S12(ICIBU))
+ ALLOCATE(ZVEC1_S21(ICIBU))
+ ALLOCATE(ZVEC1_S22(ICIBU))
+ ALLOCATE(ZVEC1_S31(ICIBU))
+ ALLOCATE(ZVEC1_S32(ICIBU))
+ ALLOCATE(ZVEC1_S41(ICIBU))
+ ALLOCATE(ZVEC1_S42(ICIBU))
+ ALLOCATE(ZVEC2_S1(ICIBU))
+ ALLOCATE(IVEC2_S1(ICIBU))
+ ALLOCATE(ZVEC2_S2(ICIBU))
+ ALLOCATE(IVEC2_S2(ICIBU))
+!
+!
+! 1.3.1.1 select the PLBDS
+!
+ ZVEC1_S(:) = PACK( PLBDS(:),MASK=GCIBU(:) )
+!
+!
+! 1.3.1.2 find the next lower indice for the PLBDS in the
+! geometrical set of Lbda_s used to tabulate some moments of the
+! incomplete gamma function, for boundary 1 (0.2 mm)
+!
+ ZVEC2_S1(1:ICIBU) = MAX( 1.0001, MIN( FLOAT(NGAMINC)-0.0001,XCIBUINTP_S &
+ * LOG( ZVEC1_S(1:ICIBU) ) + XCIBUINTP1_S ) )
+ IVEC2_S1(1:ICIBU) = INT( ZVEC2_S1(1:ICIBU) )
+ ZVEC2_S1(1:ICIBU) = ZVEC2_S1(1:ICIBU) - FLOAT( IVEC2_S1(1:ICIBU) )
+!
+!
+! 1.3.1.3 find the next lower indice for the PLBDS in the
+! geometrical set of Lbda_s used to tabulate some moments of the
+! incomplete gamma function, for boundary 2 (1 mm)
+!
+ ZVEC2_S2(1:ICIBU) = MAX( 1.0001, MIN( FLOAT(NGAMINC)-0.0001,XCIBUINTP_S &
+ * LOG( ZVEC1_S(1:ICIBU) ) + XCIBUINTP2_S ) )
+ IVEC2_S2(1:ICIBU) = INT( ZVEC2_S2(1:ICIBU) )
+ ZVEC2_S2(1:ICIBU) = ZVEC2_S2(1:ICIBU) - FLOAT( IVEC2_S2(1:ICIBU) )
+!
+!
+! 1.3.1.4 perform the linear interpolation of the
+! normalized "0"-moment of the incomplete gamma function
+!
+! For lower boundary (0.2 mm)
+ ZVEC1_S11(1:ICIBU) = XGAMINC_CIBU_S(1,IVEC2_S1(1:ICIBU)+1) * ZVEC2_S1(1:ICIBU) &
+ - XGAMINC_CIBU_S(1,IVEC2_S1(1:ICIBU)) * (ZVEC2_S1(1:ICIBU)-1.0)
+!
+! For upper boundary (1 mm)
+ ZVEC1_S12(1:ICIBU) = XGAMINC_CIBU_S(1,IVEC2_S2(1:ICIBU)+1) * ZVEC2_S2(1:ICIBU) &
+ - XGAMINC_CIBU_S(1,IVEC2_S2(1:ICIBU)) * (ZVEC2_S2(1:ICIBU)-1.0)
+!
+! Computation of spectrum from 0.2 mm to 1 mm
+ ZVEC1_S1(1:ICIBU) = ZVEC1_S12(1:ICIBU) - ZVEC1_S11(1:ICIBU)
+!
+!
+! 1.3.1.5 perform the linear interpolation of the
+! normalized "XDS"-moment of the incomplete gamma function
+!
+! For lower boundary (0.2 mm)
+ ZVEC1_S21(1:ICIBU) = XGAMINC_CIBU_S(2,IVEC2_S1(1:ICIBU)+1) * ZVEC2_S1(1:ICIBU) &
+ - XGAMINC_CIBU_S(2,IVEC2_S1(1:ICIBU)) * (ZVEC2_S1(1:ICIBU)-1.0)
+!
+! For upper boundary (1 mm)
+ ZVEC1_S22(1:ICIBU) = XGAMINC_CIBU_S(2,IVEC2_S2(1:ICIBU)+1) * ZVEC2_S2(1:ICIBU) &
+ - XGAMINC_CIBU_S(2,IVEC2_S2(1:ICIBU)) * (ZVEC2_S2(1:ICIBU)-1.0)
+!
+! From 0.2 mm to 1 mm we need
+ ZVEC1_S2(1:ICIBU) = XMOMGS_CIBU_1 * (ZVEC1_S22(1:ICIBU) - ZVEC1_S21(1:ICIBU))
+!
+! For lower boundary (0.2 mm)
+ ZVEC1_S31(1:ICIBU) = XGAMINC_CIBU_S(3,IVEC2_S1(1:ICIBU)+1) * ZVEC2_S1(1:ICIBU) &
+ - XGAMINC_CIBU_S(3,IVEC2_S1(1:ICIBU)) * (ZVEC2_S1(1:ICIBU)-1.0)
+!
+! For upper boundary (1 mm)
+ ZVEC1_S32(1:ICIBU) = XGAMINC_CIBU_S(3,IVEC2_S2(1:ICIBU)+1) * ZVEC2_S2(1:ICIBU) &
+ - XGAMINC_CIBU_S(3,IVEC2_S2(1:ICIBU)) * (ZVEC2_S2(1:ICIBU)-1.0)
+!
+! From 0.2 mm to 1 mm we need
+ ZVEC1_S3(1:ICIBU) = XMOMGS_CIBU_2 * (ZVEC1_S32(1:ICIBU) - ZVEC1_S31(1:ICIBU))
+!
+!
+! 1.3.1.6 perform the linear interpolation of the
+! normalized "XBS+XDS"-moment of the incomplete gamma function
+!
+! For lower boundary (0.2 mm)
+ ZVEC1_S41(1:ICIBU) = XGAMINC_CIBU_S(4,IVEC2_S1(1:ICIBU)+1) * ZVEC2_S1(1:ICIBU) &
+ - XGAMINC_CIBU_S(4,IVEC2_S1(1:ICIBU)) * (ZVEC2_S1(1:ICIBU)-1.0)
+!
+! For upper boundary (1 mm)
+ ZVEC1_S42(1:ICIBU) = XGAMINC_CIBU_S(4,IVEC2_S2(1:ICIBU)+1) * ZVEC2_S2(1:ICIBU) &
+ - XGAMINC_CIBU_S(4,IVEC2_S2(1:ICIBU)) * (ZVEC2_S2(1:ICIBU)-1.0)
+!
+! From 0.2 mm to 1 mm we need
+ ZVEC1_S4(1:ICIBU) = XMOMGS_CIBU_3 * (ZVEC1_S42(1:ICIBU) - ZVEC1_S41(1:ICIBU))
+!
+ ZINTG_SNOW_1(:) = UNPACK ( VECTOR=ZVEC1_S1(:),MASK=GCIBU,FIELD=0.0 )
+ ZINTG_SNOW_2(:) = UNPACK ( VECTOR=ZVEC1_S2(:),MASK=GCIBU,FIELD=0.0 )
+ ZINTG_SNOW_3(:) = UNPACK ( VECTOR=ZVEC1_S3(:),MASK=GCIBU,FIELD=0.0 )
+ ZINTG_SNOW_4(:) = UNPACK ( VECTOR=ZVEC1_S4(:),MASK=GCIBU,FIELD=0.0 )
+!
+!
+! 1.3.2 Compute the partial integration of graupel gamma function
+!
+! 1.3.2.0 allocations
+!
+ ALLOCATE(ZVEC1_G(ICIBU))
+ ALLOCATE(ZVEC1_G1(ICIBU))
+ ALLOCATE(ZVEC1_G2(ICIBU))
+ ALLOCATE(ZVEC2_G(ICIBU))
+ ALLOCATE(IVEC2_G(ICIBU))
+!
+!
+! 1.3.2.1 select the PLBDG
+!
+ ZVEC1_G(:) = PACK( PLBDG(:),MASK=GCIBU(:) )
+!
+!
+! 1.3.2.2 find the next lower indice for the PLBDG in the
+! geometrical set of Lbda_g used to tabulate some moments of the
+! incomplete gamma function, for the "2mm" boundary
+!
+ ZVEC2_G(1:ICIBU) = MAX( 1.0001, MIN( FLOAT(NGAMINC)-0.0001,XCIBUINTP_G &
+ * LOG( ZVEC1_G(1:ICIBU) ) + XCIBUINTP1_G ) )
+ IVEC2_G(1:ICIBU) = INT( ZVEC2_G(1:ICIBU) )
+ ZVEC2_G(1:ICIBU) = ZVEC2_G(1:ICIBU) - FLOAT( IVEC2_G(1:ICIBU) )
+!
+!
+! 1.3.2.3 perform the linear interpolation of the
+! normalized "2+XDG"-moment of the incomplete gamma function
+!
+ ZVEC1_G1(1:ICIBU) = XGAMINC_CIBU_G(1,IVEC2_G(1:ICIBU)+1) * ZVEC2_G(1:ICIBU) &
+ - XGAMINC_CIBU_G(1,IVEC2_G(1:ICIBU)) * (ZVEC2_G(1:ICIBU)-1.0)
+!
+! From 2 mm to infinity we need
+ ZVEC1_G1(1:ICIBU) = XMOMGG_CIBU_1 * (1.0 - ZVEC1_G1(1:ICIBU))
+!
+!
+! 1.3.2.4 perform the linear interpolation of the
+! normalized "2.0"-moment of the incomplete gamma function
+!
+ ZVEC1_G2(1:ICIBU) = XGAMINC_CIBU_G(2,IVEC2_G(1:ICIBU)+1) * ZVEC2_G(1:ICIBU) &
+ - XGAMINC_CIBU_G(2,IVEC2_G(1:ICIBU)) * (ZVEC2_G(1:ICIBU)-1.0)
+!
+! From 2 mm to infinity we need
+ ZVEC1_G2(1:ICIBU) = XMOMGG_CIBU_2 * (1.0 - ZVEC1_G2(1:ICIBU))
+!
+!
+ ZINTG_GRAUPEL_1(:) = UNPACK ( VECTOR=ZVEC1_G1(:),MASK=GCIBU,FIELD=0.0 )
+ ZINTG_GRAUPEL_2(:) = UNPACK ( VECTOR=ZVEC1_G2(:),MASK=GCIBU,FIELD=0.0 )
+!
+!
+! 1.3.3 To compute final "CIBU" contributions
+!
+ ZFRAG_CIBU(:) = UNPACK ( VECTOR=ZFRAGMENTS(:),MASK=GCIBU,FIELD=0.0 )
+ ZNI_CIBU(:) = ZFRAG_CIBU(:) * (XFACTOR_CIBU_NI * PRST(:) * PRHODREF(:) / (PRHODREF(:)**(XCEXVT-1.0))) * &
+ (XCG * ZINTG_GRAUPEL_1(:) * ZINTG_SNOW_1(:) * &
+ PLBDS(:)**(XBS) * PLBDG(:)**(XCXG-(XDG+2.0)) &
+ - XCS * ZINTG_GRAUPEL_2(:) * ZINTG_SNOW_2(:) * &
+ PLBDS(:)**(-XDS+XBS) * PLBDG(:)**(XCXG-2.0) * &
+ (1+(XFVELOS/PLBDS(:))**XALPHAS)**(-XNUS-XDS/XALPHAS) )
+
+ P_CI_CIBU(:) = MAX(ZNI_CIBU(:),0.)
+!
+ DEALLOCATE(ZFRAGMENTS)
+!
+! Max value of rs removed by CIBU
+ ZRI_CIBU(:) = (XFACTOR_CIBU_RI * PRST(:) * PRHODREF(:) / (PRHODREF(:)**(XCEXVT+1.0))) * &
+ (XCG * ZINTG_GRAUPEL_1(:) * ZINTG_SNOW_3(:) * &
+ PLBDG(:)**(XCXG-(XDG+2.0)) &
+ - XCS * ZINTG_GRAUPEL_2(:) * ZINTG_SNOW_4(:) * &
+ PLBDS(:)**(-XDS) * PLBDG(:)**(XCXG-2.0) * &
+ (1+(XFVELOS/PLBDS(:))**XALPHAS)**(-XNUS-(XBS+XDS)/XALPHAS) )
+!
+! The value of rs removed by CIBU is determined by the mean mass of pristine ice
+ WHERE( PRIT(:)>XRTMIN(4) .AND. PCIT(:)>XCTMIN(4) )
+ ZRI_CIBU(:) = MIN( ZRI_CIBU(:), ZNI_CIBU(:)*PRIT(:)/PCIT(:) )
+ ELSE WHERE
+ ZRI_CIBU(:) = MIN( ZRI_CIBU(:), MAX( ZNI_CIBU(:)*XMNU0,XRTMIN(4) ) )
+ END WHERE
+!
+ P_RI_CIBU(:) = MAX(ZRI_CIBU(:), 0.)
+!
+ DEALLOCATE(ZVEC1_S)
+ DEALLOCATE(ZVEC1_S1)
+ DEALLOCATE(ZVEC1_S2)
+ DEALLOCATE(ZVEC1_S3)
+ DEALLOCATE(ZVEC1_S4)
+ DEALLOCATE(ZVEC1_S11)
+ DEALLOCATE(ZVEC1_S12)
+ DEALLOCATE(ZVEC1_S21)
+ DEALLOCATE(ZVEC1_S22)
+ DEALLOCATE(ZVEC1_S31)
+ DEALLOCATE(ZVEC1_S32)
+ DEALLOCATE(ZVEC1_S41)
+ DEALLOCATE(ZVEC1_S42)
+ DEALLOCATE(ZVEC2_S1)
+ DEALLOCATE(IVEC2_S1)
+ DEALLOCATE(ZVEC2_S2)
+ DEALLOCATE(IVEC2_S2)
+ DEALLOCATE(ZVEC1_G)
+ DEALLOCATE(ZVEC1_G1)
+ DEALLOCATE(ZVEC1_G2)
+ DEALLOCATE(ZVEC2_G)
+ DEALLOCATE(IVEC2_G)
+END IF
+!
+!-------------------------------------------------------------------------------
+!
+END SUBROUTINE LIMA_COLLISIONAL_ICE_BREAKUP
diff --git a/src/MNH/lima_mixed_fast_processes.f90 b/src/MNH/lima_mixed_fast_processes.f90
index 662eecdd9..24e282c6b 100644
--- a/src/MNH/lima_mixed_fast_processes.f90
+++ b/src/MNH/lima_mixed_fast_processes.f90
@@ -686,11 +686,12 @@ IF (ICIBU > 0) THEN
ALLOCATE(ZFRAG_CIBU(SIZE(PZT)))
!
ZFRAG_CIBU(:) = UNPACK ( VECTOR=ZFRAGMENTS(:),MASK=GCIBU,FIELD=0.0 )
- ZNI_CIBU(:) = ZFRAG_CIBU(:) * (XFACTOR_CIBU_NI / (PRHODREF(:)**(XCEXVT-1.0))) * &
- (XCG * ZINTG_GRAUPEL_1(:) * ZINTG_SNOW_1(:) * &
- PLBDAS(:)**(XCXS) * PLBDAG(:)**(XCXG-(XDG+2.0)) &
- - XCS * ZINTG_GRAUPEL_2(:) * ZINTG_SNOW_2(:) * &
- PLBDAS(:)**(XCXS-XDS) * PLBDAG(:)**(XCXG-2.0) )
+ ZNI_CIBU(:) = ZFRAG_CIBU(:) * (XFACTOR_CIBU_NI * PRST1D(:) * PRHODREF(:) / (PRHODREF(:)**(XCEXVT-1.0))) * &
+ (XCG * ZINTG_GRAUPEL_1(:) * ZINTG_SNOW_1(:) * &
+ PLBDAS(:)**(XBS) * PLBDAG(:)**(XCXG-(XDG+2.0)) &
+ - XCS * ZINTG_GRAUPEL_2(:) * ZINTG_SNOW_2(:) * &
+ PLBDAS(:)**(-XDS+XBS) * PLBDAG(:)**(XCXG-2.0) * &
+ (1+(XFVELOS/PLBDAS(:))**XALPHAS)**(-XNUS-XDS/XALPHAS) )
PCIS1D(:) = PCIS1D(:) + MAX(ZNI_CIBU(:), 0.)
!
DEALLOCATE(ZFRAG_CIBU)
@@ -698,11 +699,12 @@ IF (ICIBU > 0) THEN
!
! Max value of rs removed by CIBU
ALLOCATE(ZRI_CIBU(SIZE(PZT)))
- ZRI_CIBU(:) = (XFACTOR_CIBU_RI / (PRHODREF(:)**(XCEXVT+1.0))) * &
- (XCG * ZINTG_GRAUPEL_1(:) * ZINTG_SNOW_3(:) * &
- PLBDAS(:)**(XCXS-XBS) * PLBDAG(:)**(XCXG-(XDG+2.0)) &
- - XCS * ZINTG_GRAUPEL_2(:) * ZINTG_SNOW_4(:) * &
- PLBDAS(:)**(XCXS-(XBS+XDS)) * PLBDAG(:)**(XCXG-2.0))
+ ZRI_CIBU(:) = (XFACTOR_CIBU_RI * PRST1D(:) * PRHODREF(:) / (PRHODREF(:)**(XCEXVT+1.0))) * &
+ (XCG * ZINTG_GRAUPEL_1(:) * ZINTG_SNOW_3(:) * &
+ PLBDAG(:)**(XCXG-(XDG+2.0)) &
+ - XCS * ZINTG_GRAUPEL_2(:) * ZINTG_SNOW_4(:) * &
+ PLBDAS(:)**(-XDS) * PLBDAG(:)**(XCXG-2.0) * &
+ (1+(XFVELOS/PLBDAS(:))**XALPHAS)**(-XNUS-(XBS+XDS)/XALPHAS) )
!
! The value of rs removed by CIBU is determined by the mean mass of pristine ice
WHERE( PRIT1D(:)>XRTMIN(4) .AND. PCIT1D(:)>XCTMIN(4) )
diff --git a/src/MNH/lima_raindrop_shattering_freezing.f90 b/src/MNH/lima_raindrop_shattering_freezing.f90
new file mode 100644
index 000000000..c64dca2e5
--- /dev/null
+++ b/src/MNH/lima_raindrop_shattering_freezing.f90
@@ -0,0 +1,152 @@
+!MNH_LIC Copyright 2018-2021 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
+!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
+!MNH_LIC for details. version 1.
+!-------------------------------------------------------------------------------
+! #############################################
+ MODULE MODI_LIMA_RAINDROP_SHATTERING_FREEZING
+! #############################################
+!
+INTERFACE
+ SUBROUTINE LIMA_RAINDROP_SHATTERING_FREEZING (LDCOMPUTE, &
+ PRHODREF, &
+ PRRT, PCRT, PRIT, PCIT, &
+ PLBDR, &
+ P_RI_RDSF, P_CI_RDSF )
+!
+LOGICAL, DIMENSION(:),INTENT(IN) :: LDCOMPUTE
+!
+REAL, DIMENSION(:), INTENT(IN) :: PRHODREF
+!
+REAL, DIMENSION(:), INTENT(IN) :: PRRT
+REAL, DIMENSION(:), INTENT(IN) :: PCRT
+REAL, DIMENSION(:), INTENT(IN) :: PRIT
+REAL, DIMENSION(:), INTENT(IN) :: PCIT
+REAL, DIMENSION(:), INTENT(IN) :: PLBDR
+!
+REAL, DIMENSION(:), INTENT(OUT) :: P_RI_RDSF
+REAL, DIMENSION(:), INTENT(OUT) :: P_CI_RDSF
+!
+END SUBROUTINE LIMA_RAINDROP_SHATTERING_FREEZING
+END INTERFACE
+END MODULE MODI_LIMA_RAINDROP_SHATTERING_FREEZING
+!
+! #######################################################################
+ SUBROUTINE LIMA_RAINDROP_SHATTERING_FREEZING (LDCOMPUTE, &
+ PRHODREF, &
+ PRRT, PCRT, PRIT, PCIT, &
+ PLBDR, &
+ P_RI_RDSF, P_CI_RDSF )
+! #######################################################################
+!
+!! PURPOSE
+!! -------
+!! Compute the raindrop shattering when freezing (secondary ice production process)
+!!
+!! AUTHOR
+!! ------
+!! J.-P. Pinty * Laboratoire d'Aerologie*
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 04/2022 duplicate from original for LIMA_SPLIT
+!
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODD_PARAM_LIMA, ONLY : XRTMIN, XCTMIN, LRDSF, XCEXVT
+USE MODD_PARAM_LIMA_COLD, ONLY : XMNU0
+USE MODD_PARAM_LIMA_WARM, ONLY : XDR
+USE MODD_PARAM_LIMA_MIXED, ONLY : NGAMINC, XGAMINC_RDSF_R, &
+ XFACTOR_RDSF_NI, XMOMGR_RDSF, XRDSFINTP1_R, XRDSFINTP_R
+!
+IMPLICIT NONE
+!
+!* 0.1 Declarations of dummy arguments :
+!
+!
+LOGICAL, DIMENSION(:),INTENT(IN) :: LDCOMPUTE
+!
+REAL, DIMENSION(:), INTENT(IN) :: PRHODREF
+!
+REAL, DIMENSION(:), INTENT(IN) :: PRRT
+REAL, DIMENSION(:), INTENT(IN) :: PCRT
+REAL, DIMENSION(:), INTENT(IN) :: PRIT
+REAL, DIMENSION(:), INTENT(IN) :: PCIT
+REAL, DIMENSION(:), INTENT(IN) :: PLBDR
+!
+REAL, DIMENSION(:), INTENT(OUT) :: P_RI_RDSF
+REAL, DIMENSION(:), INTENT(OUT) :: P_CI_RDSF
+!
+!
+!* 0.2 Declarations of local variables :
+!
+LOGICAL, DIMENSION(SIZE(PRRT)) :: GRDSF ! Test where to compute collision process
+INTEGER :: IRDSF
+REAL, DIMENSION(:), ALLOCATABLE :: ZVEC1_R ! Work vectors for rain
+REAL, DIMENSION(:), ALLOCATABLE :: ZVEC1_R1 ! Work vectors for rain
+REAL, DIMENSION(:), ALLOCATABLE :: ZVEC2_R ! Work vectors for rain
+INTEGER, DIMENSION(:), ALLOCATABLE :: IVEC2_R ! Rain indice vector
+REAL, DIMENSION(SIZE(PRRT)) :: ZINTG_RAIN ! incomplete gamma function for rain
+REAL, DIMENSION(SIZE(PRRT)) :: ZNI_RDSF,ZRI_RDSF ! RDSF rates
+!
+!-------------------------------------------------------------------------------
+GRDSF(:) = LRDSF .AND. LDCOMPUTE .AND. (PRIT(:)>XRTMIN(4)) .AND. (PRRT(:)>XRTMIN(3)) &
+ .AND. (PCIT(:)>XCTMIN(4)) .AND. (PCRT(:)>XCTMIN(3))
+IRDSF = COUNT( GRDSF(:) )
+!
+IF (IRDSF > 0) THEN
+!
+ ALLOCATE(ZVEC1_R(IRDSF))
+ ALLOCATE(ZVEC1_R1(IRDSF))
+ ALLOCATE(ZVEC2_R(IRDSF))
+ ALLOCATE(IVEC2_R(IRDSF))
+!
+!* 2.2.1 select the ZLBDAR
+!
+ ZVEC1_R(:) = PACK( PLBDR(:),MASK=GRDSF(:) )
+!
+!* 2.2.2 find the next lower indice for the ZLBDAR in the
+! geometrical set of Lbda_r used to tabulate some moments of the
+! incomplete gamma function, for the lower boundary (0.1 mm)
+!
+ ZVEC2_R(1:IRDSF) = MAX( 1.00001, MIN( FLOAT(NGAMINC)-0.00001,XRDSFINTP_R &
+ * LOG( ZVEC1_R(1:IRDSF) ) + XRDSFINTP1_R ) )
+ IVEC2_R(1:IRDSF) = INT( ZVEC2_R(1:IRDSF) )
+ ZVEC2_R(1:IRDSF) = ZVEC2_R(1:IRDSF) - FLOAT( IVEC2_R(1:IRDSF) )
+!
+!* 2.2.3 perform the linear interpolation of the
+! normalized "2+XDR"-moment of the incomplete gamma function
+!
+ ZVEC1_R1(1:IRDSF) = XGAMINC_RDSF_R(IVEC2_R(1:IRDSF)+1) * ZVEC2_R(1:IRDSF) &
+ - XGAMINC_RDSF_R(IVEC2_R(1:IRDSF)) * (ZVEC2_R(1:IRDSF) - 1.0)
+!
+! From 0.1 mm to infinity we need
+ ZVEC1_R1(1:IRDSF) = XMOMGR_RDSF * (1.0 - ZVEC1_R1(1:IRDSF))
+!
+ ZINTG_RAIN(:) = UNPACK ( VECTOR=ZVEC1_R1(:),MASK=GRDSF,FIELD=0.0 )
+!
+!* 2.2.4 To compute final "RDSF" contributions
+!
+ ZNI_RDSF(:) = (XFACTOR_RDSF_NI / (PRHODREF(:)**(XCEXVT-1.0))) * ( &
+ PCIT(:) * PCRT(:) * ZINTG_RAIN(:) * PLBDR(:)**(-(XDR+6.0)) )
+!
+ P_CI_RDSF(:) = MAX(ZNI_RDSF(:),0.)
+!
+! The value of rg removed by RDSF is determined by the mean mass of pristine ice
+ ZRI_RDSF(:) = MAX( ZNI_RDSF(:)*XMNU0,XRTMIN(5) )
+!
+ P_RI_RDSF(:) = ZRI_RDSF(:)
+!
+ DEALLOCATE(ZVEC1_R)
+ DEALLOCATE(ZVEC1_R1)
+ DEALLOCATE(ZVEC2_R)
+ DEALLOCATE(IVEC2_R)
+ !
+ENDIF
+!
+!-------------------------------------------------------------------------------
+!
+END SUBROUTINE LIMA_RAINDROP_SHATTERING_FREEZING
diff --git a/src/MNH/lima_tendencies.f90 b/src/MNH/lima_tendencies.f90
index 0fce8cb11..06e9252c4 100644
--- a/src/MNH/lima_tendencies.f90
+++ b/src/MNH/lima_tendencies.f90
@@ -29,6 +29,8 @@ MODULE MODI_LIMA_TENDENCIES
P_TH_ACC, P_RR_ACC, P_CR_ACC, P_RS_ACC, P_RG_ACC, &
P_RS_CMEL, &
P_TH_CFRZ, P_RR_CFRZ, P_CR_CFRZ, P_RI_CFRZ, P_CI_CFRZ, &
+ P_RI_CIBU, P_CI_CIBU, &
+ P_RI_RDSF, P_CI_RDSF, &
P_TH_WETG, P_RC_WETG, P_CC_WETG, P_RR_WETG, P_CR_WETG, &
P_RI_WETG, P_CI_WETG, P_RS_WETG, P_RG_WETG, P_RH_WETG, &
P_TH_DRYG, P_RC_DRYG, P_CC_DRYG, P_RR_DRYG, P_CR_DRYG, &
@@ -128,6 +130,12 @@ REAL, DIMENSION(:), INTENT(INOUT) :: P_CR_CFRZ
REAL, DIMENSION(:), INTENT(INOUT) :: P_RI_CFRZ
REAL, DIMENSION(:), INTENT(INOUT) :: P_CI_CFRZ ! rain freezing (CFRZ) : rr, Nr, ri, Ni, rg=-rr-ri, th
!
+REAL, DIMENSION(:), INTENT(INOUT) :: P_RI_CIBU
+REAL, DIMENSION(:), INTENT(INOUT) :: P_CI_CIBU ! collisional ice break-up (CIBU) : ri, Ni, rs=-ri
+!
+REAL, DIMENSION(:), INTENT(INOUT) :: P_RI_RDSF
+REAL, DIMENSION(:), INTENT(INOUT) :: P_CI_RDSF ! rain drops freezing shattering (RDSF) : ri, Ni, rg=-ri
+!
REAL, DIMENSION(:), INTENT(INOUT) :: P_TH_WETG
REAL, DIMENSION(:), INTENT(INOUT) :: P_RC_WETG
REAL, DIMENSION(:), INTENT(INOUT) :: P_CC_WETG
@@ -203,6 +211,8 @@ SUBROUTINE LIMA_TENDENCIES (PTSTEP, LDCOMPUTE,
P_TH_ACC, P_RR_ACC, P_CR_ACC, P_RS_ACC, P_RG_ACC, &
P_RS_CMEL, &
P_TH_CFRZ, P_RR_CFRZ, P_CR_CFRZ, P_RI_CFRZ, P_CI_CFRZ, &
+ P_RI_CIBU, P_CI_CIBU, &
+ P_RI_RDSF, P_CI_RDSF, &
P_TH_WETG, P_RC_WETG, P_CC_WETG, P_RR_WETG, P_CR_WETG, &
P_RI_WETG, P_CI_WETG, P_RS_WETG, P_RG_WETG, P_RH_WETG, &
P_TH_DRYG, P_RC_DRYG, P_CC_DRYG, P_RR_DRYG, P_CR_DRYG, &
@@ -242,7 +252,7 @@ SUBROUTINE LIMA_TENDENCIES (PTSTEP, LDCOMPUTE,
!
USE MODD_CST, ONLY : XP00, XRD, XRV, XMD, XMV, XCPD, XCPV, XCL, XCI, XLVTT, XLSTT, XTT, &
XALPW, XBETAW, XGAMW, XALPI, XBETAI, XGAMI
-USE MODD_PARAM_LIMA, ONLY : XRTMIN, XCTMIN, XNUS, &
+USE MODD_PARAM_LIMA, ONLY : XRTMIN, XCTMIN, XNUS, LCIBU, LRDSF, &
LCOLD, LNUCL, LSNOW, LHAIL, LWARM, LACTI, LRAIN, LKHKO, LSNOW_T, NMOM_I
USE MODD_PARAM_LIMA_WARM, ONLY : XLBC, XLBEXC, XLBR, XLBEXR
USE MODD_PARAM_LIMA_MIXED, ONLY : XLBG, XLBEXG, XLBH, XLBEXH, XLBDAG_MAX
@@ -263,6 +273,8 @@ USE MODI_LIMA_DROPLETS_RIMING_SNOW
USE MODI_LIMA_RAIN_ACCR_SNOW
USE MODI_LIMA_CONVERSION_MELTING_SNOW
USE MODI_LIMA_RAIN_FREEZING
+USE MODI_LIMA_COLLISIONAL_ICE_BREAKUP
+USE MODI_LIMA_RAINDROP_SHATTERING_FREEZING
USE MODI_LIMA_GRAUPEL
!
USE MODI_LIMA_BERGERON
@@ -355,6 +367,12 @@ REAL, DIMENSION(:), INTENT(INOUT) :: P_CR_CFRZ
REAL, DIMENSION(:), INTENT(INOUT) :: P_RI_CFRZ
REAL, DIMENSION(:), INTENT(INOUT) :: P_CI_CFRZ ! rain freezing (CFRZ) : rr, Nr, ri, Ni, rg=-rr-ri, th
!
+REAL, DIMENSION(:), INTENT(INOUT) :: P_RI_CIBU
+REAL, DIMENSION(:), INTENT(INOUT) :: P_CI_CIBU ! collisional ice break-up (CIBU) : ri, Ni, rs=-ri
+!
+REAL, DIMENSION(:), INTENT(INOUT) :: P_RI_RDSF
+REAL, DIMENSION(:), INTENT(INOUT) :: P_CI_RDSF ! rain drops freezing shattering (RDSF) : ri, Ni, rg=-ri
+!
REAL, DIMENSION(:), INTENT(INOUT) :: P_TH_WETG
REAL, DIMENSION(:), INTENT(INOUT) :: P_RC_WETG
REAL, DIMENSION(:), INTENT(INOUT) :: P_CC_WETG
@@ -784,6 +802,44 @@ IF (LWARM .AND. LRAIN .AND. LCOLD ) THEN
PA_CI(:) = PA_CI(:) + P_CI_CFRZ(:)
PA_RG(:) = PA_RG(:) - P_RR_CFRZ(:) - P_RI_CFRZ(:)
+END IF
+!
+IF (LWARM .AND. LCOLD .AND. LSNOW .AND. LCIBU) THEN
+ !
+ ! Conversion melting of snow should account for collected droplets and drops where T>0C, but does not !
+ ! Some thermodynamical computations inside, to externalize ?
+ !
+ CALL LIMA_COLLISIONAL_ICE_BREAKUP (LDCOMPUTE, & ! depends on PF (IF for fragments size)
+ PRHODREF, &
+ PRIT/ZIF1D, PRST/ZPF1D, PRGT/ZPF1D, PCIT/ZIF1D, &
+ ZLBDS, ZLBDG, &
+ P_RI_CIBU, P_CI_CIBU )
+ P_RI_CIBU(:) = P_RI_CIBU(:) * ZPF1D(:)
+ P_CI_CIBU(:) = P_CI_CIBU(:) * ZPF1D(:)
+ !
+ PA_RI(:) = PA_RI(:) + P_RI_CIBU(:)
+ PA_CI(:) = PA_CI(:) + P_CI_CIBU(:)
+ PA_RS(:) = PA_RS(:) - P_RI_CIBU(:)
+
+END IF
+!
+IF (LWARM .AND. LRAIN .AND. LCOLD .AND. LSNOW .AND. LRDSF) THEN
+ !
+ ! Conversion melting of snow should account for collected droplets and drops where T>0C, but does not !
+ ! Some thermodynamical computations inside, to externalize ?
+ !
+ CALL LIMA_RAINDROP_SHATTERING_FREEZING (LDCOMPUTE, & ! depends on PF, IF
+ PRHODREF, &
+ PRRT/ZPF1D, PCRT/ZPF1D, PRIT/ZIF1D, PCIT/ZIF1D, &
+ ZLBDR, &
+ P_RI_RDSF, P_CI_RDSF )
+ P_RI_RDSF(:) = P_RI_RDSF(:) * ZIF1D(:)
+ P_CI_RDSF(:) = P_CI_RDSF(:) * ZIF1D(:)
+ !
+ PA_RI(:) = PA_RI(:) + P_RI_RDSF(:)
+ PA_CI(:) = PA_CI(:) + P_CI_RDSF(:)
+ PA_RG(:) = PA_RG(:) - P_RI_RDSF(:)
+
END IF
!
IF (LWARM .AND. LCOLD) THEN
diff --git a/src/MNH/modd_param_lima_cold.f90 b/src/MNH/modd_param_lima_cold.f90
index dde912b42..cb23cb131 100644
--- a/src/MNH/modd_param_lima_cold.f90
+++ b/src/MNH/modd_param_lima_cold.f90
@@ -127,40 +127,6 @@ REAL,SAVE :: XCONCI_MAX ! Limitation of the pristine
! ice concentration (init and grid-nesting)
REAL,SAVE :: XFREFFI ! Factor to compute the cloud ice effective radius
!
-!
-! Constants for ice-ice collision : CIBU
-!
-REAL, SAVE :: XDCSLIM_CIBU_MIN, & ! aggregates min diam. : 0.2 mm
- XDCSLIM_CIBU_MAX, & ! aggregates max diam. : 1.0 mm
- XDCGLIM_CIBU_MIN, & ! graupel min diam. : 2 mm
- XGAMINC_BOUND_CIBU_SMIN, & ! Min val. of Lbda_s*dlim
- XGAMINC_BOUND_CIBU_SMAX, & ! Max val. of Lbda_s*dlim
- XGAMINC_BOUND_CIBU_GMIN, & ! Min val. of Lbda_g*dlim
- XGAMINC_BOUND_CIBU_GMAX, & ! Max val. of Lbda_g*dlim
- XCIBUINTP_S,XCIBUINTP1_S, & !
- XCIBUINTP2_S, & !
- XCIBUINTP_G,XCIBUINTP1_G, & !
- XFACTOR_CIBU_NI,XFACTOR_CIBU_RI, & ! Factor for final CIBU Eq.
- XMOMGG_CIBU_1,XMOMGG_CIBU_2, & ! Moment computation
- XMOMGS_CIBU_1,XMOMGS_CIBU_2, &
- XMOMGS_CIBU_3
-!
-REAL, DIMENSION(:,:), SAVE, ALLOCATABLE &
- :: XGAMINC_CIBU_S, & ! Tab.incomplete Gamma function
- XGAMINC_CIBU_G ! Tab.incomplete Gamma function
-!
-! Constants for raindrop shattering : RDSF
-!
-REAL, SAVE :: XDCRLIM_RDSF_MIN, & ! Raindrops min diam. : 0.2 mm
- XGAMINC_BOUND_RDSF_RMIN, & ! Min val. of Lbda_r*dlim
- XGAMINC_BOUND_RDSF_RMAX, & ! Max val. of Lbda_r*dlim
- XRDSFINTP_R,XRDSFINTP1_R, & !
- XFACTOR_RDSF_NI, & ! Factor for final RDSF Eq.
- XMOMGR_RDSF
-!
-REAL, DIMENSION(:), SAVE, ALLOCATABLE &
- :: XGAMINC_RDSF_R ! Tab.incomplete Gamma function
-!
!-------------------------------------------------------------------------------
!
END MODULE MODD_PARAM_LIMA_COLD
diff --git a/src/MNH/modd_param_lima_mixed.f90 b/src/MNH/modd_param_lima_mixed.f90
index 57ea4d1a5..b92e8603a 100644
--- a/src/MNH/modd_param_lima_mixed.f90
+++ b/src/MNH/modd_param_lima_mixed.f90
@@ -14,6 +14,7 @@
!! MODIFICATIONS
!! -------------
!! Original ??/??/13
+!! C. Barthe 14/03/2022 add CIBU and RDSF
! J. Wurtz 03/2022: new snow characteristics
!!
!-------------------------------------------------------------------------------
@@ -49,7 +50,44 @@ REAL,SAVE :: XALPHAH,XNUH,XLBEXH,XLBH ! Hail distribution parameters
!
!-------------------------------------------------------------------------------
!
-!* 2. MICROPHYSICAL FACTORS - Graupel
+!* 2. MICROPHYSICAL FACTORS - CIBU and RDSF
+! -------------------------------------
+!
+! Constants for ice-ice collision : CIBU
+!
+REAL, SAVE :: XDCSLIM_CIBU_MIN, & ! aggregates min diam. : 0.2 mm
+ XDCSLIM_CIBU_MAX, & ! aggregates max diam. : 1.0 mm
+ XDCGLIM_CIBU_MIN, & ! graupel min diam. : 2 mm
+ XGAMINC_BOUND_CIBU_SMIN, & ! Min val. of Lbda_s*dlim
+ XGAMINC_BOUND_CIBU_SMAX, & ! Max val. of Lbda_s*dlim
+ XGAMINC_BOUND_CIBU_GMIN, & ! Min val. of Lbda_g*dlim
+ XGAMINC_BOUND_CIBU_GMAX, & ! Max val. of Lbda_g*dlim
+ XCIBUINTP_S,XCIBUINTP1_S, & !
+ XCIBUINTP2_S, & !
+ XCIBUINTP_G,XCIBUINTP1_G, & !
+ XFACTOR_CIBU_NI,XFACTOR_CIBU_RI, & ! Factor for final CIBU Eq.
+ XMOMGG_CIBU_1,XMOMGG_CIBU_2, & ! Moment computation
+ XMOMGS_CIBU_1,XMOMGS_CIBU_2, &
+ XMOMGS_CIBU_3
+!
+REAL, DIMENSION(:,:), SAVE, ALLOCATABLE &
+ :: XGAMINC_CIBU_S, & ! Tab.incomplete Gamma function
+ XGAMINC_CIBU_G ! Tab.incomplete Gamma function
+!
+! Constants for raindrop shattering : RDSF
+!
+REAL, SAVE :: XDCRLIM_RDSF_MIN, & ! Raindrops min diam. : 0.2 mm
+ XGAMINC_BOUND_RDSF_RMIN, & ! Min val. of Lbda_r*dlim
+ XGAMINC_BOUND_RDSF_RMAX, & ! Max val. of Lbda_r*dlim
+ XRDSFINTP_R,XRDSFINTP1_R, & !
+ XFACTOR_RDSF_NI, & ! Factor for final RDSF Eq.
+ XMOMGR_RDSF
+!
+REAL, DIMENSION(:), SAVE, ALLOCATABLE &
+ :: XGAMINC_RDSF_R ! Tab.incomplete Gamma function
+!
+!
+!* 3. MICROPHYSICAL FACTORS - Graupel
! -------------------------------
!
REAL,SAVE :: XFSEDG, XEXSEDG ! Sedimentation fluxes of Graupel
@@ -134,7 +172,7 @@ REAL,DIMENSION(:,:), SAVE, ALLOCATABLE &
!
!-------------------------------------------------------------------------------
!
-!* 2. MICROPHYSICAL FACTORS - Hail
+!* 4. MICROPHYSICAL FACTORS - Hail
! ----------------------------
!
REAL,SAVE :: XFSEDH,XEXSEDH ! Constants for sedimentation
--
GitLab