text:
c     Response function factor = 0.78896
      subroutine avhrr14(z0, p0, t0, u0, ux, nlev,
     *co2, n2o, f11, f12, ch4,
     *od, wght, ntau)
      INCLUDE "avhrr.inc"
      real co2, n2o, ch4, f11, f12
      real z0(mm)
      real od(3,mm)
      real wght(3)
      integer ntau,nlev
      real wi(8),xi(8),p0(mq),p(mm),t0(mq),t(mm),bnu(mm),u0(mq),u(mm),
     *ux(mq),uo3(mm),tau(6,mm),tauo3(6,mm),delw,x(mm,8),y(mm,8),
     *sx(6,mm),sy(6,mm),f(6),fac(mm),tsx(mq),tsy(mq),fo3(6)
      data xi /.01986,.10167,.23723,.40828,
     *         .59172,.76277,.89833,.98014/
      data wi /.05062,.11119,.15685,.18134,
     *         .18134,.15685,.11119,.05062/
      tb=0
      nlyr = nlev-1
      if(nlyr.GT.mm)then
        write (*,*) 'error: increase mm to at least', nlyr
        write (*,*) '       in avhrr.inc'
        return
      end if
      if(nlev.GT.mq)then
        write (*,*) 'error: increase mq to at least', nlev
        write (*,*) '       in avhrr.inc'
        return
      end if
      omega=16950.0
      delw=700.0
      sigma=5.67e-05
      sconst=1.368e+06
      tsun=5710.0
      fsun=sigma*tsun**4
      mlv=nlyr
      do 5 m=1,nlyr
        t(m)=(t0(m)+t0(m+1))/2.0
        p(m)=(p0(m)+p0(m+1))/2.0
        call flux_14(bnu(m),omega,t(m))
        fac(m)=z0(m+1)-z0(m)
        u(m)=fac(m)*(u0(m)+u0(m+1))/2.0
        uo3(m)=fac(m)*(ux(m)+ux(m+1))/2.0
 5    continue
      do 7 m=1,nlev
        tsx(m)=0.0
        tsy(m)=0.0
 7    continue
      call flux_14(back,omega,t0(1))
      call flux_14(bsun,omega,tsun)
      bsun=bsun/fsun*sconst
      call ck_14(ni,mlv,u,delw,f,p,t,tau)
      call cko3_14(no3,mlv,uo3,delw,fo3,p,t,tauo3)
c     **fo3=1.0000**
      do 10 m=1,nlyr
        do i=1,ni
          tau(i,m)=tau(i,m)+tauo3(1,m)
c          tau(i,m)=tauo3(1,m)
        end do
 10   continue
      do 90 i=1,ni
        do 30 m=nlyr,1,-1
          do 20 j=1,8
           if(m.eq.nlyr)then
c           y(m,j)=bnu(m)*(1.0-exp(-tau(i,m)/xi(j)))
            y(m,j)=bsun*0.6*exp(-tau(i,m)/xi(j))+
     *             bnu(m)*(1.0-exp(-tau(i,m)/xi(j)))
           else
            y(m,j)=y(m+1,j)*exp(-tau(i,m)/xi(j))+
     *      bnu(m)*(1.0-exp(-tau(i,m)/xi(j)))
           end if
 20       continue
 30     continue
        do 50 m=1,nlyr
          sx(i,m)=0.0
          sy(i,m)=0.0
          do 40 j=1,8
           if(m.eq.1)then
            x(m,j)=(0.10*y(m,j)+0.90*back)*exp(-tau(i,m)/xi(j))+
     *      bnu(m)*(1.0-exp(-tau(i,m)/xi(j)))
           else
            x(m,j)=x(m-1,j)*exp(-tau(i,m)/xi(j))+
     *      bnu(m)*(1.0-exp(-tau(i,m)/xi(j)))
           end if
 40       continue
 50     continue
        do 70 m=1,nlyr
          do 60 j=1,8
            sx(i,m)=sx(i,m)+2.0*x(m,j)*wi(j)*xi(j)
            sy(i,m)=sy(i,m)+2.0*y(m,j)*wi(j)*xi(j)
 60       continue
 70     continue
        do 80 m=1,nlyr
           tsx(m+1)=tsx(m+1)+sx(i,m)*delw*f(i)
           tsy(m)=tsy(m)+sy(i,m)*delw*f(i)
 80     continue
c        tb=tb+(0.10*tsy(1))+0.90*back*delw*f(i)
        tb=tb+back*delw*f(i)
c        tsx(1)=tsx(1)+0.10*tsy(1)
c        tsy(nlev)=tsy(nlev)+bsun*0.60*delw*f(i)
 90   continue
      tb=0.10*tsy(1)+0.90*tb
      tsx(1)=tb
      tsy(nlev)=tsy(nlev)+bsun*0.60*delw
      do 95 m=1,nlev
cbm         write(6,100)m-1,tsx(m),tb-tsx(m),tsy(m)
 95   continue
      do i1=1,3
       do m=1,nlyr
        od(1+(i1-1),m) =
     *tau(i1,m)
        wght(1+(i1-1)) =
     *f(i1)
       enddo
      enddo
      ntau=3
 100  format(i3,3f15.7)
      end
c *********************************************************************
      subroutine flux_14(bnu,omega,t)
      INCLUDE "avhrr.inc"
      bnu=0.0
      z=3.14159
      bnu=z/1000.0*((2.*6.625e-27)*(2.997925e+10)**2
     *    *omega**3)/(exp(min(88.0,1.4388*omega/t))-1.0)
      return
      end
c *********************************************************************
      subroutine ck_14(ni,mlv,u,delw,f,p,t,tau)
      INCLUDE "avhrr.inc"
      real f(6),k(mm),tau(6,mm),delw,p(mm),t(mm),u(mm),coefk(6,3,19),
     *     stp(19),fkg(6,mm),pmb(mm)
      ni=3
      f(1)=0.914321
      f(2)=0.054775
      f(3)=0.030904
      k(1)=0.0024289518
      do i=2,ni
        k(i)=8.0*k(i-1)
      end do
      data ( (coefk(1,jt,jp), jp = 1, 19), jt = 1, 3)/
     *0.1013,0.1013,0.1013,0.1013,0.1013,0.1013,0.1013,0.1013,0.1013,
     *0.1235,0.1569,0.2061,0.2763,0.3743,0.5061,0.6784,0.8855,1.1379,
     *1.4178,
     *5.336e-04,5.336e-04,5.336e-04,5.336e-04,5.336e-04,5.336e-04,
     *5.336e-04,5.336e-04,5.336e-04,4.715e-04,3.877e-04,2.794e-04,
     *1.297e-04,-6.213e-05,-2.831e-04,-5.150e-04,-6.824e-04,-8.774e-04,
     *-1.020e-03,
     *-7.031e-07,-7.031e-07,-7.031e-07,-7.031e-07,-7.031e-07,-7.031e-07,
     *-7.031e-07,-7.031e-07,-7.031e-07,-8.500e-07,-1.156e-06,-1.422e-06,
     *-1.387e-06,-1.353e-06,-9.406e-07,-7.187e-07,3.781e-07,8.406e-07,
     *-1.594e-07/
      data ( (coefk(2,jt,jp), jp = 1, 19), jt = 1, 3)/
     *0.5929,0.5929,0.5929,0.5929,0.5929,0.5929,0.5929,0.5929,0.5929,
     *0.6013,0.6154,0.6377,0.6729,0.7252,0.7952,0.8817,0.9649,1.0435,
     *1.0898,
     *1.616e-03,1.616e-03,1.616e-03,1.616e-03,1.616e-03,1.616e-03,
     *1.616e-03,1.616e-03,1.616e-03,1.529e-03,1.392e-03,1.196e-03,
     *9.329e-04,6.607e-04,3.900e-04,1.035e-04,-9.275e-05,-2.953e-04,
     *-4.020e-04,
     *-1.747e-06,-1.747e-06,-1.747e-06,-1.747e-06,-1.747e-06,-1.747e-06,
     *-1.747e-06,-1.747e-06,-1.747e-06,-7.375e-07,1.281e-07,1.787e-06,
     *3.016e-06,3.037e-06,3.469e-06,2.931e-06,1.950e-06,7.625e-07,
     *1.252e-08/
      data ( (coefk(3,jt,jp), jp = 1, 19), jt = 1, 3)/
     *1.5250,1.5250,1.5250,1.5250,1.5250,1.5250,1.5250,1.5250,1.5250,
     *1.5128,1.4948,1.4677,1.4271,1.3719,1.2951,1.1977,1.0700,0.9381,
     *0.7974,
     *-1.931e-03,-1.931e-03,-1.931e-03,-1.931e-03,-1.931e-03,-1.931e-03,
     *-1.931e-03,-1.931e-03,-1.931e-03,-1.891e-03,-1.820e-03,-1.716e-03,
     *-1.613e-03,-1.465e-03,-1.284e-03,-1.130e-03,-9.167e-04,-8.091e-04,
     *-7.313e-04,
     *3.766e-06,3.766e-06,3.766e-06,3.766e-06,3.766e-06,3.766e-06,
     *3.766e-06,3.766e-06,3.766e-06,4.028e-06,3.250e-06,3.303e-06,
     *3.447e-06,2.297e-06,2.563e-06,1.969e-06,2.044e-06,2.366e-06,
     *2.619e-06/
	data stp / 0.251, 0.398, 0.631, 1.000, 1.58, 2.51,
     *	             3.98, 6.31, 10.0, 15.8, 25.1, 39.8, 63.1,
     *	             100.0, 158.0, 251.0, 398.0, 631.0, 1000.0/
      do i=1,ni
      do m=1,mlv
        ml=1
        pmb(m)=p(m)*1013.25
        if(pmb(m).lt.stp(1))then
         x1=coefk(i,1,1)+coefk(i,2,1)*(t(m)-250.0)
     *   +coefk(i,3,1)*(t(m)-250.0)**2
         fkg(i,m)=x1*pmb(m)/stp(1)
        else if (pmb(m).gt.stp(19)) then
         x1=coefk(i,1,18)+coefk(i,2,18)*(t(m)-250.0)
     *   +coefk(i,3,18)*(t(m)-250.0)**2
         x2=coefk(i,1,19)+coefk(i,2,19)*(t(m)-250.0)
     *   +coefk(i,3,19)*(t(m)-250.0)**2
         fkg(i,m)=x1+(x2-x1)/(stp(19)-stp(18))
     *  *(pmb(m)-stp(18))
        else
         do while(pmb(m).ge.stp(ml))
           ml=ml+1
         end do
         x1=coefk(i,1,ml-1)+coefk(i,2,ml-1)*(t(m)-250.0)
     *   + coefk(i,3,ml-1)*(t(m)-250.0)**2
         x2=coefk(i,1,ml)+coefk(i,2,ml)*(t(m)-250.0)
     *   + coefk(i,3,ml)*(t(m)-250.0)**2
         fkg(i,m)=x1+(x2-x1)/(stp(ml)-stp(ml-1))
     *   *(pmb(m)-stp(ml-1))
        end if
      end do
      end do
      do i=1,ni
        do m=1,mlv
          tau(i,m)=k(i)*u(m)*fkg(i,m)
c         write(6,*)m,i,u(m),k(i),fkg(i,m),tau(i,m)
        end do
      end do
      return
      end
c *********************************************************************
      subroutine cko3_14(ni,mlv,u,delw,f,p,t,tau)
      INCLUDE "avhrr.inc"
      real f(6),k(mm),tau(6,mm),delw,p(mm),t(mm),u(mm)
      ni=1
      f(1)=1.000000
      k(1)=4.55e-21*6.023e+23/47.9982
      do i=1,ni
        do m=1,mlv
          tau(i,m)=k(i)*u(m)
        end do
      end do
      return
      end


text:
  subroutine gamma (evapObj, compound, uran, beta, calcVars)

! ======================================================================
!
! GAMMA
!  Calculate decay width for each particle emission channel.
!
!  Called by: STDCAY
!
!  Calls:
! ====================================================================
! <Subroutine>
!   eye    : Calculate I0,I1,I2,I3
! ====================================================================
! <Function>
!  dost    : Calculate kp, cp, or k_alpha
!  paire   : Calculate pairing energy (NOT USED)
! ====================================================================
! <variables>
!     a    :   residual mass before emission                     (IN)
!     z    :   charge number of res nucleus before emission      (IN)
!     u    :   excitation energy of nucleus before emission      (IN)
!   uran   :   random number                                     (IN)
!     r    :   decay width                                       (OUT)
!     rr   :   decay width enchancement factor by excited-state
!                            particle emission                   (OUT)
!   ifa    :   mass # of emitted particle                        (OUT)
!   ifz    :   charge # of emitted particle                      (OUT)
!  omega   :   spin of emittor                                   (OUT)
!   gj     :   gj in eq.(39)  (~ spin multiplicity factor*mu)    (OUT)
!   q      :   Q-value                                           (OUT)
!  delta   :   pairing energy                                    (OUT)
!  smalla  :   level density parameter                           (OUT)
!   v      :   Coulomb barrier                                   (OUT)
!   exm    :   mass of an excited state [MeV]                    (OUT)
!  spin    :   spin of an excited state                          (OUT)
!  width   :   lifetime of an excited state [MeV]                (OUT)
!   gamn   :   Decay width for neutron emission                  (OUT)
!   an     :   Level density parameter for neutron emission      (OUT)
! beta,alp :   Inverse cross section params for neutron emission (OUT)
!   s      : NOT USED!  Removed 1/24/05 AJS.
!
!    "Last" change: 13-AUG-2003 by NVMokhov
!    Edited by A. J. Sierk, LANL T-16, September, 2003.
!    Modified by A. J. Sierk, LANL T-16, January, 2005.
!    Modified by A. J. Sierk, LANL T-16, March, 2005.
!    Modified by A. J. Sierk, LANL T-2, February, 2009.
!    Edited by AJS, LANL T-2, December, 2011.
!    Edited by LMK, XCP-3, July 2013 (included error protection).
!    Edited by CMJ, XCP-3, July 2018 (Evap class creation)
!
! ======================================================================

    use, intrinsic :: iso_fortran_env, only: int32, real64
    use evaporationParams, only: zro, hlf, one, two, thr, hbarc, amu, pi
    use evaporationFissionData, only: exm, spin, width, omega, ifa, &
         & ifz, pz, pn, paire0, maxSize, aMax, omfct, alpp, betap

    implicit none
    class(Evaporation),  intent(inout) :: evapObj
    type(evapCompound),  intent(in   ) :: compound
    real(real64),        intent(in   ) :: uran
    real(real64),        intent(  out) :: beta
    type(evaporationCalculation), intent(inout) :: calcVars

    integer(int32) :: i, ia, iaa, &
         & isdum, iz, izz, j, k, nn, nn0
    real(real64)   :: aa, alp, aq, bett, del, emred, gq, q1, q2, qq, &
         & rhop, rmass, rq, rrq, rrqg, temp, zz

! ======================================================================

    real(real64), parameter :: cst0 = amu / ( pi * (hbarc**2) )

! ======================================================================

!   Parent nucleus:
    ia = nint(compound%numBaryons)
    iz = nint(compound%numProtons)
    q1 = evapObj%energy (iz, ia)
    if (q1 >= 1.d10) then
       write(evapObj%io%message, 2000) ia, iz
       call evapObj%io%print(2, 2, evapObj%io%message)
       q1 = 1.d10
    endif

    if (evapObj%options%inverseParameter == zro) then
!   Set cj and kj for Dostrovsky's parameter set
       calcVars%couk(2) = dost (1, compound%numProtons-dble(ifz(2)))
       calcVars%couk(3) = calcVars%couk(2) + 6.d-2
       calcVars%couk(4) = calcVars%couk(2) + 1.2d-1
       calcVars%couk(6) = dost (2, compound%numProtons-dble(ifz(6)))
       calcVars%couk(5) = calcVars%couk(6) - 6.d-2
       calcVars%couc(2) = dost (3, compound%numProtons-dble(ifz(2)))
       calcVars%couc(3) = calcVars%couc(2)*hlf
       calcVars%couc(4) = calcVars%couc(2)/thr
!   ck is given by expression: ck = C1*log(Z) + C2*log(A) + C3
    elseif (evapObj%options%inverseParameter < zro) then
       do j = 2,6
          calcVars%couk(j) = ckcal (j, compound%numProtons - dble(ifz(j)), compound%numBaryons - dble(ifa(j)))
       end do
    endif
    nn0 = ia - iz
    calcVars%smalla0 = geta (evapObj%data%alev(), compound%kinEnergy, &
         & iz, nn0, isdum)
!  Pairing energy for the  parent nucleus:
    del = paire0(iz,ia-iz)
!  Rho_i (level density of parent))
    rhop = evapObj%rho (ia, iz, compound%kinEnergy, del, calcVars)
    if (rhop <= zro) then
       write(evapObj%io%message, 2100) rhop
       call evapObj%io%print(2, 3, evapObj%io%message)
    end if

!     Start calculating Gamma for each particle type:
    do j = 1, evapObj%options%numEvapType

!  Initialization:
       calcVars%r(j) = zro
       calcVars%rr(j) = one
       calcVars%gj(j) = zro
       if (ifa(j) == 0) go to 20
! ----------------------------------------------------------------------
!    Reduce calculation time:
! ----------------------------------------------------------------------
       if (evapObj%options%redCompTime.ne.0) then
          if (j > 6) then
             if (ia > 40 .and. uran < 0.95d0) go to 20
             if (ia > 30 .and. ia <= 40 .and. uran < 0.93d0) go to 20
             if (ia > 20 .and. ia <= 30 .and. uran < 0.7d0) go to 20
          endif
       endif
!  Daughter nucleus mass and charge:
       iaa = ia - ifa(j)
       aa = dble(iaa)
       izz = iz - ifz(j)
       zz = dble(izz)
       nn = iaa - izz
! ----------------------------------------------------------------------
!    Check the residual nuclei after the emission for reasonable Z, A
!    and avoid double counting:
! ----------------------------------------------------------------------
       if (iaa <= 0 .or. izz <= 0 .or. iaa < izz .or. nn <= 0) &
            & go to 20
       if (iaa < ifa(j) .or. izz < ifz(j)) then
          do k = 1,j-1
             if (iaa == ifa(k) .and. izz == ifz(k)) go to 20
          end do
       endif

!   Q-value
       q2 = evapObj%energy (izz, iaa)
       if (q2 >= 1.d10) then
          write(evapObj%io%message,2000) iaa, izz
          call evapObj%io%print(2, 2, evapObj%io%message)
          q2 = 1.d10
       endif
       calcVars%q(j) = q2 - q1 + evapObj%energy (ifz(j), ifa(j))

!   Pairing energy of daughter:
       calcVars%delta(j) = paire0(izz, nn)

!   ............................................................
!  Coulomb potential
       calcVars%v(j) = evapObj%vcoul (zz, aa, ifz(j), ifa(j), calcVars%couk(j), j)
!    ..Set Coulomb potential to zero if Q-value < 0, for light residual
!    ..... e.g. 8Be, 9B
       if (calcVars%q(j) <= zro .and. iaa <= 20) calcVars%v(j) = zro

!   Check whether the emission is energetically possible or not:
       if (compound%kinEnergy - calcVars%q(j) - calcVars%v(j) <= zro) go to 20

!   alpha and beta for neutrons: the precise parameter set:
       if (j == 1) then
          if (evapObj%options%inverseParameter <= zro) then
             alp = alpp(iaa)
             beta = betap(iaa)
          else
!   Simple model:
             alp = one
             beta = zro
          endif
          bett = beta
       else
          bett = -calcVars%v(j)
       endif

       if (evapObj%options%inverseParameter == -one .and. j > 6) then
!    Exact calculation for subthreshold reaction
          if (calcVars%r(1) + calcVars%r(2) + calcVars%r(3) + calcVars%r(4) + &
               & calcVars%r(5) + calcVars%r(6).ne.zro) call evapObj%calr &
               & (aa, zz, compound%kinEnergy, calcVars%q(j), calcVars%v(j), calcVars%delta(j), &
               & calcVars%r(j), calcVars)
       else
!   Level density parameter of daughter:
          calcVars%smalla(j) = geta ( evapObj%data%alev(), &
               & (compound%kinEnergy-calcVars%q(j)-calcVars%v(j)), &
               & izz, nn, isdum)
          call evapObj%eye (iaa, izz, nn, compound%kinEnergy, calcVars%q(j), calcVars%v(j), &
               & calcVars%delta(j), calcVars%smalla(j), bett, calcVars%r(j), calcVars)
       endif
!   ............................................................

!   (2Sj + 1)*muj*alpha
       temp = compound%numBaryons
       if (temp < div0Lim .and. temp > -div0Lim) then
          temp = div0Lim
          write(evapObj%io%message,1000) '308'
          call evapObj%io%print(4, 3, evapObj%io%message)
       end if
       emred = dble(ifa(j))*aa/temp
       calcVars%gj(j) = omfct(j)*emred
       if (j == 1) then
          calcVars%gj(j) = calcVars%gj(j)*alp
       else
          calcVars%gj(j) = calcVars%gj(j)*(one + calcVars%couc(j))
       endif

!  Geometric cross section... sigma R
       rmass = evapObj%rb (aa, ifa(j), zz, ifz(j), j)
       calcVars%gj(j) = rmass*rmass*calcVars%gj(j)
       if (iz==7 .and. ia==13 .and. j==2) then
!          print *, 'Geom. radius: ', rmass
!          print *, 'cp = ', calcVars%couc(2)
!          print *, 'Fragment energy: ', q(2), 'GEM2 cross section: ', gj(2)
       end if
!  Decay width:
       calcVars%r(j) = calcVars%gj(j)*calcVars%r(j)
       calcVars%r(j) = max(calcVars%r(j), zro)

!  Decay width calculation from excited states:
       if (calcVars%r(j) == zro .or. j <= 6) go to 20

!   Move outside the loop over j; AJS (02/11/09)
!  Pairing energy for the  parent nucleus:
!  Rho_i (level density of parent))
!   End of j-independent code fragment.  (02/11/09)
       rrq = calcVars%r(j)
       do i = 1,200
          if (width(j-6,i) == zro) go to 10
!   Q-value for an excited state:
          qq = calcVars%q(j) + exm(j-6,i)
!   Level density for an excited state:
          aq = geta (evapObj%data%alev(), compound%kinEnergy - qq - calcVars%v(j), izz, nn, isdum)
!   Integral part of gamma for an excited state:
          call evapObj%eye (iaa, izz, nn, compound%kinEnergy, qq, calcVars%v(j), &
               & calcVars%delta(j), aq, bett, rq, calcVars)
!   (2Sj + 1)*mu_j*alpha*sigma_R for an excited state:
          gq = (two*spin(j-6,i) + one)*emred
          gq = rmass*rmass*gq
!   gamma for an excited state [MeV]:
          temp = rhop
          if (temp < div0Lim .and. temp > -div0Lim) then
             temp = div0Lim
             write(evapObj%io%message,1000) '354'
             call evapObj%io%print(4, 3, evapObj%io%message)
          end if
          rrqg = gq*rq*cst0/temp
!   Reject an excited state if the decay width [MeV] < level width [MeV]:
          if (width(j-6,i) >= rrqg) go to 10
          rrq = rrq + gq*rq
       end do
10     temp = rrq
       if (temp < 1.0d-90 .and. temp > -1.0d-90) then
          print *, temp
          temp = 1.0d-90
          write(evapObj%io%message,1000) '370'
          call evapObj%io%print(4, 3, evapObj%io%message)
       end if
       calcVars%rr(j) = calcVars%r(j)/temp
       calcVars%r(j) = rrq
20     continue
    end do
!  End of j (ejectile index) DO loop ^
    calcVars%sigma = zro
    do j = 1, evapObj%options%numEvapType
       calcVars%sigma = calcVars%sigma + calcVars%r(j)
    end do

    return

! ======================================================================
1000 format("Divide by zero error prevented in 'gamma.f90', line ", A)
2000 format("Mass calculation error occurred in function ", &
            & "'energy' for the nucleus with ia = ", i3, &
            & " and iz = ", i3, ".")
2100 format("Negative level density parameter (", f8.2, &
          & ") ocurred in function 'rho'.")
! ======================================================================
  end subroutine gamma

max_stars_repo_name: lanl/generalized-spallation-model

max_stars_repo_licenses: Intel Unlicense

max_issues_repo_name: lanl/generalized-spallation-model

max_issues_repo_licenses: Intel Unlicense

max_forks_repo_name: lanl/generalized-spallation-model

max_forks_repo_licenses: Intel Unlicense


! { dg-do run }
! PR 49479 - this used not to print anything.
! Test case by Joost VandeVondele.
MODULE M1
  IMPLICIT NONE
  type foo
     character(len=5) :: x
  end type foo
CONTAINS
  SUBROUTINE S1(data)
    INTEGER, DIMENSION(:), INTENT(IN), &
         OPTIONAL                               :: DATA
    character(20) :: line
    IF (.not. PRESENT(data)) call abort
    write (unit=line,fmt='(I5)') size(data)
    if (line /= '    0               ') call abort
  END SUBROUTINE S1

  subroutine s_type(data)
    type(foo), dimension(:), intent(in), optional :: data
    character(20) :: line
    IF (.not. PRESENT(data)) call abort
    write (unit=line,fmt='(I5)') size(data)
    if (line /= '    0               ') call abort
  end subroutine s_type

  SUBROUTINE S2(N)
    INTEGER :: N
    INTEGER, ALLOCATABLE, DIMENSION(:, :)    :: blki
    type(foo), allocatable, dimension(:, :)  :: bar
    ALLOCATE(blki(3,N))
    allocate (bar(3,n))
    blki=0
    CALL S1(RESHAPE(blki,(/3*N/)))
    call s_type(reshape(bar, (/3*N/)))
  END SUBROUTINE S2

END MODULE M1

USE M1
CALL S2(0)
END


text:
      PROGRAM xks2d2s
C     driver for routine ks2d2s
      INTEGER NMAX
      PARAMETER(NMAX=1000)
      INTEGER idum,j,jtrial,n1,n2,ntrial
      REAL d,gasdev,prob,shrink,u,v,x1(NMAX),y1(NMAX),x2(NMAX),y2(NMAX)
1     write(*,*) 'INPUT N1,N2'
      read(*,*,END=999) n1,n2
      if (n1.gt.NMAX) then
        write(*,*) 'n1 too large.'
        goto 1
      endif
      if (n2.gt.NMAX) then
        write(*,*) 'n2 too large.'
        goto 1
      endif
      write(*,*) 'WHAT SHRINKAGE?'
      read(*,*,END=999) shrink
      write(*,*) 'HOW MANY TRIALS?'
      read(*,*,END=999) ntrial
      if (ntrial.gt.NMAX) then
        write(*,*) 'Too many trials.'
        goto 1
      endif
      idum=-287-ntrial-n1-n2
      do 13 jtrial=1,ntrial
        do 11 j=1,n1
          u=gasdev(idum)
          v=gasdev(idum)*shrink
          x1(j)=u+v
          y1(j)=u-v
11      continue
        do 12 j=1,n2
          u=gasdev(idum)*shrink
          v=gasdev(idum)
          x2(j)=u+v
          y2(j)=u-v
12      continue
        call ks2d2s(x1,y1,n1,x2,y2,n2,d,prob)
        write(*,'(1x,a7,2f12.6)') 'D,PROB= ',d,prob
13    continue
      goto 1
999   write(*,*) 'NORMAL COMPLETION'
      STOP
      END

ext: for

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


!*==slarrf.f90  processed by SPAG 7.51RB at 20:08 on  3 Mar 2022
!> \brief \b SLARRF finds a new relatively robust representation such that at least one of the eigenvalues is relatively isolated.
!
!  =========== DOCUMENTATION ===========
!
! Online html documentation available at
!            http://www.netlib.org/lapack/explore-html/
!
!> \htmlonly
!> Download SLARRF + dependencies
!>
!> [TGZ]
!>
!> [ZIP]
!>
!> [TXT]
!> \endhtmlonly
!
!  Definition:
!  ===========
!
!       SUBROUTINE SLARRF( N, D, L, LD, CLSTRT, CLEND,
!                          W, WGAP, WERR,
!                          SPDIAM, CLGAPL, CLGAPR, PIVMIN, SIGMA,
!                          DPLUS, LPLUS, WORK, INFO )
!
!       .. Scalar Arguments ..
!       INTEGER            CLSTRT, CLEND, INFO, N
!       REAL               CLGAPL, CLGAPR, PIVMIN, SIGMA, SPDIAM
!       ..
!       .. Array Arguments ..
!       REAL               D( * ), DPLUS( * ), L( * ), LD( * ),
!      $          LPLUS( * ), W( * ), WGAP( * ), WERR( * ), WORK( * )
!       ..
!
!
!> \par Purpose:
!  =============
!>
!> \verbatim
!>
!> Given the initial representation L D L^T and its cluster of close
!> eigenvalues (in a relative measure), W( CLSTRT ), W( CLSTRT+1 ), ...
!> W( CLEND ), SLARRF finds a new relatively robust representation
!> L D L^T - SIGMA I = L(+) D(+) L(+)^T such that at least one of the
!> eigenvalues of L(+) D(+) L(+)^T is relatively isolated.
!> \endverbatim
!
!  Arguments:
!  ==========
!
!> \param[in] N
!> \verbatim
!>          N is INTEGER
!>          The order of the matrix (subblock, if the matrix split).
!> \endverbatim
!>
!> \param[in] D
!> \verbatim
!>          D is REAL array, dimension (N)
!>          The N diagonal elements of the diagonal matrix D.
!> \endverbatim
!>
!> \param[in] L
!> \verbatim
!>          L is REAL array, dimension (N-1)
!>          The (N-1) subdiagonal elements of the unit bidiagonal
!>          matrix L.
!> \endverbatim
!>
!> \param[in] LD
!> \verbatim
!>          LD is REAL array, dimension (N-1)
!>          The (N-1) elements L(i)*D(i).
!> \endverbatim
!>
!> \param[in] CLSTRT
!> \verbatim
!>          CLSTRT is INTEGER
!>          The index of the first eigenvalue in the cluster.
!> \endverbatim
!>
!> \param[in] CLEND
!> \verbatim
!>          CLEND is INTEGER
!>          The index of the last eigenvalue in the cluster.
!> \endverbatim
!>
!> \param[in] W
!> \verbatim
!>          W is REAL array, dimension
!>          dimension is >=  (CLEND-CLSTRT+1)
!>          The eigenvalue APPROXIMATIONS of L D L^T in ascending order.
!>          W( CLSTRT ) through W( CLEND ) form the cluster of relatively
!>          close eigenalues.
!> \endverbatim
!>
!> \param[in,out] WGAP
!> \verbatim
!>          WGAP is REAL array, dimension
!>          dimension is >=  (CLEND-CLSTRT+1)
!>          The separation from the right neighbor eigenvalue in W.
!> \endverbatim
!>
!> \param[in] WERR
!> \verbatim
!>          WERR is REAL array, dimension
!>          dimension is >=  (CLEND-CLSTRT+1)
!>          WERR contain the semiwidth of the uncertainty
!>          interval of the corresponding eigenvalue APPROXIMATION in W
!> \endverbatim
!>
!> \param[in] SPDIAM
!> \verbatim
!>          SPDIAM is REAL
!>          estimate of the spectral diameter obtained from the
!>          Gerschgorin intervals
!> \endverbatim
!>
!> \param[in] CLGAPL
!> \verbatim
!>          CLGAPL is REAL
!> \endverbatim
!>
!> \param[in] CLGAPR
!> \verbatim
!>          CLGAPR is REAL
!>          absolute gap on each end of the cluster.
!>          Set by the calling routine to protect against shifts too close
!>          to eigenvalues outside the cluster.
!> \endverbatim
!>
!> \param[in] PIVMIN
!> \verbatim
!>          PIVMIN is REAL
!>          The minimum pivot allowed in the Sturm sequence.
!> \endverbatim
!>
!> \param[out] SIGMA
!> \verbatim
!>          SIGMA is REAL
!>          The shift used to form L(+) D(+) L(+)^T.
!> \endverbatim
!>
!> \param[out] DPLUS
!> \verbatim
!>          DPLUS is REAL array, dimension (N)
!>          The N diagonal elements of the diagonal matrix D(+).
!> \endverbatim
!>
!> \param[out] LPLUS
!> \verbatim
!>          LPLUS is REAL array, dimension (N-1)
!>          The first (N-1) elements of LPLUS contain the subdiagonal
!>          elements of the unit bidiagonal matrix L(+).
!> \endverbatim
!>
!> \param[out] WORK
!> \verbatim
!>          WORK is REAL array, dimension (2*N)
!>          Workspace.
!> \endverbatim
!>
!> \param[out] INFO
!> \verbatim
!>          INFO is INTEGER
!>          Signals processing OK (=0) or failure (=1)
!> \endverbatim
!
!  Authors:
!  ========
!
!> \author Univ. of Tennessee
!> \author Univ. of California Berkeley
!> \author Univ. of Colorado Denver
!> \author NAG Ltd.
!
!> \date June 2016
!
!> \ingroup OTHERauxiliary
!
!> \par Contributors:
!  ==================
!>
!> Beresford Parlett, University of California, Berkeley, USA \n
!> Jim Demmel, University of California, Berkeley, USA \n
!> Inderjit Dhillon, University of Texas, Austin, USA \n
!> Osni Marques, LBNL/NERSC, USA \n
!> Christof Voemel, University of California, Berkeley, USA
!
!  =====================================================================
      SUBROUTINE SLARRF(N,D,L,Ld,Clstrt,Clend,W,Wgap,Werr,Spdiam,Clgapl,&
     &                  Clgapr,Pivmin,Sigma,Dplus,Lplus,Work,Info)
      IMPLICIT NONE
!*--SLARRF195
!
!  -- LAPACK auxiliary routine (version 3.7.1) --
!  -- LAPACK is a software package provided by Univ. of Tennessee,    --
!  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
!     June 2016
!
!     .. Scalar Arguments ..
      INTEGER Clstrt , Clend , Info , N
      REAL Clgapl , Clgapr , Pivmin , Sigma , Spdiam
!     ..
!     .. Array Arguments ..
      REAL D(*) , Dplus(*) , L(*) , Ld(*) , Lplus(*) , W(*) , Wgap(*) , &
     &     Werr(*) , Work(*)
!     ..
!
!  =====================================================================
!
!     .. Parameters ..
      REAL MAXGROWTH1 , MAXGROWTH2 , ONE , QUART , TWO
      PARAMETER (ONE=1.0E0,TWO=2.0E0,QUART=0.25E0,MAXGROWTH1=8.E0,      &
     &           MAXGROWTH2=8.E0)
!     ..
!     .. Local Scalars ..
      LOGICAL dorrr1 , forcer , nofail , sawnan1 , sawnan2 , tryrrr1
      INTEGER i , indx , ktry , KTRYMAX , SLEFT , SRIGHT , shift
      PARAMETER (KTRYMAX=1,SLEFT=1,SRIGHT=2)
      REAL avgap , bestshift , clwdth , eps , fact , fail , fail2 ,     &
     &     growthbound , ldelta , ldmax , lsigma , max1 , max2 ,        &
     &     mingap , oldp , prod , rdelta , rdmax , rrr1 , rrr2 ,        &
     &     rsigma , s , smlgrowth , tmp , znm2
!     ..
!     .. External Functions ..
      LOGICAL SISNAN
      REAL SLAMCH
      EXTERNAL SISNAN , SLAMCH
!     ..
!     .. External Subroutines ..
      EXTERNAL SCOPY
!     ..
!     .. Intrinsic Functions ..
      INTRINSIC ABS
!     ..
!     .. Executable Statements ..
!
      Info = 0
!
!     Quick return if possible
!
      IF ( N<=0 ) RETURN
!
      fact = REAL(2**KTRYMAX)
      eps = SLAMCH('Precision')
      shift = 0
      forcer = .FALSE.

!     Note that we cannot guarantee that for any of the shifts tried,
!     the factorization has a small or even moderate element growth.
!     There could be Ritz values at both ends of the cluster and despite
!     backing off, there are examples where all factorizations tried
!     (in IEEE mode, allowing zero pivots & infinities) have INFINITE
!     element growth.
!     For this reason, we should use PIVMIN in this subroutine so that at
!     least the L D L^T factorization exists. It can be checked afterwards
!     whether the element growth caused bad residuals/orthogonality.

!     Decide whether the code should accept the best among all
!     representations despite large element growth or signal INFO=1
!     Setting NOFAIL to .FALSE. for quick fix for bug 113
      nofail = .FALSE.
!

!     Compute the average gap length of the cluster
      clwdth = ABS(W(Clend)-W(Clstrt)) + Werr(Clend) + Werr(Clstrt)
      avgap = clwdth/REAL(Clend-Clstrt)
      mingap = MIN(Clgapl,Clgapr)
!     Initial values for shifts to both ends of cluster
      lsigma = MIN(W(Clstrt),W(Clend)) - Werr(Clstrt)
      rsigma = MAX(W(Clstrt),W(Clend)) + Werr(Clend)

!     Use a small fudge to make sure that we really shift to the outside
      lsigma = lsigma - ABS(lsigma)*TWO*eps
      rsigma = rsigma + ABS(rsigma)*TWO*eps

!     Compute upper bounds for how much to back off the initial shifts
      ldmax = QUART*mingap + TWO*Pivmin
      rdmax = QUART*mingap + TWO*Pivmin

      ldelta = MAX(avgap,Wgap(Clstrt))/fact
      rdelta = MAX(avgap,Wgap(Clend-1))/fact
!
!     Initialize the record of the best representation found
!
      s = SLAMCH('S')
      smlgrowth = ONE/s
      fail = REAL(N-1)*mingap/(Spdiam*eps)
      fail2 = REAL(N-1)*mingap/(Spdiam*SQRT(eps))
      bestshift = lsigma
!
!     while (KTRY <= KTRYMAX)
      ktry = 0
      growthbound = MAXGROWTH1*Spdiam

 100  sawnan1 = .FALSE.
      sawnan2 = .FALSE.
!     Ensure that we do not back off too much of the initial shifts
      ldelta = MIN(ldmax,ldelta)
      rdelta = MIN(rdmax,rdelta)

!     Compute the element growth when shifting to both ends of the cluster
!     accept the shift if there is no element growth at one of the two ends

!     Left end
      s = -lsigma
      Dplus(1) = D(1) + s
      IF ( ABS(Dplus(1))<Pivmin ) THEN
         Dplus(1) = -Pivmin
!        Need to set SAWNAN1 because refined RRR test should not be used
!        in this case
         sawnan1 = .TRUE.
      ENDIF
      max1 = ABS(Dplus(1))
      DO i = 1 , N - 1
         Lplus(i) = Ld(i)/Dplus(i)
         s = s*Lplus(i)*L(i) - lsigma
         Dplus(i+1) = D(i+1) + s
         IF ( ABS(Dplus(i+1))<Pivmin ) THEN
            Dplus(i+1) = -Pivmin
!           Need to set SAWNAN1 because refined RRR test should not be used
!           in this case
            sawnan1 = .TRUE.
         ENDIF
         max1 = MAX(max1,ABS(Dplus(i+1)))
      ENDDO
      sawnan1 = sawnan1 .OR. SISNAN(max1)

      IF ( forcer .OR. (max1<=growthbound .AND. .NOT.sawnan1) ) THEN
         Sigma = lsigma
         shift = SLEFT
         GOTO 300
      ENDIF

!     Right end
      s = -rsigma
      Work(1) = D(1) + s
      IF ( ABS(Work(1))<Pivmin ) THEN
         Work(1) = -Pivmin
!        Need to set SAWNAN2 because refined RRR test should not be used
!        in this case
         sawnan2 = .TRUE.
      ENDIF
      max2 = ABS(Work(1))
      DO i = 1 , N - 1
         Work(N+i) = Ld(i)/Work(i)
         s = s*Work(N+i)*L(i) - rsigma
         Work(i+1) = D(i+1) + s
         IF ( ABS(Work(i+1))<Pivmin ) THEN
            Work(i+1) = -Pivmin
!           Need to set SAWNAN2 because refined RRR test should not be used
!           in this case
            sawnan2 = .TRUE.
         ENDIF
         max2 = MAX(max2,ABS(Work(i+1)))
      ENDDO
      sawnan2 = sawnan2 .OR. SISNAN(max2)

      IF ( forcer .OR. (max2<=growthbound .AND. .NOT.sawnan2) ) THEN
         Sigma = rsigma
         shift = SRIGHT
         GOTO 300
      ENDIF
!     If we are at this point, both shifts led to too much element growth

!     Record the better of the two shifts (provided it didn't lead to NaN)
!        both MAX1 and MAX2 are NaN
      IF ( sawnan1 .AND. sawnan2 ) GOTO 200
      IF ( .NOT.sawnan1 ) THEN
         indx = 1
         IF ( max1<=smlgrowth ) THEN
            smlgrowth = max1
            bestshift = lsigma
         ENDIF
      ENDIF
      IF ( .NOT.sawnan2 ) THEN
         IF ( sawnan1 .OR. max2<=max1 ) indx = 2
         IF ( max2<=smlgrowth ) THEN
            smlgrowth = max2
            bestshift = rsigma
         ENDIF
      ENDIF

!     If we are here, both the left and the right shift led to
!     element growth. If the element growth is moderate, then
!     we may still accept the representation, if it passes a
!     refined test for RRR. This test supposes that no NaN occurred.
!     Moreover, we use the refined RRR test only for isolated clusters.
      IF ( (clwdth<mingap/REAL(128)) .AND. (MIN(max1,max2)<fail2) .AND. &
     &     (.NOT.sawnan1) .AND. (.NOT.sawnan2) ) THEN
         dorrr1 = .TRUE.
      ELSE
         dorrr1 = .FALSE.
      ENDIF
      tryrrr1 = .TRUE.
      IF ( tryrrr1 .AND. dorrr1 ) THEN
         IF ( indx==1 ) THEN
            tmp = ABS(Dplus(N))
            znm2 = ONE
            prod = ONE
            oldp = ONE
            DO i = N - 1 , 1 , -1
               IF ( prod<=eps ) THEN
                  prod = ((Dplus(i+1)*Work(N+i+1))/(Dplus(i)*Work(N+i)))&
     &                   *oldp
               ELSE
                  prod = prod*ABS(Work(N+i))
               ENDIF
               oldp = prod
               znm2 = znm2 + prod**2
               tmp = MAX(tmp,ABS(Dplus(i)*prod))
            ENDDO
            rrr1 = tmp/(Spdiam*SQRT(znm2))
            IF ( rrr1<=MAXGROWTH2 ) THEN
               Sigma = lsigma
               shift = SLEFT
               GOTO 300
            ENDIF
         ELSEIF ( indx==2 ) THEN
            tmp = ABS(Work(N))
            znm2 = ONE
            prod = ONE
            oldp = ONE
            DO i = N - 1 , 1 , -1
               IF ( prod<=eps ) THEN
                  prod = ((Work(i+1)*Lplus(i+1))/(Work(i)*Lplus(i)))    &
     &                   *oldp
               ELSE
                  prod = prod*ABS(Lplus(i))
               ENDIF
               oldp = prod
               znm2 = znm2 + prod**2
               tmp = MAX(tmp,ABS(Work(i)*prod))
            ENDDO
            rrr2 = tmp/(Spdiam*SQRT(znm2))
            IF ( rrr2<=MAXGROWTH2 ) THEN
               Sigma = rsigma
               shift = SRIGHT
               GOTO 300
            ENDIF
         ENDIF
      ENDIF

 200  IF ( ktry<KTRYMAX ) THEN
!        If we are here, both shifts failed also the RRR test.
!        Back off to the outside
         lsigma = MAX(lsigma-ldelta,lsigma-ldmax)
         rsigma = MIN(rsigma+rdelta,rsigma+rdmax)
         ldelta = TWO*ldelta
         rdelta = TWO*rdelta
         ktry = ktry + 1
         GOTO 100
!        None of the representations investigated satisfied our
!        criteria. Take the best one we found.
      ELSEIF ( (smlgrowth<fail) .OR. nofail ) THEN
         lsigma = bestshift
         rsigma = bestshift
         forcer = .TRUE.
         GOTO 100
      ELSE
         Info = 1
         RETURN
      ENDIF

 300  IF ( shift==SLEFT ) THEN
      ELSEIF ( shift==SRIGHT ) THEN
!        store new L and D back into DPLUS, LPLUS
         CALL SCOPY(N,Work,1,Dplus,1)
         CALL SCOPY(N-1,Work(N+1),1,Lplus,1)
      ENDIF

!
!     End of SLARRF
!
      END SUBROUTINE SLARRF


text:
      SUBROUTINE wrout(fz1s, fz1c, fz3s, fz3c, srces)

c  This SUBROUTINE creates the file DKESF containing the magnetic field
c  and the final solution distributions.

C-----------------------------------------------
C   M o d u l e s
C-----------------------------------------------
      USE Vnamecl2
      USE safe_open_mod
      USE dkes_input
      USE dkes_realspace, ONLY: mvalue, nvalue, mpnt, borbi1
      IMPLICIT NONE
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      REAL(rprec), DIMENSION(mpnt,0:lalpha) :: fz1s, fz1c, fz3s, fz3c
      REAL(rprec), DIMENSION(mpnt,4,2,2) :: srces
C-----------------------------------------------
C   L o c a l   P a r a m e t e r s
C-----------------------------------------------
      INTEGER, PARAMETER :: idkes = 7, iobin = 66
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER :: m, n, mn, ll, l, iodkes = idkes
C-----------------------------------------------

      CALL safe_open(iodkes, mn, 'DKESF', 'unknown', 'formatted')
      IF (mn .ne. 0) STOP 'ERROR OPENING DKESF FILE'

      WRITE (iodkes, 100) chip, psip, nzperiod, cmul1, efield1, wtov,
     1   weov, mpolb, ntorb, lfout, mpnt, blabl(ibbi)
  100 FORMAT(/8x,'CHIP',8x,'PSIP',4x,'NZPERIOD',8x,'CMUL',6x,'EFIELD',4x
     1   ,'OMEGAT/V',4x,'OMEGAE/V',7x,'MPOLB',7x,'NTORB',7x,'LFOUT',2x,
     2   'BLOCK SIZE'//1p,2e12.4,i12,4e12.4,4i12,2/4x,'M',4x,'N',14x,a3,
     3   18x,'BTHETA',7x,'BZETA'/)

!     WRITE FIRST ROW

      DO mn = 1, mpnt
         m = mvalue(mn)
         n = nvalue(mn)
         IF (m.ne.0 .or. n.ne.0) CYCLE
         WRITE (iodkes, 200) m, n, borbi1(mn), btheta, bzeta
         EXIT
      END DO
  200 FORMAT(2i5,5x,1pe12.4,12x,2e12.4)

      DO mn = 1, mpnt
         m = mvalue(mn)
         n = nvalue(mn)
         IF (m.ne.0 .and. n.eq.0 .and. borbi1(mn).ne.zero) THEN
            WRITE (iodkes, 300) m, n, borbi1(mn)
         END IF
      END DO
  300 FORMAT(2i5,5x,1pe12.4)

      DO n = 1, ntorb
         DO mn = 1, mpnt
           IF (nvalue(mn).eq.n .and. borbi1(mn).ne.zero) THEN
               WRITE (iodkes, 300) mvalue(mn),n,borbi1(mn)
           END IF
         END DO
      END DO

      IF (ipmb .eq. 1) THEN
         fz1c(:,:lfout) = zero
         fz3c(:,:lfout) = zero
      ENDIF
      IF (ipmb .eq. 2) THEN
         fz1s(:,:lfout) = zero
         fz3s(:,:lfout) = zero
      ENDIF

      n = nzperiod
      IF (n .eq. 0) n = 1

!
!     WRITE "CONSTRAINT" PHANTOM DISTRIBUTION (l = -1, remember l is offset by 1 here)
!
      WRITE (iodkes, 400)
  400 FORMAT(/,'  Constraint Phantom Distribution',//,
     &   4x,'M',4x,'N',13x,'F11C',8x,'F11S',8x,'F13C',8x,'F13S'/)
      WRITE (iodkes, 500) (mvalue(mn),nvalue(mn),fz1s(mn,0),fz1c(mn,0),
     1   fz3s(mn,0),fz3c(mn,0),mn=1,mpnt)
  500 FORMAT(2i5,5x,1p,4e12.4)

!
!     WRITE DISTRIBUTION FUNCTION FOURIER COEFFIENTS FOR EACH LEGENDRE PITCH
!     IN NAMELIST, SET idisk:
!        idisk = 0 : DKES computes all l values.
!        idisk = 1 : DKES computes only l = 0,1,2,3,4,5 (default)
!
!     THESE ARE COEFFICIENTS OF ORTHONORMALIZED LEGENDRE FUNCTIONS and FOURIER HARMONICS
!
      OPEN (unit=iobin, file='dkes_bin', status='replace',
     &      form='unformatted')
      WRITE (iodkes, 600)
  600 FORMAT(/, '  Fourier Harmonics of Distribution Functions',//,
     1   4x,'M',4x,'N',5x,'F01C',8x,'F01S',8x,'F03C',8x,
     1   'F03S',14x,'S01C',8x,'S01S',8x,'S03C',8x,'S03S'/)

      WRITE (iobin) mpnt, lfout
      l = -1
      DO ll = 1, lfout
         l = l + 1
         WRITE (iodkes, 650) l
         WRITE (iobin) l
         WRITE (iobin) (mvalue(mn), nvalue(mn), fz1c(mn,ll),
     &               fz1s(mn,ll), fz3c(mn,ll), fz3s(mn,ll), mn=1,mpnt)
         IF (ll > 4) WRITE (iodkes, 700) (mvalue(mn),nvalue(mn),
     1      fz1c(mn,ll),fz1s(mn,ll),fz3c(mn,ll),fz3s(mn,ll),mn=1,mpnt)
         IF (ll <= 4)
     1      WRITE (iodkes, 800) (mvalue(mn),nvalue(mn),
     2      fz1c(mn,ll),fz1s(mn,ll),fz3c(mn,ll),fz3s(mn,ll),
     3      srces(mn,ll,1,2),srces(mn,ll,1,1),srces(mn,ll,2,2),
     4      srces(mn,ll,2,1),mn=1,mpnt)
      END DO

  650 FORMAT(/,2x,'PITCH L =',i5,/)
  700 FORMAT(2i5,1p,4e12.4)
  800 FORMAT(2i5,1p,4e12.4,6x,4e12.4)

      CLOSE(iodkes)
      CLOSE(iobin)

      END SUBROUTINE wrout

ext: f

max_stars_repo_name: joseluisvelasco/STELLOPT

max_issues_repo_name: joseluisvelasco/STELLOPT

max_forks_repo_name: joseluisvelasco/STELLOPT


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

I meant to say WON, college square was voted best apartments in davis for  2008
20081020 21:38:25 nbsp Welcome to the Wiki Howdy Ann, and Welcome to the Wiki!     My names Evan, pleased to meet you!  Thanks for adding your take on College Square Apartments. Users/JabberWokky Evan JabberWokky Edwards

20081021 17:04:43 nbsp Hi, I fixed your comment for College Square Apartments. You can fix typos like that by clicking the edit button at the top of the page.  Users/ElleWeber

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
The ASUCD Perestroika was a period from Fall 2004  Spring 2005, in which ASUCD underwent rapid change.

During this era of ASUCD History, there was a period of unprecedented interest in ASUCD reforms.  Persons without official positions within ASUCD are enthusiastically participating in the writing of bills to change perceived injustice in ASUCD to a degree markedly higher than any time in known ASUCD history.  Some have declared this renaissance the ASUCD wiki:wikipedia:Perestroika.  Those some are named Kris Fricke, but its a good enough name that we can use it too.

 Setting the Stage

At the beginning of the 2004 calendar year, everything looked pretty good for Student Focus  they were on top, and controlled the wheels of student government.  Well, there was that pesky History of Choice Voting in Davis Choice Voting Amendment passed in 2003, which had allowed a couple of LEAD senators to get elected, but that was it.

But for Student Focus, the perfect storm was coming.  L.E.A.D. Senator Caliph Assagai was well connected to real politicians through his parents.  And a certain Davis College Democrats DCD member, James Schwab had proven himself an excellent organizer.  The two them combined to plan an excellent campaign strategy for 2004.   But the first true sign of problems was the Senator Ackerman Scandal, in which StuFo Senator James Ackerman resigned over having stolen books from the UCD Bookstore, which has like 10 million cameras.FootNote(Of course, to this, wiki:wikipedia:True Grit Rooster Cogburn would have noted that You cant rob from a bunch of thieves, right?)

From a secret computer lab hidden on campusFootNote(Actually it was Hoagland Hall 124), the Students for an Orwellian Society, a crazy group with too much time on their hands decided: Well, if Orwellians are running the federal government, why arent we running the student government?  They formed SOSSS, a joke slate with a few serious issues, told everyone how Big Brother wanted them to vote down to the correct order, and made some silly flyers.  And they attacked Ackerman, who just happened to be the boyfriend of the sitting Elections Committee Chair, Christine Schacter.  As the election came to a vote, this spawned the Election Statement Controversy from her Committees editing of ballot statements.

Immediately after the election, Student Focus Senatorselect Sean Ruel and Nafeh Malik were discovered to have been campaigning in the dorms, leading to the aptly named Campaigning in Dorms Controversy.  While solicitation is a clear violation of Student Housing policy, the implication that candidates were helping people with their ballots went beyond the pale.  The ASUCD Court was literally powerless to do anything about either scandal.  And though the dirty pair resigned, likely due to pressure from Student Judicial Affairs, the stage was now set for a revolution.

 Perestroika

Due to the double resignation, eight seats were up for grabs in the Winter 2005 ASUCD Election, as well as the presidency, so a sea change was possible.  LEAD developed an even better organizational network, and recruited six candidates to run.

The Orwellians were not the kind of people that you wanted to push around.  Led by Users/RevChad Chad van Shoelandt and Users/BrentLaabs Brent Laabs, they were incensed at their treatment  and the treatment of democracy.  So they redoubled their efforts, and became determined to fix the problem.  Recruiting local celebrity Users/RobRoy to help and seeking the advice of Chief Justice Users/KrisFricke Fricke, they had an eightpage elections reform bill introduced that drastically reduced the power of the Elections Committee.  And they launched a brand new slate, Friends Urging Campus Kindness with the help of the brandnew Davis Wiki.  The Friends, who for some reason didnt use their acryonym, originally planned on running 23 candidates in a zerg rush attack on the Senate.  However, most of them backed out (or were ineligible for being Users/BrentLaabs a grad student), so they ended up with seven.  They endorsed Assagai, LEADs presidential candidate, due to lack of hubris on their own part.

The results of the election were split: 2 LEAD, 2 Focus, 2 Friends, 2 Independent Independents.  But Rob Roy became the first candidate ever elected in the first round, and Caliph Assagai won election.  The results were clear, and though they knew it not, the era of Student Focus dominance had ended for good.

 Reforms

Many reforms came to pass in the next few months.  Most of the legislation was fairly long, but the bylaws had been somewhat neglected since the ASUCD Federal Era (around 1995).  The ASUCD Court ASUCD Supreme Court was granted authority over some elections cases (questions of law, but not questions of fact), since the Campus Judicial Board made it clear that they thought judging internal ASUCD cases really shouldnt be part of their jobs.

The power of the ASUCD Elections Commission was clearly defined, as were behavior around polling places.  Penalties were more clearly defined.  It of course could not prevent all scandals in ASUCD, but seriously dudes, the codes were a fucking mess before then.

And further, attempts at transparency were pushed to new levels.  The ASUCD Website started to be used more, though it took years to finally get things updated regularly.  Eventually, it would become a resource for insiders as well, which is good for transparency when everyone uses the same basis for discussion.

 Problems

There was that one time where the Senate decided to resolve an elections dispute with an ex post facto law, just because they could.  That, of course, led to a constitutional amendment banning ex post facto laws.

 Legacy

Though the era was short lived, many of the key players were integrated into the government at some level.  SOSSSs #1 candidate, Jonathon Leathers became Elections Committee Chair in a rare show of irony from the Senate.

The focus on transparency, fairness, and efficiency continued into the following Harms Hegemony era, as technocrats played a much greater role in ASUCD politics during that time.

Yet for the radicals who had sought to break the slate system, the end results were disappointing.  The reforms they had pushed through made it easier for the incoming LEAD group to survive, and in fact, thrive.  Still, Choice Voting endured, with LEADs blessing, which allowed opposition Senators to endure.  Meanwhile, LEAD became exceptionally good at choosing candidates to support their broad coalition, and in training successors; they generally prospered under Choice Voting.

 Documents

The original purpose of this page was to add some wiki:wikipedia:glasnost to the perestroika in hopes of garnering yet more public involvement.  The idea was to make some of the current reform proposals available here for public perusal and contribution to the writing process.  It worked for the most part, as a few more people than the usual suspects participated.

So without further ado, and in their full but unfinished glory, are various Senate bills to reform the ASUCD government:

 Constitutional Amendment Proposals
 ASUCD Judicial Reform Act
Current Status: Passed unanimously by Senate, will appear on Winter ballot
Contacts: Senator Darnell Holloway (djholloway@ucdavis.edu) or Justice Daniel Raff
Current Version: version posted here (v2) hella out of date  will post current (v5) when I can
Other Commentary: other election reform legislation also being currently worked on intends to render this change.  As this necessitates a change to the http://asucd.ucdavis.edu/gov/docs/ ASUCD Constitution, this separate constitutional amendment is necessary.  It must be passed by the ASUCD Senate Senate in order to appear on the winter 2005 ballot, and passed by a majority of voters to take effect.

Authored by:
CoAuthored by:
Introduced by:

A bill to add this ASUCD Judicial Reform Act, a constitutional amendment, to the Winter 2005 election ballot.

Background: 	In 2001, the elections process within ASUCD was taken entirely from the jurisdiction of the ASUCD Court.  The Elections Committee was given direct jurisdiction over alleged campaign violations.  Appeal of the Committees decisions on these cases could be made to the Campus Judicial Board (Campus Judicial Board) of Student Judicial Affairs (Student Judicial Affairs), a nonASUCD body.

In effect this has shielded the Elections Committee from judicial review within ASUCD.  wiki:WikiPedia:Original jurisdiction over elections complaints fell upon a body that had not been assembled with an interest in judicial matters.  The primary purpose of the Elections Committee is the administration of elections, and no members are hired for expertise or interest in jurisprudence.  Furthermore the appeal process from the Elections Committee leads not to the Court but to SJA, a place far removed from ASUCD experience that most students are fortunate enough to know little about.  And if the CJB of SJA were to actually issue a binding opinion involving orders upon the Elections Committee, we would have the ideologically repugnant situation of ASUCD operations being ordered by a nonASUCD entity.  CJB has avoided this however, refusing to issue binding orders pursuant to its wiki:wikipedia:Appellate jurisdiction appellate jurisdiction over the Elections Committee, which creates the similarly alarming circumstance of no entity actually exercising authority over the Committee.

This maldesigned system has since 2001 led to increasing cynicism and distrust of the electoral process in ASUCD.  By the 2004 Fall election this erupted into a Election Statement Controversy major controversy, with at least five cases coming before SJA, numerous Aggie articles, and significant damage to ASUCDs reputation.

This bill seeks to rectify this situation by placing the ASUCD Court in the appellate role rather than SJA.  The Elections Committee shall retain original jurisdiction over elections complaints, and its decisions may be appealed to the ASUCD Court regarding matters of law.  This means that it shall remain the Elections Committees role to investigate allegations, and the Court shall not interfere with its decisions as to the wiki:wikipedia:Question of fact preponderance of evidence.  The Committees interpretation of or adherence to the bylaws, however, may be challenged in an appeal before the ASUCD Court.  In this manner the administration of elections and the interpretation of jurisprudence are placed with the entities most qualified in their respective areas, there is a checkandbalance upon the elections process, and ASUCD operations shall not be ordered by nonASUCD bodies.

Article VII

Section 1	The judicial authority of the ASUCD shall be vested in one ASUCD Court (hereafter referred to as the ASUCD Court) Xand the Student Judicial Affairs / Campus Judicial BoardX.

Sections 8 through 12 deleted __Section 8	The ASUCD Court shall have appellate jurisdiction over elections.  The Elections Committee shall retain original jurisdiction over potential disqualifications or violations in elections, but specific questions of law may be appealed to the ASUCD Court.__

 Countback
Countback has its Countback own page

 Open Government
It has been proposed that a bill be written requiring that any legislation heard by the Senate be available online for at least three days before it is voted on by the Senate, to prevent ASUCD Open Government stealthy questionable $33,000 allocations.

Status: The foundation to solve this problem was first created in Spring 2003 when Paul Amnuaypayoat wrote a bill to put government documents online in the first place. Unfortunately, compliance has been rather half hearted lately, and the bill Paul Amnuaypayoat Paul wrote did not address certain issues, such as the one above. This problem is currently being worked on by Jonathon Leathers. Check out his page for the actual text of the bill to solve this problem.

In the interests of an open and transparent election process, it has also been proposed that ASUCD post the underlying ballot data at the same time election results are announced.  The city of San Francisco does this as a matter of course for its ranked voting elections.  In the past, ASUCD has been either unwilling or reluctant to release its ballot data, even though this language is spelled out in the Constitution and further amplified in the Government Codes.  Several attempts have been made to accomplish this, with some gains.

 Impeachment
Currently ASUCD elected officials cannot be impeached, whereas in the socalled real world they almost invariably can be.

Status: not being worked on

 Proposed Changes to Other Bylaws
 Judicial Codes
 Legal Enfranchisement
Current Status: Introduced, will go to Internal Affairs Commission at their meeting of Monday, 02/07/05, (5:30pm, one of the outdoor rooms on the second floor of the MU), and be discussed in Senate at the following meeting (02/10/05)
Contacts: Senator Darnell Holloway (djholloway@ucdavis.edu) or Chief Justice Users/KrisFricke Kris Fricke
Current Version: this is the first version
Other Commentary: none comes to mind

Authored by:  Holloway			 Coauthored by: 	 Introduced by:	Holloway

Referred to: Internal Affairs Commission

A bill to amend the ASUCD Standing Rules to include the ASUCD judiciary in ASUCD legal matters.

BACKGROUND: A constant criticism of judicial authority in ASUCD is that it might come into conflict with legal circumstances related to the nature of ASUCD.  Many of these mystical issues of potential litigation are barely understood by anyone in ASUCD, and from those issues that do arise, nothing is permitted to permeate down to the Court. If the ASUCD Court is party to these legal discussions the Senate and Executive Office participate in, they shall gain a crucial understanding of the legal circumstances of ASUCD.  Whereas members of the Senate and the Executive Office are rarely involved in ASUCD leadership more than two years at most, the members of the Court can carry the experience gleamed from legal discussions through time, being in their positions for much longer.

Section 1. The Senate hereby amends ASUCD Rule Eight, Section 5.4.

__Section 5.4.3

The Members of the ASUCD Court shall sit with the Senate as ex officio members during closed sessions on litigation matters.  They shall take part in discussion but as long as there are more than one legal options, not advocate one legal course of action over another, and shall not vote.  Members of the Court will be notified and invited as soon as the litigation matter is added to the agenda.

Section 5.4.4

The Chief Justice and Alternate Presiding Justice of the ASUCD Court shall be invited to participate whenever the ASUCD Executive Office discusses legal matters with legal council. They shall take part in discussion but as long as there are more than one legal options, not advocate one legal course of action over another.  These members of the Court shall be notified and invited as soon as the meeting is scheduled with legal council.__

 Legislative Enfranchisement
Current Status: Introduced, will go to Internal Affairs Commission at their meeting of Monday, 02/07/05, (5:30pm, one of the outdoor rooms on the second floor of the MU), and be discussed in Senate at the following meeting (02/10/05)
Contacts: Senator Janine Fiel (jmfiel@ucdavis.edu) or Chief Justice Users/KrisFricke Kris Fricke
Current Version: an earlier version of this was passed by over twothirds of Senate, but vetoed by Users/KalenGallagher President Gallagher
Other Commentary: CJB may be substituted in place of Senate here.  Also, part in italics isnt officially in the bill yet but will be added in IAC

BACKGROUND: The ASUCD Court is in need of various reforms to improve its functionality. Members of the ASUCD court have the most working knowledge of the Court, yet they have limited ability to speak on their behalf when the Senate is reviewing legislation related to court operations. Currently, members of the Court are forced to go through a third party to author legislation, who in turn writes almost exactly what the members tell them to. If ASUCD adopts the amendment in Section 1, Court members will be able to speak during authors comments and field questions when the Senate is reviewing legislation pertaining to court operations.

Section 1. The Senate hereby amends ASUCD Judicial Code Six, Section 605.

Section 605

__Members of the ASUCD Court may author, or coauthor, legislation to amend the Judicial Codes and/or article VII of the ASUCD Constitution.

a)	Any challenges to the constitutionality of legislation written by the Court, or of official actions of the Court, shall be considered by the Senate, except in the circumstance described in subsection c

b)	The Senate shall hold a special hearing in accordance with Judicial Code Seven, with the Senators acting as justices and the Senate president pro tempore acting as Chief Justice.

c)	If the ASUCD Court never discussed the legislation as a body, or a quorum of Court members can be gathered who abstained when the Court considered the legislation, were not present, or were not on the Court at the time, those Court members shall consider challenges to the legislation.__

 Right to Counsel
Current Status: Could be introduced as early as next Thursday (03/10/05), will be discussed by the ASUCD Judiciary this evening (Monday the 7th) at 8:30pm in the ASUCD Conf. Rm., 3rd Floor MU.
Contacts: Users/KrisFricke, Chief Justice
Current Version: Second draft
Other Commentary: the idea of having persons in ASUCD whose job it is to represent parties in the ASUCD judicial system was first highlighted as a necessity when last year the Court faced several named defendants who had no studen head to defend themselves (SGAO & Creative Media).  Having someone who would defend the ASUCD govt at large would have been highly useful.  Also the Court cannot call out potential cases it sees but the proposed Plaintiff Advocate could.  I should have a full background for this legislation written out later this week.

Legislation can be found http://www.geocities.com/winged_snail/judref04b.doc HERE because Im not about to try to type all that out again.

 Judicial Accountability and Decisions on Elections (JADE) Act
This proposal has two aspects to it: (1) update the codes to match Winter 2005 ASUCD Election/Ballot Measures Constitutional Amendment #6 and (2) allow the ASUCD Court to serve as the court of first resort in cases where ASUCD is breaking state or federal law.  This bill is intended to help avoid things like the Election Statement Controversy and prevent future Lamargatetype incidents.

Current version (1.0) is available http://atm.ucdavis.edu/~bslaabs/JADEAct.doc here.

 Elections Codes
See ASUCD Elections Reform

20050122 08:07:53 nbsp Somehow the phrase Campus Judicial Board of Student Judicial Affairs was used... and now it seems to have stuck. The name is incorrect, and implies that CJB is merely an arm of SJA  which is false. I wrote the Aggie about this (because their articles make the same mistake and I felt it was potentially misleading to their readers) and they apologized, but did not publish my letter. They also continue to make the same mistake. shrug Users/ElvinLee

20050204 18:25:00 nbsp If anyone can get a copy of my passed bill from 2003, that would be a really good first step in solving the open government problem. I am willing to help if someone can go into the SGAO during its open hours to get this important legislation for me. Paul Amnuaypayoat Paul

20050205 20:21:11 nbsp Legislative Enfranchisement is probably a bad idea. I dont think the Senate should decide the constitutionality of anything, so if youre not conforatable with the court deciding the constitutionality of a bill it authored then you shouldnt let them author bills. The US Congress doesnt let the US Supreme court author legislation after all. (Actually, that last sentence may not be right. I think anyone in the US can author legislation, but the limits are on who can introduce it  and nothing prevents the US Supreme court from deciding the constitutionality of anything. Users/KenBloom

20050205 21:14:14 nbsp Yeah, Ken, Im pretty sure anyone can write anything, but it has to be introduced by a member of the House or Senate.  Theres nothing in the Constitution, anyway, which prohibits the executive or judiciary from writing laws.  The reason the Supreme Court never writes laws is due in part to judicial activism.  But the main reason is that they can just have one of their 3,652 lawyer friends write it for them.   In regards to this issue, I dont see a problem with the court reviewing the constitutionality of a bill that some members wrote.  Some members might want to recuse themselves for the sake of personal integrity. Users/BrentLaabs

20050207 01:31:13 nbsp Incidently, in the United States the Supreme Court of the United States (SCOTUS) HAS been directly responsible for writing legislation, the Federal Rules of Civil Procedure come directly to mind.  In fact, anything authored by the SCOTUS is assumed to be constitutional by its very nature, how bout THEM apples?  Anyway, the fact is that our ASUCD Senate seems reluctant to allow members of the Court the same writes of authorship everyone else enjoys without some kind of compromise  Im leaning more towards sending things to CJB now rather than the Senate.  Thanks for your interest though and please contact me if you want to discuss it further (and or continue discussion here)   Users/KrisFricke


     program demo_factor
     use M_calcomp, only : plots, plot, number, symbol, newpen
     use M_calcomp, only : nframe, factor, rect
     use M_calcomp, only : MOVE, END, DRAW
     call plots(0.0,10.0,0.0,10.0)
     call draw_car_prices()
     call nframe()
     call factor(0.5)
     call draw_car_prices()
     call plot( 5.0, 5.0,-3)
     call factor(0.5)
     call draw_car_prices()
     call plot(0.0,0.0,end)
     contains
     subroutine draw_car_prices()
        character(len=21) :: ichr6
        character(len=19) :: ichr7
        character(len=17) :: ichr8
        ichr6='CAR MODEL AGE (YEARS)'
        ichr7='CAR VALUE (DOLLARS)'
        ichr8='AVERAGE CAR VALUE'
        !     CALL TO SYMBOL USES -0.5Y, -0.8-.14  X
        !     (-.14 FOR CHARACTER HEIGHT)
        call rect(0.0,0.0,10.0,10.0,0.0,7)
        call plot(0.95,0.5,-MOVE)
        ! PLOT CAR VALUE CHART WITHOUT USING SCALE,AXIS,OR LINE
        x=1.0
        ! PLOT X-AXIS
        do i=1,7
           call plot(x-1.0,0.0,MOVE)
           call plot(x   , 0.0,DRAW)
           call plot(x   ,-0.1,DRAW)
           call number(x-.02,-0.25,0.1,x,0.0,-1)
           x=x+1.0
        enddo
        call symbol(2.0,-0.5,0.14,ichr6,inteq,0.0,21)
        ! PLOT Y-AXIS
        value=1000.0
        do i=1,6
           y=0.0015*value
           call plot(0.0,y-1.5,MOVE)
           call plot(0.0,y-.75,DRAW)
           call plot(-.1,y-.75,DRAW)
           call plot(0.0,y-.75,DRAW)
           call plot(0.0,y    ,DRAW)
           call plot(-.1,y    ,DRAW)
           call number(-0.7,y,0.14,value,0.0,-1)
           value=value+1000.0
        enddo
        call symbol(-0.8,3.1,0.14,ichr7,inteq,90.0,19)
        ! PLOT CURVES
        call newpen(2)
        do i=2000,6000,500
           value=i
           age=0.0
           call plot(age,0.0015*value,MOVE)
           do j=1,84
              value=value*0.972
              age=age+0.08333
              call plot(age,0.0015*value,DRAW)
           enddo
        enddo
        call newpen(3)
        call symbol(3.0,6.0,0.21,ichr8,inteq,0.0,17)
        end subroutine draw_car_prices
        end program demo_factor


text:
      subroutine      SWATAN (N5,U5,Y5)
c Copyright (c) 1996 California Institute of Technology, Pasadena, CA.
c ALL RIGHTS RESERVED.
c Based on Government Sponsored Research NAS7-03001.
C>> 2001-06-18 SWCOMP Krogh   Replaced ". and." with " .and."
C>> 1996-04-15 SWCOMP Krogh   Changed single common block name to SWCOM.
C>> 1995-11-21 SWCOMP Krogh   Removed multiple entries for C conversion.
C>> 1994-11-11 SWCOMP Krogh   Declared all vars.
C>> 1994-10-31 SWCOMP Krogh  Changes to use M77CON
C>> 1987-12-07 SWCOMP Lawson  Initial code.
c--S replaces "?": ?WCOMP,?WATAN,?WASIN,?WACOS,?WATN2,?WSUM,?WDIF,
c--&   ?WSQRT,?WEXP,?WSIN,?WCOS,?WTAN,?WSINH,?WCOSH,?WTANH,?WSET,?WSUM1,
c--&   ?WDIF1,?WPRO1,?WQUO1,?WLOG,?WPWRI,?WCHN,?WPRO,?WQUO,?PASCL,
c--&   ?WRCHN,?WACSI,?WCOM
c
C     File SWCOMP [- Wengert derivative COMPutation]
c     The file SWCOMP contains a set of program units to perform
c     computation on (N+1)-tuples representing the value of a quantity
c     and its first N derivatives with respect to a single variable.
c     In the comments we use T as the name of the ultimate independent
c     variable.
c          The parameter NMAX must be set in program units
c                   SWATAN, SWSUM, and SWPRO
c     to establish the largest order of derivative the package
c     will be able to handle.
c     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c     C. L. Lawson, JPL, 1971.
C     References:
c     1.   Wengert, R. E., A simple automatic derivative evaluation
c          program, Comm. ACM, 1, Aug 1964, 463-464.
c     2.   C. L. Lawson, Computing Derivatives using W-arithmetic and
c          U-arithmetic., JPL Appl Math TM 289, Sept 1971.
c     Revised by CLL for Fortran 77 in Jan 1987.  Deleted subr WSTART.
c     Put a first time flag in WPRO.  Now user does not need to make
c     any initialization call.
c     9/18/87, CLL. Added new entry names SWASIN & SWACOS.
c     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C     NMAX        Setting NMAX in the parameter statement determines
c           the highest order derivative this program unit will be able
c           to handle.  The following arrays have dimensions depending
c           on NMAX:  S1(), S2(), S3(), S4(), C()
c
c     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c                         ERROR HANDLING
c
c     All detected error conditions are essentially fatal for the
c     requested operation.  We do not stop; rather, we issue an error
c     message and return.  By use of the the MATH77 library subroutine,
c     ERMSET, the user can change the error processing action to cause a
c     STOP following the error message.
c
c     Error  Called
c      No.    name   Explanation
c
c       1    SWASIN  Infinite deriv when arg = -1 or +1
c       1    SWACOS  Infinite deriv when arg = -1 or +1
c       2    SWSQRT  Infinite deriv when arg = 0
c       3    SWQUO1  Zero divisor
c       4    SWPWRI  U**M is infinite when U = 0 and M < 0
c       5    SWPRO   Require dimension NMAX .ge. NDERIV
c       6    SWQUO   Require dimension NMAX .ge. NDERIV
c       7    SWQUO   Zero divisor
c       8    SPASCL  Require N .ge. 0
c       9    SWRCHN  Require U(2) .ne. 0.
c     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c               PROGRAM UNITS, entry NAMES, and calls
c
c     All of these subroutine and entry names may be called directly by
c     a user, however, it is expected that SWCHN, SWRCHN, and SPASCL
c     would only be used to augment the package with new functions.
c
C     Subroutine names:                       Calls to:
c
c     subroutine SWATAN  (N5,U5,Y5)           SWPWRI,SWCHN
c     subroutine SWASIN (N16,U16,Y16)         SWACSI
c     subroutine SWACOS (N16,U16,Y16)         SWACSI
c     subroutine SWACSI (N16,U16,Y16,ACOSIN)  SWSQRT,SWPWRI,SWPRO1,SWCHN
c     subroutine SWATN2 (N9,Y9,X9,A9)      SWSUM,SWDIF,SWPRO,SWQUO,SWCHN
c
c     subroutine SWSUM  (NDERIV,U,V,W)        None
c     subroutine SWDIF  (N2,U2,V2,W2)         None
c     subroutine SWSQRT (N15,U15,Y15)         SWCHN
c     subroutine SWEXP  (N4,U4,Y4)            SWCHN
c     subroutine SWSIN  (N7,U7,Y7)            SWCHN
c     subroutine SWCOS  (N7,U7,Y7)            SWCHN
c     subroutine SWTAN  (N7,U7,Y7)            SWCHN
c     subroutine SWSINH (N7,U7,Y7)            SWCHN
c     subroutine SWCOSH (N7,U7,Y7)            SWCHN
c     subroutine SWTANH (N7,U7,Y7)            SWCHN, SWPRO, SWQUO
c     subroutine SWSET  (N10,UVAL,UDER,U10)   None
c     subroutine SWSUM1 (N11,C11,U11,Y11)     None
c     subroutine SWDIF1 (N12,C12,U12,Y12)     None
c     subroutine SWPRO1 (N13,C13,U13,Y13)     None
c     subroutine SWQUO1 (N14,C14,U14,Y14)     SWQUO
c
c     subroutine SWLOG  (NDERIV,U,Y)          SWCHN
c     subroutine SWPWRI (NDERIV,MPWR,U,Y)     SWCHN
c
c     subroutine SWCHN  (NDERIV,U,F)          SWPRO
c
c     subroutine SWPRO  (NDERIV,U,V,W)        SPASCL
c     subroutine SWQUO  (NDERIV,U,V,W)        SPASCL
c
c     subroutine SPASCL (N   ,C)              None
c     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c                   Dependencies between program units
c
c     The progran units are SWATAN,SWSUM,SWLOG,SWCHN,SWPRO,SPASCL
c     SWATAN has calls into        SWSUM,SWLOG,SWCHN,SWPRO
c     SWSUM has calls into               SWLOG,SWCHN,SWPRO
c     SWLOG has calls into                     SWCHN
c     SWCHN has calls into                           SWPRO
c     SWPRO has calls into                                 SPASCL
c     SPASCL has no calls
c     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c
      integer N5
      real             U5(*),Y5(*)
      integer NMAX, NMAXP1
      parameter (NMAX=10, NMAXP1=NMAX+1)
      integer I, NS
      real             S1(NMAXP1)
c     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C
C     COMPUTE..              Y5=ATAN(U5)    and derivs w.r.t. T
C
      Y5(1)= ATAN (U5(1))
      IF(N5 .eq. 0) return
C
C     deriv of Y w.r.t. U IS  1./(1. + U**2)
C     START BY CONSTRUCTING 1.+U**2 AND ITS derivs IN S1().
C
      S1(1)=1.0E0 + U5(1)**2
      NS=N5-1
      IF(NS .eq. 0) go to 52
      S1(2)= U5(1)+U5(1)
      IF(NS .eq. 1) go to 52
      S1(3)= 2.0E0
      IF(NS .eq. 2) go to 52
          DO 50 I=3,NS
   50     S1(I+1) = 0.0E0
C
C     NOW S1() = 1.+U**2 and derivs w.r.t. U .
C     COMPUTE 1./S1 and derivs w.r.t. U.  STORE RESULT STARTING IN Y(2)
C
   52 call SWPWRI (NS,-1,S1,Y5(2))
C
C     Convert Y from derivs w.r.t. U TO derivs w.r.t. T
      call SWCHN(N5,U5,Y5)
      return
      end
C     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      subroutine SWASIN(N16,U16,Y16)
C
C     COMPUTES..        Y = ASIN(U) and derivs w.r.t. T
C
      integer N16
      real             U16(*), Y16(*)
c
      Y16(1)= ASIN (U16(1))
      if (N16 .ne. 0) call SWACSI(N16,U16,Y16,.FALSE.)
      return
      end
C     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      subroutine SWACOS(N16,U16,Y16)
C
C     COMPUTES..        Y = ACOS(U) and derivs w.r.t. T
C
      integer N16
      real             U16(*), Y16(*)
c
      Y16(1)= ACOS (U16(1))
      if (N16 .ne. 0) call SWACSI(N16,U16,Y16,.TRUE.)
      return
      end
C     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      subroutine SWACSI(N16,U16,Y16,ACOSIN)
c      Common code for SWACOS and SWASIN  -- ACOSIN - .true. for SWACOS
      integer N16
      real             U16(*), Y16(*)
      logical ACOSIN
c
      integer I, NS
      integer NMAX, NMAXP1
      parameter (NMAX=10, NMAXP1=NMAX+1)
      real             S1(NMAXP1),S2(NMAXP1)
c
C
C     Deriv of Y w.r.t. U is  (+,-)1.0/sqrt(1.0 - U**2)
C     Start by constructing 1.0 - U**2 and its derivs in S1().
C
      S1(1)=1.0E0 - U16(1)**2
      if(S1(1) .eq. 0.0E0) then
c                                     Error condition.
         if(ACOSIN) then
            call ERMSG('SWACOS',1,0,
     *      'Deriv of ACOS(x) is infinite at x = -1 or +1','.')
         else
            call ERMSG('SWASIN',1,0,
     *      'Deriv of ASIN(x) is infinite at x = -1 or +1','.')
         endif
         return
      endif
c
      NS=N16-1
      IF(NS .eq. 0) go to 66
      S1(2)= -2.0E0 * U16(1)
      IF(NS .eq. 1) go to 66
      S1(3)= -2.0E0
      IF(NS .eq. 2) go to 66
          DO 64 I=3,NS
   64     S1(I+1) = 0.0E0
C
C     Now S1() = 1.0 - U**2 and derivs w.r.t. U.
c     Compute S2() = sqrt(1.0 - U**2) and derivs w.r.t. U.
C     Then change sign if doing Arc Cosine.
c
   66 continue
      call SWSQRT(NS,S1,S2)
      if(ACOSIN) call SWPRO1(NS,-1.0E0,S2,S2)
c
C     Compute 1.0/S2 and derivs w.r.t. U.  Store result starting in Y(2)
c
      call SWPWRI (NS,-1,S2,Y16(2))
C
C     Convert Y from derivs w.r.t. U to derivs w.r.t. T
      call SWCHN(N16,U16,Y16)
      return
      end
C     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      subroutine SWATN2(N9,Y9,X9,A9)
C
C     COMPUTES..        A = ATAN2(Y,X) = ATAN(Y/X) and derivs w.r.t. T
C
      integer N9
      real             Y9(*),X9(*),A9(*)
      integer I, NM1, NP1
      real             BIG, XX, YY
      integer NMAX, NMAXP1
      parameter (NMAX=10, NMAXP1=NMAX+1)
      real             S1(NMAXP1),S2(NMAXP1),S3(NMAXP1),S4(NMAXP1)
c
      A9(1)= ATAN2(Y9(1),X9(1))
      IF(N9 - 1) 90,92,94
   90 return
C                            SPECIAL FOR N9 = 1
   92 continue
      BIG= max( ABS(X9(1)), ABS(Y9(1)))
      YY  =Y9(1)/BIG
      XX    =X9(1)/BIG
      A9(2)=(XX*(Y9(2)/BIG) - YY*(X9(2)/BIG) )/(XX**2 + YY**2)
      return
C                            GENERAL CASE    N9 .ge. 2
   94 continue
      NP1=N9+1
      BIG= max( ABS(X9(1)), ABS(Y9(1)))
          DO 96 I=1,NP1
          S1(I)=X9(I)/BIG
   96     S2(I)= Y9(I)/BIG
      NM1=N9-1
      call SWPRO (NM1,S1,S2(2),S3)
      call SWPRO (NM1,S2,S1(2),S4)
      call SWDIF  (NM1,S3,S4,S3)
      call SWPRO (NM1,S1  ,S1  ,S1  )
      call SWPRO (NM1,S2,S2,S2)
      call SWSUM  (NM1,S1  ,S2,S1 )
      call SWQUO (NM1,S3,S1  ,A9(2))
      return
      end

      subroutine SWSUM  (NDERIV,U,V,W)
C
C     COMPUTE..              W=U+V      and derivs w.r.t. T
C
      integer NDERIV
      real             U(*),V(*),W(*)
      integer I, NP1
c
      NP1=NDERIV+1
          DO 10 I=1,NP1
   10     W(I)=U(I)+V(I)
      return
      end
C
C
C
      subroutine SWDIF  (N2,U2,V2,W2)
C
C     COMPUTE..              W2=U2-V2   and derivs w.r.t. T
      integer N2
      real             U2(*),V2(*),W2(*)
      integer I, NP1
c
      NP1= N2+1
          DO 20 I=1,NP1
   20     W2(I)=U2(I)-V2(I)
      return
      end
C
C
C
      subroutine SWSQRT (N15,U15,Y15)
C
C     COMPUTE..        Y= SQRT(U)      and derivs w.r.t. T
C
      integer N15
      real             U15(*),Y15(*)
      integer I
      real              FAC, SCALE, SQRTU
      integer NMAX, NMAXP1
      parameter (NMAX=10, NMAXP1=NMAX+1)
      real             S1(NMAXP1)
c
      SQRTU = SQRT(U15(1))
      if(N15 .eq. 0) go to 152
      if(SQRTU .eq. 0.0E0) then
         call ERMSG('SWSQRT',2,0,
     *    'Deriv of sqrt(x) is infinite at x = 0','.')
         do 155 I = 2,N15+1
            Y15(I) = 0.0E0
  155    continue
         return
      endif
      if(N15 .gt. 1) go to 154
         Y15(2)=U15(2)/(SQRTU+SQRTU)
  152    Y15(1)=SQRTU
         return
C                            GENERAL CASE FOR N .ge. 2
  154 continue
      SCALE=1.0E0/U15(1)
      Y15(1)=1.0E0
      S1(1)  =1.0E0
      FAC=0.5E0
          DO 156 I=1,N15
          Y15(I+1)=Y15(I)*FAC
          FAC=FAC-1.0E0
  156     S1(I+1)=U15(I+1)*SCALE
      call SWCHN(N15,S1,Y15)
      Y15(1)=SQRTU
          DO 158 I=1,N15
  158     Y15(I+1)=Y15(I+1)*SQRTU
      return
      end
C     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C
C
      subroutine SWEXP  (N4,U4,Y4)
C
C     COMPUTE..              Y4=EXP(U4)     and derivs w.r.t. T
C
      integer N4
      real             U4(*),Y4(*)
      integer I
C
      Y4(1)=EXP(U4(1))
      IF(N4 .eq. 0) return
          DO 40 I=1,N4
   40     Y4(I+1)=Y4(I)
      call SWCHN(N4,U4,Y4)
      return
      end
C
C
C     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      subroutine SWSIN(N7,U7,Y7)
C
C     COMPUTE..              Y7=SIN(U7) and derivs w.r.t. T
C
      integer N7
      real             U7(*),Y7(*)
      integer I
C
      Y7(1) = SIN(U7(1))
      IF(N7 .eq. 0) return
      Y7(2) = COS(U7(1))
      if(N7 .eq.1) then
         Y7(2) = U7(2) * Y7(2)
      else
         DO 70 I=2,N7
            Y7(I+1)=-Y7(I-1)
   70    continue
         call SWCHN(N7,U7,Y7)
      end if
      return
      end
C
C
C
      subroutine SWCOS  (N7,U7,Y7)
C
C     COMPUTE..              Y7=COS(U7)      and derivs w.r.t. T
C
      integer N7
      real             U7(*),Y7(*)
      integer I
C
      Y7(1)= COS (U7(1))
      IF(N7 .eq. 0) return
      Y7(2)= -SIN(U7(1))
      if(N7 .eq.1) then
         Y7(2) = U7(2) * Y7(2)
      else
         DO 70 I=2,N7
            Y7(I+1)=-Y7(I-1)
   70    continue
         call SWCHN(N7,U7,Y7)
      end if
      return
      end
c     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      subroutine SWTAN  (N7,U7,Y7)
C
C     COMPUTE..              Y7=TAN(U7)      and derivs w.r.t. T
C     The first deriv is 1/cos(u)**2
c
      integer N7
      real             U7(*),Y7(*)
      integer I, N7M1
      real             CS, SN
      integer NMAX, NMAXP1
      parameter (NMAX=10, NMAXP1=NMAX+1)
      real             S1(NMAXP1)
c
      if(N7 .eq. 0) then
         Y7(1) = tan(U7(1))
         return
      end if
      SN = SIN(U7(1))
      CS = COS(U7(1))
      Y7(1)= SN/CS
c
c           Compute Cos(U) and its derivs w.r.t. U in S1(2..N+1)
c
      S1(2) = CS
      S1(3) = -SN
          DO 76 I=4,N7+1
   76     S1(I) = -S1(I-2)
      N7M1 = N7-1
c                                       Convert to derivs w.r.t. t
      call SWCHN(N7M1,U7,S1(2))
c                               Compute Cos(U)**2 & derivs w.r.t. t
      call SWPRO(N7M1, S1(2), S1(2), S1(2))
c                          Divide first deriv of U by Cos(U)**2,
c                          both with derivs w.r.t. t
      call SWQUO(N7M1, U7(2), S1(2), Y7(2))
      return
      end
c     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      subroutine SWSINH(N7,U7,Y7)
C
C     COMPUTE..              Y7=SINH(U7) and derivs w.r.t. T
C
      integer N7
      real             U7(*),Y7(*)
      integer I
C
      Y7(1) = SINH(U7(1))
      IF(N7 .eq. 0) return
      Y7(2) = COSH(U7(1))
      if(N7 .eq.1) then
         Y7(2) = Y7(2) * U7(2)
      else
         DO 80 I=2,N7
            Y7(I+1) = Y7(I-1)
   80    CONTINUE
         call SWCHN(N7,U7,Y7)
      end if
      return
      end
C
C
C
      subroutine SWCOSH  (N7,U7,Y7)
C
C     COMPUTE..              Y7=COSH(U7)      and derivs w.r.t. T
C
      integer N7
      real             U7(*),Y7(*)
      integer I
C
      Y7(1)= COSH (U7(1))
      IF(N7 .eq. 0) return
      Y7(2)= SINH(U7(1))
      if(N7 .eq.1) then
         Y7(2) = Y7(2) * U7(2)
      else
         DO 80 I=2,N7
            Y7(I+1) = Y7(I-1)
   80    CONTINUE
         call SWCHN(N7,U7,Y7)
      end if
      return
      end
C
C
C
      subroutine SWTANH  (N7,U7,Y7)
C
C     COMPUTE..              Y7=TANH(U7)      and derivs w.r.t. T
C     The first deriv is 1/cosh(u)**2
c
      integer N7
      real             U7(*),Y7(*)
      integer I, N7M1
      real             CH, SH
      integer NMAX, NMAXP1
      parameter (NMAX=10, NMAXP1=NMAX+1)
      real             S1(NMAXP1)
C
      if(N7 .eq. 0) then
         Y7(1) = tanh(U7(1))
         return
      endif
      SH = SINH(U7(1))
      CH = COSH(U7(1))
      Y7(1)= SH/CH
c
c           Compute Cosh(U) and its derivs w.r.t. U in S1(2..N+1)
c
      S1(2) = CH
      S1(3) = SH
          DO 94 I=4,N7+1
   94     S1(I) = S1(I-2)
      N7M1 = N7-1
c                                       Convert to derivs w.r.t. t
      call SWCHN(N7M1,U7,S1(2))
c                               Compute Cosh(U)**2 & derivs w.r.t. t
      call SWPRO(N7M1, S1(2), S1(2), S1(2))
c                          Divide first deriv of U by Cosh(U)**2,
c                          both with derivs w.r.t. t
      call SWQUO(N7M1, U7(2), S1(2), Y7(2))
      return
      end
C     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      subroutine SWSET  (N10,UVAL,UDER,U10)
C
C     INITIALIZE THE ARRAY U10()..  SET VALUE=UVAL, 1-ST DERIV = UDER,
C                                   AND HIGHER derivs = ZERO
      integer N10
      real             U10(*), UVAL, UDER
      integer I
C
      IF(N10-1) 104,102,100
  100     DO 101 I=2, N10
  101     U10(I+1)= 0.0E0
  102 U10(2)=UDER
  104 U10(1)=UVAL
      return
      end
C
C
C
      subroutine SWSUM1 (N11,C11,U11,Y11)
C
C     COMPUTE..        Y=C+U      WHERE C IS A SCALAR
C
      integer N11
      real             C11, U11(*), Y11(*)
      integer I
C
      Y11(1)=C11+U11(1)
          DO 110 I=1,N11
  110     Y11(I+1)=U11(I+1)
      return
      end
C
C
C
      subroutine SWDIF1 (N12,C12,U12,Y12)
C
C     COMPUTE..        Y=C-U      WHERE C IS A SCALAR
C
      integer N12
      real             C12, U12(*), Y12(*)
      integer I
C
      Y12(1)=C12-U12(1)
          DO 120 I=1,N12
  120     Y12(I+1)=-U12(I+1)
      return
      end
C
C
C
      subroutine SWPRO1(N13,C13,U13,Y13)
C
C     COMPUTE..        Y=C*U      WHERE C IS A SCALAR
C
      integer N13
      real             C13, U13(*), Y13(*)
      integer I
C
          DO 130 I=1,N13+1
  130     Y13(I) = C13*U13(I)
      return
      end
C
C
C
      subroutine SWQUO1(N14,C14,U14,Y14)
C
C     COMPUTE..        Y=C/U      WHERE C IS A SCALAR
C
      integer N14
      real             C14, U14(*), Y14(*)
      integer I
      integer NMAX, NMAXP1
      parameter (NMAX=10, NMAXP1=NMAX+1)
      real             S1(NMAXP1)
C
      if(U14(1) .eq. 0.0E0) then
         call ERMSG('SWQUO1',3,0, 'Zero divisor.','.')
         return
      endif
      S1(1) = C14
      do 160 I = 2, N14+1
         S1(I) = 0.0E0
  160 continue
      call SWQUO (N14,S1,U14,Y14)
      return
      end

      subroutine SWLOG  (NDERIV,U,Y)
C     C. L. LAWSON, JPL, 1971 NOV 19
C
C     COMPUTE..              Y= LOG BASE E OF U   and derivs w.r.t. T
C     ERROR STOP IN this subr IF U(1) .le. 0.
C
      integer         II, NDERIV
      real                         U(*),Y(*)
      real             FAC, UINV
C
      Y(1)= LOG(U(1))
      if (NDERIV .le. 1) then
         if (NDERIV .eq. 1) Y(2)= U(2)/U(1)
         return
      end if
      UINV=1.0E0/U(1)
      Y(2)=UINV
      FAC=-UINV
      DO 120 II = 2, NDERIV
         Y(II+1)=Y(II)*FAC
         FAC=FAC-UINV
  120 continue
C
C     Y() NOW CONTAINS VALUE and derivs w.r.t. U
C     USE CHAIN RULE TO convert TO derivs w.r.t. T
C
      call SWCHN(NDERIV,U,Y)
      return
      end
C
      subroutine SWPWRI (NDERIV,MPWR,U,Y)
C
C     COMPUTE..              Y = U**MPWR  and derivs w.r.t. T
C     MPWR is an integer, POS.,NEG.,OR ZERO.
C     MPWR does not depend ON T.
C     U may depend on T
C     Issues an error message if both (U(1) .eq. 0.) and (MPWR .lt. 0)
C     If MPWR .eq. 0 then set Y(1) = 1. and all derivs = 0.
C
      integer         I, II, II1, II2, J, M, MPWR, NDERIV
      real                         U(*),Y(*)
      logical          SHORT
      real             D, FAC, UINV
C
      II2=NDERIV
      SHORT = .FALSE.
      M=MPWR
      IF(M     ) 110,35,45
C                                      HERE MPWR .eq. 0
   35 Y(1)= 1.0E0
      IF(NDERIV .le. 0) return
          DO 40 J=1,NDERIV
   40     Y(J+1) = 0.0E0
      return
C                                      HERE MPWR IS POSITIVE
   45 IF(U(1) .ne. 0.0E0) go to 50
      Y(1)=0.0E0
      IF(NDERIV .eq. 0) return
      IF(M .le. NDERIV) go to 46
      II=1
      go to 130
C
   46     DO 47 I=1,NDERIV
   47     Y(I+1) = 0.0E0
      FAC=1.0E0
      D  =0.0E0
          DO 48 I=1,M
          D=D+1.0E0
   48     FAC=FAC*D
      Y(M+1) = FAC
      go to 150
C
   50 IF( M .ge. NDERIV) go to 112
      II2=M
      SHORT = .TRUE.
      go to 112
C
C                                      HERE MPWR .lt 0
  110 continue
      if(U(1) .eq. 0.0E0) then
            call IERM1('SWPWRI',4,0,
     *      'U**M is infinite when U = 0. and M < 0',
     *      'M',M,'.')
         do 111 I = 1,NDERIV+1
            Y(I) = 0.0E0
  111    continue
         return
      endif
  112 continue
      Y(1)=U(1)**M
      if(NDERIV .eq. 0) return
      UINV=1.0E0/U(1)
      FAC=UINV*M
      II1=1
          DO 120 II=II1,II2
          Y(II+1)=Y(II)*FAC
          FAC=FAC-UINV
  120     continue
      if(.NOT. SHORT) go to 150
      II=II2+1
C
C                            SET HIGHER derivs TO ZERO.
C
  130 continue
          DO 140 J=II,NDERIV
  140     Y(J+1)=0.0E0
C
C     Y() NOW CONTAINS VALUE and derivs w.r.t. U
C     USE CHAIN RULE TO convert TO derivs w.r.t. T
C
  150 call SWCHN(NDERIV,U,Y)
      return
      end

      subroutine SWCHN(NDERIV,U,Y)
C
c     Implements the chain rule.
c     Given values of y(u) and its derivs with respect to u in Y(), and
c     u(x) and its derivs with respect to x in U(), this subr replaces
c     the contents of Y() with y(u(x)) and its derivs with respect to x.
c
C     NDERIV                 HIGHEST ORDER DERIVATIVE
c
C     (U(I),I=1,NDERIV+1)    INPUT.. value of u and derivs w.r.t. x
C
C     (Y(I),I=1,NDERIV+1)    INPUT.. value of y and derivs w.r.t. u
C                            OUTPUT.. (Y(I+1),I=1,NDERIV) will be
C                                     replaced by derivs w.r.t. x
      integer        I, L, NDERIV, NP1
      real                 U(*), Y(*)
C
      if(NDERIV .eq. 0) return
      if(U(2) .ne. 1.) go to 20
C                             U(2) .eq. 1.,  See if higher derivs are 0.
      if(NDERIV .eq. 1) return
          DO 10 I=2,NDERIV
          if(U(I+1) .ne. 0.) go to 50
   10     continue
      return
C                                   Test if all derivs of U are 0.
   20     DO 30 I=1,NDERIV
          if(U(I+1) .ne. 0.) go to 50
   30     continue
C
          DO 40 I=1,NDERIV
   40     Y(I+1)=0.
      return
C
C         Here for the general case when U is not T and not constant.
C
   50     NP1=NDERIV+1
          DO 60 L=0,NDERIV-1
   60     call SWPRO (L, Y(NP1-L), U(2), Y(NP1-L))
      return
      end

      subroutine SWRCHN(NDERIV, U, Y)
C
c     Implements the "reverse" chain rule.
c     Given values of y(u(x)) and its derivs with respect to x in Y(),
c     and u(x) and its derivs with respect to x in U(), this subr replac
c     the contents of Y() with y(u) and its derivs with respect to u.
c     Requires u'(x) .ne. 0.
c
c     Note that this subr can be used to compute a representation of
c     the inverse function of a given function.  This transformation
c     is called "reversion" of a power series.  To do this let y() be
c     the inverse function of u() in a neighborhood of a point, x0.
c     Thus y(u(x)) = x for all x in a neighborhood of x0.
c     Then the value of y() and its derivs with respect to x, evaluated
c     at x0 are (x0, 1.0, 0.0, ..., 0.0).
c     If Y() contains these values on entry, and U() contains u(x) and
c     its derivs w.r.t. x, evaluated at x0, then on return Y() will
c     contain y(u) and its derivs w.r.t. u, evaluated at u0 = u(x0).
c     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C     NDERIV  [in]  Highest order derivative to be considered.
C
C     (U(I),I=1,NDERIV+1)  [in]  Contains the value of u(x) and its
c           first NDERIV derivatives with respect to x, evaluated at x0.
c           Require U(2) .ne. 0.
c
C     (Y(I),I=1,NDERIV+1)  [inout]  On entry, Y() must contain y(u(x))
c           and its first NDERIV derivatives w.r.t. x, evaluated at x0.
c           On return, Y() will contain y(u) and its first NDERIV
c           derivatives w.r.t. u, evaluated at u0 = u(x0).
c     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      integer NDERIV, NP1, L
      real              U(*), Y(*)
c     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      NP1=NDERIV+1
c
      if( U(2) .eq. 0.0E0 .and. NP1 .ge. 2) then
         call ERMSG('SWRCHN',9,0,'U(2) is zero.','.')
         return
      endif
c
      do 60 L=2, NP1
         CALL SWQUO (NP1 - L, Y(L), U(2), Y(L))
   60 continue
      return
      end

      subroutine SWPRO(NDERIV,U,V,W)
C                      COMPUTE W(1)=U(1)*V(1) AND DERIVATIVES .
C     NDERIV           HIGHEST ORDER DERIVATIVE TO BE COMPUTED
C                      NDERIV .ge. 0
C     (U(I),I=1,NDERIV+1)    INPUT ARRAY..
C                                       U(1)= VALUE
C                                       U(I+1)= I-TH DERIVATIVE OF U(1)
C     (V(I),I=1,NDERIV+1)    INPUT ARRAY..      SAME FORMAT AS U().
C
C     (W(I),I=1,NDERIV+1)    OUTPUT ARRAY..
C                                            W(1)= U(1)*V(1)
C                                            W(I+1)=Ith deriv of W(1)
c     It is permissible for W() to occupy the same storage locations as
c     U() and/or V().
c     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C     NMAX        Setting NMAX in the parameter statement determines
c           the highest order derivative this program unit will be able
c           to handle.  The dimension of the internal saved array C()
c           will be set as a function of NMAX.
C     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      integer I, ID2, IMID, IP1, J, JC, NDERIV, NP1
      real                         U(*),V(*),W(*)
      integer NMAX, NCOEF
      parameter(NMAX = 10, NCOEF = 1 + (NMAX*(NMAX+1))/2 )
      real                C(NCOEF)
      real                TEMP
      logical FIRST
      common / SWCOM / C
      save FIRST, /SWCOM/
      data FIRST /.TRUE. /
C     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      if(FIRST) then
         call SPASCL(NMAX, C)
         FIRST = .FALSE.
      endif
      if( NDERIV .gt. NMAX) then
         call IERM1('SWPRO',5,0,
     *   'Dimension NMAX in SWATAN, SWSUM, & SWPRO must be .ge. NDERIV',
     *   'NMAX',NMAX,',')
         call IERV1('NDERIV',NDERIV,'.')
         return
      endif
      NP1=NDERIV+1
c     JC = (NDERIV*(NDERIV-1))/2
c     DO 20 I=NP1,1,-1
c         TEMP=0.0E0
c         IP1=I+1
c              DO 10 J=1,I
c              TEMP = TEMP + U(J)*V(IP1-J)*C(JC+J)
c  10     continue
c         W(I) = TEMP
c         JC=JC- I+2
c  20 continue
c
c     The following code does the same as the above commented-out
c     code but reduces the number of multiplies.
c     Treats the cases of NDERIV = 0, 1, 2, 3, and 4 special for
c     efficiency, since it is expected that most usage will involve
c     small values of NDERIV.  In particular the chain rule subroutine,
c     _WCHN, will call this subr with a sequence of values of NDERIV
c     going down to 0.
c     Treats NDERIV > 4 in a general way, but takes account of
c     the fact that the C() terms in the "do 10" loop above are
c     symmetric about the middle term and the first and last term
c     are each = 1.
c
      go to (41,42,43,44,45), NP1
c
c     Fortran 77 specifies that a Computed Go To drops through is the
c     index is out of range.  Thus we arrive here if NP1 < 1 or > 5.
c     NP1 < 1 would be an error.  We choose not to use time to test for
c     it however.
c     NP1 > 5 is valid. In that case we do the following loop for I
c     running down from NP1 to 6, and then do the special statements
c     for 5, 4, 3, 2, and 1.
c
      JC = (NDERIV*(NDERIV-1))/2
      DO 20 I=NP1,6,-1
          TEMP = U(1) * V(I) + U(I) * V(1)
          IP1=I+1
          ID2 = I/2
          DO 10 J=2,ID2
               TEMP = TEMP + C(JC+J) * (U(J)*V(IP1-J) + U(IP1-J)*V(J))
   10     continue
          if(2*ID2 .eq. I) then
             W(I) = TEMP
          else
             IMID = ID2 + 1
             W(I) = TEMP + C(JC+IMID) * U(IMID) * V(IMID)
          endif
          JC=JC- I+2
   20 continue
c
   45 W(5) = U(1)*V(5) + 4.0E0*(U(2)*V(4) + U(4)*V(2)) +
     *                   6.0E0*U(3)*V(3) + U(5)*V(1)
   44 W(4) = U(1) * V(4) + 3.0E0*(U(2)*V(3) + U(3) * V(2)) + U(4)*V(1)
   43 W(3) = U(1) * V(3) + 2.0E0 * U(2) * V(2) + U(3) * V(1)
   42 W(2) = U(1) * V(2) + U(2) * V(1)
   41 W(1) = U(1) * V(1)
      return
      end
c     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      subroutine SWQUO (NDERIV, U ,V ,W )
C
C     COMPUTE..          W = U / V  and derivs w.r.t. T
c
c     V() must be distinct in storage from W(), but U() can be the same
c     array as W().
c
C     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      integer NDERIV
      real             U(*),V(*),W(*)
      integer I, ID2, IMID, IP1, J, JC, NP1
      integer NMAX, NCOEF
      parameter(NMAX = 10, NCOEF = 1 + (NMAX*(NMAX+1))/2 )
      real                C(NCOEF)
      real                FAC, TEMP
      logical FIRST
      common / SWCOM / C
      save FIRST, /SWCOM/
      data FIRST /.TRUE. /
C     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      IF(FIRST) THEN
         CALL SPASCL(NMAX, C)
         FIRST = .FALSE.
      END IF
      if( NDERIV .GT. NMAX) then
         call IERM1('SWQUO',6,0,
     *   'Dimension NMAX in SWATAN, SWSUM, & SWPRO must be .ge. NDERIV',
     *   'NMAX',NMAX,',')
         call IERV1('NDERIV',NDERIV,'.')
         return
      endif
c
      if( V(1) .eq. 0.0E0) then
         call ERMSG('SWQUO',7,0,'Zero divisor.','.')
         return
      endif
c
      FAC = 1.0E0/V(1)
      W(1) = U(1) * FAC
      NP1=NDERIV+1
      JC = 0
c         DO 40 I=2,NP1
c            TEMP = U(I)
c            IP1=I+1
c            DO 30 J=2,I
c               TEMP = TEMP - V(J)*W(IP1-J)*C(JC+J)
c  30        CONTINUE
c            W(I) = TEMP * FAC
c            JC=JC + I-1
c  40     CONTINUE
      DO 40 I= 2, NP1
          TEMP = U(I) - W(1) * V(I)
          IP1=I+1
          ID2 = I/2
          DO 30 J=2,ID2
               TEMP = TEMP - C(JC+J) * (W(J)*V(IP1-J) + W(IP1-J)*V(J))
   30     continue
          if(2*ID2 .ne. I) then
             IMID = ID2 + 1
             TEMP = TEMP - C(JC+IMID) * W(IMID) * V(IMID)
          endif
          W(I) = TEMP * FAC
          JC=JC + I-1
   40 continue
      return
      end

      subroutine SPASCL(N   ,C)
C        C.L.LAWSON,JPL, 1969 DEC 10
C
C     N                      MAXIMUM ORDER DERIVATIVE TO BE
C                            COMPUTED USING W-ARITHMETIC
C
C     (C(I),I=1,NC)          FIRST N+1 ROWS OF PASCAL TRIANGLE PACKED
C                            WITH 1-ST DIAGONAL OMITTED AFTER 1-ST ROW.
C                            NC= (N*(N+1)/2) + 1
C     EXAMPLE..    IF N=3 then TRIANGLE IS       1
C                                              1    1
C                                            1   2    1
C                                          1   3    3   1
C
C     WHICH WILL BE STORED AS.. 1,1,2,1,3,3,1
C
      integer          I, J, K, L, N
      logical          DONE
      real             C(*)
      save             DONE
      data DONE / .false. /
C     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      if (DONE) return
      DONE = .true.
      if(N .lt. 0) then
            call IERM1('SPASCL',8,0, 'Require N .ge. 0','N',N,'.')
            return
      endif
      C(1)=1.0E0
      if(N .eq. 0) return
      C(2)=1.0E0
      if(N .eq. 1) return
      K=3
      DO 520 I=1,N-1
          L=I+1
          DO 510 J=1,I
               C(K)=C(K-L)+C(K-L+1)
               K=K+1
  510     continue
          C(K)=1.0E0
          K=K+1
  520 continue
      return
      end