! Copyright (c) 2012 Joseph A. Levin
!
! Permission is hereby granted, free of charge, to any person obtaining a copy of this
! software and associated documentation files (the "Software"), to deal in the Software
! without restriction, including without limitation the rights to use, copy, modify, merge,
! publish, distribute, sublicense, and/or sell copies of the Software, and to permit
! persons to whom the Software is furnished to do so, subject to the following conditions:
!
! The above copyright notice and this permission notice shall be included in all copies or
! substantial portions of the Software.
!
! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,
! INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
! PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
! LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT
! OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
! DEALINGS IN THE SOFTWARE.

! FSON MODULE
!
! File:   fson.f95
! Author: Joseph A. Levin
!
! Created on January 8, 2013
!

program example_basic

    ! Typical usage should only require an explicit use of the fson module.
    ! The other modules will be used privatley by fson as required.
!    use fson

    ! declare a pointer variable.  Always use a pointer with fson_value.
!    type(fson_value), pointer :: value

    ! parse the json file
!    value => fson_parse("test1.json")

    ! print the parsed data to the console
!    call fson_print(value)

    ! extract data from the parsed value

    ! clean up
!    call fson_destroy(value)

end program example_basic


! { dg-do compile }
! PR fortran/66979
program p
  implicit none
  integer::i
  flush (iostat=i) ! { dg-error "UNIT number missing" }
end program p


text:
   use TestModule, only : dot
   type(dot) :: dot1
   write(*,*) dot1%<caret>x
end

max_stars_repo_name: LChernigovskaya/fortran-plugin

max_issues_repo_name: LChernigovskaya/fortran-plugin

max_forks_repo_name: LChernigovskaya/fortran-plugin


text:
! Lagrange particle tracking on one particle
! !INTERFACE:
!TODO: Needs to implement Ross & Sharples (2004)
SUBROUTINE LAGRANGE(nlev,zlev, nuh,w,zi,zp)
USE PARAM_MOD, only: dtdays, d_per_s
! dtdays: timestep in days
! d_per_s: how many seconds in one day
!
IMPLICIT NONE
!
! !INPUT PARAMETERS:
INTEGER, intent(in)                :: nlev          !Number of vertical grids
REAL,    intent(in)                :: zlev(0:nlev)  !Depth of vertical grids at w points
REAL,    intent(in)                :: nuh(0:nlev)   !Vertical diffusivity at w points
REAL,    intent(in)                :: w(0:nlev)     !Vertical advection at w points

! INPUT/OUTPUT PARAMETERS:
! the ith grid where the particle currently resides
INTEGER, intent(inout)             :: zi

! Z position of particles
REAL,    intent(inout)             :: zp
!
! !LOCAL VARIABLES:
INTEGER            :: i
REAL               :: rnd, rnd_var_inv, dt, step1, step2
REAL,parameter     :: visc_back=0.e-6,rnd_var=0.333333333
REAL               :: dz(nlev),dzn(nlev),step,zp_old
REAL               :: visc,rat,dt_inv,zloc
REAL               :: wint  !interpolated w
!EOP
!-----------------------------------------------------------------------
!BOC

dt = dtdays * d_per_s

!Judge whether the input is correct
IF (zp < zlev(0) .or. zp > zlev(nlev) .or. zi < 1 .or. zi > nlev) THEN
   write(6,*) "zp = ", zp
   write(6,*) "zi = ", zi
   stop "The particle out of the domain!"
ENDIF

dt_inv      = 1./dt
rnd_var_inv = 1./rnd_var

DO i = 1,nlev

   ! Depth of each grid
   dz(i) = zlev(i)-zlev(i-1)

   ! gradient of Kv
   dzn(i)= (nuh(i)-nuh(i-1))/dz(i)

END DO

!  local viscosity calculation
! correction suggested by Visser [1997]

rat  = (zp-zlev(zi-1))/dz(zi)
wint = rat*w(zi) +(1.-rat)*w(zi-1)  ! interpolate w

! The new position after dt/2
zloc = zp + 0.5*(dzn(zi)+wint)*dt

step = zloc - zp   ! The distance traveled

! First judge whether the particle jumps out of the domain
if (zloc > zlev(nlev)) then  !Reflective boundary at surface
   zloc = 2.*zlev(nlev) - zloc
   zi   = nlev
elseif (zloc < zlev(0)) then  !Open boundary at bottom, the particle got lost
   zp   = zloc
   zi   = -1
   return
endif

step = zloc-zp     ! The distance traveled

! The following is to find the index for the grid where the particle currently resides
do i = 1,nlev
   if (zlev(i) .gt. zloc .and. zlev(i-1) .le. zloc) then
       zi = i
       EXIT
   endif
end do   !End of diffusivity correction

rat  = (zloc-zlev(i-1))/dz(i)
visc = rat*nuh(i)+(1.-rat)*nuh(i-1)  ! interpolate Kv
wint = rat*w(i) +(1.-rat)*w(i-1)     ! interpolate w

IF (visc.lt.visc_back) visc=visc_back   !Kv background

zp_old = zp

CALL RANDOM_NUMBER(rnd)
rnd  = 2.d0*rnd - 1.d0
step1= SQRT(2d0 * rnd_var_inv * dt * visc)*rnd
step2= dt*(wint+dzn(i))
step = step1 + step2

zp = zp + step

! First judge whether the particle jumps out of the domain
IF (zp > zlev(nlev)) THEN
      zp = 2.*zlev(nlev) - zp   !Reflective boundary at surface
ELSEIF (zp < zlev(0)) THEN      ! Particle lost at bottom
      zi = -1
      return
ENDIF

!Compute new zi
do i=1,nlev
   if (zlev(i) .gt. zp .and. zlev(i-1) .le. zp) then
       zi = i
       EXIT
   endif
end do

return
END SUBROUTINE LAGRANGE

max_stars_repo_name: BingzhangChen/HOT

max_issues_repo_name: BingzhangChen/HOT

max_forks_repo_name: BingzhangChen/HOT


text:
PROGRAM VECTOR_LOOP
  CALL claw
  CALL claw_transformed
END

SUBROUTINE claw ( )

 INTEGER :: j
 INTEGER :: vec1 ( 0 : 10 )
 INTEGER :: claw_induction_0

 DO j = 0 , 10 , 1
  vec1 ( j ) = j
 END DO
!$acc data present(vec1(:))
!$acc parallel
!$acc loop gang vector
 DO claw_induction_0 = 1 , size ( vec1 , 1 )
  vec1 ( claw_induction_0 ) = vec1 ( claw_induction_0 ) + 10
 END DO
!$acc end parallel
!$acc end data
 PRINT * , vec1
END SUBROUTINE claw

SUBROUTINE claw_transformed
  INTEGER :: i = 10
  INTEGER, DIMENSION(0:10) :: vec1
  INTEGER :: claw_induc1 ! induction variable is declared

  DO j = 0, i
    vec1(j) = j
  END DO

  ! Do stmt is inserted
  DO claw_induc1=0,i
    ! index range are replaced with arrayIndex
    vec1(claw_induc1) = vec1(claw_induc1) + 10
  END DO

  PRINT*,vec1
END SUBROUTINE claw_transformed

max_stars_repo_name: clementval/claw-compiler

max_issues_repo_name: skosukhin/claw-compiler

max_forks_repo_name: skosukhin/claw-compiler


text:
program test
use iso_fortran_env
use mod_mpifw
implicit none

integer(4) :: i, a(3) , b(3)

call init_mpi()

    a = [( this_proc*i , i=1,3)]

    write(*,'(*(g0,1x))') this_proc,")", a

    call sync_mpi()
    if(this_proc==0) write(*,*) "=================="
    call sync_mpi()

    write(*,'(*(g0,1x))') this_proc,")",reduce_mpi(a, 0, "sum")

    call sync_mpi()
    if(this_proc==0) write(*,*) "=================="
    call sync_mpi()

    write(*,'(*(g0,1x))') this_proc,")",allreduce_mpi(a, "sum")

call finish_mpi()
end program

max_stars_repo_name: msvbd/Functional-MPI-Wrapper

max_issues_repo_name: msvbd/Functional-MPI-Wrapper

max_forks_repo_name: msvbd/Functional-MPI-Wrapper


text:
*=======================================================================
* File:    tmatn.f
* Author:  Alex Stivala
* Created: July 2008
*
* FORTRAN-77 implementation of tableau matching numeric (ie Omega
* matrix rather than actual tableau) using the SPSOQP implementation
* of Prof. Yinyu Ye's interior point QP solver.
*
* See detailed documentation in tmatd.f header comments.
*
* $Id: tmatn.F 2969 2009-11-22 01:29:39Z astivala $
*=======================================================================

      subroutine tmatn(omega1,ldo1,n1,omega2,ldo2,n2,
     $     ltype,lorder,dmat1,dmat2,score,soln,info)

      implicit none
*
* tmatn - find matching score of two Omega matrices using relaxed QP
*
*   The tableaux are specified here as double precision Omega matrices.
*

*
*     .. Scalar Arguments ..
      integer ldo1,n1,ldo2,n2,info
      logical ltype,lorder
      double precision score
*     ..
*     .. Array Arguments ..
      double precision   omega1(ldo1,*),omega2(ldo2,*)
      double precision   dmat1(ldo1,*),dmat2(ldo2,*)
      double precision   soln(*)
*     ..
*

*
* Arguments
* =========
*
* omega1 (input) DOUBLE PRECISION array, dimension (n1,n1)
*        Omega matrix for one structure. Symmetric.
*
* ldo1   (input) INTEGER
*        Leading dimension of omega1 array. ldo1 >= n1
*
* n1     (input) INTEGER
*        Dimension of omega1 matrix.
*
* omega2 (input) DOUBLE PRECISION array, dimension (n2,n2)
*        Omega matrix for second structure. Symmetric.
*
* ldo2   (input) INTEGER
*        Leading dimension of omega2 array. ldo2 >= n2
*
* n2     (input) INTEGER
*        Dimension of omega2 matrix.
*
* ltype  (input) LOGICAL
*        if true, disalow matches between SSEs that are not the same
*        type. This uses the entry on the main diagonal as the indicator
*        of SSE type, and penalizes matches between entries with
*        different types.
*
* lorder (input) LOGICAL
*        if true, penalize matches between SSEs not maintaining sequence
*        order between the tableaux i.e. if i < k and j >= l for i,k
*        indices in omega1 and j,l indices in omega2.
*
* dmat1  (input) DOUBLE PRECISION array, dimension (n1,n1)
*        SSE distance matrix for one structure. symmetric
*
* dmat2  (input) DOUBLE PRECISION array, dimension (n2,n2)
*        SSE distance matrix for second structure. symmetric
*
* score  (output) DOUBLE PRECISION
*        score of matching omega1 with omega2 using relaxed QP.
*
* soln   (output) DOUBLE PRECISION vector, dimension (n1*n2+n1+n2)
*        solution x vector of matching omega1 with omega2 using relaxed QP
*
* info   (output) INTEGER
*        on exit, status of the computation
*        =  0 : successful exit
*        = -1 : maximum order of Q matrix exceeded
*        = -2 : maximum number of constraints exceeded
*        = -3 : SOLVQP failed
*
*=======================================================================
*

*     .. Parameters ..

*     largest tableau (omega matrix) dimension allowed
      integer maxdim1,maxdim2
      parameter(maxdim1 = 110, maxdim2 = 110)
      integer maxbdim, maxcdim
      parameter (maxcdim = maxdim1*maxdim2+maxdim1+maxdim2)
      parameter (maxbdim = maxdim1+maxdim2)
*     Threshold for penalizing difference in distances between SSEs
*     (Angstroms) - if difference in distances between SSEs in the two
*     structures exceeds this then penalize that match
      double precision mxssed
      parameter (mxssed = 4.0d0)
*     ..
*     .. Local Scalars ..
      integer i,j,k,l,qdim,arows
*     ..
*     .. Local Arrays ..
*     Q: symmetric objective matrix
*     A: constraint left-hand matrix
*     b: constraint rhs column vector
*     c: objective column vector
*     x: solution vector, with objective value as last element
*     y: dual solution vector (Langrangian multiplier)
      double precision Q(maxcdim, maxcdim),
     $     A(maxbdim, maxcdim),
     $     b(maxbdim),
     $     c(maxcdim),
     $     x(maxcdim + 1),
     $     y(maxbdim)
*     ..
*     .. Intrinsic Functions ..
*     ..
*     .. External Subroutines and Functions ..
      external dcopy
      external solvqp
      double precision tscorn
      external tscorn

*     ..
*     .. Executable Statements ..
*
      qdim = n1*n2 + n1 + n2
      arows = n1 + n2
      if (qdim .gt. maxcdim) then
         write(*,*) 'Max order of Q matrix exceeded'
         info = -1
         return
      endif
      if (arows .gt. maxbdim) then
         write(*,*) 'Max number of constraints exceeded'
         info = -2
         return
      endif

*     zero Q and A and c. c will stay zero.
      do 20 j = 1, qdim
         c(j) = 0.0d0
         do 10 i = 1, qdim
            Q(i, j) = 0.0d0
            if (i .le. arows) then
               A(i, j) = 0.0d0
            endif
 10      continue
 20   continue

*     set b to ones
      do 30 j = 1, arows
         b(j) = 1.0d0
 30   continue

*     setup the constraint matrix A
      do 50 j = 1, n2
         do 40 i = 1, n1
            A(i, (i-1)*n2+j) = 1.0d0
 40      continue
 50   continue
      do 70 j = 1, n1
         do 60 i = 1, n2
            A(n1+i, (j-1)*n2+i) = 1.0d0
 60      continue
 70   continue
      do 80 i = 1, n1+n2
         A(i, n1*n2+i) = 1.0d0
 80   continue

C       write(*,99),( (A(i,j), j=1,maxcdim), i=1,maxbdim)
C  99   format (111(f8.4))

*     setup the objective matrix Q
*     NB values in Q negated since we are
*     minimizing not maximizing (c is zero so irrelevant there)

      do 200 i = 1, n1
         do 190 j = 1, n2
            do 180 k = 1, n1
               do 170 l = 1, n2
*                 diagonal enries are SSE type not angle
                  if (i .ne. k .and. j .ne. l) then
                     if (ltype .and.
     $                   (omega1(i,i) .ne. omega2(j,j) .or.
     $                    omega1(k,k) .ne. omega2(l,l))) then
*                       penalize matches between SSEs of different type
                        Q((i-1)*n2 + j, (k-1)*n2 + l) = 0.0d0
                     elseif (abs(dmat1(i,k) - dmat2(j,l)) .gt. mxssed)
*                       penalize matches where difference between SSE
*                       distances exceeds the threshold
     $               then
                        Q((i-1)*n2 + j, (k-1)*n2 + l) = 0.0d0
                     else
                        if (lorder .and.
     $                       ( (i .lt. k .and. j .gt. l) .or.
     $                         (i .gt. k .and. j .lt. l) )) then
*                       penalize SSEs that are out of order
                           Q((i-1)*n2 + j, (k-1)*n2 + l) = 1.0d0
                        else
                           Q((i-1)*n2 + j, (k-1)*n2 + l) = -1.0d0 *
     $                          tscorn(omega1(i,k), omega2(j,l))
                        endif
                     endif
                  endif
 170           continue
 180        continue
 190     continue
 200  continue

C       write(*,205),( (Q(i,j), j=1,maxcdim), i=1,maxcdim)
C  205  format (111(f8.4))

      call solvqp(Q,maxcdim,A,maxbdim,arows,qdim,b,c,x,y,info)
      if (info .ne. 0) then
         write(*,210) info
 210     format('SOLVQP failed, info = ',i4)
         info = -3
         return
      endif

      score = x(qdim+1)
      call dcopy(qdim, x, 1, soln, 1)
      info = 0
      return
      end

ext: f

max_stars_repo_name: stivalaa/qptabsearch

max_issues_repo_name: stivalaa/qptabsearch

max_forks_repo_name: stivalaa/qptabsearch


module compute_flux_module

  implicit none

  private

  public :: &
       godunov_flux_2d, &
       mol2ndord_flux_2d, &
       mol4thord_flux_2d_limited, &
       mol4thord_flux_2d_nolimit, &
       getmaxmin_bounds

contains

  subroutine godunov_flux_2d(lo, hi, dt, dx, &
                             phi,ph_lo,ph_hi, &
                             umac,  u_lo,  u_hi, &
                             vmac,  v_lo,  v_hi, &
                             flxx, fx_lo, fx_hi, &
                             flxy, fy_lo, fy_hi, &
                             phix_1d, phiy_1d, phix, phiy, slope, glo, ghi, nu,  uselimit)

    use slope_module, only: slopex, slopey, slopex_nolim, slopey_nolim

    integer, intent(in) :: lo(2), hi(2), glo(2), ghi(2), uselimit
    double precision, intent(in) :: dt, dx(2), nu
    integer, intent(in) :: ph_lo(2), ph_hi(2)
    integer, intent(in) ::  u_lo(2),  u_hi(2)
    integer, intent(in) ::  v_lo(2),  v_hi(2)
    integer, intent(in) :: fx_lo(2), fx_hi(2)
    integer, intent(in) :: fy_lo(2), fy_hi(2)
    double precision, intent(in   ) :: phi (ph_lo(1):ph_hi(1),ph_lo(2):ph_hi(2))
    double precision, intent(in   ) :: umac( u_lo(1): u_hi(1), u_lo(2): u_hi(2))
    double precision, intent(in   ) :: vmac( v_lo(1): v_hi(1), v_lo(2): v_hi(2))
    double precision, intent(  out) :: flxx(fx_lo(1):fx_hi(1),fx_lo(2):fx_hi(2))
    double precision, intent(  out) :: flxy(fy_lo(1):fy_hi(1),fy_lo(2):fy_hi(2))
    double precision, dimension(glo(1):ghi(1),glo(2):ghi(2)) :: &
         phix_1d, phiy_1d, phix, phiy, slope

    integer :: i, j
    double precision :: hdtdx(2)

    hdtdx = 0.5*(dt/dx)

    if(uselimit .eq. 0) then
       call slopex_nolim(glo, ghi, &
            phi, ph_lo, ph_hi, &
            slope, glo, ghi)
    else

       call slopex(glo, ghi, &
            phi, ph_lo, ph_hi, &
            slope, glo, ghi)
    endif

    ! compute phi on x faces using umac to upwind; ignore transverse terms
    do    j = lo(2)-1, hi(2)+1
       do i = lo(1)  , hi(1)+1

          if (umac(i,j) .lt. 0.d0) then
             phix_1d(i,j) = phi(i  ,j) - (0.5d0 + hdtdx(1)*umac(i,j))*slope(i  ,j)
          else
             phix_1d(i,j) = phi(i-1,j) + (0.5d0 - hdtdx(1)*umac(i,j))*slope(i-1,j)
          end if

       end do
    end do

    if(uselimit .eq. 0) then
       call slopey_nolim(glo, ghi, &
            phi, ph_lo, ph_hi, &
            slope, glo, ghi)
    else
       call slopey(glo, ghi, &
            phi, ph_lo, ph_hi, &
            slope, glo, ghi)
    endif
    ! compute phi on y faces using umac to upwind; ignore transverse terms
    do    j = lo(2)  , hi(2)+1
       do i = lo(1)-1, hi(1)+1

          if (vmac(i,j) .lt. 0.d0) then
             phiy_1d(i,j) = phi(i,j  ) - (0.5d0 + hdtdx(2)*vmac(i,j))*slope(i,j  )
          else
             phiy_1d(i,j) = phi(i,j-1) + (0.5d0 - hdtdx(2)*vmac(i,j))*slope(i,j-1)
          end if

       end do
    end do

    ! update phi on x faces by adding in y-transverse terms
    do    j = lo(2), hi(2)
       do i = lo(1), hi(1)+1

          if (umac(i,j) .lt. 0.d0) then
             phix(i,j) = phix_1d(i,j) &
                  - hdtdx(2)*( 0.5d0*(vmac(i  ,j+1)+vmac(i  ,j)) * (phiy_1d(i  ,j+1)-phiy_1d(i  ,j)) )
          else
             phix(i,j) = phix_1d(i,j) &
                  - hdtdx(2)*( 0.5d0*(vmac(i-1,j+1)+vmac(i-1,j)) * (phiy_1d(i-1,j+1)-phiy_1d(i-1,j)) )
          end if

          ! compute final x-fluxes
          ! including diffusive fluxes
          flxx(i,j) = phix(i,j)*umac(i,j)  - nu*(phi(i,j) - phi(i-1,j))/dx(1)

       end do
    end do

    ! update phi on y faces by adding in x-transverse terms
    do    j = lo(2), hi(2)+1
       do i = lo(1), hi(1)

          if (vmac(i,j) .lt. 0.d0) then
             phiy(i,j) = phiy_1d(i,j) &
                  - hdtdx(1)*( 0.5d0*(umac(i+1,j  )+umac(i,j  )) * (phix_1d(i+1,j  )-phix_1d(i,j  )) )
          else
             phiy(i,j) = phiy_1d(i,j) &
                  - hdtdx(1)*( 0.5d0*(umac(i+1,j-1)+umac(i,j-1)) * (phix_1d(i+1,j-1)-phix_1d(i,j-1)) )
          end if

          ! compute final y-fluxes
          ! including diffusive fluxes
          flxy(i,j) = phiy(i,j)*vmac(i,j)  - nu*(phi(i,j) - phi(i,j-1))/dx(2)

       end do
    end do

  end subroutine godunov_flux_2d

  subroutine mol2ndord_flux_2d(lo, hi, dt, dx, &
       phi,ph_lo,ph_hi, &
       umac,  u_lo,  u_hi, &
       vmac,  v_lo,  v_hi, &
       flxx, fx_lo, fx_hi, &
       flxy, fy_lo, fy_hi, &
       phix_1d, phiy_1d, slope, glo, ghi, nu, uselimit)

    use slope_module, only: slopex, slopey, slopex_nolim, slopey_nolim

    integer, intent(in) :: lo(2), hi(2), glo(2), ghi(2), uselimit
    double precision, intent(in) :: dt, dx(2), nu
    integer, intent(in) :: ph_lo(2), ph_hi(2)
    integer, intent(in) ::  u_lo(2),  u_hi(2)
    integer, intent(in) ::  v_lo(2),  v_hi(2)
    integer, intent(in) :: fx_lo(2), fx_hi(2)
    integer, intent(in) :: fy_lo(2), fy_hi(2)
    double precision, intent(in   ) :: phi (ph_lo(1):ph_hi(1),ph_lo(2):ph_hi(2))
    double precision, intent(in   ) :: umac( u_lo(1): u_hi(1), u_lo(2): u_hi(2))
    double precision, intent(in   ) :: vmac( v_lo(1): v_hi(1), v_lo(2): v_hi(2))
    double precision, intent(  out) :: flxx(fx_lo(1):fx_hi(1),fx_lo(2):fx_hi(2))
    double precision, intent(  out) :: flxy(fy_lo(1):fy_hi(1),fy_lo(2):fy_hi(2))
    double precision, dimension(glo(1):ghi(1),glo(2):ghi(2)) :: &
         phix_1d, phiy_1d, slope
    double precision phi1, phi2, ftemp, veltemp, diffterm

    integer :: i, j

    if(uselimit .eq. 0) then
       call slopex_nolim(glo, ghi, &
            phi, ph_lo, ph_hi, &
            slope, glo, ghi)
    else
       call slopex(glo, ghi, &
            phi, ph_lo, ph_hi, &
            slope, glo, ghi)
    endif

    ! compute phi on x faces using umac to upwind; ignore transverse terms
    do    j = lo(2)-1, hi(2)+1
       do i = lo(1)  , hi(1)+1

          phi1 =  phi(i  ,j) - 0.5d0*slope(i  ,j)
          phi2 =  phi(i-1,j) + 0.5d0*slope(i-1,j)

          if (umac(i,j) .lt. 0.d0) then
             phix_1d(i,j) = phi(i  ,j) - 0.5d0*slope(i  ,j)
          else
             phix_1d(i,j) = phi(i-1,j) + 0.5d0*slope(i-1,j)
          end if

       end do
    end do

    if(uselimit .eq. 0) then
       call slopey_nolim(glo, ghi, &
            phi, ph_lo, ph_hi, &
            slope, glo, ghi)
    else
       call slopey(glo, ghi, &
            phi, ph_lo, ph_hi, &
            slope, glo, ghi)
    endif

    ! compute phi on y faces using umac to upwind; ignore transverse terms
    do    j = lo(2)  , hi(2)+1
       do i = lo(1)-1, hi(1)+1

          if (vmac(i,j) .lt. 0.d0) then
             phiy_1d(i,j) = phi(i,j  ) - (0.5d0)*slope(i,j  )
          else
             phiy_1d(i,j) = phi(i,j-1) + (0.5d0)*slope(i,j-1)
          end if

       end do
    end do

    ! compute final x-fluxes
    do    j = lo(2), hi(2)
       do i = lo(1), hi(1)+1
          ! including diffusive fluxes
          veltemp  = umac(i,j)
          diffterm = - nu*(phi(i,j) - phi(i-1,j))/dx(1)
          phi1  = phix_1d(i,j)
          ftemp  = phix_1d(i,j)*umac(i,j)  - nu*(phi(i,j) - phi(i-1,j))/dx(1)
          flxx(i,j) = phix_1d(i,j)*umac(i,j)  - nu*(phi(i,j) - phi(i-1,j))/dx(1)
       end do
    end do

    ! compute final y-fluxes
    do    j = lo(2), hi(2)+1
       do i = lo(1), hi(1)
          ! including diffusive fluxes
          ftemp  = phiy_1d(i,j)*vmac(i,j)  - nu*(phi(i,j) - phi(i,j-1))/dx(2)

          flxy(i,j) = phiy_1d(i,j)*vmac(i,j)  - nu*(phi(i,j) - phi(i,j-1))/dx(2)
       end do
    end do

  end subroutine mol2ndord_flux_2d

  !phi coming in is assumed to be cell-averaged
  !velocity is pointwise on faces---no limiting here
  subroutine mol4thord_flux_2d_nolimit(lo, hi, dt, dx, &
                                       phi ,ph_lo,ph_hi, &
                                       umac,  u_lo,  u_hi, &
                                       vmac,  v_lo,  v_hi, &
                                       flxx, fx_lo, fx_hi, &
                                       flxy, fy_lo, fy_hi, &
                                       fluxptx, phiptx, phiavex, &
                                       fluxpty, phipty, phiavey, &
                                       phiptcc, glo, ghi, nu, &
                                       deblocell, debhicell, &
                                       hisidedebfacelo, hisidedebfacehi, &
                                       losidedebfacelo, losidedebfacehi, printstuff)

    integer, intent(in) :: lo(2), hi(2), glo(2), ghi(2), printstuff
    integer, intent(in) :: deblocell(2), debhicell(2)
    integer, intent(in) :: hisidedebfacelo(2), hisidedebfacehi(2)
    integer, intent(in) :: losidedebfacelo(2), losidedebfacehi(2)
    double precision, intent(in) :: dt, dx(2), nu
    integer, intent(in) :: ph_lo(2), ph_hi(2)
    integer, intent(in) ::  u_lo(2),  u_hi(2)
    integer, intent(in) ::  v_lo(2),  v_hi(2)
    integer, intent(in) :: fx_lo(2), fx_hi(2)
    integer, intent(in) :: fy_lo(2), fy_hi(2)
    double precision, intent(in   ) :: phi (ph_lo(1):ph_hi(1),ph_lo(2):ph_hi(2))
    double precision, intent(in   ) :: umac( u_lo(1): u_hi(1), u_lo(2): u_hi(2))
    double precision, intent(in   ) :: vmac( v_lo(1): v_hi(1), v_lo(2): v_hi(2))
    double precision, intent(  out) :: flxx(fx_lo(1):fx_hi(1),fx_lo(2):fx_hi(2))
    double precision, intent(  out) :: flxy(fy_lo(1):fy_hi(1),fy_lo(2):fy_hi(2))
    double precision, dimension(glo(1):ghi(1),glo(2):ghi(2)) :: &
         fluxptx, fluxpty,  phiptx, phipty,  phiavex, phiavey, phiptcc

    double precision :: diffflux, phicctemp, debtemp
    integer :: i, j
!    integer ::  numphi
!    double precision :: phitot

    fluxptx = 1.0d30
    fluxpty = 1.0d30
    phiptx  = 1.0d30
    phipty  = 1.0d30
    phiavex = 1.0d30
    phiavey = 1.0d30
    phiptcc = 1.0d30
    !STEP 0
    ! 2.1 get cell-centered phi so we can compute a pointwise, fourth order gradient at faces
    ! needed  for diffusive fluxes

    do    j = lo(2)-3, hi(2)+3
       do i = lo(1)-3, hi(1)+3
          phicctemp  = phi(i,j) - (1.0d0/24.d0)* &
               (phi(i+1,j  )+phi(i-1,j  )&
               +phi(i  ,j+1)+phi(i  ,j-1)&
               -4.d0*phi(i,j))
          phiptcc(i,j)  = phicctemp
       end do
    end do

    ! STEP ONE--- HYPERBOLIC FLUXES
    ! compute face average phi on x faces via eqn 17 of mccorquodale, colella

    do    j = lo(2)-2, hi(2)+2
       do i = lo(1)  , hi(1)+1

          debtemp  = &
                (7.d0/12.d0)*(phi(i  ,j) + phi(i-1,j)) &
               -(1.d0/12.d0)*(phi(i-2,j) + phi(i+1,j))

          phiavex(i,j)  = debtemp
       end do
    end do

    !same for y faces
    do    j = lo(2)  , hi(2)+1
       do i = lo(1)-2, hi(1)+2

          debtemp  = &
                (7.d0/12.d0)*(phi(i,j  ) + phi(i,j-1)) &
               -(1.d0/12.d0)*(phi(i,j-2) + phi(i,j+1))
          phiavey(i,j)  = debtemp

       end do
    end do

    !now get point valued phi at faces so we can multiply by point valued velocity
    ! phipt = phiave - (h^2/24)*(lapl^2d(phi_ave))
    ! while I am at it, multiply in pointwise velocity so we get pointwise *HYPERBOLIC* flux
    !also  get pointwise diffusive fluxes using 4th order finite differences of pointwise phi

    do    j = lo(2)-1, hi(2)+1
       do i = lo(1)  , hi(1)+1
          phiptx(i,j)  =    phiavex(i,j  ) - (1.0d0/24.d0)* &
               (phiavex(i,j+1) + phiavex(i,j-1) - 2.d0*phiavex(i,j))

          debtemp  = phiptx(i,j)

          diffflux  = (-nu/dx(1))* &
               ((27.0d0/24.0d0)*(phiptcc(i  ,j) - phiptcc(i-1,j)) &
               +( 1.0d0/24.0d0)*(phiptcc(i-2,j) - phiptcc(i+1,j)))

          fluxptx(i,j) = umac(i,j)*phiptx(i,j) + diffflux

!uncomment to just do diffusion
!          fluxptx(i,j) =  diffflux

       end do
    end do

    !same for y faces
    do    j = lo(2)  , hi(2)+1
       do i = lo(1)-1, hi(1)+1

          phipty(i,j)  = phiavey(i,j) - (1.0d0/24.d0)* &
               (phiavey(i+1,j)+phiavey(i-1,j)-2.d0*phiavey(i,j))

          debtemp  = phipty(i,j)

          diffflux  = (-nu/dx(2))* &
               ((27.0d0/24.0d0)*(phiptcc(i,j  ) - phiptcc(i,j-1)) &
               +( 1.0d0/24.0d0)*(phiptcc(i,j-2) - phiptcc(i,j+1)))

          fluxpty(i,j) = vmac(i,j)*phipty(i,j) + diffflux

!uncomment to just do diffusion
!          fluxpty(i,j) =  diffflux

       end do
    end do

    ! now transform pointwise fluxes into face-averaged fluxes
    !fluxave = fluxpt + (1/24)(Lapl2d(fluxpt))

    do    j = lo(2), hi(2)
       do i = lo(1), hi(1)+1
          debtemp  = fluxptx(i,j) + (1.0d0/24.d0)* &
               (fluxptx(i,j+1) + fluxptx(i,j-1) - 2.d0*fluxptx(i,j))

          flxx(i,j)  = debtemp

!debug   just set the flux = phiavex
!          flxx(i,j)  = phiavex(i,j)
!end debug
       end do
    end do

    do    j = lo(2), hi(2)+1
       do i = lo(1), hi(1)
          debtemp  = fluxpty(i,j) + (1.0d0/24.d0)* &
               (fluxpty(i+1,j) + fluxpty(i-1,j) - 2.d0*fluxpty(i,j))
          flxy(i,j)  = debtemp

!debug   just set the flux = phiavex
!          flxy(i,j)  = phiavey(i,j)
!end debug
       end do
    end do

 end subroutine mol4thord_flux_2d_nolimit

  subroutine getmaxmin_bounds(lo, hi, glo, ghi, phi, maxval, minval, imax, imin)

    integer, intent(in) :: lo(2), hi(2), glo(2), ghi(2)
    double precision, intent(in   ) :: phi (glo(1):ghi(1),glo(2):ghi(2))
    double precision, intent(out   ) :: maxval, minval
    integer, intent(out) :: imax(2), imin(2)
    integer :: i, j

    maxval = -1.0d30
    minval =  1.0d30

    do    j = lo(2), hi(2)
       do i = lo(1), hi(1)
          if(phi(i,j) .gt. maxval) then
             maxval = phi(i,j)
             imax(1) = i
             imax(2) = j
          endif
          if(phi(i,j) .lt. minval) then
             minval = phi(i,j)
             imin(1) = i
             imin(2) = j
          endif
       enddo
    enddo

    return
  end subroutine getmaxmin_bounds

  !phi coming in is assumed to be cell-averaged
  !velocity is pointwise on faces---this is the more complicated one that does limiting.
  subroutine mol4thord_flux_2d_limited(lo, hi, dt, dx, &
                                       phi ,ph_lo,ph_hi, &
                                       umac,  u_lo,  u_hi, &
                                       vmac,  v_lo,  v_hi, &
                                       flxx, fx_lo, fx_hi, &
                                       flxy, fy_lo, fy_hi, &
                                       fluxptx, phiptx, phiavex, &
                                       fluxpty, phipty, phiavey, &
                                       phiptcc, glo, ghi, nu, &
                                       deblocell, debhicell, &
                                       hisidedebfacelo, hisidedebfacehi, &
                                       losidedebfacelo, losidedebfacehi, printstuff)

    use amrex_mempool_module, only : bl_allocate, bl_deallocate

    integer, intent(in) :: lo(2), hi(2), glo(2), ghi(2), printstuff
    integer, intent(in) :: deblocell(2), debhicell(2)
    integer, intent(in) :: hisidedebfacelo(2), hisidedebfacehi(2)
    integer, intent(in) :: losidedebfacelo(2), losidedebfacehi(2)
    double precision, intent(in) :: dt, dx(2), nu
    integer, intent(in) :: ph_lo(2), ph_hi(2)
    integer, intent(in) ::  u_lo(2),  u_hi(2)
    integer, intent(in) ::  v_lo(2),  v_hi(2)
    integer, intent(in) :: fx_lo(2), fx_hi(2)
    integer, intent(in) :: fy_lo(2), fy_hi(2)
    double precision, intent(in   ) :: phi (ph_lo(1):ph_hi(1),ph_lo(2):ph_hi(2))
    double precision, intent(in   ) :: umac( u_lo(1): u_hi(1), u_lo(2): u_hi(2))
    double precision, intent(in   ) :: vmac( v_lo(1): v_hi(1), v_lo(2): v_hi(2))
    double precision, intent(  out) :: flxx(fx_lo(1):fx_hi(1),fx_lo(2):fx_hi(2))
    double precision, intent(  out) :: flxy(fy_lo(1):fy_hi(1),fy_lo(2):fy_hi(2))
    double precision, dimension(glo(1):ghi(1),glo(2):ghi(2)) :: &
         fluxptx, fluxpty,  phiptx, phipty,  phiavex, phiavey, phiptcc

    double precision, dimension(:,:), pointer, contiguous :: &
         dwfminux, dwfminuy, dwfplusx, dwfplusy, d2wfx, d2wfy, d2wcx, d2wcy, d3wx, d3wy, &
         wleftx, wrighx, wlefty, wrighy

    double precision :: diffflux, phicctemp, debtemp, mono1, mono2, signd2, d2wlim, &
         d2wc1, d2wc2, d2wc3, d2wf1, rho, d3wmin, d3wmax, dwthresh, maxw
    logical :: applylim
    integer :: i, j
    integer :: imax(2), imin(2), blo(2), bhi(2);
    double precision ::     maxval, minval
!    double precision :: phitot
!    integer ::  numphi

    call bl_allocate(dwfplusx,glo(1), ghi(1), glo(2), ghi(2))
    call bl_allocate(dwfminux,glo(1), ghi(1), glo(2), ghi(2))
    call bl_allocate(dwfplusy,glo(1), ghi(1), glo(2), ghi(2))
    call bl_allocate(dwfminuy,glo(1), ghi(1), glo(2), ghi(2))
    call bl_allocate(d2wfx   ,glo(1), ghi(1), glo(2), ghi(2))
    call bl_allocate(d2wcx   ,glo(1), ghi(1), glo(2), ghi(2))
    call bl_allocate(d3wx    ,glo(1), ghi(1), glo(2), ghi(2))
    call bl_allocate(d3wy    ,glo(1), ghi(1), glo(2), ghi(2))
    call bl_allocate(d2wfy   ,glo(1), ghi(1), glo(2), ghi(2))
    call bl_allocate(d2wcy   ,glo(1), ghi(1), glo(2), ghi(2))

    call bl_allocate(wleftx  ,glo(1), ghi(1), glo(2), ghi(2))
    call bl_allocate(wlefty  ,glo(1), ghi(1), glo(2), ghi(2))
    call bl_allocate(wrighx  ,glo(1), ghi(1), glo(2), ghi(2))
    call bl_allocate(wrighy  ,glo(1), ghi(1), glo(2), ghi(2))

    phiavey  = 1.0d30
    phiptcc  = 1.0d30
    dwfplusx = 1.0d30
    dwfminux = 1.0d30
    dwfplusy = 1.0d30
    dwfminuy = 1.0d30
    d2wfx    = 1.0d30
    d2wcx    = 1.0d30
    d3wx     = 1.0d30
    d2wfy    = 1.0d30
    d2wcy    = 1.0d30
    d3wy     = 1.0d30

    !STEP 0
    ! 2.1 get cell-centered phi so we can compute a pointwise, fourth order gradient at faces
    ! needed  for diffusive fluxes (eqn 16 of m + c)

    do    j = lo(2)-3, hi(2)+3
       do i = lo(1)-3, hi(1)+3
          phicctemp  = phi(i,j) - (1.0d0/24.d0)* &
               (phi(i+1,j  )+phi(i-1,j  )&
               +phi(i  ,j+1)+phi(i  ,j-1)&
               -4.d0*phi(i,j))
          phiptcc(i,j)  = phicctemp

       end do
    end do

    blo(1) = lo(1)-3
    bhi(1) = hi(1)+3
    blo(2) = lo(2)-3
    bhi(2) = hi(2)+3
    call getmaxmin_bounds(blo, bhi, glo, ghi, phiptcc, maxval, minval, imax, imin)

    ! STEP ONE--- HYPERBOLIC FLUXES
    ! compute face average phi on x faces via eqn 17 of mccorquodale, colella

    do    j = lo(2)-2, hi(2)+2
       do i = lo(1)-1, hi(1)+2

          debtemp  = &
                (7.d0/12.d0)*(phi(i  ,j) + phi(i-1,j)) &
               -(1.d0/12.d0)*(phi(i-2,j) + phi(i+1,j))

          phiavex(i,j)  = debtemp

       end do
    end do

    blo(2) = lo(2)-2
    bhi(2) = hi(2)+2
    blo(1) = lo(1)-1
    bhi(1) = hi(1)+2
    call getmaxmin_bounds(blo, bhi, glo, ghi, phiavex, maxval, minval, imax, imin)

    !same for y faces
    do    j = lo(2)-1, hi(2)+2
       do i = lo(1)-2, hi(1)+2

          debtemp  = &
                (7.d0/12.d0)*(phi(i,j  ) + phi(i,j-1)) &
               -(1.d0/12.d0)*(phi(i,j-2) + phi(i,j+1))
          phiavey(i,j)  = debtemp

       end do
    end do

    blo(2) = lo(2)-1
    bhi(2) = hi(2)+2
    blo(1) = lo(1)-2
    bhi(1) = hi(1)+2
    call getmaxmin_bounds(blo, bhi, glo, ghi, phiavey, maxval, minval, imax, imin)

    !this stuff is for limiting
    !The variable names are from mccorquodale & colella.   I could not come
    ! up with anything more sensible.   This is from section 2.4.1.
    do    j = lo(2)-2, hi(2)+2
       do i = lo(1)  , hi(1)
          dwfminux(i,j) =     phi(i  ,j) - phiavex(i,j)
          dwfplusx(i,j) = phiavex(i+1,j) -     phi(i,j)
       enddo
    enddo

    blo(2) = lo(2)-2
    bhi(2) = hi(2)+2
    blo(1) = lo(1)
    bhi(1) = hi(1)
    call getmaxmin_bounds(blo, bhi, glo, ghi, dwfminux, maxval, minval, imax, imin)
    call getmaxmin_bounds(blo, bhi, glo, ghi, dwfplusx, maxval, minval, imax, imin)

    do    j = lo(2)    , hi(2)
       do i = lo(1)-2  , hi(1)+2
          dwfminuy(i,j) =     phi(i,j  ) - phiavey(i,j)
          dwfplusy(i,j) = phiavey(i,j+1) -     phi(i,j)

       enddo
    enddo

    blo(2) = lo(2)
    bhi(2) = hi(2)
    blo(1) = lo(1)-2
    bhi(1) = hi(1)+2
    call getmaxmin_bounds(blo, bhi, glo, ghi, dwfminuy, maxval, minval, imax, imin)
    call getmaxmin_bounds(blo, bhi, glo, ghi, dwfplusy, maxval, minval, imax, imin)

    do    j = lo(2)-2, hi(2)+2
       do i = lo(1)-1, hi(1)+1
          d2wfx(i,j) = 6.0d0*(phiavex(i  ,j) + phiavex(i+1,j) - 2.0d0*phi(i,j))
          d2wcx(i,j) =       (    phi(i-1,j) +     phi(i+1,j) - 2.0d0*phi(i,j))

       enddo
    enddo

    blo(2) = lo(2)-2
    bhi(2) = hi(2)+2
    blo(1) = lo(1)-1
    bhi(1) = hi(1)+1
    call getmaxmin_bounds(blo, bhi, glo, ghi, d2wfx, maxval, minval, imax, imin)
    call getmaxmin_bounds(blo, bhi, glo, ghi, d2wcx, maxval, minval, imax, imin)

    do    j = lo(2)-1, hi(2)+1
       do i = lo(1)-2, hi(1)+2
          d2wfy(i,j) = 6.0d0*(phiavey(i,j  ) + phiavey(i,j+1) - 2.0d0*phi(i,j))
          d2wcy(i,j) =       (    phi(i,j-1) +     phi(i,j+1) - 2.0d0*phi(i,j))

       enddo
    enddo

    blo(2) = lo(2)-1
    bhi(2) = hi(2)+1
    blo(1) = lo(1)-2
    bhi(1) = hi(1)+2
    call getmaxmin_bounds(blo, bhi, glo, ghi, d2wfy, maxval, minval, imax, imin)
    call getmaxmin_bounds(blo, bhi, glo, ghi, d2wcy, maxval, minval, imax, imin)

    do    j = lo(2)-2, hi(2)+2
       do i = lo(1),   hi(1)+1
          d3wx(i,j) = d2wcx(i,j) - d2wcx(i-1,j)
       enddo
    enddo

    blo(2) = lo(2)-2
    bhi(2) = hi(2)+2
    blo(1) = lo(1)
    bhi(1) = hi(1)+1
    call getmaxmin_bounds(blo, bhi, glo, ghi, d3wx, maxval, minval, imax, imin)

    do    j = lo(2)  , hi(2)+1
       do i = lo(1)-2, hi(1)+2
          d3wy(i,j) = d2wcy(i,j) - d2wcy(i,j-1)
       enddo
    enddo

    blo(2) = lo(2)
    bhi(2) = hi(2)+1
    blo(1) = lo(1)-2
    bhi(1) = hi(1)+2
    call getmaxmin_bounds(blo, bhi, glo, ghi, d3wy, maxval, minval, imax, imin)

    !initialize left and right states
    do    j = lo(2)-2, hi(2)+2
       do i = lo(1)  , hi(1)+1
          wleftx(i,j) = phiavex(i,j)
          wrighx(i,j) = phiavex(i,j)

       enddo
    enddo
    do    j = lo(2)  , hi(2)+1
       do i = lo(1)-2, hi(1)+2
          wlefty(i,j) = phiavey(i,j)
          wrighy(i,j) = phiavey(i,j)
       enddo
    enddo

    !now for the really complicted stuff in 2.4.1
    !If there are extrema detected a cell, modify the states
    !that are extrapolate from that cell.  If monotonoic but there
    ! is a huge jump in the state, mollify that
    do    j = lo(2)-2, hi(2)+2
       do i = lo(1)  , hi(1)+1
          mono1 = dwfminux(i,j) * dwfplusx(i,j)
          mono2 = (phi(i,j)-phi(i-2,j))*(phi(i+2,j)-phi(i,j))
          if((mono1 .lt. 0.0d0).or.(mono2 .lt. 0.0d0)) then

             !see if all second derivs in sight are the same sign.   If they are,
             ! calculate d2wlim.     Otherwise, d2wlim = 0
             !this stuff comes from equation 26
             d2wlim = 0.0d0
             signd2 = 0.0d0
             d2wc1 = d2wcx(i-1,j)
             d2wc2 = d2wcx(i  ,j)
             d2wc3 = d2wcx(i+1,j)
             d2wf1 = d2wfx(i  ,j)
             d2wlim = min(abs(d2wf1), 1.25d0*abs(d2wc1))
             d2wlim = min(d2wlim    , 1.25d0*abs(d2wc2))
             d2wlim = min(d2wlim    , 1.25d0*abs(d2wc3))
             if((d2wc1.gt.0.0d0) .and. &
                (d2wc2.gt.0.0d0) .and. &
                (d2wc3.gt.0.0d0) .and. &
                (d2wf1.gt.0.0d0)) then
                signd2 = 1.0d0
                d2wlim = d2wlim*signd2
             else if((d2wc1.lt.0.0d0) .and. &
                     (d2wc2.lt.0.0d0) .and. &
                     (d2wc3.lt.0.0d0) .and. &
                     (d2wf1.lt.0.0d0)) then
                signd2 = -1.0d0
                d2wlim = d2wlim*signd2
             else
                d2wlim = 0.0d0
             endif
             !equation 27 below
             maxw = max(abs(phi(i,j)),abs(phi(i+1,j)))
             maxw = max(maxw         ,abs(phi(i-1,j)))
             maxw = max(maxw         ,abs(phi(i+2,j)))
             maxw = max(maxw         ,abs(phi(i-2,j)))
             dwthresh = 1.0d-12*(maxw)
             if(abs(d2wf1).lt.dwthresh) then
                rho = 0.0d0
             else
                rho = d2wlim/d2wf1
             endif
             if(rho .lt. (1.0d0 - 1.0d-12)) then
                applylim = .false.
             else
                !this stuff is equation 28
                d3wmin = min(d3wx(i,j),d3wx(i-1,j))
                d3wmin = min(d3wmin   ,d3wx(i+1,j))
                d3wmin = min(d3wmin   ,d3wx(i+2,j))
                d3wmax = max(d3wx(i,j),d3wx(i-1,j))
                d3wmax = max(d3wmax   ,d3wx(i+1,j))
                d3wmax = max(d3wmax   ,d3wx(i+2,j))
                dwthresh = 0.1d0*max(abs(d3wmax),abs(d3wmin))
                if((d3wmax - d3wmin).lt. dwthresh) then
                   applylim = .false.
                else
                   applylim = .true.
                endif
             endif
             !the following covered equations 29-32
             if(applylim) then
                if((dwfminux(i,j)*dwfplusx(i,j)).lt.0.0) then

                   wrighx(i,j) = phi(i,j) -  rho*dwfminux(i,j) !eqn 29 (corrected)
                   wleftx(i,j) = phi(i,j) +  rho*dwfplusx(i,j) !eqn 30

                else if(abs(dwfminux(i,j)) .ge. 2.0d0*abs(dwfplusx(i,j))) then
                   !equation 31
                   wrighx(i,j) = phi(i,j) -  2.0d0*(1.0d0-rho)*dwfplusx(i,j) - rho*dwfminux(i,j)
                else if(abs(dwfplusx(i,j)) .ge. 2.0d0*abs(dwfminux(i,j))) then
                   !equation 32
                   wleftx(i,j) = phi(i,j) +  2.0d0*(1.0d0-rho)*dwfminux(i,j) + rho*dwfplusx(i,j)
                endif
             endif
          else !monotonic but still need to check for big jumps in the state
             if(abs(dwfminux(i,j)) .ge. (2.0d0*abs(dwfplusx(i,j)))) then
                wrighx(i,j) = phi(i,j) -  2.0d0*dwfplusx(i,j) !eqn 33
             endif
             if(abs(dwfplusx(i,j)) .ge. (2.0d0*abs(dwfminux(i,j)))) then
                wleftx(i,j) = phi(i,j) +  2.0d0*dwfminux(i,j) !eqn 34
             endif
          endif
       enddo
    enddo
    !now we do it all again for the y direction.
    do    j = lo(2)  , hi(2)+1
       do i = lo(1)-2, hi(1)+2
          mono1 = dwfminuy(i,j) * dwfplusy(i,j)
          mono2 = (phi(i,j)-phi(i,j-2))*(phi(i,j+2)-phi(i,j))
          if((mono1 .lt. 0.0d0).or.(mono2 .lt. 0.0d0)) then

             !see if all second derivs in sight are the same sign.   If they are,
             ! calculate d2wlim.     Otherwise, d2wlim = 0
             !this stuff comes from equation 26
             d2wlim = 0.0d0
             signd2 = 0.0d0
             d2wc1 = d2wcy(i,j-1)
             d2wc2 = d2wcy(i,j  )
             d2wc3 = d2wcy(i,j+1)
             d2wf1 = d2wfy(i,j  )
             d2wlim = min(abs(d2wf1), 1.25d0*abs(d2wc1))
             d2wlim = min(d2wlim    , 1.25d0*abs(d2wc2))
             d2wlim = min(d2wlim    , 1.25d0*abs(d2wc3))
             if((d2wc1.gt.0.0d0) .and. &
                (d2wc2.gt.0.0d0) .and. &
                (d2wc3.gt.0.0d0) .and. &
                (d2wf1.gt.0.0d0)) then
                signd2 = 1.0d0
                d2wlim = d2wlim*signd2
             else if((d2wc1.lt.0.0d0) .and. &
                     (d2wc2.lt.0.0d0) .and. &
                     (d2wc3.lt.0.0d0) .and. &
                     (d2wf1.lt.0.0d0)) then
                signd2 = -1.0d0
                d2wlim = d2wlim*signd2
             else
                d2wlim = 0.0d0
             endif
             !equation 27 below
             maxw = max(abs(phi(i,j)),abs(phi(i,j+1)))
             maxw = max(maxw         ,abs(phi(i,j-1)))
             maxw = max(maxw         ,abs(phi(i,j+2)))
             maxw = max(maxw         ,abs(phi(i,j-2)))
             dwthresh = 1.0d-12*(maxw)
             if(abs(d2wf1).lt.dwthresh) then
                rho = 0.0d0
             else
                rho = d2wlim/d2wf1
             endif
             if(rho .lt. (1.0d0 - 1.0d-12)) then
                applylim = .false.
             else
                !this stuff is equation 28
                d3wmin = min(d3wy(i,j),d3wy(i,j-1))
                d3wmin = min(d3wmin   ,d3wy(i,j+1))
                d3wmin = min(d3wmin   ,d3wy(i,j+2))
                d3wmax = max(d3wy(i,j),d3wy(i,j-1))
                d3wmax = max(d3wmax   ,d3wy(i,j+1))
                d3wmax = max(d3wmax   ,d3wy(i,j+2))
                dwthresh = 0.1d0*max(abs(d3wmax),abs(d3wmin))
                if((d3wmax - d3wmin).lt. dwthresh) then
                   applylim = .false.
                else
                   applylim = .true.
                endif
             endif
             !the following covered equations 29-32
             if(applylim) then
                if((dwfminuy(i,j)*dwfplusy(i,j)).lt.0.0) then

                   wrighy(i,j) = phi(i,j) -  rho*dwfminuy(i,j) !eqn 29 (corrected)
                   wlefty(i,j) = phi(i,j) +  rho*dwfplusy(i,j) !eqn 30

                else if(abs(dwfminuy(i,j)) .ge. 2.0d0*abs(dwfplusy(i,j))) then
                   !equation 31
                   wrighy(i,j) = phi(i,j) -  2.0d0*(1.0d0-rho)*dwfplusy(i,j) - rho*dwfminuy(i,j)
                else if(abs(dwfplusy(i,j)) .ge. 2.0d0*abs(dwfminuy(i,j))) then
                   !equation 32
                   wlefty(i,j) = phi(i,j) +  2.0d0*(1.0d0-rho)*dwfminuy(i,j) + rho*dwfplusy(i,j)
                endif
             endif

          else !monotonic but still need to check for big jumps in the state
             if(abs(dwfminuy(i,j)) .ge. (2.0d0*abs(dwfplusy(i,j)))) then
                wrighy(i,j) = phi(i,j) -  2.0d0*dwfplusy(i,j) !eqn 33
             endif
             if(abs(dwfplusx(i,j)) .ge. (2.0d0*abs(dwfminux(i,j)))) then
                wlefty(i,j) = phi(i,j) +  2.0d0*dwfminuy(i,j) !eqn 34
             endif
          endif

       enddo
    enddo
    !now get point valued phi at faces so we can multiply by point valued velocity
    ! phipt = phiave - (h^2/24)*(lapl^2d(phi_ave))
    ! while I am at it, multiply in pointwise velocity so we get pointwise *HYPERBOLIC* flux
    !also  get pointwise diffusive fluxes using 4th order finite differences of pointwise phi

    do    j = lo(2)-1, hi(2)+1
       do i = lo(1)  , hi(1)+1
          phiptx(i,j)  =    phiavex(i,j  ) - (1.0d0/24.d0)* &
               (phiavex(i,j+1) + phiavex(i,j-1) - 2.d0*phiavex(i,j))

          debtemp  = phiptx(i,j)

          diffflux  = (-nu/dx(1))* &
               ((27.0d0/24.0d0)*(phiptcc(i  ,j) - phiptcc(i-1,j)) &
               +( 1.0d0/24.0d0)*(phiptcc(i-2,j) - phiptcc(i+1,j)))

          fluxptx(i,j) = umac(i,j)*phiptx(i,j) + diffflux

!uncomment to just do diffusion
!          fluxptx(i,j) =  diffflux

       end do
    end do

    !same for y faces
    do    j = lo(2)  , hi(2)+1
       do i = lo(1)-1, hi(1)+1

          phipty(i,j)  = phiavey(i,j) - (1.0d0/24.d0)* &
               (phiavey(i+1,j)+phiavey(i-1,j)-2.d0*phiavey(i,j))

          debtemp  = phipty(i,j)

          diffflux  = (-nu/dx(2))* &
               ((27.0d0/24.0d0)*(phiptcc(i,j  ) - phiptcc(i,j-1)) &
               +( 1.0d0/24.0d0)*(phiptcc(i,j-2) - phiptcc(i,j+1)))

          fluxpty(i,j) = vmac(i,j)*phipty(i,j) + diffflux

!uncomment to just do diffusion
!          fluxpty(i,j) =  diffflux

       end do
    end do

    ! now transform pointwise fluxes into face-averaged fluxes
    !fluxave = fluxpt + (1/24)(Lapl2d(fluxpt))

    do    j = lo(2), hi(2)
       do i = lo(1), hi(1)+1
          debtemp  = fluxptx(i,j) + (1.0d0/24.d0)* &
               (fluxptx(i,j+1) + fluxptx(i,j-1) - 2.d0*fluxptx(i,j))

          flxx(i,j)  = debtemp

!debug   just set the flux = phiavex
!          flxx(i,j)  = phiavex(i,j)
!end debug
       end do
    end do

    do    j = lo(2), hi(2)+1
       do i = lo(1), hi(1)
          debtemp  = fluxpty(i,j) + (1.0d0/24.d0)* &
               (fluxpty(i+1,j) + fluxpty(i-1,j) - 2.d0*fluxpty(i,j))
          flxy(i,j)  = debtemp

!debug   just set the flux = phiavex
!          flxy(i,j)  = phiavey(i,j)
!end debug
       end do
    end do

    call bl_deallocate(dwfplusx)
    call bl_deallocate(dwfminux)
    call bl_deallocate(dwfplusy)
    call bl_deallocate(dwfminuy)
    call bl_deallocate(d2wfx   )
    call bl_deallocate(d2wcx   )
    call bl_deallocate(d3wx    )
    call bl_deallocate(d3wy    )
    call bl_deallocate(d2wfy   )
    call bl_deallocate(d2wcy   )
    call bl_deallocate(d3wy    )
    call bl_deallocate(wleftx  )
    call bl_deallocate(wlefty  )
    call bl_deallocate(wrighx  )
    call bl_deallocate(wrighy  )

 end subroutine mol4thord_flux_2d_limited

end module compute_flux_module


text:
program helloworld
  implicit none
  write(*,*) 'Hello world!'
end program helloworld

max_stars_repo_name: raramakr/aomp

max_issues_repo_name: raramakr/aomp

max_forks_repo_name: raramakr/aomp


text:
! { dg-do compile }

CONTAINS
SUBROUTINE send_forward ()

    INTEGER, DIMENSION(3)                    :: lz, ub, uz
    REAL, ALLOCATABLE, DIMENSION(:, :, :)    :: buffer
    COMPLEX, DIMENSION ( :, :, : ), POINTER  :: cc3d

    cc3d ( lz(1):uz(1), lz(2):uz(2), lz(3):uz(3) ) = &
           CMPLX ( buffer ( lz(1):uz(1), lz(2):uz(2), lz(3):uz(3) ), &
                   KIND = SELECTED_REAL_KIND ( 14, 200 ) )

END SUBROUTINE send_forward
END

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
c	Modified from: timedecm2.f
c	  (1) stupid questions removed GA 3/96
c	  (2) made a db-routine

c  sample time-domain deconvolution by damped least-squares inversion
c  this version does gauss-jordan
c  version does multiple stations simultaneously:  normalization by pwr in H
c --- GAA, 9/19/91
c     modified by afs 1/8/92

	parameter (NOUTMAX=512, NOUTMAX1=513)
	parameter (NDATMAX=2048, NBUFMAX=25000)
	dimension amat(NDATMAX,NOUTMAX),dat(NDATMAX)
c -- undisturbed copies of amat, dat:
	dimension amat2(NDATMAX,NOUTMAX), dat2(NDATMAX)
c --- thisone for dumping "A" matrix
c	dimension indump(NOUTMAX)
c	character*22 adump
c --- gauss-jordan I/O matricies
	dimension ata(NOUTMAX,NOUTMAX)
	dimension gjdat(NOUTMAX,NOUTMAX1)
	real*8 ataji, diagmax, gtd
c  resolution, error output matricies
	dimension resmat(NOUTMAX,NOUTMAX)
	dimension varout(NOUTMAX), resout(NOUTMAX)
c  spread-functions
	dimension sprout(NOUTMAX),sprbg(NOUTMAX)
c  Data buffers (!) -- from ahiof lib
        real*4 buf1(NBUFMAX),buf2(NBUFMAX),buf3(NBUFMAX)

C  DB stuff
	integer db(4), dbw(4), dbo(4), dbz(4), dbr(4), tr(4)
	integer wfid_z, wfid_r
	character*64 subset_z, subset_r
	character*6 sta2, sta1, instype
	character*8 chan2, chan1
	character*32 dfile
	character*64 dir
	character*2 datatype
	real*8 t0_req, t1_req, t0, t1, t0r,  srate1, srate2
	real*8 dtime, calib
	real zero
	integer nout, izero
	include 'db.i'

	character argu*80, dbout*80, iyesno*1
	external fill
c--------------
c	data luah1/9/,luah2/10/,luahout/11/
c --- for dumping "A"
c	data adump/"ZYXWVUTSR-0+123456789*"/
	data calib/1.d0/, dir/"wfrfcn"/, izero/0/, izero2/0/
	data instype/"TDrfcn"/

	ndatamax=NBUFMAX

	print *,'dbtimerf: Solves G*f=H for f, G=Greens fcn (Z-comp)'

	if (iargc().lt.6) then
	  write(0,*) "Usage: dbtimerf dbin dbout tshift tout apm eps
     &subset_z subset_r [kres]"
	  write(0,*) "tshift = shift of R rel. to G"
	  write(0,*) "tout = time length for output (should be <Ndata)"
	  write(0,*) "apm = prior modl uncert.= ap1/sqrt(NOUT) [1.] "
	  write(0,*) "eps = prior thry+data uncert.= eps*RMS(R) [1.]"
	  write(0,*) "subset_z, subset_r = dbselect logic strings for",
     -		" z,r component"
	  write(0,*) "kres = 1 to o/p resolution + spread, 0 otherwise"

	  stop
	endif

	call getarg(1, argu)
	ierr =  dbopen(argu(1:lnblnk(argu)), "r", db)
	if (ierr.eq.dbINVALID) call die (0,"Could not open db ")
	call getarg(2, dbout)
	ierr =  dbopen(dbout(1:lnblnk(dbout)), "r+", dbo)
	if (ierr.eq.dbINVALID) call die (0,"Could not open dbout ")
	call getarg(3, argu)
	read(argu,*) ppsec
	call getarg(4, argu)
	read(argu,*) tout
	call getarg(5,argu)
	read(argu,*) aprior1
	call getarg(6, argu)
	read(argu,*) eps
	call getarg(7, subset_z)
	call getarg(8, subset_r)
	if (iargc().gt.8) then
	  call getarg(9, argu)
	  read(argu,*) kresmat
	else
	  kresmat = 0
	endif

c---  inversion parameters (weighty weighting stuff)
	ppsec0 = ppsec
	dmom = aprior1
	if (eps.eq.0.) eps = 1.

c
c---open input wfdiscs
c
	call dblookup(dbw, db, "\0", "wfdisc", "\0", "\0" )
	call dbquery ( dbw, dbRECORD_COUNT, n )
	if (n.le.0) call die (0, " wfdisc empty or not found")
	call dbsubset(dbz, dbw, subset_z(1:lnblnk(subset_z)), "\0" )
	call dbquery ( dbz, dbRECORD_COUNT, nz )
	if (nz.le.0) call die (0, " z-subset empty or not found")
	call dbsubset(dbr, dbw, subset_r(1:lnblnk(subset_r)), "\0" )
	call dbquery ( dbr, dbRECORD_COUNT, nr )
	if (nr.le.0) call die (0, " r-subset empty or not found")
	if (nr.ne.nz) call die(0, " R and Z subsets don't match")
	nwfin = nz

c --- pointer to output db
	call dblookup(dbo, dbo, "\0", "wfdisc", "\0", "\0" )
	call dbquery ( dbo, dbRECORD_COUNT, n )
	write(0,*) n," records already in output db ",dbout
	nrecin = 0
	ndatin = 0
	negfin = 0
	sumwt = 0.
	sumegf = 0.
	ndata1 = 0
	ndata2 = 0

c	start reading in seismograms

	iquitflag = 0
	do 1 iwf=1,nwfin
	  if (ndatin.ge.NDATMAX.or.iquitflag.eq.1) then
	    print *,"too much input data; truncating at rec",nrecin
	    go to 2
	  end if

c--- read header G (Z)

	dbz(4) = iwf - 1
	if( dbgetv ( dbz, "",
     * 	  "sta", sta1, "chan", chan1,
     *    "time", t0_req, "samprate", srate1,
     *		"nsamp", ndata1, 0) .lt. 0)
     *    call complain (0, "dbgetv error" )

	   write (0,1000) sta1(1:lnblnk(sta1)),
     -		chan1(1:lnblnk(chan1)),srate1
1000	format( "Z station ",a6," chan: ",a8," samprate ",f8.3)

	ndata1=min(ndata1, NBUFMAX)
	t1_req=t0_req + dble(ndata1-1)/srate1

c---read data G (Z)
c

c	call trgetwf ( dbz, 0, buf1, ndatamax, t0_req, t1_req,
c     *         t0, t1, npts, izero, izero2 )
	call readmywf(dbz, t0_req, t1_req, tr, t0, t1,
     *		npts, buf1)

	if ( npts .lt. 1 ) call die ( 0, "Couldnt read Zwaveform")
	write(0,1002) npts, t0
1002	format(" Z: read in ",i6," data points starting ",f14.3)
	ndata1=npts

c--- read header H (R)
	dbr(4) = iwf - 1
	if( dbgetv ( dbr, "",
     * 	  "sta", sta2, "chan", chan2,
     *    "time", t0_req, "samprate", srate2,
     *		"nsamp", ndata2, 0) .lt. 0)
     *    call complain (0, "dbgetv 2 error" )

	   write (0,1001) sta2(1:lnblnk(sta2)),
     -		chan2(1:lnblnk(chan2)),srate2
1001	format( "R station ",a6," chan: ",a8," samprate ",f8.3)

	ndata2=min(ndata2, NBUFMAX)
	t1_req=t0_req + dble(ndata2-1)/srate2
c---read data H (R)
c
c	call trgetwf ( dbr, 0, buf2, ndatamax, t0_req, t1_req,
c     *         t0r, t1, npts, izero, izero2 )
	call readmywf(dbr, t0_req, t1_req, tr, t0r, t1,
     *		npts, buf2)
	if ( npts .lt. 1 ) call die ( 0, "Couldnt read waveform")
	write(0,1003) npts, t0r
1003	format(" R: read in ",i6," data points starting ",f14.3)
	ndata2=npts

c
c--- check to see if things are similar enough to deconvolve

	if (srate2.ne.srate1) then
	  print *,'Samprates dont match; skipping'
	  print *,'delta-G = ',srate1,' delta-H = ',srate2
	  go to 1
	end if
	delta2 = 1./srate2
	if (dabs(t0-t0r).gt.delta2) then
	  print *,'Start times dont match; skipping'
	  print *,'t0-z = ',t0,' t0-r = ',t0r
	  go to 1
	end if

c--- get rms(h) and data weight
	rmsh=0.
	sumh=0.
	do i=1,ndata2
	  sumh=sumh+buf2(i)
	  rmsh=rmsh+buf2(i)*buf2(i)
	end do
	if (ndata2.lt.2) then
	  rmsh = abs(buf2(1))
	else
	  en=ndata2
	  rmsh=(rmsh-sumh*sumh/en)/(en-1.)
	  if (rmsh.lt.0.) rmsh=1.
	  rmsh=sqrt(rmsh)
	end if
	datwt = 1./(rmsh*eps)

c --- pre-pad H, if necessary
	if (ppsec.gt.0.) then
  	  npad = ppsec/delta2
	  t0r = t0r - npad*delta2
	  nout2 = ndata2+npad
	  do i=ndata2,1,-1
	    buf2(i+npad) = buf2(i)
	  end do
	  do i=1,npad
	    buf2(i) = 0.
	  end do
	  ndata2 = nout2
	end if

c --- check for too much data
	if (ndata2+ndatin.ge.NDATMAX) then
	  print *,'Too much data; truncating to ',NDATMAX
	  ndata2 = NDATMAX - ndatin-1
	  iquitflag = 1
	end if
	if (ndata1+ndatin.ge.NDATMAX) ndata1=NDATMAX-ndatin-1

c ---- set number of output points, on first seismogram
	if (nrecin.eq.0) then
	  nout = tout*srate2 + 1
	  if (nout.gt.NOUTMAX) nout=NOUTMAX
	  print *,nout,' points for output'
	end if

c --- build matrix to decompose (weight by rms-H)
c --- NOTE:  this should directly calculate aTa and gjdat; it would
c		be a lot more efficient

	do j=1,nout
	    if (j.gt.1.and.ndata2.gt.1) then
		jm1 = j - 1
	 	do i=1,jm1
		  ii = i + ndatin
		  if (ii.le.NDATMAX) amat(ii,j) = 0.
		end do
	    end if
	    n2 = j + ndata1 - 1
	    if (n2.gt.ndata2) n2 = ndata2
	    do i=j,n2
	      ii = i + ndatin
	      if (ii.le.NDATMAX) amat(ii,j) = buf1(i-j+1)*datwt
	    end do
	    if (n2.lt.ndata2) then
		do i=n2,ndata2
		  ii = i + ndatin
		  if (ii.le.NDATMAX) amat(ii,j) = 0.
		end do
	    end if
	end do
c--- get pwr in g*wt, to check stability
	rmsg = 0.
	do i=1,ndata1
	  rmsg=rmsg+buf1(i)*buf1(i)*datwt*datwt
	end do

c --- build data array
	do i=1,ndata2
	  ii = i + ndatin
	  dat(ii) = buf2(i)*datwt
	end do

	sumwt = sumwt + float(ndata2-1)*datwt*datwt
	sumegf = sumegf + rmsg
	negfin = negfin + ndata1
	ndatin = ndatin + ndata2
	nrecin = nrecin + 1
1	continue

c -- END OF DATA INPUT LOOP
2	continue
	wtavg = sqrt(sumwt/float(ndatin))
	egfavg = sqrt(sumegf/float(negfin))
	print *,"Average weight = ",wtavg," rms of weighted A-mat=",egfavg

	dmom = dmom/sqrt(real(nout))
	print *,"Est. prior uncert. in rfcn = ",dmom

c -- save input matricies
	do i=1,ndatin
	  dat2(i) = dat(i)
	  do j=1,nout
	    amat2(i,j) = amat(i,j)
	  end do
	end do
c
c--- calculate Cmm-1
	  cmm1 = 0.
	  if (dmom.gt.0.) cmm1=(1./dmom)**2
	  print *,"inverse variance diag Cmm-1 = ",cmm1,
     -		" A priori error = ",dmom," root-damper=",1/dmom

c--- do gauss-jordan elimination
	  print *,"Making AtA+Cm-1 and At[d|AtA] matricies....."
	  nosvd = 1

c--- build g-j matricies: ata = normal eq'ns, gjdat=[d|AtA] = data
	  diagmax = 0.
	  do i=1,nout
	    do j=1,i
	      ataji = 0.d0
	      do k=1,ndatin
		ataji=ataji+amat2(k,i)*amat2(k,j)
	      end do
	      if (kresmat.eq.1) then
	 	gjdat(j,i+1)=ataji
		gjdat(i,j+1) = ataji
	      end if
	      ata(j,i)=ataji
	      ata(i,j)=ataji
	      if (j.eq.i) ata(j,i) = ataji + cmm1
	      if (i.eq.j.and.ataji.gt.diagmax) diagmax=ataji
	    end do
	    gtd = 0.d0
	    do k=1,ndatin
	      gtd=gtd+dble(amat2(k,i))*dble(dat(k))
	    end do
	    gjdat(i,1) = gtd
c	    gjdat1(i) = gtd
	  end do
	  print *,'AtA: max-diag=',diagmax,' effective CN=',
     -		diagmax/dble(cmm1)+1.d0
c-- dump AtA matrix
c	open (41,file="a.dump")
c	do i=1,nout
c	  do j=1,nout
c	    ii = int((ata(i,j)/diagmax)*10.d0)+11
c	    if (ata(i,j).lt.0.d0) ii=ii-1
c	    if (ata(i,j).gt.0.d0) ii=ii+1
c	    if (ii.lt.1) ii=1
c	    if (ii.gt.22) ii=22
c	    indump(j) = ii
c	  end do
c	  write(41,"(132a1)") (adump(indump(j):indump(j)),j=1,nout)
c	end do
c	close(41)

	  nout1 = 1
	  if (kresmat.eq.1) nout1=nout+1

c--- invert using numerical recipes subroutine [SVD encountered ERRORS]
c  A note on speed:  This method should be comparable to all others as long
c  as the entire inverse matrix ata-1 is desired.  If complete resolution and error
c  info is not desired, this may not be necessary.  ALso, it will crash if singular.

	  call gaussj(ata,nout,NOUTMAX,gjdat,nout1,NOUTMAX1,infogj)

	  if (infogj.ne.0) then
	    print *,"GAUSSJ failed, status = ",info," **STOPPING**"
	    stop
	  else
	    print *,"GAUSSJ done, successful"
	  end if

c--- do gauss-jordan solution
	  sumbuf = 0.
	  do i=1,nout
	    buf3(i) = gjdat(i,1)
	    sumbuf = sumbuf + buf3(i)
	  end do
	print *,"Integrated STF = ",sumbuf

c--- get errors, possibly resolutions, from Gauss-jordan
c     using Total or a posteriori error estimate = inv(AtA+inv(Cm))
	  resmax = 0.
	  varmax = 0.
	  do i=1,nout
	    varout(i) = (sqrt(ata(i,i)))
	    if (varout(i).gt.varmax) varmax = varout(i)
	  end do
	  if (kresmat.eq.1) then
	    do i=1,nout
	      ip1 = i+1
	      resout(i) = (gjdat(i,ip1))
	      if (resout(i).gt.resmax) resmax=resout(i)
	      do j=1,nout
		resmat(j,i) = (gjdat(j,ip1))
	      end do
	    end do
	  end if

c--- calculate spread functions
	if (kresmat.eq.1) then
	  do i=1,nout
	    spread1 = 0.
	    spread2 = 0.
	    spread3 = 0.
	    sprout(i) = 0.
	    sprbg(i) = 0.
	    do k=1,nout
	      spread1=spread1+resmat(i,k)
	      ss2 = resmat(i,k)**2
	      spread2=spread2+ss2
	      spread3=spread3+ss2*float((i-k)*(i-k))
	    end do
	    if (spread2.gt.0.) then
c --- this spread-measure works well for sinx/x distributions
	      sprout(i)=spread1*spread1/spread2*delta2
c --- this is a Backus-Gilbert-type measure
	      sprbg(i)=sqrt(spread3/spread2)*delta2
	    end if
	  end do
	end if

c--- calculate residuals
	write (*,"('Calculating residuals...')")
	nd = ndatin
	ressum = 0.
	datsum = 0.
	do i=1,nd
	  resi = 0.
	  do j=1,nout
	    resi=resi+amat2(i,j)*buf3(j)
	  end do
	  resi = dat2(i)-resi
	  ressum = ressum + resi*resi
	  datsum = datsum + dat2(i)*dat2(i)
	end do
	ressum = ressum/sumwt
	datsum = datsum/sumwt
	fracout = ressum/datsum
	sigdav = sqrt(float(nd)/sumwt)
	print *," Final RMS Residual              = ",sqrt(ressum)
     	print *," var(res)/var(dat)               = ",fracout
  	print *," mean a priori data+theory error = ",sigdav

c correct sign for receiver-function and rotation conventions....
	do i=1,nout
	    buf3(i) = -buf3(i)
	enddo
c  write header

	call dblookup(dbo, dbo, "\0", "wfdisc", "\0", "\0" )
	chan2="tdrf\0"
	call e2h(t0, iyr, idy, ihr, imn, sec)
	dtime=(dble(iyr*1000+idy)*100d0+dble(ihr))*100.d0+dble(imn)
	dtime=dtime+dble(int(sec))/60.
	write(dfile,'(f14.2,".",a3,".",a4)') dtime,sta2(1:3),chan2(1:4)

	call writemywf(dbo,  sta2, chan2,
     -		t0r, srate2, nout, calib,
     -		instype, dir, dfile, buf3)

c--- write standard error
	chan2="terr\0"
	write(dfile,'(f14.2,".",a3,".",a4)') dtime,sta2(1:3),chan2(1:4)
	call writemywf(dbo,  sta2, chan2,
     -		t0r, srate2, nout, calib,
     -		instype, dir, dfile, varout)

c--- write resolution diagonals
	if (kresmat.eq.1) then
	  print *,"res-matrix output to timedec.resmat"
	  chan2="res\0"
	write(dfile,'(f14.2,".",a3,".",a3)') dtime,sta2(1:3),chan2(1:3)
	call writemywf(dbo,  sta2, chan2,
     -		t0r, srate2, nout, calib,
     -		instype, dir, dfile, resout)

c--- write resolution spread widths (sinc)
	  chan2="sprd\0"
	write(dfile,'(f14.2,".",a3,".",a3)') dtime,sta2(1:3),chan2(1:4)
	call writemywf(dbo,  sta2, chan2,
     -		t0r, srate2, nout, calib,
     -		instype, dir, dfile, sprout)

c--- write Backus-Gilbert spread
	  chan2="BGsp\0"
	write(dfile,'(f14.2,".",a3,".",a3)') dtime,sta2(1:3),chan2(1:4)
	call writemywf(dbo,  sta2, chan2,
     -		t0r, srate2, nout, calib,
     -		instype, dir, dfile, sprbg)

c -- output full resolution matrix in rmm-format (ASCII)
	  open(22,file="timedec.resmat")
	  nelem=nout*(nout+1)/2
	  write(22,"(5i6,f10.3)") 1,1,0,2,nelem,1.
	  write(22,"(10e11.3)") ((resmat(j,i),j=1,i),i=1,nout)
	  close(22)
	endif

	stop
	end

          subroutine fill ( db, data, i0, i1, imax, value )
          implicit none
          integer db(4)
          integer i, i0, i1, imax
          real data(imax)
          real value
           do i=i0+1, i1
             data(i) = value
	   end do
	end


text:
program tst
use constants
implicit none
integer, parameter :: dp = selected_real_kind(15,307)
! program to calculate tst rate
real (dp) :: deltag,rate,temp,deg
! read deltag (in kcal/mol) and temperature (in K)
read(*,*) deltag,temp,deg
rate=deg*boltz*temp/planck*exp(-deltag/r/temp)
! rate in s-1
print*, rate

end program tst

max_stars_repo_name: emartineznunez/AutoMeKin

max_issues_repo_name: emartineznunez/AutoMeKin

max_forks_repo_name: emartineznunez/AutoMeKin


text:
program main
    use omp_lib
    integer :: nteams, nthreads

    !$omp target map(nteams, nthreads)
        nteams = omp_get_num_teams()
        nthreads = omp_get_num_threads()
    !$omp end target

    write(*,*) nteams, nthreads
    if (nthreads .ne. 1) write (*,*) "Failed nthreads for gneric kernell"
end program main

max_stars_repo_name: aperdeus-amd/aomp

max_issues_repo_name: aperdeus-amd/aomp

max_forks_repo_name: aperdeus-amd/aomp


text:
      SUBROUTINE CTLOOP(NLOOPLINE,PL,M2L,RANK,RES,STABLE)
C
C     Generated by MadGraph5_aMC@NLO v. %(version)s, %(date)s
C     By the MadGraph5_aMC@NLO Development Team
C     Visit launchpad.net/madgraph5 and amcatnlo.web.cern.ch
C
C     Interface between MG5 and CutTools.
C
C     Process: g g > w- t b~ QED=1 QCD=2 [ virt = QCD ] WEIGHTED=10
C
C
C     CONSTANTS
C
      INTEGER    NEXTERNAL
      PARAMETER (NEXTERNAL=5)
      LOGICAL CHECKPCONSERVATION
      PARAMETER (CHECKPCONSERVATION=.TRUE.)
      REAL*8 NORMALIZATION
      PARAMETER (NORMALIZATION = 1.D0/(16.D0*3.14159265358979323846D0*
     $ *2))
C
C     ARGUMENTS
C
      INTEGER NLOOPLINE, RANK
      REAL*8 PL(0:3,NLOOPLINE)
      REAL*8 PCT(0:3,0:NLOOPLINE-1)
      COMPLEX*16 M2L(NLOOPLINE)
      COMPLEX*16 M2LCT(0:NLOOPLINE-1)
      COMPLEX*16 RES(3)
      LOGICAL STABLE
C
C     LOCAL VARIABLES
C
      COMPLEX*16 R1, ACC
      INTEGER I, J, K
      LOGICAL INIT
      DATA INIT /.TRUE./
      SAVE INIT
C
C     EXTERNAL FUNCTIONS
C
      EXTERNAL LOOPNUM
      EXTERNAL MPLOOPNUM
C
C     GLOBAL VARIABLES
C
      INCLUDE 'coupl.inc'
      INTEGER CTMODE
      REAL*8 LSCALE
      COMMON/CT/LSCALE,CTMODE

      INTEGER ID,R
      COMMON/LOOP/ID,R

C     ----------
C     BEGIN CODE
C     ----------

C     INITIALIZE CUTTOOLS IF NEEDED
      IF (INIT) THEN
        INIT=.FALSE.
        CALL INITCT()
      ENDIF

C     YOU CAN FIND THE DETAILS ABOUT THE DIFFERENT CTMODE AT THE
C      BEGINNING OF THE FILE CTS_CUTS.F90 IN THE CUTTOOLS DISTRIBUTION

C     CONVERT THE MASSES TO BE COMPLEX
      DO I=1,NLOOPLINE
        M2LCT(I-1)=M2L(I)
      ENDDO

C     CONVERT THE MOMENTA FLOWING IN THE LOOP LINES TO CT CONVENTIONS
      DO I=0,3
        DO J=0,(NLOOPLINE-1)
          PCT(I,J)=0.D0
        ENDDO
      ENDDO
      DO I=0,3
        DO J=1,NLOOPLINE
          PCT(I,0)=PCT(I,0)+PL(I,J)
        ENDDO
      ENDDO
      IF (CHECKPCONSERVATION) THEN
        IF (PCT(0,0).GT.1.D-6) THEN
          WRITE(*,*) 'energy is not conserved ',PCT(0,0)
          STOP 'energy is not conserved'
        ELSEIF (PCT(1,0).GT.1.D-6) THEN
          WRITE(*,*) 'px is not conserved ',PCT(1,0)
          STOP 'px is not conserved'
        ELSEIF (PCT(2,0).GT.1.D-6) THEN
          WRITE(*,*) 'py is not conserved ',PCT(2,0)
          STOP 'py is not conserved'
        ELSEIF (PCT(3,0).GT.1.D-6) THEN
          WRITE(*,*) 'pz is not conserved ',PCT(3,0)
          STOP 'pz is not conserved'
        ENDIF
      ENDIF
      DO I=0,3
        DO J=1,(NLOOPLINE-1)
          DO K=1,J
            PCT(I,J)=PCT(I,J)+PL(I,K)
          ENDDO
        ENDDO
      ENDDO

      CALL CTSXCUT(CTMODE,LSCALE,MU_R,NLOOPLINE,LOOPNUM,MPLOOPNUM,RANK
     $ ,PCT,M2LCT,RES,ACC,R1,STABLE)
      RES(1)=NORMALIZATION*2.0D0*DBLE(RES(1))
      RES(2)=NORMALIZATION*2.0D0*DBLE(RES(2))
      RES(3)=NORMALIZATION*2.0D0*DBLE(RES(3))
C     WRITE(*,*) 'Loop ID',ID,' =',RES(1),RES(2),RES(3)
      END

      SUBROUTINE INITCT()
C
C     INITIALISATION OF CUTTOOLS
C
C     LOCAL VARIABLES
C
      REAL*8 THRS
      LOGICAL EXT_NUM_FOR_R1
C
C     GLOBAL VARIABLES
C
      INCLUDE 'MadLoopParams.inc'
C     ----------
C     BEGIN CODE
C     ----------

C     DEFAULT PARAMETERS FOR CUTTOOLS
C     -------------------------------
C     THRS1 IS THE PRECISION LIMIT BELOW WHICH THE MP ROUTINES
C      ACTIVATES
      THRS=CTSTABTHRES
C     LOOPLIB SET WHAT LIBRARY CT USES
C     1 -> LOOPTOOLS
C     2 -> AVH
C     3 -> QCDLOOP
      LOOPLIB=CTLOOPLIBRARY
C     MADLOOP'S NUMERATOR IN THE OPEN LOOP IS MUCH FASTER THAN THE
C      RECONSTRUCTED ONE IN CT. SO WE BETTER USE MADLOOP ONE IN THIS
C      CASE.
      EXT_NUM_FOR_R1=.TRUE.
C     -------------------------------

C     The initialization below is for CT v1.8.+
      CALL CTSINIT(THRS,LOOPLIB,EXT_NUM_FOR_R1)
C     The initialization below is for the older stable CT v1.7, still
C      used for now in the beta release.
C     CALL CTSINIT(THRS,LOOPLIB)

      END

      SUBROUTINE LOOP_5( LID, W1, W2, W3, W4, W5, M1, M2, M3, M4, M5
     $ ,  RANK, RES, STABLE, LOOPNUM)

      INTEGER    NEXTERNAL
      PARAMETER (NEXTERNAL=5)
      INTEGER    NLOOPLINE
      PARAMETER (NLOOPLINE=5)
      INTEGER    NWAVEFUNCS
      PARAMETER (NWAVEFUNCS=28)
      INTEGER    NLOOPGROUPS
      PARAMETER (NLOOPGROUPS=77)
      INTEGER    NCOMB
      PARAMETER (NCOMB=48)
C
C     ARGUMENTS
C
      INTEGER W1, W2, W3, W4, W5
      COMPLEX*16 M1, M2, M3, M4, M5

      COMPLEX*16 RES(3)
      INTEGER LID, RANK, LSYMFACT
      INTEGER LOOPNUM
      LOGICAL STABLE
C
C     LOCAL VARIABLES
C
      REAL*8 PL(0:3,NLOOPLINE)
      COMPLEX*16 M2L(NLOOPLINE)
      INTEGER PAIRING(NLOOPLINE),WE(5)
      INTEGER I, J, K, TEMP
C
C     GLOBAL VARIABLES
C
      INTEGER ID,R
      COMMON/LOOP/ID,R

      LOGICAL CHECKPHASE, HELDOUBLECHECKED
      COMMON/INIT/CHECKPHASE, HELDOUBLECHECKED

      INTEGER HELOFFSET
      INTEGER GOODHEL(NCOMB)
      LOGICAL GOODAMP(NLOOPGROUPS)
      COMMON/FILTERS/GOODAMP,GOODHEL,HELOFFSET

      COMPLEX*16 W(20,NWAVEFUNCS)
      COMMON/W/W

C     ----------
C     BEGIN CODE
C     ----------

      IF (CHECKPHASE.OR.(.NOT.HELDOUBLECHECKED).OR.GOODAMP(LOOPNUM)
     $ ) THEN
        WE(1)=W1
        WE(2)=W2
        WE(3)=W3
        WE(4)=W4
        WE(5)=W5
        M2L(1)=M5**2
        M2L(2)=M1**2
        M2L(3)=M2**2
        M2L(4)=M3**2
        M2L(5)=M4**2
        DO I=1,NLOOPLINE
          PAIRING(I)=1
        ENDDO

        R=RANK
        ID=LID
        DO I=0,3
          TEMP=1
          DO J=1,NLOOPLINE
            PL(I,J)=0.D0
            DO K=TEMP,(TEMP+PAIRING(J)-1)
              PL(I,J)=PL(I,J)-DBLE(W(1+I,WE(K)))
            ENDDO
            TEMP=TEMP+PAIRING(J)
          ENDDO
        ENDDO
        CALL CTLOOP(NLOOPLINE,PL,M2L,RANK,RES,STABLE)
      ELSE
        RES(1)=0.0D0
        RES(2)=0.0D0
        RES(3)=0.0D0
        STABLE=.TRUE.
      ENDIF
      END

      SUBROUTINE LOOP_4_5( LID, P1, P2, P3, P4, W1, W2, W3, W4, W5, M1
     $ , M2, M3, M4,  RANK, RES, STABLE, LOOPNUM)

      INTEGER    NEXTERNAL
      PARAMETER (NEXTERNAL=5)
      INTEGER    NLOOPLINE
      PARAMETER (NLOOPLINE=4)
      INTEGER    NWAVEFUNCS
      PARAMETER (NWAVEFUNCS=28)
      INTEGER    NLOOPGROUPS
      PARAMETER (NLOOPGROUPS=77)
      INTEGER    NCOMB
      PARAMETER (NCOMB=48)
C
C     ARGUMENTS
C
      INTEGER W1, W2, W3, W4, W5
      COMPLEX*16 M1, M2, M3, M4
      INTEGER P1, P2, P3, P4
      COMPLEX*16 RES(3)
      INTEGER LID, RANK, LSYMFACT
      INTEGER LOOPNUM
      LOGICAL STABLE
C
C     LOCAL VARIABLES
C
      REAL*8 PL(0:3,NLOOPLINE)
      COMPLEX*16 M2L(NLOOPLINE)
      INTEGER PAIRING(NLOOPLINE),WE(5)
      INTEGER I, J, K, TEMP
C
C     GLOBAL VARIABLES
C
      INTEGER ID,R
      COMMON/LOOP/ID,R

      LOGICAL CHECKPHASE, HELDOUBLECHECKED
      COMMON/INIT/CHECKPHASE, HELDOUBLECHECKED

      INTEGER HELOFFSET
      INTEGER GOODHEL(NCOMB)
      LOGICAL GOODAMP(NLOOPGROUPS)
      COMMON/FILTERS/GOODAMP,GOODHEL,HELOFFSET

      COMPLEX*16 W(20,NWAVEFUNCS)
      COMMON/W/W

C     ----------
C     BEGIN CODE
C     ----------

      IF (CHECKPHASE.OR.(.NOT.HELDOUBLECHECKED).OR.GOODAMP(LOOPNUM)
     $ ) THEN
        WE(1)=W1
        WE(2)=W2
        WE(3)=W3
        WE(4)=W4
        WE(5)=W5
        M2L(1)=M4**2
        M2L(2)=M1**2
        M2L(3)=M2**2
        M2L(4)=M3**2
        PAIRING(1)=P1
        PAIRING(2)=P2
        PAIRING(3)=P3
        PAIRING(4)=P4
        R=RANK
        ID=LID
        DO I=0,3
          TEMP=1
          DO J=1,NLOOPLINE
            PL(I,J)=0.D0
            DO K=TEMP,(TEMP+PAIRING(J)-1)
              PL(I,J)=PL(I,J)-DBLE(W(1+I,WE(K)))
            ENDDO
            TEMP=TEMP+PAIRING(J)
          ENDDO
        ENDDO
        CALL CTLOOP(NLOOPLINE,PL,M2L,RANK,RES,STABLE)
      ELSE
        RES(1)=0.0D0
        RES(2)=0.0D0
        RES(3)=0.0D0
        STABLE=.TRUE.
      ENDIF
      END

      SUBROUTINE LOOP_4( LID, W1, W2, W3, W4, M1, M2, M3, M4,  RANK
     $ , RES, STABLE, LOOPNUM)

      INTEGER    NEXTERNAL
      PARAMETER (NEXTERNAL=5)
      INTEGER    NLOOPLINE
      PARAMETER (NLOOPLINE=4)
      INTEGER    NWAVEFUNCS
      PARAMETER (NWAVEFUNCS=28)
      INTEGER    NLOOPGROUPS
      PARAMETER (NLOOPGROUPS=77)
      INTEGER    NCOMB
      PARAMETER (NCOMB=48)
C
C     ARGUMENTS
C
      INTEGER W1, W2, W3, W4
      COMPLEX*16 M1, M2, M3, M4

      COMPLEX*16 RES(3)
      INTEGER LID, RANK, LSYMFACT
      INTEGER LOOPNUM
      LOGICAL STABLE
C
C     LOCAL VARIABLES
C
      REAL*8 PL(0:3,NLOOPLINE)
      COMPLEX*16 M2L(NLOOPLINE)
      INTEGER PAIRING(NLOOPLINE),WE(4)
      INTEGER I, J, K, TEMP
C
C     GLOBAL VARIABLES
C
      INTEGER ID,R
      COMMON/LOOP/ID,R

      LOGICAL CHECKPHASE, HELDOUBLECHECKED
      COMMON/INIT/CHECKPHASE, HELDOUBLECHECKED

      INTEGER HELOFFSET
      INTEGER GOODHEL(NCOMB)
      LOGICAL GOODAMP(NLOOPGROUPS)
      COMMON/FILTERS/GOODAMP,GOODHEL,HELOFFSET

      COMPLEX*16 W(20,NWAVEFUNCS)
      COMMON/W/W

C     ----------
C     BEGIN CODE
C     ----------

      IF (CHECKPHASE.OR.(.NOT.HELDOUBLECHECKED).OR.GOODAMP(LOOPNUM)
     $ ) THEN
        WE(1)=W1
        WE(2)=W2
        WE(3)=W3
        WE(4)=W4
        M2L(1)=M4**2
        M2L(2)=M1**2
        M2L(3)=M2**2
        M2L(4)=M3**2
        DO I=1,NLOOPLINE
          PAIRING(I)=1
        ENDDO

        R=RANK
        ID=LID
        DO I=0,3
          TEMP=1
          DO J=1,NLOOPLINE
            PL(I,J)=0.D0
            DO K=TEMP,(TEMP+PAIRING(J)-1)
              PL(I,J)=PL(I,J)-DBLE(W(1+I,WE(K)))
            ENDDO
            TEMP=TEMP+PAIRING(J)
          ENDDO
        ENDDO
        CALL CTLOOP(NLOOPLINE,PL,M2L,RANK,RES,STABLE)
      ELSE
        RES(1)=0.0D0
        RES(2)=0.0D0
        RES(3)=0.0D0
        STABLE=.TRUE.
      ENDIF
      END

      SUBROUTINE LOOP_2_3( LID, P1, P2, W1, W2, W3, M1, M2,  RANK, RES
     $ , STABLE, LOOPNUM)

      INTEGER    NEXTERNAL
      PARAMETER (NEXTERNAL=5)
      INTEGER    NLOOPLINE
      PARAMETER (NLOOPLINE=2)
      INTEGER    NWAVEFUNCS
      PARAMETER (NWAVEFUNCS=28)
      INTEGER    NLOOPGROUPS
      PARAMETER (NLOOPGROUPS=77)
      INTEGER    NCOMB
      PARAMETER (NCOMB=48)
C
C     ARGUMENTS
C
      INTEGER W1, W2, W3
      COMPLEX*16 M1, M2
      INTEGER P1, P2
      COMPLEX*16 RES(3)
      INTEGER LID, RANK, LSYMFACT
      INTEGER LOOPNUM
      LOGICAL STABLE
C
C     LOCAL VARIABLES
C
      REAL*8 PL(0:3,NLOOPLINE)
      COMPLEX*16 M2L(NLOOPLINE)
      INTEGER PAIRING(NLOOPLINE),WE(3)
      INTEGER I, J, K, TEMP
C
C     GLOBAL VARIABLES
C
      INTEGER ID,R
      COMMON/LOOP/ID,R

      LOGICAL CHECKPHASE, HELDOUBLECHECKED
      COMMON/INIT/CHECKPHASE, HELDOUBLECHECKED

      INTEGER HELOFFSET
      INTEGER GOODHEL(NCOMB)
      LOGICAL GOODAMP(NLOOPGROUPS)
      COMMON/FILTERS/GOODAMP,GOODHEL,HELOFFSET

      COMPLEX*16 W(20,NWAVEFUNCS)
      COMMON/W/W

C     ----------
C     BEGIN CODE
C     ----------

      IF (CHECKPHASE.OR.(.NOT.HELDOUBLECHECKED).OR.GOODAMP(LOOPNUM)
     $ ) THEN
        WE(1)=W1
        WE(2)=W2
        WE(3)=W3
        M2L(1)=M2**2
        M2L(2)=M1**2
        PAIRING(1)=P1
        PAIRING(2)=P2
        R=RANK
        ID=LID
        DO I=0,3
          TEMP=1
          DO J=1,NLOOPLINE
            PL(I,J)=0.D0
            DO K=TEMP,(TEMP+PAIRING(J)-1)
              PL(I,J)=PL(I,J)-DBLE(W(1+I,WE(K)))
            ENDDO
            TEMP=TEMP+PAIRING(J)
          ENDDO
        ENDDO
        CALL CTLOOP(NLOOPLINE,PL,M2L,RANK,RES,STABLE)
      ELSE
        RES(1)=0.0D0
        RES(2)=0.0D0
        RES(3)=0.0D0
        STABLE=.TRUE.
      ENDIF
      END

      SUBROUTINE LOOP_2( LID, W1, W2, M1, M2,  RANK, RES, STABLE
     $ , LOOPNUM)

      INTEGER    NEXTERNAL
      PARAMETER (NEXTERNAL=5)
      INTEGER    NLOOPLINE
      PARAMETER (NLOOPLINE=2)
      INTEGER    NWAVEFUNCS
      PARAMETER (NWAVEFUNCS=28)
      INTEGER    NLOOPGROUPS
      PARAMETER (NLOOPGROUPS=77)
      INTEGER    NCOMB
      PARAMETER (NCOMB=48)
C
C     ARGUMENTS
C
      INTEGER W1, W2
      COMPLEX*16 M1, M2

      COMPLEX*16 RES(3)
      INTEGER LID, RANK, LSYMFACT
      INTEGER LOOPNUM
      LOGICAL STABLE
C
C     LOCAL VARIABLES
C
      REAL*8 PL(0:3,NLOOPLINE)
      COMPLEX*16 M2L(NLOOPLINE)
      INTEGER PAIRING(NLOOPLINE),WE(2)
      INTEGER I, J, K, TEMP
C
C     GLOBAL VARIABLES
C
      INTEGER ID,R
      COMMON/LOOP/ID,R

      LOGICAL CHECKPHASE, HELDOUBLECHECKED
      COMMON/INIT/CHECKPHASE, HELDOUBLECHECKED

      INTEGER HELOFFSET
      INTEGER GOODHEL(NCOMB)
      LOGICAL GOODAMP(NLOOPGROUPS)
      COMMON/FILTERS/GOODAMP,GOODHEL,HELOFFSET

      COMPLEX*16 W(20,NWAVEFUNCS)
      COMMON/W/W

C     ----------
C     BEGIN CODE
C     ----------

      IF (CHECKPHASE.OR.(.NOT.HELDOUBLECHECKED).OR.GOODAMP(LOOPNUM)
     $ ) THEN
        WE(1)=W1
        WE(2)=W2
        M2L(1)=M2**2
        M2L(2)=M1**2
        DO I=1,NLOOPLINE
          PAIRING(I)=1
        ENDDO

        R=RANK
        ID=LID
        DO I=0,3
          TEMP=1
          DO J=1,NLOOPLINE
            PL(I,J)=0.D0
            DO K=TEMP,(TEMP+PAIRING(J)-1)
              PL(I,J)=PL(I,J)-DBLE(W(1+I,WE(K)))
            ENDDO
            TEMP=TEMP+PAIRING(J)
          ENDDO
        ENDDO
        CALL CTLOOP(NLOOPLINE,PL,M2L,RANK,RES,STABLE)
      ELSE
        RES(1)=0.0D0
        RES(2)=0.0D0
        RES(3)=0.0D0
        STABLE=.TRUE.
      ENDIF
      END

      SUBROUTINE LOOP_3_4( LID, P1, P2, P3, W1, W2, W3, W4, M1, M2, M3
     $ ,  RANK, RES, STABLE, LOOPNUM)

      INTEGER    NEXTERNAL
      PARAMETER (NEXTERNAL=5)
      INTEGER    NLOOPLINE
      PARAMETER (NLOOPLINE=3)
      INTEGER    NWAVEFUNCS
      PARAMETER (NWAVEFUNCS=28)
      INTEGER    NLOOPGROUPS
      PARAMETER (NLOOPGROUPS=77)
      INTEGER    NCOMB
      PARAMETER (NCOMB=48)
C
C     ARGUMENTS
C
      INTEGER W1, W2, W3, W4
      COMPLEX*16 M1, M2, M3
      INTEGER P1, P2, P3
      COMPLEX*16 RES(3)
      INTEGER LID, RANK, LSYMFACT
      INTEGER LOOPNUM
      LOGICAL STABLE
C
C     LOCAL VARIABLES
C
      REAL*8 PL(0:3,NLOOPLINE)
      COMPLEX*16 M2L(NLOOPLINE)
      INTEGER PAIRING(NLOOPLINE),WE(4)
      INTEGER I, J, K, TEMP
C
C     GLOBAL VARIABLES
C
      INTEGER ID,R
      COMMON/LOOP/ID,R

      LOGICAL CHECKPHASE, HELDOUBLECHECKED
      COMMON/INIT/CHECKPHASE, HELDOUBLECHECKED

      INTEGER HELOFFSET
      INTEGER GOODHEL(NCOMB)
      LOGICAL GOODAMP(NLOOPGROUPS)
      COMMON/FILTERS/GOODAMP,GOODHEL,HELOFFSET

      COMPLEX*16 W(20,NWAVEFUNCS)
      COMMON/W/W

C     ----------
C     BEGIN CODE
C     ----------

      IF (CHECKPHASE.OR.(.NOT.HELDOUBLECHECKED).OR.GOODAMP(LOOPNUM)
     $ ) THEN
        WE(1)=W1
        WE(2)=W2
        WE(3)=W3
        WE(4)=W4
        M2L(1)=M3**2
        M2L(2)=M1**2
        M2L(3)=M2**2
        PAIRING(1)=P1
        PAIRING(2)=P2
        PAIRING(3)=P3
        R=RANK
        ID=LID
        DO I=0,3
          TEMP=1
          DO J=1,NLOOPLINE
            PL(I,J)=0.D0
            DO K=TEMP,(TEMP+PAIRING(J)-1)
              PL(I,J)=PL(I,J)-DBLE(W(1+I,WE(K)))
            ENDDO
            TEMP=TEMP+PAIRING(J)
          ENDDO
        ENDDO
        CALL CTLOOP(NLOOPLINE,PL,M2L,RANK,RES,STABLE)
      ELSE
        RES(1)=0.0D0
        RES(2)=0.0D0
        RES(3)=0.0D0
        STABLE=.TRUE.
      ENDIF
      END

      SUBROUTINE LOOP_3( LID, W1, W2, W3, M1, M2, M3,  RANK, RES
     $ , STABLE, LOOPNUM)

      INTEGER    NEXTERNAL
      PARAMETER (NEXTERNAL=5)
      INTEGER    NLOOPLINE
      PARAMETER (NLOOPLINE=3)
      INTEGER    NWAVEFUNCS
      PARAMETER (NWAVEFUNCS=28)
      INTEGER    NLOOPGROUPS
      PARAMETER (NLOOPGROUPS=77)
      INTEGER    NCOMB
      PARAMETER (NCOMB=48)
C
C     ARGUMENTS
C
      INTEGER W1, W2, W3
      COMPLEX*16 M1, M2, M3

      COMPLEX*16 RES(3)
      INTEGER LID, RANK, LSYMFACT
      INTEGER LOOPNUM
      LOGICAL STABLE
C
C     LOCAL VARIABLES
C
      REAL*8 PL(0:3,NLOOPLINE)
      COMPLEX*16 M2L(NLOOPLINE)
      INTEGER PAIRING(NLOOPLINE),WE(3)
      INTEGER I, J, K, TEMP
C
C     GLOBAL VARIABLES
C
      INTEGER ID,R
      COMMON/LOOP/ID,R

      LOGICAL CHECKPHASE, HELDOUBLECHECKED
      COMMON/INIT/CHECKPHASE, HELDOUBLECHECKED

      INTEGER HELOFFSET
      INTEGER GOODHEL(NCOMB)
      LOGICAL GOODAMP(NLOOPGROUPS)
      COMMON/FILTERS/GOODAMP,GOODHEL,HELOFFSET

      COMPLEX*16 W(20,NWAVEFUNCS)
      COMMON/W/W

C     ----------
C     BEGIN CODE
C     ----------

      IF (CHECKPHASE.OR.(.NOT.HELDOUBLECHECKED).OR.GOODAMP(LOOPNUM)
     $ ) THEN
        WE(1)=W1
        WE(2)=W2
        WE(3)=W3
        M2L(1)=M3**2
        M2L(2)=M1**2
        M2L(3)=M2**2
        DO I=1,NLOOPLINE
          PAIRING(I)=1
        ENDDO

        R=RANK
        ID=LID
        DO I=0,3
          TEMP=1
          DO J=1,NLOOPLINE
            PL(I,J)=0.D0
            DO K=TEMP,(TEMP+PAIRING(J)-1)
              PL(I,J)=PL(I,J)-DBLE(W(1+I,WE(K)))
            ENDDO
            TEMP=TEMP+PAIRING(J)
          ENDDO
        ENDDO
        CALL CTLOOP(NLOOPLINE,PL,M2L,RANK,RES,STABLE)
      ELSE
        RES(1)=0.0D0
        RES(2)=0.0D0
        RES(3)=0.0D0
        STABLE=.TRUE.
      ENDIF
      END


text:
*dk,rankvolume
      subroutine rankvolume(imsgin,xmsgin,cmsgin,msgtype,nwds,
     &   ierror)
C
C
C #####################################################################
C
C     PURPOSE -
C
C        RANKVOLUME prints data on the smallest volume elements
C        (smallest 100 by default).
C
C     INPUT ARGUMENTS -
C
C         imsgin()  - Integer array of command input tokens
C         xmsgin()  - Real array of command input tokens
C         cmsgin()  - Character array of command input tokens
C         msgtype() - Integer array of command input token types
C         nwds      - Number of command input tokens
C
C     OUTPUT ARGUMENTS -
C
C         ierror - Error Return Code (==0 ==> OK, <>0 ==> Error)
C
C     CHANGE HISTORY -
C
C        $Log: rankvolume.f,v $
C        Revision 2.00  2007/11/09 20:03:59  spchu
C        Import to CVS
C
CPVCS
CPVCS       Rev 1.1   30 Jan 2001 08:12:42   dcg
CPVCS    change format to allow for larger element numbers
CPVCS
CPVCS       Rev 1.0   Tue Dec 02 19:28:58 1997   kuprat
CPVCS    Initial revision.
CPVCS
CPVCS       Rev 1.0   Tue Jun 03 12:44:34 1997   kuprat
CPVCS    Initial revision.
C
C #####################################################################
c
c   RANKVOLUME prints out the lowest volume elements from a mesh, ranked
c   in increasing order.  The default is to print out the 100 lowest
c   volume elements, but this number can be changed by specifying it as an
c   optional second argument to the command.
c   Also printed are the contents of the ITET and JTET arrays
c   corresponding to these elements.
c
c   FORMAT:
c
c      RANKVOLUME / number to rank
c
C #####################################################################

      implicit none

      include 'local_element.h'
C
      integer lenptr
      parameter (lenptr=1000000)

      integer nwds, imsgin(nwds), msgtype(nwds)
      real*8 xmsgin(nwds)
      character*(*) cmsgin(nwds)
C
      integer ierror,ierrw

      character*132 logmess
      pointer (ipitet,itet)
      pointer (ipitetoff,itetoff)
      pointer (ipjtet,jtet)
      pointer (ipjtetoff,jtetoff)
      pointer (ipitettyp,itettyp)
      pointer (ipxic,xic)
      pointer (ipyic,yic)
      pointer (ipzic,zic)
      integer itet(lenptr),itetoff(lenptr),jtet(lenptr),
     &   jtetoff(lenptr),itettyp(lenptr)
      real*8 xic(lenptr),yic(lenptr),zic(lenptr)

      pointer (ipvol,vol),(ipieltlist,ieltlist)
      real*8 vol(lenptr)
      integer ieltlist(lenptr)

      real*8 xicvol(maxnen),yicvol(maxnen),zicvol(maxnen)
C
      character*32 cmo,isubname
      integer length,icmotype,icscode,nnodes,i,j,nelements,mbndry,
     &   nrank,ierr,ielt,nbf,nif
      real*8 ascend

      isubname = 'rankvolume'
C
      ierror=0
C
C  Check that user has specified a valid mesh object.

      call cmo_get_name(cmo,ierror)
      if(ierror.ne.0) then
         write(logmess,'(a)')
     *      'RANKVOLUME: ',cmo,' not a valid mesh object'
         call writloga('default',0,logmess,0,ierrw)
         goto 9999
      endif

         call cmo_get_info('nnodes',cmo,
     *      nnodes,length,icmotype,icscode)
         call cmo_get_info('nelements',cmo,
     *      nelements,length,icmotype,icscode)
         call cmo_get_info('itet',cmo,
     *      ipitet,length,icmotype,icscode)
         call cmo_get_info('itetoff',cmo,
     *      ipitetoff,length,icmotype,icscode)
         call cmo_get_info('jtet',cmo,
     *      ipjtet,length,icmotype,icscode)
         call cmo_get_info('jtetoff',cmo,
     *      ipjtetoff,length,icmotype,icscode)
         call cmo_get_info('xic',cmo,
     *      ipxic,length,icmotype,icscode)
         call cmo_get_info('yic',cmo,
     *      ipyic,length,icmotype,icscode)
         call cmo_get_info('zic',cmo,
     *      ipzic,length,icmotype,icscode)
         call cmo_get_info('itettyp',cmo,
     *      ipitettyp,length,icmotype,icscode)
         call cmo_get_info('mbndry',cmo,
     *      mbndry,length,icmotype,icscode)

      if (nwds.ge.2) then
         call test_argument_type(1,1,2,imsgin,xmsgin,cmsgin,msgtype,
     *      nwds)
         nrank=min(nelements,imsgin(2))
      else
         nrank=min(nelements,100)
      endif
      call mmgetblk('ieltlist',isubname,ipieltlist,nelements,1,icscode)
      call mmgetblk('vol',isubname,ipvol,nelements,2,icscode)
      do i=1,nelements
         do j=1,nelmnen(itettyp(i))
            xicvol(j)=xic(itet(itetoff(i)+j))
            yicvol(j)=yic(itet(itetoff(i)+j))
            zicvol(j)=zic(itet(itetoff(i)+j))
         enddo
         call volume_element(itettyp(i),xicvol,yicvol,zicvol,vol(i))
         ieltlist(i)=i
      enddo

      ascend=1.
      call hpsort1(nelements,vol,ascend,ieltlist)

      write(logmess,9000)
 9000 format('elt. no.       volume     ',
     &   '#ext.bound.faces #int.bound.faces')
      call writloga('default',0,logmess,0,ierr)

      do i=1,nrank
         ielt=ieltlist(i)
         nbf=0
         nif=0
         do j=1,nelmnef(itettyp(ielt))
            if (jtet(jtetoff(ielt)+j).eq.mbndry) then
               nbf=nbf+1
            elseif (jtet(jtetoff(ielt)+j).gt.mbndry) then
               nif=nif+1
            endif
         enddo
         write(logmess,9010) ielt,vol(ielt),nbf,nif
 9010    format(i10,2x,e12.6,2x,i5,6x,i5)
         call writloga('default',0,logmess,0,ierr)
      enddo

 9999 continue
      call mmrelprt(isubname,icscode)

      return
      end

max_forks_repo_name: daniellivingston/LaGriT


text:
template<typename T>
chan<T> &chan<T>::operator>>(T &val) {
    this->d_impl->receive(&val);
    return *this;
}

template<typename T>
chan<T> &chan<T>::operator>>(T &&) {
    this->d_impl->receive();
    return *this;
}

max_stars_repo_name: grebnetiew/cpp-chan

max_issues_repo_name: grebnetiew/cpp-chan

max_forks_repo_name: grebnetiew/cpp-chan


text:
Alhambra Drive runs south from East Covell Boulevard before turning east to Mace Boulevard, making a mirror of the Mace curve. Both Loyola Drive and 5th Street tee in near Fred T. Korematsu Elementary School Korematsu Elementary and Mace Ranch Park.

Traveling along Alhambra you might see:

 Alhambra Apartments
 Seville at Mace Ranch Apartments

http://maps.google.com/maps?qmap+davis+ca&ll38.558855,121.695814&spn0.00641,0.021973&clientubuntu&channelfs&oeutf8&glus&z16&layerc&cbll38.558849,121.69582&panoidUQLTTwMlHMQx725yDJDlCg&cbp12,264.24,,0,5.27 Alhambra Drive Google Streetview

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      MODULE cnvg_I
      INTERFACE
!...Generated by Pacific-Sierra Research 77to90  4.4G  21:13:27  03/10/06
      SUBROUTINE cnvg (PNEW, P, P1, NITER, PL)
      USE vast_kind_param,ONLY: DOUBLE
      REAL(DOUBLE), DIMENSION(1234), INTENT(INOUT) :: PNEW
      REAL(DOUBLE), DIMENSION(1234), INTENT(INOUT) :: P
      REAL(DOUBLE), DIMENSION(1234), INTENT(INOUT) :: P1
      INTEGER, INTENT(IN) :: NITER
      REAL(DOUBLE), INTENT(OUT) :: PL
      END SUBROUTINE
      END INTERFACE
      END MODULE

max_stars_repo_name: openmopac/MOPAC-archive

max_issues_repo_name: openmopac/MOPAC-archive

max_forks_repo_name: openmopac/MOPAC-archive


text:
      SUBROUTINE G1CONC(XPOS,YPOS,THETA,ELONG)
C
C          ------------------------------------------------
C          ROUTINE NO. (1022)   VERSION (A7.7)    11:FEB:85
C          ------------------------------------------------
C
C          THIS GENERATES ELLIPTICAL ARCS CENTRED
C          ABOUT THE CURRENT PLOTTING POSITION.
C
C          THE ARGUMENTS ARE:
C
C          <XPOS>   IS THE START-POSITION X-COORDINATE,
C          <YPOS>   IS THE START-POSITION Y-COORDINATE
C                   (BOTH RELATIVE TO THE CURRENT PLOTTING POSITION),
C          <THETA>  IS THE ANGLE SUBTENDED BY THE ARC (IN RADIANS),
C          <ELONG>  IS THE RATIO OF MAJOR:MINOR AXES LENGTHS.
C
C
      COMMON /T1CPOS/ XCHAR,YCHAR
      COMMON /T1DRES/ DRESX,DRESY
      COMMON /T1PPAT/ XPLOTA,YPLOTA
      COMMON /T1PPBT/ XPLOTB,YPLOTB
      COMMON /T3CONS/ PI
C
C
C          FIRSTLY THE CURRENT POSITIONS ARE STORED.
C
      XSAVEB= XPLOTB
      YSAVEB= YPLOTB
      XSAVEA= XPLOTA
      YSAVEA= YPLOTA
      XSAVEC= XCHAR
      YSAVEC= YCHAR
C
C          THE ELLIPSE CONSTANTS ARE THEN CALCULATED,
C          AND THE REQUIRED PORTION OF THE CURVE IS DRAWN.
C
      DIST= 13.0*DRESX*SQRT(XPOS*XPOS+YPOS*YPOS*ELONG*ELONG)
      IF (DIST.LE.0.0) RETURN
C
      NPTS= MAX0(8,INT(SQRT(DIST)))
      NPTS= MIN0(500,NPTS)*(ABS(THETA)/(2.0*PI))
      IF (NPTS.LE.0) NPTS= 1
C
      DANGL= THETA/NPTS
      COSDA= COS(DANGL)
      SINDA= SIN(DANGL)
      RATIO= ABS(ELONG)
      COEFXX=  COSDA
      COEFXY= -SINDA*RATIO
      COEFYX=  SINDA/(COSDA*RATIO)
      COEFYY=  1.0/COSDA
C
      XNEW= XPOS
      YNEW= YPOS
      XARG= XNEW+XSAVEB
      YARG= YNEW+YSAVEB
      CALL G1LINE(XARG,YARG,.FALSE.)
      DO 100 IPT= 1,NPTS
        XNEW= COEFXX*XNEW+COEFXY*YNEW
        YNEW= COEFYX*XNEW+COEFYY*YNEW
        XARG= XNEW+XSAVEB
        YARG= YNEW+YSAVEB
        CALL G1LINE(XARG,YARG,.TRUE.)
  100 CONTINUE
C
C          LASTLY, THE CURRENT POSITIONS ARE RESTORED.
C
      XPLOTB= XSAVEB
      YPLOTB= YSAVEB
      XPLOTA= XSAVEA
      YPLOTA= YSAVEA
      XCHAR=  XSAVEC
      YCHAR=  YSAVEC
C
      RETURN
      END

max_stars_repo_name: ZedThree/GHOST

max_issues_repo_name: ZedThree/GHOST

max_forks_repo_name: ZedThree/GHOST


text:
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!   EVB-QMDFF - RPMD molecular dynamics and rate constant calculations on
!               black-box generated potential energy surfaces
!
!   Copyright (c) 2021 by Julien Steffen (steffen@pctc.uni-kiel.de)
!                         Stefan Grimme (grimme@thch.uni-bonn.de) (QMDFF code)
!
!   Permission is hereby granted, free of charge, to any person obtaining a
!   copy of this software and associated documentation files (the "Software"),
!   to deal in the Software without restriction, including without limitation
!   the rights to use, copy, modify, merge, publish, distribute, sublicense,
!   and/or sell copies of the Software, and to permit persons to whom the
!   Software is furnished to do so, subject to the following conditions:
!
!   The above copyright notice and this permission notice shall be included in
!   all copies or substantial portions of the Software.
!
!   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
!   IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
!   FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
!   THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
!   LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
!   FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
!   DEALINGS IN THE SOFTWARE.
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!
!     subroutine "fatal" terminates execution due to a user request, a severe
!     error or some other nonstandard condition
!
subroutine fatal
implicit none
write(*,*)
write(*,*) "An ERROR occured! EVB-QMDFF will be terminated."
write(*,*) "If the error message remains ask J.Steffen (steffen@pctc.uni-kiel.de) for further maintenance."
write(*,*)
stop

end subroutine fatal

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_name: chrinide/EVB-QMDFF

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_name: chrinide/EVB-QMDFF

max_forks_repo_licenses: MIT Unlicense