From 01f0b0e044e7ee67243b212cdbfa32bbdeb3e34c Mon Sep 17 00:00:00 2001
From: Philippe WAUTELET <philippe.wautelet@aero.obs-mip.fr>
Date: Mon, 13 Mar 2023 15:55:25 +0100
Subject: [PATCH] Philippe 13/03/2023: retrieve dyn_sources developments from
 ZSOLVER

---
 src/MNH/dyn_sources.f90 | 258 ++++++++++++++++++++++++++++++++++++----
 1 file changed, 235 insertions(+), 23 deletions(-)

diff --git a/src/MNH/dyn_sources.f90 b/src/MNH/dyn_sources.f90
index a8f3c4d69..c87ce1d18 100644
--- a/src/MNH/dyn_sources.f90
+++ b/src/MNH/dyn_sources.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 1994-2021 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 1994-2023 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.
@@ -161,10 +161,20 @@ USE MODD_DYN
 USE MODD_DYN_n,      ONLY: LOCEAN
 !
 use mode_budget,     only: Budget_store_init, Budget_store_end
+#ifdef MNH_OPENACC
+USE MODE_DEVICE
+USE MODE_MNH_ZWORK,      ONLY: MNH_MEM_GET, MNH_MEM_POSITION_PIN, MNH_MEM_RELEASE
+#endif
 USE MODE_MPPDB
 
 USE MODI_GRADIENT_M
+#ifndef MNH_OPENACC
 USE MODI_SHUMAN
+#else
+USE MODI_SHUMAN
+USE MODI_SHUMAN_DEVICE
+#endif
+
 
 IMPLICIT NONE
 !  
@@ -200,21 +210,68 @@ REAL, DIMENSION(:,:,:),   INTENT(INOUT) :: PRTHS            ! Sources of theta
 REAL       ::  ZCPD_OV_RD   ! = CPD / RD
 REAL       ::  ZG_OV_CPD    ! =-XG / XCPD
 INTEGER    ::  JWATER       ! loop index on the different types of water
-INTEGER    ::  IKU          ! array size along the k direction 
+INTEGER    ::  JIU, JJU, JKU
+INTEGER    ::  JI,JJ,JK
 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, ZRVTT, ZRUTT
+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
 !
 !-------------------------------------------------------------------------------
 !
+JIU = size( PUT, 1 )
+JJU = size( PUT, 2 )
+JKU = size( PUT, 3 )
+
+#ifdef MNH_OPENACC
+CALL MNH_MEM_POSITION_PIN( 'DYN_SOURCES' )
+
+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( 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(:,:,:))
-!
-IKU = SIZE(PUT,3)
+#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
 !
 !-------------------------------------------------------------------------------
 !
@@ -228,10 +285,18 @@ IF ((.NOT.L1D).AND.(.NOT.LCARTESIAN) )  THEN
   if ( lbudget_v ) call Budget_store_init( tbudgets(NBUDGET_V), 'CURV', prvs(:, :, :) )
 
   IF ( LTHINSHELL ) THEN           !  THINSHELL approximation
-!
-    ZWORK1(:,:,:) = SPREAD( PCURVX(:,:),DIM=3,NCOPIES=IKU ) / XRADIUS
-    ZWORK2(:,:,:) = SPREAD( PCURVY(:,:),DIM=3,NCOPIES=IKU ) / XRADIUS
-!
+#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 )
@@ -240,12 +305,54 @@ IF ((.NOT.L1D).AND.(.NOT.LCARTESIAN) )  THEN
     - 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=IKU )
-    ZWORK2(:,:,:) = SPREAD( PCURVY(:,:),DIM=3,NCOPIES=IKU )
+    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,&
@@ -258,6 +365,7 @@ IF ((.NOT.L1D).AND.(.NOT.LCARTESIAN) )  THEN
          &  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 )
@@ -268,6 +376,49 @@ IF ((.NOT.L1D).AND.(.NOT.LCARTESIAN) )  THEN
 !
     PRWS(:,:,:) = PRWS                                              &
     +MZM( ( MXF(ZRUT*PUT) + MYF(ZRVT*PVT) ) * ZWORK3 )
+#else
+!$acc kernels present_cr(ZRVTT,ZRUTT)
+    ZRVTT(:,:,:)=ZRVT(:,:,:)*PVT(:,:,:)
+    ZRUTT(:,:,:)=ZRUT(:,:,:)*PUT(:,:,:)
+!$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
+    CALL MXF_DEVICE(ZRUTT,ZTMP1_DEVICE ) !MXF(ZRUT*PUT)
+    CALL MXF_DEVICE(PUT,ZTMP2_DEVICE )   !MXF(PUT)
+!$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
+    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
@@ -285,26 +436,71 @@ END IF
 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(:, :, :) )
-
-  ZWORK3(:,:,:) = SPREAD( PCORIOZ(:,:),DIM=3,NCOPIES=IKU ) * PRHODJ(:,:,:)
-!
+!$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
 
-    ZWORK1(:,:,:) = SPREAD( PCORIOX(:,:),DIM=3,NCOPIES=IKU) * PRHODJ(:,:,:) 
-    ZWORK2(:,:,:) = SPREAD( PCORIOY(:,:),DIM=3,NCOPIES=IKU) * PRHODJ(:,:,:)
-    !
     CALL MPPDB_CHECK3DM("DYN_SOOURCES:CORIOLIS",PRECISION,&
-         & ZWORK1,ZWORK2,ZWORK3 )
+    & 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
+    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
+
+!$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
+
+    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
@@ -335,29 +531,45 @@ IF( .NOT.L1D ) THEN
     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( 'DYN_SOURCES' )
+#endif
+!$acc end data
+!$acc update self(PRUS, PRVS, PRWS,PRTHS)
 !-------------------------------------------------------------------------------
-!
 END SUBROUTINE DYN_SOURCES 
-- 
GitLab