From f0829b1a717101cf64ead847e179bfaf134066d7 Mon Sep 17 00:00:00 2001
From: Quentin Rodier <quentin.rodier@meteo.fr>
Date: Tue, 31 Jan 2023 15:53:17 +0100
Subject: [PATCH] P.Tulet 31/01/2023: add aerosols scavenging with LIMA
---
src/MNH/aer_effic.f90 | 176 +++---
src/MNH/aer_wet_dep_kmt_warm.f90 | 935 ++++++++++++++++++++-----------
src/MNH/aeroopt_get.f90 | 68 ++-
src/MNH/ch_aer_depos.f90 | 47 +-
4 files changed, 786 insertions(+), 440 deletions(-)
diff --git a/src/MNH/aer_effic.f90 b/src/MNH/aer_effic.f90
index 27064451b..8c7b7af94 100644
--- a/src/MNH/aer_effic.f90
+++ b/src/MNH/aer_effic.f90
@@ -16,7 +16,8 @@ SUBROUTINE AER_EFFIC(PRG,PVGG, & !aerosol radius/fall speed (m/s)
PRRS, & ! Rain water m.r. at time
KMODE, & ! Number of aerosol modes
PTEMP, PCOR, & ! air temp, cunningham corr factor
- PDENSITY_AER ) ! aerosol density
+ PDENSITY_AER, & ! aerosol density
+ PRR, PNT ) ! radius and number of rain drops
!
IMPLICIT NONE
REAL, DIMENSION(:,:), INTENT(IN) :: PRG, PVGG
@@ -27,6 +28,7 @@ REAL, DIMENSION(:,:), INTENT(INOUT) :: PEFC
REAL, DIMENSION(:), INTENT(IN) :: PRRS
REAL, DIMENSION(:), INTENT(IN) :: PTEMP
REAL, DIMENSION(:,:), INTENT(IN) :: PCOR
+REAL, DIMENSION(:), INTENT(IN) :: PRR, PNT
INTEGER, INTENT(IN) :: KMODE
REAL, DIMENSION(:,:), INTENT(IN) :: PDENSITY_AER
@@ -40,10 +42,11 @@ SUBROUTINE AER_EFFIC(PRG,PVGG, & !aerosol radius/fall speed (m/s)
PRHODREF, & !Air density
PMUW, PMU, & !mu water/air
PDPG, PEFC, & !diffusivity, efficiency
- PRRS, & ! Rain water m.r. at time t
+ PRRT, & ! Rain water m.r. at time t
KMODE, & ! Number of aerosol modes
PTEMP, PCOR, & ! air temp, cunningham corr factor
- PDENSITY_AER ) ! aerosol density
+ PDENSITY_AER, & ! aerosol density
+ PRR, PNT ) ! radius and number of rain drops
!! #######################################
!!**********AER_EFFIC**********
!! PURPOSE
@@ -71,19 +74,31 @@ SUBROUTINE AER_EFFIC(PRG,PVGG, & !aerosol radius/fall speed (m/s)
!! K. Crahan Kaku / P. Tulet (CNRM/GMEI)
!!
!! MODIFICATIONS
-!! -------------
-!! Philippe Wautelet 28/05/2018: corrected truncated integer division (1/12 -> 1./12.)
+!! -------------
+!! T. Hoarau (LACy) 15/05/17 add LIMA
+!! Philippe Wautelet 28/05/2018: corrected truncated integer division (1/12 -> 1./12.)
+!! P. Tulet and C. Barthe (LAERO) 15/01/22 correction for lima
!!
!-----------------------------------------------------------------
!
!* 0. DECLARATIONS
! ------------
!
-USE MODD_RAIN_ICE_PARAM
-USE MODD_RAIN_ICE_DESCR
-USE MODD_CST, ONLY : XPI, XRHOLW, XP00, XRD
-USE MODD_PARAMETERS , ONLY : JPVEXT
-USE MODD_REF, ONLY : XTHVREFZ
+USE MODD_RAIN_ICE_PARAM, ONLY : YFSEDR => XFSEDR, YEXSEDR => XEXSEDR
+!++cb++
+!++th++
+USE MODD_RAIN_ICE_DESCR, ONLY : YCCR => XCCR, YLBR => XLBR, YLBEXR => XLBEXR, &
+ YCEXVT => XCEXVT
+USE MODD_PARAM_LIMA_WARM, ONLY : WCCR => XCCR, WLBR => XLBR, WLBEXR => XLBEXR, &
+ XFSEDRR, XFSEDRC
+USE MODD_PARAM_LIMA, ONLY : WCEXVT => XCEXVT, WFSEDR => XFSEDR, WFSEDC=>XFSEDC, &
+ XRTMIN
+!--cb--
+USE MODD_PARAM_n, ONLY: CCLOUD
+!--th--
+USE MODD_CST, ONLY : XPI, XRHOLW, XP00, XRD
+USE MODD_PARAMETERS , ONLY : JPVEXT
+USE MODD_REF, ONLY : XTHVREFZ
!
IMPLICIT NONE
!
@@ -93,8 +108,9 @@ REAL, DIMENSION(:), INTENT(IN) :: PRHODREF
REAL, DIMENSION(:,:), INTENT(IN) :: PDPG
REAL, DIMENSION(:), INTENT(IN) :: PMU, PMUW
REAL, DIMENSION(:,:), INTENT(INOUT) :: PEFC
-REAL, DIMENSION(:), INTENT(IN) :: PRRS
+REAL, DIMENSION(:), INTENT(IN) :: PRRT
REAL, DIMENSION(:), INTENT(IN) :: PTEMP
+REAL, DIMENSION(:), INTENT(IN) :: PRR, PNT
REAL, DIMENSION(:,:), INTENT(IN) :: PCOR
INTEGER, INTENT(IN) :: KMODE
REAL, DIMENSION(:,:), INTENT(IN) :: PDENSITY_AER
@@ -110,9 +126,7 @@ REAL, DIMENSION(SIZE(PRG,1)) :: ZOMG, ZREY
!rain radius, m, and rain fall speed, m/s; aerosol radius (m),
REAL, DIMENSION(SIZE(PRG,1)) :: ZRR, ZVR
!lambda, number concentration according to marshall palmer,
-REAL, DIMENSION(SIZE(PRG,1)) :: ZNT, ZLBDA
-! Rain water m.r. source
-REAL, DIMENSION(SIZE(PRG,1)) :: ZRRS
+REAL, DIMENSION(SIZE(PRG,1)) :: ZNT, ZLBDA1
!RHO_dref*r_r, Rain LWC
REAL, DIMENSION(SIZE(PRG,1)) :: RLWC
! schmidts number
@@ -129,83 +143,115 @@ REAL, DIMENSION(SIZE(PRG,1),KMODE) :: ZT1, ZT2
REAL, DIMENSION(SIZE(PRG,1),KMODE) :: ZT3, ZT4
!
INTEGER :: JI,JK
+!++th++
+REAL :: KLBEXR, KLBR, KCEXVT, KCCR, ZFSEDR, ZBR, ZDR, ZEXSEDR
+!--th--
!
!-----------------------------------------------------------------
-ZRRS(:)=PRRS(:)
IKB = 1 + JPVEXT
-ZRHO00 = XP00/(XRD*XTHVREFZ(IKB))
-ZRG(:,:)=PRG(:,:)*1.E-6 !change units to meters
-!
+ZRHO00 = XP00 / (XRD * XTHVREFZ(IKB))
+ZRG(:,:) = PRG(:,:) * 1.E-6 !change units to meters
+ZVR(:) = 0.
+
+SELECT CASE(CCLOUD)
+CASE('ICE3')
+ KLBEXR = YLBEXR
+ KLBR = YLBR
+ KCEXVT = YCEXVT
+ KCCR = YCCR
+ ZFSEDR = YFSEDR
+ ZEXSEDR = YEXSEDR
+
!Fall Speed calculations
!similar to rain_ice.f90, chapter 17.3.4, MESONH Handbook
-!
-ZVR (:)= XFSEDR * ZRRS(:)**(XEXSEDR-1) * &
- PRHODREF(:)**(XEXSEDR-XCEXVT-1)
-! Drop Radius calculation in m
-!lbda = pi*No*rho(lwc)/(rho(dref)*rain rate) p.212 MESONH Handbook
-! compute the slope parameter Lbda_r
-ZLBDA(:) = XLBR*( PRHODREF(:)* ZRRS(:) )**XLBEXR
-!Number concentration NT=No/lbda p. 415 Jacobson
-ZNT (:) = XCCR/ZLBDA (:)
-!rain lwc (kg/m3) = rain m.r.(kg/kg) * rho_air(kg/m3)
-RLWC(:)=ZRRS(:)*PRHODREF(:)
-!4/3 *pi *rł*NT*rho_eau(kg/m3) =rho(lwc)=rho(air)* qc(kg/kg)
-ZRR(:) = (RLWC(:)/(XRHOLW*ZNT(:)*4./3.*XPI))**(1./3.)
-!
+ ZVR(:) = ZFSEDR * PRRT(:)**(ZEXSEDR-1) * &
+ PRHODREF(:)**(ZEXSEDR-KCEXVT)
+
+CASE('LIMA')
+ KLBEXR = WLBEXR
+ KLBR = WLBR
+ KCEXVT = WCEXVT
+ KCCR = WCCR
+ ZFSEDR = XFSEDRR
+ ZBR = 3.0
+ ZDR = 0.8
+ ZEXSEDR = (ZBR + ZDR + 1.0) / (ZBR + 1.0)
+ WHERE (PRRT(:) > XRTMIN(3) .AND. PNT(:) > 0.)
+ ZLBDA1(:) = (KLBR * PNT(:) / PRRT(:))**KLBEXR
+ ZVR(:) = XFSEDRR * PRHODREF(:)**(1.-KCEXVT) * ZLBDA1(:)**(-ZDR)
+ END WHERE
+END SELECT
+
+
!Fall speed cannot be faster than 7 m/s
-ZVR (:)=MIN(ZVR (:),7.)
+ZVR(:) = MIN(ZVR(:), 7.)
+
+KCCR = 8.E6
+
!Ref SEINFELD AND PANDIS p.1019
! Viscosity Ratio
-ZOMG(:)=PMUW(:)/PMU(:)
+ZOMG(:) = PMUW(:) / PMU(:)
!!Reynolds number
-ZREY(:)=ZRR(:)*ZVR(:)*PRHODREF(:)/PMU(:)
-ZREY(:)= MAX(ZREY(:), 1E-2)
-
+ZREY(:) = PRR(:) * ZVR(:) * PRHODREF(:) / PMU(:)
+ZREY(:) = MAX(ZREY(:), 1.E-2)
+!
!S Star
-ZSTA(:)=(1.2+1./12.*LOG(1.+ZREY(:)))/(1.+LOG(1.+ZREY(:)))
-PEFC(:,:)=0.0
-DO JI=1,KMODE
+ZSTA(:) = (1.2 + 1./12. * LOG(1.+ZREY(:))) / (1. + LOG(1.+ZREY(:)))
+
+PEFC(:,:) = 0.0
!
+DO JI = 1, KMODE
!Scmidts number
- ZSCH(:,JI)=PMU(:)/PRHODREF(:)/PDPG(:,JI)
+ ZSCH(:,JI) = PMU(:) / PRHODREF(:) / PDPG(:,JI)
+!
! Rain-Aerosol relative velocity
- ZDIFF(:) = MAX(ZVR(:)-PVGG(:,JI),0.)
+ ZDIFF(:) = MAX(ZVR(:)-PVGG(:,JI), 0.)
+!
! Relaxation time
- ZTAU(:) = (ZRG(:,JI)*2.)**2. * PDENSITY_AER(:,JI) * PCOR(:,JI) / (18.*PMU(:))
+ ZTAU(:) = (ZRG(:,JI)*2.)**2. * PDENSITY_AER(:,JI) * PCOR(:,JI) / (18. * PMU(:))
+!
! Stockes number
- ZSTO(:,JI)= ZTAU(:) * ZDIFF(:) / ZRR(:)
+ ZSTO(:,JI) = ZTAU(:) * ZDIFF(:) / PRR(:)
+!
!Ratio of diameters
- ZPHI(:,JI)=ZRG(:,JI)/ZRR(:)
- ZPHI(:,JI)=MIN(ZPHI(:,JI), 1.)
+ ZPHI(:,JI) = ZRG(:,JI) / PRR(:)
+ ZPHI(:,JI) = MIN(ZPHI(:,JI), 1.)
+!
!Term 1
- ZT1(:,JI)=4.0/ZREY(:)/ZSCH(:,JI)
+ ZT1(:,JI) = 4.0 / ZREY(:) / ZSCH(:,JI)
+!
!Term 2
- ZT2(:,JI)=1.0+0.4*ZREY(:)**(0.5)*ZSCH(:,JI)**(1./3.)+ &
- 0.16*ZREY(:)**(0.5)*ZSCH(:,JI)**(0.5)
-
+ ZT2(:,JI) = 1.0 + 0.4 * ZREY(:)**(0.5) * ZSCH(:,JI)**(1./3.) + &
+ 0.16 * ZREY(:)**(0.5) * ZSCH(:,JI)**(0.5)
+!
!Brownian diffusion
- ZT1(:,JI)= ZT1(:,JI)*ZT2(:,JI)
+ ZT1(:,JI) = ZT1(:,JI) * ZT2(:,JI)
+!
!Term 3 - interception
- ZT3(:,JI)=4.*ZPHI(:,JI)*(1./ZOMG(:)+ &
- (1.0+2.0*ZREY(:)**0.5)*ZPHI(:,JI))
-
- ZT4(:,JI)=0.0
- WHERE(ZSTO(:,JI).GT.ZSTA(:))
+ ZT3(:,JI) = 4. * ZPHI(:,JI) * (1. / ZOMG(:) + &
+ (1.0 + 2.0 * ZREY(:)**0.5) * ZPHI(:,JI))
+!
+ ZT4(:,JI) = 0.0
+!
+ WHERE(ZSTO(:,JI) .GT. ZSTA(:))
!Term 4 - impaction
- ZT4(:,JI)=((ZSTO(:,JI)-ZSTA(:))/ &
- (ZSTO(:,JI)-ZSTA(:)+2./3.))**(3./2.) &
- *(XRHOLW/PDENSITY_AER(:,JI))**(1./2.)
+ ZT4(:,JI) = ((ZSTO(:,JI) - ZSTA(:)) / &
+ (ZSTO(:,JI) - ZSTA(:) + 2. / 3.))**(3./2.) * &
+ (XRHOLW / PDENSITY_AER(:,JI))**(1./2.)
END WHERE
+!
!Collision Efficiancy
- PEFC(:,JI)=ZT1(:,JI)+ ZT3(:,JI)+ZT4(:,JI)
+ PEFC(:,JI) = ZT1(:,JI) + ZT3(:,JI) + ZT4(:,JI)
+!
! Physical radius of a rain collector droplet up than 20 um
-WHERE (ZRR(:) .LE. 20.E-6)
- PEFC(:,JI)= 0.
-END WHERE
+ WHERE (PRR(:) .LE. 20.E-6)
+ PEFC(:,JI) = 0.
+ END WHERE
ENDDO
-PEFC(:,:)=MIN(PEFC(:,:),1.0)
-PEFC(:,:)=MAX(PEFC(:,:),0.0)
+!
+PEFC(:,:) = MIN(PEFC(:,:), 1.0)
+PEFC(:,:) = MAX(PEFC(:,:), 0.0)
END SUBROUTINE AER_EFFIC
diff --git a/src/MNH/aer_wet_dep_kmt_warm.f90 b/src/MNH/aer_wet_dep_kmt_warm.f90
index 8d80ec2e9..0b3ce1984 100644
--- a/src/MNH/aer_wet_dep_kmt_warm.f90
+++ b/src/MNH/aer_wet_dep_kmt_warm.f90
@@ -124,11 +124,31 @@ END MODULE MODI_AER_WET_DEP_KMT_WARM
!
USE MODD_CST
USE MODD_RAIN_ICE_PARAM
-USE MODD_RAIN_ICE_DESCR
+!++th++ 10/05/17
+USE MODD_RAIN_ICE_DESCR, ONLY : YRTMIN => XRTMIN, YCEXVT => XCEXVT, &
+ XCONC_LAND, XCONC_SEA, XCONC_URBAN, &
+ XNUC2, XALPHAC2, XNUC, XALPHAC, &
+ YLBC => XLBC, XLBEXC, &
+ XCCR, &
+ YLBR => XLBR, YLBEXR => XLBEXR
+!--th--
USE MODD_PRECIP_n
USE MODI_AER_VELGRAV
USE MODI_AER_EFFIC
USE MODI_GAMMA
+!++th++ 10/05/17
+USE MODD_PARAM_LIMA, ONLY : XCTMIN, WRTMIN => XRTMIN, WCEXVT => XCEXVT
+USE MODD_PARAM_LIMA_WARM, ONLY : WLBR => XLBR, WLBEXR => XLBEXR, & ! for
+ XFSEDRR, XDR, XBR, & ! sedim.
+ XAUTO1, XAUTO2, XCAUTR, XITAUTR, XLAUTR, & ! for
+ XLAUTR_THRESHOLD, XITAUTR_THRESHOLD, & ! autoconv.
+ WLBC => XLBC, &
+ XACCR1, XACCR2, XACCR3, XACCR4, XACCR5, & ! for
+ XACCR_RLARGE1, XACCR_RLARGE2, & ! accr.
+ XACCR_RSMALL1, XACCR_RSMALL2
+USE MODD_PARAM_n, ONLY: CCLOUD
+!--th--
+
!
IMPLICIT NONE
!
@@ -179,8 +199,17 @@ INTEGER :: ICLOUD, IRAIN
LOGICAL, DIMENSION(SIZE(PSVT,1),SIZE(PSVT,2),SIZE(PSVT,3)) &
:: GRAIN, GCLOUD ! Test where to compute all processes
! Test where to compute the SED/EVAP processes
+!++cb++ 15/05/17
+!REAL, DIMENSION(SIZE(PSVT,1),SIZE(PSVT,2),SIZE(PSVT,3)) &
+! :: ZW, ZZW1, ZZW2, ZZW4 ! work array
REAL, DIMENSION(SIZE(PSVT,1),SIZE(PSVT,2),SIZE(PSVT,3)) &
- :: ZW, ZZW1, ZZW2, ZZW4 ! work array
+ :: ZW, ZZW1, ZZW2, ZZW4, & ! work array
+ ZZW3, ZZW5
+REAL, DIMENSION(SIZE(PSVT,1),SIZE(PSVT,2),SIZE(PSVT,3)) &
+ :: ZDIM, &
+ ZLBDC3, ZLBDC, &
+ ZLBDR3, ZLBDR
+!--cb--
REAL, DIMENSION(SIZE(PSVT,1),SIZE(PSVT,2),SIZE(PSVT,3)) &
:: ZWEVAP ! sedimentation fluxes
REAL, DIMENSION(SIZE(PSVT,1),SIZE(PSVT,2),SIZE(PSVT,3)+1) &
@@ -196,12 +225,11 @@ REAL, DIMENSION(SIZE(PSVT,1),SIZE(PSVT,2),SIZE(PSVT,3)) &
ZFSEDC
REAL, DIMENSION(SIZE(PSVT,1),SIZE(PSVT,2)) :: ZCONC_TMP
REAL, DIMENSION(SIZE(PSVT,1),SIZE(PSVT,2),SIZE(PSVT,3)) :: ZCONC
-REAL, DIMENSION(SIZE(PSVT,1),SIZE(PSVT,2),SIZE(PSVT,3)) :: ZRRS
!
REAL, DIMENSION(:,:), ALLOCATABLE :: ZSVT ! Tracer m.r. concentration
!
REAL, DIMENSION(:,:), ALLOCATABLE :: ZVGG, ZDPG !aerosol velocity [m/s], diffusivity [m2/s]
-REAL, DIMENSION(:,:), ALLOCATABLE :: ZRG !Dust R[ľm]
+REAL, DIMENSION(:,:), ALLOCATABLE :: ZRG !Dust R[\b5m]
REAL, DIMENSION(:,:), ALLOCATABLE :: ZCOR !Cunningham correction factor [unitless]
REAL, DIMENSION(:,:), ALLOCATABLE :: ZMASSMIN ! Aerosol mass minimum value
REAL, DIMENSION(:,:), ALLOCATABLE :: ZDENSITY_AER ! Aerosol density
@@ -221,9 +249,9 @@ REAL, DIMENSION(:), ALLOCATABLE &
ZFLUX, & ! Effective precipitation flux (kg.m-2.s-1)
ZCONC1D, & ! Weighted droplets concentration
ZWLBDC, & ! Slope parameter of the droplet distribution
- ZGAMMA ! scavenging coefficient
-REAL, DIMENSION(:), ALLOCATABLE :: ZW1 ! Work arrays
-REAL, DIMENSION(SIZE(XRTMIN)) :: ZRTMIN
+ ZGAMMA, & ! scavenging coefficient
+ ZLBDA ! lambda parameter for lima distribution
+REAL, DIMENSION(:), ALLOCATABLE :: ZW1 ! Work arrays
INTEGER :: JL ! and PACK intrinsics
!
@@ -231,40 +259,78 @@ INTEGER :: JKAQ, JSV
!
REAL :: A0, A1, A2, A3 ! Constants for computing viscocity
INTEGER :: IKE
-
+!
+REAL, DIMENSION(:), ALLOCATABLE :: KRTMIN
+REAL :: KCEXVT, KLBR, KLBEXR, KLBC, ZLBEXC
+REAL, DIMENSION(2) :: ZXLBC
+REAL :: ZEXSEDR, ZDR
!
!-------------------------------------------------------------------------------
!
!* 0. Initialize work array
! ---------------------
!
+!++cb++ 15/05/17 gestion des parametres redondants entre lima et ice3
+! ATTENTION : pour le moment, les autres schemas microphysiques ne sont pas geres
+! NOTE : les noms sont changes dans toute la routine X... --> K...
+SELECT CASE(CCLOUD)
+CASE('ICE3')
+ ALLOCATE(KRTMIN(SIZE(YRTMIN)))
+ KRTMIN(:) = YRTMIN(:)
+ KCEXVT = YCEXVT
+ KLBR = YLBR
+ KLBEXR = YLBEXR
+ ZXLBC(:) = YLBC(:)
+ ZLBEXC = XLBEXC
+CASE('LIMA')
+ ALLOCATE(KRTMIN(SIZE(WRTMIN)))
+ KRTMIN = WRTMIN
+ KCEXVT = WCEXVT
+ KLBR = WLBR
+ KLBEXR = WLBEXR
+ KLBC = WLBC
+ ZLBEXC = 1.0 / 3.0
+ ZDR = 0.8
+END SELECT
+!--cb--
+!
! Compute Effective cloud radius
-ZRAY(:,:,:) = 0.
-ZLBC(:,:,:) = 0.
-IF (PRESENT(PCCT)) THEN ! case KHKO, C2R2, C3R5 (two moments schemes)
- ZRAY(:,:,:) = 3.* PRCT(:,:,:) / (4.*XPI*XRHOLW*PCCT(:,:,:))
- ZRAY(:,:,:) = ZRAY(:,:,:)**(1./3.) ! Cloud mean radius in m
-ELSE IF (PRESENT(PSEA)) THEN ! Case ICE3, RAVE, KESS, ..
-ZLBC(:,:,:) = XLBC(1)
-ZFSEDC(:,:,:) = XFSEDC(1)
-ZCONC(:,:,:) = XCONC_LAND
-ZCONC_TMP(:,:)=PSEA(:,:)*XCONC_SEA+(1.-PSEA(:,:))*XCONC_LAND
-DO JK=1,SIZE(PRHODREF,3)
- ZLBC(:,:,JK) = PSEA(:,:)*XLBC(2)+(1.-PSEA(:,:))*XLBC(1)
- ZFSEDC(:,:,JK) = (PSEA(:,:)*XFSEDC(2)+(1.-PSEA(:,:))*XFSEDC(1))
- ZFSEDC(:,:,JK) = MAX(MIN(XFSEDC(1),XFSEDC(2)),ZFSEDC(:,:,JK))
- ZCONC(:,:,JK) = (1.-PTOWN(:,:))*ZCONC_TMP(:,:)+PTOWN(:,:)*XCONC_URBAN
- ZRAY(:,:,JK) = 0.5*((1.-PSEA(:,:))*GAMMA(XNUC+1.0/XALPHAC)/(GAMMA(XNUC)) + &
- PSEA(:,:)*GAMMA(XNUC2+1.0/XALPHAC2)/(GAMMA(XNUC2)))
-END DO
-ZRAY(:,:,:) = MAX(1.,ZRAY(:,:,:))
-ZLBC(:,:,:) = MAX(MIN(XLBC(1),XLBC(2)),ZLBC(:,:,:))
+ZRAY(:,:,:) = 0.
+ZLBC(:,:,:) = 0.
+!
+!++th++ 05/05/17 test thomas
+IF (PRESENT(PCCT)) THEN ! case KHKO, C2R2, C3R5, LIMA (two moments schemes)
+!
+ WHERE (PCCT(:,:,:) .GT. 0. .AND. PRCT(:,:,:) .GT. 0.)
+ ZRAY(:,:,:) = 3. * PRCT(:,:,:) / (4. * XPI * XRHOLW * PCCT(:,:,:))
+ ZRAY(:,:,:) = ZRAY(:,:,:)**(1./3.) ! Cloud mean radius in m
+ ELSEWHERE
+ ZRAY(:,:,:) = 30. ! Cloud mean radius in m
+ ENDWHERE
+!--th--
+!
+ELSE IF (PRESENT(PSEA)) THEN ! Case ICE3, REVE, KESS, ..
+ ZLBC(:,:,:) = ZXLBC(1)
+ ZFSEDC(:,:,:) = XFSEDC(1)
+ ZCONC(:,:,:) = XCONC_LAND
+ ZCONC_TMP(:,:) = PSEA(:,:) * XCONC_SEA + (1. - PSEA(:,:)) * XCONC_LAND
+!
+ DO JK = 1, SIZE(PRHODREF,3)
+ ZLBC(:,:,JK) = PSEA(:,:) * ZXLBC(2) + (1. - PSEA(:,:)) * ZXLBC(1)
+ ZFSEDC(:,:,JK) = (PSEA(:,:) * XFSEDC(2) + (1. - PSEA(:,:)) * XFSEDC(1))
+ ZFSEDC(:,:,JK) = MAX(MIN(XFSEDC(1),XFSEDC(2)),ZFSEDC(:,:,JK))
+ ZCONC(:,:,JK) = (1. - PTOWN(:,:)) * ZCONC_TMP(:,:) + PTOWN(:,:) * XCONC_URBAN
+ ZRAY(:,:,JK) = 0.5 * ((1. - PSEA(:,:)) * GAMMA(XNUC+1.0/XALPHAC) / (GAMMA(XNUC)) + &
+ PSEA(:,:) * GAMMA(XNUC2+1.0/XALPHAC2) / (GAMMA(XNUC2)))
+ END DO
+ ZRAY(:,:,:) = MAX(1., ZRAY(:,:,:))
+ ZLBC(:,:,:) = MAX(MIN(ZXLBC(1),ZXLBC(2)), ZLBC(:,:,:))
ELSE
-ZRAY(:,:,:) = 30. ! default value for cloud radius
+ ZRAY(:,:,:) = 30. ! default value for cloud radius
END IF
!
-ZNRT(:,:,:) = 0.
-IF (PRESENT(PCRT)) THEN ! case KHKO, C2R2, C3R5
+ZNRT(:,:,:) = 0.
+IF (PRESENT(PCRT)) THEN ! case KHKO, C2R2, C3R5, LIMA
! Transfert Number of rain droplets
ZNRT(:,:,:) = PCRT(:,:,:)
END IF
@@ -272,7 +338,9 @@ END IF
!
!* 1. COMPUTE THE AEROSOL/CLOUD-RAIN MASS TRANSFER
! ----------------------------------------------
+!
CALL AER_WET_MASS_TRANSFER
+!
!-------------------------------------------------------------------------------
!
!* 2. COMPUTE THE SEDIMENTATION (RS) SOURCE
@@ -281,20 +349,22 @@ CALL AER_WET_MASS_TRANSFER
CALL AER_WET_DEP_KMT_WARM_SEDIMENT
!
!-------------------------------------------------------------------------------
-!!
-!!* 3. COMPUTES THE SLOW WARM PROCESS SOURCES
-!! --------------------------------------
-!!
+!
+!* 3. COMPUTES THE SLOW WARM PROCESS SOURCES
+! --------------------------------------
+!
CALL AER_WET_DEP_KMT_ICE_WARM
-
!
!-------------------------------------------------------------------------------
-!!* 4. COMPUTES EVAPORATION PROCESS
-!! ----------------------------
-!!
+!* 4. COMPUTES EVAPORATION PROCESS
+! ----------------------------
+!
CALL AER_WET_DEP_KMT_EVAP
!
-
+!++cb++
+DEALLOCATE(KRTMIN)
+!--cb--
+!
!-------------------------------------------------------------------------------
!
!
@@ -324,208 +394,241 @@ INTEGER :: JKAQ ! counter for chemistry
! 1 Mass transfer Aerosol to cloud (Tost et al., 2006)
!
GCLOUD(:,:,:) = .FALSE.
-GCLOUD(:,:,:) = PRCS(:,:,:)>XRTMIN(2)
+!
+IF (PRESENT(PCCT)) THEN ! case KHKO, C2R2, C3R5, LIMA (2-moment schemes)
+ GCLOUD(:,:,:) = PRCS(:,:,:) > KRTMIN(2) .AND. PCCT(:,:,:) > XCTMIN(2)
+ELSE ! Case ICE3, REVE, KESS, ... (1-moment schemes)
+ GCLOUD(:,:,:) = PRCS(:,:,:) > KRTMIN(2)
+END IF
+!--cb--
+!--th--
+
ICLOUD = COUNTJV( GCLOUD(:,:,:),I1C(:),I2C(:),I3C(:))
IF( ICLOUD >= 1 ) THEN
-ALLOCATE(ZSVT(ICLOUD,KMODE*3))
-ALLOCATE(ZRHODREF(ICLOUD))
-ALLOCATE(ZTHT(ICLOUD))
-ALLOCATE(ZRC(ICLOUD))
-ALLOCATE(ZPABST(ICLOUD))
-ALLOCATE(ZRG(ICLOUD,KMODE))
-ALLOCATE(ZTEMP(ICLOUD))
-ALLOCATE(ZMU(ICLOUD))
-ALLOCATE(ZRCT(ICLOUD))
-ALLOCATE(ZVGG(ICLOUD,KMODE))
-ALLOCATE(ZDPG(ICLOUD,KMODE))
-ALLOCATE(ZGAMMA(ICLOUD))
-ALLOCATE(ZW1(ICLOUD))
-ALLOCATE(ZCOR(ICLOUD,KMODE))
-ALLOCATE(ZMASSMIN(ICLOUD,KMODE))
-ALLOCATE(ZWLBDC(ICLOUD))
-ALLOCATE(ZCONC1D(ICLOUD))
-ALLOCATE(ZDENSITY_AER(ICLOUD,KMODE))
-ZSVT(:,:) = 0.
-DO JL=1,ICLOUD
- DO JKAQ = 1, KMODE
- ZRG(JL,JKAQ) = PRGAER(I1C(JL),I2C(JL),I3C(JL),JKAQ)
- ENDDO
- DO JKAQ = 1, KMODE*3
- ZSVT(JL,JKAQ) = PSVT(I1C(JL),I2C(JL),I3C(JL),JKAQ)
- END DO
+ ALLOCATE(ZSVT(ICLOUD,KMODE*3))
+ ALLOCATE(ZRHODREF(ICLOUD))
+ ALLOCATE(ZTHT(ICLOUD))
+ ALLOCATE(ZRC(ICLOUD))
+ ALLOCATE(ZPABST(ICLOUD))
+ ALLOCATE(ZRG(ICLOUD,KMODE))
+ ALLOCATE(ZTEMP(ICLOUD))
+ ALLOCATE(ZMU(ICLOUD))
+ ALLOCATE(ZRCT(ICLOUD))
+ ALLOCATE(ZVGG(ICLOUD,KMODE))
+ ALLOCATE(ZDPG(ICLOUD,KMODE))
+ ALLOCATE(ZGAMMA(ICLOUD))
+ ALLOCATE(ZW1(ICLOUD))
+ ALLOCATE(ZCOR(ICLOUD,KMODE))
+ ALLOCATE(ZMASSMIN(ICLOUD,KMODE))
+ ALLOCATE(ZWLBDC(ICLOUD))
+ ALLOCATE(ZCONC1D(ICLOUD))
+ ALLOCATE(ZDENSITY_AER(ICLOUD,KMODE))
+!
+ ZSVT(:,:) = 0.
+!
+ DO JL = 1, ICLOUD
+ DO JKAQ = 1, KMODE
+ ZRG(JL,JKAQ) = PRGAER(I1C(JL),I2C(JL),I3C(JL),JKAQ)
+ ENDDO
+ DO JKAQ = 1, KMODE*3
+ ZSVT(JL,JKAQ) = PSVT(I1C(JL),I2C(JL),I3C(JL),JKAQ)
+ END DO
!
- ZTHT(JL) = PTHT(I1C(JL),I2C(JL),I3C(JL))
- ZRC(JL) = ZRAY(I1C(JL),I2C(JL),I3C(JL))
- ZPABST(JL) = PPABST(I1C(JL),I2C(JL),I3C(JL))
- ZRCT(JL) = PRCS(I1C(JL),I2C(JL),I3C(JL))
- ZRHODREF(JL) = PRHODREF(I1C(JL),I2C(JL),I3C(JL))
- ZMASSMIN(JL,:) = PMASSMIN(I1C(JL),I2C(JL),I3C(JL),:)
- ZWLBDC(JL) = ZLBC(I1C(JL),I2C(JL),I3C(JL))
- ZCONC1D(JL) = ZCONC(I1C(JL),I2C(JL),I3C(JL))
- ZDENSITY_AER(JL,:) = PDENSITY_AER(I1C(JL),I2C(JL),I3C(JL),:)
-END DO
-IF (ANY(ZWLBDC(:)/=0.)) THEN ! case one moments
+ ZTHT(JL) = PTHT(I1C(JL),I2C(JL),I3C(JL))
+ ZRC(JL) = ZRAY(I1C(JL),I2C(JL),I3C(JL))
+ ZPABST(JL) = PPABST(I1C(JL),I2C(JL),I3C(JL))
+ ZRCT(JL) = PRCS(I1C(JL),I2C(JL),I3C(JL))
+ ZRHODREF(JL) = PRHODREF(I1C(JL),I2C(JL),I3C(JL))
+ ZMASSMIN(JL,:) = PMASSMIN(I1C(JL),I2C(JL),I3C(JL),:)
+ ZWLBDC(JL) = ZLBC(I1C(JL),I2C(JL),I3C(JL))
+ ZCONC1D(JL) = ZCONC(I1C(JL),I2C(JL),I3C(JL))
+ ZDENSITY_AER(JL,:) = PDENSITY_AER(I1C(JL),I2C(JL),I3C(JL),:)
+ END DO
+!
+ IF (ANY(ZWLBDC(:) /= 0.)) THEN ! case one moments
! On calcule Rc a partir de M(3) car c'est le seul moment indt de alpha et nu
! Rho_air * Rc / (Pi/6 * Rho_eau * Nc) = M(3) = 1/ (Lambda**3 * rapport des
! gamma)
- ZWLBDC(:) = ZWLBDC(:) * ZCONC1D(:) / (ZRHODREF(:) * ZRCT(:))
- ZWLBDC(:) = ZWLBDC(:)**XLBEXC
- ZRC(:) = ZRC(:) / ZWLBDC(:)
-END IF
-
+ ZWLBDC(:) = ZWLBDC(:) * ZCONC1D(:) / (ZRHODREF(:) * ZRCT(:))
+ ZWLBDC(:) = ZWLBDC(:)**ZLBEXC
+ ZRC(:) = ZRC(:) / ZWLBDC(:)
+ END IF
!
! initialize temperature
- ZTEMP(:)=ZTHT(:)*(ZPABST(:)/XP00)**(XRD/XCPD)
-
+ ZTEMP(:) = ZTHT(:) * (ZPABST(:) / XP00)**(XRD/XCPD)
+!
! compute diffusion and gravitation velocity
+ CALL AER_VELGRAV(ZRG(:,:), ZPABST(:), &
+ KMODE, ZMU(:), ZVGG(:,:), &
+ ZDPG(:,:),ZTEMP(:),ZCOR(:,:), &
+ ZDENSITY_AER(:,:))
- CALL AER_VELGRAV(ZRG(:,:), ZPABST(:), &
- KMODE, ZMU(:), ZVGG(:,:), &
- ZDPG(:,:),ZTEMP(:),ZCOR(:,:), &
- ZDENSITY_AER(:,:))
-
-DO JKAQ = 1, KMODE
+ DO JKAQ = 1, KMODE
! Browninan nucleation scavenging (Pruppacher and Klett, 2000, p723)
- ZGAMMA(:) = 1.35 * ZRCT(:)*ZRHODREF(:)*1.E-3 * ZDPG(:,JKAQ) /&
- (ZRC(:)*ZRC(:))
-
- ZW1(:) = ZSVT(:,JKAQ) * EXP(-ZGAMMA(:) * PTSTEP)
- ZW1(:) = MAX(ZW1(:), ZMASSMIN(:,JKAQ))
- ZW1(:) = MIN(ZW1(:),ZSVT(:,JKAQ))
+ ZGAMMA(:) = 1.35 * ZRCT(:) * ZRHODREF(:) * 1.E-3 * ZDPG(:,JKAQ) / &
+ (ZRC(:) * ZRC(:))
+!
+ ZW1(:) = ZSVT(:,JKAQ) * EXP(-ZGAMMA(:) * PTSTEP)
+ ZW1(:) = MAX(ZW1(:), ZMASSMIN(:,JKAQ))
+! ZW1(:) = MIN(ZW1(:), ZSVT(:,JKAQ))
! Aerosol mass in cloud
- ZSVT(:,KMODE+JKAQ) = ZSVT(:,KMODE+JKAQ) + ZSVT(:,JKAQ) - ZW1(:)
+ ZSVT(:,KMODE+JKAQ) = ZSVT(:,KMODE+JKAQ) + ZSVT(:,JKAQ) - ZW1(:)
! New aerosol mass
- ZSVT(:,JKAQ) = ZW1(:)
+ ZSVT(:,JKAQ) = ZW1(:)
! Return in 3D
- PSVT(:,:,:,JKAQ) = &
- UNPACK(ZSVT(:,JKAQ),MASK=GCLOUD(:,:,:),FIELD=PSVT(:,:,:,JKAQ))
- PSVT(:,:,:,KMODE+JKAQ) = &
- UNPACK(ZSVT(:,KMODE+JKAQ),MASK=GCLOUD(:,:,:),FIELD=PSVT(:,:,:,KMODE+JKAQ))
-ENDDO
-DEALLOCATE(ZSVT)
-DEALLOCATE(ZRHODREF)
-DEALLOCATE(ZTHT)
-DEALLOCATE(ZRC)
-DEALLOCATE(ZPABST)
-DEALLOCATE(ZRG)
-DEALLOCATE(ZTEMP)
-DEALLOCATE(ZMU)
-DEALLOCATE(ZRCT)
-DEALLOCATE(ZVGG)
-DEALLOCATE(ZDPG)
-DEALLOCATE(ZGAMMA)
-DEALLOCATE(ZW1)
-DEALLOCATE(ZCOR)
-DEALLOCATE(ZMASSMIN)
-DEALLOCATE(ZWLBDC)
-DEALLOCATE(ZCONC1D)
-DEALLOCATE(ZDENSITY_AER)
+ PSVT(:,:,:,JKAQ) = &
+ UNPACK(ZSVT(:,JKAQ),MASK=GCLOUD(:,:,:),FIELD=PSVT(:,:,:,JKAQ))
+ PSVT(:,:,:,KMODE+JKAQ) = &
+ UNPACK(ZSVT(:,KMODE+JKAQ),MASK=GCLOUD(:,:,:),FIELD=PSVT(:,:,:,KMODE+JKAQ))
+ ENDDO
+!
+ DEALLOCATE(ZSVT)
+ DEALLOCATE(ZRHODREF)
+ DEALLOCATE(ZTHT)
+ DEALLOCATE(ZRC)
+ DEALLOCATE(ZPABST)
+ DEALLOCATE(ZRG)
+ DEALLOCATE(ZTEMP)
+ DEALLOCATE(ZMU)
+ DEALLOCATE(ZRCT)
+ DEALLOCATE(ZVGG)
+ DEALLOCATE(ZDPG)
+ DEALLOCATE(ZGAMMA)
+ DEALLOCATE(ZW1)
+ DEALLOCATE(ZCOR)
+ DEALLOCATE(ZMASSMIN)
+ DEALLOCATE(ZWLBDC)
+ DEALLOCATE(ZCONC1D)
+ DEALLOCATE(ZDENSITY_AER)
END IF
!
! 2 Mass transfer Aerosol to Rain (Seinfeld and Pandis, 1998, Tost et al., 2006)
!
GRAIN(:,:,:) = .FALSE.
-GRAIN(:,:,:) = PRRS(:,:,:)>XRTMIN(3)
+!
+IF (PRESENT(PCRT)) THEN ! case KHKO, C2R2, C3R5, LIMA (2-moment schemes)
+ GRAIN(:,:,:) = PRRT(:,:,:) > KRTMIN(3) .AND. PCRT(:,:,:) > XCTMIN(3)
+ELSE ! Case ICE3, REVE, KESS, ... (1-moment schemes)
+ GRAIN(:,:,:) = PRRT(:,:,:) > KRTMIN(3)
+END IF
+
IRAIN = COUNTJV( GRAIN(:,:,:),I1R(:),I2R(:),I3R(:))
IF( IRAIN >= 1 ) THEN
!
- ALLOCATE(ZRRT(IRAIN))
- ALLOCATE(ZSVT(IRAIN,3*KMODE))
- ALLOCATE(ZRHODREF(IRAIN))
- ALLOCATE(ZTHT(IRAIN))
- ALLOCATE(ZRR(IRAIN))
- ALLOCATE(ZNT(IRAIN))
- ALLOCATE(ZPABST(IRAIN))
- ALLOCATE(ZRG(IRAIN,KMODE))
- ALLOCATE(ZCOR(IRAIN,KMODE))
- ALLOCATE(ZTEMP(IRAIN))
- ALLOCATE(ZMU(IRAIN))
- ALLOCATE(ZVGG(IRAIN,KMODE))
- ALLOCATE(ZDPG(IRAIN,KMODE))
- ALLOCATE(ZMUW(IRAIN))
- ALLOCATE(ZEFC(IRAIN,KMODE))
- ALLOCATE(ZW1(IRAIN))
- ALLOCATE(ZFLUX(IRAIN))
- ALLOCATE(ZGAMMA(IRAIN))
- ALLOCATE(ZMASSMIN(IRAIN,KMODE))
- ALLOCATE(ZDENSITY_AER(IRAIN,KMODE))
-
- ZSVT(:,:) = 0.
- DO JL=1,IRAIN
- DO JKAQ = 1, KMODE
- ZRG(JL,JKAQ) = PRGAER(I1R(JL),I2R(JL),I3R(JL),JKAQ )
- ZSVT(JL,JKAQ) = PSVT(I1R(JL),I2R(JL),I3R(JL),JKAQ)
- ZSVT(JL,KMODE*2+JKAQ) = PSVT(I1R(JL),I2R(JL),I3R(JL),KMODE*2+JKAQ)
- END DO
- !
- ZTHT(JL) = PTHT(I1R(JL),I2R(JL),I3R(JL))
- ZPABST(JL) = PPABST(I1R(JL),I2R(JL),I3R(JL))
- ZRRT(JL) = PRRS(I1R(JL),I2R(JL),I3R(JL))
- ZRHODREF(JL) = PRHODREF(I1R(JL),I2R(JL),I3R(JL))
- ZMASSMIN(JL,:) = PMASSMIN(I1R(JL),I2R(JL),I3R(JL),:)
- ZNT(JL) = ZNRT(I1R(JL),I2R(JL),I3R(JL))
- ZDENSITY_AER(JL,:) = PDENSITY_AER(I1R(JL),I2R(JL),I3R(JL),:)
- ENDDO
-!
-CALL AER_WET_DEP_KMT_EFFIC
+ ALLOCATE(ZRRT(IRAIN))
+ ALLOCATE(ZSVT(IRAIN,3*KMODE))
+ ALLOCATE(ZRHODREF(IRAIN))
+ ALLOCATE(ZTHT(IRAIN))
+ ALLOCATE(ZRR(IRAIN))
+ ALLOCATE(ZNT(IRAIN))
+ ALLOCATE(ZPABST(IRAIN))
+ ALLOCATE(ZRG(IRAIN,KMODE))
+ ALLOCATE(ZCOR(IRAIN,KMODE))
+ ALLOCATE(ZTEMP(IRAIN))
+ ALLOCATE(ZMU(IRAIN))
+ ALLOCATE(ZVGG(IRAIN,KMODE))
+ ALLOCATE(ZDPG(IRAIN,KMODE))
+ ALLOCATE(ZMUW(IRAIN))
+ ALLOCATE(ZEFC(IRAIN,KMODE))
+ ALLOCATE(ZW1(IRAIN))
+ ALLOCATE(ZFLUX(IRAIN))
+ ALLOCATE(ZGAMMA(IRAIN))
+ ALLOCATE(ZMASSMIN(IRAIN,KMODE))
+ ALLOCATE(ZDENSITY_AER(IRAIN,KMODE))
+ ALLOCATE(ZLBDA(IRAIN))
+!
+ ZSVT(:,:) = 0.
+!
+ DO JL = 1, IRAIN
+ DO JKAQ = 1, KMODE
+ ZRG(JL,JKAQ) = PRGAER(I1R(JL),I2R(JL),I3R(JL),JKAQ )
+ ZSVT(JL,JKAQ) = PSVT(I1R(JL),I2R(JL),I3R(JL),JKAQ)
+ ZSVT(JL,KMODE*2+JKAQ) = PSVT(I1R(JL),I2R(JL),I3R(JL),KMODE*2+JKAQ)
+ END DO
+!
+ ZTHT(JL) = PTHT(I1R(JL),I2R(JL),I3R(JL))
+ ZPABST(JL) = PPABST(I1R(JL),I2R(JL),I3R(JL))
+ ZRRT(JL) = PRRT(I1R(JL),I2R(JL),I3R(JL))
+ ZRHODREF(JL) = PRHODREF(I1R(JL),I2R(JL),I3R(JL))
+ ZMASSMIN(JL,:) = PMASSMIN(I1R(JL),I2R(JL),I3R(JL),:)
+ ZNT(JL) = ZNRT(I1R(JL),I2R(JL),I3R(JL))
+ ZDENSITY_AER(JL,:) = PDENSITY_AER(I1R(JL),I2R(JL),I3R(JL),:)
+ ENDDO
! Compute scavenging coefficient
-ZFLUX(:) = 0.
-ZRRT(:) = MAX(ZRRT(:), 0.)
+ ZFLUX(:) = 0.
+ ZRRT(:) = MAX(ZRRT(:), 0.)
+!
! Effective precipitation flux (kg.m-2.s-1)
-ZFLUX(:) = XFSEDR * ZRRT(:)**(XEXSEDR ) &
- * ZRHODREF(:)**(XEXSEDR-XCEXVT)
-ZFLUX(:) = MAX(ZFLUX(:), 0.)
-
-IF (ALL(ZNT(:) == 0.)) THEN ! case one moments
-!Number concentration NT=No/lbda p. 415 Jacobson
-!4/3 *pi *rł*NT*rho_eau(kg/m3) =rho(lwc)=rho(air)* qc(kg/kg)
-ZNT (:) = XCCR/(XLBR*( ZRHODREF(:)* ZRRT(:) )**XLBEXR)
-END IF
+ IF (PRESENT(PCRT)) THEN ! cf lima_precip_scavenging.f90 (l. 751)
+ ZEXSEDR = (XBR + XDR + 1.0) / (XBR + 1.0)
-ZRR(:) = (ZRRT(:)*ZRHODREF(:)/(XRHOLW*ZNT(:)*4./3.*XPI))**(1./3.)
+ ZLBDA(:) = (KLBR * ZNT(:) / ZRRT(:))**KLBEXR
+ ZFLUX(:) = XFSEDRR * ZRRT(:) * ZRHODREF(:)**(1.-KCEXVT) * ZLBDA(:)**(-ZDR)
-DO JKAQ = 1, KMODE
- ! Tost et al, 2006
- ZGAMMA(:) = 0.75 * ZEFC(:,JKAQ) * ZFLUX(:) / (ZRR(:)*1E3)
- ZW1(:) = ZSVT(:,JKAQ) * EXP(-ZGAMMA(:) * PTSTEP)
- ZW1(:) = MAX(ZW1(:), ZMASSMIN(:,JKAQ))
- ZW1(:) = MIN(ZW1(:),ZSVT(:,JKAQ))
-
- ! Aerosol mass in rain
- ZSVT(:,KMODE*2+JKAQ) = ZSVT(:,KMODE*2+JKAQ) + ZSVT(:,JKAQ) - ZW1(:)
- ! New aerosol mass
- ZSVT(:,JKAQ) = ZW1(:)
-
- ! Return to 3D
- PSVT(:,:,:,JKAQ) = &
- UNPACK(ZSVT(:,JKAQ),MASK=GRAIN(:,:,:),FIELD=PSVT(:,:,:,JKAQ))
- PSVT(:,:,:,KMODE*2+JKAQ) = &
- UNPACK(ZSVT(:,KMODE*2+JKAQ),MASK=GRAIN(:,:,:),FIELD=PSVT(:,:,:,KMODE*2+JKAQ))
-ENDDO
- DEALLOCATE(ZRRT)
- DEALLOCATE(ZSVT)
- DEALLOCATE(ZRHODREF)
- DEALLOCATE(ZTHT)
- DEALLOCATE(ZRR)
- DEALLOCATE(ZNT)
- DEALLOCATE(ZPABST)
- DEALLOCATE(ZRG)
- DEALLOCATE(ZCOR)
- DEALLOCATE(ZTEMP)
- DEALLOCATE(ZMU)
- DEALLOCATE(ZVGG)
- DEALLOCATE(ZDPG)
- DEALLOCATE(ZMUW)
- DEALLOCATE(ZEFC)
- DEALLOCATE(ZW1)
- DEALLOCATE(ZFLUX)
- DEALLOCATE(ZGAMMA)
- DEALLOCATE(ZMASSMIN)
- DEALLOCATE(ZDENSITY_AER)
- END IF
+ ELSE ! cf ZWSED dans rain_ice.f90 (l. 1077)
+ ZFLUX(:) = XFSEDR * ZRRT(:)**(XEXSEDR) * ZRHODREF(:)**(XEXSEDR-KCEXVT)
+ END IF
+ ZFLUX(:) = MAX(ZFLUX(:), 0.)
+
+ IF (ALL(ZNT(:) == 0.)) THEN ! case one moments
+! Number concentration NT=No/lbda p. 415 Jacobson
+! 4/3 *pi *r\b3*NT*rho_eau(kg/m3) =rho(lwc)=rho(air)* qc(kg/kg)
+ ZNT(:) = XCCR / (KLBR * (ZRHODREF(:) * ZRRT(:))**KLBEXR)
+ END IF
+!
+ ZRR(:) = (ZRRT(:) * ZRHODREF(:) / &
+ (XRHOLW * ZNT(:) * 4. / 3. * XPI))**(1./3.)
+
+ CALL AER_WET_DEP_KMT_EFFIC
+
+ DO JKAQ = 1, KMODE
+ ! Tost et al, 2006
+ ZGAMMA(:) = 0.75 * ZEFC(:,JKAQ) * ZFLUX(:) / (ZRR(:) * 1.E3)
+
+ ZW1(:) = ZSVT(:,JKAQ) * EXP(-ZGAMMA(:)*PTSTEP)
+ ZW1(:) = MAX(ZW1(:), ZMASSMIN(:,JKAQ))
+
+ ! Aerosol mass in rain
+ ZSVT(:,KMODE*2+JKAQ) = ZSVT(:,KMODE*2+JKAQ) + ZSVT(:,JKAQ) - ZW1(:)
+
+ ! New aerosol mass
+ ZSVT(:,JKAQ) = ZW1(:)
+
+ ! Return to 3D
+ PSVT(:,:,:,JKAQ) = &
+ UNPACK(ZSVT(:,JKAQ),MASK=GRAIN(:,:,:),FIELD=PSVT(:,:,:,JKAQ))
+ PSVT(:,:,:,KMODE*2+JKAQ) = &
+ UNPACK(ZSVT(:,KMODE*2+JKAQ),MASK=GRAIN(:,:,:),FIELD=PSVT(:,:,:,KMODE*2+JKAQ))
+ ENDDO
+!
+ DEALLOCATE(ZRRT)
+ DEALLOCATE(ZSVT)
+ DEALLOCATE(ZRHODREF)
+ DEALLOCATE(ZTHT)
+ DEALLOCATE(ZRR)
+ DEALLOCATE(ZNT)
+ DEALLOCATE(ZPABST)
+ DEALLOCATE(ZRG)
+ DEALLOCATE(ZCOR)
+ DEALLOCATE(ZTEMP)
+ DEALLOCATE(ZMU)
+ DEALLOCATE(ZVGG)
+ DEALLOCATE(ZDPG)
+ DEALLOCATE(ZMUW)
+ DEALLOCATE(ZEFC)
+ DEALLOCATE(ZW1)
+ DEALLOCATE(ZFLUX)
+ DEALLOCATE(ZGAMMA)
+ DEALLOCATE(ZMASSMIN)
+ DEALLOCATE(ZDENSITY_AER)
+ DEALLOCATE(ZLBDA)
+END IF
+!
END SUBROUTINE AER_WET_MASS_TRANSFER
!
+!-------------------------------------------------------------------------------
+!
SUBROUTINE AER_WET_DEP_KMT_WARM_SEDIMENT
!
!* Sedimentation of aerosol in rain droplets
@@ -533,52 +636,150 @@ SUBROUTINE AER_WET_DEP_KMT_WARM_SEDIMENT
!* 0. DECLARATIONS
! ------------
!
+use mode_tools, only: Countjv
+!
IMPLICIT NONE
!
!* declaration of local variables
!
-!
-INTEGER :: JL ! and PACK intrinsics
-INTEGER :: JKAQ ! counter for acquous aerosols
+INTEGER :: JL ! and PACK intrinsics
+INTEGER :: JKAQ ! counter for acquous aerosols
+INTEGER :: IRAIN, ILISTLENR
+INTEGER :: ILENALLOCR
+INTEGER, SAVE :: IOLDALLOCR = 6000
+INTEGER, DIMENSION(SIZE(PZZ)) :: IR1,IR2,IR3 ! Used to replace the COUNT
+INTEGER, DIMENSION(:), ALLOCATABLE :: ILISTR
+REAL, DIMENSION(:), ALLOCATABLE :: ZLAMBDA, ZRHODREF, ZCRT, ZRRT
+REAL, DIMENSION(:,:), ALLOCATABLE :: ZSVT
!
!-------------------------------------------------------------------------------
!
!* Time splitting initialization
+!
ZTSPLITR = PTSTEP / REAL(KSPLITR)
!
ZW(:,:,:)=0.
-ZRRS(:,:,:) = MAX(PRRS(:,:,:), 0.)
-IKE = SIZE(PRCS,3)
+ZWSED(:,:,:) = 0.
+IKE = SIZE(PRCT,3)
+ILENALLOCR = 0
DO JK = 1 , SIZE(PZZ,3)-1
- ZW(:,:,JK) =ZTSPLITR/(( PZZ(:,:,JK+1)-PZZ(:,:,JK) ))
+ ZW(:,:,JK) = ZTSPLITR / ((PZZ(:,:,JK+1) - PZZ(:,:,JK)))
END DO
-WHERE (ZRRS(:,:,:)<=XRTMIN(3))
- ZW(:,:,:)=0.
-END WHERE
-!
-ZWSED(:,:,IKE+1) = 0.
+IF (PRESENT(PCRT)) THEN !two moments
+ WHERE (PRRT(:,:,:) > KRTMIN(3) .AND. PCRT(:,:,:) > XCTMIN(3))
+ ZW(:,:,:) = 0.
+ END WHERE
+ELSE ! one moment
+ WHERE (PRRT(:,:,:) <= KRTMIN(3))
+ ZW(:,:,:) = 0.
+ END WHERE
+END IF
+
+GRAIN(:,:,:) = .FALSE.
+
+IF (PRESENT(PCRT)) THEN ! case KHKO, C2R2, C3R5, LIMA (2-moment schemes)
+ GRAIN(:,:,:) = PRRT(:,:,:) > KRTMIN(3) .AND. PCRT(:,:,:) > XCTMIN(3)
+ELSE ! Case ICE3, REVE, KESS, ... (1-moment schemes)
+ GRAIN(:,:,:) = PRRT(:,:,:) > KRTMIN(3)
+END IF
+
+IRAIN = COUNTJV( GRAIN(:,:,:),IR1(:),IR2(:),IR3(:))
+
+IF( IRAIN >= 1 ) THEN
+DO JN = 1 , KSPLITR
+ IF( JN==1 ) THEN
+ DO JKAQ = 1,KMODE
+ DO JK = 1, IKE
+ PSVT(:,:,JK,KMODE*2+JKAQ) = PSVT(:,:,JK,KMODE*2+JKAQ) / FLOAT(KSPLITR)
+ END DO
+ END DO
+ END IF
+ IF ( IRAIN .GT. ILENALLOCR ) THEN
+ IF ( ILENALLOCR .GT. 0 ) THEN
+ DEALLOCATE (ILISTR,ZSVT,ZRHODREF,ZCRT,ZRRT,ZLAMBDA)
+ END IF
+ ILENALLOCR = MAX (IOLDALLOCR, 2*IRAIN )
+ IOLDALLOCR = ILENALLOCR
+ ALLOCATE(ILISTR(ILENALLOCR), ZRHODREF(ILENALLOCR), ZSVT(ILENALLOCR,3*KMODE),&
+ ZCRT(ILENALLOCR), ZRRT(ILENALLOCR), ZLAMBDA(ILENALLOCR))
+ END IF
+
+ DO JL = 1, IRAIN
+ DO JKAQ = 1, KMODE
+ ZSVT(JL,KMODE*2+JKAQ) = PSVT(IR1(JL),IR2(JL),IR3(JL),KMODE*2+JKAQ)
+ END DO
+!
+ IF (PRESENT(PCRT)) ZCRT(JL) = PCRT(IR1(JL),IR2(JL),IR2(JL))
+ ZRRT(JL) = PRRT(IR1(JL),IR2(JL),IR3(JL))
+ ZRHODREF(JL) = PRHODREF(IR1(JL),IR2(JL),IR3(JL))
+ ENDDO
+
+ ILISTLENR = 0
+ DO JL=1,IRAIN
+ IF (PRESENT(PCRT)) THEN !two moments
+ IF (ZRRT(JL) > KRTMIN(3) .AND. ZCRT(JL) > XCTMIN(3)) THEN
+ ILISTLENR = ILISTLENR + 1
+ ILISTR(ILISTLENR) = JL
+ END IF
+ ELSE ! one moment
+ IF (ZRRT(JL) > KRTMIN(3)) THEN
+ ILISTLENR = ILISTLENR + 1
+ ILISTR(ILISTLENR) = JL
+ END IF
+ END IF
+ END DO
+
+!
! Flux mass aerosol in rain droplets =
! Flux mass rain water * Mass aerosol in rain / Mass rain water
-DO JKAQ = 1,KMODE
+ DO JKAQ = 1,KMODE
+ DO JJ = 1, ILISTLENR
+ JL = ILISTR(JJ)
+ IF (PRESENT(PCRT)) THEN !two moments
+ IF (ZRRT(JL) > KRTMIN(3) .AND. ZCRT(JL) > XCTMIN(3)) THEN
+ ZLAMBDA(JL) = (KLBR * ZCRT(JL) / ZRRT(JL))**KLBEXR
- DO JN = 1 , KSPLITR
- ZWSED(:,:,1:IKE) = XFSEDR &
- * (ZRRS(:,:,:))**(XEXSEDR-1.) &
- * PRHODREF(:,:,:)**(XEXSEDR-XCEXVT) &
- * PSVT(:,:,:,KMODE*2+JKAQ)
- DO JK = 1, IKE
- PSVT(:,:,JK,KMODE*2+JKAQ)= PSVT(:,:,JK,KMODE*2+JKAQ) + &
- ZW(:,:,JK)*(ZWSED(:,:,JK+1)-ZWSED(:,:,JK))
- ! Aerosol mass in rain droplets need to be positive
- PSVT(:,:,JK,KMODE*2+JKAQ)= MAX(PSVT(:,:,JK,KMODE*2+JKAQ), 0.)
- END DO
- END DO
+ ZWSED(IR1(JL),IR2(JL),IR3(JL)) = XFSEDRR * ZRHODREF(JL)**(1.-KCEXVT) &
+ * ZLAMBDA(JL)**(-ZDR) &
+ * ZSVT(JL,KMODE*2+JKAQ)
+ END IF
+ ELSE ! one moments
+! cf rain_ice.f90 : l. 1077 (zwsed * psvt(kmode+2+jkaq) / zrrs)
+ IF (ZRRT(JL) > KRTMIN(3)) THEN
+
+ ZWSED(IR1(JL),IR2(JL),IR3(JL)) = XFSEDR &
+ * ZRRT(JL)**(XEXSEDR-1.) &
+ * ZRHODREF(JL)**(XEXSEDR-KCEXVT) &
+ * ZSVT(JL,KMODE*2+JKAQ)
+ END IF
+ END IF ! moments
+ END DO ! JJ
+
+ DO JK = 1, IKE
+ PSVT(:,:,JK,KMODE*2+JKAQ) = PSVT(:,:,JK,KMODE*2+JKAQ) + &
+ ZW(:,:,JK)*(ZWSED(:,:,JK+1)-ZWSED(:,:,JK))
+ END DO
+ END DO ! JKAQ
+
+END DO ! JN - time splitting
+
+ DO JKAQ = 1,KMODE
+! Aerosol mass in rain droplets need to be positive
+ PSVT(:,:,:,KMODE*2+JKAQ) = MAX(PSVT(:,:,:,KMODE*2+JKAQ), 0.)
+ END DO ! KKAQ
+END IF !(IRAIN)
+!
+IF (ALLOCATED(ILISTR)) DEALLOCATE(ILISTR)
+IF (ALLOCATED(ZSVT)) DEALLOCATE(ZSVT)
+IF (ALLOCATED(ZRHODREF)) DEALLOCATE(ZRHODREF)
+IF (ALLOCATED(ZCRT)) DEALLOCATE(ZCRT)
+IF (ALLOCATED(ZRRT)) DEALLOCATE(ZRRT)
+IF (ALLOCATED(ZLAMBDA)) DEALLOCATE(ZLAMBDA)
-END DO
!
- END SUBROUTINE AER_WET_DEP_KMT_WARM_SEDIMENT
+END SUBROUTINE AER_WET_DEP_KMT_WARM_SEDIMENT
!
!-------------------------------------------------------------------------------
!
@@ -587,127 +788,206 @@ END DO
!* 0. DECLARATIONS
!
IMPLICIT NONE
+!
!-------------------------------------------------------------------------------
!
!* 1. compute the autoconversion of r_c for r_r production: RCAUTR
!
-ZZW4(:,:,:)=0.0
+ZZW4(:,:,:) = 0.0
! to be sure no division by zero in case of ZZRCT = 0.
ZZRCT(:,:,:) = PRCT(:,:,:)
-ZZRCT(:,:,:) = MAX(ZZRCT(:,:,:), XRTMIN(2)/2.)
-
-WHERE( (ZZRCT(:,:,:)>XRTMIN(2)) .AND. (PRCS(:,:,:)>0.0 ) )
- ZZW4(:,:,:) = MIN( PRCS(:,:,:),XTIMAUTC* &
- MAX((ZZRCT(:,:,:)-XCRIAUTC/ PRHODREF(:,:,:)),0.0))
-END WHERE
+ZZRCT(:,:,:) = MAX(ZZRCT(:,:,:), KRTMIN(2)/2.)
+!
+IF (PRESENT(PCRT)) THEN ! 2-moment schemes
+!
+! from lima_warm_coal.f90 (AUTO)
+ ZLBDC3(:,:,:) = 1.E45
+ ZLBDC(:,:,:) = 1.E15
+ WHERE (ZZRCT(:,:,:) > KRTMIN(2) .AND. PCCT(:,:,:) > XCTMIN(2))
+ ZLBDC3(:,:,:) = KLBC * PCCT(:,:,:) / PRCT(:,:,:)
+ ZLBDC(:,:,:) = ZLBDC3(:,:,:)**ZLBEXC
+ END WHERE
+!
+ ZZW3(:,:,:) = 0.
+ WHERE (ZZRCT(:,:,:) > KRTMIN(2))
+ ZZW3(:,:,:) = MAX(0.0, XLAUTR*PRHODREF(:,:,:)*ZZRCT(:,:,:)* &
+ (XAUTO1/ZLBDC3(:,:,:)**4-XLAUTR_THRESHOLD)) ! L
+ ZZW4(:,:,:) = MIN(PRCS(:,:,:), MAX(0.0, XITAUTR*ZZW3(:,:,:)*ZZRCT(:,:,:)* &
+ (XAUTO2/ZLBDC3(:,:,:)-XITAUTR_THRESHOLD))) ! L/tau
+ END WHERE
+!
+ELSE ! 1-moment scheme
+!
+ WHERE ((ZZRCT(:,:,:) > KRTMIN(2)) .AND. (PRCS(:,:,:) > 0.0))
+ ZZW4(:,:,:) = MIN(PRCS(:,:,:), XTIMAUTC* &
+ MAX((ZZRCT(:,:,:)-XCRIAUTC/PRHODREF(:,:,:)), 0.0))
+ END WHERE
+!
+END IF
+!--cb--
DO JKAQ = 1,KMODE
- ZZW2(:,:,:) =0.0
- ZZW2(:,:,:)=ZZW4(:,:,:) * PSVT(:,:,:,KMODE+JKAQ)/ZZRCT(:,:,:) * PTSTEP
- ZZW2(:,:,:) = MAX(MIN( ZZW2(:,:,:), PSVT(:,:,:,KMODE+JKAQ)),0.0)
+ ZZW2(:,:,:) = 0.0
+ ZZW2(:,:,:) = ZZW4(:,:,:) * PSVT(:,:,:,KMODE+JKAQ) / ZZRCT(:,:,:) * PTSTEP
+ ZZW2(:,:,:) = MAX(MIN(ZZW2(:,:,:), PSVT(:,:,:,KMODE+JKAQ)), 0.0)
! For rain - Increase the aerosol conc in rain
- PSVT(:,:,:,KMODE*2+JKAQ) = &
- PSVT(:,:,:,KMODE*2+JKAQ) + ZZW2(:,:,:)
+ PSVT(:,:,:,KMODE*2+JKAQ) = PSVT(:,:,:,KMODE*2+JKAQ) + ZZW2(:,:,:)
! For Cloud Decrease the aerosol conc in cloud
- PSVT(:,:,:,KMODE+JKAQ) = &
- PSVT(:,:,:,KMODE+JKAQ) - ZZW2(:,:,:)
+ PSVT(:,:,:,KMODE+JKAQ) = PSVT(:,:,:,KMODE+JKAQ) - ZZW2(:,:,:)
ENDDO
-
-
+!
!
!* 2. compute the accretion of r_c for r_r production: RCACCR
!
-ZZW4(:,:,:)=0.0
-ZLBDAR(:,:,:)=0.0
-WHERE ( (ZZRCT(:,:,:)>XRTMIN(2)) .AND. (PRRT(:,:,:)>XRTMIN(3)) &
- .AND. (PRCS(:,:,:)> 0.0 ) )
- ZLBDAR(:,:,:) = XLBR*( PRHODREF(:,:,:)* PRRT(:,:,:) )**XLBEXR
- ZZW4(:,:,:) = MIN( PRCS(:,:,:),XFCACCR * ZZRCT(:,:,:) &
- * ZLBDAR(:,:,:)**XEXCACCR &
- * PRHODREF(:,:,:)**(-XCEXVT) )
-END WHERE
-!
-DO JKAQ = 1,KMODE
- ZZW2(:,:,:)=0.0
- ZZW2(:,:,:)=ZZW4(:,:,:) * PSVT(:,:,:,KMODE+JKAQ)/ZZRCT(:,:,:) * PTSTEP
- ZZW2(:,:,:) = MAX(MIN(ZZW2(:,:,:),PSVT(:,:,:,KMODE+JKAQ)), 0.0)
+ZZW4(:,:,:) = 0.0
+ZZW5(:,:,:) = 0.
+ZDIM(:,:,:) = 0.
+ZLBDAR(:,:,:)=0.
-!* 3. compute the new acquous aerosol mass
+!
+IF (PRESENT(PCRT)) THEN ! 2-moment schemes
+!
+! from lima_warm_coal.f90 (ACCR)
+ ZLBDR3(:,:,:) = 1.E30
+ ZLBDR(:,:,:) = 1.E10
+ WHERE (PRRT(:,:,:) > KRTMIN(3) .AND. PCRT(:,:,:) > XCTMIN(3))
+ ZLBDAR(:,:,:) = KLBR * (PRHODREF(:,:,:) * PRRT(:,:,:))**KLBEXR
+ ZLBDR3(:,:,:) = KLBR * PCRT(:,:,:) / PRRT(:,:,:)
+ ZLBDR(:,:,:) = ZLBDR3(:,:,:)**KLBEXR
+ ZZW4(:,:,:) = MIN(PRCS(:,:,:), XFCACCR * ZZRCT(:,:,:) &
+ * ZLBDAR(:,:,:)**XEXCACCR &
+ * PRHODREF(:,:,:)**(-KCEXVT) )
+ ZDIM(:,:,:) = XACCR1 / ZLBDAR(:,:,:)
+ END WHERE
+!
+! Accretion for D > 100 10-6 m
+ WHERE (PRRT(:,:,:) > KRTMIN(3) .AND. PCRT(:,:,:) > XCTMIN(3) .AND. &
+ ZZRCT(:,:,:) > KRTMIN(2) .AND. ZZW4(:,:,:) > 1.E-4 .AND. &
+ (PRRT(:,:,:) > 1.2*ZZW3(:,:,:)/PRHODREF(:,:,:) .OR. &
+ ZDIM(:,:,:) >= MAX(XACCR2,XACCR3/(XACCR4/ZLBDC(:,:,:)-XACCR5))))
+ ZZW5(:,:,:) = ZLBDC3(:,:,:) / ZLBDR3(:,:,:)
+ ZZW1(:,:,:) = (PCCT(:,:,:) * PCRT(:,:,:) / ZLBDC3(:,:,:)**2) * PRHODREF(:,:,:)
+ ZZW4(:,:,:) = MIN(ZZW1(:,:,:)*(XACCR_RLARGE1+XACCR_RLARGE2*ZZW5(:,:,:)), &
+ PRCS(:,:,:))
+ END WHERE
+! Accretion for D < 100 10-6 m
+ WHERE (PRRT(:,:,:) > KRTMIN(3) .AND. PCRT(:,:,:) > XCTMIN(3) .AND. &
+ ZZRCT(:,:,:) > KRTMIN(2) .AND. ZZW4(:,:,:) <= 1.E-4 .AND. &
+ (PRRT(:,:,:) > (1.2*ZZW2(:,:,:)/PRHODREF(:,:,:)) .OR. &
+ ZDIM(:,:,:) >= MAX(XACCR2,XACCR3/(XACCR4/ZLBDC(:,:,:)-XACCR5))))
+ ZZW5(:,:,:) = (ZLBDC3(:,:,:) / ZLBDR3(:,:,:))**2
+ ZZW1(:,:,:) = (PCCT(:,:,:) * PCRT(:,:,:) / ZLBDC3(:,:,:)**3) * PRHODREF(:,:,:)
+ ZZW4(:,:,:) = MIN(ZZW1(:,:,:)*(XACCR_RSMALL1+XACCR_RSMALL2*ZZW5(:,:,:)), &
+ PRCS(:,:,:))
+ END WHERE
+!
+ELSE ! 1-moment schemes
+!
+ ZLBDR(:,:,:) = 0.0
+ WHERE ((ZZRCT(:,:,:) > KRTMIN(2)) .AND. (PRRT(:,:,:) > KRTMIN(3)) &
+ .AND. (PRCS(:,:,:) > 0.0))
+ ZLBDR(:,:,:) = KLBR * (PRHODREF(:,:,:) * PRRT(:,:,:))**KLBEXR
+ ZZW4(:,:,:) = MIN(PRCS(:,:,:), XFCACCR * ZZRCT(:,:,:) &
+ * ZLBDR(:,:,:)**XEXCACCR &
+ * PRHODREF(:,:,:)**(-KCEXVT) )
+ END WHERE
+END IF
+!--cb--
+!
+DO JKAQ = 1, KMODE
+ ZZW2(:,:,:) = 0.0
+ ZZW2(:,:,:) = ZZW4(:,:,:) * PSVT(:,:,:,KMODE+JKAQ) / ZZRCT(:,:,:) * PTSTEP
+ ZZW2(:,:,:) = MAX(MIN(ZZW2(:,:,:), PSVT(:,:,:,KMODE+JKAQ)), 0.0)
+!
+!
+!* 3. compute the new acqueous aerosol mass
!
! For rain - Increase the aerosol conc in rain
- PSVT(:,:,:,KMODE*2+JKAQ) = PSVT(:,:,:,KMODE*2+JKAQ) + ZZW2(:,:,:)
+ PSVT(:,:,:,KMODE*2+JKAQ) = PSVT(:,:,:,KMODE*2+JKAQ) + ZZW2(:,:,:)
! For Cloud Decrease the aerosol conc in cloud
- PSVT(:,:,:,KMODE+JKAQ) = PSVT(:,:,:,KMODE+JKAQ) - ZZW2(:,:,:)
+ PSVT(:,:,:,KMODE+JKAQ) = PSVT(:,:,:,KMODE+JKAQ) - ZZW2(:,:,:)
ENDDO
-
- END SUBROUTINE AER_WET_DEP_KMT_ICE_WARM
+!
+END SUBROUTINE AER_WET_DEP_KMT_ICE_WARM
+!
!---------------------------------------------------------------------------------------
+!
SUBROUTINE AER_WET_DEP_KMT_EVAP
!
!* COMPUTES THE EVAPORATION OF CLOUD-RAIN FOR THE
!* RE-RELEASE OF AER INTO THE ENVIRONMENT
! --------------------------------------
!
-
!
!* 0. DECLARATIONS
! ------------
!
IMPLICIT NONE
+!
!* declaration of local variables
!
INTEGER :: JKAQ ! counter for aerosols
-
+!
+!-------------------------------------------------------------------------------
+!
!* 1. compute the evaporation of r_r: RREVAV
-
+!
!When partial reevaporation of precip takes place, the fraction of
!tracer precipitating form above is reevaporated is equal to
!half of the evaporation rate of water
!
! Rain water evaporated during PTSTEP in kg/kg
ZZEVAP(:,:,:) = PEVAP3D(:,:,:) * PTSTEP
+!
! Fraction of rain water evaporated
! at this stage (bulk), we consider that the flux of evaporated aerosol
! is a ratio of the evaporated rain water.
! It will interested to calculate with a two moment scheme (C2R2 or C3R5)
! the complete evaporation of rain droplet to use it for the compuation
! of the evaporated aerosol flux.
-ZWEVAP(:,:,:)=0.0
-WHERE( PRRT(:,:,:) .GT. XRTMIN(3) )
- ZWEVAP(:,:,:) = ZZEVAP(:,:,:)/(PRRT(:,:,:))
+ZWEVAP(:,:,:) = 0.0
+WHERE(PRRT(:,:,:) .GT. KRTMIN(3))
+ ZWEVAP(:,:,:) = ZZEVAP(:,:,:) / (PRRT(:,:,:))
END WHERE
-ZWEVAP(:,:,:)=MIN(ZWEVAP(:,:,:),1.0)
-ZWEVAP(:,:,:)=MAX(ZWEVAP(:,:,:),0.0)
-
+ZWEVAP(:,:,:) = MIN(ZWEVAP(:,:,:), 1.0)
+ZWEVAP(:,:,:) = MAX(ZWEVAP(:,:,:), 0.0)
+!
+!
!* 2. compute the mask of r_c evaporation : all cloud is evaporated
! no partial cloud evaporation at this stage
+!
ZMASK(:,:,:) = 0.
-WHERE( PRCS(:,:,:) .LT. XRTMIN(2) )
- ZMASK(:,:,:) = 1.
+WHERE(PRCS(:,:,:) .LT. KRTMIN(2))
+ ZMASK(:,:,:) = 1.
END WHERE
!
+DO JKAQ = 1, KMODE
+ ZZW1(:,:,:) = ZMASK(:,:,:) * PSVT(:,:,:,KMODE+JKAQ)
+ ZZW2(:,:,:) = ZWEVAP(:,:,:) * PSVT(:,:,:,KMODE*2+JKAQ)
+!
+ ZZW1(:,:,:) = MIN(ZZW1(:,:,:),PSVT(:,:,:,KMODE+JKAQ))
+ ZZW2(:,:,:) = MIN(ZZW2(:,:,:),PSVT(:,:,:,KMODE*2+JKAQ))
!
-DO JKAQ = 1,KMODE
- ZZW1(:,:,:) = ZMASK(:,:,:)*PSVT(:,:,:,KMODE+JKAQ)
-
- ZZW2(:,:,:) = ZWEVAP(:,:,:)*PSVT(:,:,:,KMODE*2+JKAQ)
-
! 3. New dry aerosol mass
!
- PSVT(:,:,:,JKAQ) = PSVT(:,:,:,JKAQ) + ZZW2(:,:,:) + ZZW1(:,:,:)
+ PSVT(:,:,:,JKAQ) = PSVT(:,:,:,JKAQ) + ZZW2(:,:,:) + ZZW1(:,:,:)
+!
+!
! 4. New cloud aerosol mass
!
- PSVT(:,:,:,KMODE+JKAQ) = PSVT(:,:,:,KMODE+JKAQ) - ZZW1(:,:,:)
-
+ PSVT(:,:,:,KMODE+JKAQ) = PSVT(:,:,:,KMODE+JKAQ) - ZZW1(:,:,:)
+!
+!
! 5. New rain aerosol mass
!
- PSVT(:,:,:,KMODE*2+JKAQ) = PSVT(:,:,:,KMODE*2+JKAQ) - ZZW2(:,:,:)
+ PSVT(:,:,:,KMODE*2+JKAQ) = PSVT(:,:,:,KMODE*2+JKAQ) - ZZW2(:,:,:)
END DO
!
+END SUBROUTINE AER_WET_DEP_KMT_EVAP
!
- END SUBROUTINE AER_WET_DEP_KMT_EVAP
!---------------------------------------------------------------------------------------
+!
SUBROUTINE AER_WET_DEP_KMT_EFFIC
!
!* COMPUTES THE EFFICIENCY FACTOR
@@ -719,6 +999,7 @@ END DO
!
IMPLICIT NONE
!
+!-------------------------------------------------------------------------------
!
!* 1. COMPUTES THE EFFICIENCY FACTOR
! --------------------------------------
@@ -726,48 +1007,46 @@ IMPLICIT NONE
!* 1.1 compute gravitational velocities
!
!initialize
- ZTEMP(:)=ZTHT(:)*(ZPABST(:)/XP00)**(XRD/XCPD)
- ZTEMP(:)=MAX(ZTEMP(:),1.e-12)
-
- CALL AER_VELGRAV(ZRG(:,:), ZPABST(:), KMODE, &
- ZMU(:), ZVGG(:,:), &
- ZDPG(:,:),ZTEMP(:), &
- ZCOR(:,:), ZDENSITY_AER(:,:))
+ZTEMP(:) = ZTHT(:) * (ZPABST(:) / XP00)**(XRD/XCPD)
+ZTEMP(:) = MAX(ZTEMP(:), 1.e-12)
+!
+CALL AER_VELGRAV(ZRG(:,:), ZPABST(:), KMODE, &
+ ZMU(:), ZVGG(:,:), &
+ ZDPG(:,:),ZTEMP(:), &
+ ZCOR(:,:), ZDENSITY_AER(:,:))
-! Above gives mu (ZMU), v(aerosol)(PVGG, m/s), diffusion (ZDPG, m2/s)
+! Above gives mu (ZMU), v(aerosol)(PVGG, m/s), diffusion (ZDPG, m2/s)
!
!* 1.2 Compute Water Viscocity in kg/m/s Prup. & Klett, p.95
!
-!
- A0=1.76
- A1=-5.5721e-2
- A2=-1.3943e-3
- A3=-4.3015e-5
- ZMUW(:)=A0*EXP(A1*(ZTEMP(:)-273.15) &
- +A2*(ZTEMP(:)-273.15) + A3*(ZTEMP(:)-273.15))*1.e-3
-
- A1=-3.5254e-2
- A2=4.7163e-4
- A3=-6.0667e-6
-
- WHERE (ZTEMP(:)>273.15)
- ZMUW(:)=A0*EXP(A1*(ZTEMP(:)-273.15) &
- +A2*(ZTEMP(:)-273.15) + A3*(ZTEMP(:)-273.15))*1.e-3
- END WHERE
- ZMUW(:)=MAX(ZMUW(:),1.e-12)
-
+A0 = 1.76
+A1 = -5.5721e-2
+A2 = -1.3943e-3
+A3 = -4.3015e-5
+ZMUW(:) = A0 * EXP(A1*(ZTEMP(:)-273.15) &
+ + A2*(ZTEMP(:)-273.15) + A3*(ZTEMP(:)-273.15)) * 1.e-3
+!
+A1 = -3.5254e-2
+A2 = 4.7163e-4
+A3 = -6.0667e-6
+WHERE (ZTEMP(:) > 273.15)
+ ZMUW(:) = A0 * EXP(A1*(ZTEMP(:)-273.15) &
+ + A2*(ZTEMP(:)-273.15) + A3*(ZTEMP(:)-273.15)) * 1.e-3
+END WHERE
+ZMUW(:) = MAX(ZMUW(:), 1.e-12)
!
!* 1.3 compute efficiency factor
!
! This gives aerosol collection efficiency by calculating Reynolds number
! schmidt number, stokes number, etc
- CALL AER_EFFIC(ZRG(:,:), ZVGG(:,:), & !aerosol radius/velocity
- ZRHODREF(:), & !Air density
- ZMUW(:), ZMU(:), & !mu water/air
- ZDPG(:,:), ZEFC(:,:), & !diffusivity, efficiency
- ZRRT(:), KMODE, & !Rain water, nb aerosols modes
- ZTEMP(:),ZCOR(:,:), & ! Temperature, Cunnimgham coeff
- ZDENSITY_AER(:,:)) ! aerosol density
+CALL AER_EFFIC(ZRG(:,:), ZVGG(:,:), & !aerosol radius/velocity
+ ZRHODREF(:), & !Air density
+ ZMUW(:), ZMU(:), & !mu water/air
+ ZDPG(:,:), ZEFC(:,:), & !diffusivity, efficiency
+ ZRRT(:), KMODE, & !Rain water, nb aerosols modes
+ ZTEMP(:),ZCOR(:,:), & ! Temperature, Cunnimgham coeff
+ ZDENSITY_AER(:,:), & ! aerosol density
+ ZRR, ZNT ) ! radius and number of rain drops
!
END SUBROUTINE AER_WET_DEP_KMT_EFFIC
!
diff --git a/src/MNH/aeroopt_get.f90 b/src/MNH/aeroopt_get.f90
index 485e87b68..d8ec58193 100644
--- a/src/MNH/aeroopt_get.f90
+++ b/src/MNH/aeroopt_get.f90
@@ -96,10 +96,10 @@
!LOCAL VARIABLES
INTEGER :: NMODE_AER
REAL, DIMENSION(SIZE(PSVTA,1),SIZE(PSVTA,2),SIZE(PSVTA,3),SIZE(PSVTA,4)) :: ZSVT
- REAL, DIMENSION(SIZE(PSVTA,1),SIZE(PSVTA,2),SIZE(PSVTA,3),16 ,2) :: ZMASS ![kg/m3] mass of aerosol mode
- REAL, DIMENSION(SIZE(PSVTA,1),SIZE(PSVTA,2),SIZE(PSVTA,3) ,2) :: ZMASSeq ![kg/m3] mass of aerosol mode
- REAL, DIMENSION(SIZE(PSVTA,1),SIZE(PSVTA,2),SIZE(PSVTA,3), 2) :: ZRADIUS ![um] number median radius ofaerosol mode
- REAL, DIMENSION(SIZE(PSVTA,1),SIZE(PSVTA,2),SIZE(PSVTA,3), 2) :: ZSIGMA ![-] dispersion coefficientaerosol mode
+ REAL, DIMENSION(SIZE(PSVTA,1),SIZE(PSVTA,2),SIZE(PSVTA,3),NSP+NCARB+NSOA ,JPMODE) :: ZMASS ![kg/m3] mass of aerosol mode
+ REAL, DIMENSION(SIZE(PSVTA,1),SIZE(PSVTA,2),SIZE(PSVTA,3) , JPMODE) :: ZMASSeq ![kg/m3] mass of aerosol mode
+ REAL, DIMENSION(SIZE(PSVTA,1),SIZE(PSVTA,2),SIZE(PSVTA,3), JPMODE) :: ZRADIUS ![um] number median radius ofaerosol mode
+ REAL, DIMENSION(SIZE(PSVTA,1),SIZE(PSVTA,2),SIZE(PSVTA,3), JPMODE) :: ZSIGMA ![-] dispersion coefficientaerosol mode
REAL, ALLOCATABLE, DIMENSION(:,:,:,:) :: ZTAU550_MDE ![-] opt.depth 550nm one mode
REAL, ALLOCATABLE, DIMENSION(:,:,:,:,:) :: ZTAU_WVL_MDE ![-] opt.depth @ wvl, one mode
REAL, ALLOCATABLE, DIMENSION(:,:,:,:,:) :: ZPIZA_WVL_MDE ![-] single scattering albedo @ wvl, one mode
@@ -133,13 +133,13 @@
COMPLEX, DIMENSION(size(ZMASS,1),size(ZMASS,2),size(ZMASS,3),6)::Req ! Equivalent refractive index
REAL, PARAMETER :: EPSILON=1.d-8 ![um] a small number used to avoid zero
- INTEGER ::JJJ
+ INTEGER ::JJJ, JI, JSV
!-------------------------------------------------------------------------
!-------------------------------------------------------------------------
!-------------------------------------------------------------------------
!
- NMODE_AER=2
+NMODE_AER=JPMODE ! in case of ORILAM
!Allocate arrays which size depend on number of modes
@@ -148,6 +148,7 @@
ALLOCATE(ZPIZA_WVL_MDE(SIZE(PTAU550,1),SIZE(PTAU550,2),SIZE(PTAU550,3),KSWB,NMODE_AER))
ALLOCATE(ZCGA_WVL_MDE(SIZE(PTAU550,1),SIZE(PTAU550,2),SIZE(PTAU550,3),KSWB,NMODE_AER))
ZSVT(:,:,:,:)=PSVTA(:,:,:,:)
+
CALL PPP2AERO( &
ZSVT & !I [moments/molec_{air}] moments of aerosol for all modes
,PRHODREFA & !I [kg/m3] air density
@@ -157,8 +158,8 @@
)
- ZMASS(:,:,:,:,:)=ZMASS(:,:,:,:,:)*1.E-9
-
+ ZMASS(:,:,:,:,:)=MAX(ZMASS(:,:,:,:,:)*1.E-9, 1E-20)
+
DO JMDE=1,NMODE_AER
Ri(1,1)=CMPLX(1.80,-7.40E-1,kind=kind(Ri(1,1)))
@@ -207,30 +208,35 @@
! Computation of the refractive index for the whole aerosol mode according to
! Maxwell-Garnett mixing rule
! IF (LDUST) THEN
-! VDDST(:,:,:)=(ZMASS(:,:,:,JP_AER_DST,JMDE))/XFAC(JP_AER_DST)
! ELSE
- VDDST(:,:,:)=0.
+! VDDST(:,:,:)=0.
! ENDIF
- VOC(:,:,:)=(ZMASS(:,:,:,5,JMDE))/XFAC(5)
- VH2O(:,:,:)=(ZMASS(:,:,:,4,JMDE))/XFAC(4)
- VAM(:,:,:)=(ZMASS(:,:,:,3,JMDE))/XFAC(3)
- VSU(:,:,:)=(ZMASS(:,:,:,1,JMDE))/XFAC(1)
- VNI(:,:,:)=(ZMASS(:,:,:,2,JMDE))/XFAC(2)
- VBC(:,:,:)=(ZMASS(:,:,:,6,JMDE))/XFAC(6)
- VSOA1(:,:,:)=(ZMASS(:,:,:,7,JMDE))/XFAC(7)
- VSOA2(:,:,:)=(ZMASS(:,:,:,8,JMDE))/XFAC(8)
- VSOA3(:,:,:)=(ZMASS(:,:,:,9,JMDE))/XFAC(9)
- VSOA4(:,:,:)=(ZMASS(:,:,:,10,JMDE))/XFAC(10)
- VSOA5(:,:,:)=(ZMASS(:,:,:,11,JMDE))/XFAC(11)
- VSOA6(:,:,:)=(ZMASS(:,:,:,12,JMDE))/XFAC(12)
- VSOA7(:,:,:)=(ZMASS(:,:,:,13,JMDE))/XFAC(13)
- VSOA8(:,:,:)=(ZMASS(:,:,:,14,JMDE))/XFAC(14)
- VSOA9(:,:,:)=(ZMASS(:,:,:,15,JMDE))/XFAC(15)
- VSOA10(:,:,:)=(ZMASS(:,:,:,16,JMDE))/XFAC(16)
+
+ VOC(:,:,:)=(ZMASS(:,:,:,JP_AER_OC,JMDE))/XFAC(JP_AER_OC)
+ VH2O(:,:,:)=(ZMASS(:,:,:,JP_AER_H2O,JMDE))/XFAC(JP_AER_H2O)
+ VAM(:,:,:)=(ZMASS(:,:,:,JP_AER_NH3,JMDE))/XFAC(JP_AER_NH3)
+ VSU(:,:,:)=(ZMASS(:,:,:,JP_AER_SO4,JMDE))/XFAC(JP_AER_SO4)
+ VNI(:,:,:)=(ZMASS(:,:,:,JP_AER_NO3,JMDE))/XFAC(JP_AER_NO3)
+ VBC(:,:,:)=(ZMASS(:,:,:,JP_AER_BC,JMDE))/XFAC(JP_AER_BC)
+ VDDST(:,:,:)=(ZMASS(:,:,:,JP_AER_DST,JMDE))/XFAC(JP_AER_DST)
+ IF (NSOA .EQ. 10) THEN
+ VSOA1(:,:,:)=(ZMASS(:,:,:,JP_AER_SOA1,JMDE))/XFAC(JP_AER_SOA1)
+ VSOA2(:,:,:)=(ZMASS(:,:,:,JP_AER_SOA2,JMDE))/XFAC(JP_AER_SOA2)
+ VSOA3(:,:,:)=(ZMASS(:,:,:,JP_AER_SOA3,JMDE))/XFAC(JP_AER_SOA3)
+ VSOA4(:,:,:)=(ZMASS(:,:,:,JP_AER_SOA4,JMDE))/XFAC(JP_AER_SOA4)
+ VSOA5(:,:,:)=(ZMASS(:,:,:,JP_AER_SOA5,JMDE))/XFAC(JP_AER_SOA5)
+ VSOA6(:,:,:)=(ZMASS(:,:,:,JP_AER_SOA6,JMDE))/XFAC(JP_AER_SOA6)
+ VSOA7(:,:,:)=(ZMASS(:,:,:,JP_AER_SOA7,JMDE))/XFAC(JP_AER_SOA7)
+ VSOA8(:,:,:)=(ZMASS(:,:,:,JP_AER_SOA8,JMDE))/XFAC(JP_AER_SOA8)
+ VSOA9(:,:,:)=(ZMASS(:,:,:,JP_AER_SOA9,JMDE))/XFAC(JP_AER_SOA9)
+ VSOA10(:,:,:)=(ZMASS(:,:,:,JP_AER_SOA10,JMDE))/XFAC(JP_AER_SOA10)
VSOA(:,:,:)=VSOA1(:,:,:)+VSOA2(:,:,:)+VSOA3(:,:,:)+VSOA4(:,:,:)+&
VSOA5(:,:,:)+VSOA6(:,:,:)+VSOA7(:,:,:)+VSOA8(:,:,:)+&
VSOA9(:,:,:)+VSOA10(:,:,:)
+ ELSE
+ VSOA(:,:,:)=0.
+ END IF
VEXTR(:,:,:)=VSOA(:,:,:)+VH2O(:,:,:)+VAM(:,:,:)+VSU(:,:,:)+VNI(:,:,:)
@@ -251,11 +257,11 @@
ENDWHERE
ENDDO
- ZMASSeq(:,:,:,JMDE)=ZMASS(:,:,:,1,JMDE)+ZMASS(:,:,:,2,JMDE)+ZMASS(:,:,:,3,JMDE)&
- +ZMASS(:,:,:,4,JMDE)+ZMASS(:,:,:,5,JMDE)+ZMASS(:,:,:,6,JMDE)+ZMASS(:,:,:,7,JMDE)&
- +ZMASS(:,:,:,8,JMDE)+ZMASS(:,:,:,9,JMDE)+ZMASS(:,:,:,10,JMDE)+ZMASS(:,:,:,11,JMDE)&
- +ZMASS(:,:,:,12,JMDE)+ZMASS(:,:,:,13,JMDE)+ZMASS(:,:,:,14,JMDE)+ZMASS(:,:,:,15,JMDE)&
- +ZMASS(:,:,:,16,JMDE)
+ ZMASSeq(:,:,:,JMDE)=0.
+ DO JSV=1,NSP+NCARB+NSOA
+ ZMASSeq(:,:,:,JMDE)=ZMASSeq(:,:,:,JMDE)+ ZMASS(:,:,:,JSV,JMDE)
+ END DO
+
PII(:,:,:,:) = aimag(CMPLX(Req(:,:,:,:),kind=kind(PII(1,1,1,1))))
PIR(:,:,:,:) = real( CMPLX(Req(:,:,:,:),kind=kind(PIR(1,1,1,1))))
!Get aerosol optical properties from look up tables
diff --git a/src/MNH/ch_aer_depos.f90 b/src/MNH/ch_aer_depos.f90
index 0e29a73db..cf5659aa3 100644
--- a/src/MNH/ch_aer_depos.f90
+++ b/src/MNH/ch_aer_depos.f90
@@ -63,7 +63,7 @@ USE MODE_AERO_PSD
USE MODD_CST, ONLY: XMD, &! Molar mass of dry air
XPI
USE MODD_NSV, ONLY : NSV_C2R2BEG, NSV_AERBEG,NSV_AEREND,NSV_AER, &
- NSV_AERDEP, NSV_AERDEPBEG
+ NSV_AERDEP, NSV_AERDEPBEG, NSV_LIMA_NC, NSV_LIMA_NR
USE MODD_PARAM_n, ONLY: CCLOUD
USE MODI_AER_WET_DEP_KMT_WARM
@@ -156,7 +156,7 @@ DO JN=1,JPMODE*2
ENDDO
DO JN=1,JPMODE
- ZSVAER(:,:,:,JN) = ZPM(:,:,:,JN*3-1) * XMD * XPI * &
+ ZSVAER(:,:,:,JN) = ZPM(:,:,:,NM3(JN)) * XMD * XPI * &
4./3. * ZRHOP(:,:,:,JN) / &
(ZMI(:,:,:,JN)*1.d15*PRHODREF(:,:,:))
ENDDO
@@ -180,6 +180,15 @@ SELECT CASE (CCLOUD)
PSEA=ZSEA, PTOWN=ZTOWN, &
PCCT=PSVT(:,:,:,NSV_C2R2BEG+1), &
PCRT=PSVT(:,:,:,NSV_C2R2BEG+2) )
+ CASE ('LIMA')
+! Two moment cloud scheme
+ CALL AER_WET_DEP_KMT_WARM (NSPLITR, PTSTEP, PZZ, PRHODREF, &
+ PRT(:,:,:,2), PRT(:,:,:,3), ZRCS, &
+ ZRRS, ZSVAER, PTHT, PPABST, ZRG, &
+ PEVAP3D, JPMODE,ZRHOP, ZMASSMIN, &
+ PSEA=ZSEA, PTOWN=ZTOWN, &
+ PCCT=PSVT(:,:,:,NSV_LIMA_NC), &
+ PCRT=PSVT(:,:,:,NSV_LIMA_NR) )
END SELECT
@@ -190,32 +199,38 @@ ENDDO
DO JN=1,JPMODE
! Return from dry mass on m3
- ZPM(:,:,:,JN*3-1) = ZSVAER(:,:,:,JN) * &
+ ! Compute m0 and m6 using m3, Rg and sigma
+
+ !ZPM(:,:,:,JN*3-1) = ZSVAER(:,:,:,JN) * &
+ ZPM(:,:,:,NM3(JN)) = ZSVAER(:,:,:,JN) * &
ZMI(:,:,:,JN)*1.d15*PRHODREF(:,:,:) / &
(XMD * XPI * 4./3. * ZRHOP(:,:,:,JN) )
-
- ! Compute m0 and m6 using m3, Rg and sigma
- ZPM(:,:,:,JN*3-2) = ZPM(:,:,:,JN*3-1)/ &
+ !ZPM(:,:,:,JN*3-2) = ZPM(:,:,:,JN*3-1)/ &
+ ZPM(:,:,:,NM0(JN)) = ZPM(:,:,:,NM3(JN)) / &
( (ZRG(:,:,:,JN)**3)*EXP(4.5 * LOG(ZSIG(:,:,:,JN))**2) )
- ZPM(:,:,:,JN*3) = ZPM(:,:,:,JN*3-2)*(ZRG(:,:,:,JN)**6) * &
+ !ZPM(:,:,:,JN*3) = ZPM(:,:,:,JN*3-2)*(ZRG(:,:,:,JN)**6) * &
+ ZPM(:,:,:,NM6(JN)) = ZPM(:,:,:,NM0(JN))*(ZRG(:,:,:,JN)**6) * &
EXP(18 *(LOG(ZSIG(:,:,:,JN)))**2)
+
ENDDO
! Filter minimum values
-!DO JN=1,JPMODE
-! WHERE ((ZPM(:,:,:,NM0(JN)) .LE. ZPMIN(NM0(JN))).OR.&
-! (ZPM(:,:,:,NM3(JN)) .LE. ZPMIN(NM3(JN))).OR.&
-! (ZPM(:,:,:,NM6(JN)) .LE. ZPMIN(NM6(JN))))
-! ZPM(:,:,:,NM0(JN)) = ZPMOLD(:,:,:,NM0(JN))
-! ZPM(:,:,:,NM3(JN)) = ZPMOLD(:,:,:,NM3(JN))
-! ZPM(:,:,:,NM6(JN)) = ZPMOLD(:,:,:,NM6(JN))
-! END WHERE
-!ENDDO
+DO JN=1,JPMODE
+ WHERE ((ZPM(:,:,:,NM0(JN)) .LE. ZPMIN(NM0(JN))).OR.&
+ (ZPM(:,:,:,NM3(JN)) .LE. ZPMIN(NM3(JN))).OR.&
+ (ZPM(:,:,:,NM6(JN)) .LE. ZPMIN(NM6(JN))))
+ ZPM(:,:,:,NM0(JN)) = ZPMOLD(:,:,:,NM0(JN))
+ ZPM(:,:,:,NM3(JN)) = ZPMOLD(:,:,:,NM3(JN))
+ ZPM(:,:,:,NM6(JN)) = ZPMOLD(:,:,:,NM6(JN))
+ END WHERE
+ENDDO
! Compute final sedimentation tendency with adding acqueous sedimentation
DO JN=1,JPIN
PSEDA(:,:,:,JN) = PSEDA(:,:,:,JN) + (ZPM(:,:,:,JN) - ZPMOLD(:,:,:,JN)) / PTSTEP
END DO
+
!
+
END SUBROUTINE CH_AER_DEPOS
--
GitLab