#include "options.inc"

c
c-----------------------------------------------------------------------
c     Time series diagnostics.
c     Linked in the code in diag.F 
c     Activated by namelist parameter enable_timeseries
c     Binary data can be converted to NetCDF with ../cv_cdf/cv_anna.F
c-----------------------------------------------------------------------
c 
c ===== local defines ======
c
c The following defines can be either switched on or off
c to specify the data to be written out
c
c
c write u,v,t,s and tracer completely
c-------------------------
#define write_uvts_complete
c
c mixed layer depth 
c-------------------------
#define write_dml
c
c temperature in 100m depth
c-------------------------
c#define write_t100
c
c U,V in 100m depth
c-------------------------
c#define write_uv100
c
c T,S,U,V in 50m depth
c-------------------------
c#define write_tsuv50
c
c T,S,U,V in 0m depth
c-------------------------
c#define write_tsuv0
c
c u/viscosity at Eq.
c -------------------------
c#define write_u0N
c#define write_Av0N
c
c u/v/w at 35W
c -------------------------
c#define write_u35W
c#define write_v35W
c#define write_w35W
c
c zonal velocity at 5W 
c -------------------------
c#define write_u5W
c
c section at 57N
c -------------------------
c#define write_T57N
c
c pot. vorticity (times f)
c -------------------------
c#define write_pot_vort
c#define write_pot_vort_35W
c#define write_pot_vort_0N
c
c Tracer at BATS station (65W, 32N)
c -------------------------
c#define write_BATS
c
c ====== end of local defines ======



      module timeserie_module
      real :: spval=-9.9e12,tt
      real,allocatable :: dmix(:,:),potv(:,:,:),bvort(:,:,:)
      integer :: k_t100,k_t50,k_t10,k_t0
      integer :: j_u0N,j_t0N
      integer :: i_u35W,i_t35W
      integer :: i_u5W
      integer :: j_t57N
      integer :: i_bats,j_bats
      integer :: nr_var
      character (len=80) :: outfile
      end module timeserie_module


      subroutine init_diag_timeserie
      use spflame_module
      use timeserie_module
      implicit none
#ifdef netcdf_diagnostics
#include "netcdf.inc"
      integer :: ncid,iret,i,j,k,n
      integer :: lat_tdim,itimedim, lon_tdim,lat_udim,lon_udim
      integer :: depth_wdim,depth_tdim,lat_uid,lon_uid
      integer :: depth_wid,depth_tid, lat_tid,itimeid,lon_tid
      integer :: dims(4), corner(4), edges(4)
      character :: name*24, unit*16, text*80
#ifdef netcdf_real4
      real (kind=4), allocatable :: v4(:)
#else
      real         , allocatable :: v4(:)
#endif
#endif
      integer :: indp,year,month,day,io,varid

      year  = get_year(get_current_time())
      month = get_month_of_year(get_current_time())
      day   = get_day_of_month(get_current_time())

#ifdef netcdf_diagnostics
      write(outfile,'("anna_",i3,"_y",i4,"m",i2,"d",i2,"_",i8,".cdf")') 
#else
      write(outfile,'("anna_",i3,"_y",i4,"m",i2,"d",i2,"_",i8,".dta")') 
#endif
     &            sub_domain,year,month,day,itt

      call replace_space_zero(outfile)

      if (my_pe == 0) then
        print*,''
        print*,' Setting up time series diagnostics:'
        print*,'  writing to file ',
     &         outfile(1:len_trim(outfile)),
#ifdef netcdf_diagnostics
     &         '  (netcdf format)'
#else
     &         '  (binary output)'
#endif
      endif

      nr_var=0

      k_t100 = indp(100.0*100.0,zt,km)
      k_t50  = indp(50.0*100.0,zt,km)
      k_t10  = indp(10.0*100.0,zt,km)
      k_t0   = indp(0.0*100.0,zt,km)

      i_t35W  = indp(360.0-35.0,xt,imt)
      i_u35W  = indp(360.0-35.0,xu,imt)
      i_u5W  = indp(360.0-5.0,xu,imt)
      j_t0N   = indp(0.0,yt,jmt)
      j_u0N  = indp(0.0,yu,jmt)
      j_t57N  = indp(0.0,yt,jmt)

#ifdef write_uvts_complete
      if (my_pe == 0) print*,'  writing u,v,t,s,... completely'
      nr_var=nr_var+1
#endif

#ifdef write_dml
      if (my_pe == 0) print*,'  writing mixed layer depth'
      allocate(dmix(is_pe:ie_pe,js_pe:je_pe) )
      dmix=0.
      nr_var=nr_var+1
#endif

#if defined write_pot_vort 
      if (my_pe == 0) print*,'  writing potential vorticity'
      nr_var=nr_var+2
#endif
#if defined write_pot_vort_35W
      if (my_pe == 0) print*,'  writing potential vorticity at 35W'
      nr_var=nr_var+2
#endif
#if defined write_pot_vort_0N
      if (my_pe == 0) print*,'  writing potential vorticity at 0N'
      nr_var=nr_var+2
#endif
#if defined write_pot_vort || defined write_pot_vort_35W || defined write_pot_vort_0N 
      allocate(potv(is_pe:ie_pe,km,js_pe:je_pe) )
      allocate(bvort(is_pe:ie_pe,km,js_pe:je_pe) )
      potv=0.0; bvort=0.0
#endif

#ifdef write_t100
      if (my_pe == 0) print*,'  writing temperature in ',
     &   zt(k_t100)/100.0,'m depth'
      nr_var=nr_var+1
#endif

#ifdef write_uv100
      if (my_pe == 0) print*,'  writing velocity in ',
     &   zt(k_t100)/100.0,'m depth'
      nr_var=nr_var+2
#endif

#ifdef write_tsuv50
      if (my_pe == 0) print*,'  writing velocity,T,S in ',
     &   zt(k_t50)/100.0,'m depth'
      nr_var=nr_var+4
#endif

#ifdef write_tsuv0
      if (my_pe == 0) print*,'  writing velocity,T,S in ',
     &   zt(k_t0)/100.0,'m depth'
      nr_var=nr_var+4
#endif

#ifdef write_u0N
      if (my_pe == 0) print*,'  writing zonal velocity at ',
     &   yu(j_u0N),' N'
      nr_var=nr_var+1
#endif

#ifdef write_Av0N
      if (my_pe == 0) print*,'  writing vert. viscosity at ',
     &   yu(j_u0N),' N'
      nr_var=nr_var+1
#endif

#ifdef write_T57N
      if (my_pe == 0) print*,'  writing temperature at ',
     &   yt(j_t57N),' N'
      nr_var=nr_var+1
#endif

#ifdef write_u35W
      if (my_pe == 0) print*,'  writing zonal velocity at ',
     &   xu(i_u35W),' W'
      nr_var=nr_var+1
#endif

#ifdef write_v35W
      if (my_pe == 0) print*,'  writing merid. velocity at ',
     &   xu(i_u35W),' W'
      nr_var=nr_var+1
#endif

#ifdef write_w35W
      if (my_pe == 0) print*,'  writing vert. velocity at ',
     &   xu(i_u35W),' W'
      nr_var=nr_var+1
#endif

#ifdef write_u5W
      if (my_pe == 0) print*,'  writing zonal velocity at ',
     &   xu(i_u5W),' W'
      nr_var=nr_var+1
#endif

#ifdef write_BATS
      i_bats  = indp(360.0-65.0,xt,imt)
      j_bats  = indp(32.0,yt,jmt)
      if (my_pe == 0) print*,'  writing tracer at station BATS '
      nr_var=nr_var+1
#endif

      if (my_pe==0) print*,'  writing altogether ',nr_var,' variables'

#ifdef netcdf_diagnostics
      if (my_pe == 0) then
      print*,'  initializing NetCDF output file ',
     &         outfile(1:len_trim(outfile))
      ncid = nccre (outfile, NCCLOB, iret)
      iret=nf_set_fill(ncid, NF_NOFILL, iret)
      call store_info_cdf(ncid)
c     dimensions
      Lat_tdim  = ncddef(ncid, 'Latitude_t',  jmt, iret)
      Lon_udim  = ncddef(ncid, 'Longitude_u', imt, iret)
      lon_tdim  = ncddef(ncid, 'Longitude_t', imt, 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)  = 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)
      dims(1)  = iTimedim
      iTimeid   = ncvdef(ncid,'Time',       NCFLOAT,1,dims,iret)

c     attributes of the grid
      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 = '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 = '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
c     define variables
c
#ifdef write_uvts_complete
      dims=(/lon_udim,lat_udim,depth_tdim,itimedim/)
      varid=ncvdef (ncid,'u', NCFLOAT,4,dims,iret)
      name = 'Zonal velocity'; unit = 'm/s'
      call dvcdf(ncid,varid,name,len_trim(name),
     &              unit,len_trim(unit),spval)
      varid=ncvdef (ncid,'v', NCFLOAT,4,dims,iret)
      name = 'Meridional velocity'; unit = 'm/s'
      call dvcdf(ncid,varid,name,len_trim(name),
     &              unit,len_trim(unit),spval)
      dims=(/lon_tdim,lat_tdim,depth_tdim,itimedim/)
      varid=ncvdef (ncid,'temp', NCFLOAT,4,dims,iret)
      name = 'Temperature'; unit = 'deg C'
      call dvcdf(ncid,varid,name,len_trim(name),
     &              unit,len_trim(unit),spval)
      varid=ncvdef (ncid,'salt', NCFLOAT,4,dims,iret)
      name = 'Salinity'; unit = 'psu'
      call dvcdf(ncid,varid,name,len_trim(name),
     &              unit,len_trim(unit),spval)
      if (nt>2) then
        do n=3,nt
          write(name,'("tracer_",i2)') n
          call replace_space_zero(name)
          varid=ncvdef (ncid,name, NCFLOAT,4,dims,iret)
          name = 'Passive tracer'; unit = ' '
          call dvcdf(ncid,varid,name,len_trim(name),
     &               unit,len_trim(unit),spval)
        enddo
      endif
#endif

#ifdef write_dml
      dims=(/lon_tdim,lat_tdim,itimedim,1/)
      varid=ncvdef (ncid,'dml', NCFLOAT,3,dims,iret)
      name = 'Mixed layer depth'; unit = 'm'
      call dvcdf(ncid,varid,name,len_trim(name),
     &              unit,len_trim(unit),spval)
#endif

#ifdef write_pot_vort
      dims=(/lon_tdim,lat_tdim,depth_tdim,itimedim/)
      varid=ncvdef (ncid,'potv', NCFLOAT,4,dims,iret)
      name = 'Potential Vorticity'; unit = '1/s^2/m'
      call dvcdf(ncid,varid,name,len_trim(name),
     &              unit,len_trim(unit),spval)
      varid=ncvdef (ncid,'bvort', NCFLOAT,4,dims,iret)
      name = 'Barotropic Vorticity'; unit = '1/s^2/m'
      call dvcdf(ncid,varid,name,len_trim(name),
     &              unit,len_trim(unit),spval)
#endif

#ifdef write_pot_vort_35W
      dims=(/lat_tdim,depth_tdim,itimedim,1/)
      varid=ncvdef (ncid,'potv35W', NCFLOAT,3,dims,iret)
      name = 'Potential Vorticity at 35W'; unit = '1/s^2/m'
      call dvcdf(ncid,varid,name,len_trim(name),
     &              unit,len_trim(unit),spval)
      varid=ncvdef (ncid,'bvort35W', NCFLOAT,3,dims,iret)
      name = 'Barotropic Vorticity at 35W'; unit = '1/s^2/m'
      call dvcdf(ncid,varid,name,len_trim(name),
     &              unit,len_trim(unit),spval)
#endif

#ifdef write_pot_vort_0N
      dims=(/lon_tdim,depth_tdim,itimedim,1/)
      varid=ncvdef (ncid,'potv0N', NCFLOAT,3,dims,iret)
      name = 'Potential Vorticity at 0N'; unit = '1/s^2/m'
      call dvcdf(ncid,varid,name,len_trim(name),
     &              unit,len_trim(unit),spval)
      varid=ncvdef (ncid,'bvort0N', NCFLOAT,3,dims,iret)
      name = 'Barotropic Vorticity at 0N'; unit = '1/s^2/m'
      call dvcdf(ncid,varid,name,len_trim(name),
     &              unit,len_trim(unit),spval)
#endif

#ifdef write_t100
      dims=(/lon_tdim,lat_tdim,itimedim,1/)
      varid =ncvdef (ncid,'t100', NCFLOAT,3,dims,iret)
      name = 'Temperature in 100m'; unit = 'deg C'
      call dvcdf(ncid,varid ,name,len_trim(name),
     &              unit,len_trim(unit),spval)
#endif

#ifdef write_uv100
      dims=(/lon_udim,lat_udim,itimedim,1/)
      varid =ncvdef (ncid,'u100', NCFLOAT,3,dims,iret)
      name = 'Zonal velocity in 100m'; unit = 'm/s'
      call dvcdf(ncid,varid ,name,len_trim(name),
     &              unit,len_trim(unit),spval)
      dims=(/lon_udim,lat_udim,itimedim,1/)
      varid =ncvdef (ncid,'v100', NCFLOAT,3,dims,iret)
      name = 'Meridional velocity in 100m'; unit = 'm/s'
      call dvcdf(ncid,varid ,name,len_trim(name),
     &              unit,len_trim(unit),spval)
#endif

#ifdef write_tsuv50
      dims=(/lon_tdim,lat_tdim,itimedim,1/)
      varid =ncvdef (ncid,'s50', NCFLOAT,3,dims,iret)
      name = 'Salinity in 50m'; unit = 'psu'
      call dvcdf(ncid,varid ,name,len_trim(name),
     &              unit,len_trim(unit),spval)
      dims=(/lon_tdim,lat_tdim,itimedim,1/)
      varid =ncvdef (ncid,'t50', NCFLOAT,3,dims,iret)
      name = 'Temperature in 50m'; unit = 'deg C'
      call dvcdf(ncid,varid ,name,len_trim(name),
     &              unit,len_trim(unit),spval)
      dims=(/lon_udim,lat_udim,itimedim,1/)
      varid =ncvdef (ncid,'u50', NCFLOAT,3,dims,iret)
      name = 'Zonal velocity in 50m'; unit = 'm/s'
      call dvcdf(ncid,varid ,name,len_trim(name),
     &              unit,len_trim(unit),spval)
      dims=(/lon_udim,lat_udim,itimedim,1/)
      varid =ncvdef (ncid,'v50', NCFLOAT,3,dims,iret)
      name = 'Meridional velocity in 50m'; unit = 'm/s'
      call dvcdf(ncid,varid ,name,len_trim(name),
     &              unit,len_trim(unit),spval)
#endif

#ifdef write_tsuv0
      dims=(/lon_tdim,lat_tdim,itimedim,1/)
      varid =ncvdef (ncid,'s0', NCFLOAT,3,dims,iret)
      name = 'Salinity in 0m'; unit = 'psu'
      call dvcdf(ncid,varid ,name,len_trim(name),
     &              unit,len_trim(unit),spval)
      dims=(/lon_tdim,lat_tdim,itimedim,1/)
      varid =ncvdef (ncid,'t0', NCFLOAT,3,dims,iret)
      name = 'Temperature in 0m'; unit = 'deg C'
      call dvcdf(ncid,varid ,name,len_trim(name),
     &              unit,len_trim(unit),spval)
      dims=(/lon_udim,lat_udim,itimedim,1/)
      varid =ncvdef (ncid,'u0', NCFLOAT,3,dims,iret)
      name = 'Zonal velocity in 0m'; unit = 'm/s'
      call dvcdf(ncid,varid ,name,len_trim(name),
     &              unit,len_trim(unit),spval)
      dims=(/lon_udim,lat_udim,itimedim,1/)
      varid =ncvdef (ncid,'v0', NCFLOAT,3,dims,iret)
      name = 'Meridional velocity in 0m'; unit = 'm/s'
      call dvcdf(ncid,varid ,name,len_trim(name),
     &              unit,len_trim(unit),spval)
#endif


#ifdef write_u0N
      dims=(/lon_udim,depth_tdim,itimedim,1/)
      varid=ncvdef (ncid,'u0N', NCFLOAT,3,dims,iret)
      name = 'Zonal velocity at Eq.'; unit = 'm/s'
      call dvcdf(ncid,varid,name,len_trim(name),
     &              unit,len_trim(unit),spval)
#endif

#ifdef write_Av0N
      dims=(/lon_udim,depth_wdim,itimedim,1/)
      varid=ncvdef (ncid,'Av0N', NCFLOAT,3,dims,iret)
      name = 'Vert. viscosity at Eq.'; unit = 'm^2/s'
      call dvcdf(ncid,varid,name,len_trim(name),
     &              unit,len_trim(unit),spval)
#endif

#ifdef write_T57N
      dims=(/lon_tdim,depth_tdim,itimedim,1/)
      varid=ncvdef (ncid,'T57N', NCFLOAT,3,dims,iret)
      name = 'Temperature at 57N'; unit = 'deg C'
      call dvcdf(ncid,varid,name,len_trim(name),
     &              unit,len_trim(unit),spval)
#endif

#ifdef write_u35W
      dims=(/lat_udim,depth_tdim,itimedim,1/)
      varid=ncvdef (ncid,'u35W', NCFLOAT,3,dims,iret)
      name = 'Zonal velocity at 35W'; unit = 'm/s'
      call dvcdf(ncid,varid,name,len_trim(name),
     &              unit,len_trim(unit),spval)
#endif

#ifdef write_v35W
      dims=(/lat_udim,depth_tdim,itimedim,1/)
      varid=ncvdef (ncid,'v35W', NCFLOAT,3,dims,iret)
      name = 'Merid. velocity at 35W'; unit = 'm/s'
      call dvcdf(ncid,varid,name,len_trim(name),
     &              unit,len_trim(unit),spval)
#endif

#ifdef write_w35W
      dims=(/lat_udim,depth_wdim,itimedim,1/)
      varid=ncvdef (ncid,'w35W', NCFLOAT,3,dims,iret)
      name = 'Vert. velocity at 35W'; unit = 'm/s'
      call dvcdf(ncid,varid,name,len_trim(name),
     &              unit,len_trim(unit),spval)
#endif

#ifdef write_u5W
      dims=(/lat_udim,depth_tdim,itimedim,1/)
      varid=ncvdef (ncid,'u5W', NCFLOAT,3,dims,iret)
      name = 'Zonal velocity at 5W'; unit = 'm/s'
      call dvcdf(ncid,varid,name,len_trim(name),
     &              unit,len_trim(unit),spval)
#endif

#ifdef write_BATS
      dims=(/depth_tdim,itimedim,1,1/)
      varid=ncvdef (ncid,'t_bats', NCFLOAT,2,dims,iret)
      name = 'Temperature at 65W, 32N'; unit = 'deg C'
      call dvcdf(ncid,varid,name,len_trim(name),
     &              unit,len_trim(unit),spval)
      varid=ncvdef (ncid,'s_bats', NCFLOAT,2,dims,iret)
      name = 'Salinity at 65W, 32N'; unit = 'psu'
      call dvcdf(ncid,varid,name,len_trim(name),
     &              unit,len_trim(unit),spval)
      do n=3,nt
       write(name,'("tr",i2,"_bats")') n
       call replace_space_zero(name)
       varid=ncvdef (ncid,name, NCFLOAT,2,dims,iret)
       write(name,'("Passive tracer ",i2)') n; unit = ' '
       call dvcdf(ncid,varid,name,len_trim(name),
     &            unit,len_trim(unit),spval)
      enddo
#endif

      call ncendf(ncid, iret)

      allocate( v4(max(imt,jmt,km)) )
      corner(1) = 1; edges(1) = jmt
      v4(1:jmt)=yt
      call ncvpt(ncid, Lat_tid, corner, edges,v4, iret)
      v4(1:jmt)=yu
      call ncvpt(ncid, Lat_uid, corner, edges,v4, iret)
      corner(1) = 1; edges(1) = imt
      v4(1:imt)=xt
      call ncvpt(ncid, lon_tid, corner, edges,v4, iret)
      v4(1:imt)=xu
      call ncvpt(ncid, Lon_uid, corner, edges,v4, iret)
      corner(1) = 1; edges(1) = km
      v4(1:km)=zt/100.
      call ncvpt(ncid, depth_tid, corner, edges,v4, iret)
      corner(1) = 1; edges(1) = km
      v4(1:km)=zw/100.
      call ncvpt(ncid, depth_wid, corner, edges,v4, iret)
      deallocate(v4)
      call ncclos (ncid, iret)
      print*,' done'
      endif ! my_pe==0
#else
c
c     write a header
c
      if (my_pe==0) then
       print*,'  initializing binary output file ',
     &         outfile(1:len_trim(outfile))
       call getunit(io,outfile,'usr ieee')
       write(io) imt,jmt,km,nt,nr_var
       write(io) xt,yt,xu,yu,zt,zw,spval
       close(io)
      endif
#endif

      if (my_pe == 0) then
        print*,' done setting up time series diagnostics'
        print*,''
      endif


      end subroutine init_diag_timeserie



      subroutine diag_timeserie
      use spflame_module
      use timeserie_module
      implicit none
      real    :: ro(is_pe:ie_pe,km)
      integer :: m1(is_pe:ie_pe,km),i,j,k
      integer :: kd(is_pe:ie_pe)
      real, parameter :: delta = 0.01/1000. ! [g/cm^3] mixed layer depth criterium
#if defined write_pot_vort || defined write_pot_vort_35W || defined write_pot_vort_0N
      real ::  fxf(is_pe-1:ie_pe+1,   js_pe-1:je_pe+1),f
      real :: sig0(is_pe-1:ie_pe+1,km,js_pe-1:je_pe+1),rx,ry
      real ::  frz(is_pe-1:ie_pe+1,km,js_pe-1:je_pe+1),uz,vz
      real ::  fuy(is_pe-1:ie_pe+1,km,js_pe-1:je_pe+1),uy
      real ::  fvx(is_pe-1:ie_pe+1,km,js_pe-1:je_pe+1),vx

#endif
      integer :: is,ie,js,je

      is=max(is_pe,2); ie=min(ie_pe,imt-1)
      js=max(2,js_pe); je=min(je_pe,jmt-1)
c
#ifdef write_dml
      dmix=0.
      do j=js,je
c      calculate sigma theta  [g/cm^3]
       do k=1,km
        call model_dens(t(is,k,j,1,tau),
     &                  t(is,k,j,2,tau),
     &                  ro(is,k),1,ie-is+1
#ifdef partial_cell
     &                       ,ztp(is,k,j)
#endif
     &                       )
       enddo
c      compute the mixed layer dmix from sigma theta (vector code here)
       m1(:,1)=1
       do k=2,km
        do i=is,ie
         m1(i,k)=1*m1(i,k-1) ! if the above is already zero do not go further
         if (ro(i,k) > ro(i,1) +delta ) m1(i,k)=0
        enddo
       enddo
       kd(:)=1
       do k=2,km
        do i=is,ie
         kd(i)=kd(i)+m1(i,k)
        enddo
       enddo
       do i=is,ie
        k = min( max(2,kd(i)),kmt(i,j) )
        dmix(i,j)=zt(k)/100.
       enddo
      enddo
#endif


#if defined write_pot_vort || defined write_pot_vort_35W || defined write_pot_vort_0N
c
c   ------------------------
c     calculate potential vorticity:
c   ------------------------
c
c     Q = 1/rho_0 ( u_z rho_y - v_z rho_x + rho_z( v_x-u_y +f) )
c
c     get potential density referenced to k=1


      do j=js,je
       do k=1,km
        call model_dens(t(is-1,k,j,1,tau),
     &                  t(is-1,k,j,2,tau),
     &                  sig0(is-1,k,j),1,ie-is+3
#ifdef partial_cell
     &                       ,ztp(is-1,k,j)
#endif
     &                       )
       enddo
      enddo
c
c     drhodz
c   ------------------------
      do j=js,je
       do k=2,km
        do i=is,ie
         frz(i,k,j) = -(sig0(i,k-1,j)-sig0(i,k,j))*dzwr(k)
     &                  *tmask(i,k,j)*tmask(i,k-1,j)
        enddo
       enddo
       frz(:,1,j)=0.0
      enddo 
c     interpolate rho_z to t grid
      do j=js,je
       do k=1,km-1
        do i=is,ie
         frz(i,k,j) = (frz(i,k,j)+frz(i,k+1,j))*0.5
     &                 *tmask(i,k,j)*tmask(i,k+1,j)
        enddo
       enddo
       frz(:,km,j)=0.0
      enddo 
c
c     dudy      
c   ------------------------
      do j=js,je
       do k=1,km
        do i=is-1,ie
         fuy(i,k,j) = (u(i,k,j,1,tau)-u(i,k,j-1,1,tau))*dytr(j)
        enddo
       enddo
      enddo 
c
c     dvdx
c   ------------------------
      do j=js-1,je
       do k=1,km
        do i=is,ie
         fvx(i,k,j) =  (u(i,k,j,2,tau)-u(i-1,k,j,2,tau))
     &        *dxtr(i)*cstr(max(1,j))
        enddo
       enddo
      enddo 
c
c     Coriolis factor
      fxf=0.0
      do j=js,je
       do i=is,ie
        fxf(i,j)=cori(i,j,1)
       enddo 
      enddo 
      call set_cyclic(fxf,1,1)
      call border_exchg(fxf,1,1)
c
c     calculate  rho_z ( v_x - u_y + f)
c   ------------------------
      do j=js,je
       do k=1,km
        do i=is,ie
         f = (fxf(i,j)+fxf(i-1,j)+fxf(i,j-1)+fxf(i-1,j-1) )*0.25
         uy= (fuy(i,k,j)+fuy(i-1,k,j))*0.5
         vx= (fvx(i,k,j)+fvx(i,k,j-1))*0.5
         potv(i,k,j)=frz(i,k,j)*( vx -uy+ f)
         bvort(i,k,j)=f*( vx -uy+ f)
        enddo 
       enddo 
      enddo 
c
c     dudz
c   ------------------------
      do j=js-1,je
       do k=2,km
        do i=is-1,ie
         frz(i,k,j) = (u(i,k-1,j,1,tau)-u(i,k,j,1,tau))
     &                 *dzwr(k)*umask(i,k,j)*umask(i,k-1,j)
        enddo
       enddo
       frz(:,1,j)=0.0
      enddo 
c     interpolate u_z to t grid vertically
      do j=js-1,je
       do k=1,km-1
        do i=is-1,ie
         frz(i,k,j) = (frz(i,k,j)+frz(i,k+1,j))*0.5
     &                 *umask(i,k,j)*umask(i,k+1,j)
        enddo
       enddo
       frz(:,km,j)=0.0
      enddo 
c
c     drhody      
c   ------------------------
      do j=js-1,je
       do k=1,km
        do i=is,ie
         fuy(i,k,j) = -(sig0(i,k,j+1)-sig0(i,k,j))*dyur(max(1,j))
        enddo
       enddo
      enddo 
c
c     add u_z rho_y to pot. vorticity
c   ------------------------
      do j=js,je
       do k=1,km
        do i=is,ie
         ry = (fuy(i,k,j)+fuy(i,k,j-1))*0.5
         uz = (frz(i  ,k,j)+frz(i  ,k,j-1)+
     &         frz(i-1,k,j)+frz(i-1,k,j-1))*0.25
         potv(i,k,j)=potv(i,k,j)+uz*ry
        enddo 
       enddo 
      enddo 
c
c     dvdz
c   ------------------------
      do j=js-1,je
       do k=2,km
        do i=is-1,ie
         frz(i,k,j) = (u(i,k-1,j,2,tau)-u(i,k,j,2,tau))
     &                 *dztr(k)*umask(i,k,j)*umask(i,k-1,j)
        enddo
       enddo
       frz(:,1,j)=0.0
      enddo 
c     interpolate v_z to t grid vertically
      do j=js-1,je
       do k=1,km-1
        do i=is-1,ie
         frz(i,k,j) = (frz(i,k,j)+frz(i,k+1,j))*0.5
     &                 *umask(i,k,j)*umask(i,k+1,j)
        enddo
       enddo
       frz(i,:,j)=0.0
      enddo 
c
c     drhodx      
c   ------------------------
      do j=js,je
       do k=1,km
        do i=is-1,ie
         fvx(i,k,j) = -(sig0(i+1,k,j)-sig0(i,k,j))
     &                    *dxtr(max(1,i))*cstr(j)
        enddo
       enddo
      enddo 
c
c     add  - v_z rho_x to pot. vorticity
c   ------------------------
      do j=js,je
       do k=1,km
        do i=is,ie
         rx = (fvx(i,k,j)+fvx(i-1,k,j))*0.5
         vz = (frz(i  ,k,j)+frz(i  ,k,j-1)+
     &         frz(i-1,k,j)+frz(i-1,k,j-1))*0.25
         f = (fxf(i,j)+fxf(i-1,j)+fxf(i,j-1)+fxf(i-1,j-1) )*0.25
c        multiply pot. vort with f
         potv(i,k,j)=(potv(i,k,j)-vz*rx)*rho0r*f
        enddo 
       enddo 
      enddo 

#endif

      call diag_timeserie_output

      end subroutine diag_timeserie




      subroutine diag_timeserie_output
      use spflame_module
      use timeserie_module
      implicit none

#ifdef netcdf_diagnostics
#include "netcdf.inc"
      integer :: npe,ncid,iret,timedim,timeid,varid
      integer :: corner(4),edges(4)
      character (len=80) :: name

#ifdef netcdf_real4
      real (kind=4), allocatable :: v4(:,:)
#else
      real         , allocatable :: v4(:,:)
#endif

#else
      integer :: io
      real         , allocatable :: buf(:,:)
      character (len=80) :: ident
#endif
      integer :: i,j,k,n,m

#ifdef netcdf_diagnostics
c
c     write to netcdf format here
c
      if (my_pe==0) print*,' --> writing time serie to file ',
     &         outfile(1:len_trim(outfile))
      do npe=0,n_pes
       call barrier
       if (my_pe==npe) then
        iret=nf_open(outfile,NF_WRITE,ncid)
        iret=nf_set_fill(ncid, NF_NOFILL, iret)
        iret=nf_inq_varid(ncid,'Time',timeid)
        iret=nf_inq_dimid(ncid,'Time',timedim)
        iret=nf_inq_dimlen(ncid, timedim,n)
        allocate( v4(1,1) )
        if (my_pe==0) then
         n=n+1
         corner(1)=n
         edges(1)=1
         call read_stamp(current_stamp,tt)
         print*,' at stamp=',current_stamp,
     &          ' (days since origin : ',tt,')',
     &          ' (time steps in file : ',n,')'
         v4(1,1)=tt
         iret= nf_put_vara_real (ncid,timeid, corner,edges,v4)
        endif
        deallocate(v4)

#ifdef write_uvts_complete
        allocate( v4(is_pe:ie_pe,js_pe:je_pe) )
        edges=(/ie_pe-is_pe+1,je_pe-js_pe+1,1,1/)
        do k=1,km
         corner=(/is_pe,js_pe,k,n/); 

         v4=u(is_pe:ie_pe,k,js_pe:je_pe,1,tau)
         where( umask(is_pe:ie_pe,k,js_pe:je_pe)==0) v4=spval
         iret=nf_inq_varid(ncid,'u',varid)
         iret= nf_put_vara_real (ncid,varid ,corner, edges,v4)

         v4=u(is_pe:ie_pe,k,js_pe:je_pe,2,tau)
         where( umask(is_pe:ie_pe,k,js_pe:je_pe)==0) v4=spval
         iret=nf_inq_varid(ncid,'v',varid)
         iret= nf_put_vara_real (ncid,varid ,corner, edges,v4)

         v4=t(is_pe:ie_pe,k,js_pe:je_pe,1,tau)
         where( tmask(is_pe:ie_pe,k,js_pe:je_pe)==0) v4=spval
         iret=nf_inq_varid(ncid,'temp',varid)
         iret= nf_put_vara_real (ncid,varid ,corner, edges,v4)

         v4=t(is_pe:ie_pe,k,js_pe:je_pe,2,tau)
         where( tmask(is_pe:ie_pe,k,js_pe:je_pe)==0) v4=spval
         iret=nf_inq_varid(ncid,'salt',varid)
         iret= nf_put_vara_real (ncid,varid ,corner, edges,v4)

         if (nt>2) then 
          do m=3,nt
           v4=t(is_pe:ie_pe,k,js_pe:je_pe,m,tau)
           where( tmask(is_pe:ie_pe,k,js_pe:je_pe)==0) v4=spval
           write(name,'("tracer_",i2)') m
           call replace_space_zero(name)
           iret=nf_inq_varid(ncid,name,varid)
           iret= nf_put_vara_real (ncid,varid ,corner, edges,v4)
          enddo
         endif

        enddo
        deallocate(v4)
#endif

#ifdef write_dml
        allocate( v4(is_pe:ie_pe,js_pe:je_pe) )
        v4=dmix(is_pe:ie_pe,js_pe:je_pe)
        where( tmask(is_pe:ie_pe,1,js_pe:je_pe)==0) v4=spval
        iret=nf_inq_varid(ncid,'dml',varid)
        corner=(/is_pe,js_pe,n,1/); 
        edges=(/ie_pe-is_pe+1,je_pe-js_pe+1,1,1/)
        iret= nf_put_vara_real (ncid,varid ,corner, edges,v4)
        deallocate(v4)
#endif

#ifdef write_pot_vort
        allocate( v4(is_pe:ie_pe,js_pe:je_pe) )
        iret=nf_inq_varid(ncid,'potv',varid)
        v4=spval
        corner=(/is_pe,js_pe,1,n/); 
        edges=(/ie_pe-is_pe+1,je_pe-js_pe+1,1,1/)
        iret= nf_put_vara_real (ncid,varid ,corner, edges,v4)
        corner=(/is_pe,js_pe,km,n/); 
        iret= nf_put_vara_real (ncid,varid ,corner, edges,v4)
        do k=2,km-1
         v4=potv(is_pe:ie_pe,k,js_pe:je_pe)
         where( tmask(is_pe:ie_pe,k,js_pe:je_pe)==0)   v4=spval
         where( tmask(is_pe:ie_pe,k-1,js_pe:je_pe)==0) v4=spval
         where( tmask(is_pe:ie_pe,k+1,js_pe:je_pe)==0) v4=spval
         corner=(/is_pe,js_pe,k,n/); 
         edges=(/ie_pe-is_pe+1,je_pe-js_pe+1,1,1/)
         iret= nf_put_vara_real (ncid,varid ,corner, edges,v4)
        enddo

        iret=nf_inq_varid(ncid,'bvort',varid)
        v4=spval
        corner=(/is_pe,js_pe,1,n/); 
        edges=(/ie_pe-is_pe+1,je_pe-js_pe+1,1,1/)
        iret= nf_put_vara_real (ncid,varid ,corner, edges,v4)
        corner=(/is_pe,js_pe,km,n/); 
        iret= nf_put_vara_real (ncid,varid ,corner, edges,v4)
        do k=2,km-1
         v4=bvort(is_pe:ie_pe,k,js_pe:je_pe)
         where( tmask(is_pe:ie_pe,k,js_pe:je_pe)==0)   v4=spval
         where( tmask(is_pe:ie_pe,k-1,js_pe:je_pe)==0) v4=spval
         where( tmask(is_pe:ie_pe,k+1,js_pe:je_pe)==0) v4=spval
         corner=(/is_pe,js_pe,k,n/); 
         edges=(/ie_pe-is_pe+1,je_pe-js_pe+1,1,1/)
         iret= nf_put_vara_real (ncid,varid ,corner, edges,v4)
        enddo
        deallocate(v4)
#endif

#ifdef write_pot_vort_35W
        if (is_pe <= i_t35W .and. ie_pe >= i_t35W ) then
         allocate( v4(js_pe:je_pe,km) )
         do j=js_pe,je_pe
          do k=1,km
           v4(j,k)=potv(i_t35W,k,j)
           if (tmask(i_t35W,k,j)==0) v4(j,k)=spval
          enddo
         enddo
         iret=nf_inq_varid(ncid,'potv35W',varid)
         corner=(/js_pe,1,n,1/); 
         edges=(/je_pe-js_pe+1,km,1,1/)
         iret= nf_put_vara_real (ncid,varid ,corner, edges,v4)
         do j=js_pe,je_pe
          do k=1,km
           v4(j,k)=bvort(i_t35W,k,j)
           if (tmask(i_t35W,k,j)==0) v4(j,k)=spval
          enddo
         enddo
         iret=nf_inq_varid(ncid,'bvort35W',varid)
         corner=(/js_pe,1,n,1/); 
         edges=(/je_pe-js_pe+1,km,1,1/)
         iret= nf_put_vara_real (ncid,varid ,corner, edges,v4)
         deallocate(v4)
        endif
#endif

#ifdef write_pot_vort_0N
        if (js_pe <= j_t0N .and. je_pe >= j_t0N ) then
         allocate( v4(is_pe:ie_pe,km) )
         do k=1,km
          do i=is_pe,ie_pe
           v4(i,k)=potv(i,k,j_t0N)
           if (tmask(i,k,j_t0N)==0) v4(i,k)=spval
          enddo
         enddo
         iret=nf_inq_varid(ncid,'potv0N',varid)
         corner=(/is_pe,1,n,1/); 
         edges=(/ie_pe-is_pe+1,km,1,1/)
         iret= nf_put_vara_real (ncid,varid ,corner, edges,v4)
         do k=1,km
          do i=is_pe,ie_pe
           v4(i,k)=bvort(i,k,j_t0N)
           if (tmask(i,k,j_t0N)==0) v4(i,k)=spval
          enddo
         enddo
         iret=nf_inq_varid(ncid,'bvort0N',varid)
         corner=(/is_pe,1,n,1/); 
         edges=(/ie_pe-is_pe+1,km,1,1/)
         iret= nf_put_vara_real (ncid,varid ,corner, edges,v4)
         deallocate(v4)
        endif
#endif

#ifdef write_t100
        allocate( v4(is_pe:ie_pe,js_pe:je_pe) )
        v4=t(is_pe:ie_pe,k_t100,js_pe:je_pe,1,tau)
        where( tmask(is_pe:ie_pe,k_t100,js_pe:je_pe)==0) v4=spval
        iret=nf_inq_varid(ncid,'t100',varid )
        corner=(/is_pe,js_pe,n,1/); 
        edges=(/ie_pe-is_pe+1,je_pe-js_pe+1,1,1/)
        iret= nf_put_vara_real (ncid,varid  ,corner, edges,v4)
        deallocate(v4)
#endif

#ifdef write_uv100
        allocate( v4(is_pe:ie_pe,js_pe:je_pe) )
        v4=u(is_pe:ie_pe,k_t100,js_pe:je_pe,1,tau)
        where( umask(is_pe:ie_pe,k_t100,js_pe:je_pe)==0) v4=spval
        iret=nf_inq_varid(ncid,'u100',varid )
        corner=(/is_pe,js_pe,n,1/); 
        edges=(/ie_pe-is_pe+1,je_pe-js_pe+1,1,1/)
        iret= nf_put_vara_real (ncid,varid  ,corner, edges,v4)
        v4=u(is_pe:ie_pe,k_t100,js_pe:je_pe,2,tau)
        where( umask(is_pe:ie_pe,k_t100,js_pe:je_pe)==0) v4=spval
        iret=nf_inq_varid(ncid,'v100',varid )
        corner=(/is_pe,js_pe,n,1/); 
        edges=(/ie_pe-is_pe+1,je_pe-js_pe+1,1,1/)
        iret= nf_put_vara_real (ncid,varid  ,corner, edges,v4)
        deallocate(v4)
#endif

#ifdef write_tsuv50
        allocate( v4(is_pe:ie_pe,js_pe:je_pe) )
        v4=t(is_pe:ie_pe,k_t50,js_pe:je_pe,2,tau)
        where( tmask(is_pe:ie_pe,k_t50,js_pe:je_pe)==0) v4=spval
        iret=nf_inq_varid(ncid,'s50',varid )
        corner=(/is_pe,js_pe,n,1/); 
        edges=(/ie_pe-is_pe+1,je_pe-js_pe+1,1,1/)
        iret= nf_put_vara_real (ncid,varid  ,corner, edges,v4)
        v4=t(is_pe:ie_pe,k_t50,js_pe:je_pe,1,tau)
        where( tmask(is_pe:ie_pe,k_t50,js_pe:je_pe)==0) v4=spval
        iret=nf_inq_varid(ncid,'t50',varid )
        corner=(/is_pe,js_pe,n,1/); 
        edges=(/ie_pe-is_pe+1,je_pe-js_pe+1,1,1/)
        iret= nf_put_vara_real (ncid,varid  ,corner, edges,v4)
        v4=u(is_pe:ie_pe,k_t50,js_pe:je_pe,1,tau)
        where( umask(is_pe:ie_pe,k_t50,js_pe:je_pe)==0) v4=spval
        iret=nf_inq_varid(ncid,'u50',varid )
        corner=(/is_pe,js_pe,n,1/); 
        edges=(/ie_pe-is_pe+1,je_pe-js_pe+1,1,1/)
        iret= nf_put_vara_real (ncid,varid  ,corner, edges,v4)
        v4=u(is_pe:ie_pe,k_t50,js_pe:je_pe,2,tau)
        where( umask(is_pe:ie_pe,k_t50,js_pe:je_pe)==0) v4=spval
        iret=nf_inq_varid(ncid,'v50',varid )
        corner=(/is_pe,js_pe,n,1/); 
        edges=(/ie_pe-is_pe+1,je_pe-js_pe+1,1,1/)
        iret= nf_put_vara_real (ncid,varid  ,corner, edges,v4)
        deallocate(v4)
#endif

#ifdef write_tsuv0
        allocate( v4(is_pe:ie_pe,js_pe:je_pe) )
        v4=t(is_pe:ie_pe,k_t0,js_pe:je_pe,2,tau)
        where( tmask(is_pe:ie_pe,k_t0,js_pe:je_pe)==0) v4=spval
        iret=nf_inq_varid(ncid,'s0',varid )
        corner=(/is_pe,js_pe,n,1/); 
        edges=(/ie_pe-is_pe+1,je_pe-js_pe+1,1,1/)
        iret= nf_put_vara_real (ncid,varid  ,corner, edges,v4)
        v4=t(is_pe:ie_pe,k_t0,js_pe:je_pe,1,tau)
        where( tmask(is_pe:ie_pe,k_t0,js_pe:je_pe)==0) v4=spval
        iret=nf_inq_varid(ncid,'t0',varid )
        corner=(/is_pe,js_pe,n,1/); 
        edges=(/ie_pe-is_pe+1,je_pe-js_pe+1,1,1/)
        iret= nf_put_vara_real (ncid,varid  ,corner, edges,v4)
        v4=u(is_pe:ie_pe,k_t0,js_pe:je_pe,1,tau)
        where( umask(is_pe:ie_pe,k_t0,js_pe:je_pe)==0) v4=spval
        iret=nf_inq_varid(ncid,'u0',varid )
        corner=(/is_pe,js_pe,n,1/); 
        edges=(/ie_pe-is_pe+1,je_pe-js_pe+1,1,1/)
        iret= nf_put_vara_real (ncid,varid  ,corner, edges,v4)
        v4=u(is_pe:ie_pe,k_t0,js_pe:je_pe,2,tau)
        where( umask(is_pe:ie_pe,k_t0,js_pe:je_pe)==0) v4=spval
        iret=nf_inq_varid(ncid,'v0',varid )
        corner=(/is_pe,js_pe,n,1/); 
        edges=(/ie_pe-is_pe+1,je_pe-js_pe+1,1,1/)
        iret= nf_put_vara_real (ncid,varid  ,corner, edges,v4)
        deallocate(v4)
#endif

#ifdef write_u0N
        if (js_pe <= j_u0N .and. je_pe >= j_u0N ) then
         allocate( v4(is_pe:ie_pe,km) )
         v4=u(is_pe:ie_pe,1:km,j_u0N,1,tau)/100.0
         where( umask(is_pe:ie_pe,1:km,j_u0N)==0) v4=spval
         iret=nf_inq_varid(ncid,'u0N',varid)
         corner=(/is_pe,1,n,1/); 
         edges=(/ie_pe-is_pe+1,km,1,1/)
         iret= nf_put_vara_real (ncid,varid ,corner, edges,v4)
         deallocate(v4)
        endif
#endif

#ifdef write_Av0N
        if (js_pe <= j_u0N .and. je_pe >= j_u0N ) then
         allocate( v4(is_pe:ie_pe,km) )
         v4=visc_cbu(is_pe:ie_pe,1:km,j_u0N)/1.0e4
         where( umask(is_pe:ie_pe,1:km,j_u0N)==0) v4=spval
         iret=nf_inq_varid(ncid,'Av0N',varid)
         corner=(/is_pe,1,n,1/); 
         edges=(/ie_pe-is_pe+1,km,1,1/)
         iret= nf_put_vara_real (ncid,varid ,corner, edges,v4)
         deallocate(v4)
        endif
#endif

#ifdef write_T57N
        if (js_pe <= j_t57N .and. je_pe >= j_t57N ) then
         allocate( v4(is_pe:ie_pe,km) )
         v4=t(is_pe:ie_pe,1:km,j_t57N,1,tau)
         where( tmask(is_pe:ie_pe,1:km,j_t57N)==0) v4=spval
         iret=nf_inq_varid(ncid,'T57N',varid)
         corner=(/is_pe,1,n,1/); 
         edges=(/ie_pe-is_pe+1,km,1,1/)
         iret= nf_put_vara_real (ncid,varid ,corner, edges,v4)
         deallocate(v4)
        endif
#endif

#ifdef write_u35W
        if (is_pe <= i_u35W .and. ie_pe >= i_u35W ) then
         allocate( v4(js_pe:je_pe,km) )
         do j=js_pe,je_pe
          do k=1,km
           v4(j,k)=u(i_u35W,k,j,1,tau)/100.0
           if (umask(i_u35W,k,j)==0) v4(j,k)=spval
          enddo
         enddo
         iret=nf_inq_varid(ncid,'u35W',varid)
         corner=(/js_pe,1,n,1/); 
         edges=(/je_pe-js_pe+1,km,1,1/)
         iret= nf_put_vara_real (ncid,varid ,corner, edges,v4)
         deallocate(v4)
        endif
#endif

#ifdef write_v35W
        if (is_pe <= i_u35W .and. ie_pe >= i_u35W ) then
         allocate( v4(js_pe:je_pe,km) )
         do j=js_pe,je_pe
          do k=1,km
           v4(j,k)=u(i_u35W,k,j,2,tau)/100.0
           if (umask(i_u35W,k,j)==0) v4(j,k)=spval
          enddo
         enddo
         iret=nf_inq_varid(ncid,'v35W',varid)
         corner=(/js_pe,1,n,1/); 
         edges=(/je_pe-js_pe+1,km,1,1/)
         iret= nf_put_vara_real (ncid,varid ,corner, edges,v4)
         deallocate(v4)
        endif
#endif

#ifdef write_w35W
        if (is_pe <= i_u35W .and. ie_pe >= i_u35W ) then
         allocate( v4(js_pe:je_pe,km) )
         do j=js_pe,je_pe
          do k=1,km
           v4(j,k)=adv_vbu(i_u35W,k,j)/100.0
           if (umask(i_u35W,k,j)==0) v4(j,k)=spval
          enddo
         enddo
         iret=nf_inq_varid(ncid,'w35W',varid)
         corner=(/js_pe,1,n,1/); 
         edges=(/je_pe-js_pe+1,km,1,1/)
         iret= nf_put_vara_real (ncid,varid ,corner, edges,v4)
         deallocate(v4)
        endif
#endif

#ifdef write_u5W
        if (is_pe <= i_u5W .and. ie_pe >= i_u5W ) then
         allocate( v4(js_pe:je_pe,km) )
         do j=js_pe,je_pe
          do k=1,km
           v4(j,k)=u(i_u5W,k,j,1,tau)/100.0
           if (umask(i_u5W,k,j)==0) v4(j,k)=spval
          enddo
         enddo
         iret=nf_inq_varid(ncid,'u5W',varid)
         corner=(/js_pe,1,n,1/); 
         edges=(/je_pe-js_pe+1,km,1,1/)
         iret= nf_put_vara_real (ncid,varid ,corner, edges,v4)
         deallocate(v4)
        endif
#endif

#ifdef write_BATS
        if (is_pe <= i_bats .and. ie_pe >= i_bats .and.
     &      js_pe <= j_bats .and. je_pe >= j_bats ) then

         allocate( v4(km,1) )
         corner=(/1,n,1,1/); edges=(/km,1,1,1/)
         do n=1,nt
           if (n==1) then
            name='t_bats'
           elseif (n==2) then
            name='s_bats'
           else
            write(name,'("tr",i2,"_bats")') n
            call replace_space_zero(name)
           endif
           iret= nf_inq_varid(ncid,name,varid)
           v4(1:km,1)  = t(i_bats,1:km,j_bats,n,tau)
           if (n==2) v4=v4*1000.0+35.0
           where( tmask(i_bats,1:km,j_bats)==0) v4(1:km,1)=spval
           iret= nf_put_vara_real (ncid,varid,corner, edges,v4)
         enddo
         deallocate(v4)
        endif
#endif
        call ncclos (ncid, iret)
       endif
       call barrier
      enddo
#else
c
c     Write unformatted binary output here
c
      if (my_pe==0) then
       print*,' --> writing time serie to file ',
     &         outfile(1:len_trim(outfile))
       call getunit(io,outfile,'usa ieee')
       write (io) current_stamp, itt
      endif

#ifdef write_uvts_complete
      ident='uvts'
      if (my_pe==0) write(io) ident
      allocate( buf(imt,jmt) )
      do k=1,km
       buf(is_pe:ie_pe,js_pe:je_pe)=u(is_pe:ie_pe,k,js_pe:je_pe,1,tau)
       where(umask(is_pe:ie_pe,k,js_pe:je_pe)==0) 
     &                      buf(is_pe:ie_pe,js_pe:je_pe)=spval
       call pe0_recv_2D(buf); if (my_pe==0) write(io) buf

       buf(is_pe:ie_pe,js_pe:je_pe)=u(is_pe:ie_pe,k,js_pe:je_pe,2,tau)
       where(umask(is_pe:ie_pe,k,js_pe:je_pe)==0) 
     &                      buf(is_pe:ie_pe,js_pe:je_pe)=spval
       call pe0_recv_2D(buf); if (my_pe==0) write(io) buf

       do n=1,nt
        buf(is_pe:ie_pe,js_pe:je_pe)=t(is_pe:ie_pe,k,js_pe:je_pe,n,tau)
        where(tmask(is_pe:ie_pe,k,js_pe:je_pe)==0) 
     &                      buf(is_pe:ie_pe,js_pe:je_pe)=spval
        call pe0_recv_2D(buf); if (my_pe==0) write(io) buf
       enddo
      enddo
      deallocate(buf)
#endif

#ifdef write_dml
      ident='dml'
      if (my_pe==0) write(io) ident
      allocate( buf(imt,jmt) )
      buf(is_pe:ie_pe,js_pe:je_pe)=dmix(is_pe:ie_pe,js_pe:je_pe)
      where(tmask(is_pe:ie_pe,1,js_pe:je_pe)==0) 
     &  buf(is_pe:ie_pe,js_pe:je_pe)=spval
      call pe0_recv_2D(buf); if (my_pe==0) write(io) buf
      deallocate(buf)
#endif

#ifdef write_pot_vort
      ident='potv'
      if (my_pe==0) write(io) ident
      allocate( buf(imt,jmt) )
      do k=1,km
       buf(is_pe:ie_pe,js_pe:je_pe)=potv(is_pe:ie_pe,k,js_pe:je_pe)
       where(tmask(is_pe:ie_pe,k,js_pe:je_pe)==0) 
     &  buf(is_pe:ie_pe,js_pe:je_pe)=spval
       call pe0_recv_2D(buf); if (my_pe==0) write(io) buf
      enddo

      ident='bvor'
      if (my_pe==0) write(io) ident
      do k=1,km
       buf(is_pe:ie_pe,js_pe:je_pe)=bvort(is_pe:ie_pe,k,js_pe:je_pe)
       where(tmask(is_pe:ie_pe,k,js_pe:je_pe)==0) 
     &  buf(is_pe:ie_pe,js_pe:je_pe)=spval
       call pe0_recv_2D(buf); if (my_pe==0) write(io) buf
      enddo
      deallocate(buf)
#endif


#ifdef write_pot_vort_35W
      ident='pv35'
      if (my_pe==0) write(io) ident
      allocate( buf(jmt,km) )
      if (is_pe <= i_t35W .and. ie_pe >= i_t35W ) then
       do j=js_pe,je_pe
        do k=1,km
          buf(j,k)=potv(i_t35W,k,j)
          if (tmask(i_t35W,k,j)==0) buf(j,k)=spval
        enddo
       enddo
      endif
      call pe0_recv_merid_vec(buf,km,i_t35W)
      if (my_pe==0) write(io) buf

      ident='bv35'
      if (my_pe==0) write(io) ident
      if (is_pe <= i_t35W .and. ie_pe >= i_t35W ) then
       do j=js_pe,je_pe
        do k=1,km
          buf(j,k)=bvort(i_t35W,k,j)
          if (tmask(i_t35W,k,j)==0) buf(j,k)=spval
        enddo
       enddo
      endif
      call pe0_recv_merid_vec(buf,km,i_t35W)
      if (my_pe==0) write(io) buf

      deallocate(buf)
#endif

#ifdef write_pot_vort_0N

      ident='pv0N'
      if (my_pe==0) write(io) ident
      allocate( buf(imt,km) )
      if (js_pe <= j_t0N .and. je_pe >= j_t0N ) then
       buf(is_pe:ie_pe,:)=potv(is_pe:ie_pe,1:km,j_t0N)
       where( tmask(is_pe:ie_pe,1:km,j_t0N)==0) 
     &         buf(is_pe:ie_pe,:)=spval
      endif
      call pe0_recv_zonal_vec(buf,km,j_t0N)
      if (my_pe==0) write(io) buf

      ident='bv0N'
      if (my_pe==0) write(io) ident
      if (js_pe <= j_t0N .and. je_pe >= j_t0N ) then
       buf(is_pe:ie_pe,:)=bvort(is_pe:ie_pe,1:km,j_t0N)
       where( tmask(is_pe:ie_pe,1:km,j_t0N)==0) 
     &         buf(is_pe:ie_pe,:)=spval
      endif
      call pe0_recv_zonal_vec(buf,km,j_t0N)
      if (my_pe==0) write(io) buf
      deallocate(buf)
#endif

#ifdef write_t100
      ident='t100'
      if (my_pe==0) write(io) ident
      allocate( buf(imt,jmt) )
      buf(is_pe:ie_pe,js_pe:je_pe)=
     &  t(is_pe:ie_pe,k_t100,js_pe:je_pe,1,tau)
      where(tmask(is_pe:ie_pe,k_t100,js_pe:je_pe)==0) 
     &  buf(is_pe:ie_pe,js_pe:je_pe)=spval
      call pe0_recv_2D(buf); if (my_pe==0) write(io) buf
      deallocate(buf)
#endif

#ifdef write_uv100
      allocate( buf(imt,jmt) )
      ident='u100'
      if (my_pe==0) write(io) ident
      buf(is_pe:ie_pe,js_pe:je_pe)=
     &  u(is_pe:ie_pe,k_t100,js_pe:je_pe,1,tau)
      where(umask(is_pe:ie_pe,k_t100,js_pe:je_pe)==0) 
     &  buf(is_pe:ie_pe,js_pe:je_pe)=spval
      call pe0_recv_2D(buf); if (my_pe==0) write(io) buf
      ident='v100'
      if (my_pe==0) write(io) ident
      buf(is_pe:ie_pe,js_pe:je_pe)=
     &  u(is_pe:ie_pe,k_t100,js_pe:je_pe,2,tau)
      where(umask(is_pe:ie_pe,k_t100,js_pe:je_pe)==0) 
     &  buf(is_pe:ie_pe,js_pe:je_pe)=spval
      call pe0_recv_2D(buf); if (my_pe==0) write(io) buf
      deallocate(buf)
#endif

#ifdef write_tsuv50
      allocate( buf(imt,jmt) )
      ident='s50'
      if (my_pe==0) write(io) ident
      buf(is_pe:ie_pe,js_pe:je_pe)=
     &  t(is_pe:ie_pe,k_t50,js_pe:je_pe,2,tau)
      where(tmask(is_pe:ie_pe,k_t50,js_pe:je_pe)==0) 
     &  buf(is_pe:ie_pe,js_pe:je_pe)=spval
      call pe0_recv_2D(buf); if (my_pe==0) write(io) buf
      ident='t50'
      if (my_pe==0) write(io) ident
      buf(is_pe:ie_pe,js_pe:je_pe)=
     &  t(is_pe:ie_pe,k_t50,js_pe:je_pe,1,tau)
      where(tmask(is_pe:ie_pe,k_t50,js_pe:je_pe)==0) 
     &  buf(is_pe:ie_pe,js_pe:je_pe)=spval
      call pe0_recv_2D(buf); if (my_pe==0) write(io) buf
      ident='u50'
      if (my_pe==0) write(io) ident
      buf(is_pe:ie_pe,js_pe:je_pe)=
     &  u(is_pe:ie_pe,k_t50,js_pe:je_pe,1,tau)
      where(umask(is_pe:ie_pe,k_t50,js_pe:je_pe)==0) 
     &  buf(is_pe:ie_pe,js_pe:je_pe)=spval
      call pe0_recv_2D(buf); if (my_pe==0) write(io) buf
      ident='v50'
      if (my_pe==0) write(io) ident
      buf(is_pe:ie_pe,js_pe:je_pe)=
     &  u(is_pe:ie_pe,k_t50,js_pe:je_pe,2,tau)
      where(umask(is_pe:ie_pe,k_t50,js_pe:je_pe)==0) 
     &  buf(is_pe:ie_pe,js_pe:je_pe)=spval
      call pe0_recv_2D(buf); if (my_pe==0) write(io) buf
      deallocate(buf)
#endif

#ifdef write_tsuv0
      allocate( buf(imt,jmt) )
      ident='s0'
      if (my_pe==0) write(io) ident
      buf(is_pe:ie_pe,js_pe:je_pe)=
     &  t(is_pe:ie_pe,k_t0,js_pe:je_pe,2,tau)
      where(tmask(is_pe:ie_pe,k_t0,js_pe:je_pe)==0) 
     &  buf(is_pe:ie_pe,js_pe:je_pe)=spval
      call pe0_recv_2D(buf); if (my_pe==0) write(io) buf
      ident='t0'
      if (my_pe==0) write(io) ident
      buf(is_pe:ie_pe,js_pe:je_pe)=
     &  t(is_pe:ie_pe,k_t0,js_pe:je_pe,1,tau)
      where(tmask(is_pe:ie_pe,k_t0,js_pe:je_pe)==0) 
     &  buf(is_pe:ie_pe,js_pe:je_pe)=spval
      call pe0_recv_2D(buf); if (my_pe==0) write(io) buf
      ident='u0'
      if (my_pe==0) write(io) ident
      buf(is_pe:ie_pe,js_pe:je_pe)=
     &  u(is_pe:ie_pe,k_t0,js_pe:je_pe,1,tau)
      where(umask(is_pe:ie_pe,k_t0,js_pe:je_pe)==0) 
     &  buf(is_pe:ie_pe,js_pe:je_pe)=spval
      call pe0_recv_2D(buf); if (my_pe==0) write(io) buf
      ident='v0'
      if (my_pe==0) write(io) ident
      buf(is_pe:ie_pe,js_pe:je_pe)=
     &  u(is_pe:ie_pe,k_t0,js_pe:je_pe,2,tau)
      where(umask(is_pe:ie_pe,k_t0,js_pe:je_pe)==0) 
     &  buf(is_pe:ie_pe,js_pe:je_pe)=spval
      call pe0_recv_2D(buf); if (my_pe==0) write(io) buf
      deallocate(buf)
#endif

#ifdef write_u0N
      ident='u0N'
      if (my_pe==0) write(io) ident
      allocate( buf(imt,km) )
      if (js_pe <= j_u0N .and. je_pe >= j_u0N ) then
       buf(is_pe:ie_pe,:)=u(is_pe:ie_pe,1:km,j_u0N,1,tau)/100.0
       where( umask(is_pe:ie_pe,1:km,j_u0N)==0) 
     &         buf(is_pe:ie_pe,:)=spval
      endif
      call pe0_recv_zonal_vec(buf,km,j_u0N)
      if (my_pe==0) write(io) buf
      deallocate(buf)
#endif

#ifdef write_Av0N
      ident='Av0N'
      if (my_pe==0) write(io) ident
      allocate( buf(imt,km) )
      if (js_pe <= j_u0N .and. je_pe >= j_u0N ) then
       buf(is_pe:ie_pe,:)=visc_cbu(is_pe:ie_pe,1:km,j_u0N)/1.0e4
       where( umask(is_pe:ie_pe,1:km,j_u0N)==0) 
     &         buf(is_pe:ie_pe,:)=spval
      endif
      call pe0_recv_zonal_vec(buf,km,j_u0N)
      if (my_pe==0) write(io) buf
      deallocate(buf)
#endif

#ifdef write_T57N
      ident='T57N'
      if (my_pe==0) write(io) ident
      allocate( buf(imt,km) )
      if (js_pe <= j_t57N .and. je_pe >= j_t57N ) then
       buf(is_pe:ie_pe,:)=t(is_pe:ie_pe,1:km,j_t57N,1,tau)
       where( tmask(is_pe:ie_pe,1:km,j_t57N)==0) 
     &         buf(is_pe:ie_pe,:)=spval
      endif
      call pe0_recv_zonal_vec(buf,km,j_t57N)
      if (my_pe==0) write(io) buf
      deallocate(buf)
#endif

#ifdef write_u35W
      ident='u35W'
      if (my_pe==0) write(io) ident
      allocate( buf(jmt,km) )
      if (is_pe <= i_u35W .and. ie_pe >= i_u35W ) then
       do j=js_pe,je_pe
        do k=1,km
          buf(j,k)=u(i_u35W,k,j,1,tau)/100.0
          if (umask(i_u35W,k,j)==0) buf(j,k)=spval
        enddo
       enddo
      endif
      call pe0_recv_merid_vec(buf,km,i_u35W)
      if (my_pe==0) write(io) buf
      deallocate(buf)
#endif

#ifdef write_v35W
      ident='v35W'
      if (my_pe==0) write(io) ident
      allocate( buf(jmt,km) )
      if (is_pe <= i_u35W .and. ie_pe >= i_u35W ) then
       do j=js_pe,je_pe
        do k=1,km
          buf(j,k)=u(i_u35W,k,j,2,tau)/100.0
          if (umask(i_u35W,k,j)==0) buf(j,k)=spval
        enddo
       enddo
      endif
      call pe0_recv_merid_vec(buf,km,i_u35W)
      if (my_pe==0) write(io) buf
      deallocate(buf)
#endif

#ifdef write_w35W
      ident='w35W'
      if (my_pe==0) write(io) ident
      allocate( buf(jmt,km) )
      if (is_pe <= i_u35W .and. ie_pe >= i_u35W ) then
       do j=js_pe,je_pe
        do k=1,km
          buf(j,k)=adv_vbu(i_u35W,k,j)/100.0
          if (umask(i_u35W,k,j)==0) buf(j,k)=spval
        enddo
       enddo
      endif
      call pe0_recv_merid_vec(buf,km,i_u35W)
      if (my_pe==0) write(io) buf
      deallocate(buf)
#endif

#ifdef write_u5W
      ident='u5W'
      if (my_pe==0) write(io) ident
      allocate( buf(jmt,km) )
      if (is_pe <= i_u5W .and. ie_pe >= i_u5W ) then
       do j=js_pe,je_pe
        do k=1,km
          buf(j,k)=u(i_u5W,k,j,1,tau)/100.0
          if (umask(i_u5W,k,j)==0) buf(j,k)=spval
        enddo
       enddo
      endif
      call pe0_recv_merid_vec(buf,km,i_u5W)
      if (my_pe==0) write(io) buf
      deallocate(buf)
#endif

#ifdef write_BATS
      ident='bats'
      if (my_pe==0) write(io) ident
      allocate( buf(km,nt) )
      if (is_pe <= i_bats .and. ie_pe >= i_bats .and.
     &    js_pe <= j_bats .and. je_pe >= j_bats ) then
       do n=1,nt
        buf(:,n) = t(i_bats,1:km,j_bats,n,tau)
        where( tmask(i_bats,1:km,j_bats)==0) buf(:,n)=spval
       enddo
      endif
      call pe0_recv_depth_vec(buf,nt,i_bats,j_bats)
      if (my_pe==0) write(io) buf
      deallocate(buf)
#endif

      if (my_pe==0) close(io)
#endif

      end subroutine diag_timeserie_output



