text:
module quadpack

  use, intrinsic :: iso_c_binding, only: c_double
  private dqelg, dqk21, dqpsrt
  public dqagse

  abstract interface
     ! f: real(c_double) --> real(c_double)
     real(c_double) function scalar_func(x)
       use, intrinsic :: iso_c_binding, only: c_double
       implicit none

       real(c_double), intent(in) :: x
     end function scalar_func
  end interface

contains

subroutine dqagse(f,a,b,epsabs,epsrel,limit,result,abserr,neval, &
     ier,alist,blist,rlist,elist,iord,last)

!*****************************************************************************80
!
!! DQAGSE estimates the integral of a function.
!
! D - double precision
! Q - quadrature
! A - adaptive
! G - General integrand (i.e. INT f(x), not weighted INT w(x) f(x))
! S - Singularities handled
! E - Extended
!
!  Modified:
!
!    11 September 2015
!
!  Author:
!
!    Robert Piessens, Elise de Doncker
!
!***purpose  the routine calculates an approximation result to a given
!      definite integral i = integral of f over (a,b),
!      hopefully satisfying following claim for accuracy
!      abs(i-result).le.max(epsabs,epsrel*abs(i)).
!
!  Parameters:
!
!   on entry
!      f      - real ( kind = 8 )
!               function subprogram defining the integrand
!               function f(x). the actual name for f needs to be
!               declared e x t e r n a l in the driver program.
!
!      a      - real ( kind = 8 )
!               lower limit of integration
!
!      b      - real ( kind = 8 )
!               upper limit of integration
!
!      epsabs - real ( kind = 8 )
!               absolute accuracy requested
!      epsrel - real ( kind = 8 )
!               relative accuracy requested
!               if  epsabs.le.0
!               and epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
!               the routine will end with ier = 6.
!
!      limit  - integer ( kind = 4 )
!               gives an upperbound on the number of subintervals
!               in the partition of (a,b)
!
!   on return
!      result - real ( kind = 8 )
!               approximation to the integral
!
!      abserr - real ( kind = 8 )
!               estimate of the modulus of the absolute error,
!               which should equal or exceed abs(i-result)
!
!      neval  - integer ( kind = 4 )
!               number of integrand evaluations
!
!      ier    - integer ( kind = 4 )
!               ier = 0 normal and reliable termination of the
!                       routine. it is assumed that the requested
!                       accuracy has been achieved.
!               ier.gt.0 abnormal termination of the routine
!                       the estimates for integral and error are
!                       less reliable. it is assumed that the
!                       requested accuracy has not been achieved.
!      error messages
!                   = 1 maximum number of subdivisions allowed
!                       has been achieved. one can allow more sub-
!                       divisions by increasing the value of limit
!                       (and taking the according dimension
!                       adjustments into account). however, if
!                       this yields no improvement it is advised
!                       to analyze the integrand in order to
!                       determine the integration difficulties. if
!                       the position of a local difficulty can be
!                       determined (e.g. singularity,
!                       discontinuity within the interval) one
!                       will probably gain from splitting up the
!                       interval at this point and calling the
!                       integrator on the subranges. if possible,
!                       an appropriate special-purpose integrator
!                       should be used, which is designed for
!                       handling the type of difficulty involved.
!                   = 2 the occurrence of roundoff error is detec-
!                       ted, which prevents the requested
!                       tolerance from being achieved.
!                       the error may be under-estimated.
!                   = 3 extremely bad integrand behaviour
!                       occurs at some points of the integration
!                       interval.
!                   = 4 the algorithm does not converge.
!                       roundoff error is detected in the
!                       extrapolation table.
!                       it is presumed that the requested
!                       tolerance cannot be achieved, and that the
!                       returned result is the best which can be
!                       obtained.
!                   = 5 the integral is probably divergent, or
!                       slowly convergent. it must be noted that
!                       divergence can occur with any other value
!                       of ier.
!                   = 6 the input is invalid, because
!                       epsabs.le.0 and
!                       epsrel.lt.max(50*rel.mach.acc.,0.5d-28).
!                       result, abserr, neval, last, rlist(1),
!                       iord(1) and elist(1) are set to zero.
!                       alist(1) and blist(1) are set to a and b
!                       respectively.
!
!      alist  - real ( kind = 8 )
!               vector of dimension at least limit, the first
!                last  elements of which are the left end points
!               of the subintervals in the partition of the
!               given integration range (a,b)
!
!      blist  - real ( kind = 8 )
!               vector of dimension at least limit, the first
!                last  elements of which are the right end points
!               of the subintervals in the partition of the given
!               integration range (a,b)
!
!      rlist  - real ( kind = 8 )
!               vector of dimension at least limit, the first
!                last  elements of which are the integral
!               approximations on the subintervals
!
!      elist  - real ( kind = 8 )
!               vector of dimension at least limit, the first
!                last  elements of which are the moduli of the
!               absolute error estimates on the subintervals
!
!      iord   - integer ( kind = 4 )
!               vector of dimension at least limit, the first k
!               elements of which are pointers to the
!               error estimates over the subintervals,
!               such that elist(iord(1)), ..., elist(iord(k))
!               form a decreasing sequence, with k = last
!               if last.le.(limit/2+2), and k = limit+1-last
!               otherwise
!
!      last   - integer ( kind = 4 )
!               number of subintervals actually produced in the
!               subdivision process
!
!  Local parameters:
!
!      the dimension of rlist2 is determined by the value of
!      limexp in routine dqelg (rlist2 should be of dimension
!      (limexp+2) at least).
!
!      list of major variables
!
!     alist     - list of left end points of all subintervals
!                 considered up to now
!     blist     - list of right end points of all subintervals
!                 considered up to now
!     rlist(i)  - approximation to the integral over
!                 (alist(i),blist(i))
!     rlist2    - array of dimension at least limexp+2 containing
!                 the part of the epsilon table which is still
!                 needed for further computations
!     elist(i)  - error estimate applying to rlist(i)
!     maxerr    - pointer to the interval with largest error
!                 estimate
!     errmax    - elist(maxerr)
!     erlast    - error on the interval currently subdivided
!                 (before that subdivision has taken place)
!     area      - sum of the integrals over the subintervals
!     errsum    - sum of the errors over the subintervals
!     errbnd    - requested accuracy max(epsabs,epsrel*
!                 abs(result))
!     *****1    - variable for the left interval
!     *****2    - variable for the right interval
!     last      - index for subdivision
!     nres      - number of calls to the extrapolation routine
!     numrl2    - number of elements currently in rlist2. if an
!                 appropriate approximation to the compounded
!                 integral has been obtained it is put in
!                 rlist2(numrl2) after numrl2 has been increased
!                 by one.
!     small     - length of the smallest interval considered up
!                 to now, multiplied by 1.5
!     erlarg    - sum of the errors over the intervals larger
!                 than the smallest interval considered up to now
!     extrap    - logical variable denoting that the routine is
!                 attempting to perform extrapolation i.e. before
!                 subdividing the smallest interval we try to
!                 decrease the value of erlarg.
!     noext     - logical variable denoting that extrapolation
!                 is no longer allowed (true value)
!
!      machine dependent constants
!
!     epmach is the largest relative spacing.
!     uflow is the smallest positive magnitude.
!     oflow is the largest positive magnitude.
!
  implicit none

  real ( kind = 8 ) a,abseps,abserr,alist,area,area1,area12,area2,a1, &
    a2,b,blist,b1,b2,correc,defabs,defab1,defab2, &
    dres,elist,epmach,epsabs,epsrel,erlarg,erlast,errbnd,errmax, &
    error1,error2,erro12,errsum,ertest,oflow,resabs,reseps,result, &
    res3la,rlist,rlist2,small,uflow
  integer ( kind = 4 ) id,ier,ierro,iord,iroff1,iroff2,iroff3,jupbnd, &
    k,ksgn,ktmin,last,limit,maxerr,neval,nres,nrmax,numrl2
  logical extrap,noext
  dimension alist(limit),blist(limit),elist(limit),iord(limit), &
   res3la(3),rlist(limit),rlist2(52)

  procedure(scalar_func) :: f

  correc = 0.0D+00
  erlarg = 0.0D+00
  ertest = 0.0D+00
  small = 0.0D+00
  epmach = epsilon ( epmach )
!
!  test on validity of parameters
!
  ier = 0
  neval = 0
  last = 0
  result = 0.0D+00
  abserr = 0.0D+00
  alist(1) = a
  blist(1) = b
  rlist(1) = 0.0D+00
  elist(1) = 0.0D+00

  if(epsabs.le.0.0D+00.and.epsrel.lt. max ( 0.5D+02*epmach,0.5d-28)) then
    ier = 6
    return
  end if
!
!  first approximation to the integral
!
  uflow = tiny ( uflow )
  oflow = huge ( oflow )
  ierro = 0
  call dqk21(f,a,b,result,abserr,defabs,resabs)
!
!  test on accuracy.
!
  dres =  abs ( result)
  errbnd =  max ( epsabs,epsrel*dres)
  last = 1
  rlist(1) = result
  elist(1) = abserr
  iord(1) = 1
  if(abserr.le.1.0D+02*epmach*defabs.and.abserr.gt.errbnd) ier = 2
  if(limit.eq.1) ier = 1
  if(ier.ne.0.or.(abserr.le.errbnd.and.abserr.ne.resabs).or. &
    abserr.eq.0.0D+00) go to 140
!
!  initialization
!
  rlist2(1) = result
  errmax = abserr
  maxerr = 1
  area = result
  errsum = abserr
  abserr = oflow
  nrmax = 1
  nres = 0
  numrl2 = 2
  ktmin = 0
  extrap = .false.
  noext = .false.
  iroff1 = 0
  iroff2 = 0
  iroff3 = 0
  ksgn = -1
  if(dres.ge.(0.1D+01-0.5D+02*epmach)*defabs) ksgn = 1
!
!  main do-loop
!
  do 90 last = 2,limit
!
!  bisect the subinterval with the nrmax-th largest error estimate.
!
    a1 = alist(maxerr)
    b1 = 0.5D+00*(alist(maxerr)+blist(maxerr))
    a2 = b1
    b2 = blist(maxerr)
    erlast = errmax
    call dqk21(f,a1,b1,area1,error1,resabs,defab1)
    call dqk21(f,a2,b2,area2,error2,resabs,defab2)
!
!  improve previous approximations to integral
!  and error and test for accuracy.
!
    area12 = area1+area2
    erro12 = error1+error2
    errsum = errsum+erro12-errmax
    area = area+area12-rlist(maxerr)
    if(defab1.eq.error1.or.defab2.eq.error2) go to 15
    if( abs ( rlist(maxerr)-area12).gt.0.1D-04* abs ( area12) &
    .or.erro12.lt.0.99D+00*errmax) go to 10
    if(extrap) iroff2 = iroff2+1
    if(.not.extrap) iroff1 = iroff1+1
   10   if(last.gt.10.and.erro12.gt.errmax) iroff3 = iroff3+1
   15   rlist(maxerr) = area1
    rlist(last) = area2
    errbnd =  max ( epsabs,epsrel* abs ( area))
!
!  test for roundoff error and eventually set error flag.
!
    if(iroff1+iroff2.ge.10.or.iroff3.ge.20) ier = 2
    if(iroff2.ge.5) ierro = 3
!
!  set error flag in the case that the number of subintervals
!  equals limit.
!
    if(last.eq.limit) ier = 1
!
!  set error flag in the case of bad integrand behaviour
!  at a point of the integration range.
!
    if( max (  abs ( a1), abs ( b2)).le.(0.1D+01+0.1D+03*epmach)* &
    ( abs ( a2)+0.1D+04*uflow)) ier = 4
!
!  append the newly-created intervals to the list.
!
    if(error2.gt.error1) go to 20
    alist(last) = a2
    blist(maxerr) = b1
    blist(last) = b2
    elist(maxerr) = error1
    elist(last) = error2
    go to 30
   20   alist(maxerr) = a2
    alist(last) = a1
    blist(last) = b1
    rlist(maxerr) = area2
    rlist(last) = area1
    elist(maxerr) = error2
    elist(last) = error1
!
!  call dqpsrt to maintain the descending ordering
!  in the list of error estimates and select the subinterval
!  with nrmax-th largest error estimate (to be bisected next).
!
   30   call dqpsrt(limit,last,maxerr,errmax,elist,iord,nrmax)
    if(errsum.le.errbnd) go to 115
    if(ier.ne.0) go to 100
    if(last.eq.2) go to 80
    if(noext) go to 90
    erlarg = erlarg-erlast
    if( abs ( b1-a1).gt.small) erlarg = erlarg+erro12
    if(extrap) go to 40
!
!  test whether the interval to be bisected next is the
!  smallest interval.
!
    if( abs ( blist(maxerr)-alist(maxerr)).gt.small) go to 90
    extrap = .true.
    nrmax = 2
   40   if(ierro.eq.3.or.erlarg.le.ertest) go to 60
!
!  the smallest interval has the largest error.
!  before bisecting decrease the sum of the errors over the
!  larger intervals (erlarg) and perform extrapolation.
!
    id = nrmax
    jupbnd = last
    if(last.gt.(2+limit/2)) jupbnd = limit+3-last
    do k = id,jupbnd
      maxerr = iord(nrmax)
      errmax = elist(maxerr)
      if( abs ( blist(maxerr)-alist(maxerr)).gt.small) go to 90
      nrmax = nrmax+1
    end do
!
!  perform extrapolation.
!
   60   numrl2 = numrl2+1
    rlist2(numrl2) = area
    call dqelg(numrl2,rlist2,reseps,abseps,res3la,nres)
    ktmin = ktmin+1
    if(ktmin.gt.5.and.abserr.lt.0.1D-02*errsum) ier = 5
    if(abseps.ge.abserr) go to 70
    ktmin = 0
    abserr = abseps
    result = reseps
    correc = erlarg
    ertest =  max ( epsabs,epsrel* abs ( reseps))
    if(abserr.le.ertest) go to 100
!
!  prepare bisection of the smallest interval.
!
   70   if(numrl2.eq.1) noext = .true.
    if(ier.eq.5) go to 100
    maxerr = iord(1)
    errmax = elist(maxerr)
    nrmax = 1
    extrap = .false.
    small = small*0.5D+00
    erlarg = errsum
    go to 90
   80   small =  abs ( b-a)*0.375D+00
    erlarg = errsum
    ertest = errbnd
    rlist2(2) = area
   90 continue
!
!  set final result and error estimate.
!
  100 if(abserr.eq.oflow) go to 115
  if(ier+ierro.eq.0) go to 110
  if(ierro.eq.3) abserr = abserr+correc
  if(ier.eq.0) ier = 3
  if(result.ne.0.0D+00.and.area.ne.0.0D+00) go to 105
  if(abserr.gt.errsum) go to 115
  if(area.eq.0.0D+00) go to 130
  go to 110
  105 if(abserr/ abs ( result).gt.errsum/ abs ( area)) go to 115
!
!  test on divergence.
!
  110 if(ksgn.eq.(-1).and. max (  abs ( result), abs ( area)).le. &
   defabs*0.1D-01) go to 130
  if(0.1D-01.gt.(result/area).or.(result/area).gt.0.1D+03 &
   .or.errsum.gt. abs ( area)) ier = 6
  go to 130
!
!  compute global integral sum.
!
  115 result = 0.0D+00
  do k = 1,last
     result = result+rlist(k)
  end do
  abserr = errsum
  130 if(ier.gt.2) ier = ier-1
  140 neval = 42*last-21

  return
end subroutine dqagse

subroutine dqelg ( n, epstab, result, abserr, res3la, nres )

!*****************************************************************************80
!
!! DQELG carries out the Epsilon extrapolation algorithm.
!
!  Modified:
!
!    11 September 2015
!
!  Author:
!
!    Robert Piessens, Elise de Doncker
!
!***purpose  the routine determines the limit of a given sequence of
!      approximations, by means of the epsilon algorithm of
!      p.wynn. an estimate of the absolute error is also given.
!      the condensed epsilon table is computed. only those
!      elements needed for the computation of the next diagonal
!      are preserved.
!
!  Parameters:
!
!        n      - integer ( kind = 4 )
!                 epstab(n) contains the new element in the
!                 first column of the epsilon table.
!
!        epstab - real ( kind = 8 )
!                 vector of dimension 52 containing the elements
!                 of the two lower diagonals of the triangular
!                 epsilon table. the elements are numbered
!                 starting at the right-hand corner of the
!                 triangle.
!
!        result - real ( kind = 8 )
!                 resulting approximation to the integral
!
!        abserr - real ( kind = 8 )
!                 estimate of the absolute error computed from
!                 result and the 3 previous results
!
!        res3la - real ( kind = 8 )
!                 vector of dimension 3 containing the last 3
!                 results
!
!        nres   - integer ( kind = 4 )
!                 number of calls to the routine
!                 (should be zero at first call)
!
!  Local Parameters:
!
!     e0     - the 4 elements on which the computation of a new
!     e1       element in the epsilon table is based
!     e2
!     e3                 e0
!                  e3    e1    new
!                        e2
!     newelm - number of elements to be computed in the new
!              diagonal
!     error  - error = abs(e1-e0)+abs(e2-e1)+abs(new-e2)
!     result - the element in the new diagonal with least value
!              of error
!
!     machine dependent constants
!
!     epmach is the largest relative spacing.
!     oflow is the largest positive magnitude.
!     limexp is the maximum number of elements the epsilon
!     table can contain. if this number is reached, the upper
!     diagonal of the epsilon table is deleted.
!
  implicit none

  real ( kind = 8 ) abserr,delta1,delta2,delta3, &
    epmach,epsinf,epstab,error,err1,err2,err3,e0,e1,e1abs,e2,e3, &
    oflow,res,result,res3la,ss,tol1,tol2,tol3
  integer ( kind = 4 ) i,ib,ib2,ie,indx,k1,k2,k3,limexp,n,newelm
  integer ( kind = 4 ) nres
  integer ( kind = 4 ) num
  dimension epstab(52),res3la(3)

  epmach = epsilon ( epmach )
  oflow = huge ( oflow )
  nres = nres+1
  abserr = oflow
  result = epstab(n)
  if(n.lt.3) go to 100
  limexp = 50
  epstab(n+2) = epstab(n)
  newelm = (n-1)/2
  epstab(n) = oflow
  num = n
  k1 = n

  do 40 i = 1,newelm

    k2 = k1-1
    k3 = k1-2
    res = epstab(k1+2)
    e0 = epstab(k3)
    e1 = epstab(k2)
    e2 = res
    e1abs =  abs ( e1)
    delta2 = e2-e1
    err2 =  abs ( delta2)
    tol2 =  max (  abs ( e2),e1abs)*epmach
    delta3 = e1 - e0
    err3 =  abs ( delta3)
    tol3 =  max ( e1abs, abs ( e0))*epmach
    if(err2.gt.tol2.or.err3.gt.tol3) go to 10
!
!  if e0, e1 and e2 are equal to machine accuracy, convergence is assumed.
!
    result = res
    abserr = err2+err3
    go to 100
   10   e3 = epstab(k1)
    epstab(k1) = e1
    delta1 = e1-e3
    err1 =  abs ( delta1)
    tol1 =  max ( e1abs, abs ( e3))*epmach
!
!  if two elements are very close to each other, omit
!  a part of the table by adjusting the value of n
!
    if(err1.le.tol1.or.err2.le.tol2.or.err3.le.tol3) go to 20
    ss = 0.1D+01/delta1+0.1D+01/delta2-0.1D+01/delta3
    epsinf =  abs ( ss*e1)
!
!  test to detect irregular behaviour in the table, and
!  eventually omit a part of the table adjusting the value
!  of n.
!
    if(epsinf.gt.0.1D-03) go to 30
   20   n = i+i-1
    go to 50
!
!  compute a new element and eventually adjust
!  the value of result.
!
   30   res = e1+0.1D+01/ss
    epstab(k1) = res
    k1 = k1-2
    error = err2 + abs ( res-e2 ) + err3

    if ( error .le. abserr ) then
      abserr = error
      result = res
    end if

   40 continue
!
!  shift the table.
!
   50 if(n.eq.limexp) n = 2*(limexp/2)-1
  ib = 1
  if((num/2)*2.eq.num) ib = 2
  ie = newelm+1
  do i=1,ie
    ib2 = ib+2
    epstab(ib) = epstab(ib2)
    ib = ib2
  end do
  if(num.eq.n) go to 80
  indx = num-n+1
  do i = 1,n
    epstab(i)= epstab(indx)
    indx = indx+1
  end do
   80 if(nres.ge.4) go to 90
  res3la(nres) = result
  abserr = oflow
  go to 100
!
!  compute error estimate
!
   90 abserr =  abs ( result-res3la(3))+ abs ( result-res3la(2)) &
    + abs ( result-res3la(1))
  res3la(1) = res3la(2)
  res3la(2) = res3la(3)
  res3la(3) = result
  100 continue

  abserr =  max ( abserr, 0.5D+01*epmach* abs ( result))

  return
end subroutine dqelg

subroutine dqk21(f,a,b,result,abserr,resabs,resasc)

!*****************************************************************************80
!
!! DQK21 carries out a 21 point Gauss-Kronrod quadrature rule.
!
!  Modified:
!
!    11 September 2015
!
!  Author:
!
!    Robert Piessens, Elise de Doncker
!
!***purpose  to compute i = integral of f over (a,b), with error
!                     estimate
!                 j = integral of abs(f) over (a,b)
!
!  Parameters:
!
!      on entry
!        f      - real ( kind = 8 )
!                 function subprogram defining the integrand
!                 function f(x). the actual name for f needs to be
!                 declared e x t e r n a l in the driver program.
!
!        a      - real ( kind = 8 )
!                 lower limit of integration
!
!        b      - real ( kind = 8 )
!                 upper limit of integration
!
!      on return
!        result - real ( kind = 8 )
!                 approximation to the integral i
!                 result is computed by applying the 21-point
!                 kronrod rule (resk) obtained by optimal addition
!                 of abscissae to the 10-point gauss rule (resg).
!
!        abserr - real ( kind = 8 )
!                 estimate of the modulus of the absolute error,
!                 which should not exceed abs(i-result)
!
!        resabs - real ( kind = 8 )
!                 approximation to the integral j
!
!        resasc - real ( kind = 8 )
!                 approximation to the integral of abs(f-i/(b-a))
!                 over (a,b)
!
!  Local Parameters:
!
!
!     the abscissae and weights are given for the interval (-1,1).
!     because of symmetry only the positive abscissae and their
!     corresponding weights are given.
!
!     xgk    - abscissae of the 21-point kronrod rule
!              xgk(2), xgk(4), ...  abscissae of the 10-point
!              gauss rule
!              xgk(1), xgk(3), ...  abscissae which are optimally
!              added to the 10-point gauss rule
!
!     wgk    - weights of the 21-point kronrod rule
!
!     wg     - weights of the 10-point gauss rule
!
!
! gauss quadrature weights and kronron quadrature abscissae and weights
! as evaluated with 80 decimal digit arithmetic by l. w. fullerton,
! bell labs, nov. 1981.
!
!     centr  - mid point of the interval
!     hlgth  - half-length of the interval
!     absc   - abscissa
!     fval*  - function value
!     resg   - result of the 10-point gauss formula
!     resk   - result of the 21-point kronrod formula
!     reskh  - approximation to the mean value of f over (a,b),
!              i.e. to i/(b-a)
!
!
!     machine dependent constants
!
!     epmach is the largest relative spacing.
!     uflow is the smallest positive magnitude.
!
  implicit none

  real ( kind = 8 ) a,absc,abserr,b,centr,dhlgth, &
    epmach,fc,fsum,fval1,fval2,fv1,fv2,hlgth,resabs,resasc, &
    resg,resk,reskh,result,uflow,wg,wgk,xgk
  integer ( kind = 4 ) j,jtw,jtwm1
  procedure(scalar_func) :: f
  dimension fv1(10),fv2(10),wg(5),wgk(11),xgk(11)

  data wg  (  1) / 0.066671344308688137593568809893332d0 /
  data wg  (  2) / 0.149451349150580593145776339657697d0 /
  data wg  (  3) / 0.219086362515982043995534934228163d0 /
  data wg  (  4) / 0.269266719309996355091226921569469d0 /
  data wg  (  5) / 0.295524224714752870173892994651338d0 /

  data xgk (  1) / 0.995657163025808080735527280689003d0 /
  data xgk (  2) / 0.973906528517171720077964012084452d0 /
  data xgk (  3) / 0.930157491355708226001207180059508d0 /
  data xgk (  4) / 0.865063366688984510732096688423493d0 /
  data xgk (  5) / 0.780817726586416897063717578345042d0 /
  data xgk (  6) / 0.679409568299024406234327365114874d0 /
  data xgk (  7) / 0.562757134668604683339000099272694d0 /
  data xgk (  8) / 0.433395394129247190799265943165784d0 /
  data xgk (  9) / 0.294392862701460198131126603103866d0 /
  data xgk ( 10) / 0.148874338981631210884826001129720d0 /
  data xgk ( 11) / 0.000000000000000000000000000000000d0 /

  data wgk (  1) / 0.011694638867371874278064396062192d0 /
  data wgk (  2) / 0.032558162307964727478818972459390d0 /
  data wgk (  3) / 0.054755896574351996031381300244580d0 /
  data wgk (  4) / 0.075039674810919952767043140916190d0 /
  data wgk (  5) / 0.093125454583697605535065465083366d0 /
  data wgk (  6) / 0.109387158802297641899210590325805d0 /
  data wgk (  7) / 0.123491976262065851077958109831074d0 /
  data wgk (  8) / 0.134709217311473325928054001771707d0 /
  data wgk (  9) / 0.142775938577060080797094273138717d0 /
  data wgk ( 10) / 0.147739104901338491374841515972068d0 /
  data wgk ( 11) / 0.149445554002916905664936468389821d0 /

  epmach = epsilon ( epmach )
  uflow = tiny ( uflow )
  centr = 0.5D+00*(a+b)
  hlgth = 0.5D+00*(b-a)
  dhlgth =  abs ( hlgth)
!
!  compute the 21-point kronrod approximation to
!  the integral, and estimate the absolute error.
!
  resg = 0.0D+00
  fc = f(centr)
  resk = wgk(11)*fc
  resabs =  abs ( resk)
  do j=1,5
    jtw = 2*j
    absc = hlgth*xgk(jtw)
    fval1 = f(centr-absc)
    fval2 = f(centr+absc)
    fv1(jtw) = fval1
    fv2(jtw) = fval2
    fsum = fval1+fval2
    resg = resg+wg(j)*fsum
    resk = resk+wgk(jtw)*fsum
    resabs = resabs+wgk(jtw)*( abs ( fval1)+ abs ( fval2))
  end do

  do j = 1,5
    jtwm1 = 2*j-1
    absc = hlgth*xgk(jtwm1)
    fval1 = f(centr-absc)
    fval2 = f(centr+absc)
    fv1(jtwm1) = fval1
    fv2(jtwm1) = fval2
    fsum = fval1+fval2
    resk = resk+wgk(jtwm1)*fsum
    resabs = resabs+wgk(jtwm1)*( abs ( fval1)+ abs ( fval2))
  end do

  reskh = resk*0.5D+00
  resasc = wgk(11)* abs ( fc-reskh)

  do j=1,10
    resasc = resasc+wgk(j)*( abs ( fv1(j)-reskh)+ abs ( fv2(j)-reskh))
  end do

  result = resk*hlgth
  resabs = resabs*dhlgth
  resasc = resasc*dhlgth
  abserr =  abs ( (resk-resg)*hlgth)
  if(resasc.ne.0.0D+00.and.abserr.ne.0.0D+00) &
    abserr = resasc*min(0.1D+01,(0.2D+03*abserr/resasc)**1.5D+00)
  if(resabs.gt.uflow/(0.5D+02*epmach)) abserr = max &
    ((epmach*0.5D+02)*resabs,abserr)

  return
end subroutine dqk21

subroutine dqpsrt ( limit, last, maxerr, ermax, elist, iord, nrmax )

!*****************************************************************************80
!
!! DQPSRT maintains the order of a list of local error estimates.
!
!  Modified:
!
!    11 September 2015
!
!  Author:
!
!    Robert Piessens, Elise de Doncker
!
!***purpose  this routine maintains the descending ordering in the
!      list of the local error estimated resulting from the
!      interval subdivision process. at each call two error
!      estimates are inserted using the sequential search
!      method, top-down for the largest error estimate and
!      bottom-up for the smallest error estimate.
!
!  Parameters:
!
!        limit  - integer ( kind = 4 )
!                 maximum number of error estimates the list
!                 can contain
!
!        last   - integer ( kind = 4 )
!                 number of error estimates currently in the list
!
!        maxerr - integer ( kind = 4 )
!                 maxerr points to the nrmax-th largest error
!                 estimate currently in the list
!
!        ermax  - real ( kind = 8 )
!                 nrmax-th largest error estimate
!                 ermax = elist(maxerr)
!
!        elist  - real ( kind = 8 )
!                 vector of dimension last containing
!                 the error estimates
!
!        iord   - integer ( kind = 4 )
!                 vector of dimension last, the first k elements
!                 of which contain pointers to the error
!                 estimates, such that
!                 elist(iord(1)),...,  elist(iord(k))
!                 form a decreasing sequence, with
!                 k = last if last.le.(limit/2+2), and
!                 k = limit+1-last otherwise
!
!        nrmax  - integer ( kind = 4 )
!                 maxerr = iord(nrmax)
!
  implicit none

  real ( kind = 8 ) elist,ermax,errmax,errmin
  integer ( kind = 4 ) i,ibeg,ido,iord,isucc,j,jbnd,jupbn,k,last
  integer ( kind = 4 ) limit
  integer ( kind = 4 ) maxerr
  integer ( kind = 4 ) nrmax
  dimension elist(last),iord(last)
!
!  check whether the list contains more than
!  two error estimates.
!
  if(last.gt.2) go to 10
  iord(1) = 1
  iord(2) = 2
  go to 90
!
!  this part of the routine is only executed if, due to a
!  difficult integrand, subdivision increased the error
!  estimate. in the normal case the insert procedure should
!  start after the nrmax-th largest error estimate.
!
   10 errmax = elist(maxerr)

  ido = nrmax-1
  do i = 1,ido
    isucc = iord(nrmax-1)
    if(errmax.le.elist(isucc)) go to 30
    iord(nrmax) = isucc
    nrmax = nrmax-1
  end do
!
!  compute the number of elements in the list to be maintained
!  in descending order. this number depends on the number of
!  subdivisions still allowed.
!
   30 jupbn = last
  if(last.gt.(limit/2+2)) jupbn = limit+3-last
  errmin = elist(last)
!
!  insert errmax by traversing the list top-down,
!  starting comparison from the element elist(iord(nrmax+1)).
!
  jbnd = jupbn-1
  ibeg = nrmax+1

  do i=ibeg,jbnd
    isucc = iord(i)
    if(errmax.ge.elist(isucc)) go to 60
    iord(i-1) = isucc
  end do

  iord(jbnd) = maxerr
  iord(jupbn) = last
  go to 90
!
!  insert errmin by traversing the list bottom-up.
!
   60 iord(i-1) = maxerr
  k = jbnd

  do j=i,jbnd
    isucc = iord(k)
    if(errmin.lt.elist(isucc)) go to 80
    iord(k+1) = isucc
    k = k-1
  end do

  iord(i) = last
  go to 90
   80 iord(k+1) = last
!
!     set maxerr and ermax.
!
   90 maxerr = iord(nrmax)
  ermax = elist(maxerr)

  return
end subroutine dqpsrt

end module quadpack

max_stars_repo_name: kianmeng/bezier

max_issues_repo_name: kianmeng/bezier

max_forks_repo_name: kianmeng/bezier


text:
      SUBROUTINE SMLEIG1(D,O,VAL)
C
C     COMPUTES EIGENVALUES AND VECTORS FOR 1X1 AND 2X2
C
      REAL    D(2),O(2),VAL(2),P,Q
      INTEGER ENTRY,XENTRY,SYSBUF,PHIA,MCB(7)
      DIMENSION VCOM(30)
C
      COMMON/SYSTEM/SYSBUF
      COMMON /GIVN / TITLE(150)
      COMMON /PACKX/IT1,IT2,II,JJ,INCR
      COMMON /UNPAKX/IT3,III,JJJ,INCR1
C
      EQUIVALENCE
     1 (MO,TITLE(2)),(MD,TITLE(3)),(ENTRY,TITLE(11)),(XENTRY,TITLE(20)),
     2 (VCOM(1),TITLE(101)),(N,VCOM(1)),(LAMA,VCOM(6)),(PHIA,VCOM(12)),
     3 (NFOUND,VCOM(10))
C
      DATA MCB/7*0/
C
C     D        ARRAY OF DIAGONALS
C     O        ARRAY OF OFF DIAGONALS
C     VAL      ARRAY OF EIGENVALUES
C     LAMA     FILE OF EIGENVALUES--HEADER,VALUES,ORDER FOUND
C     PHIA     FILE OF  VECTORS   --......,VECTORS-D.P.
C     MO       RESTART TAPE FOR MORE EIGENVALUES
C     MD       INPUT MATRIX
C     N        ORDER OF  PROBLEM
C     NFOUND   NUMBER OF EIGENVALUES/VECTOR PREVIOUSLY FOUND
C
CWKBR 2/94      IBUF1 =(KORSZ(O) - SYSBUF +1 )/2  -1
      IBUF1 =(KORSZ(O) - SYSBUF +1 )  -1
C
C     OPEN INPUT MATRIX
C
      CALL GOPEN(MD,O(IBUF1),0)
C
C     SETUP FOR UNPACK
C
      IT3 = 1
      III = 1
      JJJ = N
      INCR1= 1
      ASSIGN 101 TO ITRA
      CALL UNPACK(*1000,MD,D)
  101 IF(N .EQ. 2) GO TO 110
C
C     THE MATRIX IS A 1X1
C
      O(1) = 0.0
      VAL(1) = D(1)
      LOC = 1
      GO TO 120
C
C     THE MATRIX IS A 2X2
C
  110 O(1) = D(2)
      O(2) = 0.0
      ASSIGN 111 TO ITRA
      III = 2
      CALL UNPACK(*1000,MD,D(2))
  111 P = D(1) + D(2)
      Q = SQRT( P*P -4.0*(D(1)*D(2)- O(1)**2))
      VAL(1) =(P+Q)/2.0
      VAL(2) = (P-Q)/2.0
      LOC = 0
C
C     WRAP UP ROUTINE
C
  120 CALL CLOSE(MD,1)
C
C     COPY D,O,LOC ONTO MO FOR RESTART
C
      CALL GOPEN(MO,O(IBUF1),1)
C
C     SETUP FOR PACK
C
      IM1=1
      IT1=1
      IT2=1
      II =1
      JJ = N
      INCR =1
      CALL PACK(D,MO,MCB)
      CALL PACK(O,MO,MCB)
      CALL WRITE(MO,LOC,1,1)
      CALL CLOSE(MO,1)
      IF(N .NE. 1) GO TO 125
C
C     1X1 WRITE OUT VECTORS AND VALUES
C
      MCB(1) = PHIA
      MCB(2) = 0
      MCB(3) = 1
      MCB(4) = 2
      MCB(5) = 2
      MCB(6) = 0
      CALL GOPEN(PHIA,O(IBUF1),  1)
      JJ = 1
      CALL PACK(1.0,PHIA,MCB)
      CALL CLOSE(PHIA,1)
      CALL WRTTRL(MCB(1))
      CALL GOPEN(LAMA,O(IBUF1),1 )
      IF (NFOUND .EQ. 0) GO TO 128
      DO 126 I= 1,NFOUND
      CALL WRITE(LAMA,0.0,1,0)
  126 CONTINUE
  128 VALX = VAL(1)
      CALL WRITE(LAMA,VALX,1,1)
      IF (NFOUND .EQ. 0) GO TO 129
      DO 127 I= 1,NFOUND
      CALL WRITE(LAMA,I,1,0)
  127 CONTINUE
  129 CALL WRITE(LAMA,NFOUND+1,1,1)
      CALL CLOSE(LAMA,1)
      MCB(1) = LAMA
      CALL WRTTRL(MCB)
  125 XENTRY = -ENTRY
      RETURN
 1000 DO 1001 I =III,JJJ
      D(I) = 0.0
 1001 CONTINUE
      GO TO  ITRA,(111,101)
      END

ext: f


text:
This editor can edit this entry and tell us a bit about themselves by clicking the Edit icon.

Welcome to the Wiki Howdy, David, and Welcome to the Wiki!   I reverted your edit to Users/Gabriella because you didnt leave any reason as to why you were editing her personal entry.  Normally editing the personal entry of another editor in such a manner is frowned on.  Its possible she asked you to edit it, which would make it okay, but you didnt mention any reason when you deleted a bunch of her entry.  It might be a good idea to note the reason for your edit in the edit comment (the field just below the editor, just above the Save Changes button).  Basically, right now it is unclear if you were vandalizing her entry (possibly just kidding around), or if youre acting on her behalf.   If its something she wants, please feel free to do the edit again (and let people know that youre doing it with her consent... the reason your edit was undone was to protect the entry she created).  In general, however, welcome to the wiki, and I hope you enjoy the community!  Users/JabberWokky Evan JabberWokky Edwards

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


! { dg-do compile }
!
! PR fortran/51809
!
! Contributed by Kacper Kowalik
!
module foo
   implicit none

   type foo_t
   contains
      procedure :: func_foo
   end type foo_t

contains

   subroutine func_foo(this)
      implicit none
      class(foo_t), intent(in) :: this
   end subroutine func_foo

end module foo

module bar
   use foo,   only: foo_t

   implicit none

   type, extends(foo_t) :: bar_t
   contains
      procedure :: func_bar
   end type bar_t

contains

   subroutine func_bar(this)
      use foo,    only: foo_t     ! <--- removing this line also fixes ICE
      implicit none
      class(bar_t), intent(in) :: this
   end subroutine func_bar

end module bar

module merry_ICE
   use foo,  only: foo_t   ! <------ change order to prevent ICE
   use bar,  only: bar_t   ! <------ change order to prevent ICE
end module merry_ICE


text:
! DART software - Copyright UCAR. This open source software is provided
! by UCAR, "as is", without charge, subject to all terms of use at
! http://www.image.ucar.edu/DAReS/DART/DART_download
!
! $Id$

program convert_madis_acars

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!   convert_madis_acars - program that reads a MADIS netCDF ACARS
!                         observation file and writes a DART
!                         obs_seq file using the DART library routines.
!
!     created Dec. 2007 Ryan Torn, NCAR/MMM
!     modified Dec. 2008 Soyoung Ha and David Dowell, NCAR/MMM
!     - added dewpoint as an output variable
!     - added relative humidity as an output variable
!
!     modified to use a common set of utilities, better netcdf error checks,
!     able to insert obs with any time correctly (not only monotonically
!     increasing times)    nancy collins,  ncar/image   11 march 2010
!
!     keep original obs times, make source for all converters as similar
!     as possbile.   nancy collins,  ncar/image   26 march 2010
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

use         types_mod, only : r8, missing_r8
use      location_mod, only : VERTISPRESSURE
use     utilities_mod, only : initialize_utilities, finalize_utilities
use  netcdf_utilities_mod, only : nc_open_file_readonly, nc_close_file
use  time_manager_mod, only : time_type, set_calendar_type, set_date, operator(>=), &
                              increment_time, get_time, operator(-), GREGORIAN
use  obs_sequence_mod, only : obs_sequence_type, obs_type, read_obs_seq, &
                              static_init_obs_sequence, init_obs, write_obs_seq, &
                              init_obs_sequence, get_num_obs, &
                              set_copy_meta_data, set_qc_meta_data
use        meteor_mod, only : pres_alt_to_pres, sat_vapor_pressure, &
                              specific_humidity, wind_dirspd_to_uv
use       obs_err_mod, only : acars_wind_error, acars_temp_error, &
                              acars_rel_hum_error
use dewpoint_obs_err_mod, only : dewpt_error_from_rh_and_temp, &
                                 rh_error_from_dewpt_and_temp
use          sort_mod, only : index_sort
use      obs_kind_mod, only : ACARS_U_WIND_COMPONENT, ACARS_V_WIND_COMPONENT, &
                              ACARS_TEMPERATURE, ACARS_SPECIFIC_HUMIDITY, &
                              ACARS_DEWPOINT, ACARS_RELATIVE_HUMIDITY
use obs_utilities_mod, only : getvar_real, get_or_fill_QC, add_obs_to_seq, &
                              create_3d_obs, getvar_int, getdimlen, set_missing_name

implicit none

character(len=14),  parameter :: acars_netcdf_file = 'acars_input.nc'
character(len=129), parameter :: acars_out_file    = 'obs_seq.acars'

! the following logical parameters control which water-vapor variables appear in the output file,
! whether to use the NCEP error or Lin and Hubbard (2004) moisture error model, and if the
! input file has data quality control fields, whether to use or ignore them.
logical, parameter :: LH_err                    = .false.
logical, parameter :: include_specific_humidity = .true.
logical, parameter :: include_relative_humidity = .false.
logical, parameter :: include_dewpoint          = .false.
logical, parameter :: use_input_qc              = .true.

integer, parameter :: num_copies = 1,   &   ! number of copies in sequence
                      num_qc     = 1        ! number of QC entries

integer :: ncid, nobs, nvars, n, i, oday, osec, nused

logical  :: file_exist, first_obs

real(r8) :: palt_miss, tair_miss, relh_miss, tdew_miss, wdir_miss, wspd_miss, uwnd, &
            vwnd, qobs, qsat, oerr, pres, qc, qerr

integer,  allocatable :: tobs(:), tused(:)
real(r8), allocatable :: lat(:), lon(:), palt(:), tair(:), relh(:), tdew(:), &
                         wdir(:), wspd(:)
integer,  allocatable :: qc_palt(:), qc_tair(:), qc_relh(:), qc_tdew(:), qc_wdir(:), qc_wspd(:)
integer,  allocatable :: used(:), sorted_used(:)

type(obs_sequence_type) :: obs_seq
type(obs_type)          :: obs, prev_obs
type(time_type)         :: comp_day0, time_obs, prev_time

call initialize_utilities('convert_madis_acars')

! put the reference date into DART format
call set_calendar_type(GREGORIAN)
comp_day0 = set_date(1970, 1, 1, 0, 0, 0)

first_obs = .true.

ncid = nc_open_file_readonly(acars_netcdf_file, 'convert_madis_acars')

call getdimlen(ncid, "recNum", nobs)
call set_missing_name("missing_value")

allocate( lat(nobs))  ;  allocate( lon(nobs))
allocate(palt(nobs))  ;  allocate(tobs(nobs))
allocate(tair(nobs))  ;  allocate(relh(nobs))
allocate(wdir(nobs))  ;  allocate(wspd(nobs))
allocate(tdew(nobs))  ;  allocate(tused(nobs))
allocate(used(nobs))  ;  allocate(sorted_used(nobs))

nvars = 3
if (include_specific_humidity) nvars = nvars + 1
if (include_relative_humidity) nvars = nvars + 1
if (include_dewpoint)          nvars = nvars + 1

allocate(qc_palt(nobs)) ;  allocate(qc_relh(nobs))
allocate(qc_tair(nobs)) ;  allocate(qc_tdew(nobs))
allocate(qc_wdir(nobs)) ;  allocate(qc_wspd(nobs))

! read in the data arrays
call getvar_real(ncid, "latitude",     lat            ) ! latitudes
call getvar_real(ncid, "longitude",    lon            ) ! longitudes
call getvar_real(ncid, "altitude",     palt, palt_miss) ! pressure altitude
call getvar_real(ncid, "temperature",  tair, tair_miss) ! air temperature
call getvar_real(ncid, "downlinkedRH", relh, relh_miss) ! relative humidity
call getvar_real(ncid, "dewpoint",     tdew, tdew_miss) ! dew-point temperature
call getvar_real(ncid, "windDir",      wdir, wdir_miss) ! wind direction
call getvar_real(ncid, "windSpeed",    wspd, wspd_miss) ! wind speed
call getvar_int (ncid, "timeObs",      tobs           ) ! observation time

! if user says to use them, read in QCs if present
if (use_input_qc) then
   call get_or_fill_QC(ncid, "altitudeQCR",     qc_palt)
   call get_or_fill_QC(ncid, "temperatureQCR",  qc_tair)
   call get_or_fill_QC(ncid, "downlinkedRHQCR", qc_relh)
   call get_or_fill_QC(ncid, "dewpointQCR",     qc_tdew)
   call get_or_fill_QC(ncid, "windDirQCR",      qc_wdir)
   call get_or_fill_QC(ncid, "windSpeedQCR",    qc_wspd)
else
   qc_palt = 0 ;  qc_tair = 0
   qc_relh = 0 ;  qc_tdew = 0
   qc_wdir = 0 ;  qc_wspd = 0
endif

call nc_close_file(ncid, 'convert_madis_acars')

!  either read existing obs_seq or create a new one
call static_init_obs_sequence()
call init_obs(obs,      num_copies, num_qc)
call init_obs(prev_obs, num_copies, num_qc)

inquire(file=acars_out_file, exist=file_exist)

if ( file_exist ) then

  ! existing file found, append to it
  call read_obs_seq(acars_out_file, 0, 0, nvars*nobs, obs_seq)

else

  ! create a new one
  call init_obs_sequence(obs_seq, num_copies, num_qc, nvars*nobs)
  do i = 1, num_copies
    call set_copy_meta_data(obs_seq, i, 'MADIS observation')
  end do
  do i = 1, num_qc
    call set_qc_meta_data(obs_seq, i, 'Data QC')
  end do

endif

! Set the DART data quality control.  Be consistent with NCEP codes;
! 0 is 'must use', 1 is good, no reason not to use it.
qc = 1.0_r8

nused = 0
obsloop1: do n = 1, nobs

  ! check the lat/lon values to see if they are ok
  if ( lat(n) >  90.0_r8 .or. lat(n) <  -90.0_r8 ) cycle obsloop1
  if ( lon(n) > 180.0_r8 .or. lon(n) < -180.0_r8 ) cycle obsloop1

  ! change lon from -180 to 180 into 0-360
  if ( lon(n) < 0.0_r8 )  lon(n) = lon(n) + 360.0_r8

  ! Check for duplicate observations
  do i = 1, nused
    if ( lon(n)  == lon(used(i)) .and. &
         lat(n)  == lat(used(i)) .and. &
         tobs(n) == tobs(used(i)) .and. &
         palt(n) == palt(used(i)) ) cycle obsloop1
  end do

  nused = nused + 1
  used(nused) = n
  tused(nused) = tobs(n)

end do obsloop1

! sort obs by time
call index_sort(tused, sorted_used, nused)

obsloop2: do i = 1, nused

  ! get the next unique observation in sorted time order
  n = used(sorted_used(i))

  ! compute time of observation
  time_obs = increment_time(comp_day0, tobs(n))

  if ( palt(n) == palt_miss .or. qc_palt(n) /= 0 ) cycle obsloop2
  pres = pres_alt_to_pres(palt(n))

  ! extract actual time of observation in file into oday, osec.
  call get_time(time_obs, osec, oday)

  ! add wind component data to obs_seq
  if ( wdir(n) /= wdir_miss .and. qc_wdir(n) == 0 .and. &
       wspd(n) /= wspd_miss .and. qc_wspd(n) == 0  ) then

      call wind_dirspd_to_uv(wdir(n), wspd(n), uwnd, vwnd)
      oerr = acars_wind_error(pres * 0.01_r8)
      if ( abs(uwnd) < 150.0_r8 .and. abs(vwnd) < 150.0_r8 .and. oerr /= missing_r8) then

         call create_3d_obs(lat(n), lon(n), pres, VERTISPRESSURE, uwnd, &
                            ACARS_U_WIND_COMPONENT, oerr, oday, osec, qc, obs)
         call add_obs_to_seq(obs_seq, obs, time_obs, prev_obs, prev_time, first_obs)

         call create_3d_obs(lat(n), lon(n), pres, VERTISPRESSURE, vwnd, &
                            ACARS_V_WIND_COMPONENT, oerr, oday, osec, qc, obs)
         call add_obs_to_seq(obs_seq, obs, time_obs, prev_obs, prev_time, first_obs)

      endif

  endif

  ! add air temperature data to obs_seq
  if ( tair(n) /= tair_miss .and. qc_tair(n) == 0 ) then

      oerr = acars_temp_error(pres * 0.01_r8)
      if ( tair(n) >= 180.0_r8 .and. tair(n) <= 330.0_r8 .and. oerr /= missing_r8 ) then

         call create_3d_obs(lat(n), lon(n), pres, VERTISPRESSURE, tair(n), &
                              ACARS_TEMPERATURE, oerr, oday, osec, qc, obs)
         call add_obs_to_seq(obs_seq, obs, time_obs, prev_obs, prev_time, first_obs)

      endif

  endif

  ! if the air, dewpoint, and relh obs are all ok, then see which of
  ! the three types of moisture obs to generate.
  if ( tair(n) /= tair_miss .and. qc_tair(n) == 0 .and. &
       tdew(n) /= tdew_miss .and. qc_tdew(n) == 0 .and. &
       relh(n) /= relh_miss .and. qc_relh(n) == 0  ) then

       ! before we start computing things based on the dewpoint,
       ! make sure it isn't larger than the air temp.  if it is
       ! more than a degree larger, skip it completely.  if it is
       ! less, set them equal and continue.
       if (tdew(n) > tair(n)) then
          if (tdew(n) > tair(n) + 1.0_r8) goto 100
          tdew(n) = tair(n)
       endif

      ! add specific humidity to obs_seq
      if ( include_specific_humidity ) then

         qsat = specific_humidity(sat_vapor_pressure(tair(n)), pres)
         qobs = qsat * relh(n)
         if ( LH_err ) then
            !GSR decided to get error from tdew and tair, while passing
            !the rh from the read ob
            qerr = rh_error_from_dewpt_and_temp(tair(n), tdew(n))
         else
            qerr = acars_rel_hum_error(pres * 0.01_r8, tair(n), relh(n))
         endif
         oerr = max(qerr * qsat, 0.0001_r8)

         if ( qobs >  0.0_r8  .and. &
              qobs <= 0.07_r8 .and. qerr /= missing_r8 ) then

            call create_3d_obs(lat(n), lon(n), pres, VERTISPRESSURE, qobs, &
                               ACARS_SPECIFIC_HUMIDITY, oerr, oday, osec, qc, obs)
            call add_obs_to_seq(obs_seq, obs, time_obs, prev_obs, prev_time, first_obs)

         endif

       endif

       ! add relative humidity data to obs_seq
       if ( include_relative_humidity ) then

         if ( LH_err ) then
            !GSR decided to get error from tdew and tair, while passing the rh
            !from the read ob
            oerr = rh_error_from_dewpt_and_temp(tair(n), tdew(n))
         else
            oerr = acars_rel_hum_error(pres * 0.01_r8, tair(n), relh(n))
         endif

         if ( relh(n) >  0.0_r8 .and. &
              relh(n) <= 1.5_r8 .and. oerr /= missing_r8 ) then

           call create_3d_obs(lat(n), lon(n), pres, VERTISPRESSURE, relh(n), &
                              ACARS_RELATIVE_HUMIDITY, oerr, oday, osec, qc, obs)
           call add_obs_to_seq(obs_seq, obs, time_obs, prev_obs, prev_time, first_obs)
         endif

      endif

      ! add dew-point temperature data to obs_seq
      if ( include_dewpoint ) then

         oerr = dewpt_error_from_rh_and_temp(tair(n), relh(n))

         if ( relh(n) >  0.0_r8 .and. &
              relh(n) <= 1.5_r8 .and. oerr /= missing_r8 ) then

           call create_3d_obs(lat(n), lon(n), pres, VERTISPRESSURE, tdew(n), &
                              ACARS_DEWPOINT, oerr, oday, osec, qc, obs)
           call add_obs_to_seq(obs_seq, obs, time_obs, prev_obs, prev_time, first_obs)

         endif
      endif

  endif     ! quality control/missing check on tair, tdew, and relh

100 continue

end do obsloop2

! if we added any obs to the sequence, write it now.
if ( get_num_obs(obs_seq) > 0 )  call write_obs_seq(obs_seq, acars_out_file)

! end of main program
call finalize_utilities()

end program convert_madis_acars

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

max_stars_repo_name: hkershaw-brown/feature-preprocess

max_issues_repo_name: hkershaw-brown/feature-preprocess

max_forks_repo_name: hkershaw-brown/feature-preprocess


text:
!!  Copyright (C)  Stichting Deltares, 2012-2016.
!!
!!  This program is free software: you can redistribute it and/or modify
!!  it under the terms of the GNU General Public License version 3,
!!  as published by the Free Software Foundation.
!!
!!  This program is distributed in the hope that it will be useful,
!!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!!  GNU General Public License for more details.
!!
!!  You should have received a copy of the GNU General Public License
!!  along with this program. If not, see <http://www.gnu.org/licenses/>.
!!
!!  contact: delft3d.support@deltares.nl
!!  Stichting Deltares
!!  P.O. Box 177
!!  2600 MH Delft, The Netherlands
!!
!!  All indications and logos of, and references to registered trademarks
!!  of Stichting Deltares remain the property of Stichting Deltares. All
!!  rights reserved.

      subroutine nutrel ( pmsa   , fl     , ipoint , increm , noseg  ,
     &                    noflux , iexpnt , iknmrk , noq1   , noq2   ,
     &                    noq3   , noq4   )
!>\file
!>       Release (nutrients/detritus) by of mortality algae DYNAMO

!
!     Description of the module :
!
! Name    T   L I/O   Description                                   Unit
! ----    --- -  -    -------------------                            ---
! FL( 1)  R*4 1 O autolysis of NN4                               [gN/m3/
! FL( 2)  R*4 1 O production of N-det                            [gN/m3/
! FL( 3)  R*4 1 O autolysis of P                                 [gP/m3/
! FL( 4)  R*4 1 O production of P-det                            [gP/m3/
! FL( 5)  R*4 1 O autolysis of Si                               [gSi/m3/
! FL( 6)  R*4 1 O production of Si-det                         [gSiC/m3/
! MRT1    R*4 1 I fraction of mortality dissolved as nutrients         [
! MRT2    R*4 1 I fraction of mortality dissolved as nutrients         [
! NCRAT1  R*4 1 I Nitrogen-Carbon ratio in green-algea             [gN/g
! NCRAT2  R*4 1 I Nitrogen-Carbon ratio in diatoms                 [gN/g
! PCRAT1  R*4 1 I Phosphorus-Carbon ratio in green-algea           [gP/g
! PCRAT2  R*4 1 I Phosphorus-Carbon ratio in diatoms               [gP/g
! RESP1   R*4 1 L total respiration rate const. green-algea          [1/
! RESP2   R*4 1 L total respiration rate const. diatoms              [1/
! SCRAT3  R*4 1 I Silicate-Carbon ratio in diatoms                [gSi/g

!     Logical Units : -

!     Modules called : -

!     Name     Type   Library
!     ------   -----  ------------

      IMPLICIT REAL (A-H,J-Z)

      REAL     PMSA  ( * ) , FL    (*)
      INTEGER  IPOINT( * ) , INCREM(*) , NOSEG , NOFLUX,
     +         IEXPNT(4,*) , IKNMRK(*) , NOQ1, NOQ2, NOQ3, NOQ4

      IP1  = IPOINT( 1)
      IP2  = IPOINT( 2)
      IP3  = IPOINT( 3)
      IP4  = IPOINT( 4)
      IP5  = IPOINT( 5)
      IP6  = IPOINT( 6)
      IP7  = IPOINT( 7)
      IP8  = IPOINT( 8)
      IP9  = IPOINT( 9)
      IP10 = IPOINT(10)
      IP11 = IPOINT(11)
!
      IFLUX = 0
      DO 9000 ISEG = 1 , NOSEG
!!    CALL DHKMRK(1,IKNMRK(ISEG),IKMRK1)
!!    IF (IKMRK1.GT.0) THEN
      IF (BTEST(IKNMRK(ISEG),0)) THEN
!

      MORT1     = PMSA(IP1 )
      NCRAT1    = PMSA(IP2 )
      PCRAT1    = PMSA(IP3 )
      FMRT1A    = PMSA(IP4 )
      FMRT2A    = PMSA(IP5 )
      MORT2     = PMSA(IP6 )
      NCRAT2    = PMSA(IP7 )
      PCRAT2    = PMSA(IP8 )
      SCRAT2    = PMSA(IP9 )
      FMRT1D    = PMSA(IP10)
      FMRT2D    = PMSA(IP11)

!***********************************************************************
!**** Processes connected to the ALGEA model
!***********************************************************************

!     Calculate fractions for carbon (different from nutrient fractions)
!     no part of carbon to autolyse!
       FDCA = 0.0
       FDCD = 0.0
       IF (FMRT1A .LT. 1.0) FDCA = FMRT2A / (1-FMRT1A)
       IF (FMRT1D .LT. 1.0) FDCD = FMRT2D / (1-FMRT1D)

!@    Production of DETC
      FL ( 1 + IFLUX ) = ( MORT1 * FDCA +
     &            MORT2 * FDCD )

!@    Production of OOC
      FL ( 2 + IFLUX ) = ( MORT1 * ( 1.0 - FDCA ) +
     &            MORT2 * ( 1.0 - FDCD ) )

!@    Autolysis of NH4
      FL ( 3 + IFLUX ) = ( MORT1 * NCRAT1 * FMRT1A +
     &            MORT2 * NCRAT2 * FMRT1D )

!@    Production of DETN
      FL ( 4 + IFLUX ) = ( MORT1 * NCRAT1 * FMRT2A +
     &            MORT2 * NCRAT2 * FMRT2D )

!@    Production of OON
      FL ( 5 + IFLUX ) = ( MORT1 * NCRAT1 * ( 1.0 - FMRT1A - FMRT2A ) +
     &            MORT2 * NCRAT2 * ( 1.0 - FMRT1D - FMRT2D ) )

!@    Autolysis of PO4
      FL ( 6 + IFLUX ) = ( MORT1 * PCRAT1 * FMRT1A +
     &            MORT2 * PCRAT2 * FMRT1D )

!@    Production of DETP
      FL ( 7 + IFLUX ) = ( MORT1 * PCRAT1 * FMRT2A +
     &            MORT2 * PCRAT2 * FMRT2D )

!@    Production of OOP
      FL ( 8 + IFLUX ) = ( MORT1 * PCRAT1 * ( 1.0 - FMRT1A - FMRT2A ) +
     &            MORT2 * PCRAT2 * ( 1.0 - FMRT1D - FMRT2D ) )

!@    Autolysis of Si
      FL ( 9 + IFLUX ) =   MORT2 * SCRAT2  * FMRT1D

!@    Production of Si-det
      FL (10 + IFLUX ) =   MORT2 * SCRAT2 * FMRT2D

!@    Production of OOSI
      FL (11 + IFLUX ) =   MORT2 * SCRAT2 * ( 1.0 - FMRT1D - FMRT2D )

      ENDIF
!
      IFLUX = IFLUX + NOFLUX
      IP1   = IP1   + INCREM (  1 )
      IP2   = IP2   + INCREM (  2 )
      IP3   = IP3   + INCREM (  3 )
      IP4   = IP4   + INCREM (  4 )
      IP5   = IP5   + INCREM (  5 )
      IP6   = IP6   + INCREM (  6 )
      IP7   = IP7   + INCREM (  7 )
      IP8   = IP8   + INCREM (  8 )
      IP9   = IP9   + INCREM (  9 )
      IP10  = IP10  + INCREM ( 10 )
      IP11  = IP11  + INCREM ( 11 )
!
 9000 CONTINUE
!
      RETURN
      END

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


text:
      SUBROUTINE INTEGRATE( NSENSIT, Y, TIN, TOUT )

      INCLUDE 'KPP_ROOT_params.h'
      INCLUDE 'KPP_ROOT_global.h'

      INTEGER NSENSIT
C TIN - Start Time
      KPP_REAL TIN
C TOUT - End Time
      KPP_REAL TOUT
C Y - Concentrations and Sensitivities
      KPP_REAL Y(NVAR*(NSENSIT+1))
C ---  Note: Y contains: (1:NVAR) concentrations, followed by
C ---                   (1:NVAR) sensitivities w.r.t. first parameter, followed by
C ---                   etc.,  followed by
C ---                   (1:NVAR) sensitivities w.r.t. NSENSIT's parameter

      INTEGER    INFO(5)

      EXTERNAL FUNC_CHEM, JAC_CHEM

      INFO(1) = Autonomous

      CALL ROS2_DDM(NVAR,NSENSIT,TIN,TOUT,STEPMIN,STEPMAX,
     +                   STEPMIN,Y,ATOL,RTOL,
     +                   Info,FUNC_CHEM,JAC_CHEM)

      RETURN
      END

      SUBROUTINE ROS2_DDM(N,NSENSIT,T,Tnext,Hmin,Hmax,Hstart,
     +                   y,AbsTol,RelTol,
     +                   Info,FUNC_CHEM,JAC_CHEM)
      IMPLICIT NONE
      INCLUDE 'KPP_ROOT_params.h'
      INCLUDE 'KPP_ROOT_global.h'
      INCLUDE 'KPP_ROOT_sparse.h'
C
C  Ros2 with direct-decoupled calculation of sensitivities
C
C The global variable DDMTYPE distinguishes between:
C      DDMTYPE = 0 : sensitivities w.r.t. initial values
C      DDMTYPE = 1 : sensitivities w.r.t. parameters
C
C  INPUT ARGUMENTS:
C     y = Vector of:   (1:NVAR) concentrations, followed by
C                      (1:NVAR) sensitivities w.r.t. first parameter, followed by
C                       etc.,  followed by
C                      (1:NVAR) sensitivities w.r.t. NSENSIT's parameter
C         (y contains initial values at input, final values at output)
C     [T, Tnext] = the integration interval
C     Hmin, Hmax = lower and upper bounds for the selected step-size.
C          Note that for Step = Hmin the current computed
C          solution is unconditionally accepted by the error
C          control mechanism.
C     AbsTol, RelTol = (NVAR) dimensional vectors of
C          componentwise absolute and relative tolerances.
C     FUNC_CHEM = name of routine of derivatives. KPP syntax.
C          See the header below.
C     JAC_CHEM = name of routine that computes the Jacobian, in
C          sparse format. KPP syntax. See the header below.
C     Info(1) = 1  for  autonomous   system
C             = 0  for nonautonomous system
C
C  OUTPUT ARGUMENTS:
C     y = the values of concentrations at TEND.
C     T = equals TEND on output.
C     Info(2) = # of FUNC_CHEM calls.
C     Info(3) = # of JAC_CHEM calls.
C     Info(4) = # of accepted steps.
C     Info(5) = # of rejected steps.
C
C  Adrian Sandu, December 2001

      INTEGER NSENSIT
      KPP_REAL y(NVAR*(NSENSIT+1)), ynew(NVAR*(NSENSIT+1))
      KPP_REAL K1(NVAR*(NSENSIT+1))
      KPP_REAL K2(NVAR*(NSENSIT+1))
      KPP_REAL K3(NVAR)
      KPP_REAL DFDT(NVAR*(NSENSIT+1))
      KPP_REAL DFDP(NVAR*NSENSIT+1), DFDPDT(NVAR*NSENSIT+1)
      KPP_REAL DJDP(NVAR*NSENSIT+1)
      KPP_REAL F1(NVAR), F2(NVAR)
      KPP_REAL JAC(LU_NONZERO), AJAC(LU_NONZERO)
      KPP_REAL DJDT(LU_NONZERO)
      KPP_REAL HESS(NHESS)
      KPP_REAL Hmin,Hmax,Hnew,Hstart,ghinv,uround
      KPP_REAL AbsTol(NVAR), RelTol(NVAR)
      KPP_REAL T, Tnext, H, Hold, Tplus, e
      KPP_REAL ERR, factor, facmax, dround, elo, tau, gam

      INTEGER    n,nfcn,njac,Naccept,Nreject,i,j,ier
      INTEGER    Info(5)
      LOGICAL    IsReject,Autonomous,Embed3
      EXTERNAL   FUNC_CHEM, JAC_CHEM, HESS_CHEM

      LOGICAL    negative
      KPP_REAL gamma, m1, m2, alpha, beta, delta, theta, w
      KPP_REAL gamma3, d1, d2, d3, beta1, beta2
      KPP_REAL c31, c32, c34

c     Initialization of counters, etc.
      Autonomous = Info(1) .EQ. 1
      Embed3  = Info(2) .EQ. 1
      uround = 1.d-15
      dround = 1.0d-7 ! DSQRT(uround)
      H = DMAX1(1.d-8, Hstart)
      Tplus = T
      IsReject = .false.
      Naccept  = 0
      Nreject  = 0
      Nfcn     = 0
      Njac     = 0
      gamma = 1.d0 + 1.d0/DSQRT(2.0d0)
      c31 = -1.0D0/gamma**2*(1.0D0-7.0D0*gamma+9.0D0*gamma**2)
     &         /(-1.0D0+2.0D0*gamma)
      c32 = -1.0D0/gamma**2*(1.0D0-6.0D0*gamma+6.0D0*gamma**2)
     &          /(-1.0D0+2.0D0*gamma)/2
      gamma3 = 0.5D0 - 2*gamma
      d1 = ((-9.0D0*gamma+8.0D0*gamma**2+2.0D0)/gamma**2/
     &          (-1.0D0+2*gamma))/6.0D0
      d2 = ((-1.0D0+3.0D0*gamma)/gamma**2/
     &          (-1.0D0+2.0D0*gamma))/6.0D0
      d3 = -1.0D0/(3.0D0*gamma)
       m1 = -3.d0/(2.d0*gamma)
       m2 = -1.d0/(2.d0*gamma)

C === Starting the time loop ===
 10    CONTINUE
       Tplus = T + H
       IF ( Tplus .gt. Tnext ) THEN
          H = Tnext - T
          Tplus = Tnext
       END IF

C              Initial Function, Jacobian, and Hessian Values
       CALL FUNC_CHEM(NVAR, T, y, F1)
       CALL JAC_CHEM(NVAR, T, y, JAC)
       CALL HESS_CHEM( NVAR, T, y, HESS )
       IF (DDMTYPE .EQ. 1) THEN
          CALL DFUNDPAR(NVAR, NSENSIT, T, y, DFDP)
       END IF

C              Estimate the time derivatives in non-autonomous case
       IF (.not. Autonomous) THEN
         tau = DSIGN(dround*DMAX1( 1.0d0, DABS(T) ), T)
         CALL FUNC_CHEM(NVAR, T+tau, y, K2)
         nfcn=nfcn+1
         CALL JAC_CHEM(NVAR, T+tau, y, AJAC)
         njac=njac+1
         DO 20 j = 1,NVAR
           DFDT(j) = ( K2(j)-F1(j) )/tau
 20      CONTINUE
         DO 30 j = 1,LU_NONZERO
           DJDT(j) = ( AJAC(j)-JAC(j) )/tau
 30      CONTINUE
         DO 40 i=1,NSENSIT
	    CALL Jac_SP_Vec (DJDT,y(i*NVAR+1),DFDT(i*NVAR+1))
 40      CONTINUE
        END IF ! .not. Autonomous

       Njac = Njac+1
       ghinv = - 1.0d0/(gamma*H)
       DO 50 j=1,LU_NONZERO
         AJAC(j) = JAC(j)
 50    CONTINUE
       DO 60 j=1,NVAR
         AJAC(LU_DIAG(j)) = JAC(LU_DIAG(j)) + ghinv
 60    CONTINUE
       CALL KppDecomp (AJAC, ier)

       IF (ier.ne.0) THEN
         IF ( H.gt.Hmin) THEN
            H = 5.0d-1*H
            go to 10
         ELSE
            print *,'IER <> 0, H=',H
            stop
         END IF
       END IF

C ----- STAGE 1 -----
         delta = gamma*H
         DO 70 j = 1,NVAR
           K1(j) =  F1(j)
 70      CONTINUE
         IF (.NOT. Autonomous) THEN
           DO 80 j = 1,NVAR
             K1(j) =  K1(j) + delta*DFDT(j)
 80        CONTINUE
	 END IF
         CALL KppSolve (AJAC, K1)
C               --- If  derivative w.r.t. parameters
	 IF (DDMTYPE .EQ. 1) THEN
	    CALL DJACDPAR(NVAR, NSENSIT, T, y, K1(1), DJDP)
	 END IF
C               --- End of derivative w.r.t. parameters
	 DO 120 i=1,NSENSIT
	    CALL Jac_SP_Vec (JAC,y(i*NVAR+1),K1(i*NVAR+1))
	    CALL Hess_Vec ( HESS, K1(1), y(i*NVAR+1), F2 )
	    DO 90 j=1,NVAR
	      K1(i*NVAR+j) = K1(i*NVAR+j) + gHinv*F2(j)
 90	    CONTINUE
           IF (.NOT. Autonomous) THEN
             DO 100 j = 1,NVAR
               K1(i*NVAR+j) =  K1(i*NVAR+j) + delta*DFDT(i*NVAR+j)
 100          CONTINUE
	   END IF
C               --- If  derivative w.r.t. parameters
	  IF (DDMTYPE .EQ. 1) THEN
            DO 110 j = 1,NVAR
	       K1(i*NVAR+j) = K1(i*NVAR+j) + DFDP((i-1)*NVAR+j)
     &                           + DJDP((i-1)*NVAR+j)
 110        CONTINUE
	  END IF
C               --- End of derivative w.r.t. parameters
            CALL KppSolve (AJAC, K1(i*NVAR+1))
 120	 CONTINUE

C ----- STAGE 2   -----
       alpha = - 1.d0/gamma
       DO 130 j = 1,NVAR*(NSENSIT+1)
         ynew(j) = y(j) + alpha*K1(j)
 130    CONTINUE
       CALL FUNC_CHEM(NVAR, T+H, ynew, F1)
       IF (DDMTYPE.EQ.1) THEN
         CALL DFUNDPAR(NVAR, NSENSIT, T+H, ynew, DFDP)
       END IF
       nfcn=nfcn+1
       beta1 = 2.d0/(gamma*H)
       delta = -gamma*H
       DO 140 j = 1,NVAR
         K2(j) = F1(j) + beta1*K1(j)
 140   CONTINUE
       IF (.NOT. Autonomous) THEN
          DO 150 j = 1,NVAR
               K2(j) =  K2(j) + delta*DFDT(j)
 150      CONTINUE
       END IF
       CALL KppSolve (AJAC, K2)
C               --- If  derivative w.r.t. parameters
       IF (DDMTYPE .EQ. 1) THEN
	    CALL DJACDPAR(NVAR, NSENSIT, T, y, K2(1), DJDP)
       END IF
C               --- End of derivative w.r.t. parameters

       CALL JAC_CHEM(NVAR, T+H, Ynew, JAC)
       njac=njac+1
       DO 190 i=1,NSENSIT
	  CALL Jac_SP_Vec (JAC,ynew(i*NVAR+1),K2(i*NVAR+1))
	  CALL Jac_SP_Vec (DJDT,y(i*NVAR+1),F1)
	  CALL Hess_Vec ( HESS, K2(1), y(i*NVAR+1), F2 )
          DO 160 j = 1,NVAR
	     K2(i*NVAR+j) = K2(i*NVAR+j) + beta1*K1(i*NVAR+j)
     &                       + gHinv*F2(j)
 160      CONTINUE
          IF (.NOT. Autonomous) THEN
             DO 170 j = 1,NVAR
               K2(i*NVAR+j) =  K2(i*NVAR+j) + delta*DFDT(i*NVAR+j)
 170         CONTINUE
	  END IF
C               --- If  derivative w.r.t. parameters
	  IF (DDMTYPE .EQ. 1) THEN
            DO 180 j = 1,NVAR
	       K2(i*NVAR+j) = K2(i*NVAR+j)  + DFDP((i-1)*NVAR+j)
     &                           + DJDP((i-1)*NVAR+j)
 180        CONTINUE
	  END IF
C               --- End of derivative w.r.t. parameters
          CALL KppSolve (AJAC, K2(i*NVAR+1))
 190    CONTINUE

C ----- STAGE 3  for error control only -----
       IF (Embed3) THEN
       beta1 = -c31/H
       beta2 = -c32/H
       delta = gamma3*H
       DO 195 j = 1,NVAR
         K3(j) = F1(j) + beta1*K1(j) + beta2*K2(j)
 195   CONTINUE
       IF (.NOT. Autonomous) THEN
         DO 196 j = 1,NVAR
           K3(j) = K3(j) + delta*DFDT(j)
 196     CONTINUE
       END IF
       CALL KppSolve (AJAC, K3)
       END IF

C ---- The Solution ---
       DO 200 j = 1,NVAR*(NSENSIT+1)
         ynew(j) = y(j) + m1*K1(j) + m2*K2(j)
 200   CONTINUE

C ====== Error estimation for concentrations only; this can be easily adapted to
C                                             estimate the sensitivity error too ========

        ERR=0.d0
        DO 210 i=1,NVAR
           w = AbsTol(i) + RelTol(i)*DMAX1(DABS(y(i)),DABS(ynew(i)))
	   IF (Embed3) THEN
	     e = d1*K1(i) + d2*K2(i) + d3*K3(i)
	   ELSE
             e = (1.d0/(2.d0*gamma))*(K1(i)+K2(i))
	   END IF
           ERR = ERR + ( e/w )**2
 210    CONTINUE
        ERR = DMAX1( uround, DSQRT( ERR/NVAR ) )

C ======= Choose the stepsize ===============================

	IF (Embed3) THEN
           elo    = 3.0D0 ! estimator local order
	ELSE
	   elo    = 2.0D0
	END IF
        factor = DMAX1(2.0D-1,DMIN1(6.0D0,ERR**(1.0D0/elo)/.9D0))
        Hnew   = DMIN1(Hmax,DMAX1(Hmin, H/factor))

C ======= Rejected/Accepted Step ============================

        IF ( (ERR.gt.1).and.(H.gt.Hmin) ) THEN
          IsReject = .true.
	  H = DMIN1(H/10,Hnew)
          Nreject  = Nreject+1
        ELSE
          DO 300 i=1,NVAR*(NSENSIT+1)
             y(i)  = ynew(i)
 300      CONTINUE
          T = Tplus
	  IF (.NOT.IsReject) THEN
	      H = Hnew   ! Do not increase stepsize if previous step was rejected
	  END IF
          IsReject = .false.
          Naccept = Naccept+1
        END IF

C ======= End of the time loop ===============================
      IF ( T .lt. Tnext ) GO TO 10

C ======= Output Information =================================
      Info(2) = Nfcn
      Info(3) = Njac
      Info(4) = Naccept
      Info(5) = Nreject

      RETURN
      END

      SUBROUTINE FUNC_CHEM(N, T, Y, P)
      INCLUDE 'KPP_ROOT_params.h'
      INCLUDE 'KPP_ROOT_global.h'
      KPP_REAL   T, Told
      KPP_REAL   Y(NVAR), P(NVAR)
      Told = TIME
      TIME = T
      CALL Update_SUN()
      CALL Update_RCONST()
      CALL Fun( Y,  FIX, RCONST, P )
      TIME = Told
      RETURN
      END

      SUBROUTINE DFUNDPAR(N, NSENSIT, T, Y, P)
C ---  Computes the partial derivatives of FUNC_CHEM w.r.t. parameters
      INCLUDE 'KPP_ROOT_params.h'
      INCLUDE 'KPP_ROOT_global.h'
C ---  NCOEFF, JCOEFF useful for derivatives w.r.t. rate coefficients
      INTEGER NCOEFF, JCOEFF(NREACT)
      COMMON /DDMRCOEFF/ NCOEFF, JCOEFF

      KPP_REAL   T, Told
      KPP_REAL   Y(NVAR), P(NVAR*NSENSIT)
      Told = TIME
      TIME = T
      CALL Update_SUN()
      CALL Update_RCONST()
C
      IF (DDMTYPE .EQ. 0) THEN
C ---  Note: the values below are for sensitivities w.r.t. initial values;
C ---       they may have to be changed for other applications
        DO j=1,NSENSIT
          DO i=1,NVAR
	    P(i+NVAR*(j-1)) = 0.0D0
	  END DO
        END DO
      ELSE
C ---  Example: the call below is for sensitivities w.r.t. rate coefficients;
C ---       JCOEFF(1:NSENSIT) are the indices of the NSENSIT rate coefficients
C ---       w.r.t. which one differentiates
        CALL dFun_dRcoeff( Y,  FIX, NCOEFF, JCOEFF, P )
      END IF
      TIME = Told
      RETURN
      END

      SUBROUTINE DJACDPAR(N, NSENSIT, T, Y, U, P)
C ---  Computes the partial derivatives of JAC w.r.t. parameters times user vector U
      INCLUDE 'KPP_ROOT_params.h'
      INCLUDE 'KPP_ROOT_global.h'
C ---  NCOEFF, JCOEFF useful for derivatives w.r.t. rate coefficients
      INTEGER NCOEFF, JCOEFF(NREACT)
      COMMON /DDMRCOEFF/ NCOEFF, JCOEFF

      KPP_REAL   T, Told
      KPP_REAL   Y(NVAR), U(NVAR)
      KPP_REAL   P(NVAR*NSENSIT)
      Told = TIME
      TIME = T
      CALL Update_SUN()
      CALL Update_RCONST()
C
      IF (DDMTYPE .EQ. 0) THEN
C ---  Note: the values below are for sensitivities w.r.t. initial values;
C ---       they may have to be changed for other applications
        DO j=1,NSENSIT
          DO i=1,NVAR
	    P(i+NVAR*(j-1)) = 0.0D0
	  END DO
        END DO
      ELSE
C ---  Example: the call below is for sensitivities w.r.t. rate coefficients;
C ---       JCOEFF(1:NSENSIT) are the indices of the NSENSIT rate coefficients
C ---       w.r.t. which one differentiates
        CALL dJac_dRcoeff( Y,  FIX, U, NCOEFF, JCOEFF, P )
      END IF
      TIME = Told
      RETURN
      END

      SUBROUTINE JAC_CHEM(N, T, Y, J)
      INCLUDE 'KPP_ROOT_params.h'
      INCLUDE 'KPP_ROOT_global.h'
      INTEGER N
      KPP_REAL   Told, T
      KPP_REAL   Y(NVAR), J(LU_NONZERO)
      Told = TIME
      TIME = T
      CALL Update_SUN()
      CALL Update_RCONST()
      CALL Jac_SP( Y, FIX, RCONST, J )
      TIME = Told
      RETURN
      END

      SUBROUTINE HESS_CHEM(N, T, Y, HESS)
      INCLUDE 'KPP_ROOT_params.h'
      INCLUDE 'KPP_ROOT_global.h'
      INTEGER N
      KPP_REAL   Told, T
      KPP_REAL   Y(NVAR), HESS(NHESS)
      Told = TIME
      TIME = T
      CALL Update_SUN()
      CALL Update_RCONST()
      CALL Hessian( Y, FIX, RCONST, HESS )
      TIME = Told
      RETURN
      END

max_stars_repo_name: matzegoebel/WRF-fluxavg

max_issues_repo_name: teb-model/wrf-teb

max_forks_repo_name: teb-model/wrf-teb


text:
      subroutine putc (c)
      integer c
      call putch (c, 1)
      return
      end

max_stars_repo_name: ysBach/irafdocgen

max_issues_repo_name: kirxkirx/iraf

max_forks_repo_name: kirxkirx/iraf