diff --git a/docs/TODO b/docs/TODO
index ba3cf158895eb8e06b43759cdb5aa99cc9f1abf0..2d640ac7e0a8bd13f58b6e260afc50f7e10c3330 100644
--- a/docs/TODO
+++ b/docs/TODO
@@ -47,13 +47,12 @@ Pb identifiés à corriger plus tard:
l'option n'est pas remontée en namelist)
- th_r_from_thl_rt appelée partout, il faudrait limiter à OTEST
- sedimentation momentum non branchée
- - si possible, modifier ice4_sedimentation_split* dans le même esprit que stat
+ - si possible, modifier ice4_sedimentation_split* dans le même esprit que stat, supprimer le faux packing
- La taille du buffer utilisé pour th_r_from_thl_rt doit être mise en module et
utilisée pour déclarer le buffer dans les routines appelantes et dans th_r_from_thl_rt
- rain_ice:
- séparer l'avance temporelle du découpage en sous-blocs en créant une couche driver supplémentaire. Cette couche pourrait
avoir différentes implémentations (filtre LLMICRO seul, filtre LLMICRO + découpage en sous-blocs, filtre LLMICR
- - mettre le code des interpolations linéaires et bi-linéaires dans des routines avec deux implémentations: avec et sans packing
- shuman, turb, shuman_mf, apl_arome
- peut-on remettre à jour la partie sous le sol des variables pronostiques avant d'entrer dans EDKFi dans Méso-NH? Si oui, on devrait pouvoir
récrire les shuman_mf pour qu'ils produisent les mêmes résultats qu'aujourd'hui mais sans utiliser les valeurs sous le sol ou dans "l'espace infini"
diff --git a/src/common/micro/mode_ice4_correct_negativities.F90 b/src/common/micro/mode_ice4_correct_negativities.F90
new file mode 100644
index 0000000000000000000000000000000000000000..12f8b08aa4eed8fc9565d2fc7462d48cf92eabff
--- /dev/null
+++ b/src/common/micro/mode_ice4_correct_negativities.F90
@@ -0,0 +1,118 @@
+!MNH_LIC Copyright 1995-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 MODE_ICE4_CORRECT_NEGATIVITIES
+IMPLICIT NONE
+CONTAINS
+SUBROUTINE ICE4_CORRECT_NEGATIVITIES(D, ICED, KRR, PRV, PRC, PRR, &
+ &PRI, PRS, PRG, &
+ &PTH, PLVFACT, PLSFACT, PRH)
+!
+USE PARKIND1, ONLY : JPRB
+USE YOMHOOK , ONLY : LHOOK, DR_HOOK
+USE MODD_DIMPHYEX, ONLY: DIMPHYEX_t
+USE MODD_RAIN_ICE_DESCR, ONLY: RAIN_ICE_DESCR_t
+!
+IMPLICIT NONE
+!
+TYPE(DIMPHYEX_t), INTENT(IN) :: D
+TYPE(RAIN_ICE_DESCR_t), INTENT(IN) :: ICED
+INTEGER, INTENT(IN) :: KRR
+REAL, DIMENSION(D%NIJT, D%NKT), INTENT(INOUT) :: PRV, PRC, PRR, PRI, PRS, PRG, PTH
+REAL, DIMENSION(D%NIJT, D%NKT), INTENT(IN) :: PLVFACT, PLSFACT
+REAL, DIMENSION(D%NIJT, D%NKT), OPTIONAL, INTENT(INOUT) :: PRH
+!
+REAL :: ZW
+INTEGER :: JIJ, JK, IKTB, IKTE, IIJB, IIJE
+
+REAL(KIND=JPRB) :: ZHOOK_HANDLE
+!
+IF (LHOOK) CALL DR_HOOK('ICE4_CORRECT_NEGATIVITIES', 0, ZHOOK_HANDLE)
+!
+IKTB=D%NKTB
+IKTE=D%NKTE
+IIJB=D%NIJB
+IIJE=D%NIJE
+!
+!We correct negativities with conservation
+DO JK = IKTB, IKTE
+ DO JIJ = IIJB, IIJE
+ ! 1) deal with negative values for mixing ratio, except for vapor
+ ZW =PRC(JIJ,JK)-MAX(PRC(JIJ,JK), 0.)
+ PRV(JIJ,JK)=PRV(JIJ,JK)+ZW
+ PTH(JIJ,JK)=PTH(JIJ,JK)-ZW*PLVFACT(JIJ,JK)
+ PRC(JIJ,JK)=PRC(JIJ,JK)-ZW
+
+ ZW =PRR(JIJ,JK)-MAX(PRR(JIJ,JK), 0.)
+ PRV(JIJ,JK)=PRV(JIJ,JK)+ZW
+ PTH(JIJ,JK)=PTH(JIJ,JK)-ZW*PLVFACT(JIJ,JK)
+ PRR(JIJ,JK)=PRR(JIJ,JK)-ZW
+
+ ZW =PRI(JIJ,JK)-MAX(PRI(JIJ,JK), 0.)
+ PRV(JIJ,JK)=PRV(JIJ,JK)+ZW
+ PTH(JIJ,JK)=PTH(JIJ,JK)-ZW*PLSFACT(JIJ,JK)
+ PRI(JIJ,JK)=PRI(JIJ,JK)-ZW
+
+ ZW =PRS(JIJ,JK)-MAX(PRS(JIJ,JK), 0.)
+ PRV(JIJ,JK)=PRV(JIJ,JK)+ZW
+ PTH(JIJ,JK)=PTH(JIJ,JK)-ZW*PLSFACT(JIJ,JK)
+ PRS(JIJ,JK)=PRS(JIJ,JK)-ZW
+
+ ZW =PRG(JIJ,JK)-MAX(PRG(JIJ,JK), 0.)
+ PRV(JIJ,JK)=PRV(JIJ,JK)+ZW
+ PTH(JIJ,JK)=PTH(JIJ,JK)-ZW*PLSFACT(JIJ,JK)
+ PRG(JIJ,JK)=PRG(JIJ,JK)-ZW
+
+ IF(KRR==7) THEN
+ ZW =PRH(JIJ,JK)-MAX(PRH(JIJ,JK), 0.)
+ PRV(JIJ,JK)=PRV(JIJ,JK)+ZW
+ PTH(JIJ,JK)=PTH(JIJ,JK)-ZW*PLSFACT(JIJ,JK)
+ PRH(JIJ,JK)=PRH(JIJ,JK)-ZW
+ ENDIF
+
+ ! 2) deal with negative vapor mixing ratio
+
+ ! for rc and ri, we keep ice fraction constant
+ ZW=MIN(1., MAX(ICED%XRTMIN(1)-PRV(JIJ,JK), 0.) / &
+ &MAX(PRC(JIJ,JK)+PRI(JIJ,JK), 1.E-20)) ! Proportion of rc+ri to convert into rv
+ PTH(JIJ,JK)=PTH(JIJ,JK)-ZW* &
+ &(PRC(JIJ,JK)*PLVFACT(JIJ,JK)+PRI(JIJ,JK)*PLSFACT(JIJ,JK))
+ PRV(JIJ,JK)=PRV(JIJ,JK)+ZW*(PRC(JIJ,JK)+PRI(JIJ,JK))
+ PRC(JIJ,JK)=(1.-ZW)*PRC(JIJ,JK)
+ PRI(JIJ,JK)=(1.-ZW)*PRI(JIJ,JK)
+
+ ZW=MIN(MAX(PRR(JIJ,JK), 0.), &
+ &MAX(ICED%XRTMIN(1)-PRV(JIJ,JK), 0.)) ! Quantity of rr to convert into rv
+ PRV(JIJ,JK)=PRV(JIJ,JK)+ZW
+ PRR(JIJ,JK)=PRR(JIJ,JK)-ZW
+ PTH(JIJ,JK)=PTH(JIJ,JK)-ZW*PLVFACT(JIJ,JK)
+
+ ZW=MIN(MAX(PRS(JIJ,JK), 0.), &
+ &MAX(ICED%XRTMIN(1)-PRV(JIJ,JK), 0.)) ! Quantity of rs to convert into rv
+ PRV(JIJ,JK)=PRV(JIJ,JK)+ZW
+ PRS(JIJ,JK)=PRS(JIJ,JK)-ZW
+ PTH(JIJ,JK)=PTH(JIJ,JK)-ZW*PLSFACT(JIJ,JK)
+
+ ZW=MIN(MAX(PRG(JIJ,JK), 0.), &
+ &MAX(ICED%XRTMIN(1)-PRV(JIJ,JK), 0.)) ! Quantity of rg to convert into rv
+ PRV(JIJ,JK)=PRV(JIJ,JK)+ZW
+ PRG(JIJ,JK)=PRG(JIJ,JK)-ZW
+ PTH(JIJ,JK)=PTH(JIJ,JK)-ZW*PLSFACT(JIJ,JK)
+
+ IF(KRR==7) THEN
+ ZW=MIN(MAX(PRH(JIJ,JK), 0.), &
+ &MAX(ICED%XRTMIN(1)-PRV(JIJ,JK), 0.)) ! Quantity of rh to convert into rv
+ PRV(JIJ,JK)=PRV(JIJ,JK)+ZW
+ PRH(JIJ,JK)=PRH(JIJ,JK)-ZW
+ PTH(JIJ,JK)=PTH(JIJ,JK)-ZW*PLSFACT(JIJ,JK)
+ ENDIF
+ ENDDO
+ENDDO
+!
+IF (LHOOK) CALL DR_HOOK('ICE4_CORRECT_NEGATIVITIES', 1, ZHOOK_HANDLE)
+!
+END SUBROUTINE ICE4_CORRECT_NEGATIVITIES
+!
+END MODULE MODE_ICE4_CORRECT_NEGATIVITIES
diff --git a/src/common/micro/mode_ice4_sedimentation.F90 b/src/common/micro/mode_ice4_sedimentation.F90
new file mode 100644
index 0000000000000000000000000000000000000000..310edd5123b778896a0244f3cf990b0585c7731a
--- /dev/null
+++ b/src/common/micro/mode_ice4_sedimentation.F90
@@ -0,0 +1,236 @@
+!MNH_LIC Copyright 1994-2019 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
+!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
+!MNH_LIC for details. version 1.
+!-----------------------------------------------------------------
+MODULE MODE_ICE4_SEDIMENTATION
+IMPLICIT NONE
+CONTAINS
+SUBROUTINE ICE4_SEDIMENTATION(D, CST, ICEP, ICED, PARAMI, BUCONF, &
+ &PTSTEP, KRR, PDZZ, &
+ &PLVFACT, PLSFACT, PRHODREF, PPABST, PTHT, PT, PRHODJ, &
+ &PTHS, PRVS, PRCS, PRCT, PRRS, PRRT, PRIS, PRIT, PRSS, PRST, PRGS, PRGT,&
+ &PINPRC, PINPRR, PINPRS, PINPRG, &
+ &TBUDGETS, KBUDGETS, &
+ &PSEA, PTOWN, &
+ &PINPRH, PRHT, PRHS, PFPR)
+!!
+!!** PURPOSE
+!! -------
+!! Computes the sedimentation
+!!
+!! AUTHOR
+!! ------
+!! S. Riette code extracted from rain_ice
+!!
+!! MODIFICATIONS
+!! -------------
+!!
+!
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE PARKIND1, ONLY : JPRB
+USE YOMHOOK , ONLY : LHOOK, DR_HOOK
+USE MODD_DIMPHYEX, ONLY: DIMPHYEX_t
+USE MODD_BUDGET, ONLY: TBUDGETDATA, TBUDGETCONF_t, NBUDGET_TH, NBUDGET_RV, NBUDGET_RC, &
+ NBUDGET_RI, NBUDGET_RR, NBUDGET_RS, NBUDGET_RG, NBUDGET_RH
+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
+!
+USE MODE_MSG, ONLY: PRINT_MSG, NVERB_FATAL
+USE MODE_BUDGET, ONLY: BUDGET_STORE_INIT_PHY, BUDGET_STORE_END_PHY
+!
+USE MODE_ICE4_SEDIMENTATION_STAT, ONLY: ICE4_SEDIMENTATION_STAT
+USE MODE_ICE4_SEDIMENTATION_SPLIT, ONLY: ICE4_SEDIMENTATION_SPLIT
+USE MODE_ICE4_SEDIMENTATION_SPLIT_MOMENTUM, ONLY: ICE4_SEDIMENTATION_SPLIT_MOMENTUM
+USE MODE_ICE4_CORRECT_NEGATIVITIES, ONLY: ICE4_CORRECT_NEGATIVITIES
+!
+!
+IMPLICIT NONE
+!
+!* 0.1 Declarations of dummy arguments :
+!
+TYPE(DIMPHYEX_t), INTENT(IN) :: D !array dimensions
+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
+TYPE(TBUDGETCONF_t), INTENT(IN) :: BUCONF
+REAL, INTENT(IN) :: PTSTEP ! Double Time step (single if cold start)
+INTEGER, INTENT(IN) :: KRR ! Number of moist variable
+REAL, DIMENSION(D%NIJT,D%NKT), INTENT(IN) :: PDZZ ! Layer thikness (m)
+REAL, DIMENSION(D%NIJT,D%NKT), INTENT(IN) :: PLVFACT
+REAL, DIMENSION(D%NIJT,D%NKT), INTENT(IN) :: PLSFACT
+REAL, DIMENSION(D%NIJT,D%NKT), INTENT(IN) :: PRHODREF! Reference density
+REAL, DIMENSION(D%NIJT,D%NKT), INTENT(IN) :: PPABST ! absolute pressure at t
+REAL, DIMENSION(D%NIJT,D%NKT), INTENT(IN) :: PTHT ! Theta at time t
+REAL, DIMENSION(D%NIJT,D%NKT), INTENT(IN) :: PT ! Temperature at time t
+REAL, DIMENSION(D%NIJT,D%NKT), INTENT(IN) :: PRHODJ ! Dry density * Jacobian
+REAL, DIMENSION(D%NIJT,D%NKT), INTENT(INOUT) :: PTHS
+REAL, DIMENSION(D%NIJT,D%NKT), INTENT(INOUT) :: PRVS
+REAL, DIMENSION(D%NIJT,D%NKT), INTENT(INOUT) :: PRCS ! Cloud water m.r. source
+REAL, DIMENSION(D%NIJT,D%NKT), INTENT(IN) :: PRCT ! Cloud water m.r. at t
+REAL, DIMENSION(D%NIJT,D%NKT), INTENT(INOUT) :: PRRS ! Rain water m.r. source
+REAL, DIMENSION(D%NIJT,D%NKT), INTENT(IN) :: PRRT ! Rain water m.r. at t
+REAL, DIMENSION(D%NIJT,D%NKT), INTENT(INOUT) :: PRIS ! Pristine ice m.r. source
+REAL, DIMENSION(D%NIJT,D%NKT), INTENT(IN) :: PRIT ! Pristine ice m.r. at t
+REAL, DIMENSION(D%NIJT,D%NKT), INTENT(INOUT) :: PRSS ! Snow/aggregate m.r. source
+REAL, DIMENSION(D%NIJT,D%NKT), INTENT(IN) :: PRST ! Snow/aggregate m.r. at t
+REAL, DIMENSION(D%NIJT,D%NKT), INTENT(INOUT) :: PRGS ! Graupel m.r. source
+REAL, DIMENSION(D%NIJT,D%NKT), INTENT(IN) :: PRGT ! Graupel/hail m.r. at t
+REAL, DIMENSION(D%NIJT), INTENT(OUT) :: PINPRC ! Cloud instant precip
+REAL, DIMENSION(D%NIJT), INTENT(OUT) :: PINPRR ! Rain instant precip
+REAL, DIMENSION(D%NIJT), INTENT(OUT) :: PINPRS ! Snow instant precip
+REAL, DIMENSION(D%NIJT), INTENT(OUT) :: PINPRG ! Graupel instant precip
+TYPE(TBUDGETDATA), DIMENSION(KBUDGETS), INTENT(INOUT) :: TBUDGETS
+INTEGER, INTENT(IN) :: KBUDGETS
+REAL, DIMENSION(D%NIJT), OPTIONAL, INTENT(IN) :: PSEA ! Sea Mask
+REAL, DIMENSION(D%NIJT), OPTIONAL, INTENT(IN) :: PTOWN ! Fraction that is town
+REAL, DIMENSION(D%NIJT), OPTIONAL, INTENT(OUT) :: PINPRH ! Hail instant precip
+REAL, DIMENSION(D%NIJT,D%NKT), OPTIONAL, INTENT(IN) :: PRHT ! Hail m.r. at t
+REAL, DIMENSION(D%NIJT,D%NKT), OPTIONAL, INTENT(INOUT) :: PRHS ! Hail m.r. source
+REAL, DIMENSION(D%NIJT,D%NKT,KRR), OPTIONAL, INTENT(OUT) :: PFPR ! upper-air precipitation fluxes
+!
+!* 0.2 declaration of local variables
+!
+REAL, DIMENSION(D%NIJT,D%NKT) :: ZRCT, ZRRT, ZRIT, ZRST, ZRGT, ZRHT
+REAL, DIMENSION(D%NIJT) :: ZINPRI
+INTEGER :: JK, JIJ, IKTB, IKTE, IIJB, IIJE
+REAL(KIND=JPRB) :: ZHOOK_HANDLE
+!
+IF (LHOOK) CALL DR_HOOK('ICE4_SEDIMENTATION', 0, ZHOOK_HANDLE)
+IKTB=D%NKTB
+IKTE=D%NKTE
+IIJB=D%NIJB
+IIJE=D%NIJE
+!
+!
+!-------------------------------------------------------------------------------
+!
+!* 2. COMPUTE THE SEDIMENTATION (RS) SOURCE
+! -------------------------------------
+!
+!
+!
+IF (BUCONF%LBUDGET_RC .AND. PARAMI%LSEDIC) CALL BUDGET_STORE_INIT_PHY(D, TBUDGETS(NBUDGET_RC), 'SEDI', PRCS(:, :) * PRHODJ(:, :))
+IF (BUCONF%LBUDGET_RR) CALL BUDGET_STORE_INIT_PHY(D, TBUDGETS(NBUDGET_RR), 'SEDI', PRRS(:, :) * PRHODJ(:, :))
+IF (BUCONF%LBUDGET_RI) CALL BUDGET_STORE_INIT_PHY(D, TBUDGETS(NBUDGET_RI), 'SEDI', PRIS(:, :) * PRHODJ(:, :))
+IF (BUCONF%LBUDGET_RS) CALL BUDGET_STORE_INIT_PHY(D, TBUDGETS(NBUDGET_RS), 'SEDI', PRSS(:, :) * PRHODJ(:, :))
+IF (BUCONF%LBUDGET_RG) CALL BUDGET_STORE_INIT_PHY(D, TBUDGETS(NBUDGET_RG), 'SEDI', PRGS(:, :) * PRHODJ(:, :))
+IF (BUCONF%LBUDGET_RH .AND. KRR==7) CALL BUDGET_STORE_INIT_PHY(D, TBUDGETS(NBUDGET_RH), 'SEDI', PRHS(:, :) * PRHODJ(:, :))
+
+IF(PARAMI%CSEDIM=='STAT') THEN
+ IF (KRR==7) THEN
+ DO JK = IKTB,IKTE
+ DO JIJ = IIJB,IIJE
+ ZRCT(JIJ,JK)=PRCS(JIJ,JK)*PTSTEP
+ ZRRT(JIJ,JK)=PRRS(JIJ,JK)*PTSTEP
+ ZRIT(JIJ,JK)=PRIS(JIJ,JK)*PTSTEP
+ ZRST(JIJ,JK)=PRSS(JIJ,JK)*PTSTEP
+ ZRGT(JIJ,JK)=PRGS(JIJ,JK)*PTSTEP
+ ZRHT(JIJ,JK)=PRHS(JIJ,JK)*PTSTEP
+ ENDDO
+ ENDDO
+ CALL ICE4_SEDIMENTATION_STAT(D, CST, ICEP, ICED, PARAMI, &
+ &PTSTEP, KRR, PDZZ, &
+ &PRHODREF, PPABST, PTHT, PT, PRHODJ, &
+ &PRCS, ZRCT, PRRS, ZRRT, PRIS, ZRIT,&
+ &PRSS, ZRST, PRGS, ZRGT,&
+ &PINPRC, PINPRR, ZINPRI, PINPRS, PINPRG, &
+ &PSEA=PSEA, PTOWN=PTOWN, &
+ &PINPRH=PINPRH, PRHT=ZRHT, PRHS=PRHS, PFPR=PFPR)
+ ELSE
+ DO JK = IKTB,IKTE
+ DO JIJ = IIJB,IIJE
+ ZRCT(JIJ,JK)=PRCS(JIJ,JK)*PTSTEP
+ ZRRT(JIJ,JK)=PRRS(JIJ,JK)*PTSTEP
+ ZRIT(JIJ,JK)=PRIS(JIJ,JK)*PTSTEP
+ ZRST(JIJ,JK)=PRSS(JIJ,JK)*PTSTEP
+ ZRGT(JIJ,JK)=PRGS(JIJ,JK)*PTSTEP
+ ENDDO
+ ENDDO
+ CALL ICE4_SEDIMENTATION_STAT(D, CST, ICEP, ICED, PARAMI, &
+ &PTSTEP, KRR, PDZZ, &
+ &PRHODREF, PPABST, PTHT, PT, PRHODJ, &
+ &PRCS, ZRCT, PRRS, ZRRT, PRIS, ZRIT,&
+ &PRSS, ZRST, PRGS, ZRGT,&
+ &PINPRC, PINPRR, ZINPRI, PINPRS, PINPRG, &
+ &PSEA=PSEA, PTOWN=PTOWN, &
+ &PFPR=PFPR)
+ ENDIF
+ PINPRS(IIJB:IIJE) = PINPRS(IIJB:IIJE) + ZINPRI(IIJB:IIJE)
+ !No negativity correction here as we apply sedimentation on PR.S*PTSTEP variables
+ELSEIF(PARAMI%CSEDIM=='SPLI') THEN
+ IF(KRR==7) THEN
+ CALL ICE4_SEDIMENTATION_SPLIT(D, CST, ICEP, ICED, PARAMI, &
+ &PTSTEP, KRR, PDZZ, &
+ &PRHODREF, PPABST, PTHT, PT, PRHODJ, &
+ &PRCS, PRCT, PRRS, PRRT, PRIS, PRIT, PRSS, PRST, PRGS, PRGT,&
+ &PINPRC, PINPRR, ZINPRI, PINPRS, PINPRG, &
+ &PSEA=PSEA, PTOWN=PTOWN, &
+ &PINPRH=PINPRH, PRHT=PRHT, PRHS=PRHS, PFPR=PFPR)
+ ELSE
+ CALL ICE4_SEDIMENTATION_SPLIT(D, CST, ICEP, ICED, PARAMI, &
+ &PTSTEP, KRR, PDZZ, &
+ &PRHODREF, PPABST, PTHT, PT, PRHODJ, &
+ &PRCS, PRCT, PRRS, PRRT, PRIS, PRIT, PRSS, PRST, PRGS, PRGT,&
+ &PINPRC, PINPRR, ZINPRI, PINPRS, PINPRG, &
+ &PSEA=PSEA, PTOWN=PTOWN, &
+ &PFPR=PFPR)
+ ENDIF
+ PINPRS(IIJB:IIJE) = PINPRS(IIJB:IIJE) + ZINPRI(IIJB:IIJE)
+ !We correct negativities with conservation
+ !SPLI algorith uses a time-splitting. Inside the loop a temporary m.r. is used.
+ ! It is initialized with the m.r. at T and is modified by two tendencies:
+ ! sedimentation tendency and an external tendency which represents all other
+ ! processes (mainly advection and microphysical processes). If both tendencies
+ ! are negative, sedimentation can remove a species at a given sub-timestep. From
+ ! this point sedimentation stops for the remaining sub-timesteps but the other tendency
+ ! will be still active and will lead to negative values.
+ ! We could prevent the algorithm to not consume too much a species, instead we apply
+ ! a correction here.
+ CALL ICE4_CORRECT_NEGATIVITIES(D, ICED, KRR, PRVS, PRCS, PRRS, &
+ &PRIS, PRSS, PRGS, &
+ &PTHS, PLVFACT, PLSFACT, PRHS)
+ELSEIF(PARAMI%CSEDIM=='NONE') THEN
+ELSE
+ CALL PRINT_MSG(NVERB_FATAL, 'GEN', 'RAIN_ICE', 'no sedimentation scheme for PARAMI%CSEDIM='//PARAMI%CSEDIM)
+END IF
+
+
+
+
+
+
+!!! ajouter momentum
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+!
+!
+IF (BUCONF%LBUDGET_RC .AND. PARAMI%LSEDIC) CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_RC), 'SEDI', PRCS(:, :) * PRHODJ(:, :))
+IF (BUCONF%LBUDGET_RR) CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_RR), 'SEDI', PRRS(:, :) * PRHODJ(:, :))
+IF (BUCONF%LBUDGET_RI) CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_RI), 'SEDI', PRIS(:, :) * PRHODJ(:, :))
+IF (BUCONF%LBUDGET_RS) CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_RS), 'SEDI', PRSS(:, :) * PRHODJ(:, :))
+IF (BUCONF%LBUDGET_RG) CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_RG), 'SEDI', PRGS(:, :) * PRHODJ(:, :))
+IF (BUCONF%LBUDGET_RH .AND. KRR==7) CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_RH), 'SEDI', PRHS(:, :) * PRHODJ(:, :))
+!
+IF (LHOOK) CALL DR_HOOK('ICE4_SEDIMENTATION', 1, ZHOOK_HANDLE)
+!
+END SUBROUTINE ICE4_SEDIMENTATION
+END MODULE MODE_ICE4_SEDIMENTATION
diff --git a/src/common/micro/mode_ice4_sedimentation_split.F90 b/src/common/micro/mode_ice4_sedimentation_split.F90
index e7b90b3d4ae6dd7f2c3b4ba22908543ab8e152e1..8d77611e3d17fc0f9819df34a2a10742c46197ee 100644
--- a/src/common/micro/mode_ice4_sedimentation_split.F90
+++ b/src/common/micro/mode_ice4_sedimentation_split.F90
@@ -56,46 +56,46 @@ TYPE(RAIN_ICE_DESCR_t), INTENT(IN) :: ICED
TYPE(PARAM_ICE_t), INTENT(IN) :: PARAMI
REAL, INTENT(IN) :: PTSTEP ! Double Time step (single if cold start)
INTEGER, INTENT(IN) :: KRR ! Number of moist variable
-REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PDZZ ! Layer thikness (m)
-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
-REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(INOUT) :: PRRS ! Rain water m.r. source
-REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PRRT ! Rain water m.r. at t
-REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(INOUT) :: PRIS ! Pristine ice m.r. source
-REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PRIT ! Pristine ice m.r. at t
-REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(INOUT) :: PRSS ! Snow/aggregate m.r. source
-REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PRST ! Snow/aggregate m.r. at t
-REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(INOUT) :: PRGS ! Graupel m.r. source
-REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PRGT ! Graupel/hail m.r. at t
-REAL, DIMENSION(D%NIT,D%NJT), INTENT(OUT) :: PINPRC ! Cloud instant precip
-REAL, DIMENSION(D%NIT,D%NJT), INTENT(OUT) :: PINPRR ! Rain instant precip
-REAL, DIMENSION(D%NIT,D%NJT), INTENT(OUT) :: PINPRI ! Pristine ice instant precip
-REAL, DIMENSION(D%NIT,D%NJT), INTENT(OUT) :: PINPRS ! Snow instant precip
-REAL, DIMENSION(D%NIT,D%NJT), INTENT(OUT) :: PINPRG ! Graupel instant precip
-REAL, DIMENSION(D%NIT,D%NJT), OPTIONAL, INTENT(IN) :: PSEA ! Sea Mask
-REAL, DIMENSION(D%NIT,D%NJT), OPTIONAL, INTENT(IN) :: PTOWN ! Fraction that is town
-REAL, DIMENSION(D%NIT,D%NJT), OPTIONAL, INTENT(OUT) :: PINPRH ! Hail instant precip
-REAL, DIMENSION(D%NIT,D%NJT,D%NKT), OPTIONAL, INTENT(IN) :: PRHT ! Hail m.r. at t
-REAL, DIMENSION(D%NIT,D%NJT,D%NKT), OPTIONAL, INTENT(INOUT) :: PRHS ! Hail m.r. source
-REAL, DIMENSION(D%NIT,D%NJT,D%NKT,KRR), OPTIONAL, INTENT(OUT) :: PFPR ! upper-air precipitation fluxes
+REAL, DIMENSION(D%NIJT,D%NKT), INTENT(IN) :: PDZZ ! Layer thikness (m)
+REAL, DIMENSION(D%NIJT,D%NKT), INTENT(IN) :: PRHODREF! Reference density
+REAL, DIMENSION(D%NIJT,D%NKT), INTENT(IN) :: PPABST ! absolute pressure at t
+REAL, DIMENSION(D%NIJT,D%NKT), INTENT(IN) :: PTHT ! Theta at time t
+REAL, DIMENSION(D%NIJT,D%NKT), INTENT(IN) :: PT ! Temperature at time t
+REAL, DIMENSION(D%NIJT,D%NKT), INTENT(IN) :: PRHODJ ! Dry density * Jacobian
+REAL, DIMENSION(D%NIJT,D%NKT), INTENT(INOUT) :: PRCS ! Cloud water m.r. source
+REAL, DIMENSION(D%NIJT,D%NKT), INTENT(IN) :: PRCT ! Cloud water m.r. at t
+REAL, DIMENSION(D%NIJT,D%NKT), INTENT(INOUT) :: PRRS ! Rain water m.r. source
+REAL, DIMENSION(D%NIJT,D%NKT), INTENT(IN) :: PRRT ! Rain water m.r. at t
+REAL, DIMENSION(D%NIJT,D%NKT), INTENT(INOUT) :: PRIS ! Pristine ice m.r. source
+REAL, DIMENSION(D%NIJT,D%NKT), INTENT(IN) :: PRIT ! Pristine ice m.r. at t
+REAL, DIMENSION(D%NIJT,D%NKT), INTENT(INOUT) :: PRSS ! Snow/aggregate m.r. source
+REAL, DIMENSION(D%NIJT,D%NKT), INTENT(IN) :: PRST ! Snow/aggregate m.r. at t
+REAL, DIMENSION(D%NIJT,D%NKT), INTENT(INOUT) :: PRGS ! Graupel m.r. source
+REAL, DIMENSION(D%NIJT,D%NKT), INTENT(IN) :: PRGT ! Graupel/hail m.r. at t
+REAL, DIMENSION(D%NIJT), INTENT(OUT) :: PINPRC ! Cloud instant precip
+REAL, DIMENSION(D%NIJT), INTENT(OUT) :: PINPRR ! Rain instant precip
+REAL, DIMENSION(D%NIJT), INTENT(OUT) :: PINPRI ! Pristine ice instant precip
+REAL, DIMENSION(D%NIJT), INTENT(OUT) :: PINPRS ! Snow instant precip
+REAL, DIMENSION(D%NIJT), INTENT(OUT) :: PINPRG ! Graupel instant precip
+REAL, DIMENSION(D%NIJT), OPTIONAL, INTENT(IN) :: PSEA ! Sea Mask
+REAL, DIMENSION(D%NIJT), OPTIONAL, INTENT(IN) :: PTOWN ! Fraction that is town
+REAL, DIMENSION(D%NIJT), OPTIONAL, INTENT(OUT) :: PINPRH ! Hail instant precip
+REAL, DIMENSION(D%NIJT,D%NKT), OPTIONAL, INTENT(IN) :: PRHT ! Hail m.r. at t
+REAL, DIMENSION(D%NIJT,D%NKT), OPTIONAL, INTENT(INOUT) :: PRHS ! Hail m.r. source
+REAL, DIMENSION(D%NIJT,D%NKT,KRR), OPTIONAL, INTENT(OUT) :: PFPR ! upper-air precipitation fluxes
!
!* 0.2 declaration of local variables
!
!
-INTEGER :: JI, JJ, JK
-INTEGER :: IKTB, IKTE, IKB, IKL, IIE, IIB, IJB, IJE
+INTEGER :: JIJ, JK
+INTEGER :: IKTB, IKTE, IKB, IKL, IIJE, IIJB
INTEGER :: IRR !Workaround of PGI bug with OpenACC (at least up to 18.10 version)
LOGICAL :: GSEDIC !Workaround of PGI bug with OpenACC (at least up to 18.10 version)
LOGICAL :: GPRESENT_PFPR, GPRESENT_PSEA
REAL :: ZINVTSTEP
-REAL, DIMENSION(D%NIT, D%NJT) :: ZCONC_TMP ! Weighted concentration
-REAL, DIMENSION(D%NIT,D%NJT,D%NKTB:D%NKTE) :: ZW ! work array
-REAL, DIMENSION(D%NIT, D%NJT, D%NKT) :: ZCONC3D, & ! droplet condensation
+REAL, DIMENSION(D%NIJT) :: ZCONC_TMP ! Weighted concentration
+REAL, DIMENSION(D%NIJT,D%NKTB:D%NKTE) :: ZW ! work array
+REAL, DIMENSION(D%NIJT, D%NKT) :: ZCONC3D, & ! droplet condensation
& ZRAY, & ! Cloud Mean radius
& ZLBC, & ! XLBC weighted by sea fraction
& ZFSEDC, &
@@ -118,10 +118,8 @@ IRR = KRR
!
IKTB=D%NKTB
IKTE=D%NKTE
-IIB=D%NIB
-IIE=D%NIE
-IJB=D%NJB
-IJE=D%NJE
+IIJB=D%NIJB
+IIJE=D%NIJE
!
IF (PRESENT(PFPR)) THEN
GPRESENT_PFPR = .TRUE.
@@ -139,41 +137,41 @@ END IF
!
ZINVTSTEP=1./PTSTEP
IF (GPRESENT_PFPR) THEN
- PFPR(:,:,:,:) = 0.
+ PFPR(:,:,:) = 0.
END IF
!
!* 1. Parameters for cloud sedimentation
!
IF (GSEDIC) THEN
- ZRAY(:,:,:) = 0.
- ZLBC(:,:,:) = ICED%XLBC(1)
- ZFSEDC(:,:,:) = ICEP%XFSEDC(1)
- ZCONC3D(:,:,:)= ICED%XCONC_LAND
- ZCONC_TMP(:,:)= ICED%XCONC_LAND
+ ZRAY(:,:) = 0.
+ ZLBC(:,:) = ICED%XLBC(1)
+ ZFSEDC(:,:) = ICEP%XFSEDC(1)
+ ZCONC3D(:,:)= ICED%XCONC_LAND
+ ZCONC_TMP(:)= ICED%XCONC_LAND
IF (GPRESENT_PSEA) THEN
- DO JJ = IJB, IJE
- DO JI = IIB, IIE
- ZCONC_TMP(JI,JJ)=PSEA(JI,JJ)*ICED%XCONC_SEA+(1.-PSEA(JI,JJ))*ICED%XCONC_LAND
- ENDDO
+ DO JIJ = IIJB, IIJE
+ ZCONC_TMP(JIJ)=PSEA(JIJ)*ICED%XCONC_SEA+(1.-PSEA(JIJ))*ICED%XCONC_LAND
ENDDO
DO JK=IKTB, IKTE
- DO JJ = IJB, IJE
- DO JI = IIB, IIE
- ZLBC(JI,JJ,JK) = PSEA(JI,JJ)*ICED%XLBC(2)+(1.-PSEA(JI,JJ))*ICED%XLBC(1)
- ZFSEDC(JI,JJ,JK) = (PSEA(JI,JJ)*ICEP%XFSEDC(2)+(1.-PSEA(JI,JJ))*ICEP%XFSEDC(1))
- ZFSEDC(JI,JJ,JK) = MAX(MIN(ICEP%XFSEDC(1),ICEP%XFSEDC(2)),ZFSEDC(JI,JJ,JK))
- ZCONC3D(JI,JJ,JK)= (1.-PTOWN(JI,JJ))*ZCONC_TMP(JI,JJ)+PTOWN(JI,JJ)*ICED%XCONC_URBAN
- ZRAY(JI,JJ,JK) = 0.5*((1.-PSEA(JI,JJ))*GAMMA(ICED%XNUC+1.0/ICED%XALPHAC)/(GAMMA(ICED%XNUC)) + &
- & PSEA(JI,JJ)*GAMMA(ICED%XNUC2+1.0/ICED%XALPHAC2)/(GAMMA(ICED%XNUC2)))
- ENDDO
+ DO JIJ = IIJB, IIJE
+ ZLBC(JIJ,JK) = PSEA(JIJ)*ICED%XLBC(2)+(1.-PSEA(JIJ))*ICED%XLBC(1)
+ ZFSEDC(JIJ,JK) = (PSEA(JIJ)*ICEP%XFSEDC(2)+(1.-PSEA(JIJ))*ICEP%XFSEDC(1))
+ ZFSEDC(JIJ,JK) = MAX(MIN(ICEP%XFSEDC(1),ICEP%XFSEDC(2)),ZFSEDC(JIJ,JK))
+ ZCONC3D(JIJ,JK)= (1.-PTOWN(JIJ))*ZCONC_TMP(JIJ)+PTOWN(JIJ)*ICED%XCONC_URBAN
+ ZRAY(JIJ,JK) = 0.5*((1.-PSEA(JIJ))*GAMMA(ICED%XNUC+1.0/ICED%XALPHAC)/(GAMMA(ICED%XNUC)) + &
+ & PSEA(JIJ)*GAMMA(ICED%XNUC2+1.0/ICED%XALPHAC2)/(GAMMA(ICED%XNUC2)))
ENDDO
END DO
ELSE
- ZCONC3D(:,:,:) = ICED%XCONC_LAND
- ZRAY(:,:,:) = 0.5*(GAMMA(ICED%XNUC+1.0/ICED%XALPHAC)/(GAMMA(ICED%XNUC)))
+ ZCONC3D(:,:) = ICED%XCONC_LAND
+ ZRAY(:,:) = 0.5*(GAMMA(ICED%XNUC+1.0/ICED%XALPHAC)/(GAMMA(ICED%XNUC)))
END IF
- ZRAY(:,:,:) = MAX(1.,ZRAY(:,:,:))
- ZLBC(:,:,:) = MAX(MIN(ICED%XLBC(1),ICED%XLBC(2)),ZLBC(:,:,:))
+ DO JK=IKTB, IKTE
+ DO JIJ = IIJB, IIJE
+ ZRAY(JIJ,JK) = MAX(1.,ZRAY(JIJ,JK))
+ ZLBC(JIJ,JK) = MAX(MIN(ICED%XLBC(1),ICED%XLBC(2)),ZLBC(JIJ,JK))
+ ENDDO
+ ENDDO
ENDIF
!
!* 2. compute the fluxes
@@ -183,32 +181,30 @@ ENDIF
! For optimization we consider each variable separately
!
DO JK=IKTB, IKTE
- DO JJ = IJB, IJE
- DO JI = IIB, IIE
- ! External tendecies
- IF (GSEDIC) THEN
- ZPRCS(JI,JJ,JK) = PRCS(JI,JJ,JK)-PRCT(JI,JJ,JK)*ZINVTSTEP
- ENDIF
- ZPRRS(JI,JJ,JK) = PRRS(JI,JJ,JK)-PRRT(JI,JJ,JK)*ZINVTSTEP
- ZPRIS(JI,JJ,JK) = PRIS(JI,JJ,JK)-PRIT(JI,JJ,JK)*ZINVTSTEP
- ZPRSS(JI,JJ,JK) = PRSS(JI,JJ,JK)-PRST(JI,JJ,JK)*ZINVTSTEP
- ZPRGS(JI,JJ,JK) = PRGS(JI,JJ,JK)-PRGT(JI,JJ,JK)*ZINVTSTEP
- IF ( IRR == 7 ) THEN
- ZPRHS(JI,JJ,JK) = PRHS(JI,JJ,JK)-PRHT(JI,JJ,JK)*ZINVTSTEP
- END IF
- !
- ! mr values inside the time-splitting loop
- ZRCT(JI,JJ,JK) = PRCT(JI,JJ,JK)
- ZRRT(JI,JJ,JK) = PRRT(JI,JJ,JK)
- ZRIT(JI,JJ,JK) = PRIT(JI,JJ,JK)
- ZRST(JI,JJ,JK) = PRST(JI,JJ,JK)
- ZRGT(JI,JJ,JK) = PRGT(JI,JJ,JK)
- IF (IRR==7) THEN
- ZRHT(JI,JJ,JK) = PRHT(JI,JJ,JK)
- END IF
- !
- ZW(JI,JJ,JK) =1./(PRHODREF(JI,JJ,JK)* PDZZ(JI,JJ,JK))
- ENDDO
+ DO JIJ = IIJB, IIJE
+ ! External tendecies
+ IF (GSEDIC) THEN
+ ZPRCS(JIJ,JK) = PRCS(JIJ,JK)-PRCT(JIJ,JK)*ZINVTSTEP
+ ENDIF
+ ZPRRS(JIJ,JK) = PRRS(JIJ,JK)-PRRT(JIJ,JK)*ZINVTSTEP
+ ZPRIS(JIJ,JK) = PRIS(JIJ,JK)-PRIT(JIJ,JK)*ZINVTSTEP
+ ZPRSS(JIJ,JK) = PRSS(JIJ,JK)-PRST(JIJ,JK)*ZINVTSTEP
+ ZPRGS(JIJ,JK) = PRGS(JIJ,JK)-PRGT(JIJ,JK)*ZINVTSTEP
+ IF ( IRR == 7 ) THEN
+ ZPRHS(JIJ,JK) = PRHS(JIJ,JK)-PRHT(JIJ,JK)*ZINVTSTEP
+ END IF
+ !
+ ! mr values inside the time-splitting loop
+ ZRCT(JIJ,JK) = PRCT(JIJ,JK)
+ ZRRT(JIJ,JK) = PRRT(JIJ,JK)
+ ZRIT(JIJ,JK) = PRIT(JIJ,JK)
+ ZRST(JIJ,JK) = PRST(JIJ,JK)
+ ZRGT(JIJ,JK) = PRGT(JIJ,JK)
+ IF (IRR==7) THEN
+ ZRHT(JIJ,JK) = PRHT(JIJ,JK)
+ END IF
+ !
+ ZW(JIJ,JK) =1./(PRHODREF(JIJ,JK)* PDZZ(JIJ,JK))
ENDDO
ENDDO
!
@@ -296,38 +292,38 @@ 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, 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, DIMENSION(D%NIJT,D%NKT), INTENT(IN) :: PRHODREF ! Reference density
+REAL, DIMENSION(D%NIJT,D%NKTB:D%NKTE), INTENT(IN) :: POORHODZ ! One Over (Rhodref times delta Z)
+REAL, DIMENSION(D%NIJT,D%NKT), INTENT(IN) :: PDZZ ! layer thikness (m)
+REAL, DIMENSION(D%NIJT,D%NKT), INTENT(IN) :: PPABST
+REAL, DIMENSION(D%NIJT,D%NKT), INTENT(IN) :: PTHT
+REAL, DIMENSION(D%NIJT,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
-REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(INOUT) :: PRXS !Tendency of the specy KSPE
-REAL, DIMENSION(D%NIT,D%NJT), INTENT(OUT) :: PINPRX ! instant precip
-REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PPRXS ! external tendencie
-REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN), OPTIONAL :: PRAY, PLBC, PFSEDC, PCONC3D
-REAL, DIMENSION(D%NIT,D%NJT,D%NKT,KRR), INTENT(INOUT), OPTIONAL :: PFPR ! upper-air precipitation fluxes
+REAL, DIMENSION(D%NIJT,D%NKT), INTENT(INOUT) :: PRXT ! mr of specy X
+REAL, DIMENSION(D%NIJT,D%NKT), INTENT(INOUT) :: PRXS !Tendency of the specy KSPE
+REAL, DIMENSION(D%NIJT), INTENT(OUT) :: PINPRX ! instant precip
+REAL, DIMENSION(D%NIJT,D%NKT), INTENT(IN) :: PPRXS ! external tendencie
+REAL, DIMENSION(D%NIJT,D%NKT), INTENT(IN), OPTIONAL :: PRAY, PLBC, PFSEDC, PCONC3D
+REAL, DIMENSION(D%NIJT,D%NKT,KRR), INTENT(INOUT), OPTIONAL :: PFPR ! upper-air precipitation fluxes
!
!* 0.2 declaration of local variables
!
CHARACTER(LEN=10) :: YSPE ! String for error message
INTEGER :: IDX, ISEDIM
-INTEGER :: JI, JJ, JK, JL
-INTEGER, DIMENSION(D%NIT*D%NJT*D%NKT) :: I1,I2,I3 ! Used to replace the COUNT
+INTEGER :: JIJ, JK, JL
+INTEGER, DIMENSION(D%NIJT*D%NKT) :: I1,I2 ! 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
+REAL, DIMENSION(D%NIJT) :: ZMAX_TSTEP ! Maximum CFL in column
REAL, DIMENSION(SIZE(ICED%XRTMIN)) :: ZRSMIN
-REAL, DIMENSION(D%NIT, D%NJT) :: ZREMAINT ! Remaining time until the timestep end
-REAL, DIMENSION(D%NIT, D%NJT, 0:D%NKT+1) :: ZWSED ! Sedimentation fluxes
-INTEGER :: IKTB, IKTE, IKB, IKL, IIE, IIB, IJB, IJE
+REAL, DIMENSION(D%NIJT) :: ZREMAINT ! Remaining time until the timestep end
+REAL, DIMENSION(D%NIJT, 0:D%NKT+1) :: ZWSED ! Sedimentation fluxes
+INTEGER :: IKTB, IKTE, IKB, IKL, IIJE, IIJB
REAL(KIND=JPRB) :: ZHOOK_HANDLE
IF (LHOOK) CALL DR_HOOK('ICE4_SEDIMENTATION_SPLIT:INTERNAL_SEDIM_SPLIT', 0, ZHOOK_HANDLE)
!
@@ -335,10 +331,8 @@ IKTB=D%NKTB
IKTE=D%NKTE
IKB=D%NKB
IKL=D%NKL
-IIB=D%NIB
-IIE=D%NIE
-IJB=D%NJB
-IJE=D%NJE
+IIJB=D%NIJB
+IIJE=D%NIJE
!-------------------------------------------------------------------------------
IF (KSPE<2 .OR. KSPE>7) CALL PRINT_MSG(NVERB_FATAL,'GEN','INTERNAL_SEDIM_SPLIT','invalid species (KSPE variable)')
!
@@ -348,27 +342,24 @@ ELSE
GPRESENT_PFPR = .FALSE.
END IF
!
-PINPRX(:,:) = 0.
+PINPRX(:) = 0.
ZINVTSTEP=1./PTSTEP
ZRSMIN(:) = ICED%XRTMIN(:) * ZINVTSTEP
-ZREMAINT(:,:) = 0.
-ZREMAINT(IIB:IIE,IJB:IJE) = PTSTEP
+ZREMAINT(:) = 0.
+ZREMAINT(IIJB:IIJE) = PTSTEP
!
DO WHILE (ANY(ZREMAINT>0.))
ISEDIM = 0
DO JK = IKTB,IKTE
- DO JJ = IJB,IJE
- DO JI = IIB,IIE
- IF( (PRXT (JI,JJ,JK)>ICED%XRTMIN(KSPE) .OR. &
- PPRXS(JI,JJ,JK)>ZRSMIN(KSPE)) .AND. &
- ZREMAINT(JI,JJ)>0. ) THEN
- ISEDIM = ISEDIM + 1
- IDX = ISEDIM
- I1(IDX) = JI
- I2(IDX) = JJ
- I3(IDX) = JK
- END IF
- END DO
+ DO JIJ = IIJB,IIJE
+ IF( (PRXT (JIJ,JK)>ICED%XRTMIN(KSPE) .OR. &
+ PPRXS(JIJ,JK)>ZRSMIN(KSPE)) .AND. &
+ ZREMAINT(JIJ)>0. ) THEN
+ ISEDIM = ISEDIM + 1
+ IDX = ISEDIM
+ I1(IDX) = JIJ
+ I2(IDX) = JK
+ END IF
END DO
END DO
!
@@ -381,57 +372,54 @@ DO WHILE (ANY(ZREMAINT>0.))
!
IF(KSPE==2) THEN
!******* for cloud
- ZWSED(:,:,:) = 0.
+ ZWSED(:,:) = 0.
DO JL=1, ISEDIM
- JI=I1(JL)
- JJ=I2(JL)
- JK=I3(JL)
- IF(PRXT(JI,JJ,JK)>ICED%XRTMIN(KSPE)) THEN
- ZZWLBDC = PLBC(JI,JJ,JK) * PCONC3D(JI,JJ,JK) / &
- &(PRHODREF(JI,JJ,JK) * PRXT(JI,JJ,JK))
+ JIJ=I1(JL)
+ JK=I2(JL)
+ IF(PRXT(JIJ,JK)>ICED%XRTMIN(KSPE)) THEN
+ ZZWLBDC = PLBC(JIJ,JK) * PCONC3D(JIJ,JK) / &
+ &(PRHODREF(JIJ,JK) * PRXT(JIJ,JK))
ZZWLBDC = ZZWLBDC**ICED%XLBEXC
- ZRAY = PRAY(JI,JJ,JK) / ZZWLBDC
- ZZT = PTHT(JI,JJ,JK) * (PPABST(JI,JJ,JK)/CST%XP00)**(CST%XRD/CST%XCPD)
- ZZWLBDA = 6.6E-8*(101325./PPABST(JI,JJ,JK))*(ZZT/293.15)
+ ZRAY = PRAY(JIJ,JK) / ZZWLBDC
+ ZZT = PTHT(JIJ,JK) * (PPABST(JIJ,JK)/CST%XP00)**(CST%XRD/CST%XCPD)
+ ZZWLBDA = 6.6E-8*(101325./PPABST(JIJ,JK))*(ZZT/293.15)
ZZCC = ICED%XCC*(1.+1.26*ZZWLBDA/ZRAY)
- ZWSED(JI, JJ, JK) = PRHODREF(JI,JJ,JK)**(-ICED%XCEXVT +1 ) * &
- &ZZWLBDC**(-ICED%XDC)*ZZCC*PFSEDC(JI,JJ,JK) * PRXT(JI,JJ,JK)
+ ZWSED(JIJ, JK) = PRHODREF(JIJ,JK)**(-ICED%XCEXVT +1 ) * &
+ &ZZWLBDC**(-ICED%XDC)*ZZCC*PFSEDC(JIJ,JK) * PRXT(JIJ,JK)
ENDIF
ENDDO
ELSEIF(KSPE==4) THEN
! ******* for pristine ice
- ZWSED(:,:,:) = 0.
+ ZWSED(:,:) = 0.
DO JL=1, ISEDIM
- JI=I1(JL)
- JJ=I2(JL)
- JK=I3(JL)
- IF(PRXT(JI, JJ, JK) .GT. MAX(ICED%XRTMIN(4), 1.0E-7)) THEN
- ZWSED(JI, JJ, JK) = ICEP%XFSEDI * PRXT(JI, JJ, JK) * &
- & PRHODREF(JI,JJ,JK)**(1.-ICED%XCEXVT) * & ! McF&H
+ JIJ=I1(JL)
+ JK=I2(JL)
+ IF(PRXT(JIJ, JK) .GT. MAX(ICED%XRTMIN(4), 1.0E-7)) THEN
+ ZWSED(JIJ, JK) = ICEP%XFSEDI * PRXT(JIJ, JK) * &
+ & PRHODREF(JIJ,JK)**(1.-ICED%XCEXVT) * & ! McF&H
& MAX( 0.05E6,-0.15319E6-0.021454E6* &
- & ALOG(PRHODREF(JI,JJ,JK)*PRXT(JI,JJ,JK)) )**ICEP%XEXCSEDI
+ & ALOG(PRHODREF(JIJ,JK)*PRXT(JIJ,JK)) )**ICEP%XEXCSEDI
ENDIF
ENDDO
#if defined(REPRO48) || defined(REPRO55)
#else
ELSEIF(KSPE==5) THEN
! ******* for snow
- ZWSED(:,:,:) = 0.
+ 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
+ JIJ=I1(JL)
+ JK=I2(JL)
+ IF(PRXT(JIJ,JK)> ICED%XRTMIN(KSPE)) THEN
+ IF (PARAMI%LSNOW_T .AND. PT(JIJ,JK)>263.15) THEN
+ ZLBDA = MAX(MIN(ICED%XLBDAS_MAX, 10**(14.554-0.0423*PT(JIJ,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
+ ZLBDA = MAX(MIN(ICED%XLBDAS_MAX, 10**(6.226 -0.0106*PT(JIJ,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)
+ ZLBDA=MAX(MIN(ICED%XLBDAS_MAX, ICED%XLBS * ( PRHODREF(JIJ,JK) * PRXT(JIJ,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) * &
+ ZWSED(JIJ, JK) = ICEP%XFSEDS * &
+ & PRXT(JIJ,JK)* &
+ & PRHODREF(JIJ,JK)**(1-ICED%XCEXVT) * &
& (1 + (ICED%XFVELOS/ZLBDA)**ICED%XALPHAS)** (-ICED%XNUS+ICEP%XEXSEDS/ICED%XALPHAS) * &
& ZLBDA ** (ICED%XBS+ICEP%XEXSEDS)
@@ -461,45 +449,39 @@ DO WHILE (ANY(ZREMAINT>0.))
CALL PRINT_MSG(NVERB_FATAL, 'GEN', 'ICE4_SEDIMENTATION_SPLIT', 'no sedimentation parameter for KSPE='//TRIM(YSPE) )
END SELECT
!
- ZWSED(:,:,:) = 0.
+ ZWSED(:,:) = 0.
DO JL=1, ISEDIM
- JI=I1(JL)
- JJ=I2(JL)
- JK=I3(JL)
- IF(PRXT(JI,JJ,JK)>ICED%XRTMIN(KSPE)) THEN
- ZWSED(JI, JJ, JK) = ZFSED * PRXT(JI, JJ, JK)**ZEXSED &
- & * PRHODREF(JI, JJ, JK)**(ZEXSED-ICED%XCEXVT)
+ JIJ=I1(JL)
+ JK=I2(JL)
+ IF(PRXT(JIJ,JK)>ICED%XRTMIN(KSPE)) THEN
+ ZWSED(JIJ, JK) = ZFSED * PRXT(JIJ, JK)**ZEXSED &
+ & * PRHODREF(JIJ, JK)**(ZEXSED-ICED%XCEXVT)
ENDIF
ENDDO
ENDIF
- ZMAX_TSTEP(:,:) = ZREMAINT(:,:)
+ ZMAX_TSTEP(:) = ZREMAINT(:)
DO JL=1, ISEDIM
- JI=I1(JL)
- 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), PARAMI%XSPLIT_MAXCFL * PRHODREF(JI, JJ, JK) * &
- & PRXT(JI, JJ, JK) * PDZZ(JI, JJ, JK) / ZWSED(JI, JJ, JK))
+ JIJ=I1(JL)
+ JK=I2(JL)
+ IF(PRXT(JIJ,JK)>ICED%XRTMIN(KSPE) .AND. ZWSED(JIJ, JK)>1.E-20) THEN
+ ZMAX_TSTEP(JIJ) = MIN(ZMAX_TSTEP(JIJ), PARAMI%XSPLIT_MAXCFL * PRHODREF(JIJ, JK) * &
+ & PRXT(JIJ, JK) * PDZZ(JIJ, JK) / ZWSED(JIJ, JK))
ENDIF
ENDDO
- DO JJ = IJB, IJE
- DO JI = IIB, IIE
- ZREMAINT(JI,JJ) = ZREMAINT(JI,JJ) - ZMAX_TSTEP(JI,JJ)
- PINPRX(JI,JJ) = PINPRX(JI,JJ) + ZWSED(JI,JJ,IKB) / CST%XRHOLW * (ZMAX_TSTEP(JI,JJ) * ZINVTSTEP)
- ENDDO
+ DO JIJ = IIJB, IIJE
+ ZREMAINT(JIJ) = ZREMAINT(JIJ) - ZMAX_TSTEP(JIJ)
+ PINPRX(JIJ) = PINPRX(JIJ) + ZWSED(JIJ,IKB) / CST%XRHOLW * (ZMAX_TSTEP(JIJ) * ZINVTSTEP)
ENDDO
DO JK = IKTB , IKTE
- DO JJ = IJB, IJE
- DO JI = IIB, IIE
- ZMRCHANGE = ZMAX_TSTEP(JI,JJ) * POORHODZ(JI,JJ,JK)*(ZWSED(JI,JJ,JK+IKL)-ZWSED(JI,JJ,JK))
- PRXT(JI,JJ,JK) = PRXT(JI,JJ,JK) + ZMRCHANGE + PPRXS(JI,JJ,JK) * ZMAX_TSTEP(JI,JJ)
- PRXS(JI,JJ,JK) = PRXS(JI,JJ,JK) + ZMRCHANGE * ZINVTSTEP
- IF (GPRESENT_PFPR) THEN
- PFPR(JI,JJ,JK,KSPE) = PFPR(JI,JJ,JK,KSPE) + ZWSED(JI,JJ,JK) * (ZMAX_TSTEP(JI,JJ) * ZINVTSTEP)
- ENDIF
- ENDDO
+ DO JIJ = IIJB, IIJE
+ ZMRCHANGE = ZMAX_TSTEP(JIJ) * POORHODZ(JIJ,JK)*(ZWSED(JIJ,JK+IKL)-ZWSED(JIJ,JK))
+ PRXT(JIJ,JK) = PRXT(JIJ,JK) + ZMRCHANGE + PPRXS(JIJ,JK) * ZMAX_TSTEP(JIJ)
+ PRXS(JIJ,JK) = PRXS(JIJ,JK) + ZMRCHANGE * ZINVTSTEP
+ IF (GPRESENT_PFPR) THEN
+ PFPR(JIJ,JK,KSPE) = PFPR(JIJ,JK,KSPE) + ZWSED(JIJ,JK) * (ZMAX_TSTEP(JIJ) * ZINVTSTEP)
+ ENDIF
ENDDO
ENDDO
!
diff --git a/src/common/micro/mode_ice4_sedimentation_stat.F90 b/src/common/micro/mode_ice4_sedimentation_stat.F90
index 33e4e6e01e7ff5992ce919206038b62e0647a2e3..5d64266c1c84008bcaad0f3cbcb54baeb0f8221b 100644
--- a/src/common/micro/mode_ice4_sedimentation_stat.F90
+++ b/src/common/micro/mode_ice4_sedimentation_stat.F90
@@ -6,10 +6,9 @@
MODULE MODE_ICE4_SEDIMENTATION_STAT
IMPLICIT NONE
CONTAINS
-SUBROUTINE ICE4_SEDIMENTATION_STAT(D, CST, ICEP, ICED, &
- &PTSTEP, KRR, OSEDIC, PDZZ, &
- &PRHODREF, PPABST, PTHT, PRHODJ, &
- &PLBDAS, &
+SUBROUTINE ICE4_SEDIMENTATION_STAT(D, CST, ICEP, ICED, PARAMI, &
+ &PTSTEP, KRR, PDZZ, &
+ &PRHODREF, PPABST, PTHT, PT, PRHODJ, &
&PRCS, PRCT, PRRS, PRRT, PRIS, PRIT, &
&PRSS, PRST, PRGS, PRGT,&
&PINPRC, PINPRR, PINPRI, PINPRS, PINPRG, &
@@ -45,6 +44,7 @@ USE MODD_DIMPHYEX, ONLY: DIMPHYEX_t
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
USE MODI_GAMMA, ONLY: GAMMA
!
IMPLICIT NONE
@@ -55,46 +55,46 @@ TYPE(DIMPHYEX_t), INTENT(IN) :: D !array dimensio
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
REAL, INTENT(IN) :: PTSTEP ! Double Time step (single if cold start)
INTEGER, INTENT(IN) :: KRR ! Number of moist variable
-LOGICAL, INTENT(IN) :: OSEDIC ! Switch for droplet sedim.
-REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PDZZ ! Layer thikness (m)
-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) :: 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
-REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PRRT ! Rain water m.r. at t
-REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(INOUT) :: PRIS ! Pristine ice m.r. source
-REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PRIT ! Pristine ice m.r. at t
-REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(INOUT) :: PRSS ! Snow/aggregate m.r. source
-REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PRST ! Snow/aggregate m.r. at t
-REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(INOUT) :: PRGS ! Graupel m.r. source
-REAL, DIMENSION(D%NIT,D%NJT,D%NKT), INTENT(IN) :: PRGT ! Graupel/hail m.r. at t
-REAL, DIMENSION(D%NIT,D%NJT), INTENT(OUT) :: PINPRC ! Cloud instant precip
-REAL, DIMENSION(D%NIT,D%NJT), INTENT(OUT) :: PINPRR ! Rain instant precip
-REAL, DIMENSION(D%NIT,D%NJT), INTENT(OUT) :: PINPRI ! Pristine ice instant precip
-REAL, DIMENSION(D%NIT,D%NJT), INTENT(OUT) :: PINPRS ! Snow instant precip
-REAL, DIMENSION(D%NIT,D%NJT), INTENT(OUT) :: PINPRG ! Graupel instant precip
-REAL, DIMENSION(D%NIT,D%NJT), OPTIONAL, INTENT(IN) :: PSEA ! Sea Mask
-REAL, DIMENSION(D%NIT,D%NJT), OPTIONAL, INTENT(IN) :: PTOWN ! Fraction that is town
-REAL, DIMENSION(D%NIT,D%NJT), OPTIONAL, INTENT(OUT) :: PINPRH ! Hail instant precip
-REAL, DIMENSION(D%NIT,D%NJT,D%NKT), OPTIONAL, INTENT(IN) :: PRHT ! Hail m.r. at t
-REAL, DIMENSION(D%NIT,D%NJT,D%NKT), OPTIONAL, INTENT(INOUT) :: PRHS ! Hail m.r. source
-REAL, DIMENSION(D%NIT,D%NJT,D%NKT,KRR), OPTIONAL, INTENT(OUT) :: PFPR ! upper-air precipitation fluxes
+REAL, DIMENSION(D%NIJT,D%NKT), INTENT(IN) :: PDZZ ! Layer thikness (m)
+REAL, DIMENSION(D%NIJT,D%NKT), INTENT(IN) :: PRHODREF! Reference density
+REAL, DIMENSION(D%NIJT,D%NKT), INTENT(IN) :: PPABST ! absolute pressure at t
+REAL, DIMENSION(D%NIJT,D%NKT), INTENT(IN) :: PTHT ! Theta at time t
+REAL, DIMENSION(D%NIJT,D%NKT), INTENT(IN) :: PT ! Temperature
+REAL, DIMENSION(D%NIJT,D%NKT), INTENT(IN) :: PRHODJ ! Dry density * Jacobian
+REAL, DIMENSION(D%NIJT,D%NKT), INTENT(INOUT) :: PRCS ! Cloud water m.r. source
+REAL, DIMENSION(D%NIJT,D%NKT), INTENT(IN) :: PRCT ! Cloud water m.r. at t
+REAL, DIMENSION(D%NIJT,D%NKT), INTENT(INOUT) :: PRRS ! Rain water m.r. source
+REAL, DIMENSION(D%NIJT,D%NKT), INTENT(IN) :: PRRT ! Rain water m.r. at t
+REAL, DIMENSION(D%NIJT,D%NKT), INTENT(INOUT) :: PRIS ! Pristine ice m.r. source
+REAL, DIMENSION(D%NIJT,D%NKT), INTENT(IN) :: PRIT ! Pristine ice m.r. at t
+REAL, DIMENSION(D%NIJT,D%NKT), INTENT(INOUT) :: PRSS ! Snow/aggregate m.r. source
+REAL, DIMENSION(D%NIJT,D%NKT), INTENT(IN) :: PRST ! Snow/aggregate m.r. at t
+REAL, DIMENSION(D%NIJT,D%NKT), INTENT(INOUT) :: PRGS ! Graupel m.r. source
+REAL, DIMENSION(D%NIJT,D%NKT), INTENT(IN) :: PRGT ! Graupel/hail m.r. at t
+REAL, DIMENSION(D%NIJT), INTENT(OUT) :: PINPRC ! Cloud instant precip
+REAL, DIMENSION(D%NIJT), INTENT(OUT) :: PINPRR ! Rain instant precip
+REAL, DIMENSION(D%NIJT), INTENT(OUT) :: PINPRI ! Pristine ice instant precip
+REAL, DIMENSION(D%NIJT), INTENT(OUT) :: PINPRS ! Snow instant precip
+REAL, DIMENSION(D%NIJT), INTENT(OUT) :: PINPRG ! Graupel instant precip
+REAL, DIMENSION(D%NIJT), OPTIONAL, INTENT(IN) :: PSEA ! Sea Mask
+REAL, DIMENSION(D%NIJT), OPTIONAL, INTENT(IN) :: PTOWN ! Fraction that is town
+REAL, DIMENSION(D%NIJT), OPTIONAL, INTENT(OUT) :: PINPRH ! Hail instant precip
+REAL, DIMENSION(D%NIJT,D%NKT), OPTIONAL, INTENT(IN) :: PRHT ! Hail m.r. at t
+REAL, DIMENSION(D%NIJT,D%NKT), OPTIONAL, INTENT(INOUT) :: PRHS ! Hail m.r. source
+REAL, DIMENSION(D%NIJT,D%NKT,KRR), OPTIONAL, INTENT(OUT) :: PFPR ! upper-air precipitation fluxes
!
!* 0.2 declaration of local variables
!
LOGICAL :: LLSEA_AND_TOWN
-INTEGER :: JRR, JI, JJ, JK, IKB, IKE,IKL, IIE, IIB, IJB, IJE, IKTB, IKTE
+INTEGER :: JRR, JIJ, JK, IKB, IKE,IKL, IIJB, IIJE, IKTB, IKTE
INTEGER :: ISHIFT, IK, IKPLUS
-REAL :: ZQP, ZP1, ZINVTSTEP, ZGAC, ZGC, ZGAC2, ZGC2, ZRAYDEFO
-REAL, DIMENSION(D%NIT) :: ZWSEDW1, ZWSEDW2 ! sedimentation speed
-REAL, DIMENSION(D%NIT,D%NJT) :: ZTSORHODZ ! TimeStep Over (Rhodref times delta Z)
-REAL, DIMENSION(D%NIT,D%NJT,0:1,2:KRR) :: ZSED ! sedimentation flux array for each species and for above and current levels
+REAL :: ZQP, ZP1, ZINVTSTEP, ZGAC, ZGC, ZGAC2, ZGC2, ZRAYDEFO, ZLBDAS
+REAL, DIMENSION(D%NIJT) :: ZWSEDW1, ZWSEDW2 ! sedimentation speed
+REAL, DIMENSION(D%NIJT) :: ZTSORHODZ ! TimeStep Over (Rhodref times delta Z)
+REAL, DIMENSION(D%NIJT,0:1,2:KRR) :: ZSED ! sedimentation flux array for each species and for above and current levels
REAL :: FWSED1, FWSED2, PWSEDW, PWSEDWSUP, PINVTSTEP, PTSTEP1, PDZZ1, PRHODREF1, PRXT1
REAL(KIND=JPRB) :: ZHOOK_HANDLE
@@ -115,17 +115,15 @@ IF (LHOOK) CALL DR_HOOK('ICE4_SEDIMENTATION_STAT',0,ZHOOK_HANDLE)
IKB=D%NKB
IKE=D%NKE
IKL=D%NKL
-IIB=D%NIB
-IIE=D%NIE
-IJB=D%NJB
-IJE=D%NJE
+IIJB=D%NIJB
+IIJE=D%NIJE
IKTB=D%NKTB
IKTE=D%NKTE
!
IF ( PRESENT( PFPR ) ) THEN
!Set to 0. to avoid undefined values (in files)
- PFPR(:, :, : IKTB, :) = 0.
- PFPR(:, :, IKTE :, :) = 0.
+ PFPR(:, : IKTB, :) = 0.
+ PFPR(:, IKTE :, :) = 0.
END IF
!-------------------------------------------------------------------------------
@@ -149,16 +147,14 @@ CALL SHIFT
! Initialize vertical loop
DO JRR=2,KRR
- ZSED(:,:,IKPLUS,JRR) = 0.
+ ZSED(:,IKPLUS,JRR) = 0.
ENDDO
! calculation sedimentation flux
DO JK = IKE , IKB, -1*IKL
- DO JJ = IJB, IJE
- DO JI = IIB, IIE
- ZTSORHODZ(JI,JJ) =PTSTEP/(PRHODREF(JI,JJ,JK)*PDZZ(JI,JJ,JK))
- ENDDO
+ DO JIJ = IIJB, IIJE
+ ZTSORHODZ(JIJ) =PTSTEP/(PRHODREF(JIJ,JK)*PDZZ(JIJ,JK))
ENDDO
!
DO JRR=2,KRR
@@ -166,40 +162,40 @@ DO JK = IKE , IKB, -1*IKL
IF (JRR==2) THEN
!******* for cloud
- IF (OSEDIC) THEN
- CALL CLOUD(PRCT(:,:,JK))
+ IF (PARAMI%LSEDIC) THEN
+ CALL CLOUD(PRCT(:,JK))
ELSE
- ZSED(:,:,IK,JRR)=0.
+ ZSED(:,IK,JRR)=0.
ENDIF
ELSEIF (JRR==3) THEN
!* 2.2 for rain
- CALL OTHER_SPECIES(ICEP%XFSEDR,ICEP%XEXSEDR,PRRT(:,:,JK))
+ CALL OTHER_SPECIES(ICEP%XFSEDR,ICEP%XEXSEDR,PRRT(:,JK))
ELSEIF (JRR==4) THEN
- CALL PRISTINE_ICE(PRIT(:,:,JK))
+ CALL PRISTINE_ICE(PRIT(:,JK))
ELSEIF (JRR==5) THEN
!* 2.4 for aggregates/snow
#ifdef REPRO48
- CALL OTHER_SPECIES(ICEP%XFSEDS,ICEP%XEXSEDS,PRST(:,:,JK))
+ CALL OTHER_SPECIES(ICEP%XFSEDS,ICEP%XEXSEDS,PRST(:,JK))
#else
- CALL SNOW(PRST(:,:,JK))
+ CALL SNOW(PRST(:,JK))
#endif
ELSEIF (JRR==6) THEN
!* 2.5 for graupeln
- CALL OTHER_SPECIES(ICEP%XFSEDG,ICEP%XEXSEDG,PRGT(:,:,JK))
+ CALL OTHER_SPECIES(ICEP%XFSEDG,ICEP%XEXSEDG,PRGT(:,JK))
ELSEIF (JRR==7) THEN
!* 2.6 for hail
IF (PRESENT(PRHT)) THEN
- CALL OTHER_SPECIES(ICEP%XFSEDH,ICEP%XEXSEDH,PRHT(:,:,JK))
+ CALL OTHER_SPECIES(ICEP%XFSEDH,ICEP%XEXSEDH,PRHT(:,JK))
ENDIF
ENDIF
@@ -210,35 +206,31 @@ DO JK = IKE , IKB, -1*IKL
IF(PRESENT(PFPR)) THEN
DO JRR=2,KRR
- PFPR(:,:,JK,JRR)=ZSED(:,:,IK,JRR)
+ PFPR(:,JK,JRR)=ZSED(:,IK,JRR)
ENDDO
ENDIF
- DO JJ = IJB, IJE
- DO JI = IIB, IIE
- PRCS(JI,JJ,JK) = PRCS(JI,JJ,JK)+ZTSORHODZ(JI,JJ)*(ZSED(JI,JJ,IKPLUS,2)-ZSED(JI,JJ,IK,2))*ZINVTSTEP
- PRRS(JI,JJ,JK) = PRRS(JI,JJ,JK)+ZTSORHODZ(JI,JJ)*(ZSED(JI,JJ,IKPLUS,3)-ZSED(JI,JJ,IK,3))*ZINVTSTEP
- PRIS(JI,JJ,JK) = PRIS(JI,JJ,JK)+ZTSORHODZ(JI,JJ)*(ZSED(JI,JJ,IKPLUS,4)-ZSED(JI,JJ,IK,4))*ZINVTSTEP
- PRSS(JI,JJ,JK) = PRSS(JI,JJ,JK)+ZTSORHODZ(JI,JJ)*(ZSED(JI,JJ,IKPLUS,5)-ZSED(JI,JJ,IK,5))*ZINVTSTEP
- PRGS(JI,JJ,JK) = PRGS(JI,JJ,JK)+ZTSORHODZ(JI,JJ)*(ZSED(JI,JJ,IKPLUS,6)-ZSED(JI,JJ,IK,6))*ZINVTSTEP
- IF (PRESENT(PRHS)) THEN
- PRHS(JI,JJ,JK) = PRHS(JI,JJ,JK)+ZTSORHODZ(JI,JJ)*(ZSED(JI,JJ,IKPLUS,7)-ZSED(JI,JJ,IK,7))*ZINVTSTEP
- ENDIF
- ENDDO
+ DO JIJ = IIJB, IIJE
+ PRCS(JIJ,JK) = PRCS(JIJ,JK)+ZTSORHODZ(JIJ)*(ZSED(JIJ,IKPLUS,2)-ZSED(JIJ,IK,2))*ZINVTSTEP
+ PRRS(JIJ,JK) = PRRS(JIJ,JK)+ZTSORHODZ(JIJ)*(ZSED(JIJ,IKPLUS,3)-ZSED(JIJ,IK,3))*ZINVTSTEP
+ PRIS(JIJ,JK) = PRIS(JIJ,JK)+ZTSORHODZ(JIJ)*(ZSED(JIJ,IKPLUS,4)-ZSED(JIJ,IK,4))*ZINVTSTEP
+ PRSS(JIJ,JK) = PRSS(JIJ,JK)+ZTSORHODZ(JIJ)*(ZSED(JIJ,IKPLUS,5)-ZSED(JIJ,IK,5))*ZINVTSTEP
+ PRGS(JIJ,JK) = PRGS(JIJ,JK)+ZTSORHODZ(JIJ)*(ZSED(JIJ,IKPLUS,6)-ZSED(JIJ,IK,6))*ZINVTSTEP
+ IF (PRESENT(PRHS)) THEN
+ PRHS(JIJ,JK) = PRHS(JIJ,JK)+ZTSORHODZ(JIJ)*(ZSED(JIJ,IKPLUS,7)-ZSED(JIJ,IK,7))*ZINVTSTEP
+ ENDIF
ENDDO
IF (JK==IKB) THEN
- DO JJ = IJB, IJE
- DO JI = IIB, IIE
- IF(OSEDIC) PINPRC(JI,JJ) = ZSED(JI,JJ,IK,2)/CST%XRHOLW
- PINPRR(JI,JJ) = ZSED(JI,JJ,IK,3)/CST%XRHOLW
- PINPRI(JI,JJ) = ZSED(JI,JJ,IK,4)/CST%XRHOLW
- PINPRS(JI,JJ) = ZSED(JI,JJ,IK,5)/CST%XRHOLW
- PINPRG(JI,JJ) = ZSED(JI,JJ,IK,6)/CST%XRHOLW
- IF (PRESENT(PINPRH)) THEN
- PINPRH(JI,JJ) = ZSED(JI,JJ,IK,7)/CST%XRHOLW
- ENDIF
- ENDDO
+ DO JIJ = IIJB, IIJE
+ IF(PARAMI%LSEDIC) PINPRC(JIJ) = ZSED(JIJ,IK,2)/CST%XRHOLW
+ PINPRR(JIJ) = ZSED(JIJ,IK,3)/CST%XRHOLW
+ PINPRI(JIJ) = ZSED(JIJ,IK,4)/CST%XRHOLW
+ PINPRS(JIJ) = ZSED(JIJ,IK,5)/CST%XRHOLW
+ PINPRG(JIJ) = ZSED(JIJ,IK,6)/CST%XRHOLW
+ IF (PRESENT(PINPRH)) THEN
+ PINPRH(JIJ) = ZSED(JIJ,IK,7)/CST%XRHOLW
+ ENDIF
ENDDO
ENDIF
@@ -253,7 +245,7 @@ CONTAINS
SUBROUTINE CLOUD(PRXT)
- REAL, INTENT(IN) :: PRXT(D%NIT,D%NJT) ! mr of specy X
+ REAL, INTENT(IN) :: PRXT(D%NIJT) ! mr of specy X
REAL :: ZLBC ! XLBC weighted by sea fraction
REAL :: ZFSEDC
@@ -265,53 +257,51 @@ CONTAINS
!!IF (LHOOK) CALL DR_HOOK('ICE4_SEDIMENTATION_STAT:CLOUD',0,ZHOOK_HANDLE)
- DO JJ = IJB, IJE
- DO JI = IIB, IIE
- !estimation of q' taking into account incoming ZWSED from previous vertical level
- ZQP=ZSED(JI,JJ,IKPLUS,JRR)*ZTSORHODZ(JI,JJ)
- IF ((PRXT(JI,JJ) > ICED%XRTMIN(JRR)) .OR. (ZQP > ICED%XRTMIN(JRR))) THEN
- IF (LLSEA_AND_TOWN) THEN
- ZRAY = MAX(1.,0.5*((1.-PSEA(JI,JJ))*ZGAC/ZGC+PSEA(JI,JJ)*ZGAC2/ZGC2))
- ZLBC = MAX(MIN(ICED%XLBC(1),ICED%XLBC(2)),(PSEA(JI,JJ)*ICED%XLBC(2)+(1.-PSEA(JI,JJ))*ICED%XLBC(1)) )
- ZFSEDC = MAX(MIN(ICEP%XFSEDC(1),ICEP%XFSEDC(2)), (PSEA(JI,JJ)*ICEP%XFSEDC(2)+(1.-PSEA(JI,JJ))*ICEP%XFSEDC(1)) )
- ZCONC3D= (1.-PTOWN(JI,JJ))*(PSEA(JI,JJ)*ICED%XCONC_SEA+(1.-PSEA(JI,JJ))*ICED%XCONC_LAND) + &
- PTOWN(JI,JJ) *ICED%XCONC_URBAN
- ELSE
- ZRAY = ZRAYDEFO
- ZLBC = ICED%XLBC(1)
- ZFSEDC = ICEP%XFSEDC(1)
- ZCONC3D= ICED%XCONC_LAND
- ENDIF
- !calculation of w
- IF(PRXT(JI,JJ) > ICED%XRTMIN(JRR)) THEN
- ZZWLBDA=6.6E-8*(101325./PPABST(JI,JJ,JK))*(PTHT(JI,JJ,JK)/293.15)
- ZZWLBDC=(ZLBC*ZCONC3D/(PRHODREF(JI,JJ,JK)*PRXT(JI,JJ)))**ICED%XLBEXC
- ZZCC=ICED%XCC*(1.+1.26*ZZWLBDA*ZZWLBDC/ZRAY) !! ZCC : Fall speed
- ZWSEDW1(JI)=PRHODREF(JI,JJ,JK)**(-ICED%XCEXVT ) * ZZWLBDC**(-ICED%XDC)*ZZCC*ZFSEDC
- ELSE
- ZWSEDW1(JI)=0.
- ENDIF
- IF ( ZQP > ICED%XRTMIN(JRR) ) THEN
- ZZWLBDA=6.6E-8*(101325./PPABST(JI,JJ,JK))*(PTHT(JI,JJ,JK)/293.15)
- ZZWLBDC=(ZLBC*ZCONC3D/(PRHODREF(JI,JJ,JK)*ZQP))**ICED%XLBEXC
- ZZCC=ICED%XCC*(1.+1.26*ZZWLBDA*ZZWLBDC/ZRAY) !! ZCC : Fall speed
- ZWSEDW2(JI)=PRHODREF(JI,JJ,JK)**(-ICED%XCEXVT ) * ZZWLBDC**(-ICED%XDC)*ZZCC*ZFSEDC
- ELSE
- ZWSEDW2(JI)=0.
- ENDIF
+ DO JIJ = IIJB, IIJE
+ !estimation of q' taking into account incoming ZWSED from previous vertical level
+ ZQP=ZSED(JIJ,IKPLUS,JRR)*ZTSORHODZ(JIJ)
+ IF ((PRXT(JIJ) > ICED%XRTMIN(JRR)) .OR. (ZQP > ICED%XRTMIN(JRR))) THEN
+ IF (LLSEA_AND_TOWN) THEN
+ ZRAY = MAX(1.,0.5*((1.-PSEA(JIJ))*ZGAC/ZGC+PSEA(JIJ)*ZGAC2/ZGC2))
+ ZLBC = MAX(MIN(ICED%XLBC(1),ICED%XLBC(2)),(PSEA(JIJ)*ICED%XLBC(2)+(1.-PSEA(JIJ))*ICED%XLBC(1)) )
+ ZFSEDC = MAX(MIN(ICEP%XFSEDC(1),ICEP%XFSEDC(2)), (PSEA(JIJ)*ICEP%XFSEDC(2)+(1.-PSEA(JIJ))*ICEP%XFSEDC(1)) )
+ ZCONC3D= (1.-PTOWN(JIJ))*(PSEA(JIJ)*ICED%XCONC_SEA+(1.-PSEA(JIJ))*ICED%XCONC_LAND) + &
+ PTOWN(JIJ) *ICED%XCONC_URBAN
ELSE
- ZWSEDW1(JI)=0.
- ZWSEDW2(JI)=0.
+ ZRAY = ZRAYDEFO
+ ZLBC = ICED%XLBC(1)
+ ZFSEDC = ICEP%XFSEDC(1)
+ ZCONC3D= ICED%XCONC_LAND
ENDIF
-!- duplicated code -------------------------------------------------------------------------
- IF (ZWSEDW2(JI) /= 0.) THEN
- ZSED(JI,JJ,IK,JRR)=FWSED1(ZWSEDW1(JI),PTSTEP,PDZZ(JI,JJ,JK),PRHODREF(JI,JJ,JK),PRXT(JI,JJ),ZINVTSTEP) &
- & + FWSED2(ZWSEDW2(JI),PTSTEP,PDZZ(JI,JJ,JK),ZSED(JI,JJ,IKPLUS,JRR))
+ !calculation of w
+ IF(PRXT(JIJ) > ICED%XRTMIN(JRR)) THEN
+ ZZWLBDA=6.6E-8*(101325./PPABST(JIJ,JK))*(PTHT(JIJ,JK)/293.15)
+ ZZWLBDC=(ZLBC*ZCONC3D/(PRHODREF(JIJ,JK)*PRXT(JIJ)))**ICED%XLBEXC
+ ZZCC=ICED%XCC*(1.+1.26*ZZWLBDA*ZZWLBDC/ZRAY) !! ZCC : Fall speed
+ ZWSEDW1(JIJ)=PRHODREF(JIJ,JK)**(-ICED%XCEXVT ) * ZZWLBDC**(-ICED%XDC)*ZZCC*ZFSEDC
+ ELSE
+ ZWSEDW1(JIJ)=0.
+ ENDIF
+ IF ( ZQP > ICED%XRTMIN(JRR) ) THEN
+ ZZWLBDA=6.6E-8*(101325./PPABST(JIJ,JK))*(PTHT(JIJ,JK)/293.15)
+ ZZWLBDC=(ZLBC*ZCONC3D/(PRHODREF(JIJ,JK)*ZQP))**ICED%XLBEXC
+ ZZCC=ICED%XCC*(1.+1.26*ZZWLBDA*ZZWLBDC/ZRAY) !! ZCC : Fall speed
+ ZWSEDW2(JIJ)=PRHODREF(JIJ,JK)**(-ICED%XCEXVT ) * ZZWLBDC**(-ICED%XDC)*ZZCC*ZFSEDC
ELSE
- ZSED(JI,JJ,IK,JRR)=FWSED1(ZWSEDW1(JI),PTSTEP,PDZZ(JI,JJ,JK),PRHODREF(JI,JJ,JK),PRXT(JI,JJ),ZINVTSTEP)
+ ZWSEDW2(JIJ)=0.
ENDIF
+ ELSE
+ ZWSEDW1(JIJ)=0.
+ ZWSEDW2(JIJ)=0.
+ ENDIF
+!- duplicated code -------------------------------------------------------------------------
+ IF (ZWSEDW2(JIJ) /= 0.) THEN
+ ZSED(JIJ,IK,JRR)=FWSED1(ZWSEDW1(JIJ),PTSTEP,PDZZ(JIJ,JK),PRHODREF(JIJ,JK),PRXT(JIJ),ZINVTSTEP) &
+ & + FWSED2(ZWSEDW2(JIJ),PTSTEP,PDZZ(JIJ,JK),ZSED(JIJ,IKPLUS,JRR))
+ ELSE
+ ZSED(JIJ,IK,JRR)=FWSED1(ZWSEDW1(JIJ),PTSTEP,PDZZ(JIJ,JK),PRHODREF(JIJ,JK),PRXT(JIJ),ZINVTSTEP)
+ ENDIF
!-------------------------------------------------------------------------------------------
- ENDDO
ENDDO
!!IF (LHOOK) CALL DR_HOOK('ICE4_SEDIMENTATION_STAT:CLOUD',1,ZHOOK_HANDLE)
@@ -320,47 +310,45 @@ CONTAINS
SUBROUTINE PRISTINE_ICE(PRXT)
- REAL, INTENT(IN) :: PRXT(D%NIT,D%NJT) ! mr of specy X
+ REAL, INTENT(IN) :: PRXT(D%NIJT) ! mr of specy X
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!!IF (LHOOK) CALL DR_HOOK('ICE4_SEDIMENTATION_STAT:PRISTINE_ICE',0,ZHOOK_HANDLE)
! ******* for pristine ice
- DO JJ = IJB, IJE
- DO JI = IIB, IIE
- ZQP=ZSED(JI,JJ,IKPLUS,JRR)*ZTSORHODZ(JI,JJ)
- IF ((PRXT(JI,JJ) > ICED%XRTMIN(JRR)) .OR. (ZQP > ICED%XRTMIN(JRR))) THEN
- !calculation of w
- IF ( PRXT(JI,JJ) > MAX(ICED%XRTMIN(JRR),1.0E-7 ) ) THEN
- ZWSEDW1(JI)= ICEP%XFSEDI * &
- & PRHODREF(JI,JJ,JK)**(-ICED%XCEXVT) * & ! McF&H
- & MAX( 0.05E6,-0.15319E6-0.021454E6* &
- & ALOG(PRHODREF(JI,JJ,JK)*PRXT(JI,JJ)) )**ICEP%XEXCSEDI
- ELSE
- ZWSEDW1(JI)=0.
- ENDIF
- IF ( ZQP > MAX(ICED%XRTMIN(JRR),1.0E-7 ) ) THEN
- ZWSEDW2(JI)= ICEP%XFSEDI * &
- & PRHODREF(JI,JJ,JK)**(-ICED%XCEXVT) * & ! McF&H
- & MAX( 0.05E6,-0.15319E6-0.021454E6* &
- & ALOG(PRHODREF(JI,JJ,JK)*ZQP) )**ICEP%XEXCSEDI
- ELSE
- ZWSEDW2(JI)=0.
- ENDIF
+ DO JIJ = IIJB, IIJE
+ ZQP=ZSED(JIJ,IKPLUS,JRR)*ZTSORHODZ(JIJ)
+ IF ((PRXT(JIJ) > ICED%XRTMIN(JRR)) .OR. (ZQP > ICED%XRTMIN(JRR))) THEN
+ !calculation of w
+ IF ( PRXT(JIJ) > MAX(ICED%XRTMIN(JRR),1.0E-7 ) ) THEN
+ ZWSEDW1(JIJ)= ICEP%XFSEDI * &
+ & PRHODREF(JIJ,JK)**(-ICED%XCEXVT) * & ! McF&H
+ & MAX( 0.05E6,-0.15319E6-0.021454E6* &
+ & ALOG(PRHODREF(JIJ,JK)*PRXT(JIJ)) )**ICEP%XEXCSEDI
ELSE
- ZWSEDW1(JI)=0.
- ZWSEDW2(JI)=0.
+ ZWSEDW1(JIJ)=0.
ENDIF
-!- duplicated code -------------------------------------------------------------------------
- IF (ZWSEDW2(JI) /= 0.) THEN
- ZSED(JI,JJ,IK,JRR)=FWSED1(ZWSEDW1(JI),PTSTEP,PDZZ(JI,JJ,JK),PRHODREF(JI,JJ,JK),PRXT(JI,JJ),ZINVTSTEP) &
- & + FWSED2(ZWSEDW2(JI),PTSTEP,PDZZ(JI,JJ,JK),ZSED(JI,JJ,IKPLUS,JRR))
+ IF ( ZQP > MAX(ICED%XRTMIN(JRR),1.0E-7 ) ) THEN
+ ZWSEDW2(JIJ)= ICEP%XFSEDI * &
+ & PRHODREF(JIJ,JK)**(-ICED%XCEXVT) * & ! McF&H
+ & MAX( 0.05E6,-0.15319E6-0.021454E6* &
+ & ALOG(PRHODREF(JIJ,JK)*ZQP) )**ICEP%XEXCSEDI
ELSE
- ZSED(JI,JJ,IK,JRR)=FWSED1(ZWSEDW1(JI),PTSTEP,PDZZ(JI,JJ,JK),PRHODREF(JI,JJ,JK),PRXT(JI,JJ),ZINVTSTEP)
+ ZWSEDW2(JIJ)=0.
ENDIF
+ ELSE
+ ZWSEDW1(JIJ)=0.
+ ZWSEDW2(JIJ)=0.
+ ENDIF
+!- duplicated code -------------------------------------------------------------------------
+ IF (ZWSEDW2(JIJ) /= 0.) THEN
+ ZSED(JIJ,IK,JRR)=FWSED1(ZWSEDW1(JIJ),PTSTEP,PDZZ(JIJ,JK),PRHODREF(JIJ,JK),PRXT(JIJ),ZINVTSTEP) &
+ & + FWSED2(ZWSEDW2(JIJ),PTSTEP,PDZZ(JIJ,JK),ZSED(JIJ,IKPLUS,JRR))
+ ELSE
+ ZSED(JIJ,IK,JRR)=FWSED1(ZWSEDW1(JIJ),PTSTEP,PDZZ(JIJ,JK),PRHODREF(JIJ,JK),PRXT(JIJ),ZINVTSTEP)
+ ENDIF
!-------------------------------------------------------------------------------------------
- ENDDO
ENDDO
!!IF (LHOOK) CALL DR_HOOK('ICE4_SEDIMENTATION_STAT:PRISTINE_ICE',1,ZHOOK_HANDLE)
@@ -369,47 +357,55 @@ CONTAINS
SUBROUTINE SNOW(PRXT)
- REAL, INTENT(IN) :: PRXT(D%NIT,D%NJT) ! mr of specy X
+ REAL, INTENT(IN) :: PRXT(D%NIJT) ! mr of specy X
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!!IF (LHOOK) CALL DR_HOOK('ICE4_SEDIMENTATION_STAT:SNOW',0,ZHOOK_HANDLE)
! ******* for snow
- DO JJ = IJB, IJE
- DO JI = IIB, IIE
- ZQP=ZSED(JI,JJ,IKPLUS,JRR)*ZTSORHODZ(JI,JJ)
- IF ((PRXT(JI,JJ) > ICED%XRTMIN(JRR)) .OR. (ZQP > ICED%XRTMIN(JRR))) THEN
- !calculation of w
- IF ( PRXT(JI,JJ) > ICED%XRTMIN(JRR) ) THEN
- ZWSEDW1(JI)= ICEP%XFSEDS * &
- & PRHODREF(JI,JJ,JK)**(-ICED%XCEXVT) * &
- & (1+(ICED%XFVELOS/PLBDAS(JI,JJ,JK))**ICED%XALPHAS)**(-ICED%XNUS+ICEP%XEXSEDS/ICED%XALPHAS)* &
- & PLBDAS(JI,JJ,JK)**(ICED%XBS+ICEP%XEXSEDS)
- ELSE
- ZWSEDW1(JI)=0.
- ENDIF
- IF ( ZQP > ICED%XRTMIN(JRR) ) THEN
- ZWSEDW2(JI)= ICEP%XFSEDS * &
- & PRHODREF(JI,JJ,JK)**(-ICED%XCEXVT) * &
- & (1+(ICED%XFVELOS/PLBDAS(JI,JJ,JK))**ICED%XALPHAS)**(-ICED%XNUS+ICEP%XEXSEDS/ICED%XALPHAS)* &
- & PLBDAS(JI,JJ,JK)**(ICED%XBS+ICEP%XEXSEDS)
+ DO JIJ = IIJB, IIJE
+ ZQP=ZSED(JIJ,IKPLUS,JRR)*ZTSORHODZ(JIJ)
+ IF ((PRXT(JIJ) > ICED%XRTMIN(JRR)) .OR. (ZQP > ICED%XRTMIN(JRR))) THEN
+ !Compute lambda_snow parameter
+ IF (PARAMI%LSNOW_T) THEN
+ IF(PT(JIJ,JK)>CST%XTT-10.0) THEN
+ ZLBDAS = MAX(MIN(ICED%XLBDAS_MAX, 10**(14.554-0.0423*PT(JIJ,JK))),ICED%XLBDAS_MIN)*ICED%XTRANS_MP_GAMMAS
ELSE
- ZWSEDW2(JI)=0.
- ENDIF
+ ZLBDAS = MAX(MIN(ICED%XLBDAS_MAX, 10**(6.226-0.0106*PT(JIJ,JK))),ICED%XLBDAS_MIN)*ICED%XTRANS_MP_GAMMAS
+ END IF
ELSE
- ZWSEDW1(JI)=0.
- ZWSEDW2(JI)=0.
+ ZLBDAS = MAX(MIN(ICED%XLBDAS_MAX,ICED%XLBS*(PRHODREF(JIJ,JK)*PRXT(JIJ))**ICED%XLBEXS),ICED%XLBDAS_MIN)
+ END IF
+ !calculation of w
+ IF ( PRXT(JIJ) > ICED%XRTMIN(JRR) ) THEN
+ ZWSEDW1(JIJ)= ICEP%XFSEDS * &
+ & PRHODREF(JIJ,JK)**(-ICED%XCEXVT) * &
+ & (1+(ICED%XFVELOS/ZLBDAS)**ICED%XALPHAS)**(-ICED%XNUS+ICEP%XEXSEDS/ICED%XALPHAS)* &
+ & ZLBDAS**(ICED%XBS+ICEP%XEXSEDS)
+ ELSE
+ ZWSEDW1(JIJ)=0.
ENDIF
-!- duplicated code -------------------------------------------------------------------------
- IF (ZWSEDW2(JI) /= 0.) THEN
- ZSED(JI,JJ,IK,JRR)=FWSED1(ZWSEDW1(JI),PTSTEP,PDZZ(JI,JJ,JK),PRHODREF(JI,JJ,JK),PRXT(JI,JJ),ZINVTSTEP) &
- & + FWSED2(ZWSEDW2(JI),PTSTEP,PDZZ(JI,JJ,JK),ZSED(JI,JJ,IKPLUS,JRR))
+ IF ( ZQP > ICED%XRTMIN(JRR) ) THEN
+ ZWSEDW2(JIJ)= ICEP%XFSEDS * &
+ & PRHODREF(JIJ,JK)**(-ICED%XCEXVT) * &
+ & (1+(ICED%XFVELOS/ZLBDAS)**ICED%XALPHAS)**(-ICED%XNUS+ICEP%XEXSEDS/ICED%XALPHAS)* &
+ & ZLBDAS**(ICED%XBS+ICEP%XEXSEDS)
ELSE
- ZSED(JI,JJ,IK,JRR)=FWSED1(ZWSEDW1(JI),PTSTEP,PDZZ(JI,JJ,JK),PRHODREF(JI,JJ,JK),PRXT(JI,JJ),ZINVTSTEP)
+ ZWSEDW2(JIJ)=0.
ENDIF
+ ELSE
+ ZWSEDW1(JIJ)=0.
+ ZWSEDW2(JIJ)=0.
+ ENDIF
+!- duplicated code -------------------------------------------------------------------------
+ IF (ZWSEDW2(JIJ) /= 0.) THEN
+ ZSED(JIJ,IK,JRR)=FWSED1(ZWSEDW1(JIJ),PTSTEP,PDZZ(JIJ,JK),PRHODREF(JIJ,JK),PRXT(JIJ),ZINVTSTEP) &
+ & + FWSED2(ZWSEDW2(JIJ),PTSTEP,PDZZ(JIJ,JK),ZSED(JIJ,IKPLUS,JRR))
+ ELSE
+ ZSED(JIJ,IK,JRR)=FWSED1(ZWSEDW1(JIJ),PTSTEP,PDZZ(JIJ,JK),PRHODREF(JIJ,JK),PRXT(JIJ),ZINVTSTEP)
+ ENDIF
!-------------------------------------------------------------------------------------------
- ENDDO
ENDDO
!!IF (LHOOK) CALL DR_HOOK('ICE4_SEDIMENTATION_STAT:SNOW',1,ZHOOK_HANDLE)
@@ -420,41 +416,39 @@ CONTAINS
REAL, INTENT(IN) :: PFSED
REAL, INTENT(IN) :: PEXSED
- REAL, INTENT(IN) :: PRXT(D%NIT,D%NJT) ! mr of specy X
+ REAL, INTENT(IN) :: PRXT(D%NIJT) ! mr of specy X
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!!IF (LHOOK) CALL DR_HOOK('ICE4_SEDIMENTATION_STAT:OTHER_SPECIES',0,ZHOOK_HANDLE)
! for all but cloud and pristine ice :
- DO JJ = IJB, IJE
- DO JI = IIB, IIE
- ZQP=ZSED(JI,JJ,IKPLUS,JRR)*ZTSORHODZ(JI,JJ)
- IF ((PRXT(JI,JJ) > ICED%XRTMIN(JRR)) .OR. (ZQP > ICED%XRTMIN(JRR))) THEN
- !calculation of w
- IF ( PRXT(JI,JJ) > ICED%XRTMIN(JRR) ) THEN
- ZWSEDW1(JI)= PFSED *PRXT(JI,JJ)**(PEXSED-1)*PRHODREF(JI,JJ,JK)**(PEXSED-ICED%XCEXVT-1)
- ELSE
- ZWSEDW1(JI)=0.
- ENDIF
- IF ( ZQP > ICED%XRTMIN(JRR) ) THEN
- ZWSEDW2(JI)= PFSED *ZQP**(PEXSED-1)*PRHODREF(JI,JJ,JK)**(PEXSED-ICED%XCEXVT-1)
- ELSE
- ZWSEDW2(JI)=0.
- ENDIF
+ DO JIJ = IIJB, IIJE
+ ZQP=ZSED(JIJ,IKPLUS,JRR)*ZTSORHODZ(JIJ)
+ IF ((PRXT(JIJ) > ICED%XRTMIN(JRR)) .OR. (ZQP > ICED%XRTMIN(JRR))) THEN
+ !calculation of w
+ IF ( PRXT(JIJ) > ICED%XRTMIN(JRR) ) THEN
+ ZWSEDW1(JIJ)= PFSED *PRXT(JIJ)**(PEXSED-1)*PRHODREF(JIJ,JK)**(PEXSED-ICED%XCEXVT-1)
ELSE
- ZWSEDW1(JI)=0.
- ZWSEDW2(JI)=0.
+ ZWSEDW1(JIJ)=0.
ENDIF
-!- duplicated code -------------------------------------------------------------------------
- IF (ZWSEDW2(JI) /= 0.) THEN
- ZSED(JI,JJ,IK,JRR)=FWSED1(ZWSEDW1(JI),PTSTEP,PDZZ(JI,JJ,JK),PRHODREF(JI,JJ,JK),PRXT(JI,JJ),ZINVTSTEP) &
- & + FWSED2(ZWSEDW2(JI),PTSTEP,PDZZ(JI,JJ,JK),ZSED(JI,JJ,IKPLUS,JRR))
+ IF ( ZQP > ICED%XRTMIN(JRR) ) THEN
+ ZWSEDW2(JIJ)= PFSED *ZQP**(PEXSED-1)*PRHODREF(JIJ,JK)**(PEXSED-ICED%XCEXVT-1)
ELSE
- ZSED(JI,JJ,IK,JRR)=FWSED1(ZWSEDW1(JI),PTSTEP,PDZZ(JI,JJ,JK),PRHODREF(JI,JJ,JK),PRXT(JI,JJ),ZINVTSTEP)
+ ZWSEDW2(JIJ)=0.
ENDIF
+ ELSE
+ ZWSEDW1(JIJ)=0.
+ ZWSEDW2(JIJ)=0.
+ ENDIF
+!- duplicated code -------------------------------------------------------------------------
+ IF (ZWSEDW2(JIJ) /= 0.) THEN
+ ZSED(JIJ,IK,JRR)=FWSED1(ZWSEDW1(JIJ),PTSTEP,PDZZ(JIJ,JK),PRHODREF(JIJ,JK),PRXT(JIJ),ZINVTSTEP) &
+ & + FWSED2(ZWSEDW2(JIJ),PTSTEP,PDZZ(JIJ,JK),ZSED(JIJ,IKPLUS,JRR))
+ ELSE
+ ZSED(JIJ,IK,JRR)=FWSED1(ZWSEDW1(JIJ),PTSTEP,PDZZ(JIJ,JK),PRHODREF(JIJ,JK),PRXT(JIJ),ZINVTSTEP)
+ ENDIF
!-------------------------------------------------------------------------------------------
- ENDDO
ENDDO
!!IF (LHOOK) CALL DR_HOOK('ICE4_SEDIMENTATION_STAT:OTHER_SPECIES',1,ZHOOK_HANDLE)
diff --git a/src/common/micro/rain_ice.F90 b/src/common/micro/rain_ice.F90
index aa9b05d272540471381ce91dd65f468fb73b6be2..44edc4c0ee65413e7c9f5434da0c9080f4864b2d 100644
--- a/src/common/micro/rain_ice.F90
+++ b/src/common/micro/rain_ice.F90
@@ -202,11 +202,10 @@ USE MODE_BUDGET, ONLY: BUDGET_STORE_ADD_PHY, BUDGET_STORE_INIT_PHY, BUDG
USE MODE_MSG, ONLY: PRINT_MSG, NVERB_FATAL
USE MODE_ICE4_RAINFR_VERT, ONLY: ICE4_RAINFR_VERT
-USE MODE_ICE4_SEDIMENTATION_STAT, ONLY: ICE4_SEDIMENTATION_STAT
-USE MODE_ICE4_SEDIMENTATION_SPLIT, ONLY: ICE4_SEDIMENTATION_SPLIT
-USE MODE_ICE4_SEDIMENTATION_SPLIT_MOMENTUM, ONLY: ICE4_SEDIMENTATION_SPLIT_MOMENTUM
+USE MODE_ICE4_SEDIMENTATION, ONLY: ICE4_SEDIMENTATION
USE MODE_ICE4_TENDENCIES, ONLY: ICE4_TENDENCIES
USE MODE_ICE4_NUCLEATION, ONLY: ICE4_NUCLEATION
+USE MODE_ICE4_CORRECT_NEGATIVITIES, ONLY: ICE4_CORRECT_NEGATIVITIES
!
IMPLICIT NONE
!
@@ -302,8 +301,6 @@ 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
!
LOGICAL :: GEXT_TEND
LOGICAL :: LSOFT ! Must we really compute tendencies or only adjust them to new T variables
@@ -504,125 +501,25 @@ END DO
! -------------------------------------
!
IF(.NOT. PARAMI%LSEDIM_AFTER) THEN
- !
- !* 2.1 sedimentation
- !
- IF (BUCONF%LBUDGET_RC .AND. PARAMI%LSEDIC) CALL BUDGET_STORE_INIT_PHY(D, TBUDGETS(NBUDGET_RC), 'SEDI', PRCS(:, :) * PRHODJ(:, :))
- IF (BUCONF%LBUDGET_RR) CALL BUDGET_STORE_INIT_PHY(D, TBUDGETS(NBUDGET_RR), 'SEDI', PRRS(:, :) * PRHODJ(:, :))
- IF (BUCONF%LBUDGET_RI) CALL BUDGET_STORE_INIT_PHY(D, TBUDGETS(NBUDGET_RI), 'SEDI', PRIS(:, :) * PRHODJ(:, :))
- IF (BUCONF%LBUDGET_RS) CALL BUDGET_STORE_INIT_PHY(D, TBUDGETS(NBUDGET_RS), 'SEDI', PRSS(:, :) * PRHODJ(:, :))
- IF (BUCONF%LBUDGET_RG) CALL BUDGET_STORE_INIT_PHY(D, TBUDGETS(NBUDGET_RG), 'SEDI', PRGS(:, :) * PRHODJ(:, :))
- IF (BUCONF%LBUDGET_RH .AND. KRR==7) CALL BUDGET_STORE_INIT_PHY(D, TBUDGETS(NBUDGET_RH), 'SEDI', PRHS(:, :) * PRHODJ(:, :))
-
- IF(PARAMI%CSEDIM=='STAT') THEN
- IF(KRR==7) THEN
- DO JK = IKTB,IKTE
- DO JIJ = IIJB,IIJE
- ZRCT(JIJ,JK)=PRCS(JIJ,JK)*PTSTEP
- ZRRT(JIJ,JK)=PRRS(JIJ,JK)*PTSTEP
- ZRIT(JIJ,JK)=PRIS(JIJ,JK)*PTSTEP
- ZRST(JIJ,JK)=PRSS(JIJ,JK)*PTSTEP
- ZRGT(JIJ,JK)=PRGS(JIJ,JK)*PTSTEP
- ZRHT(JIJ,JK)=PRHS(JIJ,JK)*PTSTEP
- ENDDO
- ENDDO
- CALL ICE4_SEDIMENTATION_STAT(D, CST, ICEP, ICED, &
- &PTSTEP, KRR, PARAMI%LSEDIC, PDZZ, &
- &PRHODREF, PPABST, PTHT, PRHODJ, &
- &ZLBDAS, &
- &PRCS, ZRCT, PRRS, ZRRT, PRIS, ZRIT,&
- &PRSS, ZRST, PRGS, ZRGT,&
- &PINPRC, PINPRR, ZINPRI, PINPRS, PINPRG, &
- &PSEA=PSEA, PTOWN=PTOWN, &
- &PINPRH=PINPRH, PRHT=ZRHT, PRHS=PRHS, PFPR=PFPR)
- ELSE
- DO JK = IKTB,IKTE
- DO JIJ = IIJB,IIJE
- ZRCT(JIJ,JK)=PRCS(JIJ,JK)*PTSTEP
- ZRRT(JIJ,JK)=PRRS(JIJ,JK)*PTSTEP
- ZRIT(JIJ,JK)=PRIS(JIJ,JK)*PTSTEP
- ZRST(JIJ,JK)=PRSS(JIJ,JK)*PTSTEP
- ZRGT(JIJ,JK)=PRGS(JIJ,JK)*PTSTEP
- ENDDO
- ENDDO
- CALL ICE4_SEDIMENTATION_STAT(D, CST, ICEP, ICED, &
- &PTSTEP, KRR, PARAMI%LSEDIC, PDZZ, &
- &PRHODREF, PPABST, PTHT, PRHODJ, &
- &ZLBDAS, &
- &PRCS, ZRCT, PRRS, ZRRT, PRIS, ZRIT,&
- &PRSS, ZRST, PRGS, ZRGT,&
- &PINPRC, PINPRR, ZINPRI, PINPRS, PINPRG, &
- &PSEA=PSEA, PTOWN=PTOWN, &
- &PFPR=PFPR)
- ENDIF
- PINPRS(IIJB:IIJE) = PINPRS(IIJB:IIJE) + ZINPRI(IIJB:IIJE)
- !No negativity correction here as we apply sedimentation on PR.S*PTSTEP variables
- ELSEIF(PARAMI%CSEDIM=='SPLI') THEN
- IF(KRR==7) THEN
- CALL ICE4_SEDIMENTATION_SPLIT(D, CST, ICEP, ICED, PARAMI, &
- &PTSTEP, KRR, PDZZ, &
- &PRHODREF, PPABST, PTHT, ZT, PRHODJ, &
- &PRCS, PRCT, PRRS, PRRT, PRIS, PRIT, PRSS, PRST, PRGS, PRGT,&
- &PINPRC, PINPRR, ZINPRI, PINPRS, PINPRG, &
- &PSEA=PSEA, PTOWN=PTOWN, &
- &PINPRH=PINPRH, PRHT=PRHT, PRHS=PRHS, PFPR=PFPR)
- ELSE
- CALL ICE4_SEDIMENTATION_SPLIT(D, CST, ICEP, ICED, PARAMI, &
- &PTSTEP, KRR, PDZZ, &
- &PRHODREF, PPABST, PTHT, ZT, PRHODJ, &
- &PRCS, PRCT, PRRS, PRRT, PRIS, PRIT, PRSS, PRST, PRGS, PRGT,&
- &PINPRC, PINPRR, ZINPRI, PINPRS, PINPRG, &
- &PSEA=PSEA, PTOWN=PTOWN, &
- &PFPR=PFPR)
- ENDIF
- PINPRS(IIJB:IIJE) = PINPRS(IIJB:IIJE) + ZINPRI(IIJB:IIJE)
- !We correct negativities with conservation
- !SPLI algorith uses a time-splitting. Inside the loop a temporary m.r. is used.
- ! It is initialized with the m.r. at T and is modified by two tendencies:
- ! sedimentation tendency and an external tendency which represents all other
- ! processes (mainly advection and microphysical processes). If both tendencies
- ! are negative, sedimentation can remove a species at a given sub-timestep. From
- ! this point sedimentation stops for the remaining sub-timesteps but the other tendency
- ! will be still active and will lead to negative values.
- ! We could prevent the algorithm to not consume too much a species, instead we apply
- ! a correction here.
- CALL CORRECT_NEGATIVITIES(D, KRR, PRVS, PRCS, PRRS, &
- &PRIS, PRSS, PRGS, &
- &PTHS, ZZ_LVFACT, ZZ_LSFACT, PRHS)
- ELSEIF(PARAMI%CSEDIM=='NONE') THEN
+ IF(KRR==7) THEN
+ CALL ICE4_SEDIMENTATION(D, CST, ICEP, ICED, PARAMI, BUCONF, &
+ &PTSTEP, KRR, PDZZ, &
+ &ZZ_LVFACT, ZZ_LSFACT, PRHODREF, PPABST, PTHT, ZT, PRHODJ, &
+ &PTHS, PRVS, PRCS, PRCT, PRRS, PRRT, PRIS, PRIT, PRSS, PRST, PRGS, PRGT,&
+ &PINPRC, PINPRR, PINPRS, PINPRG, &
+ &TBUDGETS, KBUDGETS, &
+ &PSEA=PSEA, PTOWN=PTOWN, &
+ &PINPRH=PINPRH, PRHT=PRHT, PRHS=PRHS, PFPR=PFPR)
ELSE
- CALL PRINT_MSG(NVERB_FATAL, 'GEN', 'RAIN_ICE', 'no sedimentation scheme for PARAMI%CSEDIM='//PARAMI%CSEDIM)
- END IF
-
-
-
-
-
-
-!!!!! ajouter momentum
-
-
-
-
-
-
-
-
-
-
-
-
-
-
- !
- !* 2.2 budget storage
- !
- IF (BUCONF%LBUDGET_RC .AND. PARAMI%LSEDIC) CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_RC), 'SEDI', PRCS(:, :) * PRHODJ(:, :))
- IF (BUCONF%LBUDGET_RR) CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_RR), 'SEDI', PRRS(:, :) * PRHODJ(:, :))
- IF (BUCONF%LBUDGET_RI) CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_RI), 'SEDI', PRIS(:, :) * PRHODJ(:, :))
- IF (BUCONF%LBUDGET_RS) CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_RS), 'SEDI', PRSS(:, :) * PRHODJ(:, :))
- IF (BUCONF%LBUDGET_RG) CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_RG), 'SEDI', PRGS(:, :) * PRHODJ(:, :))
- IF (BUCONF%LBUDGET_RH .AND. KRR==7) CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_RH), 'SEDI', PRHS(:, :) * PRHODJ(:, :))
+ CALL ICE4_SEDIMENTATION(D, CST, ICEP, ICED, PARAMI, BUCONF, &
+ &PTSTEP, KRR, PDZZ, &
+ &ZZ_LVFACT, ZZ_LSFACT, PRHODREF, PPABST, PTHT, ZT, PRHODJ, &
+ &PTHS, PRVS, PRCS, PRCT, PRRS, PRRT, PRIS, PRIT, PRSS, PRST, PRGS, PRGT,&
+ &PINPRC, PINPRR, PINPRS, PINPRG, &
+ &TBUDGETS, KBUDGETS, &
+ &PSEA=PSEA, PTOWN=PTOWN, &
+ &PFPR=PFPR)
+ ENDIF
ENDIF
!
@@ -1571,9 +1468,9 @@ END IF
! call must only be due to the correction of negativities.
!
!We correct negativities with conservation
-CALL CORRECT_NEGATIVITIES(D, KRR, PRVS, PRCS, PRRS, &
- &PRIS, PRSS, PRGS, &
- &PTHS, ZZ_LVFACT, ZZ_LSFACT, PRHS)
+CALL ICE4_CORRECT_NEGATIVITIES(D, ICED, KRR, PRVS, PRCS, PRRS, &
+ &PRIS, PRSS, PRGS, &
+ &PTHS, ZZ_LVFACT, ZZ_LSFACT, PRHS)
IF (BUCONF%LBU_ENABLE) THEN
IF (BUCONF%LBUDGET_TH) CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_TH), 'CORR', PTHS(:, :)*PRHODJ(:, :))
@@ -1592,104 +1489,25 @@ END IF
! -------------------------------------
!
IF(PARAMI%LSEDIM_AFTER) THEN
- !
- !* 8.1 sedimentation
- !
- IF (BUCONF%LBUDGET_RC .AND. PARAMI%LSEDIC) CALL BUDGET_STORE_INIT_PHY(D, TBUDGETS(NBUDGET_RC), 'SEDI', PRCS(:, :) * PRHODJ(:, :))
- IF (BUCONF%LBUDGET_RR) CALL BUDGET_STORE_INIT_PHY(D, TBUDGETS(NBUDGET_RR), 'SEDI', PRRS(:, :) * PRHODJ(:, :))
- IF (BUCONF%LBUDGET_RI) CALL BUDGET_STORE_INIT_PHY(D, TBUDGETS(NBUDGET_RI), 'SEDI', PRIS(:, :) * PRHODJ(:, :))
- IF (BUCONF%LBUDGET_RS) CALL BUDGET_STORE_INIT_PHY(D, TBUDGETS(NBUDGET_RS), 'SEDI', PRSS(:, :) * PRHODJ(:, :))
- IF (BUCONF%LBUDGET_RG) CALL BUDGET_STORE_INIT_PHY(D, TBUDGETS(NBUDGET_RG), 'SEDI', PRGS(:, :) * PRHODJ(:, :))
- IF (BUCONF%LBUDGET_RH .AND. KRR==7) CALL BUDGET_STORE_INIT_PHY(D, TBUDGETS(NBUDGET_RH), 'SEDI', PRHS(:, :) * PRHODJ(:, :))
-
- IF(PARAMI%CSEDIM=='STAT') THEN
- IF (KRR==7) THEN
- DO JK = IKTB,IKTE
- DO JIJ = IIJB,IIJE
- ZRCT(JIJ,JK)=PRCS(JIJ,JK)*PTSTEP
- ZRRT(JIJ,JK)=PRRS(JIJ,JK)*PTSTEP
- ZRIT(JIJ,JK)=PRIS(JIJ,JK)*PTSTEP
- ZRST(JIJ,JK)=PRSS(JIJ,JK)*PTSTEP
- ZRGT(JIJ,JK)=PRGS(JIJ,JK)*PTSTEP
- ZRHT(JIJ,JK)=PRHS(JIJ,JK)*PTSTEP
- ENDDO
- ENDDO
- CALL ICE4_SEDIMENTATION_STAT(D, CST, ICEP, ICED, &
- &PTSTEP, KRR, PARAMI%LSEDIC, PDZZ, &
- &PRHODREF, PPABST, PTHT, PRHODJ, &
- &ZLBDAS, &
- &PRCS, ZRCT, PRRS, ZRRT, PRIS, ZRIT,&
- &PRSS, ZRST, PRGS, ZRGT,&
- &PINPRC, PINPRR, ZINPRI, PINPRS, PINPRG, &
- &PSEA=PSEA, PTOWN=PTOWN, &
- &PINPRH=PINPRH, PRHT=ZRHT, PRHS=PRHS, PFPR=PFPR)
- ELSE
- DO JK = IKTB,IKTE
- DO JIJ = IIJB,IIJE
- ZRCT(JIJ,JK)=PRCS(JIJ,JK)*PTSTEP
- ZRRT(JIJ,JK)=PRRS(JIJ,JK)*PTSTEP
- ZRIT(JIJ,JK)=PRIS(JIJ,JK)*PTSTEP
- ZRST(JIJ,JK)=PRSS(JIJ,JK)*PTSTEP
- ZRGT(JIJ,JK)=PRGS(JIJ,JK)*PTSTEP
- ENDDO
- ENDDO
- CALL ICE4_SEDIMENTATION_STAT(D, CST, ICEP, ICED, &
- &PTSTEP, KRR, PARAMI%LSEDIC, PDZZ, &
- &PRHODREF, PPABST, PTHT, PRHODJ, &
- &ZLBDAS, &
- &PRCS, ZRCT, PRRS, ZRRT, PRIS, ZRIT,&
- &PRSS, ZRST, PRGS, ZRGT,&
- &PINPRC, PINPRR, ZINPRI, PINPRS, PINPRG, &
- &PSEA=PSEA, PTOWN=PTOWN, &
- &PFPR=PFPR)
- ENDIF
- PINPRS(IIJB:IIJE) = PINPRS(IIJB:IIJE) + ZINPRI(IIJB:IIJE)
- !No negativity correction here as we apply sedimentation on PR.S*PTSTEP variables
- ELSEIF(PARAMI%CSEDIM=='SPLI') THEN
- !SR: It *seems* that we must have two separate calls for ifort
- IF(KRR==7) THEN
- CALL ICE4_SEDIMENTATION_SPLIT(D, CST, ICEP, ICED, PARAMI, &
- &PTSTEP, KRR, PDZZ, &
- &PRHODREF, PPABST, PTHT, ZT, PRHODJ, &
- &PRCS, PRCT, PRRS, PRRT, PRIS, PRIT, PRSS, PRST, PRGS, PRGT,&
- &PINPRC, PINPRR, ZINPRI, PINPRS, PINPRG, &
- &PSEA=PSEA, PTOWN=PTOWN, &
- &PINPRH=PINPRH, PRHT=PRHT, PRHS=PRHS, PFPR=PFPR)
- ELSE
- CALL ICE4_SEDIMENTATION_SPLIT(D, CST, ICEP, ICED, PARAMI, &
- &PTSTEP, KRR, PDZZ, &
- &PRHODREF, PPABST, PTHT, ZT, PRHODJ, &
- &PRCS, PRCT, PRRS, PRRT, PRIS, PRIT, PRSS, PRST, PRGS, PRGT,&
- &PINPRC, PINPRR, ZINPRI, PINPRS, PINPRG, &
- &PSEA=PSEA, PTOWN=PTOWN, &
- &PFPR=PFPR)
- ENDIF
- PINPRS(IIJB:IIJE) = PINPRS(IIJB:IIJE) + ZINPRI(IIJB:IIJE)
- !We correct negativities with conservation
- !SPLI algorith uses a time-splitting. Inside the loop a temporary m.r. is used.
- ! It is initialized with the m.r. at T and is modified by two tendencies:
- ! sedimentation tendency and an external tendency which represents all other
- ! processes (mainly advection and microphysical processes). If both tendencies
- ! are negative, sedimentation can remove a species at a given sub-timestep. From
- ! this point sedimentation stops for the remaining sub-timesteps but the other tendency
- ! will be still active and will lead to negative values.
- ! We could prevent the algorithm to not consume too much a species, instead we apply
- ! a correction here.
- CALL CORRECT_NEGATIVITIES(D, KRR, PRVS, PRCS, PRRS, &
- &PRIS, PRSS, PRGS, &
- &PTHS, ZZ_LVFACT, ZZ_LSFACT, PRHS)
+ IF(KRR==7) THEN
+ CALL ICE4_SEDIMENTATION(D, CST, ICEP, ICED, PARAMI, BUCONF, &
+ &PTSTEP, KRR, PDZZ, &
+ &ZZ_LVFACT, ZZ_LSFACT, PRHODREF, PPABST, PTHT, ZT, PRHODJ, &
+ &PTHS, PRVS, PRCS, PRCT, PRRS, PRRT, PRIS, PRIT, PRSS, PRST, PRGS, PRGT,&
+ &PINPRC, PINPRR, PINPRS, PINPRG, &
+ &TBUDGETS, KBUDGETS, &
+ &PSEA=PSEA, PTOWN=PTOWN, &
+ &PINPRH=PINPRH, PRHT=PRHT, PRHS=PRHS, PFPR=PFPR)
ELSE
- CALL PRINT_MSG(NVERB_FATAL, 'GEN', 'RAIN_ICE', 'no sedimentation scheme for PARAMI%CSEDIM='//PARAMI%CSEDIM)
- END IF
- !
- !* 8.2 budget storage
- !
- IF (BUCONF%LBUDGET_RC .AND. PARAMI%LSEDIC) CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_RC), 'SEDI', PRCS(:, :) * PRHODJ(:, :))
- IF (BUCONF%LBUDGET_RR) CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_RR), 'SEDI', PRRS(:, :) * PRHODJ(:, :))
- IF (BUCONF%LBUDGET_RI) CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_RI), 'SEDI', PRIS(:, :) * PRHODJ(:, :))
- IF (BUCONF%LBUDGET_RS) CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_RS), 'SEDI', PRSS(:, :) * PRHODJ(:, :))
- IF (BUCONF%LBUDGET_RG) CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_RG), 'SEDI', PRGS(:, :) * PRHODJ(:, :))
- IF (BUCONF%LBUDGET_RH .AND. KRR==7) CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_RH), 'SEDI', PRHS(:, :) * PRHODJ(:, :))
+ CALL ICE4_SEDIMENTATION(D, CST, ICEP, ICED, PARAMI, BUCONF, &
+ &PTSTEP, KRR, PDZZ, &
+ &ZZ_LVFACT, ZZ_LSFACT, PRHODREF, PPABST, PTHT, ZT, PRHODJ, &
+ &PTHS, PRVS, PRCS, PRCT, PRRS, PRRT, PRIS, PRIT, PRSS, PRST, PRGS, PRGT,&
+ &PINPRC, PINPRR, PINPRS, PINPRG, &
+ &TBUDGETS, KBUDGETS, &
+ &PSEA=PSEA, PTOWN=PTOWN, &
+ &PFPR=PFPR)
+ ENDIF
!"sedimentation" of rain fraction
IF (PRESENT(PRHS)) THEN
@@ -1724,103 +1542,4 @@ ENDIF
IF (LHOOK) CALL DR_HOOK('RAIN_ICE', 1, ZHOOK_HANDLE)
!
-CONTAINS
- !
- SUBROUTINE CORRECT_NEGATIVITIES(D, KRR, PRV, PRC, PRR, &
- &PRI, PRS, PRG, &
- &PTH, PLVFACT, PLSFACT, PRH)
- !
- IMPLICIT NONE
- !
- TYPE(DIMPHYEX_t), INTENT(IN) :: D
- INTEGER, INTENT(IN) :: KRR
- REAL, DIMENSION(D%NIJT, D%NKT), INTENT(INOUT) :: PRV, PRC, PRR, PRI, PRS, PRG, PTH
- REAL, DIMENSION(D%NIJT, D%NKT), INTENT(IN) :: PLVFACT, PLSFACT
- REAL, DIMENSION(D%NIJT, D%NKT), OPTIONAL, INTENT(INOUT) :: PRH
- !
- REAL :: ZW
- INTEGER :: JIJ, JK
- REAL(KIND=JPRB) :: ZHOOK_HANDLE
- !
- IF (LHOOK) CALL DR_HOOK('RAIN_ICE:CORRECT_NEGATIVITIES', 0, ZHOOK_HANDLE)
- !
- !We correct negativities with conservation
- DO JK = IKTB, IKTE
- DO JIJ = IIJB, IIJE
- ! 1) deal with negative values for mixing ratio, except for vapor
- ZW =PRC(JIJ,JK)-MAX(PRC(JIJ,JK), 0.)
- PRV(JIJ,JK)=PRV(JIJ,JK)+ZW
- PTH(JIJ,JK)=PTH(JIJ,JK)-ZW*PLVFACT(JIJ,JK)
- PRC(JIJ,JK)=PRC(JIJ,JK)-ZW
-
- ZW =PRR(JIJ,JK)-MAX(PRR(JIJ,JK), 0.)
- PRV(JIJ,JK)=PRV(JIJ,JK)+ZW
- PTH(JIJ,JK)=PTH(JIJ,JK)-ZW*PLVFACT(JIJ,JK)
- PRR(JIJ,JK)=PRR(JIJ,JK)-ZW
-
- ZW =PRI(JIJ,JK)-MAX(PRI(JIJ,JK), 0.)
- PRV(JIJ,JK)=PRV(JIJ,JK)+ZW
- PTH(JIJ,JK)=PTH(JIJ,JK)-ZW*PLSFACT(JIJ,JK)
- PRI(JIJ,JK)=PRI(JIJ,JK)-ZW
-
- ZW =PRS(JIJ,JK)-MAX(PRS(JIJ,JK), 0.)
- PRV(JIJ,JK)=PRV(JIJ,JK)+ZW
- PTH(JIJ,JK)=PTH(JIJ,JK)-ZW*PLSFACT(JIJ,JK)
- PRS(JIJ,JK)=PRS(JIJ,JK)-ZW
-
- ZW =PRG(JIJ,JK)-MAX(PRG(JIJ,JK), 0.)
- PRV(JIJ,JK)=PRV(JIJ,JK)+ZW
- PTH(JIJ,JK)=PTH(JIJ,JK)-ZW*PLSFACT(JIJ,JK)
- PRG(JIJ,JK)=PRG(JIJ,JK)-ZW
-
- IF(KRR==7) THEN
- ZW =PRH(JIJ,JK)-MAX(PRH(JIJ,JK), 0.)
- PRV(JIJ,JK)=PRV(JIJ,JK)+ZW
- PTH(JIJ,JK)=PTH(JIJ,JK)-ZW*PLSFACT(JIJ,JK)
- PRH(JIJ,JK)=PRH(JIJ,JK)-ZW
- ENDIF
-
- ! 2) deal with negative vapor mixing ratio
-
- ! for rc and ri, we keep ice fraction constant
- ZW=MIN(1., MAX(ICED%XRTMIN(1)-PRV(JIJ,JK), 0.) / &
- &MAX(PRC(JIJ,JK)+PRI(JIJ,JK), 1.E-20)) ! Proportion of rc+ri to convert into rv
- PTH(JIJ,JK)=PTH(JIJ,JK)-ZW* &
- &(PRC(JIJ,JK)*PLVFACT(JIJ,JK)+PRI(JIJ,JK)*PLSFACT(JIJ,JK))
- PRV(JIJ,JK)=PRV(JIJ,JK)+ZW*(PRC(JIJ,JK)+PRI(JIJ,JK))
- PRC(JIJ,JK)=(1.-ZW)*PRC(JIJ,JK)
- PRI(JIJ,JK)=(1.-ZW)*PRI(JIJ,JK)
-
- ZW=MIN(MAX(PRR(JIJ,JK), 0.), &
- &MAX(ICED%XRTMIN(1)-PRV(JIJ,JK), 0.)) ! Quantity of rr to convert into rv
- PRV(JIJ,JK)=PRV(JIJ,JK)+ZW
- PRR(JIJ,JK)=PRR(JIJ,JK)-ZW
- PTH(JIJ,JK)=PTH(JIJ,JK)-ZW*PLVFACT(JIJ,JK)
-
- ZW=MIN(MAX(PRS(JIJ,JK), 0.), &
- &MAX(ICED%XRTMIN(1)-PRV(JIJ,JK), 0.)) ! Quantity of rs to convert into rv
- PRV(JIJ,JK)=PRV(JIJ,JK)+ZW
- PRS(JIJ,JK)=PRS(JIJ,JK)-ZW
- PTH(JIJ,JK)=PTH(JIJ,JK)-ZW*PLSFACT(JIJ,JK)
-
- ZW=MIN(MAX(PRG(JIJ,JK), 0.), &
- &MAX(ICED%XRTMIN(1)-PRV(JIJ,JK), 0.)) ! Quantity of rg to convert into rv
- PRV(JIJ,JK)=PRV(JIJ,JK)+ZW
- PRG(JIJ,JK)=PRG(JIJ,JK)-ZW
- PTH(JIJ,JK)=PTH(JIJ,JK)-ZW*PLSFACT(JIJ,JK)
-
- IF(KRR==7) THEN
- ZW=MIN(MAX(PRH(JIJ,JK), 0.), &
- &MAX(ICED%XRTMIN(1)-PRV(JIJ,JK), 0.)) ! Quantity of rh to convert into rv
- PRV(JIJ,JK)=PRV(JIJ,JK)+ZW
- PRH(JIJ,JK)=PRH(JIJ,JK)-ZW
- PTH(JIJ,JK)=PTH(JIJ,JK)-ZW*PLSFACT(JIJ,JK)
- ENDIF
- ENDDO
- ENDDO
- !
- IF (LHOOK) CALL DR_HOOK('RAIN_ICE:CORRECT_NEGATIVITIES', 1, ZHOOK_HANDLE)
- !
- END SUBROUTINE CORRECT_NEGATIVITIES
-!
END SUBROUTINE RAIN_ICE