module forcing_vaues
    real*8 F_1, F_2
    real*8 D_1, D_2
end module

subroutine rhs_velocity

    use m_openmpi
    use m_io
    use m_parameters
    use m_fields
    use m_work
    use x_fftw
    use m_les

    implicit none

    integer :: i, j, k, n, nx3
    real*8  :: t1(0:6), rtmp, wnum2, rnx3, epsilon1, epsilon2

    ! The IFFT of velocities has been done earlier in rhs_scalars
    ! the velocities were kept in wrk1...wrk3, intact.

!!$  ! putting the velocity field in the wrk array
!!$  wrk(:,:,:,1:3) = fields(:,:,:,1:3)
!!$  ! performing IFFT to convert them to the X-space
!!$  call xFFT3d(-1,1)
!!$  call xFFT3d(-1,2)
!!$  call xFFT3d(-1,3)

    ! right now the velocities should be in x space

    ! getting the Courant number (on the master process only)
    wrk(:, :, :, 4) = abs(wrk(:, :, :, 1)) + abs(wrk(:, :, :, 2)) + abs(wrk(:, :, :, 3))
    rtmp = maxval(wrk(1:nx, :, :, 4))
    call MPI_REDUCE(rtmp, courant, 1, MPI_REAL8, MPI_MAX, 0, MPI_COMM_TASK, mpi_err)
    if (variable_dt) then
        count = 1
        call MPI_BCAST(courant, count, MPI_REAL8, 0, MPI_COMM_TASK, mpi_err)
    end if
    courant = courant*dt/dx

!-------------------------------------------------------------------------

!--------------------------------------------------------------------------------
!  Calculating the right-hand side for the velocities
!
!  There are two options available: the standard 2/3 rule (dealias=0) and
!  combination of phase shift and truncation (dealias=1).  The latter retains
!  more modes but requires more calculations thus slowing down the simulation.
!  These are treated separately in two different "if" blocks.  This is done in
!  order not to complicate the logic.  Also this way both blocks can be
!  optimized separately.
!--------------------------------------------------------------------------------

    two_thirds_rule: if (dealias .eq. 0) then

        ! getting all 6 products of velocities
        do k = 1, nz
            do j = 1, ny
                do i = 1, nx
                    t1(1) = wrk(i, j, k, 1)*wrk(i, j, k, 1)
                    t1(2) = wrk(i, j, k, 1)*wrk(i, j, k, 2)
                    t1(3) = wrk(i, j, k, 1)*wrk(i, j, k, 3)
                    t1(4) = wrk(i, j, k, 2)*wrk(i, j, k, 2)
                    t1(5) = wrk(i, j, k, 2)*wrk(i, j, k, 3)
                    t1(6) = wrk(i, j, k, 3)*wrk(i, j, k, 3)

                    do n = 1, 6
                        wrk(i, j, k, n) = t1(n)
                    end do

                end do
            end do
        end do

        ! converting the products to the Fourier space
        do n = 1, 6
            call xFFT3d(1, n)
        end do

        ! Building the RHS.
        ! First, put into wrk arrays the convective terms (that will be multiplyed by "i"
        ! later) and the factor that corresponds to the diffusion

        ! Do not forget that in Fourier space the indicies are (ix, iz, iy)
        do k = 1, nz
            do j = 1, ny
                do i = 1, nx + 2
                    t1(1) = -(akx(i)*wrk(i, j, k, 1) + aky(k)*wrk(i, j, k, 2) + akz(j)*wrk(i, j, k, 3))
                    t1(2) = -(akx(i)*wrk(i, j, k, 2) + aky(k)*wrk(i, j, k, 4) + akz(j)*wrk(i, j, k, 5))
                    t1(3) = -(akx(i)*wrk(i, j, k, 3) + aky(k)*wrk(i, j, k, 5) + akz(j)*wrk(i, j, k, 6))

                    ! this looks like a diffusion term
                    t1(4) = -nu*(akx(i)**2 + aky(k)**2 + akz(j)**2)

                    do n = 1, 4
                        wrk(i, j, k, n) = t1(n)
                    end do

                end do
            end do
        end do

        ! call a subroutine to handle both viscous compensation and regular compensation
        ! if the viscous compensation parameter was not used
        if (PERT == 1) then
            call update_forcing_viscous_compensation(PERTamp1, PERTamp2)
        end if

        ! now take the actual fields from fields(:,:,:,:) and calculate the RHSs

        ! at this moment the contains of wrk(:,:,:,1:3) are the convective terms in the RHS
        ! which are not yet multiplied by "i"
        ! wrk(:,:,:,4) contains the Laplace operator in Fourier space.  To get the diffusion term
        ! we need to take wrk(:,:,:,4) and multiply it by the velocity

        t1(6) = real(kmax, 8)

        do k = 1, nz
            do j = 1, ny
                do i = 1, nx+ 1, 2
                    ! If the dealiasing option is 2/3-rule (dealias=0) then we retain the modes
                    ! inside the cube described by $| k_i | \leq  k_{max}$, $i=1,2,3$.
                    ! The rest of the modes is purged

                    if (ialias(i, j, k) .gt. 0) then ! run for dealias=0
                        ! setting the Fourier components to zero
                        wrk(i, j, k, 1:3) = zip
                        wrk(i + 1, j, k, 1:3) = zip

                    else ! run for dealias=1

                        ! RHS for u, v and w
                        do n = 1, 3

                            ! this is what the code was like before MGM
                            ! https://github.com/Fluid-Dynamics-Group/hit3d/blob/92db6d9e864bdbf2248b025a71f315cba21d72fd/rhs_velocity.f90#L129-L135
                            if (skip_diffusion == 1) then ! skip diffusion calculation
                                ! ===========================MGM-Forcing=====================
                                ! Validation by inviscid flow
                                rtmp = -wrk(i + 1, j, k, n) + fcomp(i  ,j,k,n)
                                wrk(i + 1, j, k, n) = wrk(i, j, k, n) + fcomp(i+1,j,k,n)
                                wrk(i, j, k, n) = rtmp
                            else ! dont skip the diffusion calculation
                                ! ==========================================================
                                ! BROOKS: original hit3d code with MGM forcing is here (with diffusion, i think)
                                ! ==========================================================

                                ! the diffusion terms from wrk(4) here have already been calculated in real space
                                ! if there is viscous compensation - but it simplifies the code to not do anything
                                ! about it here

                                ! taking the convective term, multiply it by "i"
                                ! (see how it's done in x_fftw.f90)
                                ! and adding the diffusion term
                                rtmp =           - wrk(i+1,j,k,n) + wrk(i  ,j,k,4) * fields(i  ,j,k,n) + fcomp(i  ,j,k,n)
                                wrk(i+1,j,k,n) =   wrk(i  ,j,k,n) + wrk(i+1,j,k,4) * fields(i+1,j,k,n) + fcomp(i+1,j,k,n)
                                wrk(i, j, k, n) = rtmp

                                ! dot this wrk variable with u - can try doing truncation here based on rtmp

                            end if

                        end do

                    end if

                end do
            end do
        end do

    end if two_thirds_rule

!--------------------------------------------------------------------------------
!  The second option (dealias=1).  All pairwise products of velocities are
!  dealiased using one phase shift of (dx/2,dy/2,dz/2).
!--------------------------------------------------------------------------------
    phase_shifting: if (dealias .eq. 1) then

        ! work parameters
        wrk(:, :, :, 0) = zip

        ! getting all 6 products of velocities
        do k = 1, nz
            do j = 1, ny
                do i = 1, nx

                    t1(1) = wrk(i, j, k, 1)*wrk(i, j, k, 1)
                    t1(2) = wrk(i, j, k, 1)*wrk(i, j, k, 2)
                    t1(3) = wrk(i, j, k, 1)*wrk(i, j, k, 3)
                    t1(4) = wrk(i, j, k, 2)*wrk(i, j, k, 2)
                    t1(5) = wrk(i, j, k, 2)*wrk(i, j, k, 3)
                    t1(6) = wrk(i, j, k, 3)*wrk(i, j, k, 3)
                    do n = 1, 6
                        wrk(i, j, k, n) = t1(n)
                    end do
                end do
            end do
        end do

        ! converting the products to the Fourier space
        do n = 1, 6
            call xFFT3d(1, n)
        end do

        ! Building the RHS.
        ! First, put into wrk arrays the convectove terms (that will be multiplyed by "i"
        ! later) and the factor that corresponds to the diffusion

        ! Do not forget that in Fourier space the indicies are (ix, iz, iy)
        do k = 1, nz
            do j = 1, ny
                do i = 1, nx + 2

                    t1(1) = -(akx(i)*wrk(i, j, k, 1) + aky(k)*wrk(i, j, k, 2) + akz(j)*wrk(i, j, k, 3))
                    t1(2) = -(akx(i)*wrk(i, j, k, 2) + aky(k)*wrk(i, j, k, 4) + akz(j)*wrk(i, j, k, 5))
                    t1(3) = -(akx(i)*wrk(i, j, k, 3) + aky(k)*wrk(i, j, k, 5) + akz(j)*wrk(i, j, k, 6))

                    t1(4) = -nu*(akx(i)**2 + aky(k)**2 + akz(j)**2)

                    do n = 1, 4
                        wrk(i, j, k, n) = t1(n)
                    end do
                end do
            end do
        end do

        ! now use the actual fields from fields(:,:,:,:) to calculate the RHSs

        ! at this moment the contains of wrk(:,:,:,1:3) are the convective terms in the RHS
        ! which are not yet multiplied by "i"

        ! wrk(:,:,:,4) contains the Laplace operator in Fourier space.  To get the diffusion term
        ! we need to take wrk(:,:,:,4) and multiply it by the velocity

        do k = 1, nz
            do j = 1, ny
                do i = 1, nx + 1, 2

                    ! If the dealiasing option is (dealias=1) then we retain the modes
                    ! for which no more than one component of the k-vector is larger than nx/3.
                    ! The rest of the modes is purged.

                    if (ialias(i, j, k) .gt. 1) then ! run for dealias = 1
                        ! setting the Fourier components to zero
                        wrk(i, j, k, 1:3) = zip
                        wrk(i + 1, j, k, 1:3) = zip
                    else ! run for dealias = 0
                        if (skip_diffusion == 1) then ! we _ARE NOT_ doing diffusion calculations
                            ! RHS for u, v and w
                            do n = 1, 3
                                rtmp = -0.5d0*wrk(i + 1, j, k, n) !+ wrk(i  ,j,k,4) * fields(i  ,j,k,n)
                                wrk(i + 1, j, k, n) = 0.5d0*wrk(i, j, k, n) !+ wrk(i+1,j,k,4) * fields(i+1,j,k,n)
                                wrk(i, j, k, n) = rtmp
                            end do
                        else ! we _ARE_ doing diffusion calculations
                            ! RHS for u, v and w
                            do n = 1, 3
                                ! taking the HALF of the convective term, multiply it by "i"
                                ! and adding the diffusion term
                                rtmp = -0.5d0*wrk(i + 1, j, k, n) + wrk(i, j, k, 4)*fields(i, j, k, n)
                                wrk(i + 1, j, k, n) = 0.5d0*wrk(i, j, k, n) + wrk(i + 1, j, k, 4)*fields(i + 1, j, k, n)
                                wrk(i, j, k, n) = rtmp
                            end do

                        end if
                    end if

                end do
            end do
        end do

!--------------------------------------------------------------------------------
!  Second part of the phase shifting technique
!--------------------------------------------------------------------------------

        ! since wrk1...3 are taken by parts of RHS constructed earlier, we can use
        ! only wrk0 and wrk4...6.

        do k = 1, nz
            do j = 1, ny
                do i = 1, nx + 1, 2

                    ! computing sines and cosines for the phase shift of dx/2,dy/2,dz/2
                    ! and putting them into wrk0
                    wrk(i, j, k, 0) = cos(half*(akx(i) + aky(k) + akz(j))*dx)
                    wrk(i + 1, j, k, 0) = sin(half*(akx(i + 1) + aky(k) + akz(j))*dx)

                    ! wrk4 will have phase-shifted u
                    wrk(i, j, k, 4) = fields(i, j, k, 1)*wrk(i, j, k, 0) - fields(i + 1, j, k, 1)*wrk(i + 1, j, k, 0)
                    wrk(i + 1, j, k, 4) = fields(i + 1, j, k, 1)*wrk(i, j, k, 0) + fields(i, j, k, 1)*wrk(i + 1, j, k, 0)

                    ! wrk5 will have phase-shifted v
                    wrk(i, j, k, 5) = fields(i, j, k, 2)*wrk(i, j, k, 0) - fields(i + 1, j, k, 2)*wrk(i + 1, j, k, 0)
                    wrk(i + 1, j, k, 5) = fields(i + 1, j, k, 2)*wrk(i, j, k, 0) + fields(i, j, k, 2)*wrk(i + 1, j, k, 0)

                end do
            end do
        end do

        ! transforming u+ and v+ into X-space
        call xFFT3d(-1, 4)
        call xFFT3d(-1, 5)

        ! now wrk4 and wrk5 contain u+ and v+

        ! getting (u+)*(u+) in real space, converting it to Fourier space,
        ! phase shifting back and adding -0.5*(the results)  to the RHS for u
        wrk(:, :, :, 6) = wrk(:, :, :, 4)**2
        call xFFT3d(1, 6)
        do k = 1, nz
            do j = 1, ny
                do i = 1, nx + 1, 2
                    rtmp = wrk(i, j, k, 6)*wrk(i, j, k, 0) + wrk(i + 1, j, k, 6)*wrk(i + 1, j, k, 0)
                    wrk(i + 1, j, k, 6) = wrk(i + 1, j, k, 6)*wrk(i, j, k, 0) - wrk(i, j, k, 6)*wrk(i + 1, j, k, 0)
                    wrk(i, j, k, 6) = rtmp
                end do
            end do
        end do
        do k = 1, nz
            do j = 1, ny
                do i = 1, nx + 1, 2
                    wrk(i, j, k, 1) = wrk(i, j, k, 1) + 0.5d0*akx(i + 1)*wrk(i + 1, j, k, 6)
                    wrk(i + 1, j, k, 1) = wrk(i + 1, j, k, 1) - 0.5d0*akx(i)*wrk(i, j, k, 6)
                end do
            end do
        end do

        ! getting (u+)*(v+) in real space, converting it to Fourier space,
        ! phase shifting back and adding -0.5*(the results)  to the RHSs for u and v
        wrk(:, :, :, 6) = wrk(:, :, :, 4)*wrk(:, :, :, 5)
        call xFFT3d(1, 6)
        do k = 1, nz
            do j = 1, ny
                do i = 1, nx + 1, 2
                    rtmp = wrk(i, j, k, 6)*wrk(i, j, k, 0) + wrk(i + 1, j, k, 6)*wrk(i + 1, j, k, 0)
                    wrk(i + 1, j, k, 6) = wrk(i + 1, j, k, 6)*wrk(i, j, k, 0) - wrk(i, j, k, 6)*wrk(i + 1, j, k, 0)
                    wrk(i, j, k, 6) = rtmp
                end do
            end do
        end do
        do k = 1, nz
            do j = 1, ny
                do i = 1, nx + 1, 2
                    wrk(i, j, k, 1) = wrk(i, j, k, 1) + 0.5d0*aky(k)*wrk(i + 1, j, k, 6)
                    wrk(i + 1, j, k, 1) = wrk(i + 1, j, k, 1) - 0.5d0*aky(k)*wrk(i, j, k, 6)

                    wrk(i, j, k, 2) = wrk(i, j, k, 2) + 0.5d0*akx(i + 1)*wrk(i + 1, j, k, 6)
                    wrk(i + 1, j, k, 2) = wrk(i + 1, j, k, 2) - 0.5d0*akx(i)*wrk(i, j, k, 6)
                end do
            end do
        end do

        ! getting (v+)*(v+) in real space, converting it to Fourier space,
        ! phase shifting back and adding -0.5*(the results)  to the RHS for v
        wrk(:, :, :, 6) = wrk(:, :, :, 5)**2
        call xFFT3d(1, 6)
        do k = 1, nz
            do j = 1, ny
                do i = 1, nx + 1, 2
                    rtmp = wrk(i, j, k, 6)*wrk(i, j, k, 0) + wrk(i + 1, j, k, 6)*wrk(i + 1, j, k, 0)
                    wrk(i + 1, j, k, 6) = wrk(i + 1, j, k, 6)*wrk(i, j, k, 0) - wrk(i, j, k, 6)*wrk(i + 1, j, k, 0)
                    wrk(i, j, k, 6) = rtmp
                end do
            end do
        end do
        do k = 1, nz
            do j = 1, ny
                do i = 1, nx + 1, 2
                    wrk(i, j, k, 2) = wrk(i, j, k, 2) + 0.5d0*aky(k)*wrk(i + 1, j, k, 6)
                    wrk(i + 1, j, k, 2) = wrk(i + 1, j, k, 2) - 0.5d0*aky(k)*wrk(i, j, k, 6)
                end do
            end do
        end do

        ! now get the (w+) in wrk6
        do k = 1, nz
            do j = 1, ny
                do i = 1, nx + 1, 2
                    ! wrk6 will have phase-shifted w
                    wrk(i, j, k, 6) = fields(i, j, k, 3)*wrk(i, j, k, 0) - fields(i + 1, j, k, 3)*wrk(i + 1, j, k, 0)
                    wrk(i + 1, j, k, 6) = fields(i + 1, j, k, 3)*wrk(i, j, k, 0) + fields(i, j, k, 3)*wrk(i + 1, j, k, 0)
                end do
            end do
        end do
        ! transforming w+ into X-space
        call xFFT3d(-1, 6)

        ! at this point wrk4..6 contain (u+), (v+) and (w+) in real space.
        ! the combinations that we have not dealt with are: uw, vw and ww.
        ! we'll deal with all three of them at once.

        ! first get all three of these in wrk4...6 and
        wrk(:, :, :, 4) = wrk(:, :, :, 4)*wrk(:, :, :, 6)
        wrk(:, :, :, 5) = wrk(:, :, :, 5)*wrk(:, :, :, 6)
        wrk(:, :, :, 6) = wrk(:, :, :, 6)**2

        ! transform them into Fourier space
        call xFFT3d(1, 4)
        call xFFT3d(1, 5)
        call xFFT3d(1, 6)

        ! phase shift back to origianl grid and add to corresponding RHSs
        do n = 4, 6
            do k = 1, nz
                do j = 1, ny
                    do i = 1, nx + 1, 2
                        rtmp = wrk(i, j, k, n)*wrk(i, j, k, 0) + wrk(i + 1, j, k, n)*wrk(i + 1, j, k, 0)
                        wrk(i + 1, j, k, n) = wrk(i + 1, j, k, n)*wrk(i, j, k, 0) - wrk(i, j, k, n)*wrk(i + 1, j, k, 0)
                        wrk(i, j, k, n) = rtmp
                    end do
                end do
            end do
        end do

        ! adding to corresponding RHSs
        do k = 1, nz
            do j = 1, ny
                do i = 1, nx + 1, 2

                    ! If the dealiasing option is (dealias=1) then we retain the modes
                    ! for which no more than one component of the k-vector is larger than nx/3.
                    ! The rest of the modes is purged.

                    if (ialias(i, j, k) .lt. 2) then

                        wrk(i, j, k, 1) = wrk(i, j, k, 1) + 0.5d0*akz(j)*wrk(i + 1, j, k, 4)
                        wrk(i + 1, j, k, 1) = wrk(i + 1, j, k, 1) - 0.5d0*akz(j)*wrk(i, j, k, 4)

                        wrk(i, j, k, 2) = wrk(i, j, k, 2) + 0.5d0*akz(j)*wrk(i + 1, j, k, 5)
                        wrk(i + 1, j, k, 2) = wrk(i + 1, j, k, 2) - 0.5d0*akz(j)*wrk(i, j, k, 5)

                        wrk(i, j, k, 3) = wrk(i, j, k, 3) + 0.5d0* &
                                          (akx(i + 1)*wrk(i + 1, j, k, 4) + aky(k)*wrk(i + 1, j, k, 5) + akz(j)*wrk(i + 1, j, k, 6))
                        wrk(i + 1, j, k, 3) = wrk(i + 1, j, k, 3) - 0.5d0* &
                                              (akx(i)*wrk(i, j, k, 4) + aky(k)*wrk(i, j, k, 5) + akz(j)*wrk(i, j, k, 6))

                    else
                        wrk(i:i + 1, j, k, 1) = zip
                        wrk(i:i + 1, j, k, 2) = zip
                        wrk(i:i + 1, j, k, 3) = zip
                    end if

                end do
            end do
        end do

    end if phase_shifting

    ! if performing large eddy simulations, call LES subroutine to augment
    ! the right hand side for velocioties
    les_active: if (les) then
        call les_rhs_velocity
    end if les_active

    return
end subroutine rhs_velocity

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

subroutine test_rhs_velocity

    use m_openmpi
    use m_io
    use m_parameters
    use m_fields
    use m_work
    use x_fftw

    implicit none

    integer :: i, j, k, n
    real*8 :: a, b, c, x, y, z

    ! defining very particular velocities so the RHS can be computed analytically

    if (task .eq. 'hydro') then

        write (out, *) 'inside.'
        call flush (out)

        a = 1.d0
        b = 1.d0
        c = 1.d0

        do k = 1, nz
            do j = 1, ny
                do i = 1, nx

                    x = dx*real(i - 1)
                    y = dx*real(j - 1)
                    z = dx*real(myid*nz + k - 1)

                    wrk(i, j, k, 1) = sin(a*x)
                    wrk(i, j, k, 2) = sin(b*y)
                    wrk(i, j, k, 3) = sin(c*z)
                end do
            end do
        end do

        write (out, *) 'did work'
        call flush (out)

        do n = 1, 3
            call xFFT3d(1, n)
            fields(:, :, :, n) = wrk(:, :, :, n)
        end do

        write (out, *) 'did FFTs'
        call flush (out)

        nu = 0.d0

        call rhs_velocity

        write (out, *) 'got rhs'
        call flush (out)

        do n = 1, 3
            call xFFT3d(-1, n)
        end do

        write (out, *) 'did FFTs'
        call flush (out)

        do k = 1, nz
            do j = 1, ny
                do i = 1, nx

                    x = dx*real(i - 1)
                    y = dx*real(j - 1)
                    z = dx*real(myid*nz + k - 1)

                    ! checking u
                    wrk(i, j, k, 4) = -sin(a*x)*(two*a*cos(a*x) + b*cos(b*y) + c*cos(c*z) + nu*a**2)

                    ! checking v
                    wrk(i, j, k, 5) = -sin(b*y)*(two*b*cos(b*y) + a*cos(a*x) + c*cos(c*z) + nu*b**2)

                    ! checking w
                    wrk(i, j, k, 6) = -sin(c*z)*(two*c*cos(c*z) + b*cos(b*y) + a*cos(a*x) + nu*c**2)

                end do
            end do
        end do

!!$    do k = 1,nz
!!$      write(out,"(3e15.6)") wrk(1,1,k,3),wrk(1,1,k,5),wrk(1,1,k,4)
!!$    end do

        wrk(:, :, :, 0) = &
            abs(wrk(:, :, :, 1) - wrk(:, :, :, 4)) + &
            abs(wrk(:, :, :, 2) - wrk(:, :, :, 5)) + &
            abs(wrk(:, :, :, 3) - wrk(:, :, :, 6))

        print *, 'Maximum error is ', maxval(wrk(1:nx, :, :, 0))

        tmp4(:, :, :) = wrk(1:nx, :, :, 1) - wrk(1:nx, :, :, 4)
        fname = 'e1.arr'
        call write_tmp4

        tmp4(:, :, :) = wrk(1:nx, :, :, 2) - wrk(1:nx, :, :, 5)
        fname = 'e2.arr'
        call write_tmp4

        tmp4(:, :, :) = wrk(1:nx, :, :, 3) - wrk(1:nx, :, :, 6)
        fname = 'e3.arr'
        call write_tmp4

    end if
    return
end subroutine test_rhs_velocity

! calculates vorticity using fields (in fourier space) and sets wrk(:,:,:1-3) = omg x omgy omgz
! WARNING: This function will overwrite any values currently in wrk
subroutine calculate_vorticity()
    use m_fields
    use m_work
    use x_fftw
    use m_parameters ! dx

    implicit none

    integer :: i,j,k

    ! velocities in Fourier space
    wrk(:, :, :, 1:3) = fields(:, :, :, 1:3)
    ! Taking derivatives
    ! derivative of (1) WRT (2) -> store it in (3)
    call x_derivative(3, 'y', 6) ! dw/dy -> 6
    call x_derivative(3, 'x', 5) ! dw/dx -> 5
    call x_derivative(2, 'z', 4) ! dv/dz -> 4
    call x_derivative(2, 'x', 3) ! dv/dx -> 3
    call x_derivative(1, 'z', 2) ! du/dz -> 2
    call x_derivative(1, 'y', 1) ! du/dy -> 1

    ! getting vorticity
    wrk(:, :, :, 3) = wrk(:, :, :, 3) - wrk(:, :, :, 1)  ! omega_3 = v_x - u_y
    wrk(:, :, :, 2) = wrk(:, :, :, 2) - wrk(:, :, :, 5)  ! omega_2 = u_z - w_x
    wrk(:, :, :, 1) = wrk(:, :, :, 6) - wrk(:, :, :, 4)  ! omega_1 = w_y - v_z

    ! wrk(:, :, :, 1:3) = fields(:, :, :, 1:3)
    !
    ! ! invert to physical space
    ! do i = 1,3
    !     call xFFT3d(-1, i)
    ! end do

    ! call gradient3D(nx, ny, nz, &
    !     wrk(1:nx, 1:ny, 1:nz, 3), & ! WRT W
    !     dx, dx, dx, &    ! step sizes in each direction
    !     wrk(1:nx, 1:ny, 1:nz, 5), & ! dW / dx -> 5
    !     wrk(1:nx, 1:ny, 1:nz, 6), & ! dW / dy -> 6
    !     wrk(1:nx, 1:ny, 1:nz, 0)  & !
    ! )

    ! call gradient3D(nx, ny, nz, &
    !     wrk(1:nx, 1:ny, 1:nz, 2), & ! WRT V
    !     dx, dx, dx, &    ! step sizes in each direction
    !     wrk(1:nx, 1:ny, 1:nz, 3), & ! dV / dx -> 3
    !     wrk(1:nx, 1:ny, 1:nz, 0), & !
    !     wrk(1:nx, 1:ny, 1:nz, 4)  & ! dV / dz -> 4
    ! )

    ! call gradient3D(nx, ny, nz, &
    !     wrk(1:nx, 1:ny, 1:nz, 1), & ! WRT U
    !     dx, dx, dx, &    ! step sizes in each direction
    !     tmp_wrk(1:nx, 1:ny, 1:nz, 8), & ! this is calculated before the next section so this is actually ok
    !     tmp_wrk(1:nx, 1:ny, 1:nz, 9), & ! dU / dy -> 1
    !     wrk(1:nx, 1:ny, 1:nz, 2)  & ! dU / dZ -> 2
    ! )

    ! write(*,*) "DU / DY"
    ! do k = 1,nz
    !     write(*,*) "k = ", k -1
    !     do j = 1,ny
    !         do i = 1,nx
    !             write(*,*) "i = ", i-1, "j = ", j-1, "velocity x: ", tmp_wrk(i,j,k,9)
    !         end do
    !     end do
    ! end do

    ! wrk(1:nx, 1:ny, 1:nz, 3) = wrk(1:nx, 1:ny, 1:nz, 3) - tmp_wrk(1:nx, 1:ny, 1:nz, 9)  ! omega_3 = v_x - u_y
    ! wrk(1:nx, 1:ny, 1:nz, 2) = wrk(1:nx, 1:ny, 1:nz, 2) - wrk(1:nx, 1:ny, 1:nz, 5)  ! omega_2 = u_z - w_x
    ! wrk(1:nx, 1:ny, 1:nz, 1) = wrk(1:nx, 1:ny, 1:nz, 6) - wrk(1:nx, 1:ny, 1:nz, 4)  ! omega_1 = w_y - v_z

    ! ! transform to fourier space since thats what the rest of the code expects
    ! do i = 1,3
    !     call xFFT3d(1, i)
    ! end do
end

SUBROUTINE gradient3D(m, n, o, f, hx, hy, hz, dfdx, dfdy, dfdz)
    INTEGER, INTENT(IN) :: m, n, o
    REAL*8, DIMENSION(m, n, o), INTENT(IN) :: f
    REAL*8, INTENT(IN) :: hx, hy, hz
    REAL*8, DIMENSION(m, n, o), INTENT(OUT) :: dfdx, dfdy, dfdz
    INTEGER :: i, j, k

    dfdx = 0.0; dfdy = 0.0; dfdz = 0.0;
    !dfdx
    do k = 1, o
        do j = 1, n
    !        forward difference at start point
            dfdx(1, j, k) = (f(2, j, k) - f(1, j, k))/hx;
    !        central difference at middle region
            dfdx(2:m - 1, j, k) = (f(3:m, j, k) - f(1:m - 2, j, k))/(2.*hx);
    !        backward difference at end point
            dfdx(m, j, k) = (f(m, j, k) - f(m - 1, j, k))/hx;
        end do
    end do

    !dfdy
    do k = 1, o
        do i = 1, m
    !        forward difference at start point
            dfdy(i, 1, k) = (f(i, 2, k) - f(i, 1, k))/hy;
    !        central difference at middle region
            dfdy(i, 2:n - 1, k) = (f(i, 3:n, k) - f(i, 1:n - 2, k))/(2.*hy);
    !        backward difference at end point
            dfdy(i, n, k) = (f(i, n, k) - f(i, n - 1, k))/hy;
        end do
    end do

    !dfdz
    do j = 1, n
        do i = 1, m
!        forward difference at start point
            dfdz(i, j, 1) = (f(i, j, 2) - f(i, j, 1))/hz;
!        central difference at middle region
            dfdz(i, j, 2:o - 1) = (f(i, j, 3:o) - f(i, j, 1:o - 2))/(2.*hz);
!        backward difference at end point
            dfdz(i, j, o) = (f(i, j, o) - f(i, j, o - 1))/hz;
        end do
    end do

END SUBROUTINE gradient3D

! the input to this subroutine has the wrk array in fourier space
! the output arrays should also be in fourier space
! wrk(:,:,:,1-3) contains the convective terms
! wrk(:,:,:,4) contains the diffusion term
subroutine update_forcing_viscous_compensation(epsilon_1, epsilon_2)
    use m_work
    use m_fields
    use x_fftw
    use forcing_vaues
    use m_parameters

    implicit none

    ! F_1, D_i and d/dt (Q_1)
    real*8 :: dQ_1, dQ_2
    real*8 :: diffusion_x, diffusion_y, diffusion_z
    real*8 :: epsilon_1, epsilon_2
    real*8 :: new_epsilon_1, new_epsilon_2
    real*8 :: f_left, f_right, f_total, diffusion
    integer :: i,j,k,n

    if (ITIME == 1) then
        write(out, *) "calculating forcing components with visc compensation routine"
    end if

    F_1 = 0.
    D_1 = 0.
    dQ_1 = 0.

    F_2 = 0.
    D_2 = 0.
    dQ_2 = 0.

    ! copy the first 3 terms (used for derivatives or something) to the temp array
    ! so that we can use the first three indicies to store the results of the diffusion term
    tmp_wrk(:,:,:,1:3) = wrk(:,:,:,1:3)

    ! copy the diffusion term to the 7th slot
    tmp_wrk(:,:,:,7) = wrk(:,:,:,4)

    ! FIRST:
    ! we evaluate the diffusion term from the RHS in fourier space
    !
    ! this code is replicated (read: Copied) from the for-loop that is directly AFTER this subroutine is called
    ! copy all the old varaibles to tmp_work
    if (viscous_compensation ==1) then
        call calculate_diffusion_term()
    end if

    ! after calling calculate_diffusion_term the diffusion term WRT u,v,w is now populated in wrk(:,:,:,1-3)
    ! so we can now use it to calculate whatever derivatives we may need. Swap the diffusion terms to X-space
    ! so we can use them when calculating the forcing

    do n=1,3
        call xFFT3d(-1, n)
    end do

    ! Lets copy these diffusion terms
    ! to tmp_wrk so that we can use all 6 values of wrk() when calculating vorticity

    tmp_wrk(:,:,:,4:6) = wrk(:,:,:,1:3)

    ! so now:
    ! tmp_wrk(:,:,:,1) -> some fourier derivative term we need to return to wrk() in this position
    ! tmp_wrk(:,:,:,2) -> some fourier derivative term we need to return to wrk() in this position
    ! tmp_wrk(:,:,:,3) -> some fourier derivative term we need to return to wrk() in this position
    ! tmp_wrk(:,:,:,4) -> diffusion WRT x (X-space)
    ! tmp_wrk(:,:,:,5) -> diffusion WRT y (X-space)
    ! tmp_wrk(:,:,:,6) -> diffusion WRT z (X-space)

    ! now - we dont have anything useful in the wrk array so we can use wrk to calculate the vorticity
    ! which requires 6 slots in wrk

    call calculate_vorticity()

    ! now:
    ! wrk(:,:,:,1) -> omg x (fourier-space)
    ! wrk(:,:,:,2) -> omg y (fourier-space)
    ! wrk(:,:,:,3) -> omg z (fourier-space)

    ! copy the velocities into wrk(:,:,:,4-6) so that we can shamelessly copy-paste more code from murali

    wrk(:, :, :, 4:6) = fields(:, :, :, 1:3)

    ! convert all wrk variables to X-space
    do n = 1, 6
        call xFFT3d(-1, n)
    end do

    ! now:
    ! wrk(:,:,:,1) -> omg x (x-space)
    ! wrk(:,:,:,2) -> omg y (x-space)
    ! wrk(:,:,:,3) -> omg z (x-space)
    ! wrk(:,:,:,4) -> u (x-space)
    ! wrk(:,:,:,5) -> v (x-space)
    ! wrk(:,:,:,6) -> w (x-space)

    ! ******2. Compute forcing terms & take FFTs
    ! u \cdot omg
    fcomp(:, :, :, 0) = wrk(:, :, :, 1)*wrk(:, :, :, 4) + wrk(:, :, :, 2)*wrk(:, :, :, 5) + wrk(:, :, :, 3)*wrk(:, :, :, 6)
    ! ||omg||^2
    fcomp(:, :, :, 1) = wrk(:, :, :, 1)**2 + wrk(:, :, :, 2)**2 + wrk(:, :, :, 3)**2
    ! ||u||^2
    fcomp(:, :, :, 2) = wrk(:, :, :, 4)**2 + wrk(:, :, :, 5)**2 + wrk(:, :, :, 6)**2

    do i=1,nx
        do j=1,ny
            do k=1,nz

                ! for each n direction
                do n = 1,3
                    !
                    ! Calculate the forcing components for 1 and 2
                    !

                    ! u \cdot omg * omg_i - |omg|^2 * u_i
                    f_left = fcomp(i,j,k,0) * wrk(i,j,k,n) - fcomp(i,j,k,1) * wrk(i,j,k,3+n)
                    ! u \cdot omg * u_i - |u|^2 * omg_i
                    f_right = fcomp(i,j,k,0) * wrk(i,j,k,3+n) - fcomp(i,j,k,2) * wrk(i,j,k,n)

                    ! store fcomp left and right so they can be exported to a VTK file later
                    ! in a 3D flowfield

                    fcomp_individual(i,j,k,n) = f_left
                    fcomp_individual(i,j,k,n+3) = f_right

                    ! the total forcing
                    f_total = (f_left * epsilon_1) + (f_right * epsilon_2)

                    ! integrate all the forcing components without the corresponding
                    ! epsilons into each F term
                    !
                    ! these values are outside the if statement because we want to know F_1 and F_2
                    ! for all the time steps so that we can (sometimes) output them into the energy.csv file
                    F_1 = F_1 + (f_left * wrk(i,j,k,3+n))
                    F_2 = F_2 + (f_right * wrk(i,j,k,n))

                    diffusion = tmp_wrk(i,j,k,n+3)

                    ! D_1 = u \cdot d_u
                    D_1 = D_1 + wrk(i,j,k,n+3) * diffusion

                    ! D_2 = \omega \cdot d_u
                    D_2 = D_2 + wrk(i,j,k,n) * diffusion

                    if (viscous_compensation == 1) then

                        ! d/dt Q_1 = u \cdot (d_u + f_u)
                        dQ_1 = dQ_1 +&
                            wrk(i,j,k,n+3) * (diffusion + f_total)

                        ! d/dt Q_2 = \omega \cdot (d_u + f_u)
                        dQ_2 = dQ_2 +&
                            wrk(i,j,k,n) * (diffusion + f_total)
                    else
                        ! forcing results if no viscous compensation go in 3:5
                        fcomp(i,j,k,n+2) = f_total
                    end if
                end do

                ! end
            end do
        end do
    end do

    ! now we have evaluated the integral so we can set the forcing components to their true value
    if (viscous_compensation == 1) then

        F_1 = -1 * F_1
        F_2 = -1 * F_2

        if (viscous_compensation_validation == 1) then
            dQ_1 = 0.0
            dQ_2 = 0.0
        end if

        if (skip_diffusion == 1) then
            D_1 = 0.0
            D_2 = 0.0
        end if

        new_epsilon_1 = epsilon_1 + ((D_1 - dQ_1)/F_1)
        new_epsilon_2 = epsilon_2 + ((D_2 - dQ_2)/F_2)

        ! recalculate the forcing components with the new values
        do n = 1, 3
            fcomp(:, :, :, 2 + n) = new_epsilon_1*(fcomp(:, :, :, 0)*wrk(:, :, :, n) - fcomp(:, :, :, 1)*wrk(:, :, :, 3 + n)) + &
                                    new_epsilon_2*(fcomp(:, :, :, 0)*wrk(:, :, :, 3 + n) - fcomp(:, :, :, 2)*wrk(:, :, :, n))
        end do
    end if

    ! copy either forcing result to the wrk array (we dont care about any of the data in it anymore)
    ! so that we can go back to fourier space

    wrk(:,:,:,1:3) = fcomp(:,:,:,3:5)

    ! Do the IFFT over the 3 forcing components
    do n = 1, 3
        call xFFT3d(1, n)
    end do

    fcomp(:,:,:, 1:3) = wrk(:,:,:,1:3)

    ! copy back whatever fourier derivatives used to be at this position
    wrk(:,:,:,1:3) = tmp_wrk(:,:,:,1:3)

    ! copy back the original diffusion term (fourier space) and just recalculate diffusion components as normal
    ! This is a missed optimization in the parent routine

    wrk(:,:,:,4) = tmp_wrk(:,:,:,7)
end

! calculate and store the result of the diffusion term
! in wrk(:,:,:,4)
! in order to call this subroutine you must have ensured that the laplace
! operation has been evaluated in fourier space in the main rhs_velocity
! initialization and that it is stored in wrk(:,:,:)
!
! This subroutine does _NOT_ take care of any forcing
subroutine calculate_diffusion_term()
    use m_work      !wrk
    use m_fields  ! fields
    use x_fftw ! ialias

    implicit none

    integer :: k, j, i, n

    do k = 1, nz
        do j = 1, ny
            do i = 1, nx + 1, 2

                if (ialias(i, j, k) .gt. 0) then ! run for dealias=0
                    ! setting the Fourier components to zero
                    wrk(i, j, k, 1:3) = zip
                    wrk(i + 1, j, k, 1:3) = zip

                else
                    ! RHS for u, v and w
                    do n = 1, 3
                        ! for the inviscid case
                        if (skip_diffusion == 1) then
                            wrk(i + 1, j, k, n) = 0.
                            wrk(i, j, k, n) =    0.
                        ! for the viscouse case - we evaluate the diffusion term
                        else
                            !write(*,*) "calculating diffusion terms"
                            ! the reason these indicies look like this is the convective term has
                            ! a negative in front of it in the impt in the main rhs_velocity subroutine
                            ! that we never fixed.
                            !wrk(i+1,j,k,n)  =   wrk(i  ,j,k,n) + wrk(i+1,j,k,4) * fields(i+1,j,k,n)
                            !wrk(i, j, k, n) =- wrk(i+1,j,k,n)  + wrk(i  ,j,k,4) * fields(i  ,j,k,n)
                            wrk(i+1,j,k,n)  = wrk(i+1,j,k,4) * fields(i+1,j,k,n)
                            wrk(i, j, k, n) = wrk(i  ,j,k,4) * fields(i  ,j,k,n)
                        end if
                    end do

                end if

            end do
        end do
    end do

end


!=====================================================|
subroutine artdif
!=====================================================|
  use geometry_def, only: qqm,qqp,nxg,dx,dt,mtype,cc
  use const_def, only: cv,cp,gm,rr,pii4
  implicit none

  integer :: i,m
  real(KIND(0.d0)), parameter :: fh = 0.d0, ep = 1.d0
  real(KIND(0.d0)), dimension(nxg,mtype) :: fx,dui
  real(KIND(0.d0)) :: dul,dur,duc,mup,mlo
  real(KIND(0.d0)) :: ul,ur,du,ra,pp
  real(KIND(0.d0)) :: dro,drx,dry,drz,dbx,dby,dbz,dse

  real(KIND(0.d0)), dimension(nxg) :: pr,te,hh ! pressure, temperature, enthalpy
!-----------------------------------------------------|

  do i = 1,nxg
     qqp(i,1) = qqm(i,1)
     qqp(i,2) = qqm(i,2)
     qqp(i,3) = qqm(i,3)
     qqp(i,4) = qqm(i,4)
     qqp(i,5) = qqm(i,5)
     qqp(i,6) = qqm(i,6)
     qqp(i,7) = qqm(i,7)
     qqp(i,8) = qqm(i,8)

     pr(i) = qqp(i,1)**gm*exp(qqp(i,8)/cv)
     te(i) = pr(i)/rr/qqp(i,1)
     hh(i) = cp*te(i) + 0.5d0*(qqp(i,2)**2 + qqp(i,3)**2 + qqp(i,4)**2)
  enddo

  ! calculating artificial diffusivity flux
  do i = 2,nxg-1
  do m = 1,mtype
     dul = qqp(i  ,m) - qqp(i-1,m)
     dur = qqp(i+1,m) - qqp(i  ,m)
     duc = 0.5d0*(qqp(i+1,m) - qqp(i-1,m))

     mup = max(ep*dul,ep*dur,duc)
     mlo = min(ep*dul,ep*dur,duc)

     dui(i,m) = min(0.d0,mup) + max(0.d0,mlo)
  enddo
  enddo

  do i = 2,nxg-2
  do m = 1,mtype
     ul = qqp(i  ,m) + 0.5d0*dui(i  ,m)
     ur = qqp(i+1,m) - 0.5d0*dui(i+1,m)
     du = qqp(i+1,m) - qqp(i,m)
     du = sign(1.d0,du)*max(abs(du),1.d-10)
     ra = min((ur-ul)/du,1.d0)
     pp = (0.5d0 + sign(0.5d0,ra))*max(0.d0,1.d0 + fh*(ra-1.d0))
     fx(i,m) = -0.5d0*(ur - ul)*pp*max(cc(i),cc(i+1))
     !fx(i,m) = -0.5d0*(ur - ul)*pp*0.5d0*(cc(i)+cc(i+1))
  enddo
  enddo

  do i = 2,nxg-2
     fx(i,8) = &
          & + 0.5d0*(qqp(i,1)*te(i) + qqp(i+1,1)*te(i+1))*fx(i,8) &
          & + 0.5d0*(qqp(i,1)*qqp(i,2) + qqp(i+1,1)*qqp(i+1,2))*fx(i,2) &
          & + 0.5d0*(qqp(i,1)*qqp(i,3) + qqp(i+1,1)*qqp(i+1,3))*fx(i,3) &
          & + 0.5d0*(qqp(i,1)*qqp(i,4) + qqp(i+1,1)*qqp(i+1,4))*fx(i,4) &
          & + 0.5d0*(qqp(i,5) + qqp(i+1,5))*fx(i,5)*pii4 &
          & + 0.5d0*(qqp(i,6) + qqp(i+1,6))*fx(i,6)*pii4 &
          & + 0.5d0*(qqp(i,7) + qqp(i+1,7))*fx(i,7)*pii4

     fx(i,2) = + 0.5d0*(qqp(i,1) + qqp(i+1,1))*fx(i,2) &
          &    + 0.5d0*(qqp(i,2) + qqp(i+1,2))*fx(i,1)

     fx(i,3) = + 0.5d0*(qqp(i,1) + qqp(i+1,1))*fx(i,3) &
          &    + 0.5d0*(qqp(i,3) + qqp(i+1,3))*fx(i,1)

     fx(i,4) = + 0.5d0*(qqp(i,1) + qqp(i+1,1))*fx(i,4) &
          &    + 0.5d0*(qqp(i,4) + qqp(i+1,4))*fx(i,1)
  enddo

  do i = 3,nxg-2
     dro = -(fx(i,1) - fx(i-1,1))/dx*dt
     drx = -(fx(i,2) - fx(i-1,2))/dx*dt
     dry = -(fx(i,3) - fx(i-1,3))/dx*dt
     drz = -(fx(i,4) - fx(i-1,4))/dx*dt
     dbx = -(fx(i,5) - fx(i-1,5))/dx*dt
     dby = -(fx(i,6) - fx(i-1,6))/dx*dt
     dbz = -(fx(i,7) - fx(i-1,7))/dx*dt
     dse =  ( &
          & - (fx(i,8) - fx(i-1,8))/dx*dt &
          & - (hh(i) - qqp(i,2)**2 - qqp(i,3)**2 - qqp(i,4)**2)*dro &
          & - (qqp(i,2)*drx + qqp(i,3)*dry + qqp(i,4)*drz) &
          & - (qqp(i,5)*dbx + qqp(i,6)*dby + qqp(i,7)*dbz)*pii4 &
          & )
     qqm(i,1) = qqm(i,1) + dro
     qqm(i,2) = (qqp(i,1)*qqm(i,2) + drx)/qqm(i,1)
     qqm(i,3) = (qqp(i,1)*qqm(i,3) + dry)/qqm(i,1)
     qqm(i,4) = (qqp(i,1)*qqm(i,4) + drz)/qqm(i,1)
     qqm(i,5) = qqm(i,5) + dbx
     qqm(i,6) = qqm(i,6) + dby
     qqm(i,7) = qqm(i,7) + dbz
     qqm(i,8) = qqm(i,8) + dse/qqp(i,1)/te(i)
  enddo

  call bc(qqm)
  return
end subroutine artdif


text: Describe Users/Olivia Glass here.

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
!
!     CalculiX - A 3-dimensional finite element program
!              Copyright (C) 1998-2015 Guido Dhondt
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation(version 2);
!
!
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!     GNU General Public License for more details.
!
!     You should have received a copy of the GNU General Public License
!     along with this program; if not, write to the Free Software
!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
!
!     completing hel:
!
!     at the start of the subroutine: rhs of conservation of momentum
!                                     without pressure contribution
!     at the end of the subroutine: neighboring velocity terms subtracted
!
      subroutine complete_hel(neq,bv,hel,adv,auv,jq,irow,nzs)
!
      implicit none
!
      integer irow(*),neq,nzs,j,l,jdof1,jq(*)
      real*8 hel(3,*),bv(neq,3),auv(*),adv(*)
c      do j=1,neq
c      write(*,*) 'complete_hel_1 ',j,hel(1,j),hel(2,j),hel(3,j)
c      enddo
!
!     off-diagonal terms
!
      do j=1,neq
         do l=jq(j),jq(j+1)-1
            jdof1=irow(l)
!
!           subdiagonal terms
!
            hel(1,jdof1)=hel(1,jdof1)-auv(l)*bv(j,1)
            hel(2,jdof1)=hel(2,jdof1)-auv(l)*bv(j,2)
            hel(3,jdof1)=hel(3,jdof1)-auv(l)*bv(j,3)
!
!           superdiagonal terms
!
            hel(1,j)=hel(1,j)-auv(l+nzs)*bv(jdof1,1)
            hel(2,j)=hel(2,j)-auv(l+nzs)*bv(jdof1,2)
            hel(3,j)=hel(3,j)-auv(l+nzs)*bv(jdof1,3)
         enddo
      enddo
c      do j=1,neq
c         write(*,*) 'complete_hel ',j,bv(j,1),bv(j,2),bv(j,3)
c      enddo
!
      return
      end

max_stars_repo_name: alleindrach/calculix-desktop

max_issues_repo_name: alleindrach/calculix-desktop

max_forks_repo_name: alleindrach/calculix-desktop


text:
C    Copyright(C) 1999-2020 National Technology & Engineering Solutions
C    of Sandia, LLC (NTESS).  Under the terms of Contract DE-NA0003525 with
C    NTESS, the U.S. Government retains certain rights in this software.
C
C    See packages/seacas/LICENSE for details

      SUBROUTINE GETSBC (MXND, MXNPER, NPER, NL, ML, MAXSBC, MAXPRM,
     &   NPRM, NID, LISTL, XN, YN, NUID, LXK, KXL, NXL, LSTSBC, NPERIM,
     &   KSBC, LCON, ISBOUN, LINKL, NSPF, IFSB, LISTSB, LINKSB, LLL,
     &   BAR, ERR)
C***********************************************************************

C  SUBROUTINE GETSBC = GETS THE SIDE BOUNDARY LIST

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

      DIMENSION NID (MXNPER, MAXPRM), NPERIM (MAXPRM)
      DIMENSION LISTL (NL), XN (MXND), YN (MXND), NUID (MXND)
      DIMENSION LXK (4, MXND), KXL (2, MXND*3), NXL (2, MXND*3)
      DIMENSION LCON (3, ML), ISBOUN (ML), LINKL (2, ML)
      DIMENSION NSPF (ML), IFSB (ML), LISTSB (2, ML), LINKSB (2, ML)
      DIMENSION LSTSBC (MAXSBC), NODES (4)

      LOGICAL EXTER, ERR, CCW, BAR, ADDLNK

      ERR = .TRUE.
      CCW = .TRUE.
      ADDLNK = .FALSE.
      NPERIM (1) = NPER

      DO 110 II = 1, NPRM
         DO 100 I = 1, NPERIM(II)
            IF (BAR) THEN
               NID (I, II) = IABS (NUID (I))
            ELSE
               IF (NID (I, II) .LT. 0) NID (I, II) = - NID (I, II)
            ENDIF
  100    CONTINUE
  110 CONTINUE

C  SORT THROUGH AND PICK OFF ELEMENTS WITH SIDE BOUNDARY CONDITIONS

      DO 240 I = 1, LLL
         IF (BAR) THEN
            I1 = LXK (1, I)
            I2 = LXK (2, I)
         ELSE
            I1 = NXL (1, I)
            I2 = NXL (2, I)
         ENDIF

C  SEE IF THE LINE IS CLEARLY INTERIOR

         IF (I1 .GT. 0 .AND. I2 .GT. 0) THEN
           if ((NUID (I1) .NE. 0) .AND. (NUID (I2) .NE. 0)) THEN
            LTEST = 0
            EXTER = .FALSE.

C  CHECK AGAINST THE PERIMETER LIST TO SEE IF IT IS TRULY EXTERIOR

            DO 130 JJ  =  1, NPRM
               DO 120 J = 1, NPERIM (JJ)
                  IF (ABS (NUID (I1)) .EQ. NID (J, JJ)) THEN
                     IF (J .EQ. 1) THEN
                        J1 = J + 1
                        J2 = NPERIM(JJ)
                     ELSEIF (J .EQ. NPERIM(JJ)) THEN
                        J1 = J - 1
                        J2 = 1
                     ELSE
                        J1 = J - 1
                        J2 = J + 1
                     ENDIF
                     IF ( (ABS (NUID (I2)) .EQ. NID (J1, JJ)) .OR.
     &                  (ABS (NUID (I2)) .EQ. NID (J2, JJ)) )
     &                  EXTER = .TRUE.
                     GOTO 140
                  ENDIF
  120          CONTINUE
  130       CONTINUE
  140       CONTINUE
            IF (EXTER) THEN

C  FIND THE LINE NUMBER IT BELONGS TO

               IF (ABS (NUID (I1)) .GT. 1000000000) THEN
                  LTEST =  (ABS (NUID (I1)) - 1000000000) / 100000
               ELSEIF (ABS (NUID (I2)) .GT. 1000000000) THEN
                  LTEST =  (ABS (NUID (I2)) - 1000000000) / 100000
               ELSE
                  NSUM = ABS (NUID (I1)) + ABS (NUID (I2))
                  DO 150 J = 1, NL
                     CALL LTSORT (ML, LINKL, LISTL (J), K, ADDLNK)
                     IF ((LCON (1, K) + LCON (2, K)) .EQ. NSUM) THEN
                        IF (( (LCON (1, K) .EQ. ABS (NUID (I1))) .AND.
     +                     (LCON (2, K) .EQ. ABS (NUID (I2)))) .OR.
     +                     ((LCON (1, K) .EQ. ABS (NUID (I2))) .AND.
     +                     (LCON (2, K) .EQ. ABS (NUID (I1))))) THEN
                           LTEST = LISTL (J)
                           GOTO 160
                        ENDIF
                     ENDIF
  150             CONTINUE
  160             CONTINUE
               ENDIF

C  FIND THE ELEMENT BOUNDARY FLAG IF THERE IS ONE

               IF (LTEST.LE.0) THEN
                  CALL MESSAGE(' ERROR IN SEARCHING NXL FOR '//
     &               'ELEMENT BCC')
                  RETURN
               ELSE
                  CALL LTSORT (ML, LINKL, LTEST, J, ADDLNK)
                  IF (ISBOUN (J) .GT. 0) THEN
                     IFLAG = ISBOUN (J)

C  CHECK TO MAKE SURE LINE IS LINKED TO FLAG
C  AND GET THE NEXT LINK  (NFLAG)

                     CALL LTSORT (ML, LINKSB, IFLAG, L, ADDLNK)
                     DO 170 JJ = IFSB (L), IFSB (L) + NSPF (L) - 1
                        IF (LISTSB (1, JJ) .LT. 0) THEN
                           CALL MESSAGE('PROBLEMS WITH SIDES IN '//
     &                        'FLAG LIST IN GETSBC')
                        ELSE
                           IF (LISTSB (1, JJ) .EQ. LTEST) THEN
                              NFLAG = LISTSB (2, JJ)
                              GOTO 180
                           ENDIF
                        ENDIF
  170                CONTINUE
                     WRITE (*, 10000)IFLAG
                     RETURN
  180                CONTINUE
                     IF (BAR) THEN
                        NELEM = I
                     ELSE
                        NELEM = KXL (1, I)
                        IF (NELEM .EQ. 0)NELEM = KXL (2, I)
                     ENDIF
                     KSBC = KSBC + 1
                     LSTSBC (KSBC) = - IFLAG
                     KSBC = KSBC + 1
                     if (ksbc .gt. maxsbc) stop 'maxsbc error'
                     LSTSBC (KSBC) = NELEM

C  GET THE CORRECT ELEMENT SIDE

                     IF (BAR) THEN
                        JSIDE = 1
                     ELSE
                        CALL GNXKA (MXND, XN, YN, NELEM, NODES, AREA,
     &                     LXK, NXL, CCW)
                        DO 190 J = 1, 4
                           IF (I1 .EQ. NODES (J)) THEN
                              JP1 = J + 1
                              JM1 = J - 1
                              IF (JP1 .EQ. 5)JP1 = 1
                              IF (JM1 .EQ. 0)JM1 = 4
                              IF (I2 .EQ. NODES (JP1)) THEN
                                 JSIDE = J
                                 GOTO 200
                              ELSEIF (I2 .EQ. NODES (JM1)) THEN
                                 JSIDE = JM1
                                 GOTO 200
                              ENDIF
                           ENDIF
  190                   CONTINUE
                        WRITE (*, 10010)NELEM
                        RETURN
  200                   CONTINUE
                     ENDIF
                     KSBC = KSBC + 1
                     LSTSBC (KSBC) = JSIDE

C  SEE IF ANY MORE FLAGS ARE ATTACHED TO THIS SIDE

  210                CONTINUE
                     IF (NFLAG .GT. 0) THEN

C  CHECK TO MAKE SURE LINE IS LINKED TO FLAG
C  AND GET THE NEXT LINK  (NFLAG)

                        IFLAG = NFLAG
                        CALL LTSORT (ML, LINKSB, IFLAG, L, ADDLNK)
                        DO 220 JJ = IFSB (L), IFSB (L) + NSPF (L)
                           IF (LISTSB (1, JJ) .EQ. LTEST) THEN
                              NFLAG = LISTSB (2, JJ)
                              GOTO 230
                           ENDIF
  220                   CONTINUE
                        WRITE (*, 10000)IFLAG
                        RETURN
  230                   CONTINUE
                        KSBC = KSBC + 1
                        LSTSBC (KSBC) = - IFLAG
                        KSBC = KSBC + 1
                        LSTSBC (KSBC) = NELEM
                        KSBC = KSBC + 1
                        LSTSBC (KSBC) = JSIDE
                        GOTO 210
                     ENDIF
                  ENDIF
               ENDIF
            ENDIF
         ENDIF
       END IF
  240 CONTINUE

      ERR = .FALSE.
      RETURN

10000 FORMAT (' SIDE BOUNDARY FLAG', I5, ' IS NOT PROPERLY LINKED')
10010 FORMAT (' ERROR FINDING CORRECT BOUNDARY SIDE ON ELEMENT', I5)
      END

max_stars_repo_name: jschueller/seacas

max_issues_repo_name: jschueller/seacas

max_forks_repo_name: jschueller/seacas


text:
! { dg-do compile }
! { dg-options "-std=f95" }
abstract interface ! { dg-error "Fortran 2003: ABSTRACT INTERFACE" }
  subroutine two()
  end subroutine two
end interface ! { dg-error "Expecting END PROGRAM statement" }
end

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
      SUBROUTINE CALAVB2 (ISTL_)
C
C **  SUBROUTINE CALAV CALCULATES VERTICAL VISCOSITY AND DIFFUSIVITY
C **  USING GLAPERIN ET AL'S MODIFICATION OF THE MELLOR-YAMADA MODEL
C **  (NOTE AV, AB, AND AQ ARE ACTUALLY DIVIDED BY H)
C **  IF ISGA=1 VALUES ARE GEOMETRIC AVERAGES WITH THE PREVIOUS VALUES
C CHANGE RECORD
C
      USE GLOBAL

	IMPLICIT NONE
C
      !REAL*8 HPD(LCM),DMLD(LCM,KCM),BD(LCM,KCM),DZIGD(KCM)
      !REAL*8 CTURBB1D(LCM,KCM),CTURBB2D(LCM,KCM),QQD(LCM,KCM)
      !REAL*8 CTURBD,CTURB2BD,GPD
      !REAL*8 RIQMIN,RAVBTMP,DELBTMP,RITMP,BFUN,RIQ
      !REAL*8 TMPVAL,TMPVAL1,SBTOP,SVTOP,TMPVAL2,SVTOP2,SVBOT
      !REAL*8 SFAV,SFAB,SBBOT,TMPVAL3
      !REAL*8 ATURB1,ATURB2,TURBC1
      !REAL*8 AVTMP,ABTMP

      REAL*8 RIQMIN,DELBTMP,RITMP,BFUN,RIQ
      REAL*8 TMPVAL,TMPVAL1,SBTOP,SVTOP,TMPVAL2,SVTOP2,SVBOT
      REAL*8 SFAV,SFAB,SBBOT,TMPVAL3
      REAL*8 ATURB1,ATURB2,TURBC1
      REAL AVTMP,ABTMP,AQTMP
C
      INTEGER K,L,LS,ISTL_
C
      DATA ATURB1,ATURB2,TURBC1/0.92,0.74,0.08/
C
C   SHTOP    =      0.4939
C   SHBOT    =     34.6764
C   SMTOP1   =      0.3933
C   SMTOP2   =      7.8464
C   SMBOT1   =     34.6764
C   SMBOT2   =      6.1272
C   RLIMIT   =      0.0233
C   SHMIN    =      0.0934
C   SMMIN    =      0.1099
C   SHMAX    =      5.2073
C   SMMAX    =      4.9639
C
      RIQMIN=-0.023
C
      DO K=1,KS
        DO L=2,LA
          CTURBB1(L,K)=CTURB
          CTURBB2(L,K)=CTURB2B   ! PMC ONLY NEEDS TO BE ASSIGNED ONCE SINCE NEVER CHANGED
        ENDDO
      ENDDO
C
      DO K=1,KS    ! PMC SKIP IF BSC.LE.1.E-6
        DO L=2,LA
          DELBTMP=(B(L,K+1)-B(L,K))*DZIG(K)
          RITMP=-GP*HP(L)*DELBTMP/QQ(L,K)
          IF(RITMP.GT.0.)THEN
            RIQ=DML(L,K)*DML(L,K)*RITMP
            BFUN=EXP(-3.11*RIQ)
            CTURBB1(L,K)=CTURB/(BFUN+1.E-12)

            ! *** Original Code
            IF(BBT(L,K).GT.0.)THEN   ! *** PMC BBT is never set, so this is never used
              TMPVAL=DELBTMP*DELBTMP/(RITMP*BBT(L,K))
              CTURBB2(L,K)=CTURB2B/(1.+0.61*(1.-BFUN)*TMPVAL)
            ENDIF
          ENDIF
        ENDDO
      ENDDO

      DO K=1,KS
        DO L=2,LA
          RIQ=-GP*HP(L)*DML(L,K)*DML(L,K)*DZIG(K)
     &            *(B(L,K+1)-B(L,K))/QQ(L,K)
          RIQ=MAX(RIQ,RIQMIN)
          TMPVAL1=1.-( 6.*ATURB1/CTURBB1(L,K) )
          SBTOP=ATURB2*TMPVAL1
          SBBOT=3.*ATURB2*( 6.*ATURB1+CTURBB2(L,K) )
          SVTOP=ATURB1*( TMPVAL1-3.*TURBC1)
          TMPVAL2=TMPVAL1*( CTURBB2(L,K)-3.*ATURB2 )
          TMPVAL3=-3.*TURBC1*( CTURBB2(L,K)+6.*ATURB1)
          SVTOP2=3.*ATURB2*(TMPVAL2+TMPVAL3)/SVTOP
          SVBOT=9.*ATURB1*ATURB2
          SFAV=SVTOP*(1.+SVTOP2*RIQ)/((1.+SVBOT*RIQ)*(1.+SBBOT*RIQ))
          SFAB=SBTOP/(1.+SBBOT*RIQ)
          ABTMP=AVCON*SFAB*DML(L,K)*HP(L)*QQSQR(L,K)+ABO
          AVTMP=AVCON*SFAV*DML(L,K)*HP(L)*QQSQR(L,K)+AVO
C
          IF(ISFAVB.EQ.0)THEN
            AV(L,K)=AVTMP*HPI(L)
            AB(L,K)=SCB(L)*ABTMP*HPI(L)
          ELSEIF(ISFAVB.EQ.1)THEN
            AV(L,K)=0.5*(AV(L,K)+AVTMP*HPI(L))
            AB(L,K)=SCB(L)*0.5*(AB(L,K)+ABTMP*HPI(L))
          ELSEIF(ISFAVB.EQ.2)THEN
            AV(L,K)=SQRT(AV(L,K)*AVTMP*HPI(L))
            AB(L,K)=SCB(L)*SQRT(AB(L,K)*ABTMP*HPI(L))
          ENDIF
        ENDDO
      ENDDO

      IF(ISLOG.GE.1)THEN
        AVMAX=AVO
        ABMAX=ABO
        AVMIN=10.
        ABMIN=10.
        DO K=1,KS
          DO L=2,LA
            AVMAX=MAX(AVMAX,AV(L,K))
            ABMAX=MAX(ABMAX,AB(L,K))
            AVMIN=MIN(AVMIN,AV(L,K))
            ABMIN=MIN(ABMIN,AB(L,K))
          ENDDO
        ENDDO
      ENDIF

      ! *** NOW APPLY MAXIMUM, IF REQURIED
      IF(ISAVBMX.GE.1)THEN
        DO K=1,KS
          DO L=2,LA
            AVTMP=AVMX*HPI(L)
            ABTMP=ABMX*HPI(L)
            AV(L,K)=MIN(AV(L,K),AVTMP)
            AB(L,K)=MIN(AB(L,K),ABTMP)
          ENDDO
        ENDDO
      ENDIF

      DO K=1,KS
        DO L=2,LA
          LS=LSC(L)
          AVUI(L,K)=2./(AV(L,K)+AV(LWEST(L),K))
          AVVI(L,K)=2./(AV(L,K)+AV(LS,K))
        ENDDO
      ENDDO

      IF(ISTL_.EQ.3)THEN
        DO K=2,KS
          DO L=2,LA
            AQ(L,K)=0.205*(AV(L,K-1)+AV(L,K))
          ENDDO
        ENDDO
        DO L=2,LA
          AQ(L,1)=0.205*AV(L,1)
          AQ(L,KC)=0.205*AV(L,KS)
        ENDDO
      ELSE
        DO K=2,KS
          DO L=2,LA
            AQTMP=0.205*(AV(L,K-1)+AV(L,K))
            AQ(L,K)=AQTMP
          ENDDO
        ENDDO
        DO L=2,LA
          AQTMP=0.205*AV(L,1)
          AQ(L,1)=AQTMP
          AQTMP=0.205*AV(L,KS)
          AQ(L,KC)=AQTMP
        ENDDO
      ENDIF

      RETURN
      END


!!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!! 2018.03.07: fortran lib for python
!!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!!==============================================================================
!!==============================================================================
    SUBROUTINE coupling_prediction(N,L,S,St1,dS,M,C_inv,W2,H02,cost2)
    IMPLICIT NONE
    !! W here = W in python code (the same order: from j--> i: W(i,j))
    INTEGER,PARAMETER ::  nloop=1000
    INTEGER :: N,L,i,iloop,stop_iteration,n_iter,i0

    REAL (KIND=8) :: H_av
    REAL (KIND=8) :: S(L,N),St1(L,N),dS(L,N),W(N),H0(nloop),cost(0:nloop),H(L)
    REAL (KIND=8) :: HS_av(N),C_inv(N,N),M(N),W1(nloop,N),W2(N,N),H02(N),cost2(N)

   !f2py intent(in) :: S,St1,dS,M,C_inv
   !f2py intent(hide), depend(S) :: L=shape(S,0),N=shape(S,1)
   !f2py intent(out) :: W2,H02,cost2

    DO i0=1,N
    H(1:L)=St1(1:L,i0) ; iloop=1 ; cost(0)=100. ; stop_iteration=0
    !!--------------------------------------------------------------
    DO WHILE ((iloop<=nloop).and.(stop_iteration==0))
        H_av=sum(H(1:L))/L

        DO i=1,N
        HS_av(i)=sum((H(1:L)-H_av)*dS(1:L,i))/L
        END DO

        W(:)=matmul(HS_av(:),C_inv(:,:))
        H0(iloop)=H_av - sum(W(:)*M(:))

        H(1:L)=matmul(S(1:L,:),W(:))+H0(iloop)

        cost(iloop)=sum((St1(1:L,i0)-tanh(H(1:L)))**2)/L
        !WRITE(*,*)iloop,cost(iloop)

        IF (cost(iloop)>=cost(iloop-1)) THEN
         stop_iteration=1
        END IF

        H(:)=St1(:,i0)*H(:)/tanh(H(:))
        W1(iloop,:)=W(:)

        iloop=iloop+1
    END DO
    !!--------------------------------------------------------------
    n_iter=iloop-2

    W2(i0,:)=W1(n_iter,:)
    cost2(i0)=cost(n_iter)
    H02(i0)=H0(n_iter)

    !WRITE(*,*)i0,cost2(i0)
    END DO

    END SUBROUTINE
!!==============================================================================
!!==============================================================================
    SUBROUTINE update_hidden_ens(N1,N2,L0,W,H0,S,S_out,nens)
    IMPLICIT NONE
    !! W here = W.T in python code
    INTEGER :: N1,N2,t,i,L0,L01,nens
    REAL (KIND=8) :: W(N2,N2),H0(N2),S(L0,N2),S_out(L0,N1:N2)
    REAL (KIND=8) :: H11(N1:N2),P11(N1:N2),H12(N2),P1,P2,rdn_S(L0,N1:N2)

    !f2py intent(in) :: N1,W,H0,S,nens
    !f2py intent(hide), depend(S) :: L0=shape(S,0),N2=shape(W,0)
    !f2py intent(out) :: S_out

    L01=L0/nens
    CALL random_number(rdn_S)

    DO t=1,L0
        IF (mod(t,L01) > 1) THEN  !! t = 2, .., L01-1
            H11(N1:N2)=H0(N1:N2)+matmul(S(t-1,:),W(:,N1:N2))
            P11(N1:N2)=1/(1+exp(-2*H11(N1:N2)))

            DO i=N1,N2
                S(t,i)=1.0
                H12(:)=H0(:)+matmul(S(t,:),W(:,:))
                P1=P11(i)/product(1+exp(-2*S(t+1,:)*H12(:)))
                P2=(1-P11(i))/product(1+exp(-2*S(t+1,:)*(H12(:)-2.*W(i,:))))
                S(t,i) = sign(1.d+00,(P1/(P1+P2)-rdn_S(t,i)))
            END DO
        END IF
    END DO

    !DO t=1,L0
    !    IF (mod(t,L01) > 1) THEN  !! t = 2, .., L01-1
    !        H11(N1:N2)=H0(N1:N2)+matmul(S(t-1,:),W(:,N1:N2))
    !        P11(N1:N2)=1/(1+exp(-2*H11(N1:N2)))

     !       DO i=N1,N2
     !           S(t,i)=1.0
     !           H12(:)=H0(:)+matmul(S(t,:),W(:,:))
     !           P1=P11(i)/product(1+exp(-2*S(t+1,:)*H12(:)))
     !           P2=(1-P11(i))/product(1+exp(-2*S(t+1,:)*(H12(:)-2.*W(i,:))))
     !           S(t,i) = sign(1.d+00,(P1/(P1+P2)-rdn_S(t,i)))
      !      END DO
      !  ELSE

        !IF (mod(t,L01) == 0) THEN  !! t = L01
        !    H11(N1:N2)=H0(N1:N2)+matmul(S(t-1,:),W(:,N1:N2))
        !    P11(N1:N2)=1/(1+exp(-2*H11(N1:N2)))
        !    S(t,N1:N2) = sign(1.d+00,(P11(N1:N2)-rdn_S(t,N1:N2)))
        !ELSE  !! t=1,L01+1,...
        !    DO i=N1,N2
        !        S(t,i)=1.0
        !        H12(:)=H0(:)+matmul(S(t,:),W(:,:))
        !        P1=1/product(1+exp(-2*S(t+1,:)*H12(:)))
        !        P2=1/product(1+exp(-2*S(t+1,:)*(H12(:)-2.*W(i,:))))
        !        S(t,i) = sign(1.d+00,(P1/(P1+P2)-rdn_S(t,i)))
        !    END DO
        !END IF
        !END IF

    !END DO

    S_out(:,N1:N2)=S(:,N1:N2)

    END SUBROUTINE
!!==============================================================================
   SUBROUTINE update_hidden(N1,N2,L,W,H0,S,S_out)
   IMPLICIT NONE
   !! W here = W.T in python code

   INTEGER :: N1,N2,t,L,i
   REAL (KIND=8) :: W(N2,N2),H0(N2),S(L,N2),S_out(L,N1:N2)
   REAL (KIND=8) :: H11(N1:N2),P11(N1:N2),H12(N2),P1,P2,rdn_S(L,N1:N2)

   !f2py intent(in) :: N1,W,H0,S
   !f2py intent(hide), depend(S) :: L=shape(S,0),N2=shape(W,0)
   !f2py intent(out) :: S_out

   CALL random_number(rdn_S)

   DO t=2,L-1

      H11(N1:N2)=H0(N1:N2)+matmul(S(t-1,:),W(:,N1:N2))
      P11(N1:N2)=1/(1+exp(-2*H11(N1:N2)))

      DO i=N1,N2
         S(t,i)=1.0
         H12(:)=H0(:)+matmul(S(t,:),W(:,:))
         P1=P11(i)/product(1+exp(-2*S(t+1,:)*H12(:)))
         P2=(1-P11(i))/product(1+exp(-2*S(t+1,:)*(H12(:)-2.*W(i,:))))

         S(t,i) = sign(1.d+00,(P1/(P1+P2)-rdn_S(t,i)))

      END DO

   END DO

   S_out(:,N1:N2)=S(:,N1:N2)

   END SUBROUTINE


text:
!----------------------------------------------------------------------
! Module: module_apply_com.f90
!> Purpose: communication routines
!>

!> Author: Benjamin Menetrier
!>

!> Licensing: this code is distributed under the CeCILL-C license
!>

!> Copyright © 2017 METEO-FRANCE
!----------------------------------------------------------------------
module module_apply_com

use tools_display, only: msgerror
use type_fields, only: fldtype,alphatype,buftype
use type_sdata, only: sdatatype,sdatampitype
implicit none

private
public :: fld_com_gl,fld_com_lg
public :: alpha_com_AB,alpha_com_BA,alpha_com_CA,alpha_copy_AB,alpha_copy_AC,alpha_copy_BC

contains

!----------------------------------------------------------------------
! Subroutine: fld_com_gl
!> Purpose: communicate full field from global to local distribution
!----------------------------------------------------------------------
subroutine fld_com_gl(sdata,fld_global,fld_local)

implicit none

! Passed variables
type(sdatatype),intent(in) :: sdata                   !< Sampling data
type(fldtype),intent(inout) :: fld_global             !< Global field
type(fldtype),intent(inout) :: fld_local(sdata%nproc) !< Local field

! Local variables
integer :: ic0,iproc,ic0a

! Allocation
do iproc=1,sdata%nproc
   if (.not.allocated(fld_local(iproc)%vala)) allocate(fld_local(iproc)%vala(sdata%mpi(iproc)%nc0a,sdata%mpi(iproc)%nl0))
end do

! Loop over cells
do ic0=1,sdata%nc0
   iproc = sdata%ic0_to_iproc(ic0)
   ic0a = sdata%ic0_to_ic0a(ic0)
   fld_local(iproc)%vala(ic0a,:) = fld_global%val(ic0,:)
end do

end subroutine fld_com_gl

!----------------------------------------------------------------------
! Subroutine: fld_com_lg
!> Purpose: communicate full field from local to global distribution
!----------------------------------------------------------------------
subroutine fld_com_lg(sdata,fld_local,fld_global)

implicit none

! Passed variables
type(sdatatype),intent(in) :: sdata                   !< Sampling data
type(fldtype),intent(inout) :: fld_local(sdata%nproc) !< Local field
type(fldtype),intent(inout) :: fld_global             !< Global field

! Local variables
integer :: ic0,iproc,ic0a

! Allocation
if (.not.allocated(fld_global%val)) allocate(fld_global%val(sdata%nc0,sdata%nl0))

! Loop over cells
do ic0=1,sdata%nc0
   iproc = sdata%ic0_to_iproc(ic0)
   ic0a = sdata%ic0_to_ic0a(ic0)
   fld_global%val(ic0,:) = fld_local(iproc)%vala(ic0a,:)
end do

end subroutine fld_com_lg

!----------------------------------------------------------------------
! Subroutine: alpha_com_AB
!> Purpose: communicate reduced field from zone A to zone B
!----------------------------------------------------------------------
subroutine alpha_com_AB(sdata,alpha)

implicit none

! Passed variables
type(sdatatype),intent(in) :: sdata                 !< Sampling data
type(alphatype),intent(inout) :: alpha(sdata%nproc) !< Subgrid variable

! Local variables
integer :: iproc,jproc,ir_s,ir_e,is_s,is_e
type(buftype) :: bufs(sdata%nproc),bufr(sdata%nproc)

! Allocation
do iproc=1,sdata%nproc
   allocate(alpha(iproc)%valb(sdata%mpi(iproc)%nsb))
end do

if (sdata%nproc==1) then
   ! Copy data
   alpha(1)%valb = alpha(1)%vala
elseif (sdata%nproc>1) then
   ! Allocation
   do iproc=1,sdata%nproc
      allocate(bufs(iproc)%val(sdata%mpi(iproc)%AB%nexcl))
      allocate(bufr(iproc)%val(sdata%mpi(iproc)%AB%nhalo))
   end do

   ! Prepare buffers to send
   do iproc=1,sdata%nproc
      bufs(iproc)%val = alpha(iproc)%vala(sdata%mpi(iproc)%AB%excl)
   end do

   ! Emulate mpi_alltoallv
   do iproc=1,sdata%nproc
      do jproc=1,sdata%nproc
         ! Communication iproc->jproc
         is_s = sdata%mpi(iproc)%AB%jexcldispl(jproc)+1
         is_e = sdata%mpi(iproc)%AB%jexcldispl(jproc)+sdata%mpi(iproc)%AB%jexclcounts(jproc)
         ir_s = sdata%mpi(jproc)%AB%jhalodispl(iproc)+1
         ir_e = sdata%mpi(jproc)%AB%jhalodispl(iproc)+sdata%mpi(jproc)%AB%jhalocounts(iproc)
         if (is_e-is_s/=ir_e-ir_s) call msgerror('wrong mpi_alltoallv indices')
         if (ir_e>=ir_s) bufr(jproc)%val(ir_s:ir_e) = bufs(iproc)%val(is_s:is_e)
      end do
   end do

   do iproc=1,sdata%nproc
      ! Copy zone A into zone B
      alpha(iproc)%valb(sdata%mpi(iproc)%isa_to_isb) = alpha(iproc)%vala

      ! Copy halo into zone B
      alpha(iproc)%valb(sdata%mpi(iproc)%AB%halo) = bufr(iproc)%val
   end do
end if

! Release memory
do iproc=1,sdata%nproc
   deallocate(alpha(iproc)%vala)
end do

end subroutine alpha_com_AB

!----------------------------------------------------------------------
! Subroutine: alpha_com_BA
!> Purpose: communicate reduced field from zone B to zone A
!----------------------------------------------------------------------
subroutine alpha_com_BA(sdata,alpha)

implicit none

! Passed variables
type(sdatatype),intent(in) :: sdata                 !< Sampling data
type(alphatype),intent(inout) :: alpha(sdata%nproc) !< Subgrid variable

! Local variables
integer :: iproc,jproc,ir_s,ir_e,is_s,is_e,i
type(buftype) :: bufs(sdata%nproc),bufr(sdata%nproc)

! Allocation
do iproc=1,sdata%nproc
   allocate(alpha(iproc)%vala(sdata%mpi(iproc)%nsa))
end do

if (sdata%nproc==1) then
   ! Copy data
   alpha(1)%vala = alpha(1)%valb
elseif (sdata%nproc>1) then
   ! Allocation
   do iproc=1,sdata%nproc
      allocate(bufs(iproc)%val(sdata%mpi(iproc)%AB%nhalo))
      allocate(bufr(iproc)%val(sdata%mpi(iproc)%AB%nexcl))
   end do

   ! Prepare buffers to send
   do iproc=1,sdata%nproc
      bufs(iproc)%val = alpha(iproc)%valb(sdata%mpi(iproc)%AB%halo)
   end do

   ! Emulate mpi_alltoallv
   do iproc=1,sdata%nproc
      do jproc=1,sdata%nproc
         ! Communication iproc->jproc
         is_s = sdata%mpi(iproc)%AB%jhalodispl(jproc)+1
         is_e = sdata%mpi(iproc)%AB%jhalodispl(jproc)+sdata%mpi(iproc)%AB%jhalocounts(jproc)
         ir_s = sdata%mpi(jproc)%AB%jexcldispl(iproc)+1
         ir_e = sdata%mpi(jproc)%AB%jexcldispl(iproc)+sdata%mpi(jproc)%AB%jexclcounts(iproc)
         if (is_e-is_s/=ir_e-ir_s) call msgerror('wrong mpi_alltoallv indices')
         if (ir_e>=ir_s) bufr(jproc)%val(ir_s:ir_e) = bufs(iproc)%val(is_s:is_e)
      end do
   end do

   do iproc=1,sdata%nproc
      ! Copy zone B into zone A
      alpha(iproc)%vala = alpha(iproc)%valb(sdata%mpi(iproc)%isa_to_isb)

      ! Copy halo into zone B
      do i=1,sdata%mpi(iproc)%AB%nexcl
         alpha(iproc)%vala(sdata%mpi(iproc)%AB%excl(i)) = alpha(iproc)%vala(sdata%mpi(iproc)%AB%excl(i))+bufr(iproc)%val(i)
      end do
   end do
end if

! Release memory
do iproc=1,sdata%nproc
   deallocate(alpha(iproc)%valb)
end do

end subroutine alpha_com_BA

!----------------------------------------------------------------------
! Subroutine: alpha_com_CA
!> Purpose: communicate reduced field from zone C to zone A
!----------------------------------------------------------------------
subroutine alpha_com_CA(sdata,alpha)

implicit none

! Passed variables
type(sdatatype),intent(in) :: sdata                 !< Sampling data
type(alphatype),intent(inout) :: alpha(sdata%nproc) !< Subgrid variable

! Local variables
integer :: iproc,jproc,ir_s,ir_e,is_s,is_e,i
type(buftype) :: bufs(sdata%nproc),bufr(sdata%nproc)

! Allocation
do iproc=1,sdata%nproc
   allocate(alpha(iproc)%vala(sdata%mpi(iproc)%nsa))
end do

if (sdata%nproc==1) then
   ! Copy data
   alpha(1)%vala = alpha(1)%valc
elseif (sdata%nproc>1) then
   ! Allocation
   do iproc=1,sdata%nproc
      allocate(bufs(iproc)%val(sdata%mpi(iproc)%AC%nhalo))
      allocate(bufr(iproc)%val(sdata%mpi(iproc)%AC%nexcl))
   end do

   ! Prepare buffers to send
   do iproc=1,sdata%nproc
      bufs(iproc)%val = alpha(iproc)%valc(sdata%mpi(iproc)%AC%halo)
   end do

   ! Emulate mpi_alltoallv
   do iproc=1,sdata%nproc
      do jproc=1,sdata%nproc
         ! Communication iproc->jproc
         is_s = sdata%mpi(iproc)%AC%jhalodispl(jproc)+1
         is_e = sdata%mpi(iproc)%AC%jhalodispl(jproc)+sdata%mpi(iproc)%AC%jhalocounts(jproc)
         ir_s = sdata%mpi(jproc)%AC%jexcldispl(iproc)+1
         ir_e = sdata%mpi(jproc)%AC%jexcldispl(iproc)+sdata%mpi(jproc)%AC%jexclcounts(iproc)
         if (is_e-is_s/=ir_e-ir_s) call msgerror('wrong mpi_alltoallv indices')
         if (ir_e>=ir_s) bufr(jproc)%val(ir_s:ir_e) = bufs(iproc)%val(is_s:is_e)
      end do
   end do

   do iproc=1,sdata%nproc
      ! Copy zone C into zone A
      alpha(iproc)%vala = alpha(iproc)%valc(sdata%mpi(iproc)%isa_to_isc)

      ! Copy halo into zone C
      do i=1,sdata%mpi(iproc)%AC%nexcl
         alpha(iproc)%vala(sdata%mpi(iproc)%AC%excl(i)) = alpha(iproc)%vala(sdata%mpi(iproc)%AC%excl(i))+bufr(iproc)%val(i)
      end do
   end do
end if

! Release memory
do iproc=1,sdata%nproc
   deallocate(alpha(iproc)%valc)
end do

end subroutine alpha_com_CA

!----------------------------------------------------------------------
! Subroutine: alpha_copy_AB
!> Purpose: copy reduced field from zone A to zone B
!----------------------------------------------------------------------
subroutine alpha_copy_AB(sdatampi,alpha)

implicit none

! Passed variables
type(sdatampitype),intent(in) :: sdatampi !< Sampling data
type(alphatype),intent(inout) :: alpha    !< Subgrid variable

! Allocation
allocate(alpha%valb(sdatampi%nsb))

! Initialize
alpha%valb = 0.0

! Copy zone A into zone B
alpha%valb(sdatampi%isa_to_isb) = alpha%vala

! Release memory
deallocate(alpha%vala)

end subroutine alpha_copy_AB

!----------------------------------------------------------------------
! Subroutine: alpha_copy_AC
!> Purpose: copy reduced field from zone A to zone C
!----------------------------------------------------------------------
subroutine alpha_copy_AC(sdatampi,alpha)

implicit none

! Passed variables
type(sdatampitype),intent(in) :: sdatampi !< Sampling data
type(alphatype),intent(inout) :: alpha    !< Subgrid variable

! Allocation
allocate(alpha%valc(sdatampi%nsc))

! Initialize
alpha%valc = 0.0

! Copy zone A into zone C
alpha%valc(sdatampi%isa_to_isc) = alpha%vala

! Release memory
deallocate(alpha%vala)

end subroutine alpha_copy_AC

!----------------------------------------------------------------------
! Subroutine: alpha_copy_BC
!> Purpose: copy reduced field from zone B to zone C
!----------------------------------------------------------------------
subroutine alpha_copy_BC(sdatampi,alpha)

implicit none

! Passed variables
type(sdatampitype),intent(in) :: sdatampi !< Sampling data
type(alphatype),intent(inout) :: alpha    !< Subgrid variable

! Allocation
allocate(alpha%valc(sdatampi%nsc))

! Initialization
alpha%valc = 0.0

! Copy zone B into zone C
alpha%valc(sdatampi%isb_to_isc) = alpha%valb

! Release memory
deallocate(alpha%valb)

end subroutine alpha_copy_BC

end module module_apply_com

max_stars_repo_name: thomasauligne/NICAS

max_issues_repo_name: thomasauligne/NICAS

max_forks_repo_name: thomasauligne/NICAS


text:
C******************************************************************************
C FILE: mpi_bug6.f
C DESCRIPTION:
C   This program attempts to compare isend/irecv with send/irecv, but there is
C   a bug. Under AIX it will cause program failure. Intel Linux ignores it.
C   4 tasks are required.
C SOURCE: Blaise Barney
C LAST REVISED: 04/02/05
C*****************************************************************************/

      program mpi_bug6
      include 'mpif.h'

      integer REPS, DISP, numtasks, rank, buf, tag, i, tmp,
     +        ierr, dest, src, offset, nreqs
      parameter (DISP=100)
      parameter (REPS=1000)
      integer reqs(REPS*2), stats(MPI_STATUS_SIZE,REPS*2)
      real*8 T1, T2

      tag = 1

      call MPI_INIT(ierr)
      call MPI_COMM_RANK(MPI_COMM_WORLD, rank, ierr)
      call MPI_COMM_SIZE(MPI_COMM_WORLD, numtasks, ierr)

C     Require 4 tasks
      if (rank .EQ. 0) then
        if (numtasks .NE. 4) then
          print *, 'ERROR: Number of tasks must be 4. Quitting.'
          call MPI_ABORT(MPI_COMM_WORLD, ierr)
        endif
        print *, 'Starting isend/irecv send/irecv test...'
      endif

C     Use barriers for clean output
      call MPI_BARRIER(MPI_COMM_WORLD, ierr)
      print *, 'Task',rank,'starting...'
      call MPI_BARRIER(MPI_COMM_WORLD, ierr)

C     Start the clock
      T1 = MPI_WTIME()

C     Tasks 0 and 1 do the isend/irecv test.
C     Determine who to send/receive with. nreqs specifies how many non-blocking
C     operation request handles to capture. offset is where the task should
C     store each request as it is captured in the reqs() array.
      if (rank .LT. 2) then
        nreqs = REPS*2
        if (rank .EQ. 0) then
          src = 1
          offset = 1
        else
          src = 0
          offset = REPS
        endif
        dest = src

C       Do the non-blocking send and receive operations
        do i=1,REPS
          call MPI_ISEND(rank, 1, MPI_INTEGER, dest, tag,
     +                   MPI_COMM_WORLD, reqs(offset), ierr)
          call MPI_IRECV(buf, 1, MPI_INTEGER, src, tag,
     +                   MPI_COMM_WORLD, reqs(offset+1), ierr)
          offset = offset + 2
          tmp = MOD(i, DISP)
          if (tmp .EQ. 0) then
            print *, 'Task',rank,'has done',i,'isends/irecvs'
          endif
        enddo
      endif

C     Tasks 2 and 3 do the send/irecv test.
C     Determine who to send/receive with. nreqs specifies how many non-blocking
C     operation request handles to capture. offset is where the task should
C     store each request as it is captured in the reqs() array.
      if (rank .GT. 1) then
        nreqs = REPS

C       Task 2 does the blocking send operation
        if (rank .EQ. 2) then
          dest = 3
          do i=1,REPS
            call MPI_SEND(rank, 1, MPI_INTEGER, dest, tag,
     +                    MPI_COMM_WORLD, ierr)
            tmp = MOD(i, DISP)
            if (tmp .EQ. 0) then
              print *, 'Task',rank,'has done',i,'sends'
            endif
          enddo
        endif

C       Task 3 does the non-blocking receive operation
        if (rank .EQ. 3) then
          src = 2
          offset = 1
          do i=1,REPS
            call MPI_IRECV(buf, 1, MPI_INTEGER, src, tag,
     +                   MPI_COMM_WORLD, reqs(offset), ierr)
            offset = offset + 1
            tmp = MOD(i, DISP)
            if (tmp .EQ. 0) then
              print *, 'Task',rank,'has done',i,'irecvs'
            endif
          enddo
        endif

      endif

C     Wait for all non-blocking operations to complete and record time
      call MPI_WAITALL(nreqs, reqs, stats, ierr)
      T2 = MPI_WTIME()
      call MPI_BARRIER(MPI_COMM_WORLD, ierr)

      print *, 'Task',rank,'time(wall)=',T2-T1,' sec'

      call MPI_FINALIZE(ierr)

      end

max_stars_repo_name: hixio-mh/HPC-Tutorials

max_issues_repo_name: hixio-mh/HPC-Tutorials

max_forks_repo_name: hixio-mh/HPC-Tutorials


text:
a = "a";
b = " ";

while (length(a) < 10)
    a = a + "a";

while (length(b) < 20) {
    b = b + ".";
    print(b);
};

print(a);

max_stars_repo_name: git-bauerseb/M-Lang

max_issues_repo_name: git-bauerseb/M-Lang

max_forks_repo_name: git-bauerseb/M-Lang


text:
! computes trapezoidal method
program CompTrap
implicit none
double precision x0,xn,h,sum,f,exact,diff,xj
integer n,j
write(*,*)'enter n'
read(*,*) n
x0=0.0d0
xn=2.d0
h=(xn-x0)/dfloat(n)
sum=f(x0)+f(xn)
do j=1,n-1
  xj=x0+dfloat(j)*h
  sum=sum+2.d0*f(xj)
enddo
sum=h/2.d0*sum
exact=dlog(1.5d0)
write(*,*) 'exact=',exact
diff=dabs(sum-exact)
write(*,*) 'Trapezoidal approx=',sum
write(*,*) 'diff=',diff
stop
end program CompTrap

function f(x)
double precision x,f
f=1/(x+4)
return
end function f

max_stars_repo_name: whisker-rebellion/Fortran

max_issues_repo_name: whisker-rebellion/Fortran

max_forks_repo_name: whisker-rebellion/Fortran


text:
active component C {

  param set port paramSetOut

}

ext: fpp

max_stars_repo_name: kevin-f-ortega/fpp

max_issues_repo_name: kevin-f-ortega/fpp

max_forks_repo_name: kevin-f-ortega/fpp


text:
	SUBROUTINE IN_SKIP  ( skip, iskpcn, iskplt, iret )
C************************************************************************
C* IN_SKIP								*
C*									*
C* This subroutine decodes the skip factor which is in the form:	*
C*		skip_contour / skip_plot_x; skip_plot_y			*
C*									*
C* IN_SKIP  ( SKIP, ISKPCN, ISKPLT, IRET )				*
C*									*
C* Input parameters:							*
C*	SKIP		CHAR*		Skip factor input		*
C*									*
C* Output parameters:							*
C*	ISKPCN		INTEGER		Skip for contours		*
C*	ISKPLT (2)	INTEGER		Skip for plotting at x, y locs	*
C*	IRET		INTEGER		Return code			*
C*					  0 = normal return		*
C**									*
C* Log:									*
C* M. desJardins/NMC	12/91						*
C* L. Sager/NMC		 6/93		Stagger wind plots		*
C* S. Jacobs/NMC	 8/94		Cleaned up			*
C* S. Jacobs/NMC	 8/94		Force the user to type a / to	*
C*					   get the plot skip values	*
C* S. Jacobs/NMC	 9/94		Force neg values to be odd	*
C* S. Jacobs/NMC	 9/94		Default the Y-skip to 1/2 the	*
C*					   X-skip for staggering,	*
C*					   otherwise use the requested	*
C*					   value			*
C* D.W.Plummer/NCEP	10/96		Allow iskplt(2) to be negative.	*
C* T. Piper/GSC		11/98		Updated prolog			*
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
C*
	CHARACTER*(*)	skip
	INTEGER		iskplt (2)
C*
	CHARACTER	carr (2)*72
	INTEGER		iarr (2)
C------------------------------------------------------------------------
	iret = 0
C
C*	Break input into data and vector elements.
C
	CALL ST_CLST  ( skip, '/', ' ', 2, carr, num, ier )
C
C*	Get the contour skip value from the first string.
C
	CALL ST_ILST ( carr (1), ';', 0, 1, iskpcn, num, ier )
	IF  ( iskpcn .lt. 0 )  iskpcn = 0
C
C*	Parse the second string. Get the values for the x and y point
C*	increments.
C
	CALL ST_ILST ( carr (2), ';', IMISSD, 2, iarr, num, ier )
C
C*	Check for a missing x value.
C*	Also check for negative values which signal staggering.
C
	IF  ( iarr (1) .eq. IMISSD )  THEN
	    iskplt (1) = 0
	ELSE IF  ( iarr (1) .lt. 0 )  THEN
	    ival = ABS ( iarr (1) )
	    IF  ( MOD ( ival, 2 ) .eq. 0 )  THEN
		iskplt (1) = iarr (1) + 1
	    ELSE
		iskplt (1) = iarr (1)
	    END IF
	ELSE
	    iskplt (1) = iarr (1)
	END IF
C
C*	Check for a missing y value.
C
	IF  ( iarr (2) .eq. IMISSD )  THEN
	    IF  ( iskplt (1) .lt. 0 )  THEN
		iskplt (2) = ABS ( iskplt (1) ) / 2
	    ELSE
		iskplt (2) = iskplt (1)
	    END IF
	ELSE
	    iskplt (2) = iarr (2)
	END IF
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
SUBROUTINE TIMOSHENKO
! . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
! .                                                                   .
! .   To set up storage and call the BEAMELE element subroutine
! .   CALCULATED THE BEAM ELEMENT STIFFNESS MATRIX
!.    POSTPROCESSING FOR THE BEAM ELEMENT
!
!     LIUCHANGWU                     .
! . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  USE GLOBALS
  USE MEMALLOCATE

  IMPLICIT NONE
  INTEGER :: NUME, NUMMAT, MM, N(14)

  NUME = NPAR(2)
  NUMMAT = NPAR(3)
  NPAR(5) = 2

! Calculate the pointer to the arrays in the element group data
! N101: E(NUMMAT)
! N102: G(NUMMAT)
! N103: AREA(NUMMAT)
! N104: I_X(NUMMAT)
! N105: I_Y(NUMMAT)
! N106: I_Z(NUMMAT)
! N107: J_X(NUMMAT)
! N108: J_Y(NUMMAT)
! N109: J_Z(NUMMAT)
! N110: LM(12,NUME)
! N111: XYZ(6,NUME)
! N112: MTAP(NUME)

  N(1)=0
  N(2)=N(1)+NUMMAT*ITWO
  N(3)=N(2)+NUMMAT*ITWO
  N(4)=N(3)+NUMMAT*ITWO
  N(5)=N(4)+NUMMAT*ITWO
  N(6)=N(5)+NUMMAT*ITWO
  N(7)=N(6)+NUMMAT*ITWO
  N(8)=N(7)+NUMMAT*ITWO
  N(9)=N(8)+NUMMAT*ITWO
  N(10)=N(9)+NUMMAT*ITWO
  N(11)=N(10)+12*NUME
  N(12)=N(11)+6*NUME*ITWO
  N(13)=N(12)+NUME
  N(14)=N(13)+NPAR(5)*NPAR(2)

  MIDEST=N(14)
  if (IND .EQ. 1) then
        ! Allocate storage for element group data
        call MemAlloc(11,"ELEGP",MIDEST,1)
  end if
  NFIRST = NP(11)   ! Pointer to the first entry in the element group data array in the unit of single precision (corresponding to A)
  N(:) = N(:) + NFIRST
  NLAST=N(14)

  CALL TIMOELE (IA(NP(1)),DA(NP(2)),DA(NP(3)),DA(NP(4)),DA(NP(4)),IA(NP(5)),   &
       A(N(1)),A(N(2)),A(N(3)),A(N(4)),A(N(5)),A(N(6)),A(N(7)),A(N(8)),A(N(9)),A(N(10)),A(N(11)),A(N(12)),A(N(13)))

  RETURN

END SUBROUTINE TIMOSHENKO

SUBROUTINE TIMOELE (ID,X,Y,Z,U,MHT,E,G,AREA,I_X,I_Y,I_Z,J_X,J_Y,J_Z,LM,XYZ,MATP,Node)
! . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
! .                                                                   .
! .   TIMOSHENKO ELEMENT subroutine                                        .
! .
! .   OUTPUT: S(12,12) THE ELEMENT STIFFNESS MATRIX
! . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  USE GLOBALS
  USE MEMALLOCATE

  IMPLICIT NONE
  INTEGER :: ID(6,NUMNP),LM(12,NPAR(2)),MATP(NPAR(2)),MHT(NEQ)
  REAL(8) :: X(NUMNP),Y(NUMNP),Z(NUMNP),E(NPAR(3)),G(NPAR(3)),AREA(NPAR(3)),  &
             I_X(NPAR(3)),I_Y(NPAR(3)),I_Z(NPAR(3)),J_X(NPAR(3)),J_Y(NPAR(3)),J_Z(NPAR(3)),XYZ(6,NPAR(2)),U(NEQ)
  REAL(8) :: S(12,12),D(3),UELE(12),FORCE(12),SIGMA, M(12,12), Rho, Density(NPAR(3))   !UELE: THE DISPLACEMENT OF ELEMENT N; FORCE: THE INTERNAL FORCE AND MOMENT FOR ELEMENT N

  INTEGER :: NPAR1, NUME, NUMMAT, ND, I, J, L, N, Node(NPAR(2),NPAR(5))
  INTEGER :: MTYPE, IPRINT
  REAL(8) :: XL2, XL ,CXX,CXY,CXZ,CYX,CYY,CYZ,CZX,CZY,CZZ,T(12,12),KFACTOR,POISSON!TRANSLATION MATRIX
  REAL(8) :: StressCollection(6,NPAR(2)), GaussianCollection(3,NPAR(2))

  NPAR1  = NPAR(1)
  NUME   = NPAR(2)
  NUMMAT = NPAR(3)

  ND=12

! Read and generate element information
  IF (IND .EQ. 1) THEN

     WRITE (IOUT,"(' E L E M E N T   D E F I N I T I O N',//,  &
                   ' ELEMENT TYPE ',13(' .'),'( NPAR(1) ) . . =',I10,/,   &
                   '     EQ.1, TRUSS ELEMENTS',/,      &
                   '     EQ.2, ELEMENTS CURRENTLY',/,  &
                   '     EQ.5, 3D BEAM',//,      &
                   ' NUMBER OF ELEMENTS.',10(' .'),'( NPAR(2) ) . . =',I10,/)") NPAR1,NUME

     IF (NUMMAT.EQ.0) NUMMAT=1

     WRITE (IOUT,"(' M A T E R I A L   D E F I N I T I O N',//,  &
                   ' NUMBER OF DIFFERENT SETS OF MATERIAL',/,  &
                   ' AND CROSS-SECTIONAL  CONSTANTS ',         &
                   4 (' .'),'( NPAR(3) ) . . =',I10,/)") NUMMAT

     WRITE (IOUT,"('  SET       YOUNG''S     SHEARING     AREA     I_X     I_Y     I_Z     J_X     J_Y     J_Z',/,  &
                   ' NUMBER     MODULUS',10X,'MODULUS',/,  &
                   15 X,'E',14X,'A')")

     DO I=1,NUMMAT
        READ (IIN,'(I10,9F10.0)') N,E(N),G(N),AREA(N),I_X(N),I_Y(N),I_Z(N),J_X(N),J_Y(N),J_Z(N)  ! Read material information
        WRITE (IOUT,"(I10,1X,E12.5,1X,E12.5,1X,E12.5,1X,E12.5,1X,E12.5,1X, E12.5,1X,E12.5, 1X,E12.5,1X,E12.5)") &
        N,E(N),G(N),AREA(N),I_X(N),I_Y(N),I_Z(N),J_X(N),J_Y(N),J_Z(N)
     END DO

     WRITE (IOUT,"(//,' E L E M E N T   I N F O R M A T I O N',//,  &
                      ' ELEMENT     NODE     NODE       MATERIAL',/,   &
                      ' NUMBER-N      I        J       SET NUMBER')")

     N=0
     DO WHILE (N .NE. NUME)
        READ (IIN,'(5I10)') N,Node(N,1:NPAR(5)),MTYPE  ! Read in element information

!       Save element information
        XYZ(1:NPAR(5)*3-1:3,N)=X(Node(N,:))  ! Coordinates of the element's nodes
        XYZ(2:NPAR(5)*3  :3,N)=Y(Node(N,:))
        XYZ(3:NPAR(5)*3+1:3,N)=Z(Node(N,:))

        MATP(N)=MTYPE  ! Material type

        DO L=1,12
           LM(L,N)=0
        END DO

        DO L=1,6
           LM(L,N)=ID(L,Node(N,1))     ! Connectivity matrix
           LM(L+6,N)=ID(L,Node(N,2))
        END DO

!       Update column heights and bandwidth
        if (.NOT. PARDISODOOR) CALL COLHT (MHT,ND,LM(1,N))

        WRITE (IOUT,"(I10,6X,I10,4X,I10,7X,I10)") N,Node(N,1:NPAR(5)),MTYPE

        write (VTKNodeTmp) NPAR(5), Node(N,:)-1

     END DO

     RETURN

! Assemble stucture stiffness matrix
  ELSE IF (IND .EQ. 2) THEN

     DO N=1,NUME
        MTYPE=MATP(N)

        !CALCULATE THE LENGTH OF THE ELEMENT

        XL2=0.
        DO L=1,3
           D(L)=XYZ(L,N) - XYZ(L+3,N)
           XL2=XL2 + D(L)*D(L)
        END DO
        XL=SQRT(XL2)

        ! COORDINATE TRANSLATION MATRIX T

        DO I=1,12
            DO J=1,12
                T(I,J)=0
            ENDDO
        ENDDO

        !ASSUME THE DIRECTION FOR PRINCIPLE AXIS

        CXX=(XYZ(4,N)-XYZ(1,N))/XL
        CXY=(XYZ(5,N)-XYZ(2,N))/XL
        CXZ=(XYZ(6,N)-XYZ(3,N))/XL
        IF ((ABS(CXX-0.0)<10E-5).AND.(ABS(CXY-0.0)<10E-5))THEN
          CYX=1.0
          CYY=0.0
          CYZ=0.0
          CZX=0.0
          CZY=1.0
          CZZ=0.0

        ELSE
          CYX=-CXX*CXZ/SQRT(CXX*CXX+CXY*CXY)
          CYY=-CXY*CXZ/SQRT(CXX*CXX+CXY*CXY)
          CYZ=SQRT(CXX*CXX+CXY*CXY)
          CZX=CXY/SQRT(CXX*CXX+CXY*CXY)
          CZY=-CXX/SQRT(CXX*CXX+CXY*CXY)
          CZZ=0.0
        ENDIF

        T(1,1)=CXX
        T(1,2)=CXY
        T(1,3)=CXZ
        T(2,1)=CYX
        T(2,2)=CYY
        T(2,3)=CYZ
        T(3,1)=CZX
        T(3,2)=CZY
        T(3,3)=CZZ

        DO I=1,3
           DO J=1,3
              T(I+3,J+3)=T(I,J)
              T(I+6,J+6)=T(I,J)
              T(I+9,J+9)=T(I,J)
           ENDDO
        ENDDO

        !CALCULATE ELEMENT STIFFNESS MATRIX

        DO I=1,12
           DO J=1,12
              S(I,J)=0
           ENDDO
        ENDDO

        !MODIFICATION FACTOR FOR THE SHEARING STRESS
        POISSON=E(MTYPE)/(2.0*G(MTYPE))-1.0
        KFACTOR=10.0*(1.0+POISSON)/(12.0+11.0*POISSON)

        S(1,1)=E(MTYPE)*AREA(MTYPE)/XL
        S(2,2)=G(MTYPE)*AREA(MTYPE)*KFACTOR/XL
        S(3,3)=G(MTYPE)*AREA(MTYPE)*KFACTOR/XL
        S(4,4)=G(MTYPE)*J_Z(MTYPE)/XL
        S(5,5)=E(MTYPE)*I_Y(MTYPE)/XL+G(MTYPE)*AREA(MTYPE)*KFACTOR*XL/4.0
        S(6,6)=E(MTYPE)*I_Z(MTYPE)/XL+G(MTYPE)*AREA(MTYPE)*KFACTOR*XL/4.0
        S(7,7)=E(MTYPE)*AREA(MTYPE)/XL
        S(8,8)=G(MTYPE)*AREA(MTYPE)*KFACTOR/XL
        S(9,9)=G(MTYPE)*AREA(MTYPE)*KFACTOR/XL
        S(10,10)=G(MTYPE)*J_Z(MTYPE)/XL
        S(11,11)=E(MTYPE)*I_Y(MTYPE)/XL+G(MTYPE)*AREA(MTYPE)*KFACTOR*XL/4.0
        S(12,12)=E(MTYPE)*I_Z(MTYPE)/XL+G(MTYPE)*AREA(MTYPE)*KFACTOR*XL/4.0
        S(5,3)=-G(MTYPE)*AREA(MTYPE)*KFACTOR/2.0
        S(6,2)=G(MTYPE)*AREA(MTYPE)*KFACTOR/2.0
        S(7,1)=-E(MTYPE)*AREA(MTYPE)/XL
        S(8,2)=-G(MTYPE)*AREA(MTYPE)*KFACTOR/XL
        S(8,6)=-G(MTYPE)*AREA(MTYPE)*KFACTOR/2.0
        S(9,3)=-G(MTYPE)*AREA(MTYPE)*KFACTOR/XL
        S(9,5)=G(MTYPE)*AREA(MTYPE)*KFACTOR/2.0
        S(10,4)=-G(MTYPE)*J_X(MTYPE)/XL
        S(11,3)=-G(MTYPE)*AREA(MTYPE)*KFACTOR/2.0
        S(11,5)=-E(MTYPE)*I_Y(MTYPE)/XL+G(MTYPE)*AREA(MTYPE)*KFACTOR*XL/4.0
        S(11,9)=G(MTYPE)*AREA(MTYPE)*KFACTOR/2.0
        S(12,2)=G(MTYPE)*AREA(MTYPE)*KFACTOR/2.0
        S(12,6)=-E(MTYPE)*I_Z(MTYPE)/XL+G(MTYPE)*AREA(MTYPE)*KFACTOR*XL/4.0
        S(12,8)=-G(MTYPE)*AREA(MTYPE)*KFACTOR/2.0

        DO I=1,11
           DO J=I+1,12
               S(I,J)=S(J,I)
           ENDDO
        ENDDO

        S=MATMUL(TRANSPOSE(T),MATMUL(S,T))

        if(pardisodoor) then
            if(huge) then
                call pardiso_addban(stff,IA(NP(2)),columns,S,LM(1,N),ND)
            else
                call pardiso_addban(DA(NP(3)),IA(NP(2)),IA(NP(5)),S,LM(1,N),ND)
            end if
        else
            CALL ADDBAN (DA(NP(3)),IA(NP(2)),S,LM(1,N),ND)
        end if

     END DO

     RETURN

! Stress calculations
  ELSE IF (IND .EQ. 3) THEN

     IPRINT=0
     DO N=1,NUME
        IPRINT=IPRINT + 1
        IF (IPRINT.GT.50) IPRINT=1
        IF (IPRINT.EQ.1) WRITE (IOUT,"(//,' S T R E S S  C A L C U L A T I O N S  F O R  ',  &
              'E L E M E N T  G R O U P',I4,//,   &
              '  ELEMENT',2X,'FORCE_X1',3X,'FORCE_Y1',3X,'FORCE_Z1',3X,'MOMENT_X1',3X,'MOMENT_Y1',3X, 'MOMENT_Z1',3X, &
              'FORCE_X2',3X,'FORCE_Y2',3X,'FORCE_Z2',3X,'MOMENT_X2',3X,'MOMENT_Y2',3X,'MOMENT_Z2',/,'  NUMBER')") NG
        MTYPE=MATP(N)

        !CALCULATE THE LENGTH OF THE ELEMENT

        XL2=0.
        DO L=1,3
           D(L) = XYZ(L,N) - XYZ(L+3,N)
           XL2=XL2 + D(L)*D(L)
        END DO

         XL=SQRT(XL2)

        ! COORDINATE TRANSLATION MATRIX T

        DO I=1,12
            DO J=1,12
                T(I,J)=0
            ENDDO
        ENDDO

        !ASSUME THE DIRECTION FOR PRINCIPLE AXIS

        CXX=(XYZ(4,N)-XYZ(1,N))/XL
        CXY=(XYZ(5,N)-XYZ(2,N))/XL
        CXZ=(XYZ(6,N)-XYZ(3,N))/XL

        IF ((ABS(CXX-0.0)<10E-5).AND.(ABS(CXY-0.0)<10E-5))THEN
          CYX=1.0
          CYY=0.0
          CYZ=0.0
          CZX=0.0
          CZY=1.0
          CZZ=0.0

        ELSE
          CYX=-CXX*CXZ/SQRT(CXX*CXX+CXY*CXY)
          CYY=-CXY*CXZ/SQRT(CXX*CXX+CXY*CXY)
          CYZ=SQRT(CXX*CXX+CXY*CXY)
          CZX=CXY/SQRT(CXX*CXX+CXY*CXY)
          CZY=-CXX/SQRT(CXX*CXX+CXY*CXY)
          CZZ=0.0
        ENDIF

        T(1,1)=CXX
        T(1,2)=CXY
        T(1,3)=CXZ
        T(2,1)=CYX
        T(2,2)=CYY
        T(2,3)=CYZ
        T(3,1)=CZX
        T(3,2)=CZY
        T(3,3)=CZZ

        DO I=1,3
           DO J=1,3
              T(I+3,J+3)=T(I,J)
              T(I+6,J+6)=T(I,J)
              T(I+9,J+9)=T(I,J)
           ENDDO
        ENDDO

        !CALCULATE ELEMENT STIFFNESS MATRIX

        DO I=1,12
           DO J=1,12
              S(I,J)=0
           ENDDO
        ENDDO

        !MODIFIED FACTOR FOR SHEARING STRESS
        POISSON=E(MTYPE)/(2.0*G(MTYPE))-1.0
        KFACTOR=10.0*(1.0+POISSON)/(12.0+11.0*POISSON)

        S(1,1)=E(MTYPE)*AREA(MTYPE)/XL
        S(2,2)=G(MTYPE)*AREA(MTYPE)*KFACTOR/XL
        S(3,3)=G(MTYPE)*AREA(MTYPE)*KFACTOR/XL
        S(4,4)=G(MTYPE)*J_Z(MTYPE)/XL
        S(5,5)=E(MTYPE)*I_Y(MTYPE)/XL+G(MTYPE)*AREA(MTYPE)*KFACTOR*XL/4.0
        S(6,6)=E(MTYPE)*I_Z(MTYPE)/XL+G(MTYPE)*AREA(MTYPE)*KFACTOR*XL/4.0
        S(7,7)=E(MTYPE)*AREA(MTYPE)/XL
        S(8,8)=G(MTYPE)*AREA(MTYPE)*KFACTOR/XL
        S(9,9)=G(MTYPE)*AREA(MTYPE)*KFACTOR/XL
        S(10,10)=G(MTYPE)*J_Z(MTYPE)/XL
        S(11,11)=E(MTYPE)*I_Y(MTYPE)/XL+G(MTYPE)*AREA(MTYPE)*KFACTOR*XL/4.0
        S(12,12)=E(MTYPE)*I_Z(MTYPE)/XL+G(MTYPE)*AREA(MTYPE)*KFACTOR*XL/4.0
        S(5,3)=-G(MTYPE)*AREA(MTYPE)*KFACTOR/2.0
        S(6,2)=G(MTYPE)*AREA(MTYPE)*KFACTOR/2.0
        S(7,1)=-E(MTYPE)*AREA(MTYPE)/XL
        S(8,2)=-G(MTYPE)*AREA(MTYPE)*KFACTOR/XL
        S(8,6)=-G(MTYPE)*AREA(MTYPE)*KFACTOR/2.0
        S(9,3)=-G(MTYPE)*AREA(MTYPE)*KFACTOR/XL
        S(9,5)=G(MTYPE)*AREA(MTYPE)*KFACTOR/2.0
        S(10,4)=-G(MTYPE)*J_X(MTYPE)/XL
        S(11,3)=-G(MTYPE)*AREA(MTYPE)*KFACTOR/2.0
        S(11,5)=-E(MTYPE)*I_Y(MTYPE)/XL+G(MTYPE)*AREA(MTYPE)*KFACTOR*XL/4.0
        S(11,9)=G(MTYPE)*AREA(MTYPE)*KFACTOR/2.0
        S(12,2)=G(MTYPE)*AREA(MTYPE)*KFACTOR/2.0
        S(12,6)=-E(MTYPE)*I_Z(MTYPE)/XL+G(MTYPE)*AREA(MTYPE)*KFACTOR*XL/4.0
        S(12,8)=-G(MTYPE)*AREA(MTYPE)*KFACTOR/2.0

        DO I=1,11
           DO J=I+1,12
               S(I,J)=S(J,I)
           ENDDO
        ENDDO

        S=MATMUL(TRANSPOSE(T),MATMUL(S,T))

     !GET THE DISPLACEMENT FOR ELEMENT N
        DO L=1,12
           UELE(L)=0.0
        ENDDO
        DO L=1,12
           I=LM(L,N)
           IF (I .GT. 0)THEN
               UELE(L)=U(I)
           ENDIF
        ENDDO

      !CALCULATE FORCE AND MOMENT OF THE ELEMENT

        DO I=1,12
           SIGMA=0.0
           DO J=1,12
               SIGMA=SIGMA+UELE(J)*S(I,J)
           ENDDO
           FORCE(I)=SIGMA
        ENDDO

       !PRINT THE RESULTS FOR POSTPROCESSING

       WRITE (IOUT,"(1X,I10,1X,12(E12.6,1X))") N,(FORCE(I),I=1,12)
       GaussianCollection(:,N) = 0.5*(XYZ(4:6,N)+XYZ(1:3,N))
       StressCollection((/1,2,3,5,6,4/),N) = 0.5*(FORCE(7:12)-FORCE(1:6))
     END DO
     call PostProcessor(NPAR(1), 3, XYZ, &
                        Node, 1, GaussianCollection, StressCollection, U)
  ELSE
     STOP "*** ERROR *** Invalid IND value."
  END IF

 END SUBROUTINE TIMOELE

max_stars_repo_name: JAMMIASHOK/OpenSTAP


text:
      SUBROUTINE SBINIT
C
      IMPLICIT REAL*8(A-H,O-Z)
C
      COMMON/ACFILT/ACPARM(20),IACFLT(20)
C
      COMMON/ADSTAT/DER(150),DEP(150)
C
      COMMON/CONSTS/PI,TWOPI,RADIAN
C
      COMMON/CRASTR/DELA(3),TAUA(4,3),ANG20(3),ADD0(3),TC(3),TTAB(4,3)
C
      COMMON/CSECBD/SECM,SECI(3,3),ZBAR2(3),YI02(3),ZZ2(3,3),ZZ2Y(3,3)
     1             ,YZ2(3),YZT(3),ZB(3),YIZ02(3)
C
      COMMON/ISECBD/I2BDY,NDOF2,IRAST,NDEP2,IARST(3),IRSCY(3)
C
      COMMON/PWHEEL/XMOMIN(3),DVMOM(3),VMOM(3),VW(3)
C
      COMMON/SECBDW/GAMMA2(3,6),DELTA2(3,3),DELT2(3,3),DELT2I(3,3)
C
      COMMON/SECICS/GAM20,ALP20,BET20,B0(3,3),GAM2I,ALP2I,BET2I
     1             ,OM2I(3),SBUP(2),SBDN(2)
C
      COMMON/XIN4  /UP(150),DN(150),BNDS(22)
C
C
      IF(I2BDY.EQ.0) GO TO 100
C
      GAM=GAM20*RADIAN
      ALP=ALP20*RADIAN
      BET=BET20*RADIAN
C
      CALL DTR312(GAM,ALP,BET,B0)
C
      ZZ2(1,1)=0.5D0*(SECI(2,2)+SECI(3,3)-SECI(1,1))
     1         +SECM*ZBAR2(1)*ZBAR2(1)
      ZZ2(2,2)=0.5D0*(SECI(3,3)+SECI(1,1)-SECI(2,2))
     1         +SECM*ZBAR2(2)*ZBAR2(2)
      ZZ2(3,3)=0.5D0*(SECI(1,1)+SECI(2,2)-SECI(3,3))
     1         +SECM*ZBAR2(3)*ZBAR2(3)
      ZZ2(1,2)=SECI(1,2)+SECM*ZBAR2(1)*ZBAR2(2)
      ZZ2(1,3)=SECI(1,3)+SECM*ZBAR2(1)*ZBAR2(3)
      ZZ2(2,3)=SECI(2,3)+SECM*ZBAR2(2)*ZBAR2(3)
      ZZ2(2,1)=ZZ2(1,2)
      ZZ2(3,1)=ZZ2(1,3)
      ZZ2(3,2)=ZZ2(2,3)
C
      IF(NDOF2.EQ.0) GO TO 100
C
      GAM2=GAM2I*RADIAN
      ALP2=ALP2I*RADIAN
      BET2=BET2I*RADIAN
C
      SALP=DSIN(ALP2)
      CALP=DCOS(ALP2)
      SBET=DSIN(BET2)
      CBET=DCOS(BET2)
      GAM2D=0.0D0
C
      IF(NDOF2.EQ.1) GO TO 20
      GAM2D=(OM2I(3)*CBET-OM2I(1)*SBET)*RADIAN/CALP
C
   20 CONTINUE
C
      BET2D=OM2I(2)*RADIAN-GAM2D*SALP
      I0=NDEP2
      I1=I0+NDOF2
      DEP(I0)=BET2
      DEP(I1)=BET2D
      UP(I0)=SBUP(1)
      DN(I0)=SBDN(1)
      UP(I1)=SBUP(2)
      DN(I1)=SBDN(2)
C
      IF(NDOF2.EQ.1) GO TO 80
C
      I0=I0+1
      I1=I1+1
      DEP(I0)=GAM2
      DEP(I1)=GAM2D
      UP(I0)=SBUP(1)
      DN(I0)=SBDN(1)
      UP(I1)=SBUP(2)
      DN(I1)=SBDN(2)
C
      IF(NDOF2.EQ.2) GO TO 80
C
      ALP2D=(OM2I(1)*CBET+OM2I(3)*SBET)*RADIAN
      I0=I0+1
      I1=I1+1
      DEP(I0)=ALP2
      DEP(I1)=ALP2D
      UP(I0)=SBUP(1)
      DN(I0)=SBDN(1)
      UP(I1)=SBUP(2)
      DN(I1)=SBDN(2)
C
   80 CONTINUE
C
      DELTA2(1,1)=ZZ2(2,2)+ZZ2(3,3)
      DELTA2(2,2)=ZZ2(3,3)+ZZ2(1,1)
      DELTA2(3,3)=ZZ2(1,1)+ZZ2(2,2)
      DELTA2(1,2)=-ZZ2(1,2)
      DELTA2(1,3)=-ZZ2(1,3)
      DELTA2(2,3)=-ZZ2(2,3)
      DELTA2(2,1)=DELTA2(1,2)
      DELTA2(3,1)=DELTA2(1,3)
      DELTA2(3,2)=DELTA2(2,3)
C
C
  100 CONTINUE
C
      DO 105 I=1,3
      VW(I)=VMOM(I)
  105 CONTINUE
C
      IF(IRAST.EQ.0) GO TO 150
C
      DO 120 I=1,3
      ITEST=IARST(I)
      IF(ITEST.EQ.0) GO TO 120
      IF(ITEST.EQ.2) GO TO 110
C
C     TYPE 1 CYCLE
C
      TC(I)=2.0D0*(TAUA(2,I)+2.0D0*TAUA(3,I)+TAUA(4,I))
      TTAB(1,I)=TAUA(2,I)/2.0D0
      TTAB(2,I)=TTAB(1,I)+TAUA(3,I)
      TTAB(3,I)=TTAB(2,I)+TAUA(4,I)
      TTAB(4,I)=TTAB(3,I)+TAUA(3,I)
      ADD0(I)=DELA(I)*PI/(2.0D0*TAUA(3,I)*(TAUA(3,I)+TAUA(4,I)))
C
      GO TO 120
C
  110 CONTINUE
C
C     TYPE 2 CYCLE
C
      TC(I)=TAUA(2,I)+TAUA(3,I)
      TTAB(1,I)=TAUA(2,I)
      TTAB(2,I)=TTAB(1,I)+TAUA(3,I)
      ADD0(I)=DELA(I)*TWOPI/TAUA(2,I)**2
C
  120 CONTINUE
C
  150 CONTINUE
C
      IF(IACFLT(1).EQ.0) RETURN
C
      I0=IACFLT(6)
      I1=I0+1
      UP(I0)=ACPARM(19)
      DN(I0)=ACPARM(20)
      UP(I1)=ACPARM(19)
      DN(I1)=ACPARM(20)
      IF(IACFLT(4).EQ.0) GO TO 155
      DEP(I0)=ACPARM(17)
      DEP(I1)=ACPARM(18)
  155 CONTINUE
C
C
      RETURN
C
      END

ext: for