module general_tll2xy_mod
!$$$   module documentation block
!                .      .    .                                       .
! module:    generic_tll2xy_mod  generalized conversion from latlon to xy
!   prgmmr: parrish          org: np22                date: 2010-10-28
!
! abstract: This module contains generalized tll2xy and related routines to convert
!             earth lat lon coordinates to grid coordinates for an orthogonal
!             non-staggered grid for which the earth lat and lon is known for each grid point.
!            NOTE: all routines tested against same routines in gridmod, using region_lat, region_lon first,
!                then random x,y,ug,vg.  get exact match between gridmod routines and the corresponding
!                  general_ routines here.  See subroutine test1_general_ll2xy at end of this module.
!
! program history log:
!   2010-10-28  parrish - initial documentation
!   2013-01-23  parrish - modify so calls to ll2rpolar and rpolar2ll avoid type mismatch error
!                           when using WCOSS intel debug compile.
!
! subroutines included:
!   sub general_create_llxy_transform - initialize type(llxy_cons) for desired grid
!   sub general_tll2xy                - convert from lat-lon to grid coordinates for desired grid
!   sub general_txy2ll                - convert from grid coordinates to lat-lon for desired grid
!   sub general_rotate_wind_ll2xy     - rotate earth wind to grid wind for desired grid
!   sub general_rotate_wind_xy2ll     - rotate grid wind to earth wind for desired grid
!
! Variable Definitions:
!   def llxy_cons:                    - type variable which is used to make parameter set previously
!                                         hardwired in gridmod.f90 portable for arbitrary grids.
!
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$ end documentation block

   use kinds, only: r_kind,i_kind,i_byte

   implicit none

! set default to private
   private
! set subroutines to public
   public :: general_create_llxy_transform
   public :: general_tll2xy
   public :: general_txy2ll
   public :: general_rotate_wind_ll2xy
   public :: general_rotate_wind_xy2ll
! set passed variables to public
   public :: llxy_cons

   type llxy_cons
! The following is for the generalized transform
      real(r_kind) rlon_min_dd,rlon_max_dd,rlat_min_dd,rlat_max_dd
      real(r_kind) pihalf,sign_pole,rlambda0
      real(r_kind) atilde_x,btilde_x,atilde_y,btilde_y
      real(r_kind) btilde_xinv,btilde_yinv
      integer(i_kind) nlon,nlat,nxtilde,nytilde
      integer(i_kind),pointer::i0_tilde(:,:) => NULL()
      integer(i_kind),pointer::j0_tilde(:,:) => NULL()
      integer(i_byte),pointer::ip_tilde(:,:) => NULL()
      integer(i_byte),pointer::jp_tilde(:,:) => NULL()
      real(r_kind),pointer::xtilde0(:,:) => NULL()
      real(r_kind),pointer::ytilde0(:,:) => NULL()
      real(r_kind),pointer::cos_beta_ref(:,:) => NULL()
      real(r_kind),pointer::sin_beta_ref(:,:) => NULL()
      real(r_kind),pointer::region_lat(:,:) => NULL()
      real(r_kind),pointer::region_lon(:,:) => NULL()
      logical:: lallocated = .false.

   end type llxy_cons

!  other declarations  ...

   contains

 subroutine general_create_llxy_transform(region_lat,region_lon,nlat,nlon,gt)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    general_create_llxy_transform
!   prgmmr:  parrish
!
! abstract:  copy routines from gridmod.f90 to make a general purpose module which allows
!     conversion of earth lat lons to grid coordinates for any non-staggered rectangular 
!     orthogonal grid with known earth lat and lon of each grid point.
!     set up constants to allow conversion between earth lat lon and analysis grid units.
!     There is no need to specify details of the analysis grid projection.  All that is required
!     is the earth latitude and longitude in radians of each analysis grid point.
!
! program history log:
!   2010-10-28  parrish - initial documentation
!   2012-11-28  tong - added gt%lallocated=.true. after arrays of gt are allocated and 
!                      removed the duplicate allocation of gt%region_lat,gt%region_lon
!                      at the begining
!
!   input argument list:
!    glats       - earth latitudes of each grid point for desired grid.
!    glons       - earth longitudes of each grid point for desired grid.
!
!   output argument list:
!    gt          - variable of type llxy_cons, which contains all information needed for 
!                   conversion between earth lat lon and this grid grid units.
!
! attributes:
!   language: f90
!   machine:
!
!$$$ end documentation block

  use kinds, only: r_kind,i_kind
  use constants, only: zero,one,half,pi
  implicit none

  integer(i_kind),intent(in   ) :: nlat,nlon
  real(r_kind)   ,intent(in   ) :: region_lat(nlat,nlon),region_lon(nlat,nlon)
  type(llxy_cons),intent(inout) :: gt

  real(r_kind),parameter:: rbig =1.0e30_r_kind
  real(r_kind) xbar_min,xbar_max,ybar_min,ybar_max
  real(r_kind) clon,slon,r_of_lat,xbar,ybar
  integer(i_kind) i,j,istart0,iend,iinc,itemp,ilast,jlast
  real(r_kind) glats(nlon,nlat),glons(nlon,nlat)
  real(r_kind),allocatable:: clata(:,:),slata(:,:),clona(:,:),slona(:,:)
  real(r_kind) clat0,slat0,clon0,slon0
  real(r_kind) clat_m1,slat_m1,clon_m1,slon_m1
  real(r_kind) clat_p1,slat_p1,clon_p1,slon_p1
  real(r_kind) x,y,z,xt,yt,zt,xb,yb,zb
  real(r_kind) rlonb_m1,clonb_m1,slonb_m1
  real(r_kind) rlonb_p1,clonb_p1,slonb_p1
  real(r_kind) crot,srot

  do j=1,nlon
     do i=1,nlat
        glats(j,i)=region_lat(i,j)
        glons(j,i)=region_lon(i,j)
     end do
  end do
  gt%pihalf=half*pi
  gt%nlon=nlon
  gt%nlat=nlat
  gt%rlon_min_dd=one
  gt%rlat_min_dd=one
  gt%rlon_max_dd=nlon
  gt%rlat_max_dd=nlat

!  define xtilde, ytilde grid, transform

!      glons,glats are lons, lats of input grid points of dimension nlon,nlat
  call general_get_xytilde_domain(gt,gt%nlon,gt%nlat,glons,glats,gt%nxtilde,gt%nytilde, &
                   xbar_min,xbar_max,ybar_min,ybar_max)
  if(gt%lallocated) then
     deallocate(gt%i0_tilde,gt%j0_tilde,gt%ip_tilde,gt%jp_tilde,gt%xtilde0,gt%ytilde0)
     deallocate(gt%cos_beta_ref,gt%sin_beta_ref,gt%region_lat,gt%region_lon)
     gt%lallocated=.false.
  end if
  allocate(gt%i0_tilde(gt%nxtilde,gt%nytilde),gt%j0_tilde(gt%nxtilde,gt%nytilde))
  allocate(gt%ip_tilde(gt%nxtilde,gt%nytilde),gt%jp_tilde(gt%nxtilde,gt%nytilde))
  allocate(gt%xtilde0(gt%nlon,gt%nlat),gt%ytilde0(gt%nlon,gt%nlat))
  allocate(gt%cos_beta_ref(gt%nlon,gt%nlat),gt%sin_beta_ref(gt%nlon,gt%nlat))
  allocate(gt%region_lat(gt%nlat,gt%nlon),gt%region_lon(gt%nlat,gt%nlon))

  gt%lallocated=.true.

  do j=1,nlon
     do i=1,nlat
        gt%region_lat(i,j)=region_lat(i,j)
        gt%region_lon(i,j)=region_lon(i,j)
     end do
  end do

! define atilde_x, btilde_x, atilde_y, btilde_y

  gt%btilde_x   =(gt%nxtilde -one     )/(xbar_max-xbar_min)
  gt%btilde_xinv=(xbar_max-xbar_min)/(gt%nxtilde -one     )
  gt%atilde_x   =one-gt%btilde_x*xbar_min
  gt%btilde_y   =(gt%nytilde -one     )/(ybar_max-ybar_min)
  gt%btilde_yinv=(ybar_max-ybar_min)/(gt%nytilde -one     )
  gt%atilde_y   =one-gt%btilde_y*ybar_min

! define xtilde0,ytilde0
  do j=1,gt%nlat
     do i=1,gt%nlon
        r_of_lat=gt%pihalf+gt%sign_pole*glats(i,j)
        clon=cos(glons(i,j)+gt%rlambda0)
        slon=sin(glons(i,j)+gt%rlambda0)
        xbar=r_of_lat*clon
        ybar=r_of_lat*slon
        gt%xtilde0(i,j)=gt%atilde_x+gt%btilde_x*xbar
        gt%ytilde0(i,j)=gt%atilde_y+gt%btilde_y*ybar
     end do
  end do

!  now get i0_tilde, j0_tilde, ip_tilde,jp_tilde
  ilast=1 ; jlast=1
  istart0=gt%nxtilde
  iend=1
  iinc=-1
  do j=1,gt%nytilde
     itemp=istart0
     istart0=iend
     iend=itemp
     iinc=-iinc
     ybar=j
     do i=istart0,iend,iinc
        xbar=i
        call general_nearest_3(ilast,jlast,gt%i0_tilde(i,j),gt%j0_tilde(i,j), &
                       gt%ip_tilde(i,j),gt%jp_tilde(i,j),xbar,ybar,gt%nlon,gt%nlat,gt%xtilde0,gt%ytilde0)
     end do
  end do

!   new, more accurate and robust computation of cos_beta_ref and sin_beta_ref which is independent
!     of sign_pole and works for any orientation of grid on sphere (only restriction for now is that
!     x-y coordinate of analysis grid is right handed).
  allocate(clata(gt%nlon,gt%nlat),slata(gt%nlon,gt%nlat),clona(gt%nlon,gt%nlat),slona(gt%nlon,gt%nlat))
  do j=1,gt%nlat
     do i=1,gt%nlon
        clata(i,j)=cos(glats(i,j))
        slata(i,j)=sin(glats(i,j))
        clona(i,j)=cos(glons(i,j))
        slona(i,j)=sin(glons(i,j))
     end do
  end do
  do j=1,gt%nlat
     do i=2,gt%nlon-1

!     do all interior lon points to 2nd order accuracy

!   transform so pole is at rlat0,rlon0 and 0 meridian is tangent to earth latitude at rlat0,rlon0.

        clat0=clata(i,j) ; slat0=slata(i,j) ; clon0=clona(i,j) ; slon0=slona(i,j)

!    now obtain new coordinates for m1 and p1 points.

        clat_m1=clata(i-1,j) ; slat_m1=slata(i-1,j) ; clon_m1=clona(i-1,j) ; slon_m1=slona(i-1,j)
        clat_p1=clata(i+1,j) ; slat_p1=slata(i+1,j) ; clon_p1=clona(i+1,j) ; slon_p1=slona(i+1,j)

        x=clat_m1*clon_m1 ; y=clat_m1*slon_m1 ; z=slat_m1
        xt=x*clon0+y*slon0 ; yt=-x*slon0+y*clon0 ; zt=z
        yb=zt*clat0-xt*slat0
        xb=yt
        zb=xt*clat0+zt*slat0

        rlonb_m1=atan2(-yb,-xb)   !  the minus signs here are so line for m1 is directed same
        clonb_m1=cos(rlonb_m1)
        slonb_m1=sin(rlonb_m1)

        x=clat_p1*clon_p1 ; y=clat_p1*slon_p1 ; z=slat_p1
        xt=x*clon0+y*slon0 ; yt=-x*slon0+y*clon0 ; zt=z
        yb=zt*clat0-xt*slat0
        xb=yt
        zb=xt*clat0+zt*slat0
        rlonb_p1=atan2(yb,xb)
        clonb_p1=cos(rlonb_p1)
        slonb_p1=sin(rlonb_p1)
        crot=half*(clonb_m1+clonb_p1)
        srot=half*(slonb_m1+slonb_p1)
        gt%cos_beta_ref(i,j)=crot*clon0-srot*slon0
        gt%sin_beta_ref(i,j)=srot*clon0+crot*slon0
     end do
!               now do i=1 and i=gt%nlon at 1st order accuracy
     i=1

!   transform so pole is at rlat0,rlon0 and 0 meridian is tangent to earth latitude at rlat0,rlon0.

        clat0=clata(i,j) ; slat0=slata(i,j) ; clon0=clona(i,j) ; slon0=slona(i,j)
!    now obtain new coordinates for m1 and p1 points.

        clat_p1=clata(i+1,j) ; slat_p1=slata(i+1,j) ; clon_p1=clona(i+1,j) ; slon_p1=slona(i+1,j)

        x=clat_p1*clon_p1 ; y=clat_p1*slon_p1 ; z=slat_p1
        xt=x*clon0+y*slon0 ; yt=-x*slon0+y*clon0 ; zt=z
        yb=zt*clat0-xt*slat0
        xb=yt
        zb=xt*clat0+zt*slat0
        rlonb_p1=atan2(yb,xb)
        clonb_p1=cos(rlonb_p1)
        slonb_p1=sin(rlonb_p1)
        crot=clonb_p1
        srot=slonb_p1
        gt%cos_beta_ref(i,j)=crot*clon0-srot*slon0
        gt%sin_beta_ref(i,j)=srot*clon0+crot*slon0

     i=gt%nlon

!   transform so pole is at rlat0,rlon0 and 0 meridian is tangent to earth latitude at rlat0,rlon0.

        clat0=clata(i,j) ; slat0=slata(i,j) ; clon0=clona(i,j) ; slon0=slona(i,j)

!    now obtain new coordinates for m1 and p1 points.

        clat_m1=clata(i-1,j) ; slat_m1=slata(i-1,j) ; clon_m1=clona(i-1,j) ; slon_m1=slona(i-1,j)

        x=clat_m1*clon_m1 ; y=clat_m1*slon_m1 ; z=slat_m1
        xt=x*clon0+y*slon0 ; yt=-x*slon0+y*clon0 ; zt=z
        yb=zt*clat0-xt*slat0
        xb=yt
        zb=xt*clat0+zt*slat0

        rlonb_m1=atan2(-yb,-xb)   !  the minus signs here are so line for m1 is directed same
        clonb_m1=cos(rlonb_m1)
        slonb_m1=sin(rlonb_m1)

        crot=clonb_m1
        srot=slonb_m1
        gt%cos_beta_ref(i,j)=crot*clon0-srot*slon0
        gt%sin_beta_ref(i,j)=srot*clon0+crot*slon0
  end do
  deallocate(clata,slata,clona,slona)

end subroutine general_create_llxy_transform

subroutine general_tll2xy(gt,rlon,rlat,x,y,outside)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    general_tll2xy
!   prgmmr:  parrish
!
! abstract:  copy routines from gridmod.f90 to make a general purpose module which allows
!     conversion of earth lat lons to grid coordinates for any non-staggered rectangular 
!     orthogonal grid with known earth lat and lon of each grid point.
!     general_tll2xy converts earth lon-lat to x-y grid units of a
!           general regional rectangular domain.  Also, decides if
!           point is inside this domain.  As a result, there is
!           no restriction on type of horizontal coordinate for
!           a regional run, other than that it not have periodicity
!           or polar singularities.
!           This is done by first converting rlon, rlat to an
!           intermediate coordinate xtilde,ytilde, which has
!           precomputed pointers and constants for final conversion
!           to the desired x,y via 3 point inverse interpolation.
!           All of the information needed is derived from arrays
!           specifying earth latitude and longitude of every point
!           on the input grid.  Currently, the input x-y grid that
!           this is based on must be non-staggered.  This restriction
!           will eventually be lifted so we can run directly from
!           model grids that are staggered without first resorting
!           to interpolation of the guess to a non-staggered grid.
!
! program history log:
!   2010-10-28  parrish - initial documentation
!
!   input argument list:
!    gt          - type(llxy_cons) contains transform information to convert lat,lon to grid units x,y
!    rlon,rlat   - earth longitude and latitude (radians)
!
!   output argument list:
!    x,y         - grid coordinates in grid units.
!    outside     - if .false., then point is inside grid domain
!                    otherwise point is outside domain.
!
! attributes:
!   language: f90
!   machine:
!
!$$$ end documentation block

    use kinds, only: r_kind,i_kind
    use constants, only: one
    implicit none

    type(llxy_cons),intent(in) :: gt
    real(r_kind),intent(in   ) :: rlon
    real(r_kind),intent(in   ) :: rlat
    real(r_kind),intent(  out) :: x
    real(r_kind),intent(  out) :: y
    logical     ,intent(  out) :: outside

    real(r_kind) clon,slon,r_of_lat,xtilde,ytilde
    real(r_kind) dtilde,etilde
    real(r_kind) d1tilde,d2tilde,e1tilde,e2tilde,detinv
    integer(i_kind) itilde,jtilde
    integer(i_kind) i0,j0,ip,jp

!   first compute xtilde, ytilde

    clon=cos(rlon+gt%rlambda0)
    slon=sin(rlon+gt%rlambda0)
    r_of_lat=gt%pihalf+gt%sign_pole*rlat

    xtilde=gt%atilde_x+gt%btilde_x*r_of_lat*clon
    ytilde=gt%atilde_y+gt%btilde_y*r_of_lat*slon

!  next get interpolation information

    itilde=max(1,min(nint(xtilde),gt%nxtilde))
    jtilde=max(1,min(nint(ytilde),gt%nytilde))

    i0     =   gt%i0_tilde(itilde,jtilde)
    j0     =   gt%j0_tilde(itilde,jtilde)
    ip     =i0+gt%ip_tilde(itilde,jtilde)
    jp     =j0+gt%jp_tilde(itilde,jtilde)
    dtilde =xtilde-gt%xtilde0(i0,j0)
    etilde =ytilde-gt%ytilde0(i0,j0)
    d1tilde=(gt%xtilde0(ip,j0)-gt%xtilde0(i0,j0))*(ip-i0)
    d2tilde=(gt%xtilde0(i0,jp)-gt%xtilde0(i0,j0))*(jp-j0)
    e1tilde=(gt%ytilde0(ip,j0)-gt%ytilde0(i0,j0))*(ip-i0)
    e2tilde=(gt%ytilde0(i0,jp)-gt%ytilde0(i0,j0))*(jp-j0)
    detinv =one/(d1tilde*e2tilde-d2tilde*e1tilde)
    x = i0+detinv*(e2tilde*dtilde-d2tilde*etilde)
    y = j0+detinv*(d1tilde*etilde-e1tilde*dtilde)

    outside=x < gt%rlon_min_dd .or. x > gt%rlon_max_dd .or. &
            y < gt%rlat_min_dd .or. y > gt%rlat_max_dd

end subroutine general_tll2xy

subroutine general_txy2ll(gt,x,y,rlon,rlat)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    general_txy2ll
!   prgmmr:  parrish
!
! abstract: convert x-y grid units to earth lat-lon coordinates
!
! program history log:
!   2010-10-28 parrish - initial documentation
!
!   input argument list:
!    gt          - type(llxy_cons) contains transform information to convert x,y to lat,lon
!    x,y         - grid coordinates in grid units.
!
!   output argument list:
!    rlon,rlat   - earth longitude and latitude (radians)
!
! attributes:
!   language: f90
!   machine:
!
!$$$ end documentation block

    use kinds, only: r_kind,i_kind
    use constants, only: one
    implicit none

    type(llxy_cons),intent(in) :: gt
    real(r_kind),intent(in   ) :: x
    real(r_kind),intent(in   ) :: y
    real(r_kind),intent(  out) :: rlon
    real(r_kind),intent(  out) :: rlat

    real(r_kind) r_of_lat,xtilde,ytilde
    real(r_kind) dtilde,etilde,xbar,ybar
    real(r_kind) d1tilde,d2tilde,e1tilde,e2tilde
    integer(i_kind) i0,j0,ip,jp

    i0=nint(x)
    j0=nint(y)
    i0=max(1,min(i0,gt%nlon))
    j0=max(1,min(j0,gt%nlat))
    ip=i0+nint(sign(one,x-i0))
    jp=j0+nint(sign(one,y-j0))
    if(ip<1) then
       i0=2
       ip=1
    end if
    if(jp<1) then
       j0=2
       jp=1
    end if
    if(ip>gt%nlon) then
       i0=gt%nlon-1
       ip=gt%nlon
    end if
    if(jp>gt%nlat) then
       j0=gt%nlat-1
       jp=gt%nlat
    end if
    d1tilde=(gt%xtilde0(ip,j0)-gt%xtilde0(i0,j0))*(ip-i0)
    d2tilde=(gt%xtilde0(i0,jp)-gt%xtilde0(i0,j0))*(jp-j0)
    e1tilde=(gt%ytilde0(ip,j0)-gt%ytilde0(i0,j0))*(ip-i0)
    e2tilde=(gt%ytilde0(i0,jp)-gt%ytilde0(i0,j0))*(jp-j0)
    dtilde =d1tilde*(x-i0) +d2tilde*(y-j0)
    etilde =e1tilde*(x-i0) +e2tilde*(y-j0)
    xtilde =dtilde         +gt%xtilde0(i0,j0)
    ytilde =etilde         +gt%ytilde0(i0,j0)

    xbar=(xtilde-gt%atilde_x)*gt%btilde_xinv
    ybar=(ytilde-gt%atilde_y)*gt%btilde_yinv
    r_of_lat=sqrt(xbar**2+ybar**2)
    rlat=(r_of_lat-gt%pihalf)*gt%sign_pole
    rlon=atan2(ybar,xbar)-gt%rlambda0

end subroutine general_txy2ll

subroutine general_nearest_3(ilast,jlast,i0,j0,ip,jp,x,y,nx0,ny0,x0,y0)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    nearest_3
!   prgmmr:
!
! abstract: find closest 3 points to (x,y) on grid defined by x0,y0
!
! program history log:
!   2009-08-04  lueken - added subprogram doc block
!
!   input argument list:
!    ilast,jlast
!    nx0,ny0
!    x,y
!    x0,y0
!
!   output argument list:
!    ilast,jlast
!    i0,j0
!    ip,jp
!
! attributes:
!   language: f90
!   machine:
!
!$$$ end documentation block

  use kinds, only: r_kind,i_kind,i_byte
  implicit none

  integer(i_kind),intent(inout) :: ilast,jlast
  integer(i_kind),intent(  out) :: i0,j0
  integer(i_byte),intent(  out) :: ip,jp
  integer(i_kind),intent(in   ) :: nx0,ny0
  real(r_kind)   ,intent(in   ) :: x,y
  real(r_kind)   ,intent(in   ) :: x0(nx0,ny0),y0(nx0,ny0)
 
  real(r_kind) dista,distb,dist2,dist2min
  integer(i_kind) i,inext,j,jnext

  do
     i0=ilast
     j0=jlast
     dist2min=huge(dist2min)
     inext=0
     jnext=0
     do j=max(j0-1,1),min(j0+1,ny0)
        do i=max(i0-1,1),min(i0+1,nx0)
           dist2=(x-x0(i,j))**2+(y-y0(i,j))**2
           if(dist2<dist2min) then
              dist2min=dist2
              inext=i
              jnext=j
           end if
        end do
     end do
     if(inext==i0.and.jnext==j0) exit
     ilast=inext
     jlast=jnext
  end do

!  now find which way to go in x for second point

  ip=0
  if(i0==nx0)  ip=-1
  if(i0==1) ip=1
  if(ip==0) then
     dista=(x-x0(i0-1,j0))**2+(y-y0(i0-1,j0))**2
     distb=(x-x0(i0+1,j0))**2+(y-y0(i0+1,j0))**2
     if(distb<dista) then
        ip=1
     else
        ip=-1
     end if
  end if

!  repeat for y for 3rd point

  jp=0
  if(j0==ny0  ) jp=-1
  if(j0==1 ) jp=1
  if(jp==0) then
     dista=(x-x0(i0,j0-1))**2+(y-y0(i0,j0-1))**2
     distb=(x-x0(i0,j0+1))**2+(y-y0(i0,j0+1))**2
     if(distb<dista) then
        jp=1
     else
        jp=-1
     end if
  end if

  ilast=i0
  jlast=j0
    
end subroutine general_nearest_3

subroutine general_get_xytilde_domain(gt,nx0,ny0,rlons0,rlats0, &
                                  nx,ny,xminout,xmaxout,yminout,ymaxout)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    get_xytilde_domain
!   prgmmr:
!
! abstract:
!
! program history log:
!   2009-08-04  lueken - added subprogram doc block
!
!   input argument list:
!    nx0,ny0
!    rlons0,rlats0
!
!   output argument list:
!    nx,ny
!    xminout,xmaxout,yminout,ymaxout
!
! attributes:
!   language: f90
!   machine:
!
!$$$ end documentation block

   use kinds, only: r_kind,i_kind
   use constants, only: one, deg2rad,half,zero,r10
!  define parameters for xy domain which optimally overlays input grid

  implicit none

  type(llxy_cons),intent(inout) :: gt
  integer(i_kind),intent(in   ) :: nx0,ny0
  real(r_kind)   ,intent(in   ) :: rlons0(nx0,ny0),rlats0(nx0,ny0)

  integer(i_kind),intent(  out) :: nx,ny
  real(r_kind)   ,intent(  out) :: xminout,xmaxout,yminout,ymaxout

  real(r_kind),parameter:: r37=37.0_r_kind

  real(r_kind) area,areamax,areamin,extra,rlats0max,rlats0min,testlambda
  real(r_kind) xthis,ythis
  integer(i_kind) i,ip1,j,jp1,m

  real(r_kind) coslon0(nx0,ny0),sinlon0(nx0,ny0)
  real(r_kind) coslat0(nx0,ny0),sinlat0(nx0,ny0)
  real(r_kind) count,delbar
  real(r_kind) dx,dy,disti,distj,distmin,distmax
  real(r_kind) xmin,xmax,ymin,ymax

!  get range of lats for input grid

  rlats0max=maxval(rlats0) ; rlats0min=minval(rlats0)

!   assign hemisphere ( parameter sign_pole )

  if(rlats0min>-r37*deg2rad) gt%sign_pole=-one   !  northern hemisphere xy domain
  if(rlats0max< r37*deg2rad) gt%sign_pole= one   !  southern hemisphere xy domain


!   get optimum rotation angle rlambda0

  areamin= huge(areamin)
  areamax=-huge(areamax)
  do m=0,359
     testlambda=m*deg2rad
     xmax=-huge(xmax)
     xmin= huge(xmin)
     ymax=-huge(ymax)
     ymin= huge(ymin)
     do j=1,ny0,ny0-1
        do i=1,nx0
           xthis=(gt%pihalf+gt%sign_pole*rlats0(i,j))*cos(rlons0(i,j)+testlambda)
           ythis=(gt%pihalf+gt%sign_pole*rlats0(i,j))*sin(rlons0(i,j)+testlambda)
           xmax=max(xmax,xthis)
           ymax=max(ymax,ythis)
           xmin=min(xmin,xthis)
           ymin=min(ymin,ythis)
        end do
     end do
     do j=1,ny0
        do i=1,nx0,nx0-1
           xthis=(gt%pihalf+gt%sign_pole*rlats0(i,j))*cos(rlons0(i,j)+testlambda)
           ythis=(gt%pihalf+gt%sign_pole*rlats0(i,j))*sin(rlons0(i,j)+testlambda)
           xmax=max(xmax,xthis)
           ymax=max(ymax,ythis)
           xmin=min(xmin,xthis)
           ymin=min(ymin,ythis)
        end do
     end do
     area=(xmax-xmin)*(ymax-ymin)
     areamax=max(area,areamax)
     if(area<areamin) then
        areamin =area
        gt%rlambda0=testlambda
        xmaxout =xmax
        xminout =xmin
        ymaxout =ymax
        yminout =ymin
     end if
  end do


!   now determine resolution of input grid and choose nx,ny of xy grid accordingly
!                 (currently hard-wired at 1/2 the average input grid increment)

  do j=1,ny0
     do i=1,nx0
        coslon0(i,j)=cos(one*rlons0(i,j)) ; sinlon0(i,j)=sin(one*rlons0(i,j))
        coslat0(i,j)=cos(one*rlats0(i,j)) ; sinlat0(i,j)=sin(one*rlats0(i,j))
     end do
  end do

  delbar=zero
  count =zero
  do j=1,ny0-1
     jp1=j+1
     do i=1,nx0-1
        ip1=i+1
        disti=acos(sinlat0(i,j)*sinlat0(ip1,j)+coslat0(i,j)*coslat0(ip1,j)* &
                  (sinlon0(i,j)*sinlon0(ip1,j)+coslon0(i,j)*coslon0(ip1,j)))
        distj=acos(sinlat0(i,j)*sinlat0(i,jp1)+coslat0(i,j)*coslat0(i,jp1)* &
                  (sinlon0(i,j)*sinlon0(i,jp1)+coslon0(i,j)*coslon0(i,jp1)))
        distmax=max(disti,distj)
        distmin=min(disti,distj)
        delbar=delbar+distmax
        count=count+one
     end do
  end do
  delbar=delbar/count
  dx=half*delbar
  dy=dx

!   add extra space to computational grid to push any boundary problems away from
!     area of interest

  extra=r10*dx
  xmaxout=xmaxout+extra
  xminout=xminout-extra
  ymaxout=ymaxout+extra
  yminout=yminout-extra
  nx=1+(xmaxout-xminout)/dx
  ny=1+(ymaxout-yminout)/dy
 
end subroutine general_get_xytilde_domain

!-------------------------------------------------------------------------
!    NOAA/NCEP, National Centers for Environmental Prediction GSI        !
!-------------------------------------------------------------------------
!BOP
!
! !IROUTINE:  general_rotate_wind_ll2xy ---  Rotate earth vector wind
!
! !INTERFACE:
!
  subroutine general_rotate_wind_ll2xy(gt,u0,v0,u,v,rlon0,x,y)

! !USES:

    use kinds, only: r_kind,i_kind
    use constants, only: one,two,pi,rad2deg,one_tenth
    implicit none

! !INPUT PARAMETERS:

    type(llxy_cons),intent(in) :: gt
    real(r_kind),intent(in   ) :: u0,v0        ! earth wind component
    real(r_kind),intent(in   ) :: rlon0        ! earth   lon (radians)
    real(r_kind),intent(in   ) :: x,y          ! local x,y coordinate (grid units)

! !OUTPUT PARAMETERS:

    real(r_kind),intent(  out) :: u,v          ! rotated coordinate of winds

! !DESCRIPTION: to convert earth vector wind components to corresponding
!           local x,y coordinate
!
! !REVISION HISTORY:
!   2003-09-30  parrish
!   2004-05-13  kleist, documentation
!   2004-07-15  todling, protex-compliant prologue
!   2010-09-08  parrish, remove sign_pole variable--no longer needed, due to more accurate and
!                 robust computation of reference wind rotation angle defined by
!                 cos_beta_ref, sin_beta_ref.
!
! !REMARKS:
!   language: f90
!   machine:  ibm rs/6000 sp; SGI Origin 2000; Compaq/HP
!
! !AUTHOR:
!   parrish          org: np22                date: 2003-09-30
!
!EOP
!-------------------------------------------------------------------------

  real(r_kind) beta,delx,delxp,dely,delyp
  real(r_kind) sin_beta,cos_beta
  integer(i_kind) ix,iy

!  interpolate departure from longitude part of angle between earth positive east and local positive x

  ix=x
  iy=y
  ix=max(1,min(ix,gt%nlon-1))
  iy=max(1,min(iy,gt%nlat-1))
  delx=x-ix
  dely=y-iy
  delxp=one-delx
  delyp=one-dely
  cos_beta=gt%cos_beta_ref(ix  ,iy  )*delxp*delyp+gt%cos_beta_ref(ix+1,iy  )*delx *delyp+ &
           gt%cos_beta_ref(ix  ,iy+1)*delxp*dely +gt%cos_beta_ref(ix+1,iy+1)*delx *dely
  sin_beta=gt%sin_beta_ref(ix  ,iy  )*delxp*delyp+gt%sin_beta_ref(ix+1,iy  )*delx *delyp+ &
           gt%sin_beta_ref(ix  ,iy+1)*delxp*dely +gt%sin_beta_ref(ix+1,iy+1)*delx *dely
  beta=atan2(sin_beta,cos_beta)

!  now rotate;

  u= u0*cos(beta-rlon0)+v0*sin(beta-rlon0)
  v=-u0*sin(beta-rlon0)+v0*cos(beta-rlon0)

end subroutine general_rotate_wind_ll2xy

!-------------------------------------------------------------------------
!    NOAA/NCEP, National Centers for Environmental Prediction GSI        !
!-------------------------------------------------------------------------
!BOP
!
! !IROUTINE:  rotate_wind_xy2ll ---  Unrotate earth vector wind
!
! !INTERFACE:
!
  subroutine general_rotate_wind_xy2ll(gt,u,v,u0,v0,rlon0,x,y)

! !USES:

    use kinds, only: r_kind,i_kind
    use constants, only: one
    implicit none

! !INPUT PARAMETERS:

    type(llxy_cons),intent(in) :: gt
    real(r_kind),intent(in   ) :: u,v         ! rotated coordinate winds
    real(r_kind),intent(in   ) :: rlon0       ! earth   lon     (radians)
    real(r_kind),intent(in   ) :: x,y         ! rotated lon/lat (radians)

! !OUTPUT PARAMETERS:

    real(r_kind),intent(  out) :: u0,v0       ! earth winds

! !DESCRIPTION: rotate u,v in local x,y coordinate to u0,v0 in earth 
!           lat, lon coordinate
!
! !REVISION HISTORY:
!   2003-09-30  parrish
!   2004-05-13  kleist, documentation
!   2004-07-15  todling, protex-compliant prologue
!   2004-07-20  todling, fixed description
!   2010-09-08  parrish, remove sign_pole variable--no longer needed, due to more accurate and
!                 robust computation of reference wind rotation angle defined by
!                 cos_beta_ref, sin_beta_ref.
!
! !REMARKS:
!   language: f90
!   machine:  ibm rs/6000 sp; SGI Origin 2000; Compaq/HP
!
! !AUTHOR:
!   parrish          org: np22                date: 2003-09-30
!
!EOP
!-------------------------------------------------------------------------
  real(r_kind) beta,delx,delxp,dely,delyp
  real(r_kind) sin_beta,cos_beta
  integer(i_kind) ix,iy

!  interpolate departure from longitude part of angle between earth 
!  positive east and local positive x

  ix=x
  iy=y
  ix=max(1,min(ix,gt%nlon-1))
  iy=max(1,min(iy,gt%nlat-1))
  delx=x-ix
  dely=y-iy
  delxp=one-delx
  delyp=one-dely
  cos_beta=gt%cos_beta_ref(ix  ,iy  )*delxp*delyp+gt%cos_beta_ref(ix+1,iy  )*delx *delyp+ &
           gt%cos_beta_ref(ix  ,iy+1)*delxp*dely +gt%cos_beta_ref(ix+1,iy+1)*delx *dely
  sin_beta=gt%sin_beta_ref(ix  ,iy  )*delxp*delyp+gt%sin_beta_ref(ix+1,iy  )*delx *delyp+ &
           gt%sin_beta_ref(ix  ,iy+1)*delxp*dely +gt%sin_beta_ref(ix+1,iy+1)*delx *dely
  beta=atan2(sin_beta,cos_beta)

!  now rotate;

  u0= u*cos(beta-rlon0)-v*sin(beta-rlon0)
  v0= u*sin(beta-rlon0)+v*cos(beta-rlon0)

end subroutine general_rotate_wind_xy2ll

end module general_tll2xy_mod

subroutine merge_grid_e_to_grid_a_initialize(region_lat_e,region_lon_e,region_lat_a,region_lon_a, &
                                   nlat_e,nlon_e,nlat_a,nlon_a,nord_e2a,nord_blend,nmix,gt_e,gt_a,p_e2a)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    merge_grid_e_to_grid_a_initialize
!   prgmmr:  parrish
!
! abstract:  initialize parameters needed for routines merge_grid_e_to_grid_a and merge_vgrid_e_to_vgrid_a.
!
! program history log:
!   2010-10-28  parrish - initial documentation
!   2012-03-01  tong - modified the way to call create_egrid2points_slow. If nmix <= 0, the orginal 
!               way to call the subroutine will cause segmentation fault
!
!   input argument list:
!    region_lat_e - earth lats for e grid (radians) 
!    region_lon_e - earth lons for e grid (radians)
!    region_lat_a - earth lats for a grid (radians)
!    region_lon_a - earth lons for a grid (radians)
!    nlat_e       - "y" dimension of e grid
!    nlon_e       - "x" dimension of e grid
!    nlat_a       - "y" dimension of a grid
!    nlon_a       - "x" dimension of a grid
!    nord_e2a     - order of accuracy of interpolation from grid e to grid a (1=bilinear, 2=biquadratic, etc)
!    nord_blend   - order of continuity of blend function (1=continuous 1st derivative, etc)
!    nmix         - width of blend zone on edge of e grid in e grid units.
!
!   output argument list:
!    gt_e         - navigation structure variable for grid e
!    gt_a         - navigation structure variable for grid a
!    p_e2a        - structure variable containing interpolation info for interpolation from e grid to a grid
!
! attributes:
!   language: f90
!   machine:
!
!$$$ end documentation block

  use kinds, only: r_kind,i_kind
  use general_tll2xy_mod, only: llxy_cons,general_create_llxy_transform,general_tll2xy
  use egrid2agrid_mod, only: egrid2agrid_parm,create_egrid2points_slow
  use mpimod, only: mype
  use constants, only: zero,one
  implicit none

  type(llxy_cons),intent(inout):: gt_e,gt_a
  type(egrid2agrid_parm),intent(inout) :: p_e2a
  integer(i_kind),intent(in):: nlat_e,nlon_e,nlat_a,nlon_a,nord_e2a,nord_blend,nmix
  real(r_kind),intent(in)::region_lat_e(nlat_e,nlon_e),region_lon_e(nlat_e,nlon_e)
  real(r_kind),intent(in)::region_lat_a(nlat_a,nlon_a),region_lon_a(nlat_a,nlon_a)

  integer(i_kind) i,j,np
  real(r_kind),allocatable::xa_e(:,:),ya_e(:,:),xe(:),ye(:)
  logical outside
  real(r_kind) diffmax,range_lat,range_lon

! define navigation for e grid and a grid

  if(nlat_e == nlat_a .and. nlon_e == nlon_a) then
     range_lat=max(abs(region_lat_a(nlat_a,1)-region_lat_a(1,1)),abs(region_lat_a(nlat_a,nlon_a)-region_lat_a(1,nlon_a)))
     range_lat=max(range_lat,abs(region_lat_a(nlat_a,nlon_a/2)-region_lat_a(1,nlon_a/2)))
     range_lat=max(range_lat,abs(region_lat_e(nlat_e,1)-region_lat_e(1,1)))
     range_lat=max(range_lat,abs(region_lat_e(nlat_e,nlon_e)-region_lat_e(1,nlon_e)))
     range_lat=max(range_lat,abs(region_lat_e(nlat_e,nlon_e/2)-region_lat_e(1,nlon_e/2)))
     if(nlat_a == 1) range_lat=one
     range_lon=max(abs(region_lon_a(1,nlon_a)-region_lon_a(1,1)),abs(region_lon_a(nlat_a,nlon_a)-region_lon_a(nlat_a,1)))
     range_lon=max(range_lon,abs(region_lon_a(nlat_a/2,nlon_a)-region_lon_a(nlat_a/2,1)))
     range_lon=max(range_lon,abs(region_lon_e(1,nlon_e)-region_lon_e(1,1)))
     range_lon=max(range_lon,abs(region_lon_e(nlat_e,nlon_e)-region_lon_e(nlat_e,1)))
     range_lon=max(range_lon,abs(region_lon_e(nlat_e/2,nlon_e)-region_lon_e(nlat_e/2,1)))
     if(nlon_a == 1) range_lon=one
     diffmax=zero
     do j=1,nlon_a
        do i=1,nlat_a
           diffmax=max(diffmax,abs(region_lat_a(i,j)-region_lat_e(i,j))/range_lat)
        end do
     end do
     do j=1,nlon_a
        do i=1,nlat_a
           diffmax=max(diffmax,abs(region_lon_a(i,j)-region_lon_e(i,j))/range_lon)
        end do
     end do
     if(diffmax < .0000001_r_kind)then
        p_e2a%identity=.true.
        return
     end if
  end if

  call general_create_llxy_transform(region_lat_e,region_lon_e,nlat_e,nlon_e,gt_e)
  call general_create_llxy_transform(region_lat_a,region_lon_a,nlat_a,nlon_a,gt_a)

!    obtain xa_e, ya_e, the coordinates of the a grid points in e-grid measure.

  allocate(xa_e(nlat_a,nlon_a),ya_e(nlat_a,nlon_a))
  do j=1,nlon_a
     do i=1,nlat_a
        call general_tll2xy(gt_e,region_lon_a(i,j),region_lat_a(i,j),xa_e(i,j),ya_e(i,j),outside)
     end do
  end do

  allocate(xe(nlon_e),ye(nlat_e))
  do j=1,nlon_e
     xe(j)=j
  end do
  do i=1,nlat_e
     ye(i)=i
  end do
  np=nlat_a*nlon_a
  if(nord_blend > 0 .and. nmix > 0)then
     call create_egrid2points_slow(np,ya_e,xa_e,nlat_e,ye,nlon_e,xe,nord_e2a,p_e2a,nord_blend,nmix)
  else
     call create_egrid2points_slow(np,ya_e,xa_e,nlat_e,ye,nlon_e,xe,nord_e2a,p_e2a)
  end if

  deallocate(xe,ye,xa_e,ya_e)

end subroutine merge_grid_e_to_grid_a_initialize

subroutine merge_grid_e_to_grid_a(e_in,a_in,a_out,gt_e,gt_a,p_e2a)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    merge_grid_e_to_grid_a  merge scalar field from one grid to another
!   prgmmr: parrish          org: np22                date: 2010-11-02
!
! abstract: interpolate input scalar field e_in from e grid to a grid, blend with existing
!             field a_in, and return result in a_out.
!
! program history log:
!   2010-11-02  parrish, initial documentation
!
!   input argument list:
!      e_in:             input field on e-grid
!      a_in:             input field on a-grid
!      gt_e:             structure variable containing navigation info and dimensions for e_in
!      gt_a:             same for a_in, a_out
!      p_e2a:            structure variable containing info for interpolating from e-grid to a-grid.
!
!   output argument list:
!      a_out:            blended field with interpolated e-grid field smoothly replacing
!                           region of a-grid covered by e-grid.
!              
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$ end documentation block
  use kinds, only: r_kind,i_kind
  use general_tll2xy_mod, only: llxy_cons
  use egrid2agrid_mod, only: egrid2agrid_parm,egrid2points
  use constants, only: one
  implicit none

  type(llxy_cons),intent(in):: gt_e,gt_a
  type(egrid2agrid_parm),intent(in):: p_e2a
  real(r_kind),intent(in),dimension(gt_e%nlat,gt_e%nlon):: e_in
  real(r_kind),intent(in),dimension(gt_a%nlat,gt_a%nlon):: a_in
  real(r_kind),intent(out),dimension(gt_a%nlat,gt_a%nlon):: a_out

  integer(i_kind) i,ii,j

  call egrid2points(p_e2a,e_in,a_out)
  ii=0
  do j=1,gt_a%nlon
     do i=1,gt_a%nlat
        ii=ii+1
!   blend smoothly with a_in
        a_out(i,j)=p_e2a%blend(ii)*a_out(i,j)+(one-p_e2a%blend(ii))*a_in(i,j)
!        a_out(i,j)=blend(i,j)*a_out(i,j)+(one-blend(i,j))*a_in(i,j)
     end do
  end do

end subroutine merge_grid_e_to_grid_a

subroutine merge_vgrid_e_to_vgrid_a(eu_in,ev_in,au_in,av_in,au_out,av_out,gt_e,gt_a,p_e2a)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    merge_vgrid_e_to_vgrid_a  merge vector field from one grid to another
!   prgmmr: parrish          org: np22                date: 2010-11-02
!
! abstract: interpolate input vector field eu_in,ev_in from e grid to a grid, blend with existing
!             vector field au_in, av_in, and return result in au_out,av_out.
!
! program history log:
!   2010-11-02  parrish, initial documentation
!
!   input argument list:
!      eu_in:            input u component on e-grid
!      ev_in:            input v component on e-grid
!      au_in:            input u component on a-grid
!      av_in:            input v component on a-grid
!      gt_e:             structure variable containing navigation info and dimensions for eu_in, ev_in
!      gt_a:             same for au_in, av_in, au_out,av_out
!      p_e2a:            structure variable containing info for interpolating from e-grid to a-grid.
!
!   output argument list:
!      au_out:           blended u component with interpolated e-grid u component smoothly replacing
!                           region of a-grid covered by e-grid.
!      av_out:           same for v component
!              
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$ end documentation block
  use kinds, only: r_kind,i_kind
  use general_tll2xy_mod, only: llxy_cons,general_rotate_wind_ll2xy,general_rotate_wind_xy2ll
  use egrid2agrid_mod, only: egrid2agrid_parm,egrid2points
  use constants, only: one
  implicit none

  type(llxy_cons),intent(in):: gt_e,gt_a
  type(egrid2agrid_parm),intent(in):: p_e2a
  real(r_kind),intent(in),dimension(gt_e%nlat,gt_e%nlon):: eu_in,ev_in
  real(r_kind),intent(in),dimension(gt_a%nlat,gt_a%nlon):: au_in,av_in
  real(r_kind),intent(out),dimension(gt_a%nlat,gt_a%nlon):: au_out,av_out

  integer(i_kind) i,ii,j
  real(r_kind) xa,ya,au_earth,av_earth

  call egrid2points(p_e2a,eu_in,au_out)
  call egrid2points(p_e2a,ev_in,av_out)
  ii=0
  do j=1,gt_a%nlon
     xa=j
     do i=1,gt_a%nlat
        ya=i
        ii=ii+1
!    rotate vector from e grid orientation to earth lat-lon orientation
        call general_rotate_wind_xy2ll(gt_e,au_out(i,j),av_out(i,j),au_earth,av_earth, &
                                       gt_a%region_lon(i,j),p_e2a%xa_e(ii),p_e2a%ya_e(ii))
!    rotate vector from earth lat-lon orientation to a grid orientation
        call general_rotate_wind_ll2xy(gt_a,au_earth,av_earth,au_out(i,j),av_out(i,j), &
                                       gt_a%region_lon(i,j),xa,ya)
!   blend smoothly with au_in, av_in
        au_out(i,j)=p_e2a%blend(ii)*au_out(i,j)+(one-p_e2a%blend(ii))*au_in(i,j)
        av_out(i,j)=p_e2a%blend(ii)*av_out(i,j)+(one-p_e2a%blend(ii))*av_in(i,j)
!        au_out(i,j)=blend(i,j)*au_out(i,j)+(one-blend(i,j))*au_in(i,j)
!        av_out(i,j)=blend(i,j)*av_out(i,j)+(one-blend(i,j))*av_in(i,j)
     end do
  end do

end subroutine merge_vgrid_e_to_vgrid_a