diff --git a/src/ZSOLVER/advection_uvw.f90 b/src/ZSOLVER/advection_uvw.f90 index 9660d8993df2f27f18cb2bc2061011f3d30e7aa1..466ff36cee4a99bcbe92d014d8297557b3a18717 100644 --- a/src/ZSOLVER/advection_uvw.f90 +++ b/src/ZSOLVER/advection_uvw.f90 @@ -197,6 +197,8 @@ TYPE(LIST_ll), POINTER :: TZFIELDS0_ll ! list of fields to exchange ! INTEGER :: IIU,IJU,IKU ! +INTEGER :: JI,JJ,JK +! !------------------------------------------------------------------------------- ! !* 0. INITIALIZATION @@ -273,8 +275,8 @@ CALL MNH_MEM_GET( ZMZM_RHODJ, IIU, IJU, IKU ) #endif !$acc data present( zrut, zrvt, zrwt, zruct, zrvct, zrwct, zu, zv, zw, & -!$acc & zrus_other, zrvs_other, zrws_other, zrus_adv, zrvs_adv, zrws_adv, & -!$acc & zmxm_rhodj, zmym_rhodj, zmzm_rhodj ) +!$acc & zrus_other, zrvs_other, zrws_other, zrus_adv, zrvs_adv, zrws_adv, & +!$acc & zmxm_rhodj, zmym_rhodj, zmzm_rhodj ) IKE = SIZE(PWT,3) - JPVEXT ! @@ -435,24 +437,32 @@ DO JSPL=1,ISPLIT ! Tendencies on wind !$acc update device(ZRUS_ADV,ZRVS_ADV,ZRWS_ADV) !$acc kernels - PRUS(:,:,:) = PRUS(:,:,:) + ZRUS_ADV(:,:,:) / ISPLIT - PRVS(:,:,:) = PRVS(:,:,:) + ZRVS_ADV(:,:,:) / ISPLIT - PRWS(:,:,:) = PRWS(:,:,:) + ZRWS_ADV(:,:,:) / ISPLIT -!$acc end kernels -! +#ifdef MNH_COMPILER_NVHPC +!$acc loop independent collapse(3) +#endif +DO CONCURRENT (JI=1:IIU , JJ=1:IJU , JK=1:IKU ) + PRUS(JI,JJ,JK) = PRUS(JI,JJ,JK) + ZRUS_ADV(JI,JJ,JK) / ISPLIT + PRVS(JI,JJ,JK) = PRVS(JI,JJ,JK) + ZRVS_ADV(JI,JJ,JK) / ISPLIT + PRWS(JI,JJ,JK) = PRWS(JI,JJ,JK) + ZRWS_ADV(JI,JJ,JK) / ISPLIT +END DO IF (JSPL<ISPLIT) THEN ! ! Guesses for next time splitting loop - ! - !$acc kernels - ZU(:,:,:) = ZU(:,:,:) + ZTSTEP / ZMXM_RHODJ * & - (ZRUS_OTHER(:,:,:) + ZRUS_ADV(:,:,:)) - ZV(:,:,:) = ZV(:,:,:) + ZTSTEP / ZMYM_RHODJ * & - (ZRVS_OTHER(:,:,:) + ZRVS_ADV(:,:,:)) - ZW(:,:,:) = ZW(:,:,:) + ZTSTEP / ZMZM_RHODJ * & - (ZRWS_OTHER(:,:,:) + ZRWS_ADV(:,:,:)) - !$acc end kernels - END IF +! +! +#ifdef MNH_COMPILER_NVHPC +!$acc loop independent collapse(3) +#endif +DO CONCURRENT (JI=1:IIU , JJ=1:IJU , JK=1:IKU ) + ZU(JI,JJ,JK) = ZU(JI,JJ,JK) + ZTSTEP / ZMXM_RHODJ(JI,JJ,JK) * & + (ZRUS_OTHER(JI,JJ,JK) + ZRUS_ADV(JI,JJ,JK)) + ZV(JI,JJ,JK) = ZV(JI,JJ,JK) + ZTSTEP / ZMYM_RHODJ(JI,JJ,JK) * & + (ZRVS_OTHER(JI,JJ,JK) + ZRVS_ADV(JI,JJ,JK)) + ZW(JI,JJ,JK) = ZW(JI,JJ,JK) + ZTSTEP / ZMZM_RHODJ(JI,JJ,JK) * & + (ZRWS_OTHER(JI,JJ,JK) + ZRWS_ADV(JI,JJ,JK)) +END DO +END IF +!$acc end kernels ! ! Top and bottom Boundaries ! diff --git a/src/ZSOLVER/tke_eps_sources.f90 b/src/ZSOLVER/tke_eps_sources.f90 new file mode 100644 index 0000000000000000000000000000000000000000..15c82a0df4d52eef8274009e142bf6ee840d993d --- /dev/null +++ b/src/ZSOLVER/tke_eps_sources.f90 @@ -0,0 +1,682 @@ +!MNH_LIC Copyright 1994-2022 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_TKE_EPS_SOURCES +! ########################### +INTERFACE +! + SUBROUTINE TKE_EPS_SOURCES(KKA,KKU,KKL,KMI,PTKEM,PLM,PLEPS,PDP,PTRH, & + PRHODJ,PDZZ,PDXX,PDYY,PDZX,PDZY,PZZ, & + PTSTEP,PIMPL,PEXPL, & + HTURBLEN,HTURBDIM, & + TPFILE,OTURB_DIAG, & + PTP,PRTKES,PRTKESM, PRTHLS,PCOEF_DISS,PTR,PDISS ) +! +USE MODD_IO, ONLY: TFILEDATA +! +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) :: KMI ! model index number +REAL, DIMENSION(:,:,:), INTENT(IN) :: PTKEM ! TKE at t-deltat +REAL, DIMENSION(:,:,:), INTENT(IN) :: PLM ! mixing length +REAL, DIMENSION(:,:,:), INTENT(IN) :: PLEPS ! dissipative length +REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODJ ! density * grid volume +REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX,PDYY,PDZZ,PDZX,PDZY + ! metric coefficients +REAL, DIMENSION(:,:,:), INTENT(IN) :: PZZ ! physical height w-pt +REAL, INTENT(IN) :: PTSTEP ! Time step +REAL, INTENT(IN) :: PEXPL, PIMPL ! Coef. temporal. disc. +CHARACTER(len=4), INTENT(IN) :: HTURBDIM ! dimensionality of the + ! turbulence scheme +CHARACTER(len=4), INTENT(IN) :: HTURBLEN ! kind of mixing length +TYPE(TFILEDATA), INTENT(INOUT):: TPFILE ! Output file +LOGICAL, INTENT(IN) :: OTURB_DIAG ! switch to write some + ! diagnostic fields in the syncronous FM-file +REAL, DIMENSION(:,:,:), INTENT(INOUT):: PDP ! Dyn. prod. of TKE +REAL, DIMENSION(:,:,:), INTENT(IN) :: PTRH +REAL, DIMENSION(:,:,:), INTENT(IN) :: PTP ! Ther. prod. of TKE +REAL, DIMENSION(:,:,:), INTENT(INOUT):: PRTKES ! RHOD * Jacobian * + ! TKE at t+deltat +REAL, DIMENSION(:,:,:), INTENT(IN) :: PRTKESM ! Advection source +REAL, DIMENSION(:,:,:), INTENT(INOUT):: PRTHLS ! Source of Theta_l +REAL, DIMENSION(:,:,:), INTENT(IN) :: PCOEF_DISS ! 1/(Cph*Exner) +REAL, DIMENSION(:,:,:), INTENT(OUT) :: PTR ! Transport prod. of TKE +REAL, DIMENSION(:,:,:), INTENT(OUT) :: PDISS ! Dissipati prod. of TKE +! +! +! +END SUBROUTINE TKE_EPS_SOURCES +! +END INTERFACE +! +END MODULE MODI_TKE_EPS_SOURCES +! +! ################################################################## + SUBROUTINE TKE_EPS_SOURCES(KKA,KKU,KKL,KMI,PTKEM,PLM,PLEPS,PDP, & + PTRH,PRHODJ,PDZZ,PDXX,PDYY,PDZX,PDZY,PZZ, & + PTSTEP,PIMPL,PEXPL, & + HTURBLEN,HTURBDIM, & + TPFILE,OTURB_DIAG, & + PTP,PRTKES,PRTKESM, PRTHLS,PCOEF_DISS,PTR,PDISS ) +! ################################################################## +! +! +!!**** *TKE_EPS_SOURCES* - routine to compute the sources of the turbulent +!! evolutive variables: TKE and its dissipation when it is taken into +!! account. The contribution to the heating of tke dissipation is computed. +!! +!! PURPOSE +!! ------- +! The purpose of this routine is to compute the sources necessary for +! the evolution of the turbulent kinetic energy and its dissipation +! if necessary. +! +!!** METHOD +!! ------ +!! The vertical turbulent flux is computed in an off-centered +!! implicit scheme (a Crank-Nicholson type with coefficients different +!! than 0.5), which allows to vary the degree of implicitness of the +!! formulation. +!! In high resolution, the horizontal transport terms are also +!! calculated, but explicitly. +!! The evolution of the dissipation as a variable is made if +!! the parameter HTURBLEN is set equal to KEPS. The same reasoning +!! made for TKE applies. +!! +!! EXTERNAL +!! -------- +!! GX_U_M,GY_V_M,GZ_W_M +!! GX_M_U,GY_M_V : Cartesian vertical gradient operators +!! +!! MXF,MXM.MYF,MYM,MZF,MZM: Shuman functions (mean operators) +!! DZF : Shuman functions (difference operators) +!! +!! +!! IMPLICIT ARGUMENTS +!! ------------------ +!! Module MODD_CST : contains physical constants +!! +!! XG : gravity constant +!! +!! Module MODD_CTURB: contains the set of constants for +!! the turbulence scheme +!! +!! XCET,XCED : transport and dissipation cts. for the TKE +!! XCDP,XCDD,XCDT: constants from the parameterization of +!! the K-epsilon equation +!! XTKEMIN,XEPSMIN : minimum values for the TKE and its +!! dissipation +!! +!! Module MODD_PARAMETERS: +!! +!! JPVEXT +!! Module MODD_BUDGET: +!! NBUMOD : model in which budget is calculated +!! CBUTYPE : type of desired budget +!! 'CART' for cartesian box configuration +!! 'MASK' for budget zone defined by a mask +!! 'NONE' ' for no budget +!! LBU_RTKE : logical for budget of RTKE (turbulent kinetic energy) +!! .TRUE. = budget of RTKE +!! .FALSE. = no budget of RTKE +!! +!! +!! REFERENCE +!! --------- +!! Book 2 of documentation (routine TKE_EPS_SOURCES) +!! Book 1 of documentation (Chapter: Turbulence) +!! +!! AUTHOR +!! ------ +!! Joan Cuxart * INM and Meteo-France * +!! +!! MODIFICATIONS +!! ------------- +!! Original August 23, 1994 +!! Modifications: Feb 14, 1995 (J.Cuxart and J.Stein) +!! Doctorization and Optimization +!! June 29, 1995 (J.Stein) TKE budget +!! June 28, 1995 (J.Cuxart) Add LES tools +!! Modifications: February 29, 1996 (J. Stein) optimization +!! Modifications: May 6, 1996 (N. Wood) Extend some loops over +!! the outer points +!! Modifications: August 30, 1996 (P. Jabouille) calcul ZFLX at the +!! IKU level +!! October 10, 1996 (J.Stein) set Keff at t-deltat +!! Oct 8, 1996 (Cuxart,Sanchez) Var.LES: XETR_TF,XDISS_TF +!! December 20, 1996 (J.-P. Pinty) update the CALL BUDGET +!! November 24, 1997 (V. Masson) bug in <v'e> +!! removes the DO loops +!! Augu. 9, 1999 (J.Stein) TKE budget correction +!! Mar 07 2001 (V. Masson and J. Stein) remove the horizontal +!! turbulent transports of Tke computation +!! Nov 06, 2002 (V. Masson) LES budgets +!! July 20, 2003 (J.-P. Pinty P Jabouille) add the dissipative heating +!! May 2006 Remove KEPS +!! October 2009 (G. Tanguy) add ILENCH=LEN(YCOMMENT) after +!! change of YCOMMENT +!! 2012-02 Y. Seity, add possibility to run with reversed +!! vertical levels +!! 2015-01 (J. Escobar) missing get_halo(ZRES) for JPHEXT<> 1 +!! J.Escobar : 15/09/2015 : WENO5 & JPHEXT <> 1 +!! Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O +! P. Wautelet 20/05/2019: add name argument to ADDnFIELD_ll + new ADD4DFIELD_ll subroutine +! P. Wautelet 02/2020: use the new data structures and subroutines for budgets +! -------------------------------------------------------------------------- +! +!* 0. DECLARATIONS +! ------------ +! +USE MODD_ARGSLIST_ll, ONLY: LIST_ll +use modd_budget, only: lbudget_tke, lbudget_th, NBUDGET_TKE, NBUDGET_TH, tbudgets +USE MODD_CONF +USE MODD_CST +USE MODD_CTURB +USE MODD_DIAG_IN_RUN, ONLY: LDIAG_IN_RUN, XCURRENT_TKE_DISS +use modd_field, only: tfielddata, TYPEREAL +USE MODD_IO, ONLY: TFILEDATA +USE MODD_LES +USE MODD_PARAMETERS +! +USE MODE_ARGSLIST_ll, ONLY: ADD3DFIELD_ll, CLEANLIST_ll +use mode_budget, only: Budget_store_add, Budget_store_end, Budget_store_init +USE MODE_EXCHANGE_ll, ONLY: UPDATE_HALO_ll +USE MODE_IO_FIELD_WRITE, only: IO_Field_write +#ifdef MNH_OPENACC +USE MODE_MNH_ZWORK, ONLY: MNH_MEM_GET, MNH_MEM_POSITION_PIN, MNH_MEM_RELEASE +#endif +use mode_mppdb +! +#ifdef MNH_BITREP +USE MODI_BITREP +#endif +USE MODI_GET_HALO +USE MODI_GRADIENT_M +USE MODI_GRADIENT_U +USE MODI_GRADIENT_V +USE MODI_GRADIENT_W +USE MODI_LES_MEAN_SUBGRID +#ifndef MNH_OPENACC +USE MODI_SHUMAN +#else +USE MODI_SHUMAN_DEVICE +#endif +USE MODI_TRIDIAG_TKE +! +IMPLICIT NONE +! +! +!* 0.1 declarations of 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) :: KMI ! model index number +REAL, DIMENSION(:,:,:), INTENT(IN) :: PTKEM ! TKE at t-deltat +REAL, DIMENSION(:,:,:), INTENT(IN) :: PLM ! mixing length +REAL, DIMENSION(:,:,:), INTENT(IN) :: PLEPS ! dissipative length +REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODJ ! density * grid volume +REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX,PDYY,PDZZ,PDZX,PDZY + ! metric coefficients +REAL, DIMENSION(:,:,:), INTENT(IN) :: PZZ ! physical height w-pt +REAL, INTENT(IN) :: PTSTEP ! Time step +REAL, INTENT(IN) :: PEXPL, PIMPL ! Coef. temporal. disc. +CHARACTER(len=4), INTENT(IN) :: HTURBDIM ! dimensionality of the + ! turbulence scheme +CHARACTER(len=4), INTENT(IN) :: HTURBLEN ! kind of mixing length +TYPE(TFILEDATA), INTENT(INOUT):: TPFILE ! Output file +LOGICAL, INTENT(IN) :: OTURB_DIAG ! switch to write some + ! diagnostic fields in the syncronous FM-file +REAL, DIMENSION(:,:,:), INTENT(INOUT):: PDP ! Dyn. prod. of TKE +REAL, DIMENSION(:,:,:), INTENT(IN) :: PTRH +REAL, DIMENSION(:,:,:), INTENT(IN) :: PTP ! Ther. prod. of TKE +REAL, DIMENSION(:,:,:), INTENT(INOUT):: PRTKES ! RHOD * Jacobian * + ! TKE at t+deltat +REAL, DIMENSION(:,:,:), INTENT(INOUT):: PRTHLS ! Source of Theta_l +REAL, DIMENSION(:,:,:), INTENT(IN) :: PCOEF_DISS ! 1/(Cph*Exner) +REAL, DIMENSION(:,:,:), INTENT(IN) :: PRTKESM ! Advection source +REAL, DIMENSION(:,:,:), INTENT(OUT) :: PTR ! Transport prod. of TKE +REAL, DIMENSION(:,:,:), INTENT(OUT) :: PDISS ! Dissipati prod. of TKE +! +!* 0.2 declaration of local variables +! +REAL, DIMENSION(:,:,:), pointer , contiguous :: & + ZA, & ! under diagonal elements of the tri-diagonal matrix involved + ! in the temporal implicit scheme + ZRES, & ! treated variable at t+ deltat when the turbu- + ! lence is the only source of evolution added to the ones + ! considered in ZSOURCE. This variable is also used to + ! temporarily store some diagnostics stored in FM file + ZFLX, & ! horizontal or vertical flux of the treated variable + ZSOURCE, & ! source of evolution for the treated variable + ZKEFF ! effectif diffusion coeff = LT * SQRT( TKE ) +INTEGER :: IIB,IIE,IJB,IJE,IKB,IKE + ! Index values for the Beginning and End + ! mass points of the domain +INTEGER :: IIU,IJU,IKU ! array size in the 3 dimensions +! +TYPE(LIST_ll), POINTER :: TZFIELDDISS_ll ! list of fields to exchange +INTEGER :: IINFO_ll ! return code of parallel routine +TYPE(TFIELDDATA) :: TZFIELD +! +#ifdef MNH_OPENACC +REAL, DIMENSION(:,:,:), pointer , contiguous :: ZTMP1_DEVICE, ZTMP2_DEVICE, ZTMP3_DEVICE, ZTMP4_DEVICE +#endif +! +INTEGER :: JIU,JJU,JKU +INTEGER :: JI,JJ,JK +!---------------------------------------------------------------------------- + +!$acc data present( PTKEM, PLM, PLEPS, PDP, PTRH, & +!$acc & PRHODJ, PDZZ, PDXX, PDYY, PDZX, PDZY, PZZ, & +!$acc & PTP, PRTKES, PRTKESM, PRTHLS, PCOEF_DISS, & +!$acc & PTR, PDISS ) + +if ( mppdb_initialized ) then + !Check all in arrays + call Mppdb_check( ptkem, "Tke_eps_sources beg:ptkem" ) + call Mppdb_check( plm, "Tke_eps_sources beg:plm" ) + call Mppdb_check( pleps, "Tke_eps_sources beg:pleps" ) + call Mppdb_check( prhodj, "Tke_eps_sources beg:prhodj" ) + call Mppdb_check( pdxx, "Tke_eps_sources beg:pdxx" ) + call Mppdb_check( pdyy, "Tke_eps_sources beg:pdyy" ) + call Mppdb_check( pdzz, "Tke_eps_sources beg:pdzz" ) + call Mppdb_check( pdzx, "Tke_eps_sources beg:pdzx" ) + call Mppdb_check( pdzy, "Tke_eps_sources beg:pdzy" ) + call Mppdb_check( pzz, "Tke_eps_sources beg:pzz" ) + call Mppdb_check( ptrh, "Tke_eps_sources beg:ptrh" ) + call Mppdb_check( ptp, "Tke_eps_sources beg:ptp" ) + call Mppdb_check( pcoef_diss, "Tke_eps_sources beg:pcoef_diss" ) + call Mppdb_check( prtkesm, "Tke_eps_sources beg:prtkesm" ) + !check all inout arrays + call Mppdb_check( pdp, "Tke_eps_sources beg:pdp" ) + call Mppdb_check( prtkes, "Tke_eps_sources beg:prtkes" ) + call Mppdb_check( prthls, "Tke_eps_sources beg:prthls" ) +end if + +JIU = size(ptkem, 1 ) +JJU = size(ptkem, 2 ) +JKU = size(ptkem, 3 ) + +#ifndef MNH_OPENACC +allocate( za (JIU,JJU,JKU ) ) +allocate( zres (JIU,JJU,JKU ) ) +allocate( zflx (JIU,JJU,JKU ) ) +allocate( zsource(JIU,JJU,JKU ) ) +allocate( zkeff (JIU,JJU,JKU ) ) +#else +!Pin positions in the pools of MNH memory +CALL MNH_MEM_POSITION_PIN() + +CALL MNH_MEM_GET( za ,JIU,JJU,JKU) +CALL MNH_MEM_GET( zres ,JIU,JJU,JKU) +CALL MNH_MEM_GET( zflx ,JIU,JJU,JKU) +CALL MNH_MEM_GET( zsource,JIU,JJU,JKU) +CALL MNH_MEM_GET( zkeff ,JIU,JJU,JKU) + +CALL MNH_MEM_GET( ztmp1_device,JIU,JJU,JKU) +CALL MNH_MEM_GET( ztmp2_device,JIU,JJU,JKU) +CALL MNH_MEM_GET( ztmp3_device,JIU,JJU,JKU) +CALL MNH_MEM_GET( ztmp4_device,JIU,JJU,JKU) +#endif + +!$acc data present( ZA, ZRES, ZFLX, ZSOURCE, ZKEFF, ZTMP1_DEVICE, ZTMP2_DEVICE, ZTMP3_DEVICE, ZTMP4_DEVICE ) + +NULLIFY(TZFIELDDISS_ll) +! +!* 1. PRELIMINARY COMPUTATIONS +! ------------------------ +! +CALL GET_INDICE_ll (IIB,IJB,IIE,IJE) +!$acc kernels +IIU=SIZE(PTKEM,1) +IJU=SIZE(PTKEM,2) +IKB=KKA+JPVEXT_TURB*KKL +IKE=KKU-JPVEXT_TURB*KKL +! +! compute the effective diffusion coefficient at the mass point +#ifndef MNH_BITREP +ZKEFF(:,:,:) = PLM(:,:,:) * SQRT(PTKEM(:,:,:)) +#else +#ifdef MNH_COMPILER_NVHPC +!$acc loop independent collapse(3) +#endif +DO CONCURRENT ( JI=1:JIU,JJ=1:JJU,JK=1:JKU) + ZKEFF(JI,JJ,JK) = PLM(JI,JJ,JK) * BR_POW(PTKEM(JI,JJ,JK),0.5) +END DO +#endif +! +!$acc end kernels +if (lbudget_th) call Budget_store_init( tbudgets(NBUDGET_TH), 'DISSH', prthls(:, :, :) ) +!$acc kernels +!---------------------------------------------------------------------------- +! +!* 2. TKE EQUATION +! ------------ +! +!* 2.1 Horizontal turbulent explicit transport +! +! +! Complete the sources of TKE with the horizontal turbulent explicit transport +! +IF (HTURBDIM=='3DIM') THEN + PTR(:,:,:) = PTRH(:,:,:) +ELSE + PTR(:,:,:) = 0. +END IF +! +! +! +!* 2.2 Explicit TKE sources except horizontal turbulent transport +! +! +! extrapolate the dynamic production with a 1/Z law from its value at the +! W(IKB+1) value stored in PDP(IKB) to the mass localization tke(IKB) +PDP(:,:,IKB) = PDP(:,:,IKB) * (1. + PDZZ(:,:,IKB+KKL)/PDZZ(:,:,IKB)) +! +! Compute the source terms for TKE: ( ADVECtion + NUMerical DIFFusion + ..) +! + (Dynamical Production) + (Thermal Production) - (dissipation) +#ifndef MNH_BITREP +ZFLX(:,:,:) = XCED * SQRT(PTKEM(:,:,:)) / PLEPS(:,:,:) +#else +#ifdef MNH_COMPILER_NVHPC +!$acc loop independent collapse(3) +#endif +DO CONCURRENT ( JI=1:JIU,JJ=1:JJU,JK=1:JKU) + ZFLX(JI,JJ,JK) = XCED * BR_POW(PTKEM(JI,JJ,JK),0.5) / PLEPS(JI,JJ,JK) +END DO +#endif +#ifdef MNH_COMPILER_NVHPC +!$acc loop independent collapse(3) +#endif +DO CONCURRENT ( JI=1:JIU,JJ=1:JJU,JK=1:JKU) + ZSOURCE(JI,JJ,JK) = ( PRTKES(JI,JJ,JK) + PRTKESM(JI,JJ,JK) ) / PRHODJ(JI,JJ,JK) & + - PTKEM(JI,JJ,JK) / PTSTEP & + + PDP(JI,JJ,JK) + PTP(JI,JJ,JK) + PTR(JI,JJ,JK) - PEXPL * ZFLX(JI,JJ,JK) * PTKEM(JI,JJ,JK) +END DO !CONCURRENT +!$acc end kernels +! +!* 2.2 implicit vertical TKE transport +! +! +! Compute the vector giving the elements just under the diagonal for the +! matrix inverted in TRIDIAG +! +#ifndef MNH_OPENACC +ZA(:,:,:) = - PTSTEP * XCET * & +#ifndef MNH_BITREP + MZM(ZKEFF) * MZM(PRHODJ) / PDZZ(:,:,:)**2 +#else + MZM(ZKEFF) * MZM(PRHODJ) / BR_P2(PDZZ(:,:,:)) +#endif +#else +CALL MZM_DEVICE(ZKEFF, ZTMP1_DEVICE) !Warning: re-used later +CALL MZM_DEVICE(PRHODJ,ZTMP2_DEVICE) !Warning: re-used later +!$acc kernels +#ifndef MNH_BITREP +ZA(:,:,:) = - PTSTEP * XCET * ZTMP1_DEVICE(:,:,:) * ZTMP2_DEVICE(:,:,:) / PDZZ(:,:,:)**2 +#else +#ifdef MNH_COMPILER_NVHPC +!$acc loop independent collapse(3) +#endif +DO CONCURRENT ( JI=1:JIU,JJ=1:JJU,JK=1:JKU) + ZA(JI,JJ,JK) = - PTSTEP * XCET * ZTMP1_DEVICE(JI,JJ,JK) * ZTMP2_DEVICE(JI,JJ,JK) / BR_P2(PDZZ(JI,JJ,JK)) +END DO !CONCURRENT +#endif +!$acc end kernels +#endif +! +! Compute TKE at time t+deltat: ( stored in ZRES ) +! +#ifndef MNH_OPENACC +CALL TRIDIAG_TKE(KKA,KKU,KKL,PTKEM,ZA,PTSTEP,PEXPL,PIMPL,PRHODJ,& + & ZSOURCE,PTSTEP*ZFLX,ZRES) +CALL GET_HALO(ZRES) +#else +!$acc kernels +#ifdef MNH_COMPILER_NVHPC +!$acc loop independent collapse(3) +#endif +DO CONCURRENT ( JI=1:JIU,JJ=1:JJU,JK=1:JKU) + ZTMP3_DEVICE(JI,JJ,JK) = PTSTEP*ZFLX(JI,JJ,JK) +END DO !CONCURRENT +!$acc end kernels +CALL TRIDIAG_TKE(KKA,KKU,KKL,PTKEM,ZA,PTSTEP,PEXPL,PIMPL,PRHODJ,& + & ZSOURCE,ZTMP3_DEVICE,ZRES) +! acc update self(ZRES) +CALL GET_HALO_D(ZRES) +! acc update device(ZRES) +#endif +! +!* diagnose the dissipation +! +IF (LDIAG_IN_RUN) THEN +!$acc data copyout(XCURRENT_TKE_DISS) + XCURRENT_TKE_DISS = ZFLX(:,:,:) * PTKEM(:,:,:) & + *(PEXPL*PTKEM(:,:,:) + PIMPL*ZRES(:,:,:)) +!$acc end data + CALL ADD3DFIELD_ll( TZFIELDDISS_ll, XCURRENT_TKE_DISS, 'TKE_EPS_SOURCES::XCURRENT_TKE_DISS' ) + CALL UPDATE_HALO_ll(TZFIELDDISS_ll,IINFO_ll) + CALL CLEANLIST_ll(TZFIELDDISS_ll) +ENDIF +! +! TKE must be greater than its minimum value +! +! CL : Now done at the end of the time step in ADVECTION_METSV +!GTKENEG = ZRES <= XTKEMIN +!WHERE ( GTKENEG ) +! ZRES = XTKEMIN +!END WHERE +! +IF ( LLES_CALL .OR. & + (OTURB_DIAG .AND. tpfile%lopened) ) THEN +! +! Compute the cartesian vertical flux of TKE in ZFLX +! +#ifndef MNH_OPENACC + ZFLX(:,:,:) = - XCET * MZM(ZKEFF) * & + DZM(PIMPL * ZRES + PEXPL * PTKEM ) / PDZZ(:,:,:) +#else +!$acc kernels + ZTMP3_DEVICE(:,:,:) = PIMPL * ZRES(:,:,:) + PEXPL * PTKEM(:,:,:) +!$acc end kernels + CALL DZM_DEVICE( ZTMP3_DEVICE, ZTMP4_DEVICE ) +!$acc kernels + ZFLX(:,:,:) = - XCET * ZTMP1_DEVICE(:,:,:) * ZTMP4_DEVICE(:,:,:) / PDZZ(:,:,:) !Re-use of ZTMP1_DEVICE +#endif +! + ZFLX(:,:,IKB) = 0. + ZFLX(:,:,KKA) = 0. +! +! Compute the whole turbulent TRansport of TKE: +! +#ifndef MNH_OPENACC + PTR(:,:,:)= PTR - DZF( MZM(PRHODJ) * ZFLX / PDZZ ) /PRHODJ +#else + ZTMP1_DEVICE(:,:,:) = ZTMP2_DEVICE(:,:,:) * ZFLX(:,:,:) / PDZZ(:,:,:) !Re-use of ZTMP2_DEVICE +!$acc end kernels + CALL DZF_DEVICE( ZTMP1_DEVICE, ZTMP2_DEVICE ) +!$acc kernels + PTR(:,:,:)= PTR(:,:,:) - ZTMP2_DEVICE(:,:,:) / PRHODJ(:,:,:) +!$acc end kernels +#endif +! +! Storage in the LES configuration +! + IF (LLES_CALL) THEN +#ifndef MNH_OPENACC + CALL LES_MEAN_SUBGRID( MZF(ZFLX), X_LES_SUBGRID_WTke ) + CALL LES_MEAN_SUBGRID( -PTR, X_LES_SUBGRID_ddz_WTke ) +#else +!$acc data copy(X_LES_SUBGRID_WTke,X_LES_SUBGRID_ddz_WTke) + CALL MZF_DEVICE( ZFLX, ZTMP1_DEVICE ) + CALL LES_MEAN_SUBGRID( ZTMP1_DEVICE, X_LES_SUBGRID_WTke ) +!$acc kernels + ZTMP1_DEVICE(:,:,:) = -PTR(:,:,:) +!$acc end kernels + CALL LES_MEAN_SUBGRID( ZTMP1_DEVICE, X_LES_SUBGRID_ddz_WTke ) +!$acc end data +#endif + END IF +! +END IF +! +!* 2.4 stores the explicit sources for budget purposes +! +if (lbudget_tke) then + ! Dynamical production + call Budget_store_add( tbudgets(NBUDGET_TKE), 'DP', pdp(:, :, :) * prhodj(:, :, :) ) + ! Thermal production + call Budget_store_add( tbudgets(NBUDGET_TKE), 'TP', ptp(:, :, :) * prhodj(:, :, :) ) + ! Dissipation + call Budget_store_add( tbudgets(NBUDGET_TKE), 'DISS', -xced * sqrt( ptkem(:, :, :) ) / pleps(:, :, :) & + * ( pexpl * ptkem(:, :, :) + pimpl * zres(:, :, :) ) * prhodj(:, :, :) ) +end if +! +!* 2.5 computes the final RTKE and stores the whole turbulent transport +! with the removal of the advection part + +if (lbudget_tke) then + !Store the previous source terms in prtkes before initializing the next one + !$acc kernels + PRTKES(:,:,:) = PRTKES(:,:,:) + PRHODJ(:,:,:) * & + ( PDP(:,:,:) + PTP(:,:,:) & + - XCED * SQRT(PTKEM(:,:,:)) / PLEPS(:,:,:) * ( PEXPL*PTKEM(:,:,:) + PIMPL*ZRES(:,:,:) ) ) + !$acc end kernels + call Budget_store_init( tbudgets(NBUDGET_TKE), 'TR', prtkes(:, :, :) ) +end if + +!$acc kernels +PRTKES(:,:,:) = ZRES(:,:,:) * PRHODJ(:,:,:) / PTSTEP - PRTKESM(:,:,:) +!$acc end kernels +! +! stores the whole turbulent transport +! +if (lbudget_tke) call Budget_store_end( tbudgets(NBUDGET_TKE), 'TR', prtkes(:, :, :) ) + +!---------------------------------------------------------------------------- +! +!* 3. COMPUTE THE DISSIPATIVE HEATING +! ------------------------------- +! +!$acc kernels +#ifndef MNH_BITREP +PRTHLS(:,:,:) = PRTHLS(:,:,:) + XCED * SQRT(PTKEM(:,:,:)) / PLEPS(:,:,:) * & +#else +PRTHLS(:,:,:) = PRTHLS(:,:,:) + XCED * BR_POW(PTKEM(:,:,:),0.5) / PLEPS(:,:,:) * & +#endif + (PEXPL*PTKEM(:,:,:) + PIMPL*ZRES(:,:,:)) * PRHODJ(:,:,:) * PCOEF_DISS(:,:,:) +!$acc end kernels + +if (lbudget_th) call Budget_store_end( tbudgets(NBUDGET_TH), 'DISSH', prthls(:, :, :) ) + +!---------------------------------------------------------------------------- +! +!* 4. STORES SOME DIAGNOSTICS +! ----------------------- +! +!$acc kernels +#ifndef MNH_BITREP +PDISS(:,:,:) = -XCED * (PTKEM(:,:,:)**1.5) / PLEPS(:,:,:) +#else +PDISS(:,:,:) = -XCED * BR_POW(PTKEM(:,:,:),1.5) / PLEPS(:,:,:) +#endif +!$acc end kernels +! +IF ( OTURB_DIAG .AND. tpfile%lopened ) THEN +!$acc update self(PDP,PTP,PTR,PDISS) +! +! stores the dynamic production +! + TZFIELD%CMNHNAME = 'DP' + TZFIELD%CSTDNAME = '' + TZFIELD%CLONGNAME = 'DP' + TZFIELD%CUNITS = 'm2 s-3' + TZFIELD%CDIR = 'XY' + TZFIELD%CCOMMENT = 'X_Y_Z_DP' + TZFIELD%NGRID = 1 + TZFIELD%NTYPE = TYPEREAL + TZFIELD%NDIMS = 3 + TZFIELD%LTIMEDEP = .TRUE. + CALL IO_Field_write(TPFILE,TZFIELD,PDP) +! +! stores the thermal production +! + TZFIELD%CMNHNAME = 'TP' + TZFIELD%CSTDNAME = '' + TZFIELD%CLONGNAME = 'TP' + TZFIELD%CUNITS = 'm2 s-3' + TZFIELD%CDIR = 'XY' + TZFIELD%CCOMMENT = 'X_Y_Z_TP' + TZFIELD%NGRID = 1 + TZFIELD%NTYPE = TYPEREAL + TZFIELD%NDIMS = 3 + TZFIELD%LTIMEDEP = .TRUE. + CALL IO_Field_write(TPFILE,TZFIELD,PTP) +! +! stores the whole turbulent transport +! + TZFIELD%CMNHNAME = 'TR' + TZFIELD%CSTDNAME = '' + TZFIELD%CLONGNAME = 'TR' + TZFIELD%CUNITS = 'm2 s-3' + TZFIELD%CDIR = 'XY' + TZFIELD%CCOMMENT = 'X_Y_Z_TR' + TZFIELD%NGRID = 1 + TZFIELD%NTYPE = TYPEREAL + TZFIELD%NDIMS = 3 + TZFIELD%LTIMEDEP = .TRUE. + CALL IO_Field_write(TPFILE,TZFIELD,PTR) +! +! stores the dissipation of TKE +! + TZFIELD%CMNHNAME = 'DISS' + TZFIELD%CSTDNAME = '' + TZFIELD%CLONGNAME = 'DISS' + TZFIELD%CUNITS = 'm2 s-3' + TZFIELD%CDIR = 'XY' + TZFIELD%CCOMMENT = 'X_Y_Z_DISS' + TZFIELD%NGRID = 1 + TZFIELD%NTYPE = TYPEREAL + TZFIELD%NDIMS = 3 + TZFIELD%LTIMEDEP = .TRUE. + CALL IO_Field_write(TPFILE,TZFIELD,PDISS) +END IF +! +! Storage in the LES configuration of the Dynamic Production of TKE and +! the dissipation of TKE +! +IF (LLES_CALL ) THEN +!$acc data copy(X_LES_SUBGRID_DISS_Tke) + CALL LES_MEAN_SUBGRID( PDISS, X_LES_SUBGRID_DISS_Tke ) +!$acc end data +END IF + +if ( mppdb_initialized ) then + !check all inout arrays + call Mppdb_check( pdp, "Tke_eps_sources end:pdp" ) + call Mppdb_check( prtkes, "Tke_eps_sources end:prtkes" ) + call Mppdb_check( prthls, "Tke_eps_sources end:prthls" ) + !check all out arrays + call Mppdb_check( ptr, "Tke_eps_sources end:ptr" ) + call Mppdb_check( pdiss, "Tke_eps_sources end:pdiss" ) +end if + +!$acc end data + +#ifndef MNH_OPENACC +deallocate( ZA, ZRES, ZFLX, ZSOURCE, ZKEFF ) +#else +!Release all memory allocated with MNH_MEM_GET calls since last call to MNH_MEM_POSITION_PIN +CALL MNH_MEM_RELEASE() +#endif + +!$acc end data + +END SUBROUTINE TKE_EPS_SOURCES