!---------------------------------------------------------------------------------------
! PREP_REGIONS: read_lsm
!---------------------------------------------------------------------------------------
!  FLEXINVERT is free software: you can redistribute it and/or modify
!  it under the terms of the GNU General Public License as published by
!  the Free Software Foundation, either version 3 of the License, or
!  (at your option) any later version.
!
!  FLEXINVERT is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!  GNU General Public License for more details.
!
!  You should have received a copy of the GNU General Public License
!  along with FLEXINVERT.  If not, see <http://www.gnu.org/licenses/>.
!
!  Copyright 2017, Rona Thompson
!---------------------------------------------------------------------------------------
!
!> read_lsm
!!
!! Purpose:    Reads the land-sea mask from ncdf file.
!! 
!---------------------------------------------------------------------------------------

subroutine read_lsm(files)

  use mod_var
  use mod_settings
  use mod_ncdf,    only : check
  use netcdf

  implicit none

  type (files_t),                 intent(in)     :: files

  character(len=max_path_len)                    :: filename
  character(len=max_name_len)                    :: varname, lonname, latname
  logical                                        :: lexist
  integer                                        :: varid, lonid, latid
  integer                                        :: nlon, nlat
  integer                                        :: ncid, lon_dimid, lat_dimid
  integer                                        :: ierr
  integer                                        :: ix, jy
  real, dimension(:), allocatable                :: lon, lat
  real, dimension(:,:,:), allocatable            :: ice
  real, dimension(:,:), allocatable              :: ice_av

  ! read lsm file
  ! -------------

  allocate ( lsm(nxregrid,nyregrid), stat=ierr )
  if ( ierr.ne.0 ) stop 'ERROR read_lsm: not enough memory'

  filename = files%file_lsm
  varname = files%varname_lsm
  lonname = files%lonname_lsm
  latname = files%latname_lsm

  inquire(file=trim(filename),exist=lexist)
  if (.not.lexist) then
    write(*,*) 'ERROR: cannot find '//trim(filename)
    stop
  endif

  write(*,*) 'Reading file: '//trim(filename)

  call check( nf90_open(trim(filename),nf90_NOWRITE,ncid) )
  call check( nf90_inq_dimid(ncid,trim(lonname),lon_dimid) )
  call check( nf90_inq_dimid(ncid,trim(latname),lat_dimid) )

  call check( nf90_inquire_dimension(ncid,lon_dimid,len=nlon) )
  call check( nf90_inquire_dimension(ncid,lat_dimid,len=nlat) )

  print*, 'read_lsm: nlon, nlat =', nlon, nlat

  allocate ( lon(nlon) )
  allocate ( lat(nlat) )

  call check( nf90_inq_varid(ncid,trim(lonname),lonid) )
  call check( nf90_get_var(ncid,lonid,lon) )
  call check( nf90_inq_varid(ncid,trim(latname),latid) )
  call check( nf90_get_var(ncid,latid,lat) )

  ! check latitude direction
  if ( lat(1).gt.lat(2) ) then
    write(*,*) 'ERROR: lsm latitude must be increasing'
    stop
  endif
    
  ! indices to variable
  ix = minloc(lon, dim=1, mask=abs(lon-rdx/2.-rllx).eq.minval(abs(lon-rdx/2.-rllx)))
  jy = minloc(lat, dim=1, mask=abs(lat-rdy/2.-rlly).eq.minval(abs(lat-rdy/2.-rlly)))

  print*, 'read_lsm: lsm indices lon(ix), lat(jy) = ',lon(ix),lat(jy)

  ! check lsm extent
  if ( (ix.eq.0).or.(jy.eq.0) ) then
    write(*,*) 'ERROR: incompatible domain in '//trim(filename)
    stop
  endif
  if ( (nlon.lt.(ix+nxregrid-1)).or.(nlat.lt.(jy+nyregrid-1)) ) then
    write(*,*) 'ERROR: incompatible dimensions in '//trim(filename)
    stop
  endif

  ! read lsm
  call check( nf90_inq_varid(ncid,trim(varname),varid) )
  call check( nf90_get_var(ncid,varid,lsm,start=(/ix,jy/),count=(/nxregrid,nyregrid/)) )
  call check( nf90_close(ncid) )

  deallocate(lon)
  deallocate(lat)

  ! read ice file
  ! --------------

  ! only if file_ice is specified
  if ( len_trim(files%file_ice).gt.0 ) then

    allocate ( ice(nxregrid,nyregrid,12), stat=ierr )
    if ( ierr.ne.0 ) stop 'ERROR read_lsm: not enough memory'
    allocate ( ice_av(nxregrid,nyregrid), stat=ierr )
    if ( ierr.ne.0 ) stop 'ERROR read_lsm: not enough memory'

    filename = files%file_ice
    varname = files%varname_ice
    lonname = files%lonname_ice
    latname = files%latname_ice

    inquire(file=trim(filename),exist=lexist)
    if (.not.lexist) then
      write(*,*) 'ERROR: cannot find '//trim(filename)
      stop
    endif

    write(*,*) 'Reading file: '//trim(filename)

    call check( nf90_open(trim(filename),nf90_NOWRITE,ncid) )
    call check( nf90_inq_dimid(ncid,trim(lonname),lon_dimid) )
    call check( nf90_inq_dimid(ncid,trim(latname),lat_dimid) )

    call check( nf90_inquire_dimension(ncid,lon_dimid,len=nlon) )
    call check( nf90_inquire_dimension(ncid,lat_dimid,len=nlat) )

    print*, 'read_lsm: nlon, nlat =', nlon, nlat

    allocate ( lon(nlon) )
    allocate ( lat(nlat) )

    call check( nf90_inq_varid(ncid,trim(lonname),lonid) )
    call check( nf90_get_var(ncid,lonid,lon) )
    call check( nf90_inq_varid(ncid,trim(latname),latid) )
    call check( nf90_get_var(ncid,latid,lat) )

    ! indices to variable
    ix = minloc(lon, dim=1, mask=abs(lon-rdx/2.-rllx).eq.minval(abs(lon-rdx/2.-rllx)))
    jy = minloc(lat, dim=1, mask=abs(lat-rdy/2.-rlly).eq.minval(abs(lat-rdy/2.-rlly)))

    print*, 'read_lsm: ice indices lon(ix), lat(jy) = ',lon(ix),lat(jy)

    ! check ice extent
    if ( (ix.eq.0).or.(jy.eq.0) ) then
      write(*,*) 'ERROR: incompatible domain in '//trim(filename)
      stop
    endif
    if ( (nlon.lt.(ix+nxregrid-1)).or.(nlat.lt.(jy+nyregrid-1)) ) then
      write(*,*) 'ERROR: incompatible dimensions in '//trim(filename)
      stop
    endif

    call check( nf90_inq_varid(ncid,trim(varname),varid) )
    call check( nf90_get_var(ncid,varid,ice,start=(/ix,jy,1/),count=(/nxregrid,nyregrid,12/)) )
    call check( nf90_close(ncid) )

    ! integrate ice into lsm
    ice_av = sum(ice,dim=3)/12.
    where( ice_av(:,:).gt.0.99 ) lsm = 0.

    deallocate(lon)
    deallocate(lat)
    deallocate(ice)
    deallocate(ice_av)

  endif

end subroutine read_lsm

