module m_cvgridLookup !$$$ subprogram documentation block ! . . . . ! subprogram: module m_cvgridLookup ! prgmmr: j guo ! org: NASA/GSFC, Global Modeling and Assimilation Office, 610.3 ! date: 2015-08-18 ! ! abstract: cv (control-vector) grid utilities for obs. partitioning ! ! program history log: ! 2015-08-18 j guo - initial implementation, from m_obsNode.F90 ! 2016-06-22 j guo - added []_sdget() to support m_latlonRange ! . refined []_islocal() algorithm and related module ! variables, ilon_lbound, ilon_ubound, etc., to support ! regional/non-periodic subdomain grid types. ! ! input argument list: see Fortran 90 style document below ! ! output argument list: see Fortran 90 style document below ! ! attributes: ! language: Fortran 90 and/or above ! machine: ! !$$$ end subprogram documentation block ! module interface: use kinds, only: i_kind,r_kind use mpeu_util, only: tell,perr,die use mpeu_util, only: assert_ use mpeu_util, only: stdout implicit none private ! except public :: cvgridLookup_sdget ! get subdomain parameters public :: cvgridLookup_islocal ! a lat-lon pair is "local" on a PE. public :: cvgridLookup_isluse ! a lat-lon pair is "luse" on a PE public :: cvgridLookup_getiw ! get interpolation operator (index,weight) interface cvgridLookup_sdget ; module procedure sdget_ ; end interface interface cvgridLookup_islocal; module procedure islocal_; end interface ! ... = []_local(elat,elon,myPE) interface cvgridLookup_isluse ; module procedure isluse_ ; end interface ! ... = []_luse(elat,elon,myPE) interface cvgridLookup_getiw; module procedure & get_Dij_ , & ! call []_getiw(elat,elon,ij(1:4),wij(1:4) [, debugging arguments]) get_Dijk_, & ! call []_getiw(elat,elon,dlev,ijk(1:8),wijk(1:8) [, debugging arguments]) get_Hij_ , & ! call []_getiw(elat,elon,ij(1:4),wij(1:4)) get_Hijk_ ! call []_getiw(elat,elon,dlev,ijk(1:8),wijk(1:8)) end interface !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ character(len=*),parameter :: myname='m_cvgridLookup' logical,save:: isubdom_configured_=.false. integer(i_kind),allocatable,dimension(:),save:: ilatbox_lbound, ilatbox_ubound integer(i_kind),allocatable,dimension(:),save:: ilonbox_lbound, ilonbox_ubound real(r_kind),save:: sdlat_ref_ , sdlon_ref_ ! sub-domain reference in degree lat-lon real(r_kind),save:: sdlat_lbnd_, sdlon_lbnd_ ! sub-domain lower-left corner in degree lat-lon real(r_kind),save:: sdlat_ubnd_, sdlon_ubnd_ ! sub-domain upper-right corner in degree lat-lon real(r_kind),parameter:: DEG_IN_RAD = 4._r_kind*atan(1._r_kind)/180._r_kind real(r_kind),parameter:: RAD_IN_DEG = 1._r_kind/DEG_IN_RAD integer(i_kind),parameter:: iHalo = 1 #include "mytrace.H" #include "myassert.H" contains !---- local routine: configure sub-domains of the control-vector grid. subroutine isubdom_config_() use gridmod, only: nlon, rlons, jlon1, jstart !, periodic_s use gridmod, only: nlat, rlats, ilat1, istart use gridmod, only: regional use mpimod , only: nPEs => npe ! number of PEs in GSI_comm_world use mpimod , only: myPE ! my rank in GSI_comm_world implicit none character(len=*),parameter:: myname_=myname//"::isubdom_config_" integer(i_kind):: ipe,ip1 integer(i_kind):: ilon,ilat !integer(i_kind):: mlat_lbound, mlat_ubound ! lower/upper lat. boundaries of the local subdomain !integer(i_kind):: mlon_lbound, mlon_ubound ! lower/upper lon. boundaries of the local subdomain real(r_kind):: r,rad_to_deg rad_to_deg(r) = r*RAD_IN_DEG ! The lower and the upper grid index boundaries of the subdomains of ! iHalo==1, are defined in the same way as they are in obs_para(), ! except that the array index are confined in 0:nPEs-1, instead of 1:nPEs. allocate(ilatbox_lbound(0:nPEs-1),ilatbox_ubound(0:nPEs-1)) allocate(ilonbox_lbound(0:nPEs-1),ilonbox_ubound(0:nPEs-1)) ! It was assumed in this isubdom_ module, the halo grid is +/-1 grid ! line. The current implementation let a single iHalo be a configurable ! parameter for both latitude and longitude grids. However, no test has ! been done for iHalo>1, which would need some changes to other part of ! the GSI to run. do ipe=0,nPEs-1 ! However, the periodicity variable periodic_s(:) are not used here, ! since its true definition is unknown. For example, they are undefined ! in at least one case of global periodic longitude grid. Their ! definitions for some regional/non-periodic longitude grids (expected ! to be .false.?) are also suspecious. ip1=ipe+1 ! latitude grid box index range of the local subdomain ilatbox_lbound(ipe) = istart(ip1)-iHalo ! lower grid box boundary ilatbox_ubound(ipe) = istart(ip1)+ilat1(ip1)-1 ilatbox_ubound(ipe) = ilatbox_ubound(ipe)-1 +iHalo ! upper grid box boundary without halo ! longitude grid box index range of the local subdomain ilonbox_lbound(ipe) = jstart(ip1)-iHalo ! lower grid box boundary ilonbox_ubound(ipe) = jstart(ip1)+jlon1(ip1)-1 ilonbox_ubound(ipe) = ilonbox_ubound(ipe)-1 +iHalo ! upper grid box boundary without halo enddo ip1=myPE+1 sdlat_ref_ = rad_to_deg(rlats(istart(ip1))) sdlon_ref_ = rad_to_deg(rlons(jstart(ip1))) ! Map local grid box range into local grid point range ! First, for the local latitude grid, the lower-bound and the upper-bound ilat= ilatbox_lbound(myPE) ! the lower grid point is the lower grid box ilat= min(max(1,ilat),nlat) ! map it to a valid grid index sdlat_lbnd_= rad_to_deg(rlats(ilat)) ! get the grid value ilat= ilatbox_ubound(myPE)+1 ! the upper grid point is the outer edge of the upper grid box ilat= min(max(1,ilat),nlat) ! map it to a valid grid index sdlat_ubnd_= rad_to_deg(rlats(ilat)) ! get the grid value ! Then, for the local longitude grid, the lower-bound and the upper-bound ilon=ilonbox_lbound(myPE) ! the lower grid point is the lower grid box if(regional) then ilon=min(max(1,ilon),nlon) ! map it to a valid non-periodic grid index else if(ilon<1) ilon= nlon+ilon ! map it to a valid periodic grid index endif sdlon_lbnd_= rad_to_deg(rlons(ilon)) ! get the grid value ilon=ilonbox_ubound(myPE)+1 ! the upper grid point is the outer edge of the upper grid box if(regional) then ilon=min(max(1,ilon),nlon) ! map it to a valid non-periodic grid index else if(ilon>nlon) ilon= ilon-nlon ! map it to a valid periodic grid index endif sdlon_ubnd_= rad_to_deg(rlons(ilon)) ! get the grid value isubdom_configured_ = .true. end subroutine isubdom_config_ !---- local routine: get the grid indices on the cv-grid. subroutine isubdom_index_(elat,elon,ilat,ilon) use gridmod, only: nlat,rlats use gridmod, only: nlon,rlons use kinds, only: i_kind, r_kind use constants, only: deg2rad implicit none real(kind=r_kind),intent(in ):: elat,elon integer(kind=i_kind),intent(out):: ilat,ilon real(kind=r_kind):: dlat,dlon dlat=elat*deg2rad dlon=elon*deg2rad call grdcrd1(dlat,rlats,nlat,1) ! returns dlat:=. call grdcrd1(dlon,rlons,nlon,1) ! returns dlon:=. ilat=dlat ! truncate out to only ilon=dlon ! truncate out to only ! For non-periodic (regional) grid, given a precondition of ! elon in [rlons(1),rlons(nlon)), the returned value of index, ! floor(dlon), is in [1,nlon-1]. ! ! For periodic (global) grid, given a precondition of elon is in ! [rlons(1),rlons(1)+360.), assuming ascending, or rlons(1)=1) ASSERT(ilat<=nlat) ASSERT(ilon>=0) ASSERT(ilon<=nlon) end subroutine isubdom_index_ !---- get subdomain boundaries subroutine sdget_( sdlat_ref, sdlat_lbnd, sdlat_ubnd, & sdlon_ref, sdlon_lbnd, sdlon_ubnd ) use kinds, only: r_kind implicit none ! subdomain grid reference, lower-bound, and upper-bound values real(r_kind),optional,intent(out):: sdlat_ref, sdlat_lbnd, sdlat_ubnd ! degree latitude real(r_kind),optional,intent(out):: sdlon_ref, sdlon_lbnd, sdlon_ubnd ! degree longitude if(.not.isubdom_configured_) call isubdom_config_() if(present(sdlat_ref )) sdlat_ref =sdlat_ref_ if(present(sdlat_lbnd)) sdlat_lbnd=sdlat_lbnd_ if(present(sdlat_ubnd)) sdlat_ubnd=sdlat_ubnd_ if(present(sdlon_ref )) sdlon_ref =sdlon_ref_ if(present(sdlon_lbnd)) sdlon_lbnd=sdlon_lbnd_ if(present(sdlon_ubnd)) sdlon_ubnd=sdlon_ubnd_ end subroutine sdget_ !----- external routine: determine if (elat,elon) is local on the cv-grid function islocal_(elat,elon,iPE) use gridmod, only: nlat use gridmod, only: nlon !,periodic_s use gridmod, only: regional implicit none logical :: islocal_ real (r_kind), intent(in):: elat,elon integer(i_kind), intent(in):: iPE integer(kind=i_kind):: ilat,ilon _ENTRY_('isluse_') if(.not.isubdom_configured_) call isubdom_config_() call isubdom_index_(elat,elon,ilat,ilon) ASSERT( allocated(ilonbox_lbound )) ASSERT(iPE>=lbound(ilonbox_lbound,1)) ASSERT(iPE<=ubound(ilonbox_lbound,1)) ! The conditions defining "local" were original kept the same as they ! were in obs_para(). However, the way it uses periodic_s(:) was ! questionable, either for a global periodic grid or for a regional grid. ! ! The current implementation avoids the use of periodic_s(:), and uses ! configuration variable _regional_ instead. ilat=min(max(1,ilat),nlat) ! ilat==nlat should not happen. See grdcrd() through isubdom_index_() islocal_ = ilat>=ilatbox_lbound(iPE) .and. ilat<=ilatbox_ubound(iPE) if(islocal_) then ilon=min(max(0,ilon),nlon) ! ilon==0 should not happen. See grdcrd() through isubdom_index_() ! ilon==nlon should happend only for periodic grid islocal_ = ilon>=ilonbox_lbound(iPE) .and. ilon<=ilonbox_ubound(iPE) if(.not.islocal_) then if(.not.regional) then ! periodicity applies islocal_ = ilon-nlon >= ilonbox_lbound(iPE) .or. & ! left to grid box 1 ilon+nlon <= ilonbox_ubound(iPE) ! right to grid box nlon-1 endif endif !-- This is the original algorithm as it is defined in obs_para() --- !islocal_ = ( ilon>=ilon_west(iPE) .and. ilon<=ilon_east(iPE) ) .or. & ! ( ilon==0 .and. ilon_east(iPE)>=nlon ) .or. & ! ( ilon==nlon .and. ilon_west(iPE)<=1 ) .or. ilon_periodic(iPE) endif _EXIT_('isluse_') end function islocal_ !----- external routine: determine if (elat,elon) is luse local on the cv-grid function isluse_(elat,elon,mype) result(luse) ! a wrapper to use the subroutine as a function. implicit none logical:: luse real (r_kind), intent(in ):: elat,elon integer(i_kind), intent(in ):: myPE call setluse_(luse,elat,elon,mype) end function isluse_ subroutine setluse_(luse,elat,elon,mype) implicit none logical , intent(out):: luse real (r_kind), intent(in ):: elat,elon integer(i_kind), intent(in ):: myPE integer(kind=i_kind):: ipe _ENTRY_('setluse_') ! Similar to what it was in obs_para() (and almost the same lately), for ! a given location (elat,elon), luse is set to .true., if and only if, ! islocal(myPE) .and. .not.any(islocal_(myPE-1:0:-1)) luse=islocal_(elat,elon,mype) ! It would not be here, if it were not even islocal() on the current PE. ! So does this assertion. ASSERT(luse) do ipe=mype-1,0,-1 ! On any PE of lower rank (iPE = mype-1, mype-2, ..., 0), if this ob ! islocal(), then it is .not. luse. if(islocal_(elat,elon,iPE)) then luse=.false. exit endif enddo _EXIT_('setluse_') end subroutine setluse_ subroutine get_Hij_(elat,elon,ij,wij) use gridmod, only: nlat,rlats use gridmod, only: nlon,rlons use gridmod, only: get_ij use constants, only: deg2rad use mpimod, only: myPE implicit none real(r_kind),intent(in):: elat,elon integer(i_kind),dimension(:),intent(inout):: ij real (r_kind),dimension(:),intent(inout):: wij real (r_kind):: dlat,dlon _ENTRY_('get_Hij_') dlat=elat*deg2rad dlon=elon*deg2rad call grdcrd1(dlat,rlats,nlat,1) ! returns dlat:=. call grdcrd1(dlon,rlons,nlon,1) ! returns dlon:=. call get_ij(myPE+1,dlat,dlon,ij(1:4),wij(1:4)) _EXIT_('get_Hij_') end subroutine get_Hij_ subroutine get_Hijk_(elat,elon,dlev,ijk,wijk) use gridmod, only: nlat,rlats use gridmod, only: nlon,rlons use gridmod, only: get_ijk use constants, only: deg2rad use mpimod, only: myPE implicit none real (r_kind),intent(in):: elat,elon,dlev integer(i_kind),dimension(:),intent(inout):: ijk real (r_kind),dimension(:),intent(inout):: wijk real (r_kind):: dlat,dlon _ENTRY_('get_Hijk_') dlat=elat*deg2rad dlon=elon*deg2rad call grdcrd1(dlat,rlats,nlat,1) ! returns dlat:=. call grdcrd1(dlon,rlons,nlon,1) ! returns dlon:=. call get_ijk(myPE+1,dlat,dlon,dlev,ijk(1:8),wijk(1:8)) _EXIT_('get_Hijk_') end subroutine get_Hijk_ subroutine get_Dij_(elat,elon,dlat,dlon,ij,wij,jtype,asis) use gridmod, only: nlat,rlats use gridmod, only: nlon,rlons use gridmod, only: get_ij use constants, only: deg2rad use mpimod, only: myPE implicit none real(r_kind),intent(in):: elat,elon real(r_kind),intent(in):: dlat,dlon integer(i_kind),dimension(:),intent(inout):: ij real (r_kind),dimension(:),intent(inout):: wij integer(i_kind),intent(in):: jtype logical,optional,intent(in):: asis real (r_kind):: dlat_,dlon_ integer(i_kind),dimension(4):: ij_ real (r_kind),dimension(4):: wij_ ! integer(i_kind):: ilat,ilon ! real (r_kind):: wlat,wlon logical:: asis_ _ENTRY_('get_Dij_') asis_=.false. if(present(asis)) asis_=asis if(asis_) then _EXIT_('get_Dij_') return endif dlat_=elat*deg2rad dlon_=elon*deg2rad call grdcrd1(dlat_,rlats,nlat,1) ! returns dlat:=. call grdcrd1(dlon_,rlons,nlon,1) ! returns dlon:=. ij_(1:4)= ij(1:4) wij_(1:4)=wij_(1:4) call get_ij(myPE+1,dlat_,dlon_,ij(1:4),wij(1:4)) ! call get_ij(myPE+1,dlat_,dlon_,ij(1:4),wij(1:4), & ! jjlat=ilat,jjlon=ilon,wwlat=wlat,wwlon=wlon,verbose=.true.) if(.not.all(ij_(1:4)==ij(1:4))) then call tell('get_Dij_','all( ij_(1:4)== ij(1:4)) =',all( ij_(1:4)== ij(1:4))) call tell('get_Dij_','all(wij_(1:4)==wij(1:4)) =',all(wij_(1:4)==wij(1:4))) call tell('get_Dij_',' jtype =',jtype) call tell('get_Dij_',' elat =',elat ) call tell('get_Dij_',' dlat =',dlat_) call tell('get_Dij_',' dlat0=',dlat ) call tell('get_Dij_',' dlat-dlat0=',dlat_-dlat) call tell('get_Dij_',' elon =',elon ) call tell('get_Dij_',' dlon =',dlon_) call tell('get_Dij_',' dlon0=',dlon ) call tell('get_Dij_',' dlon-dlon0=',dlon_-dlon) ! call tell('get_Dij_',' ilat =',ilat) ! call tell('get_Dij_',' ilon =',ilon) ! call tell('get_Dij_',' wlat =',wlat) ! call tell('get_Dij_',' wlon =',wlon) write(stdout,'(a,8i8 )') 'get_Dij(): ij (1:4) =', ij (1:4) write(stdout,'(a,8i8 )') 'get_Dij(): ij~(1:4) =', ij_(1:4) write(stdout,'(a,8f8.4)') 'get_Dij(): wij (1:4) =',wij (1:4) write(stdout,'(a,8f8.4)') 'get_Dij(): wij~(1:4) =',wij_(1:4) call die('get_Dij_') endif ! ASSERT(all( ij_(1:4)== ij(1:4))) ! ASSERT(all(wij_(1:4)==wij(1:4))) _EXIT_('get_Dij_') end subroutine get_Dij_ subroutine get_Dijk_(elat,elon,dlat,dlon,dlev,ijk,wijk,jtype,asis) use gridmod, only: nlat,rlats use gridmod, only: nlon,rlons use gridmod, only: get_ijk use constants, only: deg2rad use mpimod, only: myPE implicit none real (r_kind),intent(in):: elat,elon,dlev real (r_kind),intent(in):: dlat,dlon integer(i_kind),dimension(:),intent(inout):: ijk real (r_kind),dimension(:),intent(inout):: wijk integer(i_kind) ,intent(in):: jtype logical,optional,intent(in):: asis integer(i_kind),dimension(8):: ijk_ real (r_kind),dimension(8):: wijk_ real (r_kind):: dlat_,dlon_ logical:: asis_ _ENTRY_('get_Dijk_') asis_=.false. if(present(asis)) asis_=asis if(asis_) then _EXIT_('get_Dijk_') return endif dlat_=elat*deg2rad dlon_=elon*deg2rad call grdcrd1(dlat_,rlats,nlat,1) ! returns dlat:=. call grdcrd1(dlon_,rlons,nlon,1) ! returns dlon:=. ijk_(1:8)= ijk(1:8) wijk_(1:8)=wijk(1:8) call get_ijk(myPE+1,dlat_,dlon_,dlev,ijk(1:8),wijk(1:8)) if(.not.all(ijk_(1:8)==ijk(1:8))) then call tell('get_Dijk_','all( ijk_(1:8)== ijk(1:8)) =',all( ijk_(1:8)== ijk(1:8))) call tell('get_Dijk_','all(wijk_(1:8)==wijk(1:8)) =',all(wijk_(1:8)==wijk(1:8))) call tell('get_Dijk_',' jtype =',jtype) call tell('get_DijK_',' elat =',elat ) call tell('get_Dijk_',' dlat =',dlat_) call tell('get_Dijk_',' dlat0=',dlat ) call tell('get_Dijk_',' dlat-dlat0=',dlat_-dlat) call tell('get_DijK_',' elon =',elon ) call tell('get_Dijk_',' dlon =',dlon_) call tell('get_Dijk_',' dlon0=',dlon ) call tell('get_Dijk_',' dlon-dlon0=',dlon_-dlon) call tell('get_Dijk_',' dlev =',dlev ) write(stdout,'(a,8i8 )') 'get_Dijk(): ijk (1:8) =', ijk (1:8) write(stdout,'(a,8i8 )') 'get_Dijk(): ijk~(1:8) =', ijk_(1:8) write(stdout,'(a,8f8.4)') 'get_Dijk(): wijk (1:8) =',wijk (1:8) write(stdout,'(a,8f8.4)') 'get_Dijk(): wijk~(1:8) =',wijk_(1:8) call die('get_Dijk_') endif ! ASSERT(all( ijk_(1:8)/= 0)) ! ASSERT(all( ijk_(1:8)== ijk(1:8))) ! ASSERT(all(wijk_(1:8)==wijk(1:8))) _EXIT_('get_Dijk_') end subroutine get_Dijk_ end module m_cvgridLookup