text:
    ! barth jespersen on cut cells, recentering on whole cells
      subroutine limiter_a(qp, qx, qy, mitot, mjtot, irr, nvar, hx, hy,
     . lwidth, lstgrd, xlow, ylow)
       implicit double precision(a-h,o-z)
       dimension qp(mitot,mjtot,4),qx(mitot,mjtot,4),
     &           qy(mitot,mjtot,4),irr(mitot,mjtot)

       common /order2/ ssw, primitive, nolimiter
       dimension nlist(25,2)
       logical primitive
       logical quadrecon
       include "./reconinfo.i"
       include "cirr.i"

      logical  LP_success, output_debug, tricell
      double precision neighb_x(8),neighb_y(8), neighb_u(8)
      double precision Phi_x, Phi_y, LP_Dx,LP_Dy, LS_Dx, LS_Dy
      integer row_A
      logical do_pos_constr ! are we really a cutcell?
      integer my_iters, tot_iters, max_iters   ! maximum no of iterations needed in Simplex
      double precision avg_iters    ! avg number of iterations needed in Simplex
      integer count_cells
      logical   prflag, enufNbor, outside, almost_outside
      logical dump/.false./

       dimension dumax(4),dumin(4),phimin(4)
       dimension graddot(4),alpha(4),recon(4)
       logical IS_GHOST, verbose/.true./
       IS_GHOST(i,j) = (i .le. lwidth .or. i .gt. mitot-lwidth .or.
     .                  j .le. lwidth .or. j .gt. mjtot-lwidth)

      outside(i,j) = ((i .lt. 1) .or. (j.lt. 1) .or.
     .                (i .gt. mitot) .or. (j.gt. mjtot))

       if(primitive .eqv. .false.) then
        print *, "PROBLEM: THE CURRENT LIMITER IMPLEMENTATION ONLY WORKS
     .            WHEN RECONSTRUCTING IN PRIMITIVE VARIABLES."
       endif

      ! quadratic limiting does not work right now.
      quadrecon = .false.

      ! limit using the MC limiter everywhere that you can
      do 50 i = 1, mitot
      do 50 j = 1, mjtot
        k = irr(i,j)
        if(k .eq. lstgrd .and. .not. outside(i+1,j)) then
        if(irr(i+1,j) .eq. lstgrd) then

            ! simple MC limiter
            do m = 1,nvar
                dfx = 2.d0*( qp(i+1,j,m) - qp(i,j,m) ) / hx
                qx(i,j,m) = dminmod(qx(i,j,m) ,dfx)
            enddo
        endif
        endif

        if(k .eq. lstgrd .and. .not. outside(i-1,j) ) then
        if( irr(i-1,j) .eq. lstgrd) then
            ! simple MC limiter
            do m = 1,nvar
                dbx = 2.d0*( qp(i,j,m) - qp(i-1,j,m) ) / hx
                qx(i,j,m) = dminmod(dbx, qx(i,j,m) )
            enddo
        endif
        endif

        if(k .eq. lstgrd .and. .not. outside(i,j+1)) then
        if( irr(i,j+1) .eq. lstgrd ) then

            ! simple MC limiter
            do m = 1,nvar
                dfy = 2.d0*( qp(i,j+1,m) - qp(i,j,m) ) / hy
                qy(i,j,m) = dminmod(qy(i,j,m) ,dfy)
            enddo
        endif
        endif

        if(k .eq. lstgrd .and. .not. outside(i,j-1) ) then
        if( irr(i,j-1) .eq. lstgrd ) then

            ! simple MC limiter
            do m = 1,nvar
                dby = 2.d0*( qp(i,j,m) - qp(i,j-1,m) ) / hy
                qy(i,j,m) = dminmod(dby, qy(i,j,m) )
            enddo
        endif
        endif

50    continue

      ! I still must
      ! limit the cut cells
      ! AND
      ! whole cells adjacent the cut cells

c
c  **************************************
c
c   LIMIT CUT CELLS USING LINEARITY PRESERVING BARTH-JESPERSEN OVER RECONSTRUCTION NEIGHBORHOOD
c

      do 30 j = 1, mjtot
      do 30 i = 1, mitot
          k = irr(i,j)
          if (k .eq. -1 .or. k .eq. lstgrd) goto 30 ! solid/full cells have no need for a gradient

           call getCellCentroid(lstgrd,i,j,xcent,ycent,xlow,
     .                            ylow,hx,hy,k)

          ! find max and min needed for BJ limiting
          dumax = 0.d0
          dumin = 0.d0
          do 31 joff = -jjr(i,j),jjr(i,j)
          do 31 ioff = -iir(i,j),iir(i,j)
!            if (ioff .eq. 0 .and. joff .eq. 0) go to 31
            if (outside(i+ioff,j+joff)) go to 31

            koff = irr(i+ioff,j+joff)
            if (koff .eq. -1) go to 31
            dumax = max(dumax,qp(i+ioff,j+joff,:)-qp(i,j,:))
            dumin = min(dumin,qp(i+ioff,j+joff,:)-qp(i,j,:))
  31      continue

          phimin = 1.d0
          do 32 joff = -jjr(i,j), jjr(i,j)
          do 32 ioff = -iir(i,j), iir(i,j)
!             if (ioff .eq. 0 .and. joff .eq. 0) go to 32 ! no eqn to solve
             if (outside(i+ioff,j+joff)) go to 32
             koff = irr(i+ioff,j+joff)
             if (koff .eq. -1) go to 32
             call getCellCentroid(lstgrd,i+ioff,j+joff,xc,yc,xlow,
     .                            ylow,hx,hy,koff)
              diffx = xc-xcent
              diffy = yc-ycent
              graddot  = qx(i,j,:)*diffx + qy(i,j,:)*diffy
              recon = qp(i,j,:) + graddot
                do m = 1,4
                   if (graddot(m) > 1.d-10) then
                      alpha(m) = min(1.d0, dumax(m)/graddot(m))
                   else if (graddot(m) < -1.d-10) then
                      alpha(m) = min(1.d0, dumin(m)/graddot(m))
                   else
                      alpha(m) = 1.d0
                   endif
                end do
                if (recon(1) .le. 0.d0) alpha(1) = 0.d0
                if (recon(4) .le. 0.d0) alpha(4) = 0.d0
                phimin = min(phimin, alpha)
  32      continue
          qx(i,j,:) = qx(i,j,:)*phimin(:)
          qy(i,j,:) = qy(i,j,:)*phimin(:)

      ! positivity at quadrature points is required, midpoints because p is always = 1.
         do 510 kside=1,6

               x1 = poly(kside,1,k)
               y1 = poly(kside,2,k)
               x2 = poly(kside+1,1,k)
               y2 = poly(kside+1,2,k)

               d_bxp = (x1+x2)/2.d0
               d_byp = (y1+y2)/2.d0

               diffx = d_bxp-xcent
               diffy = d_byp-ycent
               graddot  = qx(i,j,:)*diffx + qy(i,j,:)*diffy
               recon = qp(i,j,:) + graddot

              alpha = 1.d0
              if (recon(1) .le. 0.d0) alpha(1) = 0.d0
              if (recon(4) .le. 0.d0) alpha(4) = 0.d0

!                if (alpha(1) < 0.5d0) then
!                print *, "here"
!                endif

              qx(i,j,1) = qx(i,j,1)*alpha(1)
              qy(i,j,4) = qy(i,j,4)*alpha(4)
              if (poly(kside+2,1,k).eq.-11.) exit
 510      continue

!      alpha = 1.d0
!      if (recon(1) .le. 0.d0 .and. graddot(1) > 0) then
!      alpha(1) = (d_eps-qp(i,j,1)) / graddot(1)
!      alpha(1) = min(alpha(1), 1.d0)
!      endif
!
!      if (recon(4) .le. 0.d0 .and. graddot(4) > 0) then
!      alpha(4) = (d_eps-qp(i,j,4)) / graddot(4)
!      alpha(4) = min(alpha(4), 1.d0)
!      endif
!
!                if (alpha(1) < 1.d0) then
!                print *, "here"
!                endif
!
!              qx(i,j,1) = qx(i,j,1)*alpha(1)
!              qy(i,j,4) = qy(i,j,4)*alpha(4)

  30  continue

c
c  **************************************
c
c   LIMIT FLOW CELLS ADJACENT TO CUTS USING RECENTERING
c
      do 210 ix0 = 1, mitot
      do 210 iy0 = 1, mjtot
        k = irr(ix0,iy0)
        if (k .eq. -1) go to 210
        if (k .ne. lstgrd) go to 210

            x0 = xlow + (ix0-.5d0)*hx
            y0 = ylow + (iy0-.5d0)*hy

            do m = 1, nvar
c
c  for each irregular side, use recentering. adjacent cant be solid since this cell regular
       if(.not. outside(ix0+1,iy0)) then
       if (irr(ix0+1,iy0) .ne. lstgrd) then
            kp1 = irr(ix0+1,iy0)
            xcut = xcirr(kp1)
            ycut = ycirr(kp1)
            recentM =   qp(ix0+1,iy0,m) + (y0-ycut)*qy(ix0+1,iy0,m)
            cutmin= (recentM - qp(ix0,iy0,m))/(xcut-x0)

             othermin = qx(ix0,iy0,m)
             qx(ix0,iy0,m) = .5*min(abs(cutmin),abs(othermin)) *
     &                         (sign(1.d0,othermin)+sign(1.d0,cutmin))
!            endif
       endif
       endif
       if(.not. outside(ix0-1,iy0)) then
       if (irr(ix0-1,iy0) .ne. lstgrd) then
            km1 = irr(ix0-1,iy0)
            xcut = xcirr(km1)
            ycut = ycirr(km1)
            recentM =   qp(ix0-1,iy0,m) + (y0-ycut)*qy(ix0-1,iy0,m)
            cutmin= (recentM - qp(ix0,iy0,m))/(xcut-x0)

            othermin = qx(ix0,iy0,m)
            qx(ix0,iy0,m) = .5*min(abs(cutmin),abs(othermin))*
     &                     (sign(1.d0,othermin)+sign(1.d0,cutmin))
       endif
       endif
c  for each irregular side, use recentering.
       if(.not. outside(ix0,iy0+1)) then
       if (irr(ix0,iy0+1) .ne. lstgrd) then
            kp1 = irr(ix0,iy0+1)
            xcut = xcirr(kp1)
            ycut = ycirr(kp1)
            recentM =   qp(ix0,iy0+1,m) + (x0-xcut)*qx(ix0,iy0+1,m)
            cutmin= (recentM - qp(ix0,iy0,m))/(ycut-y0)

            othermin = qy(ix0,iy0,m)
            qy(ix0,iy0,m) = .5*min(abs(cutmin),abs(othermin))*
     &                          (sign(1.d0,othermin)+sign(1.d0,cutmin))
!            endif
       endif
       endif
       if(.not. outside(ix0,iy0-1)) then
       if (irr(ix0,iy0-1) .ne. lstgrd) then
            km1 = irr(ix0,iy0-1)
            xcut = xcirr(km1)
            ycut = ycirr(km1)
            recentM =   qp(ix0,iy0-1,m) + (x0-xcut)*qx(ix0,iy0-1,m)
            cutmin= (recentM - qp(ix0,iy0,m))/(ycut-y0)

            othermin = qy(ix0,iy0,m)
            qy(ix0,iy0,m) = .5* min(abs(cutmin),abs(othermin))*
     &                      (sign(1.d0,othermin)+sign(1.d0,cutmin))

       endif
       endif

       end do
 210  continue

      return
      end

      function dminmod(a,b)
      implicit double precision(a-h,o-z)
      if( a > 0 .and. b > 0 ) then
        dminmod = min(a,b)
      elseif( a < 0 .and. b < 0 ) then
        dminmod = max(a,b)
      else
        dminmod = 0.d0
      endif
      end

! barth jespersen on cut cells and on whole cells
      subroutine limiter(qp, qx, qy, mitot, mjtot, irr, nvar, hx, hy,
     . lwidth, lstgrd, xlow, ylow)
       implicit double precision(a-h,o-z)
       dimension qp(mitot,mjtot,4),qx(mitot,mjtot,4),
     &           qy(mitot,mjtot,4),irr(mitot,mjtot)

       common /order2/ ssw, primitive, nolimiter
       dimension nlist(25,2)
       logical primitive
       logical quadrecon
       include "./reconinfo.i"
       include "cirr.i"

      logical  LP_success, output_debug, tricell
      double precision neighb_x(8),neighb_y(8), neighb_u(8)
      double precision Phi_x, Phi_y, LP_Dx,LP_Dy, LS_Dx, LS_Dy
      integer row_A
      logical do_pos_constr ! are we really a cutcell?
      integer my_iters, tot_iters, max_iters   ! maximum no of iterations needed in Simplex
      double precision avg_iters    ! avg number of iterations needed in Simplex
      integer count_cells
      logical   prflag, enufNbor, outside, almost_outside
      logical dump/.false./

       dimension dumax(4),dumin(4),phimin(4)
       dimension graddot(4),alpha(4),recon(4)
       logical IS_GHOST, verbose/.true./
       IS_GHOST(i,j) = (i .le. lwidth .or. i .gt. mitot-lwidth .or.
     .                  j .le. lwidth .or. j .gt. mjtot-lwidth)

      outside(i,j) = ((i .lt. 1) .or. (j.lt. 1) .or.
     .                (i .gt. mitot) .or. (j.gt. mjtot))

       if(primitive .eqv. .false.) then
        print *, "PROBLEM: THE CURRENT LIMITER IMPLEMENTATION ONLY WORKS
     .            WHEN RECONSTRUCTING IN PRIMITIVE VARIABLES."
       endif

      ! quadratic limiting does not work right now.
      quadrecon = .false.

      ! limit using the MC limiter everywhere that you can
      do 50 i = 1, mitot
      do 50 j = 1, mjtot

        k = irr(i,j)
        if ((k .eq. -1) .or. (inuf(i,j) .eq. 0)) cycle

        do m = 1,nvar
           dfx = 2.d0*( qp(i+1,j,m) - qp(i,j,m) ) / hx
           dbx = 2.d0*( qp(i,j,m) - qp(i-1,j,m) ) / hx
           dfy = 2.d0*( qp(i,j+1,m) - qp(i,j,m) ) / hy
           dby = 2.d0*( qp(i,j,m) - qp(i,j-1,m) ) / hy

           qx(i,j,m) = dminmod(qx(i,j,m) ,dfx)
           qx(i,j,m) = dminmod(dbx, qx(i,j,m) )
           qy(i,j,m) = dminmod(qy(i,j,m) ,dfy)
           qy(i,j,m) = dminmod(dby, qy(i,j,m) )
        enddo

50    continue
c
c  **************************************
c
c   LIMIT CUT CELLS USING LINEARITY PRESERVING BARTH-JESPERSEN OVER RECONSTRUCTION NEIGHBORHOOD
c

      do 30 j = 1, mjtot
      do 30 i = 1, mitot
          k = irr(i,j)
          if ((k .eq. -1) .or. (inuf(i,j) .eq. 1)) cycle ! skips solid/regular cells

           call getCellCentroid(lstgrd,i,j,xcent,ycent,xlow,
     .                            ylow,hx,hy,k)

          ! find max and min needed for BJ limiting
          dumax = 0.d0
          dumin = 0.d0
          do 31 joff = -jjr(i,j),jjr(i,j)
          do 31 ioff = -iir(i,j),iir(i,j)
!            if (ioff .eq. 0 .and. joff .eq. 0) go to 31
            if (outside(i+ioff,j+joff)) go to 31

            koff = irr(i+ioff,j+joff)
            if (koff .eq. -1) go to 31
            dumax = max(dumax,qp(i+ioff,j+joff,:)-qp(i,j,:))
            dumin = min(dumin,qp(i+ioff,j+joff,:)-qp(i,j,:))
  31      continue

          phimin = 1.d0
          do 32 joff = -jjr(i,j), jjr(i,j)
          do 32 ioff = -iir(i,j), iir(i,j)
!             if (ioff .eq. 0 .and. joff .eq. 0) go to 32 ! no eqn to solve
             if (outside(i+ioff,j+joff)) go to 32
             koff = irr(i+ioff,j+joff)
             if (koff .eq. -1) go to 32
             call getCellCentroid(lstgrd,i+ioff,j+joff,xc,yc,xlow,
     .                            ylow,hx,hy,koff)
              diffx = xc-xcent
              diffy = yc-ycent
              graddot  = qx(i,j,:)*diffx + qy(i,j,:)*diffy
              recon = qp(i,j,:) + graddot
                do m = 1,4
                   if (graddot(m) > 1.d-10) then
                      alpha(m) = min(1.d0, dumax(m)/graddot(m))
                   else if (graddot(m) < -1.d-10) then
                      alpha(m) = min(1.d0, dumin(m)/graddot(m))
                   else
                      alpha(m) = 1.d0
                   endif
                end do
                if (recon(1) .le. 0.d0) alpha(1) = 0.d0
                if (recon(4) .le. 0.d0) alpha(4) = 0.d0
                phimin = min(phimin, alpha)
  32      continue
          qx(i,j,:) = qx(i,j,:)*phimin(:)
          qy(i,j,:) = qy(i,j,:)*phimin(:)

      ! check for positivity
      if(irr(i,j) .eq. lstgrd) then
         alpha = 1.d0

         diffx = hx/2.d0
         diffy = 0.
         graddot  = qx(i,j,:)*diffx + qy(i,j,:)*diffy
         recon = qp(i,j,:) + graddot

         if (recon(1) .le. 0.d0) alpha(1) = 0.d0
         if (recon(4) .le. 0.d0) alpha(4) = 0.d0

         diffx = -hx/2.d0
         diffy = 0.
         graddot  = qx(i,j,:)*diffx + qy(i,j,:)*diffy
         recon = qp(i,j,:) + graddot

         if (recon(1) .le. 0.d0) alpha(1) = 0.d0
         if (recon(4) .le. 0.d0) alpha(4) = 0.d0

         diffx = 0.d0
         diffy = hy/2.d0
         graddot  = qx(i,j,:)*diffx + qy(i,j,:)*diffy
         recon = qp(i,j,:) + graddot

         if (recon(1) .le. 0.d0) alpha(1) = 0.d0
         if (recon(4) .le. 0.d0) alpha(4) = 0.d0

         diffx = 0.d0
         diffy = -hy/2.d0
         graddot  = qx(i,j,:)*diffx + qy(i,j,:)*diffy
         recon = qp(i,j,:) + graddot

         if (recon(1) .le. 0.d0) alpha(1) = 0.d0
         if (recon(4) .le. 0.d0) alpha(4) = 0.d0

         qx(i,j,1) = qx(i,j,1)*alpha(1)
         qy(i,j,4) = qy(i,j,4)*alpha(4)
      else
      ! positivity at quadrature points is required, midpoints because p is always = 1.
         do 510 kside=1,6
               x1 = poly(kside,1,k)
               y1 = poly(kside,2,k)
               x2 = poly(kside+1,1,k)
               y2 = poly(kside+1,2,k)

               d_bxp = (x1+x2)/2.d0
               d_byp = (y1+y2)/2.d0

               diffx = d_bxp-xcent
               diffy = d_byp-ycent
               graddot  = qx(i,j,:)*diffx + qy(i,j,:)*diffy
               recon = qp(i,j,:) + graddot

              alpha = 1.d0
              if (recon(1) .le. 0.d0) alpha(1) = 0.d0
              if (recon(4) .le. 0.d0) alpha(4) = 0.d0

!                if (alpha(1) < 0.5d0) then
!                print *, "here"
!                endif

              qx(i,j,1) = qx(i,j,1)*alpha(1)
              qy(i,j,4) = qy(i,j,4)*alpha(4)
              if (poly(kside+2,1,k).eq.-11.) exit
 510      continue
      endif

  30  continue

      return
      end

        ! apparently recentering + Sandra's limiter does not work

      ! limit the cut cells first
!      do 110 ix0 = 1, mitot
!      do 110 iy0 = 1, mjtot
!         k = irr(ix0,iy0)
!         if (k .eq. -1) go to 110
!         if (k .eq. lstgrd) go to 110
!
!
!         x0 = xcirr(k)
!         y0 = ycirr(k)
!
!
!         do 37 kside = 1, 6
!            if (poly(kside+2,1,k).eq.-11)then
!               sidex2 = poly(kside,1,k)
!               sidey2 = poly(kside,2,k)
!               sidex1 = poly(kside+1,1,k)
!               sidey1 = poly(kside+1,2,k)
!               go to 39
!            endif
! 37      continue
! 39      bxpt = .5d0*(sidex1+sidex2)
!         bypt = .5d0*(sidey1+sidey2)
!
!
!
!
!
!
!
!         nlist(1,1) = ix0
!         nlist(1,2) = iy0
!         nst        = 1
!         nend       = 1
!         call addneighbors(irr,nlist,nst,nend,newend,mitot,mjtot,
!     .        lstgrd,quadrecon,xlow,ylow,hx,hy)
!
!c
!c convert from my routine to Sandras cut cell lp routine
!c
!          center_x = x0
!          center_y = y0
!
!          num_neighb = newend-1  ! not counting cell itself
!          output_debug = .false.
!
!
!          do m = 1, 4   ! loop over the 4 variables
!            row_A     = 2*num_neighb + 4
!            do_pos_constr = .false.
!            if (irr(ix0,iy0) .ne. lstgrd) then
!               if (m .eq. 1 .OR. m .eq. 4) then
!                  do_pos_constr = .true.
!                  row_A = row_A + 1 ! additional pos constraint for rho, p
!               endif
!            endif
!
!            do inbor = 1, num_neighb
!              ixn = nlist(inbor+1,1)
!              iyn = nlist(inbor+1,2)
!              kn =  irr(ixn,iyn)
!              if (kn .ne. lstgrd) then
!                 xn = xcirr(kn)
!                 yn = ycirr(kn)
!              else
!                 xn = xlow + (ixn-.5d0)*hx
!                 yn = ylow + (iyn-.5d0)*hy
!              endif
!
!              neighb_x(inbor) = xn ! neighboring centroid, not diffs
!              neighb_y(inbor) = yn ! neighboring centroid, not diffs
!
!              neighb_u(inbor) = qp(ixn,iyn, m)  ! neighboring u, need to loop over m
!            end do  ! do inbor = 1, num_neighb
!              center_u = qp(ix0,iy0,m)
!              LS_Dx = qx(ix0,iy0,m)
!              LS_Dy = qy(ix0,iy0,m)
!
!            LP_success = .false.
!
!
!
!             call vector_lim_allineq(num_neighb,neighb_x,neighb_y,
!     &               neighb_u,center_x,center_y,center_u,output_debug,
!     &               row_A,Phi_x,Phi_y, LP_Dx,LP_Dy,
!     &         LP_success,LS_Dx,LS_Dy,bxpt,bypt,do_pos_constr,my_iters)
!
!c           Increase the count of number of iterations
!            tot_iters = tot_iters + my_iters
!            if (my_iters .GT. max_iters) then
!               max_iters = my_iters
!            endif
!            count_cells = count_cells + 1
!
!
! 402         format("sandras lsq ",2e15.7,10x," limited ",2e15.7)
!
!             if (LP_success) then
!                  qx(ix0,iy0,m) = LP_Dx
!                  qy(ix0,iy0,m) = LP_Dy
!
! 400              format("limiter values for m = ",i3," are ",2e15.7)
!             else
!                  write(*,*)" LP failure for cell ",ix0,iy0
!c                  stop
!             endif
!
!          end do   ! do m=1,4
!
!
!
!
!
!110   continue
!
!
!c
!c  **************************************
!c
!c   LIMIT FLOW CELLS ADJACENT TO CUTS USING RECENTERING
!c
!      do 210 ix0 = 1, mitot
!      do 210 iy0 = 1, mjtot
!        k = irr(ix0,iy0)
!        if (k .eq. -1) go to 210
!        if (k .ne. lstgrd) go to 210
!
!
!            x0 = xlow + (ix0-.5d0)*hx
!            y0 = ylow + (iy0-.5d0)*hy
!
!            do m = 1, nvar
!c
!c  for each irregular side, use recentering. adjacent cant be solid since this cell regular
!       if(.not. outside(ix0+1,iy0)) then
!       if (irr(ix0+1,iy0) .ne. lstgrd) then
!            kp1 = irr(ix0+1,iy0)
!            xcut = xcirr(kp1)
!            ycut = ycirr(kp1)
!            recentM =   qp(ix0+1,iy0,m) + (y0-ycut)*qy(ix0+1,iy0,m)
!            cutmin= (recentM - qp(ix0,iy0,m))/(xcut-x0)
!
!             othermin = qx(ix0,iy0,m)
!             qx(ix0,iy0,m) = .5*min(abs(cutmin),abs(othermin)) *
!     &                         (sign(1.d0,othermin)+sign(1.d0,cutmin))
!!            endif
!       endif
!       endif
!       if(.not. outside(ix0-1,iy0)) then
!       if (irr(ix0-1,iy0) .ne. lstgrd) then
!            km1 = irr(ix0-1,iy0)
!            xcut = xcirr(km1)
!            ycut = ycirr(km1)
!            recentM =   qp(ix0-1,iy0,m) + (y0-ycut)*qy(ix0-1,iy0,m)
!            cutmin= (recentM - qp(ix0,iy0,m))/(xcut-x0)
!
!            othermin = qx(ix0,iy0,m)
!            qx(ix0,iy0,m) = .5*min(abs(cutmin),abs(othermin))*
!     &                     (sign(1.d0,othermin)+sign(1.d0,cutmin))
!       endif
!       endif
!c  for each irregular side, use recentering.
!       if(.not. outside(ix0,iy0+1)) then
!       if (irr(ix0,iy0+1) .ne. lstgrd) then
!            kp1 = irr(ix0,iy0+1)
!            xcut = xcirr(kp1)
!            ycut = ycirr(kp1)
!            recentM =   qp(ix0,iy0+1,m) + (x0-xcut)*qx(ix0,iy0+1,m)
!            cutmin= (recentM - qp(ix0,iy0,m))/(ycut-y0)
!
!            othermin = qy(ix0,iy0,m)
!            qy(ix0,iy0,m) = .5*min(abs(cutmin),abs(othermin))*
!     &                          (sign(1.d0,othermin)+sign(1.d0,cutmin))
!!            endif
!       endif
!       endif
!       if(.not. outside(ix0,iy0-1)) then
!       if (irr(ix0,iy0-1) .ne. lstgrd) then
!            km1 = irr(ix0,iy0-1)
!            xcut = xcirr(km1)
!            ycut = ycirr(km1)
!            recentM =   qp(ix0,iy0-1,m) + (x0-xcut)*qx(ix0,iy0-1,m)
!            cutmin= (recentM - qp(ix0,iy0,m))/(ycut-y0)
!
!            othermin = qy(ix0,iy0,m)
!            qy(ix0,iy0,m) = .5* min(abs(cutmin),abs(othermin))*
!     &                      (sign(1.d0,othermin)+sign(1.d0,cutmin))
!
!       endif
!       endif
!
!       end do
! 210  continue

ext: f


text:
! Auto-generated. DO NOT MODIFY!
module tensor_permute_loops_module
    use, intrinsic :: iso_fortran_env, only : &
            int64, &
            real64, &
            real32
    use :: data_api, only : &
            stream
    use :: tensor_api, only : &
            tensor, &
            tensor_fortran_converter, &
            matrix, &
            tensor_rank3, &
            tensor_rank4, &
            tensor_rank5, &
            dt_c128, &
            dt_c64, &
            dt_r64, &
            dt_r32
    use :: raw_tensor_permute_module, only : &
            raw_tensor_permute

    implicit none
    private

    public :: tensor_permute_loops

    type, extends(raw_tensor_permute) :: tensor_permute_loops
        type(tensor_fortran_converter) :: converter
    contains
        procedure :: permute_tensor => permute_tensor
        procedure :: permute_complex128_matrix => permute_complex128_matrix
        procedure :: permute_complex128_tensor_rank3 => permute_complex128_tensor_rank3
        procedure :: permute_complex128_tensor_rank4 => permute_complex128_tensor_rank4
        procedure :: permute_complex128_tensor_rank5 => permute_complex128_tensor_rank5
        procedure :: permute_complex64_matrix => permute_complex64_matrix
        procedure :: permute_complex64_tensor_rank3 => permute_complex64_tensor_rank3
        procedure :: permute_complex64_tensor_rank4 => permute_complex64_tensor_rank4
        procedure :: permute_complex64_tensor_rank5 => permute_complex64_tensor_rank5
        procedure :: permute_real64_matrix => permute_real64_matrix
        procedure :: permute_real64_tensor_rank3 => permute_real64_tensor_rank3
        procedure :: permute_real64_tensor_rank4 => permute_real64_tensor_rank4
        procedure :: permute_real64_tensor_rank5 => permute_real64_tensor_rank5
        procedure :: permute_real32_matrix => permute_real32_matrix
        procedure :: permute_real32_tensor_rank3 => permute_real32_tensor_rank3
        procedure :: permute_real32_tensor_rank4 => permute_real32_tensor_rank4
        procedure :: permute_real32_tensor_rank5 => permute_real32_tensor_rank5
        procedure :: cleanup => cleanup
        procedure :: synchronize => synchronize
    end type tensor_permute_loops

    interface tensor_permute_loops
        module procedure constructor
    end interface tensor_permute_loops
contains
    function constructor() result(this)
        type(tensor_permute_loops) :: this

        this%converter = tensor_fortran_converter()
    end function constructor

    subroutine permute_tensor(this, src, dst, perm, astream)
        class(tensor_permute_loops), intent(inout) :: this
        class(tensor), intent(in) :: src
        class(tensor), intent(inout) :: dst
        integer, dimension(:), intent(in) :: perm
        type(stream), intent(in), optional :: astream

        integer :: rank

        if ( .not. this%is_compatible(src, dst, perm) ) &
                error stop "tensor_permute_loops:Tensors are not compatible."

        select type(src)
        type is (matrix)
            select type(dst)
            type is (matrix)
                rank = 2
                select case ( dst%get_datatype() )
                case(dt_c128)
                    call this%permute_complex128_matrix(src, dst, perm(:2), astream)
                case(dt_c64)
                    call this%permute_complex64_matrix(src, dst, perm(:2), astream)
                case(dt_r64)
                    call this%permute_real64_matrix(src, dst, perm(:2), astream)
                case(dt_r32)
                    call this%permute_real32_matrix(src, dst, perm(:2), astream)
                end select
            end select
        type is (tensor_rank3)
            select type(dst)
            type is (tensor_rank3)
                rank = 3
                select case ( dst%get_datatype() )
                case(dt_c128)
                    call this%permute_complex128_tensor_rank3(src, dst, perm(:3), astream)
                case(dt_c64)
                    call this%permute_complex64_tensor_rank3(src, dst, perm(:3), astream)
                case(dt_r64)
                    call this%permute_real64_tensor_rank3(src, dst, perm(:3), astream)
                case(dt_r32)
                    call this%permute_real32_tensor_rank3(src, dst, perm(:3), astream)
                end select
            end select
        type is (tensor_rank4)
            select type(dst)
            type is (tensor_rank4)
                rank = 4
                select case ( dst%get_datatype() )
                case(dt_c128)
                    call this%permute_complex128_tensor_rank4(src, dst, perm(:4), astream)
                case(dt_c64)
                    call this%permute_complex64_tensor_rank4(src, dst, perm(:4), astream)
                case(dt_r64)
                    call this%permute_real64_tensor_rank4(src, dst, perm(:4), astream)
                case(dt_r32)
                    call this%permute_real32_tensor_rank4(src, dst, perm(:4), astream)
                end select
            end select
        type is (tensor_rank5)
            select type(dst)
            type is (tensor_rank5)
                rank = 5
                select case ( dst%get_datatype() )
                case(dt_c128)
                    call this%permute_complex128_tensor_rank5(src, dst, perm(:5), astream)
                case(dt_c64)
                    call this%permute_complex64_tensor_rank5(src, dst, perm(:5), astream)
                case(dt_r64)
                    call this%permute_real64_tensor_rank5(src, dst, perm(:5), astream)
                case(dt_r32)
                    call this%permute_real32_tensor_rank5(src, dst, perm(:5), astream)
                end select
            end select
        end select
    end subroutine permute_tensor

    subroutine permute_complex128_matrix(this, src, dst, perm, astream)
        class(tensor_permute_loops), intent(inout) :: this
        type(matrix), intent(in) :: src
        type(matrix), intent(inout) :: dst
        integer, dimension(2), intent(in) :: perm
        type(stream), intent(in), optional :: astream

        integer(int64) :: d1, d2
        complex(real64), dimension(:,:), pointer, contiguous :: tsrc, tdst
        integer(int64) :: sidx(2), didx(2)

        associate ( sdims => src%dims, ddims => dst%dims )
            call this%converter%secure_fortran_pointer(tsrc, src, astream)
            call this%converter%secure_fortran_pointer(tdst, dst, astream)

            do d2 = 1, sdims(2)
                sidx(2) = d2
                do d1 = 1, sdims(1)
                    sidx(1) = d1

                    didx = sidx(perm)
                    tdst(didx(1),didx(2)) = tsrc(sidx(1),sidx(2))
                end do
            end do

            call this%converter%update_remote_and_release_pointer(tdst, dst, astream)
            call this%converter%release_pointer(tsrc, src, astream)
        end associate
    end subroutine permute_complex128_matrix

    subroutine permute_complex128_tensor_rank3(this, src, dst, perm, astream)
        class(tensor_permute_loops), intent(inout) :: this
        type(tensor_rank3), intent(in) :: src
        type(tensor_rank3), intent(inout) :: dst
        integer, dimension(3), intent(in) :: perm
        type(stream), intent(in), optional :: astream

        integer(int64) :: d1, d2, d3
        complex(real64), dimension(:,:,:), pointer, contiguous :: tsrc, tdst
        integer(int64) :: sidx(3), didx(3)

        associate ( sdims => src%dims, ddims => dst%dims )
            call this%converter%secure_fortran_pointer(tsrc, src, astream)
            call this%converter%secure_fortran_pointer(tdst, dst, astream)

            do d3 = 1, sdims(3)
                sidx(3) = d3
                do d2 = 1, sdims(2)
                    sidx(2) = d2
                    do d1 = 1, sdims(1)
                        sidx(1) = d1

                        didx = sidx(perm)
                        tdst(didx(1),didx(2),didx(3)) = tsrc(sidx(1),sidx(2),sidx(3))
                    end do
                end do
            end do

            call this%converter%update_remote_and_release_pointer(tdst, dst, astream)
            call this%converter%release_pointer(tsrc, src, astream)
        end associate
    end subroutine permute_complex128_tensor_rank3

    subroutine permute_complex128_tensor_rank4(this, src, dst, perm, astream)
        class(tensor_permute_loops), intent(inout) :: this
        type(tensor_rank4), intent(in) :: src
        type(tensor_rank4), intent(inout) :: dst
        integer, dimension(4), intent(in) :: perm
        type(stream), intent(in), optional :: astream

        integer(int64) :: d1, d2, d3, d4
        complex(real64), dimension(:,:,:,:), pointer, contiguous :: tsrc, tdst
        integer(int64) :: sidx(4), didx(4)

        associate ( sdims => src%dims, ddims => dst%dims )
            call this%converter%secure_fortran_pointer(tsrc, src, astream)
            call this%converter%secure_fortran_pointer(tdst, dst, astream)

            do d4 = 1, sdims(4)
                sidx(4) = d4
                do d3 = 1, sdims(3)
                    sidx(3) = d3
                    do d2 = 1, sdims(2)
                        sidx(2) = d2
                        do d1 = 1, sdims(1)
                            sidx(1) = d1

                            didx = sidx(perm)
                            tdst(didx(1),didx(2),didx(3),didx(4)) = tsrc(sidx(1),sidx(2),sidx(3),sidx(4))
                        end do
                    end do
                end do
            end do

            call this%converter%update_remote_and_release_pointer(tdst, dst, astream)
            call this%converter%release_pointer(tsrc, src, astream)
        end associate
    end subroutine permute_complex128_tensor_rank4

    subroutine permute_complex128_tensor_rank5(this, src, dst, perm, astream)
        class(tensor_permute_loops), intent(inout) :: this
        type(tensor_rank5), intent(in) :: src
        type(tensor_rank5), intent(inout) :: dst
        integer, dimension(5), intent(in) :: perm
        type(stream), intent(in), optional :: astream

        integer(int64) :: d1, d2, d3, d4, d5
        complex(real64), dimension(:,:,:,:,:), pointer, contiguous :: tsrc, tdst
        integer(int64) :: sidx(5), didx(5)

        associate ( sdims => src%dims, ddims => dst%dims )
            call this%converter%secure_fortran_pointer(tsrc, src, astream)
            call this%converter%secure_fortran_pointer(tdst, dst, astream)

            do d5 = 1, sdims(5)
                sidx(5) = d5
                do d4 = 1, sdims(4)
                    sidx(4) = d4
                    do d3 = 1, sdims(3)
                        sidx(3) = d3
                        do d2 = 1, sdims(2)
                            sidx(2) = d2
                            do d1 = 1, sdims(1)
                                sidx(1) = d1

                                didx = sidx(perm)
                                tdst(didx(1),didx(2),didx(3),didx(4),didx(5)) = tsrc(sidx(1),sidx(2),sidx(3),sidx(4),sidx(5))
                            end do
                        end do
                    end do
                end do
            end do

            call this%converter%update_remote_and_release_pointer(tdst, dst, astream)
            call this%converter%release_pointer(tsrc, src, astream)
        end associate
    end subroutine permute_complex128_tensor_rank5

    subroutine permute_complex64_matrix(this, src, dst, perm, astream)
        class(tensor_permute_loops), intent(inout) :: this
        type(matrix), intent(in) :: src
        type(matrix), intent(inout) :: dst
        integer, dimension(2), intent(in) :: perm
        type(stream), intent(in), optional :: astream

        integer(int64) :: d1, d2
        complex(real32), dimension(:,:), pointer, contiguous :: tsrc, tdst
        integer(int64) :: sidx(2), didx(2)

        associate ( sdims => src%dims, ddims => dst%dims )
            call this%converter%secure_fortran_pointer(tsrc, src, astream)
            call this%converter%secure_fortran_pointer(tdst, dst, astream)

            do d2 = 1, sdims(2)
                sidx(2) = d2
                do d1 = 1, sdims(1)
                    sidx(1) = d1

                    didx = sidx(perm)
                    tdst(didx(1),didx(2)) = tsrc(sidx(1),sidx(2))
                end do
            end do

            call this%converter%update_remote_and_release_pointer(tdst, dst, astream)
            call this%converter%release_pointer(tsrc, src, astream)
        end associate
    end subroutine permute_complex64_matrix

    subroutine permute_complex64_tensor_rank3(this, src, dst, perm, astream)
        class(tensor_permute_loops), intent(inout) :: this
        type(tensor_rank3), intent(in) :: src
        type(tensor_rank3), intent(inout) :: dst
        integer, dimension(3), intent(in) :: perm
        type(stream), intent(in), optional :: astream

        integer(int64) :: d1, d2, d3
        complex(real32), dimension(:,:,:), pointer, contiguous :: tsrc, tdst
        integer(int64) :: sidx(3), didx(3)

        associate ( sdims => src%dims, ddims => dst%dims )
            call this%converter%secure_fortran_pointer(tsrc, src, astream)
            call this%converter%secure_fortran_pointer(tdst, dst, astream)

            do d3 = 1, sdims(3)
                sidx(3) = d3
                do d2 = 1, sdims(2)
                    sidx(2) = d2
                    do d1 = 1, sdims(1)
                        sidx(1) = d1

                        didx = sidx(perm)
                        tdst(didx(1),didx(2),didx(3)) = tsrc(sidx(1),sidx(2),sidx(3))
                    end do
                end do
            end do

            call this%converter%update_remote_and_release_pointer(tdst, dst, astream)
            call this%converter%release_pointer(tsrc, src, astream)
        end associate
    end subroutine permute_complex64_tensor_rank3

    subroutine permute_complex64_tensor_rank4(this, src, dst, perm, astream)
        class(tensor_permute_loops), intent(inout) :: this
        type(tensor_rank4), intent(in) :: src
        type(tensor_rank4), intent(inout) :: dst
        integer, dimension(4), intent(in) :: perm
        type(stream), intent(in), optional :: astream

        integer(int64) :: d1, d2, d3, d4
        complex(real32), dimension(:,:,:,:), pointer, contiguous :: tsrc, tdst
        integer(int64) :: sidx(4), didx(4)

        associate ( sdims => src%dims, ddims => dst%dims )
            call this%converter%secure_fortran_pointer(tsrc, src, astream)
            call this%converter%secure_fortran_pointer(tdst, dst, astream)

            do d4 = 1, sdims(4)
                sidx(4) = d4
                do d3 = 1, sdims(3)
                    sidx(3) = d3
                    do d2 = 1, sdims(2)
                        sidx(2) = d2
                        do d1 = 1, sdims(1)
                            sidx(1) = d1

                            didx = sidx(perm)
                            tdst(didx(1),didx(2),didx(3),didx(4)) = tsrc(sidx(1),sidx(2),sidx(3),sidx(4))
                        end do
                    end do
                end do
            end do

            call this%converter%update_remote_and_release_pointer(tdst, dst, astream)
            call this%converter%release_pointer(tsrc, src, astream)
        end associate
    end subroutine permute_complex64_tensor_rank4

    subroutine permute_complex64_tensor_rank5(this, src, dst, perm, astream)
        class(tensor_permute_loops), intent(inout) :: this
        type(tensor_rank5), intent(in) :: src
        type(tensor_rank5), intent(inout) :: dst
        integer, dimension(5), intent(in) :: perm
        type(stream), intent(in), optional :: astream

        integer(int64) :: d1, d2, d3, d4, d5
        complex(real32), dimension(:,:,:,:,:), pointer, contiguous :: tsrc, tdst
        integer(int64) :: sidx(5), didx(5)

        associate ( sdims => src%dims, ddims => dst%dims )
            call this%converter%secure_fortran_pointer(tsrc, src, astream)
            call this%converter%secure_fortran_pointer(tdst, dst, astream)

            do d5 = 1, sdims(5)
                sidx(5) = d5
                do d4 = 1, sdims(4)
                    sidx(4) = d4
                    do d3 = 1, sdims(3)
                        sidx(3) = d3
                        do d2 = 1, sdims(2)
                            sidx(2) = d2
                            do d1 = 1, sdims(1)
                                sidx(1) = d1

                                didx = sidx(perm)
                                tdst(didx(1),didx(2),didx(3),didx(4),didx(5)) = tsrc(sidx(1),sidx(2),sidx(3),sidx(4),sidx(5))
                            end do
                        end do
                    end do
                end do
            end do

            call this%converter%update_remote_and_release_pointer(tdst, dst, astream)
            call this%converter%release_pointer(tsrc, src, astream)
        end associate
    end subroutine permute_complex64_tensor_rank5

    subroutine permute_real64_matrix(this, src, dst, perm, astream)
        class(tensor_permute_loops), intent(inout) :: this
        type(matrix), intent(in) :: src
        type(matrix), intent(inout) :: dst
        integer, dimension(2), intent(in) :: perm
        type(stream), intent(in), optional :: astream

        integer(int64) :: d1, d2
        real(real64), dimension(:,:), pointer, contiguous :: tsrc, tdst
        integer(int64) :: sidx(2), didx(2)

        associate ( sdims => src%dims, ddims => dst%dims )
            call this%converter%secure_fortran_pointer(tsrc, src, astream)
            call this%converter%secure_fortran_pointer(tdst, dst, astream)

            do d2 = 1, sdims(2)
                sidx(2) = d2
                do d1 = 1, sdims(1)
                    sidx(1) = d1

                    didx = sidx(perm)
                    tdst(didx(1),didx(2)) = tsrc(sidx(1),sidx(2))
                end do
            end do

            call this%converter%update_remote_and_release_pointer(tdst, dst, astream)
            call this%converter%release_pointer(tsrc, src, astream)
        end associate
    end subroutine permute_real64_matrix

    subroutine permute_real64_tensor_rank3(this, src, dst, perm, astream)
        class(tensor_permute_loops), intent(inout) :: this
        type(tensor_rank3), intent(in) :: src
        type(tensor_rank3), intent(inout) :: dst
        integer, dimension(3), intent(in) :: perm
        type(stream), intent(in), optional :: astream

        integer(int64) :: d1, d2, d3
        real(real64), dimension(:,:,:), pointer, contiguous :: tsrc, tdst
        integer(int64) :: sidx(3), didx(3)

        associate ( sdims => src%dims, ddims => dst%dims )
            call this%converter%secure_fortran_pointer(tsrc, src, astream)
            call this%converter%secure_fortran_pointer(tdst, dst, astream)

            do d3 = 1, sdims(3)
                sidx(3) = d3
                do d2 = 1, sdims(2)
                    sidx(2) = d2
                    do d1 = 1, sdims(1)
                        sidx(1) = d1

                        didx = sidx(perm)
                        tdst(didx(1),didx(2),didx(3)) = tsrc(sidx(1),sidx(2),sidx(3))
                    end do
                end do
            end do

            call this%converter%update_remote_and_release_pointer(tdst, dst, astream)
            call this%converter%release_pointer(tsrc, src, astream)
        end associate
    end subroutine permute_real64_tensor_rank3

    subroutine permute_real64_tensor_rank4(this, src, dst, perm, astream)
        class(tensor_permute_loops), intent(inout) :: this
        type(tensor_rank4), intent(in) :: src
        type(tensor_rank4), intent(inout) :: dst
        integer, dimension(4), intent(in) :: perm
        type(stream), intent(in), optional :: astream

        integer(int64) :: d1, d2, d3, d4
        real(real64), dimension(:,:,:,:), pointer, contiguous :: tsrc, tdst
        integer(int64) :: sidx(4), didx(4)

        associate ( sdims => src%dims, ddims => dst%dims )
            call this%converter%secure_fortran_pointer(tsrc, src, astream)
            call this%converter%secure_fortran_pointer(tdst, dst, astream)

            do d4 = 1, sdims(4)
                sidx(4) = d4
                do d3 = 1, sdims(3)
                    sidx(3) = d3
                    do d2 = 1, sdims(2)
                        sidx(2) = d2
                        do d1 = 1, sdims(1)
                            sidx(1) = d1

                            didx = sidx(perm)
                            tdst(didx(1),didx(2),didx(3),didx(4)) = tsrc(sidx(1),sidx(2),sidx(3),sidx(4))
                        end do
                    end do
                end do
            end do

            call this%converter%update_remote_and_release_pointer(tdst, dst, astream)
            call this%converter%release_pointer(tsrc, src, astream)
        end associate
    end subroutine permute_real64_tensor_rank4

    subroutine permute_real64_tensor_rank5(this, src, dst, perm, astream)
        class(tensor_permute_loops), intent(inout) :: this
        type(tensor_rank5), intent(in) :: src
        type(tensor_rank5), intent(inout) :: dst
        integer, dimension(5), intent(in) :: perm
        type(stream), intent(in), optional :: astream

        integer(int64) :: d1, d2, d3, d4, d5
        real(real64), dimension(:,:,:,:,:), pointer, contiguous :: tsrc, tdst
        integer(int64) :: sidx(5), didx(5)

        associate ( sdims => src%dims, ddims => dst%dims )
            call this%converter%secure_fortran_pointer(tsrc, src, astream)
            call this%converter%secure_fortran_pointer(tdst, dst, astream)

            do d5 = 1, sdims(5)
                sidx(5) = d5
                do d4 = 1, sdims(4)
                    sidx(4) = d4
                    do d3 = 1, sdims(3)
                        sidx(3) = d3
                        do d2 = 1, sdims(2)
                            sidx(2) = d2
                            do d1 = 1, sdims(1)
                                sidx(1) = d1

                                didx = sidx(perm)
                                tdst(didx(1),didx(2),didx(3),didx(4),didx(5)) = tsrc(sidx(1),sidx(2),sidx(3),sidx(4),sidx(5))
                            end do
                        end do
                    end do
                end do
            end do

            call this%converter%update_remote_and_release_pointer(tdst, dst, astream)
            call this%converter%release_pointer(tsrc, src, astream)
        end associate
    end subroutine permute_real64_tensor_rank5

    subroutine permute_real32_matrix(this, src, dst, perm, astream)
        class(tensor_permute_loops), intent(inout) :: this
        type(matrix), intent(in) :: src
        type(matrix), intent(inout) :: dst
        integer, dimension(2), intent(in) :: perm
        type(stream), intent(in), optional :: astream

        integer(int64) :: d1, d2
        real(real32), dimension(:,:), pointer, contiguous :: tsrc, tdst
        integer(int64) :: sidx(2), didx(2)

        associate ( sdims => src%dims, ddims => dst%dims )
            call this%converter%secure_fortran_pointer(tsrc, src, astream)
            call this%converter%secure_fortran_pointer(tdst, dst, astream)

            do d2 = 1, sdims(2)
                sidx(2) = d2
                do d1 = 1, sdims(1)
                    sidx(1) = d1

                    didx = sidx(perm)
                    tdst(didx(1),didx(2)) = tsrc(sidx(1),sidx(2))
                end do
            end do

            call this%converter%update_remote_and_release_pointer(tdst, dst, astream)
            call this%converter%release_pointer(tsrc, src, astream)
        end associate
    end subroutine permute_real32_matrix

    subroutine permute_real32_tensor_rank3(this, src, dst, perm, astream)
        class(tensor_permute_loops), intent(inout) :: this
        type(tensor_rank3), intent(in) :: src
        type(tensor_rank3), intent(inout) :: dst
        integer, dimension(3), intent(in) :: perm
        type(stream), intent(in), optional :: astream

        integer(int64) :: d1, d2, d3
        real(real32), dimension(:,:,:), pointer, contiguous :: tsrc, tdst
        integer(int64) :: sidx(3), didx(3)

        associate ( sdims => src%dims, ddims => dst%dims )
            call this%converter%secure_fortran_pointer(tsrc, src, astream)
            call this%converter%secure_fortran_pointer(tdst, dst, astream)

            do d3 = 1, sdims(3)
                sidx(3) = d3
                do d2 = 1, sdims(2)
                    sidx(2) = d2
                    do d1 = 1, sdims(1)
                        sidx(1) = d1

                        didx = sidx(perm)
                        tdst(didx(1),didx(2),didx(3)) = tsrc(sidx(1),sidx(2),sidx(3))
                    end do
                end do
            end do

            call this%converter%update_remote_and_release_pointer(tdst, dst, astream)
            call this%converter%release_pointer(tsrc, src, astream)
        end associate
    end subroutine permute_real32_tensor_rank3

    subroutine permute_real32_tensor_rank4(this, src, dst, perm, astream)
        class(tensor_permute_loops), intent(inout) :: this
        type(tensor_rank4), intent(in) :: src
        type(tensor_rank4), intent(inout) :: dst
        integer, dimension(4), intent(in) :: perm
        type(stream), intent(in), optional :: astream

        integer(int64) :: d1, d2, d3, d4
        real(real32), dimension(:,:,:,:), pointer, contiguous :: tsrc, tdst
        integer(int64) :: sidx(4), didx(4)

        associate ( sdims => src%dims, ddims => dst%dims )
            call this%converter%secure_fortran_pointer(tsrc, src, astream)
            call this%converter%secure_fortran_pointer(tdst, dst, astream)

            do d4 = 1, sdims(4)
                sidx(4) = d4
                do d3 = 1, sdims(3)
                    sidx(3) = d3
                    do d2 = 1, sdims(2)
                        sidx(2) = d2
                        do d1 = 1, sdims(1)
                            sidx(1) = d1

                            didx = sidx(perm)
                            tdst(didx(1),didx(2),didx(3),didx(4)) = tsrc(sidx(1),sidx(2),sidx(3),sidx(4))
                        end do
                    end do
                end do
            end do

            call this%converter%update_remote_and_release_pointer(tdst, dst, astream)
            call this%converter%release_pointer(tsrc, src, astream)
        end associate
    end subroutine permute_real32_tensor_rank4

    subroutine permute_real32_tensor_rank5(this, src, dst, perm, astream)
        class(tensor_permute_loops), intent(inout) :: this
        type(tensor_rank5), intent(in) :: src
        type(tensor_rank5), intent(inout) :: dst
        integer, dimension(5), intent(in) :: perm
        type(stream), intent(in), optional :: astream

        integer(int64) :: d1, d2, d3, d4, d5
        real(real32), dimension(:,:,:,:,:), pointer, contiguous :: tsrc, tdst
        integer(int64) :: sidx(5), didx(5)

        associate ( sdims => src%dims, ddims => dst%dims )
            call this%converter%secure_fortran_pointer(tsrc, src, astream)
            call this%converter%secure_fortran_pointer(tdst, dst, astream)

            do d5 = 1, sdims(5)
                sidx(5) = d5
                do d4 = 1, sdims(4)
                    sidx(4) = d4
                    do d3 = 1, sdims(3)
                        sidx(3) = d3
                        do d2 = 1, sdims(2)
                            sidx(2) = d2
                            do d1 = 1, sdims(1)
                                sidx(1) = d1

                                didx = sidx(perm)
                                tdst(didx(1),didx(2),didx(3),didx(4),didx(5)) = tsrc(sidx(1),sidx(2),sidx(3),sidx(4),sidx(5))
                            end do
                        end do
                    end do
                end do
            end do

            call this%converter%update_remote_and_release_pointer(tdst, dst, astream)
            call this%converter%release_pointer(tsrc, src, astream)
        end associate
    end subroutine permute_real32_tensor_rank5

    subroutine synchronize(this, astream)
        class(tensor_permute_loops), intent(in) :: this
        type(stream), intent(in), optional :: astream

        continue
    end subroutine synchronize

    subroutine cleanup(this)
        class(tensor_permute_loops), intent(inout) :: this

        call this%converter%cleanup()
    end subroutine cleanup
end module tensor_permute_loops_module

max_stars_repo_name: cheshyre/ntcl-algorithms

max_issues_repo_name: cheshyre/ntcl-algorithms

max_forks_repo_name: cheshyre/ntcl-algorithms


text:
Steve Larsen was a former Davis resident (born 19700313 in Santa Clara County) who started cycling at 13 and went on to become a successful professional athlete in the areas of road bike, mountain bike, cyclocross and triathlons. Larsen began cycling in 1984 and jumped into cycling just as Greg LeMond was becoming a household name. He entered his first cyclocross race and bronzed in 1986 at the first Junior National Championship in Santa Cruz. Larsen went on to road race through the Junior National Team Program in Europe before moving up to the Amateur World Championships with the U.S. team.

Then firmly entrenched in road bike racing, Larsen got his first professional contract for the Motorola squad around 1991 and relocated to Italy where he raced professionally for almost four years and became a highly respected road racer, riding alongside teammate Lance Armstrong. By the close of 1994, he relocated to Bend, Oregon where he transitioned to mountain bike racing and won the NORBA National CrossCountry Mountainbike Champion title in 1998 and again in 2000. Despite a stellar career as Americas top male mountain bike racer, 2000 saw him barely miss a spot on the U.S. Olympic squad at the NORBA NCS event in Deer Valley, Utah.

Always looking for greater challenges, Steve began to eye triathlons. In May, 2000 he entered the California Wildflower Triathlon on a whim and surprised many by surviving the swim, moving into third after turning in the days fastest bike split and then holding on for fourth against a field that included some of the sports best. In July of the following year, Larsen entered the Isuzu Ironman USA Lake Placid triathlon and won the first Ironman he entered by scorching the bike course in 4:33 (a new course record) and setting a new overall course record by four minutes 8:33:11.

But these other triathlons were just stepping stones on Steves path to Hawaiis Ironman World Championship in October 2001. Back in Davis with his wife and two children now, after purchasing Wheelworks in April, Steve was again road racing with Pearl Izumi as his sponsor. The purchase of Wheelworks, renamed Steve Larsens Wheelworks, was especially meaningful because it was at this shop that he began riding. During this rampup period, Steve took full advantage of Davis many opportunties including its strong Davis Aquatic Masters swim program and working with former Motorola/7Eleven team doctor Massimo Testa and former 7Eleven pro Eric Heiden Dr. Eric Heiden. With his cycling skills and intensity, Larsen easily qualified for the event and was off to Kona in October. There he made an impressive showing, getting out of the water nearly 12 minutes behind the leaders then blazing through the 112mile bike leg, holding the front of the race until mile 10, ultimately finishing in 9th place but knowing exactly where I can improve for next year.

By 2003, Steve Larsen was done with professional racing and looking forward to spending more time focusing on his family than on sports. He retired in December and sold Wheelworks shortly thereafter, relocating once more to Bend, Oregon. There, he went into commercial real estate and shifted to competing on an amateur level.

Steve Larsen died 20090519 at the age of 39 in Bend, Oregon after collapsing during a running workout. A preliminary autopsy report found heart failure to be unlikely, and pointed to a possible viral causeFootnote(http://www.velonews.com/article/92268/autopsyshowsstevelarsendidnotdieofaheartattack).

 Steves Career Milestones
Steve is the only American to have ridden the Cycling World Championships in 4 different disciplines (road, mountain bike, track, and cyclocross).

1998 NORBA National CrossCountry Mountainbike Champion
2000 NORBA National CrossCountry Mountainbike Champion
2001 4th Place Wildflower Triathlon  1st American (Bike course record)
2001 1st Place HalfVineman Triathlon (Course Record)
2001 1st Place Ironman USA (Bike course record)
2002 6th Place Quelle Challenge (Ironman Distance)
2003 6th Place IM New Zealand (New bike course record 4:26)

Also see Steve Larsen Plaza.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
( useful extras )

: PLACE     ( a n a' -- )  2DUP C!  1+ SWAP MOVE ;

: ON   -1 SWAP ! ;
: OFF   0 SWAP ! ;

( comment to end of file )
: \S  BEGIN REFILL 0= UNTIL ;

( Multi-line comments )
: COMMENT  CHAR
    BEGIN  DUP DUP PARSE + C@ = NOT
    WHILE  REFILL NOT ABORT" comment?"
    REPEAT DROP ;

comment * this is a test *
comment *
this is also a test*
comment ~
this
is
typical
~
\ comment ~ missing delimiter fails...

( copied from standard )
: [ELSE] ( -- )
    1 BEGIN
        BEGIN PARSE-NAME DUP WHILE
          2DUP S" [IF]" COMPARE 0= IF
            2DROP 1+
          ELSE
            2DUP S" [ELSE]" COMPARE 0= IF
              2DROP 1- DUP IF 1+ THEN
            ELSE
              S" [THEN]" COMPARE 0= IF
                1-
              THEN
            THEN
          THEN ?DUP 0= IF EXIT THEN
        REPEAT 2DROP
    REFILL 0= UNTIL  DROP ; IMMEDIATE

: [IF]  0= IF  [COMPILE] [ELSE]  THEN ; IMMEDIATE
: [THEN] ; IMMEDIATE

: DEFINED ( -- f )  BL WORD FIND NIP ;

\ My needs are simpler than require, e.g.
\ need locals from opt/locals.f
\ need off : off 0 swap ! ;
: need  defined if [COMPILE] \ then ;
: from  include ; ( sugar )

\ Save dictionary image
: ?IOERR  ABORT" File I/O Error" ;
: SAVE ( <filename> -- ) \ format for include
    PARSE-NAME W/O CREATE-FILE ?IOERR
    HERE 0 DO
        DUP I C@ 0 <# ',' HOLD #S #> ROT
        I 15 AND 15 = IF WRITE-LINE ELSE WRITE-FILE THEN ?IOERR
    LOOP
    CLOSE-FILE ?IOERR ;
: SAVE-IMAGE ( <filename> -- )
    PARSE-NAME W/O CREATE-FILE ?IOERR
    DUP 0 HERE ROT WRITE-FILE ?IOERR
    CLOSE-FILE ?IOERR ;

max_stars_repo_name: amckewan/mcforth

max_issues_repo_name: amckewan/mcforth

max_forks_repo_name: amckewan/mcforth


text:
c
c
c     ###################################################
c     ##  COPYRIGHT (C)  1995  by  Jay William Ponder  ##
c     ##              All Rights Reserved              ##
c     ###################################################
c
c     ################################################################
c     ##                                                            ##
c     ##  subroutine diffeq  --  differential equation integration  ##
c     ##                                                            ##
c     ################################################################
c
c
c     "diffeq" performs the numerical integration of an ordinary
c     differential equation using an adaptive stepsize method to
c     solve the corresponding coupled first-order equations of the
c     general form dyi/dx = f(x,y1,...,yn) for yi = y1,...,yn
c
c     variables and parameters :
c
c     nvar      number of coupled first-order differential equations
c     y         contains the values of the dependent variables
c     x1        value of the beginning integration limit
c     x2        value of the ending integration limit
c     eps       relative accuracy required of the integration steps
c     h1        initial guess for the first integration stepsize
c     hmin      minimum allowed integration stepsize
c     nok       number of initially successful integration steps
c     nbad      number of integration steps that required retry
c
c     required external routines :
c
c     gvalue    subroutine to find the right-hand side of the
c                  first-order differential equations
c
c
      subroutine diffeq (nvar,y,x1,x2,eps,h1,hmin,nok,nbad,gvalue)
      include 'sizes.i'
      include 'iounit.i'
      real*8 tiny
      parameter (tiny=1.0d-30)
      integer i,nvar,nok,nbad
      integer nstep,maxstep
      real*8 x,x1,x2,h,h1
      real*8 eps,hnext
      real*8 hmin,hdid
      real*8 y(*)
      real*8, allocatable :: dydx(:)
      real*8, allocatable :: yscal(:)
      logical terminate
      character*7 status
      external gvalue
c
c
c     initialize starting limit, step size and status counters
c
      terminate = .false.
      x = x1
      h = sign(h1,x2-x1)
      nstep = 0
      nok = 0
      nbad = 0
      maxstep = 1000
c
c     perform dynamic allocation of some local arrays
c
      allocate (dydx(nvar))
      allocate (yscal(nvar))
c
c     perform a series of individual integration steps
c
      do while (.not. terminate)
         call gvalue (x,y,dydx)
         do i = 1, nvar
            yscal(i) = abs(y(i)) + abs(h*dydx(i)) + tiny
         end do
c
c     set the final step to stop at the integration limit
c
         if ((x+h-x2)*(x+h-x1) .gt. 0.0d0)  h = x2 - x
c
c     take a Bulirsch-Stoer integration step
c
         call bsstep (nvar,x,dydx,y,h,eps,yscal,hdid,hnext,gvalue)
c
c     mark the current step as either good or bad
c
         if (hdid .eq. h) then
            nok = nok + 1
            status = 'Success'
         else
            nbad = nbad + 1
            status = ' Retry '
         end if
c
c     update stepsize and get information about the current step
c
         h = hnext
         nstep = nstep + 1
         call gdastat (nstep,x,y,status)
c
c     test for convergence to the final integration limit
c
         if ((x-x2)*(x2-x1) .ge. 0.0d0) then
            write (iout,10)
   10       format (/,' DIFFEQ  --  Normal Termination',
     &                 ' at Integration Limit')
            terminate = .true.
         end if
c
c     test for a trial stepsize that is too small
c
         if (abs(hnext) .lt. hmin) then
            write (iout,20)
   20       format (/,' DIFFEQ  --  Incomplete Integration',
     &                 ' due to SmallStep')
            terminate = .true.
         end if
c
c     test for too many total integration steps
c
         if (nstep .ge. maxstep) then
            write (iout,30)
   30       format (/,' DIFFEQ  --  Incomplete Integration',
     &                 ' due to IterLimit')
            terminate = .true.
         end if
      end do
c
c     perform deallocation of some local arrays
c
      deallocate (dydx)
      deallocate (yscal)
      return
      end
c
c
c     ##############################################################
c     ##                                                          ##
c     ##  subroutine bsstep  --  Bulirsch-Stoer integration step  ##
c     ##                                                          ##
c     ##############################################################
c
c
c     "bsstep" takes a single Bulirsch-Stoer step with monitoring
c     of local truncation error to ensure accuracy
c
c     literature reference:
c
c     W. H. Press, S. A. Teukolsky, W. T. Vetterling and B. P.
c     Flannery, Numerical Recipes (Fortran), 2nd Ed., Cambridge
c     University Press, 1992, Section 16.4
c
c
      subroutine bsstep (nvar,x,dydx,y,htry,eps,yscal,hdid,hnext,gvalue)
      include 'sizes.i'
      include 'iounit.i'
      integer kmaxx,imax
      real*8 safe1,safe2
      real*8 redmax,redmin
      real*8 tiny,scalmx
      parameter (kmaxx=8)
      parameter (imax=kmaxx+1)
      parameter (safe1=0.25d0)
      parameter (safe2=0.7d0)
      parameter (redmax=1.0d-5)
      parameter (redmin=0.7d0)
      parameter (tiny=1.0d-30)
      parameter (scalmx=0.1d0)
      integer i,iq,k,kk,nvar
      integer km,kmax,kopt
      integer nseq(imax)
      real*8 eps,eps1,epsold
      real*8 h,hdid,hnext,htry
      real*8 errmax,fact,red
      real*8 scale,work,wrkmin
      real*8 x,xest,xnew
      real*8 dydx(*)
      real*8 y(*)
      real*8 yscal(*)
      real*8 a(imax)
      real*8 err(kmaxx)
      real*8 alf(kmaxx,kmaxx)
      real*8, allocatable :: yerr(:)
      real*8, allocatable :: ysav(:)
      real*8, allocatable :: yseq(:)
      logical first,reduct
      save a,alf,epsold,first
      save kmax,kopt,nseq,xnew
      external gvalue
      data first  / .true. /
      data epsold / -1.0d0 /
      data nseq   / 2,4,6,8,10,12,14,16,18 /
c
c
c     setup prior to the Bulirsch-Stoer integration step
c
      if (eps .ne. epsold) then
         hnext = -1.0d29
         xnew = -1.0d29
         eps1 = safe1 * eps
         a(1) = 1.0d0 + dble(nseq(1))
         do k = 1, kmaxx
            a(k+1) = a(k) + dble(nseq(k+1))
         end do
         do iq = 2, kmaxx
            do k = 1, iq-1
               alf(k,iq) = eps1**((a(k+1)-a(iq+1))/((a(iq+1)-a(1)+1.0d0)
     &                                                 *(2*k+1)))
            end do
         end do
         epsold = eps
         do kopt = 2, kmaxx-1
            kmax = kopt
            if (a(kopt+1) .gt. a(kopt)*alf(kopt-1,kopt))  goto 10
         end do
   10    continue
      end if
c
c     perform dynamic allocation of some local arrays
c
      allocate (yerr(nvar))
      allocate (ysav(nvar))
      allocate (yseq(nvar))
c
c     make an integration step using Bulirsch-Stoer method
c
      h = htry
      do i = 1, nvar
         ysav(i) = y(i)
      end do
      if (h.ne.hnext .or. x.ne.xnew) then
         first = .true.
         kopt = kmax
      end if
      reduct = .false.
   20 continue
      do k = 1, kmax
         xnew = x + h
         if (xnew .eq. x) then
            write (iout,30)
   30       format (' BSSTEP  --  Underflow of Step Size')
            call fatal
         end if
         call mmid (nseq(k),h,nvar,x,dydx,ysav,yseq,gvalue)
         xest = (h/dble(nseq(k)))**2
         call pzextr (k,nvar,xest,yseq,y,yerr)
         if (k .ne. 1) then
            errmax = tiny
            do i = 1, nvar
               errmax = max(errmax,abs(yerr(i)/yscal(i)))
            end do
            errmax = errmax / eps
            km = k - 1
            err(km) = (errmax/safe1)**(1.0d0/(2*km+1))
         end if
         if (k.ne.1 .and. (k.ge.kopt-1 .or. first)) then
            if (errmax .lt. 1.0d0)  goto 50
            if (k.eq.kmax .or. k.eq.kopt+1) then
               red = safe2 / err(km)
               goto 40
            else if (k .eq. kopt) then
               if (alf(kopt-1,kopt) .lt. err(km)) then
                  red = 1.0d0 / err(km)
                  goto 40
               end if
            else if (kopt .eq. kmax)then
               if (alf(km,kmax-1) .lt. err(km)) then
                  red = alf(km,kmax-1) * safe2 / err(km)
                  goto 40
               end if
            else if (alf(km,kopt) .lt. err(km)) then
               red = alf(km,kopt-1) / err(km)
               goto 40
            end if
         end if
      end do
   40 continue
      red = min(red,redmin)
      red = max(red,redmax)
      h = h * red
      reduct = .true.
      goto 20
   50 continue
      x = xnew
      hdid = h
      first = .false.
      wrkmin = 1.0d35
      do kk = 1, km
         fact = max(err(kk),scalmx)
         work = fact * a(kk+1)
         if (work .lt. wrkmin) then
            scale = fact
            wrkmin = work
            kopt = kk + 1
         end if
      end do
      hnext = h / scale
      if (kopt.ge.k .and. kopt.ne.kmax .and. .not.reduct) then
         fact = max(scale/alf(kopt-1,kopt),scalmx)
         if (a(kopt+1)*fact .le. wrkmin) then
            hnext = h / fact
            kopt = kopt + 1
         end if
      end if
c
c     perform deallocation of some local arrays
c
      deallocate (yerr)
      deallocate (ysav)
      deallocate (yseq)
      return
      end
c
c
c     ###########################################################
c     ##                                                       ##
c     ##  subroutine mmid  --  takes a modified midpoint step  ##
c     ##                                                       ##
c     ###########################################################
c
c
c     "mmid" implements a modified midpoint method to advance the
c     integration of a set of first order differential equations
c
c
      subroutine mmid (nstep,htot,nvar,xs,dydx,y,yout,gvalue)
      include 'sizes.i'
      integer i,k
      integer nstep,nvar
      real*8 htot,h,h2
      real*8 xs,x,temp
      real*8 dydx(*)
      real*8 y(*)
      real*8 yout(*)
      real*8, allocatable :: ym(:)
      real*8, allocatable :: yn(:)
      external gvalue
c
c
c     set substep size based on number of steps to be taken
c
      h = htot / dble(nstep)
      h2 = 2.0d0 * h
c
c     perform dynamic allocation of some local arrays
c
      allocate (ym(nvar))
      allocate (yn(nvar))
c
c     take the first substep and get values at ends of step
c
      do i = 1, nvar
         ym(i) = y(i)
         yn(i) = y(i) + h*dydx(i)
      end do
      x = xs + h
      call gvalue (x,yn,yout)
c
c     take the second and subsequent substeps
c
      do k = 2, nstep
         do i = 1, nvar
            temp = ym(i) + h2*yout(i)
            ym(i) = yn(i)
            yn(i) = temp
         end do
         x = x + h
         call gvalue (x,yn,yout)
      end do
c
c     complete the update of values for the last substep
c
      do i = 1, nvar
         yout(i) = 0.5d0 * (ym(i)+yn(i)+h*yout(i))
      end do
c
c     perform deallocation of some local arrays
c
      deallocate (ym)
      deallocate (yn)
      return
      end
c
c
c     ##############################################################
c     ##                                                          ##
c     ##  subroutine pzextr  --  polynomial extrapolation method  ##
c     ##                                                          ##
c     ##############################################################
c
c
c     "pzextr" is a polynomial extrapolation routine used during
c     Bulirsch-Stoer integration of ordinary differential equations
c
c
      subroutine pzextr (iest,nvar,xest,yest,yz,dy)
      include 'sizes.i'
      integer maxgda,imax
      parameter (maxgda=4*maxatm)
      parameter (imax=13)
      integer i,j,iest,nvar
      real*8 xest,delta
      real*8 f1,f2,q
      real*8 x(imax)
      real*8 yz(*)
      real*8 dy(*)
      real*8 yest(*)
      real*8, allocatable :: d(:)
      real*8, allocatable, save :: qcol(:,:)
      save x
c
c
c     perform dynamic allocation of some local arrays
c
      allocate (d(nvar))
      if (.not. allocated(qcol))  allocate (qcol(nvar,imax))
c
c     polynomial extrapolation needed for Bulirsch-Stoer step
c
      x(iest) = xest
      do j = 1, nvar
         dy(j) = yest(j)
         yz(j) = yest(j)
      end do
      if (iest .eq. 1) then
         do j = 1, nvar
            qcol(j,1) = yest(j)
         end do
      else
         do j = 1, nvar
            d(j) = yest(j)
         end do
         do i = 1, iest-1
            delta = 1.0d0 / (x(iest-i)-xest)
            f1 = xest * delta
            f2 = x(iest-i) * delta
            do j = 1, nvar
               q = qcol(j,i)
               qcol(j,i) = dy(j)
               delta = d(j) - q
               dy(j) = f1 * delta
               d(j) = f2 * delta
               yz(j) = yz(j) + dy(j)
            end do
         end do
         do j = 1, nvar
            qcol(j,iest) = dy(j)
         end do
      end if
c
c     perform deallocation of some local arrays
c
      deallocate (d)
      return
      end
c
c
c     ################################################################
c     ##                                                            ##
c     ##  subroutine gdastat  --  results for GDA integration step  ##
c     ##                                                            ##
c     ################################################################
c
c     "gdastat" finds the energy, radius of gyration, and average M2
c     for a GDA integration step; also saves the coordinates
c
c
      subroutine gdastat (nstep,beta,xx,status)
      implicit none
      include 'sizes.i'
      include 'atoms.i'
      include 'iounit.i'
      include 'math.i'
      include 'warp.i'
      integer i,nvar
      integer nstep
      real*8 beta
      real*8 e,energy
      real*8 rg,m2ave
      real*8 xx(*)
      character*7 status
c
c
c     translate optimization parameters to coordinates and M2's
c
      nvar = 0
      do i = 1, n
         nvar = nvar + 1
         x(i) = xx(nvar)
         nvar = nvar + 1
         y(i) = xx(nvar)
         nvar = nvar + 1
         z(i) = xx(nvar)
      end do
      do i = 1, n
         nvar = nvar + 1
         m2(i) = abs(xx(nvar))
      end do
c
c     get some info about the current integration step
c
      e = energy ()
      call gyrate (rg)
      m2ave = 0.0d0
      do i = 1, n
         m2ave = m2ave + m2(i)
      end do
      m2ave = m2ave / dble(n)
      write (iout,10)  nstep,log(beta)/logten,e,rg,
     &                 log(m2ave)/logten,status
   10 format (i6,2x,4f13.4,6x,a7)
c
c     save the current coordinates to a disk file
c
      call optsave (nstep,e,xx)
      return
      end

max_stars_repo_name: htna/HCsbLib

max_issues_repo_name: htna/HCsbLib

max_forks_repo_name: htna/HCsbLib


text:
c
c     cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine cmtrins(ma,mb,kr,kc,rst,cst,fac,a,lda,b,ldb)
c     ================================================================
c
      implicit   none
c
      integer    ma
      integer    mb
      integer    kr
      integer    kc
      integer    rst
      integer    cst
      integer    lda
      integer    ldb
      integer    i,j
c
      complex*16 fac
      complex*16 a(lda,ma)
      complex*16 b(ldb,mb)
c
c     ****************************************************************
c     a :     a matrix with dimension (lda x ma)
c     =
c
c     b :     a matrix assigned within dimension (kr x kc)
c     =
c     construct a new matrix b with dimension (ldb x mb):
c                            =
c
c             | b       |
c             | =       |
c     b   =   |         |
c     =       |   fac*a |
c             |       = |
c     ****************************************************************
c
      do j=1,kc
c        -------------------------------------------------------------
	 call zcopy(kr,a(1,j),1,b(rst+1,cst+j),1)
	 call zscal(kr,fac,b(rst+1,cst+j),1)
c        -------------------------------------------------------------
      enddo
c
      return
      end

max_stars_repo_name: rdietric/lsms

max_issues_repo_name: rdietric/lsms

max_forks_repo_name: hkershaw-brown/MuST


text:
      program prog
      character buf*24
      buf = 'abcdefghijklmnopqrstuvwx'
      write(6, '(3a)') '[', buf, ']'
      write(buf, *) 123
      write(6, '(3a)') '[', buf, ']'
      write(buf, *) 1000000, 2000000
      write(6, '(3a)') '[', buf, ']'
      end


text:
	SUBROUTINE DM_LSSF  ( iflno, sttype, ilstid, ilstnm, ilslat,
     +			      ilslon, ilselv, ilstat, ilcoun, ilstd2,
     +			      ilspri, ilswfo, ilwfo2, iret )
C************************************************************************
C* DM_LSSF								*
C*									*
C* This subroutine finds the location of the surface station keywords.	*
C* Both SLAT and SLON must be row or column keys.  The locations of 	*
C* the keywords STID, STNM, SELV, STAT, COUN, STD2, and SPRI are also	*
C* checked.  If present, they must be the same type as the SLAT and	*
C* SLON keys.  If a key is not found, the location is set to 0.		*
C*									*
C* DM_LSSF  ( IFLNO,  STTYPE, ILSTID, ILSTNM, ILSLAT, ILSLON, 		*
C*            ILSELV, ILSTAT, ILCOUN, ILSTD2, ILSPRI, ILSWFO,		*
C*	      ILWFO2, IRET )						*
C*									*
C* Input parameters:							*
C*	IFLNO		INTEGER		File number			*
C*									*
C* Output parameters:							*
C*	STTYPE		CHAR*		Location type:  ROW or COL	*
C*	ILSTID		INTEGER		Location of STID		*
C*	ILSTNM		INTEGER		Location of STNM		*
C*	ILSLAT		INTEGER		Location of SLAT		*
C*	ILSLON		INTEGER		Location of SLON		*
C*	ILSELV		INTEGER		Location of SELV		*
C*	ILSTAT		INTEGER		Location of STAT		*
C*	ILCOUN		INTEGER		Location of COUN		*
C*	ILSTD2		INTEGER		Location of extended stn ID	*
C*	ILSPRI		INTEGER		Location of priority flag	*
C*	ILSWFO		INTEGER		Location of SWFO		*
C*	ILWFO2		INTEGER		Location of extended SWFO	*
C*	IRET		INTEGER		Return code			*
C*					  0 = normal return		*
C*					 -4 = file not open		*
C*					-25 = invalid station keywords	*
C**									*
C* Log:									*
C* M. desJardins/GSFC	 5/87						*
C* K. Brill/NMC		 8/93	Created from DM_LSTN			*
C* T. Lee/GSC		10/97	Added SWFO and WFO2 keyword checks	*
C* T. Piper/GSC		11/98	Updated prolog				*
C* m. gamazaychikov/CWS 04/11   Add code for A2DB connectivity          *
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
	INCLUDE		'dmcmn.cmn'
C
	CHARACTER*(*)	sttype
	CHARACTER*4	tstid, tstnm, tslat, tslon, tselv, tstat, tcoun
	CHARACTER*4	tstd2, tspri, tswfo, twfo2
C-----------------------------------------------------------------------
C*	Initialize variables.
C
	iret   = 0
	sttype = ' '
	ilstid = 0
	ilstnm = 0
	ilslat = 0
	ilslon = 0
	ilselv = 0
	ilstat = 0
	ilcoun = 0
	ilstd2 = 0
	ilspri = 0
	ilswfo = 0
	ilwfo2 = 0
C
C*	Check that the file is open.
C
	CALL DM_CHKF  ( iflno, iret )
	IF  ( iret .ne. 0 )  RETURN
C
C*	Check for keywords.
C
	CALL DM_FKEY  ( iflno, 'STID', tstid, lstid, irstid )
	CALL DM_FKEY  ( iflno, 'STNM', tstnm, lstnm, irstnm )
	CALL DM_FKEY  ( iflno, 'SLAT', tslat, lslat, irslat )
	CALL DM_FKEY  ( iflno, 'SLON', tslon, lslon, irslon )
	CALL DM_FKEY  ( iflno, 'SELV', tselv, lselv, irselv )
	CALL DM_FKEY  ( iflno, 'STAT', tstat, lstat, irstat )
	CALL DM_FKEY  ( iflno, 'COUN', tcoun, lcoun, ircoun )
	CALL DM_FKEY  ( iflno, 'STD2', tstd2, lstd2, irstd2 )
	CALL DM_FKEY  ( iflno, 'SPRI', tspri, lspri, irspri )
	CALL DM_FKEY  ( iflno, 'SWFO', tswfo, lswfo, irswfo )
	CALL DM_FKEY  ( iflno, 'WFO2', twfo2, lwfo2, irwfo2 )
C
C*	Return error message if either latitude or longitude is missing.
C
	IF  ( ( irslat .ne. 0 )  .or.  ( irslon .ne. 0 )  .or.
     +	      ( tslat .ne. tslon ) )  THEN
	    iret = -25
	  ELSE IF  ( ( irstid .eq. 0 ) .and. ( tstid .ne. tslat ) ) THEN
	    iret = -25
	  ELSE IF  ( ( irstnm .eq. 0 ) .and. ( tstnm .ne. tslat ) ) THEN
	    iret = -25
	  ELSE IF  ( ( irselv .eq. 0 ) .and. ( tselv .ne. tslat ) ) THEN
	    iret = -25
	  ELSE IF  ( ( irstat .eq. 0 ) .and. ( tstat .ne. tslat ) ) THEN
	    iret = -25
	  ELSE IF  ( ( ircoun .eq. 0 ) .and. ( tcoun .ne. tslat ) ) THEN
	    iret = -25
	  ELSE IF  ( ( irstd2 .eq. 0 ) .and. ( tstd2 .ne. tslat ) ) THEN
	    iret = -25
	  ELSE IF  ( ( irspri .eq. 0 ) .and. ( tspri .ne. tslat ) ) THEN
	    iret = -25
	  ELSE IF  ( ( irswfo .eq. 0 ) .and. ( tswfo .ne. tslat ) ) THEN
	    iret = -25
	  ELSE IF  ( ( irwfo2 .eq. 0 ) .and. ( twfo2 .ne. tslat ) ) THEN
	    iret = -25
	  ELSE
	    ilstid = lstid
	    ilstnm = lstnm
	    ilslat = lslat
	    ilslon = lslon
	    ilselv = lselv
	    ilstat = lstat
	    ilcoun = lcoun
	    ilstd2 = lstd2
	    ilspri = lspri
	    ilswfo = lswfo
	    ilwfo2 = lwfo2
	    sttype = tslat
	END IF
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
program hello
  implicit none
  write(*,*) "hello"
  stop
end program hello

max_stars_repo_name: teigaCHENG/program

max_issues_repo_name: teigaCHENG/program

max_forks_repo_name: teigaCHENG/program


text:
C**********************************************************************
C***                                                                ***
C***  THIS PROGRAM CALCULATES                                       ***
C***  THE TOTAL CROSS-SECTION AND DISTRIBUTIONS                     ***
C***  FOR SQUARK-GLUINO PRODUCTION AT HADRON COLLIDERS              ***
C***  INCLUDING FULL SUSY-QCD CORRECTIONS                           ***
C***                                                                ***
C***  WRITTEN BY: W. BEENAKKER, R. HOPKER AND M. SPIRA              ***
C***                                                                ***
C**********************************************************************

      PROGRAM TOTALSG
      IMPLICIT REAL*8 (A-H,M-Z)
      IMPLICIT INTEGER (I,J)
      COMMON/IOUT/IPRINT
      COMMON/CONST1/S,ENERGY,ALPHAS,MS,MG,MT
      COMMON/CONST2/SCALE,SCAFAC,ICOLL,ISCAPT
      COMMON/CONST3/IPDFSET
      COMMON/CONST5/ILO,INLO,IONLYLO
      COMMON/CUT1/PTMIN,PTMAX
      COMMON/CUT2/YMIN,YMAX
      COMMON/FLAVOR/IFLAVOR,ITOTAL

C ---------------------------------------------------------------
C --- INPUT PARAMETERS, CAN BE CHANGED  -------------------------
C ---------------------------------------------------------------

C***  THE MASSES (IN GEV)

      MS = 280.D0
      MG = 200.D0
      MT = 175.D0

C***  THE COLLIDER TYPE ( P PBAR = 0, P P = 1 )

      ICOLL = 1

C***  THE CENTER OF MASS ENERGY (IN GEV)

      ENERGY = 14000D0

C***  THE SET OF PARTON DENSITIES (GRV = 0, MRSAP = 1, PDFLIB = 2)
C***  FOR PDFLIB PLEASE MAKE CHANGES IN SUBROUTINES INILO AND ININLO
C***  IN FILE INITPDF.F

      IPDFSET = 2

C***  THE INITIAL STATE
C***  ALL = 0, G G = 1, Q Q = 2, G Q = 3, MAJOR=4, MINOR=5

      IFLAVOR = 0

C***  THE CROSS-SECTION WITH CUTS (0)
C***  THE TOTAL CROSS-SECTION WITHOUT CUTS IN A FASTER WAY (1)

      ITOTAL = 1

C***  THE SCALE FOR RENORMALIZATION AND FACTORIZATION
C***  ISCAPT = 0  --> SCALE = (MS+MG)/2. * SCAFAC
C***  ISCAPT = 1  --> SCALE = SQRT(MG**2 +PT**2) * SCAFAC
C***                  ONLY FOR ITOTAL = 0
C***  DEFAULT FOR SCAFAC = 1.0

      ISCAPT = 0
      SCAFAC = 1.D0

C***  THE CUTS ON THE CROSS-SECTION IN PT (DEFAULT: 0,ENERGY )
C***  PT IS THE TRANSVERSE MOMENTUM OF THE GLUINO
C***  PT IS ONLY DEFINED FOR POSITIVE VALUES: PTMIN >= 0
C***  EQUAL LOWER AND UPPER CUT GIVES DSIGMA/DPT
C***  ITOTAL = 0 NECESSARY

      PTMIN = 0D0
      PTMAX = ENERGY

C***  THE CUTS ON THE CROSS-SECTION IN Y (DEFAULT: 0,+9.99)
C***  Y IS THE RAPIDITY OF THE GLUINO
C***  RAPIDITY IS ONLY DEFINED FOR POSITIVE VALUES: YMIN >= 0
C***  EQUAL LOWER AND UPPER CUT GIVES DSIGMA/DY
C***  ITOTAL = 0 NECESSARY

      YMIN =  0.D0
      YMAX = +9.99D0

C***  IONLYLO = 0 CALCULATES BORN AND NLO CROSS-SECTIONS
C***  IONLYLO = 1 CALCULATES ONLY THE BORN CROSS-SECTION

      IONLYLO = 0

C***  THE NUMBER OF VEGAS CALLS ( DEFAULT = 1000 )

      ILO   = 1000
      INLO  = 500

C***  PRINT VEGAS STATISTICS (10) OR NOT (0)

      IPRINT = 0

C --- PRINT THE HEADER ----------------------------------------

      CALL PRIHEADSG

C --- INITIALIZE VEGAS -----------------------------------------

      CALL RSTART(12,34,56,78)

C ---------------------------------------------------------------
C --- INTEGRATION BY VEGAS --------------------------------------
C --- CALCULATION OF CROSS-SECTIONS AND DISTRIBUTIONS -----------
C --- CAN BE CHANGED --------------------------------------------
C ---------------------------------------------------------------

      do img=160,1000,20
         do ims=160,1000,20

            mg = dble(img)
            ms = dble(ims)

C***  CALCULATE THE CROSS-SECTIONS IN LO AND NLO
            CALL INTEGSG(RESLO,ERRLO,RESNLO,ERRNLO)

C***  PRINT THE CROSS-SECTIONS IN LO AND NLO AND THEIR RELATIVE ERRORS
            CALL PRIRESSG(RESLO,ERRLO,RESNLO,ERRNLO)

         end do
      end do

c$$$C***  CALCULATE THE DIFFERENTIAL CROSS-SECTION DSIGMA/DPT
c$$$C***  WITH THE SCALE Q**2 = MG**2 + PT**2
c$$$C***  FOR PT = 50, 100, 150 GEV
c$$$      ITOTAL = 0
c$$$      ISCAPT = 1
c$$$      DO 100 I = 1,3
c$$$         PTMIN = 50.D0 * I
c$$$         PTMAX = PTMIN
c$$$         CALL INTEGSG(RESLO,ERRLO,RESNLO,ERRNLO)
c$$$         CALL PRIRESSG(RESLO,ERRLO,RESNLO,ERRNLO)
c$$$ 100  CONTINUE

      STOP
      END

C**********************************************************************

max_stars_repo_name: sliem/docker-hep

max_issues_repo_name: sliem/docker-hep

max_forks_repo_name: sliem/docker-hep


text:
!
! Format statement strings containing single apostrophes fails to parse.
!
program test_format_strings
050 format(' num''ber:')
100 FORMAT ('"')
200 FORMAT ("""")
300 FORMAT ('hellno')
 10200 format(' number of separate'' processes to treat:')
end program

max_stars_repo_name: OpenFortranProject/ofp-sdf

max_issues_repo_name: OpenFortranProject/ofp-sdf

max_forks_repo_name: OpenFortranProject/ofp-sdf


text:
!! -*- mode: F90; -*-
!!
!! Copyright (c) 2013, Georgia Tech Research Corporation
!! All rights reserved.
!!
!! Author(s): Neil T. Dantam <ntd@gatech.edu>
!! Georgia Tech Humanoid Robotics Lab
!! Under Direction of Prof. Mike Stilman <mstilman@cc.gatech.edu>
!!
!!
!! This file is provided under the following "BSD-style" License:
!!
!!
!!   Redistribution and use in source and binary forms, with or
!!   without modification, are permitted provided that the following
!!   conditions are met:
!!
!!   * Redistributions of source code must retain the above copyright
!!     notice, this list of conditions and the following disclaimer.
!!
!!   * Redistributions in binary form must reproduce the above
!!     copyright notice, this list of conditions and the following
!!     disclaimer in the documentation and/or other materials provided
!!     with the distribution.
!!
!!   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND
!!   CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,
!!   INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
!!   MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
!!   DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR
!!   CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
!!   SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
!!   LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
!!   USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED
!!   AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
!!   LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
!!   ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!!   POSSIBILITY OF SUCH DAMAGE.

!!! Kalman Filters !!!

subroutine rfx_lqg_kf_predict_work( n_x, n_u, A, x, B, u, P, V, work ) &
     bind( C, name="rfx_lqg_kf_predict_work" )
  integer(C_SIZE_T), intent(in), value :: n_x, n_u
  real(C_DOUBLE), intent(in) :: A(n_x,n_x), B(n_x, n_u), u(n_u), V(n_x,n_x)
  real(C_DOUBLE), intent(inout) :: x(n_x), P(n_x,n_x)
  real(C_DOUBLE), intent(out) :: work(n_x,n_x)

  call state(work(:,1))
  call covariance(work)

contains

  subroutine state( tmpx )
    real(C_DOUBLE), intent(out) :: tmpx(:)
    ! x = A*x + B*u
    tmpx = matmul(A,x)
    x = tmpx
    tmpx = matmul(B,u)
    x = x + tmpx
  end subroutine state

  subroutine covariance( tmpp )
    real(C_DOUBLE), intent(out) :: tmpp(:,:)
    ! P = A*P*A^T + V
    tmpp = matmul(A,P)
    P = matmul(tmpp,transpose(A))
    P = P + V
  end subroutine covariance

end subroutine rfx_lqg_kf_predict_work

subroutine rfx_lqg_kf_gain_work( n_x, n_z, C, P, W, K, work_xz, work_zz ) &
     bind( C, name="rfx_lqg_kf_gain_work" )
  integer(C_SIZE_T), intent(in), value :: n_x, n_z
  real(C_DOUBLE), intent(in) :: C(n_z,n_x), P(n_x, n_x), W(n_z,n_z)
  real(C_DOUBLE), intent(out) :: work_xz(n_x,n_z), work_zz(n_z,n_z), K(n_x, n_z)
  integer :: i
  !! K = P * C**T * (C * P * C**T + W)**-1
  work_xz = matmul(P,transpose(C))
  work_zz = matmul(C,work_xz)
  work_zz = work_zz + W
  i = aa_la_inv( n_z, work_zz )
  K = matmul(work_xz, work_zz)
end subroutine rfx_lqg_kf_gain_work

subroutine rfx_lqg_kf_innovate_work( n_x, n_z, x, K, z, C, workx, workz ) &
     bind( C, name="rfx_lqg_kf_innovate_work" )
  integer(C_SIZE_T) :: n_x, n_z
  real(C_DOUBLE), intent(inout) :: x(n_x)
  real(C_DOUBLE), intent(in) :: K(n_x, n_z), z(n_z), C(n_x,n_z)
  real(C_DOUBLE), intent(out) :: workz(n_z), workx(n_x)
  !! xh = x + K * (z - C*x)
  workz =  matmul(C,x)
  workz = z - workz
  workx = matmul(K,workz)
  x = x + workx
end subroutine rfx_lqg_kf_innovate_work

subroutine rfx_lqg_kf_covariance_work( n_x, n_z, C, P, K, work_xx2 ) &
     bind( C, name="rfx_lqg_kf_covariance_work" )
  integer(C_SIZE_T), intent(in), value :: n_x, n_z
  real(C_DOUBLE), intent(inout) :: P(n_x,n_x)
  real(C_DOUBLE), intent(in) :: C(n_z, n_x), K(n_x, n_z)
  real(C_DOUBLE), intent(out) :: work_xx2(n_x,n_x,2)
  integer(C_SIZE_T) :: i
  !! P = (I - K*C) * P
  work_xx2(:,:,1) = matmul(K, C)
  work_xx2(:,:,1) = -work_xx2(:,:,1)
  forall (i = 1:n_x)
     work_xx2(i,i,1) = 1 + work_xx2(i,i,1)
  end forall
  work_xx2(:,:,2) = P
  P = matmul( work_xx2(:,:,1), work_xx2(:,:,1) )
end subroutine rfx_lqg_kf_covariance_work

elemental subroutine rfx_lqg_kf_simple(  x, z, p, v, w, k )
  real(C_DOUBLE), intent(in) :: z,v,w
  real(C_DOUBLE), intent(inout) :: x,p
  real(C_DOUBLE), intent(out) :: k

  ! A = C = I

  ! predict
  p = p + v

  ! correct
  k = p / (p+w)
  x = x + k * (z - x)
  p = (1-k) * p
end subroutine rfx_lqg_kf_simple

subroutine rfx_lqg_kf_simple_work( n_x, x, z, P, V, W, K ) &
     bind( C, name="rfx_lqg_kf_simple_work" )
  integer(C_SIZE_T), intent(in), value :: n_x
  real(C_DOUBLE), intent(in) :: z(n_x), W(n_x), V(n_x)
  real(C_DOUBLE), intent(inout) :: x(n_x), P(n_x)
  real(C_DOUBLE), intent(out) :: K(n_x)
  call rfx_lqg_kf_simple( x, z, P, V, W, K )
end subroutine rfx_lqg_kf_simple_work

subroutine rfx_lqg_kf_predict_diag_work( n_x, A, x, P, V ) &
     bind( C, name="rfx_lqg_kf_predict_diag_work" )
  integer(C_SIZE_T), intent(in), value :: n_x
  real(C_DOUBLE), intent(in) :: A(n_x), V(n_x)
  real(C_DOUBLE), intent(inout) :: x(n_x), P(n_x)
  ! state
  x = A*x
  ! covariance
  P = A*P*A + V
end subroutine rfx_lqg_kf_predict_diag_work

subroutine rfx_lqg_kf_correct_diag_work( n_x, z, C, x, P, W, K ) &
     bind( C, name="rfx_lqg_kf_correct_diag_work" )
  integer(C_SIZE_T), intent(in), value :: n_x
  real(C_DOUBLE), intent(in) :: z(n_x), C(n_x), W(n_x)
  real(C_DOUBLE), intent(inout) :: x(n_x), P(n_x)
  real(C_DOUBLE), intent(out) :: K(n_x)

  ! kalman gain
  K = (P * C ) / (C*P*C + W)

  ! state
  x = x + K * (z - C*x)

  ! covariance
  P = (1-K*C) * P
end subroutine rfx_lqg_kf_correct_diag_work

max_stars_repo_name: golems/reflex

max_issues_repo_name: golems/reflex

max_forks_repo_name: golems/reflex


! PR middle-end/49308
! { dg-do compile }
! { dg-options "-O2 -funroll-loops -g" }

subroutine foo(n, b, d, e)
  type t
    integer :: f
  end type t
  type s
    type(t), pointer :: g
  end type s
  type u
    type(s), dimension(:), pointer :: h
  end type
  integer :: i, k, n
  type(u), pointer :: a, e
  character(len=250) :: b, c, d
  logical :: l
  do i = 1, n
    j = i - 1
    if (j/=0) c = trim(b) // adjustl(d(j))
  end do
  a => e
  do k = 1, size(a%h)
     l = (a%h(k)%g%f == a%h(1)%g%f)
     if (.not.(l)) call bar()
  enddo
end subroutine foo


text:
Living in Gloucester, UK and blogging about Gardening growing, cooking and eating food grown on my allotment
http://www.nipitinthebud.co.uk


text:
program funcao

    implicit none
    double precision :: area_circulo, volume_esfera, x, y

        write(*,*) "Digite o raio da circunferencia"
        read(*,*) x
        write(*,*) area_circulo(x)

	write(*,*) "Digite o raio da esfera"
	read(*,*) y
	write(*,*) volume_esfera(y)

end program

    double precision function area_circulo(r)
        double precision :: r
        area_circulo = 3.1415*r**2
        return
    end function area_circulo

    double precision function volume_esfera(r)
	double precision :: r
	volume_esfera = (4d0/3d0)*3.1415*r**3
	return
    end function volume_esfera

max_stars_repo_name: j-rheinheimer/Fisica-Computacional-I

max_issues_repo_name: j-rheinheimer/Fisica-Computacional-I

max_forks_repo_name: j-rheinheimer/Fisica-Computacional-I