Index: NEMO/trunk/src/OCE/DIA/dia25h.F90
===================================================================
--- NEMO/trunk/src/OCE/DIA/dia25h.F90	(revision 15234)
+++ NEMO/trunk/src/OCE/DIA/dia25h.F90	(revision 15249)
@@ -32,4 +32,6 @@
    REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   en_25h  , rmxln_25h
 
+!! * Substitutions
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -51,4 +53,5 @@
       INTEGER ::   ios                 ! Local integer output status for namelist read
       INTEGER ::   ierror              ! Local integer for memory allocation
+      INTEGER ::   ji, jj, jk
       !
       NAMELIST/nam_dia25h/ ln_dia25h
@@ -73,12 +76,12 @@
       ! ------------------- !
       !                                ! ocean arrays
-      ALLOCATE( tn_25h (jpi,jpj,jpk), sn_25h (jpi,jpj,jpk), sshn_25h(jpi,jpj)  ,     &
-         &      un_25h (jpi,jpj,jpk), vn_25h (jpi,jpj,jpk), wn_25h(jpi,jpj,jpk),     &
-         &      avt_25h(jpi,jpj,jpk), avm_25h(jpi,jpj,jpk),                      STAT=ierror )
+      ALLOCATE( tn_25h (A2D(0),jpk), sn_25h (A2D(0),jpk), sshn_25h(A2D(0))  ,     &
+         &      un_25h (A2D(0),jpk), vn_25h (A2D(0),jpk), wn_25h(A2D(0),jpk),     &
+         &      avt_25h(A2D(0),jpk), avm_25h(A2D(0),jpk),                      STAT=ierror )
       IF( ierror > 0 ) THEN
          CALL ctl_stop( 'dia_25h: unable to allocate ocean arrays' )   ;   RETURN
       ENDIF
       IF( ln_zdftke ) THEN             ! TKE physics
-         ALLOCATE( en_25h(jpi,jpj,jpk), STAT=ierror )
+         ALLOCATE( en_25h(A2D(0),jpk), STAT=ierror )
          IF( ierror > 0 ) THEN
             CALL ctl_stop( 'dia_25h: unable to allocate en_25h' )   ;   RETURN
@@ -86,5 +89,5 @@
       ENDIF
       IF( ln_zdfgls ) THEN             ! GLS physics
-         ALLOCATE( en_25h(jpi,jpj,jpk), rmxln_25h(jpi,jpj,jpk), STAT=ierror )
+         ALLOCATE( en_25h(A2D(0),jpk), rmxln_25h(A2D(0),jpk), STAT=ierror )
          IF( ierror > 0 ) THEN
             CALL ctl_stop( 'dia_25h: unable to allocate en_25h and rmxln_25h' )   ;   RETURN
@@ -94,18 +97,26 @@
       ! 2 - Assign Initial Values !
       ! ------------------------- !
-      cnt_25h = 1  ! sets the first value of sum at timestep 1 (note - should strictly be at timestep zero so before values used where possible) 
-      tn_25h  (:,:,:) = ts (:,:,:,jp_tem,Kbb)
-      sn_25h  (:,:,:) = ts (:,:,:,jp_sal,Kbb)
-      sshn_25h(:,:)   = ssh(:,:,Kbb)
-      un_25h  (:,:,:) = uu  (:,:,:,Kbb)
-      vn_25h  (:,:,:) = vv  (:,:,:,Kbb)
-      avt_25h (:,:,:) = avt (:,:,:)
-      avm_25h (:,:,:) = avm (:,:,:)
+      cnt_25h = 1  ! sets the first value of sum at timestep 1 (note - should strictly be at timestep zero so before values used where possible)
+      DO_3D( 0, 0, 0, 0, 1, jpk )
+         tn_25h (ji,jj,jk) = ts (ji,jj,jk,jp_tem,Kbb)
+         sn_25h (ji,jj,jk) = ts (ji,jj,jk,jp_sal,Kbb)
+         un_25h (ji,jj,jk) = uu (ji,jj,jk,Kbb)
+         vn_25h (ji,jj,jk) = vv (ji,jj,jk,Kbb)
+         avt_25h(ji,jj,jk) = avt(ji,jj,jk)
+         avm_25h(ji,jj,jk) = avm(ji,jj,jk)
+      END_3D
+      DO_2D( 0, 0, 0, 0 )
+         sshn_25h(ji,jj) = ssh(ji,jj,Kbb)
+      END_2D
       IF( ln_zdftke ) THEN
-         en_25h(:,:,:) = en(:,:,:)
+         DO_3D( 0, 0, 0, 0, 1, jpk )
+            en_25h(ji,jj,jk) = en(ji,jj,jk)
+         END_3D
       ENDIF
       IF( ln_zdfgls ) THEN
-         en_25h   (:,:,:) = en    (:,:,:)
-         rmxln_25h(:,:,:) = hmxl_n(:,:,:)
+         DO_3D( 0, 0, 0, 0, 1, jpk )
+            en_25h   (ji,jj,jk) = en    (ji,jj,jk)
+            rmxln_25h(ji,jj,jk) = hmxl_n(ji,jj,jk)
+         END_3D
       ENDIF
 #if defined key_si3
@@ -132,7 +143,7 @@
       REAL(wp)                         ::   zsto, zout, zmax, zjulian, zmdi   ! local scalars
       INTEGER                          ::   i_steps                           ! no of timesteps per hour
-      REAL(wp), DIMENSION(jpi,jpj    ) ::   zw2d, un_dm, vn_dm                ! workspace
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zw3d                              ! workspace
-      REAL(wp), DIMENSION(jpi,jpj,3)   ::   zwtmb                             ! workspace
+      REAL(wp), DIMENSION(A2D(0)    )  ::   zw2d, un_dm, vn_dm                ! workspace
+      REAL(wp), DIMENSION(A2D(0),jpk)  ::   zw3d                              ! workspace
+      REAL(wp), DIMENSION(A2D(0),3)    ::   zwtmb                             ! workspace
       !!----------------------------------------------------------------------
 
@@ -151,5 +162,7 @@
       ! wn_25h could not be initialised in dia_25h_init, so we do it here instead
       IF( kt == nn_it000 ) THEN
-         wn_25h(:,:,:) = ww(:,:,:)
+         DO_3D( 0, 0, 0, 0, 1, jpk )
+            wn_25h(ji,jj,jk) = ww(ji,jj,jk)
+         END_3D
       ENDIF
 
@@ -162,18 +175,26 @@
          ENDIF
 
-         tn_25h  (:,:,:)     = tn_25h  (:,:,:) + ts (:,:,:,jp_tem,Kmm)
-         sn_25h  (:,:,:)     = sn_25h  (:,:,:) + ts (:,:,:,jp_sal,Kmm)
-         sshn_25h(:,:)       = sshn_25h(:,:)   + ssh(:,:,Kmm)
-         un_25h  (:,:,:)     = un_25h  (:,:,:) + uu  (:,:,:,Kmm)
-         vn_25h  (:,:,:)     = vn_25h  (:,:,:) + vv  (:,:,:,Kmm)
-         wn_25h  (:,:,:)     = wn_25h  (:,:,:) + ww  (:,:,:)
-         avt_25h (:,:,:)     = avt_25h (:,:,:) + avt (:,:,:)
-         avm_25h (:,:,:)     = avm_25h (:,:,:) + avm (:,:,:)
+         DO_3D( 0, 0, 0, 0, 1, jpk )
+            tn_25h  (ji,jj,jk) = tn_25h  (ji,jj,jk) + ts (ji,jj,jk,jp_tem,Kmm)
+            sn_25h  (ji,jj,jk) = sn_25h  (ji,jj,jk) + ts (ji,jj,jk,jp_sal,Kmm)
+            un_25h  (ji,jj,jk) = un_25h  (ji,jj,jk) + uu (ji,jj,jk,Kmm)
+            vn_25h  (ji,jj,jk) = vn_25h  (ji,jj,jk) + vv (ji,jj,jk,Kmm)
+            wn_25h  (ji,jj,jk) = wn_25h  (ji,jj,jk) + ww (ji,jj,jk)
+            avt_25h (ji,jj,jk) = avt_25h (ji,jj,jk) + avt(ji,jj,jk)
+            avm_25h (ji,jj,jk) = avm_25h (ji,jj,jk) + avm(ji,jj,jk)
+         END_3D
+         DO_2D( 0, 0, 0, 0 )
+            sshn_25h(ji,jj)    = sshn_25h(ji,jj)    + ssh(ji,jj,Kmm)
+         END_2D
          IF( ln_zdftke ) THEN
-            en_25h(:,:,:)    = en_25h  (:,:,:) + en(:,:,:)
+            DO_3D( 0, 0, 0, 0, 1, jpk )
+               en_25h(ji,jj,jk) = en_25h(ji,jj,jk) + en(ji,jj,jk)
+            END_3D
          ENDIF
          IF( ln_zdfgls ) THEN
-            en_25h   (:,:,:) = en_25h   (:,:,:) + en    (:,:,:)
-            rmxln_25h(:,:,:) = rmxln_25h(:,:,:) + hmxl_n(:,:,:)
+            DO_3D( 0, 0, 0, 0, 1, jpk )
+               en_25h   (ji,jj,jk) = en_25h   (ji,jj,jk) + en    (ji,jj,jk)
+               rmxln_25h(ji,jj,jk) = rmxln_25h(ji,jj,jk) + hmxl_n(ji,jj,jk)
+            END_3D
          ENDIF
          cnt_25h = cnt_25h + 1
@@ -212,9 +233,15 @@
          zmdi=1.e+20 !missing data indicator for masking
          ! write tracers (instantaneous)
-         zw3d(:,:,:) = tn_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
+         DO_3D( 0, 0, 0, 0, 1, jpk )
+            zw3d(ji,jj,jk) = tn_25h(ji,jj,jk)*tmask(ji,jj,jk) + zmdi*(1.0-tmask(ji,jj,jk))
+         END_3D
          CALL iom_put("temper25h", zw3d)   ! potential temperature
-         zw3d(:,:,:) = sn_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
+         DO_3D( 0, 0, 0, 0, 1, jpk )
+            zw3d(ji,jj,jk) = sn_25h(ji,jj,jk)*tmask(ji,jj,jk) + zmdi*(1.0-tmask(ji,jj,jk))
+         END_3D
          CALL iom_put( "salin25h", zw3d  )   ! salinity
-         zw2d(:,:) = sshn_25h(:,:)*tmask(:,:,1) + zmdi*(1.0-tmask(:,:,1))
+         DO_2D( 0, 0, 0, 0 )
+            zw2d(ji,jj) = sshn_25h(ji,jj)*tmask(ji,jj,1) + zmdi*(1.0-tmask(ji,jj,1))
+         END_2D
          IF( ll_wd ) THEN
             CALL iom_put( "ssh25h", zw2d+ssh_ref )   ! sea surface 
@@ -223,41 +250,65 @@
          ENDIF
          ! Write velocities (instantaneous)
-         zw3d(:,:,:) = un_25h(:,:,:)*umask(:,:,:) + zmdi*(1.0-umask(:,:,:))
+         DO_3D( 0, 0, 0, 0, 1, jpk )
+            zw3d(ji,jj,jk) = un_25h(ji,jj,jk)*umask(ji,jj,jk) + zmdi*(1.0-umask(ji,jj,jk))
+         END_3D
          CALL iom_put("vozocrtx25h", zw3d)    ! i-current
-         zw3d(:,:,:) = vn_25h(:,:,:)*vmask(:,:,:) + zmdi*(1.0-vmask(:,:,:))
+         DO_3D( 0, 0, 0, 0, 1, jpk )
+            zw3d(ji,jj,jk) = vn_25h(ji,jj,jk)*vmask(ji,jj,jk) + zmdi*(1.0-vmask(ji,jj,jk))
+         END_3D
          CALL iom_put("vomecrty25h", zw3d  )   ! j-current
-         zw3d(:,:,:) = wn_25h(:,:,:)*wmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
+         DO_3D( 0, 0, 0, 0, 1, jpk )
+            zw3d(ji,jj,jk) = wn_25h(ji,jj,jk)*wmask(ji,jj,jk) + zmdi*(1.0-tmask(ji,jj,jk))
+         END_3D
          CALL iom_put("vovecrtz25h", zw3d )   ! k-current
          ! Write vertical physics
-         zw3d(:,:,:) = avt_25h(:,:,:)*wmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
+         DO_3D( 0, 0, 0, 0, 1, jpk )
+            zw3d(ji,jj,jk) = avt_25h(ji,jj,jk)*wmask(ji,jj,jk) + zmdi*(1.0-tmask(ji,jj,jk))
+         END_3D
          CALL iom_put("avt25h", zw3d )   ! diffusivity
-         zw3d(:,:,:) = avm_25h(:,:,:)*wmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
+         DO_3D( 0, 0, 0, 0, 1, jpk )
+            zw3d(ji,jj,jk) = avm_25h(ji,jj,jk)*wmask(ji,jj,jk) + zmdi*(1.0-tmask(ji,jj,jk))
+         END_3D
          CALL iom_put("avm25h", zw3d)   ! viscosity
          IF( ln_zdftke ) THEN
-            zw3d(:,:,:) = en_25h(:,:,:)*wmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
+            DO_3D( 0, 0, 0, 0, 1, jpk )
+               zw3d(ji,jj,jk) = en_25h(ji,jj,jk)*wmask(ji,jj,jk) + zmdi*(1.0-tmask(ji,jj,jk))
+            END_3D
             CALL iom_put("tke25h", zw3d)   ! tke
          ENDIF
          IF( ln_zdfgls ) THEN
-            zw3d(:,:,:) = en_25h(:,:,:)*wmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
+            DO_3D( 0, 0, 0, 0, 1, jpk )
+               zw3d(ji,jj,jk) = en_25h(ji,jj,jk)*wmask(ji,jj,jk) + zmdi*(1.0-tmask(ji,jj,jk))
+            END_3D
             CALL iom_put("tke25h", zw3d)   ! tke
-            zw3d(:,:,:) = rmxln_25h(:,:,:)*wmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
+            DO_3D( 0, 0, 0, 0, 1, jpk )
+               zw3d(ji,jj,jk) = rmxln_25h(ji,jj,jk)*wmask(ji,jj,jk) + zmdi*(1.0-tmask(ji,jj,jk))
+            END_3D
             CALL iom_put( "mxln25h",zw3d)
          ENDIF
          !
          ! After the write reset the values to cnt=1 and sum values equal current value 
-         tn_25h  (:,:,:) = ts (:,:,:,jp_tem,Kmm)
-         sn_25h  (:,:,:) = ts (:,:,:,jp_sal,Kmm)
-         sshn_25h(:,:)   = ssh(:,:,Kmm)
-         un_25h  (:,:,:) = uu  (:,:,:,Kmm)
-         vn_25h  (:,:,:) = vv  (:,:,:,Kmm)
-         wn_25h  (:,:,:) = ww  (:,:,:)
-         avt_25h (:,:,:) = avt (:,:,:)
-         avm_25h (:,:,:) = avm (:,:,:)
+         DO_3D( 0, 0, 0, 0, 1, jpk )
+            tn_25h  (ji,jj,jk) = ts (ji,jj,jk,jp_tem,Kmm)
+            sn_25h  (ji,jj,jk) = ts (ji,jj,jk,jp_sal,Kmm)
+            un_25h  (ji,jj,jk) = uu (ji,jj,jk,Kmm)
+            vn_25h  (ji,jj,jk) = vv (ji,jj,jk,Kmm)
+            wn_25h  (ji,jj,jk) = ww (ji,jj,jk)
+            avt_25h (ji,jj,jk) = avt(ji,jj,jk)
+            avm_25h (ji,jj,jk) = avm(ji,jj,jk)
+         END_3D
+         DO_2D( 0, 0, 0, 0 )
+            sshn_25h(ji,jj)    = ssh(ji,jj,Kmm)
+         END_2D
          IF( ln_zdftke ) THEN
-            en_25h(:,:,:) = en(:,:,:)
+            DO_3D( 0, 0, 0, 0, 1, jpk )
+               en_25h(ji,jj,jk) = en(ji,jj,jk)
+            END_3D
          ENDIF
          IF( ln_zdfgls ) THEN
-            en_25h   (:,:,:) = en    (:,:,:)
-            rmxln_25h(:,:,:) = hmxl_n(:,:,:)
+            DO_3D( 0, 0, 0, 0, 1, jpk )
+               en_25h   (ji,jj,jk) = en    (ji,jj,jk)
+               rmxln_25h(ji,jj,jk) = hmxl_n(ji,jj,jk)
+            END_3D
          ENDIF
          cnt_25h = 1
