#include "options.inc"


c
c  local CPP defines for writing NetCDF snapshots
c
c  write UV 
#define write_UV
c  wrtie T and S
#define write_TS
c  write SSH
c#define write_eta
c#define write_ubar

c  write wind stress components
#define write_tau
c  write surface heat/salt fluxes
#define write_tracer_fluxes

c  write the forcing of the barotr. mode
c#define write_barotropic_forcing
c  write the restoring b.c. for TS at the open boundaries
c#define write_restore_TS_obc
c  write blue density
c#define write_blue_density
c  write also the blue parameter
c#define write_blue_parameter
c  write vertical velocity
#define write_adv_vbt
c  write T/S climatologies and restoring time scales
c#define write_surface_restoring
c  write the forcing for the KT model
c#define write_ustar
c  write the penetration depth of wind stirring
c#define write_ktmix_dml
c  write the hor. diffusivity/viscosity
c#define write_am_scale
c#define write_ambi_scale
c#define write_ah_scale
c#define write_ahbi_scale
c  write the GM velocities
c#define write_gm_velocity
c  write some informations from convection algorythm
c#define write_ventilated_depth
c#define write_convection_below_ml
c  write vertical diffusivity
c#define write_vert_diff
c#define write_vert_visc
c#define write_tkemix
c#define write_Nsqr





      subroutine init_snap_cdf
c
c-----------------------------------------------------------------------
c     initialize NetCDF snapshot file
c-----------------------------------------------------------------------
c
      use spflame_module
      implicit none
#ifdef netcdf_diagnostics
#include "netcdf.inc"
      integer ncid,iret,i,j,k,n
      integer lon_tdim,lon_udim,lat_tdim,lat_udim,itimedim
      integer depth_wdim,depth_wid, depth_tdim,depth_tid
      integer lon_tid,lon_uid,lat_tid,lat_uid,itimeid
      integer ipsiid,itopoid,itauxid,itauyid,ihtpid
      integer rhoid,brhoid ,zu1id,zu2id,guessid,ptdid
      integer uid,vid,trid(nt),wid,balphaid
      integer tobcnid,sobcnid, tobcsid,sobcsid
      integer tobcwid,sobcwid, tobceid,sobceid
      integer stfid(nt),ivz,ivz2,ekeid
      integer ihfrestid,isfrestid,iustarid,idmlid
      integer ihfclimid,isfclimid,sshid,ubarid,vbarid
      integer iamid,iahid,iahbiid,iambiid,adv_sid
      integer ugmid,vgmid,wgmid,diffvid,viscvid,Nsqrid
      integer dims(4), corner(4), edges(4)
      character name*24, unit*16, text*80
      real :: spval=-9.9e12
#ifdef netcdf_real4
      real (kind=4), allocatable :: var4(:)
#else
      real         , allocatable :: var4(:)
#endif
       integer :: year,month,day
c
c    overiding namelist input for file name
c
       year  = get_year(get_current_time())
       month = get_month_of_year(get_current_time())
       day   = get_day_of_month(get_current_time())
       write(snap_file,'("snap_",i3,"_y",i4,"m",i2,"d",i2,".cdf")') 
     &            sub_domain,year,month,day
       call replace_space_zero(snap_file)

      if (my_pe == 0) then
      print*,' initializing NetCDF output file ',
     &        snap_file(1:len_trim(snap_file))

      ncid = nccre (snap_file, NCCLOB, iret)
      iret=nf_set_fill(ncid, NF_NOFILL, iret)
      call store_info_cdf(ncid)
c     dimensions
      lon_tdim  = ncddef(ncid, 'Longitude_t', imt, iret)
      Lon_udim  = ncddef(ncid, 'Longitude_u', imt, iret)
      Lat_tdim  = ncddef(ncid, 'Latitude_t',  jmt, iret)
      Lat_udim  = ncddef(ncid, 'Latitude_u',  jmt, iret)
      depth_wdim = ncddef(ncid, 'depth_w',  km, iret)
      depth_tdim = ncddef(ncid, 'depth_t',  km, iret)
      iTimedim  = ncddef(ncid, 'Time', nf_unlimited, iret)
c     grid variables
      dims(1)  = Lon_tdim
      Lon_tid  = ncvdef (ncid,'Longitude_t',NCFLOAT,1,dims,iret)
      dims(1)  = Lon_udim
      Lon_uid  = ncvdef (ncid,'Longitude_u',NCFLOAT,1,dims,iret)
      dims(1)  = Lat_tdim
      Lat_tid  = ncvdef (ncid,'Latitude_t', NCFLOAT,1,dims,iret)
      dims(1)  = Lat_udim
      Lat_uid  = ncvdef (ncid,'Latitude_u', NCFLOAT,1,dims,iret)
      dims(1)  = iTimedim
      iTimeid   = ncvdef(ncid,'Time',       NCFLOAT,1,dims,iret)
      dims(1)  = depth_wdim
      depth_wid = ncvdef (ncid,'depth_w', NCFLOAT,1,dims,iret)
      dims(1)  = depth_tdim
      depth_tid = ncvdef (ncid,'depth_t', NCFLOAT,1,dims,iret)
c     2 dim variables on t grid
      dims(1) = Lon_tdim; dims(2) = Lat_tdim; dims(3) = iTimedim
      ipsiid  = ncvdef (ncid,'psi', NCFLOAT,3,dims,iret)
      itopoid = ncvdef (ncid,'topo', NCFLOAT,2,dims,iret)
#ifdef partial_cell
      ihtpid = ncvdef (ncid,'htp', NCFLOAT,2,dims,iret)
#endif
#ifdef write_tracer_fluxes
      do n=1,nt
       name='heat_fl';if (n==2) name='salt_fl'
       if (n>2) then
         write(name, '("tracer_",i2,"_sflux")') n
         call replace_space_zero(name)
       endif
       stfid(n) = ncvdef (ncid,name, NCFLOAT,3,dims,iret)
      enddo
#endif
#ifdef write_surface_restoring
      ihfrestid = ncvdef (ncid,'T_rest', NCFLOAT,3,dims,iret)
      ihfclimid = ncvdef (ncid,'T_clim', NCFLOAT,3,dims,iret)
      isfrestid = ncvdef (ncid,'S_rest', NCFLOAT,3,dims,iret)
      isfclimid = ncvdef (ncid,'S_clim', NCFLOAT,3,dims,iret)
#endif
#if defined write_eta 
      if (enable_freesurf) 
     &      sshid = ncvdef (ncid,'eta', NCFLOAT,3,dims,iret)
#endif
#if defined write_ubar
      if (enable_freesurf)  then
         ubarid = ncvdef (ncid,'ubar', NCFLOAT,3,dims,iret)
         vbarid = ncvdef (ncid,'vbar', NCFLOAT,3,dims,iret)
      endif
#endif

      if (enable_ktmix.or.enable_tkemix) then
#if defined write_ustar 
       iustarid = ncvdef (ncid,'ustar', NCFLOAT,3,dims,iret)
#endif
      endif
      if (enable_ktmix) then
#ifdef write_ktmix_dml
       idmlid = ncvdef (ncid,'dml', NCFLOAT,3,dims,iret)
#endif
      endif
      if (enable_expl_convection) then
#ifdef write_ventilated_depth
       ivz = ncvdef (ncid,'convect', NCFLOAT,3,dims,iret)
#endif
#ifdef write_convection_below_ml
       ivz2 = ncvdef (ncid,'convect2', NCFLOAT,3,dims,iret)
#endif
      endif
c     3 dim variables
      dims(1) = Lon_udim; dims(2) = Lat_udim
      dims(3) = depth_tdim; dims(4) = iTimedim
#ifdef write_UV
      uid=ncvdef (ncid,'u', NCFLOAT,4,dims,iret)
      vid=ncvdef (ncid,'v', NCFLOAT,4,dims,iret)
#endif
#if defined write_gm_velocity 
      if (enable_diffusion_isoneutral.or.
     &    enable_diffusion_isopycnic) then
       dims=(/lon_udim,lat_tdim,depth_tdim,itimedim/)
       ugmid=ncvdef (ncid,'u_gm', NCFLOAT,4,dims,iret)
       dims=(/lon_tdim,lat_udim,depth_tdim,itimedim/)
       vgmid=ncvdef (ncid,'v_gm', NCFLOAT,4,dims,iret)
      endif
#endif
c     2 dim variables on u grid
      dims(1) = Lon_udim; dims(2) = Lat_udim; dims(3) = iTimedim
#ifdef write_tau
      itauxid=ncvdef (ncid,'taux', NCFLOAT,3,dims,iret)
      itauyid=ncvdef (ncid,'tauy', NCFLOAT,3,dims,iret)
#endif
#ifdef write_barotropic_forcing
      zu1id=ncvdef (ncid,'zu1', NCFLOAT,3,dims,iret)
      zu2id=ncvdef (ncid,'zu2', NCFLOAT,3,dims,iret)
      dims(1) = Lon_tdim; dims(2) = Lat_tdim; dims(3) = iTimedim
      guessid=ncvdef (ncid,'guess', NCFLOAT,3,dims,iret)
      ptdid=ncvdef (ncid,'ptd', NCFLOAT,3,dims,iret)
      dims(1) = Lon_udim; dims(2) = Lat_udim; dims(3) = iTimedim
#endif
#if defined write_am_scale 
      iamid=ncvdef (ncid,'am', NCFLOAT,2,dims,iret)
#endif
#if defined write_ambi_scale 
      iambiid=ncvdef (ncid,'ambi', NCFLOAT,2,dims,iret)
#endif
#if defined write_ah_scale 
      iahid=ncvdef (ncid,'ah', NCFLOAT,2,dims,iret)
#endif
#if defined write_ahbi_scale 
      iahbiid=ncvdef (ncid,'ahbi', NCFLOAT,2,dims,iret)
#endif
      dims(1) = Lon_tdim; dims(2) = Lat_tdim; dims(3) = depth_tdim
      dims(4) = iTimedim
#ifdef write_TS
      do n=1,nt
       name='temp';if (n==2) name='salt'
       if (n>2) then
         write(name, '("tracer_",i2)') n
         call replace_space_zero(name)
       endif
       trid(n)=ncvdef (ncid,name, NCFLOAT,4,dims,iret)
      enddo
#endif
      if (enable_blue) then
#if defined write_blue_density 
       rhoid=ncvdef (ncid,'rho', NCFLOAT,4,dims,iret)
       brhoid=ncvdef (ncid,'blue_rho', NCFLOAT,4,dims,iret)
#endif
#if defined write_blue_parameter
       balphaid=ncvdef (ncid,'blue_alpha', NCFLOAT,3,dims,iret)
#endif
      endif
#ifdef write_adv_vbt
      dims=(/lon_tdim,lat_tdim,depth_wdim,itimedim/)
      wid=ncvdef (ncid,'w', NCFLOAT,4,dims,iret)
#endif
#if defined write_gm_velocity 
      if (enable_diffusion_isoneutral.or.
     &    enable_diffusion_isopycnic) then
       dims=(/lon_tdim,lat_tdim,depth_wdim,itimedim/)
       wgmid=ncvdef (ncid,'w_gm', NCFLOAT,4,dims,iret)
      endif
#endif
      dims(1) = Lon_tdim; dims(2) = Lat_tdim; dims(3) = depth_wdim
      dims(4) = iTimedim
#ifdef write_vert_diff
      diffvid=ncvdef (ncid,'diff_cbt', NCFLOAT,4,dims,iret)
#endif
#ifdef write_vert_visc
      dims(1) = Lon_udim; dims(2) = Lat_udim; dims(3) = depth_wdim
      dims(4) = iTimedim
      viscvid=ncvdef (ncid,'visc_cbu', NCFLOAT,4,dims,iret)
#endif
#ifdef write_Nsqr
      dims(1) = Lon_tdim; dims(2) = Lat_tdim; dims(3) = depth_wdim
      dims(4) = iTimedim
      Nsqrid=ncvdef (ncid,'Nsqr', NCFLOAT,4,dims,iret)
#endif

#ifdef write_tkemix
      dims(1) = Lon_tdim; dims(2) = Lat_tdim; dims(3) = depth_wdim
      dims(4) = iTimedim
      if (enable_tkemix) then
       ekeid=ncvdef (ncid,'tke', NCFLOAT,4,dims,iret)
      endif
#endif

#ifdef write_restore_TS_obc
      dims(1) = Lon_tdim; dims(2) = depth_tdim; dims(3) = iTimedim
      tobcnid=ncvdef (ncid,'t_obc_north', NCFLOAT,3,dims,iret)
      sobcnid=ncvdef (ncid,'s_obc_north', NCFLOAT,3,dims,iret)
      tobcsid=ncvdef (ncid,'t_obc_south', NCFLOAT,3,dims,iret)
      sobcsid=ncvdef (ncid,'s_obc_south', NCFLOAT,3,dims,iret)
      dims(1) = Lat_tdim
      tobcwid=ncvdef (ncid,'t_obc_west', NCFLOAT,3,dims,iret)
      sobcwid=ncvdef (ncid,'s_obc_west', NCFLOAT,3,dims,iret)
      tobceid=ncvdef (ncid,'t_obc_east', NCFLOAT,3,dims,iret)
      sobceid=ncvdef (ncid,'s_obc_east', NCFLOAT,3,dims,iret)
#endif
c     attributes of the grid
      name = 'Longitude on T grid     '; unit = 'degrees_W       '
      call ncaptc(ncid, Lon_tid, 'long_name', NCCHAR, 24, name, iret) 
      call ncaptc(ncid, Lon_tid, 'units',     NCCHAR, 16, unit, iret) 
      name = 'Longitude on U grid     '; unit = 'degrees_W       '
      call ncaptc(ncid, Lon_uid, 'long_name', NCCHAR, 24, name, iret) 
      call ncaptc(ncid, Lon_uid, 'units',     NCCHAR, 16, unit, iret) 
      name = 'Latitude on T grid      '; unit = 'degrees_N       '
      call ncaptc(ncid, Lat_tid, 'long_name', NCCHAR, 24, name, iret) 
      call ncaptc(ncid, Lat_tid, 'units',     NCCHAR, 16, unit, iret) 
      name = 'Latitude on U grid      '; unit = 'degrees_N       '
      call ncaptc(ncid, Lat_uid, 'long_name', NCCHAR, 24, name, iret) 
      call ncaptc(ncid, Lat_uid, 'units',     NCCHAR, 16, unit, iret) 
      name = 'Depth of T Grid points  '; unit = 'm '
      call ncaptc(ncid, depth_tid, 'long_name', NCCHAR, 24, name, iret) 
      call ncaptc(ncid, depth_tid, 'units',     NCCHAR, 16, unit, iret) 
      name = 'Depth of W Grid points  '; unit = 'm '
      call ncaptc(ncid, depth_wid, 'long_name', NCCHAR, 24, name, iret) 
      call ncaptc(ncid, depth_wid, 'units',     NCCHAR, 16, unit, iret) 
      name = 'Time                    '; unit = 'days            '
      call ncaptc(ncid, iTimeid, 'long_name', NCCHAR, 24, name, iret) 
      call ncaptc(ncid, iTimeid, 'units',     NCCHAR, 16, unit, iret) 
      call ncaptc(ncid, iTimeid,'time_origin',NCCHAR, 20,
     &  '31-DEC-1899 00:00:00', iret)
c     attributes of variables
      name = 'Topography            '; unit = 'm    '
      call dvcdf(ncid,itopoid,name,24,unit,16,spval)
#ifdef partial_cell
      name = 'Depth of partial cells'; unit = 'm    '
      call dvcdf(ncid,ihtpid,name,24,unit,16,spval)
#endif
      name = 'STREAMFUNCTION        '; unit = 'Sv    '
      call dvcdf(ncid,ipsiid,name,24,unit,16,spval)
#ifdef write_TS
      do n=1,nt
       name='Temperature';if (n==2) name='Salinity'
       if (n>2) write(name,'("Tracer Nr.",i2)') n
       unit='deg C'; if (n==2) unit='PSU'
       if (n>2) unit='  '
       call dvcdf(ncid,trid(n),name,24,unit,16,spval)
      enddo
#endif
#ifdef write_UV
      name = 'ZONAL VELOCITY     '; unit = 'm/s'
      call dvcdf(ncid,uid,name,24,unit,16,spval)
      name = 'MERIDIONAL VELOCITY     '; unit = 'm/s'
      call dvcdf(ncid,vid,name,24,unit,16,spval)
#endif
#if defined write_gm_velocity 
      if (enable_diffusion_isoneutral.or.
     &    enable_diffusion_isopycnic) then
       name = 'ZONAL BOLUS VELOCITY '; unit = 'm/s'
       call dvcdf(ncid,ugmid,name,24,unit,16,spval)
       name = 'MERIDIONAL BOLUS VELOCITY'; unit = 'm/s'
       call dvcdf(ncid,vgmid,name,24,unit,16,spval)
       name = 'VERTICAL BOLUS VELOCITY'; unit = 'm/s'
       call dvcdf(ncid,wgmid,name,24,unit,16,spval)
      endif
#endif
      if (enable_blue) then
#if defined write_blue_density 
       name = 'Density     '; unit = 'kg/m^3'
       call dvcdf(ncid,rhoid,name,24,unit,16,spval)
       name = 'Density from sub domain '; unit = 'kg/m^3'
       call dvcdf(ncid,brhoid,name,24,unit,16,spval)
#endif
#if defined write_blue_parameter
       name = 'Blue parameter'; unit = '--'
       call dvcdf(ncid,balphaid,name,24,unit,16,spval)
#endif
      endif
#ifdef write_adv_vbt
      name = 'Vertical Velocity  '; unit = 'm/s'
      call dvcdf(ncid,wid,name,24,unit,16,spval)
#endif
#ifdef write_vert_diff
      name = 'Vertical Diffusivity'; unit = 'cm^2/s'
      call dvcdf(ncid,diffvid,name,24,unit,16,spval)
#endif
#ifdef write_vert_visc
      name = 'Vertical Viscosity'; unit = 'cm^2/s'
      call dvcdf(ncid,viscvid,name,24,unit,16,spval)
#endif
#ifdef write_Nsqr
      name = 'Brunt-Vaisala frequency'; unit = '1/s^2'
      call dvcdf(ncid,Nsqrid,name,24,unit,16,spval)
#endif
#ifdef write_tkemix
      if (enable_tkemix) then
       name = 'Turbulent kinetic energy'; unit = 'm^2/s^2'
       call dvcdf(ncid,ekeid,name,24,unit,16,spval)
      endif
#endif
#ifdef write_tracer_fluxes
      do n=1,nt
       name = 'Surface heat flux'; unit = 'W/m^2'
       if (n==2) then
        name = 'Surface salt flux'; unit = 'kg/m**2/s'
       endif
       if (n>2) then
        write(name,'("Surface flux tracer ",i2)') n
        call replace_space_zero(name)
        unit = '  '
       endif
       call dvcdf(ncid,stfid(n),name,24,unit,16,spval)
      enddo
#endif
#ifdef write_surface_restoring
      name = 'Temp. restoring coeff.'; unit = 'W/m^2/K'
      call dvcdf(ncid,ihfrestid,name,24,unit,16,spval)
      name = 'Sal. restoring coeff.'; unit = 'kg/m**2/s'
      call dvcdf(ncid,isfrestid,name,24,unit,16,spval)
      name = 'Temp. restoring clim.'; unit = 'deg C'
      call dvcdf(ncid,ihfclimid,name,24,unit,16,spval)
      name = 'Sal. restoring clim.'; unit = 'psu'
      call dvcdf(ncid,isfclimid,name,24,unit,16,spval)
#endif
#if defined  write_eta 
      if (enable_freesurf) then
       name = 'Sea surface height'; unit = 'm'
       call dvcdf(ncid,sshid,name,24,unit,16,spval)
      endif
#endif
#if defined  write_ubar
      if (enable_freesurf) then
       name = 'Barotropic velocity'; unit = 'm/s'
       call dvcdf(ncid,ubarid,name,24,unit,16,spval)
       name = 'Barotropic velocity'; unit = 'm/s'
       call dvcdf(ncid,vbarid,name,24,unit,16,spval)
      endif
#endif
      if (enable_ktmix.or.enable_tkemix) then
#if defined write_ustar
       name = 'U^star'; unit = '(cm/s)^3'
       call dvcdf(ncid,iustarid,name,24,unit,16,spval)
#endif
      endif
      if (enable_ktmix) then
#ifdef write_ktmix_dml
       name = 'Mixed layer depth due to wind stirring'; unit = 'm'
       call dvcdf(ncid,idmlid,name,24,unit,16,spval)
#endif
      endif
      if (enable_expl_convection) then
#ifdef write_ventilated_depth
       name = 'Ventilated depth'; unit = 'm'
       call dvcdf(ncid,ivz,name,24,unit,16,spval)
#endif
#ifdef write_convection_below_ml
       name = 'Convection below Mixel Layer'; unit = '# of lev.'
       call dvcdf(ncid,ivz2,name,24,unit,16,spval)
#endif
      endif
#ifdef write_tau
      name = 'ZONAL WIND STRESS'; unit = 'N/m^2'
      call dvcdf(ncid,itauxid,name,24,unit,16,spval)
      name = 'MERIDIONAL WIND STRESS'; unit = 'N/m^2'
      call dvcdf(ncid,itauyid,name,24,unit,16,spval)
#endif
#ifdef write_barotropic_forcing
      name = 'Zonal comp. of barotr. forcing'; unit = '*'
      call dvcdf(ncid,zu1id,name,24,unit,16,spval)
      name = 'Mer. comp. of barotr. forcing'; unit = '*'
      call dvcdf(ncid,zu2id,name,24,unit,16,spval)
      name = 'Last change in streamfct'; unit = 'Sv*dt'
      call dvcdf(ncid,guessid,name,24,unit,16,spval)
      name = 'Actual change in streamfct'; unit = 'Sv*dt'
      call dvcdf(ncid,ptdid,name,24,unit,16,spval)
#endif
#if defined write_am_scale 
      name = 'Hor. harm. viscosity'; unit = 'm^2/s'
      call dvcdf(ncid,iamid,name,24,unit,16,spval)
#endif
#if defined write_ambi_scale 
      name = 'Hor. biharm. viscosity'; unit = 'm^4/s'
      call dvcdf(ncid,iambiid,name,24,unit,16,spval)
#endif
#if defined write_ah_scale
      name = 'Hor. harm. Diffusivity'; unit = 'm^2/s'
      call dvcdf(ncid,iahid,name,24,unit,16,spval)
#endif
#if defined write_ahbi_scale
      name = 'Hor. biharm. Diffusivity'; unit = 'm^4/s'
      call dvcdf(ncid,iahbiid,name,24,unit,16,spval)
#endif
#ifdef write_restore_TS_obc
      name = 'Restoring Temp. at northern obc'; unit = 'deg C'
      call dvcdf(ncid,tobcnid,name,24,unit,16,spval)
      name = 'Restoring Sal. at northern obc'; unit = 'psu'
      call dvcdf(ncid,sobcnid,name,24,unit,16,spval)
      name = 'Restoring Temp. at southern obc'; unit = 'deg C'
      call dvcdf(ncid,tobcsid,name,24,unit,16,spval)
      name = 'Restoring Sal. at southern obc'; unit = 'psu'
      call dvcdf(ncid,sobcsid,name,24,unit,16,spval)
      name = 'Restoring Temp. at western obc'; unit = 'deg C'
      call dvcdf(ncid,tobcwid,name,24,unit,16,spval)
      name = 'Restoring Sal. at western obc'; unit = 'psu'
      call dvcdf(ncid,sobcwid,name,24,unit,16,spval)
      name = 'Restoring Temp. at eastern obc'; unit = 'deg C'
      call dvcdf(ncid,tobceid,name,24,unit,16,spval)
      name = 'Restoring Sal. at eastern obc'; unit = 'psu'
      call dvcdf(ncid,sobceid,name,24,unit,16,spval)
#endif
      call ncendf(ncid, iret)

      Corner=0; edges=0; Corner(1) = 1; edges(1) = imt
      allocate( var4(max(imt,jmt,km)) )
      var4(1:imt)=xt
      call ncvpt(ncid, Lon_tid, corner, edges,var4, iret)
      var4(1:imt)=xu
      call ncvpt(ncid, Lon_uid, corner, edges,var4, iret)
      corner(1) = 1; edges(1) = jmt
      var4(1:jmt)=yt
      call ncvpt(ncid, Lat_tid, corner, edges,var4, iret)
      var4(1:jmt)=yu
      call ncvpt(ncid, Lat_uid, corner, edges,var4, iret)
      corner(1) = 1; edges(1) = km
      var4(1:km)=zt/100.
      call ncvpt(ncid, depth_tid, corner, edges,var4, iret)
      corner(1) = 1; edges(1) = km
      var4(1:km)=zw/100.
      call ncvpt(ncid, depth_wid, corner, edges,var4, iret)
      deallocate(var4)
      call ncclos (ncid, iret)
      print*,' done'
      endif
#endif
      end subroutine init_snap_cdf



      subroutine diag_snap
c
c-----------------------------------------------------------------------
c     write to NetCDF snapshot file
c-----------------------------------------------------------------------
c
      use spflame_module
#if defined write_gm_velocity 
      use isoneutral_module
      use isopycnic_module
#endif
      use blue_module
      use freesurf_module
      implicit none
#ifdef netcdf_diagnostics
#include "netcdf.inc"
      integer ncid,iret,n, corner(4), edges(4)
      real :: spval=-9.9e12
      integer ipsiid,itdimid,ilen,itopoid,ihtpid
      integer itauxid,itauyid, uid,vid,wid,trid(nt)
#ifdef netcdf_real4
      real (kind=4) :: var(is_pe:ie_pe,js_pe:je_pe)
      real (kind=4) , allocatable :: v2(:)
#else
      real var(is_pe:ie_pe,js_pe:je_pe)
      real, allocatable :: v2(:)
#endif
      integer i,j,k,is,ie,js,je
      integer rhoid,brhoid ,zu1id,zu2id,guessid,ptdid
      integer npe,ierr,itimeid, stfid(nt),ivz,ivz2
      integer tobcnid,sobcnid,tobcsid,sobcsid
      integer tobcwid,sobcwid, tobceid,sobceid
      integer ihfrestid,isfrestid,iustarid,idmlid
      integer ihfclimid,isfclimid, balphaid
      integer iamid,iahid,iambiid,iahbiid,adv_sid,sshid,ubarid,vbarid
      integer ugmid,vgmid,wgmid,diffvid,ekeid,viscvid,Nsqrid
      character (len=80) :: name
      real, allocatable :: u_gm(:,:,:),v_gm(:,:,:)
      real, allocatable :: w_gm(:,:,:)
      real :: fxa
#ifdef write_Nsqr
      real :: Nsqr(is_pe:ie_pe,km,js_pe:je_pe)  
      real, dimension(is_pe:ie_pe)    :: ro1,ro2
#endif

       if (my_pe==0)
     &  print*,' --> writing a snapshot of sub_domain to file ',
     &        snap_file(1:len_trim(snap_file))

#if defined write_gm_velocity 
      if (enable_diffusion_isoneutral.or.
     &    enable_diffusion_isopycnic) then
       allocate( u_gm(is_pe-1:ie_pe,  km,js_pe  :je_pe) ); u_gm=0.
       allocate( v_gm(is_pe  :ie_pe,  km,js_pe-1:je_pe) ); v_gm=0.
       allocate( w_gm(is_pe  :ie_pe,0:km,js_pe  :je_pe) ); w_gm=0.
       if (enable_diffusion_isoneutral) then
        call isoneutral_adv(u_gm,v_gm,w_gm)
       else
        u_gm=adv_vetiso
        v_gm=adv_vntiso
        w_gm=adv_vbtiso
       endif
      endif
#endif
#ifdef write_Nsqr
      Nsqr=0.
      do j=js_pe,je_pe
       do k=1,km-1
        call model_dens(t(is_pe,k  ,j,1,tau),t(is_pe,k  ,j,2,tau),
     &                  ro1(is_pe),k, ie_pe-is_pe+1
#ifdef partial_cell
     &                       ,ztp(is_pe,k,j)
#endif
     &                       )
        call model_dens(t(is_pe,k+1,j,1,tau),t(is_pe,k+1,j,2,tau),
     &                  ro2(is_pe),k, ie_pe-is_pe+1
#ifdef partial_cell
     &                       ,ztp(is_pe,k,j)
#endif
     &                       )
        do i=is_pe,ie_pe
#ifndef partial_cell
         fxa = -grav/rho0*dzwr(k)
#endif
#ifdef partial_cell
         fxa = -grav/(rho0*dhwt(i,k,j))
#endif
         Nsqr(i,k,j) =fxa*(ro1(i)-ro2(i))*tmask(i,k+1,j)
        enddo
       enddo
      enddo
#endif

       do npe=0,n_pes
       call barrier

       if (my_pe==npe) then

       iret=nf_open(snap_file,NF_WRITE,ncid)
       iret=nf_set_fill(ncid, NF_NOFILL, iret)

       iret=nf_inq_varid(ncid,'psi',ipsiid)
       iret=nf_inq_varid(ncid,'topo',itopoid)
#ifdef partial_cell
       iret=nf_inq_varid(ncid,'htp',ihtpid)
#endif
#ifdef write_UV
       iret=nf_inq_varid(ncid,'u',uid)
       iret=nf_inq_varid(ncid,'v',vid)
#endif
#ifdef write_adv_vbt
       iret=nf_inq_varid(ncid,'w',wid)
#endif
#if defined write_gm_velocity 
       if (enable_diffusion_isoneutral.or.
     &     enable_diffusion_isopycnic) then
        iret=nf_inq_varid(ncid,'u_gm',ugmid)
        iret=nf_inq_varid(ncid,'v_gm',vgmid)
        iret=nf_inq_varid(ncid,'w_gm',wgmid)
       endif
#endif
#ifdef write_vert_diff
       iret=nf_inq_varid(ncid,'diff_cbt',diffvid)
#endif
#ifdef write_vert_visc
       iret=nf_inq_varid(ncid,'visc_cbu',viscvid)
#endif
#ifdef write_Nsqr
       iret=nf_inq_varid(ncid,'Nsqr',Nsqrid)
#endif
#ifdef write_tkemix
       if (enable_tkemix) then
        iret=nf_inq_varid(ncid,'tke',ekeid)
       endif
#endif
#ifdef write_TS
       do n=1,nt
        name='temp';if (n==2) name='salt'
        if (n>2) then
          write(name, '("tracer_",i2)') n
          call replace_space_zero(name)
        endif
        iret=nf_inq_varid(ncid,name,trid(n))
       enddo
#endif
       if (enable_blue) then
#if defined write_blue_density 
        iret=nf_inq_varid(ncid,'rho',rhoid)
        iret=nf_inq_varid(ncid,'blue_rho',brhoid)
#endif
#if defined write_blue_parameter
        iret=nf_inq_varid(ncid,'blue_alpha',balphaid)
#endif
       endif
#ifdef write_tau
       iret=nf_inq_varid(ncid,'taux',itauxid)
       iret=nf_inq_varid(ncid,'tauy',itauyid)
#endif
#ifdef write_tracer_fluxes
       do n=1,nt
        name='heat_fl';if (n==2) name='salt_fl'
        if (n>2) then
          write(name, '("tracer_",i2,"_sflux")') n
          call replace_space_zero(name)
        endif
        iret=nf_inq_varid(ncid,name,stfid(n))
       enddo
#endif
#ifdef write_surface_restoring
       iret=nf_inq_varid(ncid,'T_rest',ihfrestid)
       iret=nf_inq_varid(ncid,'S_rest',isfrestid)
       iret=nf_inq_varid(ncid,'T_clim',ihfclimid)
       iret=nf_inq_varid(ncid,'S_clim',isfclimid)
#endif
#if defined write_eta 
       if (enable_freesurf) then
        iret=nf_inq_varid(ncid,'eta',sshid)
       endif
#endif
#if defined write_ubar
       if (enable_freesurf) then
        iret=nf_inq_varid(ncid,'ubar',ubarid)
        iret=nf_inq_varid(ncid,'vbar',vbarid)
       endif
#endif
       if (enable_ktmix.or.enable_tkemix) then
#if defined write_ustar 
        iret=nf_inq_varid(ncid,'ustar',iustarid)
#endif
       endif
       if (enable_ktmix) then
#ifdef write_ktmix_dml
        iret=nf_inq_varid(ncid,'dml',idmlid)
#endif
       endif
       if (enable_expl_convection) then
#ifdef write_ventilated_depth
        iret=nf_inq_varid(ncid,'convect',ivz)
#endif
#ifdef write_convection_below_ml
        iret=nf_inq_varid(ncid,'convect2',ivz2)
#endif
       endif
#ifdef write_barotropic_forcing
       iret=nf_inq_varid(ncid,'zu1',zu1id)
       iret=nf_inq_varid(ncid,'zu2',zu2id)
       iret=nf_inq_varid(ncid,'guess',guessid)
       iret=nf_inq_varid(ncid,'ptd',ptdid)
#endif
#if defined write_am_scale 
       iret=nf_inq_varid(ncid,'am',iamid)
#endif
#if defined write_ambi_scale 
       iret=nf_inq_varid(ncid,'ambi',iambiid)
#endif
#if defined write_ah_scale 
       iret=nf_inq_varid(ncid,'ah',iahid)
#endif
#if defined write_ahbi_scale 
       iret=nf_inq_varid(ncid,'ahbi',iahbiid)
#endif

       iret=nf_inq_varid(ncid,'Time',itimeid)
#ifdef write_restore_TS_obc
       iret=nf_inq_varid(ncid,'t_obc_north',tobcnid)
       iret=nf_inq_varid(ncid,'s_obc_north',sobcnid)
       iret=nf_inq_varid(ncid,'t_obc_south',tobcsid)
       iret=nf_inq_varid(ncid,'s_obc_south',sobcsid)
       iret=nf_inq_varid(ncid,'t_obc_west',tobcwid)
       iret=nf_inq_varid(ncid,'s_obc_west',sobcwid)
       iret=nf_inq_varid(ncid,'t_obc_east',tobceid)
       iret=nf_inq_varid(ncid,'s_obc_east',sobceid)
#endif

       iret=nf_inq_dimid(ncid,'Time',itdimid)
       iret=nf_inq_dimlen(ncid, itdimid,ilen)

       if (my_pe==0) then
         ilen=ilen+1
         corner(1)=ilen
         edges(1)=1
         call read_stamp(current_stamp,fxa)
         print*,' at stamp=',current_stamp,
     &          ' (days since origin : ',fxa,')',
     &          ' (time steps in file : ',ilen,')'
         var(is_pe,js_pe)=fxa
         iret= nf_put_vara_real (ncid,itimeid , corner, edges,var)
       endif

       Corner = (/is_pe,js_pe,ilen,0/)
       edges  = (/ie_pe-is_pe+1,je_pe-js_pe+1,1,0/)
       var=psi(is_pe:ie_pe,js_pe:je_pe,1)/1e12
       where( tmask(is_pe:ie_pe,1,js_pe:je_pe) == 0. ) var = spval
       iret= nf_put_vara_real (ncid,ipsiid , corner, edges,var)

#ifdef write_tau
       var=smf(is_pe:ie_pe,js_pe:je_pe,1)/10.
       where( umask(is_pe:ie_pe,1,js_pe:je_pe) == 0. ) var = spval
       iret= nf_put_vara_real (ncid,itauxid , corner, edges,var)

       var=smf(is_pe:ie_pe,js_pe:je_pe,2)/10.
       where( umask(is_pe:ie_pe,1,js_pe:je_pe) == 0. ) var = spval
       iret= nf_put_vara_real (ncid,itauyid , corner, edges,var)
#endif

       is=max(2,is_pe); ie=min(imt-1,ie_pe)
       js=max(2,js_pe); je=min(jmt-1,je_pe)

#ifdef write_tracer_fluxes
       do n=1,nt
        var=spval
        var(is:ie,js:je)=stf(is:ie,js:je,n)
        if (n==1) var(is:ie,js:je)=var(is:ie,js:je)*41868
        if (n==2) var(is:ie,js:je)=var(is:ie,js:je)*10
        where( tmask(is:ie,1,js:je) == 0. ) var(is:ie,js:je)= spval
        iret= nf_put_vara_real (ncid,stfid(n) , corner, edges,var)
       enddo
#endif
#ifdef write_surface_restoring
       var=spval
       var(is:ie,js:je)=stf_rest(is:ie,js:je,1)*41868.  
       where( tmask(is:ie,1,js:je) == 0. ) var(is:ie,js:je)= spval
       iret= nf_put_vara_real (ncid,ihfrestid , corner, edges,var)

       var=spval
       var(is:ie,js:je)=stf_rest(is:ie,js:je,2)   ! do not known what to do
       where( tmask(is:ie,1,js:je) == 0. ) var(is:ie,js:je)= spval
       iret= nf_put_vara_real (ncid,isfrestid , corner, edges,var)

       var=spval
       var(is:ie,js:je)=stf_clim(is:ie,js:je,1)
       where( tmask(is:ie,1,js:je) == 0. ) var(is:ie,js:je)= spval
       iret= nf_put_vara_real (ncid,ihfclimid , corner, edges,var)

       var=spval
       var(is:ie,js:je)=stf_clim(is:ie,js:je,2)*1000.+35.
       where( tmask(is:ie,1,js:je) == 0. ) var(is:ie,js:je)= spval
       iret= nf_put_vara_real (ncid,isfclimid , corner, edges,var)
#endif
#if defined  write_eta 
       if (enable_freesurf) then
        var=spval
        var(is:ie,js:je)=etat(is:ie,js:je,3)/100.
        where( tmask(is:ie,1,js:je) == 0. ) var(is:ie,js:je)= spval
        iret= nf_put_vara_real (ncid,sshid , corner, edges,var)
       endif
#endif
#if defined  write_ubar
       if (enable_freesurf) then
        var=spval
        var(is:ie,js:je)=ubar(is:ie,js:je,1)/100.
        where( umask(is:ie,1,js:je) == 0. ) var(is:ie,js:je)= spval
        iret= nf_put_vara_real (ncid,ubarid , corner, edges,var)
        var(is:ie,js:je)=ubar(is:ie,js:je,2)/100.
        where( umask(is:ie,1,js:je) == 0. ) var(is:ie,js:je)= spval
        iret= nf_put_vara_real (ncid,vbarid , corner, edges,var)
       endif
#endif


       if (enable_ktmix.or.enable_tkemix) then
#if defined write_ustar 
        var=spval; var(is:ie,js:je)=ustar(is:ie,js:je)
        where( tmask(is:ie,1,js:je) == 0. ) var(is:ie,js:je)= spval
        iret= nf_put_vara_real (ncid,iustarid , corner, edges,var)
#endif
       endif
       if (enable_ktmix) then
#ifdef write_ktmix_dml
        var=spval; var(is:ie,js:je)=dml(is:ie,js:je)/100.
        where( tmask(is:ie,1,js:je) == 0. ) var(is:ie,js:je)= spval
        iret= nf_put_vara_real (ncid,idmlid , corner, edges,var)
#endif
       endif
       if (enable_expl_convection) then
#ifdef write_ventilated_depth
        var=spval; 
        do j=js,je
         do i=is,ie
          var(i,j)=zt( max(1, min(km,lcven(i,j)) )  )/100.
         enddo
        enddo
        iret= nf_put_vara_real (ncid,ivz, corner, edges,var)
#endif
#ifdef write_convection_below_ml
        var(is:ie,js:je)=lctot(is:ie,js:je)-lcven(is:ie,js:je)
        where( tmask(is_pe:ie_pe,1,js_pe:je_pe) == 0. ) var= spval
        iret= nf_put_vara_real (ncid,ivz2 , corner, edges,var)
#endif
       endif

#ifdef write_barotropic_forcing
       var=zu(is_pe:ie_pe,js_pe:je_pe,1)
c       where( umask(is_pe:ie_pe,1,js_pe:je_pe) == 0. ) var = spval
       iret= nf_put_vara_real (ncid,zu1id , corner, edges,var)
       var=zu(is_pe:ie_pe,js_pe:je_pe,2)
c       where( umask(is_pe:ie_pe,1,js_pe:je_pe) == 0. ) var = spval
       iret= nf_put_vara_real (ncid,zu2id , corner, edges,var)
       var=guess(is_pe:ie_pe,js_pe:je_pe)
       where( tmask(is_pe:ie_pe,1,js_pe:je_pe) == 0. ) var = spval
       iret= nf_put_vara_real (ncid,guessid , corner, edges,var)
       var=ptd(is_pe:ie_pe,js_pe:je_pe)
       where( tmask(is_pe:ie_pe,1,js_pe:je_pe) == 0. ) var = spval
       iret= nf_put_vara_real (ncid,ptdid , corner, edges,var)
#endif

       if (ilen==1) then
        var=0.
        do j=js_pe,je_pe 
         do k=1,km
          do i=is_pe,ie_pe
           var(i,j)=var(i,j)+tmask(i,k,j)*dzt(k)/100.
          enddo
         enddo
        enddo
        iret= nf_put_vara_real(ncid,itopoid,corner, edges,var)
        if (iret.ne.0) print*,nf_strerror(iret),my_pe
#ifdef partial_cell
        var=htp(is_pe:ie_pe,js_pe:je_pe)/100.
        iret= nf_put_vara_real(ncid,ihtpid,corner, edges,var)
#endif
#if defined write_am_scale 
        var=am_scale(is_pe:ie_pe,js_pe:je_pe)*am/1e4
        where( umask(is_pe:ie_pe,1,js_pe:je_pe) == 0. ) var = spval
        iret= nf_put_vara_real (ncid,iamid , corner, edges,var)
#endif
#if defined write_ambi_scale 
        var=ambi_scale(is_pe:ie_pe,js_pe:je_pe)*ambi/1e8
        where( umask(is_pe:ie_pe,1,js_pe:je_pe) == 0. ) var = spval
        iret= nf_put_vara_real (ncid,iambiid , corner, edges,var)
#endif
#if defined write_ah_scale 
        var=ah_scale(is_pe:ie_pe,js_pe:je_pe)*ah/1e4
        where( tmask(is_pe:ie_pe,1,js_pe:je_pe) == 0. ) var = spval
        iret= nf_put_vara_real (ncid,iahid , corner, edges,var)
#endif
#if defined write_ahbi_scale 
        var=ahbi_scale(is_pe:ie_pe,js_pe:je_pe)*ahbi/1e8
        where( tmask(is_pe:ie_pe,1,js_pe:je_pe) == 0. ) var = spval
        iret= nf_put_vara_real (ncid,iahbiid , corner, edges,var)
#endif
#if defined write_blue_parameter
        if (enable_blue) then
         do k=1,km
          Corner = (/is_pe,js_pe,k,ilen/)
          edges  = (/ie_pe-is_pe+1,je_pe-js_pe+1,1,1/)
          var=blue_alpha(is_pe:ie_pe,k,js_pe:je_pe)
          where( tmask(is_pe:ie_pe,k,js_pe:je_pe) == 0. ) var = spval
          iret= nf_put_vara_real (ncid,balphaid,corner,edges,var)
         enddo
        endif
#endif
       endif

       do k=1,km
        Corner = (/is_pe,js_pe,k,ilen/)
        edges  = (/ie_pe-is_pe+1,je_pe-js_pe+1,1,1/)
#ifdef write_UV
        var=u(is_pe:ie_pe,k,js_pe:je_pe,1,tau)/100.
        where( umask(is_pe:ie_pe,k,js_pe:je_pe) == 0. ) var = spval
        iret= nf_put_vara_real (ncid,uid , corner, edges,var)
        var=u(is_pe:ie_pe,k,js_pe:je_pe,2,tau)/100.
        where( umask(is_pe:ie_pe,k,js_pe:je_pe) == 0. ) var = spval
        iret= nf_put_vara_real (ncid,vid , corner, edges,var)
#endif
#ifdef write_TS
        do n=1,nt
         var=t(is_pe:ie_pe,k,js_pe:je_pe,n,tau)
         if (n==2) var=var*1000+35.
         where( tmask(is_pe:ie_pe,k,js_pe:je_pe) == 0. ) var = spval
         iret= nf_put_vara_real (ncid,trid(n) , corner, edges,var)
        enddo
#endif
#if defined write_blue_density 
        if (enable_blue) then
         var=rho(is_pe:ie_pe,k,js_pe:je_pe)*1000.
         where( tmask(is_pe:ie_pe,k,js_pe:je_pe) == 0. ) var = spval
         iret= nf_put_vara_real (ncid,rhoid , corner, edges,var)
         var = (blue_rho(is_pe:ie_pe,k,js_pe:je_pe,blue_p1)*blue_f1+
     &          blue_rho(is_pe:ie_pe,k,js_pe:je_pe,blue_p2)*blue_f2)
     &         *1000.
         where( tmask(is_pe:ie_pe,k,js_pe:je_pe) == 0. ) var = spval
         iret= nf_put_vara_real (ncid,brhoid , corner, edges,var)
        endif
#endif

#ifdef write_adv_vbt
        var=spval
        var(is:ie,js:je)=adv_vbt(is:ie,k,js:je)/100.
        where( tmask(is_pe:ie_pe,k,js_pe:je_pe) == 0. ) var = spval
        iret= nf_put_vara_real (ncid,wid , corner, edges,var)
#endif
#ifdef write_vert_diff
        var=spval
        var(is:ie,js:je)=diff_cbt(is:ie,k,js:je)
        where( tmask(is_pe:ie_pe,k,js_pe:je_pe) == 0. ) var = spval
        iret= nf_put_vara_real (ncid,diffvid , corner, edges,var)
#endif
#ifdef write_vert_visc
        var=spval
        var(is:ie,js:je)=visc_cbu(is:ie,k,js:je)
        where( umask(is_pe:ie_pe,k,js_pe:je_pe) == 0. ) var = spval
        iret= nf_put_vara_real (ncid,viscvid , corner, edges,var)
#endif
#ifdef write_Nsqr
        var=spval
        var(is:ie,js:je)=Nsqr(is:ie,k,js:je)
        where( tmask(is_pe:ie_pe,k,js_pe:je_pe) == 0. ) var = spval
        iret= nf_put_vara_real (ncid,Nsqrid , corner, edges,var)
#endif
#ifdef write_tkemix
        if (enable_tkemix) then
         var=spval
         var(is:ie,js:je)=eke(is:ie,k,js:je,tau)/1.e4
         where( tmask(is_pe:ie_pe,k,js_pe:je_pe) == 0. ) var = spval
         iret= nf_put_vara_real (ncid,ekeid , corner, edges,var)
        endif
#endif
       enddo


#ifdef write_restore_TS_obc
       do k=1,km
        Corner = (/is_pe,k,ilen,1/)
        edges  = (/ie_pe-is_pe+1,1,1,0/)

        allocate(v2(is_pe:ie_pe) )

        if (my_blk_j == n_pes_j ) then
         v2=spval
         if (enable_obc_north .and. restore_TS_obc_north) 
     *       v2=ts_obc_north(is_pe:ie_pe,k,1)
         where( tmask(is_pe:ie_pe,k,je_pe) == 0. ) v2 = spval
         iret= nf_put_vara_real (ncid,tobcnid , corner, edges,v2)

         v2=spval
         if (enable_obc_north .and. restore_TS_obc_north) 
     &      v2=ts_obc_north(is_pe:ie_pe,k,2)*1000.+35.
         where( tmask(is_pe:ie_pe,k,je_pe) == 0. ) v2 = spval
         iret= nf_put_vara_real (ncid,sobcnid , corner, edges,v2)
        endif

        if (my_blk_j == 1) then
         v2=spval
         if (enable_obc_south .and. restore_TS_obc_south) 
     &         v2=ts_obc_south(is_pe:ie_pe,k,1)
         where( tmask(is_pe:ie_pe,k,js_pe) == 0. ) v2 = spval
         iret= nf_put_vara_real (ncid,tobcsid , corner, edges,v2)

         v2=spval
         if (enable_obc_south .and. restore_TS_obc_south) 
     &         v2=ts_obc_south(is_pe:ie_pe,k,2)*1000.+35.
         where( tmask(is_pe:ie_pe,k,js_pe) == 0. ) v2 = spval
         iret= nf_put_vara_real (ncid,sobcsid , corner, edges,v2)
        endif

        deallocate( v2 )

        Corner = (/js_pe,k,ilen,0/)
        edges  = (/je_pe-js_pe+1,1,1,0/)

        allocate(v2(js_pe:je_pe) )

        if (my_blk_i == n_pes_i) then
         v2=spval
         if (enable_obc_east .and. restore_TS_obc_east) 
     &       v2=ts_obc_east(js_pe:je_pe,k,1)
         where( tmask(ie_pe,k,js_pe:je_pe) == 0. ) v2 = spval
         iret= nf_put_vara_real (ncid,tobceid , corner, edges,v2)

         v2=spval
         if (enable_obc_east .and. restore_TS_obc_east) 
     &       v2=ts_obc_east(js_pe:je_pe,k,2)*1000.+35.
         where( tmask(ie_pe,k,js_pe:je_pe) == 0. ) v2 = spval
         iret= nf_put_vara_real (ncid,sobceid , corner, edges,v2)

        endif

        if (my_blk_i == 1) then
         v2=spval
         if (enable_obc_west .and. restore_TS_obc_west) 
     &       v2=ts_obc_west(js_pe:je_pe,k,1)
         where( tmask(is_pe,k,js_pe:je_pe) == 0. ) v2 = spval
         iret= nf_put_vara_real (ncid,tobcwid , corner, edges,v2)
c
         v2=spval
         if (enable_obc_west .and. restore_TS_obc_west) 
     &       v2=ts_obc_west(js_pe:je_pe,k,2)*1000.+35.
         where( tmask(is_pe,k,js_pe:je_pe) == 0. ) v2 = spval
         iret= nf_put_vara_real (ncid,sobcwid , corner, edges,v2)

        endif

        deallocate(v2)

       enddo
#endif

#if defined write_gm_velocity 
       if (enable_diffusion_isoneutral.or.
     &     enable_diffusion_isopycnic) then
        do k=1,km
         Corner = (/is_pe,js_pe,k,ilen/)
         edges  = (/ie_pe-is_pe+1,je_pe-js_pe+1,1,1/)
         do j=js_pe,je_pe
          do i=is_pe,ie_pe
           var(i,j) = u_gm(i,k,j)/100.
          enddo
         enddo
         where( tmask(is_pe:ie_pe,k,js_pe:je_pe) == 0. ) var = spval
         iret= nf_put_vara_real (ncid,ugmid , corner, edges,var)
         do j=js_pe,je_pe
          do i=is_pe,ie_pe 
           var(i,j) = v_gm(i,k,j)/csu(j)/100.
          enddo
         enddo
         where( tmask(is_pe:ie_pe,k,js_pe:je_pe) == 0. ) var = spval
         iret= nf_put_vara_real (ncid,vgmid , corner, edges,var)
         do j=js_pe,je_pe
          do i=is_pe,ie_pe
           var(i,j) = w_gm(i,k,j)   /100.
          enddo
         enddo
         where( tmask(is_pe:ie_pe,k,js_pe:je_pe) == 0. ) var = spval
         iret= nf_put_vara_real (ncid,wgmid , corner, edges,var)
        enddo
       endif
#endif

       call ncclos (ncid, iret)
       endif
       call barrier
      enddo

#if defined write_gm_velocity 
      if (enable_diffusion_isoneutral.or.
     &    enable_diffusion_isopycnic) then
       deallocate( u_gm,v_gm,w_gm )
      endif
#endif
#else
c
c     write to binary files here
c

#endif
      end subroutine diag_snap

