Index: /branches/UKMO/AMM15_v3_6_STABLE_package_collate/NEMOGCM/NEMO/OPA_SRC/DIA/diaopfoam.F90
===================================================================
--- /branches/UKMO/AMM15_v3_6_STABLE_package_collate/NEMOGCM/NEMO/OPA_SRC/DIA/diaopfoam.F90	(revision 12569)
+++ /branches/UKMO/AMM15_v3_6_STABLE_package_collate/NEMOGCM/NEMO/OPA_SRC/DIA/diaopfoam.F90	(revision 12570)
@@ -15,4 +15,5 @@
    USE diurnal_bulk
    USE cool_skin
+   USE ioipsl
 
 
@@ -21,4 +22,5 @@
 
    LOGICAL , PUBLIC ::   ln_diaopfoam     !: Diaopfoam output 
+   LOGICAL , PUBLIC ::   ln_diaopfoam_Tzero     !: Diaopfoam first time step output
    PUBLIC   dia_diaopfoam_init            ! routine called by nemogcm.F90
    PUBLIC   dia_diaopfoam                 ! routine called by diawri.F90
@@ -49,9 +51,10 @@
       INTEGER            ::   ierror   ! local integer
       !!
-      NAMELIST/nam_diadiaop/ ln_diaopfoam
+      NAMELIST/nam_diadiaop/ ln_diaopfoam,ln_diaopfoam_Tzero
       !!
       !!----------------------------------------------------------------------
       !
       ln_diaopfoam = .false.         ! default value for diaopfoam stream
+      ln_diaopfoam_Tzero = .false.         ! default value for diaopfoam Tzero stream
       REWIND ( numnam_ref )              ! Read Namelist nam_diadiaop in reference namelist : 3D hourly diagnostics
       READ   ( numnam_ref, nam_diadiaop, IOSTAT=ios, ERR= 901 )
@@ -68,9 +71,10 @@
          WRITE(numout,*) '~~~~~~~~~~~~'
          WRITE(numout,*) '   Namelist nam_diadiaop : set diaopfoam outputs '
-         WRITE(numout,*) '      Switch for diaopfoam diagnostics (T) or not (F)  ln_diaopfoam  = ', ln_diaopfoam
+         WRITE(numout,*) '      Switch for diaopfoam diagnostics (T) or not (F)                 ln_diaopfoam        = ', ln_diaopfoam
+         WRITE(numout,*) '      Switch for diaopfoam first timestep diagnostics (T) or not (F)  ln_diaopfoam_Tzero  = ', ln_diaopfoam_Tzero
       ENDIF
     END SUBROUTINE dia_diaopfoam_init
 
-    SUBROUTINE dia_diaopfoam
+    SUBROUTINE dia_diaopfoam( kt )
       !!----------------------------------------------------------------------
       !!                 ***  ROUTINE dia_diaopfoam  ***
@@ -82,4 +86,8 @@
       !!         
       !!--------------------------------------------------------------------
+      IMPLICIT NONE
+
+      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
+
       REAL(wp), DIMENSION(jpi,jpj)   :: zw2d       ! 2D workspace
       REAL(wp), DIMENSION(jpi,jpj,jpk) :: zw3d       ! 3D workspace
@@ -94,6 +102,12 @@
 
       zmdi=1.e+20                               !  missing data indicator for masking
+
       ! Diaopfoam stream if needed
       IF (ln_diaopfoam) THEN
+         IF ( kt .eq. nn_it000 .AND. ln_diaopfoam_Tzero ) THEN
+            IF(lwp) WRITE(numout,*) 'diaopfoam: writing T0 at kt = ', kt
+            CALL dia_diaopfoam_zero( 'Tzero', kt )
+         ENDIF
+         
          CALL theta2t
          CALL iom_put( "insitut_op" , rinsitu_t(:,:,:)                      )    ! insitu temperature
@@ -101,6 +115,6 @@
          CALL iom_put( "soce_op"   , tsn(:,:,:,jp_sal)                     )    ! salinity
          IF (ln_diurnal) THEN
-            CALL iom_put( "sst_wl_op"   , x_dsst                    )    ! warm layer 
-            CALL iom_put( "sst_cs_op"   , x_csdsst                  )    ! cool skin 
+            CALL iom_put( "sst_wl_op"   , x_dsst                             )    ! warm layer 
+            CALL iom_put( "sst_cs_op"   , x_csdsst                           )    ! cool skin 
          ENDIF
          zw2d(:,:)=sshn(:,:)*tmask(:,:,1) + zmdi*(1.0-tmask(:,:,1))
@@ -175,3 +189,189 @@
    END SUBROUTINE calc_max_cur
 
+   SUBROUTINE dia_diaopfoam_zero( cdfile_name, kt )
+      !!---------------------------------------------------------------------
+      !!                 ***  ROUTINE dia_diaopfoam_zero  ***
+      !!        
+      !! ** Purpose :   create a NetCDF file named cdfile_name which contains 
+      !!      the instantaneous ocean state at the first time step.
+      !!
+      !! ** Method  :   NetCDF files using ioipsl
+      !!----------------------------------------------------------------------
+      CHARACTER (len=* ), INTENT( in ) ::   cdfile_name      ! name of the file created
+      INTEGER           , INTENT( in ) ::   kt               ! ocean time-step index
+      !! 
+      CHARACTER (len=32) :: clhstnam
+      CHARACTER (len=40) :: clop
+      INTEGER  ::   iimi, iima, ipk, ijmi, ijma          ! local integers
+      INTEGER  ::   nid_T, nz_T, nh_T, ndim_T, ndim_hT   ! 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  ::   id_i , nz_i, nh_i       
+      INTEGER, DIMENSION(1) ::   idex                    ! local workspace
+      INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hT, ndex_hU, ndex_hV
+      INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_T, ndex_U, ndex_V
+      INTEGER :: ierr
+      INTEGER :: jkbot, jj, ji
+      REAL(wp) :: zsto, zout, zmax, zjulian, zdt
+      REAL(wp) :: zmdi
+      REAL(wp), POINTER, DIMENSION(:,:)   :: zwu
+      REAL(wp), POINTER, DIMENSION(:,:)   :: zwv
+      REAL(wp), POINTER, DIMENSION(:,:)   :: zwz
+      REAL(wp), DIMENSION(jpi,jpj) :: zw2d           ! 2D workspace
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zw3d       ! 3D workspace
+      REAL(wp), DIMENSION(jpi,jpj) :: z2d
+
+      !!----------------------------------------------------------------------
+
+      ! -----------------------------------------------------------------
+      ! 0. Allocations
+      ! -----------------------------------------------------------------
+      ierr = 0
+      ALLOCATE( ndex_hT(jpi*jpj) , ndex_T(jpi*jpj*jpk) ,     &
+         &      ndex_hU(jpi*jpj) , ndex_U(jpi*jpj*jpk) ,     &
+         &      ndex_hV(jpi*jpj) , ndex_V(jpi*jpj*jpk) , STAT=ierr )
+      IF( lk_mpp )   CALL mpp_sum( ierr )
+         IF( ierr /= 0 ) THEN
+            CALL ctl_stop('dia_diaopfoam_zero: failed to allocate arrays')
+            RETURN
+      ENDIF
+      CALL wrk_alloc( jpi , jpj      , zwu )
+      CALL wrk_alloc( jpi , jpj      , zwv )
+      CALL wrk_alloc( jpi , jpj      , zwz )
+
+      zmdi=1.e+20                               !  missing data indicator for masking
+
+      ! -----------------------------------------------------------------
+      ! 1. Define NETCDF files and fields at beginning of first time step
+      ! -----------------------------------------------------------------
+
+      ! Define indices of the horizontal output zoom and vertical limit storage
+      iimi = 1      ;      iima = jpi
+      ijmi = 1      ;      ijma = jpj
+      ipk = jpk
+
+      ! Define frequency of output and means
+      zdt  = rdt
+      zsto = rdt
+      clop = "inst(x)"           ! no use of the mask value (require less cpu time)
+      zout = rdt
+      zmax = ( nitend - nit000 + 1 ) * zdt
+
+      ! Compute julian date from starting date of the run
+      CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian )         ! time axis 
+      zjulian = zjulian - adatrj   !   set calendar origin to the beginning of the experiment
+
+      ! Define the T grid FILE ( nid_T )
+      clhstnam = TRIM(cdfile_name)//".grid_T"
+      IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam    ! filename
+      CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,           &  ! Horizontal grid: glamt and gphit
+         &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
+         &          nit000-1, zjulian, zdt, nh_T, nid_T, domain_id=nidom, snc4chunks=snc4set )
+      CALL histvert( nid_T, "deptht", "Vertical T levels",      &  ! Vertical grid: gdept
+         &           "m", ipk, gdept_1d, nz_T, "down" )
+      !                                                            ! Index of ocean points
+      CALL wheneq( jpi*jpj*ipk, tmask, 1, 1., ndex_T , ndim_T  )      ! volume
+      CALL wheneq( jpi*jpj    , tmask, 1, 1., ndex_hT, ndim_hT )      ! surface
+
+      ! Define the U grid FILE ( nid_U )
+      clhstnam = TRIM(cdfile_name)//".grid_U"
+      IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam    ! filename
+      CALL histbeg( clhstnam, jpi, glamu, jpj, gphiu,           &  ! Horizontal grid: glamu and gphiu
+         &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
+         &          nit000-1, zjulian, zdt, nh_U, nid_U, domain_id=nidom, snc4chunks=snc4set )
+      CALL histvert( nid_U, "depthu", "Vertical U levels",      &  ! Vertical grid: gdepu
+         &           "m", ipk, gdept_1d, nz_U, "down" )
+      !                                                            ! Index of ocean points
+      CALL wheneq( jpi*jpj*ipk, umask, 1, 1., ndex_U , ndim_U  )      ! volume
+      CALL wheneq( jpi*jpj    , umask, 1, 1., ndex_hU, ndim_hU )      ! surface
+
+      ! Define the V grid FILE ( nid_V )
+      clhstnam = TRIM(cdfile_name)//".grid_V"
+      IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
+      CALL histbeg( clhstnam, jpi, glamv, jpj, gphiv,           &  ! Horizontal grid: glamv and gphiv
+         &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
+         &          nit000-1, zjulian, zdt, nh_V, nid_V, domain_id=nidom, snc4chunks=snc4set )
+      CALL histvert( nid_V, "depthv", "Vertical V levels",      &  ! Vertical grid : gdepv
+         &          "m", ipk, gdept_1d, nz_V, "down" )
+      !                                                            ! Index of ocean points
+      CALL wheneq( jpi*jpj*ipk, vmask, 1, 1., ndex_V , ndim_V  )      ! volume
+      CALL wheneq( jpi*jpj    , vmask, 1, 1., ndex_hV, ndim_hV )      ! surface
+
+
+      ! -----------------------------------------------------------------
+      ! 2. Declare all the output fields as NETCDF variables
+      ! -----------------------------------------------------------------
+
+      !                                                                                     !!! nid_T : 3D
+      !CALL histdef( nid_T, "votempis", "Insitu Temperature" , "C"      ,       &  ! 
+      !   &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
+      CALL histdef( nid_T, "votemper", "Temperature"         , "C"     ,       &  ! tn
+         &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
+      CALL histdef( nid_T, "vosaline", "Salinity"            , "PSU"   ,       & ! sn
+         &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
+      CALL histdef( nid_T, "sossheig", "Sea Surface Height"  , "m"     ,       &  ! sshn
+         &          jpi, jpj, nh_T, 1  , 1, 1  , nz_T, 32, clop, zsto, zout )
+      CALL histdef( nid_T, "votempis", "Insitu Temperature"  , "C"     ,       &  ! rinsitu_t
+         &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
+      CALL histdef( nid_T, "maxu"    , "Max Zonal Current"   , "m/s"     ,       &  ! zwu
+         &          jpi, jpj, nh_T, 1  , 1, 1  , nz_T, 32, clop, zsto, zout )
+      CALL histdef( nid_T, "maxv"    , "Max Meridional Current" , "m/s"     ,    &  ! zwv
+         &          jpi, jpj, nh_T, 1  , 1, 1  , nz_T, 32, clop, zsto, zout )
+      CALL histdef( nid_T, "maxz"    , "Max Current Depth"   , "m/s"     ,       &  ! zwz
+         &          jpi, jpj, nh_T, 1  , 1, 1  , nz_T, 32, clop, zsto, zout )
+      CALL histdef( nid_T, "sbt"       , "Bottom Temperature"  , "C"     ,       &  ! sbt
+         &          jpi, jpj, nh_T, 1  , 1, 1  , nz_T, 32, clop, zsto, zout )
+      CALL histend( nid_T, snc4chunks=snc4set )
+
+      !                                                                                     !!! nid_U : 3D
+      CALL histdef( nid_U, "vozocrtx", "Zonal Current"                      , "m/s"    ,   &  ! un
+         &          jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
+      CALL histend( nid_U, snc4chunks=snc4set )
+
+      !                                                                                     !!! nid_V : 3D
+      CALL histdef( nid_V, "vomecrty", "Meridional Current"                 , "m/s"    ,   &  ! vn
+         &          jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
+      CALL histend( nid_V, snc4chunks=snc4set )
+
+
+      ! -----------------------------------------------------------------
+      ! 3. Write the data
+      ! -----------------------------------------------------------------
+
+      idex(1) = 1
+      CALL histwrite( nid_T, "votemper", kt, tsn(:,:,:,jp_tem), jpi*jpj*jpk, idex )    ! now temperature
+      CALL histwrite( nid_T, "vosaline", kt, tsn(:,:,:,jp_sal), jpi*jpj*jpk, idex )    ! now salinity
+      CALL histwrite( nid_T, "sossheig", kt, sshn             , jpi*jpj    , idex )    ! sea surface height
+      CALL theta2t
+      CALL histwrite( nid_T, "votempis", kt, rinsitu_t(:,:,:) , jpi*jpj*jpk, idex )    ! now insitu-temperature
+      CALL calc_max_cur(zwu,zwv,zwz,zmdi)
+      CALL lbc_lnk( zwu, 'T', 1. )
+      CALL lbc_lnk( zwv, 'T', 1. )
+      CALL lbc_lnk( zwz, 'T', 1. )
+      CALL histwrite( nid_T, "maxu"    , kt, zwu              , jpi*jpj    , idex )    ! max u-current
+      CALL histwrite( nid_T, "maxv"    , kt, zwv              , jpi*jpj    , idex )    ! max v-current
+      CALL histwrite( nid_T, "maxz"    , kt, zwz              , jpi*jpj    , idex )    ! max current depth
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            jkbot = mbkt(ji,jj)
+            z2d(ji,jj) = tsn(ji,jj,jkbot,jp_tem)
+         END DO
+      END DO
+      CALL histwrite( nid_T, "sbt"     , kt, z2d              , jpi*jpj    , idex )    ! sbt
+      ! U file
+      CALL histwrite( nid_U, "vozocrtx", kt, un              , jpi*jpj*jpk, idex )     ! now i-velocity
+      ! V file
+      CALL histwrite( nid_V, "vomecrty", kt, vn              , jpi*jpj*jpk, idex )     ! now j-velocity
+
+
+      ! -----------------------------------------------------------------
+      ! 4. Close the files
+      ! -----------------------------------------------------------------
+
+      CALL histclo( nid_T )
+      CALL histclo( nid_U )
+      CALL histclo( nid_V )
+
+   END SUBROUTINE dia_diaopfoam_zero
+
 END MODULE diaopfoam
Index: /branches/UKMO/AMM15_v3_6_STABLE_package_collate/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90
===================================================================
--- /branches/UKMO/AMM15_v3_6_STABLE_package_collate/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90	(revision 12569)
+++ /branches/UKMO/AMM15_v3_6_STABLE_package_collate/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90	(revision 12570)
@@ -405,5 +405,5 @@
       ENDIF
       IF (ln_diaopfoam) THEN
-         CALL dia_diaopfoam
+         CALL dia_diaopfoam( kt )
       ENDIF
       !
