diff --git a/src/ZSOLVER/ice_adjust.f90 b/src/ZSOLVER/ice_adjust.f90
deleted file mode 100644
index d19aa6b93f9c21a966943bdc63c3e4ab4343f274..0000000000000000000000000000000000000000
--- a/src/ZSOLVER/ice_adjust.f90
+++ /dev/null
@@ -1,738 +0,0 @@
-!MNH_LIC Copyright 1996-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_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
-!
-#ifndef MNH_OPENACC
-LOGICAL,DIMENSION(:,:,:), allocatable :: GTEMP
-!
-REAL, DIMENSION(:,:,:), allocatable :: ZSIGS,ZSRCS
-REAL, DIMENSION(:,:,:), allocatable &
- :: 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
-#else
-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
-REAL, DIMENSION(:,:,:), pointer , contiguous :: ZTEMP_BUD
-#endif
-!
-LOGICAL :: GPOUT_RV,GPOUT_RC,GPOUT_RI,GPOUT_TH
-!-------------------------------------------------------------------------------
-!
-GPOUT_RV = PRESENT(POUT_RV)
-GPOUT_RC = PRESENT(POUT_RC)
-GPOUT_RI = PRESENT(POUT_RI)
-GPOUT_TH = PRESENT(POUT_TH)
-!
-! 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 ) )
-allocate( ZTEMP_BUD (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 )
-CALL MNH_MEM_GET( ZTEMP_BUD , 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 ) then
- !$acc kernels
- ZTEMP_BUD(:,:,:) = pths(:, :, :) * prhodj(:, :, :)
- !$acc end kernels
- call Budget_store_init( tbudgets(NBUDGET_TH), trim( hbuname ), ZTEMP_BUD(:,:,:) )
-end if
-if ( lbudget_rv ) then
- !$acc kernels
- ZTEMP_BUD(:,:,:) = prvs(:, :, :) * prhodj(:, :, :)
- !$acc end kernels
- call Budget_store_init( tbudgets(NBUDGET_RV), trim( hbuname ), ZTEMP_BUD(:,:,:) )
-end if
-if ( lbudget_rc ) then
- !$acc kernels
- ZTEMP_BUD(:,:,:) = prcs(:, :, :) * prhodj(:, :, :)
- !$acc end kernels
- call Budget_store_init( tbudgets(NBUDGET_RC), trim( hbuname ), ZTEMP_BUD(:,:,:) )
-end if
-if ( lbudget_ri ) then
- !$acc kernels
- ZTEMP_BUD(:,:,:) = pris(:, :, :) * prhodj(:, :, :)
- !$acc end kernels
- call Budget_store_init( tbudgets(NBUDGET_RI), trim( hbuname ), ZTEMP_BUD(:,:,:) )
-end if
-!
-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 present_cr(ZLV,ZLS,ZCPH)
- 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 present_cr(GTEMP)
-DO CONCURRENT (JI=1:IIU,JJ=1:IJU,JK=1:IKU)
- GTEMP(JI,JJ,JK) = PRCS(JI,JJ,JK) + PRIS(JI,JJ,JK) > 1.E-12 / PTSTEP
- IF ( GTEMP(JI,JJ,JK) )THEN
- PCLDFR(JI,JJ,JK) = 1.
- ELSE
- PCLDFR(JI,JJ,JK) = 0.
- ENDIF
-ENDDO
- 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(GPOUT_RV .OR. GPOUT_RC .OR. &
- &GPOUT_RI .OR. GPOUT_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(GPOUT_RV) POUT_RV=ZRV
-IF(GPOUT_RC) POUT_RC=ZRC
-IF(GPOUT_RI) POUT_RI=ZRI
-IF(GPOUT_TH) POUT_TH=ZT / PEXN(:,:,:)
-!$acc end kernels
-!
-!
-!* 6. STORE THE BUDGET TERMS
-! ----------------------
-!
-if ( lbudget_th ) then
- !$acc kernels
- ZTEMP_BUD(:,:,:) = pths(:, :, :) * prhodj(:, :, :)
- !$acc end kernels
- call Budget_store_end( tbudgets(NBUDGET_TH), trim( hbuname ), ZTEMP_BUD(:,:,:) )
-end if
-if ( lbudget_rv ) then
- !$acc kernels
- ZTEMP_BUD(:,:,:) = prvs(:, :, :) * prhodj(:, :, :)
- !$acc end kernels
- call Budget_store_end( tbudgets(NBUDGET_RV), trim( hbuname ), ZTEMP_BUD(:,:,:) )
-end if
-if ( lbudget_rc ) then
- !$acc kernels
- ZTEMP_BUD(:,:,:) = prcs(:, :, :) * prhodj(:, :, :)
- !$acc end kernels
- call Budget_store_end( tbudgets(NBUDGET_RC), trim( hbuname ), ZTEMP_BUD(:,:,:) )
-end if
-if ( lbudget_ri ) then
- !$acc kernels
- ZTEMP_BUD(:,:,:) = pris(:, :, :) * prhodj(:, :, :)
- !$acc end kernels
- call Budget_store_end( tbudgets(NBUDGET_RI), trim( hbuname ), ZTEMP_BUD(:,:,:) )
-end if
-
-!$acc end data
-
-#ifndef MNH_OPENACC
-deallocate( gtemp, zsigs, zsrcs, zt, zrv, zrc, zri, zcph, zlv, zls, zw1, zw2, zcriaut, zhcf, zhr , ZTEMP_BUD )
-#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 (GPOUT_RV) CALL MPPDB_CHECK3D(POUT_RV,"ICE_ADJUST end:POUT_RV")
-IF (GPOUT_RC) CALL MPPDB_CHECK3D(POUT_RC,"ICE_ADJUST end:POUT_RC")
-IF (GPOUT_RI) CALL MPPDB_CHECK3D(POUT_RI,"ICE_ADJUST end:POUT_RI")
-IF (GPOUT_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
diff --git a/src/ZSOLVER/sources_neg_correct.f90 b/src/ZSOLVER/sources_neg_correct.f90
deleted file mode 100644
index 5602b75ccf615817e299c332cf382beb9053c7f6..0000000000000000000000000000000000000000
--- a/src/ZSOLVER/sources_neg_correct.f90
+++ /dev/null
@@ -1,580 +0,0 @@
-!MNH_LIC Copyright 2020-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.
-!-----------------------------------------------------------------
-! Author: P. Wautelet 25/06/2020 (deduplication of code from advection_metsv, resolved_cloud and turb)
-! Modifications:
-! P. Wautelet 30/06/2020: remove non-local corrections in resolved_cloud for NEGA => new local corrections here
-! J. Escobar 21/07/2020: bug <-> array of size(:,:,:,0) => return if krr=0
-! P. Wautelet 10/02/2021: budgets: add missing sources for NSV_C2R2BEG+3 budget
-!-----------------------------------------------------------------
-module mode_sources_neg_correct
-
-implicit none
-
-private
-
-public :: Sources_neg_correct
-
-contains
-
-subroutine Sources_neg_correct( hcloud, hbudname, krr, ptstep, ppabst, ptht, prt, prths, prrs, prsvs, prhodj )
-
-use modd_budget, only: lbudget_th, lbudget_rv, lbudget_rc, lbudget_rr, lbudget_ri, &
- lbudget_rs, lbudget_rg, lbudget_rh, lbudget_sv, &
- NBUDGET_TH, NBUDGET_RV, NBUDGET_RC, NBUDGET_RR, NBUDGET_RI, &
- NBUDGET_RS, NBUDGET_RG, NBUDGET_RH, NBUDGET_SV1, &
- tbudgets
-use modd_cst, only: xci, xcl, xcpd, xcpv, xlstt, xlvtt, xp00, xrd, xtt
-use modd_nsv, only: nsv_c2r2beg, nsv_c2r2end, nsv_lima_beg, nsv_lima_end, nsv_lima_nc, nsv_lima_nr, nsv_lima_ni
-use modd_param_lima, only: lcold_lima => lcold, lrain_lima => lrain, lspro_lima => lspro, lwarm_lima => lwarm, &
- xctmin_lima => xctmin, xrtmin_lima => xrtmin
-
-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
-use mode_msg
-
-#ifdef MNH_BITREP
-use modi_bitrep
-#endif
-
-implicit none
-
-character(len=*), intent(in) :: hcloud ! Kind of cloud parameterization
-character(len=*), intent(in) :: hbudname ! Budget name
-integer, intent(in) :: krr ! Number of moist variables
-real, intent(in) :: ptstep ! Timestep
-real, dimension(:, :, :), intent(in) :: ppabst ! Absolute pressure at time t
-real, dimension(:, :, :), intent(in) :: ptht ! Theta at time t
-real, dimension(:, :, :, :), intent(in) :: prt ! Moist variables at time t
-real, dimension(:, :, :), intent(inout) :: prths ! Source terms
-real, dimension(:, :, :, :), intent(inout) :: prrs ! Source terms
-real, dimension(:, :, :, :), intent(inout) :: prsvs ! Source terms
-real, dimension(:, :, :), intent(in), optional :: prhodj ! Dry density * jacobian
-
-integer :: jiu, jju, jku
-integer :: ji, jj, jk
-integer :: jr
-integer :: jrmax
-integer :: jsv
-integer :: isv_lima_end
-real, dimension(:, :, :), pointer, contiguous :: zt, zexn, zlv, zls, zcph, zcor
-REAL, DIMENSION(:,:,:), pointer , contiguous :: ZTEMP_BUD
-
-if ( krr == 0 ) return
-
-zcor => null()
-
-jiu = Size(prths, 1 )
-jju = Size(prths, 2 )
-jku = Size(prths, 3 )
-
-#ifndef MNH_OPENACC
-allocate( zt ( jiu, jju, jku ) )
-allocate( zexn( jiu, jju, jku ) )
-allocate( zlv ( jiu, jju, jku ) )
-allocate( zcph( jiu, jju, jku ) )
-if ( hcloud == 'ICE3' .or. hcloud == 'ICE4' .or. hcloud == 'LIMA' ) then
- allocate( zls( jiu, jju, jku ) )
- if ( krr > 3 ) then
- allocate( zcor( jiu, jju, jku ) )
- end if
-else
- allocate( zls(0, 0, 0) )
-end if
-if ( .not. Associated( zcor ) ) Allocate( zcor(0, 0, 0) )
-ALLOCATE(ZTEMP_BUD( jiu, jju, jku ))
-#else
-!Pin positions in the pools of MNH memory
-call Mnh_mem_position_pin()
-
-call Mnh_mem_get( zt, jiu, jju, jku )
-call Mnh_mem_get( zexn, jiu, jju, jku )
-call Mnh_mem_get( zlv, jiu, jju, jku )
-call Mnh_mem_get( zcph, jiu, jju, jku )
-if ( hcloud == 'ICE3' .or. hcloud == 'ICE4' .or. hcloud == 'LIMA' ) then
- call Mnh_mem_get( zls, jiu, jju, jku )
- if ( krr > 3 ) then
- call Mnh_mem_get( zcor, jiu, jju, jku )
- else
- call Mnh_mem_get( zcor, 0, 0, 0 )
- end if
-else
- call Mnh_mem_get( zls, 0, 0, 0 )
- call Mnh_mem_get( zcor, 0, 0, 0 )
-end if
-call MNH_MEM_GET( ZTEMP_BUD, jiu, jju, jku )
-#endif
-
-!$acc data present( ppabst, ptht, prt, prths, prrs, prsvs, prhodj )
-
-if ( mppdb_initialized ) then
- !Check all IN arrays
- call Mppdb_check( ppabst, "Sources_neg_correct beg:ppabst")
- call Mppdb_check( ptht, "Sources_neg_correct beg:ptht")
- call Mppdb_check( prt, "Sources_neg_correct beg:prt")
- if ( Present( prhodj ) ) call Mppdb_check( prhodj, "Sources_neg_correct beg:prhodj")
- !Check all INOUT arrays
- call Mppdb_check( prths, "Sources_neg_correct beg:prths")
- call Mppdb_check( prrs, "Sources_neg_correct beg:prrs")
- call Mppdb_check( prsvs, "Sources_neg_correct beg:prsvs")
-end if
-
-if ( hbudname /= 'NEADV' .and. hbudname /= 'NECON' .and. hbudname /= 'NEGA' .and. hbudname /= 'NETUR' ) &
- call Print_msg( NVERB_WARNING, 'GEN', 'Sources_neg_correct', 'budget '//hbudname//' not yet tested' )
-
-if ( hcloud == 'LIMA' ) then
- ! The negativity correction does not apply to the SPRO (supersaturation) variable which may be naturally negative
- if ( lspro_lima ) then
- isv_lima_end = nsv_lima_end - 1
- else
- isv_lima_end = nsv_lima_end
- end if
-end if
-
-if ( hbudname /= 'NECON' .and. hbudname /= 'NEGA' ) then
- if ( hcloud == 'KESS' .or. hcloud == 'ICE3' .or. hcloud == 'ICE4' .or. &
- hcloud == 'KHKO' .or. hcloud == 'C2R2' .or. hcloud == 'LIMA' ) then
- if ( lbudget_th ) call Budget_store_init( tbudgets(NBUDGET_TH), Trim( hbudname ), prths(:, :, :) )
- if ( lbudget_rv ) call Budget_store_init( tbudgets(NBUDGET_RV), Trim( hbudname ), prrs (:, :, :, 1) )
- if ( lbudget_rc ) call Budget_store_init( tbudgets(NBUDGET_RC), Trim( hbudname ), prrs (:, :, :, 2) )
- if ( lbudget_rr .and. &
- ( hbudname /= 'NETUR' .or. &
- ( hbudname == 'NETUR' .and. ( hcloud == 'C2R2' .or. hcloud == 'KHKO' .or. hcloud == 'LIMA' ) ) ) ) &
- call Budget_store_init( tbudgets(NBUDGET_RR), Trim( hbudname ), prrs (:, :, :, 3) )
- IF (lbudget_ri .and. &
- ( hbudname /= 'NETUR' .or. &
- ( hbudname == 'NETUR' .and. ( hcloud == 'ICE3' .or. hcloud == 'ICE4' .or. hcloud == 'LIMA' ) ) ) ) &
- call Budget_store_init( tbudgets(NBUDGET_RI), Trim( hbudname ), prrs (:, :, :, 4) )
- if ( lbudget_rs .and. hbudname /= 'NETUR' ) call Budget_store_init( tbudgets(NBUDGET_RS), Trim( hbudname ), prrs (:, :, :, 5) )
- if ( lbudget_rg .and. hbudname /= 'NETUR' ) call Budget_store_init( tbudgets(NBUDGET_RG), Trim( hbudname ), prrs (:, :, :, 6) )
- if ( lbudget_rh .and. hbudname /= 'NETUR' ) call Budget_store_init( tbudgets(NBUDGET_RH), Trim( hbudname ), prrs (:, :, :, 7) )
- end if
-
- if ( lbudget_sv .and. ( hcloud == 'C2R2' .or. hcloud == 'KHKO' ) ) then
- do ji = nsv_c2r2beg, nsv_c2r2end
- call Budget_store_init( tbudgets(NBUDGET_SV1 - 1 + ji), Trim( hbudname ), prsvs(:, :, :, ji) )
- end do
- end if
- if ( lbudget_sv .and. hcloud == 'LIMA' ) then
- do ji = nsv_lima_beg, isv_lima_end
- call Budget_store_init( tbudgets(NBUDGET_SV1 - 1 + ji), Trim( hbudname ), prsvs(:, :, :, ji) )
- end do
- end if
-else !NECON + NEGA
- if ( .not. present( prhodj ) ) &
- call Print_msg( NVERB_FATAL, 'GEN', 'Sources_neg_correct', 'optional argument prhodj not present' )
-
- if ( hcloud == 'KESS' .or. hcloud == 'ICE3' .or. hcloud == 'ICE4' .or. &
- hcloud == 'KHKO' .or. hcloud == 'C2R2' .or. hcloud == 'LIMA' ) then
- if ( lbudget_th) then
- !$acc kernels present(ZTEMP_BUD)
- ZTEMP_BUD(:,:,:) = prths(:, :, :) * prhodj(:, :, :)
- !$acc end kernels
- call Budget_store_init( tbudgets(NBUDGET_TH), Trim( hbudname ), ZTEMP_BUD(:,:,:) )
- end if
- if ( lbudget_rv) then
- !$acc kernels present(ZTEMP_BUD)
- ZTEMP_BUD(:,:,:) = prrs (:, :, :, 1) * prhodj(:, :, :)
- !$acc end kernels
- call Budget_store_init( tbudgets(NBUDGET_RV), Trim( hbudname ), ZTEMP_BUD(:,:,:) )
- end if
- if ( lbudget_rc) then
- !$acc kernels present(ZTEMP_BUD)
- ZTEMP_BUD(:,:,:) = prrs (:, :, :, 2) * prhodj(:, :, :)
- !$acc end kernels
- call Budget_store_init( tbudgets(NBUDGET_RC), Trim( hbudname ), ZTEMP_BUD(:,:,:) )
- end if
- if ( lbudget_rr) then
- !$acc kernels present(ZTEMP_BUD)
- ZTEMP_BUD(:,:,:) = prrs (:, :, :, 3) * prhodj(:, :, :)
- !$acc end kernels
- call Budget_store_init( tbudgets(NBUDGET_RR), Trim( hbudname ), ZTEMP_BUD(:,:,:) )
- end if
- if ( lbudget_ri) then
- !$acc kernels present(ZTEMP_BUD)
- ZTEMP_BUD(:,:,:) = prrs (:, :, :, 4) * prhodj(:, :, :)
- !$acc end kernels
- call Budget_store_init( tbudgets(NBUDGET_RI), Trim( hbudname ), ZTEMP_BUD(:,:,:) )
- end if
- if ( lbudget_rs) then
- !$acc kernels present(ZTEMP_BUD)
- ZTEMP_BUD(:,:,:) = prrs (:, :, :, 5) * prhodj(:, :, :)
- !$acc end kernels
- call Budget_store_init( tbudgets(NBUDGET_RS), Trim( hbudname ), ZTEMP_BUD(:,:,:) )
- end if
- if ( lbudget_rg) then
- !$acc kernels present(ZTEMP_BUD)
- ZTEMP_BUD(:,:,:) = prrs (:, :, :, 6) * prhodj(:, :, :)
- !$acc end kernels
- call Budget_store_init( tbudgets(NBUDGET_RG), Trim( hbudname ), ZTEMP_BUD(:,:,:) )
- end if
- if ( lbudget_rh) then
- !$acc kernels present(ZTEMP_BUD)
- ZTEMP_BUD(:,:,:) = prrs (:, :, :, 7) * prhodj(:, :, :)
- !$acc end kernels
- call Budget_store_init( tbudgets(NBUDGET_RH), Trim( hbudname ), ZTEMP_BUD(:,:,:) )
- end if
- end if
-
- if ( lbudget_sv .and. ( hcloud == 'C2R2' .or. hcloud == 'KHKO' ) ) then
- do ji = nsv_c2r2beg, nsv_c2r2end
- !$acc kernels present(ZTEMP_BUD)
- ZTEMP_BUD(:,:,:) = prsvs(:, :, :, ji) * prhodj(:, :, :)
- !$acc end kernels
- call Budget_store_init( tbudgets(NBUDGET_SV1 - 1 + ji), Trim( hbudname ), ZTEMP_BUD(:,:,:) )
- end do
- end if
- if ( lbudget_sv .and. hcloud == 'LIMA' ) then
- do ji = nsv_lima_beg, isv_lima_end
- !$acc kernels present(ZTEMP_BUD)
- ZTEMP_BUD(:,:,:) = prsvs(:, :, :, ji) * prhodj(:, :, :)
- !$acc end kernels
- call Budget_store_init( tbudgets(NBUDGET_SV1 - 1 + ji), Trim( hbudname ), ZTEMP_BUD(:,:,:) )
- end do
- end if
-end if
-
-
-
-!$acc data present( zt, zexn, zlv, zcph, zls, zcor )
-
-!$acc kernels present_cr(zexn,zt,zlv)
-#ifndef MNH_BITREP
-zexn(:, :, :) = ( ppabst(:, :, :) / xp00 ) ** (xrd / xcpd )
-#else
-zexn(:, :, :) = Br_pow( ppabst(:, :, :) / xp00, xrd / xcpd )
-#endif
-zt (:, :, :) = ptht(:, :, :) * zexn(:, :, :)
-zlv (:, :, :) = xlvtt + ( xcpv - xcl ) * ( zt(:, :, :) - xtt )
-!$acc end kernels
-if ( hcloud == 'ICE3' .or. hcloud == 'ICE4' .or. hcloud == 'LIMA' ) then
-!$acc kernels present_cr(zls)
- zls(:, :, :) = xlstt + ( xcpv - xci ) * ( zt(:, :, :) - xtt )
-!$acc end kernels
-end if
-!$acc kernels
-zcph(:, :, :) = xcpd + xcpv * prt(:, :, :, 1)
-!$acc end kernels
-
-#ifndef MNH_OPENACC
-deallocate( zt )
-#endif
-
-CLOUD: select case ( hcloud )
- case ( 'KESS' )
- jrmax = Size( prrs, 4 )
- do jr = 2, jrmax
-!PW: kernels directive inside do loop on jr because compiler bug... (NVHPC 21.7)
-!$acc kernels
- where ( prrs(:, :, :, jr) < 0. )
- prrs(:, :, :, 1) = prrs(:, :, :, 1) + prrs(:, :, :, jr)
- prths(:, :, :) = prths(:, :, :) - prrs(:, :, :, jr) * zlv(:, :, :) / &
- ( zcph(:, :, :) * zexn(:, :, :) )
- prrs(:, :, :, jr) = 0.
- end where
-!$acc end kernels
- end do
-
-!$acc kernels
- where ( prrs(:, :, :, 1) < 0. .and. prrs(:, :, :, 2) > 0. )
- prrs(:, :, :, 1) = prrs(:, :, :, 1) + prrs(:, :, :, 2)
- prths(:, :, :) = prths(:, :, :) - prrs(:, :, :, 2) * zlv(:, :, :) / &
- ( zcph(:, :, :) * zexn(:, :, :) )
- prrs(:, :, :, 2) = 0.
- end where
-!$acc end kernels
-
-
- case( 'ICE3', 'ICE4' )
- if ( hbudname == 'NETUR' ) then
- jrmax = 4
- else
- jrmax = Size( prrs, 4 )
- end if
-!$acc kernels
- do jr = 4, jrmax
- do concurrent( ji = 1 : jiu, jj = 1 : jju, jk = 1 : jku )
- if ( prrs(ji, jj, jk, jr) < 0. ) then
- prrs(ji, jj, jk, 1) = prrs(ji, jj, jk, 1) + prrs(ji, jj, jk, jr)
- prths(ji, jj, jk) = prths(ji, jj, jk) - prrs(ji, jj, jk, jr) * zls(ji, jj, jk) / &
- ( zcph(ji, jj, jk) * zexn(ji, jj, jk) )
- prrs(ji, jj, jk, jr) = 0.
- end if
- end do
- end do
-!$acc end kernels
-!
-! cloud
- if ( hbudname == 'NETUR' ) then
- jrmax = 2
- else
- jrmax = 3
- end if
-!$acc kernels
- do jr = 2, jrmax
- do concurrent( ji = 1 : jiu, jj = 1 : jju, jk = 1 : jku )
- if ( prrs(ji, jj, jk, jr) < 0. ) then
- prrs(ji, jj, jk, 1) = prrs(ji, jj, jk, 1) + prrs(ji, jj, jk, jr)
- prths(ji, jj, jk) = prths(ji, jj, jk) - prrs(ji, jj, jk, jr) * zlv(ji, jj, jk) / &
- ( zcph(ji, jj, jk) * zexn(ji, jj, jk) )
- prrs(ji, jj, jk, jr) = 0.
- end if
- end do
- end do
-!
-! if rc or ri are positive, we can correct negative rv
-! cloud
- do concurrent( ji = 1 : jiu, jj = 1 : jju, jk = 1 : jku )
- if ( prrs(ji, jj, jk, 1) < 0. .and. prrs(ji, jj, jk, 2) > 0. ) then
- prrs(ji, jj, jk, 1) = prrs(ji, jj, jk, 1) + prrs(ji, jj, jk, 2)
- prths(ji, jj, jk) = prths(ji, jj, jk) - prrs(ji, jj, jk, 2) * zlv(ji, jj, jk) / &
- ( zcph(ji, jj, jk) * zexn(ji, jj, jk) )
- prrs(ji, jj, jk, 2) = 0.
- end if
- end do
-! ice
- if ( krr > 3 ) then
-#ifdef MNH_COMPILER_NVHPC
-!$acc loop independent collapse(3)
-#endif
- do concurrent( ji = 1 : jiu, jj = 1 : jju, jk = 1 : jku )
- if ( prrs(ji, jj, jk, 1) < 0. .and. prrs(ji, jj, jk, 4) > 0. ) then
- zcor(ji, jj, jk) = Min( -prrs(ji, jj, jk, 1), prrs(ji, jj, jk, 4) )
- prrs(ji, jj, jk, 1) = prrs(ji, jj, jk, 1) + zcor(ji, jj, jk)
- prths(ji, jj, jk) = prths(ji, jj, jk) - zcor(ji, jj, jk) * zls(ji, jj, jk) / &
- ( zcph(ji, jj, jk) * zexn(ji, jj, jk) )
- prrs(ji, jj, jk, 4) = prrs(ji, jj, jk, 4) - zcor(ji, jj, jk)
- end if
- end do
- end if
-!$acc end kernels
-!
-!
- case( 'C2R2', 'KHKO' )
-!$acc kernels
- where ( prrs(:, :, :, 2) < 0. .or. prsvs(:, :, :, nsv_c2r2beg + 1) < 0. )
- prsvs(:, :, :, nsv_c2r2beg) = 0.
- end where
-!$acc end kernels
- do jsv = 2, 3
-!$acc kernels
-!PW: kernels directive inside do loop on jr because compiler bug... (NVHPC 21.7)
- where ( prrs(:, :, :, jsv) < 0. .or. prsvs(:, :, :, nsv_c2r2beg - 1 + jsv) < 0. )
- prrs(:, :, :, 1) = prrs(:, :, :, 1) + prrs(:, :, :, jsv)
- prths(:, :, :) = prths(:, :, :) - prrs(:, :, :, jsv) * zlv(:, :, :) / &
- ( zcph(:, :, :) * zexn(:, :, :) )
- prrs(:, :, :, jsv) = 0.
- prsvs(:, :, :, nsv_c2r2beg - 1 + jsv) = 0.
- end where
-!$acc end kernels
- end do
-!$acc kernels
- where ( prrs(:, :, :, 1) < 0. .and. prrs(:, :, :, 2) > 0. )
- prrs(:, :, :, 1) = prrs(:, :, :, 1) + prrs(:, :, :, 2)
- prths(:, :, :) = prths(:, :, :) - prrs(:, :, :, 2) * zlv(:, :, :) / &
- ( zcph(:, :, :) * zexn(:, :, :) )
- prrs(:, :, :, 2) = 0.
- prsvs(:, :, :, nsv_c2r2beg + 1) = 0.
- end where
-!$acc end kernels
-!
-!
- case( 'LIMA' )
-!$acc kernels
-! Correction where rc<0 or Nc<0
- if ( lwarm_lima ) then
- where ( prrs(:, :, :, 2) < xrtmin_lima(2) / ptstep .or. prsvs(:, :, :, nsv_lima_nc) < xctmin_lima(2) / ptstep )
- prrs(:, :, :, 1) = prrs(:, :, :, 1) + prrs(:, :, :, 2)
- prths(:, :, :) = prths(:, :, :) - prrs(:, :, :, 2) * zlv(:, :, :) / &
- ( zcph(:, :, :) * zexn(:, :, :) )
- prrs(:, :, :, 2) = 0.
- prsvs(:, :, :, nsv_lima_nc) = 0.
- end where
- where ( prrs(:, :, :, 1) < 0. .and. prrs(:, :, :, 2) > 0. )
- prrs(:, :, :, 1) = prrs(:, :, :, 1) + prrs(:, :, :, 2)
- prths(:, :, :) = prths(:, :, :) - prrs(:, :, :, 2) * zlv(:, :, :) / &
- ( zcph(:, :, :) * zexn(:, :, :) )
- prrs(:, :, :, 2) = 0.
- prsvs(:, :, :, nsv_lima_nc) = 0.
- end where
- end if
-! Correction where rr<0 or Nr<0
- if ( lwarm_lima .and. lrain_lima ) then
- where ( prrs(:, :, :, 3) < xrtmin_lima(3) / ptstep .or. prsvs(:, :, :, nsv_lima_nr) < xctmin_lima(3) / ptstep )
- prrs(:, :, :, 1) = prrs(:, :, :, 1) + prrs(:, :, :, 3)
- prths(:, :, :) = prths(:, :, :) - prrs(:, :, :, 3) * zlv(:, :, :) / &
- ( zcph(:, :, :) * zexn(:, :, :) )
- prrs(:, :, :, 3) = 0.
- prsvs(:, :, :, nsv_lima_nr) = 0.
- end where
- end if
-!$acc end kernels
-! Correction where ri<0 or Ni<0
- if ( lcold_lima ) then
-!$acc kernels
- where ( prrs(:, :, :, 4) < xrtmin_lima(4) / ptstep .or. prsvs(:, :, :, nsv_lima_ni) < xctmin_lima(4) / ptstep )
- prrs(:, :, :, 1) = prrs(:, :, :, 1) + prrs(:, :, :, 4)
- prths(:, :, :) = prths(:, :, :) - prrs(:, :, :, 4) * zls(:, :, :) / &
- ( zcph(:, :, :) * zexn(:, :, :) )
- prrs(:, :, :, 4) = 0.
- prsvs(:, :, :, nsv_lima_ni) = 0.
- end where
-!$acc end kernels
- if ( hbudname /= 'NETUR' ) then
- do jr = 5, Size( prrs, 4 )
-!PW: kernels directive inside do loop on jr because compiler bug... (NVHPC 21.7)
-!$acc kernels
- where ( prrs(:, :, :, jr) < 0. )
- prrs(:, :, :, 1) = prrs(:, :, :, 1) + prrs(:, :, :, jr)
- prths(:, :, :) = prths(:, :, :) - prrs(:, :, :, jr) * zls(:, :, :) / &
- ( zcph(:, :, :) * zexn(:, :, :) )
- prrs(:, :, :, jr) = 0.
- end where
-!$acc end kernels
- end do
- end if
- if(krr > 3) then
-!$acc kernels
- where ( prrs(:, :, :, 1) < 0. .and. prrs(:, :, :, 4) > 0. )
- zcor(:, :, :) = Min( -prrs(:, :, :, 1), prrs(:, :, :, 4) )
- prrs(:, :, :, 1) = prrs(:, :, :, 1) + zcor(:, :, :)
- prths(:, :, :) = prths(:, :, :) - zcor(:, :, :) * zls(:, :, :) / &
- ( zcph(:, :, :) * zexn(:, :, :) )
- prrs(:, :, :, 4) = prrs(:, :, :, 4) - zcor(:, :, :)
- end where
-!$acc end kernels
- end if
- end if
-
-!$acc kernels
- prsvs(:, :, :, nsv_lima_beg : isv_lima_end) = Max( 0.0, prsvs(:, :, :, nsv_lima_beg : isv_lima_end) )
-!$acc end kernels
-
-end select CLOUD
-
-!$acc end data
-
-#ifndef MNH_OPENACC
-deallocate( zexn, zlv, zcph, zls, zcor , ZTEMP_BUD)
-#else
-!Release all memory allocated with Mnh_mem_get calls since last call to Mnh_mem_position_pin
-call Mnh_mem_release()
-#endif
-
-if ( hbudname /= 'NECON' .and. hbudname /= 'NEGA' ) then
- if ( hcloud == 'KESS' .or. hcloud == 'ICE3' .or. hcloud == 'ICE4' .or. &
- hcloud == 'KHKO' .or. hcloud == 'C2R2' .or. hcloud == 'LIMA' ) then
- if ( lbudget_th ) call Budget_store_end( tbudgets(NBUDGET_TH), Trim( hbudname ), prths(:, :, :) )
- if ( lbudget_rv ) call Budget_store_end( tbudgets(NBUDGET_RV), Trim( hbudname ), prrs (:, :, :, 1) )
- if ( lbudget_rc ) call Budget_store_end( tbudgets(NBUDGET_RC), Trim( hbudname ), prrs (:, :, :, 2) )
- if ( lbudget_rr .and. &
- ( hbudname /= 'NETUR' .or. &
- ( hbudname == 'NETUR' .and. ( hcloud == 'C2R2' .or. hcloud == 'KHKO' .or. hcloud == 'LIMA' ) ) ) ) &
- call Budget_store_end( tbudgets(NBUDGET_RR), Trim( hbudname ), prrs (:, :, :, 3) )
- IF (lbudget_ri .and. &
- ( hbudname /= 'NETUR' .or. &
- ( hbudname == 'NETUR' .and. ( hcloud == 'ICE3' .or. hcloud == 'ICE4' .or. hcloud == 'LIMA' ) ) ) ) &
- call Budget_store_end( tbudgets(NBUDGET_RI), Trim( hbudname ), prrs (:, :, :, 4) )
- if ( lbudget_rs .and. hbudname /= 'NETUR' ) call Budget_store_end( tbudgets(NBUDGET_RS), Trim( hbudname ), prrs (:, :, :, 5) )
- if ( lbudget_rg .and. hbudname /= 'NETUR' ) call Budget_store_end( tbudgets(NBUDGET_RG), Trim( hbudname ), prrs (:, :, :, 6) )
- if ( lbudget_rh .and. hbudname /= 'NETUR' ) call Budget_store_end( tbudgets(NBUDGET_RH), Trim( hbudname ), prrs (:, :, :, 7) )
- end if
-
- if ( lbudget_sv .and. ( hcloud == 'C2R2' .or. hcloud == 'KHKO' ) ) then
- do ji = nsv_c2r2beg, nsv_c2r2end
- call Budget_store_end( tbudgets(NBUDGET_SV1 - 1 + ji), Trim( hbudname ), prsvs(:, :, :, ji) )
- end do
- end if
- if ( lbudget_sv .and. hcloud == 'LIMA' ) then
- do ji = nsv_lima_beg, isv_lima_end
- call Budget_store_end( tbudgets(NBUDGET_SV1 - 1 + ji), Trim( hbudname ), prsvs(:, :, :, ji) )
- end do
- end if
-else !NECON + NEGA
- if ( hcloud == 'KESS' .or. hcloud == 'ICE3' .or. hcloud == 'ICE4' .or. &
- hcloud == 'KHKO' .or. hcloud == 'C2R2' .or. hcloud == 'LIMA' ) then
- if ( lbudget_th) then
- !$acc kernels present(ZTEMP_BUD)
- ZTEMP_BUD(:,:,:) = prths(:, :, :) * prhodj(:, :, :)
- !$acc end kernels
- call Budget_store_end( tbudgets(NBUDGET_TH), Trim( hbudname ), ZTEMP_BUD(:,:,:) )
- end if
- if ( lbudget_rv) then
- !$acc kernels present(ZTEMP_BUD)
- ZTEMP_BUD(:,:,:) = prrs (:, :, :, 1) * prhodj(:, :, :)
- !$acc end kernels
- call Budget_store_end( tbudgets(NBUDGET_RV), Trim( hbudname ), ZTEMP_BUD(:,:,:) )
- end if
- if ( lbudget_rc) then
- !$acc kernels present(ZTEMP_BUD)
- ZTEMP_BUD(:,:,:) = prrs (:, :, :, 2) * prhodj(:, :, :)
- !$acc end kernels
- call Budget_store_end( tbudgets(NBUDGET_RC), Trim( hbudname ), ZTEMP_BUD(:,:,:) )
- end if
- if ( lbudget_rr) then
- !$acc kernels present(ZTEMP_BUD)
- ZTEMP_BUD(:,:,:) = prrs (:, :, :, 3) * prhodj(:, :, :)
- !$acc end kernels
- call Budget_store_end( tbudgets(NBUDGET_RR), Trim( hbudname ), ZTEMP_BUD(:,:,:) )
- end if
- if ( lbudget_ri) then
- !$acc kernels present(ZTEMP_BUD)
- ZTEMP_BUD(:,:,:) = prrs (:, :, :, 4) * prhodj(:, :, :)
- !$acc end kernels
- call Budget_store_end( tbudgets(NBUDGET_RI), Trim( hbudname ), ZTEMP_BUD(:,:,:) )
- end if
- if ( lbudget_rs) then
- !$acc kernels present(ZTEMP_BUD)
- ZTEMP_BUD(:,:,:) = prrs (:, :, :, 5) * prhodj(:, :, :)
- !$acc end kernels
- call Budget_store_end( tbudgets(NBUDGET_RS), Trim( hbudname ), ZTEMP_BUD(:,:,:) )
- end if
- if ( lbudget_rg) then
- !$acc kernels present(ZTEMP_BUD)
- ZTEMP_BUD(:,:,:) = prrs (:, :, :, 6) * prhodj(:, :, :)
- !$acc end kernels
- call Budget_store_end( tbudgets(NBUDGET_RG), Trim( hbudname ), ZTEMP_BUD(:,:,:) )
- end if
- if ( lbudget_rh) then
- !$acc kernels present(ZTEMP_BUD)
- ZTEMP_BUD(:,:,:) = prrs (:, :, :, 7) * prhodj(:, :, :)
- !$acc end kernels
- call Budget_store_end( tbudgets(NBUDGET_RH), Trim( hbudname ), ZTEMP_BUD(:,:,:) )
- end if
- end if
-
- if ( lbudget_sv .and. ( hcloud == 'C2R2' .or. hcloud == 'KHKO' ) ) then
- do ji = nsv_c2r2beg, nsv_c2r2end
- !$acc kernels present(ZTEMP_BUD)
- ZTEMP_BUD(:,:,:) = prsvs(:, :, :, ji) * prhodj(:, :, :)
- !$acc end kernels
- call Budget_store_end( tbudgets(NBUDGET_SV1 - 1 + ji), Trim( hbudname ), ZTEMP_BUD(:,:,:) )
- end do
- end if
- if ( lbudget_sv .and. hcloud == 'LIMA' ) then
- do ji = nsv_lima_beg, isv_lima_end
- !$acc kernels present(ZTEMP_BUD)
- ZTEMP_BUD(:,:,:) = prsvs(:, :, :, ji) * prhodj(:, :, :)
- !$acc end kernels
- call Budget_store_end( tbudgets(NBUDGET_SV1 - 1 + ji), Trim( hbudname ), ZTEMP_BUD(:,:,:) )
- end do
- end if
-end if
-
-if ( mppdb_initialized ) then
- !Check all INOUT arrays
- call Mppdb_check( prths, "Sources_neg_correct end:prths")
- call Mppdb_check( prrs, "Sources_neg_correct end:prrs")
- call Mppdb_check( prsvs, "Sources_neg_correct end:prsvs")
-end if
-
-!$acc end data
-
-end subroutine Sources_neg_correct
-
-end module mode_sources_neg_correct
diff --git a/src/ZSOLVER/tridiag_thermo.f90 b/src/ZSOLVER/tridiag_thermo.f90
deleted file mode 100644
index c96d99382707f1623ae9460b973256a7855ce63b..0000000000000000000000000000000000000000
--- a/src/ZSOLVER/tridiag_thermo.f90
+++ /dev/null
@@ -1,477 +0,0 @@
-!MNH_LIC Copyright 2003-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_TRIDIAG_THERMO
-! ###################
-!
-INTERFACE
-!
- SUBROUTINE TRIDIAG_THERMO(KKA,KKU,KKL,PVARM,PF,PDFDDTDZ,PTSTEP,PIMPL, &
- PDZZ,PRHODJ,PVARP )
-!
-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=AR
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PVARM ! variable at t-1 at mass point
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PF ! flux in dT/dt=-dF/dz at flux point
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PDFDDTDZ! dF/d(dT/dz) at flux point
-REAL, INTENT(IN) :: PTSTEP ! Double time step
-REAL, INTENT(IN) :: PIMPL ! implicit weight
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ ! Dz at flux point
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODJ ! (dry rho)*J at mass point
-!
-REAL, DIMENSION(:,:,:), INTENT(OUT):: PVARP ! variable at t+1 at mass point
-!
-END SUBROUTINE TRIDIAG_THERMO
-!
-END INTERFACE
-!
-END MODULE MODI_TRIDIAG_THERMO
-!
-! #################################################
- SUBROUTINE TRIDIAG_THERMO(KKA,KKU,KKL,PVARM,PF,PDFDDTDZ,PTSTEP,PIMPL, &
- PDZZ,PRHODJ,PVARP )
-! #################################################
-!
-!
-!!**** *TRIDIAG_THERMO* - routine to solve a time implicit scheme
-!!
-!!
-!! PURPOSE
-!! -------
-! The purpose of this routine is to give a field PVARP at t+1, by
-! solving an implicit TRIDIAGonal system obtained by the
-! discretization of the vertical turbulent diffusion. It should be noted
-! that the degree of implicitness can be varied (PIMPL parameter) and that
-! the function of F(dT/dz) must have been linearized.
-! PVARP is localized at a mass point.
-!
-!!** METHOD
-!! ------
-!!
-!! [T(+) - T(-)]/2Dt = -d{ F + dF/d(dT/dz) * [impl*dT/dz(+) + expl* dT/dz(-)] }/dz
-!!
-!! It is discretized as follows:
-!!
-!! PRHODJ(k)*PVARP(k)/PTSTEP
-!! =
-!! PRHODJ(k)*PVARM(k)/PTSTEP
-!! - (PRHODJ(k+1)+PRHODJ(k) )/2. * PF(k+1)/PDZZ(k+1)
-!! + (PRHODJ(k) +PRHODJ(k-1))/2. * PF(k) /PDZZ(k)
-!! + (PRHODJ(k+1)+PRHODJ(k) )/2. * ZEXPL* PDFDDTDZ(k+1) * PVARM(k+1)/PDZZ(k+1)**2
-!! - (PRHODJ(k+1)+PRHODJ(k) )/2. * PIMPL* PDFDDTDZ(k+1) * PVARP(k+1)/PDZZ(k+1)**2
-!! - (PRHODJ(k+1)+PRHODJ(k) )/2. * ZEXPL* PDFDDTDZ(k+1) * PVARM(k) /PDZZ(k+1)**2
-!! + (PRHODJ(k+1)+PRHODJ(k) )/2. * PIMPL* PDFDDTDZ(k+1) * PVARP(k) /PDZZ(k+1)**2
-!! - (PRHODJ(k) +PRHODJ(k-1))/2. * ZEXPL* PDFDDTDZ(k) * PVARM(k) /PDZZ(k)**2
-!! + (PRHODJ(k) +PRHODJ(k-1))/2. * PIMPL* PDFDDTDZ(k) * PVARP(k) /PDZZ(k)**2
-!! + (PRHODJ(k) +PRHODJ(k-1))/2. * ZEXPL* PDFDDTDZ(k) * PVARM(k-1)/PDZZ(k)**2
-!! - (PRHODJ(k) +PRHODJ(k-1))/2. * PIMPL* PDFDDTDZ(k) * PVARP(k-1)/PDZZ(k)**2
-!!
-!!
-!! The system to solve is:
-!!
-!! A*PVARP(k-1) + B*PVARP(k) + C*PVARP(k+1) = Y(k)
-!!
-!!
-!! The RHS of the linear system in PVARP writes:
-!!
-!! y(k) = PRHODJ(k)*PVARM(k)/PTSTEP
-!! - (PRHODJ(k+1)+PRHODJ(k) )/2. * PF(k+1)/PDZZ(k+1)
-!! + (PRHODJ(k) +PRHODJ(k-1))/2. * PF(k) /PDZZ(k)
-!! + (PRHODJ(k+1)+PRHODJ(k) )/2. * ZEXPL* PDFDDTDZ(k+1) * PVARM(k+1)/PDZZ(k+1)**2
-!! - (PRHODJ(k+1)+PRHODJ(k) )/2. * ZEXPL* PDFDDTDZ(k+1) * PVARM(k) /PDZZ(k+1)**2
-!! - (PRHODJ(k) +PRHODJ(k-1))/2. * ZEXPL* PDFDDTDZ(k) * PVARM(k) /PDZZ(k)**2
-!! + (PRHODJ(k) +PRHODJ(k-1))/2. * ZEXPL* PDFDDTDZ(k) * PVARM(k-1)/PDZZ(k)**2
-!!
-!!
-!! Then, the classical TRIDIAGonal algorithm is used to invert the
-!! implicit operator. Its matrix is given by:
-!!
-!! ( b(ikb) c(ikb) 0 0 0 0 0 0 )
-!! ( 0 a(ikb+1) b(ikb+1) c(ikb+1) 0 ... 0 0 0 )
-!! ( 0 0 a(ikb+2) b(ikb+2) c(ikb+2). 0 0 0 )
-!! .......................................................................
-!! ( 0 ... 0 a(k) b(k) c(k) 0 ... 0 0 )
-!! .......................................................................
-!! ( 0 0 0 0 0 ...a(ike-1) b(ike-1) c(ike-1))
-!! ( 0 0 0 0 0 ... 0 a(ike) b(ike) )
-!!
-!! ikb and ike represent the first and the last inner mass levels of the
-!! model. The coefficients are:
-!!
-!! a(k) = + (PRHODJ(k) +PRHODJ(k-1))/2. * PIMPL* PDFDDTDZ(k) /PDZZ(k)**2
-!! b(k) = PRHODJ(k) / PTSTEP
-!! - (PRHODJ(k+1)+PRHODJ(k) )/2. * PIMPL* PDFDDTDZ(k+1)/PDZZ(k+1)**2
-!! - (PRHODJ(k) +PRHODJ(k-1))/2. * PIMPL* PDFDDTDZ(k) /PDZZ(k)**2
-!! c(k) = + (PRHODJ(k+1)+PRHODJ(k) )/2. * PIMPL* PDFDDTDZ(k+1)/PDZZ(k+1)**2
-!!
-!! for all k /= ikb or ike
-!!
-!!
-!! b(ikb) = PRHODJ(ikb) / PTSTEP
-!! -(PRHODJ(ikb+1)+PRHODJ(ikb))/2.*PIMPL*PDFDDTDZ(ikb+1)/PDZZ(ikb+1)**2
-!! c(ikb) = +(PRHODJ(ikb+1)+PRHODJ(ikb))/2.*PIMPL*PDFDDTDZ(ikb+1)/PDZZ(ikb+1)**2
-!!
-!! b(ike) = PRHODJ(ike) / PTSTEP
-!! -(PRHODJ(ike)+PRHODJ(ike-1))/2.*PIMPL*PDFDDTDZ(ike)/PDZZ(ike)**2
-!! a(ike) = +(PRHODJ(ike)+PRHODJ(ike-1))/2.*PIMPL*PDFDDTDZ(ike)/PDZZ(ike)**2
-!!
-!!
-!! EXTERNAL
-!! --------
-!!
-!! NONE
-!!
-!! IMPLICIT ARGUMENTS
-!! ------------------
-!!
-!! REFERENCE
-!! ---------
-!! Press et al: Numerical recipes (1986) Cambridge Univ. Press
-!!
-!! AUTHOR
-!! ------
-!! V. Masson * Meteo-France *
-!!
-!! MODIFICATIONS
-!! -------------
-!! Original 04/2003 (from tridiag.f90)
-!! M.Moge 04/2016 Use openACC directives to port the TURB part of Meso-NH on GPU
-!! ---------------------------------------------------------------------
-!
-!* 0. DECLARATIONS
-!
-USE MODD_PARAMETERS, ONLY : JPVEXT_TURB
-
-use mode_mppdb
-
-#ifndef MNH_OPENACC
-USE MODI_SHUMAN
-#else
-USE MODI_SHUMAN_DEVICE
-#endif
-#ifdef MNH_BITREP
-USE MODI_BITREP
-#endif
-!
-#ifdef MNH_OPENACC
-USE MODE_MNH_ZWORK, ONLY: MNH_MEM_GET, MNH_MEM_POSITION_PIN, MNH_MEM_RELEASE
-#endif
-!
-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
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PVARM ! variable at t-1 at mass point
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PF ! flux in dT/dt=-dF/dz at flux point
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PDFDDTDZ! dF/d(dT/dz) at flux point
-REAL, INTENT(IN) :: PTSTEP ! Double time step
-REAL, INTENT(IN) :: PIMPL ! implicit weight
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ ! Dz at flux point
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODJ ! (dry rho)*J at mass point
-!
-REAL, DIMENSION(:,:,:), INTENT(OUT):: PVARP ! variable at t+1 at mass point
-!
-!* 0.2 declarations of local variables
-!
-REAL, DIMENSION(:,:,:), pointer , contiguous :: ZRHODJ_DFDDTDZ_O_DZ2
-REAL, DIMENSION(:,:,:), pointer , contiguous :: ZMZM_RHODJ
-REAL, DIMENSION(:,:,:), pointer , contiguous :: ZA, ZB, ZC
-REAL, DIMENSION(:,:,:), pointer , contiguous :: ZY ,ZGAM
- ! RHS of the equation, 3D work array
-REAL, DIMENSION(:,:), pointer , contiguous :: ZBET
- ! 2D work array
-INTEGER :: JI,JJ,JK ! loop counter
-INTEGER :: IKB,IKE ! inner vertical limits
-INTEGER :: IKT ! array size in k direction
-INTEGER :: IKTB,IKTE ! start, end of k loops in physical domain
-!
-!
-#ifdef MNH_OPENACC
-REAL, DIMENSION(:,:,:), pointer , contiguous :: ZTMP1_DEVICE
-#endif
-
-INTEGER :: JIU,JJU,JKU
-
-! ---------------------------------------------------------------------------
-
-!$acc data present( PVARM, PF, PDFDDTDZ, PDZZ, PRHODJ, PVARP )
-
-if ( mppdb_initialized ) then
- !Check all in arrays
- call Mppdb_check( pvarm, "Tridiag_thermo beg:pvarm" )
- call Mppdb_check( pf, "Tridiag_thermo beg:pf" )
- call Mppdb_check( pdfddtdz, "Tridiag_thermo beg:pdfddtdz" )
- call Mppdb_check( pdzz, "Tridiag_thermo beg:pdzz" )
- call Mppdb_check( prhodj, "Tridiag_thermo beg:prhodj" )
-end if
-
-JIU = size( pvarm, 1 )
-JJU = size( pvarm, 2 )
-JKU = size( pvarm, 3 )
-
-#ifndef MNH_OPENACC
-allocate( zrhodj_dfddtdz_o_dz2(JIU,JJU,JKU ) )
-allocate( zmzm_rhodj (JIU,JJU,JKU ) )
-allocate( za (JIU,JJU,JKU ) )
-allocate( zb (JIU,JJU,JKU ) )
-allocate( zc (JIU,JJU,JKU ) )
-allocate( zy (JIU,JJU,JKU ) )
-allocate( zgam (JIU,JJU,JKU ) )
-allocate( zbet (JIU,JJU ) )
-#else
-!Pin positions in the pools of MNH memory
-CALL MNH_MEM_POSITION_PIN()
-
-CALL MNH_MEM_GET( zrhodj_dfddtdz_o_dz2, JIU, JJU, JKU )
-CALL MNH_MEM_GET( zmzm_rhodj , JIU, JJU, JKU )
-CALL MNH_MEM_GET( za , JIU, JJU, JKU )
-CALL MNH_MEM_GET( zb , JIU, JJU, JKU )
-CALL MNH_MEM_GET( zc , JIU, JJU, JKU )
-CALL MNH_MEM_GET( zy , JIU, JJU, JKU )
-CALL MNH_MEM_GET( zgam , JIU, JJU, JKU )
-CALL MNH_MEM_GET( zbet , JIU, JJU )
-
-CALL MNH_MEM_GET( ztmp1_device, JIU, JJU, JKU )
-#endif
-
-!$acc data present( zrhodj_dfddtdz_o_dz2, zmzm_rhodj, za, zb, zc, zy, zgam, zbet, ztmp1_device )
-!
-!* 1. Preliminaries
-! -------------
-!
-IKT=SIZE(PVARM,3)
-IKTB=1+JPVEXT_TURB
-IKTE=IKT-JPVEXT_TURB
-IKB=KKA+JPVEXT_TURB*KKL
-IKE=KKU-JPVEXT_TURB*KKL
-!
-#ifndef MNH_OPENACC
-ZMZM_RHODJ = MZM(PRHODJ)
-#else
-CALL MZM_DEVICE(PRHODJ,ZMZM_RHODJ)
-#endif
-!$acc kernels present_cr(ZRHODJ_DFDDTDZ_O_DZ2) ! async
-#ifndef MNH_BITREP
-ZRHODJ_DFDDTDZ_O_DZ2(:,:,:) = ZMZM_RHODJ(:,:,:)*PDFDDTDZ(:,:,:)/PDZZ(:,:,:)**2
-#else
-!$acc_nv loop independent collapse(3)
-DO CONCURRENT ( JI=1:JIU,JJ=1:JJU,JK=1:JKU)
- ZRHODJ_DFDDTDZ_O_DZ2(JI,JJ,JK) = ZMZM_RHODJ(JI,JJ,JK)*PDFDDTDZ(JI,JJ,JK)/BR_P2(PDZZ(JI,JJ,JK))
-END DO !CONCURRENT
-#endif
-!$acc end kernels
-!
-!$acc kernels ! async
-ZA=0.
-ZB=0.
-ZC=0.
-ZY=0.
-!$acc end kernels
-! acc wait
-!
-!
-!* 2. COMPUTE THE RIGHT HAND SIDE
-! ---------------------------
-!
-!$acc kernels ! async
-#ifdef MNH_COMPILER_NVHPC
-!$acc loop independent collapse(2)
-#endif
-DO CONCURRENT ( JI=1:JIU,JJ=1:JJU)
-ZY(JI,JJ,IKB) = PRHODJ(JI,JJ,IKB)*PVARM(JI,JJ,IKB)/PTSTEP &
- - ZMZM_RHODJ(JI,JJ,IKB+KKL) * PF(JI,JJ,IKB+KKL)/PDZZ(JI,JJ,IKB+KKL) &
- + ZMZM_RHODJ(JI,JJ,IKB ) * PF(JI,JJ,IKB )/PDZZ(JI,JJ,IKB ) &
- + ZRHODJ_DFDDTDZ_O_DZ2(JI,JJ,IKB+KKL) * PIMPL * PVARM(JI,JJ,IKB+KKL) &
- - ZRHODJ_DFDDTDZ_O_DZ2(JI,JJ,IKB+KKL) * PIMPL * PVARM(JI,JJ,IKB )
-END DO !CONCURRENT
-!$acc end kernels
-!
-!$acc kernels ! async
-#ifdef MNH_COMPILER_NVHPC
-!$acc loop independent collapse(3)
-#endif
-DO CONCURRENT ( JI=1:JIU,JJ=1:JJU,JK=IKTB+1:IKTE-1)
- ZY(JI,JJ,JK) = PRHODJ(JI,JJ,JK)*PVARM(JI,JJ,JK)/PTSTEP &
- - ZMZM_RHODJ(JI,JJ,JK+KKL) * PF(JI,JJ,JK+KKL)/PDZZ(JI,JJ,JK+KKL) &
- + ZMZM_RHODJ(JI,JJ,JK ) * PF(JI,JJ,JK )/PDZZ(JI,JJ,JK ) &
- + ZRHODJ_DFDDTDZ_O_DZ2(JI,JJ,JK+KKL) * PIMPL * PVARM(JI,JJ,JK+KKL) &
- - ZRHODJ_DFDDTDZ_O_DZ2(JI,JJ,JK+KKL) * PIMPL * PVARM(JI,JJ,JK ) &
- - ZRHODJ_DFDDTDZ_O_DZ2(JI,JJ,JK ) * PIMPL * PVARM(JI,JJ,JK ) &
- + ZRHODJ_DFDDTDZ_O_DZ2(JI,JJ,JK ) * PIMPL * PVARM(JI,JJ,JK-KKL)
-END DO !CONCURRENT
-!$acc end kernels
-!
-!$acc kernels ! async
-#ifdef MNH_COMPILER_NVHPC
-!$acc loop independent collapse(2)
-#endif
-DO CONCURRENT ( JI=1:JIU,JJ=1:JJU)
-ZY(JI,JJ,IKE) = PRHODJ(JI,JJ,IKE)*PVARM(JI,JJ,IKE)/PTSTEP &
- - ZMZM_RHODJ(JI,JJ,IKE+KKL) * PF(JI,JJ,IKE+KKL)/PDZZ(JI,JJ,IKE+KKL) &
- + ZMZM_RHODJ(JI,JJ,IKE ) * PF(JI,JJ,IKE )/PDZZ(JI,JJ,IKE ) &
- - ZRHODJ_DFDDTDZ_O_DZ2(JI,JJ,IKE ) * PIMPL * PVARM(JI,JJ,IKE ) &
- + ZRHODJ_DFDDTDZ_O_DZ2(JI,JJ,IKE ) * PIMPL * PVARM(JI,JJ,IKE-KKL)
-END DO !CONCURRENT
-!$acc end kernels
-!
-! acc wait
-!
-!* 3. INVERSION OF THE TRIDIAGONAL SYSTEM
-! -----------------------------------
-!
-IF ( PIMPL > 1.E-10 ) THEN
-!
-!* 3.1 arrays A, B, C
-! --------------
-!
-!$acc kernels ! async
-#ifdef MNH_COMPILER_NVHPC
-!$acc loop independent collapse(2)
-#endif
-DO CONCURRENT ( JI=1:JIU,JJ=1:JJU)
- ZB(JI,JJ,IKB) = PRHODJ(JI,JJ,IKB)/PTSTEP &
- - ZRHODJ_DFDDTDZ_O_DZ2(JI,JJ,IKB+KKL) * PIMPL
-END DO !CONCURRENT
-!$acc end kernels
-!
-!$acc kernels ! async
-#ifdef MNH_COMPILER_NVHPC
-!$acc loop independent collapse(2)
-#endif
-DO CONCURRENT ( JI=1:JIU,JJ=1:JJU)
- ZC(JI,JJ,IKB) = ZRHODJ_DFDDTDZ_O_DZ2(JI,JJ,IKB+KKL) * PIMPL
-END DO !CONCURRENT
-!$acc end kernels
-!
-!$acc kernels ! async
-#ifdef MNH_COMPILER_NVHPC
-!$acc loop independent collapse(3)
-#endif
-DO CONCURRENT ( JI=1:JIU,JJ=1:JJU,JK=IKTB+1:IKTE-1)
- ZA(JI,JJ,JK) = ZRHODJ_DFDDTDZ_O_DZ2(JI,JJ,JK) * PIMPL
- ZB(JI,JJ,JK) = PRHODJ(JI,JJ,JK)/PTSTEP &
- - ZRHODJ_DFDDTDZ_O_DZ2(JI,JJ,JK+KKL) * PIMPL &
- - ZRHODJ_DFDDTDZ_O_DZ2(JI,JJ,JK) * PIMPL
- ZC(JI,JJ,JK) = ZRHODJ_DFDDTDZ_O_DZ2(JI,JJ,JK+KKL) * PIMPL
-END DO !CONCURRENT
-!$acc end kernels
-!
-!$acc kernels ! async
-#ifdef MNH_COMPILER_NVHPC
-!$acc loop independent collapse(2)
-#endif
-DO CONCURRENT ( JI=1:JIU,JJ=1:JJU)
- ZA(JI,JJ,IKE) = ZRHODJ_DFDDTDZ_O_DZ2(JI,JJ,IKE ) * PIMPL
- ZB(JI,JJ,IKE) = PRHODJ(JI,JJ,IKE)/PTSTEP &
- - ZRHODJ_DFDDTDZ_O_DZ2(JI,JJ,IKE ) * PIMPL
-END DO !CONCURRENT
-!$acc end kernels
-!
-! acc wait
-!
-!
-!* 3.2 going up
-! --------
-!
-!$acc kernels
-#ifdef MNH_COMPILER_NVHPC
-!$acc loop independent collapse(2)
-#endif
-DO CONCURRENT ( JI=1:JIU,JJ=1:JJU)
- ZBET(JI,JJ) = ZB(JI,JJ,IKB) ! bet = b(ikb)
- PVARP(JI,JJ,IKB) = ZY(JI,JJ,IKB) / ZBET(JI,JJ)
-END DO !CONCURRENT
-!
-!$acc loop seq
-DO JK = IKB+KKL,IKE-KKL,KKL
- ! gang+vector needed or parallisation vector only
- !$acc_nv loop independent gang, vector collapse(2)
- DO CONCURRENT ( JI=1:JIU,JJ=1:JJU)
- ZGAM(JI,JJ,JK) = ZC(JI,JJ,JK-KKL) / ZBET(JI,JJ)
- ! gam(k) = c(k-1) / bet
- ZBET(JI,JJ) = ZB(JI,JJ,JK) - ZA(JI,JJ,JK) * ZGAM(JI,JJ,JK)
- ! bet = b(k) - a(k)* gam(k)
- PVARP(JI,JJ,JK)= ( ZY(JI,JJ,JK) - ZA(JI,JJ,JK) * PVARP(JI,JJ,JK-KKL) ) / ZBET(JI,JJ)
- ! res(k) = (y(k) -a(k)*res(k-1))/ bet
- END DO !CONCURRENT
-END DO
-! special treatment for the last level
-#ifdef MNH_COMPILER_NVHPC
-!$acc loop independent collapse(2)
-#endif
-DO CONCURRENT ( JI=1:JIU,JJ=1:JJU)
- ZGAM(JI,JJ,IKE) = ZC(JI,JJ,IKE-KKL) / ZBET(JI,JJ)
- ! gam(k) = c(k-1) / bet
- ZBET(JI,JJ) = ZB(JI,JJ,IKE) - ZA(JI,JJ,IKE) * ZGAM(JI,JJ,IKE)
- ! bet = b(k) - a(k)* gam(k)
- PVARP(JI,JJ,IKE)= ( ZY(JI,JJ,IKE) - ZA(JI,JJ,IKE) * PVARP(JI,JJ,IKE-KKL) ) / ZBET(JI,JJ)
- ! res(k) = (y(k) -a(k)*res(k-1))/ bet
-END DO !CONCURRENT
-!
-!* 3.3 going down
-! ----------
-!
-!$acc loop seq
-DO JK = IKE-KKL,IKB,-1*KKL
- ! gang+vector needed or parallisation vector only
- !$acc_nv loop independent gang, vector collapse(2)
- DO CONCURRENT ( JI=1:JIU,JJ=1:JJU)
- PVARP(JI,JJ,JK) = PVARP(JI,JJ,JK) - ZGAM(JI,JJ,JK+KKL) * PVARP(JI,JJ,JK+KKL)
- END DO !CONCURRENT
-END DO
-!$acc end kernels
-!
-ELSE
-!
-!$acc kernels
-#ifdef MNH_COMPILER_NVHPC
-!$acc loop independent collapse(3)
-#endif
-DO CONCURRENT ( JI=1:JIU,JJ=1:JJU,JK=IKTB:IKTE)
- PVARP(JI,JJ,JK) = ZY(JI,JJ,JK) * PTSTEP / PRHODJ(JI,JJ,JK)
-END DO !CONCURRENT
-!$acc end kernels
-!
-END IF
-!
-!
-!* 4. FILL THE UPPER AND LOWER EXTERNAL VALUES
-! ----------------------------------------
-!
-!$acc kernels
-#ifdef MNH_COMPILER_NVHPC
-!$acc loop independent collapse(2)
-#endif
-DO CONCURRENT ( JI=1:JIU,JJ=1:JJU)
- PVARP(JI,JJ,KKA)=PVARP(JI,JJ,IKB)
- PVARP(JI,JJ,KKU)=PVARP(JI,JJ,IKE)
-END DO !CONCURRENT
-!$acc end kernels
-
-if ( mppdb_initialized ) then
- !Check all out arrays
- call Mppdb_check( pvarp, "Tridiag_thermo end:pvarp" )
-end if
-
-!$acc end data
-
-#ifndef MNH_OPENACC
-deallocate (zrhodj_dfddtdz_o_dz2,zmzm_rhodj,za,zb,zc,zy,zgam,zbet)
-#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 TRIDIAG_THERMO
diff --git a/src/ZSOLVER/turb_hor_thermo_flux.f90 b/src/ZSOLVER/turb_hor_thermo_flux.f90
deleted file mode 100644
index 52cff51886aedfc51234d2c944b6cf074f25af2b..0000000000000000000000000000000000000000
--- a/src/ZSOLVER/turb_hor_thermo_flux.f90
+++ /dev/null
@@ -1,1802 +0,0 @@
-!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_TURB_HOR_THERMO_FLUX
-! ################################
-!
-INTERFACE
-!
- SUBROUTINE TURB_HOR_THERMO_FLUX(KSPLT, KRR, KRRL, KRRI, &
- OTURB_FLX,OSUBG_COND, &
- TPFILE, &
- PK,PINV_PDXX,PINV_PDYY,PINV_PDZZ,PMZM_PRHODJ, &
- PDXX,PDYY,PDZZ,PDZX,PDZY, &
- PDIRCOSXW,PDIRCOSYW, &
- PRHODJ, &
- PSFTHM,PSFRM, &
- PWM,PTHLM,PRM, &
- PATHETA,PAMOIST,PSRCM,PFRAC_ICE, &
- PRTHLS,PRRS )
-!
-USE MODD_IO, ONLY: TFILEDATA
-!
-INTEGER, INTENT(IN) :: KSPLT ! split process index
-INTEGER, INTENT(IN) :: KRR ! number of moist var.
-INTEGER, INTENT(IN) :: KRRL ! number of liquid water var.
-INTEGER, INTENT(IN) :: KRRI ! number of ice water var.
-LOGICAL, INTENT(IN) :: OTURB_FLX ! switch to write the
- ! turbulent fluxes in the syncronous FM-file
-LOGICAL, INTENT(IN) :: OSUBG_COND ! Switch for sub-grid
-! condensation
-TYPE(TFILEDATA), INTENT(INOUT) :: TPFILE ! Output file
-!
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PK ! Turbulent diffusion doef.
- ! PK = PLM * SQRT(PTKEM)
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PINV_PDXX ! 1./PDXX
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PINV_PDYY ! 1./PDYY
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PINV_PDZZ ! 1./PDZZ
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PMZM_PRHODJ ! MZM(PRHODJ)
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX, PDYY, PDZZ, PDZX, PDZY
- ! Metric coefficients
-REAL, DIMENSION(:,:), INTENT(IN) :: PDIRCOSXW, PDIRCOSYW
-! Director Cosinus along x, y and z directions at surface w-point
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODJ ! density * grid volume
-!
-REAL, DIMENSION(:,:), INTENT(IN) :: PSFTHM,PSFRM
-!
-! Variables at t-1
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PWM
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHLM
-REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PRM ! mixing ratios at t-1,
- ! where PRM(:,:,:,1) = conservative mixing ratio
-!
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PATHETA ! coefficients between
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PAMOIST ! s and Thetal and Rnp
-
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PSRCM
- ! normalized 2nd-order flux
- ! s'r'c/2Sigma_s2 at t-1 multiplied by Lambda_3
-!
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PFRAC_ICE ! ri fraction of rc+ri
-!
-REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PRTHLS
-REAL, DIMENSION(:,:,:,:), INTENT(INOUT) :: PRRS ! var. at t+1 -split-
-!
-END SUBROUTINE TURB_HOR_THERMO_FLUX
-!
-END INTERFACE
-!
-END MODULE MODI_TURB_HOR_THERMO_FLUX
-! ################################################################
- SUBROUTINE TURB_HOR_THERMO_FLUX(KSPLT, KRR, KRRL, KRRI, &
- OTURB_FLX,OSUBG_COND, &
- TPFILE, &
- PK,PINV_PDXX,PINV_PDYY,PINV_PDZZ,PMZM_PRHODJ, &
- PDXX,PDYY,PDZZ,PDZX,PDZY, &
- PDIRCOSXW,PDIRCOSYW, &
- PRHODJ, &
- PSFTHM,PSFRM, &
- PWM,PTHLM,PRM, &
- PATHETA,PAMOIST,PSRCM,PFRAC_ICE, &
- PRTHLS,PRRS )
-! ################################################################
-!
-!
-!!**** *TURB_HOR* -routine to compute the source terms in the meso-NH
-!! model equations due to the non-vertical turbulent fluxes.
-!!
-!! PURPOSE
-!! -------
-!!
-!! see TURB_HOR
-!!
-!!** METHOD
-!! ------
-!!
-!! EXTERNAL
-!! --------
-!!
-!! IMPLICIT ARGUMENTS
-!! ------------------
-!!
-!! REFERENCE
-!! ---------
-!!
-!! AUTHOR
-!! ------
-!!
-!! Joan Cuxart * INM and Meteo-France *
-!!
-!! MODIFICATIONS
-!! -------------
-!! Aug , 1997 (V. Saravane) spliting of TURB_HOR
-!! Nov 27, 1997 (V. Masson) clearing of the routine
-!! Feb. 18, 1998 (J. Stein) bug for v'RC'
-!! Oct 18, 2000 (V. Masson) LES computations + LFLAT switch
-!! Nov 06, 2002 (V. Masson) LES budgets
-!! Feb 20, 2003 (JP Pinty) Add PFRAC_ICE
-!! October 2009 (G. Tanguy) add ILENCH=LEN(YCOMMENT) after
-!! change of YCOMMENT
-!! 04/2016 (M.Moge) Use openACC directives to port the TURB part of Meso-NH on GPU
-!! Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
-!! --------------------------------------------------------------------------
-!
-!* 0. DECLARATIONS
-! ------------
-!
-USE MODD_CST
-USE MODD_CONF
-USE MODD_CTURB
-use modd_field, only: tfielddata, TYPEREAL
-USE MODD_IO, ONLY: TFILEDATA
-USE MODD_PARAMETERS
-USE MODD_LES
-!
-USE MODE_IO_FIELD_WRITE, only: IO_Field_write
-use mode_mppdb
-!
-USE MODI_GRADIENT_M
-USE MODI_GRADIENT_U
-USE MODI_GRADIENT_V
-USE MODI_GRADIENT_W
-#ifndef MNH_OPENACC
-USE MODI_SHUMAN
-#else
-USE MODI_SHUMAN_DEVICE
-#endif
-USE MODI_LES_MEAN_SUBGRID
-!
-USE MODI_SECOND_MNH
-!
-#ifdef MNH_OPENACC
-USE MODE_MNH_ZWORK, ONLY: MNH_MEM_GET, MNH_MEM_POSITION_PIN, MNH_MEM_RELEASE
-#endif
-!
-IMPLICIT NONE
-!
-!
-!* 0.1 declaration of arguments
-!
-!
-!
-INTEGER, INTENT(IN) :: KSPLT ! split process index
-INTEGER, INTENT(IN) :: KRR ! number of moist var.
-INTEGER, INTENT(IN) :: KRRL ! number of liquid water var.
-INTEGER, INTENT(IN) :: KRRI ! number of ice water var.
-LOGICAL, INTENT(IN) :: OTURB_FLX ! switch to write the
- ! turbulent fluxes in the syncronous FM-file
-LOGICAL, INTENT(IN) :: OSUBG_COND ! Switch for sub-grid
-! condensation
-TYPE(TFILEDATA), INTENT(INOUT) :: TPFILE ! Output file
-!
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PK ! Turbulent diffusion doef.
- ! PK = PLM * SQRT(PTKEM)
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PINV_PDXX ! 1./PDXX
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PINV_PDYY ! 1./PDYY
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PINV_PDZZ ! 1./PDZZ
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PMZM_PRHODJ ! MZM(PRHODJ)
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX, PDYY, PDZZ, PDZX, PDZY
- ! Metric coefficients
-REAL, DIMENSION(:,:), INTENT(IN) :: PDIRCOSXW, PDIRCOSYW
-! Director Cosinus along x, y and z directions at surface w-point
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODJ ! density * grid volume
-!
-REAL, DIMENSION(:,:), INTENT(IN) :: PSFTHM,PSFRM
-!
-! Variables at t-1
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PWM
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHLM
-REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PRM ! mixing ratios at t-1,
- ! where PRM(:,:,:,1) = conservative mixing ratio
-!
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PATHETA ! coefficients between
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PAMOIST ! s and Thetal and Rnp
-
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PSRCM
- ! normalized 2nd-order flux
- ! s'r'c/2Sigma_s2 at t-1 multiplied by Lambda_3
-!
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PFRAC_ICE ! ri fraction of rc+ri
-!
-REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PRTHLS
-REAL, DIMENSION(:,:,:,:), INTENT(INOUT) :: PRRS ! var. at t+1 -split-
-!
-!
-!* 0.2 declaration of local variables
-!
-REAL, DIMENSION(:,:,:), pointer , contiguous :: ZFLX,ZFLXC ! work arrays
-!! REAL, DIMENSION(:,:,:), pointer , contiguous :: ZVPTV
-INTEGER :: IKB,IKE,IKU
- ! Index values for the Beginning and End
- ! mass points of the domain
-REAL, DIMENSION(:,:,:), pointer , contiguous :: ZCOEFF
- ! coefficients for the uncentred gradient
- ! computation near the ground
-!
-REAL :: ZTIME1, ZTIME2
-!
-#ifdef MNH_OPENACC
-REAL, DIMENSION(:,:,:), pointer , contiguous :: ZTMP1_DEVICE, ZTMP2_DEVICE, ZTMP3_DEVICE, ZTMP4_DEVICE
-REAL, DIMENSION(:,:,:), pointer , contiguous :: ZTMP5_DEVICE, ZTMP6_DEVICE, ZTMP7_DEVICE, ZTMP8_DEVICE
-#endif
-!
-TYPE(TFIELDDATA) :: TZFIELD
-!
-INTEGER :: JIU,JJU,JKU
-INTEGER :: JI,JJ,JK
-!
-! ---------------------------------------------------------------------------
-
-!$acc data present( PK, PINV_PDXX, PINV_PDYY, PINV_PDZZ, PMZM_PRHODJ, &
-!$acc & PDXX, PDYY, PDZZ, PDZX, PDZY, &
-!$acc & PDIRCOSXW, PDIRCOSYW, &
-!$acc & PRHODJ, &
-!$acc & PSFTHM, PSFRM, &
-!$acc & PWM, PTHLM, PRM, &
-!$acc & PATHETA, PAMOIST, PSRCM, PFRAC_ICE, &
-!$acc & PRTHLS, PRRS )
-
-if ( mppdb_initialized ) then
- !Check all in arrays
- call Mppdb_check( pinv_pdxx, "Turb_hor_thermo_flux beg:pinv_pdxx" )
- call Mppdb_check( pinv_pdyy, "Turb_hor_thermo_flux beg:pinv_pdyy" )
- call Mppdb_check( pinv_pdzz, "Turb_hor_thermo_flux beg:pinv_pdzz" )
- call Mppdb_check( pmzm_prhodj, "Turb_hor_thermo_flux beg:pmzm_prhodj" )
- call Mppdb_check( pdxx, "Turb_hor_thermo_flux beg:pdxx" )
- call Mppdb_check( pdyy, "Turb_hor_thermo_flux beg:pdyy" )
- call Mppdb_check( pdzz, "Turb_hor_thermo_flux beg:pdzz" )
- call Mppdb_check( pdzx, "Turb_hor_thermo_flux beg:pdzx" )
- call Mppdb_check( pdzy, "Turb_hor_thermo_flux beg:pdzy" )
- call Mppdb_check( pdircosxw, "Turb_hor_thermo_flux beg:pdircosxw" )
- call Mppdb_check( pdircosyw, "Turb_hor_thermo_flux beg:pdircosyw" )
- call Mppdb_check( prhodj, "Turb_hor_thermo_flux beg:prhodj" )
- call Mppdb_check( psfthm, "Turb_hor_thermo_flux beg:psfthm" )
- call Mppdb_check( psfrm, "Turb_hor_thermo_flux beg:psfrm" )
- call Mppdb_check( pwm, "Turb_hor_thermo_flux beg:pwm" )
- call Mppdb_check( pthlm, "Turb_hor_thermo_flux beg:pthlm" )
- call Mppdb_check( prm, "Turb_hor_thermo_flux beg:prm" )
- call Mppdb_check( patheta, "Turb_hor_thermo_flux beg:patheta" )
- call Mppdb_check( pamoist, "Turb_hor_thermo_flux beg:pamoist" )
- call Mppdb_check( psrcm, "Turb_hor_thermo_flux beg:psrcm" )
- call Mppdb_check( pfrac_ice, "Turb_hor_thermo_flux beg:pfrac_ice" )
- !Check all inout arrays
- call Mppdb_check( prthls, "Turb_hor_thermo_flux beg:prthls" )
- call Mppdb_check( prrs, "Turb_hor_thermo_flux beg:prrs" )
-end if
-
-JIU = size(pthlm, 1 )
-JJU = size(pthlm, 2 )
-JKU = size(pthlm, 3 )
-
-#ifndef MNH_OPENACC
-allocate( zflx (JIU,JJU,JKU) )
-allocate( zflxc(JIU,JJU,JKU) )
-! allocate( zvptv(JIU,JJU,JKU) )
-
-allocate( zcoeff(JIU,JJU, 1 + jpvext : 3 + jpvext ) )
-#else
-!Pin positions in the pools of MNH memory
-CALL MNH_MEM_POSITION_PIN()
-
-CALL MNH_MEM_GET( zflx, JIU, JJU, JKU )
-CALL MNH_MEM_GET( zflxc, JIU, JJU, JKU )
-! CALL MNH_MEM_GET( zvptv, JIU, JJU, JKU )
-
-CALL MNH_MEM_GET( zcoeff, 1, JIU, 1, JJU, 1 + jpvext, 3 + jpvext )
-
-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 )
-CALL MNH_MEM_GET( ztmp5_device, JIU, JJU, JKU )
-CALL MNH_MEM_GET( ztmp6_device, JIU, JJU, JKU )
-CALL MNH_MEM_GET( ztmp7_device, JIU, JJU, JKU )
-CALL MNH_MEM_GET( ztmp8_device, JIU, JJU, JKU )
-#endif
-
-!$acc data present( ZFLX, ZFLXC, ZCOEFF, &
-!$acc & ZTMP1_DEVICE, ZTMP2_DEVICE, ZTMP3_DEVICE, ZTMP4_DEVICE, &
-!$acc & ZTMP5_DEVICE, ZTMP6_DEVICE, ZTMP7_DEVICE, ZTMP8_DEVICE )
-
-!
-!* 1. PRELIMINARY COMPUTATIONS
-! ------------------------
-!
-IKB = 1+JPVEXT
-IKE = SIZE(PTHLM,3)-JPVEXT
-IKU = SIZE(PTHLM,3)
-!
-!
-! compute the coefficients for the uncentred gradient computation near the
-! ground
-!$acc kernels
-ZCOEFF(:,:,IKB+2)= - PDZZ(:,:,IKB+1) / &
- ( (PDZZ(:,:,IKB+2)+PDZZ(:,:,IKB+1)) * PDZZ(:,:,IKB+2) )
-ZCOEFF(:,:,IKB+1)= (PDZZ(:,:,IKB+2)+PDZZ(:,:,IKB+1)) / &
- ( PDZZ(:,:,IKB+1) * PDZZ(:,:,IKB+2) )
-ZCOEFF(:,:,IKB)= - (PDZZ(:,:,IKB+2)+2.*PDZZ(:,:,IKB+1)) / &
- ( (PDZZ(:,:,IKB+2)+PDZZ(:,:,IKB+1)) * PDZZ(:,:,IKB+1) )
-!$acc end kernels
-!
-!* 2. < U' THETA'l >
-! --------------
-!
-!
-#ifndef MNH_OPENACC
-ZFLX(:,:,:) = -XCSHF * MXM( PK ) * GX_M_U(1,IKU,1,PTHLM,PDXX,PDZZ,PDZX)
-ZFLX(:,:,IKE+1) = ZFLX(:,:,IKE)
-#else
-CALL MXM_DEVICE( PK, ZTMP1_DEVICE )
-CALL GX_M_U_DEVICE(1,IKU,1,PTHLM,PDXX,PDZZ,PDZX,ZTMP2_DEVICE)
-!$acc kernels
-!$acc_nv loop independent collapse(3)
-DO CONCURRENT ( JI=1:JIU,JJ=1:JJU,JK=1:JKU)
- ZFLX(JI,JJ,JK) = -XCSHF * ZTMP1_DEVICE(JI,JJ,JK) * ZTMP2_DEVICE(JI,JJ,JK)
-END DO
-ZFLX(:,:,IKE+1) = ZFLX(:,:,IKE)
-!$acc end kernels
-#endif
-!
-! Compute the flux at the first inner U-point with an uncentred vertical
-! gradient
-#ifndef MNH_OPENACC
-ZFLX(:,:,IKB:IKB) = -XCSHF * MXM( PK(:,:,IKB:IKB) ) * &
- ( DXM(PTHLM(:,:,IKB:IKB)) * PINV_PDXX(:,:,IKB:IKB) &
- -MXM( ZCOEFF(:,:,IKB+2:IKB+2)*PTHLM(:,:,IKB+2:IKB+2) &
- +ZCOEFF(:,:,IKB+1:IKB+1)*PTHLM(:,:,IKB+1:IKB+1) &
- +ZCOEFF(:,:,IKB :IKB )*PTHLM(:,:,IKB :IKB )) &
- *0.5* ( PDZX(:,:,IKB+1:IKB+1)+PDZX(:,:,IKB:IKB)) &
- * PINV_PDXX(:,:,IKB:IKB) )
-#else
-CALL MXM_DEVICE( PK(:,:,IKB:IKB), ZTMP1_DEVICE(:,:,1:1) )
-CALL DXM_DEVICE( PTHLM(:,:,IKB:IKB), ZTMP2_DEVICE(:,:,1:1) )
-!$acc kernels present_cr(ZTMP3_DEVICE)
-ZTMP3_DEVICE(:,:,1) = ZCOEFF(:,:,IKB+2)*PTHLM(:,:,IKB+2) &
- +ZCOEFF(:,:,IKB+1)*PTHLM(:,:,IKB+1) &
- +ZCOEFF(:,:,IKB )*PTHLM(:,:,IKB )
-!$acc end kernels
-CALL MXM_DEVICE( ZTMP3_DEVICE(:,:,1:1), ZTMP4_DEVICE(:,:,1:1))
-!$acc kernels present_cr(ZFLX,ZTMP1_DEVICE)
-ZFLX(:,:,IKB) = -XCSHF * ZTMP1_DEVICE(:,:,1) * &
- ( ZTMP2_DEVICE(:,:,1) * PINV_PDXX(:,:,IKB) - ZTMP4_DEVICE(:,:,1) &
- *0.5* ( PDZX(:,:,IKB+1)+PDZX(:,:,IKB)) &
- * PINV_PDXX(:,:,IKB) )
-! acc end kernels
-#endif
-! extrapolates the flux under the ground so that the vertical average with
-! the IKB flux gives the ground value ( warning the tangential surface
-! flux has been set to 0 for the moment !! to be improved )
-#ifndef MNH_OPENACC
-ZFLX(:,:,IKB-1:IKB-1) = 2. * MXM( SPREAD( PSFTHM(:,:)* PDIRCOSXW(:,:), 3,1) ) &
- - ZFLX(:,:,IKB:IKB)
-#else
-! acc kernels
-DO CONCURRENT ( JI=1:JIU,JJ=1:JJU )
- ZTMP1_DEVICE(JI,JJ,1) = PSFTHM(JI,JJ)* PDIRCOSXW(JI,JJ)
-END DO
-!$acc end kernels
- CALL MXM_DEVICE( ZTMP1_DEVICE(:,:,1:1), ZTMP2_DEVICE(:,:,1:1) )
-!$acc kernels present_cr(ZFLX)
- ZFLX(:,:,IKB-1) = 2. * ZTMP2_DEVICE(:,:,1) - ZFLX(:,:,IKB)
-!$acc end kernels
-#endif
-!
-! Add this source to the Theta_l sources
-!
-#ifndef MNH_OPENACC
-IF (.NOT. LFLAT) THEN
- PRTHLS(:,:,:) = PRTHLS(:,:,:) &
- - DXF( MXM(PRHODJ) * ZFLX(:,:,:) * PINV_PDXX(:,:,:) ) &
- + DZF( PMZM_PRHODJ(:,:,:) *MXF(PDZX*(MZM(ZFLX(:,:,:) * PINV_PDXX(:,:,:)))) * PINV_PDZZ(:,:,:) )
-ELSE
- PRTHLS(:,:,:) = PRTHLS(:,:,:) - DXF( MXM(PRHODJ) * ZFLX(:,:,:) * PINV_PDXX(:,:,:) )
-END IF
-#else
-IF (.NOT. LFLAT) THEN
- CALL MXM_DEVICE(PRHODJ, ZTMP1_DEVICE)
- !$acc kernels
- !$acc_nv loop independent collapse(3)
- DO CONCURRENT ( JI=1:JIU,JJ=1:JJU,JK=1:JKU)
- ZTMP2_DEVICE(JI,JJ,JK) = ZTMP1_DEVICE(JI,JJ,JK) * ZFLX(JI,JJ,JK) * PINV_PDXX(JI,JJ,JK)
- END DO
- !$acc end kernels
- CALL DXF_DEVICE(ZTMP2_DEVICE, ZTMP3_DEVICE)
- !$acc kernels
- !$acc_nv loop independent collapse(3)
- DO CONCURRENT ( JI=1:JIU,JJ=1:JJU,JK=1:JKU)
- ZTMP2_DEVICE(JI,JJ,JK) = ZFLX(JI,JJ,JK) * PINV_PDXX(JI,JJ,JK)
- END DO
- !$acc end kernels
- CALL MZM_DEVICE(ZTMP2_DEVICE,ZTMP4_DEVICE)
- !$acc kernels
- !$acc_nv loop independent collapse(3)
- DO CONCURRENT ( JI=1:JIU,JJ=1:JJU,JK=1:JKU)
- ZTMP2_DEVICE(JI,JJ,JK) = PDZX(JI,JJ,JK)*ZTMP4_DEVICE(JI,JJ,JK)
- END DO
-!$acc end kernels
- CALL MXF_DEVICE(ZTMP2_DEVICE, ZTMP4_DEVICE)
-!$acc kernels
-!$acc_nv loop independent collapse(3)
-DO CONCURRENT ( JI=1:JIU,JJ=1:JJU,JK=1:JKU)
- ZTMP2_DEVICE(JI,JJ,JK) = PMZM_PRHODJ(JI,JJ,JK) * ZTMP4_DEVICE(JI,JJ,JK) * PINV_PDZZ(JI,JJ,JK)
-END DO
-!$acc end kernels
- CALL DZF_DEVICE( ZTMP2_DEVICE, ZTMP4_DEVICE )
-!$acc kernels
- PRTHLS(:,:,:) = PRTHLS(:,:,:) - ZTMP3_DEVICE(:,:,:) + ZTMP4_DEVICE(:,:,:)
-!$acc end kernels
-ELSE
- CALL MXM_DEVICE(PRHODJ, ZTMP1_DEVICE)
-!$acc kernels
- ZTMP2_DEVICE(:,:,:) = ZTMP1_DEVICE(:,:,:) * ZFLX(:,:,:) * PINV_PDXX(:,:,:)
-!$acc end kernels
- CALL DXF_DEVICE(ZTMP2_DEVICE, ZTMP3_DEVICE)
-!$acc kernels
- PRTHLS(:,:,:) = PRTHLS(:,:,:) - ZTMP3_DEVICE(:,:,:)
-!$acc end kernels
-END IF
-#endif
-!
-! Compute the equivalent tendancy for Rc and Ri
-!
-#ifndef MNH_OPENACC
-IF ( KRRL >= 1 ) THEN
- IF (.NOT. LFLAT) THEN
- ZFLXC(:,:,:) = 2.*( MXF( MXM( PRHODJ(:,:,:)*PATHETA(:,:,:)*PSRCM(:,:,:) )*ZFLX(:,:,:) ) &
- +MZF( MZM( PRHODJ(:,:,:)*PATHETA(:,:,:)*PSRCM(:,:,:) )*MXF( &
- PDZX*(MZM( ZFLX(:,:,:)*PINV_PDXX(:,:,:) )) ) )&
- )
- IF ( KRRI >= 1 ) THEN
- PRRS(:,:,:,2) = PRRS(:,:,:,2) + 2. * &
- (- DXF( MXM( PRHODJ(:,:,:)*PATHETA(:,:,:)*PSRCM(:,:,:) )*ZFLX(:,:,:)*PINV_PDXX(:,:,:) ) &
- + DZF( MZM( PRHODJ(:,:,:)*PATHETA(:,:,:)*PSRCM(:,:,:) )* &
- MXF( PDZX*(MZM( ZFLX(:,:,:)*PINV_PDXX(:,:,:) )) )&
- *PINV_PDZZ(:,:,:) ) &
- )*(1.0-PFRAC_ICE(:,:,:))
- PRRS(:,:,:,4) = PRRS(:,:,:,4) + 2. * &
- (- DXF( MXM( PRHODJ(:,:,:)*PATHETA(:,:,:)*PSRCM(:,:,:) )*ZFLX(:,:,:)*PINV_PDXX(:,:,:) ) &
- + DZF( MZM( PRHODJ(:,:,:)*PATHETA(:,:,:)*PSRCM(:,:,:) )* &
- MXF( PDZX*(MZM( ZFLX(:,:,:)*PINV_PDXX(:,:,:) )) )&
- *PINV_PDZZ(:,:,:) ) &
- )*PFRAC_ICE(:,:,:)
- ELSE
- PRRS(:,:,:,2) = PRRS(:,:,:,2) + 2. * &
- (- DXF( MXM( PRHODJ(:,:,:)*PATHETA(:,:,:)*PSRCM(:,:,:) )*ZFLX(:,:,:)*PINV_PDXX(:,:,:) ) &
- + DZF( MZM( PRHODJ(:,:,:)*PATHETA(:,:,:)*PSRCM(:,:,:) )* &
- MXF( PDZX*(MZM( ZFLX(:,:,:)*PINV_PDXX(:,:,:) )) )&
- *PINV_PDZZ(:,:,:) ) &
- )
- END IF
- ELSE
- ZFLXC(:,:,:) = 2.*MXF( MXM( PRHODJ(:,:,:)*PATHETA(:,:,:)*PSRCM(:,:,:) )*ZFLX )
- IF ( KRRI >= 1 ) THEN
- PRRS(:,:,:,2) = PRRS(:,:,:,2) - 2. * &
- DXF( MXM( PRHODJ(:,:,:)*PATHETA(:,:,:)*PSRCM(:,:,:) )*ZFLX(:,:,:)*PINV_PDXX(:,:,:) )*(1.0-PFRAC_ICE(:,:,:))
- PRRS(:,:,:,4) = PRRS(:,:,:,4) - 2. * &
- DXF( MXM( PRHODJ(:,:,:)*PATHETA(:,:,:)*PSRCM(:,:,:) )*ZFLX(:,:,:)*PINV_PDXX(:,:,:) )*PFRAC_ICE(:,:,:)
- ELSE
- PRRS(:,:,:,2) = PRRS(:,:,:,2) - 2. * &
- DXF( MXM( PRHODJ(:,:,:)*PATHETA(:,:,:)*PSRCM(:,:,:) )*ZFLX(:,:,:)*PINV_PDXX(:,:,:) )
- END IF
- END IF
-END IF
-#else
-IF ( KRRL >= 1 ) THEN
- IF (.NOT. LFLAT) THEN
- !$acc kernels
- ZTMP1_DEVICE(:,:,:) = PRHODJ(:,:,:)*PATHETA(:,:,:)*PSRCM(:,:,:)
- !$acc end kernels
- CALL MZM_DEVICE( ZTMP1_DEVICE, ZTMP4_DEVICE )
- CALL MXM_DEVICE( ZTMP1_DEVICE, ZTMP2_DEVICE )
- !$acc kernels present_cr(ZTMP1_DEVICE)
- ZTMP1_DEVICE(:,:,:) = ZTMP2_DEVICE(:,:,:) *ZFLX(:,:,:)
- !$acc end kernels
- CALL MXF_DEVICE( ZTMP1_DEVICE, ZTMP2_DEVICE)
- !$acc kernels present_cr(ZTMP1_DEVICE)
- ZTMP1_DEVICE(:,:,:) = ZFLX(:,:,:)*PINV_PDXX(:,:,:)
- !$acc end kernels
- CALL MZM_DEVICE( ZTMP1_DEVICE, ZTMP5_DEVICE )
- !$acc kernels present_cr(ZTMP6_DEVICE)
- ZTMP6_DEVICE(:,:,:) = PDZX(:,:,:)*ZTMP5_DEVICE(:,:,:)
- !$acc end kernels
- CALL MXF_DEVICE( ZTMP6_DEVICE, ZTMP5_DEVICE )
- !$acc kernels present_cr(ZTMP6_DEVICE)
- ZTMP6_DEVICE(:,:,:) = ZTMP4_DEVICE(:,:,:)*ZTMP5_DEVICE(:,:,:)
- !$acc end kernels
- CALL MZF_DEVICE( ZTMP6_DEVICE, ZTMP7_DEVICE )
- !$acc kernels present_cr(ZFLXC)
- ZFLXC(:,:,:) = 2.*( ZTMP2_DEVICE(:,:,:) +ZTMP7_DEVICE(:,:,:) )
- !$acc end kernels
- IF ( KRRI >= 1 ) THEN
- !$acc kernels
- ZTMP1_DEVICE(:,:,:) = PRHODJ(:,:,:)*PATHETA(:,:,:)*PSRCM(:,:,:)
- !$acc end kernels
- CALL MXM_DEVICE( ZTMP1_DEVICE, ZTMP2_DEVICE )
- !$acc kernels present_cr(ZTMP6_DEVICE)
- ZTMP6_DEVICE(:,:,:) = ZTMP2_DEVICE(:,:,:)*ZFLX(:,:,:)*PINV_PDXX(:,:,:)
- !$acc end kernels
- CALL DXF_DEVICE( ZTMP6_DEVICE, ZTMP2_DEVICE)
- !$acc kernels present_cr(ZTMP3_DEVICE)
- ZTMP3_DEVICE(:,:,:) = ZTMP4_DEVICE(:,:,:)*ZTMP5_DEVICE(:,:,:)*PINV_PDZZ(:,:,:)
- !$acc end kernels
- CALL DZF_DEVICE( ZTMP3_DEVICE, ZTMP4_DEVICE )
-!$acc kernels
- PRRS(:,:,:,2) = PRRS(:,:,:,2) + 2. * (- ZTMP2_DEVICE(:,:,:) + ZTMP4_DEVICE(:,:,:) )*(1.0-PFRAC_ICE(:,:,:))
- PRRS(:,:,:,4) = PRRS(:,:,:,4) + 2. * (- ZTMP2_DEVICE(:,:,:) + ZTMP4_DEVICE(:,:,:) )*PFRAC_ICE(:,:,:)
-!$acc end kernels
- ELSE
- !$acc kernels
- ZTMP1_DEVICE(:,:,:) = PRHODJ(:,:,:)*PATHETA(:,:,:)*PSRCM(:,:,:)
- !$acc end kernels
- CALL MXM_DEVICE( ZTMP1_DEVICE, ZTMP2_DEVICE )
- !$acc kernels
- ZTMP1_DEVICE(:,:,:) = ZTMP2_DEVICE(:,:,:) *ZFLX(:,:,:)*PINV_PDXX(:,:,:)
- !$acc end kernels
- CALL DXF_DEVICE( ZTMP6_DEVICE, ZTMP2_DEVICE)
- !$acc kernels
- ZTMP3_DEVICE(:,:,:) = ZTMP4_DEVICE(:,:,:)*ZTMP5_DEVICE(:,:,:)*PINV_PDZZ(:,:,:)
- !$acc end kernels
- CALL DZF_DEVICE( ZTMP3_DEVICE, ZTMP4_DEVICE )
-!$acc kernels
- PRRS(:,:,:,2) = PRRS(:,:,:,2) + 2. * (- ZTMP2_DEVICE(:,:,:) + ZTMP4_DEVICE(:,:,:) )
-!$acc end kernels
- END IF
- ELSE
- !$acc kernels
- ZTMP1_DEVICE(:,:,:) = PRHODJ(:,:,:)*PATHETA(:,:,:)*PSRCM(:,:,:)
- !$acc end kernels
- CALL MXM_DEVICE( ZTMP1_DEVICE,ZTMP2_DEVICE )
- !$acc kernels
- ZTMP3_DEVICE(:,:,:) = ZTMP2_DEVICE(:,:,:)*ZFLX(:,:,:)
- !$acc end kernels
- CALL MXF_DEVICE( ZTMP3_DEVICE, ZTMP4_DEVICE )
- !$acc kernels
- ZFLXC(:,:,:) = 2.*ZTMP4_DEVICE(:,:,:)
- !$acc end kernels
- IF ( KRRI >= 1 ) THEN
- !$acc kernels
- ZTMP1_DEVICE(:,:,:) = ZTMP2_DEVICE(:,:,:)*ZFLX(:,:,:)*PINV_PDXX(:,:,:)
- !$acc end kernels
- CALL DXF_DEVICE( ZTMP1_DEVICE, ZTMP2_DEVICE )
- !$acc kernels
- PRRS(:,:,:,2) = PRRS(:,:,:,2) - 2. * ZTMP2_DEVICE(:,:,:)*(1.0-PFRAC_ICE(:,:,:))
- PRRS(:,:,:,4) = PRRS(:,:,:,4) - 2. * ZTMP2_DEVICE(:,:,:)*PFRAC_ICE(:,:,:)
- !$acc end kernels
- ELSE
- !$acc kernels
- ZTMP1_DEVICE(:,:,:) = ZTMP2_DEVICE(:,:,:)*ZFLX(:,:,:)*PINV_PDXX(:,:,:)
- !$acc end kernels
- CALL DXF_DEVICE( ZTMP1_DEVICE, ZTMP2_DEVICE )
- !$acc kernels
- PRRS(:,:,:,2) = PRRS(:,:,:,2) - 2. * ZTMP2_DEVICE(:,:,:)
- !$acc end kernels
- END IF
- END IF
-END IF
-#endif
-!
-!! stores this flux in ZWORK to compute later <U' VPT'>
-!!ZWORK(:,:,:) = ZFLX(:,:,:)
-!
-! stores the horizontal <U THl>
-IF ( tpfile%lopened .AND. OTURB_FLX ) THEN
- TZFIELD%CMNHNAME = 'UTHL_FLX'
- TZFIELD%CSTDNAME = ''
- TZFIELD%CLONGNAME = 'UTHL_FLX'
- TZFIELD%CUNITS = 'K m s-1'
- TZFIELD%CDIR = 'XY'
- TZFIELD%CCOMMENT = 'X_Y_Z_UTHL_FLX'
- TZFIELD%NGRID = 2
- TZFIELD%NTYPE = TYPEREAL
- TZFIELD%NDIMS = 3
- TZFIELD%LTIMEDEP = .TRUE.
-!$acc update self(ZFLX)
- CALL IO_Field_write(TPFILE,TZFIELD,ZFLX(:,:,:))
-END IF
-!
-IF (KSPLT==1 .AND. LLES_CALL) THEN
- CALL SECOND_MNH(ZTIME1)
-#ifndef MNH_OPENACC
- CALL LES_MEAN_SUBGRID( MXF(ZFLX), X_LES_SUBGRID_UThl )
- CALL LES_MEAN_SUBGRID( MZF(MXF(GX_W_UW(PWM,PDXX,PDZZ,PDZX)*MZM(ZFLX))),&
- X_LES_RES_ddxa_W_SBG_UaThl , .TRUE. )
- CALL LES_MEAN_SUBGRID( GX_M_M(PTHLM,PDXX,PDZZ,PDZX)*MXF(ZFLX),&
- X_LES_RES_ddxa_Thl_SBG_UaThl , .TRUE. )
- IF (KRR>=1) THEN
- CALL LES_MEAN_SUBGRID( GX_M_M(PRM(:,:,:,1),PDXX,PDZZ,PDZX)*MXF(ZFLX), &
- X_LES_RES_ddxa_Rt_SBG_UaThl , .TRUE. )
- END IF
-#else
-!$acc data copy(X_LES_SUBGRID_UThl,X_LES_RES_ddxa_W_SBG_UaThl, &
-!$acc & X_LES_RES_ddxa_Thl_SBG_UaThl,X_LES_RES_ddxa_Rt_SBG_UaThl)
- !
- CALL MXF_DEVICE(ZFLX,ZTMP1_DEVICE)
- CALL LES_MEAN_SUBGRID( ZTMP1_DEVICE, X_LES_SUBGRID_UThl )
- !
- CALL GX_W_UW_DEVICE(PWM,PDXX,PDZZ,PDZX,ZTMP1_DEVICE)
- CALL MZM_DEVICE(ZFLX,ZTMP2_DEVICE)
- !$acc kernels
- ZTMP3_DEVICE(:,:,:) = ZTMP1_DEVICE(:,:,:)*ZTMP2_DEVICE(:,:,:)
- !$acc end kernels
- CALL MXF_DEVICE(ZTMP3_DEVICE,ZTMP1_DEVICE)
- CALL MZF_DEVICE(ZTMP1_DEVICE,ZTMP2_DEVICE)
- CALL LES_MEAN_SUBGRID( ZTMP2_DEVICE,X_LES_RES_ddxa_W_SBG_UaThl , .TRUE. )
- !
- CALL GX_M_M_DEVICE(PTHLM,PDXX,PDZZ,PDZX,ZTMP1_DEVICE)
- CALL MXF_DEVICE(ZFLX,ZTMP2_DEVICE)
- !$acc kernels
- ZTMP3_DEVICE(:,:,:) = ZTMP1_DEVICE(:,:,:) * ZTMP2_DEVICE(:,:,:)
- !$acc end kernels
- CALL LES_MEAN_SUBGRID( ZTMP3_DEVICE,X_LES_RES_ddxa_Thl_SBG_UaThl , .TRUE. )
- !
- IF (KRR>=1) THEN
- CALL GX_M_M_DEVICE(PRM(:,:,:,1),PDXX,PDZZ,PDZX,ZTMP1_DEVICE)
- !$acc kernels
- ZTMP3_DEVICE(:,:,:) = ZTMP1_DEVICE(:,:,:) * ZTMP2_DEVICE(:,:,:)
- !$acc end kernels
- CALL LES_MEAN_SUBGRID( ZTMP3_DEVICE,X_LES_RES_ddxa_Rt_SBG_UaThl , .TRUE. )
- END IF
-!$acc end data
-
-#endif
-
- CALL SECOND_MNH(ZTIME2)
- XTIME_LES = XTIME_LES + ZTIME2 - ZTIME1
-END IF
-!
-!* 3. < U' R'np >
-! -----------
-IF (KRR/=0) THEN
- !
-#ifndef MNH_OPENACC
- ZFLX(:,:,:) = -XCHF * MXM( PK ) * GX_M_U(1,IKU,1,PRM(:,:,:,1),PDXX,PDZZ,PDZX)
- ZFLX(:,:,IKE+1) = ZFLX(:,:,IKE)
-!
-! Compute the flux at the first inner U-point with an uncentred vertical
-! gradient
- ZFLX(:,:,IKB:IKB) = -XCHF * MXM( PK(:,:,IKB:IKB) ) * &
- ( DXM(PRM(:,:,IKB:IKB,1)) * PINV_PDXX(:,:,IKB:IKB) &
- -MXM( ZCOEFF(:,:,IKB+2:IKB+2)*PRM(:,:,IKB+2:IKB+2,1) &
- +ZCOEFF(:,:,IKB+1:IKB+1)*PRM(:,:,IKB+1:IKB+1,1) &
- +ZCOEFF(:,:,IKB :IKB )*PRM(:,:,IKB :IKB ,1)) &
- *0.5* ( PDZX(:,:,IKB+1:IKB+1)+PDZX(:,:,IKB:IKB)) &
- * PINV_PDXX(:,:,IKB:IKB) )
-! extrapolates the flux under the ground so that the vertical average with
-! the IKB flux gives the ground value ( warning the tangential surface
-! flux has been set to 0 for the moment !! to be improved )
- ZFLX(:,:,IKB-1:IKB-1) = 2. * MXM( SPREAD( PSFRM(:,:)* PDIRCOSXW(:,:), 3,1) ) &
- - ZFLX(:,:,IKB:IKB)
- !
- ! Add this source to the conservative mixing ratio sources
- !
- IF (.NOT. LFLAT) THEN
- PRRS(:,:,:,1) = PRRS(:,:,:,1) &
- - DXF( MXM(PRHODJ) * ZFLX(:,:,:) * PINV_PDXX(:,:,:) ) &
- + DZF( PMZM_PRHODJ(:,:,:) *MXF(PDZX*(MZM(ZFLX * PINV_PDXX(:,:,:)))) * PINV_PDZZ(:,:,:) )
- ELSE
- PRRS(:,:,:,1) = PRRS(:,:,:,1) - DXF( MXM(PRHODJ) * ZFLX(:,:,:) * PINV_PDXX(:,:,:) )
- END IF
- !
- ! Compute the equivalent tendancy for Rc and Ri
- !
- IF ( KRRL >= 1 ) THEN
- IF (.NOT. LFLAT) THEN
- ZFLXC(:,:,:) = ZFLXC(:,:,:) &
- + 2.*( MXF( MXM( PRHODJ(:,:,:)*PAMOIST(:,:,:)*PSRCM(:,:,:) )*ZFLX(:,:,:) ) &
- +MZF( MZM( PRHODJ(:,:,:)*PAMOIST(:,:,:)*PSRCM(:,:,:) )*MXF( &
- PDZX(:,:,:)*(MZM( ZFLX(:,:,:)*PINV_PDXX(:,:,:) )) ) )&
- )
- IF ( KRRI >= 1 ) THEN
- PRRS(:,:,:,2) = PRRS(:,:,:,2) + 2. * &
- (- DXF( MXM( PRHODJ(:,:,:)*PAMOIST(:,:,:)*PSRCM(:,:,:) )*ZFLX(:,:,:)*PINV_PDXX(:,:,:) ) &
- + DZF( MZM( PRHODJ(:,:,:)*PAMOIST(:,:,:)*PSRCM(:,:,:) )*MXF( PDZX(:,:,:)* &
- (MZM( ZFLX(:,:,:)*PINV_PDXX(:,:,:) )) )&
- *PINV_PDZZ(:,:,:) ) &
- )*(1.0-PFRAC_ICE(:,:,:))
- PRRS(:,:,:,2) = PRRS(:,:,:,2) + 2. * &
- (- DXF( MXM( PRHODJ(:,:,:)*PAMOIST(:,:,:)*PSRCM(:,:,:) )*ZFLX(:,:,:)*PINV_PDXX(:,:,:) ) &
- + DZF( MZM( PRHODJ(:,:,:)*PAMOIST(:,:,:)*PSRCM(:,:,:) )*MXF( PDZX(:,:,:)* &
- (MZM( ZFLX(:,:,:)*PINV_PDXX(:,:,:) )) )&
- *PINV_PDZZ(:,:,:) ) &
- )*PFRAC_ICE(:,:,:)
- ELSE
- PRRS(:,:,:,2) = PRRS(:,:,:,2) + 2. * &
- (- DXF( MXM( PRHODJ(:,:,:)*PAMOIST(:,:,:)*PSRCM(:,:,:) )*ZFLX(:,:,:)*PINV_PDXX(:,:,:) ) &
- + DZF( MZM( PRHODJ(:,:,:)*PAMOIST(:,:,:)*PSRCM(:,:,:) )*MXF( PDZX(:,:,:)* &
- (MZM( ZFLX(:,:,:)*PINV_PDXX(:,:,:) )) )&
- *PINV_PDZZ(:,:,:) ) &
- )
- END IF
- ELSE
- ZFLXC(:,:,:) = ZFLXC(:,:,:) + 2.*MXF( MXM( PRHODJ(:,:,:)*PAMOIST(:,:,:)*PSRCM(:,:,:) )*ZFLX(:,:,:) )
- IF ( KRRI >= 1 ) THEN
- PRRS(:,:,:,2) = PRRS(:,:,:,2) - 2. * &
- DXF( MXM( PRHODJ(:,:,:)*PAMOIST(:,:,:)*PSRCM(:,:,:) )*ZFLX(:,:,:)*PINV_PDXX(:,:,:) )*(1.0-PFRAC_ICE(:,:,:))
- PRRS(:,:,:,4) = PRRS(:,:,:,4) - 2. * &
- DXF( MXM( PRHODJ(:,:,:)*PAMOIST(:,:,:)*PSRCM(:,:,:) )*ZFLX(:,:,:)*PINV_PDXX(:,:,:) )*PFRAC_ICE(:,:,:)
- ELSE
- PRRS(:,:,:,2) = PRRS(:,:,:,2) - 2. * &
- DXF( MXM( PRHODJ(:,:,:)*PAMOIST(:,:,:)*PSRCM(:,:,:) )*ZFLX(:,:,:)*PINV_PDXX(:,:,:) )
- END IF
- END IF
- END IF
- !
- ! stores the horizontal <U Rnp>
- IF ( tpfile%lopened .AND. OTURB_FLX ) THEN
- TZFIELD%CMNHNAME = 'UR_FLX'
- TZFIELD%CSTDNAME = ''
- TZFIELD%CLONGNAME = 'UR_FLX'
- TZFIELD%CUNITS = 'kg kg-1 m s-1'
- TZFIELD%CDIR = 'XY'
- TZFIELD%CCOMMENT = 'X_Y_Z_UR_FLX'
- TZFIELD%NGRID = 2
- TZFIELD%NTYPE = TYPEREAL
- TZFIELD%NDIMS = 3
- TZFIELD%LTIMEDEP = .TRUE.
- CALL IO_Field_write(TPFILE,TZFIELD,ZFLX(:,:,:))
- END IF
- !
- IF (KSPLT==1 .AND. LLES_CALL) THEN
- CALL SECOND_MNH(ZTIME1)
- CALL LES_MEAN_SUBGRID( MXF(ZFLX), X_LES_SUBGRID_URt )
- CALL LES_MEAN_SUBGRID( MZF(MXF(GX_W_UW(PWM,PDXX,PDZZ,PDZX)*MZM(ZFLX))),&
- X_LES_RES_ddxa_W_SBG_UaRt , .TRUE. )
- CALL LES_MEAN_SUBGRID( GX_M_M(PTHLM,PDXX,PDZZ,PDZX)*MXF(ZFLX),&
- X_LES_RES_ddxa_Thl_SBG_UaRt , .TRUE. )
- CALL LES_MEAN_SUBGRID( GX_M_M(PRM(:,:,:,1),PDXX,PDZZ,PDZX)*MXF(ZFLX),&
- X_LES_RES_ddxa_Rt_SBG_UaRt , .TRUE. )
- CALL SECOND_MNH(ZTIME2)
- XTIME_LES = XTIME_LES + ZTIME2 - ZTIME1
- END IF
-!
- !
- IF (KRRL>0 .AND. KSPLT==1 .AND. LLES_CALL) THEN
- CALL SECOND_MNH(ZTIME1)
- CALL LES_MEAN_SUBGRID(MXF(ZFLXC), X_LES_SUBGRID_URc )
- CALL SECOND_MNH(ZTIME2)
- XTIME_LES = XTIME_LES + ZTIME2 - ZTIME1
- END IF
-!
-END IF
-#else
- CALL MXM_DEVICE( PK, ZTMP1_DEVICE )
- CALL GX_M_U_DEVICE(1,IKU,1,PRM(:,:,:,1),PDXX,PDZZ,PDZX,ZTMP2_DEVICE)
-!$acc kernels
-!$acc_nv loop independent collapse(3)
-DO CONCURRENT ( JI=1:JIU,JJ=1:JJU,JK=1:JKU)
- ZFLX(JI,JJ,JK) = -XCHF * ZTMP1_DEVICE(JI,JJ,JK) * ZTMP2_DEVICE(JI,JJ,JK)
-END DO
- ZFLX(:,:,IKE+1) = ZFLX(:,:,IKE)
-!$acc end kernels
-!
-! Compute the flux at the first inner U-point with an uncentred vertical
-! gradient
- CALL MXM_DEVICE( PK(:,:,IKB:IKB), ZTMP1_DEVICE(:,:,1:1) )
- CALL DXM_DEVICE(PRM(:,:,IKB:IKB,1), ZTMP2_DEVICE(:,:,1:1))
-!$acc kernels present_cr(ZTMP3_DEVICE)
- ZTMP3_DEVICE(:,:,1) = ZCOEFF(:,:,IKB+2)*PRM(:,:,IKB+2,1) &
- +ZCOEFF(:,:,IKB+1)*PRM(:,:,IKB+1,1) &
- +ZCOEFF(:,:,IKB )*PRM(:,:,IKB ,1)
-!$acc end kernels
- CALL MXM_DEVICE(ZTMP3_DEVICE(:,:,1:1),ZTMP4_DEVICE(:,:,1:1))
-!$acc kernels present_cr(ZFLX)
- ZFLX(:,:,IKB) = -XCHF * ZTMP1_DEVICE(:,:,1) * &
- ( ZTMP2_DEVICE(:,:,1) * PINV_PDXX(:,:,IKB) &
- -ZTMP4_DEVICE(:,:,1) &
- *0.5* ( PDZX(:,:,IKB+1)+PDZX(:,:,IKB)) &
- * PINV_PDXX(:,:,IKB) )
-! extrapolates the flux under the ground so that the vertical average with
-! the IKB flux gives the ground value ( warning the tangential surface
-! flux has been set to 0 for the moment !! to be improved )
- ZTMP1_DEVICE(:,:,1) = PSFRM(:,:)* PDIRCOSXW(:,:)
-!$acc end kernels
- CALL MXM_DEVICE(ZTMP1_DEVICE(:,:,1:1),ZTMP2_DEVICE(:,:,1:1))
- !$acc kernels present_cr(ZFLX)
- ZFLX(:,:,IKB-1) = 2. * ZTMP2_DEVICE(:,:,1) - ZFLX(:,:,IKB)
- !$acc end kernels
-
- !
- ! Add this source to the conservative mixing ratio sources
- !
- IF (.NOT. LFLAT) THEN
- CALL MXM_DEVICE(PRHODJ,ZTMP1_DEVICE)
- !$acc kernels
- !$acc_nv loop independent collapse(3)
- DO CONCURRENT ( JI=1:JIU,JJ=1:JJU,JK=1:JKU)
- ZTMP2_DEVICE(JI,JJ,JK) = ZTMP1_DEVICE(JI,JJ,JK) * ZFLX(JI,JJ,JK) * PINV_PDXX(JI,JJ,JK)
- END DO
- !$acc end kernels
- CALL DXF_DEVICE( ZTMP2_DEVICE, ZTMP3_DEVICE )
- !$acc kernels
- !$acc_nv loop independent collapse(3)
- DO CONCURRENT ( JI=1:JIU,JJ=1:JJU,JK=1:JKU)
- ZTMP2_DEVICE(JI,JJ,JK) = ZFLX(JI,JJ,JK) * PINV_PDXX(JI,JJ,JK)
- END DO
- !$acc end kernels
- CALL MZM_DEVICE(ZTMP2_DEVICE,ZTMP4_DEVICE)
- !$acc kernels
- !$acc_nv loop independent collapse(3)
- DO CONCURRENT ( JI=1:JIU,JJ=1:JJU,JK=1:JKU)
- ZTMP2_DEVICE(JI,JJ,JK) = PDZX(JI,JJ,JK)*ZTMP4_DEVICE(JI,JJ,JK)
- END DO
- !$acc end kernels
- CALL MXF_DEVICE(ZTMP2_DEVICE,ZTMP4_DEVICE)
- !$acc kernels
- !$acc_nv loop independent collapse(3)
- DO CONCURRENT ( JI=1:JIU,JJ=1:JJU,JK=1:JKU)
- ZTMP2_DEVICE(JI,JJ,JK) = PMZM_PRHODJ(JI,JJ,JK) * ZTMP4_DEVICE(JI,JJ,JK) * PINV_PDZZ(JI,JJ,JK)
- END DO
- !$acc end kernels
- CALL DZF_DEVICE(ZTMP2_DEVICE,ZTMP4_DEVICE)
- !$acc kernels
- !$acc_nv loop independent collapse(3)
- DO CONCURRENT ( JI=1:JIU,JJ=1:JJU,JK=1:JKU)
- PRRS(JI,JJ,JK,1) = PRRS(JI,JJ,JK,1) - ZTMP3_DEVICE(JI,JJ,JK) + ZTMP4_DEVICE(JI,JJ,JK)
- END DO
- !$acc end kernels
- ELSE
- CALL MXM_DEVICE(PRHODJ,ZTMP1_DEVICE)
- !$acc kernels
- ZTMP2_DEVICE(:,:,:) = ZTMP1_DEVICE(:,:,:) * ZFLX(:,:,:) * PINV_PDXX(:,:,:)
- !$acc end kernels
- CALL DXF_DEVICE( ZTMP2_DEVICE, ZTMP3_DEVICE )
- !$acc kernels
- PRRS(:,:,:,1) = PRRS(:,:,:,1) - ZTMP3_DEVICE(:,:,:)
- !$acc end kernels
- END IF
- !
- ! Compute the equivalent tendancy for Rc and Ri
- !
- IF ( KRRL >= 1 ) THEN
- !$acc kernels present_cr(ZTMP2_DEVICE)
- ZTMP1_DEVICE(:,:,:) = PRHODJ(:,:,:)*PAMOIST(:,:,:)*PSRCM(:,:,:)
- ZTMP2_DEVICE(:,:,:) = ZFLX(:,:,:)*PINV_PDXX(:,:,:)
- !$acc end kernels
- CALL MXM_DEVICE( ZTMP1_DEVICE, ZTMP8_DEVICE )
- IF (.NOT. LFLAT) THEN
- !$acc kernels present_cr(ZTMP4_DEVICE)
- ZTMP4_DEVICE(:,:,:) = ZTMP8_DEVICE(:,:,:) * ZFLX(:,:,:)
- !$acc end kernels
- CALL MXF_DEVICE( ZTMP4_DEVICE, ZTMP3_DEVICE )
- CALL MZM_DEVICE( ZTMP1_DEVICE, ZTMP4_DEVICE )
- CALL MZM_DEVICE( ZTMP2_DEVICE, ZTMP5_DEVICE )
- !$acc kernels present_cr(ZTMP6_DEVICE)
- ZTMP6_DEVICE(:,:,:) = PDZX(:,:,:)*ZTMP5_DEVICE(:,:,:)
- !$acc end kernels
- CALL MXF_DEVICE( ZTMP6_DEVICE, ZTMP5_DEVICE )
- !$acc kernels present_cr(ZTMP6_DEVICE)
- ZTMP6_DEVICE(:,:,:) = ZTMP4_DEVICE(:,:,:)*ZTMP5_DEVICE(:,:,:)
- !$acc end kernels
- CALL MZF_DEVICE( ZTMP6_DEVICE, ZTMP7_DEVICE )
- !$acc kernels present_cr(ZFLXC,ZTMP6_DEVICE)
- ZFLXC(:,:,:) = ZFLXC(:,:,:) + 2.*( ZTMP3_DEVICE(:,:,:) + ZTMP7_DEVICE(:,:,:) )
- !
- ZTMP6_DEVICE(:,:,:) = ZTMP4_DEVICE(:,:,:)*ZTMP5_DEVICE(:,:,:)*PINV_PDZZ(:,:,:)
- !$acc end kernels
- CALL DZF_DEVICE( ZTMP6_DEVICE, ZTMP3_DEVICE )
- !$acc kernels present_cr(ZTMP4_DEVICE)
- ZTMP4_DEVICE(:,:,:) = ZTMP8_DEVICE(:,:,:) * ZFLX(:,:,:)*PINV_PDXX(:,:,:)
- !$acc end kernels
- CALL DXF_DEVICE(ZTMP4_DEVICE, ZTMP5_DEVICE)
- !
- IF ( KRRI >= 1 ) THEN
- !$acc kernels
- PRRS(:,:,:,2) = PRRS(:,:,:,2) + 2. * (- ZTMP5_DEVICE(:,:,:)+ ZTMP3_DEVICE(:,:,:))*(1.0-PFRAC_ICE(:,:,:))
- PRRS(:,:,:,2) = PRRS(:,:,:,2) + 2. * (- ZTMP5_DEVICE(:,:,:)+ ZTMP3_DEVICE(:,:,:))*PFRAC_ICE(:,:,:)
- !$acc end kernels
- ELSE
- !$acc kernels
- PRRS(:,:,:,2) = PRRS(:,:,:,2) + 2. * (- ZTMP5_DEVICE(:,:,:) + ZTMP3_DEVICE(:,:,:))
- !$acc end kernels
- END IF
- ELSE
- !$acc kernels
- ZTMP4_DEVICE(:,:,:) = ZTMP8_DEVICE(:,:,:)*ZTMP2_DEVICE(:,:,:)
- !$acc end kernels
- CALL DXF_DEVICE(ZTMP4_DEVICE, ZTMP5_DEVICE)
- !$acc kernels
- ZTMP3_DEVICE(:,:,:) = ZTMP8_DEVICE(:,:,:)*ZFLX(:,:,:)
- !$acc end kernels
- CALL MXF_DEVICE( ZTMP3_DEVICE, ZTMP4_DEVICE )
- !$acc kernels
- ZFLXC(:,:,:) = ZFLXC(:,:,:) + 2.*ZTMP4_DEVICE(:,:,:)
- !$acc end kernels
- IF ( KRRI >= 1 ) THEN
- !$acc kernels
- PRRS(:,:,:,2) = PRRS(:,:,:,2) - 2. * ZTMP5_DEVICE(:,:,:)*(1.0-PFRAC_ICE(:,:,:))
- PRRS(:,:,:,4) = PRRS(:,:,:,4) - 2. * ZTMP5_DEVICE(:,:,:)*PFRAC_ICE(:,:,:)
- !$acc end kernels
- ELSE
- !$acc kernels
- PRRS(:,:,:,2) = PRRS(:,:,:,2) - 2. * ZTMP5_DEVICE(:,:,:)
- !$acc end kernels
- END IF
- END IF
- END IF
- !
- ! stores the horizontal <U Rnp>
- IF ( tpfile%lopened .AND. OTURB_FLX ) THEN
-!$acc update self(ZFLX)
- TZFIELD%CMNHNAME = 'UR_FLX'
- TZFIELD%CSTDNAME = ''
- TZFIELD%CLONGNAME = 'UR_FLX'
- TZFIELD%CUNITS = 'kg kg-1 m s-1'
- TZFIELD%CDIR = 'XY'
- TZFIELD%CCOMMENT = 'X_Y_Z_UR_FLX'
- TZFIELD%NGRID = 2
- TZFIELD%NTYPE = TYPEREAL
- TZFIELD%NDIMS = 3
- TZFIELD%LTIMEDEP = .TRUE.
- CALL IO_Field_write(TPFILE,TZFIELD,ZFLX(:,:,:))
- END IF
- !
- IF (KSPLT==1 .AND. LLES_CALL) THEN
- CALL SECOND_MNH(ZTIME1)
- !
-!$acc data copy(X_LES_SUBGRID_URt,X_LES_RES_ddxa_W_SBG_UaRt, &
-!$acc & X_LES_RES_ddxa_Thl_SBG_UaRt,X_LES_RES_ddxa_Rt_SBG_UaRt)
- !
- CALL MXF_DEVICE(ZFLX,ZTMP1_DEVICE)
- CALL LES_MEAN_SUBGRID( ZTMP1_DEVICE, X_LES_SUBGRID_URt )
- !
- CALL GX_W_UW_DEVICE(PWM,PDXX,PDZZ,PDZX,ZTMP1_DEVICE)
- CALL MZM_DEVICE(ZFLX,ZTMP2_DEVICE)
- !$acc kernels
- ZTMP3_DEVICE(:,:,:) = ZTMP1_DEVICE(:,:,:)*ZTMP2_DEVICE(:,:,:)
- !$acc end kernels
- CALL MXF_DEVICE(ZTMP3_DEVICE,ZTMP4_DEVICE)
- CALL MZF_DEVICE( ZTMP4_DEVICE, ZTMP3_DEVICE )
- CALL LES_MEAN_SUBGRID( ZTMP3_DEVICE, X_LES_RES_ddxa_W_SBG_UaRt , .TRUE. )
- !
- CALL GX_M_M_DEVICE(PTHLM,PDXX,PDZZ,PDZX,ZTMP1_DEVICE)
- CALL MXF_DEVICE(ZFLX,ZTMP2_DEVICE)
- !$acc kernels
- ZTMP3_DEVICE(:,:,:) = ZTMP1_DEVICE(:,:,:)*ZTMP2_DEVICE(:,:,:)
- !$acc end kernels
- CALL LES_MEAN_SUBGRID( ZTMP3_DEVICE, X_LES_RES_ddxa_Thl_SBG_UaRt , .TRUE. )
- !
- CALL GX_M_M_DEVICE(PRM(:,:,:,1),PDXX,PDZZ,PDZX,ZTMP1_DEVICE)
- CALL MXF_DEVICE(ZFLX,ZTMP2_DEVICE)
- !$acc kernels
- ZTMP3_DEVICE(:,:,:) = ZTMP1_DEVICE(:,:,:)*ZTMP2_DEVICE(:,:,:)
- !$acc end kernels
- CALL LES_MEAN_SUBGRID( ZTMP3_DEVICE, X_LES_RES_ddxa_Rt_SBG_UaRt , .TRUE. )
- !
-!$acc end data
- !
- CALL SECOND_MNH(ZTIME2)
- XTIME_LES = XTIME_LES + ZTIME2 - ZTIME1
- END IF
-!
- !
- IF (KRRL>0 .AND. KSPLT==1 .AND. LLES_CALL) THEN
- CALL SECOND_MNH(ZTIME1)
- !
- !$acc data copy(X_LES_SUBGRID_URc)
- !
- CALL MXF_DEVICE(ZFLXC,ZTMP1_DEVICE)
- CALL LES_MEAN_SUBGRID(ZTMP1_DEVICE, X_LES_SUBGRID_URc )
- !
- !$acc end data
- !
- CALL SECOND_MNH(ZTIME2)
- XTIME_LES = XTIME_LES + ZTIME2 - ZTIME1
- END IF
-!
-END IF
-#endif
-!
-!* 4. < U' TPV' >
-! -----------
-!
-!! to be tested later
-!!IF (KRR/=0) THEN
-!! ! here ZFLX= <U'Rnp'> and ZWORK= <U'Thetal'>
-!! !
-!! ZVPTU(:,:,:) = &
-!! ZWORK(:,:,:)*MXM(ETHETA(KRR,KRRI,PTHLT,PEXNREF,PRT,PLOCPT,PSRCM)) + &
-!! ZFLX(:,:,:)*MXM(EMOIST(KRR,KRRI,PTHLT,PEXNREF,PRT,PLOCPT,PSRCM))
-!! !
-!! ! stores the horizontal <U VPT>
-!! IF ( tpfile%lopened .AND. OTURB_FLX ) THEN
-!! TZFIELD%CMNHNAME = 'UVPT_FLX'
-!! TZFIELD%CSTDNAME = ''
-!! TZFIELD%CLONGNAME = 'UVPT_FLX'
-!! TZFIELD%CUNITS = 'K m s-1'
-!! TZFIELD%CDIR = 'XY'
-!! TZFIELD%CCOMMENT = 'X_Y_Z_UVPT_FLX'
-!! TZFIELD%NGRID = 2
-!! TZFIELD%NTYPE = TYPEREAL
-!! TZFIELD%NDIMS = 3
-!! TZFIELD%LTIMEDEP = .TRUE.
-!! CALL IO_Field_write(TPFILE,TZFIELD,ZVPTU)
-!! END IF
-!!!
-!!ELSE
-!! ZVPTU(:,:,:)=ZWORK(:,:,:)
-!!END IF
-!
-!
-!* 5. < V' THETA'l >
-! --------------
-!
-!
-IF (.NOT. L2D) THEN
-#ifndef MNH_OPENACC
- ZFLX(:,:,:) = -XCSHF * MYM( PK ) * GY_M_V(1,IKU,1,PTHLM,PDYY,PDZZ,PDZY)
- ZFLX(:,:,IKE+1) = ZFLX(:,:,IKE)
-ELSE
- ZFLX(:,:,:) = 0.
-END IF
-!
-!
-! Compute the flux at the first inner U-point with an uncentred vertical
-! gradient
-ZFLX(:,:,IKB:IKB) = -XCSHF * MYM( PK(:,:,IKB:IKB) ) * &
- ( DYM(PTHLM(:,:,IKB:IKB)) * PINV_PDYY(:,:,IKB:IKB) &
- -MYM( ZCOEFF(:,:,IKB+2:IKB+2)*PTHLM(:,:,IKB+2:IKB+2) &
- +ZCOEFF(:,:,IKB+1:IKB+1)*PTHLM(:,:,IKB+1:IKB+1) &
- +ZCOEFF(:,:,IKB :IKB )*PTHLM(:,:,IKB :IKB ) ) &
- *0.5* ( PDZY(:,:,IKB+1:IKB+1)+PDZY(:,:,IKB:IKB)) &
- * PINV_PDYY(:,:,IKB:IKB) )
-! extrapolates the flux under the ground so that the vertical average with
-! the IKB flux gives the ground value ( warning the tangential surface
-! flux has been set to 0 for the moment !! to be improved )
-ZFLX(:,:,IKB-1:IKB-1) = 2. * MYM( SPREAD( PSFTHM(:,:)* PDIRCOSYW(:,:), 3,1) ) &
- - ZFLX(:,:,IKB:IKB)
-!
-! Add this source to the Theta_l sources
-!
-IF (.NOT. L2D) THEN
- IF (.NOT. LFLAT) THEN
- PRTHLS(:,:,:) = PRTHLS(:,:,:) &
- - DYF( MYM(PRHODJ) * ZFLX(:,:,:) * PINV_PDYY(:,:,:) ) &
- + DZF( PMZM_PRHODJ *MYF(PDZY(:,:,:)*(MZM(ZFLX(:,:,:) * PINV_PDYY(:,:,:)))) * PINV_PDZZ(:,:,:) )
- ELSE
- PRTHLS(:,:,:) = PRTHLS(:,:,:) - DYF( MYM(PRHODJ) * ZFLX(:,:,:) * PINV_PDYY(:,:,:) )
- END IF
-END IF
-!
-! Compute the equivalent tendancy for Rc and Ri
-!
-!IF ( OSUBG_COND .AND. KRRL > 0 .AND. .NOT. L2D) THEN
-IF ( KRRL >= 1 .AND. .NOT. L2D) THEN
- IF (.NOT. LFLAT) THEN
- ZFLXC(:,:,:) = 2.*( MYF( MYM( PRHODJ(:,:,:)*PATHETA(:,:,:)*PSRCM(:,:,:) )*ZFLX(:,:,:) ) &
- +MZF( MZM( PRHODJ(:,:,:)*PATHETA(:,:,:)*PSRCM(:,:,:) )*MYF( &
- PDZY(:,:,:)*(MZM( ZFLX(:,:,:)*PINV_PDYY(:,:,:) )) ) )&
- )
- IF ( KRRI >= 1 ) THEN
- PRRS(:,:,:,2) = PRRS(:,:,:,2) + 2. * &
- (- DYF( MYM( PRHODJ(:,:,:)*PATHETA(:,:,:)*PSRCM(:,:,:) )*ZFLX(:,:,:)*PINV_PDYY(:,:,:) ) &
- + DZF( MZM( PRHODJ(:,:,:)*PATHETA(:,:,:)*PSRCM(:,:,:) )*MYF( PDZY(:,:,:)* &
- (MZM( ZFLX(:,:,:)*PINV_PDYY(:,:,:) )) )&
- *PINV_PDZZ(:,:,:) ) &
- )*(1.0-PFRAC_ICE(:,:,:))
- PRRS(:,:,:,4) = PRRS(:,:,:,4) + 2. * &
- (- DYF( MYM( PRHODJ(:,:,:)*PATHETA(:,:,:)*PSRCM(:,:,:) )*ZFLX(:,:,:)*PINV_PDYY(:,:,:) ) &
- + DZF( MZM( PRHODJ(:,:,:)*PATHETA(:,:,:)*PSRCM(:,:,:) )*MYF( PDZY(:,:,:)* &
- (MZM( ZFLX(:,:,:)*PINV_PDYY(:,:,:) )) )&
- *PINV_PDZZ(:,:,:) ) &
- )*PFRAC_ICE(:,:,:)
- ELSE
- PRRS(:,:,:,2) = PRRS(:,:,:,2) + 2. * &
- (- DYF( MYM( PRHODJ(:,:,:)*PATHETA(:,:,:)*PSRCM(:,:,:) )*ZFLX(:,:,:)*PINV_PDYY(:,:,:) ) &
- + DZF( MZM( PRHODJ(:,:,:)*PATHETA(:,:,:)*PSRCM(:,:,:) )*MYF( PDZY(:,:,:)* &
- (MZM( ZFLX(:,:,:)*PINV_PDYY(:,:,:) )) )&
- *PINV_PDZZ(:,:,:) ) &
- )
- END IF
- ELSE
- ZFLXC(:,:,:) = 2.*MYF( MYM( PRHODJ(:,:,:)*PATHETA(:,:,:)*PSRCM(:,:,:) )*ZFLX(:,:,:) )
- IF ( KRRI >= 1 ) THEN
- PRRS(:,:,:,2) = PRRS(:,:,:,2) - 2. * &
- DYF( MYM( PRHODJ(:,:,:)*PATHETA(:,:,:)*PSRCM(:,:,:) )*ZFLX(:,:,:)*PINV_PDYY(:,:,:) )*(1.0-PFRAC_ICE(:,:,:))
- PRRS(:,:,:,4) = PRRS(:,:,:,4) - 2. * &
- DYF( MYM( PRHODJ(:,:,:)*PATHETA(:,:,:)*PSRCM(:,:,:) )*ZFLX(:,:,:)*PINV_PDYY(:,:,:) )*PFRAC_ICE(:,:,:)
- ELSE
- PRRS(:,:,:,2) = PRRS(:,:,:,2) - 2. * &
- DYF( MYM( PRHODJ(:,:,:)*PATHETA(:,:,:)*PSRCM(:,:,:) )*ZFLX(:,:,:)*PINV_PDYY(:,:,:) )
- END IF
- END IF
-END IF
-!
-!! stores this flux in ZWORK to compute later <V' VPT'>
-!!ZWORK(:,:,:) = ZFLX(:,:,:)
-!
-! stores the horizontal <V THl>
-IF ( tpfile%lopened .AND. OTURB_FLX ) THEN
- TZFIELD%CMNHNAME = 'VTHL_FLX'
- TZFIELD%CSTDNAME = ''
- TZFIELD%CLONGNAME = 'VTHL_FLX'
- TZFIELD%CUNITS = 'K m s-1'
- TZFIELD%CDIR = 'XY'
- TZFIELD%CCOMMENT = 'X_Y_Z_VTHL_FLX'
- TZFIELD%NGRID = 3
- TZFIELD%NTYPE = TYPEREAL
- TZFIELD%NDIMS = 3
- TZFIELD%LTIMEDEP = .TRUE.
- CALL IO_Field_write(TPFILE,TZFIELD,ZFLX(:,:,:))
-END IF
-!
-IF (KSPLT==1 .AND. LLES_CALL) THEN
- CALL SECOND_MNH(ZTIME1)
- CALL LES_MEAN_SUBGRID( MYF(ZFLX), X_LES_SUBGRID_VThl )
- CALL LES_MEAN_SUBGRID( MZF(MYF(GY_W_VW(PWM,PDYY,PDZZ,PDZY)*MZM(ZFLX))),&
- X_LES_RES_ddxa_W_SBG_UaThl , .TRUE. )
- CALL LES_MEAN_SUBGRID( GY_M_M(PTHLM,PDYY,PDZZ,PDZY)*MYF(ZFLX),&
- X_LES_RES_ddxa_Thl_SBG_UaThl , .TRUE. )
- IF (KRR>=1) THEN
- CALL LES_MEAN_SUBGRID( GY_M_M(PRM(:,:,:,1),PDYY,PDZZ,PDZY)*MYF(ZFLX),&
- X_LES_RES_ddxa_Rt_SBG_UaThl , .TRUE. )
- END IF
- CALL SECOND_MNH(ZTIME2)
- XTIME_LES = XTIME_LES + ZTIME2 - ZTIME1
-END IF
-#else
- CALL MYM_DEVICE( PK, ZTMP1_DEVICE )
- CALL GY_M_V_DEVICE(1,IKU,1,PTHLM,PDYY,PDZZ,PDZY,ZTMP2_DEVICE)
- !$acc kernels
- !$acc_nv loop independent collapse(3)
- DO CONCURRENT ( JI=1:JIU,JJ=1:JJU,JK=1:JKU)
- ZFLX(JI,JJ,JK) = -XCSHF * ZTMP1_DEVICE(JI,JJ,JK) * ZTMP2_DEVICE(JI,JJ,JK)
- END DO
- ZFLX(:,:,IKE+1) = ZFLX(:,:,IKE)
- !$acc end kernels
-ELSE
- !$acc kernels
- ZFLX(:,:,:) = 0.
- !$acc end kernels
-END IF
-!
-!
-! Compute the flux at the first inner U-point with an uncentred vertical
-! gradient
-CALL MYM_DEVICE( PK(:,:,IKB:IKB), ZTMP1_DEVICE(:,:,1:1) )
-CALL DYM_DEVICE(PTHLM(:,:,IKB:IKB), ZTMP2_DEVICE(:,:,1:1) )
-!$acc kernels present_cr(ZTMP3_DEVICE)
-ZTMP3_DEVICE(:,:,1) = ZCOEFF(:,:,IKB+2)*PTHLM(:,:,IKB+2) &
- +ZCOEFF(:,:,IKB+1)*PTHLM(:,:,IKB+1) &
- +ZCOEFF(:,:,IKB )*PTHLM(:,:,IKB )
-!$acc end kernels
-CALL MYM_DEVICE( ZTMP3_DEVICE(:,:,1:1), ZTMP4_DEVICE(:,:,1:1) )
-!$acc kernels present_cr(ZFLX)
-ZFLX(:,:,IKB) = -XCSHF * ZTMP1_DEVICE(:,:,1) * &
- ( ZTMP2_DEVICE(:,:,1) * PINV_PDYY(:,:,IKB) &
- - ZTMP4_DEVICE(:,:,1) &
- *0.5* ( PDZY(:,:,IKB+1)+PDZY(:,:,IKB)) &
- * PINV_PDYY(:,:,IKB) )
-!$acc end kernels
-! extrapolates the flux under the ground so that the vertical average with
-! the IKB flux gives the ground value ( warning the tangential surface
-! flux has been set to 0 for the moment !! to be improved )
-!$acc kernels
-ZTMP1_DEVICE(:,:,1) = PSFTHM(:,:)* PDIRCOSYW(:,:)
-!$acc end kernels
-CALL MYM_DEVICE( ZTMP1_DEVICE(:,:,1:1), ZTMP2_DEVICE(:,:,1:1) )
-!$acc kernels present_cr(ZFLX)
-ZFLX(:,:,IKB-1) = 2. * ZTMP2_DEVICE(:,:,1) - ZFLX(:,:,IKB)
-!$acc end kernels
-!
-! Add this source to the Theta_l sources
-!
-IF (.NOT. L2D) THEN
- IF (.NOT. LFLAT) THEN
- CALL MYM_DEVICE(PRHODJ, ZTMP1_DEVICE)
- !$acc kernels
- !$acc_nv loop independent collapse(3)
- DO CONCURRENT ( JI=1:JIU,JJ=1:JJU,JK=1:JKU)
- ZTMP2_DEVICE(JI,JJ,JK) = ZTMP1_DEVICE(JI,JJ,JK) * ZFLX(JI,JJ,JK) * PINV_PDYY(JI,JJ,JK)
- END DO
- !$acc end kernels
- CALL DYF_DEVICE( ZTMP2_DEVICE, ZTMP3_DEVICE )
- !$acc kernels
- !$acc_nv loop independent collapse(3)
- DO CONCURRENT ( JI=1:JIU,JJ=1:JJU,JK=1:JKU)
- ZTMP1_DEVICE(JI,JJ,JK) = ZFLX(JI,JJ,JK) * PINV_PDYY(JI,JJ,JK)
- END DO
- !$acc end kernels
- CALL MZM_DEVICE(ZTMP1_DEVICE, ZTMP2_DEVICE)
- !$acc kernels
- !$acc_nv loop independent collapse(3)
- DO CONCURRENT ( JI=1:JIU,JJ=1:JJU,JK=1:JKU)
- ZTMP1_DEVICE(JI,JJ,JK) = PDZY(JI,JJ,JK)*ZTMP2_DEVICE(JI,JJ,JK)
- END DO
- !$acc end kernels
- CALL MYF_DEVICE(ZTMP1_DEVICE, ZTMP2_DEVICE)
- !$acc kernels
- !$acc_nv loop independent collapse(3)
- DO CONCURRENT ( JI=1:JIU,JJ=1:JJU,JK=1:JKU)
- ZTMP1_DEVICE(JI,JJ,JK) = PMZM_PRHODJ(JI,JJ,JK) * ZTMP2_DEVICE(JI,JJ,JK) * PINV_PDZZ(JI,JJ,JK)
- END DO
- !$acc end kernels
- CALL DZF_DEVICE( ZTMP1_DEVICE, ZTMP2_DEVICE )
- !$acc kernels
- !$acc_nv loop independent collapse(3)
- DO CONCURRENT ( JI=1:JIU,JJ=1:JJU,JK=1:JKU)
- PRTHLS(JI,JJ,JK) = PRTHLS(JI,JJ,JK) - ZTMP3_DEVICE(JI,JJ,JK) + ZTMP2_DEVICE(JI,JJ,JK)
- END DO
- !$acc end kernels
- ELSE
- CALL MYM_DEVICE(PRHODJ, ZTMP1_DEVICE)
- !$acc kernels
- ZTMP2_DEVICE(:,:,:) = ZTMP1_DEVICE(:,:,:) * ZFLX(:,:,:) * PINV_PDYY(:,:,:)
- !$acc end kernels
- CALL DYF_DEVICE( ZTMP2_DEVICE, ZTMP3_DEVICE )
- !$acc kernels
- PRTHLS(:,:,:) = PRTHLS(:,:,:) - ZTMP3_DEVICE(:,:,:)
- !$acc end kernels
- END IF
-END IF
-!
-! Compute the equivalent tendancy for Rc and Ri
-!
-!IF ( OSUBG_COND .AND. KRRL > 0 .AND. .NOT. L2D) THEN
-IF ( KRRL >= 1 .AND. .NOT. L2D) THEN
- IF (.NOT. LFLAT) THEN
- !$acc kernels
- ZTMP1_DEVICE(:,:,:) = PRHODJ(:,:,:)*PATHETA(:,:,:)*PSRCM(:,:,:)
- !$acc end kernels
- CALL MYM_DEVICE( ZTMP1_DEVICE, ZTMP2_DEVICE )
- !$acc kernels present_cr(ZTMP4_DEVICE)
- ZTMP4_DEVICE(:,:,:) = ZTMP2_DEVICE(:,:,:)*ZFLX(:,:,:)
- !$acc end kernels
- CALL MYF_DEVICE( ZTMP4_DEVICE, ZTMP3_DEVICE )
- !$acc kernels present_cr(ZTMP4_DEVICE)
- ZTMP4_DEVICE(:,:,:) = ZFLX(:,:,:)*PINV_PDYY(:,:,:)
- !$acc end kernels
- CALL MZM_DEVICE( ZTMP4_DEVICE, ZTMP5_DEVICE )
- !$acc kernels present_cr(ZTMP4_DEVICE)
- ZTMP4_DEVICE(:,:,:) = PDZY(:,:,:)*ZTMP5_DEVICE(:,:,:)
- !$acc end kernels
- CALL MYF_DEVICE( ZTMP4_DEVICE, ZTMP5_DEVICE)
- CALL MZM_DEVICE( ZTMP1_DEVICE, ZTMP4_DEVICE )
- !$acc kernels present_cr(ZTMP6_DEVICE)
- ZTMP6_DEVICE(:,:,:) = ZTMP4_DEVICE(:,:,:)*ZTMP5_DEVICE(:,:,:)
- !$acc end kernels
- CALL MZF_DEVICE( ZTMP6_DEVICE, ZTMP4_DEVICE )
- !$acc kernels present_cr(ZFLXC)
- ZFLXC(:,:,:) = 2.*( ZTMP3_DEVICE(:,:,:) + ZTMP4_DEVICE(:,:,:) )
- !$acc end kernels
- IF ( KRRI >= 1 ) THEN
- !$acc kernels present_cr(ZTMP3_DEVICE)
- ZTMP3_DEVICE(:,:,:) = ZTMP2_DEVICE(:,:,:)*ZFLX(:,:,:)*PINV_PDYY(:,:,:)
- !$acc end kernels
- CALL DYF_DEVICE( ZTMP3_DEVICE, ZTMP4_DEVICE )
- !$acc kernels present_cr(ZTMP3_DEVICE)
- ZTMP3_DEVICE(:,:,:) = ZTMP6_DEVICE(:,:,:)*PINV_PDZZ(:,:,:)
- !$acc end kernels
- CALL DZF_DEVICE( ZTMP3_DEVICE, ZTMP5_DEVICE )
- !$acc kernels
- PRRS(:,:,:,2) = PRRS(:,:,:,2) + 2. * (- ZTMP4_DEVICE(:,:,:) + ZTMP5_DEVICE(:,:,:) )*(1.0-PFRAC_ICE(:,:,:))
- PRRS(:,:,:,4) = PRRS(:,:,:,4) + 2. * (- ZTMP4_DEVICE(:,:,:) + ZTMP5_DEVICE(:,:,:) )*PFRAC_ICE(:,:,:)
- !$acc end kernels
- ELSE
- !$acc kernels
- ZTMP3_DEVICE(:,:,:) = ZTMP2_DEVICE(:,:,:)*ZFLX(:,:,:)*PINV_PDYY(:,:,:)
- !$acc end kernels
- CALL DYF_DEVICE( ZTMP3_DEVICE, ZTMP4_DEVICE )
- !$acc kernels
- ZTMP3_DEVICE(:,:,:) = ZTMP6_DEVICE(:,:,:)*PINV_PDZZ(:,:,:)
- !$acc end kernels
- CALL DZF_DEVICE( ZTMP3_DEVICE, ZTMP5_DEVICE )
- !$acc kernels
- PRRS(:,:,:,2) = PRRS(:,:,:,2) + 2. * (- ZTMP4_DEVICE(:,:,:) + ZTMP5_DEVICE(:,:,:) )
- !$acc end kernels
- END IF
- ELSE
- !$acc kernels
- ZTMP1_DEVICE(:,:,:) = PRHODJ(:,:,:)*PATHETA(:,:,:)*PSRCM(:,:,:)
- !$acc end kernels
- CALL MYM_DEVICE( ZTMP1_DEVICE, ZTMP2_DEVICE )
- !$acc kernels
- ZTMP1_DEVICE(:,:,:) = ZTMP2_DEVICE(:,:,:)*ZFLX(:,:,:)
- !$acc end kernels
- CALL MYF_DEVICE( ZTMP1_DEVICE, ZTMP3_DEVICE )
- !$acc kernels
- ZFLXC(:,:,:) = 2.*ZTMP3_DEVICE(:,:,:)
- !$acc end kernels
- !
- !$acc kernels
- ZTMP1_DEVICE(:,:,:) = ZTMP2_DEVICE(:,:,:)*ZFLX(:,:,:)*PINV_PDYY(:,:,:)
- !$acc end kernels
- CALL DYF_DEVICE( ZTMP1_DEVICE, ZTMP2_DEVICE )
- IF ( KRRI >= 1 ) THEN
- !$acc kernels
- PRRS(:,:,:,2) = PRRS(:,:,:,2) - 2. * ZTMP2_DEVICE(:,:,:)*(1.0-PFRAC_ICE(:,:,:))
- PRRS(:,:,:,4) = PRRS(:,:,:,4) - 2. * ZTMP2_DEVICE(:,:,:)*PFRAC_ICE(:,:,:)
- !$acc end kernels
- ELSE
- !$acc kernels
- PRRS(:,:,:,2) = PRRS(:,:,:,2) - 2. * ZTMP2_DEVICE(:,:,:)
- !$acc end kernels
- END IF
- END IF
-END IF
-!! stores this flux in ZWORK to compute later <V' VPT'>
-!!ZWORK(:,:,:) = ZFLX(:,:,:)
-!
-! stores the horizontal <V THl>
-IF ( tpfile%lopened .AND. OTURB_FLX ) THEN
-!$acc update self(ZFLX)
- TZFIELD%CMNHNAME = 'VTHL_FLX'
- TZFIELD%CSTDNAME = ''
- TZFIELD%CLONGNAME = 'VTHL_FLX'
- TZFIELD%CUNITS = 'K m s-1'
- TZFIELD%CDIR = 'XY'
- TZFIELD%CCOMMENT = 'X_Y_Z_VTHL_FLX'
- TZFIELD%NGRID = 3
- TZFIELD%NTYPE = TYPEREAL
- TZFIELD%NDIMS = 3
- TZFIELD%LTIMEDEP = .TRUE.
- CALL IO_Field_write(TPFILE,TZFIELD,ZFLX(:,:,:))
-END IF
-!
-IF (KSPLT==1 .AND. LLES_CALL) THEN
- CALL SECOND_MNH(ZTIME1)
- !
-!$acc data copy(X_LES_SUBGRID_VThl,X_LES_RES_ddxa_W_SBG_UaThl,X_LES_RES_ddxa_Thl_SBG_UaThl)
- !
- CALL MYF_DEVICE(ZFLX, ZTMP1_DEVICE)
- CALL LES_MEAN_SUBGRID( ZTMP1_DEVICE, X_LES_SUBGRID_VThl )
- !
- CALL GY_W_VW_DEVICE(PWM,PDYY,PDZZ,PDZY,ZTMP1_DEVICE)
- CALL MZM_DEVICE(ZFLX,ZTMP2_DEVICE)
- !$acc kernels
- ZTMP3_DEVICE(:,:,:) = ZTMP1_DEVICE(:,:,:) * ZTMP2_DEVICE(:,:,:)
- !$acc end kernels
- CALL MYF_DEVICE(ZTMP3_DEVICE, ZTMP4_DEVICE)
- CALL MZF_DEVICE( ZTMP4_DEVICE, ZTMP1_DEVICE )
- CALL LES_MEAN_SUBGRID( ZTMP1_DEVICE, X_LES_RES_ddxa_W_SBG_UaThl , .TRUE. )
- !
- CALL GY_M_M_DEVICE(PTHLM,PDYY,PDZZ,PDZY,ZTMP1_DEVICE)
- CALL MYF_DEVICE(ZFLX,ZTMP2_DEVICE)
- !$acc kernels
- ZTMP3_DEVICE(:,:,:) = ZTMP1_DEVICE(:,:,:) * ZTMP2_DEVICE(:,:,:)
- !$acc end kernels
- CALL LES_MEAN_SUBGRID( ZTMP3_DEVICE, X_LES_RES_ddxa_Thl_SBG_UaThl , .TRUE. )
- !
-!$acc end data
- !
- IF (KRR>=1) THEN
-!$acc data copy(X_LES_RES_ddxa_Rt_SBG_UaThl)
- CALL GY_M_M_DEVICE(PRM(:,:,:,1),PDYY,PDZZ,PDZY,ZTMP1_DEVICE)
- CALL MYF_DEVICE(ZFLX,ZTMP2_DEVICE)
- !$acc kernels
- ZTMP3_DEVICE(:,:,:) = ZTMP1_DEVICE(:,:,:) * ZTMP2_DEVICE(:,:,:)
- !$acc end kernels
- CALL LES_MEAN_SUBGRID( ZTMP3_DEVICE,X_LES_RES_ddxa_Rt_SBG_UaThl , .TRUE. )
-!$acc end data
- END IF
- !
- CALL SECOND_MNH(ZTIME2)
- XTIME_LES = XTIME_LES + ZTIME2 - ZTIME1
-END IF
-#endif
-!
-!
-!* 6. < V' R'np >
-! -----------
-!
-#ifndef MNH_OPENACC
-IF (KRR/=0) THEN
- !
- IF (.NOT. L2D) THEN
- ZFLX(:,:,:) = -XCHF * MYM( PK ) * GY_M_V(1,IKU,1,PRM(:,:,:,1),PDYY,PDZZ,PDZY)
- ZFLX(:,:,IKE+1) = ZFLX(:,:,IKE)
- ELSE
- ZFLX(:,:,:) = 0.
- END IF
-!
-! Compute the flux at the first inner U-point with an uncentred vertical
-! gradient
- ZFLX(:,:,IKB:IKB) = -XCHF * MYM( PK(:,:,IKB:IKB) ) * &
- ( DYM(PRM(:,:,IKB:IKB,1)) * PINV_PDYY(:,:,IKB:IKB) &
- -MYM( ZCOEFF(:,:,IKB+2:IKB+2)*PRM(:,:,IKB+2:IKB+2,1) &
- +ZCOEFF(:,:,IKB+1:IKB+1)*PRM(:,:,IKB+1:IKB+1,1) &
- +ZCOEFF(:,:,IKB :IKB )*PRM(:,:,IKB :IKB ,1) ) &
- *0.5* ( PDZY(:,:,IKB+1:IKB+1)+PDZY(:,:,IKB:IKB)) &
- * PINV_PDYY(:,:,IKB:IKB) )
-! extrapolates the flux under the ground so that the vertical average with
-! the IKB flux gives the ground value ( warning the tangential surface
-! flux has been set to 0 for the moment !! to be improved )
- ZFLX(:,:,IKB-1:IKB-1) = 2. * MYM( SPREAD( PSFRM(:,:)* PDIRCOSYW(:,:), 3,1) ) &
- - ZFLX(:,:,IKB:IKB)
- !
- ! Add this source to the conservative mixing ratio sources
- !
- IF (.NOT. L2D) THEN
- IF (.NOT. LFLAT) THEN
- PRRS(:,:,:,1) = PRRS(:,:,:,1) &
- - DYF( MYM(PRHODJ) * ZFLX(:,:,:) * PINV_PDYY(:,:,:) ) &
- + DZF( PMZM_PRHODJ(:,:,:) *MYF(PDZY(:,:,:)* &
- (MZM(ZFLX(:,:,:) * PINV_PDYY(:,:,:)))) * PINV_PDZZ(:,:,:) )
- ELSE
- PRRS(:,:,:,1) = PRRS(:,:,:,1) - DYF( MYM(PRHODJ) * ZFLX(:,:,:) * PINV_PDYY(:,:,:) )
- END IF
- END IF
- !
- ! Compute the equivalent tendancy for Rc and Ri
- !
- IF ( KRRL >= 1 .AND. .NOT. L2D) THEN ! Sub-grid condensation
- IF (.NOT. LFLAT) THEN
- ZFLXC(:,:,:) = ZFLXC(:,:,:) &
- + 2.*( MXF( MYM( PRHODJ(:,:,:)*PAMOIST(:,:,:)*PSRCM(:,:,:) )*ZFLX(:,:,:) ) &
- + MZF( MZM( PRHODJ(:,:,:)*PAMOIST(:,:,:)*PSRCM(:,:,:) )*MYF( &
- PDZY(:,:,:)*(MZM( ZFLX(:,:,:)*PINV_PDYY(:,:,:) )) ) ) &
- )
- IF ( KRRI >= 1 ) THEN
- PRRS(:,:,:,2) = PRRS(:,:,:,2) + 2. * &
- (- DYF( MYM( PRHODJ(:,:,:)*PAMOIST(:,:,:)*PSRCM(:,:,:) )*ZFLX(:,:,:)*PINV_PDYY ) &
- + DZF( MZM( PRHODJ(:,:,:)*PAMOIST(:,:,:)*PSRCM(:,:,:) )*MYF( PDZY(:,:,:)* &
- (MZM( ZFLX(:,:,:)*PINV_PDYY(:,:,:) )) ) &
- * PINV_PDZZ(:,:,:) ) &
- )*(1.0-PFRAC_ICE(:,:,:))
- PRRS(:,:,:,4) = PRRS(:,:,:,4) + 2. * &
- (- DYF( MYM( PRHODJ(:,:,:)*PAMOIST(:,:,:)*PSRCM(:,:,:) )*ZFLX(:,:,:)*PINV_PDYY ) &
- + DZF( MZM( PRHODJ(:,:,:)*PAMOIST(:,:,:)*PSRCM(:,:,:) )*MYF( PDZY(:,:,:)* &
- (MZM( ZFLX(:,:,:)*PINV_PDYY(:,:,:) )) ) &
- * PINV_PDZZ(:,:,:) ) &
- )*PFRAC_ICE(:,:,:)
- ELSE
- PRRS(:,:,:,2) = PRRS(:,:,:,2) + 2. * &
- (- DYF( MYM( PRHODJ(:,:,:)*PAMOIST(:,:,:)*PSRCM(:,:,:) )*ZFLX(:,:,:)*PINV_PDYY ) &
- + DZF( MZM( PRHODJ(:,:,:)*PAMOIST(:,:,:)*PSRCM(:,:,:) )*MYF( PDZY(:,:,:)* &
- (MZM( ZFLX(:,:,:)*PINV_PDYY(:,:,:) )) ) &
- * PINV_PDZZ(:,:,:) ) &
- )
- END IF
- ELSE
- ZFLXC(:,:,:) = ZFLXC(:,:,:) + 2.*MXF( MYM( PRHODJ(:,:,:)*PAMOIST(:,:,:)*PSRCM(:,:,:) )*ZFLX )
- IF ( KRRI >= 1 ) THEN
- PRRS(:,:,:,2) = PRRS(:,:,:,2) - 2. * &
- DYF( MYM( PRHODJ(:,:,:)*PAMOIST(:,:,:)*PSRCM(:,:,:) )*ZFLX(:,:,:)*PINV_PDYY )*(1.0-PFRAC_ICE(:,:,:))
- PRRS(:,:,:,4) = PRRS(:,:,:,4) - 2. * &
- DYF( MYM( PRHODJ(:,:,:)*PAMOIST(:,:,:)*PSRCM(:,:,:) )*ZFLX(:,:,:)*PINV_PDYY )*PFRAC_ICE(:,:,:)
- ELSE
- PRRS(:,:,:,2) = PRRS(:,:,:,2) - 2. * &
- DYF( MYM( PRHODJ(:,:,:)*PAMOIST(:,:,:)*PSRCM(:,:,:) )*ZFLX(:,:,:)*PINV_PDYY )
- END IF
- END IF
- END IF
- !
- ! stores the horizontal <V Rnp>
- IF ( tpfile%lopened .AND. OTURB_FLX ) THEN
- TZFIELD%CMNHNAME = 'VR_FLX'
- TZFIELD%CSTDNAME = ''
- TZFIELD%CLONGNAME = 'VR_FLX'
- TZFIELD%CUNITS = 'kg kg-1 m s-1'
- TZFIELD%CDIR = 'XY'
- TZFIELD%CCOMMENT = 'X_Y_Z_VR_FLX'
- TZFIELD%NGRID = 3
- TZFIELD%NTYPE = TYPEREAL
- TZFIELD%NDIMS = 3
- TZFIELD%LTIMEDEP = .TRUE.
- CALL IO_Field_write(TPFILE,TZFIELD,ZFLX(:,:,:))
- END IF
- !
- IF (KSPLT==1 .AND. LLES_CALL) THEN
- CALL SECOND_MNH(ZTIME1)
- CALL LES_MEAN_SUBGRID( MYF(ZFLX), X_LES_SUBGRID_VRt )
- CALL LES_MEAN_SUBGRID( MZF(MYF(GY_W_VW(PWM,PDYY,PDZZ,PDZY)*MZM(ZFLX))),&
- X_LES_RES_ddxa_W_SBG_UaRt , .TRUE. )
- CALL LES_MEAN_SUBGRID( GY_M_M(PTHLM,PDYY,PDZZ,PDZY)*MYF(ZFLX), &
- X_LES_RES_ddxa_Thl_SBG_UaRt , .TRUE. )
- CALL LES_MEAN_SUBGRID( GY_M_M(PRM(:,:,:,1),PDYY,PDZZ,PDZY)*MYF(ZFLX), &
- X_LES_RES_ddxa_Rt_SBG_UaRt , .TRUE. )
- CALL SECOND_MNH(ZTIME2)
- XTIME_LES = XTIME_LES + ZTIME2 - ZTIME1
- END IF
-!
- !
- IF (KRRL>0 .AND. KSPLT==1 .AND. LLES_CALL) THEN
- CALL SECOND_MNH(ZTIME1)
- CALL LES_MEAN_SUBGRID(MYF(ZFLXC), X_LES_SUBGRID_VRc )
- CALL SECOND_MNH(ZTIME2)
- XTIME_LES = XTIME_LES + ZTIME2 - ZTIME1
- END IF
- !
-END IF
-#else
-IF (KRR/=0) THEN
- !
- IF (.NOT. L2D) THEN
- CALL MYM_DEVICE( PK, ZTMP1_DEVICE )
- CALL GY_M_V_DEVICE(1,IKU,1,PRM(:,:,:,1),PDYY,PDZZ,PDZY, ZTMP2_DEVICE)
- !$acc kernels
- !$acc_nv loop independent collapse(3)
- DO CONCURRENT ( JI=1:JIU,JJ=1:JJU,JK=1:JKU)
- ZFLX(JI,JJ,JK) = -XCHF * ZTMP1_DEVICE(JI,JJ,JK) * ZTMP2_DEVICE(JI,JJ,JK)
- END DO !CONCURRENT
- ZFLX(:,:,IKE+1) = ZFLX(:,:,IKE)
- !$acc end kernels
- ELSE
- !$acc kernels
- ZFLX(:,:,:) = 0.
- !$acc end kernels
- END IF
-!
-! Compute the flux at the first inner U-point with an uncentred vertical
-! gradient
- CALL MYM_DEVICE( PK(:,:,IKB:IKB), ZTMP1_DEVICE(:,:,1:1) )
- CALL DYM_DEVICE(PRM(:,:,IKB:IKB,1), ZTMP2_DEVICE(:,:,1:1))
- !$acc kernels present_cr(ZTMP3_DEVICE)
- ZTMP3_DEVICE(:,:,1) = ZCOEFF(:,:,IKB+2)*PRM(:,:,IKB+2,1) &
- +ZCOEFF(:,:,IKB+1)*PRM(:,:,IKB+1,1) &
- +ZCOEFF(:,:,IKB )*PRM(:,:,IKB ,1)
- !$acc end kernels
- CALL MYM_DEVICE( ZTMP3_DEVICE(:,:,1:1), ZTMP4_DEVICE(:,:,1:1) )
- !$acc kernels present_cr(ZFLX)
- ZFLX(:,:,IKB) = -XCHF * ZTMP1_DEVICE(:,:,1) * &
- ( ZTMP2_DEVICE(:,:,1) * PINV_PDYY(:,:,IKB) &
- - ZTMP4_DEVICE(:,:,1) &
- *0.5* ( PDZY(:,:,IKB+1)+PDZY(:,:,IKB)) &
- * PINV_PDYY(:,:,IKB) )
- !$acc end kernels
-! extrapolates the flux under the ground so that the vertical average with
-! the IKB flux gives the ground value ( warning the tangential surface
-! flux has been set to 0 for the moment !! to be improved )
- !$acc kernels
- ZTMP1_DEVICE(:,:,1) = PSFRM(:,:)* PDIRCOSYW(:,:)
- !$acc end kernels
- CALL MYM_DEVICE( ZTMP1_DEVICE(:,:,1:1), ZTMP2_DEVICE(:,:,1:1) )
- !$acc kernels present_cr(ZFLX)
- ZFLX(:,:,IKB-1) = 2. * ZTMP2_DEVICE(:,:,1) - ZFLX(:,:,IKB)
- !$acc end kernels
- !
- ! Add this source to the conservative mixing ratio sources
- !
- IF (.NOT. L2D) THEN
- IF (.NOT. LFLAT) THEN
- CALL MYM_DEVICE(PRHODJ, ZTMP1_DEVICE)
- !$acc kernels
- !$acc_nv loop independent collapse(3)
- DO CONCURRENT ( JI=1:JIU,JJ=1:JJU,JK=1:JKU)
- ZTMP2_DEVICE(JI,JJ,JK) = ZTMP1_DEVICE(JI,JJ,JK) * ZFLX(JI,JJ,JK) * PINV_PDYY(JI,JJ,JK)
- END DO
- !$acc end kernels
- CALL DYF_DEVICE( ZTMP2_DEVICE, ZTMP3_DEVICE )
- !
- !$acc kernels
- !$acc_nv loop independent collapse(3)
- DO CONCURRENT ( JI=1:JIU,JJ=1:JJU,JK=1:JKU)
- ZTMP1_DEVICE(JI,JJ,JK) = ZFLX(JI,JJ,JK) * PINV_PDYY(JI,JJ,JK)
- END DO
- !$acc end kernels
- CALL MZM_DEVICE(ZTMP1_DEVICE,ZTMP2_DEVICE)
- !$acc kernels
- !$acc_nv loop independent collapse(3)
- DO CONCURRENT ( JI=1:JIU,JJ=1:JJU,JK=1:JKU)
- ZTMP1_DEVICE(JI,JJ,JK) = PDZY(JI,JJ,JK)*ZTMP2_DEVICE(JI,JJ,JK)
- END DO
- !$acc end kernels
- CALL MYF_DEVICE(ZTMP1_DEVICE,ZTMP2_DEVICE)
- !$acc kernels
- !$acc_nv loop independent collapse(3)
- DO CONCURRENT ( JI=1:JIU,JJ=1:JJU,JK=1:JKU)
- ZTMP1_DEVICE(JI,JJ,JK) = PMZM_PRHODJ(JI,JJ,JK) * ZTMP2_DEVICE(JI,JJ,JK) * PINV_PDZZ(JI,JJ,JK)
- END DO
- !$acc end kernels
- CALL DZF_DEVICE( ZTMP1_DEVICE, ZTMP2_DEVICE )
- !
- !$acc kernels
- !$acc_nv loop independent collapse(3)
- DO CONCURRENT ( JI=1:JIU,JJ=1:JJU,JK=1:JKU)
- PRRS(JI,JJ,JK,1) = PRRS(JI,JJ,JK,1) - ZTMP3_DEVICE(JI,JJ,JK) + ZTMP2_DEVICE(JI,JJ,JK)
- END DO
- !$acc end kernels
- ELSE
- CALL MYM_DEVICE(PRHODJ, ZTMP1_DEVICE)
- !$acc kernels
- ZTMP2_DEVICE(:,:,:) = ZTMP1_DEVICE(:,:,:) * ZFLX(:,:,:) * PINV_PDYY(:,:,:)
- !$acc end kernels
- CALL DYF_DEVICE( ZTMP2_DEVICE, ZTMP3_DEVICE )
- !$acc kernels
- PRRS(:,:,:,1) = PRRS(:,:,:,1) - ZTMP3_DEVICE(:,:,:)
- !$acc end kernels
- END IF
- END IF
- !
- ! Compute the equivalent tendancy for Rc and Ri
- !
- IF ( KRRL >= 1 .AND. .NOT. L2D) THEN ! Sub-grid condensation
- IF (.NOT. LFLAT) THEN
- !$acc kernels
- ZTMP1_DEVICE(:,:,:) = PRHODJ(:,:,:)*PAMOIST(:,:,:)*PSRCM(:,:,:)
- !$acc end kernels
- CALL MYM_DEVICE( ZTMP1_DEVICE, ZTMP2_DEVICE )
- !$acc kernels present_cr(ZTMP3_DEVICE)
- ZTMP3_DEVICE(:,:,:) = ZTMP2_DEVICE(:,:,:)*ZFLX(:,:,:)
- !$acc end kernels
- CALL MXF_DEVICE( ZTMP3_DEVICE, ZTMP4_DEVICE )
- CALL MZM_DEVICE( ZTMP1_DEVICE, ZTMP5_DEVICE )
- !$acc kernels present_cr(ZTMP1_DEVICE)
- ZTMP1_DEVICE(:,:,:) = ZFLX(:,:,:)*PINV_PDYY(:,:,:)
- !$acc end kernels
- CALL MZM_DEVICE( ZTMP1_DEVICE, ZTMP2_DEVICE )
- !$acc kernels present_cr(ZTMP1_DEVICE)
- ZTMP1_DEVICE(:,:,:) = PDZY(:,:,:)*ZTMP2_DEVICE(:,:,:)
- !$acc end kernels
- CALL MYF_DEVICE( ZTMP1_DEVICE, ZTMP2_DEVICE )
- !$acc kernels present_cr(ZTMP1_DEVICE)
- ZTMP1_DEVICE(:,:,:) = ZTMP5_DEVICE(:,:,:)*ZTMP2_DEVICE(:,:,:)
- !$acc end kernels
- CALL MZF_DEVICE( ZTMP1_DEVICE, ZTMP2_DEVICE )
- !$acc kernels present_cr(ZFLXC)
- ZFLXC(:,:,:) = ZFLXC(:,:,:) + 2.*( ZTMP4_DEVICE(:,:,:) + ZTMP2_DEVICE(:,:,:) )
- !$acc end kernels
- IF ( KRRI >= 1 ) THEN
- !$acc kernels present_cr(ZTMP2_DEVICE)
- ZTMP2_DEVICE(:,:,:) = ZTMP3_DEVICE(:,:,:) * PINV_PDYY(:,:,:)
- !$acc end kernels
- CALL DYF_DEVICE( ZTMP2_DEVICE, ZTMP3_DEVICE )
- !$acc kernels present_cr(ZTMP2_DEVICE)
- ZTMP2_DEVICE(:,:,:) = ZTMP1_DEVICE(:,:,:)* PINV_PDZZ(:,:,:)
- !$acc end kernels
- CALL DZF_DEVICE( ZTMP2_DEVICE, ZTMP4_DEVICE )
- !$acc kernels
- PRRS(:,:,:,2) = PRRS(:,:,:,2) + 2. * (- ZTMP3_DEVICE(:,:,:) + ZTMP4_DEVICE(:,:,:) )*(1.0-PFRAC_ICE(:,:,:))
- PRRS(:,:,:,4) = PRRS(:,:,:,4) + 2. * (- ZTMP3_DEVICE(:,:,:) + ZTMP4_DEVICE(:,:,:) )*PFRAC_ICE(:,:,:)
- !$acc end kernels
- ELSE
- !$acc kernels
- PRRS(:,:,:,2) = PRRS(:,:,:,2) + 2. * (- ZTMP3_DEVICE(:,:,:) + ZTMP4_DEVICE(:,:,:) )
- !$acc end kernels
- END IF
- ELSE
- !$acc kernels
- ZTMP1_DEVICE(:,:,:) = PRHODJ(:,:,:)*PAMOIST(:,:,:)*PSRCM(:,:,:)
- !$acc end kernels
- CALL MYM_DEVICE( ZTMP1_DEVICE, ZTMP2_DEVICE )
- !$acc kernels
- ZTMP3_DEVICE(:,:,:) = ZTMP2_DEVICE(:,:,:)*ZFLX(:,:,:)
- !$acc end kernels
- CALL MXF_DEVICE( ZTMP3_DEVICE, ZTMP4_DEVICE )
- !$acc kernels
- ZFLXC(:,:,:) = ZFLXC(:,:,:) + 2.*ZTMP4_DEVICE(:,:,:)
- !$acc end kernels
- !
- !$acc kernels
- ZTMP1_DEVICE(:,:,:) = ZTMP3_DEVICE(:,:,:)* PINV_PDYY(:,:,:)
- !$acc end kernels
- CALL DYF_DEVICE( ZTMP1_DEVICE, ZTMP2_DEVICE )
- IF ( KRRI >= 1 ) THEN
- !$acc kernels
- PRRS(:,:,:,2) = PRRS(:,:,:,2) - 2. * ZTMP2_DEVICE(:,:,:)*(1.0-PFRAC_ICE(:,:,:))
- PRRS(:,:,:,4) = PRRS(:,:,:,4) - 2. * ZTMP2_DEVICE(:,:,:)*PFRAC_ICE(:,:,:)
- !$acc end kernels
- ELSE
- !$acc kernels
- PRRS(:,:,:,2) = PRRS(:,:,:,2) - 2. * ZTMP2_DEVICE(:,:,:)
- !$acc end kernels
- END IF
- END IF
- END IF
- !
- ! stores the horizontal <V Rnp>
- IF ( tpfile%lopened .AND. OTURB_FLX ) THEN
- !$acc update self(ZFLX)
- TZFIELD%CMNHNAME = 'VR_FLX'
- TZFIELD%CSTDNAME = ''
- TZFIELD%CLONGNAME = 'VR_FLX'
- TZFIELD%CUNITS = 'kg kg-1 m s-1'
- TZFIELD%CDIR = 'XY'
- TZFIELD%CCOMMENT = 'X_Y_Z_VR_FLX'
- TZFIELD%NGRID = 3
- TZFIELD%NTYPE = TYPEREAL
- TZFIELD%NDIMS = 3
- TZFIELD%LTIMEDEP = .TRUE.
- CALL IO_Field_write(TPFILE,TZFIELD,ZFLX(:,:,:))
- END IF
- !
- IF (KSPLT==1 .AND. LLES_CALL) THEN
- CALL SECOND_MNH(ZTIME1)
- !
-!$acc data copy(X_LES_SUBGRID_VRt,X_LES_RES_ddxa_W_SBG_UaRt, &
-!$acc & X_LES_RES_ddxa_Thl_SBG_UaRt,X_LES_RES_ddxa_Rt_SBG_UaRt)
- !
- CALL MYF_DEVICE(ZFLX,ZTMP1_DEVICE)
- CALL LES_MEAN_SUBGRID( ZTMP1_DEVICE, X_LES_SUBGRID_VRt )
- !
- CALL GY_W_VW_DEVICE(PWM,PDYY,PDZZ,PDZY,ZTMP1_DEVICE)
- CALL MZM_DEVICE(ZFLX,ZTMP2_DEVICE)
- !$acc kernels
- ZTMP3_DEVICE(:,:,:) = ZTMP1_DEVICE(:,:,:) * ZTMP2_DEVICE(:,:,:)
- !$acc end kernels
- CALL MYF_DEVICE(ZTMP3_DEVICE, ZTMP4_DEVICE)
- CALL MZF_DEVICE( ZTMP4_DEVICE, ZTMP1_DEVICE )
- CALL LES_MEAN_SUBGRID( ZTMP1_DEVICE,X_LES_RES_ddxa_W_SBG_UaRt , .TRUE. )
- !
- CALL GY_M_M_DEVICE(PTHLM,PDYY,PDZZ,PDZY,ZTMP1_DEVICE)
- CALL MYF_DEVICE(ZFLX,ZTMP2_DEVICE)
- !$acc kernels
- ZTMP3_DEVICE(:,:,:) = ZTMP1_DEVICE(:,:,:) * ZTMP2_DEVICE(:,:,:)
- !$acc end kernels
- CALL LES_MEAN_SUBGRID( ZTMP3_DEVICE, X_LES_RES_ddxa_Thl_SBG_UaRt , .TRUE. )
- !
- CALL GY_M_M_DEVICE(PRM(:,:,:,1),PDYY,PDZZ,PDZY,ZTMP1_DEVICE)
- CALL MYF_DEVICE(ZFLX,ZTMP2_DEVICE)
- !$acc kernels
- ZTMP3_DEVICE(:,:,:) = ZTMP1_DEVICE(:,:,:) * ZTMP2_DEVICE(:,:,:)
- !$acc end kernels
- CALL LES_MEAN_SUBGRID( ZTMP3_DEVICE, X_LES_RES_ddxa_Rt_SBG_UaRt , .TRUE. )
- !
-!$acc end data
- !
- CALL SECOND_MNH(ZTIME2)
- XTIME_LES = XTIME_LES + ZTIME2 - ZTIME1
- END IF
-!
- !
- IF (KRRL>0 .AND. KSPLT==1 .AND. LLES_CALL) THEN
- CALL SECOND_MNH(ZTIME1)
- !
-!$acc data copy(X_LES_SUBGRID_VRc)
- !
- CALL MYF_DEVICE(ZFLXC,ZTMP1_DEVICE)
- CALL LES_MEAN_SUBGRID(ZTMP1_DEVICE, X_LES_SUBGRID_VRc )
- !
-!$acc end data
- !
- CALL SECOND_MNH(ZTIME2)
- XTIME_LES = XTIME_LES + ZTIME2 - ZTIME1
- END IF
- !
-END IF
-#endif
-!
-!* 7. < V' TPV' >
-! -----------
-!
-!! to be tested later
-!!IF (KRR/=0) THEN
-!! ! here ZFLX= <V'R'np> and ZWORK= <V'Theta'l>
-!! !
-!! IF (.NOT. L2D) THEN &
-!! ZVPTV(:,:,:) = &
-!! ZWORK(:,:,:)*MYM(ETHETA(KRR,KRRI,PTHLT,PEXNREF,PRT,PLOCPT,PSRCM)) + &
-!! ZFLX(:,:,:)*MYM(EMOIST(KRR,KRRI,PTHLT,PEXNREF,PRT,PLOCPT,PSRCM))
-!! ELSE
-!! ZVPTV(:,:,:) = 0.
-!! END IF
-!! !
-!! ! stores the horizontal <V VPT>
-!! IF ( tpfile%lopened .AND. OTURB_FLX ) THEN
-!! TZFIELD%CMNHNAME = 'VVPT_FLX'
-!! TZFIELD%CSTDNAME = ''
-!! TZFIELD%CLONGNAME = 'VVPT_FLX'
-!! TZFIELD%CUNITS = 'K m s-1'
-!! TZFIELD%CDIR = 'XY'
-!! TZFIELD%CCOMMENT = 'X_Y_Z_VVPT_FLX'
-!! TZFIELD%NGRID = 3
-!! TZFIELD%NTYPE = TYPEREAL
-!! TZFIELD%NDIMS = 3
-!! TZFIELD%LTIMEDEP = .TRUE.
-!! CALL IO_Field_write(TPFILE,TZFIELD,ZVPTV)
-!! END IF
-!!!
-!!ELSE
-!! ZVPTV(:,:,:)=ZWORK(:,:,:)
-!!END IF
-
-if ( mppdb_initialized ) then
- !Check all inout arrays
- call Mppdb_check( prthls, "Turb_hor_thermo_flux end:prthls" )
- call Mppdb_check( prrs, "Turb_hor_thermo_flux end:prrs" )
-end if
-
-!$acc end data
-
-#ifndef MNH_OPENACC
-deallocate (zflx,zflxc,zcoeff)
-#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 TURB_HOR_THERMO_FLUX