Index: NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/BENCH/MY_SRC/diawri.F90
===================================================================
--- NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/BENCH/MY_SRC/diawri.F90	(revision 12486)
+++ 	(revision )
@@ -1,504 +1,0 @@
-MODULE diawri
-   !!======================================================================
-   !!                     ***  MODULE  diawri  ***
-   !! Ocean diagnostics :  write ocean output files
-   !!=====================================================================
-   !! History :  OPA  ! 1991-03  (M.-A. Foujols)  Original code
-   !!            4.0  ! 1991-11  (G. Madec)
-   !!                 ! 1992-06  (M. Imbard)  correction restart file
-   !!                 ! 1992-07  (M. Imbard)  split into diawri and rstwri
-   !!                 ! 1993-03  (M. Imbard)  suppress writibm
-   !!                 ! 1998-01  (C. Levy)  NETCDF format using ioipsl INTERFACE
-   !!                 ! 1999-02  (E. Guilyardi)  name of netCDF files + variables
-   !!            8.2  ! 2000-06  (M. Imbard)  Original code (diabort.F)
-   !!   NEMO     1.0  ! 2002-06  (A.Bozec, E. Durand)  Original code (diainit.F)
-   !!             -   ! 2002-09  (G. Madec)  F90: Free form and module
-   !!             -   ! 2002-12  (G. Madec)  merge of diabort and diainit, F90
-   !!                 ! 2005-11  (V. Garnier) Surface pressure gradient organization
-   !!            3.2  ! 2008-11  (B. Lemaire) creation from old diawri
-   !!            3.7  ! 2014-01  (G. Madec) remove eddy induced velocity from no-IOM output
-   !!                 !                     change name of output variables in dia_wri_state
-   !!----------------------------------------------------------------------
-
-   !!----------------------------------------------------------------------
-   !!   dia_wri       : create the standart output files
-   !!   dia_wri_state : create an output NetCDF file for a single instantaeous ocean state and forcing fields
-   !!----------------------------------------------------------------------
-   USE oce            ! ocean dynamics and tracers 
-   USE dom_oce        ! ocean space and time domain
-   USE phycst         ! physical constants
-   USE dianam         ! build name of file (routine)
-   USE diahth         ! thermocline diagnostics
-   USE dynadv   , ONLY: ln_dynadv_vec
-   USE icb_oce        ! Icebergs
-   USE icbdia         ! Iceberg budgets
-   USE ldftra         ! lateral physics: eddy diffusivity coef.
-   USE ldfdyn         ! lateral physics: eddy viscosity   coef.
-   USE sbc_oce        ! Surface boundary condition: ocean fields
-   USE sbc_ice        ! Surface boundary condition: ice fields
-   USE sbcssr         ! restoring term toward SST/SSS climatology
-   USE sbcwave        ! wave parameters
-   USE wet_dry        ! wetting and drying
-   USE zdf_oce        ! ocean vertical physics
-   USE zdfdrg         ! ocean vertical physics: top/bottom friction
-   USE zdfmxl         ! mixed layer
-   !
-   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
-   USE in_out_manager ! I/O manager
-   USE diatmb         ! Top,middle,bottom output
-   USE dia25h         ! 25h Mean output
-   USE iom            ! 
-   USE ioipsl         ! 
-
-#if defined key_si3
-   USE ice 
-   USE icewri 
-#endif
-   USE lib_mpp         ! MPP library
-   USE timing          ! preformance summary
-   USE diu_bulk        ! diurnal warm layer
-   USE diu_coolskin    ! Cool skin
-
-   IMPLICIT NONE
-   PRIVATE
-
-   PUBLIC   dia_wri                 ! routines called by step.F90
-   PUBLIC   dia_wri_state
-   PUBLIC   dia_wri_alloc           ! Called by nemogcm module
-
-   INTEGER ::   nid_T, nz_T, nh_T, ndim_T, ndim_hT   ! grid_T file
-   INTEGER ::          nb_T              , ndim_bT   ! grid_T file
-   INTEGER ::   nid_U, nz_U, nh_U, ndim_U, ndim_hU   ! grid_U file
-   INTEGER ::   nid_V, nz_V, nh_V, ndim_V, ndim_hV   ! grid_V file
-   INTEGER ::   nid_W, nz_W, nh_W                    ! grid_W file
-   INTEGER ::   ndex(1)                              ! ???
-   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hT, ndex_hU, ndex_hV
-   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_T, ndex_U, ndex_V
-   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_bT
-
-   !!----------------------------------------------------------------------
-   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
-   !! $Id: diawri.F90 9598 2018-05-15 22:47:16Z nicolasmartin $
-   !! Software governed by the CeCILL licence     (./LICENSE)
-   !!----------------------------------------------------------------------
-CONTAINS
-
-#if defined key_iomput
-   !!----------------------------------------------------------------------
-   !!   'key_iomput'                                        use IOM library
-   !!----------------------------------------------------------------------
-   INTEGER FUNCTION dia_wri_alloc()
-      !
-      dia_wri_alloc = 0
-      !
-   END FUNCTION dia_wri_alloc
-
-   
-   SUBROUTINE dia_wri( kt, Kmm )
-      !!---------------------------------------------------------------------
-      !!                  ***  ROUTINE dia_wri  ***
-      !!                   
-      !! ** Purpose :   Standard output of opa: dynamics and tracer fields 
-      !!      NETCDF format is used by default 
-      !!
-      !! ** Method  :  use iom_put
-      !!----------------------------------------------------------------------
-      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
-      INTEGER, INTENT( in ) ::   Kmm     ! ocean time level index
-      !!
-      INTEGER ::   ji, jj, jk       ! dummy loop indices
-      INTEGER ::   ikbot            ! local integer
-      REAL(wp)::   zztmp , zztmpx   ! local scalar
-      REAL(wp)::   zztmp2, zztmpy   !   -      -
-      REAL(wp), DIMENSION(jpi,jpj)     ::   z2d   ! 2D workspace
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   z3d   ! 3D workspace
-      !!----------------------------------------------------------------------
-      ! 
-      IF( ln_timing )   CALL timing_start('dia_wri')
-      ! 
-      ! Output the initial state and forcings
-      IF( ninist == 1 ) THEN                       
-         CALL dia_wri_state( Kmm, 'output.init' )
-         ninist = 0
-      ENDIF
-
-      ! Output of initial vertical scale factor
-      CALL iom_put("e3t_0", e3t_0(:,:,:) )
-      CALL iom_put("e3u_0", e3u_0(:,:,:) )
-      CALL iom_put("e3v_0", e3v_0(:,:,:) )
-      !
-      CALL iom_put( "e3t" , e3t(:,:,:,Kmm) )
-      CALL iom_put( "e3u" , e3u(:,:,:,Kmm) )
-      CALL iom_put( "e3v" , e3v(:,:,:,Kmm) )
-      CALL iom_put( "e3w" , e3w(:,:,:,Kmm) )
-      IF( iom_use("e3tdef") )   &
-         CALL iom_put( "e3tdef"  , ( ( e3t(:,:,:,Kmm) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 )
-
-      IF( ll_wd ) THEN
-         CALL iom_put( "ssh" , (ssh(:,:,Kmm)+ssh_ref)*tmask(:,:,1) )   ! sea surface height (brought back to the reference used for wetting and drying)
-      ELSE
-         CALL iom_put( "ssh" , ssh(:,:,Kmm) )              ! sea surface height
-      ENDIF
-
-      IF( iom_use("wetdep") )   &                  ! wet depth
-         CALL iom_put( "wetdep" , ht_0(:,:) + ssh(:,:,Kmm) )
-      
-      CALL iom_put( "toce", ts(:,:,:,jp_tem,Kmm) )    ! 3D temperature
-      CALL iom_put(  "sst", ts(:,:,1,jp_tem,Kmm) )    ! surface temperature
-      IF ( iom_use("sbt") ) THEN
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               ikbot = mbkt(ji,jj)
-               z2d(ji,jj) = ts(ji,jj,ikbot,jp_tem,Kmm)
-            END DO
-         END DO
-         CALL iom_put( "sbt", z2d )                ! bottom temperature
-      ENDIF
-      
-      CALL iom_put( "soce", ts(:,:,:,jp_sal,Kmm) )    ! 3D salinity
-      CALL iom_put(  "sss", ts(:,:,1,jp_sal,Kmm) )    ! surface salinity
-      IF ( iom_use("sbs") ) THEN
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               ikbot = mbkt(ji,jj)
-               z2d(ji,jj) = ts(ji,jj,ikbot,jp_sal,Kmm)
-            END DO
-         END DO
-         CALL iom_put( "sbs", z2d )                ! bottom salinity
-      ENDIF
-
-      IF ( iom_use("taubot") ) THEN                ! bottom stress
-         zztmp = rho0 * 0.25
-         z2d(:,:) = 0._wp
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               zztmp2 = (  ( rCdU_bot(ji+1,jj)+rCdU_bot(ji  ,jj) ) * uu(ji  ,jj,mbku(ji  ,jj),Kmm)  )**2   &
-                  &   + (  ( rCdU_bot(ji  ,jj)+rCdU_bot(ji-1,jj) ) * uu(ji-1,jj,mbku(ji-1,jj),Kmm)  )**2   &
-                  &   + (  ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj  ) ) * vv(ji,jj  ,mbkv(ji,jj  ),Kmm)  )**2   &
-                  &   + (  ( rCdU_bot(ji,jj  )+rCdU_bot(ji,jj-1) ) * vv(ji,jj-1,mbkv(ji,jj-1),Kmm)  )**2
-               z2d(ji,jj) = zztmp * SQRT( zztmp2 ) * tmask(ji,jj,1) 
-               !
-            END DO
-         END DO
-         CALL lbc_lnk( 'diawri', z2d, 'T', 1. )
-         CALL iom_put( "taubot", z2d )           
-      ENDIF
-         
-      CALL iom_put( "uoce", uu(:,:,:,Kmm) )            ! 3D i-current
-      CALL iom_put(  "ssu", uu(:,:,1,Kmm) )            ! surface i-current
-      IF ( iom_use("sbu") ) THEN
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               ikbot = mbku(ji,jj)
-               z2d(ji,jj) = uu(ji,jj,ikbot,Kmm)
-            END DO
-         END DO
-         CALL iom_put( "sbu", z2d )                ! bottom i-current
-      ENDIF
-      
-      CALL iom_put( "voce", vv(:,:,:,Kmm) )            ! 3D j-current
-      CALL iom_put(  "ssv", vv(:,:,1,Kmm) )            ! surface j-current
-      IF ( iom_use("sbv") ) THEN
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               ikbot = mbkv(ji,jj)
-               z2d(ji,jj) = vv(ji,jj,ikbot,Kmm)
-            END DO
-         END DO
-         CALL iom_put( "sbv", z2d )                ! bottom j-current
-      ENDIF
-
-      IF( ln_zad_Aimp ) ww = ww + wi               ! Recombine explicit and implicit parts of vertical velocity for diagnostic output
-      !
-      CALL iom_put( "woce", ww )                   ! vertical velocity
-      IF( iom_use('w_masstr') .OR. iom_use('w_masstr2') ) THEN   ! vertical mass transport & its square value
-         ! Caution: in the VVL case, it only correponds to the baroclinic mass transport.
-         z2d(:,:) = rho0 * e1e2t(:,:)
-         DO jk = 1, jpk
-            z3d(:,:,jk) = ww(:,:,jk) * z2d(:,:)
-         END DO
-         CALL iom_put( "w_masstr" , z3d )  
-         IF( iom_use('w_masstr2') )   CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) )
-      ENDIF
-      !
-      IF( ln_zad_Aimp ) ww = ww - wi               ! Remove implicit part of vertical velocity that was added for diagnostic output
-
-      CALL iom_put( "avt" , avt )                  ! T vert. eddy diff. coef.
-      CALL iom_put( "avs" , avs )                  ! S vert. eddy diff. coef.
-      CALL iom_put( "avm" , avm )                  ! T vert. eddy visc. coef.
-
-      IF( iom_use('logavt') )   CALL iom_put( "logavt", LOG( MAX( 1.e-20_wp, avt(:,:,:) ) ) )
-      IF( iom_use('logavs') )   CALL iom_put( "logavs", LOG( MAX( 1.e-20_wp, avs(:,:,:) ) ) )
-
-      IF ( iom_use("sstgrad") .OR. iom_use("sstgrad2") ) THEN
-         DO jj = 2, jpjm1                                    ! sst gradient
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               zztmp  = ts(ji,jj,1,jp_tem,Kmm)
-               zztmpx = ( ts(ji+1,jj,1,jp_tem,Kmm) - zztmp ) * r1_e1u(ji,jj) + ( zztmp - ts(ji-1,jj  ,1,jp_tem,Kmm) ) * r1_e1u(ji-1,jj)
-               zztmpy = ( ts(ji,jj+1,1,jp_tem,Kmm) - zztmp ) * r1_e2v(ji,jj) + ( zztmp - ts(ji  ,jj-1,1,jp_tem,Kmm) ) * r1_e2v(ji,jj-1)
-               z2d(ji,jj) = 0.25 * ( zztmpx * zztmpx + zztmpy * zztmpy )   &
-                  &              * umask(ji,jj,1) * umask(ji-1,jj,1) * vmask(ji,jj,1) * umask(ji,jj-1,1)
-            END DO
-         END DO
-         CALL lbc_lnk( 'diawri', z2d, 'T', 1. )
-         CALL iom_put( "sstgrad2",  z2d )          ! square of module of sst gradient
-         z2d(:,:) = SQRT( z2d(:,:) )
-         CALL iom_put( "sstgrad" ,  z2d )          ! module of sst gradient
-      ENDIF
-         
-      ! heat and salt contents
-      IF( iom_use("heatc") ) THEN
-         z2d(:,:)  = 0._wp 
-         DO jk = 1, jpkm1
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_tem,Kmm) * tmask(ji,jj,jk)
-               END DO
-            END DO
-         END DO
-         CALL iom_put( "heatc", rho0_rcp * z2d )   ! vertically integrated heat content (J/m2)
-      ENDIF
-
-      IF( iom_use("saltc") ) THEN
-         z2d(:,:)  = 0._wp 
-         DO jk = 1, jpkm1
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) * tmask(ji,jj,jk)
-               END DO
-            END DO
-         END DO
-         CALL iom_put( "saltc", rho0 * z2d )          ! vertically integrated salt content (PSU*kg/m2)
-      ENDIF
-      !
-      IF ( iom_use("eken") ) THEN
-         z3d(:,:,jpk) = 0._wp 
-         DO jk = 1, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zztmp  = 0.25_wp * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
-                  z3d(ji,jj,jk) = zztmp * (  uu(ji-1,jj,jk,Kmm)**2 * e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm)   &
-                     &                     + uu(ji  ,jj,jk,Kmm)**2 * e2u(ji  ,jj) * e3u(ji  ,jj,jk,Kmm)   &
-                     &                     + vv(ji,jj-1,jk,Kmm)**2 * e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm)   &
-                     &                     + vv(ji,jj  ,jk,Kmm)**2 * e1v(ji,jj  ) * e3v(ji,jj  ,jk,Kmm)   )
-               END DO
-            END DO
-         END DO
-         CALL lbc_lnk( 'diawri', z3d, 'T', 1. )
-         CALL iom_put( "eken", z3d )                 ! kinetic energy
-      ENDIF
-      !
-      CALL iom_put( "hdiv", hdiv )                  ! Horizontal divergence
-      !
-      IF( iom_use("u_masstr") .OR. iom_use("u_masstr_vint") .OR. iom_use("u_heattr") .OR. iom_use("u_salttr") ) THEN
-         z3d(:,:,jpk) = 0.e0
-         z2d(:,:) = 0.e0
-         DO jk = 1, jpkm1
-            z3d(:,:,jk) = rho0 * uu(:,:,jk,Kmm) * e2u(:,:) * e3u(:,:,jk,Kmm) * umask(:,:,jk)
-            z2d(:,:) = z2d(:,:) + z3d(:,:,jk)
-         END DO
-         CALL iom_put( "u_masstr"     , z3d )         ! mass transport in i-direction
-         CALL iom_put( "u_masstr_vint", z2d )         ! mass transport in i-direction vertical sum
-      ENDIF
-      
-      IF( iom_use("u_heattr") ) THEN
-         z2d(:,:) = 0._wp 
-         DO jk = 1, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( ts(ji,jj,jk,jp_tem,Kmm) + ts(ji+1,jj,jk,jp_tem,Kmm) )
-               END DO
-            END DO
-         END DO
-         CALL lbc_lnk( 'diawri', z2d, 'U', -1. )
-         CALL iom_put( "u_heattr", 0.5*rcp * z2d )    ! heat transport in i-direction
-      ENDIF
-
-      IF( iom_use("u_salttr") ) THEN
-         z2d(:,:) = 0.e0 
-         DO jk = 1, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( ts(ji,jj,jk,jp_sal,Kmm) + ts(ji+1,jj,jk,jp_sal,Kmm) )
-               END DO
-            END DO
-         END DO
-         CALL lbc_lnk( 'diawri', z2d, 'U', -1. )
-         CALL iom_put( "u_salttr", 0.5 * z2d )        ! heat transport in i-direction
-      ENDIF
-
-      
-      IF( iom_use("v_masstr") .OR. iom_use("v_heattr") .OR. iom_use("v_salttr") ) THEN
-         z3d(:,:,jpk) = 0.e0
-         DO jk = 1, jpkm1
-            z3d(:,:,jk) = rho0 * vv(:,:,jk,Kmm) * e1v(:,:) * e3v(:,:,jk,Kmm) * vmask(:,:,jk)
-         END DO
-         CALL iom_put( "v_masstr", z3d )              ! mass transport in j-direction
-      ENDIF
-      
-      IF( iom_use("v_heattr") ) THEN
-         z2d(:,:) = 0.e0 
-         DO jk = 1, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( ts(ji,jj,jk,jp_tem,Kmm) + ts(ji,jj+1,jk,jp_tem,Kmm) )
-               END DO
-            END DO
-         END DO
-         CALL lbc_lnk( 'diawri', z2d, 'V', -1. )
-         CALL iom_put( "v_heattr", 0.5*rcp * z2d )    !  heat transport in j-direction
-      ENDIF
-
-      IF( iom_use("v_salttr") ) THEN
-         z2d(:,:) = 0._wp 
-         DO jk = 1, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( ts(ji,jj,jk,jp_sal,Kmm) + ts(ji,jj+1,jk,jp_sal,Kmm) )
-               END DO
-            END DO
-         END DO
-         CALL lbc_lnk( 'diawri', z2d, 'V', -1. )
-         CALL iom_put( "v_salttr", 0.5 * z2d )        !  heat transport in j-direction
-      ENDIF
-
-      IF( iom_use("tosmint") ) THEN
-         z2d(:,:) = 0._wp
-         DO jk = 1, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) *  ts(ji,jj,jk,jp_tem,Kmm)
-               END DO
-            END DO
-         END DO
-         CALL lbc_lnk( 'diawri', z2d, 'T', -1. )
-         CALL iom_put( "tosmint", rho0 * z2d )        ! Vertical integral of temperature
-      ENDIF
-      IF( iom_use("somint") ) THEN
-         z2d(:,:)=0._wp
-         DO jk = 1, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm)
-               END DO
-            END DO
-         END DO
-         CALL lbc_lnk( 'diawri', z2d, 'T', -1. )
-         CALL iom_put( "somint", rho0 * z2d )         ! Vertical integral of salinity
-      ENDIF
-
-      CALL iom_put( "bn2", rn2 )                      ! Brunt-Vaisala buoyancy frequency (N^2)
-      !
-
-      IF (ln_diatmb)   CALL dia_tmb( Kmm )            ! tmb values 
-          
-      IF (ln_dia25h)   CALL dia_25h( kt, Kmm )        ! 25h averaging
-
-      IF( ln_timing )   CALL timing_stop('dia_wri')
-      !
-   END SUBROUTINE dia_wri
-
-#else
-   !!----------------------------------------------------------------------
-   !!   Default option                                  use IOIPSL  library
-   !!----------------------------------------------------------------------
-
-   INTEGER FUNCTION dia_wri_alloc()
-      !
-      dia_wri_alloc = 0
-      !
-   END FUNCTION dia_wri_alloc
-
-   
-   SUBROUTINE dia_wri( kt, Kmm )
-
-      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
-      INTEGER, INTENT( in ) ::   Kmm  ! ocean time level index
-      IF( ninist == 1 ) THEN     !==  Output the initial state and forcings  ==!
-         CALL dia_wri_state( Kmm, 'output.init' )
-         ninist = 0
-      ENDIF
-      !
-      ! 0. Initialisation
-      ! -----------------
-
-   END SUBROUTINE dia_wri
-
-#endif
-
-   SUBROUTINE dia_wri_state( Kmm, cdfile_name )
-      !!---------------------------------------------------------------------
-      !!                 ***  ROUTINE dia_wri_state  ***
-      !!        
-      !! ** Purpose :   create a NetCDF file named cdfile_name which contains 
-      !!      the instantaneous ocean state and forcing fields.
-      !!        Used to find errors in the initial state or save the last
-      !!      ocean state in case of abnormal end of a simulation
-      !!
-      !! ** Method  :   NetCDF files using ioipsl
-      !!      File 'output.init.nc'  is created if ninist = 1 (namelist)
-      !!      File 'output.abort.nc' is created in case of abnormal job end
-      !!----------------------------------------------------------------------
-      INTEGER           , INTENT( in ) ::   Kmm              ! time level index
-      CHARACTER (len=* ), INTENT( in ) ::   cdfile_name      ! name of the file created
-      !!
-      INTEGER :: inum
-      !!----------------------------------------------------------------------
-      ! 
-      IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) 'dia_wri_state : single instantaneous ocean state'
-      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~   and forcing fields file created '
-      IF(lwp) WRITE(numout,*) '                and named :', cdfile_name, '...nc'
-
-#if defined key_si3
-     CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE., kdlev = jpl )
-#else
-     CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE. )
-#endif
-
-      CALL iom_rstput( 0, 0, inum, 'votemper', ts(:,:,:,jp_tem,Kmm) )    ! now temperature
-      CALL iom_rstput( 0, 0, inum, 'vosaline', ts(:,:,:,jp_sal,Kmm) )    ! now salinity
-      CALL iom_rstput( 0, 0, inum, 'sossheig', ssh(:,:,Kmm)              )    ! sea surface height
-      CALL iom_rstput( 0, 0, inum, 'vozocrtx', uu(:,:,:,Kmm)                )    ! now i-velocity
-      CALL iom_rstput( 0, 0, inum, 'vomecrty', vv(:,:,:,Kmm)                )    ! now j-velocity
-      IF( ln_zad_Aimp ) THEN
-         CALL iom_rstput( 0, 0, inum, 'vovecrtz', ww + wi        )    ! now k-velocity
-      ELSE
-         CALL iom_rstput( 0, 0, inum, 'vovecrtz', ww             )    ! now k-velocity
-      ENDIF
-      IF( ALLOCATED(ahtu) ) THEN
-         CALL iom_rstput( 0, 0, inum,  'ahtu', ahtu              )    ! aht at u-point
-         CALL iom_rstput( 0, 0, inum,  'ahtv', ahtv              )    ! aht at v-point
-      ENDIF
-      IF( ALLOCATED(ahmt) ) THEN
-         CALL iom_rstput( 0, 0, inum,  'ahmt', ahmt              )    ! ahmt at u-point
-         CALL iom_rstput( 0, 0, inum,  'ahmf', ahmf              )    ! ahmf at v-point
-      ENDIF
-      CALL iom_rstput( 0, 0, inum, 'sowaflup', emp - rnf         )    ! freshwater budget
-      CALL iom_rstput( 0, 0, inum, 'sohefldo', qsr + qns         )    ! total heat flux
-      CALL iom_rstput( 0, 0, inum, 'soshfldo', qsr               )    ! solar heat flux
-      CALL iom_rstput( 0, 0, inum, 'soicecov', fr_i              )    ! ice fraction
-      CALL iom_rstput( 0, 0, inum, 'sozotaux', utau              )    ! i-wind stress
-      CALL iom_rstput( 0, 0, inum, 'sometauy', vtau              )    ! j-wind stress
-      IF(  .NOT.ln_linssh  ) THEN             
-         CALL iom_rstput( 0, 0, inum, 'vovvldep', gdept(:,:,:,Kmm)        )    !  T-cell depth 
-         CALL iom_rstput( 0, 0, inum, 'vovvle3t', e3t(:,:,:,Kmm)          )    !  T-cell thickness  
-      END IF
-      IF( ln_wave .AND. ln_sdw ) THEN
-         CALL iom_rstput( 0, 0, inum, 'sdzocrtx', usd            )    ! now StokesDrift i-velocity
-         CALL iom_rstput( 0, 0, inum, 'sdmecrty', vsd            )    ! now StokesDrift j-velocity
-         CALL iom_rstput( 0, 0, inum, 'sdvecrtz', wsd            )    ! now StokesDrift k-velocity
-      ENDIF
- 
-#if defined key_si3
-      IF( nn_ice == 2 ) THEN   ! condition needed in case agrif + ice-model but no-ice in child grid
-         CALL ice_wri_state( inum )
-      ENDIF
-#endif
-      !
-      CALL iom_close( inum )
-      ! 
-   END SUBROUTINE dia_wri_state
-
-   !!======================================================================
-END MODULE diawri
Index: NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/BENCH/MY_SRC/zdfiwm.F90
===================================================================
--- NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/BENCH/MY_SRC/zdfiwm.F90	(revision 12486)
+++ NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/BENCH/MY_SRC/zdfiwm.F90	(revision 12488)
@@ -48,8 +48,10 @@
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hcri_iwm   ! decay scale for low-mode critical slope dissipation (m)
 
+   !! * Substitutions
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
-   !! $Id: zdfiwm.F90 8093 2017-05-30 08:13:14Z gm $
-   !! Software governed by the CeCILL licence     (./LICENSE)
+   !! $Id: zdfiwm.F90 12377 2020-02-12 14:39:06Z acc $
+   !! Software governed by the CeCILL license (see ./LICENSE)
    !!----------------------------------------------------------------------
 CONTAINS
@@ -121,5 +123,5 @@
       !
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
-      REAL(wp) ::   zztmp        ! scalar workspace
+      REAL(wp) ::   zztmp, ztmp1, ztmp2        ! scalar workspace
       REAL(wp), DIMENSION(jpi,jpj)     ::   zfact       ! Used for vertical structure
       REAL(wp), DIMENSION(jpi,jpj)     ::   zhdep       ! Ocean depth
@@ -147,21 +149,22 @@
       !                       !* Critical slope mixing: distribute energy over the time-varying ocean depth,
       !                                                 using an exponential decay from the seafloor.
-      DO jj = 1, jpj                ! part independent of the level
-         DO ji = 1, jpi
-            zhdep(ji,jj) = gdepw_0(ji,jj,mbkt(ji,jj)+1)       ! depth of the ocean
-            zfact(ji,jj) = rho0 * (  1._wp - EXP( -zhdep(ji,jj) / hcri_iwm(ji,jj) )  )
-            IF( zfact(ji,jj) /= 0._wp )   zfact(ji,jj) = ecri_iwm(ji,jj) / zfact(ji,jj)
-         END DO
-      END DO
-!!gm gde3w ==>>>  check for ssh taken into account.... seem OK gde3w_n=gdept(Kmm) - ssh(Kmm)
-      DO jk = 2, jpkm1              ! complete with the level-dependent part
-         zemx_iwm(:,:,jk) = zfact(:,:) * (  EXP( ( gde3w(:,:,jk  ) - zhdep(:,:) ) / hcri_iwm(:,:) )                      &
-            &                             - EXP( ( gde3w(:,:,jk-1) - zhdep(:,:) ) / hcri_iwm(:,:) )  ) * wmask(:,:,jk)   &
-            &                          / ( gde3w(:,:,jk) - gde3w(:,:,jk-1) )
-
-!!gm delta(gde3w) = e3t(Kmm)  !!  Please verify the grid-point position w versus t-point
+      DO_2D_11_11
+         zhdep(ji,jj) = gdepw_0(ji,jj,mbkt(ji,jj)+1)       ! depth of the ocean
+         zfact(ji,jj) = rho0 * (  1._wp - EXP( -zhdep(ji,jj) / hcri_iwm(ji,jj) )  )
+         IF( zfact(ji,jj) /= 0._wp )   zfact(ji,jj) = ecri_iwm(ji,jj) / zfact(ji,jj)
+      END_2D
+!!gm gde3w ==>>>  check for ssh taken into account.... seem OK gde3w_n=gdept(:,:,:,Kmm) - ssh(:,:,Kmm)
+      DO_3D_11_11( 2, jpkm1 )
+         IF ( zfact(ji,jj) == 0._wp .OR. wmask(ji,jj,jk) == 0._wp ) THEN   ! optimization
+            zemx_iwm(ji,jj,jk) = 0._wp
+         ELSE
+            zemx_iwm(ji,jj,jk) = zfact(ji,jj) * (  EXP( ( gde3w(ji,jj,jk  ) - zhdep(ji,jj) ) / hcri_iwm(ji,jj) )     &
+                 &                               - EXP( ( gde3w(ji,jj,jk-1) - zhdep(ji,jj) ) / hcri_iwm(ji,jj) ) )   &
+                 &                            / ( gde3w(ji,jj,jk) - gde3w(ji,jj,jk-1) )
+         ENDIF
+      END_3D
+!!gm delta(gde3w) = e3t(:,:,:,Kmm)  !!  Please verify the grid-point position w versus t-point
 !!gm it seems to me that only 1/hcri_iwm  is used ==>  compute it one for all
 
-      END DO
 
       !                        !* Pycnocline-intensified mixing: distribute energy over the time-varying 
@@ -177,9 +180,7 @@
          END DO
          !
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = epyc_iwm(ji,jj) / ( rho0 * zfact(ji,jj) )
-            END DO
-         END DO
+         DO_2D_11_11
+            IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = epyc_iwm(ji,jj) / ( rho0 * zfact(ji,jj) )
+         END_2D
          !
          DO jk = 2, jpkm1              ! complete with the level-dependent part
@@ -194,9 +195,7 @@
          END DO
          !
-         DO jj= 1, jpj
-            DO ji = 1, jpi
-               IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = epyc_iwm(ji,jj) / ( rho0 * zfact(ji,jj) )
-            END DO
-         END DO
+         DO_2D_11_11
+            IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = epyc_iwm(ji,jj) / ( rho0 * zfact(ji,jj) )
+         END_2D
          !
          DO jk = 2, jpkm1              ! complete with the level-dependent part
@@ -223,19 +222,18 @@
 !!gm
       !
-      DO jk = 2, jpkm1
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               IF( zfact(ji,jj) /= 0 )   zwkb(ji,jj,jk) = zhdep(ji,jj) * ( zfact(ji,jj) - zwkb(ji,jj,jk) )   &
-                  &                                     * wmask(ji,jj,jk) / zfact(ji,jj)
-            END DO
-         END DO
-      END DO
+      DO_3D_11_11( 2, jpkm1 )
+         IF( zfact(ji,jj) /= 0 )   zwkb(ji,jj,jk) = zhdep(ji,jj) * ( zfact(ji,jj) - zwkb(ji,jj,jk) )   &
+            &                                     * wmask(ji,jj,jk) / zfact(ji,jj)
+      END_3D
       zwkb(:,:,1) = zhdep(:,:) * wmask(:,:,1)
       !
-      zweight(:,:,:) = 0._wp
-      DO jk = 2, jpkm1
-         zweight(:,:,jk) = MAX( 0._wp, rn2(:,:,jk) ) * hbot_iwm(:,:) * wmask(:,:,jk)                    &
-            &   * (  EXP( -zwkb(:,:,jk) / hbot_iwm(:,:) ) - EXP( -zwkb(:,:,jk-1) / hbot_iwm(:,:) )  )
-      END DO
+      DO_3D_11_11( 2, jpkm1 )
+         IF ( rn2(ji,jj,jk) <= 0._wp .OR. wmask(ji,jj,jk) == 0._wp ) THEN   ! optimization
+            zweight(ji,jj,jk) = 0._wp
+         ELSE
+            zweight(ji,jj,jk) = rn2(ji,jj,jk) * hbot_iwm(ji,jj)    &
+               &   * (  EXP( -zwkb(ji,jj,jk) / hbot_iwm(ji,jj) ) - EXP( -zwkb(ji,jj,jk-1) / hbot_iwm(ji,jj) )  )
+         ENDIF
+      END_3D
       !
       zfact(:,:) = 0._wp
@@ -244,9 +242,7 @@
       END DO
       !
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = ebot_iwm(ji,jj) / ( rho0 * zfact(ji,jj) )
-         END DO
-      END DO
+      DO_2D_11_11
+         IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = ebot_iwm(ji,jj) / ( rho0 * zfact(ji,jj) )
+      END_2D
       !
       DO jk = 2, jpkm1              ! complete with the level-dependent part
@@ -276,15 +272,11 @@
       !
       IF( ln_mevar ) THEN              ! Variable mixing efficiency case : modify zav_wave in the
-         DO jk = 2, jpkm1              ! energetic (Reb > 480) and buoyancy-controlled (Reb <10.224 ) regimes
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  IF( zReb(ji,jj,jk) > 480.00_wp ) THEN
-                     zav_wave(ji,jj,jk) = 3.6515_wp * znu_w(ji,jj,jk) * SQRT( zReb(ji,jj,jk) )
-                  ELSEIF( zReb(ji,jj,jk) < 10.224_wp ) THEN
-                     zav_wave(ji,jj,jk) = 0.052125_wp * znu_w(ji,jj,jk) * zReb(ji,jj,jk) * SQRT( zReb(ji,jj,jk) )
-                  ENDIF
-               END DO
-            END DO
-         END DO
+         DO_3D_11_11( 2, jpkm1 )
+            IF( zReb(ji,jj,jk) > 480.00_wp ) THEN
+               zav_wave(ji,jj,jk) = 3.6515_wp * znu_w(ji,jj,jk) * SQRT( zReb(ji,jj,jk) )
+            ELSEIF( zReb(ji,jj,jk) < 10.224_wp ) THEN
+               zav_wave(ji,jj,jk) = 0.052125_wp * znu_w(ji,jj,jk) * zReb(ji,jj,jk) * SQRT( zReb(ji,jj,jk) )
+            ENDIF
+         END_3D
       ENDIF
       !
@@ -296,12 +288,8 @@
          zztmp = 0._wp
 !!gm used of glosum 3D....
-         DO jk = 2, jpkm1
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zztmp = zztmp + e3w(ji,jj,jk,Kmm) * e1e2t(ji,jj)   &
-                     &          * MAX( 0._wp, rn2(ji,jj,jk) ) * zav_wave(ji,jj,jk) * wmask(ji,jj,jk) * tmask_i(ji,jj)
-               END DO
-            END DO
-         END DO
+         DO_3D_11_11( 2, jpkm1 )
+            zztmp = zztmp + e3w(ji,jj,jk,Kmm) * e1e2t(ji,jj)   &
+               &          * MAX( 0._wp, rn2(ji,jj,jk) ) * zav_wave(ji,jj,jk) * wmask(ji,jj,jk) * tmask_i(ji,jj)
+         END_3D
          CALL mpp_sum( 'zdfiwm', zztmp )
          zztmp = rho0 * zztmp ! Global integral of rauo * Kz * N^2 = power contributing to mixing 
@@ -321,13 +309,13 @@
       !      
       IF( ln_tsdiff ) THEN          !* Option for differential mixing of salinity and temperature
-         DO jk = 2, jpkm1              ! Calculate S/T diffusivity ratio as a function of Reb
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zav_ratio(ji,jj,jk) = ( 0.505_wp + 0.495_wp *                                                                  &
-                      &   TANH(    0.92_wp * (   LOG10(  MAX( 1.e-20_wp, zReb(ji,jj,jk) * 5._wp * r1_6 )  ) - 0.60_wp   )    )   &
-                      &                 ) * wmask(ji,jj,jk)
-               END DO
-            END DO
-         END DO
+         ztmp1 = 0.505_wp + 0.495_wp * TANH( 0.92_wp * ( LOG10( 1.e-20_wp ) - 0.60_wp ) )
+         DO_3D_11_11( 2, jpkm1 )
+            ztmp2 = zReb(ji,jj,jk) * 5._wp * r1_6
+            IF ( ztmp2 > 1.e-20_wp .AND. wmask(ji,jj,jk) == 1._wp ) THEN
+               zav_ratio(ji,jj,jk) = 0.505_wp + 0.495_wp * TANH( 0.92_wp * ( LOG10(ztmp2) - 0.60_wp ) )
+            ELSE
+               zav_ratio(ji,jj,jk) = ztmp1 * wmask(ji,jj,jk)
+            ENDIF
+         END_3D
          CALL iom_put( "av_ratio", zav_ratio )
          DO jk = 2, jpkm1           !* update momentum & tracer diffusivity with wave-driven mixing
@@ -364,5 +352,5 @@
       CALL iom_put( "emix_iwm", zemx_iwm )
       
-      IF(ln_ctl)   CALL prt_ctl(tab3d_1=zav_wave , clinfo1=' iwm - av_wave: ', tab3d_2=avt, clinfo2=' avt: ', kdim=jpk)
+      IF(sn_cfctl%l_prtctl)   CALL prt_ctl(tab3d_1=zav_wave , clinfo1=' iwm - av_wave: ', tab3d_2=avt, clinfo2=' avt: ', kdim=jpk)
       !
    END SUBROUTINE zdf_iwm
@@ -395,5 +383,4 @@
       !!              de Lavergne et al. in prep., 2017
       !!----------------------------------------------------------------------
-      INTEGER  ::   ji, jj, jk   ! dummy loop indices
       INTEGER  ::   inum         ! local integer
       INTEGER  ::   ios
