diff --git a/src/ZSOLVER/dyn_sources.f90 b/src/ZSOLVER/dyn_sources.f90
new file mode 100644
index 0000000000000000000000000000000000000000..aa8feb2fc99f4836c19fd04492f442826b0457a1
--- /dev/null
+++ b/src/ZSOLVER/dyn_sources.f90
@@ -0,0 +1,607 @@
+!MNH_LIC Copyright 1994-2021 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
+!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
+!MNH_LIC for details. version 1.
+!-----------------------------------------------------------------
+! #######################
+ MODULE MODI_DYN_SOURCES
+! #######################
+!
+INTERFACE
+!
+ SUBROUTINE DYN_SOURCES( KRR,KRRL, KRRI, &
+ PUT, PVT, PWT, PTHT, PRT, &
+ PCORIOX, PCORIOY, PCORIOZ, PCURVX, PCURVY, &
+ PRHODJ, PZZ, PTHVREF, PEXNREF, &
+ PRUS, PRVS, PRWS, PRTHS )
+!
+INTEGER, INTENT(IN) :: KRR ! Total number of water var.
+INTEGER, INTENT(IN) :: KRRL ! Number of liquid water var.
+INTEGER, INTENT(IN) :: KRRI ! Number of ice water var.
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PUT, PVT, PWT ! variables
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHT ! at
+REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PRT ! time t
+!
+REAL, DIMENSION(:,:), INTENT(IN) :: PCORIOX ! 2D horizontal
+REAL, DIMENSION(:,:), INTENT(IN) :: PCORIOY ! Coriolis
+REAL, DIMENSION(:,:), INTENT(IN) :: PCORIOZ ! factors
+REAL, DIMENSION(:,:), INTENT(IN) :: PCURVX ! 2D horizontal
+REAL, DIMENSION(:,:), INTENT(IN) :: PCURVY ! curvature factors
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODJ ! Density * Jacobian
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PZZ ! Height (z)
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHVREF ! Virtual Temperature
+ ! of the reference state
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PEXNREF ! Exner function
+ ! of the reference state
+!
+REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PRUS, PRVS, PRWS ! Sources of Momentum
+REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PRTHS ! Sources of theta
+!
+END SUBROUTINE DYN_SOURCES
+!
+END INTERFACE
+!
+END MODULE MODI_DYN_SOURCES
+! ######################################################################
+ SUBROUTINE DYN_SOURCES( KRR,KRRL, KRRI, &
+ PUT, PVT, PWT, PTHT, PRT, &
+ PCORIOX, PCORIOY, PCORIOZ, PCURVX, PCURVY, &
+ PRHODJ, PZZ, PTHVREF, PEXNREF, &
+ PRUS, PRVS, PRWS, PRTHS )
+! ######################################################################
+!
+!!**** *DYN_SOURCES * - routine to compute the curvature, coriolis and gravity terms
+!!
+!! PURPOSE
+!! -------
+!! The purpose of this routine is to compute the dynamical sources
+!! (curvature, coriolis and gravity terms) for each component of the
+!! momentum. The curvature terms arise in case of non-cartesian geometry
+!! only. For both curvature and coriolis terms, the "thin shell"
+!! approximation can be assumed or not. Gravity affects only the vertical
+!! component of the momentum.
+!! This routine also adds the Lagrangien derivative of the pressure as a
+!! source for the theta evolution.
+!!
+!!
+!!** METHOD
+!! ------
+!! The horizontal curvature and coriolis field arrays are expanded in
+!! the vertical direction with the SPREAD intrinsics to match with the
+!! Shuman operators. The source term for theta due to the advection of the
+!! Exner function is treated in an advective form (or non-conservative form).
+!! Values of the calorific capacity of the melting and of the potential
+!! virtual temperature are recovered taking care of the number of water
+!! variables. The different sources terms are stored for the budget
+!! computations.
+!!
+!! EXTERNAL
+!! --------
+!! MXM,MYM,MZM : Shuman functions (mean operators)
+!! MXF,MYF,MZF : Shuman functions (mean operators)
+!! GZ_M_W : projection along the vertical direction of the gradient
+!! vector. It acts on a field localized in mass point and
+!! the result is localized in w point.
+!! BUDGET : Stores the different budget components
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!! Module MODD_CONF : contains configuration variables for all models
+!! LTHINSHELL = .TRUE. if the THINSHELL approximation is made
+!! LCARTESIAN = .TRUE. if the CARTESIAN approximation is made
+!! Module MODD_CST :
+!! XG gravity acceleration
+!! XRADIUS Earth radius
+!! XRD,XRV Gas constant for dry air and wator vapor
+!! XCPD,XCPV Cp for dry air and wator vapor
+!! XCL,XCI C (calorific capacity) for liquid and solid water
+!! Module MODD_DYN : contains the parameters for the dynamics
+!! LCORIO = .FALSE. if the earth rotation is neglected
+!!
+!! Module MODD_BUDGET:
+!! NBUMOD : model in which budget is calculated
+!! CBUTYPE : type of desired budget
+!! 'CART' for cartesian box configuration
+!! 'MASK' for budget zone defined by a mask
+!! 'NONE' ' for no budget
+!! LBU_RTH : logical for budget of RTH (potential temperature)
+!! .TRUE. = budget of RTH
+!! .FALSE. = no budget of RTH
+!! LBU_RU : logical for budget of RU (wind component along x)
+!! .TRUE. = budget of RU
+!! .FALSE. = no budget of RU
+!! LBU_RV : logical for budget of RV (wind component along y)
+!! .TRUE. = budget of RV
+!! .FALSE. = no budget of RV
+!! LBU_RW : logical for budget of RW (wind component along z)
+!! .TRUE. = budget of RW
+!! .FALSE. = no budget of RW
+!!
+!! REFERENCE
+!! ---------
+!! Book2 of documentation ( routine DYN_SOURCE )
+!!
+!! AUTHOR
+!! ------
+!! J.-P. Pinty * Laboratoire d'Aerologie*
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 06/06/94
+!! Corrections 06/08/94 (J.-P. Lafore)
+!! Corrections 17/10/94 (Stein) For LCORIO
+!! Corrections 22/12/94 (Stein) add the pressure term in the theta evolution
+!! Corrections 30/12/94 (J.P. Lafore) bug corrections for the pressure term
+!! Corrections 16/03/95 (Stein) remove R from the historical variables and
+!! correction of the pressure term
+!! Corrections 14/04/95 (Masson Stein) bugs in the curvatureand Coriolis
+!! terms for PRUS,PRVS, PRWS
+!! Corrections 01/04/95 (Ph. Hereil J. Nicolau) add the budget computation
+!! Corrections 16/10/95 (J. Stein) change the budget calls
+!! Corrections 16/02/96 (J.-P. Pinty) Introduce L1D switches
+!! Corrections 19/12/96 (J.-P. Pinty) Update the CALL BUDGET
+!! Corrections 03/12/02 (P. Jabouille) add no thinshell condition
+!! Correction 06/10 (C.Lac) Exclude L1D for Coriolis term
+!! Modification 03/11 (C.Lac) Split the gravity term due to buoyancy
+!! J.Escobar : 15/09/2015 : WENO5 & JPHEXT <> 1
+! P. Wautelet 02/2020: use the new data structures and subroutines for budgets
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+use modd_budget, only: lbudget_u, lbudget_v, lbudget_w, lbudget_th, &
+ NBUDGET_U, NBUDGET_V, NBUDGET_W, NBUDGET_TH, &
+ tbudgets
+USE MODD_CONF
+USE MODD_CST
+USE MODD_DYN
+USE MODD_DYN_n, ONLY: LOCEAN
+!
+use mode_budget, only: Budget_store_init, Budget_store_end
+USE MODE_MPPDB
+
+USE MODI_GRADIENT_M
+USE MODI_SHUMAN
+#ifdef MNH_OPENACC
+USE MODI_SHUMAN_DEVICE
+#endif
+#ifdef MNH_OPENACC
+USE MODE_DEVICE
+USE MODE_MNH_ZWORK, ONLY: MNH_MEM_GET, MNH_MEM_POSITION_PIN, MNH_MEM_RELEASE
+#endif
+
+
+IMPLICIT NONE
+!
+!* 0.1 Declarations of dummy arguments :
+!
+INTEGER, INTENT(IN) :: KRR ! Total number of water var.
+INTEGER, INTENT(IN) :: KRRL ! Number of liquid water var.
+INTEGER, INTENT(IN) :: KRRI ! Number of ice water var.
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PUT, PVT, PWT ! variables
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHT ! at
+REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PRT ! time t
+!
+REAL, DIMENSION(:,:), INTENT(IN) :: PCORIOX ! 2D horizontal
+REAL, DIMENSION(:,:), INTENT(IN) :: PCORIOY ! Coriolis
+REAL, DIMENSION(:,:), INTENT(IN) :: PCORIOZ ! factors
+REAL, DIMENSION(:,:), INTENT(IN) :: PCURVX ! 2D horizontal
+REAL, DIMENSION(:,:), INTENT(IN) :: PCURVY ! curvature factors
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODJ ! dry Density * Jacobian
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PZZ ! Height (z)
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHVREF ! Virtual Temperature
+ ! of the reference state
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PEXNREF ! Exner function
+ ! of the reference state
+!
+REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PRUS, PRVS, PRWS ! Sources of Momentum
+REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PRTHS ! Sources of theta
+!
+!
+!* 0.2 Declarations of local variables :
+!
+REAL :: ZCPD_OV_RD ! = CPD / RD
+REAL :: ZG_OV_CPD ! =-XG / XCPD
+INTEGER :: JWATER ! loop index on the different types of water
+INTEGER :: JIU, JJU, JKU
+INTEGER :: JI,JJ,JK
+INTEGER :: IKU ! array size along the k direction
+REAL :: ZD1 ! DELTA1 (switch 0/1) for thinshell approximation
+#ifndef MNH_OPENACC
+REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3)) :: &
+ ZWORK1, ZWORK2, ZWORK3, ZRUT, ZRVT
+#else
+REAL, DIMENSION(:,:,:), pointer , contiguous ::ZWORK1, ZWORK2, ZWORK3, ZRUT, ZRVT
+#endif
+#ifdef MNH_OPENACC
+REAL, DIMENSION(:,:,:), pointer , contiguous :: ZTMP1_DEVICE,ZTMP2_DEVICE,ZTMP3_DEVICE,ZTMP4_DEVICE,ZRVTT,ZRUTT
+REAL, DIMENSION(:,:,:), pointer , contiguous :: ZTMP5_DEVICE,ZTMP6_DEVICE,ZTMP7_DEVICE,ZTMP8_DEVICE
+!!!REAL, DIMENSION(:,:,:), pointer , contiguous :: ZTMP9_DEVICE,ZTMP10_DEVICE,ZTMP11_DEVICE,ZTMP12_DEVICE
+#endif
+!
+!-------------------------------------------------------------------------------
+!
+
+JIU = size( PUT, 1 )
+JJU = size( PUT, 2 )
+JKU = size( PUT, 3 )
+
+#ifdef MNH_OPENACC
+CALL MNH_MEM_POSITION_PIN()
+
+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 )
+!CALL MNH_MEM_GET( ztmp9_device, JIU, JJU, JKU )
+!CALL MNH_MEM_GET( ztmp10_device, JIU, JJU, JKU )
+!CALL MNH_MEM_GET( ztmp11_device, JIU, JJU, JKU )
+!CALL MNH_MEM_GET( ztmp12_device, JIU, JJU, JKU )
+CALL MNH_MEM_GET( ZRVTT, JIU, JJU, JKU )
+CALL MNH_MEM_GET( ZRUTT, JIU, JJU, JKU )
+CALL MNH_MEM_GET( ZWORK1, JIU, JJU, JKU )
+CALL MNH_MEM_GET( ZWORK2, JIU, JJU, JKU )
+CALL MNH_MEM_GET( ZWORK3, JIU, JJU, JKU )
+CALL MNH_MEM_GET( ZRUT, JIU, JJU, JKU )
+CALL MNH_MEM_GET( ZRVT, JIU, JJU, JKU )
+
+#endif
+
+!
+!* 1. COMPUTES THE TRUE VELOCITY COMPONENTS
+! -------------------------------------
+!
+
+!$acc data copyin ( PUT, PVT, PWT,PTHT, PRT, &
+!$acc & PCORIOX, PCORIOY, PCORIOZ, PCURVX, PCURVY, &
+!$acc & PZZ, PTHVREF, PEXNREF) &
+!$acc & present (PRUS,PRVS, PRWS, PRTHS,PRHODJ)
+!$acc update device(PRUS, PRVS, PRWS,PRTHS,PRHODJ)
+!$acc update device ( PUT, PVT, PWT,PTHT, PRT, &
+!$acc & PCORIOX, PCORIOY, PCORIOZ, PCURVX, PCURVY, &
+!$acc & PZZ, PTHVREF, PEXNREF)
+
+
+#ifndef MNH_OPENACC
+ZRUT(:,:,:) = PUT(:,:,:) * MXM(PRHODJ(:,:,:))
+ZRVT(:,:,:) = PVT(:,:,:) * MYM(PRHODJ(:,:,:))
+
+print *,"PUT(1,1,1)",PUT(1,1,1)
+#else
+CALL MXM_DEVICE(PRHODJ(:,:,:), ZTMP1_DEVICE (:,:,:))
+CALL MYM_DEVICE(PRHODJ(:,:,:), ZTMP2_DEVICE (:,:,:))
+!$acc kernels present(ZRUT,ZRVT,PUT,PVT,ZTMP1_DEVICE,ZTMP2_DEVICE)
+ZRUT(:,:,:) = PUT(:,:,:) * ZTMP1_DEVICE(:,:,:)
+ZRVT(:,:,:) = PVT(:,:,:) * ZTMP2_DEVICE(:,:,:)
+!$acc end kernels
+#endif
+!
+IKU = SIZE(PUT,3)
+!
+!-------------------------------------------------------------------------------
+!
+!* 2. COMPUTES THE CURVATURE TERMS
+! ----------------------------
+!
+! Only when earth rotation is considered but not in 1D and CARTESIAN cases
+!
+IF ((.NOT.L1D).AND.(.NOT.LCARTESIAN) ) THEN
+ if ( lbudget_u ) call Budget_store_init( tbudgets(NBUDGET_U), 'CURV', prus(:, :, :) )
+ if ( lbudget_v ) call Budget_store_init( tbudgets(NBUDGET_V), 'CURV', prvs(:, :, :) )
+
+ IF ( LTHINSHELL ) THEN ! THINSHELL approximation
+#ifndef MNH_OPENACC
+ ZWORK1(:,:,:) = SPREAD( PCURVX(:,:),DIM=3,NCOPIES=JKU ) / XRADIUS
+ ZWORK2(:,:,:) = SPREAD( PCURVY(:,:),DIM=3,NCOPIES=JKU ) / XRADIUS
+#else
+ !$acc kernels present_cr(ZWORK1,ZWORK2)
+ !$mnh_expand_array(JI=1:JIU , JJ=1:JJU , JK=1:JKU)
+ ZWORK1(:,:,:) = PCURVX(:,:) / XRADIUS
+ ZWORK2(:,:,:) = PCURVY(:,:) / XRADIUS
+ !$mnh_end_expand_array()
+!$acc end kernels
+#endif
+#ifndef MNH_OPENACC
+ PRUS(:,:,:) = PRUS &
+ + ZRUT * MXM( MYF(PVT) * ZWORK1 ) &
+ + MXM( MYF(ZRVT*PVT) * ZWORK2 )
+!
+ PRVS(:,:,:) = PRVS &
+ - MYM( MXF(ZRUT*PUT) * ZWORK1 ) &
+ - ZRVT * MYM( MXF(PUT) * ZWORK2 )
+!
+#else
+!$acc kernels
+ZRVTT(:,:,:)=ZRVT(:,:,:)*PVT(:,:,:)
+ZRUTT(:,:,:)=ZRUT(:,:,:)*PUT(:,:,:)
+!$acc end kernels
+CALL MYF_DEVICE(PVT,ZTMP1_DEVICE ) !MYF(PVT)
+CALL MYF_DEVICE(ZRVTT,ZTMP2_DEVICE ) !MYF(ZRVT*PVT)
+CALL MXF_DEVICE(ZRUTT,ZTMP3_DEVICE ) !MXF(ZRUT*PUT)
+CALL MXF_DEVICE(PUT,ZTMP4_DEVICE ) !MXF(PUT)
+
+!$acc kernels present_cr(ZTMP1_DEVICE,ZTMP1_DEVICE,ZTMP1_DEVICE,ZTMP1_DEVICE)
+ZTMP1_DEVICE(:,:,:)=ZTMP1_DEVICE(:,:,:)* ZWORK1(:,:,:) !MYF(PVT) * ZWORK1
+ZTMP2_DEVICE(:,:,:)=ZTMP2_DEVICE(:,:,:)* ZWORK2(:,:,:) !MYF(ZRVT*PVT) * ZWORK2
+ZTMP3_DEVICE(:,:,:)=ZTMP3_DEVICE(:,:,:)* ZWORK1(:,:,:) !MXF(ZRUT*PUT) * ZWORK1
+ZTMP4_DEVICE(:,:,:)=ZTMP4_DEVICE(:,:,:) * ZWORK2(:,:,:) !MXF(PUT) * ZWORK2
+!$acc end kernels
+
+CALL MXM_DEVICE(ZTMP1_DEVICE,ZTMP5_DEVICE) !MXM( MYF(PVT) * ZWORK1 )
+CALL MXM_DEVICE(ZTMP2_DEVICE,ZTMP6_DEVICE) !MXM( MYF(ZRVT*PVT) * ZWORK2 )
+
+CALL MYM_DEVICE(ZTMP3_DEVICE,ZTMP7_DEVICE) !MYM( MXF(ZRUT*PUT) * ZWORK1 )
+CALL MYM_DEVICE(ZTMP4_DEVICE,ZTMP8_DEVICE) !MYM( MXF(PUT) * ZWORK2 )
+
+!$acc kernels present_cr(PRUS,PRVS)
+PRUS(:,:,:) = PRUS(:,:,:) &
++ ZRUT(:,:,:) * ZTMP5_DEVICE(:,:,:) &
+ + ZTMP6_DEVICE(:,:,:)
+
+PRVS(:,:,:) = PRVS(:,:,:) &
+ - ZTMP7_DEVICE(:,:,:) &
+ - ZRVT(:,:,:) * ZTMP8_DEVICE(:,:,:)
+
+!$acc end kernels
+#endif
+
+ ELSE ! NO THINSHELL approximation
+ if ( lbudget_w ) call Budget_store_init( tbudgets(NBUDGET_W), 'CURV', prws(:, :, :) )
+#ifndef MNH_OPENACC
+ ZWORK3(:,:,:) = 1.0 / ( XRADIUS + MZF(PZZ(:,:,:)) )
+ ZWORK1(:,:,:) = SPREAD( PCURVX(:,:),DIM=3,NCOPIES=JKU )
+ ZWORK2(:,:,:) = SPREAD( PCURVY(:,:),DIM=3,NCOPIES=JKU )
+#else
+CALL MZF_DEVICE(PZZ,ZWORK3) !MZF(PZZ(:,:,:))
+!$acc kernels present_cr(ZWORK3,ZWORK1,ZWORK2)
+ ZWORK3(:,:,:) = 1.0 / ( XRADIUS + ZWORK3(:,:,:) )
+ ZWORK1(:,:,:) = SPREAD( PCURVX(:,:),DIM=3,NCOPIES=JKU )
+ ZWORK2(:,:,:) = SPREAD( PCURVY(:,:),DIM=3,NCOPIES=JKU )
+!$acc end kernels
+#endif
+ IF (MPPDB_INITIALIZED) THEN
+ CALL MPPDB_CHECK3DM("DYN_SOURCES:ZWORK3,ZWORK1,ZWORK2",PRECISION,&
+ & ZWORK3,ZWORK1,ZWORK2,&
+ & MXM( MYF(ZRVT*PVT) * ZWORK2 * ZWORK3 ) , &
+ & MXM( ( MYF(PVT) * ZWORK1 - MZF(PWT) ) * ZWORK3 ) ,&
+ & MYF(PVT) * ZWORK1 - MZF(PWT) , &
+ & MYF(PVT) , MZF(PWT) , MXM(PWT) , MYM(PWT) )
+ CALL MPPDB_CHECK3DM("DYN_SOOURCES:SUITE",PRECISION,&
+ & MXM(ZRVT),MXM(PVT),MXM(PWT),MXM(ZWORK1),MXM(ZWORK2),MXM(ZWORK3),&
+ & ZRUT,ZRVT,PRUS,PRVS,PRWS )
+ ENDIF
+!
+#ifndef MNH_OPENACC
+ PRUS(:,:,:) = PRUS &
+ + MXM( MYF(ZRVT*PVT) * ZWORK2 * ZWORK3 ) &
+ + ZRUT * MXM( ( MYF(PVT) * ZWORK1 - MZF(PWT) ) * ZWORK3 )
+!
+ PRVS(:,:,:) = PRVS &
+ - MYM( MXF(ZRUT*PUT) * ZWORK1 * ZWORK3 ) &
+ - ZRVT * MYM( (MXF(PUT) * ZWORK2 + MZF(PWT) ) * ZWORK3 )
+!
+ PRWS(:,:,:) = PRWS &
+ +MZM( ( MXF(ZRUT*PUT) + MYF(ZRVT*PVT) ) * ZWORK3 )
+#else
+!$acc kernels present_cr(ZRVTT,ZRUTT)
+ZRVTT(:,:,:)=ZRVT(:,:,:)*PVT(:,:,:)
+ZRUTT(:,:,:)=ZRUT(:,:,:)*PUT(:,:,:)
+
+!PRUS(:,:,:) = PRUS &
+! + MXM( MYF(ZRVT*PVT) * ZWORK2 * ZWORK3 ) &
+! + ZRUT * MXM( ( MYF(PVT) * ZWORK1 - MZF(PWT) ) * ZWORK3 )
+!$acc end kernels
+CALL MYF_DEVICE(ZRVTT,ZTMP1_DEVICE ) !MYF(ZRVT*PVT)
+CALL MYF_DEVICE(PVT,ZTMP2_DEVICE ) !MYF(PVT)
+CALL MZF_DEVICE(PWT,ZTMP3_DEVICE ) !MZF(PWT)
+!$acc kernels present_cr(ZTMP1_DEVICE,ZTMP4_DEVICE)
+ZTMP1_DEVICE(:,:,:)=ZTMP1_DEVICE(:,:,:)* ZWORK2(:,:,:) * ZWORK3(:,:,:) ! MYF(ZRVT*PVT) * ZWORK2 * ZWORK3
+ZTMP4_DEVICE(:,:,:)=(ZTMP2_DEVICE(:,:,:)* ZWORK1(:,:,:) - ZTMP3_DEVICE(:,:,:)) * ZWORK3(:,:,:) !( MYF(PVT) * ZWORK1 - MZF(PWT) ) * ZWORK3
+!$acc end kernels
+CALL MXM_DEVICE(ZTMP1_DEVICE,ZTMP2_DEVICE) ! MXM( MYF(ZRVT*PVT) * ZWORK2 * ZWORK3 )
+CALL MXM_DEVICE(ZTMP4_DEVICE,ZTMP5_DEVICE) ! MXM( ( MYF(PVT) * ZWORK1 - MZF(PWT) ) * ZWORK3 )
+!$acc kernels present_cr(PRUS)
+ PRUS(:,:,:) = PRUS(:,:,:) &
+ +ZTMP2_DEVICE(:,:,:) &
+ + ZRUT(:,:,:) * ZTMP5_DEVICE(:,:,:)
+!$acc end kernels
+
+! PRVS(:,:,:) = PRVS &
+! - MYM( MXF(ZRUT*PUT) * ZWORK1 * ZWORK3 ) &
+! - ZRVT * MYM( (MXF(PUT) * ZWORK2 + MZF(PWT) ) * ZWORK3 )
+CALL MXF_DEVICE(ZRUTT,ZTMP1_DEVICE ) !MXF(ZRUT*PUT)
+CALL MXF_DEVICE(PUT,ZTMP2_DEVICE ) !MXF(PUT)
+!CALL MZF_DEVICE(PWT,ZTMP3_DEVICE ) !MZF(PWT)
+!$acc kernels present_cr(ZTMP1_DEVICE,ZTMP4_DEVICE)
+ZTMP1_DEVICE(:,:,:)=ZTMP1_DEVICE(:,:,:)* ZWORK1(:,:,:) * ZWORK3(:,:,:) !MXF(ZRUT*PUT) * ZWORK1 * ZWORK3
+ZTMP4_DEVICE(:,:,:)=(ZTMP2_DEVICE(:,:,:)* ZWORK2(:,:,:) + ZTMP3_DEVICE(:,:,:)) * ZWORK3(:,:,:) !(MXF(PUT) * ZWORK2 + MZF(PWT) ) * ZWORK3
+!$acc end kernels
+CALL MYM_DEVICE(ZTMP1_DEVICE,ZTMP2_DEVICE) !MYM( MXF(ZRUT*PUT) * ZWORK1 * ZWORK3 )
+CALL MYM_DEVICE(ZTMP4_DEVICE,ZTMP5_DEVICE) !MYM( (MXF(PUT) * ZWORK2 + MZF(PWT) ) * ZWORK3 )
+!$acc kernels present_cr(PRVS)
+ PRVS(:,:,:) = PRVS(:,:,:) &
+ - ZTMP2_DEVICE(:,:,:) &
+ - ZRVT(:,:,:) * ZTMP5_DEVICE(:,:,:)
+!$acc end kernels
+
+! PRWS(:,:,:) = PRWS &
+! +MZM( ( MXF(ZRUT*PUT) + MYF(ZRVT*PVT) ) * ZWORK3 )
+CALL MYF_DEVICE(ZRVTT,ZTMP1_DEVICE ) !MXF(ZRUT*PUT)
+CALL MXF_DEVICE(ZRUTT,ZTMP2_DEVICE ) !MYF(ZRVT*PVT)
+!$acc kernels present_cr(ZTMP3_DEVICE)
+ZTMP3_DEVICE(:,:,:)=(ZTMP2_DEVICE(:,:,:)+ZTMP1_DEVICE(:,:,:)) * ZWORK3(:,:,:) ! (MXF(ZRUT*PUT) + MYF(ZRVT*PVT) ) * ZWORK3
+!$acc end kernels
+CALL MZM_DEVICE(ZTMP3_DEVICE,ZTMP1_DEVICE) ! MZM( ( MXF(ZRUT*PUT) + MYF(ZRVT*PVT) ) * ZWORK3 )
+!$acc kernels present_cr(PRWS)
+PRWS(:,:,:) = PRWS(:,:,:)+ZTMP1_DEVICE(:,:,:)
+!$acc end kernels
+#endif
+
+
+ if ( lbudget_w ) call Budget_store_end( tbudgets(NBUDGET_W), 'CURV', prws(:, :, :) )
+ END IF
+
+ if ( lbudget_u ) call Budget_store_end( tbudgets(NBUDGET_U), 'CURV', prus(:, :, :) )
+ if ( lbudget_v ) call Budget_store_end( tbudgets(NBUDGET_V), 'CURV', prvs(:, :, :) )
+END IF
+!
+!
+!-------------------------------------------------------------------------------
+!
+!* 3. COMPUTES THE CORIOLIS TERMS
+! ---------------------------
+!
+IF (LCORIO) THEN
+ if ( lbudget_u ) call Budget_store_init( tbudgets(NBUDGET_U), 'COR', prus(:, :, :) )
+ if ( lbudget_v ) call Budget_store_init( tbudgets(NBUDGET_V), 'COR', prvs(:, :, :) )
+!$acc kernels present_cr(ZWORK3)
+ ZWORK3(:,:,:) = SPREAD( PCORIOZ(:,:),DIM=3,NCOPIES=JKU ) * PRHODJ(:,:,:)
+!$acc end kernels
+#ifndef MNH_OPENACC
+PRUS(:,:,:) = PRUS + MXM( ZWORK3 * MYF(PVT) )
+PRVS(:,:,:) = PRVS - MYM( ZWORK3 * MXF(PUT) )
+#else
+CALL MYF_DEVICE(PVT,ZTMP1_DEVICE) !MYF(PVT)
+CALL MXF_DEVICE(PUT,ZTMP2_DEVICE) !MXF(PUT)
+!$acc kernels present_cr(ZTMP1_DEVICE,ZTMP2_DEVICE)
+ZTMP1_DEVICE(:,:,:)=ZWORK3(:,:,:) * ZTMP1_DEVICE(:,:,:)
+ZTMP2_DEVICE(:,:,:)=ZWORK3(:,:,:) * ZTMP2_DEVICE(:,:,:)
+!$acc end kernels
+CALL MXM_DEVICE(ZTMP1_DEVICE, ZTMP3_DEVICE) ! MXM( ZWORK3 * MYF(PVT)
+CALL MYM_DEVICE(ZTMP2_DEVICE, ZTMP4_DEVICE) ! MYM( ZWORK3 * MXF(PUT)
+!$acc kernels present_cr(PRUS,PRVS)
+PRUS(:,:,:) = PRUS(:,:,:) + ZTMP3_DEVICE(:,:,:)
+PRVS(:,:,:) = PRVS(:,:,:) - ZTMP4_DEVICE(:,:,:)
+!$acc end kernels
+#endif
+!
+
+IF ((.NOT.LTHINSHELL) .AND. (.NOT.L1D)) THEN
+if ( lbudget_w ) call Budget_store_init( tbudgets(NBUDGET_W), 'COR', prws(:, :, :) )
+!$acc kernels present_cr(ZWORK1,ZWORK2)
+ZWORK1(:,:,:) = SPREAD( PCORIOX(:,:),DIM=3,NCOPIES=JKU) * PRHODJ(:,:,:)
+ZWORK2(:,:,:) = SPREAD( PCORIOY(:,:),DIM=3,NCOPIES=JKU) * PRHODJ(:,:,:)
+!$acc end kernels
+
+CALL MPPDB_CHECK3DM("DYN_SOOURCES:CORIOLIS",PRECISION,&
+ & ZWORK1,ZWORK2,ZWORK3 )
+!
+#ifndef MNH_OPENACC
+ PRUS(:,:,:) = PRUS - MXM( ZWORK2 * MZF(PWT) )
+!
+ PRVS(:,:,:) = PRVS - MYM( ZWORK1 * MZF(PWT) )
+!
+ PRWS(:,:,:) = PRWS + MZM( ZWORK2 * MXF(PUT) + ZWORK1 * MYF(PVT) )
+#else
+!***PRUS(:,:,:) = PRUS - MXM( ZWORK2 * MZF(PWT) )
+CALL MZF_DEVICE(PWT,ZTMP1_DEVICE) !MZF(PWT)
+!$acc kernels present_cr(ZTMP2_DEVICE)
+ZTMP2_DEVICE(:,:,:)=ZWORK2(:,:,:) * ZTMP1_DEVICE(:,:,:) !ZWORK2 * MZF(PWT)
+!$acc end kernels
+CALL MXM_DEVICE(ZTMP2_DEVICE,ZTMP3_DEVICE) !MXM( ZWORK2 * MZF(PWT) )
+!$acc kernels present_cr(PRUS)
+ PRUS(:,:,:) = PRUS(:,:,:) - ZTMP3_DEVICE(:,:,:)
+!$acc end kernels
+
+!*** PRVS(:,:,:) = PRVS - MYM( ZWORK1 * MZF(PWT) )
+!$acc kernels present_cr(ZTMP2_DEVICE)
+ZTMP2_DEVICE(:,:,:)=ZWORK1(:,:,:) * ZTMP1_DEVICE(:,:,:) !ZWORK1 * MZF(PWT)
+!$acc end kernels
+CALL MYM_DEVICE(ZTMP2_DEVICE,ZTMP3_DEVICE) !MYM( ZWORK1 * MZF(PWT) )
+!$acc kernels present_cr(PRVS)
+ PRVS(:,:,:) = PRVS(:,:,:) - ZTMP2_DEVICE(:,:,:)
+!$acc end kernels
+
+!*** PRWS(:,:,:) = PRWS + MZM( ZWORK2 * MXF(PUT) + ZWORK1 * MYF(PVT) )
+CALL MXF_DEVICE(PUT,ZTMP1_DEVICE) !MXF(PUT)
+CALL MYF_DEVICE(PVT,ZTMP2_DEVICE) !MYF(PVT)
+!$acc kernels present_cr(ZTMP3_DEVICE)
+ZTMP3_DEVICE(:,:,:)= ZWORK2(:,:,:) * ZTMP1_DEVICE(:,:,:)+ ZWORK1(:,:,:) * ZTMP2_DEVICE(:,:,:) !ZWORK2 * MXF(PUT) + ZWORK1 * MYF(PVT)
+!$acc end kernels
+CALL MZM_DEVICE(ZTMP3_DEVICE,ZTMP1_DEVICE) !MZM( ZWORK2 * MXF(PUT) + ZWORK1 * MYF(PVT) )
+!$acc kernels present_cr(PRWS)
+ PRWS(:,:,:) = PRWS(:,:,:) + ZTMP1_DEVICE(:,:,:)
+!$acc end kernels
+#endif
+
+ if ( lbudget_w ) call Budget_store_end( tbudgets(NBUDGET_W), 'COR', prws(:, :, :) )
+ END IF
+
+ if ( lbudget_u ) call Budget_store_end( tbudgets(NBUDGET_U), 'COR', prus(:, :, :) )
+ if ( lbudget_v ) call Budget_store_end( tbudgets(NBUDGET_V), 'COR', prvs(:, :, :) )
+END IF
+!!!!$acc end data
+!
+!-------------------------------------------------------------------------------
+!
+!* 4. COMPUTES THE THETA SOURCE TERM DUE TO THE REFERENCE PRESSURE
+! ------------------------------------------------------------
+!
+IF (LCARTESIAN .OR. LTHINSHELL) THEN
+ ZD1=0.
+ELSE
+ ZD1=1.
+ENDIF
+!
+IF( .NOT.L1D ) THEN
+ IF (.NOT. LOCEAN) THEN
+!
+ IF (KRR > 0) THEN
+ if ( lbudget_th ) call Budget_store_init( tbudgets(NBUDGET_TH), 'PREF', prths(:, :, :) )
+
+
+ ZCPD_OV_RD = XCPD / XRD
+ ZG_OV_CPD = -XG / XCPD
+!
+! stores the specific heat capacity (Cph) in ZWORK1
+
+!$acc kernels present_cr( ZWORK1)
+ ZWORK1(:,:,:) = XCPD + XCPV * PRT(:,:,:,1) ! gas mixing
+!$acc loop seq
+ DO JWATER = 2,1+KRRL ! loop on the liquid components
+ ZWORK1(:,:,:) = ZWORK1(:,:,:) + XCL * PRT(:,:,:,JWATER)
+ END DO
+!$acc loop seq
+ DO JWATER = 2+KRRL,1+KRRL+KRRI ! loop on the solid components
+ ZWORK1(:,:,:) = ZWORK1(:,:,:) + XCI * PRT(:,:,:,JWATER)
+ END DO
+!$acc end kernels
+! computes the source term
+!
+#ifndef MNH_OPENACC
+ PRTHS(:,:,:) = PRTHS(:,:,:) + PRHODJ(:,:,:) &
+ * ( ( XRD + XRV * PRT(:,:,:,1) ) * ZCPD_OV_RD / ZWORK1(:,:,:) - 1. ) &
+ * PTHT(:,:,:)/PEXNREF(:,:,:)*MZF(PWT(:,:,:))*(ZG_OV_CPD/PTHVREF(:,:,:) &
+ -ZD1*4./7.*PEXNREF(:,:,:)/( XRADIUS+MZF(PZZ(:,:,:)) ))
+#else
+CALL MZF_DEVICE(PWT,ZTMP1_DEVICE) !MZF(PWT(:,:,:))
+CALL MZF_DEVICE(PZZ(:,:,:),ZTMP2_DEVICE(:,:,:)) !MZF(PZZ(:,:,:))
+!$acc kernels present_cr(PRTHS)
+PRTHS(:,:,:) = PRTHS(:,:,:) + PRHODJ(:,:,:) &
+ * ( ( XRD + XRV * PRT(:,:,:,1) ) * ZCPD_OV_RD / ZWORK1(:,:,:) - 1. ) &
+ * PTHT(:,:,:)/PEXNREF(:,:,:)*ZTMP1_DEVICE*(ZG_OV_CPD/PTHVREF(:,:,:) &
+ -ZD1*4./7.*PEXNREF(:,:,:)/( XRADIUS+ZTMP2_DEVICE(:,:,:) ))
+!$acc end kernels
+#endif
+ if ( lbudget_th ) call Budget_store_end( tbudgets(NBUDGET_TH), 'PREF', prths(:, :, :) )
+ END IF
+ END IF
+!
+END IF
+!
+!-------------------------------------------------------------------------------
+#ifdef MNH_OPENACC
+CALL MNH_MEM_RELEASE()
+#endif
+!$acc end data
+!$acc update self(PRUS, PRVS, PRWS,PRTHS)
+END SUBROUTINE DYN_SOURCES