PROGRAM test_spreading
! Tests spreading module by:
!
!  1. Reads a coupler history file, cplhist.nc, to get O_Runoff
!  2. Modifies runoff using runoff_spreading subroutine
!  3. Writes runoff_out.nc with the modified runoff.
!
    use spreading
    use netcdf
    implicit none
    
    ! Define variables for netcdf
    character (len = *), parameter :: INFILE = "cplhist.nc"   ! Coupler history file (input)
    character (len = *), parameter :: OFILE = "runoff_out.nc" ! Modified runoff file (output)
        
    ! Netcdf file dimensions
    integer, parameter :: NX = 362, NY = 292, NT = 240 
    real (dp), DIMENSION(NX, NY, NT) ::  runoff_in, runoff_out
    real (dp), DIMENSION(NX, NY)     ::  mask_ocn, ocn_grid_cell_area, rin, rout
    character(30) :: date
    ! This will be the netCDF ID for the file and data variable.
    integer :: ncid, varid, varid_rout
    integer :: x_dimid, y_dimid, t_dimid, dimids(3)
    ! Loop indexes
    integer :: l

! 1. READ Coupled history file data
    ! Open the file. NF90_NOWRITE tells netCDF we want read-only access to
    ! the file.
    call check( nf90_open(INFILE, NF90_NOWRITE, ncid))
    ! Get the varid of the data variable, based on its name.
    call check( nf90_inq_varid(ncid, "O_Runoff", varid))
    call check( nf90_get_var(ncid, varid, runoff_in))
!    call check( nf90_get_att(ncid, varid, '_FillValue', mv))  
    
    call check( nf90_inq_varid(ncid, "mask_ocn", varid))
    call check( nf90_get_var(ncid, varid, mask_ocn))
    
    call check( nf90_inq_varid(ncid, "ocn_grid_cell_area", varid))
    call check( nf90_get_var(ncid, varid, ocn_grid_cell_area))  
    
    call check( nf90_close(ncid))
    

! 2. Do the spreading
    DO l=1,nt
        print*, 'time index ', l
        rin=runoff_in(:,:,l)
        call runoff_spreading(rin, mask_ocn, ocn_grid_cell_area, rout)
        runoff_out(:,:,l) = rout
    ENDDO
!

! 3. Write modified runoff to NetCDF

    ! Create the output netCDF file. The nf90_clobber parameter tells netCDF to 
    ! overwrite this file, if it already exists.
    call check( nf90_create(OFILE, NF90_CLOBBER, ncid)) 

    ! Define the dimensions. NetCDF will hand back an ID for each.
    !ff= nf90_def_dim(ncid,  "time", NT, t_dimid) 
    call check( nf90_def_dim(ncid, "j", NY, y_dimid))    
    call check( nf90_def_dim(ncid, "i", NX, x_dimid)) 
    call check( nf90_def_dim(ncid, "time", NF90_UNLIMITED, t_dimid) )
    ! The dimids array is used to pass the IDs of the dimensions of
    ! the variables. Note that in fortran arrays are stored in
    ! column-major format.
    dimids =  (/ x_dimid, y_dimid, t_dimid /)

    ! define variables
    call check( nf90_def_var(ncid, "O_Runoff", NF90_DOUBLE, dimids, varid_rout))
    call check( nf90_put_att(ncid, varid_rout, "standard_name", "Runoff"))
    call check( nf90_put_att(ncid, varid_rout, "units", "kg m-2 s-1"))

    ! Some global attributes
    call check( nf90_put_att(ncid, NF90_GLOBAL, "source", "Runoff modified &
      from couopler history files."))
    call check( nf90_put_att(ncid, NF90_GLOBAL, "author", "Neil Swart &
      <Neil.Swart@canada.ca>"))      
    call fdate(date)     
    call check( nf90_put_att(ncid, NF90_GLOBAL, "history", &
      date//" ./river_spreading.exe\n"))

    ! End define mode. This tells netCDF we are done defining metadata.
    call check(nf90_enddef(ncid)) 

    ! Write the pretend data to the file. Although netCDF supports
    ! reading and writing subsets of data, in this case we write all the
    ! data in one operation.
    call check( nf90_put_var(ncid, varid_rout, runoff_out)) 

    ! Close the file. This frees up any internal netCDF resources
    ! associated with the file, and flushes any buffers.
    call check( nf90_close(ncid)) 

    print *, "*** SUCCESS writing runoff_modified nc file! "

CONTAINS
  subroutine check(status)
    integer, intent ( in) :: status
    
    if(status /= nf90_noerr) then 
      print *, trim(nf90_strerror(status))
      stop "Stopped"
    end if
  end subroutine check  


END PROGRAM test_spreading
