module m_get_mod_grid

! The reading of depths is not needed for assimilation of sea-ice data
! but left here for potential future use. And to retain the calling format.

  use netcdf
#if defined (QMPI)
  use qmpi
#else
  use qmpi_fake
#endif

  private handle_err

contains 
  subroutine get_mod_grid(modlon,modlat,depths,mindx,meandx,nx,ny)

  implicit none

  ! In/out
  integer,                intent(in)  :: nx,ny  ! 182, resp. 149
  real, dimension(nx,ny), intent(out) :: modlon, modlat, depths
  real,                   intent(out) :: mindx, meandx

  ! NetCDF vars
  integer ncid, varID, error, ncid_mask, ncid_zgr
  character(len=80), parameter :: maskfile = 'mask.nc' !hc!
  character(len=80), parameter :: meshfile_hgr = 'mesh_hgr.nc' !hc!
  character(len=80), parameter :: meshfile_zgr = 'mesh_zgr.nc' !hc!
    
  logical ex

  ! Variables for mindx   ! uncomment whatever is needed
  real, allocatable, dimension(:,:) :: e1t, e2t!, e1u, e2u, e1v, e2v, e1f, e2f
  real, allocatable, dimension(:,:) :: tmask!, umask, vmask, fmask


  ! check the netCDF file exists
  inquire(file=meshfile_hgr, exist=ex)
  if (.not.ex) then
     if (master) print *, '(get_mod_grid): file does not exist: '//trim(meshfile_hgr)
     call stop_mpi()
  end if

  ! open the netCDF file
  error = nf90_open(trim(maskfile),nf90_NoWrite,ncid_mask)
  if (error.ne.nf90_noerr) call handle_err(error, "opening")
  error = nf90_open(trim(meshfile_hgr),nf90_NoWrite,ncid)
  if (error.ne.nf90_noerr) call handle_err(error, "opening")
  error = nf90_open(trim(meshfile_zgr),nf90_NoWrite,ncid_zgr)
  if (error.ne.nf90_noerr) call handle_err(error, "opening")

  ! Longitude
  ! Find VarID
  error = nf90_inq_varid(ncid, 'nav_lon', varID)
  if (error.ne.nf90_noerr) call handle_err(error, "inquiring varID modlon")
  ! Get values of variable
  error = nf90_get_var(ncid, varID, modlon)
  if (error.ne.nf90_noerr) call handle_err(error, "getting variable modlon")

  ! Latitude
  ! Find VarID
  error = nf90_inq_varid(ncid, 'nav_lat', varID)
  if (error.ne.nf90_noerr) call handle_err(error, "inquiring varID modlat")
  ! Get values of variable
  error = nf90_get_var(ncid, varID, modlat)
  if (error.ne.nf90_noerr) call handle_err(error, "getting variable modlat")

  ! Depths:
  ! Find VarID
  error = nf90_inq_varid(ncid_zgr, 'gdept_0', varID)
  if (error.ne.nf90_noerr) call handle_err(error, "inquiring varID depths")
  ! Get values of variable
  error = nf90_get_var(ncid_zgr, varID, depths)
  if (error.ne.nf90_noerr) call handle_err(error, "getting variable depths")

  ! mindx: Smallest horizontal grid spacing. Requires some 'math'.
  ! Load grid spacing and corresponding masks
  allocate( e1t(nx,ny), e2t(nx,ny) )
  !allocate( e1u(nx,ny), e2u(nx,ny) ) ! In case those variables are ...
  !allocate( e1v(nx,ny), e2v(nx,ny) ) ! ... needed, feel free ...
  !allocate( e1f(nx,ny), e2f(nx,ny) ) ! ... to uncomment.
  allocate( tmask(nx,ny) ) ! umask(nx,ny), vmask(nx,ny), fmask(nx,ny) 

  ! Get e1t, e1u, e1v, e1f, e2t, e2u, e2v, and e2f
  error = nf90_inq_varid(ncid, 'e1t', varID) ; if (error.ne.nf90_noerr) call handle_err(error, "inquiring varID")
  error = nf90_get_var(ncid, varID, e1t) ;     if (error.ne.nf90_noerr) call handle_err(error, "getting variable")
  error = nf90_inq_varid(ncid, 'e2t', varID) ; if (error.ne.nf90_noerr) call handle_err(error, "inquiring varID")
  error = nf90_get_var(ncid, varID, e2t) ;     if (error.ne.nf90_noerr) call handle_err(error, "getting variable")
  !error = nf90_inq_varid(ncid, 'e1u', varID) ; if (error.ne.nf90_noerr) call handle_err(error, "inquiring varID")
  !error = nf90_get_var(ncid, varID, e1u) ;     if (error.ne.nf90_noerr) call handle_err(error, "getting variable")
  !error = nf90_inq_varid(ncid, 'e2u', varID) ; if (error.ne.nf90_noerr) call handle_err(error, "inquiring varID")
  !error = nf90_get_var(ncid, varID, e2u) ;     if (error.ne.nf90_noerr) call handle_err(error, "getting variable")
  !error = nf90_inq_varid(ncid, 'e1v', varID) ; if (error.ne.nf90_noerr) call handle_err(error, "inquiring varID")
  !error = nf90_get_var(ncid, varID, e1v) ;     if (error.ne.nf90_noerr) call handle_err(error, "getting variable")
  !error = nf90_inq_varid(ncid, 'e2v', varID) ; if (error.ne.nf90_noerr) call handle_err(error, "inquiring varID")
  !error = nf90_get_var(ncid, varID, e2v) ;     if (error.ne.nf90_noerr) call handle_err(error, "getting variable")
  !error = nf90_inq_varid(ncid, 'e1f', varID) ; if (error.ne.nf90_noerr) call handle_err(error, "inquiring varID")
  !error = nf90_get_var(ncid, varID, e1f) ;     if (error.ne.nf90_noerr) call handle_err(error, "getting variable")
  !error = nf90_inq_varid(ncid, 'e2f', varID) ; if (error.ne.nf90_noerr) call handle_err(error, "inquiring varID")
  !error = nf90_get_var(ncid, varID, e2f) ;     if (error.ne.nf90_noerr) call handle_err(error, "getting variable")
  ! Get tmask, umask, vmask, fmask   !!! only first level of 3d-nc-var fits into local var. It's all we need.
  error = nf90_inq_varid(ncid_mask, 'tmask', varID) ; if (error.ne.nf90_noerr) call handle_err(error, "inquiring varID")
  error = nf90_get_var(ncid_mask, varID, tmask) ;     if (error.ne.nf90_noerr) call handle_err(error, "getting variable")
  !error = nf90_inq_varid(ncid, 'umask', varID) ; if (error.ne.nf90_noerr) call handle_err(error, "inquiring varID")
  !error = nf90_get_var(ncid, varID, umask) ;     if (error.ne.nf90_noerr) call handle_err(error, "getting variable")
  !error = nf90_inq_varid(ncid, 'vmask', varID) ; if (error.ne.nf90_noerr) call handle_err(error, "inquiring varID")
  !error = nf90_get_var(ncid, varID, vmask) ;     if (error.ne.nf90_noerr) call handle_err(error, "getting variable")
  !error = nf90_inq_varid(ncid, 'fmask', varID) ; if (error.ne.nf90_noerr) call handle_err(error, "inquiring varID")
  !error = nf90_get_var(ncid, varID, fmask) ;     if (error.ne.nf90_noerr) call handle_err(error, "getting variable")

  ! Smart use of min/maxval

  ! Find absolute minimum
  mindx = min(minval(e1t, tmask>0.5), minval(e2t, tmask>0.5))!, &
            ! minval(e1u, umask>0.5), minval(e2u, umask>0.5), &
            ! minval(e1v, vmask>0.5), minval(e2v, vmask>0.5), &
            ! minval(e1f, fmask>0.5), minval(e2f, fmask>0.5) )
  if (master) then
     print *,'(get_mod_grid) MINIMUM grid size from mesh_mask : ', mindx
  end if

  ! Find mean horizontal distance
  meandx = (sum(e1t,mask=tmask>0.5) + sum(e2t,mask=tmask>0.5) ) &
       / count(tmask>0.5)
  if (master) then
     print *,'(get_mod_grid) MEAN grid size from mesh_mask: ', meandx
  end if


  ! Safety check ..        inherited from KAL
  if (mindx<2000.) then
     if (master) print *,'(get_mod_grid) min grid size lower than safety threshold - fix if you want'
     call stop_mpi()
  end if

  ! Safety check .. This one is not that critical so the value is set high
  if (mindx>500000.) then
     if (master) print *,'(get_mod_grid) min grid size higher than safety threshold - fix if you want'
     call stop_mpi()
  end if

  ! Close file
  error = nf90_close(ncid)              ! close netCDF dataset
  if (error.ne.nf90_noerr) call handle_err(error, "closing")
  error = nf90_close(ncid_mask)              ! close netCDF dataset
  if (error.ne.nf90_noerr) call handle_err(error, "closing")
  error = nf90_close(ncid_zgr)              ! close netCDF dataset
  if (error.ne.nf90_noerr) call handle_err(error, "closing")

end subroutine  get_mod_grid

subroutine handle_err(status, infomsg) 
  integer,            intent ( in) :: status 
  character(len = *), intent ( in), optional :: infomsg
  if(status /= nf90_noerr) then
     if (master) then
        if (present(infomsg)) then
           print *, 'Error while '//infomsg//' - '//trim(nf90_strerror(status)) 
        else
           print *, trim(nf90_strerror(status)) 
        endif ! opt arg
        print *,'(get_mod_grid)'
     endif ! only master outputs
     call stop_mpi()
  end if ! check error status
end subroutine handle_err

end module  m_get_mod_grid