module m_latlonRange
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    module m_latlonRange
!   prgmmr:      j guo <jguo@nasa.gov>
!      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
!.