From 99e457cc1b04afd0a1abc54a0b78a5a9412ad054 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Beno=C3=AEt=20Vi=C3=A9?= <benoit.vie@meteo.fr>
Date: Thu, 24 Mar 2022 15:38:38 +0100
Subject: [PATCH] LIMA_KHKO compilation bugfix
---
src/MNH/lima_warm_coal.f90 | 109 ++++++++++++++++++++++++-------------
src/MNH/lima_warm_evap.f90 | 81 ++++++++++++++++++++-------
2 files changed, 131 insertions(+), 59 deletions(-)
diff --git a/src/MNH/lima_warm_coal.f90 b/src/MNH/lima_warm_coal.f90
index 4ec69ac58..23e991643 100644
--- a/src/MNH/lima_warm_coal.f90
+++ b/src/MNH/lima_warm_coal.f90
@@ -104,6 +104,7 @@ END MODULE MODI_LIMA_WARM_COAL
! ------------
!
use modd_budget, only: lbudget_rc, lbudget_rr, lbudget_sv, NBUDGET_RC, NBUDGET_RR, NBUDGET_SV1, tbudgets
+USE MODD_CST, ONLY: XPI, XRHOLW
USE MODD_NSV, ONLY: NSV_LIMA_NC, NSV_LIMA_NR
USE MODD_PARAMETERS, ONLY: JPHEXT, JPVEXT
USE MODD_PARAM_LIMA
@@ -252,7 +253,7 @@ IF (LRAIN) THEN
GSELF(:) = ZCCT(:)>XCTMIN(2)
ISELF = COUNT(GSELF(:))
- IF( ISELF>0 ) THEN
+ IF( ISELF>0 .AND. .NOT.LKHKO) THEN
ZZW1(:) = XSELFC*(ZCCT(:)/ZLBDC3(:))**2 * ZRHODREF(:) ! analytical integration
WHERE( GSELF(:) )
ZCCS(:) = ZCCS(:) - MIN( ZCCS(:),ZZW1(:) )
@@ -273,38 +274,54 @@ IF (LRAIN) THEN
if ( lbudget_rc ) call Budget_store_init( tbudgets(NBUDGET_RC), 'AUTO', prcs(:, :, :) * prhodj(:, :, :) )
if ( lbudget_rr ) call Budget_store_init( tbudgets(NBUDGET_RR), 'AUTO', prrs(:, :, :) * prhodj(:, :, :) )
if ( lbudget_sv ) then
- !call Budget_store_init( tbudgets(NBUDGET_SV1 - 1 + nsv_lima_nc), 'AUTO', pccs(:, :, :) * prhodj(:, :, :) )
+ call Budget_store_init( tbudgets(NBUDGET_SV1 - 1 + nsv_lima_nc), 'AUTO', pccs(:, :, :) * prhodj(:, :, :) )
call Budget_store_init( tbudgets(NBUDGET_SV1 - 1 + nsv_lima_nr), 'AUTO', pcrs(:, :, :) * prhodj(:, :, :) )
end if
ZZW2(:) = 0.0
ZZW1(:) = 0.0
- WHERE( ZRCT(:)>XRTMIN(2) )
- ZZW2(:) = MAX( 0.0,XLAUTR*ZRHODREF(:)*ZRCT(:)* &
+ IF (LKHKO) THEN
+ WHERE ( ZRCT(:) .GT. XRTMIN(2) .AND. ZCCT(:) .GT. XCTMIN(2) &
+ .AND. (ZRCS(:) .GT. 0.0) .AND. (ZCCS(:) .GT. 0.0))
+!
+ ZZW1(:)= 1350.0 * ZRCT(:)**(2.47) * (ZCCT(:)/1.0E6)**(-1.79) ! ZCCT in cm-3
+ ZZW1(:) = min (ZRCS(:), ZZW1(:))
+ ZRCS(:) = ZRCS(:) - ZZW1(:)
+ ZRRS(:) = ZRRS(:) + ZZW1(:)
+!
+ ZCRS(:) = ZCRS(:) + ZZW1(:) * 3. * ZRHODREF(:)/(4.*XPI*XRHOLW*(XR0)**(3.))
+!
+ ZZW1(:) = min ( ZCCS(:),ZZW1(:) * ZCCT(:) / ZRCT(:))
+ ZCCS(:) = ZCCS(:) - ZZW1(:)
+!
+ END WHERE
+ ELSE
+ WHERE( ZRCT(:)>XRTMIN(2) )
+ ZZW2(:) = MAX( 0.0,XLAUTR*ZRHODREF(:)*ZRCT(:)* &
(XAUTO1/min(ZLBDC(:),1.e9)**4-XLAUTR_THRESHOLD) ) ! L
!
- ZZW3(:) = MIN( ZRCS(:), MAX( 0.0,XITAUTR*ZZW2(:)*ZRCT(:)* &
+ ZZW3(:) = MIN( ZRCS(:), MAX( 0.0,XITAUTR*ZZW2(:)*ZRCT(:)* &
(XAUTO2/ZLBDC(:)-XITAUTR_THRESHOLD) ) ) ! L/tau
!
- ZRCS(:) = ZRCS(:) - ZZW3(:)
- ZRRS(:) = ZRRS(:) + ZZW3(:)
+ ZRCS(:) = ZRCS(:) - ZZW3(:)
+ ZRRS(:) = ZRRS(:) + ZZW3(:)
!
- ZZW1(:) = MIN( MIN( 1.2E4,(XACCR4/ZLBDC(:)-XACCR5)/XACCR3), &
+ ZZW1(:) = MIN( MIN( 1.2E4,(XACCR4/ZLBDC(:)-XACCR5)/XACCR3), &
ZLBDR(:)/XACCR1 ) ! D**-1 threshold diameter for
! switching the autoconversion regimes
! min (80 microns, D_h, D_r)
- ZZW3(:) = ZZW3(:) * MAX( 0.0,ZZW1(:) )**3 / XAC
- ZCRS(:) = ZCRS(:) + ZZW3(:)
- END WHERE
-
+ ZZW3(:) = ZZW3(:) * MAX( 0.0,ZZW1(:) )**3 / XAC
+ ZCRS(:) = ZCRS(:) + ZZW3(:)
+ END WHERE
+ END IF
if ( lbudget_rc ) call Budget_store_end( tbudgets(NBUDGET_RC), 'AUTO', &
Unpack( zrcs(:), mask = gmicro(:, :, :), field = prcs(:, :, :) ) * prhodj(:, :, :) )
if ( lbudget_rr ) call Budget_store_end( tbudgets(NBUDGET_RR), 'AUTO', &
Unpack( zrrs(:), mask = gmicro(:, :, :), field = prrs(:, :, :) ) * prhodj(:, :, :) )
if ( lbudget_sv ) then
!This budget is = 0 for nsv_lima_nc => not necessary to call it (ZCCS is not modified in this part)
- !call Budget_store_end( tbudgets(NBUDGET_SV1 - 1 + nsv_lima_nc), 'AUTO', &
- ! Unpack( zccs(:), mask = gmicro(:, :, :), field = pccs(:, :, :) ) * prhodj(:, :, :) )
+ call Budget_store_end( tbudgets(NBUDGET_SV1 - 1 + nsv_lima_nc), 'AUTO', &
+ Unpack( zccs(:), mask = gmicro(:, :, :), field = pccs(:, :, :) ) * prhodj(:, :, :) )
call Budget_store_end( tbudgets(NBUDGET_SV1 - 1 + nsv_lima_nr), 'AUTO', &
Unpack( zcrs(:), mask = gmicro(:, :, :), field = pcrs(:, :, :) ) * prhodj(:, :, :) )
end if
@@ -334,31 +351,45 @@ IF (LRAIN) THEN
Unpack( zrrs(:), mask = gmicro(:, :, :), field = prrs(:, :, :) ) * prhodj(:, :, :) )
if ( lbudget_sv ) call Budget_store_init( tbudgets(NBUDGET_SV1 - 1 + nsv_lima_nc), 'ACCR', &
Unpack( zccs(:), mask = gmicro(:, :, :), field = pccs(:, :, :) ) * prhodj(:, :, :) )
- WHERE( GACCR(:).AND.(ZZW4(:)>1.E-4) ) ! Accretion for D>100 10-6 m
- ZZW3(:) = ZLBDC3(:) / ZLBDR3(:)
- ZZW1(:) = ( ZCCT(:)*ZCRT(:) / ZLBDC3(:) )*ZRHODREF(:)
- ZZW2(:) = MIN( ZZW1(:)*(XACCR_CLARGE1+XACCR_CLARGE2*ZZW3(:)),ZCCS(:) )
- ZCCS(:) = ZCCS(:) - ZZW2(:)
-!
- ZZW1(:) = ( ZZW1(:) / ZLBDC3(:) )
- ZZW2(:) = MIN( ZZW1(:)*(XACCR_RLARGE1+XACCR_RLARGE2*ZZW3(:)),ZRCS(:) )
- ZRCS(:) = ZRCS(:) - ZZW2(:)
- ZRRS(:) = ZRRS(:) + ZZW2(:)
- END WHERE
- WHERE( GACCR(:).AND.(ZZW4(:)<=1.E-4) ) ! Accretion for D<100 10-6 m
- ZZW3(:) = MIN(ZLBDC3(:) / ZLBDR3(:), 1.E8)
- ZZW1(:) = ( ZCCT(:)*ZCRT(:) / ZLBDC3(:) )*ZRHODREF(:)
- ZZW1(:) = ZZW1(:) / ZLBDC3(:)
- ZZW3(:) = ZZW3(:)**2
- ZZW2(:) = MIN( ZZW1(:)*(XACCR_CSMALL1+XACCR_CSMALL2*ZZW3(:)),ZCCS(:) )
- ZCCS(:) = ZCCS(:) - ZZW2(:)
-!
- ZZW1(:) = ZZW1(:) / ZLBDC3(:)
- ZZW2(:) = MIN( ZZW1(:)*(XACCR_RSMALL1+XACCR_RSMALL2*ZZW3(:)) &
+ IF (LKHKO) THEN
+ WHERE ( (ZRCT(:) .GT. XRTMIN(2)) .AND. (ZRRT(:) .GT. XRTMIN(3)) &
+ .AND. (ZRCS(:) .GT. 0.0) .AND. (ZCCS(:) .GT. 0.0))
+ ZZW1(:) = 67.0 * ( ZRCT(:) * ZRRT(:) )**1.15
+ ZZW1(:) = MIN (ZRCS(:),ZZW1(:))
+ ZRCS(:) = ZRCS(:) - ZZW1(:)
+ ZRRS(:) = ZRRS(:) + ZZW1(:)
+!
+ ZZW1(:) = MIN (ZCCS(:),ZZW1(:) * ZCCT(:) / ZRCT(:))
+ ZCCS(:) = ZCCS(:) - ZZW1(:)
+!
+ END WHERE
+ ELSE
+ WHERE( GACCR(:).AND.(ZZW4(:)>1.E-4) ) ! Accretion for D>100 10-6 m
+ ZZW3(:) = ZLBDC3(:) / ZLBDR3(:)
+ ZZW1(:) = ( ZCCT(:)*ZCRT(:) / ZLBDC3(:) )*ZRHODREF(:)
+ ZZW2(:) = MIN( ZZW1(:)*(XACCR_CLARGE1+XACCR_CLARGE2*ZZW3(:)),ZCCS(:) )
+ ZCCS(:) = ZCCS(:) - ZZW2(:)
+!
+ ZZW1(:) = ( ZZW1(:) / ZLBDC3(:) )
+ ZZW2(:) = MIN( ZZW1(:)*(XACCR_RLARGE1+XACCR_RLARGE2*ZZW3(:)),ZRCS(:) )
+ ZRCS(:) = ZRCS(:) - ZZW2(:)
+ ZRRS(:) = ZRRS(:) + ZZW2(:)
+ END WHERE
+ WHERE( GACCR(:).AND.(ZZW4(:)<=1.E-4) ) ! Accretion for D<100 10-6 m
+ ZZW3(:) = MIN(ZLBDC3(:) / ZLBDR3(:), 1.E8)
+ ZZW1(:) = ( ZCCT(:)*ZCRT(:) / ZLBDC3(:) )*ZRHODREF(:)
+ ZZW1(:) = ZZW1(:) / ZLBDC3(:)
+ ZZW3(:) = ZZW3(:)**2
+ ZZW2(:) = MIN( ZZW1(:)*(XACCR_CSMALL1+XACCR_CSMALL2*ZZW3(:)),ZCCS(:) )
+ ZCCS(:) = ZCCS(:) - ZZW2(:)
+!
+ ZZW1(:) = ZZW1(:) / ZLBDC3(:)
+ ZZW2(:) = MIN( ZZW1(:)*(XACCR_RSMALL1+XACCR_RSMALL2*ZZW3(:)) &
,ZRCS(:) )
- ZRCS(:) = ZRCS(:) - ZZW2(:)
- ZRRS(:) = ZRRS(:) + ZZW2(:)
- END WHERE
+ ZRCS(:) = ZRCS(:) - ZZW2(:)
+ ZRRS(:) = ZRRS(:) + ZZW2(:)
+ END WHERE
+ END IF
if ( lbudget_rc ) call Budget_store_end( tbudgets(NBUDGET_RC), 'ACCR', &
Unpack( zrcs(:), mask = gmicro(:, :, :), field = prcs(:, :, :) ) * prhodj(:, :, :) )
@@ -380,7 +411,7 @@ IF (LRAIN) THEN
ELSE
ISCBU = 0.0
END IF
- IF( ISCBU>0 ) THEN
+ IF( ISCBU>0 .AND. .NOT.LKHKO) THEN
if ( lbudget_sv ) call Budget_store_init( tbudgets(NBUDGET_SV1 - 1 + nsv_lima_nr), 'SCBU', &
Unpack( zcrs(:), mask = gmicro(:, :, :), field = pcrs(:, :, :) ) * prhodj(:, :, :) )
!
diff --git a/src/MNH/lima_warm_evap.f90 b/src/MNH/lima_warm_evap.f90
index 9a67a4b82..dd0cd055e 100644
--- a/src/MNH/lima_warm_evap.f90
+++ b/src/MNH/lima_warm_evap.f90
@@ -137,6 +137,7 @@ REAL, DIMENSION(:) , ALLOCATABLE :: ZCRT ! rain conc. at t
!
REAL, DIMENSION(:) , ALLOCATABLE :: ZRVS ! Water vapor m.r. source
REAL, DIMENSION(:) , ALLOCATABLE :: ZRRS ! Rain water m.r. source
+REAL, DIMENSION(:) , ALLOCATABLE :: ZCRS ! Rain water m.r. source
REAL, DIMENSION(:) , ALLOCATABLE :: ZTHS ! Theta source
!
! Other packed variables
@@ -151,7 +152,7 @@ REAL, DIMENSION(:), ALLOCATABLE :: ZZW1, ZZW2, ZZW3, &
ZZLV ! Latent heat of vaporization at T
!
REAL, DIMENSION(SIZE(PRHODREF,1),SIZE(PRHODREF,2),SIZE(PRHODREF,3)) &
- :: ZW, ZW2, ZRVSAT, ZDR
+ :: ZW, ZW2, ZRVSAT, ZDR, ZLV
!
!
REAL :: ZEPS, ZFACT
@@ -179,11 +180,14 @@ ZCTMIN(:) = XCTMIN(:) / PTSTEP
ZEPS= XMV / XMD
ZRVSAT(:,:,:) = ZEPS / (PPABST(:,:,:) * &
EXP(-XALPW+XBETAW/ZT(:,:,:)+XGAMW*ALOG(ZT(:,:,:))) - 1.0)
-
+ZLV(:,:,:) = XLVTT + (XCPV-XCL)*(ZT(:,:,:)-XTT)
!
GEVAP(:,:,:) = .FALSE.
GEVAP(IIB:IIE,IJB:IJE,IKB:IKE) = &
- PRRS(IIB:IIE,IJB:IJE,IKB:IKE)>ZRTMIN(3) .AND. &
+ PRRS(IIB:IIE,IJB:IJE,IKB:IKE)>ZRTMIN(3) .AND. &
+ PCRS(IIB:IIE,IJB:IJE,IKB:IKE)>ZCTMIN(3) .AND. &
+ PRRT(IIB:IIE,IJB:IJE,IKB:IKE)>ZRTMIN(3) .AND. &
+ PCRT(IIB:IIE,IJB:IJE,IKB:IKE)>ZCTMIN(3) .AND. &
PRVT(IIB:IIE,IJB:IJE,IKB:IKE)<ZRVSAT(IIB:IIE,IJB:IJE,IKB:IKE)
!
IEVAP = COUNTJV( GEVAP(:,:,:),I1(:),I2(:),I3(:))
@@ -196,6 +200,7 @@ IF( IEVAP >= 1 ) THEN
!
ALLOCATE(ZRVS(IEVAP))
ALLOCATE(ZRRS(IEVAP))
+ ALLOCATE(ZCRS(IEVAP))
ALLOCATE(ZTHS(IEVAP))
!
ALLOCATE(ZLBDR(IEVAP))
@@ -212,6 +217,7 @@ IF( IEVAP >= 1 ) THEN
ZRRT(JL) = PRRT(I1(JL),I2(JL),I3(JL))
ZCRT(JL) = PCRT(I1(JL),I2(JL),I3(JL))
ZRRS(JL) = PRRS(I1(JL),I2(JL),I3(JL))
+ ZCRS(JL) = PCRS(I1(JL),I2(JL),I3(JL))
ZRVS(JL) = PRVS(I1(JL),I2(JL),I3(JL))
ZTHS(JL) = PTHS(I1(JL),I2(JL),I3(JL))
ZZT(JL) = ZT(I1(JL),I2(JL),I3(JL))
@@ -219,8 +225,8 @@ IF( IEVAP >= 1 ) THEN
ZLBDR(JL) = ZWLBDR(I1(JL),I2(JL),I3(JL))
ZRHODREF(JL) = PRHODREF(I1(JL),I2(JL),I3(JL))
ZEXNREF(JL) = PEXNREF(I1(JL),I2(JL),I3(JL))
+ ZZLV(JL) = ZLV(I1(JL),I2(JL),I3(JL))
END DO
- ZZLV(:) = XLVTT + (XCPV-XCL)*(ZZT(:)-XTT)
!
ALLOCATE(ZZW2(IEVAP))
ALLOCATE(ZZW3(IEVAP))
@@ -242,10 +248,16 @@ IF( IEVAP >= 1 ) THEN
!
! Compute the evaporation tendency
!
- ZZW2(:) = MIN( ZZW2(:) * ZZW3(:) * ZRRT(:) * &
+ IF (LKHKO) THEN
+ ZZW2(:) = 3.0 * XCEVAP * ZZW2(:) * (4.*XPI*XRHOLW/(3.*ZRHODREF(:)))**(2./3.) * &
+ (ZRRT(:))**(1./3.) * (ZCRT(:))**(2./3.) * ZZW3(:)
+ ZZW2(:) = MIN(ZZW2(:),ZRRS(:))
+ ELSE
+ ZZW2(:) = MIN( ZZW2(:) * ZZW3(:) * ZRRT(:) * &
(X0EVAR*ZLBDR(:)**XEX0EVAR + X1EVAR*ZRHODREF(:)**XEX2EVAR* &
ZLBDR(:)**XEX1EVAR),ZRRS(:) )
- ZZW2(:) = MAX(ZZW2(:),0.0)
+ ZZW2(:) = MAX(ZZW2(:),0.0)
+ END IF
!
! Adjust sources
!
@@ -271,12 +283,20 @@ IF( IEVAP >= 1 ) THEN
ZW(:,:,:)= PEVAP3D(:,:,:)
PEVAP3D(:,:,:) = UNPACK( ZZW2(:),MASK=GEVAP(:,:,:),FIELD=ZW(:,:,:) )
!
+ IF (LKHKO) THEN
+ ZZW2(:) = MIN(ZZW2(:) * ZCRT(:)/ZRRT(:),ZCRS(:))
+ ZCRS(:) = ZCRS(:) - ZZW2(:)
+ ZW(:,:,:) = PCRS(:,:,:)
+ PCRS(:,:,:) = UNPACK( ZCRS(:),MASK=GEVAP(:,:,:),FIELD=ZW(:,:,:) )
+ ENDIF
+
DEALLOCATE(ZRCT)
DEALLOCATE(ZRRT)
DEALLOCATE(ZRVT)
DEALLOCATE(ZCRT)
DEALLOCATE(ZRVS)
DEALLOCATE(ZRRS)
+ DEALLOCATE(ZCRS)
DEALLOCATE(ZTHS)
DEALLOCATE(ZZLV)
DEALLOCATE(ZZT)
@@ -295,21 +315,40 @@ IF( IEVAP >= 1 ) THEN
! ---------------------------------------
!
!
- GEVAP(:,:,:) = PRRS(:,:,:)>ZRTMIN(3) .AND. PCRS(:,:,:)>ZCTMIN(3)
- ZDR(:,:,:) = 9999.
- WHERE (GEVAP(:,:,:))
- ZDR(:,:,:)=(6.*PRRS(:,:,:)/XPI/XRHOLW/PCRS(:,:,:))**0.33
- ZWLBDR3(:,:,:) = XLBR * PCRS(:,:,:) / PRRS(:,:,:)
- ZWLBDR(:,:,:) = ZWLBDR3(:,:,:)**XLBEXR
- END WHERE
+ IF (LKHKO) THEN
+!* correct negative values for rain
+! --------------------------------
+!
+ WHERE (PRRS(:,:,:)<0.)
+ PRCS(:,:,:) = PRCS(:,:,:)+PRRS(:,:,:)
+ PRRS(:,:,:) = 0.
+ PCRS(:,:,:) = 0.
+ END WHERE
+!
+!* REMOVES NON-PHYSICAL LOW VALUES
+ GEVAP(:,:,:) = PRRS(:,:,:)<ZRTMIN(3) .AND. PCRS(:,:,:)< ZCTMIN(3)
+ WHERE (GEVAP(:,:,:))
+ PRVS(:,:,:) = PRVS(:,:,:) + PRRS(:,:,:)
+ PTHS(:,:,:) = PTHS(:,:,:) - PRRS(:,:,:) * ZLV(:,:,:) / &
+ ( PEXNREF(:,:,:)*(XCPD + XCPV*PRVT(:,:,:) + XCL*(PRCT(:,:,:) + PRRT(:,:,:)) ) )
+ PCRS(:,:,:) = 0.0
+ PRRS(:,:,:) = 0.0
+ END WHERE
+ ELSE
+ GEVAP(:,:,:) = PRRS(:,:,:)>ZRTMIN(3) .AND. PCRS(:,:,:)>ZCTMIN(3)
+ ZDR(:,:,:) = 9999.
+ WHERE (GEVAP(:,:,:))
+ ZDR(:,:,:)=(6.*PRRS(:,:,:)/XPI/XRHOLW/PCRS(:,:,:))**0.33
+ ZWLBDR3(:,:,:) = XLBR * PCRS(:,:,:) / PRRS(:,:,:)
+ ZWLBDR(:,:,:) = ZWLBDR3(:,:,:)**XLBEXR
+ END WHERE
!
- WHERE (GEVAP(:,:,:) .AND. ZDR(:,:,:).LT.82.E-6)
- PCCS(:,:,:) = PCCS(:,:,:)+PCRS(:,:,:)
- PCRS(:,:,:) = 0.
- PRCS(:,:,:) = PRCS(:,:,:)+PRRS(:,:,:)
- PRRS(:,:,:) = 0.
- END WHERE
-
+ WHERE (GEVAP(:,:,:) .AND. ZDR(:,:,:).LT.82.E-6)
+ PCCS(:,:,:) = PCCS(:,:,:)+PCRS(:,:,:)
+ PCRS(:,:,:) = 0.
+ PRCS(:,:,:) = PRCS(:,:,:)+PRRS(:,:,:)
+ PRRS(:,:,:) = 0.
+ END WHERE
!!$ GMICRO(:,:,:) = GEVAP(:,:,:) .AND. ZWLBDR(:,:,:)/XACCR1>ZWLBDC3(:,:,:)
!!$ ! the raindrops are too small, that is lower than D_h
!!$ ZFACT = 1.2E4*XACCR1
@@ -341,6 +380,8 @@ IF( IEVAP >= 1 ) THEN
!!$ PRRS(:,:,:) = 0.0
!!$ END WHERE
!
+ END IF ! LKHKO
+ !
END IF ! IEVAP
!
!++cb++
--
GitLab