!===============================================================================
! SVN $Id: shr_map_mod.F90 22686 2010-05-03 16:30:51Z tcraig $
! SVN $URL: https://svn-ccsm-models.cgd.ucar.edu/csm_share/branch_tags/cesm1_0_rel_tags/cesm1_0_rel01_share3_100616/shr/shr_map_mod.F90 $
!===============================================================================
!BOP ===========================================================================
!
! !MODULE: shr_map_mod -- generic map data type and associated methods
!
! !DESCRIPTION:
! Generic map data type and associated methods
! \newline
! This module supports mapping of fields from one grid to another.
! A general datatype, shr\_map\_mapType, stores the mapping information
! set in shr\_map\_mapSet. shr\_map\_mapData then allows this mapping
! to be applied to an input array to generate the output array.
! \newline
! The mapType has several flags that give the user various options
! for setting the mapping
! type: [remap,fill]
! remap - mapping of data between different grids, primarily
! for the active grid area
! fill - mapping of data on the same grid, primarily to fill missing
! areas, copy data, or set the array to a spval.
! algo: [copy,bilinear,nn,nnoni,nnonj,spval]
! copy - copy data from one array to another using indexing
! bilinear - bilinear remapping using 4 corner points
! nn - nearest neighbor, set value to nn value
! nnoni - nearest neighbor using i, search for nearest neighbor in the
! i direction first, then j
! nnonj - nearest neighbor using j, search for nearest neighbor in the
! j direction first, then i
! spval - set values to the spval
! mask: [srcmask,dstmask,nomask,bothmask]
! srcmask - use only src points with mask = true in mapping
! dstmask - map only to dst points where mask = true
! nomask - ignore both src and dst mask in mapping
! bothmask - use both src and dst mask in mapping (srcmask and dstmask)
! vect: [scalar,vector]
! scalar - fields are scalar type (default)
! vector - fields are vector type, operates only on 2 fields to 2 fields
! NOTE: Not all combinatations are unique and not all combinations are valid
! \newline
! The above settings are put into the maptype using shr\_map\_put. Public
! parameters are available to users to set the switches. The first three
! switches must be set then the mapSet method can be called. After the
! mapSet method is called, the mapData method can be used.
! \newline
! A Note on Subroutine Arguments:
! Lat, lon, and mask arguments in these routines are 2d (nx,ny)
! Array arguments are 2d (nf,nxy), number of fields by grid point
! \newline
! General Usage:
! type(shr\_map\_mapType) :: mymap
! call shr\_map\_put(mymap,'type','remap')
! call shr\_map\_put(mymap,shr\_map\_fs\_algo,shr\_map\_fs\_bilinear)
! call shr\_map\_put(mymap,shr\_map\_fs\_mask,'bothmask')
! call shr\_map\_put(mymap,shr\_map\_fs\_vect,'scalar')
! call shr\_map\_mapSet(mymap,Xsrc,Ysrc,Msrc,Xdst,Ydst,Mdst,rc=rCode)
! call shr\_map\_mapData(Asrc,Adst,mymap)
! \newline
! call shr\_map\_mapSet(mymap,Xs,Ys,Ms,Xd,Yd,Md,name='fillnnoni',type='fill',algo='nnoni',mask='dstmask',rc=rc)
! call shr\_map\_mapData(Asrc,Adst,mymap)
! \newline
! call shr\_map\_mapData(Ain,Aout,Xs,Ys,Ms,Xd,Yd,Md,type='remap',algo='nn',mask='dstmask',rc)
!
! !REMARKS:
! nn needs a faster algorithm
!
! !REVISION HISTORY:
! 2005-Mar-27 - T. Craig - first version
!
! !INTERFACE: ------------------------------------------------------------------
module shr_map_mod 5,6
! !USES:
use shr_const_mod
use shr_kind_mod
use shr_sys_mod
use shr_timer_mod
use shr_log_mod, only: s_loglev => shr_log_Level
use shr_log_mod, only: s_logunit => shr_log_Unit
implicit none
private
! !PUBLIC TYPES:
public :: shr_map_maptype ! shr_map datatype
type shr_map_mapType ! like mct sparsematrix datatype
private
character(SHR_KIND_CS) :: name
character(SHR_KIND_CS) :: type
character(SHR_KIND_CS) :: algo
character(SHR_KIND_CS) :: mask
character(SHR_KIND_CS) :: vect
integer(SHR_KIND_IN) :: nsrc ! grid size or src
integer(SHR_KIND_IN) :: ndst ! grid size of dst
integer(SHR_KIND_IN) :: nwts ! number of total weights
real(SHR_KIND_R8) ,pointer :: xsrc(:) ! longitude, for vector, rad
real(SHR_KIND_R8) ,pointer :: ysrc(:) ! latitude , for vector, rad
real(SHR_KIND_R8) ,pointer :: xdst(:) ! longitude, for vector, rad
real(SHR_KIND_R8) ,pointer :: ydst(:) ! latitude , for vector, rad
real(SHR_KIND_R8) ,pointer :: wgts(:) ! weights
integer(SHR_KIND_IN),pointer :: isrc(:) ! input grid index
integer(SHR_KIND_IN),pointer :: idst(:) ! output grid index
character(SHR_KIND_CS) :: fill ! string to check if filled
character(SHR_KIND_CS) :: init ! initialization of dst array
end type shr_map_mapType
! PUBLIC MEMBER FUNCTIONS:
public :: shr_map_checkInit ! check whether map type is set
public :: shr_map_checkFilled ! check whether map wts are set
public :: shr_map_put ! put stuff into the datatype
public :: shr_map_get ! get stuff out of the datatype
public :: shr_map_mapSet ! compute weights in map
public :: shr_map_mapData ! map data
public :: shr_map_listValidOpts ! list valid options
public :: shr_map_print ! print map datatype info
public :: shr_map_clean ! clean map datatype
public :: shr_map_setAbort ! set abort flag for shr_map
public :: shr_map_setDebug ! set debug level for shr_map
public :: shr_map_setDopole ! set dopole flag
! PUBLIC DATA MEMBERS:
!--- Field Strings (fldStr) ---
character(SHR_KIND_CS),public,parameter :: shr_map_fs_name = 'name'
character(SHR_KIND_CS),public,parameter :: shr_map_fs_type = 'type'
character(SHR_KIND_CS),public,parameter :: shr_map_fs_algo = 'algo'
character(SHR_KIND_CS),public,parameter :: shr_map_fs_mask = 'mask'
character(SHR_KIND_CS),public,parameter :: shr_map_fs_vect = 'vect'
character(SHR_KIND_CS),public,parameter :: shr_map_fs_nwts = 'nwts'
character(SHR_KIND_CS),public,parameter :: shr_map_fs_nsrc = 'nsrc'
character(SHR_KIND_CS),public,parameter :: shr_map_fs_ndst = 'ndst'
!--- "type" options ---
character(len=*),public,parameter :: shr_map_fs_fill = 'fill '
character(len=*),public,parameter :: shr_map_fs_cfill = 'cfill '
character(len=*),public,parameter :: shr_map_fs_remap = 'remap '
!--- "algorithm" options ---
character(len=*),public,parameter :: shr_map_fs_copy = 'copy '
character(len=*),public,parameter :: shr_map_fs_bilinear = 'bilinear'
character(len=*),public,parameter :: shr_map_fs_nn = 'nn '
character(len=*),public,parameter :: shr_map_fs_nnoni = 'nnoni '
character(len=*),public,parameter :: shr_map_fs_nnonj = 'nnonj '
character(len=*),public,parameter :: shr_map_fs_spval = 'spval '
!--- "mask" options ---
character(len=*),public,parameter :: shr_map_fs_srcmask = 'srcmask '
character(len=*),public,parameter :: shr_map_fs_dstmask = 'dstmask '
character(len=*),public,parameter :: shr_map_fs_nomask = 'nomask '
character(len=*),public,parameter :: shr_map_fs_bothmask = 'bothmask'
!--- "vect" options ---
character(len=*),public,parameter :: shr_map_fs_scalar = 'scalar '
character(len=*),public,parameter :: shr_map_fs_vector = 'vector '
!--- other public parameters ---
character(SHR_KIND_CS),public,parameter :: shr_map_setTru = 'TRUE map'
character(SHR_KIND_CS),public,parameter :: shr_map_setFal = 'FALSE m '
integer(SHR_KIND_IN) ,public,parameter :: shr_map_ispval = -99
real(SHR_KIND_R8) ,public,parameter :: shr_map_spval = shr_const_spval
!EOP
!--- Must update these if anything above changes ---
integer(SHR_KIND_IN),public,parameter :: shr_map_fs_ntype = 3
character(len=*),public,parameter :: &
shr_map_fs_types(shr_map_fs_ntype) = (/shr_map_fs_fill, &
shr_map_fs_cfill, &
shr_map_fs_remap /)
integer(SHR_KIND_IN),public,parameter :: shr_map_fs_nalgo = 6
character(len=*),public,parameter :: &
shr_map_fs_algos(shr_map_fs_nalgo) = (/shr_map_fs_copy, &
shr_map_fs_bilinear, &
shr_map_fs_nn, &
shr_map_fs_nnoni, &
shr_map_fs_nnonj, &
shr_map_fs_spval /)
integer(SHR_KIND_IN),public,parameter :: shr_map_fs_nmask = 4
character(len=*),public,parameter :: &
shr_map_fs_masks(shr_map_fs_nmask) = (/shr_map_fs_srcmask, &
shr_map_fs_dstmask, &
shr_map_fs_nomask , &
shr_map_fs_bothmask /)
integer(SHR_KIND_IN),public,parameter :: shr_map_fs_nvect = 2
character(len=*),public,parameter :: &
shr_map_fs_vects(shr_map_fs_nvect) = (/shr_map_fs_scalar, &
shr_map_fs_vector /)
interface shr_map_put ; module procedure &
shr_map_putCS, &
shr_map_putR8, &
shr_map_putIN
end interface
interface shr_map_get ; module procedure &
shr_map_getCS, &
shr_map_getR8, &
shr_map_getIN, &
shr_map_getAR
end interface
interface shr_map_mapSet ; module procedure &
shr_map_mapSet_global, &
shr_map_mapSet_dest
end interface
interface shr_map_mapData ; module procedure &
shr_map_mapDatam, &
shr_map_mapDatanm
end interface
logical,save :: doabort = .true.
logical,save :: dopole = .true. ! for bilinear
integer(SHR_KIND_IN),save :: debug = 0
character(SHR_KIND_CS),parameter :: fillstring = 'mapisfilled'
character(SHR_KIND_CS),parameter :: inispval = 'spval'
character(SHR_KIND_CS),parameter :: initcopy = 'copy'
real(SHR_KIND_R8) ,parameter :: c0 = 0._SHR_KIND_R8
real(SHR_KIND_R8) ,parameter :: c1 = 1._SHR_KIND_R8
real(SHR_KIND_R8) ,parameter :: c2 = 2._SHR_KIND_R8
real(SHR_KIND_R8) ,parameter :: eps = 1.0e-12_SHR_KIND_R8
real(SHR_KIND_R8) ,parameter :: pi = shr_const_pi
!===============================================================================
contains
!===============================================================================
!BOP ===========================================================================
!
! !IROUTINE: shr_map_checkInit -- returns init state of map
!
! !DESCRIPTION:
! Returns init state of map. shr\_map\_checkInit is true
! if the type, algo, and mask are set to valid values.
! \newline
! test = shr\_map\_checkInit(map)
!
! !REVISION HISTORY:
! 2005-Mar-27 - T. Craig - first version
!
! !INTERFACE: ------------------------------------------------------------------
logical function shr_map_checkInit(map) 1,3
implicit none
! !INPUT/OUTPUT PARAMETERS:
type(shr_map_mapType),intent(in) :: map
!EOP
!--- local ---
!--- formats ---
character(*),parameter :: subName = "('shr_map_checkInit') "
!-------------------------------------------------------------------------------
if (shr_map_checkFldStrOpt(shr_map_fs_type,map%type) .and. &
shr_map_checkFldStrOpt(shr_map_fs_algo,map%algo) .and. &
shr_map_checkFldStrOpt(shr_map_fs_mask,map%mask)) then
shr_map_checkInit = .true.
else
shr_map_checkInit = .false.
endif
end function shr_map_checkInit
!===============================================================================
!BOP ===========================================================================
!
! !IROUTINE: shr_map_checkFilled -- returns fill state of map
!
! !DESCRIPTION:
! Returns fill state of map. shr\_map\_checkFilled is true
! if the number of weights are greater than zero in map
! and if the wgts, isrc, and idst arrays have been allocated to
! that size.
! \newline
! test = shr\_map\_checkFilled(map)
!
! !REVISION HISTORY:
! 2005-Mar-27 - T. Craig - first version
!
! !INTERFACE: ------------------------------------------------------------------
logical function shr_map_checkFilled(map) 1
implicit none
! !INPUT/OUTPUT PARAMETERS:
type(shr_map_mapType),intent(in) :: map
!EOP
!--- local ---
integer(SHR_KIND_IN) :: nwts
!--- formats ---
character(*),parameter :: subName = "('shr_map_checkFilled') "
!-------------------------------------------------------------------------------
shr_map_checkFilled = .false.
nwts = map%nwts
if (map%fill == fillstring .and. nwts >= 0) then
if (size(map%wgts) == nwts .and. size(map%isrc) == nwts &
.and. size(map%idst) == nwts ) then
shr_map_checkFilled = .true.
endif
endif
end function shr_map_checkFilled
!===============================================================================
!XXBOP ===========================================================================
!
! !IROUTINE: shr_map_checkFldStr -- checks fldstr for validity
!
! !DESCRIPTION:
! Returns true if fldstr is valid (ie. 'type','algo','mask')
! \newline
! test = shr\_map\_checkFldStr('type')
!
! !REVISION HISTORY:
! 2005-Mar-27 - T. Craig - first version
!
! !INTERFACE: ------------------------------------------------------------------
logical function shr_map_checkFldStr(fldStr)
implicit none
! !INPUT/OUTPUT PARAMETERS:
character(*) :: fldStr
!XXEOP
!--- local ---
!--- formats ---
character(*),parameter :: subName = "('shr_map_checkFldStr') "
!-------------------------------------------------------------------------------
shr_map_checkFldStr = .false.
if (trim(fldStr) == trim(shr_map_fs_type).or. &
trim(fldStr) == trim(shr_map_fs_name).or. &
trim(fldStr) == trim(shr_map_fs_algo).or. &
trim(fldStr) == trim(shr_map_fs_mask).or. &
trim(fldStr) == trim(shr_map_fs_vect).or. &
trim(fldStr) == trim(shr_map_fs_nsrc).or. &
trim(fldStr) == trim(shr_map_fs_ndst).or. &
trim(fldStr) == trim(shr_map_fs_nwts)) then
shr_map_checkFldStr = .true.
endif
end function shr_map_checkFldStr
!===============================================================================
!XXBOP ===========================================================================
!
! !IROUTINE: shr_map_checkFldStrOpt -- checks cval for validity with fldstr
!
! !DESCRIPTION:
! Returns true if cval is valid for fldstr (ie. 'type,remap','algo,bilinear',
! 'mask,srcmask')
! \newline
! test = shr\_map\_checkFldStrOpt(shr_map_fs_algo,'bilinear')
!
! !REVISION HISTORY:
! 2005-Mar-27 - T. Craig - first version
!
! !INTERFACE: ------------------------------------------------------------------
logical function shr_map_checkFldStrOpt(fldStr,cval) 3
implicit none
! !INPUT/OUTPUT PARAMETERS:
character(*),intent(in) :: fldStr
character(*),intent(in) :: cval
!XXEOP
!--- local ---
integer(SHR_KIND_IN) :: n
!--- formats ---
character(*),parameter :: subName = "('shr_map_checkFldStrOpt') "
!-------------------------------------------------------------------------------
shr_map_checkFldStrOpt = .false.
if (.not.shr_map_checkFldStr(fldStr)) return
if (trim(fldStr) == trim(shr_map_fs_name)) then
shr_map_checkFldStrOpt = .true.
elseif (trim(fldStr) == trim(shr_map_fs_type)) then
do n = 1,shr_map_fs_ntype
if (trim(cval) == trim(shr_map_fs_types(n))) shr_map_checkFldStrOpt = .true.
enddo
elseif (trim(fldStr) == trim(shr_map_fs_algo)) then
do n = 1,shr_map_fs_nalgo
if (trim(cval) == trim(shr_map_fs_algos(n))) shr_map_checkFldStrOpt = .true.
enddo
elseif (trim(fldStr) == trim(shr_map_fs_mask)) then
do n = 1,shr_map_fs_nmask
if (trim(cval) == trim(shr_map_fs_masks(n))) shr_map_checkFldStrOpt = .true.
enddo
elseif (trim(fldStr) == trim(shr_map_fs_vect)) then
do n = 1,shr_map_fs_nvect
if (trim(cval) == trim(shr_map_fs_vects(n))) shr_map_checkFldStrOpt = .true.
enddo
endif
end function shr_map_checkFldStrOpt
!===============================================================================
!BOP ===========================================================================
!
! !IROUTINE: shr_map_getCS -- get string from map
!
! !DESCRIPTION:
! one of the shr\_map\_get methods for chars
! returns value cval for input fldstr in map
! \newline
! call shr\_map\_get(mymap,shr\_map\_fs\_type,cval)
!
! !REVISION HISTORY:
! 2005-Mar-27 - T. Craig - first version
!
! !INTERFACE: ------------------------------------------------------------------
subroutine shr_map_getCS(map,fldStr,cval),2
implicit none
! !INPUT/OUTPUT PARAMETERS:
type(shr_map_mapType) ,intent(in) :: map
character(*) ,intent(in) :: fldStr
character(*) ,intent(out):: cval
!EOP
!--- local ---
!--- formats ---
character(*),parameter :: subName = "('shr_map_getCS') "
!-------------------------------------------------------------------------------
cval = shr_map_setFal
if (.not.shr_map_checkFldStr(fldStr)) then
call shr_map_abort(subName//' ERROR illegal fldStr '//trim(fldStr))
return
endif
if (trim(fldStr) == trim(shr_map_fs_name)) then
cval = map%name
elseif (trim(fldStr) == trim(shr_map_fs_type)) then
cval = map%type
elseif (trim(fldStr) == trim(shr_map_fs_algo)) then
cval = map%algo
elseif (trim(fldStr) == trim(shr_map_fs_mask)) then
cval = map%mask
elseif (trim(fldStr) == trim(shr_map_fs_vect)) then
cval = map%vect
else
call shr_map_abort(subName//' ERROR illegal fldStr '//trim(fldStr))
endif
end subroutine shr_map_getCS
!===============================================================================
!BOP ===========================================================================
!
! !IROUTINE: shr_map_getIN -- get integer from map
!
! !DESCRIPTION:
! one of the shr\_map\_get methods for integers
! returns value ival for input fldstr in map
! \newline
! call shr\_map\_get(mymap,shr\_map\_fs\_nwts,ival)
!
! !REVISION HISTORY:
! 2005-Mar-27 - T. Craig - first version
!
! !INTERFACE: ------------------------------------------------------------------
subroutine shr_map_getIN(map,fldStr,ival),2
implicit none
! !INPUT/OUTPUT PARAMETERS:
type(shr_map_mapType) ,intent(in) :: map
character(*) ,intent(in) :: fldStr
integer(SHR_KIND_IN) ,intent(out):: ival
!EOP
!--- local ---
!--- formats ---
character(*),parameter :: subName = "('shr_map_getIN') "
!-------------------------------------------------------------------------------
ival = shr_map_ispval
if (.not.shr_map_checkFldStr(fldStr)) then
call shr_map_abort(subName//' ERROR illegal fldStr '//trim(fldStr))
return
endif
if (trim(fldStr) == trim(shr_map_fs_nwts)) then
ival = map%nwts
elseif (trim(fldStr) == trim(shr_map_fs_nsrc)) then
ival = map%nsrc
elseif (trim(fldStr) == trim(shr_map_fs_ndst)) then
ival = map%ndst
else
call shr_map_abort(subName//' ERROR illegal fldStr '//trim(fldStr))
endif
end subroutine shr_map_getIN
!===============================================================================
!BOP ===========================================================================
!
! !IROUTINE: shr_map_getR8 -- get real from map
!
! !DESCRIPTION:
! one of the shr\_map\_get methods for reals
! returns value rval for input fldstr in map
!
! !REVISION HISTORY:
! 2005-Mar-27 - T. Craig - first version
!
! !INTERFACE: ------------------------------------------------------------------
subroutine shr_map_getR8(map,fldStr,rval),2
implicit none
! !INPUT/OUTPUT PARAMETERS:
type(shr_map_mapType) ,intent(in) :: map
character(*) ,intent(in) :: fldStr
real(SHR_KIND_R8) ,intent(out):: rval
!EOP
!--- local ---
!--- formats ---
character(*),parameter :: subName = "('shr_map_getR8') "
!-------------------------------------------------------------------------------
rval = shr_map_spval
if (.not.shr_map_checkFldStr(fldStr)) then
call shr_map_abort(subName//' ERROR illegal fldStr '//trim(fldStr))
return
endif
call shr_map_abort(subName//' ERROR illegal fldStr '//trim(fldStr))
end subroutine shr_map_getR8
!===============================================================================
!BOP ===========================================================================
!
! !IROUTINE: shr_map_getAR -- get arrays from map
!
! !DESCRIPTION:
! one of the shr\_map\_get methods for arrays
! returns value ival for input fldstr in map
! \newline
! call shr\_map\_get(mymap,idst,isrc,wgts)
!
! !REVISION HISTORY:
! 2009-Jul-27 - T. Craig - first version
!
! !INTERFACE: ------------------------------------------------------------------
subroutine shr_map_getAR(map,isrc,idst,wgts),3
implicit none
! !INPUT/OUTPUT PARAMETERS:
type(shr_map_mapType) ,intent(in) :: map
integer(SHR_KIND_IN),pointer,optional :: isrc(:)
integer(SHR_KIND_IN),pointer,optional :: idst(:)
real (SHR_KIND_R8),pointer,optional :: wgts(:)
!EOP
!--- local ---
integer(SHR_KIND_IN) :: nwts
!--- formats ---
character(*),parameter :: subName = "('shr_map_getAR') "
!-------------------------------------------------------------------------------
nwts = map%nwts
if (present(isrc)) then
if (size(isrc) < nwts) then
call shr_sys_abort(subName//' ERROR is isrc size')
endif
isrc(1:nwts) = map%isrc(1:nwts)
endif
if (present(idst)) then
if (size(idst) < nwts) then
call shr_sys_abort(subName//' ERROR is idst size')
endif
idst(1:nwts) = map%idst(1:nwts)
endif
if (present(wgts)) then
if (size(wgts) < nwts) then
call shr_sys_abort(subName//' ERROR is wgts size')
endif
wgts(1:nwts) = map%wgts(1:nwts)
endif
end subroutine shr_map_getAR
!===============================================================================
!BOP ===========================================================================
!
! !IROUTINE: shr_map_putCS -- put char to map
!
! !DESCRIPTION:
! one of the shr\_map\_put methods for chars
! puts value cval for input fldstr in map
! verify is optional argument that check validity and will
! call abort if cval is not valid option for fldstr.
! \newline
! call shr\_map\_put(mymap,shr\_map\_fs\_algo,shr\_map\_fs\_bilinear)
!
! !REVISION HISTORY:
! 2005-Mar-27 - T. Craig - first version
!
! !INTERFACE: ------------------------------------------------------------------
subroutine shr_map_putCS(map,fldStr,cval,verify),2
implicit none
! !INPUT/OUTPUT PARAMETERS:
type(shr_map_mapType) ,intent(inout):: map
character(*) ,intent(in) :: fldStr
character(*) ,intent(in) :: cval
logical,optional ,intent(in) :: verify ! check if string is valid
!EOP
!--- local ---
logical :: lverify
!--- formats ---
character(*),parameter :: subName = "('shr_map_putCS') "
!-------------------------------------------------------------------------------
lverify = .true.
if (present(verify)) lverify = verify
if (lverify .and. .not.shr_map_checkFldStrOpt(fldStr,cval)) then
call shr_map_abort(subName//' ERROR illegal fldStr '//trim(fldStr)//' '//trim(cval))
return
endif
if (trim(fldStr) == trim(shr_map_fs_name)) then
map%name = cval
elseif (trim(fldStr) == trim(shr_map_fs_type)) then
map%type = cval
elseif (trim(fldStr) == trim(shr_map_fs_algo)) then
map%algo = cval
elseif (trim(fldStr) == trim(shr_map_fs_mask)) then
map%mask = cval
elseif (trim(fldStr) == trim(shr_map_fs_vect)) then
map%vect = cval
else
call shr_map_abort(subName//' ERROR illegal fldStr '//trim(fldStr))
endif
end subroutine shr_map_putCS
!===============================================================================
!BOP ===========================================================================
!
! !IROUTINE: shr_map_putIN -- put integer to map
!
! !DESCRIPTION:
! one of the shr\_map\_put methods for integers
! puts value ival for input fldstr in map
! verify is optional argument that check validity and will
! call abort if ival is not valid option for fldstr.
! \newline
! call shr\_map\_put(mymap,shr\_map\_fs\_nwts,-1)
!
! !REVISION HISTORY:
! 2005-Mar-27 - T. Craig - first version
!
! !INTERFACE: ------------------------------------------------------------------
subroutine shr_map_putIN(map,fldStr,ival,verify),2
implicit none
! !INPUT/OUTPUT PARAMETERS:
type(shr_map_mapType) ,intent(inout):: map
character(*) ,intent(in) :: fldStr
integer(SHR_KIND_IN) ,intent(in) :: ival
logical,optional ,intent(in) :: verify ! check if string is valid
!EOP
!--- local ---
logical :: lverify
!--- formats ---
character(*),parameter :: subName = "('shr_map_putIN') "
character(*),parameter :: F01 = "('(shr_map_putIN) ',a,i8) "
!-------------------------------------------------------------------------------
lverify = .true.
if (present(verify)) lverify = verify
if (lverify .and. .not.shr_map_checkFldStr(fldStr)) then
call shr_map_abort(subName//' ERROR illegal fldStr '//trim(fldStr))
return
endif
if (trim(fldStr) == trim(shr_map_fs_nwts)) then
map%nwts = ival
elseif (trim(fldStr) == trim(shr_map_fs_nsrc)) then
map%nsrc = ival
elseif (trim(fldStr) == trim(shr_map_fs_ndst)) then
map%ndst = ival
else
call shr_map_abort(subName//' ERROR illegal fldStr '//trim(fldStr))
endif
end subroutine shr_map_putIN
!===============================================================================
!BOP ===========================================================================
!
! !IROUTINE: shr_map_putR8 -- put real to map
!
! !DESCRIPTION:
! one of the shr\_map\_put methods for reals
! puts value rval for input fldstr in map
! verify is optional argument that check validity and will
! call abort if rval is not valid option for fldstr.
!
! !REVISION HISTORY:
! 2005-Mar-27 - T. Craig - first version
!
! !INTERFACE: ------------------------------------------------------------------
subroutine shr_map_putR8(map,fldStr,rval,verify),2
implicit none
! !INPUT/OUTPUT PARAMETERS:
type(shr_map_mapType) ,intent(inout):: map
character(*) ,intent(in) :: fldStr
real(SHR_KIND_R8) ,intent(in) :: rval
logical,optional ,intent(in) :: verify ! check if string is valid
!EOP
!--- local ---
logical :: lverify
!--- formats ---
character(*),parameter :: subName = "('shr_map_putR8') "
!-------------------------------------------------------------------------------
lverify = .true.
if (present(verify)) lverify = verify
if (lverify .and. .not.shr_map_checkFldStr(fldStr)) then
call shr_map_abort(subName//' ERROR illegal fldStr '//trim(fldStr))
return
endif
call shr_map_abort(subName//' ERROR illegal fldStr '//trim(fldStr))
end subroutine shr_map_putR8
!===============================================================================
!BOP ===========================================================================
!
! !IROUTINE: shr_map_print -- write map to stdout
!
! !DESCRIPTION:
! Write map info to stdout
! \newline
! call shr\_map\_print(mymap)
!
! !REVISION HISTORY:
! 2005-Mar-27 - T. Craig - first version
!
! !INTERFACE: ------------------------------------------------------------------
subroutine shr_map_print(map) 3,3
implicit none
! !INPUT/OUTPUT PARAMETERS:
type(shr_map_mapType) ,intent(in) :: map
!EOP
!--- local ---
!--- formats ---
character(*),parameter :: subName = "('shr_map_print') "
character(*),parameter :: F00 = "('(shr_map_print) ',a) "
character(*),parameter :: F01 = "('(shr_map_print) ',a,2l2) "
character(*),parameter :: F02 = "('(shr_map_print) ',a,i8) "
character(*),parameter :: F03 = "('(shr_map_print) ',a,3i8) "
character(*),parameter :: F04 = "('(shr_map_print) ',a,2i8) "
character(*),parameter :: F05 = "('(shr_map_print) ',a,2e20.13) "
if (s_loglev > 0) then
write(s_logunit,*) ' '
write(s_logunit,F01) ' name : '//trim(map%name),shr_map_checkInit(map),shr_map_checkFilled(map)
write(s_logunit,F00) ' type : '//trim(map%type)
write(s_logunit,F00) ' algo : '//trim(map%algo)
write(s_logunit,F00) ' mask : '//trim(map%mask)
write(s_logunit,F00) ' vect : '//trim(map%vect)
write(s_logunit,F04) ' gsiz : ',map%nsrc,map%ndst
write(s_logunit,F05) ' xsrc : ',minval(map%xsrc),maxval(map%xsrc)
write(s_logunit,F05) ' ysrc : ',minval(map%ysrc),maxval(map%ysrc)
write(s_logunit,F05) ' xdst : ',minval(map%xdst),maxval(map%xdst)
write(s_logunit,F05) ' ydst : ',minval(map%ydst),maxval(map%ydst)
write(s_logunit,F02) ' nwts : ',map%nwts
write(s_logunit,F03) ' wsiz : ',size(map%wgts),size(map%isrc),size(map%idst)
write(s_logunit,F00) ' init : '//trim(map%init)
call shr_sys_flush(s_logunit)
endif
end subroutine shr_map_print
!===============================================================================
!BOP ===========================================================================
!
! !IROUTINE: shr_map_listValidOpts -- list the valid switches for map
!
! !DESCRIPTION:
! Lists the valid switches for map, informational only
! \newline
! call shr\_map\_listValidOpts()
!
! !REVISION HISTORY:
! 2005-Mar-27 - T. Craig - first version
!
! !INTERFACE: ------------------------------------------------------------------
subroutine shr_map_listValidOpts(),1
implicit none
! !INPUT/OUTPUT PARAMETERS:
!EOP
!--- local ---
integer(SHR_KIND_IN) :: n
!--- formats ---
character(*),parameter :: subName = "('shr_map_listValidOpts') "
character(*),parameter :: F00 = "('(shr_map_listValidOpts) ',a) "
!-------------------------------------------------------------------------------
if (s_loglev > 0) then
write(s_logunit,F00) ':'
write(s_logunit,F00) ' '//trim(shr_map_fs_name)//' : any character string'
do n = 1,shr_map_fs_ntype
write(s_logunit,F00) ' '//trim(shr_map_fs_type)//' : '//trim(shr_map_fs_types(n))
enddo
do n = 1,shr_map_fs_nalgo
write(s_logunit,F00) ' '//trim(shr_map_fs_algo)//' : '//trim(shr_map_fs_algos(n))
enddo
do n = 1,shr_map_fs_nmask
write(s_logunit,F00) ' '//trim(shr_map_fs_mask)//' : '//trim(shr_map_fs_masks(n))
enddo
do n = 1,shr_map_fs_nvect
write(s_logunit,F00) ' '//trim(shr_map_fs_vect)//' : '//trim(shr_map_fs_vects(n))
enddo
call shr_sys_flush(s_logunit)
endif
end subroutine shr_map_listValidOpts
!===============================================================================
!BOP ===========================================================================
!
! !IROUTINE: shr_map_clean -- cleans map
!
! !DESCRIPTION:
! Cleans map by resetting switches, deallocating arrays
! \newline
! call shr\_map\_clean(mymap)
!
! !REVISION HISTORY:
! 2005-Mar-27 - T. Craig - first version
!
! !INTERFACE: ------------------------------------------------------------------
subroutine shr_map_clean(map) 3
implicit none
! !INPUT/OUTPUT PARAMETERS:
type(shr_map_mapType) ,intent(inout):: map
!EOP
!--- local ---
integer :: rc
!--- formats ---
character(*),parameter :: subName = "('shr_map_clean') "
character(*),parameter :: F00 = "('(shr_map_clean) ',a) "
!-------------------------------------------------------------------------------
map%fill = ' '
map%init = ' '
call shr_map_put(map,shr_map_fs_name,shr_map_setFal,verify=.false.)
call shr_map_put(map,shr_map_fs_type,shr_map_setFal,verify=.false.)
call shr_map_put(map,shr_map_fs_algo,shr_map_setFal,verify=.false.)
call shr_map_put(map,shr_map_fs_mask,shr_map_setFal,verify=.false.)
call shr_map_put(map,shr_map_fs_mask,shr_map_setFal,verify=.false.)
call shr_map_put(map,shr_map_fs_nwts,shr_map_ispval)
call shr_map_put(map,shr_map_fs_nsrc,shr_map_ispval)
call shr_map_put(map,shr_map_fs_ndst,shr_map_ispval)
deallocate(map%xsrc,stat=rc)
if (rc > 0.and.debug > 1 .and. s_loglev > 0) write(s_logunit,F00) 'Warning: unable to deallocate map wgts'
deallocate(map%ysrc,stat=rc)
if (rc > 0.and.debug > 1 .and. s_loglev > 0) write(s_logunit,F00) 'Warning: unable to deallocate map wgts'
deallocate(map%xdst,stat=rc)
if (rc > 0.and.debug > 1 .and. s_loglev > 0) write(s_logunit,F00) 'Warning: unable to deallocate map wgts'
deallocate(map%ydst,stat=rc)
if (rc > 0.and.debug > 1 .and. s_loglev > 0) write(s_logunit,F00) 'Warning: unable to deallocate map wgts'
deallocate(map%wgts,stat=rc)
if (rc > 0.and.debug > 1 .and. s_loglev > 0) write(s_logunit,F00) 'Warning: unable to deallocate map wgts'
deallocate(map%isrc,stat=rc)
if (rc > 0.and.debug > 1 .and. s_loglev > 0) write(s_logunit,F00) 'Warning: unable to deallocate map isrc'
deallocate(map%idst,stat=rc)
if (rc > 0.and.debug > 1 .and. s_loglev > 0) write(s_logunit,F00) 'Warning: unable to deallocate map idst'
end subroutine shr_map_clean
!===============================================================================
!BOP ===========================================================================
!
! !IROUTINE: shr_map_mapSet_global -- Compute mapping weights
!
! !DESCRIPTION:
! Compute mapping weights based on setting in map. Fill the
! weights in the map. Currently supported maps and action:
! fill :copy = copy array by index, mask switch used
! fill :spval = copy array, fill with spval, mask switch not used
! fill :nn* = copy array, fill with nnval, mask switch not used
! remap:copy = copy array by index, mask switch used
! remap:spval = sets array to spval, mask switch used
! remap:bil* = bilinear interpolation, mask switch used
! remap:nn* = sets array to nnval, mask switch used
! \newline
! Requirements for input grids:
! Xsrc,Ysrc must be regular lat/lon grid, monotonically increasing,
! can be degrees or radians
! Xdst,Ydst are arbitrary list of lats/lons, must be same units as src
! Msrc,Mdst have nonzero for active grid point, zero for non-active
! src and dst must be the grid for type = fill
! Grids are check for validity
! \newline
! call shr\_map\_mapSet(mymap,Xs,Ys,Ms,Xd,Yd,Md)
! \newline
! call shr\_map\_mapSet(mymap,Xs,Ys,Ms,Xd,Yd,Md,algo='bilinear')
!
! !REMARKS
! If bothmask or srcmask is used with remap and some algorithms, active
! dst grid points can have invalid values. A report is produced after
! weights are calculated and this information will be detailed.
!
! !REVISION HISTORY:
! 2005-Mar-27 - T. Craig - first version
!
! !INTERFACE: ------------------------------------------------------------------
subroutine shr_map_mapSet_global(map,Xsrc,Ysrc,Msrc,Xdst_in,Ydst,Mdst,name,type,algo,mask,vect,rc),35
implicit none
! !INPUT/OUTPUT PARAMETERS:
type(shr_map_mapType) ,intent(inout):: map ! map
real(SHR_KIND_R8) ,intent(in) :: Xsrc(:,:) ! lon of src grid
real(SHR_KIND_R8) ,intent(in) :: Ysrc(:,:) ! lat of src grid
integer(SHR_KIND_IN) ,intent(in) :: Msrc(:,:) ! mask of src grid
real(SHR_KIND_R8) ,intent(in) :: Xdst_in(:,:) ! lon of dst grid
real(SHR_KIND_R8) ,intent(in) :: Ydst(:,:) ! lat of dst grid
integer(SHR_KIND_IN) ,intent(in) :: Mdst(:,:) ! mask of dst grid
character(*) ,optional,intent(in) :: name ! name
character(*) ,optional,intent(in) :: type ! type
character(*) ,optional,intent(in) :: algo ! algo
character(*) ,optional,intent(in) :: mask ! mask
character(*) ,optional,intent(in) :: vect ! vect
integer(SHR_KIND_IN),optional,intent(out) :: rc ! error code
!EOP
!--- local ---
integer(SHR_KIND_IN) :: nis,njs,nid,njd
integer(SHR_KIND_IN) :: nwts,n,n1,n2,ncnt,i,j,inn,jnn
integer(SHR_KIND_IN) :: irc,lrc
real(SHR_KIND_R8) :: rmin,rmax ! min/max value
real(SHR_KIND_R8) :: cang ! circle angle, deg or rad
real(SHR_KIND_R8),allocatable :: Xdst(:,:) ! lon of dst grid, wrapped as needed
integer(SHR_KIND_IN) :: pmax ! max num of wgts in pti...
integer(SHR_KIND_IN) :: ptot,ptot2 ! max num of wgts in lis...
integer(SHR_KIND_IN) :: pnum ! num of wgts set in pti...
integer(SHR_KIND_IN),allocatable :: pti(:) ! i index for wgts
integer(SHR_KIND_IN),allocatable :: ptj(:) ! j index for wgts
real(SHR_KIND_R8) ,allocatable :: ptw(:) ! weights for pti,ptj
integer(SHR_KIND_IN),allocatable :: lis(:) ! tmp src/dst index
integer(SHR_KIND_IN),allocatable :: lid(:) ! tmp src/dst index
real(SHR_KIND_R8) ,allocatable :: lwt(:) ! tmp wgt array
real(SHR_KIND_R8) ,allocatable :: sum(:) ! tmp sum array
integer(SHR_KIND_IN),allocatable :: ltmp(:) ! tmp src/dst index, for resize
real(SHR_KIND_R8) ,allocatable :: lwtmp(:) ! tmp wgt array, for resize
character(len=8) :: units ! radians or degrees
logical :: masksrc ! local var to turn on masking using src mask
logical :: maskdst ! local var to turn on masking using dst mask
logical :: maskdstbysrc ! local var to turn on masking using src mask for
! dst array, especially for fill
logical :: renorm ! local var to turn on renormalization
!--- formats ---
character(*),parameter :: subName = "('shr_map_mapSet_global') "
character(*),parameter :: F00 = "('(shr_map_mapSet_global) ',a) "
character(*),parameter :: F01 = "('(shr_map_mapSet_global) ',a,l2) "
character(*),parameter :: F02 = "('(shr_map_mapSet_global) ',a,2i8) "
character(*),parameter :: F03 = "('(shr_map_mapSet_global) ',a,2e20.13) "
!-------------------------------------------------------------------------------
lrc = 0
if (present(rc)) rc = lrc
if (present(name)) call shr_map_put(map,shr_map_fs_name,name)
if (present(type)) call shr_map_put(map,shr_map_fs_type,type,verify=.true.)
if (present(algo)) call shr_map_put(map,shr_map_fs_algo,algo,verify=.true.)
if (present(mask)) call shr_map_put(map,shr_map_fs_mask,mask,verify=.true.)
if (present(vect)) call shr_map_put(map,shr_map_fs_vect,vect,verify=.true.)
map%init = inispval
if (.NOT.shr_map_checkInit(map)) then
call shr_map_abort(subName//' ERROR map not initialized')
endif
!--- is lat/lon degrees or radians? ---
cang = 360._SHR_KIND_R8
units = 'degrees'
if (shr_map_checkRad(Ysrc)) then
cang=c2*pi
units = 'radians'
endif
nis = size(Xsrc,1)
njs = size(Xsrc,2)
nid = size(Xdst_in,1)
njd = size(Xdst_in,2)
!--- shift Xdst by 2pi to range of Xsrc as needed ---
allocate(Xdst(nid,njd))
rmin = minval(Xsrc)
rmax = maxval(Xsrc)
do j=1,njd
do i=1,nid
Xdst(i,j) = Xdst_in(i,j)
do while ((Xdst(i,j) < rmin .and. Xdst(i,j)+cang <= rmax).or. &
(Xdst(i,j) > rmax .and. Xdst(i,j)-cang >= rmin))
if (Xdst(i,j) < rmin) then
Xdst(i,j) = Xdst(i,j) + cang
elseif (Xdst(i,j) > rmax) then
Xdst(i,j) = Xdst(i,j) - cang
else
call shr_sys_abort(subName//' ERROR in Xdst wrap')
endif
enddo
enddo
enddo
call shr_map_checkGrids_global(Xsrc,Ysrc,Msrc,Xdst,Ydst,Mdst,map,lrc)
map%nwts = 0
map%nsrc = nis*njs
map%ndst = nid*njd
! deallocate(map%xsrc,stat=irc) ! this used to be a safe way to delloc when necessary,
! deallocate(map%ysrc,stat=irc) ! but do nothing when pointers were undefined or
! deallocate(map%xdst,stat=irc) ! un-associated, in Oct 2005, undefined ptrs started
! deallocate(map%ydst,stat=irc) ! causing seg-faults on bluesky (B. Kauffman)
allocate(map%xsrc(nis*njs))
allocate(map%ysrc(nis*njs))
allocate(map%xdst(nid*njd))
allocate(map%ydst(nid*njd))
do j=1,njs
do i=1,nis
call shr_map_2dto1d(n1,nis,njs,i,j)
map%xsrc(n1) = Xsrc(i,j)*c2*pi/cang
map%ysrc(n1) = Ysrc(i,j)*c2*pi/cang
enddo
enddo
do j=1,njd
do i=1,nid
call shr_map_2dto1d(n1,nid,njd,i,j)
map%xdst(n1) = Xdst(i,j)*c2*pi/cang
map%ydst(n1) = Ydst(i,j)*c2*pi/cang
enddo
enddo
masksrc = .false.
maskdstbysrc = .false.
maskdst = .false.
renorm = .true.
if (trim(map%type) /= trim(shr_map_fs_fill) .and. &
trim(map%type) /= trim(shr_map_fs_cfill)) then
if (trim(map%mask) == trim(shr_map_fs_bothmask) .or. &
trim(map%mask) == trim(shr_map_fs_srcmask)) masksrc = .true.
if (trim(map%mask) == trim(shr_map_fs_bothmask) .or. &
trim(map%mask) == trim(shr_map_fs_dstmask)) maskdst = .true.
endif
if (trim(map%algo) == trim(shr_map_fs_spval)) then
masksrc = .false.
renorm = .false.
endif
if (debug > 1) then
if (s_loglev > 0) write(s_logunit,*) ' '
call shr_map_print(map)
endif
if (lrc /= 0) then
if (present(rc)) rc = lrc
return
endif
if (trim(map%algo) == trim(shr_map_fs_bilinear)) then
if (dopole) then
pmax = nis+2 ! possible for high lat points
ptot = 4*nid*njd ! start with bilinear estimate
else
pmax = 4 ! bilinear with 4 wts/map
ptot = 4*nid*njd
endif
else
pmax = 1 ! nn with 1 wts/map
ptot = 1*nid*njd
endif
allocate(lis(ptot))
allocate(lid(ptot))
allocate(lwt(ptot))
allocate(pti(pmax))
allocate(ptj(pmax))
allocate(ptw(pmax))
!--- full array copy is default ---
nwts = nid*njd
do n=1,nwts
lid(n) = n
lis(n) = mod(n-1,nis*njs)+1
lwt(n) = c1
enddo
!--- index copy anytime algo = copy ---
if (trim(map%algo) == trim(shr_map_fs_copy)) then
map%init = initcopy
! just use copy default
!--- for fill ---
elseif (trim(map%type) == trim(shr_map_fs_fill) .or. &
trim(map%type) == trim(shr_map_fs_cfill)) then
map%init = initcopy
if (trim(map%algo) == trim(shr_map_fs_spval)) then
maskdstbysrc = .true.
elseif (trim(map%algo) == trim(shr_map_fs_nn)) then
do n=1,nwts
call shr_map_1dto2d(lis(n),nis,njs,i,j)
if (Msrc(i,j) == 0) then
call shr_map_findnn(Xsrc(i,j),Ysrc(i,j),Xsrc,Ysrc,Msrc,inn,jnn)
call shr_map_2dto1d(lis(n),nis,njs,inn,jnn)
endif
enddo
elseif (trim(map%algo) == trim(shr_map_fs_nnoni)) then
do n=1,nwts
call shr_map_1dto2d(lis(n),nis,njs,i,j)
if (Msrc(i,j) == 0) then
call shr_map_findnnon('i',Xsrc(i,j),Ysrc(i,j),Xsrc,Ysrc,Msrc,inn,jnn)
call shr_map_2dto1d(lis(n),nis,njs,inn,jnn)
endif
enddo
elseif (trim(map%algo) == trim(shr_map_fs_nnonj)) then
do n=1,nwts
call shr_map_1dto2d(lis(n),nis,njs,i,j)
if (Msrc(i,j) == 0) then
call shr_map_findnnon('j',Xsrc(i,j),Ysrc(i,j),Xsrc,Ysrc,Msrc,inn,jnn)
call shr_map_2dto1d(lis(n),nis,njs,inn,jnn)
endif
enddo
else
call shr_map_abort(subName//' ERROR: unsupported map option combo')
endif
!--- for remap ---
elseif (trim(map%type) == trim(shr_map_fs_remap)) then
map%init = inispval
if (trim(map%algo) == trim(shr_map_fs_spval)) then
nwts = 0
elseif (trim(map%algo) == trim(shr_map_fs_nn)) then
do n=1,nwts
call shr_map_1dto2d(lid(n),nid,njd,i,j)
call shr_map_findnn(Xdst(i,j),Ydst(i,j),Xsrc,Ysrc,Msrc,inn,jnn)
call shr_map_2dto1d(lis(n),nis,njs,inn,jnn)
enddo
elseif (trim(map%algo) == trim(shr_map_fs_nnoni)) then
do n=1,nwts
call shr_map_1dto2d(lid(n),nid,njd,i,j)
call shr_map_findnnon('i',Xdst(i,j),Ydst(i,j),Xsrc,Ysrc,Msrc,inn,jnn)
call shr_map_2dto1d(lis(n),nis,njs,inn,jnn)
enddo
elseif (trim(map%algo) == trim(shr_map_fs_nnonj)) then
do n=1,nwts
call shr_map_1dto2d(lid(n),nid,njd,i,j)
call shr_map_findnnon('j',Xdst(i,j),Ydst(i,j),Xsrc,Ysrc,Msrc,inn,jnn)
call shr_map_2dto1d(lis(n),nis,njs,inn,jnn)
enddo
elseif (trim(map%algo) == trim(shr_map_fs_bilinear)) then
nwts = 0
do n=1,nid*njd
call shr_map_1dto2d(n,nid,njd,i,j)
call shr_map_getWts(Xdst(i,j),Ydst(i,j),Xsrc,Ysrc,pti,ptj,ptw,pnum,units)
if (nwts + pnum > size(lwt)) then
!--- resize lis, lid, lwt. ptot is old size, ptot2 is new size
ptot = size(lwt)
ptot2 = ptot + max(ptot/2,pnum*10)
if (debug > 1 .and. s_loglev > 0) write(s_logunit,F02) 'resize wts ',ptot,ptot2
allocate(ltmp(ptot))
ltmp(1:nwts) = lis(1:nwts)
deallocate(lis)
allocate(lis(ptot2))
lis(1:nwts) = ltmp(1:nwts)
ltmp(1:nwts) = lid(1:nwts)
deallocate(lid)
allocate(lid(ptot2))
lid(1:nwts) = ltmp(1:nwts)
deallocate(ltmp)
allocate(lwtmp(ptot))
lwtmp(1:nwts) = lwt(1:nwts)
deallocate(lwt)
allocate(lwt(ptot2))
lwt(1:nwts) = lwtmp(1:nwts)
deallocate(lwtmp)
endif
do n1 = 1,pnum
nwts = nwts + 1
lid(nwts) = n
call shr_map_2dto1d(lis(nwts),nis,njs,pti(n1),ptj(n1))
lwt(nwts) = ptw(n1)
enddo
enddo
else
call shr_map_abort(subName//' ERROR: unsupported map option combo')
if (present(rc)) rc = 1
return
endif
else
call shr_map_abort(subName//' ERROR: unsupported map option combo')
if (present(rc)) rc = 1
return
endif
!--- compress weights and copy to map ---
!--- remove 1:1 copies if initcopy
if (debug > 1 .and. s_loglev > 0) write(s_logunit,F01) 'init: ',map%init
if (map%init == initcopy .and. &
trim(map%type) /= trim(shr_map_fs_cfill)) then
ncnt = 0
do n=1,nwts
if (lid(n) == lis(n) .and. abs(lwt(n)-c1) < eps) then
! skipit
else
ncnt = ncnt+1
lid(ncnt) = lid(n)
lis(ncnt) = lis(n)
lwt(ncnt) = lwt(n)
endif
enddo
if (debug > 1 .and. s_loglev > 0) write(s_logunit,F02) ' rm dst grid points, nwts old/new = ',nwts,ncnt
nwts = ncnt
endif
!--- remove dst grid points ---
if (debug > 1 .and. s_loglev > 0) write(s_logunit,F01) 'maskdst: ',maskdst
if (maskdst) then
ncnt = 0
do n=1,nwts
call shr_map_1dto2d(lid(n),nid,njd,i,j)
if (Mdst(i,j) /= 0) then
ncnt = ncnt+1
lid(ncnt) = lid(n)
lis(ncnt) = lis(n)
lwt(ncnt) = lwt(n)
endif
enddo
if (debug > 1 .and. s_loglev > 0) write(s_logunit,F02) ' rm dst grid points, nwts old/new = ',nwts,ncnt
nwts = ncnt
endif
!--- remove dst grid points based on src mask---
if (debug > 1 .and. s_loglev > 0) write(s_logunit,F01) 'maskdstbysrc: ',maskdstbysrc
if (maskdstbysrc) then
ncnt = 0
do n=1,nwts
call shr_map_1dto2d(lid(n),nid,njd,i,j)
if (Msrc(i,j) /= 0) then
ncnt = ncnt+1
lid(ncnt) = lid(n)
lis(ncnt) = lis(n)
lwt(ncnt) = lwt(n)
endif
enddo
if (debug > 1 .and. s_loglev > 0) write(s_logunit,F02) ' rm dst grid points by src, nwts old/new = ',nwts,ncnt
nwts = ncnt
endif
!--- remove src grid points ---
if (debug > 1 .and. s_loglev > 0) write(s_logunit,F01) 'masksrc: ',masksrc
if (masksrc) then
ncnt = 0
do n=1,nwts
call shr_map_1dto2d(lis(n),nis,njs,i,j)
if (Msrc(i,j) /= 0) then
ncnt = ncnt+1
lid(ncnt) = lid(n)
lis(ncnt) = lis(n)
lwt(ncnt) = lwt(n)
endif
enddo
if (debug > 1 .and. s_loglev > 0) write(s_logunit,F02) ' rm src grid points, nwts old/new = ',nwts,ncnt
nwts = ncnt
endif
!--- renormalize wgts to 1.0 ---
allocate(sum(nid*njd))
!--- sum weights for dst grid points ---
sum(:) = c0
do n=1,nwts
sum(lid(n)) = sum(lid(n)) + lwt(n)
enddo
!--- print min/max sum ---
rmin = maxval(sum)
rmax = minval(sum)
do n=1,nid*njd
if (sum(n) > eps) then
rmin = min(rmin,sum(n))
rmax = max(rmax,sum(n))
endif
enddo
if (debug > 1 .and. s_loglev > 0) write(s_logunit,F03) 'sum wts min/max ',rmin,rmax
!--- renormalize so sum on destination is always 1.0 for active dst points
if (debug > 1 .and. s_loglev > 0) write(s_logunit,F01) 'renorm: ',renorm
if (renorm) then
do n=1,nwts
if (sum(lid(n)) > eps) then
lwt(n) = lwt(n) / sum(lid(n))
endif
enddo
!--- sum weights for dst grid points ---
sum(:) = c0
do n=1,nwts
sum(lid(n)) = sum(lid(n)) + lwt(n)
enddo
!--- print min/max sum ---
rmin = maxval(sum)
rmax = minval(sum)
do n=1,nid*njd
if (sum(n) > eps) then
rmin = min(rmin,sum(n))
rmax = max(rmax,sum(n))
endif
enddo
if (debug > 1 .and. s_loglev > 0) write(s_logunit,F03) 'sum wts min/max ',rmin,rmax
endif
map%nwts = nwts
! deallocate(map%idst,stat=irc)
! deallocate(map%isrc,stat=irc)
! deallocate(map%wgts,stat=irc)
allocate(map%idst(nwts))
allocate(map%isrc(nwts))
allocate(map%wgts(nwts))
do n=1,nwts
map%idst(n) = lid(n)
map%isrc(n) = lis(n)
map%wgts(n) = lwt(n)
enddo
deallocate(Xdst)
deallocate(lis)
deallocate(lid)
deallocate(lwt)
deallocate(sum)
deallocate(pti)
deallocate(ptj)
deallocate(ptw)
map%fill = fillstring
call shr_map_checkWgts_global(Msrc,Mdst,map)
if (present(rc)) rc = lrc
end subroutine shr_map_mapSet_global
!===============================================================================
!BOP ===========================================================================
!
! !IROUTINE: shr_map_mapSet_dest -- Compute mapping weights
!
! !DESCRIPTION:
! Compute mapping weights based on setting in map. Fill the
! weights in the map. Currently supported maps and action:
! fill :copy = copy array by index, mask switch used
! fill :spval = copy array, fill with spval, mask switch not used
! fill :nn* = copy array, fill with nnval, mask switch not used
! remap:copy = copy array by index, mask switch used
! remap:spval = sets array to spval, mask switch used
! remap:bil* = bilinear interpolation, mask switch used
! remap:nn* = sets array to nnval, mask switch used
! \newline
! Requirements for input grids:
! Xsrc,Ysrc must be regular lat/lon grid, monotonically increasing,
! can be degrees or radians
! Xdst,Ydst are arbitrary list of lats/lons, must be same units as src
! Msrc,Mdst have nonzero for active grid point, zero for non-active
! src and dst must be the grid for type = fill
! Grids are check for validity
! \newline
! call shr\_map\_mapSet(mymap,Xs,Ys,Ms,Xd,Yd,Md)
! \newline
! call shr\_map\_mapSet(mymap,Xs,Ys,Ms,Xd,Yd,Md,algo='bilinear')
!
! !REMARKS
! If bothmask or srcmask is used with remap and some algorithms, active
! dst grid points can have invalid values. A report is produced after
! weights are calculated and this information will be detailed.
!
! !REVISION HISTORY:
! 2005-Mar-27 - T. Craig - first version
!
! !INTERFACE: ------------------------------------------------------------------
subroutine shr_map_mapSet_dest(map,Xsrc,Ysrc,Msrc,Xdst_in,Ydst,Mdst,ndst,Idst,name,type,algo,mask,vect,rc),26
implicit none
! !INPUT/OUTPUT PARAMETERS:
type(shr_map_mapType) ,intent(inout):: map ! map
real(SHR_KIND_R8) ,intent(in) :: Xsrc(:,:) ! lon of src grid
real(SHR_KIND_R8) ,intent(in) :: Ysrc(:,:) ! lat of src grid
integer(SHR_KIND_IN) ,intent(in) :: Msrc(:,:) ! mask of src grid
real(SHR_KIND_R8) ,intent(in) :: Xdst_in(:) ! lon of dst grid
real(SHR_KIND_R8) ,intent(in) :: Ydst(:) ! lat of dst grid
integer(SHR_KIND_IN) ,intent(in) :: Mdst(:) ! mask of dst grid
integer(SHR_KIND_IN) ,intent(in) :: ndst ! global size of dst
integer(SHR_KIND_IN) ,intent(in) :: Idst(:) ! global index of dst grid
character(*) ,optional,intent(in) :: name ! name
character(*) ,optional,intent(in) :: type ! type
character(*) ,optional,intent(in) :: algo ! algo
character(*) ,optional,intent(in) :: mask ! mask
character(*) ,optional,intent(in) :: vect ! vect
integer(SHR_KIND_IN),optional,intent(out) :: rc ! error code
!EOP
!--- local ---
integer(SHR_KIND_IN) :: nis,njs,nid,njd
integer(SHR_KIND_IN) :: nwts,n,n1,n2,ncnt,i,j,inn,jnn
integer(SHR_KIND_IN) :: irc,lrc
real(SHR_KIND_R8) :: rmin,rmax ! min/max value
real(SHR_KIND_R8) :: cang ! circle angle, deg or rad
real(SHR_KIND_R8),allocatable :: Xdst(:) ! lon of dst grid, wrapped as needed
integer(SHR_KIND_IN) :: pmax ! max num of wgts in pti...
integer(SHR_KIND_IN) :: ptot,ptot2 ! max num of wgts in lis...
integer(SHR_KIND_IN) :: pnum ! num of wgts set in pti...
integer(SHR_KIND_IN),allocatable :: pti(:) ! i index for wgts
integer(SHR_KIND_IN),allocatable :: ptj(:) ! j index for wgts
real(SHR_KIND_R8) ,allocatable :: ptw(:) ! weights for pti,ptj
integer(SHR_KIND_IN),allocatable :: lis(:) ! tmp src/dst index
integer(SHR_KIND_IN),allocatable :: lid(:) ! tmp src/dst index
real(SHR_KIND_R8) ,allocatable :: lwt(:) ! tmp wgt array
real(SHR_KIND_R8) ,allocatable :: sum(:) ! tmp sum array
integer(SHR_KIND_IN),allocatable :: ltmp(:) ! tmp src/dst index, for resize
real(SHR_KIND_R8) ,allocatable :: lwtmp(:) ! tmp wgt array, for resize
character(len=8) :: units ! radians or degrees
logical :: masksrc ! local var to turn on masking using src mask
logical :: maskdst ! local var to turn on masking using dst mask
logical :: maskdstbysrc ! local var to turn on masking using src mask for
! dst array, especially for fill
logical :: renorm ! local var to turn on renormalization
!--- formats ---
character(*),parameter :: subName = "('shr_map_mapSet_dest') "
character(*),parameter :: F00 = "('(shr_map_mapSet_dest) ',a) "
character(*),parameter :: F01 = "('(shr_map_mapSet_dest) ',a,l2) "
character(*),parameter :: F02 = "('(shr_map_mapSet_dest) ',a,2i8) "
character(*),parameter :: F03 = "('(shr_map_mapSet_dest) ',a,2e20.13) "
!-------------------------------------------------------------------------------
write(s_logunit,F00) 'ERROR this routine is not validated'
call shr_sys_abort(subName//' ERROR subroutine not validated')
lrc = 0
if (present(rc)) rc = lrc
if (present(name)) call shr_map_put(map,shr_map_fs_name,name)
if (present(type)) call shr_map_put(map,shr_map_fs_type,type,verify=.true.)
if (present(algo)) call shr_map_put(map,shr_map_fs_algo,algo,verify=.true.)
if (present(mask)) call shr_map_put(map,shr_map_fs_mask,mask,verify=.true.)
if (present(vect)) call shr_map_put(map,shr_map_fs_vect,vect,verify=.true.)
map%init = inispval
if (.NOT.shr_map_checkInit(map)) then
call shr_map_abort(subName//' ERROR map not initialized')
endif
!--- is lat/lon degrees or radians? ---
cang = 360._SHR_KIND_R8
units = 'degrees'
if (shr_map_checkRad(Ysrc)) then
cang=c2*pi
units = 'radians'
endif
nis = size(Xsrc,1)
njs = size(Xsrc,2)
nid = size(Xdst_in,1)
njd = 1
!--- shift Xdst by 2pi to range of Xsrc as needed ---
allocate(Xdst(nid))
rmin = minval(Xsrc)
rmax = maxval(Xsrc)
do i=1,nid
Xdst(i) = Xdst_in(i)
do while ((Xdst(i) < rmin .and. Xdst(i)+cang <= rmax).or. &
(Xdst(i) > rmax .and. Xdst(i)-cang >= rmin))
if (Xdst(i) < rmin) then
Xdst(i) = Xdst(i) + cang
elseif (Xdst(i) > rmax) then
Xdst(i) = Xdst(i) - cang
else
call shr_sys_abort(subName//' ERROR in Xdst wrap')
endif
enddo
enddo
call shr_map_checkGrids_dest(Xsrc,Ysrc,Msrc,Xdst,Ydst,Mdst,map,lrc)
map%nwts = 0
map%nsrc = nis*njs
map%ndst = ndst
! deallocate(map%xsrc,stat=irc) ! this used to be a safe way to delloc when necessary,
! deallocate(map%ysrc,stat=irc) ! but do nothing when pointers were undefined or
! deallocate(map%xdst,stat=irc) ! un-associated, in Oct 2005, undefined ptrs started
! deallocate(map%ydst,stat=irc) ! causing seg-faults on bluesky (B. Kauffman)
allocate(map%xsrc(nis*njs))
allocate(map%ysrc(nis*njs))
allocate(map%xdst(nid*njd))
allocate(map%ydst(nid*njd))
do j=1,njs
do i=1,nis
call shr_map_2dto1d(n1,nis,njs,i,j)
map%xsrc(n1) = Xsrc(i,j)*c2*pi/cang
map%ysrc(n1) = Ysrc(i,j)*c2*pi/cang
enddo
enddo
do i=1,nid
map%xdst(i) = Xdst(i)*c2*pi/cang
map%ydst(i) = Ydst(i)*c2*pi/cang
enddo
masksrc = .false.
maskdstbysrc = .false.
maskdst = .false.
renorm = .true.
if (trim(map%type) /= trim(shr_map_fs_fill) .and. &
trim(map%type) /= trim(shr_map_fs_cfill)) then
if (trim(map%mask) == trim(shr_map_fs_bothmask) .or. &
trim(map%mask) == trim(shr_map_fs_srcmask)) masksrc = .true.
if (trim(map%mask) == trim(shr_map_fs_bothmask) .or. &
trim(map%mask) == trim(shr_map_fs_dstmask)) maskdst = .true.
endif
if (trim(map%algo) == trim(shr_map_fs_spval)) then
masksrc = .false.
renorm = .false.
endif
if (debug > 1) then
if (s_loglev > 0) write(s_logunit,*) ' '
call shr_map_print(map)
endif
if (lrc /= 0) then
if (present(rc)) rc = lrc
return
endif
if (trim(map%algo) == trim(shr_map_fs_bilinear)) then
if (dopole) then
pmax = nis+2 ! possible for high lat points
ptot = 4*nid*njd ! start with bilinear estimate
else
pmax = 4 ! bilinear with 4 wts/map
ptot = 4*nid*njd
endif
else
pmax = 1 ! nn with 1 wts/map
ptot = 1*nid*njd
endif
allocate(lis(ptot))
allocate(lid(ptot))
allocate(lwt(ptot))
allocate(pti(pmax))
allocate(ptj(pmax))
allocate(ptw(pmax))
!--- full array copy is default ---
nwts = nid*njd
do n=1,nwts
lid(n) = Idst(n)
lis(n) = Idst(n)
lwt(n) = c1
enddo
!--- index copy anytime algo = copy ---
if (trim(map%algo) == trim(shr_map_fs_copy)) then
map%init = initcopy
! just use copy default
!--- for fill ---
elseif (trim(map%type) == trim(shr_map_fs_fill) .or. &
trim(map%type) == trim(shr_map_fs_cfill)) then
map%init = initcopy
if (trim(map%algo) == trim(shr_map_fs_spval)) then
maskdstbysrc = .true.
elseif (trim(map%algo) == trim(shr_map_fs_nn)) then
do n=1,nwts
if (Mdst(n) == 0) then
call shr_map_findnn(Xdst(n),Ydst(n),Xsrc,Ysrc,Msrc,inn,jnn)
call shr_map_2dto1d(lis(n),nis,njs,inn,jnn)
endif
enddo
elseif (trim(map%algo) == trim(shr_map_fs_nnoni)) then
do n=1,nwts
if (Mdst(n) == 0) then
call shr_map_findnnon('i',Xdst(n),Ydst(n),Xsrc,Ysrc,Msrc,inn,jnn)
call shr_map_2dto1d(lis(n),nis,njs,inn,jnn)
endif
enddo
elseif (trim(map%algo) == trim(shr_map_fs_nnonj)) then
do n=1,nwts
if (Mdst(n) == 0) then
call shr_map_findnnon('j',Xdst(n),Ydst(n),Xsrc,Ysrc,Msrc,inn,jnn)
call shr_map_2dto1d(lis(n),nis,njs,inn,jnn)
endif
enddo
else
call shr_map_abort(subName//' ERROR: unsupported map option combo')
endif
!--- for remap ---
elseif (trim(map%type) == trim(shr_map_fs_remap)) then
map%init = inispval
if (trim(map%algo) == trim(shr_map_fs_spval)) then
nwts = 0
elseif (trim(map%algo) == trim(shr_map_fs_nn)) then
do n=1,nwts
call shr_map_findnn(Xdst(n),Ydst(n),Xsrc,Ysrc,Msrc,inn,jnn)
call shr_map_2dto1d(lis(n),nis,njs,inn,jnn)
enddo
elseif (trim(map%algo) == trim(shr_map_fs_nnoni)) then
do n=1,nwts
call shr_map_findnnon('i',Xdst(n),Ydst(n),Xsrc,Ysrc,Msrc,inn,jnn)
call shr_map_2dto1d(lis(n),nis,njs,inn,jnn)
enddo
elseif (trim(map%algo) == trim(shr_map_fs_nnonj)) then
do n=1,nwts
call shr_map_findnnon('j',Xdst(n),Ydst(n),Xsrc,Ysrc,Msrc,inn,jnn)
call shr_map_2dto1d(lis(n),nis,njs,inn,jnn)
enddo
elseif (trim(map%algo) == trim(shr_map_fs_bilinear)) then
nwts = 0
do n=1,nid*njd
call shr_map_getWts(Xdst(n),Ydst(n),Xsrc,Ysrc,pti,ptj,ptw,pnum,units)
if (nwts + pnum > size(lwt)) then
!--- resize lis, lid, lwt. ptot is old size, ptot2 is new size
ptot = size(lwt)
ptot2 = ptot + max(ptot/2,pnum*10)
if (debug > 1 .and. s_loglev > 0) write(s_logunit,F02) 'resize wts ',ptot,ptot2
allocate(ltmp(ptot))
ltmp(1:nwts) = lis(1:nwts)
deallocate(lis)
allocate(lis(ptot2))
lis(1:nwts) = ltmp(1:nwts)
ltmp(1:nwts) = lid(1:nwts)
deallocate(lid)
allocate(lid(ptot2))
lid(1:nwts) = ltmp(1:nwts)
deallocate(ltmp)
allocate(lwtmp(ptot))
lwtmp(1:nwts) = lwt(1:nwts)
deallocate(lwt)
allocate(lwt(ptot2))
lwt(1:nwts) = lwtmp(1:nwts)
deallocate(lwtmp)
endif
do n1 = 1,pnum
nwts = nwts + 1
lid(nwts) = Idst(n)
call shr_map_2dto1d(lis(nwts),nis,njs,pti(n1),ptj(n1))
lwt(nwts) = ptw(n1)
enddo
enddo
else
call shr_map_abort(subName//' ERROR: unsupported map option combo')
if (present(rc)) rc = 1
return
endif
else
call shr_map_abort(subName//' ERROR: unsupported map option combo')
if (present(rc)) rc = 1
return
endif
!--- compress weights and copy to map ---
!--- remove 1:1 copies if initcopy
if (debug > 1 .and. s_loglev > 0) write(s_logunit,F01) 'init: ',map%init
if (map%init == initcopy .and. &
trim(map%type) /= trim(shr_map_fs_cfill)) then
ncnt = 0
do n=1,nwts
if (lid(n) == lis(n) .and. abs(lwt(n)-c1) < eps) then
! skipit
else
ncnt = ncnt+1
lid(ncnt) = lid(n)
lis(ncnt) = lis(n)
lwt(ncnt) = lwt(n)
endif
enddo
if (debug > 1 .and. s_loglev > 0) write(s_logunit,F02) ' rm dst grid points, nwts old/new = ',nwts,ncnt
nwts = ncnt
endif
!--- remove dst grid points ---
if (debug > 1 .and. s_loglev > 0) write(s_logunit,F01) 'maskdst: ',maskdst
if (maskdst) then
ncnt = 0
do n=1,nwts
if (Mdst(n) /= 0) then
ncnt = ncnt+1
lid(ncnt) = lid(n)
lis(ncnt) = lis(n)
lwt(ncnt) = lwt(n)
endif
enddo
if (debug > 1 .and. s_loglev > 0) write(s_logunit,F02) ' rm dst grid points, nwts old/new = ',nwts,ncnt
nwts = ncnt
endif
!--- remove dst grid points based on src mask---
if (debug > 1 .and. s_loglev > 0) write(s_logunit,F01) 'maskdstbysrc: ',maskdstbysrc
if (maskdstbysrc) then
ncnt = 0
do n=1,nwts
call shr_map_1dto2d(lid(n),nis,njs,i,j)
if (Msrc(i,j) /= 0) then
ncnt = ncnt+1
lid(ncnt) = lid(n)
lis(ncnt) = lis(n)
lwt(ncnt) = lwt(n)
endif
enddo
if (debug > 1 .and. s_loglev > 0) write(s_logunit,F02) ' rm dst grid points by src, nwts old/new = ',nwts,ncnt
nwts = ncnt
endif
!--- remove src grid points ---
if (debug > 1 .and. s_loglev > 0) write(s_logunit,F01) 'masksrc: ',masksrc
if (masksrc) then
ncnt = 0
do n=1,nwts
call shr_map_1dto2d(lis(n),nis,njs,i,j)
if (Msrc(i,j) /= 0) then
ncnt = ncnt+1
lid(ncnt) = lid(n)
lis(ncnt) = lis(n)
lwt(ncnt) = lwt(n)
endif
enddo
if (debug > 1 .and. s_loglev > 0) write(s_logunit,F02) ' rm src grid points, nwts old/new = ',nwts,ncnt
nwts = ncnt
endif
!--- renormalize wgts to 1.0 ---
allocate(sum(ndst))
!--- sum weights for dst grid points ---
sum(:) = c0
do n=1,nwts
sum(lid(n)) = sum(lid(n)) + lwt(n)
enddo
!--- print min/max sum ---
rmin = maxval(sum)
rmax = minval(sum)
do n=1,ndst
if (sum(n) > eps) then
rmin = min(rmin,sum(n))
rmax = max(rmax,sum(n))
endif
enddo
if (debug > 1 .and. s_loglev > 0) write(s_logunit,F03) 'sum wts min/max ',rmin,rmax
!--- renormalize so sum on destination is always 1.0 for active dst points
if (debug > 1 .and. s_loglev > 0) write(s_logunit,F01) 'renorm: ',renorm
if (renorm) then
do n=1,nwts
if (sum(lid(n)) > eps) then
lwt(n) = lwt(n) / sum(lid(n))
endif
enddo
!--- sum weights for dst grid points ---
sum(:) = c0
do n=1,nwts
sum(lid(n)) = sum(lid(n)) + lwt(n)
enddo
!--- print min/max sum ---
rmin = maxval(sum)
rmax = minval(sum)
do n=1,nid*njd
if (sum(n) > eps) then
rmin = min(rmin,sum(n))
rmax = max(rmax,sum(n))
endif
enddo
if (debug > 1 .and. s_loglev > 0) write(s_logunit,F03) 'sum wts min/max ',rmin,rmax
endif
map%nwts = nwts
! deallocate(map%idst,stat=irc)
! deallocate(map%isrc,stat=irc)
! deallocate(map%wgts,stat=irc)
allocate(map%idst(nwts))
allocate(map%isrc(nwts))
allocate(map%wgts(nwts))
do n=1,nwts
map%idst(n) = lid(n)
map%isrc(n) = lis(n)
map%wgts(n) = lwt(n)
enddo
deallocate(Xdst)
deallocate(lis)
deallocate(lid)
deallocate(lwt)
deallocate(sum)
deallocate(pti)
deallocate(ptj)
deallocate(ptw)
map%fill = fillstring
!! call shr_map_checkWgts_dest(Msrc,Mdst,map)
if (present(rc)) rc = lrc
end subroutine shr_map_mapSet_dest
!===============================================================================
!BOP ===========================================================================
!
! !IROUTINE: shr_map_mapDatam -- maps arrays using input map
!
! !DESCRIPTION:
! Maps arrays using preset map
! \newline
! call shr\_map\_mapData(Ain,Aout,mymap)
!
! !REVISION HISTORY:
! 2005-Mar-27 - T. Craig - first version
!
! !INTERFACE: ------------------------------------------------------------------
subroutine shr_map_mapDatam(arrsrc,arrdst,map),22
!--- map arrsrc to arrdst, each array is dimension (fields,grid index) ---
implicit none
! !INPUT/OUTPUT PARAMETERS:
real(SHR_KIND_R8) ,intent(in) :: arrsrc(:,:) ! src array(fields,grid)
real(SHR_KIND_R8) ,intent(out):: arrdst(:,:) ! dst array(fields,grid)
type(shr_map_mapType) ,intent(in) :: map ! map
!EOP
!--- local ---
integer(SHR_KIND_IN) :: n,n2 ! counters
integer(SHR_KIND_IN) :: indi,indo ! array indices, in/out
real(SHR_KIND_R8) :: wgt ! value of weight
integer(SHR_KIND_IN) :: nfi,nfo ! number of fields in array, in/out
integer(SHR_KIND_IN) :: nsi,nso ! size of grid in array, in/out
real(SHR_KIND_R8) :: theta ! angle difference
integer(SHR_KIND_IN),save :: t0=-1,t1,t2,t3,t4,t5 ! timers
integer(SHR_KIND_IN),parameter :: timing=0 ! turn timers off/on (0/1)
logical,pointer :: initnew(:) ! mask for initialization
!--- formats ---
character(*),parameter :: subName = "('shr_map_mapDatam') "
character(*),parameter :: F00 = "('(shr_map_mapDatam) ',a) "
character(*),parameter :: F01 = "('(shr_map_mapDatam) ',a,2i8) "
!-------------------------------------------------------------------------------
if (timing>0 .and. t0 == -1) then
call shr_timer_get(t0,subName//"everything")
call shr_timer_get(t1,subName//"initial checks")
call shr_timer_get(t2,subName//"dst to spval")
call shr_timer_get(t4,subName//"map vector")
call shr_timer_get(t5,subName//"map scalar")
end if
if (timing>0) call shr_timer_start(t0)
if (timing>0) call shr_timer_start(t1)
!--- get number of fields ---
nfi = size(arrsrc,1)
nfo = size(arrdst,1)
!--- check number of fields ---
if (nfi /= nfo) then
write(s_logunit,F01) ' field numbers dont match ',nfi,nfo
call shr_map_abort(subName//' ERROR number of fields')
endif
!--- check two fields for vector ---
if (trim(map%vect) == trim(shr_map_fs_vector).and.(nfi /= 2)) then
write(s_logunit,F01) ' vector mapping, must map only two fields',nfi,nfo
call shr_map_abort(subName//' ERROR vector mapping fields not two')
endif
!--- check that map is set ---
if (.not.shr_map_checkFilled(map)) then
write(s_logunit,F00) ' map is not filled'
call shr_map_abort(subName//' ERROR map is not filled')
endif
!--- get size of grid ---
nsi = size(arrsrc,2)
nso = size(arrdst,2)
!--- check size of grid ---
if (nsi /= map%nsrc) then
write(s_logunit,F01) ' src grid size doesnt match ',nsi,map%nsrc
call shr_map_abort(subName//' ERROR src grid size')
endif
if (nso /= map%ndst) then
write(s_logunit,F01) ' dst grid size doesnt match ',nso,map%ndst
call shr_map_abort(subName//' ERROR dst grid size')
endif
if (timing>0) call shr_timer_stop(t1)
if (timing>0) call shr_timer_start(t2)
allocate(initnew(1:nso))
initnew = .true.
!--- set arrdst to spval, all points, default ---
if (map%init == inispval) then
arrdst = shr_map_spval
elseif (map%init == initcopy) then
if (nsi /= nso) then
write(s_logunit,F01) ' initcopy has nsi ne nso ',nsi,nso
call shr_map_abort(subName//' ERROR initcopy size')
else
do n = 1,nsi
do n2 = 1,nfo
arrdst(n2,n) = arrsrc(n2,n)
enddo
enddo
endif
else
write(s_logunit,F00) ' map%init illegal '//trim(map%init)
call shr_map_abort(subName//' ERROR map init')
endif
if (timing>0) call shr_timer_stop(t2)
!--- generate output array ---
if (trim(map%vect) == trim(shr_map_fs_vector)) then
if (timing>0) call shr_timer_start(t4)
do n=1,map%nwts
indi = map%isrc(n)
indo = map%idst(n)
wgt = map%wgts(n)
theta = map%xdst(indo) - map%xsrc(indi)
if (initnew(indo)) then
initnew(indo) = .false.
arrdst(1,indo) = wgt*( arrsrc(1,indi)*cos(theta) &
+arrsrc(2,indi)*sin(theta))
arrdst(2,indo) = wgt*(-arrsrc(1,indi)*sin(theta) &
+arrsrc(2,indi)*cos(theta))
else
arrdst(1,indo) = arrdst(1,indo) + wgt*( arrsrc(1,indi)*cos(theta) &
+arrsrc(2,indi)*sin(theta))
arrdst(2,indo) = arrdst(2,indo) + wgt*(-arrsrc(1,indi)*sin(theta) &
+arrsrc(2,indi)*cos(theta))
endif
enddo
if (timing>0) call shr_timer_stop(t4)
else
if (timing>0) call shr_timer_start(t5)
do n=1,map%nwts
indi = map%isrc(n)
indo = map%idst(n)
wgt = map%wgts(n)
if (initnew(indo)) then
initnew(indo) = .false.
do n2 = 1,nfo
arrdst(n2,indo) = arrsrc(n2,indi)*wgt
enddo
else
do n2 = 1,nfo
arrdst(n2,indo) = arrdst(n2,indo) + arrsrc(n2,indi)*wgt
enddo
endif
enddo
if (timing>0) call shr_timer_stop(t5)
endif
deallocate(initnew)
if (timing>0) call shr_timer_stop(t0)
end subroutine shr_map_mapDatam
!===============================================================================
!BOP ===========================================================================
!
! !IROUTINE: shr_map_mapDatanm -- maps arrays without map
!
! !DESCRIPTION:
! Maps arrays, don't save the map
! \newline
! call shr\_map\_mapData(Ain,Aout,Xs,Ys,Ms,Xd,Yd,Md,name,type,algo,vect,rc)
!
! !REVISION HISTORY:
! 2005-Mar-27 - T. Craig - first version
!
! !INTERFACE: ------------------------------------------------------------------
subroutine shr_map_mapDatanm(arrsrc,arrdst,Xsrc,Ysrc,Msrc,Xdst,Ydst,Mdst,name,type,algo,mask,vect,rc),1
implicit none
! !INPUT/OUTPUT PARAMETERS:
!--- map arrsrc to arrdst, each array is dimension (fields,grid index) ---
real(SHR_KIND_R8) ,intent(in) :: arrsrc(:,:) ! src array(fields,grid)
real(SHR_KIND_R8) ,intent(out):: arrdst(:,:) ! dst array(fields,grid)
real(SHR_KIND_R8) ,intent(in) :: Xsrc(:,:) ! lon of src grid
real(SHR_KIND_R8) ,intent(in) :: Ysrc(:,:) ! lat of src grid
integer(SHR_KIND_IN) ,intent(in) :: Msrc(:,:) ! mask of src grid
real(SHR_KIND_R8) ,intent(in) :: Xdst(:,:) ! lon of dst grid
real(SHR_KIND_R8) ,intent(in) :: Ydst(:,:) ! lat of dst grid
integer(SHR_KIND_IN) ,intent(in) :: Mdst(:,:) ! mask of dst grid
character(*) ,intent(in) :: name ! name
character(*) ,intent(in) :: type ! type
character(*) ,intent(in) :: algo ! algo
character(*) ,intent(in) :: mask ! mask
character(*) ,optional,intent(in) :: vect ! vect
integer(SHR_KIND_IN),optional,intent(out) :: rc ! error code
!EOP
!--- local ---
type(shr_map_mapType) :: map
integer(SHR_KIND_IN) :: lrc
!--- formats ---
character(*),parameter :: subName = "('shr_map_mapDatanm') "
character(*),parameter :: F00 = "('(shr_map_mapDatanm) ',a) "
!-------------------------------------------------------------------------------
lrc = 0
call shr_map_put(map,shr_map_fs_name,name,verify=.false.)
call shr_map_put(map,shr_map_fs_type,type,verify=.true.)
call shr_map_put(map,shr_map_fs_algo,algo,verify=.true.)
call shr_map_put(map,shr_map_fs_mask,mask,verify=.true.)
if (present(vect)) then
call shr_map_put(map,shr_map_fs_vect,vect,verify=.true.)
else
call shr_map_put(map,shr_map_fs_vect,'scalar',verify=.true.)
endif
call shr_map_mapSet(map,Xsrc,Ysrc,Msrc,Xdst,Ydst,Mdst,rc=lrc)
call shr_map_mapData(arrsrc,arrdst,map)
call shr_map_clean(map)
if (present(rc)) rc = lrc
end subroutine shr_map_mapDatanm
!===============================================================================
!BOP ===========================================================================
!
! !IROUTINE: shr_map_setAbort -- local interface to abort routine
!
! !DESCRIPTION:
! Set doabort flag for shr_map methods, true = call shr\_sys\_abort,
! false = write error message and continue
! \newline
! call shr\_map\_abort(subName//' ERROR: illegal option')
!
! !REVISION HISTORY:
! 2005-Mar-27 - T. Craig - first version
!
! !INTERFACE: ------------------------------------------------------------------
subroutine shr_map_setAbort(flag)
implicit none
! !INPUT/OUTPUT PARAMETERS:
logical,intent(in) :: flag
!EOP
!--- local ---
!--- formats ---
character(*),parameter :: subName = "('shr_map_setAbort') "
character(*),parameter :: F00 = "('(shr_map_setAbort) ',a) "
!-------------------------------------------------------------------------------
doabort = flag
end subroutine shr_map_setAbort
!===============================================================================
!BOP ===========================================================================
!
! !IROUTINE: shr_map_setDebug -- set local debug level
!
! !DESCRIPTION:
! Set debug level for shr_map methods, 0 = production
! \newline
! call shr\_map\_setDebug(2)
!
! !REVISION HISTORY:
! 2005-Apr-15 - T. Craig - first version
!
! !INTERFACE: ------------------------------------------------------------------
subroutine shr_map_setDebug(iflag)
implicit none
! !INPUT/OUTPUT PARAMETERS:
integer,intent(in) :: iflag
!EOP
!--- local ---
!--- formats ---
character(*),parameter :: subName = "('shr_map_setDebug') "
character(*),parameter :: F00 = "('(shr_map_setDebug) ',a) "
!-------------------------------------------------------------------------------
debug = iflag
end subroutine shr_map_setDebug
!===============================================================================
!BOP ===========================================================================
!
! !IROUTINE: shr_map_setDopole -- set dopole flag
!
! !DESCRIPTION:
! set dopole flag
! \newline
! call shr\_map\_setDopole(flag)
!
! !REVISION HISTORY:
! 2009-Jun-22 - T. Craig - first version
!
! !INTERFACE: ------------------------------------------------------------------
subroutine shr_map_setDopole(flag) 1
implicit none
! !INPUT/OUTPUT PARAMETERS:
logical, intent(in) :: flag
!EOP
!--- local ---
!--- formats ---
character(*),parameter :: subName = "('shr_map_setDopole') "
character(*),parameter :: F00 = "('(shr_map_setDopole) ',a) "
dopole = flag
end subroutine shr_map_setDopole
!===============================================================================
!XXBOP ===========================================================================
!
! !IROUTINE: shr_map_abort -- local interface to abort routine
!
! !DESCRIPTION:
! Local interface to abort routine. Depending on local flag, abort,
! either calls shr\_sys\_abort or writes abort message and continues.
! \newline
! call shr\_map\_abort(subName//' ERROR: illegal option')
!
! !REVISION HISTORY:
! 2005-Mar-27 - T. Craig - first version
!
! !INTERFACE: ------------------------------------------------------------------
subroutine shr_map_abort(string) 37,1
implicit none
! !INPUT/OUTPUT PARAMETERS:
character(*),optional,intent(in) :: string
!XXEOP
!--- local ---
character(shr_kind_CL) :: lstring
!--- formats ---
character(*),parameter :: subName = "('shr_map_abort') "
character(*),parameter :: F00 = "('(shr_map_abort) ',a) "
!-------------------------------------------------------------------------------
lstring = ''
if (present(string)) lstring = string
if (doabort) then
call shr_sys_abort(lstring)
else
write(s_logunit,F00) trim(lstring)
endif
end subroutine shr_map_abort
!===============================================================================
!XXBOP ===========================================================================
!
! !IROUTINE: shr_map_checkGrids_global -- local routine to check mapSet grids
!
! !DESCRIPTION:
! Local method to check grid arguments in shr\_map\_mapSet
! \newline
! call shr\_map\_checkGrids_global(Xs,Ys,Ms,Xd,Yd,Md,mymap)
!
! !REVISION HISTORY:
! 2005-Mar-27 - T. Craig - first version
!
! !INTERFACE: ------------------------------------------------------------------
subroutine shr_map_checkGrids_global(Xsrc,Ysrc,Msrc,Xdst,Ydst,Mdst,map,rc) 1,1
implicit none
! !INPUT/OUTPUT PARAMETERS:
real(SHR_KIND_R8) ,intent(in) :: Xsrc(:,:) ! src lat
real(SHR_KIND_R8) ,intent(in) :: Ysrc(:,:) ! src lon
integer(SHR_KIND_IN) ,intent(in) :: Msrc(:,:) ! src mask
real(SHR_KIND_R8) ,intent(in) :: Xdst(:,:) ! dst lat
real(SHR_KIND_R8) ,intent(in) :: Ydst(:,:) ! dst lon
integer(SHR_KIND_IN) ,intent(in) :: Mdst(:,:) ! dst mask
type(shr_map_mapType),intent(in) :: map ! map
integer(SHR_KIND_IN),optional,intent(out) :: rc ! error code
!XXEOP
!--- local ---
integer(SHR_KIND_IN) :: i,j,nis,njs,nid,njd,ncnt
logical :: error,flag
!--- formats ---
character(*),parameter :: subName = "('shr_map_checkGrids_global') "
character(*),parameter :: F00 = "('(shr_map_checkGrids_global) ',a) "
character(*),parameter :: F01 = "('(shr_map_checkGrids_global) ',a,2i8) "
character(*),parameter :: F02 = "('(shr_map_checkGrids_global) ',a,4i8) "
character(*),parameter :: F03 = "('(shr_map_checkGrids_global) ',a,2g20.13) "
character(*),parameter :: F04 = "('(shr_map_checkGrids_global) ',a,i8,a,i8) "
character(*),parameter :: F05 = "('(shr_map_checkGrids_global) ',a,i8,2g20.13) "
character(*),parameter :: F06 = "('(shr_map_checkGrids_global) ',a,2i8,2g20.13) "
!-------------------------------------------------------------------------------
error = .false.
if (present(rc)) rc = 0
!--- get size of X arrays
nis = size(Xsrc,1)
njs = size(Xsrc,2)
nid = size(Xdst,1)
njd = size(Xdst,2)
!--- check array size consistency for src and dst
if (size(Ysrc,1) /= nis) then
write(s_logunit,F01) 'ERROR Xsrc,Ysrc i-dim mismatch',nis,size(Ysrc,1)
error = .true.
endif
if (size(Ysrc,2) /= njs) then
write(s_logunit,F01) 'ERROR Xsrc,Ysrc j-dim mismatch',njs,size(Ysrc,2)
error = .true.
endif
if (size(Msrc,1) /= nis) then
write(s_logunit,F01) 'ERROR Xsrc,Msrc i-dim mismatch',nis,size(Msrc,1)
error = .true.
endif
if (size(Msrc,2) /= njs) then
write(s_logunit,F01) 'ERROR Xsrc,Msrc j-dim mismatch',njs,size(Msrc,2)
error = .true.
endif
if (size(Ydst,1) /= nid) then
write(s_logunit,F01) 'ERROR Xdst,Ydst i-dim mismatch',nid,size(Ydst,1)
error = .true.
endif
if (size(Ydst,2) /= njd) then
write(s_logunit,F01) 'ERROR Xdst,Ydst j-dim mismatch',njd,size(Ydst,2)
error = .true.
endif
if (size(Mdst,1) /= nid) then
write(s_logunit,F01) 'ERROR Xdst,Mdst i-dim mismatch',nid,size(Mdst,1)
error = .true.
endif
if (size(Mdst,2) /= njd) then
write(s_logunit,F01) 'ERROR Xdst,Mdst j-dim mismatch',njd,size(Mdst,2)
error = .true.
endif
!--- fill type must have same grid size on src and dst ---
if (trim(map%type) == trim(shr_map_fs_fill) .or. &
trim(map%type) == trim(shr_map_fs_cfill)) then
if (nis*njs /= nid*njd) then
write(s_logunit,F02) 'ERROR: fill type, src/dst sizes ',nis*njs,nid*njd
error = .true.
endif
endif
!--- write min/max or X, Y and M count ---
if (debug > 1 .and. s_loglev > 0) then
write(s_logunit,F03) ' Xsrc min/max ',minval(Xsrc),maxval(Xsrc)
write(s_logunit,F03) ' Ysrc min/max ',minval(Ysrc),maxval(Ysrc)
write(s_logunit,F03) ' Xdst min/max ',minval(Xdst),maxval(Xdst)
write(s_logunit,F03) ' Ydst min/max ',minval(Ydst),maxval(Ydst)
endif
ncnt = 0
do j=1,njs
do i=1,nis
if (Msrc(i,j) == 0) ncnt = ncnt + 1
enddo
enddo
if (debug > 1 .and. s_loglev > 0) write(s_logunit,F04) ' Msrc mask T ',nis*njs-ncnt,' of ',nis*njs
ncnt = 0
do j=1,njd
do i=1,nid
if (Mdst(i,j) == 0) ncnt = ncnt + 1
enddo
enddo
if (debug > 1 .and. s_loglev > 0) write(s_logunit,F04) ' Mdst mask T ',nid*njd-ncnt,' of ',nid*njd
if (trim(map%algo) == trim(shr_map_fs_bilinear)) then
!--- check that Xsrc is monotonically increasing for bilinear ---
flag = .false.
i = 1
do while (i < nis .and. .not.flag)
if (Xsrc(i+1,1) <= Xsrc(i,1)) then
write(s_logunit,F05) 'ERROR Xsrc not increasing ',i,Xsrc(i+1,1),Xsrc(i,1)
flag = .true.
error = .true.
endif
i = i+1
enddo
!--- check that Ysrc is monotonically increasing for bilinear ---
flag = .false.
j = 1
do while (j < njs .and. .not.flag)
if (Ysrc(1,j+1) <= Ysrc(1,j)) then
write(s_logunit,F05) 'ERROR Ysrc not increasing ',i,Ysrc(1,j+1),Ysrc(1,j)
flag = .true.
error = .true.
endif
j = j+1
enddo
!--- check that Xsrc and Ysrc are regular lat/lon grids for bilinear
flag = .false.
i = 1
do while (i < nis .and. .not.flag)
j = 2
do while (j < njs .and. .not.flag)
if (abs(Xsrc(i,j)-Xsrc(i,1)) > eps) then
write(s_logunit,F06) ' ERROR Xsrc not regular lat,lon ',i,j, &
Xsrc(i,j),Xsrc(1,j)
flag = .true.
error = .true.
endif
j = j+1
enddo
i = i+1
enddo
flag = .false.
j = 1
do while (j < njs .and. .not.flag)
i = 2
do while (i < nis .and. .not.flag)
if (abs(Ysrc(i,j)-Ysrc(1,j)) > eps) then
write(s_logunit,F06) ' ERROR Ysrc not regular lat,lon ',i,j, &
Ysrc(i,j),Ysrc(1,j)
flag = .true.
error = .true.
endif
i = i+1
enddo
j = j+1
enddo
endif
if (error) then
call shr_map_abort(subName//' ERROR ')
if (present(rc)) rc = 1
endif
end subroutine shr_map_checkGrids_global
!===============================================================================
!XXBOP ===========================================================================
!
! !IROUTINE: shr_map_checkGrids_dest -- local routine to check mapSet grids
!
! !DESCRIPTION:
! Local method to check grid arguments in shr\_map\_mapSet
! \newline
! call shr\_map\_checkGrids_dest(Xs,Ys,Ms,Xd,Yd,Md,mymap)
!
! !REVISION HISTORY:
! 2005-Mar-27 - T. Craig - first version
!
! !INTERFACE: ------------------------------------------------------------------
subroutine shr_map_checkGrids_dest(Xsrc,Ysrc,Msrc,Xdst,Ydst,Mdst,map,rc) 1,1
implicit none
! !INPUT/OUTPUT PARAMETERS:
real(SHR_KIND_R8) ,intent(in) :: Xsrc(:,:) ! src lat
real(SHR_KIND_R8) ,intent(in) :: Ysrc(:,:) ! src lon
integer(SHR_KIND_IN) ,intent(in) :: Msrc(:,:) ! src mask
real(SHR_KIND_R8) ,intent(in) :: Xdst(:) ! dst lat
real(SHR_KIND_R8) ,intent(in) :: Ydst(:) ! dst lon
integer(SHR_KIND_IN) ,intent(in) :: Mdst(:) ! dst mask
type(shr_map_mapType),intent(in) :: map ! map
integer(SHR_KIND_IN),optional,intent(out) :: rc ! error code
!XXEOP
!--- local ---
integer(SHR_KIND_IN) :: i,j,nis,njs,nid,njd,ncnt
logical :: error,flag
!--- formats ---
character(*),parameter :: subName = "('shr_map_checkGrids_dest') "
character(*),parameter :: F00 = "('(shr_map_checkGrids_dest) ',a) "
character(*),parameter :: F01 = "('(shr_map_checkGrids_dest) ',a,2i8) "
character(*),parameter :: F02 = "('(shr_map_checkGrids_dest) ',a,4i8) "
character(*),parameter :: F03 = "('(shr_map_checkGrids_dest) ',a,2g20.13) "
character(*),parameter :: F04 = "('(shr_map_checkGrids_dest) ',a,i8,a,i8) "
character(*),parameter :: F05 = "('(shr_map_checkGrids_dest) ',a,i8,2g20.13) "
character(*),parameter :: F06 = "('(shr_map_checkGrids_dest) ',a,2i8,2g20.13) "
!-------------------------------------------------------------------------------
error = .false.
if (present(rc)) rc = 0
!--- get size of X arrays
nis = size(Xsrc,1)
njs = size(Xsrc,2)
nid = size(Xdst,1)
njd = 1
!--- check array size consistency for src and dst
if (size(Ysrc,1) /= nis) then
write(s_logunit,F01) 'ERROR Xsrc,Ysrc i-dim mismatch',nis,size(Ysrc,1)
error = .true.
endif
if (size(Ysrc,2) /= njs) then
write(s_logunit,F01) 'ERROR Xsrc,Ysrc j-dim mismatch',njs,size(Ysrc,2)
error = .true.
endif
if (size(Msrc,1) /= nis) then
write(s_logunit,F01) 'ERROR Xsrc,Msrc i-dim mismatch',nis,size(Msrc,1)
error = .true.
endif
if (size(Msrc,2) /= njs) then
write(s_logunit,F01) 'ERROR Xsrc,Msrc j-dim mismatch',njs,size(Msrc,2)
error = .true.
endif
if (size(Ydst,1) /= nid) then
write(s_logunit,F01) 'ERROR Xdst,Ydst i-dim mismatch',nid,size(Ydst,1)
error = .true.
endif
if (size(Mdst,1) /= nid) then
write(s_logunit,F01) 'ERROR Xdst,Mdst i-dim mismatch',nid,size(Mdst,1)
error = .true.
endif
!--- tcraig, can't check this with dest mapset ---
! !--- fill type must have same grid size on src and dst ---
! if (trim(map%type) == trim(shr_map_fs_fill) .or. &
! trim(map%type) == trim(shr_map_fs_cfill)) then
! if (nis*njs /= nid*njd) then
! write(s_logunit,F02) 'ERROR: fill type, src/dst sizes ',nis*njs,nid*njd
! error = .true.
! endif
! endif
!--- write min/max or X, Y and M count ---
if (debug > 1 .and. s_loglev > 0) then
write(s_logunit,F03) ' Xsrc min/max ',minval(Xsrc),maxval(Xsrc)
write(s_logunit,F03) ' Ysrc min/max ',minval(Ysrc),maxval(Ysrc)
write(s_logunit,F03) ' Xdst min/max ',minval(Xdst),maxval(Xdst)
write(s_logunit,F03) ' Ydst min/max ',minval(Ydst),maxval(Ydst)
endif
ncnt = 0
do j=1,njs
do i=1,nis
if (Msrc(i,j) == 0) ncnt = ncnt + 1
enddo
enddo
if (debug > 1 .and. s_loglev > 0) write(s_logunit,F04) ' Msrc mask T ',nis*njs-ncnt,' of ',nis*njs
ncnt = 0
do i=1,nid
if (Mdst(i) == 0) ncnt = ncnt + 1
enddo
if (debug > 1 .and. s_loglev > 0) write(s_logunit,F04) ' Mdst mask T ',nid*njd-ncnt,' of ',nid*njd
if (trim(map%algo) == trim(shr_map_fs_bilinear)) then
!--- check that Xsrc is monotonically increasing for bilinear ---
flag = .false.
i = 1
do while (i < nis .and. .not.flag)
if (Xsrc(i+1,1) <= Xsrc(i,1)) then
write(s_logunit,F05) 'ERROR Xsrc not increasing ',i,Xsrc(i+1,1),Xsrc(i,1)
flag = .true.
error = .true.
endif
i = i+1
enddo
!--- check that Ysrc is monotonically increasing for bilinear ---
flag = .false.
j = 1
do while (j < njs .and. .not.flag)
if (Ysrc(1,j+1) <= Ysrc(1,j)) then
write(s_logunit,F05) 'ERROR Ysrc not increasing ',i,Ysrc(1,j+1),Ysrc(1,j)
flag = .true.
error = .true.
endif
j = j+1
enddo
!--- check that Xsrc and Ysrc are regular lat/lon grids for bilinear
flag = .false.
i = 1
do while (i < nis .and. .not.flag)
j = 2
do while (j < njs .and. .not.flag)
if (abs(Xsrc(i,j)-Xsrc(i,1)) > eps) then
write(s_logunit,F06) ' ERROR Xsrc not regular lat,lon ',i,j, &
Xsrc(i,j),Xsrc(1,j)
flag = .true.
error = .true.
endif
j = j+1
enddo
i = i+1
enddo
flag = .false.
j = 1
do while (j < njs .and. .not.flag)
i = 2
do while (i < nis .and. .not.flag)
if (abs(Ysrc(i,j)-Ysrc(1,j)) > eps) then
write(s_logunit,F06) ' ERROR Ysrc not regular lat,lon ',i,j, &
Ysrc(i,j),Ysrc(1,j)
flag = .true.
error = .true.
endif
i = i+1
enddo
j = j+1
enddo
endif
if (error) then
call shr_map_abort(subName//' ERROR ')
if (present(rc)) rc = 1
endif
end subroutine shr_map_checkGrids_dest
!===============================================================================
!XXBOP ===========================================================================
!
! !IROUTINE: shr_map_checkWgts_global -- checks weights
!
! !DESCRIPTION:
! Checks weights in map for validity
! \newline
! call shr\_map\_checkWgts_global(Ms,Md,mymap)
!
! !REVISION HISTORY:
! 2005-Mar-27 - T. Craig - first version
!
! !INTERFACE: ------------------------------------------------------------------
subroutine shr_map_checkWgts_global(Msrc,Mdst,map) 1,5
implicit none
! !INPUT/OUTPUT PARAMETERS:
integer(SHR_KIND_IN) ,intent(in) :: Msrc(:,:) ! src mask
integer(SHR_KIND_IN) ,intent(in) :: Mdst(:,:) ! dst mask
type(shr_map_mapType),intent(in) :: map ! map
!XXEOP
!--- local ---
integer(SHR_KIND_IN) :: i,j,nis,njs,nid,njd,n
integer(SHR_KIND_IN) :: ic1,ic2,ic3,ic4,ic5 ! counters
logical :: error
real(SHR_KIND_R8),allocatable :: Csrc(:,:)
real(SHR_KIND_R8),allocatable :: Cdst(:,:)
!--- formats ---
character(*),parameter :: subName = "('shr_map_checkWgts_global') "
character(*),parameter :: F00 = "('(shr_map_checkWgts_global) ',a) "
character(*),parameter :: F01 = "('(shr_map_checkWgts_global) ',a,i8) "
character(*),parameter :: F02 = "('(shr_map_checkWgts_global) ',a,3i8) "
character(*),parameter :: F03 = "('(shr_map_checkWgts_global) ',a,i8,a) "
!-------------------------------------------------------------------------------
error = .false.
if (debug > 0) call shr_map_print(map)
if (map%nwts < 1) then
if (s_loglev > 0) write(s_logunit,F00) 'WARNING map size is zero'
endif
if (size(map%wgts) /= map%nwts .or. &
size(map%isrc) /= map%nwts .or. &
size(map%idst) /= map%nwts) then
call shr_map_abort(subName//'ERROR sizes inconsistent')
endif
!--- get size of X arrays
nis = size(Msrc,1)
njs = size(Msrc,2)
nid = size(Mdst,1)
njd = size(Mdst,2)
allocate(Csrc(nis,njs))
allocate(Cdst(nid,njd))
Csrc = c0
Cdst = c0
do n = 1,map%nwts
call shr_map_1dto2d(map%isrc(n),nis,njs,i,j)
Csrc(i,j) = c1
call shr_map_1dto2d(map%idst(n),nid,njd,i,j)
Cdst(i,j) = Cdst(i,j) + map%wgts(n)
enddo
ic1 = 0
ic2 = 0
ic3 = 0
ic4 = 0
ic5 = 0
do j=1,njs
do i=1,nis
if (Msrc(i,j) /= 0) then ! live src pt
if (abs(Csrc(i,j)-c1) < eps) then
ic1 = ic1 + 1 ! in use
else
ic2 = ic2 + 1 ! not used
endif
else ! dead src pt
if (abs(Csrc(i,j)-c1) < eps) then
ic3 = ic3 + 1 ! in use
else
ic5 = ic5 + 1 ! not used
endif
endif
enddo
enddo
! if (ic3 > 0) error = .true.
if (debug > 0 .and. s_loglev > 0) then
write(s_logunit,F01) ' total number of SRC points : ',nis*njs
write(s_logunit,F01) ' wgts from SRC TRUE points; used : ',ic1
write(s_logunit,F01) ' wgts from SRC TRUE points; not used : ',ic2
write(s_logunit,F01) ' wgts from SRC FALSE points; used : ',ic3
write(s_logunit,F01) ' wgts from SRC FALSE points; not used : ',ic5
endif
ic1 = 0
ic2 = 0
ic3 = 0
ic4 = 0
ic5 = 0
do j=1,njd
do i=1,nid
if (Mdst(i,j) /= 0) then ! wgts should sum to one
if (abs(Cdst(i,j)-c1) < eps) then
ic1 = ic1 + 1 ! wgts sum to one
else
ic2 = ic2 + 1 ! invalid wgts
endif
else ! wgts should sum to one or zero
if (abs(Cdst(i,j)-c1) < eps) then
ic3 = ic3 + 1 ! wgts sum to one
elseif (abs(Cdst(i,j)) < eps) then
ic4 = ic4 + 1 ! wgts sum to zero
else
ic5 = ic5 + 1 ! invalid wgts
endif
endif
enddo
enddo
! if (ic2 > 0) error = .true.
! if (ic5 > 0) error = .true.
if (debug > 0 .and. s_loglev > 0) then
write(s_logunit,F01) ' total number of DST points : ',nid*njd
write(s_logunit,F01) ' sum wgts for DST TRUE points; one : ',ic1
if (ic2 > 0) then
write(s_logunit,F03) ' sum wgts for DST TRUE points; not : ',ic2,' **-WARNING-**'
else
write(s_logunit,F01) ' sum wgts for DST TRUE points; not : ',ic2
endif
write(s_logunit,F01) ' sum wgts for DST FALSE points; one : ',ic3
write(s_logunit,F01) ' sum wgts for DST FALSE points; zero : ',ic4
write(s_logunit,F01) ' sum wgts for DST FALSE points; not : ',ic5
endif
deallocate(Csrc)
deallocate(Cdst)
if (error) call shr_map_abort(subName//' ERROR invalid weights')
end subroutine shr_map_checkWgts_global
!===============================================================================
!XXBOP ===========================================================================
!
! !IROUTINE: shr_map_getWts -- local code that sets weights for a point
!
! !DESCRIPTION:
! Local code that sets weights for a point. Executes searches
! and computes weights. For bilinear remap for example.
!
! !REMARKS:
! Assumes Xsrc,Ysrc are regular lat/lon grids, monotonicallly increasing
! on constant latitude and longitude lines.
! Assumes Xdst,Ydst,Xsrc,Ysrc are all either radians or degrees
!
! !REVISION HISTORY:
! 2005-Mar-27 - T. Craig - first version
!
! !INTERFACE: ------------------------------------------------------------------
subroutine shr_map_getWts(Xdst,Ydst,Xsrc,Ysrc,pti,ptj,ptw,pnum,units) 2,7
implicit none
! !INPUT/OUTPUT PARAMETERS:
!XXEOP
real(SHR_KIND_R8) ,intent(in) :: Xdst,Ydst
real(SHR_KIND_R8) ,intent(in) :: Xsrc(:,:),Ysrc(:,:)
integer(SHR_KIND_IN),intent(out):: pti(:),ptj(:)
real(SHR_KIND_R8) ,intent(out):: ptw(:)
integer(SHR_KIND_IN),intent(out):: pnum
character(len=*),optional,intent(in) :: units
!--- local ---
integer(SHR_KIND_IN) :: isize,jsize ! array sizes
integer(SHR_KIND_IN) :: i,j ! indices
integer(SHR_KIND_IN) :: n ! do loop counter
integer(SHR_KIND_IN) :: il,ir ! index of i left/mid/right
integer(SHR_KIND_IN) :: jl,ju ! index of j lower/mid/upper
integer(SHR_KIND_IN) :: pmax ! size of pti,ptj,ptw
real(SHR_KIND_R8) :: xsl,xsr ! value of Xsrc, left/right
real(SHR_KIND_R8) :: ysl,ysu ! value of Ysrc, left/right
real(SHR_KIND_R8) :: xd,yd ! value of Xdst,Ydst
real(SHR_KIND_R8) :: dx,dy,dx1,dy1 ! some d_lengths for weights calc
real(SHR_KIND_R8) :: csize ! circle angle/radians
real(SHR_KIND_R8) :: rmin,rmax ! min/max
real(SHR_KIND_R8) :: cpole ! the r8 lat value of the pole
integer(SHR_KIND_IN) :: pole ! 0=no, 1=north, 2=south
!--- formats ---
character(*),parameter :: subName = "('shr_map_getWts') "
character(*),parameter :: F00 = "('(shr_map_getWts) ',a) "
character(*),parameter :: F02 = "('(shr_map_getWts) ',a,4g20.13) "
character(*),parameter :: F03 = "('(shr_map_getWts) ',a,2g20.13) "
character(*),parameter :: F04 = "('(shr_map_getWts) ',a,4i8) "
character(*),parameter :: F05 = "('(shr_map_getWts) ',a,3g20.13) "
!-------------------------------------------------------------------------------
pmax = size(pti,1)
!--- is lat/lon degrees or radians? needed for X wraparound ---
if (present(units)) then
if (trim(units) == 'degrees') then
csize = 360._SHR_KIND_R8
elseif (trim(units) == 'radians') then
csize = c2*pi
else
call shr_sys_abort(subName//' ERROR in optional units = '//trim(units))
endif
else
csize = 360._SHR_KIND_R8
if (shr_map_checkRad(Ysrc)) csize = c2*pi
endif
isize = size(Xsrc,1)
jsize = size(Xsrc,2)
pti = 0
ptj = 0
ptw = c0
cpole = csize/(c2*c2)
xd = Xdst
yd = Ydst
if (yd > cpole + 1.0e-3 .or. &
yd < -cpole - 1.0e-3) then
write(s_logunit,*) trim(subname),' ERROR: yd outside bounds ',yd
call shr_map_abort(subName//' ERROR yd outside 90 degree bounds')
endif
if (yd > cpole) yd = cpole
if (yd < -cpole) yd = -cpole
call shr_map_find4corners(Xdst,yd,Xsrc,Ysrc,il,ir,jl,ju)
!--- bilinear ---
pnum = 4
pole = 0
xsl = Xsrc(il,1)
xsr = Xsrc(ir,1)
ysl = Ysrc(1,jl)
ysu = Ysrc(1,ju)
if (Xdst < Xsrc(1,1) .or. Xdst > Xsrc(isize,1)) then
xsl = mod(Xsrc(il,1),csize)
xsr = mod(Xsrc(ir,1),csize)
xd = mod(Xdst ,csize)
if (xsl > xd) xsl = xsl - csize
if (xsr < xd) xsr = xsr + csize
endif
if (yd > Ysrc(1,jsize)) then
if (dopole) then
pnum = isize+2
pole = 1
endif
ysu = cpole
elseif (yd < Ysrc(1,1)) then
if (dopole) then
pnum = isize+2
pole = 2
endif
ysl = -cpole
endif
!--- compute dx1,dy1; distance from src(1) to dst
dx = (xsr-xsl)
dy = (ysu-ysl)
dx1 = ( xd-xsl)
dy1 = ( yd-Ysl)
if (dx1 > dx .and. dx1-dx < 1.0e-7 ) dx1 = dx
if (dy1 > dy .and. dy1-dy < 1.0e-7 ) dy1 = dy
if (dx <= c0 .or. dy <= c0 .or. dx1 > dx .or. dy1 > dy) then
write(s_logunit,*) ' '
write(s_logunit,F02) 'ERROR in dx,dy: ',dx1,dx,dy1,dy
write(s_logunit,F03) ' dst: ',Xdst,Ydst
write(s_logunit,F04) ' ind: ',il,ir,jl,ju
write(s_logunit,F02) ' dis: ',dx1,dx,dy1,dy
write(s_logunit,F05) ' x3 : ',xsl,xd,xsr
write(s_logunit,F05) ' y3 : ',ysl,yd,ysu
write(s_logunit,*) ' '
call shr_map_abort(subName//' ERROR in dx,dy calc')
stop
return
endif
dx1 = dx1 / dx
dy1 = dy1 / dy
if (pnum > pmax) then
call shr_sys_abort(subName//' ERROR pti not big enough')
endif
if (pole == 0) then ! bilinear
pti(1) = il
pti(2) = ir
pti(3) = il
pti(4) = ir
ptj(1) = jl
ptj(2) = jl
ptj(3) = ju
ptj(4) = ju
ptw(1) = (c1-dx1)*(c1-dy1)
ptw(2) = ( dx1)*(c1-dy1)
ptw(3) = (c1-dx1)*( dy1)
ptw(4) = ( dx1)*( dy1)
elseif (pole == 1) then ! north pole
pti(1) = il
pti(2) = ir
ptj(1) = jl
ptj(2) = jl
ptw(1) = (c1-dx1)*(c1-dy1)
ptw(2) = ( dx1)*(c1-dy1)
do n=1,isize
pti(2+n) = n
ptj(2+n) = ju
ptw(2+n) = (dy1)/real(isize,SHR_KIND_R8)
enddo
elseif (pole == 2) then ! south pole
pti(1) = il
pti(2) = ir
ptj(1) = ju
ptj(2) = ju
ptw(1) = (c1-dx1)*( dy1)
ptw(2) = ( dx1)*( dy1)
do n=1,isize
pti(2+n) = n
ptj(2+n) = jl
ptw(2+n) = (c1-dy1)/real(isize,SHR_KIND_R8)
enddo
else
write(s_logunit,F00) ' ERROR illegal pnum situation '
call shr_map_abort(subName//' ERROR illegal pnum situation')
endif
end subroutine shr_map_getWts
!===============================================================================
subroutine shr_map_find4corners(Xdst,Ydst,Xsrc,Ysrc,il,ir,jl,ju) 2
! finds 4 corner points surrounding dst in src
! returns left, right, lower, and upper i and j index
implicit none
! !INPUT/OUTPUT PARAMETERS:
real(SHR_KIND_R8) ,intent(in) :: Xdst,Ydst
real(SHR_KIND_R8) ,intent(in) :: Xsrc(:,:),Ysrc(:,:)
integer(SHR_KIND_IN),intent(out):: il,ir,jl,ju
!--- local ---
integer(SHR_KIND_IN) :: isize,jsize
integer(SHR_KIND_IN) :: im,jm
!--- formats ---
character(*),parameter :: subName = "('shr_map_find4corners') "
character(*),parameter :: F00 = "('(shr_map_find4corners) ',a,2i8) "
!-------------------------------------------------------------------------------
isize = size(Xsrc,1)
jsize = size(Xsrc,2)
if (Xdst < Xsrc(1,1) .or. Xdst > Xsrc(isize,1)) then
il = isize
ir = 1
else
!--- find i index where Xsrc(i) <= Xdst < Xsrc(i+1) ---
il = 1
ir = isize
do while (ir-il > 1)
im = (ir+il)/2
if (Xdst >= Xsrc(im,1)) then
il = im
else
ir = im
endif
enddo
endif
if (Ydst > Ysrc(1,jsize)) then
jl = jsize
ju = jsize
elseif (Ydst < Ysrc(1,1)) then
jl = 1
ju = 1
else
!--- find j index where Ysrc(j) <= Ydst < Ysrc(j+1) ---
jl = 1
ju = jsize
do while (ju-jl > 1)
jm = (ju+jl)/2
if (Ydst >= Ysrc(1,jm)) then
jl = jm
else
ju = jm
endif
enddo
endif
end subroutine shr_map_find4corners
!===============================================================================
subroutine shr_map_findnn(Xdst,Ydst,Xsrc,Ysrc,Msrc,inn,jnn) 4,1
! finds point in src nearest to dst, returns inn,jnn src index
! searches using Msrc active points only
implicit none
! !INPUT/OUTPUT PARAMETERS:
real(SHR_KIND_R8) ,intent(in) :: Xdst,Ydst
real(SHR_KIND_R8) ,intent(in) :: Xsrc(:,:),Ysrc(:,:)
integer(SHR_KIND_IN),intent(in) :: Msrc(:,:)
integer(SHR_KIND_IN),intent(out):: inn,jnn
!--- local ---
integer(SHR_KIND_IN) :: isize,jsize
integer(SHR_KIND_IN) :: i,j
real(SHR_KIND_R8) :: dnn,dist
!--- formats ---
character(*),parameter :: subName = "('shr_map_findnn') "
character(*),parameter :: F00 = "('(shr_map_findnn) ',a,2i8) "
!-------------------------------------------------------------------------------
isize = size(Xsrc,1)
jsize = size(Xsrc,2)
inn = -1
jnn = -1
dnn = -1._SHR_KIND_R8
do j=1,jsize
do i=1,isize
if (Msrc(i,j) /= 0) then
dist = shr_map_finddist(Xdst,Ydst,Xsrc(i,j),Ysrc(i,j))
if (dist < dnn .or. inn < 0) then
dnn = dist
inn = i
jnn = j
endif
endif
enddo
enddo
end subroutine shr_map_findnn
!===============================================================================
subroutine shr_map_findnnon(dir,Xdst,Ydst,Xsrc,Ysrc,Msrc,inn,jnn) 8,8
! finds point in src nearest to dst searching i dir first
! returns inn,jnn src index
! searches using Msrc active points only
implicit none
! !INPUT/OUTPUT PARAMETERS:
character(*) ,intent(in) :: dir
real(SHR_KIND_R8) ,intent(in) :: Xdst,Ydst
real(SHR_KIND_R8) ,intent(in) :: Xsrc(:,:),Ysrc(:,:)
integer(SHR_KIND_IN),intent(in) :: Msrc(:,:)
integer(SHR_KIND_IN),intent(out):: inn,jnn
!--- local ---
integer(SHR_KIND_IN) :: isize,jsize
integer(SHR_KIND_IN) :: il,ir,jl,ju
integer(SHR_KIND_IN) :: n,i,j
integer(SHR_KIND_IN) :: is,js
integer(SHR_KIND_IN) :: i2,j2
real(SHR_KIND_R8) :: dnn,dist,ds
!--- formats ---
character(*),parameter :: subName = "('shr_map_findnnon') "
character(*),parameter :: F00 = "('(shr_map_findnnon) ',a,2i8) "
!-------------------------------------------------------------------------------
isize = size(Xsrc,1)
jsize = size(Xsrc,2)
!--- find 4 corner points
call shr_map_find4corners(Xdst,Ydst,Xsrc,Ysrc,il,ir,jl,ju)
!--- find closest of 4 corner points to dst, set that to is,js
is = il
js = jl
ds = shr_map_finddist(Xdst,Ydst,Xsrc(il,jl),Ysrc(il,jl))
dist = shr_map_finddist(Xdst,Ydst,Xsrc(ir,jl),Ysrc(ir,jl))
if (dist < ds) then
is = ir
js = jl
ds = dist
endif
dist = shr_map_finddist(Xdst,Ydst,Xsrc(il,ju),Ysrc(il,ju))
if (dist < ds) then
is = il
js = ju
ds = dist
endif
dist = shr_map_finddist(Xdst,Ydst,Xsrc(ir,ju),Ysrc(ir,ju))
if (dist < ds) then
is = ir
js = ju
ds = dist
endif
inn = -1
jnn = -1
dnn = -1._SHR_KIND_R8
i2 = 0
j2 = 0
if (trim(dir) == 'i') then
!--- search biased over i ---
do while (inn < 0 .and. j2 < jsize)
do n=1,2
if (n == 1) j = min(js + j2,jsize)
if (n == 2) j = max(js - j2,1)
do i=1,isize
if (Msrc(i,j) /= 0) then
dist = shr_map_finddist(Xdst,Ydst,Xsrc(i,j),Ysrc(i,j))
if (dist < dnn .or. inn < 0) then
dnn = dist
inn = i
jnn = j
endif
endif
enddo
enddo
j2 = j2 + 1
enddo
elseif (trim(dir) == 'j') then
!--- search biased over j ---
do while (inn < 0 .and. i2 < isize)
do n=1,2
if (n == 1) i = min(is + i2,isize)
if (n == 2) i = max(is - i2,1)
do j=1,jsize
if (Msrc(i,j) /= 0) then
dist = shr_map_finddist(Xdst,Ydst,Xsrc(i,j),Ysrc(i,j))
if (dist < dnn .or. inn < 0) then
dnn = dist
inn = i
jnn = j
endif
endif
enddo
enddo
i2 = i2 + 1
enddo
else
call shr_map_abort(subName//' ERROR illegal dir '//trim(dir))
endif
end subroutine shr_map_findnnon
!===============================================================================
real(SHR_KIND_R8) function shr_map_finddist(Xdst,Ydst,Xsrc,Ysrc) 7
! x,y distance computation
implicit none
real(SHR_KIND_R8),intent(in) :: Xdst,Ydst,Xsrc,Ysrc
character(*),parameter :: subName = "('shr_map_finddist') "
!-------------------------------------------------------------------------------
shr_map_finddist = sqrt((Ydst-Ysrc)**2 + (Xdst-Xsrc)**2)
end function shr_map_finddist
!===============================================================================
logical function shr_map_checkRad(Grid) 3
! check if grid is rad or degree
implicit none
real(SHR_KIND_R8),intent(in) :: Grid(:,:)
character(*),parameter :: subName = "('shr_map_checkRad') "
real(SHR_KIND_R8) :: rmin,rmax
!-------------------------------------------------------------------------------
shr_map_checkRad = .false.
rmin = minval(Grid)
rmax = maxval(Grid)
if ((rmax - rmin) < 1.01_SHR_KIND_R8*c2*pi) shr_map_checkRad = .true.
end function shr_map_checkRad
!===============================================================================
subroutine shr_map_1dto2d(gid,ni,nj,i,j) 14,1
! convert from a 1d index system to a 2d index system
! gid is 1d index; ni,nj are 2d grid size; i,j are local 2d index
implicit none
integer(SHR_KIND_IN),intent(in) :: gid,ni,nj
integer(SHR_KIND_IN),intent(out):: i,j
character(*),parameter :: subName = "('shr_map_1dto2d') "
character(*),parameter :: F01 = "('(shr_map_1dto2d) ',a,3i8)"
!-------------------------------------------------------------------------------
if (gid < 1 .or. gid > ni*nj) then
write(s_logunit,F01) 'ERROR: illegal gid ',gid,ni,nj
call shr_map_abort(subName//' ERROR')
endif
j = (gid-1)/ni+1
i = mod(gid-1,ni)+1
end subroutine shr_map_1dto2d
!===============================================================================
subroutine shr_map_2dto1d(gid,ni,nj,i,j) 17,1
! convert from a 2d index system to a 1d index system
! gid is 1d index; ni,nj are 2d grid size; i,j are local 2d index
implicit none
integer(SHR_KIND_IN),intent(in) :: ni,nj,i,j
integer(SHR_KIND_IN),intent(out):: gid
character(*),parameter :: subName = "('shr_map_2dto1d') "
character(*),parameter :: F01 = "('(shr_map_2dto1d) ',a,4i8)"
!-------------------------------------------------------------------------------
if (i < 1 .or. i > ni .or. j < 1 .or. j > nj) then
write(s_logunit,F01) 'ERROR: illegal i,j ',i,ni,j,nj
call shr_map_abort(subName//' ERROR')
endif
gid = (j-1)*ni + i
end subroutine shr_map_2dto1d
!===============================================================================
!===============================================================================
end module shr_map_mod