c     ------------------
c     make_cdf_ARCS_02.f
c     ------------------

      program converting_to_netcdf
c     ----------------------------      
      
c This program provides a model for converting a data set to netCDF.
c The basic strategy used in this program is to open an existing netCDF
c file, query the file for the ID's of the variables it contains, and
c then write the data to those variables.
      
c The output netCDF file must be created **before** this program is run.
c The simplest way to do this is to cd to your scratch directory and
c    % cp $FER_DIR/doc/converting_to_netcdf.basic converting_to_netcdf.cdl
c and then edit converting_to_netcdf.cdl (an ASCII file) to describe YOUR data
c set.  If your data set requires unequally spaced axes, climatological time
c axes, staggered grids, etc. then converting_to_netcdf.supplement may be
c a better starting point then the "basic" file used above.
c After you edit converting_to_netcdf.cdl then create the netCDF file with
c the command
c    % ncgen -o converting_to_netcdf.cdf converting_to_netcdf.cdl
      
c Now we will read in **your** data (gridded oceanic temperature and
c salt in this example) and write it out into the netCDF file 
c converting_to_netcdf.cdf.  Note that the axis coordinates can be written
c out exactly the same methodology - including time step values (as below).
*************************************************************************
c An alternative to modifying this program is to use the command:
      
c     ncgen -f converting_to_netcdf.cdl
      
c This will create a large source code to which select lines can 
c be added to write out your data.
*************************************************************************
c To compile and link converting_to_netcdf.f, use:
      
c     f77 -o converting_to_netcdf converting_to_netcdf.f -lnetcdf
      
*************************************************************************
c include file necessary for netCDF

c     update for your location
      
      include '/usr/local/netcdf/netcdf-3.5.0/include/netcdf.inc'    
c     may be found in $FER_DIR/fmt/cmn
*************************************************************************
c parameters relevant to the data being read in
c THESE NORMALLY MATCH THE DIMENSIONS IN THE CDL FILE
c (except nt which may be "unlimited")
      
      integer imt, jmt, km, nt, lnew, inlun
      parameter (imt=360, jmt=61, km=1, nt=240)
      
c imt is longitude, jmt latitude, km depth, and nt number of time steps 
*************************************************************************
c variable declaration
      
      real datain(imt)
      real temp(imt,jmt), taux(imt,jmt), tauy(imt,jmt)
      real hflx(imt,jmt), hc500(imt,jmt)
      real z20(imt,jmt), sst(imt,jmt), u(imt,jmt), v(imt,jmt)
      real hmxl(imt,jmt), eta(imt,jmt)

      real sum_temp(imt,jmt,12) 
      real sum_taux(imt,jmt,12)
      real sum_tauy(imt,jmt,12)
      real sum_hflx(imt,jmt,12)
      real sum_hmxl(imt,jmt,12)
      real sum_eta(imt,jmt,12)
      real sum_hc500(imt,jmt,12)
      real sum_z20(imt,jmt,12)
      real sum_u(imt,jmt,12)
      real sum_v(imt,jmt,12)
      integer int_z20(imt,jmt,12)
       
      real tmp1,tmp2,tmp3     

      real miss
      integer cdfid, rcode, index

c     netcdf field ids
      integer sstid,tauxid,tauyid,hflxid, timeaxid
      integer hc500id, z20id, uid, vid
      integer hmxlid, etaid

      integer itime,i,j,k,l
*************************************************************************
c dimension corner and step for defining size of gridded data

      integer corner1(4)
      integer step1(4)
      integer corner2(3)
      integer step2(3)
   
c corner and step are used to define the size of the gridded data 
c to be written out.  Since temp and salt are four dimensional arrays,
c corner and step must be four dimensions as well.  In each output
c to my_data.cdf within the do loop, the entire array of data (160 long. 
c pts, 100 lat. pts., 27 depth pts.) will be written for one time step.
c Corner tells netCDF where to start, and step indicates how many steps
c in each dimension to take.
      
 
      data corner1/1, 1, 1, -1/

      data corner2/1, 1, -1/
      
      data step1/imt, jmt, km, 1/
      
      data step2/imt, jmt,  1/

      real dayofdat(12)

      data dayofdat / 14,  45,  74, 105, 135, 166,
     &                196, 227, 258, 288, 319, 349 /


      
c ***NOTE*** Since Fortran and C access data differently, the order of
c the variables in the Fortran code must be opposite that in the cdl
c file.  In Fortran, the first variable varies fastest while in C, the
c last variables varies fastest.
**************************************************************************
c initialize cdfid by using ncopn
      
      cdfid = ncopn('make_cdf_61x360_grid_clim_01.cdf', ncwrite, rcode)
      if (rcode.ne.ncnoerr) stop 'error with ncopn'
      
**************************************************************************

      miss = -9.99e12
c get variable id's by using ncvid
c THE VARIABLE NAMES MUST MATCH THE CDL FILE (case sensitive)

      sstid = ncvid(cdfid, 'sst', rcode)
      if (rcode.ne.ncnoerr) stop 'error with tempid'  
      tauxid = ncvid(cdfid, 'taux', rcode)
      if (rcode.ne.ncnoerr) stop 'error with tauxid'      
      tauyid = ncvid(cdfid, 'tauy', rcode)
      if (rcode.ne.ncnoerr) stop 'error with tauyid'  
      hflxid = ncvid(cdfid, 'hflx', rcode)
      if (rcode.ne.ncnoerr) stop 'error with hflxid'  
      hmxlid = ncvid(cdfid, 'hmxl', rcode)
      if (rcode.ne.ncnoerr) stop 'error with hmxlid'  
      etaid = ncvid(cdfid, 'eta', rcode)
      if (rcode.ne.ncnoerr) stop 'error with etaid' 
      hc500id = ncvid(cdfid, 'hc500', rcode)
      if (rcode.ne.ncnoerr) stop 'error with hc500id' 
      z20id = ncvid(cdfid, 'z20', rcode)
      if (rcode.ne.ncnoerr) stop 'error with z20id' 
      uid = ncvid(cdfid, 'u', rcode)
      if (rcode.ne.ncnoerr) stop 'error with uid' 
      vid = ncvid(cdfid, 'v', rcode)
      if (rcode.ne.ncnoerr) stop 'error with vid' 

       
      timeaxid = ncvid(cdfid, 'time', rcode)
      if (rcode.ne.ncnoerr) stop 'error with timeaxid'   


**************************************************************************
c this is a good place to open your data file
      open(unit=32,
     *  file=
     *  'make_cdf_61x360_temp.dat',
     *  status='old',form='unformatted')
      open(unit=33,
     *  file=
     *  'make_cdf_61x360_taux.dat',
     *  status='old',form='unformatted')
      open(unit=34,
     *  file=
     *  'make_cdf_61x360_tauy.dat',
     *  status='old',form='unformatted')
      open(unit=37,
     *  file=
     *  'make_cdf_61x360_hc500.dat',
     *  status='old',form='unformatted')
      open(unit=38,
     *  file=
     *  'make_cdf_61x360_z20.dat',
     *  status='old',form='unformatted')
      open(unit=39,
     *  file=
     *  'make_cdf_61x360_u.dat',
     *  status='old',form='unformatted')
      open(unit=40,
     *  file=
     *  'make_cdf_61x360_v.dat',
     *  status='old',form='unformatted')
      open(unit=41,
     *  file=
     *  'make_cdf_61x360_hflx.dat',
     *  status='old',form='unformatted')
      open(unit=42,
     *  file=
     *  'make_cdf_61x360_hmxl.dat',
     *  status='old',form='unformatted')
      open(unit=43,
     *  file=
     *  'make_cdf_61x360_eta.dat',
     *  status='old',form='unformatted')


**************************************************************************
 
     
c begin do loop.  Each step will read in one time step of data
c and then write it out to my_data.cdf. 

c     zero the climatology fields
      do l=1,12
        do i=1,imt
          do j=1,jmt
            sum_taux(i,j,l) = 0.0
            sum_tauy(i,j,l) = 0.0
            sum_hflx(i,j,l) = 0.0
            sum_hmxl(i,j,l) = 0.0
            sum_eta(i,j,l) = 0.0
            sum_hc500(i,j,l) = 0.0
            sum_z20(i,j,l) = 0.0
            sum_u(i,j,l) = 0.0
            sum_v(i,j,l) = 0.0
            int_z20(i,j,l) = 0
            sum_temp(i,j,l) = 0.0
          end do
        end do
      end do

      index = 0
     
      do itime = 1, nt

        write(6,*) itime
      
        index = index +1
        if(index .eq. 13)then
          index = 1
        end if
      
c       read data
        do j=1, jmt
          read(32) datain
          do i=1, imt
            if(datain(i) .lt. -1.0e30)then
             sum_temp(i,j,index) = miss  
            else
              sum_temp(i,j,index) = sum_temp(i,j,index)
     &                              + datain(i)
            end if
          end do
        end do

        do j=1, jmt
          read(33) datain
          do i=1, imt
            if(datain(i) .lt. -1.0e30)then
              sum_taux(i,j,index) = miss
            else
              sum_taux(i,j,index) = sum_taux(i,j,index)
     &                            + datain(i)/10.0
            end if
          end do
        end do

        do j=1, jmt
          read(34) datain
          do i=1, imt
            if(datain(i) .lt. -1.0e30)then
              sum_tauy(i,j,index) = miss
            else
              sum_tauy(i,j,index) = sum_tauy(i,j,index)
     &                            + datain(i)/10.0
            end if
          end do
        end do

        do j=1, jmt
          read(41) datain
          do i=1, imt
            if(datain(i) .lt. -1.0e30)then
              sum_hflx(i,j,index) = miss
            else
              sum_hflx(i,j,index) = sum_hflx(i,j,index)
     &                            + datain(i)*4.1868e4
            end if
          end do
        end do
        do j=1, jmt
          read(42) datain
          do i=1, imt
            if(datain(i) .lt. -1.0e30)then
              sum_hmxl(i,j,index) = miss
            else
              sum_hmxl(i,j,index) = sum_hmxl(i,j,index)
     &                            + datain(i)/100.0
            end if
          end do
        end do
        do j=1, jmt
          read(43) datain
          do i=1, imt
            if(datain(i) .lt. -1.0e30)then
              sum_eta(i,j,index) = miss
            else
              sum_eta(i,j,index) = sum_eta(i,j,index)
     &                            + datain(i)
            end if
          end do
        end do
        do j=1, jmt
          read(37) datain
          do i=1, imt
            if(datain(i) .lt. -1.0e30)then
              sum_hc500(i,j,index) = miss
            else
              sum_hc500(i,j,index) = sum_hc500(i,j,index)
     &                            + datain(i)
            end if
          end do
        end do
        do j=1, jmt
          read(38) datain
          do i=1, imt
            if(datain(i) .ge. -1.0e30)then
              int_z20(i,j,index) = int_z20(i,j,index) + 1
              sum_z20(i,j,index) = sum_z20(i,j,index)
     &                            + datain(i)
            end if
          end do
        end do
        do j=1, jmt
          read(39) datain
          do i=1, imt
            if(datain(i) .lt. -1.0e30)then
              sum_u(i,j,index) = miss
            else
              u(i,j) = datain(i)
              sum_u(i,j,index) = sum_u(i,j,index)
     &                            + datain(i)/100.0
            end if
          end do
        end do
        do j=1, jmt
          read(40) datain
          do i=1, imt
            if(datain(i) .lt. -1.0e30)then
              sum_v(i,j,index) = miss
            else
              sum_v(i,j,index) = sum_v(i,j,index)
     &                            + datain(i)/100.0
            end if
          end do
        end do
      end do 

c     compute number of years
      tmp1 = float(nt)/12.0
 
      write(6,*) 'the number of years is ', tmp1
      do l=1,12
        do i=1,imt
          do j=1,jmt
            if(sum_taux(i,j,l) .ne. miss)then
              sum_taux(i,j,l)  = sum_taux(i,j,l)/tmp1
            end if
            if(sum_tauy(i,j,l) .ne. miss)then
              sum_tauy(i,j,l)  = sum_tauy(i,j,l)/tmp1
            end if
            if(sum_hflx(i,j,l) .ne. miss)then
              sum_hflx(i,j,l)  = sum_hflx(i,j,l)/tmp1
            end if
            if(sum_hmxl(i,j,l) .ne. miss)then
              sum_hmxl(i,j,l)  = sum_hmxl(i,j,l)/tmp1
            end if
            if(sum_eta(i,j,l) .ne. miss)then
              sum_eta(i,j,l)  =  sum_eta(i,j,l)/tmp1
            end if
            if(sum_hc500(i,j,l) .ne. miss)then
              sum_hc500(i,j,l) = sum_hc500(i,j,l)/tmp1
            end if
            if(sum_u(i,j,l) .ne. miss)then
              sum_u(i,j,l)     = sum_u(i,j,l)/tmp1
            end if
            if(sum_v(i,j,l) .ne. miss)then
              sum_v(i,j,l)     = sum_v(i,j,l)/tmp1
            end if
            if(int_z20(i,j,l) .ne. 0)then
              if(sum_z20(i,j,l) .ne. miss)then
                sum_z20(i,j,l) = 
     &            sum_z20(i,j,l)/tmp1
              end if
            end if
            if(sum_temp(i,j,l) .ne. miss)then
              sum_temp(i,j,l) = sum_temp(i,j,l)/tmp1
            end if
          end do
        end do
      end do

      close(32)
      close(33)
      close(34)
      close(36)
      close(37)
      close(38)
      close(39)
      close(40)
      close(41)
      close(42)
      close(43)
  
**************************************************************************
 
c begin do loop.  Each step will read in one time step of data
c and then write it out to my_data.cdf. 
       
      write(6,*) 'Data writing routine'

      index = 0      

      do itime = 1, 12

        time_step = dayofdat(itime)

        index = index +1
        if(index .eq. 13)then
          index = 1
        end if
     
c       initialize time step in corner
        corner2(3) = itime
        corner1(4) = itime
      
        do i=1,imt
          do j=1,jmt
            sst(i,j) = sum_temp(i,j,itime)
          end do
        end do

        do i=1,imt
          do j=1,jmt
            taux(i,j) = sum_taux(i,j,itime)
          end do
        end do

        do i=1,imt
          do j=1,jmt
            tauy(i,j) = sum_tauy(i,j,itime)
          end do
        end do

        do i=1,imt
          do j=1,jmt
            hflx(i,j) = sum_hflx(i,j,itime)
          end do
        end do

        do i=1,imt
          do j=1,jmt
            hmxl(i,j) = sum_hmxl(i,j,itime)
          end do
        end do

        do i=1,imt
          do j=1,jmt
            eta(i,j) = sum_eta(i,j,itime)
          end do
        end do

        do i=1,imt
          do j=1,jmt
            hc500(i,j) = sum_hc500(i,j,itime)
          end do
        end do

        do i=1,imt
          do j=1,jmt
c           use the sst mask for land since z20 has missing
c           values when it shoals
            if(sum_temp(i,j,itime) .eq. miss)then 
              z20(i,j) = miss
            else
              z20(i,j) = sum_z20(i,j,itime)
            end if
          end do
        end do

        do i=1,imt
          do j=1,jmt
            u(i,j) = sum_u(i,j,itime)
          end do
        end do

        do i=1,imt
          do j=1,jmt
            v(i,j) = sum_v(i,j,itime)
          end do
        end do      

c       writing the data to the netCDF file

        write (6,*) 'writing time step: ',itime, time_step 


        write (6,*) 'writing sst'
        call ncvpt(cdfid,sstid,corner2,step2,sst(1,1),rcode)
        if (rcode.ne.ncnoerr) stop 'error with sst-put'

        write (6,*) 'writing taux'
        call ncvpt(cdfid,tauxid,corner2,step2,taux(1,1),rcode) 
        if (rcode.ne.ncnoerr) stop 'error with taux-put'

        write (6,*) 'writing tauy'
        call ncvpt(cdfid,tauyid,corner2,step2,tauy(1,1),rcode) 
        if (rcode.ne.ncnoerr) stop 'error with tauy-put'

        write (6,*) 'writing hflx'
        call ncvpt(cdfid,hflxid,corner2,step2,hflx(1,1),rcode) 
        if (rcode.ne.ncnoerr) stop 'error with hflx-put'

        write (6,*) 'writing hmxl'
        call ncvpt(cdfid,hmxlid,corner2,step2,hmxl(1,1),rcode) 
        if (rcode.ne.ncnoerr) stop 'error with hmxl-put'

        write (6,*) 'writing eta'
        call ncvpt(cdfid,etaid,corner2,step2,eta(1,1),rcode) 
        if (rcode.ne.ncnoerr) stop 'error with eta-put'

        write (6,*) 'writing hc500'
        call ncvpt(cdfid,hc500id,corner2,step2,hc500(1,1),rcode) 
        if (rcode.ne.ncnoerr) stop 'error with hc500-put'

        write (6,*) 'writing z20'
        call ncvpt(cdfid,z20id,corner2,step2,z20(1,1),rcode) 
        if (rcode.ne.ncnoerr) stop 'error with z20-put'

        write (6,*) 'writing u'
        call ncvpt(cdfid,uid,corner2,step2,u(1,1),rcode) 
        if (rcode.ne.ncnoerr) stop 'error with u-put'

        write (6,*) 'writing v'
        call ncvpt(cdfid,vid,corner2,step2,v(1,1),rcode) 
        if (rcode.ne.ncnoerr) stop 'error with v-put'

        call ncvpt1(cdfid,timeaxid,itime,time_step,rcode) 
        if (rcode.ne.ncnoerr) stop 'error with timax-put'

      end do

      close(32)
      close(33)
      close(34)
      close(36)
      close(37)
      close(38)
      close(39)
      close(40)
      close(41)
      close(42)
      close(43)
 
**************************************************************************
c close my_data.cdf using ncclos
      call ncclos(cdfid, rcode)
      if (rcode.ne.ncnoerr) stop 'error with ncclos'
**************************************************************************
      stop
      end