From 18b976b420cd5950e24915cf40ed2ecc7a2c67c7 Mon Sep 17 00:00:00 2001
From: Gaelle DELAUTIER <gaelle.delautier@meteo.fr>
Date: Mon, 18 Sep 2017 16:09:03 +0200
Subject: [PATCH] G.Delautier 18/09/2017 : modif R.Honnert compute_updraft (add
RAHA)
---
src/MNH/compute_bl89_ml.f90 | 128 ++++--
src/MNH/compute_entr_detr.f90 | 684 ++++++++++++++---------------
src/MNH/compute_updraft.f90 | 108 +++--
src/MNH/compute_updraft_hrio.f90 | 7 +-
src/MNH/compute_updraft_raha.f90 | 666 ++++++++++++++++++++++++++++
src/MNH/compute_updraft_rhcj10.f90 | 4 +-
src/MNH/default_desfmn.f90 | 3 +
src/MNH/shallow_mf.f90 | 47 +-
8 files changed, 1210 insertions(+), 437 deletions(-)
create mode 100644 src/MNH/compute_updraft_raha.f90
diff --git a/src/MNH/compute_bl89_ml.f90 b/src/MNH/compute_bl89_ml.f90
index 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..b7346fad9 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,JLOOP
!----------------------------------------------------------------------------------
@@ -220,269 +223,260 @@ 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)
-
- ! Compute pressure at flux level KK+KKL
+ ZPRE(:)=PPRE_MINUS_HALF(:)
+ ZMIXTHL(:)=0.1
+ ZMIXRT(:)=0.1
+
+! 1.4 Estimation of PPART_DRY
+ DO JLOOP=1,SIZE(OTEST)
+ IF(OTEST(JLOOP) .AND. OTESTLCL(JLOOP)) THEN
+ !No dry part when condensation level is reached
+ PPART_DRY(JLOOP)=0.
+ ZDZ_STOP(JLOOP)=0.
+ ZPRE(JLOOP)=PPRE_MINUS_HALF(JLOOP)
+ ELSE IF (OTEST(JLOOP) .AND. .NOT. OTESTLCL(JLOOP)) THEN
+ !Temperature at flux level KK
+ ZT(JLOOP)=PTH_UP(JLOOP)*(PPRE_MINUS_HALF(JLOOP)/XP00) ** (XRD/XCPD)
+ !Saturating vapor pressure at flux level KK
+ ZFOESW(JLOOP) = MIN(EXP( XALPW - XBETAW/ZT(JLOOP) - XGAMW*LOG(ZT(JLOOP)) ), 0.99*PPRE_MINUS_HALF(JLOOP))
+ ZFOESI(JLOOP) = MIN(EXP( XALPI - XBETAI/ZT(JLOOP) - XGAMI*LOG(ZT(JLOOP)) ), 0.99*PPRE_MINUS_HALF(JLOOP))
+ !Computation of d.Rsat / dP (partial derivations with respect to P and T
+ !and use of T=Theta*(P/P0)**(R/Cp) to transform dT into dP with theta_up
+ !constant at the vertical)
+ ZDRSATODP(JLOOP)=(XBETAW/ZT(JLOOP)-XGAMW)*(1-ZFRAC_ICE(JLOOP))+(XBETAI/ZT(JLOOP)-XGAMI)*ZFRAC_ICE(JLOOP)
+ ZDRSATODP(JLOOP)=((XRD/XCPD)*ZDRSATODP(JLOOP)-1.)*PRSAT_UP(JLOOP)/ &
+ &(PPRE_MINUS_HALF(JLOOP)-(ZFOESW(JLOOP)*(1-ZFRAC_ICE(JLOOP)) + ZFOESI(JLOOP)*ZFRAC_ICE(JLOOP)))
+ !Use of d.Rsat / dP and pressure at flux level KK to find pressure (ZPRE)
+ !where Rsat is equal to PRT_UP
+ ZPRE(JLOOP)=PPRE_MINUS_HALF(JLOOP)+(PRT_UP(JLOOP)-PRSAT_UP(JLOOP))/ZDRSATODP(JLOOP)
+ !Fraction of dry part (computed with pressure and used with heights, no
+ !impact found when using log function here and for pressure on flux levels
+ !computation)
+ PPART_DRY(JLOOP)=MAX(0., MIN(1., (PPRE_MINUS_HALF(JLOOP)-ZPRE(JLOOP))/(PPRE_MINUS_HALF(JLOOP)-PPRE_PLUS_HALF(JLOOP))))
+ !Height above flux level KK of the cloudy part
+ ZDZ_STOP(JLOOP) = (PZZ(JLOOP,KK+KKL)-PZZ(JLOOP,KK))*PPART_DRY(JLOOP)
+ END IF
+ END DO
+
+! 1.5 Gradient and flux values of thetav
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
+! --------------------------------------------------------------------
+
+DO JLOOP=1,SIZE(OTEST)
+ IF (OTEST(JLOOP) .AND. PPART_DRY(JLOOP)>0.) THEN
+ ZDZ(JLOOP)=MIN(ZDZ_STOP(JLOOP),(PZZ(JLOOP,KK+KKL)-PZZ(JLOOP,KK))*0.5)
+ PBUO_INTEG_DRY(JLOOP) = ZG_O_THVREF(JLOOP)*ZDZ(JLOOP)*&
+ (0.5 * ( - ZCOEFF_MINUS_HALF(JLOOP))*ZDZ(JLOOP) &
+ - ZTHV_MINUS_HALF(JLOOP) + PTHV_UP(JLOOP) )
+
+ !Between mass flux KK and bottom of cloudy part (if above mass flux)
+ ZDZ(JLOOP)=MAX(0., ZDZ_STOP(JLOOP)-(PZZ(JLOOP,KK+KKL)-PZZ(JLOOP,KK))*0.5)
+ PBUO_INTEG_DRY(JLOOP) = PBUO_INTEG_DRY(JLOOP) + ZG_O_THVREF(JLOOP)*ZDZ(JLOOP)*&
+ (0.5 * ( - ZCOEFF_PLUS_HALF(JLOOP))*ZDZ(JLOOP) &
+ - PTHVM(JLOOP,KK) + PTHV_UP(JLOOP) )
+ IF (PBUO_INTEG_DRY(JLOOP)>=0.) THEN
+ PENTR(JLOOP) = 0.5/(XABUO-XBENTR*XENTR_DRY)*&
+ LOG(1.+ (2.*(XABUO-XBENTR*XENTR_DRY)/PW_UP2(JLOOP,KK))* &
+ PBUO_INTEG_DRY(JLOOP))
+ PDETR(JLOOP) = 0.
+ ELSE
+ PENTR(JLOOP) = 0.
+ PDETR(JLOOP) = 0.5/(XABUO)*&
+ LOG(1.+ (2.*(XABUO)/PW_UP2(JLOOP,KK))* &
+ (-PBUO_INTEG_DRY(JLOOP)))
+ ENDIF
+ PENTR(JLOOP) = XENTR_DRY*PENTR(JLOOP)/(PZZ(JLOOP,KK+KKL)-PZZ(JLOOP,KK))
+ PDETR(JLOOP) = XDETR_DRY*PDETR(JLOOP)/(PZZ(JLOOP,KK+KKL)-PZZ(JLOOP,KK))
+ !Minimum value of detrainment
+ ZWK(JLOOP)=PLUP(JLOOP)-0.5*(PZZ(JLOOP,KK)+PZZ(JLOOP,KK+KKL))
+ ZWK(JLOOP)=SIGN(MAX(1., ABS(ZWK(JLOOP))), ZWK(JLOOP)) ! ZWK must not be zero
+ PDETR(JLOOP) = MAX(PPART_DRY(JLOOP)*XDETR_LUP/ZWK(JLOOP), PDETR(JLOOP))
+ ELSE
+ !No dry part, consation reached (OTESTLCL)
+ PBUO_INTEG_DRY(JLOOP) = 0.
+ PENTR(JLOOP)=0.
+ PDETR(JLOOP)=0.
+ END IF
+ENDDO
+
+
+! 3 Wet part computation
+! -----------------------
+
+! 3.1 Integral buoyancy for cloudy part
- ! Compute 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.
-! -----------------------------------------------------
-
-
-! 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(:))
-
- ! Compute theta_v of updraft at flux level KK+KKL
- ZRCMIX_F2=PRC_UP
- ZRIMIX_F2=PRI_UP
+ ZTHV_UP_F2(:) = ZTHMIX(:)*(1.+ZRVORD*ZRVMIX(:))/(1.+PRT_UP(:))
+
+ ! Integral buoyancy for cloudy part
+ DO JLOOP=1,SIZE(OTEST)
+ IF(OTEST(JLOOP) .AND. PPART_DRY(JLOOP)<1.) THEN
+ !Gradient of Theta V updraft over the cloudy part, assuming that thetaV updraft don't change
+ !between flux level KK and bottom of cloudy part
+ ZCOTHVU(JLOOP)=(ZTHV_UP_F2(JLOOP)-PTHV_UP(JLOOP))/((PZZ(JLOOP,KK+KKL)-PZZ(JLOOP,KK))*(1-PPART_DRY(JLOOP)))
+
+ !Computation in two parts to use change of gradient of theta v of environment
+ !Between bottom of cloudy part (if under mass level) and mass level KK
+ ZDZ(JLOOP)=MAX(0., 0.5*(PZZ(JLOOP,KK+KKL)-PZZ(JLOOP,KK))-ZDZ_STOP(JLOOP))
+ PBUO_INTEG_CLD(JLOOP) = ZG_O_THVREF(JLOOP)*ZDZ(JLOOP)*&
+ (0.5*( ZCOTHVU(JLOOP) - ZCOEFF_MINUS_HALF(JLOOP))*ZDZ(JLOOP) &
+ - (PTHVM(JLOOP,KK)-ZDZ(JLOOP)*ZCOEFF_MINUS_HALF(JLOOP)) + PTHV_UP(JLOOP) )
+
+ !Between max(mass level, bottom of cloudy part) and flux level KK+KKL
+ ZDZ(JLOOP)=(PZZ(JLOOP,KK+KKL)-PZZ(JLOOP,KK))-MAX(ZDZ_STOP(JLOOP),0.5*(PZZ(JLOOP,KK+KKL)-PZZ(JLOOP,KK)))
+ PBUO_INTEG_CLD(JLOOP) = PBUO_INTEG_CLD(JLOOP)+ZG_O_THVREF(JLOOP)*ZDZ(JLOOP)*&
+ (0.5*( ZCOTHVU(JLOOP) - ZCOEFF_PLUS_HALF(JLOOP))*ZDZ(JLOOP)&
+ - (PTHVM(JLOOP,KK)+(0.5*((PZZ(JLOOP,KK+KKL)-PZZ(JLOOP,KK)))-ZDZ(JLOOP))*ZCOEFF_PLUS_HALF(JLOOP)) +&
+ PTHV_UP(JLOOP) )
+
+ ELSE
+ !No cloudy part
+ PBUO_INTEG_CLD(JLOOP)=0.
+ END IF
+ END DO
+
+! 3.2 Critical mixed fraction for KK+KKL flux level (ZKIC_F2) and
+! for bottom of cloudy part (ZKIC), then a mean for the cloudy part
+! (put also in ZKIC)
+!
+! computation by estimating unknown
+! T^mix r_c^mix and r_i^mix from enthalpy^mix and r_w^mix
+! We determine the zero crossing of the linear curve
+! evaluating the derivative using ZMIXF=0.1
+
+ ZKIC_INIT=0.1 ! starting value for critical mixed fraction for CLoudy Part
+
+ ! Compute thetaV of environment at the bottom of cloudy part
+ ! and cons then non cons. var. of mixture at the bottom of cloudy part
+
+ ! JI computed to avoid KKL(KK-KKL) being < KKL*KKB
+ JI=KKL*MAX(KKL*(KK-KKL),KKL*KKB)
+ DO JLOOP=1,SIZE(OTEST)
+ IF(OTEST(JLOOP) .AND. PPART_DRY(JLOOP)>0.5) THEN
+ ZDZ(JLOOP)=ZDZ_STOP(JLOOP)-0.5*(PZZ(JLOOP,KK+KKL)-PZZ(JLOOP,KK))
+ ZTHV(JLOOP)= PTHVM(JLOOP,KK)+ZCOEFF_PLUS_HALF(JLOOP)*ZDZ(JLOOP)
+ ZMIXTHL(JLOOP) = ZKIC_INIT * &
+ (PTHLM(JLOOP,KK)+ZDZ(JLOOP)*(PTHLM(JLOOP,KK+KKL)-PTHLM(JLOOP,KK))/PDZZ(JLOOP,KK+KKL)) + &
+ (1. - ZKIC_INIT)*PTHL_UP(JLOOP)
+ ZMIXRT(JLOOP) = ZKIC_INIT * &
+ (PRTM(JLOOP,KK)+ZDZ(JLOOP)*(PRTM(JLOOP,KK+KKL)-PRTM(JLOOP,KK))/PDZZ(JLOOP,KK+KKL)) + &
+ (1. - ZKIC_INIT)*PRT_UP(JLOOP)
+ ELSEIF(OTEST(JLOOP)) THEN
+ ZDZ(JLOOP)=0.5*(PZZ(JLOOP,KK+KKL)-PZZ(JLOOP,KK))-ZDZ_STOP(JLOOP)
+ ZTHV(JLOOP)= PTHVM(JLOOP,KK)-ZCOEFF_MINUS_HALF(JLOOP)*ZDZ(JLOOP)
+ ZMIXTHL(JLOOP) = ZKIC_INIT * &
+ (PTHLM(JLOOP,KK)-ZDZ(JLOOP)*(PTHLM(JLOOP,KK)-PTHLM(JLOOP,JI))/PDZZ(JLOOP,KK)) + &
+ (1. - ZKIC_INIT)*PTHL_UP(JLOOP)
+ ZMIXRT(JLOOP) = ZKIC_INIT * &
+ (PRTM(JLOOP,KK)-ZDZ(JLOOP)*(PRTM(JLOOP,KK)-PRTM(JLOOP,JI))/PDZZ(JLOOP,KK)) + &
+ (1. - ZKIC_INIT)*PRT_UP(JLOOP)
+ ENDIF
+ ENDDO
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
- ENDWHERE
-
- IF(KK/=KKB)THEN
- ZCOEFF_MINUS_HALF = ((PTHVM(:,KK)-PTHVM(:,KK-KKL))/PDZZ(:,KK))
- ELSE
- ZCOEFF_MINUS_HALF = 0.
- ENDIF
- ZCOEFF_PLUS_HALF = ((PTHVM(:,KK+KKL)-PTHVM(:,KK))/PDZZ(:,KK+KKL))
-
- ZCOTHVU_MINUS_HALF = (ZTHV_UP(:)-ZTHV_UP_F1(:))/ZDZ_HALF(:)
- ZCOTHVU_PLUS_HALF = (ZTHV_UP_F2(:)-ZTHV_UP(:))/ZDZ_HALF(:)
-
- IF(KK/=KKB)THEN
- ZTHV_MINUS_HALF = ZCOEFF_MINUS_HALF*0.5*(PZZ(:,KK)-PZZ(:,KK-KKL))+PTHVM(:,KK-KKL)
- ZTHV_PLUS_HALF = ZCOEFF_PLUS_HALF*0.5*(PZZ(:,KK)-PZZ(:,KK-KKL))+ ZTHV_MINUS_HALF
- ELSE
- ZTHV_MINUS_HALF = PTHVM(:,KK)
- ZTHV_PLUS_HALF = ZCOEFF_PLUS_HALF*0.5*(PZZ(:,KK+KKL)-PZZ(:,KK))+ ZTHV_MINUS_HALF !according to PZZ computation at KKB-KKL
- ENDIF
-
- ! Integral Buoyancy between flux level KK and mass level KK
- PBUO_INTEG = ZG_O_THVREF(:,KK)*ZDZ_HALF(:)*&
- (0.5*( ZCOTHVU_MINUS_HALF - ZCOEFF_MINUS_HALF)*ZDZ_HALF(:) &
- - ZTHV_MINUS_HALF + ZTHV_UP_F1(:) )
-
- WHERE ((OTEST).AND.(.NOT.OTESTLCL))
- PENTR=0.
- PDETR=0.
-
- ZBUO_INTEG = ZG_O_THVREF(:,KK)*ZDZ_STOP(:)*&
- (0.5 * ( - ZCOEFF_MINUS_HALF)* ZDZ_STOP(:) &
- - ZTHV_MINUS_HALF + ZTHV_UP_F1(:) )
- WHERE (ZBUO_INTEG(:)>=0.)
- PENTR = 0.5/(XABUO-XBENTR*XENTR_DRY)*&
- LOG(1.+ (2.*(XABUO-XBENTR*XENTR_DRY)/PW_UP2(:,KK))* &
- ZBUO_INTEG)
- PDETR = 0.
-
- ZW2_HALF = PW_UP2(:,KK) + 2*(XABUO-XBENTR*XENTR_DRY)*(ZBUO_INTEG)
- ELSEWHERE
- PENTR = 0.
- PDETR = 0.5/(XABUO)*&
- LOG(1.+ (2.*(XABUO)/PW_UP2(:,KK))* &
- MAX(0.,-ZBUO_INTEG))
-
- ZW2_HALF = PW_UP2(:,KK) + 2*(XABUO)*(ZBUO_INTEG)
- ENDWHERE
- ENDWHERE
-
-
- ZDZ_STOP(:) = ZDZ_HALF(:)
-
-! total Integral Buoyancy between flux level KK and flux level KK+KKL
- PBUO_INTEG = PBUO_INTEG + ZG_O_THVREF(:,KK)*ZDZ_HALF(:)*&
- (0.5*(ZCOTHVU_PLUS_HALF - ZCOEFF_PLUS_HALF)* ZDZ_HALF(:) - &
- PTHVM(:,KK) + ZTHV_UP(:) )
-
- IF(KK*KKL<(KKE-KKL)*KKL) THEN !Computation only if we are strictly below KKE-KKL
- WHERE ((((ZRCMIX_F2(:)+ZRIMIX_F2(:)>0.).AND.(ZRCMIX(:)+ZRIMIX(:)<=0.)).AND.(.NOT.OTESTLCL)).AND.(.NOT.GTEST_LOCAL_LCL(:)))
- ! If rc and/or ri is found between mass level KK and flux level KK+KKL
- ! a part of dry entrainment is defined
- ! the exact height of LCL is also determined
- ZCOEFF_QSAT(:) = ((ZRCMIX_M2(:)+ZRIMIX_M2(:)) - (ZRCMIX_F2(:)+ZRIMIX_F2(:))) / &
- & (0.5* (PZZ(:,KK+2*KKL)-PZZ(:,KK+KKL)))
- !old formulation without ice (and perhaps with errors)
- !ZCOEFF_QSAT(:) = (PRT_UP(:) - &
- ! QSAT(ZTHMIX_F2(:)*((PPABSM(:,KK+KKL)/XP00)**(XRD/XCPD)),&
- ! PPABSM(:,KK+KKL)) - &
- ! ZRCMIX(:))/ (0.5* (PZZ(:,KK+2*KKL)-PZZ(:,KK+KKL)))
- WHERE ((ZCOEFF_QSAT(:)>0.) .OR. (ZCOEFF_QSAT(:)<0.))
- ZRC_ORD(:) = ZRCMIX_F2(:)+ZRIMIX_F2(:) - ZCOEFF_QSAT(:) * ZDZ_HALF(:)
- ZDZ_STOP = (- ZRC_ORD(:)/ZCOEFF_QSAT(:))
- ZPART_DRY(:) = 0.5+MAX(MIN(ZDZ_STOP / (PZZ(:,KK+KKL)-PZZ(:,KK)),0.5),0.)
- GTEST_LOCAL_LCL2(:)=.TRUE.
- ENDWHERE
- ENDWHERE
- ENDIF
-
- WHERE (((OTEST).AND.(.NOT.OTESTLCL)).AND.(.NOT.GTEST_LOCAL_LCL(:)))
- ZBUO_INTEG = ZG_O_THVREF(:,KK)*ZDZ_STOP(:)*&
- (0.5*( - ZCOEFF_PLUS_HALF)* ZDZ_STOP(:)&
- - PTHVM(:,KK) + ZTHV_UP(:) )
-
- WHERE (ZW2_HALF>0.)
- WHERE (ZBUO_INTEG(:)>=0.)
- PENTR = PENTR + 0.5/(XABUO-XBENTR*XENTR_DRY)* &
- LOG(1.+ (2.*(XABUO-XBENTR*XENTR_DRY)/ZW2_HALF(:)) * ZBUO_INTEG)
-
- PDETR = PDETR
- ELSEWHERE
- PENTR = PENTR
- PDETR = PDETR + 0.5/(XABUO)* &
- LOG(1.+ (2.*(XABUO)/ZW2_HALF(:)) * &
- MAX(-ZBUO_INTEG,0.))
- ENDWHERE
- ELSEWHERE
- ! if w**2<0 the updraft is stopped
- OTEST=.FALSE.
- PENTR = PENTR
- PDETR = PDETR
- ENDWHERE
- ENDWHERE
- PENTR = XENTR_DRY*PENTR/(PZZ(:,KK+KKL)-PZZ(:,KK))
- PDETR = XDETR_DRY*PDETR/(PZZ(:,KK+KKL)-PZZ(:,KK))
-
- ZWK(:)=PLUP(:) - 0.5*(PZZ(:,KK)+PZZ(:,KK+KKL))
- ZWK(:)=SIGN(MAX(1., ABS(ZWK(:))), ZWK(:))
- PDETR(:)=MAX(ZPART_DRY(:)*XDETR_LUP/ZWK(:),PDETR(:))
-
-! compute final value of critical mixed fraction using theta_v
-! of mixture, grid-scale and updraft in cloud
- WHERE ((OTEST).AND.(OTESTLCL))
- ZKIC(:) = MAX(0.,ZTHV_UP(:)-PTHVM(:,KK))*ZKIC(:) / &
- (ZTHV_UP(:)-ZTHVMIX_1(:)+XMNH_EPSILON)
-
- ZKIC(:) = MAX(0., MIN(1., ZKIC(:)))
-
- ZEPSI(:) = ZKIC(:) **2.
- ZDELTA(:) = (1.-ZKIC(:))**2.
- ZEPSI_CLOUD=MIN(ZDELTA,ZEPSI)
- ZCOEFFMF_CLOUD(:)=XENTR_MF * XG / XCRAD_MF
- PENTR(:) = ZCOEFFMF_CLOUD(:)*ZEPSI_CLOUD(:)
- PDETR(:) = ZCOEFFMF_CLOUD(:)*ZDELTA(:)
- ENDWHERE
-
-! compute final value of critical mixed fraction using theta_v
-! of mixture, grid-scale and updraft in cloud
- WHERE (((OTEST).AND.(.NOT.(OTESTLCL))).AND.((GTEST_LOCAL_LCL(:).OR.GTEST_LOCAL_LCL2(:))))
- ZKIC(:) = MAX(0.,ZTHV_UP_F2(:)-ZTHV_PLUS_HALF)*ZKIC(:) / &
- (ZTHV_UP_F2(:)-ZTHVMIX_2(:)+XMNH_EPSILON)
- ZKIC(:) = MAX(0., MIN(1., ZKIC(:)))
- ZEPSI(:) = ZKIC(:) **2.
- ZDELTA(:) = (1.-ZKIC(:))**2.
- ZEPSI_CLOUD=MIN(ZDELTA,ZEPSI)
- ZCOEFFMF_CLOUD(:)=XENTR_MF * XG / XCRAD_MF
- PENTR(:) = PENTR+(1.-ZPART_DRY(:))*ZCOEFFMF_CLOUD(:)*ZEPSI_CLOUD(:)
- PDETR(:) = PDETR+(1.-ZPART_DRY(:))*ZCOEFFMF_CLOUD(:)*ZDELTA(:)
- ENDWHERE
+ ZTHVMIX_F2(:) = ZTHMIX(:)*(1.+ZRVORD*ZRVMIX(:))/(1.+ZMIXRT(:))
+
+ !Computation of mean ZKIC over the cloudy part
+ DO JLOOP=1,SIZE(OTEST)
+ IF (OTEST(JLOOP)) THEN
+ ! Compute ZKIC at the bottom of cloudy part
+ ! Thetav_up at bottom is equal to Thetav_up at flux level KK
+ IF (ABS(PTHV_UP(JLOOP)-ZTHVMIX(JLOOP))<1.E-10) THEN
+ ZKIC(JLOOP)=1.
+ ELSE
+ ZKIC(JLOOP) = MAX(0.,PTHV_UP(JLOOP)-ZTHV(JLOOP))*ZKIC_INIT / &
+ (PTHV_UP(JLOOP)-ZTHVMIX(JLOOP))
+ END IF
+ ! Compute ZKIC_F2 at flux level KK+KKL
+ IF (ABS(ZTHV_UP_F2(JLOOP)-ZTHVMIX_F2(JLOOP))<1.E-10) THEN
+ ZKIC_F2(JLOOP)=1.
+ ELSE
+ ZKIC_F2(JLOOP) = MAX(0.,ZTHV_UP_F2(JLOOP)-ZTHV_PLUS_HALF(JLOOP))*ZKIC_INIT / &
+ (ZTHV_UP_F2(JLOOP)-ZTHVMIX_F2(JLOOP))
+ END IF
+ !Mean ZKIC over the cloudy part
+ ZKIC(JLOOP)=MAX(MIN(0.5*(ZKIC(JLOOP)+ZKIC_F2(JLOOP)),1.),0.)
+ END IF
+ END DO
+
+
+! 3.3 Integration of PDF
+! According to Kain and Fritsch (1990), we replace delta Mt
+! in eq. (7) and (8) using eq. (5). Here we compute the ratio
+! of integrals without computing delta Me
+
+ !Constant PDF
+ !For this PDF, eq. (5) is delta Me=0.5*delta Mt
+ DO JLOOP=1,SIZE(OTEST)
+ IF(OTEST(JLOOP)) THEN
+ ZEPSI(JLOOP) = ZKIC(JLOOP)**2. !integration multiplied by 2
+ ZDELTA(JLOOP) = (1.-ZKIC(JLOOP))**2. !idem
+ ENDIF
+ ENDDO
+
+ !Triangular PDF
+ !Calculus must be verified before activating this part, but in this state,
+ !results on ARM case are almost identical
+ !For this PDF, eq. (5) is also delta Me=0.5*delta Mt
+ !WHERE(OTEST)
+ ! !Integration multiplied by 2
+ ! WHERE(ZKIC<0.5)
+ ! ZEPSI(:)=8.*ZKIC(:)**3/3.
+ ! ZDELTA(:)=1.-4.*ZKIC(:)**2+8.*ZKIC(:)**3/3.
+ ! ELSEWHERE
+ ! ZEPSI(:)=5./3.-4*ZKIC(:)**2+8.*ZKIC(:)**3/3.
+ ! ZDELTA(:)=8.*(1.-ZKIC(:))**3/3.
+ ! ENDWHERE
+ !ENDWHERE
+
+! 3.4 Computation of PENTR and PDETR
+ DO JLOOP=1,SIZE(OTEST)
+ IF(OTEST(JLOOP)) THEN
+ ZEPSI_CLOUD=MIN(ZDELTA,ZEPSI)
+ PENTR_CLD(JLOOP) = (1.-PPART_DRY(JLOOP))*ZCOEFFMF_CLOUD*PRHODREF(JLOOP)*ZEPSI_CLOUD(JLOOP)
+ PDETR_CLD(JLOOP) = (1.-PPART_DRY(JLOOP))*ZCOEFFMF_CLOUD*PRHODREF(JLOOP)*ZDELTA(JLOOP)
+ PENTR(JLOOP) = PENTR(JLOOP)+PENTR_CLD(JLOOP)
+ PDETR(JLOOP) = PDETR(JLOOP)+PDETR_CLD(JLOOP)
+ ELSE
+ PENTR_CLD(JLOOP) = 0.
+ PDETR_CLD(JLOOP) = 0.
+ ENDIF
+ ENDDO
END SUBROUTINE COMPUTE_ENTR_DETR
diff --git a/src/MNH/compute_updraft.f90 b/src/MNH/compute_updraft.f90
index 1d0e2a821..e41658775 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,6 +162,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
@@ -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
@@ -239,14 +248,17 @@ LOGICAL :: GLMIX
LOGICAL, DIMENSION(SIZE(PTHM,1)) :: GWORK1
LOGICAL, DIMENSION(SIZE(PTHM,1),SIZE(PTHM,2)) :: GWORK2
-INTEGER :: ITEST
+INTEGER :: ITEST,JLOOP
-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 +275,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)
@@ -367,7 +378,7 @@ 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,KKB,GLMIX,.FALSE.,ZLUP)
ZLUP(:)=MAX(ZLUP(:),1.E-10)
! Compute Buoyancy flux at the ground
@@ -375,8 +386,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.
@@ -400,7 +417,6 @@ GTESTLCL(:)=.FALSE.
GTESTETL(:)=.FALSE.
! Loop on vertical level
-
DO JK=KKB,KKE-KKL,KKL
! IF the updraft top is reached for all column, stop the loop on levels
@@ -418,7 +434,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 +441,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
@@ -457,15 +478,24 @@ 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) !&
+! WHERE(GTEST)
+ DO JLOOP=1,SIZE(GTEST)
+ IF (GTEST(JLOOP) ) THEN
+ ZMIX2(JLOOP) = (PZZ(JLOOP,JK+KKL)-PZZ(JLOOP,JK))*PENTR(JLOOP,JK) !&
+ ZMIX3_CLD(JLOOP) = (PZZ(JLOOP,JK+KKL)-PZZ(JLOOP,JK))*(1.-ZPART_DRY(JLOOP))*ZDETR_CLD(JLOOP,JK) !&
+ ZMIX2_CLD(JLOOP) = (PZZ(JLOOP,JK+KKL)-PZZ(JLOOP,JK))*(1.-ZPART_DRY(JLOOP))*ZENTR_CLD(JLOOP,JK)
- PTHL_UP(:,JK+KKL)=(PTHL_UP(:,JK)*(1.-0.5*ZMIX2(:)) + PTHLM(:,JK)*ZMIX2(:)) &
- /(1.+0.5*ZMIX2(:))
- PRT_UP(:,JK+KKL) =(PRT_UP (:,JK)*(1.-0.5*ZMIX2(:)) + PRTM(:,JK)*ZMIX2(:)) &
- /(1.+0.5*ZMIX2(:))
- ENDWHERE
+ !PTHL_UP(JLOOP,JK+KKL)=(PTHL_UP(JLOOP,JK)*(1.-0.5*ZMIX2(JLOOP)) + PTHLM(JLOOP,JK)*ZMIX2(JLOOP)) &
+ ! /(1.+0.5*ZMIX2(JLOOP))
+ !PRT_UP(JLOOP,JK+KKL) =(PRT_UP (JLOOP,JK)*(1.-0.5*ZMIX2(JLOOP)) + PRTM(JLOOP,JK)*ZMIX2(JLOOP)) &
+ ! /(1.+0.5*ZMIX2(JLOOP))
+
+ PTHL_UP(JLOOP,JK+KKL)=PTHL_UP(JLOOP,JK)*EXP(-ZMIX2(JLOOP)) + PTHLM(JLOOP,JK)*(1-EXP(-ZMIX2(JLOOP)))
+ PRT_UP(JLOOP,JK+KKL) =PRT_UP (JLOOP,JK)*EXP(-ZMIX2(JLOOP)) + PRTM(JLOOP,JK)*(1-EXP(-ZMIX2(JLOOP)))
+
+ END IF
+ END DO
+! ENDWHERE
IF(OMIXUV) THEN
@@ -522,20 +552,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 +627,13 @@ IF(OENTR_DETR) THEN
PDEPTH(JI) = MAX(0., PZZ(JI,KKCTL(JI)) - PZZ(JI,KKLCL(JI)) )
END DO
+ GWORK1(:)= (GTESTLCL(:) .AND. (PDEPTH(:) > ZDEPTH_MAX1) )
+ GWORK2(:,:) = SPREAD( GWORK1(:), DIM=2, NCOPIES=MAX(KKU,KKA) )
+ ZCOEF(:,:) = SPREAD( (1.-(PDEPTH(:)-ZDEPTH_MAX1)/(ZDEPTH_MAX2-ZDEPTH_MAX1)), DIM=2, NCOPIES=SIZE(ZCOEF,2))
+ ZCOEF=MIN(MAX(ZCOEF,0.),1.)
+ WHERE (GWORK2)
+ PEMF(:,:) = PEMF(:,:) * ZCOEF(:,:)
+ PFRAC_UP(:,:) = PFRAC_UP(:,:) * ZCOEF(:,:)
+ ENDWHERE
ENDIF
-
-GWORK1(:)= (GTESTLCL(:) .AND. (PDEPTH(:) > ZDEPTH_MAX1) )
-GWORK2(:,:) = SPREAD( GWORK1(:), DIM=2, NCOPIES=MAX(KKU,KKA) )
-ZCOEF(:,:) = SPREAD( (1.-(PDEPTH(:)-ZDEPTH_MAX1)/(ZDEPTH_MAX2-ZDEPTH_MAX1)), DIM=2, NCOPIES=SIZE(ZCOEF,2))
-
END SUBROUTINE COMPUTE_UPDRAFT
diff --git a/src/MNH/compute_updraft_hrio.f90 b/src/MNH/compute_updraft_hrio.f90
index 9fa0ad5f0..346644141 100644
--- a/src/MNH/compute_updraft_hrio.f90
+++ b/src/MNH/compute_updraft_hrio.f90
@@ -390,6 +390,9 @@ 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(:,:))
@@ -444,7 +447,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
@@ -845,5 +849,6 @@ 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..553a95de0
--- /dev/null
+++ b/src/MNH/compute_updraft_raha.f90
@@ -0,0 +1,666 @@
+!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
+!! --------------------------------------------------------------------------
+!
+!* 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..f062dcf0a 100644
--- a/src/MNH/compute_updraft_rhcj10.f90
+++ b/src/MNH/compute_updraft_rhcj10.f90
@@ -390,8 +390,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/default_desfmn.f90 b/src/MNH/default_desfmn.f90
index 618f2611a..e687622a5 100644
--- a/src/MNH/default_desfmn.f90
+++ b/src/MNH/default_desfmn.f90
@@ -1058,6 +1058,8 @@ XA1 = 2./3.
XB = 0.002
XC = 0.012
XBETA1 = 0.9
+XR = 2.
+XLAMBDA= 0.
!
!-------------------------------------------------------------------------------
!
@@ -1167,6 +1169,7 @@ LUSECHEM = .FALSE.
LUSECHAQ = .FALSE.
LUSECHIC = .FALSE.
LCH_INIT_FIELD = .FALSE.
+!LCH_SURFACE_FLUX = .FALSE.
LCH_CONV_SCAV = .FALSE.
LCH_CONV_LINOX = .FALSE.
LCH_PH = .FALSE.
diff --git a/src/MNH/shallow_mf.f90 b/src/MNH/shallow_mf.f90
index 1bf2aaa96..a294f5859 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,9 @@ 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, &
@@ -341,6 +348,21 @@ ELSEIF (HMF_UPDRAFT == 'RHCJ') 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 == '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 +407,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 +420,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 +453,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