module m_latlonRange !$$$ subprogram documentation block ! . . . . ! subprogram: module m_latlonRange ! prgmmr: j guo ! org: NASA/GSFC, Global Modeling and Assimilation Office, 610.3 ! date: 2016-02-09 ! ! abstract: lat-lon range information for subdomains and observation clusters ! ! program history log: ! 2016-02-09 j guo - initial code ! 2016-02-09 j guo - added this document block ! 2016-06-22 j guo - refined the definition for local subdomain ranges ! . added text-dump routines for diagnosis ! . added "pure" qualifier to some routines. ! ! 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_mpif, only: MPI_iKIND use mpeu_util, only: assert_ use mpeu_util, only: tell,perr,die implicit none private ! except public :: latlonRange ! data structure public :: latlonRange_reset ! the initializer or the finalizer public :: latlonRange_enclose ! updating the range by including a (elat,elon) public :: latlonRange_islocal ! falls (partially) in the "local" subdomain. public :: latlonRange_set ! set the range for a known lat-lon domain public :: latlonRange_overlaps ! []_overlaps(ar, br) is "(br) overlaps (ar)" public :: latlonRange_gatherWrite ! gather then root-write, for all-ranges public :: latlonRange_readBcast ! root-read then bcast, for all-ranges ! Interfaces for debugging purposes public :: latlonRange_alldump ! all-text-dump a local latlonRange array. interface latlonRange_alldump; module procedure alldump_; end interface public :: latlonRange_gatherdump ! gather local scalar latlonRanges, then root-text-dump all. interface latlonRange_gatherdump; module procedure & gatherdump_local_, & ! gather-dump local cvgrid lat-lon-range gatherdump_; end interface ! gather-dump user provided lat-lon-range type latlonRange private integer(kind=i_kind):: nuse=0 ! > 0; count of luse==.true., or ! 0; for no data (so no range defined), or ! -1; for a subdomain grid range. ! latitude range value defined within an absolute domain, ! [-90.,+90.] real(kind=r_kind):: alat_min=0. ! the lower bound of lat real(kind=r_kind):: alat_max=0. ! the upper bound of lat ! longitude range value is defined within a relative domain, ! [alon_ref-180.,alon_ref+180.] real(kind=r_kind):: alon_ref=0. ! a reference lon used for all data real(kind=r_kind):: alon_min=0. ! the lower bound of lon real(kind=r_kind):: alon_max=0. ! the upper bound of lon end type latlonRange interface latlonRange_reset; module procedure & reset_, & ! rank-0 reset_r1_; end interface ! rank-1 interface latlonRange_enclose ; module procedure enclose_ ; end interface interface latlonRange_islocal ; module procedure islocal_ ; end interface interface latlonRange_set ; module procedure set_ ; end interface interface latlonRange_overlaps ; module procedure overlaps_; end interface interface latlonRange_gatherWrite; module procedure gatherWrite_; end interface interface latlonRange_readBcast ; module procedure readBcast_ ; end interface ! Usecase I: Creat a file of all-ranges, for meta info of the all-obsdiags outputs. ! ! use m_latlonRange, only: latlonRange ! use m_latlonRange, only: latlonRange_reset ! use m_latlonRange, only: latlonRange_enclose ! use m_latlonRange, only: latlonRange_gatherWrite() ! type(latlonRange):: oRange ! ! call latlonRange_reset(oRange) ! do .. in an obstypes node entries ! if(obsNode%isluse()) & ! call latlonRange_enclose(oRange, obsNode%elat,obsNode%elon) ! enddo ! ! call latlonRange_gatherWrite(oRange,'obsdiags.01',root,comm) ! call latlonRange_reset(oRange) ! ! Usecase II: Read in a file of all-ranges, to identify the locality before the ! read of any obsdiags file. ! ! type(latlonRange):: myRange ! type(latlonRange),allocatable,dimension(:):: allRanges ! ! allocatge(allRanges(0:mPEs-1) ! call latlonRange_reset(allRanges(:)) ! call latlonRange_readBcast('obsdiags.01.headers',allRanges(:),comm,root) ! call latlonRange_set(myRange, ...) ! ! do iPE=lPE,uPE ! if(latlonRange_overlaps(myRange,allRanges(iPE)) then ! call read_(...) ! endif ! enddo ! ! call latlonRange_reset(allRanges) ! call latlonRange_reset( myRange ) ! deallocate(allRanges) real(r_kind),parameter:: alat_NP=+90._r_kind real(r_kind),parameter:: alat_SP=-90._r_kind 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 real(r_kind),parameter:: R360DEG =360._r_kind integer(i_kind),parameter:: ROOT=0 integer(kind=MPI_iKIND),save:: mpi_type_=0 logical,save:: mpi_type_defined_=.false. type(latlonRange),save:: localRange_ = latlonRange() logical,save:: localRange_defined_ =.false. logical,parameter:: alwaysLocal_ =.true. #include "myassert.H" #define _TIMER_ON_ #ifdef _TIMER_ON_ #undef _TIMER_ON_ #undef _TIMER_OFF_ #undef _TIMER_USE_ #define _TIMER_ON_(id) call timer_ini(id) #define _TIMER_OFF_(id) call timer_fnl(id) #define _TIMER_USE_ use timermod, only: timer_ini,timer_fnl #else #define _TIMER_ON_(id) #define _TIMER_OFF_(id) #define _TIMER_USE_ #endif !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ character(len=*),parameter :: myname='m_latlonRange' contains pure function lonDiff_(clon,glon) result(dlon) ! lonDiff_ = dlon := (glon.minus.clon +180.) .modulo. 360. -180. ! Note: it is NOT clon.minus.glon, BUT glon.minus.clon. implicit none real(kind=r_kind):: dlon ! is defined in range [-180.,180.]. ! The actual result for an argument at an end point could ! be either -180 or +180, depending on the value of glon ! relative to clon. real(kind=r_kind),intent(in):: clon,glon real(kind=kind(dlon)):: rad_clon,cos_clon,sin_clon real(kind=kind(dlon)):: rad_glon,cos_glon,sin_glon !print'(a,2f10.3)','DEG_in_RAD=',DEG_in_RAD,3.1415926/180. !print'(a,2f10.3)','RAD_in_DEG=',RAD_in_DEG,180.00/3.1415926 rad_clon = clon *DEG_in_RAD cos_clon = cos(rad_clon) sin_clon = sin(rad_clon) rad_glon = glon *DEG_in_RAD cos_glon = cos(rad_glon) sin_glon = sin(rad_glon) ! dlon = glon .minus. clon, in a trignomitry based algorithm. dlon = atan2( cos_clon*sin_glon - sin_clon*cos_glon, & ! sin(g-c) = sin(g)*cos(c)-cos(g)*sin(c) cos_clon*cos_glon + sin_clon*sin_glon) ! cos(g-c) = cos(g)*cos(c)+sin(g)*sin(c) dlon = dlon *RAD_in_DEG ! d := "g-c" = atan2(sin(g-c),cos(g-c)) end function lonDiff_ pure subroutine reset_(llrange) ! set llrange to null. implicit none type(latlonRange),intent(out):: llrange llrange=latlonRange() end subroutine reset_ pure subroutine reset_r1_(llranges) implicit none type(latlonRange),dimension(:),intent(out):: llranges llranges(:)=latlonRange() end subroutine reset_r1_ pure subroutine set_(llrange,alat1,alon1,alat2,alon2) ! set llrange for the local lat-lon subdomain implicit none type(latlonRange),intent(out):: llrange real(kind=r_kind),intent( in):: alat1,alon1 ! corner #1 real(kind=r_kind),intent( in):: alat2,alon2 ! corner #2 llrange=latlonRange() ! reinitialize to zero call enclose_(llrange,alat1,alon1) call enclose_(llrange,alat2,alon2) llrange%nuse=-1 ! it suggests this range is for a grid. end subroutine set_ pure subroutine enclose_(llrange,alat,alon) implicit none type(latlonRange),intent(inout):: llrange real(kind=r_kind),intent(in ):: alat,alon real(kind=kind(alat)):: alat_,alon_ real(kind=kind(alat)),parameter:: alon_HalfInterval=360._r_kind*epsilon(1._r_kind) if(llrange%nuse==0) then llrange%alat_min=alat llrange%alat_max=alat llrange%alon_ref=alon llrange%alon_min=alon -alon_HalfInterval llrange%alon_max=alon +alon_HalfInterval else alat_=max(alat_SP,min(alat,alat_NP)) llrange%alat_min=min(llrange%alat_min,alat_) llrange%alat_max=max(llrange%alat_max,alat_) alon_ = llrange%alon_ref + lonDiff_(llrange%alon_ref,alon) llrange%alon_min=min(llrange%alon_min,alon_-alon_HalfInterval) llrange%alon_max=max(llrange%alon_max,alon_+alon_HalfInterval) endif llrange%nuse=llrange%nuse+1 end subroutine enclose_ subroutine localRange_config_() !-- Config localRange_ for its lat-lon range. use m_cvgridLookup, only: cvgridLookup_sdget implicit none real(r_kind):: sdlat_ref ,sdlon_ref ! a reference lat-lon grid point real(r_kind):: sdlat_lbnd,sdlon_lbnd ! lower-left corner real(r_kind):: sdlat_ubnd,sdlon_ubnd ! upper-right corner call cvgridLookup_sdget( sdlat_ref =sdlat_ref , sdlon_ref =sdlon_ref , & sdlat_lbnd=sdlat_lbnd, sdlon_lbnd=sdlon_lbnd, & sdlat_ubnd=sdlat_ubnd, sdlon_ubnd=sdlon_ubnd ) localRange_=latlonRange() ! reinitialize call enclose_(localRange_,sdlat_ref ,sdlon_ref ) ! include a reference grid point call enclose_(localRange_,sdlat_lbnd,sdlon_lbnd) ! include the lower-left corner call enclose_(localRange_,sdlat_ubnd,sdlon_ubnd) ! include the upper-right corner ! This algorithm does not assume sdlat_lbnd < sdlat_ubnd or ! sdlon_lbnd < sdlon_ubnd. localRange_defined_=.true. end subroutine localRange_config_ function islocal_(r) !-- is the given latlonRange::r local, against *localRange_*? implicit none logical:: islocal_ type(latlonRange),intent(in):: r if(.not.localRange_defined_) call localRange_config_() islocal_=overlaps_(localRange_,r) !-- Or for debugging purposes (see module parameter alwaysLocal_ defined on !-- the top), let ! islocal_ = alwaysLocal_.or.overlaps_(localRange_,r) end function islocal_ function overlaps_(aRange,bRange) ! overlaps_(ar,br) := br .overlaps. ar implicit none logical:: overlaps_ type(latlonRange),intent(in):: aRange type(latlonRange),intent(in):: bRange ! An overlaping status of a pair of ranges, requires neither to be null overlaps_ = aRange%nuse /= 0 .and. bRange%nuse /= 0 ! In some cases, both boundaries of the larger range may be outside ! the smaler range. So we need to swap the two ranges to find the ! right result. if(overlaps_) overlaps_ = lat_overlaps_(aRange,bRange) .or. & lat_overlaps_(bRange,aRange) ! an overlaping status is both overlaping-in-lat and overlaping-in-lon. if(overlaps_) overlaps_ = lon_overlaps_(aRange,bRange) .or. & lon_overlaps_(bRange,aRange) contains pure function lat_overlaps_(ar,br) result(overlaps_) ! at least one of boundary points of br, falls in the range of ar implicit none logical:: overlaps_ type(latlonRange),intent(in):: ar,br real(kind=r_kind):: alat_ alat_ = br%alat_min overlaps_ = alat_>=ar%alat_min .and. & alat_<=ar%alat_max if(.not.overlaps_) then alat_ = br%alat_max overlaps_ = alat_>=ar%alat_min .and. & alat_<=ar%alat_max endif end function lat_overlaps_ pure function lon_overlaps_(ar,br) result(overlaps_) ! at least one of boundary points of br, falls in the range of ar implicit none logical:: overlaps_ type(latlonRange),intent(in):: ar,br real(kind=r_kind):: aref_,alon_ aref_ = ar%alon_ref alon_ = aref_ + lonDiff_(aref_,br%alon_min) overlaps_ = alon_ >= ar%alon_min .and. & alon_ <= ar%alon_max if(.not.overlaps_) then alon_ = aref_ + lonDiff_(aref_,br%alon_max) overlaps_ = alon_ >= ar%alon_min .and. & alon_ <= ar%alon_max endif end function lon_overlaps_ end function overlaps_ subroutine gatherWrite_(llrange,hdfile,root,comm) ! gather all llranges from all PEs of this comm, then write them all to ! a "header" file (hdfile="obsdiags.##.headers") by the root PE. use mpeu_mpif, only: MPI_type use gsi_unformatted, only: unformatted_open use mpimod, only: nPEs => nPE use mpimod, only: myPE _TIMER_USE_ implicit none type(latlonRange ),intent(in):: llrange ! a laglonRange to write character(len=* ),intent(in):: hdfile ! a filename to write to integer(kind=i_kind),intent(in):: root integer(kind=i_kind),intent(in):: comm character(len=*),parameter:: myname_=myname//"::gatherWrite_" integer(kind=i_kind):: ier,lu integer(kind=i_kind):: lsize,isize,rsize,mtype integer(kind=i_kind):: irec,nrec integer(kind=i_kind), dimension(1 ):: isend integer(kind=i_kind),allocatable,dimension(:,:):: irecv real (kind=r_kind), dimension(5 ):: rsend real (kind=r_kind),allocatable,dimension(:,:):: rrecv _TIMER_ON_(myname_) lsize=0 if(myPE==root) lsize=nPEs isize =size(isend) rsize =size(rsend) allocate( irecv(isize,0:lsize-1), & rrecv(rsize,0:lsize-1) ) isend(1)=llrange%nuse mtype =MPI_type(isend) call MPI_gather(isend,isize,mtype, irecv,isize,mtype, root,comm,ier) ASSERT(ier==0) rsend(:) = (/ llrange%alat_min, & llrange%alat_max, & llrange%alon_ref, & llrange%alon_min, & llrange%alon_max /) mtype =MPI_type(rsend) call MPI_gather(rsend,rsize,mtype, rrecv,rsize,mtype, root,comm,ier) ASSERT(ier==0) if(myPE==root) then call unformatted_open(unit=lu, & file=trim(hdfile), & class='.obsdiags.', & action='write', & status='unknown', & newunit=.true., & ! with newunit=.true., unit returns a value assigned by Fortran. iostat=ier,silent=.true.) if(ier/=0) then call perr(myname_,'unformatted_open(), file =',trim(hdfile)) call perr(myname_,' newunit =',lu) call perr(myname_,' iostat =',ier) call die(myname_) endif nrec=size(irecv,2) write(lu) nrec do irec=0,nrec-1 write(lu) irec, irecv(1:isize,irec),rrecv(1:rsize,irec) enddo close(lu) endif deallocate(irecv,rrecv) _TIMER_OFF_(myname_) end subroutine gatherWrite_ subroutine readBcast_(hdfile,allRanges,root,comm) use gsi_unformatted, only: unformatted_open use mpeu_mpif, only: MPI_type use mpimod, only: myPE _TIMER_USE_ implicit none character(len=* ),intent(in):: hdfile ! input file type(latlonRange),dimension(0:),intent(out):: allranges ! data received by all PEs. integer(kind=i_kind),intent(in):: root ! the root PE does the read(). integer(kind=i_kind),intent(in):: comm ! the communicator of my "world". character(len=*),parameter:: myname_=myname//"::readBcast_" integer(kind=i_kind):: ier,lu,irec,jrec,nrec integer(kind=i_kind),allocatable,dimension(:,:):: ibufr real (kind=r_kind),allocatable,dimension(:,:):: rbufr _TIMER_ON_(myname_) nrec=0 if(myPE==root) then call unformatted_open(unit=lu, & file=trim(hdfile), & class='.obsdiags.', & action='read', & status='old', & newunit=.true., & ! with newunit=.true., unit returns a value assigned by Fortran. iostat=ier,silent=.true.) if(ier/=0) then call perr(myname_,'unformatted_open(), file =',trim(hdfile)) call perr(myname_,' newunit =',lu) call perr(myname_,' iostat =',ier) call die(myname_) endif read(lu) nrec endif call MPI_bcast(nrec,1,mpi_type(nrec),root,comm,ier) ASSERT(ier==0) allocate( ibufr(1,0:nrec-1), & rbufr(5,0:nrec-1) ) if(myPE==root) then do irec=0,nrec-1 read(lu) jrec, ibufr(1,irec),rbufr(1:5,irec) ASSERT(irec==jrec) enddo close(lu) endif call MPI_bcast(ibufr,size(ibufr),mpi_type(ibufr),root,comm,ier) ASSERT(ier==0) call MPI_bcast(rbufr,size(rbufr),mpi_type(rbufr),root,comm,ier) ASSERT(ier==0) ASSERT(nrec==size(allRanges)) do irec=0,nrec-1 allRanges(irec)=latlonRange( nuse=ibufr(1,irec), & alat_min=rbufr(1,irec), & alat_max=rbufr(2,irec), & alon_ref=rbufr(3,irec), & alon_min=rbufr(4,irec), & alon_max=rbufr(5,irec) ) enddo deallocate(ibufr,rbufr) _TIMER_OFF_(myname_) end subroutine readBcast_ subroutine alldump_(allRanges,varname) !-- simply text-dump a local allRanges to stdout use mpeu_util, only: stdout use mpeu_util, only: stdout_lead use mpeu_util, only: stdout_open use mpeu_util, only: stdout_close _TIMER_USE_ implicit none type(latlonRange), dimension(0:), intent(in):: allRanges character(len=* ), intent(in):: varname character(len=*),parameter:: myname_=myname//"::alldump_" integer(i_kind):: irec character(len=:), allocatable:: varlead_ _TIMER_ON_(myname_) varlead_=stdout_lead(varname) write (stdout,'(a,i4,l2,i8,5f12.4)') varlead_, -1,islocal_(localRange_ ),localRange_ write (stdout,'(a,1x,a,i8)') varlead_,'NPES =',size(allranges) do irec=0,size(allRanges)-1 write(stdout,'(a,i4,l2,i8,5f12.4)') varlead_, irec,islocal_( allRanges(irec)),allRanges(irec) enddo write (stdout,'(a,1x,a)') varlead_,'========' _TIMER_OFF_(myname_) end subroutine alldump_ subroutine gatherdump_local_(varname,root,comm) !-- gahter-dump the internal lat-lon-range of the cvgrid, localRange_ _TIMER_USE_ implicit none character(len=* ),intent(in):: varname ! identity of the range integer(kind=i_kind),intent(in):: root integer(kind=i_kind),intent(in):: comm character(len=*),parameter:: myname_=myname//"::gatherdump_local_" _TIMER_ON_(myname_) if(.not.localRange_defined_) call localRange_config_() call gatherdump_(localRange_,varname,root,comm) _TIMER_OFF_(myname_) end subroutine gatherdump_local_ subroutine gatherdump_(llrange,varname,root,comm) !-- gather all llranges from all PEs of this comm, then text dump the gathered !-- array of llranges to stdout. use mpeu_mpif, only: MPI_type use mpeu_util, only: stdout use mpeu_util, only: stdout_lead use mpimod, only: nPEs => nPE use mpimod, only: myPE _TIMER_USE_ implicit none type(latlonRange ),intent(in):: llrange ! a laglonRange to write character(len=* ),intent(in):: varname ! identity of the range integer(kind=i_kind),intent(in):: root integer(kind=i_kind),intent(in):: comm character(len=*),parameter:: myname_=myname//"::gatherdump_" integer(kind=i_kind):: ier integer(kind=i_kind):: lsize,isize,rsize,mtype integer(kind=i_kind):: irec,nrec integer(kind=i_kind), dimension(1 ):: isend integer(kind=i_kind),allocatable,dimension(:,:):: irecv real (kind=r_kind), dimension(5 ):: rsend real (kind=r_kind),allocatable,dimension(:,:):: rrecv character(len=:),allocatable:: varlead_ _TIMER_ON_(myname_) lsize=0 if(myPE==root) lsize=nPEs isize =size(isend) rsize =size(rsend) allocate( irecv(isize,0:lsize-1), & rrecv(rsize,0:lsize-1) ) isend(1)=llrange%nuse mtype =MPI_type(isend) call MPI_gather(isend,isize,mtype, irecv,isize,mtype, root,comm,ier) ASSERT(ier==0) rsend(:) = (/ llrange%alat_min, & llrange%alat_max, & llrange%alon_ref, & llrange%alon_min, & llrange%alon_max /) mtype =MPI_type(rsend) call MPI_gather(rsend,rsize,mtype, rrecv,rsize,mtype, root,comm,ier) ASSERT(ier==0) varlead_=stdout_lead(myname,varname) if(myPE==root) then nrec=size(irecv,2) write (stdout,'(a,1x,a,i8)') varlead_,'nrec =',nrec do irec=0,nrec-1 write(stdout,'(a,2i8,5f12.4)') varlead_,irec, irecv(1:isize,irec),rrecv(1:rsize,irec) enddo write (stdout,'(a,1x,a)') varlead_,'==========' endif deallocate(irecv,rrecv) _TIMER_OFF_(myname_) end subroutine gatherdump_ end module m_latlonRange !.