Index: NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/DIA/diaar5.F90
===================================================================
--- NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/DIA/diaar5.F90	(revision 12529)
+++ NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/DIA/diaar5.F90	(revision 13151)
@@ -32,5 +32,4 @@
    REAL(wp)                         ::   vol0         ! ocean volume (interior domain)
    REAL(wp)                         ::   area_tot     ! total ocean surface (interior domain)
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:  ) ::   area         ! cell surface (interior domain)
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:  ) ::   thick0       ! ocean thickness (interior domain)
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sn0          ! initial salinity
@@ -40,4 +39,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -54,5 +54,5 @@
       !!----------------------------------------------------------------------
       !
-      ALLOCATE( area(jpi,jpj), thick0(jpi,jpj) , sn0(jpi,jpj,jpk) , STAT=dia_ar5_alloc )
+      ALLOCATE( thick0(jpi,jpj) , sn0(jpi,jpj,jpk) , STAT=dia_ar5_alloc )
       !
       CALL mpp_sum ( 'diaar5', dia_ar5_alloc )
@@ -77,6 +77,6 @@
       !
       REAL(wp), ALLOCATABLE, DIMENSION(:,:)     :: zarea_ssh , zbotpres       ! 2D workspace 
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:)     :: zpe, z2d                   ! 2D workspace 
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   :: zrhd , zrhop, ztpot   ! 3D workspace
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:)     :: z2d, zpe                   ! 2D workspace 
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   :: z3d, zrhd , zrhop, ztpot, zgdept   ! 3D workspace (zgdept: needed to use the substitute)
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) :: ztsn                       ! 4D workspace
 
@@ -90,18 +90,21 @@
          ALLOCATE( zrhd(jpi,jpj,jpk) , zrhop(jpi,jpj,jpk) )
          ALLOCATE( ztsn(jpi,jpj,jpk,jpts) )
-         zarea_ssh(:,:) = area(:,:) * ssh(:,:,Kmm)
-      ENDIF
-      !
-      CALL iom_put( 'e2u'      , e2u (:,:) )
-      CALL iom_put( 'e1v'      , e1v (:,:) )
-      CALL iom_put( 'areacello', area(:,:) )
+         zarea_ssh(:,:) = e1e2t(:,:) * ssh(:,:,Kmm)
+      ENDIF
+      !
+      CALL iom_put( 'e2u'      , e2u  (:,:) )
+      CALL iom_put( 'e1v'      , e1v  (:,:) )
+      CALL iom_put( 'areacello', e1e2t(:,:) )
       !
       IF( iom_use( 'volcello' ) .OR. iom_use( 'masscello' )  ) THEN  
          zrhd(:,:,jpk) = 0._wp        ! ocean volume ; rhd is used as workspace
          DO jk = 1, jpkm1
-            zrhd(:,:,jk) = area(:,:) * e3t(:,:,jk,Kmm) * tmask(:,:,jk)
+            zrhd(:,:,jk) = e1e2t(:,:) * e3t(:,:,jk,Kmm) * tmask(:,:,jk)
          END DO
+         DO jk = 1, jpk
+            z3d(:,:,jk) =  rho0 * e3t(:,:,jk,Kmm) * tmask(:,:,jk)
+         END DO 
          CALL iom_put( 'volcello'  , zrhd(:,:,:)  )  ! WARNING not consistent with CMIP DR where volcello is at ca. 2000
-         CALL iom_put( 'masscello' , rho0 * e3t(:,:,:,Kmm) * tmask(:,:,:) )  ! ocean mass
+         CALL iom_put( 'masscello' , z3d (:,:,:) )   ! ocean mass
       ENDIF 
       !
@@ -129,5 +132,8 @@
          ztsn(:,:,:,jp_tem) = ts(:,:,:,jp_tem,Kmm)                    ! thermosteric ssh
          ztsn(:,:,:,jp_sal) = sn0(:,:,:)
-         CALL eos( ztsn, zrhd, gdept(:,:,:,Kmm) )                       ! now in situ density using initial salinity
+         DO jk = 1, jpk
+            zgdept(:,:,jk) = gdept(:,:,jk,Kmm)
+         END DO
+         CALL eos( ztsn, zrhd, zgdept)                       ! now in situ density using initial salinity
          !
          zbotpres(:,:) = 0._wp                        ! no atmospheric surface pressure, levitating sea-ice
@@ -151,10 +157,10 @@
          END IF
          !                                         
-         zarho = glob_sum( 'diaar5', area(:,:) * zbotpres(:,:) ) 
+         zarho = glob_sum( 'diaar5', e1e2t(:,:) * zbotpres(:,:) ) 
          zssh_steric = - zarho / area_tot
          CALL iom_put( 'sshthster', zssh_steric )
       
          !                                         ! steric sea surface height
-         CALL eos( ts(:,:,:,:,Kmm), zrhd, zrhop, gdept(:,:,:,Kmm) )                 ! now in situ and potential density
+         CALL eos( ts(:,:,:,:,Kmm), zrhd, zrhop, zgdept )                 ! now in situ and potential density
          zrhop(:,:,jpk) = 0._wp
          CALL iom_put( 'rhop', zrhop )
@@ -177,5 +183,5 @@
          END IF
          !    
-         zarho = glob_sum( 'diaar5', area(:,:) * zbotpres(:,:) ) 
+         zarho = glob_sum( 'diaar5', e1e2t(:,:) * zbotpres(:,:) ) 
          zssh_steric = - zarho / area_tot
          CALL iom_put( 'sshsteric', zssh_steric )
@@ -191,5 +197,5 @@
           ztsn(:,:,:,:) = 0._wp                    ! ztsn(:,:,1,jp_tem/sal) is used here as 2D Workspace for temperature & salinity
           DO_3D_11_11( 1, jpkm1 )
-             zztmp = area(ji,jj) * e3t(ji,jj,jk,Kmm)
+             zztmp = e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm)
              ztsn(ji,jj,1,jp_tem) = ztsn(ji,jj,1,jp_tem) + zztmp * ts(ji,jj,jk,jp_tem,Kmm)
              ztsn(ji,jj,1,jp_sal) = ztsn(ji,jj,1,jp_sal) + zztmp * ts(ji,jj,jk,jp_sal,Kmm)
@@ -237,5 +243,5 @@
                z2d(:,:) = 0._wp
                DO jk = 1, jpkm1
-                 z2d(:,:) = z2d(:,:) + area(:,:) * e3t(:,:,jk,Kmm) * ztpot(:,:,jk)
+                 z2d(:,:) = z2d(:,:) + e1e2t(:,:) * e3t(:,:,jk,Kmm) * ztpot(:,:,jk)
                END DO
                ztemp = glob_sum( 'diaar5', z2d(:,:)  ) 
@@ -244,5 +250,5 @@
              !
              IF( iom_use( 'ssttot' ) ) THEN   ! Output potential temperature in case we use TEOS-10
-               zsst = glob_sum( 'diaar5',  area(:,:) * ztpot(:,:,1)  ) 
+               zsst = glob_sum( 'diaar5',  e1e2t(:,:) * ztpot(:,:,1)  ) 
                CALL iom_put( 'ssttot', zsst / area_tot )
              ENDIF
@@ -259,5 +265,5 @@
       ELSE       
          IF( iom_use('ssttot') ) THEN   ! Output sst in case we use EOS-80
-            zsst  = glob_sum( 'diaar5', area(:,:) * ts(:,:,1,jp_tem,Kmm) )
+            zsst  = glob_sum( 'diaar5', e1e2t(:,:) * ts(:,:,1,jp_tem,Kmm) )
             CALL iom_put('ssttot', zsst / area_tot )
          ENDIF
@@ -375,6 +381,5 @@
          IF( dia_ar5_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dia_ar5_init : unable to allocate arrays' )
 
-         area(:,:) = e1e2t(:,:)
-         area_tot  = glob_sum( 'diaar5', area(:,:) )
+         area_tot  = glob_sum( 'diaar5', e1e2t(:,:) )
 
          ALLOCATE( zvol0(jpi,jpj) )
@@ -383,5 +388,5 @@
          DO_3D_11_11( 1, jpkm1 )
             idep = tmask(ji,jj,jk) * e3t_0(ji,jj,jk)
-            zvol0 (ji,jj) = zvol0 (ji,jj) +  idep * area(ji,jj)
+            zvol0 (ji,jj) = zvol0 (ji,jj) +  idep * e1e2t(ji,jj)
             thick0(ji,jj) = thick0(ji,jj) +  idep    
          END_3D
