text:
*...................................................................

      subroutine ODS_NCVPTR ( ncid,   varid,
     .                        start,  count,
     .                        values, ierr )

      implicit NONE

!-------------------------------------------------------------------------
!         NASA/GSFC, Data Assimilation Office, Code 910.3, GEOS/DAS      !
!-------------------------------------------------------------------------
!
! !ROUTINE: ODS_NCVPTR
!
! !DESCRIPTION:
!     This routine writes an array of native real numbers to a
!     file using NetCDF interface routines. The routine invokes
!     other modules to convert the values from native format to
!     a format required by the NetCDF routines, and then writes
!     the data to the file.  If the attributes, scale_factor and
!     add_offset, are defined in the NetCDF file, then the values
!     are also scaled using these attributes.

!     NOTE: This routine does perform checks to verify whether
!     the range of values on input is consistent with the
!     internal variable type in the file. For example, an
!     floating point variable here could be stored as a two
!     byte real variable on file.  Therefore, all elements in
!     the array to be stored as two byte integers should have
!     values between -32,767 and 32,767.  Be careful.  If the
!     software does detect a number inconsistent with the
!     variable type, then the routine returns with an error
!     message without saving at least some of the  values.
!     This check prevents overflows from occurring which would
!     produce unexpected results and probable abnormal program
!     termination.  The range of values for the internal
!     variable types are as follows:
!
!       type                      minimum value    maximum value
!       ----                      -------------    -------------
!       byte                               -128              255
!       short (2 byte) integer          -32,767          -32,767
!       long  (4 byte) integer   -2,147,483,643    2,147,483,643
!       real  (4 bytes )               -3.40e38          3.40e38
!
!     A value of 256 is added to any negative integer greater than
!     -128 that is to be stored as a byte.
!
! !INTERFACE:
!     call ODS_NCVPTR ( ncid, varid, start, count, values, ierr )
!
      integer   ncid         ! NetCDF file id
      integer   varid        ! NetCDF variable id
      integer   start  ( * ) ! NetCDF file indicies specifying the
                             !   location of the corner of the
                             !   hyperslab where the first of the
                             !   data values will be written.  A
                             !   hyperslab is a multidimensional
                             !   block of data within a NetCDF
                             !   file.
      integer   count  ( * ) ! The edge lengths of the NetCDF
                             !   file hyperslab.
      real      values ( * ) ! Values to be written
!
! !OUTPUT PARAMETER:
      integer   ierr         ! Return error code
!
! !LIBRARIES ACCESSED:
!     NetCDF
!
! !FILES USED:
!     netcdf.inc, a header file, for defining NetCDF library
!            parameters
!     ods_worksp.h, a header file, for defining hardwired
!            constants and defining global variables and
!            setting up data structures for work space
!     ods_stdio.h, a header file, for defining standard input/output
!            unit numbers
!
! !REVISION HISTORY:
!     05Jun1996   C. Redder   Origional version
!     13Sep2002   C. Redder   Made changes in calling ods_defwsp in order
!                             process and store one-byte integers in a
!                             multi-character string rather than in an
!                             array of one-character strings
!
!-------------------------------------------------------------------------

      include 'netcdf.inc'
      include 'ods_worksp.h'
      include 'ods_stdio.h'

*     NetCDF Variable information
*     ---------------------------
      character VarNam * ( MaxNCNam ) ! name of variable
      integer   NC_VarType            ! type of NetCDF variable
      integer   NVDims                ! number of NetCDF dimensions
      integer   VDims    ( MaxVDims ) ! NetCDF dimensions
      integer   NVAtts                ! number of attributes

*     Parameters determining the usage of work space
*     ----------------------------------------------
      integer   nblocks               ! number of blocks of values
                                      !   that are required to
                                      !   write all data via the
                                      !   work space.  This
                                      !   parameter is also the
                                      !   number of NetCDF I/O
                                      !   routines calls
      integer   iblock                ! block index number for do
                                      !   loop
      integer   nval                  ! number of values in the block
      integer   ival                  ! value index number for do loop
      integer   iworksp               ! index number for work space

*     Default value for error code
*     ----------------------------
      ierr    = NCNoErr

*     Get information about the NetCDF variable
*     -----------------------------------------
      call ncvinq ( ncid,   varid,
     .              VarNam, NC_VarType,
     .              NVDims, VDims,
     .              NVAtts, ierr )
      if ( ierr .ne. NCNoErr ) return

*     The data is being partitioned into subslabs ( or blocks or
*     subsets).  Each subslab is processed using the work space
*     as defined in the header file, ods_worksp.h.  Set up work
*     space and determine NBlocks for ...
*     ----------------------------------------------------------
      if ( NC_VarType .eq. NCByte ) then   ! ... one byte integers
         call ODS_DefWSp ( NVDims, start, count, Max_CWorkSp, nblocks )

      else                                 ! ... for all other types
         call ODS_DefWSp ( NVDims, start, count, Max_WorkSp,  nblocks )

      end if

*     For each block of values...
*     ---------------------------
      do 20, iblock = 1, nblocks

*        Get indecies defining block in 1-D array and NetCDF hyperslab
*        -------------------------------------------------------------
         call ODS_ValWSp ( iblock,         ival_beg,    ival_end    )
         call ODS_NCVWSp ( iblock, NVDims, start_block, count_block )

*        Determine the number of values in the block
*        -------------------------------------------
         nval = ival_end - ival_beg + 1

*        Transfer all values in the block to the work space
*        --------------------------------------------------
         do 10, ival = ival_beg, ival_end
            iworksp           = ival - ival_beg + 1
            R_Val ( iworksp ) = values ( ival )
 10      continue

*        Check to determine if the values are within the range
*        as specified by the NetCDF file attributes
*        -----------------------------------------------------
         call ODS_CheckR    ( ncid, varid, nval, R_Val, ierr )
         if ( ierr .ne. NCNoErr ) return

*        Scale the array of values
*        -------------------------
         call ODS_ScaleRRev ( ncid, varid, nval, R_Val, ierr )
         if ( ierr .ne. NCNoErr ) return

*        Convert and write array of values to
*        file if the NetCDF variable type is ...
*        ---------------------------------------

*        one byte integer
*        -----------------
         if      ( NC_VarType .eq. NCByte  ) then

            call ODS_NCP_RtoI1 ( ncid,        varid,
     .                           start_block, count_block,
     .                           nval,        R_Val,
     .                           ierr )

*        two byte integers
*        -----------------
         else if ( NC_VarType .eq. NCShort ) then

            call ODS_NCP_RtoI2 ( ncid,        varid,
     .                           start_block, count_block,
     .                           nval,        R_Val,
     .                           ierr )

*        four byte integers
*        ------------------
         else if ( NC_VarType .eq. NCLong  ) then

            call ODS_NCP_RtoI4 ( ncid,        varid,
     .                           start_block, count_block,
     .                           nval,        R_Val,
     .                           ierr )

*        four byte real numbers
*        ----------------------
         else if ( NC_VarType .eq. NCFloat ) then

            call ODS_NCP_RtoR4 ( ncid,        varid,
     .                           start_block, count_block,
     .                           nval,        R_Val,
     .                           ierr )

         else

            write ( stderr, 901 )
            ierr = NCEBadTy

         end if

*        If there is an error, print out variable name
*        ---------------------------------------------
         if ( ierr .ne. NCNoErr ) then
            write ( stderr, 902 ) VarNam
            return
         end if

 20   continue

*     Insure that disk copy of file is updated
*     ----------------------------------------
*      call NCSNC ( ncid, ierr )
*      if ( ierr .ne. NCNoErr ) return

      return
*     ------

 901  format ( /, ' ODS_NCVPTR : not a recognized NetCDF ',
     .                          'variable type ' )
 902  format ( /, ' ODS_NCVPTR : variable name is ', a )

      end

max_stars_repo_name: GEOS-ESM/AeroApps

max_forks_repo_name: GEOS-ESM/AeroApps


text:
c     Subroutine to get quality flag vector
c     AJ_Kettle, 09Jan2019

      SUBROUTINE get_quality_flag_vector2(f_ndflag,l_channel,
     +  f_prcp_datavalue_mm,f_tmin_datavalue_c,
     +  f_tmax_datavalue_c,f_tavg_datavalue_c,
     +  f_snwd_datavalue_mm,f_snow_datavalue_mm,
     +  f_awnd_datavalue_ms,f_awdr_datavalue_deg,
     +  f_wesd_datavalue_mm,
     +  s_prcp_qflag,s_tmin_qflag,
     +  s_tmax_qflag,s_tavg_qflag,
     +  s_snwd_qflag,s_snow_qflag,
     +  s_awnd_qflag,s_awdr_qflag,
     +  s_wesd_qflag,
     +  s_vec_qc_code)

      IMPLICIT NONE
c************************************************************************
c     Declare variables passed into program

      REAL                :: f_ndflag
      INTEGER             :: l_channel

      REAL                :: f_prcp_datavalue_mm
      REAL                :: f_tmin_datavalue_c
      REAL                :: f_tmax_datavalue_c
      REAL                :: f_tavg_datavalue_c
      REAL                :: f_snwd_datavalue_mm
      REAL                :: f_snow_datavalue_mm
      REAL                :: f_awnd_datavalue_ms
      REAL                :: f_awdr_datavalue_deg
      REAL                :: f_wesd_datavalue_mm

      CHARACTER(LEN=*)    :: s_prcp_qflag
      CHARACTER(LEN=*)    :: s_tmin_qflag
      CHARACTER(LEN=*)    :: s_tmax_qflag
      CHARACTER(LEN=*)    :: s_tavg_qflag
      CHARACTER(LEN=*)    :: s_snwd_qflag
      CHARACTER(LEN=*)    :: s_snow_qflag
      CHARACTER(LEN=*)    :: s_awnd_qflag
      CHARACTER(LEN=*)    :: s_awdr_qflag
      CHARACTER(LEN=*)    :: s_wesd_qflag

      CHARACTER(LEN=*)    :: s_vec_qc_code(l_channel)
c*****
c     Variables used within program

      INTEGER             :: i,j,k,ii,jj,kk
      CHARACTER(LEN=1)    :: s_badflag(14)
      INTEGER             :: l_badflag

      INTEGER             :: i_index

      CHARACTER(LEN=1)    :: s_single_qflag
c************************************************************************
c      print*,'just entered get_quality_flag_vector'

c      print*,'f_ndflag=',f_ndflag
c      print*,'l_channel=',l_channel
c************************************************************************
c     Define list of bad flags
      s_badflag(1) ='D'
      s_badflag(2) ='G'
      s_badflag(3) ='I'
      s_badflag(4) ='K'
      s_badflag(5) ='L'
      s_badflag(6) ='M'
      s_badflag(7) ='N'
      s_badflag(8) ='O'
      s_badflag(9) ='R'
      s_badflag(10)='S'
      s_badflag(11)='T'
      s_badflag(12)='W'
      s_badflag(13)='X'
      s_badflag(14)='Z'

      l_badflag=14
c************************************************************************
c     Initialize the QC variable to missing=3
      DO i=1,l_channel
       s_vec_qc_code(i)='3'
      ENDDO

c     Call subroutine to determine 1 qc_vector element
      i_index=1
      s_single_qflag=s_prcp_qflag
      CALL qc_single_variable2(l_channel,s_vec_qc_code,i_index,
     +  l_badflag,s_badflag,
     +  s_single_qflag)

      i_index=2
      s_single_qflag=s_tmin_qflag
      CALL qc_single_variable2(l_channel,s_vec_qc_code,i_index,
     +  l_badflag,s_badflag,
     +  s_single_qflag)

      i_index=3
      s_single_qflag=s_tmax_qflag
      CALL qc_single_variable2(l_channel,s_vec_qc_code,i_index,
     +  l_badflag,s_badflag,
     +  s_single_qflag)

      i_index=4
      s_single_qflag=s_tavg_qflag
      CALL qc_single_variable2(l_channel,s_vec_qc_code,i_index,
     +  l_badflag,s_badflag,
     +  s_single_qflag)

      i_index=5
      s_single_qflag=s_snwd_qflag
      CALL qc_single_variable2(l_channel,s_vec_qc_code,i_index,
     +  l_badflag,s_badflag,
     +  s_single_qflag)

      i_index=6
      s_single_qflag=s_snow_qflag
      CALL qc_single_variable2(l_channel,s_vec_qc_code,i_index,
     +  l_badflag,s_badflag,
     +  s_single_qflag)

      i_index=7
      s_single_qflag=s_awnd_qflag
      CALL qc_single_variable2(l_channel,s_vec_qc_code,i_index,
     +  l_badflag,s_badflag,
     +  s_single_qflag)

      i_index=8
      s_single_qflag=s_awdr_qflag
      CALL qc_single_variable2(l_channel,s_vec_qc_code,i_index,
     +  l_badflag,s_badflag,
     +  s_single_qflag)

      i_index=9
      s_single_qflag=s_wesd_qflag
      CALL qc_single_variable2(l_channel,s_vec_qc_code,i_index,
     +  l_badflag,s_badflag,
     +  s_single_qflag)

c      print*,'s_vec_qc_code=',(s_vec_qc_code(i),i=1,l_channel)
c      print*,'just leaving get_quality_flag_vector'

      RETURN
      END

ext: f

max_stars_repo_name: glamod/glamod-nuim

max_issues_repo_name: glamod/glamod-nuim

max_forks_repo_name: glamod/glamod-nuim


text:
Hibachi Grill and Take Out was a Japanese and Korean Restaurants restaurant which featured teriyaki, spicy pork, and bulgogi. They also served Korean BBQ. Do not confuse them with HibachiSan.

20050706 07:05:40 This place offers simple, tasty rice bowls, Chinese chicken salad and even smoothies. They have sushi, but it can be a little funky. Id avoid it, and go for the good old chicken teriyaki bowl.  Users/SiennaGrass

20060319 18:08:47 nbsp The owner is one of the nicest owners Ive ever met.  He makes friends with all his clients.  I highly reccommend the spicy pork.  Spicy Pork bowls, while not on the menu, are also available if you ask for it.  Their shakes are also superb. Users/KevinRollins

20060516 09:00:51 nbsp The only thing Ive ever tried on the menu is the Chicken Tariyaki rice bowl, but its always been great.   I never liked teriyaki very much until I tried Hibachis.  Lately, Ive been going here for lunch several times a month. Users/AdamSchneider

20060710 14:36:30 nbsp This place is a little gem.  So far, the restaurant is immaculately clean and the food ingredients seem fresh, and the dishes are tasty. (Ive had the chicken teriyaki and the hibap).  I concur with previous comments: the owner IS very nice.  According to him, the place has been open since about January 2006.  Check it out. Users/AnnieSirrah

20060712 23:29:54 nbsp The owner is very, very nice.  The food is good.  Perhaps the bowls are a bit small.  I had the BBQ dinner tonight and it was quite good. Users/MisterProfessor

20061018 16:12:04 nbsp For the amount that you pay, you dont get a lot food. I paid $4.60 plus tax for a dinky chicken salad, and it was mostly lettuce and sauce; barely any chicken in it. Users/MattLow

20061206 09:19:53 nbsp Im coming here again... Users/CarlosOverstreet

20070210 17:31:00 nbsp This place has been closed for a week or more.  Is it going to open again?  I hope the owners are okay. Users/MisterProfessor

20070211 12:57:03 nbsp I tried the chicken teriyaki here for the first time about a month ago  it was pretty good.  Id return if I was in the area and was craving some teriyaki. Users/AlexPomeranz

20070228 13:06:52 nbsp I think that this place is closed. Everything was taken down when i tired to eat here on feb. 28th Users/DonGibson

20070818 08:09:14 nbsp do you know what is opening up in it now.  there has been a lot of construction. Users/MattHh

20070909 13:26:32 nbsp It looks like its being replaced by Burgers and Brew owned and operated by Crepeville (or so the sign on the front says). Wonder how this will compete with the new 3rd & U Cafe seems like the same kinda food is going to be offered at both locations. Users/ArmyDecoy

20111017 14:02:42 nbsp Ate here so many times from 9498, great marinade on the chicken Users/Jeffrey9

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
c---------------------------------------------------------------------
c---------------------------------------------------------------------

      subroutine x_solve

c---------------------------------------------------------------------
c---------------------------------------------------------------------

c---------------------------------------------------------------------
c
c     Performs line solves in X direction by first factoring
c     the block-tridiagonal matrix into an upper triangular matrix,
c     and then performing back substitution to solve for the unknow
c     vectors of each line.
c
c     Make sure we treat elements zero to cell_size in the direction
c     of the sweep.
c
c---------------------------------------------------------------------

      include 'header.h'

      integer i,j,k,m,n,isize

c---------------------------------------------------------------------
c---------------------------------------------------------------------

      if (timeron) call timer_start(t_xsolve)

c---------------------------------------------------------------------
c---------------------------------------------------------------------

c---------------------------------------------------------------------
c     This function computes the left hand side in the xi-direction
c---------------------------------------------------------------------

      isize = grid_points(1)-1

c---------------------------------------------------------------------
c     determine a (labeled f) and n jacobians
c---------------------------------------------------------------------
!$omp parallel do default(shared) shared(isize)
!$omp& private(i,j,k,m,n)
      do k = 1, grid_points(3)-2
         do j = 1, grid_points(2)-2
            do i = 0, isize

               tmp1 = rho_i(i,j,k)
               tmp2 = tmp1 * tmp1
               tmp3 = tmp1 * tmp2
c---------------------------------------------------------------------
c
c---------------------------------------------------------------------
               fjac(1,1,i) = 0.0d+00
               fjac(1,2,i) = 1.0d+00
               fjac(1,3,i) = 0.0d+00
               fjac(1,4,i) = 0.0d+00
               fjac(1,5,i) = 0.0d+00

               fjac(2,1,i) = -(u(2,i,j,k) * tmp2 *
     >              u(2,i,j,k))
     >              + c2 * qs(i,j,k)
               fjac(2,2,i) = ( 2.0d+00 - c2 )
     >              * ( u(2,i,j,k) / u(1,i,j,k) )
               fjac(2,3,i) = - c2 * ( u(3,i,j,k) * tmp1 )
               fjac(2,4,i) = - c2 * ( u(4,i,j,k) * tmp1 )
               fjac(2,5,i) = c2

               fjac(3,1,i) = - ( u(2,i,j,k)*u(3,i,j,k) ) * tmp2
               fjac(3,2,i) = u(3,i,j,k) * tmp1
               fjac(3,3,i) = u(2,i,j,k) * tmp1
               fjac(3,4,i) = 0.0d+00
               fjac(3,5,i) = 0.0d+00

               fjac(4,1,i) = - ( u(2,i,j,k)*u(4,i,j,k) ) * tmp2
               fjac(4,2,i) = u(4,i,j,k) * tmp1
               fjac(4,3,i) = 0.0d+00
               fjac(4,4,i) = u(2,i,j,k) * tmp1
               fjac(4,5,i) = 0.0d+00

               fjac(5,1,i) = ( c2 * 2.0d0 * square(i,j,k)
     >              - c1 * u(5,i,j,k) )
     >              * ( u(2,i,j,k) * tmp2 )
               fjac(5,2,i) = c1 *  u(5,i,j,k) * tmp1
     >              - c2
     >              * ( u(2,i,j,k)*u(2,i,j,k) * tmp2
     >              + qs(i,j,k) )
               fjac(5,3,i) = - c2 * ( u(3,i,j,k)*u(2,i,j,k) )
     >              * tmp2
               fjac(5,4,i) = - c2 * ( u(4,i,j,k)*u(2,i,j,k) )
     >              * tmp2
               fjac(5,5,i) = c1 * ( u(2,i,j,k) * tmp1 )

               njac(1,1,i) = 0.0d+00
               njac(1,2,i) = 0.0d+00
               njac(1,3,i) = 0.0d+00
               njac(1,4,i) = 0.0d+00
               njac(1,5,i) = 0.0d+00

               njac(2,1,i) = - con43 * c3c4 * tmp2 * u(2,i,j,k)
               njac(2,2,i) =   con43 * c3c4 * tmp1
               njac(2,3,i) =   0.0d+00
               njac(2,4,i) =   0.0d+00
               njac(2,5,i) =   0.0d+00

               njac(3,1,i) = - c3c4 * tmp2 * u(3,i,j,k)
               njac(3,2,i) =   0.0d+00
               njac(3,3,i) =   c3c4 * tmp1
               njac(3,4,i) =   0.0d+00
               njac(3,5,i) =   0.0d+00

               njac(4,1,i) = - c3c4 * tmp2 * u(4,i,j,k)
               njac(4,2,i) =   0.0d+00
               njac(4,3,i) =   0.0d+00
               njac(4,4,i) =   c3c4 * tmp1
               njac(4,5,i) =   0.0d+00

               njac(5,1,i) = - ( con43 * c3c4
     >              - c1345 ) * tmp3 * (u(2,i,j,k)**2)
     >              - ( c3c4 - c1345 ) * tmp3 * (u(3,i,j,k)**2)
     >              - ( c3c4 - c1345 ) * tmp3 * (u(4,i,j,k)**2)
     >              - c1345 * tmp2 * u(5,i,j,k)

               njac(5,2,i) = ( con43 * c3c4
     >              - c1345 ) * tmp2 * u(2,i,j,k)
               njac(5,3,i) = ( c3c4 - c1345 ) * tmp2 * u(3,i,j,k)
               njac(5,4,i) = ( c3c4 - c1345 ) * tmp2 * u(4,i,j,k)
               njac(5,5,i) = ( c1345 ) * tmp1

            enddo
c---------------------------------------------------------------------
c     now jacobians set, so form left hand side in x direction
c---------------------------------------------------------------------
            call lhsinit(lhs, isize)
            do i = 1, isize-1

               tmp1 = dt * tx1
               tmp2 = dt * tx2

               lhs(1,1,aa,i) = - tmp2 * fjac(1,1,i-1)
     >              - tmp1 * njac(1,1,i-1)
     >              - tmp1 * dx1
               lhs(1,2,aa,i) = - tmp2 * fjac(1,2,i-1)
     >              - tmp1 * njac(1,2,i-1)
               lhs(1,3,aa,i) = - tmp2 * fjac(1,3,i-1)
     >              - tmp1 * njac(1,3,i-1)
               lhs(1,4,aa,i) = - tmp2 * fjac(1,4,i-1)
     >              - tmp1 * njac(1,4,i-1)
               lhs(1,5,aa,i) = - tmp2 * fjac(1,5,i-1)
     >              - tmp1 * njac(1,5,i-1)

               lhs(2,1,aa,i) = - tmp2 * fjac(2,1,i-1)
     >              - tmp1 * njac(2,1,i-1)
               lhs(2,2,aa,i) = - tmp2 * fjac(2,2,i-1)
     >              - tmp1 * njac(2,2,i-1)
     >              - tmp1 * dx2
               lhs(2,3,aa,i) = - tmp2 * fjac(2,3,i-1)
     >              - tmp1 * njac(2,3,i-1)
               lhs(2,4,aa,i) = - tmp2 * fjac(2,4,i-1)
     >              - tmp1 * njac(2,4,i-1)
               lhs(2,5,aa,i) = - tmp2 * fjac(2,5,i-1)
     >              - tmp1 * njac(2,5,i-1)

               lhs(3,1,aa,i) = - tmp2 * fjac(3,1,i-1)
     >              - tmp1 * njac(3,1,i-1)
               lhs(3,2,aa,i) = - tmp2 * fjac(3,2,i-1)
     >              - tmp1 * njac(3,2,i-1)
               lhs(3,3,aa,i) = - tmp2 * fjac(3,3,i-1)
     >              - tmp1 * njac(3,3,i-1)
     >              - tmp1 * dx3
               lhs(3,4,aa,i) = - tmp2 * fjac(3,4,i-1)
     >              - tmp1 * njac(3,4,i-1)
               lhs(3,5,aa,i) = - tmp2 * fjac(3,5,i-1)
     >              - tmp1 * njac(3,5,i-1)

               lhs(4,1,aa,i) = - tmp2 * fjac(4,1,i-1)
     >              - tmp1 * njac(4,1,i-1)
               lhs(4,2,aa,i) = - tmp2 * fjac(4,2,i-1)
     >              - tmp1 * njac(4,2,i-1)
               lhs(4,3,aa,i) = - tmp2 * fjac(4,3,i-1)
     >              - tmp1 * njac(4,3,i-1)
               lhs(4,4,aa,i) = - tmp2 * fjac(4,4,i-1)
     >              - tmp1 * njac(4,4,i-1)
     >              - tmp1 * dx4
               lhs(4,5,aa,i) = - tmp2 * fjac(4,5,i-1)
     >              - tmp1 * njac(4,5,i-1)

               lhs(5,1,aa,i) = - tmp2 * fjac(5,1,i-1)
     >              - tmp1 * njac(5,1,i-1)
               lhs(5,2,aa,i) = - tmp2 * fjac(5,2,i-1)
     >              - tmp1 * njac(5,2,i-1)
               lhs(5,3,aa,i) = - tmp2 * fjac(5,3,i-1)
     >              - tmp1 * njac(5,3,i-1)
               lhs(5,4,aa,i) = - tmp2 * fjac(5,4,i-1)
     >              - tmp1 * njac(5,4,i-1)
               lhs(5,5,aa,i) = - tmp2 * fjac(5,5,i-1)
     >              - tmp1 * njac(5,5,i-1)
     >              - tmp1 * dx5

               lhs(1,1,bb,i) = 1.0d+00
     >              + tmp1 * 2.0d+00 * njac(1,1,i)
     >              + tmp1 * 2.0d+00 * dx1
               lhs(1,2,bb,i) = tmp1 * 2.0d+00 * njac(1,2,i)
               lhs(1,3,bb,i) = tmp1 * 2.0d+00 * njac(1,3,i)
               lhs(1,4,bb,i) = tmp1 * 2.0d+00 * njac(1,4,i)
               lhs(1,5,bb,i) = tmp1 * 2.0d+00 * njac(1,5,i)

               lhs(2,1,bb,i) = tmp1 * 2.0d+00 * njac(2,1,i)
               lhs(2,2,bb,i) = 1.0d+00
     >              + tmp1 * 2.0d+00 * njac(2,2,i)
     >              + tmp1 * 2.0d+00 * dx2
               lhs(2,3,bb,i) = tmp1 * 2.0d+00 * njac(2,3,i)
               lhs(2,4,bb,i) = tmp1 * 2.0d+00 * njac(2,4,i)
               lhs(2,5,bb,i) = tmp1 * 2.0d+00 * njac(2,5,i)

               lhs(3,1,bb,i) = tmp1 * 2.0d+00 * njac(3,1,i)
               lhs(3,2,bb,i) = tmp1 * 2.0d+00 * njac(3,2,i)
               lhs(3,3,bb,i) = 1.0d+00
     >              + tmp1 * 2.0d+00 * njac(3,3,i)
     >              + tmp1 * 2.0d+00 * dx3
               lhs(3,4,bb,i) = tmp1 * 2.0d+00 * njac(3,4,i)
               lhs(3,5,bb,i) = tmp1 * 2.0d+00 * njac(3,5,i)

               lhs(4,1,bb,i) = tmp1 * 2.0d+00 * njac(4,1,i)
               lhs(4,2,bb,i) = tmp1 * 2.0d+00 * njac(4,2,i)
               lhs(4,3,bb,i) = tmp1 * 2.0d+00 * njac(4,3,i)
               lhs(4,4,bb,i) = 1.0d+00
     >              + tmp1 * 2.0d+00 * njac(4,4,i)
     >              + tmp1 * 2.0d+00 * dx4
               lhs(4,5,bb,i) = tmp1 * 2.0d+00 * njac(4,5,i)

               lhs(5,1,bb,i) = tmp1 * 2.0d+00 * njac(5,1,i)
               lhs(5,2,bb,i) = tmp1 * 2.0d+00 * njac(5,2,i)
               lhs(5,3,bb,i) = tmp1 * 2.0d+00 * njac(5,3,i)
               lhs(5,4,bb,i) = tmp1 * 2.0d+00 * njac(5,4,i)
               lhs(5,5,bb,i) = 1.0d+00
     >              + tmp1 * 2.0d+00 * njac(5,5,i)
     >              + tmp1 * 2.0d+00 * dx5

               lhs(1,1,cc,i) =  tmp2 * fjac(1,1,i+1)
     >              - tmp1 * njac(1,1,i+1)
     >              - tmp1 * dx1
               lhs(1,2,cc,i) =  tmp2 * fjac(1,2,i+1)
     >              - tmp1 * njac(1,2,i+1)
               lhs(1,3,cc,i) =  tmp2 * fjac(1,3,i+1)
     >              - tmp1 * njac(1,3,i+1)
               lhs(1,4,cc,i) =  tmp2 * fjac(1,4,i+1)
     >              - tmp1 * njac(1,4,i+1)
               lhs(1,5,cc,i) =  tmp2 * fjac(1,5,i+1)
     >              - tmp1 * njac(1,5,i+1)

               lhs(2,1,cc,i) =  tmp2 * fjac(2,1,i+1)
     >              - tmp1 * njac(2,1,i+1)
               lhs(2,2,cc,i) =  tmp2 * fjac(2,2,i+1)
     >              - tmp1 * njac(2,2,i+1)
     >              - tmp1 * dx2
               lhs(2,3,cc,i) =  tmp2 * fjac(2,3,i+1)
     >              - tmp1 * njac(2,3,i+1)
               lhs(2,4,cc,i) =  tmp2 * fjac(2,4,i+1)
     >              - tmp1 * njac(2,4,i+1)
               lhs(2,5,cc,i) =  tmp2 * fjac(2,5,i+1)
     >              - tmp1 * njac(2,5,i+1)

               lhs(3,1,cc,i) =  tmp2 * fjac(3,1,i+1)
     >              - tmp1 * njac(3,1,i+1)
               lhs(3,2,cc,i) =  tmp2 * fjac(3,2,i+1)
     >              - tmp1 * njac(3,2,i+1)
               lhs(3,3,cc,i) =  tmp2 * fjac(3,3,i+1)
     >              - tmp1 * njac(3,3,i+1)
     >              - tmp1 * dx3
               lhs(3,4,cc,i) =  tmp2 * fjac(3,4,i+1)
     >              - tmp1 * njac(3,4,i+1)
               lhs(3,5,cc,i) =  tmp2 * fjac(3,5,i+1)
     >              - tmp1 * njac(3,5,i+1)

               lhs(4,1,cc,i) =  tmp2 * fjac(4,1,i+1)
     >              - tmp1 * njac(4,1,i+1)
               lhs(4,2,cc,i) =  tmp2 * fjac(4,2,i+1)
     >              - tmp1 * njac(4,2,i+1)
               lhs(4,3,cc,i) =  tmp2 * fjac(4,3,i+1)
     >              - tmp1 * njac(4,3,i+1)
               lhs(4,4,cc,i) =  tmp2 * fjac(4,4,i+1)
     >              - tmp1 * njac(4,4,i+1)
     >              - tmp1 * dx4
               lhs(4,5,cc,i) =  tmp2 * fjac(4,5,i+1)
     >              - tmp1 * njac(4,5,i+1)

               lhs(5,1,cc,i) =  tmp2 * fjac(5,1,i+1)
     >              - tmp1 * njac(5,1,i+1)
               lhs(5,2,cc,i) =  tmp2 * fjac(5,2,i+1)
     >              - tmp1 * njac(5,2,i+1)
               lhs(5,3,cc,i) =  tmp2 * fjac(5,3,i+1)
     >              - tmp1 * njac(5,3,i+1)
               lhs(5,4,cc,i) =  tmp2 * fjac(5,4,i+1)
     >              - tmp1 * njac(5,4,i+1)
               lhs(5,5,cc,i) =  tmp2 * fjac(5,5,i+1)
     >              - tmp1 * njac(5,5,i+1)
     >              - tmp1 * dx5

            enddo

c---------------------------------------------------------------------
c---------------------------------------------------------------------

c---------------------------------------------------------------------
c     performs guaussian elimination on this cell.
c
c     assumes that unpacking routines for non-first cells
c     preload C' and rhs' from previous cell.
c
c     assumed send happens outside this routine, but that
c     c'(IMAX) and rhs'(IMAX) will be sent to next cell
c---------------------------------------------------------------------

c---------------------------------------------------------------------
c     outer most do loops - sweeping in i direction
c---------------------------------------------------------------------

c---------------------------------------------------------------------
c     multiply c(0,j,k) by b_inverse and copy back to c
c     multiply rhs(0) by b_inverse(0) and copy to rhs
c---------------------------------------------------------------------
            call binvcrhs( lhs(1,1,bb,0),
     >                        lhs(1,1,cc,0),
     >                        rhs(1,0,j,k) )

c---------------------------------------------------------------------
c     begin inner most do loop
c     do all the elements of the cell unless last
c---------------------------------------------------------------------
            do i=1,isize-1

c---------------------------------------------------------------------
c     rhs(i) = rhs(i) - A*rhs(i-1)
c---------------------------------------------------------------------
               call matvec_sub(lhs(1,1,aa,i),
     >                         rhs(1,i-1,j,k),rhs(1,i,j,k))

c---------------------------------------------------------------------
c     B(i) = B(i) - C(i-1)*A(i)
c---------------------------------------------------------------------
               call matmul_sub(lhs(1,1,aa,i),
     >                         lhs(1,1,cc,i-1),
     >                         lhs(1,1,bb,i))

c---------------------------------------------------------------------
c     multiply c(i,j,k) by b_inverse and copy back to c
c     multiply rhs(1,j,k) by b_inverse(1,j,k) and copy to rhs
c---------------------------------------------------------------------
               call binvcrhs( lhs(1,1,bb,i),
     >                        lhs(1,1,cc,i),
     >                        rhs(1,i,j,k) )

            enddo

c---------------------------------------------------------------------
c     rhs(isize) = rhs(isize) - A*rhs(isize-1)
c---------------------------------------------------------------------
            call matvec_sub(lhs(1,1,aa,isize),
     >                         rhs(1,isize-1,j,k),rhs(1,isize,j,k))

c---------------------------------------------------------------------
c     B(isize) = B(isize) - C(isize-1)*A(isize)
c---------------------------------------------------------------------
            call matmul_sub(lhs(1,1,aa,isize),
     >                         lhs(1,1,cc,isize-1),
     >                         lhs(1,1,bb,isize))

c---------------------------------------------------------------------
c     multiply rhs() by b_inverse() and copy to rhs
c---------------------------------------------------------------------
            call binvrhs( lhs(1,1,bb,isize),
     >                       rhs(1,isize,j,k) )

c---------------------------------------------------------------------
c     back solve: if last cell, then generate U(isize)=rhs(isize)
c     else assume U(isize) is loaded in un pack backsub_info
c     so just use it
c     after call u(istart) will be sent to next cell
c---------------------------------------------------------------------

            do i=isize-1,0,-1
               do m=1,BLOCK_SIZE
                  do n=1,BLOCK_SIZE
                     rhs(m,i,j,k) = rhs(m,i,j,k)
     >                    - lhs(m,n,cc,i)*rhs(n,i+1,j,k)
                  enddo
               enddo
            enddo

         enddo
      enddo
      if (timeron) call timer_stop(t_xsolve)

      return
      end

c---------------------------------------------------------------------
c---------------------------------------------------------------------

      subroutine matvec_sub(ablock,avec,bvec)

c---------------------------------------------------------------------
c---------------------------------------------------------------------

c---------------------------------------------------------------------
c     subtracts bvec=bvec - ablock*avec
c---------------------------------------------------------------------

      implicit none

      double precision ablock,avec,bvec
      dimension ablock(5,5),avec(5),bvec(5)
      integer i

      do i=1,5
c---------------------------------------------------------------------
c            rhs(i,ic,jc,kc) = rhs(i,ic,jc,kc)
c     $           - lhs(i,1,ablock,ia)*
c---------------------------------------------------------------------
         bvec(i) = bvec(i) - ablock(i,1)*avec(1)
     >                     - ablock(i,2)*avec(2)
     >                     - ablock(i,3)*avec(3)
     >                     - ablock(i,4)*avec(4)
     >                     - ablock(i,5)*avec(5)
      enddo

      return
      end

c---------------------------------------------------------------------
c---------------------------------------------------------------------

      subroutine matmul_sub(ablock, bblock, cblock)

c---------------------------------------------------------------------
c---------------------------------------------------------------------

c---------------------------------------------------------------------
c     subtracts a(i,j,k) X b(i,j,k) from c(i,j,k)
c---------------------------------------------------------------------

      implicit none

      double precision ablock, bblock, cblock
      dimension ablock(5,5), bblock(5,5), cblock(5,5)
      integer j

      do j=1,5
         cblock(1,j) = cblock(1,j) - ablock(1,1)*bblock(1,j)
     >                             - ablock(1,2)*bblock(2,j)
     >                             - ablock(1,3)*bblock(3,j)
     >                             - ablock(1,4)*bblock(4,j)
     >                             - ablock(1,5)*bblock(5,j)
         cblock(2,j) = cblock(2,j) - ablock(2,1)*bblock(1,j)
     >                             - ablock(2,2)*bblock(2,j)
     >                             - ablock(2,3)*bblock(3,j)
     >                             - ablock(2,4)*bblock(4,j)
     >                             - ablock(2,5)*bblock(5,j)
         cblock(3,j) = cblock(3,j) - ablock(3,1)*bblock(1,j)
     >                             - ablock(3,2)*bblock(2,j)
     >                             - ablock(3,3)*bblock(3,j)
     >                             - ablock(3,4)*bblock(4,j)
     >                             - ablock(3,5)*bblock(5,j)
         cblock(4,j) = cblock(4,j) - ablock(4,1)*bblock(1,j)
     >                             - ablock(4,2)*bblock(2,j)
     >                             - ablock(4,3)*bblock(3,j)
     >                             - ablock(4,4)*bblock(4,j)
     >                             - ablock(4,5)*bblock(5,j)
         cblock(5,j) = cblock(5,j) - ablock(5,1)*bblock(1,j)
     >                             - ablock(5,2)*bblock(2,j)
     >                             - ablock(5,3)*bblock(3,j)
     >                             - ablock(5,4)*bblock(4,j)
     >                             - ablock(5,5)*bblock(5,j)
      enddo

      return
      end

c---------------------------------------------------------------------
c---------------------------------------------------------------------

      subroutine binvcrhs( lhs,c,r )

c---------------------------------------------------------------------
c---------------------------------------------------------------------

c---------------------------------------------------------------------
c
c---------------------------------------------------------------------

      implicit none

      double precision pivot, coeff, lhs
      dimension lhs(5,5)
      double precision c(5,5), r(5)

c---------------------------------------------------------------------
c
c---------------------------------------------------------------------

      pivot = 1.00d0/lhs(1,1)
      lhs(1,2) = lhs(1,2)*pivot
      lhs(1,3) = lhs(1,3)*pivot
      lhs(1,4) = lhs(1,4)*pivot
      lhs(1,5) = lhs(1,5)*pivot
      c(1,1) = c(1,1)*pivot
      c(1,2) = c(1,2)*pivot
      c(1,3) = c(1,3)*pivot
      c(1,4) = c(1,4)*pivot
      c(1,5) = c(1,5)*pivot
      r(1)   = r(1)  *pivot

      coeff = lhs(2,1)
      lhs(2,2)= lhs(2,2) - coeff*lhs(1,2)
      lhs(2,3)= lhs(2,3) - coeff*lhs(1,3)
      lhs(2,4)= lhs(2,4) - coeff*lhs(1,4)
      lhs(2,5)= lhs(2,5) - coeff*lhs(1,5)
      c(2,1) = c(2,1) - coeff*c(1,1)
      c(2,2) = c(2,2) - coeff*c(1,2)
      c(2,3) = c(2,3) - coeff*c(1,3)
      c(2,4) = c(2,4) - coeff*c(1,4)
      c(2,5) = c(2,5) - coeff*c(1,5)
      r(2)   = r(2)   - coeff*r(1)

      coeff = lhs(3,1)
      lhs(3,2)= lhs(3,2) - coeff*lhs(1,2)
      lhs(3,3)= lhs(3,3) - coeff*lhs(1,3)
      lhs(3,4)= lhs(3,4) - coeff*lhs(1,4)
      lhs(3,5)= lhs(3,5) - coeff*lhs(1,5)
      c(3,1) = c(3,1) - coeff*c(1,1)
      c(3,2) = c(3,2) - coeff*c(1,2)
      c(3,3) = c(3,3) - coeff*c(1,3)
      c(3,4) = c(3,4) - coeff*c(1,4)
      c(3,5) = c(3,5) - coeff*c(1,5)
      r(3)   = r(3)   - coeff*r(1)

      coeff = lhs(4,1)
      lhs(4,2)= lhs(4,2) - coeff*lhs(1,2)
      lhs(4,3)= lhs(4,3) - coeff*lhs(1,3)
      lhs(4,4)= lhs(4,4) - coeff*lhs(1,4)
      lhs(4,5)= lhs(4,5) - coeff*lhs(1,5)
      c(4,1) = c(4,1) - coeff*c(1,1)
      c(4,2) = c(4,2) - coeff*c(1,2)
      c(4,3) = c(4,3) - coeff*c(1,3)
      c(4,4) = c(4,4) - coeff*c(1,4)
      c(4,5) = c(4,5) - coeff*c(1,5)
      r(4)   = r(4)   - coeff*r(1)

      coeff = lhs(5,1)
      lhs(5,2)= lhs(5,2) - coeff*lhs(1,2)
      lhs(5,3)= lhs(5,3) - coeff*lhs(1,3)
      lhs(5,4)= lhs(5,4) - coeff*lhs(1,4)
      lhs(5,5)= lhs(5,5) - coeff*lhs(1,5)
      c(5,1) = c(5,1) - coeff*c(1,1)
      c(5,2) = c(5,2) - coeff*c(1,2)
      c(5,3) = c(5,3) - coeff*c(1,3)
      c(5,4) = c(5,4) - coeff*c(1,4)
      c(5,5) = c(5,5) - coeff*c(1,5)
      r(5)   = r(5)   - coeff*r(1)

      pivot = 1.00d0/lhs(2,2)
      lhs(2,3) = lhs(2,3)*pivot
      lhs(2,4) = lhs(2,4)*pivot
      lhs(2,5) = lhs(2,5)*pivot
      c(2,1) = c(2,1)*pivot
      c(2,2) = c(2,2)*pivot
      c(2,3) = c(2,3)*pivot
      c(2,4) = c(2,4)*pivot
      c(2,5) = c(2,5)*pivot
      r(2)   = r(2)  *pivot

      coeff = lhs(1,2)
      lhs(1,3)= lhs(1,3) - coeff*lhs(2,3)
      lhs(1,4)= lhs(1,4) - coeff*lhs(2,4)
      lhs(1,5)= lhs(1,5) - coeff*lhs(2,5)
      c(1,1) = c(1,1) - coeff*c(2,1)
      c(1,2) = c(1,2) - coeff*c(2,2)
      c(1,3) = c(1,3) - coeff*c(2,3)
      c(1,4) = c(1,4) - coeff*c(2,4)
      c(1,5) = c(1,5) - coeff*c(2,5)
      r(1)   = r(1)   - coeff*r(2)

      coeff = lhs(3,2)
      lhs(3,3)= lhs(3,3) - coeff*lhs(2,3)
      lhs(3,4)= lhs(3,4) - coeff*lhs(2,4)
      lhs(3,5)= lhs(3,5) - coeff*lhs(2,5)
      c(3,1) = c(3,1) - coeff*c(2,1)
      c(3,2) = c(3,2) - coeff*c(2,2)
      c(3,3) = c(3,3) - coeff*c(2,3)
      c(3,4) = c(3,4) - coeff*c(2,4)
      c(3,5) = c(3,5) - coeff*c(2,5)
      r(3)   = r(3)   - coeff*r(2)

      coeff = lhs(4,2)
      lhs(4,3)= lhs(4,3) - coeff*lhs(2,3)
      lhs(4,4)= lhs(4,4) - coeff*lhs(2,4)
      lhs(4,5)= lhs(4,5) - coeff*lhs(2,5)
      c(4,1) = c(4,1) - coeff*c(2,1)
      c(4,2) = c(4,2) - coeff*c(2,2)
      c(4,3) = c(4,3) - coeff*c(2,3)
      c(4,4) = c(4,4) - coeff*c(2,4)
      c(4,5) = c(4,5) - coeff*c(2,5)
      r(4)   = r(4)   - coeff*r(2)

      coeff = lhs(5,2)
      lhs(5,3)= lhs(5,3) - coeff*lhs(2,3)
      lhs(5,4)= lhs(5,4) - coeff*lhs(2,4)
      lhs(5,5)= lhs(5,5) - coeff*lhs(2,5)
      c(5,1) = c(5,1) - coeff*c(2,1)
      c(5,2) = c(5,2) - coeff*c(2,2)
      c(5,3) = c(5,3) - coeff*c(2,3)
      c(5,4) = c(5,4) - coeff*c(2,4)
      c(5,5) = c(5,5) - coeff*c(2,5)
      r(5)   = r(5)   - coeff*r(2)

      pivot = 1.00d0/lhs(3,3)
      lhs(3,4) = lhs(3,4)*pivot
      lhs(3,5) = lhs(3,5)*pivot
      c(3,1) = c(3,1)*pivot
      c(3,2) = c(3,2)*pivot
      c(3,3) = c(3,3)*pivot
      c(3,4) = c(3,4)*pivot
      c(3,5) = c(3,5)*pivot
      r(3)   = r(3)  *pivot

      coeff = lhs(1,3)
      lhs(1,4)= lhs(1,4) - coeff*lhs(3,4)
      lhs(1,5)= lhs(1,5) - coeff*lhs(3,5)
      c(1,1) = c(1,1) - coeff*c(3,1)
      c(1,2) = c(1,2) - coeff*c(3,2)
      c(1,3) = c(1,3) - coeff*c(3,3)
      c(1,4) = c(1,4) - coeff*c(3,4)
      c(1,5) = c(1,5) - coeff*c(3,5)
      r(1)   = r(1)   - coeff*r(3)

      coeff = lhs(2,3)
      lhs(2,4)= lhs(2,4) - coeff*lhs(3,4)
      lhs(2,5)= lhs(2,5) - coeff*lhs(3,5)
      c(2,1) = c(2,1) - coeff*c(3,1)
      c(2,2) = c(2,2) - coeff*c(3,2)
      c(2,3) = c(2,3) - coeff*c(3,3)
      c(2,4) = c(2,4) - coeff*c(3,4)
      c(2,5) = c(2,5) - coeff*c(3,5)
      r(2)   = r(2)   - coeff*r(3)

      coeff = lhs(4,3)
      lhs(4,4)= lhs(4,4) - coeff*lhs(3,4)
      lhs(4,5)= lhs(4,5) - coeff*lhs(3,5)
      c(4,1) = c(4,1) - coeff*c(3,1)
      c(4,2) = c(4,2) - coeff*c(3,2)
      c(4,3) = c(4,3) - coeff*c(3,3)
      c(4,4) = c(4,4) - coeff*c(3,4)
      c(4,5) = c(4,5) - coeff*c(3,5)
      r(4)   = r(4)   - coeff*r(3)

      coeff = lhs(5,3)
      lhs(5,4)= lhs(5,4) - coeff*lhs(3,4)
      lhs(5,5)= lhs(5,5) - coeff*lhs(3,5)
      c(5,1) = c(5,1) - coeff*c(3,1)
      c(5,2) = c(5,2) - coeff*c(3,2)
      c(5,3) = c(5,3) - coeff*c(3,3)
      c(5,4) = c(5,4) - coeff*c(3,4)
      c(5,5) = c(5,5) - coeff*c(3,5)
      r(5)   = r(5)   - coeff*r(3)

      pivot = 1.00d0/lhs(4,4)
      lhs(4,5) = lhs(4,5)*pivot
      c(4,1) = c(4,1)*pivot
      c(4,2) = c(4,2)*pivot
      c(4,3) = c(4,3)*pivot
      c(4,4) = c(4,4)*pivot
      c(4,5) = c(4,5)*pivot
      r(4)   = r(4)  *pivot

      coeff = lhs(1,4)
      lhs(1,5)= lhs(1,5) - coeff*lhs(4,5)
      c(1,1) = c(1,1) - coeff*c(4,1)
      c(1,2) = c(1,2) - coeff*c(4,2)
      c(1,3) = c(1,3) - coeff*c(4,3)
      c(1,4) = c(1,4) - coeff*c(4,4)
      c(1,5) = c(1,5) - coeff*c(4,5)
      r(1)   = r(1)   - coeff*r(4)

      coeff = lhs(2,4)
      lhs(2,5)= lhs(2,5) - coeff*lhs(4,5)
      c(2,1) = c(2,1) - coeff*c(4,1)
      c(2,2) = c(2,2) - coeff*c(4,2)
      c(2,3) = c(2,3) - coeff*c(4,3)
      c(2,4) = c(2,4) - coeff*c(4,4)
      c(2,5) = c(2,5) - coeff*c(4,5)
      r(2)   = r(2)   - coeff*r(4)

      coeff = lhs(3,4)
      lhs(3,5)= lhs(3,5) - coeff*lhs(4,5)
      c(3,1) = c(3,1) - coeff*c(4,1)
      c(3,2) = c(3,2) - coeff*c(4,2)
      c(3,3) = c(3,3) - coeff*c(4,3)
      c(3,4) = c(3,4) - coeff*c(4,4)
      c(3,5) = c(3,5) - coeff*c(4,5)
      r(3)   = r(3)   - coeff*r(4)

      coeff = lhs(5,4)
      lhs(5,5)= lhs(5,5) - coeff*lhs(4,5)
      c(5,1) = c(5,1) - coeff*c(4,1)
      c(5,2) = c(5,2) - coeff*c(4,2)
      c(5,3) = c(5,3) - coeff*c(4,3)
      c(5,4) = c(5,4) - coeff*c(4,4)
      c(5,5) = c(5,5) - coeff*c(4,5)
      r(5)   = r(5)   - coeff*r(4)

      pivot = 1.00d0/lhs(5,5)
      c(5,1) = c(5,1)*pivot
      c(5,2) = c(5,2)*pivot
      c(5,3) = c(5,3)*pivot
      c(5,4) = c(5,4)*pivot
      c(5,5) = c(5,5)*pivot
      r(5)   = r(5)  *pivot

      coeff = lhs(1,5)
      c(1,1) = c(1,1) - coeff*c(5,1)
      c(1,2) = c(1,2) - coeff*c(5,2)
      c(1,3) = c(1,3) - coeff*c(5,3)
      c(1,4) = c(1,4) - coeff*c(5,4)
      c(1,5) = c(1,5) - coeff*c(5,5)
      r(1)   = r(1)   - coeff*r(5)

      coeff = lhs(2,5)
      c(2,1) = c(2,1) - coeff*c(5,1)
      c(2,2) = c(2,2) - coeff*c(5,2)
      c(2,3) = c(2,3) - coeff*c(5,3)
      c(2,4) = c(2,4) - coeff*c(5,4)
      c(2,5) = c(2,5) - coeff*c(5,5)
      r(2)   = r(2)   - coeff*r(5)

      coeff = lhs(3,5)
      c(3,1) = c(3,1) - coeff*c(5,1)
      c(3,2) = c(3,2) - coeff*c(5,2)
      c(3,3) = c(3,3) - coeff*c(5,3)
      c(3,4) = c(3,4) - coeff*c(5,4)
      c(3,5) = c(3,5) - coeff*c(5,5)
      r(3)   = r(3)   - coeff*r(5)

      coeff = lhs(4,5)
      c(4,1) = c(4,1) - coeff*c(5,1)
      c(4,2) = c(4,2) - coeff*c(5,2)
      c(4,3) = c(4,3) - coeff*c(5,3)
      c(4,4) = c(4,4) - coeff*c(5,4)
      c(4,5) = c(4,5) - coeff*c(5,5)
      r(4)   = r(4)   - coeff*r(5)

      return
      end

c---------------------------------------------------------------------
c---------------------------------------------------------------------

      subroutine binvrhs( lhs,r )

c---------------------------------------------------------------------
c---------------------------------------------------------------------

c---------------------------------------------------------------------
c
c---------------------------------------------------------------------

      implicit none

      double precision pivot, coeff, lhs
      dimension lhs(5,5)
      double precision r(5)

c---------------------------------------------------------------------
c
c---------------------------------------------------------------------

      pivot = 1.00d0/lhs(1,1)
      lhs(1,2) = lhs(1,2)*pivot
      lhs(1,3) = lhs(1,3)*pivot
      lhs(1,4) = lhs(1,4)*pivot
      lhs(1,5) = lhs(1,5)*pivot
      r(1)   = r(1)  *pivot

      coeff = lhs(2,1)
      lhs(2,2)= lhs(2,2) - coeff*lhs(1,2)
      lhs(2,3)= lhs(2,3) - coeff*lhs(1,3)
      lhs(2,4)= lhs(2,4) - coeff*lhs(1,4)
      lhs(2,5)= lhs(2,5) - coeff*lhs(1,5)
      r(2)   = r(2)   - coeff*r(1)

      coeff = lhs(3,1)
      lhs(3,2)= lhs(3,2) - coeff*lhs(1,2)
      lhs(3,3)= lhs(3,3) - coeff*lhs(1,3)
      lhs(3,4)= lhs(3,4) - coeff*lhs(1,4)
      lhs(3,5)= lhs(3,5) - coeff*lhs(1,5)
      r(3)   = r(3)   - coeff*r(1)

      coeff = lhs(4,1)
      lhs(4,2)= lhs(4,2) - coeff*lhs(1,2)
      lhs(4,3)= lhs(4,3) - coeff*lhs(1,3)
      lhs(4,4)= lhs(4,4) - coeff*lhs(1,4)
      lhs(4,5)= lhs(4,5) - coeff*lhs(1,5)
      r(4)   = r(4)   - coeff*r(1)

      coeff = lhs(5,1)
      lhs(5,2)= lhs(5,2) - coeff*lhs(1,2)
      lhs(5,3)= lhs(5,3) - coeff*lhs(1,3)
      lhs(5,4)= lhs(5,4) - coeff*lhs(1,4)
      lhs(5,5)= lhs(5,5) - coeff*lhs(1,5)
      r(5)   = r(5)   - coeff*r(1)

      pivot = 1.00d0/lhs(2,2)
      lhs(2,3) = lhs(2,3)*pivot
      lhs(2,4) = lhs(2,4)*pivot
      lhs(2,5) = lhs(2,5)*pivot
      r(2)   = r(2)  *pivot

      coeff = lhs(1,2)
      lhs(1,3)= lhs(1,3) - coeff*lhs(2,3)
      lhs(1,4)= lhs(1,4) - coeff*lhs(2,4)
      lhs(1,5)= lhs(1,5) - coeff*lhs(2,5)
      r(1)   = r(1)   - coeff*r(2)

      coeff = lhs(3,2)
      lhs(3,3)= lhs(3,3) - coeff*lhs(2,3)
      lhs(3,4)= lhs(3,4) - coeff*lhs(2,4)
      lhs(3,5)= lhs(3,5) - coeff*lhs(2,5)
      r(3)   = r(3)   - coeff*r(2)

      coeff = lhs(4,2)
      lhs(4,3)= lhs(4,3) - coeff*lhs(2,3)
      lhs(4,4)= lhs(4,4) - coeff*lhs(2,4)
      lhs(4,5)= lhs(4,5) - coeff*lhs(2,5)
      r(4)   = r(4)   - coeff*r(2)

      coeff = lhs(5,2)
      lhs(5,3)= lhs(5,3) - coeff*lhs(2,3)
      lhs(5,4)= lhs(5,4) - coeff*lhs(2,4)
      lhs(5,5)= lhs(5,5) - coeff*lhs(2,5)
      r(5)   = r(5)   - coeff*r(2)

      pivot = 1.00d0/lhs(3,3)
      lhs(3,4) = lhs(3,4)*pivot
      lhs(3,5) = lhs(3,5)*pivot
      r(3)   = r(3)  *pivot

      coeff = lhs(1,3)
      lhs(1,4)= lhs(1,4) - coeff*lhs(3,4)
      lhs(1,5)= lhs(1,5) - coeff*lhs(3,5)
      r(1)   = r(1)   - coeff*r(3)

      coeff = lhs(2,3)
      lhs(2,4)= lhs(2,4) - coeff*lhs(3,4)
      lhs(2,5)= lhs(2,5) - coeff*lhs(3,5)
      r(2)   = r(2)   - coeff*r(3)

      coeff = lhs(4,3)
      lhs(4,4)= lhs(4,4) - coeff*lhs(3,4)
      lhs(4,5)= lhs(4,5) - coeff*lhs(3,5)
      r(4)   = r(4)   - coeff*r(3)

      coeff = lhs(5,3)
      lhs(5,4)= lhs(5,4) - coeff*lhs(3,4)
      lhs(5,5)= lhs(5,5) - coeff*lhs(3,5)
      r(5)   = r(5)   - coeff*r(3)

      pivot = 1.00d0/lhs(4,4)
      lhs(4,5) = lhs(4,5)*pivot
      r(4)   = r(4)  *pivot

      coeff = lhs(1,4)
      lhs(1,5)= lhs(1,5) - coeff*lhs(4,5)
      r(1)   = r(1)   - coeff*r(4)

      coeff = lhs(2,4)
      lhs(2,5)= lhs(2,5) - coeff*lhs(4,5)
      r(2)   = r(2)   - coeff*r(4)

      coeff = lhs(3,4)
      lhs(3,5)= lhs(3,5) - coeff*lhs(4,5)
      r(3)   = r(3)   - coeff*r(4)

      coeff = lhs(5,4)
      lhs(5,5)= lhs(5,5) - coeff*lhs(4,5)
      r(5)   = r(5)   - coeff*r(4)

      pivot = 1.00d0/lhs(5,5)
      r(5)   = r(5)  *pivot

      coeff = lhs(1,5)
      r(1)   = r(1)   - coeff*r(5)

      coeff = lhs(2,5)
      r(2)   = r(2)   - coeff*r(5)

      coeff = lhs(3,5)
      r(3)   = r(3)   - coeff*r(5)

      coeff = lhs(4,5)
      r(4)   = r(4)   - coeff*r(5)

      return
      end


!----------------------------------------------------------------------
!----------------------------------------------------------------------
!   int :    An ODE with boundary and integral constraints
!----------------------------------------------------------------------
!----------------------------------------------------------------------

      SUBROUTINE FUNC(NDIM,U,ICP,PAR,IJAC,F,DFDU,DFDP)
!     ---------- ----

      IMPLICIT NONE
      INTEGER, INTENT(IN) :: NDIM, ICP(*), IJAC
      DOUBLE PRECISION, INTENT(IN) :: U(NDIM), PAR(*)
      DOUBLE PRECISION, INTENT(OUT) :: F(NDIM)
      DOUBLE PRECISION, INTENT(INOUT) :: DFDU(NDIM,NDIM), DFDP(NDIM,*)

      DOUBLE PRECISION E

       E=EXP(U(1))
       F(1)=U(2)
       F(2)=-PAR(1)*E

      IF(IJAC.EQ.0)RETURN

       DFDU(1,1)=0.0
       DFDU(1,2)=1
       DFDU(2,1)=-PAR(1)*E
       DFDU(2,2)=0.0

      IF(IJAC.EQ.1)RETURN

!      *Parameter derivatives
       DFDP(1,1)=0.0
       DFDP(2,1)=-E

      END SUBROUTINE FUNC

      SUBROUTINE STPNT(NDIM,U,PAR,T)
!     ---------- -----

      IMPLICIT NONE
      INTEGER, INTENT(IN) :: NDIM
      DOUBLE PRECISION, INTENT(INOUT) :: U(NDIM),PAR(*)
      DOUBLE PRECISION, INTENT(IN) :: T

       PAR(1)=0
       PAR(2)=0
       PAR(3)=0

       U(1)=0.0
       U(2)=0.0

      END SUBROUTINE STPNT

      SUBROUTINE BCND(NDIM,PAR,ICP,NBC,U0,U1,FB,IJAC,DBC)
!     ---------- ----

      IMPLICIT NONE
      INTEGER, INTENT(IN) :: NDIM, ICP(*), NBC, IJAC
      DOUBLE PRECISION, INTENT(IN) :: PAR(*), U0(NDIM), U1(NDIM)
      DOUBLE PRECISION, INTENT(OUT) :: FB(NBC)
      DOUBLE PRECISION, INTENT(INOUT) :: DBC(NBC,*)

       FB(1)=U0(1)-U1(1)-PAR(2)

      IF(IJAC.EQ.0)RETURN

       DBC(1,1)=1.0
       DBC(1,2)=0.0

       DBC(1,3)=-1.0
       DBC(1,4)=0.0

      IF(IJAC.EQ.1)RETURN

!      *Parameter derivatives
       DBC(1,5)=0.0
       DBC(1,6)=-1.0
       DBC(1,7)=0.0

      END SUBROUTINE BCND

      SUBROUTINE ICND(NDIM,PAR,ICP,NINT,U,UOLD,UDOT,UPOLD,FI,IJAC,DINT)
!     ---------- ----

      IMPLICIT NONE
      INTEGER, INTENT(IN) :: NDIM, ICP(*), NINT, IJAC
      DOUBLE PRECISION, INTENT(IN) :: PAR(*)
      DOUBLE PRECISION, INTENT(IN) :: U(NDIM), UOLD(NDIM), UDOT(NDIM), UPOLD(NDIM)
      DOUBLE PRECISION, INTENT(OUT) :: FI(NINT)
      DOUBLE PRECISION, INTENT(INOUT) :: DINT(NINT,*)

       FI(1)=U(1)-PAR(3)

      IF(IJAC.EQ.0)RETURN

       DINT(1,1)=1.0
       DINT(1,2)=0.0

      IF(IJAC.EQ.1)RETURN

!      *Parameter derivatives
       DINT(1,3)=0.0
       DINT(1,4)=0.0
       DINT(1,5)=-1.0

      END SUBROUTINE ICND

      SUBROUTINE FOPT
      END SUBROUTINE FOPT

      SUBROUTINE PVLS
      END SUBROUTINE PVLS


!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||

 module forcing_tools

!BOP
! !MODULE: forcing_tools
! !DESCRIPTION:
! Contains common variables, parameters and routines that are used by
! all of the individual forcing modules.
!
! !REVISION HISTORY:
! CVS:$Id: forcing_tools.F90,v 1.10 2003/01/10 00:11:25 jfd Exp $
! CVS:$Name: POP_2_0_1 $
!
! !USES:

   use kinds_mod
   use blocks
   use domain
   use constants
   use io
   use io_types
   use grid
   use time_management
   use exit_mod

   implicit none
   private
   save

! !PUBLIC MEMBER FUNCTIONS:

   public :: find_interp_time, &
             get_forcing_filename, &
             find_forcing_times, &
             interpolate_forcing, &
             update_forcing_data, &
             echo_forcing_options

! !PUBLIC DATA MEMBERS:

   real (r8), parameter, public :: &
      never = 1.e20_r8, &! value to signify never
      always = -1.e20_r8 ! value to signify always

   !*** common data shared by forcing routines

!EOP
!BOC
!-----------------------------------------------------------------------
!
! internal module variables
!
!-----------------------------------------------------------------------

   real (r8), dimension(12) :: &
      thour00_midmonth, &
      thour00_endmonth

!EOC
!***********************************************************************

 contains

!***********************************************************************
!BOP
! !IROUTINE: find_interp_time
! !INTERFACE:

 subroutine find_interp_time(forcing_interp_inc, forcing_interp_next)

! !DESCRIPTION:
! Finds the next time at which temporal interpolation is needed
! in the forcing.
!
! !REVISION HISTORY:
! same as module

! !INPUT PARAMETERS:

   real (r8), intent(in) :: &
      forcing_interp_inc ! increment (hours) for intepolation time

! !OUTPUT PARAMETERS:

   real (r8), intent(out) :: &
      forcing_interp_next ! next interpolation time

!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------

   integer (int_kind) :: nstart, n

   real (r8) :: forcing_time_test

!-----------------------------------------------------------------------
!
! find the first interpolation time that is greater than the current
! time, which then becomes the next interpolation time.
!
!-----------------------------------------------------------------------

   nstart = int(thour00/forcing_interp_inc) - 2

   do n = nstart, nstart + 5
      forcing_time_test = (n + p5)*forcing_interp_inc
      if (forcing_time_test > thour00) exit
   enddo

   forcing_interp_next = forcing_time_test

!-----------------------------------------------------------------------
!EOC

 end subroutine find_interp_time

!***********************************************************************
!BOP
! !IROUTINE: get_forcing_filename
! !INTERFACE:

 subroutine get_forcing_filename(forcing_filename, forcing_infile, &
                                 forcing_time, forcing_data_inc)

! !DESCRIPTION:
! Finds the name of an n-hour file given the forcing time needed.
! Files can be labeled by the date of the beginning of the
! forcing interval or by the date of the middle in the forcing
! interval (which corresponds to the value of the variable
! forcing\_time) depending on the value of
! filename\_at\_begin\_of\_interval (T or F).
!
! File names are asumed to be of the form:
! fileroot.year.day.hour (fileroot.yyyy.ddd.hh)
! where fileroot is specified by namelist input, year is
! 4 digits, day is 3 digits, and hour is 2 digits.
! for example, Jan 1 of year 22 with a root of 'ws' would be
! ws.0022.001.00 (labeled by beginning of interval)
! ws.0022.001.12 (labeled by middle of 1-day interval)
!
! Remember, forcing times are relative to 01-01-0000 so results
! may not be what you expect. for example, if the forcing
! increment is 2 days, then even numbered years will have
! mid-interval-referenced files that are on even numbered
! days, and odd numbered years will reference odd numbered
! days (assuming no leap years).
!
! !REVISION HISTORY:
! same as module

! !INPUT PARAMETERS:

   real (r8), intent(in) :: &
      forcing_time, &! time at begin of forcing interval
      forcing_data_inc ! increment in hours to read data

   character (*), intent(in) :: &
      forcing_infile ! root filename

! !OUTPUT PARAMETERS:

   character (char_len), intent(out) :: &
      forcing_filename ! output forcing filename

!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------

   integer (int_kind) :: &
      file_year, &! year for file name
      file_day, &! day for file name
      file_hour, &! hour for file name
      n, &! dummy loop index
      days_in_prev_year, &! days in previous year
      cindx ! index for character strings

   real (r8) :: &
      file_time, &!
      file_time00 !

   character (len = 5) :: file_cyear ! character version of the
   character (len = 4) :: file_cday ! above integer fields
   character (len = 3) :: file_chour

   logical (log_kind) :: &
      filename_at_begin_of_interval ! flag for position of data file
                                    ! in forcing interval (begin,middle)

!-----------------------------------------------------------------------
!
! filename_at_begin_of_interval = .true. for filenames defined at
! the beginning of the forcing interval.
! filename_at_begin_of_interval = .false. for filenames defined at
! the middle of the forcing interval.
!
!-----------------------------------------------------------------------

   filename_at_begin_of_interval = .false.

!-----------------------------------------------------------------------
!
! find the time associated with the beginning or middle of the
! forcing interval.
!
!-----------------------------------------------------------------------

   if (filename_at_begin_of_interval) then
      file_time = forcing_time - p5*forcing_data_inc + eps
   else ! filename defined at forcing time (middle of forcing interval)
      file_time = forcing_time + eps
   endif

!-----------------------------------------------------------------------
!
! if the file time is negative, find the first appropriate
! positive value and use it.
!
!-----------------------------------------------------------------------

   if (file_time < c0 ) then
      if (my_task == master_task) then
         write(stdout,'(a52)') &
            'WARNING: apparently need a forcing file earlier than'
         write(stdout,'(a54)') &
            '01-01-0000 will attempt to find first appropriate file'
      endif
      do n = 1,10
         if (filename_at_begin_of_interval) then
            file_time = forcing_time + (n-p5)*forcing_data_inc + eps
         else
            file_time = forcing_time + n*forcing_data_inc + eps
         endif
         if(file_time >= c0 ) exit
      enddo
   endif
   if (file_time < c0 ) then
      call exit_POP(sigAbort,'get_forcing_filename confused')
   endif

!-----------------------------------------------------------------------
!
! convert file time to time relative to the beginning of the year.
! then find the day and hour corresponding to that time.
!
! note that the use of 'eps' above is important here due to the
! 'int' and 'nint' functions, especially if file_time00 < 0.
!
!-----------------------------------------------------------------------

   thour00_begin_this_year = thour00 - &
                             (elapsed_days_this_year + frac_day)*24.0_r8
   file_time00 = file_time - thour00_begin_this_year
   file_day = int(abs(file_time00)/24.0_r8)
   file_hour = nint(abs(file_time00) - 24*file_day)

!-----------------------------------------------------------------------
!
! if file_time00 < 0, then the time is in the previous year.
!
! if file_time00 > the number of hours in a year, then the time is
! in the next year.
!
!-----------------------------------------------------------------------

   if (file_time00 < c0 )then
      file_hour = 24 - file_hour ! counting backwards
      file_year = iyear - 1
      days_in_prev_year = days_in_norm_year
      if (allow_leapyear .and. mod(file_year,4) == 0) then ! check if previous
         if (mod(file_year,100) /= 0) & ! year was leapyear
            days_in_prev_year = days_in_leap_year
         if (mod(file_year,400) == 0) &
            days_in_prev_year = days_in_leap_year
      endif
      file_day = days_in_prev_year - file_day

   elseif (file_time00 > 24*days_in_year) then
      file_day = file_day - days_in_year + 1
      file_year = iyear + 1
   else
      file_day = file_day + 1 ! day counting starts at 1, not 0
      file_year = iyear
   endif

!-----------------------------------------------------------------------
!
! write year, day, and hour into character variables and use them
! to create filename.
!
!-----------------------------------------------------------------------

   write(file_chour,'(i3)') 100 + file_hour
   write(file_cday, '(i4)') 1000 + file_day
   write(file_cyear,'(i5)') 10000 + file_year

   forcing_filename = char_blank
   cindx = len_trim(forcing_infile)
   forcing_filename(1:cindx) = forcing_infile(1:cindx)
   cindx = cindx + 1
   forcing_filename(cindx:cindx) = '.'
   cindx = cindx + 1
   forcing_filename(cindx:cindx+3) = file_cyear(2:5)
   cindx = cindx + 4
   forcing_filename(cindx:cindx) = '.'
   cindx = cindx + 1
   forcing_filename(cindx:cindx+2) = file_cday(2:4)
   cindx = cindx + 3
   forcing_filename(cindx:cindx) = '.'
   cindx = cindx + 1
   forcing_filename(cindx:cindx+1) = file_chour(2:3)

!-----------------------------------------------------------------------
!EOC

 end subroutine get_forcing_filename

!***********************************************************************
!BOP
! !IROUTINE: find_forcing_times
! !INTERFACE:

 subroutine find_forcing_times(forcing_time, forcing_data_inc, &
                            forcing_interp_type, forcing_time_next, &
                            forcing_time_min_loc, forcing_data_update, &
                            forcing_data_type)

! !DESCRIPTION:
! Find the times associated with forcing data needed for
! interpolation given the current time and time increment
! between forcing data fields. forcing times are assumed
! to be centered; for example, with daily forcing, the value
! is assumed to occur at noon (hour = 12).
!
! !REVISION HISTORY:
! same as module

! !INPUT PARAMETERS:

   character (char_len), intent(in) :: &
      forcing_data_type, &! frequency of forcing data
      forcing_interp_type ! type of interpolation to be used

   real (r8), intent(in) :: &
      forcing_data_inc ! time increment between forcing data

! !INPUT/OUTPUT PARAMETERS:

   real (r8), dimension(12), intent(inout) :: &
      forcing_time

! !OUTPUT PARAMETERS:

   integer (int_kind), intent(out) :: &
      forcing_time_min_loc ! index to first location (in time)
                            ! of forcing data array

   real (r8), intent(out) :: &
      forcing_time_next, &! next time for updating forcing
      forcing_data_update ! next time to update forcing data

!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------

   integer (int_kind) :: &
      nstart, n, &
      second, third, fourth

   real (r8) :: &
      forcing_time_test, &
      forcing_time_min

!-----------------------------------------------------------------------
!
! set mid_month and end_month times for monthly options.
!
!-----------------------------------------------------------------------

   if (forcing_data_type == 'monthly-calendar') then
      thour00_midmonth = thour00_midmonth_calendar
      thour00_endmonth = thour00_endmonth_calendar
   else if (forcing_data_type == 'monthly-equal') then
      thour00_midmonth = thour00_midmonth_equal
      thour00_endmonth = thour00_endmonth_equal
   endif

!-----------------------------------------------------------------------
!
! for 'monthly-equal' and 'monthly-calendar' set the forcing times to
! be the mid-month times plus the time at the beginning of the year.
! then find the first forcing time greater than the current time
! and exit. note that the loop will not exit if the current time
! is in the second half of December in which case n will have the
! value 13.
! for 'n-hour' forcing, start at a time that is definitely less
! than the first forcing time greater than the current time and
! just increment the test time by the forcing increment until we
! find a forcing time greater than the current time.
!
!-----------------------------------------------------------------------

   select case (forcing_data_type)

   case ('monthly-equal','monthly-calendar')

      forcing_time(1:12) = thour00_begin_this_year + &
                           thour00_midmonth(1:12)
      do n = 1,12
         if (forcing_time(n) > thour00) exit
      enddo

   case ('n-hour')

      nstart = int(thour00/forcing_data_inc) - 3
      do n = nstart, nstart + 6
         forcing_time_test = (n - p5)*forcing_data_inc
         if (forcing_time_test > thour00) exit
      enddo

   end select

!-----------------------------------------------------------------------
!
! Set forcing times depending on the type of temporal interpolation.
! Also set forcing_time_min_loc, which is the temporal index of the
! minimum forcing time to be used for the interpolation on the forcing
! data arrays. for 'monthly' data, forcing_time_min_loc varies between
! 1 and 12. For example, if a 4point interpolation is required and
! the months that are needed are 3,4,5,6 then forcing_time_min_loc = 3.
! For 'n-hour' data, forcing_time_min_loc = 1 always because
! initially the data is read chronologically 1 file at a time.
! Also set the times for the next interpolation and next time
! that new data will be needed.
!
!-----------------------------------------------------------------------

   select case(forcing_interp_type)

   case ('nearest') ! select nearest forcing time

      select case (forcing_data_type)

      case ('monthly-equal','monthly-calendar')

         do n = 1,12
            if (thour00 <= &
                thour00_begin_this_year + thour00_endmonth(n)) exit
         enddo
         if (n == 13) call exit_POP(sigAbort, &
            'Error finding current month in find_forcing_times')

         forcing_time_min_loc = n
         forcing_time_min_loc = mod(forcing_time_min_loc,12)
         if (forcing_time_min_loc <= 0) &
            forcing_time_min_loc = forcing_time_min_loc + 12
         if (n == 12) then
            forcing_time_next = forcing_time(1) + hours_in_year ! next Jan
         else
            forcing_time_next = forcing_time(forcing_time_min_loc + 1)
         endif
         forcing_data_update = thour00_begin_this_year + &
                               thour00_endmonth(forcing_time_min_loc)

      case ('n-hour')

         if ((forcing_time_test-thour00) <= p5*forcing_data_inc ) then
            forcing_time_min = (n - p5)*forcing_data_inc
         else
            forcing_time_min = (n - 1.5_r8)*forcing_data_inc
         endif
         forcing_time_min_loc = 1
         forcing_time(forcing_time_min_loc) = forcing_time_min
         forcing_time_next = forcing_time(forcing_time_min_loc) + &
                               forcing_data_inc
         forcing_data_update = forcing_time(forcing_time_min_loc) + &
                               p5*forcing_data_inc

      end select

   case ('linear') ! linear interpolation

      select case(forcing_data_type)

      case ('monthly-equal','monthly-calendar')
         forcing_time_min_loc = n - 1
         forcing_time_min_loc = mod(forcing_time_min_loc,12)
         if (forcing_time_min_loc <= 0) &
            forcing_time_min_loc = forcing_time_min_loc + 12
         second = mod(forcing_time_min_loc+1,12)
         if (second == 0) second = 12
         if (n == 1) then
            forcing_time(forcing_time_min_loc) = &
            forcing_time(forcing_time_min_loc) - hours_in_year ! previous Dec
         elseif (n >= 3) then
            forcing_time(1:n-2) = forcing_time(1:n-2) + hours_in_year
         endif
         if (n == 12) then
            forcing_time_next = forcing_time(1) ! next Jan
         else
            forcing_time_next = forcing_time(second+1) ! next month
         endif

      case ('n-hour')
         forcing_time_min_loc = 1
         second = 2
         forcing_time(forcing_time_min_loc) = (n - 1.5_r8)* &
                                              forcing_data_inc
         forcing_time(second) = forcing_time(forcing_time_min_loc) + &
                                forcing_data_inc
         forcing_time_next = forcing_time(forcing_time_min_loc) + &
                                c2*forcing_data_inc

      end select

      forcing_data_update = forcing_time(second)

   case ('4point') ! 4-point interpolation

      select case(forcing_data_type)

      case ('monthly-equal','monthly-calendar')

         forcing_time_min_loc = n - 2
         forcing_time_min_loc = mod(forcing_time_min_loc,12)
         if (forcing_time_min_loc <= 0) &
            forcing_time_min_loc = forcing_time_min_loc + 12
         second = mod(forcing_time_min_loc+1,12)
         third = mod(forcing_time_min_loc+2,12)
         fourth = mod(forcing_time_min_loc+3,12)
         if (second == 0) second = 12
         if (third == 0) third = 12
         if (fourth == 0) fourth = 12
         if (n == 1) then
            forcing_time(forcing_time_min_loc) = &
            forcing_time(forcing_time_min_loc) - hours_in_year ! previous Nov
            forcing_time(second) = &
            forcing_time(second) - hours_in_year ! previous Dec
         else if (n == 2) then
            forcing_time(forcing_time_min_loc) = &
            forcing_time(forcing_time_min_loc) - hours_in_year ! previous Dec
         else if (n >= 4) then
            forcing_time(1:n-3) = forcing_time(1:n-3) + hours_in_year
         endif
         if (n == 11) then
            forcing_time_next = forcing_time(1) ! next Jan
         else
            forcing_time_next = forcing_time(fourth+1) ! next month
         endif

      case ('n-hour')

         forcing_time_min_loc = 1
         second = 2
         third = 3
         fourth = 4
         forcing_time(forcing_time_min_loc) = (n - 2.5_r8)* &
                                              forcing_data_inc
         forcing_time(second) = forcing_time(forcing_time_min_loc) + &
                                forcing_data_inc
         forcing_time(third) = forcing_time(forcing_time_min_loc) + &
                                c2*forcing_data_inc
         forcing_time(fourth) = forcing_time(forcing_time_min_loc) + &
                                c3*forcing_data_inc
         forcing_time_next = forcing_time(forcing_time_min_loc) + &
                                c4*forcing_data_inc

      end select

      forcing_data_update = forcing_time(third)

   end select

!-----------------------------------------------------------------------
!EOC

 end subroutine find_forcing_times

!***********************************************************************
!BOC
! !IROUTINE: interpolate_forcing
! !INTERFACE:

   subroutine interpolate_forcing(INTERP, FIELD, forcing_time, &
         forcing_interp_type, forcing_time_min_loc, &
         forcing_interp_freq, forcing_interp_inc, forcing_interp_next, &
         forcing_interp_last, nsteps_run_check)

! !DESCRIPTION:
! Temporally interpolates field based on input interpolation options
! and indices.
!
! !REVISION HISTORY:
! same as module

! !INPUT PARAMETERS:

   character (char_len), intent(in) :: &
      forcing_interp_type, &! type of interpolation to use
      forcing_interp_freq ! frequency for interpolating forcing

   integer (int_kind), intent(in) :: &
      nsteps_run_check, &! check for handling first time step
      forcing_time_min_loc ! first location of data to interpolate

   real (r8), dimension(:,:,:,:,:), intent(in) :: &
      FIELD ! forcing field to be interpolated

   real (r8), intent(in) :: &
      forcing_interp_inc, &! forcing increment
      forcing_interp_next

   real (r8), dimension(12), intent(in) :: &
      forcing_time ! times at which forcing located for interpolation

! !INPUT/OUTPUT PARAMETERS:

   real (r8), intent(inout) :: &
      forcing_interp_last ! time when last interpolation done

! !OUTPUT PARAMETERS:

   real (r8), dimension(:,:,:,:), intent(out) :: &
     INTERP ! result of interpolation

!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------

   integer (int_kind) :: &
      k,n, &! loop indices
      iblock, &! dummy block index
      second, third, &! time indices for data to be interpolated
      fourth, &!
      ntime ! number of time slices in data array

   real (r8) :: &
      weight ! linear interpolation weight

!-----------------------------------------------------------------------
!
! if interpolation is done every timestep, then interpolate to
! the current time.
! if interpolation is done every n-hours, then interpolate to
! the current interpolation time (which is >= current time)
! or if it is the first step of the current sub-run, the value
! of forcing_interp_last was read in from the restart file.
!
!-----------------------------------------------------------------------

   ntime = size(FIELD, dim=5)

   if (forcing_interp_freq == 'every-timestep') then
      forcing_interp_last = thour00
   else if (forcing_interp_freq == 'n-hour') then
      if (nsteps_run_check /= 0) forcing_interp_last = &
                                 forcing_interp_next
   else
      call exit_POP(sigAbort,'Error in interpolate_forcing')
   endif

   select case(forcing_interp_type)

   case ('nearest')

      INTERP = FIELD(:,:,:,:,forcing_time_min_loc)

   case ('linear')

      second = mod(forcing_time_min_loc+1,ntime)
      if (second == 0) second = ntime

      weight = (forcing_time(second) - forcing_interp_last)/ &
               (forcing_time(second) - &
                forcing_time(forcing_time_min_loc))

      !*** 3d interior forcing and surface forcing have different
      !*** order of indices so must be treated differently

      if (size(FIELD,dim=3) == km) then ! 3d interior forcing

         !$OMP PARALLEL DO
         do iblock=1,nblocks_clinic
            INTERP(:,:,:,iblock) = &
                    weight *FIELD(:,:,:,iblock,forcing_time_min_loc) + &
              (c1 - weight)*FIELD(:,:,:,iblock,second)
         end do
         !$OMP END PARALLEL DO

      else ! 2d surface forcing with 1 or more fields

         !$OMP PARALLEL DO PRIVATE(iblock, n)
         do iblock=1,nblocks_clinic
         do n=1,size(FIELD,dim=4)
            INTERP(:,:,iblock,n) = &
                   weight *FIELD(:,:,iblock,n,forcing_time_min_loc) + &
             (c1 - weight)*FIELD(:,:,iblock,n,second)
         end do
         end do
         !$OMP END PARALLEL DO

      endif

   case ('4point')

      second = mod(forcing_time_min_loc+1,ntime)
      third = mod(forcing_time_min_loc+2,ntime)
      fourth = mod(forcing_time_min_loc+3,ntime)
      if (second == 0) second = ntime
      if (third == 0) third = ntime
      if (fourth == 0) fourth = ntime

      !*** 3d interior forcing and surface forcing have different
      !*** order of indices so must be treated differently

      if (size(FIELD,dim=3) == km) then ! 3d interior forcing

         !$OMP PARALLEL DO PRIVATE(iblock, k)
         do iblock=1,nblocks_clinic
         do k=1,km

            call interp_4pt(INTERP(:,:,k,iblock), &
                      FIELD(:,:,k,iblock,forcing_time_min_loc), &
                      FIELD(:,:,k,iblock,second), &
                      FIELD(:,:,k,iblock,third), &
                      FIELD(:,:,k,iblock,fourth), &
                      forcing_time(forcing_time_min_loc), &
                      forcing_time(second),forcing_time(third), &
                      forcing_time(fourth),forcing_interp_last)

         end do
         end do
         !$OMP END PARALLEL DO

      else ! 2d surface forcing with 1 or more fields

         !$OMP PARALLEL DO PRIVATE(iblock, n)
         do iblock=1,nblocks_clinic
         do n=1,size(FIELD,dim=4)
            call interp_4pt(INTERP(:,:,iblock,n), &
                      FIELD(:,:,iblock,n,forcing_time_min_loc), &
                      FIELD(:,:,iblock,n,second), &
                      FIELD(:,:,iblock,n,third), &
                      FIELD(:,:,iblock,n,fourth), &
                      forcing_time(forcing_time_min_loc), &
                      forcing_time(second),forcing_time(third), &
                      forcing_time(fourth),forcing_interp_last)

         end do
         end do
         !$OMP END PARALLEL DO
      endif

   end select

!-----------------------------------------------------------------------
!EOC

 end subroutine interpolate_forcing

!***********************************************************************
!BOP
! !IROUTINE:
! !INTERFACE:

 subroutine update_forcing_data(forcing_time, forcing_time_min_loc, &
                            forcing_interp_type, forcing_data_next, &
                            forcing_data_update, forcing_data_type, &
                            forcing_data_inc, FIELD, &
                            forcing_data_rescale, &
                            forcing_data_label, forcing_data_names, &
                            forcing_bndy_loc, forcing_bndy_type, &
                            forcing_infile, forcing_infile_fmt)

! !DESCRIPTION:
! Updates the data to be used in the next temporal interpolation.
! If data is monthly this is a shuffling of indices since
! all 12 months are in memory. if data is n-hour, new data
! needs to be read in as well as shuffling of indices.
!
! !REVISION HISTORY:
! same as module

! !INPUT PARAMETERS:

   character (*), intent(in) :: &
      forcing_infile, &! file containing forcing data
      forcing_infile_fmt, &! format (bin or netcdf) for input file
      forcing_data_type, &! type of forcing data
      forcing_interp_type, &! interpolation to use for forcing
      forcing_data_label ! description of forcing

   character (char_len), dimension(:), intent(in) :: &
      forcing_data_names ! short names for input forcing fields

   integer (int_kind), dimension(:), intent(in) :: &
      forcing_bndy_loc, &! location and field type for ghost
      forcing_bndy_type ! cell updates

   real (r8), dimension(20), intent(in) :: &
      forcing_data_rescale ! factors for converting to model units

   real (r8), intent(in) :: &
      forcing_data_inc ! time increment between forcing data

! !INPUT/OUTPUT PARAMETERS:

   integer (int_kind), intent(inout) :: &
      forcing_time_min_loc ! first temporal index for interpolation

   real (r8), dimension(12), intent(inout) :: &
      forcing_time !

   real (r8), intent(inout) :: &
      forcing_data_next, &! time next data required
      forcing_data_update ! time data to be updated

   real (r8), dimension(:,:,:,:,:), intent(inout) :: &
      FIELD ! forcing field to be read if necessary

!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------

   integer (int_kind) :: &
      n, &! dummy index
      forcing_time_max_loc ! time index for last forcing time

   character (char_len) :: &
      forcing_filename ! name of forcing file

   type (datafile) :: &
      in_forcing_file ! data file type for input forcing file

   type (io_field_desc) :: &
      forcing_field ! io type for input forcing field

   type (io_dim) :: &
      i_dim, j_dim, &! dimension descriptors for horiz dims
      k_dim ! dimension descriptor for vertical levels

!-----------------------------------------------------------------------
!
! determine data by forcing type
!
!-----------------------------------------------------------------------

   select case(forcing_data_type)

!-----------------------------------------------------------------------
!
! for monthly forcing, all twelve months are stored - simply compute
! proper indices into forcing data array
!
!-----------------------------------------------------------------------

   case ('monthly-equal','monthly-calendar')

      !*** compute first forcing time

      forcing_time(forcing_time_min_loc) = &
      forcing_time(forcing_time_min_loc) + hours_in_year
      forcing_time_min_loc = forcing_time_min_loc + 1
      if (forcing_time_min_loc > 12) forcing_time_min_loc = 1

      !*** compute times for other times for use in interpolation

      select case(forcing_interp_type)

      case ('nearest')
         forcing_time_max_loc = forcing_time_min_loc
         forcing_data_update = forcing_data_update + &
                        24.0_r8*days_in_month(forcing_time_max_loc)

      case ('linear')
         forcing_time_max_loc = mod(forcing_time_min_loc + 1,12)
         if (forcing_time_max_loc == 0) forcing_time_max_loc = 12
         forcing_data_update = forcing_time(forcing_time_max_loc)

      case ('4point')
         forcing_time_max_loc = mod(forcing_time_min_loc + 3,12)
         if (forcing_time_max_loc == 0) forcing_time_max_loc = 12
         n = forcing_time_max_loc - 1
         if (n == 0) n = 12
         forcing_data_update = forcing_time(n)

      end select

      !*** compute next forcing time

      forcing_time(forcing_time_max_loc) = forcing_data_next
      n = mod(forcing_time_max_loc + 1, 12)
      if (n == 0) n = 12
      forcing_data_next = forcing_time(n)

!-----------------------------------------------------------------------
!
! for n-hour option, read in required forcing data from forcing file
!
!-----------------------------------------------------------------------

   case ('n-hour')

      call get_forcing_filename(forcing_filename, forcing_infile, &
                                forcing_data_next, forcing_data_inc)

      in_forcing_file = construct_file(forcing_infile_fmt, &
                                    full_name=trim(forcing_filename), &
                                    record_length=rec_type_dbl, &
                                    recl_words=nx_global*ny_global)

      call data_set(in_forcing_file, 'open_read')

      i_dim = construct_io_dim('i',nx_global)
      j_dim = construct_io_dim('j',ny_global)

      if (size(FIELD,dim=3) == km) then ! 3d interior forcing field

         k_dim = construct_io_dim('k',km)

         forcing_field = construct_io_field( &
                         trim(forcing_data_names(1)), &
                         i_dim, j_dim, dim3=k_dim, &
                         field_loc = forcing_bndy_loc(1), &
                         field_type = forcing_bndy_type(1), &
                         d3d_array=FIELD(:,:,:,:,forcing_time_min_loc))

         call data_set (in_forcing_file, 'define', forcing_field)
         call data_set (in_forcing_file, 'read', forcing_field)
         call destroy_io_field(forcing_field)

         !*** renormalize if necessary to compensate for different
         !*** units.

         if (forcing_data_rescale(1) /= c1) &
            FIELD(:,:,:,:,forcing_time_min_loc) = &
            FIELD(:,:,:,:,forcing_time_min_loc)*forcing_data_rescale(1)

      else ! forcing field(s) for surface forcing

         do n = 1,size(FIELD,dim=4)

            forcing_field = construct_io_field( &
                         trim(forcing_data_names(n)), &
                         i_dim, j_dim, &
                         field_loc = forcing_bndy_loc(n), &
                         field_type = forcing_bndy_type(n), &
                         d2d_array=FIELD(:,:,:,n,forcing_time_min_loc))

            call data_set (in_forcing_file, 'define', forcing_field)
            call data_set (in_forcing_file, 'read', forcing_field)
            call destroy_io_field(forcing_field)

            !*** renormalize if necessary to compensate for different
            !*** units.

            if (forcing_data_rescale(n) /= c1) &
               FIELD(:,:,:,n,forcing_time_min_loc) = &
               FIELD(:,:,:,n,forcing_time_min_loc)*forcing_data_rescale(n)

         enddo
      endif

      call data_set (in_forcing_file, 'close')
      call destroy_file(in_forcing_file)

      if (my_task == master_task) then
         write(stdout,blank_fmt)
         write(stdout,'(a9,f12.3,1x,a,a12,a)') 'tday00 = ', tday00, &
            trim(forcing_data_label),' file read: ', &
            trim(forcing_filename)
      endif

      forcing_time(forcing_time_min_loc) = forcing_data_next

      forcing_time_min_loc = forcing_time_min_loc + 1
      if (forcing_time_min_loc > size(FIELD,dim=5)) &
          forcing_time_min_loc = 1
      forcing_data_next = forcing_data_next + forcing_data_inc
      forcing_data_update = forcing_data_update + forcing_data_inc

!-----------------------------------------------------------------------

   end select

!-----------------------------------------------------------------------
!EOC

 end subroutine update_forcing_data

!***********************************************************************
!BOP
! !IROUTINE: echo_forcing_options
! !INTERFACE:

 subroutine echo_forcing_options(forcing_data_type, &
                                 forcing_formulation, &
                                 forcing_data_inc, &
                                 forcing_interp_freq, &
                                 forcing_interp_type, &
                                 forcing_interp_inc, &
                                 forcing_label)

! !DESCRIPTION:
! Writes out forcing options to stdout (or log file).
!
! !REVISION HISTORY:
! same as module

! !INPUT PARAMETERS:

   real (r8), intent(in) :: &
      forcing_data_inc, &
      forcing_interp_inc

   character (char_len), intent(in) :: &
      forcing_data_type, &
      forcing_interp_freq, &
      forcing_interp_type, &
      forcing_label, &
      forcing_formulation

!EOP
!BOC
!-----------------------------------------------------------------------
!
! echo forcing selections to stdout. only write from master task.
!
!-----------------------------------------------------------------------

   if (my_task == master_task) then

!-----------------------------------------------------------------------
!
! no forcing
!
!-----------------------------------------------------------------------

      if (forcing_data_type == 'none') then

         write(stdout,'(a3,a)') 'No ', trim(forcing_label)

!-----------------------------------------------------------------------
!
! if there is some kind of forcing, write out the type (eg, annual)
! and the formulation (eg, restoring).
! then write increment between forcing data if appropriate.
!
!-----------------------------------------------------------------------

      else

         if (len_trim(forcing_formulation) == 0) then
            ! Work around current Linux runtime error
            write(stdout,'(a,a1,a)') trim(forcing_label), ':', &
                                     trim(forcing_data_type)
         else
            write(stdout,'(a,a1,a,a1,a)') trim(forcing_label), ':', &
                                          trim(forcing_data_type),' ', &
                                          trim(forcing_formulation)
         end if

         if (forcing_data_type == 'n-hour') then
            write(stdout,'(a,a29,f7.3)') trim(forcing_label), &
                                  ' Data Time Increment (days): ', &
                                  forcing_data_inc/24.0_r8
         endif

!-----------------------------------------------------------------------
!
! write the frequency that the forcing is updated.
!
!-----------------------------------------------------------------------

         select case(forcing_interp_freq)

         case ('never')
            write(stdout,'(a,a27)') trim(forcing_label), &
                                    ' data is never interpolated'

         case ('n-hour')
            write(stdout,'(a,a25,f7.3,a5)') trim(forcing_label), &
                        ' data interpolated every ', &
                        forcing_interp_inc/24.0_r8,' days'

         case ('every-timestep')
            write(stdout,'(a,a33)') trim(forcing_label), &
                                    ' data interpolated every timestep'

         end select

!-----------------------------------------------------------------------
!
! write order of interpolation if appropriate.
!
!-----------------------------------------------------------------------

         if (forcing_data_type /= 'none' .and. &
             forcing_data_type /= 'analytic' .and. &
             forcing_data_type /= 'annual') &
            write(stdout,'(a,a21,a)') trim(forcing_label), &
              ' Interpolation type: ',trim(forcing_interp_type)

      endif ! forcing type /= none

      write(stdout,blank_fmt)

!-----------------------------------------------------------------------
!
! end of output
!
!-----------------------------------------------------------------------

   endif ! master task

!-----------------------------------------------------------------------
!EOC

 end subroutine echo_forcing_options

!***********************************************************************
!BOP
! !IROUTINE: interp_4pt
! !INTERFACE:

 subroutine interp_4pt(INTERP, DATA1, DATA2, DATA3, DATA4, &
                       tdata1, tdata2, tdata3, tdata4, time)

! !DESCRIPTION:
! Interpolates data at four consecutive (in time) points to a
! particular time using a 3rd degree polynomial fit. This is an
! implementation of Nevilles algorithm as described in Numerical
! Recipes.
!
! !REVISION HISTORY:
! same as module

! !INPUT PARAMETERS:

   real (r8), dimension(nx_block,ny_block), intent(in) :: &
      DATA1, DATA2, DATA3, DATA4 ! data for each grid point in
                                   ! a block for 4 consecutive
                                   ! instants in time

   real (r8), intent(in) :: &
      time, &! time for which interpolated data required
      tdata1, tdata2, &! time points at which above data is located
      tdata3, tdata4

! !OUTPUT PARAMETERS:

   real (r8), dimension(nx_block,ny_block), intent(out) :: &
      INTERP ! data interpolated to input time

!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------

   real (r8), dimension(nx_block,ny_block) :: &
      WORK1,WORK2,WORK3,WORK4 ! local work space

!-----------------------------------------------------------------------

   WORK1 = c0
   WORK2 = c0
   WORK3 = c0
   WORK4 = c0
   call det(WORK1,DATA1,DATA2,tdata1,tdata2,time) ! P12
   call det(WORK2,DATA2,DATA3,tdata2,tdata3,time) ! P23
   call det(WORK3,DATA3,DATA4,tdata3,tdata4,time) ! P34

   call det(WORK4,WORK1,WORK2,tdata1,tdata3,time) ! P123
   call det(WORK1,WORK2,WORK3,tdata2,tdata4,time) ! P234

   call det(INTERP,WORK4,WORK1,tdata1,tdata4,time) ! P1234

!-----------------------------------------------------------------------
!EOC

 end subroutine interp_4pt

!***********************************************************************
!BOP
! !IROUTINE: det
! !INTERFACE:

 subroutine det(C,A,B,y,z,x)

! !DESCRIPTION:
! Determinant routine used by Nevilles algorithm in the 4-point
! time interpolation.
!
! !REVISION HISTORY:
! same as module

! !INPUT PARAMETERS:

   real (r8), intent(in) :: &
      y,z,x

   real (r8), dimension(nx_block,ny_block), intent(in) :: &
      A,B

! !OUTPUT PARAMETERS:

   real (r8), dimension(nx_block,ny_block), intent(out) :: &
      C ! output determinant

!-----------------------------------------------------------------------

   C = ( A*(z-x) - B*(y-x) )/(z-y)

!-----------------------------------------------------------------------
!EOC

 end subroutine det

!***********************************************************************

 end module forcing_tools

!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||


text:
      SUBROUTINE CHSTAB(ITERM)
*
*
*       Chain stability test.
*       ---------------------
*
      INCLUDE 'commonc.h'
      INCLUDE 'common2.h'
      COMMON/CHREG/  TIMEC,TMAX,RMAXC,CM(10),NAMEC(6),NSTEP1,KZ27,KZ30
      COMMON/CPERT/  RGRAV,GPERT,IPERT,NPERT
      REAL*8  M,MB,MB1,R2(NMX,NMX),XCM(3),VCM(3),XX(3,3),VV(3,3),
     &        A1(3),A2(3),XREL(3),VREL(3),EI(3),HI(3),HO(3)
      INTEGER  IJ(NMX)
*
*
*       Sort particle separations (I1 & I2 form closest pair).
      CALL R2SORT(IJ,R2)
      I1 = IJ(1)
      I2 = IJ(2)
      I3 = IJ(3)
      MB = M(I1) + M(I2)
      MB1 = MB + M(I3)
*
*       Form output diagnostics.
      VREL2 = 0.0D0
      VREL21 = 0.0D0
      RDOT = 0.0D0
      RDOT3 = 0.0D0
      R3 = 0.0D0
      DO 5 K = 1,3
          J1 = 3*(I1-1) + K
          J2 = 3*(I2-1) + K
          J3 = 3*(I3-1) + K
          VREL2 = VREL2 + (V(J1) - V(J2))**2
          RDOT = RDOT + (X(J1) - X(J2))*(V(J1) - V(J2))
          XCM(K) = (M(I1)*X(J1) + M(I2)*X(J2))/MB
          VCM(K) = (M(I1)*V(J1) + M(I2)*V(J2))/MB
          R3 = R3 + (X(J3) - XCM(K))**2
          VREL21 = VREL21 + (V(J3) - VCM(K))**2
          RDOT3 = RDOT3 + (X(J3) - XCM(K))*(V(J3) - VCM(K))
          XX(K,1) = X(J1)
          XX(K,2) = X(J2)
          XX(K,3) = X(J3)
          VV(K,1) = V(J1)
          VV(K,2) = V(J2)
          VV(K,3) = V(J3)
    5 CONTINUE
*
*       Evaluate orbital elements for inner and outer motion.
      RB = SQRT(R2(I1,I2))
      R3 = SQRT(R3)
      SEMI = 2.0D0/RB - VREL2/MB
      SEMI = 1.0/SEMI
      ECC = SQRT((1.0D0 - RB/SEMI)**2 + RDOT**2/(SEMI*MB))
      SEMI1 = 2.0/R3 - VREL21/MB1
      SEMI1 = 1.0/SEMI1
      ECC1 = SQRT((1.0D0 - R3/SEMI1)**2 + RDOT3**2/(SEMI1*MB1))
      PMIN = SEMI1*(1.0D0 - ECC1)
*
*       Form hierarchical stability ratio (Eggleton & Kiseleva 1995).
*     QL = MB/M(I3)
*     Q1 = MAX(M(I2)/M(I1),M(I1)/M(I2))
*     Q3 = QL**0.33333
*     Q13 = Q1**0.33333
*     AR = 1.0 + 3.7/Q3 - 2.2/(1.0 + Q3) + 1.4/Q13*(Q3 - 1.0)/(Q3 + 1.0)
*     EK = AR*SEMI*(1.0D0 + ECC)
*
*       Obtain the inclination.
      CALL INCLIN(XX,VV,XCM,VCM,ALPHA)
*
*       Replace the EK criterion by the MA 1999 stability expression.
*     PCRIT = stability(M(I1),M(I2),M(I3),ECC,ECC1,ALPHA)*SEMI
*       Add 1% for perturbation to avoid repeated switching.
*     PCRIT = 1.01*PCRIT
*
*       Evaluate the general stability function (Mardling, Cambody 2008).
      IF (ECC1.LT.1.0.AND.ECC.LT.1.0) THEN
          NST = NSTAB(SEMI,SEMI1,ECC,ECC1,ALPHA,M(I1),M(I2),M(I3))
          IF (NST.EQ.0) THEN
              PCRIT = 0.99*PMIN
          ELSE
              PCRIT = 1.01*PMIN
          END IF
      ELSE
*       Set nominal failed value for SEMI1 > 0 test.
          PCRIT = 1.01*PMIN
      END IF
*
*       Prepare evaluation of maximum eccentricity (see routine HIGROW).
      DO 10 K = 1,3
          XREL(K) = XX(K,1) - XX(K,2)
          VREL(K) = VV(K,1) - VV(K,2)
   10 CONTINUE
      A12 = 0.0
      A22 = 0.0
      A1A2 = 0.0
      RI2 = 0.0
      VI2 = 0.0
      RVI = 0.0
      DO 12 K = 1,3
          K1 = K + 1
          IF (K1.GT.3) K1 = 1
          K2 = K1 + 1
          IF (K2.GT.3) K2 = 1
          A1(K) = XREL(K1)*VREL(K2) - XREL(K2)*VREL(K1)
          A2(K) = (XX(K1,3) - XCM(K1))*(VV(K2,3) - VCM(K2))
     &          - (XX(K2,3) - XCM(K2))*(VV(K1,3) - VCM(K1))
          A12 = A12 + A1(K)**2
          A22 = A22 + A2(K)**2
          A1A2 = A1A2 + A1(K)*A2(K)
          RI2 = RI2 + XREL(K)**2
          VI2 = VI2 + VREL(K)**2
          RVI = RVI + XREL(K)*VREL(K)
   12 CONTINUE
*
*       Construct the Runge-Lenz vector (Heggie & Rasio 1995, Eq.(5)).
      EI2 = 0.0
      DO 15 K = 1,3
          EI(K) = (VI2*XREL(K) - RVI*VREL(K))/MB - XREL(K)/SQRT(RI2)
          EI2 = EI2 + EI(K)**2
   15 CONTINUE
      EI2 = MIN(EI2,0.9999D0)
*
*       Define unit vectors for inner eccentricity and angular momenta.
      COSJ = 0.0
      SJSG = 0.0
      DO 18 K = 1,3
          EI(K) = EI(K)/SQRT(EI2)
          HI(K) = A1(K)/SQRT(A12)
          HO(K) = A2(K)/SQRT(A22)
          COSJ = COSJ + HI(K)*HO(K)
          SJSG = SJSG + EI(K)*HO(K)
   18 CONTINUE
*
*       Evaluate the expressions A & Z.
      A = COSJ*SQRT(1.0 - EI2)
      Z = (1.0 - EI2)*(2.0 - COSJ**2) + 5.0*EI2*SJSG**2
*
*       Obtain maximum inner eccentricity (Douglas Heggie, Sept. 1995).
      Z2 = Z**2 + 25.0 + 16.0*A**4 - 10.0*Z - 20.0*A**2 - 8.0*A**2*Z
      Z2 = MAX(Z2,0.0D0)
      EMAX = (Z + 1.0 - 4.0*A**2 + SQRT(Z2))/6.0
      EMAX = MAX(EMAX,0.0001D0)
      EMAX = SQRT(EMAX)
*
*       Check hierarchical stability condition for bound close pair (RB > a).
      ITERM = 0
      ALPHA = 180.0*ALPHA/3.1415
      IF (PMIN.GT.PCRIT.AND.SEMI.GT.0.0.AND.SEMI1.GT.0.0.AND.
     &    RB.GT.SEMI) THEN
*       Delay termination until perturbation < 10^{-6} (includes CHAIN B-B).
          IF (R3.GT.SEMI1.AND.GPERT.LT.1.0D-06) THEN
          ITERM = -1
          WRITE (6,20)  NAMEC(I1), NAMEC(I2), NAMEC(I3), ECC, EMAX,
     &                  ECC1, SEMI, SEMI1, PMIN, PCRIT, ALPHA
   20     FORMAT (' NEW HIARCH    NM =',3I6,'  E =',F6.3,'  EX =',F7.4,
     &                         '  E1 =',F6.3,'  A =',1P,E8.1,
     &                         '  A1 =',E8.1,'  PM =',E9.2,
     &                         '  PC =',E9.2,'  IN =',0P,F7.1)
          RI = SQRT(CM(1)**2 + CM(2)**2 + CM(3)**2)
          Q0 = M(I3)/MB
          WRITE (81,30)  TIMEC, RI, NAMEC(I3), Q0, ECC, EMAX, ECC1,
     &                   SEMI, SEMI1, PCRIT/PMIN, ALPHA
   30     FORMAT (2F8.4,I6,F6.2,3F6.3,1P,2E10.2,0P,F5.2,F8.1)
          CALL FLUSH(81)
          END IF
*       Include termination test for wide triple system (exclude ECC1 > 0.9).
      ELSE IF (PMIN.GT.3.0*SEMI*(1.0 + ECC).AND.SEMI.GT.0.0.AND.
     &         ECC1.LT.0.9) THEN
*       Wait for favourable configuration (R > SEMI, R3 > SEMI1 & RDOT3 > 0.
          IF (RB.GT.SEMI.AND.R3.GT.SEMI1.AND.RDOT3.GT.0.0) THEN
              APO = SEMI*(1.0 + ECC)
              WRITE (6,40)  ECC, ECC1, ALPHA, RB, R3, PCRIT, PMIN, APO
   40         FORMAT (' WIDE CHAIN    E E1 IN RB R3 PC PM APO ',
     &                                2F7.3,F7.1,1P,5E10.2)
              CALL FLUSH(6)
              ITERM = -1
          END IF
*       Enforce termination of long-lived configuration using two limits.
      ELSE IF (NSTEP1.GT.50000.AND.PMIN.GT.0.9*PCRIT.AND.
     &    R3.GT.SEMI1) THEN
          ITERM = -1
      ELSE IF (NSTEP1.GT.500000.AND.R3.GT.SEMI1) THEN
          ITERM = -1
      END IF
*
      RETURN
*
      END