From 2fd76d72fd2690ee76815a4beb11f12f1a8d8ed5 Mon Sep 17 00:00:00 2001
From: ESCOBAR Juan <escj@nuwa>
Date: Thu, 21 Nov 2013 00:42:50 +0100
Subject: [PATCH] Juan 20/11/2013: OK advecuvw , manque optimisation

---
 MNH/advec_4th_order_aux.f90 | 113 +++++++++++++++------------
 MNH/advection.f90           |   1 +
 MNH/advecuvw_4th.f90        | 147 +++++++++++++++++++++++++++++++-----
 3 files changed, 195 insertions(+), 66 deletions(-)

diff --git a/MNH/advec_4th_order_aux.f90 b/MNH/advec_4th_order_aux.f90
index 96440f2cd..559928dec 100644
--- a/MNH/advec_4th_order_aux.f90
+++ b/MNH/advec_4th_order_aux.f90
@@ -120,7 +120,10 @@ CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCX ! X direction LBC type
 CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCY ! Y direction LBC type
 !
 REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PMEANX, PMEANY ! fluxes
+!$acc declare pcopy(PMEANX,PMEANY)
+!acc declare present(PMEANX,PMEANY)
 REAL, DIMENSION(:,:,:), INTENT(IN)    :: PFIELDT  ! variable at t
+!$acc declare pcopyin(PFIELDT)
 INTEGER,                INTENT(IN)    :: KGRID    ! C grid localisation
 !
 TYPE(HALO2_ll), OPTIONAL, POINTER :: TPHALO2      ! halo2 for the field at t
@@ -134,20 +137,36 @@ INTEGER:: IIE,IJE        ! End useful area in x,y directions
 !
 INTEGER:: ILUOUT,IRESP   ! for prints
 !
+! JUAN ACC
+LOGICAL                                               :: GWEST , GEAST
+LOGICAL                                               :: GSOUTH , GNORTH
+REAL, DIMENSION(SIZE(PFIELDT,2),SIZE(PFIELDT,3))      :: ZHALO2_WEST,ZHALO2_EAST
+REAL, DIMENSION(SIZE(PFIELDT,1),SIZE(PFIELDT,3))      :: ZHALO2_SOUTH,ZHALO2_NORTH
+!$acc declare create (ZHALO2_WEST,ZHALO2_EAST,ZHALO2_SOUTH,ZHALO2_NORTH)
+!
 !-------------------------------------------------------------------------------
 !
 !*       0.3.     COMPUTES THE DOMAIN DIMENSIONS
 !                 ------------------------------
 !
+!$acc data  
+!
 CALL GET_INDICE_ll(IIB,IJB,IIE,IJE)
 !
+GWEST = LWEST_ll()
+GEAST = LEAST_ll()
+GSOUTH=LSOUTH_ll()
+GNORTH=LNORTH_ll()
+!
 !-------------------------------------------------------------------------------
 !
 !*       0.4.   INITIALIZE THE FIELDS 
 !               ---------------------
 !
-PMEANX(:,:,:) = 0.0
-PMEANY(:,:,:) = 0.0
+!!$!$acc kernels
+!!$PMEANX(:,:,:) = 0.0
+!!$PMEANY(:,:,:) = 0.0
+!!$!$acc end kernels
 !
 !-------------------------------------------------------------------------------
 !
@@ -161,6 +180,11 @@ SELECT CASE ( HLBCX(1) ) ! X direction LBC type: (1) for left side
 !
 CASE ('CYCL')          ! In that case one must have HLBCX(1) == HLBCX(2)
 !
+ZHALO2_WEST(:,:) = TPHALO2%WEST(:,:)
+ZHALO2_EAST(:,:) = TPHALO2%EAST(:,:)
+!$acc update device (ZHALO2_WEST,ZHALO2_EAST)
+!
+!$acc kernels
   IF(NHALO == 1) THEN
     IW=IIB+1
     IE=IIE
@@ -187,33 +211,26 @@ CASE ('CYCL')          ! In that case one must have HLBCX(1) == HLBCX(2)
 !
 !* lateral boundary conditions
   PMEANX(IWF-1,:,:) = (7.0*( PFIELDT(IW-1,:,:)+PFIELDT(IW-2,:,:) ) -  &
-                           ( PFIELDT(IW,:,:)+TPHALO2%WEST(:,:) ) )/12.0
+                           ( PFIELDT(IW,:,:)+ZHALO2_WEST(:,:) ) )/12.0
 !
   PMEANX(IEF+1,:,:) = (7.0*( PFIELDT(IE+1,:,:)+PFIELDT(IE,:,:) ) -  &
-                       ( TPHALO2%EAST(:,:)+PFIELDT(IE-1,:,:) ) )/12.0
+                       ( ZHALO2_EAST(:,:)+PFIELDT(IE-1,:,:) ) )/12.0
 !
 !* inner domain
   PMEANX(IWF:IEF,:,:) = (7.0*( PFIELDT(IW:IE,:,:)+PFIELDT(IW-1:IE-1,:,:) ) -  &
                        ( PFIELDT(IW+1:IE+1,:,:)+PFIELDT(IW-2:IE-2,:,:) ) )/12.0
+!$acc end kernels
 !
-!!$!
-!!$
-!!$  IF(NHALO == 1) THEN
-!!$    PMEANX(IWF-1,:,:) = (7.0*( PFIELDT(IW-1,:,:)+PFIELDT(IW-2,:,:) ) -  &
-!!$                             ( PFIELDT(IW,:,:)+TPHALO2%WEST(:,:) ) )/12.0
-!!$!
-!!$    PMEANX(IEF+1,:,:) = (7.0*( PFIELDT(IE+1,:,:)+PFIELDT(IE,:,:) ) -  &
-!!$                         ( TPHALO2%EAST(:,:)+PFIELDT(IE-1,:,:) ) )/12.0
-!!$  ENDIF
-!!$!
-!!$  PMEANX(IWF:IEF,:,:) = (7.0*( PFIELDT(IW:IE,:,:)+PFIELDT(IW-1:IE-1,:,:) ) -  &
-!!$                       ( PFIELDT(IW+1:IE+1,:,:)+PFIELDT(IW-2:IE-2,:,:) ) )/12.0
-!!$!
 !*       1.2    NON CYCLIC CASE IN THE X DIRECTION 
 !
 CASE ('OPEN','WALL','NEST')
 !
-  IF (LWEST_ll()) THEN
+ZHALO2_WEST(:,:) = TPHALO2%WEST(:,:)
+ZHALO2_EAST(:,:) = TPHALO2%EAST(:,:)
+!$acc update device (ZHALO2_WEST,ZHALO2_EAST)
+!
+!$acc kernels
+  IF (GWEST) THEN
     IF(KGRID == 2) THEN
       IW=IIB+2          ! special case of C grid
     ELSE
@@ -226,7 +243,7 @@ CASE ('OPEN','WALL','NEST')
       IW=IIB
     ENDIF
   ENDIF
-  IF (LEAST_ll() .OR. NHALO == 1) THEN
+  IF (GEAST .OR. NHALO == 1) THEN
 ! T. Maric
 !    IE=IIE-1 ! original
     IE=IIE
@@ -248,7 +265,7 @@ CASE ('OPEN','WALL','NEST')
 !
 !*             Use a second order scheme at the physical border
 !
-  IF (LWEST_ll()) THEN
+  IF (GWEST) THEN
     PMEANX(IWF-1,:,:) = 0.5*( PFIELDT(IW-1,:,:)+PFIELDT(IW-2,:,:) )
     ! T. Maric
     ! PMEANX(1,:,:) = PMEANX(IWF-1,:,:)
@@ -256,20 +273,21 @@ CASE ('OPEN','WALL','NEST')
     !PMEANX(1,:,:) = 0.5*(3.0*PFIELDT(1,:,:) - PFIELDT(2,:,:))
   ELSEIF (NHALO == 1) THEN
     PMEANX(IWF-1,:,:) = (7.0*( PFIELDT(IW-1,:,:)+PFIELDT(IW-2,:,:) ) -  &
-                             ( PFIELDT(IW,:,:)+TPHALO2%WEST(:,:) ) )/12.0
+                             ( PFIELDT(IW,:,:)+ZHALO2_WEST(:,:) ) )/12.0
   ENDIF
 !
-  IF (LEAST_ll()) THEN
+  IF (GEAST) THEN
     PMEANX(IEF+1,:,:) = 0.5*( PFIELDT(IE+1,:,:)+PFIELDT(IE,:,:) )
   ELSEIF (NHALO == 1) THEN
     PMEANX(IEF+1,:,:) = (7.0*( PFIELDT(IE+1,:,:)+PFIELDT(IE,:,:) ) -  &
-                         ( TPHALO2%EAST(:,:)+PFIELDT(IE-1,:,:) ) )/12.0
+                         ( ZHALO2_EAST(:,:)+PFIELDT(IE-1,:,:) ) )/12.0
   ENDIF
 !
 !*             Use a fourth order scheme elsewhere 
 !
   PMEANX(IWF:IEF,:,:) = (7.0*( PFIELDT(IW:IE,:,:)+PFIELDT(IW-1:IE-1,:,:) ) -  &
                        ( PFIELDT(IW+1:IE+1,:,:)+PFIELDT(IW-2:IE-2,:,:) ) )/12.0
+!$acc end kernels
 END SELECT
 !
 !-------------------------------------------------------------------------------
@@ -284,6 +302,12 @@ IF ( .NOT. L2D ) THEN
 !
   CASE ('CYCL')          ! In that case one must have HLBCY(1) == HLBCY(2)
 !
+ZHALO2_SOUTH(:,:) = TPHALO2%SOUTH(:,:) 
+ZHALO2_NORTH(:,:) = TPHALO2%NORTH(:,:)
+!$acc update device (ZHALO2_SOUTH,ZHALO2_NORTH)
+!
+!$acc kernels
+!
 !
     IF(NHALO == 1) THEN
       IS=IJB+1
@@ -311,31 +335,26 @@ IF ( .NOT. L2D ) THEN
 !
 !* lateral boundary conditions
     PMEANY(:,ISF-1,:) = (7.0*( PFIELDT(:,IS-1,:)+PFIELDT(:,IS-2,:) ) -  &
-                        ( PFIELDT(:,IS,:)+TPHALO2%SOUTH(:,:) ) )/12.0
+                        ( PFIELDT(:,IS,:)+ZHALO2_SOUTH(:,:) ) )/12.0
 !
     PMEANY(:,INF+1,:) = (7.0*( PFIELDT(:,IN+1,:)+PFIELDT(:,IN,:) ) -  &
-                        ( TPHALO2%NORTH(:,:)+PFIELDT(:,IN-1,:) ) )/12.0
+                        ( ZHALO2_NORTH(:,:)+PFIELDT(:,IN-1,:) ) )/12.0
 !
 !* inner domain
     PMEANY(:,ISF:INF,:) = (7.0*( PFIELDT(:,IS:IN,:)+PFIELDT(:,IS-1:IN-1,:)) -  &
                          ( PFIELDT(:,IS+1:IN+1,:)+PFIELDT(:,IS-2:IN-2,:) ))/12.0
-!!$!
-!!$    IF(NHALO == 1) THEN
-!!$      PMEANY(:,ISF-1,:) = (7.0*( PFIELDT(:,IS,:)+PFIELDT(:,IS-1,:) ) -  &
-!!$                          ( PFIELDT(:,IS+1,:)+TPHALO2%SOUTH(:,:) ) )/12.0
-!!$!
-!!$      PMEANY(:,ISF+1,:) = (7.0*( PFIELDT(:,IS,:)+PFIELDT(:,IS-1,:) ) -  &
-!!$                          ( TPHALO2%NORTH(:,:)+PFIELDT(:,IS-2,:) ) )/12.0
-!!$    ENDIF
-!!$!
-!!$    PMEANY(:,ISF:INF,:) = (7.0*( PFIELDT(:,IS:IN,:)+PFIELDT(:,IS-1:IN-1,:)) -  &
-!!$                         ( PFIELDT(:,IS+1:IN+1,:)+PFIELDT(:,IS-2:IN-2,:) ))/12.0
-!!$!
+!$acc end kernels
+!
 !*       2.2    NON CYCLIC CASE IN THE Y DIRECTION
 !
   CASE ('OPEN','WALL','NEST')
 !
-    IF (LSOUTH_ll()) THEN
+ZHALO2_SOUTH(:,:) = TPHALO2%SOUTH(:,:)
+ZHALO2_NORTH(:,:) = TPHALO2%NORTH(:,:)
+!$acc update device (ZHALO2_SOUTH,ZHALO2_NORTH)
+!
+!$acc kernels
+    IF (GSOUTH) THEN
       IF(KGRID == 3) THEN
         IS=IJB+2          ! special case of C grid
       ELSE
@@ -348,7 +367,7 @@ IF ( .NOT. L2D ) THEN
         IS=IJB
       ENDIF
     ENDIF
-    IF (LNORTH_ll() .OR. NHALO == 1) THEN
+    IF (GNORTH .OR. NHALO == 1) THEN
 ! T. Maric
 !      IN=IJE-1  ! original
       IN=IJE
@@ -366,38 +385,38 @@ IF ( .NOT. L2D ) THEN
 !
 !*             Use a second order scheme at the physical border
 !
-    IF (LSOUTH_ll()) THEN
+    IF (GSOUTH) THEN
       PMEANY(:,ISF-1,:) = 0.5*( PFIELDT(:,IS-1,:)+PFIELDT(:,IS-2,:) )
       ! T. Maric
       ! PMEANY(:,1,:) = PMEANY(:,ISF-1,:)
       ! extrapolate
       !PMEANY(:,1,:) = 0.5*(3.0*PFIELDT(:,1,:) - PFIELDT(:,2,:))
     ELSEIF (NHALO == 1) THEN
-!!$      PMEANY(:,ISF-1,:) = (7.0*( PFIELDT(:,IS,:)+PFIELDT(:,IS-1,:)) -  &
-!!$                          ( PFIELDT(:,IS+1,:)+TPHALO2%SOUTH(:,:) ))/12.0
        PMEANY(:,ISF-1,:) = (7.0*( PFIELDT(:,IS-1,:)+PFIELDT(:,IS-2,:)) -  &
-                           ( PFIELDT(:,IS,:)+TPHALO2%SOUTH(:,:) ))/12.0
+                           ( PFIELDT(:,IS,:)+ZHALO2_SOUTH(:,:) ))/12.0
     ENDIF
 !
-    IF (LNORTH_ll()) THEN
+    IF (GNORTH) THEN
       PMEANY(:,INF+1,:) = 0.5*( PFIELDT(:,IN+1,:)+PFIELDT(:,IN,:) )
     ELSEIF (NHALO == 1) THEN
-!!$      PMEANY(:,INF+1,:) = (7.0*( PFIELDT(:,IN,:)+PFIELDT(:,IN-1,:)) -  &
-!!$                          ( TPHALO2%NORTH(:,:)+PFIELDT(:,IN-2,:) ))/12.0
        PMEANY(:,INF+1,:) = (7.0*( PFIELDT(:,IN+1,:)+PFIELDT(:,IN,:)) -  &
-                           ( TPHALO2%NORTH(:,:)+PFIELDT(:,IN-1,:) ))/12.0
+                           ( ZHALO2_NORTH(:,:)+PFIELDT(:,IN-1,:) ))/12.0
     ENDIF
 !
 !*             Use a fourth order scheme elsewhere 
 !
     PMEANY(:,ISF:INF,:) = (7.0*( PFIELDT(:,IS:IN,:)+PFIELDT(:,IS-1:IN-1,:)) -  &
                          ( PFIELDT(:,IS+1:IN+1,:)+PFIELDT(:,IS-2:IN-2,:) ))/12.0
+!$acc end kernels
 !
   END SELECT
 ELSE
+!$acc kernels
   PMEANY(:,:,:) = 0.0
+!$acc end kernels
 ENDIF
 !
+!$acc end data
 !-------------------------------------------------------------------------------
 !
 END SUBROUTINE ADVEC_4TH_ORDER_ALGO
diff --git a/MNH/advection.f90 b/MNH/advection.f90
index adfe7e0ac..694756cba 100644
--- a/MNH/advection.f90
+++ b/MNH/advection.f90
@@ -352,6 +352,7 @@ IF(NHALO == 1) THEN
   CALL UPDATE_HALO_ll(TZFIELDS_ll,IINFO_ll)
   CALL CLEANLIST_ll(TZFIELDS_ll)
 END IF
+!$acc update device(ZRUCT,ZRVCT,ZRWCT)
 !
 !-------------------------------------------------------------------------------
 !
diff --git a/MNH/advecuvw_4th.f90 b/MNH/advecuvw_4th.f90
index a8ef52ae1..359e6fb28 100644
--- a/MNH/advecuvw_4th.f90
+++ b/MNH/advecuvw_4th.f90
@@ -171,10 +171,13 @@ CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCY ! Y direction LBC type
 REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PRUCT  ! contravariant
 REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PRVCT  !  components
 REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PRWCT  ! of momentum
+!$acc declare present(PRUCT,PRVCT,PRWCT)
 !
 REAL, DIMENSION(:,:,:),   INTENT(IN) :: PUT, PVT, PWT        ! Variables at t
+!$acc declare pcopyin(PUT,PVT,PWT)
 !
 REAL, DIMENSION(:,:,:),   INTENT(INOUT) :: PRUS, PRVS, PRWS     ! Source terms
+!$acc declare pcopy(PRUS, PRVS, PRWS)
 !
 TYPE(HALO2LIST_ll), POINTER :: TPHALO2LIST ! list for diffusion
 !
@@ -187,13 +190,42 @@ TYPE(HALO2LIST_ll), POINTER :: TZHALO2LIST
 !
 INTEGER :: IGRID ! localisation on the model grid
 REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3)) :: ZMEANX, ZMEANY ! fluxes
+!acc declare create(ZMEANX, ZMEANY)
+!
+REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3)) :: ZTEMP1,ZTEMP2,ZTEMP3,ZTEMP4
+!$acc declare create(ZTEMP1,ZTEMP2,ZTEMP3,ZTEMP4) 
+
+INTEGER                                              :: IIU,IJU,IKU
+INTEGER                                              :: II
+
+!
+#define dxm(PDXM,PA) PDXM(2:IIU,:,:)   = PA(2:IIU,:,:) - PA(1:IIU-1,:,:)       ; PDXM(1,:,:) = PDXM(IIU-2*JPHEXT+1,:,:) ! DXM(PDXM,PA)
+#define mxf(PMXF,PA) PMXF(1:IIU-1,:,:) = 0.5*( PA(2:IIU,:,:)+PA(1:IIU-1,:,:) ) ; PMXF(IIU,:,:) = PMXF(2*JPHEXT,:,:) ! MXF(PMXF,PA)
+#define mxm(PMXM,PA) PMXM(2:IIU,:,:) = 0.5*( PA(2:IIU,:,:)+PA(1:IIU-1,:,:) ) ; PMXM(1,:,:) = PMXM(IIU-2*JPHEXT+1,:,:) ! MXM(PMXM,PA)
+#define dyf(PDYF,PA) PDYF(:,1:IJU-1,:) = PA(:,2:IJU,:) - PA(:,1:IJU-1,:) ; PDYF(:,IJU,:) = PDYF(:,2*JPHEXT,:) ! DYF(PDYF,PA)
+#define dzf(PDZF,PA) PDZF(:,:,1:IKU-1) = PA(:,:,2:IKU) - PA(:,:,1:IKU-1) ; PDZF(:,:,IKU) = -999. ! DZF(PDZF,PA)
+#define mzm4(PMZM4,PA) PMZM4(:,:,3:IKU-1) = (7.0*( PA(:,:,3:IKU-1)+PA(:,:,2:IKU-2) ) - (PA(:,:,4:IKU)+PA(:,:,1:IKU-3) ) )/12.0 ; \
+  PMZM4(:,:,2) = 0.5*( PA(:,:,2)+PA(:,:,1) ) ; PMZM4(:,:,IKU) = 0.5*( PA(:,:,IKU)+PA(:,:,IKU-1) ) ; PMZM4(:,:,1) = -999.
+#define mym(PMYM,PA) PMYM(:,2:IJU,:) = 0.5*( PA(:,2:IJU,:)+PA(:,1:IJU-1,:) ) ; PMYM(:,1,:) = PMYM(:,IJU-2*JPHEXT+1,:) !  MYM(PMYM,PA)
+#define dxf(PDXF,PA) PDXF(1:IIU-1,:,:) = PA(2:IIU,:,:) - PA(1:IIU-1,:,:) ; PDXF(IIU,:,:)    = PDXF(2*JPHEXT,:,:) ! DXF(PDXF,PA)
+#define myf(PMYF,PA) PMYF(:,1:IJU-1,:) = 0.5*( PA(:,1:IJU-1,:)+PA(:,2:IJU,:) ) ; PMYF(:,IJU,:) = PMYF(:,2*JPHEXT,:) ! MYF(PMYF,PA)
+#define dym(PDYM,PA) PDYM(:,2:IJU,:) = PA(:,2:IJU,:) - PA(:,1:IJU-1,:) ; PDYM(:,1,:) = PDYM(:,IJU-2*JPHEXT+1,:) ! DYM(PDYM,PA)
+#define mzm(PMZM,PA) PMZM(:,:,2:IKU) = 0.5*( PA(:,:,2:IKU)+PA(:,:,1:IKU-1) ) ; PMZM(:,:,1)    = -999. !  MZM(PMZM,PA)
+#define mzf(PMZF,PA) PMZF(:,:,1:IKU-1) = 0.5*( PA(:,:,1:IKU-1)+PA(:,:,2:IKU) ) ; PMZF(:,:,IKU) = -999. ! MZF(PMZF,PA)
+#define dzm(PDZM,PA) PDZM(:,:,2:IKU) = PA(:,:,2:IKU) - PA(:,:,1:IKU-1) ; PDZM(:,:,1) = -999. !  DZM(PDZM,PA)
+#define mzf4(PMZF4,PA) PMZF4(:,:,2:IKU-2) = (7.0*( PA(:,:,3:IKU-1)+PA(:,:,2:IKU-2) ) - (PA(:,:,4:IKU)+PA(:,:,1:IKU-3) ) )/12.0 ; \
+  PMZF4(:,:,1) = 0.5*( PA(:,:,2)+PA(:,:,1) ) ; PMZF4(:,:,IKU-1) = 0.5*( PA(:,:,IKU)+PA(:,:,IKU-1) ) ; PMZF4(:,:,IKU) = -999.
 !
 !-------------------------------------------------------------------------------
 !
 !*       1.     COMPUTES THE DOMAIN DIMENSIONS
+!
 !               ------------------------------
+!$acc data  create(ZMEANX, ZMEANY)
 !
 CALL GET_INDICE_ll(IIB,IJB,IIE,IJE)
+CALL GET_DIM_EXT_ll('B',IIU,IJU)     
+IKU=SIZE(PUT,3)
 !
 !-------------------------------------------------------------------------------
 !
@@ -209,16 +241,43 @@ ELSE
   CALL ADVEC_4TH_ORDER_ALGO(HLBCX, HLBCY, PUT, IGRID, ZMEANX, ZMEANY)
 ENDIF
 !
-PRUS(:,:,:) = PRUS(:,:,:)                          &
-             -DXM( MXF(PRUCT(:,:,:))*ZMEANX(:,:,:) ) 
+! pcopy(prus) pcopyin(pruct,ZMEANX) create(ZTEMP1,ZTEMP2,ZTEMP3)
+!!$PRUS(:,:,:) = PRUS(:,:,:)                          &
+!!$             -DXM( MXF(PRUCT(:,:,:))*ZMEANX(:,:,:) ) 
+
+!$acc kernels   pcopyin(ZMEANX)       
+mxf(ZTEMP1,PRUCT)         
+ZTEMP2 = ZTEMP1 * ZMEANX
+dxm(ZTEMP3,ZTEMP2)
+PRUS(:,:,:) = PRUS(:,:,:) -  ZTEMP3
+!$acc end kernels   
+         
 IF (LBUDGET_U)  CALL BUDGET (PRUS,1,'ADVX_BU_RU')
-!
-PRUS(:,:,:) = PRUS(:,:,:)                          &
-             -DYF( MXM(PRVCT(:,:,:))*ZMEANY(:,:,:) ) 
+
+!
+!!$PRUS(:,:,:) = PRUS(:,:,:)                          &
+!!$             -DYF( MXM(PRVCT(:,:,:))*ZMEANY(:,:,:) )
+
+!$acc kernels   pcopyin(ZMEANY) 
+mxm(ZTEMP1,PRVCT)
+ZTEMP2 = ZTEMP1 * ZMEANY
+dyf(ZTEMP3,ZTEMP2)
+PRUS(:,:,:) = PRUS(:,:,:) - ZTEMP3
+!$acc end kernels            
+
 IF (LBUDGET_U)  CALL BUDGET (PRUS,1,'ADVY_BU_RU')
 !
-PRUS(:,:,:) = PRUS(:,:,:)                             &
-             -DZF( MXM(PRWCT(:,:,:))*MZM4(PUT(:,:,:)) )
+!!$PRUS(:,:,:) = PRUS(:,:,:)                             &
+!!$             -DZF( MXM(PRWCT(:,:,:))*MZM4(PUT(:,:,:)) )
+
+!$acc kernels            
+mzm4(ZTEMP1,PUT)
+mxm(ZTEMP2,PRWCT)            
+ZTEMP3 = ZTEMP1 * ZTEMP2
+dzf(ZTEMP4,ZTEMP3)
+PRUS(:,:,:) = PRUS(:,:,:) - ZTEMP4
+!$acc end kernels            
+
 IF (LBUDGET_U)  CALL BUDGET (PRUS,1,'ADVZ_BU_RU')
 !
 !
@@ -231,16 +290,41 @@ ELSE
   CALL ADVEC_4TH_ORDER_ALGO(HLBCX, HLBCY, PVT, IGRID, ZMEANX, ZMEANY)
 ENDIF
 !
-PRVS(:,:,:) = PRVS(:,:,:)                          &
-             -DXF( MYM(PRUCT(:,:,:))*ZMEANX(:,:,:) ) 
+!!$PRVS(:,:,:) = PRVS(:,:,:)                          &
+!!$             -DXF( MYM(PRUCT(:,:,:))*ZMEANX(:,:,:) ) 
+
+!$acc kernels  pcopyin(ZMEANX)          
+mym(ZTEMP1,PRUCT)
+ZTEMP2 = ZTEMP1 * ZMEANX
+dxf(ZTEMP3,ZTEMP2) 
+PRVS(:,:,:) = PRVS(:,:,:) -  ZTEMP3   
+!$acc end kernels   
+                     
 IF (LBUDGET_V)  CALL BUDGET (PRVS,2,'ADVX_BU_RV')
 !
-PRVS(:,:,:) = PRVS(:,:,:)                          &
-             -DYM( MYF(PRVCT(:,:,:))*ZMEANY(:,:,:) )  
+!!$PRVS(:,:,:) = PRVS(:,:,:)                          &
+!!$             -DYM( MYF(PRVCT(:,:,:))*ZMEANY(:,:,:) )
+
+!$acc kernels   pcopyin(ZMEANY)    
+myf(ZTEMP1,PRVCT)
+ZTEMP2 = ZTEMP1 * ZMEANY
+dym(ZTEMP3,ZTEMP2)
+PRVS(:,:,:) = PRVS(:,:,:) - ZTEMP3
+!$acc end kernels   
+
 IF (LBUDGET_V)  CALL BUDGET (PRVS,2,'ADVY_BU_RV')
 !
-PRVS(:,:,:) = PRVS(:,:,:)                             &
-             -DZF( MYM(PRWCT(:,:,:))*MZM4(PVT(:,:,:)) )
+!!$PRVS(:,:,:) = PRVS(:,:,:)                             &
+!!$             -DZF( MYM(PRWCT(:,:,:))*MZM4(PVT(:,:,:)) )
+
+!$acc kernels    
+mym(ZTEMP1,PRWCT)
+mzm4(ZTEMP2,PVT)
+ZTEMP3 = ZTEMP1 * ZTEMP2
+dzf(ZTEMP4,ZTEMP3)
+PRVS(:,:,:) = PRVS(:,:,:) - ZTEMP4
+!$acc end kernels
+
 IF (LBUDGET_V)  CALL BUDGET (PRVS,2,'ADVZ_BU_RV')
 !
 !
@@ -254,18 +338,43 @@ ELSE
   CALL ADVEC_4TH_ORDER_ALGO(HLBCX, HLBCY, PWT, IGRID, ZMEANX, ZMEANY)
 ENDIF
 !
-PRWS(:,:,:) = PRWS(:,:,:)                          &
-             -DXF( MZM(PRUCT(:,:,:))*ZMEANX(:,:,:) ) 
+!!$PRWS(:,:,:) = PRWS(:,:,:)                          &
+!!$             -DXF( MZM(PRUCT(:,:,:))*ZMEANX(:,:,:) ) 
+
+!$acc kernels    pcopyin(ZMEANX) 
+mzm(ZTEMP1,PRUCT)
+ZTEMP2 = ZTEMP1 * ZMEANX
+dxf(ZTEMP3,ZTEMP2)
+PRWS(:,:,:) = PRWS(:,:,:)  - ZTEMP3
+!$acc end kernels
+
 IF (LBUDGET_W)  CALL BUDGET (PRWS,3,'ADVX_BU_RW')
 !
-PRWS(:,:,:) = PRWS(:,:,:)                          &
-             -DYF( MZM(PRVCT(:,:,:))*ZMEANY(:,:,:) ) 
+!!$PRWS(:,:,:) = PRWS(:,:,:)                          &
+!!$             -DYF( MZM(PRVCT(:,:,:))*ZMEANY(:,:,:) )
+
+!$acc kernels   pcopyin(ZMEANY)  
+mzm(ZTEMP1,PRVCT)
+ZTEMP2 = ZTEMP1 * ZMEANY
+dyf(ZTEMP3,ZTEMP2)
+PRWS(:,:,:) = PRWS(:,:,:) - ZTEMP3
+!$acc end kernels
+
 IF (LBUDGET_W)  CALL BUDGET (PRWS,3,'ADVY_BU_RW')
 !
-PRWS(:,:,:) = PRWS(:,:,:)                             &
-             -DZM( MZF(PRWCT(:,:,:))*MZF4(PWT(:,:,:)) )
+!!$PRWS(:,:,:) = PRWS(:,:,:)                             &
+!!$             -DZM( MZF(PRWCT(:,:,:))*MZF4(PWT(:,:,:)) )
+
+!$acc kernels
+mzf(ZTEMP1,PRWCT)
+mzf4(ZTEMP2,PWT)
+ZTEMP1 = ZTEMP1 * ZTEMP2 
+dzm(ZTEMP4,ZTEMP1)
+PRWS(:,:,:) = PRWS(:,:,:) - ZTEMP4
+!$acc end kernels
 IF (LBUDGET_W)  CALL BUDGET (PRWS,3,'ADVZ_BU_RW')
 !
+!$acc end data 
 !-------------------------------------------------------------------------------
 !
 END SUBROUTINE ADVECUVW_4TH
-- 
GitLab