From 79d81425112e7287f522e9efcdbdb883bf87b528 Mon Sep 17 00:00:00 2001 From: Alexandre MARY Date: Thu, 12 Sep 2024 16:10:59 +0000 Subject: [PATCH 1/3] Cf. following IAL commits to 49T1 21d95b428e Merge remote-tracking branch 'antoinettealiasMF/alias_CY49T0_gmgec' into accord_CY49T0_to_T1 e1f45ecf4c Workaround for NVHPC 23.11 4ef4a377db Merge remote-tracking branch 'BertvanUlft/ulft_CY49T0_HIRLAM_physics_aerosol' into accord_CY49T0_to_T1 302034816c Merge remote-tracking branch 'pmarguinaud/rapsgpu' into accord_CY49T0_to_T1 5a2cc52214 add a comment b1dc74233b re-instate broken vectorization 79925fce8f Merge remote-tracking branch 'romick-knmi/romick-knmi_CY49T0_HIRLAM_tech_contrib_3' into accord_CY49T0_to_T1 bdab329b73 Merge branch 'honnert_CY49_TURB3D' into accord_CY49T0_to_T1 0115fffb54 turb optimizations 74b59a5680 Merge remote-tracking branch 'SebastienRietteMTO/riette_CY49T1PHYEX_bf' into accord_CY49T0_to_T1 0a2bccab62 Merge remote-tracking branch 'RachelHonnert/honnert_CY49_LS18' into accord_CY49T0_to_T1 3d39e3c59d Merge remote-tracking branch 'YannSeity/seity_CY49_from48t1op' into accord_CY49T0_to_T1 23276ac69b remove LMNHLEONARD 8b0cdd0bb0 Merge branch 'accord_CY49T0_bf' into accord_CY49T0_to_T1 be9be2a97c S. Riette comments 1e52eb3fc4 rename LBL89ARP into LBL89TOP and optimize associated code ee310cefe5 Merge remote-tracking branch 'gh-SebastienRietteMTO/riette_CY49T0PHYEX_optim' into accord_CY49T0_to_T1 64b53509ed bug in fypp 93a4e77a6c bug interface structure 81124b3226 Set-up modifications 5c4ba693bb compilation bugs cdf9bcc9cf TURB3D modifications with PHYEX and GPU 07301cbd08 LGZ already in accord_CY49T0_to_T1/ change comment cb4d2aef56 modifs compilation et ajout LBL89ARP dans NAM_TURBn et bl89.F90 bb470178ce merge with the following 48t1_op arome branches : seity_CY48T1_op0.04%physArome seity_CY48T1_op1.02_cheks_apl_arome seity_CY48T1_op1.03_bf_STAT_arome seity_CY48T1_op1.04_LITOTC seity_CY48T1_op1.06_bfLITOTC seity_CY48T1_op1.05_LimitTendQ 3c73227c02 Merge remote-tracking branch 'gh-SebastienRietteMTO/riette_CY49T0_PHYEX_integration' into mf_CY49T0_raps --- aux/ini_phyex.F90 | 5 +- aux/modd_aerosol_prop.F90 | 85 ++++ aux/modd_nrt_aerosols.F90 | 61 +++ aux/modd_nsv.F90 | 8 +- aux/modi_aermr_nc.F90 | 30 ++ aux/modi_aerosol_process.F90 | 38 ++ aux/modi_ini_phyex.F90 | 3 +- conv/convect_chem_transport.F90 | 11 +- conv/convect_closure.F90 | 26 +- conv/convect_closure_adjust_shal.F90 | 14 +- conv/convect_closure_shal.F90 | 396 ++++++++------- conv/convect_closure_thrvlcl.F90 | 54 ++- conv/convect_downdraft.F90 | 18 +- conv/convect_mixing_funct.F90 | 30 +- conv/convect_satmixratio.F90 | 95 +--- conv/convect_satmixratio.h | 99 +++- conv/convect_trigger_funct.F90 | 14 +- conv/convect_trigger_shal.F90 | 149 +++--- conv/convect_updraft.F90 | 27 +- conv/convect_updraft_shal.F90 | 93 ++-- conv/ini_convpar.h | 93 +--- micro/aermr_nc.F90 | 291 +++++++++++ micro/aerosol_process.F90 | 354 ++++++++++++++ micro/ice_adjust.F90 | 22 +- micro/ini_aerosols_cams.F90 | 637 +++++++++++++++++++++++++ micro/mode_ini_snow.F90 | 2 +- micro/mode_rain_ice_old_nucleation.F90 | 7 + micro/modi_rain_ice_old.F90 | 6 + micro/rain_ice_old.F90 | 23 +- turb/modd_turbn.F90 | 17 +- turb/mode_bl89.F90 | 22 +- turb/mode_compute_updraft.F90 | 2 +- turb/modi_turb.F90 | 11 +- turb/turb.F90 | 90 +++- 34 files changed, 2223 insertions(+), 610 deletions(-) create mode 100644 aux/modd_aerosol_prop.F90 create mode 100644 aux/modd_nrt_aerosols.F90 create mode 100644 aux/modi_aermr_nc.F90 create mode 100644 aux/modi_aerosol_process.F90 create mode 100644 micro/aermr_nc.F90 create mode 100644 micro/aerosol_process.F90 create mode 100644 micro/ini_aerosols_cams.F90 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_aerosol_prop.F90 b/aux/modd_aerosol_prop.F90 new file mode 100644 index 00000000..33bb750a --- /dev/null +++ b/aux/modd_aerosol_prop.F90 @@ -0,0 +1,85 @@ +! ######spl + MODULE MODD_AEROSOL_PROP +! ########################## +! +!!**** *MODD_AEROSOL_PROP* - declaration of aerosol properties. +!! +!! PURPOSE +!! ------- +! The purpose of this declarative module is to declare ... +! +!! +!!** IMPLICIT ARGUMENTS +!! ------------------ +!! None +!! +!! REFERENCE +!! --------- +!! +!! AUTHOR +!! ------ +!! +!! MODIFICATIONS +!! ------------- +!! Original 04/12/17 +!! +!------------------------------------------------------------------------------- +! +!* 0. DECLARATIONS +! ------------ +! +IMPLICIT NONE + +! CAMS Aerosols : 14 species +! MOCAGE Aerosols : ??? + +INTEGER, PARAMETER :: NAEROMAX=16 + +REAL,SAVE :: XCONC_MIN ! Minimun droplet concentration + +INTEGER, SAVE :: NCCN +INTEGER, SAVE :: NIFN +INTEGER, SAVE :: NCOAR +INTEGER, SAVE :: NSEASALT +INTEGER, SAVE :: NCOARSEASALT +INTEGER, SAVE :: NDUST +INTEGER, SAVE :: NDRYDEP + +! Cloud condensation nuclei +INTEGER, DIMENSION(NAEROMAX), SAVE :: XACCNI +! Ice freezing Nuclei +INTEGER, DIMENSION(NAEROMAX), SAVE :: XIFNI +! Coarse modes for dry sedimentation +INTEGER, DIMENSION(NAEROMAX), SAVE :: XCOARSE +! Sea salt species +INTEGER, DIMENSION(NAEROMAX), SAVE :: XSEASALT +INTEGER, DIMENSION(NAEROMAX), SAVE :: XCOARSEASALT +! Desert dust +INTEGER, DIMENSION(NAEROMAX), SAVE :: XDUST +! Dry deposition active +INTEGER, DIMENSION(NAEROMAX), SAVE :: XDRYDEP + +! Aerosol distribution: log normal size distribution +! number mode radius(rm) +! geometric standard deviation(s) +REAL, DIMENSION(NAEROMAX), SAVE :: XMRAE, XSDAE +! Variables related with the bin limits +REAL, DIMENSION(NAEROMAX), SAVE :: XBINDOWN,XBINUP +REAL, DIMENSION(NAEROMAX), SAVE :: XERFDOWN,XERFUP,XR3CORR + +! Physical properties +REAL, DIMENSION(NAEROMAX), SAVE :: XKHYGROS ! hygroscopic factor +REAL, DIMENSION(NAEROMAX), SAVE :: XRHOAE ! aerosol density (Kg m-3) +! Optical properties +! Mass extinction +REAL, DIMENSION(NAEROMAX,12), SAVE :: XMEXT ! mass extinction + +REAL, DIMENSION(NAEROMAX,2), SAVE :: XVDRY ! Dry Deposition Velocity + ! over sea and land +REAL, DIMENSION(NAEROMAX), SAVE :: XFRICS ! Fracion of aerosols + ! in aquaeous phase + ! for in-cloud scavenging + +REAL, DIMENSION(NAEROMAX), SAVE :: XINTSS ! Sea salt emision + +END MODULE MODD_AEROSOL_PROP diff --git a/aux/modd_nrt_aerosols.F90 b/aux/modd_nrt_aerosols.F90 new file mode 100644 index 00000000..ea0542ec --- /dev/null +++ b/aux/modd_nrt_aerosols.F90 @@ -0,0 +1,61 @@ +! ######spl + MODULE MODD_NRT_AEROSOLS +! ##################### +! +!!**** *MODD_NRT_AEROSOLS* - declaration of the control parameters for +!! the use of near real time aerosols +!! +!! PURPOSE +!! ------- +!! The purpose of this declarative module is to define the set of +!! parameters for the use of nrt aerosols in the microphysics. +!! +!! +!!** IMPLICIT ARGUMENTS +!! ------------------ +!! None +!! +!! REFERENCE +!! --------- +!! +!! +!! AUTHOR +!! ------ +!! D. Martin-Perez +!! +!! MODIFICATIONS +!! ------------- +!! Original 05/08/2022 +!! +!------------------------------------------------------------------------------- +! +!* 0. DECLARATIONS +! ------------ +! +IMPLICIT NONE +! +LOGICAL, SAVE :: LAEIFN ! switch to activate ice nuclei +LOGICAL, SAVE :: LAECCN2CLDR ! switch to convert CCN into Cloud + ! droplets +LOGICAL, SAVE :: LAERDRDEP ! Switch for dry deposition activation +LOGICAL, SAVE :: LAERSSEM ! Switch for sea salt emission. +REAL, SAVE :: XSSMINLO ! Minimum supersaturation at lowest level +REAL, SAVE :: XSSMINUP ! Minimum supersaturation at high levels +REAL, SAVE :: XSSMAX ! Maximum supersaturation +REAL, SAVE :: XSSHEIGHT ! Height limit for SSMIN separation +REAL, SAVE :: XSSFACVV ! Factor for SS dependence with vert. velo. +REAL, SAVE :: XSSFACSS ! Factor for SS dependence with coarse sea salt. +REAL, SAVE :: XPANMIN ! Minimum particle concentration number + ! (m-3) +REAL, SAVE :: XCCNMIN ! Minimum cloud condensation number + ! concentration (m-3) +REAL, SAVE :: XCLDROPMIN ! Minimum cloud droplet number + ! concentration (m-3) +REAL, SAVE :: XIFNMINSIZE ! Minimum size of aerosol to be considered + ! ice nuclei (micrometers) + +! +!------------------------------------------------------------------------------- +! +END MODULE MODD_NRT_AEROSOLS + 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_aermr_nc.F90 b/aux/modi_aermr_nc.F90 new file mode 100644 index 00000000..83833e83 --- /dev/null +++ b/aux/modi_aermr_nc.F90 @@ -0,0 +1,30 @@ +! ######spl + MODULE MODI_AERMR_NC +! #################### +! +INTERFACE + SUBROUTINE AERMR_NC(PHYEX, PTSTEP, PRHODREF, PEXNREF, PPABST, & + & PCLDFR, PTHT, PRVT, PRCT, KNAERO, & + & PAER_MR, PSSAT, PAER_NC, PCCN_NC, PIFN_NC) +! +USE MODD_PHYEX, ONLY: PHYEX_t +! +TYPE(PHYEX_t), INTENT(IN) :: PHYEX +REAL, INTENT(IN) :: PTSTEP +INTEGER, INTENT(IN) :: KNAERO +REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODREF +REAL, DIMENSION(:,:,:), INTENT(IN) :: PEXNREF +REAL, DIMENSION(:,:,:), INTENT(IN) :: PPABST +REAL, DIMENSION(:,:,:), INTENT(IN) :: PCLDFR +REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHT +REAL, DIMENSION(:,:,:), INTENT(IN) :: PRVT +REAL, DIMENSION(:,:,:), INTENT(IN) :: PRCT +REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PAER_MR +REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PSSAT +REAL, DIMENSION(:,:,:,:), INTENT( OUT) :: PAER_NC +REAL, DIMENSION(:,:,:,:), INTENT( OUT) :: PCCN_NC +REAL, DIMENSION(:,:,:), INTENT( OUT) :: PIFN_NC +! +END SUBROUTINE AERMR_NC +END INTERFACE +END MODULE MODI_AERMR_NC diff --git a/aux/modi_aerosol_process.F90 b/aux/modi_aerosol_process.F90 new file mode 100644 index 00000000..4c252119 --- /dev/null +++ b/aux/modi_aerosol_process.F90 @@ -0,0 +1,38 @@ +! ######spl + MODULE MODI_AEROSOL_PROCESS +! #################### +! +INTERFACE + SUBROUTINE AEROSOL_PROCESS(PHYEX, KKA,KKU,KKL,KNAERO,PTSTEP, & + & PPABST,PTHT,PDZZ,PRHODREF,PRCT,PRIT, & + & PCLDFR,PFPR, & + & PLSM,PUGST,PVGST, & + & PAERMR,PAERTEND) +! +USE MODD_PHYEX, ONLY: PHYEX_t +! +TYPE(PHYEX_t), INTENT(IN) :: PHYEX +INTEGER, INTENT(IN) :: KKA +INTEGER, INTENT(IN) :: KKU +INTEGER, INTENT(IN) :: KKL +INTEGER, INTENT(IN) :: KNAERO +REAL, INTENT(IN) :: PTSTEP ! Double Time step + ! (single if cold + ! start) +REAL, DIMENSION(:,:,:), INTENT(IN) :: PPABST +REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHT +REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ ! Layer thikness (m) +REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODREF! Reference density +REAL, DIMENSION(:,:,:), INTENT(IN) :: PRCT ! Cloud water m.r. at t +REAL, DIMENSION(:,:,:), INTENT(IN) :: PRIT ! Pristine ice m.r. at t +REAL, DIMENSION(:,:,:), INTENT(IN) :: PCLDFR ! Convective Mass Flux Cloud fraction +REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PFPR ! upper-air precipitation fluxes +REAL, DIMENSION(:), INTENT(IN) :: PLSM ! Land/Sea Mask +REAL, DIMENSION(:), INTENT(IN) :: PUGST ! 10m u-wind gust +REAL, DIMENSION(:), INTENT(IN) :: PVGST ! 10m v-wind gust +REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PAERMR ! Aerosol mixing ratio +REAL, DIMENSION(:,:,:,:), INTENT( OUT) :: PAERTEND! Aerosol tendencies: scavenging + emission +END SUBROUTINE AEROSOL_PROCESS +END INTERFACE +END MODULE MODI_AEROSOL_PROCESS + 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/aermr_nc.F90 b/micro/aermr_nc.F90 new file mode 100644 index 00000000..99804021 --- /dev/null +++ b/micro/aermr_nc.F90 @@ -0,0 +1,291 @@ +! ######spl + SUBROUTINE AERMR_NC ( PHYEX, PTSTEP, PRHODREF, PEXNREF,PPABST, & + PCLDFR, PTHT, PRVT, PRCT, KNAERO, & + PAER_MR, PSSAT, PAER_NC, PCCN_NC, PIFN_NC) + + USE YOMHOOK, ONLY : LHOOK, DR_HOOK, JPHOOK +! ########################################################################## +! +!!**** * - Calculate the number of condensation nuclei +!! +!! PURPOSE +!! ------- +!! Calculate the number of nuclei from the mixing ratio +!! +!!** METHOD +!! ------ +!! Considering a log normal distribution for the particles of every specie +!! the mean cubic radius is used for the caltulation. +!! +!! EXTERNAL +!! -------- +!! +!! IMPLICIT ARGUMENTS +!! ------------------ +!! Module MODD_PARAMETERS +!! JPHEXT : Horizontal external points number +!! JPVEXT : Vertical external points number +!! Module MODD_CST +!! XPI +!! XP00 ! Reference pressure +!! XRD,XRV ! Gaz constant for dry air, vapor +!! XMD,XMV ! Molecular weight for dry air, vapor +!! XCPD ! Cpd (dry air) +!! XCL ! Cl (liquid) +!! XCI ! Ci (solid) +!! XTT ! Triple point temperature +!! XLVTT ! Vaporization heat constant +!! XALPW,XBETAW,XGAMW ! Constants for saturation vapor pressure +!! function over liquid water +!! XALPI,XBETAI,XGAMI ! Constants for saturation vapor pressure +!! function over solid ice +!! AUTHOR +!! ------ +!! D. Martin-Perez +!! +!! MODIFICATIONS +!! ------------- +!! Original 03.2018 +!! 02.2021: Inclusion of IFN +!! 28.06.2022: Modification of SSat to include effect of coarse Sea Salt +!! +!------------------------------------------------------------------------------- +! +!* 0. DECLARATIONS +! ------------ +! +USE MODD_CST +USE MODD_PARAMETERS +USE MODD_PHYEX, ONLY: PHYEX_t +USE MODD_AEROSOL_PROP +USE MODD_NRT_AEROSOLS, ONLY: XSSMINLO,XSSFACSS,XIFNMINSIZE,XCLDROPMIN +USE MODI_GAMMA +! +IMPLICIT NONE +! +!* 0.1 Declarations of dummy arguments : +! +TYPE(PHYEX_t), INTENT(IN) :: PHYEX +REAL, INTENT(IN) :: PTSTEP ! Time step +REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODREF! Reference density +REAL, DIMENSION(:,:,:), INTENT(IN) :: PEXNREF ! Reference Exner function +REAL, DIMENSION(:,:,:), INTENT(IN) :: PPABST ! absolute pressure at t +REAL, DIMENSION(:,:,:), INTENT(IN) :: PCLDFR ! Convective Mass Flux Cloud fraction +REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHT ! Theta at time t +REAL, DIMENSION(:,:,:), INTENT(IN) :: PRVT ! Water vapor m.r. at t +REAL, DIMENSION(:,:,:), INTENT(IN) :: PRCT ! Cloud water m.r. at t +! +INTEGER, INTENT(IN) :: KNAERO ! Number of aerosol species +REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PAER_MR ! Aerosol mixing ratio (kg/kg) +REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PSSAT ! Supersaturation +REAL, DIMENSION(:,:,:,:), INTENT( OUT) :: PAER_NC ! Aerosol number concentration (m-3) +REAL, DIMENSION(:,:,:,:), INTENT( OUT) :: PCCN_NC ! Cloud condensation nuclei number concentration (m-3) per species +REAL, DIMENSION(:,:,:), INTENT( OUT) :: PIFN_NC ! Ice freezing nuclei number concentration (m-3) + +! ****** Used for aerosol . **** +! Weighted concentration +REAL, DIMENSION(SIZE(PEXNREF,1),SIZE(PEXNREF,2),SIZE(PEXNREF,3),KNAERO) :: & + ZCONC_AER ! Weighted concentration +REAL, DIMENSION(SIZE(PEXNREF,1),SIZE(PEXNREF,2),SIZE(PEXNREF,3)) :: & + ZT, & ! Temperature + ZAER_NCCOARSS, ZAER_NCSS, ZAER_NCDU +! +! Working variables +! +REAL :: ZZT ! Temperature +REAL :: ZZZSSW ! Supersaturation +REAL :: ZACTAE ! Activated nuclei +REAL :: ZPPABST ! pressure +! +REAL :: ZDKA, ZDKB ! Kohler parameters (a,b) +REAL :: ZDKR ! Critical radius +REAL :: ZDKRMIN ! Minimum aerosol radius for saturation +REAL :: ZRTINCLOUD + +! ****** end aerosol ***** +! +!* 0.2 Declarations of local variables : +! +REAL :: ZCCNAE ! Condensation nuclei +REAL :: ZAERPR, ZR3AE +REAL :: ZERFRAE, ZERFIFN +REAL :: ZWVDIF,ZRA,ZRD,ZW1,ZSSC,ZDELTASS + +INTEGER :: IC,JC +INTEGER :: JI,JJ,JK +INTEGER :: IIB,IIE,IJE,IJB +INTEGER :: IKT,IKTB,IKTE +INTEGER :: JCCN + +! +!------------------------------------------------------------------------------- +! +!* 1.1 COMPUTE THE LOOP BOUNDS +! ----------------------- +! +REAL(KIND=JPHOOK) :: ZHOOK_HANDLE +IF (LHOOK) CALL DR_HOOK('AERMR_NC',0,ZHOOK_HANDLE) + +IIB=1+JPHEXT +IIE=SIZE(PRHODREF,1) - JPHEXT +IJB=1+JPHEXT +IJE=SIZE(PRHODREF,2) - JPHEXT +IKT=SIZE(PRHODREF,3) +IKTB=1+JPVEXT +IKTE=IKT-JPVEXT + +! +!* 1.2 COMPUTE SOME CONSTANT PARAMETERS +! + +! Initialize Variables + +ZDKRMIN=2.0E-8 ! Minimum aerosol radius to be activated (m) +ZERFIFN=0.0 + +PAER_NC(:,:,:,:)=0.0 +PCCN_NC(:,:,:,:)=0.0 +PIFN_NC(:,:,:)=0.0 + +ZCONC_AER(:,:,:,:)=0.0 + +ZRTINCLOUD=PHYEX%RAIN_ICE_DESCRN%XRTMIN(2) + +! Temperature +ZT(:,:,:) = PTHT(:,:,:) * ( PPABST(:,:,:) / XP00 ) ** (XRD/XCPD) + +! Compute the number of aerosol nuclei from the mixing ratio. +DO JC=1,KNAERO + ! rm**3 exp(4.5 (ln s)**2 ) + ZR3AE=XR3CORR(JC)*1.0E-18 + ZAERPR=XRHOAE(JC)*ZR3AE + PAER_NC(:,:,:,JC)=(0.75*PAER_MR(:,:,:,JC)*PRHODREF(:,:,:))/(XPI*ZAERPR) +ENDDO + +! Total number of aerosol nuclei per specie +ZCONC_AER(:,:,:,:)=MAX(PAER_NC(:,:,:,:),0.0) + +! Estimation of the reduction of the supersaturation due to the presence coarse sea salt +ZWVDIF=0.211E-4 ! Water vapor diffusivity (m**2/s) +ZRA=5.0E-6 ! reference aerosol radius (m) +ZAER_NCCOARSS(:,:,:)=0.0 +ZDELTASS=0.03E-2 +IF ( NCOARSEASALT > 0 ) THEN + DO JK = IKTB, IKTE + DO JI = IIB,IIE + DO JJ = IJE,IJB + IF (PSSAT(JI,JJ,JK)>0.0) THEN + ZZT = ZT(JI,JJ,JK) + ZDKA=0.32538E-6/ZZT + DO IC=1,NCOARSEASALT + JC=XCOARSEASALT(IC) + ZAER_NCCOARSS(JI,JJ,JK)=ZAER_NCCOARSS(JI,JJ,JK)+ & + & ZCONC_AER(JI,JJ,JK,JC) + ZDKB=XKHYGROS(JC)*(ZRA)**3.0 + ENDDO + ! Critical radius and supersaturation + ZRD=(3.0*ZDKB/ZDKA)**0.5 + ZSSC=(4.0*ZDKA**3.0/(27.0*ZDKB))**0.5 + ZW1=MAX(0.0,4.0*XPI*ZWVDIF*(XSSFACSS*PTSTEP)*ZRD*((XSSMINLO-ZDELTASS)-ZSSC)) + PSSAT(JI,JJ,JK)=PSSAT(JI,JJ,JK)-MIN(ZDELTASS, & + & ZW1*(MIN(XCLDROPMIN,ZAER_NCCOARSS(JI,JJ,JK)))) + ENDIF + ENDDO + ENDDO + ENDDO +ENDIF + +! ***************************************************** +! Compute total number of IFN + +DO IC=1,NIFN + JC=XIFNI(IC) + IF ( JC > KNAERO ) CYCLE + ! Number of c.n. bigger than XIFNMINSIZE micrometers + ZERFIFN=ERFPDF(XIFNMINSIZE,XMRAE(JC),XSDAE(JC)) + PIFN_NC(:,:,:)=PIFN_NC(:,:,:)+PAER_NC(:,:,:,JC)*(XERFUP(JC)-ZERFIFN)/(XERFUP(JC)-XERFDOWN(JC)) +ENDDO + + +! ***************************************************** +! Compute Active CCN +! + +! Number of aerosol species to become CCN + + +! 1. Consider grid points with Cloud water content greater than ZRTINCLOUD + +DO JK = IKTB,IKTE + DO JI = IIB,IIE + DO JJ = IJE,IJB + IF ( PRCT(JI,JJ,JK)>ZRTINCLOUD ) THEN + ZZT = ZT(JI,JJ,JK) + ZZZSSW=PSSAT(JI,JJ,JK) + ZPPABST=PPABST(JI,JJ,JK) + + ZDKA=0.32538/ZZT + +! 2. Calculate activated cloud condensation nuclei + DO JC=1,NCCN + ZACTAE=0.0 + JCCN=XACCNI(JC) + IF ( JCCN > KNAERO ) CYCLE + ZCCNAE=0.0 + ZERFRAE=0.0 + ZDKR=0.0 + ZDKB=XKHYGROS(JCCN) + ZDKR=MAX(ZDKA/(6.75*ZDKB)**(0.33)/& + & (ZZZSSW)**(0.66),ZDKRMIN) + ZCCNAE=ZCONC_AER(JI,JJ,JK,JCCN)/(XERFUP(JCCN)-XERFDOWN(JCCN)) + ZERFRAE=MAX(XERFDOWN(JCCN),0.5*& + & erf(0.7071*LOG(ZDKR/XMRAE(JCCN))/LOG(XSDAE(JCCN)))) + ZACTAE=MAX(ZCCNAE*(XERFUP(JCCN)-ZERFRAE),0.0) +! Total cloud condensation nuclei + PCCN_NC(JI,JJ,JK,JC)=MAX(ZACTAE,0.0) + ENDDO + ENDIF + ENDDO + ENDDO +ENDDO + +! Total aerosol nuclei +PAER_NC(:,:,:,:)=MAX(PAER_NC(:,:,:,:),0.0) + +! End of calculation of Active CCN +!------------------------------------------------------------------------------- +! +IF (LHOOK) CALL DR_HOOK('AERMR_NC',1,ZHOOK_HANDLE) +! +! +!------------------------------------------------------------------------------- +! +! +CONTAINS + +! +! + +FUNCTION ERFPDF(PRAD,PMODRAD,PSIGMA) RESULT (PERFPDF) +! +! 0.5*erf(log(PRAD/Rg)/(sqrt(2.0))*log(sigma) +! +! + IMPLICIT NONE +! + REAL :: PRAD + REAL :: PMODRAD + REAL :: PSIGMA + REAL :: PERFPDF +! + REAL(KIND=JPHOOK) :: ZHOOK_HANDLE + IF (LHOOK) CALL DR_HOOK('AERMR_NC:ERFPDF',0,ZHOOK_HANDLE) +! + PERFPDF = 0.5*erf(0.7071*LOG(PRAD/PMODRAD)/LOG(PSIGMA)) +! + IF (LHOOK) CALL DR_HOOK('AERMR_NC:ERFPDF',1,ZHOOK_HANDLE) +! +END FUNCTION ERFPDF +! +! +END SUBROUTINE AERMR_NC diff --git a/micro/aerosol_process.F90 b/micro/aerosol_process.F90 new file mode 100644 index 00000000..b9a485ac --- /dev/null +++ b/micro/aerosol_process.F90 @@ -0,0 +1,354 @@ +! ######spl + SUBROUTINE AEROSOL_PROCESS(PHYEX,KKA,KKU,KKL,KNAERO,PTSTEP,& + & PPABST,PTHT,PDZZ,PRHODREF,PRCT,PRIT,& + & PCLDFR,PFPR, & + & PLSM,PUGST,PVGST, & + & PAERMR,PAERTEND) + + USE YOMHOOK , ONLY : LHOOK, DR_HOOK, JPHOOK +! ########################################################################## +! +!! Aerosol removal processes +!! +!! PURPOSE +!! ------- +!! Calculate the scavenging and deposition of aerosols. +!! +!! IMPLICIT ARGUMENTS +!! ------------------ +!! Module MODD_PARAMETERS +!! JPHEXT : Horizontal external points number +!! JPVEXT : Vertical external points number +!! Module MODD_CST +!! XPI +!! XG : Gravity constant +!! XP00 : Reference pressure +!! XRD : Gaz constant for dry air +!! XCPD : Cpd (dry air) +!! XRHOWL : density of liquid water +!! XRHOLI : density of ice water +!! Module MODD_RAIN_ICE_DESCR +!! XRTMIN : Min values allowed for mixing ratios +!! Module MODD_AEROSOL_PROP +!! XVDRY : Dry Deposition Velocity +!! XFRICS : Fracion of aerosols in aquaeous phase +!! : for in-cloud scavenging +!! +!! REFERENCE +!! --------- +!! Mocrette et al. 2009 and Redy, 2019. +!! +!! AUTHOR +!! ------ +!! Daniel Martin Perez (AEMET) +!! +!! MODIFICATIONS +!! ------------- +!! Original: 2018 +!! D. Martin : 27.03.2023 : Adition of sea salt emission. +!! D. Martin : 28.06.2022 : Modification of scavenging, addition of +!! snow fluxes. +!! +!------------------------------------------------------------------------- +!! +!* 0. DECLARATIONS +! ------------ + +USE MODD_PARAMETERS, ONLY : JPHEXT, JPVEXT +USE MODD_CST, ONLY : XG, XPI, XRHOLW, XRHOLI, XP00, XRD, XCPD +USE MODD_PHYEX, ONLY: PHYEX_t +USE MODD_AEROSOL_PROP, ONLY : XVDRY, XFRICS, NCOAR, XCOARSE, XSDAE, XRHOAE, XMRAE, & + & XERFUP, XERFDOWN, XKHYGROS, NDRYDEP, XDRYDEP, NSEASALT, XSEASALT, & + & NDUST, XDUST, XINTSS +USE MODD_NRT_AEROSOLS, ONLY : LAERDRDEP, LAERSSEM +! +IMPLICIT NONE +! +TYPE(PHYEX_t), INTENT(IN) :: PHYEX +INTEGER, INTENT(IN) :: KKA ! near ground array index +INTEGER, INTENT(IN) :: KKU ! uppest atmosphere array index +INTEGER, INTENT(IN) :: KKL ! vert. levels type 1=MNH -1=ARO +INTEGER, INTENT(IN) :: KNAERO ! number of aerosol species +REAL, INTENT(IN) :: PTSTEP ! Double Time step + ! (single if cold start) +REAL, DIMENSION(:,:,:), INTENT(IN) :: PPABST ! Pressure +REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHT ! Theta +REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ ! Layer thickness (m) +REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODREF ! Reference density +REAL, DIMENSION(:,:,:), INTENT(IN) :: PRCT ! Cloud water m.r. at t +REAL, DIMENSION(:,:,:), INTENT(IN) :: PRIT ! Pristine ice m.r. at t +REAL, DIMENSION(:,:,:), INTENT(IN) :: PCLDFR ! Convective Mass Flux Cloud fraction +REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PFPR ! upper-air precipitation fluxes +REAL, DIMENSION(:), INTENT(IN) :: PLSM ! Land/Sea Mask +REAL, DIMENSION(:), INTENT(IN) :: PUGST ! 10m u-wind gust +REAL, DIMENSION(:), INTENT(IN) :: PVGST ! 10m v-wind gust +REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PAERMR +REAL, DIMENSION(:,:,:,:), INTENT( OUT) :: PAERTEND ! Aerosol tendencies: emission + scavenging coefficient rate +! +!* 0.2 declaration of local variables +! +INTEGER :: IIB ! Define the domain where is +INTEGER :: IIE ! the microphysical sources have to be computed +INTEGER :: IJB ! +INTEGER :: IJE ! +INTEGER :: IKB ! +INTEGER :: IKE ! +INTEGER :: JC, IC +INTEGER :: JK, JI, JJ + +REAL, DIMENSION(SIZE(PRHODREF,1),SIZE(PRHODREF,2),SIZE(PRHODREF,3),KNAERO) :: & + ZW1, ZW2, ZW3, ZW4, & ! work arrays + ZAER_MR, & + ZAER_DEP ! Aerosol deposition +REAL, DIMENSION(SIZE(PRHODREF,1),SIZE(PRHODREF,2),SIZE(PRHODREF,3)) :: & + ZT, & ! Temperature + ZVELSED, & ! Velocity sedimentation + ZMFREEPATH, & ! mean free path + ZCUNN, & ! Cunningham factor for sedimentation + ZPFPRR, ZPFPRS + +REAL, DIMENSION(SIZE(PRHODREF,1),3) :: & + ZWS_SS ! work array + +! Sink and source at surface of aerosol in kg/m2 +REAL, DIMENSION(SIZE(PRHODREF,1)) :: & + ZAERSNKDU, & ! Desert dust sink + ZAERSNKSS, & ! Sea salt sink + ZAERSRCSS ! Sea salt source + +! XRTMIN = Minimum value for the mixing ratio +REAL :: ZINVTSTEP +REAL :: ZPFMIN ! Minimum Flux precipitation +REAL :: ZEF_RAIN, ZEF_SNOW +REAL :: ZRAY_RAIN, ZRAY_SNOW +REAL :: ZCLFK +REAL :: ZVDRY +REAL :: ZVISC +REAL :: ZAIRMAS +REAL :: ZBOLTZMANN +REAL :: ZAER_R +REAL :: Z10WIND + +REAL :: ZDKA, ZDKB, ZDKR, ZDKRMIN, ZFRICS, ZERFRAE + +!------------------------------------------------------------------------------- +! +!* 1. Parameters for aerosol wet deposition +REAL(KIND=JPHOOK) :: ZHOOK_HANDLE +IF (LHOOK) CALL DR_HOOK('AEROSOL_PROCESS',0,ZHOOK_HANDLE) + + +IIB=1+JPHEXT +IIE=SIZE(PDZZ,1) - JPHEXT +IJB=1+JPHEXT +IJE=SIZE(PDZZ,2) - JPHEXT +IKB=KKA+JPVEXT*KKL +IKE=KKU-JPVEXT*KKL +! +!* 1.2 COMPUTE SOME CONSTANT PARAMETERS +! +ZINVTSTEP=1./PTSTEP + +! +!* 2. compute the fluxes +! + +!ZPFMIN=1.0E-12 +ZPFMIN=0.0 +ZEF_RAIN=0.001 +ZEF_SNOW=0.01 +ZRAY_RAIN=1.0E-3 +ZRAY_SNOW=1.0E-3 +ZVISC=1.78E-5 ! Air viscosity (Pa*s) +ZAIRMAS=28.97E-3 ! Air mean molar mass(Kg*mol-1) +ZBOLTZMANN=8.31445 ! (J*K-1*mol-1) +ZDKRMIN=2.0E-8 + +ZAER_MR(:,:,:,:)=MAX(PAERMR(:,:,:,:),0.0) + +ZW1(:,:,:,:) = 0. +ZW2(:,:,:,:) = 0. +ZW3(:,:,:,:) = 0. +ZW4(:,:,:,:) = 0. + +ZWS_SS(:,:)=0.0 + +ZT(:,:,:) = PTHT(:,:,:) * ( PPABST(:,:,:) / XP00 ) ** (XRD/XCPD) + +ZAERSNKDU(:)=0.0 +ZAERSNKSS(:)=0.0 +ZAERSRCSS(:)=0.0 + +! ************************************************ +! AEROSOL REMOVAL PROCESSES +! ************************************************ + +! Rain and snow fluxes initialized with sedimentation fluxes + +DO JI = IIB,IIE + DO JJ = IJE,IJB + DO JK = IKE+1 , IKB, -1*KKL + ZPFPRR(JI,JJ,JK)=MAX(PFPR(JI,JJ,JK,3),0.0) + ZPFPRS(JI,JJ,JK)=MAX(PFPR(JI,JJ,JK,5),0.0) + ENDDO + ENDDO +ENDDO + + +! Aerosol wet deposition +! (Loop from top to bottom KKL=-1 ARO) +DO JK = IKE+1 , IKB, -1*KKL + DO JI = IIB,IIE + DO JJ = IJE,IJB + ! inside cloud + IF (ZPFPRR(JI,JJ,JK) > ZPFMIN .OR. & + & ZPFPRS(JI,JJ,JK) > ZPFMIN) THEN + ! inside cloud + ZCLFK=PCLDFR(JI,JJ,JK) + ZDKA=0.32538/ZT(JI,JJ,JK) + IF ((PRCT(JI,JJ,JK) > PHYEX%RAIN_ICE_DESCRN%XRTMIN(2) .OR. & + & PRIT(JI,JJ,JK) > PHYEX%RAIN_ICE_DESCRN%XRTMIN(4)) .AND. ZCLFK > 0.1) THEN + DO JC=1,KNAERO + ZW1(JI,JJ,JK,JC)=((ZPFPRR(JI,JJ,JK+KKL)-ZPFPRR(JI,JJ,JK))+& + & (ZPFPRS(JI,JJ,JK+KKL)-ZPFPRS(JI,JJ,JK)))/ & + & (PRCT(JI,JJ,JK)+PRIT(JI,JJ,JK))/ & + & (PRHODREF(JI,JJ,JK)*PDZZ(JI,JJ,JK)*ZCLFK) + ZDKB=XKHYGROS(JC) + IF ( ZDKB == 0.0 ) THEN + ZFRICS=XFRICS(JC) + ELSE + ZDKR=MAX(ZDKA/(6.75*ZDKB)**(0.33)/(0.05E-2)**(0.66),ZDKRMIN) + ZERFRAE=MAX(XERFDOWN(JC),0.5*erf(0.7071*LOG(ZDKR/XMRAE(JC))/LOG(XSDAE(JC)))) + ZFRICS=MIN((XERFUP(JC)-ZERFRAE)/(XERFUP(JC)-XERFDOWN(JC)),XFRICS(JC)) + ENDIF + ZW1(JI,JJ,JK,JC)=ZCLFK*ZFRICS*ZAER_MR(JI,JJ,JK,JC)*MIN(MAX(-1.0,ZW1(JI,JJ,JK,JC)),0.0) + ENDDO + ELSE + ! below cloud (rain) + DO JC=1,KNAERO + ZW2(JI,JJ,JK,JC)=MIN(-0.75*ZPFPRR(JI,JJ,JK)*ZEF_RAIN/& + & (ZRAY_RAIN*XRHOLW)- & + & 0.75*ZPFPRS(JI,JJ,JK)*ZEF_SNOW/(ZRAY_SNOW*XRHOLI),0.0)* & + & ZAER_MR(JI,JJ,JK,JC) + ENDDO + ENDIF + ENDIF + ENDDO + ENDDO +ENDDO + +! Aerosol dry deposition +! Surface layer +! +IF (LAERDRDEP) THEN + DO JI = IIB,IIE + DO JJ = IJE,IJB + IF (PRCT(JI,JJ,IKB) < PHYEX%RAIN_ICE_DESCRN%XRTMIN(2) .OR. & + & (ZPFPRR(JI,JJ,IKB) < ZPFMIN .AND. & + & ZPFPRS(JI,JJ,IKB) < ZPFMIN)) THEN + DO IC=1,NDRYDEP ! only for those aerosol with parametrized emission + JC=XDRYDEP(IC) + IF ( JC > KNAERO ) CYCLE + ZVDRY=(1.0-PLSM(JI))*XVDRY(JC,1)+PLSM(JI)*XVDRY(JC,2) + ZW3(JI,JJ,IKB,JC)=-0.01*ZVDRY*ZAER_MR(JI,JJ,IKB,JC)/& + & PDZZ(JI,JJ,IKB) + ENDDO + ENDIF + ENDDO + ENDDO +ENDIF + +ZAER_DEP(:,:,:,:)=ZW1(:,:,:,:)+ZW2(:,:,:,:)+ZW3(:,:,:,:) + + +! Dry sedimentation only for coarse modes (Redy, 2019) +ZMFREEPATH(:,:,:)=ZVISC/PPABST(:,:,:)*SQRT(0.5*XPI*ZBOLTZMANN*ZT(:,:,:)/ZAIRMAS) + +DO IC=1,NCOAR + JC=XCOARSE(IC) + IF ( JC > KNAERO ) CYCLE + ZAER_R=XMRAE(JC)*EXP(0.5*(LOG(XSDAE(JC))**2))*1.0E-6 + ZCUNN(:,:,:)=1.0+ZMFREEPATH(:,:,:)/ZAER_R*& + & (1.257+0.4*EXP(-1.1*ZAER_R/ZMFREEPATH(:,:,:))) + ZVELSED(:,:,:)=MIN(0.22*(ZAER_R**2)*XRHOAE(JC)*XG*ZCUNN(:,:,:)/ZVISC, & + & PDZZ(:,:,:)*ZINVTSTEP) + + DO JK = IKE+1 , IKB, -1*KKL + ZW4(:,:,JK,JC)=(ZAER_MR(:,:,JK+KKL,JC)*ZVELSED(:,:,JK+KKL)-& + & ZAER_MR(:,:,JK,JC)*ZVELSED(:,:,JK))/PDZZ(:,:,JK) + ENDDO + +ENDDO + +! Sum all the aerosol removal contributions + +PAERTEND(:,:,:,:)=MAX(ZAER_DEP(:,:,:,:)+ZW4(:,:,:,:),-0.8*ZAER_MR(:,:,:,:)*ZINVTSTEP) + +! Aerosol sinks on the surface (kg/m2/s) + +DO JI = IIB,IIE + ZAERSNKDU(JI)=0.0 + DO JJ = IJE,IJB + DO JK = IKE+1 , IKB, -1*KKL + DO IC=1,NSEASALT ! Sea Salt species + JC=XSEASALT(IC) + ! diagnostics of all aerosol deposition in kg/m2/s + ZAERSNKSS(JI)=ZAERSNKSS(JI)+(MIN(ZW1(JI,JJ,JK,JC),0.0)+ & + & MIN(ZW2(JI,JJ,JK,JC),0.0))*PRHODREF(JI,JJ,JK)*& + & PDZZ(JI,JJ,JK) + ENDDO + DO IC=1,NDUST ! Dust species + JC=XDUST(IC) + ! diagnostics of all aerosol deposition in kg/m2/s + ZAERSNKDU(JI)=ZAERSNKDU(JI)+(MIN(ZW1(JI,JJ,JK,JC),0.0)+ & + & MIN(ZW2(JI,JJ,JK,JC),0.0))*PRHODREF(JI,JJ,JK)*& + & PDZZ(JI,JJ,JK) + ENDDO + ENDDO + DO IC=1,NSEASALT ! Sea Salt species + JC=XSEASALT(IC) + ZAERSNKSS(JI)=ZAERSNKSS(JI)+(MIN(ZW3(JI,JJ,IKB,JC),0.0)+& + & MIN(ZW4(JI,JJ,IKB,JC),0.0))*& + & PRHODREF(JI,JJ,IKB)*PDZZ(JI,JJ,IKB) + ENDDO + DO IC=1,NDUST ! Dust species + JC=XDUST(IC) + ZAERSNKDU(JI)=ZAERSNKDU(JI)+(MIN(ZW3(JI,JJ,IKB,JC),0.0)+& + & MIN(ZW4(JI,JJ,IKB,JC),0.0))*& + & PRHODREF(JI,JJ,IKB)*PDZZ(JI,JJ,IKB) + ENDDO + ENDDO + ZAERSNKDU(JI)=MIN(ZAERSNKDU(JI),0.0) + ZAERSNKSS(JI)=MIN(ZAERSNKSS(JI),0.0) +ENDDO + + +!********************************* +! AEROSOL EMISSIONS +!********************************* +! +! Sea salt production (Monahan,1986) +! +IF (LAERSSEM) THEN + DO JI = IIB,IIE + ZAERSRCSS(JI)=0.0 + IF ( PLSM(JI)<0.5 ) THEN + Z10WIND=SQRT(PUGST(JI)**2+PVGST(JI)**2) + DO IC=1,NSEASALT ! Sea Salt species + JC=XSEASALT(IC) + ZWS_SS(JI,JC)=1.373*Z10WIND**(3.41)*XINTSS(JC)* & + & (XRHOAE(JC)/8.0)*(4.0*XPI/3.0)*1.0E-18 + ZAERSRCSS(JI)=ZAERSRCSS(JI)+ZWS_SS(JI,JC) + DO JJ = IJE,IJB + PAERTEND(JI,JJ,IKB,JC)=PAERTEND(JI,JJ,IKB,JC)+ZWS_SS(JI,JC)/& + & (PRHODREF(JI,JJ,IKB)*PDZZ(JI,JJ,IKB)) + ENDDO + ENDDO + ENDIF + ENDDO +ENDIF + + +IF (LHOOK) CALL DR_HOOK('AEROSOL_PROCESS',1,ZHOOK_HANDLE) +! +END SUBROUTINE AEROSOL_PROCESS 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/ini_aerosols_cams.F90 b/micro/ini_aerosols_cams.F90 new file mode 100644 index 00000000..8382b45d --- /dev/null +++ b/micro/ini_aerosols_cams.F90 @@ -0,0 +1,637 @@ +! ######spl + SUBROUTINE INI_AEROSOLS_CAMS + USE YOMHOOK, ONLY : LHOOK, DR_HOOK, JPHOOK +! ########################################################### +! +!!**** *INI_AEROSOLS_CAMS* - initialize the constants necessary +!! for cams aerosols +!! +!! PURPOSE +!! ------- +!! The purpose of this routine is to initialize the constants used to +!! resolve the mixed phase microphysical scheme. +!! AEROSOLS FROM CAMS: A total number of 14 species. +!! +!!** METHOD +!! ------ +!! The constants are initialized to their numerical values and the number +!! +!! EXTERNAL +!! -------- +!! +!! +!! IMPLICIT ARGUMENTS +!! ------------------ +!! Module MODD_AEROSOL_PROP +!! +!! REFERENCE +!! --------- +!! +!! AUTHOR +!! ------ +!! +!! MODIFICATIONS +!! ------------- +!! Original D. Martin-Perez 11/02/2019 +!! +!------------------------------------------------------------------------------- +! +!* 0. DECLARATIONS +! ------------ +! +USE MODD_AEROSOL_PROP +! +IMPLICIT NONE +! +!* 0.1 Declarations of dummy arguments : +! +! +! +!* 0.2 Declarations of local variables : +! +! +LOGICAL :: GFLAG ! Logical flag for printing the constatnts on the output + ! listing +!------------------------------------------------------------------------------- +! +! +REAL(KIND=JPHOOK) :: ZHOOK_HANDLE +IF (LHOOK) CALL DR_HOOK('INI_AEROSOLS_CAMS',0,ZHOOK_HANDLE) +! +! +!------------------------------------------------------------------------------- +! +!* 1. CLOUD CONDENSATION NUMBER +! ------------------------- +! +! XACCNI: Index of aerosol species considered condensation nuclei + + +XACCNI(:)=0 +XIFNI(:)=0 +XCOARSE(:)=0 +XSEASALT(:)=0 +XDUST(:)=0 +XDRYDEP(:)=0 +XINTSS(:)=0 + +! Cloud condensation nuclei species +NCCN=9 ! Number of species becoming CCN +XACCNI(1)=1 ! Sea salt 1 +XACCNI(2)=2 ! Sea salt 2 +XACCNI(3)=3 ! Sea salt 3 +XACCNI(4)=8 ! Hydrophilic Organic Matter +XACCNI(5)=10 ! Hydrophilic Black Carbon +XACCNI(6)=11 ! Sulfate +XACCNI(7)=12 ! Nitrate fine mode +XACCNI(8)=13 ! Nitrate coarse mode +XACCNI(9)=14 ! Ammonium + +! Ice nuclei species +NIFN=5 ! Number of species becoming IFN +XIFNI(1)=4 ! Desert dust 1 +XIFNI(2)=5 ! Desert dust 2 +XIFNI(3)=6 ! Desert dust 3 +XIFNI(4)=7 ! Hydrophobic OM +XIFNI(5)=9 ! Hydrophobic BC + +! For supercoarse sea salt and dust (Remy 2019) +NCOAR=2 +XCOARSE(1)=3 +XCOARSE(2)=6 + +! Sea salt +NSEASALT=3 +XSEASALT(1)=1 +XSEASALT(2)=2 +XSEASALT(3)=3 +! Coarse Sea salt +NCOARSEASALT=2 +XCOARSEASALT(1)=2 +XCOARSEASALT(2)=3 + +! Desert Dust +NDUST=3 +XDUST(1)=4 +XDUST(2)=5 +XDUST(3)=6 + +! Dry Deposition activated +NDRYDEP=6 +XDRYDEP(1)=1 +XDRYDEP(2)=2 +XDRYDEP(3)=3 +XDRYDEP(4)=4 +XDRYDEP(5)=5 +XDRYDEP(6)=6 + + +! +!* 2. AEROSOL PHYSICAL PROPERTIES +! --------------------------- +! + +! +! Optical properties +! +! XMEXT: Mass Extinction (m2 kg-1) +! +! Aerosol distribution: log normal size distribution +! +! XMRAE: Modal Radius(micrometers),Rg +! XSDAE: geometric Standard Ddeviation,sigma +! +! Limits of the error function +! +! XBINDOWN: Lower bin limit (micrometers) +! XBINUP: Upper bin limit (micrometers) +! +! XERFDOWN(JC)=0.5*erf(log({DOWNlimit}/Rg)/(sqrt(2.0))*log(sigma) +! XERFUP(JC)=0.5*erf(log({UPlimit}/Rg)/(sqrt(2.0))*log(sigma) +! +! Mean cubic radius of the distribution: +! = Rg**3*((3.0/sqrt(2.0)))*log(sigma))**2 +! X = (3.0/sqrt(2.0)))*log(sigma) +! {ERF3DOWN,ERF3UP} = +! 0.5*erf(log({DOWNlimit,Uplimit}-X)/Rg)/(sqrt(2.0))*log(sigma) +! Correction of the mean cubic radius when the bin size limits are considered +! +! XR3CORR(JC)= *(ERF3UP-ERF3DOWN)/(XERFUP(JC)-XERFDOWN(JC) +! +! Physical properties +! +! XRHOAE: mass density (Kg m-3) +! XKHYGROS: Hygroscopic parameter +! +! XMEXT: Mass extinction +! +! Variables for removal processes +! +! XVDRY: Dry Deposition Velocity (cm s-1) (Remy et al., 2019) +! XFRICS: Fraction of aerosols in aqueous phase for in-cloud scavenging +! (Remy et al., 2019) + +XMEXT(:,:)=0.0 +XMRAE(:)=0.0 +XSDAE(:)=0.0 +XBINDOWN(:)=0.0 +XBINUP(:)=0.0 +XERFDOWN(:)=0.0 +XERFUP(:)=0.0 +XR3CORR(:)=0.0 +XRHOAE(:)=0.0 +XKHYGROS(:)=0.0 +XVDRY(:,:)=0.0 ! XVDRY= Dry deposition over (sea,land) +XFRICS(:)=0.0 + +!------------------------------------------------- +! 2.1 Sea salt +!------------------------------------------------- +! +! Sea salt aerosol density (Kg m-3) +! 2160 (dry particles) +! 1182 (sea salt production asuming 80% humidity, Mocrette et al. 2009) + +! Factor for Sea Salt emision +! Monahan, 1986 +! INTEGRAL (1+0.057*r**1.05)10**(1.19*exp(-B**2)) +! B=(0.380-log(r))/0.650 +! r: aerosol radius in um + +XINTSS(1)=2.09913 ! r: 2x0.03 - 2x0.5 um +XINTSS(2)=29.6649 ! r: 2x0.5 - 2x5.0 um +XINTSS(3)=80.3799 ! r: 2x5.0 - 2x20.0 um + +! +! Sea salt, size bin limits: 0.03 - 0.5 microm. +! Size distribution parameters: +XMRAE(1)=0.1992 +XSDAE(1)= 1.92 +XBINDOWN(1)=0.03 +XBINUP(1)=0.5 +XERFDOWN(1)=ERFPDF(XBINDOWN(1),XMRAE(1),XSDAE(1)) ! -0.4981 +XERFUP(1)=ERFPDF(XBINUP(1),XMRAE(1),XSDAE(1)) ! 0.4208 +! Mean cubic radius of the distribution: =5.3640e-02 +! Corrected mean cubic radius +XR3CORR(1)=1.7071e-02 +! Density (kg/m3): +XRHOAE(1)= 2160.0 + +XKHYGROS(1)=1.28 + +XVDRY(1,:)=(/0.05,0.2/) +XFRICS(1)=0.7 + +! Mass extinction +XMEXT(1,:)=(/827.98,755.72,684.17,614.08,2197.81,2612.42,3223.86,4048.34,5381.25,6506.83,8497.96,13330.20/) + +! +! Sea salt, size bin limits: 0.5 - 5.0 microm. +! + +! Size distribution parameters: +XMRAE(2)=1.9920 +XSDAE(2)= 2.00 +XBINDOWN(2)=0.5 +XBINUP(2)=5.0 +XERFDOWN(2)=ERFPDF(XBINDOWN(2),XMRAE(2),XSDAE(2)) ! -0.4769 +XERFUP(2)=ERFPDF(XBINUP(2),XMRAE(2),XSDAE(2)) ! 0.4079 +! Mean cubic radius of the distribution: =6.8680e+01 +! Corrected mean cubic radius +XR3CORR(2)=1.7549e+01 +! Density (kg/m3): +XRHOAE(2)=2160.0 + +XKHYGROS(2)=1.28 + +XVDRY(2,:)=(/0.05,0.5/) +XFRICS(2)=0.9 + +! Mass extinction +XMEXT(2,:)=(/143.79,143.36,143.30,143.72,285.42,330.74,376.77,432.98,520.54,590.13,713.12,1030.59/) + +! +! Sea salt, size bin limits: 5.0 - 20.0 microm. +! +! Size distribution parameters: +XMRAE(3)=1.9920 +XSDAE(3)=2.00 +XBINDOWN(3)=5.0 +XBINUP(3)=20.0 +XERFDOWN(3)=ERFPDF(XBINDOWN(3),XMRAE(3),XSDAE(3)) ! 0.4079 +XERFUP(3)=ERFPDF(XBINUP(3),XMRAE(3),XSDAE(3)) ! 0.4996 +! Mean cubic radius of the distribution: =6.8680e+01 +! Corrected mean cubic radius +XR3CORR(3)=5.0026e+02 +! Density (kg/m3): +XRHOAE(3)=2160.0 + +XKHYGROS(3)=1.28 + +XVDRY(3,:)=(/1.20,1.50/) +XFRICS(3)=0.9 + +! Mass extinction +XMEXT(3,:)=(/38.65,38.85,38.85,38.63,79.27,92.17,105.23,122.77,149.16,171.24,209.44,309.58/) + +!------------------------------------------------- +! 2.2 Desert Dust +!------------------------------------------------- + +! +! Desert Dust, size bin limits: 0.03 - 0.55 microm. +! +! Size distribution parameters: +XMRAE(4)=0.2900 +XSDAE(4)=2.00 +XBINDOWN(4)=0.03 +XBINUP(4)=0.55 +XERFDOWN(4)=ERFPDF(XBINDOWN(4),XMRAE(4),XSDAE(4)) ! -0.4995 +XERFUP(4)=ERFPDF(XBINUP(4),XMRAE(4),XSDAE(4)) ! 0.3221 +! Mean cubic radius of the distribution: =2.1191e-01 +! Corrected mean cubic radius +XR3CORR(4)=3.1940e-02 +! Density (kg/m3): +XRHOAE(4)=2610.0 + +XKHYGROS(4)=0.0 + +XVDRY(4,:)=(/0.02,0.02/) +XFRICS(4)=0.3 + +! Mass extinction +XMEXT(4,:)=2496.68 + +! +! Desert Dust, size bin limits: 0.55 - 0.9 microm. +! + +! Size distribution parameters: +XMRAE(5)=0.2900 +XSDAE(5)=2.00 +XBINDOWN(5)=0.55 +XBINUP(5)=0.9 +XERFDOWN(5)=ERFPDF(XBINDOWN(5),XMRAE(5),XSDAE(5)) ! 0.3221 +XERFUP(5)=ERFPDF(XBINUP(5),XMRAE(5),XSDAE(5)) ! 0.4489 +! Mean cubic radius of the distribution: =2.1191e-01 +! Corrected mean cubic radius +XR3CORR(5)=3.4124e-01 +! Density (kg/m3): +XRHOAE(5)=2610.0 + +XKHYGROS(5)=0.0 + +XVDRY(5,:)=(/0.2,0.2/) +XFRICS(5)=0.4 + +! Mass extinction +XMEXT(5,:)=955.078 + +! +! Desert Dust, size bin limits: 0.9 - 20.0 microm. +! + +! Size distribution parameters: +XMRAE(6)=0.2900 +XSDAE(6)=2.00 +XBINDOWN(6)=0.9 +XBINUP(6)=20.0 +XERFDOWN(6)=ERFPDF(XBINDOWN(6),XMRAE(6),XSDAE(6)) ! 0.4489 +XERFUP(6)=ERFPDF(XBINUP(6),XMRAE(6),XSDAE(6)) ! 0.5000 +! Mean cubic radius of the distribution: =2.1191e-01 +! Corrected mean cubic radius +XR3CORR(6)=2.7845e+00 +! Density (kg/m3): +XRHOAE(6)=2610.0 + +XKHYGROS(6)=0.0 + +XVDRY(6,:)=(/1.2,1.2/) +XFRICS(6)=0.4 + +! Mass extinction +XMEXT(6,:)=406.533 + +!------------------------------------------------- +! 2.3 Organic matter +!------------------------------------------------- +! + +! +! 7:Hydrophobic Organic Matter +! + +! Size distribution parameters: +XMRAE(7)=0.02 +XSDAE(7)=2.2 +XBINDOWN(7)=0.005 +XBINUP(7)=20.0 +XERFDOWN(7)=ERFPDF(XBINDOWN(7),XMRAE(7),XSDAE(7)) ! -0.4996 +XERFUP(7)=ERFPDF(XBINUP(7),XMRAE(7),XSDAE(7)) ! 0.4996 +! Mean cubic radius of the distribution: =1.0861e-03 +! Corrected mean cubic radius +XR3CORR(7)=9.7073e-04 +! Density (kg/m3): +XRHOAE(7)=1825.0 + +XKHYGROS(7)=0.0 + +XVDRY(7,:)=(/0.1,0.1/) +XFRICS(7)=0.0 + +! Mass extinction (RH=80) +XMEXT(7,:)=2321.03 + +! +! 8:Hydrophilic Organic Matter +! +! Size distribution parameters: +XMRAE(8)=0.02 +XSDAE(8)=2.2 +XBINDOWN(8)=0.005 +XBINUP(8)=20.0 +XERFDOWN(8)=ERFPDF(XBINDOWN(8),XMRAE(8),XSDAE(8)) ! -0.4996 +XERFUP(8)=ERFPDF(XBINUP(8),XMRAE(8),XSDAE(8)) ! 0.4996 +! Mean cubic radius of the distribution: =1.0861e-03 +! Corrected mean cubic radius +XR3CORR(8)=9.7073e-04 +! Density (kg/m3): +XRHOAE(8)=1825.0 + +XKHYGROS(8)=0.3 + +XVDRY(8,:)=(/0.1,0.1/) +XFRICS(8)=0.7 + +! Mass extinction +XMEXT(8,:)=(/1942.13,2131.58,2321.03,2510.48,2699.93,2889.38,3185.61,3481.84,4119.88,4909.06,5698.24,8206.75/) + +! +!------------------------------------------------- +! 2.4 Black Carbon +!------------------------------------------------- +! + +! Hydrophobic Black Carbon size bin limits 0.005 - 0.500 micrometers. + +! Size distribution (Oshima et al., 2009) +! XMRAE(9)=0.053; XSDAE(9)=1.5; XRHOAE(9)=1800. +! Size distribution parameters: +XMRAE(9)=0.0118 +XSDAE(9)=2.00 +XBINDOWN(9)=0.005 +XBINUP(9)=0.5 +XERFDOWN(9)=ERFPDF(XBINDOWN(9),XMRAE(9),XSDAE(9)) ! -0.3923 +XERFUP(9)=ERFPDF(XBINUP(9),XMRAE(9),XSDAE(9)) ! 0.5000 +! Mean cubic radius of the distribution: =1.4276e-05 +! Corrected mean cubic radius +XR3CORR(9)=1.5985e-05 +! Density (kg/m3): +XRHOAE(9)=1000.0 + +XKHYGROS(9)=0.0 + +XVDRY(9,:)=(/0.1,0.1/) +XFRICS(9)=0.0 + +XMEXT(9,:)=13487.8 + +! +! Hydrophilic Black Carbon size bin limits 0.005 - 0.500 micrometers. +! + + +! Size distribution parameters: +XMRAE(10)=0.0118 +XSDAE(10)=2.00 +XBINDOWN(10)=0.005 +XBINUP(10)=0.5 +XERFDOWN(10)=ERFPDF(XBINDOWN(10),XMRAE(10),XSDAE(10)) ! -0.3923 +XERFUP(10)=ERFPDF(XBINUP(10),XMRAE(10),XSDAE(10)) ! 0.5000 +! Mean cubic radius of the distribution: =1.4276e-05 +! Corrected mean cubic radius +XR3CORR(10)=1.5985e-05 +! Density (kg/m3): +XRHOAE(10)=1000.0 + +XKHYGROS(10)=0.1 + +XVDRY(10,:)=(/0.1,0.1/) +XFRICS(10)=0.7 + +! Mass extinction +XMEXT(10,:)=13487.8 + +! +!------------------------------------------------- +! 2.5. Sulphate properties (0.005-20 micrometers) +!------------------------------------------------- +! +! Sulfate (Vd=0.05 for oceans; Vd=0.25 all other surfaces ) + +! Size distribution parameters: +XMRAE(11)=0.0355 +XSDAE(11)=2.00 +XBINDOWN(11)=0.005 +XBINUP(11)=20.0 +XERFDOWN(11)=ERFPDF(XBINDOWN(11),XMRAE(11),XSDAE(11)) ! -0.4977 +XERFUP(11)=ERFPDF(XBINUP(11),XMRAE(11),XSDAE(11)) ! 0.5000 +! Mean cubic radius of the distribution: =3.8873e-04 +! Corrected mean cubic radius +XR3CORR(11)=3.8964e-04 +! Density (kg/m3): +XRHOAE(11)=1760.0 + +XKHYGROS(11)=0.6 + +XVDRY(11,:)=(/0.15,0.25/) +XFRICS(11)=0.7 + +! Mass extinction +XMEXT(11,:)=(/3119.62,3510.56,3901.51,4292.45,4683.40,5074.35,5685.64,6296.94,7613.58,9242.11,10870.64,16047.13/) + +! +!------------------------------------------------- +! 2.6. Nitrate properties (0.005-0.9 micrometers) +!------------------------------------------------- + +! +! Nitrate fine mode (NH4NO3), size bin limits: 0.005 - 0.9 microm. +! + +! Size distribution parameters: +XMRAE(12)=0.0355 +XSDAE(12)=2.0 +XBINDOWN(12)=0.005 +XBINUP(12)=0.9 +XERFDOWN(12)=ERFPDF(XBINDOWN(12),XMRAE(12),XSDAE(12)) ! -0.4977 +XERFUP(12)=ERFPDF(XBINUP(12),XMRAE(12),XSDAE(12)) ! 0.5000 +! Corrected mean cubic radius +XR3CORR(12)=3.8774e-04 +! Density (kg/m3): +XRHOAE(12)=1730.0 + +XKHYGROS(12)=0.64 + +XVDRY(12,:)=(/0.15,0.15/) +XFRICS(12)=0.4 + +XMEXT(12,:)=(/3501.201,3501.201,3501.201,3501.201,4757.47,5354.868,6161.467,7361.848,13555.21,15649.4,14113.68,23958.8/) + +! +! Nitrate coarse mode ( NaNO3+Ca(NO3)2 ), size bin limits: 0.9 - 20.0 microm. +! + +! Size distribution parameters: +XMRAE(13)=1.992 +XSDAE(13)=2.0 +XBINDOWN(13)=0.9 +XBINUP(13)=20.0 +XERFDOWN(13)=ERFPDF(XBINDOWN(13),XMRAE(13),XSDAE(13)) ! -0.3741 +XERFUP(13)=ERFPDF(XBINUP(13),XMRAE(13),XSDAE(13)) ! 0.4996 +! Corrected mean cubic radius +XR3CORR(13)=7.0228e+01 +! Density (kg/m3): +XRHOAE(13)=1400.0 + +XKHYGROS(13)=0.9 + +XVDRY(13,:)=(/0.15,0.15/) +XFRICS(13)=0.4 + +! Mass extinction +XMEXT(13,:)=(/4295.583,4295.583,4295.583,4295.583,5150.058,6185.985,6780.04,7425.708,8129.261,10575.19,14709.95,26338.47/) + +! +!------------------------------------------------- +! 2.7. Ammonium properties (0.005-20 micrometers) +!------------------------------------------------- +! (Ammonium sulfate for Amonia) +! Ammonium (Vd=0.05 for oceans; Vd=0.25 all other surfaces ) + +! Size distribution parameters: +XMRAE(14)=0.0355 +XSDAE(14)=2.00 +XBINDOWN(14)=0.005 +XBINUP(14)=20.0 +XERFDOWN(14)=ERFPDF(XBINDOWN(14),XMRAE(14),XSDAE(14)) ! -0.4977 +XERFUP(14)=ERFPDF(XBINUP(14),XMRAE(14),XSDAE(14)) ! 0.5000 +! Corrected mean cubic radius +XR3CORR(14)=3.8964e-04 +! Density (kg/m3): +XRHOAE(14)=1760.0 + +XKHYGROS(14)=0.6 + +XVDRY(14,:)=(/0.15,0.15/) +XFRICS(14)=0.4 + +! Mass extinction +XMEXT(14,:)=(/190.3313,232.0307,274.9777,321.6228,346.3498,371.914,425.8234,483.4848,544.9835,610.332,751.631,906.7416/) + + +! +IF (LHOOK) CALL DR_HOOK('INI_AEROSOL_CAMS',1,ZHOOK_HANDLE) + +CONTAINS +! +!---------------------------------------- +! + +FUNCTION ERFPDF(PRAD,PMODRAD,PSIGMA) RESULT (PERFPDF) +! +! 0.5*erf(log(PRAD/Rg)/(sqrt(2.0)*log(sigma))) +! +! + IMPLICIT NONE +! + REAL :: PRAD + REAL :: PMODRAD + REAL :: PSIGMA + REAL :: PERFPDF +! + REAL(KIND=JPHOOK) :: ZHOOK_HANDLE + IF (LHOOK) CALL DR_HOOK('INI_AEROSOLS_CAMS:ERFPDF',0,ZHOOK_HANDLE) +! + PERFPDF = 0.5*erf(0.7071*LOG(PRAD/PMODRAD)/LOG(PSIGMA)) +! + IF (LHOOK) CALL DR_HOOK('INI_AEROSOLS_CAMS:ERFPDF',1,ZHOOK_HANDLE) +! +END FUNCTION ERFPDF + + +FUNCTION RAD3CORR(PMODRAD,PSIGMA,PBINDOWN,PBINUP) RESULT (PR3CORR) +! Mean cubic radius of the distribution: +! = Rg**3*((3.0/sqrt(2.0)))*log(sigma))**2 +! X = (3.0/sqrt(2.0)))*log(sigma) +! {ERF3DOWN,ERF3UP} = +! 0.5*erf(log({DOWNlimit,Uplimit}-X)/Rg)/(sqrt(2.0))*log(sigma) +! Correction of the mean cubic radius when the bin size limits are considered +! +! XR3CORR(JC)= *(ERF3UP-ERF3DOWN)/(XERFUP(JC)-XERFDOWN(JC) +! + IMPLICIT NONE +! + REAL :: PMODRAD,PSIGMA,PBINDOWN,PBINUP + REAL :: PR3CORR + REAL :: ZR3, ZX,ZWUP,ZWDOWN,ZW3UP,ZW3DOWN +! + REAL(KIND=JPHOOK) :: ZHOOK_HANDLE + IF (LHOOK) CALL DR_HOOK('INI_AEROSOLS_CAMS:RAD3CORR',0,ZHOOK_HANDLE) + + ZX=2.12132*LOG(PSIGMA) + + ZR3=(PMODRAD**3)*ZX**2 + + ZWUP=ERFPDF(PBINUP,PMODRAD,PSIGMA) + ZWDOWN=ERFPDF(PBINDOWN,PMODRAD,PSIGMA) + ZW3UP=ERFPDF(PBINUP-ZX,PMODRAD,PSIGMA) + ZW3DOWN=ERFPDF(PBINDOWN-ZX,PMODRAD,PSIGMA) + + PR3CORR=ZR3*(ZW3UP-ZW3DOWN)/(ZWUP-ZWDOWN) + + IF (LHOOK) CALL DR_HOOK('INI_AEROSOLS_CAMS:RAD3CORR',1,ZHOOK_HANDLE) + +END FUNCTION RAD3CORR + +END SUBROUTINE INI_AEROSOLS_CAMS 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..ac9fae31 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, PIFNNC, & PTHS, PRVS, PRIS, PCIT, & PICENU, PT, PZZZ, & PRHT) @@ -21,6 +22,7 @@ SUBROUTINE RAIN_ICE_OLD_NUCLEATION(D, CST, ICEP, KSIZE, OCND2, LMODICEDEP, KRR, USE MODD_CST, ONLY: CST_T USE MODE_TIWMX, ONLY: ESATI, ESATW, AM3, REDIN USE MODD_RAIN_ICE_PARAM_N,ONLY: RAIN_ICE_PARAM_T + USE MODD_NRT_AEROSOLS, ONLY: LAEIFN ! !* 0. DECLARATIONS ! ------------ @@ -52,6 +54,8 @@ 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 + 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.LAEIFN) 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..e28f162f 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, PCLDROP, PIFNNC, & TBUDGETS, KBUDGETS, & PICENU, PKGN_ACON, PKGN_SBGR, & PRHT, PRHS, PINPRH, PFPR) @@ -98,6 +99,11 @@ 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 +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..c6b9be1b 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, 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 @@ -181,6 +183,7 @@ SUBROUTINE RAIN_ICE_OLD (D, CST, PARAMI, ICEP, ICED, BUCONF, USE MODI_GAMMA, ONLY: GAMMA USE MODE_TIWMX, ONLY: ESATI, ESATW, AA2, BB3, AA2W, BB3W USE MODE_TIWMX_TAB, ONLY: TIWMX_TAB +USE MODD_NRT_AEROSOLS, ONLY : LAEIFN ! USE MODE_RAIN_ICE_OLD_NUCLEATION, ONLY: RAIN_ICE_OLD_NUCLEATION USE MODE_RAIN_ICE_OLD_SEDIMENTATION_STAT, ONLY: RAIN_ICE_OLD_SEDIMENTATION_STAT @@ -268,6 +271,11 @@ 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 +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, 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..c0516fc3 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 ! @@ -2052,7 +2080,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 +2133,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 +2159,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 From 9b5a2237de3585d80654a12826ad87e572e9c691 Mon Sep 17 00:00:00 2001 From: Alexandre MARY Date: Wed, 25 Sep 2024 13:29:48 +0000 Subject: [PATCH 2/3] remove aerosols stuff from phyex (back to mpa/), cf. commits 7ec04ee1 & 8d1b4852 of IAL --- aux/modd_aerosol_prop.F90 | 85 ---- aux/modd_nrt_aerosols.F90 | 61 --- aux/modi_aermr_nc.F90 | 30 -- aux/modi_aerosol_process.F90 | 38 -- micro/aermr_nc.F90 | 291 ----------- micro/aerosol_process.F90 | 354 -------------- micro/ini_aerosols_cams.F90 | 637 ------------------------- micro/mode_rain_ice_old_nucleation.F90 | 6 +- micro/modi_rain_ice_old.F90 | 3 +- micro/rain_ice_old.F90 | 6 +- 10 files changed, 8 insertions(+), 1503 deletions(-) delete mode 100644 aux/modd_aerosol_prop.F90 delete mode 100644 aux/modd_nrt_aerosols.F90 delete mode 100644 aux/modi_aermr_nc.F90 delete mode 100644 aux/modi_aerosol_process.F90 delete mode 100644 micro/aermr_nc.F90 delete mode 100644 micro/aerosol_process.F90 delete mode 100644 micro/ini_aerosols_cams.F90 diff --git a/aux/modd_aerosol_prop.F90 b/aux/modd_aerosol_prop.F90 deleted file mode 100644 index 33bb750a..00000000 --- a/aux/modd_aerosol_prop.F90 +++ /dev/null @@ -1,85 +0,0 @@ -! ######spl - MODULE MODD_AEROSOL_PROP -! ########################## -! -!!**** *MODD_AEROSOL_PROP* - declaration of aerosol properties. -!! -!! PURPOSE -!! ------- -! The purpose of this declarative module is to declare ... -! -!! -!!** IMPLICIT ARGUMENTS -!! ------------------ -!! None -!! -!! REFERENCE -!! --------- -!! -!! AUTHOR -!! ------ -!! -!! MODIFICATIONS -!! ------------- -!! Original 04/12/17 -!! -!------------------------------------------------------------------------------- -! -!* 0. DECLARATIONS -! ------------ -! -IMPLICIT NONE - -! CAMS Aerosols : 14 species -! MOCAGE Aerosols : ??? - -INTEGER, PARAMETER :: NAEROMAX=16 - -REAL,SAVE :: XCONC_MIN ! Minimun droplet concentration - -INTEGER, SAVE :: NCCN -INTEGER, SAVE :: NIFN -INTEGER, SAVE :: NCOAR -INTEGER, SAVE :: NSEASALT -INTEGER, SAVE :: NCOARSEASALT -INTEGER, SAVE :: NDUST -INTEGER, SAVE :: NDRYDEP - -! Cloud condensation nuclei -INTEGER, DIMENSION(NAEROMAX), SAVE :: XACCNI -! Ice freezing Nuclei -INTEGER, DIMENSION(NAEROMAX), SAVE :: XIFNI -! Coarse modes for dry sedimentation -INTEGER, DIMENSION(NAEROMAX), SAVE :: XCOARSE -! Sea salt species -INTEGER, DIMENSION(NAEROMAX), SAVE :: XSEASALT -INTEGER, DIMENSION(NAEROMAX), SAVE :: XCOARSEASALT -! Desert dust -INTEGER, DIMENSION(NAEROMAX), SAVE :: XDUST -! Dry deposition active -INTEGER, DIMENSION(NAEROMAX), SAVE :: XDRYDEP - -! Aerosol distribution: log normal size distribution -! number mode radius(rm) -! geometric standard deviation(s) -REAL, DIMENSION(NAEROMAX), SAVE :: XMRAE, XSDAE -! Variables related with the bin limits -REAL, DIMENSION(NAEROMAX), SAVE :: XBINDOWN,XBINUP -REAL, DIMENSION(NAEROMAX), SAVE :: XERFDOWN,XERFUP,XR3CORR - -! Physical properties -REAL, DIMENSION(NAEROMAX), SAVE :: XKHYGROS ! hygroscopic factor -REAL, DIMENSION(NAEROMAX), SAVE :: XRHOAE ! aerosol density (Kg m-3) -! Optical properties -! Mass extinction -REAL, DIMENSION(NAEROMAX,12), SAVE :: XMEXT ! mass extinction - -REAL, DIMENSION(NAEROMAX,2), SAVE :: XVDRY ! Dry Deposition Velocity - ! over sea and land -REAL, DIMENSION(NAEROMAX), SAVE :: XFRICS ! Fracion of aerosols - ! in aquaeous phase - ! for in-cloud scavenging - -REAL, DIMENSION(NAEROMAX), SAVE :: XINTSS ! Sea salt emision - -END MODULE MODD_AEROSOL_PROP diff --git a/aux/modd_nrt_aerosols.F90 b/aux/modd_nrt_aerosols.F90 deleted file mode 100644 index ea0542ec..00000000 --- a/aux/modd_nrt_aerosols.F90 +++ /dev/null @@ -1,61 +0,0 @@ -! ######spl - MODULE MODD_NRT_AEROSOLS -! ##################### -! -!!**** *MODD_NRT_AEROSOLS* - declaration of the control parameters for -!! the use of near real time aerosols -!! -!! PURPOSE -!! ------- -!! The purpose of this declarative module is to define the set of -!! parameters for the use of nrt aerosols in the microphysics. -!! -!! -!!** IMPLICIT ARGUMENTS -!! ------------------ -!! None -!! -!! REFERENCE -!! --------- -!! -!! -!! AUTHOR -!! ------ -!! D. Martin-Perez -!! -!! MODIFICATIONS -!! ------------- -!! Original 05/08/2022 -!! -!------------------------------------------------------------------------------- -! -!* 0. DECLARATIONS -! ------------ -! -IMPLICIT NONE -! -LOGICAL, SAVE :: LAEIFN ! switch to activate ice nuclei -LOGICAL, SAVE :: LAECCN2CLDR ! switch to convert CCN into Cloud - ! droplets -LOGICAL, SAVE :: LAERDRDEP ! Switch for dry deposition activation -LOGICAL, SAVE :: LAERSSEM ! Switch for sea salt emission. -REAL, SAVE :: XSSMINLO ! Minimum supersaturation at lowest level -REAL, SAVE :: XSSMINUP ! Minimum supersaturation at high levels -REAL, SAVE :: XSSMAX ! Maximum supersaturation -REAL, SAVE :: XSSHEIGHT ! Height limit for SSMIN separation -REAL, SAVE :: XSSFACVV ! Factor for SS dependence with vert. velo. -REAL, SAVE :: XSSFACSS ! Factor for SS dependence with coarse sea salt. -REAL, SAVE :: XPANMIN ! Minimum particle concentration number - ! (m-3) -REAL, SAVE :: XCCNMIN ! Minimum cloud condensation number - ! concentration (m-3) -REAL, SAVE :: XCLDROPMIN ! Minimum cloud droplet number - ! concentration (m-3) -REAL, SAVE :: XIFNMINSIZE ! Minimum size of aerosol to be considered - ! ice nuclei (micrometers) - -! -!------------------------------------------------------------------------------- -! -END MODULE MODD_NRT_AEROSOLS - diff --git a/aux/modi_aermr_nc.F90 b/aux/modi_aermr_nc.F90 deleted file mode 100644 index 83833e83..00000000 --- a/aux/modi_aermr_nc.F90 +++ /dev/null @@ -1,30 +0,0 @@ -! ######spl - MODULE MODI_AERMR_NC -! #################### -! -INTERFACE - SUBROUTINE AERMR_NC(PHYEX, PTSTEP, PRHODREF, PEXNREF, PPABST, & - & PCLDFR, PTHT, PRVT, PRCT, KNAERO, & - & PAER_MR, PSSAT, PAER_NC, PCCN_NC, PIFN_NC) -! -USE MODD_PHYEX, ONLY: PHYEX_t -! -TYPE(PHYEX_t), INTENT(IN) :: PHYEX -REAL, INTENT(IN) :: PTSTEP -INTEGER, INTENT(IN) :: KNAERO -REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODREF -REAL, DIMENSION(:,:,:), INTENT(IN) :: PEXNREF -REAL, DIMENSION(:,:,:), INTENT(IN) :: PPABST -REAL, DIMENSION(:,:,:), INTENT(IN) :: PCLDFR -REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHT -REAL, DIMENSION(:,:,:), INTENT(IN) :: PRVT -REAL, DIMENSION(:,:,:), INTENT(IN) :: PRCT -REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PAER_MR -REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PSSAT -REAL, DIMENSION(:,:,:,:), INTENT( OUT) :: PAER_NC -REAL, DIMENSION(:,:,:,:), INTENT( OUT) :: PCCN_NC -REAL, DIMENSION(:,:,:), INTENT( OUT) :: PIFN_NC -! -END SUBROUTINE AERMR_NC -END INTERFACE -END MODULE MODI_AERMR_NC diff --git a/aux/modi_aerosol_process.F90 b/aux/modi_aerosol_process.F90 deleted file mode 100644 index 4c252119..00000000 --- a/aux/modi_aerosol_process.F90 +++ /dev/null @@ -1,38 +0,0 @@ -! ######spl - MODULE MODI_AEROSOL_PROCESS -! #################### -! -INTERFACE - SUBROUTINE AEROSOL_PROCESS(PHYEX, KKA,KKU,KKL,KNAERO,PTSTEP, & - & PPABST,PTHT,PDZZ,PRHODREF,PRCT,PRIT, & - & PCLDFR,PFPR, & - & PLSM,PUGST,PVGST, & - & PAERMR,PAERTEND) -! -USE MODD_PHYEX, ONLY: PHYEX_t -! -TYPE(PHYEX_t), INTENT(IN) :: PHYEX -INTEGER, INTENT(IN) :: KKA -INTEGER, INTENT(IN) :: KKU -INTEGER, INTENT(IN) :: KKL -INTEGER, INTENT(IN) :: KNAERO -REAL, INTENT(IN) :: PTSTEP ! Double Time step - ! (single if cold - ! start) -REAL, DIMENSION(:,:,:), INTENT(IN) :: PPABST -REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHT -REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ ! Layer thikness (m) -REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODREF! Reference density -REAL, DIMENSION(:,:,:), INTENT(IN) :: PRCT ! Cloud water m.r. at t -REAL, DIMENSION(:,:,:), INTENT(IN) :: PRIT ! Pristine ice m.r. at t -REAL, DIMENSION(:,:,:), INTENT(IN) :: PCLDFR ! Convective Mass Flux Cloud fraction -REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PFPR ! upper-air precipitation fluxes -REAL, DIMENSION(:), INTENT(IN) :: PLSM ! Land/Sea Mask -REAL, DIMENSION(:), INTENT(IN) :: PUGST ! 10m u-wind gust -REAL, DIMENSION(:), INTENT(IN) :: PVGST ! 10m v-wind gust -REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PAERMR ! Aerosol mixing ratio -REAL, DIMENSION(:,:,:,:), INTENT( OUT) :: PAERTEND! Aerosol tendencies: scavenging + emission -END SUBROUTINE AEROSOL_PROCESS -END INTERFACE -END MODULE MODI_AEROSOL_PROCESS - diff --git a/micro/aermr_nc.F90 b/micro/aermr_nc.F90 deleted file mode 100644 index 99804021..00000000 --- a/micro/aermr_nc.F90 +++ /dev/null @@ -1,291 +0,0 @@ -! ######spl - SUBROUTINE AERMR_NC ( PHYEX, PTSTEP, PRHODREF, PEXNREF,PPABST, & - PCLDFR, PTHT, PRVT, PRCT, KNAERO, & - PAER_MR, PSSAT, PAER_NC, PCCN_NC, PIFN_NC) - - USE YOMHOOK, ONLY : LHOOK, DR_HOOK, JPHOOK -! ########################################################################## -! -!!**** * - Calculate the number of condensation nuclei -!! -!! PURPOSE -!! ------- -!! Calculate the number of nuclei from the mixing ratio -!! -!!** METHOD -!! ------ -!! Considering a log normal distribution for the particles of every specie -!! the mean cubic radius is used for the caltulation. -!! -!! EXTERNAL -!! -------- -!! -!! IMPLICIT ARGUMENTS -!! ------------------ -!! Module MODD_PARAMETERS -!! JPHEXT : Horizontal external points number -!! JPVEXT : Vertical external points number -!! Module MODD_CST -!! XPI -!! XP00 ! Reference pressure -!! XRD,XRV ! Gaz constant for dry air, vapor -!! XMD,XMV ! Molecular weight for dry air, vapor -!! XCPD ! Cpd (dry air) -!! XCL ! Cl (liquid) -!! XCI ! Ci (solid) -!! XTT ! Triple point temperature -!! XLVTT ! Vaporization heat constant -!! XALPW,XBETAW,XGAMW ! Constants for saturation vapor pressure -!! function over liquid water -!! XALPI,XBETAI,XGAMI ! Constants for saturation vapor pressure -!! function over solid ice -!! AUTHOR -!! ------ -!! D. Martin-Perez -!! -!! MODIFICATIONS -!! ------------- -!! Original 03.2018 -!! 02.2021: Inclusion of IFN -!! 28.06.2022: Modification of SSat to include effect of coarse Sea Salt -!! -!------------------------------------------------------------------------------- -! -!* 0. DECLARATIONS -! ------------ -! -USE MODD_CST -USE MODD_PARAMETERS -USE MODD_PHYEX, ONLY: PHYEX_t -USE MODD_AEROSOL_PROP -USE MODD_NRT_AEROSOLS, ONLY: XSSMINLO,XSSFACSS,XIFNMINSIZE,XCLDROPMIN -USE MODI_GAMMA -! -IMPLICIT NONE -! -!* 0.1 Declarations of dummy arguments : -! -TYPE(PHYEX_t), INTENT(IN) :: PHYEX -REAL, INTENT(IN) :: PTSTEP ! Time step -REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODREF! Reference density -REAL, DIMENSION(:,:,:), INTENT(IN) :: PEXNREF ! Reference Exner function -REAL, DIMENSION(:,:,:), INTENT(IN) :: PPABST ! absolute pressure at t -REAL, DIMENSION(:,:,:), INTENT(IN) :: PCLDFR ! Convective Mass Flux Cloud fraction -REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHT ! Theta at time t -REAL, DIMENSION(:,:,:), INTENT(IN) :: PRVT ! Water vapor m.r. at t -REAL, DIMENSION(:,:,:), INTENT(IN) :: PRCT ! Cloud water m.r. at t -! -INTEGER, INTENT(IN) :: KNAERO ! Number of aerosol species -REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PAER_MR ! Aerosol mixing ratio (kg/kg) -REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PSSAT ! Supersaturation -REAL, DIMENSION(:,:,:,:), INTENT( OUT) :: PAER_NC ! Aerosol number concentration (m-3) -REAL, DIMENSION(:,:,:,:), INTENT( OUT) :: PCCN_NC ! Cloud condensation nuclei number concentration (m-3) per species -REAL, DIMENSION(:,:,:), INTENT( OUT) :: PIFN_NC ! Ice freezing nuclei number concentration (m-3) - -! ****** Used for aerosol . **** -! Weighted concentration -REAL, DIMENSION(SIZE(PEXNREF,1),SIZE(PEXNREF,2),SIZE(PEXNREF,3),KNAERO) :: & - ZCONC_AER ! Weighted concentration -REAL, DIMENSION(SIZE(PEXNREF,1),SIZE(PEXNREF,2),SIZE(PEXNREF,3)) :: & - ZT, & ! Temperature - ZAER_NCCOARSS, ZAER_NCSS, ZAER_NCDU -! -! Working variables -! -REAL :: ZZT ! Temperature -REAL :: ZZZSSW ! Supersaturation -REAL :: ZACTAE ! Activated nuclei -REAL :: ZPPABST ! pressure -! -REAL :: ZDKA, ZDKB ! Kohler parameters (a,b) -REAL :: ZDKR ! Critical radius -REAL :: ZDKRMIN ! Minimum aerosol radius for saturation -REAL :: ZRTINCLOUD - -! ****** end aerosol ***** -! -!* 0.2 Declarations of local variables : -! -REAL :: ZCCNAE ! Condensation nuclei -REAL :: ZAERPR, ZR3AE -REAL :: ZERFRAE, ZERFIFN -REAL :: ZWVDIF,ZRA,ZRD,ZW1,ZSSC,ZDELTASS - -INTEGER :: IC,JC -INTEGER :: JI,JJ,JK -INTEGER :: IIB,IIE,IJE,IJB -INTEGER :: IKT,IKTB,IKTE -INTEGER :: JCCN - -! -!------------------------------------------------------------------------------- -! -!* 1.1 COMPUTE THE LOOP BOUNDS -! ----------------------- -! -REAL(KIND=JPHOOK) :: ZHOOK_HANDLE -IF (LHOOK) CALL DR_HOOK('AERMR_NC',0,ZHOOK_HANDLE) - -IIB=1+JPHEXT -IIE=SIZE(PRHODREF,1) - JPHEXT -IJB=1+JPHEXT -IJE=SIZE(PRHODREF,2) - JPHEXT -IKT=SIZE(PRHODREF,3) -IKTB=1+JPVEXT -IKTE=IKT-JPVEXT - -! -!* 1.2 COMPUTE SOME CONSTANT PARAMETERS -! - -! Initialize Variables - -ZDKRMIN=2.0E-8 ! Minimum aerosol radius to be activated (m) -ZERFIFN=0.0 - -PAER_NC(:,:,:,:)=0.0 -PCCN_NC(:,:,:,:)=0.0 -PIFN_NC(:,:,:)=0.0 - -ZCONC_AER(:,:,:,:)=0.0 - -ZRTINCLOUD=PHYEX%RAIN_ICE_DESCRN%XRTMIN(2) - -! Temperature -ZT(:,:,:) = PTHT(:,:,:) * ( PPABST(:,:,:) / XP00 ) ** (XRD/XCPD) - -! Compute the number of aerosol nuclei from the mixing ratio. -DO JC=1,KNAERO - ! rm**3 exp(4.5 (ln s)**2 ) - ZR3AE=XR3CORR(JC)*1.0E-18 - ZAERPR=XRHOAE(JC)*ZR3AE - PAER_NC(:,:,:,JC)=(0.75*PAER_MR(:,:,:,JC)*PRHODREF(:,:,:))/(XPI*ZAERPR) -ENDDO - -! Total number of aerosol nuclei per specie -ZCONC_AER(:,:,:,:)=MAX(PAER_NC(:,:,:,:),0.0) - -! Estimation of the reduction of the supersaturation due to the presence coarse sea salt -ZWVDIF=0.211E-4 ! Water vapor diffusivity (m**2/s) -ZRA=5.0E-6 ! reference aerosol radius (m) -ZAER_NCCOARSS(:,:,:)=0.0 -ZDELTASS=0.03E-2 -IF ( NCOARSEASALT > 0 ) THEN - DO JK = IKTB, IKTE - DO JI = IIB,IIE - DO JJ = IJE,IJB - IF (PSSAT(JI,JJ,JK)>0.0) THEN - ZZT = ZT(JI,JJ,JK) - ZDKA=0.32538E-6/ZZT - DO IC=1,NCOARSEASALT - JC=XCOARSEASALT(IC) - ZAER_NCCOARSS(JI,JJ,JK)=ZAER_NCCOARSS(JI,JJ,JK)+ & - & ZCONC_AER(JI,JJ,JK,JC) - ZDKB=XKHYGROS(JC)*(ZRA)**3.0 - ENDDO - ! Critical radius and supersaturation - ZRD=(3.0*ZDKB/ZDKA)**0.5 - ZSSC=(4.0*ZDKA**3.0/(27.0*ZDKB))**0.5 - ZW1=MAX(0.0,4.0*XPI*ZWVDIF*(XSSFACSS*PTSTEP)*ZRD*((XSSMINLO-ZDELTASS)-ZSSC)) - PSSAT(JI,JJ,JK)=PSSAT(JI,JJ,JK)-MIN(ZDELTASS, & - & ZW1*(MIN(XCLDROPMIN,ZAER_NCCOARSS(JI,JJ,JK)))) - ENDIF - ENDDO - ENDDO - ENDDO -ENDIF - -! ***************************************************** -! Compute total number of IFN - -DO IC=1,NIFN - JC=XIFNI(IC) - IF ( JC > KNAERO ) CYCLE - ! Number of c.n. bigger than XIFNMINSIZE micrometers - ZERFIFN=ERFPDF(XIFNMINSIZE,XMRAE(JC),XSDAE(JC)) - PIFN_NC(:,:,:)=PIFN_NC(:,:,:)+PAER_NC(:,:,:,JC)*(XERFUP(JC)-ZERFIFN)/(XERFUP(JC)-XERFDOWN(JC)) -ENDDO - - -! ***************************************************** -! Compute Active CCN -! - -! Number of aerosol species to become CCN - - -! 1. Consider grid points with Cloud water content greater than ZRTINCLOUD - -DO JK = IKTB,IKTE - DO JI = IIB,IIE - DO JJ = IJE,IJB - IF ( PRCT(JI,JJ,JK)>ZRTINCLOUD ) THEN - ZZT = ZT(JI,JJ,JK) - ZZZSSW=PSSAT(JI,JJ,JK) - ZPPABST=PPABST(JI,JJ,JK) - - ZDKA=0.32538/ZZT - -! 2. Calculate activated cloud condensation nuclei - DO JC=1,NCCN - ZACTAE=0.0 - JCCN=XACCNI(JC) - IF ( JCCN > KNAERO ) CYCLE - ZCCNAE=0.0 - ZERFRAE=0.0 - ZDKR=0.0 - ZDKB=XKHYGROS(JCCN) - ZDKR=MAX(ZDKA/(6.75*ZDKB)**(0.33)/& - & (ZZZSSW)**(0.66),ZDKRMIN) - ZCCNAE=ZCONC_AER(JI,JJ,JK,JCCN)/(XERFUP(JCCN)-XERFDOWN(JCCN)) - ZERFRAE=MAX(XERFDOWN(JCCN),0.5*& - & erf(0.7071*LOG(ZDKR/XMRAE(JCCN))/LOG(XSDAE(JCCN)))) - ZACTAE=MAX(ZCCNAE*(XERFUP(JCCN)-ZERFRAE),0.0) -! Total cloud condensation nuclei - PCCN_NC(JI,JJ,JK,JC)=MAX(ZACTAE,0.0) - ENDDO - ENDIF - ENDDO - ENDDO -ENDDO - -! Total aerosol nuclei -PAER_NC(:,:,:,:)=MAX(PAER_NC(:,:,:,:),0.0) - -! End of calculation of Active CCN -!------------------------------------------------------------------------------- -! -IF (LHOOK) CALL DR_HOOK('AERMR_NC',1,ZHOOK_HANDLE) -! -! -!------------------------------------------------------------------------------- -! -! -CONTAINS - -! -! - -FUNCTION ERFPDF(PRAD,PMODRAD,PSIGMA) RESULT (PERFPDF) -! -! 0.5*erf(log(PRAD/Rg)/(sqrt(2.0))*log(sigma) -! -! - IMPLICIT NONE -! - REAL :: PRAD - REAL :: PMODRAD - REAL :: PSIGMA - REAL :: PERFPDF -! - REAL(KIND=JPHOOK) :: ZHOOK_HANDLE - IF (LHOOK) CALL DR_HOOK('AERMR_NC:ERFPDF',0,ZHOOK_HANDLE) -! - PERFPDF = 0.5*erf(0.7071*LOG(PRAD/PMODRAD)/LOG(PSIGMA)) -! - IF (LHOOK) CALL DR_HOOK('AERMR_NC:ERFPDF',1,ZHOOK_HANDLE) -! -END FUNCTION ERFPDF -! -! -END SUBROUTINE AERMR_NC diff --git a/micro/aerosol_process.F90 b/micro/aerosol_process.F90 deleted file mode 100644 index b9a485ac..00000000 --- a/micro/aerosol_process.F90 +++ /dev/null @@ -1,354 +0,0 @@ -! ######spl - SUBROUTINE AEROSOL_PROCESS(PHYEX,KKA,KKU,KKL,KNAERO,PTSTEP,& - & PPABST,PTHT,PDZZ,PRHODREF,PRCT,PRIT,& - & PCLDFR,PFPR, & - & PLSM,PUGST,PVGST, & - & PAERMR,PAERTEND) - - USE YOMHOOK , ONLY : LHOOK, DR_HOOK, JPHOOK -! ########################################################################## -! -!! Aerosol removal processes -!! -!! PURPOSE -!! ------- -!! Calculate the scavenging and deposition of aerosols. -!! -!! IMPLICIT ARGUMENTS -!! ------------------ -!! Module MODD_PARAMETERS -!! JPHEXT : Horizontal external points number -!! JPVEXT : Vertical external points number -!! Module MODD_CST -!! XPI -!! XG : Gravity constant -!! XP00 : Reference pressure -!! XRD : Gaz constant for dry air -!! XCPD : Cpd (dry air) -!! XRHOWL : density of liquid water -!! XRHOLI : density of ice water -!! Module MODD_RAIN_ICE_DESCR -!! XRTMIN : Min values allowed for mixing ratios -!! Module MODD_AEROSOL_PROP -!! XVDRY : Dry Deposition Velocity -!! XFRICS : Fracion of aerosols in aquaeous phase -!! : for in-cloud scavenging -!! -!! REFERENCE -!! --------- -!! Mocrette et al. 2009 and Redy, 2019. -!! -!! AUTHOR -!! ------ -!! Daniel Martin Perez (AEMET) -!! -!! MODIFICATIONS -!! ------------- -!! Original: 2018 -!! D. Martin : 27.03.2023 : Adition of sea salt emission. -!! D. Martin : 28.06.2022 : Modification of scavenging, addition of -!! snow fluxes. -!! -!------------------------------------------------------------------------- -!! -!* 0. DECLARATIONS -! ------------ - -USE MODD_PARAMETERS, ONLY : JPHEXT, JPVEXT -USE MODD_CST, ONLY : XG, XPI, XRHOLW, XRHOLI, XP00, XRD, XCPD -USE MODD_PHYEX, ONLY: PHYEX_t -USE MODD_AEROSOL_PROP, ONLY : XVDRY, XFRICS, NCOAR, XCOARSE, XSDAE, XRHOAE, XMRAE, & - & XERFUP, XERFDOWN, XKHYGROS, NDRYDEP, XDRYDEP, NSEASALT, XSEASALT, & - & NDUST, XDUST, XINTSS -USE MODD_NRT_AEROSOLS, ONLY : LAERDRDEP, LAERSSEM -! -IMPLICIT NONE -! -TYPE(PHYEX_t), INTENT(IN) :: PHYEX -INTEGER, INTENT(IN) :: KKA ! near ground array index -INTEGER, INTENT(IN) :: KKU ! uppest atmosphere array index -INTEGER, INTENT(IN) :: KKL ! vert. levels type 1=MNH -1=ARO -INTEGER, INTENT(IN) :: KNAERO ! number of aerosol species -REAL, INTENT(IN) :: PTSTEP ! Double Time step - ! (single if cold start) -REAL, DIMENSION(:,:,:), INTENT(IN) :: PPABST ! Pressure -REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHT ! Theta -REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ ! Layer thickness (m) -REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODREF ! Reference density -REAL, DIMENSION(:,:,:), INTENT(IN) :: PRCT ! Cloud water m.r. at t -REAL, DIMENSION(:,:,:), INTENT(IN) :: PRIT ! Pristine ice m.r. at t -REAL, DIMENSION(:,:,:), INTENT(IN) :: PCLDFR ! Convective Mass Flux Cloud fraction -REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PFPR ! upper-air precipitation fluxes -REAL, DIMENSION(:), INTENT(IN) :: PLSM ! Land/Sea Mask -REAL, DIMENSION(:), INTENT(IN) :: PUGST ! 10m u-wind gust -REAL, DIMENSION(:), INTENT(IN) :: PVGST ! 10m v-wind gust -REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PAERMR -REAL, DIMENSION(:,:,:,:), INTENT( OUT) :: PAERTEND ! Aerosol tendencies: emission + scavenging coefficient rate -! -!* 0.2 declaration of local variables -! -INTEGER :: IIB ! Define the domain where is -INTEGER :: IIE ! the microphysical sources have to be computed -INTEGER :: IJB ! -INTEGER :: IJE ! -INTEGER :: IKB ! -INTEGER :: IKE ! -INTEGER :: JC, IC -INTEGER :: JK, JI, JJ - -REAL, DIMENSION(SIZE(PRHODREF,1),SIZE(PRHODREF,2),SIZE(PRHODREF,3),KNAERO) :: & - ZW1, ZW2, ZW3, ZW4, & ! work arrays - ZAER_MR, & - ZAER_DEP ! Aerosol deposition -REAL, DIMENSION(SIZE(PRHODREF,1),SIZE(PRHODREF,2),SIZE(PRHODREF,3)) :: & - ZT, & ! Temperature - ZVELSED, & ! Velocity sedimentation - ZMFREEPATH, & ! mean free path - ZCUNN, & ! Cunningham factor for sedimentation - ZPFPRR, ZPFPRS - -REAL, DIMENSION(SIZE(PRHODREF,1),3) :: & - ZWS_SS ! work array - -! Sink and source at surface of aerosol in kg/m2 -REAL, DIMENSION(SIZE(PRHODREF,1)) :: & - ZAERSNKDU, & ! Desert dust sink - ZAERSNKSS, & ! Sea salt sink - ZAERSRCSS ! Sea salt source - -! XRTMIN = Minimum value for the mixing ratio -REAL :: ZINVTSTEP -REAL :: ZPFMIN ! Minimum Flux precipitation -REAL :: ZEF_RAIN, ZEF_SNOW -REAL :: ZRAY_RAIN, ZRAY_SNOW -REAL :: ZCLFK -REAL :: ZVDRY -REAL :: ZVISC -REAL :: ZAIRMAS -REAL :: ZBOLTZMANN -REAL :: ZAER_R -REAL :: Z10WIND - -REAL :: ZDKA, ZDKB, ZDKR, ZDKRMIN, ZFRICS, ZERFRAE - -!------------------------------------------------------------------------------- -! -!* 1. Parameters for aerosol wet deposition -REAL(KIND=JPHOOK) :: ZHOOK_HANDLE -IF (LHOOK) CALL DR_HOOK('AEROSOL_PROCESS',0,ZHOOK_HANDLE) - - -IIB=1+JPHEXT -IIE=SIZE(PDZZ,1) - JPHEXT -IJB=1+JPHEXT -IJE=SIZE(PDZZ,2) - JPHEXT -IKB=KKA+JPVEXT*KKL -IKE=KKU-JPVEXT*KKL -! -!* 1.2 COMPUTE SOME CONSTANT PARAMETERS -! -ZINVTSTEP=1./PTSTEP - -! -!* 2. compute the fluxes -! - -!ZPFMIN=1.0E-12 -ZPFMIN=0.0 -ZEF_RAIN=0.001 -ZEF_SNOW=0.01 -ZRAY_RAIN=1.0E-3 -ZRAY_SNOW=1.0E-3 -ZVISC=1.78E-5 ! Air viscosity (Pa*s) -ZAIRMAS=28.97E-3 ! Air mean molar mass(Kg*mol-1) -ZBOLTZMANN=8.31445 ! (J*K-1*mol-1) -ZDKRMIN=2.0E-8 - -ZAER_MR(:,:,:,:)=MAX(PAERMR(:,:,:,:),0.0) - -ZW1(:,:,:,:) = 0. -ZW2(:,:,:,:) = 0. -ZW3(:,:,:,:) = 0. -ZW4(:,:,:,:) = 0. - -ZWS_SS(:,:)=0.0 - -ZT(:,:,:) = PTHT(:,:,:) * ( PPABST(:,:,:) / XP00 ) ** (XRD/XCPD) - -ZAERSNKDU(:)=0.0 -ZAERSNKSS(:)=0.0 -ZAERSRCSS(:)=0.0 - -! ************************************************ -! AEROSOL REMOVAL PROCESSES -! ************************************************ - -! Rain and snow fluxes initialized with sedimentation fluxes - -DO JI = IIB,IIE - DO JJ = IJE,IJB - DO JK = IKE+1 , IKB, -1*KKL - ZPFPRR(JI,JJ,JK)=MAX(PFPR(JI,JJ,JK,3),0.0) - ZPFPRS(JI,JJ,JK)=MAX(PFPR(JI,JJ,JK,5),0.0) - ENDDO - ENDDO -ENDDO - - -! Aerosol wet deposition -! (Loop from top to bottom KKL=-1 ARO) -DO JK = IKE+1 , IKB, -1*KKL - DO JI = IIB,IIE - DO JJ = IJE,IJB - ! inside cloud - IF (ZPFPRR(JI,JJ,JK) > ZPFMIN .OR. & - & ZPFPRS(JI,JJ,JK) > ZPFMIN) THEN - ! inside cloud - ZCLFK=PCLDFR(JI,JJ,JK) - ZDKA=0.32538/ZT(JI,JJ,JK) - IF ((PRCT(JI,JJ,JK) > PHYEX%RAIN_ICE_DESCRN%XRTMIN(2) .OR. & - & PRIT(JI,JJ,JK) > PHYEX%RAIN_ICE_DESCRN%XRTMIN(4)) .AND. ZCLFK > 0.1) THEN - DO JC=1,KNAERO - ZW1(JI,JJ,JK,JC)=((ZPFPRR(JI,JJ,JK+KKL)-ZPFPRR(JI,JJ,JK))+& - & (ZPFPRS(JI,JJ,JK+KKL)-ZPFPRS(JI,JJ,JK)))/ & - & (PRCT(JI,JJ,JK)+PRIT(JI,JJ,JK))/ & - & (PRHODREF(JI,JJ,JK)*PDZZ(JI,JJ,JK)*ZCLFK) - ZDKB=XKHYGROS(JC) - IF ( ZDKB == 0.0 ) THEN - ZFRICS=XFRICS(JC) - ELSE - ZDKR=MAX(ZDKA/(6.75*ZDKB)**(0.33)/(0.05E-2)**(0.66),ZDKRMIN) - ZERFRAE=MAX(XERFDOWN(JC),0.5*erf(0.7071*LOG(ZDKR/XMRAE(JC))/LOG(XSDAE(JC)))) - ZFRICS=MIN((XERFUP(JC)-ZERFRAE)/(XERFUP(JC)-XERFDOWN(JC)),XFRICS(JC)) - ENDIF - ZW1(JI,JJ,JK,JC)=ZCLFK*ZFRICS*ZAER_MR(JI,JJ,JK,JC)*MIN(MAX(-1.0,ZW1(JI,JJ,JK,JC)),0.0) - ENDDO - ELSE - ! below cloud (rain) - DO JC=1,KNAERO - ZW2(JI,JJ,JK,JC)=MIN(-0.75*ZPFPRR(JI,JJ,JK)*ZEF_RAIN/& - & (ZRAY_RAIN*XRHOLW)- & - & 0.75*ZPFPRS(JI,JJ,JK)*ZEF_SNOW/(ZRAY_SNOW*XRHOLI),0.0)* & - & ZAER_MR(JI,JJ,JK,JC) - ENDDO - ENDIF - ENDIF - ENDDO - ENDDO -ENDDO - -! Aerosol dry deposition -! Surface layer -! -IF (LAERDRDEP) THEN - DO JI = IIB,IIE - DO JJ = IJE,IJB - IF (PRCT(JI,JJ,IKB) < PHYEX%RAIN_ICE_DESCRN%XRTMIN(2) .OR. & - & (ZPFPRR(JI,JJ,IKB) < ZPFMIN .AND. & - & ZPFPRS(JI,JJ,IKB) < ZPFMIN)) THEN - DO IC=1,NDRYDEP ! only for those aerosol with parametrized emission - JC=XDRYDEP(IC) - IF ( JC > KNAERO ) CYCLE - ZVDRY=(1.0-PLSM(JI))*XVDRY(JC,1)+PLSM(JI)*XVDRY(JC,2) - ZW3(JI,JJ,IKB,JC)=-0.01*ZVDRY*ZAER_MR(JI,JJ,IKB,JC)/& - & PDZZ(JI,JJ,IKB) - ENDDO - ENDIF - ENDDO - ENDDO -ENDIF - -ZAER_DEP(:,:,:,:)=ZW1(:,:,:,:)+ZW2(:,:,:,:)+ZW3(:,:,:,:) - - -! Dry sedimentation only for coarse modes (Redy, 2019) -ZMFREEPATH(:,:,:)=ZVISC/PPABST(:,:,:)*SQRT(0.5*XPI*ZBOLTZMANN*ZT(:,:,:)/ZAIRMAS) - -DO IC=1,NCOAR - JC=XCOARSE(IC) - IF ( JC > KNAERO ) CYCLE - ZAER_R=XMRAE(JC)*EXP(0.5*(LOG(XSDAE(JC))**2))*1.0E-6 - ZCUNN(:,:,:)=1.0+ZMFREEPATH(:,:,:)/ZAER_R*& - & (1.257+0.4*EXP(-1.1*ZAER_R/ZMFREEPATH(:,:,:))) - ZVELSED(:,:,:)=MIN(0.22*(ZAER_R**2)*XRHOAE(JC)*XG*ZCUNN(:,:,:)/ZVISC, & - & PDZZ(:,:,:)*ZINVTSTEP) - - DO JK = IKE+1 , IKB, -1*KKL - ZW4(:,:,JK,JC)=(ZAER_MR(:,:,JK+KKL,JC)*ZVELSED(:,:,JK+KKL)-& - & ZAER_MR(:,:,JK,JC)*ZVELSED(:,:,JK))/PDZZ(:,:,JK) - ENDDO - -ENDDO - -! Sum all the aerosol removal contributions - -PAERTEND(:,:,:,:)=MAX(ZAER_DEP(:,:,:,:)+ZW4(:,:,:,:),-0.8*ZAER_MR(:,:,:,:)*ZINVTSTEP) - -! Aerosol sinks on the surface (kg/m2/s) - -DO JI = IIB,IIE - ZAERSNKDU(JI)=0.0 - DO JJ = IJE,IJB - DO JK = IKE+1 , IKB, -1*KKL - DO IC=1,NSEASALT ! Sea Salt species - JC=XSEASALT(IC) - ! diagnostics of all aerosol deposition in kg/m2/s - ZAERSNKSS(JI)=ZAERSNKSS(JI)+(MIN(ZW1(JI,JJ,JK,JC),0.0)+ & - & MIN(ZW2(JI,JJ,JK,JC),0.0))*PRHODREF(JI,JJ,JK)*& - & PDZZ(JI,JJ,JK) - ENDDO - DO IC=1,NDUST ! Dust species - JC=XDUST(IC) - ! diagnostics of all aerosol deposition in kg/m2/s - ZAERSNKDU(JI)=ZAERSNKDU(JI)+(MIN(ZW1(JI,JJ,JK,JC),0.0)+ & - & MIN(ZW2(JI,JJ,JK,JC),0.0))*PRHODREF(JI,JJ,JK)*& - & PDZZ(JI,JJ,JK) - ENDDO - ENDDO - DO IC=1,NSEASALT ! Sea Salt species - JC=XSEASALT(IC) - ZAERSNKSS(JI)=ZAERSNKSS(JI)+(MIN(ZW3(JI,JJ,IKB,JC),0.0)+& - & MIN(ZW4(JI,JJ,IKB,JC),0.0))*& - & PRHODREF(JI,JJ,IKB)*PDZZ(JI,JJ,IKB) - ENDDO - DO IC=1,NDUST ! Dust species - JC=XDUST(IC) - ZAERSNKDU(JI)=ZAERSNKDU(JI)+(MIN(ZW3(JI,JJ,IKB,JC),0.0)+& - & MIN(ZW4(JI,JJ,IKB,JC),0.0))*& - & PRHODREF(JI,JJ,IKB)*PDZZ(JI,JJ,IKB) - ENDDO - ENDDO - ZAERSNKDU(JI)=MIN(ZAERSNKDU(JI),0.0) - ZAERSNKSS(JI)=MIN(ZAERSNKSS(JI),0.0) -ENDDO - - -!********************************* -! AEROSOL EMISSIONS -!********************************* -! -! Sea salt production (Monahan,1986) -! -IF (LAERSSEM) THEN - DO JI = IIB,IIE - ZAERSRCSS(JI)=0.0 - IF ( PLSM(JI)<0.5 ) THEN - Z10WIND=SQRT(PUGST(JI)**2+PVGST(JI)**2) - DO IC=1,NSEASALT ! Sea Salt species - JC=XSEASALT(IC) - ZWS_SS(JI,JC)=1.373*Z10WIND**(3.41)*XINTSS(JC)* & - & (XRHOAE(JC)/8.0)*(4.0*XPI/3.0)*1.0E-18 - ZAERSRCSS(JI)=ZAERSRCSS(JI)+ZWS_SS(JI,JC) - DO JJ = IJE,IJB - PAERTEND(JI,JJ,IKB,JC)=PAERTEND(JI,JJ,IKB,JC)+ZWS_SS(JI,JC)/& - & (PRHODREF(JI,JJ,IKB)*PDZZ(JI,JJ,IKB)) - ENDDO - ENDDO - ENDIF - ENDDO -ENDIF - - -IF (LHOOK) CALL DR_HOOK('AEROSOL_PROCESS',1,ZHOOK_HANDLE) -! -END SUBROUTINE AEROSOL_PROCESS diff --git a/micro/ini_aerosols_cams.F90 b/micro/ini_aerosols_cams.F90 deleted file mode 100644 index 8382b45d..00000000 --- a/micro/ini_aerosols_cams.F90 +++ /dev/null @@ -1,637 +0,0 @@ -! ######spl - SUBROUTINE INI_AEROSOLS_CAMS - USE YOMHOOK, ONLY : LHOOK, DR_HOOK, JPHOOK -! ########################################################### -! -!!**** *INI_AEROSOLS_CAMS* - initialize the constants necessary -!! for cams aerosols -!! -!! PURPOSE -!! ------- -!! The purpose of this routine is to initialize the constants used to -!! resolve the mixed phase microphysical scheme. -!! AEROSOLS FROM CAMS: A total number of 14 species. -!! -!!** METHOD -!! ------ -!! The constants are initialized to their numerical values and the number -!! -!! EXTERNAL -!! -------- -!! -!! -!! IMPLICIT ARGUMENTS -!! ------------------ -!! Module MODD_AEROSOL_PROP -!! -!! REFERENCE -!! --------- -!! -!! AUTHOR -!! ------ -!! -!! MODIFICATIONS -!! ------------- -!! Original D. Martin-Perez 11/02/2019 -!! -!------------------------------------------------------------------------------- -! -!* 0. DECLARATIONS -! ------------ -! -USE MODD_AEROSOL_PROP -! -IMPLICIT NONE -! -!* 0.1 Declarations of dummy arguments : -! -! -! -!* 0.2 Declarations of local variables : -! -! -LOGICAL :: GFLAG ! Logical flag for printing the constatnts on the output - ! listing -!------------------------------------------------------------------------------- -! -! -REAL(KIND=JPHOOK) :: ZHOOK_HANDLE -IF (LHOOK) CALL DR_HOOK('INI_AEROSOLS_CAMS',0,ZHOOK_HANDLE) -! -! -!------------------------------------------------------------------------------- -! -!* 1. CLOUD CONDENSATION NUMBER -! ------------------------- -! -! XACCNI: Index of aerosol species considered condensation nuclei - - -XACCNI(:)=0 -XIFNI(:)=0 -XCOARSE(:)=0 -XSEASALT(:)=0 -XDUST(:)=0 -XDRYDEP(:)=0 -XINTSS(:)=0 - -! Cloud condensation nuclei species -NCCN=9 ! Number of species becoming CCN -XACCNI(1)=1 ! Sea salt 1 -XACCNI(2)=2 ! Sea salt 2 -XACCNI(3)=3 ! Sea salt 3 -XACCNI(4)=8 ! Hydrophilic Organic Matter -XACCNI(5)=10 ! Hydrophilic Black Carbon -XACCNI(6)=11 ! Sulfate -XACCNI(7)=12 ! Nitrate fine mode -XACCNI(8)=13 ! Nitrate coarse mode -XACCNI(9)=14 ! Ammonium - -! Ice nuclei species -NIFN=5 ! Number of species becoming IFN -XIFNI(1)=4 ! Desert dust 1 -XIFNI(2)=5 ! Desert dust 2 -XIFNI(3)=6 ! Desert dust 3 -XIFNI(4)=7 ! Hydrophobic OM -XIFNI(5)=9 ! Hydrophobic BC - -! For supercoarse sea salt and dust (Remy 2019) -NCOAR=2 -XCOARSE(1)=3 -XCOARSE(2)=6 - -! Sea salt -NSEASALT=3 -XSEASALT(1)=1 -XSEASALT(2)=2 -XSEASALT(3)=3 -! Coarse Sea salt -NCOARSEASALT=2 -XCOARSEASALT(1)=2 -XCOARSEASALT(2)=3 - -! Desert Dust -NDUST=3 -XDUST(1)=4 -XDUST(2)=5 -XDUST(3)=6 - -! Dry Deposition activated -NDRYDEP=6 -XDRYDEP(1)=1 -XDRYDEP(2)=2 -XDRYDEP(3)=3 -XDRYDEP(4)=4 -XDRYDEP(5)=5 -XDRYDEP(6)=6 - - -! -!* 2. AEROSOL PHYSICAL PROPERTIES -! --------------------------- -! - -! -! Optical properties -! -! XMEXT: Mass Extinction (m2 kg-1) -! -! Aerosol distribution: log normal size distribution -! -! XMRAE: Modal Radius(micrometers),Rg -! XSDAE: geometric Standard Ddeviation,sigma -! -! Limits of the error function -! -! XBINDOWN: Lower bin limit (micrometers) -! XBINUP: Upper bin limit (micrometers) -! -! XERFDOWN(JC)=0.5*erf(log({DOWNlimit}/Rg)/(sqrt(2.0))*log(sigma) -! XERFUP(JC)=0.5*erf(log({UPlimit}/Rg)/(sqrt(2.0))*log(sigma) -! -! Mean cubic radius of the distribution: -! = Rg**3*((3.0/sqrt(2.0)))*log(sigma))**2 -! X = (3.0/sqrt(2.0)))*log(sigma) -! {ERF3DOWN,ERF3UP} = -! 0.5*erf(log({DOWNlimit,Uplimit}-X)/Rg)/(sqrt(2.0))*log(sigma) -! Correction of the mean cubic radius when the bin size limits are considered -! -! XR3CORR(JC)= *(ERF3UP-ERF3DOWN)/(XERFUP(JC)-XERFDOWN(JC) -! -! Physical properties -! -! XRHOAE: mass density (Kg m-3) -! XKHYGROS: Hygroscopic parameter -! -! XMEXT: Mass extinction -! -! Variables for removal processes -! -! XVDRY: Dry Deposition Velocity (cm s-1) (Remy et al., 2019) -! XFRICS: Fraction of aerosols in aqueous phase for in-cloud scavenging -! (Remy et al., 2019) - -XMEXT(:,:)=0.0 -XMRAE(:)=0.0 -XSDAE(:)=0.0 -XBINDOWN(:)=0.0 -XBINUP(:)=0.0 -XERFDOWN(:)=0.0 -XERFUP(:)=0.0 -XR3CORR(:)=0.0 -XRHOAE(:)=0.0 -XKHYGROS(:)=0.0 -XVDRY(:,:)=0.0 ! XVDRY= Dry deposition over (sea,land) -XFRICS(:)=0.0 - -!------------------------------------------------- -! 2.1 Sea salt -!------------------------------------------------- -! -! Sea salt aerosol density (Kg m-3) -! 2160 (dry particles) -! 1182 (sea salt production asuming 80% humidity, Mocrette et al. 2009) - -! Factor for Sea Salt emision -! Monahan, 1986 -! INTEGRAL (1+0.057*r**1.05)10**(1.19*exp(-B**2)) -! B=(0.380-log(r))/0.650 -! r: aerosol radius in um - -XINTSS(1)=2.09913 ! r: 2x0.03 - 2x0.5 um -XINTSS(2)=29.6649 ! r: 2x0.5 - 2x5.0 um -XINTSS(3)=80.3799 ! r: 2x5.0 - 2x20.0 um - -! -! Sea salt, size bin limits: 0.03 - 0.5 microm. -! Size distribution parameters: -XMRAE(1)=0.1992 -XSDAE(1)= 1.92 -XBINDOWN(1)=0.03 -XBINUP(1)=0.5 -XERFDOWN(1)=ERFPDF(XBINDOWN(1),XMRAE(1),XSDAE(1)) ! -0.4981 -XERFUP(1)=ERFPDF(XBINUP(1),XMRAE(1),XSDAE(1)) ! 0.4208 -! Mean cubic radius of the distribution: =5.3640e-02 -! Corrected mean cubic radius -XR3CORR(1)=1.7071e-02 -! Density (kg/m3): -XRHOAE(1)= 2160.0 - -XKHYGROS(1)=1.28 - -XVDRY(1,:)=(/0.05,0.2/) -XFRICS(1)=0.7 - -! Mass extinction -XMEXT(1,:)=(/827.98,755.72,684.17,614.08,2197.81,2612.42,3223.86,4048.34,5381.25,6506.83,8497.96,13330.20/) - -! -! Sea salt, size bin limits: 0.5 - 5.0 microm. -! - -! Size distribution parameters: -XMRAE(2)=1.9920 -XSDAE(2)= 2.00 -XBINDOWN(2)=0.5 -XBINUP(2)=5.0 -XERFDOWN(2)=ERFPDF(XBINDOWN(2),XMRAE(2),XSDAE(2)) ! -0.4769 -XERFUP(2)=ERFPDF(XBINUP(2),XMRAE(2),XSDAE(2)) ! 0.4079 -! Mean cubic radius of the distribution: =6.8680e+01 -! Corrected mean cubic radius -XR3CORR(2)=1.7549e+01 -! Density (kg/m3): -XRHOAE(2)=2160.0 - -XKHYGROS(2)=1.28 - -XVDRY(2,:)=(/0.05,0.5/) -XFRICS(2)=0.9 - -! Mass extinction -XMEXT(2,:)=(/143.79,143.36,143.30,143.72,285.42,330.74,376.77,432.98,520.54,590.13,713.12,1030.59/) - -! -! Sea salt, size bin limits: 5.0 - 20.0 microm. -! -! Size distribution parameters: -XMRAE(3)=1.9920 -XSDAE(3)=2.00 -XBINDOWN(3)=5.0 -XBINUP(3)=20.0 -XERFDOWN(3)=ERFPDF(XBINDOWN(3),XMRAE(3),XSDAE(3)) ! 0.4079 -XERFUP(3)=ERFPDF(XBINUP(3),XMRAE(3),XSDAE(3)) ! 0.4996 -! Mean cubic radius of the distribution: =6.8680e+01 -! Corrected mean cubic radius -XR3CORR(3)=5.0026e+02 -! Density (kg/m3): -XRHOAE(3)=2160.0 - -XKHYGROS(3)=1.28 - -XVDRY(3,:)=(/1.20,1.50/) -XFRICS(3)=0.9 - -! Mass extinction -XMEXT(3,:)=(/38.65,38.85,38.85,38.63,79.27,92.17,105.23,122.77,149.16,171.24,209.44,309.58/) - -!------------------------------------------------- -! 2.2 Desert Dust -!------------------------------------------------- - -! -! Desert Dust, size bin limits: 0.03 - 0.55 microm. -! -! Size distribution parameters: -XMRAE(4)=0.2900 -XSDAE(4)=2.00 -XBINDOWN(4)=0.03 -XBINUP(4)=0.55 -XERFDOWN(4)=ERFPDF(XBINDOWN(4),XMRAE(4),XSDAE(4)) ! -0.4995 -XERFUP(4)=ERFPDF(XBINUP(4),XMRAE(4),XSDAE(4)) ! 0.3221 -! Mean cubic radius of the distribution: =2.1191e-01 -! Corrected mean cubic radius -XR3CORR(4)=3.1940e-02 -! Density (kg/m3): -XRHOAE(4)=2610.0 - -XKHYGROS(4)=0.0 - -XVDRY(4,:)=(/0.02,0.02/) -XFRICS(4)=0.3 - -! Mass extinction -XMEXT(4,:)=2496.68 - -! -! Desert Dust, size bin limits: 0.55 - 0.9 microm. -! - -! Size distribution parameters: -XMRAE(5)=0.2900 -XSDAE(5)=2.00 -XBINDOWN(5)=0.55 -XBINUP(5)=0.9 -XERFDOWN(5)=ERFPDF(XBINDOWN(5),XMRAE(5),XSDAE(5)) ! 0.3221 -XERFUP(5)=ERFPDF(XBINUP(5),XMRAE(5),XSDAE(5)) ! 0.4489 -! Mean cubic radius of the distribution: =2.1191e-01 -! Corrected mean cubic radius -XR3CORR(5)=3.4124e-01 -! Density (kg/m3): -XRHOAE(5)=2610.0 - -XKHYGROS(5)=0.0 - -XVDRY(5,:)=(/0.2,0.2/) -XFRICS(5)=0.4 - -! Mass extinction -XMEXT(5,:)=955.078 - -! -! Desert Dust, size bin limits: 0.9 - 20.0 microm. -! - -! Size distribution parameters: -XMRAE(6)=0.2900 -XSDAE(6)=2.00 -XBINDOWN(6)=0.9 -XBINUP(6)=20.0 -XERFDOWN(6)=ERFPDF(XBINDOWN(6),XMRAE(6),XSDAE(6)) ! 0.4489 -XERFUP(6)=ERFPDF(XBINUP(6),XMRAE(6),XSDAE(6)) ! 0.5000 -! Mean cubic radius of the distribution: =2.1191e-01 -! Corrected mean cubic radius -XR3CORR(6)=2.7845e+00 -! Density (kg/m3): -XRHOAE(6)=2610.0 - -XKHYGROS(6)=0.0 - -XVDRY(6,:)=(/1.2,1.2/) -XFRICS(6)=0.4 - -! Mass extinction -XMEXT(6,:)=406.533 - -!------------------------------------------------- -! 2.3 Organic matter -!------------------------------------------------- -! - -! -! 7:Hydrophobic Organic Matter -! - -! Size distribution parameters: -XMRAE(7)=0.02 -XSDAE(7)=2.2 -XBINDOWN(7)=0.005 -XBINUP(7)=20.0 -XERFDOWN(7)=ERFPDF(XBINDOWN(7),XMRAE(7),XSDAE(7)) ! -0.4996 -XERFUP(7)=ERFPDF(XBINUP(7),XMRAE(7),XSDAE(7)) ! 0.4996 -! Mean cubic radius of the distribution: =1.0861e-03 -! Corrected mean cubic radius -XR3CORR(7)=9.7073e-04 -! Density (kg/m3): -XRHOAE(7)=1825.0 - -XKHYGROS(7)=0.0 - -XVDRY(7,:)=(/0.1,0.1/) -XFRICS(7)=0.0 - -! Mass extinction (RH=80) -XMEXT(7,:)=2321.03 - -! -! 8:Hydrophilic Organic Matter -! -! Size distribution parameters: -XMRAE(8)=0.02 -XSDAE(8)=2.2 -XBINDOWN(8)=0.005 -XBINUP(8)=20.0 -XERFDOWN(8)=ERFPDF(XBINDOWN(8),XMRAE(8),XSDAE(8)) ! -0.4996 -XERFUP(8)=ERFPDF(XBINUP(8),XMRAE(8),XSDAE(8)) ! 0.4996 -! Mean cubic radius of the distribution: =1.0861e-03 -! Corrected mean cubic radius -XR3CORR(8)=9.7073e-04 -! Density (kg/m3): -XRHOAE(8)=1825.0 - -XKHYGROS(8)=0.3 - -XVDRY(8,:)=(/0.1,0.1/) -XFRICS(8)=0.7 - -! Mass extinction -XMEXT(8,:)=(/1942.13,2131.58,2321.03,2510.48,2699.93,2889.38,3185.61,3481.84,4119.88,4909.06,5698.24,8206.75/) - -! -!------------------------------------------------- -! 2.4 Black Carbon -!------------------------------------------------- -! - -! Hydrophobic Black Carbon size bin limits 0.005 - 0.500 micrometers. - -! Size distribution (Oshima et al., 2009) -! XMRAE(9)=0.053; XSDAE(9)=1.5; XRHOAE(9)=1800. -! Size distribution parameters: -XMRAE(9)=0.0118 -XSDAE(9)=2.00 -XBINDOWN(9)=0.005 -XBINUP(9)=0.5 -XERFDOWN(9)=ERFPDF(XBINDOWN(9),XMRAE(9),XSDAE(9)) ! -0.3923 -XERFUP(9)=ERFPDF(XBINUP(9),XMRAE(9),XSDAE(9)) ! 0.5000 -! Mean cubic radius of the distribution: =1.4276e-05 -! Corrected mean cubic radius -XR3CORR(9)=1.5985e-05 -! Density (kg/m3): -XRHOAE(9)=1000.0 - -XKHYGROS(9)=0.0 - -XVDRY(9,:)=(/0.1,0.1/) -XFRICS(9)=0.0 - -XMEXT(9,:)=13487.8 - -! -! Hydrophilic Black Carbon size bin limits 0.005 - 0.500 micrometers. -! - - -! Size distribution parameters: -XMRAE(10)=0.0118 -XSDAE(10)=2.00 -XBINDOWN(10)=0.005 -XBINUP(10)=0.5 -XERFDOWN(10)=ERFPDF(XBINDOWN(10),XMRAE(10),XSDAE(10)) ! -0.3923 -XERFUP(10)=ERFPDF(XBINUP(10),XMRAE(10),XSDAE(10)) ! 0.5000 -! Mean cubic radius of the distribution: =1.4276e-05 -! Corrected mean cubic radius -XR3CORR(10)=1.5985e-05 -! Density (kg/m3): -XRHOAE(10)=1000.0 - -XKHYGROS(10)=0.1 - -XVDRY(10,:)=(/0.1,0.1/) -XFRICS(10)=0.7 - -! Mass extinction -XMEXT(10,:)=13487.8 - -! -!------------------------------------------------- -! 2.5. Sulphate properties (0.005-20 micrometers) -!------------------------------------------------- -! -! Sulfate (Vd=0.05 for oceans; Vd=0.25 all other surfaces ) - -! Size distribution parameters: -XMRAE(11)=0.0355 -XSDAE(11)=2.00 -XBINDOWN(11)=0.005 -XBINUP(11)=20.0 -XERFDOWN(11)=ERFPDF(XBINDOWN(11),XMRAE(11),XSDAE(11)) ! -0.4977 -XERFUP(11)=ERFPDF(XBINUP(11),XMRAE(11),XSDAE(11)) ! 0.5000 -! Mean cubic radius of the distribution: =3.8873e-04 -! Corrected mean cubic radius -XR3CORR(11)=3.8964e-04 -! Density (kg/m3): -XRHOAE(11)=1760.0 - -XKHYGROS(11)=0.6 - -XVDRY(11,:)=(/0.15,0.25/) -XFRICS(11)=0.7 - -! Mass extinction -XMEXT(11,:)=(/3119.62,3510.56,3901.51,4292.45,4683.40,5074.35,5685.64,6296.94,7613.58,9242.11,10870.64,16047.13/) - -! -!------------------------------------------------- -! 2.6. Nitrate properties (0.005-0.9 micrometers) -!------------------------------------------------- - -! -! Nitrate fine mode (NH4NO3), size bin limits: 0.005 - 0.9 microm. -! - -! Size distribution parameters: -XMRAE(12)=0.0355 -XSDAE(12)=2.0 -XBINDOWN(12)=0.005 -XBINUP(12)=0.9 -XERFDOWN(12)=ERFPDF(XBINDOWN(12),XMRAE(12),XSDAE(12)) ! -0.4977 -XERFUP(12)=ERFPDF(XBINUP(12),XMRAE(12),XSDAE(12)) ! 0.5000 -! Corrected mean cubic radius -XR3CORR(12)=3.8774e-04 -! Density (kg/m3): -XRHOAE(12)=1730.0 - -XKHYGROS(12)=0.64 - -XVDRY(12,:)=(/0.15,0.15/) -XFRICS(12)=0.4 - -XMEXT(12,:)=(/3501.201,3501.201,3501.201,3501.201,4757.47,5354.868,6161.467,7361.848,13555.21,15649.4,14113.68,23958.8/) - -! -! Nitrate coarse mode ( NaNO3+Ca(NO3)2 ), size bin limits: 0.9 - 20.0 microm. -! - -! Size distribution parameters: -XMRAE(13)=1.992 -XSDAE(13)=2.0 -XBINDOWN(13)=0.9 -XBINUP(13)=20.0 -XERFDOWN(13)=ERFPDF(XBINDOWN(13),XMRAE(13),XSDAE(13)) ! -0.3741 -XERFUP(13)=ERFPDF(XBINUP(13),XMRAE(13),XSDAE(13)) ! 0.4996 -! Corrected mean cubic radius -XR3CORR(13)=7.0228e+01 -! Density (kg/m3): -XRHOAE(13)=1400.0 - -XKHYGROS(13)=0.9 - -XVDRY(13,:)=(/0.15,0.15/) -XFRICS(13)=0.4 - -! Mass extinction -XMEXT(13,:)=(/4295.583,4295.583,4295.583,4295.583,5150.058,6185.985,6780.04,7425.708,8129.261,10575.19,14709.95,26338.47/) - -! -!------------------------------------------------- -! 2.7. Ammonium properties (0.005-20 micrometers) -!------------------------------------------------- -! (Ammonium sulfate for Amonia) -! Ammonium (Vd=0.05 for oceans; Vd=0.25 all other surfaces ) - -! Size distribution parameters: -XMRAE(14)=0.0355 -XSDAE(14)=2.00 -XBINDOWN(14)=0.005 -XBINUP(14)=20.0 -XERFDOWN(14)=ERFPDF(XBINDOWN(14),XMRAE(14),XSDAE(14)) ! -0.4977 -XERFUP(14)=ERFPDF(XBINUP(14),XMRAE(14),XSDAE(14)) ! 0.5000 -! Corrected mean cubic radius -XR3CORR(14)=3.8964e-04 -! Density (kg/m3): -XRHOAE(14)=1760.0 - -XKHYGROS(14)=0.6 - -XVDRY(14,:)=(/0.15,0.15/) -XFRICS(14)=0.4 - -! Mass extinction -XMEXT(14,:)=(/190.3313,232.0307,274.9777,321.6228,346.3498,371.914,425.8234,483.4848,544.9835,610.332,751.631,906.7416/) - - -! -IF (LHOOK) CALL DR_HOOK('INI_AEROSOL_CAMS',1,ZHOOK_HANDLE) - -CONTAINS -! -!---------------------------------------- -! - -FUNCTION ERFPDF(PRAD,PMODRAD,PSIGMA) RESULT (PERFPDF) -! -! 0.5*erf(log(PRAD/Rg)/(sqrt(2.0)*log(sigma))) -! -! - IMPLICIT NONE -! - REAL :: PRAD - REAL :: PMODRAD - REAL :: PSIGMA - REAL :: PERFPDF -! - REAL(KIND=JPHOOK) :: ZHOOK_HANDLE - IF (LHOOK) CALL DR_HOOK('INI_AEROSOLS_CAMS:ERFPDF',0,ZHOOK_HANDLE) -! - PERFPDF = 0.5*erf(0.7071*LOG(PRAD/PMODRAD)/LOG(PSIGMA)) -! - IF (LHOOK) CALL DR_HOOK('INI_AEROSOLS_CAMS:ERFPDF',1,ZHOOK_HANDLE) -! -END FUNCTION ERFPDF - - -FUNCTION RAD3CORR(PMODRAD,PSIGMA,PBINDOWN,PBINUP) RESULT (PR3CORR) -! Mean cubic radius of the distribution: -! = Rg**3*((3.0/sqrt(2.0)))*log(sigma))**2 -! X = (3.0/sqrt(2.0)))*log(sigma) -! {ERF3DOWN,ERF3UP} = -! 0.5*erf(log({DOWNlimit,Uplimit}-X)/Rg)/(sqrt(2.0))*log(sigma) -! Correction of the mean cubic radius when the bin size limits are considered -! -! XR3CORR(JC)= *(ERF3UP-ERF3DOWN)/(XERFUP(JC)-XERFDOWN(JC) -! - IMPLICIT NONE -! - REAL :: PMODRAD,PSIGMA,PBINDOWN,PBINUP - REAL :: PR3CORR - REAL :: ZR3, ZX,ZWUP,ZWDOWN,ZW3UP,ZW3DOWN -! - REAL(KIND=JPHOOK) :: ZHOOK_HANDLE - IF (LHOOK) CALL DR_HOOK('INI_AEROSOLS_CAMS:RAD3CORR',0,ZHOOK_HANDLE) - - ZX=2.12132*LOG(PSIGMA) - - ZR3=(PMODRAD**3)*ZX**2 - - ZWUP=ERFPDF(PBINUP,PMODRAD,PSIGMA) - ZWDOWN=ERFPDF(PBINDOWN,PMODRAD,PSIGMA) - ZW3UP=ERFPDF(PBINUP-ZX,PMODRAD,PSIGMA) - ZW3DOWN=ERFPDF(PBINDOWN-ZX,PMODRAD,PSIGMA) - - PR3CORR=ZR3*(ZW3UP-ZW3DOWN)/(ZWUP-ZWDOWN) - - IF (LHOOK) CALL DR_HOOK('INI_AEROSOLS_CAMS:RAD3CORR',1,ZHOOK_HANDLE) - -END FUNCTION RAD3CORR - -END SUBROUTINE INI_AEROSOLS_CAMS diff --git a/micro/mode_rain_ice_old_nucleation.F90 b/micro/mode_rain_ice_old_nucleation.F90 index ac9fae31..7c240de3 100644 --- a/micro/mode_rain_ice_old_nucleation.F90 +++ b/micro/mode_rain_ice_old_nucleation.F90 @@ -12,7 +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, PIFNNC, & + OAERONRT, OAEIFN, PIFNNC, & PTHS, PRVS, PRIS, PCIT, & PICENU, PT, PZZZ, & PRHT) @@ -22,7 +22,6 @@ SUBROUTINE RAIN_ICE_OLD_NUCLEATION(D, CST, ICEP, KSIZE, OCND2, LMODICEDEP, KRR, USE MODD_CST, ONLY: CST_T USE MODE_TIWMX, ONLY: ESATI, ESATW, AM3, REDIN USE MODD_RAIN_ICE_PARAM_N,ONLY: RAIN_ICE_PARAM_T - USE MODD_NRT_AEROSOLS, ONLY: LAEIFN ! !* 0. DECLARATIONS ! ------------ @@ -55,6 +54,7 @@ SUBROUTINE RAIN_ICE_OLD_NUCLEATION(D, CST, ICEP, KSIZE, OCND2, LMODICEDEP, KRR, 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 @@ -132,7 +132,7 @@ 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.LAEIFN) THEN + 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 diff --git a/micro/modi_rain_ice_old.F90 b/micro/modi_rain_ice_old.F90 index e28f162f..ad6e69e0 100644 --- a/micro/modi_rain_ice_old.F90 +++ b/micro/modi_rain_ice_old.F90 @@ -15,7 +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, PCLDROP, PIFNNC, & + OAERONRT, OAEIFN, PCLDROP, PIFNNC, & TBUDGETS, KBUDGETS, & PICENU, PKGN_ACON, PKGN_SBGR, & PRHT, PRHS, PINPRH, PFPR) @@ -101,6 +101,7 @@ SUBROUTINE RAIN_ICE_OLD (D, CST, PARAMI, ICEP, ICED, BUCONF, 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 ! diff --git a/micro/rain_ice_old.F90 b/micro/rain_ice_old.F90 index c6b9be1b..68a8a2ab 100644 --- a/micro/rain_ice_old.F90 +++ b/micro/rain_ice_old.F90 @@ -10,7 +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, PCLDROP, PIFNNC, & + OAERONRT, OAEIFN, PCLDROP, PIFNNC, & TBUDGETS, KBUDGETS, & PICENU, PKGN_ACON, PKGN_SBGR, & PRHT, PRHS, PINPRH, PFPR) @@ -183,7 +183,6 @@ SUBROUTINE RAIN_ICE_OLD (D, CST, PARAMI, ICEP, ICED, BUCONF, USE MODI_GAMMA, ONLY: GAMMA USE MODE_TIWMX, ONLY: ESATI, ESATW, AA2, BB3, AA2W, BB3W USE MODE_TIWMX_TAB, ONLY: TIWMX_TAB -USE MODD_NRT_AEROSOLS, ONLY : LAEIFN ! USE MODE_RAIN_ICE_OLD_NUCLEATION, ONLY: RAIN_ICE_OLD_NUCLEATION USE MODE_RAIN_ICE_OLD_SEDIMENTATION_STAT, ONLY: RAIN_ICE_OLD_SEDIMENTATION_STAT @@ -273,6 +272,7 @@ SUBROUTINE RAIN_ICE_OLD (D, CST, PARAMI, ICEP, ICED, BUCONF, 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 ! @@ -522,7 +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, PIFNNC, & + OAERONRT, OAEIFN, PIFNNC, & PTHS, PRVS, PRIS, PCIT, & PICENU, ZT, ZZZZ, & PRHT) From b57e83aca5ff0eb334338e3ed882b351e334b9a3 Mon Sep 17 00:00:00 2001 From: Alexandre MARY Date: Wed, 25 Sep 2024 13:30:45 +0000 Subject: [PATCH 3/3] bf turb single precision (commit c1fa35bb of IAL) --- turb/turb.F90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/turb/turb.F90 b/turb/turb.F90 index c0516fc3..5e6257ea 100644 --- a/turb/turb.F90 +++ b/turb/turb.F90 @@ -1545,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') !