    module grid_setup
!
!  SYNOPSIS
!    use grid_setup
!
!  DESCRIPTION
!   GRID_SETUP contains the initialization and indexing
!   routines for managing the model grids.
!
!  REQUIRED PACKAGES
!   grid_types, covar_setup - in main source directory
!   ferror, scalings - in support directory
!
!  ROUTINES IMPLEMENTED
!
!   pure subroutine zdiam(loc,low,high)
!   pure subroutine model_grid_boundary(south,north)
!   pure subroutine patch_radius(loc,lon,lat,z)
!   pure integer function seqid()
!   pure subroutine assimilation_boundary(south,north)
!   subroutine set_model_grid_boundary(south,north,psouth,pnorth)
!   subroutine set_patch_radius(lon,lat)
!   subroutine vertical_localization_setup(zfile)
!   subroutine set_sequence_id(id)
!
!  AUTHOR
!   Eric Kostelich, kostelich@asu.edu, 480-965-5006
!
!  MODIFICATION HISTORY
!   10/12/04 - first version derived from old grid_manager code.
!
    use covar_setup
    use ferror
!
!  Patch radius defaults.
!
    integer,parameter,private::XWIN_DEFAULT=3
    integer,parameter,private::YWIN_DEFAULT=3
!
!  The prefix X refers to quantities that index the longitudinal
!  ("horizontal") direction; Y for latitudinal.
!  Data are assimilated in rectangular patches of radius XWIN and YWIN.
!  Vertical radii are level-dependent.
!
    integer,save,private::xwin=XWIN_DEFAULT
    integer,save,private::ywin=YWIN_DEFAULT
!
!  We assimilate data in latitudinal strips.  When the program is run
!  on multiple processors, each instance gets a grid and an observation
!  file with header information that identifies the range ot longitudes
!  to which the data refer.  YSOUTH and YNORTH define the boundaries of
!  the model grid, which by default is the whole globe.  ASOUTH and ANORTH
!  identify the southern- and northernmost latitudes at which patches may
!  be centered for the data assimilation procedure.
!
    integer,save,private::ysouth=1
    integer,save,private::ynorth=LATG
    integer,save,private::asouth=1
    integer,save,private::anorth=LATG
!
!  The MRF code outputs values in arrays whose leading dimension is LONFD,
!  which is greater than LONF (the size of the longitude grid).  We
!  allocate, read, and write global arrays whose longitudinal dimension
!  is at most 1:LONF.
!
!  All grid arrays are indexed by (latitude,longitude,bred_vector).
!  The leading dimension is max(1,YLEFT-NPAD):min(LATG,YRIGHT+NPAD)
!  in all cases (including the global case).
!  The service routine LEADING_DIMENSIONS returns these values.
!
!  The convenience routine XY_ADJUST handles longitude wrapping at the
!  0-degree line and at the poles, when patches encompass grid points
!  opposite the pole.
!
!  To help keep track of the correspondence between grid and observation
!  arrays, the splitting routines insert "sequence ids" into the header
!  portion of the data files.  We store it here.  A zero value means
!  the the sequence id is undefined or will not be used.  If nonzero, we
!  usually use a date field in the form yyyyddmmhh.
!
    integer,save,private::sequence_id=0
!
!  The vertical thickness of each patch is determined by zrad, whose values
!  are read in by read_vert_info_file.
!
    integer,save,private::zrad(LEVS,LATG)=LEVS ! default: no vert localization
    contains
! -------------------
    pure subroutine zdiam(loc,low,high)
!  ZDIAM - return the lower and upper indices for the vertical levels
!  contained in the patch LOC.
!
    implicit none
    type(patch_index),intent(in)::loc
    integer,intent(out)::low,high
!
    integer::z
!
    z=zrad(loc%z,loc%y)
    low=max(loc%z-z,1)
    high=min(loc%z+z,LEVS)
    return
    end subroutine zdiam
! -------------------
    pure subroutine model_grid_boundary(south,north)
!  MODEL_GRID_BOUNDARY - return the latitude boundaries of the model grids.
!
    integer,intent(out)::south,north
!
    south=ysouth
    north=ynorth
    return
    end subroutine model_grid_boundary
! -------------------
    pure subroutine assimilation_boundary(south,north)
!  ASSIMILATION_BOUNDARY - return the starting and stopping latitudes at which
!  patches may be centered for the data assimilation.  These boundaries
!  must be in the closed interior of the model and observation
!  grid latitude boundaries.
!  Note: SET_MODEL_GRID_BOUNDARY must be called first.
!
    integer,intent(out)::south,north
!
    south=asouth
    north=anorth
    return
    end subroutine assimilation_boundary
! -------------------
    pure subroutine patch_radius(loc,lon,lat,z)
!  PATCH_RADIUS - return the patch radius in longitudinal (X),
!  latitudinal (Y) coordinates, and/or vertical (Z) coordinates
!  for a patch centered at LOC.  Right now we assume a constant patch radius,
!  but in the future we can allow for patch sizes that vary with location.
!  Service routine.
!  Errors signaled by 0 return values for lon,lat,z.
!
    implicit none
    type(patch_index),intent(in),optional::loc
    integer,intent(out),optional::lon,lat,z
!
    if(present(lon)) lon=xwin
    if(present(lat)) lat=ywin
    if(present(z)) then
       if(present(loc)) then
          z=zrad(loc%z,loc%y)
       else
          z=0
       endif
    endif
    return
    end subroutine patch_radius
! -------------------
    pure integer function seqid()
!  SEQID - return the sequence id of a split grid file; 0 if none appears
!  in the file or if the file contains the whole globe.
!
    seqid=sequence_id
    return
    end function seqid
! -------------------
    subroutine set_model_grid_boundary(south,north,psouth,pnorth)
!  SET_MODEL_GRID_BOUNDARY - define the southern- and northernmost latitudes
!  contained in the model grid.  Cannot be called more than once.
!  SOUTH, NORTH : define the latitude strip that contains all model grid points.
!  PSOUTH, PNORTH : the latitude strip at which patches may be centered
!    for the data assimilation procedure.
!
    integer,intent(in)::south,north,psouth,pnorth
!
    logical,save::grid_defined=.false.
!
    if(grid_defined) then
       call fabort('set_strip_boundary: cannot reset boundary')
    else if(south.lt.1.or.north.gt.LATG) then
       call fabort('set_strip_boundary: illegal parameter value')
    endif
    ysouth=south
    ynorth=north
    asouth=psouth
    anorth=pnorth
    grid_defined=.true.
    return
    end subroutine set_model_grid_boundary
! -------------------
    subroutine set_patch_radius(lon,lat)
!  SET_PATCH_RADIUS - set the patch radius in longitudinal (X) and
!  latitudinal (Y) coordinates. 
!  For the moment, this is constant for the entire globe.  In the future,
!  we might consider radii that vary with location; these might be set
!  in the vert_info file.
!
    implicit none
    integer,intent(in),optional::lon,lat
!
    if(present(lon)) xwin=lon
    if(present(lat)) ywin=lat
    if(xwin.lt.1.or.ywin.lt.1) call fabort('set_patch_radius: illegal argument')
    return
    end subroutine set_patch_radius
! -------------------
    subroutine vertical_localization_setup(zfile)
    character(*),intent(in)::zfile
!
    call read_vert_info_file(zfile,zrad)
    return
    end subroutine vertical_localization_setup
! -------------------
    subroutine set_sequence_id(id)
!  SET_SEQUENCE_ID - If ID is negative, then create an ID.  Otherwise, set the
!  sequence id to the value of ID.
!
    integer,intent(in)::id
!
    logical,save::called=.false.
!
    if(called) call fabort('set_sequence_id: cannot be called twice')
    called=.true.
    if(id.ge.0) then
       sequence_id=id
    else
       call system_clock(sequence_id)
    endif
    return
    end subroutine set_sequence_id
! -------------------
    end module grid_setup
