diff --git a/src/ZSOLVER/turb_hor_thermo_flux.f90 b/src/ZSOLVER/turb_hor_thermo_flux.f90
new file mode 100644
index 0000000000000000000000000000000000000000..14cba9ab43501f46ef6dc434816b66e72cd8ba64
--- /dev/null
+++ b/src/ZSOLVER/turb_hor_thermo_flux.f90
@@ -0,0 +1,1818 @@
+!MNH_LIC Copyright 1994-2020 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, &
+ OCLOSE_OUT,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) :: OCLOSE_OUT ! switch for syncronous
+ ! file opening
+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(IN) :: 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, &
+ OCLOSE_OUT,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_IO, ONLY: TFILEDATA
+USE MODD_PARAMETERS
+USE MODD_LES
+!
+USE MODE_FIELD, ONLY: TFIELDDATA, TYPEREAL
+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_ALLOCATE_ZT3D , MNH_REL_ZT3D, MNH_ALLOCATE_ZT3DP , MNH_ALLOCATE_ZT2D, &
+ MNH_ALLOCATE_ZT4D , MNH_REL_ZT4D, &
+ MNH_CHECK_IN_ZT3D,MNH_CHECK_OUT_ZT3D
+#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) :: OCLOSE_OUT ! switch for syncronous
+ ! file opening
+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(IN) :: 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
+INTEGER :: IZFLX,IZFLXC,IZCOEFF
+!
+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
+INTEGER :: IZTMP1_DEVICE, IZTMP2_DEVICE, IZTMP3_DEVICE, IZTMP4_DEVICE, &
+ IZTMP5_DEVICE, IZTMP6_DEVICE, IZTMP7_DEVICE, IZTMP8_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
+CALL MNH_CHECK_IN_ZT3D("TURB_HOR_THERMO_FLUX")
+izflx = MNH_ALLOCATE_ZT3D( zflx ,JIU,JJU,JKU )
+izflxc = MNH_ALLOCATE_ZT3D( zflxc,JIU,JJU,JKU )
+! izvptv= MNH_ALLOCATE_ZT3D( zvptv,JIU,JJU,JKU )
+
+izcoeff= MNH_ALLOCATE_ZT3DP( zcoeff,JIU,JJU, 1 + jpvext , 3 + jpvext )
+
+#endif
+
+#ifdef MNH_OPENACC
+iztmp1_device= MNH_ALLOCATE_ZT3D( ztmp1_device,JIU,JJU,JKU )
+iztmp2_device= MNH_ALLOCATE_ZT3D( ztmp2_device,JIU,JJU,JKU )
+iztmp3_device= MNH_ALLOCATE_ZT3D( ztmp3_device,JIU,JJU,JKU )
+iztmp4_device= MNH_ALLOCATE_ZT3D( ztmp4_device,JIU,JJU,JKU )
+iztmp5_device= MNH_ALLOCATE_ZT3D( ztmp5_device,JIU,JJU,JKU )
+iztmp6_device= MNH_ALLOCATE_ZT3D( ztmp6_device,JIU,JJU,JKU )
+iztmp7_device= MNH_ALLOCATE_ZT3D( ztmp7_device,JIU,JJU,JKU )
+iztmp8_device= MNH_ALLOCATE_ZT3D( 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 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
+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
+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
+!$acc loop independent collapse(2)
+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
+ 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 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 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 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 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(1,IKU,1,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
+ ZTMP1_DEVICE(:,:,:) = ZTMP2_DEVICE(:,:,:) *ZFLX(:,:,:)
+ !$acc end kernels
+ CALL MXF_DEVICE( ZTMP1_DEVICE, ZTMP2_DEVICE)
+ !$acc kernels
+ ZTMP1_DEVICE(:,:,:) = ZFLX(:,:,:)*PINV_PDXX(:,:,:)
+ !$acc end kernels
+ CALL MZM_DEVICE( ZTMP1_DEVICE, ZTMP5_DEVICE )
+ !$acc kernels
+ ZTMP6_DEVICE(:,:,:) = PDZX(:,:,:)*ZTMP5_DEVICE(:,:,:)
+ !$acc end kernels
+ CALL MXF_DEVICE( ZTMP6_DEVICE, ZTMP5_DEVICE )
+ !$acc kernels
+ ZTMP6_DEVICE(:,:,:) = ZTMP4_DEVICE(:,:,:)*ZTMP5_DEVICE(:,:,:)
+ !$acc end kernels
+ CALL MZF_DEVICE(1,IKU,1, ZTMP6_DEVICE,ZTMP7_DEVICE )
+ !$acc kernels
+ 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
+ ZTMP6_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(1,IKU,1, 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(1,IKU,1, 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 ( OCLOSE_OUT .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(1,IKU,1,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(1,IKU,1,ZTMP1_DEVICE,ZTMP2_DEVICE)
+ CALL LES_MEAN_SUBGRID( ZTMP2_DEVICE,X_LES_RES_ddxa_W_SBG_UaThl , .TRUE. )
+ !
+ CALL GX_M_M_DEVICE(1,IKU,1,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(1,IKU,1,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 ( OCLOSE_OUT .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 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
+ 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
+ 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
+ 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 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 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 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 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(1,IKU,1, ZTMP2_DEVICE, ZTMP4_DEVICE)
+ !$acc kernels
+ !$acc 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
+ 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
+ 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
+ ZTMP6_DEVICE(:,:,:) = PDZX(:,:,:)*ZTMP5_DEVICE(:,:,:)
+ !$acc end kernels
+ CALL MXF_DEVICE( ZTMP6_DEVICE, ZTMP5_DEVICE )
+ !$acc kernels
+ ZTMP6_DEVICE(:,:,:) = ZTMP4_DEVICE(:,:,:)*ZTMP5_DEVICE(:,:,:)
+ !$acc end kernels
+ CALL MZF_DEVICE(1,IKU,1, ZTMP6_DEVICE, ZTMP7_DEVICE )
+ !$acc kernels
+ ZFLXC(:,:,:) = ZFLXC(:,:,:) + 2.*( ZTMP3_DEVICE(:,:,:) + ZTMP7_DEVICE(:,:,:) )
+ !$acc end kernels
+ !
+ !
+ !$acc kernels
+ ZTMP6_DEVICE(:,:,:) = ZTMP4_DEVICE(:,:,:)*ZTMP5_DEVICE(:,:,:)*PINV_PDZZ(:,:,:)
+ !$acc end kernels
+ CALL DZF_DEVICE(1,IKU,1, ZTMP6_DEVICE, ZTMP3_DEVICE )
+ !$acc kernels
+ 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 ( OCLOSE_OUT .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(1,IKU,1,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(1,IKU,1,ZTMP4_DEVICE,ZTMP3_DEVICE)
+ CALL LES_MEAN_SUBGRID( ZTMP3_DEVICE, X_LES_RES_ddxa_W_SBG_UaRt , .TRUE. )
+ !
+ CALL GX_M_M_DEVICE(1,IKU,1,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(1,IKU,1,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 ( OCLOSE_OUT .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 ( OCLOSE_OUT .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 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
+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
+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
+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 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 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 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 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(1,IKU,1, ZTMP1_DEVICE, ZTMP2_DEVICE )
+ !$acc kernels
+ !$acc 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
+ ZTMP4_DEVICE(:,:,:) = ZTMP2_DEVICE(:,:,:)*ZFLX(:,:,:)
+ !$acc end kernels
+ CALL MYF_DEVICE( ZTMP4_DEVICE, ZTMP3_DEVICE )
+ !$acc kernels
+ ZTMP4_DEVICE(:,:,:) = ZFLX(:,:,:)*PINV_PDYY(:,:,:)
+ !$acc end kernels
+ CALL MZM_DEVICE( ZTMP4_DEVICE, ZTMP5_DEVICE )
+ !$acc kernels
+ ZTMP4_DEVICE(:,:,:) = PDZY(:,:,:)*ZTMP5_DEVICE(:,:,:)
+ !$acc end kernels
+ CALL MYF_DEVICE( ZTMP4_DEVICE, ZTMP5_DEVICE)
+ CALL MZM_DEVICE( ZTMP1_DEVICE, ZTMP4_DEVICE )
+ !$acc kernels
+ ZTMP6_DEVICE(:,:,:) = ZTMP4_DEVICE(:,:,:)*ZTMP5_DEVICE(:,:,:)
+ !$acc end kernels
+ CALL MZF_DEVICE(1,IKU,1, ZTMP6_DEVICE, ZTMP4_DEVICE )
+ !$acc kernels
+ ZFLXC(:,:,:) = 2.*( ZTMP3_DEVICE(:,:,:) + ZTMP4_DEVICE(:,:,:) )
+ !$acc end kernels
+ IF ( KRRI >= 1 ) THEN
+ !$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(1,IKU,1, 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(1,IKU,1, 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 ( OCLOSE_OUT .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(1,IKU,1,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(1,IKU,1,ZTMP4_DEVICE, ZTMP1_DEVICE)
+ CALL LES_MEAN_SUBGRID( ZTMP1_DEVICE, X_LES_RES_ddxa_W_SBG_UaThl , .TRUE. )
+ !
+ CALL GY_M_M_DEVICE(1,IKU,1,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(1,IKU,1,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(:,:,:)/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(:,:,:)/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(:,:,:)/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(:,:,:)/PDYY )*(1.0-PFRAC_ICE(:,:,:))
+ PRRS(:,:,:,4) = PRRS(:,:,:,4) - 2. * &
+ DYF( MYM( PRHODJ(:,:,:)*PAMOIST(:,:,:)*PSRCM(:,:,:) )*ZFLX(:,:,:)/PDYY )*PFRAC_ICE(:,:,:)
+ ELSE
+ PRRS(:,:,:,2) = PRRS(:,:,:,2) - 2. * &
+ DYF( MYM( PRHODJ(:,:,:)*PAMOIST(:,:,:)*PSRCM(:,:,:) )*ZFLX(:,:,:)/PDYY )
+ END IF
+ END IF
+ END IF
+ !
+ ! stores the horizontal <V Rnp>
+ IF ( OCLOSE_OUT .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 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
+ 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
+ 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
+ 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 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 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 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 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(1,IKU,1,ZTMP1_DEVICE,ZTMP2_DEVICE )
+ !
+ !$acc kernels
+ !$acc 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
+ ZTMP3_DEVICE(:,:,:) = ZTMP2_DEVICE(:,:,:)*ZFLX(:,:,:)
+ !$acc end kernels
+ CALL MXF_DEVICE( ZTMP3_DEVICE, ZTMP4_DEVICE )
+ CALL MZM_DEVICE( ZTMP1_DEVICE, ZTMP5_DEVICE )
+ !$acc kernels
+ ZTMP1_DEVICE(:,:,:) = ZFLX(:,:,:)*PINV_PDYY(:,:,:)
+ !$acc end kernels
+ CALL MZM_DEVICE( ZTMP1_DEVICE, ZTMP2_DEVICE )
+ !$acc kernels
+ ZTMP1_DEVICE(:,:,:) = PDZY(:,:,:)*ZTMP2_DEVICE(:,:,:)
+ !$acc end kernels
+ CALL MYF_DEVICE( ZTMP1_DEVICE, ZTMP2_DEVICE )
+ !$acc kernels
+ ZTMP1_DEVICE(:,:,:) = ZTMP5_DEVICE(:,:,:)*ZTMP2_DEVICE(:,:,:)
+ !$acc end kernels
+ CALL MZF_DEVICE(1,IKU,1, ZTMP1_DEVICE, ZTMP2_DEVICE )
+ !$acc kernels
+ ZFLXC(:,:,:) = ZFLXC(:,:,:) + 2.*( ZTMP4_DEVICE(:,:,:) + ZTMP2_DEVICE(:,:,:) )
+ !$acc end kernels
+ IF ( KRRI >= 1 ) THEN
+ !$acc kernels
+ ZTMP2_DEVICE(:,:,:) = ZTMP3_DEVICE(:,:,:) * PINV_PDYY(:,:,:)
+ !$acc end kernels
+ CALL DYF_DEVICE( ZTMP2_DEVICE, ZTMP3_DEVICE )
+ !$acc kernels
+ ZTMP2_DEVICE(:,:,:) = ZTMP1_DEVICE(:,:,:)* PINV_PDZZ(:,:,:)
+ !$acc end kernels
+ CALL DZF_DEVICE(1,IKU,1, 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 ( OCLOSE_OUT .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(1,IKU,1,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(1,IKU,1,ZTMP4_DEVICE,ZTMP1_DEVICE)
+ CALL LES_MEAN_SUBGRID( ZTMP1_DEVICE,X_LES_RES_ddxa_W_SBG_UaRt , .TRUE. )
+ !
+ CALL GY_M_M_DEVICE(1,IKU,1,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(1,IKU,1,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 ( OCLOSE_OUT .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
+CALL MNH_REL_ZT3D ( IZFLX, IZFLXC, IZCOEFF, &
+ IZTMP1_DEVICE, IZTMP2_DEVICE, IZTMP3_DEVICE, IZTMP4_DEVICE, &
+ IZTMP5_DEVICE, IZTMP6_DEVICE, IZTMP7_DEVICE, IZTMP8_DEVICE )
+CALL MNH_CHECK_OUT_ZT3D("TURB_HOR_THERMO_FLUX")
+#endif
+
+!$acc end data
+
+END SUBROUTINE TURB_HOR_THERMO_FLUX