Forked from
Méso-NH / Méso-NH code
2331 commits behind, 928 commits ahead of the upstream repository.
-
ESCOBAR MUNOZ Juan authoredJuan 15/04/2022:MNH/ZSOLVER/turb_hor_thermo_flux.f90, Cray GPU Opt/Bug bypass, add present_cr + acc_nv , where needeed
ESCOBAR MUNOZ Juan authoredJuan 15/04/2022:MNH/ZSOLVER/turb_hor_thermo_flux.f90, Cray GPU Opt/Bug bypass, add present_cr + acc_nv , where needeed
turb_hor_thermo_flux.f90 71.06 KiB
!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