MODULE depwri
   !!======================================================================
   !!                     ***  MODULE  depwri  ***
   !! Ocean diagnostics :  write ocean output files
   !!=====================================================================

   !!----------------------------------------------------------------------
   !! * Modules used
   USE dom_oce         ! ocean space and time domain
   USE in_out_manager
   USE daymod
   USE ioipsl
   USE lib_mpp         ! MPP library

   IMPLICIT NONE
   PRIVATE

   !! * Accessibility
   PUBLIC dep_wri                 ! routine called by step.F90

   !! * Module variables
   INTEGER ::  nid_T, nz_T, nh_T, ndim_T, ndim_hT   ! grid_T file
   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_T

CONTAINS

   INTEGER FUNCTION dep_wri_alloc()
      !!----------------------------------------------------------------------
      INTEGER, DIMENSION(2) :: ierr
      !!----------------------------------------------------------------------
      !
      ierr = 0
      !
      ALLOCATE( ndex_T(jpi*jpj*jpk), STAT=ierr(1) )
         !
      dep_wri_alloc = MAXVAL(ierr)
      IF( lk_mpp )   CALL mpp_sum( dep_wri_alloc )
      !
   END FUNCTION dep_wri_alloc

   SUBROUTINE dep_wri( kt )
      !!---------------------------------------------------------------------
      !!                  ***  ROUTINE dep_wri  ***
      !!                   
      !! ** Purpose :  Write out depths
      !!
      !! * Arguments
      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index

      !! * Local declarations
      CHARACTER (len=40) ::           &
         clop, clmx            ! temporary names
      INTEGER ::   &
         iimi, iima, ipk, it,         &  ! temporary integers
         ijmi, ijma                      !    "          "
      REAL(wp) ::   &
         zsto, zout, zmax,            &  ! temporary scalars
         zjulian, zdt                    !    "         "
      !!----------------------------------------------------------------------

      ! 0. Initialisation
      ! -----------------
      
      ! Define frequency of output and means
      zdt = rdt
      clop = "once"
      zsto = nitend * zdt
      zout = nitend * zdt

      ! Define indices of the horizontal output zoom and vertical limit storage
      iimi = 1      ;      iima = jpi
      ijmi = 1      ;      ijma = jpj
      ipk = jpk

      ! define time axis
      it = kt - nit000 + 1


      ! 1. Define NETCDF files and fields at beginning of first time step
      ! -----------------------------------------------------------------

      IF( kt == nit000 ) THEN

         ! Define the NETCDF
         
         ! Compute julian date from starting date of the run
         CALL ymds2ju( nyear, nmonth, nday, 0.e0, zjulian )
         IF(lwp)WRITE(numout,*)
         IF(lwp)WRITE(numout,*) 'Date 0 used :', nit000, ' YEAR ', nyear,   &
            &                    ' MONTH ', nmonth, ' DAY ', nday, 'Julian day : ', zjulian
         IF(lwp)WRITE(numout,*) ' indexes of zoom = ', iimi, iima, ijmi, ijma,   &
                                 ' limit storage in depth = ', ipk
         
         CALL histbeg( 'depths', jpi, glamt, jpj, gphit,           &  ! Horizontal grid: glamt and gphit
            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
            &          0, zjulian, zdt, nh_T, nid_T, domain_id=nidom )
         CALL histvert( nid_T, "deptht", "Vertical T levels",      &  ! Vertical grid: gdept_0
            &           "m", ipk, gdept_0, nz_T, pdirect="down" )

         ! Ignore that u, v, f are not on T grid - stick em all in one file

         CALL histdef( nid_T, "gdept", "Ocean depth (masked)", "m" ,   &  ! gdept
            &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
         CALL histdef( nid_T, "tmask", "Ocean T grid mask", "1" ,   &  ! tmask
            &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
         CALL histdef( nid_T, "umask", "Ocean U grid mask", "1" ,   &  ! umask
            &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
         CALL histdef( nid_T, "vmask", "Ocean V grid mask", "1" ,   &  ! vmask
            &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
         CALL histdef( nid_T, "fmask", "Ocean F grid mask", "1" ,   &  ! fmask
            &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )

         CALL histend( nid_T )

      ENDIF ! kt == nit000

      ! 2. Start writing data
      ! ---------------------

      ! Write fields on T grid
      CALL histwrite( nid_T, "gdept", it, gdept*tmask, ndim_T , ndex_T  )
      CALL histwrite( nid_T, "tmask", it, tmask, ndim_T , ndex_T  )
      CALL histwrite( nid_T, "umask", it, umask, ndim_T , ndex_T  )
      CALL histwrite( nid_T, "vmask", it, vmask, ndim_T , ndex_T  )
      CALL histwrite( nid_T, "fmask", it, fmask, ndim_T , ndex_T  )

      ! 3. Synchronise and close all files
      ! ---------------------------------------
      IF( kt == nitend ) THEN
         CALL histsync( nid_T )
         CALL histclo( nid_T )
      ENDIF

   END SUBROUTINE dep_wri

   !!======================================================================
END MODULE depwri