module kinds

  implicit none
  integer, parameter :: r8 = selected_real_kind(15, 307)

end module kinds

! This program reads a genesis grid file and adds a sizing field to it.
! It uses the new f90 interface to netCDF, although there is a cheat
! using the f77 interface to write text into the file. Please let me
! know if you fix this.
!
! The NEW flag here is to tell it whether or not this is the first time
! we are running add_size on the grid file. The first time through,
! we add the new variables, while the second time, we just have to
! write the new sizing function.
!
!                       Kate Hedstrom (kate@imcs.rutgers.edu)
!                       16 May 2000

#define NEW

program add_size
   use netcdf
   use kinds
   implicit none
#include "netcdf.inc"


   integer :: i, nnodes, times, ndims, varid, lenstr, nvars
   integer :: fileID, timedimID, nnodeID, ndimID, vardimID
   integer :: stringID
   real(r8) :: xtmp, ytmp, dist, dist2
   real(r8), allocatable :: coord(:,:)
   real(r8), allocatable :: size(:)
   real(r8)        size1, size2, size3, size4
   parameter (     size1 = 20, size2 = 40, size3 = 80, size4 = 120 )
   real(r8)  :: c_ell, b1_ell, b2_ell, b3_ell
   parameter  ( b1_ell = 500., b2_ell = 620., b3_ell = 900. )
   real(r8)  :: a1_ell, a2_ell, a3_ell
   parameter  ( a1_ell = 700., a2_ell = 790., a3_ell = 1025. )
   parameter  ( c_ell = 490.  )
   real(r8)  :: x_off, y_off
   parameter  ( x_off = 140., y_off = 200. )
   real(r8)  :: c_ell2, a1_ell_2, a2_ell_2
   parameter  ( c_ell2 = 400., a1_ell_2 = 425., a2_ell_2 = 640.  )

! open the file
   call check(nf90_open('ccs11_size.nc', nf90_write, fileID))

! read in the dimensions
   call check(nf90_inq_dimid(fileID, 'num_nodes', nnodeID))
   call check(nf90_inquire_dimension(fileID, nnodeID, len = nnodes))
   call check(nf90_inq_dimid(fileID, 'num_dim', ndimID))
   call check(nf90_inquire_dimension(fileID, ndimID, len = ndims))
   call check(nf90_inq_dimid(fileID, 'len_string', stringID))
   call check(nf90_inquire_dimension(fileID, stringID, len = lenstr))
   call check(nf90_inq_dimid(fileID, 'time_step', timedimID))
   call check(nf90_inquire_dimension(fileID, timedimID, len = times))

#ifdef NEW
! add new stuff to the file
   call check(nf90_redef(fileID))
   call check(nf90_def_dim(fileID, 'num_nod_var', 1, vardimID))
   call check(nf90_def_var(fileID, 'vals_nod_var', nf90_double, &
       (/ nnodeID, vardimID, timedimID /), varid))
   call check(nf90_def_var(fileID, 'name_nod_var', nf90_char, &
       (/ stringID, vardimID /), varid))
   call check(nf90_enddef(fileID))
#else
   call check(nf90_inq_dimid(fileID, 'num_nod_var', vardimID))
   call check(nf90_inquire_dimension(fileID, vardimID, len = nvars))
#endif /* NEW */
   print *, 'numbers = ', nnodes, ndims, lenstr, times, nvars

! compute sizing function
   allocate (coord(nnodes, ndims), size(nnodes))

   call check(nf90_inq_varid(fileID, 'coord', varid))
   call check(nf90_get_var(fileID, varid, coord))

! prepare the element size field
      do i=1,nnodes
          xtmp = coord(i,1)
          ytmp = coord(i,2)
          dist = sqrt((xtmp - x_off)**2 + (ytmp - y_off + c_ell)**2) &
               + sqrt((xtmp - x_off)**2 + (ytmp - y_off - c_ell)**2)
          dist2 = sqrt((xtmp - x_off)**2 + (ytmp - y_off + c_ell2)**2) &
                + sqrt((xtmp - x_off)**2 + (ytmp - y_off - c_ell2)**2)
          if (dist2 .le. 2*a1_ell_2) then
            size(i) = size1 
          else if (dist2 .le. 2*a2_ell_2) then
            size(i) = size1 + (size2 - size1)* &
     &             (dist2 - 2*a1_ell_2)/(2.*(a2_ell_2 - a1_ell_2))
          else if (dist .le. 2*a1_ell) then
            size(i) = size2
          else if (dist .le. 2*a2_ell) then
            size(i) = size2 + (size3 - size2)* &
     &             (dist - 2*a1_ell)/(2.*(a2_ell - a1_ell))
          else if (dist .le. 2*a3_ell) then
            size(i) = size3 + (size4 - size3)* &
     &             (dist - 2*a2_ell)/(2.*(a3_ell - a2_ell))
          else
            size(i) = size4
          endif
      enddo

#ifdef NEW
   call check(nf90_inq_varid(fileID, 'time_whole', varid))
   call check(nf90_put_var(fileID, varid, (/ 0.d0 /), &
         start = (/ 1 /), count = (/ 1 /) ))
#endif /* NEW */
   call check(nf90_inq_varid(fileID, 'vals_nod_var', varid))
   call check(nf90_put_var(fileID, varid, size, &
      start = (/ 1, 1, 1 /), count = (/ nnodes, 1, 1 /) ))
#ifdef NEW
#endif /* NEW */
   call check(nf90_inq_varid(fileID, 'name_nod_var', varid))
!  call check(nf90_put_var(fileID, varid, 'ELSIZE', &
!     start = (/ 1, 1 /), count = (/ 6, 1 /) ))
   call check(nf_put_vara_text(fileID, varid, &
      (/ 1, 1 /), (/ 6, 1 /), 'ELSIZE' ))
   call check(nf90_close(fileID))
 
 contains

  subroutine check(status)
    integer, intent(in) :: status
    if (status /= nf90_noerr) then
      print *, nf90_strerror(status)
      call exit(1)
    endif
  end subroutine check
end program add_size