subroutine compute_Mder_legacy(s,n,der,I,J,F)
  real*8, intent(in) :: s
  integer*8, intent(in) :: n
  integer*8, intent(in) :: der
  integer*8, intent(out), dimension(3*n):: I
  integer*8, intent(out), dimension(3*n):: J
  real*8, intent(out), dimension(3*n):: F
  integer*8 :: p
  real*8 :: factor
  if (der==0) then
     factor=1;
  else
     factor=0;
  end if
  do p = 1, n
     I(p) = p
     J(p) = p
     F(p) = 2.0*factor;
  end do
  do p = 1, n-1
     I(n+p) = p
     J(n+p) = p+1
     F(n+p) = -1.0*factor;
     I(2*n+p) = p+1
     J(2*n+p) = p
     F(2*n+p) = -1.0*factor;
  end do
  I(2*n)=n;
  J(2*n)=1;
  if (der == 0) then
     F(2*n)=s*s*s;
  else if (der == 1) then
     F(2*n)=3*s*s;
  else if (der == 2) then
     F(2*n)=3*2*s;
  else if (der == 3) then
     F(2*n)=3*2;
  else
     F(2*n)=0;
  end if
  I(3*n)=1;
  J(3*n)=n;
  F(3*n)=-exp(s);
end subroutine compute_Mder_legacy


text:
      MODULE rotmol_I
      INTERFACE
!...Generated by Pacific-Sierra Research 77to90  4.4G  08:15:45  03/09/06
      SUBROUTINE rotmol (NUMAT, COORD, SINA, COSA, I, J, R)
      USE vast_kind_param,ONLY: DOUBLE
      INTEGER :: NUMAT
      REAL(DOUBLE), DIMENSION(3,NUMAT) :: COORD
      REAL(DOUBLE), INTENT(IN) :: SINA
      REAL(DOUBLE), INTENT(IN) :: COSA
      INTEGER, INTENT(IN) :: I
      INTEGER, INTENT(IN) :: J
      REAL(DOUBLE), DIMENSION(3,3), INTENT(INOUT) :: R
      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:
module mod_mom

 use mod_param

 implicit none

 private
 public momxad,momyad,momzad,momsad, momxp,momyp,momzp, mom_main

 contains

   ! #0
   subroutine mom_main(RHSu,RHSv,RHSw)

     use mod_common
     use mod_gfm

     real, dimension(1:imax,1:jmax,1:kmax) :: RCu,RCv,RCw
     real, dimension(1:imax,1:jmax,1:kmax) :: RDu,RDv,RDw
     real, dimension(1:imax,1:jmax,1:kmax) :: RGu,RGv,RGw

     real, dimension(1:imax,1:jmax,1:kmax), intent(out) :: RHSu,RHSv,RHSw

     !-- Convection --!

     select case (NS_conv)
     case('CNTR2')
        call CNTR2_conv(RCu,RCv,RCw, unew,vnew,wnew)

     case('WENO5')
        call WENO5_conv('x-dir',RCu)
        call WENO5_conv('y-dir',RCv)
        call WENO5_conv('z-dir',RCw)
     end select

     !-- Diffusion --!

     select case (NS_diff_discretization)
     case('del')
        call delta_diff(RDu,RDv,RDw, unew,vnew,wnew, miu, rho_u,rho_v,rho_w)

     case('GFM')
        call gfm4diff(RDu,RDv,RDw, unew,vnew,wnew, rho_u,rho_v,rho_w)
     end select

     !-- Gravity --!

     call NS_grav(RGu,RGv,RGw)

     !-- Output --!

     RHSu = RCu + RDu*visc + RGu
     RHSv = RCv + RDv*visc + RGv
     RHSw = RCw + RDw*visc + RGw

     return
   end subroutine mom_main

   ! #1
   subroutine CNTR2_conv(RCu,RCv,RCw, u,v,w)

     integer :: i,j,k, ip,jp,kp, im,jm,km

     real :: uuip,uuim, uvjp,uvjm, uwkp,uwkm
     real :: uvip,uvim, vvjp,vvjm, wvkp,wvkm
     real :: uwip,uwim, vwjp,vwjm, wwkp,wwkm

     real, dimension(0:,0:,0:), intent(in) :: u,v,w
     real, dimension(1:imax,1:jmax,1:kmax), intent(out) :: RCu,RCv,RCw

     do k = 1,kmax
        do j = 1,jmax
           do i = 1,imax

              ip = i + 1
              jp = j + 1
              kp = k + 1
              im = i - 1
              jm = j - 1
              km = k - 1

              !-- x direction

              uuip  = 0.25 *( u(ip,j ,k ) +u(i,j,k) )*( u(ip,j ,k ) +u(i,j ,k ) )
              uuim  = 0.25 *( u(im,j ,k ) +u(i,j,k) )*( u(im,j ,k ) +u(i,j ,k ) )
              uvjp  = 0.25 *( u(i ,jp,k ) +u(i,j,k) )*( v(ip,j ,k ) +v(i,j ,k ) )
              uvjm  = 0.25 *( u(i ,jm,k ) +u(i,j,k) )*( v(ip,jm,k ) +v(i,jm,k ) )
              uwkp  = 0.25 *( u(i ,j ,kp) +u(i,j,k) )*( w(ip,j ,k ) +w(i,j ,k ) )
              uwkm  = 0.25 *( u(i ,j ,km) +u(i,j,k) )*( w(ip,j ,km) +w(i,j ,km) )

              RCu(i,j,k) = (-uuip +uuim)/dx + (-uvjp +uvjm)/dy + (-uwkp +uwkm)/dz

              !-- y direction

              uvip  = 0.25 *(u(i ,j,k ) +u(i ,jp,k ) ) *(v(i,j,k ) +v(ip,j ,k) )
              uvim  = 0.25 *(u(im,j,k ) +u(im,jp,k ) ) *(v(i,j,k ) +v(im,j ,k) )
              vvjp  = 0.25 *(v(i ,j,k ) +v(i ,jp,k ) ) *(v(i,j,k ) +v(i ,jp,k) )
              vvjm  = 0.25 *(v(i ,j,k ) +v(i ,jm,k ) ) *(v(i,j,k ) +v(i ,jm,k) )
              wvkp  = 0.25 *(w(i ,j,k ) +w(i ,jp,k ) ) *(v(i,j,kp) +v(i ,j ,k) )
              wvkm  = 0.25 *(w(i ,j,km) +w(i ,jp,km) ) *(v(i,j,km) +v(i ,j ,k) )

              RCv(i,j,k) = (-uvip +uvim)/dx + (-vvjp +vvjm)/dy + (-wvkp +wvkm)/dz

              !-- z direction

              uwip  = 0.25 *(w(i,j,k) +w(ip,j ,k ) ) *(u(i ,j ,k) +u(i ,j ,kp) )
              uwim  = 0.25 *(w(i,j,k) +w(im,j ,k ) ) *(u(im,j ,k) +u(im,j ,kp) )
              vwjp  = 0.25 *(w(i,j,k) +w(i ,jp,k ) ) *(v(i ,j ,k) +v(i ,j ,kp) )
              vwjm  = 0.25 *(w(i,j,k) +w(i ,jm,k ) ) *(v(i ,jm,k) +v(i ,jm,kp) )
              wwkp  = 0.25 *(w(i,j,k) +w(i ,j ,kp) ) *(w(i ,j ,k) +w(i ,j ,kp ))
              wwkm  = 0.25 *(w(i,j,k) +w(i ,j ,km) ) *(w(i ,j ,k) +w(i ,j ,km ))

              RCw(i,j,k) = (-uwip +uwim)/dx + (-vwjp +vwjm)/dy + (-wwkp +wwkm)/dz

           enddo
        enddo
     enddo

     if (TwoD) RCu = 0.

     return
   end subroutine CNTR2_conv

   ! #2
   subroutine delta_diff(RDu,RDv,RDw, u,v,w, miu, rho_u,rho_v,rho_w)

     real, dimension(0:,0:,0:), intent(in) :: u,v,w, miu
     real, dimension(1:imax,1:jmax,1:kmax), intent(in) :: rho_u,rho_v,rho_w
     real, dimension(1:imax,1:jmax,1:kmax), intent(out) :: RDu,RDv,RDw

     integer :: i,j,k, im,ip, jm,jp, km,kp

     real :: dudxp,dudxm,dudyp,dudym,dudzp,dudzm
     real :: dvdxp,dvdxm,dvdyp,dvdym,dvdzp,dvdzm
     real :: dwdxp,dwdxm,dwdyp,dwdym,dwdzp,dwdzm
     real :: ma1,ma2,ma3,ma4,ma5,ma6,ma7,ma8,ma9
     real :: rhou,rhov,rhow
     real :: RD1,RD2,RD3,RD4,RD5

     do k = 1,kmax
        do j = 1,jmax
           do i = 1,imax

              ip = i + 1
              jp = j + 1
              kp = k + 1
              im = i - 1
              jm = j - 1
              km = k - 1

              !-- u-component
              !
              ! first-derivatives at u-grid

              dudxp = (u(ip,j,k)-u(i ,j,k))/dx
              dudxm = (u(i ,j,k)-u(im,j,k))/dx

              dudyp = (u(i,jp,k)-u(i,j ,k))/dy
              dudym = (u(i,j ,k)-u(i,jm,k))/dy

              dvdxp = (v(ip,j ,k)-v(i,j ,k))/dx
              dvdxm = (v(ip,jm,k)-v(i,jm,k))/dx

              dudzp = (u(i,j,kp)-u(i,j,k ))/dz
              dudzm = (u(i,j,k )-u(i,j,km))/dz

              dwdxp = (w(ip,j,k )-w(i,j,k ))/dx
              dwdxm = (w(ip,j,km)-w(i,j,km))/dx

              ! averaging at face centers

              ma1  = 0.25 *( miu(i,j,k)+miu(i,jm,k )+miu(ip,jm,k )+miu(ip,j,k) )   ! miu(i+0.5,j-0.5,k)
              ma2  = 0.25 *( miu(i,j,k)+miu(i,jp,k )+miu(ip,jp,k )+miu(ip,j,k) )   ! miu(i+0.5,j+0.5,k)
              ma3  = 0.25 *( miu(i,j,k)+miu(i,j ,km)+miu(ip,j ,km)+miu(ip,j,k) )   ! miu(i+0.5,j,k-0.5)
              ma4  = 0.25 *( miu(i,j,k)+miu(i,j ,kp)+miu(ip, j,kp)+miu(ip,j,k) )   ! miu(i+0.5,j,k+0.5)

              RD1 = 2.*(dudxp*miu(ip,j,k)-dudxm*miu(i,j,k))/dx  ! d/dx(2*miu*du/dx)
              RD2 = (dudyp*ma2-dudym*ma1)/dy                    ! d/dy(miu*du/dy)
              RD3 = (dvdxp*ma2-dvdxm*ma1)/dy                    ! d/dy(miu*dv/dx)
              RD4 = (dudzp*ma4-dudzm*ma3)/dz                    ! d/dz(miu*du/dz)
              RD5 = (dwdxp*ma4-dwdxm*ma3)/dz                    ! d/dz(miu*dw/dx)

              RDu(i,j,k) = (RD1+RD2+RD3+RD4+RD5)/rho_u(i,j,k)

              !-- v-component
              !
              ! first-derivatives at v-grid

              dvdxp = (v(ip,j,k)-v(i ,j,k))/dx
              dvdxm = (v(i ,j,k)-v(im,j,k))/dx

              dudyp = (u(i ,jp,k)-u(i ,j,k))/dy
              dudym = (u(im,jp,k)-u(im,j,k))/dy

              dvdyp = (v(i,jp,k)-v(i,j ,k))/dy
              dvdym = (v(i,j ,k)-v(i,jm,k))/dy

              dvdzp = (v(i,j,kp)-v(i,j,k ))/dz
              dvdzm = (v(i,j,k )-v(i,j,km))/dz

              dwdyp = (w(i,jp,k )-w(i,j,k ))/dy
              dwdym = (w(i,jp,km)-w(i,j,km))/dy

              ! averaging at face centers

              ma2  = 0.25 *( miu(i,j,k)+miu(i,jp,k )+miu(ip,jp,k )+miu(ip,j ,k) )   ! miu(i+0.5,j+0.5,k)
              ma5  = 0.25 *( miu(i,j,k)+miu(i,jp,k )+miu(im,jp,k )+miu(im,j ,k) )   ! miu(i-0.5,j+0.5,k)
              ma6  = 0.25 *( miu(i,j,k)+miu(i,j ,km)+miu(i ,jp,km)+miu(i ,jp,k) )   ! miu(i,j+0.5,k-0.5)
              ma7  = 0.25 *( miu(i,j,k)+miu(i,j ,kp)+miu(i ,jp,kp)+miu(i ,jp,k) )   ! miu(i,j+0.5,k+0.5)

              RD1 = (dvdxp*ma2-dvdxm*ma5)/dx                       ! d/dx(miu*dv/dx)
              RD2 = (dudyp*ma2-dudym*ma5)/dx                       ! d/dx(miu*du/dy)
              RD3 = 2.*(dvdyp*miu(i,jp,k)-dvdym*miu(i,j,k))/dy     ! d/dy(2*miu*dv/dy)
              RD4 = (dvdzp*ma7-dvdzm*ma6)/dz                       ! d/dz(miu*dv/dz)
              RD5 = (dwdyp*ma7-dwdym*ma6)/dz                       ! d/dz(miu*dw/dy)

              RDv(i,j,k) = (RD1+RD2+RD3+RD4+RD5)/rho_v(i,j,k)

              !-- w-component
              !
              ! first-derivatives at w-grid

              dwdxp = (w(ip,j,k)-w(i ,j,k))/dx
              dwdxm = (w(i ,j,k)-w(im,j,k))/dx

              dudzp = (u(i ,j,kp)-u(i ,j,k))/dz
              dudzm = (u(im,j,kp)-u(im,j,k))/dz

              dwdyp = (w(i,jp,k)-w(i,j ,k))/dy
              dwdym = (w(i,j ,k)-w(i,jm,k))/dy

              dvdzp = (v(i,j ,kp)-v(i,j ,k))/dz
              dvdzm = (v(i,jm,kp)-v(i,jm,k))/dz

              dwdzp = (w(i,j,kp)-w(i,j,k ))/dz
              dwdzm = (w(i,j,k )-w(i,j,km))/dz

              ! averaging at cell faces

              ma4  = 0.25 *(miu(i,j,k)+miu(i,j,kp)+miu(ip,j ,kp)+miu(ip,j ,k) )   ! miu(i+0.5,j,k+0.5)
              ma8  = 0.25 *(miu(i,j,k)+miu(i,j,kp)+miu(im,j ,kp)+miu(im,j ,k) )   ! miu(i-0.5,j,k+0.5)
              ma7  = 0.25 *(miu(i,j,k)+miu(i,j,kp)+miu(i ,jp,kp)+miu(i ,jp,k) )   ! miu(i,j+0.5,k+0.5)
              ma9  = 0.25 *(miu(i,j,k)+miu(i,j,kp)+miu(i ,jm,kp)+miu(i ,jm,k) )   ! miu(i,j-0.5,k+0.5)

              RD1 = (dwdxp*ma4-dwdxm*ma8)/dx                    ! d/dx(miu*dw/dx)
              RD2 = (dudzp*ma4-dudzm*ma8)/dx                    ! d/dx(miu*du/dx)
              RD3 = (dwdyp*ma7-dwdym*ma9)/dy                    ! d/dy(miu*dw/dy)
              RD4 = (dvdzp*ma7-dvdzm*ma9)/dy                    ! d/dy(miu*dv/dy)
              RD5 = 2.*(dwdzp*miu(i,j,kp)-dwdzm*miu(i,j,k))/dz  ! d/dz(2*miu*dw/dz)

              RDw(i,j,k) = (RD1+RD2+RD3+RD4+RD5)/rho_w(i,j,k)

           enddo
        enddo
     enddo

     if (TwoD) RDu = 0.

     return
   end subroutine delta_diff

   ! #3
   subroutine NS_grav(RGu,RGv,RGw)

     use mod_common

     integer :: i,j,k

     real, dimension(1:imax,1:jmax,1:kmax), intent(out) :: RGu,RGv,RGw

     do k = 1,kmax
        do j = 1,jmax
           do i = 1,imax

              RGu(i,j,k) = iFr*(1.-rho_avr/rho_u(i,j,k))*g_unit(1)
              RGv(i,j,k) = iFr*(1.-rho_avr/rho_v(i,j,k))*g_unit(2)
              RGw(i,j,k) = iFr*(1.-rho_avr/rho_w(i,j,k))*g_unit(3)

           enddo
        enddo
     enddo

     if (TwoD) RGu = 0.

     return
   end subroutine NS_grav

   ! #1+
   subroutine WENO5_conv(dir,Lconv)
   !                                              !
   ! Compute the NS convective terms using WENO5. !
   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!

   use mod_common
   use mod_bound

   character(len=5), intent(in) :: dir

   real, dimension(1:imax,1:jmax,1:kmax), intent(out) :: Lconv
   real, dimension(1:imax,1:jmax,1:kmax) :: wind_u, wind_v, wind_w
   real, dimension(1:imax,1:jmax,1:kmax) :: convx, convy, convz

   !-- Buffer flux matrix

   real, dimension(-2:i1+2,-2:j1+2,-2:k1+2) :: fx,fy,fz, fdir

   !-- Parameters

   real :: eps_weno
   parameter ( eps_weno = 1.e-6 )
   real :: sigma_1, sigma_2, sigma_3
   parameter ( sigma_1 = 0.1 )
   parameter ( sigma_2 = 0.6 )
   parameter ( sigma_3 = 0.3 )
   integer :: i, j, k
   real :: v1, v2, v3, v4, v5
   real :: S1, S2, S3
   real :: a1, a2, a3
   real :: w1, w2, w3
   real :: vel_vertice, f_x,f_y,f_z

    !-- Initialize fluxes --!

    fx = 0.
    fy = 0.
    fz = 0.

    fx(0:i1,0:j1,0:k1) = unew
    fy(0:i1,0:j1,0:k1) = vnew
    fz(0:i1,0:j1,0:k1) = wnew

    call boundfu(fx)
    call boundfv(fy)
    call boundfw(fz)

    select case (dir)

      case('x-dir')

        if (TwoD) then
           Lconv = 0.
           return
        endif

        !-- Direction flux

        fdir = fx

        !-- Face-velocity

        do k = 1,kmax
          do j = 1,jmax
            do i = 1,imax

              !-- U at unew(i,j,k)

              wind_u(i,j,k) = fx(i,j,k)

              !-- V at unew(i,j,k)

              vel_vertice = (fy(i,j  ,k) + fy(i+1,j  ,k) + &
                             fy(i,j-1,k) + fy(i+1,j-1,k) )/4.

              wind_v(i,j,k) = vel_vertice

              !-- W at unew(i,j,k)

              vel_vertice = (fz(i,j,k  ) + fz(i+1,j,k  ) + &
                             fz(i,j,k-1) + fz(i+1,j,k-1) )/4.

              wind_w(i,j,k) = vel_vertice

            enddo
          enddo
        enddo

      case('y-dir')

        !-- Direction flux

        fdir = fy

        !-- Face-velocity

        do k = 1,kmax
          do j = 1,jmax
            do i = 1,imax

              !-- U at vnew(i,j,k)

              vel_vertice = (fx(i  ,j,k) + fx(i  ,j+1,k) + &
                             fx(i-1,j,k) + fx(i-1,j+1,k) )/4.

              wind_u(i,j,k) = vel_vertice

              !-- V at vnew(i,j,k)

              wind_v(i,j,k) = fy(i,j,k)

              !-- W at vnew(i,j,k)

              vel_vertice = (fz(i,j,k  ) + fz(i,j+1,k  ) + &
                             fz(i,j,k-1) + fz(i,j+1,k-1) )/4.

              wind_w(i,j,k) = vel_vertice

            enddo
          enddo
        enddo

      case('z-dir')

        !-- Direction flux

        fdir = fz

        !-- Face-velocity

        do k = 1,kmax
          do j = 1,jmax
            do i = 1,imax

              !-- U at wnew(i,j,k)

              vel_vertice = (fx(i  ,j,k) + fx(i  ,j,k+1) + &
                             fx(i-1,j,k) + fx(i-1,j,k+1) )/4.

              wind_u(i,j,k) = vel_vertice

              !-- V at wnew(i,j,k)

              vel_vertice = (fy(i,j  ,k) + fy(i,j  ,k+1) + &
                             fy(i,j-1,k) + fy(i,j-1,k+1) )/4.

              wind_v(i,j,k) = vel_vertice

              !-- W at wnew(i,j,k)

              wind_w(i,j,k) = fz(i,j,k)

            enddo
          enddo
        enddo

    end select

    do k = 1,kmax
      do j = 1,jmax
        do i = 1,imax

          !-- X-derivative

          if (wind_u(i,j,k) .gt. 0.) then

            v1 = (fdir(i-2,j,k) - fdir(i-3,j,k))/dx
            v2 = (fdir(i-1,j,k) - fdir(i-2,j,k))/dx
            v3 = (fdir(i  ,j,k) - fdir(i-1,j,k))/dx
            v4 = (fdir(i+1,j,k) - fdir(i  ,j,k))/dx
            v5 = (fdir(i+2,j,k) - fdir(i+1,j,k))/dx

            !-- Smoothness indicators

            S1 = (13./12.) * ( v1 - 2.*v2 + v3 )**2 &
                 + (1./4.) * ( v1 - 4.*v2 + 3.*v3)**2

            S2 = (13./12.) * (v2 - 2.*v3 + v4)**2 &
                 + (1./4.) * (v2 - v4)**2

            S3 = (13./12.) * (v3 - 2.*v4 + v5)**2 &
                 + (1./4.) * (3.*v3 - 4.*v4 + v5)**2

            !-- Nonlinear weights

            a1 = sigma_1 / (eps_weno + S1)**2
            a2 = sigma_2 / (eps_weno + S2)**2
            a3 = sigma_3 / (eps_weno + S3)**2

            w1 = a1 / (a1 + a2 + a3)
            w2 = a2 / (a1 + a2 + a3)
            w3 = a3 / (a1 + a2 + a3)

            !-- Minus flux

            f_x = w1 * ( (1./3.)  * v1   &
                       - (7./6.)  * v2   &
                       + (11./6.) * v3 ) &
                + w2 * ( (-1./6.) * v2   &
                       + (5./6.)  * v3   &
                       + (1./3.)  * v4 ) &
                + w3 * ( (1./3.)  * v3   &
                       + (5./6.)  * v4   &
                       - (1./6.)  * v5 )
          else

            v1 = (fdir(i+3,j,k) - fdir(i+2,j,k))/dx
            v2 = (fdir(i+2,j,k) - fdir(i+1,j,k))/dx
            v3 = (fdir(i+1,j,k) - fdir(i  ,j,k))/dx
            v4 = (fdir(i  ,j,k) - fdir(i-1,j,k))/dx
            v5 = (fdir(i-1,j,k) - fdir(i-2,j,k))/dx

            !++ Smoothness indicators

            S1 = (13./12.) * ( v1 - 2.*v2 + v3 )**2 &
                 + (1./4.) * ( v1 - 4.*v2 + 3.*v3)**2

            S2 = (13./12.) * (v2 - 2.*v3 + v4)**2 &
                         + (1./4.) * (v2 - v4)**2

            S3 = (13./12.) * (v3 - 2.*v4 + v5)**2 &
                 + (1./4.) * (3.*v3 - 4.*v4 + v5)**2

            !++ Nonlinear weights

            a1 = sigma_1 / (eps_weno + S1)**2
            a2 = sigma_2 / (eps_weno + S2)**2
            a3 = sigma_3 / (eps_weno + S3)**2

            w1 = a1 / (a1 + a2 + a3)
            w2 = a2 / (a1 + a2 + a3)
            w3 = a3 / (a1 + a2 + a3)

            !++ Plus flux

            f_x = w1 * ( (1./3.)  * v1   &
                       - (7./6.)  * v2   &
                       + (11./6.) * v3 ) &
                + w2 * ( (-1./6.) * v2   &
                       + (5./6.)  * v3   &
                       + (1./3.)  * v4 ) &
                + w3 * ( (1./3.)  * v3   &
                       + (5./6.)  * v4   &
                       - (1./6.)  * v5 )

          endif

          ! Time derivative (x-component)

          convx(i,j,k) = wind_u(i,j,k)*f_x

          !-- Y-derivative

          if (wind_v(i,j,k) .gt. 0.) then

            v1 = (fdir(i,j-2,k) - fdir(i,j-3,k))/dy
            v2 = (fdir(i,j-1,k) - fdir(i,j-2,k))/dy
            v3 = (fdir(i,j  ,k) - fdir(i,j-1,k))/dy
            v4 = (fdir(i,j+1,k) - fdir(i,j  ,k))/dy
            v5 = (fdir(i,j+2,k) - fdir(i,j+1,k))/dy

            !-- Smoothness indicators

            S1 = (13./12.) * ( v1 - 2.*v2 + v3 )**2 &
                 + (1./4.) * ( v1 - 4.*v2 + 3.*v3)**2

            S2 = (13./12.) * (v2 - 2.*v3 + v4)**2 &
                           + (1./4.) * (v2 - v4)**2

            S3 = (13./12.) * (v3 - 2.*v4 + v5)**2 &
                 + (1./4.) * (3.*v3 - 4.*v4 + v5)**2

            !-- Nonlinear weights

            a1 = sigma_1 / (eps_weno + S1)**2
            a2 = sigma_2 / (eps_weno + S2)**2
            a3 = sigma_3 / (eps_weno + S3)**2

            w1 = a1 / (a1 + a2 + a3)
            w2 = a2 / (a1 + a2 + a3)
            w3 = a3 / (a1 + a2 + a3)

            !-- Minus flux

            f_y = w1 * ( (1./3.)  * v1   &
                       - (7./6.)  * v2   &
                       + (11./6.) * v3 ) &
                + w2 * ( (-1./6.) * v2   &
                       + (5./6.)  * v3   &
                       + (1./3.)  * v4 ) &
                + w3 * ( (1./3.)  * v3   &
                       + (5./6.)  * v4   &
                       - (1./6.)  * v5 )

          else

            v1 = (fdir(i,j+3,k) - fdir(i,j+2,k))/dy
            v2 = (fdir(i,j+2,k) - fdir(i,j+1,k))/dy
            v3 = (fdir(i,j+1,k) - fdir(i,j  ,k))/dy
            v4 = (fdir(i,j  ,k) - fdir(i,j-1,k))/dy
            v5 = (fdir(i,j-1,k) - fdir(i,j-2,k))/dy

            !++ Smoothness indicators

            S1 = (13./12.) * ( v1 - 2.*v2 + v3 )**2 &
                 + (1./4.) * ( v1 - 4.*v2 + 3.*v3)**2

            S2 = (13./12.) * (v2 - 2.*v3 + v4)**2 &
                           + (1./4.) * (v2 - v4)**2

            S3 = (13./12.) * (v3 - 2.*v4 + v5)**2 &
                 + (1./4.) * (3.*v3 - 4.*v4 + v5)**2

            !++ Nonlinear weights

            a1 = sigma_1 / (eps_weno + S1)**2
            a2 = sigma_2 / (eps_weno + S2)**2
            a3 = sigma_3 / (eps_weno + S3)**2

            w1 = a1 / (a1 + a2 + a3)
            w2 = a2 / (a1 + a2 + a3)
            w3 = a3 / (a1 + a2 + a3)

            !++ Plus flux

            f_y = w1 * ( (1./3.)  * v1   &
                       - (7./6.)  * v2   &
                       + (11./6.) * v3 ) &
                + w2 * ( (-1./6.) * v2   &
                       + (5./6.)  * v3   &
                       + (1./3.)  * v4 ) &
                + w3 * ( (1./3.)  * v3   &
                       + (5./6.)  * v4   &
                       - (1./6.)  * v5 )

          endif

          ! Time derivative (y-component)

          convy(i,j,k) = wind_v(i,j,k)*f_y

          !-- Z-derivative

          if (wind_w(i,j,k) .gt. 0.) then

            v1 = (fdir(i,j,k-2) - fdir(i,j,k-3))/dz
            v2 = (fdir(i,j,k-1) - fdir(i,j,k-2))/dz
            v3 = (fdir(i,j,k  ) - fdir(i,j,k-1))/dz
            v4 = (fdir(i,j,k+1) - fdir(i,j,k  ))/dz
            v5 = (fdir(i,j,k+2) - fdir(i,j,k+1))/dz

            !-- Smoothness indicators

            S1 = (13./12.) * ( v1 - 2.*v2 + v3 )**2 &
                 + (1./4.) * ( v1 - 4.*v2 + 3.*v3)**2

            S2 = (13./12.) * (v2 - 2.*v3 + v4)**2 &
                 + (1./4.) * (v2 - v4)**2

            S3 = (13./12.) * (v3 - 2.*v4 + v5)**2 &
                 + (1./4.) * (3.*v3 - 4.*v4 + v5)**2

            !-- Nonlinear weights

            a1 = sigma_1 / (eps_weno + S1)**2
            a2 = sigma_2 / (eps_weno + S2)**2
            a3 = sigma_3 / (eps_weno + S3)**2

            w1 = a1 / (a1 + a2 + a3)
            w2 = a2 / (a1 + a2 + a3)
            w3 = a3 / (a1 + a2 + a3)

            !-- Minus flux

            f_z = w1 * ( (1./3.)  * v1   &
                       - (7./6.)  * v2   &
                       + (11./6.) * v3 ) &
                + w2 * ( (-1./6.) * v2   &
                       + (5./6.)  * v3   &
                       + (1./3.)  * v4 ) &
                + w3 * ( (1./3.)  * v3   &
                       + (5./6.)  * v4   &
                       - (1./6.)  * v5 )

          else !++

            v1 = (fdir(i,j,k+3) - fdir(i,j,k+2))/dz
            v2 = (fdir(i,j,k+2) - fdir(i,j,k+1))/dz
            v3 = (fdir(i,j,k+1) - fdir(i,j,k  ))/dz
            v4 = (fdir(i,j,k  ) - fdir(i,j,k-1))/dz
            v5 = (fdir(i,j,k-1) - fdir(i,j,k-2))/dz

            !++ Smoothness indicators

            S1 = (13./12.) * ( v1 - 2.*v2 + v3 )**2 &
                 + (1./4.) * ( v1 - 4.*v2 + 3.*v3)**2

            S2 = (13./12.) * (v2 - 2.*v3 + v4)**2 &
                 + (1./4.) * (v2 - v4)**2

            S3 = (13./12.) * (v3 - 2.*v4 + v5)**2 &
                 + (1./4.) * (3.*v3 - 4.*v4 + v5)**2

            !++ Nonlinear weights

            a1 = sigma_1 / (eps_weno + S1)**2
            a2 = sigma_2 / (eps_weno + S2)**2
            a3 = sigma_3 / (eps_weno + S3)**2

            w1 = a1 / (a1 + a2 + a3)
            w2 = a2 / (a1 + a2 + a3)
            w3 = a3 / (a1 + a2 + a3)

            !++ Plus flux

            f_z = w1 * ( (1./3.)  * v1   &
                       - (7./6.)  * v2   &
                       + (11./6.) * v3 ) &
                + w2 * ( (-1./6.) * v2   &
                       + (5./6.)  * v3   &
                       + (1./3.)  * v4 ) &
                + w3 * ( (1./3.)  * v3   &
                       + (5./6.)  * v4   &
                       - (1./6.)  * v5 )
          endif

          ! Time derivative (z-component)

          convz(i,j,k) = wind_w(i,j,k)*f_z

        enddo
      enddo
    enddo

    !-- Output convective fluxes --!

    Lconv = - convx - convy - convz

    return
   end subroutine WENO5_conv

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!

  !----------------
  !
  !     Subroutines mom?p compute pressure corection terms split from the pressure gradient in ? direction
  !     as used in Dodd & Ferrante (the rest computed with FFT)
  !
  !     INPUT  : Predicted pressure and density at the next time level
  !     OUTPUT : "Acceleration" at an intermediate time level (AB2)
  !
  !     Key internal variables
  !     ----------------------
  !     rho?         : Averaged density at face centers
  !     rho0         : Min density of phase 1 and 2
  !
  !-----

  subroutine momxp(dudt,p,p_j,rho_u)

    real, dimension(0:,0:,0:), intent(in) :: p,p_j
    real, dimension(1:,1:,1:), intent(in) :: rho_u
    real, dimension(1:,1:,1:), intent(out) :: dudt
    integer :: i,j,k
    integer :: ip

    do k=1,kmax
       do j=1,jmax
          do i=1,imax

             ip = i + 1
             dudt(i,j,k) = - (1./rho_u(i,j,k)-1./rho0)*dxi*( p(ip,j,k)-p_j(i,j,k) -p(i,j,k) )

          enddo
       enddo
    enddo

    return
  end subroutine momxp

  subroutine momyp(dvdt,p,p_j,rho_v)

    real, dimension(0:,0:,0:), intent(in) :: p,p_j
    real, dimension(1:,1:,1:), intent(in) :: rho_v
    real, dimension(1:,1:,1:), intent(out) :: dvdt
    integer :: i,j,k
    integer :: jp

    do k=1,kmax
       do j=1,jmax
          do i=1,imax

             jp = j + 1
             dvdt(i,j,k) = - (1./rho_v(i,j,k)-1./rho0)*dyi*( p(i,jp,k)-p_j(i,j,k) -p(i,j,k) )

          enddo
       enddo
    enddo

    return
  end subroutine momyp

  subroutine momzp(dwdt,p,p_j,rho_w)

    real, dimension(0:,0:,0:), intent(in) :: p,p_j
    real, dimension(1:,1:,1:), intent(in) :: rho_w
    real, dimension(1:,1:,1:), intent(out) :: dwdt
    integer :: kp
    integer :: i,j,k

    do k=1,kmax
       do j=1,jmax
          do i=1,imax

             kp = k + 1
             dwdt(i,j,k) = - (1./rho_w(i,j,k)-1./rho0)*dzi*( p(i,j,kp)-p_j(i,j,k) -p(i,j,k) )

          enddo
       enddo
    enddo

    return
  end subroutine momzp

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! Not used below

 !
 ! Subroutines mom?ad compute convection, diffusion, and gravity terms on the RHS of NS in ? direction
 ! using a second-order accurate central difference scheme.
 !
 ! INPUT  : Velocity, viscosity, and density at a known time level
 ! OUTPUT : "Acceleration" at an intermediate time level (AB2) or at the previous time level (restarting)
 !
 ! Key internal variables
 ! ----------------------
 ! ma?          : Averaged viscosity at face centers
 ! rho?         : Averaged density at face centers
 ! RC           : Negative convection terms
 ! sum of RD?   : Diffusion terms
 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

   ! #1

   subroutine momxad(dudt, u,v,w,miu,rho,g_unit)

     real, intent(in) :: g_unit(1)
     real, dimension(0:,0:,0:), intent(in) :: u,v,w,miu,rho
     real, dimension(0:,0:,0:), intent(out) :: dudt

     integer :: im,ip,jm,jp,km,kp,i,j,k

     real :: uuip,uuim,uvjp,uvjm,uwkp,uwkm
     real :: dudxp,dudxm,dudyp,dudym,dvdxp,dvdxm,dudzp,dudzm,dwdxp,dwdxm
     real :: ma1,ma2,ma3,ma4
     real :: rhou  !#
     real :: RD1,RD2,RD3,RD4,RD5

     real, dimension(1:imax,1:jmax,1:kmax) :: RC,RD,RG

     !$omp parallel default(shared) &
     !$omp&private(i,j,k,ip,jp,kp,im,jm,km) &
     !$omp&private(uuip,uuim,uvjp,uvjm,uwkp,uwkm,dudxp,dudxm,dudyp,dudym,dudzp,dudzm)
     !$omp do

     if (TwoD) then
        dudt = 0.
        return
     endif

     !-- Convection --!

     select case (NS_conv)

     case('CNTR2')

        do k = 1,kmax
           do j = 1,jmax
              do i = 1,imax

                 ip = i + 1
                 jp = j + 1
                 kp = k + 1
                 im = i - 1
                 jm = j - 1
                 km = k - 1

                 uuip  = 0.25 *( u(ip,j ,k ) +u(i,j,k) )*( u(ip,j ,k ) +u(i,j ,k ) )
                 uuim  = 0.25 *( u(im,j ,k ) +u(i,j,k) )*( u(im,j ,k ) +u(i,j ,k ) )
                 uvjp  = 0.25 *( u(i ,jp,k ) +u(i,j,k) )*( v(ip,j ,k ) +v(i,j ,k ) )
                 uvjm  = 0.25 *( u(i ,jm,k ) +u(i,j,k) )*( v(ip,jm,k ) +v(i,jm,k ) )
                 uwkp  = 0.25 *( u(i ,j ,kp) +u(i,j,k) )*( w(ip,j ,k ) +w(i,j ,k ) )
                 uwkm  = 0.25 *( u(i ,j ,km) +u(i,j,k) )*( w(ip,j ,km) +w(i,j ,km) )

                 RC(i,j,k) = (-uuip +uuim)/dx + (-uvjp +uvjm)/dy + (-uwkp +uwkm)/dz

              enddo
           enddo
        enddo

     case('WENO5')

        call WENO5_conv('x-dir',RC)

     end select

     !-- Diffusion --!

     do k = 1,kmax
        do j = 1,jmax
           do i = 1,imax

              ip = i + 1
              jp = j + 1
              kp = k + 1
              im = i - 1
              jm = j - 1
              km = k - 1

              ! velocity gradients

              dudxp = (u(ip,j,k)-u(i ,j,k))/dx
              dudxm = (u(i ,j,k)-u(im,j,k))/dx

              dudyp = (u(i,jp,k)-u(i,j ,k))/dy
              dudym = (u(i,j ,k)-u(i,jm,k))/dy

              dvdxp = (v(ip,j ,k)-v(i,j ,k))/dx
              dvdxm = (v(ip,jm,k)-v(i,jm,k))/dx

              dudzp = (u(i,j,kp)-u(i,j,k ))/dz
              dudzm = (u(i,j,k )-u(i,j,km))/dz

              dwdxp = (w(ip,j,k )-w(i,j,k ))/dx
              dwdxm = (w(ip,j,km)-w(i,j,km))/dx

              ! averaging at face centers

              ma1  = 0.25 *( miu(i,j,k)+miu(i,jm,k )+miu(ip,jm,k )+miu(ip,j,k) )   ! miu(i+0.5,j-0.5,k)
              ma2  = 0.25 *( miu(i,j,k)+miu(i,jp,k )+miu(ip,jp,k )+miu(ip,j,k) )   ! miu(i+0.5,j+0.5,k)
              ma3  = 0.25 *( miu(i,j,k)+miu(i,j ,km)+miu(ip,j ,km)+miu(ip,j,k) )   ! miu(i+0.5,j,k-0.5)
              ma4  = 0.25 *( miu(i,j,k)+miu(i,j ,kp)+miu(ip, j,kp)+miu(ip,j,k) )   ! miu(i+0.5,j,k+0.5)

              RD1 = 2.*(dudxp*miu(ip,j,k)-dudxm*miu(i,j,k))/dx  ! d/dx(2*miu*du/dx)
              RD2 = (dudyp*ma2-dudym*ma1)/dy                    ! d/dy(miu*du/dy)
              RD3 = (dvdxp*ma2-dvdxm*ma1)/dy                    ! d/dy(miu*dv/dx)
              RD4 = (dudzp*ma4-dudzm*ma3)/dz                    ! d/dz(miu*du/dz)
              RD5 = (dwdxp*ma4-dwdxm*ma3)/dz                    ! d/dz(miu*dw/dx)

              rhou = (rho(i,j,k)+rho(ip,j,k))/2.

              RD(i,j,k) = (RD1+RD2+RD3+RD4+RD5)/rhou

           enddo
        enddo
     enddo

     !-- Gravity --!

     do k = 1,kmax
        do j = 1,jmax
           do i = 1,imax

              rhou = (rho(i,j,k)+rho(i+1,j,k))/2.

              RG(i,j,k) = iFr*(1.-rho_avr/rhou)*g_unit(1)

           enddo
        enddo
     enddo

     !-- Output  --!

     dudt(1:imax,1:jmax,1:kmax) = RC + RD*visc + RG

     !$omp end parallel

     return
   end subroutine momxad

   ! #2

   subroutine momyad(dvdt, u,v,w,miu,rho,g_unit)

     real, intent(in) :: g_unit(2)
     real, dimension(0:,0:,0:), intent(in) :: u,v,w,miu,rho
     real, dimension(0:,0:,0:), intent(out) :: dvdt

     integer :: im,ip,jm,jp,km,kp,i,j,k

     real :: uvip,uvim,vvjp,vvjm,wvkp,wvkm
     real :: dvdxp,dvdxm,dudyp,dudym,dvdyp,dvdym,dvdzp,dvdzm,dwdyp,dwdym
     real :: ma2,ma5,ma6,ma7
     real :: rhov  !#
     real :: RD1,RD2,RD3,RD4,RD5

     real, dimension(1:imax,1:jmax,1:kmax) :: RC,RD,RG

     !$omp parallel default(shared) &
     !$omp&private(i,j,k,ip,jp,kp,im,jm,km) &
     !$omp&private(uvip,uvim,vvjp,vvjm,wvkp,wvkm,dvdxp,dvdxm,dvdyp,dvdym,dvdzp,dvdzm)
     !$omp do

     !-- Convection --!

     select case (NS_conv)

     case('CNTR2')

        do k = 1,kmax
           do j = 1,jmax
              do i = 1,imax

                 ip = i + 1
                 jp = j + 1
                 kp = k + 1
                 im = i - 1
                 jm = j - 1
                 km = k - 1

                 uvip  = 0.25 *(u(i ,j,k ) +u(i ,jp,k ) ) *(v(i,j,k ) +v(ip,j ,k) )
                 uvim  = 0.25 *(u(im,j,k ) +u(im,jp,k ) ) *(v(i,j,k ) +v(im,j ,k) )
                 vvjp  = 0.25 *(v(i ,j,k ) +v(i ,jp,k ) ) *(v(i,j,k ) +v(i ,jp,k) )
                 vvjm  = 0.25 *(v(i ,j,k ) +v(i ,jm,k ) ) *(v(i,j,k ) +v(i ,jm,k) )
                 wvkp  = 0.25 *(w(i ,j,k ) +w(i ,jp,k ) ) *(v(i,j,kp) +v(i ,j ,k) )
                 wvkm  = 0.25 *(w(i ,j,km) +w(i ,jp,km) ) *(v(i,j,km) +v(i ,j ,k) )

                 RC(i,j,k) = (-uvip +uvim)/dx + (-vvjp +vvjm)/dy + (-wvkp +wvkm)/dz

              enddo
           enddo
        enddo

     case('WENO5')

        call WENO5_conv('y-dir',RC)

     end select

     !-- Diffusion --!

     do k = 1,kmax
        do j = 1,jmax
           do i = 1,imax

              ip = i + 1
              jp = j + 1
              kp = k + 1
              im = i - 1
              jm = j - 1
              km = k - 1

              ! velocity gradients

              dvdxp = (v(ip,j,k)-v(i ,j,k))/dx
              dvdxm = (v(i ,j,k)-v(im,j,k))/dx

              dudyp = (u(i ,jp,k)-u(i ,j,k))/dy
              dudym = (u(im,jp,k)-u(im,j,k))/dy

              dvdyp = (v(i,jp,k)-v(i,j ,k))/dy
              dvdym = (v(i,j ,k)-v(i,jm,k))/dy

              dvdzp = (v(i,j,kp)-v(i,j,k ))/dz
              dvdzm = (v(i,j,k )-v(i,j,km))/dz

              dwdyp = (w(i,jp,k )-w(i,j,k ))/dy
              dwdym = (w(i,jp,km)-w(i,j,km))/dy

              ! averaging at face centers

              ma2  = 0.25 *( miu(i,j,k)+miu(i,jp,k )+miu(ip,jp,k )+miu(ip,j ,k) )   ! miu(i+0.5,j+0.5,k)
              ma5  = 0.25 *( miu(i,j,k)+miu(i,jp,k )+miu(im,jp,k )+miu(im,j ,k) )   ! miu(i-0.5,j+0.5,k)
              ma6  = 0.25 *( miu(i,j,k)+miu(i,j ,km)+miu(i ,jp,km)+miu(i ,jp,k) )   ! miu(i,j+0.5,k-0.5)
              ma7  = 0.25 *( miu(i,j,k)+miu(i,j ,kp)+miu(i ,jp,kp)+miu(i ,jp,k) )   ! miu(i,j+0.5,k+0.5)

              RD1 = (dvdxp*ma2-dvdxm*ma5)/dx                       ! d/dx(miu*dv/dx)
              RD2 = (dudyp*ma2-dudym*ma5)/dx                       ! d/dx(miu*du/dx)
              RD3 = 2.*(dvdyp*miu(i,jp,k)-dvdym*miu(i,j,k))/dy     ! d/dy(2*miu*dv/dy)
              RD4 = (dvdzp*ma7-dvdzm*ma6)/dz                       ! d/dz(miu*dv/dz)
              RD5 = (dwdyp*ma7-dwdym*ma6)/dz                       ! d/dz(miu*dw/dy)

              rhov = (rho(i,j,k)+rho(i,jp,k))/2.

              RD(i,j,k) = (RD1+RD2+RD3+RD4+RD5)/rhov

           enddo
        enddo
     enddo

     !-- Gravity --!

     do k = 1,kmax
        do j = 1,jmax
           do i = 1,imax

              rhov = (rho(i,j,k)+rho(i,j+1,k))/2.

              RG(i,j,k) = iFr*(1.-rho_avr/rhov)*g_unit(2)

           enddo
        enddo
     enddo
     !$omp end parallel

     !-- Output --!

     dvdt(1:imax,1:jmax,1:kmax) = RC + RD*visc + RG

     return
   end subroutine momyad

   ! #3

   subroutine momzad(dwdt, u,v,w,miu,rho,g_unit)

     real, intent(in) :: g_unit(3)
     real, dimension(0:,0:,0:), intent(in) :: u,v,w,miu,rho
     real, dimension(0:,0:,0:), intent(out) :: dwdt

     integer :: im,ip,jm,jp,km,kp,i,j,k

     real :: uwip,uwim,vwjp,vwjm,wwkp,wwkm
     real :: dwdxp,dwdxm,dudzp,dudzm,dwdyp,dwdym,dvdzp,dvdzm,dwdzp,dwdzm
     real :: ma4,ma8,ma7,ma9
     real :: rhow  !#
     real :: RD1,RD2,RD3,RD4,RD5

     real, dimension(1:imax,1:jmax,1:kmax) :: RC,RD,RG

     !$omp parallel default(shared) &
     !$omp&private(i,j,k,ip,jp,kp,im,jm,km) &
     !$omp&private(uwip,uwim,vwjp,vwjm,wwkp,wwkm,dwdxp,dwdxm,dwdyp,dwdym,dwdzp,dwdzm)
     !$omp do

     !-- Convection --!

     select case (NS_conv)

     case('CNTR2')

        do k = 1,kmax
           do j = 1,jmax
              do i = 1,imax

                 ip = i + 1
                 jp = j + 1
                 kp = k + 1
                 im = i - 1
                 jm = j - 1
                 km = k - 1

                 uwip  = 0.25 *(w(i,j,k) +w(ip,j ,k ) ) *(u(i ,j ,k) +u(i ,j ,kp) )
                 uwim  = 0.25 *(w(i,j,k) +w(im,j ,k ) ) *(u(im,j ,k) +u(im,j ,kp) )
                 vwjp  = 0.25 *(w(i,j,k) +w(i ,jp,k ) ) *(v(i ,j ,k) +v(i ,j ,kp) )
                 vwjm  = 0.25 *(w(i,j,k) +w(i ,jm,k ) ) *(v(i ,jm,k) +v(i ,jm,kp) )
                 wwkp  = 0.25 *(w(i,j,k) +w(i ,j ,kp) ) *(w(i ,j ,k) +w(i ,j ,kp ))
                 wwkm  = 0.25 *(w(i,j,k) +w(i ,j ,km) ) *(w(i ,j ,k) +w(i ,j ,km ))

                 RC(i,j,k) = (-uwip +uwim)/dx + (-vwjp +vwjm)/dy + (-wwkp +wwkm)/dz

              enddo
           enddo
        enddo

     case('WENO5')

        call WENO5_conv('z-dir',RC)

     end select

     !-- Diffusion --!

     do k = 1,kmax
        do j = 1,jmax
           do i = 1,imax

              ip = i + 1
              jp = j + 1
              kp = k + 1
              im = i - 1
              jm = j - 1
              km = k - 1

              ! velocity gradients

              dwdxp = (w(ip,j,k)-w(i ,j,k))/dx
              dwdxm = (w(i ,j,k)-w(im,j,k))/dx

              dudzp = (u(i ,j,kp)-u(i ,j,k))/dz
              dudzm = (u(im,j,kp)-u(im,j,k))/dz

              dwdyp = (w(i,jp,k)-w(i,j ,k))/dy
              dwdym = (w(i,j ,k)-w(i,jm,k))/dy

              dvdzp = (v(i,j ,kp)-v(i,j ,k))/dz
              dvdzm = (v(i,jm,kp)-v(i,jm,k))/dz

              dwdzp = (w(i,j,kp)-w(i,j,k ))/dz
              dwdzm = (w(i,j,k )-w(i,j,km))/dz

              ! averaging at cell faces

              ma4  = 0.25 *(miu(i,j,k)+miu(i,j,kp)+miu(ip,j ,kp)+miu(ip,j ,k) )   ! miu(i+0.5,j,k+0.5)
              ma8  = 0.25 *(miu(i,j,k)+miu(i,j,kp)+miu(im,j ,kp)+miu(im,j ,k) )   ! miu(i-0.5,j,k+0.5)
              ma7  = 0.25 *(miu(i,j,k)+miu(i,j,kp)+miu(i ,jp,kp)+miu(i ,jp,k) )   ! miu(i,j+0.5,k+0.5)
              ma9  = 0.25 *(miu(i,j,k)+miu(i,j,kp)+miu(i ,jm,kp)+miu(i ,jm,k) )   ! miu(i,j-0.5,k+0.5)

              RD1 = (dwdxp*ma4-dwdxm*ma8)/dx                    ! d/dx(miu*dw/dx)
              RD2 = (dudzp*ma4-dudzm*ma8)/dx                    ! d/dx(miu*du/dx)
              RD3 = (dwdyp*ma7-dwdym*ma9)/dy                    ! d/dy(miu*dw/dy)
              RD4 = (dvdzp*ma7-dvdzm*ma9)/dy                    ! d/dy(miu*dv/dy)
              RD5 = 2.*(dwdzp*miu(i,j,kp)-dwdzm*miu(i,j,k))/dz  ! d/dz(2*miu*dw/dz)

              rhow = (rho(i,j,k)+rho(i,j,kp))/2.

              RD(i,j,k) = (RD1+RD2+RD3+RD4+RD5)/rhow

           enddo
        enddo
     enddo

     !-- Gravity --!

     do k = 1,kmax
        do j = 1,jmax
           do i = 1,imax

              rhow = (rho(i,j,k)+rho(i,j,k+1))/2.

              RG(i,j,k) = iFr*(1.-rho_avr/rhow)*g_unit(3)

           enddo
        enddo
     enddo
     !$omp end parallel

     !-- Output --!

     dwdt(1:imax,1:jmax,1:kmax) = RC + RD*visc + RG

     return
   end subroutine momzad

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

!# convection of a scalar
subroutine momsad(dydt,u,v,w,y)
implicit none

real, dimension(0:,0:,0:), intent(in) :: u,v,w
real, dimension(-2:,-2:,-2:), intent(in) :: y
real, dimension(-2:,-2:,-2:), intent(out) :: dydt
integer :: im,ip,jm,jp,km,kp,i,j,k
real :: yuip,yuim,yvjp,yvjm,ywkp,ywkm
real :: dydxp,dydxm,dydyp,dydym,dydzp,dydzm

!$omp parallel default(shared) &
!$omp&private(i,j,k,ip,jp,kp,im,jm,km) &
!$omp&private(uwip,uwim,vwjp,vwjm,wwkp,wwkm,dwdxp,dwdxm,dwdyp,dwdym,dwdzp,dwdzm)
!$omp do
 do k=1,kmax
   do j=1,jmax
     do i=1,imax
       ip = i + 1
       jp = j + 1
       kp = k + 1
       im = i - 1
       jm = j - 1
       km = k - 1
       yuip  = 0.5 * (Y(i,j,k)+Y(ip,j,k))*U(i,j,k)
       yuim  = 0.5 * (Y(i,j,k)+Y(im,j,k))*U(im,j,k)
       yvjp  = 0.5 * (Y(i,j,k)+Y(i,jp,k))*V(i,j,k)
       yvjm  = 0.5 * (Y(i,j,k)+Y(i,jm,k))*V(i,jm,k)
       ywkp  = 0.5 * (Y(i,j,k)+Y(i,j,kp))*W(i,j,k)
       ywkm  = 0.5 * (Y(i,j,k)+Y(i,j,km))*W(i,j,km)

       ! Momentum balance

       dydt(i,j,k) = dxi*( -yuip + yuim ) + &
                     dyi*( -yvjp + yvjm ) + &
                     dzi*( -ywkp + ywkm )
     enddo
   enddo
 enddo
!$omp end parallel

return
end subroutine momsad

end module mod_mom

max_stars_repo_name: GeZhouyang/ICLS-release

max_issues_repo_name: GeZhouyang/ICLS-release

max_forks_repo_name: GeZhouyang/ICLS-release


text:
      SUBROUTINE zcofil6 (IFTOLD, IFTNEW, IARRAY, KARRAY, IBUFF, KBUFF,
     * LRCDEL, LRETAGX)
C
C
C     Z-Copy-File.
C     Copies all records in file opened with IFTNEW into
C     file opened with IFTOLD.
C     Records that have been deleted (but have not yet been
C     squeezed out) may be restored by setting LRCDEL to true.
C     New tags will be generated if LRETAG is set to true.
C
C     zcofil6 copies via a brute-force approach;  Thus it
C     will copy all vaild data from a damaged file
C     (restoring it as best as possible).
C
C     Written by Bill Charley at HEC, 1988.
C
      INTEGER IFTNEW(*), IFTOLD(*), IBUFF(KBUFF), INFOCO(30)
      INTEGER IARRAY(KARRAY)
      LOGICAL LRCDEL, LRETAG, LFOUND, LPROT2, LTSCOPY, LRETAGX, LDUP
      CHARACTER CPATH*400
C
C
      INCLUDE 'zdsskz.h'
C
      INCLUDE 'zdsscz.h'
C
      INCLUDE 'zdssiz.h'
C
      INCLUDE 'zdsslz.h'
C
      INCLUDE 'zdssmz.h'
C
      INCLUDE 'zdssnz.h'
C
      COMMON /ZSTATUS/ TOTAL_NUMB,  CURRENT_NUMB,
     *                 INTERRUPT, NERROR, MAXERROR
      INTEGER TOTAL_NUMB,  CURRENT_NUMB,  INTERRUPT
C
C
C
C     The follow is if another thread interrupts us to stop the process
      INTERRUPT = 0
      LRETAG = .FALSE.
      IF (MLEVEL.GE.11) WRITE (MUNIT,20) IFTOLD(KUNIT)
 20   FORMAT (T6,'-----DSS---Debug:  Enter zcofil6;  Unit:',I5)
C
C     Make Sure that we have enough space to work with
      IF ((KARRAY.LT.NBSIZE).OR.(KBUFF.LT.NBSIZE)) GO TO 900
C
C     Be sure the files are same type endian!
      IF (IFTOLD(KSWAP).NE.IFTNEW(KSWAP)) GO TO 930
C
C     Make sure that we are copying the same version
      IF ((IFTOLD(1).NE.6).OR.(IFTNEW(1).NE.6)) GO TO 910
      LPROT2 = LPROTC
      LCOFIL = .TRUE.
C
C     Are we in a read only state?
      IF (IFTNEW(KREADO).EQ.1) GO TO 920
C
C     Print any squeeze status message
C     IF (LSQSTA) CALL CHRWT (MUNIT, CHAR(13)//CHAR(10)//'  0% Complete'
C    * // CHAR(13), 16)
C     IF (LSQSTA) CALL CHRWT (MUNIT, ' 0% Complete' // CHAR(13), 14)    u
      !IF (LSQSTA) CALL CHRWT (MUNIT, '+ 0% Complete' // CHAR(13), 14)   Md
C
C     Get the number of physical records in file (not data recs)
      IADD = IFTOLD(KFSIZE) - 1
      CALL zgetrw6 (IADD, NRECS, IW)
C
      TOTAL_NUMB = IFTOLD(KNRECS)
      IFTNEW(KCOLL) = IFTOLD(KCOLL)
      CURRENT_NUMB = 0
C
      DO 200 IREC=1,NRECS
C
C     If another thead set the interrupt flag, return
      IF (INTERRUPT.NE.0) GO TO 840
C
C     Should the status (% complete) be printed
      IF (LSQSTA) THEN
      N = MOD (IREC,10)
      IF (N.EQ.0) THEN
      IP = ((IREC-1)*100) / NRECS
      WRITE (CPATH, 25) IP
C25   FORMAT (I3,'% Complete')                                          u
 25   FORMAT ('+',I3,'% Complete')                                      Md
      !CALL CHRWT (MUNIT, CPATH(1:14)//CHAR(13), 15)
      ENDIF
      ENDIF
C
C     Read physical record IREC from file
      CALL zgetad6 (IADD, IREC, 1)
      CALL zgtrec6 (IFTOLD, IARRAY, NBSIZE, IADD, .FALSE.)
C
C     Did we hit end of file?
      IF (IADD.LT.0) THEN
      WRITE (MUNIT, 28) IREC
 28   FORMAT (//, ' **** Caution:  Reached End of File at Physical',
     * ' Record:',I8,/,' **** Probable Incomplete File Copy',/)
      GO TO 800
      ENDIF
C
C     Search through this record, looking for pathname flags
C
      DO 200 IWRD=1,NBSIZE
C
      IF (IARRAY(IWRD).EQ.NPFLAG) THEN
C
C     Found a flag - Get the first three words of the information
C     block to see if this indeed is the start of a data record
      CALL zgetad6 (IADD, IREC, IWRD)
      CALL zgtrec6 ( IFTOLD, INFO, 3, IADD, .FALSE.)
      IF (IFTOLD(KSTAT).NE.0) GO TO 940
C
      IF (MLEVEL.GE.10) WRITE (MUNIT,30)IADD, INFO(KISTAT), INFO(KINPAT)
 30   FORMAT (T5,'Found a INFO Flag at address',I13,/,
     * T5,'Status:',I5,';  Pathname Length:',I5)
      IPSTAT = INFO(KISTAT)
      IF (IPSTAT.GT.10) IPSTAT = IPSTAT - 10
C     Is this an undelete of a renamed or deleted record?
      IF ((IPSTAT.EQ.2).AND.(LRCDEL)) INFO(KISTAT) = INFO(KISTAT) - 1
C     Check for a valid status flag (1)
      IF (IPSTAT.EQ.1) THEN
C     Check for a valid pathname length
      NPATH = INFO(KINPAT)
      IF ((NPATH.GT.0).AND.(NPATH.LE.392)) THEN
C
C     Passed, therefore a valid record - Get the full information block
      NPPWRD = (NPATH - 1)/NCPW + 1
      NPMWRD = (NPATH - 1)/NCMW + 1
      CALL zgtrec6 ( IFTOLD, INFO, NINFO+NPPWRD, IADD, .FALSE.)
      IF (IFTOLD(KSTAT).NE.0) GO TO 940
C
C     Get the pathname
      CPATH = ' '
      CALL HOLCHR (INFO(KIPATH), 1, NPATH, CPATH, 1)
      CURRENT_NUMB = CURRENT_NUMB + 1
C
C
C     If this data is time-series, and we need to use the compression
C     method of the new file, or we are copying 15 or 30 min data,
C     call zrrtsx6, then zsrtsx6 via zcorec6
      LTSCOPY = .FALSE.
      IF (INFO(NPPWRD+KITYPE).EQ.100) THEN
      IF ((IFTOLD(KVERNO).LT.67300).OR.(IFTNEW(KVERNO).LT.67300)) THEN
      IF (INDEX(CPATH,'15MIN').GT.0) LTSCOPY = .TRUE.
      IF (INDEX(CPATH,'20MIN').GT.0) LTSCOPY = .TRUE.
      IF (INDEX(CPATH,'30MIN').GT.0) LTSCOPY = .TRUE.
      ENDIF
      IF ((LTSCMP).AND.(INFO(NPPWRD+KIQUAL).EQ.0)) THEN
      IF  ((IFTNEW(KCOMPN).GT.0).OR.(INFO(NPPWRD+KICOMP).GT.0)) THEN
      LTSCOPY = .TRUE.
      ENDIF
      ENDIF
      ENDIF
C
C     If converting from a previous version with a possible double
C     swap miss-match, use zcorec6 also
      IF (IFTOLD(KDSWAP).NE.IFTNEW(KDSWAP)) LTSCOPY = .TRUE.
C
      IF (LTSCOPY) THEN
      IF (LRETAG) THEN
      IST = -5
      ELSE
      IST = 0
      ENDIF
      LDUP = .FALSE.
      CALL zcorec6 (IFTOLD, IFTNEW, CPATH(1:NPATH), CPATH(1:NPATH),
     * IARRAY, KARRAY, IBUFF, KBUFF, LDUP, IST)
      IF (IFTNEW(KSTAT).NE.0) GO TO 940
      LPROTC = LPROT2
C     Re-read the record that we just destroyed
      CALL zgetad6 (IADD, IREC, 1)
      CALL zgtrec6 (IFTOLD, IARRAY, NBSIZE, IADD, .FALSE.)
      GO TO 200
C
      ENDIF
C
C
C     Save pertainent info from the info block
      NIHEAD = INFO(NPPWRD+KINIHE)
      NCHEAD = INFO(NPPWRD+KINCHE)
      NUHEAD = INFO(NPPWRD+KINUHE)
      NDATA  = INFO(NPPWRD+KINDAT)
      IIHADD = INFO(NPPWRD+KIAIHE)
      ICHADD = INFO(NPPWRD+KIACHE)
      IUHADD = INFO(NPPWRD+KIAUHE)
      IDADD  = INFO(NPPWRD+KIADAT)
      ITYPE  = INFO(NPPWRD+KITYPE)
      NLDATA = INFO(NPPWRD+KILNDA)
      IF (LRETAG) THEN
      CTAG = ' '
      ELSE
      CALL HOLCHR (INFO(NPPWRD+KITAG), 1, NTAGC, CTAG, 1)
      ENDIF
C
      NLEN = KIQUAL - KILNDA + 1
      DO 40 I=1,NLEN
      J = I + NPPWRD + KILNDA - 1
      INFOCO(I) = INFO(J)
 40   CONTINUE
      DO 45 I=1,NPASS
      J = I + NPPWRD + KIPASS - 1
      IPASS(I) = INFO(J)
 45   CONTINUE
C
C
C     Check for pseudo regular time series data
      IF ((ITYPE.GE.110).AND.(ITYPE.LE.120)) THEN
        IADD = IIHADD + 4
        CALL ZGTREC6 (IFTOLD, INTL_PSEUDO, 1, IADD, .FALSE.)
      ENDIF
C
C
C     Check for addresses out of range.
      IF ((IDADD.LT.100).OR.(IDADD.GE.IFTOLD(KFSIZE))) THEN
      IF (MLEVEL.GE.0) WRITE (MUNIT, 50) 'Data address',
     *                                    IADD, IDADD, CPATH(1:NPATH)
 50   FORMAT(/,' -----DSS*** zcofil6;  Error: ',A,' out of Range',
     * /,' Unable to Copy Record;',/,' Info Block Address:',I12,
     * ',   Value:',I12,/,' Pathname: ',A)
      GO TO 200
      ENDIF
      IF ((NIHEAD.GT.0).AND.
     *   ((IIHADD.GE.IFTOLD(KFSIZE)).OR.(IIHADD.LE.0))) THEN
      IF (MLEVEL.GE.0) WRITE (MUNIT, 50) 'Internal header address',
     *                                    IADD, IIHADD, CPATH(1:NPATH)
      GO TO 200
      ENDIF
      IF ((NCHEAD.GT.0).AND.
     *   ((ICHADD.GE.IFTOLD(KFSIZE)).OR.(ICHADD.LE.0))) THEN
      IF (MLEVEL.GE.0) WRITE (MUNIT, 50) 'Compression header address',
     *                                    IADD, ICHADD, CPATH(1:NPATH)
      GO TO 200
      ENDIF
      IF ((NUHEAD.GT.0).AND.
     *   ((IUHADD.GE.IFTOLD(KFSIZE)).OR.(IUHADD.LE.0))) THEN
      IF (MLEVEL.GE.0) WRITE (MUNIT, 50) 'User header address',
     *                                    IADD, IUHADD, CPATH(1:NPATH)
      GO TO 200
      ENDIF
      IF (NDATA.GT.10000) THEN
C        Large data set, probably an error.  For now, only check a
C        frequent UNET error
      IF (NDATA.EQ.IDADD) THEN
      IF (MLEVEL.GE.0) WRITE (MUNIT, 50) 'Number of data values',
     *                                    IADD, NDATA, CPATH(1:NPATH)
      GO TO 200
         ENDIF
      ENDIF
C
C     Do a sanity check on our numbers
      IF ((NIHEAD.GT.100000).OR.(NIHEAD.LT.0)) THEN
      IF (MLEVEL.GE.0) WRITE (MUNIT, 50) 'Internal Header values',
     *                                    IADD, NIHEAD, CPATH(1:NPATH)
      GO TO 200
      ENDIF
      IF ((NCHEAD.GT.100000).OR.(NCHEAD.LT.0)) THEN
      IF (MLEVEL.GE.0) WRITE (MUNIT, 50) 'Compression Header values',
     *                                    IADD, NCHEAD, CPATH(1:NPATH)
      GO TO 200
      ENDIF
      IF ((NUHEAD.GT.100000).OR.(NUHEAD.LT.0)) THEN
      IF (MLEVEL.GE.0) WRITE (MUNIT, 50) 'User Header values',
     *                                    IADD, NUHEAD, CPATH(1:NPATH)
      GO TO 200
      ENDIF
C
C     Write the data record to the DSS file w/ IFTNEW
C
C     Get multiple user access
      LWRITE = .FALSE.
      CALL zmultu6 ( IFTNEW, .TRUE., .TRUE.)
C
C     Check if new record exists
      CALL zcheck6 (IFTNEW, CPATH(1:NPATH), NPATH, JHEAD, JDATA, LFOUND)
      IF (IFTNEW(KSTAT).NE.0) GO TO 940
C
C
C     If the pathname was not found by zcheck6 write new pointers
C
      IF (.NOT.LFOUND) THEN
      CALL znwrit6 (IFTNEW, CPATH(1:NPATH), NPATH, NIHEAD, NCHEAD,
     * NUHEAD, NDATA)
      ELSE
      IF (LPROT2) THEN
      IF (MLEVEL.GE.2) WRITE (MUNIT, 60) CPATH(1:NPATH)
 60   FORMAT (' -----DSS---zcofil6:  Write Protection for Existing',
     * ' Record (no data written)',/,
     * ' Pathname: ',A)
      CALL zmultu6 ( IFTNEW, .FALSE., .TRUE.)
      GO TO 200
      ENDIF
      CALL zowrit6 (IFTNEW, CPATH(1:NPATH), NPATH, NIHEAD, NCHEAD,
     * NUHEAD, NDATA)
      ENDIF
C
C
C     Update the information block to contain what the old one had
C     But don't copy the tag
C
      NLEN = KITIME - KILNDA + 1
      DO 80 I=1,NLEN
      J = I + NPPWRD + KILNDA - 1
      INFO(J) = INFOCO(I)
 80   CONTINUE
C
      NLEN = KIQUAL - KILNDA + 1
      IBEG = KITYPE - KILNDA + 1
      DO 90 I=IBEG,NLEN
      J = I + NPPWRD + KILNDA - 1
      INFO(J) = INFOCO(I)
 90   CONTINUE
      DO 100 I=1,NPASS
      J = I + NPPWRD + KIPASS - 1
      INFO(J) = IPASS(I)
 100  CONTINUE
C
      ISIZE = NPPWRD + NINFO
      CALL zptrec6 (IFTNEW, INFO, ISIZE, IPNBIN(JPNBIN+NPPWRD+KBAINF),
     * .FALSE.)
      IF (IFTNEW(KSTAT).NE.0) GO TO 940
C
C     Now copy the headers and the data array
C
C     Internal header
      IF (NIHEAD.GT.0) THEN
      N = MIN0(NIHEAD,KBUFF)
      NTOT = 0
      IADD = INFO(NPPWRD+KIAIHE)
 120  CONTINUE
      CALL zgtrec6 (IFTOLD, IBUFF, N, IIHADD, .FALSE.)
      CALL zptrec6 (IFTNEW, IBUFF, N, IADD,   .FALSE.)
      NTOT = NTOT + N
      IF (NTOT.LT.NIHEAD) THEN
      IIHADD = IIHADD + N
      IADD = IADD + N
      N = NIHEAD - NTOT
      N = MIN0(N,KBUFF)
      GO TO 120
      ENDIF
      ENDIF
C
C     Compression header
      IF (NCHEAD.GT.0) THEN
      N = MIN0(NCHEAD,KBUFF)
      NTOT = 0
      IADD = INFO(NPPWRD+KIACHE)
 140  CONTINUE
      CALL zgtrec6 (IFTOLD, IBUFF, N, ICHADD, .FALSE.)
      CALL zptrec6 (IFTNEW, IBUFF, N, IADD,   .FALSE.)
      NTOT = NTOT + N
      IF (NTOT.LT.NCHEAD) THEN
      ICHADD = ICHADD + N
      IADD = IADD + N
      N = NCHEAD - NTOT
      N = MIN0(N,KBUFF)
      GO TO 140
      ENDIF
      ENDIF
C
C     User Header
      IF (NUHEAD.GT.0) THEN
      N = MIN0(NUHEAD,KBUFF)
      NTOT = 0
      IADD = INFO(NPPWRD+KIAUHE)
 160  CONTINUE
      CALL zgtrec6 (IFTOLD, IBUFF, N, IUHADD, .FALSE.)
      CALL zptrec6 (IFTNEW, IBUFF, N, IADD,   .FALSE.)
      NTOT = NTOT + N
      IF (NTOT.LT.NUHEAD) THEN
      IUHADD = IUHADD + N
      IADD = IADD + N
      N = NUHEAD - NTOT
      N = MIN0(N,KBUFF)
      GO TO 160
      ENDIF
      ENDIF
C
C     Data Array
      IF (NDATA.GT.0) THEN
      N = MIN0(NDATA,KBUFF)
      NTOT = 0
      IADD = INFO(NPPWRD+KIADAT)
 180  CONTINUE
      CALL zgtrec6 (IFTOLD, IBUFF, N, IDADD, .FALSE.)
      CALL zptrec6 (IFTNEW, IBUFF, N, IADD,   .FALSE.)
      NTOT = NTOT + N
      IF (NTOT.LT.NDATA) THEN
      IDADD = IDADD + N
      IADD = IADD + N
      N = NDATA - NTOT
      N = MIN0(N,KBUFF)
      GO TO 180
      ENDIF
      ENDIF
C
      IF (MLEVEL.GE.3) WRITE (MUNIT,190) CPATH(1:NPATH)
 190  FORMAT (' -----DSS--- zcofil6;  Record: ',A)
C
C     Release multiple user access
      CALL zmultu6 ( IFTNEW, .FALSE., .TRUE.)
      LWRITE = .FALSE.
C
      ENDIF
      ENDIF
      ENDIF
C
      IF (IFTOLD(KSTAT).NE.0) GO TO 940
      IF (IFTNEW(KSTAT).NE.0) GO TO 940
C
 200  CONTINUE
C
C
C     IF (LSQSTA) CALL CHRWT (MUNIT, '100' // CHAR(13) // CHAR(10), 6)  u
      !IF (LSQSTA) CALL CHRWT (MUNIT, '+100' // CHAR(13) // CHAR(10), 6) Md
C
 800  CONTINUE
      LSQSTA = .FALSE.
      LTSCMP = .FALSE.
      LPROTC = .FALSE.
      LCOFIL = .FALSE.
      IF (MLEVEL.GE.11) WRITE (MUNIT,820)
 820  FORMAT (T6,'-----DSS---Debug:  Exit zcofil6')
      !CALL FLUSH(MUNIT)                                                 Mu
C
      RETURN
C
C
 840  CONTINUE
      IF (MLEVEL.GE.4) WRITE (MUNIT, 850)
 850  FORMAT (' ----- DSS --- zcofil6:  Copy interrupted.')
      GO TO 800
C
C
 900  CONTINUE
      WRITE (MUNIT, 901) KARRAY, KBUFF, NBSIZE
 901  FORMAT (/,' ----- DSS --- zcofil6:  ERROR;  Buffer size(s) not',
     * ' Large enough to copy file',/,' Sizes provided:',2I8,
     * ';   Size Required:',I5,/)
      CALL zabort6 (IFTOLD, 331, 'zcofil6 ', KBUFF, NBSIZE,
     * 'Buffer sizes not large enough')
      RETURN
C
 910  CONTINUE
      WRITE (MUNIT, 911) IFTOLD(1), IFTNEW(1)
 911  FORMAT (/,' ----- DSS --- zcofil6:  ERROR;  DSS File(s) not',
     * ' Version 6.',/,' Versions:',2I8,/)
      RETURN
C
C
 920  CONTINUE
      WRITE (MUNIT, 921)
 921  FORMAT(/' -----DSS---zcofil6: ERROR;  File has Read Access Only')
      RETURN
C
C
 930  CONTINUE
      WRITE (MUNIT, 931)
 931  FORMAT (/' -----DSS---zcofil6:  ERROR;  Files are of different',
     * ' endian type.'/' Can only copy files with the same byte order')
      RETURN
C
C
 940  CONTINUE
      WRITE (MUNIT, 941) IFTNEW(KSTAT), IFTOLD(KSTAT)
 941  FORMAT (/,' *****DSS*** zcofil6:  ERROR  - UNABLE TO ',
     * ' COPY DATA',/,' Status: ',2I8)
      RETURN
C
      END

ext: f

max_stars_repo_name: sharkcooper/hec-dss

max_issues_repo_name: sharkcooper/hec-dss

max_forks_repo_name: sharkcooper/hec-dss


text:
*> \brief \b IZAMAX
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*  Definition:
*  ===========
*
*       INTEGER FUNCTION IZAMAX(N,ZX,INCX)
*
*       .. Scalar Arguments ..
*       INTEGER INCX,N
*       ..
*       .. Array Arguments ..
*       COMPLEX*16 ZX(*)
*       ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*>    IZAMAX finds the index of element having max. absolute value.
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \date November 2011
*
*> \ingroup aux_blas
*
*> \par Further Details:
*  =====================
*>
*> \verbatim
*>
*>     jack dongarra, 1/15/85.
*>     modified 3/93 to return if incx .le. 0.
*>     modified 12/3/93, array(1) declarations changed to array(*)
*> \endverbatim
*>
*  =====================================================================
      INTEGER FUNCTION IZAMAX(N,ZX,INCX)
*
*  -- Reference BLAS level1 routine (version 3.4.0) --
*  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*     November 2011
*
*     .. Scalar Arguments ..
      INTEGER INCX,N
*     ..
*     .. Array Arguments ..
      COMPLEX*16 ZX(*)
*     ..
*
*  =====================================================================
*
*     .. Local Scalars ..
      DOUBLE PRECISION DMAX
      INTEGER I,IX
*     ..
*     .. External Functions ..
      DOUBLE PRECISION DCABS1
      EXTERNAL DCABS1
*     ..
      IZAMAX = 0
      IF (N.LT.1 .OR. INCX.LE.0) RETURN
      IZAMAX = 1
      IF (N.EQ.1) RETURN
      IF (INCX.EQ.1) THEN
*
*        code for increment equal to 1
*
         DMAX = DCABS1(ZX(1))
         DO I = 2,N
            IF (DCABS1(ZX(I)).GT.DMAX) THEN
               IZAMAX = I
               DMAX = DCABS1(ZX(I))
            END IF
         END DO
      ELSE
*
*        code for increment not equal to 1
*
         IX = 1
         DMAX = DCABS1(ZX(1))
         IX = IX + INCX
         DO I = 2,N
            IF (DCABS1(ZX(IX)).GT.DMAX) THEN
               IZAMAX = I
               DMAX = DCABS1(ZX(IX))
            END IF
            IX = IX + INCX
         END DO
      END IF
      RETURN
      END
c $Id$

max_stars_repo_name: dinisAbranches/nwchem

max_issues_repo_name: dinisAbranches/nwchem

max_forks_repo_name: dinisAbranches/nwchem


text:
      subroutine ZZmbc(st,j1,j2,j3,j4,j5,j6,za,zb,bcoeff)
      implicit none
c----Bubble coefficients extracted from BDK 11.5, 11.8
c----after performing a transformation
      integer j1,j2,j3,j4,j5,j6
      include 'constants.f'
      include 'zprods_decl.f'
      include 'sprods_com.f'
      include 'blabels.f'
      logical ggZZuse6d
      character*9 st
      double complex bcoeff(8),zab2,funcb34pp,funcb134pp
      double complex funcb134mp,funcb34mp,
     & oldfuncb134mp,oldfuncb34mp
      double precision IDelta
      common/ggZZuse6d/ggZZuse6d
!$omp threadprivate(/ggZZuse6d/)

C----statement functions
      zab2(j1,j2,j3,j4)=+za(j1,j2)*zb(j2,j4)+za(j1,j3)*zb(j3,j4)
      funcb34pp(j1,j2,j3,j4,j5,j6)=za(j3,j4)/(za(j1,j2)**2*za(j5,j6))
     & *((za(j1,j5)*zb(j1,j4)/(s(j1,j3)+s(j1,j4)))**2
     &  +(za(j2,j5)*zb(j2,j4)/(s(j2,j3)+s(j2,j4)))**2)
     & +2d0*za(j1,j5)*za(j2,j5)/(za(j1,j2)**3*za(j5,j6))
     & *(za(j3,j1)*zb(j4,j1)/(s(j1,j3)+s(j1,j4))
     &  +za(j2,j3)*zb(j4,j2)/(s(j2,j3)+s(j2,j4)))

      funcb134pp(j1,j2,j3,j4,j5,j6)=(2d0*za(j1,j3)*za(j2,j5)
     & *(za(j2,j3)*zb(j2,j6)*za(j5,j6)/(s(j2,j5)+s(j2,j6))
     &  -za(j1,j5)*zb(j1,j4)*za(j3,j4)/(s(j1,j3)+s(j1,j4)))
     & -za(j1,j2)
     & *((za(j3,j4)*za(j1,j5)*zb(j1,j4)/(s(j1,j3)+s(j1,j4)))**2
     &  +(za(j2,j3)*za(j5,j6)*zb(j2,j6)/(s(j2,j5)+s(j2,j6)))**2)
     & -(za(j1,j3)*za(j2,j5)+za(j2,j3)*za(j1,j5))*za(j3,j5))
     & /(za(j1,j2)**3*za(j3,j4)*za(j5,j6))
c---end statement functions

      IDelta=s(j1,j2)**2+s(j3,j4)**2+s(j5,j6)**2
     & -2d0*s(j1,j2)*s(j3,j4)
     & -2d0*s(j3,j4)*s(j5,j6)
     & -2d0*s(j1,j2)*s(j5,j6)
      IDelta=1d0/IDelta

      bcoeff(:)=czip

      if     (st.eq.'q+qb-g+g+') then

      bcoeff(b34)=funcb34pp(j1,j2,j3,j4,j5,j6)
      bcoeff(b56)=funcb34pp(j1,j2,j5,j6,j3,j4)
      bcoeff(b134)=funcb134pp(j1,j2,j3,j4,j5,j6)
c      bcoeff(b234)=funcb134pp(j2,j1,j3,j4,j5,j6)
      bcoeff(b234)=-bcoeff(b34)-bcoeff(b56)-bcoeff(b134)

      bcoeff(rat)=(
     &  (za(j1,j5)**2*zb(j1,j4)**2/(s(j1,j3)+s(j1,j4))
     &  +za(j2,j5)**2*zb(j2,j4)**2/(s(j2,j3)+s(j2,j4)))
     & /(za(j5,j6)*zb(j3,j4))
     & +(za(j3,j1)**2*zb(j1,j6)**2/(s(j1,j5)+s(j1,j6))
     &  +za(j3,j2)**2*zb(j2,j6)**2/(s(j2,j5)+s(j2,j6)))
     & /(zb(j5,j6)*za(j3,j4))
     & -za(j3,j5)**2/(za(j3,j4)*za(j5,j6))
     & +zb(j4,j6)**2/(zb(j3,j4)*zb(j5,j6)))/za(j1,j2)**2

      elseif (st.eq.'q+qb-g-g+') then

c      bcoeff(b12zm)= + IDelta**2 * ( 3.D0*zba2(j4,j1,j2,j3)*zba2(j2,j3,
C     &    j4,j1)*zba2(j6,j1,j2,j5)/zba2(j1,j3,j4,j2)*s(j3,j4) + 3.D0*
C     &    zba2(j4,j1,j2,j3)*zba2(j2,j4,j3,j1)*zba2(j6,j1,j2,j5)/zba2(
C     &    j1,j3,j4,j2)*s(j1,j2) - 3.D0*zba2(j4,j1,j2,j3)*zba2(j2,j3,j4,
C     &    j1)*zba2(j6,j1,j2,j5)/zba2(j1,j3,j4,j2)*s(j5,j6) - 3.D0*
C     &    zba2(j4,j1,j2,j3)*zba2(j2,j5,j6,j1)*zba2(j6,j1,j2,j5)/zba2(
C     &    j1,j5,j6,j2)*s(j3,j4) + 3.D0*zba2(j4,j1,j2,j3)*zba2(j2,j5,j6,
C     &    j1)*zba2(j6,j1,j2,j5)/zba2(j1,j5,j6,j2)*s(j1,j2) + 3.D0*
C     &    zba2(j4,j1,j2,j3)*zba2(j2,j5,j6,j1)*zba2(j6,j1,j2,j5)/zba2(
C     &    j1,j5,j6,j2)*s(j5,j6) + 3.D0*zba2(j4,j1,j2,j3)*zba2(j2,j4,j3,
C     &    j1)*zba2(j6,j1,j2,j5)/zab2(j2,j3,j4,j1)*s(j3,j4) + 3.D0*
C     &    zba2(j4,j1,j2,j3)*zba2(j2,j4,j3,j1)*zba2(j6,j1,j2,j5)/zba2(
C     &    j1,j4,j3,j2)*s(j1,j2) - 3.D0*zba2(j4,j1,j2,j3)*zba2(j2,j4,j3,
C     &    j1)*zba2(j6,j1,j2,j5)/zab2(j2,j3,j4,j1)*s(j5,j6) - 3.D0*
C     &    zba2(j4,j1,j2,j3)*zba2(j2,j6,j5,j1)*zba2(j6,j1,j2,j5)/zba2(
C     &    j1,j6,j5,j2)*s(j3,j4) )
C     & + IDelta**2 * ( 3.D0*zba2(j4,j2,j1,
C     &    j3)*zba2(j2,j6,j5,j1)*zba2(j6,j1,j2,j5)/zba2(j1,j6,j5,j2)*s(
C     &    j1,j2) + 3.D0*zba2(j4,j1,j2,j3)*zba2(j2,j6,j5,j1)*zba2(j6,j2,
C     &    j1,j5)/zba2(j1,j6,j5,j2)*s(j5,j6) )
C     & + IDelta * (  - 2.D0*zb(j3,j4)*zb(
C     &    j4,j1)*zb(j6,j5)*za(j3,j4)*za(j3,j5)*za(j2,j5)/zba2(j1,j3,j4
C     &    ,j2)**3*t(j3,j4,j1) + 2.D0*zb(j3,j4)*zb(j4,j1)*zb(j6,j5)*za(
C     &    j3,j4)*za(j3,j5)*za(j2,j5)/zba2(j1,j3,j4,j2)**3*t(j3,j4,j2)
C     &     - 2.D0*zb(j3,j4)*zb(j4,j6)*zb(j1,j6)*za(j3,j4)*za(j3,j2)*za(
C     &    j6,j5)/zab2(j2,j3,j4,j1)**3*t(j4,j3,j1) + 2.D0*zb(j3,j4)*zb(
C     &    j4,j6)*zb(j1,j6)*za(j3,j4)*za(j3,j2)*za(j6,j5)/zba2(j1,j4,j3
C     &    ,j2)**3*t(j4,j3,j2) + zb(j3,j4)*zb(j4,j6)*za(j3,j4)*za(j3,j5)
C     &    /zba2(j1,j3,j4,j2)**2*s(j3,j4) - zb(j3,j4)*zb(j4,j6)*za(j3,
C     &    j4)*za(j3,j5)/zba2(j1,j3,j4,j2)**2*s(j1,j2) - zb(j3,j4)*zb(
C     &    j4,j6)*za(j3,j4)*za(j3,j5)/zba2(j1,j3,j4,j2)**2*s(j5,j6) +
C     &    zb(j3,j4)*zb(j4,j6)*za(j3,j4)*za(j3,j5)/zab2(j2,j3,j4,j1)**2
C     &    *s(j3,j4) - zb(j3,j4)*zb(j4,j6)*za(j3,j4)*za(j3,j5)/zba2(j1,
C     &    j4,j3,j2)**2*s(j1,j2) - zb(j3,j4)*zb(j4,j6)*za(j3,j4)*za(j3,
C     &    j5)/zab2(j2,j3,j4,j1)**2*s(j5,j6) - 2.D0*zb(j3,j4)*zb(j1,j6)
C     &    *zb(j6,j5)*za(j3,j2)*za(j3,j5)*za(j6,j5)/zba2(j1,j5,j6,j2)**
C     &    3*t(j5,j6,j1) )
C     & + IDelta * ( 2.D0*zb(j3,j4)*zb(j1,
C     &    j6)*zb(j6,j5)*za(j3,j2)*za(j3,j5)*za(j6,j5)/zba2(j1,j5,j6,j2
C     &    )**3*t(j5,j6,j2) - 2.D0*zb(j3,j4)*zb(j1,j6)*za(j3,j2)*za(j3,
C     &    j5)/zab2(j2,j3,j4,j1)**3*t(j4,j3,j1)*t(j4,j3,j2) + 2.D0*zb(
C     &    j3,j4)*zb(j1,j6)*za(j3,j2)*za(j3,j5)/zab2(j2,j3,j4,j1)**3*t(
C     &    j4,j3,j2)**2 + 2.D0*zb(j3,j4)*zb(j2,j6)*za(j3,j2)*za(j3,j5)*
C     &    izab2(j2,j3,j4,j1)**2*t(j4,j3,j1) - 2.D0*zb(j3,j4)*zb(j2,j6)*
C     &    za(j3,j2)*za(j3,j5)/zab2(j2,j3,j4,j1)**2*t(j4,j3,j2) - 1d0/
C     &    2.D0*zb(j3,j4)*za(j3,j5)**2/za(j6,j5)*zba2(j2,j4,j3,j1)*
C     &    izba2(j1,j3,j4,j2) - 1d0/2.D0*zb(j3,j4)*za(j3,j5)**2/za(j6,
C     &    j5)*zba2(j2,j4,j3,j1)/zab2(j2,j3,j4,j1) + zb(j3,j4)*za(j3,j5
C     &    )**2/za(j6,j5)/zab2(j2,j3,j4,j1)**2*s(j3,j4)*t(j4,j3,j2) -
C     &    zb(j3,j4)*za(j3,j5)**2/za(j6,j5)/zab2(j2,j3,j4,j1)**2*s(j1,
C     &    j2)*t(j4,j3,j2) - zb(j3,j4)*za(j3,j5)**2/za(j6,j5)/zba2(j1,
C     &    j4,j3,j2)**2*s(j5,j6)*t(j4,j3,j2) - 2.D0*zb(j4,j1)*zb(j4,j6)*
C     &    zb(j6,j5)*za(j3,j4)*za(j2,j5)*za(j6,j5)/zba2(j1,j6,j5,j2)**3
C     &    *t(j6,j5,j1) )
C     & + IDelta * ( 2.D0*zb(j4,j1)*zb(j4,
C     &    j6)*zb(j6,j5)*za(j3,j4)*za(j2,j5)*za(j6,j5)/zba2(j1,j6,j5,j2
C     &    )**3*t(j6,j5,j2) + 2.D0*zb(j4,j1)*zb(j4,j6)*za(j3,j4)*za(j1,
C     &    j5)/zba2(j1,j3,j4,j2)**2*t(j3,j4,j1) - 2.D0*zb(j4,j1)*zb(j4,
C     &    j6)*za(j3,j4)*za(j1,j5)/zba2(j1,j3,j4,j2)**2*t(j3,j4,j2) - 2.
C     &    D0*zb(j4,j1)*zb(j4,j6)*za(j3,j4)*za(j2,j5)/zba2(j1,j3,j4,j2)
C     &    **3*t(j3,j4,j1)**2 + 2.D0*zb(j4,j1)*zb(j4,j6)*za(j3,j4)*za(j2
C     &    ,j5)/zba2(j1,j3,j4,j2)**3*t(j3,j4,j1)*t(j3,j4,j2) - 2.D0*zb(
C     &    j4,j1)*zb(j6,j5)*za(j3,j5)*za(j2,j5)/zba2(j1,j6,j5,j2)**3*t(
C     &    j6,j5,j1)*t(j6,j5,j2) + 2.D0*zb(j4,j1)*zb(j6,j5)*za(j3,j5)*
C     &    za(j2,j5)/zba2(j1,j6,j5,j2)**3*t(j6,j5,j2)**2 + 2.D0*zb(j4,
C     &    j2)*zb(j6,j5)*za(j3,j5)*za(j2,j5)/zba2(j1,j6,j5,j2)**2*t(j6,
C     &    j5,j1) - 2.D0*zb(j4,j2)*zb(j6,j5)*za(j3,j5)*za(j2,j5)/zba2(
C     &    j1,j6,j5,j2)**2*t(j6,j5,j2) - 1d0/2.D0*zb(j4,j6)**2*za(j3,j4
C     &    )/zb(j6,j5)*zba2(j2,j4,j3,j1)/zba2(j1,j3,j4,j2) - 1d0/2.D0
C     &    *zb(j4,j6)**2*za(j3,j4)/zb(j6,j5)*zba2(j2,j4,j3,j1)/zba2(j1
C     &    ,j4,j3,j2) )
C     & + IDelta * ( zb(j4,j6)**2*za(j3,j4)
C     &    /zb(j6,j5)/zba2(j1,j3,j4,j2)**2*s(j3,j4)*t(j3,j4,j1) - zb(
C     &    j4,j6)**2*za(j3,j4)/zb(j6,j5)/zba2(j1,j3,j4,j2)**2*s(j1,j2)
C     &    *t(j3,j4,j1) - zb(j4,j6)**2*za(j3,j4)/zb(j6,j5)/zba2(j1,j3,
C     &    j4,j2)**2*s(j5,j6)*t(j3,j4,j1) - 1d0/2.D0*zb(j4,j6)**2*za(j6
C     &    ,j5)/zb(j3,j4)*zba2(j2,j6,j5,j1)/zba2(j1,j6,j5,j2) - 1d0/2.
C     &    D0*zb(j4,j6)**2*za(j6,j5)/zb(j3,j4)*zba2(j2,j5,j6,j1)/zba2(
C     &    j1,j5,j6,j2) - zb(j4,j6)**2*za(j6,j5)/zb(j3,j4)/zba2(j1,j5,
C     &    j6,j2)**2*s(j3,j4)*t(j5,j6,j1) - zb(j4,j6)**2*za(j6,j5)/zb(
C     &    j3,j4)/zba2(j1,j5,j6,j2)**2*s(j1,j2)*t(j5,j6,j1) + zb(j4,j6)
C     &    **2*za(j6,j5)/zb(j3,j4)/zba2(j1,j5,j6,j2)**2*s(j5,j6)*t(j5,
C     &    j6,j1) + 2.D0*zb(j4,j6)*zb(j1,j6)*za(j3,j1)*za(j6,j5)/zba2(
C     &    j1,j5,j6,j2)**2*t(j5,j6,j1) - 2.D0*zb(j4,j6)*zb(j1,j6)*za(j3,
C     &    j1)*za(j6,j5)/zba2(j1,j5,j6,j2)**2*t(j5,j6,j2) - 2.D0*zb(j4,
C     &    j6)*zb(j1,j6)*za(j3,j2)*za(j6,j5)/zba2(j1,j5,j6,j2)**3*t(j5,
C     &    j6,j1)**2 )
C     & + IDelta * ( 2.D0*zb(j4,j6)*zb(j1,
C     &    j6)*za(j3,j2)*za(j6,j5)/zba2(j1,j5,j6,j2)**3*t(j5,j6,j1)*t(
C     &    j5,j6,j2) - zb(j4,j6)*zb(j6,j5)*za(j3,j5)*za(j6,j5)/zba2(j1,
C     &    j6,j5,j2)**2*s(j3,j4) - zb(j4,j6)*zb(j6,j5)*za(j3,j5)*za(j6,
C     &    j5)/zba2(j1,j6,j5,j2)**2*s(j1,j2) + zb(j4,j6)*zb(j6,j5)*za(
C     &    j3,j5)*za(j6,j5)/zba2(j1,j6,j5,j2)**2*s(j5,j6) - zb(j4,j6)*
C     &    zb(j6,j5)*za(j3,j5)*za(j6,j5)/zba2(j1,j5,j6,j2)**2*s(j3,j4)
C     &     - zb(j4,j6)*zb(j6,j5)*za(j3,j5)*za(j6,j5)/zba2(j1,j5,j6,j2)
C     &    **2*s(j1,j2) + zb(j4,j6)*zb(j6,j5)*za(j3,j5)*za(j6,j5)/zba2(
C     &    j1,j5,j6,j2)**2*s(j5,j6) + zb(j4,j6)*za(j3,j5)*zba2(j2,j3,j4,
C     &    j1)/zba2(j1,j3,j4,j2) + zb(j4,j6)*za(j3,j5)*zba2(j2,j4,j3,j1
C     &    )/zab2(j2,j3,j4,j1) + zb(j4,j6)*za(j3,j5)*zba2(j2,j6,j5,j1)*
C     &    izba2(j1,j6,j5,j2) + zb(j4,j6)*za(j3,j5)*zba2(j2,j5,j6,j1)*
C     &    izba2(j1,j5,j6,j2) - 1d0/2.D0*zb(j6,j5)*za(j3,j5)**2/za(j3,
C     &    j4)*zba2(j2,j6,j5,j1)/zba2(j1,j6,j5,j2) - 1d0/2.D0*zb(j6,j5
C     &    )*za(j3,j5)**2/za(j3,j4)*zba2(j2,j5,j6,j1)/zba2(j1,j5,j6,j2
C     &    ) )
C     & + IDelta * (  - zb(j6,j5)*za(j3,j5)
C     &    **2/za(j3,j4)/zba2(j1,j6,j5,j2)**2*s(j3,j4)*t(j6,j5,j2) -
C     &    zb(j6,j5)*za(j3,j5)**2/za(j3,j4)/zba2(j1,j6,j5,j2)**2*s(j1,
C     &    j2)*t(j6,j5,j2) + zb(j6,j5)*za(j3,j5)**2/za(j3,j4)/zba2(j1,
C     &    j6,j5,j2)**2*s(j5,j6)*t(j6,j5,j2) )
C     & - 2.D0*zb(j4,j1)*zb(j4,j2)*za(j1,j5
C     & )*za(j2,j5)/zb(j3,j4)/za(j6,j5)/zba2(j1,j3,j4,j2)**2 + 2.D0*
C     &    zb(j4,j1)*za(j2,j5)/zb(j3,j4)/za(j6,j5)*zba2(j4,j3,j2,j5)*
C     &    izba2(j1,j3,j4,j2)**3*t(j3,j4,j1) - 2.D0*zb(j1,j6)*zb(j2,j6)*
C     &    za(j3,j1)*za(j3,j2)/zb(j6,j5)/za(j3,j4)/zab2(j2,j3,j4,j1)
C     &    **2 + 2.D0*zb(j1,j6)*za(j3,j2)/zb(j6,j5)/za(j3,j4)*zba2(j6,
C     &    j4,j1,j3)/zab2(j2,j3,j4,j1)**3*t(j4,j3,j2) - izb(j3,j4)/za(
C     &    j6,j5)*zba2(j4,j3,j2,j5)**2/zba2(j1,j3,j4,j2)**2 - izb(j6,j5
C     &    )/za(j3,j4)*zba2(j6,j4,j1,j3)**2/zab2(j2,j3,j4,j1)**2

c      write(6,*) 'old:b12zm',bcoeff(b12zm)

c      bcoeff(b12zm)= + IDelta**2
c     & * ( 3.D0*zba2(j4,j1,j2,j3)*zab2(j1,j3,j4,j2)*zba2(j6,j1,j2,j5)
c     & /zba2(j1,j3,j4,j2)*s(j3,j4)
c     &  + 3.D0*zba2(j4,j1,j2,j3)*zba2(j2,j4,j3,j1)*zba2(j6,j1,j2,j5)
c     & /zba2(j1,j3,j4,j2)*s(j1,j2)
c     &  -3.D0*zba2(j4,j1,j2,j3)*zab2(j1,j3,j4,j2)*zba2(j6,j1,j2,j5)
c     & /zba2(j1,j3,j4,j2)*s(j5,j6)
c     & - 3.D0*zba2(j4,j1,j2,j3)*zba2(j2,j5,j6,j1)*zba2(j6,j1,j2,j5)
c     & /zba2(j1,j5,j6,j2)*s(j3,j4)
c     &  + 3.D0*zba2(j4,j1,j2,j3)*zba2(j2,j5,j6,j1)*zba2(j6,j1,j2,j5)
c     & /zba2(j1,j5,j6,j2)*s(j1,j2)
c     &  + 3.D0*zba2(j4,j1,j2,j3)*zba2(j2,j5,j6,j1)*zba2(j6,j1,j2,j5)
c     & /zba2(j1,j5,j6,j2)*s(j5,j6)
c     &  + 3.D0*zba2(j4,j1,j2,j3)*zba2(j2,j4,j3,j1)*zba2(j6,j1,j2,j5)
c     & /zab2(j2,j3,j4,j1)*s(j3,j4)
c     & + 3.D0*zba2(j4,j1,j2,j3)*zba2(j2,j4,j3,j1)*zba2(j6,j1,j2,j5)
c     & /zab2(j2,j3,j4,j1)*s(j1,j2)
c     & - 3.D0*zba2(j4,j1,j2,j3)*zba2(j2,j4,j3,j1)*zba2(j6,j1,j2,j5)
c     & /zab2(j2,j3,j4,j1)*s(j5,j6)
c     &  - 3.D0*zba2(j4,j1,j2,j3)*zba2(j2,j6,j5,j1)*zba2(j6,j1,j2,j5)
c     & /zba2(j1,j6,j5,j2)*s(j3,j4))

c     & + IDelta**2
c     & * ( 3.D0*zba2(j4,j2,j1,j3)*zba2(j2,j6,j5,j1)*zba2(j6,j1,j2,j5)
c     & /zba2(j1,j6,j5,j2)*s(j1,j2)
c     &  + 3.D0*zba2(j4,j1,j2,j3)*zba2(j2,j6,j5,j1)*zba2(j6,j2,j1,j5)
c     & /zba2(j1,j6,j5,j2)*s(j5,j6))

c     & + IDelta
c     &  * (-2.D0*zb(j3,j4)*zb(j4,j1)*zb(j6,j5)*za(j3,j4)*za(j3,j5)
c     & *za(j2,j5)/zba2(j1,j3,j4,j2)**3*t(j3,j4,j1)
c     &  + 2.D0*zb(j3,j4)*zb(j4,j1)*zb(j6,j5)*za(j3,j4)*za(j3,j5)
c     & *za(j2,j5)/zba2(j1,j3,j4,j2)**3*t(j3,j4,j2)
c     & -2.D0*zb(j3,j4)*zb(j4,j6)*zb(j1,j6)*za(j3,j4)*za(j3,j2)
c     & *za(j6,j5)/zab2(j2,j3,j4,j1)**3*t(j4,j3,j1)
c     &  +2.D0*zb(j3,j4)*zb(j4,j6)*zb(j1,j6)*za(j3,j4)*za(j3,j2)
c     &  *za(j6,j5)/zab2(j2,j3,j4,j1)**3*t(j4,j3,j2)
c     &  + zb(j3,j4)*zb(j4,j6)*za(j3,j4)*za(j3,j5)
c     & /zba2(j1,j3,j4,j2)**2*s(j3,j4)
c     & - zb(j3,j4)*zb(j4,j6)*za(j3,j4)*za(j3,j5)
c     & /zba2(j1,j3,j4,j2)**2*s(j1,j2)
c     &  - zb(j3,j4)*zb(j4,j6)*za(j3,j4)*za(j3,j5)
c     & /zba2(j1,j3,j4,j2)**2*s(j5,j6)
c     &  + zb(j3,j4)*zb(j4,j6)*za(j3,j4)*za(j3,j5)
c     & /zab2(j2,j3,j4,j1)**2*s(j3,j4)
c     &  -zb(j3,j4)*zb(j4,j6)*za(j3,j4)*za(j3,j5)
c     & /zab2(j2,j3,j4,j1)**2*s(j1,j2)
c     &  -zb(j3,j4)*zb(j4,j6)*za(j3,j4)*za(j3,j5)
c     & /zab2(j2,j3,j4,j1)**2*s(j5,j6)
c     &  -2.D0*zb(j3,j4)*zb(j1,j6)*zb(j6,j5)*za(j3,j2)*za(j3,j5)
c     & *za(j6,j5)/zba2(j1,j5,j6,j2)**3*t(j5,j6,j1))

c     & + IDelta * (
c     &  2.D0*zb(j3,j4)*zb(j1,j6)*zb(j6,j5)*za(j3,j2)*za(j3,j5)
c     &  *za(j6,j5)/zba2(j1,j5,j6,j2)**3*t(j5,j6,j2)
c     &  -2.D0*zb(j3,j4)*zb(j1,j6)*za(j3,j2)*za(j3,j5)
c     & /zab2(j2,j3,j4,j1)**3*t(j4,j3,j1)*t(j4,j3,j2)
c     & + 2.D0*zb(j3,j4)*zb(j1,j6)*za(j3,j2)*za(j3,j5)
c     & /zab2(j2,j3,j4,j1)**3*t(j4,j3,j2)**2
c     &  + 2.D0*zb(j3,j4)*zb(j2,j6)*za(j3,j2)*za(j3,j5)*
c     &    izab2(j2,j3,j4,j1)**2*t(j4,j3,j1)
c     &  -2.D0*zb(j3,j4)*zb(j2,j6)*za(j3,j2)*za(j3,j5)
c     & /zab2(j2,j3,j4,j1)**2*t(j4,j3,j2)
c     &  - 1d0/2.D0*zb(j3,j4)*za(j3,j5)**2/za(j6,j5)*zba2(j2,j4,j3,j1)*
c     &    izba2(j1,j3,j4,j2)
c     & - 1d0/2.D0*zb(j3,j4)*za(j3,j5)**2/za(j6,j5)*zba2(j2,j4,j3,j1)
c     &  /zab2(j2,j3,j4,j1)
c     & + zb(j3,j4)*za(j3,j5)**2/za(j6,j5)
c     & /zab2(j2,j3,j4,j1)**2*s(j3,j4)*t(j4,j3,j2)
c     & -zb(j3,j4)*za(j3,j5)**2
c     & /za(j6,j5)/zab2(j2,j3,j4,j1)**2*s(j1,j2)*t(j4,j3,j2)
c     &  -zb(j3,j4)*za(j3,j5)**2/za(j6,j5)/zab2(j2,j3,j4,j1)**2
c     & *s(j5,j6)*t(j4,j3,j2)
c     &  -2.D0*zb(j4,j1)*zb(j4,j6)*zb(j6,j5)*za(j3,j4)*za(j2,j5)
c     & *za(j6,j5)/zba2(j1,j6,j5,j2)**3*t(j6,j5,j1))

c     & + IDelta
c     & * ( 2.D0*zb(j4,j1)*zb(j4,j6)*zb(j6,j5)*za(j3,j4)*za(j2,j5)
c     &   *za(j6,j5)/zba2(j1,j6,j5,j2)**3*t(j6,j5,j2)
c     &  + 2.D0*zb(j4,j1)*zb(j4,j6)*za(j3,j4)*za(j1,j5)
c     &   /zba2(j1,j3,j4,j2)**2*t(j3,j4,j1)
c     &  - 2.D0*zb(j4,j1)*zb(j4,j6)*za(j3,j4)*za(j1,j5)
c     &   /zba2(j1,j3,j4,j2)**2*t(j3,j4,j2)
c     &  - 2.D0*zb(j4,j1)*zb(j4,j6)*za(j3,j4)*za(j2,j5)
c     &  /zba2(j1,j3,j4,j2)**3*t(j3,j4,j1)**2
c     &  +2.D0*zb(j4,j1)*zb(j4,j6)*za(j3,j4)*za(j2,j5)
c     & /zba2(j1,j3,j4,j2)**3*t(j3,j4,j1)*t(j3,j4,j2)
c     &  - 2.D0*zb(j4,j1)*zb(j6,j5)*za(j3,j5)*za(j2,j5)
c     & /zba2(j1,j6,j5,j2)**3*t(j6,j5,j1)*t(j6,j5,j2)
c     &  + 2.D0*zb(j4,j1)*zb(j6,j5)*za(j3,j5)*za(j2,j5)
c     & /zba2(j1,j6,j5,j2)**3*t(j6,j5,j2)**2
c     &  + 2.D0*zb(j4,j2)*zb(j6,j5)*za(j3,j5)*za(j2,j5)
c     & /zba2(j1,j6,j5,j2)**2*t(j6,j5,j1)
c     &  - 2.D0*zb(j4,j2)*zb(j6,j5)*za(j3,j5)*za(j2,j5)
c     & /zba2(j1,j6,j5,j2)**2*t(j6,j5,j2)
c     &  - 1d0/2.D0*zb(j4,j6)**2*za(j3,j4)/zb(j6,j5)*zba2(j2,j4,j3,j1)
c     &  /zba2(j1,j3,j4,j2)
c     &  - 1d0/2.D0*zb(j4,j6)**2*za(j3,j4)/zb(j6,j5)*zba2(j2,j4,j3,j1)
c     & /zab2(j2,j3,j4,j1))
c     & + IDelta
c     & * ( zb(j4,j6)**2*za(j3,j4)/zb(j6,j5)/zba2(j1,j3,j4,j2)**2
c     & *s(j3,j4)*t(j3,j4,j1)
c     &  -zb(j4,j6)**2*za(j3,j4)/zb(j6,j5)
c     & /zba2(j1,j3,j4,j2)**2*s(j1,j2)*t(j3,j4,j1)
c     &  - zb(j4,j6)**2*za(j3,j4)/zb(j6,j5)
c     & /zba2(j1,j3,j4,j2)**2*s(j5,j6)*t(j3,j4,j1)
c     &  -1d0/2.D0*zb(j4,j6)**2*za(j6,j5)/zb(j3,j4)*zba2(j2,j6,j5,j1)
c     & /zba2(j1,j6,j5,j2)
c     &  - 1d0/2.D0*zb(j4,j6)**2*za(j6,j5)/zb(j3,j4)*zba2(j2,j5,j6,j1)
c     & /zba2(j1,j5,j6,j2)
c     &  - zb(j4,j6)**2*za(j6,j5)/zb(j3,j4)
c     & /zba2(j1,j5,j6,j2)**2*s(j3,j4)*t(j5,j6,j1)
c     &  -zb(j4,j6)**2*za(j6,j5)/zb(j3,j4)
c     & /zba2(j1,j5,j6,j2)**2*s(j1,j2)*t(j5,j6,j1)
c     & + zb(j4,j6)**2*za(j6,j5)/zb(j3,j4)
c     & /zba2(j1,j5,j6,j2)**2*s(j5,j6)*t(j5,j6,j1)
c     & + 2.D0*zb(j4,j6)*zb(j1,j6)*za(j3,j1)*za(j6,j5)
c     &  /zba2(j1,j5,j6,j2)**2*t(j5,j6,j1)
c     &  - 2.D0*zb(j4,j6)*zb(j1,j6)*za(j3,j1)*za(j6,j5)
c     &  /zba2(j1,j5,j6,j2)**2*t(j5,j6,j2)
c     &  - 2.D0*zb(j4,j6)*zb(j1,j6)*za(j3,j2)*za(j6,j5)
c     & /zba2(j1,j5,j6,j2)**3*t(j5,j6,j1)**2)
c     & + IDelta
c     &  * ( 2.D0*zb(j4,j6)*zb(j1,j6)*za(j3,j2)*za(j6,j5)
c     &  /zba2(j1,j5,j6,j2)**3*t(j5,j6,j1)*t(j5,j6,j2)
c     &  -zb(j4,j6)*zb(j6,j5)*za(j3,j5)*za(j6,j5)
c     &  /zba2(j1,j6,j5,j2)**2*s(j3,j4)
c     &  - zb(j4,j6)*zb(j6,j5)*za(j3,j5)*za(j6,j5)
c     &  /zba2(j1,j6,j5,j2)**2*s(j1,j2)
c     &  + zb(j4,j6)*zb(j6,j5)*za(j3,j5)*za(j6,j5)
c     &  /zba2(j1,j6,j5,j2)**2*s(j5,j6)
c     & - zb(j4,j6)*zb(j6,j5)*za(j3,j5)*za(j6,j5)
c     & /zba2(j1,j5,j6,j2)**2*s(j3,j4)
c     &  -zb(j4,j6)*zb(j6,j5)*za(j3,j5)*za(j6,j5)
c     &  /zba2(j1,j5,j6,j2)**2*s(j1,j2)
c     &  +zb(j4,j6)*zb(j6,j5)*za(j3,j5)*za(j6,j5)
c     & /zba2(j1,j5,j6,j2)**2*s(j5,j6)
c     &  +zb(j4,j6)*za(j3,j5)*zab2(j1,j3,j4,j2)/zba2(j1,j3,j4,j2)
c     &  +zb(j4,j6)*za(j3,j5)*zba2(j2,j4,j3,j1)/zab2(j2,j3,j4,j1)
c     &  +zb(j4,j6)*za(j3,j5)*zba2(j2,j6,j5,j1)*izba2(j1,j6,j5,j2)
c     &  +zb(j4,j6)*za(j3,j5)*zba2(j2,j5,j6,j1)*izba2(j1,j5,j6,j2)
c     &  - 1d0/2.D0*zb(j6,j5)*za(j3,j5)**2/za(j3,j4)*zba2(j2,j6,j5,j1)
c     & /zba2(j1,j6,j5,j2)
c     &  - 1d0/2.D0*zb(j6,j5)*za(j3,j5)**2/za(j3,j4)*zba2(j2,j5,j6,j1)
c     & /zba2(j1,j5,j6,j2))

c     & + IDelta
c     & * (-zb(j6,j5)*za(j3,j5)**2
c     & /za(j3,j4)/zba2(j1,j6,j5,j2)**2*s(j3,j4)*t(j6,j5,j2)
c     & -zb(j6,j5)*za(j3,j5)**2
c     & /za(j3,j4)/zba2(j1,j6,j5,j2)**2*s(j1,j2)*t(j6,j5,j2)
c     &  +zb(j6,j5)*za(j3,j5)**2
c     & /za(j3,j4)/zba2(j1,j6,j5,j2)**2*s(j5,j6)*t(j6,j5,j2))

c     & - 2.D0*zb(j4,j1)*zb(j4,j2)*za(j1,j5)*za(j2,j5)
c     &  /zb(j3,j4)/za(j6,j5)/zba2(j1,j3,j4,j2)**2
c     &  + 2.D0*zb(j4,j1)*za(j2,j5)
c     & /zb(j3,j4)/za(j6,j5)*zba2(j4,j3,j2,j5)
c     & /zba2(j1,j3,j4,j2)**3*t(j3,j4,j1)
c     &  -2.D0*zb(j1,j6)*zb(j2,j6)*za(j3,j1)*za(j3,j2)
c     & /zb(j6,j5)/za(j3,j4)/zab2(j2,j3,j4,j1)**2
c     & + 2.D0*zb(j1,j6)*za(j3,j2)/zb(j6,j5)/za(j3,j4)
c     & *zba2(j6,j4,j1,j3)/zab2(j2,j3,j4,j1)**3*t(j4,j3,j2)
c     &  -izb(j3,j4)/za(j6,j5)*zba2(j4,j3,j2,j5)**2/zba2(j1,j3,j4,j2)**2
c     &   -izb(j6,j5)/za(j3,j4)*zba2(j6,j4,j1,j3)**2/zab2(j2,j3,j4,j1)**2

c      write(6,*) 'new:b12zm',bcoeff(b12zm)
c      pause

c      write(6,*) 'old b134',oldfuncb134mp(j1,j2,j3,j4,j5,j6,za,zb)
c      write(6,*) 'NEW b134',funcb134mp(j1,j2,j3,j4,j5,j6,za,zb)
c      write(6,*) 'ratio',funcb134mp(j1,j2,j3,j4,j5,j6,za,zb)
c     &                  /oldfuncb134mp(j1,j2,j3,j4,j5,j6,za,zb)
c      write(6,*)
c      write(6,*) 'old b234',oldfuncb134mp(j2,j1,j4,j3,j6,j5,zb,za)
c      write(6,*) 'NEW b234',funcb134mp(j2,j1,j4,j3,j6,j5,zb,za)
c      write(6,*) 'ratio',funcb134mp(j2,j1,j4,j3,j6,j5,zb,za)
c     &                  /oldfuncb134mp(j2,j1,j4,j3,j6,j5,zb,za)
c      pause

c      write(6,*) 'old b34',oldfuncb34mp(j1,j2,j3,j4,j5,j6,za,zb)
c      write(6,*) 'NEW b34',funcb34mp(j1,j2,j3,j4,j5,j6,za,zb)
c      write(6,*) 'ratio',funcb34mp(j1,j2,j3,j4,j5,j6,za,zb)
c     &                  /oldfuncb34mp(j1,j2,j3,j4,j5,j6,za,zb)
c      write(6,*)
c      write(6,*) 'old b56',oldfuncb34mp(j2,j1,j4,j3,j6,j5,zb,za)
c      write(6,*) 'NEW b56',funcb34mp(j2,j1,j4,j3,j6,j5,zb,za)
c      write(6,*) 'ratio',funcb34mp(j2,j1,j4,j3,j6,j5,zb,za)
c     &                  /oldfuncb34mp(j2,j1,j4,j3,j6,j5,zb,za)
c      pause

      if (ggZZuse6d) then
      bcoeff(b34)=funcb34mp(j1,j2,j3,j4,j5,j6,za,zb)
      bcoeff(b56)=funcb34mp(j1,j2,j5,j6,j3,j4,za,zb)
      bcoeff(b134)=funcb134mp(j1,j2,j3,j4,j5,j6,za,zb)
      bcoeff(b234)=funcb134mp(j2,j1,j4,j3,j6,j5,zb,za)
      else
      bcoeff(b34)=oldfuncb34mp(j1,j2,j3,j4,j5,j6,za,zb)
      bcoeff(b56)=oldfuncb34mp(j1,j2,j5,j6,j3,j4,za,zb)
      bcoeff(b134)=oldfuncb134mp(j1,j2,j3,j4,j5,j6,za,zb)
      bcoeff(b234)=oldfuncb134mp(j2,j1,j4,j3,j6,j5,zb,za)
      endif
      bcoeff(b12)=-bcoeff(b34)-bcoeff(b56)
     & -bcoeff(b134)-bcoeff(b234)

      bcoeff(rat)=
     & ((zb(j1,j4)*za(j1,j5))**2*za(j4,j3)
     & /(za(j6,j5)*(s(j1,j4)+s(j1,j3)))
     & +(zb(j2,j6)*za(j2,j3))**2*zb(j4,j3)
     & /(zb(j6,j5)*(s(j2,j4)+s(j2,j3)))
     & +(zb(j1,j6)*za(j1,j3))**2*za(j6,j5)
     & /(za(j4,j3)*(s(j1,j6)+s(j1,j5)))
     & +(zb(j2,j4)*za(j2,j5))**2*zb(j6,j5)
     & /(zb(j4,j3)*(s(j2,j6)+s(j2,j5)))
     & -zab2(j5,j2,j3,j4)**2/zb(j4,j3)/za(j6,j5)
     & -zab2(j3,j1,j4,j6)**2/za(j4,j3)/zb(j6,j5)
     & -2d0*zb(j4,j6)*za(j3,j5))/zab2(j2,j3,j4,j1)**2
     & -IDelta*zab2(j1,j3,j4,j2)/zab2(j2,j3,j4,j1)*(
     & +4d0*zb(j4,j6)*za(j3,j5)
     & +(s(j1,j2)-s(j3,j4)-s(j5,j6))
     & *(zb(j4,j6)**2/zb(j4,j3)/zb(j6,j5)
     &  +za(j3,j5)**2/za(j4,j3)/za(j6,j5)))

      endif

      return
      end

c--- This form has been massaged to remove the explicit <2|3+4|1]^3 poles
      double complex function funcb34mp(j1,j2,j3,j4,j5,j6,za,zb)
      implicit none
      include 'constants.f'
      include 'zprods_decl.f'
      include 'sprods_com.f'
      integer j1,j2,j3,j4,j5,j6
      double precision t,sbar134,sbar234,del,s134,s234
      double complex zab2
      double precision IDelta
c---statement functions
      zab2(j1,j2,j3,j4)=za(j1,j2)*zb(j2,j4)+za(j1,j3)*zb(j3,j4)
      t(j1,j2,j3)=s(j1,j2)+s(j2,j3)+s(j1,j3)
      del(j1,j2,j3,j4,j5,j6)=s(j1,j2)-s(j3,j4)-s(j5,j6)
c---end statement functions

      IDelta=s(j1,j2)**2+s(j3,j4)**2+s(j5,j6)**2
     & -2d0*s(j1,j2)*s(j3,j4)
     & -2d0*s(j3,j4)*s(j5,j6)
     & -2d0*s(j1,j2)*s(j5,j6)

      IDelta=1d0/IDelta
      s134=t(j1,j3,j4)
      s234=t(j2,j3,j4)
      sbar134=s134-s(j3,j4)
      sbar234=s234-s(j3,j4)

      funcb34mp=
     & -za(j4,j3)*za(j1,j5)**2*zb(j1,j4)**2*s(j3,j4)
     & /sbar134**2/za(j6,j5)/zab2(j2,j3,j4,j1)**2

     & -zb(j4,j3)*zb(j2,j6)**2*za(j2,j3)**2*s(j3,j4)
     & /sbar234**2/zb(j6,j5)/zab2(j2,j3,j4,j1)**2

     & -2d0*za(j3,j2)*zb(j1,j4)*s(j3,j4)
     & /s(j5,j6)/zab2(j2,j3,j4,j1)**3/sbar134/sbar234
     & *(s(j3,j4)*(za(j2,j5)*zb(j1,j6)*zab2(j1,j3,j4,j2)
     &            -za(j1,j5)*zb(j2,j6)*zab2(j2,j3,j4,j1))
     &  -za(j1,j5)*zb(j5,j6)*zab2(j5,j3,j4,j2)*zab2(j2,j3,j4,j1)
     &  +za(j5,j6)*zb(j2,j6)*zab2(j1,j3,j4,j6)*zab2(j2,j3,j4,j1))

     & +6d0*IDelta**2
     & *zab2(j3,j1,j2,j4)*zab2(j5,j1,j2,j6)*del(j5,j6,j1,j2,j3,j4)
     & *zab2(j1,j3,j4,j2)/zab2(j2,j3,j4,j1)

     & + IDelta *(
     & +2d0*(s134-s234)/s(j5,j6)*zab2(j1,j3,j4,j2)/zab2(j2,j3,j4,j1)*(
     &  +s234*zb(j1,j4)*zb(j1,j6)*za(j2,j3)*za(j2,j5)
     &       /zab2(j2,j3,j4,j1)**2
     &  +s234*zb(j1,j4)*za(j1,j5)*zab2(j3,j2,j5,j6)
     &       /zab2(j1,j3,j4,j2)/zab2(j2,j3,j4,j1)
     &  +zb(j4,j6)*za(j2,j5)*zab2(j3,j2,j4,j1)/zab2(j2,j3,j4,j1)
     &  -s134*za(j2,j3)*za(j2,j5)*zb(j1,j4)*zb(j1,j6)
     &       /zab2(j2,j3,j4,j1)**2
     &  -s134*za(j2,j3)*zb(j2,j6)*zab2(j5,j1,j6,j4)
     &       /zab2(j1,j3,j4,j2)/zab2(j2,j3,j4,j1)
     &  -za(j3,j5)*zb(j1,j6)*zab2(j2,j1,j3,j4)/zab2(j2,j3,j4,j1) )

     & -8d0*za(j1,j5)*za(j2,j5)*zb(j1,j4)**2*za(j3,j4)
     &  /za(j5,j6)/zab2(j2,j3,j4,j1)**2*zab2(j1,j3,j4,j2)

     & -8d0*za(j2,j3)**2*zb(j3,j4)*zb(j2,j6)*zb(j1,j6)
     &  /zb(j5,j6)/zab2(j2,j3,j4,j1)**2*zab2(j1,j3,j4,j2)

     & +2d0*za(j3,j4)*za(j1,j5)*zb(j4,j1)*zb(j4,j6)*(s234-s134)
     & /zab2(j2,j3,j4,j1)**2
     & +2d0*za(j3,j2)*za(j3,j5)*zb(j3,j4)*zb(j2,j6)*(s234-s134)
     & /zab2(j2,j3,j4,j1)**2
     & +2d0*za(j3,j5)*zb(j4,j6)*s(j3,j4)*del(j3,j4,j1,j2,j5,j6)
     & /zab2(j2,j3,j4,j1)**2
     & -za(j3,j5)**2*zb(j3,j4)*del(j3,j4,j1,j2,j5,j6)*s234
     & /za(j6,j5)/zab2(j2,j3,j4,j1)**2
     & -za(j3,j4)*zb(j4,j6)**2*del(j3,j4,j1,j2,j5,j6)*s134
     & /zb(j6,j5)/zab2(j2,j3,j4,j1)**2

     & +za(j3,j5)**2*zb(j3,j4)*zab2(j1,j3,j4,j2)
     & /za(j6,j5)/zab2(j2,j3,j4,j1)

     & +za(j3,j4)*zb(j4,j6)**2*zab2(j1,j3,j4,j2)
     & /zb(j6,j5)/zab2(j2,j3,j4,j1)

     & -2d0*za(j3,j5)*zb(j4,j6)*zb(j6,j5)*zab2(j1,j3,j4,j2)
     & /zb(j6,j5)/zab2(j2,j3,j4,j1)
     & )

      return
      end

c--- This form is extracted from BDK
      double complex function oldfuncb34mp(j1,j2,j3,j4,j5,j6,za,zb)
      implicit none
      include 'constants.f'
      include 'zprods_decl.f'
      include 'sprods_com.f'
      integer j1,j2,j3,j4,j5,j6
      double precision t,sbar134,sbar234,del,s134,s234
      double complex zab2
      double precision IDelta
c---statement functions
      zab2(j1,j2,j3,j4)=za(j1,j2)*zb(j2,j4)+za(j1,j3)*zb(j3,j4)
      t(j1,j2,j3)=s(j1,j2)+s(j2,j3)+s(j1,j3)
      del(j1,j2,j3,j4,j5,j6)=s(j1,j2)-s(j3,j4)-s(j5,j6)
c---end statement functions

      IDelta=s(j1,j2)**2+s(j3,j4)**2+s(j5,j6)**2
     & -2d0*s(j1,j2)*s(j3,j4)
     & -2d0*s(j3,j4)*s(j5,j6)
     & -2d0*s(j1,j2)*s(j5,j6)

      IDelta=1d0/IDelta
      s134=t(j1,j3,j4)
      s234=t(j2,j3,j4)
      sbar134=s134-s(j3,j4)
      sbar234=s234-s(j3,j4)

      oldfuncb34mp=
     & -za(j4,j3)*za(j1,j5)**2*zb(j1,j4)**2*s(j3,j4)
     & /sbar134**2/za(j6,j5)/zab2(j2,j3,j4,j1)**2

     & -zb(j4,j3)*zb(j2,j6)**2*za(j2,j3)**2*s(j3,j4)
     & /sbar234**2/zb(j6,j5)/zab2(j2,j3,j4,j1)**2

     & -2d0*za(j1,j5)*za(j3,j2)*zb(j4,j1)*zab2(j5,j3,j4,j1)*s(j3,j4)
     & /sbar134/za(j6,j5)/zab2(j2,j3,j4,j1)**3

     & -2d0*za(j3,j2)*zb(j2,j6)*zb(j1,j4)*zab2(j2,j3,j4,j6)*s(j3,j4)
     &  /sbar234/zb(j6,j5)/zab2(j2,j3,j4,j1)**3

     & -2d0*za(j1,j5)*za(j2,j5)*zb(j1,j4)**2*za(j4,j3)
     &  /za(j6,j5)/zab2(j2,j3,j4,j1)**3

     & -2d0*za(j2,j3)**2*zb(j4,j3)*zb(j2,j6)*zb(j1,j6)
     &  /zb(j6,j5)/zab2(j2,j3,j4,j1)**3

     & +6d0*IDelta**2
     & *zab2(j3,j1,j2,j4)*zab2(j5,j1,j2,j6)*del(j5,j6,j1,j2,j3,j4)
     & *zab2(j1,j3,j4,j2)/zab2(j2,j3,j4,j1)

     & + IDelta *(
     & +2d0*za(j3,j5)*za(j2,j5)*zb(j4,j1)*zb(j6,j5)*s(j3,j4)*(s234-s134)
     & /zab2(j2,j3,j4,j1)**3
     & -2d0*za(j3,j4)*za(j2,j5)*zb(j4,j1)*zb(j4,j6)*s134*(s234-s134)
     & /zab2(j2,j3,j4,j1)**3
     & +2d0*za(j3,j2)*za(j6,j5)*zb(j4,j6)*zb(j1,j6)*s(j3,j4)*(s234-s134)
     & /zab2(j2,j3,j4,j1)**3
     & -2d0*za(j3,j2)*za(j3,j5)*zb(j3,j4)*zb(j1,j6)*(s234-s134)*s234
     & /zab2(j2,j3,j4,j1)**3

     & +2d0*za(j3,j4)*za(j1,j5)*zb(j4,j1)*zb(j4,j6)*(s234-s134)
     & /zab2(j2,j3,j4,j1)**2
     & +2d0*za(j3,j2)*za(j3,j5)*zb(j3,j4)*zb(j2,j6)*(s234-s134)
     & /zab2(j2,j3,j4,j1)**2
     & +2d0*za(j3,j5)*zb(j4,j6)*s(j3,j4)*del(j3,j4,j1,j2,j5,j6)
     & /zab2(j2,j3,j4,j1)**2
     & -za(j3,j5)**2*zb(j3,j4)*del(j3,j4,j1,j2,j5,j6)*s234
     & /za(j6,j5)/zab2(j2,j3,j4,j1)**2
     & -za(j3,j4)*zb(j4,j6)**2*del(j3,j4,j1,j2,j5,j6)*s134
     & /zb(j6,j5)/zab2(j2,j3,j4,j1)**2

     & +za(j3,j5)**2*zb(j3,j4)*zab2(j1,j3,j4,j2)
     & /za(j6,j5)/zab2(j2,j3,j4,j1)

     & +za(j3,j4)*zb(j4,j6)**2*zab2(j1,j3,j4,j2)
     & /zb(j6,j5)/zab2(j2,j3,j4,j1)

     & -2d0*za(j3,j5)*zb(j4,j6)*zb(j6,j5)*zab2(j1,j3,j4,j2)
     & /zb(j6,j5)/zab2(j2,j3,j4,j1)
     & )

      return
      end

c--- This form has been massaged to remove the explicit <2|3+4|1]^3 poles
      double complex function funcb134mp(j1,j2,j3,j4,j5,j6,za,zb)
      implicit none
      include 'constants.f'
      include 'zprods_decl.f'
      include 'sprods_com.f'
      integer j1,j2,j4,j3,j6,j5
      double complex zab2,zbb22
c---statement functions
      zab2(j1,j2,j3,j4)=za(j1,j2)*zb(j2,j4)+za(j1,j3)*zb(j3,j4)
      zbb22(j1,j2,j3,j4,j5,j6)
     & =zb(j1,j2)*zab2(j2,j4,j5,j6)+zb(j1,j3)*zab2(j3,j4,j5,j6)
c---end statement functions

      funcb134mp=

     &+ zab2(j1,j3,j4,j2)/zab2(j2,j3,j4,j1)**3 * (
     & + 2d0*za(j2,j5)**2*zb(j1,j4)/za(j1,j2)/za(j5,j6)/zb(j1,j2)
     &   /zb(j3,j4)*zbb22(j4,j5,j6,j3,j4,j1))

     &+ 1d0/zab2(j2,j3,j4,j1)**2 * (
     & + s(j3,j4)*za(j1,j5)**2*za(j3,j4)*zb(j1,j4)**2
     &   /za(j5,j6)/zab2(j1,j3,j4,j1)**2
     & + s(j5,j6)*za(j2,j5)**2*zb(j2,j4)**2*zb(j5,j6)
     &   /zb(j3,j4)/zab2(j2,j5,j6,j2)**2
     & - 2d0*s(j3,j4)*za(j1,j5)**2*za(j2,j3)*zb(j1,j4)
     &   /za(j1,j2)/za(j5,j6)/zab2(j1,j3,j4,j1)
     & - 2d0*s(j5,j6)*za(j2,j5)*zb(j1,j6)*zb(j2,j4)**2
     &   /zb(j1,j2)/zb(j3,j4)/zab2(j2,j5,j6,j2)
     & + 2d0*za(j1,j5)*za(j2,j5)*zb(j1,j4)*zbb22(j2,j3,j4,j5,j6,j4)
     &   /za(j1,j2)/zb(j1,j2)/zb(j3,j4)/za(j5,j6)
     & - 1d0/za(j5,j6)/zb(j3,j4)*zab2(j5,j2,j3,j4)**2)

      return
      end

c--- This form is extracted from BDK
      double complex function oldfuncb134mp(j1,j2,j3,j4,j5,j6,za,zb)
      implicit none
      include 'constants.f'
      include 'zprods_decl.f'
      include 'sprods_com.f'
      integer j1,j2,j4,j3,j6,j5
      double precision t
      double complex zab2
c---statement functions
      zab2(j1,j2,j3,j4)=za(j1,j2)*zb(j2,j4)+za(j1,j3)*zb(j3,j4)
      t(j1,j2,j4)=s(j1,j2)+s(j2,j4)+s(j1,j4)
c---end statement functions

      oldfuncb134mp=
     &  +za(j3,j4)*za(j1,j5)**2*zb(j1,j4)**2*s(j3,j4)
     & /((s(j1,j4)+s(j1,j3))**2*za(j5,j6)*zab2(j2,j3,j4,j1)**2)

     & +za(j2,j5)**2*zb(j2,j4)**2*zb(j5,j6)*s(j5,j6)
     & /((s(j2,j6)+s(j2,j5))**2*zb(j3,j4)*zab2(j2,j5,j6,j1)**2)

     & -2d0*za(j2,j3)*zb(j1,j4)*za(j1,j5)*zab2(j5,j3,j4,j1)*s(j3,j4)
     & /((s(j1,j4)+s(j1,j3))*za(j5,j6)*zab2(j2,j3,j4,j1)**3)

     & +2d0*za(j2,j5)*zb(j2,j4)*zb(j1,j6)*zab2(j2,j5,j6,j4)*s(j5,j6)
     & /((s(j2,j6)+s(j2,j5))*zb(j3,j4)*zab2(j2,j5,j6,j1)**3)

     & +(2d0*zb(j1,j4)*za(j2,j5)*(
     & +zb(j1,j4)*za(j1,j5)*t(j2,j4,j3)
     & -zb(j2,j4)*za(j2,j5)*s(j5,j6)
     & +zb(j2,j1)*za(j1,j5)*zb(j4,j3)*za(j3,j2)
     & +t(j1,j4,j3)*zab2(j5,j2,j3,j4))
     & +zab2(j5,j2,j3,j4)**2*zab2(j2,j3,j4,j1))
     & /zb(j3,j4)/za(j6,j5)/zab2(j2,j3,j4,j1)**3

      return
      end

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text:
      program unstruc11

C     Check a scc which is path coverable but not a cycle

      k = 0
      i = 0
      j = 0
      l = 0

 100  continue
      if(x.gt.0.) go to 300
      l = l + 4
 200  continue
      k = k + 3
      i = i + 1
      if(y.gt.0.) go to 100
      j = j + 2
      go to 200

 300  continue
      print *, i, j, k, l

      end