text:
      subroutine jgibbs(aamatr,afcnst,affp,cdrs,cscale,csts,delvec,
     $ eps100,gmmatr,iindx1,iodb,ipivot,ipndx1,jpflag,kbt,kdim,kmax,
     $ km1,kmt,kx1,kxt,mtb,nbasp,nbtmax,ndrs,ndrsmx,ndrsr,nodbmx,
     $ noutpt,npadd,npdel,nptmax,nstmax,nsts,nstsmx,nstsr,nttyo,
     $ rhsvec,uphase,uspec,xlks)
c
c     This subroutine determines which one of a set of phases should be
c     removed from the equilibrium system when there is a mineralogic
c     phase rule violation. It has not yet been extended to handle
c     solid solutions.
c
c     It is assumed heres that a mineralogical phase rule violation
c     does in fact exist. This subroutine does not make an independent
c     confirmation.
c
c     This subroutine is called by:
c
c       EQ6/eqcalc.f.
c
c-----------------------------------------------------------------------
c
c     Principal input:
c
c       npadd  = index of the last phase added to the equilibrium
c                  phase assemblage
c
c     Principal output:
c
c       npdel = index of the phase which should be deleted from the
c                 equilibrium phase assemblage in order to resolve
c                 the violation
c
c
c-----------------------------------------------------------------------
c
      implicit none
c
c-----------------------------------------------------------------------
c
c     Calling sequence variable declarations.
c
      integer kmax,nbtmax,ndrsmx,nodbmx,nptmax,nstmax,nstsmx
c
      integer noutpt,nttyo
c
      integer iindx1(kmax),iodb(nodbmx),ipivot(kmax),ipndx1(kmax),
     $ jpflag(nptmax),nbasp(nbtmax),ndrs(ndrsmx),ndrsr(2,nstmax),
     $ nsts(nstsmx),nstsr(2,nstmax)
c
      integer kbt,kdim,km1,kmt,kx1,kxt,npadd,npdel
c
      character*48 uspec(nstmax)
      character*24 uphase(nptmax)
c
      real*8 aamatr(kmax,kmax),affp(nptmax),cdrs(ndrsmx),
     $ cscale(nstmax),csts(nstsmx),delvec(kmax),gmmatr(kmax,kmax),
     $ mtb(nbtmax),rhsvec(kmax),xlks(nstmax)
c
      real*8 afcnst,eps100
c
c-----------------------------------------------------------------------
c
c     Local variable declarations with global dimensions.
c
      integer isv_kmax
c
      SAVE isv_kmax
c
      integer, dimension(:), allocatable :: kmprv,ksprv
c
      SAVE kmprv,ksprv
c
      real(8), dimension(:), allocatable :: cscx
c
      SAVE cscx
c
c-----------------------------------------------------------------------
c
c     Local variable declarations.
c
      integer i,idim,ier,irow1,irow2,j,jcol1,jcol2,jdim,j2,j3,j4,k,
     $ kcol,km,kmm,kp,kpadd,kpp,kpt,n,nb,np,npp,nr1,nr2,ns,nse,nss
c
      integer ilnobl
c
      logical qldep,qpr
c
      real*8 acx,acx2,affc,affmin,afscal,cx,mxx,mxxe,mxxm
c
      real*8 coefdr
c
c-----------------------------------------------------------------------
c
      data qpr    /.false./
c
c-----------------------------------------------------------------------
c
c     Allocate or reallocate local work arrays as needed.
c
      if (.not.ALLOCATED(kmprv)) then
c
c       Local work arrays are not allocated. Zero the saved
c       array size variables. Note that only one array is tested
c       to see if it is allocated. It is assumed that all local
c       work arrays are either allocated or not.
c
        isv_kmax = 0
      else
c
c       Local work arrays are allocated. Check to see if any of the
c       array size variables have changed. If so, deallocate
c       the corresponding local work arrays and zero the corresponding
c       saved size variables.
c
        if (kmax .ne. isv_kmax) then
          DEALLOCATE(kmprv)
          DEALLOCATE(ksprv)
          DEALLOCATE(cscx)
          isv_kmax = 0
        endif
      endif
c
c     At this point, the saved array size values are zero if the
c     corresponding arrays need to be allocated.
c
      if (isv_kmax .eq. 0) then
        ALLOCATE(kmprv(kmax),ksprv(kmax))
        ALLOCATE(cscx(kmax))
        isv_kmax = kmax
      endif
c
c     Zero the contents of the local work arrays.
c
      do k = 1,kmax
        kmprv(k) = 0
        ksprv(k) = 0
      enddo
c
      do k = 1,kmax
        cscx(k) = 0.
      enddo
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
      if (iodb(1) .gt. 0) then
        if (kxt .ge. kx1) then
          write (noutpt,1000)
          write (nttyo,1000)
 1000     format(/' * Note - (EQ6/jgibbs) This subroutine finds the',
     $    ' phases or phases',/7x,'involved in a violation of the',
     $    ' mineralogic phase rule, but it',/7x,'presently considers',
     $    ' only pure minerals in the analysis.')
        endif
      endif
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
c     Clear the Jacobian matrix for use as workspace.
c
      do i = 1,kdim
        do j = 1,kdim
          aamatr(i,j) = 0.
        enddo
      enddo
c
      do i = 1,kbt + 1
        kmprv(i) = 0
      enddo
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
c     Find the phases which are directly involved in the mineralogic
c     phase rule violation. For each possible phase, test to see if
c     the violation disappears when the phase is removed from the
c     tested assemblage.
c
      kpt = 0
      do km = km1,kmt
        ns = iindx1(km)
c
c       For each of the other minerals involved in the current
c       assemblage, set up a corresponding row which contains the
c       reaction coefficients of the basis species.
c
        idim = 0
        do kmm = km1,kmt
          if (kmm .ne. km) then
            nss = iindx1(kmm)
            idim = idim + 1
            do j = 1,kbt
              nb = iindx1(j)
              nse = nbasp(nb)
c
c             Calling sequence substitutions:
c               nss for ns
c
              cx = coefdr(cdrs,ndrs,ndrsmx,ndrsr,nse,nss,nstmax)
              aamatr(idim,j) = cx
            enddo
          endif
        enddo
c
c       Test for linear dependence.
c
        irow1 = 1
        irow2 = idim
        jcol1 = 1
        jcol2 = kbt
        call lindep(aamatr,eps100,irow1,irow2,jcol1,jcol2,kmax,qldep)
c
        if (.not.qldep) then
c
c         Have not obtained linear dependence. The phase being tested
c         is involved in the violation. Note its matrix index in the
c         kmprv array.
c
          kpt = kpt + 1
          kmprv(kpt) = km
        endif
      enddo
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
      if (kpt .ge. 1) then
        write (noutpt,1050)
        write (nttyo,1050)
 1050   format(/' * Note - (EQ6/jgibbs) The phases involved in the',
     $  ' mineralogic',/7x,'phase rule violation are:',/)
c
        do kp = 1,kpt
          km = kmprv(kp)
          ns = iindx1(km)
          np = ipndx1(km)
          j2 = ilnobl(uphase(np))
          write (noutpt,1060) np,uphase(np)(1:j2)
          write (nttyo,1060) np,uphase(np)(1:j2)
 1060     format(2x,i5,2x,a)
        enddo
c
        if (kpt .eq. 1) then
          write (noutpt,1070)
          write (nttyo,1070)
 1070     format(/' * Error (EQ6/jgibbs) Programming error trap:',
     $    ' Have a supposed',/7x,'mineralogical phase rule violation',
     $    ' involving only one phase.')
          stop
        endif
      else
        write (noutpt,1080)
        write (nttyo,1080)
 1080   format(/' * Error (EQ6/jgibbs) Programming error trap:',
     $  ' Have a supposed',/7x,'mineralogical phase rule violation',
     $  ' involving no phases.')
        stop
      endif
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
c     Now determine the affinity of each involved phase under the
c     assumption that the solution is saturated with respect to
c     the others.
c
      do kp = 1,kpt
        km = kmprv(kp)
        ns = iindx1(km)
        np = ipndx1(km)
        if (iodb(1) .gt. 0) then
c
          j2 = ilnobl(uphase(np))
          write (noutpt,1090) uphase(np)(1:j2)
 1090     format(/1x,'Hypothetical affinity calculation for ',a,':')
c
          write (noutpt,1100)
 1100     format(/3x,'Cognate minerals:',/)
          do kpp = 1,kpt
            if (kpp .ne. kp) then
              kmm = kmprv(kpp)
              npp = ipndx1(kmm)
              j3 = ilnobl(uphase(npp))
              write (noutpt,1110) uphase(npp)(1:j3)
 1110         format(5x,a)
            endif
          enddo
        endif
c
c       Once again, clear the Jacobian matrix for use as workspace.
c
        do i = 1,kdim
          do j = 1,kdim
            aamatr(i,j) = 0.
          enddo
        enddo
c
c       For each of the other minerals involved in the violation,
c       set up a corresponding row which contains the reaction
c       coefficients of the basis species. The objective is to
c       find a set of such coefficients which does not exhibit linear
c       dependence. Such dependence is possible because the number of
c       basis species appearing in these reactions may exceed the
c       number of "other" minerals. For example, if three ferric iron
c       minerals are present in the violating assemblage, Fe+++ is
c       a dependent species, and Fe++ is the strict basis species for
c       iron, then the reactions will be written in terms of O2(g), H+,
c       H2O, and Fe++ instead of in terms of Fe+++. Depending on what
c       is exactly in a given set, one or more of O2(g), H+, and H2O
c       may be redundant in terms of relating the current mineral to
c       to the current set of other minerals involved in the violation.
c
        do i = 1,kbt
          ksprv(i) = 0
        enddo
c
        idim = 0
        do kpp = 1,kpt
          if (kpp .ne. kp) then
            kmm = kmprv(kpp)
            nss = iindx1(kmm)
            idim = idim + 1
            do j = 1,kbt
              nb = iindx1(j)
              nse = nbasp(nb)
c
c             Calling sequence substitutions:
c               nss for ns
c
              cx = coefdr(cdrs,ndrs,ndrsmx,ndrsr,nse,nss,nstmax)
              aamatr(idim,j) = cx
            enddo
          endif
        enddo
c
c       Factor the rows.
c
        irow1 = 1
        irow2 = idim
        jcol1 = 1
        jcol2 = kbt
        call lindep(aamatr,eps100,irow1,irow2,jcol1,jcol2,kmax,qldep)
c
        if (qldep) then
          j2 = ilnobl(uphase(np))
          write (noutpt,1120) uphase(np)(1:j2)
          write (nttyo,1120) uphase(np)(1:j2)
 1120     format(/' * Error (EQ6/jgibbs) Programming error trap:',
     $    " Couldn't find a",/7x,'set of independent basis species',
     $    ' reaction coefficients needed',/7x,'to calculate the',
     $    ' hypothetical affinity of ',a,',',/7x,'which is involved',
     $    ' in a violation of the mineralogic phase rule.',/7x,
     $    'Previous testing should guarantee the existence of such',
     $    ' a set.',/7x,'Factoring shows linear dependence.')
          stop
        endif
c
        jdim = 0
        do i = 1,idim
          do kcol = jdim + 1,kbt
            if (abs(aamatr(i,kcol)) .ne. 0.) then
              jdim = jdim + 1
              ksprv(jdim) = kcol
              go to 100
            endif
          enddo
  100     continue
        enddo
c
        if (jdim .ne. idim) then
          j2 = ilnobl(uphase(np))
          write (noutpt,1130) uphase(np)(1:j2),jdim,idim
          write (nttyo,1130) uphase(np)(1:j2),jdim,idim
 1130     format(/' * Error (EQ6/jgibbs) Programming error trap:',
     $    " Couldn't find a",/7x,'set of independent basis species',
     $    ' reaction coefficients needed',/7x,'to calculate the',
     $    ' hypothetical affinity of ',a,',',/7x,'which is involved',
     $    ' in a violation of the mineralogic phase rule.',/7x,
     $    'Previous testing should guarantee the existence of such',
     $    ' a set.',/7x,'Found ',i3,' independent species for ',i3,
     $    ' "other" phases.')
          stop
        endif
c
        if (iodb(1) .gt. 0) then
          write (noutpt,1140)
 1140     format(/3x,'Cognate independent basis species:',/)
          do j = 1,jdim
            kcol = ksprv(j)
            nb = iindx1(kcol)
            nse = nbasp(nb)
            j4 = ilnobl(uspec(nse))
            write (noutpt,1110) uspec(nse)(1:j4)
          enddo
        endif
c
c       Write the matrix for calculating the hypothetical affinity.
c
        j = 0
        do kpp = 1,kpt
          if (kpp .ne. kp) then
            kmm = kmprv(kpp)
            nss = iindx1(kmm)
            j = j + 1
            do i = 1,idim
              kcol = ksprv(i)
              nb = iindx1(kcol)
              nse = nbasp(nb)
c
c             Calling sequence substitutions:
c               nss for ns
c
              cx = coefdr(cdrs,ndrs,ndrsmx,ndrsr,nse,nss,nstmax)
              aamatr(i,j) = cx
            enddo
          endif
        enddo
c
c       Write the right-hand-side vector for calculating the
c       hypothetical affinity.
c
        do j = 1,jdim
          kcol = ksprv(j)
          nb = iindx1(kcol)
          nse = nbasp(nb)
          cx = coefdr(cdrs,ndrs,ndrsmx,ndrsr,nse,ns,nstmax)
          rhsvec(j) = cx
        enddo
c
        if (iodb(1) .gt. 0) then
          write (noutpt,1150)
 1150     format(/3x,'Matrix:',/)
          do i = 1,idim
            write (noutpt,1160) (aamatr(i,j), j = 1,jdim)
 1160       format(5x,5(1pe10.3,2x),e10.3)
          enddo
c
          write (noutpt,1170)
 1170     format(/3x,'RHS vector:',/)
          do j = 1,jdim
            write (noutpt,1180) rhsvec(j)
 1180       format(5x,1pe10.3)
          enddo
        endif
c
c       Calling sequence substitutions:
c         idim for kdim
c
        call msolvr(aamatr,delvec,gmmatr,ier,ipivot,idim,kmax,
     $  noutpt,nttyo,qpr,rhsvec)
c
        if (ier .gt. 0) then
          j2 = ilnobl(uphase(np))
          write (noutpt,1220) uphase(np)(1:j2)
          write (nttyo,1220) uphase(np)(1:j2)
 1220     format(/" * Error (EQ6/jgibbs) Can't resolve a violation of",
     $    ' the mineralogical',/7x,'phase rule involving ',a,' because',
     $    ' of a failure',/7x,'in the calculation of the hypothetical',
     $    ' affinity of this phase',/7x,'in the presence of the',
     $    ' other phases involved in the violation.',/7x,"Couldn't",
     $    ' solve the requisite matrix equation.')
          stop
        endif
c
        if (iodb(1) .gt. 0) then
          write (noutpt,1230)
 1230     format(/3x,'The interphase reaction is:',/)
          cx = 1.0
          j2 = ilnobl(uphase(np))
          write (noutpt,1240) cx,uphase(np)(1:j2)
 1240     format(4x,f7.3,2x,a)
          j = 0
          do kpp = 1,kpt
            if (kpp .ne. kp) then
              kmm = kmprv(kpp)
              npp = ipndx1(kmm)
              j = j + 1
              cx = delvec(j)
              if (cx .lt. 0.) then
                acx = -cx
                j3 = ilnobl(uphase(npp))
                write (noutpt,1250) acx,uphase(npp)(1:j3)
 1250           format(2x,'+ ',f7.3,2x,a)
              endif
            endif
          enddo
          write (noutpt,1260)
 1260     format(10x,'==')
          j = 0
          do kpp = 1,kpt
            if (kpp .ne. kp) then
              kmm = kmprv(kpp)
              npp = ipndx1(kmm)
              j = j + 1
              cx = delvec(j)
              if (cx .gt. 0.) then
                j3 = ilnobl(uphase(npp))
                write (noutpt,1250) cx,uphase(npp)(1:j3)
              endif
            endif
          enddo
        endif
c
        affc = -xlks(ns)
        j = 0
        do kpp = 1,kpt
          if (kpp .ne. kp) then
            j = j + 1
            kmm = kmprv(kpp)
            nss = iindx1(kmm)
            affc = affc + xlks(nss)*delvec(j)
          endif
        enddo
        affp(np) = afcnst*affc
c
        if (iodb(1) .gt. 0) then
          j2 = ilnobl(uphase(np))
          write (noutpt,1270) uphase(np)(1:j2),affp(np)
 1270     format(/3x,'The hypothetical affinity for ',a,' is ',f10.5,/)
        endif
      enddo
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
c     Only phases with negative affinities are candidates to be
c     dropped. It is guaranteed that the last phase added will
c     not be one of these. The best approach to affinity scaling
c     is to find the phase which is limiting with respect to how
c     many moles can be created of the last phase added.
c
      kpadd = 0
      do kp = 1,kpt
        km = kmprv(kp)
        np = ipndx1(km)
        if (np .eq. npadd) then
          kpadd = kp
          go to 110
        endif
      enddo
c
      write (noutpt,1280)
      write (nttyo,1280)
 1280 format(/" * Note (EQ6/jgibbs) Don't know which phase was the",
     $ ' last one added'/7x,'to the ES phase assemblage. This',
     $ ' datum is needed to generate',/7x,'the optimum choice of',
     $ ' phase to delete from the assemblage',/7x,'to resolve the',
     $ ' present mineralogical phase rule violation.',/7x,'Will',
     $ ' generate a phase to delete using the available data.')
c
  110 continue
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
      if (kpadd .eq. 0) then
c
c       Use the regular affinity scaling factors.
c
        do kp = 1,kpt
          km = kmprv(kp)
          ns = iindx1(km)
          cscx(kp) = cscale(ns)
        enddo
      else
c
c       Normalize the inter-phase reaction so that the last phase
c       added has a reaction coefficient of -1.0. Note that the
c       reaction at this point is set up for the last phase involved
c       in the violation.
c
        delvec(kpt) = -1.0
        cx = delvec(kpadd)
        if (cx .lt. 0.) cx = -cx
        do kp = 1,kpt
          delvec(kp) = delvec(kp)/cx
        enddo
c
c       Calculate the maximum possible mass of each phase. Define
c       the scaling factor as equivalent to the maximum number
c       of mole equivalents of the last phase added.
c
        do kp = 1,kpt
          cscx(kp) = 1.
          km = kmprv(kp)
          ns = iindx1(km)
          np = ipndx1(km)
          mxxm = 1.e+38
          nr1 = nstsr(1,ns)
          nr2 = nstsr(2,ns)
          do n = nr1,nr2
            nb = nsts(n)
            cx = csts(n)
            mxx = mtb(nb)
            if (mxx.gt.0. .and. cx.gt.0.) then
              mxxe = mxx/cx
              if (mxxe .lt. mxxm) mxxm = mxxe
            endif
          enddo
          acx = abs(delvec(kp))
          acx2 = acx*acx
c
c         Note: dividing mxxm by acx once gives the maximum number
c         of mole equivalents of the last phase addded. Dividing
c         mxmm by acx a second time is the equivalent of multiplying
c         the hypothetical affinity by the same factor, which normalizes
c         the affinities with respect to the corresponding reaction
c         stoichiometries. Such normalization causes each of the
c         hypothetical affinities to take on the same magnitude.
c
          cscx(kp) = mxxm/acx2
        enddo
      endif
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
      if (iodb(1) .gt. 0) then
        if (kpadd .gt. 0) then
c
c         Write an example of the linking inter-phase reaction such
c         that the last phase added is on the right hand side. It is
c         presently on the left hand side, so reverse the reaction
c         coefficients.
c
          do kp = 1,kpt
            delvec(kp) = -delvec(kp)
          enddo
c
          write (noutpt,1340)
 1340     format(//3x,'The phases are linked by the reaction:',/)
c
          j = 0
          do kp = 1,kpt
            km = kmprv(kp)
            np = ipndx1(km)
            cx = delvec(kp)
            if (cx .lt. 0.) then
              acx = -cx
              j3 = ilnobl(uphase(np))
              if (j .eq. 0) then
                write (noutpt,1350) acx,uphase(np)(1:j3)
 1350           format(4x,f7.3,2x,a)
                j = j + 1
              else
                write (noutpt,1360) acx,uphase(np)(1:j3)
 1360           format(2x,'+ ',f7.3,2x,a)
              endif
            endif
          enddo
          write (noutpt,1370)
 1370     format(10x,'==')
          do kp = 1,kpt
            km = kmprv(kp)
            np = ipndx1(km)
            cx = delvec(kp)
            if (cx .gt. 0.) then
              j3 = ilnobl(uphase(np))
              write (noutpt,1360) cx,uphase(np)(1:j3)
            endif
          enddo
          write (noutpt,1380)
 1380     format(1x)
        endif
      endif
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
      if (iodb(1) .gt. 0) then
c
c       Write a table of the results.
c
        write (noutpt,1420)
 1420   format(/8x,'Phase',22x,'Affinity',2x,'Scale Factor',2x,
     $  'Scaled Affinity',/)
        do kp = 1,kpt
          km = kmprv(kp)
          np = ipndx1(km)
          afscal = affp(np)/cscx(kp)
          write (noutpt,1430) uphase(np),affp(np),cscx(kp),afscal
 1430     format(6x,a24,3x,f10.5,3x,1pe10.3,3x,e10.3)
        enddo
        if (kpadd .eq. 0) then
          write (noutpt,1440)
 1440     format(/3x,'Scaling employs the regular affinity scaling',
     $    ' factors.',/)
        else
          j2 = ilnobl(uphase(npadd))
          write (noutpt,1450) uphase(npadd)(1:j2)
 1450     format(/3x,'Scaling is based on stoichiometric normalization',
     $    ' of the reactions',/3x,'from which the hypothetical',
     $    ' affinities are calculated and on the',/3x,'maximum number',
     $    ' of mole equivalents of each phase in terms of',/3x,a,
     $    ', the last phase added to the ES phase assemblage.',/)
        endif
      endif
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
c     Determine the candidate to drop from the equilibrium phase
c     assemblage.
c
      npdel = 0
      affmin = 0.
      do kp = 1,kpt
        km = kmprv(kp)
        np = ipndx1(km)
        afscal = affp(np)/cscx(kp)
        if (afscal .lt. affmin) then
          if (affmin.ge.0. .or. npadd.ne.np) then
            affmin = afscal
            npdel = np
          endif
        endif
      enddo
c
      if (npdel .le. 0) then
        j2 = ilnobl(uphase(npdel))
        write (noutpt,1480) uphase(npdel)(1:j2)
        write (nttyo,1480) uphase(npdel)(1:j2)
 1480   format(/' * Error (EQ6/jgibbs) Programming error trap:',
     $  " Couldn't find",/7x,'a candidate phase to delete to',
     $  ' resolve the mineralogical',/7x,'phase rule violation.')
        stop
      endif
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
  999 continue
      end


text:
less_toxic(y1,h1).
less_toxic(hh1,w1).
less_toxic(ee1,b1).
less_toxic(m1,cc1).
less_toxic(bb1,z1).
less_toxic(ff1,v1).
less_toxic(ll1,b1).
less_toxic(o1,jj1).
less_toxic(j1,dd1).
less_toxic(n1,cc1).
less_toxic(w1,aa1).
less_toxic(q1,l1).
less_toxic(m1,l1).
less_toxic(bb1,aa1).
less_toxic(q1,h1).
less_toxic(b1,i1).
less_toxic(ee1,l1).
less_toxic(u1,p1).
less_toxic(ee1,jj1).
less_toxic(m1,w1).
less_toxic(q1,d1).
less_toxic(cc1,l1).
less_toxic(kk1,v1).
less_toxic(n1,t1).
less_toxic(b1,d1).
less_toxic(o1,a1).
less_toxic(b1,aa1).
less_toxic(f1,d1).
less_toxic(i1,d1).
less_toxic(ll1,t1).
less_toxic(j1,i1).
less_toxic(h1,p1).
less_toxic(n1,aa1).
less_toxic(n1,u1).
less_toxic(o1,i1).
less_toxic(y1,i1).
less_toxic(o1,e1).
less_toxic(w1,v1).
less_toxic(ff1,t1).
less_toxic(g1,cc1).
less_toxic(cc1,z1).
less_toxic(g1,l1).
less_toxic(o1,l1).
less_toxic(y1,p1).

ext: f

max_stars_repo_name: manoelfranca/cilppp

max_issues_repo_name: manoelfranca/cilppp

max_forks_repo_name: manoelfranca/cilppp


text:
c###outmuf.for
      SUBROUTINE OUTMUF
C--------------------------------
C.....THIS ROUTINE TABULATES THE "MUFS" USING A VECTORED FORMAT SCHEME
C.....WHICH PRINTS A "DASH" IF "MUFS" .LE. ONE AND THE VALUE OTHERWISE
      COMMON / FILES / LUI, LUO, LU2, LU5, LU6, LU15, LU16, LU20, LU25,
     A LU26, LU35
      COMMON / ION / IHRE, IHRO, IHRS, LUFP, METHOD, NOISE, NPSL
      COMMON/MUFS/EMUF(24),F1MUF(24),F2MUF(24),ESMUF(24),ALLMUF(24)
     1,FOT(24),XLUF(24),HPF(24),ANGMUF(24),MODMUF,SIGL(4),SIGU(4),DELMU
     2F(4),HPMUF(4),HTMUF(4),FVMUF(4),AFMUF(4),NHOPMF(4),YFOT(4),YHPF(4)
     3 ,YMUF(4)
      COMMON / OUTLAB / LABEL(11), LAYTYP(5), IEAST, INORTH, ISOUTH,
     A IWEST, LABLI, LABLJ, LABLK
      CHARACTER IEAST*1, ISOUTH*1, INORTH*1, IWEST*1, LABEL*5, LAYTYP*2,
     A LABLI*5, LABLJ*5, LABLK*5
      COMMON / OUTPRT / LINBOT(26), LINBD(14), LINTOP(15), LINTP(14),
     A GRPTYP, JOUT, LINBYP, LINES, LINMAX, LINTYP, LPAGES, NLINE
      COMMON / TIME / IT, GMT, UTIME(24), GMTR, XLMT(24), ITIM, JTX
      DIMENSION XFMT(5)
      CHARACTER NFMT(4)*6, IFMT(14)*6, IGMT*6, ILMT*6
      DATA NFMT/',5X,  ', '2H-   ', ',F7.2 ', '      '/
      DATA IFMT(1)/'(1H , '/, IFMT(2)/'2X,   '/, IFMT(3)/'2F6.1 '/,
     A IFMT(14)/'     )'/, IGMT/'   GMT'/, ILMT/'   LMT'/
      LINES = 0
C.....CALL SUBROUTINE OUTTOP TO OUTPUT THE HEADER LINES
      CALL OUTTOP
C.....PRINT THE TABULATION HEADING
      IF(METHOD - 26) 130, 120, 130
  120 WRITE(LUO,901) IGMT, ILMT, LABEL(4), LABEL(7), LABEL(3),
     1 LABEL(10), LABEL(8)
       GO TO 135
  130 WRITE(LUO,902) IGMT, ILMT, LABEL(4), LABEL(7), LABEL(3), LABEL(10)
  135 CONTINUE
C.....BEGIN HOUR LOOP
      DO 400 JT = 1,24
C.....INCREMENT HOUR IF GMT . LT. 0
      IF(UTIME(JT)) 400, 140, 140
C.....INITIALIZE LINE TO PRINT DASHES (CHANGE LATER WHERE NEEDED)
  140 KOUNT = 0
      DO 150 J = 4, 12, 2
  150 IFMT(J) = NFMT(1)
      DO 160 J = 5, 13, 2
  160 IFMT(J) = NFMT(2)
C.....SET TO TABULATE FOT
      KOUNT = KOUNT + 1
      IFMT(4) = NFMT(3)
      IFMT(5) = NFMT(4)
      XFMT(KOUNT) = FOT(JT)
C.....SET TO TABULATE HPF
      KOUNT = KOUNT + 1
      IFMT(6) = NFMT(3)
      IFMT(7) = NFMT(4)
      XFMT(KOUNT) = HPF(JT)
C.....SET TO TABULATE ESMUF
      KOUNT = KOUNT + 1
      IFMT(8) = NFMT(3)
      IFMT(9) = NFMT(4)
      XFMT(KOUNT) = ESMUF(JT)
C.....SET TO TABULATE ALLMUF
      KOUNT = KOUNT + 1
      IFMT(10) = NFMT(3)
      IFMT(11) = NFMT(4)
      XFMT(KOUNT) = ALLMUF(JT)
C.....SET TO TABULATE LUFS (IF METHOD = 26)
      IF(METHOD - 26) 250, 260, 250
  250 IFMT(12) = NFMT(4)
      IFMT(13) = NFMT(4)
      GO TO 280
  260 CONTINUE
      KOUNT = KOUNT + 1
      IFMT(12) = NFMT(3)
      IFMT(13) = NFMT(4)
      XFMT(KOUNT) = XLUF(JT)
  280 CONTINUE
      WRITE(LUO,IFMT) UTIME(JT), XLMT(JT), (XFMT(K),K=1,KOUNT)
  400 CONTINUE
C.....END OF HOUR LOOP
      RETURN
  901 FORMAT(//,' ',2X,2A6,5(2X,A5),/)
  902 FORMAT(//,' ',2X,2A6,4(2X,A5),/)
      END
C--------------------------------

ext: for

max_stars_repo_name: jawatson/voacapl

max_issues_repo_name: jawatson/voacapl

max_forks_repo_name: jawatson/voacapl


! NOTE: this example is unsupported by the implementation
program monte_carlo_pi
  use iso_fortran_env, only : STAT_FAILED_IMAGE
  implicit none
  integer, parameter :: NFAIL=2, NSAMPLES=10**6
  integer :: i, status, n_i, n[*], active_images = 1
  double precision :: x, y

  call random_init(repeatable=.true., image_distinct=.true.)

  n = 0

  do i = 1, NSAMPLES
    call random_number(x); call random_number(y)
    if (hypot(x, y) <= 1) n = n + 1
  end do

  ! simulate failure in last NFAIL images
  if (num_images() - this_image() < NFAIL) fail image

  if (this_image() == 1) then
    do i = 2, num_images()
       sync images(i, stat=status)
       n_i = n[i, stat=status]
       if (status /= STAT_FAILED_IMAGE) then
          n = n + n_i
          active_images = active_images + 1
       end if
    end do
    write(*,*) (4.0d0 * n / NSAMPLES) / active_images
  else
     sync images(1)
  end if
end program monte_carlo_pi


text:
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c      written by the UFO converter
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      DOUBLE COMPLEX FUNCTION COND(CONDITION,TRUECASE,FALSECASE)
      IMPLICIT NONE
      DOUBLE COMPLEX CONDITION,TRUECASE,FALSECASE
      IF(CONDITION.EQ.(0.0D0,0.0D0)) THEN
        COND=TRUECASE
      ELSE
        COND=FALSECASE
      ENDIF
      END

      DOUBLE COMPLEX FUNCTION CONDIF(CONDITION,TRUECASE,FALSECASE)
      IMPLICIT NONE
      LOGICAL CONDITION
      DOUBLE COMPLEX TRUECASE,FALSECASE
      IF(CONDITION) THEN
        CONDIF=TRUECASE
      ELSE
        CONDIF=FALSECASE
      ENDIF
      END

      DOUBLE COMPLEX FUNCTION RECMS(CONDITION,EXPR)
      IMPLICIT NONE
      LOGICAL CONDITION
      DOUBLE COMPLEX EXPR
      IF(CONDITION)THEN
        RECMS=EXPR
      ELSE
        RECMS=DCMPLX(DBLE(EXPR))
      ENDIF
      END

      DOUBLE COMPLEX FUNCTION REGLOG(ARG)
      IMPLICIT NONE
      DOUBLE COMPLEX TWOPII
      PARAMETER (TWOPII=2.0D0*3.1415926535897932D0*(0.0D0,1.0D0))
      DOUBLE COMPLEX ARG
      IF(ARG.EQ.(0.0D0,0.0D0)) THEN
        REGLOG=(0.0D0,0.0D0)
      ELSE
        REGLOG=LOG(ARG)
      ENDIF
      END

      DOUBLE COMPLEX FUNCTION REGLOGP(ARG)
      IMPLICIT NONE
      DOUBLE COMPLEX TWOPII
      PARAMETER (TWOPII=2.0D0*3.1415926535897932D0*(0.0D0,1.0D0))
      DOUBLE COMPLEX ARG
      IF(ARG.EQ.(0.0D0,0.0D0))THEN
        REGLOGP=(0.0D0,0.0D0)
      ELSE
        IF(DBLE(ARG).LT.0.0D0.AND.DIMAG(ARG).LT.0.0D0)THEN
          REGLOGP=LOG(ARG) + TWOPII
        ELSE
          REGLOGP=LOG(ARG)
        ENDIF
      ENDIF
      END

      DOUBLE COMPLEX FUNCTION REGLOGM(ARG)
      IMPLICIT NONE
      DOUBLE COMPLEX TWOPII
      PARAMETER (TWOPII=2.0D0*3.1415926535897932D0*(0.0D0,1.0D0))
      DOUBLE COMPLEX ARG
      IF(ARG.EQ.(0.0D0,0.0D0))THEN
        REGLOGM=(0.0D0,0.0D0)
      ELSE
        IF(DBLE(ARG).LT.0.0D0.AND.DIMAG(ARG).GT.0.0D0)THEN
          REGLOGM=LOG(ARG) - TWOPII
        ELSE
          REGLOGM=LOG(ARG)
        ENDIF
      ENDIF
      END

      DOUBLE COMPLEX FUNCTION REGSQRT(ARG_IN)
      IMPLICIT NONE
      DOUBLE COMPLEX ARG_IN
      DOUBLE COMPLEX ARG
      ARG=ARG_IN
      IF(DABS(DIMAG(ARG)).EQ.0.0D0)THEN
        ARG=DCMPLX(DBLE(ARG),0.0D0)
      ENDIF
      IF(DABS(DBLE(ARG)).EQ.0.0D0)THEN
        ARG=DCMPLX(0.0D0,DIMAG(ARG))
      ENDIF
      REGSQRT=SQRT(ARG)
      END

      DOUBLE COMPLEX FUNCTION GRREGLOG(LOGSW,EXPR1_IN,EXPR2_IN)
      IMPLICIT NONE
      DOUBLE COMPLEX TWOPII
      PARAMETER (TWOPII=2.0D0*3.1415926535897932D0*(0.0D0,1.0D0))
      DOUBLE COMPLEX EXPR1_IN,EXPR2_IN
      DOUBLE COMPLEX EXPR1,EXPR2
      DOUBLE PRECISION LOGSW
      DOUBLE PRECISION IMAGEXPR
      LOGICAL FIRSTSHEET
      EXPR1=EXPR1_IN
      EXPR2=EXPR2_IN
      IF(DABS(DIMAG(EXPR1)).EQ.0.0D0)THEN
        EXPR1=DCMPLX(DBLE(EXPR1),0.0D0)
      ENDIF
      IF(DABS(DBLE(EXPR1)).EQ.0.0D0)THEN
        EXPR1=DCMPLX(0.0D0,DIMAG(EXPR1))
      ENDIF
      IF(DABS(DIMAG(EXPR2)).EQ.0.0D0)THEN
        EXPR2=DCMPLX(DBLE(EXPR2),0.0D0)
      ENDIF
      IF(DABS(DBLE(EXPR2)).EQ.0.0D0)THEN
        EXPR2=DCMPLX(0.0D0,DIMAG(EXPR2))
      ENDIF
      IF(EXPR1.EQ.(0.0D0,0.0D0))THEN
        GRREGLOG=(0.0D0,0.0D0)
      ELSE
        IMAGEXPR=DIMAG(EXPR1)*DIMAG(EXPR2)
        FIRSTSHEET=IMAGEXPR.GE.0.0D0
        FIRSTSHEET=FIRSTSHEET.OR.DBLE(EXPR1).GE.0.0D0
        FIRSTSHEET=FIRSTSHEET.OR.DBLE(EXPR2).GE.0.0D0
        IF(FIRSTSHEET)THEN
          GRREGLOG=LOG(EXPR1)
        ELSE
          IF(DIMAG(EXPR1).GT.0.0D0)THEN
            GRREGLOG=LOG(EXPR1) - LOGSW*TWOPII
          ELSE
            GRREGLOG=LOG(EXPR1) + LOGSW*TWOPII
          ENDIF
        ENDIF
      ENDIF
      END

      DOUBLE COMPLEX FUNCTION ARG(COMNUM)
      IMPLICIT NONE
      DOUBLE COMPLEX COMNUM
      DOUBLE COMPLEX IIM
      IIM = (0.0D0,1.0D0)
      IF(COMNUM.EQ.(0.0D0,0.0D0)) THEN
        ARG=(0.0D0,0.0D0)
      ELSE
        ARG=LOG(COMNUM/ABS(COMNUM))/IIM
      ENDIF
      END

      COMPLEX*32 FUNCTION MP_COND(CONDITION,TRUECASE,FALSECASE)
      IMPLICIT NONE
      COMPLEX*32 CONDITION,TRUECASE,FALSECASE
      IF(CONDITION.EQ.(0.0E0_16,0.0E0_16)) THEN
        MP_COND=TRUECASE
      ELSE
        MP_COND=FALSECASE
      ENDIF
      END

      COMPLEX*32 FUNCTION MP_CONDIF(CONDITION,TRUECASE,FALSECASE)
      IMPLICIT NONE
      LOGICAL CONDITION
      COMPLEX*32 TRUECASE,FALSECASE
      IF(CONDITION) THEN
        MP_CONDIF=TRUECASE
      ELSE
        MP_CONDIF=FALSECASE
      ENDIF
      END

      COMPLEX*32 FUNCTION MP_RECMS(CONDITION,EXPR)
      IMPLICIT NONE
      LOGICAL CONDITION
      COMPLEX*32 EXPR
      IF(CONDITION)THEN
        MP_RECMS=EXPR
      ELSE
        MP_RECMS=CMPLX(REAL(EXPR),KIND=16)
      ENDIF
      END

      COMPLEX*32 FUNCTION MP_REGLOG(ARG)
      IMPLICIT NONE
      COMPLEX*32 TWOPII
      PARAMETER (TWOPII=2.0E0_16
     $ *3.14169258478796109557151794433593750E0_16*(0.0E0_16,1.0E0_16))
      COMPLEX*32 ARG
      IF(ARG.EQ.(0.0E0_16,0.0E0_16)) THEN
        MP_REGLOG=(0.0E0_16,0.0E0_16)
      ELSE
        MP_REGLOG=LOG(ARG)
      ENDIF
      END

      COMPLEX*32 FUNCTION MP_REGLOGP(ARG)
      IMPLICIT NONE
      COMPLEX*32 TWOPII
      PARAMETER (TWOPII=2.0E0_16
     $ *3.14169258478796109557151794433593750E0_16*(0.0E0_16,1.0E0_16))
      COMPLEX*32 ARG
      IF(ARG.EQ.(0.0E0_16,0.0E0_16))THEN
        MP_REGLOGP=(0.0E0_16,0.0E0_16)
      ELSE
        IF(REAL(ARG,KIND=16).LT.0.0E0_16.AND.IMAGPART(ARG).LT.0.0E0_16)
     $   THEN
          MP_REGLOGP=LOG(ARG) + TWOPII
        ELSE
          MP_REGLOGP=LOG(ARG)
        ENDIF
      ENDIF
      END

      COMPLEX*32 FUNCTION MP_REGLOGM(ARG)
      IMPLICIT NONE
      COMPLEX*32 TWOPII
      PARAMETER (TWOPII=2.0E0_16
     $ *3.14169258478796109557151794433593750E0_16*(0.0E0_16,1.0E0_16))
      COMPLEX*32 ARG
      IF(ARG.EQ.(0.0E0_16,0.0E0_16))THEN
        MP_REGLOGM=(0.0E0_16,0.0E0_16)
      ELSE
        IF(REAL(ARG,KIND=16).LT.0.0E0_16.AND.IMAGPART(ARG).GT.0.0E0_16)
     $   THEN
          MP_REGLOGM=LOG(ARG) - TWOPII
        ELSE
          MP_REGLOGM=LOG(ARG)
        ENDIF
      ENDIF
      END

      COMPLEX*32 FUNCTION MP_REGSQRT(ARG_IN)
      IMPLICIT NONE
      COMPLEX*32 ARG_IN
      COMPLEX*32 ARG
      ARG=ARG_IN
      IF(ABS(IMAGPART(ARG)).EQ.0.0E0_16)THEN
        ARG=CMPLX(REAL(ARG,KIND=16),0.0E0_16)
      ENDIF
      IF(ABS(REAL(ARG,KIND=16)).EQ.0.0E0_16)THEN
        ARG=CMPLX(0.0E0_16,IMAGPART(ARG))
      ENDIF
      MP_REGSQRT=SQRT(ARG)
      END

      COMPLEX*32 FUNCTION MP_GRREGLOG(LOGSW,EXPR1_IN,EXPR2_IN)
      IMPLICIT NONE
      COMPLEX*32 TWOPII
      PARAMETER (TWOPII=2.0E0_16
     $ *3.14169258478796109557151794433593750E0_16*(0.0E0_16,1.0E0_16))
      COMPLEX*32 EXPR1_IN,EXPR2_IN
      COMPLEX*32 EXPR1,EXPR2
      REAL*16 LOGSW
      REAL*16 IMAGEXPR
      LOGICAL FIRSTSHEET
      EXPR1=EXPR1_IN
      EXPR2=EXPR2_IN
      IF(ABS(IMAGPART(EXPR1)).EQ.0.0E0_16)THEN
        EXPR1=CMPLX(REAL(EXPR1,KIND=16),0.0E0_16)
      ENDIF
      IF(ABS(REAL(EXPR1,KIND=16)).EQ.0.0E0_16)THEN
        EXPR1=CMPLX(0.0E0_16,IMAGPART(EXPR1))
      ENDIF
      IF(ABS(IMAGPART(EXPR2)).EQ.0.0E0_16)THEN
        EXPR2=CMPLX(REAL(EXPR2,KIND=16),0.0E0_16)
      ENDIF
      IF(ABS(REAL(EXPR2,KIND=16)).EQ.0.0E0_16)THEN
        EXPR2=CMPLX(0.0E0_16,IMAGPART(EXPR2))
      ENDIF
      IF(EXPR1.EQ.(0.0E0_16,0.0E0_16))THEN
        MP_GRREGLOG=(0.0E0_16,0.0E0_16)
      ELSE
        IMAGEXPR=IMAGPART(EXPR1)*IMAGPART(EXPR2)
        FIRSTSHEET=IMAGEXPR.GE.0.0E0_16
        FIRSTSHEET=FIRSTSHEET.OR.REAL(EXPR1,KIND=16).GE.0.0E0_16
        FIRSTSHEET=FIRSTSHEET.OR.REAL(EXPR2,KIND=16).GE.0.0E0_16
        IF(FIRSTSHEET)THEN
          MP_GRREGLOG=LOG(EXPR1)
        ELSE
          IF(IMAGPART(EXPR1).GT.0.0E0_16)THEN
            MP_GRREGLOG=LOG(EXPR1) - LOGSW*TWOPII
          ELSE
            MP_GRREGLOG=LOG(EXPR1) + LOGSW*TWOPII
          ENDIF
        ENDIF
      ENDIF
      END

      COMPLEX*32 FUNCTION MP_ARG(COMNUM)
      IMPLICIT NONE
      COMPLEX*32 COMNUM
      COMPLEX*32 IMM
      IMM = (0.0E0_16,1.0E0_16)
      IF(COMNUM.EQ.(0.0E0_16,0.0E0_16)) THEN
        MP_ARG=(0.0E0_16,0.0E0_16)
      ELSE
        MP_ARG=LOG(COMNUM/ABS(COMNUM))/IMM
      ENDIF
      END


text:
!====================================================
SUBROUTINE splper(n,h,y,bi,ci,di)

! Makes a cubic spline of periodic function y(x)
!
! Input:  n                   number of values in y
!         h                   step size in x (aequidistant)
!         y(n)                y-values
! Output: bi(n),ci(n),di(n)   Spline parameters

  USE neo_precision

  IMPLICIT NONE

  INTEGER,                     INTENT(in)  :: n
  REAL(kind=dp),               INTENT(in)  :: h
  REAL(kind=dp), DIMENSION(n), INTENT(in)  :: y
  REAL(kind=dp), DIMENSION(n), INTENT(out) :: bi, ci, di

  REAL(kind=dp)                            :: psi, ss
  REAL(kind=dp), DIMENSION(:), ALLOCATABLE :: bmx, yl
  REAL(kind=dp), DIMENSION(:), ALLOCATABLE :: amx1, amx2, amx3
  INTEGER                                  :: nmx, n1, n2, i, i1

  ALLOCATE ( bmx(n), yl(n), amx1(n), amx2(n), amx3(n) )

  bmx(1) = 1.e30_dp

  nmx=n-1
  n1=nmx-1
  n2=nmx-2
  psi=3.0_dp/h/h

  CALL spfper(n,amx1,amx2,amx3)

  bmx(nmx) = (y(nmx+1)-2*y(nmx)+y(nmx-1))*psi
  bmx(1)   =(y(2)-y(1)-y(nmx+1)+y(nmx))*psi
  DO i = 3,nmx
     bmx(i-1) = (y(i)-2*y(i-1)+y(i-2))*psi
  END DO
  yl(1) = bmx(1)/amx1(1)
  DO i = 2,n1
     i1 = i-1
     yl(i) = (bmx(i)-yl(i1)*amx2(i1))/amx1(i)
  END DO
  ss = 0
  DO i = 1,n1
     ss = ss+yl(i)*amx3(i)
  END DO
  yl(nmx) = (bmx(nmx)-ss)/amx1(nmx)
  bmx(nmx) = yl(nmx)/amx1(nmx)
  bmx(n1) = (yl(n1)-amx2(n1)*bmx(nmx))/amx1(n1)
  DO i = n2,1,-1
     bmx(i) = (yl(i)-amx3(i)*bmx(nmx)-amx2(i)*bmx(i+1))/amx1(i)
  END DO
  DO i = 1,nmx
     ci(i) = bmx(i)
  END DO

  DO i = 1,n1
     bi(i) = (y(i+1)-y(i))/h-h*(ci(i+1)+2*ci(i))/3
     di(i) = (ci(i+1)-ci(i))/h/3
  END DO
  bi(nmx) = (y(n)-y(n-1))/h-h*(ci(1)+2*ci(nmx))/3
  di(nmx) = (ci(1)-ci(nmx))/h/3
!
! Fix of problems at upper periodicity boundary
!
  bi(n) = bi(1)
  ci(n) = ci(1)
  di(n) = di(1)

  DEALLOCATE ( bmx, yl, amx1, amx2, amx3 )

  RETURN
END SUBROUTINE splper

max_stars_repo_name: joseluisvelasco/STELLOPT

max_issues_repo_name: joseluisvelasco/STELLOPT

max_forks_repo_name: joseluisvelasco/STELLOPT


text:
      SUBROUTINE READLC(LUNIT,CHR,STR)

C$$$  SUBPROGRAM DOCUMENTATION BLOCK
C
C SUBPROGRAM:    READLC
C   PRGMMR: WOOLLEN          ORG: NP20       DATE: 2003-11-04
C
C ABSTRACT: THIS SUBROUTINE RETURNS A CHARACTER DATA ELEMENT ASSOCIATED
C   WITH A PARTICULAR SUBSET MNEMONIC FROM THE INTERNAL MESSAGE BUFFER
C   (ARRAY MBAY IN COMMON BLOCK /BITBUF/).  IT IS DESIGNED TO BE USED
C   TO RETURN CHARACTER ELEMENTS GREATER THAN THE USUAL LENGTH OF EIGHT
C   BYTES.
C
C PROGRAM HISTORY LOG:
C 2003-11-04  J. WOOLLEN -- ORIGINAL AUTHOR
C 2003-11-04  D. KEYSER  -- UNIFIED/PORTABLE FOR WRF; ADDED
C                           DOCUMENTATION; OUTPUTS MORE COMPLETE
C                           DIAGNOSTIC INFO WHEN ROUTINE TERMINATES
C                           ABNORMALLY OR UNUSUAL THINGS HAPPEN
C 2004-08-09  J. ATOR    -- MAXIMUM MESSAGE LENGTH INCREASED FROM
C                           20,000 TO 50,000 BYTES
C 2007-01-19  J. ATOR    -- REPLACED CALL TO PARSEQ WITH CALL TO PARSTR
C 2009-03-23  J. ATOR    -- ADDED CAPABILITY FOR COMPRESSED MESSAGES;
C                           ADDED CHECK FOR OVERFLOW OF CHR; ADDED '#'
C                           OPTION FOR MORE THAN ONE OCCURRENCE OF STR
C 2009-04-21  J. ATOR    -- USE ERRWRT
C 2012-12-07  J. ATOR    -- ALLOW STR MNEMONIC LENGTH OF UP TO 14 CHARS
C                           WHEN USED WITH '#' OCCURRENCE CODE
C
C USAGE:    CALL READLC (LUNIT, CHR, STR)
C   INPUT ARGUMENT LIST:
C     LUNIT    - INTEGER: FORTRAN LOGICAL UNIT NUMBER FOR BUFR FILE
C     STR      - CHARACTER*(*): STRING (I.E., MNEMONIC)
C
C   OUTPUT ARGUMENT LIST:
C     CHR      - CHARACTER*(*): UNPACKED CHARACTER STRING (I.E.,
C                CHARACTER DATA ELEMENT GREATER THAN EIGHT BYTES)
C
C REMARKS:
C    THIS ROUTINE CALLS:        BORT     ERRWRT   PARSTR   PARUTG
C                               STATUS   UPC
C    THIS ROUTINE IS CALLED BY: UFBDMP   UFDUMP   WRTREE
C                               Also called by application programs.
C
C ATTRIBUTES:
C   LANGUAGE: FORTRAN 77
C   MACHINE:  PORTABLE TO ALL PLATFORMS
C
C$$$

      INCLUDE 'bufrlib.prm'

      COMMON /BITBUF/ MAXBYT,IBIT,IBAY(MXMSGLD4),MBYT(NFILES),
     .                MBAY(MXMSGLD4,NFILES)
      COMMON /TABLES/ MAXTAB,NTAB,TAG(MAXJL),TYP(MAXJL),KNT(MAXJL),
     .                JUMP(MAXJL),LINK(MAXJL),JMPB(MAXJL),
     .                IBT(MAXJL),IRF(MAXJL),ISC(MAXJL),
     .                ITP(MAXJL),VALI(MAXJL),KNTI(MAXJL),
     .                ISEQ(MAXJL,2),JSEQ(MAXJL)
      COMMON /USRINT/ NVAL(NFILES),INV(MAXSS,NFILES),VAL(MAXSS,NFILES)
      COMMON /RLCCMN/ NRST,IRNCH(MXRST),IRBIT(MXRST),CRTAG(MXRST)
      COMMON /USRBIT/ NBIT(MAXSS),MBIT(MAXSS)
      COMMON /UNPTYP/ MSGUNP(NFILES)
      COMMON /QUIET / IPRT

      CHARACTER*(*) CHR,STR
      CHARACTER*128 BORT_STR,ERRSTR
      CHARACTER*10  TAG,CTAG,CRTAG
      CHARACTER*14  TGS(10)
      CHARACTER*3   TYP
      REAL*8        VAL

      DATA MAXTG /10/

C-----------------------------------------------------------------------
C-----------------------------------------------------------------------

      CHR = ' '

C  CHECK THE FILE STATUS
C  ---------------------

      CALL STATUS(LUNIT,LUN,IL,IM)
      IF(IL.EQ.0) GOTO 900
      IF(IL.GT.0) GOTO 901
      IF(IM.EQ.0) GOTO 902

C  CHECK FOR TAGS (MNEMONICS) IN INPUT STRING (THERE CAN ONLY BE ONE)
C  ------------------------------------------------------------------

      CALL PARSTR(STR,TGS,MAXTG,NTG,' ',.TRUE.)
      IF(NTG.GT.1) GOTO 903

C     Check if a specific occurrence of the input string was requested;
C     if not, then the default is to return the first occurrence.

      CALL PARUTG(LUN,0,TGS(1),NNOD,KON,ROID)
      IF(KON.EQ.6) THEN
         IOID=NINT(ROID)
         IF(IOID.LE.0) IOID = 1
         CTAG = ' '
         II = 1
         DO WHILE((II.LE.10).AND.(TGS(1)(II:II).NE.'#'))
            CTAG(II:II)=TGS(1)(II:II)
            II = II + 1
         ENDDO
      ELSE
         IOID = 1
         CTAG = TGS(1)(1:10)
      ENDIF

C  LOCATE AND DECODE THE LONG CHARACTER STRING
C  -------------------------------------------

      IF(MSGUNP(LUN).EQ.0.OR.MSGUNP(LUN).EQ.1) THEN

C        The message is uncompressed

         ITAGCT = 0
         DO N=1,NVAL(LUN)
           NOD = INV(N,LUN)
           IF(CTAG.EQ.TAG(NOD)) THEN
             ITAGCT = ITAGCT + 1
             IF(ITAGCT.EQ.IOID) THEN
               IF(ITP(NOD).NE.3) GOTO 904
               NCHR = NBIT(N)/8
               IF(NCHR.GT.LEN(CHR)) GOTO 905
               KBIT = MBIT(N)
               CALL UPC(CHR,NCHR,MBAY(1,LUN),KBIT)
               GOTO 100
             ENDIF
           ENDIF
         ENDDO
      ELSEIF(MSGUNP(LUN).EQ.2) THEN

C        The message is compressed

         IF(NRST.GT.0) THEN
           ITAGCT = 0
           DO II=1,NRST
             IF(CTAG.EQ.CRTAG(II)) THEN
               ITAGCT = ITAGCT + 1
               IF(ITAGCT.EQ.IOID) THEN
                 NCHR = IRNCH(II)
                 IF(NCHR.GT.LEN(CHR)) GOTO 905
                 KBIT = IRBIT(II)
                 CALL UPC(CHR,NCHR,MBAY(1,LUN),KBIT)
                 GOTO 100
               ENDIF
             ENDIF
           ENDDO
         ENDIF
      ELSE
         GOTO 906
      ENDIF

C     If we made it here, then we couldn't find the requested string.

      IF(IPRT.GE.0) THEN
      CALL ERRWRT('++++++++++++++BUFR ARCHIVE LIBRARY+++++++++++++++++')
      ERRSTR = 'BUFRLIB: READLC - MNEMONIC ' // TGS(1) //
     .   ' NOT LOCATED IN REPORT SUBSET - RETURN WITH BLANK' //
     .   ' STRING FOR CHARACTER DATA ELEMENT'
      CALL ERRWRT(ERRSTR)
      CALL ERRWRT('++++++++++++++BUFR ARCHIVE LIBRARY+++++++++++++++++')
      CALL ERRWRT(' ')
      ENDIF

C  EXITS
C  -----

100   RETURN
900   CALL BORT('BUFRLIB: READLC - INPUT BUFR FILE IS CLOSED, IT MUST'//
     . ' BE OPEN FOR INPUT')
901   CALL BORT('BUFRLIB: READLC - INPUT BUFR FILE IS OPEN FOR '//
     . 'OUTPUT, IT MUST BE OPEN FOR INPUT')
902   CALL BORT('BUFRLIB: READLC - A MESSAGE MUST BE OPEN IN INPUT '//
     . 'BUFR FILE, NONE ARE')
903   WRITE(BORT_STR,'("BUFRLIB: READLC - THERE CANNOT BE MORE THAN '//
     . 'ONE MNEMONIC IN THE INPUT STRING (",A,") (HERE THERE ARE ",'//
     . 'I3,")")') STR,NTG
      CALL BORT(BORT_STR)
904   WRITE(BORT_STR,'("BUFRLIB: READLC - MNEMONIC ",A," DOES NOT '//
     . 'REPRESENT A CHARACTER ELEMENT (ITP=",I2,")")') TGS(1),ITP(NOD)
      CALL BORT(BORT_STR)
905   WRITE(BORT_STR,'("BUFRLIB: READLC - MNEMONIC ",A," IS A '//
     . 'CHARACTER STRING OF LENGTH",I4," BUT SPACE WAS PROVIDED '//
     . 'FOR ONLY",I4, " CHARACTERS")') TGS(1),NCHR,LEN(CHR)
      CALL BORT(BORT_STR)
906   WRITE(BORT_STR,'("BUFRLIB: READLC - MESSAGE UNPACK TYPE",I3,'//
     . '" IS NOT RECOGNIZED")') MSGUNP
      CALL BORT(BORT_STR)
      END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
      SUBROUTINE FIXPQF(N,Y,IFLAG,ARCRE,ARCAE,ANSRE,ANSAE,TRACE,A,
     $     NFE,ARCLEN,YP,YOLD,YPOLD,QT,R,F0,F1,Z0,DZ,W,T,YSAV,
     $     SSPAR,PAR,IPAR)
C
C SUBROUTINE  FIXPQF  FINDS A FIXED POINT OR ZERO OF THE
C N-DIMENSIONAL VECTOR FUNCTION  F(X), OR TRACKS A ZERO CURVE OF A
C GENERAL HOMOTOPY MAP  RHO(A,LAMBDA,X).  FOR THE FIXED POINT PROBLEM
C F(X) IS ASSUMED TO BE A C2 MAP OF SOME BALL INTO ITSELF.  THE
C EQUATION  X=F(X)  IS SOLVED BY FOLLOWING THE ZERO CURVE OF THE
C HOMOTOPY MAP
C
C  LAMBDA*(X - F(X)) + (1 - LAMBDA)*(X - A),
C
C STARTING FROM  LAMBDA = 0, X = A.   THE CURVE IS PARAMETERIZED
C BY ARC LENGTH  S, AND IS FOLLOWED BY SOLVING THE ORDINARY
C DIFFERENTIAL EQUATION  D(HOMOTOPY MAP)/DS = 0  FOR
C Y(S) = (LAMBDA(S), X(S)).  THIS IS DONE BY USING A HERMITE CUBIC
C PREDICTOR AND A CORRECTOR WHICH RETURNS TO THE ZERO CURVE IN A
C HYPERPLANE PERPENDICULAR TO THE TANGENT TO THE ZERO CURVE AT THE
C MOST RECENT POINT.
C
C FOR THE ZERO FINDING PROBLEM  F(X)  IS ASSUMED TO BE A C2 MAP
C SUCH THAT FOR SOME  R > 0,  X*F(X) >= 0  WHENEVER  NORM(X) = R.
C THE EQUATION  F(X) = 0  IS SOLVED BY FOLLOWING THE ZERO CURVE OF
C THE HOMOTOPY MAP
C
C  LAMBDA*F(X) + (1 - LAMBDA)*(X - A)
C
C EMANATING FROM  LAMBDA = 0, X = A.
C
C A  MUST BE AN INTERIOR POINT OF THE ABOVE MENTIONED BALLS.
C
C FOR THE CURVE TRACKING PROBLEM RHO(A,LAMBDA,X) IS ASSUMED TO
C BE A C2 MAP FROM  E**M X [0,1) X E**N  INTO  E**N, WHICH FOR
C ALMOST ALL PARAMETER VECTORS  A  IN SOME NONEMPTY OPEN SUBSET
C OF E**M SATISFIES
C
C  RANK [D RHO(A,LAMBDA,X)/D LAMBDA, D RHO(A,LAMBDA,X)/DX] = N
C
C FOR ALL POINTS  (LAMBDA,X)  SUCH THAT  RHO(A,LAMBDA,X) = 0.  IT IS
C FURTHER ASSUMED THAT
C
C         RANK [ D RHO(A,0,X0)/DX ] = N.
C
C WITH  A  FIXED, THE ZERO CURVE OF  RHO(A,LAMBDA,X)  EMANATING FROM
C LAMBDA = 0, X = X0  IS TRACKED UNTIL  LAMBDA = 1  BY SOLVING THE
C ORDINARY DIFFERENTIAL EQUATION    D RHO(A,LAMBDA(S),X(S))/DS = 0
C FOR  Y(S) = (LAMBDA(S), X(S)), WHERE  S  IS ARC LENGTH ALONG THE
C ZERO CURVE.  ALSO THE HOMOTOPY MAP  RHO(A,LAMBDA,X)  IS ASSUMED TO
C BE CONSTRUCTED SUCH THAT
C
C         D LAMBDA(0)/DS > 0.
C
C FOR THE FIXED POINT AND ZERO FINDING PROBLEMS, THE USER MUST SUPPLY
C A SUBROUTINE  F(X,V)  WHICH EVALUATES  F(X)  AT  X  AND RETURNS THE
C VECTOR F(X) IN  V, AND A SUBROUTINE  FJAC(X,V,K)  WHICH RETURNS IN  V
C THE KTH COLUMN OF THE JACOBIAN MATRIX OF F(X) EVALUATED AT X.  FOR
C THE CURVE TRACKING PROBLEM, THE USER MUST SUPPLY A SUBROUTINE
C RHO(A,LAMBDA,X,V,PAR,IPAR)  WHICH EVALUATES THE HOMOTOPY MAP  RHO AT
C (A,LAMBDA,X)  AND RETURNS THE VECTOR  RHO(A,LAMBDA,X)  IN  V, AND
C A SUBROUTINE  RHOJAC(A,LAMBDA,X,V,K,PAR,IPAR)  WHICH RETURNS IN  V
C THE KTH COLUMN OF THE  N X (N+1)  JACOBIAN MATRIX
C [D RHO/D LAMBDA, D RHO/DX]  EVALUATED AT  (A,LAMBDA,X).  FIXPQF
C DIRECTLY OR INDIRECTLY USES THE SUBROUTINES  D1MACH, F (OR RHO),
C FJAC (OR RHOJAC), QRFAQF, QRSLQF, ROOT, ROOTQF, STEPQF, TANGQF,
C UPQRQF  AND THE BLAS ROUTINES  DAXPY, DCOPY, DDOT, DNRM2, AND  DSCAL.
C ONLY  D1MACH  CONTAINS MACHINE DEPENDENT CONSTANTS.  NO OTHER
C MODIFICATIONS BY THE USER ARE REQUIRED.
C
C
C ON INPUT:
C
C N  IS THE DIMENSION OF X, F(X), AND RHO(A,LAMBDA,X).
C
C Y(1:N+1)  CONTAINS THE STARTING POINT FOR TRACKING THE HOMOTOPY MAP.
C    (Y(2),...,Y(N+1)) = A  FOR THE FIXED POINT AND ZERO FINDING
C    PROBLEMS.  (Y(2),...,Y(N+1)) = X0  FOR THE CURVE TRACKING PROBLEM.
C    Y(1)  NEED NOT BE DEFINED BY THE USER.
C
C IFLAG CAN BE -2, -1, 0, 2, OR 3.  IFLAG SHOULD BE 0 ON THE FIRST
C    CALL TO  FIXPQF  FOR THE PROBLEM  X=F(X), -1 FOR THE PROBLEM
C    F(X)=0, AND -2 FOR THE PROBLEM  RHO(A,LAMBDA,X)=0.   IN CERTAIN
C    SITUATIONS  IFLAG  IS SET TO 2 OR 3 BY  FIXPQF, AND  FIXPQF  CAN
C    BE CALLED AGAIN WITHOUT CHANGING  IFLAG.
C
C ARCRE, ARCAE  ARE THE RELATIVE AND ABSOLUTE ERRORS, RESPECTIVELY,
C    ALLOWED THE QUASI-NEWTON ITERATION ALONG THE ZERO CURVE.  IF
C    ARC?E .LE. 0.0  ON INPUT, IT IS RESET TO  .5*SQRT(ANS?E).
C    NORMALLY  ARC?E  SHOULD BE CONSIDERABLY LARGER THAN  ANS?E.
C
C ANSRE, ANSAE  ARE THE RELATIVE AND ABSOLUTE ERROR VALUES USED FOR
C    THE ANSWER AT  LAMBDA = 1.  THE ACCEPTED ANSWER  Y = (LAMBDA, X)
C    SATISFIES
C
C      |Y(1) - 1| .LE. ANSRE + ANSAE      .AND.
C
C      ||DZ|| .LE. ANSRE*||Y|| + ANSAE      WHERE
C
C      DZ IS THE QUASI-NEWTON STEP TO Y.
C
C TRACE  IS AN INTEGER SPECIFYING THE LOGICAL I/O UNIT FOR
C    INTERMEDIATE OUTPUT.  IF  TRACE .GT. 0  THE POINTS COMPUTED ON
C    THE ZERO CURVE ARE WRITTEN TO I/O UNIT  TRACE .
C
C A(1:*)  CONTAINS THE PARAMETER VECTOR  A.  FOR THE FIXED POINT
C    AND ZERO FINDING PROBLEMS,  A  NEED NOT BE INITIALIZED BY THE
C    USER, AND IS ASSUMED TO HAVE LENGTH  N.  FOR THE CURVE
C    TRACKING PROBLEM,  A  MUST BE INITIALIZED BY THE USER.
C
C YP(1:N+1)  IS A WORK ARRAY CONTAINING THE TANGENT VECTOR TO THE
C    ZERO CURVE AT THE CURRENT POINT  Y.
C
C YOLD(1:N+1) IS A WORK ARRAY CONTAINING THE PREVIOUS POINT FOUND
C    ON THE ZERO CURVE.
C
C YPOLD(1:N+1) IS A WORK ARRAY CONTAINING THE TANGENT VECTOR TO
C    THE ZERO CURVE AT  YOLD.
C
C QT(1:N+1,1:N+1), R((N+1)*(N+2)/2), F0(1:N+1), F1(1:N+1), Z0(1:N+1),
C    DZ(1:N+1), W(1:N+1), T(1:N+1), YSAV(1:N+1)  ARE ALL WORK ARRAYS
C    USED BY  STEPQF, TANGQF AND ROOTQF TO CALCULATE THE TANGENT
C    VECTORS AND QUASI-NEWTON STEPS.
C
C SSPAR(1:4) =  (HMIN, HMAX, BMIN, BMAX)  IS A VECTOR OF PARAMETERS
C    USED FOR THE OPTIMAL STEP SIZE ESTIMATION.  A DEFAULT VALUE
C    CAN BE SPECIFIED FOR ANY OF THESE FOUR PARAMETERS BY SETTING IT
C    .LE. 0.0  ON INPUT.  SEE THE COMMENTS IN  STEPQF  FOR MORE
C    INFORMATION ABOUT THESE PARAMETERS.
C
C PAR(1:*)  AND  IPAR(1:*)  ARE ARRAYS FOR (OPTIONAL) USER PARAMETERS,
C    WHICH ARE SIMPLY PASSED THROUGH TO THE USER WRITTEN SUBROUTINES
C    RHO, RHOJAC.
C
C
C ON OUTPUT:
C
C N , TRACE , A  ARE UNCHANGED.
C
C Y(1) = LAMBDA, (Y(2),...,Y(N+1)) = X, AND  Y  IS AN APPROXIMATE
C    ZERO OF THE HOMOTOPY MAP.  NORMALLY  LAMBDA = 1  AND  X  IS A
C    FIXED POINT OR ZERO OF  F(X).   IN ABNORMAL SITUATIONS,  LAMBDA
C    MAY ONLY BE NEAR 1 AND  X  NEAR A FIXED POINT OR ZERO.
C
C IFLAG =
C
C   1   NORMAL RETURN
C
C   2   SPECIFIED ERROR TOLERANCE CANNOT BE MET.  SOME OR ALL OF
C       ARCRE, ARCAE, ANSRE, ANSAE  HAVE BEEN INCREASED TO
C       SUITABLE VALUES.  TO CONTINUE, JUST CALL  FIXPQF  AGAIN
C       WITHOUT CHANGING ANY PARAMETERS.
C
C   3   STEPQF  HAS BEEN CALLED 1000 TIMES.  TO CONTINUE, CALL
C       FIXPQF  AGAIN WITHOUT CHANGING ANY PARAMETERS.
C
C   4   JACOBIAN MATRIX DOES NOT HAVE FULL RANK.  THE ALGORITHM
C       HAS FAILED (THE ZERO CURVE OF THE HOMOTOPY MAP CANNOT BE
C       FOLLOWED ANY FURTHER).
C
C   5   THE TRACKING ALGORITHM HAS LOST THE ZERO CURVE OF THE
C       HOMOTOPY MAP AND IS NOT MAKING PROGRESS.  THE ERROR
C       TOLERANCES  ARC?E  AND  ANS?E  WERE TOO LENIENT.  THE PROBLEM
C       SHOULD BE RESTRARTED BY CALLING  FIXPQF  WITH SMALLER ERROR
C       TOLERANCES AND  IFLAG = 0 (-1, -2).
C
C   6   THE QUASI-NEWTON ITERATION IN  STEPQF  OR  ROOTQF  FAILED TO
C       CONVERGE.  THE ERROR TOLERANCES  ANS?E  MAY BE TOO STRINGENT.
C
C   7   ILLEGAL INPUT PARAMETERS, A FATAL ERROR.
C
C ARCRE, ARCAE, ANSRE, ANSAE  ARE UNCHANGED AFTER A NORMAL RETURN
C    (IFLAG = 1).  THEY ARE INCREASED TO APPROPRIATE VALUES ON THE
C    RETURN  IFLAG = 2.
C
C NFE  IS THE NUMBER OF JACOBIAN EVALUATIONS.
C
C ARCLEN  IS THE APPROXIMATE LENGTH OF THE ZERO CURVE.
C
C ***** DECLARATIONS *****
C
C     FUNCTION DECLARATIONS
C
        DOUBLE PRECISION D1MACH, DNRM2
C
C     LOCAL VARIABLES
C
        DOUBLE PRECISION ABSERR, H, HOLD, RELERR, S, WK
        INTEGER IFLAGC, ITER, JW, LIMITD, LIMIT, NP1
        LOGICAL CRASH, START
C
C     SCALAR ARGUMENTS
C
        DOUBLE PRECISION ARCRE, ARCAE, ANSRE, ANSAE, ARCLEN
        INTEGER N,IFLAG,TRACE,NFE
C
C     ARRAY DECLARATIONS
C
        DOUBLE PRECISION A(N), Y(N+1), YP(N+1), YOLD(N+1), YPOLD(N+1),
     $    QT(N+1,N+1), R((N+1)*(N+2)/2), F0(N+1), F1(N+1), Z0(N+1),
     $    DZ(N+1), W(N+1), T(N+1), YSAV(N+1), SSPAR(4), PAR(1)
        INTEGER IPAR(1)
C
        SAVE
C
C ***** END OF DECLARATIONS *****
C
C LIMITD IS AN UPPER BOUND ON THE NUMBER OF STEPS.  IT MAY BE
C CHANGED BY CHANGING THE FOLLOWING PARAMETER STATEMENT:
        PARAMETER (LIMITD =1000)
C
C ***** FIRST EXECUTABLE STATEMENT *****
C
C CHECK IFLAG
C
        IF (N .LE. 0 .OR. ANSRE .LE. 0.0 .OR. ANSAE .LT. 0.0)
     $    IFLAG = 7
        IF (IFLAG .GE. -2 .AND. IFLAG .LE. 0) GO TO 10
        IF (IFLAG .EQ. 2) GO TO 50
        IF (IFLAG .EQ. 3) GO TO 40
C
C ONLY VALID INPUT FOR IFLAG IS -2, -1, 0, 2, 3.
C
        IFLAG = 7
        RETURN
C
C ***** INITIALIZATION BLOCK  *****
C
  10    ARCLEN = 0.0
        IF (ARCRE .LE. 0.0) ARCRE = .5*SQRT(ANSRE)
        IF (ARCAE .LE. 0.0) ARCAE = .5*SQRT(ANSAE)
        NFE=0
        IFLAGC = IFLAG
        NP1=N+1
C
C SET INITIAL CONDITIONS FOR FIRST CALL TO STEPQF.
C
        START=.TRUE.
        CRASH=.FALSE.
        RELERR = ARCRE
        ABSERR = ARCAE
        HOLD=1.0
        H=0.1
        S=0.0
        YPOLD(1) = 1.0
        Y(1) = 0.0
        DO 20 JW=2,NP1
          YPOLD(JW)=0.0
  20    CONTINUE
C
C SET OPTIMAL STEP SIZE ESTIMATION PARAMETERS.
C
C     MINIMUM STEP SIZE HMIN
        IF (SSPAR(1) .LE. 0.0) SSPAR(1)= (SQRT(N+1.0)+4.0)*D1MACH(4)
C     MAXIMUM STEP SIZE HMAX
        IF (SSPAR(2) .LE. 0.0) SSPAR(2)= 1.0
C     MINIMUM STEP REDUCTION FACTOR BMIN
        IF (SSPAR(3) .LE. 0.0) SSPAR(3)= 0.1
C     MAXIMUM STEP EXPANSION FACTOR BMAX
        IF (SSPAR(4) .LE. 0.0) SSPAR(4)= 7.0
C
C LOAD  A  FOR THE FIXED POINT AND ZERO FINDING PROBLEMS.
C
        IF (IFLAGC .GE. -1) THEN
          CALL DCOPY(N,Y(2),1,A,1)
        ENDIF
C
  40    LIMIT=LIMITD
C
C ***** END OF INITIALIZATION BLOCK. *****
C
C ***** MAIN LOOP. *****
C
  50    DO 400 ITER=1,LIMIT
        IF (Y(1) .LT. 0.0) THEN
          ARCLEN = S
          IFLAG = 5
          RETURN
        END IF
C
C TAKE A STEP ALONG THE CURVE.
C
        CALL STEPQF(N,NFE,IFLAGC,START,CRASH,HOLD,H,WK,
     $    RELERR,ABSERR,S,Y,YP,YOLD,YPOLD,A,QT,R,F0,F1,Z0,DZ,
     $    W,T,SSPAR,PAR,IPAR)
C
C PRINT LATEST POINT ON CURVE IF REQUESTED.
C
      IF (TRACE .GT. 0) THEN
        WRITE (TRACE,217) ITER,NFE,S,Y(1),(Y(JW),JW=2,NP1)
217     FORMAT(/' STEP',I5,3X,'NFE =',I5,3X,'ARC LENGTH =',F9.4,3X,
     $  'LAMBDA =',F7.4,5X,'X vector:'/1P,(1X,6E12.4))
      ENDIF
C
C CHECK IF THE STEP WAS SUCCESSFUL.
C
        IF (IFLAGC .GT. 0) THEN
          ARCLEN=S
          IFLAG=IFLAGC
          RETURN
        END IF
C
        IF (CRASH) THEN
C
C         RETURN CODE FOR ERROR TOLERANCE TOO SMALL.
C
          IFLAG=2
C
C         CHANGE ERROR TOLERANCES.
C
          IF (ARCRE .LT. RELERR) THEN
            ARCRE=RELERR
            ANSRE=RELERR
          END IF
          IF (ARCAE .LT. ABSERR) ARCAE=ABSERR
C
C         CHANGE LIMIT ON NUMBER OF ITERATIONS.
C
          LIMIT = LIMIT - ITER
          RETURN
        END IF
C
C IF LAMBDA >= 1.0, USE ROOTQF TO FIND SOLUTION.
C
        IF (Y(1) .GE. 1.0) GOTO 500
C
  400   CONTINUE
C
C ***** END OF MAIN LOOP *****
C
C DID NOT CONVERGE IN  LIMIT  ITERATIONS, SET  IFLAG  AND RETURN.
C
        ARCLEN = S
        IFLAG = 3
        RETURN
C
C ***** FINAL STEP -- FIND SOLUTION AT LAMBDA=1 *****
C
C SAVE  YOLD  FOR ARC LENGTH CALCULATION LATER.
C
  500   CALL DCOPY(NP1,YOLD,1,YSAV,1)
C
C FIND SOLUTION.
C
        CALL ROOTQF(N,NFE,IFLAGC,ANSRE,ANSAE,Y,YP,YOLD,
     $    YPOLD,A,QT,R,DZ,Z0,W,T,F0,F1,PAR,IPAR)
C
C CHECK IF SOLUTION WAS FOUND AND SET  IFLAG  ACCORDINGLY.
C
        IFLAG=1
C
C     SET ERROR FLAG IF ROOTQF COULD NOT GET THE POINT ON THE ZERO
C     CURVE AT  LAMBDA = 1.0.
C
        IF (IFLAGC .GT. 0) IFLAG=IFLAGC
C
C CALCULATE FINAL ARC LENGTH.
C
        CALL DCOPY(NP1,Y,1,DZ,1)
        WK=-1.0
        CALL DAXPY(NP1,WK,YSAV,1,DZ,1)
        ARCLEN = S - HOLD + DNRM2(NP1,DZ,1)
C
C ***** END OF FINAL STEP *****
C
        RETURN
C
C ***** END OF SUBROUTINE FIXPQF *****
        END

max_stars_repo_name: sbwhitecap/clocc-hg

max_issues_repo_name: sbwhitecap/clocc-hg

max_forks_repo_name: sbwhitecap/clocc-hg


text:
      SUBROUTINE UFBSTP(LUNIO,USR,I1,I2,IRET,STR)

C$$$  SUBPROGRAM DOCUMENTATION BLOCK
C
C SUBPROGRAM:    UFBSTP (docblock incomplete)
C   PRGMMR: WOOLLEN          ORG: NP20       DATE: 1999-11-18
C
C ABSTRACT: THIS SUBROUTINE WRITES OR READS SPECIFIED VALUES TO OR FROM
C   THE CURRENT BUFR DATA SUBSET WITHIN INTERNAL ARRAYS, WITH THE
C   DIRECTION OF THE DATA TRANSFER DETERMINED BY THE CONTEXT OF
C   ABS(LUNIO) (I.E., IF ABS(LUNIO) POINTS TO A BUFR FILE THAT IS OPEN
C   FOR INPUT, THEN DATA VALUES ARE READ FROM THE INTERNAL DATA SUBSET;
C   OTHERWISE, DATA VALUES ARE WRITTEN TO THE INTERNAL DATA SUBSET.
C   THE DATA VALUES CORRESPOND TO MNEMONICS WHICH ARE
C   ....
C   IF UFBSTP IS READING VALUES, THEN EITHER BUFR ARCHIVE LIBRARY
C   SUBROUTINE READSB OR READNS MUST HAVE BEEN PREVIOUSLY
C   CALLED TO READ THE SUBSET FROM UNIT ABS(LUNIO) INTO INTERNAL
C   MEMORY.  IF IT IS WRITING VALUES, THEN EITHER BUFR ARCHIVE LIBRARY
C   SUBROUTINE OPENMG OR OPENMB MUST HAVE BEEN PREVIOUSLY CALLED TO
C   OPEN AND INITIALIZE A BUFR MESSAGE WITHIN MEMORY FOR THIS ABS(LUNIO).
C
C PROGRAM HISTORY LOG:
C 1999-11-18  J. WOOLLEN -- ORIGINAL AUTHOR
C 2003-11-04  S. BENDER  -- ADDED REMARKS/BUFRLIB ROUTINE
C                           INTERDEPENDENCIES
C 2003-11-04  D. KEYSER  -- MAXJL (MAXIMUM NUMBER OF JUMP/LINK ENTRIES)
C                           INCREASED FROM 15000 TO 16000 (WAS IN
C                           VERIFICATION VERSION); UNIFIED/PORTABLE FOR
C                           WRF; ADDED DOCUMENTATION (INCLUDING
C                           HISTORY) (INCOMPLETE); OUTPUTS MORE
C                           COMPLETE DIAGNOSTIC INFO WHEN ROUTINE
C                           TERMINATES ABNORMALLY OR UNUSUAL THINGS
C                           HAPPEN; CHANGED CALL FROM BORT TO BORT2 IN
C                           SOME CASES
C 2004-08-18  J. ATOR    -- ADDED SAVE FOR IFIRST1 AND IFIRST2 FLAGS
C
C USAGE:    CALL UFBSTP (LUNIO, USR, I1, I2, IRET, STR)
C   INPUT ARGUMENT LIST:
C     LUNIO    - INTEGER: ABSOLUTE VALUE IS FORTRAN LOGICAL UNIT
C                NUMBER FOR BUFR FILE
C                  - IF BUFR FILE OPEN FOR OUTPUT AND LUNIO IS LESS
C                    THAN ZERO, UFBSTP TREATS THE BUFR FILE AS THOUGH
C                    IT WERE OPEN FOR INPUT
C     USR      - ONLY IF BUFR FILE OPEN FOR OUTPUT:
C                   REAL*8: (I1,I2) STARTING ADDRESS OF DATA VALUES
C                   WRITTEN TO DATA SUBSET
C     I1       - INTEGER: LENGTH OF FIRST DIMENSION OF USR OR THE
C                NUMBER OF BLANK-SEPARATED MNEMONICS IN STR (FORMER
C                MUST BE AT LEAST AS LARGE AS LATTER)
C     I2       - INTEGER:
C                  - IF BUFR FILE OPEN FOR INPUT:  LENGTH OF SECOND
C                    DIMENSION OF USR
C                  - IF BUFR FILE OPEN FOR OUTPUT: NUMBER OF "LEVELS
C                    OF DATA VALUES TO BE WRITTEN TO DATA SUBSET
C                    (MAXIMUM VALUE IS 255)
C     STR      - CHARACTER*(*): STRING ....
C
C   OUTPUT ARGUMENT LIST:
C     USR      - ONLY IF BUFR FILE OPEN FOR INPUT:
C                   REAL*8: (I1,I2) STARTING ADDRESS OF DATA VALUES
C                   READ FROM DATA SUBSET
C     IRET     - INTEGER:
C                  - IF BUFR FILE OPEN FOR INPUT: NUMBER OF "LEVELS" OF
C                    DATA VALUES READ FROM DATA SUBSET (MUST BE NO
C                    LARGER THAN I2)
C                  - IF BUFR FILE OPEN FOR OUTPUT: NUMBER OF "LEVELS"
C                    OF DATA VALUES WRITTEN TO DATA SUBSET (SHOULD BE
C                    SAME AS I2)
C
C   OUTPUT FILES:
C     UNIT 06  - STANDARD OUTPUT PRINT
C
C REMARKS:
C    THIS ROUTINE CALLS:        BORT     BORT2    STATUS   STRING
C                               UFBSP
C    THIS ROUTINE IS CALLED BY: None
C                               Normally called only by application
C                               programs.
C
C ATTRIBUTES:
C   LANGUAGE: FORTRAN 77
C   MACHINE:  PORTABLE TO ALL PLATFORMS
C
C$$$

      INCLUDE 'bufrlib.prm'

      COMMON /MSGCWD/ NMSG(NFILES),NSUB(NFILES),MSUB(NFILES),
     .                INODE(NFILES),IDATE(NFILES)
      COMMON /USRINT/ NVAL(NFILES),INV(MAXJL,NFILES),VAL(MAXJL,NFILES)
      COMMON /QUIET / IPRT

      CHARACTER*(*) STR
      CHARACTER*128 BORT_STR1,BORT_STR2
      REAL*8        USR(I1,I2),VAL

      DATA IFIRST1/0/,IFIRST2/0/

      SAVE IFIRST1, IFIRST2

C----------------------------------------------------------------------
C----------------------------------------------------------------------

      IRET = 0

C  CHECK THE FILE STATUS AND I-NODE
C  --------------------------------

      LUNIT = ABS(LUNIO)
      CALL STATUS(LUNIT,LUN,IL,IM)
      IF(IL.EQ.0) GOTO 900
      IF(IM.EQ.0) GOTO 901
      IF(INODE(LUN).NE.INV(1,LUN)) GOTO 902

      IO = MIN(MAX(0,IL),1)
      IF(LUNIO.NE.LUNIT) IO = 0

      IF(I1.LE.0) THEN
         IF(IPRT.GE.0) THEN
      PRINT*
      PRINT*,'+++++++++++++++++++++++WARNING+++++++++++++++++++++++++'
         PRINT*,'BUFRLIB: UFBSTP - THIRD ARGUMENT (INPUT) IS .LE. 0',
     .    ' -  RETURN WITH FIFTH ARGUMENT (IRET) = 0'
         PRINT*,'STR = ',STR
      PRINT*,'+++++++++++++++++++++++WARNING+++++++++++++++++++++++++'
      PRINT*
         ENDIF
         GOTO 100
      ELSEIF(I2.LE.0) THEN
         IF(IPRT.EQ.-1)  IFIRST1 = 1
         IF(IO.EQ.0 .OR. IFIRST1.EQ.0 .OR. IPRT.GE.1)  THEN
      PRINT*
      PRINT*,'+++++++++++++++++++++++WARNING+++++++++++++++++++++++++'
            PRINT*,'BUFRLIB: UFBSTP - FOURTH ARGUMENT (INPUT) IS .LE. ',
     .       '0 -  RETURN WITH FIFTH ARGUMENT (IRET) = 0'
            PRINT*,'STR = ',STR
            IF(IPRT.EQ.0 .AND. IO.EQ.1)  PRINT 101
101   FORMAT('Note: Only the first occurrence of this WARNING message ',
     . 'is printed, there may be more.  To output'/6X,'ALL WARNING ',
     . 'messages, modify your application program to add ',
     . '"CALL OPENBF(0,''QUIET'',1)" prior'/6X,'to the first call to a',
     . ' BUFRLIB routine.')
      PRINT*,'+++++++++++++++++++++++WARNING+++++++++++++++++++++++++'
      PRINT*
            IFIRST1 = 1
         ENDIF
         GOTO 100
      ENDIF

C  INITIALIZE USR ARRAY PRECEEDING AN INPUT OPERATION
C  --------------------------------------------------

      IF(IO.EQ.0) THEN
         DO J=1,I2
         DO I=1,I1
         USR(I,J) = BMISS
         ENDDO
         ENDDO
      ENDIF

C  PARSE OR RECALL THE INPUT STRING - READ/WRITE VALUES
C  ----------------------------------------------------

      CALL STRING(STR,LUN,I1,IO)

C  CALL THE MNEMONIC READER/WRITER
C  -------------------------------

      CALL UFBSP(LUN,USR,I1,I2,IO,IRET)

      IF(IO.EQ.1 .AND. IRET.NE.I2) GOTO 903

      IF(IRET.EQ.0)  THEN
         IF(IO.EQ.0) THEN
            IF(IPRT.GE.1)  THEN
      PRINT*
      PRINT*,'+++++++++++++++++++++++WARNING+++++++++++++++++++++++++'
               PRINT*,'BUFRLIB: UFBSTP - NO SPECIFIED VALUES READ IN',
     .          ' -  RETURN WITH FIFTH ARGUMENT (IRET) = 0'
               PRINT*,'STR = ',STR
      PRINT*,'+++++++++++++++++++++++WARNING+++++++++++++++++++++++++'
      PRINT*
            ENDIF
         ELSE
            IF(IPRT.EQ.-1)  IFIRST2 = 1
            IF(IFIRST2.EQ.0 .OR. IPRT.GE.1)  THEN
      PRINT*
      PRINT*,'+++++++++++++++++++++++WARNING+++++++++++++++++++++++++'
               PRINT*,'BUFRLIB: UFBSTP - NO SPECIFIED VALUES WRITTEN ',
     .          'OUT -  RETURN WITH FIFTH ARGUMENT (IRET) = 0'
               PRINT*,'STR = ',STR,' MAY NOT BE IN THE BUFR TABLE(?)'
               IF(IPRT.EQ.0)  PRINT 101
      PRINT*,'+++++++++++++++++++++++WARNING+++++++++++++++++++++++++'
      PRINT*
               IFIRST2 = 1
            ENDIF
         ENDIF
      ENDIF

C  EXITS
C  -----

100   RETURN
900   CALL BORT('BUFRLIB: UFBSTP - BUFR FILE IS CLOSED, IT MUST BE'//
     . ' OPEN')
901   CALL BORT('BUFRLIB: UFBSTP - A MESSAGE MUST BE OPEN IN BUFR '//
     . 'FILE, NONE ARE')
902   CALL BORT('BUFRLIB: UFBSTP - LOCATION OF INTERNAL TABLE FOR '//
     . 'BUFR FILE DOES NOT AGREE WITH EXPECTED LOCATION IN INTERNAL '//
     . 'SUBSET ARRAY')
903   WRITE(BORT_STR1,'("BUFRLIB: UFBSTP - MNEMONIC STRING READ IN IS'//
     . ': ",A)') STR
      WRITE(BORT_STR2,'(18X,"THE NUMBER OF ''LEVELS'' ACTUALLY '//
     . 'WRITTEN (",I3,") DOES NOT EQUAL THE NUMBER REQUESTED (",I3,")'//
     . ' - INCOMPLETE WRITE")')  IRET,I2
      CALL BORT2(BORT_STR1,BORT_STR2)
      END


! Erick Juarez
! CPSC 479 Sec 1
! HOMEWORK 6 - 4/20/20
! tested using GNU Fortran (GCC) 4.8.5
program fortran_matrix
implicit none

! Initialize 3 matrices with 4 rows & 4 columns
integer, dimension(4, 4) :: A, B, C
integer:: i, j

A = 4 ! set values of A to 4
do i = 1, 4 ! set values of B(i,j) = i + j + 1
    do j = 1, 4
        B(i, j) = i + j + 1
    end do
end do

C = A + B ! set values of C

! Change values of 8 to 16 in C
where(C == 8)
        C = 16
end where

! print results
do i = 1,4
    print '(a, i2, a, i2, 3i5)', 'A (row', i, '): ', A(i, :)
end do

do i = 1,4
    print '(a, i2, a, i2, 3i5)', 'B (row', i, '): ', B(i, :)
end do

do i = 1,4
    print '(a, i2, a, i2, 3i5)', 'C (row', i, '): ', C(i, :)
end do

end