Index: NEMO/branches/2020/dev_r13747_HPC-10_HPDAonline_DiagGPU/tests/DIA_GPU/DIA_GPU.fcm
===================================================================
--- NEMO/branches/2020/dev_r13747_HPC-10_HPDAonline_DiagGPU/tests/DIA_GPU/DIA_GPU.fcm	(revision 14091)
+++ NEMO/branches/2020/dev_r13747_HPC-10_HPDAonline_DiagGPU/tests/DIA_GPU/DIA_GPU.fcm	(revision 14091)
@@ -0,0 +1,1 @@
+bld::tool::fppkeys   key_si3 key_top key_iomput key_mpp_mpi key_gpu
Index: NEMO/branches/2020/dev_r13747_HPC-10_HPDAonline_DiagGPU/tests/DIA_GPU/EXPREF/axis_def_nemo.xml
===================================================================
--- NEMO/branches/2020/dev_r13747_HPC-10_HPDAonline_DiagGPU/tests/DIA_GPU/EXPREF/axis_def_nemo.xml	(revision 14091)
+++ NEMO/branches/2020/dev_r13747_HPC-10_HPDAonline_DiagGPU/tests/DIA_GPU/EXPREF/axis_def_nemo.xml	(revision 14091)
@@ -0,0 +1,1 @@
+link ../../../cfgs/SHARED/axis_def_nemo.xml
Index: NEMO/branches/2020/dev_r13747_HPC-10_HPDAonline_DiagGPU/tests/DIA_GPU/EXPREF/domain_def_nemo.xml
===================================================================
--- NEMO/branches/2020/dev_r13747_HPC-10_HPDAonline_DiagGPU/tests/DIA_GPU/EXPREF/domain_def_nemo.xml	(revision 14091)
+++ NEMO/branches/2020/dev_r13747_HPC-10_HPDAonline_DiagGPU/tests/DIA_GPU/EXPREF/domain_def_nemo.xml	(revision 14091)
@@ -0,0 +1,1 @@
+link ../../../cfgs/SHARED/domain_def_nemo.xml
Index: NEMO/branches/2020/dev_r13747_HPC-10_HPDAonline_DiagGPU/tests/DIA_GPU/EXPREF/field_def_nemo-ice.xml
===================================================================
--- NEMO/branches/2020/dev_r13747_HPC-10_HPDAonline_DiagGPU/tests/DIA_GPU/EXPREF/field_def_nemo-ice.xml	(revision 14091)
+++ NEMO/branches/2020/dev_r13747_HPC-10_HPDAonline_DiagGPU/tests/DIA_GPU/EXPREF/field_def_nemo-ice.xml	(revision 14091)
@@ -0,0 +1,1 @@
+link ../../../cfgs/SHARED/field_def_nemo-ice.xml
Index: NEMO/branches/2020/dev_r13747_HPC-10_HPDAonline_DiagGPU/tests/DIA_GPU/EXPREF/field_def_nemo-innerttrc.xml
===================================================================
--- NEMO/branches/2020/dev_r13747_HPC-10_HPDAonline_DiagGPU/tests/DIA_GPU/EXPREF/field_def_nemo-innerttrc.xml	(revision 14091)
+++ NEMO/branches/2020/dev_r13747_HPC-10_HPDAonline_DiagGPU/tests/DIA_GPU/EXPREF/field_def_nemo-innerttrc.xml	(revision 14091)
@@ -0,0 +1,1 @@
+link ../../../cfgs/SHARED/field_def_nemo-innerttrc.xml
Index: NEMO/branches/2020/dev_r13747_HPC-10_HPDAonline_DiagGPU/tests/DIA_GPU/EXPREF/field_def_nemo-oce.xml
===================================================================
--- NEMO/branches/2020/dev_r13747_HPC-10_HPDAonline_DiagGPU/tests/DIA_GPU/EXPREF/field_def_nemo-oce.xml	(revision 14091)
+++ NEMO/branches/2020/dev_r13747_HPC-10_HPDAonline_DiagGPU/tests/DIA_GPU/EXPREF/field_def_nemo-oce.xml	(revision 14091)
@@ -0,0 +1,1 @@
+link ../../../cfgs/SHARED/field_def_nemo-oce.xml
Index: NEMO/branches/2020/dev_r13747_HPC-10_HPDAonline_DiagGPU/tests/DIA_GPU/EXPREF/field_def_nemo-pisces.xml
===================================================================
--- NEMO/branches/2020/dev_r13747_HPC-10_HPDAonline_DiagGPU/tests/DIA_GPU/EXPREF/field_def_nemo-pisces.xml	(revision 14091)
+++ NEMO/branches/2020/dev_r13747_HPC-10_HPDAonline_DiagGPU/tests/DIA_GPU/EXPREF/field_def_nemo-pisces.xml	(revision 14091)
@@ -0,0 +1,1 @@
+link ../../../cfgs/SHARED/field_def_nemo-pisces.xml
Index: NEMO/branches/2020/dev_r13747_HPC-10_HPDAonline_DiagGPU/tests/DIA_GPU/EXPREF/grid_def_nemo.xml
===================================================================
--- NEMO/branches/2020/dev_r13747_HPC-10_HPDAonline_DiagGPU/tests/DIA_GPU/EXPREF/grid_def_nemo.xml	(revision 14091)
+++ NEMO/branches/2020/dev_r13747_HPC-10_HPDAonline_DiagGPU/tests/DIA_GPU/EXPREF/grid_def_nemo.xml	(revision 14091)
@@ -0,0 +1,1 @@
+link ../../../cfgs/SHARED/grid_def_nemo.xml
Index: NEMO/branches/2020/dev_r13747_HPC-10_HPDAonline_DiagGPU/tests/DIA_GPU/EXPREF/namelist_ice_ref
===================================================================
--- NEMO/branches/2020/dev_r13747_HPC-10_HPDAonline_DiagGPU/tests/DIA_GPU/EXPREF/namelist_ice_ref	(revision 14091)
+++ NEMO/branches/2020/dev_r13747_HPC-10_HPDAonline_DiagGPU/tests/DIA_GPU/EXPREF/namelist_ice_ref	(revision 14091)
@@ -0,0 +1,1 @@
+link ../../../cfgs/SHARED/namelist_ice_ref
Index: NEMO/branches/2020/dev_r13747_HPC-10_HPDAonline_DiagGPU/tests/DIA_GPU/EXPREF/namelist_pisces_ref
===================================================================
--- NEMO/branches/2020/dev_r13747_HPC-10_HPDAonline_DiagGPU/tests/DIA_GPU/EXPREF/namelist_pisces_ref	(revision 14091)
+++ NEMO/branches/2020/dev_r13747_HPC-10_HPDAonline_DiagGPU/tests/DIA_GPU/EXPREF/namelist_pisces_ref	(revision 14091)
@@ -0,0 +1,1 @@
+link ../../../cfgs/SHARED/namelist_pisces_ref
Index: NEMO/branches/2020/dev_r13747_HPC-10_HPDAonline_DiagGPU/tests/DIA_GPU/EXPREF/namelist_ref
===================================================================
--- NEMO/branches/2020/dev_r13747_HPC-10_HPDAonline_DiagGPU/tests/DIA_GPU/EXPREF/namelist_ref	(revision 14091)
+++ NEMO/branches/2020/dev_r13747_HPC-10_HPDAonline_DiagGPU/tests/DIA_GPU/EXPREF/namelist_ref	(revision 14091)
@@ -0,0 +1,1 @@
+link ../../../cfgs/SHARED/namelist_ref
Index: NEMO/branches/2020/dev_r13747_HPC-10_HPDAonline_DiagGPU/tests/DIA_GPU/EXPREF/namelist_top_ref
===================================================================
--- NEMO/branches/2020/dev_r13747_HPC-10_HPDAonline_DiagGPU/tests/DIA_GPU/EXPREF/namelist_top_ref	(revision 14091)
+++ NEMO/branches/2020/dev_r13747_HPC-10_HPDAonline_DiagGPU/tests/DIA_GPU/EXPREF/namelist_top_ref	(revision 14091)
@@ -0,0 +1,1 @@
+link ../../../cfgs/SHARED/namelist_top_ref
Index: NEMO/branches/2020/dev_r13747_HPC-10_HPDAonline_DiagGPU/tests/DIA_GPU/EXPREF/namelist_trc_ref
===================================================================
--- NEMO/branches/2020/dev_r13747_HPC-10_HPDAonline_DiagGPU/tests/DIA_GPU/EXPREF/namelist_trc_ref	(revision 14091)
+++ NEMO/branches/2020/dev_r13747_HPC-10_HPDAonline_DiagGPU/tests/DIA_GPU/EXPREF/namelist_trc_ref	(revision 14091)
@@ -0,0 +1,1 @@
+link ../../../cfgs/SHARED/namelist_trc_ref
Index: NEMO/branches/2020/dev_r13747_HPC-10_HPDAonline_DiagGPU/tests/DIA_GPU/MY_SRC/cuda_fortran.F90
===================================================================
--- NEMO/branches/2020/dev_r13747_HPC-10_HPDAonline_DiagGPU/tests/DIA_GPU/MY_SRC/cuda_fortran.F90	(revision 14091)
+++ NEMO/branches/2020/dev_r13747_HPC-10_HPDAonline_DiagGPU/tests/DIA_GPU/MY_SRC/cuda_fortran.F90	(revision 14091)
@@ -0,0 +1,118 @@
+#ifdef key_gpu
+MODULE cuda_fortran
+    USE cudafor
+    USE dom_oce
+    USE lib_mpp
+    CONTAINS
+      ATTRIBUTES(global) &
+      SUBROUTINE dia_hsb_kernel(surf, e3t, surf_ini, e3t_ini, ts, hc_loc_ini, sc_loc_ini,     &
+                            & tmask, tmask_ini, zwrkv, zwrkh, zwrks, zwrk, jpi, jpj, jpk, jpt, Kmm)
+            !
+            INTEGER, VALUE :: jpi, jpj, jpk, jpt                                                                         ! MPI process sub-domain and stream tile
+            INTEGER, VALUE :: Kmm                                                                                   !ocean time level indice
+            REAL(8)        :: surf(jpi, jpj), e3t(jpi, jpj, jpk, jpt), surf_ini(jpi, jpj), e3t_ini(jpi, jpj, jpk), ts(jpi, jpj, jpk, 2, jpt)
+            REAL(8)        :: hc_loc_ini(jpi, jpj, jpk), sc_loc_ini(jpi, jpj, jpk), tmask(jpi, jpj, jpk), tmask_ini(jpi, jpj, jpk)
+            REAL(8)        :: zwrkv(jpi, jpj, jpk), zwrkh(jpi, jpj, jpk), zwrks(jpi, jpj, jpk), zwrk(jpi, jpj, jpk)
+            !
+
+            INTEGER        :: i, j, k                       ! dummy indexes
+
+            !
+            i = blockDim%x * (blockIdx%x -1) + threadIdx%x
+            j = blockDim%y * (blockIdx%y -1) + threadIdx%y
+            k = blockDim%z * (blockIdx%z -1) + threadIdx%z
+            !
+            IF ( (i .le. jpi ) .AND. (j .le. jpj) .AND. (k .le. jpk-1) ) THEN
+                zwrkv(i, j, k) =   surf    (i,j) * e3t    (i,j,k,Kmm) * tmask    (i,j,k)      &
+                &                - surf_ini(i,j) * e3t_ini(i,j,k)     * tmask_ini(i,j,k)
+                
+                zwrkh(i, j, k) = ( surf(i, j) * e3t(i, j, k, Kmm) * ts(i, j, k, 1, Kmm) - surf_ini(i, j) * hc_loc_ini(i, j, k) )
+                
+                zwrks(i, j, k) = ( surf(i, j) * e3t(i, j, k, Kmm) * ts(i, j, k, 2, Kmm) - surf_ini(i, j) * sc_loc_ini(i, j, k) )
+                
+                zwrk(i, j, k)  =   surf(i, j) * e3t(i, j, k, Kmm) * tmask(i, j, k)
+            END IF
+      END SUBROUTINE dia_hsb_kernel
+
+!     ATTRIBUTES(global) &
+!     SUBROUTINE dia_hsb_kernel1d(surf, e3t_n, surf_ini, e3t_ini, &
+!        &  tsn, hc_loc_ini, sc_loc_ini, tmask, zwrkv, zwrkh, zwrks, zwrk, jpi, jpj, jpk)
+!        IMPLICIT NONE
+!        REAL(kind=8) :: surf(:, :), e3t_n(:, :, :), surf_ini(:, :), e3t_ini(:, :, :), tsn(:, :, :, :), &
+!                            & hc_loc_ini(:, :, :), sc_loc_ini(:, :, :), tmask(:, :, :), zwrkv(:, :, :),      &
+!                            & zwrkh(:, :, :), zwrks(:, :, :), zwrk(:, :, :)
+!        REAL(kind=8), SHARED       :: sdata(*)
+!        INTEGER, VALUE             :: jpi, jpj, jpk
+!        INTEGER                    :: i, si, ti, globsize
+!
+!        globsize = jpi*jpj*jpk
+!
+!        i  = blockDim%x * (blockIdx%x-1) + threadIdx%x
+!        ti = threadIdx%x
+!        si = MOD(i, jpk)+1 !Shared memory indexing
+!
+!        !Compute volume (zwrk), volume variation (zwrkv), heat deviation (zwrkh)
+!        !and salinity deviation (zwrks)
+!        IF ( i .le. globsize ) THEN
+!            sdata(ti) = surf(si) * e3t_n(i)
+!            zwrkv(i) = ( sdata(ti) - surf_ini(si) * e3t_ini(i) ) * tmask(i) * surf(si)
+!            zwrkh(i) = ( sdata(ti) * tsn(i) - surf_ini(si) * hc_loc_ini(i) ) * tmask(i) * surf(si)
+!            zwrks(i) = ( sdata(ti) * tsn(globsize + i) - surf_ini(si) * sc_loc_ini(i) ) &
+!                    * tmask(i) * surf(si)
+!            zwrk(i) = sdata(ti) * tmask(i) * surf(si)
+!        END IF
+!     END SUBROUTINE dia_hsb_kernel1d
+
+
+       ATTRIBUTES(global) &
+       SUBROUTINE filter_cuda(ptab, mask, jpi, jpj, jpk)
+           !
+           REAL(kind=8)         :: ptab(jpi,jpj,jpk) , mask(jpi,jpj)
+           !
+           INTEGER, VALUE  :: jpi, jpj, jpk                 ! MPI process sub-domain
+           INTEGER         :: i, j, k                 ! dummy indexes
+           !
+           i = blockDim%x * (blockIdx%x -1) + threadIdx%x
+           j = blockDim%y * (blockIdx%y -1) + threadIdx%y
+           k = blockDim%z * (blockIdx%z -1) + threadIdx%z
+           !
+           tile = 1
+           IF ( (i .le. jpi ) .AND. (j .le. jpj) .AND. (k .le. jpk-1) ) THEN
+               ptab(i, j, k)  = ptab(i, j, k) * mask(i, j)
+           END IF
+       END SUBROUTINE filter_cuda
+
+       ATTRIBUTES(global) SUBROUTINE array3dto1d(d_inp,d_out,ipi, ipj, ipk)
+           !
+           REAL(kind=8)                 :: d_inp(:,:,:)
+           REAL(kind=8), intent(out)    :: d_out(:)
+           !
+           INTEGER, VALUE  :: ipi, ipj, ipk                 ! MPI process sub-domain
+           INTEGER         :: i, j, k                       ! dummy indexes
+           !
+           i = blockDim%x * (blockIdx%x -1) + threadIdx%x
+           j = blockDim%y * (blockIdx%y -1) + threadIdx%y
+           k = blockDim%z * (blockIdx%z -1) + threadIdx%z
+           !
+           IF ( (i .le. ipi ) .AND. (j .le. ipj) .AND. (k .le. ipk) ) THEN
+               d_out(i + ipi*(j-1) + ipi*ipj*(k-1) )  = d_inp(i, j, k)
+           END IF
+       END SUBROUTINE array3dto1d
+
+       !Knuth's trick
+       ATTRIBUTES(device) &
+       SUBROUTINE DDPDD_d(ydda, yddb)
+           IMPLICIT NONE
+           COMPLEX(kind=8), INTENT(in   ) :: ydda              !Scalar to add
+           COMPLEX(kind=8), INTENT(inout) :: yddb              !Total sum
+           REAL(kind=8)                   :: zerr, zt1, zt2
+       
+           zt1  = REAL(ydda) + REAL(yddb)
+           zerr = zt1 - REAL(ydda)
+           zt2  = ( (REAL(yddb) - zerr) + (REAL(ydda) - (zt1 - zerr)) ) &
+                 & + AIMAG(ydda) + AIMAG(yddb)
+           !The result is t1 + t2 after normalization
+           yddb = CMPLX( zt1 + zt2, zt2 - ((zt1 + zt2) - zt1), 8 )
+       END SUBROUTINE DDPDD_d       
+END MODULE
+#endif
Index: NEMO/branches/2020/dev_r13747_HPC-10_HPDAonline_DiagGPU/tests/DIA_GPU/MY_SRC/diahsb.F90
===================================================================
--- NEMO/branches/2020/dev_r13747_HPC-10_HPDAonline_DiagGPU/tests/DIA_GPU/MY_SRC/diahsb.F90	(revision 14091)
+++ NEMO/branches/2020/dev_r13747_HPC-10_HPDAonline_DiagGPU/tests/DIA_GPU/MY_SRC/diahsb.F90	(revision 14091)
@@ -0,0 +1,438 @@
+#ifdef key_gpu
+#include "diahsb_gpu.h90"
+#else
+MODULE diahsb
+   !!======================================================================
+   !!                       ***  MODULE  diahsb  ***
+   !! Ocean diagnostics: Heat, salt and volume budgets
+   !!======================================================================
+   !! History :  3.3  ! 2010-09  (M. Leclair)  Original code 
+   !!                 ! 2012-10  (C. Rousset)  add iom_put
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   dia_hsb       : Diagnose the conservation of ocean heat and salt contents, and volume
+   !!   dia_hsb_rst   : Read or write DIA file in restart file
+   !!   dia_hsb_init  : Initialization of the conservation diagnostic
+   !!----------------------------------------------------------------------
+   USE oce            ! ocean dynamics and tracers
+   USE dom_oce        ! ocean space and time domain
+   USE phycst         ! physical constants
+   USE sbc_oce        ! surface thermohaline fluxes
+   USE isf_oce        ! ice shelf fluxes
+   USE sbcrnf         ! river runoff
+   USE domvvl         ! vertical scale factors
+   USE traqsr         ! penetrative solar radiation
+   USE trabbc         ! bottom boundary condition 
+   USE trabbc         ! bottom boundary condition
+   USE restart        ! ocean restart
+   USE bdy_oce , ONLY : ln_bdy
+   !
+   USE iom            ! I/O manager
+   USE in_out_manager ! I/O manager
+   USE lib_fortran    ! glob_sum
+   USE lib_mpp        ! distributed memory computing library
+   USE timing         ! preformance summary
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   dia_hsb        ! routine called by step.F90
+   PUBLIC   dia_hsb_init   ! routine called by nemogcm.F90
+
+   LOGICAL, PUBLIC ::   ln_diahsb   !: check the heat and salt budgets
+
+   REAL(wp) ::   surf_tot              ! ocean surface
+   REAL(wp) ::   frc_t, frc_s, frc_v   ! global forcing trends
+   REAL(wp) ::   frc_wn_t, frc_wn_s    ! global forcing trends
+   !
+   REAL(wp), DIMENSION(:,:)  , ALLOCATABLE ::   surf 
+   REAL(wp), DIMENSION(:,:)  , ALLOCATABLE ::   surf_ini      , ssh_ini          !
+   REAL(wp), DIMENSION(:,:)  , ALLOCATABLE ::   ssh_hc_loc_ini, ssh_sc_loc_ini   !
+   REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   hc_loc_ini, sc_loc_ini, e3t_ini  !
+   REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   tmask_ini
+
+   !! * Substitutions
+#  include "domzgr_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+   !! $Id$
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE dia_hsb( kt, Kbb, Kmm )
+      !!---------------------------------------------------------------------------
+      !!                  ***  ROUTINE dia_hsb  ***
+      !!     
+      !! ** Purpose: Compute the ocean global heat content, salt content and volume conservation
+      !!	
+      !! ** Method : - Compute the deviation of heat content, salt content and volume
+      !!	            at the current time step from their values at nit000
+      !!	            - Compute the contribution of forcing and remove it from these deviations
+      !!
+      !!---------------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt         ! ocean time-step index
+      INTEGER, INTENT(in) ::   Kbb, Kmm   ! ocean time level indices
+      !
+      INTEGER    ::   ji, jj, jk                  ! dummy loop indice
+      REAL(wp)   ::   zdiff_hc    , zdiff_sc      ! heat and salt content variations
+      REAL(wp)   ::   zdiff_hc1   , zdiff_sc1     !  -         -     -        - 
+      REAL(wp)   ::   zdiff_v1    , zdiff_v2      ! volume variation
+      REAL(wp)   ::   zerr_hc1    , zerr_sc1      ! heat and salt content misfit
+      REAL(wp)   ::   zvol_tot                    ! volume
+      REAL(wp)   ::   z_frc_trd_t , z_frc_trd_s   !    -     -
+      REAL(wp)   ::   z_frc_trd_v                 !    -     -
+      REAL(wp)   ::   z_wn_trd_t , z_wn_trd_s     !    -     -
+      REAL(wp)   ::   z_ssh_hc , z_ssh_sc         !    -     -
+      REAL(wp), DIMENSION(jpi,jpj)       ::   z2d0, z2d1   ! 2D workspace
+      REAL(wp), DIMENSION(jpi,jpj,jpkm1) ::   zwrk         ! 3D workspace
+      !!---------------------------------------------------------------------------
+      IF( ln_timing )   CALL timing_start('dia_hsb')      
+      !
+      ts(:,:,:,1,Kmm) = ts(:,:,:,1,Kmm) * tmask(:,:,:) ; ts(:,:,:,1,Kbb) = ts(:,:,:,1,Kbb) * tmask(:,:,:) ;
+      ts(:,:,:,2,Kmm) = ts(:,:,:,2,Kmm) * tmask(:,:,:) ; ts(:,:,:,2,Kbb) = ts(:,:,:,2,Kbb) * tmask(:,:,:) ;
+      ! ------------------------- !
+      ! 1 - Trends due to forcing !
+      ! ------------------------- !
+      z_frc_trd_v = r1_rho0 * glob_sum( 'diahsb', - ( emp(:,:) - rnf(:,:) + fwfisf_cav(:,:) + fwfisf_par(:,:) ) * surf(:,:) )   ! volume fluxes
+      z_frc_trd_t =           glob_sum( 'diahsb', sbc_tsc(:,:,jp_tem) * surf(:,:) )                       ! heat fluxes
+      z_frc_trd_s =           glob_sum( 'diahsb', sbc_tsc(:,:,jp_sal) * surf(:,:) )                       ! salt fluxes
+      !                    !  Add runoff    heat & salt input
+      IF( ln_rnf    )   z_frc_trd_t = z_frc_trd_t + glob_sum( 'diahsb', rnf_tsc(:,:,jp_tem) * surf(:,:) )
+      IF( ln_rnf_sal)   z_frc_trd_s = z_frc_trd_s + glob_sum( 'diahsb', rnf_tsc(:,:,jp_sal) * surf(:,:) )
+      !                    ! Add ice shelf heat & salt input
+      IF( ln_isf    )   z_frc_trd_t = z_frc_trd_t &
+         &                          + glob_sum( 'diahsb', ( risf_cav_tsc(:,:,jp_tem) + risf_par_tsc(:,:,jp_tem) ) * surf(:,:) )
+      !                    ! Add penetrative solar radiation
+      IF( ln_traqsr )   z_frc_trd_t = z_frc_trd_t + r1_rho0_rcp * glob_sum( 'diahsb', qsr     (:,:) * surf(:,:) )
+      !                    ! Add geothermal heat flux
+      IF( ln_trabbc )   z_frc_trd_t = z_frc_trd_t +               glob_sum( 'diahsb', qgh_trd0(:,:) * surf(:,:) )
+      !
+      IF( ln_linssh ) THEN
+         IF( ln_isfcav ) THEN
+            DO ji=1,jpi
+               DO jj=1,jpj
+                  z2d0(ji,jj) = surf(ji,jj) * ww(ji,jj,mikt(ji,jj)) * ts(ji,jj,mikt(ji,jj),jp_tem,Kbb)
+                  z2d1(ji,jj) = surf(ji,jj) * ww(ji,jj,mikt(ji,jj)) * ts(ji,jj,mikt(ji,jj),jp_sal,Kbb)
+               END DO
+            END DO
+         ELSE
+            z2d0(:,:) = surf(:,:) * ww(:,:,1) * ts(:,:,1,jp_tem,Kbb)
+            z2d1(:,:) = surf(:,:) * ww(:,:,1) * ts(:,:,1,jp_sal,Kbb)
+         END IF
+         z_wn_trd_t = - glob_sum( 'diahsb', z2d0 ) 
+         z_wn_trd_s = - glob_sum( 'diahsb', z2d1 )
+      ENDIF
+
+      frc_v = frc_v + z_frc_trd_v * rn_Dt
+      frc_t = frc_t + z_frc_trd_t * rn_Dt
+      frc_s = frc_s + z_frc_trd_s * rn_Dt
+      !                                          ! Advection flux through fixed surface (z=0)
+      IF( ln_linssh ) THEN
+         frc_wn_t = frc_wn_t + z_wn_trd_t * rn_Dt
+         frc_wn_s = frc_wn_s + z_wn_trd_s * rn_Dt
+      ENDIF
+
+      ! ------------------------ !
+      ! 2 -  Content variations  !
+      ! ------------------------ !
+      ! glob_sum_full is needed because you keep the full interior domain to compute the sum (iscpl)
+
+      !                    ! volume variation (calculated with ssh)
+      zdiff_v1 = glob_sum_full( 'diahsb', surf(:,:)*ssh(:,:,Kmm) - surf_ini(:,:)*ssh_ini(:,:) )
+
+      !                    ! heat & salt content variation (associated with ssh)
+      IF( ln_linssh ) THEN       ! linear free surface case
+         IF( ln_isfcav ) THEN          ! ISF case
+            DO ji = 1, jpi
+               DO jj = 1, jpj
+                  z2d0(ji,jj) = surf(ji,jj) * ( ts(ji,jj,mikt(ji,jj),jp_tem,Kmm) * ssh(ji,jj,Kmm) - ssh_hc_loc_ini(ji,jj) ) 
+                  z2d1(ji,jj) = surf(ji,jj) * ( ts(ji,jj,mikt(ji,jj),jp_sal,Kmm) * ssh(ji,jj,Kmm) - ssh_sc_loc_ini(ji,jj) ) 
+               END DO
+            END DO
+         ELSE                          ! no under ice-shelf seas
+            z2d0(:,:) = surf(:,:) * ( ts(:,:,1,jp_tem,Kmm) * ssh(:,:,Kmm) - ssh_hc_loc_ini(:,:) ) 
+            z2d1(:,:) = surf(:,:) * ( ts(:,:,1,jp_sal,Kmm) * ssh(:,:,Kmm) - ssh_sc_loc_ini(:,:) ) 
+         END IF
+         z_ssh_hc = glob_sum_full( 'diahsb', z2d0 ) 
+         z_ssh_sc = glob_sum_full( 'diahsb', z2d1 ) 
+      ENDIF
+      !
+      DO jk = 1, jpkm1           ! volume variation (calculated with scale factors)
+         zwrk(:,:,jk) =   surf    (:,:) * e3t    (:,:,jk,Kmm)*tmask    (:,:,jk)   &
+            &           - surf_ini(:,:) * e3t_ini(:,:,jk    )*tmask_ini(:,:,jk)
+      END DO
+      zdiff_v2 = glob_sum_full( 'diahsb', zwrk(:,:,:) )     ! glob_sum_full needed as tmask and tmask_ini could be different
+      DO jk = 1, jpkm1           ! heat content variation
+         zwrk(:,:,jk) = ( surf    (:,:) * e3t(:,:,jk,Kmm)*ts(:,:,jk,jp_tem,Kmm)   &
+            &           - surf_ini(:,:) *         hc_loc_ini(:,:,jk) )
+      END DO
+      zdiff_hc = glob_sum_full( 'diahsb', zwrk(:,:,:) )
+      DO jk = 1, jpkm1           ! salt content variation
+         zwrk(:,:,jk) = ( surf    (:,:) * e3t(:,:,jk,Kmm)*ts(:,:,jk,jp_sal,Kmm)   &
+            &           - surf_ini(:,:) *         sc_loc_ini(:,:,jk) )
+      END DO
+      zdiff_sc = glob_sum_full( 'diahsb', zwrk(:,:,:) )
+
+      ! ------------------------ !
+      ! 3 -  Drifts              !
+      ! ------------------------ !
+      zdiff_v1 = zdiff_v1 - frc_v
+      IF( .NOT.ln_linssh )   zdiff_v2 = zdiff_v2 - frc_v
+      zdiff_hc = zdiff_hc - frc_t
+      zdiff_sc = zdiff_sc - frc_s
+      IF( ln_linssh ) THEN
+         zdiff_hc1 = zdiff_hc + z_ssh_hc 
+         zdiff_sc1 = zdiff_sc + z_ssh_sc
+         zerr_hc1  = z_ssh_hc - frc_wn_t
+         zerr_sc1  = z_ssh_sc - frc_wn_s
+      ENDIF
+
+      ! ----------------------- !
+      ! 4 - Diagnostics writing !
+      ! ----------------------- !
+      DO jk = 1, jpkm1           ! total ocean volume (calculated with scale factors)
+         zwrk(:,:,jk) = surf(:,:) * e3t(:,:,jk,Kmm) * tmask(:,:,jk)
+      END DO
+      zvol_tot = glob_sum( 'diahsb', zwrk(:,:,:) )
+
+!!gm to be added ?
+!      IF( ln_linssh ) THEN            ! fixed volume, add the ssh contribution
+!        zvol_tot = zvol_tot + glob_sum( 'diahsb', surf(:,:) * ssh(:,:,Kmm) )
+!      ENDIF
+!!gm end
+
+      CALL iom_put(   'bgfrcvol' , frc_v    * 1.e-9    )              ! vol - surface forcing (km3) 
+      CALL iom_put(   'bgfrctem' , frc_t    * rho0 * rcp * 1.e-20 )   ! hc  - surface forcing (1.e20 J) 
+      CALL iom_put(   'bgfrchfx' , frc_t    * rho0 * rcp /  &         ! hc  - surface forcing (W/m2) 
+         &                       ( surf_tot * kt * rn_Dt )        )
+      CALL iom_put(   'bgfrcsal' , frc_s    * 1.e-9    )              ! sc  - surface forcing (psu*km3) 
+
+      IF( .NOT. ln_linssh ) THEN
+         CALL iom_put( 'bgtemper' , zdiff_hc / zvol_tot )              ! Temperature drift     (C) 
+         CALL iom_put( 'bgsaline' , zdiff_sc / zvol_tot )              ! Salinity    drift     (PSU)
+         CALL iom_put( 'bgheatco' , zdiff_hc * 1.e-20 * rho0 * rcp )   ! Heat content drift    (1.e20 J) 
+         CALL iom_put( 'bgheatfx' , zdiff_hc * rho0 * rcp /  &         ! Heat flux drift       (W/m2) 
+            &                       ( surf_tot * kt * rn_Dt )        )
+         CALL iom_put( 'bgsaltco' , zdiff_sc * 1.e-9    )              ! Salt content drift    (psu*km3)
+         CALL iom_put( 'bgvolssh' , zdiff_v1 * 1.e-9    )              ! volume ssh drift      (km3)  
+         CALL iom_put( 'bgvole3t' , zdiff_v2 * 1.e-9    )              ! volume e3t drift      (km3)  
+         !
+         IF( kt == nitend .AND. lwp ) THEN
+            WRITE(numout,*)
+            WRITE(numout,*) 'dia_hsb : last time step hsb diagnostics: at it= ', kt,' date= ', ndastp
+            WRITE(numout,*) '~~~~~~~'
+            WRITE(numout,*) '   Temperature drift = ', zdiff_hc / zvol_tot, ' C'
+            WRITE(numout,*) '   Salinity    drift = ', zdiff_sc / zvol_tot, ' PSU'
+            WRITE(numout,*) '   volume ssh  drift = ', zdiff_v1 * 1.e-9   , ' km^3'
+            WRITE(numout,*) '   volume e3t  drift = ', zdiff_v2 * 1.e-9   , ' km^3'
+         ENDIF
+         !
+      ELSE
+         CALL iom_put( 'bgtemper' , zdiff_hc1 / zvol_tot)              ! Heat content drift    (C) 
+         CALL iom_put( 'bgsaline' , zdiff_sc1 / zvol_tot)              ! Salt content drift    (PSU)
+         CALL iom_put( 'bgheatco' , zdiff_hc1 * 1.e-20 * rho0 * rcp )  ! Heat content drift    (1.e20 J) 
+         CALL iom_put( 'bgheatfx' , zdiff_hc1 * rho0 * rcp /  &        ! Heat flux drift       (W/m2) 
+            &                       ( surf_tot * kt * rn_Dt )         )
+         CALL iom_put( 'bgsaltco' , zdiff_sc1 * 1.e-9    )             ! Salt content drift    (psu*km3)
+         CALL iom_put( 'bgvolssh' , zdiff_v1 * 1.e-9    )              ! volume ssh drift      (km3)  
+         CALL iom_put( 'bgmistem' , zerr_hc1 / zvol_tot )              ! hc  - error due to free surface (C)
+         CALL iom_put( 'bgmissal' , zerr_sc1 / zvol_tot )              ! sc  - error due to free surface (psu)
+      ENDIF
+      !
+      IF( lrst_oce )   CALL dia_hsb_rst( kt, Kmm, 'WRITE' )
+      !
+      IF( ln_timing )   CALL timing_stop('dia_hsb')
+      !
+   END SUBROUTINE dia_hsb
+
+
+   SUBROUTINE dia_hsb_rst( kt, Kmm, cdrw )
+      !!---------------------------------------------------------------------
+      !!                   ***  ROUTINE dia_hsb_rst  ***
+      !!                     
+      !! ** Purpose :   Read or write DIA file in restart file
+      !!
+      !! ** Method  :   use of IOM library
+      !!----------------------------------------------------------------------
+      INTEGER         , INTENT(in) ::   kt     ! ocean time-step
+      INTEGER         , INTENT(in) ::   Kmm    ! ocean time level index
+      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
+      !
+      INTEGER ::   ji, jj, jk   ! dummy loop indices
+      !!----------------------------------------------------------------------
+      !
+      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise 
+         IF( ln_rstart ) THEN                   !* Read the restart file
+            !
+            IF(lwp) WRITE(numout,*)
+            IF(lwp) WRITE(numout,*) '   dia_hsb_rst : read hsb restart at it= ', kt,' date= ', ndastp
+            IF(lwp) WRITE(numout,*)
+            CALL iom_get( numror, 'frc_v', frc_v, ldxios = lrxios )
+            CALL iom_get( numror, 'frc_t', frc_t, ldxios = lrxios )
+            CALL iom_get( numror, 'frc_s', frc_s, ldxios = lrxios )
+            IF( ln_linssh ) THEN
+               CALL iom_get( numror, 'frc_wn_t', frc_wn_t, ldxios = lrxios )
+               CALL iom_get( numror, 'frc_wn_s', frc_wn_s, ldxios = lrxios )
+            ENDIF
+            CALL iom_get( numror, jpdom_auto, 'surf_ini'  , surf_ini  , ldxios = lrxios ) ! ice sheet coupling
+            CALL iom_get( numror, jpdom_auto, 'ssh_ini'   , ssh_ini   , ldxios = lrxios )
+            CALL iom_get( numror, jpdom_auto, 'e3t_ini'   , e3t_ini   , ldxios = lrxios )
+            CALL iom_get( numror, jpdom_auto, 'tmask_ini' , tmask_ini , ldxios = lrxios )
+            CALL iom_get( numror, jpdom_auto, 'hc_loc_ini', hc_loc_ini, ldxios = lrxios )
+            CALL iom_get( numror, jpdom_auto, 'sc_loc_ini', sc_loc_ini, ldxios = lrxios )
+            IF( ln_linssh ) THEN
+               CALL iom_get( numror, jpdom_auto, 'ssh_hc_loc_ini', ssh_hc_loc_ini, ldxios = lrxios )
+               CALL iom_get( numror, jpdom_auto, 'ssh_sc_loc_ini', ssh_sc_loc_ini, ldxios = lrxios )
+            ENDIF
+         ELSE
+            IF(lwp) WRITE(numout,*)
+            IF(lwp) WRITE(numout,*) '   dia_hsb_rst : initialise hsb at initial state '
+            IF(lwp) WRITE(numout,*)
+            surf_ini(:,:) = e1e2t(:,:) * tmask_i(:,:)         ! initial ocean surface
+            ssh_ini(:,:) = ssh(:,:,Kmm)                          ! initial ssh
+            DO jk = 1, jpk
+              ! if ice sheet/oceqn coupling, need to mask ini variables here (mask could change at the next NEMO instance).
+               e3t_ini   (:,:,jk) = e3t(:,:,jk,Kmm)                      * tmask(:,:,jk)  ! initial vertical scale factors
+               tmask_ini (:,:,jk) = tmask(:,:,jk)                                       ! initial mask
+               hc_loc_ini(:,:,jk) = ts(:,:,jk,jp_tem,Kmm) * e3t(:,:,jk,Kmm) * tmask(:,:,jk)  ! initial heat content
+               sc_loc_ini(:,:,jk) = ts(:,:,jk,jp_sal,Kmm) * e3t(:,:,jk,Kmm) * tmask(:,:,jk)  ! initial salt content
+            END DO
+            frc_v = 0._wp                                           ! volume       trend due to forcing
+            frc_t = 0._wp                                           ! heat content   -    -   -    -   
+            frc_s = 0._wp                                           ! salt content   -    -   -    -        
+            IF( ln_linssh ) THEN
+               IF( ln_isfcav ) THEN
+                  DO ji = 1, jpi
+                     DO jj = 1, jpj
+                        ssh_hc_loc_ini(ji,jj) = ts(ji,jj,mikt(ji,jj),jp_tem,Kmm) * ssh(ji,jj,Kmm)   ! initial heat content in ssh
+                        ssh_sc_loc_ini(ji,jj) = ts(ji,jj,mikt(ji,jj),jp_sal,Kmm) * ssh(ji,jj,Kmm)   ! initial salt content in ssh
+                     END DO
+                   END DO
+                ELSE
+                  ssh_hc_loc_ini(:,:) = ts(:,:,1,jp_tem,Kmm) * ssh(:,:,Kmm)   ! initial heat content in ssh
+                  ssh_sc_loc_ini(:,:) = ts(:,:,1,jp_sal,Kmm) * ssh(:,:,Kmm)   ! initial salt content in ssh
+               END IF
+               frc_wn_t = 0._wp                                       ! initial heat content misfit due to free surface
+               frc_wn_s = 0._wp                                       ! initial salt content misfit due to free surface
+            ENDIF
+         ENDIF
+         !
+      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
+         !                                   ! -------------------
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) '   dia_hsb_rst : write restart at it= ', kt,' date= ', ndastp
+         IF(lwp) WRITE(numout,*)
+         !
+         IF( lwxios ) CALL iom_swap(      cwxios_context          )
+         CALL iom_rstput( kt, nitrst, numrow, 'frc_v', frc_v, ldxios = lwxios )
+         CALL iom_rstput( kt, nitrst, numrow, 'frc_t', frc_t, ldxios = lwxios )
+         CALL iom_rstput( kt, nitrst, numrow, 'frc_s', frc_s, ldxios = lwxios )
+         IF( ln_linssh ) THEN
+            CALL iom_rstput( kt, nitrst, numrow, 'frc_wn_t', frc_wn_t, ldxios = lwxios )
+            CALL iom_rstput( kt, nitrst, numrow, 'frc_wn_s', frc_wn_s, ldxios = lwxios )
+         ENDIF
+         CALL iom_rstput( kt, nitrst, numrow, 'surf_ini'  , surf_ini  , ldxios = lwxios )      ! ice sheet coupling
+         CALL iom_rstput( kt, nitrst, numrow, 'ssh_ini'   , ssh_ini   , ldxios = lwxios )
+         CALL iom_rstput( kt, nitrst, numrow, 'e3t_ini'   , e3t_ini   , ldxios = lwxios )
+         CALL iom_rstput( kt, nitrst, numrow, 'tmask_ini' , tmask_ini , ldxios = lwxios )
+         CALL iom_rstput( kt, nitrst, numrow, 'hc_loc_ini', hc_loc_ini, ldxios = lwxios )
+         CALL iom_rstput( kt, nitrst, numrow, 'sc_loc_ini', sc_loc_ini, ldxios = lwxios )
+         IF( ln_linssh ) THEN
+            CALL iom_rstput( kt, nitrst, numrow, 'ssh_hc_loc_ini', ssh_hc_loc_ini, ldxios = lwxios )
+            CALL iom_rstput( kt, nitrst, numrow, 'ssh_sc_loc_ini', ssh_sc_loc_ini, ldxios = lwxios )
+         ENDIF
+         IF( lwxios ) CALL iom_swap(      cxios_context          )
+         !
+      ENDIF
+      !
+   END SUBROUTINE dia_hsb_rst
+
+
+   SUBROUTINE dia_hsb_init( Kmm )
+      !!---------------------------------------------------------------------------
+      !!                  ***  ROUTINE dia_hsb  ***
+      !!     
+      !! ** Purpose: Initialization for the heat salt volume budgets
+      !!	
+      !! ** Method : Compute initial heat content, salt content and volume
+      !!
+      !! ** Action : - Compute initial heat content, salt content and volume
+      !!             - Initialize forcing trends
+      !!             - Compute coefficients for conversion
+      !!---------------------------------------------------------------------------
+      INTEGER, INTENT(in) :: Kmm ! time level index
+      !
+      INTEGER ::   ierror, ios   ! local integer
+      !!
+      NAMELIST/namhsb/ ln_diahsb
+      !!----------------------------------------------------------------------
+      !
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'dia_hsb_init : heat and salt budgets diagnostics'
+         WRITE(numout,*) '~~~~~~~~~~~~ '
+      ENDIF
+      READ  ( numnam_ref, namhsb, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namhsb in reference namelist' )
+      READ  ( numnam_cfg, namhsb, IOSTAT = ios, ERR = 902 )
+902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namhsb in configuration namelist' )
+      IF(lwm) WRITE( numond, namhsb )
+
+      IF(lwp) THEN
+         WRITE(numout,*) '   Namelist  namhsb :'
+         WRITE(numout,*) '      check the heat and salt budgets (T) or not (F)       ln_diahsb = ', ln_diahsb
+      ENDIF
+      !
+      IF( .NOT. ln_diahsb )   RETURN
+
+      IF(lwxios) THEN
+! define variables in restart file when writing with XIOS
+        CALL iom_set_rstw_var_active('frc_v')
+        CALL iom_set_rstw_var_active('frc_t')
+        CALL iom_set_rstw_var_active('frc_s')
+        CALL iom_set_rstw_var_active('surf_ini')
+        CALL iom_set_rstw_var_active('ssh_ini')
+        CALL iom_set_rstw_var_active('e3t_ini')
+        CALL iom_set_rstw_var_active('hc_loc_ini')
+        CALL iom_set_rstw_var_active('sc_loc_ini')
+        IF( ln_linssh ) THEN
+           CALL iom_set_rstw_var_active('ssh_hc_loc_ini')
+           CALL iom_set_rstw_var_active('ssh_sc_loc_ini')
+           CALL iom_set_rstw_var_active('frc_wn_t')
+           CALL iom_set_rstw_var_active('frc_wn_s')
+        ENDIF
+      ENDIF
+      ! ------------------- !
+      ! 1 - Allocate memory !
+      ! ------------------- !
+      ALLOCATE( hc_loc_ini(jpi,jpj,jpk), sc_loc_ini(jpi,jpj,jpk), surf_ini(jpi,jpj), &
+         &      e3t_ini(jpi,jpj,jpk), surf(jpi,jpj),  ssh_ini(jpi,jpj), tmask_ini(jpi,jpj,jpk),STAT=ierror  )
+      IF( ierror > 0 ) THEN
+         CALL ctl_stop( 'dia_hsb_init: unable to allocate hc_loc_ini' )   ;   RETURN
+      ENDIF
+
+      IF( ln_linssh )   ALLOCATE( ssh_hc_loc_ini(jpi,jpj), ssh_sc_loc_ini(jpi,jpj),STAT=ierror )
+      IF( ierror > 0 ) THEN
+         CALL ctl_stop( 'dia_hsb: unable to allocate ssh_hc_loc_ini' )   ;   RETURN
+      ENDIF
+
+      ! ----------------------------------------------- !
+      ! 2 - Time independant variables and file opening !
+      ! ----------------------------------------------- !
+      surf(:,:) = e1e2t(:,:) * tmask_i(:,:)               ! masked surface grid cell area
+      surf_tot  = glob_sum( 'diahsb', surf(:,:) )         ! total ocean surface area
+
+      IF( ln_bdy ) CALL ctl_warn( 'dia_hsb_init: heat/salt budget does not consider open boundary fluxes' )         
+      !
+      ! ---------------------------------- !
+      ! 4 - initial conservation variables !
+      ! ---------------------------------- !
+      CALL dia_hsb_rst( nit000, Kmm, 'READ' )  !* read or initialize all required files
+      !
+   END SUBROUTINE dia_hsb_init
+
+   !!======================================================================
+END MODULE diahsb
+#endif
Index: NEMO/branches/2020/dev_r13747_HPC-10_HPDAonline_DiagGPU/tests/DIA_GPU/MY_SRC/diahsb_gpu.h90
===================================================================
--- NEMO/branches/2020/dev_r13747_HPC-10_HPDAonline_DiagGPU/tests/DIA_GPU/MY_SRC/diahsb_gpu.h90	(revision 14091)
+++ NEMO/branches/2020/dev_r13747_HPC-10_HPDAonline_DiagGPU/tests/DIA_GPU/MY_SRC/diahsb_gpu.h90	(revision 14091)
@@ -0,0 +1,728 @@
+MODULE diahsb
+   !!======================================================================
+   !! *** MODULE diahsb ***
+   !! Ocean diagnostics: Heat, salt and volume budgets
+   !!======================================================================
+   !! History : 3.3 ! 2010-09 (M. Leclair) Original code
+   !! ! 2012-10 (C. Rousset) add iom_put
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !! dia_hsb : Diagnose the conservation of ocean heat and salt contents, and volume
+   !! dia_hsb_rst : Read or write DIA file in restart file
+   !! dia_hsb_init : Initialization of the conservation diagnostic
+   !!----------------------------------------------------------------------
+   USE oce ! ocean dynamics and tracers
+   USE dom_oce ! ocean space and time domain
+   USE phycst ! physical constants
+   USE sbc_oce ! surface thermohaline fluxes
+   USE isf_oce ! ice shelf fluxes
+   USE sbcrnf ! river runoff
+   USE domvvl ! vertical scale factors
+   USE traqsr ! penetrative solar radiation
+   USE trabbc ! bottom boundary condition
+   USE trabbc ! bottom boundary condition
+   USE restart ! ocean restart
+   USE bdy_oce , ONLY : ln_bdy
+   !
+   USE iom ! I/O manager
+   USE in_out_manager ! I/O manager
+
+    USE gpu_manager ! GPU manager
+    USE cudafor ! CUDA toolkit libs
+    USE cuda_fortran ! CUDA routines
+    !USE nvtx ! CUDA profiling/DEGUG tools
+
+   USE lib_fortran ! glob_sum
+   USE lib_mpp ! distributed memory computing library
+   USE timing ! preformance summary
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC dia_hsb ! routine called by step.F90
+   PUBLIC dia_hsb_init ! routine called by nemogcm.F90
+
+   LOGICAL, PUBLIC :: ln_diahsb !: check the heat and salt budgets
+
+   REAL(wp) :: surf_tot ! ocean surface
+
+
+
+
+   REAL(wp) , DIMENSION(2), SAVE :: frc_t, frc_s, frc_v ! global forcing trends
+   REAL(wp) , DIMENSION(2), SAVE :: frc_wn_t, frc_wn_s ! global forcing trends
+
+   REAL(wp), DIMENSION(:,:) , ALLOCATABLE :: surf
+   REAL(wp), DIMENSION(:,:) , ALLOCATABLE :: surf_ini , ssh_ini !
+   REAL(wp), DIMENSION(:,:) , ALLOCATABLE :: ssh_hc_loc_ini, ssh_sc_loc_ini !
+
+
+
+   REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, PINNED :: hc_loc_ini, sc_loc_ini !
+   REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: e3t_ini !
+   REAL(wp), DIMENSION(:) , ALLOCATABLE, PINNED, SAVE :: h_ztmpv, h_ztmph, h_ztmps, h_ztmp !
+
+   REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: tmask_ini
+
+
+    !Device data associate to PUBLIC arrays
+    REAL(8), DIMENSION(:,:,:,:) , ALLOCATABLE, DEVICE :: d_e3t !
+    REAL(8), DIMENSION(:,:,:) , ALLOCATABLE, DEVICE :: d_tmask !
+    REAL(8), DIMENSION(:,:) , ALLOCATABLE, DEVICE :: d_tmask_h !
+    REAL(8), DIMENSION(:,:,:) , ALLOCATABLE, DEVICE :: d_tmask_ini !
+    REAL(8), DIMENSION(:,:,:,:,:), ALLOCATABLE, DEVICE :: d_ts !
+    !Device data associate to LOCAL/DEVICE arrays
+    REAL(8), DEVICE , DIMENSION(:,:) , ALLOCATABLE :: d_surf !
+    REAL(8), DEVICE , DIMENSION(:,:) , ALLOCATABLE :: d_surf_ini !
+    REAL(8), DEVICE , DIMENSION(:,:,:) , ALLOCATABLE :: d_hc_loc_ini !
+    REAL(8), DEVICE , DIMENSION(:,:,:) , ALLOCATABLE :: d_sc_loc_ini !
+    REAL(8), DEVICE , DIMENSION(:,:,:) , ALLOCATABLE :: d_e3t_ini !
+    REAL(8), DEVICE , DIMENSION(:,:,:) , ALLOCATABLE :: d_zwrkv, d_zwrkh, d_zwrks, d_zwrk ! 3D GPU workspace
+    REAL(8), DEVICE :: ztmpv, ztmph, ztmps, ztmp ! Device Reduction
+    !
+    INTEGER :: globsize ! 3D workspace size
+    type(dim3) :: dimGrid, dimBlock ! cuda parameters
+    INTEGER, parameter :: nstreams = 3 ! Streams Number
+    INTEGER(kind=cuda_stream_kind) :: stream(nstreams), str ! Stream ID
+    !DEBUG
+    !REAL(8) , save , DIMENSION(:,:,:) , ALLOCATABLE :: prev_3d
+    !REAL(8) :: accum
+
+
+
+
+   !! * Substitutions
+!# include "domzgr_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+   !! $Id$
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE dia_hsb( kt, Kbb, Kmm )
+      !!---------------------------------------------------------------------------
+      !! *** ROUTINE dia_hsb ***
+      !!
+      !! ** Purpose: Compute the ocean global heat content, salt content and volume conservation
+      !!
+      !! ** Method : - Compute the deviation of heat content, salt content and volume
+      !! at the current time step from their values at nit000
+      !! - Compute the contribution of forcing and remove it from these deviations
+      !!
+      !!---------------------------------------------------------------------------
+      INTEGER, INTENT(in) :: kt ! ocean time-step index
+      INTEGER, INTENT(in) :: Kbb, Kmm ! ocean time level indices
+      !
+
+      INTEGER, VALUE :: ji, jj, jk, kts ! dummy loop indice
+      INTEGER, VALUE :: localsize ! jpi * jpj * jpk
+      INTEGER :: istat ! CUDA error check
+      COMPLEX :: ctmp ! dummy complex number
+      INTEGER(kind=cuda_stream_kind) :: str ! dummy kernel stream
+      INTEGER :: tile_n, tile_b ! tile indexe. _n now, _b before
+      REAL(wp) , DIMENSION(2), SAVE :: zdiff_hc1, zdiff_sc1 ! heat and salt content variations
+      REAL(wp) , DIMENSION(2), SAVE :: zdiff_hc, zdiff_sc ! - - - -
+      REAL(wp) , DIMENSION(2), SAVE :: zdiff_v2 ! volume variation
+      REAL(wp) , DIMENSION(2), SAVE :: zdiff_v1 ! volume variation
+      REAL(wp) , DIMENSION(2), SAVE :: zerr_hc1, zerr_sc1 ! heat and salt content misfit
+      REAL(wp) , DIMENSION(2), SAVE :: zvol_tot ! volume
+      REAL(wp) , DIMENSION(2), SAVE :: z_frc_trd_t, z_frc_trd_s ! - -
+      REAL(wp) , DIMENSION(2), SAVE :: z_frc_trd_v ! - -
+      REAL(wp) , DIMENSION(2), SAVE :: z_wn_trd_t, z_wn_trd_s ! - -
+      REAL(wp) , DIMENSION(2), SAVE :: z_ssh_hc, z_ssh_sc ! - -
+# 147 "diahsb_new.F90"
+      REAL(wp), DIMENSION(jpi,jpj) :: z2d0, z2d1 ! 2D workspace
+      REAL(wp), DIMENSION(jpi,jpj,jpkm1) :: zwrk ! 3D workspace
+      !!---------------------------------------------------------------------------
+      IF( ln_timing ) CALL timing_start('dia_hsb')
+
+      localsize = jpi * jpj * jpk
+      kts = kt
+      IF (kts == 1) THEN
+          tile_n = 1
+          tile_b = 1
+      ELSE
+          IF( MOD(kts,2) == 0) THEN
+              tile_n = 2
+              tile_b = 1
+          ELSE IF( MOD(kts,2) == 1 ) THEN
+              tile_n = 1
+              tile_b = 2
+          END IF
+      END IF
+
+      !
+      ts(:,:,:,1,Kmm) = ts(:,:,:,1,Kmm) * tmask(:,:,:) ; ts(:,:,:,1,Kbb) = ts(:,:,:,1,Kbb) * tmask(:,:,:) ;
+      ts(:,:,:,2,Kmm) = ts(:,:,:,2,Kmm) * tmask(:,:,:) ; ts(:,:,:,2,Kbb) = ts(:,:,:,2,Kbb) * tmask(:,:,:) ;
+      ! ------------------------- !
+      ! 1 - Trends due to forcing !
+      ! ------------------------- !
+
+      z_frc_trd_v(tile_n) = r1_rho0 * glob_sum( 'diahsb', - ( emp(:,:) - rnf(:,:) + fwfisf_cav(:,:) + fwfisf_par(:,:) ) * surf(:,:) )! volume fluxes
+      z_frc_trd_t(tile_n) = glob_sum( 'diahsb', sbc_tsc(:,:,jp_tem) * surf(:,:) ) ! heat fluxes
+      z_frc_trd_s(tile_n) = glob_sum( 'diahsb', sbc_tsc(:,:,jp_sal) * surf(:,:) ) ! salt fluxes
+      ! ! Add runoff heat & salt input
+      IF( ln_rnf ) z_frc_trd_t(tile_n) = z_frc_trd_t(tile_n) + glob_sum( 'diahsb', rnf_tsc(:,:,jp_tem) * surf(:,:) )
+      IF( ln_rnf_sal) z_frc_trd_s(tile_n) = z_frc_trd_s(tile_n) + glob_sum( 'diahsb', rnf_tsc(:,:,jp_sal) * surf(:,:) ) ! Add ice shelf heat & salt input
+      IF( ln_isf ) z_frc_trd_t(tile_n) = z_frc_trd_t(tile_n) &
+                             & + glob_sum( 'diahsb', ( risf_cav_tsc(:,:,jp_tem) + risf_par_tsc(:,:,jp_tem) ) * surf(:,:) ) ! Add penetrative solar radiation
+      IF( ln_traqsr ) z_frc_trd_t(tile_n) = z_frc_trd_t(tile_n) + r1_rho0_rcp * glob_sum( 'diahsb', qsr (:,:) * surf(:,:) ) ! Add geothermal heat flux
+      IF( ln_trabbc ) z_frc_trd_t(tile_n) = z_frc_trd_t(tile_n) + glob_sum( 'diahsb', qgh_trd0(:,:) * surf(:,:) )
+      !
+      IF( ln_linssh ) THEN
+         IF( ln_isfcav ) THEN
+            DO ji=1,jpi
+               DO jj=1,jpj
+                  z2d0(ji,jj) = surf(ji,jj) * ww(ji,jj,mikt(ji,jj)) * ts(ji,jj,mikt(ji,jj),jp_tem,Kbb)
+                  z2d1(ji,jj) = surf(ji,jj) * ww(ji,jj,mikt(ji,jj)) * ts(ji,jj,mikt(ji,jj),jp_sal,Kbb)
+               END DO
+            END DO
+         ELSE
+            z2d0(:,:) = surf(:,:) * ww(:,:,1) * ts(:,:,1,jp_tem,Kbb)
+            z2d1(:,:) = surf(:,:) * ww(:,:,1) * ts(:,:,1,jp_sal,Kbb)
+         END IF
+
+         z_wn_trd_t(tile_n) = - glob_sum( 'diahsb', z2d0 )
+         z_wn_trd_s(tile_n) = - glob_sum( 'diahsb', z2d1 )
+
+
+
+
+      ENDIF
+
+      IF (kts>1) THEN
+         frc_v(tile_n) = frc_v(tile_b)
+         frc_t(tile_n) = frc_t(tile_b)
+         frc_s(tile_n) = frc_s(tile_b)
+         frc_wn_t(tile_n) = frc_wn_t(tile_b)
+         frc_wn_s(tile_n) = frc_wn_s(tile_b)
+      END IF
+      frc_v(tile_n) = frc_v(tile_n) + z_frc_trd_v(tile_n) * rn_Dt
+      frc_t(tile_n) = frc_t(tile_n) + z_frc_trd_t(tile_n) * rn_Dt
+      frc_s(tile_n) = frc_s(tile_n) + z_frc_trd_s(tile_n) * rn_Dt
+      ! ! Advection flux through fixed surface (z=0)
+      IF( ln_linssh ) THEN
+          frc_wn_t(tile_n) = frc_wn_t(tile_n) + z_wn_trd_t(tile_n) * rn_Dt
+          frc_wn_s(tile_n) = frc_wn_s(tile_n) + z_wn_trd_s(tile_n) * rn_Dt
+      ENDIF
+      ! ------------------------ !
+      ! 2 - Content variations !
+      ! ------------------------ !
+      ! glob_sum_full is needed because you keep the full interior domain to compute the sum (iscpl)
+
+      ! ! volume variation (calculated with ssh)
+
+      zdiff_v1(tile_n) = glob_sum_full( 'diahsb', surf(:,:)*ssh(:,:,Kmm) - surf_ini(:,:)*ssh_ini(:,:) )
+
+
+
+      ! ! heat & salt content variation (associated with ssh)
+      IF( ln_linssh ) THEN ! linear free surface case
+         IF( ln_isfcav ) THEN ! ISF case
+            DO ji = 1, jpi
+               DO jj = 1, jpj
+                  z2d0(ji,jj) = surf(ji,jj) * ( ts(ji,jj,mikt(ji,jj),jp_tem,Kmm) * ssh(ji,jj,Kmm) - ssh_hc_loc_ini(ji,jj) )
+                  z2d1(ji,jj) = surf(ji,jj) * ( ts(ji,jj,mikt(ji,jj),jp_sal,Kmm) * ssh(ji,jj,Kmm) - ssh_sc_loc_ini(ji,jj) )
+               END DO
+            END DO
+         ELSE ! no under ice-shelf seas
+            z2d0(:,:) = surf(:,:) * ( ts(:,:,1,jp_tem,Kmm) * ssh(:,:,Kmm) - ssh_hc_loc_ini(:,:) )
+            z2d1(:,:) = surf(:,:) * ( ts(:,:,1,jp_sal,Kmm) * ssh(:,:,Kmm) - ssh_sc_loc_ini(:,:) )
+         END IF
+
+         z_ssh_hc(tile_n) = glob_sum_full( 'diahsb', z2d0 )
+         z_ssh_sc(tile_n) = glob_sum_full( 'diahsb', z2d1 )
+
+
+
+
+      ENDIF
+
+      str = stream(tile_n)
+      istat = 0
+      istat = cudaMemcpyAsync( d_e3t, e3t , jpi*jpj*jpk*jpt , str ) + istat
+      istat = cudaMemcpyAsync(d_hc_loc_ini, hc_loc_ini, jpi*jpj*jpk , str ) + istat
+      istat = cudaMemcpyAsync(d_sc_loc_ini, sc_loc_ini, jpi*jpj*jpk , str ) + istat
+      istat = cudaMemcpyAsync( d_ts, ts , jpi*jpj*jpk*2*jpt, str ) + istat
+      IF( istat /= 0 ) THEN
+         CALL ctl_stop( 'dia_hsb: unable to async GPU copy H2D' ) ; RETURN
+      ENDIF
+      dimBlock = dim3(4,4,4)
+      dimGrid = dim3( ceiling( real( jpi ) / dimBlock%x ) , ceiling( real( jpj ) / dimBlock%y ) , &
+                                                             ceiling( real( jpkm1 ) / dimBlock%z ) )
+      !
+      CALL dia_hsb_kernel<<<dimGrid, dimBlock, 0, str>>> (d_surf , d_e3t, d_surf_ini, d_e3t_ini, &
+      & d_ts, d_hc_loc_ini, d_sc_loc_ini, d_tmask, d_tmask_ini, d_zwrkv, d_zwrkh, d_zwrks, d_zwrk, jpi, jpj, jpk, jpt, Kmm)
+
+
+      CALL filter_cuda<<<dimGrid, dimBlock, 0, str>>>(d_zwrkv , d_tmask_h , jpi, jpj, jpk)
+      CALL filter_cuda<<<dimGrid, dimBlock, 0, str>>>(d_zwrkh , d_tmask_h , jpi, jpj, jpk)
+      CALL filter_cuda<<<dimGrid, dimBlock, 0, str>>>(d_zwrks , d_tmask_h , jpi, jpj, jpk)
+      CALL filter_cuda<<<dimGrid, dimBlock, 0, str>>>(d_zwrk , d_tmask_h , jpi, jpj, jpk)
+
+      ztmpv = 0.e0
+      ztmph = 0.e0
+      ztmps = 0.e0
+      ztmp = 0.e0
+
+      istat = 0
+      !$cuf kernel do <<< *, *, stream=str >>>
+      do ji = 1, localsize
+        ztmpv = ztmpv + d_zwrkv(ji)
+      end do
+      istat = cudaMemcpyAsync( h_ztmpv(tile_n) , ztmpv , 1 , str ) + istat
+      !$cuf kernel do <<< *, *, stream=str >>>
+      do ji = 1, localsize
+        ztmph = ztmph + d_zwrkh(ji)
+      end do
+      istat = cudaMemcpyAsync( h_ztmph(tile_n) , ztmph , 1 , str ) + istat
+      !$cuf kernel do <<< *, *, stream=str >>>
+      do ji = 1, localsize
+        ztmps = ztmps + d_zwrks(ji)
+      end do
+      istat = cudaMemcpyAsync( h_ztmps(tile_n) , ztmps , 1 , str ) + istat
+      !$cuf kernel do <<< *, *, stream=str >>>
+      do ji = 1, localsize
+        ztmp = ztmp + d_zwrk(ji)
+      end do
+      istat = cudaMemcpyAsync( h_ztmp (tile_n) , ztmp , 1 , str ) + istat
+      !
+      IF( istat /= 0 ) THEN
+         CALL ctl_stop( 'dia_hsb: unable to async GPU copy D2H' ) ; RETURN
+      ENDIF
+      !
+      istat = cudaStreamSynchronize(stream(tile_b))
+      !
+      IF( istat /= 0 ) THEN
+         CALL ctl_stop( 'dia_hsb: unable to stream synchronize' ) ; RETURN
+      ENDIF
+      !
+      ctmp = CMPLX( h_ztmpv(tile_b) , 0.e0, 8 )
+      CALL mpp_sum('diahsb', ctmp )
+      zdiff_v2(tile_b) = REAL( ctmp, 8 )
+
+      ctmp = CMPLX( h_ztmph(tile_b) , 0.e0, 8 )
+      CALL mpp_sum('diahsb', ctmp )
+      zdiff_hc(tile_b) = REAL( ctmp, 8 )
+
+      ctmp = CMPLX( h_ztmps(tile_b) , 0.e0, 8 )
+      CALL mpp_sum('diahsb', ctmp )
+      zdiff_sc(tile_b) = REAL( ctmp, 8 )
+
+      ctmp = CMPLX( h_ztmp(tile_b) , 0.e0, 8 )
+      CALL mpp_sum('diahsb', ctmp )
+      zvol_tot(tile_b) = REAL( ctmp, 8 )
+
+      IF ( kt == nitend ) THEN
+         !
+         istat = cudaStreamSynchronize(stream(tile_n))
+         IF( istat /= 0 ) THEN
+         CALL ctl_stop( 'dia_hsb: unable to stream synchronize' ) ; RETURN
+         ENDIF
+         !
+         ctmp = CMPLX( h_ztmpv(tile_n) , 0.e0, 8 )
+         CALL mpp_sum('diahsb', ctmp )
+         zdiff_v2(tile_n) = REAL( ctmp, 8 )
+
+         ctmp = CMPLX( h_ztmph(tile_n) , 0.e0, 8 )
+         CALL mpp_sum('diahsb', ctmp )
+         zdiff_hc(tile_n) = REAL( ctmp, 8 )
+
+         ctmp = CMPLX( h_ztmps(tile_n) , 0.e0, 8 )
+         CALL mpp_sum('diahsb', ctmp )
+         zdiff_sc(tile_n) = REAL( ctmp, 8 )
+
+         ctmp = CMPLX( h_ztmp(tile_n) , 0.e0, 8 )
+         CALL mpp_sum('diahsb', ctmp )
+         zvol_tot(tile_n) = REAL( ctmp, 8 )
+      ENDIF
+      ! ------------------------ !
+      ! 3 - Drifts !
+      ! ------------------------ !
+
+      IF ( kt > 1 ) THEN
+            kts = kts - 1
+            zdiff_v1(tile_b) = zdiff_v1(tile_b) - frc_v(tile_b)
+            IF( .NOT.ln_linssh ) zdiff_v2(tile_b) = zdiff_v2(tile_b) - frc_v(tile_b)
+            zdiff_hc(tile_b) = zdiff_hc(tile_b) - frc_t(tile_b)
+            zdiff_sc(tile_b) = zdiff_sc(tile_b) - frc_s(tile_b)
+            IF( ln_linssh ) THEN
+                zdiff_hc1(tile_b) = zdiff_hc (tile_b) + z_ssh_hc(tile_b)
+                zdiff_sc1(tile_b) = zdiff_sc (tile_b) + z_ssh_sc(tile_b)
+                zerr_hc1 (tile_b) = z_ssh_hc(tile_b) - frc_wn_t(tile_b)
+                zerr_sc1 (tile_b) = z_ssh_sc(tile_b) - frc_wn_s(tile_b)
+            ENDIF
+            !!gm to be added ?
+            ! IF( ln_linssh ) THEN ! fixed volume, add the ssh contribution
+            ! zvol_tot = zvol_tot + glob_sum( 'diahsb', surf(:,:) * sshn(:,:) )
+            ! ENDIF
+            !!gm end
+
+            CALL iom_put( 'bgfrcvol' , frc_v(tile_b) * 1.e-9 ) ! vol - surface forcing (km3)
+            CALL iom_put( 'bgfrctem' , frc_t(tile_b) * rho0 * rcp * 1.e-20 ) ! hc - surface forcing (1.e20 J)
+            CALL iom_put( 'bgfrchfx' , frc_t(tile_b) * rho0 * rcp / & ! hc - surface forcing (W/m2)
+                & ( surf_tot * kts * rn_Dt ) )
+            CALL iom_put( 'bgfrcsal' , frc_s(tile_b) * 1.e-9 ) ! sc - surface forcing (psu*km3)
+            IF( .NOT. ln_linssh ) THEN
+                CALL iom_put( 'bgtemper' , zdiff_hc(tile_b) / zvol_tot(tile_b) ) ! Temperature drift (C)
+                CALL iom_put( 'bgsaline' , zdiff_sc(tile_b) / zvol_tot(tile_b) ) ! Salinity drift (PSU)
+                CALL iom_put( 'bgheatco' , zdiff_hc(tile_b) * 1.e-20 * rho0 * rcp ) ! Heat content drift (1.e20 J)
+                CALL iom_put( 'bgheatfx' , zdiff_hc(tile_b) * rho0 * rcp / & ! Heat flux drift (W/m2)
+                    & ( surf_tot * kts * rn_Dt ) )
+                CALL iom_put( 'bgsaltco' , zdiff_sc(tile_b) * 1.e-9 ) ! Salt content drift (psu*km3)
+                CALL iom_put( 'bgvolssh' , zdiff_v1(tile_b) * 1.e-9 ) ! volume ssh drift (km3)
+                CALL iom_put( 'bgvole3t' , zdiff_v2(tile_b) * 1.e-9 ) ! volume e3t drift (km3)
+                !
+! IF( lwp ) THEN
+! WRITE(numout,*)
+! WRITE(numout,*) 'dia_hsb : last time step hsb diagnostics: at it= ', kt,' date= ', ndastp
+! WRITE(numout,*) '~~~~~~~'
+! WRITE(numout,*) '   Temperature drift = ', zdiff_hc(tile_b) / zvol_tot(tile_b), ' C'
+! WRITE(numout,*) '   Salinity12  drift = ', zdiff_sc(tile_b) / zvol_tot(tile_b), ' PSU'
+! WRITE(numout,*) '   volume ssh  drift = ', zdiff_v1(tile_b) * 1.e-9 , ' km^3'
+! WRITE(numout,*) '   volume e3t  drift = ', zdiff_v2(tile_b) * 1.e-9 , ' km^3'
+! ENDIF
+            ELSE
+                CALL iom_put( 'bgtemper' , zdiff_hc1(tile_b) / zvol_tot(tile_b)) ! Heat content drift (C)
+                CALL iom_put( 'bgsaline' , zdiff_sc1(tile_b) / zvol_tot(tile_b)) ! Salt content drift (PSU)
+                CALL iom_put( 'bgheatco' , zdiff_hc1(tile_b) * 1.e-20 * rho0 * rcp ) ! Heat content drift (1.e20 J)
+                CALL iom_put( 'bgheatfx' , zdiff_hc1(tile_b) * rho0 * rcp / & ! Heat flux drift (W/m2)
+                    & ( surf_tot * kts * rn_Dt ) )
+                CALL iom_put( 'bgsaltco' , zdiff_sc1(tile_b) * 1.e-9 ) ! Salt content drift (psu*km3)
+                CALL iom_put( 'bgvolssh' , zdiff_v1(tile_b) * 1.e-9 ) ! volume ssh drift (km3)
+                CALL iom_put( 'bgmistem' , zerr_hc1(tile_b) / zvol_tot(tile_b) ) ! hc - error due to free surface (C)
+                CALL iom_put( 'bgmissal' , zerr_sc1(tile_b) / zvol_tot(tile_b) ) ! sc - error due to free surface (psu)
+            ENDIF
+            !
+            IF( lrst_oce ) CALL dia_hsb_rst( kts, Kmm, tile_n, 'WRITE' )
+        !
+        END IF
+        IF ( kt == nitend ) THEN
+
+            zdiff_v1(tile_n) = zdiff_v1(tile_n) - frc_v(tile_n)
+            IF( .NOT.ln_linssh ) zdiff_v2(tile_n) = zdiff_v2(tile_n) - frc_v(tile_n)
+            zdiff_hc(tile_n) = zdiff_hc(tile_n) - frc_t(tile_n)
+            zdiff_sc(tile_n) = zdiff_sc(tile_n) - frc_s(tile_n)
+            IF( ln_linssh ) THEN
+                zdiff_hc1(tile_n) = zdiff_hc (tile_n) + z_ssh_hc(tile_n)
+                zdiff_sc1(tile_n) = zdiff_sc (tile_n) + z_ssh_sc(tile_n)
+                zerr_hc1 (tile_n) = z_ssh_hc(tile_n) - frc_wn_t(tile_n)
+                zerr_sc1 (tile_n) = z_ssh_sc(tile_n) - frc_wn_s(tile_n)
+            ENDIF
+            !!gm to be added ?
+            ! IF( ln_linssh ) THEN ! fixed volume, add the ssh contribution
+            ! zvol_tot = zvol_tot + glob_sum( 'diahsb', surf(:,:) * sshn(:,:) )
+            ! ENDIF
+            !!gm end
+
+            CALL iom_put( 'bgfrcvol' , frc_v(tile_n) * 1.e-9 ) ! vol - surface forcing (km3)
+            CALL iom_put( 'bgfrctem' , frc_t(tile_n) * rho0 * rcp * 1.e-20 ) ! hc - surface forcing (1.e20 J)
+            CALL iom_put( 'bgfrchfx' , frc_t(tile_n) * rho0 * rcp / & ! hc - surface forcing (W/m2)
+                & ( surf_tot * kt * rn_Dt ) )
+            CALL iom_put( 'bgfrcsal' , frc_s(tile_n) * 1.e-9 ) ! sc - surface forcing (psu*km3)
+            IF( .NOT. ln_linssh ) THEN
+                CALL iom_put( 'bgtemper' , zdiff_hc(tile_n) / zvol_tot(tile_n) ) ! Temperature drift (C)
+                CALL iom_put( 'bgsaline' , zdiff_sc(tile_n) / zvol_tot(tile_n) ) ! Salinity drift (PSU)
+                CALL iom_put( 'bgheatco' , zdiff_hc(tile_n) * 1.e-20 * rho0 * rcp ) ! Heat content drift (1.e20 J)
+                CALL iom_put( 'bgheatfx' , zdiff_hc(tile_n) * rho0 * rcp / & ! Heat flux drift (W/m2)
+                    & ( surf_tot * kt * rn_Dt ) )
+                CALL iom_put( 'bgsaltco' , zdiff_sc(tile_n) * 1.e-9 ) ! Salt content drift (psu*km3)
+                CALL iom_put( 'bgvolssh' , zdiff_v1(tile_n) * 1.e-9 ) ! volume ssh drift (km3)
+                CALL iom_put( 'bgvole3t' , zdiff_v2(tile_n) * 1.e-9 ) ! volume e3t drift (km3)
+                !
+                IF( kt == nitend .AND. lwp ) THEN
+                    WRITE(numout,*)
+                    WRITE(numout,*) 'dia_hsb : last time step hsb diagnostics: at it= ', kt,' date= ', ndastp
+                    WRITE(numout,*) '~~~~~~~'
+                    WRITE(numout,*) '   Temperature drift = ', zdiff_hc(tile_n) / zvol_tot(tile_n), ' C'
+                    WRITE(numout,*) '   Salinity  drift   = ', zdiff_sc(tile_n) / zvol_tot(tile_n), ' PSU'
+                    WRITE(numout,*) '   volume ssh  drift = ', zdiff_v1(tile_n) * 1.e-9 , ' km^3'
+                    WRITE(numout,*) '   volume e3t  drift = ', zdiff_v2(tile_n) * 1.e-9 , ' km^3'
+                !
+                ENDIF
+               !
+            ELSE
+                CALL iom_put( 'bgtemper' , zdiff_hc1(tile_n) / zvol_tot(tile_n)) ! Heat content drift (C)
+                CALL iom_put( 'bgsaline' , zdiff_sc1(tile_n) / zvol_tot(tile_n)) ! Salt content drift (PSU)
+                CALL iom_put( 'bgheatco' , zdiff_hc1(tile_n) * 1.e-20 * rho0 * rcp ) ! Heat content drift (1.e20 J)
+                CALL iom_put( 'bgheatfx' , zdiff_hc1(tile_n) * rho0 * rcp / & ! Heat flux drift (W/m2)
+                    & ( surf_tot * kt * rn_Dt ) )
+                CALL iom_put( 'bgsaltco' , zdiff_sc1(tile_n) * 1.e-9 ) ! Salt content drift (psu*km3)
+                CALL iom_put( 'bgvolssh' , zdiff_v1(tile_n) * 1.e-9 ) ! volume ssh drift (km3)
+                CALL iom_put( 'bgmistem' , zerr_hc1(tile_n) / zvol_tot(tile_n) ) ! hc - error due to free surface (C)
+                CALL iom_put( 'bgmissal' , zerr_sc1(tile_n) / zvol_tot(tile_n) ) ! sc - error due to free surface (psu)
+            ENDIF
+            !
+            !Last step, don't need restart
+            !IF( lrst_oce ) CALL dia_hsb_rst( kts, Kmm, tile_n, 'WRITE' )
+        !
+        END IF
+      IF( ln_timing ) CALL timing_stop('dia_hsb')
+      !
+   END SUBROUTINE dia_hsb
+
+
+   SUBROUTINE dia_hsb_rst( kt, Kmm, tile, cdrw )
+
+
+
+      !!---------------------------------------------------------------------
+      !! *** ROUTINE dia_hsb_rst ***
+      !!
+      !! ** Purpose : Read or write DIA file in restart file
+      !!
+      !! ** Method : use of IOM library
+      !!----------------------------------------------------------------------
+      INTEGER , INTENT(in) :: kt ! ocean time-step
+      INTEGER , INTENT(in) :: Kmm ! ocean time level index
+
+      INTEGER , INTENT(in) :: tile ! host tile
+
+      CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag
+      !
+      INTEGER :: ji, jj, jk ! dummy loop indices
+      !!----------------------------------------------------------------------
+      !
+      IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise
+         IF( ln_rstart ) THEN !* Read the restart file
+            !
+            IF(lwp) WRITE(numout,*)
+            IF(lwp) WRITE(numout,*) '   dia_hsb_rst : read hsb restart at it= ', kt,' date= ', ndastp
+            IF(lwp) WRITE(numout,*)
+
+            CALL iom_get( numror, 'frc_v', frc_v(tile), ldxios = lrxios )
+            CALL iom_get( numror, 'frc_t', frc_t(tile), ldxios = lrxios )
+            CALL iom_get( numror, 'frc_s', frc_s(tile), ldxios = lrxios )
+            IF( ln_linssh ) THEN
+               CALL iom_get( numror, 'frc_wn_t', frc_wn_t(tile), ldxios = lrxios )
+               CALL iom_get( numror, 'frc_wn_s', frc_wn_s(tile), ldxios = lrxios )
+            ENDIF
+            CALL iom_get( numror, jpdom_auto, 'surf_ini' , surf_ini , ldxios = lrxios ) ! ice sheet coupling
+            CALL iom_get( numror, jpdom_auto, 'ssh_ini' , ssh_ini , ldxios = lrxios )
+            CALL iom_get( numror, jpdom_auto, 'e3t_ini' , e3t_ini , ldxios = lrxios )
+            CALL iom_get( numror, jpdom_auto, 'tmask_ini' , tmask_ini , ldxios = lrxios )
+            CALL iom_get( numror, jpdom_auto, 'hc_loc_ini', hc_loc_ini, ldxios = lrxios )
+            CALL iom_get( numror, jpdom_auto, 'sc_loc_ini', sc_loc_ini, ldxios = lrxios )
+            IF( ln_linssh ) THEN
+               CALL iom_get( numror, jpdom_auto, 'ssh_hc_loc_ini', ssh_hc_loc_ini, ldxios = lrxios )
+               CALL iom_get( numror, jpdom_auto, 'ssh_sc_loc_ini', ssh_sc_loc_ini, ldxios = lrxios )
+            ENDIF
+         ELSE
+            IF(lwp) WRITE(numout,*)
+            IF(lwp) WRITE(numout,*) '   dia_hsb_rst : initialise hsb at initial state '
+            IF(lwp) WRITE(numout,*)
+            surf_ini(:,:) = e1e2t(:,:) * tmask_i(:,:) ! initial ocean surface
+            ssh_ini(:,:) = ssh(:,:,Kmm) ! initial ssh
+            DO jk = 1, jpk
+              ! if ice sheet/oceqn coupling, need to mask ini variables here (mask could change at the next NEMO instance).
+               e3t_ini (:,:,jk) = e3t(:,:,jk,Kmm) * tmask(:,:,jk) ! initial vertical scale factors
+               tmask_ini (:,:,jk) = tmask(:,:,jk) ! initial mask
+               hc_loc_ini(:,:,jk) = ts(:,:,jk,jp_tem,Kmm) * e3t(:,:,jk,Kmm) * tmask(:,:,jk) ! initial heat content
+               sc_loc_ini(:,:,jk) = ts(:,:,jk,jp_sal,Kmm) * e3t(:,:,jk,Kmm) * tmask(:,:,jk) ! initial salt content
+            END DO
+
+            d_surf_ini = surf_ini
+            d_e3t_ini = e3t_ini
+            d_tmask_ini = tmask_ini
+            d_hc_loc_ini = hc_loc_ini
+            d_sc_loc_ini = sc_loc_ini
+            frc_v(tile) = 0._wp ! volume trend due to forcing
+            frc_t(tile) = 0._wp ! heat content - - - -
+            frc_s(tile) = 0._wp ! salt content - - - -
+
+
+
+
+
+            IF( ln_linssh ) THEN
+               IF( ln_isfcav ) THEN
+                  DO ji = 1, jpi
+                     DO jj = 1, jpj
+                        ssh_hc_loc_ini(ji,jj) = ts(ji,jj,mikt(ji,jj),jp_tem,Kmm) * ssh(ji,jj,Kmm) ! initial heat content in ssh
+                        ssh_sc_loc_ini(ji,jj) = ts(ji,jj,mikt(ji,jj),jp_sal,Kmm) * ssh(ji,jj,Kmm) ! initial salt content in ssh
+                     END DO
+                   END DO
+                ELSE
+                  ssh_hc_loc_ini(:,:) = ts(:,:,1,jp_tem,Kmm) * ssh(:,:,Kmm) ! initial heat content in ssh
+                  ssh_sc_loc_ini(:,:) = ts(:,:,1,jp_sal,Kmm) * ssh(:,:,Kmm) ! initial salt content in ssh
+               END IF
+
+               frc_wn_t(tile) = 0._wp ! initial heat content misfit due to free surface
+               frc_wn_s(tile) = 0._wp ! initial salt content misfit due to free surface
+
+
+
+
+            ENDIF
+         ENDIF
+         !
+      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file
+         ! ! -------------------
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) '   dia_hsb_rst : write restart at it= ', kt,' date= ', ndastp
+         IF(lwp) WRITE(numout,*)
+         !
+
+         IF( lwxios ) CALL iom_swap( cwxios_context )
+         CALL iom_rstput( kt, nitrst, numrow, 'frc_v', frc_v(tile), ldxios = lwxios )
+         CALL iom_rstput( kt, nitrst, numrow, 'frc_t', frc_t(tile), ldxios = lwxios )
+         CALL iom_rstput( kt, nitrst, numrow, 'frc_s', frc_s(tile), ldxios = lwxios )
+         IF( ln_linssh ) THEN
+            CALL iom_rstput( kt, nitrst, numrow, 'frc_wn_t', frc_wn_t(tile), ldxios = lwxios )
+            CALL iom_rstput( kt, nitrst, numrow, 'frc_wn_s', frc_wn_s(tile), ldxios = lwxios )
+         ENDIF
+         CALL iom_rstput( kt, nitrst, numrow, 'surf_ini' , surf_ini , ldxios = lwxios ) ! ice sheet coupling
+         CALL iom_rstput( kt, nitrst, numrow, 'ssh_ini' , ssh_ini , ldxios = lwxios )
+         CALL iom_rstput( kt, nitrst, numrow, 'e3t_ini' , e3t_ini , ldxios = lwxios )
+         CALL iom_rstput( kt, nitrst, numrow, 'tmask_ini' , tmask_ini , ldxios = lwxios )
+         CALL iom_rstput( kt, nitrst, numrow, 'hc_loc_ini', hc_loc_ini, ldxios = lwxios )
+         CALL iom_rstput( kt, nitrst, numrow, 'sc_loc_ini', sc_loc_ini, ldxios = lwxios )
+         IF( ln_linssh ) THEN
+            CALL iom_rstput( kt, nitrst, numrow, 'ssh_hc_loc_ini', ssh_hc_loc_ini, ldxios = lwxios )
+            CALL iom_rstput( kt, nitrst, numrow, 'ssh_sc_loc_ini', ssh_sc_loc_ini, ldxios = lwxios )
+         ENDIF
+         IF( lwxios ) CALL iom_swap( cxios_context )
+         !
+      ENDIF
+      !
+   END SUBROUTINE dia_hsb_rst
+
+
+   SUBROUTINE dia_hsb_init( Kmm )
+      !!---------------------------------------------------------------------------
+      !! *** ROUTINE dia_hsb ***
+      !!
+      !! ** Purpose: Initialization for the heat salt volume budgets
+      !!
+      !! ** Method : Compute initial heat content, salt content and volume
+      !!
+      !! ** Action : - Compute initial heat content, salt content and volume
+      !! - Initialize forcing trends
+      !! - Compute coefficients for conversion
+      !!---------------------------------------------------------------------------
+      INTEGER, INTENT(in) :: Kmm ! time level index
+      !
+      INTEGER :: ierror, ios ! local integer
+
+      INTEGER :: i, istat ! local integer
+
+      !!
+      NAMELIST/namhsb/ ln_diahsb
+      !!----------------------------------------------------------------------
+      !
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'dia_hsb_init : heat and salt budgets diagnostics'
+         WRITE(numout,*) '~~~~~~~~~~~~ '
+      ENDIF
+      READ ( numnam_ref, namhsb, IOSTAT = ios, ERR = 901)
+901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namhsb in reference namelist' )
+      READ ( numnam_cfg, namhsb, IOSTAT = ios, ERR = 902 )
+902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namhsb in configuration namelist' )
+
+      IF(lwm) WRITE( numond, namhsb )
+
+      IF(lwp) THEN
+         WRITE(numout,*) '   Namelist  namhsb :'
+         WRITE(numout,*) '      check the heat and salt budgets (T) or not (F)       ln_diahsb = ', ln_diahsb
+      ENDIF
+      !
+      IF( .NOT. ln_diahsb ) RETURN
+
+      IF(lwxios) THEN
+! define variables in restart file when writing with XIOS
+        CALL iom_set_rstw_var_active('frc_v')
+        CALL iom_set_rstw_var_active('frc_t')
+        CALL iom_set_rstw_var_active('frc_s')
+        CALL iom_set_rstw_var_active('surf_ini')
+        CALL iom_set_rstw_var_active('ssh_ini')
+        CALL iom_set_rstw_var_active('e3t_ini')
+        CALL iom_set_rstw_var_active('hc_loc_ini')
+        CALL iom_set_rstw_var_active('sc_loc_ini')
+        IF( ln_linssh ) THEN
+           CALL iom_set_rstw_var_active('ssh_hc_loc_ini')
+           CALL iom_set_rstw_var_active('ssh_sc_loc_ini')
+           CALL iom_set_rstw_var_active('frc_wn_t')
+           CALL iom_set_rstw_var_active('frc_wn_s')
+        ENDIF
+      ENDIF
+      ! ------------------- !
+      ! 1 - Allocate memory !
+      ! ------------------- !
+
+      CALL setdevice()
+      !Device data associate to PUBLIC arrays
+      ALLOCATE(d_e3t (jpi,jpj,jpk,jpt) ) !
+      ALLOCATE(d_tmask (jpi,jpj,jpk) ) !
+      ALLOCATE(d_tmask_ini (jpi,jpj,jpk) ) !
+      ALLOCATE(d_tmask_h (jpi,jpj) ) !
+      ALLOCATE(d_ts (jpi,jpj,jpk,2,jpj) ) !
+      !Device data associate to LOCAL/DEVICE arrays !
+      ALLOCATE(d_surf (jpi,jpj) ) !
+      ALLOCATE(d_surf_ini (jpi,jpj) ) !
+      ALLOCATE(d_hc_loc_ini (jpi,jpj,jpk) ) !
+      ALLOCATE(d_sc_loc_ini (jpi,jpj,jpk) ) !
+      ALLOCATE(d_e3t_ini (jpi,jpj,jpk) ) !
+      ALLOCATE(d_zwrkv (jpi,jpj,jpkm1) ) !
+      ALLOCATE(d_zwrkh (jpi,jpj,jpkm1) ) !
+      ALLOCATE(d_zwrks (jpi,jpj,jpkm1) ) !
+      ALLOCATE(d_zwrk (jpi,jpj,jpkm1) ) !
+      ALLOCATE(h_ztmpv(2),h_ztmph(2),h_ztmps(2),h_ztmp(2)) !
+
+      DO i = 1, nstreams !Create Streams
+          istat = cudaStreamCreate(stream(i))
+         IF( istat /= 0 ) THEN
+         CALL ctl_stop( 'dia_hsb_init: error in Stream creation' ) ; RETURN
+         ENDIF
+      END DO
+            !
+      !Pinned reallocation step non constant
+      istat = cudaHostRegister(C_LOC(ts ), sizeof(ts ), cudaHostRegisterMapped)
+      istat = cudaHostRegister(C_LOC(e3t), sizeof(e3t), cudaHostRegisterMapped)
+      IF( istat /= 0 ) THEN
+         CALL ctl_stop( 'dia_hsb_init: unable to pin host memory to GPU' ) ; RETURN
+      ENDIF
+
+      ALLOCATE( hc_loc_ini(jpi,jpj,jpk), sc_loc_ini(jpi,jpj,jpk), surf_ini(jpi,jpj), &
+         & e3t_ini(jpi,jpj,jpk), surf(jpi,jpj), ssh_ini(jpi,jpj), tmask_ini(jpi,jpj,jpk),STAT=ierror )
+
+      IF( ierror > 0 ) THEN
+         CALL ctl_stop( 'dia_hsb_init: unable to allocate hc_loc_ini' ) ; RETURN
+      ENDIF
+
+      IF( ln_linssh ) ALLOCATE( ssh_hc_loc_ini(jpi,jpj), ssh_sc_loc_ini(jpi,jpj),STAT=ierror )
+      IF( ierror > 0 ) THEN
+         CALL ctl_stop( 'dia_hsb: unable to allocate ssh_hc_loc_ini' ) ; RETURN
+      ENDIF
+
+      ! ----------------------------------------------- !
+      ! 2 - Time independant variables and file opening !
+      ! ----------------------------------------------- !
+      surf(:,:) = e1e2t(:,:) * tmask_i(:,:) ! masked surface grid cell area
+      surf_tot = glob_sum( 'diahsb', surf(:,:) ) ! total ocean surface area
+
+       d_surf = surf
+       d_surf_ini = surf_ini
+       d_e3t_ini = e3t_ini
+       d_tmask = tmask
+       d_tmask_ini = tmask_ini
+       d_tmask_h = tmask_h
+       h_ztmp = 0.0
+
+      IF( ln_bdy ) CALL ctl_warn( 'dia_hsb_init: heat/salt budget does not consider open boundary fluxes' )
+      !
+      ! ---------------------------------- !
+      ! 4 - initial conservation variables !
+      ! ---------------------------------- !
+
+      CALL dia_hsb_rst( nit000, Kmm, 1, 'READ' ) !* read or initialize all required files
+
+
+
+      !
+   END SUBROUTINE dia_hsb_init
+
+   !!======================================================================
+END MODULE diahsb
Index: NEMO/branches/2020/dev_r13747_HPC-10_HPDAonline_DiagGPU/tests/DIA_GPU/MY_SRC/gpu_manager.F90
===================================================================
--- NEMO/branches/2020/dev_r13747_HPC-10_HPDAonline_DiagGPU/tests/DIA_GPU/MY_SRC/gpu_manager.F90	(revision 14091)
+++ NEMO/branches/2020/dev_r13747_HPC-10_HPDAonline_DiagGPU/tests/DIA_GPU/MY_SRC/gpu_manager.F90	(revision 14091)
@@ -0,0 +1,34 @@
+#ifdef key_gpu
+MODULE gpu_manager
+    USE lib_mpp        ! distributed memory computing library
+    USE cudafor
+    USE par_oce
+    !!======================================================================
+    !!                       ***  MODULE  gpu_manager  ***
+    !! GPU CUDA utilities : Defines run parameters for GPGPU routines
+    !!=====================================================================
+    !! History :   1.0  !  2019-10  (M. Faria)   original code
+    !!----------------------------------------------------------------------
+
+    !!----------------------------------------------------------------------
+    IMPLICIT NONE
+    PUBLIC
+    !!----------------------------------------------------------------------
+    !!                   gpu configuration parameters
+    !!----------------------------------------------------------------------
+    INTEGER          :: cudaistat, usecudev
+    INTEGER          :: nodesize
+    INTEGER          :: cudadev
+    INTEGER          :: cusocket
+CONTAINS
+    SUBROUTINE setdevice()
+        nodesize  = 1                                                    !Node number of MPI processes, CTE-Power default 40
+        cusocket  = 1                                                    !Number of NVlink or PCI-E nodes, CTE-Power default 2
+        cudadev   = 1                                                    !Number of GPUs on a node, CTE-Power default 4
+        usecudev  = MOD( mpprank /  ( nodesize / cusocket )  , cudadev ) !Bind MPI RANK to GPU, consider balanced NVlink and constant number of devices per node 
+        cudaistat = cudaSetDevice( usecudev  )                           !Set which GPU use. Ex for 4 GPUs. cpus_per_socket / sockets_per_node for "bind-to core"
+        print *, 'CUDA istat ', cudaistat, 'rank', mpprank, 'cuda device', cudadev 
+    END SUBROUTINE setdevice
+
+END MODULE gpu_manager
+#endif
Index: NEMO/branches/2020/dev_r13747_HPC-10_HPDAonline_DiagGPU/tests/DIA_GPU/MY_SRC/nvtx.F90
===================================================================
--- NEMO/branches/2020/dev_r13747_HPC-10_HPDAonline_DiagGPU/tests/DIA_GPU/MY_SRC/nvtx.F90	(revision 14091)
+++ NEMO/branches/2020/dev_r13747_HPC-10_HPDAonline_DiagGPU/tests/DIA_GPU/MY_SRC/nvtx.F90	(revision 14091)
@@ -0,0 +1,66 @@
+#ifdef key_gpu
+module nvtx
+
+use iso_c_binding
+implicit none
+
+integer,private :: col(7) = [ Z'0000ff00', Z'000000ff', Z'00ffff00', Z'00ff00ff', Z'0000ffff', Z'00ff0000', Z'00ffffff']
+character(len=256),private :: tempName
+
+type, bind(C):: nvtxEventAttributes
+  integer(C_INT16_T):: version=1
+  integer(C_INT16_T):: size=48 !
+  integer(C_INT):: category=0
+  integer(C_INT):: colorType=1 ! NVTX_COLOR_ARGB = 1
+  integer(C_INT):: color
+  integer(C_INT):: payloadType=0 ! NVTX_PAYLOAD_UNKNOWN = 0
+  integer(C_INT):: reserved0
+  integer(C_INT64_T):: payload   ! union uint,int,double
+  integer(C_INT):: messageType=1  ! NVTX_MESSAGE_TYPE_ASCII     = 1 
+  type(C_PTR):: message  ! ascii char
+end type
+
+interface nvtxRangePush
+  ! push range with custom label and standard color
+  subroutine nvtxRangePushA(name) bind(C, name='nvtxRangePushA')
+  use iso_c_binding
+  character(kind=C_CHAR,len=*) :: name
+  end subroutine
+
+  ! push range with custom label and custom color
+  subroutine nvtxRangePushEx(event) bind(C, name='nvtxRangePushEx')
+  use iso_c_binding
+  import:: nvtxEventAttributes
+  type(nvtxEventAttributes):: event
+  end subroutine
+end interface
+
+interface nvtxRangePop
+  subroutine nvtxRangePop() bind(C, name='nvtxRangePop')
+  end subroutine
+end interface
+
+contains
+
+subroutine nvtxStartRange(name,id)
+  character(kind=c_char,len=*) :: name
+  integer, optional:: id
+  type(nvtxEventAttributes):: event
+
+  tempName=trim(name)//c_null_char
+
+  if ( .not. present(id)) then
+    call nvtxRangePush(tempName)
+  else
+    event%color=col(mod(id,7)+1)
+    event%message=c_loc(tempName)
+    call nvtxRangePushEx(event)
+  end if
+end subroutine
+
+subroutine nvtxEndRange
+  call nvtxRangePop
+end subroutine
+
+end module nvtx
+#endif
Index: NEMO/branches/2020/dev_r13747_HPC-10_HPDAonline_DiagGPU/tests/DIA_GPU/MY_SRC/stpctl.F90
===================================================================
--- NEMO/branches/2020/dev_r13747_HPC-10_HPDAonline_DiagGPU/tests/DIA_GPU/MY_SRC/stpctl.F90	(revision 14091)
+++ NEMO/branches/2020/dev_r13747_HPC-10_HPDAonline_DiagGPU/tests/DIA_GPU/MY_SRC/stpctl.F90	(revision 14091)
@@ -0,0 +1,321 @@
+MODULE stpctl
+   !!======================================================================
+   !!                       ***  MODULE  stpctl  ***
+   !! Ocean run control :  gross check of the ocean time stepping
+   !!======================================================================
+   !! History :  OPA  ! 1991-03  (G. Madec) Original code
+   !!            6.0  ! 1992-06  (M. Imbard)
+   !!            8.0  ! 1997-06  (A.M. Treguier)
+   !!   NEMO     1.0  ! 2002-06  (G. Madec)  F90: Free form and module
+   !!            2.0  ! 2009-07  (G. Madec)  Add statistic for time-spliting
+   !!            3.7  ! 2016-09  (G. Madec)  Remove solver
+   !!            4.0  ! 2017-04  (G. Madec)  regroup global communications
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   stp_ctl      : Control the run
+   !!----------------------------------------------------------------------
+   USE oce             ! ocean dynamics and tracers variables
+   USE dom_oce         ! ocean space and time domain variables 
+   USE c1d             ! 1D vertical configuration
+   USE zdf_oce ,  ONLY : ln_zad_Aimp       ! ocean vertical physics variables
+   USE wet_dry,   ONLY : ll_wd, ssh_ref    ! reference depth for negative bathy
+   !  
+   USE diawri          ! Standard run outputs       (dia_wri_state routine)
+   USE in_out_manager  ! I/O manager
+   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
+   USE lib_mpp         ! distributed memory computing
+   !
+   USE netcdf          ! NetCDF library
+   USE, INTRINSIC :: IEEE_ARITHMETIC, ONLY : IEEE_IS_NAN
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC stp_ctl           ! routine called by step.F90
+
+   INTEGER                ::   nrunid   ! netcdf file id
+   INTEGER, DIMENSION(8)  ::   nvarid   ! netcdf variable id
+   !!----------------------------------------------------------------------
+   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+   !! $Id$
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE stp_ctl( kt, Kmm )
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE stp_ctl  ***
+      !!                     
+      !! ** Purpose :   Control the run
+      !!
+      !! ** Method  : - Save the time step in numstp
+      !!              - Stop the run IF problem encountered by setting nstop > 0
+      !!                Problems checked: |ssh| maximum larger than 10 m
+      !!                                  |U|   maximum larger than 10 m/s 
+      !!                                  negative sea surface salinity
+      !!
+      !! ** Actions :   "time.step" file = last ocean time-step
+      !!                "run.stat"  file = run statistics
+      !!                 nstop indicator sheared among all local domain
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in   ) ::   kt       ! ocean time-step index
+      INTEGER, INTENT(in   ) ::   Kmm      ! ocean time level index
+      !!
+      INTEGER                         ::   ji                                    ! dummy loop indices
+      INTEGER                         ::   idtime, istatus
+      INTEGER , DIMENSION(9)          ::   iareasum, iareamin, iareamax
+      INTEGER , DIMENSION(3,4)        ::   iloc                                  ! min/max loc indices
+      REAL(wp)                        ::   zzz                                   ! local real 
+      REAL(wp), DIMENSION(9)          ::   zmax, zmaxlocal
+      LOGICAL                         ::   ll_wrtstp, ll_colruns, ll_wrtruns, ll_0oce
+      LOGICAL, DIMENSION(jpi,jpj,jpk) ::   llmsk
+      CHARACTER(len=20)               ::   clname
+      !!----------------------------------------------------------------------
+      IF( nstop > 0 .AND. ngrdstop > -1 )   RETURN   !   stpctl was already called by a child grid
+      !
+      ll_wrtstp  = ( MOD( kt-nit000, sn_cfctl%ptimincr ) == 0 ) .OR. ( kt == nitend )
+      ll_colruns = ll_wrtstp .AND. sn_cfctl%l_runstat .AND. jpnij > 1 
+      ll_wrtruns = ( ll_colruns .OR. jpnij == 1 ) .AND. lwm
+      !
+      IF( kt == nit000 ) THEN
+         !
+         IF( lwp ) THEN
+            WRITE(numout,*)
+            WRITE(numout,*) 'stp_ctl : time-stepping control'
+            WRITE(numout,*) '~~~~~~~'
+         ENDIF
+         !                                ! open time.step    ascii file, done only by 1st subdomain
+         IF( lwm )   CALL ctl_opn( numstp, 'time.step', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
+         !
+         IF( ll_wrtruns ) THEN
+            !                             ! open run.stat     ascii file, done only by 1st subdomain
+            CALL ctl_opn( numrun, 'run.stat', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
+            !                             ! open run.stat.nc netcdf file, done only by 1st subdomain
+            clname = 'run.stat.nc'
+            IF( .NOT. Agrif_Root() )   clname = TRIM(Agrif_CFixed())//"_"//TRIM(clname)
+            istatus = NF90_CREATE( TRIM(clname), NF90_CLOBBER, nrunid )
+            istatus = NF90_DEF_DIM( nrunid, 'time', NF90_UNLIMITED, idtime )
+            istatus = NF90_DEF_VAR( nrunid, 'abs_ssh_max', NF90_DOUBLE, (/ idtime /), nvarid(1) )
+            istatus = NF90_DEF_VAR( nrunid,   'abs_u_max', NF90_DOUBLE, (/ idtime /), nvarid(2) )
+            istatus = NF90_DEF_VAR( nrunid,       's_min', NF90_DOUBLE, (/ idtime /), nvarid(3) )
+            istatus = NF90_DEF_VAR( nrunid,       's_max', NF90_DOUBLE, (/ idtime /), nvarid(4) )
+            istatus = NF90_DEF_VAR( nrunid,       't_min', NF90_DOUBLE, (/ idtime /), nvarid(5) )
+            istatus = NF90_DEF_VAR( nrunid,       't_max', NF90_DOUBLE, (/ idtime /), nvarid(6) )
+            IF( ln_zad_Aimp ) THEN
+               istatus = NF90_DEF_VAR( nrunid,   'Cf_max', NF90_DOUBLE, (/ idtime /), nvarid(7) )
+               istatus = NF90_DEF_VAR( nrunid,'abs_wi_max',NF90_DOUBLE, (/ idtime /), nvarid(8) )
+            ENDIF
+            istatus = NF90_ENDDEF(nrunid)
+         ENDIF
+         !    
+      ENDIF
+      !
+      !                                   !==              write current time step              ==!
+      !                                   !==  done only by 1st subdomain at writting timestep  ==!
+      IF( lwm .AND. ll_wrtstp ) THEN
+         WRITE ( numstp, '(1x, i8)' )   kt
+         REWIND( numstp )
+      ENDIF
+      !                                   !==            test of local extrema           ==!
+      !                                   !==  done by all processes at every time step  ==!
+      !
+      llmsk(   1:Nis1,:,:) = .FALSE.                                              ! exclude halos from the checked region
+      llmsk(Nie1: jpi,:,:) = .FALSE.
+      llmsk(:,   1:Njs1,:) = .FALSE.
+      llmsk(:,Nje1: jpj,:) = .FALSE.
+      !
+      llmsk(Nis0:Nie0,Njs0:Nje0,1) = ssmask(Nis0:Nie0,Njs0:Nje0) == 1._wp         ! define only the inner domain
+      !
+      ll_0oce = .NOT. ANY( llmsk(:,:,1) )                                         ! no ocean point in the inner domain?
+      !
+      IF( ll_wd ) THEN
+         zmax(1) = MAXVAL( ABS( ssh(:,:,Kmm) + ssh_ref ), mask = llmsk(:,:,1) )   ! ssh max
+      ELSE
+         zmax(1) = MAXVAL( ABS( ssh(:,:,Kmm)           ), mask = llmsk(:,:,1) )   ! ssh max
+      ENDIF
+      llmsk(Nis0:Nie0,Njs0:Nje0,:) = umask(Nis0:Nie0,Njs0:Nje0,:) == 1._wp        ! define only the inner domain
+      zmax(2) = MAXVAL(  ABS( uu(:,:,:,Kmm) ), mask = llmsk )                     ! velocity max (zonal only)
+      llmsk(Nis0:Nie0,Njs0:Nje0,:) = tmask(Nis0:Nie0,Njs0:Nje0,:) == 1._wp        ! define only the inner domain
+      zmax(3) = MAXVAL( -ts(:,:,:,jp_sal,Kmm), mask = llmsk )                     ! minus salinity max
+      zmax(4) = MAXVAL(  ts(:,:,:,jp_sal,Kmm), mask = llmsk )                     !       salinity max
+      IF( ll_colruns .OR. jpnij == 1 ) THEN     ! following variables are used only in the netcdf file
+         zmax(5) = MAXVAL( -ts(:,:,:,jp_tem,Kmm), mask = llmsk )                  ! minus temperature max
+         zmax(6) = MAXVAL(  ts(:,:,:,jp_tem,Kmm), mask = llmsk )                  !       temperature max
+         IF( ln_zad_Aimp ) THEN
+            zmax(7) = MAXVAL(   Cu_adv(:,:,:)   , mask = llmsk )                  ! partitioning coeff. max
+            llmsk(:,:,:) = wmask(:,:,:) == 1._wp
+            zmax(8) = MAXVAL(  ABS( wi(:,:,:) ) , mask = llmsk )                  ! implicit vertical vel. max
+         ELSE
+            zmax(7:8) = 0._wp
+         ENDIF
+      ELSE
+         zmax(5:8) = 0._wp
+      ENDIF
+      zmax(9) = REAL( nstop, wp )                                                 ! stop indicator
+      !
+      !                                   !==               get global extrema             ==!
+      !                                   !==  done by all processes if writting run.stat  ==!
+      IF( ll_colruns ) THEN
+         zmaxlocal(:) = zmax(:)
+         CALL mpp_max( "stpctl", zmax )          ! max over the global domain: ok even of ll_0oce = .true. 
+         nstop = NINT( zmax(9) )                 ! update nstop indicator (now sheared among all local domains)
+      ELSE
+         ! if no ocean point: MAXVAL returns -HUGE => we must overwrite this value to avoid error handling bellow.
+         IF( ll_0oce )   zmax(1:4) = (/ 0._wp, 0._wp, -1._wp, 1._wp /)   ! default "valid" values...
+      ENDIF
+      !
+      zmax(3) = -zmax(3)                         ! move back from max(-zz) to min(zz) : easier to manage! 
+      zmax(5) = -zmax(5)                         ! move back from max(-zz) to min(zz) : easier to manage!
+      IF( ll_colruns ) THEN
+         zmaxlocal(3) = -zmaxlocal(3)            ! move back from max(-zz) to min(zz) : easier to manage! 
+         zmaxlocal(5) = -zmaxlocal(5)            ! move back from max(-zz) to min(zz) : easier to manage!
+      ENDIF
+      !
+      !                                   !==              write "run.stat" files              ==!
+      !                                   !==  done only by 1st subdomain at writting timestep  ==!
+      IF( ll_wrtruns ) THEN
+         WRITE(numrun,9500) kt, zmax(1), zmax(2), zmax(3), zmax(4)
+         DO ji = 1, 6 + 2 * COUNT( (/ln_zad_Aimp/) )
+            istatus = NF90_PUT_VAR( nrunid, nvarid(ji), (/zmax(ji)/), (/kt/), (/1/) )
+         END DO
+         IF( kt == nitend )   istatus = NF90_CLOSE(nrunid)
+      END IF
+      !                                   !==               error handling               ==!
+      !                                   !==  done by all processes at every time step  ==!
+      !
+      IF(  zmax(1) >   20._wp .OR.   &                   ! too large sea surface height ( > 20 m )
+         & zmax(2) >   10._wp .OR.   &                   ! too large velocity ( > 10 m/s)
+         & zmax(3) <=   0._wp .OR.   &                   ! negative or zero sea surface salinity
+         & zmax(4) >= 100._wp .OR.   &                   ! too large sea surface salinity ( > 100 )
+         & zmax(4) <    0._wp .OR.   &                   ! too large sea surface salinity (keep this line for sea-ice)
+         & IEEE_IS_NAN( zmax(1) + zmax(2) + zmax(3) ) .OR.   &               ! NaN encounter in the tests
+         & ABS(   zmax(1) + zmax(2) + zmax(3) ) > HUGE(1._wp) ) THEN   ! Infinity encounter in the tests
+         !
+         iloc(:,:) = 0
+         IF( ll_colruns ) THEN   ! zmax is global, so it is the same on all subdomains -> no dead lock with mpp_maxloc
+            ! first: close the netcdf file, so we can read it
+            IF( lwm .AND. kt /= nitend )   istatus = NF90_CLOSE(nrunid)
+            ! get global loc on the min/max
+            llmsk(Nis0:Nie0,Njs0:Nje0,1) = ssmask(Nis0:Nie0,Njs0:Nje0 ) == 1._wp         ! define only the inner domain
+            CALL mpp_maxloc( 'stpctl', ABS(ssh(:,:,         Kmm)), llmsk(:,:,1), zzz, iloc(1:2,1) )   ! mpp_maxloc ok if mask = F 
+            llmsk(Nis0:Nie0,Njs0:Nje0,:) = umask(Nis0:Nie0,Njs0:Nje0,:) == 1._wp        ! define only the inner domain
+            CALL mpp_maxloc( 'stpctl', ABS( uu(:,:,:,       Kmm)), llmsk(:,:,:), zzz, iloc(1:3,2) )
+            llmsk(Nis0:Nie0,Njs0:Nje0,:) = tmask(Nis0:Nie0,Njs0:Nje0,:) == 1._wp        ! define only the inner domain
+            CALL mpp_minloc( 'stpctl',      ts(:,:,:,jp_sal,Kmm) , llmsk(:,:,:), zzz, iloc(1:3,3) )
+            CALL mpp_maxloc( 'stpctl',      ts(:,:,:,jp_sal,Kmm) , llmsk(:,:,:), zzz, iloc(1:3,4) )
+            ! find which subdomain has the max.
+            iareamin(:) = jpnij+1   ;   iareamax(:) = 0   ;   iareasum(:) = 0
+            DO ji = 1, 9
+               IF( zmaxlocal(ji) == zmax(ji) ) THEN
+                  iareamin(ji) = narea   ;   iareamax(ji) = narea   ;   iareasum(ji) = 1
+               ENDIF
+            END DO
+            CALL mpp_min( "stpctl", iareamin )         ! min over the global domain
+            CALL mpp_max( "stpctl", iareamax )         ! max over the global domain
+            CALL mpp_sum( "stpctl", iareasum )         ! sum over the global domain
+         ELSE                    ! find local min and max locations:
+            ! if we are here, this means that the subdomain contains some oce points -> no need to test the mask used in maxloc
+            llmsk(Nis0:Nie0,Njs0:Nje0,1) = ssmask(Nis0:Nie0,Njs0:Nje0 ) == 1._wp        ! define only the inner domain
+            iloc(1:2,1) = MAXLOC( ABS( ssh(:,:,         Kmm)), mask = llmsk(:,:,1) )
+            llmsk(Nis0:Nie0,Njs0:Nje0,:) = umask(Nis0:Nie0,Njs0:Nje0,:) == 1._wp        ! define only the inner domain
+            iloc(1:3,2) = MAXLOC( ABS(  uu(:,:,:,       Kmm)), mask = llmsk(:,:,:) )
+            llmsk(Nis0:Nie0,Njs0:Nje0,:) = tmask(Nis0:Nie0,Njs0:Nje0,:) == 1._wp        ! define only the inner domain
+            iloc(1:3,3) = MINLOC(       ts(:,:,:,jp_sal,Kmm) , mask = llmsk(:,:,:) )
+            iloc(1:3,4) = MAXLOC(       ts(:,:,:,jp_sal,Kmm) , mask = llmsk(:,:,:) )
+            DO ji = 1, 4   ! local domain indices ==> global domain indices, excluding halos
+               iloc(1:2,ji) = (/ mig0(iloc(1,ji)), mjg0(iloc(2,ji)) /)
+            END DO
+            iareamin(:) = narea   ;   iareamax(:) = narea   ;   iareasum(:) = 1         ! this is local information
+         ENDIF
+         !
+         WRITE(ctmp1,*) ' stp_ctl: |ssh| > 20 m  or  |U| > 10 m/s  or  S <= 0  or  S >= 100  or  NaN encounter in the tests'
+         CALL wrt_line( ctmp2, kt, '|ssh| max', zmax(1), iloc(:,1), iareasum(1), iareamin(1), iareamax(1) )
+         CALL wrt_line( ctmp3, kt, '|U|   max', zmax(2), iloc(:,2), iareasum(2), iareamin(2), iareamax(2) )
+         CALL wrt_line( ctmp4, kt, 'Sal   min', zmax(3), iloc(:,3), iareasum(3), iareamin(3), iareamax(3) )
+         CALL wrt_line( ctmp5, kt, 'Sal   max', zmax(4), iloc(:,4), iareasum(4), iareamin(4), iareamax(4) )
+         IF( Agrif_Root() ) THEN
+            WRITE(ctmp6,*) '      ===> output of last computed fields in output.abort* files'
+         ELSE
+            WRITE(ctmp6,*) '      ===> output of last computed fields in '//TRIM(Agrif_CFixed())//'_output.abort* files'
+         ENDIF
+         !
+         CALL dia_wri_state( Kmm, 'output.abort' )     ! create an output.abort file
+         !
+         IF( ll_colruns .or. jpnij == 1 ) THEN   ! all processes synchronized -> use lwp to print in opened ocean.output files
+            IF(lwp) THEN   ;   CALL ctl_stop( ctmp1, ' ', ctmp2, ctmp3, ctmp4, ctmp5, ' ', ctmp6 )
+            ELSE           ;   nstop = MAX(1, nstop)   ! make sure nstop > 0 (automatically done when calling ctl_stop)
+            ENDIF
+         ELSE                                    ! only mpi subdomains with errors are here -> STOP now
+            CALL ctl_stop( 'STOP', ctmp1, ' ', ctmp2, ctmp3, ctmp4, ctmp5, ' ', ctmp6 )
+         ENDIF
+         !
+      ENDIF
+      !
+      IF( nstop > 0 ) THEN                                                  ! an error was detected and we did not abort yet...
+         ngrdstop = Agrif_Fixed()                                           ! store which grid got this error
+         IF( .NOT. ll_colruns .AND. jpnij > 1 )   CALL ctl_stop( 'STOP' )   ! we must abort here to avoid MPI deadlock
+      ENDIF
+      !
+9500  FORMAT(' it :', i8, '    |ssh|_max: ', D23.16, ' |U|_max: ', D23.16,' S_min: ', D23.16,' S_max: ', D23.16)
+      !
+   END SUBROUTINE stp_ctl
+
+
+   SUBROUTINE wrt_line( cdline, kt, cdprefix, pval, kloc, ksum, kmin, kmax )
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE wrt_line  ***
+      !!
+      !! ** Purpose :   write information line
+      !!
+      !!----------------------------------------------------------------------
+      CHARACTER(len=*),      INTENT(  out) ::   cdline
+      CHARACTER(len=*),      INTENT(in   ) ::   cdprefix
+      REAL(wp),              INTENT(in   ) ::   pval
+      INTEGER, DIMENSION(3), INTENT(in   ) ::   kloc
+      INTEGER,               INTENT(in   ) ::   kt, ksum, kmin, kmax
+      !
+      CHARACTER(len=80) ::   clsuff
+      CHARACTER(len=9 ) ::   clkt, clsum, clmin, clmax
+      CHARACTER(len=9 ) ::   cli, clj, clk
+      CHARACTER(len=1 ) ::   clfmt
+      CHARACTER(len=4 ) ::   cl4   ! needed to be able to compile with Agrif, I don't know why
+      INTEGER           ::   ifmtk
+      !!----------------------------------------------------------------------
+      WRITE(clkt , '(i9)') kt
+      
+      WRITE(clfmt, '(i1)') INT(LOG10(REAL(jpnij  ,wp))) + 1     ! how many digits to we need to write ? (we decide max = 9)
+      !!! WRITE(clsum, '(i'//clfmt//')') ksum                   ! this is creating a compilation error with AGRIF
+      cl4 = '(i'//clfmt//')'   ;   WRITE(clsum, cl4) ksum
+      WRITE(clfmt, '(i1)') INT(LOG10(REAL(MAX(1,jpnij-1),wp))) + 1    ! how many digits to we need to write ? (we decide max = 9)
+      cl4 = '(i'//clfmt//')'   ;   WRITE(clmin, cl4) kmin-1
+                                   WRITE(clmax, cl4) kmax-1
+      !
+      WRITE(clfmt, '(i1)') INT(LOG10(REAL(jpiglo,wp))) + 1      ! how many digits to we need to write jpiglo? (we decide max = 9)
+      cl4 = '(i'//clfmt//')'   ;   WRITE(cli, cl4) kloc(1)      ! this is ok with AGRIF
+      WRITE(clfmt, '(i1)') INT(LOG10(REAL(jpjglo,wp))) + 1      ! how many digits to we need to write jpjglo? (we decide max = 9)
+      cl4 = '(i'//clfmt//')'   ;   WRITE(clj, cl4) kloc(2)      ! this is ok with AGRIF
+      !
+      IF( ksum == 1 ) THEN   ;   WRITE(clsuff,9100) TRIM(clmin)
+      ELSE                   ;   WRITE(clsuff,9200) TRIM(clsum), TRIM(clmin), TRIM(clmax)
+      ENDIF
+      IF(kloc(3) == 0) THEN
+         ifmtk = INT(LOG10(REAL(jpk,wp))) + 1                   ! how many digits to we need to write jpk? (we decide max = 9)
+         clk = REPEAT(' ', ifmtk)                               ! create the equivalent in blank string
+         WRITE(cdline,9300) TRIM(ADJUSTL(clkt)), TRIM(ADJUSTL(cdprefix)), pval, TRIM(cli), TRIM(clj), clk(1:ifmtk), TRIM(clsuff)
+      ELSE
+         WRITE(clfmt, '(i1)') INT(LOG10(REAL(jpk,wp))) + 1      ! how many digits to we need to write jpk? (we decide max = 9)
+         !!! WRITE(clk, '(i'//clfmt//')') kloc(3)               ! this is creating a compilation error with AGRIF
+         cl4 = '(i'//clfmt//')'   ;   WRITE(clk, cl4) kloc(3)   ! this is ok with AGRIF
+         WRITE(cdline,9400) TRIM(ADJUSTL(clkt)), TRIM(ADJUSTL(cdprefix)), pval, TRIM(cli), TRIM(clj),    TRIM(clk), TRIM(clsuff)
+      ENDIF
+      !
+9100  FORMAT('MPI rank ', a)
+9200  FORMAT('found in ', a, ' MPI tasks, spread out among ranks ', a, ' to ', a)
+9300  FORMAT('kt ', a, ' ', a, ' ', 1pg11.4, ' at i j   ', a, ' ', a, ' ', a, ' ', a)
+9400  FORMAT('kt ', a, ' ', a, ' ', 1pg11.4, ' at i j k ', a, ' ', a, ' ', a, ' ', a)
+      !
+   END SUBROUTINE wrt_line
+
+
+   !!======================================================================
+END MODULE stpctl
Index: NEMO/branches/2020/dev_r13747_HPC-10_HPDAonline_DiagGPU/tests/README.rst
===================================================================
--- NEMO/branches/2020/dev_r13747_HPC-10_HPDAonline_DiagGPU/tests/README.rst	(revision 14085)
+++ NEMO/branches/2020/dev_r13747_HPC-10_HPDAonline_DiagGPU/tests/README.rst	(revision 14091)
@@ -210,2 +210,9 @@
 | This test case checks the OASIS interface in OCE/SBC, allowing to set up 
 a coupled configuration through OASIS. See CPL_OASIS/README.md for more information.
+
+DIA_GPU
+---------
+| This is a demonstrator of diagnostic DIAHSB ported to GPU using CUDA Fortran. 
+Memory communications between host and device are asynchronous given the device has that capability. 
+This experiment is target for ORCA2_ICE_PISCES
+
Index: NEMO/branches/2020/dev_r13747_HPC-10_HPDAonline_DiagGPU/tests/demo_cfgs.txt
===================================================================
--- NEMO/branches/2020/dev_r13747_HPC-10_HPDAonline_DiagGPU/tests/demo_cfgs.txt	(revision 14085)
+++ NEMO/branches/2020/dev_r13747_HPC-10_HPDAonline_DiagGPU/tests/demo_cfgs.txt	(revision 14091)
@@ -12,2 +12,3 @@
 STATION_ASF OCE
 CPL_OASIS  OCE TOP ICE NST
+DIA_GPU OCE TOP ICE NST ABL
