#include "options.inc"

c
c-----------------------------------------------------------------------
c
c     convert NetCDF files containing the blue correction
c     to forcing files for the model run
c
c-----------------------------------------------------------------------
c
      subroutine halt_stop(string)
      character (len=*) :: string
      print *,string
      stop
      end 

      program prep_blue_const
      implicit none
      integer :: itm,imt,jmt,km,len
      real (kind=4), allocatable :: buf(:,:,:),aprec(:)
      integer :: n,io,i,k,j
      character(len=32):: stamp
      integer :: 
     &days_in_month(12) = (/31,28,31,30,31,30,31,31,30,31,30,31/)
#include "netcdf.inc"
      integer :: iret,ncid,blueid
      integer   start(4),count(4),dims(4)
      real (kind=4) :: spval
      character (len=80) :: infile,outfile

      print*,' Converting blue_const forcing file ',
     &       ' from NetCDF to binary '

      call get_commandline_arg(1,infile,iret)
      if (iret/=0) goto 99
      print*,' reading from file ',infile(1:len_trim(infile))
      call get_commandline_arg(2,outfile,iret)
      if (iret/=0) goto 99
      print*,' writing to file ',outfile(1:len_trim(outfile))

      iret=nf_open(infile, NF_nowrite,ncid)

      iret=nf_inq_varid(ncid,'blue_corr',blueid)

      if (iret /= 0 ) then
        print*,' cannotfind variable blue_corr in file'
        call halt_stop(' ERROR ')
      endif

      iret=NF_INQ_VARDIMID (NCID,blueid,dims)
      iret= nf_get_att_real (ncid,blueid ,'missing_value',spval)

      iret=NF_INQ_DIMLEN (NCID, dims(1), imt)
      iret=NF_INQ_DIMLEN (NCID, dims(2), jmt)
      iret=NF_INQ_DIMLEN (NCID, dims(3), km)
      iret=NF_INQ_DIMLEN (NCID, dims(4), len)

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

      if (len==1) then
        print*,' found an annual mean correction '
        itm=2
      elseif (len==12) then
        print*,' found a monthly mean correction '
        itm=12
      else
         print*,' do not know what average period '
         call halt_stop(' in prep_const')
      endif

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

      stamp ='m/d/y= 1/ 1/1900, h:m:s= 0: 0: 0' 

      if (len==1) then
        aprec(:)=365.
      elseif (len==12) then
        do n=1,itm
         aprec(n)=days_in_month(n)
        enddo
      else
         print*,' do not know what average period '
         call halt_stop(' in prep_const')
      endif
       
      call getunit(io,outfile,'usr ieee')
      print*,' writing blue correction '
      write (io) itm,stamp,aprec
      do n=1,itm
        print*,' time step #',n,' (well, really #',min(n,len),')'
        start=(/1,1,1,min(n,len)/); count=(/imt,jmt,km,1/); 
        iret= nf_get_vara_real(ncid,blueid,start,count,buf)
        if (iret /= 0 ) then
         print*,' problem reading variable blue_corr from file'
         call halt_stop(' ERROR ')
        endif
        where( buf == spval ) buf = 0.
        buf = buf /1000. ! convert to cgs units
        do j=1,jmt
         write(io) buf(:,j,:)
        enddo
      enddo
      close(io)
      call ncclos (ncid, iret)

      print*,' all done'

      stop

 99   print*,' program needs command line parameter'
      print*,' 1. input NetCDF file '
      print*,' 2. output binary file '

      end program prep_blue_const
