text:
! This code may (or may not) be part of the COAMPS distribution,
! So it is not protected by the DART copyright agreement.
!
! DART $Id$

module coamps_grid_mod

!------------------------------
! MODULE:       coamps_grid_mod
! AUTHOR:       T. R. Whitcomb
!               Naval Research Laboratory
! DART VERSION: Jamaica
!
! Module containing the definition of a COAMPS grid data structure
! and routines dealing with COAMPS grid information.
!------------------------------

  use coamps_intrinsic_mod, only : ij2ll,                &
                                   ll2ij

  use coamps_util_mod,      only : C_REAL,               &
                                   check_io_status,      &
                                   check_alloc_status,   &
                                   check_dealloc_status, &
                                   fix_for_platform

  use location_mod,  only : get_location,  &
                            location_type, &
                            set_location

  use utilities_mod, only : do_output,     &
                            E_ERR,         &
                            E_WARN,        &
                            error_handler, &
                            get_unit

  use types_mod,     only : r8

  implicit none

  private

  !------------------------------
  ! BEGIN PUBLIC INTERFACE
  !------------------------------

  ! COAMPS grid data structure
  public :: coamps_grid
  public :: initialize_grid

  ! Location conversion functions
  public :: gridpt_to_latlon
  public :: latlon_to_gridpt
  public :: gridpt_to_location
  public :: location_to_gridpt
  public :: grid_ij_to_vector_index

  ! Information about the grid size
  public :: get_grid_dims
  public :: get_grid_field_size
  public :: get_grid_num_levels
  public :: check_ij_within_grid

  ! Information about horizontal coordinates
  public :: get_grid_delta_x
  public :: get_grid_delta_y

  ! Information about vertical coordinates
  public :: get_grid_msigma
  public :: get_grid_wsigma
  public :: get_grid_dsigmaw
  public :: get_terrain_height_at_points

  ! Domain decomposition tools
  public :: initialize_decomposition
  public :: decompose_domain
  public :: get_grid_iminf
  public :: get_grid_imaxf
  public :: get_grid_jminf
  public :: get_grid_jmaxf

  ! Debugging output
  public :: dump_grid_info

  !------------------------------
  ! END PUBLIC INTERFACE
  !------------------------------

  !------------------------------
  ! BEGIN EXTERNAL INTERFACES
  !------------------------------
  ! [none]
  !------------------------------
  ! END EXTERNAL INTERFACES
  !------------------------------

  !------------------------------
  ! BEGIN TYPES AND CONSTANTS
  !------------------------------

  ! Derived data type for the COAMPS grid information - this is all
  ! contained in the datahd_*_infofld files that the COAMPS analysis
  ! produces, so we can just read that rather than processing a
  ! namelist.  Most of these are pretty straighforward, but:
  !  anchor_i - the i-grid point where the grid is anchored to parent
  !  anchor_j - the j-grid point where the grid is anchored to parent
  !  iref     - the i-grid point where the grid is anchored to Earth
  !  jref     - the j-grid point where the grid is anchored to Earth
  !  dsigma   - the spacing between mass sigma levels
  !  msigma   - actual sigma values at mass levels
  !  terrain  - terrain height
  !  [ij]minf - the minimum grid point for process N
  !  [ij]maxf - the maximum grid point for process N
  ! NOTE: For now, "nests" is currently unused and the pts, anchor,
  ! and [ij]ref variables are scalars - when multiple grid nests are
  ! supported, these will be vectors.  Currently only a single grid
  ! is supported to make my life much, much, much easier.
  type :: coamps_grid
     real(kind=r8)                              :: pres_lvls
     real(kind=r8)                              :: sigm_lvls
     real(kind=r8)                              :: map_proj
     real(kind=r8)                              :: stdlat1
     real(kind=r8)                              :: stdlat2
     real(kind=r8)                              :: stdlon
     real(kind=r8)                              :: reflat
     real(kind=r8)                              :: reflon
     real(kind=r8)                              :: delta_x
     real(kind=r8)                              :: delta_y
     real(kind=r8)                              :: nests
     real(kind=r8)                              :: pts_x
     real(kind=r8)                              :: pts_y
     real(kind=r8)                              :: anchor_i
     real(kind=r8)                              :: anchor_j
     real(kind=r8)                              :: iref
     real(kind=r8)                              :: jref
     real(kind=r8), dimension(:),   allocatable :: dsigma
     real(kind=r8), dimension(:),   allocatable :: dsigmaw
     real(kind=r8), dimension(:),   allocatable :: msigma
     real(kind=r8), dimension(:),   allocatable :: wsigma
     real(kind=r8), dimension(:,:), allocatable :: terrain

     ! These fields are used for dealing with domain decomposition
     integer, dimension(:), allocatable :: iminf
     integer, dimension(:), allocatable :: imaxf
     integer, dimension(:), allocatable :: jminf
     integer, dimension(:), allocatable :: jmaxf
  end type coamps_grid

  integer, parameter :: SINGLE_POINT = 1

  ! Model vertical coordinate indexes top-down
  integer, parameter :: TOP_LEVEL    = 1

  ! Location of entries in the datahd file
  integer, parameter :: DATAHD_NUM_P_LEVELS   = 1
  integer, parameter :: DATAHD_NUM_S_LEVELS   = 2
  integer, parameter :: DATAHD_MAP_PROJECTION = 3
  integer, parameter :: DATAHD_STANDARD_LAT_1 = 4
  integer, parameter :: DATAHD_STANDARD_LAT_2 = 5
  integer, parameter :: DATAHD_STANDARD_LON   = 6
  integer, parameter :: DATAHD_REFERENCE_LAT  = 7
  integer, parameter :: DATAHD_REFERENCE_LON  = 8
  integer, parameter :: DATAHD_DELTA_X        = 9
  integer, parameter :: DATAHD_DELTA_Y        = 10
  integer, parameter :: DATAHD_NUM_NESTS      = 11
  integer, parameter :: DATAHD_NUM_X_POINTS   = 30
  integer, parameter :: DATAHD_NUM_Y_POINTS   = 31
  integer, parameter :: DATAHD_ANCHOR_POINT_I = 32
  integer, parameter :: DATAHD_ANCHOR_POINT_J = 33
  integer, parameter :: DATAHD_REFERENCE_I    = 34
  integer, parameter :: DATAHD_REFERENCE_J    = 35
  integer, parameter :: DATAHD_DSIGMA_OFFSET  = 500
  integer, parameter :: DATAHD_MSIGMA_OFFSET  = 800

  !------------------------------
  ! END TYPES AND CONSTANTS
  !------------------------------

  !------------------------------
  ! BEGIN MODULE VARIABLES
  !------------------------------

! version controlled file description for error handling, do not edit
character(len=256), parameter :: source   = &
   "$URL$"
character(len=32 ), parameter :: revision = "$Revision$"
character(len=128), parameter :: revdate  = "$Date$"

  logical :: module_initialized = .false.

  !------------------------------
  ! END MODULE VARIABLES
  !------------------------------

contains

  !------------------------------
  ! BEGIN PUBLIC ROUTINES
  !------------------------------

  ! initialize_grid
  ! ---------------
  ! Populates a coamps_grid data structure using the COAMPS domain
  ! data (datahd) file and the COAMPS terrain height (terrht) file
  ! based on a given date-time group.
  ! N.B. "grid" parameter is marked as intent(inout) to ensure that
  !      grid will not be undefined on exit if we don't need to run
  !      this routine
  !  PARAMETERS
  !   IN  dtg               base date-time group for model run
  ! INOUT grid              coamps_grid structure to initialize
  subroutine initialize_grid(dtg, grid)
    character(len=10), intent(in)    :: dtg
    type(coamps_grid), intent(inout) :: grid

    character(len=64) :: datahd_filename
    integer           :: datahd_unit
    character(len=64) :: terrht_filename
    integer           :: terrht_unit
    integer           :: terrht_record_len

    ! Terrain data is stored in a flat file, so it uses the COAMPS
    ! real number size.  Single dimension makes byteswapping easier
    real(kind=C_REAL), dimension(:), allocatable :: coamps_terrain

    integer :: field_size

    ! Error-checking info
    character(len=*), parameter :: routine = 'initialize_grid'
    integer :: io_status
    integer :: alloc_status
    integer :: dealloc_status

    if (module_initialized) return

    ! Initialize the domain data section of the coamps_grid structure
    call generate_datahd_filename(dtg, datahd_filename)
    datahd_unit = get_unit()
    open(unit=datahd_unit, file=datahd_filename, status='old', &
         access='sequential', action='read', form='formatted', &
         iostat=io_status)
    call check_io_status(io_status, routine, source, revision, &
                         revdate, 'Opening datahd file')
    call read_datahd_file(datahd_unit, grid)
    close(datahd_unit)

    ! Need to set up permanent and temporary storage for terrain data
    allocate( grid%terrain(int(grid%pts_x), int(grid%pts_y)), &
              stat = alloc_status )
    call check_alloc_status(alloc_status, routine, source,    &
                            revision, revdate, 'grid%terrain')

    call get_grid_field_size(grid, field_size)
    allocate (coamps_terrain(field_size), stat = alloc_status)
    call check_alloc_status(alloc_status, routine, source, revision, &
                            revdate, 'coamps_terrain (temporary)')

    ! Make sure when checking the record length to use what *COAMPS*
    ! wrote it as, not what we want to store it as
    call generate_terrht_filename(dtg, grid, terrht_filename)
    terrht_unit = get_unit()
    inquire(iolength=terrht_record_len) coamps_terrain

    open(unit = terrht_unit, file = terrht_filename, status = 'old',&
         access = 'direct', action = 'read', form = 'unformatted',  &
         recl = terrht_record_len, iostat = io_status)
    call check_io_status(io_status, routine, source, revision, &
                         revdate, 'Opening terrain height file')
    call read_terrht_file(terrht_unit, coamps_terrain)
    close(terrht_unit)

    ! COAMPS writes big-endian flat files no matter what platform is
    ! being used, so we may need to change the data we read.
    call fix_for_platform(coamps_terrain, field_size, C_REAL)

    grid%terrain = reshape(real(coamps_terrain, kind=r8),        &
                           (/ int(grid%pts_x), int(grid%pts_y) /))

    deallocate(coamps_terrain, stat=dealloc_status)
    call check_dealloc_status(dealloc_status, routine, source, revision, &
                              revdate, 'coamps_terrain')

    module_initialized = .true.
  end subroutine initialize_grid

  ! gridpt_to_latlon
  ! ----------------
  ! Given a COAMPS grid structure, convert an (i,j) representation
  ! of a point on that grid to a (lat, lon) representation.
  ! Note that this routine will happily accept and convert (i,j)
  ! points that are beyond the bounds of the domain.  For this
  ! direction of conversion, that's not too big a problem.
  !  PARAMETERS
  !   IN  grid              coamps_grid to base conversion on
  !   IN  ii                point's x-direction index
  !   IN  jj                point's y-direction index
  !   OUT lat               point's latitude
  !   OUT lon               point's longitude
  subroutine gridpt_to_latlon(grid, ii, jj, lat, lon)
    type(coamps_grid), intent(in) :: grid
    real(kind=r8), intent(in)     :: ii
    real(kind=r8), intent(in)     :: jj
    real(kind=r8), intent(out)    :: lat
    real(kind=r8), intent(out)    :: lon

    ! Work-around for not being able to pass in scalars as arrays
    ! of length 1
    real(kind=r8), dimension(SINGLE_POINT) :: ii_array,  jj_array
    real(kind=r8), dimension(SINGLE_POINT) :: lat_array, lon_array

    ii_array(SINGLE_POINT) = ii
    jj_array(SINGLE_POINT) = jj

    call ij2ll(int(grid%map_proj), grid%reflat,    grid%reflon,  &
               int(grid%iref),     int(grid%jref), grid%stdlat1, &
               grid%stdlat2,       grid%stdlon,    grid%delta_x, &
               grid%delta_y,       ii_array,       jj_array,     &
               SINGLE_POINT,       lat_array,      lon_array)

    lat = lat_array(SINGLE_POINT)
    lon = lon_array(SINGLE_POINT)
  end subroutine gridpt_to_latlon

  ! latlon_to_gridpt
  ! ----------------
  ! Given a COAMPS grid structure, convert a (lat,lon)
  ! representation of a point on that grid to an (i, j)
  ! representation.
  ! Note that this routine will happily accept and convert
  ! (lat,lon) points that are beyond the bounds of the domain.
  ! This could happen in several ways, but especially if an obs
  ! gets included in the sequence that for some reason isn't in
  ! the domain - in this case, warn the user.
  !  PARAMETERS
  !   IN  grid              coamps_grid to base conversion on
  !   IN  lat               point's latitude
  !   IN  lon               point's longitude
  !   OUT ii                point's x-direction index
  !   OUT jj                point's y-direction index
  subroutine latlon_to_gridpt(grid, lat, lon, ii, jj)
    type(coamps_grid), intent(in) :: grid
    real(kind=r8), intent(in)     :: lat
    real(kind=r8), intent(in)     :: lon
    real(kind=r8), intent(out)    :: ii
    real(kind=r8), intent(out)    :: jj

    ! Work-around for not being able to pass in scalars as arrays
    ! of length 1
    real(kind=r8), dimension(SINGLE_POINT) :: ii_array,  jj_array
    real(kind=r8), dimension(SINGLE_POINT) :: lat_array, lon_array

    ! If anything goes wrong in the conversion, pre-define values
    ! that should set off alarm bells
    ii = -999.0
    jj = -999.0

    lat_array(SINGLE_POINT) = lat
    lon_array(SINGLE_POINT) = lon

    call ll2ij(int(grid%map_proj), grid%reflat,    grid%reflon,  &
               int(grid%iref),     int(grid%jref), grid%stdlat1, &
               grid%stdlat2,       grid%stdlon,    grid%delta_x, &
               grid%delta_y,       ii_array,       jj_array,     &
               SINGLE_POINT,       lat_array,      lon_array)

    ii = ii_array(SINGLE_POINT)
    jj = jj_array(SINGLE_POINT)

    if ( ii .lt. 1 .or. ii .gt. grid%pts_x .or. &
         jj .lt. 1 .or. jj .gt. grid%pts_y) then
        call error_handler(E_WARN, 'latlon_to_gridpt', source, &
                           revision, revdate, 'Point is' //    &
                           'outside domain')
    end if
  end subroutine latlon_to_gridpt

  ! gridpt_to_location
  ! ------------------
  ! Given a COAMPS grid, an (i,j) horizontal location in that grid,
  ! and a vertical location and vertical coordinate type, returns
  ! a DART location structure for that point
  !  PARAMETERS
  !   IN  grid              coamps_grid for lat/lon conversion
  !   IN  ii                point's x-direction index
  !   IN  jj                point's y-direction index
  !   IN  vert_loc          point's vertical direction location
  !   IN  vert_coord        point's vertical direction coordinate
  !   OUT loc               point's location as a DART structure
  subroutine gridpt_to_location(grid, ii, jj, vert_loc, vert_coord,&
                                loc)
    type(coamps_grid), intent(in)    :: grid
    integer, intent(in)              :: ii
    integer, intent(in)              :: jj
    real(kind=r8), intent(in)        :: vert_loc
    integer, intent(in)              :: vert_coord
    type(location_type), intent(out) :: loc

    real(kind=r8) :: real_ii, real_jj
    real(kind=r8) :: point_lat, point_lon

    ! Get things squared away in 2D first
    real_ii = real(ii, kind=r8)
    real_jj = real(jj, kind=r8)
    call gridpt_to_latlon(grid, real_ii, real_jj, point_lat, &
                          point_lon)

    ! Do this with keyword arguments because I can't keep the order
    ! straight
    loc = set_location( lon        = point_lon, &
                        lat        = point_lat, &
                        vert_loc   = vert_loc,  &
                        which_vert = vert_coord )
  end subroutine gridpt_to_location

  ! location_to_gridpt
  ! ------------------
  ! Given a COAMPS grid and a DART location, return the location
  ! of a point as an (i, j) coordinate and a vertical coordinate.
  ! Note that DART's get_location function that's used to break
  ! this down does *not* return any information on what that
  ! vertical coordinate is, so no information is returned to the
  ! caller.
  !  PARAMETERS
  !   IN  grid              coamps_grid for lat/loon conversion
  !   IN  loc               DART location structure to unpack
  !   OUT ii                point's x-direction index
  !   OUT jj                point's y-direction index
  !   OUT vert_loc          point's vertical location
  subroutine location_to_gridpt(grid, loc, ii, jj, vert_loc)
    type(coamps_grid), intent(in)   :: grid
    type(location_type), intent(in) :: loc
    real(kind=r8), intent(out)      :: ii
    real(kind=r8), intent(out)      :: jj
    real(kind=r8), intent(out)      :: vert_loc

    ! The array returned by get_location is (lon, lat, vertical)
    real(kind=r8), dimension(3) :: loc_array

    integer, parameter :: DART_LOC_LON  = 1
    integer, parameter :: DART_LOC_LAT  = 2
    integer, parameter :: DART_LOC_VERT = 3

    loc_array = get_location(loc)

    ! Deal with the 2D component separately from the vertical
    ! Again, be explicit since I apparently like to mix these up
    call latlon_to_gridpt( grid = grid,                    &
                           lat  = loc_array(DART_LOC_LAT), &
                           lon  = loc_array(DART_LOC_LON), &
                           ii   = ii,                      &
                           jj   = jj )
    vert_loc = loc_array(DART_LOC_VERT)
  end subroutine location_to_gridpt

  ! grid_ij_to_vector_index
  ! -----------------------
  ! Given a COAMPS grid, converts an horizontal (ii,jj) location
  ! to an index in the vector representation
  !  PARAMETERS
  !   IN  grid              coamps_grid to use for the conversion
  !   IN  ii                zonal index component
  !   IN  jj                meridional index component
  !   OUT index             position that the point (ii,jj) has in
  !                         the field re-arranged in vector form
  subroutine grid_ij_to_vector_index(grid, ii, jj, index)
    type(coamps_grid), intent(in) :: grid
    integer, intent(in)           :: ii
    integer, intent(in)           :: jj
    integer, intent(out)          :: index

    index = (jj - 1) * grid%pts_x + ii
  end subroutine grid_ij_to_vector_index

  ! get_grid_dims
  ! -------------
  ! Given a COAMPS grid structure, calculate the i (i.e. x) and
  ! j (i.e. y) dimensions
  !  PARAMETERS
  !   IN  grid              coamps_grid to get limits from
  !   OUT grid_i            number of points on the i dimension
  !   OUT grid_j            number of points on the j dimension
  subroutine get_grid_dims(grid, grid_i, grid_j)
    type(coamps_grid), intent(in) :: grid
    integer, intent(out)          :: grid_i
    integer, intent(out)          :: grid_j

    grid_i = int(grid%pts_x)
    grid_j = int(grid%pts_y)
  end subroutine get_grid_dims

  ! get_grid_field_size
  ! -------------------
  ! Given a coamps grid structure, calculate the field size (the
  ! total number of points in the i and j directions)
  !  PARAMETERS
  !   IN  grid              coamps_grid to get size from
  !   OUT field_size        number of grid points on a single level
  subroutine get_grid_field_size(grid, field_size)
    type(coamps_grid), intent(in) :: grid
    integer, intent(out)          :: field_size

    field_size = int(grid%pts_x) * int(grid%pts_y)
  end subroutine get_grid_field_size

  ! get_grid_num_levels
  ! -------------------
  ! Returns the number of mass sigma levels on a given COAMPS grid
  !  PARAMETERS
  !   IN  grid              coamps_grid to get levels from
  !   OUT num_sig_lvls      number of sigma levels on the grid
  subroutine get_grid_num_levels(grid, num_sig_lvls)
    type(coamps_grid), intent(in) :: grid
    integer, intent(out)          :: num_sig_lvls

    num_sig_lvls = grid%sigm_lvls
  end subroutine get_grid_num_levels

  ! check_ij_within_grid
  ! --------------------
  ! Calculates if a given (i,j) index is within a given COAMPS grid
  !  PARAMETERS
  !   IN  grid              coamps_grid to use for checking
  !   IN  ii                zonal index
  !   IN  jj                meridional index
  !   OUT is_within_grid    True if the (i,j) location listed is
  !                         inside the given grid
  subroutine check_ij_within_grid(grid, ii, jj, is_within_grid)
    type(coamps_grid), intent(in) :: grid
    integer, intent(in)           :: ii
    integer, intent(in)           :: jj
    logical, intent(out)          :: is_within_grid

    is_within_grid = .true.
    if ( (ii < 1) .or. (ii > grid%pts_x) .or. &
         (jj < 1) .or. (jj > grid%pts_y) ) then
       is_within_grid = .false.
    end if
  end subroutine check_ij_within_grid

  ! get_grid_delta_x
  ! ----------------
  ! Returns the x-direction grid spacing (in meters) for the given grid
  !  PARAMETERS
  !   IN  grid              coamps_grid to pull data from
  !   OUT grid_delta_x      x-direction grid spacing
  subroutine get_grid_delta_x(grid, grid_delta_x)
    type(coamps_grid), intent(in)  :: grid
    real(kind=r8),     intent(out) :: grid_delta_x

    grid_delta_x = grid%delta_x
  end subroutine get_grid_delta_x

  ! get_grid_delta_y
  ! ----------------
  ! Returns the y-direction grid spacing (in meters) for the given grid
  !  PARAMETERS
  !   IN  grid              coamps_grid to pull data from
  !   OUT grid_delta_y      y-direction grid spacing
  subroutine get_grid_delta_y(grid, grid_delta_y)
    type(coamps_grid), intent(in)  :: grid
    real(kind=r8),     intent(out) :: grid_delta_y

    grid_delta_y = grid%delta_y
  end subroutine get_grid_delta_y

  ! get_grid_msigma
  ! ---------------
  ! Returns the mass sigma levels for a given COAMPS grid
  !  PARAMETERS
  !   IN  grid              coamps_grid to pull msigma out of
  !   OUT msigma            mass sigma levels for the grid
  subroutine get_grid_msigma(grid, msigma)
    type(coamps_grid), intent(in)            :: grid
    real(kind=r8), dimension(:), intent(out) :: msigma

    integer :: kka

    ! Don't attempt a copy if the data in the grid won't fit
    if ( size(grid%msigma) > size(msigma) ) then
       call error_handler(E_ERR, 'get_grid_msigma', 'msigma    &
                          &passed into function is too small', &
                          &source, revision, revdate)
    end if

    do kka=1,size(grid%msigma)
       msigma(kka) = grid%msigma(kka)
    end do
  end subroutine get_grid_msigma

  ! get_grid_wsigma
  ! ---------------
  ! Returns the w sigma levels for a given COAMPS grid
  !  PARAMETERS
  !   IN  grid              coamps_grid to pull wsigma out of
  !   OUT wsigma            w sigma levels for the grid
  subroutine get_grid_wsigma(grid, wsigma)
    type(coamps_grid), intent(in)            :: grid
    real(kind=r8), dimension(:), intent(out) :: wsigma

    integer :: kkw

    ! Bail if the data in the grid won't fit
    if ( size(grid%wsigma) > size(wsigma) ) then
       call error_handler(E_ERR, 'get_grid_wsigma', 'wsigma ' // &
                          'passed into function is too small',   &
                          source, revision, revdate)
    end if

    do kkw=1,size(grid%wsigma)
       wsigma(kkw) = grid%wsigma(kkw)
    end do
  end subroutine get_grid_wsigma

  ! get_grid_dsigmaw
  ! ----------------
  ! Returns the distance between w sigma levels for a given COAMPS
  ! grid
  !  PARAMETERS
  !   IN  grid              coamps_grid to pull dsigmaw out of
  !   OUT dsigmaw           distance between w sigma levels
  subroutine get_grid_dsigmaw(grid, dsigmaw)
    type(coamps_grid), intent(in)            :: grid
    real(kind=r8), dimension(:), intent(out) :: dsigmaw

    integer :: kkw

    ! Bail if the data in the grid won't fit
    if ( size(grid%dsigmaw) > size(dsigmaw) ) then
       call error_handler(E_ERR, 'get_grid_dsigmaw', 'dsigmaw ' // &
                          'passed into function is too small',     &
                          source, revision, revdate)
    end if

    do kkw=1,size(grid%dsigmaw)
       dsigmaw(kkw) = grid%dsigmaw(kkw)
    end do
  end subroutine get_grid_dsigmaw

  ! get_terrain_height_at_points
  ! ----------------------------
  ! Given a COAMPS grid and a vector representation of a set of
  ! (i,j) points, returns the terrain height for those points
  !  PARAMETERS
  !   IN  grid              coamps_grid to pull terrht from
  !   IN  i_points          points's x-direction index
  !   IN  j_points          points's y-direction index
  !   OUT terrain_height    points's terrain height
  subroutine get_terrain_height_at_points(grid, i_points, j_points,&
                                          terrain_height)
    type(coamps_grid), intent(in)            :: grid
    integer, dimension(:), intent(in)        :: i_points
    integer, dimension(:), intent(in)        :: j_points
    real(kind=r8), dimension(:), intent(out) :: terrain_height

    integer :: nn, num_points

    ! Sanity check
    if (size(i_points) .ne. size(j_points)) then
      call error_handler(E_ERR, 'get_terrain_height_at_points',   &
                         'Error in getting terrain height - ' //  &
                         'i/j point arrays are different sizes!', &
                         source, revision, revdate)
    end if

    num_points = size(i_points)
    do nn=1,num_points
       terrain_height(nn) = grid%terrain(i_points(nn), j_points(nn))
    end do
  end subroutine get_terrain_height_at_points

  ! initialize_decomposition
  ! ------------------------
  ! Given a COAMPS grid the total number of non-I/O processors for
  ! that grid, allocate space for each processor's [ij](min|max)f
  ! extent
  ! N.B. The grid is not defined as intent(out) here since that
  !      could mean that the other components of the coamps_grid
  !      are undefined after this routine call.  This is avoided
  !      by using intent(inout).
  !  PARAMETERS
  ! INOUT grid              coamps_grid to allocate space in
  !   IN  totalprocs        total number of non-I/O processors
  subroutine initialize_decomposition(grid, totalprocs)
    type(coamps_grid), intent(inout) :: grid
    integer, intent(in)              :: totalprocs

    character(len=*), parameter :: routine = 'initialize_decomposition'
    integer :: alloc_status

    allocate( grid%iminf(totalprocs), stat=alloc_status )
    call check_alloc_status(alloc_status, routine, source, revision, &
                            revdate, 'grid%iminf')
    allocate( grid%imaxf(totalprocs), stat=alloc_status )
    call check_alloc_status(alloc_status, routine, source, revision, &
                            revdate, 'grid%imaxf')
    allocate( grid%jminf(totalprocs), stat=alloc_status )
    call check_alloc_status(alloc_status, routine, source, revision, &
                            revdate, 'grid%jminf')
    allocate( grid%jmaxf(totalprocs), stat=alloc_status )
    call check_alloc_status(alloc_status, routine, source, revision, &
                            revdate, 'grid%jmaxf')
  end subroutine initialize_decomposition

  ! decompose_domain
  ! ----------------
  ! Given the grid, the number of I/O processors, the number of
  ! processors in x and y, and the number of boundary points (nbnam
  ! in the COAMPS namelist), mimics the domain decomposition
  ! performed by COAMPS routine domdec.F (the math is taken from that
  ! module and variables renamed for use here).  The maximum and
  ! minimum [ij] points for each processor are given in arrays
  ! [ij](min|max)f within the coamps_grid structure.  This routine
  ! populates those arrays.
  !  PARAMETERS
  ! INOUT grid              coamps_grid to decompose
  !   IN  num_io_procs      number of I/O processors (0 or 1)
  !   IN  dom_proc_x        number of processors in the x direction
  !   IN  dom_proc_y        number of processors in the y direction
  !   IN  nbound            number of halo boundary points
  subroutine decompose_domain(grid, num_io_procs, dom_proc_x, &
                              dom_proc_y, nbound)
    type(coamps_grid), intent(inout) :: grid
    integer, intent(in)              :: num_io_procs
    integer, intent(in)              :: dom_proc_x
    integer, intent(in)              :: dom_proc_y
    integer, intent(in)              :: nbound

    ! Total number of processors we are iterating over, along with
    ! the current processor and some counters
    integer :: nproc_x, nproc_y
    integer :: totalprocs, curproc, x_count, y_count

    ! Number of points per processor (not including the boundary)
    integer :: xperproc, yperproc
    integer :: xextra, yextra

    ! The initial boundaries for the domain (not counting the
    ! halo points)
    integer, allocatable, dimension(:) :: begin_i
    integer, allocatable, dimension(:) :: begin_j
    integer, allocatable, dimension(:) :: end_i
    integer, allocatable, dimension(:) :: end_j

    ! Error handling
    character(len=*), parameter :: routine = 'decompose_domain'
    integer :: alloc_status, dealloc_status

    ! If we are only doing single processor I/O, it doesn't matter
    ! how many processors there are - we treat it as a single
    ! decomposition
    if (num_io_procs .gt. 0) then
       nproc_x = 1
       nproc_y = 1
    else
       nproc_x = dom_proc_x
       nproc_y = dom_proc_y
    endif
    totalprocs = nproc_x * nproc_y

    ! This sequence of allocations could potentially be offloaded to
    ! the initialize_decomposition subroutine but as the arrays in
    ! question are local and temporary, I've kept them here.
    allocate( begin_i(totalprocs), stat=alloc_status )
    call check_alloc_status(alloc_status, routine, source, revision, &
                            revdate, 'begin_i')
    allocate( begin_j(totalprocs), stat=alloc_status )
    call check_alloc_status(alloc_status, routine, source, revision, &
                            revdate, 'begin_j')
    allocate( end_i(totalprocs), stat=alloc_status )
    call check_alloc_status(alloc_status, routine, source, revision, &
                            revdate, 'end_i')
    allocate( end_j(totalprocs), stat=alloc_status )
    call check_alloc_status(alloc_status, routine, source, revision, &
                            revdate, 'end_j')

    ! Get the processor-point counts
    xperproc = (grid%pts_x - 2) / nproc_x
    yperproc = (grid%pts_y - 2) / nproc_y
    xextra   = (grid%pts_x - 2) - (xperproc * nproc_x)
    yextra   = (grid%pts_y - 2) - (yperproc * nproc_y)

    do y_count = 1,nproc_y
       do x_count = 1,nproc_x

          ! 2D grid -> 1D number
          curproc = (y_count - 1) * nproc_x + x_count

          ! X points
          if (x_count .eq. 1) then
             begin_i(curproc) = 2
          else
             begin_i(curproc) = end_i(curproc-1) + 1
          endif

          end_i(curproc) = begin_i(curproc) + xperproc - 1
          if (x_count .gt. nproc_x-xextra) then
             end_i(curproc) = end_i(curproc) + 1
          endif

          ! Y Points
          if (y_count .eq. 1) then
             begin_j(curproc) = 2
          else
             begin_j(curproc) = end_j(curproc-nproc_x) + 1
          endif
          end_j(curproc) = begin_j(curproc) + yperproc - 1
          if (y_count .gt. nproc_y-yextra) then
             end_j(curproc) = end_j(curproc) + 1
          endif

          ! Take the boundaries into consideration
          grid%iminf(curproc) = begin_i(curproc) - nbound
          grid%imaxf(curproc) = end_i(curproc)   + nbound
          grid%jminf(curproc) = begin_j(curproc) - nbound
          grid%jmaxf(curproc) = end_j(curproc)   + nbound
       enddo
    enddo

    deallocate( begin_i, begin_j, end_i, end_j, stat=dealloc_status )
    call check_dealloc_status(dealloc_status, routine, source, &
                              revision, revdate, '[ij](min|max)f')
  end subroutine decompose_domain

  ! get_grid_iminf
  ! --------------
  ! Accessor method for the index of a subgrid's left extent
  !  PARAMETERS
  !   IN  grid              coamps_grid to get data from
  !   IN  subgrid_num       decomposed subgrid number
  !   OUT iminf             i-index in large domain for the
  !                         left edge of this subgrid
  subroutine get_grid_iminf(grid, subgrid_num, iminf)
    type(coamps_grid), intent(in)  :: grid
    integer,           intent(in)  :: subgrid_num
    integer,           intent(out) :: iminf

    iminf = grid%iminf(subgrid_num)
  end subroutine

  ! get_grid_imaxf
  ! --------------
  ! Accessor method for the index of a subgrid's right extent
  !  PARAMETERS
  !   IN  grid              coamps_grid to get data from
  !   IN  subgrid_num       decomposed subgrid number
  !   OUT imaxf             i-index in large domain for the
  !                         right edge of this subgrid
  subroutine get_grid_imaxf(grid, subgrid_num, imaxf)
    type(coamps_grid), intent(in)  :: grid
    integer,           intent(in)  :: subgrid_num
    integer,           intent(out) :: imaxf

    imaxf = grid%imaxf(subgrid_num)
  end subroutine

  ! get_grid_jminf
  ! --------------
  ! Accessor method for the index of a subgrid's lower extent
  !  PARAMETERS
  !   IN  grid              coamps_grid to get data from
  !   IN  subgrid_num       decomposed subgrid number
  !   OUT jminf             j-index in large domain for the
  !                         lower edge of this subgrid
  subroutine get_grid_jminf(grid, subgrid_num, jminf)
    type(coamps_grid), intent(in)  :: grid
    integer,           intent(in)  :: subgrid_num
    integer,           intent(out) :: jminf

    jminf = grid%jminf(subgrid_num)
  end subroutine

  ! get_grid_jmaxf
  ! --------------
  ! Accessor method for the index of a subgrid's upper extent
  !  PARAMETERS
  !   IN  grid              coamps_grid to get data from
  !   IN  subgrid_num       decomposed subgrid number
  !   OUT jmaxf             j-index in large domain for the
  !                         upper edge of this subgrid
  subroutine get_grid_jmaxf(grid, subgrid_num, jmaxf)
    type(coamps_grid), intent(in)  :: grid
    integer,           intent(in)  :: subgrid_num
    integer,           intent(out) :: jmaxf

    jmaxf = grid%jmaxf(subgrid_num)
  end subroutine

  ! dump_grid_info
  ! --------------
  ! Dumps a human-readable summary of the supplied grid
  !  PARAMETERS
  !   IN  grid              coamps_grid to print out
  subroutine dump_grid_info(grid)
    type(coamps_grid), intent(in) :: grid

    if (do_output()) then
       write (*,*) "**** COAMPS GRID INFORMATION ****"
       write (*,*) "Number of pressure levels:", grid%pres_lvls
       write (*,*) "Number of sigma levels:   ", grid%sigm_lvls
       write (*,*) "Map projection type:      ", grid%map_proj
       write (*,*) "Latitude intersection 1:  ", grid%stdlat1
       write (*,*) "Latitude intersection 2:  ", grid%stdlat2
       write (*,*) "Longitude intersection:   ", grid%stdlon
       write (*,*) "Reference latitude:       ", grid%reflat
       write (*,*) "Reference longitude:      ", grid%reflon
       write (*,*) "X-direction spacing (m):  ", grid%delta_x
       write (*,*) "Y-direction spacing (m):  ", grid%delta_y
       write (*,*) "Number of nests:          ", grid%nests
       write (*,*) "Number of X grid points:  ", grid%pts_x
       write (*,*) "Number of Y grid points:  ", grid%pts_y
       write (*,*) "i point on parent grid:   ", grid%anchor_i
       write (*,*) "j point on parent grid:   ", grid%anchor_j
       write (*,*) "i point at ref. lat/lon:  ", grid%iref
       write (*,*) "j point at ref. lat/lon:  ", grid%jref
    end if
  end subroutine dump_grid_info

  !------------------------------
  ! END PUBLIC ROUTINES
  !------------------------------

  !------------------------------
  ! BEGIN PRIVATE ROUTINES
  !------------------------------

  ! generate_datahd_filename
  ! ------------------------
  ! Generates the COAMPS domain information file name for the first
  ! nest at a given date-time group.  This does *not* generate any
  ! path information - it only returns the file name.
  !  PARAMETERS
  !   IN  dtg               base date-time group for model run
  !   OUT datahd_filename   name of domain information file
  subroutine generate_datahd_filename(dtg, datahd_filename)
    character(len=10), intent(in)  :: dtg
    character(len=64), intent(out) :: datahd_filename

    ! The format of the datahd filename is fixed except for the date
    ! -time group: mimics a 2000x1 flat file with no level
    ! information or tau information
    call generate_flat_file_name( var_name   = 'datahd',      &
                                  level_type = 'sfc',         &
                                  level1     = 0,             &
                                  level2     = 0,             &
                                  gridnum    = 1,             &
                                  aoflag     = 'a',           &
                                  xpts       = 2000,          &
                                  ypts       = 1,             &
                                  dtg        = dtg,           &
                                  tau_dd     = 0,             &
                                  tau_hh     = 0,             &
                                  tau_mm     = 0,             &
                                  tau_ss     = 0,             &
                                  field_type = 'infofld',     &
                                  file_name  = datahd_filename )
  end subroutine generate_datahd_filename

  ! generate_terrht_filename
  ! ------------------------
  ! Generates the COAMPS terrain height flat file name for the first
  ! nest of a grid at a given date-time group.  This does *not*
  ! generate any path information - it only returns the file name.
  !
  ! Assumes that the coamps_grid structure has already been populated
  ! with the domain information, since we need the number of x and y
  ! grid points to generate the file name
  !  PARAMETERS
  !   IN  dtg               base date-time group for model run
  !   IN  grid              coamps_grid that we need terrain info for
  !   OUT terrht_filename   name of terrain height flat file
  subroutine generate_terrht_filename(dtg, grid, terrht_filename)
    character(len=10), intent(in)  :: dtg
    type(coamps_grid), intent(in)  :: grid
    character(len=64), intent(out) :: terrht_filename

    integer :: x_points
    integer :: y_points

    x_points = int(grid%pts_x)
    y_points = int(grid%pts_y)

    ! The format of the terrain height file is fixed except for
    ! the horizontal grid size and date-time group
    call generate_flat_file_name( var_name   = 'terrht',      &
                                  level_type = 'sfc',         &
                                  level1     = 0,             &
                                  level2     = 0,             &
                                  gridnum    = 1,             &
                                  aoflag     = 'a',           &
                                  xpts       = x_points,      &
                                  ypts       = y_points,      &
                                  dtg        = dtg,           &
                                  tau_dd     = 0,             &
                                  tau_hh     = 0,             &
                                  tau_mm     = 0,             &
                                  tau_ss     = 0,             &
                                  field_type = 'fcstfld',     &
                                  file_name  = terrht_filename )
  end subroutine generate_terrht_filename

  ! generate_flat_file_name
  ! -----------------------
  ! Given field, level, and grid information, generate the properly
  ! formatted 64-character COAMPS flat file name.  Note that this
  ! does *not* generate any path information - it only returns the
  ! file name.
  !  PARAMETERS
  !   IN  var_name          the field the file contains
  !   IN  level_type        vertical level type (height/pressure/etc)
  !   IN  level1            lowest vertical level in the file
  !   IN  level2            highest vertical level in the file
  !                         (for files for a single level, level1 is
  !                          that level and level2 is left to 0)
  !   IN  gridnum           nest number (only 1 supported for now)
  !   IN  aoflag            field type: (a)tmosphere or (o)cean
  !   IN  xpts              number of points in the x direction
  !   IN  ypts              number of points in the y direction
  !   IN  dtg               base date-time group
  !   IN  tau_dd            forecast lead time - day component
  !   IN  tau_hh            forecast lead time - hour component
  !   IN  tau_mm            forecast lead time - minute component
  !   IN  tau_ss            forecast lead time - second component
  !   IN  field_type        type of field (e.g. fcstfld, infofld)
  !   OUT file_name         COAMPS flat file name
  subroutine generate_flat_file_name(var_name, level_type, level1, &
                                     level2, gridnum, aoflag, xpts,&
                                     ypts, dtg, tau_dd, tau_hh,    &
                                     tau_mm, tau_ss, field_type,   &
                                     file_name)
    character(len=6), intent(in)   :: var_name
    character(len=3), intent(in)   :: level_type
    integer, intent(in)            :: level1
    integer, intent(in)            :: level2
    integer, intent(in)            :: gridnum
    character(len=1), intent(in)   :: aoflag
    integer, intent(in)            :: xpts
    integer, intent(in)            :: ypts
    character(len=10), intent(in)  :: dtg
    integer, intent(in)            :: tau_dd
    integer, intent(in)            :: tau_hh
    integer, intent(in)            :: tau_mm
    integer, intent(in)            :: tau_ss
    character(len=7), intent(in)   :: field_type
    character(len=64), intent(out) :: file_name

    write(file_name, 100) var_name, level_type, level1, level2,&
         & gridnum, aoflag, xpts, ypts, dtg, tau_dd, tau_hh, tau_mm,&
         & tau_ss, field_type

100 format(A6'_'A3'_'I6.6'_'I6.6'_'I1A1I4.4'x'I4.4'_'A10'_'I2.2I2.2&
         I2.2I2.2'_'A7)
  end subroutine generate_flat_file_name

  ! read_datahd_file
  ! ----------------
  ! Given the unit number of an *open* COAMPS domain information file,
  ! populate the portion of a COAMPS grid structure that corresponds
  ! to the data contained in the file
  !  PARAMETERS
  !   IN  datahd_unit       unit number of open datahd file
  !   OUT grid              coamps_grid structure to be filled
  subroutine read_datahd_file(datahd_unit, grid)
    integer, intent(in)            :: datahd_unit
    type(coamps_grid), intent(out) :: grid

    ! Domain information from file
    real(kind=r8), dimension(2000) :: domain_data

    ! Model sigma levels
    integer :: kka
    integer :: cur_lev

    ! True if a the datahd file is opened & if memory's allocated
    logical :: is_opened, is_allocated

    ! Error checking
    character(len=*), parameter :: routine = 'read_datahd_file'
    integer :: io_status, alloc_status

    ! This is really an "assert" statement rather than error capture
    ! since the user should not have to worry about this
    inquire(unit=datahd_unit, opened=is_opened)
    if ( .not. is_opened ) then
      call error_handler(E_ERR,'read_datahd_file','Trying to read a &
                         &datahd file that is not open yet', source,&
                         &revision, revdate)
    end if

    ! Probably shouldn't do this as free-form I/O, but it works so
    ! stick with it for now - maybe add a real format statement later.
    read(unit=datahd_unit, fmt=*,iostat=io_status) domain_data
    call check_io_status(io_status, routine, source, revision, &
                         revdate, 'Reading datahd file')

    ! Populate the grid structure - do the scalar variables first
    ! This information can be found around line 875 in the COAMPS
    ! atmos/libsrc/aalib/coama.F file
    grid%pres_lvls = domain_data(DATAHD_NUM_P_LEVELS)
    grid%sigm_lvls = domain_data(DATAHD_NUM_S_LEVELS)
    grid%map_proj  = domain_data(DATAHD_MAP_PROJECTION)
    grid%stdlat1   = domain_data(DATAHD_STANDARD_LAT_1)
    grid%stdlat2   = domain_data(DATAHD_STANDARD_LAT_2)
    grid%stdlon    = domain_data(DATAHD_STANDARD_LON)
    grid%reflat    = domain_data(DATAHD_REFERENCE_LAT)
    grid%reflon    = domain_data(DATAHD_REFERENCE_LON)
    grid%delta_x   = domain_data(DATAHD_DELTA_X)
    grid%delta_y   = domain_data(DATAHD_DELTA_Y)
    grid%nests     = domain_data(DATAHD_NUM_NESTS)

    ! Individual nest information - currently only 1 next so
    ! explicitly specify indices.  This would be a loop later if I
    ! allowed nested grids.
    grid%pts_x    = domain_data(DATAHD_NUM_X_POINTS)
    grid%pts_y    = domain_data(DATAHD_NUM_Y_POINTS)
    grid%anchor_i = domain_data(DATAHD_ANCHOR_POINT_I)
    grid%anchor_j = domain_data(DATAHD_ANCHOR_POINT_J)
    grid%iref     = domain_data(DATAHD_REFERENCE_I)
    grid%jref     = domain_data(DATAHD_REFERENCE_J)

    ! Again, this is an assert more than an error trap
    is_allocated = allocated(grid%dsigma)  .and. &
                   allocated(grid%dsigmaw) .and. &
                   allocated(grid%msigma)  .and. &
                   allocated(grid%wsigma)
    if (is_allocated) then
      call error_handler(E_ERR, routine, 'Grid sigma arrays are ' // &
                         'already allocated', source, revision, revdate)
    end if

    kka = grid%sigm_lvls
    allocate( grid%dsigma(kka),      stat=alloc_status)
    call check_alloc_status(alloc_status, routine, source, revision, &
                            revdate, 'grid%dsigma')
    allocate( grid%dsigmaw(kka + 1), stat=alloc_status)
    call check_alloc_status(alloc_status, routine, source, revision, &
                            revdate, 'grid%dsigmaw')
    allocate( grid%msigma(kka),      stat=alloc_status)
    call check_alloc_status(alloc_status, routine, source, revision, &
                            revdate, 'grid%msigma')
    allocate( grid%wsigma(kka + 1),  stat=alloc_status)
    call check_alloc_status(alloc_status, routine, source, revision, &
                            revdate, 'grid%wsigma')

    ! Now that we have the space, information about mass sigma
    ! levels can just be read in from the file...
    do cur_lev=1,kka
       grid%dsigma(cur_lev) = domain_data(DATAHD_DSIGMA_OFFSET + cur_lev)
       grid%msigma(cur_lev) = domain_data(DATAHD_MSIGMA_OFFSET + cur_lev)
    end do

    ! ... but we have to calculate the w sigma level information
    ! (this is based on the COAMPS utility package source code)
    grid%dsigmaw(TOP_LEVEL) = grid%dsigma(TOP_LEVEL)   * 0.5
    grid%dsigmaw(kka + 1)   = grid%dsigma(kka)         * 0.5
    do cur_lev = 1, (kka - 1)
       grid%dsigmaw(cur_lev+1) = grid%msigma(cur_lev) - grid%msigma(cur_lev+1)
    end do
    grid%wsigma(kka + 1) = 0.0
    do cur_lev=1,kka
       grid%wsigma(kka+1 - cur_lev) = grid%wsigma(kka+1 - cur_lev + 1) + &
                                      grid%dsigma(kka+1 - cur_lev)
    end do
  end subroutine read_datahd_file

  ! read_terrht_file
  ! ----------------
  ! Given the unit number of an *open* COAMPS terrain height flat
  ! file, read the terrain heights into an array.
  !  PARAMETERS
  !   IN  terrht_unit       unit number of an open terrht file
  !   OUT terrain_array     coamps_grid structure to be filled
  subroutine read_terrht_file(terrht_unit, terrain_array)
    integer, intent(in)                          :: terrht_unit
    real(kind=C_REAL), dimension(:), intent(out) :: terrain_array

    character(len=*), parameter :: routine = 'read_terrht_file'
    integer :: io_status

    ! Read in the data - COAMPS writes this out as size C_REAL
    read(unit=terrht_unit, iostat=io_status) terrain_array
    call check_io_status(io_status, routine, source, revision, &
                         revdate, 'Reading terrain height')
  end subroutine read_terrht_file

  !------------------------------
  ! END PRIVATE ROUTINES
  !------------------------------

end module coamps_grid_mod

! <next few lines under version control, do not edit>
! $URL$
! $Id$
! $Revision$
! $Date$

max_stars_repo_name: hkershaw-brown/feature-preprocess

max_issues_repo_name: hkershaw-brown/feature-preprocess

max_forks_repo_name: hkershaw-brown/feature-preprocess


text:
	SUBROUTINE GH_SAVE ( iret )
C************************************************************************
C* GH_SAVE								*
C*									*
C* This subroutine stores the current settings for lines, color, 	*
C* text and special markers. 						*
C*									*
C* GH_SAVE ( IRET )  							*
C*									*
C* Input parameters:							*
C*									*
C* Output parameters:							*
C*	IRET		INTEGER		Return code			*
C*									*
C**									*
C* Log:									*
C* A. Hardy/GSC		 6/01   					*
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
        INCLUDE		'ghcmn.cmn'
C-----------------------------------------------------------------------
	iret = 0
C
C*	Query current settings.
C
	CALL GQSPCL ( szghspcl, jghspwd, ier )
        CALL GQLINE ( jghltyp, jghlthw, jghwdth, jghwhw, ier )
	CALL GQTEXT ( jghtxfn, jghtxhw, ghsztx, jghtxwd, jghbrdr,
     +		      jghrotn, jghjust, ier )
C*
	RETURN
	END

ext: f

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
interface

   double precision function ppnd (p, ifault)
     implicit none
     double precision p
     integer ifault
   end function ppnd

   subroutine normp(x, p, q, pdf)
     implicit none
     double precision x,p,q,pdf
   end subroutine normp

   subroutine sortrx(n,data,index)
     implicit none
     integer n
     double precision, intent(in) :: data(n)
     integer, intent(out) :: index(n)
   end subroutine sortrx

   double precision function betain(x, p, q, beta, ifault)
     implicit none
     double precision x,p,q,beta,alpha
     integer ifault
   end function betain

   double precision function xinbta(p,q,beta,alpha,ifault)
     implicit none
     double precision p,q,beta,alpha
     integer ifault
   end function xinbta

   double precision function alngam(xvalue, ifault)
     implicit none
     double precision xvalue
     integer ifault
   end function alngam

end interface

max_forks_repo_name: GarciaLab/mRNADynamics


text:
      SUBROUTINE HMDOT(J,IMERGE,M1,KW,MC,DMS,RNEW,ITERM)
*
*
*       Mass loss from inner hierarchical binary.
*       -----------------------------------------
*
      INCLUDE 'common6.h'
      COMMON/BINARY/  CM(4,MMAX),XREL(3,MMAX),VREL(3,MMAX),
     &                HM(MMAX),UM(4,MMAX),UMDOT(4,MMAX),TMDIS(MMAX),
     &                NAMEM(MMAX),NAMEG(MMAX),KSTARM(MMAX),IFLAGM(MMAX)
      REAL*8  BODYI(2),W(2)
      REAL*8  M1,MC
      LOGICAL IQUIT
*
*
*       Define merger termination index and relevant KS/c.m. indices.
      ITERM = 0
      IPAIR = KSPAIR
      I1 = 2*IPAIR - 1
      ICM = N + IPAIR
*
*       Include quitting conditions to avoid large changes at end-point.
      IQUIT = .FALSE.
*     IF (RADIUS(J)*SU.GT.500.0) IQUIT = .TRUE.
*     IF (BODY0(J)*ZMBAR.GT.15.0.AND.KSTAR(J).GE.4) IQUIT = .TRUE.
*     IF (ABS(M1 - MC).LT.0.1*M1.OR.MC.GT.1.4) IQUIT = .TRUE.
*     IF (BODY0(J)*ZMBAR - M1.GT.0.5*M1) IQUIT = .TRUE.
*       Note that tidal dissipation needs full KS for rectification in MDOT.
*     IF (KSTARM(IMERGE).LT.0.OR.BODY(ICM).EQ.0.0D0) IQUIT = .TRUE.
      IF (BODY(ICM).EQ.0.0D0) IQUIT = .TRUE.
*
*       Quit for mis-identification, advanced evolution or double merger.
      IF (J.LE.0.OR.IQUIT) THEN
          ITERM = 1
          GO TO 50
      END IF
*
*       Update radius for #J and copy #I1 for MDOT & HCORR in case of ghost.
      RSTAR = RNEW
      IF (J.GT.I1) THEN
          RADIUS(J) = RNEW
          RSTAR = RADIUS(I1)
*       Return RNEW for copying back to RADIUS(I1) on ghost & DM/M < 0.01.
      END IF
*
*       Skip further modifications on zero mass loss of same type.
      IF (ABS(DMS/M1).EQ.0.0D0.AND.KW.EQ.KSTAR(J)) GO TO 50
*
*       Set two-body elements for outer orbit and inner semi-major axis.
      SEMI = -0.5*BODY(ICM)/H(IPAIR)
      ECC2 = (1.0 - R(IPAIR)/SEMI)**2 + TDOT2(IPAIR)**2/(BODY(ICM)*SEMI)
      ECC1 = SQRT(ECC2)
      ZMB0 = CM(1,IMERGE) + CM(2,IMERGE)
      SEMI0 = -0.5*ZMB0/HM(IMERGE)
*
*       Form modified mass ratio and semi-major axes from M*A = const.
      DM = DMS/ZMBAR
      ZMB = ZMB0 - DM
      SEMI1 = SEMI*(ZMB0 + BODY(2*IPAIR))/(ZMB + BODY(2*IPAIR))
      SEMI2 = SEMI0*ZMB0/ZMB
      PMIN = SEMI1*(1.0 - ECC1)
*
*       Evaluate old separation, square regularized velocity, t'' & ECC.
      RI = 0.0D0
      V20 = 0.0
      TD2 = 0.0
      DO 10 K = 1,4
          RI = RI + UM(K,IMERGE)**2
          V20 = V20 + UMDOT(K,IMERGE)**2
          TD2 = TD2 + 2.0*UM(K,IMERGE)*UMDOT(K,IMERGE)
   10 CONTINUE
      ECC2 = (1.0 - RI/SEMI0)**2 + TD2**2/(SEMI0*ZMB0)
      ECC = SQRT(ECC2)
*
*       Determine inclination (use #I1 as first KS component).
      CALL HIMAX(I1,IMERGE,ECC,SEMI0,EMAX,EMIN,ZI,TG,EDAV)
*
*       Evaluate the general stability function (Mardling 2008).
      IF (ECC1.LT.1.0) THEN
          EOUT = ECC1
*       Increase tolerance near sensitive stability boundary (RM 10/2008).
          IF (EOUT.GT.0.9) THEN
              DE = 0.5*(1.0 - EOUT)
              DE = MIN(DE,0.01D0)
*       Add extra amount 0.011 to avoid switching.
              IF (ECC1.GT.0.9) DE = DE + 0.011
              EOUT = EOUT - DE
              PMIN = SEMI1*(1.0 - EOUT)
          END IF
          NST = NSTAB(SEMI2,SEMI1,ECC,EOUT,ZI,CM(1,IMERGE),
     &                                     CM(2,IMERGE),BODY(2*IPAIR))
          IF (NST.EQ.0) THEN
              PCRIT = 0.99*PMIN
              PCR = stability(CM(1,IMERGE),CM(2,IMERGE),BODY(2*IPAIR),
     &                                          ECC,ECC1,ZI)*SEMI2
          ELSE
              PCRIT = 1.01*PMIN
          END IF
      ELSE
          PCRIT = stability(CM(1,IMERGE),CM(2,IMERGE),BODY(2*IPAIR),
     &                                          ECC,ECC1,ZI)*SEMI2
      END IF
*
*       Update pericentre distance on successful stability test or exit.
      IF (PMIN.GT.PCRIT) THEN
          R0(IPAIR) = PCRIT
      ELSE
          ITERM = 1
          GO TO 50
      END IF
*
*       Check for possible circularization (tidal or sequential).
      RP = SEMI2*(1.0 - ECC)
      IF (KZ(27).GT.0.AND.RP.LT.4.0*RADIUS(J).AND.ECC.GT.0.002.AND.
     &    KSTARM(IMERGE).GE.0) THEN
          TC = 1.0D+10
          IF (KZ(27).EQ.2) THEN
              I2 = 0
*       Identify the ghost by searching single bodies.
              DO 12 K = IFIRST,N
                  IF (BODY(K).EQ.0.0D0.AND.
     &                NAME(K).EQ.NAMEG(IMERGE)) THEN
                      I2 = K
                  END IF
   12         CONTINUE
              IF (I2.GT.0) THEN
                  DO 14 K = 1,2
                      BODYI(K) = CM(K,IMERGE)
   14             CONTINUE
                  TG = 1.0
*       Evaluate the circularization time.
                  CALL HICIRC(RP,ECC,I1,I2,BODYI,TG,TC,EC,ED,W)
                  IF (TC.GT.2000) GO TO 18
              END IF
          END IF
          if(rank.eq.0)
     &    WRITE (6,15)  NAME(J), KSTAR(J), KSTARM(IMERGE), ECC, DMS,
     &                  RP, RADIUS(J), TC
   15     FORMAT (' HMDOT TERM    NAM K* E DM RP R* TC ',
     &                            I6,2I4,F8.4,F7.3,1P,3E10.2)
          ITERM = 1
          GO TO 50
      ELSE IF (ECC.GT.0.002) THEN
          IF (RP.LT.2.0*RADIUS(J)) THEN
              if(rank.eq.0)
     &        WRITE (6,15)  NAME(J), KSTAR(J), KSTARM(IMERGE), ECC, DMS,
     &                      RP, RADIUS(J)
              ITERM = 1
              GO TO 50
          END IF
      END IF
*
*       Obtain energy change from M*A = const and H = -M/(2*A) (DCH 8/96).
   18 DH = DM/SEMI0*(1.0 - 0.5*DM/ZMB0)
      HM0 = HM(IMERGE)
      HM(IMERGE) = HM(IMERGE) + DH
*
*       Form KS coordinate & velocity scaling factors (general point is OK).
      SEMI2 = -0.5*ZMB/HM(IMERGE)
      C2 = SQRT(SEMI2/SEMI0)
      V2 = 0.5*(ZMB + HM(IMERGE)*RI*(SEMI2/SEMI0))
      C1 = SQRT(V2/V20)
*
*       Re-scale KS variables to new energy (H < 0: constant eccentricity).
      DO 20 K = 1,4
          UM(K,IMERGE) = C2*UM(K,IMERGE)
          UMDOT(K,IMERGE) = C1*UMDOT(K,IMERGE)
*       Note no need to update XREL & VREL for RESTART (c.m. error cancels).
   20 CONTINUE
*
*       Reduce mass of relevant component (ghost is original second member).
      ZMU0 = CM(1,IMERGE)*CM(2,IMERGE)/ZMB0
      KM = 1
      IF (J.GE.IFIRST) KM = 2
      CM(KM,IMERGE) = CM(KM,IMERGE) - DM
*
*       Include corrections to EMERGE & EMDOT (consistency; no net effect!).
      ZMU = CM(1,IMERGE)*CM(2,IMERGE)/ZMB
      DECORR = ZMU*HM(IMERGE) - ZMU0*HM0
      EMERGE = EMERGE + DECORR
      EMDOT = EMDOT - DECORR
*       Compensate for EMERGE contribution in DEGRAV (cf. EVENTS).
      EGRAV = EGRAV - DECORR
*
*       Correct outer orbit for mass loss (use #I1 in case #J is a ghost).
      CALL HCORR(I1,DM,RSTAR)
*
*       Print some diagnostics on non-zero mass loss.
      IF (rank.eq.0.and.DMS.GT.1.0D-03) THEN
          WRITE (6,30)  NAME(J), KSTAR(J), IMERGE, KM, M1, DMS, PMIN,
     &                  PCRIT, SEMI2, SEMI1, DECORR
   30     FORMAT (' HMDOT    NAM K* IM KM M1 DM PM PC A A1 DE ',
     &                       I6,3I4,F6.2,1P,5E10.2,0P,F10.6)
      END IF
*
   50 RETURN
*
      END

max_stars_repo_name: rknop/amuse

max_issues_repo_name: rknop/amuse

max_forks_repo_name: rieder/amuse


text:
The Yolo County Animal Shelter is run by the Animal Services unit of Yolo County Sheriffs Department. It is not the same as the Yolo County SPCA  the SPCA is a privatelyrun organization that also serves animals in Yolo County.  The shelter primarily takes in Dogs dogs and Cats cats; the shelter also provides pet adoptions, redemptions, rescues, and licensing. They investigate noise complaints, inspect kennels, pick up loose animals including livestock, respond to animal attacks, rent Town Wildlife traps, and provide animal welfare checks  some of these services may require a fee. If your pet has died, you can bring it to the shelter (or to your local veterinarian for aftercare), however the Deceased Animal Deposit no longer exists. You can come to the office and let them know, and they will assist you.

It is easy to adopt an animal from the shelter. The adoption process includes a brief interview about your home environment, including questions about previous pet ownership, yard containment, etc.  If you live in a rental property, you must provide proof that your landlord will allow you to keep pets, which an Animal Services officer will verify.  All shelter animals are spayed or neutered before going to their new homes.

All animals over the age of four months must be licensed and vaccinated for rabies. Once your pet is vaccinated, your veterinarian will report this information to Animal Services, initiating you in the licensing system. The shelter will mail or fax you a license application renewable annually or you can complete the application in person at either the shelter or Davis Police Department Davis PD. Proof of a current rabies vaccination is required for all dog licenses. Licensing fees for altered cats and dogs is $15, while unaltered licensing fees are $25.

The Animal Shelter offers low cost rabies clinics on a regular basis (and low cost vaccinations are also available through privately run mobile vaccine clinics and veterinary clinics in Yolo County).  The cost of the rabies vaccination is $8.00.  Please call the office for an appointment. (5306685287)

If your pet has been impounded at the shelter, you need to show a photo ID, proof of ownership (such as a picture of your pet) and proof of current vaccinations plus spay/neuter status. You will also have to pay fees (redemption, boarding, vaccination, license, etc.) before they will release your pet. Some Wiki animals have spent time at the Yolo County Animal Shelter.

The shelter also has a volunteer program which includes providing foster homes, cat socialization, grooming, and dog walking. Also, if you have blankets, towels, pet medications or other animal items, please consider donating them to the shelter. They also can always use liquid laundry detergent! If you are interested in volunteering contact Mailto(ycasrescue@hotmail.com), for more information.

Remember, the SPCA and the Shelter are not the same. The Yolo County Animal Shelter is a division of the Sheriffs Department. It runs the shelter and you can adopt directly from them with no go between. The SPCA is a volunteer organization that works to help place homeless animals and has their own rules and procedures for placing animals.

Is anyone else appalled that Yolo County has a kill shelter? I mean, kill shelters are horrible and inhumane, and Yolo County of all places should eliminate this disgusting practice.  Users/EliseKane
 Why Yolo more than anywhere else?  Lets face it, there arent enough people adopting to keep nokill shelters in the black without massive private donations.   Users/AllisonEriksen

 Are you kidding me? Where are they going to put the dozens of animals that come in each day? Who is going to pay for a huge sanctuary for animals that cannot be adopted? I have a brilliant idea  you want to stop the killing? Donate all your spare cash to the shelter. And all your friends too. And all of Yolo County. Life aint no fairy tale and nationwide 6070% of all animals dropped off at shelters are put to sleep. But dont you dare criticize the shelter workers for dedicating their lives to animals that no one else will take responsibility for.   Users/dizzyditz

20061009 16:19:51 nbsp The Yolo County SPCA, despite how much I hate them, gets a dozen kittens A DAY. Theres just too many. It sucks, but thats life. Users/SandyKnowles

20061009 16:23:18 nbsp So I really just gotta put in my two cents about the Yolo County SPCA. I adopted (more like bought, since they charge a gajillion dollars to adopt) my cat, Fizzle, from the Yolo County SPCA in 2003. Well, she was a stray and when I adopted her she was less than 2 pounds and had a minor upper respiratory infection. The girl who answered my questions and took me through the adoption bureacracy told me about the infection, and Im not an idiot, so I wasnt concerned. They did, however, INSIST on spaying her before I could take her home. Lets recap here: shes barely 6  8 weeks old, shes less than 2 lbs (underweight), and she has an upper respiratory infection. WHAT IDIOT OPERATES ON A SICK ANIMAL?! But they did, and it nearly killed her. It cost me about $500 dollars in vet bills from postop opportunistic infections before she was well again. Therefore, I have no respect for the shelter.  Users/SandyKnowles

20061230 20:28:58 nbsp I think that its REALLY important to point out two things: 1) The animal shelter in Woodland IS NOT run by the SPCA...it is run by Yolo County. The Yolo County SPCA is a totally SEPARATE nonprofit group that simply assists in the adoption process at the shelter. They RESCUE the animals that are in need of removal from the shelter to SAVE them from being euthanized.  The animals are then placed with families throughout the county as fosters until a good home is found for them and they can be adopted.  2)It is not the SPCA nor the Countys fault that you had to have your animal spayed or neutered prior to adoption.  IT IS THE LAW THAT ALL ANIMALS BE SPAYED OR NEUTERED or they cannot be adopted out.  It is to protect the animals from overpopulation and crowding the shelter.  Otherwise its pure irresponsibility on our part.   Users/DarcyMcNie

20070317 14:38:43 nbsp I adopted a cat from the SPCA (not the shelter) and I was told that the cat would become antisocial.  My friend petsat him at her home and she said he was the most affectionate cat and was out roaming around all the time (for 1 yr and 8 months). We realized that the cat wasnt antisocial, instead he just didnt like my dogs.  In fact, he constantly sprayed in areas that my dogs like to sleep in (my clothes, bed, etc).  We realized we couldnt force him to be happy so he needed a home with no dogs so we took him back to the spca.  The people at the spca were so rude.  First, they DEMANDED a monetary donation which defeats the whole point of calling it a donation.  They said that after 2 years, you have to make a donation which is bull b/c the contract we signed never stated so.  So not only did we pay $90 for a cat that we had to take back to the spca, but they wanted additional money because we took care of the cat for 2 years?  It does not make sense to me at all. Users/TiffanyLee

20070317 14:44:20 nbsp They threatened us and told us that we couldnt adopt from spca again.  They made it sound like we abused the cat, tossed him in the trash and kicked him around for a bit.  We were responsible owners who saw our cat in distress and made a hard but thoughtful decision to give the cat back so he can find a home where he will be Happy! They literally just pust us on a list with a bunch of animal killers and abusers! They seemed like they expected us to keep the cat in an environment where he was miserable! isnt the point of the spca to encourage people to adopt and not buy from breeders or pet stores???  Their rudeness literally DRIVES people towards pet stores!   Users/TiffanyLee

20070317 17:42:20 nbsp I adopted a 4 month old puppy from the Animal Shelter today. It cost me $15.00 for the adoption fee, $9.50 for the license, and $18.50 for the first round of vaccinations for a grand total of $48.00.  I am also paying for the dog to be neutered which is $60 to a vet in Davis.  So the grand total of the adoption will be $108.00.  I think this is a reasonable cost since I would have to pay for vaccinations and neutering with any puppy.  I see that the only real cost of the adoption from the Yolo County Shelter was the $15.00 adoption fee.   Users/SharlaDaly

20070317 22:28:53 nbsp We adopted squeaks from here last April. It cost $90 and she was spayed and had been given her first shots, and had been given antibiotics when she first came in. We got a free vet visit as well to make sure she was all right and to establish us with a vet. We got Boo from a family in Vacaville for free and had to pay $300 in vets visits for her shots and to get her spayed. I think its a much better deal to go through the shelter and pay a little bit of money up front than to get a free animal and pay tons in vet visits. Plus, they were fine with me coming in several times before adopting and spending time with the cats before we were ready to adopt.  Users/ElleWeber

20080609 12:20:40 nbsp Our last three dogs came from the SPCA, and all three were/are fantastic dogs that I would never part from.  Yes, adoption prices have gone up, but everyone has a bottom line. I couldnt imagine any of our dogs being put down because of overcrowding. Weve never purchased an animal from a pet store, and never will. Users/jazzybatgirl

20080609 18:34:29 nbsp Remember, the SPCA and the Shelter are not one in the same.  The Yolo County Animal Shelter is a division of the Sheriffs Department.  It runs the shelter and you can adopt directly from them with no go between.  The SPCA is a volunteer organization that works to help place homeless animals and has their own rules and proceedures for placing animals.   Users/SharlaDaly

20090217 15:51:09 nbsp Unfortunately, there are just too many stray animals taken to the shelter that they can all be adopted; its a terrible shame, but I do understand the necessity to euthanize as a last resort.  Because I  knew this, after my cat died in 2007, I made a point of adopting two cats from the shelter.  They cost about $100 each in total.  They were spayed/neutered, had their shots and a clean bill of health, and the shelter gave us a bag of food for each cat as well as the free vet visit.  An SPCA representative also interviewed me and went over a questionnaire Id filled out about pet ownership before the adoption was approved.  I found out later from another visitor that hed wanted to adopt one of the cats I adopted, but they wouldnt allow it because he intended for the cat to be kept outdoors at all times.  Im glad that they do screen prospective owners this way, in order to provide the best chance possible for the animals.  Needless to say...my two adorable kitties are one of the bright spots in my life, and theyll have a home with me and my family all their lives.  Users/amay

20090511 22:42:17 nbsp I currently work for Yolo County Animal Services, and I have been there 2 years this June 2009. After reading all the posts on this place, I feel the need to update to currrent. Yes, this is a shelter that uses euthansia, however, we try ou ABSOLUTE best to place the animal elsewhere before coming to this final decision. We will find a foster home, a rescue, etc before coming to the last straw. The procedure is completely humane, it only takes a matter of seconds and is not at all painful. A simple injection and the animal simply falls into a peaceful sleep. It is definatley something we all wish we didnt have to do, but if there were more responsible pet owners out there, we wouldnt have to do any of that.

As far as the sick cat insisted on being spayed before it was well, I have never heard of that issue, although I dont doubt it. Since I have started working here and fostering kittens myself, I know they have to be at least 1.5 to 2.0 lbs before any such operation can occur. The cat has to be completely healthy, and adoptable.

I absolutely love working here, and taking care of the animals. I like to make sure eachone is comfortable, and let them know I love them. Because I do. All those poor animals looking sadly at me as I walk by their kennel, needs to know someone loves them, and as I rotate them from their kennels, I make sure I pet each one and make it happy and have it wag its tail. Its a happy feeling for me!

 Users/animallover

20090512 13:21:23 nbsp I adopted a beautiful dog there 17 years ago. A great experience the staff was very helpful. There were a few people wanting to adopt her, but we were the lucky family. Sadly we had to say goodbye to her a week ago. Mahina was an amazing dog and will be missed.  Users/agerould

20111108 22:41:02 nbsp for folks wanting to really understand that there ARE enough homes for ALL the animals in shelters as noted by numerous studies done by the HSUS and maddies fund. Yolo co animal control receives about 6,000 animals a year total. Stray animals, lost animals, dead animals, and those truly without homes or dropped off by owners no longer wanting/able to care for their pets. So when Elise Kane says of all places Yolo county should be no kill she is correct. This is a low volume animal shelter, there are lots of missing programs being offered to the public. WE should be ashamed that we as a community are not doing something that many many places have achieved with many more animals a year and much less money.  here visit this link if you want more info. and do some research yourself and you can see that OPEN DOOR shelters across the country have achieved no kill. Yes, animals are still euthanized... the animals so ill that they are not treatable and are suffering, and animals that are a danger or risk to the public.However NO animal is euthanized that is healthy and adoptable.   http://www.rescue50.org/politicaladvocacy.html  Users/CayceWallace

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
C Copyright 1981-2012 ECMWF.
C
C This software is licensed under the terms of the Apache Licence
C Version 2.0 which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
C
C In applying this licence, ECMWF does not waive the privileges and immunities
C granted to it by virtue of its status as an intergovernmental organisation
C nor does it submit to any jurisdiction.
C
      SUBROUTINE ROWINA3( P, KO, KI, PW, KCODE, PMSVAL, KRET,
     X                    OMISNG, OPERIO, OVEGGY)
C
C---->
C**** ROWINA3 - Interpolation of row of values.
C
C     Purpose.
C     --------
C
C     Interpolate a row of values.
C
C
C**   Interface.
C     ----------
C
C     CALL ROWINA3( P, KO, KI, PW, KCODE, PMSVAL, KRET, OMISNG, OPERIO)
C
C
C     Input Parameters.
C     -----------------
C
C     P      - Row of values to be interpolated.
C              Dimension must be at least KO.
C
C     KO     - Number of values required.
C
C     KI     - Number of values in P on input.
C
C     PW     - Working array.
C              Dimension must be at least (0:KO+2,3).
C
C     KCODE  - Interpolation required.
C              1 , linear.
C              3 , cubic.
C
C     PMSVAL - Value used for missing data indicator.
C
C     OMISNG - True if missing values are present in field.
C
C     OPERIO - True if input field is periodic.
C
C     OVEGGY     - True if 'nearest neighbour' processing must be used
C                  for interpolation
C
C     Output Parameters.
C     ------------------
C
C     P     - Now contains KO values.
C     KRET  - Return code
C             0, OK
C             Non-zero, error
C
C
C     Method.
C     -------
C
C     Linear or cubic interpolation performed as required.
C
C
C     Externals.
C     ----------
C
C     SCM0
C
C
C     Reference.
C     ----------
C
C     None.
C
C
C     Comments.
C     ---------
C
C     This is a version of ROWINA which allows for missing data
C     values and hence for bitmapped fields.
C
C
C     Author.
C     -------
C
C     J.D.Chambers    ECMWF     22.07.94
C
C
C     Modifications.
C     --------------
C
C     J.D.Chambers    ECMWF     13.09.94
C     Add return code KRET and remove calls to ABORT.
C
C     J. Clochard, Meteo France, for ECMWF - January 1998.
C     Addition of OMISNG and OPERIO arguments.
C
C
C     -----------------------------------------------------------------
C----<
C
C*******************************************************************
C     Section 0.  Declarations.
C*******************************************************************
C
      IMPLICIT NONE
C
C#include "grprs.h"
C
C     Subroutine arguments
      REAL P, PW, PMSVAL
      INTEGER KO,KI,KCODE, KRET
      DIMENSION P(KO),PW(0:KO+2,3)
C
C     Local variables
      INTEGER JL, IP
      REAL ZRDI, ZDO, ZPOS, ZWT, ZWT1
C
      LOGICAL OMISNG, OPERIO, OVEGGY
C
C*******************************************************************
C     Section 1.  Linear interpolation ..
C*******************************************************************
C
  100 CONTINUE
C
      KRET = 0
C
      IF( KCODE.EQ.1 ) THEN
C
C        Move input values to work array
         DO 110 JL = 1, KI
            PW(JL,1) = P(JL)
  110    CONTINUE
C
         IF (OPERIO) THEN
C
C        Arrange wrap-around value in work array
           PW(KI+1,1) = P(1)
C
C          Set up constants to be used to figure out weighting for
C          values in interpolation.
           ZRDI = REAL(KI)
           ZDO  = 1.0 / REAL(KO)
         ELSE
C
C        Repeat last value, to cope with "implicit truncation" below
           PW(KI+1,1) = P(KI)
C
C          Set up constants to be used to figure out weighting for
C          values in interpolation.
           ZRDI = REAL(KI-1)
           ZDO  = 1.0 / REAL(KO-1)
         ENDIF
C
C        Loop through the output points
         IF ( .NOT. OMISNG ) THEN
C
            DO 120 JL = 1, KO
C
C           Calculate weight from the start of row
            ZPOS = (JL-1) * ZDO
            ZWT  = ZPOS * ZRDI
C
C           Get the current array position(minus 1) from the weight -
C           note the implicit truncation.
            IP   = ZWT
C
C           Adjust the weight to range (0.0 to 1.0)
            ZWT  = ZWT - IP
C
C           If 'nearest neighbour' processing must be used
C
            IF( OVEGGY ) THEN
C
              IF( ZWT.LT.0.5 ) THEN
                P(JL) = PW(IP+1,1)
              ELSE
                P(JL) = PW(IP+2,1)
              ENDIF
            ELSE
C
C             Interpolate using the weighted values on either side
C             of the output point position
              P(JL) = (1.0-ZWT) * PW(IP+1,1) + ZWT * PW(IP+2,1)
C
            ENDIF
C
  120       CONTINUE
C
         ELSE
C
         DO 130 JL = 1, KO
C
C           Calculate weight from the start of row
            ZPOS = (JL-1) * ZDO
            ZWT  = ZPOS * ZRDI
C
C           If 'nearest neighbour' processing must be used
C
            IF( OVEGGY ) THEN
C
              IF( ZWT.LT.0.5 ) THEN
                P(JL) = PW(IP+1,1)
              ELSE
                P(JL) = PW(IP+2,1)
              ENDIF
            ELSE
C
C             Get the current array position(minus 1) from the weight -
C             note the implicit truncation.
              IP   = ZWT
C
C             If the left value is missing, use the right value
              IF ( PW(IP+1,1) .EQ. PMSVAL ) THEN
                 P(JL) = PW(IP+2,1)
C
C             If the right value is missing, use the left value
              ELSE IF ( PW(IP+2,1) .EQ. PMSVAL ) THEN
                 P(JL) = PW(IP+1,1)
C
C             If neither missing, interpolate ...
              ELSE
C
C               Adjust the weight to range (0.0 to 1.0)
                ZWT  = ZWT - IP
C
C               Interpolate using the weighted values on either side
C               of the output point position
                P(JL) = (1.0-ZWT) * PW(IP+1,1) + ZWT * PW(IP+2,1)
              ENDIF
C
            ENDIF
C
  130    CONTINUE
C
       ENDIF
C
C*******************************************************************
C     Section 2.  Cubic interpolation ..
C*******************************************************************
C
  200 CONTINUE
C
      ELSEIF(KCODE.EQ.3) THEN
C
         IF ( OMISNG ) THEN
C            WRITE(GRPRSM,*)
            WRITE(*,*)
     X        'ROWINA3: Cubic interpolation not supported'
C            WRITE(GRPRSM,*)
            WRITE(*,*)
     X        'ROWINA3: for fields containing missing data.'
C            WRITE(GRPRSM,*) ' Sorry!'
            WRITE(*,*)
            KRET = 1
            GOTO 900
         ENDIF
C
         DO 210 JL = 1,KI
            IF ( P(JL) .EQ. PMSVAL ) THEN
C              WRITE(GRPRSM,*)
              WRITE(*,*)
     X          'ROWINA3: Cubic interpolation not supported'
C              WRITE(GRPRSM,*)
              WRITE(*,*)
     X          'ROWINA3: for fields containing missing data.'
C              WRITE(GRPRSM,*) ' Sorry!'
              WRITE(*,*) ' Sorry!'
              KRET = 1
              GOTO 900
            ENDIF
            PW(JL,1) = P(JL)
  210    CONTINUE
         PW(0,1) = P(KI)
         PW(KI+1,1) = P(1)
         PW(KI+2,1) = P(2)
         DO 220 JL = 1,KI
            PW(JL,2) =  - PW(JL-1,1)/3.0 - 0.5*PW(JL,1)
     1             + PW(JL+1,1)    - PW(JL+2,1)/6.0
            PW(JL+1,3) =    PW(JL-1,1)/6.0 - PW(JL,1)
     1             + 0.5*PW(JL+1,1) + PW(JL+2,1)/3.0
  220    CONTINUE
         CALL SCM0(PW(1,2),PW(2,3),PW(1,1),PW(2,1),KI)
         ZRDI = REAL(KI)
         ZDO = 1.0/REAL(KO)
         DO 230 JL = 1,KO
            ZPOS = (JL-1)*ZDO
            ZWT = ZPOS*ZRDI
            IP = ZWT+1
            ZWT = ZWT+1.0-IP
            ZWT1  =  1.0 - ZWT
            P(JL) = ((3.0-2.0*ZWT1)*PW(IP,1) + ZWT*PW(IP,2))*ZWT1*ZWT1
     1       + ((3.0-2.0*ZWT) *PW(IP+1,1) - ZWT1*PW(IP+1,3))*ZWT*ZWT
  230    CONTINUE
C
      ELSE
C
C*******************************************************************
C     Section 3.  Invalid interpolation code ..
C*******************************************************************
C
  300   CONTINUE
C
C        WRITE(GRPRSM,9001) KCODE
        WRITE(*,9001) KCODE
        KRET = 2
      ENDIF
C
C*******************************************************************
C     Section 9.  Closedown.
C*******************************************************************
C
900   CONTINUE
      RETURN
C
 9001 FORMAT (1H ,'ROWINA3 : Invalid interpolation code = ',I4)
C
      END

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


text:
      SUBROUTINE HQ9(CODE)	!Implement the rather odd HQ9+ instruction set.
       CHARACTER*(*) CODE	!One operation code per character.
       INTEGER I,B	!Steppers.
       INTEGER A	!An accumulator.
        A = 0		!Initialised.
        DO I = 1,LEN(CODE)	!Step through the code.
          SELECT CASE(CODE(I:I))!Inspect the operation code.
           CASE(" ")		!Might as well do nothing.
           CASE("+")		!Increment the accumulator.
            A = A + 1		!Thus. Though, nothing refers to it...
           CASE("h","H")	!Might as well allow upper or lower case.
            WRITE (6,*) "Hello, world!"	!Hi there!
           CASE("q","Q")	!Show the (rather questionable) code.
            WRITE (6,*) CODE	!Thus.
           CASE("9")		!Recite "99 bottles of beer"...
            DO B = 99,2,-1	!Grammar is to be upheld, so the singular case is special.
              WRITE (6,1) B,"bottles"," on the wall,",B,"bottles","."	!Two lots: number, text, text.
    1         FORMAT (I2,1X,A," of beer",A)	!Exhausted by the first triplet, so a new line for the second.
              WRITE (6,2)			!Now for the reduction.
    2         FORMAT ("Take one down, pass it around,")	!Announce.
              IF (B.GT.2) WRITE (6,1) B - 1,"bottles"," on the wall."	!But, not for the singular state.
            END DO		!Recite the next stanza.
            WRITE (6,1) 1,"bottle"," on the wall,",1,"bottle","."	!The singular case. No longer "bottles".
            WRITE (6,2)						!There's nothing so lonesome, morbid or drear,
            WRITE (6,*) "No bottles of beer on the wall."	!Than to stand at the bar of a pub with no beer.
            WRITE (6,*) "Go to the store, buy some more."	!Take action.
           CASE DEFAULT		!Denounce any unknown operation codes.
            WRITE (6,*) "Unrecognised code:",CODE(I:I)	!This is why a space is treated separately.
          END SELECT		!So much for that operation code.
        END DO			!On to the next.
      END SUBROUTINE HQ9	!That was odd.

      PROGRAM POKE
      CALL HQ9("hq9")
      END

max_issues_repo_name: seanwallawalla-forks/RosettaCodeData

max_forks_repo_name: seanwallawalla-forks/RosettaCodeData


text:
American Graffiti Departed Businesses was a small tattoo and piercing parlor located in the basement complex on E st.

Ive always had very good experiences with the folks here Users/AmeliaCarlson aec

The people that work there are always really nice, but they dont have a huge selection of jewelry  Users/SharonZimmerman Sharon Zimmerman

Clean, Professional, no complaints here!

Im not sure what the deal is, but I went there today (7/25/06) and the shop was no longer at this location. Users/NilaGurusinghe Nila


text:
C--------------------------------------------------------------------------
        subroutine ftikys(ounit,keywrd,strval,comm,status)

C       insert a string keyword into the header at the current position
C
C       ounit   i  fortran output unit number
C       keywrd  c  keyword name    ( 8 characters, cols.  1- 8)
C       strval  c  keyword value
C       comm    c  keyword comment
C       OUTPUT PARAMETERS
C       status  i  output error status (0 = ok)
C
C       written by Wm Pence, HEASARC/GSFC, March 1993
C       Modifed 9/94 to call FTPKLS, supporting the OGIP long string convention

        character*(*) keywrd,comm,strval
        integer ounit,status

C-------COMMON BLOCK DEFINITIONS:--------------------------------------------
        integer nb,ne
        parameter (nb = 20)
        parameter (ne = 200)
        integer bufnum,chdu,hdutyp,maxhdu,hdstrt,hdend,nxthdr,dtstrt
        integer nxtfld
        logical wrmode
        common/ft0001/bufnum(199),chdu(nb),hdutyp(nb),maxhdu(nb),
     &  wrmode(nb),hdstrt(nb,ne),hdend(nb),nxthdr(nb),dtstrt(nb),nxtfld
C-------END OF COMMON BLOCK DEFINITIONS:------- -----------------------------

        integer lenval,length,i,nspace,ibuff,nexthd,endhd,nkeys,keypos

        if (status .gt. 0)return

C       find how many keywords are required to write the string, in case it
C       cannot fit onto one keyword and has to be continued on multiple lines.

        lenval=len(strval)
        length=0
        do 10 i=lenval,1,-1
                if (strval(i:i) .ne. ' ')then
                        length=i
                        go to 20
                end if
10      continue
20      nspace=max(1,(length-2)/67+1)

C       save current pointer values
        ibuff=bufnum(ounit)
        endhd=hdend(ibuff)
        nexthd=nxthdr(ibuff)

C       insert enough spaces in the header at the current location
        call ftghps(ounit,nkeys,keypos,status)

        do 30 i=1,nspace
            call ftirec(ounit,keypos,' ',status)
30      continue

C       temporarily reset position of the end of header to force keyword
C       to be written at the current header position.
        hdend(ibuff)=nexthd

C       write the keyword (supporting the OGIP long string convention)
        call ftpkls(ounit,keywrd,strval,comm,status)

C       reset the next keyword pointer to follow the inserted keyword
        nxthdr(ibuff)=nexthd+80*nspace

C       reset the end-of-header pointer to its real location
        hdend(ibuff)=endhd+80*nspace
        end

max_stars_repo_name: iraf-community/stsdas

max_issues_repo_name: kirxkirx/iraf

max_forks_repo_name: kirxkirx/iraf


SUBROUTINE sweep_ahotn_nefd(g)

!-------------------------------------------------------------
!
!  Sweeps across the 3-D matrix
!   Starts at top, far, right corner (mu, eta, xi < 0), then sweeps
!   down all planes and rows, accounting for reflection if necessary. Then
!   sweeps up all planes and rows for xi>0
!
!-------------------------------------------------------------

USE invar
USE solvar
use kernel_module
use precision_module, only: dp

IMPLICIT NONE
INTEGER, INTENT(IN) :: g
INTEGER :: xs, xe, ys, ye, zs, ze, incx, incy, incz, ord, nfy, nfz
INTEGER :: i, j, k, t, u, v, m, n, ydir, xdir, zdir
INTEGER :: ix,iy,iz,jx,jy,jz,indx

REAL(kind=dp), DIMENSION(0:lambda, 0:lambda) :: fx
REAL(kind=dp), DIMENSION(0:lambda, 0:lambda,nx, 2) :: fy
REAL(kind=dp), DIMENSION(0:lambda, 0:lambda,nx, ny, 2, 2) :: fz
REAL(kind=dp), DIMENSION(ordcb) :: b
REAL(kind=dp) :: sig, mu, eta, xi, x, y, z, c, sgn

! Initialize the flux solution to zero
f=0.0d0

! Start with loop over all angles
DO n = 1, apo
   ! Set up the angles
   mu  = ang(n,1)
   eta = ang(n,2)
   xi  = ang(n,3)

   ! Loop over xi<0 then xi>0
   DO zdir = 1, 2
      IF (zdir == 1) THEN
         zs = nz
         ze = 1
         incz = -1
         ! Set boundary conditions
         IF (zebc==0) THEN
            fz=0.0
         ELSE IF (zebc==2) THEN
            !
            DO iy=1,ny
               DO ix=1,nx
                  DO jy=0,lambda
                     DO jx=0,lambda
                        fz(jx,jy,ix,iy,2,2)=tobc(jx,jy,ix,iy,n,1)  ! mu>0,eta>0
                        fz(jx,jy,ix,iy,2,1)=tobc(jx,jy,ix,iy,n,2)  ! mu>0,eta<0
                        fz(jx,jy,ix,iy,1,2)=tobc(jx,jy,ix,iy,n,3)  ! mu<0,eta>0
                        fz(jx,jy,ix,iy,1,1)=tobc(jx,jy,ix,iy,n,4)  ! mu<0,eta<0
                     END DO
                  END DO
               END DO
            END DO
            !
         END IF
      ELSE IF (zdir == 2) THEN
         zs = 1
         ze = nz
         incz = 1
         ! Set boundary conditions
         IF (zsbc==0) THEN
            fz=0.0
         ELSE IF (zsbc==1) THEN
            fz = zsbc*fz
         ELSE IF (zsbc==2) THEN
            !
            DO iy=1,ny
               DO ix=1,nx
                  DO jy=0,lambda
                     DO jx=0,lambda
                        fz(jx,jy,ix,iy,2,2)=bobc(jx,jy,ix,iy,n,1)  ! mu>0,eta>0
                        fz(jx,jy,ix,iy,2,1)=bobc(jx,jy,ix,iy,n,2)  ! mu>0,eta<0
                        fz(jx,jy,ix,iy,1,2)=bobc(jx,jy,ix,iy,n,3)  ! mu<0,eta>0
                        fz(jx,jy,ix,iy,1,1)=bobc(jx,jy,ix,iy,n,4)  ! mu<0,eta<0
                     END DO
                  END DO
               END DO
            END DO
            !
         END IF
         !
      END IF

   ! Start the loop in the negative z-direction, then do positive z-direction
   DO k = zs, ze, incz
      z = dz(k)

      ! Loop over eta<0 then eta>0
      DO ydir = 1, 2
         IF (ydir == 1) THEN
            ys = ny
            ye = 1
            incy = -1
            nfz = 1
            ! Set back face boundary conditions
            IF (yebc==0) THEN
              fy=0.0
            ELSE IF (yebc==2 .and. incz>0 ) THEN
              !
              DO ix=1,nx
                 DO jz=0,lambda
                    DO jx=0,lambda
                      fy(jx,jz,ix,1)=babc(jx,jz,ix,k,n,2) ! xi>0, mu<0
                      fy(jx,jz,ix,2)=babc(jx,jz,ix,k,n,1) ! xi>0, mu>0
                    END DO
                 END DO
              END DO
              !
            ELSE IF (yebc==2 .and. incz<0 ) THEN
              !
              DO ix=1,nx
                 DO jz=0,lambda
                    DO jx=0,lambda
                      fy(jx,jz,ix,1)=babc(jx,jz,ix,k,n,4) ! xi<0, mu<0
                      fy(jx,jz,ix,2)=babc(jx,jz,ix,k,n,3) ! xi<0, mu>0
                    END DO
                 END DO
              END DO
              !
            END IF
            !
         ELSE IF (ydir == 2) THEN
            ys = 1
            ye = ny
            incy = 1
            nfz = 2
            ! Set front face boundary conditions
            IF (ysbc==0) THEN
              fy=0.0
            ELSE IF (ysbc==1) THEN
              fy = fy
            ELSE IF (ysbc==2 .and. incz>0 ) THEN
              !
              DO ix=1,nx
                 DO jz=0,lambda
                    DO jx=0,lambda
                      fy(jx,jz,ix,1)=frbc(jx,jz,ix,k,n,2) ! xi>0, mu<0
                      fy(jx,jz,ix,2)=frbc(jx,jz,ix,k,n,1) ! xi>0, mu>0
                    END DO
                 END DO
              END DO
              !
            ELSE IF (ysbc==2 .and. incz<0 ) THEN
              !
              DO ix=1,nx
                 DO jz=0,lambda
                    DO jx=0,lambda
                      fy(jx,jz,ix,1)=frbc(jx,jz,ix,k,n,4) ! xi<0, mu<0
                      fy(jx,jz,ix,2)=frbc(jx,jz,ix,k,n,3) ! xi<0, mu>0
                    END DO
                 END DO
              END DO
              !
            END IF
         END IF

         ! Start the loop in negative y-direction, then do positve y-direction
      DO j = ys, ye, incy
         y = dy(j)

         ! Perform two loops, one in negative x-direction, then positive
         DO xdir = 1, 2
            IF (xdir == 1) THEN
               xs = nx
               xe = 1
               incx = -1
               nfy = 1
               IF(xebc==0) THEN
                  fx=0.0
               ELSE IF(xebc==2 .and. incy>0 .and. incz>0) THEN
                  fx=ribc(:,:,j,k,n,1)
               ELSE IF(xebc==2 .and. incy>0 .and. incz<0) THEN
                  fx=ribc(:,:,j,k,n,2)
               ELSE IF(xebc==2 .and. incy<0 .and. incz>0) THEN
                  fx=ribc(:,:,j,k,n,3)
               ELSE IF(xebc==2 .and. incy<0 .and. incz<0) THEN
                  fx=ribc(:,:,j,k,n,4)
               END IF
            ELSE IF (xdir == 2) THEN
               xs = 1
               xe = nx
               incx = 1
               nfy = 2
               ! Reset the incoming x-flux for the x-lo bc
               IF(xsbc==0) THEN
                  fx=0.0
               ELSE IF(xsbc==1) THEN
                  fx=fx
               ELSE IF(xsbc==2 .and. incy>0 .and. incz>0) THEN
                  fx=lebc(:,:,j,k,n,1)
               ELSE IF(xsbc==2 .and. incy>0 .and. incz<0) THEN
                  fx=lebc(:,:,j,k,n,2)
               ELSE IF(xsbc==2 .and. incy<0 .and. incz>0) THEN
                  fx=lebc(:,:,j,k,n,3)
               ELSE IF(xsbc==2 .and. incy<0 .and. incz<0) THEN
                  fx=lebc(:,:,j,k,n,4)
               END IF
            END IF

         ! Start the loop in the negative x-direction
         DO i = xs, xe, incx

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

            x = dx(i)
            m = mat(i,j,k)
            sig = sigt(m,g)
            ord = lambda
            c = sigs(m,g,g)/sig     ! Scattering ratio

            ! Prepare source
            DO v = 0, lambda
               DO u = 0, lambda
                  DO t = 0, lambda
!                     cell_source(t,u,v)=c*e(t,u,v,i,j,k) + s(t,u,v,i,j,k,g)/sig
                      indx = ordsq*t + order*u + v + 1
                      b(indx)=c*e(t,u,v,i,j,k) + s(t,u,v,i,j,k,g)/sig
                  END DO
               END DO
            END DO

            ! call AHOTN NEFD kernel
            call  ahotn_nefd_kernel(x,y,z,mu,eta,xi,incx,incy,incz,sig,c,fx,fy(:,:,i,nfy),fz(:,:,i,j,nfy,nfz),b,lambda,ordcb)

            ! Update the scalar flux solution
            DO v = 0, lambda
               DO u = 0, lambda
                  DO t = 0, lambda
                     indx = ordsq*t + order*u + v + 1
                     f(t,u,v,i,j,k,g) = f(t,u,v,i,j,k,g) + w(n)*b(indx)
                  END DO
               END DO
            END DO

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

         ! End loop over x cells
         END DO
         ! End loop over negative and positive x-directions
         END DO

      ! End loop over y cells
      END DO
      ! End loop over negative and positve y-directions
      END DO

   ! End loop over z cells
   END DO
   ! End loop over negative and positive z-directions
   END DO

! End loop over angles
END DO

RETURN
END SUBROUTINE sweep_ahotn_nefd