From 0545f84d6b6d590a99d48e7f18d28595e605518d Mon Sep 17 00:00:00 2001
From: Gaelle DELAUTIER <gaelle.delautier@meteo.fr>
Date: Thu, 15 Dec 2016 16:01:43 +0100
Subject: [PATCH] Christine 15/12/16 : retour arriere compute_updraft
---
src/MNH/compute_bl89_ml.f90 | 128 ++----
src/MNH/compute_entr_detr.f90 | 684 +++++++++++++++--------------
src/MNH/compute_updraft.f90 | 91 ++--
src/MNH/compute_updraft_hrio.f90 | 4 +-
src/MNH/compute_updraft_raha.f90 | 667 ----------------------------
src/MNH/compute_updraft_rhcj10.f90 | 5 +-
src/MNH/default_desfmn.f90 | 4 +-
src/MNH/shallow_mf.f90 | 71 +--
8 files changed, 443 insertions(+), 1211 deletions(-)
delete mode 100644 src/MNH/compute_updraft_raha.f90
diff --git a/src/MNH/compute_bl89_ml.f90 b/src/MNH/compute_bl89_ml.f90
index 4c5124602..4baf16acd 100644
--- a/src/MNH/compute_bl89_ml.f90
+++ b/src/MNH/compute_bl89_ml.f90
@@ -10,7 +10,7 @@ INTERFACE
! ###################################################################
SUBROUTINE COMPUTE_BL89_ML(KKA,KKB,KKE,KKU,KKL,PDZZ2D, &
- PTKEM_DEP,PG_O_THVREF,PVPT,KK,OUPORDN,OFLUX,PLWORK)
+ PTKEM2D,PG_O_THVREF2D,PVPT,KK,OUPORDN,PLWORK)
! ###################################################################
!* 1.1 Declaration of Arguments
@@ -20,25 +20,25 @@ INTEGER, INTENT(IN) :: KKB ! near ground physical inde
INTEGER, INTENT(IN) :: KKE ! uppest atmosphere physical index
INTEGER, INTENT(IN) :: KKU ! uppest atmosphere array index
INTEGER, INTENT(IN) :: KKL ! +1 if grid goes from ground to atmosphere top, -1 otherwise
-REAL, DIMENSION(:,:), INTENT(IN) :: PDZZ2D ! height difference between two mass levels
-REAL, DIMENSION(:), INTENT(IN) :: PTKEM_DEP ! TKE to consume
-REAL, DIMENSION(:), INTENT(IN) :: PG_O_THVREF ! g/ThetaVRef at the departure point
-REAL, DIMENSION(:,:), INTENT(IN) :: PVPT ! ThetaV on mass levels
-INTEGER, INTENT(IN) :: KK ! index of departure level
-LOGICAL, INTENT(IN) :: OUPORDN ! switch to compute upward (true) or
- ! downward (false) mixing length
-LOGICAL, INTENT(IN) :: OFLUX ! Computation must be done from flux level
-REAL, DIMENSION(:), INTENT(OUT) :: PLWORK ! Resulting mixing length
+REAL, DIMENSION(:,:), INTENT(IN) :: PDZZ2D
+REAL, DIMENSION(:,:), INTENT(IN) :: PTKEM2D
+REAL, DIMENSION(:,:), INTENT(IN) :: PG_O_THVREF2D
+REAL, DIMENSION(:,:), INTENT(IN) :: PVPT
+INTEGER, INTENT(IN) :: KK
+LOGICAL, INTENT(IN) :: OUPORDN
+REAL, DIMENSION(:), INTENT(OUT) :: PLWORK
END SUBROUTINE COMPUTE_BL89_ML
END INTERFACE
!
END MODULE MODI_COMPUTE_BL89_ML
-! ######spl
- SUBROUTINE COMPUTE_BL89_ML(KKA,KKB,KKE,KKU,KKL,PDZZ2D, &
- PTKEM_DEP,PG_O_THVREF,PVPT,KK,OUPORDN,OFLUX,PLWORK)
+
+
+! ###################################################################
+ SUBROUTINE COMPUTE_BL89_ML(KKA,KKB,KKE,KKU,KKL,PDZZ2D, &
+ PTKEM2D,PG_O_THVREF2D,PVPT,KK,OUPORDN,PLWORK)
! ###################################################################
!!
!! COMPUTE_BL89_ML routine to:
@@ -52,7 +52,6 @@ END MODULE MODI_COMPUTE_BL89_ML
!! -------------
!! Original 19/01/06
!! S. Riette Jan 2012: support for both order of vertical levels and cleaning
-!! R.Honnert Oct 2016 : Update with AROME
!!
!-------------------------------------------------------------------------------
!
@@ -83,36 +82,35 @@ INTEGER, INTENT(IN) :: KKB ! near ground physical inde
INTEGER, INTENT(IN) :: KKE ! uppest atmosphere physical index
INTEGER, INTENT(IN) :: KKU ! uppest atmosphere array index
INTEGER, INTENT(IN) :: KKL ! +1 if grid goes from ground to atmosphere top, -1 otherwise
-REAL, DIMENSION(:,:), INTENT(IN) :: PDZZ2D ! height difference between two mass levels
-REAL, DIMENSION(:), INTENT(IN) :: PTKEM_DEP ! TKE to consume
-REAL, DIMENSION(:), INTENT(IN) :: PG_O_THVREF ! g/ThetaVRef at the departure point
-REAL, DIMENSION(:,:), INTENT(IN) :: PVPT ! ThetaV on mass levels
-INTEGER, INTENT(IN) :: KK ! index of departure level
-LOGICAL, INTENT(IN) :: OUPORDN ! switch to compute upward (true) or
- ! downward (false) mixing length
-LOGICAL, INTENT(IN) :: OFLUX ! Computation must be done from flux level
-REAL, DIMENSION(:), INTENT(OUT) :: PLWORK ! Resulting mixing length
+REAL, DIMENSION(:,:), INTENT(IN) :: PDZZ2D
+REAL, DIMENSION(:,:), INTENT(IN) :: PTKEM2D
+REAL, DIMENSION(:,:), INTENT(IN) :: PG_O_THVREF2D
+REAL, DIMENSION(:,:), INTENT(IN) :: PVPT
+INTEGER, INTENT(IN) :: KK
+LOGICAL, INTENT(IN) :: OUPORDN
+REAL, DIMENSION(:), INTENT(OUT) :: PLWORK
! 0.2 Local variable
!
-REAL, DIMENSION(SIZE(PVPT,1)) :: ZLWORK1,ZLWORK2 ! Temporary mixing length
-REAL, DIMENSION(SIZE(PVPT,1)) :: ZINTE,ZPOTE ! TKE and potential energy
- ! between 2 levels
-REAL, DIMENSION(SIZE(PVPT,1)) :: ZVPT_DEP ! Thetav on departure point
+REAL, DIMENSION(SIZE(PTKEM2D,1)) :: ZLWORK1,ZLWORK2 ! Temporary mixing length
+REAL, DIMENSION(SIZE(PTKEM2D,1)) :: ZINTE,ZPOTE ! TKE and potential energy
+ ! between 2 levels
!
-REAL, DIMENSION(SIZE(PVPT,1),SIZE(PVPT,2)) :: ZDELTVPT,ZHLVPT
+REAL, DIMENSION(SIZE(PTKEM2D,1),SIZE(PTKEM2D,2)) :: ZDELTVPT,ZHLVPT
!Virtual Potential Temp at Half level and DeltaThv between
- !2 mass levels
+ !2 levels
INTEGER :: IIJU !Internal Domain
INTEGER :: J1D !horizontal loop counter
INTEGER :: JKK !loop counters
+INTEGER :: JRR !moist loop counter
+INTEGER :: JIJK !loop counters
REAL :: ZTEST,ZTEST0,ZTESTM !test for vectorization
!-------------------------------------------------------------------------------------
!
!* 1. INITIALISATION
! --------------
-IIJU=SIZE(PVPT,1)
+IIJU=SIZE(PTKEM2D,1)
!
ZDELTVPT(:,:)=DZM_MF(KKA,KKU,KKL,PVPT(:,:))
ZDELTVPT(:,KKA)=0.
@@ -122,11 +120,6 @@ END WHERE
!
ZHLVPT(:,:)=MZM_MF(KKA,KKU,KKL,PVPT(:,:))
!
-!We consider that gradient between mass levels KKB and KKB+KKL is the same as
-!the gradient between flux level KKB and mass level KKB
-ZDELTVPT(:,KKB)=PDZZ2D(:,KKB)*ZDELTVPT(:,KKB+KKL)/PDZZ2D(:,KKB+KKL)
-ZHLVPT(:,KKB)=PVPT(:,KKB)-ZDELTVPT(:,KKB)*0.5
-!
!
!
!* 2. CALCULATION OF THE UPWARD MIXING LENGTH
@@ -134,58 +127,27 @@ ZHLVPT(:,KKB)=PVPT(:,KKB)-ZDELTVPT(:,KKB)*0.5
!
IF (OUPORDN.EQV..TRUE.) THEN
- ZINTE(:)=PTKEM_DEP(:)
+ ZINTE(:)=PTKEM2D(:,KK)
PLWORK=0.
ZTESTM=1.
- IF(OFLUX)THEN
- ZVPT_DEP(:)=ZHLVPT(:,KK) ! departure point is on flux level
- !We must compute what happens between flux level KK and mass level KK
- DO J1D=1,IIJU
- ZTEST0=0.5+SIGN(0.5,ZINTE(J1D)) ! test if there's energy to consume
- ! Energy consumed if parcel cross the entire layer
- ZPOTE(J1D) = ZTEST0*(PG_O_THVREF(J1D) * &
- (0.5*(ZHLVPT(J1D,KK)+ PVPT(J1D,KK)) - ZVPT_DEP(J1D))) * &
- PDZZ2D(J1D,KK)*0.5
- ! Test if it rests some energy to consume
- ZTEST =0.5+SIGN(0.5,ZINTE(J1D)-ZPOTE(J1D))
- ! Length travelled by parcel if it rests energy to consume
- ZLWORK1(J1D)=PDZZ2D(J1D,KK)*0.5
- ! Lenght travelled by parcel to nullify energy
- ZLWORK2(J1D)= ( - PG_O_THVREF(J1D) * &
- ( ZHLVPT(J1D,KK) - ZVPT_DEP(J1D) ) &
- + SQRT (ABS( &
- ( PG_O_THVREF(J1D) * (ZHLVPT(J1D,KK) - ZVPT_DEP(J1D)) )**2 &
- + 2. * ZINTE(J1D) * PG_O_THVREF(J1D) &
- * ZDELTVPT(J1D,KK) / PDZZ2D(J1D,KK) )) ) / &
- ( PG_O_THVREF(J1D) * ZDELTVPT(J1D,KK) / PDZZ2D(J1D,KK) )
- ! Effective length travelled by parcel
- PLWORK(J1D)=PLWORK(J1D)+ZTEST0*(ZTEST*ZLWORK1(J1D)+ &
- (1-ZTEST)*ZLWORK2(J1D))
- ! Rest of energy to consume
- ZINTE(J1D) = ZINTE(J1D) - ZPOTE(J1D)
- ENDDO
- ELSE
- ZVPT_DEP(:)=PVPT(:,KK) ! departure point is on mass level
- ENDIF
-
DO JKK=KK+KKL,KKE,KKL
IF(ZTESTM > 0.) THEN
ZTESTM=0
DO J1D=1,IIJU
ZTEST0=0.5+SIGN(0.5,ZINTE(J1D))
- ZPOTE(J1D) = ZTEST0*(PG_O_THVREF(J1D) * &
- (ZHLVPT(J1D,JKK) - ZVPT_DEP(J1D))) * PDZZ2D(J1D,JKK) !particle keeps its temperature
+ ZPOTE(J1D) = ZTEST0*(PG_O_THVREF2D(J1D,KK) * &
+ (ZHLVPT(J1D,JKK) - PVPT(J1D,KK))) * PDZZ2D(J1D,JKK) !particle keeps its temperature
ZTEST =0.5+SIGN(0.5,ZINTE(J1D)-ZPOTE(J1D))
ZTESTM=ZTESTM+ZTEST0
ZLWORK1(J1D)=PDZZ2D(J1D,JKK)
!ZLWORK2 jump of the last reached level
- ZLWORK2(J1D)= ( - PG_O_THVREF(J1D) * &
- ( PVPT(J1D,JKK-KKL) - ZVPT_DEP(J1D) ) &
+ ZLWORK2(J1D)= ( - PG_O_THVREF2D(J1D,KK) * &
+ ( PVPT(J1D,JKK-KKL) - PVPT(J1D,KK) ) &
+ SQRT (ABS( &
- ( PG_O_THVREF(J1D) * (PVPT(J1D,JKK-KKL) - ZVPT_DEP(J1D)) )**2 &
- + 2. * ZINTE(J1D) * PG_O_THVREF(J1D) &
+ ( PG_O_THVREF2D(J1D,KK) * (PVPT(J1D,JKK-KKL) - PVPT(J1D,KK)) )**2 &
+ + 2. * ZINTE(J1D) * PG_O_THVREF2D(J1D,KK) &
* ZDELTVPT(J1D,JKK) / PDZZ2D(J1D,JKK) )) ) / &
- ( PG_O_THVREF(J1D) * ZDELTVPT(J1D,JKK) / PDZZ2D(J1D,JKK) )
+ ( PG_O_THVREF2D(J1D,KK) * ZDELTVPT(J1D,JKK) / PDZZ2D(J1D,JKK) )
!
PLWORK(J1D)=PLWORK(J1D)+ZTEST0*(ZTEST*ZLWORK1(J1D)+ &
(1-ZTEST)*ZLWORK2(J1D))
@@ -200,13 +162,7 @@ ENDIF
!
IF (OUPORDN.EQV..FALSE.) THEN
- IF(OFLUX) THEN
- WRITE(*,*) ' STOP'
- WRITE(*,*) ' OFLUX OPTION NOT CODED FOR DOWNWARD MIXING LENGTH'
- CALL ABORT
- STOP
- ENDIF
- ZINTE(:)=PTKEM_DEP(:)
+ ZINTE(:)=PTKEM2D(:,KK)
PLWORK=0.
ZTESTM=1.
DO JKK=KK,KKB,-KKL
@@ -214,18 +170,18 @@ IF (OUPORDN.EQV..FALSE.) THEN
ZTESTM=0
DO J1D=1,IIJU
ZTEST0=0.5+SIGN(0.5,ZINTE(J1D))
- ZPOTE(J1D) = -ZTEST0*(PG_O_THVREF(J1D) * &
+ ZPOTE(J1D) = -ZTEST0*(PG_O_THVREF2D(J1D,KK) * &
(ZHLVPT(J1D,JKK) - PVPT(J1D,KK))) * PDZZ2D(J1D,JKK) !particle keeps its temperature
ZTEST =0.5+SIGN(0.5,ZINTE(J1D)-ZPOTE(J1D))
ZTESTM=ZTESTM+ZTEST0
ZLWORK1(J1D)=PDZZ2D(J1D,JKK)
- ZLWORK2(J1D)= ( + PG_O_THVREF(J1D) * &
+ ZLWORK2(J1D)= ( + PG_O_THVREF2D(J1D,KK) * &
( PVPT(J1D,JKK) - PVPT(J1D,KK) ) &
+ SQRT (ABS( &
- ( PG_O_THVREF(J1D) * (PVPT(J1D,JKK) - PVPT(J1D,KK)) )**2 &
- + 2. * ZINTE(J1D) * PG_O_THVREF(J1D) &
+ ( PG_O_THVREF2D(J1D,KK) * (PVPT(J1D,JKK) - PVPT(J1D,KK)) )**2 &
+ + 2. * ZINTE(J1D) * PG_O_THVREF2D(J1D,KK) &
* ZDELTVPT(J1D,JKK) / PDZZ2D(J1D,JKK) )) ) / &
- ( PG_O_THVREF(J1D) * ZDELTVPT(J1D,JKK) / PDZZ2D(J1D,JKK) )
+ ( PG_O_THVREF2D(J1D,KK) * ZDELTVPT(J1D,JKK) / PDZZ2D(J1D,JKK) )
!
PLWORK(J1D)=PLWORK(J1D)+ZTEST0*(ZTEST*ZLWORK1(J1D)+ &
(1-ZTEST)*ZLWORK2(J1D))
diff --git a/src/MNH/compute_entr_detr.f90 b/src/MNH/compute_entr_detr.f90
index b7346fad9..fda4d89f4 100644
--- a/src/MNH/compute_entr_detr.f90
+++ b/src/MNH/compute_entr_detr.f90
@@ -8,56 +8,45 @@
!
INTERFACE
!
- SUBROUTINE COMPUTE_ENTR_DETR(KK,KKB,KKE,KKL,OTEST,OTESTLCL,&
- HFRAC_ICE,PFRAC_ICE,PRHODREF,&
- PPRE_MINUS_HALF,&
- PPRE_PLUS_HALF,PZZ,PDZZ,&
- PTHVM,PTHLM,PRTM,PW_UP2,PTH_UP,&
- PTHL_UP,PRT_UP,PLUP,&
- PRC_UP,PRI_UP,PTHV_UP,&
- PRSAT_UP,PRC_MIX,PRI_MIX, &
- PENTR,PDETR,PENTR_CLD,PDETR_CLD,&
- PBUO_INTEG_DRY,PBUO_INTEG_CLD,&
- PPART_DRY)
+ SUBROUTINE COMPUTE_ENTR_DETR(KK,KKB,KKE,KKL,OTEST,OTESTLCL,HFRAC_ICE, &
+ PFRAC_ICE,PPABSM,PZZ,PDZZ,&
+ PTHVM,PTHLM,PRTM,PW_UP2,&
+ PTHL_UP,PRT_UP,PLUP,&
+ PRC_UP,PRI_UP,PRC_MIX,PRI_MIX, &
+ PENTR,PDETR,PBUO_INTEG)
-!INTEGER, INTENT(IN) :: KK
+!
+!
+!
+INTEGER, INTENT(IN) :: KK ! near ground physical index
INTEGER, INTENT(IN) :: KKB ! near ground physical index
INTEGER, INTENT(IN) :: KKE ! uppest atmosphere physical index
INTEGER, INTENT(IN) :: KKL ! +1 if grid goes from ground to atmosphere top, -1 otherwise
-LOGICAL,DIMENSION(:), INTENT(IN) :: OTEST ! test to see if updraft is running
-LOGICAL,DIMENSION(:), INTENT(IN) :: OTESTLCL !test of condensation
-CHARACTER*1, INTENT(IN) :: HFRAC_ICE ! frac_ice can be compute using
- ! Temperature (T) or prescribed
- ! (Y)
-REAL, DIMENSION(:), INTENT(IN) :: PFRAC_ICE ! fraction of ice
+LOGICAL,DIMENSION(:),INTENT(INOUT) :: OTEST
+LOGICAL,DIMENSION(:),INTENT(INOUT) :: OTESTLCL !test of condensation
+CHARACTER*1,INTENT(IN) :: HFRAC_ICE
+REAL, DIMENSION(:) ,INTENT(IN) :: PFRAC_ICE
+
!
! prognostic variables at t- deltat
-!
-REAL, DIMENSION(:), INTENT(IN) :: PRHODREF !rhodref
-REAL, DIMENSION(:), INTENT(IN) :: PPRE_MINUS_HALF ! Pressure at flux level KK
-REAL, DIMENSION(:), INTENT(IN) :: PPRE_PLUS_HALF ! Pressure at flux level KK+KKL
-REAL, DIMENSION(:,:), INTENT(IN) :: PZZ ! Height at the flux point
-REAL, DIMENSION(:,:), INTENT(IN) :: PDZZ ! metrics coefficient
-REAL, DIMENSION(:,:), INTENT(IN) :: PTHVM ! ThetaV environment
+REAL, DIMENSION(:,:), INTENT(IN) :: PPABSM ! Pressure
+REAL, DIMENSION(:,:), INTENT(IN) :: PZZ ! Height at the flux point
+REAL, DIMENSION(:,:), INTENT(IN) :: PDZZ ! metrics coefficient
+REAL, DIMENSION(:,:), INTENT(IN) :: PTHVM ! ThetaV environment
!
-! thermodynamical variables which are transformed in conservative var.
!
-REAL, DIMENSION(:,:), INTENT(IN) :: PTHLM ! Thetal
-REAL, DIMENSION(:,:), INTENT(IN) :: PRTM ! total mixing ratio
-REAL, DIMENSION(:,:), INTENT(IN) :: PW_UP2 ! Vertical velocity^2
-REAL, DIMENSION(:), INTENT(IN) :: PTH_UP,PTHL_UP,PRT_UP ! updraft properties
-REAL, DIMENSION(:), INTENT(IN) :: PLUP ! LUP compute from the ground
+! thermodynamical variables which are transformed in conservative var.
+REAL, DIMENSION(:), INTENT(IN) :: PTHLM ! Thetal
+REAL, DIMENSION(:), INTENT(IN) :: PRTM ! total mixing ratio
+REAL, DIMENSION(:,:), INTENT(INOUT) :: PW_UP2 ! Vertical velocity^2
+REAL, DIMENSION(:), INTENT(IN) :: PTHL_UP,PRT_UP ! updraft properties
+REAL, DIMENSION(:), INTENT(IN) :: PLUP ! LUP compute from the ground
REAL, DIMENSION(:), INTENT(IN) :: PRC_UP,PRI_UP ! Updraft cloud content
-REAL, DIMENSION(:), INTENT(IN) :: PTHV_UP ! Thetav of updraft
-REAL, DIMENSION(:), INTENT(IN) :: PRSAT_UP ! Mixing ratio at saturation in updraft
REAL, DIMENSION(:), INTENT(INOUT) :: PRC_MIX, PRI_MIX ! Mixture cloud content
-REAL, DIMENSION(:), INTENT(OUT) :: PENTR ! Mass flux entrainment of the updraft
-REAL, DIMENSION(:), INTENT(OUT) :: PDETR ! Mass flux detrainment of the updraft
-REAL, DIMENSION(:), INTENT(OUT) :: PENTR_CLD ! Mass flux entrainment of the updraft in cloudy part
-REAL, DIMENSION(:), INTENT(OUT) :: PDETR_CLD ! Mass flux detrainment of the updraft in cloudy part
-REAL, DIMENSION(:), INTENT(OUT) :: PBUO_INTEG_DRY, PBUO_INTEG_CLD! Integral Buoyancy
-REAL, DIMENSION(:), INTENT(OUT) :: PPART_DRY ! ratio of dry part at the transition level
+REAL, DIMENSION(:), INTENT(INOUT) :: PENTR ! Mass flux entrainment of the updraft
+REAL, DIMENSION(:), INTENT(INOUT) :: PDETR ! Mass flux detrainment of the updraft
+REAL, DIMENSION(:), INTENT(INOUT) :: PBUO_INTEG ! Integrated Buoyancy
!
!
END SUBROUTINE COMPUTE_ENTR_DETR
@@ -67,16 +56,11 @@ END INTERFACE
END MODULE MODI_COMPUTE_ENTR_DETR
! ######spl
SUBROUTINE COMPUTE_ENTR_DETR(KK,KKB,KKE,KKL,OTEST,OTESTLCL,&
- HFRAC_ICE,PFRAC_ICE,PRHODREF,&
- PPRE_MINUS_HALF,&
- PPRE_PLUS_HALF,PZZ,PDZZ,&
- PTHVM,PTHLM,PRTM,PW_UP2,PTH_UP,&
+ HFRAC_ICE,PFRAC_ICE,PPABSM,PZZ,PDZZ,&
+ PTHVM,PTHLM,PRTM,PW_UP2,&
PTHL_UP,PRT_UP,PLUP,&
- PRC_UP,PRI_UP,PTHV_UP,&
- PRSAT_UP,PRC_MIX,PRI_MIX, &
- PENTR,PDETR,PENTR_CLD,PDETR_CLD,&
- PBUO_INTEG_DRY,PBUO_INTEG_CLD,&
- PPART_DRY)
+ PRC_UP,PRI_UP,PRC_MIX,PRI_MIX, &
+ PENTR,PDETR,PBUO_INTEG)
! #############################################################
!!
@@ -114,12 +98,6 @@ END MODULE MODI_COMPUTE_ENTR_DETR
!! protection against too big ZPART_DRY, interface modified
!! S. Riette Jan 2012: support for both order of vertical levels
!! S. Riette & J. Escobar (11/2013) : remove div by 0 on real*4 case
-!! P.Marguinaud Jun 2012: fix uninitialized variable
-!! P.Marguinaud Nov 2012: fix gfortran bug
-!! S. Riette Apr 2013: bugs correction, rewriting (for optimisation) and
-!! improvement of continuity at the condensation level
-!! S. Riette Nov 2013: protection against zero divide for min value of dry PDETR
-!! R.Honnert Oct 2016 : Update with AROME
!! --------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -143,18 +121,16 @@ INTEGER, INTENT(IN) :: KK
INTEGER, INTENT(IN) :: KKB ! near ground physical index
INTEGER, INTENT(IN) :: KKE ! uppest atmosphere physical index
INTEGER, INTENT(IN) :: KKL ! +1 if grid goes from ground to atmosphere top, -1 otherwise
-LOGICAL,DIMENSION(:), INTENT(IN) :: OTEST ! test to see if updraft is running
-LOGICAL,DIMENSION(:), INTENT(IN) :: OTESTLCL !test of condensation
-CHARACTER*1, INTENT(IN) :: HFRAC_ICE ! frac_ice can be compute using
+LOGICAL,DIMENSION(:),INTENT(INOUT) :: OTEST ! test to see if updraft is running
+LOGICAL,DIMENSION(:),INTENT(INOUT) :: OTESTLCL !test of condensation
+CHARACTER*1,INTENT(IN) :: HFRAC_ICE ! frac_ice can be compute using
! Temperature (T) or prescribed
! (Y)
REAL, DIMENSION(:), INTENT(IN) :: PFRAC_ICE ! fraction of ice
!
! prognostic variables at t- deltat
!
-REAL, DIMENSION(:), INTENT(IN) :: PRHODREF !rhodref
-REAL, DIMENSION(:), INTENT(IN) :: PPRE_MINUS_HALF ! Pressure at flux level KK
-REAL, DIMENSION(:), INTENT(IN) :: PPRE_PLUS_HALF ! Pressure at flux level KK+KKL
+REAL, DIMENSION(:,:), INTENT(IN) :: PPABSM ! Pressure at time t-1
REAL, DIMENSION(:,:), INTENT(IN) :: PZZ ! Height at the flux point
REAL, DIMENSION(:,:), INTENT(IN) :: PDZZ ! metrics coefficient
REAL, DIMENSION(:,:), INTENT(IN) :: PTHVM ! ThetaV environment
@@ -162,21 +138,16 @@ REAL, DIMENSION(:,:), INTENT(IN) :: PTHVM ! ThetaV environment
!
! thermodynamical variables which are transformed in conservative var.
!
-REAL, DIMENSION(:,:), INTENT(IN) :: PTHLM ! Thetal
-REAL, DIMENSION(:,:), INTENT(IN) :: PRTM ! total mixing ratio
-REAL, DIMENSION(:,:), INTENT(IN) :: PW_UP2 ! Vertical velocity^2
-REAL, DIMENSION(:), INTENT(IN) :: PTH_UP,PTHL_UP,PRT_UP ! updraft properties
+REAL, DIMENSION(:), INTENT(IN) :: PTHLM ! Thetal
+REAL, DIMENSION(:), INTENT(IN) :: PRTM ! total mixing ratio
+REAL, DIMENSION(:,:), INTENT(INOUT) :: PW_UP2 ! Vertical velocity^2
+REAL, DIMENSION(:), INTENT(IN) :: PTHL_UP,PRT_UP ! updraft properties
REAL, DIMENSION(:), INTENT(IN) :: PLUP ! LUP compute from the ground
REAL, DIMENSION(:), INTENT(IN) :: PRC_UP,PRI_UP ! Updraft cloud content
-REAL, DIMENSION(:), INTENT(IN) :: PTHV_UP ! Thetav of updraft
-REAL, DIMENSION(:), INTENT(IN) :: PRSAT_UP ! Mixing ratio at saturation in updraft
REAL, DIMENSION(:), INTENT(INOUT) :: PRC_MIX, PRI_MIX ! Mixture cloud content
-REAL, DIMENSION(:), INTENT(OUT) :: PENTR ! Mass flux entrainment of the updraft
-REAL, DIMENSION(:), INTENT(OUT) :: PDETR ! Mass flux detrainment of the updraft
-REAL, DIMENSION(:), INTENT(OUT) :: PENTR_CLD ! Mass flux entrainment of the updraft in cloudy part
-REAL, DIMENSION(:), INTENT(OUT) :: PDETR_CLD ! Mass flux detrainment of the updraft in cloudy part
-REAL, DIMENSION(:), INTENT(OUT) :: PBUO_INTEG_DRY, PBUO_INTEG_CLD! Integral Buoyancy
-REAL, DIMENSION(:), INTENT(OUT) :: PPART_DRY ! ratio of dry part at the transition level
+REAL, DIMENSION(:), INTENT(INOUT) :: PENTR ! Mass flux entrainment of the updraft
+REAL, DIMENSION(:), INTENT(INOUT) :: PDETR ! Mass flux detrainment of the updraft
+REAL, DIMENSION(:), INTENT(INOUT) :: PBUO_INTEG! Integral Buoyancy
!
!
! 1.2 Declaration of local variables
@@ -184,38 +155,64 @@ REAL, DIMENSION(:), INTENT(OUT) :: PPART_DRY ! ratio of dry part at the tr
!
! Variables for cloudy part
-REAL, DIMENSION(SIZE(PTHLM,1)) :: ZKIC, ZKIC_F2 ! fraction of env. mass in the muxtures
-REAL, DIMENSION(SIZE(PTHLM,1)) :: ZEPSI,ZDELTA ! factor entrainment detrainment
-REAL, DIMENSION(SIZE(PTHLM,1)) :: ZEPSI_CLOUD ! factor entrainment detrainment
-REAL :: ZCOEFFMF_CLOUD ! factor for compputing entr. detr.
-REAL, DIMENSION(SIZE(PTHLM,1)) :: ZMIXTHL,ZMIXRT ! Thetal and rt in the mixtures
-REAL, DIMENSION(SIZE(PTHLM,1)) :: ZTHMIX ! Theta and Thetav of mixtures
-REAL, DIMENSION(SIZE(PTHLM,1)) :: ZRVMIX,ZRCMIX,ZRIMIX ! mixing ratios in mixtures
-REAL, DIMENSION(SIZE(PTHLM,1)) :: ZTHVMIX, ZTHVMIX_F2 ! Theta and Thetav of mixtures
-REAL, DIMENSION(SIZE(PTHLM,1)) :: ZTHV_UP_F2 ! thv_up at flux point kk+kkl
-REAL, DIMENSION(SIZE(PTHLM,1)) :: ZRSATW, ZRSATI ! working arrays (mixing ratio at saturation)
-REAL, DIMENSION(SIZE(PTHLM,1)) :: ZTHV ! theta V of environment at the bottom of cloudy part
-REAL :: ZKIC_INIT !Initial value of ZKIC
-REAL, DIMENSION(SIZE(PTHLM,1)) :: ZCOTHVU ! Variation of Thvup between bottom and top of cloudy part
+
+REAL, DIMENSION(SIZE(PTHLM)) :: ZKIC ! fraction of env. mass in the muxtures
+REAL, DIMENSION(SIZE(PTHLM)) :: ZEPSI,ZDELTA ! factor entrainment detrainment
+REAL, DIMENSION(SIZE(PTHLM)) :: ZEPSI_CLOUD ! factor entrainment detrainment
+REAL, DIMENSION(SIZE(PTHLM)) :: ZCOEFFMF_CLOUD ! factor for compputing entr. detr.
+
+REAL, DIMENSION(SIZE(PTHLM)) :: ZMIXTHL,ZMIXRT ! Thetal and rt in the mixtures
+!
+REAL, DIMENSION(SIZE(PTHLM)) :: ZTHMIX ! Theta and Thetav of mixtures
+REAL, DIMENSION(SIZE(PTHLM)) :: ZRVMIX,ZRCMIX,ZRIMIX ! mixing ratios in mixtures
+
+REAL, DIMENSION(SIZE(PTHLM)) :: ZTHMIX_F2 ! Theta and Thetav of mixtures
+REAL, DIMENSION(SIZE(PTHLM)) :: ZRVMIX_F2,ZRCMIX_F2,ZRIMIX_F2 ! mixing ratios in mixtures
+
+REAL, DIMENSION(SIZE(PTHLM)) :: ZTHMIX_M2
+REAL, DIMENSION(SIZE(PTHLM)) :: ZRVMIX_M2, ZRCMIX_M2, ZRIMIX_M2
+
+REAL, DIMENSION(SIZE(PTHLM)) :: ZTHV_UP ! thvup at mass point kk
+
+
+REAL, DIMENSION(SIZE(PTHLM)) :: ZTHVMIX_1,ZTHVMIX_2 ! Theta and Thetav of mixtures
+
! Variables for dry part
-REAL, DIMENSION(SIZE(PTHLM,1)) :: ZFOESW, ZFOESI ! saturating vapor pressure
-REAL, DIMENSION(SIZE(PTHLM,1)) :: ZDRSATODP ! d.Rsat/dP
-REAL, DIMENSION(SIZE(PTHLM,1)) :: ZT ! Temperature
-REAL, DIMENSION(SIZE(PTHLM,1)) :: ZWK ! Work array
-
-! Variables for dry and cloudy parts
-REAL, DIMENSION(SIZE(PTHLM,1)) :: ZCOEFF_MINUS_HALF,& ! Variation of Thv between mass points kk-kkl and kk
- ZCOEFF_PLUS_HALF ! Variation of Thv between mass points kk and kk+kkl
-REAL, DIMENSION(SIZE(PTHLM,1)) :: ZPRE ! pressure at the bottom of the cloudy part
-REAL, DIMENSION(SIZE(PTHVM,1)) :: ZG_O_THVREF
-REAL, DIMENSION(SIZE(PTHLM,1)) :: ZFRAC_ICE ! fraction of ice
-REAL :: ZRVORD ! RV/RD
-REAL, DIMENSION(SIZE(PTHLM,1)) :: ZDZ_STOP,& ! Exact Height of the LCL above flux level KK
+
+REAL, DIMENSION(SIZE(PTHLM)) :: ZBUO_INTEG,& ! Temporary integral Buoyancy
+ ZDZ_HALF,& ! half-DeltaZ between 2 flux points
+ ZDZ_STOP,& ! Exact Height of the LCL
ZTHV_MINUS_HALF,& ! Thv at flux point(kk)
ZTHV_PLUS_HALF,& ! Thv at flux point(kk+kkl)
- ZDZ ! Delta Z used in computations
-INTEGER :: JI,JLOOP
+ ZCOEFF_MINUS_HALF,& ! Variation of Thv between mass points kk-kkl and kk
+ ZCOEFF_PLUS_HALF,& ! Variation of Thv between mass points kk and kk+kkl
+ ZCOTHVU_MINUS_HALF,& ! Variation of Thvup between flux point kk and mass point kk
+ ZCOTHVU_PLUS_HALF,& ! Variation of Thvup between mass point kk and flux point kk+kkl
+ ZW2_HALF,& ! w**2 at mass point KK
+ ZWK ! temp correction for Lup - z
+
+REAL, DIMENSION(SIZE(PTHLM)) :: ZCOPRE_MINUS_HALF,& ! Variation of pressure between mass points kk-kkl and kk
+ ZCOPRE_PLUS_HALF,& ! Variation of pressure between mass points kk and kk+kkl
+ ZPRE_MINUS_HALF,& ! pressure at flux point kk
+ ZPRE_PLUS_HALF,& ! pressure at flux point kk+kkl
+ ZTHV_UP_F1,& ! thv_up at flux point kk
+ ZTHV_UP_F2 ! thv_up at flux point kk+kkl
+REAL, DIMENSION(SIZE(PTHLM)) :: ZCOEFF_QSAT,& ! variation of Qsat at the transition between dry part and cloudy part
+ ZRC_ORD,& !
+ ZPART_DRY ! part of dry part at the transition level
+!
+REAL, DIMENSION(SIZE(PTHVM,1),SIZE(PTHVM,2)) ::ZG_O_THVREF
+!
+REAL, DIMENSION(SIZE(PTHLM)) :: ZFRAC_ICE ! fraction of ice
+REAL, DIMENSION(SIZE(PTHLM)) :: ZRSATW, ZRSATI
+!
+LOGICAL, DIMENSION(SIZE(OTEST,1)) :: GTEST_LOCAL_LCL,& ! true if LCL found between flux point KK and mass point KK
+ GTEST_LOCAL_LCL2 ! true if LCL found between mass point KK and flux point KK+KKL
+!
+REAL :: ZRDORV ! RD/RV
+REAL :: ZRVORD ! RV/RD
+
!----------------------------------------------------------------------------------
@@ -223,260 +220,269 @@ INTEGER :: JI,JLOOP
! ------------------
+ ZRDORV = XRD / XRV !=0.622
ZRVORD = XRV / XRD !=1.607
- ZG_O_THVREF(:)=XG/PTHVM(:,KK)
- ZCOEFFMF_CLOUD=XENTR_MF * XG / XCRAD_MF
+ ZG_O_THVREF=XG/PTHVM
+ ZCOEFF_QSAT=0.
+ ZRC_ORD=0.
+ ZPART_DRY=1.
+ GTEST_LOCAL_LCL=.FALSE.
+ ZDZ_HALF(:) = (PZZ(:,KK+KKL)-PZZ(:,KK))/2.
+ ZDZ_STOP(:) = ZDZ_HALF(:)
+
ZFRAC_ICE(:)=PFRAC_ICE(:) ! to not modify fraction of ice
- ZPRE(:)=PPRE_MINUS_HALF(:)
- ZMIXTHL(:)=0.1
- ZMIXRT(:)=0.1
-
-! 1.4 Estimation of PPART_DRY
- DO JLOOP=1,SIZE(OTEST)
- IF(OTEST(JLOOP) .AND. OTESTLCL(JLOOP)) THEN
- !No dry part when condensation level is reached
- PPART_DRY(JLOOP)=0.
- ZDZ_STOP(JLOOP)=0.
- ZPRE(JLOOP)=PPRE_MINUS_HALF(JLOOP)
- ELSE IF (OTEST(JLOOP) .AND. .NOT. OTESTLCL(JLOOP)) THEN
- !Temperature at flux level KK
- ZT(JLOOP)=PTH_UP(JLOOP)*(PPRE_MINUS_HALF(JLOOP)/XP00) ** (XRD/XCPD)
- !Saturating vapor pressure at flux level KK
- ZFOESW(JLOOP) = MIN(EXP( XALPW - XBETAW/ZT(JLOOP) - XGAMW*LOG(ZT(JLOOP)) ), 0.99*PPRE_MINUS_HALF(JLOOP))
- ZFOESI(JLOOP) = MIN(EXP( XALPI - XBETAI/ZT(JLOOP) - XGAMI*LOG(ZT(JLOOP)) ), 0.99*PPRE_MINUS_HALF(JLOOP))
- !Computation of d.Rsat / dP (partial derivations with respect to P and T
- !and use of T=Theta*(P/P0)**(R/Cp) to transform dT into dP with theta_up
- !constant at the vertical)
- ZDRSATODP(JLOOP)=(XBETAW/ZT(JLOOP)-XGAMW)*(1-ZFRAC_ICE(JLOOP))+(XBETAI/ZT(JLOOP)-XGAMI)*ZFRAC_ICE(JLOOP)
- ZDRSATODP(JLOOP)=((XRD/XCPD)*ZDRSATODP(JLOOP)-1.)*PRSAT_UP(JLOOP)/ &
- &(PPRE_MINUS_HALF(JLOOP)-(ZFOESW(JLOOP)*(1-ZFRAC_ICE(JLOOP)) + ZFOESI(JLOOP)*ZFRAC_ICE(JLOOP)))
- !Use of d.Rsat / dP and pressure at flux level KK to find pressure (ZPRE)
- !where Rsat is equal to PRT_UP
- ZPRE(JLOOP)=PPRE_MINUS_HALF(JLOOP)+(PRT_UP(JLOOP)-PRSAT_UP(JLOOP))/ZDRSATODP(JLOOP)
- !Fraction of dry part (computed with pressure and used with heights, no
- !impact found when using log function here and for pressure on flux levels
- !computation)
- PPART_DRY(JLOOP)=MAX(0., MIN(1., (PPRE_MINUS_HALF(JLOOP)-ZPRE(JLOOP))/(PPRE_MINUS_HALF(JLOOP)-PPRE_PLUS_HALF(JLOOP))))
- !Height above flux level KK of the cloudy part
- ZDZ_STOP(JLOOP) = (PZZ(JLOOP,KK+KKL)-PZZ(JLOOP,KK))*PPART_DRY(JLOOP)
- END IF
- END DO
-
-! 1.5 Gradient and flux values of thetav
+ ZKIC(:)=0.1 ! starting value for critical mixed fraction for CLoudy Part
+
+
+! Computation of KIC
+! ---------------------
+
+! 2.1 Compute critical mixed fraction by estimating unknown
+! T^mix r_c^mix and r_i^mix from thl^mix and r_t^mix
+! We determine the zero crossing of the linear curve
+! evaluating the derivative using ZMIXF=0.1.
+! -----------------------------------------------------
+
+ ZMIXTHL(:) = ZKIC(:) * PTHLM(:)+(1. - ZKIC(:))*PTHL_UP(:)
+ ZMIXRT(:) = ZKIC(:) * PRTM(:)+(1. - ZKIC(:))*PRT_UP(:)
+
+ ! MIXTURE FOR CLOUDY PART
+ ! Compute pressure at flux level KK
+ ZCOPRE_PLUS_HALF(:) = ((PPABSM(:,KK+KKL)-PPABSM(:,KK))/PDZZ(:,KK+KKL))
+ ZPRE_PLUS_HALF(:) = ZCOPRE_PLUS_HALF*0.5*(PZZ(:,KK+KKL)-PZZ(:,KK))+PPABSM(:,KK)
+
+ ! Compute pressure at flux level KK+KKL
IF(KK/=KKB)THEN
- ZCOEFF_MINUS_HALF(:)=((PTHVM(:,KK)-PTHVM(:,KK-KKL))/PDZZ(:,KK))
- ZTHV_MINUS_HALF(:) = PTHVM(:,KK) - ZCOEFF_MINUS_HALF(:)*0.5*(PZZ(:,KK+KKL)-PZZ(:,KK))
+ ZCOPRE_MINUS_HALF(:) = ((PPABSM(:,KK)-PPABSM(:,KK-KKL))/PDZZ(:,KK))
+ ZPRE_MINUS_HALF(:)= ZCOPRE_MINUS_HALF*0.5*(PZZ(:,KK)-PZZ(:,KK-KKL))+PPABSM(:,KK-KKL)
ELSE
- ZCOEFF_MINUS_HALF(:)=0.
- ZTHV_MINUS_HALF(:) = PTHVM(:,KK)
+ ZPRE_MINUS_HALF(:)=PPABSM(:,KK)
ENDIF
- ZCOEFF_PLUS_HALF(:) = ((PTHVM(:,KK+KKL)-PTHVM(:,KK))/PDZZ(:,KK+KKL))
- ZTHV_PLUS_HALF(:) = PTHVM(:,KK) + ZCOEFF_PLUS_HALF(:)*0.5*(PZZ(:,KK+KKL)-PZZ(:,KK))
-
-! 2 Dry part computation:
-! Integral buoyancy and computation of PENTR and PDETR for dry part
-! --------------------------------------------------------------------
-
-DO JLOOP=1,SIZE(OTEST)
- IF (OTEST(JLOOP) .AND. PPART_DRY(JLOOP)>0.) THEN
- ZDZ(JLOOP)=MIN(ZDZ_STOP(JLOOP),(PZZ(JLOOP,KK+KKL)-PZZ(JLOOP,KK))*0.5)
- PBUO_INTEG_DRY(JLOOP) = ZG_O_THVREF(JLOOP)*ZDZ(JLOOP)*&
- (0.5 * ( - ZCOEFF_MINUS_HALF(JLOOP))*ZDZ(JLOOP) &
- - ZTHV_MINUS_HALF(JLOOP) + PTHV_UP(JLOOP) )
-
- !Between mass flux KK and bottom of cloudy part (if above mass flux)
- ZDZ(JLOOP)=MAX(0., ZDZ_STOP(JLOOP)-(PZZ(JLOOP,KK+KKL)-PZZ(JLOOP,KK))*0.5)
- PBUO_INTEG_DRY(JLOOP) = PBUO_INTEG_DRY(JLOOP) + ZG_O_THVREF(JLOOP)*ZDZ(JLOOP)*&
- (0.5 * ( - ZCOEFF_PLUS_HALF(JLOOP))*ZDZ(JLOOP) &
- - PTHVM(JLOOP,KK) + PTHV_UP(JLOOP) )
- IF (PBUO_INTEG_DRY(JLOOP)>=0.) THEN
- PENTR(JLOOP) = 0.5/(XABUO-XBENTR*XENTR_DRY)*&
- LOG(1.+ (2.*(XABUO-XBENTR*XENTR_DRY)/PW_UP2(JLOOP,KK))* &
- PBUO_INTEG_DRY(JLOOP))
- PDETR(JLOOP) = 0.
- ELSE
- PENTR(JLOOP) = 0.
- PDETR(JLOOP) = 0.5/(XABUO)*&
- LOG(1.+ (2.*(XABUO)/PW_UP2(JLOOP,KK))* &
- (-PBUO_INTEG_DRY(JLOOP)))
- ENDIF
- PENTR(JLOOP) = XENTR_DRY*PENTR(JLOOP)/(PZZ(JLOOP,KK+KKL)-PZZ(JLOOP,KK))
- PDETR(JLOOP) = XDETR_DRY*PDETR(JLOOP)/(PZZ(JLOOP,KK+KKL)-PZZ(JLOOP,KK))
- !Minimum value of detrainment
- ZWK(JLOOP)=PLUP(JLOOP)-0.5*(PZZ(JLOOP,KK)+PZZ(JLOOP,KK+KKL))
- ZWK(JLOOP)=SIGN(MAX(1., ABS(ZWK(JLOOP))), ZWK(JLOOP)) ! ZWK must not be zero
- PDETR(JLOOP) = MAX(PPART_DRY(JLOOP)*XDETR_LUP/ZWK(JLOOP), PDETR(JLOOP))
- ELSE
- !No dry part, consation reached (OTESTLCL)
- PBUO_INTEG_DRY(JLOOP) = 0.
- PENTR(JLOOP)=0.
- PDETR(JLOOP)=0.
- END IF
-ENDDO
-
-
-! 3 Wet part computation
-! -----------------------
-
-! 3.1 Integral buoyancy for cloudy part
- ! Compute theta_v of updraft at flux level KK+KKL
- !MIX variables are used to avoid declaring new variables
- !but we are dealing with updraft and not mixture
- ZRCMIX(:)=PRC_UP(:)
- ZRIMIX(:)=PRI_UP(:)
+ ! Compute non cons. var. of mixture at the mass level
CALL TH_R_FROM_THL_RT_1D(HFRAC_ICE,ZFRAC_ICE,&
- PPRE_PLUS_HALF,PTHL_UP,PRT_UP,&
+ PPABSM(:,KK),ZMIXTHL,ZMIXRT,&
+ ZTHMIX,ZRVMIX,PRC_MIX,PRI_MIX,&
+ ZRSATW, ZRSATI)
+
+ ! Compute theta_v of mixture at mass level KK for KF90
+ ZTHVMIX_1(:) = ZTHMIX(:)*(1.+ZRVORD*ZRVMIX(:))/(1.+ZMIXRT(:))
+
+ ! Compute non cons. var. of mixture at the flux level KK+KKL
+ ZRCMIX=PRC_MIX
+ ZRIMIX=PRI_MIX
+ CALL TH_R_FROM_THL_RT_1D(HFRAC_ICE,ZFRAC_ICE,&
+ ZPRE_PLUS_HALF,ZMIXTHL,ZMIXRT,&
+ ZTHMIX,ZRVMIX,ZRCMIX,ZRIMIX,&
+ ZRSATW, ZRSATI)
+
+
+ ! compute theta_v of mixture at the flux level KK+KKL for KF90
+ ZTHVMIX_2(:) = ZTHMIX(:)*(1.+ZRVORD*ZRVMIX(:))/(1.+ZMIXRT(:))
+
+
+! 2.1 Compute critical mixed fraction by estimating unknown
+! T^mix r_c^mix and r_i^mix from thl^mix and r_t^mix
+! We determine the zero crossing of the linear curve
+! evaluating the derivative using ZMIXF=0.1.
+! -----------------------------------------------------
+
+
+! THV_UP FOR DRY PART
+ ! Compute theta_v of updraft at flux level KK
+ ZRCMIX=PRC_UP
+ ZRIMIX=PRI_UP
+ CALL TH_R_FROM_THL_RT_1D(HFRAC_ICE,ZFRAC_ICE,&
+ ZPRE_MINUS_HALF,PTHL_UP,PRT_UP,&
ZTHMIX,ZRVMIX,ZRCMIX,ZRIMIX,&
ZRSATW, ZRSATI)
- ZTHV_UP_F2(:) = ZTHMIX(:)*(1.+ZRVORD*ZRVMIX(:))/(1.+PRT_UP(:))
-
- ! Integral buoyancy for cloudy part
- DO JLOOP=1,SIZE(OTEST)
- IF(OTEST(JLOOP) .AND. PPART_DRY(JLOOP)<1.) THEN
- !Gradient of Theta V updraft over the cloudy part, assuming that thetaV updraft don't change
- !between flux level KK and bottom of cloudy part
- ZCOTHVU(JLOOP)=(ZTHV_UP_F2(JLOOP)-PTHV_UP(JLOOP))/((PZZ(JLOOP,KK+KKL)-PZZ(JLOOP,KK))*(1-PPART_DRY(JLOOP)))
-
- !Computation in two parts to use change of gradient of theta v of environment
- !Between bottom of cloudy part (if under mass level) and mass level KK
- ZDZ(JLOOP)=MAX(0., 0.5*(PZZ(JLOOP,KK+KKL)-PZZ(JLOOP,KK))-ZDZ_STOP(JLOOP))
- PBUO_INTEG_CLD(JLOOP) = ZG_O_THVREF(JLOOP)*ZDZ(JLOOP)*&
- (0.5*( ZCOTHVU(JLOOP) - ZCOEFF_MINUS_HALF(JLOOP))*ZDZ(JLOOP) &
- - (PTHVM(JLOOP,KK)-ZDZ(JLOOP)*ZCOEFF_MINUS_HALF(JLOOP)) + PTHV_UP(JLOOP) )
-
- !Between max(mass level, bottom of cloudy part) and flux level KK+KKL
- ZDZ(JLOOP)=(PZZ(JLOOP,KK+KKL)-PZZ(JLOOP,KK))-MAX(ZDZ_STOP(JLOOP),0.5*(PZZ(JLOOP,KK+KKL)-PZZ(JLOOP,KK)))
- PBUO_INTEG_CLD(JLOOP) = PBUO_INTEG_CLD(JLOOP)+ZG_O_THVREF(JLOOP)*ZDZ(JLOOP)*&
- (0.5*( ZCOTHVU(JLOOP) - ZCOEFF_PLUS_HALF(JLOOP))*ZDZ(JLOOP)&
- - (PTHVM(JLOOP,KK)+(0.5*((PZZ(JLOOP,KK+KKL)-PZZ(JLOOP,KK)))-ZDZ(JLOOP))*ZCOEFF_PLUS_HALF(JLOOP)) +&
- PTHV_UP(JLOOP) )
-
- ELSE
- !No cloudy part
- PBUO_INTEG_CLD(JLOOP)=0.
- END IF
- END DO
-
-! 3.2 Critical mixed fraction for KK+KKL flux level (ZKIC_F2) and
-! for bottom of cloudy part (ZKIC), then a mean for the cloudy part
-! (put also in ZKIC)
-!
-! computation by estimating unknown
-! T^mix r_c^mix and r_i^mix from enthalpy^mix and r_w^mix
-! We determine the zero crossing of the linear curve
-! evaluating the derivative using ZMIXF=0.1
-
- ZKIC_INIT=0.1 ! starting value for critical mixed fraction for CLoudy Part
-
- ! Compute thetaV of environment at the bottom of cloudy part
- ! and cons then non cons. var. of mixture at the bottom of cloudy part
-
- ! JI computed to avoid KKL(KK-KKL) being < KKL*KKB
- JI=KKL*MAX(KKL*(KK-KKL),KKL*KKB)
- DO JLOOP=1,SIZE(OTEST)
- IF(OTEST(JLOOP) .AND. PPART_DRY(JLOOP)>0.5) THEN
- ZDZ(JLOOP)=ZDZ_STOP(JLOOP)-0.5*(PZZ(JLOOP,KK+KKL)-PZZ(JLOOP,KK))
- ZTHV(JLOOP)= PTHVM(JLOOP,KK)+ZCOEFF_PLUS_HALF(JLOOP)*ZDZ(JLOOP)
- ZMIXTHL(JLOOP) = ZKIC_INIT * &
- (PTHLM(JLOOP,KK)+ZDZ(JLOOP)*(PTHLM(JLOOP,KK+KKL)-PTHLM(JLOOP,KK))/PDZZ(JLOOP,KK+KKL)) + &
- (1. - ZKIC_INIT)*PTHL_UP(JLOOP)
- ZMIXRT(JLOOP) = ZKIC_INIT * &
- (PRTM(JLOOP,KK)+ZDZ(JLOOP)*(PRTM(JLOOP,KK+KKL)-PRTM(JLOOP,KK))/PDZZ(JLOOP,KK+KKL)) + &
- (1. - ZKIC_INIT)*PRT_UP(JLOOP)
- ELSEIF(OTEST(JLOOP)) THEN
- ZDZ(JLOOP)=0.5*(PZZ(JLOOP,KK+KKL)-PZZ(JLOOP,KK))-ZDZ_STOP(JLOOP)
- ZTHV(JLOOP)= PTHVM(JLOOP,KK)-ZCOEFF_MINUS_HALF(JLOOP)*ZDZ(JLOOP)
- ZMIXTHL(JLOOP) = ZKIC_INIT * &
- (PTHLM(JLOOP,KK)-ZDZ(JLOOP)*(PTHLM(JLOOP,KK)-PTHLM(JLOOP,JI))/PDZZ(JLOOP,KK)) + &
- (1. - ZKIC_INIT)*PTHL_UP(JLOOP)
- ZMIXRT(JLOOP) = ZKIC_INIT * &
- (PRTM(JLOOP,KK)-ZDZ(JLOOP)*(PRTM(JLOOP,KK)-PRTM(JLOOP,JI))/PDZZ(JLOOP,KK)) + &
- (1. - ZKIC_INIT)*PRT_UP(JLOOP)
- ENDIF
- ENDDO
+ ZTHV_UP_F1(:) = ZTHMIX(:)*(1.+ZRVORD*ZRVMIX(:))/(1.+PRT_UP(:))
+
+ ! Compute theta_v of updraft at mass level KK
+ ZRCMIX=PRC_UP
+ ZRIMIX=PRI_UP
CALL TH_R_FROM_THL_RT_1D(HFRAC_ICE,ZFRAC_ICE,&
- ZPRE,ZMIXTHL,ZMIXRT,&
- ZTHMIX,ZRVMIX,PRC_MIX,PRI_MIX,&
+ PPABSM(:,KK),PTHL_UP,PRT_UP,&
+ ZTHMIX,ZRVMIX,ZRCMIX,ZRIMIX,&
ZRSATW, ZRSATI)
- ZTHVMIX(:) = ZTHMIX(:)*(1.+ZRVORD*ZRVMIX(:))/(1.+ZMIXRT(:))
+ ZTHV_UP(:) = ZTHMIX(:)*(1.+ZRVORD*ZRVMIX(:))/(1.+PRT_UP(:))
- ! Compute cons then non cons. var. of mixture at the flux level KK+KKL with initial ZKIC
- ZMIXTHL(:) = ZKIC_INIT * 0.5*(PTHLM(:,KK)+PTHLM(:,KK+KKL))+(1. - ZKIC_INIT)*PTHL_UP(:)
- ZMIXRT(:) = ZKIC_INIT * 0.5*(PRTM(:,KK)+PRTM(:,KK+KKL))+(1. - ZKIC_INIT)*PRT_UP(:)
+ ! Compute theta_v of updraft at flux level KK+KKL
+ ZRCMIX_F2=PRC_UP
+ ZRIMIX_F2=PRI_UP
CALL TH_R_FROM_THL_RT_1D(HFRAC_ICE,ZFRAC_ICE,&
- PPRE_PLUS_HALF,ZMIXTHL,ZMIXRT,&
- ZTHMIX,ZRVMIX,PRC_MIX,PRI_MIX,&
+ ZPRE_PLUS_HALF,PTHL_UP,PRT_UP,&
+ ZTHMIX_F2,ZRVMIX_F2,ZRCMIX_F2,ZRIMIX_F2,&
ZRSATW, ZRSATI)
- ZTHVMIX_F2(:) = ZTHMIX(:)*(1.+ZRVORD*ZRVMIX(:))/(1.+ZMIXRT(:))
-
- !Computation of mean ZKIC over the cloudy part
- DO JLOOP=1,SIZE(OTEST)
- IF (OTEST(JLOOP)) THEN
- ! Compute ZKIC at the bottom of cloudy part
- ! Thetav_up at bottom is equal to Thetav_up at flux level KK
- IF (ABS(PTHV_UP(JLOOP)-ZTHVMIX(JLOOP))<1.E-10) THEN
- ZKIC(JLOOP)=1.
- ELSE
- ZKIC(JLOOP) = MAX(0.,PTHV_UP(JLOOP)-ZTHV(JLOOP))*ZKIC_INIT / &
- (PTHV_UP(JLOOP)-ZTHVMIX(JLOOP))
- END IF
- ! Compute ZKIC_F2 at flux level KK+KKL
- IF (ABS(ZTHV_UP_F2(JLOOP)-ZTHVMIX_F2(JLOOP))<1.E-10) THEN
- ZKIC_F2(JLOOP)=1.
- ELSE
- ZKIC_F2(JLOOP) = MAX(0.,ZTHV_UP_F2(JLOOP)-ZTHV_PLUS_HALF(JLOOP))*ZKIC_INIT / &
- (ZTHV_UP_F2(JLOOP)-ZTHVMIX_F2(JLOOP))
- END IF
- !Mean ZKIC over the cloudy part
- ZKIC(JLOOP)=MAX(MIN(0.5*(ZKIC(JLOOP)+ZKIC_F2(JLOOP)),1.),0.)
- END IF
- END DO
-
-
-! 3.3 Integration of PDF
-! According to Kain and Fritsch (1990), we replace delta Mt
-! in eq. (7) and (8) using eq. (5). Here we compute the ratio
-! of integrals without computing delta Me
-
- !Constant PDF
- !For this PDF, eq. (5) is delta Me=0.5*delta Mt
- DO JLOOP=1,SIZE(OTEST)
- IF(OTEST(JLOOP)) THEN
- ZEPSI(JLOOP) = ZKIC(JLOOP)**2. !integration multiplied by 2
- ZDELTA(JLOOP) = (1.-ZKIC(JLOOP))**2. !idem
- ENDIF
- ENDDO
-
- !Triangular PDF
- !Calculus must be verified before activating this part, but in this state,
- !results on ARM case are almost identical
- !For this PDF, eq. (5) is also delta Me=0.5*delta Mt
- !WHERE(OTEST)
- ! !Integration multiplied by 2
- ! WHERE(ZKIC<0.5)
- ! ZEPSI(:)=8.*ZKIC(:)**3/3.
- ! ZDELTA(:)=1.-4.*ZKIC(:)**2+8.*ZKIC(:)**3/3.
- ! ELSEWHERE
- ! ZEPSI(:)=5./3.-4*ZKIC(:)**2+8.*ZKIC(:)**3/3.
- ! ZDELTA(:)=8.*(1.-ZKIC(:))**3/3.
- ! ENDWHERE
- !ENDWHERE
-
-! 3.4 Computation of PENTR and PDETR
- DO JLOOP=1,SIZE(OTEST)
- IF(OTEST(JLOOP)) THEN
- ZEPSI_CLOUD=MIN(ZDELTA,ZEPSI)
- PENTR_CLD(JLOOP) = (1.-PPART_DRY(JLOOP))*ZCOEFFMF_CLOUD*PRHODREF(JLOOP)*ZEPSI_CLOUD(JLOOP)
- PDETR_CLD(JLOOP) = (1.-PPART_DRY(JLOOP))*ZCOEFFMF_CLOUD*PRHODREF(JLOOP)*ZDELTA(JLOOP)
- PENTR(JLOOP) = PENTR(JLOOP)+PENTR_CLD(JLOOP)
- PDETR(JLOOP) = PDETR(JLOOP)+PDETR_CLD(JLOOP)
- ELSE
- PENTR_CLD(JLOOP) = 0.
- PDETR_CLD(JLOOP) = 0.
- ENDIF
- ENDDO
+ ZTHV_UP_F2(:) = ZTHMIX_F2(:)*(1.+ZRVORD*ZRVMIX_F2(:))/(1.+PRT_UP(:))
+
+ ! Computation of RC and RI on mass point KK+KKL
+ ZRCMIX_M2=PRC_UP
+ ZRIMIX_M2=PRI_UP
+ CALL TH_R_FROM_THL_RT_1D(HFRAC_ICE,ZFRAC_ICE,&
+ PPABSM(:,KK+KKL),PTHL_UP,PRT_UP,&
+ ZTHMIX_M2,ZRVMIX_M2,ZRCMIX_M2,ZRIMIX_M2,&
+ ZRSATW, ZRSATI)
+
+!
+!* 2.2 Compute final values for entr. and detr.
+! ----------------------------------------
+!
+! Dry PART
+
+ ! Computation of integral entrainment and detrainment between flux level KK
+ ! and mass level KK
+
+ WHERE ((ZRCMIX(:)+ZRIMIX(:)>0.).AND.(.NOT.OTESTLCL))
+! If rc and/or ri is found between flux level KK and mass level KK
+! a part of dry entrainment/detrainment is defined
+! the exact height of LCL is also determined
+ ZCOEFF_QSAT(:) = ((ZRCMIX_F2(:)+ZRIMIX_F2(:)) - (ZRCMIX(:)+ZRIMIX(:)))/ ZDZ_HALF(:)
+ WHERE ((ZCOEFF_QSAT(:)>0.) .OR. (ZCOEFF_QSAT(:)<0.))
+ ZRC_ORD(:) = (ZRCMIX(:)+ZRIMIX(:)) - ZCOEFF_QSAT(:) * ZDZ_HALF(:)
+ ZDZ_STOP = (- ZRC_ORD(:)/ZCOEFF_QSAT(:))
+ ZPART_DRY(:) = MAX(MIN(ZDZ_STOP / (PZZ(:,KK+KKL)-PZZ(:,KK)),0.5),0.)
+ GTEST_LOCAL_LCL(:)=.TRUE.
+ ENDWHERE
+ ENDWHERE
+
+ IF(KK/=KKB)THEN
+ ZCOEFF_MINUS_HALF = ((PTHVM(:,KK)-PTHVM(:,KK-KKL))/PDZZ(:,KK))
+ ELSE
+ ZCOEFF_MINUS_HALF = 0.
+ ENDIF
+ ZCOEFF_PLUS_HALF = ((PTHVM(:,KK+KKL)-PTHVM(:,KK))/PDZZ(:,KK+KKL))
+
+ ZCOTHVU_MINUS_HALF = (ZTHV_UP(:)-ZTHV_UP_F1(:))/ZDZ_HALF(:)
+ ZCOTHVU_PLUS_HALF = (ZTHV_UP_F2(:)-ZTHV_UP(:))/ZDZ_HALF(:)
+
+ IF(KK/=KKB)THEN
+ ZTHV_MINUS_HALF = ZCOEFF_MINUS_HALF*0.5*(PZZ(:,KK)-PZZ(:,KK-KKL))+PTHVM(:,KK-KKL)
+ ZTHV_PLUS_HALF = ZCOEFF_PLUS_HALF*0.5*(PZZ(:,KK)-PZZ(:,KK-KKL))+ ZTHV_MINUS_HALF
+ ELSE
+ ZTHV_MINUS_HALF = PTHVM(:,KK)
+ ZTHV_PLUS_HALF = ZCOEFF_PLUS_HALF*0.5*(PZZ(:,KK+KKL)-PZZ(:,KK))+ ZTHV_MINUS_HALF !according to PZZ computation at KKB-KKL
+ ENDIF
+
+ ! Integral Buoyancy between flux level KK and mass level KK
+ PBUO_INTEG = ZG_O_THVREF(:,KK)*ZDZ_HALF(:)*&
+ (0.5*( ZCOTHVU_MINUS_HALF - ZCOEFF_MINUS_HALF)*ZDZ_HALF(:) &
+ - ZTHV_MINUS_HALF + ZTHV_UP_F1(:) )
+
+ WHERE ((OTEST).AND.(.NOT.OTESTLCL))
+ PENTR=0.
+ PDETR=0.
+
+ ZBUO_INTEG = ZG_O_THVREF(:,KK)*ZDZ_STOP(:)*&
+ (0.5 * ( - ZCOEFF_MINUS_HALF)* ZDZ_STOP(:) &
+ - ZTHV_MINUS_HALF + ZTHV_UP_F1(:) )
+ WHERE (ZBUO_INTEG(:)>=0.)
+ PENTR = 0.5/(XABUO-XBENTR*XENTR_DRY)*&
+ LOG(1.+ (2.*(XABUO-XBENTR*XENTR_DRY)/PW_UP2(:,KK))* &
+ ZBUO_INTEG)
+ PDETR = 0.
+
+ ZW2_HALF = PW_UP2(:,KK) + 2*(XABUO-XBENTR*XENTR_DRY)*(ZBUO_INTEG)
+ ELSEWHERE
+ PENTR = 0.
+ PDETR = 0.5/(XABUO)*&
+ LOG(1.+ (2.*(XABUO)/PW_UP2(:,KK))* &
+ MAX(0.,-ZBUO_INTEG))
+
+ ZW2_HALF = PW_UP2(:,KK) + 2*(XABUO)*(ZBUO_INTEG)
+ ENDWHERE
+ ENDWHERE
+
+
+ ZDZ_STOP(:) = ZDZ_HALF(:)
+
+! total Integral Buoyancy between flux level KK and flux level KK+KKL
+ PBUO_INTEG = PBUO_INTEG + ZG_O_THVREF(:,KK)*ZDZ_HALF(:)*&
+ (0.5*(ZCOTHVU_PLUS_HALF - ZCOEFF_PLUS_HALF)* ZDZ_HALF(:) - &
+ PTHVM(:,KK) + ZTHV_UP(:) )
+
+ IF(KK*KKL<(KKE-KKL)*KKL) THEN !Computation only if we are strictly below KKE-KKL
+ WHERE ((((ZRCMIX_F2(:)+ZRIMIX_F2(:)>0.).AND.(ZRCMIX(:)+ZRIMIX(:)<=0.)).AND.(.NOT.OTESTLCL)).AND.(.NOT.GTEST_LOCAL_LCL(:)))
+ ! If rc and/or ri is found between mass level KK and flux level KK+KKL
+ ! a part of dry entrainment is defined
+ ! the exact height of LCL is also determined
+ ZCOEFF_QSAT(:) = ((ZRCMIX_M2(:)+ZRIMIX_M2(:)) - (ZRCMIX_F2(:)+ZRIMIX_F2(:))) / &
+ & (0.5* (PZZ(:,KK+2*KKL)-PZZ(:,KK+KKL)))
+ !old formulation without ice (and perhaps with errors)
+ !ZCOEFF_QSAT(:) = (PRT_UP(:) - &
+ ! QSAT(ZTHMIX_F2(:)*((PPABSM(:,KK+KKL)/XP00)**(XRD/XCPD)),&
+ ! PPABSM(:,KK+KKL)) - &
+ ! ZRCMIX(:))/ (0.5* (PZZ(:,KK+2*KKL)-PZZ(:,KK+KKL)))
+ WHERE ((ZCOEFF_QSAT(:)>0.) .OR. (ZCOEFF_QSAT(:)<0.))
+ ZRC_ORD(:) = ZRCMIX_F2(:)+ZRIMIX_F2(:) - ZCOEFF_QSAT(:) * ZDZ_HALF(:)
+ ZDZ_STOP = (- ZRC_ORD(:)/ZCOEFF_QSAT(:))
+ ZPART_DRY(:) = 0.5+MAX(MIN(ZDZ_STOP / (PZZ(:,KK+KKL)-PZZ(:,KK)),0.5),0.)
+ GTEST_LOCAL_LCL2(:)=.TRUE.
+ ENDWHERE
+ ENDWHERE
+ ENDIF
+
+ WHERE (((OTEST).AND.(.NOT.OTESTLCL)).AND.(.NOT.GTEST_LOCAL_LCL(:)))
+ ZBUO_INTEG = ZG_O_THVREF(:,KK)*ZDZ_STOP(:)*&
+ (0.5*( - ZCOEFF_PLUS_HALF)* ZDZ_STOP(:)&
+ - PTHVM(:,KK) + ZTHV_UP(:) )
+
+ WHERE (ZW2_HALF>0.)
+ WHERE (ZBUO_INTEG(:)>=0.)
+ PENTR = PENTR + 0.5/(XABUO-XBENTR*XENTR_DRY)* &
+ LOG(1.+ (2.*(XABUO-XBENTR*XENTR_DRY)/ZW2_HALF(:)) * ZBUO_INTEG)
+
+ PDETR = PDETR
+ ELSEWHERE
+ PENTR = PENTR
+ PDETR = PDETR + 0.5/(XABUO)* &
+ LOG(1.+ (2.*(XABUO)/ZW2_HALF(:)) * &
+ MAX(-ZBUO_INTEG,0.))
+ ENDWHERE
+ ELSEWHERE
+ ! if w**2<0 the updraft is stopped
+ OTEST=.FALSE.
+ PENTR = PENTR
+ PDETR = PDETR
+ ENDWHERE
+ ENDWHERE
+ PENTR = XENTR_DRY*PENTR/(PZZ(:,KK+KKL)-PZZ(:,KK))
+ PDETR = XDETR_DRY*PDETR/(PZZ(:,KK+KKL)-PZZ(:,KK))
+
+ ZWK(:)=PLUP(:) - 0.5*(PZZ(:,KK)+PZZ(:,KK+KKL))
+ ZWK(:)=SIGN(MAX(1., ABS(ZWK(:))), ZWK(:))
+ PDETR(:)=MAX(ZPART_DRY(:)*XDETR_LUP/ZWK(:),PDETR(:))
+
+! compute final value of critical mixed fraction using theta_v
+! of mixture, grid-scale and updraft in cloud
+ WHERE ((OTEST).AND.(OTESTLCL))
+ ZKIC(:) = MAX(0.,ZTHV_UP(:)-PTHVM(:,KK))*ZKIC(:) / &
+ (ZTHV_UP(:)-ZTHVMIX_1(:)+XMNH_EPSILON)
+
+ ZKIC(:) = MAX(0., MIN(1., ZKIC(:)))
+
+ ZEPSI(:) = ZKIC(:) **2.
+ ZDELTA(:) = (1.-ZKIC(:))**2.
+ ZEPSI_CLOUD=MIN(ZDELTA,ZEPSI)
+ ZCOEFFMF_CLOUD(:)=XENTR_MF * XG / XCRAD_MF
+ PENTR(:) = ZCOEFFMF_CLOUD(:)*ZEPSI_CLOUD(:)
+ PDETR(:) = ZCOEFFMF_CLOUD(:)*ZDELTA(:)
+ ENDWHERE
+
+! compute final value of critical mixed fraction using theta_v
+! of mixture, grid-scale and updraft in cloud
+ WHERE (((OTEST).AND.(.NOT.(OTESTLCL))).AND.((GTEST_LOCAL_LCL(:).OR.GTEST_LOCAL_LCL2(:))))
+ ZKIC(:) = MAX(0.,ZTHV_UP_F2(:)-ZTHV_PLUS_HALF)*ZKIC(:) / &
+ (ZTHV_UP_F2(:)-ZTHVMIX_2(:)+XMNH_EPSILON)
+ ZKIC(:) = MAX(0., MIN(1., ZKIC(:)))
+ ZEPSI(:) = ZKIC(:) **2.
+ ZDELTA(:) = (1.-ZKIC(:))**2.
+ ZEPSI_CLOUD=MIN(ZDELTA,ZEPSI)
+ ZCOEFFMF_CLOUD(:)=XENTR_MF * XG / XCRAD_MF
+ PENTR(:) = PENTR+(1.-ZPART_DRY(:))*ZCOEFFMF_CLOUD(:)*ZEPSI_CLOUD(:)
+ PDETR(:) = PDETR+(1.-ZPART_DRY(:))*ZCOEFFMF_CLOUD(:)*ZDELTA(:)
+ ENDWHERE
END SUBROUTINE COMPUTE_ENTR_DETR
diff --git a/src/MNH/compute_updraft.f90 b/src/MNH/compute_updraft.f90
index 714f53422..1d0e2a821 100644
--- a/src/MNH/compute_updraft.f90
+++ b/src/MNH/compute_updraft.f90
@@ -10,7 +10,6 @@ INTERFACE
!
! #################################################################
SUBROUTINE COMPUTE_UPDRAFT(KKA,KKB,KKE,KKU,KKL, HFRAC_ICE, &
- HMF_UPDRAFT, &
OENTR_DETR,OMIXUV, &
ONOMIXLG,KSV_LGBEG,KSV_LGEND, &
PZZ,PDZZ, &
@@ -36,7 +35,6 @@ INTEGER, INTENT(IN) :: KKE ! uppest atmosphere physica
INTEGER, INTENT(IN) :: KKU ! uppest atmosphere array index
INTEGER, INTENT(IN) :: KKL ! +1 if grid goes from ground to atmosphere top, -1 otherwise
CHARACTER*1, INTENT(IN) :: HFRAC_ICE ! partition liquid/ice scheme
-CHARACTER (LEN=4), INTENT(IN) :: HMF_UPDRAFT ! Type of Mass Flux Scheme
LOGICAL, INTENT(IN) :: OENTR_DETR! flag to recompute entrainment, detrainment and mass flux
LOGICAL, INTENT(IN) :: OMIXUV ! True if mixing of momentum
LOGICAL, INTENT(IN) :: ONOMIXLG ! False if mixing of lagrangian tracer
@@ -85,7 +83,6 @@ END INTERFACE
END MODULE MODI_COMPUTE_UPDRAFT
! ######spl
SUBROUTINE COMPUTE_UPDRAFT(KKA,KKB,KKE,KKU,KKL,HFRAC_ICE, &
- HMF_UPDRAFT, &
OENTR_DETR,OMIXUV, &
ONOMIXLG,KSV_LGBEG,KSV_LGEND, &
PZZ,PDZZ, &
@@ -132,8 +129,6 @@ END MODULE MODI_COMPUTE_UPDRAFT
!! S. Riette may 2011: ice added, interface modified
!! S. Riette Jan 2012: support for both order of vertical levels
!! V.Masson, C.Lac : 02/2011 : SV_UP initialized by a non-zero value
-!! S. Riette Apr 2013: improvement of continuity at the condensation level
-!! R.Honnert Oct 2016 : Add ZSURF and Update with AROME
!! --------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -147,7 +142,6 @@ USE MODI_TH_R_FROM_THL_RT_1D
USE MODI_SHUMAN_MF
USE MODI_COMPUTE_BL89_ML
-USE MODD_GRID_n, ONLY : XDXHAT, XDYHAT
IMPLICIT NONE
@@ -162,9 +156,8 @@ INTEGER, INTENT(IN) :: KKE ! uppest atmosphere physica
INTEGER, INTENT(IN) :: KKU ! uppest atmosphere array index
INTEGER, INTENT(IN) :: KKL ! +1 if grid goes from ground to atmosphere top, -1 otherwise
CHARACTER*1, INTENT(IN) :: HFRAC_ICE ! partition liquid/ice scheme
-CHARACTER (LEN=4), INTENT(IN) :: HMF_UPDRAFT ! Type of Mass Flux Scheme
-LOGICAL, INTENT(IN) :: OENTR_DETR! flag to recompute entrainment, detrainment and mass flux
-LOGICAL, INTENT(IN) :: OMIXUV ! True if mixing of momentum
+LOGICAL, INTENT(IN) :: OENTR_DETR! flag to recompute entrainment, detrainment and mass flux
+LOGICAL, INTENT(IN) :: OMIXUV ! True if mixing of momentum
LOGICAL, INTENT(IN) :: ONOMIXLG ! False if mixing of lagrangian tracer
INTEGER, INTENT(IN) :: KSV_LGBEG ! first index of lag. tracer
INTEGER, INTENT(IN) :: KSV_LGEND ! last index of lag. tracer
@@ -214,9 +207,7 @@ REAL, DIMENSION(SIZE(PTHM,1),SIZE(PTHM,2)) :: &
ZUM_F,ZVM_F,ZRHO_F, & ! density,momentum
ZPRES_F,ZTHVM_F,ZTHVM, & ! interpolated at the flux point
ZG_O_THVREF, & ! g*ThetaV ref
- ZW_UP2, & ! w**2 of the updraft
- ZBUO_INTEG_DRY, ZBUO_INTEG_CLD,&! Integrated Buoyancy
- ZENTR_CLD,ZDETR_CLD ! wet entrainment and detrainment
+ ZW_UP2 ! w**2 of the updraft
REAL, DIMENSION(SIZE(PSVM,1),SIZE(PTHM,2),SIZE(PSVM,3)) :: &
ZSVM_F ! scalar variables
@@ -234,7 +225,7 @@ REAL :: ZRDORV ! RD/RV
REAL :: ZRVORD ! RV/RD
-REAL, DIMENSION(SIZE(PTHM,1)) :: ZMIX1,ZMIX2,ZMIX3_CLD,ZMIX2_CLD
+REAL, DIMENSION(SIZE(PTHM,1)) :: ZMIX1,ZMIX2,ZMIX3
REAL, DIMENSION(SIZE(PTHM,1)) :: ZLUP ! Upward Mixing length from the ground
@@ -250,16 +241,12 @@ LOGICAL, DIMENSION(SIZE(PTHM,1),SIZE(PTHM,2)) :: GWORK2
INTEGER :: ITEST
-REAL, DIMENSION(SIZE(PTHM,1)) :: ZRC_UP, ZRI_UP, ZRV_UP,&
- ZRSATW, ZRSATI,&
- ZPART_DRY
+REAL, DIMENSION(SIZE(PTHM,1)) :: ZRC_UP, ZRI_UP, ZRV_UP, ZRSATW, ZRSATI
REAL :: ZDEPTH_MAX1, ZDEPTH_MAX2 ! control auto-extinction process
REAL :: ZTMAX,ZRMAX ! control value
-REAL, DIMENSION(SIZE(PTHM,1)) :: ZSURF
-
! Thresholds for the perturbation of
! theta_l and r_t at the first level of the updraft
ZTMAX=2.0
@@ -276,6 +263,7 @@ ZDEPTH_MAX1=3000. ! clouds with depth inferior to this value are keeped untouche
ZDEPTH_MAX2=4000. ! clouds with depth superior to this value are suppressed
! Local variables, internal domain
+
!number of scalar variables
ISV=SIZE(PSVM,3)
@@ -373,13 +361,13 @@ IF (OENTR_DETR) THEN
! Closure assumption for mass flux at KKB level
!
- ZG_O_THVREF(:,:)=XG/ZTHVM_F(:,:)
+ ZG_O_THVREF=XG/ZTHVM_F
! compute L_up
GLMIX=.TRUE.
ZTKEM_F(:,KKB)=0.
- CALL COMPUTE_BL89_ML(KKA,KKB,KKE,KKU,KKL,PDZZ,ZTKEM_F(:,KKB),ZG_O_THVREF(:,KKB),ZTHVM,KKB,GLMIX,.TRUE.,ZLUP)
+ CALL COMPUTE_BL89_ML(KKA,KKB,KKE,KKU,KKL,PDZZ,ZTKEM_F,ZG_O_THVREF,ZTHVM_F,KKB,GLMIX,ZLUP)
ZLUP(:)=MAX(ZLUP(:),1.E-10)
! Compute Buoyancy flux at the ground
@@ -387,14 +375,8 @@ IF (OENTR_DETR) THEN
(0.61*ZTHM_F(:,KKB))*PSFRV(:)
! Mass flux at KKB level (updraft triggered if PSFTH>0.)
- IF (HMF_UPDRAFT=='SURF') THEN
- ZSURF(:)=TANH(1.83*SQRT(XDXHAT(1)*XDYHAT(1))/ZLUP)
- ELSE
- ZSURF(:)=1.
- END IF
WHERE (ZWTHVSURF(:)>0.)
- PEMF(:,KKB) = XCMF * ZSURF(:) * ZRHO_F(:,KKB) * &
- ((ZG_O_THVREF(:,KKB))*ZWTHVSURF*ZLUP)**(1./3.)
+ PEMF(:,KKB) = XCMF * ZRHO_F(:,KKB) * ((ZG_O_THVREF(:,KKB))*ZWTHVSURF*ZLUP)**(1./3.)
PFRAC_UP(:,KKB)=MIN(PEMF(:,KKB)/(SQRT(ZW_UP2(:,KKB))*ZRHO_F(:,KKB)),XFRAC_UP_MAX)
ZW_UP2(:,KKB)=(PEMF(:,KKB)/(PFRAC_UP(:,KKB)*ZRHO_F(:,KKB)))**2
GTEST(:)=.TRUE.
@@ -436,6 +418,7 @@ DO JK=KKB,KKE-KKL,KKL
GTESTLCL(:)=.TRUE.
ENDWHERE
+
! COMPUTE PENTR and PDETR at mass level JK
IF (OENTR_DETR) THEN
IF(JK/=KKB) THEN
@@ -443,20 +426,14 @@ DO JK=KKB,KKE-KKL,KKL
ZRI_MIX(:,JK) = ZRI_MIX(:,JK-KKL) ! guess of Ri of mixture
ENDIF
CALL COMPUTE_ENTR_DETR(JK,KKB,KKE,KKL,GTEST,GTESTLCL,HFRAC_ICE,PFRAC_ICE_UP(:,JK),&
- PRHODREF(:,JK),ZPRES_F(:,JK),ZPRES_F(:,JK+KKL),&
- PZZ(:,:),PDZZ(:,:),ZTHVM(:,:), &
- PTHLM(:,:),PRTM(:,:),ZW_UP2(:,:),ZTH_UP(:,JK), &
+ PPABSM(:,:),PZZ(:,:),PDZZ(:,:),ZTHVM(:,:), &
+ PTHLM(:,JK),PRTM(:,JK),ZW_UP2(:,:), &
PTHL_UP(:,JK),PRT_UP(:,JK),ZLUP(:), &
- PRC_UP(:,JK),PRI_UP(:,JK),PTHV_UP(:,JK),&
- PRSAT_UP(:,JK),ZRC_MIX(:,JK),ZRI_MIX(:,JK), &
- PENTR(:,JK),PDETR(:,JK),ZENTR_CLD(:,JK),ZDETR_CLD(:,JK),&
- ZBUO_INTEG_DRY(:,JK), ZBUO_INTEG_CLD(:,JK), &
- ZPART_DRY(:) )
- PBUO_INTEG(:,JK)=ZBUO_INTEG_DRY(:,JK)+ZBUO_INTEG_CLD(:,JK)
+ PRC_UP(:,JK),PRI_UP(:,JK),ZRC_MIX(:,JK),ZRI_MIX(:,JK), &
+ PENTR(:,JK),PDETR(:,JK),PBUO_INTEG(:,JK) )
IF (JK==KKB) THEN
PDETR(:,JK)=0.
- ZDETR_CLD(:,JK)=0.
ENDIF
! Computation of updraft characteristics at level JK+KKL
@@ -482,8 +459,7 @@ DO JK=KKB,KKE-KKL,KKL
! If the updraft did not stop, compute cons updraft characteritics at jk+KKL
WHERE(GTEST)
ZMIX2(:) = (PZZ(:,JK+KKL)-PZZ(:,JK))*PENTR(:,JK) !&
- ZMIX3_CLD(:) = (PZZ(:,JK+KKL)-PZZ(:,JK))*(1.-ZPART_DRY(:))*ZDETR_CLD(:,JK) !&
- ZMIX2_CLD(:) = (PZZ(:,JK+KKL)-PZZ(:,JK))*(1.-ZPART_DRY(:))*ZENTR_CLD(:,JK)
+ ZMIX3(:) = (PZZ(:,JK+KKL)-PZZ(:,JK))*PDETR(:,JK) !&
PTHL_UP(:,JK+KKL)=(PTHL_UP(:,JK)*(1.-0.5*ZMIX2(:)) + PTHLM(:,JK)*ZMIX2(:)) &
/(1.+0.5*ZMIX2(:))
@@ -546,15 +522,20 @@ DO JK=KKB,KKE-KKL,KKL
! Compute the updraft theta_v, buoyancy and w**2 for level JK+KKL
WHERE(GTEST)
- PTHV_UP(:,JK+KKL) = ZTH_UP(:,JK+KKL)*((1+ZRVORD*PRV_UP(:,JK+KKL))/(1+PRT_UP(:,JK+KKL)))
- WHERE (ZBUO_INTEG_DRY(:,JK)>0.)
- ZW_UP2(:,JK+KKL) = ZW_UP2(:,JK) + 2.*(XABUO-XBENTR*XENTR_DRY)* ZBUO_INTEG_DRY(:,JK)
- ELSEWHERE
- ZW_UP2(:,JK+KKL) = ZW_UP2(:,JK) + 2.*XABUO* ZBUO_INTEG_DRY(:,JK)
- ENDWHERE
- ZW_UP2(:,JK+KKL) = ZW_UP2(:,JK+KKL)*(1.-(XBDETR*ZMIX3_CLD(:)+XBENTR*ZMIX2_CLD(:)))&
- /(1.+(XBDETR*ZMIX3_CLD(:)+XBENTR*ZMIX2_CLD(:))) &
- +2.*(XABUO)*ZBUO_INTEG_CLD(:,JK)/(1.+(XBDETR*ZMIX3_CLD(:)+XBENTR*ZMIX2_CLD(:)))
+ PTHV_UP(:,JK+KKL) = ZTH_UP(:,JK+KKL)*((1+ZRVORD*PRV_UP(:,JK+KKL))/(1+PRT_UP(:,JK+KKL)))
+ WHERE (.NOT.(GTESTLCL))
+ WHERE (PBUO_INTEG(:,JK)>0.)
+ ZW_UP2(:,JK+KKL) = ZW_UP2(:,JK) + 2.*(XABUO-XBENTR*XENTR_DRY)* PBUO_INTEG(:,JK)
+ ENDWHERE
+ WHERE (PBUO_INTEG(:,JK)<=0.)
+ ZW_UP2(:,JK+KKL) = ZW_UP2(:,JK) + 2.*XABUO* PBUO_INTEG(:,JK)
+ ENDWHERE
+ ENDWHERE
+ WHERE (GTESTLCL)
+ ZW_UP2(:,JK+KKL) = ZW_UP2(:,JK)*(1.-(XBDETR*ZMIX3(:)+XBENTR*ZMIX2(:)))&
+ /(1.+(XBDETR*ZMIX3(:)+XBENTR*ZMIX2(:))) &
+ +2.*(XABUO)*PBUO_INTEG(:,JK)/(1.+(XBDETR*ZMIX3(:)+XBENTR*ZMIX2(:)))
+ ENDWHERE
ENDWHERE
@@ -621,14 +602,10 @@ IF(OENTR_DETR) THEN
PDEPTH(JI) = MAX(0., PZZ(JI,KKCTL(JI)) - PZZ(JI,KKLCL(JI)) )
END DO
- GWORK1(:)= (GTESTLCL(:) .AND. (PDEPTH(:) > ZDEPTH_MAX1) )
- GWORK2(:,:) = SPREAD( GWORK1(:), DIM=2, NCOPIES=MAX(KKU,KKA) )
- ZCOEF(:,:) = SPREAD( (1.-(PDEPTH(:)-ZDEPTH_MAX1)/(ZDEPTH_MAX2-ZDEPTH_MAX1)), DIM=2, NCOPIES=SIZE(ZCOEF,2))
- ZCOEF=MIN(MAX(ZCOEF,0.),1.)
- WHERE (GWORK2)
- PEMF(:,:) = PEMF(:,:) * ZCOEF(:,:)
- PFRAC_UP(:,:) = PFRAC_UP(:,:) * ZCOEF(:,:)
- ENDWHERE
-ENDIF
+ENDIF
+
+GWORK1(:)= (GTESTLCL(:) .AND. (PDEPTH(:) > ZDEPTH_MAX1) )
+GWORK2(:,:) = SPREAD( GWORK1(:), DIM=2, NCOPIES=MAX(KKU,KKA) )
+ZCOEF(:,:) = SPREAD( (1.-(PDEPTH(:)-ZDEPTH_MAX1)/(ZDEPTH_MAX2-ZDEPTH_MAX1)), DIM=2, NCOPIES=SIZE(ZCOEF,2))
END SUBROUTINE COMPUTE_UPDRAFT
diff --git a/src/MNH/compute_updraft_hrio.f90 b/src/MNH/compute_updraft_hrio.f90
index c2e0bae75..9fa0ad5f0 100644
--- a/src/MNH/compute_updraft_hrio.f90
+++ b/src/MNH/compute_updraft_hrio.f90
@@ -134,7 +134,6 @@ END MODULE MODI_COMPUTE_UPDRAFT_HRIO! ######spl
!! S. Riette may 2011: ice added, interface modified
!! S. Riette Jan 2012: support for both order of vertical levels
!! V.Masson, C.Lac : 02/2011 : SV_UP initialized by a non-zero value
-!! 10/2016 (R.Honnert and S.Riette) : Improvement of EDKF and adaptation to the grey zone
!! --------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -445,8 +444,7 @@ IF (OENTR_DETR) THEN
GLMIX=.TRUE.
ZTKEM_F(:,KKB)=0.
- CALL COMPUTE_BL89_ML(KKA,KKB,KKE,KKU,KKL,PDZZ,ZTKEM_F(:,KKB),ZG_O_THVREF(:,KKB), &
- ZTHVM_F,KKB,GLMIX,.TRUE.,ZLUP)
+ CALL COMPUTE_BL89_ML(KKA,KKB,KKE,KKU,KKL,PDZZ,ZTKEM_F,ZG_O_THVREF,ZTHVM_F,KKB,GLMIX,ZLUP)
ZLUP(:)=MAX(ZLUP(:),1.E-10)
! Compute Buoyancy flux at the ground
diff --git a/src/MNH/compute_updraft_raha.f90 b/src/MNH/compute_updraft_raha.f90
deleted file mode 100644
index 4ddeb8938..000000000
--- a/src/MNH/compute_updraft_raha.f90
+++ /dev/null
@@ -1,667 +0,0 @@
-!MNH_LIC Copyright 1994-2014 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.
-! ######spl
-! ######spl
- MODULE MODI_COMPUTE_UPDRAFT_RAHA
-! ###########################
-!
-INTERFACE
-!
-! #################################################################
- SUBROUTINE COMPUTE_UPDRAFT_RAHA(KKA,KKB,KKE,KKU,KKL,HFRAC_ICE, &
- OENTR_DETR,OMIXUV, &
- ONOMIXLG,KSV_LGBEG,KSV_LGEND, &
- PZZ,PDZZ, &
- PSFTH,PSFRV, &
- PPABSM,PRHODREF,PUM,PVM, PTKEM, &
- PEXNM,PTHM,PRVM,PTHLM,PRTM, &
- PSVM,PTHL_UP,PRT_UP, &
- PRV_UP,PRC_UP,PRI_UP,PTHV_UP, &
- PW_UP,PU_UP, PV_UP, PSV_UP, &
- PFRAC_UP,PFRAC_ICE_UP,PRSAT_UP, &
- PEMF,PDETR,PENTR, &
- PBUO_INTEG,KKLCL,KKETL,KKCTL, &
- PDEPTH )
-! #################################################################
-!
-!* 1.1 Declaration of Arguments
-!
-!
-INTEGER, INTENT(IN) :: KKA ! near ground array index
-INTEGER, INTENT(IN) :: KKB ! near ground physical index
-INTEGER, INTENT(IN) :: KKE ! uppest atmosphere physical index
-INTEGER, INTENT(IN) :: KKU ! uppest atmosphere array index
-INTEGER, INTENT(IN) :: KKL ! +1 if grid goes from ground to atmosphere top, -1 otherwise
-CHARACTER*1, INTENT(IN) :: HFRAC_ICE ! partition liquid/ice scheme
-LOGICAL, INTENT(IN) :: OENTR_DETR! flag to recompute entrainment, detrainment and mass flux
-LOGICAL, INTENT(IN) :: OMIXUV ! True if mixing of momentum
-LOGICAL, INTENT(IN) :: ONOMIXLG ! False if mixing of lagrangian tracer
-INTEGER, INTENT(IN) :: KSV_LGBEG ! first index of lag. tracer
-INTEGER, INTENT(IN) :: KSV_LGEND ! last index of lag. tracer
-REAL, DIMENSION(:,:), INTENT(IN) :: PZZ ! Height at the flux point
-REAL, DIMENSION(:,:), INTENT(IN) :: PDZZ ! Metrics coefficient
-
-REAL, DIMENSION(:), INTENT(IN) :: PSFTH,PSFRV
-! normal surface fluxes of theta,rv,(u,v) parallel to the orography
-!
-REAL, DIMENSION(:,:), INTENT(IN) :: PPABSM ! Pressure at t-dt
-REAL, DIMENSION(:,:), INTENT(IN) :: PRHODREF ! dry density of the
- ! reference state
-REAL, DIMENSION(:,:), INTENT(IN) :: PUM ! u mean wind
-REAL, DIMENSION(:,:), INTENT(IN) :: PVM ! v mean wind
-REAL, DIMENSION(:,:), INTENT(IN) :: PTKEM ! TKE at t-dt
-
-REAL, DIMENSION(:,:), INTENT(IN) :: PEXNM ! Exner function at t-dt
-REAL, DIMENSION(:,:), INTENT(IN) :: PTHM ! liquid pot. temp. at t-dt
-REAL, DIMENSION(:,:), INTENT(IN) :: PRVM ! vapor mixing ratio at t-dt
-REAL, DIMENSION(:,:), INTENT(IN) :: PTHLM,PRTM ! cons. var. at t-dt
-
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PSVM ! scalar var. at t-dt
-
-REAL, DIMENSION(:,:), INTENT(OUT) :: PTHL_UP,PRT_UP ! updraft properties
-REAL, DIMENSION(:,:), INTENT(OUT) :: PU_UP, PV_UP ! updraft wind components
-REAL, DIMENSION(:,:), INTENT(INOUT):: PRV_UP,PRC_UP ! updraft rv, rc
-REAL, DIMENSION(:,:), INTENT(INOUT):: PRI_UP,PTHV_UP ! updraft ri, THv
-REAL, DIMENSION(:,:), INTENT(INOUT):: PW_UP,PFRAC_UP ! updraft w, fraction
-REAL, DIMENSION(:,:), INTENT(INOUT):: PFRAC_ICE_UP ! liquid/solid fraction in updraft
-REAL, DIMENSION(:,:), INTENT(INOUT):: PRSAT_UP ! Rsat
-
-REAL, DIMENSION(:,:,:), INTENT(OUT) :: PSV_UP ! updraft scalar var.
-
-REAL, DIMENSION(:,:), INTENT(INOUT):: PEMF,PDETR,PENTR ! Mass_flux,
- ! detrainment,entrainment
-REAL, DIMENSION(:,:), INTENT(INOUT) :: PBUO_INTEG ! Integrated Buoyancy
-INTEGER, DIMENSION(:), INTENT(INOUT):: KKLCL,KKETL,KKCTL! LCL, ETL, CTL
-REAL, DIMENSION(:), INTENT(OUT) :: PDEPTH ! Deepness of cloud
-
-
-END SUBROUTINE COMPUTE_UPDRAFT_RAHA
-
-END INTERFACE
-!
-END MODULE MODI_COMPUTE_UPDRAFT_RAHA
-!
-! ######spl
- SUBROUTINE COMPUTE_UPDRAFT_RAHA(KKA,KKB,KKE,KKU,KKL,HFRAC_ICE, &
- OENTR_DETR,OMIXUV, &
- ONOMIXLG,KSV_LGBEG,KSV_LGEND, &
- PZZ,PDZZ, &
- PSFTH,PSFRV, &
- PPABSM,PRHODREF,PUM,PVM, PTKEM, &
- PEXNM,PTHM,PRVM,PTHLM,PRTM, &
- PSVM,PTHL_UP,PRT_UP, &
- PRV_UP,PRC_UP,PRI_UP,PTHV_UP, &
- PW_UP,PU_UP, PV_UP, PSV_UP, &
- PFRAC_UP,PFRAC_ICE_UP,PRSAT_UP, &
- PEMF,PDETR,PENTR, &
- PBUO_INTEG,KKLCL,KKETL,KKCTL, &
- PDEPTH )
-
-! #################################################################
-!!
-!!**** *COMPUTE_UPDRAF_RAHA* - calculates caracteristics of the updraft
-!!
-!!
-!! PURPOSE
-!! -------
-!!**** The purpose of this routine is to build the updraft following Rio et al (2010)
-!! Same as compute_updraft_rhcj10 exept the use of Hourdin et al closure
-!!
-!
-!!** METHOD
-!! ------
-!!
-!! EXTERNAL
-!! --------
-!!
-!! IMPLICIT ARGUMENTS
-!! ------------------
-!!
-!! !! REFERENCE
-!! ---------
-!! Rio et al (2010) (Boundary Layer Meteorol 135:469-483)
-!! Hourdin et al (xxxx)
-!!
-!! AUTHOR
-!! ------
-!! Y. Bouteloup (2012)
-!! R. Honnert Janv 2013 ==> corection of some coding bugs
-!! Y. Bouteloup Janv 2014 ==> Allow the use of loops in the both direction
-!! 10/2016 (R.Honnert and S.Riette) : Improvement of EDKF and adaptation to the grey zone
-!! --------------------------------------------------------------------------
-!
-!* 0. DECLARATIONS
-! ------------
-
-USE MODD_CST
-USE MODD_PARAM_MFSHALL_n
-
-USE MODI_TH_R_FROM_THL_RT_1D
-USE MODI_SHUMAN_MF
-
-IMPLICIT NONE
-
-!* 1.1 Declaration of Arguments
-!
-!
-!
-INTEGER, INTENT(IN) :: KKA ! near ground array index
-INTEGER, INTENT(IN) :: KKB ! near ground physical index
-INTEGER, INTENT(IN) :: KKE ! uppest atmosphere physical index
-INTEGER, INTENT(IN) :: KKU ! uppest atmosphere array index
-INTEGER, INTENT(IN) :: KKL ! +1 if grid goes from ground to atmosphere top, -1 otherwise
-CHARACTER*1, INTENT(IN) :: HFRAC_ICE ! partition liquid/ice scheme
-LOGICAL, INTENT(IN) :: OENTR_DETR! flag to recompute entrainment, detrainment and mass flux
-LOGICAL, INTENT(IN) :: OMIXUV ! True if mixing of momentum
-LOGICAL, INTENT(IN) :: ONOMIXLG ! False if mixing of lagrangian tracer
-INTEGER, INTENT(IN) :: KSV_LGBEG ! first index of lag. tracer
-INTEGER, INTENT(IN) :: KSV_LGEND ! last index of lag. tracer
-REAL, DIMENSION(:,:), INTENT(IN) :: PZZ ! Height at the flux point
-REAL, DIMENSION(:,:), INTENT(IN) :: PDZZ ! Metrics coefficient
-
-REAL, DIMENSION(:), INTENT(IN) :: PSFTH,PSFRV
-! normal surface fluxes of theta,rv,(u,v) parallel to the orography
-!
-REAL, DIMENSION(:,:), INTENT(IN) :: PPABSM ! Pressure at t-dt
-REAL, DIMENSION(:,:), INTENT(IN) :: PRHODREF ! dry density of the
- ! reference state
-REAL, DIMENSION(:,:), INTENT(IN) :: PUM ! u mean wind
-REAL, DIMENSION(:,:), INTENT(IN) :: PVM ! v mean wind
-REAL, DIMENSION(:,:), INTENT(IN) :: PTKEM ! TKE at t-dt
-
-REAL, DIMENSION(:,:), INTENT(IN) :: PEXNM ! Exner function at t-dt
-REAL, DIMENSION(:,:), INTENT(IN) :: PTHM ! liquid pot. temp. at t-dt
-REAL, DIMENSION(:,:), INTENT(IN) :: PRVM ! vapor mixing ratio at t-dt
-REAL, DIMENSION(:,:), INTENT(IN) :: PTHLM,PRTM ! cons. var. at t-dt
-
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PSVM ! scalar var. at t-dt
-
-REAL, DIMENSION(:,:), INTENT(OUT) :: PTHL_UP,PRT_UP ! updraft properties
-REAL, DIMENSION(:,:), INTENT(OUT) :: PU_UP, PV_UP ! updraft wind components
-REAL, DIMENSION(:,:), INTENT(INOUT):: PRV_UP,PRC_UP ! updraft rv, rc
-REAL, DIMENSION(:,:), INTENT(INOUT):: PRI_UP,PTHV_UP ! updraft ri, THv
-REAL, DIMENSION(:,:), INTENT(INOUT):: PW_UP,PFRAC_UP ! updraft w, fraction
-REAL, DIMENSION(:,:), INTENT(INOUT):: PFRAC_ICE_UP ! liquid/solid fraction in updraft
-REAL, DIMENSION(:,:), INTENT(INOUT):: PRSAT_UP ! Rsat
-
-REAL, DIMENSION(:,:,:), INTENT(OUT) :: PSV_UP ! updraft scalar var.
-
-REAL, DIMENSION(:,:), INTENT(INOUT):: PEMF,PDETR,PENTR ! Mass_flux,
- ! detrainment,entrainment
-REAL, DIMENSION(:,:), INTENT(INOUT) :: PBUO_INTEG ! Integrated Buoyancy
-INTEGER, DIMENSION(:), INTENT(INOUT):: KKLCL,KKETL,KKCTL! LCL, ETL, CTL
-REAL, DIMENSION(:), INTENT(OUT) :: PDEPTH ! Deepness of cloud
-! 1.2 Declaration of local variables
-!
-!
-! Mean environment variables at t-dt at flux point
-REAL, DIMENSION(SIZE(PTHM,1),SIZE(PTHM,2)) :: ZTHM_F,ZRVM_F,ZRCM_F ! Theta,rv of
- ! updraft environnement
-REAL, DIMENSION(SIZE(PTHM,1),SIZE(PTHM,2)) :: ZRTM_F, ZTHLM_F, ZTKEM_F ! rt, thetal,TKE,pressure,
-REAL, DIMENSION(SIZE(PTHM,1),SIZE(PTHM,2)) :: ZUM_F,ZVM_F,ZRHO_F ! density,momentum
-REAL, DIMENSION(SIZE(PTHM,1),SIZE(PTHM,2)) :: ZPRES_F,ZTHVM_F,ZTHVM ! interpolated at the flux point
-REAL, DIMENSION(SIZE(PTHM,1),SIZE(PTHM,2)) :: ZG_O_THVREF ! g*ThetaV ref
-REAL, DIMENSION(SIZE(PTHM,1),SIZE(PTHM,2)) :: ZW_UP2 ! w**2 of the updraft
-
-REAL, DIMENSION(SIZE(PSVM,1),SIZE(PTHM,2),SIZE(PSVM,3)) :: ZSVM_F ! scalar variables
-
-
-
-REAL, DIMENSION(SIZE(PTHM,1),SIZE(PTHM,2)) :: ZTH_UP ! updraft THETA
-REAL, DIMENSION(SIZE(PTHM,1)) :: ZT_UP ! updraft T
-REAL, DIMENSION(SIZE(PTHM,1)) :: ZLVOCPEXN ! updraft L
-REAL, DIMENSION(SIZE(PTHM,1)) :: ZCP ! updraft cp
-REAL, DIMENSION(SIZE(PTHM,1),SIZE(PTHM,2)) :: ZBUO ! Buoyancy
-REAL, DIMENSION(SIZE(PTHM,1),SIZE(PTHM,2)) :: ZTHS_UP,ZTHSM
-
-REAL, DIMENSION(SIZE(PTHM,1),SIZE(PTHM,2)) :: ZCOEF ! diminution coefficient for too high clouds
-
-REAL, DIMENSION(SIZE(PSFTH,1) ) :: ZWTHVSURF ! Surface w'thetav'
-
-REAL :: ZRDORV ! RD/RV
-REAL :: ZRVORD ! RV/RD
-
-
-REAL, DIMENSION(SIZE(PTHM,1)) :: ZMIX1,ZMIX2,ZMIX3
-
-REAL, DIMENSION(SIZE(PTHM,1)) :: ZLUP ! Upward Mixing length from the ground
-
-REAL, DIMENSION(SIZE(PTHM,1)) :: ZDEPTH ! Deepness limit for cloud
-
-INTEGER :: ISV ! Number of scalar variables
-INTEGER :: IKU,IIJU ! array size in k
-INTEGER :: JK,JI,JJ,JSV ! loop counters
-
-LOGICAL, DIMENSION(SIZE(PTHM,1)) :: GTEST,GTESTLCL,GTESTETL
- ! Test if the ascent continue, if LCL or ETL is reached
-LOGICAL :: GLMIX
- ! To choose upward or downward mixing length
-LOGICAL, DIMENSION(SIZE(PTHM,1)) :: GWORK1
-LOGICAL, DIMENSION(SIZE(PTHM,1),SIZE(PTHM,2)) :: GWORK2
-
-
-INTEGER :: ITEST
-
-REAL, DIMENSION(SIZE(PTHM,1)) :: ZRC_UP, ZRI_UP, ZRV_UP, ZWP2, ZRSATW, ZRSATI
-
-LOGICAL, DIMENSION(SIZE(PTHM,1)) :: GTEST_FER
-REAL, DIMENSION(SIZE(PTHM,1)) :: ZPHI,ZALIM_STAR_TOT
-REAL, DIMENSION(SIZE(PTHM,1),SIZE(PTHM,2)) :: ZDTHETASDZ,ZALIM_STAR,ZZDZ,ZZZ
-INTEGER, DIMENSION(SIZE(PTHM,1)) :: IALIM
-
-REAL, DIMENSION(SIZE(PTHM,1)) :: ZTEST,ZDZ,ZWUP_MEAN !
-REAL, DIMENSION(SIZE(PTHM,1)) :: ZCOE,ZWCOE,ZBUCOE
-REAL, DIMENSION(SIZE(PTHM,1)) :: ZDETR_BUO, ZDETR_RT
-REAL, DIMENSION(SIZE(PTHM,1)) :: ZW_MAX ! w**2 max of the updraft
-REAL, DIMENSION(SIZE(PTHM,1)) :: ZZTOP ! Top of the updraft
-REAL, DIMENSION(SIZE(PTHM,1)) :: ZA,ZB,ZQTM,ZQT_UP
-
-REAL :: ZDEPTH_MAX1, ZDEPTH_MAX2 ! control auto-extinction process
-
-REAL :: ZTMAX,ZRMAX, ZEPS ! control value
-
-
-! Thresholds for the perturbation of
-! theta_l and r_t at the first level of the updraft
-
-ZTMAX=2.0
-ZRMAX=1.E-3
-ZEPS=1.E-15
-!------------------------------------------------------------------------
-! INITIALISATION
-
-! Initialisation of the constants
-ZRDORV = XRD / XRV !=0.622
-ZRVORD = (XRV / XRD)
-
-ZDEPTH_MAX1=4500. ! clouds with depth infeRIOr to this value are keeped untouched
-ZDEPTH_MAX2=5000. ! clouds with depth superior to this value are suppressed
-
-! Local variables, internal domain
-! Internal Domain
-
-IKU=SIZE(PTHM,2)
-IIJU =SIZE(PTHM,1)
-!number of scalar variables
-ISV=SIZE(PSVM,3)
-
-! Initialisation of intersesting Level :LCL,ETL,CTL
-KKLCL(:)=KKE
-KKETL(:)=KKE
-KKCTL(:)=KKE
-
-!
-! Initialisation
-!* udraft governing variables
-PEMF(:,:)=0.
-PDETR(:,:)=0.
-PENTR(:,:)=0.
-
-! Initialisation
-!* updraft core variables
-PRV_UP(:,:)=0.
-PRC_UP(:,:)=0.
-
-PW_UP(:,:)=0.
-ZTH_UP(:,:)=0.
-PFRAC_UP(:,:)=0.
-PTHV_UP(:,:)=0.
-
-PBUO_INTEG=0.
-ZBUO =0.
-
-!no ice cloud coded yet
-PRI_UP(:,:)=0.
-PFRAC_ICE_UP(:,:)=0.
-PRSAT_UP(:,:)=PRVM(:,:) ! should be initialised correctly but is (normaly) not used
-
-! Initialisation of environment variables at t-dt
-
-! variables at flux level
-ZTHLM_F(:,:) = MZM_MF(KKA,KKU,KKL,PTHLM(:,:))
-ZRTM_F (:,:) = MZM_MF(KKA,KKU,KKL,PRTM(:,:))
-ZUM_F (:,:) = MZM_MF(KKA,KKU,KKL,PUM(:,:))
-ZVM_F (:,:) = MZM_MF(KKA,KKU,KKL,PVM(:,:))
-ZTKEM_F(:,:) = MZM_MF(KKA,KKU,KKL,PTKEM(:,:))
-
-!DO JSV=1,ISV
-! IF (ONOMIXLG .AND. JSV >= KSV_LGBEG .AND. JSV<= KSV_LGEND) CYCLE
-! ZSVM_F(:,KKB:IKU,JSV) = 0.5*(PSVM(:,KKB:IKU,JSV)+PSVM(:,1:IKU-1,JSV))
-! ZSVM_F(:,1,JSV) = ZSVM_F(:,KKB,JSV)
-!END DO
-
-! Initialisation of updraft characteristics
-PTHL_UP(:,:)=ZTHLM_F(:,:)
-PRT_UP(:,:)=ZRTM_F(:,:)
-PU_UP(:,:)=ZUM_F(:,:)
-PV_UP(:,:)=ZVM_F(:,:)
-PSV_UP(:,:,:)=0.
-!IF (ONOMIXLG .AND. JSV >= KSV_LGBEG .AND. JSV<= KSV_LGEND) then
-! PSV_UP(:,:,:)=ZSVM_F(:,:,:)
-!ENDIF
-
-! Computation or initialisation of updraft characteristics at the KKB level
-! thetal_up,rt_up,thetaV_up, w�,Buoyancy term and mass flux (PEMF)
-
-PTHL_UP(:,KKB)= ZTHLM_F(:,KKB)+MAX(0.,MIN(ZTMAX,(PSFTH(:)/SQRT(ZTKEM_F(:,KKB)))*XALP_PERT))
-PRT_UP(:,KKB) = ZRTM_F(:,KKB)+MAX(0.,MIN(ZRMAX,(PSFRV(:)/SQRT(ZTKEM_F(:,KKB)))*XALP_PERT))
-
-ZQT_UP(:) = PRT_UP(:,KKB)/(1.+PRT_UP(:,KKB))
-ZTHS_UP(:,KKB)=PTHL_UP(:,KKB)*(1.+XLAMBDA*ZQT_UP(:))
-
-ZTHM_F (:,:) = MZM_MF(KKA,KKU,KKL,PTHM (:,:))
-ZPRES_F(:,:) = MZM_MF(KKA,KKU,KKL,PPABSM(:,:))
-ZRHO_F (:,:) = MZM_MF(KKA,KKU,KKL,PRHODREF(:,:))
-ZRVM_F (:,:) = MZM_MF(KKA,KKU,KKL,PRVM(:,:))
-
-! thetav at mass and flux levels
-ZTHVM_F(:,:)=ZTHM_F(:,:)*((1.+ZRVORD*ZRVM_F(:,:))/(1.+ZRTM_F(:,:)))
-ZTHVM(:,:)=PTHM(:,:)*((1.+ZRVORD*PRVM(:,:))/(1.+PRTM(:,:)))
-
-PTHV_UP(:,:)= ZTHVM_F(:,:)
-PRV_UP (:,:)= ZRVM_F (:,:)
-
-ZW_UP2(:,:)=ZEPS
-ZW_UP2(:,KKB) = MAX(0.0001,(1./6.)*ZTKEM_F(:,KKB))
-GTEST = (ZW_UP2(:,KKB) > ZEPS)
-
-! Computation of non conservative variable for the KKB level of the updraft
-! (all or nothing ajustement)
-PRC_UP(:,KKB)=0.
-PRI_UP(:,KKB)=0.
-
-CALL TH_R_FROM_THL_RT_1D(HFRAC_ICE,PFRAC_ICE_UP(:,KKB),ZPRES_F(:,KKB), &
- PTHL_UP(:,KKB),PRT_UP(:,KKB),ZTH_UP(:,KKB), &
- PRV_UP(:,KKB),PRC_UP(:,KKB),PRI_UP(:,KKB),ZRSATW(:),ZRSATI(:))
-
-! compute updraft thevav and buoyancy term at KKB level
-PTHV_UP(:,KKB) = ZTH_UP(:,KKB)*((1+ZRVORD*PRV_UP(:,KKB))/(1+PRT_UP(:,KKB)))
-! compute mean rsat in updraft
-PRSAT_UP(:,KKB) = ZRSATW(:)*(1-PFRAC_ICE_UP(:,KKB)) + ZRSATI(:)*PFRAC_ICE_UP(:,KKB)
-
-!Tout est commente pour tester dans un premier temps la s�paration en deux de la
-! boucle verticale, une pour w et une pour PEMF
-
-ZG_O_THVREF=XG/ZTHVM_F
-
-
-! Definition de l'alimentation au sens de la fermeture de Hourdin et al
-
-ZALIM_STAR(:,:) = 0.
-ZALIM_STAR_TOT(:) = 0. ! <== Normalization of ZALIM_STAR
-IALIM(:) = KKB ! <== Top level of the alimentation layer
-
-DO JK=KKB,KKE-KKL,KKL ! Vertical loop
- ZZDZ(:,JK) = MAX(ZEPS,PZZ(:,JK+KKL)-PZZ(:,JK)) ! <== Delta Z between two flux level
- ZZZ(:,JK) = MAX(0.,0.5*(PZZ(:,JK+KKL)+PZZ(:,JK)) ) ! <== Hight of mass levels
- ZDTHETASDZ(:,JK) = (ZTHVM_F(:,JK)-ZTHVM_F(:,JK+KKL)) ! <== Delta theta_v
-
- WHERE ((ZTHVM_F(:,JK+KKL)<ZTHVM_F(:,JK)) .AND. (ZTHVM_F(:,KKB)>=ZTHVM_F(:,JK)))
- ZALIM_STAR(:,JK) = SQRT(ZZZ(:,JK))*ZDTHETASDZ(:,JK)/ZZDZ(:,JK)
- ZALIM_STAR_TOT(:) = ZALIM_STAR_TOT(:)+ZALIM_STAR(:,JK)*ZZDZ(:,JK)
- IALIM(:) = JK
- ENDWHERE
-ENDDO
-
-! Normalization of ZALIM_STAR
-DO JK=KKB,KKE-KKL,KKL ! Vertical loop
- WHERE (ZALIM_STAR_TOT > ZEPS)
- ZALIM_STAR(:,JK) = ZALIM_STAR(:,JK)/ZALIM_STAR_TOT(:)
- ENDWHERE
-ENDDO
-ZALIM_STAR_TOT(:) = 0.
-
-
-! --------- END of alimentation calculation ---------------------------------------
-
-
-!--------------------------------------------------------------------------
-
-! 3. Vertical ascending loop
-! -----------------------
-!
-! If GTEST = T the updraft starts from the KKB level and stops when GTEST becomes F
-!
-!
-GTESTLCL(:)=.FALSE.
-GTESTETL(:)=.FALSE.
-
-! Loop on vertical level to compute W
-
-ZW_MAX(:) = 0.
-ZZTOP(:) = 0.
-ZPHI(:) = 0.
-
-
-DO JK=KKB,KKE-KKL,KKL
-
-! IF the updraft top is reached for all column, stop the loop on levels
-
-! ITEST=COUNT(GTEST)
-! IF (ITEST==0) CYCLE
-
-! Computation of entrainment and detrainment with KF90
-! parameterization in clouds and LR01 in subcloud layer
-
-
-! to find the LCL (check if JK is LCL or not)
-
- WHERE ((PRC_UP(:,JK)+PRI_UP(:,JK)>0.).AND.(.NOT.(GTESTLCL)))
- KKLCL(:) = JK
- GTESTLCL(:)=.TRUE.
- ENDWHERE
-
-
-! COMPUTE PENTR and PDETR at mass level JK
-
-
-! Buoyancy is computed on "flux" levels where updraft variables are known
-
- ! Compute theta_v of updraft at flux level JK
-
- ZRC_UP(:) = PRC_UP(:,JK)
- ZRI_UP(:) = PRI_UP(:,JK) ! guess
- ZRV_UP(:) = PRV_UP(:,JK)
- ZBUO (:,JK) = ZG_O_THVREF(:,JK)*(PTHV_UP(:,JK) - ZTHVM_F(:,JK))
- PBUO_INTEG(:,JK) = ZBUO(:,JK)*(PZZ(:,JK+KKL)-PZZ(:,JK))
-
- ZDZ(:) = MAX(ZEPS,PZZ(:,JK+KKL)-PZZ(:,JK))
- ZTEST(:) = XA1*ZBUO(:,JK) - XB*ZW_UP2(:,JK)
-
- ZCOE(:) = ZDZ(:)
- WHERE (ZTEST(:)>0.)
- ZCOE(:) = ZDZ(:)/(1.+ XBETA1)
- ENDWHERE
-
-! Calcul de la vitesse
-
- ZWCOE(:) = (1.-XB*ZCOE(:))/(1.+XB*ZCOE(:))
- ZBUCOE(:) = 2.*ZCOE(:)/(1.+XB*ZCOE(:))
-
- ZW_UP2(:,JK+KKL) = MAX(ZEPS,ZW_UP2(:,JK)*ZWCOE(:) + XA1*ZBUO(:,JK)*ZBUCOE(:) )
- ZW_MAX(:) = MAX(ZW_MAX(:), SQRT(ZW_UP2(:,JK+KKL)))
- ZWUP_MEAN(:) = MAX(ZEPS,0.5*(ZW_UP2(:,JK+KKL)+ZW_UP2(:,JK)))
-
-! Entrainement et detrainement
-
- PENTR(:,JK) = MAX(0.,(XBETA1/(1.+XBETA1))*(XA1*ZBUO(:,JK)/ZWUP_MEAN(:)-XB))
-
- ZDETR_BUO(:) = MAX(0., -(XBETA1/(1.+XBETA1))*XA1*ZBUO(:,JK)/ZWUP_MEAN(:))
- ZDETR_RT(:) = XC*SQRT(MAX(0.,(PRT_UP(:,JK) - ZRTM_F(:,JK))) / MAX(ZEPS,ZRTM_F(:,JK)) / ZWUP_MEAN(:))
- PDETR(:,JK) = ZDETR_RT(:)+ZDETR_BUO(:)
-
-
-! If the updraft did not stop, compute cons updraft characteritics at jk+1
- WHERE(GTEST)
- ZZTOP(:) = MAX(ZZTOP(:),PZZ(:,JK+KKL))
- ZMIX2(:) = (PZZ(:,JK+KKL)-PZZ(:,JK))*PENTR(:,JK) !&
- ZMIX3(:) = (PZZ(:,JK+KKL)-PZZ(:,JK))*PDETR(:,JK) !&
-
- ZQTM(:) = PRTM(:,JK)/(1.+PRTM(:,JK))
- ZTHSM(:,JK) = PTHLM(:,JK)*(1.+XLAMBDA*ZQTM(:))
- ZTHS_UP(:,JK+KKL)=(ZTHS_UP(:,JK)*(1.-0.5*ZMIX2(:)) + ZTHSM(:,JK)*ZMIX2(:)) &
- /(1.+0.5*ZMIX2(:))
- PRT_UP(:,JK+KKL)=(PRT_UP (:,JK)*(1.-0.5*ZMIX2(:)) + PRTM(:,JK)*ZMIX2(:)) &
- /(1.+0.5*ZMIX2(:))
- ZQT_UP(:) = PRT_UP(:,JK+KKL)/(1.+PRT_UP(:,JK+KKL))
- PTHL_UP(:,JK+KKL)=ZTHS_UP(:,JK+KKL)/(1.+XLAMBDA*ZQT_UP(:))
- ENDWHERE
-
-
- IF(OMIXUV) THEN
- IF(JK/=KKB) THEN
- WHERE(GTEST)
- PU_UP(:,JK+KKL) = (PU_UP (:,JK)*(1-0.5*ZMIX2(:)) + PUM(:,JK)*ZMIX2(:)+ &
- 0.5*XPRES_UV*(PZZ(:,JK+KKL)-PZZ(:,JK))*&
- ((PUM(:,JK+KKL)-PUM(:,JK))/PDZZ(:,JK+KKL)+&
- (PUM(:,JK)-PUM(:,JK-KKL))/PDZZ(:,JK)) ) &
- /(1+0.5*ZMIX2(:))
- PV_UP(:,JK+KKL) = (PV_UP (:,JK)*(1-0.5*ZMIX2(:)) + PVM(:,JK)*ZMIX2(:)+ &
- 0.5*XPRES_UV*(PZZ(:,JK+KKL)-PZZ(:,JK))*&
- ((PVM(:,JK+KKL)-PVM(:,JK))/PDZZ(:,JK+KKL)+&
- (PVM(:,JK)-PVM(:,JK-KKL))/PDZZ(:,JK)) ) &
- /(1+0.5*ZMIX2(:))
- ENDWHERE
- ELSE
- WHERE(GTEST)
- PU_UP(:,JK+KKL) = (PU_UP (:,JK)*(1-0.5*ZMIX2(:)) + PUM(:,JK)*ZMIX2(:)+ &
- 0.5*XPRES_UV*(PZZ(:,JK+KKL)-PZZ(:,JK))*&
- ((PUM(:,JK+KKL)-PUM(:,JK))/PDZZ(:,JK+KKL)) ) &
- /(1+0.5*ZMIX2(:))
- PV_UP(:,JK+KKL) = (PV_UP (:,JK)*(1-0.5*ZMIX2(:)) + PVM(:,JK)*ZMIX2(:)+ &
- 0.5*XPRES_UV*(PZZ(:,JK+KKL)-PZZ(:,JK))*&
- ((PVM(:,JK+KKL)-PVM(:,JK))/PDZZ(:,JK+KKL)) ) &
- /(1+0.5*ZMIX2(:))
- ENDWHERE
-
- ENDIF
- ENDIF
-! DO JSV=1,ISV
-! IF (ONOMIXLG .AND. JSV >= KSV_LGBEG .AND. JSV<= KSV_LGEND) CYCLE
-! WHERE(GTEST)
-! PSV_UP(:,JK+KKL,JSV) = (PSV_UP (:,JK,JSV)*(1-0.5*ZMIX2(:)) + &
-! PSVM(:,JK,JSV)*ZMIX2(:)) /(1+0.5*ZMIX2(:))
-! ENDWHERE
-! ENDDO
-
-
-! Compute non cons. var. at level JK+KKL
- ZRC_UP(:)=PRC_UP(:,JK) ! guess = level just below
- ZRI_UP(:)=PRI_UP(:,JK) ! guess = level just below
- ZRV_UP(:)=PRV_UP(:,JK)
- CALL TH_R_FROM_THL_RT_1D(HFRAC_ICE,PFRAC_ICE_UP(:,JK+KKL),ZPRES_F(:,JK+KKL), &
- PTHL_UP(:,JK+KKL),PRT_UP(:,JK+KKL),ZTH_UP(:,JK+KKL), &
- ZRV_UP(:),ZRC_UP(:),ZRI_UP(:),ZRSATW(:),ZRSATI(:))
- WHERE(GTEST)
- ZT_UP(:) = ZTH_UP(:,JK+KKL)*PEXNM(:,JK+KKL)
- ZCP(:) = XCPD + XCL * ZRC_UP(:)
- ZLVOCPEXN(:)=(XLVTT + (XCPV-XCL) * (ZT_UP(:)-XTT) ) / ZCP(:) / PEXNM(:,JK+KKL)
- PRC_UP(:,JK+KKL)=MIN(0.5E-3,ZRC_UP(:)) ! On ne peut depasser 0.5 g/kg (autoconversion donc elimination !)
- PTHL_UP(:,JK+KKL) = PTHL_UP(:,JK+KKL)+ZLVOCPEXN(:)*(ZRC_UP(:)-PRC_UP(:,JK+KKL))
- PRV_UP(:,JK+KKL)=ZRV_UP(:)
- PRI_UP(:,JK+KKL)=ZRI_UP(:)
- PRT_UP(:,JK+KKL) = PRC_UP(:,JK+KKL) + PRV_UP(:,JK+KKL)
- PRSAT_UP(:,JK+KKL) = ZRSATW(:)*(1-PFRAC_ICE_UP(:,JK+KKL)) + ZRSATI(:)*PFRAC_ICE_UP(:,JK+KKL)
- ENDWHERE
-
-
-! Compute the updraft theta_v, buoyancy and w**2 for level JK+1
- WHERE(GTEST)
-! PTHV_UP(:,JK+KKL) = ZTH_UP(:,JK+KKL)*((1+ZRVORD*PRV_UP(:,JK+KKL))/(1+PRT_UP(:,JK+KKL)))
- PTHV_UP(:,JK+KKL) = ZTH_UP(:,JK+KKL)*(1.+0.608*PRV_UP(:,JK+KKL) - PRC_UP(:,JK+KKL))
- ENDWHERE
-
-
-! Test if the updraft has reach the ETL
- GTESTETL(:)=.FALSE.
- WHERE (GTEST.AND.(PBUO_INTEG(:,JK)<=0.))
- KKETL(:) = JK+KKL
- GTESTETL(:)=.TRUE.
- ENDWHERE
-
-! Test is we have reached the top of the updraft
-
- WHERE (GTEST.AND.((ZW_UP2(:,JK+KKL)<=ZEPS)))
- ZW_UP2(:,JK+KKL)=ZEPS
- GTEST(:)=.FALSE.
- PTHL_UP(:,JK+KKL)=ZTHLM_F(:,JK+KKL)
- PRT_UP(:,JK+KKL)=ZRTM_F(:,JK+KKL)
- PRC_UP(:,JK+KKL)=0.
- PRI_UP(:,JK+KKL)=0.
- PRV_UP(:,JK+KKL)=0.
- PTHV_UP(:,JK+KKL)=ZTHVM_F(:,JK+KKL)
- PFRAC_UP(:,JK+KKL)=0.
- KKCTL(:)=JK+KKL
- ENDWHERE
-
-ENDDO
-
-! Closure assumption for mass flux at KKB+1 level (Mass flux is supposed to be 0 at KKB level !)
-! Hourdin et al 2002 formulation
-
-
-ZZTOP(:) = MAX(ZZTOP(:),ZEPS)
-
-DO JK=KKB+KKL,KKE-KKL,KKL ! Vertical loop
- WHERE(JK<=IALIM)
- ZALIM_STAR_TOT(:) = ZALIM_STAR_TOT(:) + ZALIM_STAR(:,JK)*ZALIM_STAR(:,JK)*ZZDZ(:,JK)/PRHODREF(:,JK)
- ENDWHERE
-ENDDO
-
-WHERE (ZALIM_STAR_TOT*ZZTOP > ZEPS)
- ZPHI(:) = ZW_MAX(:)/(XR*ZZTOP(:)*ZALIM_STAR_TOT(:))
-ENDWHERE
-
-GTEST(:) = .TRUE.
-PEMF(:,KKB+KKL) = ZPHI(:)*ZZDZ(:,KKB)*ZALIM_STAR(:,KKB)
-! Updraft fraction must be smaller than XFRAC_UP_MAX
-PFRAC_UP(:,KKB+KKL)=PEMF(:,KKB+KKL)/(SQRT(ZW_UP2(:,KKB+KKL))*ZRHO_F(:,KKB+KKL))
-PFRAC_UP(:,KKB+KKL)=MIN(XFRAC_UP_MAX,PFRAC_UP(:,KKB+KKL))
-PEMF(:,KKB+KKL) = ZRHO_F(:,KKB+KKL)*PFRAC_UP(:,KKB+KKL)*SQRT(ZW_UP2(:,KKB+KKL))
-
-DO JK=KKB+KKL,KKE-KKL,KKL ! Vertical loop
-
- GTEST = (ZW_UP2(:,JK) > ZEPS)
-
- WHERE (GTEST)
- WHERE(JK<IALIM)
- PEMF(:,JK+KKL) = MAX(0.,PEMF(:,JK) + ZPHI(:)*ZZDZ(:,JK)*(PENTR(:,JK) - PDETR(:,JK)))
- ELSEWHERE
- ZMIX1(:)=ZZDZ(:,JK)*(PENTR(:,JK)-PDETR(:,JK))
- PEMF(:,JK+KKL)=PEMF(:,JK)*EXP(ZMIX1(:))
- ENDWHERE
-
-! Updraft fraction must be smaller than XFRAC_UP_MAX
- PFRAC_UP(:,JK+KKL)=PEMF(:,JK+KKL)/(SQRT(ZW_UP2(:,JK+KKL))*ZRHO_F(:,JK+KKL))
- PFRAC_UP(:,JK+KKL)=MIN(XFRAC_UP_MAX,PFRAC_UP(:,JK+KKL))
- PEMF(:,JK+KKL) = ZRHO_F(:,JK+KKL)*PFRAC_UP(:,JK+KKL)*SQRT(ZW_UP2(:,JK+KKL))
- ENDWHERE
-
-ENDDO
-
-PW_UP(:,:)=SQRT(ZW_UP2(:,:))
-PEMF(:,KKB) =0.
-
-! Limits the shallow convection scheme when cloud heigth is higher than 3000m.
-! To do this, mass flux is multiplied by a coefficient decreasing linearly
-! from 1 (for clouds of 3000m of depth) to 0 (for clouds of 4000m of depth).
-! This way, all MF fluxes are diminished by this amount.
-! Diagnosed cloud fraction is also multiplied by the same coefficient.
-!
-DO JI=1,SIZE(PTHM,1)
- PDEPTH(JI) = MAX(0., PZZ(JI,KKCTL(JI)) - PZZ(JI,KKLCL(JI)) )
-END DO
-
-GWORK1(:)= (GTESTLCL(:) .AND. (PDEPTH(:) > ZDEPTH_MAX1) )
-GWORK2(:,:) = SPREAD( GWORK1(:), DIM=2, NCOPIES=IKU )
-ZCOEF(:,:) = SPREAD( (1.-(PDEPTH(:)-ZDEPTH_MAX1)/(ZDEPTH_MAX2-ZDEPTH_MAX1)), DIM=2, NCOPIES=IKU)
-ZCOEF=MIN(MAX(ZCOEF,0.),1.)
-WHERE (GWORK2)
- PEMF(:,:) = PEMF(:,:) * ZCOEF(:,:)
- PFRAC_UP(:,:) = PFRAC_UP(:,:) * ZCOEF(:,:)
-ENDWHERE
-
-
-END SUBROUTINE COMPUTE_UPDRAFT_RAHA
diff --git a/src/MNH/compute_updraft_rhcj10.f90 b/src/MNH/compute_updraft_rhcj10.f90
index 3a9614ec7..8b3f6a5e2 100644
--- a/src/MNH/compute_updraft_rhcj10.f90
+++ b/src/MNH/compute_updraft_rhcj10.f90
@@ -124,7 +124,6 @@ SUBROUTINE COMPUTE_UPDRAFT_RHCJ10(KKA,KKB,KKE,KKU,KKL,HFRAC_ICE, &
!! ------
!! Y. Bouteloup (2012)
!! R. Honert Janv 2013 ==> corection of some bugs
-!! 10/2016 (R.Honnert and S.Riette) : Improvement of EDKF and adaptation to the grey zone
!! --------------------------------------------------------------------------
! WARNING ==> This updraft is not yet ready to use scalar variables
@@ -391,8 +390,8 @@ ENDDO
GLMIX=.TRUE.
ZTKEM_F(:,KKB)=0.
-CALL COMPUTE_BL89_ML(KKA,KKB,KKE,KKU,KKL,PDZZ,ZTKEM_F(:,KKB),ZG_O_THVREF(:,KKB), &
- ZTHVM_F,KKB,GLMIX,.TRUE.,ZLUP)
+
+CALL COMPUTE_BL89_ML(KKA,KKB,KKE,KKU,KKL,PDZZ,ZTKEM_F,ZG_O_THVREF,ZTHVM_F,KKB,GLMIX,ZLUP)
ZLUP(:)=MAX(ZLUP(:),1.E-10)
! Compute Buoyancy flux at the ground
diff --git a/src/MNH/default_desfmn.f90 b/src/MNH/default_desfmn.f90
index 91c259277..6ec7ad809 100644
--- a/src/MNH/default_desfmn.f90
+++ b/src/MNH/default_desfmn.f90
@@ -999,6 +999,8 @@ IF (KMI == 1) THEN
LSEDIC = .FALSE.
LCONVHG = .FALSE.
CSEDIM = 'SPLI'
+ LDEPOSC = .FALSE.
+ XVDEPOSC= 0.02 ! 2 cm/s
END IF
!
!-------------------------------------------------------------------------------
@@ -1052,8 +1054,6 @@ XA1 = 2./3.
XB = 0.002
XC = 0.012
XBETA1 = 0.9
-XR = 2.
-XLAMBDA= 0.
!
!-------------------------------------------------------------------------------
!
diff --git a/src/MNH/shallow_mf.f90 b/src/MNH/shallow_mf.f90
index f30b6b2b0..1bf2aaa96 100644
--- a/src/MNH/shallow_mf.f90
+++ b/src/MNH/shallow_mf.f90
@@ -161,9 +161,7 @@ END MODULE MODI_SHALLOW_MF
!! S.Riette DUAL case
!! S. Riette Jan 2012: support for both order of vertical levels
!! R.Honnert 07/2012 : elemnts of Rio according to Bouteloup
-!! R.Honnert 07/2012 : MF gray zone
-!! R.Honnert 10/2016 : SURF=gray zone initilisation + EDKF
-!! R.Honnert 10/2016 : Update with Arome
+!! R.Honnert 07/2012 : EDKF gray zone
!! --------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -176,10 +174,8 @@ USE MODD_PARAM_MFSHALL_n
USE MODI_THL_RT_FROM_TH_R_MF
USE MODI_COMPUTE_UPDRAFT
USE MODI_COMPUTE_UPDRAFT_RHCJ10
-USE MODI_COMPUTE_UPDRAFT_RAHA
USE MODI_COMPUTE_UPDRAFT_HRIO
USE MODI_MF_TURB
-USE MODI_MF_TURB_EXPL
USE MODI_MF_TURB_GREYZONE
USE MODI_COMPUTE_MF_CLOUD
USE MODI_COMPUTE_FRAC_ICE
@@ -288,8 +284,6 @@ INTEGER :: IKE ! uppest atmosphere physical index
REAL, DIMENSION(SIZE(PTHM,1),SIZE(PTHM,2)) :: ZG_O_THVREF,PTHVREF
REAL, DIMENSION(SIZE(PTHM,1)) :: ZRESOL_NORM, ZRESOL_GRID,& ! normalized grid
ZLUP, ZPLAW
- ! Test if the ascent continue, if LCL or ETL is reached
-LOGICAL :: GLMIX
INTEGER :: JI,JJ,JK ! loop counter
!------------------------------------------------------------------------
@@ -300,9 +294,8 @@ IKB=KKA+KKL*JPVEXT
IKE=KKU-KKL*JPVEXT
! updraft governing variables
-IF (HMF_UPDRAFT == 'EDKF' .OR. HMF_UPDRAFT == 'HRIO' .OR. &
- HMF_UPDRAFT == 'RHCJ' .OR. HMF_UPDRAFT == 'BOUT' .OR. &
- HMF_UPDRAFT == 'SURF' ) THEN
+IF (HMF_UPDRAFT == 'EDKF'.OR. HMF_UPDRAFT == 'HRIO' .OR. &
+ HMF_UPDRAFT == 'RHCJ'.OR. HMF_UPDRAFT == 'BOUT' ) THEN
PENTR = 1.E20
PDETR = 1.E20
PEMF = 1.E20
@@ -324,10 +317,9 @@ ZTHVM(:,:) = PTHM(:,:)*((1.+XRV / XRD *PRM(:,:,1))/(1.+ZRTM(:,:)))
!!! 2. Compute updraft
!!! ---------------
!
-IF (HMF_UPDRAFT == 'EDKF' .OR. HMF_UPDRAFT == 'BOUT' .OR. HMF_UPDRAFT == 'SURF') THEN
+IF (HMF_UPDRAFT == 'EDKF' .OR. HMF_UPDRAFT == 'BOUT') THEN
GENTR_DETR = .TRUE.
- CALL COMPUTE_UPDRAFT(KKA,IKB,IKE,KKU,KKL,HFRAC_ICE,HMF_UPDRAFT,&
- GENTR_DETR, OMIXUV, &
+ CALL COMPUTE_UPDRAFT(KKA,IKB,IKE,KKU,KKL,HFRAC_ICE,GENTR_DETR,OMIXUV,&
ONOMIXLG,KSV_LGBEG,KSV_LGEND, &
PZZ,PDZZ, &
PSFTH,PSFRV,PPABSM,PRHODREF, &
@@ -337,36 +329,18 @@ IF (HMF_UPDRAFT == 'EDKF' .OR. HMF_UPDRAFT == 'BOUT' .OR. HMF_UPDRAFT == 'SURF')
PTHV_UP, PW_UP, PU_UP, PV_UP, ZSV_UP, &
PFRAC_UP,ZFRAC_ICE_UP,ZRSAT_UP,PEMF,PDETR,&
PENTR,ZBUO_INTEG,KKLCL,KKETL,KKCTL,ZDEPTH )
-ELSEIF (HMF_UPDRAFT == 'RHCJ') THEN
- CALL COMPUTE_UPDRAFT_RHCJ10(KKA,IKB,IKE,KKU,KKL,HFRAC_ICE, &
- GENTR_DETR,OMIXUV, &
- ONOMIXLG,KSV_LGBEG,KSV_LGEND, &
- PZZ,PDZZ, &
- PSFTH,PSFRV, &
- PPABSM,PRHODREF,PUM,PVM,PTKEM, &
- PTHM,PRM(:,:,1),ZTHLM,ZRTM, &
- PSVM, PTHL_UP,PRT_UP, &
- PRV_UP,PRC_UP,PRI_UP,PTHV_UP, &
- PW_UP, PU_UP, PV_UP, ZSV_UP, &
- PFRAC_UP,ZFRAC_ICE_UP,ZRSAT_UP, &
- PEMF,PDETR,PENTR, &
- ZBUO_INTEG,KKLCL,KKETL,KKCTL, &
- ZDEPTH )
-ELSEIF (HMF_UPDRAFT == 'RAHA') THEN
- CALL COMPUTE_UPDRAFT_RAHA(KKA,IKB,IKE,KKU,KKL,HFRAC_ICE, &
- GENTR_DETR,OMIXUV, &
+ELSEIF (HMF_UPDRAFT == 'RHCJ') THEN
+ GENTR_DETR = .TRUE.
+ CALL COMPUTE_UPDRAFT_RHCJ10(KKA,IKB,IKE,KKU,KKL,HFRAC_ICE,GENTR_DETR,OMIXUV,&
ONOMIXLG,KSV_LGBEG,KSV_LGEND, &
PZZ,PDZZ, &
- PSFTH,PSFRV, &
- PPABSM,PRHODREF,PUM,PVM,PTKEM, &
- PEXNM,PTHM,PRM(:,:,1),ZTHLM,ZRTM, &
- PSVM,PTHL_UP,PRT_UP, &
- PRV_UP,PRC_UP,PRI_UP, PTHV_UP, &
- PW_UP, PU_UP, PV_UP, ZSV_UP, &
- PFRAC_UP,ZFRAC_ICE_UP,ZRSAT_UP, &
- PEMF,PDETR,PENTR, &
- ZBUO_INTEG,KKLCL,KKETL,KKCTL, &
- ZDEPTH )
+ PSFTH,PSFRV,PPABSM,PRHODREF, &
+ PUM,PVM,PTKEM, &
+ PTHM,PRM(:,:,1),ZTHLM,ZRTM,PSVM, &
+ PTHL_UP,PRT_UP,PRV_UP,PRC_UP,PRI_UP, &
+ PTHV_UP, PW_UP, PU_UP, PV_UP, ZSV_UP, &
+ PFRAC_UP,ZFRAC_ICE_UP,ZRSAT_UP,PEMF,PDETR,&
+ PENTR,ZBUO_INTEG,KKLCL,KKETL,KKCTL,ZDEPTH )
ELSEIF (HMF_UPDRAFT == 'DUAL') THEN
!Updraft characteristics are already computed and received by interface
ELSEIF (HMF_UPDRAFT == 'HRIO') THEN
@@ -411,9 +385,7 @@ CALL COMPUTE_MF_CLOUD(KKA,IKB,IKE,KKU,KKL,KRR,KRRL,KRRI,&
!!! ------------------------------------------------------------------------
!
ZEMF_O_RHODREF=PEMF/PRHODREF
-IF(HMF_UPDRAFT == 'EDKF' .OR. HMF_UPDRAFT == 'RHCJ'.OR. HMF_UPDRAFT == 'BOUT' &
- .OR. HMF_UPDRAFT == 'SURF') THEN
- IF ( PIMPL_MF > 1.E-10 ) THEN
+ IF(HMF_UPDRAFT == 'EDKF' .OR. HMF_UPDRAFT == 'RHCJ'.OR. HMF_UPDRAFT == 'BOUT') THEN
CALL MF_TURB(KKA, IKB, IKE, KKU, KKL, OMIXUV, &
ONOMIXLG,KSV_LGBEG,KSV_LGEND, &
PIMPL_MF, PTSTEP, &
@@ -424,14 +396,6 @@ IF(HMF_UPDRAFT == 'EDKF' .OR. HMF_UPDRAFT == 'RHCJ'.OR. HMF_UPDRAFT == 'BOUT' &
ZEMF_O_RHODREF,PTHL_UP,PTHV_UP,PRT_UP,PU_UP,PV_UP,ZSV_UP,&
PFLXZTHMF,PFLXZTHVMF,PFLXZRMF,PFLXZUMF,PFLXZVMF, &
ZFLXZSVMF )
-ELSE
- CALL MF_TURB_EXPL(KKA, IKB, IKE, KKU, KKL, OMIXUV, &
- PRHODJ, &
- ZTHLM,ZTHVM,ZRTM,PUM,PVM, &
- PDTHLDT_MF,PDRTDT_MF,PDUDT_MF,PDVDT_MF, &
- ZEMF_O_RHODREF,PTHL_UP,PTHV_UP,PRT_UP,PU_UP,PV_UP, &
- PFLXZTHMF,PFLXZTHVMF,PFLXZRMF,PFLXZUMF,PFLXZVMF)
-ENDIF
ELSEIF (HMF_UPDRAFT == 'HRIO') THEN
CALL MF_TURB_GREYZONE(KKA, IKB, IKE, KKU, KKL,OMIXUV, &
ONOMIXLG,KSV_LGBEG,KSV_LGEND, &
@@ -457,8 +421,7 @@ ENDIF
PTHVREF(:,JK)=RESHAPE(XTHVREF(:,:,JK),(/SIZE(PTHM,1)*SIZE(PTHM,2)/) )
ENDDO
ZG_O_THVREF=XG/PTHVREF
- GLMIX=.TRUE.
- CALL COMPUTE_BL89_ML(KKA,IKB,IKE,KKU,KKL,PDZZ,PTKEM(:,IKB) ,ZG_O_THVREF(:,IKB),ZTHVM,IKB,GLMIX,.TRUE.,ZLUP)
+ CALL COMPUTE_BL89_ML(KKA,IKB,IKE,KKU,KKL,PDZZ,PTKEM,ZG_O_THVREF,ZTHVM,IKB,.TRUE.,ZLUP)
!! calcul de Dx/(h+hc)
DO JI=1,SIZE(XDXHAT)
DO JJ=1,SIZE(XDYHAT)
--
GitLab