diff --git a/src/ZSOLVER/ice_adjust.f90 b/src/ZSOLVER/ice_adjust.f90 new file mode 100644 index 0000000000000000000000000000000000000000..8da04e62e39528cbc163f6ecb66329684dcf72e6 --- /dev/null +++ b/src/ZSOLVER/ice_adjust.f90 @@ -0,0 +1,684 @@ +!MNH_LIC Copyright 1996-2021 CNRS, Meteo-France and Universite Paul Sabatier +!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence +!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt +!MNH_LIC for details. version 1. +!----------------------------------------------------------------- +! ###################### + MODULE MODI_ICE_ADJUST +! ###################### +! +INTERFACE +! + SUBROUTINE ICE_ADJUST (KKA, KKU, KKL, KRR, HFRAC_ICE, HCONDENS, HLAMBDA3,& + HBUNAME, OSUBG_COND, OSIGMAS, HSUBG_MF_PDF, & + PTSTEP, PSIGQSAT, & + PRHODJ, PEXNREF, PRHODREF, PSIGS, PMFCONV, & + PPABST, PZZ, & + PEXN, PCF_MF, PRC_MF, PRI_MF, & + PRV, PRC, PRVS, PRCS, PTH, PTHS, PSRCS, PCLDFR, & + PRR, PRI, PRIS, PRS, PRG, PRH, & + POUT_RV, POUT_RC, POUT_RI, POUT_TH, & + PHLC_HRC, PHLC_HCF, PHLI_HRI, PHLI_HCF ) +! +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) :: KRR ! Number of moist variables +CHARACTER(len=1), INTENT(IN) :: HFRAC_ICE +CHARACTER(len=80), INTENT(IN) :: HCONDENS +CHARACTER(len=4), INTENT(IN) :: HLAMBDA3 ! formulation for lambda3 coeff +CHARACTER(len=*), INTENT(IN) :: HBUNAME ! Name of the budget +LOGICAL, INTENT(IN) :: OSUBG_COND ! Switch for Subgrid + ! Condensation +LOGICAL :: OSIGMAS ! Switch for Sigma_s: + ! use values computed in CONDENSATION + ! or that from turbulence scheme +CHARACTER(len=*), INTENT(IN) :: HSUBG_MF_PDF +REAL, INTENT(IN) :: PTSTEP ! Double Time step + ! (single if cold start) +REAL, INTENT(IN) :: PSIGQSAT ! coeff applied to qsat variance contribution +! +REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODJ ! Dry density * Jacobian +REAL, DIMENSION(:,:,:), INTENT(IN) :: PEXNREF ! Reference Exner function +REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODREF +! +REAL, DIMENSION(:,:,:), INTENT(IN) :: PSIGS ! Sigma_s at time t +REAL, DIMENSION(:,:,:), INTENT(IN) :: PMFCONV ! convective mass flux +REAL, DIMENSION(:,:,:), INTENT(IN) :: PPABST ! Absolute Pressure at t +REAL, DIMENSION(:,:,:), INTENT(IN) :: PZZ ! height of model layer +REAL, DIMENSION(:,:,:), INTENT(IN) :: PEXN ! Exner function +! +REAL, DIMENSION(:,:,:), INTENT(IN) :: PCF_MF! Convective Mass Flux Cloud fraction +REAL, DIMENSION(:,:,:), INTENT(IN) :: PRI_MF! Convective Mass Flux ice mixing ratio +REAL, DIMENSION(:,:,:), INTENT(IN) :: PRC_MF! Convective Mass Flux liquid mixing ratio +! +REAL, DIMENSION(:,:,:), INTENT(IN) :: PRV ! Water vapor m.r. to adjust +REAL, DIMENSION(:,:,:), INTENT(IN) :: PRC ! Cloud water m.r. to adjust +REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PRVS ! Water vapor m.r. source +REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PRCS ! Cloud water m.r. source +REAL, DIMENSION(:,:,:), INTENT(IN) :: PTH ! Theta to adjust +REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PTHS ! Theta source +REAL, DIMENSION(:,:,:), INTENT(OUT) :: PSRCS ! Second-order flux + ! s'rc'/2Sigma_s2 at time t+1 + ! multiplied by Lambda_3 +REAL, DIMENSION(:,:,:), INTENT(OUT) :: PCLDFR ! Cloud fraction +REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PRIS ! Cloud ice m.r. at t+1 +REAL, DIMENSION(:,:,:), INTENT(IN) :: PRR ! Rain water m.r. to adjust +REAL, DIMENSION(:,:,:), INTENT(IN) :: PRI ! Cloud ice m.r. to adjust +REAL, DIMENSION(:,:,:), INTENT(IN) :: PRS ! Aggregate m.r. to adjust +REAL, DIMENSION(:,:,:), INTENT(IN) :: PRG ! Graupel m.r. to adjust +REAL, DIMENSION(:,:,:), OPTIONAL, INTENT(IN) :: PRH ! Hail m.r. to adjust +REAL, DIMENSION(:,:,:), OPTIONAL, INTENT(OUT) :: POUT_RV ! Adjusted value +REAL, DIMENSION(:,:,:), OPTIONAL, INTENT(OUT) :: POUT_RC ! Adjusted value +REAL, DIMENSION(:,:,:), OPTIONAL, INTENT(OUT) :: POUT_RI ! Adjusted value +REAL, DIMENSION(:,:,:), OPTIONAL, INTENT(OUT) :: POUT_TH ! Adjusted value +REAL, DIMENSION(:,:,:), OPTIONAL, INTENT(OUT) :: PHLC_HRC +REAL, DIMENSION(:,:,:), OPTIONAL, INTENT(OUT) :: PHLC_HCF +REAL, DIMENSION(:,:,:), OPTIONAL, INTENT(OUT) :: PHLI_HRI +REAL, DIMENSION(:,:,:), OPTIONAL, INTENT(OUT) :: PHLI_HCF +! +! +END SUBROUTINE ICE_ADJUST +! +END INTERFACE +! +END MODULE MODI_ICE_ADJUST + +! ########################################################################## + SUBROUTINE ICE_ADJUST (KKA, KKU, KKL, KRR, HFRAC_ICE, HCONDENS, HLAMBDA3,& + HBUNAME, OSUBG_COND, OSIGMAS, HSUBG_MF_PDF, & + PTSTEP, PSIGQSAT, & + PRHODJ, PEXNREF, PRHODREF, PSIGS, PMFCONV, & + PPABST, PZZ, & + PEXN, PCF_MF, PRC_MF, PRI_MF, & + PRV, PRC, PRVS, PRCS, PTH, PTHS, PSRCS, PCLDFR, & + PRR, PRI, PRIS, PRS, PRG, PRH, & + POUT_RV, POUT_RC, POUT_RI, POUT_TH, & + PHLC_HRC, PHLC_HCF, PHLI_HRI, PHLI_HCF ) +! ######################################################################### +! +!!**** *ICE_ADJUST* - compute the ajustment of water vapor in mixed-phase +!! clouds +!! +!! PURPOSE +!! ------- +!! The purpose of this routine is to compute the fast microphysical sources +!! through a saturation ajustement procedure in case of mixed-phase clouds. +!! +!! +!!** METHOD +!! ------ +!! Langlois, Tellus, 1973 for the cloudless version. +!! When cloud water is taken into account, refer to book 1 of the +!! documentation. +!! +!! +!! +!! EXTERNAL +!! -------- +!! None +!! +!! +!! IMPLICIT ARGUMENTS +!! ------------------ +!! Module MODD_CST +!! XP00 ! Reference pressure +!! XMD,XMV ! Molar mass of dry air and molar mass of vapor +!! XRD,XRV ! Gaz constant for dry air, gaz constant for vapor +!! XCPD,XCPV ! Cpd (dry air), Cpv (vapor) +!! XCL ! Cl (liquid) +!! XCI ! Ci (ice) +!! XTT ! Triple point temperature +!! XLVTT ! Vaporization heat constant +!! XLSTT ! Sublimation heat constant +!! XALPW,XBETAW,XGAMW ! Constants for saturation vapor over liquid +!! ! pressure function +!! XALPI,XBETAI,XGAMI ! Constants for saturation vapor over ice +!! ! pressure function +!! Module MODD_CONF +!! CCONF +!! Module MODD_BUDGET: +!! NBUMOD +!! CBUTYPE +!! LBU_RTH +!! LBU_RRV +!! LBU_RRC +!! LBU_RRI +!! +!! +!! REFERENCE +!! --------- +!! Book 1 and Book2 of documentation ( routine ICE_ADJUST ) +!! Langlois, Tellus, 1973 +!! +!! AUTHOR +!! ------ +!! J.-P. Pinty * Laboratoire d'Aerologie* +!! +!! +!! MODIFICATIONS +!! ------------- +!! Original 06/12/96 +!! M. Tomasini 27/11/00 Change CND and DEP fct of the T instead of rc and ri +!! Avoid the sub- and super-saturation before the ajustment +!! Avoid rc>0 if T<T00 before the ajustment +!! P Bechtold 12/02/02 change subgrid condensation +!! JP Pinty 29/11/02 add ICE2 and IC4 cases +!! (P. Jabouille) 27/05/04 safety test for case where esw/i(T)> pabs (~Z>40km) +!! J.Pergaud and S.Malardel Add EDKF case +!! S. Riette ice for EDKF +!! 2012-02 Y. Seity, add possibility to run with reversed vertical levels +!! J.Escobar : 15/09/2015 : WENO5 & JPHEXT <> 1 +!! 2016-07 S. Riette: adjustement is now realized on state variables (PRV, PRC, PRI, PTH) +!! whereas tendencies are still applied on S variables. +!! This modification allows to call ice_adjust on T variable +!! or to call it on S variables +!! 2016-11 S. Riette: all-or-nothing adjustment now uses condensation +! P. Wautelet 05/2016-04/2018: new data structures and calls for I/O +! P. Wautelet 02/2020: use the new data structures and subroutines for budgets +!------------------------------------------------------------------------------- +! +!* 0. DECLARATIONS +! ------------ +! +use modd_budget, only: lbudget_th, lbudget_rv, lbudget_rc, lbudget_ri, & + NBUDGET_TH, NBUDGET_RV, NBUDGET_RC, NBUDGET_RI, & + tbudgets +USE MODD_CONF +USE MODD_CST +USE MODD_PARAMETERS +USE MODD_RAIN_ICE_PARAM, ONLY : XCRIAUTC, XCRIAUTI, XACRIAUTI, XBCRIAUTI + +use mode_budget, only: Budget_store_init, Budget_store_end +#ifdef MNH_OPENACC +USE MODE_MNH_ZWORK, ONLY: MNH_MEM_GET, MNH_MEM_POSITION_PIN, MNH_MEM_RELEASE +#endif +USE MODE_MPPDB +#ifdef MNH_OPENACC +use mode_msg +#endif +use mode_tools_ll, only: GET_INDICE_ll + +#ifdef MNH_BITREP +USE MODI_BITREP +#endif +USE MODI_CONDENSATION +USE MODI_GET_HALO + +IMPLICIT NONE +! +! +!* 0.1 Declarations of dummy arguments : +! +! +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) :: KRR ! Number of moist variables +CHARACTER(len=1), INTENT(IN) :: HFRAC_ICE +CHARACTER(len=80), INTENT(IN) :: HCONDENS +CHARACTER(len=4), INTENT(IN) :: HLAMBDA3 ! formulation for lambda3 coeff +CHARACTER(len=*), INTENT(IN) :: HBUNAME ! Name of the budget +LOGICAL, INTENT(IN) :: OSUBG_COND ! Switch for Subgrid + ! Condensation +LOGICAL :: OSIGMAS ! Switch for Sigma_s: + ! use values computed in CONDENSATION + ! or that from turbulence scheme +CHARACTER(len=*), INTENT(IN) :: HSUBG_MF_PDF +REAL, INTENT(IN) :: PTSTEP ! Double Time step + ! (single if cold start) +REAL, INTENT(IN) :: PSIGQSAT ! coeff applied to qsat variance contribution +! +REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODJ ! Dry density * Jacobian +REAL, DIMENSION(:,:,:), INTENT(IN) :: PEXNREF ! Reference Exner function +REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODREF +! +REAL, DIMENSION(:,:,:), INTENT(IN) :: PSIGS ! Sigma_s at time t +REAL, DIMENSION(:,:,:), INTENT(IN) :: PMFCONV ! convective mass flux +REAL, DIMENSION(:,:,:), INTENT(IN) :: PPABST ! Absolute Pressure at t +REAL, DIMENSION(:,:,:), INTENT(IN) :: PZZ ! height of model layer +REAL, DIMENSION(:,:,:), INTENT(IN) :: PEXN ! Exner function +! +REAL, DIMENSION(:,:,:), INTENT(IN) :: PCF_MF! Convective Mass Flux Cloud fraction +REAL, DIMENSION(:,:,:), INTENT(IN) :: PRC_MF! Convective Mass Flux liquid mixing ratio +REAL, DIMENSION(:,:,:), INTENT(IN) :: PRI_MF! Convective Mass Flux ice mixing ratio +! +REAL, DIMENSION(:,:,:), INTENT(IN) :: PRV ! Water vapor m.r. to adjust +REAL, DIMENSION(:,:,:), INTENT(IN) :: PRC ! Cloud water m.r. to adjust +REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PRVS ! Water vapor m.r. source +REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PRCS ! Cloud water m.r. source +REAL, DIMENSION(:,:,:), INTENT(IN) :: PTH ! Theta to adjust +REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PTHS ! Theta source +REAL, DIMENSION(:,:,:), INTENT(OUT) :: PSRCS ! Second-order flux + ! s'rc'/2Sigma_s2 at time t+1 + ! multiplied by Lambda_3 +REAL, DIMENSION(:,:,:), INTENT(OUT) :: PCLDFR ! Cloud fraction +! +REAL, DIMENSION(:,:,:), INTENT(INOUT):: PRIS ! Cloud ice m.r. at t+1 +REAL, DIMENSION(:,:,:), INTENT(IN) :: PRR ! Rain water m.r. to adjust +REAL, DIMENSION(:,:,:), INTENT(IN) :: PRI ! Cloud ice m.r. to adjust +REAL, DIMENSION(:,:,:), INTENT(IN) :: PRS ! Aggregate m.r. to adjust +REAL, DIMENSION(:,:,:), INTENT(IN) :: PRG ! Graupel m.r. to adjust +REAL, DIMENSION(:,:,:), OPTIONAL, INTENT(IN) :: PRH ! Hail m.r. to adjust +REAL, DIMENSION(:,:,:), OPTIONAL, INTENT(OUT) :: POUT_RV ! Adjusted value +REAL, DIMENSION(:,:,:), OPTIONAL, INTENT(OUT) :: POUT_RC ! Adjusted value +REAL, DIMENSION(:,:,:), OPTIONAL, INTENT(OUT) :: POUT_RI ! Adjusted value +REAL, DIMENSION(:,:,:), OPTIONAL, INTENT(OUT) :: POUT_TH ! Adjusted value +REAL, DIMENSION(:,:,:), OPTIONAL, INTENT(OUT) :: PHLC_HRC +REAL, DIMENSION(:,:,:), OPTIONAL, INTENT(OUT) :: PHLC_HCF +REAL, DIMENSION(:,:,:), OPTIONAL, INTENT(OUT) :: PHLI_HRI +REAL, DIMENSION(:,:,:), OPTIONAL, INTENT(OUT) :: PHLI_HCF +! +!* 0.2 Declarations of local variables : +! +! +INTEGER :: IIU,IJU,IKU! dimensions of dummy arrays +INTEGER :: IIB,IJB ! Horz index values of the first inner mass points +INTEGER :: IIE,IJE ! Horz index values of the last inner mass points +INTEGER :: IKB ! K index value of the first inner mass point +INTEGER :: IKE ! K index value of the last inner mass point +INTEGER :: JITER,ITERMAX ! iterative loop for first order adjustment +INTEGER :: JI,JJ,JK +! +LOGICAL,DIMENSION(:,:,:), POINTER, CONTIGUOUS :: GTEMP +! +REAL, DIMENSION(:,:,:), POINTER, CONTIGUOUS :: ZSIGS,ZSRCS +REAL, DIMENSION(:,:,:), POINTER, CONTIGUOUS & + :: ZT, & ! adjusted temperature + ZRV, ZRC, ZRI, & ! adjusted state + ZCPH, & ! guess of the CPh for the mixing + ZLV, & ! guess of the Lv at t+1 + ZLS, & ! guess of the Ls at t+1 + ZW1,ZW2, & ! Work arrays for intermediate fields + ZCRIAUT, & ! Autoconversion thresholds + ZHCF, ZHR +! +!------------------------------------------------------------------------------- +! +! IN variables +! +!$acc data present( PRHODJ, PEXNREF, PSIGS, PMFCONV, PPABST, PZZ, & +!$acc & PEXN, PCF_MF, PRC_MF, PRI_MF, PRV, PRC, & +!$acc & PTH, PRR, PRI, PRS, PRG, PRH, & +! +! INOUT variables +! +!$acc & PRVS, PRCS, PRIS, PTHS, & +! +! OUT variables +! +!$acc & PSRCS, PCLDFR, POUT_RV, POUT_RC, POUT_RI, POUT_TH ) +! +! +!* 1. PRELIMINARIES +! ------------- +! +!Check all IN arrays +CALL MPPDB_CHECK3D(PRHODJ,"ICE_ADJUST beg:PRHODJ",PRECISION) +CALL MPPDB_CHECK3D(PEXNREF,"ICE_ADJUST beg:PEXNREF",PRECISION) +CALL MPPDB_CHECK3D(PRHODREF,"ICE_ADJUST beg:PRHODREF") +CALL MPPDB_CHECK3D(PSIGS,"ICE_ADJUST beg:PSIGS",PRECISION) +CALL MPPDB_CHECK3D(PMFCONV,"ICE_ADJUST beg:PMFCONV",PRECISION) +CALL MPPDB_CHECK3D(PPABST,"ICE_ADJUST beg:PPABST",PRECISION) +CALL MPPDB_CHECK3D(PZZ,"ICE_ADJUST beg:PZZ",PRECISION) +CALL MPPDB_CHECK3D(PEXN,"ICE_ADJUST beg:PEXN",PRECISION) +CALL MPPDB_CHECK3D(PCF_MF,"ICE_ADJUST beg:PCF_MF",PRECISION) +CALL MPPDB_CHECK3D(PRC_MF,"ICE_ADJUST beg:PRC_MF",PRECISION) +CALL MPPDB_CHECK3D(PRI_MF,"ICE_ADJUST begPRI_MF:",PRECISION) +CALL MPPDB_CHECK3D(PRV,"ICE_ADJUST beg:PRV",PRECISION) +CALL MPPDB_CHECK3D(PRC,"ICE_ADJUST beg:PRC",PRECISION) +CALL MPPDB_CHECK3D(PTH,"ICE_ADJUST beg:PTH",PRECISION) +CALL MPPDB_CHECK3D(PRR,"ICE_ADJUST beg:PRR",PRECISION) +CALL MPPDB_CHECK3D(PRI,"ICE_ADJUST beg:PRI",PRECISION) +CALL MPPDB_CHECK3D(PRS,"ICE_ADJUST beg:PRS",PRECISION) +CALL MPPDB_CHECK3D(PRG,"ICE_ADJUST beg:PRG",PRECISION) +IF (PRESENT(PRH)) CALL MPPDB_CHECK3D(PRH,"ICE_ADJUST beg:PRH",PRECISION) +!Check all INOUT arrays +CALL MPPDB_CHECK3D(PRVS,"ICE_ADJUST beg:PRVS",PRECISION) +CALL MPPDB_CHECK3D(PRCS,"ICE_ADJUST beg:PRCS",PRECISION) +CALL MPPDB_CHECK3D(PTHS,"ICE_ADJUST beg:PTHS",PRECISION) +CALL MPPDB_CHECK3D(PRIS,"ICE_ADJUST beg:PRIS",PRECISION) + +IIU = SIZE(PEXNREF,1) +IJU = SIZE(PEXNREF,2) +IKU = SIZE(PEXNREF,3) +CALL GET_INDICE_ll (IIB,IJB,IIE,IJE) +IKB=KKA+JPVEXT*KKL +IKE=KKU-JPVEXT*KKL + +#ifndef MNH_OPENACC +allocate( gtemp (IIU, IJU, IKU ) ) +allocate( zsigs (IIU, IJU, IKU ) ) +allocate( zsrcs (IIU, IJU, IKU ) ) +allocate( zt (IIU, IJU, IKU ) ) +allocate( zrv (IIU, IJU, IKU ) ) +allocate( zrc (IIU, IJU, IKU ) ) +allocate( zri (IIU, IJU, IKU ) ) +allocate( zcph (IIU, IJU, IKU ) ) +allocate( zlv (IIU, IJU, IKU ) ) +allocate( zls (IIU, IJU, IKU ) ) +allocate( zw1 (IIU, IJU, IKU ) ) +allocate( zw2 (IIU, IJU, IKU ) ) +allocate( zcriaut(IIU, IJU, IKU ) ) +allocate( zhcf (IIU, IJU, IKU ) ) +allocate( zhr (IIU, IJU, IKU ) ) +#else +!Pin positions in the pools of MNH memory +CALL MNH_MEM_POSITION_PIN() + +CALL MNH_MEM_GET( gtemp , IIU, IJU, IKU ) +CALL MNH_MEM_GET( zsigs , IIU, IJU, IKU ) +CALL MNH_MEM_GET( zsrcs , IIU, IJU, IKU ) +CALL MNH_MEM_GET( zt , IIU, IJU, IKU ) +CALL MNH_MEM_GET( zrv , IIU, IJU, IKU ) +CALL MNH_MEM_GET( zrc , IIU, IJU, IKU ) +CALL MNH_MEM_GET( zri , IIU, IJU, IKU ) +CALL MNH_MEM_GET( zcph , IIU, IJU, IKU ) +CALL MNH_MEM_GET( zlv , IIU, IJU, IKU ) +CALL MNH_MEM_GET( zls , IIU, IJU, IKU ) +CALL MNH_MEM_GET( zw1 , IIU, IJU, IKU ) +CALL MNH_MEM_GET( zw2 , IIU, IJU, IKU ) +CALL MNH_MEM_GET( zcriaut, IIU, IJU, IKU ) +CALL MNH_MEM_GET( zhcf , IIU, IJU, IKU ) +CALL MNH_MEM_GET( zhr , IIU, IJU, IKU ) + +!$acc data present( gtemp, zsigs, zsrcs, zt, zrv, zrc, zri, zcph, zlv, zls, zw1, zw2, zcriaut, zhcf, zhr ) +#endif + +if ( lbudget_th ) call Budget_store_init( tbudgets(NBUDGET_TH), trim( hbuname ), pths(:, :, :) * prhodj(:, :, :) ) +if ( lbudget_rv ) call Budget_store_init( tbudgets(NBUDGET_RV), trim( hbuname ), prvs(:, :, :) * prhodj(:, :, :) ) +if ( lbudget_rc ) call Budget_store_init( tbudgets(NBUDGET_RC), trim( hbuname ), prcs(:, :, :) * prhodj(:, :, :) ) +if ( lbudget_ri ) call Budget_store_init( tbudgets(NBUDGET_RI), trim( hbuname ), pris(:, :, :) * prhodj(:, :, :) ) +! +ITERMAX=1 +! +!------------------------------------------------------------------------------- +! +!* 2. COMPUTE QUANTITIES WITH THE GUESS OF THE FUTURE INSTANT +! ------------------------------------------------------- +! +! +! beginning of the iterative loop (to compute the adjusted state) +!$acc kernels +ZRV(:,:,:)=PRV(:,:,:) +ZRC(:,:,:)=PRC(:,:,:) +ZRI(:,:,:)=PRI(:,:,:) +ZT(:,:,:)=PTH(:,:,:) * PEXN(:,:,:) +!$acc end kernels +! +DO JITER =1,ITERMAX + ! + !* 2.3 compute the latent heat of vaporization Lv(T*) at t+1 + ! and the latent heat of sublimation Ls(T*) at t+1 + ! +!$acc kernels + ZLV(:,:,:) = XLVTT + ( XCPV - XCL ) * ( ZT(:,:,:) -XTT ) + ZLS(:,:,:) = XLSTT + ( XCPV - XCI ) * ( ZT(:,:,:) -XTT ) + ! + !* 2.4 compute the specific heat for moist air (Cph) at t+1 + ! + IF ( KRR == 7 ) THEN + ZCPH(:,:,:) = XCPD + XCPV * ZRV(:,:,:) & + + XCL * (ZRC(:,:,:) + PRR(:,:,:)) & + + XCI * (ZRI(:,:,:) + PRS(:,:,:) + PRG(:,:,:) + PRH(:,:,:)) + ELSE IF( KRR == 6 ) THEN + ZCPH(:,:,:) = XCPD + XCPV * ZRV(:,:,:) & + + XCL * (ZRC(:,:,:) + PRR(:,:,:)) & + + XCI * (ZRI(:,:,:) + PRS(:,:,:) + PRG(:,:,:)) + ELSE IF( KRR == 5 ) THEN + ZCPH(:,:,:) = XCPD + XCPV * ZRV(:,:,:) & + + XCL * (ZRC(:,:,:) + PRR(:,:,:)) & + + XCI * (ZRI(:,:,:) + PRS(:,:,:)) + ELSE IF( KRR == 3 ) THEN + ZCPH(:,:,:) = XCPD + XCPV * ZRV(:,:,:) & + + XCL * (ZRC(:,:,:) + PRR(:,:,:)) + ELSE IF( KRR == 2 ) THEN + ZCPH(:,:,:) = XCPD + XCPV * ZRV(:,:,:) & + + XCL * ZRC(:,:,:) + END IF +!$acc end kernels + ! + IF ( OSUBG_COND ) THEN + ! + !* 3. SUBGRID CONDENSATION SCHEME + ! --------------------------- + ! + ! PSRC= s'rci'/Sigma_s^2 + ! ZT, ZRV, ZRC and ZRI are INOUT + CALL CONDENSATION(IIU, IJU, IKU, IIB, IIE, IJB, IJE, IKB, IKE, KKL, & + HFRAC_ICE, HCONDENS, HLAMBDA3, & + PPABST, PZZ, PRHODREF, ZT, ZRV, ZRC, ZRI, PRS, PRG, PSIGS, PMFCONV, PCLDFR, & + PSRCS, .TRUE., OSIGMAS, & + PSIGQSAT, PLV=ZLV, PLS=ZLS, PCPH=ZCPH, PHLC_HRC=PHLC_HRC, PHLC_HCF=PHLC_HCF, PHLI_HRI=PHLI_HRI, PHLI_HCF=PHLI_HCF) + ELSE + ! + !* 4. ALL OR NOTHING CONDENSATION SCHEME + ! FOR MIXED-PHASE CLOUD + ! ----------------------------------------------- + ! + ! + ! ZT, ZRV, ZRC and ZRI are INOUT + ! + !CALL ADJUST_LANGLOIS(IIU, IJU, IKU, IIB, IIE, IJB, IJE, IKB, IKE, KKL, & + ! PPABST, ZT, ZRV, ZRC, ZRI, ZLV, ZLS, ZCPH) HFRAC_ICE must be implemented in Langlois before using it again +!PW: TODO: ZSRCS not used !!! (OUT variable from CONDENSATION) +!PW:TODO: autres optis: ZSIGS=0 and PSIGQSAT=0 -> less computation... +!$acc kernels + ZSIGS=0. +!$acc end kernels + CALL CONDENSATION(IIU, IJU, IKU, IIB, IIE, IJB, IJE, IKB, IKE, KKL, & + HFRAC_ICE, HCONDENS, HLAMBDA3, & + PPABST, PZZ, PRHODREF, ZT, ZRV, ZRC, ZRI, PRS, PRG, ZSIGS, PMFCONV, PCLDFR, & + ZSRCS, .TRUE., OSIGMAS=.TRUE., & + PSIGQSAT=0., PLV=ZLV, PLS=ZLS, PHLC_HRC=PHLC_HRC, PHLC_HCF=PHLC_HCF, PHLI_HRI=PHLI_HRI, PHLI_HCF=PHLI_HCF ) + END IF +ENDDO ! end of the iterative loop +! +!* 5. COMPUTE THE SOURCES AND STORES THE CLOUD FRACTION +! ------------------------------------------------- +! +!$acc kernels +#ifdef MNH_COMPILER_NVHPC +!$acc loop independent collapse(3) +#endif +DO CONCURRENT (JI=1:IIU,JJ=1:IJU,JK=1:IKU) + +! +!* 5.0 compute the variation of mixing ratio +! + ! Rc - Rc* +ZW1(JI,JJ,JK) = (ZRC(JI,JJ,JK) - PRC(JI,JJ,JK)) / PTSTEP ! Pcon = ---------- + ! 2 Delta t + +ZW2(JI,JJ,JK) = (ZRI(JI,JJ,JK) - PRI(JI,JJ,JK)) / PTSTEP ! idem ZW1 but for Ri +! +!* 5.1 compute the sources +! + +GTEMP(JI,JJ,JK) = ZW1(JI,JJ,JK) < 0.0 +IF( GTEMP(JI,JJ,JK) )THEN + ZW1(JI,JJ,JK) = MAX ( ZW1(JI,JJ,JK), -PRCS(JI,JJ,JK) ) +ELSE + ZW1(JI,JJ,JK) = MIN ( ZW1(JI,JJ,JK), PRVS(JI,JJ,JK) ) +ENDIF +PRVS(JI,JJ,JK) = PRVS(JI,JJ,JK) - ZW1(JI,JJ,JK) +PRCS(JI,JJ,JK) = PRCS(JI,JJ,JK) + ZW1(JI,JJ,JK) +PTHS(JI,JJ,JK) = PTHS(JI,JJ,JK) + & + ZW1(JI,JJ,JK) * ZLV(JI,JJ,JK) / (ZCPH(JI,JJ,JK) * PEXNREF(JI,JJ,JK)) +! +GTEMP(JI,JJ,JK) = ZW2(JI,JJ,JK) < 0.0 +IF( GTEMP(JI,JJ,JK) )THEN + ZW2(JI,JJ,JK) = MAX ( ZW2(JI,JJ,JK), -PRIS(JI,JJ,JK) ) +ELSE + ZW2(JI,JJ,JK) = MIN ( ZW2(JI,JJ,JK), PRVS(JI,JJ,JK) ) +ENDIF +PRVS(JI,JJ,JK) = PRVS(JI,JJ,JK) - ZW2(JI,JJ,JK) +PRIS(JI,JJ,JK) = PRIS(JI,JJ,JK) + ZW2(JI,JJ,JK) +PTHS(JI,JJ,JK) = PTHS(JI,JJ,JK) + & + ZW2(JI,JJ,JK) * ZLS(JI,JJ,JK) / (ZCPH(JI,JJ,JK) * PEXNREF(JI,JJ,JK)) +ENDDO + +!$acc end kernels +! +! +!* 5.2 compute the cloud fraction PCLDFR +! +IF ( .NOT. OSUBG_COND ) THEN +!$acc kernels + GTEMP(:,:,:) = PRCS(:,:,:) + PRIS(:,:,:) > 1.E-12 / PTSTEP + WHERE ( GTEMP(:,:,:) ) + PCLDFR(:,:,:) = 1. + ELSEWHERE + PCLDFR(:,:,:) = 0. + ENDWHERE + IF ( SIZE(PSRCS,3) /= 0 ) THEN + PSRCS(:,:,:) = PCLDFR(:,:,:) + END IF +!$acc end kernels +ELSE +#ifdef MNH_OPENACC + call Print_msg( NVERB_FATAL, 'GEN', 'ICE_ADJUST', 'OpenACC: OSUBG_COND=.T. not yet implemented' ) +#endif + !We limit PRC_MF+PRI_MF to PRVS*PTSTEP to avoid negative humidity + ZW1(:,:,:)=PRC_MF(:,:,:)/PTSTEP + ZW2(:,:,:)=PRI_MF(:,:,:)/PTSTEP + GTEMP(:,:,:) = ZW1(:,:,:)+ZW2(:,:,:)>PRVS(:,:,:) + WHERE(GTEMP(:,:,:)) + ZW1(:,:,:)=ZW1(:,:,:)*PRVS(:,:,:)/(ZW1(:,:,:)+ZW2(:,:,:)) + ZW2(:,:,:)=PRVS(:,:,:)-ZW1(:,:,:) + ENDWHERE + IF(PRESENT(PHLC_HRC) .AND. PRESENT(PHLC_HCF)) THEN + ZCRIAUT(:,:,:)=XCRIAUTC/PRHODREF + IF(HSUBG_MF_PDF=='NONE')THEN + WHERE(ZW1(:,:,:)*PTSTEP>PCF_MF * ZCRIAUT) + PHLC_HRC(:,:,:)=PHLC_HRC(:,:,:)+ZW1(:,:,:)*PTSTEP + PHLC_HCF(:,:,:)=MIN(1.,PHLC_HCF(:,:,:)+PCF_MF(:,:,:)) + ENDWHERE + ELSEIF(HSUBG_MF_PDF=='TRIANGLE')THEN + !ZHCF is the precipitating part of the *cloud* and not of the grid cell + WHERE(ZW1(:,:,:)*PTSTEP>PCF_MF*ZCRIAUT(:,:,:)) + ZHCF(:,:,:)=1.-.5*(ZCRIAUT(:,:,:)*PCF_MF(:,:,:) / MAX(1.E-20, ZW1(:,:,:)*PTSTEP))**2 + ZHR(:,:,:)=ZW1(:,:,:)*PTSTEP-(ZCRIAUT(:,:,:)*PCF_MF(:,:,:))**3 / & + &(3*MAX(1.E-20, ZW1(:,:,:)*PTSTEP)**2) + ELSEWHERE(2.*ZW1(:,:,:)*PTSTEP<=PCF_MF * ZCRIAUT(:,:,:)) + ZHCF(:,:,:)=0. + ZHR(:,:,:)=0. + ELSEWHERE + ZHCF(:,:,:)=(2.*ZW1(:,:,:)*PTSTEP-ZCRIAUT(:,:,:)*PCF_MF(:,:,:))**2 / & + &(2.*MAX(1.E-20, ZW1(:,:,:)*PTSTEP)**2) + ZHR(:,:,:)=(4.*(ZW1(:,:,:)*PTSTEP)**3-3.*ZW1(:,:,:)*PTSTEP*(ZCRIAUT(:,:,:)*PCF_MF(:,:,:))**2+& + (ZCRIAUT(:,:,:)*PCF_MF(:,:,:))**3) / & + &(3*MAX(1.E-20, ZW1(:,:,:)*PTSTEP)**2) + ENDWHERE + ZHCF(:,:,:)=ZHCF(:,:,:)*PCF_MF(:,:,:) !to retrieve the part of the grid cell + PHLC_HCF(:,:,:)=MIN(1.,PHLC_HCF(:,:,:)+ZHCF(:,:,:)) !total part of the grid cell that is precipitating + PHLC_HRC(:,:,:)=PHLC_HRC(:,:,:)+ZHR(:,:,:) + ENDIF + ENDIF + IF(PRESENT(PHLI_HRI) .AND. PRESENT(PHLI_HCF)) THEN + ZCRIAUT(:,:,:)=MIN(XCRIAUTI,10**(XACRIAUTI*(ZT(:,:,:)-XTT)+XBCRIAUTI)) + IF(HSUBG_MF_PDF=='NONE')THEN + WHERE(ZW2(:,:,:)*PTSTEP>PCF_MF * ZCRIAUT(:,:,:)) + PHLI_HRI(:,:,:)=PHLI_HRI(:,:,:)+ZW2(:,:,:)*PTSTEP + PHLI_HCF(:,:,:)=MIN(1.,PHLI_HCF(:,:,:)+PCF_MF(:,:,:)) + ENDWHERE + ELSEIF(HSUBG_MF_PDF=='TRIANGLE')THEN + !ZHCF is the precipitating part of the *cloud* and not of the grid cell + WHERE(ZW2(:,:,:)*PTSTEP>PCF_MF*ZCRIAUT) + ZHCF(:,:,:)=1.-.5*(ZCRIAUT*PCF_MF(:,:,:) / (ZW2(:,:,:)*PTSTEP))**2 + ZHR(:,:,:)=ZW2(:,:,:)*PTSTEP-(ZCRIAUT*PCF_MF(:,:,:))**3/(3*(ZW2(:,:,:)*PTSTEP)**2) + ELSEWHERE(2.*ZW2(:,:,:)*PTSTEP<=PCF_MF * ZCRIAUT) + ZHCF(:,:,:)=0. + ZHR(:,:,:)=0. + ELSEWHERE + ZHCF(:,:,:)=(2.*ZW2(:,:,:)*PTSTEP-ZCRIAUT*PCF_MF(:,:,:))**2 / (2.*(ZW2(:,:,:)*PTSTEP)**2) + ZHR(:,:,:)=(4.*(ZW2(:,:,:)*PTSTEP)**3-3.*ZW2(:,:,:)*PTSTEP*(ZCRIAUT*PCF_MF(:,:,:))**2+& + (ZCRIAUT*PCF_MF(:,:,:))**3)/(3*(ZW2(:,:,:)*PTSTEP)**2) + ENDWHERE + ZHCF(:,:,:)=ZHCF(:,:,:)*PCF_MF(:,:,:) !to retrieve the part of the grid cell + PHLI_HCF(:,:,:)=MIN(1.,PHLI_HCF(:,:,:)+ZHCF(:,:,:)) !total part of the grid cell that is precipitating + PHLI_HRI(:,:,:)=PHLI_HRI(:,:,:)+ZHR(:,:,:) + ENDIF + ENDIF + PCLDFR(:,:,:)=MIN(1.,PCLDFR(:,:,:)+PCF_MF(:,:,:)) + PRCS(:,:,:)=PRCS(:,:,:)+ZW1(:,:,:) + PRIS(:,:,:)=PRIS(:,:,:)+ZW2(:,:,:) + PRVS(:,:,:)=PRVS(:,:,:)-(ZW1(:,:,:)+ZW2(:,:,:)) + PTHS(:,:,:) = PTHS(:,:,:) + & + (ZW1 * ZLV(:,:,:) + ZW2 * ZLS(:,:,:)) / ZCPH(:,:,:) & + / PEXNREF(:,:,:) + IF(PRESENT(POUT_RV) .OR. PRESENT(POUT_RC) .OR. & + &PRESENT(POUT_RI) .OR. PRESENT(POUT_TH)) THEN + ZW1(:,:,:)=PRC_MF(:,:,:) + ZW2(:,:,:)=PRI_MF(:,:,:) + GTEMP(:,:,:) = ZW1(:,:,:)+ZW2(:,:,:)>ZRV(:,:,:) + WHERE(GTEMP(:,:,:)) + ZW1(:,:,:)=ZW1(:,:,:)*ZRV(:,:,:)/(ZW1(:,:,:)+ZW2(:,:,:)) + ZW2(:,:,:)=ZRV(:,:,:)-ZW1(:,:,:) + ENDWHERE + ZRC(:,:,:)=ZRC(:,:,:)+ZW1(:,:,:) + ZRI(:,:,:)=ZRI(:,:,:)+ZW2(:,:,:) + ZRV(:,:,:)=ZRV(:,:,:)-(ZW1(:,:,:)+ZW2(:,:,:)) + ZT(:,:,:) = ZT(:,:,:) + & + (ZW1 * ZLV(:,:,:) + ZW2 * ZLS(:,:,:)) / ZCPH(:,:,:) + ENDIF +ENDIF +! +!$acc kernels +IF(PRESENT(POUT_RV)) POUT_RV=ZRV +IF(PRESENT(POUT_RC)) POUT_RC=ZRC +IF(PRESENT(POUT_RI)) POUT_RI=ZRI +IF(PRESENT(POUT_TH)) POUT_TH=ZT / PEXN(:,:,:) +!$acc end kernels +! +! +!* 6. STORE THE BUDGET TERMS +! ---------------------- +! +if ( lbudget_th ) then +!$acc update self(pths) + call Budget_store_end( tbudgets(NBUDGET_TH), trim( hbuname ), pths(:, :, :) * prhodj(:, :, :) ) +end if +if ( lbudget_rv ) then +!$acc update self(prvs) + call Budget_store_end( tbudgets(NBUDGET_RV), trim( hbuname ), prvs(:, :, :) * prhodj(:, :, :) ) +end if +if ( lbudget_rc ) then +!$acc update self(prcs) + call Budget_store_end( tbudgets(NBUDGET_RC), trim( hbuname ), prcs(:, :, :) * prhodj(:, :, :) ) +end if +if ( lbudget_ri ) then +!$acc update self(pris) + call Budget_store_end( tbudgets(NBUDGET_RI), trim( hbuname ), pris(:, :, :) * prhodj(:, :, :) ) +end if + +!$acc end data + +#ifndef MNH_OPENACC +deallocate( gtemp, zsigs, zsrcs, zt, zrv, zrc, zri, zcph, zlv, zls, zw1, zw2, zcriaut, zhcf, zhr ) +#else +!Release all memory allocated with MNH_MEM_GET calls since last call to MNH_MEM_POSITION_PIN +CALL MNH_MEM_RELEASE() +#endif +!------------------------------------------------------------------------------ +! +!Check all INOUT arrays +CALL MPPDB_CHECK3D(PRVS,"ICE_ADJUST end:PRVS",PRECISION) +CALL MPPDB_CHECK3D(PRCS,"ICE_ADJUST end:PRCS",PRECISION) +CALL MPPDB_CHECK3D(PTHS,"ICE_ADJUST end:PTHS",PRECISION) +CALL MPPDB_CHECK3D(PRIS,"ICE_ADJUST end:PRIS",PRECISION) +!Check all OUT arrays +CALL MPPDB_CHECK3D(PSRCS,"ICE_ADJUST end:PSRCS",PRECISION) +CALL MPPDB_CHECK3D(PCLDFR,"ICE_ADJUST end:PCLDFR",PRECISION) +IF (PRESENT(POUT_RV)) CALL MPPDB_CHECK3D(POUT_RV,"ICE_ADJUST end:POUT_RV") +IF (PRESENT(POUT_RC)) CALL MPPDB_CHECK3D(POUT_RC,"ICE_ADJUST end:POUT_RC") +IF (PRESENT(POUT_RI)) CALL MPPDB_CHECK3D(POUT_RI,"ICE_ADJUST end:POUT_RI") +IF (PRESENT(POUT_TH)) CALL MPPDB_CHECK3D(POUT_TH,"ICE_ADJUST end:POUT_TH") +IF (PRESENT(PHLC_HRC)) CALL MPPDB_CHECK3D(PHLC_HRC,"ICE_ADJUST end:PHLC_HRC") +IF (PRESENT(PHLC_HCF)) CALL MPPDB_CHECK3D(PHLC_HCF,"ICE_ADJUST end:PHLC_HCF") +IF (PRESENT(PHLI_HRI)) CALL MPPDB_CHECK3D(PHLI_HRI,"ICE_ADJUST end:PHLI_HRI") +IF (PRESENT(PHLI_HCF)) CALL MPPDB_CHECK3D(PHLI_HCF,"ICE_ADJUST end:PHLI_HCF") + +!$acc end data + +END SUBROUTINE ICE_ADJUST