!---------------------------------------------------------------------------------------
! 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