#include "options.inc"

c
c-----------------------------------------------------------------------
c     convert the nice netcdf forcing file to binary files
c     so that we do not have to use netcdf library for IO
c
c     output is grid.dta/ kmt.dta in the old FLAME format
c          temp.mom.ic/salt.mom.ic also in that old format
c     Byte length for real depends on this cpp directive
#define real_8
c
c
c     taux.mom,tauy.mom is a different format 
c     sst_clim.mom sst_rest.mom sss_clim.mom sss_rest.mom
c     and ustar.mom also in special format
c
c     writes also blue T/S files if present in netcdf file
c
c     to do: write file for glat 
c
#define new_cdf_file_names
c     use new outfile file naming for netcdf files
c#define apply_topo_mask_to_sbc
c     apply topographic mask to sbc to zero out land values
c     (this was the default up to Apr 2005).
c-----------------------------------------------------------------------
c

      module kind_mod
#ifdef real_8
      integer, parameter::  r_knd  = selected_real_kind(13)
#else
      integer, parameter::  r_knd  = selected_real_kind(6)
#endif
      end module kind_mod


      module prep_to_bin_mod
      use kind_mod
      use time_manager_module
      implicit none
      integer :: imt,jmt,km
      integer , allocatable :: kmt(:,:),kmu(:,:)
      real (kind=r_knd), allocatable, dimension(:) :: xt,yt,zt
      character (len=80) :: 
#ifdef new_cdf_file_names
     &     grid_file    = 'forcing_folded.cdf',
     &     ic_file      = 'forcing_folded_ic.cdf',
     &     ic_tr_file   = 'forcing_folded_ic_tracer.cdf',
     &     wind_file    = 'forcing_folded_wind.cdf',
     &     sflx_file    = 'forcing_folded_sflx.cdf',
     &     sflx_tr_file = 'forcing_folded_sflx_tracer.cdf',
     &     obc_file     = 'forcing_folded_obc.cdf',
     &     sponge_file  = 'forcing_folded_sponge.cdf'
#else
     &     grid_file    = 'forcing_folded.cdf',
     &     ic_file      = 'forcing_folded.cdf',
     &     ic_tr_file   = 'forcing_folded.cdf',
     &     wind_file    = 'forcing_folded.cdf',
     &     sflx_file    = 'forcing_folded.cdf',
     &     sflx_tr_file = 'forcing_folded.cdf',
     &     obc_file     = 'forcing_folded.cdf',
     &     sponge_file  = 'forcing_folded.cdf'
#endif
      end module prep_to_bin_mod


      program prep_to_bin
      use prep_to_bin_mod
      implicit none
      character (len=32) :: stamp 
      integer :: m,ret
      character (len=80) :: vname,fname

#include "netcdf.inc"
      integer   ncid,iret, xtid,ytid,ztid,dims(4)

      stamp ='m/d/y= 1/ 1/1900, h:m:s= 0: 0: 0' 
      call init_time_manager(3600.,5.0,5.0,stamp,0)

      iret=nf_open(grid_file,NF_noWRITE,ncid)

      iret=nf_inq_varid(ncid,'xt',xtid)
      iret=nf_inq_varid(ncid,'yt',ytid)
      iret=nf_inq_varid(ncid,'zt',ztid)
      iret=NF_INQ_VARDIMID (ncid,xtid,dims)
      iret=NF_INQ_DIMLEN   (ncid, dims(1), imt)
      iret=NF_INQ_VARDIMID (ncid,ytid,dims)
      iret=NF_INQ_DIMLEN   (ncid, dims(1), jmt)
      iret=NF_INQ_VARDIMID (ncid,ztid,dims)
      iret=NF_INQ_DIMLEN   (ncid, dims(1), km)
      call ncclos (ncid, iret)

      print*,'imt=',imt,' jmt=',jmt,' km=',km

      allocate( xt(imt), yt(jmt),  zt(km))
      allocate( kmt(imt,jmt), kmu(imt,jmt))


      call write_grid()
      call write_topo()
      call write_ic()
      call write_blue()
      call write_blue_mean()
      call write_sbc('taux','taux.mom',.true.,ret)
      call write_sbc('tauy','tauy.mom',.true.,ret)
      call write_sbc('ustar','ustar.mom',.false.,ret)

      call write_sbc('sst_clim','sst_clim.mom',.false.,ret)
      call write_sbc('sst_rest','sst_rest.mom',.false.,ret)
      call write_sbc('sst_flux','sst_flux.mom',.false.,ret)

      call write_sbc('sss_clim','sss_clim.mom',.false.,ret)
      call write_sbc('sss_rest','sss_rest.mom',.false.,ret)
      call write_sbc('sss_flux','sss_flux.mom',.false.,ret)

      m=3
 20   continue
       write(vname,'("str_",i2,"_clim")') m
       call replace_space_zero(vname)
       write(fname,'("str_",i2,"_clim.mom")') m
       call replace_space_zero(fname)
       call write_sbc(vname,fname,.false.,ret)

       write(vname,'("str_",i2,"_rest")') m
       call replace_space_zero(vname)
       write(fname,'("str_",i2,"_rest.mom")') m
       call replace_space_zero(fname)
       call write_sbc(vname,fname,.false.,ret)

       write(vname,'("str_",i2,"_flux")') m
       call replace_space_zero(vname)
       write(fname,'("str_",i2,"_flux.mom")') m
       call replace_space_zero(fname)
       call write_sbc(vname,fname,.false.,ret)

      if (ret==0) then; m=m+1; goto 20 ; endif

      call write_obc('south','obc_south.mom',1000)
      call write_obc('north','obc_north.mom',1000)
      call write_obc('east','obc_east.mom',1000)
      call write_obc('west','obc_west.mom',1000)

      call write_sponge(1000)

      call write_sbc('tauxBB','taux_BB.mom',.true.,ret)
      call write_sbc('tauyBB','tauy_BB.mom',.true.,ret)
      call write_sbc('hflxBB','hflx_BB.mom',.false.,ret)
      call write_sbc('sstm_BB','sst_mean_BB.mom',.false.,ret)
      call write_sbc('sstf_BB','sst_feed_BB.mom',.false.,ret)

c      call ncclos (ncid, iret)

      end program prep_to_bin





      subroutine write_sponge(ntr_max)
      use prep_to_bin_mod
      implicit none
      integer :: ntr_max
      character (len=80)::  vname,name,fname
      character (len=32):: stamp
      integer :: itm,n,id,tid,k,nspg,i,kk
      integer :: is,ie,js,je,io
      type( time_type) :: time
      character(len=2) :: number
      character(len=4) :: tr_name
      integer :: m,ntr

      real (kind=r_knd), allocatable :: trbuf(:,:,:),aprec(:)
      real (kind=r_knd), allocatable :: tsbuf(:,:,:)
      real (kind=r_knd) :: tstart

#include "netcdf.inc"
      integer   ncid,iret, start(4),count(4),dims(4)

      iret=nf_open(sponge_file,NF_noWRITE,ncid)

      if (iret /=0) then
       print*,' cannot find Sponge netcdf file '
       return
      endif 

      iret= nf_get_att_int(ncid,nf_global,'number_of_sponge_layers',
     &                     nspg  )

      do k=1,nspg

       print*,' processing sponge layer #',k

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

       vname='temp_spg_'//number
       iret=nf_inq_varid(ncid,vname,id)
       if (iret /=0) then
        print*,' cannot find variable ',vname ,' in netcdf file'
        return
       endif

       m=3
 20    continue
        write(tr_name,'("tr",i2)') m
        call replace_space_zero(tr_name)
        iret=nf_inq_varid(ncid,tr_name//'_spg_'//number,id)
       if (iret==0) then; m=m+1; goto 20; endif
       ntr=m-1
       print*,' found ',ntr,' variables for sponge layer in netcdf file'
       if (ntr_max < ntr) then
        print*,' However, restricting number to ',ntr_max
        ntr=ntr_max
       endif

       iret=nf_inq_varid(ncid,'temp_spg_'//number,id)
       iret=NF_INQ_VARDIMID (NCID,id,dims)
       iret=NF_INQ_DIMLEN (NCID, dims(4), itm)
       name='';iret=NF_INQ_DIMname (NCID, dims(4), name)
       iret=nf_inq_varid(ncid,name,tid)

       print*,' time steps for sponge layer=',itm

       allocate( aprec(itm) )

       name='';iret= nf_get_att_text(ncid,tid ,'ave_period',name)
       iret=nf_inq_varid(ncid,name,id)
#ifdef real_8
       iret= nf_get_vara_double(ncid,id ,1,itm,aprec)
#else
       iret= nf_get_vara_real(ncid,id ,1,itm,aprec)
#endif

       name='';iret= nf_get_att_text(ncid,tid ,'start_ave_period',name)
       iret=nf_inq_varid(ncid,name,id)

#ifdef real_8
       iret= nf_get_vara_double(ncid,id,1,1,tstart);
#else
       iret= nf_get_vara_real(ncid,id,1,1,tstart);
#endif
       call set_time(time,
     &                int((tstart-int(tstart))*24*60*60),int(tstart))
       call set_stamp (stamp, time)

       iret=nf_inq_varid(ncid,vname,id)
       iret= nf_get_att_int(ncid,id,'range',dims )
       is=dims(1);ie=dims(2);js=dims(3);je=dims(4)

       print*,' range = ',is,ie,js,je
       allocate( trbuf(is:ie,js:je,ntr), tsbuf(is:ie,js:je,ntr)  )

       write(name, '("sponge_",i2,".mom")') k
       do i=1,len_trim(name)
         if (name(i:i) == ' ') name(i:i)='0'
       enddo

       print*,' writing sponge layer # ',k,' to file ',
     &   name(1:len_trim(name))

       call getunit(io,name,'usr ieee')
       write (io) nspg,itm,is,ie,js,je,stamp,aprec

       start=(/1,1,1,1/); count=(/ie-is+1,je-js+1,1,1/); 
       do n=1,itm
        start(4)=n
        do kk=1,km
         start(3)=kk
         do m=1,ntr

          if (m==1) then
           tr_name = 'temp'
          elseif (m==2) then
           tr_name = 'salt'
          else
           write(tr_name,'("tr",i2)') m
           call replace_space_zero(tr_name)
          endif
         
          iret=nf_inq_varid(ncid,tr_name//'_spg_'//number,id)
          if (iret/=0) then; print*,' cannot find var.';stop;endif
#ifdef real_8
          iret= nf_get_vara_double(ncid,id,start,count,trbuf(:,:,m))
#else
          iret= nf_get_vara_real(ncid,id,start,count,trbuf(:,:,m))
#endif
          iret=nf_inq_varid(ncid,tr_name//'_tscl_spg_'//number,id)
          if (iret/=0) then; print*,' cannot find var.';stop;endif
#ifdef real_8
          iret= nf_get_vara_double(ncid,id,start,count,tsbuf(:,:,m))
#else
          iret= nf_get_vara_real(ncid,id,start,count,tsbuf(:,:,m))
#endif
         enddo
         write(io) trbuf,tsbuf
        enddo
       enddo
       close(io)
       deallocate(aprec, trbuf,tsbuf)
      enddo
      call ncclos (ncid, iret)
      end subroutine write_sponge


      subroutine write_obc(orient,fname,ntr_max)
      use prep_to_bin_mod
      implicit none
      integer :: ntr_max
      character (len=*) :: orient,fname
      character (len=80)::  vname,name,tstr
      character (len=4)::  tr_name
      character (len=32):: stamp
      integer :: itm,n,id,tid,fdim,io,m,ntr
      type( time_type) :: time

      real (kind=r_knd), allocatable :: trbuf(:,:,:),aprec(:)
      real (kind=r_knd), allocatable :: pbuf(:)
      real (kind=r_knd) :: tstart

#include "netcdf.inc"
      integer   ncid,iret, start(4),count(4),dims(4)

      iret=nf_open(obc_file,NF_noWRITE,ncid)

      if (iret /=0) then
       print*,' cannot find OBC netcdf file '
       return
      endif 

      print*,' writing OBC ',orient,' to file ',
     &   fname(1:len_trim(fname))

      tstr=orient
      vname = 'temp_obc_'//tstr(1:len_trim(tstr))
      
      iret=nf_inq_varid(ncid,vname,id)
      if (iret /=0) then
        print*,' cannot find variable ',vname(1:len_trim(vname)) ,
     &             ' in netcdf file'
        return
      endif
      iret=NF_INQ_VARDIMID (NCID,id,dims)
      iret=NF_INQ_DIMLEN (NCID, dims(3), itm)
      name='';iret=NF_INQ_DIMname (NCID, dims(3), name)
      iret=nf_inq_varid(ncid,name,tid)

      allocate( aprec(itm) )
      name='';iret= nf_get_att_text(ncid,tid ,'ave_period',name)
      iret=nf_inq_varid(ncid,name,id)
#ifdef real_8
      iret= nf_get_vara_double(ncid,id ,1,itm,aprec)
#else
      iret= nf_get_vara_real(ncid,id ,1,itm,aprec)
#endif

      name='';iret= nf_get_att_text(ncid,tid ,'start_ave_period',name)
      iret=nf_inq_varid(ncid,name,id)

#ifdef real_8
      iret= nf_get_vara_double(ncid,id,1,1,tstart);
#else
      iret= nf_get_vara_real(ncid,id,1,1,tstart);
#endif
      call set_time(time,
     &                int((tstart-int(tstart))*24*60*60),int(tstart))
      call set_stamp (stamp, time)

      fdim=jmt
      if (orient == 'south' .or. orient == 'north' ) fdim=imt
      print*,' orientation ',orient,' fdim=',fdim

      m=3
 20   continue
        write(tr_name,'("tr",i2)') m
        call replace_space_zero(tr_name)
        iret=nf_inq_varid(ncid,
     &     tr_name//'_obc_'//tstr(1:len_trim(tstr)),id)
      if (iret==0) then; m=m+1; goto 20; endif
      ntr=m-1
      print*,' found ',ntr,' variables for tracers in netcdf file'
      if (ntr_max < ntr) then
        print*,' However, restricting number to ',ntr_max
        ntr=ntr_max
      endif

      allocate( trbuf(fdim,km,ntr), pbuf(fdim)  )

      call getunit(io,fname,'usr ieee')

      write (io) itm,stamp,aprec
      do n=1,itm

       do m=1,ntr

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

        start=(/1,1,n,1/); count=(/fdim,km,1,1/); 
        iret=nf_inq_varid(ncid,
     &        tr_name//'_obc_'//tstr(1:len_trim(tstr)),id)
        if (iret/=0) then; print*,' cannot find var.';stop;endif
#ifdef real_8
        iret= nf_get_vara_double(ncid,id,start,count,trbuf(:,:,m))
#else
        iret= nf_get_vara_real(ncid,id,start,count,trbuf(:,:,m))
#endif
       enddo

       start=(/1,n,1,1/); count=(/fdim,1,1,1/); 
       iret=nf_inq_varid(ncid,'psi_obc_'//tstr(1:len_trim(tstr)),id)
          if (iret/=0) then; print*,' cannot find var.';stop;endif
#ifdef real_8
       iret= nf_get_vara_double(ncid,id,start,count,pbuf)
#else
       iret= nf_get_vara_real(ncid,id,start,count,pbuf)
#endif
       write(io) trbuf,pbuf
      enddo
      close(io)
      deallocate(aprec, trbuf,pbuf)
      call ncclos (ncid, iret)

      end subroutine write_obc




      subroutine write_sbc(vname,fname,ugrid,ret)
      use prep_to_bin_mod
      implicit none
      character(len=*) :: vname,fname
      logical :: ugrid
      integer :: ret
      integer :: itm,n,id,tid,io
      character(len=80):: name
      character(len=32):: stamp
      type( time_type) :: time
      real (kind=r_knd), allocatable :: buf(:,:),aprec(:)
      real (kind=r_knd) :: tstart

#include "netcdf.inc"
      integer   ncid,iret, start(4),count(4),dims(4)

      if (vname=='taux') then
       iret=nf_open(wind_file(1:len_trim(wind_file)),NF_noWRITE,ncid)
      elseif (vname=='tauy') then
       iret=nf_open(wind_file,NF_noWRITE,ncid)
      elseif (vname=='ustar') then
       iret=nf_open(wind_file,NF_noWRITE,ncid)
      elseif (vname=='sst_clim') then
       iret=nf_open(sflx_file,NF_noWRITE,ncid)
      elseif (vname=='sst_rest') then
       iret=nf_open(sflx_file,NF_noWRITE,ncid)
      elseif (vname=='sst_flux') then
       iret=nf_open(sflx_file,NF_noWRITE,ncid)
      elseif (vname=='sss_clim') then
       iret=nf_open(sflx_file,NF_noWRITE,ncid)
      elseif (vname=='sss_rest') then
       iret=nf_open(sflx_file,NF_noWRITE,ncid)
      elseif (vname=='sss_flux') then
       iret=nf_open(sflx_file,NF_noWRITE,ncid)
      else
       iret=nf_open(sflx_tr_file,NF_noWRITE,ncid)
      endif

      ret=0
      if (iret /=0) then
        print*,' cannot open netcdf file for variable',vname
        ret=-1
        return
      endif

      iret=nf_inq_varid(ncid,vname,id)
      if (iret /=0) then
        print*,' cannot find variable ',vname(1:len_trim(vname)) ,
     &         ' in netcdf file'
        ret=-1
        return
      endif

       iret=NF_INQ_VARDIMID (NCID,id,dims)
       iret=NF_INQ_DIMLEN (NCID, dims(3), itm)
       name='';iret=NF_INQ_DIMname (NCID, dims(3), name)
       iret=nf_inq_varid(ncid,name,tid)

       allocate( aprec(itm) , buf(imt,jmt) )

       name='';iret= nf_get_att_text(ncid,tid ,'ave_period',name)
       iret=nf_inq_varid(ncid,name,id)
#ifdef real_8
       iret= nf_get_vara_double(ncid,id ,1,itm,aprec)
#else
       iret= nf_get_vara_real(ncid,id ,1,itm,aprec)
#endif

       name='';iret= nf_get_att_text(ncid,tid ,'start_ave_period',name)
       iret=nf_inq_varid(ncid,name,id)

#ifdef real_8
       iret= nf_get_vara_double(ncid,id,1,1,tstart)
#else
       iret= nf_get_vara_real(ncid,id,1,1,tstart)
#endif
       call set_time(time,
     &                int((tstart-int(tstart))*24*60*60),int(tstart))
       call set_stamp (stamp, time)

       iret=nf_inq_varid(ncid,vname,id)

       print*,' writing SBC ',vname(1:len_trim(vname)),
     &        ' to file ',fname(1:len_trim(fname))

       call getunit(io,fname,'usr ieee')
       write (io) itm,stamp,aprec
       do n=1,itm
        start=(/1,1,n,1/); count=(/imt,jmt,1,1/); 
#ifdef real_8
        iret= nf_get_vara_double(ncid,id,start,count,buf)
#else
        iret= nf_get_vara_real(ncid,id,start,count,buf)
#endif

#ifdef apply_topo_mask_to_sbc
        if (ugrid) then
         where( kmu == 0 ) buf = 0.
        else
         where( kmt == 0 ) buf = 0.
        endif
#endif
        write(io) buf
       enddo
       close(io)
       deallocate(aprec,buf)
      call ncclos (ncid, iret)
      end subroutine write_sbc






      subroutine write_blue()
      use prep_to_bin_mod
      implicit none
      integer :: itm,n,tempid,saltid,tid,j,id,io1,io2
      character(len=80):: name
      character(len=32):: stamp
      type( time_type) :: time
      real (kind=r_knd), allocatable :: buf(:,:,:),aprec(:)
      real (kind=r_knd) :: tstart

#include "netcdf.inc"
      integer   ncid,iret, start(4),count(4),dims(4)

      iret=nf_open(grid_file,NF_noWRITE,ncid)
       iret=nf_inq_varid(ncid,'temp_blue',tempid)
       if (iret /=0) then
        print*,' cannot find variable temp_blue in netcdf file'
        return
       endif
       iret=nf_inq_varid(ncid,'salt_blue',saltid)
       if (iret /=0) then
        print*,' cannot find variable salt_blue in netcdf file'
        return
       endif

       iret=NF_INQ_VARDIMID (NCID,tempid,dims)
       iret=NF_INQ_DIMLEN (NCID, dims(4), itm)
       name='';iret=NF_INQ_DIMname (NCID, dims(4), name)
       iret=nf_inq_varid(ncid,name,tid)

       allocate( aprec(itm) , buf(imt,jmt,km) )

       name='';iret= nf_get_att_text(ncid,tid ,'ave_period',name)
       iret=nf_inq_varid(ncid,name,id)
#ifdef real_8
       iret= nf_get_vara_double(ncid,id ,1,itm,aprec)
#else
       iret= nf_get_vara_real(ncid,id ,1,itm,aprec)
#endif

       name='';iret= nf_get_att_text(ncid,tid ,'start_ave_period',name)
       iret=nf_inq_varid(ncid,name,id)

#ifdef real_8
       iret= nf_get_vara_double(ncid,id,1,1,tstart)
#else
       iret= nf_get_vara_real(ncid,id,1,1,tstart)
#endif
       call set_time(time,
     &                int((tstart-int(tstart))*24*60*60),int(tstart))
       call set_stamp (stamp, time)
       
       call getunit(io1,'blue_temp.mom','usr ieee')
       print*,' writing temperature to file blue_temp.mom'
       call getunit(io2,'blue_salt.mom','usr ieee')
       print*,' writing salinity to file blue_temp.mom'
       write (io1) itm,stamp,aprec
       write (io2) itm,stamp,aprec
       do n=1,itm
        start=(/1,1,1,n/); count=(/imt,jmt,km,1/); 
#ifdef real_8
        iret= nf_get_vara_double(ncid,tempid,start,count,buf)
#else
        iret= nf_get_vara_real(ncid,tempid,start,count,buf)
#endif
        do j=1,km
         where( kmt < j ) buf(:,:,j) = 0.
        enddo
        do j=1,jmt
         write(io1) buf(:,j,:)
        enddo
#ifdef real_8
        iret= nf_get_vara_double(ncid,saltid,start,count,buf)
#else
        iret= nf_get_vara_real(ncid,saltid,start,count,buf)
#endif
        do j=1,km
         where( kmt < j ) buf(:,:,j) = 0.
        enddo
        do j=1,jmt
         write(io2) buf(:,j,:)
        enddo
       enddo
       close(io1); close(io2)
       deallocate(aprec,buf)
      call ncclos (ncid, iret)
      end subroutine write_blue



      subroutine write_blue_mean()
      use prep_to_bin_mod
      implicit none
      integer :: tempid,saltid,j,io
      real (kind=r_knd), allocatable :: buf2(:,:,:), buf1(:,:,:)

#include "netcdf.inc"
      integer   ncid,iret, start(4),count(4),dims(4)

      iret=nf_open(grid_file,NF_noWRITE,ncid)

       iret=nf_inq_varid(ncid,'temp_blue_mean',tempid)
       if (iret /=0) then
        print*,' cannot find variable temp_blue_mean in netcdf file'
        return
       endif
       iret=nf_inq_varid(ncid,'salt_blue_mean',saltid)
       if (iret /=0) then
        print*,' cannot find variable salt_blue_mean in netcdf file'
        return
       endif

       allocate( buf1(imt,jmt,km) , buf2(imt,jmt,km) )
       
       call getunit(io,'blue_mean.mom','usr ieee')
       print*,' writing TS to file blue_mean.mom'
 
       start=(/1,1,1,1/); count=(/imt,jmt,km,1/); 
#ifdef real_8
       iret= nf_get_vara_double(ncid,tempid,start,count,buf1)
       iret= nf_get_vara_double(ncid,saltid,start,count,buf2)
#else
       iret= nf_get_vara_real(ncid,tempid,start,count,buf1)
       iret= nf_get_vara_real(ncid,saltid,start,count,buf2)
#endif
        do j=1,km
         where( kmt < j ) buf1(:,:,j) = 0.
         where( kmt < j ) buf2(:,:,j) = 0.
        enddo
        do j=1,jmt
         write(io) buf1(:,j,:),buf2(:,j,:)
        enddo
       close(io)
      deallocate(buf2,buf1)
      call ncclos (ncid, iret)
      end subroutine write_blue_mean





      subroutine write_ic()
      use prep_to_bin_mod
      implicit none
      character(len=80) :: vname,fname
      character(len=80) :: FLAME_ID ='no flame id '
      character(len=32+80+60) :: iotext='no stamp, no iotext '
      character(len=32) :: stamp='no stamp '
      integer :: n,j,id,i,k, mm=1, io
      real (kind=r_knd) :: dpm = 30.
      real(kind=r_knd) , allocatable :: buf(:,:)

#include "netcdf.inc"
      integer   ncid,iret, start(4),count(4),dims(4)


      allocate( buf(imt,km) )

      do n=1,99
       if (n==1) fname='temp.mom.ic'
       if (n==2) fname='salt.mom.ic'
       if (n==1) vname='temp_ic'
       if (n==2) vname='salt_ic'
       if (n>2) then
        write(fname, '("tracer_",i2,".mom.ic")') n
        do i=1,len_trim(fname)
         if (fname(i:i) == ' ') fname(i:i)='0'
        enddo
        write(vname, '("tracer_",i2,"_ic")') n
        do i=1,len_trim(vname)
         if (vname(i:i) == ' ') vname(i:i)='0'
        enddo
       endif

       if (n>2) then
        iret=nf_open(ic_tr_file,NF_noWRITE,ncid)
       else
        iret=nf_open(ic_file,NF_noWRITE,ncid)
       endif

       if (iret ==0) then
        iret=nf_inq_varid(ncid,vname,id)
        if (iret ==0) then
         print*,' writing initial conditions for tracer ',
     &          vname(1:len_trim(vname)),' to file ',
     &     fname(1:len_trim(fname))
         call getunit(io,fname,'usr ieee')
         write (io) flame_id
         do j=1,jmt
          write(io) iotext
          start=(/1,j,1,1/)
          count=(/imt,1,km,0/); 
#ifdef real_8
          iret= nf_get_vara_double(ncid,id,start,count,buf)
#else
          iret= nf_get_vara_real(ncid,id,start,count,buf)
#endif
c        apply maybe land mask
          write(io) stamp, dpm, imt,km,j, yt(j),mm,
     &         (xt(i),i=1,imt),(zt(k),k=1,km), 
     &         ((buf(i,k),i=1,imt),k=1,km)
         enddo
         close(io)
        endif
        call ncclos (ncid, iret)
       endif

      enddo ! n
      deallocate( buf )
      end subroutine write_ic





      subroutine write_topo()
      use prep_to_bin_mod
      implicit none
      character(len=80) :: FLAME_ID ='no flame id '
      character(len=32+80+60) :: iotext='no stamp, no iotext '
      integer i,j,io
      real (kind=r_knd), allocatable :: htp(:,:) 

#include "netcdf.inc"
      integer   ncid,iret, start(4),count(4),dims(4),kmtid,htpid

      iret=nf_open(grid_file,NF_noWRITE,ncid)

      iret=nf_inq_varid(ncid,'kmt',kmtid)
#ifdef partial_cell
      iret=nf_inq_varid(ncid,'htp',htpid)
#endif

      start=1;count(1)=imt;count(2)=jmt
      iret= nf_get_vara_int (ncid,kmtid,start,count,kmt)
      if (iret==0) then
      print*,' writing topography to file kmt.dta'
      call getunit(io,'kmt.dta','usr ieee')
      write (io) flame_id ! no checks
      write (io) iotext
      write (io) imt, jmt, km
      write (io) iotext
      write (io) kmt
      close (io)
      else
       print*,' cannot find topography in netcdf file'
       stop
      endif
#ifdef partial_cell
      allocate(htp(imt,jmt))
#ifdef real_8
      iret= nf_get_vara_double(ncid,htpid,start,count,htp)
#else
      iret= nf_get_vara_real(ncid,htpid,start,count,htp)
#endif
c      htp=htp_

      if (iret==0) then
      print*,' writing partial cells to file htp.dta'
      call getunit(io,'htp.dta','usr ieee')
      write (io) flame_id ! no checks
      write (io) iotext
      write (io) imt, jmt, km
      write (io) iotext
      write (io) htp
      close (io)
      deallocate(htp)
      else
       print*,' cannot find partial cells in netcdf file'
       stop
      endif
#endif
      kmu=0
      do j=1,jmt-1
        do i=1,imt-1
           kmu(i,j) = min (kmt(i,j), kmt(i+1,j), 
     &                         kmt(i,j+1), kmt(i+1,j+1))
        enddo
      enddo
      call ncclos (ncid, iret)
      end subroutine write_topo



 
      subroutine write_grid()
      use prep_to_bin_mod
      implicit none
      integer  :: i,j,k,io,iret2=0
      character(len=80) :: FLAME_ID ='no flame id '
      real (kind=8), allocatable :: v8(:)
      real (kind=r_knd), allocatable, dimension(:) :: 
     & xu,yu,zw,dzt,dzw,dxtdeg,dxudeg,dytdeg,dyudeg

#include "netcdf.inc"
      integer   ncid,iret
      integer   xtid,ytid,ztid,xuid,yuid,zwid
      integer   dxtid,dytid,dxuid,dyuid,dztid,dzwid
      integer   start(4),count(4),dims(4)

      call getunit(io,'grid.dta','usr ieee')
      print*,' writing grid informations to file grid.dta'

      write (io) flame_id   ! this the the flame ID check
      write (io) imt, jmt, km
      allocate( dxtdeg(imt), dxudeg(imt),  xu(imt))
      allocate( dytdeg(jmt), dyudeg(jmt),  yu(jmt))
      allocate( dzt(km), dzw(0:km), zw(km))
      allocate(v8(imt) )

      iret=nf_open(grid_file,NF_noWRITE,ncid)
      iret=nf_inq_varid(ncid,'xt',xtid)
      iret=nf_inq_varid(ncid,'yt',ytid)
      iret=nf_inq_varid(ncid,'zt',ztid)
      iret=nf_inq_varid(ncid,'xu',xuid)
      iret=nf_inq_varid(ncid,'yu',yuid)
      iret=nf_inq_varid(ncid,'zw',zwid)
      iret=nf_inq_varid(ncid,'dxtdeg',dxtid)
      iret=nf_inq_varid(ncid,'dytdeg',dytid)
      iret=nf_inq_varid(ncid,'dxudeg',dxuid)
      iret=nf_inq_varid(ncid,'dyudeg',dyuid)
      iret=nf_inq_varid(ncid,'dzt',dztid)
      iret=nf_inq_varid(ncid,'dzw',dzwid)

      iret= nf_get_vara_double(ncid,xtid ,1,imt,v8);xt=v8
      if (iret/=0) iret2=iret
      iret= nf_get_vara_double(ncid,xuid ,1,imt,v8);xu=v8
      if (iret/=0) iret2=iret
      iret= nf_get_vara_double(ncid,dxtid ,1,imt,v8);dxtdeg=v8
      if (iret/=0) iret2=iret
      iret= nf_get_vara_double(ncid,dxuid ,1,imt,v8);dxudeg=v8
      if (iret/=0) iret2=iret
      deallocate(v8); allocate(v8(jmt) )
      iret= nf_get_vara_double(ncid,ytid ,1,jmt,v8);yt=v8
      if (iret/=0) iret2=iret
      iret= nf_get_vara_double(ncid,yuid ,1,jmt,v8);yu=v8
      if (iret/=0) iret2=iret
      iret= nf_get_vara_double(ncid,dytid ,1,jmt,v8);dytdeg=v8
      if (iret/=0) iret2=iret
      iret= nf_get_vara_double(ncid,dyuid ,1,jmt,v8);dyudeg=v8
      if (iret/=0) iret2=iret
      deallocate(v8); allocate(v8(km) )
      iret= nf_get_vara_double(ncid,ztid ,1,km,v8);zt=v8
      if (iret/=0) iret2=iret
      iret= nf_get_vara_double(ncid,dztid ,1,km,v8);dzt=v8
      if (iret/=0) iret2=iret
      deallocate(v8); allocate(v8(km+1) )
      iret= nf_get_vara_double(ncid,zwid ,1,km+1,v8);zw=v8(2:km+1)
      if (iret/=0) iret2=iret
      iret= nf_get_vara_double(ncid,dzwid ,1,km+1,v8);
      if (iret/=0) iret2=iret
      if (iret2/=0) then
       print*,' Problem reading grid infos from  netcdf file'
       stop
      endif
      dzw(0:km)=v8(1:km+1)
      deallocate(v8)
      write (io)(dxtdeg(i),i=1,imt), (dytdeg(j),j=1,jmt)
     &,         (dxudeg(i),i=1,imt),  (dyudeg(j),j=1,jmt)
     &,         (dzt(k),k=1,km),      (dzw(k),k=0,km)
     &,         (xt(i),i=1,imt),      (xu(i),i=1,imt)
     &,         (yt(j),j=1,jmt),      (yu(j),j=1,jmt)
     &,         (zt(k),k=1,km),       (zw(k),k=1,km)
       close (io)
       print*,'done'
      call ncclos (ncid, iret)
      end subroutine write_grid

      subroutine halt_stop(string)
      character (len=*) :: string
      print *,string
      stop
      end 
