diff --git a/src/MNH/compute_updraft_hrio.f90 b/src/MNH/compute_updraft_hrio.f90 deleted file mode 100644 index adccc8ca9fe3a2469fb42f78e4b15aa2268aef7a..0000000000000000000000000000000000000000 --- a/src/MNH/compute_updraft_hrio.f90 +++ /dev/null @@ -1,868 +0,0 @@ -!MNH_LIC Copyright 2004-2019 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 - MODULE MODI_COMPUTE_UPDRAFT_HRIO -! ########################### -! -INTERFACE -! -! ################################################################# - SUBROUTINE COMPUTE_UPDRAFT_HRIO(KKA,KKB,KKE,KKU,KKL, HFRAC_ICE, & - OENTR_DETR,OMIXUV, & - ONOMIXLG,KSV_LGBEG,KSV_LGEND, & - PZZ,PDZZ, & - PSFTH,PSFRV, & - PPABSM,PRHODREF,PUM,PVM,PTKEM,PWM,& - 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, & - PTHL_DO, PTHV_DO, PRT_DO, & - PU_DO, PV_DO, PSV_DO, & - 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(len=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) :: PWM ! w mean wind -! -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 - PRI_UP,PTHV_UP,& ! updraft ri, THv - PW_UP,PFRAC_UP,& ! updraft w, fraction - PFRAC_ICE_UP,& ! liquid/solid fraction in updraft - PRSAT_UP ! Rsat -REAL, DIMENSION(:,:), INTENT(INOUT):: PTHL_DO,PTHV_DO,PRT_DO,PU_DO,PV_DO - -REAL, DIMENSION(:,:,:), INTENT(OUT) :: PSV_UP ! updraft scalar var. -REAL, DIMENSION(:,:,:), INTENT(OUT) :: PSV_DO ! downdraft scalar var. - -REAL, DIMENSION(:,:), INTENT(INOUT):: PEMF,PDETR,PENTR ! Mass_flux, - ! entrainment, detrainment -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_HRIO - -END INTERFACE -! -END MODULE MODI_COMPUTE_UPDRAFT_HRIO! ######spl - SUBROUTINE COMPUTE_UPDRAFT_HRIO(KKA,KKB,KKE,KKU,KKL,HFRAC_ICE, & - OENTR_DETR,OMIXUV, & - ONOMIXLG,KSV_LGBEG,KSV_LGEND, & - PZZ,PDZZ, & - PSFTH,PSFRV, & - PPABSM,PRHODREF,PUM,PVM,PTKEM,PWM, & - 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, & - PTHL_DO, PTHV_DO, PRT_DO, & - PU_DO, PV_DO, PSV_DO, & - PEMF,PDETR,PENTR, & - PBUO_INTEG,KKLCL,KKETL,KKCTL, & - PDEPTH ) -! ################################################################# -!! -!!**** *COMPUTE_UPDRAFT_HRIO* - calculates caracteristics of the updraft -!! -!! -!! PURPOSE -!! ------- -!!**** The purpose of this routine is to build the updraft model -!! -! -!!** METHOD -!! ------ -!! -!! EXTERNAL -!! -------- -!! -!! IMPLICIT ARGUMENTS -!! ------------------ -!! -!! !! REFERENCE -!! --------- -!! Book 1 of Meso-NH documentation (chapter Turbulence) -!! Soares et al. 2004 QJ -!! -!! AUTHOR -!! ------ -!! J.Pergaud -!! V.Masson : Optimization 07/2010 -!! S. Riette : 07/2010 : modification for reproducibility -!! 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 -!! Q.Rodier 01/2019 : support RM17 mixing length -! P. Wautelet 12/04/2019: replace ABORT and STOP calls by Print_msg -!! -------------------------------------------------------------------------- -! -!* 0. DECLARATIONS -! ------------ -! -USE MODD_CST -USE MODD_PARAM_MFSHALL_n, ONLY : XPRES_UV,XALP_PERT,XCMF,XFRAC_UP_MAX,XA1,XB,& - XC,XBETA1 -!USE MODD_CMFSHALL, ONLY : XPRES_UV,XALP_PERT,XCMF,XFRAC_UP_MAX,XA1,XALPHA_SEUIL,LNORM_RESOL,& -! XCOEF1,XCOEF2,XCOEF3,NBUOY,XB,& -! XC,XBETA1 -USE MODD_GRID_n, ONLY : XDXHAT, XDYHAT -USE MODD_BLANK -USE MODD_TURB_n, ONLY :CTURBLEN - -use mode_msg - -!USE MODI_COMPUTE_ENTR_DETR -USE MODI_TH_R_FROM_THL_RT_1D -USE MODI_SHUMAN_MF - -USE MODI_COMPUTE_BL89_ML - -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(len=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 -! FAUX dans AROME, mais pas grave dans la param�trisation -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) :: PWM ! w mean wind -REAL, DIMENSION(:,:), INTENT(IN) :: PTKEM ! TKE 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 - PRI_UP,PTHV_UP,& ! updraft ri, THv - PW_UP,PFRAC_UP,& ! updraft w, fraction - PFRAC_ICE_UP,& ! liquid/solid fraction in updraft - PRSAT_UP ! Rsat -REAL, DIMENSION(:,:), INTENT(INOUT):: PTHL_DO,PTHV_DO,PRT_DO,PU_DO,PV_DO ! environment var. - -REAL, DIMENSION(:,:,:), INTENT(OUT) :: PSV_UP ! updraft scalar var. -REAL, DIMENSION(:,:,:), INTENT(OUT) :: PSV_DO ! environment 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 ! Theta,rv of - ! updraft environnement -REAL, DIMENSION(SIZE(PTHM,1),SIZE(PTHM,2)) :: & - ZRTM_F, ZTHLM_F, ZTKEM_F,& ! rt, thetal,TKE,pressure, - 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 -!================================================================================== -REAL, DIMENSION(SIZE(PTHM,1),SIZE(PTHM,2)) :: ZW_UP ! w of the updraft -REAL, DIMENSION(SIZE(PTHM,1),SIZE(PTHM,2)) :: ZWM_M ! w at mass levels - -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 - ZRC_MIX, ZRI_MIX ! guess of Rc and Ri for KF mixture -!================================================================================== -REAL, DIMENSION(SIZE(PTHM,1),SIZE(PTHM,2)) :: PTHVREF ! THv de r�f�rence -REAL, DIMENSION(SIZE(PTHM,1),SIZE(PTHM,2)) :: ZBUO ! Buoyancy -REAL, DIMENSION(SIZE(PTHM,1),SIZE(PTHM,2)) :: ZZDZ ! Dz -!================================================================================== - -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 - -REAL, DIMENSION(SIZE(PTHM,1)) :: ZLUP ! Upward Mixing length from the ground - -INTEGER :: ISV ! Number of scalar variables -INTEGER :: JK,JI,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, ZRSATW, ZRSATI -!================================================================================== -REAL, DIMENSION(SIZE(PTHM,1)) :: ZWUP_MEAN ! -REAL, DIMENSION(SIZE(PTHM,1)) :: ZCOE,ZWCOE,ZBUCOE -REAL, DIMENSION(SIZE(PTHM,1)) :: ZDETR_BUO, ZDETR_RT -!================================================================================== - -REAL :: ZDEPTH_MAX1, ZDEPTH_MAX2 ! control auto-extinction process -REAL :: XFRAC_LIM ! surface maximale du thermique - -REAL :: ZTMAX,ZRMAX, ZEPS ! control value - -REAL, DIMENSION(SIZE(PTHM,1),SIZE(PTHM,2)) :: ZSHEAR,ZDUDZ,ZDVDZ ! vertical wind shear -! pour le calcul de la resolution normalis�e -REAL, DIMENSION(SIZE(PTHM,1)) :: ZA1, ZRESOL_NORM -REAL :: ZRESOL_GRID -REAL, DIMENSION(SIZE(PTHM,1)) :: MODIF -! - -INTEGER :: KTCOUNT_MF ! current model time-step -REAL, DIMENSION (:), ALLOCATABLE :: ZWORK -! -!* 0.3 Declaration of namelists -! ------------------------ -! -!---------------------------------------------------------------------------- - -! 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 -XFRAC_LIM=0.5 -! -!------------------------------------------------------------------------ - -! INITIALISATION - -! Initialisation of the constants -ZRDORV = XRD / XRV !=0.622 -ZRVORD = (XRV / XRD) - -ZDEPTH_MAX1=3000. ! clouds with depth inferior to this value are keeped untouched -ZDEPTH_MAX2=4000. ! clouds with depth superior to this value are suppressed - -! Local variables, internal domain -! Internal Domain - -!number of scalar variables -ISV=SIZE(PSVM,3) - -IF (OENTR_DETR) THEN - ! si on prend en compte la r�solution - ! dans le calcul de l'entrainement et du d�trainement - !IF(LNORM_RESOL) THEN - !grid resolution in the AROME CASE - ZRESOL_GRID=sqrt(XDXHAT(1)*XDYHAT(1)) - !ENDIF - -! 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. - PRI_UP(:,:)=0. - PW_UP(:,:)=0. - ZTH_UP(:,:)=0. - PFRAC_UP(:,:)=0. - PTHV_UP(:,:)=0. - - PBUO_INTEG=0. - ZBUO =0. - - PFRAC_ICE_UP(:,:)=0. - PRSAT_UP(:,:)=PRVM(:,:) ! should be initialised correctly but is (normaly) not used - - !cloud/dry air mixture cloud content - ZRC_MIX = 0. - ZRI_MIX = 0. - -END IF - -! 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(:,:,JSV) = MZM_MF(KKA,KKU,KKL,PSVM(:,:,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(:,:,:)=ZSVM_F(:,:,:) -PSV_DO(:,:,:)=ZSVM_F(:,:,:) -PTHL_DO(:,:)=ZTHLM_F(:,:) -PRT_DO(:,:)=ZRTM_F(:,:) -PU_DO(:,:)=ZUM_F(:,:) -PV_DO(:,:)=ZVM_F(:,:) -PSV_DO(:,:,:)=0. - - ! initiation de l'�quation de la dynamique - !vertical velocity at mass level - ZWM_M(:,:)=MZF_MF(KKA,KKU,KKL,PWM(:,:)) -! Computation or initialisation of updraft characteristics at the KKB level -! thetal_up,rt_up,thetaV_up, w2,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)) -!------------------------ -! print*,OENTR_DETR -!------------------------ -IF (OENTR_DETR) THEN - 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(:,:) - - ZW_UP2(:,:)=ZEPS - ZW_UP2(:,KKB) = MAX(0.0001,(2./3.)*ZTKEM_F(:,KKB)) - - ! initiation de l'�quation de la dynamique - ! initialisation du vent de l'updraft pour la zone grise - ZW_UP(:,:)=SQRT(ZW_UP2) - - ! 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) - - ! Closure assumption for mass flux at KKB level - ! - ! calcul diff�rent de la flottabilit� - PTHVREF=300. - ! c'est la meilleure flottabilit� ! - !NBUOY=0. - !IF(NBUOY == 0.) THEN - ZG_O_THVREF=XG/ZTHVM_F - !ELSE - ! ! in AROME XTHVREF does not exist ici 300 - ! ZG_O_THVREF=XG/PTHVREF(1,1) ! on revient � l'�tat de r�f�rence et pas ZTHVM_F - !ENDIF - ! Calcul de la fermeture de Julien Pergaut comme limite max de PHY - - DO JK=KKB,KKE-KKL,KKL ! Vertical loop - ZZDZ(:,JK) = MAX(ZEPS,PZZ(:,JK+KKL)-PZZ(:,JK)) ! <== Delta Z between two flux level - ENDDO - - ! compute L_up - GLMIX=.TRUE. - ZTKEM_F(:,KKB)=0. - IF(CTURBLEN=='RM17') THEN - ZDUDZ = MZF_MF(KKA,KKU,KKL,GZ_M_W_MF(KKA,KKU,KKL,PUM,PDZZ)) - ZDVDZ = MZF_MF(KKA,KKU,KKL,GZ_M_W_MF(KKA,KKU,KKL,PVM,PDZZ)) - ZSHEAR = SQRT(ZDUDZ*ZDUDZ + ZDVDZ*ZDVDZ) - ELSE - ZSHEAR = 0. !no shear in bl89 mixing length - END IF - CALL COMPUTE_BL89_ML(KKA,KKB,KKE,KKU,KKL,PDZZ,ZTKEM_F(:,KKB),ZG_O_THVREF(:,KKB), & - ZTHVM_F,KKB,GLMIX,.TRUE.,ZSHEAR,ZLUP) - ZLUP(:)=MAX(ZLUP(:),1.E-10) - - ! Compute Buoyancy flux at the ground - ZWTHVSURF(:) = (ZTHVM_F(:,KKB)/ZTHM_F(:,KKB))*PSFTH(:)+ & - (0.61*ZTHM_F(:,KKB))*PSFRV(:) - - ! Mass flux at KKB level (updraft triggered if PSFTH>0.) - !elefant> - MODIF(:)=tanh(1.83*sqrt(XDXHAT(1)*XDYHAT(1))/ZLUP) - WHERE (ZWTHVSURF(:)>0.) - PEMF(:,KKB) = XCMF * MODIF(:) * 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. - ELSEWHERE - PEMF(:,KKB) =0. - GTEST(:)=.FALSE. - ENDWHERE - !elefant< -ELSE - GTEST(:)=PEMF(:,KKB+KKL)>0. -END IF - -!-------------------------------------------------------------------------- - -! 3. Vertical ascending loop -! ----------------------- -! -! If GTEST = T the updraft starts from the KKB level and stops when GTEST becomes F -! -! -GTESTLCL(:)=.FALSE. -GTESTETL(:)=.FALSE. - - ! CALCUL DE LA RESOLUTION NORMALISEE - - !IF(NLES_DTCOUNT .EQ. 0.)THEN - ! print*,"PROBLEME NLES_DTCOUNT" - ! print*,"NLES_DTCOUNT doit �tre 1." - !STOP - !ENDIF - - WHERE (GTEST .AND. ZLUP(:)>0. ) - ! hauteur des thermiques du pas de temps precedent - ! dans cette version ZLUP depend du flux de masse dans chaque - ! colonne, il serait peut-�tre plus simple de faire une variable - ! globale sur le domaine - ZRESOL_NORM(:)=ZRESOL_GRID/(ZLUP(:)) - ELSEWHERE - ZRESOL_NORM(:)=0.5 ! ARBITRAIRE - ENDWHERE - -! Loop on vertical level -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) print*,"cycle" - 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 - -! Buoyancy is computed on "flux" levels where updraft variables are known - !================================================================================ - ! CALCUL DE LA FLOTTABILITE - !================================================================================ - - ! Compute theta_v of updraft at flux level JK - - ZRC_UP(:) =PRC_UP(:,JK) ! guess - ZRI_UP(:) =PRI_UP(:,JK) ! guess - ZRV_UP(:) =PRV_UP(:,JK) - - CALL TH_R_FROM_THL_RT_1D(HFRAC_ICE,PFRAC_ICE_UP(:,JK),& - PPABSM(:,JK),PTHL_UP(:,JK),PRT_UP(:,JK),& - ZTH_UP(:,JK),ZRV_UP,ZRC_UP,ZRI_UP,ZRSATW(:),ZRSATI(:)) - - WHERE (GTEST) - PTHV_UP (:,JK) = ZTH_UP(:,JK)*(1.+ZRVORD*ZRV_UP(:))/(1.+PRT_UP(:,JK)) - ENDWHERE ! fin temporaire de GTEST - - ! test sur le calcul de la flottabilit� de 2 mani�res - ! le calcul de la flottabilit� est super important - ! c'est ce qui decide de Wup - !IF(NBUOY==1) THEN - ! !WHERE (GTEST) ! comme dans EDKF, mais c'est toujours positif - ! ! du coup le thermique ne s'arr�te pas - ! ZBUO (:,JK) = ZG_O_THVREF(:,JK)*(PTHV_UP(:,JK) - PTHVREF(:,JK)) - ! !ENDWHERE ! fin temporaire de GTEST - ! IF(LDUMMY1) THEN - ! print*,"ZBUO_1(:,",JK,")=",minval(ZBUO(:,JK)),maxval(ZBUO(:,JK)) - !END IF - !ELSEIF(NBUOY==2) THEN - ! !WHERE (GTEST)! comme ca devrait �tre mais trop restrictif, le thermique - ! !s'arrete d'embl�e - ! ZBUO (:,JK) = ZG_O_THVREF(:,JK)*(PTHV_UP(:,JK) -2.*ZTHVM_F(:,JK) + PTHVREF(:,JK)) - ! !ENDWHERE ! fin temporaire de GTEST - ! IF(LDUMMY1) THEN - ! print*,"ZBUO_2(:,",JK,")=",minval(ZBUO(:,JK)),maxval(ZBUO(:,JK)) - !END IF - !ELSEIF(NBUOY==3) THEN - ! !WHERE (GTEST) ! un peu plus rigoureux que celle du dessus, mais trop restrictif, le thermique - ! !s'arrete d'embl�e - ! ZBUO (:,JK) = XG/ZTHVM_F(:,JK)*(PTHV_UP(:,JK)-ZTHVM_F(:,JK))-& - ! XG/PTHVREF(:,JK)*(ZTHVM_F(:,JK)-PTHVREF(:,JK)) - ! !ENDWHERE ! fin temporaire de GTEST - !IF(LDUMMY1) THEN - ! print*,"ZBUO_3(:,",JK,")=",minval(ZBUO(:,JK)),maxval(ZBUO(:,JK)) - !END IF - !ELSEIF(NBUOY==4) THEN - ! devrait �tre plus proche de EDKF que NBUOY==0, mais ce n'est pas le cas. - ! IF(JK/=KKB) THEN - ! ZRC_MIX(:,JK) = ZRC_MIX(:,JK-KKL) ! guess of Rc of mixture - ! 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(:,:), & - ! PTHL_UP(:,JK),PRT_UP(:,JK),ZLUP(:), & - ! PRC_UP(:,JK),PRI_UP(:,JK),ZRC_MIX(:,JK),ZRI_MIX(:,JK), & - ! PENTR(:,JK),PDETR(:,JK),ZBUO(:,JK) ) - ! - - !IF(LDUMMY1) THEN - ! print*,"ZBUO_4(:,",JK,")=",minval(PBUO_INTEG(:,JK)),maxval(PBUO_INTEG(:,JK)) - !END IF - - !ELSE - !WHERE (GTEST) ! une derni�re possibilit� - ZBUO (:,JK) = ZG_O_THVREF(:,JK)*(PTHV_UP(:,JK) - ZTHVM_F(:,JK)) - !ENDWHERE ! fin temporaire de GTEST - !IF(LDUMMY1) THEN - ! print*,"ZBUO_0(:,",JK,")=",minval(ZBUO(:,JK)),maxval(ZBUO(:,JK)) - !END IF - !ENDIF - - !WHERE (GTEST) - ! anomalie de flottabilit� * DZ - PBUO_INTEG(:,JK) = ZBUO(:,JK)*(PZZ(:,JK+KKL)-PZZ(:,JK)) - ! uniquement pour le calcul de l'ancienne vitesse verticale - !ENDWHERE ! fin temporaire de GTEST - - !================================================================================ - ! CALCUL DE ZA1 - !================================================================================ - - ! si on a besoin d'un calcul sp�cifique de ZA1 vs XA1 (constante) - ! - LNORM_RESOL = calcul de ZA1 en prenant en compte la r�solution normalis�e - - ZA1(:)=XA1 - !IF (LNORM_RESOL) THEN - ! WHERE (GTEST) - ! ZA1(:)=XCOEF1+XCOEF2*ZRESOL_NORM(:)**(-1*XCOEF3) - ! ENDWHERE - !ELSEIF (.NOT. XDUMMY1 == 0. ) THEN - ! WHERE (GTEST) - ! ZA1(:)=XDUMMY1 - ! ENDWHERE - !END IF - ! - - !================================================================================ - ! EQUATION DE LA DYNAMIQUE - !================================================================================ - ! 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))) - ZWORK(:)=0. - !WHERE (GTEST .AND. PFRAC_UP(:,JK)>0 .AND. PFRAC_UP(:,JK)<=XFRAC_LIM) - WHERE (GTEST .AND. PFRAC_UP(:,JK)>0 ) - !hypothse => alpha(k)=alpha(k+1) - ZWCOE(:) = (1.-PFRAC_UP(:,JK))/(1.-PFRAC_UP(:,JK)+XBETA1) - ZBUCOE(:) = 2.*ZWCOE(:)*ZA1 - ZWORK(:)=(ZW_UP(:,JK)-(1.-ZWCOE(:))*PWM(:,JK)-ZWCOE(:)*ZWM_M(:,JK))*& - (ZW_UP(:,JK)-(1.-ZWCOE(:))*PWM(:,JK)-ZWCOE(:)*ZWM_M(:,JK))+& - ZBUCOE(:)*PBUO_INTEG(:,JK) - ELSEWHERE - ZWORK(:)=0 - ENDWHERE - ! - WHERE (GTEST) - WHERE(ZWORK>=0.) - ! il s'agit bien de Wu et non pas de Wu-Wm - ZW_UP(:,JK+KKL)=(1.-ZWCOE(:))*PWM(:,JK+KKL)+ZWCOE(:)*ZWM_M(:,JK)+SQRT(ZWORK(:)) - ELSEWHERE - ZW_UP(:,JK+KKL)=0. - ENDWHERE - ZW_UP2(:,JK+KKL) = MAX(0.,ZW_UP(:,JK+KKL)*ZW_UP(:,JK+KKL)) - ZWUP_MEAN(:) = MAX(ZEPS,0.5*(ZW_UP(:,JK+KKL)+ZW_UP(:,JK))-ZWM_M(:,JK)) - ENDWHERE - DEALLOCATE(ZWORK) - ! - !================================================================================ - ! CALCUL DE L ENTRAINEMENT ET DU DETRAINEMENT - !================================================================================ - ! Hypoth�se alpha jk= alpha point de masse - !WHERE (GTEST .AND. PFRAC_UP(:,JK)>0. .AND. PFRAC_UP(:,JK)<XFRAC_LIM ) - WHERE (GTEST .AND. PFRAC_UP(:,JK)>0. .AND. ZWUP_MEAN(:)>0 ) - ! ZWUP_MEAN est Wu-Wm au point de masse - PENTR(:,JK) = MAX(0., ((1.-PFRAC_UP(:,JK))*XBETA1)/(1.-PFRAC_UP(:,JK)+XBETA1)*(ZA1*ZBUO(:,JK)/& - (ZWUP_MEAN(:)*ZWUP_MEAN(:))& - - XB & - - 1./(ZWUP_MEAN(:))*(PWM(:,JK+KKL)-PWM(:,JK))/(PZZ(:,JK+KKL)-PZZ(:,JK)))& - ) - ZDETR_BUO(:) = MAX(0., -((1.-PFRAC_UP(:,JK))*XBETA1)/(1.-PFRAC_UP(:,JK)+XBETA1)*(ZA1*ZBUO(:,JK)/& - (ZWUP_MEAN(:)*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(:) - ENDWHERE -! - !================================================================================ - ! VARIABLES THERMODYNAMIQUES - !================================================================================ - ! computation of the updraft characteritics at jk+1 - ! WHERE (GTEST .AND. PFRAC_UP(:,JK)>0. .AND. PFRAC_UP(:,JK)<XFRAC_LIM) - WHERE (GTEST .AND. PFRAC_UP(:,JK)>0.) - ZMIX2(:) = (PZZ(:,JK+KKL)-PZZ(:,JK))*PENTR(:,JK)/(1.-PFRAC_UP(:,JK)) !& - ELSEWHERE - ZMIX2(:) = 0. !& - ENDWHERE ! GTEST - - -! If the updraft did not stop, compute cons updraft characteritics at jk+KKL - WHERE (GTEST) - 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 - - - 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 - END DO - -! 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 - 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) - PRC_UP(:,JK+KKL)=ZRC_UP(:) - PRV_UP(:,JK+KKL)=ZRV_UP(:) - PRI_UP(:,JK+KKL)=ZRI_UP(:) - 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+KKL - WHERE(GTEST) - PTHV_UP(:,JK+KKL) = ZTH_UP(:,JK+KKL)*((1+ZRVORD*PRV_UP(:,JK+KKL))/(1+PRT_UP(:,JK+KKL))) - ENDWHERE - !================================================================================ - ! CALCUL DU FLUX DE MASSE - !================================================================================ - ZMIX1(:)=0. - WHERE(GTEST) - ! identique dans Rio, EDKF ou dans la zone grise - ZMIX1(:)=ZZDZ(:,JK)*(PENTR(:,JK)-PDETR(:,JK)) - PEMF(:,JK+KKL)=PEMF(:,JK)*EXP(ZMIX1(:)) - ENDWHERE - !================================================================================ - ! FRACTION DE THERMIQUE - !================================================================================ - ! on cherche � savoir s'il y a des vitesses verticales non d�finies - ! je n'utilise que ZW_UP2 pour pouvoir avoir une valeur si ZW_UP - ! n'est pas d�fini -IF (maxval(ZW_UP2(:,JK+KKL)) .NE. maxval(ZW_UP2(:,JK+KKL))) & - call Print_msg( NVERB_FATAL, 'GEN', 'COMPUTE_UPDRAFT_HRIO', 'maxval(ZW_UP2(:,JK+KKL)) /= maxval(ZW_UP2(:,JK+KKL))' ) -! si on est dans la zone grise la d�finition du flux de masse change -! donc celle de alpha aussi -WHERE(GTEST) - !attention je ne suis pas au point de masse -WHERE ((SQRT(ZW_UP2(:,JK+KKL))-PWM(:,JK+KKL))>ZEPS) - PFRAC_UP(:,JK+KKL)=PEMF(:,JK+KKL)/((SQRT(ZW_UP2(:,JK+KKL))-PWM(:,JK+KKL))*ZRHO_F(:,JK+KKL)) -ELSEWHERE - PFRAC_UP(:,JK+KKL)=0. -ENDWHERE - PFRAC_UP(:,JK+KKL)=MIN(PFRAC_UP(:,JK+KKL),XFRAC_LIM) -ENDWHERE - - -! calcul des termes environmentaux au point de flux - WHERE(GTEST) - !WHERE(PFRAC_UP(:,JK+KKL)>0 .AND. PFRAC_UP(:,JK+KKL)< XFRAC_LIM) - WHERE( PFRAC_UP(:,JK+KKL)>0 ) - PTHL_DO(:,JK+KKL)=((PTHLM(:,JK)+PTHLM(:,JK+KKL))*0.5-PFRAC_UP(:,JK+KKL)*PTHL_UP(:,JK+KKL)) & - /(1.-PFRAC_UP(:,JK+KKL)) - PRT_DO(:,JK+KKL) =((PRTM(:,JK)+PRTM(:,JK+KKL))*0.5-PFRAC_UP(:,JK+KKL)*PRT_UP(:,JK+KKL)) & - /(1.-PFRAC_UP(:,JK+KKL)) - PU_DO(:,JK+KKL) = ((PUM(:,JK)+PUM(:,JK+KKL))*0.5-PFRAC_UP(:,JK+KKL)*PU_UP(:,JK+KKL)) & - /(1.-PFRAC_UP(:,JK+KKL)) - PV_DO(:,JK+KKL) = ((PVM(:,JK)+PVM(:,JK+KKL))*0.5-PFRAC_UP(:,JK+KKL)*PV_UP(:,JK+KKL)) & - /(1.-PFRAC_UP(:,JK+KKL)) - PTHV_DO(:,JK+KKL)=(ZTHVM_F(:,JK+KKL)-PFRAC_UP(:,JK+KKL)*PTHV_UP(:,JK+KKL))& - /(1.-PFRAC_UP(:,JK+KKL)) - ENDWHERE - ENDWHERE - DO JSV=1,ISV - IF (ONOMIXLG .AND. JSV >= KSV_LGBEG .AND. JSV<= KSV_LGEND) CYCLE - WHERE(GTEST) - !WHERE(PFRAC_UP(:,JK+KKL)>0 .AND. PFRAC_UP(:,JK+KKL)<= XFRAC_LIM) - WHERE(PFRAC_UP(:,JK+KKL)>0 ) - PSV_DO(:,JK+KKL,JSV) = ((PSVM(:,JK,JSV)+PSVM(:,JK+1,JSV))*0.5-PFRAC_UP(:,JK+KKL)*PSV_UP(:,JK+KKL,JSV))& - /(1.-PFRAC_UP(:,JK+KKL)) - ENDWHERE - ENDWHERE - ENDDO - - WHERE(GTEST) - PFRAC_UP(:,JK+KKL)=MIN(XFRAC_UP_MAX,PFRAC_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).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)) ) - ZW_UP2 (:,JK+KKL)=0. - PEMF (:,JK+KKL)=0. - 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)=ZRVM_F (:,JK+KKL) - PTHV_UP (:,JK+KKL)=ZTHVM_F(:,JK+KKL) - - PFRAC_UP (:,JK+KKL)=0. - KKCTL (:) =JK+KKL - PTHL_DO (:,JK+KKL)=ZTHLM_F(:,JK+KKL) - PRT_DO (:,JK+KKL)=ZRTM_F(:,JK+KKL) - PTHV_DO (:,JK+KKL)=ZTHVM_F(:,JK+KKL) - ENDWHERE - -ENDDO ! boucle JK -!STOP - 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 ZDEPTH_MAX1 m of depth) to 0 (for clouds of ZDEPTH_MAX2 m 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=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/read_exsegn.f90 b/src/MNH/read_exsegn.f90 index cd9d6ea125770025a819d9386f65d5fb2eba176b..9489c0986a3600a68939f1f2384a3e19edba72d7 100644 --- a/src/MNH/read_exsegn.f90 +++ b/src/MNH/read_exsegn.f90 @@ -298,6 +298,7 @@ END MODULE MODI_READ_EXSEG_n ! P. Wautelet 09/03/2021: simplify allocation of scalar variable names ! P. Wautelet 09/03/2021: move some chemistry initializations to ini_nsv ! P. Wautelet 10/03/2021: move scalar variable name initializations to ini_nsv +! R. Honnert 23/04/2021: add ADAP mixing length and delete HRIO and BOUT from CMF_UPDRAFT !------------------------------------------------------------------------------ ! !* 0. DECLARATIONS @@ -856,7 +857,7 @@ CALL TEST_NAM_VAR(ILUOUT,'CLBCY(1)',CLBCY(1),'CYCL','WALL','OPEN') CALL TEST_NAM_VAR(ILUOUT,'CLBCY(2)',CLBCY(2),'CYCL','WALL','OPEN') ! CALL TEST_NAM_VAR(ILUOUT,'CTURBDIM',CTURBDIM,'1DIM','3DIM') -CALL TEST_NAM_VAR(ILUOUT,'CTURBLEN',CTURBLEN,'DELT','BL89','RM17','DEAR','BLKR') +CALL TEST_NAM_VAR(ILUOUT,'CTURBLEN',CTURBLEN,'DELT','BL89','RM17','DEAR','BLKR','ADAP') CALL TEST_NAM_VAR(ILUOUT,'CTOM',CTOM,'NONE','TM06') CALL TEST_NAM_VAR(ILUOUT,'CSUBG_AUCV',CSUBG_AUCV,'NONE','CLFR','SIGM','PDF') ! @@ -877,8 +878,7 @@ CALL TEST_NAM_VAR(ILUOUT,'CTURBLEN_CLOUD',CTURBLEN_CLOUD,'NONE','DEAR','DELT','B ! ! The test on the mass flux scheme for shallow convection ! -CALL TEST_NAM_VAR(ILUOUT,'CMF_UPDRAFT',CMF_UPDRAFT,'NONE','EDKF','RHCJ',& - 'HRIO','BOUT') +CALL TEST_NAM_VAR(ILUOUT,'CMF_UPDRAFT',CMF_UPDRAFT,'NONE','EDKF','RHCJ') CALL TEST_NAM_VAR(ILUOUT,'CMF_CLOUD',CMF_CLOUD,'NONE','STAT','DIRE') ! ! The test on the CSOLVER name is made elsewhere @@ -1566,7 +1566,7 @@ ELSE END IF END IF ! -IF(CTURBLEN=='RM17') THEN +IF(CTURBLEN=='RM17' .OR. CTURBLEN=='ADAP') THEN XCEDIS=0.34 ELSE XCEDIS=0.84 diff --git a/src/MNH/shallow_mf.f90 b/src/MNH/shallow_mf.f90 index 4017b49d7540651b21d3d4040afd68f57418de68..26b53e19d7778582162b2f6d076be2f350d2e890 100644 --- a/src/MNH/shallow_mf.f90 +++ b/src/MNH/shallow_mf.f90 @@ -169,6 +169,7 @@ END MODULE MODI_SHALLOW_MF !! Q.Rodier 01/2019 : support RM17 mixing length !! R.Honnert 1/2019 : remove SURF ! P. Wautelet 10/04/2019: replace ABORT and STOP calls by Print_msg +! R.Honnert 04/2021: remove HRIO and BOUT schemes !! -------------------------------------------------------------------------- ! !* 0. DECLARATIONS @@ -183,7 +184,6 @@ 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 @@ -308,8 +308,7 @@ 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 == 'RHCJ') THEN PENTR = 1.E20 PDETR = 1.E20 PEMF = 1.E20 @@ -337,7 +336,7 @@ 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') THEN GENTR_DETR = .TRUE. CALL COMPUTE_UPDRAFT(KKA,IKB,IKE,KKU,KKL,HFRAC_ICE,GENTR_DETR,OMIXUV,& ONOMIXLG,KSV_LGBEG,KSV_LGEND, & @@ -378,22 +377,6 @@ ELSEIF (HMF_UPDRAFT == 'RAHA') THEN ZDEPTH ) ELSEIF (HMF_UPDRAFT == 'DUAL') THEN !Updraft characteristics are already computed and received by interface -ELSEIF (HMF_UPDRAFT == 'HRIO') THEN - GENTR_DETR = .TRUE. - ! ma version avec l'entrainement de Rio et al. - CALL COMPUTE_UPDRAFT_HRIO(KKA,IKB,IKE,KKU,KKL,HFRAC_ICE,GENTR_DETR,OMIXUV, & - ONOMIXLG,KSV_LGBEG,KSV_LGEND, & - PZZ,PDZZ, & - PSFTH,PSFRV,PPABSM,PRHODREF, & - PUM,PVM,PTKEM,PWM, & - 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, & - PTHL_DO, PTHV_DO, PRT_DO, & - PU_DO, PV_DO, ZSV_DO, & - PEMF,PDETR, & - PENTR,ZBUO_INTEG,KKLCL,KKETL,KKCTL,ZDEPTH ) ELSE call Print_msg( NVERB_FATAL, 'GEN', 'SHALLOW_MF', 'no updraft model for EDKF: CMF_UPDRAFT='//trim(HMF_UPDRAFT) ) ENDIF @@ -417,7 +400,7 @@ CALL COMPUTE_MF_CLOUD(KKA,IKB,IKE,KKU,KKL,KRR,KRRL,KRRI,& !!! ------------------------------------------------------------------------ ! ZEMF_O_RHODREF=PEMF/PRHODREF - IF(HMF_UPDRAFT == 'EDKF' .OR. HMF_UPDRAFT == 'RHCJ'.OR. HMF_UPDRAFT == 'BOUT') THEN +IF(HMF_UPDRAFT == 'EDKF' .OR. HMF_UPDRAFT == 'RHCJ') THEN IF ( PIMPL_MF > 1.E-10 ) THEN CALL MF_TURB(KKA, IKB, IKE, KKU, KKL, OMIXUV, & ONOMIXLG,KSV_LGBEG,KSV_LGEND, & @@ -429,66 +412,17 @@ 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, & + 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, & - PIMPL_MF, PTSTEP, & - PDZZ, & - PRHODJ, & - ZTHLM,ZTHVM,ZRTM,PUM,PVM,PSVM, & - PDTHLDT_MF,PDRTDT_MF,PDUDT_MF,PDVDT_MF,PDSVDT_MF, & - ZEMF_O_RHODREF,PTHL_UP,PTHV_UP,PRT_UP,PU_UP,PV_UP,ZSV_UP,& - PTHL_DO,PTHV_DO,PRT_DO,PU_DO,PV_DO,ZSV_DO, & - PFLXZTHMF,PFLXZTHVMF,PFLXZRMF,PFLXZUMF,PFLXZVMF, & - ZFLXZSVMF ) - ELSE - call Print_msg( NVERB_FATAL, 'GEN', 'SHALLOW_MF', 'no updraft model for EDKF: CMF_UPDRAFT='//trim(HMF_UPDRAFT) ) - END IF - - IF (HMF_UPDRAFT == 'BOUT') THEN - !! calcul de la hauteur de la couche limite ou de L_up - DO JK=1,IKE-KKL - PTHVREF(:,JK)=RESHAPE(XTHVREF(:,:,JK),(/SIZE(PTHM,1)*SIZE(PTHM,2)/) ) - ENDDO - ZG_O_THVREF=XG/PTHVREF - GLMIX=.TRUE. - IF(CTURBLEN=='RM17') THEN - ZDUDZ = MZF_MF(KKA,KKU,KKL,GZ_M_W_MF(KKA,KKU,KKL,PUM,PDZZ)) - ZDVDZ = MZF_MF(KKA,KKU,KKL,GZ_M_W_MF(KKA,KKU,KKL,PVM,PDZZ)) - ZSHEAR = SQRT(ZDUDZ*ZDUDZ + ZDVDZ*ZDVDZ) - ELSE - ZSHEAR = 0. !no shear in bl89 mixing length - END IF - CALL COMPUTE_BL89_ML(KKA,IKB,IKE,KKU,KKL,PDZZ,PTKEM(:,IKB) ,ZG_O_THVREF(:,IKB),ZTHVM,IKB,GLMIX,.TRUE.,ZSHEAR,ZLUP) - !! calcul de Dx/(h+hc) - DO JI=1,SIZE(XDXHAT) - DO JJ=1,SIZE(XDYHAT) - ZRESOL_GRID((JJ-1)*SIZE(XDXHAT)+JI)=SQRT(XDXHAT(JI)*XDYHAT(JJ)) - ENDDO - ENDDO - ZRESOL_NORM=ZRESOL_GRID/ZLUP - !! P=loi pour MF, on utilise la même loi à chaque fois - ZPLAW=(ZRESOL_NORM*ZRESOL_NORM+0.19*ZRESOL_NORM**(2./3.))/ & - (ZRESOL_NORM*ZRESOL_NORM+0.15*ZRESOL_NORM**(2./3.)+0.33) - !! reduction des flux a posteriori - !! MF=P*MF en première approximation, on oublie w'f' (Kgrad) et w'f'resol (nul avec ce flux) - ! - DO JK=1,IKE-KKL - PFLXZTHMF(:,JK)=PFLXZTHMF(:,JK)*ZPLAW - PFLXZTHVMF(:,JK)=PFLXZTHVMF(:,JK)*ZPLAW - PFLXZRMF(:,JK)=PFLXZRMF(:,JK)*ZPLAW - PFLXZUMF(:,JK)=PFLXZUMF(:,JK)*ZPLAW - PFLXZVMF(:,JK)=PFLXZVMF(:,JK)*ZPLAW - ENDDO - END IF + ENDIF +ELSE + call Print_msg( NVERB_FATAL, 'GEN', 'SHALLOW_MF', 'no updraft model for EDKF: CMF_UPDRAFT='//trim(HMF_UPDRAFT) ) +END IF ! security in the case HMF_UPDRAFT = 'DUAL' ! to be modified if 'DUAL' is evolving (momentum mixing for example) diff --git a/src/MNH/turb.f90 b/src/MNH/turb.f90 index 8969a939e0c888375df7f9eaa0773d7fa229f7cc..b46c4559326f0a23a007b7692440ab071fc78314 100644 --- a/src/MNH/turb.f90 +++ b/src/MNH/turb.f90 @@ -344,6 +344,7 @@ END MODULE MODI_TURB ! P. Wautelet 11/06/2020: bugfix: correct PRSVS array indices ! P. Wautelet + Benoit Vié 06/2020: improve removal of negative scalar variables + adapt the corresponding budgets ! P. Wautelet 30/06/2020: move removal of negative scalar variables to Sources_neg_correct +!! 02/2021 (R. Honnert/V. Masson) New mixing length in the grey zone ! -------------------------------------------------------------------------- ! !* 0. DECLARATIONS @@ -502,6 +503,7 @@ REAL, ALLOCATABLE, DIMENSION(:,:,:) ::& ZEXN, & ! EXN at t-1 ZT, & ! T at t-1 ZLOCPEXNM, & ! Lv/Cp/EXNREF at t-1 + ZLMW, & ! Turbulent mixing length (work array) ZLEPS, & ! Dissipative length ZTRH, & ! Dynamic and Thermal Production of TKE ZATHETA,ZAMOIST, & ! coefficients for s = f (Thetal,Rnp) @@ -541,8 +543,10 @@ INTEGER :: IKTB,IKTE ! start, end of k loops in physical domain INTEGER :: JRR,JK,JSV ! loop counters INTEGER :: JI,JJ ! loop counters REAL :: ZL0 ! Max. Mixing Length in Blakadar formula -REAL :: ZALPHA ! proportionnality constant between Dz/2 and -! ! BL89 mixing length near the surface +REAL :: ZALPHA ! work coefficient : + ! - proportionnality constant between Dz/2 and +! ! BL89 mixing length near the surface + ! - and coefficient to reduce DELT in ADAP ! REAL :: ZTIME1, ZTIME2 REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3)):: ZTT,ZEXNE,ZLV,ZLS,ZCPH,ZCOR @@ -555,6 +559,7 @@ ALLOCATE ( & ZEXN(SIZE(PTHLT,1),SIZE(PTHLT,2),SIZE(PTHLT,3)), & ZT(SIZE(PTHLT,1),SIZE(PTHLT,2),SIZE(PTHLT,3)), & ZLOCPEXNM(SIZE(PTHLT,1),SIZE(PTHLT,2),SIZE(PTHLT,3)), & + ZLMW(SIZE(PTHLT,1),SIZE(PTHLT,2),SIZE(PTHLT,3)), & ZLEPS(SIZE(PTHLT,1),SIZE(PTHLT,2),SIZE(PTHLT,3)), & ZTRH(SIZE(PTHLT,1),SIZE(PTHLT,2),SIZE(PTHLT,3)), & ZATHETA(SIZE(PTHLT,1),SIZE(PTHLT,2),SIZE(PTHLT,3)), & @@ -759,19 +764,37 @@ SELECT CASE (HTURBLEN) ZSHEAR = SQRT(ZDUDZ*ZDUDZ + ZDVDZ*ZDVDZ) CALL BL89(KKA,KKU,KKL,PZZ,PDZZ,PTHVREF,ZTHLM,KRR,ZRM,PTKET,ZSHEAR,PLEM) ! -!* 3.3 Delta mixing length +!* 3.3 Grey-zone combined RM17 & Deardorff mixing lengths +! -------------------------------------------------- + + CASE ('ADAP') + ZDUDZ = MXF(MZF(GZ_U_UW(PUT,PDZZ))) + ZDVDZ = MYF(MZF(GZ_V_VW(PVT,PDZZ))) + ZSHEAR = SQRT(ZDUDZ*ZDUDZ + ZDVDZ*ZDVDZ) + CALL BL89(KKA,KKU,KKL,PZZ,PDZZ,PTHVREF,ZTHLM,KRR,ZRM,PTKET,ZSHEAR,PLEM) + + CALL DELT(ZLMW,ODZ=.FALSE.) + ! The minimum mixing length is chosen between Horizontal grid mesh (not taking into account the vertical grid mesh) and RM17. + ! For large horizontal grid meshes, this is equal to RM17 + ! For LES grid meshes, this is equivalent to Deardorff : the base mixing lentgh is the horizontal grid mesh, + ! and it is limited by a stability-based length (RM17), as was done in Deardorff length (but taking into account shear as well) + ! For grid meshes in the grey zone, then this is the smaller of the two. + ZALPHA=0.50 + PLEM = MIN(PLEM,ZALPHA*ZLMW) +! +!* 3.4 Delta mixing length ! ------------------- ! CASE ('DELT') - CALL DELT(PLEM) + CALL DELT(PLEM,ODZ=.TRUE.) ! -!* 3.4 Deardorff mixing length +!* 3.5 Deardorff mixing length ! ----------------------- ! CASE ('DEAR') CALL DEAR(PLEM) ! -!* 3.5 Blackadar mixing length +!* 3.6 Blackadar mixing length ! ----------------------- ! CASE ('BLKR') @@ -1411,7 +1434,7 @@ REAL, DIMENSION(SIZE(PEXN,1),SIZE(PEXN,2),SIZE(PEXN,3)) :: ZDRVSATDT END SUBROUTINE COMPUTE_FUNCTION_THERMO ! ! #################### - SUBROUTINE DELT(PLM) + SUBROUTINE DELT(PLM,ODZ) ! #################### !! !!**** *DELT* routine to compute mixing length for DELT case @@ -1433,6 +1456,7 @@ END SUBROUTINE COMPUTE_FUNCTION_THERMO !* 0.1 Declarations of dummy arguments ! REAL, DIMENSION(:,:,:), INTENT(OUT) :: PLM +LOGICAL, OPTIONAL, INTENT(IN) :: ODZ ! !* 0.2 Declarations of local variables ! @@ -1440,18 +1464,32 @@ REAL :: ZD ! distance to the surface ! !------------------------------------------------------------------------------- ! -DO JK = IKTB,IKTE ! 1D turbulence scheme - PLM(:,:,JK) = PZZ(:,:,JK+KKL) - PZZ(:,:,JK) -END DO -PLM(:,:,KKU) = PLM(:,:,IKE) -PLM(:,:,KKA) = PZZ(:,:,IKB) - PZZ(:,:,KKA) -IF ( HTURBDIM /= '1DIM' ) THEN ! 3D turbulence scheme - IF ( L2D) THEN - PLM(:,:,:) = SQRT( PLM(:,:,:)*MXF(PDXX(:,:,:)) ) - ELSE - PLM(:,:,:) = (PLM(:,:,:)*MXF(PDXX(:,:,:))*MYF(PDYY(:,:,:)) ) ** (1./3.) +IF (ODZ) THEN + ! Dz is take into account in the computation + DO JK = IKTB,IKTE ! 1D turbulence scheme + PLM(:,:,JK) = PZZ(:,:,JK+KKL) - PZZ(:,:,JK) + END DO + PLM(:,:,KKU) = PLM(:,:,IKE) + PLM(:,:,KKA) = PZZ(:,:,IKB) - PZZ(:,:,KKA) + IF ( HTURBDIM /= '1DIM' ) THEN ! 3D turbulence scheme + IF ( L2D) THEN + PLM(:,:,:) = SQRT( PLM(:,:,:)*MXF(PDXX(:,:,:)) ) + ELSE + PLM(:,:,:) = (PLM(:,:,:)*MXF(PDXX(:,:,:))*MYF(PDYY(:,:,:)) ) ** (1./3.) + END IF + END IF +ELSE + ! Dz not taken into account in computation to assure invariability with vertical grid mesh + PLM=1.E10 + IF ( HTURBDIM /= '1DIM' ) THEN ! 3D turbulence scheme + IF ( L2D) THEN + PLM(:,:,:) = MXF(PDXX(:,:,:)) + ELSE + PLM(:,:,:) = (MXF(PDXX(:,:,:))*MYF(PDYY(:,:,:)) ) ** (1./2.) + END IF END IF END IF + ! ! mixing length limited by the distance normal to the surface ! (with the same factor as for BL89) @@ -1483,7 +1521,7 @@ END SUBROUTINE DELT SUBROUTINE DEAR(PLM) ! #################### !! -!!**** *DELT* routine to compute mixing length for DEARdorff case +!!**** *DEAR* routine to compute mixing length for DEARdorff case ! !! AUTHOR !! ------ @@ -1708,14 +1746,14 @@ ELSE ! !* 3.1 BL89 mixing length ! ------------------ - CASE ('BL89','RM17') + CASE ('BL89','RM17','ADAP') ZSHEAR=0. CALL BL89(KKA,KKU,KKL,PZZ,PDZZ,PTHVREF,ZTHLM,KRR,ZRM,PTKET,ZSHEAR,ZLM_CLOUD) ! !* 3.2 Delta mixing length ! ------------------- CASE ('DELT') - CALL DELT(ZLM_CLOUD) + CALL DELT(ZLM_CLOUD,ODZ=.TRUE.) ! !* 3.3 Deardorff mixing length ! -----------------------