#include "options.inc"
c
c--------------------------------------------------------
c     some netcdf output stuff
c--------------------------------------------------------
c



      subroutine init_forcing_file(forcing_file)
      use prep_module
      implicit none
#include "netcdf.inc"
      integer   ncid,iret,n
      integer lon_tdim,lat_tdim,itimedim,depth_tdim,depth_wid
      integer lon_tid ,lat_tid ,itimeid, depth_tid,depth_wdim
      integer lon_udim,lat_udim, lon_uid ,lat_uid
      integer dims(4), corner(4), edges(4),i
#ifdef netcdf_real4
      real (kind=4), allocatable :: var(:,:)
#else
      real , allocatable :: var(:,:)
#endif
      character name*80, unit*80, text*80
      character (len=*) :: forcing_file

      ncid = nccre (forcing_file , NCCLOB, iret)
      iret=nf_set_fill(ncid, NF_NOFILL, iret)
      i=imt; lon_tdim  = ncddef(ncid, 'Longitude_t', i, iret)
      i=jmt; Lat_tdim  = ncddef(ncid, 'Latitude_t',  i, iret)
      i=imt; lon_udim  = ncddef(ncid, 'Longitude_u', i, iret)
      i=jmt; Lat_udim  = ncddef(ncid, 'Latitude_u',  i, iret)
      i=km;  depth_tdim= ncddef(ncid, 'depth_t',  i, iret)
      i=km+1;depth_wdim= ncddef(ncid, 'depth_w',  i, iret)

      dims(1)  = Lon_tdim
      i=1
      Lon_tid  = ncvdef (ncid,'Longitude_t',NCFLOAT,i,dims,iret)
      dims(1)  = Lat_tdim
      Lat_tid  = ncvdef (ncid,'Latitude_t', NCFLOAT,i,dims,iret)
      dims(1)  = Lon_udim
      Lon_uid  = ncvdef (ncid,'Longitude_u',NCFLOAT,i,dims,iret)
      dims(1)  = Lat_udim
      Lat_uid  = ncvdef (ncid,'Latitude_u', NCFLOAT,i,dims,iret)
      dims(1)  = depth_tdim
      depth_tid = ncvdef (ncid,'depth_t', NCFLOAT,i,dims,iret)
      dims(1)  = depth_wdim
      depth_wid = ncvdef (ncid,'depth_w', NCFLOAT,i,dims,iret)
c
      name = 'Longitude on rotated T grid'
      unit = 'degrees_W       '
      i=len_trim(name)
      call ncaptc(ncid, Lon_tid, 'long_name', NCCHAR, i, name, iret) 
      i=len_trim(unit)
      call ncaptc(ncid, Lon_tid, 'units',     NCCHAR, i, unit, iret) 

      name = 'Longitude on rotated U grid'
      unit = 'degrees_W       '
      i=len_trim(name)
      call ncaptc(ncid, Lon_uid, 'long_name', NCCHAR, i, name, iret) 
      i=len_trim(unit)
      call ncaptc(ncid, Lon_uid, 'units',     NCCHAR, i, unit, iret) 

      name = 'Latitude on rotated T grid'
      unit = 'degrees_N       '
      i=len_trim(name)
      call ncaptc(ncid, Lat_tid, 'long_name', NCCHAR, i, name, iret) 
      i=len_trim(unit)
      call ncaptc(ncid, Lat_tid, 'units',     NCCHAR, i, unit, iret) 

      name = 'Latitude on rotated U grid'
      unit = 'degrees_N       '
      i=len_trim(name)
      call ncaptc(ncid, Lat_uid, 'long_name', NCCHAR, i, name, iret) 
      i=len_trim(unit)
      call ncaptc(ncid, Lat_uid, 'units',     NCCHAR, i, unit, iret) 

      name = 'Depth of T Grid points  '
      unit = 'm '
      i=len_trim(name)
      call ncaptc(ncid, depth_tid, 'long_name', NCCHAR, i, name, iret) 
      i=len_trim(unit)
      call ncaptc(ncid, depth_tid, 'units',     NCCHAR, i, unit, iret) 
c
      name = 'Depth of W Grid points  '
      unit = 'm '
      i=len_trim(name)
      call ncaptc(ncid, depth_wid, 'long_name', NCCHAR, i, name, iret) 
      i=len_trim(unit)
      call ncaptc(ncid, depth_wid, 'units',     NCCHAR, i, unit, iret) 

      call ncendf(ncid, iret)

      allocate( var(imt,1) )
      Corner=0; edges=0
      Corner(1) = 1; edges(1) = imt
      var(1:imt,1)=xt
      call ncvpt(ncid, Lon_tid, corner, edges,var, iret)
      Corner(1) = 1; edges(1) = imt
      var(1:imt,1)=xu
      call ncvpt(ncid, Lon_uid, corner, edges,var, iret)
      deallocate(var);allocate(var(jmt,1))
      corner(1) = 1; edges(1) = jmt
      var(1:jmt,1)=yt
      call ncvpt(ncid, Lat_tid, corner, edges,var, iret)
      corner(1) = 1; edges(1) = jmt
      var(1:jmt,1)=yu
      call ncvpt(ncid, Lat_uid, corner, edges,var, iret)
      deallocate(var);allocate(var(km,1))
      corner(1) = 1; edges(1) = km
      var(1:km,1)=zt/100.
      call ncvpt(ncid, depth_tid, corner, edges,var, iret)
      deallocate(var);allocate(var(km+1,1))
      var(2:km+1,1)=zw(1:km)/100.;var(1,1)=0.
      corner(1) = 1; edges(1) = km+1
      call ncvpt(ncid, depth_wid, corner, edges,var, iret)

      deallocate(var)
c
      call ncclos (ncid, iret)
      end subroutine init_forcing_file





      subroutine init_out_cdf()
      use prep_module
      implicit none
#include "netcdf.inc"
      integer   ncid,iret,n
      integer lon_tdim,lat_tdim,itimedim,depth_tdim,depth_wid
      integer lon_tid ,lat_tid ,itimeid, depth_tid,depth_wdim
      integer lon_udim,lat_udim
      integer lon_uid ,lat_uid
      integer ixtid,iytid,ixuid,iyuid,iztid,izwid
      integer dxtid,dytid,dxuid,dyuid,dztid,dzwid
      integer ikmtid, iroseid,iraw_roseid,iglatid
      integer itauxid,itauyid, ihtpid,topoid,rosebud,hupid
      integer dims(4), corner(4), edges(4),i
#ifdef netcdf_real4
      real (kind=4), allocatable :: var(:,:)
#else
      real , allocatable :: var(:,:)
#endif
      real (kind=8), allocatable :: v8(:)

      character name*80, unit*80, text*80

      ncid = nccre ('forcing.cdf', NCCLOB, iret)
      iret=nf_set_fill(ncid, NF_NOFILL, iret)
      i=imt; lon_tdim  = ncddef(ncid, 'Longitude_t', i, iret)
      i=jmt; Lat_tdim  = ncddef(ncid, 'Latitude_t',  i, iret)
      i=imt; lon_udim  = ncddef(ncid, 'Longitude_u', i, iret)
      i=jmt; Lat_udim  = ncddef(ncid, 'Latitude_u',  i, iret)
      i=km;  depth_tdim= ncddef(ncid, 'depth_t',  i, iret)
      i=km+1;depth_wdim= ncddef(ncid, 'depth_w',  i, iret)

      dims(1)  = Lon_tdim
      i=1
      Lon_tid  = ncvdef (ncid,'Longitude_t',NCFLOAT,i,dims,iret)
      dims(1)  = Lat_tdim
      Lat_tid  = ncvdef (ncid,'Latitude_t', NCFLOAT,i,dims,iret)
      dims(1)  = Lon_udim
      Lon_uid  = ncvdef (ncid,'Longitude_u',NCFLOAT,i,dims,iret)
      dims(1)  = Lat_udim
      Lat_uid  = ncvdef (ncid,'Latitude_u', NCFLOAT,i,dims,iret)
      dims(1)  = depth_tdim
      depth_tid = ncvdef (ncid,'depth_t', NCFLOAT,i,dims,iret)
      dims(1)  = depth_wdim
      depth_wid = ncvdef (ncid,'depth_w', NCFLOAT,i,dims,iret)
c
      name = 'Longitude on rotated T grid'
      unit = 'degrees_W       '
      i=len_trim(name)
      call ncaptc(ncid, Lon_tid, 'long_name', NCCHAR, i, name, iret) 
      i=len_trim(unit)
      call ncaptc(ncid, Lon_tid, 'units',     NCCHAR, i, unit, iret) 

      name = 'Longitude on rotated U grid'
      unit = 'degrees_W       '
      i=len_trim(name)
      call ncaptc(ncid, Lon_uid, 'long_name', NCCHAR, i, name, iret) 
      i=len_trim(unit)
      call ncaptc(ncid, Lon_uid, 'units',     NCCHAR, i, unit, iret) 

      name = 'Latitude on rotated T grid'
      unit = 'degrees_N       '
      i=len_trim(name)
      call ncaptc(ncid, Lat_tid, 'long_name', NCCHAR, i, name, iret) 
      i=len_trim(unit)
      call ncaptc(ncid, Lat_tid, 'units',     NCCHAR, i, unit, iret) 

      name = 'Latitude on rotated U grid'
      unit = 'degrees_N       '
      i=len_trim(name)
      call ncaptc(ncid, Lat_uid, 'long_name', NCCHAR, i, name, iret) 
      i=len_trim(unit)
      call ncaptc(ncid, Lat_uid, 'units',     NCCHAR, i, unit, iret) 

      name = 'Depth of T Grid points  '
      unit = 'm '
      i=len_trim(name)
      call ncaptc(ncid, depth_tid, 'long_name', NCCHAR, i, name, iret) 
      i=len_trim(unit)
      call ncaptc(ncid, depth_tid, 'units',     NCCHAR, i, unit, iret) 
c
      name = 'Depth of W Grid points  '
      unit = 'm '
      i=len_trim(name)
      call ncaptc(ncid, depth_wid, 'long_name', NCCHAR, i, name, iret) 
      i=len_trim(unit)
      call ncaptc(ncid, depth_wid, 'units',     NCCHAR, i, unit, iret) 


c     write also grid related quantities
c
      i=1
      ixtid   = ncvdef(ncid,'xt',NCdouble,i,lon_tdim,iret)
      ixuid   = ncvdef(ncid,'xu',NCdouble,i,lon_udim,iret)
      iytid   = ncvdef(ncid,'yt',NCdouble,i,lat_tdim,iret)
      iyuid   = ncvdef(ncid,'yu',NCdouble,i,lat_udim,iret)
      iztid   = ncvdef(ncid,'zt',NCdouble,i,Depth_tdim,iret)
      izwid   = ncvdef(ncid,'zw',NCdouble,i,Depth_wdim,iret)
      dxtid   = ncvdef(ncid,'dxtdeg',NCdouble,i,lon_tdim,iret)
      dxuid   = ncvdef(ncid,'dxudeg',NCdouble,i,lon_udim,iret)
      dytid   = ncvdef(ncid,'dytdeg',NCdouble,i,lat_tdim,iret)
      dyuid   = ncvdef(ncid,'dyudeg',NCdouble,i,lat_udim,iret)
      dztid   = ncvdef(ncid,'dzt',NCdouble,i,Depth_tdim,iret)
      dzwid   = ncvdef(ncid,'dzw',NCdouble,i,Depth_wdim,iret)
c
      i=2
      dims(1) = Lon_tdim
      dims(2) = Lat_tdim
      ikmtid    = ncvdef (ncid,'kmt',   NClong,i,dims,iret)
      iroseid   = ncvdef (ncid,'rose', NCFLOAT,i,dims,iret)
      iraw_roseid   = ncvdef (ncid,'raw_rose', NCFLOAT,i,dims,iret)
      iglatid   = ncvdef (ncid,'glat', NCFLOAT,i,dims,iret)
      topoid   = ncvdef (ncid,'topo', NCFLOAT,i,dims,iret)
#ifdef partial_cell
      ihtpid   = ncvdef (ncid,'htp', NCFLOAT,i,dims,iret)
      rosebud  = ncvdef (ncid,'rose_digged', NCFLOAT,i,dims,iret)
      hupid   = ncvdef (ncid,'hup', NCFLOAT,i,dims,iret)
#endif

      name = 'Model topography'
      unit = 'vertical index'
      call dvcdf(ncid,ikmtid,name,len_trim(name),unit,16,spval)

      name = 'Bottom of deepest wet T-grid point'
      unit = 'cm'
      call dvcdf(ncid,topoid,name,len_trim(name),unit,2,spval)

      name = 'Topography before disretisation'
      unit = 'm'
      call dvcdf(ncid,iroseid,name,len_trim(name),unit,1,spval)

      name = 'Topography before smoothing,etc'
      unit = 'm'
      call dvcdf(ncid,iraw_roseid,name,len_trim(name),unit,1,spval)

      name = 'Latitude in unrotated grid'
      unit = 'deg'
      call dvcdf(ncid,iglatid,name,len_trim(name),unit,3,spval)

#ifdef partial_cell
      name = 'Depth of partial bottom t cells'
      unit = 'cm'
      call dvcdf(ncid,ihtpid,name,len_trim(name),unit,2,spval)
      name = 'Depth of topography before discretisation and digging'
      unit = 'm'
      call dvcdf(ncid,rosebud,name,len_trim(name),unit,2,spval)
      name = 'Depth of partial bottom u-cells'
      unit = 'cm'
      call dvcdf(ncid,hupid,name,len_trim(name),unit,2,spval)
#endif

      call ncendf(ncid, iret)

      allocate( var(imt,1) )
      Corner=0; edges=0
      Corner(1) = 1; edges(1) = imt
      var(1:imt,1)=xt
      call ncvpt(ncid, Lon_tid, corner, edges,var, iret)
      Corner(1) = 1; edges(1) = imt
      var(1:imt,1)=xu
      call ncvpt(ncid, Lon_uid, corner, edges,var, iret)
      deallocate(var);allocate(var(jmt,1))
      corner(1) = 1; edges(1) = jmt
      var(1:jmt,1)=yt
      call ncvpt(ncid, Lat_tid, corner, edges,var, iret)
      corner(1) = 1; edges(1) = jmt
      var(1:jmt,1)=yu
      call ncvpt(ncid, Lat_uid, corner, edges,var, iret)
      deallocate(var);allocate(var(km,1))
      corner(1) = 1; edges(1) = km
      var(1:km,1)=zt/100.
      call ncvpt(ncid, depth_tid, corner, edges,var, iret)
      deallocate(var);allocate(var(km+1,1))
      var(2:km+1,1)=zw(1:km)/100.;var(1,1)=0.
      corner(1) = 1; edges(1) = km+1
      call ncvpt(ncid, depth_wid, corner, edges,var, iret)


      deallocate(var)
c
c
      allocate(v8(imt))
      corner(1) = 1; edges(1) = imt; v8=xt
      call ncvpt(ncid, ixtid, corner, edges,v8, iret)
      corner(1) = 1; edges(1) = imt; v8=xu
      call ncvpt(ncid, ixuid, corner, edges,v8, iret)
      corner(1) = 1; edges(1) = imt; v8=dxtdeg
      call ncvpt(ncid, dxtid, corner, edges,v8, iret)
      corner(1) = 1; edges(1) = imt; v8=dxudeg
      call ncvpt(ncid, dxuid, corner, edges,v8, iret)

      deallocate(v8);allocate(v8(jmt))
      corner(1) = 1; edges(1) = jmt; v8=yt
      call ncvpt(ncid, iytid, corner, edges,v8, iret)
      corner(1) = 1; edges(1) = jmt; v8=yu
      call ncvpt(ncid, iyuid, corner, edges,v8, iret)
      corner(1) = 1; edges(1) = jmt; v8=dytdeg
      call ncvpt(ncid, dytid, corner, edges,v8, iret)
      corner(1) = 1; edges(1) = jmt; v8=dyudeg
      call ncvpt(ncid, dyuid, corner, edges,v8, iret)

      deallocate(v8);allocate(v8(km))
      corner(1) = 1; edges(1) = km; v8=zt
      call ncvpt(ncid, iztid, corner, edges,v8, iret)
      corner(1) = 1; edges(1) = km; v8=dzt
      call ncvpt(ncid, dztid, corner, edges,v8, iret)
      deallocate(v8);allocate(v8(km+1))
      corner(1) = 1; edges(1) = km+1; v8(2:km+1)=zw; v8(1)=0.
      call ncvpt(ncid, izwid, corner, edges,v8, iret)
      corner(1) = 1; edges(1) = km+1; v8(1:km+1)=dzw(0:km)
      call ncvpt(ncid, dzwid, corner, edges,v8, iret)

      deallocate(v8)
c
      call ncclos (ncid, iret)

      end


      subroutine write_kmt_to_cdf(nrose,srose
#ifdef partial_cell
     & ,rose_digged
#endif
     &     )
      use prep_module
      implicit none
#include "netcdf.inc"
      real nrose(imt,jmt),srose(imt,jmt)
#ifdef partial_cell
      real rose_digged(imt,jmt)
#endif
      integer  ncid,iret,htpid,topoid,i,j
      integer  ikmtid, iroseid,iraw_roseid
      integer  corner(4), edges(4)
      integer , allocatable :: ivar(:,:)
#ifdef netcdf_real4
      real (kind=4), allocatable :: var(:,:)
#else
      real , allocatable :: var(:,:)
#endif

      iret=nf_open('forcing.cdf',NF_WRITE,ncid)
      iret=nf_set_fill(ncid, NF_NOFILL, iret)
      iret=nf_inq_varid(ncid,'kmt',ikmtid)
      iret=nf_inq_varid(ncid,'topo',topoid)
      iret=nf_inq_varid(ncid,'rose',iroseid)
      iret=nf_inq_varid(ncid,'raw_rose',iraw_roseid)
      corner=1; edges(1)=imt; edges(2)=jmt
      allocate(ivar(imt,jmt) )
      ivar=kmt
      call ncvpt(ncid, ikmtid, corner, edges,ivar, iret)
      deallocate(ivar);allocate(var(imt,jmt) )
      var=srose/100.
      call ncvpt(ncid, iroseid, corner, edges,var, iret)
      var=nrose/100.
      call ncvpt(ncid, iraw_roseid, corner, edges,var, iret)
#ifdef partial_cell
      iret=nf_inq_varid(ncid,'htp',htpid)
      var=htp
      call ncvpt(ncid, htpid, corner, edges,var, iret)
      iret=nf_inq_varid(ncid,'rose_digged',htpid)
      var=rose_digged
      call ncvpt(ncid, htpid, corner, edges,var, iret)
      iret=nf_inq_varid(ncid,'hup',htpid)
      var=hup
      call ncvpt(ncid, htpid, corner, edges,var, iret)
#endif
      do i=1,imt
       do j=1,jmt
        if (kmt(i,j) == 0 ) then 
         var(i,j)=0.0
        else
         var(i,j)=zw(kmt(i,j))
        endif
       enddo
      enddo
      call ncvpt(ncid, topoid, corner, edges,var, iret)
      call ncclos (ncid, iret)
      deallocate(var)
      end


      subroutine write_glat_to_cdf(glat)
      use prep_module
      implicit none
#include "netcdf.inc"
      real glat(imt,jmt)
      integer  ncid,iret
      integer iglatid
      integer corner(4), edges(4)
#ifdef netcdf_real4
      real (kind=4) , allocatable :: var(:,:)
#else
      real , allocatable :: var(:,:)
#endif
      iret=nf_open('forcing.cdf',NF_WRITE,ncid)
      iret=nf_set_fill(ncid, NF_NOFILL, iret)
      iret=nf_inq_varid(ncid,'glat',iglatid)
      corner=1; edges(1)=imt; edges(2)=jmt
      allocate(var(imt,jmt) ); var=glat
      call ncvpt(ncid, iglatid, corner, edges,var, iret)
      call ncclos (ncid, iret)
      deallocate(var)
      end


      subroutine write_tau_to_cdf(tau,nt,k,time,toffset)
!     write the k.th time level
!     if k=1 then define variables
      use prep_module
      implicit none
#include "netcdf.inc"
      integer :: nt,n,k
      real :: tau(imt,jmt,2),time(nt),toffset

      integer  ncid,iret,i
      integer  itauxid,itauyid,timedim,timeid,tstrtid,taveid
      integer lon_udim,lat_udim
      integer  corner(4), edges(4),dims(4)
#ifdef netcdf_real4
      real (kind=4), allocatable :: var(:,:)
#else
      real , allocatable :: var(:,:)
#endif
      integer :: nn
c      integer :: nn,dpm(12) = (/31,28,31,30,31,30,31,31,30,31,30,31/)
c      integer :: toffset = 0*365


      character(len=80) :: name,unit

      iret=nf_open('forcing_wind.cdf',NF_WRITE,ncid)
      iret=nf_set_fill(ncid, NF_NOFILL, iret)

      if (k==1) then

      print*,' defining variable for taux,tauy in NetCDF file'
      iret=nf_redef(ncid)

      timedim  = ncddef(ncid, 'time_tau', nt, iret)

      dims(1)=timedim
      timeid = ncvdef (ncid,'time_tau', NCFLOAT,1,dims,iret)
      tstrtid = ncvdef (ncid,'tstrt_tau', NCFLOAT,1,dims,iret)
      taveid = ncvdef (ncid,'tave_tau', NCFLOAT,1,dims,iret)

      name = 'Time'; i=len_trim(name)
      call ncaptc(ncid, timeid, 'long_name', NCCHAR, i, name, iret) 
      name = 'Start of averaging period'; i=len_trim(name)
      call ncaptc(ncid, tstrtid, 'long_name', NCCHAR, i, name, iret) 
      name = 'averaging period'; i=len_trim(name)
      call ncaptc(ncid, taveid, 'long_name', NCCHAR, i, name, iret) 
      name = 'days'; i=len_trim(name)
      call ncaptc(ncid, timeid, 'units', NCCHAR, i, name, iret) 
      call ncaptc(ncid, tstrtid, 'units', NCCHAR, i, name, iret) 
      call ncaptc(ncid, taveid, 'units', NCCHAR, i, name, iret) 
      name = '01-JAN-0001 00:00:00'; i=len_trim(name)
      call ncaptc(ncid, timeid, 'time_origin', NCCHAR, i, name, iret) 
      call ncaptc(ncid, tstrtid, 'time_origin', NCCHAR, i, name, iret) 

      name = 'noleap'; i=len_trim(name)
      call ncaptc(ncid, timeid, 'calendar',NCCHAR,i,name,iret) 
      name = 'tstrt_tau'; i=len_trim(name)
      call ncaptc(ncid, timeid, 'start_ave_period',NCCHAR,i,name,iret) 
      name = 'tave_tau'; i=len_trim(name)
      call ncaptc(ncid, timeid, 'ave_period',NCCHAR,i,name,iret) 
c
      iret=nf_inq_dimid(ncid,'Longitude_u',lon_udim)
      iret=nf_inq_dimid(ncid,'Latitude_u',lat_udim)
      dims(1) = Lon_udim; dims(2) = Lat_udim; dims(3) = timedim
      itauxid   = ncvdef (ncid,'taux', NCFLOAT,3,dims,iret)
      itauyid   = ncvdef (ncid,'tauy', NCFLOAT,3,dims,iret)
      name = 'Zonal Wind Stress     '
      unit = 'dyn/cm^2'
      call dvcdf(ncid,itauxid,name,
     &     len_trim(name),unit,len_trim(unit),spval)
      name = 'Meridional Wind Stress     '
      unit = 'dyn/cm^2'
      call dvcdf(ncid,itauyid,name,
     &     len_trim(name),unit,len_trim(unit),spval)
      iret=nf_enddef(ncid)

      allocate(var(nt,1)); var(1,1)=time(1)+toffset
      do n=2,nt
       var(n,1)=var(n-1,1)+time(n)
      enddo
      do n=1,nt
       var(n,1)=var(n,1)-time(n)/2.
      enddo
      call ncvpt(ncid, timeid,1,nt,var, iret) ! mid points

      var(1,1)=toffset
      do n=2,nt
       var(n,1)=var(n-1,1)+time(n-1)
      enddo
      call ncvpt(ncid, tstrtid,1, nt,var, iret) ! start of avering period

      var(1:nt,1)=time(1:nt)
      call ncvpt(ncid, taveid, 1, nt,var, iret) ! averaging period
 
      deallocate(var)
      endif !k==1


      iret=nf_inq_varid(ncid,'taux',itauxid)
      iret=nf_inq_varid(ncid,'tauy',itauyid)

      edges(1)=imt; edges(2)=jmt; corner=1
      allocate(var(imt,jmt) )
      corner(3)=k; edges(3)=1
      var=tau(:,:,1)
      where( var == spval ) var = spval+100.
      where( kmu==0 ) var=spval
      call ncvpt(ncid, itauxid, corner, edges,var, iret)
      var=tau(:,:,2)
      where( var == spval ) var = spval+100.
      where( kmu==0 ) var=spval
      call ncvpt(ncid, itauyid, corner, edges,var, iret)
      call ncclos (ncid, iret)
      deallocate(var)
      end subroutine



      subroutine write_tr_to_cdf(tr_clim,tr_rest,tr_flux,
     &                                      time,toffset,nt,m,k)
!     write tracer restoring time scales, climatologies and 
!     net fluxes to netcdf file

      use prep_module
      implicit none
#include "netcdf.inc"
      integer nt,n,m ! m: number of tracer, nt : nr of all time steps
      integer :: k ! number of time step to be written
      real tr_clim(imt,jmt)
      real tr_rest(imt,jmt)
      real tr_flux(imt,jmt)
      real time(nt),toffset
      integer  ncid,iret,i
      integer  timedim,timeid,tstrtid,taveid
      integer  trrid,trfid,trcid
      integer  lon_tdim,lat_tdim
      integer  corner(4), edges(4),dims(4)
      character(len=80) :: name,unit
#ifdef netcdf_real4
      real (kind=4), allocatable :: var(:,:)
#else
      real , allocatable :: var(:,:)
#endif
      integer :: nn
      character(len=80) :: tr_name
c
c     determine name for tracer
c
      if     (m==1) then; tr_name = 'sst'
      elseif (m==2) then; tr_name = 'sss'
      else
       write(tr_name,'("str_",i2)') m
       call replace_space_zero(tr_name)
      endif
      if (m>2) then
       iret=nf_open('forcing_sflx_tracer.cdf',NF_WRITE,ncid)
      else
       iret=nf_open('forcing_sflx.cdf',NF_WRITE,ncid)
      endif
      iret=nf_set_fill(ncid, NF_NOFILL, iret)
      if (k==1) then ! define variable in netcdf file
       print*,' writing variable ',tr_name,' to NetCDF file'
       iret=nf_redef(ncid)
c
c     define time related quantities
       timedim  = ncddef(ncid,'time_'//tr_name, nt, iret)
       timeid = ncvdef (ncid,'time_'//tr_name, NCFLOAT,1,timedim,iret)
       tstrtid= ncvdef (ncid,'tstrt_'//tr_name,NCFLOAT,1,timedim,iret)
       taveid = ncvdef (ncid,'tave_'//tr_name, NCFLOAT,1,timedim,iret)
       name = 'Time'; i=len_trim(name)
       call ncaptc(ncid, timeid, 'long_name', NCCHAR, i, name, iret) 
       name = 'Start of averaging period'; i=len_trim(name)
       call ncaptc(ncid, tstrtid, 'long_name', NCCHAR, i, name, iret) 
       name = 'averaging period'; i=len_trim(name)
       call ncaptc(ncid, taveid, 'long_name', NCCHAR, i, name, iret) 
       name = 'days'; i=len_trim(name)
       call ncaptc(ncid, timeid, 'units', NCCHAR, i, name, iret) 
       call ncaptc(ncid, tstrtid,'units', NCCHAR, i, name, iret) 
       call ncaptc(ncid, taveid, 'units', NCCHAR, i, name, iret) 
       name = '01-JAN-0001 00:00:00'; i=len_trim(name)
       call ncaptc(ncid, timeid, 'time_origin', NCCHAR, i, name, iret) 
       call ncaptc(ncid, tstrtid,'time_origin', NCCHAR, i, name, iret) 
       name = 'noleap'; i=len_trim(name)
       call ncaptc(ncid, timeid, 'calendar',NCCHAR,i,name,iret) 
       name = 'tstrt_'//tr_name; i=len_trim(name)
       call ncaptc(ncid, timeid, 'start_ave_period',NCCHAR,i,name,iret) 
       name = 'tave_'//tr_name; i=len_trim(name)
       call ncaptc(ncid, timeid, 'ave_period',NCCHAR,i,name,iret) 
c
c    define tracer flux related variables
       iret=nf_inq_dimid(ncid,'Longitude_t',lon_tdim)
       iret=nf_inq_dimid(ncid,'Latitude_t',lat_tdim)
       dims(1) = Lon_tdim; dims(2) = Lat_tdim; dims(3) = timedim
       name=tr_name(1:len_trim(tr_name))//'_clim'
       trcid   = ncvdef (ncid,name, NCFLOAT,3,dims,iret)
       name=tr_name(1:len_trim(tr_name))//'_rest'
       trrid   = ncvdef (ncid,name, NCFLOAT,3,dims,iret)
       name=tr_name(1:len_trim(tr_name))//'_flux'
       trfid   = ncvdef (ncid,name, NCFLOAT,3,dims,iret)

       if (m==1) then
        name = 'Appearant atmospheric temperature'; unit = 'deg C'
        call dvcdf(ncid,trcid,name,
     &     len_trim(name),unit,len_trim(unit),spval)
        name = 'Haney coefficientf or SST'; unit = 'cal/s/cm**2/K'
        call dvcdf(ncid,trrid,name,
     &     len_trim(name),unit,len_trim(unit),spval)
        name = 'Net heat flux'; unit = 'cal/s/cm**2'
        call dvcdf(ncid,trfid,name,
     &     len_trim(name),unit,len_trim(unit),spval)
       elseif (m==2) then
        name = 'Climatological salinity'; unit = 'psu'
        call dvcdf(ncid,trcid,name,
     &     len_trim(name),unit,len_trim(unit),spval)
        name = 'inverse restoring time scale '//
     &         'for SSS times 1. level thickness'; unit = 'cm/s'
        call dvcdf(ncid,trrid,name,
     &     len_trim(name),unit,len_trim(unit),spval)
        name = 'Net salt flux'; unit = 'psu/s/cm**2'
        call dvcdf(ncid,trfid,name,
     &     len_trim(name),unit,len_trim(unit),spval)
       else
        name = 'Climatological distribution'; unit = 'unknown'
        call dvcdf(ncid,trcid,name,
     &     len_trim(name),unit,len_trim(unit),spval)
        name = 'Restoring time scale ' ; unit = 'cm/s'
        call dvcdf(ncid,trrid,name,
     &     len_trim(name),unit,len_trim(unit),spval)
        name = 'Net tracer flux'; unit = 'unkown/s/cm**2'
        call dvcdf(ncid,trfid,name,
     &     len_trim(name),unit,len_trim(unit),spval)
       endif ! m==1
       iret=nf_enddef(ncid)
c
c    store time related stuff
       allocate(var(nt,1)); var(1,1)=time(1)+toffset
       do n=2,nt; var(n,1)=var(n-1,1)+time(n); enddo
       do n=1,nt; var(n,1)=var(n,1)-time(n)/2.; enddo
       call ncvpt(ncid, timeid, 1, nt,var, iret)
       var(1,1)=toffset
       do n=2,nt; var(n,1)=var(n-1,1)+time(n-1); enddo
       call ncvpt(ncid, tstrtid, 1, nt,var, iret)
       var(1:nt,1)=time(1:nt)
       call ncvpt(ncid, taveid, 1, nt,var, iret)
       deallocate(var);
      endif ! k==1
c
c     store tracer flux
      name=tr_name(1:len_trim(tr_name))//'_clim'
      iret=nf_inq_varid(ncid,name,trcid)
      name=tr_name(1:len_trim(tr_name))//'_rest'
      iret=nf_inq_varid(ncid,name,trrid)
      name=tr_name(1:len_trim(tr_name))//'_flux'
      iret=nf_inq_varid(ncid,name,trfid)
      allocate(var(imt,jmt))
      var=tr_clim(:,:)
      call ncvpt(ncid,trcid, (/1,1,k/),(/imt,jmt,1/) ,var, iret)
      var=tr_rest(:,:)
      call ncvpt(ncid,trrid, (/1,1,k/),(/imt,jmt,1/) ,var, iret)
      var=tr_flux(:,:)
      call ncvpt(ncid,trfid, (/1,1,k/),(/imt,jmt,1/) ,var, iret)
      call ncclos (ncid, iret)
      deallocate(var)
      end subroutine



      subroutine write_ustar_to_cdf(ustar,nt,k,time,toffset)
!     write the k.th time level
!     if k=1 then define variables
      use prep_module
      implicit none
#include "netcdf.inc"
      integer :: nt,n,k, ncid,iret,i
      real :: ustar(imt,jmt),time(nt),toffset
      integer  ustarid,timedim,timeid,tstrtid,taveid
      integer lon_tdim,lat_tdim, corner(4), edges(4),dims(4)
#ifdef netcdf_real4
      real (kind=4), allocatable :: var(:,:)
#else
      real , allocatable :: var(:,:)
#endif
      integer :: nn

      character(len=80) :: name,unit

      iret=nf_open('forcing_wind.cdf',NF_WRITE,ncid)
      iret=nf_set_fill(ncid, NF_NOFILL, iret)

      if (k==1) then

      print*,' defining variable for ustar in NetCDF file'
      iret=nf_redef(ncid)

      timedim  = ncddef(ncid, 'time_ustar', nt, iret)

      dims(1)=timedim
      timeid = ncvdef (ncid,'time_ustar', NCFLOAT,1,dims,iret)
      tstrtid = ncvdef (ncid,'tstrt_ustar', NCFLOAT,1,dims,iret)
      taveid = ncvdef (ncid,'tave_ustar', NCFLOAT,1,dims,iret)

      name = 'Time'; i=len_trim(name)
      call ncaptc(ncid, timeid, 'long_name', NCCHAR, i, name, iret) 
      name = 'Start of averaging period'; i=len_trim(name)
      call ncaptc(ncid, tstrtid, 'long_name', NCCHAR, i, name, iret) 
      name = 'averaging period'; i=len_trim(name)
      call ncaptc(ncid, taveid, 'long_name', NCCHAR, i, name, iret) 
      name = 'days'; i=len_trim(name)
      call ncaptc(ncid, timeid, 'units', NCCHAR, i, name, iret) 
      call ncaptc(ncid, tstrtid, 'units', NCCHAR, i, name, iret) 
      call ncaptc(ncid, taveid, 'units', NCCHAR, i, name, iret) 
      name = '01-JAN-0001 00:00:00'; i=len_trim(name)
      call ncaptc(ncid, timeid, 'time_origin', NCCHAR, i, name, iret) 
      call ncaptc(ncid, tstrtid, 'time_origin', NCCHAR, i, name, iret) 

      name = 'noleap'; i=len_trim(name)
      call ncaptc(ncid, timeid, 'calendar',NCCHAR,i,name,iret) 
      name = 'tstrt_ustar'; i=len_trim(name)
      call ncaptc(ncid, timeid, 'start_ave_period',NCCHAR,i,name,iret) 
      name = 'tave_ustar'; i=len_trim(name)
      call ncaptc(ncid, timeid, 'ave_period',NCCHAR,i,name,iret) 
c
      iret=nf_inq_dimid(ncid,'Longitude_t',lon_tdim)
      iret=nf_inq_dimid(ncid,'Latitude_t',lat_tdim)
      dims(1) = Lon_tdim; dims(2) = Lat_tdim; dims(3) = timedim

      ustarid   = ncvdef (ncid,'ustar', NCFLOAT,3,dims,iret)
      name = 'Surface friction velocity (**3)'
      unit = '(cm/s)^3'
      call dvcdf(ncid,ustarid,name,len_trim(name),
     &                 unit,len_trim(unit),spval)

      iret=nf_enddef(ncid)

      allocate(var(nt,1)); var(1,1)=time(1)+toffset
      do n=2,nt
       var(n,1)=var(n-1,1)+time(n)
      enddo
      do n=1,nt
       var(n,1)=var(n,1)-time(n)/2.
      enddo
      call ncvpt(ncid, timeid,1,nt,var, iret) ! mid points

      var(1,1)=toffset
      do n=2,nt
       var(n,1)=var(n-1,1)+time(n-1)
      enddo
      call ncvpt(ncid, tstrtid,1, nt,var, iret) ! start of avering period

      var(1:nt,1)=time(1:nt)
      call ncvpt(ncid, taveid, 1, nt,var, iret) ! averaging period
 
      deallocate(var)
      endif !k==1

      iret=nf_inq_varid(ncid,'ustar',ustarid)

      edges(1)=imt; edges(2)=jmt; corner=1
      allocate(var(imt,jmt) )
      corner(3)=k; edges(3)=1
      var=ustar
      call ncvpt(ncid,ustarid, corner, edges,var, iret)
      call ncclos (ncid, iret)
      deallocate(var)
      end subroutine




      subroutine write_blue_to_cdf(temp,salt,nt,month)
      use prep_module
      implicit none
#include "netcdf.inc"
      integer nt,n,month
      real temp(imt,jmt,km),salt(imt,jmt,km)

      integer ncid,iret,i
      integer tid,sid,timedim,timeid,tstrtid,taveid
      integer lon_tdim,lat_tdim,depth_tdim
      integer corner(4), edges(4),dims(4)
#ifdef netcdf_real4
      real (kind=4), allocatable :: var(:,:)
#else
      real , allocatable :: var(:,:)
#endif
      integer :: nn,dpm(12) = (/31,28,31,30,31,30,31,31,30,31,30,31/)
      integer :: toffset = 0*365
      character(len=80) :: name,unit

      iret=nf_open('forcing_ic.cdf',NF_WRITE,ncid)
      iret=nf_set_fill(ncid, NF_NOFILL, iret)

      if (month==1) then

       print*,' defining variables for blue T/S'

      iret=nf_redef(ncid)

      i=nt; timedim  = ncddef(ncid, 'time_blue', i, iret)

      dims(1)=timedim;i=1
      timeid = ncvdef (ncid,'time_blue', NCFLOAT,i,dims,iret)
      tstrtid = ncvdef (ncid,'tstrt_blue', NCFLOAT,i,dims,iret)
      taveid = ncvdef (ncid,'tave_blue', NCFLOAT,i,dims,iret)

      name = 'Time'; i=len_trim(name)
      call ncaptc(ncid, timeid, 'long_name', NCCHAR, i, name, iret) 

      name = 'Start of averaging period'; i=len_trim(name)
      call ncaptc(ncid, tstrtid, 'long_name', NCCHAR, i, name, iret) 

      name = 'averaging period'; i=len_trim(name)
      call ncaptc(ncid, taveid, 'long_name', NCCHAR, i, name, iret) 

      name = 'days'; i=len_trim(name)
      call ncaptc(ncid, timeid, 'units', NCCHAR, i, name, iret) 
      call ncaptc(ncid, tstrtid, 'units', NCCHAR, i, name, iret) 
      call ncaptc(ncid, taveid, 'units', NCCHAR, i, name, iret) 
      name = '01-JAN-0001 00:00:00'; i=len_trim(name)
      call ncaptc(ncid, timeid, 'time_origin', NCCHAR, i, name, iret) 
      call ncaptc(ncid, tstrtid, 'time_origin', NCCHAR, i, name, iret) 

      name = 'noleap'; i=len_trim(name)
      call ncaptc(ncid, timeid, 'calendar',NCCHAR,i,name,iret) 
      name = 'tstrt_blue'; i=len_trim(name)
      call ncaptc(ncid, timeid, 'start_ave_period',NCCHAR,i,name,iret) 
      name = 'tave_blue'; i=len_trim(name)
      call ncaptc(ncid, timeid, 'ave_period',NCCHAR,i,name,iret) 
c
      iret=nf_inq_dimid(ncid,'Longitude_t',lon_tdim)
      iret=nf_inq_dimid(ncid,'Latitude_t',lat_tdim)
      iret=nf_inq_dimid(ncid,'depth_t',depth_tdim)
      dims = (/Lon_tdim, Lat_tdim, depth_tdim, timedim/); i=4
      tid   = ncvdef (ncid,'temp_blue', NCFLOAT,i,dims,iret)
      sid   = ncvdef (ncid,'salt_blue', NCFLOAT,i,dims,iret)

      name = 'BLUE temperature'
      unit = 'degC'
      call dvcdf(ncid,tid,name,
     &     len_trim(name),unit,len_trim(unit),spval)
      name = 'BLUE salinity'
      unit = '(psu-35)/1000'
      call dvcdf(ncid,sid,name,
     &     len_trim(name),unit,len_trim(unit),spval)

      iret=nf_enddef(ncid)

      allocate(var(nt,1)); var(1,1)=dpm(1)+toffset
      do n=2,nt
       nn=mod(n-1,12)+1; var(n,1)=var(n-1,1)+dpm(nn)
      enddo
      do n=1,nt
       nn=mod(n-1,12)+1; var(n,1)=var(n,1)-dpm(nn)/2.
      enddo
      corner=1;edges=nt
      call ncvpt(ncid, timeid, corner, edges,var, iret)

      var(1,1)=toffset
      do n=2,nt
       nn=mod(n-2,12)+1; 
       var(n,1)=var(n-1,1)+dpm(nn)
      enddo
      corner=1;edges=nt
      call ncvpt(ncid, tstrtid, corner, edges,var, iret)

      do n=1,nt
       nn=mod(n-1,12)+1; 
       var(n,1)=dpm(nn)
      enddo
      corner=1;edges=nt
      call ncvpt(ncid, taveid, corner, edges,var, iret)
      deallocate(var)

       print*,' done '
      else
       iret=nf_inq_varid(ncid,'temp_blue',tid)
       iret=nf_inq_varid(ncid,'salt_blue',sid)
      endif !month == 1

      print*,' writing T/S for BLUE and month=',month
      allocate(var(imt,jmt) )
      do n=1,km
       edges = (/imt,jmt,1,1/)
       corner = (/1,1,n,month/)
       var=temp(:,:,n)
       call ncvpt(ncid,tid, corner, edges,var, iret)
       var=salt(:,:,n)
       call ncvpt(ncid,sid, corner, edges,var, iret)
      enddo
      call ncclos (ncid, iret)
      deallocate(var)

      end



      subroutine write_ic_to_cdf(tr,key)
!     write 3-D tracer fields
!     to netcdf file which is defined in init_out_cdf
!     either a certain month as initial conditions or annual values

      use prep_module
      implicit none
#include "netcdf.inc"
      character (len=*),intent(in) :: key 
      integer n,k,m
      real :: tr(imt,jmt,km,number_tr)
      integer  ncid1,iret,i,ncid2
      integer  trid(number_tr)
      integer  lon_tdim,lat_tdim,depth_dim
      integer  corner(4), edges(4),dims(4)
#ifdef netcdf_real4
      real (kind=4) , allocatable :: var(:,:)
#else
      real , allocatable :: var(:,:)
#endif
      character(len=80) :: name,unit

      if (key == 'initial') then
       print*,' writing initial conditions to NetCDF file '
      elseif (key == 'mean') then
       print*,' writing mean tracer to NetCDF file '
      else
       print*,' ERROR: wrong key in write_ic_to_cdf'
       print*,' got ',key
       print*,' but should be either mean or initial'
       call halt_stop('in write_ic_to_cdf')
      endif

      iret=nf_open('forcing_ic.cdf',NF_WRITE,ncid1)
      iret=nf_set_fill(ncid1, NF_NOFILL, iret)
      iret=nf_redef(ncid1)
c
      iret=nf_inq_dimid(ncid1,'Longitude_t',lon_tdim)
      iret=nf_inq_dimid(ncid1,'Latitude_t',lat_tdim)
      iret=nf_inq_dimid(ncid1,'depth_t',depth_dim)

      if (number_tr>2) then
       iret=nf_open('forcing_ic_tracer.cdf',NF_WRITE,ncid2)
       iret=nf_set_fill(ncid2, NF_NOFILL, iret)
       iret=nf_redef(ncid2)
c
       iret=nf_inq_dimid(ncid2,'Longitude_t',lon_tdim)
       iret=nf_inq_dimid(ncid2,'Latitude_t',lat_tdim)
       iret=nf_inq_dimid(ncid2,'depth_t',depth_dim)
      endif

      dims(1) = Lon_tdim; dims(2) = Lat_tdim; dims(3) = depth_dim

      if (key == 'initial') then
       trid(1)   = ncvdef (ncid1,'temp_ic', NCFLOAT,3,dims,iret)
       trid(2)   = ncvdef (ncid1,'salt_ic', NCFLOAT,3,dims,iret)

       name = 'Initial condition temperature'
       unit = 'deg C'
       call dvcdf(ncid1,trid(1),name,len_trim(name),
     &            unit,len_trim(unit),spval)

       name = 'Initial condition salinity'
       unit = 'psu'
       call dvcdf(ncid1,trid(2),name,len_trim(name),
     &            unit,len_trim(unit),spval)

       do m=3,number_tr

        write(name,'("tracer_",i2,"_ic")') m
        call replace_space_zero(name)
        trid(m)  = ncvdef (ncid2,name, NCFLOAT,3,dims,iret)

        write(name,'("tracer_",i2)') m
        call replace_space_zero(name)
        name = 'Initial condition for '//name
        unit = 'unknown'
        call dvcdf(ncid2,trid(m),name,len_trim(name),
     &             unit,len_trim(unit),spval)

       enddo

      elseif (key == 'mean') then

       trid(1)   = ncvdef (ncid1,'temp_mean', NCFLOAT,3,dims,iret)
       trid(2)   = ncvdef (ncid1,'salt_mean', NCFLOAT,3,dims,iret)

       name = 'Annual mean temperature'
       unit = 'deg C'
       call dvcdf(ncid1,trid(1),name,len_trim(name),
     &            unit,len_trim(unit),spval)

       name = 'Annual mean salinity'
       unit = 'psu'
       call dvcdf(ncid1,trid(2),name,len_trim(name),
     &            unit,len_trim(unit),spval)

       do m=3,number_tr

        write(name,'("tracer_",i2,"_mean")') m
        call replace_space_zero(name)
        trid(m)  = ncvdef (ncid2,name, NCFLOAT,3,dims,iret)

        write(name,'("tracer_",i2)') m
        call replace_space_zero(name)
        name = 'Annual mean '//name
        unit = 'unknown'
        call dvcdf(ncid2,trid(m),name,len_trim(name),
     &            unit,len_trim(unit),spval)

       enddo

      endif

      iret=nf_enddef(ncid1)

      edges(1)=imt; edges(2)=jmt; edges(3)=km; corner=1
      
      allocate(var(imt,jmt) )
      do k=1,km
       corner(3)=k; edges(3)=1
c       do m=1,number_tr
       do m=1,2
         var=tr(:,:,k,m)
         call ncvpt(ncid1,trid(m), corner, edges,var, iret)
       enddo
      enddo
      call ncclos (ncid1, iret)

      if (number_tr>2) then
       iret=nf_enddef(ncid2)
       do k=1,km
        corner(3)=k; edges(3)=1
        do m=3,number_tr
         var=tr(:,:,k,m)
         call ncvpt(ncid2,trid(m), corner, edges,var, iret)
        enddo
       enddo
       call ncclos (ncid2, iret)
      endif

      deallocate(var)
      end subroutine





      subroutine write_blue_mean_to_cdf(temp,salt)
!     write 3-D temperature and salinity fields
!     to netcdf file which is defined in init_out_cdf
!     here the annual values for BLUE mean

      use prep_module
      implicit none
#include "netcdf.inc"
      integer n,k
      real :: temp(imt,jmt,km),salt(imt,jmt,km)
      integer  ncid,iret,i
      integer  itid,isid
      integer  lon_tdim,lat_tdim,depth_dim
      integer  corner(4), edges(4),dims(4)
#ifdef netcdf_real4
      real (kind=4) , allocatable :: var(:,:)
#else
      real , allocatable :: var(:,:)
#endif
      character(len=80) :: name,unit

      print*,' writing mean T/S for BLUE mean to NetCDF file '

      iret=nf_open('forcing_ic.cdf',NF_WRITE,ncid)
      iret=nf_set_fill(ncid, NF_NOFILL, iret)
      iret=nf_redef(ncid)
c
      iret=nf_inq_dimid(ncid,'Longitude_t',lon_tdim)
      iret=nf_inq_dimid(ncid,'Latitude_t',lat_tdim)
      iret=nf_inq_dimid(ncid,'depth_t',depth_dim)
      dims(1) = Lon_tdim; dims(2) = Lat_tdim; dims(3) = depth_dim; i=3


      itid   = ncvdef (ncid,'temp_blue_mean', NCFLOAT,i,dims,iret)
      isid   = ncvdef (ncid,'salt_blue_mean', NCFLOAT,i,dims,iret)

      name = 'Annual mean correction temperature'
      unit = 'deg C'
      call dvcdf(ncid,itid,name,
     &     len_trim(name),unit,len_trim(unit),spval)

      name = 'Annual mean correction salinity'
      unit = 'psu'
      call dvcdf(ncid,isid,name,
     &     len_trim(name),unit,len_trim(unit),spval)

      iret=nf_enddef(ncid)

      edges(1)=imt; edges(2)=jmt; edges(3)=km; corner=1
      
      allocate(var(imt,jmt) )
      do k=1,km
       corner(3)=k; edges(3)=1
       var=temp(:,:,k)
       call ncvpt(ncid,itid, corner, edges,var, iret)
       var=salt(:,:,k)
       call ncvpt(ncid,isid, corner, edges,var, iret)
      enddo
      call ncclos (ncid, iret)
      deallocate(var)
      end subroutine


      subroutine write_ts_to_file(n,tr)
c
c     write all tracer to temporary file for month n
c
      use prep_module
      implicit none
      integer n,k,nn,m
      real :: tr(imt,jmt,km,number_tr)
      integer :: io
      print*,' writing TS for time step ',n,' to file ts.dta'

      call getunit(io,'ts.dta','usr ieee')
      do nn=1,n-1
       do k=1,km
        do m=1,number_tr
         read(io)
        enddo
       enddo
      enddo
      do k=1,km
       do m=1,number_tr
         write(io) tr(1:imt,1:jmt,k,m)
       enddo
      enddo
      close(io)
      end subroutine

      subroutine read_ts_from_file(n,tr)
c
c     read all tracer from temporary file for month n
c
      use prep_module
      implicit none
      integer n,k,nn,m
      real :: tr(imt,jmt,km,number_tr)
      integer :: io
      print*,' reading TS for time step ',n,' from file ts.dta'

      call getunit(io,'ts.dta','usr ieee')
      do nn=1,n-1
       do k=1,km
        do m=1,number_tr
         read(io)
        enddo
       enddo
      enddo
      do k=1,km
       do m=1,number_tr
         read(io) tr(1:imt,1:jmt,k,m)
       enddo
      enddo
      close(io)
      end subroutine

      subroutine read_atr_from_file(n,ntr,tr)
c
c     read a single tracer from temporary file for month n
c
      use prep_module
      implicit none
      integer n,k,nn,m,ntr
      real :: tr(imt,jmt,km)
      integer :: io
      print*,' reading TS for time step ',n,' from file ts.dta'

      call getunit(io,'ts.dta','usr ieee')
      do nn=1,n-1
       do k=1,km
        do m=1,number_tr
         read(io)
        enddo
       enddo
      enddo
      do k=1,km
       do m=1,ntr-1
         read(io) 
       enddo
       read(io) tr(1:imt,1:jmt,k)
       do m=ntr+1,number_tr
         read(io) 
       enddo
      enddo
      close(io)
      end subroutine


       subroutine write_obc_to_cdf(psi,orient,nt)
c
c      take psi profile as input for OBC and write
c      to netcdf file.
c
c      Tracer data are read from ts.dta
c
       use prep_module
       implicit none
       character (len=1) ::  orient  ! is 'n' for north, 's','w' or 'e'
       integer ::nt
c       real :: psi(nt,max(imt,jmt))
       real :: psi(nt,*)

#include "netcdf.inc"
       integer iret,ncid,timedim,dims(3),timeid,tstrtid,taveid
       integer zdim,xdim,ydim,trid(number_tr),psiid
       character (len=80) :: border,name,unit
       character (len=4) :: tr_name(number_tr)
#ifdef netcdf_real4
       real (kind=4), allocatable :: var(:,:)
#else
       real , allocatable :: var(:,:)
#endif
       integer edges(3),corner(3)
       integer :: dpm(12) = (/31,28,31,30,31,30,31,31,30,31,30,31/)
       integer :: toffset = 0*365,n,k,nn,i
       integer :: io,m
       real :: ts(imt,jmt)

       if (orient=='s') border='south'
       if (orient=='n') border='north'
       if (orient=='w') border='west'
       if (orient=='e') border='east'

       iret=nf_open('forcing_obc.cdf',NF_WRITE,ncid)
       iret=nf_set_fill(ncid, NF_NOFILL, iret)
       iret=nf_inq_dimid(ncid,'Longitude_t',xdim)
       iret=nf_inq_dimid(ncid,'Latitude_t',ydim)
       iret=nf_inq_dimid(ncid,'depth_t',zdim)

       iret=nf_redef(ncid)
       timedim  = ncddef(ncid, 'time_obc_'//border, nt, iret)
       dims(1)=timedim
       timeid = ncvdef (ncid,'time_obc_'//border, NCFLOAT,1,dims,iret)
       tstrtid= ncvdef(ncid,'tstrt_obc_'//border, NCFLOAT,1,dims,iret)
       taveid = ncvdef (ncid,'tave_obc_'//border, NCFLOAT,1,dims,iret)


       name = 'Time'; i=len_trim(name)
       call ncaptc(ncid, timeid, 'long_name', NCCHAR, i, name, iret) 
       name = 'Start of averaging period'; i=len_trim(name)
       call ncaptc(ncid, tstrtid, 'long_name', NCCHAR, i, name, iret) 
       name = 'averaging period'; i=len_trim(name)
       call ncaptc(ncid, taveid, 'long_name', NCCHAR, i, name, iret) 

       name = 'days'; i=len_trim(name)
       call ncaptc(ncid, timeid, 'units', NCCHAR, i, name, iret) 
       call ncaptc(ncid, tstrtid,'units', NCCHAR, i, name, iret) 
       call ncaptc(ncid, taveid, 'units', NCCHAR, i, name, iret) 
       name = '01-JAN-0001 00:00:00'; i=len_trim(name)
       call ncaptc(ncid, timeid, 'time_origin', NCCHAR, i, name, iret) 
       call ncaptc(ncid, tstrtid,'time_origin', NCCHAR, i, name, iret) 

       name = 'noleap'; i=len_trim(name)
       call ncaptc(ncid, timeid, 'calendar',NCCHAR,i,name,iret) 
       name = 'tstrt_obc_'//border; i=len_trim(name)
       call ncaptc(ncid, timeid, 'start_ave_period',NCCHAR,i,name,iret) 
       name = 'tave_obc_'//border; i=len_trim(name)
       call ncaptc(ncid, timeid, 'ave_period',NCCHAR,i,name,iret) 
c
       dims=(/xdim,zdim,timedim/)
       if (orient=='w'.or.orient=='e') dims(1)=ydim

       do m=1,number_tr

        if (m==1) then
         tr_name(m) = 'temp'
        elseif (m==2) then
         tr_name(m) = 'salt'
        else
         write(tr_name(m),'("tr",i2)') m
         call replace_space_zero(tr_name(m))
        endif

        trid(m)= ncvdef (ncid,tr_name(m)//'_obc_'//border,
     &                   NCFLOAT,3,dims,iret)
       enddo


       dims(2)=timedim
       psiid  = ncvdef (ncid,'psi_obc_'//border,NCFLOAT,2,dims,iret)

       name = 'Temperature for open boundary at the '//border
       unit = 'deg C'
       call dvcdf(ncid,trid(1),name,len_trim(name),unit,
     &            len_trim(unit),spval)

       name = 'Salinity for open boundary at the '//border
       unit = 'psu'
       call dvcdf(ncid,trid(2),name,len_trim(name),unit,
     &            len_trim(unit),spval)

       do m=1,number_tr
        name = tr_name(m)//' for open boundary at the '//border
        unit = 'psu'
        call dvcdf(ncid,trid(m),name,len_trim(name),unit,
     &             len_trim(unit),spval)
       enddo


       name = 'Streamfunction for open boundary at the '//border
       unit = '1^12 Sv'
       call dvcdf(ncid,psiid,name,len_trim(name),unit,
     &            len_trim(unit),spval)

       iret=nf_enddef(ncid)

c write time variables
       allocate(var(nt,1)); var(1,1)=dpm(1)+toffset
       do n=2,nt
        nn=mod(n-1,12)+1; var(n,1)=var(n-1,1)+dpm(nn)
       enddo
       do n=1,nt
        nn=mod(n-1,12)+1; var(n,1)=var(n,1)-dpm(nn)/2.
       enddo
       corner=1;edges=nt
       call ncvpt(ncid, timeid, corner, edges,var, iret)

       var(1,1)=toffset
       do n=2,nt
        nn=mod(n-2,12)+1; 
        var(n,1)=var(n-1,1)+dpm(nn)
       enddo
       corner=1;edges=nt
       call ncvpt(ncid, tstrtid, corner, edges,var, iret)

       do n=1,nt
        nn=mod(n-1,12)+1; 
        var(n,1)=dpm(nn)
       enddo
       corner=1;edges=nt
       call ncvpt(ncid, taveid, corner, edges,var, iret)
       deallocate(var)

c write psi

       corner=1; edges(2)=1
       if (orient=='w'.or.orient=='e') then
         allocate(var(jmt,1)) 
         edges(1)=jmt
       else
         allocate(var(imt,1)) 
         edges(1)=imt
       endif
       do n=1,nt
        corner(2)=n
        if (orient=='n') var(:,1)=psi(n,1:imt)
        if (orient=='s') var(:,1)=psi(n,1:imt)
        if (orient=='e') var(:,1)=psi(n,1:jmt)
        if (orient=='w') var(:,1)=psi(n,1:jmt)
        call ncvpt(ncid, psiid, corner, edges,var, iret)
       enddo
       deallocate(var)

c write temp/sal
      
       corner=1; edges(2:3)=(/1,1/)
       if (orient=='w'.or.orient=='e') then
         allocate(var(jmt,1)) 
         edges(1)=jmt
       else
         allocate(var(imt,1)) 
         edges(1)=imt
       endif

       call getunit(io,'ts.dta','usr ieee')
       do n=1,nt
        corner(3)=n
        do k=1,km
         corner(2)=k
         do m=1,number_tr
          read(io) ts(:,:)
          if (orient=='n') var(:,1)=ts(:,jmt-1)
          if (orient=='s') var(:,1)=ts(:,2)
          if (orient=='e') var(:,1)=ts(imt-1,:)
          if (orient=='w') var(:,1)=ts(2,:)
          call ncvpt(ncid, trid(m), corner, edges,var, iret)
         enddo ! m
        enddo ! k
       enddo ! n
       close(io)
       deallocate(var)

       call ncclos (ncid, iret)

       end subroutine write_obc_to_cdf



      subroutine write_spg_to_cdf(tscl,tr,is,ie,js,je,nt,nl,ntr)
      use prep_module
      implicit none
      integer , intent(in) :: nl,ntr
      integer is,ie,js,je,nt
      real tscl(is:ie,js:je,km)
      real tr(is:ie,js:je,km)
#include "netcdf.inc"
      integer ncid,iret,nsp,i
      character(len=2):: number
      character(len=4):: tr_name
      character(len=80):: name,unit
      integer londim,latdim,depthdim,lonid,latid,trid,tsclid
      integer timedim,timeid,tstrtid,taveid
      integer dims(4),edges(4),corner(4)
#ifdef netcdf_real4
      real (kind=4), allocatable :: var(:,:)
#else
      real , allocatable :: var(:,:)
#endif
      integer :: dpm(12) = (/31,28,31,30,31,30,31,31,30,31,30,31/)
      integer :: toffset = 0*365,nn,n,k
      
c see how many sponge layers are already defined
      iret=nf_open('forcing_sponge.cdf',NF_noWRITE,ncid)
      iret= nf_get_att_int(ncid,nf_global,'number_of_sponge_layers',nsp)
      if (iret /= 0) nsp = 0
      call ncclos (ncid, iret)

      if (nl==1 .and. ntr==1) then
       nsp = nsp+1
       print*,' defining sponge layer # ',nsp
       print*,' for tracer #',ntr
      elseif (nl==1 .and. ntr>1) then
       print*,' defining tracer #',ntr,' in sponge layer #',nsp
      else
       print*,' writing data for sponge layer # ',nsp
       print*,' for tracer #',ntr
       print*,' for month = ',nl
      endif ! nl==1

      write(number,'(i2)') nsp
      if (number(1:1) == ' ') number(1:1)='0'

      if (ntr==1) then
         tr_name = 'temp'
      elseif (ntr==2) then
         tr_name = 'salt'
      else
       write(tr_name,'("tr",i2)') ntr
       call replace_space_zero(tr_name)
      endif

      iret=nf_open('forcing_sponge.cdf',NF_WRITE,ncid)
      iret=nf_set_fill(ncid, NF_NOFILL, iret)
  
      if (nl==1) then

      iret=nf_redef(ncid)
c redefine number of sponge layers
      iret= nf_put_att_int(ncid,nf_global,'number_of_sponge_layers',
     &       nf_int,1,nsp)

c define time variables for the sponge layer
      name = 'time_spg_'//number//'_'//tr_name
      timedim  = ncddef(ncid, name, nt, iret)

      dims(1)=timedim
      name = 'time_spg_'//number//'_'//tr_name
      timeid = ncvdef (ncid,name, NCFLOAT,1,dims,iret)
      name = 'tstrt_spg_'//number//'_'//tr_name
      tstrtid= ncvdef(ncid,name, NCFLOAT,1,dims,iret)
      name = 'tave_spg_'//number//'_'//tr_name
      taveid = ncvdef (ncid,name, NCFLOAT,1,dims,iret)

      name = 'Time'; i=len_trim(name)
      call ncaptc(ncid, timeid, 'long_name', NCCHAR, i, name, iret) 
      name = 'Start of averaging period'; i=len_trim(name)
      call ncaptc(ncid, tstrtid, 'long_name', NCCHAR, i, name, iret) 
      name = 'averaging period'; i=len_trim(name)
      call ncaptc(ncid, taveid, 'long_name', NCCHAR, i, name, iret) 

      name = 'days'; i=len_trim(name)
      call ncaptc(ncid, timeid, 'units', NCCHAR, i, name, iret) 
      call ncaptc(ncid, tstrtid,'units', NCCHAR, i, name, iret) 
      call ncaptc(ncid, taveid, 'units', NCCHAR, i, name, iret) 
      name = '01-JAN-0001 00:00:00'; i=len_trim(name)
      call ncaptc(ncid, timeid, 'time_origin', NCCHAR, i, name, iret) 
      call ncaptc(ncid, tstrtid,'time_origin', NCCHAR, i, name, iret) 

      name = 'noleap'; i=len_trim(name)
      call ncaptc(ncid, timeid, 'calendar',NCCHAR,i,name,iret) 
      name = 'tstrt_spg_'//number//'_'//tr_name; i=len_trim(name)
      call ncaptc(ncid, timeid, 'start_ave_period',NCCHAR,i,name,iret) 
      name = 'tave_spg_'//number//'_'//tr_name; i=len_trim(name)
      call ncaptc(ncid, timeid, 'ave_period',NCCHAR,i,name,iret) 

      if (ntr==1) then
c
c define grid variables for the sponge layer
c
       i=ie-is+1
       londim  = ncddef(ncid, 'Lon_spg_'//number,i, iret)
       i=je-js+1
       latdim  = ncddef(ncid, 'Lat_spg_'//number,i, iret)
       iret=nf_inq_dimid(ncid,'depth_t',depthdim)

       dims(1)=londim
       Lonid  = ncvdef (ncid,'Lon_spg_'//number,NCFLOAT,1,dims,iret)
       dims(1)=latdim
       Latid  = ncvdef (ncid,'Lat_spg_'//number,NCFLOAT,1,dims,iret)

       name = 'Longitude of sponge layer #'//number
       unit = 'degrees_W       '
       i=len_trim(name)
       call ncaptc(ncid, Lonid, 'long_name', NCCHAR, i, name, iret) 
       i=len_trim(unit)
       call ncaptc(ncid, Lonid, 'units',     NCCHAR, i, unit, iret) 

       name = 'Latitude of sponge layer #'//number
       unit = 'degrees_N       '
       i=len_trim(name)
       call ncaptc(ncid, Latid, 'long_name', NCCHAR, i, name, iret) 
       i=len_trim(unit)
       call ncaptc(ncid, Latid, 'units',     NCCHAR, i, unit, iret) 
      else
       iret=nf_inq_dimid(ncid,'Lon_spg_'//number,londim)
       iret=nf_inq_dimid(ncid,'Lat_spg_'//number,latdim)
       iret=nf_inq_dimid(ncid,'depth_t',depthdim)
      endif
c
c define variable for tracers and time scale
c
      dims(1)=londim; dims(2)=latdim; dims(3)=depthdim; dims(4)=timedim
      trid = ncvdef (ncid,tr_name//'_spg_'//number,NCFLOAT,4,dims,iret)


      dims(1)=londim; dims(2)=latdim; dims(3)=depthdim; dims(4)=timedim
      tsclid  = ncvdef (ncid,
     &          tr_name//'_tscl_spg_'//number,NCFLOAT,4,dims,iret)

      if (ntr==1) then
      name = 'Temperature for sponge layer #'//number
      unit = 'deg C'
      elseif (ntr==2) then
      name = 'Salinity for sponge layer #'//number
      unit = 'psu'
      else
       name = tr_name//' for sponge layer #'//number
       unit = 'unknown'
      endif
      call dvcdf(ncid,trid,name,len_trim(name),unit,
     &           len_trim(unit),spval)
      dims=(/is,ie,js,je/)
      iret= nf_put_att_int(ncid,trid,'range',nf_int,4,dims)

      name = 'Inverse time scale for '//tr_name
     &         //' of sponge layer #'//number
      unit = '1/s'
      call dvcdf(ncid,tsclid,name,len_trim(name),unit,
     &           len_trim(unit),spval)
      iret= nf_put_att_int(ncid,tsclid,'range',nf_int,4,dims)

      iret=nf_enddef(ncid)

      if (ntr==1) then
c write grid variables

       allocate(var(is:ie,1)); var(is:ie,1)=xt(is:ie)
       corner=1;edges=ie-is+1
       call ncvpt(ncid, lonid, corner, edges,var(is,1), iret)
       deallocate(var)

       allocate(var(js:je,1)); var(js:je,1)=yt(js:je)
       corner=1;edges=je-js+1
       call ncvpt(ncid, latid, corner, edges,var(js,1), iret)
       deallocate(var)

      endif

c write time variables
      allocate(var(nt,1)); var(1,1)=dpm(1)+toffset
      do n=2,nt
       nn=mod(n-1,12)+1; var(n,1)=var(n-1,1)+dpm(nn)
      enddo
      do n=1,nt
       nn=mod(n-1,12)+1; var(n,1)=var(n,1)-dpm(nn)/2.
      enddo
      corner=1;edges=nt
      call ncvpt(ncid, timeid, corner, edges,var, iret)

      var(1,1)=toffset
      do n=2,nt
       nn=mod(n-2,12)+1; 
       var(n,1)=var(n-1,1)+dpm(nn)
      enddo
      corner=1;edges=nt
      call ncvpt(ncid, tstrtid, corner, edges,var, iret)

      do n=1,nt
       nn=mod(n-1,12)+1; 
       var(n,1)=dpm(nn)
      enddo
      corner=1;edges=nt
      call ncvpt(ncid, taveid, corner, edges,var, iret)
      deallocate(var)

      endif ! nl==1

c write temp/sal
      iret=nf_inq_varid(ncid,tr_name//'_spg_'//number,trid)
      iret=nf_inq_varid(ncid,tr_name//'_tscl_spg_'//number,tsclid)

      allocate(var(is:ie,js:je)) 
      corner=1;
      edges(1)=ie-is+1;edges(2)=je-js+1; 
      edges(3)=1; edges(4)=1
      corner(4)=nl
      do k=1,km
        corner(3)=k
        var(is:ie,js:je)=tr(is:ie,js:je,k)
        call ncvpt(ncid, trid, corner, edges,var(is,js), iret)
        var(is:ie,js:je)=tscl(is:ie,js:je,k)
        call ncvpt(ncid, tsclid, corner, edges,var(is,js), iret)
      enddo

      deallocate(var)
      call ncclos (ncid, iret)

      print*,'done'

      end subroutine



      subroutine dvcdf(ncid,ivarid,name,iname,unit,iunit,spval)
c-----------------------------------------------------------------------
c     define some standard attributes of variable ivarid in NetCDF file ncid 
c-----------------------------------------------------------------------
      implicit none
      integer ncid,ivarid,iname,iunit,iret
      character (len=*) :: name, unit
      real :: spval
#ifdef netcdf_real4
      real (kind=4) :: vv
#else
      real :: vv
#endif
#include "netcdf.inc"
      vv=spval
      call ncaptc(ncid,ivarid, 'long_name', NCCHAR,iname , name, iret) 
        if (iret.ne.0) print*,nf_strerror(iret)
      call ncaptc(ncid,ivarid, 'units',     NCCHAR,iunit, unit, iret) 
        if (iret.ne.0) print*,nf_strerror(iret)
      call ncapt (ncid,ivarid, 'missing_value',NCFLOAT,1,vv,iret)
        if (iret.ne.0) print*,nf_strerror(iret)
      call ncapt (ncid,ivarid, '_FillValue', NCFLOAT, 1,vv, iret)
        if (iret.ne.0) print*,nf_strerror(iret)
      end subroutine dvcdf


