text:
module iterator_product_test_module
    use :: assert_module, only : assert
    use :: iterator_product_module, only : iterator_product
    use :: iterator_module, only : iterator
    use :: iterator_constraint_wrapper_module, only : iterator_constraint_wrapper
    use :: pleq_constraint_module, only : pleq_constraint
    implicit none
    private

    public :: iterator_product_test

    type :: iterator_product_test
    contains
        procedure :: run => run
        procedure :: cleanup => cleanup
        procedure :: clear => clear
    end type iterator_product_test

    interface iterator_product_test
        module procedure constructor
    end interface iterator_product_test

contains
    function constructor() result(this)
        type(iterator_product_test) :: this

        call this%clear()
    end function constructor

    subroutine run(this, assertion)
        class(iterator_product_test), intent(in) :: this
        type(assert), intent(inout) :: assertion

        integer, dimension(:,:), allocatable :: combined, res
        type(iterator) :: aniterator
        type(iterator), dimension(:), allocatable :: iterators
        type(iterator_constraint_wrapper), dimension(1) :: constraints
        integer :: rank, number_of_states

        combined = get_combined([iterator(1,3)])
        allocate(res(3,1))
        res(:,1) = [1,2,3]
        call assertion%is_equal("iterator_product::rank1_size1", 3, size(combined))
        call assertion%equal("iterator_product::rank1_combined1", all(combined == res) )
        deallocate(res)

        combined = get_combined([iterator(2,3)])
        allocate(res(2,1))
        res(:,1) = [2,3]
        call assertion%is_equal("iterator_product::rank1_size2", 2, size(combined))
        call assertion%equal("iterator_product::rank1_combined2", all(combined == res) )
        deallocate(res)

        combined = get_combined([iterator(1,2), iterator(3,4)])
        allocate(res(4,2))
        res(:,1) = [1, 1, 2, 2]
        res(:,2) = [3, 4, 3, 4]
        call assertion%is_equal("iterator_product::rank2", 2, size(combined, 2))
        call assertion%is_equal("iterator_product::rank2_size1", 4, size(combined, 1))
        call assertion%equal("iterator_product::rank2_combined1", all(combined == res) )
        deallocate(res)

        combined = get_combined_column_major([iterator(1,2), iterator(3,4)])
        allocate(res(4,2))
        res(:,1) = [1, 2, 1, 2]
        res(:,2) = [3, 3, 4, 4]
        call assertion%is_equal("iterator_product::rank2 column major", 2, size(combined, 2))
        call assertion%is_equal("iterator_product::rank2_size1 column major", 4, size(combined, 1))
        call assertion%equal("iterator_product::rank2_combined1 column major", all(combined == res) )
        deallocate(res)

        number_of_states = 12
        aniterator = iterator(1, number_of_states)
        constraints = [ iterator_constraint_wrapper(pleq_constraint()) ]

        rank = 1
        iterators = [aniterator]
        combined = get_combined(iterators)
        call assertion%is_equal("iterator_product::p1::number_of_configs",  &
            size(combined, 1), number_of_states**rank)
        combined = get_combined(iterators, constraints)
        call assertion%is_equal("iterator_product::p1_pleq::number_of_configs",  &
            size(combined, 1), get_triangular_number(number_of_states,rank))

        rank = 2
        iterators = [aniterator, aniterator]

        combined = get_combined(iterators)
        call assertion%is_equal("iterator_product::p2::number_of_configs",  &
            size(combined, 1), number_of_states**rank)
        combined = get_combined(iterators, constraints)
        call assertion%is_equal("iterator_product::p2_pleq::number_of_configs",  &
            size(combined, 1), get_triangular_number(number_of_states,rank))

        rank = 3
        iterators = [aniterator, aniterator, aniterator]

        combined = get_combined(iterators)
        call assertion%is_equal("iterator_product::p3::number_of_configs",  &
            size(combined, 1), number_of_states**rank)
        combined = get_combined(iterators, constraints)
        call assertion%is_equal("iterator_product::p3_pleq::number_of_configs",  &
            size(combined, 1), get_triangular_number(number_of_states,rank))

        rank = 4
        iterators = [aniterator, aniterator, aniterator, aniterator]

        combined = get_combined(iterators)
        call assertion%is_equal("iterator_product::p4::number_of_configs",  &
            size(combined, 1), number_of_states**rank)
        combined = get_combined(iterators, constraints)
        call assertion%is_equal("iterator_product::p4_pleq::number_of_configs",  &
            size(combined, 1), get_triangular_number(number_of_states,rank))

        rank = 5
        iterators = [aniterator, aniterator, aniterator, aniterator, aniterator]

        combined = get_combined(iterators)
        call assertion%is_equal("iterator_product::p5::number_of_configs",  &
            size(combined, 1), number_of_states**rank)
        combined = get_combined(iterators, constraints)
        call assertion%is_equal("iterator_product::p5_pleq::number_of_configs",  &
            size(combined, 1), get_triangular_number(number_of_states,rank))
    end subroutine run

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

        call this%clear()
    end subroutine cleanup

    subroutine clear(this)
        class(iterator_product_test), intent(inout) :: this
    end subroutine clear

    function get_combined(iterators, constraints) result(combined)
        type(iterator), dimension(:) :: iterators
        type(iterator_constraint_wrapper), dimension(:), intent(in), optional :: constraints
        integer, dimension(:,:), allocatable :: combined

        type(iterator_product) :: iproduct
        type(iterator) :: aniterator

        iproduct = iterator_product(constraints)
        aniterator = iproduct%get_product_row_major(iterators)
        combined = aniterator%get_all_items()
    end function get_combined

    function get_combined_column_major(iterators, constraints) result(combined)
        type(iterator), dimension(:) :: iterators
        type(iterator_constraint_wrapper), dimension(:), intent(in), optional :: constraints
        integer, dimension(:,:), allocatable :: combined

        type(iterator_product) :: iproduct
        type(iterator) :: aniterator

        iproduct = iterator_product(constraints)
        aniterator = iproduct%get_product_column_major(iterators)
        combined = aniterator%get_all_items()
    end function get_combined_column_major

    function get_triangular_number(nstates, nranks) result(triangular_number)
        integer, intent(in) :: nstates, nranks
        integer :: triangular_number

        integer :: idx

        triangular_number = nstates
        do idx = 1, nranks-1
            triangular_number = triangular_number*(nstates+idx)
        end do

        do idx = nranks, 2, -1
            triangular_number = triangular_number/idx
        end do
    end function get_triangular_number

end module iterator_product_test_module

max_stars_repo_name: cheshyre/ntcl-util

max_issues_repo_name: cheshyre/ntcl-util

max_forks_repo_name: cheshyre/ntcl-util


text:
c     ****************************************************************
c     *                                                              *
c     *                      subroutine ouext2                       *
c     *                                                              *
c     *                       written by : kck                       *
c     *                                                              *
c     *                   last modified : 6/30/2019 rhd              *
c     *                                                              *
c     *     this subroutine computes derived stress/strain values    *
c     *     after the primary values have been averaged              *
c     *                                                              *
c     ****************************************************************
c
      subroutine ouext2 ( results, nrowd, numele, do_stresses )
      implicit none
c
      integer :: nrowd, numele
      double precision :: results(nrowd,*)
      logical :: do_stresses
c
      if( do_stresses ) then
         call princ_inv_stress ( results, nrowd, numele )
         call princ_stress ( results, nrowd, numele )
         call yield_function( results, nrowd, numele )
      else
         call princ_inv_strain ( results, nrowd, numele )
         call princ_strain ( results, nrowd, numele )
         call equiv_strain( results, nrowd, numele )
      end if
c
      return
      end

c     ****************************************************************
c     *                                                              *
c     *                      subroutine princ_inv_strain             *
c     *                                                              *
c     *                       written by : kck                       *
c     *                                                              *
c     *                   last modified : 10/04/94                   *
c     *                                   03/05/95 kck               *
c     *                                   06/11/97 rhd               *
c     *                                                              *
c     *     this subroutine computes principal invariants            *
c     *     at the nodes/elem after the primary values have          *
c     *     been averaged                                            *
c     *                                                              *
c     ****************************************************************
c
      subroutine princ_inv_strain( results, nrowd, numele )
      implicit integer (a-z)
      double precision ::results(nrowd,*), half, t1, t2, t3
      data  half  / 0.5d0  /
c
       do i = 1, numele
         t1 =  half * results(i,4)
         t2 =  half * results(i,5)
         t3 =  half * results(i,6)
         results(i,8) = results(i,1) +  results(i,2) +
     &                        results(i,3)
         results(i,9) = - t1 * t1 - t2 * t2 - t3 * t3 +
     &                        results(i,1) * results(i,2) +
     &                        results(i,2) * results(i,3) +
     &                        results(i,1) * results(i,3)
         results(i,10) =
     &     results(i,1) * ( results(i,2) * results(i,3) -
     &              t2 * t2 ) -  t1 * ( t1 * results(i,3) -
     &              t2 * t3 ) + t3 * ( t1 * t2 -
     &              results(i,2) * t3 )
       end do
c
       return
       end
c     ****************************************************************
c     *                                                              *
c     *                      subroutine princ_inv_stress             *
c     *                                                              *
c     *                       written by : kck                       *
c     *                                                              *
c     *                   last modified : 10/04/94                   *
c     *                                   03/05/95 kck               *
c     *                                   06/15/97 rhd               *
c     *                                                              *
c     *     this subroutine computes principal invariants            *
c     *     at the nodes/elements after the primary values have      *
c     *     been averaged                                            *
c     *                                                              *
c     ****************************************************************
c
      subroutine princ_inv_stress( results, mxvl, numele )
      implicit integer (a-z)
      double precision :: results(mxvl,*)

       do i = 1, numele
         results(i,12) = results(i,1) +  results(i,2) +
     &                        results(i,3)
         results(i,13) = results(i,1) * results(i,2) +
     &                        results(i,2) * results(i,3) +
     &                        results(i,1) * results(i,3) -
     &                        results(i,4) * results(i,4) -
     &                        results(i,5) * results(i,5) -
     &                        results(i,6) * results(i,6)
         results(i,14) = results(i,1) *
     &            ( results(i,2) * results(i,3) -
     &              results(i,5) * results(i,5) ) -
     &                        results(i,4) *
     &            ( results(i,4) * results(i,3) -
     &              results(i,5) * results(i,6) ) +
     &                        results(i,6) *
     &            ( results(i,4) * results(i,5) -
     &              results(i,2) * results(i,6) )
       end do
c
       return
       end

c     ****************************************************************
c     *                                                              *
c     *                      subroutine princ_strain                 *
c     *                                                              *
c     *                       written by : kck                       *
c     *                                                              *
c     *                   last modified : 10/04/94                   *
c     *                                   03/05/95 kck               *
c     *                                   06/10/97 rhd               *
c     *                                                              *
c     *     this subroutine computes principal strain values         *
c     *     at the nodes/elem after the primary values have          *
c     *     been averaged                                            *
c     *                                                              *
c     ****************************************************************
c
      subroutine princ_strain( results, nrowd, numele )
      implicit integer (a-z)
      include 'param_def'
      double precision :: results(nrowd,*)
c
c                    locally allocated
c
      double precision
     &  temp_strain(nstr), wk(ndim), ev(nstr), evec(ndim,ndim), half
c
      data  half / 0.5d0 /
c
c
c        calculate the principal strains and there direction cosines by
c        passing off strains into dummy array, then calling an eigenvalue
c        eigenvector routine, and then passing the values from this routine
c        back to the appropriate places within results
c
       do i = 1, numele
          temp_strain(1) = results(i,1)
          temp_strain(2) = results(i,4) * half
          temp_strain(3) = results(i,2)
          temp_strain(4) = results(i,6) * half
          temp_strain(5) = results(i,5) * half
          temp_strain(6) = results(i,3)
          call ou3dpr( temp_strain, ndim, 1, ev, evec, ndim, wk, ier )
          results(i,11) = ev(1)
          results(i,12) = ev(2)
          results(i,13) = ev(3)
          results(i,14) = evec(1,1)
          results(i,15) = evec(2,1)
          results(i,16) = evec(3,1)
          results(i,17) = evec(1,2)
          results(i,18) = evec(2,2)
          results(i,19) = evec(3,2)
          results(i,20) = evec(1,3)
          results(i,21) = evec(2,3)
          results(i,22) = evec(3,3)
      end do
      return
      end

c     ****************************************************************
c     *                                                              *
c     *                      subroutine princ_stress                 *
c     *                                                              *
c     *                       written by : kck                       *
c     *                                                              *
c     *                   last modified : 10/04/94                   *
c     *                                   03/05/95 kck               *
c     *                                   06/10/97 rhd               *
c     *                                                              *
c     *     this subroutine computes principal stress values         *
c     *     at the nodes/elem after the primary values have          *
c     *     been averaged                                            *
c     *                                                              *
c     ****************************************************************
c
      subroutine princ_stress( results, nrowd, numele )
      implicit integer (a-z)
      include 'param_def'
      double precision :: results(nrowd,*)
c
c                    locally allocated
c
      double precision  ::
     &  temp_stress(nstr), wk(ndim), ev(nstr), evec(ndim,ndim)
c
c
c        calculate the principal stresses and there direction cosines by
c        passing off stresses into dummy array, then calling an eigenvalue
c        eigenvector routine, and then passing the values from this routine
c        back to the appropriate places within results
c
       do i = 1, numele
          temp_stress(1) = results(i,1)
          temp_stress(2) = results(i,4)
          temp_stress(3) = results(i,2)
          temp_stress(4) = results(i,6)
          temp_stress(5) = results(i,5)
          temp_stress(6) = results(i,3)
          call ou3dpr( temp_stress, ndim, 1, ev, evec, ndim, wk, ier )
          results(i,15) = ev(1)
          results(i,16) = ev(2)
          results(i,17) = ev(3)
          results(i,18) = evec(1,1)
          results(i,19) = evec(2,1)
          results(i,20) = evec(3,1)
          results(i,21) = evec(1,2)
          results(i,22) = evec(2,2)
          results(i,23) = evec(3,2)
          results(i,24) = evec(1,3)
          results(i,25) = evec(2,3)
          results(i,26) = evec(3,3)
      end do
      return
      end

c     ****************************************************************
c     *                                                              *
c     *                      subroutine yield_function               *
c     *                                                              *
c     *                       written by : kck                       *
c     *                                                              *
c     *                   last modified : 10/04/94                   *
c     *                                   03/05/95 kck               *
c     *                                                              *
c     *     this subroutine computes mises yield function value      *
c     *     at the nodes/elements after the primary values have      *
c     *     been averaged                                            *
c     *                                                              *
c     ****************************************************************
c
      subroutine yield_function ( results, maxnum, numele )
      implicit integer (a-z)
      double precision :: results(maxnum,*), iroot2, six
      data iroot2, six / 0.70711d0, 6.0d0 /
c
      do i = 1, numele
        results(i,8) = sqrt(
     &     ( results(i,1) - results(i,2) ) ** 2 +
     &     ( results(i,2) - results(i,3) ) ** 2 +
     &     ( results(i,1) - results(i,3) ) ** 2 +
     &   six*( results(i,4)**2 + results(i,5)**2+
     &         results(i,6)**2 ) )*iroot2
c
      end do
c
      return
      end
c     ****************************************************************
c     *                                                              *
c     *                      subroutine equiv_strain                 *
c     *                                                              *
c     *                       written by : kck                       *
c     *                                                              *
c     *                   last modified : 10/04/94                   *
c     *                                   06/02/95 kck               *
c     *                                                              *
c     *     this subroutine computes mises yield function value      *
c     *     at the nodes/elem after the primary values have          *
c     *     been averaged                                            *
c     *                                                              *
c     ****************************************************************
c
      subroutine equiv_strain ( results, maxnum, numele )
      implicit integer (a-z)
      double precision
     &     results(maxnum,*),  root23, onep5
      data root23, onep5 / 0.471404d0, 1.5d0 /
c
      do i = 1, numele
        results(i,7) = root23 * sqrt(
     &     ( results(i,1) - results(i,2) ) ** 2 +
     &     ( results(i,2) - results(i,3) ) ** 2 +
     &     ( results(i,1) - results(i,3) ) ** 2 +
     &   onep5 * ( results(i,4)**2 + results(i,5)**2+
     &         results(i,6)**2 ) )
      end do
c
      return
      end
c *******************************************************************
c *                                                                 *
c *      element output service routine -- ou3dpr                   *
c *                                                                 *
c *******************************************************************
c
c
      subroutine ou3dpr( a,n,ijob,d,z,iz,wk,ier )
      implicit double precision (a-h,o-z)
c   purpose             - eigenvalues and (optionally) eigenvectors of
c                           a real symmetric matrix in symmetric
c                           storage mode
c
c   arguments    a      - input real symmetric matrix of order n,
c                           stored in symmetric storage mode,
c                           whose eigenvalues and eigenvectors
c                           are to be computed. input a is
c                           destroyed if ijob is equal to 0 or 1.
c                n      - input order of the matrix a.
c                ijob   - input option parameter, when
c                           ijob = 0, compute eigenvalues only
c                           ijob = 1, compute eigenvalues and eigen-
c                             vectors.
c                           ijob = 2, compute eigenvalues, eigenvectors
c                             and performance index.
c                           ijob = 3, compute performance index only.
c                           if the performance index is computed, it is
c                           returned in wk(1). the routines have
c                           performed (well, satisfactorily, poorly) if
c                           wk(1) is (less than 1, between 1 and 100,
c                           greater than 100).
c                d      - output vector of length n,
c                           containing the eigenvalues of a.
c                z      - output n by n matrix containing
c                           the eigenvectors of a.
c                           the eigenvector in column j of z corres-
c                           ponds to the eigenvalue d(j).
c                           if ijob = 0, z is not used.
c                iz     - input row dimension of matrix z exactly as
c                           specified in the dimension statement in the
c                           calling program.
c                wk     - work area, the length of wk depends
c                           on the value of ijob, when
c                           ijob = 0, the length of wk is at least n.
c                           ijob = 1, the length of wk is at least n.
c                           ijob = 2, the length of wk is at least
c                             n(n+1)/2+n.
c                           ijob = 3, the length of wk is at least 1.
c                ier    - error parameter (output)
c                         terminal error
c                           ier = 128+j, indicates that ourt2s failed
c                             to converge on eigenvalue j. eigenvalues
c                             and eigenvectors 1,...,j-1 have been
c                             computed correctly, but the eigenvalues
c                             are unordered. the performance index
c                             is set to 1000.0
c                         warning error (with fix)
c                           ier = 66, indicates ijob is less than 0 or
c                             ijob is greater than 3. ijob set to 1.
c                           ier = 67, indicates ijob is not equal to
c                             zero, and iz is less than the order of
c                             matrix a. ijob is set to zero.
c
c                                  specifications for arguments
      integer            n,ijob,iz,ier
      double precision
     &                   a(*),d(*),wk(*),z(iz,*)
c                                  specifications for local variables
      integer            jer,na,nd,iiz,ibeg,il,kk,lk,i,j,k,l
      double precision
     &                   anorm,asum,pi,sumz,sumr,an,s,ten,rdelp,zero,
     &                   one,thous
      data               rdelp /0.745058d-08/
      data               zero,one/0.d0,1.d0/,ten/10.d0/,thous/1000.d0/
c                                  initialize error parameters
c                                  first executable statement
      ier = 0
      jer = 0
      if (ijob .ge. 0 .and. ijob .le. 3) go to 5
c                                  warning error - ijob is not in the
c                                    range
      ier = 66
      ijob = 1
      go to 10
    5 if (ijob .eq. 0) go to 20
   10 if (iz .ge. n) go to 15
c                                  warning error - iz is less than n
c                                    eigenvectors can not be computed,
c                                    ijob set to zero
      ier = 67
      ijob = 0
   15 if (ijob .eq. 3) go to 65
   20 na = (n*(n+1))/2
      if (ijob .ne. 2) go to 35
      do 30 i=1,na
         wk(i) = a(i)
   30 continue
c                                  save input a if ijob = 2
   35 nd = 1
      if (ijob .eq. 2) nd = na+1
c                                  reduce a to symmetric tridiagonal
c                                    form
      call ouhous (a,n,d,wk(nd),wk(nd))
      iiz = 1
      if (ijob .eq. 0) go to 50
      iiz = iz
c                                  set z to the identity matrix
      do 45 i=1,n
         do 40 j=1,n
            z(i,j) = zero
   40    continue
         z(i,i) = one
   45 continue
c                                  compute eigenvalues and eigenvectors
   50 call ourt2s (d,wk(nd),n,z,iiz,jer)
      if (ijob .eq. 0) go to 9000
      if (jer .gt. 128) go to 55
c                                  back transform eigenvectors
      call ouobks (a,n,1,n,z,iz)
   55 if (ijob .le. 1) go to 9000
c                                  move input matrix back to a
      do 60 i=1,na
         a(i) = wk(i)
   60 continue
      wk(1) = thous
      if (jer .ne. 0) go to 9000
c                                  compute 1 - norm of a
   65 anorm = zero
      ibeg = 1
      do 75 i=1,n
         asum = zero
         il = ibeg
         kk = 1
         do 70 l=1,n
            asum =asum+abs(a(il))
            if (l .ge. i) kk = l
            il = il+kk
   70    continue
         anorm = dmax1(anorm,asum)
         ibeg = ibeg+i
   75 continue
      if (anorm .eq. zero) anorm = one
c                                  compute performance index
      pi = zero
      do 90 i=1,n
         ibeg = 1
         s = zero
         sumz = zero
         do 85 l=1,n
            lk = ibeg
            kk = 1
            sumz = sumz+abs(z(l,i))
            sumr = -d(i)*z(l,i)
            do 80 k=1,n
               sumr = sumr+a(lk)*z(k,i)
               if (k .ge. l) kk = k
               lk = lk+kk
   80       continue
            s = s+abs(sumr)
            ibeg = ibeg+l
   85    continue
         if (sumz .eq. zero) go to 90
         pi = dmax1(pi,s/sumz)
   90 continue
      an = n
      pi = pi/(anorm*ten*an*rdelp)
      wk(1) = pi
 9000 continue
      if ( ier .ne. 0 ) return
      if (jer .eq. 0) go to 9005
      ier = jer
      return
 9005 return
      end
c
c
c *******************************************************************
c *                                                                 *
c *      element output service routine -- ouobks                   *
c *                                                                 *
c *******************************************************************
c
c
      subroutine ouobks( a,n,m1,m2,z,iz )
      implicit double precision (a-h,o-z)
c
c   purpose             - back transformation to form the eigenvectors
c                           of the original symmetric matrix from the
c                           eigenvectors of the tridiagonal matrix
c
c   arguments    a      - the array contains the details of the house-
c                           holder reduction of the original matrix a
c                           as generated by routine ouhous.(input)
c                n      - order of the real symmetric matrix.(input)
c                m1     - m1 and m2 are two input scalars such that
c                           eigenvectors m1 to m2 of the tridiagonal
c                           matrix a have been found and normalized
c                           according to the euclidean norm.
c                m2     - see above - m1
c                z      - a two dimensional array of size n x (m2-m1+1)
c                           which contains eigenvectors m1 to m2 of
c                           tridiagonal matrix t, normalized according
c                           to euclidean norm. input z can be produced
c                           by routine ourt2s, the resultant
c                           matrix overwrites the input z.(input/output)
c                iz     - row dimension of matrix z exactly as
c                           specified in the dimension statement in the
c                           calling program. (input)
c
c
      double precision
     &                   zero
      dimension          a(*),z(iz,*)
      data zero /0.0/

c                                  first executable statement
      if (n.eq.1) go to 30
      do 25 i=2,n
         l = i-1
         ia = (i*l)/2
         h = a(ia+i)
         if (h.eq.zero) go to 25
c                                  derives eigenvectors m1 to m2 of
c                                  the original matrix from eigenvectors
c                                  m1 to m2 of the symmetric
c                                  tridiagonal matrix
         do 20 j = m1,m2
            s = zero
            do 10 k = 1,l
               s = s+a(ia+k)*z(k,j)
   10       continue
            s = s/h
            do 15 k=1,l
               z(k,j) = z(k,j)-s*a(ia+k)
   15       continue
   20    continue
   25 continue
   30 return
      end
c *******************************************************************
c *                                                                 *
c *      element output service routine -- ouhous                   *
c *                                                                 *
c *******************************************************************
c
c
      subroutine ouhous( a,n,d,e,e2 )
      implicit double precision (a-h,o-z)
c
c   purpose             - reduction of a symmetric matrix to symmetric
c                           tridiagonal form using a householder
c                           reduction
c
c   arguments    a      - the given n x n, real symmetric matrix a,
c                           where a is stored in symmetric storage mode.
c                           the input a is replaced by the details of
c                           the householder reduction of a.
c                n      - input order of a and the length of d, e, and
c                           e2.
c                d      - the output array of length n, giving the
c                           diagonal elements of the tridiagonal matrix.
c                e      - the output array of length n, giving the sub-
c                           diagonal in the last (n-1) elements, e(1) is
c                           set to zero.
c                e2     - output array of length n.  e2(i) = e(i)**2.
c
c
      dimension          a(*),d(n),e(n),e2(n)
      double precision
     &                   a,d,e,e2,zero,h,scale,one,scale1,f,g,hh
      data               zero/0.d0/,one/1.d0/
c                                  first executable statement
      np1 = n+1
      nn = (n*np1)/2-1
      nbeg = nn+1-n
      do 70 ii = 1,n
         i = np1-ii
         l = i-1
         h = zero
         scale = zero
         if (l .lt. 1) go to 10
c                                  scale row (algol tol then not needed)
         nk = nn
         do 5 k = 1,l
            scale = scale+abs(a(nk))
            nk = nk-1
    5    continue
         if (scale .ne. zero) go to 15
   10    e(i) = zero
         e2(i) = zero
         go to 65
   15    nk = nn
         scale1 = one/scale
         do 20 k = 1,l
            a(nk) = a(nk)*scale1
            h = h+a(nk)*a(nk)
            nk = nk-1
   20    continue
         e2(i) = scale*scale*h
         f = a(nn)
         g = -sign(sqrt(h),f)
         e(i) = scale*g
         h = h-f*g
         a(nn) = f-g
         if (l .eq. 1) go to 55
         f = zero
         jk1 = 1
         do 40 j = 1,l
            g = zero
            ik = nbeg+1
            jk = jk1
c                                  form element of a*u
            do 25 k = 1,j
               g = g+a(jk)*a(ik)
               jk = jk+1
               ik = ik+1
   25       continue
            jp1 = j+1
            if (l .lt. jp1) go to 35
            jk = jk+j-1
            do 30 k = jp1,l
               g = g+a(jk)*a(ik)
               jk = jk+k
               ik = ik+1
   30       continue
c                                  form element of p
   35       e(j) = g/h
            f = f+e(j)*a(nbeg+j)
            jk1 = jk1+j
   40    continue
         hh = f/(h+h)
c                                  form reduced a
         jk = 1
         do 50 j = 1,l
            f = a(nbeg+j)
            g = e(j)-hh*f
            e(j) = g
            do k = 1,j
               a(jk) = a(jk)-f*e(k)-g*a(nbeg+k)
               jk = jk+1
            enddo
   50    continue
   55    do 60 k = 1,l
            a(nbeg+k) = scale*a(nbeg+k)
   60    continue
   65    d(i) = a(nbeg+i)
         a(nbeg+i) = h*scale*scale
         nbeg = nbeg-i+1
         nn = nn-i
   70 continue
      return
      end
c *******************************************************************
c *                                                                 *
c *      element output service routine -- ourt2s
c *                                                                 *
c *******************************************************************
c
c
      subroutine ourt2s(d,e,n,z,iz,ier)
      implicit double precision (a-h,o-z)
c   purpose             - eigenvalues and (optionally) eigenvectors of
c                           a symmetric tridiagonal matrix using the
c                           ql method.
c
c   arguments    d      - on input, the vector d of length n contains
c                           the diagonal elements of the symmetric
c                           tridiagonal matrix t.
c                           on output, d contains the eigenvalues of
c                           t in ascending order.
c                e      - on input, the vector e of length n contains
c                           the sub-diagonal elements of t in position
c                           2,...,n. on output, e is destroyed.
c                n      - order of tridiagonal matrix t.(input)
c                z      - on input, z contains the identity matrix of
c                           order n.
c                           on output, z contains the eigenvectors
c                           of t. the eigenvector in column j of z
c                           corresponds to the eigenvalue d(j).
c                iz     - input row dimension of matrix z exactly as
c                           specified in the dimension statement in the
c                           calling program. if iz is less than n, the
c                           eigenvectors are not computed. in this case
c                           z is not used.
c                ier    - error parameter
c                         terminal error
c                           ier = 128+j, indicates that ourt2s failed
c                             to converge on eigenvalue j. eigenvalues
c                             and eigenvectors 1,...,j-1 have been
c                             computed correctly, but the eigenvalues
c                             are unordered.
c
c
      dimension          d(*),e(*),z(iz,*)
      data               rdelp /0.745058e-08/
      data               zero,one/0.0,1.0/
c                                  move the last n-1 elements
c                                  of e into the first n-1 locations
c                                  first executable statement
      ier  = 0
      if (n .eq. 1) go to 9005
      do 5  i=2,n
         e(i-1) = e(i)
    5 continue
      e(n) = zero
      b = zero
      f = zero
      do  60  l=1,n
         j = 0
         h = rdelp*(abs(d(l))+abs(e(l)))
         if (b.lt.h) b = h
c                                  look for small sub-diagonal element
         do 10  m=l,n
            k=m
            if (abs(e(k)) .le. b) go to 15
   10    continue
   15    m = k
         if (m.eq.l) go to 55
   20    if (j .eq. 30) go to 85
         j = j+1
         l1 = l+1
         g = d(l)
         p = (d(l1)-g)/(e(l)+e(l))
         r = sqrt(p*p+one)
         d(l) = e(l)/(p+sign(r,p))
         h = g-d(l)
         do 25 i = l1,n
            d(i) = d(i)-h
   25    continue
         f = f+h
c                                  ql transformation
         p = d(m)
         c = one
         s = zero
         mm1 = m-1
         mm1pl = mm1+l
         if (l.gt.mm1) go to 50
         do 45 ii=l,mm1
            i = mm1pl-ii
            G = C*E(I)
            H = C*P
            if (abs(p).lt.abs(e(i))) go to 30
            c = e(i)/p
            r = sqrt(c*c+one)
            e(i+1) = s*p*r
            s = c/r
            c = one/r
            go to 35
   30       c = p/e(i)
            r = sqrt(c*c+one)
            e(i+1) = s*e(i)*r
            s = one/r
            c = c*s
   35       p = c*d(i)-s*g
            d(i+1) = h+s*(c*g+s*d(i))
            if (iz .lt. n) go to 45
c                                  form vector
            do 40 k=1,n
               h = z(k,i+1)
               z(k,i+1) = s*z(k,i)+c*h
               z(k,i) = c*z(k,i)-s*h
   40       continue
   45    continue
   50    e(l) = s*p
         d(l) = c*p
         if ( abs(e(l)) .gt.b) go to 20
   55    d(l) = d(l) + f
   60 continue
c                                  order eigenvalues and eigenvectors
      do  80  i=1,n
         k = i
         p = d(i)
         ip1 = i+1
         if (ip1.gt.n) go to 70
         do 65  j=ip1,n
            if (d(j) .ge. p) go to 65
            k = j
            p = d(j)
   65    continue
   70    if (k.eq.i) go to 80
         d(k) = d(i)
         d(i) = p
         if (iz .lt. n) go to 80
         do 75 j = 1,n
            p = z(j,i)
            z(j,i) = z(j,k)
            z(j,k) = p
   75    continue
   80 continue
      go to 9005
   85 ier = 128+l
 9000 continue
 9005 return
      end

ext: f


program OOPlaneMestTest

use NumberKindsModule
use LoggerModule
use OutputWriterModule
use UtilitiesModule
use PolyMeshOOModule
use FieldOOModule

implicit none

type(PolyMesh2d), pointer :: mesh=>null()
type(Field), pointer :: scalarField => null()
type(Field), pointer :: vectorField => null()
integer(kint) ::  initNest, i, j, doOutput
character(len=56) :: mesh_type, vtkFile

real(kreal), parameter :: intCosBump = 16.0_kreal / PI**2
real(kreal), parameter :: triIntExact = 1.4741613930426135683_kreal
real(kreal), parameter :: dx = 2.0_kreal / 50.0_kreal
real(kreal), dimension(3), parameter :: outsidePoint = [3.0_kreal, 3.0_kreal, dzero]
real(kreal), dimension(3), parameter :: insidePoint = [0.5_kreal, -0.333_kreal, dzero]

type(Logger) :: exeLog
character(len=56) :: logkey = "OOPlaneMestTest"
character(len=MAX_STRING_LENGTH) :: logstring

call New(exeLog, DEBUG_LOGGING_LEVEL)
call getInput(initNest, doOutput)

call LogMessage(exeLog, TRACE_LOGGING_LEVEL,trim(logkey)//" part 1: build each planar mesh to initNest = ", initNest)
call LogMessage(exeLog, TRACE_LOGGING_LEVEL,trim(logkey), " part 2: define scalar and vector fields on mesh.")
call LogMessage(exeLog, TRACE_LOGGING_LEVEL,trim(logkey), " part 3: demonstrate quadrature on mesh.")
call LogMessage(exeLog, TRACE_LOGGING_LEVEL,trim(logkey), " part 4: demonstrate interpolation on a cubic mesh.")
call LogMessage(exeLog, TRACE_LOGGING_LEVEL,trim(logKey), " part 5: output to vtk xml file.")

!-------------------------------------
!   Test Planar mesh, linear triangles
!-------------------------------------
mesh_type = "planar_tri"
write(vtkFile,'(A,I1,A)') trim(mesh_type), initNest, '.vtp'
call StartSection(exeLog,'Planar Tri Test.')
call runTest(mesh, mesh_type, initNest, vtkFile)
call EndSection(exeLog)

!-------------------------------------
!   Test Planar mesh, linear quads
!-------------------------------------
mesh_type = "planar_quad"
write(vtkFile,'(A,I1,A)') trim(mesh_type), initNest, '.vtp'
call StartSection(exeLog,'Planar Quad Test.')
call runTest(mesh, mesh_type, initNest, vtkFile)
call EndSection(exeLog)

!-------------------------------------
!   Test Planar mesh, cubic quads
!-------------------------------------

mesh_type = "planar_cubic_quad"
write(vtkFile,'(A,I1,A)') trim(mesh_type), initNest, '.vtp'
call StartSection(exeLog,'Planar Cubic Quad Test.')
call runTest(mesh, mesh_type, initNest, vtkFile)
call EndSection(exeLog)

if (testPass) call LogMessage(exeLog, TRACE_LOGGING_LEVEL, trim(logkey)//" result:", " PASSED.")

contains

subroutine runTest(mesh_ptr, mesh_type, initNest, oname)
    type(PolyMesh2d), pointer, intent(inout) :: mesh_ptr
    character(len=*), intent(in) :: mesh_type
    integer(kint), intent(in) :: initNest
    character(len=*), intent(in) :: oname
    !
    real(kreal) :: computedIntegral
    real(kreal), dimension(3) :: interpTargetPt
    integer(kint) :: i, j
    real(kreal) :: interpErr, lagErr
    real(kreal), dimension(51) :: unifX, unifY
    real(kreal), dimension(51,51) :: interp, exact, interpX0, interpY0
    real(kreal), dimension(3) :: xyz0
    type(PolyMesh2d) :: copiedMesh

    allocate(mesh_ptr)
    allocate(scalarField)
    allocate(vectorField)
    call mesh_ptr%init(mesh_type, initNest, initNest, 0, 1.0_kreal)
    call mesh_ptr%logStats(exeLog)

    call copiedMesh%init(mesh_type, initNest, initNest, 0, 1.0_kreal)
    call copiedMesh%copy(mesh_ptr)
    call StartSection(exeLog, "mesh copy")
    call copiedMesh%logStats(exeLog)
    call EndSection(exeLog)

    write(logstring,'(2A)') "point (3,3) is outside mesh; mesh%pointIsOutsideMesh((3,3)) = ", &
        merge("true ", "false", copiedMesh%pointIsOutsideMesh(outsidePoint))
    call LogMessage(exeLog, TRACE_LOGGING_LEVEL, "outside mesh test: ", trim(logstring))
    write(logstring,'(2A)') "point (1/2,-1/3) is inside mesh; mesh%pointIsOutsideMesh((1/2,-1/3)) = ", &
        merge("true ", "false", copiedMesh%pointIsOutsideMesh(insidePoint))
    call LogMessage(exeLog, TRACE_LOGGING_LEVEL, "outside mesh test: ", trim(logstring))

    call scalarField%init(1, mesh_ptr%particles%N, "scalar")
    do i=1, mesh_ptr%particles%N
        call scalarField%insertScalar(cosineBump([mesh_ptr%particles%x(i), mesh_ptr%particles%y(i)]))
    enddo

    call vectorField%init(2, mesh_ptr%particles%N, "scalarGradient")
    do i=1, mesh_ptr%particles%n
        call vectorField%insertVector(bumpGrad([mesh_ptr%particles%x(i), mesh_ptr%particles%y(i)]))
    enddo

    computedIntegral = mesh_ptr%integrateScalar(scalarField)
    interpErr = dzero
    if (mesh_ptr%faceKind == QUAD_CUBIC_PANEL) then
        do j=1,51
            unifx(j) = -one + dx*(j-1)
        enddo
        do i=1, 51
            unify(i) = -one + dx*(i-1)
        enddo
        do j=1,51
            do i=1,51
                exact(i,j) = cosineBump([unifX(j), unifY(i)])
            enddo
        enddo
        interpTargetPt = dzero
        lagErr = dzero
        do j=1,51
            do i=1,51
                interpTargetPt(1) = unifx(j)
                interpTargetPt(2) = unify(i)
                interp(i,j) = mesh_ptr%interpolateScalar(interpTargetPt, scalarField)
                xyz0 = mesh_ptr%interpolateLagCoord(interpTargetPt)
                interpX0(i,j) = xyz0(1)
                interpY0(i,j) = xyz0(2)

                !write(logString,'(2(A,2G15.8),A)') "(x0,y0) = (", unifx(j), unify(i), &
                !    "), interp. (x0, y0) = (", xyz0(1), xyz0(2), ")"
                !call LogMessage(exeLog, TRACE_LOGGING_LEVEL, "lagCoordInterp: ", logstring)

                if (sqrt((unifx(j)-xyz0(1))**2 + (unify(i)-xyz0(2))**2) > lagErr) &
                    lagErr = sqrt((unifx(j)-xyz0(1))**2 + (unify(i)-xyz0(2))**2)
            enddo
        enddo
        interpErr = maxval(abs(exact-interp))
        call LogMessage(exeLog, TRACE_LOGGING_LEVEL, "scalar interp error = ", interpErr)
        call LogMessage(exeLog, TRACE_LOGGING_LEVEL, "lagCoord interp error = ", lagErr)
    endif

    call LogMessage(exeLog, TRACE_LOGGING_LEVEL, "scalar integral = ", computedIntegral)
    if (mesh_ptr%faceKind == TRI_PANEL) then
        call LogMessage(exeLog, TRACE_LOGGING_LEVEL, "mesh size = ", mesh_ptr%edges%avgLength(mesh_ptr%particles))
        call LogMessage(exeLog, TRACE_LOGGING_LEVEL, "abs(integralErr) = ", abs(computedIntegral-triIntExact))
    elseif (mesh_ptr%faceKind == QUAD_PANEL .or. mesh_ptr%faceKind==QUAD_CUBIC_PANEL) then
        call LogMessage(exeLog, TRACE_LOGGING_LEVEL, "mesh size = ", mesh_ptr%edges%avgLength(mesh_ptr%particles))
        call LogMessage(exeLog, TRACE_LOGGING_LEVEL, "abs(integralErr) = ", abs(computedIntegral-intCosBump))
    endif

    if (doOutput>0) then
        open(unit=WRITE_UNIT_1, file=vtkFile, action='write', status='replace')
        call mesh_ptr%writeVTKSerialStartXML(WRITE_UNIT_1)
        write(WRITE_UNIT_1,'(A)') '     <PointData Scalars="sineWave" Vectors="gradSineWave">'
        call scalarField%writeVtkPointDataXML(WRITE_UNIT_1)
        call vectorField%writeVtkPointDataXML(WRITE_UNIT_1)
        write(WRITE_UNIT_1,'(A)') '     </PointData>'
        call mesh_ptr%writeVTKSerialEndXML(WRITE_UNIT_1)
        close(WRITE_UNIT_1)
        write(vtkFile,'(A,I1,A)') trim(mesh_type), initNest, '.m'

        if (mesh_ptr%faceKind == QUAD_CUBIC_PANEL) then
            open(unit=WRITE_UNIT_1, file=vtkFile, action='write', status='replace')
            call mesh%writeMatlab(WRITE_UNIT_1)
            call scalarField%writeMatlab(WRITE_UNIT_1)
            call vectorField%writeMatlab(WRITE_UNIT_1)
            call writeToMatlab(unifX, WRITE_UNIT_1, "unifx")
            call writeToMatlab(unify, WRITE_UNIT_1, "unify")
            call writeToMatlab(exact, WRITE_UNIT_1, "exact")
            call writeToMatlab(interp, WRITE_UNIT_1, "interp")
            call writeToMatlab(interpX0, WRITE_UNIT_1, "x0interp")
            call writeToMatlab(interpY0, WRITE_UNIT_1, "y0interp")
            close(WRITE_UNIT_1)
        endif
    endif

    deallocate(vectorField)
    deallocate(scalarField)
    deallocate(mesh)
    nullify(mesh)
    nullify(vectorField)
    nullify(scalarField)
end subroutine

subroutine getInput(initNest, doOutput)
    integer(kint), intent(out) :: initNest, doOutput
    !
    integer(kint) :: argc
    character(len=MAX_STRING_LENGTH) :: argv
    argc = IARGC()
    if (argc == 0) then
        initNest = 0
        doOutput = 1
    elseif (argc==1) then
        call GETARG(1,argv)
        read(argv,*) initNest
        doOutput = 0
    else
        call GETARG(1,argv)
        read(argv,*) initNest
        call GETARG(2,argv)
        read(argv,*) doOutput
    endif
end subroutine

pure function cosineBump(xy)
    real(kreal) :: cosineBump
    real(kreal), intent(in) :: xy(2)
    cosineBump = cos(0.5_kreal * PI * xy(1)) * cos(0.5_kreal * PI * xy(2))
end function

pure function bumpGrad(xy)
    real(kreal) :: bumpGrad(2)
    real(kreal), intent(in) :: xy(2)
    bumpGrad(1) = -0.5_kreal * PI * sin(0.5_kreal*PI*xy(1)) * sin(0.5_kreal*PI*xy(2))
    bumpGrad(2) = -0.5_kreal *PI * cos(0.5_kreal*PI*xy(1))*sin(0.5_kreal*PI*xy(2))
end function

pure function sineWave(xy)
    real(kreal) :: sineWave
    real(kreal), dimension(2), intent(in) :: xy
    sineWave = sin(2.0_kreal*PI* xy(1)*1.0_kreal) * sin(2.0_kreal*PI*xy(2)*1.0_kreal)
end function

pure function sineWaveGrad(xy)
    real(kreal), dimension(2) :: sineWaveGrad
    real(kreal), dimension(2), intent(in) :: xy

    sinewaveGrad(1) = 2.0_kreal*PI*1.0_kreal * cos(2.0_kreal*PI*xy(1)*1.0_kreal) * sin(2.0_kreal*PI*xy(2)*1.0_kreal)
    sinewaveGrad(2) = 2.0_kreal*PI*1.0_kreal * sin(2.0_kreal*PI*xy(1)*1.0_kreal) * cos(2.0_kreal*PI*xy(2)*1.0_kreal)
end function

end program


text:
!
!     (c) 2019 Guide Star Engineering, LLC
!     This Software was developed for the US Nuclear Regulatory Commission (US NRC)
!     under contract "Multi-Dimensional Physics Implementation into Fuel Analysis under
!     Steady-state and Transients (FAST)", contract # NRC-HQ-60-17-C-0007
!
!
!    NEMO - Numerical Engine (for) Multiphysics Operators
! Copyright (c) 2007, Stefano Toninel
!                     Gian Marco Bianchi  University of Bologna
!              David P. Schmidt    University of Massachusetts - Amherst
!              Salvatore Filippone University of Rome Tor Vergata
! All rights reserved.
!
! Redistribution and use in source and binary forms, with or without modification,
! are permitted provided that the following conditions are met:
!
!     1. Redistributions of source code must retain the above copyright notice,
!        this list of conditions and the following disclaimer.
!     2. Redistributions in binary form must reproduce the above copyright notice,
!        this list of conditions and the following disclaimer in the documentation
!        and/or other materials provided with the distribution.
!     3. Neither the name of the NEMO project nor the names of its contributors
!        may be used to endorse or promote products derived from this software
!        without specific prior written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
! ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
! WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
! DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
! ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
! (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
! ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
!
!---------------------------------------------------------------------------------
!
! $Id: tools_part.f90 8157 2014-10-09 13:02:44Z sfilippo $
!
! Description:
!
SUBMODULE(tools_part) tools_part_procedures

    IMPLICIT NONE

    !INTEGER, ALLOCATABLE, SAVE :: part_cells(:)
    !TYPE(connectivity), SAVE :: c2v, c2f

    ! PART_CELLS is created in CMP_MESH_PART
    ! C2V and C2F are created in CMP_MESH_DESC

    ! PART_CELLS, C2V and C2F are destroyed in CMP_MESH_DESC

CONTAINS

    MODULE PROCEDURE part_verts
        USE tools_psblas
        !
        INTEGER :: i, ic, n, pid
        INTEGER :: ppv(nprocs)
        INTEGER, POINTER :: ic2v(:) => NULL()

        ppv(:) = -nverts

        CALL c2v%get_ith_conn(ic2v,iv)
        n = SIZE(ic2v)

        DO i = 1, n
            ic = ic2v(i)
            pid = part_cells(ic)
            ppv(pid+1) = pid
        END DO

        nv       = COUNT(ppv >= 0)
        ppv(1:nv) = PACK(ppv,ppv >= 0)
        IF (nv <= 0) THEN
            WRITE(0,*) 'part_verts: orphan vertex',iv,nverts,n
        END IF
        NULLIFY(ic2v)

        pv(1:nv) = ppv(1:nv)

    END PROCEDURE part_verts

    MODULE PROCEDURE part_faces
        !
        INTEGER :: i, ic, n, pid
        INTEGER, POINTER :: ic2f(:) => NULL()
        INTEGER :: ppv(nprocs)

        ppv(:) = -nfaces

        CALL c2f%get_ith_conn(ic2f,IF)
        n = SIZE(ic2f)

        DO i = 1, n
            ic = ic2f(i)
            pid = part_cells(ic)
            ppv(pid+1) = pid
        END DO

        nv        = COUNT(ppv >= 0)
        ppv(1:nv) = PACK(ppv,ppv >= 0)

        NULLIFY(ic2f)
        pv(1:nv) = ppv(1:nv)

    END PROCEDURE part_faces

    ! REMARK-1: the interface of PART_VERTS and PART_FACES, needed for
    ! allocating the PSBLAS descriptors, is fixed by PARTS.FH, in the
    ! PSBLAS library. Hence, the only way for carrying the C2V and C2F
    ! variables inside the two procedures is to handle them as public
    ! shared variables.

    ! REMARK-2: PV could be initialized equal to any negative value.
    ! The choice of NVERTS and NFACES is only for avoiding the detection
    ! of an unused variable.

END SUBMODULE tools_part_procedures

max_stars_repo_name: sourceryinstitute/MORFEUS-Source

max_issues_repo_name: sourceryinstitute/MORFEUS-Source

max_forks_repo_name: sourceryinstitute/MORFEUS-Source


    Subroutine open_bin(IBI,CSTAT,Cfile,IRECL)
    implicit none
!   Open binary file
    INTEGER IBI,IRECL
    CHARACTER*(*) Cfile,CSTAT
!
!   Open requested file
!
    OPEN (IBI,FORM='UNFORMATTED',FILE=Cfile,RECL=IRECL, &
         ACCESS='DIRECT',STATUS=CSTAT,ERR=902)
!
     RETURN
!*********************************************************
902 WRITE (*,*) Cfile(1:INDEX(Cfile,' '))
    WRITE (*,*) 'Error opening requested file',' Program stopped'
    STOP
    END


text:
C << MAIN PROGRAM>>
      IMPLICIT double precision (A-H,O-Z)
C *** Step1 input data ***
      write(6,*) 'input init value x1 and x2='
      read(5,*) x1, x2
      write(6,*) 'criteria of convergence epsilon='
      read(5,*) e
C *** Step 2 calculate bisection method ***
 100  x=(x1+x2)/2.0
      if(f(x1)*f(x).lt.0) then
          x2=x
      else
          x1=x
      end if

      if(abs(x1-x2).ge.e) then
          goto 100
      end if
C *** Step 3 output result
      write(6,200) 'answer = ',x
 200  format(1h,a6,f10.6)
      stop
      end

c << define nonlinear equation >>
      double precision function f(x)
c      x : value
      implicit double precision (A-H,O-Z)
      f=x**3-1
      return
      end

max_stars_repo_name: terasakisatoshi/Fortran

max_issues_repo_name: terasakisatoshi/Fortran

max_forks_repo_name: terasakisatoshi/Fortran


text:
subroutine matgen(imat, a, n, info)
  implicit none

  real(kind(1.d0)), allocatable :: a(:,:)
  integer :: imat, n, info

  real(kind(1.d0)), allocatable :: work(:), d(:)
  character(len=3) :: path
  character :: type, dist
  integer :: kl, ku, mode
  real(kind(1.d0)) :: anorm, cndnum, rcondc
  integer :: iseed(4)
  real(kind(1.d0)), parameter :: done=1.d0, dmone=-1.d0

  !ehouarn
  integer :: i,j
  double precision :: eps, norme_f
  intrinsic epsilon
  !ehouarn

  iseed = (/1988, 1989, 1990, 1991/)
  info  = 0
  path(1:1) = 'Double precision'
  path(2:3) = 'GE'

  select case(imat)
  case(1:11)
     call dlatb4( path, imat, n, n, type, kl, ku, anorm, mode, cndnum, dist )
     rcondc = 1.d0 / cndnum

     write(*,'("=======================")')
     ! write(*,'("Type  : ",a1)')type
     ! write(*,'("KL    : ",i4)')kl
     ! write(*,'("KL    : ",i4)')ku
     write(*,'("Anorm : ",es10.3)')anorm
     write(*,'("Cond  : ",es10.3)')cndnum
  write(*,'(" ")')

     allocate(a(n,n), work(3*n), d(n), stat=info)
     if(info .ne. 0) return

     call dlatms( n, n, dist, iseed, type, &
          & d, mode, cndnum, anorm, kl, ku, &
          & 'no packing', a, n, work, info )
     if(info .ne. 0) return

     deallocate(d, work, stat=info)
  case(12)
      !write(*,'("Generate matrix with bad determinant")')
      !info = 1  ! to remove when the generator is implemented
      write(*,'("Generate numerically singular matrix ")')
      call dlatb4( path, 8, n, n, type, kl, ku, anorm, mode, cndnum, dist )
      rcondc = 1.d0 / cndnum

      write(*,'("=======================")')
      ! write(*,'("Type  : ",a1)')type
      ! write(*,'("KL    : ",i4)')kl
      ! write(*,'("KL    : ",i4)')ku
      write(*,'("Anorm : ",es10.3)')anorm
      write(*,'("Cond  : ",es10.3)')cndnum
      write(*,'(" ")')

      allocate(a(n,n), work(3*n), d(n), stat=info)
      if(info .ne. 0) return

      call dlatms( n, n, dist, iseed, type, &
          & d, mode, cndnum, anorm, kl, ku, &
          & 'no packing', a, n, work, info )
      if(info .ne. 0) return

      deallocate(d, work, stat=info)

      eps=epsilon(eps)
      a(:,1)=eps*a(:,1)
      if(n.gt.5)then
        i=(n+1)/2
        a(:,i)=eps*a(:,i)
      endif
      a(:,n)=eps*a(:,n)

  case default
     write(*,'("Wrong matrix id")')
     info = 1
  end select

  return
end subroutine matgen

max_stars_repo_name: MOUDDENEHamza/ENSEEIHT

max_issues_repo_name: Hathoute/ENSEEIHT

max_forks_repo_name: Hathoute/ENSEEIHT


text:
      subroutine inhyco
!***********************************************************************
!  Copyright, 2004,  The  Regents  of the  University of California.
!  This program was prepared by the Regents of the University of
!  California at Los Alamos National Laboratory (the University) under
!  contract No. W-7405-ENG-36 with the U.S. Department of Energy (DOE).
!  All rights in the program are reserved by the DOE and the University.
!  Permission is granted to the public to copy and use this software
!  without charge, provided that this Notice and any statement of
!  authorship are reproduced on all copies. Neither the U.S. Government
!  nor the University makes any warranty, express or implied, or
!  assumes any liability or responsibility for the use of this software.
!***********************************************************************
!D1
!D1 PURPOSE
!D1
!D1 Read hydraulic conductivity data.
!D1
!***********************************************************************
!D2
!D2 REVISION HISTORY
!D2
!D2 FEHM Version 2.20
!D2
!D2 Initial implementation: 09-FEB-01, Programmer: George Zyvoloski
!D2
!D2 $Log:   /pvcs.config/fehm90/src/inhyco.f_a  $
!D2
!D2    Rev 2.5   06 Jan 2004 10:43:18   pvcs
!D2 FEHM Version 2.21, STN 10086-2.21-00, Qualified October 2003
!D2
!D2    Rev 2.4   29 Jan 2003 09:08:54   pvcs
!D2 FEHM Version 2.20, STN 10086-2.20-00
!D2
!***********************************************************************
!D3
!D3 REQUIREMENTS TRACEABILITY
!D3
!D3 2.6 Provide Input/Output Data Files
!D3 3.0 INPUT AND OUTPUT REQUIREMENTS
!D3
!***********************************************************************
!D4
!D4 SPECIAL COMMENTS AND REFERENCES
!D4
!D4 Requirements from SDN: 10086-RD-2.20-00
!D4   SOFTWARE REQUIREMENTS DOCUMENT (RD) for the
!D4   FEHM Application Version 2.20
!D4
!***********************************************************************

      use comai
      use comdi
      use comdti
      use comki
      implicit none

      integer i
      real*8 conv
      real*8, allocatable ::  hycox(:)
      real*8, allocatable ::  hycoy(:)
      real*8, allocatable ::  hycoz(:)
      parameter (conv=9.89d-13/9.66d-6)
      allocate (hycox(n0))
      allocate (hycoy(n0))
      allocate (hycoz(n0))
      do i = 1,n0
         hycox(i) = pnx(i)/conv
	 hycoy(i) = pny(i)/conv
	 hycoz(i) = pnz(i)/conv
      enddo
      macro = 'hyco'
c**** read hydraulic conductivity data (in m/sec)
      narrays = 3
      itype(1) = 8
      itype(2) = 8
      itype(3) = 8
      default(1) = zero_t
      default(2) = zero_t
      default(3) = zero_t
      igroup = 1
c     convert to intrinsic perms at 20 C and 0.1 Mpa

      call initdata2 (inpt, ischk, n0, narrays, itype,
     1     default, macroread(14), macro, igroup, ireturn,
     2     r8_1=hycox(1:n0),r8_2=hycoy(1:n0),r8_3=hycoz(1:n0))

      do i=1,n0
         if(hycox(i).ge.0.0d00) then
            hycox(i) = max (zero_t, hycox(i))
         else
            hycox(i) = 1.0d00/10.0d00**(abs(hycox(i)))
         endif
         pnx(i) = conv*hycox(i)
         if(hycoy(i).ge.0.0d00) then
            hycoy(i) = max (zero_t, hycoy(i))
         else
            hycoy(i) = 1.0d00/10.0d00**(abs(hycoy(i)))
         endif
         pny(i) = conv*hycoy(i)
         if(hycoz(i).ge.0.0d00) then
            hycoz(i) = max (zero_t, hycoz(i))
         else
            hycoz(i) = 1.0d00/10.0d00**(abs(hycoz(i)))
         endif
         pnz(i) = conv*hycoz(i)
      end do

      deallocate (hycox)
      deallocate (hycoy)
      deallocate (hycoz)

      macroread(14) = .TRUE.

      end

max_stars_repo_name: satkarra/FEHM

max_issues_repo_name: satkarra/FEHM

max_forks_repo_name: satkarra/FEHM


text:
*..............................................................

!-------------------------------------------------------------------------
!         NASA/GSFC, Data Assimilation Office, Code 910.3, GEOS/DAS      !
!-------------------------------------------------------------------------
!BOP
!
! !ROUTINE: ODS_PutR() --- Writes data in native real format to file
!
! !DESCRIPTION:
!    \label{ODS:PutR}
!     Stores the contents of a {\tt real} variable for a full
!     synoptic hour on an opened ODS file.\\
!
!  \noindent {\bf Notes: }\\
!  \begin{enumerate}
!  \item For a list of error codes, see Table~\ref{tab:errors}.
!  \item This routine does perform checks to verify whether
!     the range of values on input is consistent with the internal
!     variable type in the ODS file.  For example, a real variable
!     here could be stored as a two byte variable on file.
!     Therefore, all elements in the array to be stored as two byte
!     integers should have values between -32,767 and 32,767.
!     If the software does detect a number inconsistent with the
!     variable type, then the routine returns with an error
!     message without saving at least some of the values.  This
!     check prevents overflows from occurring which would produce
!     unexpected results and probably abnormal program termination.
!
!  \item The input values should be within the range specified in
!     Table~\ref{tab:ObAtt} or by the ODS file annotated header of
!     the ODS file produced with the MFHDF utility, ncdump.  If no
!     range is specified, then the maximum and minimum values for
!     the internal variable types are as follows:
!
!  \begin{verbatim}
!       type                      minimum value    maximum value
!       ----                      -------------    -------------
!       byte                                  0              255
!       short (2 byte) integer          -32,767           32,767
!       long  (4 byte) integer   -2,147,483,643    2,147,483,643
!       real  (4 bytes )               -3.40e38          3.40e38
!  \end{verbatim}
!
!  \item Once data from the specified date and hour have been
!     written, the routine will allow overwriting of data
!     from a previous date and hour.  The only restriction
!     is that the number to be written does not exceed the
!     space already allocated in the ODS file.  Only the first
!     {\tt nval} values will be overwritten.
!
!  \end{enumerate}
!
! !INTERFACE:
!
      subroutine ODS_PutR ( id, VarName, julian_day, syn_hour,
     .                          nval,    values,     ierr )
!
! !INPUT PARAMETERS:
      implicit        NONE
      integer         id               ! ODS file handle as returned
                                       !  from ODS_Create() or
                                       !  ODS_Open().
      character * (*) VarName          ! name of variable
      integer         julian_day       ! Julian Day.  Use the function
                                       !  ODS_Julian to obtain this
                                       !  number.
      integer         syn_hour         ! hour of synoptic time since
                                       !  0:00 GMT (e.g., 0, 6, 12, 18)
      integer         nval             ! number of values for this
                                       !  synoptic time
      real            values ( nval )  ! values to be written
!
! !OUTPUT PARAMETERS:
      integer         ierr             ! Error code. If non-zero, an
                                       !  error has occurred. For a list
                                       !  of possible values, see the
                                       !  description section of this
                                       !  prologue.
!
! !SEE ALSO:
!     ODS_PutC()    write a list of character  strings to   an ODS file
!     ODS_PutI()    write a list of integer    values  to   an ODS file
!     ODS_GetC()    get   a list of characters strings from an ODS file
!     ODS_GetI()    get   a list of integer    values  from an ODS file
!     ODS_GetR()    get   a list of real       values  from an ODS file
!     ODS_Julian()  convert the "calendar" date to the Julian day
!     ods_hdf.h     include file defining internal constants
!                    and global variables, and setting up data
!                    structures
!
! !REVISION HISTORY:
!     02Oct1995   da Silva    Specification.
!     17May1996   Redder      Original code
!     13Apr1998   C. Redder   Updated documentation to reflect changes
!                             made in other routines to create ODS
!                             version 2.00
!     26May1999   C. Redder   Fixed bug in handling the case when nval
!                             is zero or less.
!     19Nov1999   C. Redder   Added a latex label in and moved the
!                             subroutine statement into the prologue.
!                             Modified the comments for the return status
!                             code.  Removed references to the old version
!                             of the ODS documentation.
!     06Dec1999   C. Redder   Corrections to the documentation in the
!                             prologue.
!     30Nov2001   C. Redder   Fixed bug that prevents last batch from
!                             being processed twice when nval = 0
!     16Sep2002   C. Redder   Updated prologue and added error handling
!                             for call to NetCDF routine, NCDINQ
!     17Oct2002   C. Redder   Fixed bug that can causes an array index
!                             error when only part of the data for the
!                             latest date and time on file is being
!                             overwritten.
!
!EOP
!-------------------------------------------------------------------------

      include 'netcdf.inc'
      include 'ods_hdf.h'
      include 'ods_stdio.h'

*     Functions referenced
*     --------------------
      character * ( 10 )
     .         ODS_VarType ! recommended FORTRAN variable
                           !   type for data
      integer  ODS_JulHr   ! Julian hour using, as the reference,
                           !   the synoptic hour and julian day
                           !   that corresponds to the last block
                           !   of data written to file

*     Arguments for NetCDF library routines
*     -------------------------------------
      integer  nc_id       ! file id
      integer  varid       ! variable id
      integer  start  ( 3 )! starting file index for hyperslab
      integer  count  ( 3 )! length of edges in hyperslab
      character * ( MAXNCNAM )
     .         Var_Name    ! variable name ( returned from NetCDF
                           !   routine NCVINQ )
      integer  Var_Type    ! variable type
      integer  NVDims      ! number of dimensions
      integer  DimIDs ( 3 )! dimension ids
      integer  DimID       ! dimension id
      integer  NVAtts      ! number of attributes
      character * ( MAXNCNAM )
     .         DimName     ! dimension name ( returned from NetCDF
                           !   routine NCDINQ )

*     Indices defining the set of observation data
*     corresponding to a synoptic time in the ODS file
*     ------------------------------------------------
      integer  synbeg      ! starting index for the data block
      integer  synlen      ! number of observation reports in
                           !   the data block

*     Parameters for defining the observation
*     data batches to be written to the ODS file
*     ------------------------------------------
      integer  batch_len   ! number of values in each batch
      integer  first_batch ! index for the first batch
      integer  last_batch  ! index for the last batch
      integer  first_ival  ! data value index of the first
                           !   value in the first batch
      integer  first_nval  ! number of values in the first batch
      integer  last_nval   ! number of values in the last batch
      integer  ibatch      ! batch number index for do loop

*     Other variables
*     ---------------
      integer  ival        ! index value for the input array values
      character * ( 10 )
     .         VarType     ! recommended FORTRAN variable type for
                           !   data
      integer JulHr        ! Julian hour corresponding to julian_day
                           !   and syn_hour
      integer JulHr_Latest ! Julian hour corresponding to the synoptic
                           !   date and hour of the last block of data
                           !   that was written to file
      parameter
     .      ( JulHr_Latest = 0 )
                           ! Set to zero for clarity and convenience

*     Set ierr code to valid
*     ----------------------
      ierr    = NCNoErr

*     Check to determine if the file handle id is valid
*     -------------------------------------------------
      if ( id            .lt. 1      .or.
     .     id            .gt. id_max .or.
     .     IOMode ( id ) .eq. CLOSED ) then
         write ( stderr, 901 )
         ierr = NCEBadID
         return
      end if

*     If IO mode is not set to write,
*     then exit with an error message
*     -------------------------------
      if ( IOMode ( id ) .ne. NCWRITE ) then
         write ( stderr, 902 ) filenames ( id )
         ierr = NCEPerm
         return

      end if

*     Check to determine if variable is appropriate for this routine
*     --------------------------------------------------------------
      VarType    = ODS_VarType ( VarName )
      if ( VarType ( :4  ) .ne. 'real'     .and.
     .     VarType ( :10 ) .ne. 'indefinite' ) then
         write ( stderr, 903 ) VarName
         ierr = ODS_BadInt
         return

      end if

*     Return with an error message if the number of values
*     to be written is a negative integer
*     ----------------------------------------------------
      if ( nval .lt. 0 ) then
         ierr = ODS_DimErr
         write ( stderr, 904 ) VarName
         return

      end if

*     Calculate the Julian hour ( using, as the reference,
*     the date and hour of the last block of data that
*     was written to file )
*     ---------------------------------------------------
      JulHr = ODS_JulHr ( id, julian_day, syn_hour, ierr )

*     For synoptic times earlier than the synoptic time
*     for the last data block written to file ...
*     ------------------------------------------------------
      if      ( JulHr .lt. JulHr_Latest ) then

*        Get pointer data from tables stored in common
*        ---------------------------------------------
         call ODS_GetP    ( id,        'write',
     .                      julian_day, syn_hour,
     .                      synbeg,     synlen, ierr )
         if ( ierr .ne. NCNoErr ) return

*        Return with an error message if the number of values
*        exceeds the space already allocated in the ODS file
*        ----------------------------------------------------
         if ( nval .gt. synlen ) then
            write ( stderr, 905 )
            ierr  = NCEInVal
            return

         end if

         synlen  = nval ! Set syn_len to the user specified value

*     For synoptic times at the same time for the last data
*     block written to file ...
*     -----------------------------------------------------
      else if ( JulHr .eq. JulHr_Latest ) then

*        Get pointer data from tables stored in common
*        ---------------------------------------------
         call ODS_GetP    ( id,        'write',
     .                      julian_day, syn_hour,
     .                      synbeg,     synlen, ierr )
         if ( ierr .ne. NCNoErr ) return

*        Return with an error message if the number of values
*        exceeds the space already allocated in the ODS file
*        ----------------------------------------------------
         if ( nval .gt. synlen ) then
            write ( stderr, 906 )
            ierr  = NCEInVal
            return

         end if

         synlen  = nval ! Set syn_len to the user specified value

*     For synoptic times later than the time for the
*     last data block written to file ...
*     -------------------------------------------------
      else if ( JulHr .gt. JulHr_Latest ) then

*        Update pointer tables ( tables stored in common )
*        -------------------------------------------------
         call ODS_UpdateP ( id,        'no_option',
     .                      julian_day, syn_hour,
     .                      nval,       ierr )
         if ( ierr .ne. NCNoErr ) return

*        Get pointer data from tables stored in common
*        ---------------------------------------------
         call ODS_GetP    ( id,        'write',
     .                      julian_day, syn_hour,
     .                      synbeg,     synlen, ierr )
         if ( ierr .ne. NCNoErr ) return

      end if

*     Nothing more to do if no values are to be written
*     -------------------------------------------------
      if ( nval .le. 0 ) return

*     Get NetCDF file and variable ids
*     --------------------------------
      nc_id = ncid ( id )
      varid = NCVID ( nc_id, VarName, ierr )
      if ( ierr .ne. NCNoErr ) return

*     Get the length of each batch of observation data
*     in an ODS file
*     ------------------------------------------------
      call NCVINQ ( nc_id,    varid,
     .              Var_Name, Var_Type,
     .              NVDims,   DimIDs,
     .              NVatts,   ierr )
      if ( ierr .ne. NCNoErr ) return
      DimID       = DimIDs ( 1 )
      call NCDINQ ( nc_id,    DimID,
     .              DimName,  batch_len, ierr )
      if ( ierr .ne. NCNoErr ) return

*     The remaining code uses the ODS IO routines to write
*     the data in batches.  Each batch contains a subset
*     of the values to be written.  The maximum size of
*     each subset is batchlen as defined in the header
*     file, ods_hdf.h.  The values are being processed in
*     this manner in order to prevent the execution of
*     the NetCDF library routines from being
*     prohibitively slow.
*     ---------------------------------------------------

*     Determine the indices of the first and last batches
*     of data to be written to the ODS file
*     ---------------------------------------------------
      first_batch = ( synbeg - 1 ) / batch_len + 1
      last_batch  = ( synbeg + synlen - 2 ) / batch_len + 1

*     Write the first batch of values
*     -------------------------------
      first_ival  =   mod ( synbeg - 1, batch_len ) + 1
      first_nval  =   min ( synlen, batch_len - first_ival + 1 )
      start ( 1 ) =   first_ival
      start ( 2 ) =   first_batch
      count ( 1 ) =   first_nval
      count ( 2 ) =   1
      call ODS_NCVPTR ( nc_id, varid, start, count, values, ierr )
      if ( ierr .ne. NCNoErr ) return

*     Write data not in the first and last batch of values,
*     if necessary
*     -----------------------------------------------------
      ival = first_nval + 1
      do 10, ibatch = first_batch + 1, last_batch - 1
         start ( 1 ) = 1
         start ( 2 ) = ibatch
         count ( 1 ) = batch_len
         count ( 2 ) = 1
         call ODS_NCVPTR ( nc_id, varid, start, count,
     .                     values ( ival ), ierr )
         if ( ierr .ne. NCNoErr ) return
         ival        = ival + batch_len
 10   continue

*     Write the last batch of values, if necessary
*     --------------------------------------------
      start ( 1 ) =   1
      start ( 2 ) =   last_batch
      count ( 1 ) =   mod ( synbeg + synlen - 2, batch_len ) + 1
      count ( 2 ) =   1
      if ( last_batch .gt. first_batch )
     .   call ODS_NCVPTR ( nc_id, varid, start, count,
     .                     values ( ival ), ierr )
      if ( ierr .ne. NCNoErr ) return

      return
*     ------

 901  format ( /, ' ODS_PutR: File handle id number does not ',
     .         /, '           correspond to an opened ODS file' )
 902  format ( /, ' ODS_PutR: Cannot write to file unless mode ',
     .         /, '           is set to write.  Filename is: ', a )
 903  format ( /, ' ODS_PutR: Inappropriate variable for this ',
     .         /, '           routine.  Variable name is: ', a )
 904  format ( /, ' ODS_PutR: The number of values to be written ',
     .         /, '           must be a non-negative integer. ',
     .         /, '           Variable name is: ', a )
 905  format ( /, ' ODS_PutR: Cannot increase the size of a ',
     .         /, '           block of data that corresponds ',
     .         /, '           to a date and time earlier than ',
     .         /, '           the latest date and time for ',
     .         /, '           which there is data. ' )
 906  format ( /, ' ODS_PutR: Must use ODS_Append to increase the '
     .         /, '           size of a block of data that corre-'
     .         /, '           sponds to the latest date and time '
     .         /, '           for which there is data. ' )

      end

max_stars_repo_name: GEOS-ESM/AeroApps

max_forks_repo_name: GEOS-ESM/AeroApps


text:
module IO
! Description:
!    Contains I/O procedures and file units.
!
! Author:
!    Davide Amato
!    Space Dynamics Group - Technical University of Madrid
!    d.amato@upm.es
!
! ==============================================================================

! MODULE
use KINDS, only: dk,qk

! VARIABLES
implicit none
! Input file unit
integer,parameter  ::  iin = 10, id_set = 11
integer,parameter  ::  id_ares = 12, id_rres = 13, id_enres = 14
integer,parameter  ::  id_stats = 15
integer,parameter  ::  id_sma = 16

contains

subroutine READ_SETTINGS()
! Description:
!    Reads the file "input.txt", containing the settings for the ensemble of
!    propagations. The file must be in the same directory as the executable.
!
! ==============================================================================

! MODULES
use SETTINGS

! VARIABLES
implicit none
character(len=13)   ::  title

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

write(*,'(/,a)') 'Reading input file...'

open(unit=iin,file='input.txt',status='old',action='read')

write(*,'(a,/)') 'Done.'

write(*,'(a)') 'SIMULATION SETTINGS:'

call SKIP(10)

! Read parameter space intervals
read(iin,'(a13,2(e23.15))') title, d_int
read(iin,'(a13,2(e23.15))') title, e_int
read(iin,'(a13,2(e23.15))') title, th_int
read(iin,'(a13,2(e23.15))') title, RSw_int
write(*,'(a16,g14.7,'', '',g14.7,a)') 'd interval: ', d_int, ' Earth radii'
write(*,'(a16,g14.7,'', '',g14.7)')   'e interval: ', e_int
write(*,'(a16,g14.7,'', '',g14.7,a)') 'theta interval: ', th_int, ' deg'
write(*,'(a16,g14.7,'', '',g14.7,a)') 'R_sw interval: ', RSW_int, ' AU'

call SKIP(7)

! Read mesh settings
read(iin,'(a13,i4)') title, n_d
read(iin,'(a13,i4)') title, n_e
read(iin,'(a13,i4)') title, n_th
read(iin,'(a13,i4)') title, n_rs
read(iin,'(a13,i4)') title, espacing

write(*,'(a16,i4)') 'd points: ',n_d
write(*,'(a16,i4)') 'e points: ', n_e
write(*,'(a16,i4)') 'theta points: ', n_th
write(*,'(a16,i4)') 'R_Sw points: ', n_rs
write(*,'(a16,i10,/)') 'Total: ',n_d*n_e*n_th*maxval([n_rs,1])
if (n_rs == 0) then
  write(*,'(a16,a)') 'Mode: ','HELIO'
else if (n_rs > 0) then
  write(*,'(a16,a)') 'Prop mode: ','OTM'
end if
select case(espacing)
case(1)
    write(*,'(a16,a)') 'e spacing: ','  Linear'
case(2)
    write(*,'(a16,a)') 'e spacing: ','  Log'
end select

call SKIP(8)

! Read integrator settings
read(iin,'(a13,i2)') title, integ
read(iin,'(a13,e22.15)') title, tolref
read(iin,'(a13,e22.15)') title, tol
read(iin,'(a13,e22.15)') title, dt_H
read(iin,'(a13,e22.15)') title, dt_CE
read(iin,'(a13,i9)')     title, mxsteps

select case (integ)
case(1)
  write(*,'(a16,a)') 'Integrator: ','  LSODAR'
case(2)
  write(*,'(a16,a)') 'Integrator ','  Radau'
end select
write(*,'(a16,g14.7)') 'Ref. tol: ', tolref
write(*,'(a16,g14.7)') 'Test tol: ', tol
write(*,'(a16,g14.7,a)') 'dt_H: ',dt_H,' days'
write(*,'(a16,g14.7,a)') 'dt_CE: ',dt_CE,' days'
write(*,'(a16,i9)')    'Max steps: ',mxsteps

call SKIP(5)

! Read choice of equations of motion (heliocentric)
read(iin,'(a13,i2)') title, eqs_H
select case(eqs_H)
case default
   write(*,'(a,i2)') 'eqs_H = ',eqs_H
   write(*,'(a)') 'Invalid choice of the equations of motion for the heliocentric phase.'
   write(*,'(a)') 'Stopping...'
   stop
case(-1)
   write(*,'(a16,a)') 'Equations (H): ', '  Cowell (II)'
case(1)
   write(*,'(a16,a)') 'Equations (H): ', '  Cowell (I)'
case(2)
   write(*,'(a16,a)') 'Equations (H): ', '  K-S'
case(3)
   write(*,'(a16,a)') 'Equations (H): ', '  EDromo'
end select
! Read choice of equations of motion (close encounter)
read(iin,'(a13,i2)') title, eqs_CE
select case(eqs_CE)
case default
   write(*,'(a,i2)') 'eqs_CE = ',eqs_CE
   write(*,'(a)') 'Invalid choice of the equations of motion for the close encounter phase.'
   write(*,'(a)') 'Stopping...'
   stop
case(-1)
   write(*,'(a16,a)') 'Equations (CE): ', '  Cowell (II)'
case(1)
   write(*,'(a16,a)') 'Equations (CE): ', '  Cowell (I)'
case(2)
   write(*,'(a16,a)') 'Equations (CE): ', '  K-S'
case(4)
   write(*,'(a16,a)') 'Equations (CE): ', '  GDromo'
end select

call SKIP(6)

! Read regularized formulation settings
read(iin,'(a13,i2)') title, flag_time_EDr
select case(flag_time_EDr)
case(0)
    write(*,'(a16,a)') 'EDromo TE: ', '  None'
case(1)
    write(*,'(a16,a)') 'EDromo TE: ', '  Constant'
case(2)
    write(*,'(a16,a)') 'EDromo TE: ', '  Linear'
end select
read(iin,'(a13,i2)') title, flag_time_GDr
select case(flag_time_GDr)
case(0)
    write(*,'(a16,a)') 'GDromo TE: ', '  None'
case(1)
    write(*,'(a16,a)') 'GDromo TE: ', '  Constant'
case(2)
    write(*,'(a16,a)') 'GDromo TE: ', '  Linear'
end select

call SKIP(2)

! Read output path
read(iin,'(a13,a)') title, path
write(*,'(a16,a)') 'Out path: ', '  '//trim(path)

close(iin)

contains

subroutine SKIP(n)
integer,intent(in)  ::  n
integer  ::  i
character(len=120)  ::  dummy
do i=1,n
    read(iin,'(a)') dummy
end do
end subroutine SKIP

end subroutine READ_SETTINGS

subroutine WRITE_SETTINGS(pathname,runID)
! Description:
!    Saves settings of the ensemble propagation to disk in the output directory.
!
! ==============================================================================

use SETTINGS
implicit none
character(len=*),intent(in)  ::  pathname,runID

open(id_set,file=pathname,status='replace',action='write')

write(id_set,'(120(''#''))')
write(id_set,'(2(''#''),55x,a7,54x,2(''#''))') 'NAPLES'
write(id_set,'(2(''#''),49x,a19,48x,2(''#''))') 'SIMULATION SETTINGS'
write(id_set,'(120(''#''),/)')

write(id_set,'(a16,a)') 'Run ID: ', '  '//runID
write(id_set,'(a16,g14.7,'', '',g14.7,a)') 'd interval: ', d_int, ' Earth radii'
write(id_set,'(a16,g14.7,'', '',g14.7)')   'e interval: ', e_int
write(id_set,'(a16,g14.7,'', '',g14.7,a)') 'theta interval: ', th_int, ' deg'
write(id_set,'(a16,g14.7,'', '',g14.7,a)') 'RSw interval: ', RSw_int, ' AU'
write(id_set,'(a16,i4)') 'd points: ',n_d
write(id_set,'(a16,i4)') 'e points: ', n_e
write(id_set,'(a16,i4)') 'theta points: ', n_th
write(id_set,'(a16,i4)') 'RSw points: ', n_rs
select case(espacing)
case(1)
    write(id_set,'(a16,a)') 'e spacing: ','  Linear'
case(2)
    write(id_set,'(a16,a)') 'e spacing: ','  Log'
end select
select case (integ)
case(1)
  write(id_set,'(a16,a)') 'Integrator: ','  LSODAR'
case(2)
  write(id_set,'(a16,a)') 'Integrator ','  Radau'
end select
if (n_rs == 0) then
  write(id_set,'(a16,a)') 'Prop. mode: ','HELIO'
else if (n_rs > 0) then
  write(id_set,'(a16,a)') 'Prop. mode: ','OTM'
end if
write(id_set,'(a16,g14.7)') 'Ref. tol: ', tolref
write(id_set,'(a16,g14.7)') 'Test tol: ', tol
write(id_set,'(a16,g14.7,a)') 'dt_H: ', dt_H,' days'
write(id_set,'(a16,g14.7,a)') 'dt_CE: ', dt_CE,' days'
write(id_set,'(a16,i9)')    'Max steps: ',mxsteps
select case(eqs_H)
case(1,-1)
   write(id_set,'(a16,a)') 'Equations (H): ', '  Cowell'
case(2)
   write(id_set,'(a16,a)') 'Equations (H): ', '  K-S'
case(3)
   write(id_set,'(a16,a)') 'Equations (H): ', '  EDromo'
end select
select case(eqs_CE)
case(1,-1)
   write(id_set,'(a16,a)') 'Equations (CE): ', '  Cowell'
case(2)
   write(id_set,'(a16,a)') 'Equations (CE): ', '  K-S'
case(4)
   write(id_set,'(a16,a)') 'Equations (CE): ', '  GDromo'
end select
select case(flag_time_EDr)
case(0)
    write(id_set,'(a16,a)') 'EDromo TE: ', '  None'
case(1)
    write(id_set,'(a16,a)') 'EDromo TE: ', '  Constant'
case(2)
    write(id_set,'(a16,a)') 'EDromo TE: ', '  Linear'
end select
select case(flag_time_GDr)
case(0)
    write(id_set,'(a16,a)') 'GDromo TE: ', '  None'
case(1)
    write(id_set,'(a16,a)') 'GDromo TE: ', '  Constant'
case(2)
    write(id_set,'(a16,a)') 'GDromo TE: ', '  Linear'
end select
write(id_set,'(a16,a)') 'Out path: ', '  '//trim(path)

close(id_set)

end subroutine WRITE_SETTINGS

subroutine CREATE_ENRES(pathname)
! Description:
!    Creates the file containing the residuals in energy.
!
! ==============================================================================

implicit none
! Arguments IN
character(len=*),intent(in)    ::  pathname

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

open(unit=id_enres,file=pathname,status='replace',action='write')
write(id_enres,'(120(''#''))')
write(id_enres,'(2(''#''),55x,a7,54x,2(''#''))') 'NAPLES'
write(id_enres,'(2(''#''),50x,a16,50x,2(''#''))') 'ENERGY RESIDUALS'
write(id_enres,'(120(''#''),/,''#'')')

write(id_enres,'(''#'',a2,1x,3(a3,1x),7(a14,1x),a2,1x,a3)') &
&'i','j','k','l','dmin [R_E]','ecc [-]','theta [deg]', 'Rsw [AU]',&
&'dEN_IN [-]','dEn_OUT [-]','dEn_END [-]','MF','XC'
write(id_enres,'(127(''#''))')

end subroutine CREATE_ENRES

subroutine CREATE_ABSRES(pathname)
! Description:
!    Creates the file containing the absolute residuals in position.
!
! ==============================================================================

implicit none
! Arguments IN
character(len=*),intent(in)    ::  pathname

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

open(unit=id_ares,file=pathname,status='replace',action='write')
write(id_ares,'(120(''#''))')
write(id_ares,'(2(''#''),55x,a7,54x,2(''#''))') 'NAPLES'
write(id_ares,'(2(''#''),38x,a40,38x,2(''#''))') 'ABSOLUTE POSITION AND VELOCITY RESIDUALS'
write(id_ares,'(120(''#''),/,''#'')')

write(id_ares,'(''#'',a2,1x,3(a3,1x),10(a14,1x),a2,1x,a3)') &
&'i','j','k','l','dmin [R_E]','ecc [-]','theta [deg]','RSw [AU]',&
&'dR_IN [AU]','dR_OUT [AU]','dR_END [AU]','dV_IN [km/s]', 'dV_OUT [km/s]',&
&'dV_END [km/s]','MF','XC'
write(id_ares,'(172(''#''))')

end subroutine CREATE_ABSRES

subroutine CREATE_RELRES(pathname)
! Description:
!    Creates the file containing the relative residuals in position.
!
! ==============================================================================

implicit none
! Arguments IN
character(len=*),intent(in)    ::  pathname

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

open(unit=id_rres,file=pathname,status='replace',action='write')
write(id_rres,'(120(''#''))')
write(id_rres,'(2(''#''),55x,a7,54x,2(''#''))') 'NAPLES'
write(id_rres,'(2(''#''),38x,a40,38x,2(''#''))') 'RELATIVE POSITION AND VELOCITY RESIDUALS'
write(id_rres,'(120(''#''),/,''#'')')

write(id_rres,'(''#'',a2,1x,3(a3,1x),10(a14,1x),a2,1x,a3)') &
&'i','j','k','l','dmin [R_E]','ecc [-]','theta [deg]','RSw [AU]',&
&'dR_IN [-]','dR_OUT [-]','dR_END [-]','dV_IN [-]','dV_OUT [-]','dV_END [-]',&
&'MF','XC'
write(id_rres,'(172(''#''))')

end subroutine CREATE_RELRES

subroutine CREATE_SMARES(pathname)
! Description:
!    Creates the file containing the absolute residuals in semi-major axis.
!
! ==============================================================================

implicit none
! Arguments IN
character(len=*),intent(in)    ::  pathname

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

open(unit=id_sma,file=pathname,status='replace',action='write')
write(id_sma,'(120(''#''))')
write(id_sma,'(2(''#''),55x,a7,54x,2(''#''))') 'NAPLES'
write(id_sma,'(2(''#''),38x,a40,38x,2(''#''))') 'SEMI-MAJOR AXIS ABSOLUTE RESIDUALS'
write(id_sma,'(120(''#''),/,''#'')')

write(id_sma,'(''#'',a2,1x,3(a3,1x),7(a14,1x),a2,1x,a3)') &
&'i','j','k','l','dmin [R_E]','ecc [-]','theta [deg]','RSw [AU]',&
&'dSMA_IN [AU]','dSMA_OUT [AU]','dSMA_END [AU]','MF','XC'
write(id_sma,'(128(''#''))')

end subroutine CREATE_SMARES

subroutine CREATE_STATS(pathname)
! Description:
!    Creates the file containing the integration statistics.
!
! ==============================================================================

implicit none
! Arguments IN
character(len=*),intent(in)    ::  pathname

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

open(unit=id_stats,file=pathname,status='replace',action='write')
write(id_stats,'(120(''#''))')
write(id_stats,'(2(''#''),55x,a7,54x,2(''#''))') 'NAPLES'
write(id_stats,'(2(''#''),47x,a22,47x,2(''#''))') 'INTEGRATION STATISTICS'
write(id_stats,'(120(''#''),/,''#'')')

write(id_stats,'(''#'',a2,1x,3(a3,1x),4(a14,1x),3(a9,1x),3(a14,1x),3(a10,1x),a14,1x,a3,1x,a3)')&
&'i','j','k','l','dmin [R_E]','ecc [-]','theta [deg]','RSw [AU]',&
&'FEVAL_H1','FEVAL_CE','FEVAL_END','MORD_H1','MORD_CE','MORD_END',&
&'ISTEPS_H1','ISTEPS_CE','ISTEPS_END','RSw_diff','MF','XC'
write(id_stats,'(206(''#''))')

end subroutine CREATE_STATS

subroutine WRITE_RESULTS(i,j,k,l,dmin,ecc,theta,RSw,&
&dR_abs,dR_rel,dV_abs,dV_rel,dEn_rel,dSMA_abs,RSw_diff,MF,XC,callsArr,istepsArr,&
&mordArr,calls_end,mord_end,isteps_end)
! Description:
!    Writes the results from the ensemble of propagations to disk. The following
!    results are saved: relative energy residuals, absolute and relative
!    position residuals, absolute semi-major axis residuals, integration
!    statistics.
!
! ==============================================================================

use SETTINGS, only: n_rs
implicit none
! Arguments IN
integer,intent(in)   ::  i,j,k,l
integer,intent(in)   ::  MF(1:3),XC,callsArr(1:2),calls_end,istepsArr(1:2),isteps_end
real(qk),intent(in)  ::  dmin,ecc,theta,RSw
real(dk),intent(in)  ::  RSw_diff,dR_abs(1:3),dR_rel(1:3),dV_abs(1:3),dV_rel(1:3)
real(dk),intent(in)  ::  dEn_rel(1:3),dSMA_abs(1:3),mordArr(1:2),mord_end
! Format strings
character(len=*),parameter  :: fmt_en = '(4(i3,'',''),4(g14.7,'',''),3(es14.7,'',''),3i1,'','',i3.3)'
character(len=*),parameter  :: fmt_ar = '(4(i3,'',''),4(g14.7,'',''),6(es14.7,'',''),3i1,'','',i3.3)'
character(len=*),parameter  :: fmt_st = '(4(i3,'',''),4(g14.7,'',''),3(i9,'',''),3(g14.7,'',''),'//&
&'3(i10,'',''),g14.7,'','',3i1,'','',i3.3)'
character(len=*),parameter  :: fmt_enH = '(4(i3,'',''),4(g14.7,'',''),2(a14,'',''),es14.7,'','',3i1,'','',i3.3)'
character(len=*),parameter  :: fmt_arH = '(4(i3,'',''),4(g14.7,'',''),2(2(a14,'',''),es14.7,'',''),3i1,'','',i3.3)'
character(len=*),parameter  :: fmt_stH = '(4(i3,'',''),4(g14.7,'',''),2(a9,'',''),i9,'','',2(a14,'',''),g14.7,'','','//&
&'2(a10,'',''),i10,'','',g14.7,'','',3i1,'','',i3.3)'
character(len=3),parameter  :: nans(1:2) = ['NaN','NaN']

if (n_rs /= 0) then
  write(id_enres,fmt_en) i,j,k,l,dmin,ecc,theta,RSw,dEn_rel,MF,XC
  write(id_ares,fmt_ar)  i,j,k,l,dmin,ecc,theta,RSw,dR_abs,dV_abs,MF,XC
  write(id_rres,fmt_ar)  i,j,k,l,dmin,ecc,theta,RSw,dR_rel,dV_rel,MF,XC
  write(id_sma,fmt_en)   i,j,k,l,dmin,ecc,theta,RSw,dSMA_abs,MF,XC
  write(id_stats,fmt_st)&
  &i,j,k,l,dmin,ecc,theta,RSw,callsArr,calls_end,mordArr,mord_end,istepsArr,isteps_end,RSw_diff,MF,XC

else
  write(id_enres,fmt_enH) i,j,k,l,dmin,ecc,theta,RSw,nans,dEn_rel(3),MF,XC
  write(id_ares,fmt_arH)  i,j,k,l,dmin,ecc,theta,RSw,nans,dR_abs(3),nans,dV_abs(3),MF,XC
  write(id_rres,fmt_arH)  i,j,k,l,dmin,ecc,theta,RSw,nans,dR_rel(3),nans,dV_rel(3),MF,XC
  write(id_stats,fmt_stH)&
  &i,j,k,l,dmin,ecc,theta,RSw,nans,calls_end,nans,mord_end,nans,isteps_end,RSw_diff,MF,XC

end if

end subroutine WRITE_RESULTS

subroutine WRITE_EXCEPTION(i,j,k,l,dmin,ecc,theta,RSw,MF,XC)
! Description:
!    When an exception takes place in the main loops, this subroutine fills the
!    corresponding entry in the output files with NaNs and the error codes of
!    the exception.
!
! ==============================================================================

implicit none
! Arguments IN
integer,intent(in)  ::  i,j,k,l,MF(1:3),XC
real(qk),intent(in) ::  dmin,ecc,theta,RSw
! Locals

write(id_enres,'(4(i3,'',''),4(g14.7,'',''),3(a14,'',''),3i1,'','',i3.3)')&
& i,j,k,l,dmin,ecc,theta,RSw,'NaN','NaN','NaN',MF,XC
write(id_ares,'(4(i3,'',''),4(g14.7,'',''),6(a14,'',''),3i1,'','',i3.3)')&
& i,j,k,l,dmin,ecc,theta,RSw,'NaN','NaN','NaN','NaN','NaN','NaN',MF,XC
write(id_rres,'(4(i3,'',''),4(g14.7,'',''),6(a14,'',''),3i1,'','',i3.3)')&
& i,j,k,l,dmin,ecc,theta,RSw,'NaN','NaN','NaN','NaN','NaN','NaN',MF,XC
write&
&(id_stats,'(4(i3,'',''),4(g14.7,'',''),3(a9,'',''),3(a14,'',''),3(a10,'',''),a14,'','',3i1,'','',i3.3)')&
&i,j,k,l,dmin,ecc,theta,RSw,'NaN','NaN','NaN','NaN','NaN','NaN','NaN','NaN','NaN','NaN',MF,XC

end subroutine WRITE_EXCEPTION

end module IO

max_stars_repo_name: dkamato/NAPLES

max_issues_repo_name: dkamato/NAPLES

max_forks_repo_name: dkamato/NAPLES