    module grid_indexing
!
!  SYNOPSIS
!    use grid_indexing
!
!  DESCRIPTION
!
!  GRID_INDEXING contains utility indexing routines for managing the
!  model grids.  Wherever possible, these routines are declared PURE
!  to make them candidates for inlining by optimizing compilers.
!
!  REQUIRED PACKAGES
!   assim_vars, grid_setup, interpolate - in source code directory
!
!  AUTHOR
!   Eric Kostelich, kostelich@asu.edu, 480-065-5006
!
!  MODIFICATION HISTORY
!  10/13/05 - added optional OVER_POLE flag to XY_ADJUST so we can know when
!             an increment takes us over one of the poles.
!  07/30/05 - corrected XY_ADJUST so that patches over the North Pole are
!             indexed correctly.
!
    use grid_setup
    use resolution,only:LONSPERLAT,LONF,LATG,RLATS
    implicit none
    contains
!-----------------------------------------------------------------------------
    pure integer function patchsize(loc)
!  PATCHSIZE - return the number of dynamical variables that are assimilated
!  into a patch centered at LOC.  
!
    type(patch_index),intent(in)::loc
!
    integer::offset,low,high,ix,iy
    type(patch_index)::newloc
!
!  PATCH_OFFSET returns the number of variables PLUS 1 that appear
!  before the upper right corner of the patch.  See comments in PATCH_OFFSET
!  regarding the treatment of surface pressure.
!
    call patch_radius(loc,ix,iy)
    offset=patch_offset(loc,ix,iy)-1
    call xy_adjust(loc,ix,iy,newloc,low,high)
    offset=offset+(NAVARS-1)*(high-low+1)+1  ! and add pressure back
    patchsize=offset
    return
    end function patchsize
!-----------------------------------------------------------------------------
    pure subroutine xy_adjust(loc,inclon,inclat,newloc,low,high,over_pole)
!  XY_ADJUST - convenience routine for indexing an offset (INCLON,INCLAT)
!  from a grid point at LOC.  This version handles points near
!  the poles and the 0-degree line.
!  LOC    :  (LON,LAT,Z)
!  NEWLOC := (LON+INCLON,LAT+INCLAT,Z), with appropriate adjustments near poles
!            and 0-degree line.
!  LOW, HIGH := beginning and ending vertical levels centered at NEWLOC.
!  OVER_POLE := if present, set to TRUE if the incrementation moves us over
!            one of the poles.
!  We provide a simple modulus, assuming that |LON| and |LAT| are
!  not too large.
!
    implicit none
    type(patch_index),intent(in)::loc
    type(patch_index),intent(out)::newloc
    integer,intent(in)::inclon,inclat
    integer,intent(out),optional::low,high
    logical,intent(out),optional::over_pole
!
    integer::zlow,zhigh
!
    newloc%x=loc%x+inclon
    newloc%y=loc%y+inclat
    newloc%z=loc%z
!
    if(newloc%y.lt.1) then     ! adjust at south pole: 0 => 1, -1 => 2, etc.
       newloc%y=1-newloc%y
       newloc%x=newloc%x+LONF/2  ! go across pole 180 degrees
       if(present(over_pole)) over_pole=.true.
    else if(newloc%y.gt.LATG) then ! adjust at north pole: LATG+1 => LATG etc.
       newloc%y=2*LATG-newloc%y+1
       newloc%x=newloc%x+LONF/2  ! go across pole 180 degrees
       if(present(over_pole)) over_pole=.true.
    else
       if(present(over_pole)) over_pole=.false.
    endif
!
    if(newloc%x.lt.1) then
       newloc%x=newloc%x+LONF
    else if(newloc%x.gt.LONF) then
       newloc%x=newloc%x-LONF
    endif
    if(present(low).or.present(high)) then
       call zdiam(newloc,zlow,zhigh)
       if(present(low)) low=zlow
       if(present(high)) high=zhigh
    endif
    return
    end subroutine xy_adjust
!-----------------------------------------------------------------------------
    pure integer function patch_offset(loc,dx,dy)
!  PATCH_OFFSET - For an offset (dx,dy) from the center of a patch
!  centered at (IX,IY), find the index that corresponds to the start of
!  the dynamical variables for the associated grid point.
!  We allow for the possibility that the height of each cube varies by
!  coordinate of the patch array.
!  both longitude and latitude.
!  Arguments:
!  LOC     :  (lon,lat,z) coords of center of patch
!  DX, DY  :  (lon,lat) offset from center
!  OFFSET  := points to the index of the surface pressure component
!             at (lon+dx,lat+dy,z)
!
    implicit none
    type(patch_index),intent(in)::loc
    integer,intent(in)::dx,dy
!
    integer::j,k,low,high,nlev,offset,xrad,yrad
    type(patch_index)::newloc
!
!  Return a value that's sure to cause a crash if arguments are out of range
!
    call patch_radius(loc,xrad,yrad)
    if(abs(dx).gt.xrad.or.abs(dy).gt.yrad) then
       offset=huge(1)
       return
    endif
!
!  Surface pressure is a special case: we always include it in every patch,
!  even those that don't contain the surface.  All the other dynamical
!  variables included in the assimilation appear at each level of the
!  atmosphere.  NAVARS is the number of such dynamical variables,
!  including surface pressure, written out to the model grid by preanalysis.
!  In the computations below, we include the surface pressure once (at the
!  "lowest" level of the patch) followed by all other dynamical variables
!  at each level contained in the patch.
!
!  First sum up across all columns to the left of that defined by DX...
!
    offset=0
    do j=-xrad,dx-1
       do k=-yrad,yrad
          call xy_adjust(loc,j,k,newloc,low,high)
          nlev=high-low+1
          offset=offset+(NAVARS-1)*nlev+1  ! include pressure => +1
       enddo
    enddo
!
!  ... and all rows below that defined by DY.
!
    do k=-yrad,dy-1
       call xy_adjust(loc,dx,k,newloc,low,high)
       nlev=high-low+1
       offset=offset+(NAVARS-1)*nlev+1  ! include pressure => +1
    enddo
    patch_offset=offset+1   ! so we point to first component
    return
    end function patch_offset
!-----------------------------------------------------------------------------
    pure subroutine center_column_height(loc,height,offset)
!  CENTER_COLUMN_HEIGHT - return the number of variables in any array
!  returned by copy_from_grid corresponding to the center 
!  (longitude,latitude) grid point.
!  OFFSET  :=  Optional argument giving the offset corresponding to the first
!  variable in the center column of the patch.
!
    implicit none
    type(patch_index),intent(in)::loc
    integer,intent(out)::height
    integer,intent(out),optional::offset
!
    integer::zlow,zhigh
!
    if(present(offset)) offset=patch_offset(loc,0,0)
    call zdiam(loc,zlow,zhigh)
    height=(NAVARS-1)*(zhigh-zlow+1)+1
    return
    end subroutine center_column_height
!-----------------------------------------------------------------------------
    elemental logical function valid_dynamical_type(t)
!   VALID_DYNAMICAL_TYPE - return TRUE if T presents a dynamical variable
!   that is stored in the model grid.
!   We'll need to update this routine as more dynamical variables are
!   assimilated.
!
    integer,intent(in)::t
!
    valid_dynamical_type=t.ge.1.and.t.le.VWIND_TYPE
    return
    end function valid_dynamical_type
!-----------------------------------------------------------------------------
    elemental integer function signum(x)
!  SIGNUM - the usual signum function: return -1 if X is negative, 1 if
!  positive, and 0 otherwise.  The latter condition is important: we don't
!  want to interpolate measurements that coincide with a model grid point.
!
    real(R),intent(in)::x
!
    if(x.eq.0.0) then
       signum=0
    else
       signum=merge(1,-1,x.gt.0.0)
    endif
    return
    end function signum
!-----------------------------------------------------------------------------
    elemental real function xgrid(lambda)
!  XGRID - Given a longitude location in degrees (0 <= LAMBDA < 360),
!  return the same location relative to the model grid as an
!  interpolated subscript value.  That is, if LAMBDA is halfway between
!  the model longitude points J and J+1, then we return J+0.5 (as a real number).
!  The model longitude coordinates in this version are equally spaced,
!  with 360/LONF degrees per grid point.  Our convention in this code
!  is the 0 degrees longitude corresponds to the model grid location x=1;
!  longitudes in [360-360/LONF, 360) map to [LONF,LONF+1) in the model grid.
!  If LAMBDA is not in [0, 360), then return -1.
!
    use interpolate
    real,intent(in)::lambda
!
    if(lambda.lt.0.0.or.lambda.ge.360.0) then
       xgrid=-1.0   ! flag invalid longitude
    else
       xgrid=lambda*(LONF/360.0)+1.0  ! offset: 0 degrees -> x=1
    endif
    return
    end function xgrid
!-----------------------------------------------------------------------------
    elemental real function ygrid(phi)
!  YGRID - Given a latitude location in degrees (-90 <= PHI <= 90),
!  return the same location relative to the model grid as an
!  interpolated subscript value.  That is, if LAMBDA is halfway between
!  the model latitude points J and J+1, then we return J+0.5 (as a real number).
!  The model latitude coordinates in this version appear in a Gaussian grid
!  whose spacings depend on the model resolution; see the module resolution
!  in the support directory.
!  The GFS model grid does not go all the way to the poles.  The convention
!  in this code is that the south pole corresponds to YGRID=0 and the north
!  pole to YGRID=LATG+1.  Otherwise, if PHI < -90 or PHI > 90, then return -1.
!
    use interpolate
    real,intent(in)::phi
!
    real,parameter::SBDY(2)=(/0.0,RLATS(1)/)  ! latitudes to south pole
    real,parameter::NBDY(2)=(/RLATS(LATG),90.0/) ! latitudes to north pole.
    real::y
!
    if(phi.lt.-90.0.or.phi.gt.90.0) then
       ygrid=-1.0   ! flag invalid latitude
       return
    endif
!
    y=index_interp(phi,RLATS,LATG)
    if(y.lt.1.0) then  ! interpolate points S of the southernmost grid bdy
       y=index_interp(phi,SBDY,2)-1.0
    else if(y.gt.LATG) then  ! interpolate points N of northernmost grid bdy
       y=index_interp(phi,NBDY,2)+(LATG-1)
    endif
    ygrid=y
    return
    end function ygrid
!-----------------------------------------------------------------------------
    end module grid_indexing
