From b4983148019eddd3ea106cd61e54960ef4a0ed85 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Beno=C3=AEt=20Vi=C3=A9?= <benoit.vie@meteo.fr>
Date: Wed, 18 May 2022 13:34:37 +0200
Subject: [PATCH] bugfix aer2lima (call only if ORILAM or DUST or SALT) &
rain_accr_snow (rhodref)
---
src/MNH/aer2lima.f90 | 73 ++++++++++---------
src/MNH/ini_lima_cold_mixed.f90 | 2 +-
src/MNH/lima_ccn_activation.f90 | 4 +-
src/MNH/lima_init_ccn_activation_spectrum.f90 | 4 +-
src/MNH/lima_mixed_fast_processes.f90 | 8 +-
src/MNH/lima_rain_accr_snow.f90 | 19 +++--
src/MNH/lima_warm_nucl.f90 | 4 +-
src/MNH/resolved_cloud.f90 | 4 +-
8 files changed, 60 insertions(+), 58 deletions(-)
diff --git a/src/MNH/aer2lima.f90 b/src/MNH/aer2lima.f90
index 885dc0dc0..eaf03fa5b 100644
--- a/src/MNH/aer2lima.f90
+++ b/src/MNH/aer2lima.f90
@@ -278,29 +278,31 @@ END DO
ELSE ! keep lima class intiatialization
IF (CACTCCN=="ABRK") THEN
! only one CCN_FREE mode (activation is not performed upon aerosol class but by physical paramters)
- IF (NMOD_CCN .GE. 2) &
- ZCCN_SUM(:,:,:,1) = ZCCN_SUM(:,:,:,1) + &
- PSVT(:,:,:,NSV_LIMA_CCN_FREE+1) + PSVT(:,:,:,NSV_LIMA_CCN_ACTI+1)
- IF (NMOD_CCN .GE. 3) &
- ZCCN_SUM(:,:,:,1) = ZCCN_SUM(:,:,:,1) + &
- PSVT(:,:,:,NSV_LIMA_CCN_FREE+2) + PSVT(:,:,:,NSV_LIMA_CCN_ACTI+2)
-
+ IF (NMOD_CCN .GE. 2) THEN
+ DO JI = 2, NMOD_CCN
+ ZCCN_SUM(:,:,:,1) = ZCCN_SUM(:,:,:,1) + &
+ PSVT(:,:,:,NSV_LIMA_CCN_FREE+JI-1) + PSVT(:,:,:,NSV_LIMA_CCN_ACTI+JI-1)
+ END DO
+ END IF
ELSE
- IF (NMOD_CCN .GE. 2) &
- ZCCN_SUM(:,:,:,2) = PSVT(:,:,:,NSV_LIMA_CCN_FREE+1) + PSVT(:,:,:,NSV_LIMA_CCN_ACTI+1)
-
- IF (NMOD_CCN .GE. 3) &
- ZCCN_SUM(:,:,:,3) = PSVT(:,:,:,NSV_LIMA_CCN_FREE+2) + PSVT(:,:,:,NSV_LIMA_CCN_ACTI+2)
+ IF (NMOD_CCN .GE. 2) THEN
+ DO JI = 2, NMOD_CCN
+ ZCCN_SUM(:,:,:,JI) = PSVT(:,:,:,NSV_LIMA_CCN_FREE+JI-1) + PSVT(:,:,:,NSV_LIMA_CCN_ACTI+JI-1)
+ END DO
+ END IF
END IF
- IF (.NOT.(LDUST)) &
- ZIFN_SUM(:,:,:,1) = PSVT(:,:,:,NSV_LIMA_IFN_FREE) + PSVT(:,:,:,NSV_LIMA_IFN_NUCL)
-
- IF (NMOD_IFN .GE. 2) &
- ZIFN_SUM(:,:,:,2) = PSVT(:,:,:,NSV_LIMA_IFN_FREE+1) + PSVT(:,:,:,NSV_LIMA_IFN_NUCL+1)
+ IF (.NOT.(LDUST) .AND. NMOD_IFN.GE.1) &
+ ZIFN_SUM(:,:,:,1) = PSVT(:,:,:,NSV_LIMA_IFN_FREE) + PSVT(:,:,:,NSV_LIMA_IFN_NUCL)
+ IF (NMOD_IFN .GE. 2) THEN
+ DO JI = 2, NMOD_IFN
+ ZIFN_SUM(:,:,:,JI) = PSVT(:,:,:,NSV_LIMA_IFN_FREE+JI-1) + PSVT(:,:,:,NSV_LIMA_IFN_NUCL+JI-1)
+ END DO
+ END IF
END IF ! end if sur LORILAM
-
+!
+!
! Sea Salt part
IF (LSALT) THEN
!
@@ -329,7 +331,8 @@ ELSE ! keep lima class intiatialization for sea salt + ccn from orilam
ZCCN_SUM(:,:,:,1) = PSVT(:,:,:,NSV_LIMA_CCN_FREE) + PSVT(:,:,:,NSV_LIMA_CCN_ACTI)
END IF ! end if sur LSALT
-
+!
+!
! Dust part
IF (LDUST) THEN
! initatialization of dust if not macc
@@ -352,24 +355,24 @@ IF (LDUST) THEN
END DO
ELSE ! keep lima class intiatialization
-
- ZIFN_SUM(:,:,:,1) = PSVT(:,:,:,NSV_LIMA_IFN_FREE) + PSVT(:,:,:,NSV_LIMA_IFN_NUCL)
+ IF (NMOD_IFN.GE.1) &
+ ZIFN_SUM(:,:,:,1) = PSVT(:,:,:,NSV_LIMA_IFN_FREE) + PSVT(:,:,:,NSV_LIMA_IFN_NUCL)
END IF ! endif sur LDUST
-
-PSVT(:,:,:,NSV_LIMA_CCN_FREE) = MAX(ZCCN_SUM(:,:,:,1) - PSVT(:,:,:,NSV_LIMA_CCN_ACTI), 0.)
-
-IF (NMOD_CCN .GE. 2) &
-PSVT(:,:,:,NSV_LIMA_CCN_FREE+1) = MAX(ZCCN_SUM(:,:,:,2) - PSVT(:,:,:,NSV_LIMA_CCN_ACTI+1), 0.)
-
-
-IF (NMOD_CCN .GE. 3) &
-PSVT(:,:,:,NSV_LIMA_CCN_FREE+2) = MAX(ZCCN_SUM(:,:,:,3) - PSVT(:,:,:,NSV_LIMA_CCN_ACTI+2), 0.)
-
-PSVT(:,:,:,NSV_LIMA_IFN_FREE) = MAX(ZIFN_SUM(:,:,:,1) - PSVT(:,:,:,NSV_LIMA_IFN_NUCL), 0.)
-IF (NMOD_IFN .GE. 2) &
-PSVT(:,:,:,NSV_LIMA_IFN_FREE+1) = MAX(ZIFN_SUM(:,:,:,2) - PSVT(:,:,:,NSV_LIMA_IFN_NUCL+1), 0.)
-
+!
+!
+!
+IF (NMOD_CCN .GE. 1) THEN
+ DO JI=1,NMOD_CCN
+ PSVT(:,:,:,NSV_LIMA_CCN_FREE+JI-1) = MAX(ZCCN_SUM(:,:,:,JI) - PSVT(:,:,:,NSV_LIMA_CCN_ACTI+JI-1), 0.)
+ END DO
+END IF
+!
+IF (NMOD_IFN .GE. 1) THEN
+ DO JI=1,NMOD_IFN
+ PSVT(:,:,:,NSV_LIMA_IFN_FREE+JI-1) = MAX(ZIFN_SUM(:,:,:,JI) - PSVT(:,:,:,NSV_LIMA_IFN_NUCL+JI-1), 0.)
+ END DO
+END IF
!
!
END SUBROUTINE AER2LIMA
diff --git a/src/MNH/ini_lima_cold_mixed.f90 b/src/MNH/ini_lima_cold_mixed.f90
index e5cc933aa..2a59420e0 100644
--- a/src/MNH/ini_lima_cold_mixed.f90
+++ b/src/MNH/ini_lima_cold_mixed.f90
@@ -800,7 +800,7 @@ XHMLINTP2 = 1.0 + XHMLINTP1*LOG( 25.E-6/(XGAMINC_HMC_BOUND_MIN)**(1.0/XALPHAC) )
!
!* 7.2 Constants for the accretion of raindrops onto aggregates
!
-XFRACCSS = XNS*((XPI**2)/24.0)*XRHOLW*(ZRHO00**XCEXVT)
+XFRACCSS = XNS*XPI/4.0*XAR*(ZRHO00**XCEXVT)
!
XLBRACCS1 = MOMG(XALPHAS,XNUS,2.)*MOMG(XALPHAR,XNUR,3.)
XLBRACCS2 = 2.*MOMG(XALPHAS,XNUS,1.)*MOMG(XALPHAR,XNUR,4.)
diff --git a/src/MNH/lima_ccn_activation.f90 b/src/MNH/lima_ccn_activation.f90
index 22f5b399a..fa0a276eb 100644
--- a/src/MNH/lima_ccn_activation.f90
+++ b/src/MNH/lima_ccn_activation.f90
@@ -652,7 +652,7 @@ DO JL = 1, NPTS
fl(JL)=fnew(JL)
else if (PX2 .lt. 0.05) then
PX2 = PX2 + 1.0E-2
- PRINT*, 'PX2 ALWAYS too small, we put a greater one : PX2 =',PX2
+! PRINT*, 'PX2 ALWAYS too small, we put a greater one : PX2 =',PX2
fh(JL) = SINGL_FUNCSMAX(PX2,PZZW3(JL),PZZW6(JL),JL)
go to 100
end if
@@ -673,7 +673,7 @@ DO JL = 1, NPTS
PZRIDDR(JL)=PX2
else if (PX2 .lt. 0.05) then
PX2 = PX2 + 1.0E-2
- PRINT*, 'PX2 too small, we put a greater one : PX2 =',PX2
+! PRINT*, 'PX2 too small, we put a greater one : PX2 =',PX2
fh(JL) = SINGL_FUNCSMAX(PX2,PZZW3(JL),PZZW6(JL),JL)
go to 100
else
diff --git a/src/MNH/lima_init_ccn_activation_spectrum.f90 b/src/MNH/lima_init_ccn_activation_spectrum.f90
index 16a561740..4403f9702 100644
--- a/src/MNH/lima_init_ccn_activation_spectrum.f90
+++ b/src/MNH/lima_init_ccn_activation_spectrum.f90
@@ -287,7 +287,7 @@ fh = DSDD(PX2,XDDRY,XKAPPA,XT)
fl=fnew
else if (PX2 .lt. 0.05) then
PX2 = PX2 + 1.0E-2
- PRINT*, 'PX2 ALWAYS too small, we put a greater one : PX2 =',PX2
+! PRINT*, 'PX2 ALWAYS too small, we put a greater one : PX2 =',PX2
fh = DSDD(PX2,XDDRY,XKAPPA,XT)
go to 100
STOP
@@ -303,7 +303,7 @@ fh = DSDD(PX2,XDDRY,XKAPPA,XT)
PZRIDDR=PX2
else if (PX2 .lt. 0.05) then
PX2 = PX2 + 1.0E-2
- PRINT*, 'PX2 too small, we put a greater one : PX2 =',PX2
+! PRINT*, 'PX2 too small, we put a greater one : PX2 =',PX2
fh = DSDD(PX2,XDDRY,XKAPPA,XT)
go to 100
else
diff --git a/src/MNH/lima_mixed_fast_processes.f90 b/src/MNH/lima_mixed_fast_processes.f90
index b8a9fcc32..4ce1a8093 100644
--- a/src/MNH/lima_mixed_fast_processes.f90
+++ b/src/MNH/lima_mixed_fast_processes.f90
@@ -160,7 +160,7 @@ USE MODD_NSV
USE MODD_PARAM_LIMA
USE MODD_PARAM_LIMA_COLD
USE MODD_PARAM_LIMA_MIXED
-USE MODD_PARAM_LIMA_WARM, ONLY : XDR
+USE MODD_PARAM_LIMA_WARM, ONLY : XBR, XDR
use mode_budget, only: Budget_store_init, Budget_store_end
@@ -826,10 +826,10 @@ IF( IGACC>0 .AND. LRAIN) THEN
!
WHERE ( GACC(:) )
ZZW1(:,2) = PCRT1D(:) * & !! coef of RRACCS
- XFRACCSS*( PRST1D(:)*PLBDAS(:)**XBS )*( PRHODREF(:)**(-XCEXVT) ) &
+ XFRACCSS*( PRST1D(:)*PLBDAS(:)**XBS )*( PRHODREF(:)**(1-XCEXVT) ) &
*( XLBRACCS1/((PLBDAS(:)**2) ) + &
XLBRACCS2/( PLBDAS(:) * PLBDAR(:) ) + &
- XLBRACCS3/( (PLBDAR(:)**2)) )/PLBDAR(:)**3
+ XLBRACCS3/( (PLBDAR(:)**2)) )/PLBDAR(:)**XBR
ZZW1(:,4) = MIN( PRRS1D(:),ZZW1(:,2)*ZZW(:) ) ! RRACCSS
PRRS1D(:) = PRRS1D(:) - ZZW1(:,4)
PRSS1D(:) = PRSS1D(:) + ZZW1(:,4)
@@ -870,7 +870,7 @@ IF( IGACC>0 .AND. LRAIN) THEN
WHERE ( GACC(:) .AND. (PRSS1D(:)>XRTMIN(5)/PTSTEP) )
ZZW1(:,2) = MAX( MIN( PRRS1D(:),ZZW1(:,2)-ZZW1(:,4) ) , 0. ) ! RRACCSG
ZZW1(:,3) = MIN( PRSS1D(:),XFSACCRG*ZZW(:)* & ! RSACCRG
- ( PRST1D(:) )*( PRHODREF(:)**(-XCEXVT) ) &
+ ( PRST1D(:) )*( PRHODREF(:)**(1-XCEXVT) ) &
*( XLBSACCR1/((PLBDAR(:)**2) ) + &
XLBSACCR2/( PLBDAR(:) * PLBDAS(:) ) + &
XLBSACCR3/( (PLBDAS(:)**2)) ) )
diff --git a/src/MNH/lima_rain_accr_snow.f90 b/src/MNH/lima_rain_accr_snow.f90
index 264c98bb9..5ece49ac9 100644
--- a/src/MNH/lima_rain_accr_snow.f90
+++ b/src/MNH/lima_rain_accr_snow.f90
@@ -69,6 +69,7 @@ END MODULE MODI_LIMA_RAIN_ACCR_SNOW
!
USE MODD_CST, ONLY : XTT
USE MODD_PARAM_LIMA, ONLY : XRTMIN, XCEXVT
+USE MODD_PARAM_LIMA_WARM, ONLY : XBR
USE MODD_PARAM_LIMA_COLD, ONLY : XBS, XTRANS_MP_GAMMAS
USE MODD_PARAM_LIMA_MIXED, ONLY : NACCLBDAS, XACCINTP1S, XACCINTP2S, &
NACCLBDAR, XACCINTP1R, XACCINTP2R, &
@@ -211,21 +212,19 @@ WHERE( GACC )
!
! 1.3.4 raindrop accretion on the small sized aggregates
!
-! BVIE manque PCRT ???????????????????????????????????
-! ZZW4(:) = & !! coef of RRACCS and RRACCS
- ZZW4(:) = PCRT(:) & !! coef of RRACCS and RRACCS
- * XFRACCSS *( PRST(:)*PLBDS(:)**XBS )*( PRHODREF(:)**(-XCEXVT) ) &
- *( XLBRACCS1/( PLBDS(:)**2 ) + &
- XLBRACCS2/( PLBDS(:) * PLBDR(:) ) + &
- XLBRACCS3/( PLBDR(:)**2 ) ) / PLBDR(:)**3
+ ZZW4(:) = PCRT(:) * & !! coef of RRACCS and RRACCS
+ XFRACCSS *( PRST(:)*PLBDS(:)**XBS )*( PRHODREF(:)**(1-XCEXVT) ) &
+ *( XLBRACCS1/( PLBDS(:)**2 ) + &
+ XLBRACCS2/( PLBDS(:) * PLBDR(:) ) + &
+ XLBRACCS3/( PLBDR(:)**2 ) ) / PLBDR(:)**XBR
!
! 1.3.6 raindrop accretion-conversion of the large sized aggregates
! into graupeln
!
- ZZW5(:) = XFSACCRG*ZZW3(:) * & ! RSACCRG
- ( PRST(:) )*( PRHODREF(:)**(-XCEXVT) ) &
- *( XLBSACCR1/( PLBDR(:)**2 ) + &
+ ZZW5(:) = XFSACCRG * ZZW3(:) * PCRT(:) * & ! RSACCRG
+ ( PRST(:) )*( PRHODREF(:)**(1-XCEXVT) ) * &
+ ( XLBSACCR1/( PLBDR(:)**2 ) + &
XLBSACCR2/( PLBDR(:) * PLBDS(:) ) + &
XLBSACCR3/( PLBDS(:)**2 ) )
!
diff --git a/src/MNH/lima_warm_nucl.f90 b/src/MNH/lima_warm_nucl.f90
index 549a5fc84..d62e4f790 100644
--- a/src/MNH/lima_warm_nucl.f90
+++ b/src/MNH/lima_warm_nucl.f90
@@ -662,7 +662,7 @@ DO JL = 1, NPTS
fl(JL)=fnew(JL)
else if (PX2 .lt. 0.05) then
PX2 = PX2 + 1.0E-2
- PRINT*, 'PX2 ALWAYS too small, we put a greater one : PX2 =',PX2
+! PRINT*, 'PX2 ALWAYS too small, we put a greater one : PX2 =',PX2
fh(JL) = SINGL_FUNCSMAX(PX2,PZZW3(JL),PZZW6(JL),JL)
go to 100
end if
@@ -683,7 +683,7 @@ DO JL = 1, NPTS
PZRIDDR(JL)=PX2
else if (PX2 .lt. 0.05) then
PX2 = PX2 + 1.0E-2
- PRINT*, 'PX2 too small, we put a greater one : PX2 =',PX2
+! PRINT*, 'PX2 too small, we put a greater one : PX2 =',PX2
fh(JL) = SINGL_FUNCSMAX(PX2,PZZW3(JL),PZZW6(JL),JL)
go to 100
else
diff --git a/src/MNH/resolved_cloud.f90 b/src/MNH/resolved_cloud.f90
index 937990ac6..de986a0cd 100644
--- a/src/MNH/resolved_cloud.f90
+++ b/src/MNH/resolved_cloud.f90
@@ -534,8 +534,8 @@ END IF
!
!* 1. From ORILAM to LIMA:
!
-IF (HCLOUD == 'LIMA') THEN
-!IF (HCLOUD == 'LIMA' .AND. ((LORILAM).OR.(LDUST).OR.(LSALT))) THEN
+!IF (HCLOUD == 'LIMA') THEN
+IF (HCLOUD == 'LIMA' .AND. ((LORILAM).OR.(LDUST).OR.(LSALT))) THEN
! ORILAM : tendance s --> variable instant t
ALLOCATE(ZSVT(SIZE(PZZ,1),SIZE(PZZ,2),SIZE(PZZ,3),NSV))
DO II = 1, NSV
--
GitLab