From e086c6c9f09b19e4744b98773cb91d968fdb321b Mon Sep 17 00:00:00 2001
From: Gaelle TANGUY <gaelle.tanguy@meteo.fr>
Date: Tue, 15 Nov 2016 10:50:09 +0100
Subject: [PATCH] R.Honnert 15/11/2016 : Improvement of EDKF and adaptation to
the grey zone
---
src/MNH/compute_bl89_ml.f90 | 128 ++++--
src/MNH/compute_entr_detr.f90 | 659 ++++++++++++++--------------
src/MNH/compute_updraft.f90 | 91 ++--
src/MNH/compute_updraft_hrio.f90 | 78 +---
src/MNH/compute_updraft_raha.f90 | 667 +++++++++++++++++++++++++++++
src/MNH/compute_updraft_rhcj10.f90 | 5 +-
src/MNH/mf_turb_expl.f90 | 227 ++++++++++
src/MNH/modd_param_mfshalln.f90 | 15 +-
src/MNH/modn_param_mfshalln.f90 | 19 +-
src/MNH/shallow_mf.f90 | 71 ++-
10 files changed, 1450 insertions(+), 510 deletions(-)
create mode 100644 src/MNH/compute_updraft_raha.f90
create mode 100644 src/MNH/mf_turb_expl.f90
diff --git a/src/MNH/compute_bl89_ml.f90 b/src/MNH/compute_bl89_ml.f90
index 4baf16acd..4c5124602 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, &
- PTKEM2D,PG_O_THVREF2D,PVPT,KK,OUPORDN,PLWORK)
+ PTKEM_DEP,PG_O_THVREF,PVPT,KK,OUPORDN,OFLUX,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
-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
+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
END SUBROUTINE COMPUTE_BL89_ML
END INTERFACE
!
END MODULE MODI_COMPUTE_BL89_ML
-
-
-
-! ###################################################################
+! ######spl
SUBROUTINE COMPUTE_BL89_ML(KKA,KKB,KKE,KKU,KKL,PDZZ2D, &
- PTKEM2D,PG_O_THVREF2D,PVPT,KK,OUPORDN,PLWORK)
+ PTKEM_DEP,PG_O_THVREF,PVPT,KK,OUPORDN,OFLUX,PLWORK)
+
! ###################################################################
!!
!! COMPUTE_BL89_ML routine to:
@@ -52,6 +52,7 @@ 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
!!
!-------------------------------------------------------------------------------
!
@@ -82,35 +83,36 @@ 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
-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
+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
! 0.2 Local variable
!
-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)) :: 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),SIZE(PTKEM2D,2)) :: ZDELTVPT,ZHLVPT
+REAL, DIMENSION(SIZE(PVPT,1),SIZE(PVPT,2)) :: ZDELTVPT,ZHLVPT
!Virtual Potential Temp at Half level and DeltaThv between
- !2 levels
+ !2 mass 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(PTKEM2D,1)
+IIJU=SIZE(PVPT,1)
!
ZDELTVPT(:,:)=DZM_MF(KKA,KKU,KKL,PVPT(:,:))
ZDELTVPT(:,KKA)=0.
@@ -120,6 +122,11 @@ 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
@@ -127,27 +134,58 @@ ZHLVPT(:,:)=MZM_MF(KKA,KKU,KKL,PVPT(:,:))
!
IF (OUPORDN.EQV..TRUE.) THEN
- ZINTE(:)=PTKEM2D(:,KK)
+ ZINTE(:)=PTKEM_DEP(:)
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_THVREF2D(J1D,KK) * &
- (ZHLVPT(J1D,JKK) - PVPT(J1D,KK))) * PDZZ2D(J1D,JKK) !particle keeps its temperature
+ ZPOTE(J1D) = ZTEST0*(PG_O_THVREF(J1D) * &
+ (ZHLVPT(J1D,JKK) - ZVPT_DEP(J1D))) * 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_THVREF2D(J1D,KK) * &
- ( PVPT(J1D,JKK-KKL) - PVPT(J1D,KK) ) &
+ ZLWORK2(J1D)= ( - PG_O_THVREF(J1D) * &
+ ( PVPT(J1D,JKK-KKL) - ZVPT_DEP(J1D) ) &
+ SQRT (ABS( &
- ( PG_O_THVREF2D(J1D,KK) * (PVPT(J1D,JKK-KKL) - PVPT(J1D,KK)) )**2 &
- + 2. * ZINTE(J1D) * PG_O_THVREF2D(J1D,KK) &
+ ( PG_O_THVREF(J1D) * (PVPT(J1D,JKK-KKL) - ZVPT_DEP(J1D)) )**2 &
+ + 2. * ZINTE(J1D) * PG_O_THVREF(J1D) &
* ZDELTVPT(J1D,JKK) / PDZZ2D(J1D,JKK) )) ) / &
- ( PG_O_THVREF2D(J1D,KK) * ZDELTVPT(J1D,JKK) / PDZZ2D(J1D,JKK) )
+ ( PG_O_THVREF(J1D) * ZDELTVPT(J1D,JKK) / PDZZ2D(J1D,JKK) )
!
PLWORK(J1D)=PLWORK(J1D)+ZTEST0*(ZTEST*ZLWORK1(J1D)+ &
(1-ZTEST)*ZLWORK2(J1D))
@@ -162,7 +200,13 @@ ENDIF
!
IF (OUPORDN.EQV..FALSE.) THEN
- ZINTE(:)=PTKEM2D(:,KK)
+ IF(OFLUX) THEN
+ WRITE(*,*) ' STOP'
+ WRITE(*,*) ' OFLUX OPTION NOT CODED FOR DOWNWARD MIXING LENGTH'
+ CALL ABORT
+ STOP
+ ENDIF
+ ZINTE(:)=PTKEM_DEP(:)
PLWORK=0.
ZTESTM=1.
DO JKK=KK,KKB,-KKL
@@ -170,18 +214,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_THVREF2D(J1D,KK) * &
+ ZPOTE(J1D) = -ZTEST0*(PG_O_THVREF(J1D) * &
(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_THVREF2D(J1D,KK) * &
+ ZLWORK2(J1D)= ( + PG_O_THVREF(J1D) * &
( PVPT(J1D,JKK) - PVPT(J1D,KK) ) &
+ SQRT (ABS( &
- ( PG_O_THVREF2D(J1D,KK) * (PVPT(J1D,JKK) - PVPT(J1D,KK)) )**2 &
- + 2. * ZINTE(J1D) * PG_O_THVREF2D(J1D,KK) &
+ ( PG_O_THVREF(J1D) * (PVPT(J1D,JKK) - PVPT(J1D,KK)) )**2 &
+ + 2. * ZINTE(J1D) * PG_O_THVREF(J1D) &
* ZDELTVPT(J1D,JKK) / PDZZ2D(J1D,JKK) )) ) / &
- ( PG_O_THVREF2D(J1D,KK) * ZDELTVPT(J1D,JKK) / PDZZ2D(J1D,JKK) )
+ ( PG_O_THVREF(J1D) * 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 fda4d89f4..11d5fce5f 100644
--- a/src/MNH/compute_entr_detr.f90
+++ b/src/MNH/compute_entr_detr.f90
@@ -8,45 +8,56 @@
!
INTERFACE
!
- 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)
+ 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)
-!
-!
-!
-INTEGER, INTENT(IN) :: KK ! near ground physical index
+!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(INOUT) :: OTEST
-LOGICAL,DIMENSION(:),INTENT(INOUT) :: OTESTLCL !test of condensation
-CHARACTER*1,INTENT(IN) :: HFRAC_ICE
-REAL, DIMENSION(:) ,INTENT(IN) :: PFRAC_ICE
-
+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
!
! prognostic variables at t- deltat
-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
-
!
+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
+
!
! 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) :: 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
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(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
+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
!
!
END SUBROUTINE COMPUTE_ENTR_DETR
@@ -56,11 +67,16 @@ END INTERFACE
END MODULE MODI_COMPUTE_ENTR_DETR
! ######spl
SUBROUTINE COMPUTE_ENTR_DETR(KK,KKB,KKE,KKL,OTEST,OTESTLCL,&
- HFRAC_ICE,PFRAC_ICE,PPABSM,PZZ,PDZZ,&
- PTHVM,PTHLM,PRTM,PW_UP2,&
+ 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,PRC_MIX,PRI_MIX, &
- PENTR,PDETR,PBUO_INTEG)
+ 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)
! #############################################################
!!
@@ -98,6 +114,12 @@ 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
@@ -121,16 +143,18 @@ 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(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
+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
!
! prognostic variables at t- deltat
!
-REAL, DIMENSION(:,:), INTENT(IN) :: PPABSM ! Pressure at time t-1
+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
@@ -138,16 +162,21 @@ 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(INOUT) :: PW_UP2 ! Vertical velocity^2
-REAL, DIMENSION(:), INTENT(IN) :: PTHL_UP,PRT_UP ! updraft properties
+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
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(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
+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
!
!
! 1.2 Declaration of local variables
@@ -155,64 +184,38 @@ REAL, DIMENSION(:), INTENT(INOUT) :: PBUO_INTEG! Integral Buoyancy
!
! Variables for 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
-
+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
! Variables for dry part
-
-REAL, DIMENSION(SIZE(PTHLM)) :: ZBUO_INTEG,& ! Temporary integral Buoyancy
- ZDZ_HALF,& ! half-DeltaZ between 2 flux points
- ZDZ_STOP,& ! Exact Height of the LCL
+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
ZTHV_MINUS_HALF,& ! Thv at flux point(kk)
ZTHV_PLUS_HALF,& ! Thv at flux point(kk+kkl)
- 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
-
+ ZDZ ! Delta Z used in computations
+INTEGER :: JI
!----------------------------------------------------------------------------------
@@ -220,269 +223,251 @@ REAL :: ZRVORD ! RV/RD
! ------------------
- ZRDORV = XRD / XRV !=0.622
ZRVORD = XRV / XRD !=1.607
- ZG_O_THVREF=XG/PTHVM
+ ZG_O_THVREF(:)=XG/PTHVM(:,KK)
+ ZCOEFFMF_CLOUD=XENTR_MF * XG / XCRAD_MF
- 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
- 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)
+ ZPRE(:)=PPRE_MINUS_HALF(:)
+ ZMIXTHL(:)=0.1
+ ZMIXRT(:)=0.1
+
+! 1.4 Estimation of PPART_DRY
+ WHERE(OTEST)
+ WHERE(OTESTLCL)
+ !No dry part when condensation level is reached
+ PPART_DRY(:)=0.
+ ZDZ_STOP(:)=0.
+ ZPRE(:)=PPRE_MINUS_HALF(:)
+ ELSEWHERE
+ !Temperature at flux level KK
+ ZT(:)=PTH_UP(:)*(PPRE_MINUS_HALF(:)/XP00) ** (XRD/XCPD)
+ !Saturating vapor pressure at flux level KK
+ ZFOESW(:) = MIN(EXP( XALPW - XBETAW/ZT(:) - XGAMW*LOG(ZT(:)) ), 0.99*PPRE_MINUS_HALF(:))
+ ZFOESI(:) = MIN(EXP( XALPI - XBETAI/ZT(:) - XGAMI*LOG(ZT(:)) ), 0.99*PPRE_MINUS_HALF(:))
+ !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(:)=(XBETAW/ZT(:)-XGAMW)*(1-ZFRAC_ICE(:))+(XBETAI/ZT(:)-XGAMI)*ZFRAC_ICE(:)
+ ZDRSATODP(:)=((XRD/XCPD)*ZDRSATODP(:)-1.)*PRSAT_UP(:)/ &
+ &(PPRE_MINUS_HALF(:)-(ZFOESW(:)*(1-ZFRAC_ICE(:)) + ZFOESI(:)*ZFRAC_ICE(:)))
+ !Use of d.Rsat / dP and pressure at flux level KK to find pressure (ZPRE)
+ !where Rsat is equal to PRT_UP
+ ZPRE(:)=PPRE_MINUS_HALF(:)+(PRT_UP(:)-PRSAT_UP(:))/ZDRSATODP(:)
+ !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(:)=MAX(0., MIN(1., (PPRE_MINUS_HALF(:)-ZPRE(:))/(PPRE_MINUS_HALF(:)-PPRE_PLUS_HALF(:))))
+ !Height above flux level KK of the cloudy part
+ ZDZ_STOP(:) = (PZZ(:,KK+KKL)-PZZ(:,KK))*PPART_DRY(:)
+ ENDWHERE
+ ENDWHERE
- ! Compute pressure at flux level KK+KKL
+! 1.5 Gradient and flux values of thetav
IF(KK/=KKB)THEN
- 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)
+ 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))
ELSE
- ZPRE_MINUS_HALF(:)=PPABSM(:,KK)
+ ZCOEFF_MINUS_HALF(:)=0.
+ ZTHV_MINUS_HALF(:) = PTHVM(:,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
+! --------------------------------------------------------------------
+
+ WHERE(OTEST .AND. PPART_DRY(:)>0.)
+ !Buoyancy computation in two parts to use change of gradient of theta v of environment
+ !Between flux level KK and min(mass level, bottom of cloudy part)
+ ZDZ(:)=MIN(ZDZ_STOP(:),(PZZ(:,KK+KKL)-PZZ(:,KK))*0.5)
+ PBUO_INTEG_DRY(:) = ZG_O_THVREF(:)*ZDZ(:)*&
+ (0.5 * ( - ZCOEFF_MINUS_HALF(:))*ZDZ(:) &
+ - ZTHV_MINUS_HALF(:) + PTHV_UP(:) )
+
+ !Between mass flux KK and bottom of cloudy part (if above mass flux)
+ ZDZ(:)=MAX(0., ZDZ_STOP(:)-(PZZ(:,KK+KKL)-PZZ(:,KK))*0.5)
+ PBUO_INTEG_DRY(:) = PBUO_INTEG_DRY(:) + ZG_O_THVREF(:)*ZDZ(:)*&
+ (0.5 * ( - ZCOEFF_PLUS_HALF(:))*ZDZ(:) &
+ - PTHVM(:,KK) + PTHV_UP(:) )
+
+ !Entr//Detr. computation
+ WHERE (PBUO_INTEG_DRY(:)>=0.)
+ PENTR(:) = 0.5/(XABUO-XBENTR*XENTR_DRY)*&
+ LOG(1.+ (2.*(XABUO-XBENTR*XENTR_DRY)/PW_UP2(:,KK))* &
+ PBUO_INTEG_DRY(:))
+ PDETR(:) = 0.
+ ELSEWHERE
+ PENTR(:) = 0.
+ PDETR(:) = 0.5/(XABUO)*&
+ LOG(1.+ (2.*(XABUO)/PW_UP2(:,KK))* &
+ (-PBUO_INTEG_DRY(:)))
+ ENDWHERE
+ PENTR(:) = XENTR_DRY*PENTR(:)/(PZZ(:,KK+KKL)-PZZ(:,KK))
+ PDETR(:) = XDETR_DRY*PDETR(:)/(PZZ(:,KK+KKL)-PZZ(:,KK))
+ !Minimum value of detrainment
+ ZWK(:)=PLUP(:)-0.5*(PZZ(:,KK)+PZZ(:,KK+KKL))
+ ZWK(:)=SIGN(MAX(1., ABS(ZWK(:))), ZWK(:)) ! ZWK must not be zero
+ PDETR(:) = MAX(PPART_DRY(:)*XDETR_LUP/ZWK(:), PDETR(:))
+ ELSEWHERE
+ !No dry part, consation reached (OTESTLCL)
+ PBUO_INTEG_DRY(:) = 0.
+ PENTR(:)=0.
+ PDETR(:)=0.
+ ENDWHERE
- ! Compute non cons. var. of mixture at the mass level
- CALL TH_R_FROM_THL_RT_1D(HFRAC_ICE,ZFRAC_ICE,&
- 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.
-! -----------------------------------------------------
+! 3 Wet part computation
+! -----------------------
+! 3.1 Integral buoyancy for cloudy part
-! 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_F1(:) = ZTHMIX(:)*(1.+ZRVORD*ZRVMIX(:))/(1.+PRT_UP(:))
-
- ! Compute theta_v of updraft at mass level KK
- ZRCMIX=PRC_UP
- ZRIMIX=PRI_UP
+ ! 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(:)
CALL TH_R_FROM_THL_RT_1D(HFRAC_ICE,ZFRAC_ICE,&
- PPABSM(:,KK),PTHL_UP,PRT_UP,&
+ PPRE_PLUS_HALF,PTHL_UP,PRT_UP,&
ZTHMIX,ZRVMIX,ZRCMIX,ZRIMIX,&
ZRSATW, ZRSATI)
- ZTHV_UP(:) = ZTHMIX(:)*(1.+ZRVORD*ZRVMIX(:))/(1.+PRT_UP(:))
+ ZTHV_UP_F2(:) = ZTHMIX(:)*(1.+ZRVORD*ZRVMIX(:))/(1.+PRT_UP(:))
+
+ ! Integral buoyancy for cloudy part
+ WHERE(OTEST .AND. PPART_DRY(:)<1.)
+ !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(:)=(ZTHV_UP_F2(:)-PTHV_UP(:))/((PZZ(:,KK+KKL)-PZZ(:,KK))*(1-PPART_DRY(:)))
+
+ !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(:)=MAX(0., 0.5*(PZZ(:,KK+KKL)-PZZ(:,KK))-ZDZ_STOP(:))
+ PBUO_INTEG_CLD(:) = ZG_O_THVREF(:)*ZDZ(:)*&
+ (0.5*( ZCOTHVU(:) - ZCOEFF_MINUS_HALF(:))*ZDZ(:) &
+ - (PTHVM(:,KK)-ZDZ(:)*ZCOEFF_MINUS_HALF(:)) + PTHV_UP(:) )
+
+ !Between max(mass level, bottom of cloudy part) and flux level KK+KKL
+ ZDZ(:)=(PZZ(:,KK+KKL)-PZZ(:,KK))-MAX(ZDZ_STOP(:),0.5*(PZZ(:,KK+KKL)-PZZ(:,KK)))
+ PBUO_INTEG_CLD(:) = PBUO_INTEG_CLD(:)+ZG_O_THVREF(:)*ZDZ(:)*&
+ (0.5*( ZCOTHVU(:) - ZCOEFF_PLUS_HALF(:))*ZDZ(:)&
+ - (PTHVM(:,KK)+(0.5*((PZZ(:,KK+KKL)-PZZ(:,KK)))-ZDZ(:))*ZCOEFF_PLUS_HALF(:)) +&
+ PTHV_UP(:) )
+
+ ELSEWHERE
+ !No cloudy part
+ PBUO_INTEG_CLD(:)=0.
+ ENDWHERE
- ! Compute theta_v of updraft at flux level KK+KKL
- ZRCMIX_F2=PRC_UP
- ZRIMIX_F2=PRI_UP
+! 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)
+
+ WHERE(OTEST .AND. PPART_DRY(:)>0.5)
+ ZDZ(:)=ZDZ_STOP(:)-0.5*(PZZ(:,KK+KKL)-PZZ(:,KK))
+ ZTHV(:)= PTHVM(:,KK)+ZCOEFF_PLUS_HALF(:)*ZDZ(:)
+ ZMIXTHL(:) = ZKIC_INIT * &
+ (PTHLM(:,KK)+ZDZ(:)*(PTHLM(:,KK+KKL)-PTHLM(:,KK))/PDZZ(:,KK+KKL)) + &
+ (1. - ZKIC_INIT)*PTHL_UP(:)
+ ZMIXRT(:) = ZKIC_INIT * &
+ (PRTM(:,KK)+ZDZ(:)*(PRTM(:,KK+KKL)-PRTM(:,KK))/PDZZ(:,KK+KKL)) + &
+ (1. - ZKIC_INIT)*PRT_UP(:)
+ ELSEWHERE(OTEST)
+ ZDZ(:)=0.5*(PZZ(:,KK+KKL)-PZZ(:,KK))-ZDZ_STOP(:)
+ ZTHV(:)= PTHVM(:,KK)-ZCOEFF_MINUS_HALF(:)*ZDZ(:)
+ ZMIXTHL(:) = ZKIC_INIT * &
+ (PTHLM(:,KK)-ZDZ(:)*(PTHLM(:,KK)-PTHLM(:,JI))/PDZZ(:,KK)) + &
+ (1. - ZKIC_INIT)*PTHL_UP(:)
+ ZMIXRT(:) = ZKIC_INIT * &
+ (PRTM(:,KK)-ZDZ(:)*(PRTM(:,KK)-PRTM(:,JI))/PDZZ(:,KK)) + &
+ (1. - ZKIC_INIT)*PRT_UP(:)
+ ENDWHERE
CALL TH_R_FROM_THL_RT_1D(HFRAC_ICE,ZFRAC_ICE,&
- ZPRE_PLUS_HALF,PTHL_UP,PRT_UP,&
- ZTHMIX_F2,ZRVMIX_F2,ZRCMIX_F2,ZRIMIX_F2,&
+ ZPRE,ZMIXTHL,ZMIXRT,&
+ ZTHMIX,ZRVMIX,PRC_MIX,PRI_MIX,&
ZRSATW, ZRSATI)
- ZTHV_UP_F2(:) = ZTHMIX_F2(:)*(1.+ZRVORD*ZRVMIX_F2(:))/(1.+PRT_UP(:))
+ ZTHVMIX(:) = ZTHMIX(:)*(1.+ZRVORD*ZRVMIX(:))/(1.+ZMIXRT(:))
- ! Computation of RC and RI on mass point KK+KKL
- ZRCMIX_M2=PRC_UP
- ZRIMIX_M2=PRI_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(:)
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,&
+ PPRE_PLUS_HALF,ZMIXTHL,ZMIXRT,&
+ ZTHMIX,ZRVMIX,PRC_MIX,PRI_MIX,&
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
+ ZTHVMIX_F2(:) = ZTHMIX(:)*(1.+ZRVORD*ZRVMIX(:))/(1.+ZMIXRT(:))
+
+ !Computation of mean ZKIC over the cloudy part
+ WHERE (OTEST)
+ ! Compute ZKIC at the bottom of cloudy part
+ ! Thetav_up at bottom is equal to Thetav_up at flux level KK
+ WHERE (ABS(PTHV_UP(:)-ZTHVMIX(:))<1.E-10)
+ ZKIC(:)=1.
+ ELSEWHERE
+ ZKIC(:) = MAX(0.,PTHV_UP(:)-ZTHV(:))*ZKIC_INIT / &
+ (PTHV_UP(:)-ZTHVMIX(:))
+ ENDWHERE
+ ! Compute ZKIC_F2 at flux level KK+KKL
+ WHERE (ABS(ZTHV_UP_F2(:)-ZTHVMIX_F2(:))<1.E-10)
+ ZKIC_F2(:)=1.
+ ELSEWHERE
+ ZKIC_F2(:) = MAX(0.,ZTHV_UP_F2(:)-ZTHV_PLUS_HALF(:))*ZKIC_INIT / &
+ (ZTHV_UP_F2(:)-ZTHVMIX_F2(:))
+ ENDWHERE
+ !Mean ZKIC over the cloudy part
+ ZKIC(:)=MAX(MIN(0.5*(ZKIC(:)+ZKIC_F2(:)),1.),0.)
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))
+! 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
- ZW2_HALF = PW_UP2(:,KK) + 2*(XABUO)*(ZBUO_INTEG)
- ENDWHERE
- ENDWHERE
+ !Constant PDF
+ !For this PDF, eq. (5) is delta Me=0.5*delta Mt
+ WHERE(OTEST)
+ ZEPSI(:) = ZKIC(:)**2. !integration multiplied by 2
+ ZDELTA(:) = (1.-ZKIC(:))**2. !idem
+ 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
+ !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
+ WHERE (OTEST)
+ ZEPSI_CLOUD=MIN(ZDELTA,ZEPSI)
+ PENTR_CLD(:) = (1.-PPART_DRY(:))*ZCOEFFMF_CLOUD*PRHODREF(:)*ZEPSI_CLOUD(:)
+ PDETR_CLD(:) = (1.-PPART_DRY(:))*ZCOEFFMF_CLOUD*PRHODREF(:)*ZDELTA(:)
+ PENTR(:) = PENTR(:)+PENTR_CLD(:)
+ PDETR(:) = PDETR(:)+PDETR_CLD(:)
+ ELSEWHERE
+ PENTR_CLD(:) = 0.
+ PDETR_CLD(:) = 0.
+ ENDWHERE
END SUBROUTINE COMPUTE_ENTR_DETR
diff --git a/src/MNH/compute_updraft.f90 b/src/MNH/compute_updraft.f90
index 1d0e2a821..714f53422 100644
--- a/src/MNH/compute_updraft.f90
+++ b/src/MNH/compute_updraft.f90
@@ -10,6 +10,7 @@ INTERFACE
!
! #################################################################
SUBROUTINE COMPUTE_UPDRAFT(KKA,KKB,KKE,KKU,KKL, HFRAC_ICE, &
+ HMF_UPDRAFT, &
OENTR_DETR,OMIXUV, &
ONOMIXLG,KSV_LGBEG,KSV_LGEND, &
PZZ,PDZZ, &
@@ -35,6 +36,7 @@ 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
@@ -83,6 +85,7 @@ 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, &
@@ -129,6 +132,8 @@ 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
@@ -142,6 +147,7 @@ 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
@@ -156,8 +162,9 @@ 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
-LOGICAL, INTENT(IN) :: OENTR_DETR! flag to recompute entrainment, detrainment and mass flux
-LOGICAL, INTENT(IN) :: OMIXUV ! True if mixing of momentum
+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
INTEGER, INTENT(IN) :: KSV_LGBEG ! first index of lag. tracer
INTEGER, INTENT(IN) :: KSV_LGEND ! last index of lag. tracer
@@ -207,7 +214,9 @@ 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
+ ZW_UP2, & ! w**2 of the updraft
+ ZBUO_INTEG_DRY, ZBUO_INTEG_CLD,&! Integrated Buoyancy
+ ZENTR_CLD,ZDETR_CLD ! wet entrainment and detrainment
REAL, DIMENSION(SIZE(PSVM,1),SIZE(PTHM,2),SIZE(PSVM,3)) :: &
ZSVM_F ! scalar variables
@@ -225,7 +234,7 @@ REAL :: ZRDORV ! RD/RV
REAL :: ZRVORD ! RV/RD
-REAL, DIMENSION(SIZE(PTHM,1)) :: ZMIX1,ZMIX2,ZMIX3
+REAL, DIMENSION(SIZE(PTHM,1)) :: ZMIX1,ZMIX2,ZMIX3_CLD,ZMIX2_CLD
REAL, DIMENSION(SIZE(PTHM,1)) :: ZLUP ! Upward Mixing length from the ground
@@ -241,12 +250,16 @@ LOGICAL, DIMENSION(SIZE(PTHM,1),SIZE(PTHM,2)) :: GWORK2
INTEGER :: ITEST
-REAL, DIMENSION(SIZE(PTHM,1)) :: ZRC_UP, ZRI_UP, ZRV_UP, ZRSATW, ZRSATI
+REAL, DIMENSION(SIZE(PTHM,1)) :: ZRC_UP, ZRI_UP, ZRV_UP,&
+ ZRSATW, ZRSATI,&
+ ZPART_DRY
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
@@ -263,7 +276,6 @@ 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)
@@ -361,13 +373,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,ZG_O_THVREF,ZTHVM_F,KKB,GLMIX,ZLUP)
+ CALL COMPUTE_BL89_ML(KKA,KKB,KKE,KKU,KKL,PDZZ,ZTKEM_F(:,KKB),ZG_O_THVREF(:,KKB),ZTHVM,KKB,GLMIX,.TRUE.,ZLUP)
ZLUP(:)=MAX(ZLUP(:),1.E-10)
! Compute Buoyancy flux at the ground
@@ -375,8 +387,14 @@ 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 * ZRHO_F(:,KKB) * ((ZG_O_THVREF(:,KKB))*ZWTHVSURF*ZLUP)**(1./3.)
+ PEMF(:,KKB) = XCMF * ZSURF(:) * 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.
@@ -418,7 +436,6 @@ 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
@@ -426,14 +443,20 @@ 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),&
- PPABSM(:,:),PZZ(:,:),PDZZ(:,:),ZTHVM(:,:), &
- PTHLM(:,JK),PRTM(:,JK),ZW_UP2(:,:), &
+ PRHODREF(:,JK),ZPRES_F(:,JK),ZPRES_F(:,JK+KKL),&
+ PZZ(:,:),PDZZ(:,:),ZTHVM(:,:), &
+ PTHLM(:,:),PRTM(:,:),ZW_UP2(:,:),ZTH_UP(:,JK), &
PTHL_UP(:,JK),PRT_UP(:,JK),ZLUP(:), &
- PRC_UP(:,JK),PRI_UP(:,JK),ZRC_MIX(:,JK),ZRI_MIX(:,JK), &
- PENTR(:,JK),PDETR(:,JK),PBUO_INTEG(:,JK) )
+ 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)
IF (JK==KKB) THEN
PDETR(:,JK)=0.
+ ZDETR_CLD(:,JK)=0.
ENDIF
! Computation of updraft characteristics at level JK+KKL
@@ -459,7 +482,8 @@ 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(:) = (PZZ(:,JK+KKL)-PZZ(:,JK))*PDETR(:,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)
PTHL_UP(:,JK+KKL)=(PTHL_UP(:,JK)*(1.-0.5*ZMIX2(:)) + PTHLM(:,JK)*ZMIX2(:)) &
/(1.+0.5*ZMIX2(:))
@@ -522,20 +546,15 @@ 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 (.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
+ 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(:)))
ENDWHERE
@@ -602,10 +621,14 @@ IF(OENTR_DETR) THEN
PDEPTH(JI) = MAX(0., PZZ(JI,KKCTL(JI)) - PZZ(JI,KKLCL(JI)) )
END DO
-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))
+ 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
END SUBROUTINE COMPUTE_UPDRAFT
diff --git a/src/MNH/compute_updraft_hrio.f90 b/src/MNH/compute_updraft_hrio.f90
index 022f1b3fd..c2e0bae75 100644
--- a/src/MNH/compute_updraft_hrio.f90
+++ b/src/MNH/compute_updraft_hrio.f90
@@ -134,6 +134,7 @@ 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
@@ -297,16 +298,6 @@ ZRMAX=1.E-3
ZEPS=1.E-15
XFRAC_LIM=0.5
!
-if(LDUMMY1)THEN
-print*,"XPRES_UV=",XPRES_UV
-print*,"XALP_PERT=",XALP_PERT
-print*,"XCMF=",XCMF
-print*,"XFRAC_UP_MAX=",XFRAC_UP_MAX
-print*,"XA1=",XA1
-print*,"XB=",XB
-print*,"XC=",XC
-print*,"XBETA1=",XBETA1
-END IF
!------------------------------------------------------------------------
! INITIALISATION
@@ -400,9 +391,6 @@ PSV_DO(:,:,:)=0.
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))
-!------------------------
-print*,OENTR_DETR
-!------------------------
IF (OENTR_DETR) THEN
ZTHM_F (:,:) = MZM_MF(KKA,KKU,KKL,PTHM (:,:))
ZPRES_F(:,:) = MZM_MF(KKA,KKU,KKL,PPABSM(:,:))
@@ -457,7 +445,8 @@ IF (OENTR_DETR) THEN
GLMIX=.TRUE.
ZTKEM_F(:,KKB)=0.
- CALL COMPUTE_BL89_ML(KKA,KKB,KKE,KKU,KKL,PDZZ,ZTKEM_F,ZG_O_THVREF,ZTHVM_F,KKB,GLMIX,ZLUP)
+ CALL COMPUTE_BL89_ML(KKA,KKB,KKE,KKU,KKL,PDZZ,ZTKEM_F(:,KKB),ZG_O_THVREF(:,KKB), &
+ ZTHVM_F,KKB,GLMIX,.TRUE.,ZLUP)
ZLUP(:)=MAX(ZLUP(:),1.E-10)
! Compute Buoyancy flux at the ground
@@ -481,13 +470,6 @@ IF (OENTR_DETR) THEN
ELSE
GTEST(:)=PEMF(:,KKB+KKL)>0.
END IF
- !=================
- IF(LDUMMY1)THEN
- print*,"PEMF(:,",KKB,")=", minval(PEMF(:,KKB)),maxval(PEMF(:,KKB))
- print*,"PFRAC_UP(:,",KKB,")=", minval(PFRAC_UP(:,KKB)),maxval(PFRAC_UP(:,KKB))
- print*,"ZW_UP2(:,",KKB,")=", minval(ZW_UP2(:,KKB)),maxval(ZW_UP2(:,KKB))
- END IF
- !=================
!--------------------------------------------------------------------------
@@ -507,9 +489,6 @@ GTESTETL(:)=.FALSE.
! print*,"NLES_DTCOUNT doit être 1."
!STOP
!ENDIF
-IF(LDUMMY1) THEN
-print*,"GTEST=",GTEST
-END IF
WHERE (GTEST .AND. ZLUP(:)>0. )
! hauteur des thermiques du pas de temps precedent
@@ -521,9 +500,6 @@ END IF
ZRESOL_NORM(:)=0.5 ! ARBITRAIRE
ENDWHERE
-IF(LDUMMY1) THEN
-print*,"ZRESOL_NORM=",minval(ZRESOL_NORM),maxval(ZRESOL_NORM)
-END IF
! Loop on vertical level
DO JK=KKB,KKE-KKL,KKL
@@ -531,12 +507,6 @@ DO JK=KKB,KKE-KKL,KKL
ITEST=COUNT(GTEST)
!IF (ITEST==0) print*,"cycle"
IF (ITEST==0) CYCLE
-IF(LDUMMY1) THEN
- print*,"JK=",JK
- print*,"alt=",minval(PZZ(:,JK)),maxval(PZZ(:,JK))
- print*,"ITEST=",ITEST
- IF (ITEST==0) STOP
-END IF
! Computation of entrainment and detrainment with KF90
! parameterization in clouds and LR01 in subcloud layer
@@ -568,11 +538,6 @@ END IF
ENDWHERE ! fin temporaire de GTEST
! test sur le calcul de la flottabilité de 2 manières
- IF(LDUMMY1) THEN
- print*,"PTHV_UP(:,",JK,")=",minval(PTHV_UP(:,JK)),maxval(PTHV_UP(:,JK))
- print*,"ZTHVM_F(:,",JK,")=",minval(ZTHVM_F(:,JK)),maxval(ZTHVM_F(:,JK))
- print*,"PTHVREF(:,",JK,")=",minval(PTHVREF(:,JK)),maxval(PTHVREF(:,JK))
- END IF
! le calcul de la flottabilité est super important
! c'est ce qui decide de Wup
!IF(NBUOY==1) THEN
@@ -626,9 +591,6 @@ END IF
! print*,"ZBUO_0(:,",JK,")=",minval(ZBUO(:,JK)),maxval(ZBUO(:,JK))
!END IF
!ENDIF
- IF(LDUMMY1) THEN
- print*,"ZBUO(:,",JK,")=",minval(ZBUO(:,JK)),maxval(ZBUO(:,JK))
- END IF
!WHERE (GTEST)
! anomalie de flottabilité * DZ
@@ -655,17 +617,9 @@ END IF
!END IF
!
- IF(LDUMMY1) THEN
- print*,"ZA1=",minval(ZA1),maxval(ZA1)
- END IF
!================================================================================
! EQUATION DE LA DYNAMIQUE
!================================================================================
- IF(LDUMMY1) THEN
- print*,"PBUO_INTEG(:,",JK,")=",minval(PBUO_INTEG(:,JK)),maxval(PBUO_INTEG(:,JK))
- print*,"PWM(:,",JK,")=",minval(PWM(:,JK)),maxval(PWM(:,JK))
- print*,"ZWM_M(:,",JK,")=",minval(ZWM_M(:,JK)),maxval(ZWM_M(:,JK))
- END IF
! on calcule le vent vertical du thermique dans la zone grise
! ZA1 depend de la méthode utilisée
ALLOCATE(ZWORK(SIZE(ZW_UP,1)))
@@ -694,10 +648,6 @@ END IF
ENDWHERE
DEALLOCATE(ZWORK)
!
- IF(LDUMMY1) THEN
- print*,"ZW_UP(:,",JK+KKL,")=",minval(ZW_UP(:,JK+KKL)),maxval(ZW_UP(:,JK+KKL))
- print*,"ZW_UP2(:,",JK+KKL,")=",minval(ZW_UP2(:,JK+KKL)),maxval(ZW_UP2(:,JK+KKL))
- END IF
!================================================================================
! CALCUL DE L ENTRAINEMENT ET DU DETRAINEMENT
!================================================================================
@@ -716,10 +666,6 @@ END IF
PDETR(:,JK) = ZDETR_RT(:)+ZDETR_BUO(:)
ENDWHERE
!
- IF(LDUMMY1) THEN
- print*,"PENTR(:,",JK,")=",minval(PENTR(:,JK)),maxval(PENTR(:,JK))
- print*,"PDETR(:,",JK,")=",minval(PDETR(:,JK)),maxval(PDETR(:,JK))
- END IF
!================================================================================
! VARIABLES THERMODYNAMIQUES
!================================================================================
@@ -795,10 +741,6 @@ END IF
WHERE(GTEST)
PTHV_UP(:,JK+KKL) = ZTH_UP(:,JK+KKL)*((1+ZRVORD*PRV_UP(:,JK+KKL))/(1+PRT_UP(:,JK+KKL)))
ENDWHERE
- !IF(LDUMMY1) THEN
- !print*,"PTHV_UP(:,",JK+KKL,")=",minval(PTHV_UP(:,JK+KKL)),maxval(PTHV_UP(:,JK+KKL))
- !print*,"ZTHVM_F(:,",JK+KKL,")=",minval(ZTHVM_F(:,JK+KKL)),maxval(ZTHVM_F(:,JK+KKL))
- !END IF
!================================================================================
! CALCUL DU FLUX DE MASSE
!================================================================================
@@ -808,9 +750,6 @@ END IF
ZMIX1(:)=ZZDZ(:,JK)*(PENTR(:,JK)-PDETR(:,JK))
PEMF(:,JK+KKL)=PEMF(:,JK)*EXP(ZMIX1(:))
ENDWHERE
- IF(LDUMMY1) THEN
- print*,"PEMF(:,",JK+KKL,")=",minval(PEMF(:,JK+KKL)),maxval(PEMF(:,JK+KKL))
- END IF
!================================================================================
! FRACTION DE THERMIQUE
!================================================================================
@@ -831,9 +770,6 @@ ENDWHERE
ENDWHERE
- IF(LDUMMY1)THEN
- print*,"PFRAC_UP(:,",JK+KKL,")=",minval(PFRAC_UP(:,JK+KKL)),maxval(PFRAC_UP(:,JK+KKL))
- END IF
! calcul des termes environmentaux au point de flux
WHERE(GTEST)
!WHERE(PFRAC_UP(:,JK+KKL)>0 .AND. PFRAC_UP(:,JK+KKL)< XFRAC_LIM)
@@ -870,15 +806,8 @@ ENDWHERE
KKETL(:) = JK+KKL
GTESTETL(:)=.TRUE.
ENDWHERE
-IF(LDUMMY1)THEN
-print*,"GTESTETL=",GTESTETL
-END IF
! Test is we have reached the top of the updraft
-IF(LDUMMY1)THEN
-print*,"ZW_UP2(:,",JK+KKL,")=",minval(ZW_UP2(:,JK+KKL)),maxval(ZW_UP2(:,JK+KKL))
-print*,"ZEPS=",ZEPS
-END IF
! WHERE (GTEST.AND.((ZW_UP2(:,JK+KKL)<=ZEPS).OR.(PEMF(:,JK+KKL)<=ZEPS) .OR. PFRAC_UP(:,JK+KKL)<= XALPHA_SEUIL))
WHERE ( GTEST .AND. ( (ZW_UP2(:,JK+KKL)<=10E-5) .OR. (PEMF(:,JK+KKL)<=10E-5)) )
@@ -918,6 +847,5 @@ ENDDO ! boucle JK
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))
-print*,"je sors de compute_updraft"
END SUBROUTINE COMPUTE_UPDRAFT_HRIO
diff --git a/src/MNH/compute_updraft_raha.f90 b/src/MNH/compute_updraft_raha.f90
new file mode 100644
index 000000000..4ddeb8938
--- /dev/null
+++ b/src/MNH/compute_updraft_raha.f90
@@ -0,0 +1,667 @@
+!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 8b3f6a5e2..3a9614ec7 100644
--- a/src/MNH/compute_updraft_rhcj10.f90
+++ b/src/MNH/compute_updraft_rhcj10.f90
@@ -124,6 +124,7 @@ 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
@@ -390,8 +391,8 @@ ENDDO
GLMIX=.TRUE.
ZTKEM_F(:,KKB)=0.
-
-CALL COMPUTE_BL89_ML(KKA,KKB,KKE,KKU,KKL,PDZZ,ZTKEM_F,ZG_O_THVREF,ZTHVM_F,KKB,GLMIX,ZLUP)
+CALL COMPUTE_BL89_ML(KKA,KKB,KKE,KKU,KKL,PDZZ,ZTKEM_F(:,KKB),ZG_O_THVREF(:,KKB), &
+ ZTHVM_F,KKB,GLMIX,.TRUE.,ZLUP)
ZLUP(:)=MAX(ZLUP(:),1.E-10)
! Compute Buoyancy flux at the ground
diff --git a/src/MNH/mf_turb_expl.f90 b/src/MNH/mf_turb_expl.f90
new file mode 100644
index 000000000..cec383e9e
--- /dev/null
+++ b/src/MNH/mf_turb_expl.f90
@@ -0,0 +1,227 @@
+!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.
+! ######################
+ MODULE MODI_MF_TURB_EXPL
+! ######################
+!
+INTERFACE
+!
+! #################################################################
+ SUBROUTINE MF_TURB_EXPL(KKA,KKB,KKE,KKU,KKL,OMIXUV, &
+ PRHODJ, &
+ PTHLM,PTHVM,PRTM,PUM,PVM, &
+ PTHLDT,PRTDT,PUDT,PVDT, &
+ PEMF,PTHL_UP,PTHV_UP,PRT_UP,PU_UP,PV_UP, &
+ PFLXZTHLMF,PFLXZTHVMF,PFLXZRMF,PFLXZUMF,PFLXZVMF)
+! #################################################################
+!
+!* 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
+LOGICAL, INTENT(IN) :: OMIXUV ! True if mixing of momentum
+
+REAL, DIMENSION(:,:), INTENT(IN) :: PRHODJ ! dry density * Grid size
+
+! Conservative var. at t-dt
+REAL, DIMENSION(:,:), INTENT(IN) :: PTHLM ! conservative pot. temp.
+REAL, DIMENSION(:,:), INTENT(IN) :: PRTM ! water var. where
+
+! Virtual potential temperature at t-dt
+REAL, DIMENSION(:,:), INTENT(IN) :: PTHVM
+! Momentum at t-dt
+REAL, DIMENSION(:,:), INTENT(IN) :: PUM
+REAL, DIMENSION(:,:), INTENT(IN) :: PVM
+!
+! Tendencies of conservative variables
+REAL, DIMENSION(:,:), INTENT(OUT) :: PTHLDT
+
+REAL, DIMENSION(:,:), INTENT(OUT) :: PRTDT
+
+! Tendencies of momentum
+REAL, DIMENSION(:,:), INTENT(OUT) :: PUDT
+REAL, DIMENSION(:,:), INTENT(OUT) :: PVDT
+
+! Updraft characteritics
+REAL, DIMENSION(:,:), INTENT(IN) :: PEMF,PTHL_UP,PTHV_UP,PRT_UP,PU_UP,PV_UP
+
+! Fluxes
+REAL, DIMENSION(:,:), INTENT(OUT) :: PFLXZTHLMF,PFLXZTHVMF,PFLXZRMF,PFLXZUMF,PFLXZVMF
+
+END SUBROUTINE MF_TURB_EXPL
+
+END INTERFACE
+!
+END MODULE MODI_MF_TURB_EXPL
+!
+
+! ######spl
+ SUBROUTINE MF_TURB_EXPL(KKA,KKB,KKE,KKU,KKL,OMIXUV, &
+ PRHODJ, &
+ PTHLM,PTHVM,PRTM,PUM,PVM, &
+ PTHLDT,PRTDT,PUDT,PVDT, &
+ PEMF,PTHL_UP,PTHV_UP,PRT_UP,PU_UP,PV_UP, &
+ PFLXZTHLMF,PFLXZTHVMF,PFLXZRMF,PFLXZUMF,PFLXZVMF)
+
+! #################################################################
+!
+!
+!!**** *MF_TURB_EXPL* - computes the MF_turbulent source terms for the prognostic
+!! variables (when PIMPL=0)
+!!
+!! PURPOSE
+!! -------
+!!**** The purpose of this routine is to compute the source terms in
+!! the evolution equations due to the MF turbulent mixing.
+!! The source term is computed as the divergence of the turbulent fluxes.
+!
+!!** METHOD
+!! ------
+!!
+!!
+!! EXTERNAL
+!! --------
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!!
+!!
+!! REFERENCE
+!! ---------
+!!
+!! AUTHOR
+!! ------
+!!
+!!
+!! MODIFICATIONS
+!! -------------
+!!
+!! --------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+
+USE MODD_PARAM_MFSHALL_n
+USE MODI_SHUMAN_MF
+
+IMPLICIT NONE
+
+
+!* 0.1 declarations 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
+LOGICAL, INTENT(IN) :: OMIXUV ! True if mixing of momentum
+
+REAL, DIMENSION(:,:), INTENT(IN) :: PRHODJ ! dry density * Grid size
+
+! Conservative var. at t-dt
+REAL, DIMENSION(:,:), INTENT(IN) :: PTHLM ! conservative pot. temp.
+REAL, DIMENSION(:,:), INTENT(IN) :: PRTM ! water var. where
+
+! Virtual potential temperature at t-dt
+REAL, DIMENSION(:,:), INTENT(IN) :: PTHVM
+! Momentum at t-dt
+REAL, DIMENSION(:,:), INTENT(IN) :: PUM
+REAL, DIMENSION(:,:), INTENT(IN) :: PVM
+!
+! Tendencies of conservative variables
+REAL, DIMENSION(:,:), INTENT(OUT) :: PTHLDT
+
+REAL, DIMENSION(:,:), INTENT(OUT) :: PRTDT
+
+! Tendencies of momentum
+REAL, DIMENSION(:,:), INTENT(OUT) :: PUDT
+REAL, DIMENSION(:,:), INTENT(OUT) :: PVDT
+
+! Updraft characteritics
+REAL, DIMENSION(:,:), INTENT(IN) :: PEMF,PTHL_UP,PTHV_UP,PRT_UP,PU_UP,PV_UP
+
+! Fluxes
+REAL, DIMENSION(:,:), INTENT(OUT) :: PFLXZTHLMF,PFLXZTHVMF,PFLXZRMF,PFLXZUMF,PFLXZVMF
+
+REAL, DIMENSION(SIZE(PFLXZTHLMF,1),SIZE(PFLXZTHLMF,2)) :: ZFLXZTHSMF,ZTHS_UP,ZTHSM ! Theta S flux
+REAL, DIMENSION(SIZE(PFLXZTHLMF,1),SIZE(PFLXZTHLMF,2)) :: ZQT_UP,ZQTM,ZTHSDT,ZQTDT
+REAL, DIMENSION(SIZE(PFLXZTHLMF,1),SIZE(PFLXZTHLMF,2)) :: ZTHLM_F,ZRTM_F
+
+INTEGER :: JK ! loop counter
+
+!----------------------------------------------------------------------------
+!
+!* 1.PRELIMINARIES
+! -------------
+
+PFLXZRMF = 0.
+PFLXZTHVMF = 0.
+PFLXZTHLMF = 0.
+PFLXZUMF = 0.
+PFLXZVMF = 0.
+PTHLDT = 0.
+PRTDT = 0.
+PUDT = 0.
+PVDT = 0.
+
+!
+!----------------------------------------------------------------------------
+!
+!* 2. COMPUTE THE MEAN FLUX OF CONSERVATIVE VARIABLES at time t-dt
+! (equation (3) of Soares et al)
+! + THE MEAN FLUX OF THETA_V (buoyancy flux)
+! -----------------------------------------------
+! ( Resulting fluxes are in flux level (w-point) as PEMF and PTHL_UP )
+
+ZRTM_F (:,:) = MZM_MF(KKA,KKU,KKL,PRTM (:,:))
+ZTHLM_F(:,:) = MZM_MF(KKA,KKU,KKL,PTHLM(:,:))
+ZQTM (:,:) = ZRTM_F (:,:)/(1.+ZRTM_F (:,:))
+ZQT_UP (:,:) = PRT_UP (:,:)/(1.+PRT_UP (:,:))
+ZTHS_UP(:,:) = PTHL_UP(:,:)*(1.+XLAMBDA*ZQT_UP(:,:))
+ZTHSM (:,:) = ZTHLM_F(:,:)*(1.+XLAMBDA*ZQTM(:,:))
+
+PFLXZTHLMF(:,:) = PEMF(:,:)*(PTHL_UP(:,:)-MZM_MF(KKA,KKU,KKL,PTHLM(:,:))) ! ThetaL
+PFLXZRMF(:,:) = PEMF(:,:)*(PRT_UP (:,:)-MZM_MF(KKA,KKU,KKL,PRTM (:,:))) ! Rt
+PFLXZTHVMF(:,:) = PEMF(:,:)*(PTHV_UP(:,:)-MZM_MF(KKA,KKU,KKL,PTHVM(:,:))) ! ThetaV
+
+ZFLXZTHSMF(:,:) = PEMF(:,:)*(ZTHS_UP(:,:)-ZTHSM(:,:)) ! Theta S flux
+
+IF (OMIXUV) THEN
+ PFLXZUMF(:,:) = PEMF(:,:)*(PU_UP(:,:)-MZM_MF(KKA,KKU,KKL,PUM(:,:))) ! U
+ PFLXZVMF(:,:) = PEMF(:,:)*(PV_UP(:,:)-MZM_MF(KKA,KKU,KKL,PVM(:,:))) ! V
+ELSE
+ PFLXZUMF(:,:) = 0.
+ PFLXZVMF(:,:) = 0.
+ENDIF
+
+
+!----------------------------------------------------------------------------
+!
+!* 3. COMPUTE TENDENCIES OF CONSERVATIVE VARIABLES (or treated as such...)
+! (explicit formulation)
+! --------------------------------------------
+
+DO JK=KKB,KKE-KKL,KKL
+! PTHLDT(:,JK) = (PFLXZTHLMF(:,JK ) - PFLXZTHLMF(:,JK+KKL)) / PRHODJ(:,JK)
+ PRTDT (:,JK) = (PFLXZRMF (:,JK ) - PFLXZRMF (:,JK+KKL)) / PRHODJ(:,JK)
+ ZQTDT (:,JK) = PRTDT (:,JK)/(1.+ ZRTM_F (:,JK)*ZRTM_F (:,JK))
+ ZTHSDT(:,JK) = (ZFLXZTHSMF(:,JK ) - ZFLXZTHSMF(:,JK+KKL)) / PRHODJ(:,JK)
+ PTHLDT(:,JK) = ZTHSDT(:,JK)/(1.+XLAMBDA*ZQTM(:,JK)) - ZTHLM_F(:,JK)*XLAMBDA*ZQTDT(:,JK)
+END DO
+
+IF (OMIXUV) THEN
+ DO JK=KKB,KKE-KKL,KKL
+ PUDT(:,JK) = (PFLXZUMF(:,JK ) - PFLXZUMF(:,JK+KKL)) / PRHODJ(:,JK)
+ PVDT(:,JK) = (PFLXZVMF(:,JK ) - PFLXZVMF(:,JK+KKL)) / PRHODJ(:,JK)
+ END DO
+ENDIF
+
+
+END SUBROUTINE MF_TURB_EXPL
diff --git a/src/MNH/modd_param_mfshalln.f90 b/src/MNH/modd_param_mfshalln.f90
index 59e90ae3b..d501f7c6f 100644
--- a/src/MNH/modd_param_mfshalln.f90
+++ b/src/MNH/modd_param_mfshalln.f90
@@ -29,7 +29,8 @@
!!
!! MODIFICATIONS
!! -------------
-!! Original 01/02/07
+!! Original 01/02/07
+!! 10/16 R.Honnert Update with AROME
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -77,6 +78,14 @@ REAL :: XA1
REAL :: XB
REAL :: XC
REAL :: XBETA1
+! Parameters for closure assumption of Hourdin et al 2002
+
+REAL :: XR ! Aspect ratio of updraft
+
+! Thermodynamic parameter
+
+REAL :: XLAMBDA ! Lambda to compute ThetaS1 from ThetaL
+
END TYPE PARAM_MFSHALL_t
TYPE(PARAM_MFSHALL_t), DIMENSION(JPMODELMAX), TARGET, SAVE :: PARAM_MFSHALL_MODEL
@@ -108,6 +117,8 @@ REAL, POINTER :: XA1=>NULL()
REAL, POINTER :: XB=>NULL()
REAL, POINTER :: XC=>NULL()
REAL, POINTER :: XBETA1=>NULL()
+REAL, POINTER :: XR=>NULL()
+REAL, POINTER :: XLAMBDA=>NULL()
CONTAINS
@@ -144,6 +155,8 @@ XA1=>PARAM_MFSHALL_MODEL(KTO)%XA1
XB=>PARAM_MFSHALL_MODEL(KTO)%XB
XC=>PARAM_MFSHALL_MODEL(KTO)%XC
XBETA1=>PARAM_MFSHALL_MODEL(KTO)%XBETA1
+XR=>PARAM_MFSHALL_MODEL(KTO)%XR
+XLAMBDA=>PARAM_MFSHALL_MODEL(KTO)%XLAMBDA
END SUBROUTINE PARAM_MFSHALL_GOTO_MODEL
diff --git a/src/MNH/modn_param_mfshalln.f90 b/src/MNH/modn_param_mfshalln.f90
index 8722b0074..5d6b83205 100644
--- a/src/MNH/modn_param_mfshalln.f90
+++ b/src/MNH/modn_param_mfshalln.f90
@@ -3,6 +3,9 @@
!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
!MNH_LIC for details. version 1.
!-----------------------------------------------------------------
+!! MODIFICATIONS
+!! -------------
+!! 10/16 R.Honnert Update with AROME
!-----------------------------------------------------------------
! ###########################
MODULE MODN_PARAM_MFSHALL_n
@@ -39,7 +42,9 @@ USE MODD_PARAM_MFSHALL_n, ONLY: &
XA1_n => XA1, &
XB_n => XB, &
XC_n => XC, &
- XBETA1_n => XBETA1
+ XBETA1_n => XBETA1, &
+ XR_n => XR, &
+ XLAMBDA_n => XLAMBDA
!
IMPLICIT NONE
!
@@ -81,11 +86,16 @@ REAL,SAVE :: XB
REAL,SAVE :: XC
REAL,SAVE :: XBETA1
+! Tuning variables for RAHA updraft :
+
+REAL,SAVE :: XR
+REAL,SAVE :: XLAMBDA
+
NAMELIST/NAM_PARAM_MFSHALLn/XIMPL_MF,CMF_UPDRAFT,CMF_CLOUD,LMIXUV,LMF_FLX,&
XALP_PERT,XABUO,XBENTR,XBDETR,XCMF,XENTR_MF,&
XCRAD_MF,XENTR_DRY,XDETR_DRY,XDETR_LUP,XKCF_MF,&
XKRC_MF,XTAUSIGMF,XPRES_UV,XALPHA_MF,XSIGMA_MF,&
- XFRAC_UP_MAX,XA1,XB,XC,XBETA1
+ XFRAC_UP_MAX,XA1,XB,XC,XBETA1,XR,XLAMBDA
!
@@ -118,6 +128,9 @@ SUBROUTINE INIT_NAM_PARAM_MFSHALLn
XB = XB_n
XC = XC_n
XBETA1 = XBETA1_n
+ XR = XR_n
+ XLAMBDA = XLAMBDA_n
+
END SUBROUTINE INIT_NAM_PARAM_MFSHALLn
SUBROUTINE UPDATE_NAM_PARAM_MFSHALLn
@@ -147,6 +160,8 @@ SUBROUTINE UPDATE_NAM_PARAM_MFSHALLn
XB_n = XB
XC_n = XC
XBETA1_n = XBETA1
+ XR_n = XR
+ XLAMBDA_n = XLAMBDA
END SUBROUTINE UPDATE_NAM_PARAM_MFSHALLn
END MODULE MODN_PARAM_MFSHALL_n
diff --git a/src/MNH/shallow_mf.f90 b/src/MNH/shallow_mf.f90
index 1bf2aaa96..f30b6b2b0 100644
--- a/src/MNH/shallow_mf.f90
+++ b/src/MNH/shallow_mf.f90
@@ -161,7 +161,9 @@ 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 : EDKF gray zone
+!! R.Honnert 07/2012 : MF gray zone
+!! R.Honnert 10/2016 : SURF=gray zone initilisation + EDKF
+!! R.Honnert 10/2016 : Update with Arome
!! --------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -174,8 +176,10 @@ 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
@@ -284,6 +288,8 @@ 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
!------------------------------------------------------------------------
@@ -294,8 +300,9 @@ 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' ) THEN
+IF (HMF_UPDRAFT == 'EDKF' .OR. HMF_UPDRAFT == 'HRIO' .OR. &
+ HMF_UPDRAFT == 'RHCJ' .OR. HMF_UPDRAFT == 'BOUT' .OR. &
+ HMF_UPDRAFT == 'SURF' ) THEN
PENTR = 1.E20
PDETR = 1.E20
PEMF = 1.E20
@@ -317,9 +324,10 @@ ZTHVM(:,:) = PTHM(:,:)*((1.+XRV / XRD *PRM(:,:,1))/(1.+ZRTM(:,:)))
!!! 2. Compute updraft
!!! ---------------
!
-IF (HMF_UPDRAFT == 'EDKF' .OR. HMF_UPDRAFT == 'BOUT') THEN
+IF (HMF_UPDRAFT == 'EDKF' .OR. HMF_UPDRAFT == 'BOUT' .OR. HMF_UPDRAFT == 'SURF') THEN
GENTR_DETR = .TRUE.
- CALL COMPUTE_UPDRAFT(KKA,IKB,IKE,KKU,KKL,HFRAC_ICE,GENTR_DETR,OMIXUV,&
+ CALL COMPUTE_UPDRAFT(KKA,IKB,IKE,KKU,KKL,HFRAC_ICE,HMF_UPDRAFT,&
+ GENTR_DETR, OMIXUV, &
ONOMIXLG,KSV_LGBEG,KSV_LGEND, &
PZZ,PDZZ, &
PSFTH,PSFRV,PPABSM,PRHODREF, &
@@ -329,18 +337,36 @@ IF (HMF_UPDRAFT == 'EDKF' .OR. HMF_UPDRAFT == 'BOUT') THEN
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
- GENTR_DETR = .TRUE.
- CALL COMPUTE_UPDRAFT_RHCJ10(KKA,IKB,IKE,KKU,KKL,HFRAC_ICE,GENTR_DETR,OMIXUV,&
+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 )
+ 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, &
+ 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 )
ELSEIF (HMF_UPDRAFT == 'DUAL') THEN
!Updraft characteristics are already computed and received by interface
ELSEIF (HMF_UPDRAFT == 'HRIO') THEN
@@ -385,7 +411,9 @@ 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') THEN
+IF(HMF_UPDRAFT == 'EDKF' .OR. HMF_UPDRAFT == 'RHCJ'.OR. HMF_UPDRAFT == 'BOUT' &
+ .OR. HMF_UPDRAFT == 'SURF') THEN
+ IF ( PIMPL_MF > 1.E-10 ) THEN
CALL MF_TURB(KKA, IKB, IKE, KKU, KKL, OMIXUV, &
ONOMIXLG,KSV_LGBEG,KSV_LGEND, &
PIMPL_MF, PTSTEP, &
@@ -396,6 +424,14 @@ ZEMF_O_RHODREF=PEMF/PRHODREF
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, &
@@ -421,7 +457,8 @@ ZEMF_O_RHODREF=PEMF/PRHODREF
PTHVREF(:,JK)=RESHAPE(XTHVREF(:,:,JK),(/SIZE(PTHM,1)*SIZE(PTHM,2)/) )
ENDDO
ZG_O_THVREF=XG/PTHVREF
- CALL COMPUTE_BL89_ML(KKA,IKB,IKE,KKU,KKL,PDZZ,PTKEM,ZG_O_THVREF,ZTHVM,IKB,.TRUE.,ZLUP)
+ GLMIX=.TRUE.
+ CALL COMPUTE_BL89_ML(KKA,IKB,IKE,KKU,KKL,PDZZ,PTKEM(:,IKB) ,ZG_O_THVREF(:,IKB),ZTHVM,IKB,GLMIX,.TRUE.,ZLUP)
!! calcul de Dx/(h+hc)
DO JI=1,SIZE(XDXHAT)
DO JJ=1,SIZE(XDYHAT)
--
GitLab