diff --git a/aux/ini_phyex.F90 b/aux/ini_phyex.F90 index 79c583f0..e0e4e4c0 100644 --- a/aux/ini_phyex.F90 +++ b/aux/ini_phyex.F90 @@ -1,7 +1,8 @@ SUBROUTINE INI_PHYEX(HPROGRAM, KUNITNML, LDNEEDNAM, KLUOUT, KFROM, KTO, & &PTSTEP, PDZMIN, & &CMICRO, CSCONV, CTURB, & - &LDCHANGEMODEL, LDDEFAULTVAL, LDREADNAM, LDCHECK, KPRINT, LDINIT, & + &LDCHANGEMODEL, LDDEFAULTVAL, LDREADNAM, LDCHECK, & + &KPRINT, LDINIT, & &PHYEX_IN, PHYEX_OUT) ! USE MODD_PHYEX, ONLY: PHYEX_t @@ -228,7 +229,7 @@ SUBROUTINE INI_PHYEX(HPROGRAM, KUNITNML, LDNEEDNAM, KLUOUT, KFROM, KTO, & ENDIF CALL TURBN_INIT(HPROGRAM, KUNITNML, LDNEEDNAM, KLUOUT, & - &LDDEFAULTVAL, LDREADNAM, LDCHECK, KPRINT) + &LDDEFAULTVAL=LDDEFAULTVAL, LDREADNAM=LDREADNAM, LDCHECK=LDCHECK,KPRINT=KPRINT) IF(LLINIT) THEN CALL INI_TURB(HPROGRAM) ENDIF diff --git a/aux/modd_nsv.F90 b/aux/modd_nsv.F90 index 958d1f4d..7e6afa10 100644 --- a/aux/modd_nsv.F90 +++ b/aux/modd_nsv.F90 @@ -52,7 +52,7 @@ MODULE MODD_NSV LOGICAL :: LINI_NSV(JPMODELMAX) = .FALSE. ! becomes True when routine INI_NSV is called ! CHARACTER(LEN=NMNHNAMELGTMAX), DIMENSION(:,:), ALLOCATABLE :: CSV_CHEM_LIST_A !Names of all the chemical variables -CHARACTER(LEN=6), DIMENSION(:,:), ALLOCATABLE :: CSV_A !Names of the scalar variables +CHARACTER(LEN=16), DIMENSION(:,:), ALLOCATABLE :: CSV_A !Names of the scalar variables TYPE(tfieldmetadata), DIMENSION(:,:), ALLOCATABLE :: TSVLIST_A !Metadata of all the scalar variables INTEGER,DIMENSION(JPMODELMAX)::NSV_A = 0 ! total number of scalar variables @@ -167,7 +167,7 @@ MODULE MODD_NSV ! variables updated for the current model ! CHARACTER(LEN=NMNHNAMELGTMAX), DIMENSION(:), POINTER :: CSV_CHEM_LIST !Names of all the chemical variables -CHARACTER(LEN=6), DIMENSION(:), POINTER :: CSV !Names of the scalar variables +CHARACTER(LEN=16), DIMENSION(:), POINTER :: CSV !Names of the scalar variables TYPE(tfieldmetadata), DIMENSION(:), POINTER :: TSVLIST !Metadata of all the scalar variables @@ -289,7 +289,7 @@ MODULE MODD_NSV LOGICAL, POINTER :: LINI_NSV(:) => NULL() ! CHARACTER(LEN=NMNHNAMELGTMAX), DIMENSION(:,:), POINTER :: CSV_CHEM_LIST_A => NULL() -CHARACTER(LEN=6), DIMENSION(:,:), POINTER :: CSV_A => NULL() +CHARACTER(LEN=16), DIMENSION(:,:), POINTER :: CSV_A => NULL() TYPE(tfieldmetadata), DIMENSION(:,:), POINTER :: TSVLIST_A => NULL() INTEGER, DIMENSION(:), POINTER ::NSV_A => NULL(), & @@ -376,7 +376,7 @@ MODULE MODD_NSV NSV_SNWEND_A => NULL() CHARACTER(LEN=NMNHNAMELGTMAX), DIMENSION(:), POINTER :: CSV_CHEM_LIST => NULL() -CHARACTER(LEN=6), DIMENSION(:), POINTER :: CSV => NULL() +CHARACTER(LEN=16), DIMENSION(:), POINTER :: CSV => NULL() TYPE(tfieldmetadata), DIMENSION(:), POINTER :: TSVLIST => NULL() diff --git a/aux/modi_ini_phyex.F90 b/aux/modi_ini_phyex.F90 index 70e2c41f..5ddb5a18 100644 --- a/aux/modi_ini_phyex.F90 +++ b/aux/modi_ini_phyex.F90 @@ -4,7 +4,8 @@ MODULE MODI_INI_PHYEX SUBROUTINE INI_PHYEX(HPROGRAM, KUNITNML, LDNEEDNAM, KLUOUT, KFROM, KTO, & &PTSTEP, PDZMIN, & &CMICRO, CSCONV, CTURB, & - &LDCHANGEMODEL, LDDEFAULTVAL, LDREADNAM, LDCHECK, KPRINT, LDINIT, & + &LDCHANGEMODEL, LDDEFAULTVAL, LDREADNAM, LDCHECK,& + &KPRINT, LDINIT, & &PHYEX_IN, PHYEX_OUT) ! USE MODD_PHYEX, ONLY: PHYEX_t diff --git a/conv/convect_chem_transport.F90 b/conv/convect_chem_transport.F90 index b2386853..80060d95 100644 --- a/conv/convect_chem_transport.F90 +++ b/conv/convect_chem_transport.F90 @@ -95,6 +95,7 @@ SUBROUTINE CONVECT_CHEM_TRANSPORT( CVPEXT, D, NSV, KCH, PCH1, PCH1C, & INTEGER :: JI ! horizontal loop index INTEGER :: JK, JKP ! vertical loop index INTEGER :: JN ! chemical tracer loop index +INTEGER :: JCH INTEGER :: JSTEP ! fractional time loop index INTEGER :: JKLD, JKLP, JKMIN, JKMAX, JKMAX2 ! loop index for levels ! @@ -212,7 +213,9 @@ SUBROUTINE CONVECT_CHEM_TRANSPORT( CVPEXT, D, NSV, KCH, PCH1, PCH1C, & !* 4. Final closure (environmental) computations ! ------------------------------------------ ! -PCH1C(D%NIB:D%NIE,IKB:IKE,1:KCH) = PCH1(D%NIB:D%NIE,IKB:IKE,1:KCH) ! initialize adjusted envir. values +DO JCH = 1, KCH + PCH1C(:,IKB:IKE,JCH) = PCH1(:,IKB:IKE,JCH) ! initialize adjusted envir. values +ENDDO ! DO JK = IKB, IKE DO JI=D%NIB,D%NIE @@ -233,8 +236,10 @@ SUBROUTINE CONVECT_CHEM_TRANSPORT( CVPEXT, D, NSV, KCH, PCH1, PCH1C, & ENDDO ENDDO ! -ZCH1MFIN(D%NIB:D%NIE,1:D%NKT,1:KCH) = 0. -ZCH1MFOUT(D%NIB:D%NIE,1:D%NKT,1:KCH) = 0. +DO JCH = 1, KCH + ZCH1MFIN(:,1:D%NKT,JCH) = 0. + ZCH1MFOUT(:,1:D%NKT,JCH) = 0. +ENDDO ! DO JSTEP = 1, KFTSTEPS ! Enter the fractional time step loop ! diff --git a/conv/convect_closure.F90 b/conv/convect_closure.F90 index 7912c2bf..729190f6 100644 --- a/conv/convect_closure.F90 +++ b/conv/convect_closure.F90 @@ -146,6 +146,7 @@ SUBROUTINE CONVECT_CLOSURE( KLON, KLEV, & REAL, DIMENSION(KLON,KLEV), INTENT(INOUT):: PDTEVRF! downdraft evaporation rate REAL, DIMENSION(KLON,KLEV), INTENT(OUT) :: PPRLFLX! liquid precip flux REAL, DIMENSION(KLON,KLEV), INTENT(OUT) :: PPRSFLX! solid precip flux + ! !* 0.2 Declarations of local variables : ! @@ -156,6 +157,7 @@ SUBROUTINE CONVECT_CLOSURE( KLON, KLEV, & INTEGER :: JITER ! iteration loop index INTEGER :: JSTEP ! fractional time loop index REAL :: ZCPORD, ZRDOCP ! C_pd / R_d, R_d / C_pd +REAL :: ZEPS ! REAL, DIMENSION(KLON,KLEV) :: ZTHLC ! convectively adjusted ! grid scale enthalpy @@ -203,14 +205,15 @@ SUBROUTINE CONVECT_CLOSURE( KLON, KLEV, & LOGICAL, DIMENSION(KLON) :: GWORK1, GWORK3! work arrays LOGICAL, DIMENSION(KLON,KLEV) :: GWORK4 ! work array ! +REAL(KIND=JPHOOK) :: ZHOOK_HANDLE ! +#include "convect_closure_thrvlcl.h" !------------------------------------------------------------------------------- ! !* 0.2 Initialize local variables ! ---------------------------- ! ! -REAL(KIND=JPHOOK) :: ZHOOK_HANDLE IF (LHOOK) CALL DR_HOOK('CONVECT_CLOSURE',0,ZHOOK_HANDLE) PSPR(:) = 0. ZTIMC(:,:) = 0. @@ -225,6 +228,7 @@ SUBROUTINE CONVECT_CLOSURE( KLON, KLEV, & GWORK4(:,:) = .FALSE. ILCL(:) = KLCL(:) ! +ZEPS = XRD / XRV ZCPORD = XCPD / XRD ZRDOCP = XRD / XCPD ! @@ -545,10 +549,11 @@ SUBROUTINE CONVECT_CLOSURE( KLON, KLEV, & ! that in routine TRIGGER_FUNCT ! --------------------------------------------- ! - CALL CONVECT_CLOSURE_THRVLCL( KLON, KLEV, & - PPRES, PTHC, PRWC, PZ, GWORK1, & - ZTHLCL, ZRVLCL, ZZLCL, ZTLCL, ZTELCL, & - ILCL, KDPL, KPBL ) +CALL ABOR1('FIXME: THE INTERFACE IS WRONG') + !CALL CONVECT_CLOSURE_THRVLCL( KLON, KLEV, & + !PPRES, PTHC, PRWC, PZ, GWORK1, & + !ZTHLCL, ZRVLCL, ZZLCL, ZTLCL, ZTELCL, & + !ILCL, KDPL, KPBL ) ! ! ZTLCL(:) = MAX( 230., MIN( 335., ZTLCL(:) ) ) ! set some overflow bounds @@ -565,7 +570,9 @@ SUBROUTINE CONVECT_CLOSURE( KLON, KLEV, & ZPI(:) = MAX( 0.95, MIN( 1.5, ZPI(:) ) ) ZWORK1(:) = XP00 / ZPI(:) ** ZCPORD ! pressure at LCL ! - CALL CONVECT_SATMIXRATIO( KLON, ZWORK1, ZTELCL, ZWORK3, ZLV, ZLS, ZCPH ) + DO JI = 1, IIE + CALL CONVECT_SATMIXRATIO( ZWORK1(JI), ZTELCL(JI), ZEPS, ZWORK3(JI), ZLV(JI), ZLS(JI), ZCPH(JI) ) + END DO ZWORK3(:) = MIN( .1, MAX( 0., ZWORK3(:) ) ) ! ! compute theta_e updraft undilute @@ -593,7 +600,9 @@ SUBROUTINE CONVECT_CLOSURE( KLON, KLEV, & ZWORK2(JI) = PTHC(JI,JK) / ZPI(JI) END DO ! - CALL CONVECT_SATMIXRATIO( KLON, PPRES(:,JK), ZWORK2, ZWORK3, ZLV, ZLS, ZCPH ) + DO JI = 1, IIE + CALL CONVECT_SATMIXRATIO( PPRES(JI,JK), ZWORK2(JI), ZEPS, ZWORK3(JI), ZLV(JI), ZLS(JI), ZCPH(JI) ) + END DO ! ! DO JI = 1, IIE @@ -659,4 +668,7 @@ SUBROUTINE CONVECT_CLOSURE( KLON, KLEV, & ! ! IF (LHOOK) CALL DR_HOOK('CONVECT_CLOSURE',1,ZHOOK_HANDLE) +CONTAINS +INCLUDE "convect_satmixratio.h" +! END SUBROUTINE CONVECT_CLOSURE diff --git a/conv/convect_closure_adjust_shal.F90 b/conv/convect_closure_adjust_shal.F90 index 04fc6d91..23a3506d 100644 --- a/conv/convect_closure_adjust_shal.F90 +++ b/conv/convect_closure_adjust_shal.F90 @@ -93,13 +93,13 @@ SUBROUTINE CONVECT_CLOSURE_ADJUST_SHAL( CVPEXT, D, PADJ, & ! specified degree of stabilization ! ---------------------------------------------------- ! - DO JK = IKB + 1, IKE - DO JI = D%NIB, D%NIE - PUMF(JI,JK) = PZUMF(JI,JK) * PADJ(JI) - PUER(JI,JK) = PZUER(JI,JK) * PADJ(JI) - PUDR(JI,JK) = PZUDR(JI,JK) * PADJ(JI) - ENDDO - END DO +DO JK = IKB + 1, IKE + DO JI = D%NIB, D%NIE + PUMF(JI,JK) = PZUMF(JI,JK) * PADJ(JI) + PUER(JI,JK) = PZUER(JI,JK) * PADJ(JI) + PUDR(JI,JK) = PZUDR(JI,JK) * PADJ(JI) + ENDDO +END DO ! IF (LHOOK) CALL DR_HOOK('CONVECT_CLOSURE_ADJUST_SHAL',1,ZHOOK_HANDLE) END SUBROUTINE CONVECT_CLOSURE_ADJUST_SHAL diff --git a/conv/convect_closure_shal.F90 b/conv/convect_closure_shal.F90 index 2416a57d..c3ab3db1 100644 --- a/conv/convect_closure_shal.F90 +++ b/conv/convect_closure_shal.F90 @@ -85,48 +85,47 @@ SUBROUTINE CONVECT_CLOSURE_SHAL( CVP_SHAL, CVPEXT, CST, D, & ! !* 0.1 Declarations of dummy arguments : ! -TYPE(CONVPAR_SHAL) ,INTENT(IN) :: CVP_SHAL -TYPE(CONVPAREXT) ,INTENT(IN) :: CVPEXT -TYPE(CST_T) ,INTENT(IN) :: CST -TYPE(DIMPHYEX_T) ,INTENT(IN) :: D -REAL ,DIMENSION(D%NIT,D%NKT) ,INTENT(IN) :: PPRES ! pressure (P) -REAL ,DIMENSION(D%NIT,D%NKT) ,INTENT(IN) :: PDPRES ! pressure difference between - ! bottom and top of layer (Pa) -REAL ,DIMENSION(D%NIT,D%NKT) ,INTENT(IN) :: PZ ! height of model layer (m) -REAL ,DIMENSION(D%NIT,D%NKT) ,INTENT(IN) :: PLMASS ! mass of model layer (kg) -REAL ,DIMENSION(D%NIT,D%NKT) ,INTENT(IN) :: PTHL ! grid scale enthalpy (J/kg) -REAL ,DIMENSION(D%NIT,D%NKT) ,INTENT(IN) :: PTH ! grid scale theta -REAL ,DIMENSION(D%NIT,D%NKT) ,INTENT(IN) :: PRW ! grid scale total water +TYPE(CONVPAR_SHAL), INTENT(IN) :: CVP_SHAL +TYPE(CONVPAREXT), INTENT(IN) :: CVPEXT +TYPE(CST_T), INTENT(IN) :: CST +TYPE(DIMPHYEX_T), INTENT(IN) :: D +INTEGER, DIMENSION(D%NIT), INTENT(IN) :: KLCL ! index lifting condens. level +INTEGER, DIMENSION(D%NIT), INTENT(IN) :: KCTL ! index for cloud top level +INTEGER, DIMENSION(D%NIT), INTENT(IN) :: KDPL ! index for departure level +INTEGER, DIMENSION(D%NIT), INTENT(IN) :: KPBL ! index for top of source layer +REAL, DIMENSION(D%NIT), INTENT(INOUT) :: PTIMEC ! convection time step +REAL, DIMENSION(D%NIT,D%NKT),INTENT(IN) :: PTHL ! grid scale enthalpy (J/kg) +REAL, DIMENSION(D%NIT,D%NKT),INTENT(IN) :: PTH ! grid scale theta +REAL, DIMENSION(D%NIT,D%NKT),INTENT(IN) :: PRW ! grid scale total water ! mixing ratio -REAL ,DIMENSION(D%NIT,D%NKT) ,INTENT(IN) :: PRC ! grid scale r_c -REAL ,DIMENSION(D%NIT,D%NKT) ,INTENT(IN) :: PRI ! grid scale r_i -LOGICAL ,DIMENSION(D%NIT) ,INTENT(IN) :: OTRIG1 ! logical to keep trace of -REAL ,DIMENSION(D%NIT,D%NKT) ,INTENT(OUT) :: PTHC ! conv. adj. grid scale theta -REAL ,DIMENSION(D%NIT,D%NKT) ,INTENT(OUT) :: PRWC ! conv. adj. grid scale r_w -REAL ,DIMENSION(D%NIT,D%NKT) ,INTENT(OUT) :: PRCC ! conv. adj. grid scale r_c -REAL ,DIMENSION(D%NIT,D%NKT) ,INTENT(OUT) :: PRIC ! conv. adj. grid scale r_i -REAL ,DIMENSION(D%NIT,D%NKT) ,INTENT(OUT) :: PWSUB ! envir. compensating subsidence(Pa/s) -INTEGER ,DIMENSION(D%NIT) ,INTENT(IN) :: KLCL ! index lifting condens. level -INTEGER ,DIMENSION(D%NIT) ,INTENT(IN) :: KDPL ! index for departure level -INTEGER ,DIMENSION(D%NIT) ,INTENT(IN) :: KPBL ! index for top of source layer -INTEGER ,DIMENSION(D%NIT) ,INTENT(IN) :: KCTL ! index for cloud top level -REAL ,DIMENSION(D%NIT,D%NKT) ,INTENT(INOUT) :: PUMF ! updraft mass flux (kg/s) -REAL ,DIMENSION(D%NIT,D%NKT) ,INTENT(INOUT) :: PUER ! updraft entrainment (kg/s) -REAL ,DIMENSION(D%NIT,D%NKT) ,INTENT(INOUT) :: PUDR ! updraft detrainment (kg/s) -REAL ,DIMENSION(D%NIT,D%NKT) ,INTENT(IN) :: PUTHL ! updraft enthalpy (J/kg) -REAL ,DIMENSION(D%NIT,D%NKT) ,INTENT(IN) :: PURW ! updraft total water (kg/kg) -REAL ,DIMENSION(D%NIT,D%NKT) ,INTENT(IN) :: PURC ! updraft cloud water (kg/kg) -REAL ,DIMENSION(D%NIT,D%NKT) ,INTENT(IN) :: PURI ! updraft cloud ice (kg/kg) -REAL ,DIMENSION(D%NIT) ,INTENT(IN) :: PCAPE ! available potent. energy -REAL ,DIMENSION(D%NIT) ,INTENT(INOUT) :: PTIMEC ! convection time step +REAL, DIMENSION(D%NIT,D%NKT),INTENT(IN) :: PRC ! grid scale r_c +REAL, DIMENSION(D%NIT,D%NKT),INTENT(IN) :: PRI ! grid scale r_i +LOGICAL, DIMENSION(D%NIT), INTENT(IN) :: OTRIG1 ! logical to keep trace of ! convective arrays modified in UPDRAFT ! ! -INTEGER ,INTENT(OUT) :: KFTSTEPS! maximum of fract time steps +REAL, DIMENSION(D%NIT,D%NKT), INTENT(IN) :: PPRES ! pressure (P) +REAL, DIMENSION(D%NIT,D%NKT), INTENT(IN) :: PDPRES ! pressure difference between + ! bottom and top of layer (Pa) +REAL, DIMENSION(D%NIT,D%NKT), INTENT(IN) :: PLMASS ! mass of model layer (kg) +REAL, DIMENSION(D%NIT,D%NKT), INTENT(IN) :: PZ ! height of model layer (m) +REAL, DIMENSION(D%NIT), INTENT(IN) :: PCAPE ! available potent. energy +INTEGER, INTENT(OUT) :: KFTSTEPS! maximum of fract time steps ! only used for chemical tracers ! -! -! + +REAL, DIMENSION(D%NIT,D%NKT), INTENT(INOUT):: PUMF ! updraft mass flux (kg/s) +REAL, DIMENSION(D%NIT,D%NKT), INTENT(INOUT):: PUER ! updraft entrainment (kg/s) +REAL, DIMENSION(D%NIT,D%NKT), INTENT(INOUT):: PUDR ! updraft detrainment (kg/s) +REAL, DIMENSION(D%NIT,D%NKT), INTENT(IN) :: PUTHL ! updraft enthalpy (J/kg) +REAL, DIMENSION(D%NIT,D%NKT), INTENT(IN) :: PURW ! updraft total water (kg/kg) +REAL, DIMENSION(D%NIT,D%NKT), INTENT(IN) :: PURC ! updraft cloud water (kg/kg) +REAL, DIMENSION(D%NIT,D%NKT), INTENT(IN) :: PURI ! updraft cloud ice (kg/kg) +REAL, DIMENSION(D%NIT,D%NKT), INTENT(OUT) :: PTHC ! conv. adj. grid scale theta +REAL, DIMENSION(D%NIT,D%NKT), INTENT(OUT) :: PRWC ! conv. adj. grid scale r_w +REAL, DIMENSION(D%NIT,D%NKT), INTENT(OUT) :: PRCC ! conv. adj. grid scale r_c +REAL, DIMENSION(D%NIT,D%NKT), INTENT(OUT) :: PRIC ! conv. adj. grid scale r_i +REAL, DIMENSION(D%NIT,D%NKT), INTENT(OUT) :: PWSUB ! envir. compensating subsidence(Pa/s) ! !* 0.2 Declarations of local variables : ! @@ -137,6 +136,7 @@ SUBROUTINE CONVECT_CLOSURE_SHAL( CVP_SHAL, CVPEXT, CST, D, & INTEGER :: JITER ! iteration loop index INTEGER :: JSTEP ! fractional time loop index REAL :: ZCPORD, ZRDOCP ! C_pd / R_d, R_d / C_pd +REAL :: ZEPS ! REAL, DIMENSION(D%NIT,D%NKT) :: ZTHLC ! convectively adjusted ! grid scale enthalpy @@ -168,20 +168,21 @@ SUBROUTINE CONVECT_CLOSURE_SHAL( CVP_SHAL, CVPEXT, CST, D, & INTEGER, DIMENSION(D%NIT) :: ICOUNT ! timestep counter INTEGER, DIMENSION(D%NIT) :: ILCL ! index lifting condens. level INTEGER, DIMENSION(D%NIT) :: IWORK1 ! work array -REAL, DIMENSION(D%NIT) :: ZWORK1, ZWORK2, ZWORK3, ZWORK4, ZWORK5 +REAL, DIMENSION(D%NIT) :: ZWORK1, ZWORK2, ZWORK3, ZWORK5 LOGICAL, DIMENSION(D%NIT) :: GWORK1, GWORK3! work arrays LOGICAL, DIMENSION(D%NIT,D%NKT) :: GWORK4 ! work array ! +REAL(KIND=JPHOOK) :: ZHOOK_HANDLE ! +#include "convect_closure_adjust_shal.h" +#include "convect_closure_thrvlcl.h" !------------------------------------------------------------------------------- ! !* 0.2 Initialize local variables ! ---------------------------- ! ! -REAL(KIND=JPHOOK) :: ZHOOK_HANDLE -#include "convect_closure_adjust_shal.h" IF (LHOOK) CALL DR_HOOK('CONVECT_CLOSURE_SHAL',0,ZHOOK_HANDLE) ZTIMC(:,:) = 0. @@ -189,12 +190,12 @@ SUBROUTINE CONVECT_CLOSURE_SHAL( CVP_SHAL, CVPEXT, CST, D, & ZWORK1(:) = 0. ZWORK2(:) = 0. ZWORK3(:) = 0. -ZWORK4(:) = 0. GWORK1(:) = .FALSE. GWORK3(:) = .FALSE. GWORK4(:,:) = .FALSE. ILCL(:) = KLCL(:) ! +ZEPS = CST%XRD / CST%XRV ZCPORD = CST%XCPD / CST%XRD ZRDOCP = CST%XRD / CST%XCPD ! @@ -244,7 +245,7 @@ SUBROUTINE CONVECT_CLOSURE_SHAL( CVP_SHAL, CVPEXT, CST, D, & END DO ! ! -GWORK1(D%NIB:D%NIE) = OTRIG1(D%NIB:D%NIE) ! logical array to limit adjustment to not definitively +GWORK1(:) = OTRIG1(:) ! logical array to limit adjustment to not definitively ! adjusted columns ! DO JK = IKB, IKE @@ -262,18 +263,16 @@ SUBROUTINE CONVECT_CLOSURE_SHAL( CVP_SHAL, CVPEXT, CST, D, & DO JITER = 1, 4 ! Enter adjustment loop to assure that all CAPE is ! removed within the advective time interval TIMEC ! - ZTIMEC(D%NIB:D%NIE) = PTIMEC(D%NIB:D%NIE) - DO JI=1, D%NIE + ZTIMEC(:) = PTIMEC(:) DO JK=1, IKS - GWORK4(JI,JK) = GWORK1(JI) - ENDDO + GWORK4(1:D%NIE,JK) = GWORK1(1:D%NIE) ENDDO DO JK = IKB, IKE DO JI=D%NIB, D%NIE - IF(GWORK4(JI,JK))PWSUB(JI,JK) = 0. + IF(GWORK4(JI,JK)) PWSUB(JI,JK) = 0. ENDDO END DO - ZOMG(D%NIB:D%NIE,1:D%NKT)=0. + ZOMG(:,1:D%NKT)=0. ! DO JK = IKB + 1, JKMAX JKP = MAX( IKB + 1, JK - 1 ) @@ -341,22 +340,18 @@ SUBROUTINE CONVECT_CLOSURE_SHAL( CVP_SHAL, CVPEXT, CST, D, & DO JK = IKB, IKE PWSUB(:,JK) = PWSUB(:,JK) * ZWORK5(:) END DO - GWORK4(D%NIB:D%NIE,1:IKB) = .FALSE. - GWORK4(D%NIB:D%NIE,IKE:IKS) = .FALSE. + GWORK4(:,1:IKB) = .FALSE. + GWORK4(:,IKE:IKS) = .FALSE. ! DO JI=D%NIB,D%NIE ITSTEP(JI) = INT( PTIMEC(JI) / ZTIMEC(JI) ) + 1 - ENDDO - DO JI=D%NIB,D%NIE ZTIMEC(JI) = PTIMEC(JI) / REAL( ITSTEP(JI) ) ! adjust fractional time step ENDDO ! to be an integer multiple of PTIMEC - DO JI=1, D%NIE DO JK=1, IKS - ZTIMC(JI,JK) = ZTIMEC(JI) - ENDDO + ZTIMC(1:D%NIE,JK) = ZTIMEC(1:D%NIE) ENDDO - ICOUNT(D%NIB:D%NIE) = 0 + ICOUNT(:) = 0 ! ! ! @@ -366,63 +361,60 @@ SUBROUTINE CONVECT_CLOSURE_SHAL( CVP_SHAL, CVPEXT, CST, D, & ENDDO DO JSTEP = 1, KFTSTEPS ! Enter the fractional time step loop here ! - ICOUNT(D%NIB:D%NIE) = ICOUNT(D%NIB:D%NIE) + 1 -! - GWORK3(D%NIB:D%NIE) = ITSTEP(D%NIB:D%NIE) >= ICOUNT(D%NIB:D%NIE) .AND. GWORK1(D%NIB:D%NIE) + DO JI=D%NIB,D%NIE + ICOUNT(JI) = ICOUNT(JI) + 1 + GWORK3(JI) = ITSTEP(JI) >= ICOUNT(JI) .AND. GWORK1(JI) + ENDDO ! ! !* 7. Assign enthalpy and r_w values at the top and bottom of each ! layer based on the sign of w ! ------------------------------------------------------------ ! - ZTHMFIN(:,:) = 0. - ZRWMFIN(:,:) = 0. - ZRCMFIN(:,:) = 0. - ZRIMFIN(:,:) = 0. - ZTHMFOUT(:,:) = 0. - ZRWMFOUT(:,:) = 0. - ZRCMFOUT(:,:) = 0. - ZRIMFOUT(:,:) = 0. -! - DO JK = IKB + 1, JKMAX - DO JI = D%NIB, D%NIE - GWORK4(JI,JK) = GWORK3(JI) .AND. JK <= KCTL(JI) - END DO - JKP = MAX( IKB + 1, JK - 1 ) - DO JI = D%NIB, D%NIE - IF ( GWORK3(JI) ) THEN -! - ZWORK1(JI) = SIGN( 1., ZOMG(JI,JK) ) - ZWORK2(JI) = 0.5 * ( 1. + ZWORK1(JI) ) - ZWORK1(JI) = 0.5 * ( 1. - ZWORK1(JI) ) - ZTHMFIN(JI,JK) = - ZOMG(JI,JK) * ZTHLC(JI,JKP) * ZWORK1(JI) - ZTHMFOUT(JI,JK) = ZOMG(JI,JK) * ZTHLC(JI,JK) * ZWORK2(JI) - ZRWMFIN(JI,JK) = - ZOMG(JI,JK) * PRWC(JI,JKP) * ZWORK1(JI) - ZRWMFOUT(JI,JK) = ZOMG(JI,JK) * PRWC(JI,JK) * ZWORK2(JI) - ZRCMFIN(JI,JK) = - ZOMG(JI,JK) * PRCC(JI,JKP) * ZWORK1(JI) - ZRCMFOUT(JI,JK) = ZOMG(JI,JK) * PRCC(JI,JK) * ZWORK2(JI) - ZRIMFIN(JI,JK) = - ZOMG(JI,JK) * PRIC(JI,JKP) * ZWORK1(JI) - ZRIMFOUT(JI,JK) = ZOMG(JI,JK) * PRIC(JI,JK) * ZWORK2(JI) - END IF - END DO - DO JI = D%NIB, D%NIE - IF ( GWORK3(JI) ) THEN - ZTHMFIN(JI,JKP) = ZTHMFIN(JI,JKP) + ZTHMFOUT(JI,JK) * ZWORK2(JI) - ZTHMFOUT(JI,JKP) = ZTHMFOUT(JI,JKP) + ZTHMFIN(JI,JK) * ZWORK1(JI) - ZRWMFIN(JI,JKP) = ZRWMFIN(JI,JKP) + ZRWMFOUT(JI,JK) * ZWORK2(JI) - ZRWMFOUT(JI,JKP) = ZRWMFOUT(JI,JKP) + ZRWMFIN(JI,JK) * ZWORK1(JI) - ZRCMFIN(JI,JKP) = ZRCMFIN(JI,JKP) + ZRCMFOUT(JI,JK) * ZWORK2(JI) - ZRCMFOUT(JI,JKP) = ZRCMFOUT(JI,JKP) + ZRCMFIN(JI,JK) * ZWORK1(JI) - ZRIMFIN(JI,JKP) = ZRIMFIN(JI,JKP) + ZRIMFOUT(JI,JK) * ZWORK2(JI) - ZRIMFOUT(JI,JKP) = ZRIMFOUT(JI,JKP) + ZRIMFIN(JI,JK) * ZWORK1(JI) -! - END IF - END DO - END DO -! - DO JK = IKB, IKE - DO JI=D%NIB, D%NIE - IF(GWORK4(JI,JK)) THEN + DO JK = 1, D%NKT + ZTHMFIN(:,JK) = 0. + ZRWMFIN(:,JK) = 0. + ZRCMFIN(:,JK) = 0. + ZRIMFIN(:,JK) = 0. + ZTHMFOUT(:,JK) = 0. + ZRWMFOUT(:,JK) = 0. + ZRCMFOUT(:,JK) = 0. + ZRIMFOUT(:,JK) = 0. + ENDDO +! + DO JK = IKB + 1, JKMAX + DO JI = D%NIB, D%NIE + GWORK4(JI,JK) = GWORK3(JI) .AND. JK <= KCTL(JI) + ENDDO + JKP = MAX( IKB + 1, JK - 1 ) + DO JI = D%NIB, D%NIE + IF ( GWORK3(JI) ) THEN + ZWORK1(JI) = SIGN( 1., ZOMG(JI,JK) ) + ZWORK2(JI) = 0.5 * ( 1. + ZWORK1(JI) ) + ZWORK1(JI) = 0.5 * ( 1. - ZWORK1(JI) ) + ZTHMFIN(JI,JK) = - ZOMG(JI,JK) * ZTHLC(JI,JKP) * ZWORK1(JI) + ZTHMFOUT(JI,JK) = ZOMG(JI,JK) * ZTHLC(JI,JK) * ZWORK2(JI) + ZRWMFIN(JI,JK) = - ZOMG(JI,JK) * PRWC(JI,JKP) * ZWORK1(JI) + ZRWMFOUT(JI,JK) = ZOMG(JI,JK) * PRWC(JI,JK) * ZWORK2(JI) + ZRCMFIN(JI,JK) = - ZOMG(JI,JK) * PRCC(JI,JKP) * ZWORK1(JI) + ZRCMFOUT(JI,JK) = ZOMG(JI,JK) * PRCC(JI,JK) * ZWORK2(JI) + ZRIMFIN(JI,JK) = - ZOMG(JI,JK) * PRIC(JI,JKP) * ZWORK1(JI) + ZRIMFOUT(JI,JK) = ZOMG(JI,JK) * PRIC(JI,JK) * ZWORK2(JI) + ENDIF + ENDDO + DO JI = D%NIB, D%NIE + IF ( GWORK3(JI) ) THEN + ZTHMFIN(JI,JKP) = ZTHMFIN(JI,JKP) + ZTHMFOUT(JI,JK) * ZWORK2(JI) + ZTHMFOUT(JI,JKP) = ZTHMFOUT(JI,JKP) + ZTHMFIN(JI,JK) * ZWORK1(JI) + ZRWMFIN(JI,JKP) = ZRWMFIN(JI,JKP) + ZRWMFOUT(JI,JK) * ZWORK2(JI) + ZRWMFOUT(JI,JKP) = ZRWMFOUT(JI,JKP) + ZRWMFIN(JI,JK) * ZWORK1(JI) + ZRCMFIN(JI,JKP) = ZRCMFIN(JI,JKP) + ZRCMFOUT(JI,JK) * ZWORK2(JI) + ZRCMFOUT(JI,JKP) = ZRCMFOUT(JI,JKP) + ZRCMFIN(JI,JK) * ZWORK1(JI) + ZRIMFIN(JI,JKP) = ZRIMFIN(JI,JKP) + ZRIMFOUT(JI,JK) * ZWORK2(JI) + ZRIMFOUT(JI,JKP) = ZRIMFOUT(JI,JKP) + ZRIMFIN(JI,JK) * ZWORK1(JI) + ENDIF + ENDDO + ENDDO ! !****************************************************************************** ! @@ -430,6 +422,10 @@ SUBROUTINE CONVECT_CLOSURE_SHAL( CVP_SHAL, CVPEXT, CST, D, & ! NOTA: These are the MAIN EQUATIONS of the scheme ! ----------------------------------------------------------------- ! +! + DO JK = IKB, IKE + DO JI=D%NIB, D%NIE + IF(GWORK4(JI,JK)) THEN ! ZTHLC(JI,JK) = ZTHLC(JI,JK) + ZTIMC(JI,JK) / PLMASS(JI,JK) * ( & ZTHMFIN(JI,JK) + PUDR(JI,JK) * PUTHL(JI,JK) & @@ -443,36 +439,36 @@ SUBROUTINE CONVECT_CLOSURE_SHAL( CVP_SHAL, CVPEXT, CST, D, & PRIC(JI,JK) = PRIC(JI,JK) + ZTIMC(JI,JK) / PLMASS(JI,JK) * ( & ZRIMFIN(JI,JK) + PUDR(JI,JK) * PURI(JI,JK) - ZRIMFOUT(JI,JK) - & PUER(JI,JK) * MAX(0., PRI(JI,JK)) ) + ENDIF + ENDDO + ENDDO ! ! !****************************************************************************** ! - ENDIF - ENDDO - ENDDO ! - END DO ! Exit the fractional time step loop + ENDDO ! Exit the fractional time step loop ! ! !* 10. Compute final linearized value of theta envir. ! ---------------------------------------------- ! - DO JK = IKB + 1, JKMAX - DO JI = D%NIB, D%NIE - IF( GWORK1(JI) .AND. JK <= KCTL(JI) ) THEN - ZPI(JI) = ( CST%XP00 / PPRES(JI,JK) ) ** ZRDOCP - ZCPH(JI) = CST%XCPD + PRWC(JI,JK) * CST%XCPV - ZWORK2(JI) = PTH(JI,JK) / ZPI(JI) ! first temperature estimate - ZLV(JI) = CST%XLVTT + ( CST%XCPV - CST%XCL ) * ( ZWORK2(JI) - CST%XTT ) - ZLS(JI) = CST%XLVTT + ( CST%XCPV - CST%XCI ) * ( ZWORK2(JI) - CST%XTT ) - ! final linearized temperature - ZWORK2(JI) = ( ZTHLC(JI,JK) + ZLV(JI) * PRCC(JI,JK) + ZLS(JI) * PRIC(JI,JK) & - - (1. + PRWC(JI,JK) ) * CST%XG * PZ(JI,JK) ) / ZCPH(JI) - ZWORK2(JI) = MAX( 180., MIN( 340., ZWORK2(JI) ) ) - PTHC(JI,JK)= ZWORK2(JI) * ZPI(JI) ! final adjusted envir. theta - END IF - END DO - END DO + DO JK = IKB + 1, JKMAX + DO JI = D%NIB, D%NIE + IF( GWORK1(JI) .AND. JK <= KCTL(JI) ) THEN + ZPI(JI) = ( CST%XP00 / PPRES(JI,JK) ) ** ZRDOCP + ZCPH(JI) = CST%XCPD + PRWC(JI,JK) * CST%XCPV + ZWORK2(JI) = PTH(JI,JK) / ZPI(JI) ! first temperature estimate + ZLV(JI) = CST%XLVTT + ( CST%XCPV - CST%XCL ) * ( ZWORK2(JI) - CST%XTT ) + ZLS(JI) = CST%XLVTT + ( CST%XCPV - CST%XCI ) * ( ZWORK2(JI) - CST%XTT ) + ! final linearized temperature + ZWORK2(JI) = ( ZTHLC(JI,JK) + ZLV(JI) * PRCC(JI,JK) + ZLS(JI) * PRIC(JI,JK) & + - (1. + PRWC(JI,JK) ) * CST%XG * PZ(JI,JK) ) / ZCPH(JI) + ZWORK2(JI) = MAX( 180., MIN( 340., ZWORK2(JI) ) ) + PTHC(JI,JK)= ZWORK2(JI) * ZPI(JI) ! final adjusted envir. theta + ENDIF + ENDDO + ENDDO ! ! !* 11. Compute new cloud ( properties at new LCL ) @@ -480,105 +476,99 @@ SUBROUTINE CONVECT_CLOSURE_SHAL( CVP_SHAL, CVPEXT, CST, D, & ! that in routine TRIGGER_FUNCT ! --------------------------------------------- ! - CALL CONVECT_CLOSURE_THRVLCL( CVPEXT, CST, D, & + CALL CONVECT_CLOSURE_THRVLCL( CVPEXT, CST, D, & PPRES, PTHC, PRWC, PZ, GWORK1, & ZTHLCL, ZRVLCL, ZZLCL, ZTLCL, ZTELCL, & ILCL, KDPL, KPBL ) ! ! - ZTLCL(D%NIB:D%NIE) = MAX( 230., MIN( 335., ZTLCL(D%NIB:D%NIE) ) ) ! set some overflow bounds - ZTELCL(D%NIB:D%NIE) = MAX( 230., MIN( 335., ZTELCL(D%NIB:D%NIE) ) ) - ZTHLCL(D%NIB:D%NIE) = MAX( 230., MIN( 345., ZTHLCL(D%NIB:D%NIE) ) ) - ZRVLCL(D%NIB:D%NIE) = MAX( 0., MIN( 1., ZRVLCL(D%NIB:D%NIE) ) ) -! +!DIR$ IVDEP + DO JI = D%NIB, D%NIE + ZTLCL (JI) = MAX( 230., MIN( 335., ZTLCL(JI) ) ) ! set some overflow bounds + ZTELCL(JI) = MAX( 230., MIN( 335., ZTELCL(JI) ) ) + ZTHLCL(JI) = MAX( 230., MIN( 345., ZTHLCL(JI) ) ) + ZRVLCL(JI) = MAX( 0., MIN( 1., ZRVLCL(JI) ) ) + ZCAPE (JI) = 0. + ZPI (JI) = MAX( 0.95, MIN( 1.5, ZTHLCL(JI)/ZTLCL(JI) ) ) + ZWORK1(JI) = CST%XP00 / ZPI(JI) ** ZCPORD ! pressure at LCL + CALL CONVECT_SATMIXRATIO( ZWORK1(JI), ZTELCL(JI), ZEPS, ZWORK3(JI), ZLV(JI), ZLS(JI), ZCPH(JI) ) ! !* 12. Compute adjusted CAPE ! --------------------- + ZWORK3(JI) = MIN( .1, MAX( 0., ZWORK3(JI) ) ) + + + ! compute theta_e updraft undilute + ZTHEUL(JI) = ZTLCL(JI) * ZPI(JI) ** ( 1. - 0.28 * ZRVLCL(JI) ) & + & * EXP( ( 3374.6525 / ZTLCL(JI) - 2.5403 ) * ZRVLCL(JI) * ( 1. + 0.81 * ZRVLCL(JI) ) ) + ! compute theta_e saturated environment at LCL + ZTHES1(JI) = ZTELCL(JI) * ZPI(JI) ** ( 1. - 0.28 * ZWORK3(JI) ) & + & * EXP( ( 3374.6525 / ZTELCL(JI) - 2.5403 ) * ZWORK3(JI) * ( 1. + 0.81 * ZWORK3(JI) ) ) + ENDDO ! - ZCAPE(D%NIB:D%NIE) = 0. - ZPI(D%NIB:D%NIE) = ZTHLCL(D%NIB:D%NIE) / ZTLCL(D%NIB:D%NIE) - ZPI(D%NIB:D%NIE) = MAX( 0.95, MIN( 1.5, ZPI(D%NIB:D%NIE) ) ) - ZWORK1(D%NIB:D%NIE) = CST%XP00 / ZPI(D%NIB:D%NIE) ** ZCPORD ! pressure at LCL -! - CALL CONVECT_SATMIXRATIO( CST, D, ZWORK1, ZTELCL, ZWORK3, ZLV, ZLS, ZCPH ) - ZWORK3(D%NIB:D%NIE) = MIN( .1, MAX( 0., ZWORK3(D%NIB:D%NIE) ) ) -! - ! compute theta_e updraft undilute - ZTHEUL(D%NIB:D%NIE) = ZTLCL(D%NIB:D%NIE) * ZPI(D%NIB:D%NIE) ** ( 1. - 0.28 * ZRVLCL(D%NIB:D%NIE) ) & - * EXP( ( 3374.6525 / ZTLCL(D%NIB:D%NIE) - 2.5403 ) & - * ZRVLCL(D%NIB:D%NIE) * ( 1. + 0.81 * ZRVLCL(D%NIB:D%NIE) ) ) -! - ! compute theta_e saturated environment at LCL - ZTHES1(D%NIB:D%NIE) = ZTELCL(D%NIB:D%NIE) * ZPI(D%NIB:D%NIE) ** ( 1. - 0.28 * ZWORK3(D%NIB:D%NIE) ) & - * EXP( ( 3374.6525 / ZTELCL(D%NIB:D%NIE) - 2.5403 ) & - * ZWORK3(D%NIB:D%NIE) * ( 1. + 0.81 * ZWORK3(D%NIB:D%NIE) ) ) -! - DO JK = IKB, JKMAX - JKP = JK - 1 - DO JI = D%NIB, D%NIE - ZWORK4(JI) = 1. - IF ( JK == ILCL(JI) ) ZWORK4(JI) = 0. -! - ! compute theta_e saturated environment and adjusted values - ! of theta -! - GWORK3(JI) = JK >= ILCL(JI) .AND. JK <= KCTL(JI) .AND. GWORK1(JI) -! - ZPI(JI) = ( CST%XP00 / PPRES(JI,JK) ) ** ZRDOCP - ZWORK2(JI) = PTHC(JI,JK) / ZPI(JI) - END DO -! - CALL CONVECT_SATMIXRATIO( CST, D, PPRES(:,JK), ZWORK2, ZWORK3, ZLV, ZLS, ZCPH ) -! -! - DO JI = D%NIB, D%NIE - IF ( GWORK3(JI) ) THEN - ZTHES2(JI) = ZWORK2(JI) * ZPI(JI) ** ( 1. - 0.28 * ZWORK3(JI) ) & - * EXP( ( 3374.6525 / ZWORK2(JI) - 2.5403 ) & - * ZWORK3(JI) * ( 1. + 0.81 * ZWORK3(JI) ) ) -! - ZWORK3(JI) = PZ(JI,JK) - PZ(JI,JKP) * ZWORK4(JI) - & - ( 1. - ZWORK4(JI) ) * ZZLCL(JI) ! level thickness - ZWORK1(JI) = ( 2. * ZTHEUL(JI) ) / ( ZTHES1(JI) + ZTHES2(JI) ) - 1. - ZCAPE(JI) = ZCAPE(JI) + CST%XG * ZWORK3(JI) * MAX( 0., ZWORK1(JI) ) - ZTHES1(JI) = ZTHES2(JI) - END IF - END DO + DO JK = IKB, JKMAX + + JKP = JK - 1 +!DIR$ IVDEP + DO JI = D%NIB, D%NIE + ! compute theta_e saturated environment and adjusted values of theta + GWORK3(JI) = JK >= ILCL(JI) .AND. JK <= KCTL(JI) .AND. GWORK1(JI) + ZPI(JI) = ( CST%XP00 / PPRES(JI,JK) ) ** ZRDOCP + ZWORK2(JI) = PTHC(JI,JK) / ZPI(JI) + CALL CONVECT_SATMIXRATIO( PPRES(JI,JK), ZWORK2(JI), ZEPS, ZWORK3(JI), ZLV(JI), ZLS(JI), ZCPH(JI) ) END DO +! + DO JI = D%NIB, D%NIE + IF ( GWORK3(JI) ) THEN + ZTHES2(JI) = ZWORK2(JI) * ZPI(JI) ** ( 1. - 0.28 * ZWORK3(JI) ) & + & * EXP( ( 3374.6525 / ZWORK2(JI) - 2.5403 ) * ZWORK3(JI) * ( 1. + 0.81 * ZWORK3(JI) ) ) + IF ( JK == ILCL(JI) ) THEN + ZWORK3(JI) = PZ(JI,JK) - ZZLCL(JI) ! level thickness + ELSE + ZWORK3(JI) = PZ(JI,JK) - PZ(JI,JKP) ! level thickness + ENDIF + ZWORK1(JI) = ( 2. * ZTHEUL(JI) ) / ( ZTHES1(JI) + ZTHES2(JI) ) - 1. + ZCAPE(JI) = ZCAPE(JI) + CST%XG * ZWORK3(JI) * MAX( 0., ZWORK1(JI) ) + ZTHES1(JI) = ZTHES2(JI) + ENDIF + ENDDO + + ENDDO ! ! !* 13. Determine mass adjustment factor knowing how much ! CAPE has been removed. ! ------------------------------------------------- ! - DO JI=D%NIB,D%NIE - IF ( GWORK1(JI) ) THEN - ZWORK1(JI) = MAX( PCAPE(JI) - ZCAPE(JI), 0.2 * PCAPE(JI) ) - ZWORK2(JI) = ZCAPE(JI) / ( PCAPE(JI) + 1.E-8 ) + DO JI=D%NIB,D%NIE + IF ( GWORK1(JI) ) THEN + ZWORK1(JI) = MAX( PCAPE(JI) - ZCAPE(JI), 0.2 * PCAPE(JI) ) + ZWORK2(JI) = ZCAPE(JI) / ( PCAPE(JI) + 1.E-8 ) ! - GWORK1(JI) = ZWORK2(JI) > 0.2 .OR. ZCAPE(JI) == 0. ! mask for adjustment - END IF - ENDDO + GWORK1(JI) = ZWORK2(JI) > 0.2 .OR. ZCAPE(JI) == 0. ! mask for adjustment + ENDIF + ENDDO ! - DO JI=D%NIB,D%NIE - IF( ZCAPE(JI) == 0. .AND. GWORK1(JI) ) ZADJ(JI) = ZADJ(JI) * 0.5 - ENDDO - DO JI=D%NIB,D%NIE - IF ( ZCAPE(JI) /= 0. .AND. GWORK1(JI) ) THEN - ZADJ(JI) = ZADJ(JI) * CVP_SHAL%XSTABC * PCAPE(JI) / ( ZWORK1(JI) + 1.E-8 ) - ENDIF - ENDDO - DO JI=D%NIB,D%NIE - ZADJ(JI) = MIN( ZADJ(JI), ZADJMAX(JI) ) - ENDDO + DO JI=D%NIB,D%NIE + IF( ZCAPE(JI) == 0. .AND. GWORK1(JI) ) THEN + ZADJ(JI) = ZADJ(JI) * 0.5 + ENDIF + ENDDO + DO JI=D%NIB,D%NIE + IF ( ZCAPE(JI) /= 0. .AND. GWORK1(JI) ) THEN + ZADJ(JI) = ZADJ(JI) * CVP_SHAL%XSTABC * PCAPE(JI) / ( ZWORK1(JI) + 1.E-8 ) + ENDIF + ENDDO + DO JI=D%NIB,D%NIE + ZADJ(JI) = MIN( ZADJ(JI), ZADJMAX(JI) ) + ENDDO ! ! !* 13. Adjust mass flux by the factor ZADJ to converge to ! specified degree of stabilization ! ---------------------------------------------------- ! - CALL CONVECT_CLOSURE_ADJUST_SHAL( CVPEXT, D, ZADJ,& - PUMF, ZUMF, PUER, ZUER, PUDR, ZUDR ) + CALL CONVECT_CLOSURE_ADJUST_SHAL( CVPEXT, D, ZADJ, PUMF, ZUMF, PUER, ZUER, PUDR, ZUDR ) ! ! !IF ( COUNT( GWORK1(:) ) == 0 ) EXIT ! exit big adjustment iteration loop @@ -597,5 +587,7 @@ SUBROUTINE CONVECT_CLOSURE_SHAL( CVP_SHAL, CVPEXT, CST, D, & ! ! IF (LHOOK) CALL DR_HOOK('CONVECT_CLOSURE_SHAL',1,ZHOOK_HANDLE) +CONTAINS +INCLUDE "convect_satmixratio.h" +! END SUBROUTINE CONVECT_CLOSURE_SHAL - diff --git a/conv/convect_closure_thrvlcl.F90 b/conv/convect_closure_thrvlcl.F90 index 860a26fc..064f43bf 100644 --- a/conv/convect_closure_thrvlcl.F90 +++ b/conv/convect_closure_thrvlcl.F90 @@ -76,23 +76,25 @@ SUBROUTINE CONVECT_CLOSURE_THRVLCL( CVPEXT, CST, D, & ! !* 0.1 Declarations of dummy arguments : ! -TYPE(CONVPAREXT) ,INTENT(IN) :: CVPEXT -TYPE(CST_T) ,INTENT(IN) :: CST -TYPE(DIMPHYEX_T) ,INTENT(IN) :: D -REAL ,DIMENSION(D%NIT,D%NKT) ,INTENT(IN) :: PPRES ! pressure -REAL ,DIMENSION(D%NIT,D%NKT) ,INTENT(IN) :: PTH ! theta -REAL ,DIMENSION(D%NIT,D%NKT) ,INTENT(IN) :: PRV ! vapor mixing ratio -REAL ,DIMENSION(D%NIT,D%NKT) ,INTENT(IN) :: PZ ! height of grid point (m) -LOGICAL ,DIMENSION(D%NIT) ,INTENT(IN) :: OWORK1! logical mask -REAL ,DIMENSION(D%NIT) ,INTENT(OUT) :: PTHLCL ! theta at LCL -REAL ,DIMENSION(D%NIT) ,INTENT(OUT) :: PRVLCL ! vapor mixing ratio at LCL -REAL ,DIMENSION(D%NIT) ,INTENT(OUT) :: PZLCL ! height at LCL (m) -REAL ,DIMENSION(D%NIT) ,INTENT(OUT) :: PTLCL ! temperature at LCL (m) -REAL ,DIMENSION(D%NIT) ,INTENT(OUT) :: PTELCL ! environm. temp. at LCL (K) -INTEGER ,DIMENSION(D%NIT) ,INTENT(OUT) :: KLCL ! contains vert. index of LCL -INTEGER ,DIMENSION(D%NIT) ,INTENT(IN) :: KDPL ! contains vert. index of DPL -INTEGER ,DIMENSION(D%NIT) ,INTENT(IN) :: KPBL ! " vert. index of source layer top -! +TYPE(CONVPAREXT), INTENT(IN) :: CVPEXT +TYPE(CST_T), INTENT(IN) :: CST +TYPE(DIMPHYEX_T), INTENT(IN) :: D +REAL, DIMENSION(D%NIT,D%NKT), INTENT(IN) :: PTH ! theta +REAL, DIMENSION(D%NIT,D%NKT), INTENT(IN) :: PRV ! vapor mixing ratio +REAL, DIMENSION(D%NIT,D%NKT), INTENT(IN) :: PPRES ! pressure +REAL, DIMENSION(D%NIT,D%NKT), INTENT(IN) :: PZ ! height of grid point (m) +INTEGER, DIMENSION(D%NIT), INTENT(IN) :: KDPL ! contains vert. index of DPL +INTEGER, DIMENSION(D%NIT), INTENT(IN) :: KPBL ! " vert. index of source layer top +LOGICAL, DIMENSION(D%NIT), INTENT(IN) :: OWORK1! logical mask +! +REAL, DIMENSION(D%NIT), INTENT(OUT):: PTHLCL ! theta at LCL +REAL, DIMENSION(D%NIT), INTENT(OUT):: PRVLCL ! vapor mixing ratio at LCL +REAL, DIMENSION(D%NIT), INTENT(OUT):: PZLCL ! height at LCL (m) +REAL, DIMENSION(D%NIT), INTENT(OUT):: PTLCL ! temperature at LCL (m) +REAL, DIMENSION(D%NIT), INTENT(OUT):: PTELCL ! environm. temp. at LCL (K) +INTEGER, DIMENSION(D%NIT), INTENT(OUT):: KLCL ! contains vert. index of LCL + + ! !* 0.2 Declarations of local variables : ! @@ -110,13 +112,13 @@ SUBROUTINE CONVECT_CLOSURE_THRVLCL( CVPEXT, CST, D, & REAL, DIMENSION(D%NIT) :: ZDP ! pressure between LCL and model layer REAL, DIMENSION(D%NIT) :: ZWORK1, ZWORK2 ! work arrays ! +REAL(KIND=JPHOOK) :: ZHOOK_HANDLE ! !------------------------------------------------------------------------------- ! !* 0.3 Compute array bounds ! -------------------- ! -REAL(KIND=JPHOOK) :: ZHOOK_HANDLE IF (LHOOK) CALL DR_HOOK('CONVECT_CLOSURE_THRVLCL',0,ZHOOK_HANDLE) IKB = 1 + CVPEXT%JCVEXB IKE = D%NKT - CVPEXT%JCVEXT @@ -191,14 +193,18 @@ SUBROUTINE CONVECT_CLOSURE_THRVLCL( CVPEXT, CST, D, & END IF ENDDO ! - ZPLCL(D%NIB:D%NIE) = MIN( 2.E5, MAX( 10., ZPLCL(D%NIB:D%NIE) ) ) ! bound to avoid overflow +DO JI = D%NIB,D%NIE + ZPLCL(JI) = MIN( 2.E5, MAX( 10., ZPLCL(JI) ) ) ! bound to avoid overflow +ENDDO ! ! !* 3.2 Correct PTLCL in order to be completely consistent ! with MNH saturation formula ! -------------------------------------------------- ! - CALL CONVECT_SATMIXRATIO( CST, D, ZPLCL, PTLCL, ZWORK1, ZLV, ZWORK2, ZCPH ) + DO JI=D%NIB,D%NIE + CALL CONVECT_SATMIXRATIO( ZPLCL(JI), PTLCL(JI), ZEPS, ZWORK1(JI), ZLV(JI), ZWORK2(JI), ZCPH(JI) ) + ENDDO DO JI=D%NIB,D%NIE IF( OWORK1(JI) ) THEN ZWORK2(JI) = ZWORK1(JI) / PTLCL(JI) * ( CST%XBETAW / PTLCL(JI) - CST%XGAMW ) ! dr_sat/dT @@ -213,7 +219,9 @@ SUBROUTINE CONVECT_CLOSURE_THRVLCL( CVPEXT, CST, D, & ! to saturation values. ! ------------------------------------------------------- ! - CALL CONVECT_SATMIXRATIO( CST, D, ZPRESMIX, ZTMIX, ZWORK1, ZLV, ZWORK2, ZCPH ) + DO JI=D%NIB,D%NIE + CALL CONVECT_SATMIXRATIO( ZPRESMIX(JI), ZTMIX(JI), ZEPS, ZWORK1(JI), ZLV(JI), ZWORK2(JI), ZCPH(JI) ) + ENDDO DO JI=D%NIB,D%NIE IF( OWORK1(JI) .AND. PRVLCL(JI) > ZWORK1(JI) ) THEN ZWORK2(JI) = ZWORK1(JI) / ZTMIX(JI) * ( CST%XBETAW / ZTMIX(JI) - CST%XGAMW ) ! dr_sat/dT @@ -265,5 +273,7 @@ SUBROUTINE CONVECT_CLOSURE_THRVLCL( CVPEXT, CST, D, & ! ! IF (LHOOK) CALL DR_HOOK('CONVECT_CLOSURE_THRVLCL',1,ZHOOK_HANDLE) +CONTAINS +INCLUDE "convect_satmixratio.h" +! END SUBROUTINE CONVECT_CLOSURE_THRVLCL - diff --git a/conv/convect_downdraft.F90 b/conv/convect_downdraft.F90 index 1193a6ae..af7b3f83 100644 --- a/conv/convect_downdraft.F90 +++ b/conv/convect_downdraft.F90 @@ -119,6 +119,7 @@ SUBROUTINE CONVECT_DOWNDRAFT( KLON, KLEV, & REAL, DIMENSION(KLON,KLEV), INTENT(OUT):: PDTEVRF! downdraft evaporation rate INTEGER, DIMENSION(KLON), INTENT(OUT):: KLFS ! contains vert. index of LFS INTEGER, DIMENSION(KLON), INTENT(OUT):: KDBL ! contains vert. index of DBL + ! !* 0.2 Declarations of local variables : ! @@ -139,13 +140,13 @@ SUBROUTINE CONVECT_DOWNDRAFT( KLON, KLEV, & REAL, DIMENSION(KLON) :: ZWORK1, ZWORK2, ZWORK3, ZWORK4 ! work arrays LOGICAL, DIMENSION(KLON) :: GWORK1 ! work array ! +REAL(KIND=JPHOOK) :: ZHOOK_HANDLE ! !------------------------------------------------------------------------------- ! ! 0.3 Set loop bounds ! --------------- ! -REAL(KIND=JPHOOK) :: ZHOOK_HANDLE IF (LHOOK) CALL DR_HOOK('CONVECT_DOWNDRAFT',0,ZHOOK_HANDLE) IIE = KLON IKB = 1 + JCVEXB @@ -263,7 +264,9 @@ SUBROUTINE CONVECT_DOWNDRAFT( KLON, KLEV, & ! --------------------------------------------------------- ! ! compute satur. mixing ratio for melting corrected temperature -CALL CONVECT_SATMIXRATIO( KLON, ZWORK4, ZDT, ZWORK3, ZLV, ZLS, ZCPH ) + DO JI = 1, IIE + CALL CONVECT_SATMIXRATIO( ZWORK4(JI), ZDT(JI), ZEPS, ZWORK3(JI), ZLV(JI), ZLS(JI), ZCPH(JI) ) + END DO ! ! compute envir. saturated theta_e for melting corrected temperature ZWORK1(:) = MIN( ZWORK2(:), ZWORK3(:) ) @@ -382,7 +385,9 @@ SUBROUTINE CONVECT_DOWNDRAFT( KLON, KLEV, & ! -------------------------------------------------- ! DO JITER = 1, 4 - CALL CONVECT_SATMIXRATIO( KLON, PPRES(:,JK), ZDT, ZWORK1, ZLV, ZLS, ZCPH ) + DO JI = 1, IIE + CALL CONVECT_SATMIXRATIO( PPRES(JI,JK), ZDT(JI), ZEPS, ZWORK1(JI), ZLV(JI), ZLS(JI), ZCPH(JI) ) + END DO ZDTP(:) = PDTHL(:,JK) / ( ZPI(:) ** ( 1. - 0.28 * ZWORK1(:) ) & * EXP( ( 3374.6525 / ZDT(:) - 2.5403 ) & * ZWORK1(:) * ( 1. + 0.81 * ZWORK1(:) ) ) ) @@ -399,7 +404,9 @@ SUBROUTINE CONVECT_DOWNDRAFT( KLON, KLEV, & ( 1. + ZLV(:) / ZCPH(:) * ZWORK2(:) ) ! temperature perturb ! due to evaporation ZDT(:) = ZDT(:) + ZWORK2(:) ! - CALL CONVECT_SATMIXRATIO( KLON, PPRES(:,JK), ZDT, ZWORK3, ZLV, ZLS, ZCPH ) + DO JI = 1, IIE + CALL CONVECT_SATMIXRATIO( PPRES(JI,JK), ZDT(JI), ZEPS, ZWORK3(JI), ZLV(JI), ZLS(JI), ZCPH(JI) ) + END DO ! ZWORK3(:) = ZWORK3(:) * XRHDBC ZWORK1(:) = MAX( 0., ZWORK3(:) - PDRW(:,JK) ) @@ -434,4 +441,7 @@ SUBROUTINE CONVECT_DOWNDRAFT( KLON, KLEV, & END DO ! IF (LHOOK) CALL DR_HOOK('CONVECT_DOWNDRAFT',1,ZHOOK_HANDLE) +CONTAINS +INCLUDE "convect_satmixratio.h" +! END SUBROUTINE CONVECT_DOWNDRAFT diff --git a/conv/convect_mixing_funct.F90 b/conv/convect_mixing_funct.F90 index 8c9cec2e..51661eee 100644 --- a/conv/convect_mixing_funct.F90 +++ b/conv/convect_mixing_funct.F90 @@ -66,12 +66,12 @@ SUBROUTINE CONVECT_MIXING_FUNCT( D, & ! !* 0.2 Declarations of local variables : ! -REAL :: ZSIGMA = 0.166666667 ! standard deviation -REAL :: ZFE = 4.931813949 ! integral normalization -REAL :: ZSQRTP = 2.506628, ZP = 0.33267 ! constants -REAL :: ZA1 = 0.4361836, ZA2 =-0.1201676 ! constants -REAL :: ZA3 = 0.9372980, ZT1 = 0.500498 ! constants -REAL :: ZE45 = 0.01111 ! constant +REAL, PARAMETER :: ZSIGMA = 0.166666667 ! standard deviation +REAL, PARAMETER :: ZFE = 4.931813949 ! integral normalization +REAL, PARAMETER :: ZSQRTP = 2.506628, ZP = 0.33267 ! constants +REAL, PARAMETER :: ZA1 = 0.4361836, ZA2 =-0.1201676 ! constants +REAL, PARAMETER :: ZA3 = 0.9372980, ZT1 = 0.500498 ! constants +REAL, PARAMETER :: ZE45 = 0.01111 ! constant ! REAL, DIMENSION(D%NIT) :: ZX, ZY, ZW1, ZW2 ! work variables REAL :: ZW11 @@ -87,12 +87,14 @@ SUBROUTINE CONVECT_MIXING_FUNCT( D, & IF (LHOOK) CALL DR_HOOK('CONVECT_MIXING_FUNCT',0,ZHOOK_HANDLE) IF( KMF == 1 ) THEN ! ZX(:) = ( PMIXC(:) - 0.5 ) / ZSIGMA - ZX(D%NIB:D%NIE) = 6. * PMIXC(D%NIB:D%NIE) - 3. - ZW1(D%NIB:D%NIE) = 1. / ( 1.+ ZP * ABS ( ZX(D%NIB:D%NIE) ) ) - ZY(D%NIB:D%NIE) = EXP( -0.5 * ZX(D%NIB:D%NIE) * ZX(D%NIB:D%NIE) ) - ZW2(D%NIB:D%NIE) = ZA1 * ZW1(D%NIB:D%NIE) + ZA2 * ZW1(D%NIB:D%NIE) * ZW1(D%NIB:D%NIE) + & - ZA3 * ZW1(D%NIB:D%NIE) * ZW1(D%NIB:D%NIE) * ZW1(D%NIB:D%NIE) + DO JI = D%NIB, D%NIE + ZX(JI) = 6. * PMIXC(JI) - 3. + ZW1(JI) = 1. / ( 1.+ ZP * ABS ( ZX(JI) ) ) + ZY(JI) = EXP( -0.5 * ZX(JI) * ZX(JI) ) + ZW2(JI) = ZA1 * ZW1(JI) + ZA2 * ZW1(JI) * ZW1(JI) + & + ZA3 * ZW1(JI) * ZW1(JI) * ZW1(JI) ZW11 = ZA1 * ZT1 + ZA2 * ZT1 * ZT1 + ZA3 * ZT1 * ZT1 * ZT1 + ENDDO ENDIF ! DO JI=D%NIB, D%NIE @@ -117,8 +119,10 @@ SUBROUTINE CONVECT_MIXING_FUNCT( D, & ENDDO ! - PER(D%NIB:D%NIE) = PER(D%NIB:D%NIE) * ZFE - PDR(D%NIB:D%NIE) = PDR(D%NIB:D%NIE) * ZFE +DO JI = D%NIB, D%NIE + PER(JI) = PER(JI) * ZFE + PDR(JI) = PDR(JI) * ZFE +ENDDO ! ! ! 2. Use triangular function KMF=2 diff --git a/conv/convect_satmixratio.F90 b/conv/convect_satmixratio.F90 index 53c7ff88..3889ac16 100644 --- a/conv/convect_satmixratio.F90 +++ b/conv/convect_satmixratio.F90 @@ -1,94 +1,3 @@ -! ######spl - SUBROUTINE CONVECT_SATMIXRATIO(CST, D, PPRES, PT, PEW, PLV, PLS, PCPH) - USE YOMHOOK , ONLY : LHOOK, DR_HOOK, JPHOOK -! ################################################################ -! -!!**** Compute vapor saturation mixing ratio over liquid water -!! -!! -!! PDRPOSE -!! ------- -!! The purpose of this routine is to determine saturation mixing ratio -!! and to return values for L_v L_s and C_ph -!! -!! -!!** METHOD -!! ------ -!! -!! -!! EXTERNAL -!! -------- -!! None -!! -!! -!! IMPLICIT ARGUMENTS -!! ------------------ -!! Module MODD_CST -!! XALPW, XBETAW, XGAMW ! constants for water saturation pressure -!! XRD, XRV ! gaz constants for dry air and water vapor -!! XCPD, XCPV ! specific heat for dry air and water vapor -!! XCL, XCI ! specific heat for liquid water and ice -!! XTT ! triple point temperature -!! XLVTT, XLSTT ! vaporization, sublimation heat constant -!! -!! -!! REFERENCE -!! --------- -!! -!! Book1,2 of documentation ( routine CONVECT_SATMIXRATIO) -!! -!! AUTHOR -!! ------ -!! P. BECHTOLD * Laboratoire d'Aerologie * -!! -!! MODIFICATIONS -!! ------------- -!! Original 07/11/95 -!! Last modified 04/10/97 -!------------------------- ------------------------------------------------------ -! -!* 0. DECLARATIONS -! ------------ -! -USE MODD_CST, ONLY : CST_T -USE MODD_DIMPHYEX, ONLY: DIMPHYEX_T -! -! -IMPLICIT NONE -! -!* 0.1 Declarations of dummy arguments : -! -! -TYPE(CST_T), INTENT(IN) :: CST -TYPE(DIMPHYEX_T), INTENT(IN) :: D -REAL, DIMENSION(D%NIT), INTENT(IN) :: PPRES ! pressure -REAL, DIMENSION(D%NIT), INTENT(IN) :: PT ! temperature -! -REAL, DIMENSION(D%NIT), INTENT(OUT):: PEW ! vapor saturation mixing ratio -REAL, DIMENSION(D%NIT), INTENT(OUT):: PLV ! latent heat L_v -REAL, DIMENSION(D%NIT), INTENT(OUT):: PLS ! latent heat L_s -REAL, DIMENSION(D%NIT), INTENT(OUT):: PCPH ! specific heat C_ph -! -!* 0.2 Declarations of local variables : -! -REAL, DIMENSION(D%NIT) :: ZT ! temperature -REAL :: ZEPS ! R_d / R_v -! -! -!------------------------------------------------------------------------------- -! - REAL(KIND=JPHOOK) :: ZHOOK_HANDLE - IF (LHOOK) CALL DR_HOOK('CONVECT_SATMIXRATIO',0,ZHOOK_HANDLE) - ZEPS = CST%XRD / CST%XRV -! - ZT(D%NIB:D%NIE) = MIN( 400., MAX( PT(D%NIB:D%NIE), 10. ) ) ! overflow bound - PEW(D%NIB:D%NIE) = EXP( CST%XALPW - CST%XBETAW / ZT(D%NIB:D%NIE) - CST%XGAMW * ALOG( ZT(D%NIB:D%NIE) ) ) - PEW(D%NIB:D%NIE) = ZEPS * PEW(D%NIB:D%NIE) / ( PPRES(D%NIB:D%NIE) - PEW(D%NIB:D%NIE) ) -! - PLV(D%NIB:D%NIE) = CST%XLVTT + ( CST%XCPV - CST%XCL ) * ( ZT(D%NIB:D%NIE) - CST%XTT ) ! compute L_v - PLS(D%NIB:D%NIE) = CST%XLSTT + ( CST%XCPV - CST%XCI ) * ( ZT(D%NIB:D%NIE) - CST%XTT ) ! compute L_i -! - PCPH(D%NIB:D%NIE) = CST%XCPD + CST%XCPV * PEW(D%NIB:D%NIE) ! compute C_ph -! -IF (LHOOK) CALL DR_HOOK('CONVECT_SATMIXRATIO',1,ZHOOK_HANDLE) +SUBROUTINE CONVECT_SATMIXRATIO + ! dead code, replaced by include file END SUBROUTINE CONVECT_SATMIXRATIO diff --git a/conv/convect_satmixratio.h b/conv/convect_satmixratio.h index 72263a7a..444b0a71 100644 --- a/conv/convect_satmixratio.h +++ b/conv/convect_satmixratio.h @@ -1,17 +1,86 @@ -INTERFACE + ELEMENTAL SUBROUTINE CONVECT_SATMIXRATIO(PPRES, PT, PEPS, PEW, PLV, PLS, PCPH) -SUBROUTINE CONVECT_SATMIXRATIO(CST, D, PPRES, PT, PEW, PLV, PLS, PCPH) -USE YOMHOOK , ONLY : LHOOK, DR_HOOK -USE MODD_CST, ONLY : CST_T -USE MODD_DIMPHYEX, ONLY: DIMPHYEX_T -TYPE(CST_T), INTENT(IN) :: CST -TYPE(DIMPHYEX_T), INTENT(IN) :: D -REAL, DIMENSION(D%NIT), INTENT(IN) :: PPRES -REAL, DIMENSION(D%NIT), INTENT(IN) :: PT -REAL, DIMENSION(D%NIT), INTENT(OUT):: PEW -REAL, DIMENSION(D%NIT), INTENT(OUT):: PLV -REAL, DIMENSION(D%NIT), INTENT(OUT):: PLS -REAL, DIMENSION(D%NIT), INTENT(OUT):: PCPH -END SUBROUTINE CONVECT_SATMIXRATIO +! ******* TO BE INCLUDED IN THE *CONTAINS* OF A SUBROUTINE, IN ORDER TO EASE AUTOMATIC INLINING ****** +! => Don't use drHook !!! -END INTERFACE +! ################################################################ +! +!!**** Compute vapor saturation mixing ratio over liquid water +!! +!! +!! PDRPOSE +!! ------- +!! The purpose of this routine is to determine saturation mixing ratio +!! and to return values for L_v L_s and C_ph +!! +!! +!!** METHOD +!! ------ +!! +!! +!! EXTERNAL +!! -------- +!! None +!! +!! +!! IMPLICIT ARGUMENTS +!! ------------------ +!! Module MODD_CST +!! XALPW, XBETAW, XGAMW ! constants for water saturation pressure +!! XRD, XRV ! gaz constants for dry air and water vapor +!! XCPD, XCPV ! specific heat for dry air and water vapor +!! XCL, XCI ! specific heat for liquid water and ice +!! XTT ! triple point temperature +!! XLVTT, XLSTT ! vaporization, sublimation heat constant +!! +!! +!! REFERENCE +!! --------- +!! +!! Book1,2 of documentation ( routine CONVECT_SATMIXRATIO) +!! +!! AUTHOR +!! ------ +!! P. BECHTOLD * Laboratoire d'Aerologie * +!! +!! MODIFICATIONS +!! ------------- +!! Original 07/11/95 +!! Last modified 04/10/97 +!! R. El Khatib 01-Jun-2023 written as a include file +!------------------------- ------------------------------------------------------ +! +!* 0. DECLARATIONS +! ------------ +! +! +IMPLICIT NONE +! +!* 0.1 Declarations of dummy arguments : +! +! +REAL, INTENT(IN) :: PPRES ! pressure +REAL, INTENT(IN) :: PT ! temperature +REAL, INTENT(IN) :: PEPS ! CST%XRD / CST%XRV (ideally pre-computed in CST%) +! +REAL, INTENT(OUT):: PEW ! vapor saturation mixing ratio +REAL, INTENT(OUT):: PLV ! latent heat L_v +REAL, INTENT(OUT):: PLS ! latent heat L_s +REAL, INTENT(OUT):: PCPH ! specific heat C_ph +! +!* 0.2 Declarations of local variables : +! +REAL :: ZT ! temperature +! +!------------------------------------------------------------------------------- +! + ZT = MIN( 400., MAX( PT, 10. ) ) ! overflow bound + PEW = EXP( CST%XALPW - CST%XBETAW / ZT - CST%XGAMW * ALOG( ZT ) ) + PEW = PEPS * PEW / ( PPRES - PEW ) +! + PLV = CST%XLVTT + ( CST%XCPV - CST%XCL ) * ( ZT - CST%XTT ) ! compute L_v + PLS = CST%XLSTT + ( CST%XCPV - CST%XCI ) * ( ZT - CST%XTT ) ! compute L_i +! + PCPH = CST%XCPD + CST%XCPV * PEW ! compute C_ph +! +END SUBROUTINE CONVECT_SATMIXRATIO diff --git a/conv/convect_trigger_funct.F90 b/conv/convect_trigger_funct.F90 index 8b408bc8..5e28fa43 100644 --- a/conv/convect_trigger_funct.F90 +++ b/conv/convect_trigger_funct.F90 @@ -104,6 +104,7 @@ SUBROUTINE CONVECT_TRIGGER_FUNCT( KLON, KLEV, & INTEGER, DIMENSION(KLON), INTENT(INOUT):: KDPL ! contains vert. index of DPL INTEGER, DIMENSION(KLON), INTENT(INOUT):: KPBL ! contains index of source layer top REAL, DIMENSION(KLON), INTENT(OUT):: PCAPE ! CAPE (J/kg) for diagnostics + ! !* 0.2 Declarations of local variables : ! @@ -132,13 +133,13 @@ SUBROUTINE CONVECT_TRIGGER_FUNCT( KLON, KLEV, & LOGICAL, DIMENSION(KLON) :: GTRIG, GTRIG2 ! local arrays for OTRIG LOGICAL, DIMENSION(KLON) :: GWORK1 ! work array ! +REAL(KIND=JPHOOK) :: ZHOOK_HANDLE ! !------------------------------------------------------------------------------- ! !* 0.3 Compute array bounds ! -------------------- ! -REAL(KIND=JPHOOK) :: ZHOOK_HANDLE IF (LHOOK) CALL DR_HOOK('CONVECT_TRIGGER_FUNCT',0,ZHOOK_HANDLE) IIE = KLON IKB = 1 + JCVEXB @@ -252,7 +253,9 @@ SUBROUTINE CONVECT_TRIGGER_FUNCT( KLON, KLEV, & ! with MNH saturation formula ! --------------------------------------------- ! - CALL CONVECT_SATMIXRATIO( KLON, ZPLCL, ZTLCL, ZWORK1, ZLV, ZWORK2, ZCPH ) + DO JI = 1, IIE + CALL CONVECT_SATMIXRATIO( ZPLCL(JI), ZTLCL(JI), ZEPS, ZWORK1(JI), ZLV(JI), ZWORK2(JI), ZCPH(JI) ) + END DO WHERE( GWORK1(:) ) ZWORK2(:) = ZWORK1(:) / ZTLCL(:) * ( XBETAW / ZTLCL(:) - XGAMW ) ! dr_sat/dT ZWORK2(:) = ( ZWORK1(:) - ZRVLCL(:) ) / & @@ -266,7 +269,9 @@ SUBROUTINE CONVECT_TRIGGER_FUNCT( KLON, KLEV, & ! and temperature to saturation values. ! --------------------------------------------- ! - CALL CONVECT_SATMIXRATIO( KLON, ZPRESMIX, ZTMIX, ZWORK1, ZLV, ZWORK2, ZCPH ) + DO JI = 1, IIE + CALL CONVECT_SATMIXRATIO( ZPRESMIX(JI), ZTMIX(JI), ZEPS, ZWORK1(JI), ZLV(JI), ZWORK2(JI), ZCPH(JI) ) + END DO WHERE( GWORK1(:) .AND. ZRVLCL(:) > ZWORK1(:) ) ZWORK2(:) = ZWORK1(:) / ZTMIX(:) * ( XBETAW / ZTMIX(:) - XGAMW ) ! dr_sat/dT ZWORK2(:) = ( ZWORK1(:) - ZRVLCL(:) ) / & @@ -404,4 +409,7 @@ SUBROUTINE CONVECT_TRIGGER_FUNCT( KLON, KLEV, & ! ! IF (LHOOK) CALL DR_HOOK('CONVECT_TRIGGER_FUNCT',1,ZHOOK_HANDLE) +CONTAINS +INCLUDE "convect_satmixratio.h" +! END SUBROUTINE CONVECT_TRIGGER_FUNCT diff --git a/conv/convect_trigger_shal.F90 b/conv/convect_trigger_shal.F90 index 8be34936..58e91560 100644 --- a/conv/convect_trigger_shal.F90 +++ b/conv/convect_trigger_shal.F90 @@ -89,28 +89,28 @@ SUBROUTINE CONVECT_TRIGGER_SHAL( CVP_SHAL, CVPEXT, CST, D, & ! !* 0.1 Declarations of dummy arguments : ! -TYPE(CONVPAR_SHAL) ,INTENT(IN) :: CVP_SHAL -TYPE(CONVPAREXT) ,INTENT(IN) :: CVPEXT -TYPE(CST_T) ,INTENT(IN) :: CST -TYPE(DIMPHYEX_T) ,INTENT(IN) :: D -REAL ,DIMENSION(D%NIT,D%NKT) ,INTENT(IN) :: PPRES ! pressure -REAL ,DIMENSION(D%NIT,D%NKT) ,INTENT(IN) :: PTH,PTHV ! theta, theta_v -REAL ,DIMENSION(D%NIT,D%NKT) ,INTENT(IN) :: PTHES ! envir. satur. theta_e -REAL ,DIMENSION(D%NIT,D%NKT) ,INTENT(IN) :: PRV ! vapor mixing ratio -REAL ,DIMENSION(D%NIT,D%NKT) ,INTENT(IN) :: PW ! vertical velocity -REAL ,DIMENSION(D%NIT,D%NKT) ,INTENT(IN) :: PZ ! height of grid point (m) -REAL ,DIMENSION(D%NIT) ,INTENT(IN) :: PTKECLS ! TKE CLS -! -REAL ,DIMENSION(D%NIT) ,INTENT(OUT) :: PTHLCL ! theta at LCL -REAL ,DIMENSION(D%NIT) ,INTENT(OUT) :: PTLCL ! temp. at LCL -REAL ,DIMENSION(D%NIT) ,INTENT(OUT) :: PRVLCL ! vapor mixing ratio at LCL -REAL ,DIMENSION(D%NIT) ,INTENT(OUT) :: PWLCL ! parcel velocity at LCL -REAL ,DIMENSION(D%NIT) ,INTENT(OUT) :: PZLCL ! height at LCL (m) -REAL ,DIMENSION(D%NIT) ,INTENT(OUT) :: PTHVELCL ! environm. theta_v at LCL (K) -INTEGER ,DIMENSION(D%NIT) ,INTENT(INOUT) :: KLCL ! contains vert. index of LCL -INTEGER ,DIMENSION(D%NIT) ,INTENT(INOUT) :: KDPL ! contains vert. index of DPL -INTEGER ,DIMENSION(D%NIT) ,INTENT(INOUT) :: KPBL ! contains index of source layer top -LOGICAL ,DIMENSION(D%NIT) ,INTENT(OUT) :: OTRIG ! logical mask for convection +TYPE(CONVPAR_SHAL), INTENT(IN) :: CVP_SHAL +TYPE(CONVPAREXT), INTENT(IN) :: CVPEXT +TYPE(CST_T), INTENT(IN) :: CST +TYPE(DIMPHYEX_T), INTENT(IN) :: D +REAL, DIMENSION(D%NIT), INTENT(IN) :: PTKECLS ! TKE CLS +REAL, DIMENSION(D%NIT,D%NKT),INTENT(IN) :: PTH, PTHV ! theta, theta_v +REAL, DIMENSION(D%NIT,D%NKT),INTENT(IN) :: PTHES ! envir. satur. theta_e +REAL, DIMENSION(D%NIT,D%NKT),INTENT(IN) :: PRV ! vapor mixing ratio +REAL, DIMENSION(D%NIT,D%NKT),INTENT(IN) :: PPRES ! pressure +REAL, DIMENSION(D%NIT,D%NKT),INTENT(IN) :: PZ ! height of grid point (m) +REAL, DIMENSION(D%NIT,D%NKT),INTENT(IN) :: PW ! vertical velocity +! +REAL, DIMENSION(D%NIT), INTENT(OUT):: PTHLCL ! theta at LCL +REAL, DIMENSION(D%NIT), INTENT(OUT):: PTLCL ! temp. at LCL +REAL, DIMENSION(D%NIT), INTENT(OUT):: PRVLCL ! vapor mixing ratio at LCL +REAL, DIMENSION(D%NIT), INTENT(OUT):: PWLCL ! parcel velocity at LCL +REAL, DIMENSION(D%NIT), INTENT(OUT):: PZLCL ! height at LCL (m) +REAL, DIMENSION(D%NIT), INTENT(OUT):: PTHVELCL ! environm. theta_v at LCL (K) +LOGICAL, DIMENSION(D%NIT), INTENT(OUT):: OTRIG ! logical mask for convection +INTEGER, DIMENSION(D%NIT), INTENT(INOUT):: KLCL ! contains vert. index of LCL +INTEGER, DIMENSION(D%NIT), INTENT(INOUT):: KDPL ! contains vert. index of DPL +INTEGER, DIMENSION(D%NIT), INTENT(INOUT):: KPBL ! contains index of source layer top ! !* 0.2 Declarations of local variables : ! @@ -119,6 +119,7 @@ SUBROUTINE CONVECT_TRIGGER_SHAL( CVP_SHAL, CVPEXT, CST, D, & INTEGER :: IKB, IKE ! horizontal + vertical loop bounds REAL :: ZEPS, ZEPSA ! R_d / R_v, R_v / R_d REAL :: ZCPORD, ZRDOCP ! C_pd / R_d, R_d / C_pd +REAL :: ZX1, ZX2 ! REAL, DIMENSION(D%NIT) :: ZTHLCL, ZTLCL, ZRVLCL, & ! locals for PTHLCL,PTLCL ZWLCL, ZZLCL, ZTHVELCL ! PRVLCL, .... @@ -132,22 +133,23 @@ SUBROUTINE CONVECT_TRIGGER_SHAL( CVP_SHAL, CVPEXT, CST, D, & REAL, DIMENSION(D%NIT) :: ZCAPE ! convective available energy (m^2/s^2/g) REAL, DIMENSION(D%NIT) :: ZCAP ! pseudo fro CAPE REAL, DIMENSION(D%NIT) :: ZTHEUL ! updraft equiv. pot. temperature (K) -REAL, DIMENSION(D%NIT) :: ZLV, ZCPH! specific heats of vaporisation, dry air +REAL, DIMENSION(D%NIT) :: ZLVA, ZCPHA! specific heats of vaporisation, dry air +REAL, DIMENSION(D%NIT) :: ZLVB, ZCPHB! specific heats of vaporisation, dry air +REAL, DIMENSION(D%NIT) :: ZEWA, ZEWB ! vapor saturation mixing ratios +REAL, DIMENSION(D%NIT) :: ZLSA, ZLSB ! latent heat L_s REAL, DIMENSION(D%NIT) :: ZDP ! pressure between LCL and model layer REAL, DIMENSION(D%NIT) :: ZTOP,ZTOPP ! estimated cloud top (m) REAL, DIMENSION(D%NIT) :: ZWORK1, ZWORK2, ZWORK3 ! work arrays LOGICAL, DIMENSION(D%NIT) :: GTRIG2 ! local arrays for OTRIG LOGICAL, DIMENSION(D%NIT) :: GWORK1 ! work array ! +REAL(KIND=JPHOOK) :: ZHOOK_HANDLE ! !------------------------------------------------------------------------------- ! !* 0.3 Compute array bounds ! -------------------- ! -REAL(KIND=JPHOOK) :: ZHOOK_HANDLE - -#include "convect_satmixratio.h" IF (LHOOK) CALL DR_HOOK('CONVECT_TRIGGER_SHAL',0,ZHOOK_HANDLE) IKB = 1 + CVPEXT%JCVEXB @@ -188,7 +190,9 @@ SUBROUTINE CONVECT_TRIGGER_SHAL( CVP_SHAL, CVPEXT, CST, D, & ! DO JKK = IKB + 1, IKE - 2 ! - GWORK1(D%NIB:D%NIE) = ZZDPL(D%NIB:D%NIE) - PZ(D%NIB:D%NIE,IKB) < CVP_SHAL%XZLCL + DO JI=D%NIB, D%NIE + GWORK1(JI) = ZZDPL(JI) - PZ(JI,IKB) < CVP_SHAL%XZLCL + ENDDO ! we exit the trigger test when the center of the mixed layer is more ! than 1500 m above soil level. DO JI=D%NIB, D%NIE @@ -211,11 +215,11 @@ SUBROUTINE CONVECT_TRIGGER_SHAL( CVP_SHAL, CVPEXT, CST, D, & DO JI = D%NIB, D%NIE IF ( GWORK1(JI) .AND. ZDPTHMIX(JI) < CVP_SHAL%XZPBL ) THEN IPBL(JI) = JK - ZWORK1(JI) = PPRES(JI,JK) - PPRES(JI,JKM) - ZDPTHMIX(JI) = ZDPTHMIX(JI) + ZWORK1(JI) - ZPRESMIX(JI) = ZPRESMIX(JI) + PPRES(JI,JK) * ZWORK1(JI) - ZTHLCL(JI) = ZTHLCL(JI) + PTH(JI,JK) * ZWORK1(JI) - ZRVLCL(JI) = ZRVLCL(JI) + MAX(0., PRV(JI,JK)) * ZWORK1(JI) + ZX1 = PPRES(JI,JK) - PPRES(JI,JKM) + ZDPTHMIX(JI) = ZDPTHMIX(JI) + ZX1 + ZPRESMIX(JI) = ZPRESMIX(JI) + PPRES(JI,JK) * ZX1 + ZTHLCL(JI) = ZTHLCL(JI) + PTH(JI,JK) * ZX1 + ZRVLCL(JI) = ZRVLCL(JI) + MAX(0., PRV(JI,JK)) * ZX1 END IF END DO END DO @@ -241,30 +245,33 @@ SUBROUTINE CONVECT_TRIGGER_SHAL( CVP_SHAL, CVPEXT, CST, D, & ZTMIX(JI) = ZTHLCL(JI) * ( ZPRESMIX(JI) / CST%XP00 ) ** ZRDOCP ZEVMIX(JI) = ZRVLCL(JI) * ZPRESMIX(JI) / ( ZRVLCL(JI) + ZEPS ) ZEVMIX(JI) = MAX( 1.E-8, ZEVMIX(JI) ) - ZWORK1(JI) = LOG( ZEVMIX(JI) / 613.3 ) + ZX1 = LOG( ZEVMIX(JI) / 613.3 ) ! dewpoint temperature - ZWORK1(JI) = ( 4780.8 - 32.19 * ZWORK1(JI) ) / ( 17.502 - ZWORK1(JI) ) + ZX1 = ( 4780.8 - 32.19 * ZX1 ) / ( 17.502 - ZX1 ) ! adiabatic saturation temperature - ZTLCL(JI) = ZWORK1(JI) - ( .212 + 1.571E-3 * ( ZWORK1(JI) - CST%XTT ) & - - 4.36E-4 * ( ZTMIX(JI) - CST%XTT ) ) * ( ZTMIX(JI) - ZWORK1(JI) ) + ZTLCL(JI) = ZX1 - ( .212 + 1.571E-3 * ( ZX1 - CST%XTT ) & + - 4.36E-4 * ( ZTMIX(JI) - CST%XTT ) ) * ( ZTMIX(JI) - ZX1 ) ZTLCL(JI) = MIN( ZTLCL(JI), ZTMIX(JI) ) ZPLCL(JI) = CST%XP00 * ( ZTLCL(JI) / ZTHLCL(JI) ) ** ZCPORD ! END IF ENDDO ! + DO JI=D%NIB, D%NIE + CALL CONVECT_SATMIXRATIO( ZPLCL(JI), ZTLCL(JI), ZEPS, ZEWA(JI), ZLVA(JI), ZLSA(JI), ZCPHA(JI) ) + CALL CONVECT_SATMIXRATIO( ZPRESMIX(JI), ZTMIX(JI), ZEPS, ZEWB(JI), ZLVB(JI), ZLSB(JI), ZCPHB(JI) ) + ENDDO ! !* 4.2 Correct ZTLCL in order to be completely consistent ! with MNH saturation formula ! --------------------------------------------- ! - CALL CONVECT_SATMIXRATIO( CST, D, ZPLCL, ZTLCL, ZWORK1, ZLV, ZWORK2, ZCPH ) DO JI=D%NIB, D%NIE IF( GWORK1(JI) ) THEN - ZWORK2(JI) = ZWORK1(JI) / ZTLCL(JI) * ( CST%XBETAW / ZTLCL(JI) - CST%XGAMW ) ! dr_sat/dT - ZWORK2(JI) = ( ZWORK1(JI) - ZRVLCL(JI) ) / & - ( 1. + ZLV(JI) / ZCPH(JI) * ZWORK2(JI) ) - ZTLCL(JI) = ZTLCL(JI) - ZLV(JI) / ZCPH(JI) * ZWORK2(JI) + ZLSA(JI) = ZEWA(JI) / ZTLCL(JI) * ( CST%XBETAW / ZTLCL(JI) - CST%XGAMW ) ! dr_sat/dT + ZLSA(JI) = ( ZEWA(JI) - ZRVLCL(JI) ) / & + ( 1. + ZLVA(JI) / ZCPHA(JI) * ZLSA(JI) ) + ZTLCL(JI) = ZTLCL(JI) - ZLVA(JI) / ZCPHA(JI) * ZLSA(JI) ! END IF ENDDO @@ -274,14 +281,13 @@ SUBROUTINE CONVECT_TRIGGER_SHAL( CVP_SHAL, CVPEXT, CST, D, & ! and temperature to saturation values. ! --------------------------------------------- ! - CALL CONVECT_SATMIXRATIO( CST, D, ZPRESMIX, ZTMIX, ZWORK1, ZLV, ZWORK2, ZCPH ) DO JI=D%NIB, D%NIE - IF( GWORK1(JI) .AND. ZRVLCL(JI) > ZWORK1(JI) ) THEN - ZWORK2(JI) = ZWORK1(JI) / ZTMIX(JI) * ( CST%XBETAW / ZTMIX(JI) - CST%XGAMW ) ! dr_sat/dT - ZWORK2(JI) = ( ZWORK1(JI) - ZRVLCL(JI) ) / & - ( 1. + ZLV(JI) / ZCPH(JI) * ZWORK2(JI) ) - ZTLCL(JI) = ZTMIX(JI) - ZLV(JI) / ZCPH(JI) * ZWORK2(JI) - ZRVLCL(JI) = ZRVLCL(JI) - ZWORK2(JI) + IF( GWORK1(JI) .AND. ZRVLCL(JI) > ZEWB(JI) ) THEN + ZLSB(JI) = ZEWB(JI) / ZTMIX(JI) * ( CST%XBETAW / ZTMIX(JI) - CST%XGAMW ) ! dr_sat/dT + ZLSB(JI) = ( ZEWB(JI) - ZRVLCL(JI) ) / & + ( 1. + ZLVB(JI) / ZCPHB(JI) * ZLSB(JI) ) + ZTLCL(JI) = ZTMIX(JI) - ZLVB(JI) / ZCPHB(JI) * ZLSB(JI) + ZRVLCL(JI) = ZRVLCL(JI) - ZLSB(JI) ZPLCL(JI) = ZPRESMIX(JI) ZTHLCL(JI) = ZTLCL(JI) * ( CST%XP00 / ZPLCL(JI) ) ** ZRDOCP ZTHVLCL(JI)= ZTHLCL(JI) * ( 1. + ZEPSA * ZRVLCL(JI) ) & @@ -353,7 +359,8 @@ SUBROUTINE CONVECT_TRIGGER_SHAL( CVP_SHAL, CVPEXT, CST, D, & ! GTRIG(:) = ZTHVLCL(:) - ZTHVELCL(:) + ZWORK1(:) > 0. .AND. & ! ZWLCL(:) > 0. ! END WHERE - ZWLCL(D%NIB:D%NIE) = CVP_SHAL%XAW * MAX(0.,PW(D%NIB:D%NIE,IKB)) + CVP_SHAL%XBW + DO JI = D%NIB, D%NIE + ZWLCL(JI) = CVP_SHAL%XAW * MAX(0.,PW(JI,IKB)) + CVP_SHAL%XBW ! ! !* 6.3 Look for parcel that produces sufficient cloud depth. @@ -361,11 +368,15 @@ SUBROUTINE CONVECT_TRIGGER_SHAL( CVP_SHAL, CVPEXT, CST, D, & ! is smaller than a given value (based on vertical velocity eq.) ! -------------------------------------------------------------- ! - ZTHEUL(D%NIB:D%NIE) = ZTLCL(D%NIB:D%NIE) * ( ZTHLCL(D%NIB:D%NIE) / ZTLCL(D%NIB:D%NIE) ) & - ** ( 1. - 0.28 * ZRVLCL(D%NIB:D%NIE) ) & - * EXP( ( 3374.6525 / ZTLCL(D%NIB:D%NIE) - 2.5403 ) * & - ZRVLCL(D%NIB:D%NIE) * ( 1. + 0.81 * ZRVLCL(D%NIB:D%NIE) ) ) + ZTHEUL(JI) = ZTLCL(JI) * ( ZTHLCL(JI) / ZTLCL(JI) ) & + ** ( 1. - 0.28 * ZRVLCL(JI) ) & + * EXP( ( 3374.6525 / ZTLCL(JI) - 2.5403 ) * & + ZRVLCL(JI) * ( 1. + 0.81 * ZRVLCL(JI) ) ) + END DO + ! +!$acc enter data create(zwork3(1:144116), ztopp(1:144116), zcap(1:144116),zcape(1:144116),ztop(1:144116)) +!$acc kernels present(zwork3(1:144116), ztopp(1:144116), zcap(1:144116), zcape(1:144116),ztop(1:144116)) ZCAPE(D%NIB:D%NIE) = 0. ZCAP(D%NIB:D%NIE) = 0. ZTOP(D%NIB:D%NIE) = 0. @@ -374,33 +385,33 @@ SUBROUTINE CONVECT_TRIGGER_SHAL( CVP_SHAL, CVPEXT, CST, D, & JKM = IKB DO JL = JKM, JT JK = JL + 1 +!REK : should try if's instead of sign & max DO JI = D%NIB, D%NIE - ZWORK1(JI) = ( 2. * ZTHEUL(JI) / & + ZX1 = ( 2. * ZTHEUL(JI) / & ( PTHES(JI,JK) + PTHES(JI,JL) ) - 1. ) * ( PZ(JI,JK) - PZ(JI,JL) ) - IF ( JL < ILCL(JI) ) ZWORK1(JI) = 0. - ZCAPE(JI) = ZCAPE(JI) + CST%XG * MAX( 1., ZWORK1(JI) ) - ZCAP(JI) = ZCAP(JI) + ZWORK1(JI) - ZWORK2(JI) = CVP_SHAL%XNHGAM * CST%XG * ZCAP(JI) + 1.05 * ZWLCL(JI) * ZWLCL(JI) + IF ( JL < ILCL(JI) ) ZX1 = 0. + ZCAPE(JI) = ZCAPE(JI) + CST%XG * MAX( 1., ZX1 ) + ZCAP(JI) = ZCAP(JI) + ZX1 ! the factor 1.05 takes entrainment into account - ZWORK2(JI) = SIGN( 1., ZWORK2(JI) ) - ZWORK3(JI) = ZWORK3(JI) + MIN(0., ZWORK2(JI) ) - ZWORK3(JI) = MAX( -1., ZWORK3(JI) ) - ! Nota, the factors ZWORK2 and ZWORK3 are only used to avoid + ZX2 = SIGN( 1., ( CVP_SHAL%XNHGAM * CST%XG * ZCAP(JI) + 1.05 * ZWLCL(JI) * ZWLCL(JI)) ) + ZWORK3(JI) = MAX( -1., (ZWORK3(JI) + MIN(0.,ZX2)) ) + ! Nota, the factors ZX2 and ZWORK3 are only used to avoid ! if and goto statements, the difficulty is to extract only ! the level where the criterium is first fullfilled ZTOPP(JI)=ZTOP(JI) - ZTOP(JI) = PZ(JI,JL) * .5 * ( 1. + ZWORK2(JI) ) * ( 1. + ZWORK3(JI) ) + & - ZTOP(JI) * .5 * ( 1. - ZWORK2(JI) ) + ZTOP(JI) = PZ(JI,JL) * .5 * ( 1. + ZX2 ) * ( 1. + ZWORK3(JI) ) + & + ZTOP(JI) * .5 * ( 1. - ZX2 ) ZTOP(JI)=MAX(ZTOP(JI),ZTOPP(JI)) ZTOPP(JI)=ZTOP(JI) END DO END DO + !$acc end kernels + !$acc exit data delete(zwork3, ztopp, zcap, zcape, ztop) ! ! - ZWORK2(D%NIB:D%NIE) = ZTOP(D%NIB:D%NIE) - ZZLCL(D%NIB:D%NIE) - ! WHERE( ZWORK2(:) .GE. XCDEPTH .AND. ZWORK2(:) < XCDEPTH_D .AND. GTRIG2(:) & + DO JI=D%NIB, D%NIE - IF( ZWORK2(JI) .GE. CVP_SHAL%XCDEPTH .AND. GTRIG2(JI) .AND. ZCAPE(JI) > 10. )THEN + IF( (ZTOP(JI)-ZZLCL(JI)) .GE. CVP_SHAL%XCDEPTH .AND. GTRIG2(JI) .AND. ZCAPE(JI) > 10. )THEN GTRIG2(JI) = .FALSE. OTRIG(JI) = .TRUE. ! OTRIG(JI) = GTRIG(JI) ! we select the first departure level @@ -420,5 +431,7 @@ SUBROUTINE CONVECT_TRIGGER_SHAL( CVP_SHAL, CVPEXT, CST, D, & ! ! IF (LHOOK) CALL DR_HOOK('CONVECT_TRIGGER_SHAL',1,ZHOOK_HANDLE) +CONTAINS +INCLUDE "convect_satmixratio.h" +! END SUBROUTINE CONVECT_TRIGGER_SHAL - diff --git a/conv/convect_updraft.F90 b/conv/convect_updraft.F90 index ce59208d..d785c510 100644 --- a/conv/convect_updraft.F90 +++ b/conv/convect_updraft.F90 @@ -79,6 +79,7 @@ SUBROUTINE CONVECT_UPDRAFT( KLON, KLEV, & USE MODD_CST USE MODD_CONVPAR USE MODD_CONVPAREXT +USE MODD_DIMPHYEX ! ! IMPLICIT NONE @@ -160,13 +161,18 @@ SUBROUTINE CONVECT_UPDRAFT( KLON, KLEV, & ! work arrays LOGICAL, DIMENSION(KLON,KLEV) :: GWORK6 ! work array ! -! +TYPE(DIMPHYEX_T) :: D +TYPE(CONVPAR_T) :: CONVPAR + +REAL(KIND=JPHOOK) :: ZHOOK_HANDLE + +#include "convect_condens.h" +#include "convect_mixing_funct.h" !------------------------------------------------------------------------------- ! ! 0.3 Set loop bounds ! --------------- ! -REAL(KIND=JPHOOK) :: ZHOOK_HANDLE IF (LHOOK) CALL DR_HOOK('CONVECT_UPDRAFT',0,ZHOOK_HANDLE) IKB = 1 + JCVEXB IKE = KLEV - JCVEXT @@ -208,6 +214,12 @@ SUBROUTINE CONVECT_UPDRAFT( KLON, KLEV, & GWORK1(:) = .FALSE. GWORK4(:) = .FALSE. ! +CONVPAR%XTFRZ1=XTFRZ1 +CONVPAR%XTFRZ2=XTFRZ2 +D%NIT=KLON +D%NIB=1 +D%NIE=KLON + ! !* 1.1 Compute undilute updraft theta_e for CAPE computations ! Bolton (1980) formula. @@ -262,8 +274,8 @@ SUBROUTINE CONVECT_UPDRAFT( KLON, KLEV, & ! ZWORK1(:) = PURC(:,JK) + PURR(:,JK) ZWORK2(:) = PURI(:,JK) + PURS(:,JK) - CALL CONVECT_CONDENS( KLON, KICE, PPRES(:,JKP), PUTHL(:,JK), PURW(:,JK),& - ZWORK1, ZWORK2, PZ(:,JKP), GWORK1, ZUT, ZURV, & + CALL CONVECT_CONDENS( CST, D, CONVPAR, KICE, PPRES(:,JKP), PUTHL(:,JK), PURW(:,JK),& + ZWORK1, ZWORK2, PZ(:,JKP), ZUT, ZURV, & PURC(:,JKP), PURI(:,JKP), ZLV, ZLS, ZCPH ) ! ! @@ -337,8 +349,8 @@ SUBROUTINE CONVECT_UPDRAFT( KLON, KLEV, & ZWORK2(:) = ZMIXF(:) * PRW(:,JKP) & + ( 1. - ZMIXF(:) ) * PURW(:,JKP) ! mixed r_w ! - CALL CONVECT_CONDENS( KLON, KICE, PPRES(:,JKP), ZWORK1, ZWORK2, & - PURC(:,JKP), PURI(:,JKP), PZ(:,JKP), GWORK1, ZUT,& + CALL CONVECT_CONDENS( CST, D, CONVPAR, KICE, PPRES(:,JKP), ZWORK1, ZWORK2, & + PURC(:,JKP), PURI(:,JKP), PZ(:,JKP), ZUT,& ZWORK3, ZWORK4, ZWORK5, ZLV, ZLS, ZCPH ) ! put in enthalpy and r_w and get T r_c, r_i (ZUT, ZWORK4-5) ! @@ -357,7 +369,8 @@ SUBROUTINE CONVECT_UPDRAFT( KLON, KLEV, & ! ------------------------------------------------------- ! ! - CALL CONVECT_MIXING_FUNCT ( KLON, ZMIXF, 1, ZE2, ZD2 ) +CALL ABOR1('FIXME : THE INTERFACE IS WRONG') + !CALL CONVECT_MIXING_FUNCT ( KLON, ZMIXF, 1, ZE2, ZD2 ) ! Note: routine MIXING_FUNCT returns fractional entrainm/detrainm. rates ! ! ZWORK1(:) = XENTR * PMFLCL(:) * PDPRES(:,JKP) / XCRAD ! rate of env. inflow diff --git a/conv/convect_updraft_shal.F90 b/conv/convect_updraft_shal.F90 index e60e2ec9..6a12ba24 100644 --- a/conv/convect_updraft_shal.F90 +++ b/conv/convect_updraft_shal.F90 @@ -6,6 +6,7 @@ SUBROUTINE CONVECT_UPDRAFT_SHAL( CVP_SHAL, CVPEXT, CST, D, CONVPAR, PUMF, PUER, PUDR, PUTHL, PUTHV, PURW, & PURC, PURI, PCAPE, KCTL, KETL,GTRIG1 ) USE YOMHOOK , ONLY : LHOOK, DR_HOOK, JPHOOK + ! ############################################################################### ! !!**** Compute updraft properties from DPL to CTL. @@ -208,13 +209,16 @@ SUBROUTINE CONVECT_UPDRAFT_SHAL( CVP_SHAL, CVPEXT, CST, D, CONVPAR, ! Define accurate enthalpy for updraft ! ----------------------------------------------------- ! -ZTHEUL(D%NIB:D%NIE) = PTLCL(D%NIB:D%NIE) * ( PTHLCL(D%NIB:D%NIE) / PTLCL(D%NIB:D%NIE) ) ** ( 1. - 0.28 * PRVLCL(D%NIB:D%NIE) ) & - * EXP( ( 3374.6525 / PTLCL(D%NIB:D%NIE) - 2.5403 ) * & - PRVLCL(D%NIB:D%NIE) * ( 1. + 0.81 * PRVLCL(D%NIB:D%NIE) ) ) + +DO JI = D%NIB, D%NIE + ZTHEUL(JI) = PTLCL(JI) * ( PTHLCL(JI) / PTLCL(JI) ) ** ( 1. - 0.28 * PRVLCL(JI) ) & + * EXP( ( 3374.6525 / PTLCL(JI) - 2.5403 ) * & + PRVLCL(JI) * ( 1. + 0.81 * PRVLCL(JI) ) ) ! ! -ZWORK1(D%NIB:D%NIE) = ( CST%XCPD + PRVLCL(D%NIB:D%NIE) * CST%XCPV ) * PTLCL(D%NIB:D%NIE) & - + ( 1. + PRVLCL(D%NIB:D%NIE) ) * CST%XG * PZLCL(D%NIB:D%NIE) + ZWORK1(JI) = ( CST%XCPD + PRVLCL(JI) * CST%XCPV ) * PTLCL(JI) & + + ( 1. + PRVLCL(JI) ) * CST%XG * PZLCL(JI) +ENDDO ! ! !* 2. Set updraft properties between DPL and LCL @@ -242,10 +246,13 @@ SUBROUTINE CONVECT_UPDRAFT_SHAL( CVP_SHAL, CVPEXT, CST, D, CONVPAR, DO JK = IKB + 1, IKE - 1 ZWORK6(:) = 1. JKP = JK + 1 + ! - GWORK4(D%NIB:D%NIE) = JK >= KLCL(D%NIB:D%NIE) - 1 - GWORK1(D%NIB:D%NIE) = GWORK4(D%NIB:D%NIE) .AND. GWORK2(D%NIB:D%NIE) ! this mask is used to confine - ! updraft computations between the LCL and the CTL + DO JI = D%NIB, D%NIE + GWORK4(JI) = JK >= KLCL(JI) - 1 + GWORK1(JI) = GWORK4(JI) .AND. GWORK2(JI) ! this mask is used to confine + ! updraft computations between the LCL and the CTL + ENDDO ! DO JI=D%NIB, D%NIE IF( JK == KLCL(JI) - 1 ) ZWORK6(JI) = 0. ! factor that is used in buoyancy @@ -255,16 +262,19 @@ SUBROUTINE CONVECT_UPDRAFT_SHAL( CVP_SHAL, CVPEXT, CST, D, CONVPAR, !* 4. Estimate condensate, L_v L_i, Cph and theta_v at level k+1 ! ---------------------------------------------------------- ! - ZWORK1(D%NIB:D%NIE) = PURC(D%NIB:D%NIE,JK) - ZWORK2(D%NIB:D%NIE) = PURI(D%NIB:D%NIE,JK) - CALL CONVECT_CONDENS(CST, D, CONVPAR, KICE, PPRES(D%NIB:D%NIE,JKP),& - PUTHL(D%NIB:D%NIE,JK), PURW(D%NIB:D%NIE,JK), & - ZWORK1, ZWORK2, PZ(D%NIB:D%NIE,JKP), ZUT,ZURV,& - PURC(D%NIB:D%NIE,JKP), PURI(D%NIB:D%NIE,JKP), & - ZLV, ZLS, ZCPH ) + ZWORK1(:) = PURC(:,JK) + ZWORK2(:) = PURI(:,JK) + CALL CONVECT_CONDENS(CST, D, CONVPAR, KICE, PPRES(:,JKP),& + PUTHL(:,JK), PURW(:,JK), & + ZWORK1, ZWORK2, PZ(:,JKP), ZUT,ZURV,& + PURC(:,JKP), PURI(:,JKP), & + ZLV, ZLS, ZCPH ) ! ! - ZPI(D%NIB:D%NIE) = ( CST%XP00 / PPRES(D%NIB:D%NIE,JKP) ) ** ZRDOCP + DO JI = D%NIB, D%NIE + ZPI(JI) = ( CST%XP00 / PPRES(JI,JKP) ) ** ZRDOCP + ENDDO + DO JI=D%NIB, D%NIE IF ( GWORK1(JI) ) THEN ! @@ -320,26 +330,30 @@ SUBROUTINE CONVECT_UPDRAFT_SHAL( CVP_SHAL, CVPEXT, CST, D, CONVPAR, ! evaluating the derivative using ZMIXF=0.1. ! ----------------------------------------------------- ! - ZMIXF(D%NIB:D%NIE) = 0.1 ! starting value for critical mixed fraction - ZWORK1(D%NIB:D%NIE) = ZMIXF(D%NIB:D%NIE) * PTHL(D%NIB:D%NIE,JKP) & - + ( 1. - ZMIXF(D%NIB:D%NIE) ) * PUTHL(D%NIB:D%NIE,JKP) ! mixed enthalpy - ZWORK2(D%NIB:D%NIE) = ZMIXF(D%NIB:D%NIE) * PRW(D%NIB:D%NIE,JKP) & - + ( 1. - ZMIXF(D%NIB:D%NIE) ) * PURW(D%NIB:D%NIE,JKP) ! mixed r_w + DO JI = D%NIB, D%NIE + ZMIXF(JI) = 0.1 ! starting value for critical mixed fraction + ZWORK1(JI) = ZMIXF(JI) * PTHL(JI,JKP) & + + ( 1. - ZMIXF(JI) ) * PUTHL(JI,JKP) ! mixed enthalpy + ZWORK2(JI) = ZMIXF(JI) * PRW(JI,JKP) & + + ( 1. - ZMIXF(JI) ) * PURW(JI,JKP) ! mixed r_w + ENDDO ! - CALL CONVECT_CONDENS(CST, D, CONVPAR, KICE, PPRES(D%NIB:D%NIE,JKP),& - ZWORK1, ZWORK2, PURC(D%NIB:D%NIE,JKP), & - PURI(D%NIB:D%NIE,JKP), PZ(D%NIB:D%NIE,JKP), & - ZUT, ZWORK3, ZWORK4, ZWORK5, ZLV, ZLS, ZCPH) + CALL CONVECT_CONDENS(CST, D, CONVPAR, KICE, PPRES(:,JKP),& + ZWORK1, ZWORK2, PURC(:,JKP), & + PURI(:,JKP), PZ(:,JKP), & + ZUT, ZWORK3, ZWORK4, ZWORK5, ZLV, ZLS, ZCPH) ! put in enthalpy and r_w and get T r_c, r_i (ZUT, ZWORK4-5) ! ! compute theta_v of mixture - ZWORK3(D%NIB:D%NIE) = ZUT(D%NIB:D%NIE) * ZPI(D%NIB:D%NIE) * ( 1. + ZEPSA * ( & - ZWORK2(D%NIB:D%NIE) - ZWORK4(D%NIB:D%NIE) - ZWORK5(D%NIB:D%NIE) ) ) / ( 1. + ZWORK2(D%NIB:D%NIE) ) + DO JI = D%NIB, D%NIE + ZWORK3(JI) = ZUT(JI) * ZPI(JI) * ( 1. + ZEPSA * ( & + ZWORK2(JI) - ZWORK4(JI) - ZWORK5(JI) ) ) / ( 1. + ZWORK2(JI) ) ! compute final value of critical mixed fraction using theta_v ! of mixture, grid-scale and updraft - ZMIXF(D%NIB:D%NIE) = MAX( 0., PUTHV(D%NIB:D%NIE,JKP) - PTHV(D%NIB:D%NIE,JKP) ) * ZMIXF(D%NIB:D%NIE) / & - ( PUTHV(D%NIB:D%NIE,JKP) - ZWORK3(D%NIB:D%NIE) + 1.E-10 ) - ZMIXF(D%NIB:D%NIE) = MAX( 0., MIN( 1., ZMIXF(D%NIB:D%NIE) ) ) + ZMIXF(JI) = MAX( 0., PUTHV(JI,JKP) - PTHV(JI,JKP) ) * ZMIXF(JI) / & + ( PUTHV(JI,JKP) - ZWORK3(JI) + 1.E-10 ) + ZMIXF(JI) = MAX( 0., MIN( 1., ZMIXF(JI) ) ) + ENDDO ! ! !* 8.2 Compute final midlevel values for entr. and detrainment @@ -354,7 +368,9 @@ SUBROUTINE CONVECT_UPDRAFT_SHAL( CVP_SHAL, CVPEXT, CST, D, CONVPAR, ! ! ZWORK1(D%NIB:D%NIE) = XENTR * PMFLCL * PDPRES(D%NIB:D%NIE,JKP) / XCRAD ! rate of env. inflow !*MOD - zwork1(D%NIB:D%NIE) = CVP_SHAL%xentr * CST%xg / CVP_SHAL%xcrad * pumf(D%NIB:D%NIE,jk) * ( pz(D%NIB:D%NIE,jkp) - pz(D%NIB:D%NIE,jk) ) + DO JI = D%NIB, D%NIE + zwork1(JI) = CVP_SHAL%xentr * CST%xg / CVP_SHAL%xcrad * pumf(JI,jk) * ( pz(JI,jkp) - pz(JI,jk) ) + ENDDO ! ZWORK1(D%NIB:D%NIE) = XENTR * pumf(D%NIB:D%NIE,jk) * PDPRES(D%NIB:D%NIE,JKP) / XCRAD ! rate of env. inflow !*MOD ZWORK2(:) = 0. @@ -394,7 +410,10 @@ SUBROUTINE CONVECT_UPDRAFT_SHAL( CVP_SHAL, CVPEXT, CST, D, CONVPAR, DO JI=D%NIB, D%NIE IF ( GWORK2(JI) ) KCTL(JI) = JKP ! cloud top level ENDDO - GWORK1(D%NIB:D%NIE) = GWORK2(D%NIB:D%NIE) .AND. GWORK4(D%NIB:D%NIE) + + DO JI=D%NIB, D%NIE + GWORK1(JI) = GWORK2(JI) .AND. GWORK4(JI) + ENDDO ! !IF ( COUNT( GWORK2(:) ) == 0 ) EXIT ! @@ -453,8 +472,12 @@ SUBROUTINE CONVECT_UPDRAFT_SHAL( CVP_SHAL, CVPEXT, CST, D, CONVPAR, DO JI = D%NIB, D%NIE IF( .NOT. OTRIG(JI) ) KCTL(JI) = IKB ENDDO -KETL(D%NIB:D%NIE) = MAX( KETL(D%NIB:D%NIE), KLCL(D%NIB:D%NIE) + 2 ) -KETL(D%NIB:D%NIE) = MIN( KETL(D%NIB:D%NIE), KCTL(D%NIB:D%NIE) ) + + +DO JI = D%NIB, D%NIE + KETL(JI) = MAX( KETL(JI), KLCL(JI) + 2 ) + KETL(JI) = MIN( KETL(JI), KCTL(JI) ) +ENDDO ! ! !* 12.2 If the ETL and CTL are the same detrain updraft mass @@ -520,7 +543,7 @@ SUBROUTINE CONVECT_UPDRAFT_SHAL( CVP_SHAL, CVPEXT, CST, D, CONVPAR, ! ------------------------------------------------------- ! !IWORK(:) = MIN( KPBL(:), KLCL(:) - 1 ) -IWORK(D%NIB:D%NIE) = KPBL(D%NIB:D%NIE) +IWORK(:) = KPBL(:) DO JI = D%NIB, D%NIE JK = KDPL(JI) JKP = IWORK(JI) diff --git a/conv/ini_convpar.h b/conv/ini_convpar.h index d8766f2c..42ba82a2 100644 --- a/conv/ini_convpar.h +++ b/conv/ini_convpar.h @@ -1,92 +1 @@ -! ######spl - SUBROUTINE INI_CONVPAR(CONVPAR) - USE YOMHOOK , ONLY : LHOOK, DR_HOOK, JPHOOK -! ###################### -! -!!**** *INI_CONVPAR * - routine to initialize the constants modules -!! -!! PURPOSE -!! ------- -! The purpose of this routine is to initialize the constants -! stored in modules MODD_CONVPAR, MODD_CST, MODD_CONVPAREXT. -! -! -!!** METHOD -!! ------ -!! The deep convection constants are set to their numerical values -!! -!! -!! EXTERNAL -!! -------- -!! -!! IMPLICIT ARGUMENTS -!! ------------------ -!! Module MODD_CONVPAR : contains deep convection constants -!! -!! REFERENCE -!! --------- -!! Book2 of the documentation (module MODD_CONVPAR, routine INI_CONVPAR) -!! -!! -!! AUTHOR -!! ------ -!! P. BECHTOLD * Laboratoire d'Aerologie * -!! -!! MODIFICATIONS -!! ------------- -!! Original 26/03/96 -!! Last modified 15/04/98 adapted for ARPEGE -!------------------------------------------------------------------------------- -! -!* 0. DECLARATIONS -! ------------ -! -USE MODD_CONVPAR, ONLY : CONVPAR_T -IMPLICIT NONE -TYPE(CONVPAR_T), INTENT(INOUT) :: CONVPAR - -! -! -!------------------------------------------------------------------------------- -! -!* 1. Set the thermodynamical and numerical constants for -! the deep convection parameterization -! --------------------------------------------------- -! -! -REAL(KIND=JPHOOK) :: ZHOOK_HANDLE -IF (LHOOK) CALL DR_HOOK('INI_CONVPAR',0,ZHOOK_HANDLE) -CONVPAR%XA25 = 625.E6 ! 25 km x 25 km reference grid area -! -CONVPAR%XCRAD = 1500. ! cloud radius -CONVPAR%XCDEPTH = 2.5E3 ! minimum necessary cloud depth -CONVPAR%XENTR = 0.03 ! entrainment constant (m/Pa) = 0.2 (m) -! -CONVPAR%XZLCL = 3.5E3 ! maximum allowed allowed height - ! difference between the surface and the LCL -CONVPAR%XZPBL = 60.E2 ! minimum mixed layer depth to sustain convection -CONVPAR%XWTRIG = 6.00 ! constant in vertical velocity trigger -! -! -CONVPAR%XNHGAM = 1.3333 ! accounts for non-hydrost. pressure - ! in buoyancy term of w equation - ! = 2 / (1+gamma) -CONVPAR%XTFRZ1 = 268.16 ! begin of freezing interval -CONVPAR%XTFRZ2 = 248.16 ! end of freezing interval -! -CONVPAR%XRHDBC = 0.9 ! relative humidity below cloud in downdraft - -CONVPAR%XRCONV = 0.015 ! constant in precipitation conversion -CONVPAR%XSTABT = 0.75 ! factor to assure stability in fractional time - ! integration, routine CONVECT_CLOSURE -CONVPAR%XSTABC = 0.95 ! factor to assure stability in CAPE adjustment, - ! routine CONVECT_CLOSURE -CONVPAR%XUSRDPTH = 165.E2 ! pressure thickness used to compute updraft - ! moisture supply rate for downdraft -CONVPAR%XMELDPTH = 100.E2 ! layer (Pa) through which precipitation melt is - ! allowed below downdraft -CONVPAR%XUVDP = 0.7 ! constant for pressure perturb in momentum transport -! -! -IF (LHOOK) CALL DR_HOOK('INI_CONVPAR',1,ZHOOK_HANDLE) -END SUBROUTINE INI_CONVPAR +! apparently dead code. REK diff --git a/micro/ice_adjust.F90 b/micro/ice_adjust.F90 index b480da71..611c3b17 100644 --- a/micro/ice_adjust.F90 +++ b/micro/ice_adjust.F90 @@ -103,6 +103,7 @@ SUBROUTINE ICE_ADJUST (D, CST, ICEP, NEBN, TURBN, PARAMI, BUCONF, KRR, & ! P. Wautelet 02/2020: use the new data structures and subroutines for budgets !! 2020-12 U. Andrae : Introduce SPP for HARMONIE-AROME !! R. El Khatib 24-Aug-2021 Optimizations +!! R. El Khatib 24-Oct-2023 Re-vectorize ;-) !! !------------------------------------------------------------------------------- ! @@ -216,6 +217,8 @@ SUBROUTINE ICE_ADJUST (D, CST, ICEP, NEBN, TURBN, PARAMI, BUCONF, KRR, & ! REAL, DIMENSION(D%NIJT,D%NKT) :: ZSIGS, ZSRCS REAL, DIMENSION(D%NIJT) :: ZSIGQSAT +LOGICAL :: LLNONE, LLTRIANGLE, LLHLC_H, LLHLI_H + REAL(KIND=JPHOOK) :: ZHOOK_HANDLE ! !------------------------------------------------------------------------------- @@ -319,6 +322,13 @@ SUBROUTINE ICE_ADJUST (D, CST, ICEP, NEBN, TURBN, PARAMI, BUCONF, KRR, & END IF ENDDO ELSE !NEBN%LSUBG_COND case + ! Tests on characters strings can break the vectorization, or at least they would + ! slow down considerably the performance of a vector loop. One should use tests on + ! reals, integers or booleans only. REK. + LLNONE=PARAMI%CSUBG_MF_PDF=='NONE' + LLTRIANGLE=PARAMI%CSUBG_MF_PDF=='TRIANGLE' + LLHLC_H=PRESENT(PHLC_HRC).AND.PRESENT(PHLC_HCF) + LLHLI_H=PRESENT(PHLI_HRI).AND.PRESENT(PHLI_HCF) DO JIJ=IIJB,IIJE !We limit PRC_MF+PRI_MF to PRVS*PTSTEP to avoid negative humidity ZW1=PRC_MF(JIJ,JK)/PTSTEP @@ -334,14 +344,14 @@ SUBROUTINE ICE_ADJUST (D, CST, ICEP, NEBN, TURBN, PARAMI, BUCONF, KRR, & PTHS(JIJ,JK) = PTHS(JIJ,JK) + & (ZW1 * ZLV(JIJ,JK) + ZW2 * ZLS(JIJ,JK)) / ZCPH(JIJ,JK) / PEXNREF(JIJ,JK) ! - IF(PRESENT(PHLC_HRC) .AND. PRESENT(PHLC_HCF)) THEN + IF(LLHLC_H) THEN ZCRIAUT=ICEP%XCRIAUTC/PRHODREF(JIJ,JK) - IF(PARAMI%CSUBG_MF_PDF=='NONE')THEN + IF(LLNONE)THEN IF(ZW1*PTSTEP>PCF_MF(JIJ,JK) * ZCRIAUT) THEN PHLC_HRC(JIJ,JK)=PHLC_HRC(JIJ,JK)+ZW1*PTSTEP PHLC_HCF(JIJ,JK)=MIN(1.,PHLC_HCF(JIJ,JK)+PCF_MF(JIJ,JK)) ENDIF - ELSEIF(PARAMI%CSUBG_MF_PDF=='TRIANGLE')THEN + ELSEIF(LLTRIANGLE)THEN !ZHCF is the precipitating part of the *cloud* and not of the grid cell IF(ZW1*PTSTEP>PCF_MF(JIJ,JK)*ZCRIAUT) THEN ZHCF=1.-.5*(ZCRIAUT*PCF_MF(JIJ,JK) / MAX(1.E-20, ZW1*PTSTEP))**2 @@ -362,14 +372,14 @@ SUBROUTINE ICE_ADJUST (D, CST, ICEP, NEBN, TURBN, PARAMI, BUCONF, KRR, & PHLC_HRC(JIJ,JK)=PHLC_HRC(JIJ,JK)+ZHR ENDIF ENDIF - IF(PRESENT(PHLI_HRI) .AND. PRESENT(PHLI_HCF)) THEN + IF(LLHLI_H) THEN ZCRIAUT=MIN(ICEP%XCRIAUTI,10**(ICEP%XACRIAUTI*(ZT(JIJ,JK)-CST%XTT)+ICEP%XBCRIAUTI)) - IF(PARAMI%CSUBG_MF_PDF=='NONE')THEN + IF(LLNONE)THEN IF(ZW2*PTSTEP>PCF_MF(JIJ,JK) * ZCRIAUT) THEN PHLI_HRI(JIJ,JK)=PHLI_HRI(JIJ,JK)+ZW2*PTSTEP PHLI_HCF(JIJ,JK)=MIN(1.,PHLI_HCF(JIJ,JK)+PCF_MF(JIJ,JK)) ENDIF - ELSEIF(PARAMI%CSUBG_MF_PDF=='TRIANGLE')THEN + ELSEIF(LLTRIANGLE)THEN !ZHCF is the precipitating part of the *cloud* and not of the grid cell IF(ZW2*PTSTEP>PCF_MF(JIJ,JK)*ZCRIAUT) THEN ZHCF=1.-.5*(ZCRIAUT*PCF_MF(JIJ,JK) / (ZW2*PTSTEP))**2 diff --git a/micro/mode_ini_snow.F90 b/micro/mode_ini_snow.F90 index 313659cb..a6d6ebd0 100644 --- a/micro/mode_ini_snow.F90 +++ b/micro/mode_ini_snow.F90 @@ -79,7 +79,7 @@ SUBROUTINE INI_SNOW ( KLUOUT ) REAL(KIND=JPHOOK) :: ZHOOK_HANDLE -IF (LHOOK) CALL DR_HOOK('INI_RAIN_ICE',0,ZHOOK_HANDLE) +IF (LHOOK) CALL DR_HOOK('INI_SNOW',0,ZHOOK_HANDLE) XCCS = XFRMIN(16) diff --git a/micro/mode_rain_ice_old_nucleation.F90 b/micro/mode_rain_ice_old_nucleation.F90 index 563a3c15..7c240de3 100644 --- a/micro/mode_rain_ice_old_nucleation.F90 +++ b/micro/mode_rain_ice_old_nucleation.F90 @@ -12,6 +12,7 @@ MODULE MODE_RAIN_ICE_OLD_NUCLEATION SUBROUTINE RAIN_ICE_OLD_NUCLEATION(D, CST, ICEP, KSIZE, OCND2, LMODICEDEP, KRR, PTSTEP, & PTHT, PPABST, PEXNREF, PICLDFR, PRHODJ, PRHODREF, & PRVT, PRCT, PRRT, PRIT, PRST, PRGT, & + OAERONRT, OAEIFN, PIFNNC, & PTHS, PRVS, PRIS, PCIT, & PICENU, PT, PZZZ, & PRHT) @@ -52,6 +53,9 @@ SUBROUTINE RAIN_ICE_OLD_NUCLEATION(D, CST, ICEP, KSIZE, OCND2, LMODICEDEP, KRR, REAL, DIMENSION(D%NIT, D%NKT), INTENT(IN) :: PRIT ! Pristine ice m.r. at t REAL, DIMENSION(D%NIT, D%NKT), INTENT(IN) :: PRST ! Snow/aggregate m.r. at t REAL, DIMENSION(D%NIT, D%NKT), INTENT(IN) :: PRGT ! Graupel/hail m.r. at t + LOGICAL, INTENT(IN) :: OAERONRT ! Switch for nrt aerosols + LOGICAL, INTENT(IN) :: OAEIFN ! Switch to activate ice nuclei + REAL, DIMENSION(D%NIT, D%NKT), INTENT(IN) :: PIFNNC ! Ice freezing nuclei concentration ! REAL, DIMENSION(D%NIT, D%NKT), INTENT(INOUT) :: PTHS ! Theta source REAL, DIMENSION(D%NIT, D%NKT), INTENT(INOUT) :: PRVS ! Water vapor m.r. source @@ -128,6 +132,9 @@ SUBROUTINE RAIN_ICE_OLD_NUCLEATION(D, CST, ICEP, KSIZE, OCND2, LMODICEDEP, KRR, ZREDIN(JL) = REDIN(ZZT(JL)) ZSIFRC(JL) = PICLDFR(I1(JL),I2(JL)) ENDIF + IF (OAERONRT.AND.OAEIFN) THEN + ZZICENU(JL)=MAX(0.01,PIFNNC(I1(JL),I2(JL))*1.0E-6) ! convert m-3 to cm-3 + ENDIF ENDDO diff --git a/micro/modi_rain_ice_old.F90 b/micro/modi_rain_ice_old.F90 index ba188918..ad6e69e0 100644 --- a/micro/modi_rain_ice_old.F90 +++ b/micro/modi_rain_ice_old.F90 @@ -15,6 +15,7 @@ SUBROUTINE RAIN_ICE_OLD (D, CST, PARAMI, ICEP, ICED, BUCONF, PRGT, PTHS, PRVS, PRCS, PRRS, PRIS, PRSS, PRGS, & PINPRC, PINPRR, PEVAP3D, & PINPRS, PINPRG, PSIGS, PSEA, PTOWN, & + OAERONRT, OAEIFN, PCLDROP, PIFNNC, & TBUDGETS, KBUDGETS, & PICENU, PKGN_ACON, PKGN_SBGR, & PRHT, PRHS, PINPRH, PFPR) @@ -98,6 +99,12 @@ SUBROUTINE RAIN_ICE_OLD (D, CST, PARAMI, ICEP, ICED, BUCONF, REAL, DIMENSION(D%NIT), INTENT(OUT) :: PINPRG! Graupel instant precip REAL, DIMENSION(D%NIT), INTENT(IN) :: PSEA ! Sea Mask REAL, DIMENSION(D%NIT), INTENT(IN) :: PTOWN! Fraction that is town +! nrt aerosol +LOGICAL, INTENT(IN) :: OAERONRT ! Switch for nrt aerosols +LOGICAL, INTENT(IN) :: OAEIFN ! Switch to activate ice nuclei +REAL, DIMENSION(D%NIT,D%NKT), INTENT(IN) :: PCLDROP ! Activated Condensation nuclei (CCN) +REAL, DIMENSION(D%NIT,D%NKT), INTENT(IN) :: PIFNNC ! Ice freezing nuclei concentration +! TYPE(TBUDGETDATA), DIMENSION(KBUDGETS), INTENT(INOUT) :: TBUDGETS INTEGER, INTENT(IN) :: KBUDGETS REAL, DIMENSION(D%NIT), INTENT(IN) :: PICENU, PKGN_ACON, PKGN_SBGR diff --git a/micro/rain_ice_old.F90 b/micro/rain_ice_old.F90 index 71174a7f..68a8a2ab 100644 --- a/micro/rain_ice_old.F90 +++ b/micro/rain_ice_old.F90 @@ -10,6 +10,7 @@ SUBROUTINE RAIN_ICE_OLD (D, CST, PARAMI, ICEP, ICED, BUCONF, PRGT, PTHS, PRVS, PRCS, PRRS, PRIS, PRSS, PRGS, & PINPRC, PINPRR, PEVAP3D, & PINPRS, PINPRG, PSIGS, PSEA, PTOWN, & + OAERONRT, OAEIFN, PCLDROP, PIFNNC, & TBUDGETS, KBUDGETS, & PICENU, PKGN_ACON, PKGN_SBGR, & PRHT, PRHS, PINPRH, PFPR) @@ -167,6 +168,7 @@ SUBROUTINE RAIN_ICE_OLD (D, CST, PARAMI, ICEP, ICED, BUCONF, !! Sedimented ice should be preciptation !! (U. Andrae Dec 2020) Introduce SPP for HARMONIE-AROME !! (C. Wittmann Jan 2021) Introduce sublimation factor tuning +!! (D. Martin-Perez, 2021) nrt Aerosol ! ! !* 0. DECLARATIONS @@ -268,6 +270,12 @@ SUBROUTINE RAIN_ICE_OLD (D, CST, PARAMI, ICEP, ICED, BUCONF, REAL, DIMENSION(D%NIT), INTENT(OUT) :: PINPRG! Graupel instant precip REAL, DIMENSION(D%NIT), INTENT(IN) :: PSEA ! Sea Mask REAL, DIMENSION(D%NIT), INTENT(IN) :: PTOWN! Fraction that is town +! nrt aerosol +LOGICAL, INTENT(IN) :: OAERONRT ! Switch for nrt aerosols +LOGICAL, INTENT(IN) :: OAEIFN ! Switch to activate ice nuclei +REAL, DIMENSION(D%NIT,D%NKT), INTENT(IN) :: PCLDROP ! Activated Condensation nuclei (CCN) +REAL, DIMENSION(D%NIT,D%NKT), INTENT(IN) :: PIFNNC ! Ice freezing nuclei concentration +! TYPE(TBUDGETDATA), DIMENSION(KBUDGETS), INTENT(INOUT) :: TBUDGETS INTEGER, INTENT(IN) :: KBUDGETS REAL, DIMENSION(D%NIT), INTENT(IN) :: PICENU, PKGN_ACON, PKGN_SBGR @@ -481,7 +489,19 @@ SUBROUTINE RAIN_ICE_OLD (D, CST, PARAMI, ICEP, ICED, BUCONF, ENDDO ENDIF - ZCONC3D(:,D%NKTE)= ZCONC3D(:,D%NKTE)*MAX(0.001,ICEP%XFRMIN(22)) + ZCONC3D(:,D%NKTE) = ZCONC3D(:,D%NKTE)*MAX(0.001,ICEP%XFRMIN(22)) + + !Consider CCN obtained from near real time aerosol mixing ratio fields + IF (OAERONRT) THEN + DO JK=D%NKTB,D%NKTE + ZCONC3D(:,JK) = PCLDROP(:,JK) + ZLBC(:,JK) = 0.5* (ICED%XLBC(2)+ICED%XLBC(1)) ! Assume "average" distr. func + ZFSEDC(:,JK) = 0.5* (ICEP%XFSEDC(2)+ICEP%XFSEDC(1)) + ZFSEDC(:,JK) = MAX(MIN(ICEP%XFSEDC(1),ICEP%XFSEDC(2)),ZFSEDC(:,JK)) + ZRAY(:,JK) = 0.5*(0.5*GAMMA(ICED%XNUC+1.0/ICED%XALPHAC)/(GAMMA(ICED%XNUC)) + & + 0.5*GAMMA(ICED%XNUC2+1.0/ICED%XALPHAC2)/(GAMMA(ICED%XNUC2))) + ENDDO + ENDIF ZRAY(:,:) = MAX(1.,ZRAY(:,:)) ZLBC(:,:) = MAX(MIN(ICED%XLBC(1),ICED%XLBC(2)),ZLBC(:,:)) @@ -502,6 +522,7 @@ SUBROUTINE RAIN_ICE_OLD (D, CST, PARAMI, ICEP, ICED, BUCONF, OCND2, LMODICEDEP, KRR, PTSTEP, & PTHT, PPABST, PEXNREF, PICLDFR, PRHODJ, PRHODREF, & PRVT, PRCT, PRRT, PRIT, PRST, PRGT, & + OAERONRT, OAEIFN, PIFNNC, & PTHS, PRVS, PRIS, PCIT, & PICENU, ZT, ZZZZ, & PRHT) diff --git a/turb/modd_turbn.F90 b/turb/modd_turbn.F90 index 74a4d559..c61f67fd 100644 --- a/turb/modd_turbn.F90 +++ b/turb/modd_turbn.F90 @@ -77,6 +77,7 @@ MODULE MODD_TURB_n LOGICAL :: LSIG_CONV !< Switch for computing Sigma_s due to convection ! LOGICAL :: LHARAT !< if true RACMO turbulence is used + LOGICAL :: LBL89TOP !< if true modification in BL89 at PBL top LOGICAL :: LRMC01 !< Switch for computing separate mixing and dissipative length in the SBL !! according to Redelsperger, Mahe & Carlotti 2001 CHARACTER(LEN=4) :: CTOM !< type of Third Order Moments: @@ -99,6 +100,8 @@ MODULE MODD_TURB_n REAL :: XCOEFHGRADTHL !< coeff applied to thl contribution REAL :: XCOEFHGRADRM !< coeff applied to mixing ratio contribution REAL :: XALTHGRAD !< altitude from which to apply the Leonard terms + LOGICAL :: LGOGER ! < logical switch for the computation of the Goger Terms + REAL :: XSMAG ! < dimensionless Smagorinsky constant REAL :: XCLDTHOLD !< cloud threshold to apply the Leonard terms: !! negative value to apply everywhere; !! 0.000001 applied only inside the clouds ri+rc > 10**-6 kg/kg @@ -134,6 +137,7 @@ MODULE MODD_TURB_n LOGICAL, POINTER :: LSIG_CONV=>NULL() LOGICAL, POINTER :: LRMC01=>NULL() LOGICAL, POINTER :: LHARAT=>NULL() +LOGICAL, POINTER :: LBL89TOP=>NULL() CHARACTER(LEN=4),POINTER :: CTOM=>NULL() REAL, DIMENSION(:,:), POINTER :: XBL_DEPTH=>NULL() REAL, DIMENSION(:,:), POINTER :: XSBL_DEPTH=>NULL() @@ -151,6 +155,8 @@ MODULE MODD_TURB_n REAL, POINTER :: XCOEFHGRADTHL=>NULL() REAL, POINTER :: XCOEFHGRADRM=>NULL() REAL, POINTER :: XALTHGRAD=>NULL() +LOGICAL, POINTER :: LGOGER=>NULL() +REAL, POINTER :: XSMAG=>NULL() REAL, POINTER :: XCLDTHOLD=>NULL() REAL, POINTER :: XLINI=>NULL() LOGICAL, POINTER :: LROTATE_WIND=>NULL() @@ -164,8 +170,9 @@ MODULE MODD_TURB_n LSIG_CONV,LRMC01,CTOM,& XTKEMIN,XCED,XCTP,XCADAP,& LLEONARD,XCOEFHGRADTHL, XCOEFHGRADRM, & - XALTHGRAD, XCLDTHOLD, XLINI, LHARAT, & - LPROJQITURB, LSMOOTH_PRANDTL, XMINSIGS + XALTHGRAD, LGOGER, XSMAG, XCLDTHOLD, XLINI, LHARAT, & + LPROJQITURB, LSMOOTH_PRANDTL, XMINSIGS, & + LBL89TOP ! !------------------------------------------------------------------------------- ! @@ -216,6 +223,7 @@ SUBROUTINE TURB_GOTO_MODEL(KFROM, KTO) CTURBDIM=>TURB_MODEL(KTO)%CTURBDIM LTURB_FLX=>TURB_MODEL(KTO)%LTURB_FLX LHARAT=>TURB_MODEL(KTO)%LHARAT +LBL89TOP=>TURB_MODEL(KTO)%LBL89TOP LTURB_DIAG=>TURB_MODEL(KTO)%LTURB_DIAG LSIG_CONV=>TURB_MODEL(KTO)%LSIG_CONV LRMC01=>TURB_MODEL(KTO)%LRMC01 @@ -236,6 +244,8 @@ SUBROUTINE TURB_GOTO_MODEL(KFROM, KTO) XCOEFHGRADTHL=>TURB_MODEL(KTO)%XCOEFHGRADTHL XCOEFHGRADRM=>TURB_MODEL(KTO)%XCOEFHGRADRM XALTHGRAD=>TURB_MODEL(KTO)%XALTHGRAD +LGOGER=>TURB_MODEL(KTO)%LGOGER +XSMAG=>TURB_MODEL(KTO)%XSMAG XCLDTHOLD=>TURB_MODEL(KTO)%XCLDTHOLD XLINI=>TURB_MODEL(KTO)%XLINI LROTATE_WIND=>TURB_MODEL(KTO)%LROTATE_WIND @@ -339,9 +349,12 @@ SUBROUTINE TURBN_INIT(HPROGRAM, KUNITNML, LDNEEDNAM, KLUOUT, & XCOEFHGRADTHL = 1.0 XCOEFHGRADRM = 1.0 XALTHGRAD = 2000.0 + LGOGER=.FALSE. + XSMAG=0.20 XCLDTHOLD = -1.0 XLINI=0.1 !old value: 10. LHARAT=.FALSE. + LBL89TOP=.FALSE. LROTATE_WIND=.FALSE. LTKEMINTURB=.TRUE. LPROJQITURB=.TRUE. diff --git a/turb/mode_bl89.F90 b/turb/mode_bl89.F90 index 512ae2d5..a2dbe04a 100644 --- a/turb/mode_bl89.F90 +++ b/turb/mode_bl89.F90 @@ -93,7 +93,8 @@ SUBROUTINE BL89(D,CST,CSTURB,TURBN,PZZ,PDZZ,PTHVREF,PTHLM,KRR,PRM,PTKEM,PSHEAR,P ! Increment of Virtual Potential Temp between two following levels REAL, DIMENSION(D%NIJT,D%NKT) :: ZHLVPT ! Virtual Potential Temp at half levels -REAL, DIMENSION(D%NIJT) :: ZLWORK,ZINTE +REAL, DIMENSION(D%NIJT,D%NKT) :: ZLWORK +REAL, DIMENSION(D%NIJT) :: ZINTE ! ! downwards then upwards vertical displacement, ! ! residual internal energy, ! ! residual potential energy @@ -253,7 +254,7 @@ SUBROUTINE BL89(D,CST,CSTURB,TURBN,PZZ,PDZZ,PTHVREF,PTHLM,KRR,PRM,PTKEM,PSHEAR,P 2. * ZINTE(JIJ) * & (ZG_O_THVREF(JIJ,JK) * ZDELTVPT(JIJ,JKK)/ PDZZ(JIJ,JKK))))) / & (ZG_O_THVREF(JIJ,JK) * ZDELTVPT(JIJ,JKK) / PDZZ(JIJ,JKK)) - ZLWORK(JIJ)=ZLWORK(JIJ)+ZTEST0*(ZTEST*ZLWORK1+(1-ZTEST)*ZLWORK2) + ZLWORK(JIJ,JK)=ZLWORK(JIJ,JK)+ZTEST0*(ZTEST*ZLWORK1+(1-ZTEST)*ZLWORK2) ZINTE(JIJ) = ZINTE(JIJ) - ZPOTE END DO ENDIF @@ -264,7 +265,7 @@ SUBROUTINE BL89(D,CST,CSTURB,TURBN,PZZ,PDZZ,PTHVREF,PTHLM,KRR,PRM,PTKEM,PSHEAR,P ! ----------------------------------------------- ! DO JIJ=IIJB,IIJE - PLMDN(JIJ,JK)=MIN(ZLWORK(JIJ),0.5*(PZZ(JIJ,JK)+PZZ(JIJ,JK+IKL))-PZZ(JIJ,IKB)) + PLMDN(JIJ,JK)=MIN(ZLWORK(JIJ,JK),0.5*(PZZ(JIJ,JK)+PZZ(JIJ,JK+IKL))-PZZ(JIJ,IKB)) END DO ! !------------------------------------------------------------------------------- @@ -273,7 +274,7 @@ SUBROUTINE BL89(D,CST,CSTURB,TURBN,PZZ,PDZZ,PTHVREF,PTHLM,KRR,PRM,PTKEM,PSHEAR,P ! ----------------------------------------- ! ZINTE(IIJB:IIJE)=PTKEM(IIJB:IIJE,JK) - ZLWORK(IIJB:IIJE)=0. + ZLWORK(IIJB:IIJE,JK)=0. ZTESTM=1. ! DO JKK=JK+IKL,IKE,IKL @@ -298,19 +299,28 @@ SUBROUTINE BL89(D,CST,CSTURB,TURBN,PZZ,PDZZ,PTHVREF,PTHLM,KRR,PRM,PTKEM,PSHEAR,P + 2. * ZINTE(JIJ) * & (ZG_O_THVREF(JIJ,JK)* ZDELTVPT(JIJ,JKK)/PDZZ(JIJ,JKK))))) / & (ZG_O_THVREF(JIJ,JK) * ZDELTVPT(JIJ,JKK) / PDZZ(JIJ,JKK)) - ZLWORK(JIJ)=ZLWORK(JIJ)+ZTEST0*(ZTEST*ZLWORK1+(1-ZTEST)*ZLWORK2) + ZLWORK(JIJ,JK)=ZLWORK(JIJ,JK)+ZTEST0*(ZTEST*ZLWORK1+(1-ZTEST)*ZLWORK2) ZINTE(JIJ) = ZINTE(JIJ) - ZPOTE END DO ENDIF END DO ! +!* Maximal length between JK and the levels below (as in ARPEGE) +! + IF (TURBN%LBL89TOP) THEN + DO JIJ=IIJB,IIJE + ZLWORK(JIJ,JK) = MAX(ZLWORK(JIJ,JK),ZLWORK(JIJ,JK-IKL)-PDZZ(JIJ,JK)) + ENDDO + END IF + +! !------------------------------------------------------------------------------- ! !* 7. final mixing length ! DO JIJ=IIJB,IIJE ZLWORK1=MAX(PLMDN(JIJ,JK),1.E-10_MNHREAL) - ZLWORK2=MAX(ZLWORK(JIJ),1.E-10_MNHREAL) + ZLWORK2=MAX(ZLWORK(JIJ,JK),1.E-10_MNHREAL) ZPOTE = ZLWORK1 / ZLWORK2 ZLWORK2=1.d0 + ZPOTE**TURBN%XBL89EXP PLM(JIJ,JK) = ZLWORK1*(2./ZLWORK2)**TURBN%XUSRBL89 diff --git a/turb/mode_compute_updraft.F90 b/turb/mode_compute_updraft.F90 index f074794e..c067c702 100644 --- a/turb/mode_compute_updraft.F90 +++ b/turb/mode_compute_updraft.F90 @@ -60,7 +60,7 @@ SUBROUTINE COMPUTE_UPDRAFT(D,CST,NEBN,PARAMMF,TURBN,CSTURB, & !! S. Riette Apr 2013: improvement of continuity at the condensation level !! R.Honnert Oct 2016 : Add ZSURF and Update with AROME !! Q.Rodier 01/2019 : support RM17 mixing length -!! R.Honnert 01/2019 : add LGZ (reduction of the mass-flux surface closure with the resolution) +!! R.Honnert 01/2019 : remove ZSURF add LGZ (reduction of the MF surf closure with resolution) !! S. Riette 06/2022: compute_entr_detr is inlined !! -------------------------------------------------------------------------- ! diff --git a/turb/modi_turb.F90 b/turb/modi_turb.F90 index da938a34..e1ed0675 100644 --- a/turb/modi_turb.F90 +++ b/turb/modi_turb.F90 @@ -6,7 +6,8 @@ MODULE MODI_TURB INTERFACE ! SUBROUTINE TURB(CST,CSTURB,BUCONF,TURBN,NEBN,D,TLES, & - & KMI,KRR,KRRL,KRRI,HLBCX,HLBCY,KGRADIENTS,KHALO, & + & KMI,KRR,KRRL,KRRI,HLBCX,HLBCY,KGRADIENTSLEO, & + & KGRADIENTSGOG,KHALO, & & KSPLIT,KMODEL_CL,KSV,KSV_LGBEG,KSV_LGEND, & & KSV_LIMA_NR, KSV_LIMA_NS, KSV_LIMA_NG, KSV_LIMA_NH, & & O2D,ONOMIXLG,OFLAT,OCOUPLES,OBLOWSNOW,OIBM,OFLYER, & @@ -16,7 +17,7 @@ SUBROUTINE TURB(CST,CSTURB,BUCONF,TURBN,NEBN,D,TLES, & & PTSTEP,TPFILE, & & PDXX,PDYY,PDZZ,PDZX,PDZY,PZZ, & & PDIRCOSXW,PDIRCOSYW,PDIRCOSZW,PCOSSLOPE,PSINSLOPE, & - & PRHODJ,PTHVREF,PHGRAD,PZS, & + & PRHODJ,PTHVREF,PHGRADLEO,PHGRADGOG,PZS, & & PSFTH,PSFRV,PSFSV,PSFU,PSFV, & & PPABST,PUT,PVT,PWT,PTKET,PSVT,PSRCT, & & PLENGTHM,PLENGTHH,MFMOIST, & @@ -51,7 +52,8 @@ SUBROUTINE TURB(CST,CSTURB,BUCONF,TURBN,NEBN,D,TLES, & TYPE(TURB_t), INTENT(IN) :: TURBN ! modn_turbn (turb namelist) structure TYPE(NEB_t), INTENT(IN) :: NEBN ! modd_nebn structure TYPE(TLES_t), INTENT(INOUT) :: TLES ! modd_les structure -INTEGER, INTENT(IN) :: KGRADIENTS ! Number of stored horizontal gradients +INTEGER, INTENT(IN) :: KGRADIENTSLEO ! Number of stored horizontal gradients +INTEGER, INTENT(IN) :: KGRADIENTSGOG ! Number of stored horizontal gradients INTEGER, INTENT(IN) :: KMI ! model index number INTEGER, INTENT(IN) :: KRR ! number of moist var. INTEGER, INTENT(IN) :: KRRL ! number of liquid water var. @@ -92,7 +94,8 @@ SUBROUTINE TURB(CST,CSTURB,BUCONF,TURBN,NEBN,D,TLES, & REAL, DIMENSION(D%NIJT,D%NKT), INTENT(IN) :: MFMOIST ! moist mass flux dual scheme REAL, DIMENSION(D%NIJT,D%NKT), INTENT(IN) :: PTHVREF ! Virtual Potential ! Temperature of the reference state -REAL, DIMENSION(D%NIJT,D%NKT,KGRADIENTS), INTENT(IN) :: PHGRAD ! horizontal gradients +REAL, DIMENSION(D%NIJT,D%NKT,KGRADIENTSLEO), INTENT(IN) :: PHGRADLEO ! horizontal gradients +REAL, DIMENSION(D%NIJT,D%NKT,KGRADIENTSGOG), INTENT(IN) :: PHGRADGOG ! horizontal gradients ! REAL, DIMENSION(D%NIJT), INTENT(IN) :: PSFTH,PSFRV, & ! normal surface fluxes of theta and Rv diff --git a/turb/turb.F90 b/turb/turb.F90 index e738ffdb..5e6257ea 100644 --- a/turb/turb.F90 +++ b/turb/turb.F90 @@ -4,7 +4,8 @@ !MNH_LIC for details. version 1. !----------------------------------------------------------------- SUBROUTINE TURB(CST,CSTURB,BUCONF,TURBN,NEBN,D,TLES, & - & KMI,KRR,KRRL,KRRI,HLBCX,HLBCY,KGRADIENTS,KHALO, & + & KMI,KRR,KRRL,KRRI,HLBCX,HLBCY,KGRADIENTSLEO, & + & KGRADIENTSGOG,KHALO, & & KSPLIT,KMODEL_CL,KSV,KSV_LGBEG,KSV_LGEND, & & KSV_LIMA_NR, KSV_LIMA_NS, KSV_LIMA_NG, KSV_LIMA_NH, & & O2D,ONOMIXLG,OFLAT,OCOUPLES,OBLOWSNOW,OIBM,OFLYER, & @@ -14,7 +15,7 @@ SUBROUTINE TURB(CST,CSTURB,BUCONF,TURBN,NEBN,D,TLES, & & PTSTEP,TPFILE, & & PDXX,PDYY,PDZZ,PDZX,PDZY,PZZ, & & PDIRCOSXW,PDIRCOSYW,PDIRCOSZW,PCOSSLOPE,PSINSLOPE, & - & PRHODJ,PTHVREF,PHGRAD,PZS, & + & PRHODJ,PTHVREF,PHGRADLEO,PHGRADGOG,PZS, & & PSFTH,PSFRV,PSFSV,PSFU,PSFV, & & PPABST,PUT,PVT,PWT,PTKET,PSVT,PSRCT, & & PLENGTHM,PLENGTHH,MFMOIST, & @@ -248,7 +249,7 @@ SUBROUTINE TURB(CST,CSTURB,BUCONF,TURBN,NEBN,D,TLES, & USE MODD_CST, ONLY: CST_t USE MODD_CTURB, ONLY: CSTURB_t USE MODD_DIMPHYEX, ONLY: DIMPHYEX_t -USE MODD_FIELD, ONLY: TFIELDMETADATA, TYPEREAL +USE MODD_FIELD, ONLY: TFIELDMETADATA, TYPEREAL, NMNHDIM_UNKNOWN USE MODD_IO, ONLY: TFILEDATA USE MODD_LES, ONLY: TLES_t USE MODD_PARAMETERS, ONLY: JPVEXT_TURB, XUNDEF @@ -292,7 +293,8 @@ SUBROUTINE TURB(CST,CSTURB,BUCONF,TURBN,NEBN,D,TLES, & TYPE(TURB_t), INTENT(IN) :: TURBN ! modn_turbn (turb namelist) structure TYPE(NEB_t), INTENT(IN) :: NEBN ! modd_nebn structure TYPE(TLES_t), INTENT(INOUT) :: TLES ! modd_les structure -INTEGER, INTENT(IN) :: KGRADIENTS ! Number of stored horizontal gradients +INTEGER, INTENT(IN) :: KGRADIENTSLEO ! Number of stored horizontal gradients for Moeng scheme +INTEGER, INTENT(IN) :: KGRADIENTSGOG ! Number of stored horizontal gradients for Goger scheme INTEGER, INTENT(IN) :: KMI ! model index number INTEGER, INTENT(IN) :: KRR ! number of moist var. INTEGER, INTENT(IN) :: KRRL ! number of liquid water var. @@ -333,7 +335,8 @@ SUBROUTINE TURB(CST,CSTURB,BUCONF,TURBN,NEBN,D,TLES, & REAL, DIMENSION(D%NIJT,D%NKT), INTENT(IN) :: MFMOIST ! moist mass flux dual scheme REAL, DIMENSION(D%NIJT,D%NKT), INTENT(IN) :: PTHVREF ! Virtual Potential ! Temperature of the reference state -REAL, DIMENSION(D%NIJT,D%NKT,KGRADIENTS), INTENT(IN) :: PHGRAD ! horizontal gradients +REAL, DIMENSION(D%NIJT,D%NKT,KGRADIENTSLEO), INTENT(IN) :: PHGRADLEO ! horizontal gradients in Moeng +REAL, DIMENSION(D%NIJT,D%NKT,KGRADIENTSGOG), INTENT(IN) :: PHGRADGOG ! horizontal gradients in Goger ! REAL, DIMENSION(D%NIJT), INTENT(IN) :: PSFTH,PSFRV, & ! normal surface fluxes of theta and Rv @@ -498,6 +501,7 @@ SUBROUTINE TURB(CST,CSTURB,BUCONF,TURBN,NEBN,D,TLES, & ! ! BL89 mixing length near the surface ! REAL :: ZTIME1, ZTIME2 +REAL :: ZDUDX, ZDVDY, ZDUDY, ZDVDX, ZK ! work values for Göger TYPE(TFIELDMETADATA) :: TZFIELD ! !* 1.PRELIMINARIES @@ -1034,7 +1038,7 @@ SUBROUTINE TURB(CST,CSTURB,BUCONF,TURBN,NEBN,D,TLES, & END IF CALL TURB_VER(D,CST,CSTURB,TURBN,NEBN,TLES, & - KRR,KRRL,KRRI,KGRADIENTS, & + KRR,KRRL,KRRI,KGRADIENTSLEO, & OOCEAN, ODEEPOC, OCOMPUTE_SRC, & KSV,KSV_LGBEG,KSV_LGEND, & ZEXPL, O2D, ONOMIXLG, OFLAT, & @@ -1049,7 +1053,7 @@ SUBROUTINE TURB(CST,CSTURB,BUCONF,TURBN,NEBN,D,TLES, & PTKET,ZLM,PLENGTHM,PLENGTHH,ZLEPS,MFMOIST, & ZLOCPEXNM,ZATHETA,ZAMOIST,PSRCT,ZFRAC_ICE, & ZFWTH,ZFWR,ZFTH2,ZFR2,ZFTHR,PBL_DEPTH, & - PSBL_DEPTH,ZLMO,PHGRAD,PZS, & + PSBL_DEPTH,ZLMO,PHGRADLEO,PZS, & PRUS,PRVS,PRWS,PRTHLS,PRRS,PRSVS, & PDP,PTP,PSIGS,PWTH,PWRC,PWSV, & PSSTFL, PSSTFL_C, PSSRFL_C,PSSUFL_C,PSSVFL_C, & @@ -1211,7 +1215,31 @@ SUBROUTINE TURB(CST,CSTURB,BUCONF,TURBN,NEBN,D,TLES, & ENDDO ENDDO END IF -! 6.2 TKE evolution equation + +! 6.2 Horizontal gradients as in Göger et al. (2016) + +IF (TURBN%LGOGER) THEN + ! Add horizontal terms from Göger et al. (2018) + ! Increase the Dyn. Prod. + DO JK=1,IKT + DO JIJ=IIJB,IIJE + !* Computation of the horizontal mixing length + ZK=TURBN%XSMAG**2*PDXX(JIJ,JK)*PDYY(JIJ,JK) + !* Add horizontal terms + ! DUDX=PHGRADGOG(JIJ,JK,1) + ! DUDY=PHGRADGOG(JIJ,JK,2) + ! DVDX=PHGRADGOG(JIJ,JK,3) + ! DVDY=PHGRADGOG(JIJ,JK,4) + PDP(JIJ,JK)=PDP(JIJ,JK)+ZK*& + &(PHGRADGOG(JIJ,JK,1)*PHGRADGOG(JIJ,JK,1) & + &+PHGRADGOG(JIJ,JK,4)*PHGRADGOG(JIJ,JK,4) & + &+0.5*(PHGRADGOG(JIJ,JK,2)+PHGRADGOG(JIJ,JK,3)) & + &*(PHGRADGOG(JIJ,JK,2)+PHGRADGOG(JIJ,JK,3)))**(3./2.) + ENDDO + ENDDO +ENDIF + +! 6.3 TKE evolution equation IF (.NOT. TURBN%LHARAT) THEN ! @@ -1517,8 +1545,9 @@ SUBROUTINE COMPUTE_FUNCTION_THERMO(PALP,PBETA,PGAM,PLTT,PC,PT,PEXN,PCP,& ! !* 1.3 saturation mixing ratio at t ! +!YS Added protection (AROME 2024-03-12 crashs) ZRVSAT(JIJ,JK) = ZRVSAT(JIJ,JK) & - * ZEPS / ( PPABST(JIJ,JK) - ZRVSAT(JIJ,JK) ) + * ZEPS / MAX(1.E-3, PPABST(JIJ,JK) - ZRVSAT(JIJ,JK) ) ! !* 1.4 compute the saturation mixing ratio derivative (rvs') ! @@ -2052,7 +2081,21 @@ SUBROUTINE CLOUD_MODIF_LM NGRID = 1, & NTYPE = TYPEREAL, & NDIMS = 3, & - LTIMEDEP = .TRUE. ) + LTIMEDEP = .TRUE., & + !Filling the other fields with the default values, + !because of a bug in nvidia's compiler (NVHPC 23.11) + !That should be removed when the bug is fixed + NDIMLIST = NMNHDIM_UNKNOWN, & + NFILLVALUE = -2147483647, & + XFILLVALUE = 9.9692099683868690e+36,& + NVALIDMIN = -2147483646, & + NVALIDMAX = 2147483647, & + XVALIDMIN = -1.E36, & + XVALIDMAX = 1.E36, & + CLBTYPE = 'NONE' & + !End of default arguments + ) + CALL IO_FIELD_WRITE_PHY(D,TPFILE,TZFIELD,ZLM) ENDIF ! @@ -2091,7 +2134,20 @@ SUBROUTINE CLOUD_MODIF_LM NGRID = 1, & NTYPE = TYPEREAL, & NDIMS = 3, & - LTIMEDEP = .TRUE. ) + LTIMEDEP = .TRUE., & + !Filling the other fields with the default values, + !because of a bug in nvidia's compiler (NVHPC 23.11) + !That should be removed when the bug is fixed + NDIMLIST = NMNHDIM_UNKNOWN, & + NFILLVALUE = -2147483647, & + XFILLVALUE = 9.9692099683868690e+36,& + NVALIDMIN = -2147483646, & + NVALIDMAX = 2147483647, & + XVALIDMIN = -1.E36, & + XVALIDMAX = 1.E36, & + CLBTYPE = 'NONE' & + !End of default arguments + ) CALL IO_FIELD_WRITE_PHY(D,TPFILE,TZFIELD,ZCOEF_AMPL) ! TZFIELD = TFIELDMETADATA( & @@ -2104,7 +2160,20 @@ SUBROUTINE CLOUD_MODIF_LM NGRID = 1, & NTYPE = TYPEREAL, & NDIMS = 3, & - LTIMEDEP = .TRUE. ) + LTIMEDEP = .TRUE., & + !Filling the other fields with the default values, + !because of a bug in nvidia's compiler (NVHPC 23.11) + !That should be removed when the bug is fixed + NDIMLIST = NMNHDIM_UNKNOWN, & + NFILLVALUE = -2147483647, & + XFILLVALUE = 9.9692099683868690e+36,& + NVALIDMIN = -2147483646, & + NVALIDMAX = 2147483647, & + XVALIDMIN = -1.E36, & + XVALIDMAX = 1.E36, & + CLBTYPE = 'NONE' & + !End of default arguments + ) CALL IO_FIELD_WRITE_PHY(D,TPFILE,TZFIELD,ZLM_CLOUD) ! ENDIF