module GwtGwtExchangeModule
  use KindModule, only: I4B, LGP
  use ConstantsModule, only: LINELENGTH
  use ListModule, only: ListType
  use InputOutputModule, only: getunit, openfile
  use SimModule, only: store_error
  use MemoryManagerModule, only: mem_allocate, mem_deallocate
  use MemoryHelperModule, only: create_mem_path
  use BaseExchangeModule, only: BaseExchangeType, AddBaseExchangeToList
  use DisConnExchangeModule
  use ListsModule, only: baseexchangelist, basemodellist
  use BaseModelModule, only: BaseModelType, GetBaseModelFromList
  use GwtModule, only: GwtModelType

  implicit none
  private

  public :: gwtexchange_create
  public :: GwtExchangeType
  public :: GetGwtExchangeFromList
  public :: CastAsGwtExchange

  !> Data structure to hold the spatial connection data
  !! between two GWT models. The logic to get the coefficients
  !! to the solution will be in the GwtGwtConnection.
  !<
  type, extends(DisConnExchangeType) :: GwtExchangeType

    class(GwtModelType), pointer :: gwtmodel1 !< one transport model in the exchange
    class(GwtModelType), pointer :: gwtmodel2 !< the other transport model in the exchange

    integer(I4B), pointer :: iAdvScheme       !< the advection scheme at the interface:
                                              !! 0 = upstream, 1 = central, 2 = TVD

  contains
    procedure :: exg_df => gwtgwt_df
    procedure :: exg_da => gwtgwt_da
    procedure :: read_options
    procedure :: parse_option
    procedure :: connects_model => gwtgwt_connects_model
    procedure :: use_interface_model
    procedure :: allocate_scalars

  end type GwtExchangeType

contains

!> @brief Create the GWT-GWT exchange data structure
!<
subroutine gwtexchange_create(filename, id, m1id, m2id)
  character(len=*),intent(in) :: filename !< the input file for this exchange
  integer(I4B), intent(in) :: id          !< the numbered id for this exchange
  integer(I4B), intent(in) :: m1id        !< model number one
  integer(I4B), intent(in) :: m2id        !< model number two
  ! -- local
  type(GwtExchangeType), pointer :: exchange
  class(BaseExchangeType), pointer :: baseExchange
  class(BaseModelType), pointer :: mb
  character(len=20) :: cint

  allocate(exchange)
  baseExchange => exchange
  call AddBaseExchangeToList(baseexchangelist, baseExchange)

  ! assign id and name
  exchange%id = id
  write(cint, '(i0)') id
  exchange%name = 'GWT-GWT_' // trim(adjustl(cint))
  exchange%memoryPath = create_mem_path(exchange%name)

  ! -- allocate scalars and set defaults
  call exchange%allocate_scalars()
  exchange%filename = filename
  exchange%typename = 'GWT-GWT'
  !
  ! set models
  mb => GetBaseModelFromList(basemodellist, m1id)
  select type (mb)
  type is (GwtModelType)
    exchange%model1 => mb
    exchange%gwtmodel1 => mb
  end select
  mb => GetBaseModelFromList(basemodellist, m2id)
  select type (mb)
  type is (GwtModelType)
    exchange%model2 => mb
    exchange%gwtmodel2 => mb
  end select

end subroutine gwtexchange_create

!> @brief Return true when this exchange provides matrix
!< coefficients for solving @param model
function gwtgwt_connects_model(this, model) result(is_connected)
  class(GwtExchangeType) :: this                      !< the instance of the GWT-GWT exchange
  class(BaseModelType), pointer, intent(in) :: model  !< the model to which the exchange might hold a connection
  logical(LGP) :: is_connected                        !< true, when connected

  is_connected = .false.

  select type(model)
  class is (GwtModelType)
  if (associated(this%gwtmodel1, model)) then
    is_connected = .true.
  else if (associated(this%gwtmodel2, model)) then
    is_connected = .true.
  end if
  end select

end function gwtgwt_connects_model

!> @brief Should interface model be used for this exchange
!<
function use_interface_model(this) result(useIM)
  class(GwtExchangeType) :: this !< instance of exchange object
  logical(LGP) :: useIM          !< true when interface model should be used

  ! transport always uses IM for this coupling
  useIM = .true.

end function

!> @brief allocate the scalar variables for this exchange
!<
subroutine allocate_scalars(this)
  class(GwtExchangeType) :: this !< this GWT exchange

  call this%DisConnExchangeType%allocate_scalars()
  call mem_allocate(this%iAdvScheme, 'IADVSCHEME', this%memoryPath)

  this%iAdvScheme = 0
  this%ixt3d = 1

end subroutine allocate_scalars

!> @brief read the data from the input file
!< and store
subroutine gwtgwt_df(this)
  use SimVariablesModule, only: iout
  class(GwtExchangeType) :: this
  ! local
  integer(I4B) :: inunit

  inunit = getunit()
  write(iout,'(/a,a)') ' Creating exchange: ', this%name
  call openfile(inunit, iout, this%filename, 'GWT-GWT')

  call this%parser%Initialize(inunit, iout)

  ! check that models are in same solution
  if(this%gwtmodel1%idsoln /= this%gwtmodel2%idsoln) then
    call store_error('ERROR.  TWO MODELS ARE CONNECTED ' //                  &
      'IN A GWT EXCHANGE BUT THEY ARE IN DIFFERENT SOLUTIONS. ' //           &
      'GWT MODELS MUST BE IN SAME SOLUTION: ' //                             &
      trim(this%gwtmodel1%name) // ' ' // trim(this%gwtmodel2%name) )
    call this%parser%StoreErrorUnit()
  endif

  call this%read_options(iout)

  call this%read_dimensions(iout)

  call this%DisConnExchangeType%allocate_arrays()

  call this%read_data(iout)

  ! increase_edge_count equivalent?

  close(inunit)

end subroutine gwtgwt_df

!> @brief Read options block from input file
!<
subroutine read_options(this, iout)
  use ConstantsModule, only: LINELENGTH, LENAUXNAME, DEM6
  use MemoryManagerModule, only: mem_allocate
  use SimVariablesModule, only: errmsg
  class(GwtExchangeType) :: this   !< this exchange
  integer(I4B), intent(in) :: iout !< the output file unit
  ! -- local
  character(len=LINELENGTH) :: keyword
  logical :: isfound
  logical :: endOfBlock
  integer(I4B) :: ierr

  call this%parser%GetBlock('OPTIONS', isfound, ierr,                        &
    supportOpenClose=.true., blockRequired=.false.)
  !
  ! -- parse options block if detected
  if (isfound) then
    write(iout,'(1x,a)')'PROCESSING GWT-GWT EXCHANGE OPTIONS'
    do
      call this%parser%GetNextLine(endOfBlock)
      if (endOfBlock) then
        exit
      end if
      call this%parser%GetStringCaps(keyword)

      ! first parse option in base
      if (this%DisConnExchangeType%parse_option(keyword, iout)) then
        cycle
      end if

      ! it's probably ours
      if (this%parse_option(keyword, iout)) then
        cycle
      end if

      ! unknown option
      errmsg = "Unknown GWT-GWT exchange option '" // trim(keyword) // "'."
      call store_error(errmsg)
      call this%parser%StoreErrorUnit()
    end do

    write(iout,'(1x,a)') 'END OF GWT-GWT EXCHANGE OPTIONS'
  end if

  return
end subroutine read_options

!> @brief parse option from exchange file
  !<
function parse_option(this, keyword, iout) result(parsed)
  use SimVariablesModule, only: errmsg
  class(GwtExchangeType) :: this                   !< instance of exchange object
  character(len=LINELENGTH), intent(in) :: keyword !< the option name
  integer(I4B), intent(in) :: iout                 !< for logging
  logical(LGP) :: parsed                           !< true when parsed
  ! local
  character(len=LINELENGTH) :: subkey

  parsed = .true.

  select case (keyword)
  case ('ADVSCHEME')
    call this%parser%GetStringCaps(subkey)
    select case(subkey)
    case('UPSTREAM')
      this%iAdvScheme = 0
    case('CENTRAL')
      this%iAdvScheme = 1
    case('TVD')
      this%iAdvScheme = 2
    case default
      errmsg = "Unknown weighting method for advection: '" // trim(subkey) // "'."
      call store_error(errmsg)
      call this%parser%StoreErrorUnit()
    end select

    write(iout,'(4x,a,a)')                                                      &
      'CELL AVERAGING METHOD HAS BEEN SET TO: ', trim(subkey)

  case ('XT3D_OFF')
    this%ixt3d = 0
    write(iout, '(4x,a)') 'XT3D FORMULATION HAS BEEN SHUT OFF.'
  case ('XT3D_RHS')
    this%ixt3d = 2
    write(iout, '(4x,a)') 'XT3D RIGHT-HAND SIDE FORMULATION IS SELECTED.'
  case default
    parsed = .false.
  end select

end function parse_option

!> @brief deallocate members of the exchange object
!<
subroutine gwtgwt_da(this)
  class(GwtExchangeType) :: this

  call mem_deallocate(this%iAdvScheme)

  ! deallocate base
  call this%DisConnExchangeType%disconnex_da()

end subroutine gwtgwt_da

function CastAsGwtExchange(obj) result (res)
  class(*), pointer, intent(inout) :: obj
  class(GwtExchangeType), pointer :: res

  res => null()
  if (.not. associated(obj)) return

  select type (obj)
  class is (GwtExchangeType)
    res => obj
  end select
  return
end function CastAsGwtExchange

function GetGwtExchangeFromList(list, idx) result (res)
  type(ListType), intent(inout) :: list
  integer(I4B), intent(in)    :: idx
  class(GwtExchangeType), pointer    :: res
  ! -- local
  class(*), pointer :: obj

  obj => list%GetItem(idx)
  res => CastAsGwtExchange(obj)

  return
end function GetGwtExchangeFromList

end module GwtGwtExchangeModule


text:
      subroutine rngint(ncrop,nowcrp)
c
c     NOTE: BECAUSE "GDAY" IS USED TO DETERMINE WHICH EQUATION IS USED
c           FOR DEFAULT VALUES OF "TLIVE", "TILLAGE" WHICH READS GDAY
c           *MUST* ALWAYS BE CALLED BEFORE "RNGINT".
c
c   ***************************************************************
c   * NOTE: ALL SOURCE FILES THAT REFERENCE "CRINPT1.INC" MUST BE *
c   *       COMPILED WITH THE *NEW* CRINPT1.INC TO WORK WITH THIS *
c   *       FILE.                                                 *
c   ***************************************************************
c
c   ***************************************************************
c   * NOTE: ALL SOURCE FILES THAT REFERENCE "CRINPT3.INC" MUST BE *
c   *       COMPILED WITH THE *NEW* CRINPT3.INC TO WORK WITH THIS *
c   *       FILE.                                                 *
c   ***************************************************************
c
c   ***************************************************************
c   * NOTE: ALL SOURCE FILES THAT REFERENCE "CROUT.INC" MUST BE   *
c   *       COMPILED WITH THE *NEW* CROUT.INC TO WORK WITH THIS   *
c   *       FILE.                                                 *
c   ***************************************************************
c
c     Note: Need definitions of PSCDST & PSCD2S in include file CRINPT2.INC.
c           CRM -- 4/09/93.
c
c     + + + PURPOSE + + +
c     Initializes variables used in RANGE, when land use is rangeland.
c     Called by OFE from CONTIN.
c
c     Called from CONTIN
c     Author(s): Weltz, Meyer
c     Reference in User Guide: Chapter 8
c
c     Changes:
c          1) Deleted references to common block UPDATE.
c          2) Deleted local variables SHEEP & NOWRES.
c          3) Added SPAI, & TPAI to CRINPT3.INC permit single initial
c             calculation of SPAI & TPAI in RNGINT, rather than daily
c             computation in PATRIB.  This required a change in the
c             logic: a loop to initialize each plant type.  Deleted
c             the appropriate code in PATRIB.
c          4) Deleted generic SAVE of all variables.
c          5) Dimensioned rooty(MXNSL) to rooty(MXNSL,MXPLAN) in CROUT.INC.
c             Re: 5/6/93 conversation w/Weltz. -- CRM
c          6) Subscripts were reversed:
c                 if(gday(i,1,iplane) .gt. 0) then
c             changed to:
c                 if(gday(1,i,iplane) .gt. 0) then
c          7) Did a major overhaul relative to the variable PLIVE.  In
c             its original state, PLIVE would not handle more than one
c             OFE correctly since it is modified if burning or herbicide
c             are used.  It also would not supply a default value for a
c             plant that was grazed on one OFE and not grazed on another.
c             Since modules need to be able to swap between plant types
c             on an OFE during a simulation, PLIVE(NTYPE) was redimen-
c             sioned to PLIVE(NTYPE,IPLANE).  This permits an OFE to be
c             altered without affecting the other OFE's.  CRINPT1.INC
c             was altered to accommodate this.  -- CRM -- 5/13/93.
c          8) PROOT was altered analogous to PLIVE. -- CRM -- 5/14/93.
c          9) ROOT(MXPLAN) was removed from common CROUT.INC and made
c             a local variable only - since it is computed and used now
c             only in this routine.  dcf  5/19/93
c         10) statement setting rdf(plant) to 0.43 moved to top of
c             code - otherwise it is undefined at first point of use.
c         11) variable rdf was a constant value of 0.43, and was only used
c             in this subroutine.  Thus it was removed from common CRINPT1.INC
c             and uses in RNGINT were changed to values of 0.43.  dcf  5/21/93
c
c     Version: This module recoded from WEPP Version 92.25.
c     Date recoded: 04/07/93 - 4/09/93.
c                   05/06/93 - 5/13/93.
c     Recoded by: Charles R. Meyer.
c
c     + + + KEYWORDS + + +
c
c     + + + PARAMETERS + + +
      include 'pntype.inc'
      include 'pmxnsl.inc'
      include 'ptilty.inc'
      include 'pmxres.inc'
      include 'pmxpln.inc'
      include 'pmxgrz.inc'
      include 'pmxcrp.inc'
      include 'pmxtls.inc'
      include 'pmxtil.inc'
      include 'pmxhil.inc'
c
c     + + + ARGUMENT DECLARATIONS + + +
      integer ncrop, nowcrp
c
c     + + + ARGUMENT DEFINITIONS + + +
c     ncrop  - number of input crops (rangeland plant communities)
c     nowcrp - current crop
c
c     + + + COMMON BLOCKS + + +
      include 'ccrpgro.inc'
c      modify: ytn
c       write: daymin
c
      include 'ccrpvr1.inc'
c       write: rtm(mxres,mxplan)
c
      include 'ccrpvr2.inc'
c     modify: vdmt(mxplan), cf(ntype)
c
      include 'ccrpvr3.inc'
c      modify: bbb(ntype)
c
      include 'ccover.inc'
c       write: canhgt(mxplan),inrcov(mxplan),rilcov(mxplan)
c      modify: cancov(mxplan),gcover(mxplan)
c
      include 'ccrpout.inc'
c       write: lai(mxplan), rtd(mxplan)
c      modify: rtmass(mxplan),
c
      include 'ccrpprm.inc'
c        read: itype(mxcrop,mxplan)
c
      include 'cends.inc'
c      modify: rspace(mxplan) (?)
c
      include 'crinpt1.inc'
c      modify: ffk(ntype),plive(ntype,mxplan),proot(ntype,mxplan),
c              yield(mxplan)
c       write: pyield(mxplan)
c
c     include 'crinpt1a.inc'
c
      include 'crinpt2.inc'
c
      include 'crinpt3.inc'
c      modify: spai(ntype), tpai(ntype), ghgt(ntype), basden
c
      include 'crinpt4.inc'
c      modify: dwood(mxplan)
c       write: ptlive(mxplan),first(ntype),xlive(mxplan)
c
      include 'crinpt5.inc'
c        read: tdiam(ntype), tcoeff(ntype), tpop(ntype), thgt(ntype)
c
      include 'crout.inc'
c      modify: tlive(mxplan),slive
c       write: droot(mxnsl),rooty(mxnsl,mxplan)
c
      include 'crinpt6.inc'
c      modify: bcover(mxplan), rufcov(mxplan)
c
      include 'cstruc.inc'
c        read: iplane
c
      include 'cwater.inc'
c       write: watstr(mxplan)
c
      include 'cwint.inc'
c        read: deglat
c
c     + + + LOCAL VARIABLES + + +
      real gpai, totpai, tarea, prevrt, ch, h
      real tmpvr2, tmpvr3, tmpvr4, tmpvr5, tmpvr6, root(mxnsl)
      integer plant, i, j
c
c     + + + LOCAL DEFINITIONS + + +
c     gpai   - total projected area from the side, intercepted for grass &
c              herbaceous plants, per 100 meter transect; ie, vertical area
c              "blocked" by plants, measured perpendicular to the slope.
c              Used for snowdrift calculations.
c     totpai - total projected area from the side, for all plants combined
c     tarea  - area of the vertical rectangle circumscribed around the
c              largest plant type present (trees, shrubs, & herbs).
c     prevrt - root mass in the previous layer (one directly above)
c     ch     -
c     h      -
c
c     + + + END SPECIFICATIONS + + +
c
c
      plant = itype(nowcrp,iplane)
c
c     ---- root distribution coefficient for mass by depth
c     rdf(plant) = 0.43
c     NOTE - problems with use of RDF in code below - and since it is
c     always set to a constant value of 0.43 - change from use
c     of an array to a numerical constant in the appropriate
c     equations.    dcf  5/21/93
c
c     ---- Julian date peak standing crop occurs, for first growth peak
      pscdst(plant) = pscday(plant)
c     ---- Julian date peak standing crop occurs, for second growth peak
      pscd2s(plant) = scday2(plant)
c
c     From a user-supplied value of PLIVE for each plant, generate
c     an initial value of PLIVE for each plant and each OFE.  This
c     value may be modified later by RHERB or RBURN, and the plant
c     type may change on this OFE during the simulation.  Also, if
c     the user entered a zero for the current plant, estimate a
c     starting value of PLIVE which depends on whether the OFE is
c     grazed in the first year or not.
c
c     *** L1 IF ***
c     ---- Execute this code the first time it is encountered (OFE #1).
      if (iplane.eq.1) then
c       *** BEGIN L2 DO-LOOP ***
c       Initialize PLIVE values for each plant by OFE.
        do 30 i = 1, ncrop
c         Next line added because NCROP is not readily available to RNGINT.
          if (plive(i,1).ge.0.0) then
c
c           *** L3 IF ***
c           -------- user supplied an initial value
            if (plive(i,1).ne.0.0) then
              do 10 j = 1, mxplan
c               Next line added because NELMT is not readily available to RNGINT.
                if (plive(i,1).gt.0.0) then
c                 ------------ Set PLIVE for each OFE to user-input value.
                  plive(i,j) = plive(i,1)
c                 ------------ maximum potential root biomass
c                 (WEPP Equation 8.5.25)
c                 proot(i,j) = root10(i) / (10.0**rdf(i))
c                 NOTE - assuming that value for rdf=0.43 - using this to simplify equation
c                 dcf 5/21/93
c
                  proot(i,j) = root10(i) / 2.6915348
                end if
   10         continue
c
c           *** L3 ELSE ***
c           -------- model must generate initial value...
            else
c             ---------- set PLIVE for this plant for *all* OFE's
              do 20 j = 1, mxplan
c               Next line added because NELMT is not readily available to RNGINT.
                if (plive(i,1).gt.0.0) then
c                 ---------- check for grazing on this OFE the first year.
                  if (gday(1,i,j).gt.0) then
c                   -------------- grazed areas
c                   (WEPP Equation 8.5.6)
                    plive(i,j) = 0.3569 + 0.0036 * pptg(i)
c                 ------------ non-grazed areas
                  else
c                   (WEPP Equation 8.5.7)
                    plive(i,j) = 0.7723 + 0.0030 * pptg(i)
                  end if
c                 ------------ maximum potential root biomass
c                 (WEPP Equation 8.5.25)
c                 proot(i,j) = root10(i) / (10.0**rdf(i))
c                 NOTE - assuming that value for rdf=0.43 - using this to simplify equation
c                 dcf 5/21/93
c
                  proot(i,j) = root10(i) / 2.6915348
                end if
   20         continue
c           *** L3 ENDIF ***
            end if
          end if
c
c       *** END L2 DO-LOOP ***
   30   continue
c     *** L1 ENDIF ***
      end if
c
c
c     ---- Assume rooting depth (RTD) equals depth to bottom of last soil layer
      rtd(iplane) = solthk(nsl(iplane),iplane)
c
c
c     *******************************************
c     * Calculate root distribution by depth. *
c     *******************************************
c
c     ---- maximum potential root biomass
c     (WEPP Equation 8.5.25)
c     proot(plant) = root10(plant) / (10.0**rdf(plant))
c
      prevrt = 0.0
      do 40 i = 1, nsl(iplane)
c       ------ cumulative total roots to bottom of current layer
c       (WEPP Equation 8.5.24)
        rtmass(iplane) = rootf(plant) * proot(plant,iplane) * (100.0*
     1      solthk(i,iplane)) ** 0.43
c
c       1                         (100.0*solthk(i,iplane))**rdf(plant)
c       NOTE - since RDF is only used in SR RNGINT - removed from common,
c       and since value is always 0.43 for all plant types - use
c       this.  Problems in do 8 loop with values of rdf(i)
c       dcf 5/21/93
c
c       ------ root mass in current soil layer today
        root(i) = rtmass(iplane) - prevrt
c       ------ save root mass to bottom of this layer (top of next layer)
        prevrt = rtmass(iplane)
   40 continue
c
      if (rtmass(iplane).gt.0.0) then
        do 50 i = 1, nsl(iplane)
c         -------- fraction of rootmass in each soil layer (constant throughout sim'n.)
c         (WEPP Equation 8.5.26)
          droot(i) = root(i) / rtmass(iplane)
c         -------- initialize root mass yesterday
          rooty(i,iplane) = root(i)
   50   continue
      else
c       ------ ensure non-negative rootmass
        rtmass(iplane) = 0.0
      end if
c
c     ---- fraction of growth permitted by water stress
      watstr(iplane) = 1.0
c     ---- new plant growth today (min_fraction_live_biomass * max_amt_live_biomass)
      slive = rgcmin(plant) * plive(plant,iplane)
c     ---- total live plant material today
      tlive(iplane) = slive
c     ---- previous day's total live plant material
      ptlive(iplane) = tlive(iplane)
c     ---- total above ground plant production for a simulation year
      yield(iplane) = slive
c     ---- daily net primary plant production
      pyield(iplane) = yield(iplane)
c     ---- previous year's date of initial growth
      first(plant) = strrgc(plant)
c     ----  standing dead biomass left after burning or herbicide
      dwood(iplane) = wood(plant) * rmagt(iplane)
c     ---- evergreen leafy component
      xlive(iplane) = slive
c
c     *************************************************************
c     * NOTE: New equations to calculate  canhgt, cancov, rcov,   *
c     *       rilcov, lai for first day of simulation, should be  *
c     *       added to this section when plant growth subroutine  *
c     *       is updated.                                         *
c     *************************************************************
c
c     ---- total above ground dry weight
      vdmt(iplane) = slive + rmagt(iplane)
c     ---- leaf area index
      lai(iplane) = tlive(iplane) * aleaf(plant)
c
c     --- If this is the first OFE, compute SPAI & TPAI for ALL plant groups.
      if (iplane.eq.1) then
        do 60, i = 1, ncrop
c         -------- projected area for shrubs (remains constant throughout simulation)
c         (WEPP Equation 8.5.17)
          if (sdiam(i).gt.0.0) then
            tmpvr4 = sdiam(i) * scoeff(i) * spop(i)
c           ---------- prevent horizontal measure from exceeding 100 m total
            if (tmpvr4.gt.100.0) tmpvr4 = 100.0
            spai(i) = shgt(i) * tmpvr4
          else
            spai(i) = 0.0
          end if
c
c         -------- projected area for trees (remains constant throughout simulation)
c         (WEPP Equation 8.5.17)
          if (tdiam(i).gt.0.0) then
            tmpvr5 = tdiam(i) * tcoeff(i) * tpop(i)
            if (tmpvr5.gt.100.0) tmpvr5 = 100.0
            tpai(i) = thgt(i) * tmpvr5
          else
            tpai(i) = 0.0
          end if
   60   continue
      end if
c
c     ---- projected area of herbaceous plants (changes daily)
c     (WEPP Equation 8.5.17)
      if (gdiam(plant).gt.0.0) then
        tmpvr6 = gdiam(plant) * gcoeff(plant) * gpop(plant)
        if (tmpvr6.gt.100.0) tmpvr6 = 100.0
        gpai = ghgt(plant) * tmpvr6
      else
        gpai = 0.0
      end if
c
c     **********************************************************************
c     * NOTE: We assume NO OVERLAP between the grass, shrubs, and trees.   *
c     *       Since the current purpose of TOTPAI is solely to compute     *
c     *       a weighted average for canopy height, adding the complexity  *
c     *       of considering overlap did not seem worth the cost.          *
c     *       CRM -- 4/7/93.                                               *
c     **********************************************************************
c     ---- total projected area for vegetation
c     (WEPP Equation 8.5.16)
      totpai = spai(plant) + gpai + tpai(plant)
c
c     ---- vertical area required for largest plant type, using 100 meter transect.
c     (WEPP Equation 8.5.18)
      tarea = max(thgt(plant),shgt(plant),ghgt(plant)) * 100.0
c
c     ---- basal density of plants; ie, the fraction of the "vertical rectangle"
c     (re: TAREA) covered by plants
c     (WEPP Equation 8.5.15)
      if (tarea.gt.0.0) then
        basden = totpai / tarea
      else
        basden = 0.0
      end if
c
c     ---- canopy height, a weighted mean.
c     (WEPP Equation 8.5.13)
      if (basden.gt.0.0) then
        canhgt(iplane) = (ghgt(plant)*gpai+shgt(plant)*spai(plant)+
     1      thgt(plant)*tpai(plant)) / totpai
      else
        canhgt(iplane) = 0.0
      end if
c
c     ... calculate rill spacing
c     (WEPP Equation 8.5.19)
c     rspace(iplane)=100.0/(spop(plant)+gpop(plant)+tpop(plant)+1)
c     if(rspace(iplane) .gt. 5.0) rspace(iplane)=5.0
c     if(rspace(iplane) .lt. 0.5) rspace(iplane)=0.5
c
c     ******************************************************************
c     * NOTE: To calculate grass height requires height coefficient    *
c     *       "bbb", which is estimated from the following regression  *
c     *       equations.                                               *
c     ******************************************************************
c     ---- vegetative dry matter at maturity
      vdmmax(iplane) = plive(plant,iplane) + rmagt(iplane) -
     1    dwood(iplane)
c     ----  parameter value for canopy height equation
      if (bbb(plant).ge.0.0) then
        if (bbb(plant).lt.1.5) then
          bbb(plant) = 2.0 * (vdmmax(iplane)**(-1.0)) + 0.3031455e-06
        else if (bbb(plant).lt.2.5) then
c         Original Code:
c         bbb(plant)=4.0*vdmmax(iplane)**(-1.0)+.23155e-05
          bbb(plant) = 4.0 / vdmmax(iplane) + 0.23155e-05
        else
c         Original Code:
c         bbb(plant)=6.0*vdmmax(iplane)**(-1.0)-.1694262e-04
          bbb(plant) = 6.0 / vdmmax(iplane) - 0.1694262e-04
        end if
      end if
c
c     ---- calculate residue cover from residue mass
c     (WEPP Equation 8.5.20)
c     XXX - commented out by Kidwell/Weltz  3/95   dcf
c     tmpvr2 = 1.085 * (1.0+1.069*exp(-12.61*rmogt(1,iplane))) - 0.5838
c     cf(plant) = 10.0 ** (tmpvr2)
c     (WEPP Equation 8.5.20)
c     rescov(iplane) = 1.0 - exp(-cf(plant)*rmogt(1,iplane))
c
c     **Added by Kidwell 5/25/95**
      if (rescof(iplane) .le. 0.0001) then
        if (rescov(iplane).le.0.0 .or. rmogt(1,iplane).le.0.0) then
                cf(plant) = -63.9
                rescov(iplane) = 0.0
                rmogt(1,iplane) = 0.0
        else
                tmpvr2 = log(1.0 - rescov(iplane))
                cf(plant) = tmpvr2/rmogt(1,iplane)
        end if
      else if (rescof(iplane) .gt. 0.0) then
        if(rescov(iplane).le.0.0 .or. rmogt(1,iplane).le.0.0) then
                rescov(iplane) = 0.0
                rmogt(1,iplane) = 0.0
        else
                cf(plant) = rescof(iplane)
        end if
      end if
c
c     ---- calculate canopy cover from plant mass corresponding to 100% cover (COLD)
c     (WEPP Equation 8.5.22)
c     Original Code:
c     ffk(plant) =21.39148-54.90758*cold(plant)+61.11016*cold(plant)
c     1  **2-30.44471*cold(plant)**3+5.561994*cold(plant)**4
c     XXX - commented out next 3 lines (Kidwell/Weltz) 3/95  dcf
c     tmpvr3 = cold(plant) ** 2
c     ffk(plant) = 21.39148 - 54.90758 * cold(plant) + 61.11016 * tmpvr3
c    1    - 30.44471 * cold(plant) * tmpvr3 + 5.561994 * tmpvr3 ** 2
c     -- XXX -- Need a second equation number in the User Doc.
c     (WEPP Equation 8.5.22)
c     XXX Next 2 lines commented out by Kidwell/Weltz  3/95  dcf
c     cancov(iplane) = 1.0 - exp(-ffk(plant)*vdmt(iplane))
c     if (cancov(iplane).lt.0.0) cancov(iplane) = 0.0
c     **Added by Kidwell 5/25/95**
      if (cancof(iplane) .le. 0.0001) then
        if (cancov(iplane).le.0.0 .or. vdmt(iplane).le.0.0) then
                ffk(plant) = -31.5
                cancov(iplane) = 0.0
                vdmt(iplane) = 0.0
        else
                tmpvr3 = log(1.0 - cancov(iplane))
                ffk(plant) = tmpvr3/vdmt(iplane)
        end if
      else if (cancof(iplane) .gt. 0.0) then
        if(cancov(iplane).le.0.0 .or. vdmt(iplane).le.0.0) then
                cancov(iplane) = 0.0
                vdmt(iplane) = 0.0
        else
                ffk(plant) = cancof(iplane)
        end if
      end if

c
c     -- XXX -- Can't find this equation in the User Doc. -- CRM -- 4/07/93.
c     ---- basal plant cover (ground area occupied by stalks, stems, trunks)
c     XXX Next 5 lines commented out by Kidwell/Weltz  3/95   dcf
c     if (wood(plant).gt.0.0) then
c       bcover(iplane) = 0.3354242 * cancov(iplane)
c     else
c       bcover(iplane) = 0.429118 * cancov(iplane)
c     end if
c
c     -- XXX -- Can't find this equation in the User Doc. -- CRM -- 4/07/93.
c     ---- calculate hydraulic roughness cover from basal cover and moss & lichens
c     XXX Next 2 lines changed by Kidwell/Weltz  3/95   dcf
c     rufcov(iplane) = bcover(iplane) + crypto(iplane)
c     if (rufcov(iplane).lt.0.0) rufcov(iplane) = 0.0
      rufcov(iplane) = bascov(iplane)*fbasr(iplane) +
     1                 crycov(iplane)*fcryr(iplane)
c
c     ---- total ground cover (residue, rocks & gravel, stems and moss & lichens)
c     XXX Next line changed by Kidwell/Weltz  3/95   dcf
c     gcover(iplane) = rescov(iplane) + wcf(iplane) + rufcov(iplane)
      gcover(iplane) = rescov(iplane) + rokcov(iplane) + bascov(iplane)
     1                 + crycov(iplane)
      if (gcover(iplane).gt.1.0) gcover(iplane) = 1.0
c
c     -- XXX -- Need an equation number in the User Doc.
c     Following changed by Mary Kidwell 3/95  dcf
c     rilcov(iplane) = rescov(iplane) + wcf(iplane)
c
c     -- XXX -- Need an equation number in the User Doc.
c     inrcov(iplane) = gcover(iplane)
c
c     --- total rill cover
c
      rilcov(iplane) = rescov(iplane)*fresr(iplane) +
     1                 rokcov(iplane)*frokr(iplane) +
     1                 bascov(iplane)*fbasr(iplane) +
     1                 crycov(iplane)*fcryr(iplane)
c
c     ---- total interrill cover
c
      inrcov(iplane) = rescov(iplane)*fresi(iplane) +
     1                 rokcov(iplane)*froki(iplane) +
     1                 bascov(iplane)*fbasr(iplane) +
     1                 crycov(iplane)*fcryi(iplane)
c
c     ---- initial value of dead root mass for day 1 and when management changes
      rtm(1,iplane) = rtmass(iplane)
c
c     -- XXX -- Is DAYMIN needed here?  It is computed in WATBAL before calling
c     RANGE. -- CRM -- 5/14/93.
c     calculate minimum day length for the input latitude
c
c     original code 10/1/91
c     ytn = tan(deglat / 58.09)
c
c     code changed because of problem with Lahey compiler and
c     tangent function when running on AT&T 6300 machine  dcf
c
c     ---- function of latitude - used in day length calculation (EPIC 2.194)
      ytn = (sin(deglat/58.09)) / (cos(deglat/58.09))
      ch = 0.4349 * abs(ytn)
      if (ch.ge.1.0) then
        h = 0.0
      else
        h = acos(ch)
      end if
c     ---- minimum day length
      daymin = 7.72 * h
c
      return
      end

ext: for

max_stars_repo_name: jarad/dep

max_issues_repo_name: jarad/dep

max_forks_repo_name: akrherz/idep


text: #redirect Department of Music

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


	SUBROUTINE sort(arr,length)
	USE nrtype; USE nrutil, ONLY : swap,nrerror
	IMPLICIT NONE
    integer,intent(in)::length
	REAL(SP), INTENT(INOUT) :: arr(length)
	INTEGER(I4B), PARAMETER :: NN=15, NSTACK=50
	REAL(SP) :: a
	INTEGER(I4B) :: n,k,i,j,jstack,l,r
	INTEGER(I4B), DIMENSION(NSTACK) :: istack
	n=size(arr)
	jstack=0
	l=1
	r=n
	do
		if (r-l < NN) then
			do j=l+1,r
				a=arr(j)
				do i=j-1,l,-1
					if (arr(i) <= a) exit
					arr(i+1)=arr(i)
				end do
				arr(i+1)=a
			end do
			if (jstack == 0) RETURN
			r=istack(jstack)
			l=istack(jstack-1)
			jstack=jstack-2
		else
			k=(l+r)/2
			call swap(arr(k),arr(l+1))
			call swap(arr(l),arr(r),arr(l)>arr(r))
			call swap(arr(l+1),arr(r),arr(l+1)>arr(r))
			call swap(arr(l),arr(l+1),arr(l)>arr(l+1))
			i=l+1
			j=r
			a=arr(l+1)
			do
				do
					i=i+1
					if (arr(i) >= a) exit
				end do
				do
					j=j-1
					if (arr(j) <= a) exit
				end do
				if (j < i) exit
				call swap(arr(i),arr(j))
			end do
			arr(l+1)=arr(j)
			arr(j)=a
			jstack=jstack+2
			if (jstack > NSTACK) call nrerror('sort: NSTACK too small')
			if (r-i+1 >= j-l) then
				istack(jstack)=r
				istack(jstack-1)=i
				r=j-1
			else
				istack(jstack)=j-1
				istack(jstack-1)=l
				l=i
			end if
		end if
	end do
	END SUBROUTINE sort


text:
!! https://gcc.gnu.org/bugzilla/show_bug.cgi?id=79072
!! fixed in 254966 (2017-11-20)
!!
!! $ gfortran --version
!! GNU Fortran (GCC) 7.2.1 20171028
!!
!! $ gfortran gfortran-bug-20170812a.f90
!! $ ./a.out
!! b="" (expect "foo")

character(3), target :: a = 'foo'
class(*), pointer :: b
b => ptr()
select type (b)
type is (character(*))
  print '(3a)', 'b="', b, '" (expect "foo")'
end select
contains
  function ptr()
    class(*), pointer :: ptr
    ptr => a
  end function
end

max_stars_repo_name: dhnza/fortran-compiler-tests

max_issues_repo_name: dhnza/fortran-compiler-tests

max_forks_repo_name: dhnza/fortran-compiler-tests


program test_comm

	use comm_br_mod
	integer::handle
	call comm_br_initialize_object('/users/benabed/Desktop/sigma_dx9_delta_WMAP_n128_v1.fits','/users/benabed/Desktop/test_cls_v4.dat', 2, 32, 1, 20, 20, 100, 1, handle)

end


text:
! { dg-do run }
! pr35940
      program FA1031
      implicit none
      integer I
      INTEGER IDA1(10)
      LOGICAL GDA1(10)
      INTEGER RSLT(10)
      DATA RSLT /4,1,4,1,4,1,4,1,4,1/
      IDA1   = 0
      gda1 = (/ (i/2*2 .ne. I, i=1,10) /)

      IDA1 = INDEX ( 'DEFDEF' , 'DEF', GDA1 )    !fails
      do I = 1, 10
         if (IDA1(i).NE.RSLT(i)) call abort
      end do
      IDA1 = INDEX ( (/ ('DEFDEF',i=1,10) /) , 'DEF', GDA1 )    !works
      do I = 1, 10
         if (IDA1(i).NE.RSLT(i)) call abort
      end do

      END

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
Yoohoo is a Chocolate chocolatey drink that is beloved by many. Its a nice alternative to chocolate milk, although it is not dairy free.

Places in Davis to get Yoohoo.

 Anderson Road Chevron Chevron on Anderson
 Nugget Market in the cold drink section next to specialty soda
 7Eleven 711 carries it

Its a common misconception that Yoohoo is some kind of chemical concoction, but it dates back to the 1920s and the key ingredient is actually a really old and classic one: whey.  As in curds and whey, the stuff Miss Muffet was eating when sitting on her tuffet in the nursery poetry rhyme.  Its also got real cocoa.   So it is basically the stuff left over from the first step of making cheese, sweetened (like most drinks it switched to highfructose corn syrup sometime in the last century), and chocolate flavored.  Since whey is more stable than milk, it can sit on a shelf rather than having to be refrigerated.  Note: even though this drink isnt a chemical enigma, its pretty unhealthy; each serving contains more sugar than the recommended daily maximum.


text:
! SPDX-Identifier: MIT

module test_logicalloc
  use stdlib_array, only : trueloc, falseloc
  use stdlib_kinds, only : dp, i8 => int64
  use stdlib_strings, only : to_string
  use testdrive, only : new_unittest, unittest_type, error_type, check
  implicit none
  private

  public :: collect_logicalloc

contains

  !> Collect all exported unit tests
  subroutine collect_logicalloc(testsuite)
    !> Collection of tests
    type(unittest_type), allocatable, intent(out) :: testsuite(:)

    testsuite = [ &
      new_unittest("trueloc-empty", test_trueloc_empty), &
      new_unittest("trueloc-all", test_trueloc_all), &
      new_unittest("trueloc-where", test_trueloc_where), &
      new_unittest("trueloc-merge", test_trueloc_merge), &
      new_unittest("trueloc-pack", test_trueloc_pack), &
      new_unittest("falseloc-empty", test_falseloc_empty), &
      new_unittest("falseloc-all", test_falseloc_all), &
      new_unittest("falseloc-where", test_falseloc_where), &
      new_unittest("falseloc-merge", test_falseloc_merge), &
      new_unittest("falseloc-pack", test_falseloc_pack) &
      ]
  end subroutine collect_logicalloc

  subroutine test_trueloc_empty(error)
    !> Error handling
    type(error_type), allocatable, intent(out) :: error

    integer :: ndim
    real, allocatable :: avec(:), bvec(:)

    do ndim = 100, 12000, 100
      allocate(avec(ndim))

      call random_number(avec)

      bvec = avec
      bvec(trueloc(bvec < 0)) = 0.0

      call check(error, all(bvec == avec))
      deallocate(avec, bvec)
      if (allocated(error)) exit
    end do
  end subroutine test_trueloc_empty

  subroutine test_trueloc_all(error)
    !> Error handling
    type(error_type), allocatable, intent(out) :: error

    integer :: ndim
    real, allocatable :: avec(:)

    do ndim = 100, 12000, 100
      allocate(avec(-ndim/2:ndim))

      call random_number(avec)

      avec(trueloc(avec > 0, lbound(avec, 1))) = 0.0

      call check(error, all(avec == 0.0))
      deallocate(avec)
      if (allocated(error)) exit
    end do
  end subroutine test_trueloc_all

  subroutine test_trueloc_where(error)
    !> Error handling
    type(error_type), allocatable, intent(out) :: error

    integer :: ndim
    real, allocatable :: avec(:), bvec(:), cvec(:)
    real(dp) :: tl, tw

    tl = 0.0_dp
    tw = 0.0_dp
    do ndim = 100, 12000, 100
      allocate(avec(ndim))

      call random_number(avec)
      avec(:) = avec - 0.5

      bvec = avec
      tl = tl - timing()
      bvec(trueloc(bvec > 0)) = 0.0
      tl = tl + timing()

      cvec = avec
      tw = tw - timing()
      where(cvec > 0) cvec = 0.0
      tw = tw + timing()

      call check(error, all(bvec == cvec))
      deallocate(avec, bvec, cvec)
      if (allocated(error)) exit
    end do
    call report("trueloc", tl, "where", tw)
  end subroutine test_trueloc_where

  subroutine test_trueloc_merge(error)
    !> Error handling
    type(error_type), allocatable, intent(out) :: error

    integer :: ndim
    real, allocatable :: avec(:), bvec(:), cvec(:)
    real(dp) :: tl, tm

    tl = 0.0_dp
    tm = 0.0_dp
    do ndim = 100, 12000, 100
      allocate(avec(ndim))

      call random_number(avec)
      avec(:) = avec - 0.5

      bvec = avec
      tl = tl - timing()
      bvec(trueloc(bvec > 0)) = 0.0
      tl = tl + timing()

      cvec = avec
      tm = tm - timing()
      cvec(:) = merge(0.0, cvec, cvec > 0)
      tm = tm + timing()

      call check(error, all(bvec == cvec))
      deallocate(avec, bvec, cvec)
      if (allocated(error)) exit
    end do
    call report("trueloc", tl, "merge", tm)
  end subroutine test_trueloc_merge

  subroutine test_trueloc_pack(error)
    !> Error handling
    type(error_type), allocatable, intent(out) :: error

    integer :: ndim
    real, allocatable :: avec(:), bvec(:), cvec(:)
    real(dp) :: tl, tp

    tl = 0.0_dp
    tp = 0.0_dp
    do ndim = 100, 12000, 100
      allocate(avec(ndim))

      call random_number(avec)
      avec(:) = avec - 0.5

      bvec = avec
      tl = tl - timing()
      bvec(trueloc(bvec > 0)) = 0.0
      tl = tl + timing()

      cvec = avec
      tp = tp - timing()
      block
        integer :: i
        cvec(pack([(i, i=1, size(cvec))], cvec > 0)) = 0.0
      end block
      tp = tp + timing()

      call check(error, all(bvec == cvec))
      deallocate(avec, bvec, cvec)
      if (allocated(error)) exit
    end do
    call report("trueloc", tl, "pack", tp)
  end subroutine test_trueloc_pack

  subroutine test_falseloc_empty(error)
    !> Error handling
    type(error_type), allocatable, intent(out) :: error

    integer :: ndim
    real, allocatable :: avec(:), bvec(:)

    do ndim = 100, 12000, 100
      allocate(avec(ndim))

      call random_number(avec)

      bvec = avec
      bvec(falseloc(bvec > 0)) = 0.0

      call check(error, all(bvec == avec))
      deallocate(avec, bvec)
      if (allocated(error)) exit
    end do
  end subroutine test_falseloc_empty

  subroutine test_falseloc_all(error)
    !> Error handling
    type(error_type), allocatable, intent(out) :: error

    integer :: ndim
    real, allocatable :: avec(:)

    do ndim = 100, 12000, 100
      allocate(avec(-ndim/2:ndim))

      call random_number(avec)

      avec(falseloc(avec < 0, lbound(avec, 1))) = 0.0

      call check(error, all(avec == 0.0))
      deallocate(avec)
      if (allocated(error)) exit
    end do
  end subroutine test_falseloc_all

  subroutine test_falseloc_where(error)
    !> Error handling
    type(error_type), allocatable, intent(out) :: error

    integer :: ndim
    real, allocatable :: avec(:), bvec(:), cvec(:)
    real(dp) :: tl, tw

    tl = 0.0_dp
    tw = 0.0_dp
    do ndim = 100, 12000, 100
      allocate(avec(ndim))

      call random_number(avec)
      avec(:) = avec - 0.5

      bvec = avec
      tl = tl - timing()
      bvec(falseloc(bvec > 0)) = 0.0
      tl = tl + timing()

      cvec = avec
      tw = tw - timing()
      where(.not.(cvec > 0)) cvec = 0.0
      tw = tw + timing()

      call check(error, all(bvec == cvec))
      deallocate(avec, bvec, cvec)
      if (allocated(error)) exit
    end do
    call report("falseloc", tl, "where", tw)
  end subroutine test_falseloc_where

  subroutine test_falseloc_merge(error)
    !> Error handling
    type(error_type), allocatable, intent(out) :: error

    integer :: ndim
    real, allocatable :: avec(:), bvec(:), cvec(:)
    real(dp) :: tl, tm

    tl = 0.0_dp
    tm = 0.0_dp
    do ndim = 100, 12000, 100
      allocate(avec(ndim))

      call random_number(avec)
      avec(:) = avec - 0.5

      bvec = avec
      tl = tl - timing()
      bvec(falseloc(bvec > 0)) = 0.0
      tl = tl + timing()

      cvec = avec
      tm = tm - timing()
      cvec(:) = merge(cvec, 0.0, cvec > 0)
      tm = tm + timing()

      call check(error, all(bvec == cvec))
      deallocate(avec, bvec, cvec)
      if (allocated(error)) exit
    end do
    call report("falseloc", tl, "merge", tm)
  end subroutine test_falseloc_merge

  subroutine test_falseloc_pack(error)
    !> Error handling
    type(error_type), allocatable, intent(out) :: error

    integer :: ndim
    real, allocatable :: avec(:), bvec(:), cvec(:)
    real(dp) :: tl, tp

    tl = 0.0_dp
    tp = 0.0_dp
    do ndim = 100, 12000, 100
      allocate(avec(ndim))

      call random_number(avec)
      avec(:) = avec - 0.5

      bvec = avec
      tl = tl - timing()
      bvec(falseloc(bvec > 0)) = 0.0
      tl = tl + timing()

      cvec = avec
      tp = tp - timing()
      block
        integer :: i
        cvec(pack([(i, i=1, size(cvec))], cvec < 0)) = 0.0
      end block
      tp = tp + timing()

      call check(error, all(bvec == cvec))
      deallocate(avec, bvec, cvec)
      if (allocated(error)) exit
    end do
    call report("falseloc", tl, "pack", tp)
  end subroutine test_falseloc_pack

  subroutine report(l1, t1, l2, t2)
    character(len=*), intent(in) :: l1, l2
    real(dp), intent(in) :: t1, t2
    character(len=*), parameter :: fmt = "f6.4"

    !$omp critical
    print '(2x, "[Timing]", *(1x, g0))', &
      l1//":", to_string(t1, fmt)//"s", &
      l2//":", to_string(t2, fmt)//"s", &
      "ratio:", to_string(t1/t2, "f4.1")
    !$omp end critical
  end subroutine report

  function timing() result(time)
    real(dp) :: time

    integer(i8) :: time_count, time_rate, time_max
    call system_clock(time_count, time_rate, time_max)
    time = real(time_count, dp)/real(time_rate, dp)
  end function timing

end module test_logicalloc

program tester
  use, intrinsic :: iso_fortran_env, only : error_unit
  use testdrive, only : run_testsuite, new_testsuite, testsuite_type
  use test_logicalloc, only : collect_logicalloc
  implicit none
  integer :: stat, is
  type(testsuite_type), allocatable :: testsuites(:)
  character(len=*), parameter :: fmt = '("#", *(1x, a))'

  stat = 0

  testsuites = [ &
    new_testsuite("logicalloc", collect_logicalloc) &
    ]

  do is = 1, size(testsuites)
    write(error_unit, fmt) "Testing:", testsuites(is)%name
    call run_testsuite(testsuites(is)%collect, error_unit, stat)
  end do

  if (stat > 0) then
    write(error_unit, '(i0, 1x, a)') stat, "test(s) failed!"
    error stop
  end if
end program

max_stars_repo_name: freevryheid/stdlib

max_issues_repo_name: freevryheid/stdlib

max_forks_repo_name: freevryheid/stdlib


text:
      COMPLEX*16 FUNCTION javsubU0(Qs,Us,mtp2,mbt2,mv2)
*------------------------------------------------------
      REAL*8 Qs,Us,zet2,rmv2,eps
      COMPLEX*16 mtp2,mbt2,mv2,ieps,XSPENZ
      PARAMETER (eps=1D-20)
      PARAMETER (zet2=1.6449340668482264364724151666460251D0)
      ieps=DCMPLX(0d0,eps)
      rmv2=DREAL(mv2)
*
      javsubU0=
     &+1d0/Qs*(
     &   +CDLOG(-Qs/(Qs-ieps))*CDLOG((Qs+mv2)/rmv2)
     &   -1d0/2*CDLOG(-mv2/(Qs+mv2))**2
     &   +1d0/2*CDLOG(+mv2/(Qs+mv2))**2
     &   -3*zet2)
     &-1d0/(mtp2+Us)*(
     &   +CDLOG(Qs*mtp2/(mtp2+Us)**2)*CDLOG((Qs+mv2)/mv2)
     &   -1d0/2*CDLOG(mv2/(Qs+mv2))**2
     &   -XSPENZ((Us+mtp2-mv2)/(-mv2))
     &   -XSPENZ(mv2/(Qs+mv2))
     &   +2*zet2     )
*---------------------
      RETURN
      END


text:
      PROGRAM xdfpmin
C     driver for routine dfpmin
      INTEGER NDIM
      REAL GTOL
      PARAMETER(NDIM=2,GTOL=1.0E-4)
      COMMON /stats/ nfunc,ndfunc
      INTEGER iter,ndfunc,nfunc
      REAL fret,p(NDIM)
      EXTERNAL func,dfunc
      write(*,'(/1x,a)') 'True minimum is at (-2.0,+-0.89442719)'
      nfunc=0
      ndfunc=0
      p(1)=.1
      p(2)=4.2
      write(*,'(/1x,a,2(f7.4,a))') 'Starting vector: (',
     *     p(1),',',p(2),')'
      call dfpmin(p,NDIM,GTOL,iter,fret,func,dfunc)
      write(*,'(1x,a,i3)') 'Iterations:',iter
      write(*,'(1x,a,i3)') 'Func. evals:',nfunc
      write(*,'(1x,a,i3)') 'Deriv. evals:',ndfunc
      write(*,'(1x,a,2(f9.6,a))') 'Solution vector: (',
     *     p(1),',',p(2),')'
      write(*,'(1x,a,e14.6)') 'Func. value at solution',fret
      END
      REAL FUNCTION func(x)
      INTEGER ndfunc,nfunc
      COMMON /stats/ nfunc,ndfunc
      REAL x(*)
      nfunc=nfunc+1
      func=10.*(x(2)**2*(3.-x(1))-x(1)**2*(3.+x(1)))**2+(2.+x(1))**2/
     *     (1.+(2.+x(1))**2)
      END
      SUBROUTINE dfunc(x,df)
      INTEGER NMAX
      PARAMETER (NMAX=50)
      INTEGER ndfunc,nfunc
      COMMON /stats/ nfunc,ndfunc
      REAL x(*),df(NMAX)
      ndfunc=ndfunc+1
      df(1)=20.*(x(2)**2*(3.-x(1))-x(1)**2*(3.+x(1)))*(-x(2)**2-6.*
     *     x(1)-3.*x(1)**2)+2.*(2.+x(1))/(1.+(2.+x(1))**2)-
     *     2.*(2.+x(1))**3/(1.+(2.+x(1))**2)**2
      df(2)=40.*(x(2)**2*(3.-x(1))-x(1)**2*(3.+x(1)))*x(2)*(3.-x(1))
      return
      END

max_stars_repo_path: NR-Functions/Numerical Recipes- Example & Functions/Examples/xdfpmin.for


program driver

  ! driver for gfs dynamical core

  use init_mod, only: init
  use run_mod, only: run
  use finalize_mod, only: finalize

  implicit none

  ! read in namelist, initial conditions,
  ! initialize physics and dynamics, allocate arrays.
  call init()
  ! run model, write out at specified intervals.
  call run()
  ! deallocate arrays.
  call finalize()

end program driver


   PROGRAM TEST_CSEM_WaterMW_SfcOptics

      USE CSEM_Type_Kinds, ONLY: fp    => CSEM_fp, &
                                 Short => CSEM_Short
      USE CSEM_WaterMW_SfcOptics ,ONLY: yy_type=>ivar_type, &
                              CSEM_Compute_WaterMW_SfcOptics, &
                              CSEM_Compute_WaterMW_SfcOptics_TL,&
                              CSEM_Compute_WaterMW_SfcOptics_AD
      USE CSEM_Define
      USE CSEM_Model_Manager
      USE FASTEM_Coeff_Reader
      USE CSEM_LifeCycle

      IMPLICIT NONE

      INTEGER, PARAMETER  :: N_Angles = 1
      INTEGER, PARAMETER  :: N_STOKES = 4
      INTEGER, PARAMETER  :: SUCCESS  = 0

      TYPE(CSEM_Water_Surface)    :: Surface
      TYPE(CSEM_SfcOptics_Type)   :: SfcOptics
      TYPE(CSEM_Options_Type)     :: Options
      TYPE(yy_type)               :: iVar

      TYPE(CSEM_Model_ID)         :: CSEM_ALG

      INTEGER(SHORT) :: Error_Status

      !Arguments
      REAL(fp) :: Sensor_Azimuth_Angle
      REAL(fp) :: Transmittance

      ! REAL(fp) :: Angles(N_Angles) =  (/15.0_fp,30.0_fp,45.0_fp,60.0_fp,75.0_fp/)
      REAL(fp) :: Angles(N_Angles) =  (/45.0_fp/)

      INTEGER(SHORT) ::j
      Error_Status =CSEM_INIT("alg.list")
      CSEM_ALG = inq_model_option("MW_WATER")
      !CSEM_ALG%NAME = "RTTOV_FASTEM_V6"
      !CALL set_model_option(CSEM_ALG)

      CALL SfcOptics%Init(N_Angles=N_Angles)

      SfcOptics%Frequency          =   88.997_fp !23.800904_fp
      SfcOptics%Angle              =   Angles
      Surface%Wind_Speed           =   11.0_fp
      Surface%Water_Temperature    =   275.0_fp
      Surface%Salinity             =   33.00_fp
      Surface%Wind_Direction       =   50.0_fp
      Sensor_Azimuth_Angle         =   0.0_fp
      Transmittance                =   0.6_fp

      SfcOptics%Sensor_Azimuth_Angle = Sensor_Azimuth_Angle
      Options%Atmos%Transmittance = Transmittance

      ! forward
      Error_Status = CSEM_Compute_WaterMW_SfcOptics( &
                 Surface                  ,&  ! Input
                 SfcOptics                ,&  ! Input
                 Options, ivar)

      WRITE(*,*)TRIM(CSEM_ALG%NAME)
      WRITE(*,*)'Frequency:', SfcOptics%Frequency
      WRITE(*,*)'Geometry:',  SfcOptics%Angle
      WRITE(*,*)'Emissivity(1):', SfcOptics%Emissivity(:,1)
      WRITE(*,*)'Reflectivity(1):', (SfcOptics%Reflectivity(j,1,j,1), j=1,N_Angles)
      CALL CSEM_Destroy

  END PROGRAM TEST_CSEM_WaterMW_SfcOptics


text:
      subroutine read_lga(file_lga, mmax  , nmax  , nolay , nosegl, &
                          noq1    , noq2  , noq3  , lgrid )

      ! function : read a lga file and check dimensions

      ! global declarations

      use filmod                   ! module contains everything for the files
      implicit none

      ! declaration of the arguments

      type(t_dlwqfile)                       :: file_lga               ! aggregation-file
      integer                                :: mmax                   ! grid cells m direction
      integer                                :: nmax                   ! grid cells n direction
      integer                                :: nolay                  ! nolay
      integer                                :: nosegl                 ! nosegl
      integer                                :: noq1                   ! noq1
      integer                                :: noq2                   ! noq2
      integer                                :: noq3                   ! noq3
      integer                                :: lgrid(nmax,mmax)       ! active grid table

      ! local declarations

      integer                                :: mmaxd                  ! grid cells m direction from lga file
      integer                                :: nmaxd                  ! grid cells n direction from lga file
      integer                                :: ioerr                  ! error on file
      integer                                :: m                      ! loop counter
      integer                                :: n                      ! loop counter
      integer                                :: lunrep                 ! unit number report file

      call getmlu(lunrep)

      call dlwqfile_open(file_lga)
      read(file_lga%unit_nr,iostat=ioerr) nmaxd, mmaxd, nosegl, nolay, noq1, noq2, noq3
      if ( ioerr .ne. 0 ) then
         write(lunrep,*) ' error reading lga file'
         call srstop(1)
      endif

      if ( nmaxd.ne.nmax .or. mmaxd.ne.mmax ) then
         write(lunrep,*) ' dimensions lga file differ from input hydrodynamics'
         write(lunrep,*) ' mmax hyd file:',mmax
         write(lunrep,*) ' nmax hyd file:',nmax
         write(lunrep,*) ' mmax lga file:',mmaxd
         write(lunrep,*) ' nmax lga file:',nmaxd
         call srstop(1)
      endif

      read(file_lga%unit_nr,iostat=ioerr) ((lgrid(n,m),n=1,nmax),m=1,mmax)
      if ( ioerr .ne. 0 ) then
         write(lunrep,*) ' error reading lga file'
         call srstop(1)
      endif

      return
      end

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


module global_module

  implicit none

  ! Time integration

  real(kind=8), parameter :: t_end      = 1.0d-1 ! [s] total simulated time (from 0 to t_end)
  real(kind=8), parameter :: CFL_target = 0.25

  ! Specify domain and discretization. This version includes a flat plate in the domain,
  ! therefore the domain must contain it.

  real(kind=8), parameter :: x_min = -1.75 ! [m] needs to be NEGATIVE
  real(kind=8), parameter :: x_max = 2.75  ! [m] needs to be LARGER THAN PLATE LENGTH L_plate
  real(kind=8), parameter :: y_min = -0.7 ! [m] needs to be NEGATIVE
  real(kind=8), parameter :: y_max = 2.1  ! [m] needs to be POSITIVE

  ! Total number of cells (including 2 ghost cells for each side)
  ! (using 2 ghost cells makes it easy to reach second order accuracy)

  integer, parameter :: Nx = 100
  integer, parameter :: Ny = 100

  real(kind=8), parameter :: dx = (x_max-x_min)/(Nx-4)
  real(kind=8), parameter :: dy = (y_max-y_min)/(Ny-4)

  ! Flat plate inside the domain.
  ! The plate is placed with the leading edge at (x,y) = (0,0).
  ! Make sure the domain is chosen correspondingly.

  real(kind=8), parameter :: L_plate = 1 ! [m] length of the flat plate

  integer, parameter :: j_plate   = floor(Ny/(y_max - y_min)*(0 - y_min)) ! find j-index of plate
  integer, parameter :: i_1_plate = floor(Nx/(x_max - x_min)*(0 - x_min)) ! starting coordinate of plate
  integer, parameter :: i_2_plate = floor(Nx/(x_max - x_min)*(L_plate - x_min)) ! starting coordinate of plate

  ! Free stream

  real(kind=8), parameter :: rho0 = 1.225d0   ! [kg/m3]
  real(kind=8), parameter :: ux0  = 0.0     ! [m/s]
  real(kind=8), parameter :: uy0  = 2867.2  ! [m/s]
  real(kind=8), parameter :: T0   = 300.0   ! [K]

  ! Utilities
  real(kind=8) :: ws_maxabs

end module


program main
	use mpi

	implicit none

	character *(MPI_MAX_PROCESSOR_NAME) processor_name
	integer myid, numprocs, namelen, rc, ierr

	call MPI_INIT(ierr)
	call MPI_COMM_RANK(MPI_COMM_WORLD, myid, ierr)
	call MPI_COMM_SIZE(MPI_COMM_WORLD, numprocs,ierr)
	call MPI_GET_PROCESSOR_NAME(processor_name, namelen, ierr)

	write(*,100) myid,numprocs,processor_name

	call MPI_FINALIZE(ierr)
	stop

100 FORMAT('Hello World! Process,', I2, ' of ', I1, ' on ', 20A)
end program main


text:
      PROGRAM PROFILEI
      IMPLICIT NONE
C
C  hycom_profile2zi - Usage: hycom_profile2zi archv.txt zi.txt archz.txt [itype]
C
C                 converts an HYCOM isopycnal text profile file to
C                 a finite volume z-level text profile file.
C
C   archv.txt is assumed to be an HYCOM archive text profile file
C   zi.txt    is assumed to contain a list of interface depths, one per line
C               the first layer is assumed to be from the surface.
C               the last  layer need not reach the bottom, so set the
C               last interface depth very deep to guarentee full coverage.
C               z-levels are nominally half way between interface depths.
C   archz.txt will be the output, new cell, text profile file
C   itype     is the input interpolation type (default 1)
C                =0; piecewise constant  method (PCM) or donor cell
C                =1; piecewise linear    method (PLM) or VanLeer
C                =2; piecewise parabolic method (PPM)
C
C  all input interpolation schemes conserve the original cell average,
C  and the output is the average of the interpolation profile across
C  each new cell.  So the total depth average is conserved (provided
C  the last interface is at or below the total depth).
C
C  If some or all of the target coordinate system is in density space,
C  use hycom_profile_remap to both choose the interface depths and
C  perform the averaging over the new cells.
C
C  this version for "serial" Unix systems.
C
C  Tim Campbell,       Mississippi State University
C  Alan J. Wallcraft,  Naval Research Laboratory,  Sep. 2002 (Aug. 2005).
C
      INTEGER      IARGC
      INTEGER      NARG
      CHARACTER*240 CARG
C
      CHARACTER*240 CFILEA,CFILEB,CFILEC
      CHARACTER*240 CLINE
      REAL          THK,DUM5(5),FLAG
      REAL          PO(9999)
      REAL          SIT(5),SOT(5)
      INTEGER       IOS,ITYPE,K,KK,KDM,KO
      INTEGER       I,KMAX
C
      REAL, ALLOCATABLE :: SI(:,:),PI(:),SO(:,:),ZO(:)
C
      INTEGER       NSAMP
      REAL          PMAX
      PARAMETER(NSAMP=1,PMAX=500.0)
C
C     READ ARGUMENTS.
C
      NARG = IARGC()
C
      IF     (NARG.EQ.4) THEN
        CALL GETARG(1,CFILEA)
        CALL GETARG(2,CFILEB)
        CALL GETARG(3,CFILEC)
        CALL GETARG(4,CLINE)
        READ(CLINE,*) ITYPE
      ELSEIF (NARG.EQ.3) THEN
        CALL GETARG(1,CFILEA)
        CALL GETARG(2,CFILEB)
        CALL GETARG(3,CFILEC)
        ITYPE = 1
      ELSE
        WRITE(6,*)
     +    'Usage: hycom_profile2zi archv.txt zi.txt archz.txt [itype]'
        CALL EXIT(1)
      ENDIF
C
C     OPEN ALL FILES.
C
      OPEN(UNIT=11, FILE=CFILEA, FORM='FORMATTED', STATUS='OLD',
     +     IOSTAT=IOS)
      IF     (IOS.NE.0) THEN
        WRITE(6,*) 'Error: can''t open ',TRIM(CFILEA)
        WRITE(6,*) 'ios   = ',ios
        CALL EXIT(3)
      ENDIF
      OPEN(UNIT=12, FILE=CFILEB, FORM='FORMATTED', STATUS='OLD',
     +     IOSTAT=IOS)
      IF     (IOS.NE.0) THEN
        WRITE(6,*) 'Error: can''t open ',TRIM(CFILEB)
        WRITE(6,*) 'ios   = ',ios
        CALL EXIT(4)
      ENDIF
      OPEN(UNIT=21, FILE=CFILEC, FORM='FORMATTED', STATUS='NEW',
     +     IOSTAT=IOS)
      IF     (IOS.NE.0) THEN
        WRITE(6,*) 'Error: can''t open ',TRIM(CFILEC)
        WRITE(6,*) 'ios   = ',ios
        CALL EXIT(5)
      ENDIF
C
C     COPY PROFILE HEADER TO OUTPUT.
C
      DO K= 1,4
        READ( 11,'(a)')      CLINE
        WRITE(21,'(a)') TRIM(CLINE)
      ENDDO
      READ( 11,'(a)') CLINE
C
C     READ THE ISOPYCNAL PROFILE, TO GET KDM.
C
      KDM   = -1
      DO K= 1,99
        READ(11,'(a)',IOSTAT=IOS) CLINE
        IF     (IOS.NE.0) THEN
          IF     (K.NE.KDM+1) THEN
            WRITE(6,*) 'Error: inconsistent input profile'
            CALL EXIT(6)
          ENDIF
          EXIT
        ENDIF
        READ(CLINE,*) KDM
      ENDDO
C
C     RE-READ THE ISOPYCNAL PROFILE.
C
      ALLOCATE( PI(KDM+1), SI(KDM,5) )
C
      REWIND(11)
      DO K= 1,5
        READ(11,*)
      ENDDO
      PI(1) =  0.0
      DO K= 1,KDM
        READ(11,'(a)',IOSTAT=IOS) CLINE
        IF     (IOS.NE.0) THEN
          WRITE(6,*) 'Error: inconsistent input profile'
          CALL EXIT(6)
        ENDIF
        READ(CLINE,*) KDM,SI(K,1),SI(K,2),SI(K,3),SI(K,4),SI(K,5),THK
        PI(K+1) = PI(K) + THK
        IF     (THK.EQ.0.0) THEN
          DO KK= 1,5
            SI(K,KK)=SI(K-1,KK)
          ENDDO !kk
        ENDIF
      ENDDO
      CLOSE(11)
C
C     STORE COPY OF INPUT AVERAGES AND
C     COMPUTE INPUT TOTAL COLUMN AVERAGES
C
      KMAX=0
      DO KK= 1,5
        SIT(KK)=0.0
      ENDDO !kk
      DO K=1,KDM
        IF (PI(K).LE.PMAX) KMAX=KMAX+1
        THK=PI(K+1)-PI(K)
        DO KK= 1,5
          SIT(KK)=SIT(KK)+THK*SI(K,KK)
        ENDDO !kk
      ENDDO
      DO KK= 1,5
        SIT(KK)=SIT(KK)/PI(KDM+1)
      ENDDO !kk
C
C     READ THE ZI FILE.
C
      PO(1) = 0.0
      DO K= 2,99999
        READ(12,'(a)',IOSTAT=IOS) CLINE
        IF     (IOS.NE.0) THEN
          KO = K-2
          EXIT
        ENDIF
        READ(CLINE,*) PO(K)
        IF     (K.EQ.2 .AND. PO(K).LE.0.0) THEN
          WRITE(6,*) 'Error: inconsistent output profile'
          WRITE(6,*) '1st value must be > 0.0'
          WRITE(6,*) PO(2)
          CALL EXIT(7)
        ELSEIF (PO(K).LT.PO(K-1)) THEN
          WRITE(6,*) 'Error: inconsistent output profile'
          WRITE(6,*) k-1,'-th value < previous value'
          WRITE(6,*) PO(2:K)
          CALL EXIT(8)
        ENDIF
        PO(K) =MIN( PO(K), PI(KDM+1) )
      ENDDO
      CLOSE(12)
C
      ALLOCATE( ZO(KO), SO(KO,5) )
C
      DO K= 1,KO
        ZO(K) = 0.5*(PO(K)+PO(K+1))
      ENDDO
C
C     FORM NEW CELL AVERAGES
C
      FLAG = 9999.99
      IF     (ITYPE.EQ.0) THEN
        CALL LAYER2ZI_PCM(SI,PI,KDM,5,
     &                    SO,PO,KO,   FLAG)
      ELSEIF (ITYPE.EQ.1) THEN
        CALL LAYER2ZI_PLM(SI,PI,KDM,5,
     &                    SO,PO,KO,   FLAG)
      ELSEIF (ITYPE.EQ.2) THEN
        CALL LAYER2ZI_PPM(SI,PI,KDM,5,
     &                    SO,PO,KO,   FLAG)
      ELSE
        WRITE(6,*)
     +    'Usage: hycom_profile2zi archv.txt zi.txt archz.txt [itype]'
        WRITE(6,*) 'unsupported itype value'
        CALL EXIT(1)
      ENDIF
C
C     COMPUTE AND PRINT NEW TOTAL COLUMN AVERAGE AND DEVIATION METRIC
C
      DO KK= 1,5
        SOT(KK)=0.0
        DO K=1,KO
          THK=PO(K+1)-PO(K)
          SOT(KK)=SOT(KK)+THK*SO(K,KK)
        ENDDO !k
        SOT(KK)=SOT(KK)/PO(KO+1)
      ENDDO !kk
      WRITE(*,'(/A)')       '        LABEL:  U  V  T  S  R'
      WRITE(*,'(A,5e14.6)') ' INPUT TOTALS: ',SIT(:)
      WRITE(*,'(A,5e14.6)') 'OUTPUT TOTALS: ',SOT(:)
C
C     OUTPUT.
C
      WRITE(21,'(a,a)')
     &  '#   k',
     &  '    utot    vtot    temp    saln    dens    thkns      dpth'
      DO K= 1,KO
        WRITE(21,'(i4,1x,2f8.2,3f8.3,f9.3,f10.3)')
     &    K,
     &    SO(K,1),             !cm/s
     &    SO(K,2),             !cm/s
     &    SO(K,3),             !degC
     &    SO(K,4),             !psu
     &    SO(K,5),             !SigmaT
     &    PO(K+1)-PO(K),       !m
     &    ZO(K)                !m
      ENDDO
      END

      subroutine layer2zi_pcm(si,pi,ki,ks,
     &                        so,po,ko,   flag)
      implicit none
c
      integer ki,ks,ko
      real    si(ki,ks),pi(ki+1),
     &        so(ko,ks),po(ko+1),flag
c
c**********
c*
c  1) remap from one set of vertical cells to another.
c     method: piecewise constant across each input cell
c             the output is the average of the interpolation
c             profile across each output cell.
c
c  2) input arguments:
c       si    - scalar fields in pi-layer space
c       pi    - layer interface depths (non-negative m)
c                 pi(   1) is the surface
c                 pi(ki+1) is the bathymetry
c       ki    - 1st dimension of si     (number of  input layers)
c       ks    - 2nd dimension of si,so  (number of fields)
c       po    - target interface depths (non-negative m)
c                 po(k+1) >= po(k)
c       ko    - 1st dimension of so     (number of output layers)
c       flag  - data void (land) marker
c
c  3) output arguments:
c       so    - scalar fields in po-layer space
c
c  4) except at data voids, must have:
c           pi(   1) == zero (surface)
c           pi( l+1) >= pi(l)
c           pi(ki+1) == bathymetry
c           0 <= po(k) <= po(k+1)
c      output layers completely below the bathymetry inherit values
c      from the layer above.
c
c  5) Alan J. Wallcraft,  Naval Research Laboratory,  Sep. 2002 (Aug. 2005).
c*
c**********
c
      real       thin
      parameter (thin=1.e-6)  ! minimum layer thickness (no division by 0.0)
c
      integer i,k,l,lf
      real    q,zb,zt,sok(ks)
c
      if     (si(1,1).eq.flag) then
        do k= 1,ko
          do i= 1,ks
            so(k,i) = flag  !land
          enddo !i
        enddo !k
      else
        lf=1
        zb=po(1)
        do k= 1,ko
          zt = zb
          zb = po(k+1)
*         WRITE(6,*) 'k,zt,zb = ',k,zt,zb
          if     (zb-zt.lt.thin .or. zt.ge.pi(ki+1)) then
c
c ---       thin or bottomed layer, values taken from layer above
c
            do i= 1,ks
              so(k,i) = so(k-1,i)
            enddo !i
          else
c
c           form layer averages.
c
            if     (pi(lf).gt.zt) then
              WRITE(6,*) 'bad lf = ',lf
              stop
            endif
            do i= 1,ks
              sok(i) = 0.0
            enddo !i
            do l= lf,ki
              if     (pi(l).gt.zb) then
*               WRITE(6,*) 'l,lf= ',l,lf,l-1
                lf = l-1
                exit
              elseif (pi(l).ge.zt .and. pi(l+1).le.zb) then
c
c               the input layer is completely inside the output layer
c
                q   = max(pi(l+1)-pi(l),thin)/(zb-zt)
                do i= 1,ks
                  sok(i) = sok(i) + q*si(l,i)
                enddo !i
*               WRITE(6,*) 'L,q = ',l,q
              else
c
c               the input layer is partially inside the output layer
c
                q   = max(min(pi(l+1),zb)-max(pi(l),zt),thin)/(zb-zt)
                do i= 1,ks
                  sok(i) = sok(i) + q*si(l,i)
                enddo !i
*               WRITE(6,*) 'l,q = ',l,q
              endif
            enddo !l
            do i= 1,ks
              so(k,i) = sok(i)
            enddo !i
          endif
        enddo !k
      endif
      return
      end subroutine layer2zi_pcm

      subroutine layer2zi_plm(si,pi,ki,ks,
     &                        so,po,ko,   flag)
      implicit none
c
      integer ki,ks,ko
      real    si(ki,ks),pi(ki+1),
     &        so(ko,ks),po(ko+1),flag
c
c**********
c*
c  1) remap from one set of vertical cells to another.
c     method: piecewise linear across each input cell
c             the output is the average of the interpolation
c             profile across each output cell.
c
c  2) input arguments:
c       si    - scalar fields in pi-layer space
c       pi    - layer interface depths (non-negative m)
c                 pi(   1) is the surface
c                 pi(ki+1) is the bathymetry
c       ki    - 1st dimension of si     (number of  input layers)
c       ks    - 2nd dimension of si,so  (number of fields)
c       po    - target interface depths (non-negative m)
c                 po(k+1) >= po(k)
c       ko    - 1st dimension of so     (number of output layers)
c       flag  - data void (land) marker
c
c  3) output arguments:
c       so    - scalar fields in po-layer space
c
c  4) except at data voids, must have:
c           pi(   1) == zero (surface)
c           pi( l+1) >= pi(l)
c           pi(ki+1) == bathymetry
c           0 <= po(k) <= po(k+1)
c      output layers completely below the bathymetry inherit values
c      from the layer above.
c
c  5) Tim Campbell, Mississippi State University, October 2002.
C     Alan J. Wallcraft,  Naval Research Laboratory,  Aug. 2005.
c*
c**********
c
      real,parameter :: thin=1.e-6  !minimum layer thickness
c
      integer i,k,l,lf
      real    q,qc,zb,zc,zt,sok(ks)
      real    sis(ki,ks),pit(ki+1)
c
      if     (si(1,1).eq.flag) then
        do k= 1,ko
          do i= 1,ks
            so(k,i) = flag  !land
          enddo !i
        enddo
      else
c ---   compute PLM slopes for input layers
        do k=1,ki
          pit(k)=max(pi(k+1)-pi(k),thin)
        enddo
        call plm(pit,si,sis,ki,ks)
c ---   compute output layer averages
        lf=1
        zb=po(1)
        do k= 1,ko
          zt = zb
          zb = po(k+1)
*         WRITE(6,*) 'k,zt,zb = ',k,zt,zb
          if     (zb-zt.lt.thin .or. zt.ge.pi(ki+1)) then
c
c ---       thin or bottomed layer, values taken from layer above
c
            do i= 1,ks
              so(k,i) = so(k-1,i)
            enddo !i
          else
c
c           form layer averages.
c
            if     (pi(lf).gt.zt) then
              WRITE(6,*) 'bad lf = ',lf
              stop
            endif
            do i= 1,ks
              sok(i) = 0.0
            enddo !i
            do l= lf,ki
              if     (pi(l).gt.zb) then
*               WRITE(6,*) 'l,lf= ',l,lf,l-1
                lf = l-1
                exit
              elseif (pi(l).ge.zt .and. pi(l+1).le.zb) then
c
c               the input layer is completely inside the output layer
c
                q   = max(pi(l+1)-pi(l),thin)/(zb-zt)
                do i= 1,ks
                  sok(i) = sok(i) + q*si(l,i)
                enddo !i
*               WRITE(6,*) 'L,q = ',l,q
              else
c
c               the input layer is partially inside the output layer
c               average of linear profile is its center value
c
                q   = max( min(pi(l+1),zb)-max(pi(l),zt), thin )/(zb-zt)
                zc  = 0.5*(min(pi(l+1),zb)+max(pi(l),zt))
                qc  = (zc-pi(l))/pit(l) - 0.5
                do i= 1,ks
                  sok(i) = sok(i) + q*(si(l,i) + qc*sis(l,i))
                enddo !i
*               WRITE(6,*) 'l,q,qc = ',l,q,qc
              endif
            enddo !l
            do i= 1,ks
              so(k,i) = sok(i)
            enddo !i
          endif
        enddo !k
      endif
      return
      end subroutine layer2zi_plm

      subroutine plm(pt, s,ss,ki,ks)
      implicit none
c
      integer ki,ks
      real    pt(ki+1),s(ki,ks),ss(ki,ks)
c
c**********
c*
c  1) generate a monotonic PLM interpolation of a layered field
c
c  2) input arguments:
c       pt    - layer interface thicknesses (non-zero)
c       s     - scalar fields in layer space
c       ki    - 1st dimension of s (number of layers)
c       ks    - 2nd dimension of s (number of fields)
c
c  3) output arguments:
c       ss    - scalar field slopes for PLM interpolation
c
c  4) except at data voids, must have:
c           pi(   1) == zero (surface)
c           pi( l+1) >= pi(:,:,l)
c           pi(ki+1) == bathymetry
c
c  5) Tim Campbell, Mississippi State University, September 2002.
c*
c**********
c
      integer l
      real    ql(ki),qc(ki),qr(ki)
c
      !compute grid spacing ratios for slope computations
      ql(1)=0.0
      qc(1)=0.0
      qr(1)=0.0
      do l=2,ki-1
        ql(l)=2.0*pt(l)/(pt(l-1)+pt(l))
        qc(l)=2.0*pt(l)/(pt(l-1)+2.0*pt(l)+pt(l+1))
        qr(l)=2.0*pt(l)/(pt(l)+pt(l+1))
      enddo
      ql(ki)=0.0
      qc(ki)=0.0
      qr(ki)=0.0
      !compute normalized layer slopes
      do l=1,ks
        call slope(ql,qc,qr,s(1,l),ss(1,l),ki)
      enddo
      return
      end subroutine plm

      subroutine slope(rl,rc,rr,a,s,n)
      implicit none
c
      integer,intent(in)  :: n
      real,   intent(in)  :: rl(n),rc(n),rr(n),a(n)
      real,   intent(out) :: s(n)
c
c**********
c*
c  1) generate slopes for monotonic piecewise linear distribution
c
c  2) input arguments:
c       rl   - left grid spacing ratio
c       rc   - center grid spacing ratio
c       rr   - right grid spacing ratio
c       a    - scalar field zone averages
c       n    - number of zones
c
c  3) output arguments:
c       s    - zone slopes
c
c  4) Tim Campbell, Mississippi State University, September 2002.
c*
c**********
c
      integer,parameter :: ic=2, im=1, imax=100
      real,parameter :: fracmin=1e-6, dfac=0.5
c
      integer i,j
      real    sl,sc,sr
      real    dnp,dnn,dl,dr,ds,frac
c
c Compute zone slopes
c Campbell Eq(15) -- nonuniform grid
c
      s(1)=0.0
      do j=2,n-1
        sl=rl(j)*(a(j)-a(j-1))
        sr=rr(j)*(a(j+1)-a(j))
        if (sl*sr.gt.0.) then
          s(j)=sign(min(abs(sl),abs(sr)),sl)
        else
          s(j)=0.0
        endif
      enddo
      s(n)=0.0
c
c Minimize discontinuities between zones
c Apply single pass discontinuity minimization: Campbell Eq(19)
c
      do j=2,n-1
        if(s(j).ne.0.0) then
          dl=-0.5*(s(j)+s(j-1))+a(j)-a(j-1)
          dr=-0.5*(s(j+1)+s(j))+a(j+1)-a(j)
          ds=sign(min(abs(dl),abs(dr)),dl)
          s(j)=s(j)+2.0*ds
        endif
      enddo
      return
      end subroutine slope

      subroutine layer2zi_ppm(si,pi,ki,ks,
     &                        so,po,ko,   flag)
      implicit none
c
      integer ki,ks,ko
      real    si(ki,ks),pi(ki+1),
     &        so(ko,ks),po(ko+1),flag
c
c**********
c*
c  1) remap from one set of vertical cells to another.
c     method: piecewise parabolic method across each input cell
c             the output is the average of the interpolation
c             profile across each output cell.
c
c  2) input arguments:
c       si    - scalar fields in pi-layer space
c       pi    - layer interface depths (non-negative m)
c                 pi(   1) is the surface
c                 pi(ki+1) is the bathymetry
c       ki    - 1st dimension of si     (number of  input layers)
c       ks    - 2nd dimension of si,so  (number of fields)
c       po    - target interface depths (non-negative m)
c                 po(k+1) >= po(k)
c       ko    - 1st dimension of so     (number of output layers)
c       flag  - data void (land) marker
c
c  3) output arguments:
c       so    - scalar fields in po-layer space
c
c  4) except at data voids, must have:
c           pi(   1) == zero (surface)
c           pi( l+1) >= pi(l)
c           pi(ki+1) == bathymetry
c           0 <= po(k) <= po(k+1)
c      output layers completely below the bathymetry inherit values
c      from the layer above.
c
c  5) Tim Campbell, Mississippi State University, October 2002.
C     Alan J. Wallcraft,  Naval Research Laboratory,  Aug. 2005.
c*
c**********
c
      real,parameter :: thin=1.e-6  !minimum layer thickness
c
      integer i,k,l,lb,lt
      real    sz,xb,xt,zb,zt
      real    sic(ki,ks,3),pit(ki+1)
c
      if     (si(1,1).eq.flag) then
        do k= 1,ko
          do i= 1,ks
            so(k,i) = flag  !land
          enddo !i
        enddo
      else
c ---   compute PPM coefficients for input layers
        do k=1,ki
          pit(k)=max(pi(k+1)-pi(k),thin)
        enddo
        call ppm(pit,si,sic,ki,ks)
c ---   compute output layer averages
        lb=1
        zb=po(1)
        do k= 1,ko
          zt = zb
          zb = po(k+1)
*         WRITE(6,*) 'k,zt,zb = ',k,zt,zb
          if     (zb-zt.lt.thin .or. zt.ge.pi(ki+1)) then
c
c ---       thin or bottomed layer, values taken from layer above
c
            do i= 1,ks
              so(k,i) = so(k-1,i)
            enddo !i
          else
c
c           form layer averages.
c
            if     (pi(lb).gt.zt) then
              WRITE(6,*) 'bad lb = ',lb
              stop
            endif
            lt=lb !top will always correspond to bottom of previous
            lb=lt !find input layer containing bottom output interface
            do while (pi(lb+1).lt.zb.and.lb.lt.ki)
              lb=lb+1
            enddo
            xt=(zt-pi(lt))/pit(lt)
            xb=(zb-pi(lb))/pit(lb)
            do i= 1,ks
              if     (lt.ne.lb) then
                sz=   pit(lt)*(      sic(lt,i,1)  *(1.-xt)
     &                         +0.5*(sic(lt,i,2)+
     &                               sic(lt,i,3) )*(1.-xt**2)
     &                              -sic(lt,i,3)  *(1.-xt**3)/3.0)
                do l=lt+1,lb-1
                  sz=sz+pit(l)*si(l,i)
                enddo !l
                sz=sz+pit(lb)*(      sic(lb,i,1)  *    xb
     &                         +0.5*(sic(lb,i,2)+
     &                               sic(lb,i,3) )*    xb**2
     &                              -sic(lb,i,3)  *    xb**3 /3.0)
              else
                sz=pit(lt)*(      sic(lt,i,1)  *(xb-xt)
     &                      +0.5*(sic(lt,i,2)+
     &                            sic(lt,i,3) )*(xb**2-xt**2)
     &                           -sic(lt,i,3)  *(xb**3-xt**3)/3.0)
              endif
              so(k,i) = sz/(zb-zt)
            enddo !i
          endif !thin:std layer
        enddo !k
      endif
      return
      end subroutine layer2zi_ppm

      subroutine ppm(pt, s,sc,ki,ks)
      implicit none
c
      integer ki,ks
      real    pt(ki+1),s(ki,ks),sc(ki,ks,3)
c
c**********
c*
c  1) generate a monotonic PPM interpolation of a layered field:
c     Colella, P. & P.R. Woodward, 1984, J. Comp. Phys., 54, 174-201.
c
c  2) input arguments:
c       pt    - layer interface thicknesses (non-zero)
c       s     - scalar fields in layer space
c       ki    - 1st dimension of s (number of layers)
c       ks    - 2nd dimension of s (number of fields)
c
c  3) output arguments:
c       sc    - scalar field coefficients for PPM interpolation
c
c  4) except at data voids, must have:
c           pi(   1) == zero (surface)
c           pi( l+1) >= pi(:,:,l)
c           pi(ki+1) == bathymetry
c
c  5) Tim Campbell, Mississippi State University, September 2002;
C     Alan J. Wallcraft,  Naval Research Laboratory,  August 2007.
c*
c**********
c
      integer j,k,l
      real    da,a6,slj,scj,srj
      real    as(ki),al(ki),ar(ki)
      real     ptjp(ki), pt2jp(ki), ptj2p(ki),
     &        qptjp(ki),qpt2jp(ki),qptj2p(ki),ptq3(ki),qpt4(ki)
c
      !compute grid metrics
      do j=1,ki
         ptq3( j) = pt(j)/(pt(j-1)+pt(j)+pt(j+1))
         ptjp( j) = pt(j)   + pt(j+1)
         pt2jp(j) = pt(j)   + ptjp(j)
         ptj2p(j) = ptjp(j) + pt(j+1)
        qptjp( j) = 1.0/ptjp( j)
        qpt2jp(j) = 1.0/pt2jp(j)
        qptj2p(j) = 1.0/ptj2p(j)
      enddo !j
         ptq3(2) = pt(2)/(pt(1)+ptjp(2))
      do j=3,ki-1
         ptq3(j) = pt(j)/(pt(j-1)+ptjp(j))  !pt(j)/      (pt(j-1)+pt(j)+pt(j+1))
         qpt4(j) = 1.0/(ptjp(j-2)+ptjp(j))  !1.0/(pt(j-2)+pt(j-1)+pt(j)+pt(j+1))
      enddo !j
c
      do l= 1,ks
        !Compute average slopes: Colella, Eq. (1.8)
        as(1)=0.
        do j=2,ki-1
          slj=s(j,  l)-s(j-1,l)
          srj=s(j+1,l)-s(j,  l)
          if (slj*srj.gt.0.) then
            scj=ptq3(j)*( pt2jp(j-1)*srj*qptjp(j)
     &                   +ptj2p(j)  *slj*qptjp(j-1) )
            as(j)=sign(min(abs(2.0*slj),abs(scj),abs(2.0*srj)),scj)
          else
            as(j)=0.
          endif
        enddo !j
        as(ki)=0.
        !Compute "first guess" edge values: Colella, Eq. (1.6)
        al(1)=s(1,l)  !1st layer PCM
        ar(1)=s(1,l)  !1st layer PCM
        al(2)=s(1,l)  !1st layer PCM
        do j=3,ki-1
          al(j)=s(j-1,l)+pt(j-1)*(s(j,l)-s(j-1,l))*qptjp(j-1)
     &         +qpt4(j)*(
     &            2.*pt(j)*pt(j-1)*qptjp(j-1)*(s(j,l)-s(j-1,l))*
     &            ( ptjp(j-2)*qpt2jp(j-1)
     &             -ptjp(j)  *qptj2p(j-1) )
     &            -pt(j-1)*as(j)  *ptjp(j-2)*qpt2jp(j-1)
     &            +pt(j)  *as(j-1)*ptjp(j)  *qptj2p(j-1)
     &              )
          ar(j-1)=al(j)
        enddo !j
        ar(ki-1)=s(ki,l)  !last layer PCM
        al(ki)  =s(ki,l)  !last layer PCM
        ar(ki)  =s(ki,l)  !last layer PCM
        !Impose monotonicity: Colella, Eq. (1.10)
        do j=2,ki-1
          if ((s(j+1,l)-s(j,l))*(s(j,l)-s(j-1,l)).le.0.) then !local extremum
            al(j)=s(j,l)
            ar(j)=s(j,l)
          else
            da=ar(j)-al(j)
            a6=6.0*s(j,l)-3.0*(al(j)+ar(j))
            if     (da*a6 .gt.  da*da) then !peak in right half of zone
              al(j)=3.0*s(j,l)-2.0*ar(j)
            elseif (da*a6 .lt. -da*da) then !peak in left half of zone
              ar(j)=3.0*s(j,l)-2.0*al(j)
            endif
          endif
        enddo !j
        !Set coefficients
        do j=1,ki
          sc(j,l,1)=al(j)
          sc(j,l,2)=ar(j)-al(j)
          sc(j,l,3)=6.0*s(j,l)-3.0*(al(j)+ar(j))
        enddo !j
      enddo !l
      return
      end subroutine ppm

max_stars_repo_name: TillRasmussen/HYCOM-tools

max_issues_repo_name: TillRasmussen/HYCOM-tools

max_forks_repo_name: TillRasmussen/HYCOM-tools


text:
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C     *                                                               *
C     *                  copyright (c) 2011 by UCAR                   *
C     *                                                               *
C     *       University Corporation for Atmospheric Research         *
C     *                                                               *
C     *                      all rights reserved                      *
C     *                                                               *
C     *                     FFTPACK  version 5.1                      *
C     *                                                               *
C     *                 A Fortran Package of Fast Fourier             *
C     *                                                               *
C     *                Subroutines and Example Programs               *
C     *                                                               *
C     *                             by                                *
C     *                                                               *
C     *               Paul Swarztrauber and Dick Valent               *
C     *                                                               *
C     *                             of                                *
C     *                                                               *
C     *         the National Center for Atmospheric Research          *
C     *                                                               *
C     *                Boulder, Colorado  (80307)  U.S.A.             *
C     *                                                               *
C     *                   which is sponsored by                       *
C     *                                                               *
C     *              the National Science Foundation                  *
C     *                                                               *
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C
      SUBROUTINE TABLES (IDO,IP,WA)
      REAL  WA(IDO,IP-1,2)
C
      TPI = 8.*ATAN(1.)
      ARGZ = TPI/REAL(IP)
      ARG1 = TPI/REAL(IDO*IP)
      DO 110 J=2,IP
         ARG2 = REAL(J-1)*ARG1
         DO 100 I=1,IDO
            ARG3 = REAL(I-1)*ARG2
            WA(I,J-1,1) = COS(ARG3)
            WA(I,J-1,2) = SIN(ARG3)
  100    CONTINUE
         IF (IP .LE. 5) GO TO 110
         ARG4 = REAL(J-1)*ARGZ
         WA(1,J-1,1) = COS(ARG4)
         WA(1,J-1,2) = SIN(ARG4)
  110 CONTINUE
      RETURN
      END

max_stars_repo_name: thkrz/fstpack

max_issues_repo_name: thkrz/fstpack

max_forks_repo_name: thkrz/fstpack


SUBROUTINE jima(g)

!-------------------------------------------------------------
!
!  Directs the solutiong with the New Solution algorithm
!   ahot-n-ns. Given an energy group, it sweeps over each
!   angle and then cells constructing A and B matrices to
!   make the Gamma matrix. Then forms the iteration Jacobian
!   matrix. Finally divides by scattering ratio to form the
!   matrix used to solve for the fully converged scalar flux.
!
!-------------------------------------------------------------

USE invar
USE solvar
USE timevar
IMPLICIT NONE
INTEGER, INTENT(IN) :: g
INTEGER :: i, j, k, l, m, n, ieq, quad, jeq, neq, ord, ndx
INTEGER :: xs, xe, incx, ys, ye, incy, sgm, sge, kk, ii, jj, info
REAL*8 :: c, x, y, mu, eta, sig, tempx, tempy, symm, symm2
REAL*8, DIMENSION(order,ordsq,nx,ny) :: xold, yold
REAL*8, DIMENSION((ordsq*nx*ny)) :: src, q, sol

! Size amat, bmat, gmat
neq = ordsq + 2*order
ALLOCATE(amat(neq,neq), bmat(neq,neq), gmat(neq,neq))

! Allocate X-xmat and Y-ymat
ALLOCATE(xmat(nx,order,ordsq,nx,ny), ymat(order,ordsq,nx,ny))

! Allocate all the gamma submatrices sizes
ALLOCATE(gaa(ordsq,ordsq), gax(ordsq,order), gay(ordsq,order))
ALLOCATE(gxa(order,ordsq), gxx(order,order), gxy(order,order))
ALLOCATE(gya(order,ordsq), gyx(order,order), gyy(order,order))

! Allocate the Iteration Jacobian
neq = ordsq*nx*ny
ALLOCATE(jmat(neq,neq))
jmat = 0.0

! Start the loop in all direction
DO n = 1, apo
   mu  = ang(n,1)
   eta = ang(n,2)

   ! Start the loops in the quadrants
   DO quad = 1, 4
      IF (quad == 1) THEN
         xs = 1
         xe = nx
         incx = 1
         ys = 1
         ye = ny
         incy = 1
      ELSE IF (quad == 2) THEN
         xs = nx
         xe = 1
         incx = -1
         ys = 1
         ye = ny
         incy = 1
      ELSE IF (quad == 3) THEN
         xs = nx
         xe = 1
         incx = -1
         ys = ny
         ye = 1
         incy = -1
      ELSE IF (quad == 4) THEN
         xs = 1
         xe = nx
         incx = 1
         ys = ny
         ye = 1
         incy = -1
      END IF

      ! Reset the X matrix at the beginning of a new sweep
      xmat = 0.0

   ! Start the loops
   DO j = ys, ye, incy
      sge = incy
      y = dy(j)

      ! Reset the Y matrix at beginning of new row
      ymat = 0.0

      ! Start the row sweeps
      DO i = xs, xe, incx
         sgm = incx
         ! Get the spatial weights
         x = dx(i)
         m = mat(i,j)
         sig = sigt(m,g)
         ord = lambda
         c = sigs(m,g,g)/sig     ! Scattering ratio
         ! Call for the calculation of the spatial weights
         CALL weitn(ord,x,y,sig,mu,eta)

         ! Construct amat and bmat, A and B
         CALL conab(sgm,sge)
         ! Construct the gamma matrix
         CALL cong
         ! Unload the gamma matrix into useful sub-matrices
         CALL gammas

         ! Save the X and Y matrices for updates
         IF (i /= xs) THEN
            yold = ymat
         END IF
         IF (j /= ys) THEN
            DO jj = ys, ye, incy
               DO ii = xs, xe, incx
                  DO kk = 1, ordsq
                     DO k = 1, order
                        xold(k,kk,ii,jj) = xmat(i,k,kk,ii,jj)
                     END DO
                  END DO
               END DO
            END DO
         END IF

         ! Update X and Y matrices with the old Y matrix
         IF (i /= xs) THEN
            DO jj = ys, j, incy
               DO ii = xs, i, incx
                  DO k = 1, order
                     DO kk = 1, ordsq
                        tempy = 0.0
                        tempx = 0.0
                        DO l = 1, order
                           IF (i /= xe) THEN
                              tempy = tempy + gyy(k,l)*yold(l,kk,ii,jj)
                           END IF
                           IF (j /= ye) THEN
                              tempx = tempx + gxy(k,l)*yold(l,kk,ii,jj)
                           END IF
                        END DO
                        ymat(k,kk,ii,jj) = tempy
                        xmat(i,k,kk,ii,jj) = tempx
                     END DO
                  END DO
               END DO
            END DO
         END IF

         ! Update X and Y matrices with the old X matrix
         IF (j /= ys) THEN
            DO jj = ys, j, incy
               DO ii = xs, i, incx
                  DO k = 1, order
                     DO kk = 1, ordsq
                        tempx = 0.0
                        tempy = 0.0
                        DO l = 1, order
                           IF (i /= xe) THEN
                              tempy = tempy + gyx(k,l)*xold(l,kk,ii,jj)
                           END IF
                           IF (j /= ye) THEN
                              tempx = tempx + gxx(k,l)*xold(l,kk,ii,jj)
                           END IF
                        END DO
                        IF (i == xs) THEN
                           IF (i /= xe) THEN
                              ymat(k,kk,ii,jj) = tempy
                           END IF
                           IF (j /= ye) THEN
                              xmat(i,k,kk,ii,jj) = tempx
                           END IF
                        ELSE
                           IF (i /= xe) THEN
                              ymat(k,kk,ii,jj) = ymat(k,kk,ii,jj) + tempy
                           END IF
                           IF (j /= ye) THEN
                              xmat(i,k,kk,ii,jj) = xmat(i,k,kk,ii,jj) + tempx
                           END IF
                        END IF
                     END DO
                  END DO
               END DO
            END DO
         END IF

         ! Append X and Y matrices
         DO kk = 1, ordsq
            DO k = 1, order
               IF (j /= ye) THEN
                  xmat(i,k,kk,i,j) = c*gxa(k,kk)
               END IF
               IF (i /= xe) THEN
                  ymat(k,kk,i,j) = c*gya(k,kk)
               END IF
            END DO
         END DO

         ! Start to formulate Jacobian matrix, diagonal blocks
         ieq = ((i + (j-1)*nx) - 1)*ordsq
         DO l = 1, ordsq
            DO k = 1, ordsq
               jmat(ieq+k,ieq+l) = jmat(ieq+k,ieq+l) + w(n)*c*gaa(k,l)
            END DO
         END DO
         ! Off-diagonal contribution from Y
         IF (i /= xs) THEN
            DO jj = ys, j, incy
               DO ii = xs, i, incx
                  jeq = ((ii + (jj-1)*nx) - 1)*ordsq
                  DO k = 1, ordsq
                     DO kk = 1, ordsq
                        tempy = 0.0
                        DO l = 1, order
                           tempy = tempy + w(n)*gay(k,l)*yold(l,kk,ii,jj)
                        END DO
                        jmat(ieq+k,jeq+kk) = jmat(ieq+k,jeq+kk) + tempy
                     END DO
                  END DO
               END DO
            END DO
         END IF
         ! Off-diagonal contribution from X
         IF (j /= ys) THEN
            DO jj = ys, j, incy
               DO ii = xs, i, incx
                  jeq = ((ii + (jj-1)*nx) - 1)*ordsq
                  DO k = 1, ordsq
                     DO kk = 1, ordsq
                        tempx = 0.0
                        DO l = 1, order
                           tempx = tempx + w(n)*gax(k,l)*xold(l,kk,ii,jj)
                        END DO
                        jmat(ieq+k,jeq+kk) = jmat(ieq+k,jeq+kk) + tempx
                     END DO
                  END DO
               END DO
            END DO
         END IF

      ! End of the mesh sweep
      END DO   ! Rows
   END DO   ! Columns
   END DO   ! Quadrants
! End the loop over all angles
END DO

! Begin constructing the RHS
DO k = 0, lambda
   DO l = 0, lambda
      DO j = 1, nx
         DO i = 1, ny
            ieq = ((i-1) + (j-1)*nx)*ordsq + (l+1) + k*order
            src(ieq) = s(i,j,k,l,g)/sigs((mat(i,j)),g,g)
         END DO
      END DO
   END DO
END DO
! Multiply the source vector by the Jacobi iteration matrix
q = MATMUL(jmat,src)

! Form the LHS matrix
jmat = -jmat
DO k = 1, neq
   jmat(k,k) = 1.0 + jmat(k,k)
END DO

! Symmetrize the matrix
DO j = 1, ny
   DO i = 1, nx
      m = mat(i,j)
      symm = sigs(m,g,g)*dx(i)*dy(j)
      DO k = 0, lambda
         DO l = 0, lambda
            IF (order == 1) symm2 = symm
            IF (order /= 1) symm2 = symm*(2.0*(k+1)-1.0)*(2.0*(l+1)-1.0)
            ndx = ((j-1)*nx + (i-1))*ordsq + (l+1) + k*order
            q(ndx) = symm2*q(ndx)
            DO jj = 1, neq
               jmat(ndx,jj) = symm2*jmat(ndx,jj)
            END DO
         END DO
      END DO
   END DO
END DO

! Determine if the matrix will be put into a file. If so, do it.
READ(7,*) matrix
IF (matrix == 1 .AND. g == 1) THEN
   OPEN (UNIT = 14, FILE = "jmat")
   WRITE(14,*) "! File containing the matrix, RHS, and solution -- GROUP 1 ONLY"
   WRITE(14,*) "! Problem scope: I J lambda, Max #its, Conv. criterion"
   WRITE(14,'(2X,I5,1X,I5,1X,I2,1X,I5,1X,ES9.3)') nx, ny, lambda, itmx, err
   WRITE(14,*)
   WRITE(14,*) "! JMAT-symmetric printed row-wise (i.e., j varies faster than i)"
   DO i = 1, neq
      WRITE(14,116) (jmat(i,j), j = 1, neq)
   END DO
   WRITE(14,*)
   WRITE(14,*) "! RHS vector printed first to last"
   WRITE(14,116) (q(i), i = 1, neq)
   WRITE(14,*)
   WRITE(14,*) "! Solution printed first to last in vector"
END IF
116 FORMAT(2X,8ES14.6)

! Get the time to construct the matrix
! Time will be higher if required to copy matrix to file
CALL CPU_TIME(tjmat)

READ(7,*) itmflag
IF (itmflag == 1) THEN
   ! Try a conjugate gradient solver
   CALL cgsolve(g,q,sol)
   ! Properly place the solution in the f-matrix
   DO k = 0, lambda
      DO l = 0, lambda
         DO j = 1, ny
            DO i = 1, nx
               ieq = ((i-1) + (j-1)*nx)*ordsq + (l+1) + k*order
               f(i,j,k,l,g) = sol(ieq)
            END DO
         END DO
      END DO
   END DO

   ! Write the solution to separate file if requested
   IF (matrix == 1 .AND. g == 1) WRITE(14,116) (sol(i), i = 1, neq)

ELSE
   IF (itmflag /= 0) THEN
      WRITE (8,'(/,1X,A)') "WARNING: Flag for ITM solution type not 0 or 1, direct solver used by default."
      warn = warn + 1
   END IF
   ! Solve the matrix system for the converged scalar flux soluton
   ! Asymmetric solver
     ! CALL dgesv(neq,1,jmat,neq,ipiv,q,neq,info)
   ! Symmetric solver
   CALL dposv('U',neq,1,jmat,neq,q,neq,info)
   IF (info /= 0) THEN
      WRITE (8, '(//,1X,A)') "ERROR: matrix either has illegal value or is singular."
      STOP
   END IF

      ! Use the LINPACK routines for solving the system
      !CALL dgeco(jmat,neq,neq,ipiv,rcond,wrk)
      !IF (abs(rcond) < 1.0e-10) THEN
      !   WRITE (8,'(//,1X,A)') "WARNING: rcond very small, large condition number"
      !   warn = warn + 1
      !END IF
      !info = 0
      !CALL dgesl(jmat,neq,neq,ipiv,q,info)

   cnvf(g) = 1
   ! Properly place the solution in the f-matrix
   DO k = 0, lambda
      DO l = 0, lambda
         DO j = 1, ny
            DO i = 1, nx
               ieq = ((i-1) + (j-1)*nx)*ordsq + (l+1) + k*order
               f(i,j,k,l,g) = q(ieq)
            END DO
         END DO
      END DO
   END DO

   ! Write the solution to separate file if requested
   IF (matrix == 1) WRITE(14,116) (q(i), i = 1, neq)

END IF

! Set the convergence flag so output subroutine will print solution
! cnvf(g) = 1

RETURN
END SUBROUTINE jima


text:
      SUBROUTINE ztssrch (IFLTAB, CINPATH, IUNIT, COUTPATH, NFOUND)
C
      implicit none
      integer zdssversion
C
C     Time Series Search.
C     Given a DSS file IFLTAB and a time series pathname,
C     search the database for data sets with the same pathname,
C     except for the "D" (Date) part.  Either the first pathname
C     found (same except for the D part) can be returned (FAST),
C     or a list of all pathnames that match can be written to
C     a file (complete but slow, depending on the file size).
C     This routine is typically used when you have parts of the pathname
C     but you do not know when the data occurs.
C
C     IFLTAB - (input) Array from the zopen6 call.
C     CINPATH - (input) Character string containing the pathname to
C               search.  The D part is ignored (usually ts blank).
C     IUNIT - (input)  If a complete list is to be found, this should
C               be the unit number of the file to write the list to.
C               If only the first pathname is to be returned, set this
C               to zero.
C     COUTPATH - (output) Character string containing the pathname
C               of the first matching path, if IUNIT is zero.
C               If IUNIT is positive, this will be blank.
C     NFOUND - (output)  The number of pathnames found.  -1 if an error
C               occurred, zero if none.
C
C
C     Argument Dimensions
      CHARACTER CINPATH*(*), COUTPATH*(*)
      INTEGER IFLTAB(*)
      INTEGER IUNIT, NFOUND
C
      IF (zdssVersion(IFLTAB).EQ.6) THEN
         CALL ztssrch6 (IFLTAB, CINPATH, IUNIT, COUTPATH, NFOUND)
      ELSE
         CALL ztssrch7 (IFLTAB, CINPATH, IUNIT, COUTPATH, NFOUND)
      ENDIF
C
      RETURN
      END

max_stars_repo_name: HydrologicEngineeringCenter/heclib

max_issues_repo_name: HydrologicEngineeringCenter/heclib

max_forks_repo_name: HydrologicEngineeringCenter/heclib