text:
c ***********************************************
c    tree.f copy from pv.32.17
c  18 Jan. 2011   written by D. Kawata
c ***********************************************
c70*******************************************************************

c ***********************************************
c    Definition of function about octtree
c  This program build octtree and
c    compute the center of mass and total mass and
c ***********************************************

      subroutine treebuild(np)
      include 'define.f'
      integer np
      integer i,pn,nd,level,npn
c * for not-finished particle (nfp) *
c * number of not-finished particle *
      integer numnfp
c * for subcell (sc) *
c * start subcell in tree *
      integer stsc
c * number of subcells *
      integer numsc,clist(0:MNS-1)
c * diameter of subcell *
      double precision l_sc
c * for work *
      integer numtmp
c *** for calculating hm_tr ***
      integer d,num,npare,pnd,j

c *** for Tree ***
      tcx(0) = 1.0d0
      tcx(1) = 1.0d0
      tcx(2) = 1.0d0
      tcx(3) = 1.0d0
      tcx(4) = -1.0d0
      tcx(5) = -1.0d0
      tcx(6) = -1.0d0
      tcx(7) = -1.0d0
      tcy(0) = 1.0d0
      tcy(1) = 1.0d0
      tcy(4) = 1.0d0
      tcy(5) = 1.0d0
      tcy(2) = -1.0d0
      tcy(3) = -1.0d0
      tcy(6) = -1.0d0
      tcy(7) = -1.0d0
      tcz(0) = 1.0d0
      tcz(2) = 1.0d0
      tcz(4) = 1.0d0
      tcz(6) = 1.0d0
      tcz(1) = -1.0d0
      tcz(3) = -1.0d0
      tcz(5) = -1.0d0
      tcz(7) = -1.0d0

c *** Make root ***
      call makeroot(np)
      next_tr(0)=0
      pare_tr(0)=-1
      num_tr=1
      if(np.eq.0) then
        np_tr(0)=0
        num_tr=0
        return
      endif

c * initialization *
      numnfp = np

      do i=0,np-1
        pn_nfp(i)=i
        nd_nfp(i)=0
        label_nfp(i)=0
      enddo
      numsc=8
      stsc=1
      l_sc=l_tr(0)*0.5d0
      do i=0,7
        np_sc(i)=0
        pare_sc(i)=0
        c_sc(i)=i
      enddo
      flag_pc(0)=0
c * in subcell *
      daughter_tr(0) = 0
      level = 0
      list(0) = 0
c *****   start iteration in level *****
   77 if(numnfp.le.1) then
        goto 99
      endif
c      write(6,*) numnfp,level,l_sc
      if(level.gt.MAXNODE) then
        write(6,*) ' Error in tree(): failure in treebuild'
        stop
      endif
c * initialization for the case np_tr < 8 *
      do i=list(level),num_tr-1
        nd_sc(i)=0
      enddo
      level=level+1
      list(level) = num_tr
      l_sc = 0.5d0*l_sc
c * find out which subcell it is in *
      do i=0,numnfp-1
        pn=pn_nfp(i)
        nd=nd_nfp(i)
        tx(i)=x_p(pn)-cx_tr(nd)
        ty(i)=y_p(pn)-cy_tr(nd)
        tz(i)=z_p(pn)-cz_tr(nd)
      enddo
      do i=0,numnfp-1
        c_nfp(i)=0
        if(tx(i).lt.0.0d0) then
          c_nfp(i)=c_nfp(i)+4
        endif
        if(ty(i).lt.0.0d0) then
          c_nfp(i)=c_nfp(i)+2
        endif
        if(tz(i).lt.0.0d0) then
          c_nfp(i)=c_nfp(i)+1
        endif
      enddo
      npn = 0
      do i=0,numnfp-1
        if(np_tr(nd_nfp(i)).gt.8) then
          nd_nfp(i) = daughter_tr(nd_nfp(i))+c_nfp(i)
        else
          nalist(npn)=nd_nfp(i)
          talist(npn)=i
          npn = npn+1
        endif
      enddo

!ocl novrec (nd_nfp)
!cdir nodep (nd_nfp)
      do i=0,npn-1
        nd_nfp(talist(i))=daughter_tr(nalist(i))+nd_sc(nalist(i))
        nd_sc(nalist(i))=nd_sc(nalist(i))+1
      enddo

      do i=0,npn-1
        c_sc(nd_nfp(talist(i)))=c_nfp(talist(i))
      enddo

c * update info of subcell *
      do i=0,numnfp-1
        np_sc(nd_nfp(i))=np_sc(nd_nfp(i))+1
      enddo

      if(num_tr.gt.MAXNODE-numnfp) then
        write(6,*) ' Error in buildtree():Node is overflow!'
        write(6,*) ' This level is ',level,'numnfp=',numnfp
        stop
      endif

c *** subcell is conected to tree ***
      npn = 0
      do i=0,numsc-1
        if(np_sc(i).ge.1) then
          nalist(npn)=i
          talist(npn)=num_tr
          node(npn)=pare_sc(i)
          clist(npn)=c_sc(i)
          num_tr=num_tr+1
          npn=npn+1
        endif
      enddo

      do i=0,npn-1
        if(flag_pc(node(i)).eq.0) then
          daughter_tr(node(i))=talist(i)
          flag_pc(node(i))=1
        endif
      enddo
!ocl novrec (cx_tr,cy_tr,cz_tr)
!cdir nodep (cx_tr,cy_tr,cz_tr)
      do i=0,npn-1
        nd_sc(nalist(i))=talist(i)
c *****  Set subcell in tree   *****
        np_tr(talist(i))=np_sc(nalist(i))
        l_tr(talist(i))=2.0d0*l_sc
        pare_tr(talist(i))=node(i)
        cx_tr(talist(i))=tcx(clist(i))*l_sc+cx_tr(node(i))
        cy_tr(talist(i))=tcy(clist(i))*l_sc+cy_tr(node(i))
        cz_tr(talist(i))=tcz(clist(i))*l_sc+cz_tr(node(i))
      enddo

c *** Set label not-finished particle ***
      do i=0,numnfp-1
        nd_nfp(i)=nd_sc(nd_nfp(i))
        pn_tr(nd_nfp(i)) = pn_nfp(i)
c *  this node is leaf *
        if(np_tr(nd_nfp(i)).eq.1) then
          label_nfp(i)=1
        endif
      enddo
c *** rebuild not finished particle list ***
      numtmp = 0
!ocl novrec (label_nfp,pn_nfp,nd_nfp)
!cdir nodep (label_nfp,pn_nfp,nd_nfp)
      do i=0,numnfp-1
        if(label_nfp(i).eq.0) then
          pn_nfp(numtmp)=pn_nfp(i)
          nd_nfp(numtmp)=nd_nfp(i)
          label_nfp(numtmp)=0
          numtmp=numtmp+1
        endif
      enddo
      numnfp=numtmp
c *** rebuild subcell ***
      numtmp = num_tr-stsc
      numsc = 0
!ocl novrec (np_sc,nd_sc,pare_sc,c_sc)
!cdir nodep (np_sc,nd_sc,pare_sc,c_sc)
      do i=0,numtmp-1
        if(np_tr(stsc).ge.2) then
          daughter_tr(stsc)=numsc
          flag_pc(stsc) = 0
          np_sc(numsc)=0
          nd_sc(numsc)=0
          pare_sc(numsc)=stsc
          c_sc(numsc) = 0
          np_sc(numsc+1)=0
          nd_sc(numsc+1)=0
          pare_sc(numsc+1)=stsc
          c_sc(numsc+1) = 1
          np_sc(numsc+2)=0
          nd_sc(numsc+2)=0
          pare_sc(numsc+2)=stsc
          c_sc(numsc+2) = 2
          np_sc(numsc+3)=0
          nd_sc(numsc+3)=0
          pare_sc(numsc+3)=stsc
          c_sc(numsc+3) = 3
          np_sc(numsc+4)=0
          nd_sc(numsc+4)=0
          pare_sc(numsc+4)=stsc
          c_sc(numsc+4) = 4
          np_sc(numsc+5)=0
          nd_sc(numsc+5)=0
          pare_sc(numsc+5)=stsc
          c_sc(numsc+5) = 5
          np_sc(numsc+6)=0
          nd_sc(numsc+6)=0
          pare_sc(numsc+6)=stsc
          c_sc(numsc+6) = 6
          np_sc(numsc+7)=0
          nd_sc(numsc+7)=0
          pare_sc(numsc+7)=stsc
          c_sc(numsc+7) = 7
          numsc=numsc+8
        else if(np_tr(stsc).eq.1) then
          daughter_tr(stsc)=-1
        endif
        stsc=stsc+1
      enddo
      stsc = num_tr
      if(numsc.gt.MAXNODE) then
        write(6,*) ' Error in buildtree() : '
        write(6,*) '  Subcell is overflow	!'
        write(6,*) '  This level is ',level
        stop
      endif
      goto 77
c *** set next node ***
   99 next_tr(0)=0
      do i=0,num_tr-2
        pare_sc(i)=pare_tr(i+1)
      enddo
      do i=1,num_tr-2
        if(pare_sc(i).eq.pare_tr(i)) then
          next_tr(i)=i+1
        else
          next_tr(i) = next_tr(pare_tr(i))
        endif
      enddo
      next_tr(num_tr-1)=next_tr(pare_tr(num_tr-1))

c *** compute mass and center of mass ***
      call compute_mass(level)
      return
      end

c *** Definition of makeroot() ***
      subroutine makeroot(np)
      include 'define.f'
      integer i,np
c * max coordinate *
      double precision max_x,max_y,max_z
c * min coordinate <0 *
      double precision min_x,min_y,min_z
c * max,temp length *
      double precision maxl,tl

c *** Define root node ***
      max_x=-INF
      max_y=-INF
      max_z=-INF
      min_x=INF
      min_y=INF
      min_z=INF
      do i=0,np-1
        if(x_p(i).lt.min_x) then
          min_x = x_p(i)
        endif
        if(y_p(i).lt.min_y) then
          min_y = y_p(i)
        endif
        if(z_p(i).lt.min_z) then
          min_z = z_p(i)
        endif
        if(x_p(i).gt.max_x) then
          max_x = x_p(i)
        endif
        if(y_p(i).gt.max_y) then
          max_y = y_p(i)
        endif
        if(z_p(i).gt.max_z) then
          max_z = z_p(i)
        endif
      enddo
c *** get the maximum and minimum for all the particles ***
c *** maximum ***
c      tdvs(0)=max_x
c      tdvs(1)=max_y
c      tdvs(2)=max_z
c      call MPI_ALLREDUCE(tdvs,tdvr,3,MPI_DOUBLE_PRECISION
c     &   ,MPI_MAX,MPI_COMM_WORLD,ierr)
c      max_x=tdvr(0)
c      max_y=tdvr(1)
c      max_z=tdvr(2)
c *** minimum ***
c      tdvs(0)=min_x
c      tdvs(1)=min_y
c      tdvs(2)=min_z
c      call MPI_ALLREDUCE(tdvs,tdvr,3,MPI_DOUBLE_PRECISION
c     &   ,MPI_MIN,MPI_COMM_WORLD,ierr)
c      min_x=tdvr(0)
c      min_y=tdvr(1)
c      min_z=tdvr(2)
c *** check x range ***
      tl=max_x - min_x
      if(tl.lt.0.0d0) then
        tl = -tl
      endif
      maxl=tl
c *** check y range ***
      tl=max_y - min_y
      if(tl.lt.0.0d0) then
        tl = -tl
      endif
      if(tl.gt.maxl) then
        maxl = tl
      endif
c *** check z range ***
      tl=max_z-min_z
      if(tl.lt.0.0d0) then
        tl = -tl
      endif
      if(tl.gt.maxl) then
        maxl = tl
      endif
c *** Set root node ***
      l_tr(0)=MGROOT*maxl
      cx_tr(0)=(max_x+min_x)*0.5d0
      cy_tr(0)=(max_y+min_y)*0.5d0
      cz_tr(0)=(max_z+min_z)*0.5d0
      np_tr(0)=np

c      write(6,*) ' cx=',cx_tr(0),cy_tr(0),cz_tr(0)
c      write(6,*) ' l =',l_tr(0)
      return
      end

c *** Definition of compute_mass() ***
c * This function compute the center of mass and total mass *
c * tree, center of mass, total mass, node No. *
      subroutine compute_mass(level)
      include 'define.f'
      integer i,level
      integer j,d,npare,num,nd,pnd
      double precision xsd,ysd,zsd

c * compute total mass in node *
      do i=0,num_tr-1
        if(np_tr(i).eq.1) then
          mass_tr(i)=m_p(pn_tr(i))
          cmx_tr(i)=m_p(pn_tr(i))*x_p(pn_tr(i))
          cmy_tr(i)=m_p(pn_tr(i))*y_p(pn_tr(i))
          cmz_tr(i)=m_p(pn_tr(i))*z_p(pn_tr(i))
        else
          mass_tr(i)=0.0d0
          cmx_tr(i)=0.0d0
          cmy_tr(i)=0.0d0
          cmz_tr(i)=0.0d0
        endif
      enddo

c *** initialization for set hm_tr ***
      do i=0,num_tr-1
        if(np_tr(i).eq.1) then
          hm_tr(i)=h_p(pn_tr(i))
        else
          hm_tr(i)=0.0d0
        endif
      enddo

      num=num_tr-1
      do j=level,1,-1
        npare = 0
        do i=num,list(j),-1
          if(pare_tr(i).ne.pare_tr(i-1)) then
            pare_sc(npare) = pare_tr(i)
            nd_sc(npare) = i
            npare=npare+1
          endif
        enddo
        do d=0,7
!ocl novrec (mass_tr,cmx_tr,cmy_tr,cmz_tr)
!cdir nodep (mass_tr,cmx_tr,cmy_tr,cmz_tr)
          do i=0,npare-1
            nd=nd_sc(i)+d
            pnd=pare_sc(i)
            if(pnd.eq.pare_tr(nd).and.nd.le.num) then
              mass_tr(pnd)=mass_tr(pnd)+mass_tr(nd)
              cmx_tr(pnd)=cmx_tr(pnd)+cmx_tr(nd)
              cmy_tr(pnd)=cmy_tr(pnd)+cmy_tr(nd)
              cmz_tr(pnd)=cmz_tr(pnd)+cmz_tr(nd)
              if(hm_tr(nd).gt.hm_tr(pnd)) then
                hm_tr(pnd)=hm_tr(nd)
              endif
            endif
          enddo
        enddo
        num=list(j)-1
      enddo
c * compute center of mass *
      do i=0,num_tr-1
        cmx_tr(i)=cmx_tr(i)/mass_tr(i)
        cmy_tr(i)=cmy_tr(i)/mass_tr(i)
        cmz_tr(i)=cmz_tr(i)/mass_tr(i)
      enddo
c *** compute delta ***
!cdir nodep
      do i=0,num_tr-1
        delta_tr(i)=dsqrt((cmx_tr(i)-cx_tr(i))*(cmx_tr(i)-cx_tr(i))
     &    +(cmy_tr(i)-cy_tr(i))*(cmy_tr(i)-cy_tr(i))
     &    +(cmz_tr(i)-cz_tr(i))*(cmz_tr(i)-cz_tr(i)))
      enddo
c * compute Multipole Momentum *
c * initialization *
!cdir nodep
      do i=0,num_tr-1
        mx_tr(i)=0.0d0
        my_tr(i)=0.0d0
        mz_tr(i)=0.0d0
        mxx_tr(i)=0.0d0
        myy_tr(i)=0.0d0
        mzz_tr(i)=0.0d0
        mxy_tr(i)=0.0d0
        myz_tr(i)=0.0d0
        mzx_tr(i)=0.0d0
      enddo
      num = num_tr-1
      do j=level,1,-1
        npare = 0
        do i=num,list(j),-1
          if(pare_tr(i).ne.pare_tr(i-1)) then
            pare_sc(npare) = pare_tr(i)
            nd_sc(npare) = i
            npare=npare+1
          endif
        enddo
        do d=0,7
!ocl novrec(mxx_tr,myy_tr,mzz_tr,mxy_tr,myz_tr,mzx_tr)
!cdir nodep (mxx_tr,myy_tr,mzz_tr,mxy_tr,myz_tr,mzx_tr)
          do i=0,npare-1
            pnd=pare_sc(i)
            nd=nd_sc(i)+d
            if(pnd.eq.pare_tr(nd).and.nd.le.num) then
              xsd=cmx_tr(nd)-cmx_tr(pnd)
              ysd=cmy_tr(nd)-cmy_tr(pnd)
              zsd=cmz_tr(nd)-cmz_tr(pnd)
              mx_tr(pnd)=mx_tr(pnd)
     &          +(mx_tr(nd)-xsd*mass_tr(nd))
              my_tr(pnd)=my_tr(pnd)
     &          +(my_tr(nd)-ysd*mass_tr(nd))
              mz_tr(pnd)=mz_tr(pnd)
     &          +(mz_tr(nd)-zsd*mass_tr(nd))
              mxx_tr(pnd)=mxx_tr(pnd)
     &          +(mxx_tr(nd)-2.0d0*xsd*mx_tr(nd)
     &          +xsd*xsd*mass_tr(nd))
              myy_tr(pnd)=myy_tr(pnd)
     &          +(myy_tr(nd)-2.0d0*ysd*my_tr(nd)
     &           +ysd*ysd*mass_tr(nd))
              mzz_tr(pnd)=mzz_tr(pnd)
     &          +(mzz_tr(nd)-2.0d0*zsd*mz_tr(nd)
     &           +zsd*zsd*mass_tr(nd))
              mxy_tr(pnd)=mxy_tr(pnd)
     &          +(mxy_tr(nd)-xsd*my_tr(nd)-ysd*mx_tr(nd)
     &           +xsd*ysd*mass_tr(nd))
              myz_tr(pnd)=myz_tr(pnd)
     &          +(myz_tr(nd)-ysd*mz_tr(nd)-zsd*my_tr(nd)
     &           +ysd*zsd*mass_tr(nd))
              mzx_tr(pnd)=mzx_tr(pnd)
     &          +(mzx_tr(nd)-zsd*mx_tr(nd)-xsd*mz_tr(nd)
     &           +zsd*xsd*mass_tr(nd))
            endif
          enddo
        enddo
        num = list(j)-1
      enddo
c      write(6,*) ' num_tr = ',num_tr
c      write(6,*) ' cmx',cmx_tr(120)
c      write(6,*) ' mx,my,mz',mx_tr(120),my_tr(120),mz_tr(120)
c      write(6,*) ' mxx,myy,mzz',mxx_tr(120),myy_tr(120),
c     &  mzz_tr(120)
c      write(6,*) ' mxy,myz,mzx',mxy_tr(120),myz_tr(120),
c     &  mzx_tr(120)
      return
      end

max_stars_repo_name: kawatadaisuke/diskm-buildini

max_issues_repo_name: kawatadaisuke/diskm-buildini

max_forks_repo_name: kawatadaisuke/diskm-buildini


text:
module assert_equals_double_precision_tensor_test
    use iso_varying_string, only: var_str
    use vegetables, only: &
            result_t, &
            test_item_t, &
            assert_equals, &
            assert_not, &
            assert_that, &
            describe, &
            it

    implicit none
    private
    public :: test_assert_eq_double_mat

    character(len=*), parameter :: BOTH_MESSAGE = "Both Message"
    character(len=*), parameter :: SUCCESS_MESSAGE = "Success Message"
    character(len=*), parameter :: FAILURE_MESSAGE = "Failure Message"
contains
    function test_assert_eq_double_mat() result(tests)
        type(test_item_t) :: tests

        tests = describe( &
                "assert_equals with double precision tensors", &
                [ it( &
                        "passes with the same tensor", &
                        check_pass_for_same_tensor) &
                , it( &
                        "fails for tensors with different numbers", &
                        check_fail_for_different_numbers) &
                , it( &
                        "fails for tensors with different sizes", &
                        check_fail_for_different_sizes) &
                ])
    end function

    pure function check_pass_for_same_tensor() result(result_)
        type(result_t) :: result_

        type(result_t) :: example_result
        type(result_t) :: example_result_c
        type(result_t) :: example_result_s
        type(result_t) :: example_result_cc
        type(result_t) :: example_result_cs
        type(result_t) :: example_result_sc
        type(result_t) :: example_result_ss

        associate(example => reshape([1.0d0, 5.0d0, 3.0d0, 7.0d0, 2.0d0, 6.0d0, 4.0d0, 8.0d0], [2, 2, 2]))
            example_result = assert_equals(example, example)
            example_result_c = assert_equals(example, example, BOTH_MESSAGE)
            example_result_s = assert_equals(example, example, var_str(BOTH_MESSAGE))
            example_result_cc = assert_equals( &
                    example, example, SUCCESS_MESSAGE, FAILURE_MESSAGE)
            example_result_cs = assert_equals( &
                    example, example, SUCCESS_MESSAGE, var_str(FAILURE_MESSAGE))
            example_result_sc = assert_equals( &
                    example, example, var_str(SUCCESS_MESSAGE), FAILURE_MESSAGE)
            example_result_ss = assert_equals( &
                    example, example, var_str(SUCCESS_MESSAGE), var_str(FAILURE_MESSAGE))
        end associate

        result_ = &
                assert_that( &
                        example_result%passed(), &
                        example_result%verbose_description(.false.)) &
                .and.assert_that( &
                        example_result_c%passed(), &
                        example_result_c%verbose_description(.false.)) &
                .and.assert_that( &
                        example_result_s%passed(), &
                        example_result_s%verbose_description(.false.)) &
                .and.assert_that( &
                        example_result_cc%passed(), &
                        example_result_cc%verbose_description(.false.)) &
                .and.assert_that( &
                        example_result_cs%passed(), &
                        example_result_cs%verbose_description(.false.)) &
                .and.assert_that( &
                        example_result_sc%passed(), &
                        example_result_sc%verbose_description(.false.)) &
                .and.assert_that( &
                        example_result_ss%passed(), &
                        example_result_ss%verbose_description(.false.))
    end function

    pure function check_fail_for_different_numbers() result(result_)
        type(result_t) :: result_

        type(result_t) :: example_result
        type(result_t) :: example_result_c
        type(result_t) :: example_result_s
        type(result_t) :: example_result_cc
        type(result_t) :: example_result_cs
        type(result_t) :: example_result_sc
        type(result_t) :: example_result_ss

        associate( &
                example1 => reshape([1.0d0, 5.0d0, 3.0d0, 7.0d0, 2.0d0, 6.0d0, 4.0d0, 8.0d0], [2, 2, 2]), &
                example2 => reshape([9.0d0, 13.0d0, 11.0d0, 15.0d0, 10.0d0, 14.0d0, 12.0d0, 16.0d0], [2, 2, 2]))
            example_result = assert_equals(example1, example2)
            example_result_c = assert_equals(example1, example2, BOTH_MESSAGE)
            example_result_s = assert_equals(example1, example2, var_str(BOTH_MESSAGE))
            example_result_cc = assert_equals( &
                    example1, example2, SUCCESS_MESSAGE, FAILURE_MESSAGE)
            example_result_cs = assert_equals( &
                    example1, example2, SUCCESS_MESSAGE, var_str(FAILURE_MESSAGE))
            example_result_sc = assert_equals( &
                    example1, example2, var_str(SUCCESS_MESSAGE), FAILURE_MESSAGE)
            example_result_ss = assert_equals( &
                    example1, example2, var_str(SUCCESS_MESSAGE), var_str(FAILURE_MESSAGE))
        end associate

        result_ = &
                assert_not( &
                        example_result%passed(), &
                        example_result%verbose_description(.false.)) &
                .and.assert_not( &
                        example_result_c%passed(), &
                        example_result_c%verbose_description(.false.)) &
                .and.assert_not( &
                        example_result_s%passed(), &
                        example_result_s%verbose_description(.false.)) &
                .and.assert_not( &
                        example_result_cc%passed(), &
                        example_result_cc%verbose_description(.false.)) &
                .and.assert_not( &
                        example_result_cs%passed(), &
                        example_result_cs%verbose_description(.false.)) &
                .and.assert_not( &
                        example_result_sc%passed(), &
                        example_result_sc%verbose_description(.false.)) &
                .and.assert_not( &
                        example_result_ss%passed(), &
                        example_result_ss%verbose_description(.false.))
    end function

    pure function check_fail_for_different_sizes() result(result_)
        type(result_t) :: result_

        type(result_t) :: example_result
        type(result_t) :: example_result_c
        type(result_t) :: example_result_s
        type(result_t) :: example_result_cc
        type(result_t) :: example_result_cs
        type(result_t) :: example_result_sc
        type(result_t) :: example_result_ss

        associate( &
                example1 => reshape([1.0d0, 5.0d0, 3.0d0, 7.0d0, 2.0d0, 6.0d0, 4.0d0, 8.0d0], [2, 2, 2]), &
                example2 => reshape( &
                    [1.0d0, 5.0d0, 9.0d0, &
                    3.0d0, 7.0d0, 11.0d0, &
                    2.0d0, 6.0d0, 10.0d0, &
                    4.0d0, 8.0d0, 12.0d0], [3, 2, 2]))
            example_result = assert_equals(example1, example2)
            example_result_c = assert_equals(example1, example2, BOTH_MESSAGE)
            example_result_s = assert_equals(example1, example2, var_str(BOTH_MESSAGE))
            example_result_cc = assert_equals( &
                    example1, example2, SUCCESS_MESSAGE, FAILURE_MESSAGE)
            example_result_cs = assert_equals( &
                    example1, example2, SUCCESS_MESSAGE, var_str(FAILURE_MESSAGE))
            example_result_sc = assert_equals( &
                    example1, example2, var_str(SUCCESS_MESSAGE), FAILURE_MESSAGE)
            example_result_ss = assert_equals( &
                    example1, example2, var_str(SUCCESS_MESSAGE), var_str(FAILURE_MESSAGE))
        end associate

        result_ = &
                assert_not( &
                        example_result%passed(), &
                        example_result%verbose_description(.false.)) &
                .and.assert_not( &
                        example_result_c%passed(), &
                        example_result_c%verbose_description(.false.)) &
                .and.assert_not( &
                        example_result_s%passed(), &
                        example_result_s%verbose_description(.false.)) &
                .and.assert_not( &
                        example_result_cc%passed(), &
                        example_result_cc%verbose_description(.false.)) &
                .and.assert_not( &
                        example_result_cs%passed(), &
                        example_result_cs%verbose_description(.false.)) &
                .and.assert_not( &
                        example_result_sc%passed(), &
                        example_result_sc%verbose_description(.false.)) &
                .and.assert_not( &
                        example_result_ss%passed(), &
                        example_result_ss%verbose_description(.false.))
    end function
end module

max_stars_repo_name: everythingfunctional/vegetables

max_issues_repo_name: everythingfunctional/vegetables

max_forks_repo_name: everythingfunctional/vegetables


text:
      double precision function DLNREL (X)
c Copyright (c) 1996 California Institute of Technology, Pasadena, CA.
c ALL RIGHTS RESERVED.
c Based on Government Sponsored Research NAS7-03001.
c>>   1994-10-20 DLNREL Krogh  Changes to use M77CON
c>>   1990-11-27 DLNREL WV Snyder Initial Code.
c
c     Compute LOG(1+X).
c
c     This procedure uses approximation 2707 from Hart, et. al. when
c     1/sqrt(2) < 1+x < sqrt(2), and otherwise uses the Fortran
c     intrinsic logarithm routine.
c
c--D replaces "?": ?LNREL
      double precision X
      double precision P0, P1, P2, P3, P4, Q0, Q1, Q2, Q3, Z, ZSQ
      double precision BOT, TOP
      parameter (P0 = +0.75040 94990 77712 22174 55611 007 D+2 )
      parameter (P1 = -0.13456 69115 05043 02353 18253 537 D+3 )
      parameter (P2 = +0.74137 19213 24860 25127 79336 47  D+2 )
      parameter (P3 = -0.12772 49755 01233 08199 84385     D+2 )
      parameter (P4 = +0.33271 08381 08768 69381 44        D00 )
      parameter (Q0 = +0.37520 47495 38856 11087 27775 374 D+2 )
      parameter (Q1 = -0.79790 28073 71500 48794 39951 583 D+2 )
      parameter (Q2 = +0.56161 26132 11825 72920 58560 36  D+2 )
      parameter (Q3 = -0.14508 68091 85808 26853 62325     D+2 )
c     parameter (Q4 = +1.0                                 D00 )
      parameter (BOT = .70710 67811 86547 52440 D0 - 1.0D0 )
      parameter (TOP = 1.4142 13562 37309 50488 D0 - 1.0D0 )
c
      if (x .lt. bot .or. x .gt. top) then
         dlnrel = log(x+1.0d0)
      else
         z = x / (x+2.0d0)
         zsq = z*z
         dlnrel = z * ((((p4*zsq+p3)*zsq+p2)*zsq+p1)*zsq+p0)/
     1                ((((   zsq+q3)*zsq+q2)*zsq+q1)*zsq+q0)
      end if
      return
      end

ext: f


!-----------------------------------------------------------------------
!
! MODULE: mkba_sweep_module
!> @brief
!> This module contains the 2D and 3D mesh sweep logic using mini-KBA.
!
!-----------------------------------------------------------------------

MODULE mkba_sweep_module

  USE global_module, ONLY: i_knd, r_knd, zero, two, one, half

  USE plib_module, ONLY: ichunk, firsty, lasty, firstz, lastz,         &
    nnested

  USE geom_module, ONLY: nx, hi, hj, hk, ndimen, ny, nz, ndiag, diag, nc

  USE sn_module, ONLY: cmom, nang, mu, eta, xi, w, noct

  USE data_module, ONLY: src_opt, qim

  USE control_module, ONLY: fixup, tolr

  USE thrd_comm_module, ONLY: sweep_recv_bdry, sweep_send_bdry

  IMPLICIT NONE

  PUBLIC :: mkba_sweep

  CONTAINS

  SUBROUTINE mkba_sweep ( ich, id, d1, d2, d3, d4, jd, kd, oct, g, iop,&
    t, nthrd, reqs, szreq, psii, psij, psik, qtot, ec, vdelt, ptr_in,  &
    ptr_out, dinv, flux0, fluxm, jb_in, jb_out, kb_in, kb_out, wmu,    &
    weta, wxi, flkx, flky, flkz, t_xs, fmin, fmax )

!-----------------------------------------------------------------------
!
! 3-D slab mesh sweeper.
!
!-----------------------------------------------------------------------

    INTEGER(i_knd), INTENT(IN) :: ich, id, d1, d2, d3, d4, jd, kd, oct,&
      g, iop, t, nthrd, szreq

    INTEGER(i_knd), DIMENSION(szreq), INTENT(INOUT) :: reqs

    REAL(r_knd), INTENT(IN) :: vdelt

    REAL(r_knd), INTENT(INOUT) :: fmin, fmax

    REAL(r_knd), DIMENSION(nang), INTENT(IN) :: wmu, weta, wxi

    REAL(r_knd), DIMENSION(nang,cmom), INTENT(IN) :: ec

    REAL(r_knd), DIMENSION(nang,ny,nz), INTENT(INOUT) :: psii

    REAL(r_knd), DIMENSION(nang,ichunk,nz), INTENT(INOUT) :: psij,     &
      jb_in, jb_out

    REAL(r_knd), DIMENSION(nang,ichunk,ny), INTENT(INOUT) :: psik,     &
      kb_in, kb_out

    REAL(r_knd), DIMENSION(nx,ny,nz), INTENT(IN) :: t_xs

    REAL(r_knd), DIMENSION(nx,ny,nz), INTENT(INOUT) :: flux0

    REAL(r_knd), DIMENSION(nx+1,ny,nz), INTENT(INOUT) :: flkx

    REAL(r_knd), DIMENSION(nx,ny+1,nz), INTENT(INOUT) :: flky

    REAL(r_knd), DIMENSION(nx,ny,nz+1), INTENT(INOUT) :: flkz

    REAL(r_knd), DIMENSION(nang,ichunk,ny,nz), INTENT(IN) :: dinv

    REAL(r_knd), DIMENSION(cmom-1,nx,ny,nz), INTENT(INOUT) :: fluxm

    REAL(r_knd), DIMENSION(cmom,ichunk,ny,nz), INTENT(IN) :: qtot

    REAL(r_knd), DIMENSION(d1,d2,d3,d4), INTENT(IN) :: ptr_in

    REAL(r_knd), DIMENSION(d1,d2,d3,d4), INTENT(OUT) :: ptr_out
!_______________________________________________________________________
!
!   Local variables
!_______________________________________________________________________

    INTEGER(i_knd) :: ist, jst, jlo, jhi, kst, klo, khi, d, n, ic, i,  &
      j, k, l, ibl, ibr, ibb, ibt, ibf, ibk

    REAL(r_knd) :: sum_hv

    REAL(r_knd), DIMENSION(nang) :: psi, pc, den

    REAL(r_knd), DIMENSION(nang,4) :: hv, fxhv
!_______________________________________________________________________
!
!   Set up the sweep order given octant info.
!_______________________________________________________________________

    IF ( id == 1 ) THEN
      ist = -1
    ELSE
      ist = +1
    END IF

    IF ( jd == 1 ) THEN
      jst = -1; jlo = ny; jhi = 1
    ELSE
      jst = +1; jlo = 1; jhi = ny
    END IF

    IF ( kd == 1 ) THEN
      kst = -1; klo = nz; khi = 1
    ELSE
      kst = +1; klo = 1; khi = nz
    END IF
!_______________________________________________________________________
!
!   Initialize the fixup counter
!_______________________________________________________________________

    fxhv = zero
!_______________________________________________________________________
!
!   Process boundaries in y- and z-directions come from preset boundary
!   conditions (vacuum) or from upstream neighbors. Must call to receive
!   any data. If upstream is boundary, no communication.
!_______________________________________________________________________

    CALL sweep_recv_bdry ( g, jd, kd, iop, t, nthrd, reqs, szreq, nc,  &
      nang, ichunk, ny, nz, jb_in, kb_in )
!_______________________________________________________________________
!
!   Loop over cells along the diagonals. When only 1 diagonal, it's
!   normal sweep order. Otherwise, nested threading performs mini-KBA.
!_______________________________________________________________________

  !$OMP PARALLEL NUM_THREADS(nnested) IF(nnested>1) DEFAULT(SHARED)    &
  !$OMP& FIRSTPRIVATE(fxhv)
    diagonal_loop: DO d = 1, ndiag

  !$OMP DO SCHEDULE(STATIC,1) PRIVATE(n,ic,i,j,k,l,psi,pc,sum_hv,hv,den)
      line_loop: DO n = 1, diag(d)%len
!_______________________________________________________________________
!
!       Set the indices chunk indices and the global i index
!_______________________________________________________________________

        ic = diag(d)%cell_id(n)%ic
        IF ( ist < 0 ) ic = ichunk - ic + 1

        i = (ich-1)*ichunk + ic

        j = diag(d)%cell_id(n)%j
        IF ( jst < 0 ) j = ny - j + 1

        k = diag(d)%cell_id(n)%k
        IF ( kst < 0 ) k = nz - k + 1
!_______________________________________________________________________
!
!       Dummy operation to match real transport code where nx is not
!       always even divided by ichunk
!_______________________________________________________________________

        IF ( i > nx ) THEN
          jb_out(:,ic,k) = zero
          kb_out(:,ic,j) = zero
          CYCLE line_loop
        END IF
!_______________________________________________________________________
!
!       Compute the angular source. Add the MMS contribution if
!       necessary.
!_______________________________________________________________________

        psi = qtot(1,ic,j,k)

        DO l = 2, cmom
          psi = psi + ec(:,l)*qtot(l,ic,j,k)
        END DO

        IF ( src_opt == 3 ) psi = psi + qim(:,i,j,k,oct,g)
!_______________________________________________________________________
!
!       Left/right boundary conditions, always vacuum. Dummy operations
!       mayb be ignored.
!_______________________________________________________________________

        ibl = 0
        ibr = 0
        IF ( i==nx .AND. ist==-1 ) THEN
          psii(:,j,k) = zero
        ELSE IF ( i==1 .AND. ist==1 ) THEN
          SELECT CASE ( ibl )
            CASE ( 0 )
              psii(:,j,k) = zero
            CASE ( 1 )
              psii(:,j,k) = zero
          END SELECT
        END IF
!_______________________________________________________________________
!
!       Top/bottom boundary condtions for chunk. Either received from
!       upstream neighbor or set to zero (always vacuum) if no upstream
!       neighbor.
!_______________________________________________________________________

        ibb = 0
        ibt = 0
        IF ( j == jlo ) THEN
          IF ( jd==1 .AND. lasty ) THEN
            psij(:,ic,k) = zero
          ELSE IF ( jd==2 .AND. firsty ) THEN
            SELECT CASE ( ibb )
              CASE ( 0 )
                psij(:,ic,k) = zero
              CASE ( 1 )
                psij(:,ic,k) = zero
            END SELECT
          ELSE
            psij(:,ic,k) = jb_in(:,ic,k)
          END IF
        END IF
!_______________________________________________________________________
!
!       Front/back boundary condtions. Vacuum at global boundaries, but
!       set to some incoming flux from neighboring proc.
!_______________________________________________________________________

        ibf = 0
        ibk = 0
        IF ( k == klo ) THEN
          IF ( (kd==1 .AND. lastz) .OR. ndimen<3 ) THEN
            psik(:,ic,j) = zero
          ELSE IF ( kd==2 .AND. firstz ) THEN
            SELECT CASE ( ibf )
              CASE ( 0 )
                psik(:,ic,j) = zero
              CASE ( 1 )
                psik(:,ic,j) = zero
            END SELECT
          ELSE
            psik(:,ic,j) = kb_in(:,ic,j)
          END IF
        END IF
!_______________________________________________________________________
!
!       Compute the initial solution
!_______________________________________________________________________

        IF ( vdelt /= zero ) THEN
          pc = ( psi + psii(:,j,k)*mu*hi + psij(:,ic,k)*eta*hj +       &
                 psik(:,ic,j)*xi*hk + ptr_in(:,i,j,k)*vdelt ) *        &
               dinv(:,ic,j,k)
        ELSE
          pc = ( psi + psii(:,j,k)*mu*hi + psij(:,ic,k)*eta*hj +       &
            psik(:,ic,j)*xi*hk ) * dinv(:,ic,j,k)
        END IF
!_______________________________________________________________________
!
!       Compute outgoing edges with diamond difference, no negative flux
!       fixup
!_______________________________________________________________________

        IF ( fixup == 0 ) THEN

          psi = pc

          psii(:,j,k) = two*psi - psii(:,j,k)
          psij(:,ic,k) = two*psi - psij(:,ic,k)
          IF ( ndimen == 3 ) psik(:,ic,j) = two*psi - psik(:,ic,j)
          IF ( vdelt /= zero )                                         &
            ptr_out(:,i,j,k) = two*psi - ptr_in(:,i,j,k)

        ELSE
!_______________________________________________________________________
!
!         Use negative flux fixup. Compute outgoing edges. If negative,
!         set to zero and rebalance. Multi-pass until all negativities
!         eliminated. Initialize counters and determine angles that will
!         need fixup first.
!_______________________________________________________________________

          hv = one
          sum_hv = SUM( hv )

          fixup_loop: DO

            fxhv(:,1) = two*pc - psii(:,j,k)
            fxhv(:,2) = two*pc - psij(:,ic,k)
            IF ( ndimen == 3 ) fxhv(:,3) = two*pc - psik(:,ic,j)
            IF ( vdelt /= zero ) fxhv(:,4) = two*pc - ptr_in(:,i,j,k)

            WHERE ( fxhv < zero ) hv = zero
!_______________________________________________________________________
!
!           Exit loop when all angles are fixed up, i.e., no hv change
!_______________________________________________________________________

            IF ( sum_hv == SUM( hv ) ) EXIT fixup_loop
            sum_hv = SUM( hv )
!_______________________________________________________________________
!
!           Recompute balance equation numerator and denominator and get
!           new cell average flux
!_______________________________________________________________________

            IF ( vdelt /= zero ) THEN
              pc = psi + half * ( psii(:,j,k)*mu*hi*(one+hv(:,1)) +    &
                                  psij(:,ic,k)*eta*hj*(one+hv(:,2)) +  &
                                  psik(:,ic,j)*xi*hk*(one+hv(:,3)) +   &
                                  ptr_in(:,i,j,k)*vdelt*(one+hv(:,4)) )
              den = t_xs(i,j,k) + mu*hi*hv(:,1) + eta*hj*hv(:,2) +     &
                xi*hk*hv(:,3) + vdelt*hv(:,4)
            ELSE
              pc = psi + half * ( psii(:,j,k)*mu*hi*(one+hv(:,1)) +    &
                                  psij(:,ic,k)*eta*hj*(one+hv(:,2)) +  &
                                  psik(:,ic,j)*xi*hk*(one+hv(:,3)) )
              den = t_xs(i,j,k) + mu*hi*hv(:,1) + eta*hj*hv(:,2) +     &
                xi*hk*hv(:,3)
            END IF

            WHERE( den > tolr )
              pc = pc/den
            ELSEWHERE
              pc = zero
            END WHERE

          END DO fixup_loop
!_______________________________________________________________________
!
!         Fixup done, compute edges with resolved center
!_______________________________________________________________________

          psi = pc

          psii(:,j,k) = fxhv(:,1) * hv(:,1)
          psij(:,ic,k) = fxhv(:,2) * hv(:,2)
          IF ( ndimen == 3 ) psik(:,ic,j) = fxhv(:,3) * hv(:,3)
          IF ( vdelt /= zero ) ptr_out(:,i,j,k) = fxhv(:,4) * hv(:,4)

        END IF
!_______________________________________________________________________
!
!       Save edge fluxes (dummy if checks for unused non-vacuum BCs)
!_______________________________________________________________________

        IF ( j == jhi ) THEN
          IF ( jd==2 .AND. lasty ) THEN
            CONTINUE
          ELSE IF ( jd==1 .AND. firsty ) THEN
            IF ( ibb == 1 ) CONTINUE
          ELSE
            jb_out(:,ic,k) = psij(:,ic,k)
          END IF
        END IF

        IF ( k == khi ) THEN
          IF ( kd==2 .AND. lastz ) THEN
            CONTINUE
          ELSE IF ( kd==1 .AND. firstz ) THEN
            IF ( ibf == 1 ) CONTINUE
          ELSE
            kb_out(:,ic,j) = psik(:,ic,j)
          END IF
        END IF
!_______________________________________________________________________
!
!       Compute dummy leakages (not used elsewhere currently)
!_______________________________________________________________________

        flkx(i+id-1,j,k) = flkx(i+id-1,j,k) + ist*SUM( wmu*psii(:,j,k) )
        flky(i,j+jd-1,k) = flky(i,j+jd-1,k) +                          &
                           jst*SUM( weta*psij(:,ic,k) )
        flkz(i,j,k+kd-1) = flkz(i,j,k+kd-1) +                          &
                           kst*SUM( wxi*psik(:,ic,j) )
!_______________________________________________________________________
!
!       Compute the flux moments
!_______________________________________________________________________

        psi = w*psi

        IF ( oct == 1 ) THEN
          flux0(i,j,k) = SUM( psi )
          DO l = 1, cmom-1
            fluxm(l,i,j,k) = SUM( ec(:,l+1)*psi )
          END DO
        ELSE
          flux0(i,j,k) = flux0(i,j,k) + SUM( psi )
          DO l = 1, cmom-1
            fluxm(l,i,j,k) = fluxm(l,i,j,k) + SUM( ec(:,l+1)*psi )
          END DO
        END IF
!_______________________________________________________________________
!
!       Calculate dummy min and max scalar fluxes (not used elsewhere
!       currently)
!_______________________________________________________________________

        IF ( oct == noct ) THEN
          fmin = MIN( fmin, flux0(i,j,k) )
          fmax = MAX( fmax, flux0(i,j,k) )
        END IF
!_______________________________________________________________________
!
!       Finish the loops
!_______________________________________________________________________

      END DO line_loop
  !$OMP END DO

    END DO diagonal_loop

  !$OMP END PARALLEL
!_______________________________________________________________________
!
!   Send data to downstream neighbors
!_______________________________________________________________________

    CALL sweep_send_bdry ( g, jd, kd, iop, t, nthrd, reqs, szreq, nc,  &
      nang, ichunk, ny, nz, jb_out, kb_out )
!_______________________________________________________________________
!_______________________________________________________________________

  END SUBROUTINE mkba_sweep

END MODULE mkba_sweep_module


text:
      SUBROUTINE EBCASC( STRIN, NCHARS, STROUT )
C   EBCASC CONVERTS AN EBCDIC STRING TO AND ASCII STRING.  UNKNOWN CHARACTERS
C ARE REPLACED WITH A QUESTION MARK (?).
C
C  ARGUMENTS:
C     STRIN  - THE STRING OF EBCDIC CHARACTERS TO BE CONVERTED.  STRIN must be
C              DECLARED  CHARACTER.
C     NCHARS - THE NUMBER OF CHARACTERS TO CONVERT.
C     STROUT - THE OUTPUT ARRAY OF ASCII CHARACTERS.  ISTROU MAY BE THE SAME AS
C              STRIN.   STROUT must be type CHARACTER
C
C  WRITTEN AND COPYRIGHTED (C) BY:
C  PAUL HENKART, SCRIPPS INSTITUTION OF OCEANOGRAPHY, 4 MAY 1984
C  ALL RIGHTS ARE RESERVED BY THE AUTHOR.  PERMISSION TO COPY OR REPRODUCE THIS
C  SUBROUTINE, BY COMPUTER OR OTHER MEANS, MAY BE OBTAINED ONLY FROM THE AUTHOR.
C
      INTEGER IEBC(256)
      CHARACTER*(*) strin, strout
      DATA IEBC/
     *  1,2,3,63,9,63,16,63,63,63,11,12,13,14,15,16,17,18,19,
     *  63,63,08,63,24,25,63,63,63,63,63,63,63,63,28,63,63,10,23,27,
     *  63,63,63,63,63,05,06,07,00,00,22,00,63,30,63,04,63,63,63,63,
     *  20,21,00,26,32,63,63,63,63,63,63,63,63,63,63,46,60,40,43,124,
     *  38,63,63,63,63,63,63,63,63,63,33,36,42,41,59,95,45,47,63,63,
     *  63,63,63,63,63,63,63,44,37,95,62,63,63,63,63,63,63,63,63,63,
     *  63,96,58,35,64,39,61,34,63,97,98,99,100,101,102,103,104,105,
     *  63,63,63,63,63,63,63,106,107,108,109,110,
     *  111,112,113,114,63,63,63,63,63,63,
     *  63,126,115,116,117,118,119,120,121,122,
     *  63,63,63,91,63,63,63,63,63,63,63,63,63,63,63,63,63,63,63,93,
     *  63,63,123,65,66,67,68,69,70,71,72,73,63,63,63,63,63,63,125,74,
     *  75,76,77,78,79,80,81,82,63,63,63,63,63,63,92,63,83,84,85,86,
     *  87,88,89,90,63,63,63,63,63,63,48,49,50,51,52,53,54,55,56,57,
     *  124,63,63,63,63,63,63/
C
      i = 1
   10 ichr = IAND(ICHAR(strin(i:i)),255)
      strout(i:i) = CHAR(iebc(ichr))
      i = i + 1
      IF( i .GT. nchars ) RETURN
      GOTO 10
      END

max_stars_repo_name: apthorpe/sioseis

max_issues_repo_name: henkart/sioseis

max_forks_repo_name: henkart/sioseis


text:
!***********************************************************************
!                                                                      *
      REAL(KIND(0.0D0)) FUNCTION FUNL (X, K)
!                                                                      *
!   This  function  evaluates the LK(X) functions using the analytic   *
!   functions defined  in table 5  and equations  (20) and  (21)  of   *
!   Fullerton and Rinker.                                              *
!                                                                      *
!                                           Last update: 09 Oct 1992   *
!                                                                      *
!***********************************************************************
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:04:58   1/ 3/07
!...Modified by Charlotte Froese Fischer
!                     Gediminas Gaigalas  10/05/17
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  DOUBLE
      IMPLICIT NONE
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      INTEGER, INTENT(IN) :: K
      REAL(DOUBLE), INTENT(IN) :: X
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      INTEGER :: K1, I
      REAL(DOUBLE), DIMENSION(6,2) :: F
      REAL(DOUBLE), DIMENSION(3,2) :: G
      REAL(DOUBLE), DIMENSION(2,2) :: H
      REAL(DOUBLE) :: A, B, SUM, XN, X2, SUMG, SUMH, XM
!-----------------------------------------------
!
!
      DATA F/ 2.008188D00, -2.397605D00, 1.046471D00, -3.670660D-01, &
         6.374000D-02, -3.705800D-02, 1.646407D00, -2.092942D00, 9.623100D-01, &
         -2.549600D-01, 1.644040D-01, 0.0D00/
!
      DATA G/ 7.51198D-01, 1.38889D-01, 2.0886D-02, 1.37691D-01, -4.16667D-01, &
          - 9.7486D-02/
!
      DATA H/ -4.44444D-01, -3.472D-03, 4.44444D-01, 1.7361D-02/
!
      DATA A, B/ 2.2D00,  - 1.72D00/
!
      IF (K>=0 .AND. K<=1) THEN
         IF (X <= 2.0D00) THEN
            IF (X == 0.0D00) GO TO 6
!
!   Use rational approximation for X < 2
!
            K1 = K + 1
            SUM = 0.0D00
            XN = 1.0D00
            DO I = 1, 6
               SUM = SUM + XN*F(I,K1)
               XN = XN*X
            END DO
            X2 = X*X
            SUMG = G(1,K1) + X2*(G(2,K1)+X2*G(3,K1))
            SUMH = H(1,K1) + X2*X2*H(2,K1)
            XN = DLOG(X)
            SUMG = XN*(SUMG + XN*SUMH)
            IF (K /= 0) THEN
               SUM = SUM + SUMG
               GO TO 7
            ENDIF
            SUM = SUM + X*SUMG
            GO TO 7
         ENDIF
!
         SUM = A + B/X
         IF (K /= 0) SUM = SUM + (SUM + B/X)/X
         SUM = SUM/X
         XM = -X
         SUM = SUM*DEXP(XM)
         GO TO 7
    6    CONTINUE
         IF (K == 1) GO TO 98
         SUM = F(1,1)
    7    CONTINUE
         FUNL = SUM
         RETURN
!
!   Error section
!
   98    CONTINUE
         WRITE (*, 302)
         STOP
      ENDIF
      WRITE (*, 301)
      STOP
!
  301 FORMAT(/,'FUNL: K must be either 0 or 1')
  302 FORMAT(/,'FUNL: Attempt to calculate function for'/,&
         ' zero argument and K value of 1')
      RETURN
!
      END FUNCTION FUNL

max_stars_repo_name: sylas/grasp-continuum

max_issues_repo_name: compas/grasp-rci-qed

max_forks_repo_name: compas/grasp-rci-qed


text:
  ! ----------------------------------------------------------------------
  ! F11 normalization with fixed 180
  ! ----------------------------------------------------------------------
  ! Check the normalization of F11
  tot  = 0d0; tot2 = 0d0
  do j=1,n_ang
     !                  F11                      (sin phi)                      d phi        2pi
     tot  = tot +  p%F(ilam)%F11(j) * (sin(pi*(real(j)-0.5)/real(n_ang))) * (pi/180.d0) * (2.d0*pi)
     tot2 = tot2 + sin(pi*(real(j)-0.5)/real(n_ang))*pi/180.d0*2.d0*pi
  enddo
  ! print *,'theta integration gives: ',tot,tot/4.d0/pi

  ! do it again, using mu integration
  tot  = 0d0; tot2 = 0d0
  do j=1,n_ang-1
     !                  F11                      (sin phi)                      d phi        2pi
     tot  = tot +  2.d0*pi*p%F(ilam)%F11(j) * (cos(pi*(real(j)/real(n_ang))) - cos(pi*(real(j-1)/real(n_ang))))
  enddo
  ! print *,'mu    integration gives: ',tot,tot/4.d0/pi

     ! ----------------------------------------------------------------------
     ! F11 normalization after nang variable
     ! ----------------------------------------------------------------------
     ! Check the normalization of F11
     tot  = 0d0; tot2 = 0d0
     do j=1,nang
        !                  F11                      (sin phi)                      d phi        2pi
        tot  = tot  +  p%F(ilam)%F11(j) * (sin(pi*(real(j)-0.5)/real(nang))) * (pi/dble(nang)) * (2.d0*pi)
        tot2 = tot2 +                      sin(pi*(real(j)-0.5)/real(nang))  * (pi/dble(nang)) * (2.d0*pi)
     enddo
     print *,'theta integration gives: ',tot,tot/4.d0/pi

     ! do it again, using mu integration
     tot  = 0d0; tot2 = 0d0
     do j=1,nang-1
        !               2pi          F11                                     d mu
        tot  = tot +  2.d0*pi*p%F(ilam)%F11(j) * (cos(pi*(real(j)/real(nang))) - cos(pi*(real(j-1)/real(nang))))
     enddo
     print *,'mu    integration gives: ',tot,tot/4.d0/pi

max_stars_repo_name: cdominik/optool

max_issues_repo_name: cdominik/optool

max_forks_repo_name: cdominik/optool


text:
Carol’s Book Store, previously owned and operated by Carol McNeal in several South Sacramento locations, was purchased by Sharon Wright in late 2006 and will continue the bookstores legacy of supporting AfricanAmericans black culture in the Sacramento area. The stores red brick building has about 2,000square feet of retail space offering literature by and for African Americans and addressing such subjects as children, selfhelp, biography, special interests, fiction and local authors. Wright will continue the stores legacy of book signing events, spoken word performances and youth activities; the bookstore also offers a computer kiosk to allow customers to order books online and will carry literature in various languages to address its diverse clientele.

20110222 11:48:54 nbsp Hallo, can I send you some italian books of poetry and novel from my publishing house RUPE MUTEVOLE EDIZIONI? Contact me to s.denti1@tin.it
Thamk you so much, Silvia Users/divinafollia

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
SUBROUTINE beamdis(loads,nf,ratmax,interp,nels,ell,argv,nlen,ips)
!
! This subroutine produces a PostScript output file "*.dis" displaying
! the deformed 1-D finite element mesh.
!
 IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 REAL(iwp),INTENT(IN)::loads(0:),ratmax,ell(:)
 INTEGER,INTENT(IN)::ips,nf(:,:),nels,interp,nlen
 REAL(iwp)::width,height,scale=72,sxy,xo,yo,x,y,dismag,vmax
 REAL(iwp)::xmin,xmax,ymin,ymax,xnow,dmax,zero=0.0_iwp,pt5=0.5_iwp,       &
   opt5=1.5_iwp,fpt5=5.5_iwp,d8=8.0_iwp,ept5=8.5_iwp,d11=11.0_iwp,        &
   one=1.0_iwp,two=2.0_iwp,thr=3.0_iwp,                                   &
   localx,globalx,na,nb,nc,nd,ll,wold,wnew
 INTEGER::i,j,nn
 CHARACTER(LEN=15)::argv
 REAL(iwp),ALLOCATABLE::xcoord(:)
 OPEN(ips,FILE=argv(1:nlen)//'.dis')
!
 nn=nels+1
 ALLOCATE(xcoord(nn))

 xmin=zero
 xmax=zero
 xnow=zero
 ymin=zero
 ymax=zero

 xcoord(1)=xnow

 DO i=2,nn
   xnow=xnow+ell(i-1)
   xcoord(i)=xnow
 END DO

 xmax=xcoord(nn)

 DO i=1,nn
    IF(loads(nf(1,i))<ymin)ymin=loads(nf(1,i))
    IF(loads(nf(1,i))>ymax)ymax=loads(nf(1,i))
 END DO

 width=xmax-xmin
 height=ymax-ymin
 dmax=ratmax*width
 IF(height>width)dmax=ratmax*height
!
 vmax=zero
 DO i=1,nn
     IF(ABS(loads(nf(1,i)))>vmax)vmax=ABS(loads(nf(1,i)))
 END DO
 dismag=dmax/vmax
!
 ymin=zero
 ymax=zero

 DO i=1,nn
   IF(dismag*loads(nf(1,i))<ymin)                            &
     ymin=dismag*loads(nf(1,i))
   IF(dismag*loads(nf(1,i))>ymax)                            &
     ymax=dismag*loads(nf(1,i))
!
   IF(loads(nf(1,i))<ymin)ymin=loads(nf(1,i))
   IF(loads(nf(1,i))>ymax)ymax=loads(nf(1,i))
 END DO
!
 width =xmax-xmin
 height=ymax-ymin
!
!                       allow 1.5" margin minimum on each side of figure
!
!                       portrait mode
!
 IF(height.GE.d11/ept5*width)THEN
!
!                       height governs the scale
!
   sxy=scale*d8/height
   xo=scale*pt5*(ept5-d8*width/height)
   yo=scale*opt5
 ELSE
!
!                       width governs the scale
!
   sxy=scale*fpt5/width
   xo=scale*opt5
   yo=scale*pt5*(d11-fpt5*height/width)
 END IF
!
!
!                       start PostScript output
!
 WRITE(ips,'(a)')'%!PS-Adobe-1.0'
 WRITE(ips,'(a)')'%%DocumentFonts: none'
 WRITE(ips,'(a)')'%%Pages: 1'
 WRITE(ips,'(a)')'%%EndComments'
 WRITE(ips,'(a)')'/m {moveto} def'
 WRITE(ips,'(a)')'/l {lineto} def'
 WRITE(ips,'(a)')'/s {stroke} def'
 WRITE(ips,'(a)')'/c {closepath} def'
 WRITE(ips,'(a)')'%%EndProlog'
 WRITE(ips,'(a)')'%%Page: 0 1'
 WRITE(ips,'(a)')'gsave'
!
!                       draw the deformed mesh
!
 WRITE(ips,'(2f9.2,a)') xo, yo, ' translate'
 WRITE(ips,'(f9.2,a)') 0.5, ' setlinewidth'

 wnew=loads(nf(1,1))

 DO i=1,nels
    DO j=1,interp

    wold=wnew
    localx=j*ell(i)/interp
    globalx=localx+xcoord(i)
    ll=ell(i)
    na=(one/(ll**3))*(ll**3-thr*ll*localx**2+two*localx**3)
    nb=(one/(ll**2))*(localx*ll**2-two*ll*localx**2+localx**3)
    nc=(one/(ll**3))*(thr*ll*localx**2-two*localx**3)
    nd=(one/(ll**2))*(localx**3-ll*localx**2)
    wnew=na*loads(nf(1,i))+nb*loads(nf(2,i))+nc*loads(nf(1,1+i))+nd*loads(nf(2,i+1))

    x=sxy*((globalx-ell(i)/interp)-xmin)
    y=sxy*(dismag*wold-ymin)
    WRITE(ips,'(2f9.2,a)') x, y,' m'
    x=sxy*(globalx-xmin)
    y=sxy*(dismag*wnew-ymin)
    WRITE(ips,'(2f9.2,a)') x, y,' l'
    WRITE(ips,'(a)')'c s'

    END DO

 END DO
!
 WRITE(ips,'(a)')'grestore'
 WRITE(ips,'(a)')'showpage'
 CLOSE(ips)
!
RETURN
END SUBROUTINE beamdis

max_stars_repo_name: leemargetts/PoreFEM

max_issues_repo_name: leemargetts/PoreFEM

max_forks_repo_name: leemargetts/PoreFEM


text:
!  *********************************************************************
!  *                                                                   *
!  *                        subroutine thin1d                          *
!  *                                                                   *
!  *********************************************************************
!  Mixed Precision Version 1.11
!  Written by Gordon A. Fenton, TUNS, Dec. 8, 1993
!  Latest Update: Jun 9, 1999
!
!  PURPOSE  creates the parameter matrices required for the interior cell
!           subdivisions of LAS2G for the case k1 or k2 = 1 and for LAS3G
!           for k1=k2=1, k1=k3=1, or k2=k3=1.
!
!  Requires:
!   1) from libGAFsim:	DSIFA, DSISL, DCHOL2, DAXPY, DSWAP, DDOT, IDAMAX
!
!  i1, i2, and i3 are the indexes into the global covariance matrix R to be
!  extracted to form the local covariance matrix.
!  n is the number of cells in the subdivision less 1. This is equal to
!  3 for LAS2G and to 7 for LAS3G.
!
!  REVISION HISTORY:
!  1.11	replaced dummy dimensions with a (*) for GNU's compiler (Jun 9/99)
!-----------------------------------------------------------------------------
      subroutine thin1d(R,ir,B,ib,S,is,AS,ias,AI,iai,CS,CI,i1,i2,i3,n,iout,tol)
      real AS(ias,n,*), AI(iai,*), CS(*), CI(*)
      real*8 R(ir,*), B(ib,*), S(is,*)
      real*8 RT(2,2), RI(3,3), DA(3), B1(7,7), B2(7,7)
      integer indT(2), indI(3)

   1  format(a,a,a)
   2  format(a,e13.4)
!							extract R
      RT(1,1) = R(i2,i2)
      RT(1,2) = R(i2,i3)
      RT(2,2) = R(i3,i3)

      RI(1,1) = R(i1,i1)
      RI(1,2) = R(i1,i2)
      RI(1,3) = R(i1,i3)
      RI(2,2) = R(i2,i2)
      RI(2,3) = R(i2,i3)
      RI(3,3) = R(i3,i3)
!							factorize RT (edges)
      call dsifa(RT,2,2,indT,ierr)
      if( ierr .ne. 0 ) then
         write(iout,1)'Error: unable to factorize covariance matrix RT in THIN1D.'
         stop
      endif
!							factorize RI (interior)
      call dsifa(RI,3,3,indI,ierr)
      if( ierr .ne. 0 ) then
         write(iout,1)'Error: unable to factorize covariance matrix RI in THIN1D.'
         stop
      endif

      do 30 j = 1, n
         DA(1) = S(i2,j)
         DA(2) = S(i3,j)
         call dsisl( RT, 2, 2, indT, DA )
         AS(1,j,1) = DA(1)
         AS(2,j,1) = DA(2)
         do 10 i = 1, j
            B1(i,j) = B(i,j) - S(i2,i)*DA(1) - S(i3,i)*DA(2)
  10     continue

         DA(1) = S(i1,j)
         DA(2) = S(i2,j)
         DA(3) = S(i3,j)
         call dsisl( RI, 3, 3, indI, DA )
         AI(1,j) = DA(1)
         AI(2,j) = DA(2)
         AI(3,j) = DA(3)
         do 20 i = 1, j
            B2(i,j) = B(i,j) - S(i1,i)*DA(1) - S(i2,i)*DA(2) - S(i3,i)*DA(3)
  20     continue
  30  continue
!							Cholesky Decomposition
      call dchol2( B1, 7, n, rerr )
      if( rerr .gt. tol ) then
         write(iout,1)'Warning: Cholesky decomposition of THIN1D covariance matrix B1'
         write(iout,2)'         has maximum relative error of ',rerr
      endif

      call dchol2( B2, 7, n, rerr )
      if( rerr .gt. tol ) then
         write(iout,1)'Warning: Cholesky decomposition of THIN1D covariance matrix B2'
         write(iout,2)'         has maximum relative error of ',rerr
      endif
!							store in real*4
      ii = 0
      do 40 j = 1, n
      do 40 i = 1, j
         ii = ii + 1
         CS(ii) = B1(i,j)
         CI(ii) = B2(i,j)
  40  continue

      return
      end


text:
      module penta_functions_mod
	use penta_kind_mod
	implicit none
	contains

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c   Miscellaneous functions
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c----------------------------------------------------------------------
c   This function is used in conjunction with 'zeroin' to find
c    the solution to gamma_e=sum(Z*gamma_i).  Returns difference
c    in fluxes at Er_test, using spline parameters from Er_fit_pass.
c   JL 7/2009
c----------------------------------------------------------------------
c
      real(rknd) function rad_flux(Er_test)
      use penta_kind_mod
      use Er_fit_pass
      implicit none
      real(rknd) :: Er_test, Er_tmp, err
      call pol_int(Er_fit,diff_qg,size(Er_fit),Er_test,Er_tmp,err)
      rad_flux=Er_tmp
      return
      end function rad_flux
c
c-----------------------------------------------------------------------
c  Function zeroin
c  a zero of the function  f(x)  is computed in the interval ax,bx
c
c  input..
c
c  ax     left endpoint of initial interval
c  bx     right endpoint of initial interval
c  f      function subprogram which evaluates f(x) for any x in
c         the interval  ax,bx
c  tol    desired length of the interval of uncertainty of the
c         final result ( .ge. 0.0)
c
c
c  output..
c
c  zeroin abcissa approximating a zero of  f  in the interval ax,bx
c
c
c      it is assumed  that   f(ax)   and   f(bx)   have  opposite  signs
c  without  a  check.  zeroin  returns a zero  x  in the given interval
c  ax,bx  to within a tolerance  4*macheps*abs(x) + tol, where macheps
c  is the relative machine precision.
c      this function subprogram is a slightly  modified  translation  of
c  the algol 60 procedure  zero  given in  richard brent, algorithms for
c  minimization without derivatives, prentice - hall, inc. (1973).
c
c   modified by JL for kind precision and iteration limit 7/09
c
c----------------------------------------------------------------------------
c
      real(rknd) function zeroin(ax,bx,f,tol)
      use penta_kind_mod
      implicit none
      real(rknd) ax, bx, f, tol
      integer(iknd),save :: jfirst
      integer(iknd) :: iter
      integer(iknd), parameter :: max_iter = 1000_iknd
      real(rknd) :: a, b, c, d, e, eps, fa, fb, fc, tol1,
     1   xm, p, q, r, s
      external :: f

      !compute eps, the relative machine precision
      iter=0_iknd
      eps = epsilon(1._rknd)
      tol1 = 1._rknd + eps

      ! initialization
      a = ax
      b = bx
      fa = f(a)
      fb = f(b)

      ! begin step
   20 c = a
      fc = fa
      d = b - a
      e = d
   30 if (dabs(fc) .ge. dabs(fb)) go to 40
      a = b
      b = c
      c = a
      fa = fb
      fb = fc
      fc = fa

      ! convergence test
   40 tol1 = 2._rknd*eps*dabs(b) + 0.5_rknd*tol
      xm = .5_rknd*(c - b)
      iter=iter+1_iknd
      if (iter .gt. max_iter) then
        write(*,*) 'Maximum iterations exceeded in function zeroin:',
     1    max_iter
        stop
      endif
      if (dabs(xm) .le. tol1) go to 90
      if (fb .eq. 0._rknd) go to 90

      !is bisection necessary
      if (dabs(e) .lt. tol1) go to 70
      if (dabs(fa) .le. dabs(fb)) go to 70

      ! is quadratic interpolation possible
      if (a .ne. c) go to 50

      ! linear interpolation
      s = fb/fa
      p = 2._rknd*xm*s
      q = 1._rknd - s
      go to 60

      ! inverse quadratic interpolation
   50 q = fa/fc
      r = fb/fc
      s = fb/fa
      p = s*(2._rknd*xm*q*(q - r) - (b - a)*(r - 1._rknd))
      q = (q - 1._rknd)*(r - 1._rknd)*(s - 1._rknd)

      ! adjust signs
   60 if (p .gt. 0._rknd) q = -q
      p = dabs(p)

      ! is interpolation acceptable
      if ((2._rknd*p) .ge. (3._rknd*xm*q - dabs(tol1*q))) go to 70
      if (p .ge. dabs(0.5_rknd*e*q)) go to 70
      e = d
      d = p/q
      go to 80

      ! bisection
   70 d = xm
      e = d

      !complete step
   80 a = b
      fa = fb
      if (dabs(d) .gt. tol1) b = b + d
      if (dabs(d) .le. tol1) b = b + dsign(tol1, xm)
      fb = f(b)
      if ((fb*(fc/dabs(fc))) .gt. 0._rknd) go to 20
      go to 30

      ! done
   90 zeroin = b
      return
      end function zeroin
c
c-----------------------------------------------------------------------
c   This function calculates the error function for input xin
c-----------------------------------------------------------------------
c
      real(rknd) function derf(xin)
      use penta_kind_mod
      implicit none
      integer(iknd) :: k
      real(rknd) :: eps, pi, xin2, r, er, c0
      real(rknd), intent(in) :: xin

      eps=1.e-15_rknd
      pi=4._rknd*datan(1._rknd)
      xin2=xin*xin
      if (dabs(xin).lt. 3.5_rknd) then
        er=1._rknd
        r=1._rknd
        do k=1,50
            r=r*xin2/(k+0.5_rknd)
            er=er+r
            if (dabs(r).le.dabs(er)*eps) exit
        enddo
        c0=2._rknd/dsqrt(pi)*xin*dexp(-xin2)
        derf=c0*er
      else
        er=1._rknd
        r=1._rknd
        do k=1,12
            r=-r*(k-0.5_rknd)/xin2
            er=er+r
        enddo
        c0=dexp(-xin2)/(dabs(xin)*dsqrt(pi))
        derf=1.0D0-c0*er
        if (xin.lt.0.0) derf=-derf
      endif
      end function derf

c
c-----------------------------------------------------------------------
c   This function returns the integrand of the energy integral.
c       Collision frequencies are calculated using calc_perp_coll_freq.
c     Note the value set below for requested monoenergetic coefficients
c       out of range.
c   JL 7/2009 modified from don's fun123
c-----------------------------------------------------------------------
c
      real(rknd) function intfun(Ka)
      use penta_kind_mod
      use intfun_pass
      use thermal_pass
      use bspline
      use parameter_pass
      implicit none
      !dummy variable
      real(rknd) :: Ka
      !local variables
      real(rknd) :: efield, xa, va, nu_perp_aa, xb, Kb,
     1  nu_perp_a, cmul_K, ff, enrm
      real(rknd), allocatable :: nu_perp_ab(:)
      integer(iknd) :: ispec, svtb

      !normalized velocity and velocity
      xa=dsqrt(Ka)
      va=xa*vta

      !define |Er|/v and coefficient range
      efield=abs_Er/va

      !Define energy dependent collision frequencies
      svtb=size(vtb)
      call calc_perp_coll_freq(va,xa,qa,qa,ma,na,Ka,loglam,
     1  vta,Ka,nu_perp_aa)
      allocate(nu_perp_ab(svtb))

      !loop over field species to get collision frequency
      do ispec=1,svtb
        ! Field species norm. energy and velocity
        xb=xa*(vta/vtb(ispec))
        Kb=xb**2_iknd

        !Define energy dependent collision frequency
        call calc_perp_coll_freq(va,xb,qa,qb(ispec),ma,nb(ispec),Kb,
     1    loglam,vta,Ka,nu_perp_ab(ispec))
      enddo

      !define collisionality
      nu_perp_a=nu_perp_aa + sum(nu_perp_ab)
      cmul_K = nu_perp_a/va
      deallocate(nu_perp_ab)

      !use log interp or linear
      if ( log_interp ) then
        efield = dlog(efield)
        enrm = (efield - emin_pass)
     1    /(emax_pass - emin_pass)
        cmul_K=dlog(cmul_K)
      endif

      !Check for cmul or efield values out of range.  If efield is below the smallest database
      !  efield then use the smallest database efield.
      !Then look up coefficient at efield and cmul values from fits.  If the request
      ! is out of range use the value set below.
      if(efield .lt. emin_pass) efield = emin_pass                   !check for efield too small
      if(cmul_K .ge. cmin_pass .and. cmul_K .le. cmax_pass .and.     !check for out of range
     1       efield .ge. emin_pass .and. efield .le. emax_pass) then

        !interpolate database
        ff = dbs2vl(cmul_K,enrm,kcord,keord,xt_c_pass,xt_e_pass,
     1     num_c_pass,num_e_pass,c_spl_pass)

        !handle log of coefficient
        if( log_coeff ) then
          intfun=dexp(ff - Ka)*Ka**2*(Ka - 2.5_rknd)  !the K^2 instead of sqrt(K) accounts for K^1.5 dependance of the coefficients
     1      **(jval-1_iknd)
        else
          intfun=ff*dexp(-Ka)*Ka**2*(Ka - 2.5_rknd)
     1      **(jval-1_iknd)
        endif  !log coefficient
      else
        !handle out of range requests
        intfun = 0._rknd
      endif !range check
      end function intfun

      end module penta_functions_mod

max_stars_repo_name: joseluisvelasco/STELLOPT

max_issues_repo_name: joseluisvelasco/STELLOPT

max_forks_repo_name: joseluisvelasco/STELLOPT


text:
       SUBROUTINE W3FI47(ILABEL,NLABEL)
C$$$   SUBPROGRAM  DOCUMENTATION  BLOCK
C
C SUBPROGRAM: W3FI47         CONVERT LABEL TO OFF. NO. 85 FORMAT (CRAY)
C   AUTHOR: JONES, R.E.      ORG: W342       DATE: 85-07-31
C
C ABSTRACT: CONVERTS A OFFICE NOTE 85 LABEL IN IBM370 FORMAT
C   TO OFFICE NOTE 85 CRAY FORMAT. ALL EBCDIC CHARACTERS ARE
C   CONVERTED TO ASCII. CONVERTS BINARY OR CODED LABEL.
C
C PROGRAM HISTORY LOG:
C   85-07-31  R.E.JONES
C   89-10-24  R.E.JONES   CONVERT TO CRAY CFT77 FORTRAN
C   02-10-15  VUONG       REPLACED FUNCTION ICHAR WITH MOVA2I
C
C USAGE:  CALL W3FI47(ILABEL,NLABEL)
C
C   INPUT VARIABLES:
C     NAMES  INTERFACE DESCRIPTION OF VARIABLES AND TYPES
C     ------ --------- -----------------------------------------------
C     ILABEL ARG LIST  4 WORDS (32 BYTES) CHARACTERS ARE IN EBCDIC OR
C                      BINARY.
C
C   OUTPUT VARIABLES:
C     NAMES  INTERFACE DESCRIPTION OF VARIABLES AND TYPES
C     ------ --------- -----------------------------------------------
C     NLABEL ARG LIST  4 WORDS (32 BYTES), CHARACTERS ARE IN ASCII OR
C                      BINARY.
C
C ATTRIBUTES:
C   LANGUAGE: CRAY CFT77 FORTRAN
C   MACHINE:  CRAY Y-MP8/832
C
C$$$
C
       CHARACTER*1    ILABEL(32)
       CHARACTER*1    NLABEL(32)
C
C          TEST FOR CODED LABEL, IF SO, CONVERT ALL CHARACTERS
C          TEST FOR EBCDIC C, 195 IN DECIMAL
C
           IF (MOVA2I(ILABEL(7)).EQ.195) THEN
C
             CALL AEA(NLABEL(1),ILABEL(1),32)
C
           ELSE
C
C          BINARY LABEL, CONVERT BYTES 1-8, 21-30 TO ASCII
C
             CALL AEA(NLABEL(1),ILABEL(1),8)
C
C          MOVE BYTES 9 TO 20
C
             DO 10 I = 9,20
               NLABEL(I) = ILABEL(I)
 10          CONTINUE
C
C          CONVERT WASHINGTON TO ASCII
C
             CALL AEA(NLABEL(21),ILABEL(21),10)
C
C          TEST BYTES 31 AND 32 FOR BINARY ZERO, IF NOT ZERO
C          CONVERT TO ASCII
C
             IF (MOVA2I(ILABEL(31)).EQ.0) THEN
               NLABEL(31) = CHAR(0)
             ELSE
               CALL AEA(NLABEL(31),ILABEL(31),1)
             ENDIF
C
             IF (MOVA2I(ILABEL(32)).EQ.0) THEN
               NLABEL(32) = CHAR(0)
             ELSE
               CALL AEA(NLABEL(32),ILABEL(32),1)
             ENDIF
C
           ENDIF
C
         RETURN
       END

max_stars_repo_name: GMAO-SI-Team/jedi-stack

max_issues_repo_name: GMAO-SI-Team/jedi-stack

max_forks_repo_name: GMAO-SI-Team/jedi-stack


text:
Kickballpalooza is an annual kick ball tournament put on by UCD Intramural Sports & Sport Clubs UC Davis Intramural Sports. In 2007 it was scheduled for October 26^th^ and 27^th^.

  The California Aggie http://media.www.californiaaggie.com/media/storage/paper981/news/2007/10/12/CampusNews/Kickballpalooza.Registration.Begins.Monday3029101.shtml article

20081121 12:55:28 nbsp o_O  I missed it?! Users/SunjeetBaadkar


text:
*DECK RADB3
      SUBROUTINE RADB3 (IDO, L1, CC, CH, WA1, WA2)
C***BEGIN PROLOGUE  RADB3
C***SUBSIDIARY
C***PURPOSE  Calculate the fast Fourier transform of subvectors of
C            length three.
C***LIBRARY   SLATEC (FFTPACK)
C***TYPE      SINGLE PRECISION (RADB3-S)
C***AUTHOR  Swarztrauber, P. N., (NCAR)
C***ROUTINES CALLED  (NONE)
C***REVISION HISTORY  (YYMMDD)
C   790601  DATE WRITTEN
C   830401  Modified to use SLATEC library source file format.
C   860115  Modified by Ron Boisvert to adhere to Fortran 77 by
C           (a) changing dummy array size declarations (1) to (*),
C           (b) changing definition of variable TAUI by using
C               FORTRAN intrinsic function SQRT instead of a DATA
C               statement.
C   881128  Modified by Dick Valent to meet prologue standards.
C   890831  Modified array declarations.  (WRB)
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900402  Added TYPE section.  (WRB)
C***END PROLOGUE  RADB3
      DIMENSION CC(IDO,3,*), CH(IDO,L1,3), WA1(*), WA2(*)
C***FIRST EXECUTABLE STATEMENT  RADB3
      TAUR = -.5
      TAUI = .5*SQRT(3.)
      DO 101 K=1,L1
         TR2 = CC(IDO,2,K)+CC(IDO,2,K)
         CR2 = CC(1,1,K)+TAUR*TR2
         CH(1,K,1) = CC(1,1,K)+TR2
         CI3 = TAUI*(CC(1,3,K)+CC(1,3,K))
         CH(1,K,2) = CR2-CI3
         CH(1,K,3) = CR2+CI3
  101 CONTINUE
      IF (IDO .EQ. 1) RETURN
      IDP2 = IDO+2
      IF((IDO-1)/2.LT.L1) GO TO 104
      DO 103 K=1,L1
CDIR$ IVDEP
         DO 102 I=3,IDO,2
            IC = IDP2-I
            TR2 = CC(I-1,3,K)+CC(IC-1,2,K)
            CR2 = CC(I-1,1,K)+TAUR*TR2
            CH(I-1,K,1) = CC(I-1,1,K)+TR2
            TI2 = CC(I,3,K)-CC(IC,2,K)
            CI2 = CC(I,1,K)+TAUR*TI2
            CH(I,K,1) = CC(I,1,K)+TI2
            CR3 = TAUI*(CC(I-1,3,K)-CC(IC-1,2,K))
            CI3 = TAUI*(CC(I,3,K)+CC(IC,2,K))
            DR2 = CR2-CI3
            DR3 = CR2+CI3
            DI2 = CI2+CR3
            DI3 = CI2-CR3
            CH(I-1,K,2) = WA1(I-2)*DR2-WA1(I-1)*DI2
            CH(I,K,2) = WA1(I-2)*DI2+WA1(I-1)*DR2
            CH(I-1,K,3) = WA2(I-2)*DR3-WA2(I-1)*DI3
            CH(I,K,3) = WA2(I-2)*DI3+WA2(I-1)*DR3
  102    CONTINUE
  103 CONTINUE
      RETURN
  104 DO 106 I=3,IDO,2
         IC = IDP2-I
CDIR$ IVDEP
         DO 105 K=1,L1
            TR2 = CC(I-1,3,K)+CC(IC-1,2,K)
            CR2 = CC(I-1,1,K)+TAUR*TR2
            CH(I-1,K,1) = CC(I-1,1,K)+TR2
            TI2 = CC(I,3,K)-CC(IC,2,K)
            CI2 = CC(I,1,K)+TAUR*TI2
            CH(I,K,1) = CC(I,1,K)+TI2
            CR3 = TAUI*(CC(I-1,3,K)-CC(IC-1,2,K))
            CI3 = TAUI*(CC(I,3,K)+CC(IC,2,K))
            DR2 = CR2-CI3
            DR3 = CR2+CI3
            DI2 = CI2+CR3
            DI3 = CI2-CR3
            CH(I-1,K,2) = WA1(I-2)*DR2-WA1(I-1)*DI2
            CH(I,K,2) = WA1(I-2)*DI2+WA1(I-1)*DR2
            CH(I-1,K,3) = WA2(I-2)*DR3-WA2(I-1)*DI3
            CH(I,K,3) = WA2(I-2)*DI3+WA2(I-1)*DR3
  105    CONTINUE
  106 CONTINUE
      RETURN
      END

max_stars_repo_name: homiak/jabaws

max_issues_repo_name: homiak/jabaws

max_forks_repo_name: homiak/jabaws


text:
!
!     CalculiX - A 3-dimensional finite element program
!              Copyright (C) 1998-2021 Guido Dhondt
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation(version 2);
!
!
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!     GNU General Public License for more details.
!
!     You should have received a copy of the GNU General Public License
!     along with this program; if not, write to the Free Software
!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
!
      subroutine nodalthicknesss(inpc,textpart,set,istartset,iendset,
     &  ialset,nset,thickn,nk,istep,istat,n,iline,ipol,inl,ipoinp,
     &  inp,iaxial,ipoinpc,ier)
!
!     reading the input deck: *NODAL THICKNESS
!
      implicit none
!
      character*1 inpc(*)
      character*81 set(*),noset
      character*132 textpart(16)
!
      integer istartset(*),iendset(*),ialset(*),nset,nk,istep,istat,n,
     &  key,i,j,k,l,ipos,iline,ipol,inl,ipoinp(2,*),inp(3,*),iaxial,
     &  ipoinpc(0:*),ier,id
!
      real*8 thickn(2,*),thickness1,thickness2
!
      if(istep.gt.0) then
         write(*,*) '*ERROR reading *NODAL THICKNESS: *NODAL THICKNESS'
         write(*,*) '      should be placed before all step definitions'
         ier=1
         return
      endif
!
      do i=2,n
         write(*,*)
     &    '*WARNING reading *NODAL THICKNESS: parameter not recognized:'
         write(*,*) '         ',
     &        textpart(i)(1:index(textpart(i),' ')-1)
         call inputwarning(inpc,ipoinpc,iline,
     &"*NODAL THICKNESS%")
      enddo
!
      do
         call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
     &        ipoinp,inp,ipoinpc)
         if((istat.lt.0).or.(key.eq.1)) return
         read(textpart(2)(1:20),'(f20.0)',iostat=istat) thickness1
         if(istat.gt.0) then
            call inputerror(inpc,ipoinpc,iline,
     &           "*NODAL THICKNESS%",ier)
            return
         endif
         if(iaxial.eq.180) thickness1=thickness1/iaxial
         if(n.eq.2) then
            thickness2=0.d0
         else
            read(textpart(3)(1:20),'(f20.0)',iostat=istat) thickness2
            if(istat.gt.0) then
               call inputerror(inpc,ipoinpc,iline,
     &              "*NODAL THICKNESS%",ier)
               return
            endif
         endif
         read(textpart(1)(1:10),'(i10)',iostat=istat) l
         if(istat.eq.0) then
            thickn(1,l)=thickness1
            thickn(2,l)=thickness2
         else
            read(textpart(1)(1:80),'(a80)',iostat=istat) noset
            noset(81:81)=' '
            ipos=index(noset,' ')
            noset(ipos:ipos)='N'
c            do i=1,nset
c               if(set(i).eq.noset) exit
c            enddo
            call cident81(set,noset,nset,id)
            i=nset+1
            if(id.gt.0) then
              if(noset.eq.set(id)) then
                i=id
              endif
            endif
            if(i.gt.nset) then
               noset(ipos:ipos)=' '
               write(*,*) '*ERROR reading *NODAL THICKNESS: node set ',
     &                noset
               write(*,*) '  has not yet been defined. '
               call inputerror(inpc,ipoinpc,iline,
     &              "*NODAL THICKNESS%",ier)
               return
            endif
            do j=istartset(i),iendset(i)
               if(ialset(j).gt.0) then
                  thickn(1,ialset(j))=thickness1
                  thickn(2,ialset(j))=thickness2
               else
                  k=ialset(j-2)
                  do
                     k=k-ialset(j)
                     if(k.ge.ialset(j-1)) exit
                     thickn(1,k)=thickness1
                     thickn(2,k)=thickness2
                  enddo
               endif
            enddo
         endif
      enddo
!
      return
      end

max_stars_repo_name: alleindrach/calculix-desktop

max_issues_repo_name: alleindrach/calculix-desktop

max_forks_repo_name: alleindrach/calculix-desktop


text:
      subroutine ckfrsr(cbsr,csts,jcode,jflag,nbaspd,nbtd,nbtmax,
     $ nbt1mx,noutpt,nrct,nrctmx,nrndex,nsrtmx,nstmax,nsts,nstsmx,
     $ nstsr,nttyo,ureac,uspec)
c
c     This subroutine checks the reactions for special reactants. As
c     needed, it rewrites these reactions to eliminate any basis
c     species for which jflag = 30. Such species are not in the active
c     basis set, and there are no corresonding mass balances.
c
c     This subroutine is called by:
c
c       EQ6/eq6.f
c
c-----------------------------------------------------------------------
c
c     Principal input:
c
c       cbsr   = array of coefficients for reactions for special
c                  reactants
c       nbtmax = the maximum number of basis species
c       nbt1mx = the maximum number of basis species plus 1
c       nsrtmx = the maximum number of special reactants
c       ureac  = array of names of reactants
c       uspec  = array of names of species
c
c     Principal output:
c
c       cbsr   = array of coefficients for reactions for special
c                  reactants (modified if necessary)
c
c-----------------------------------------------------------------------
c
      implicit none
c
c-----------------------------------------------------------------------
c
c     Calling sequence variable declarations.
c
      integer nbtmax,nbt1mx,nrctmx,nsrtmx,nstmax,nstsmx
c
      integer noutpt,nttyo
c
      integer jcode(nrctmx),jflag(nstmax),nbaspd(nbtmax),nrndex(nrctmx),
     $ nsts(nstsmx),nstsr(2,nstmax)
c
      integer nbtd,nrct
c
      character*48 uspec(nstmax)
      character*24 ureac(nrctmx)
c
      real*8 cbsr(nbt1mx,nsrtmx),csts(nstsmx)
c
c-----------------------------------------------------------------------
c
c     Local variable declarations.
c
      integer nf
c
      integer jlen,j2,n,nb,nbb,nrc,nr1,nr2,ns,nsr
c
      integer ilnobl
c
      logical qcaught
c
      character*56 uspn56
c
      real*8 cx
c
c-----------------------------------------------------------------------
c
      do nrc = 1,nrct
        if (jcode(nrc) .eq. 2) then
          nsr = nrndex(nrc)
          qcaught = .false.
          do nb = 1,nbtd
            ns = nbaspd(nb)
            cx = cbsr(nb,nsr)
            if (cx.ne.0. .and. jflag(ns).eq.30) then
              if (.not.qcaught) then
c
c               Write the original reaction.
c
                j2 = ilnobl(ureac(nrc))
                write (noutpt,1000) ureac(nrc)(1:j2)
                write (nttyo,1000) ureac(nrc)(1:j2)
 1000           format(/1x,'The reaction for special reactant ',a,
     $          ' is written in terms of',/1x,'one or more species',
     $          ' which are not in the active basis set. The',/1x,
     $          'currently written reaction is:')
c
                nf = noutpt
                call prrsr(cbsr,jcode,nbaspd,nbtd,nbtmax,nbt1mx,
     $          nf,noutpt,nrc,nrctmx,nrndex,nsrtmx,nstmax,nttyo,
     $          ureac,uspec)
                write (noutpt,1010)
 1010           format(/1x,'The following species are not in the',
     $          ' active basis set and will be eliminated',
     $          /1x,'from the reaction:',/)
c
                nf = nttyo
                call prrsr(cbsr,jcode,nbaspd,nbtd,nbtmax,nbt1mx,
     $          nf,noutpt,nrc,nrctmx,nrndex,nsrtmx,nstmax,nttyo,
     $          ureac,uspec)
                write (nttyo,1010)
c
                qcaught = .true.
              endif
c
c             Calling sequence substitutions:
c               uspec(ns) for unam48
c
              call fmspnx(jlen,uspec(ns),uspn56)
c
              write (noutpt,1020) uspn56(1:jlen)
              write (nttyo,1020) uspn56(1:jlen)
 1020         format(3x,a)
c
c             The reaction contains a non-zero reaction coefficient for
c             a species which is not in the active basis set. Rewrite
c             the reaction so that this species does not appear.
c
              nr1 = nstsr(1,ns)
              nr2 = nstsr(2,ns)
              do n = nr1,nr2
                nbb = nsts(n)
                cbsr(nbb,nsr) = cbsr(nbb,nsr) + csts(n)*cx
              enddo
              cbsr(nb,nsr) = 0.
            endif
          enddo
c
          if (qcaught) then
c
c           Write the modified reaction.
c
            write (noutpt,1070)
            write (nttyo,1070)
 1070       format(//1x,'The modified reaction is:')
c
            nf = noutpt
            call prrsr(cbsr,jcode,nbaspd,nbtd,nbtmax,nbt1mx,
     $      nf,noutpt,nrc,nrctmx,nrndex,nsrtmx,nstmax,nttyo,
     $      ureac,uspec)
            write (noutpt,1030)
 1030       format(/1x)
c
            nf = nttyo
            call prrsr(cbsr,jcode,nbaspd,nbtd,nbtmax,nbt1mx,
     $      nf,noutpt,nrc,nrctmx,nrndex,nsrtmx,nstmax,nttyo,
     $      ureac,uspec)
            write (nttyo,1030)
          endif
c
        endif
      enddo
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
      end


text:
C# Make the comment for :
C# CALL CPUTIM, CALL ELAPSE
C#
      PROGRAM TRFD
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      REAL CPU,WALL,CPU1,WALL1,TIM1,TIM2
      LOGICAL OUT
      COMMON/IOFILE/IR,IW
      COMMON/IJPAIR/IA(255)
      COMMON/MEMORY/X(1000000)
      COMMON/TRFPAR/NUM,NORB
      DATA MAX /1000000/
C
C#      CALL CPUTIM(CPU)
C#      CALL ELAPSE(WALL)
      OPEN(1,FILE='TIV')
      OPEN(6,FILE='TIO6')
      WRITE(1,1)
1     FORMAT(1X,'OUTPUT FOR PERFECT CLUB BENCHMARK: TRFD'/
     *       1X,'$Revision: 1.5 $ $Author: kipp $'/)

      IR=5
      IW=6
      OUT=.FALSE.
C     CALL XUFLOW(0)
C
C     ----- SET MO'S AND AO INTEGRALS -----
C
      DO 30 NUM=10,40,5
C
      DO 10 I=1,NUM
   10 IA(I)=(I*(I-1))/2
      NORB=NUM
C
      NP=NUM
      NPQ=(NUM*(NUM+1))/2
      NRS=(NUM*(NUM+1))/2
      I00=1
      I10=I00+NUM*NUM
      I20=I10+NUM
      I30=I20+NORB*NUM
      I40=I30+NPQ*NRS
      LAST=I40
      NEED=LAST-1
      IF(NEED.GT.MAX) WRITE(IW,9998) NUM,NEED,MAX
 9998 FORMAT(' NOT ENOUGH CORE. NUM,NEED,MAX = ',I4,2I10)
      IF(NEED.GT.MAX) STOP
C
C     ----- CALCULATE -V- AND (IJ//KL) -----
C
      CALL INTGRL(NP,NPQ,X(I00),X(I30))
      IF(OUT) CALL TRFPRT(NP,NPQ,IA,X(I30),0)
C
C     ----- TRANSFORM INTEGRALS -----
C
      TIME=0.0D+00
      NTRF=1
      DO 20 ITRF=1,NTRF
C#      CALL CPUTIM(TIM1)
      CALL OLDA(X(I00),X(I30),X(I30),X(I30),X(I30),
     1                 X(I20),X(I20),X(I10),X(I10),
     2                 NORB,NORB,NUM,NUM)
C#      CALL CPUTIM(TIM2)
      WRITE(IW,*) '  AFTER OLDA  ',X(I30)
      TIME=TIME+(TIM2-TIM1)
   20 CONTINUE
      NIJ=(NORB*(NORB+1))/2
      IF(OUT) CALL TRFPRT(NORB,NIJ,IA,X(I30),1)
C
C---- MFLOP=NFLOP(NUM)
C---- FLOPS=MFLOP/TIME
      FN=FLOAT(NUM)
      FLOPS=(25.0*FN**5 + 58.0*FN**4 + 47.0*FN**3 + 14.0*FN**2 )/24.0
      FLOPS=FLOPS*2.0
      FLOPS=FLOPS/TIME
      FLOPS=FLOPS/1.0D+06
      WRITE(IW,9999) NUM,TIME,FLOPS,MFLOP
 9999 FORMAT(' NUM = ',I3,' TIME = ',F12.3,
     1       ' (SEC.) . (MFLOPS=',F5.2,' NFLOP=',I15,')')
C
   30 CONTINUE
      WRITE(1,2)
      WRITE(1,5) X(I30)
      IF (X(I30) .GT. 106352.   .AND.   X(I30) .LT. 106353) THEN
        WRITE(1,3)
      ELSE
        WRITE(1,4)
      ENDIF
C
C#      CALL CPUTIM (CPU1)
C#      CALL ELAPSE (WALL1)
      WRITE (IW,9997) CPU1-CPU,MOD(WALL1-WALL+1000000.0,1000000.0)
 9997 FORMAT (
     1 ' TOTAL ELAPSED CPU TIME = ',F15.5,' SECONDS. '/
     2 ' TOTAL ELAPSED WALL CLOCK TIME = ',F15.5,' SECONDS. ')
      WRITE(1,*)
      WRITE(1,*) 'ELAPSED CPU TIME IN SECONDS: ',CPU1-CPU
C MFLOPS = MFLOP 1 CPU CRAY X-MP/CPU TIME
      WRITE(1,*) 'MFLOPS RATE: ', 430.877/(CPU1-CPU)
      WRITE(1,*) 'ELAPSED WALLCLOCK TIME IN SECONDS: ',
     *            MOD(WALL1-WALL+1000000.0,1000000.0)
2     FORMAT(1X,'VALIDATION PARAMETERS:'/)
3     FORMAT(//1X,'RESULTS FOR THIS RUN ARE:   VALID')
4     FORMAT(//1X,'RESULTS FOR THIS RUN ARE: INVALID')
5     FORMAT(1X,E12.5)
      STOP
      END
      SUBROUTINE INTGRL(NP,NPQ,V,X)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C     COMMON/IJPAIR/IA(1)  ***COMMON BLOCK SIZE MUST BE THE SAME
C                          ***IN ALL SUBPROGRAMS (ANSI) 3/7/89
      COMMON/IJPAIR/IA(255)
      DIMENSION V(NP,1),X(NPQ,1)
C
      DO 20 J=1,NP
      DUM1=FLOAT(J)
      DO 10 I=1,NP
      DUM2=FLOAT(I)
      V(I,J)=1.0D+00-1.0D+00/(DUM1+DUM2)
   10 CONTINUE
   20 CONTINUE
C
      DO 140 I=1,NP
      DO 130 J=1,I
      IJ=IA(I)+J
      A=1.0D+00/(I+J)
      DO 120 K=1,I
      MAXL=K
      IF(K.EQ.I) MAXL=J
      DO 110 L=1,MAXL
      KL=IA(K)+L
      B=1.0D+00/(K+L)
      VAL=A+B
      IF(I.EQ.J) VAL=VAL*0.5D+00
      IF(K.EQ.L) VAL=VAL*0.5D+00
      X(IJ,KL)=VAL
      X(KL,IJ)=VAL
  110 CONTINUE
  120 CONTINUE
  130 CONTINUE
  140 CONTINUE
C
      RETURN
      END
      SUBROUTINE TRFA(V,XRSPQ,XRSIJ,XIJRS,XIJKL,
     1 XRSIQ,XIJKS,XIJ,XKL,NORB,MORB,NUM,NDIM)
C
C     ----- THIS IS THE ORIGINAL -TRFCOR- FROM -HONDO- -----
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION V(NDIM,1)
      DIMENSION XRSPQ(1),XRSIJ(1),XIJRS(1),XIJKL(1)
      DIMENSION XRSIQ(MORB,1),XIJKS(MORB,1),XIJ(1),XKL(1)
      DATA ZERO /0.0D+00/
C
C     ----- TRANSFORM -P- AND -Q- -----
C
      NRS=(NUM*(NUM+1))/2
      NP=NUM
      NQ=NUM
      MRSIJ0=0
      MRSPQ=0
C
      DO 100 MRS=1,NRS
C
      DO 10 MQ=1,NQ
      DO 10 MI=1,MORB
   10 XRSIQ(MI,MQ)=ZERO
      DO 40 MP=1,NP
      DO 30 MQ=1,MP
      MRSPQ=MRSPQ+1
      VAL=XRSPQ(MRSPQ)
C---- IF(VAL.EQ.ZERO) GO TO 30
      DO 20 MI=1,MORB
      XRSIQ(MI,MQ)=XRSIQ(MI,MQ)+VAL*V(MP,MI)
      XRSIQ(MI,MP)=XRSIQ(MI,MP)+VAL*V(MQ,MI)
   20 CONTINUE
   30 CONTINUE
   40 CONTINUE
C
      MRSIJ=MRSIJ0
      DO 90 MI=1,MORB
      DO 50 MJ=1,MI
   50 XIJ(MJ)=ZERO
      DO 70 MQ=1,NQ
      VAL=XRSIQ(MI,MQ)
C---- IF(VAL.EQ.ZERO) GO TO 70
      DO 60 MJ=1,MI
   60 XIJ(MJ)=XIJ(MJ)+VAL*V(MQ,MJ)
   70 CONTINUE
      DO 80 MJ=1,MI
      MRSIJ=MRSIJ+1
   80 XRSIJ(MRSIJ)=XIJ(MJ)
   90 CONTINUE
C
      MRSIJ0=MRSIJ0+NRS
  100 CONTINUE
C
C     ----- TRANSPOSE (R,S//I,J) -----
C
      NIJ=(MORB*(MORB+1))/2
      MRSIJ=0
      DO 120 MRS=1,NRS
      MIJRS=0
      MAX=MRS
      IF(MAX.GT.NIJ) MAX=NIJ
      DO 110 MIJ=1,MAX
      DUM=XRSIJ(MRSIJ+MIJ)
      XIJRS(MRSIJ+MIJ)=XRSIJ(MIJRS+MRS)
      XIJRS(MIJRS+MRS)=DUM
  110 MIJRS=MIJRS+NRS
  120 MRSIJ=MRSIJ+NRS
C
C     ----- TRANSFORM -R- AND -S- -----
C
      NR=NUM
      NS=NUM
      MIJKL=0
      MIJRS=0
      MIJ=0
      MLEFT=NRS-NIJ
C
      DO 300 MI=1,MORB
      DO 300 MJ=1,MI
      MIJ=MIJ+1
C
      DO 210 MS=1,NS
      DO 210 MK=MI,MORB
  210 XIJKS(MK,MS)=ZERO
      DO 240 MR=1,NR
      DO 230 MS=1,MR
      MIJRS=MIJRS+1
      VAL=XIJRS(MIJRS)
C---- IF(VAL.EQ.ZERO) GO TO 230
      DO 220 MK=MI,MORB
      XIJKS(MK,MS)=XIJKS(MK,MS)+VAL*V(MR,MK)
      XIJKS(MK,MR)=XIJKS(MK,MR)+VAL*V(MS,MK)
  220 CONTINUE
  230 CONTINUE
  240 CONTINUE
C
      LMIN=MJ
      LMAX=MI
      DO 290 MK=MI,MORB
      DO 250 ML=LMIN,LMAX
  250 XKL(ML)=ZERO
      DO 270 MS=1,NS
      VAL=XIJKS(MK,MS)
C---- IF(VAL.EQ.ZERO) GO TO 270
      DO 260 ML=LMIN,LMAX
  260 XKL(ML)=XKL(ML)+VAL*V(MS,ML)
  270 CONTINUE
      DO 280 ML=LMIN,LMAX
      MIJKL=MIJKL+1
  280 XIJKL(MIJKL)=XKL(ML)
      LMIN=1
      LMAX=MK+1
  290 CONTINUE
      MIJKL=MIJKL+MIJ+MLEFT
C
  300 CONTINUE
      RETURN
      END
      SUBROUTINE OLDA(V,XRSPQ,XRSIJ,XIJRS,XIJKL,
     1 XRSIQ,XIJKS,XIJ,XKL,NORB,MORB,NUM,NDIM)
C
C     ----- THIS IS THE ORIGINAL -TRFCOR- FROM -HONDO- -----
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION V(NDIM,1)
      DIMENSION XRSPQ(1),XRSIJ(1),XIJRS(1),XIJKL(1)
      DIMENSION XRSIQ(MORB,1),XIJKS(MORB,1),XIJ(1),XKL(1)
      DATA ZERO /0.0D+00/
C
C     ----- TRANSFORM -P- AND -Q- -----
C
      NRS=(NUM*(NUM+1))/2
      NP=NUM
      NQ=NUM
      MRSIJ0=0
      MRSPQ=0
C
      DO 100 MRS=1,NRS
C
      DO 10 MQ=1,NQ
      DO 10 MI=1,MORB
   10 XRSIQ(MI,MQ)=ZERO
      DO 40 MP=1,NP
      DO 30 MQ=1,MP
      MRSPQ=MRSPQ+1
      VAL=XRSPQ(MRSPQ)
      IF(VAL.EQ.ZERO) GO TO 30
      DO 20 MI=1,MORB
      XRSIQ(MI,MQ)=XRSIQ(MI,MQ)+VAL*V(MP,MI)
      XRSIQ(MI,MP)=XRSIQ(MI,MP)+VAL*V(MQ,MI)
   20 CONTINUE
   30 CONTINUE
   40 CONTINUE
C
      MRSIJ=MRSIJ0
      DO 90 MI=1,MORB
      DO 50 MJ=1,MI
   50 XIJ(MJ)=ZERO
      DO 70 MQ=1,NQ
      VAL=XRSIQ(MI,MQ)
      IF(VAL.EQ.ZERO) GO TO 70
      DO 60 MJ=1,MI
   60 XIJ(MJ)=XIJ(MJ)+VAL*V(MQ,MJ)
   70 CONTINUE
      DO 80 MJ=1,MI
      MRSIJ=MRSIJ+1
   80 XRSIJ(MRSIJ)=XIJ(MJ)
   90 CONTINUE
C
      MRSIJ0=MRSIJ0+NRS
  100 CONTINUE
C
C     ----- TRANSPOSE (R,S//I,J) -----
C
      NIJ=(MORB*(MORB+1))/2
      MRSIJ=0
      DO 120 MRS=1,NRS
      MIJRS=0
      MAX=MRS
      IF(MAX.GT.NIJ) MAX=NIJ
      DO 110 MIJ=1,MAX
      DUM=XRSIJ(MRSIJ+MIJ)
      XIJRS(MRSIJ+MIJ)=XRSIJ(MIJRS+MRS)
      XIJRS(MIJRS+MRS)=DUM
  110 MIJRS=MIJRS+NRS
  120 MRSIJ=MRSIJ+NRS
C
C     ----- TRANSFORM -R- AND -S- -----
C
      NR=NUM
      NS=NUM
      MIJKL=0
      MIJRS=0
      MIJ=0
      MLEFT=NRS-NIJ
C
      DO 300 MI=1,MORB
      DO 300 MJ=1,MI
      MIJ=MIJ+1
C
      DO 210 MS=1,NS
      DO 210 MK=MI,MORB
  210 XIJKS(MK,MS)=ZERO
      DO 240 MR=1,NR
      DO 230 MS=1,MR
      MIJRS=MIJRS+1
      VAL=XIJRS(MIJRS)
      IF(VAL.EQ.ZERO) GO TO 230
      DO 220 MK=MI,MORB
      XIJKS(MK,MS)=XIJKS(MK,MS)+VAL*V(MR,MK)
      XIJKS(MK,MR)=XIJKS(MK,MR)+VAL*V(MS,MK)
  220 CONTINUE
  230 CONTINUE
  240 CONTINUE
C
      LMIN=MJ
      LMAX=MI
      DO 290 MK=MI,MORB
      DO 250 ML=LMIN,LMAX
  250 XKL(ML)=ZERO
      DO 270 MS=1,NS
      VAL=XIJKS(MK,MS)
      IF(VAL.EQ.ZERO) GO TO 270
      DO 260 ML=LMIN,LMAX
  260 XKL(ML)=XKL(ML)+VAL*V(MS,ML)
  270 CONTINUE
      DO 280 ML=LMIN,LMAX
      MIJKL=MIJKL+1
  280 XIJKL(MIJKL)=XKL(ML)
      LMIN=1
      LMAX=MK+1
  290 CONTINUE
      MIJKL=MIJKL+MIJ+MLEFT
C
  300 CONTINUE
      RETURN
      END
      SUBROUTINE TRFPRT(NORB,NIJ,IA,TX,JFLAG)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON/IOFILE/IR,IW
      COMMON/TR2PRT/VALINT(3),JC,JCINT(15)
      DIMENSION TX(NIJ,1),IA(1)
C
 1000 FORMAT(' AO 2E-INTEGRALS : ')
 1010 FORMAT(' TRANSFORMED 2E-INTEGRALS : ')
 1020 FORMAT(' TOTAL NUMBER OF AO 2E-INTEGRALS  = ',I10)
 1030 FORMAT(' TOTAL NUMBER OF TRANSFORMED 2E-INTEGRALS  = ',I10)
C
      IF(JFLAG.EQ.0) WRITE(IW,1000)
      IF(JFLAG.EQ.1) WRITE(IW,1010)
C
      JC=0
      NX=0
      DO 140 I=1,NORB
      DO 130 J=1,I
      DO 120 K=1,I
      LMAX=K
      IF(K.EQ.I) LMAX=J
      DO 110 L=1,LMAX
      NX=NX+1
      IJ=IA(I)+J
      KL=IA(K)+L
      VAL=TX(IJ,KL)
      CALL TRFOUT(I,J,K,L,NX,VAL,0)
  110 CONTINUE
  120 CONTINUE
  130 CONTINUE
  140 CONTINUE
      IF(JFLAG.EQ.0) WRITE(IW,1020) NX
      IF(JFLAG.EQ.1) WRITE(IW,1030) NX
      RETURN
      END
      SUBROUTINE TRFOUT(I,J,K,L,N,VAL,NCALL)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON/IOFILE/IR,IW
      COMMON/TR2PRT/V(3),JC,J1(3),J2(3),J3(3),J4(3),N1(3)
      IF(NCALL.NE.0) GO TO 10
      JC=JC+1
      J1(JC)=I
      J2(JC)=J
      J3(JC)=K
      J4(JC)=L
      N1(JC)=N
      V(JC)=VAL
      IF(JC.LT.3) RETURN
   10 CONTINUE
      IF(JC.EQ.0) RETURN
      WRITE(IW,9999) (J1(M),J2(M),J3(M),J4(M),N1(M),V(M),M=1,JC)
      JC=0
 9999 FORMAT(3(4I3,I5,E20.12))
      RETURN
      END
      FUNCTION NFLOP(N)
      NFLP=0
      MORB=N
      NRS=(N*(N+1))/2
      NP=N
      NQ=N
      DO 100 MRS=1,NRS
      DO 40 MP=1,NP
      DO 30 MQ=1,MP
      DO 20 MI=1,MORB
      NFLP=NFLP+(2+2)
   20 CONTINUE
   30 CONTINUE
   40 CONTINUE
      DO 90 MI=1,MORB
      DO 70 MQ=1,NQ
      DO 60 MJ=1,MI
      NFLP=NFLP+(1+1)
   60 CONTINUE
   70 CONTINUE
   90 CONTINUE
  100 CONTINUE
C
      NR=N
      NS=N
      DO 300 MI=1,MORB
      DO 300 MJ=1,MI
      DO 240 MR=1,NR
      DO 230 MS=1,MR
      DO 220 MK=MI,MORB
      NFLP=NFLP+(2+2)
  220 CONTINUE
  230 CONTINUE
  240 CONTINUE
      LMIN=MJ
      LMAX=MI
      DO 290 MK=MI,MORB
      DO 270 MS=1,NS
      DO 260 L=LMIN,LMAX
      NFLP=NFLP+(1+1)
  260 CONTINUE
  270 CONTINUE
      LMIN=1
      LMAX=MK+1
  290 CONTINUE
  300 CONTINUE
C
      NFLOP=NFLP
      RETURN
      END


text:
Module IO
implicit none

private

character(len=9), parameter :: Euler_file    = "Euler.out"
integer,          parameter ::    Euler_unit = 10
integer,          parameter :: Particle_unit = 11

public :: create_Eulerian_file, create_Particle_file
public :: save_Eulerian, save_particles

contains

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine create_Eulerian_file
use GRID, only : Z_r
implicit none

! Create data out file
open (unit=Euler_unit, file = Euler_file, status = 'replace')
write(Euler_unit, 101) 'Type','Timestep','Day', 'Hour', Z_r
close(Euler_unit)

101 format(4(A10,1x), 200(2x,F12.5))
end subroutine

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine save_Eulerian
use state_variables, only : Labelout, Varout, Nout
use Time_setting,    only : it, current_day, current_hour
implicit none
integer :: i

! save data into the Eulerian file
open (unit=Euler_unit, file = Euler_file, status = 'old', action='write', position='append')
do i = 1, Nout
   write(Euler_unit, 103) trim(Labelout(i)), it, current_day, current_hour, Varout(i,:)
enddo
close(Euler_unit)

103 format(A5, 1x, I7, 1x, 2(I0, 1x), 200(1pe12.3, 1x))
end subroutine

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine create_Particle_file(fname)
implicit none
character(len=*), intent(in) :: fname  !Need to create multiple files for storing particles to avoid huge files

open (unit=Particle_unit, file = fname, status = 'replace')
write(Particle_unit, 102) 'Timestep','Day','Hour','ID','Grid', &
                          'Z', 'C', 'N', 'Chl', 'NO3', 'PAR',  &
                          'Temp', 'Num', 'Topt', 'CDiv', 'Iopt'
close(Particle_unit)
return

102 format(16(A8, 2x))
END subroutine create_Particle_file
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine save_particles(fname)
use state_variables, only : p_PHY
use Time_setting,    only : it, current_day, current_hour
implicit none
character(len=*), intent(in) :: fname  !The file for saving particles
integer :: i

! save data into the particle file
open (unit=Particle_unit, file=fname, status='old', action='write', position='append')
do i = 1, size(p_PHY)
   write(Particle_unit, 104) it,                             &
                  current_day, current_hour, p_PHY(i)%ID, &
                  p_PHY(i)%iz, p_PHY(i)%rz,  p_PHY(i)%C,  &
                  p_PHY(i)%N,  p_PHY(i)%Chl, p_PHY(i)%NO3,&
                  p_PHY(i)%PAR,p_PHY(i)%temp,p_PHY(i)%num, &
                  p_PHY(i)%Topt, p_PHY(i)%CDiv, p_PHY(i)%lnIopt
enddo
close(Particle_unit)

104 format(I0, 1x, 2(I8, 1x), I0, 1x, I3, 1x, F12.3, 1x, 10(1pe12.3, 1x))
end subroutine
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
END MODULE IO

max_stars_repo_name: BingzhangChen/IBM

max_issues_repo_name: BingzhangChen/IBM

max_forks_repo_name: BingzhangChen/IBM


! ##################################################################################################################################
! Begin MIT license text.
! _______________________________________________________________________________________________________

! Copyright 2019 Dr William R Case, Jr (dbcase29@gmail.com)

! Permission is hereby granted, free of charge, to any person obtaining a copy of this software and
! associated documentation files (the "Software"), to deal in the Software without restriction, including
! without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to
! the following conditions:

! The above copyright notice and this permission notice shall be included in all copies or substantial
! portions of the Software and documentation.

! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
! OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
! THE SOFTWARE.
! _______________________________________________________________________________________________________

! End MIT license text.

      SUBROUTINE GET_GRID_NUM_COMPS ( GRID_NUM, NUM_COMPS, CALLING_SUBR )

! Gets the number of components for a "grid" number from array GRID by testing GRID(I,6)
! If GRID(I,6) is 6 then this is a physical grid with 6 comps of displ and if GRID(I,6) is 1 then this is an SPOINT.

      USE PENTIUM_II_KIND, ONLY       :  BYTE, LONG, DOUBLE
      USE IOUNT1, ONLY                :  WRT_LOG, ERR, F04, F06
      USE SCONTR, ONLY                :  BLNK_SUB_NAM, FATAL_ERR, NGRID
      USE TIMDAT, ONLY                :  TSEC
      USE MODEL_STUF, ONLY            :  GRID
      USE SUBR_BEGEND_LEVELS, ONLY    :  GET_GRID_NUM_COMPS_BEGEND

      USE GET_GRID_NUM_COMPS_USE_IFs

      IMPLICIT NONE

      CHARACTER(LEN=LEN(BLNK_SUB_NAM)):: SUBR_NAME = 'GET_GRID_NUM_COMPS'
      CHARACTER(LEN=*), INTENT(IN)    :: CALLING_SUBR      ! Subr that called this one

      INTEGER(LONG), INTENT(IN)       :: GRID_NUM          ! A grid number (calling subr checked that it is an actual grid)
      INTEGER(LONG), INTENT(OUT)      :: NUM_COMPS         ! 6 if GRID_NUM is an physical grid, 1 if an SPOINT
      INTEGER(LONG)                   :: I                 ! DO loop index
      INTEGER(LONG), PARAMETER        :: SUBR_BEGEND = GET_GRID_NUM_COMPS_BEGEND

! **********************************************************************************************************************************
      IF (WRT_LOG >= SUBR_BEGEND) THEN
         CALL OURTIM
         WRITE(F04,9001) SUBR_NAME,TSEC
 9001    FORMAT(1X,A,' BEGN ',F10.3)
      ENDIF

! **********************************************************************************************************************************
! Initialize outputs

      NUM_COMPS = 0

! Calc NUM_COMPS

      DO I=1,NGRID
         IF (GRID(I,1) == GRID_NUM) THEN
            NUM_COMPS = GRID(I,6)
            EXIT
         ENDIF
      ENDDO

! Error if NUM_COMPS not 1 or 6

      IF ((NUM_COMPS /= 1) .AND. (NUM_COMPS /= 6)) THEN
         FATAL_ERR = FATAL_ERR + 1
         WRITE(ERR,950) SUBR_NAME, CALLING_SUBR, GRID_NUM, NUM_COMPS
         WRITE(F06,950) SUBR_NAME, CALLING_SUBR, GRID_NUM, NUM_COMPS
         CALL OUTA_HERE ( 'Y' )
      ENDIF

! **********************************************************************************************************************************
      IF (WRT_LOG >= SUBR_BEGEND) THEN
         CALL OURTIM
         WRITE(F04,9002) SUBR_NAME,TSEC
 9002    FORMAT(1X,A,' END  ',F10.3)
      ENDIF

      RETURN

! **********************************************************************************************************************************
  950 FORMAT(' *ERROR   945: PROGRAMMING ERROR IN SUBROUTINE ',A,' CALLED BY SUBROUTINE ',A                                        &
                    ,/,14X,' FOR GRID ',I8,' THE NUMBER OF DISPL COMPONENTS IS ',I8,' BUT MUST BE EITHER 1 (SPOINT) OR 6 (GRID)')

! **********************************************************************************************************************************

      END SUBROUTINE GET_GRID_NUM_COMPS