text:
      subroutine create_graph(imsgin,xmsgin,cmsgin,msgtype,nwds,
     *                ierror_return)
C
C
C#######################################################################
C
C      PURPOSE - Create the adjacency graph for a mesh.
C      If node option is selected, the graph of node adjacency
C      is created, if dual hption is selected, the graph of
C      element adjacency (dual graph) is created. See METIS
C      documentation for description of graph format.
C
C     Command syntax
C
C    create_graph / method / graph_type / offset_vector / adjancy_vector / action
C
C    method = metis - use METIS functions to create graph
C           = lagrit - use LAGriT function to create graph (not implemented)
C
C    graph_type = node - node adjacency graph is created
C               = dual - dual (element) adjacency graph is created
C
C    offset_vector = name of integer vector which holds offset pointers into adjancy vector
C                    default name = nxadj
C
C    adjancy_vector = name of integer vector which holds adjancy information
C                     default name = nadjncy
C
C    action = create - create graph arrays
C           = delete - delete graph arrays
C
C    create_graph / {metis,lagrit]  / [node,dual] / nxadj / nadjncy / create,delete
CC
C    lagrit option not implemented. It is a place holder for future development.
C    The graph for a hybrid dual graph can be created from jtet array.
C
C   NOTE: create_graph options will not work on a hybrid mesh
C         and the only element types supported are tri, tet, quad, hex
C
C      INPUT ARGUMENTS -
C
C         imsgin()  - Integer array of command input tokens
C         xmsgin()  - Real array of command input tokens
C         cmsgin()  - Character array of command input tokens
C         msgtype() - Integer array of command input token types
C         nwds      - Number of command input tokens
C
C      OUTPUT ARGUMENTS -
C
C         ierror_return - Error Return Code (==0 ==> OK, <>0 ==> Error)
C
C      CHANGE HISTORY -
C
C      $Log: create_graph.f,v $
C      Revision 2.00  2007/11/05 19:45:51  spchu
C      Import to CVS
C
CPVCS
CPVCS       Rev 1.0   01 Nov 2002 13:04:32   gable
CPVCS    Initial revision.
C
C#######################################################################
C
      implicit none
C
C#######################################################################
C
C  Subroutine argument variables
C
      integer nwds, imsgin(nwds), msgtype(nwds)
      REAL*8 xmsgin(nwds)
      character*(*) cmsgin(nwds)
      integer ierror_return
C
C#######################################################################
C
      character*132 logmess
      character*256 cmdmess
      character*12   char_number, char_number2
      integer ierr, len, itype, iflag, indx
      integer ifactor
      integer if_node_1_dual_2, n_vertices_graph,numflag
      integer if_create_1_delete_2
      integer isize_nxadj,     isize_nadjncy,     irank_nadjncy
      integer isize_nxadj_set, isize_nadjncy_set, irank_nadjncy_set
C
C#######################################################################
C
      integer icharlnf
C
      character*32 isubname
C
C     Character variables for adjacency arrays
C
      character*32 ch_nadjncy,ch_nxadj,ch_isize_nxadj
      character*32 ch_isize_nadjncy, ch_irank_nadjncy
C
C     Metis  Element Types triangle=1 tetrahedra=2 hexahedra=3 quadrilateral=4
C     LaGriT Element Types triangle=3 tetrahedra=5 hexahedra=8 quadrilateral=4
C
C     Make a mapping between LaGriT element type and Metis element type.
C
      integer metis_element_type(8)
C
C#######################################################################
C
C     Mesh Object Arrays and Variables
C
      character*32 cmo_name
C
      integer nnode, nelem
C
      pointer (ipitet, itet)
      integer itet(10000000)
      pointer (ipitetoff, itetoff)
      integer itetoff(10000000)
      pointer (ipnxadj, nxadj)
      integer nxadj(10000000)
      pointer (ipnadjncy, nadjncy)
      integer nadjncy(10000000)
      pointer (ipitettyp, itettyp)
      integer itettyp(10000000)
C
C#######################################################################
C
C     Set some defaults and constants
C
C     Mapping between LAGriT element types and element types.
C
      data metis_element_type / 0, 0, 1, 4, 2, 0, 0, 3 /
C
C     METIS flag to indicate Fortran-style arrays start at 1
C
      data numflag / 1 /
C
C#######################################################################
C
C     Begin create_graph
C
C#######################################################################
C
      isubname = 'create_graph'
C
C     Get name of current MO
C
      call cmo_get_name(cmo_name,ierr)
      if(ierr .ne. 0)then
         logmess = 'CREATE_GRAPH ERROR: No active MO found'
         call writloga('default',0,logmess,0,ierr)
         ierror_return = -1
         go to 9999
      endif
C
C#######################################################################
C
C     Parse and check the second argument
C
      indx = 2
      if(cmsgin(indx)(1:icharlnf(cmsgin(indx))) .eq. 'metis')then

      else
        logmess='CREATE_GRAPH ERROR: 2th argument requires /metis/'
        call writloga('default',0,logmess,0,ierr)
        ierror_return = -2
        go to 9999
      endif
C
C#######################################################################
C     Parse and check the third argument
C
      indx = 3
      if(cmsgin(indx)(1:icharlnf(cmsgin(indx))) .eq. 'node')then
        if_node_1_dual_2 = 1
      elseif(cmsgin(indx)(1:icharlnf(cmsgin(indx))) .eq. 'dual')then
        if_node_1_dual_2 = 2
      else
        logmess='CREATE_GRAPH ERROR: '
        call writloga('default',0,logmess,0,ierr)
        logmess=' ERROR: 3rd argument options /node/ or /dual/'
        call writloga('default',0,logmess,0,ierr)
        ierror_return = -2
        go to 9999
      endif
C
C#######################################################################
C
C    Parse and check the fourth and fifth argument
C
C     Create the appropriate attributes names for output variables.
C
C***********************************************************************
C     create_graph needs nxadj, nadjncy
C
      if(if_node_1_dual_2 .eq. 1)then
          ch_isize_nxadj   = 'isize_nxadj_n'
	  ch_irank_nadjncy = 'irank_nadjncy_n'
	  ch_isize_nadjncy = 'nnodes'
      else
         ch_isize_nxadj   = 'isize_nxadj_e'
	 ch_irank_nadjncy = 'irank_nadjncy_e'
	 ch_isize_nadjncy = 'nelements'
      endif
C
      indx = 4
      if((nwds .ge. indx) .and. (msgtype(indx) .ne. 3))then
         logmess = 'CREATE_GRAPH ERROR: Invalid 4th argument'
         call writloga('default',0,logmess,0,ierr)
         logmess = 'CREATE_GRAPH ERROR: Requires character variable'
         call writloga('default',0,logmess,0,ierr)
         ierror_return = 3
         go to 9999
      endif
      if((nwds .lt. indx) .or.
     1   (cmsgin(indx)(1:icharlnf(cmsgin(indx))) .eq. '-def-')) then
         if(if_node_1_dual_2 .eq. 1)then
           ch_nxadj = 'nxadj_n'
	 else
	   ch_nxadj = 'nxadj_e'
	 endif
      else
         ch_nxadj = cmsgin(indx)(1:icharlnf(cmsgin(indx)))
      endif
C
      indx = 5
       if((nwds .ge. indx) .and. (msgtype(indx) .ne. 3))then
         logmess = 'CREATE_GRAPH ERROR: Invalid 5th argument'
         call writloga('default',0,logmess,0,ierr)
         logmess = 'CREATE_GRAPH ERROR: Requires character variable'
         call writloga('default',0,logmess,0,ierr)
         ierror_return = 3
         go to 9999
      endif
      if((nwds .lt. indx) .or.
     1   (cmsgin(indx)(1:icharlnf(cmsgin(indx))) .eq. '-def-')) then
         if(if_node_1_dual_2 .eq. 1)then
           ch_nadjncy = 'nadjncy_n'
         else
           ch_nadjncy = 'nadjncy_e'
         endif
      else
         ch_nadjncy = cmsgin(indx)(1:icharlnf(cmsgin(indx)))
      endif
C
C#######################################################################
C
C    Parse and check the sixth argument
C
      indx = 6
      if((nwds .ge. indx) .and. (msgtype(indx) .ne. 3))then
         logmess = 'CREATE_GRAPH ERROR: Invalid 6th argument'
         call writloga('default',0,logmess,0,ierr)
         logmess = 'CREATE_GRAPH ERROR: options are create or delete'
         call writloga('default',0,logmess,0,ierr)
         ierror_return = 3
         go to 9999
      endif
      if((nwds .lt. indx) .or.
     1   (cmsgin(indx)(1:icharlnf(cmsgin(indx))) .eq. '-def-')) then
	   if_create_1_delete_2 = 1
      elseif(cmsgin(indx)(1:icharlnf(cmsgin(indx))) .eq. 'create')then
	   if_create_1_delete_2 = 1
      elseif(cmsgin(indx)(1:icharlnf(cmsgin(indx))) .eq. 'delete')then
	   if_create_1_delete_2 = 2
      endif
C
C
      if(if_create_1_delete_2 .eq. 2)then
       cmdmess = 'cmo/delatt/ /'//ch_nxadj//' ; finish'
       call dotask(cmdmess,ierr)
       cmdmess = 'cmo/delatt/ /'//ch_nadjncy//' ; finish'
       call dotask(cmdmess,ierr)
       cmdmess =
     1   'cmo/delatt/'//cmo_name//'/'//ch_isize_nxadj//' ; finish'
       call dotask(cmdmess,ierr)
       cmdmess =
     1   'cmo/delatt/'//cmo_name//'/'//ch_irank_nadjncy//' ; finish'
       call dotask(cmdmess,ierr)
C
C      Exit
C
       go to 9999
       endif
C
C#######################################################################
C
C     Get MO pointers
C
      call cmo_get_info('nnodes',cmo_name,nnode,len,itype,ierr)
      call cmo_get_info('nelements',cmo_name,nelem,len,itype,ierr)
      call cmo_get_info('itettyp',cmo_name,ipitettyp,len,itype,ierr)
      call cmo_get_info('itetoff',cmo_name,ipitetoff,len,itype,ierr)
      call cmo_get_info('itet',cmo_name,ipitet,len,itype,ierr)
C
C#######################################################################
C
C    Make some checks that MO is of a type METIS can handle.
C
C     Check that it is an element type Metis supports.
C
      if(metis_element_type(itettyp(1)) .eq. 0)then
        logmess = 'CREATE_GRAPH ERROR: Element type not supported'
        call writloga('default',0,logmess,0,ierr)
        return
      endif
C
C     Check that mesh is not hybrid by checking that all values of itettyp
C     are equal.
C
      call check_int_vector_equal(nelem,itettyp,iflag)
      if(iflag .ne. 0)then
C
C       If iflag is not zero then it must be a hybrid mesh.
C       METIS will not handle a hybrid mesh.
C
        logmess = 'CREATE_GRAPH ERROR: Hybrid mesh not supported'
        call writloga('default',0,logmess,0,ierr)
        return
      endif

C
C#######################################################################
C
C     Check for existance of and allocate if necessary the nxadj, nadjncy arrays.
C
C     Figure out the size various arrays will require.
C
        if(if_node_1_dual_2 .eq. 1)then
          isize_nxadj_set = nnode + 1
	  write(char_number2,'(i12)')isize_nxadj_set
       elseif(if_node_1_dual_2 .eq. 2)then
          isize_nxadj_set = nelem + 1
	  write(char_number2,'(i12)')isize_nxadj_set
       endif
C
C     Figure out the size of the nadjncy array based on
C     element type (tri, quad, tet, hex)
C     graph type (node, dual)
C
      ifactor = 2
C
      if    (itettyp(1) .eq. 3 .and. if_node_1_dual_2 .eq. 1)then
         isize_nadjncy_set = ifactor*6*nnode  ! triangles, node mesh
	 irank_nadjncy_set = ifactor*6
	 write(char_number,'(f6.2)')6.0
      elseif(itettyp(1) .eq. 3 .and. if_node_1_dual_2 .eq. 2)then
         isize_nadjncy_set = 3*nelem  ! triangles, dual mesh
	 irank_nadjncy_set = 3
	 write(char_number,'(f6.2)')3.0

      elseif(itettyp(1) .eq. 4 .and. if_node_1_dual_2 .eq. 1)then
         isize_nadjncy_set = ifactor*4*nnode  ! quads, node mesh
	 irank_nadjncy_set = ifactor*4
	 write(char_number,'(f6.2)')4.0
      elseif(itettyp(1) .eq. 4 .and. if_node_1_dual_2 .eq. 2)then
         isize_nadjncy_set = 4*nelem  ! quads, dual mesh
	 irank_nadjncy_set = 4
	 write(char_number,'(f6.2)')4.0

      elseif(itettyp(1) .eq. 5 .and. if_node_1_dual_2 .eq. 1)then
         isize_nadjncy_set = ifactor*15*nnode  ! tets, node mesh
	 irank_nadjncy_set = ifactor*15
	 write(char_number,'(f6.2)')15.0
      elseif(itettyp(1) .eq. 5 .and. if_node_1_dual_2 .eq. 2)then
         isize_nadjncy_set = 4*nelem  ! tets, dual mesh
	 irank_nadjncy_set = 4
	 write(char_number,'(f6.2)')4.0

      elseif(itettyp(1) .eq. 8 .and. if_node_1_dual_2 .eq. 1)then
         isize_nadjncy_set = ifactor*6*nnode  ! hex, node mesh
	 irank_nadjncy_set = ifactor*6
	 write(char_number,'(f6.2)')12.0
      elseif(itettyp(1) .eq. 8 .and. if_node_1_dual_2 .eq. 2)then
         isize_nadjncy_set = 6*nelem  ! hex, dual mesh
	 irank_nadjncy_set = 6
	 write(char_number,'(f6.2)')12.0
      endif
C
C     Right now the 'node mesh' case is just a guess. There are clearly
C     cases where the sizes set above are too small, such as a case where
C     you have N triangles incident on a single point and N is a large number.
C
C     Check to see if the adjacency array exists,
C     if it doesn't, create it.
C
        call cmo_get_info
     1    (ch_isize_nxadj,cmo_name,isize_nxadj,len,itype,ierr)
C
C     NOTE: Since there are dotask calls below, one can no longer use
C           the subroutine argument variables. They are changed by the
C           dotask calls.
C
C     If the isize_nxadj is different than isize_nxadj_set then deallocate
C     and reallocate the nxadj array.
C
      if((isize_nxadj .ne. isize_nxadj_set) .and. (ierr .eq. 0))then
C
C        If you are here, it means isize_nxadj already exists but is the
C        wrong size. Just deallocate arrays and reallocate arrays of the
C        correct size below.
C
C       NOTE: should be able to delete with two calls to DELATT
C             the rank and size parameters, but that appears to
C             not be working. Use lower case delatt
C
       cmdmess = 'cmo/delatt/ /'//ch_nxadj//' ; finish'
       call dotask(cmdmess,ierr)
       cmdmess = 'cmo/delatt/ /'//ch_nadjncy//' ; finish'
       call dotask(cmdmess,ierr)
       cmdmess =
     1   'cmo/delatt/'//cmo_name//'/'//ch_isize_nxadj//' ; finish'
       call dotask(cmdmess,ierr)
       cmdmess =
     1   'cmo/delatt/'//cmo_name//'/'//ch_irank_nadjncy//' ; finish'
       call dotask(cmdmess,ierr)
C
C        Set ierr to non-zero value so allocation will occur below.
C
         ierr = -1
      endif
C
        call cmo_get_info
     1    (ch_irank_nadjncy,cmo_name,irank_nadjncy,len,itype,ierr)
C
C     If the irank_nadjncy is different than irank_nadjncy_set then deallocate
C     and reallocate the nxadj array.
C
      if((irank_nadjncy .ne. irank_nadjncy_set) .and. (ierr .eq. 0))then
C
C        If you are here, it means irank_nxadj already exists but is the
C        wrong size. Just deallocate arrays and reallocate arrays of the
C        correct size below.
C
       cmdmess = 'cmo/delatt/ /'//ch_nxadj//' ; finish'
       call dotask(cmdmess,ierr)
       cmdmess = 'cmo/delatt/ /'//ch_nadjncy//' ; finish'
       call dotask(cmdmess,ierr)
       cmdmess =
     1   'cmo/delatt/'//cmo_name//'/'//ch_isize_nxadj//' ; finish'
       call dotask(cmdmess,ierr)
       cmdmess =
     1   'cmo/delatt/'//cmo_name//'/'//ch_irank_nadjncy//' ; finish'
       call dotask(cmdmess,ierr)
C
C        Set ierr to non-zero value so allocation will occur below.
C
         ierr = -1
      endif

      if (ierr.ne.0) then
         cmdmess ='cmo/addatt//'//ch_isize_nxadj//'/INT' //
     &      '/scalar/scalar/constant/temporary/l/'// char_number2 //
     &      ' ; finish'
         call dotask(cmdmess,ierr)
         call cmo_get_info
     1    (ch_isize_nxadj,cmo_name,isize_nxadj,len,itype,ierr)
       endif

      call cmo_get_info(ch_nxadj,cmo_name,ipnxadj,len,itype,ierr)
      if(ierr .ne. 0)then
          cmdmess='cmo/addatt/' //
     &      '/' //
     &      ch_nxadj(1:icharlnf(ch_nxadj)) //
     &      '/VINT' //
     &      '/scalar/'//ch_isize_nxadj//'//temporary/l/0.0' //
     &      ' ; finish'
       ierr = 0
       call dotask(cmdmess,ierr)
       call cmo_get_info(ch_nxadj,cmo_name,ipnxadj,len,itype,ierr)
      endif

      call cmo_get_info
     1    (ch_irank_nadjncy,cmo_name,irank_nadjncy,len,itype,ierr)
      if (ierr.ne.0) then
         cmdmess ='cmo/addatt//'//ch_irank_nadjncy//'/INT' //
     &      '/scalar/scalar/constant/temporary/l/'// char_number //
     &      ' ; finish'
         call dotask(cmdmess,ierr)
         call cmo_get_info
     1    (ch_irank_nadjncy,cmo_name,irank_nadjncy,len,itype,ierr)
       endif

      call cmo_get_info(ch_nadjncy,cmo_name,ipnadjncy,len,itype,ierr)
      if (ierr.ne.0) then
         cmdmess='cmo/addatt//' //
     &      ch_nadjncy(1:icharlnf(ch_nadjncy)) //
     &      '/VINT/'//ch_irank_nadjncy//'/'//ch_isize_nadjncy//
     &      '//temporary/l/0.0'//
     &      ' ; finish'

       call dotask(cmdmess,ierr)
       call cmo_get_info(ch_nadjncy,cmo_name,ipnadjncy,len,itype,ierr)
      endif
C
C     *****************************************************************
C
        if(if_node_1_dual_2 .eq. 1)then
C
C METIS_MeshToNodal
	  call metis_meshtonodal(nelem,nnode,itet,
     1         metis_element_type(itettyp(1)),numflag,nxadj,nadjncy)
          n_vertices_graph = nnode
C
	elseif(if_node_1_dual_2 .eq. 2)then
C METIS_MeshToDual
	  call metis_meshtodual(nelem,nnode,itet,
     1         metis_element_type(itettyp(1)),numflag,nxadj,nadjncy)
          n_vertices_graph = nelem
	endif
C     *****************************************************************
C
C     Release any temporary memory associated with partition isubname
C
 9999 continue
      call mmrelprt(isubname,ierr)

      return
      end

ext: f

max_forks_repo_name: daniellivingston/LaGriT


SUBROUTINE linbcg_sp(b,x,itol,tol,itmax,iter,err)
  USE mo_kind, only: sp, i4
  USE mo_nrutil, ONLY : assert_eq, nrerror
  use mo_nr, ONLY : atimes, asolve, snrm
  IMPLICIT NONE
  REAL(SP), DIMENSION(:), INTENT(IN) :: b
  REAL(SP), DIMENSION(:), INTENT(INOUT) :: x
  INTEGER(I4), INTENT(IN) :: itol,itmax
  REAL(SP), INTENT(IN) :: tol
  INTEGER(I4), INTENT(OUT) :: iter
  REAL(SP), INTENT(OUT) :: err
  REAL(SP), PARAMETER :: EPS=1.0e-9_sp
  INTEGER(I4) :: n
  REAL(SP) :: ak,akden,bk,bkden,bknum,bnrm,dxnrm,xnrm,zm1nrm,znrm
  REAL(SP), DIMENSION(size(b)) :: p,pp,r,rr,z,zz
  n=assert_eq(size(b),size(x),'linbcg_sp')
  iter=0
  call atimes(x,r,0)
  r=b-r
  rr=r
  !	call atimes(r,rr,0)
  select case(itol)
  case(1)
     bnrm=snrm(b,itol)
     call asolve(r,z,0)
  case(2)
     call asolve(b,z,0)
     bnrm=snrm(z,itol)
     call asolve(r,z,0)
  case(3:4)
     call asolve(b,z,0)
     bnrm=snrm(z,itol)
     call asolve(r,z,0)
     znrm=snrm(z,itol)
  case default
     call nrerror('illegal itol in linbcg_sp')
  end select
  do
     if (iter > itmax) exit
     iter=iter+1
     call asolve(rr,zz,1)
     bknum=dot_product(z,rr)
     if (iter == 1) then
        p=z
        pp=zz
     else
        bk=bknum/bkden
        p=bk*p+z
        pp=bk*pp+zz
     end if
     bkden=bknum
     call atimes(p,z,0)
     akden=dot_product(z,pp)
     ak=bknum/akden
     call atimes(pp,zz,1)
     x=x+ak*p
     r=r-ak*z
     rr=rr-ak*zz
     call asolve(r,z,0)
     select case(itol)
     case(1)
        err=snrm(r,itol)/bnrm
     case(2)
        err=snrm(z,itol)/bnrm
     case(3:4)
        zm1nrm=znrm
        znrm=snrm(z,itol)
        if (abs(zm1nrm-znrm) > EPS*znrm) then
           dxnrm=abs(ak)*snrm(p,itol)
           err=znrm/abs(zm1nrm-znrm)*dxnrm
        else
           err=znrm/bnrm
           cycle
        end if
        xnrm=snrm(x,itol)
        if (err <= 0.5_sp*xnrm) then
           err=err/xnrm
        else
           err=znrm/bnrm
           cycle
        end if
     end select
     !write (*,*) ' iter=',iter,' err=',err
     if (err <= tol) exit
  end do
END SUBROUTINE linbcg_sp

SUBROUTINE linbcg_dp(b,x,itol,tol,itmax,iter,err)
  USE mo_kind, only: dp, i4
  USE mo_nrutil, ONLY : assert_eq, nrerror
  use mo_nr, ONLY : atimes, asolve, snrm
  IMPLICIT NONE
  REAL(DP), DIMENSION(:), INTENT(IN) :: b
  REAL(DP), DIMENSION(:), INTENT(INOUT) :: x
  INTEGER(I4), INTENT(IN) :: itol,itmax
  REAL(DP), INTENT(IN) :: tol
  INTEGER(I4), INTENT(OUT) :: iter
  REAL(DP), INTENT(OUT) :: err
  REAL(DP), PARAMETER :: EPS=1.0e-14_dp
  INTEGER(I4) :: n
  REAL(DP) :: ak,akden,bk,bkden,bknum,bnrm,dxnrm,xnrm,zm1nrm,znrm
  REAL(DP), DIMENSION(size(b)) :: p,pp,r,rr,z,zz
  n=assert_eq(size(b),size(x),'linbcg_dp')
  iter=0
  call atimes(x,r,0)
  r=b-r
  rr=r
  !	call atimes(r,rr,0)
  select case(itol)
  case(1)
     bnrm=snrm(b,itol)
     call asolve(r,z,0)
  case(2)
     call asolve(b,z,0)
     bnrm=snrm(z,itol)
     call asolve(r,z,0)
  case(3:4)
     call asolve(b,z,0)
     bnrm=snrm(z,itol)
     call asolve(r,z,0)
     znrm=snrm(z,itol)
  case default
     call nrerror('illegal itol in linbcg_dp')
  end select
  do
     if (iter > itmax) exit
     iter=iter+1
     call asolve(rr,zz,1)
     bknum=dot_product(z,rr)
     if (iter == 1) then
        p=z
        pp=zz
     else
        bk=bknum/bkden
        p=bk*p+z
        pp=bk*pp+zz
     end if
     bkden=bknum
     call atimes(p,z,0)
     akden=dot_product(z,pp)
     ak=bknum/akden
     call atimes(pp,zz,1)
     x=x+ak*p
     r=r-ak*z
     rr=rr-ak*zz
     call asolve(r,z,0)
     select case(itol)
     case(1)
        err=snrm(r,itol)/bnrm
     case(2)
        err=snrm(z,itol)/bnrm
     case(3:4)
        zm1nrm=znrm
        znrm=snrm(z,itol)
        if (abs(zm1nrm-znrm) > EPS*znrm) then
           dxnrm=abs(ak)*snrm(p,itol)
           err=znrm/abs(zm1nrm-znrm)*dxnrm
        else
           err=znrm/bnrm
           cycle
        end if
        xnrm=snrm(x,itol)
        if (err <= 0.5_dp*xnrm) then
           err=err/xnrm
        else
           err=znrm/bnrm
           cycle
        end if
     end select
     !write (*,*) ' iter=',iter,' err=',err
     if (err <= tol) exit
  end do
END SUBROUTINE linbcg_dp


text:
program t
  implicit none
! inquire stmt 'delim' specifier (apostrophe)
  character*10::a
  open (56, status='new', file='tmpfile', delim='apostrophe')
  inquire (56, delim=a)
  print *,a
  close (56,status='delete')
   endprogram t

max_stars_repo_name: maddenp/ppp

max_issues_repo_name: maddenp/ppp

max_forks_repo_name: maddenp/ppp


text:
C    @(#)evalfn.f	20.7 5/27/98
      real function evalfn (string, ndef, defstr, defval, defind,
     1                      error, mode)
      character string*(*), defstr(*)*(*)
      integer defind(*), error, mode
      real defval(*)
C
C     This function evaluates symbolic expressions such as
C
C     2.0 * B1.BS * (B1.PG + B1.QG) ** 2
C
c     MODE defines the source: 0 (default) = base case
c                              1           = alternate base case
c
      include 'ipfinc/parametr.inc'
      include 'ipfinc/prt.inc'

      character capital * 132, symbol * 12, fn_code * 4,
     1          word(30) * 30, newstr * 132, cfmt * 7, tempc * 30
      integer o1, s1, s2, type, itor
      real rtoi
      equivalence (itor, rtoi)
      complex zc

      error = 0
      newstr = capital (string)
C
C     Search for symbols
C
      lc = len (string)
      call xscan (newstr(1:lc), word(2), nwrd)
C
C     Append level 0 parenthesis
C
      word(1) = '('
      word(nwrd+2) = ')'
      nwrd = nwrd + 2
      do 140 iw = 2, nwrd-1
         i1 = index ('ABCDEFGHIJKLMNOPQRSTUVWXYZ',word(iw)(1:1))
         if (i1 .gt. 0) then
C
C             Determine whether symbol is ordinary or indexed.
C
           jc = lastch (word(iw))
           if (index (word(iw), '.') .ne. 0) then
              i2 = index (word(iw), '.')
              symbol = word(iw)(1:i2-1)
              fn_code = word(iw)(i2+1:jc)
           else
              symbol = word(iw)(1:jc)
              fn_code = ' '
           endif
C
C          Find symbol value
C
           lc = lastch (symbol)
           do j = 1, ndef
              if (defstr(j) .eq. symbol(1:lc)) then
                 if (fn_code .ne. ' ' .and. defind(j) .ne. 0) then
c
c                   real-to-integer conversion averted
c
                    rtoi = defval(j)
                    ind = itor
                    type = defind(j)
                    if (mode. eq. 0) then
                       value = getval (ind, type, fn_code, error)
                    else
                       value = xgetval (ind, type, fn_code, error)
                    endif
                    if (error .ne. 0) then
                       write (errbuf(1), 102) symbol(1:lc)
     1                    // '.' // fn_code
  102                  format (' indexed SYMBOL (', a,
     &                         ') is not installed.')
                       last = min0 (80, len(string))
                       write (errbuf(2), 104) string(1:last)
  104                  format (' Text (', a, ')')
                       call prterx ('W', 2)
                       write (word(iw), 110) 0.0
                       error = 1
                       go to 140
                    endif
                 else if (defind(j) .ne. 0) then
                    write (errbuf(1), 106) symbol(1:lc)
     1                                     // '.' // fn_code
  106               format (' SYMBOL (', a,
     &                      ') has "INDEX" missing in >LET command.')
                    last = min0 (80, len(string))
                    write (errbuf(2), 104) string(1:last)
                    call prterx ('W', 2)
                    write (word(iw), 110) 0.0
                    error = 1
                    go to 140
                 else if (fn_code .ne. ' ') then
                    write (errbuf(1), 108) symbol(1:lc)
  108               format (' indexed SYMBOL (', a,
     1                      ') is missing "TYPE" (e.g. .PL)')
                    last = min0 (80, len(string))
                    write (errbuf(2), 104) string(1:last)
                    call prterx ('W', 2)
                    write (word(iw), 110) 0.0
                    go to 140
                 else
                    value = defval(j)
                 endif
                 write (word(iw), 110) value
  110            format (e15.8)
                 go to 140
              endif
           enddo
C
C          Test for special functions
C
           if (symbol(1:lc) .eq. 'SIN' .or.
     &         symbol(1:lc) .eq. 'COS' .or.
     &         symbol(1:lc) .eq. 'TAN' .or.
     &         symbol(1:lc) .eq. 'ARCSIN' .or.
     &         symbol(1:lc) .eq. 'ARCCOS' .or.
     &         symbol(1:lc) .eq. 'ARCTAN' .or.
     &         symbol(1:lc) .eq. 'ABS') then
           else
              write (errbuf(1), 130) symbol(1:lc)
  130         format (' Symbol (', a, ') is not defined')
              last = min0 (80, len(string))
              write (errbuf(2), 104) string(1:last)
              call prterx ('W', 2)
              write (word(iw), 110) 0.0
              error = 1
           endif
        endif
  140 continue
      if (error .ne. 0) go to 900
C
C     Evaluate level of innermost parenthesis
C
  160 maxlev = 0
      s1 = 0
      s2 = 0
      level = 0
      do iw = 1, nwrd
         if (word(iw) .eq. '(') then
            level = level + 1
            if (level .ge. maxlev) then
               maxlev = level
               s1 = iw
            endif
         else if (word(iw) .eq. ')') then
            if (level .eq. maxlev) then
               s2 = iw
            endif
            level = level - 1
         endif
      enddo

      if ((s1 .gt. 0 .and. s2 .eq. 0) .or.
     &    (s1 .eq. 0 .and. s2 .gt. 0)) then
         write (errbuf(1), 171) string
  171    format (' Unbalanced parenthesis in string (', a, ')')
         call prterx ('W', 2)
         error = 1
         go to 900
      endif
      if (maxlev .gt. 0) then
C
C        Find highest operator within string
C
         o1 = 0
         do iw = s1+1, s2-1
            if (word(iw) .eq. '**' .or. word(iw) .eq. '^') o1 = iw
         enddo
         if (o1 .eq. 0) then
            iw = s1+1
            do while (iw .le. s2-1 .and. o1 .eq. 0)
               jc = lastch (word(iw))
               if (word(iw)(1:jc) .eq. '/') then
                  o1 = iw
               endif
               iw = iw + 1
            enddo
            if (o1 .eq. 0) then
               iw = s1+1
               do while (iw .le. s2-1 .and. o1 .eq. 0)
                  jc = lastch (word(iw))
                  if (word(iw)(1:jc) .eq. '*') then
                     o1 = iw
                  endif
                  iw = iw + 1
               enddo
               if (o1 .eq. 0) then
                  iw = s1+1
                  do while (iw .le. s2-1 .and. o1 .eq. 0)
                     jc = lastch (word(iw))
                     if (index ('+-', word(iw)(1:jc)) .ne. 0) then
                        o1 = iw
                     endif
                     iw = iw + 1
                  enddo
                  if (o1 .eq. 0) then
                     iw = s1+1
                     do while (iw .le. s2-1 .and. o1 .eq. 0)
                        jc = lastch (word(iw))
                        if (index ('><', word(iw)(1:jc)) .ne. 0) then
                           o1 = iw
                        endif
                        iw = iw + 1
                     enddo
                  endif
               endif
            endif
         endif
         if (o1 .gt. 0) then
            if (o1 .gt. 2) then
C
C              Decode left operand
C
               xl = ftn_atof(word(o1-1))
            else
               jc = lastch (word(o1))
               write (errbuf(1), 206) word(o1)(1:jc)
  206          format (' Operator (', a, ') is missing left operand')
               call prterx ('W', 1)
               error = 1
               go to 900
            endif
C
C           Find right operand.
C
  208       if (o1 .lt. nwrd - 1) then
C
C              Decode right operand
C
               xr = ftn_atof(word(o1+1))
            else
               jc = lastch (word(o1))
               write (errbuf(1), 214) word(o1)(1:jc)
  214          format (' Operator (', a, ') is missing right operand')
               call prterx ('W', 1)
               error = 1
               go to 900
            endif
C
C           Perform operation.
C
  216       if (word(o1) .eq. '**' .or. word(o1) .eq. '^') then
               if (xl .ne. 0.0) then
                  zc = cmplx (xr,0.0) * clog (cmplx(xl,0.0))
                  zc = cexp (zc)
                  x = real (zc)
               else
                  x = 0.0
               endif
            else if (word(o1) .eq. '*') then
               x = xl * xr
            else if (word(o1) .eq. '/') then
               if (xr .ne. 0.0) then
                  x = xl / xr
               else
                  jc = lastch (word(o1+1))
                  write (errbuf(1), 220) word(o1+1)(1:jc)
  220             format (' Error division by zero (', a, ').')
                  call prterx ('W', 1)
                  error = 1
               endif
            else if (word(o1) .eq. '+') then
               x = xl + xr
            else if (word(o1) .eq. '-') then
               x = xl - xr
            else if (word(o1) .eq. '>') then
               x = dim (xl, xr)
            else if (word(o1) .eq. '<') then
               x = -dim (xl, xr)
            else
               jc = lastch (word(o1))
               write (errbuf(1), 230) word(o1)(1:jc)
  230          format (' Unknown operator (', a, ')')
               call prterx ('W', 1)
               error = 0
               go to 900
            endif
            write (word(o1-1), 240) x
  240       format (e15.8)
            do i = o1, nwrd-2
               word(i) = word(i+2)
            enddo
            nwrd = nwrd - 2
            go to 160
         else

C           No operator: remove inner level of parenthesis.
C
            do i = s2, nwrd-1
               word(i) = word(i+1)
            enddo
            nwrd = nwrd - 1
            do i = s1, nwrd-1
               word(i) = word(i+1)
            enddo
            nwrd = nwrd - 1
C
C           Determine if intrinsic function.
C
            if (s1 .gt. 1) then
               symbol = word(s1-1)
               lc = lastch(symbol)
               if (symbol(1:lc) .eq. 'SIN') then
                  if (word(s1) .ne. ' ') then
                     do while (lastch (word(s1)) .lt. 10)
                        tempc = ' ' // word(s1)
                        word(s1) = tempc
                     enddo
                  endif
	 	  x = ftn_atof(word(s1))
                  write (word(s1-1), 240) sin (x)
                  do i = s1, nwrd
                     word(i) = word(i+1)
                  enddo
                  nwrd = nwrd - 1
               else if (symbol(1:lc) .eq. 'COS') then
                  if (word(s1) .ne. ' ') then
                     do while (lastch (word(s1)) .lt. 10)
                        tempc = ' ' // word(s1)
                        word(s1) = tempc
                     enddo
                  endif
	 	  x = ftn_atof(word(s1))
                  write (word(s1-1), 240) cos (x)
                  do i = s1, nwrd
                     word(i) = word(i+1)
                  enddo
                  nwrd = nwrd - 1
               else if (symbol(1:lc) .eq. 'TAN') then
                  if (word(s1) .ne. ' ') then
                     do while (lastch (word(s1)) .lt. 10)
                        tempc = ' ' // word(s1)
                        word(s1) = tempc
                     enddo
                  endif
	 	  x = ftn_atof(word(s1))
                  write (word(s1-1), 240) tan (x)
                  do i = s1, nwrd
                     word(i) = word(i+1)
                  enddo
                  nwrd = nwrd - 1
               else if (symbol(1:lc) .eq. 'ARCSIN') then
                  if (word(s1) .ne. ' ') then
                     do while (lastch (word(s1)) .lt. 10)
                        tempc = ' ' // word(s1)
                        word(s1) = tempc
                     enddo
                  endif
	 	  x = ftn_atof(word(s1))
                  write (word(s1-1), 240) asin (x)
                  do i = s1, nwrd
                     word(i) = word(i+1)
                  enddo
                  nwrd = nwrd - 1
               else if (symbol(1:lc) .eq. 'ARCCOS') then
                  if (word(s1) .ne. ' ') then
                     do while (lastch (word(s1)) .lt. 10)
                        tempc = ' ' // word(s1)
                        word(s1) = tempc
                     enddo
                  endif
	 	  x = ftn_atof(word(s1))
                  write (word(s1-1), 240) acos (x)
                  do i = s1, nwrd
                     word(i) = word(i+1)
                  enddo
                  nwrd = nwrd - 1
               else if (symbol(1:lc) .eq. 'ARCTAN') then
                  if (word(s1) .ne. ' ') then
                     do while (lastch (word(s1)) .lt. 10)
                        tempc = ' ' // word(s1)
                        word(s1) = tempc
                     enddo
                  endif
	 	  x = ftn_atof(word(s1))
                  write (word(s1-1), 240) atan (x)
                  do i = s1, nwrd
                     word(i) = word(i+1)
                  enddo
                  nwrd = nwrd - 1
               else if (symbol(1:lc) .eq. 'ABS') then
                  if (word(s1) .ne. ' ') then
                     do while (lastch (word(s1)) .lt. 10)
                        tempc = ' ' // word(s1)
                        word(s1) = tempc
                     enddo
                  endif
	 	  x = ftn_atof(word(s1))
                  write (word(s1-1), 240) abs (x)
                  do i = s1, nwrd
                     word(i) = word(i+1)
                  enddo
                  nwrd = nwrd - 1
               endif
            endif
            go to 160
         endif

      else
C
C        Evaluation complete. Decode results.
C
	 evalfn = ftn_atof(word(1))
      endif

  900 continue
      return
      end

max_stars_repo_name: mbheinen/bpa-ipf-tsp

max_issues_repo_name: cuihantao/bpa-ipf-tsp

max_forks_repo_name: mbheinen/bpa-ipf-tsp


text:
Function pyspen(xrein, ximin, ireim)
  Common /ludat1/mstu(200), paru(200), mstj(200), parj(200)
  Save /ludat1/
  Dimension b(0:14)
  Data b/1.000000E+00, -5.000000E-01, 1.666667E-01, 0.000000E+00, -3.333333E-02, 0.000000E+00, 2.380952E-02, 0.000000E+00, -3.333333E-02, 0.000000E+00, 7.575757E-02, 0.000000E+00, -2.531135E-01, 0.000000E+00, 1.166667E+00/
  xre = xrein
  xim = ximin
  If (abs(1.-xre)<1.E-6 .And. abs(xim)<1.E-6) Then
    If (ireim==1) pyspen = paru(1)**2/6.
    If (ireim==2) pyspen = 0.
    Return
  End If
  xmod = sqrt(xre**2+xim**2)
  If (xmod<1.E-6) Then
    If (ireim==1) pyspen = 0.
    If (ireim==2) pyspen = 0.
    Return
  End If
  xarg = sign(acos(xre/xmod), xim)
  sp0re = 0.
  sp0im = 0.
  sgn = 1.
  If (xmod>1.) Then
    algxre = log(xmod)
    algxim = xarg - sign(paru(1), xarg)
    sp0re = -paru(1)**2/6. - (algxre**2-algxim**2)/2.
    sp0im = -algxre*algxim
    sgn = -1.
    xmod = 1./xmod
    xarg = -xarg
    xre = xmod*cos(xarg)
    xim = xmod*sin(xarg)
  End If
  If (xre>0.5) Then
    algxre = log(xmod)
    algxim = xarg
    xre = 1. - xre
    xim = -xim
    xmod = sqrt(xre**2+xim**2)
    xarg = sign(acos(xre/xmod), xim)
    algyre = log(xmod)
    algyim = xarg
    sp0re = sp0re + sgn*(paru(1)**2/6.-(algxre*algyre-algxim*algyim))
    sp0im = sp0im - sgn*(algxre*algyim+algxim*algyre)
    sgn = -sgn
  End If
  xre = 1. - xre
  xim = -xim
  xmod = sqrt(xre**2+xim**2)
  xarg = sign(acos(xre/xmod), xim)
  zre = -log(xmod)
  zim = -xarg
  spre = 0.
  spim = 0.
  savere = 1.
  saveim = 0.
  Do i = 0, 14
    termre = (savere*zre-saveim*zim)/float(i+1)
    termim = (savere*zim+saveim*zre)/float(i+1)
    savere = termre
    saveim = termim
    spre = spre + b(i)*termre
    spim = spim + b(i)*termim
  End Do
  If (ireim==1) pyspen = sp0re + sgn*spre
  If (ireim==2) pyspen = sp0im + sgn*spim
  Return
End Function pyspen

max_stars_repo_name: xiaohaijin/AMPT

max_issues_repo_name: xiaohaijin/AMPT

max_forks_repo_name: xiaohaijin/AMPT


!>
!! Evaluate the linear incident wave field  from
!! linear theory (eta, phi and derivatives)
!!
!!    Evaluation on the free surface
!!    Evaluation on boundaries
!<
!
subroutine incident_linear_wf_finite_standing(wf_type, Wavefield, & !ramp_type,&
	 nfx,xp,nfy,yp,nfz,zp,bottom,tp,wavenum,grav,H, d)
!
USE PRECISION
USE GlobalVariables, ONLY: GhostGridX,GhostGridY,GhostGridZ, swenseTransientTime, &
     West_refl, East_refl, North_refl, South_refl, curvilinearONOFF, SFSOL
USE DataTypes
USE Constants
!
IMPLICIT NONE
INTERFACE
	SUBROUTINE phi_linear_3D_standing(H,angle,Lx,Ly,x,y,z,d,omega,time,grav,phix,phiy,phiz,phi,phit)
    USE PRECISION
   	REAL(KIND=long), INTENT(IN)  :: H,Lx,Ly,x,y,z,d,omega,time,grav,angle
   	REAL(KIND=long), INTENT(OUT) :: phix,phiy,phiz
   	REAL(KIND=long), INTENT(OUT), OPTIONAL :: phi,phit
   	END SUBROUTINE phi_linear_3D_standing
END INTERFACE
!
INTEGER nfx, nfy, nfz, nbp, wf_type
REAL(KIND=long) tp, wavenum, grav, H, d
! For curvilinear
REAL(KIND=long), DIMENSION(nfx,nfy) :: xp,yp
REAL(KIND=long), DIMENSION(nfz) :: zp
REAL(KIND=long), DIMENSION(nfx,nfy) :: bottom

!
TYPE(Wavefield_FS) :: Wavefield
! Local variables
INTEGER i, j, k, tmpIdx, Gidx3
REAL(KIND=long) temp1,omega,angle
REAL(KIND=long) zbp
REAL(KIND=long) FAC
!
!
! Time ramp and its derivatives
!time_ramp(:)=ramp_value(ramp_type,tp)
IF (swenseTransientTime/=zero) THEN
  IF(tp==zero) THEN
    FAC = zero
  ELSEIF(tp<=swenseTransientTime) THEN
    !FAC = one/swenseTransientTime*tp
    FAC = -two/(swenseTransientTime**3)*tp**3+three/(swenseTransientTime**2)*tp**2
  ELSE
    FAC = one
  ENDIF
ELSE
  FAC=one
ENDIF
!
omega = sqrt(grav*wavenum*tanh(wavenum*d))
omega=wf_type/abs(wf_type)*omega
!
SFSOL%k=SQRT((two*pi/SFSOL%L)**2+(zero*two*pi/SFSOL%L)**2)
SFSOL%T=SQRT((two*pi)**2/(grav*SFSOL%k*TANH(SFSOL%k*d)) );
SFSOL%c=SQRT(grav/SFSOL%k*TANH(SFSOL%k*d))
      !
      IF (curvilinearONOFF.EQ.0) THEN
        IF((nfx.EQ.1).OR.(nfy.EQ.1)) THEN
          write(*,*) 'This is a 3d case !!'
          stop
        ENDIF
      ELSE
        ! FIXME: 3D only for now
		IF((nfx.EQ.1).OR.(nfy.EQ.1)) THEN
          write(*,*) 'The 2D case in curvilinear SEENSE has to be done...'
          stop
        ENDIF
      ENDIF
      !
      IF (curvilinearONOFF.EQ.1) THEN
      	angle = 20.0_long*PI/180.0_long
      ELSE
        angle= zero
      ENDIF
	  !
      DO j=1,nfx
        DO k=1,nfy
            CALL eta_linear_3D_standing(H,angle,SFsol%L,SFsol%L,xp(j,k),yp(j,k),omega,tp,Wavefield%E_I(j,k),Wavefield%Ex_I(j,k), &
            	 Wavefield%Ey_I(j,k),Wavefield%Exx_I(j,k),Wavefield%Eyy_I(j,k),Wavefield%Et_I(j,k))
        	!
        	! Here phi, is computed on the free-surface z=0
        	!
        	CALL phi_linear_3D_standing(H,angle,SFsol%L,SFsol%L,xp(j,k),yp(j,k),zero,d,omega,tp,grav, &
            	Wavefield%Px_I_s(j,k),Wavefield%Py_I_s(j,k),Wavefield%Pz_I_s(j,k),Wavefield%P_I_s(j,k),Wavefield%Pt_I_s(j,k))
        ENDDO
      END DO
      !
	!
    ! Determine the incident wavefield on the boundary points...
    !
	!
    Wavefield%SourceEx = -FAC*Wavefield%Ex_I
    Wavefield%SourcePx = -FAC*Wavefield%Px_I_s
    Wavefield%SourceEy = -FAC*Wavefield%Ey_I
    Wavefield%SourcePy = -FAC*Wavefield%Py_I_s
    !
  	! FIXME: Assumption here is that the plane is aligned with the x- and y-axes. Extend to
  	!        general curvilinear coordinates.
    ! FIXME: Optimization of loops
    Gidx3 = 1
    !
    ! Bottom boundary
    !
    DO j = 1+GhostGridY, nfy-GhostGridY
		DO i = 1+GhostGridX, nfx-GhostGridX
            ! Wave elevation + derivative (Dimensional)
       		Wavefield%E_I_bp(Gidx3)  = FAC*Wavefield%E_I(i,j)
        	Wavefield%Ex_I_bp(Gidx3) = FAC*Wavefield%Ex_I(i,j)
            Wavefield%Ey_I_bp(Gidx3) = FAC*Wavefield%Ey_I(i,j)
            ! transform on the real z-space... Dimensional
            zbp=bottom(i,j)*zp(1+GhostGridZ)-bottom(i,j)
            !zbp=d*zp(2)-d
            CALL phi_linear_3D_standing(FAC*H,angle,SFsol%L,SFsol%L,xp(i,j),yp(i,j),zbp,d,omega,tp,grav, &
        			Wavefield%Px_I_bp(Gidx3),Wavefield%Py_I_bp(Gidx3),Wavefield%Pz_I_bp(Gidx3))
            ! Global index for Boundary Points, which option ? for now cf NormalX,Y,Z
            Wavefield%GidxTableBP(1,i,j) = Gidx3
            !Wavefield%GidxTableBP(1+GhostGridZ,i,j) = Gidx3
			! Increment global index
            Gidx3 = Gidx3+1
		END DO
	END DO
    !
   	! West  boundary
	!
    i = 1
	!
    DO j = 1+GhostGridY, nfy-GhostGridY
        DO k = 1, nfz
            ! Wave elevation + derivative
            Wavefield%E_I_bp(Gidx3)  = FAC*Wavefield%E_I(i+GhostGridX,j)
            Wavefield%Ex_I_bp(Gidx3) = FAC*Wavefield%Ex_I(i+GhostGridX,j)
            Wavefield%Ey_I_bp(Gidx3) = FAC*Wavefield%Ey_I(i+GhostGridX,j)
            IF (West_refl==0) THEN
				Wavefield%Px_I_bp(Gidx3)=zero; Wavefield%Py_I_bp(Gidx3)=zero; Wavefield%Pz_I_bp(Gidx3)=zero
            ELSE
                ! transform on the real z-space... Dimensional
                zbp=bottom(i+GhostGridX,j)*zp(k)-bottom(i+GhostGridX,j)
                !zbp=d*zp(k)-d
                CALL phi_linear_3D_standing(FAC*H,angle,SFsol%L,SFsol%L,xp(i+GhostGridX,j),yp(i+GhostGridX,j),zbp,d,omega,tp,grav, &
                	Wavefield%Px_I_bp(Gidx3),Wavefield%Py_I_bp(Gidx3),Wavefield%Pz_I_bp(Gidx3))
            ENDIF
            ! Global index for Boundary Points, which option ? for now cf NormalX,Y,Z
            Wavefield%GidxTableBP(k,i,j) = Gidx3
            !Wavefield%GidxTableBP(k,i+GhostGridX,j) = Gidx3
            ! Increment global index
            Gidx3 = Gidx3+1
        END DO
    END DO
	! East  boundary
   	i = nfx
    ! index of the first element (=Nz*(Ny-2*GhostY))
    tmpIdx = Gidx3
	!
    DO j = 1+GhostGridY, nfy-GhostGridY
        DO k = 1, nfz
            ! Wave elevation + derivative
            Wavefield%E_I_bp(Gidx3)  = FAC*Wavefield%E_I(i-GhostGridX,j)
            Wavefield%Ex_I_bp(Gidx3) = FAC*Wavefield%Ex_I(i-GhostGridX,j)
            Wavefield%Ey_I_bp(Gidx3) = FAC*Wavefield%Ey_I(i-GhostGridX,j)
            IF (East_refl==0) THEN
				Wavefield%Px_I_bp(Gidx3)=zero; Wavefield%Py_I_bp(Gidx3)=zero; Wavefield%Pz_I_bp(Gidx3)=zero
           	ELSE
                ! transform on the real z-space... Dimensional
                zbp=bottom(i-GhostGridX,j)*zp(k)-bottom(i-GhostGridX,j)
                !zbp=d*zp(k)-d
                CALL phi_linear_3D_standing(FAC*H,angle,SFsol%L,SFsol%L,xp(i-GhostGridX,j),yp(i-GhostGridX,j),zbp,d,omega,tp,grav, &
                    Wavefield%Px_I_bp(Gidx3),Wavefield%Py_I_bp(Gidx3),Wavefield%Pz_I_bp(Gidx3))
            ENDIF
            ! Global index for Boundary Points, which option ? for now cf NormalX,Y,Z
            Wavefield%GidxTableBP(k,i,j) = Gidx3
            !Wavefield%GidxTableBP(k,i-GhostGridX,j) = Gidx3
            ! Increment global index
            Gidx3 = Gidx3+1
        END DO
    END DO
    !
    ! South boundary
    !
    j = 1
    DO i = 1+GhostGridX, nfx-GhostGridX
        ! index of the first element (=2*Nz*(Ny-2*GhostY)+i)
        tmpIdx = Gidx3
        ! Wave elevation + derivative (Dimensional)
        Wavefield%E_I_bp(Gidx3)  = FAC*Wavefield%E_I(i,j+GhostGridY)
        Wavefield%Ex_I_bp(Gidx3) = FAC*Wavefield%Ex_I(i,j+GhostGridY)
        Wavefield%Ey_I_bp(Gidx3) = FAC*Wavefield%Ey_I(i,j+GhostGridY)
        DO k = 1, nfz
            ! Wave elevation + derivative
            Wavefield%E_I_bp(Gidx3)  = Wavefield%E_I_bp(tmpIdx)
            Wavefield%Ex_I_bp(Gidx3) = Wavefield%Ex_I_bp(tmpIdx)
            Wavefield%Ey_I_bp(Gidx3) = Wavefield%Ey_I_bp(tmpIdx)
            IF (South_refl==0) THEN
                Wavefield%Px_I_bp(Gidx3)=zero; Wavefield%Py_I_bp(Gidx3)=zero; Wavefield%Pz_I_bp(Gidx3)=zero
            ELSE
                ! transform on the real z-space... Dimensional
                zbp=bottom(i,j+GhostGridY)*zp(k)-bottom(i,j+GhostGridY)
                !zbp=d*zp(k)-d
                CALL phi_linear_3D_standing(FAC*H,angle,SFsol%L,SFsol%L,xp(i,j+GhostGridY),yp(i,j+GhostGridY),zbp,d,omega,tp,grav, &
                    Wavefield%Px_I_bp(Gidx3),Wavefield%Py_I_bp(Gidx3),Wavefield%Pz_I_bp(Gidx3))
            ENDIF
            ! Global index for Boundary Points, which option ? for now cf NormalX,Y,Z
            Wavefield%GidxTableBP(k,i,j) = Gidx3
            !Wavefield%GidxTableBP(k,i,j+GhostGridY) = Gidx3
            ! Increment global index
            Gidx3 = Gidx3+1
        END DO
    END DO
    !
    ! North boundary
    j = nfy
    DO i = 1+GhostGridX, nfx-GhostGridX
        ! index of the first element (=2*Nz*(Ny-2*GhostY)+...+i)
        tmpIdx = Gidx3
        ! Wave elevation + derivative (Dimensional)
        Wavefield%E_I_bp(Gidx3)  = FAC*Wavefield%E_I(i,j-GhostGridY)
        Wavefield%Ex_I_bp(Gidx3) = FAC*Wavefield%Ex_I(i,j-GhostGridY)
        Wavefield%Ey_I_bp(Gidx3) = FAC*Wavefield%Ey_I(i,j-GhostGridY)
        DO k = 1, nfz
            ! Wave elevation + derivative
            Wavefield%E_I_bp(Gidx3)  = Wavefield%E_I_bp(tmpIdx)
            Wavefield%Ex_I_bp(Gidx3) = Wavefield%Ex_I_bp(tmpIdx)
            Wavefield%Ey_I_bp(Gidx3) = Wavefield%Ey_I_bp(tmpIdx)
            IF (North_refl==0) THEN
                Wavefield%Px_I_bp(Gidx3)=zero; Wavefield%Py_I_bp(Gidx3)=zero; Wavefield%Pz_I_bp(Gidx3)=zero
            ELSE
              ! transform on the real z-space... Dimensional
              zbp=bottom(i,j-GhostGridY)*zp(k)-bottom(i,j-GhostGridY)
              !zbp=d*zp(k)-d
              CALL phi_linear_3D_standing(FAC*H,angle,SFsol%L,SFsol%L,xp(i,j-GhostGridY),yp(i,j-GhostGridY),zbp,d,omega,tp,grav, &
                  Wavefield%Px_I_bp(Gidx3),Wavefield%Py_I_bp(Gidx3),Wavefield%Pz_I_bp(Gidx3))
            ENDIF
            ! Global index for Boundary Points, which option ? for now cf NormalX,Y,Z
            Wavefield%GidxTableBP(k,i,j) = Gidx3
            !Wavefield%GidxTableBP(k,i,j-GhostGridY) = Gidx3
            ! Increment global index
            Gidx3 = Gidx3+1
        END DO
    END DO
    ! FIXME: How to treat the corners ?
    IF (West_refl==0) THEN
      Wavefield%SourceEx(1+GhostGridX,2+GhostGridY:nfy-1-GhostGridY) = zero
      Wavefield%SourcePx(1+GhostGridX,2+GhostGridY:nfy-1-GhostGridY) = zero
      Wavefield%SourceEy(1+GhostGridX,2+GhostGridY:nfy-1-GhostGridY) = zero
      Wavefield%SourcePy(1+GhostGridX,2+GhostGridY:nfy-1-GhostGridY) = zero
    ENDIF
    IF (East_refl==0) THEN
      Wavefield%SourceEx(nfx-GhostGridX,2+GhostGridY:nfy-1-GhostGridY) = zero
      Wavefield%SourcePx(nfx-GhostGridX,2+GhostGridY:nfy-1-GhostGridY) = zero
      Wavefield%SourceEy(nfx-GhostGridX,2+GhostGridY:nfy-1-GhostGridY) = zero
      Wavefield%SourcePy(nfx-GhostGridX,2+GhostGridY:nfy-1-GhostGridY) = zero
    ENDIF
    IF (South_refl==0) THEN
      Wavefield%SourceEx(2+GhostGridX:nfx-1-GhostGridX,1+GhostGridY) = zero
      Wavefield%SourcePx(2+GhostGridX:nfx-1-GhostGridX,1+GhostGridY) = zero
      Wavefield%SourceEy(2+GhostGridX:nfx-1-GhostGridX,1+GhostGridY) = zero
      Wavefield%SourcePy(2+GhostGridX:nfx-1-GhostGridX,1+GhostGridY) = zero
    ENDIF
    IF (North_refl==0) THEN
      Wavefield%SourceEx(2+GhostGridX:nfx-1-GhostGridX,nfy-GhostGridY) = zero
      Wavefield%SourcePx(2+GhostGridX:nfx-1-GhostGridX,nfy-GhostGridY) = zero
      Wavefield%SourceEy(2+GhostGridX:nfx-1-GhostGridX,nfy-GhostGridY) = zero
      Wavefield%SourcePy(2+GhostGridX:nfx-1-GhostGridX,nfy-GhostGridY) = zero
    ENDIF
    ! Corners are set to zero if both reflexions are zero...
    IF (West_refl==0) THEN
      IF (South_refl==0) THEN
      	Wavefield%SourceEx(1+GhostGridX,1+GhostGridY) = zero
      	Wavefield%SourcePx(1+GhostGridX,1+GhostGridY) = zero
      	Wavefield%SourceEy(1+GhostGridX,1+GhostGridY) = zero
      	Wavefield%SourcePy(1+GhostGridX,1+GhostGridY) = zero
      ENDIF
      IF (North_refl==0) THEN
      	Wavefield%SourceEx(1+GhostGridX,nfy-GhostGridY) = zero
      	Wavefield%SourcePx(1+GhostGridX,nfy-GhostGridY) = zero
      	Wavefield%SourceEy(1+GhostGridX,nfy-GhostGridY) = zero
      	Wavefield%SourcePy(1+GhostGridX,nfy-GhostGridY) = zero
      ENDIF
    ENDIF
    IF (East_refl==0) THEN
      IF (South_refl==0) THEN
      	Wavefield%SourceEx(nfx-GhostGridX,1+GhostGridY) = zero
      	Wavefield%SourcePx(nfx-GhostGridX,1+GhostGridY) = zero
      	Wavefield%SourceEy(nfx-GhostGridX,1+GhostGridY) = zero
      	Wavefield%SourcePy(nfx-GhostGridX,1+GhostGridY) = zero
      ENDIF
      IF (North_refl==0) THEN
      	Wavefield%SourceEx(nfx-GhostGridX,nfy-GhostGridY) = zero
      	Wavefield%SourcePx(nfx-GhostGridX,nfy-GhostGridY) = zero
      	Wavefield%SourceEy(nfx-GhostGridX,nfy-GhostGridY) = zero
      	Wavefield%SourcePy(nfx-GhostGridX,nfy-GhostGridY) = zero
      ENDIF
    ENDIF

RETURN
END SUBROUTINE incident_linear_wf_finite_standing


text:
!> @ingroup #<sort><name>#
!>
!> @brief
!> Run #<sort><docRef># on a vector of #<data><full># values.
!>
!> @details
!> Run #<sort><docRef># on a vector of #<data><full># values to sort them in a specified order. #IF<map><true>THEN<The vector of intergers #<map><name># is permuted in the same way that the vector vec is permuted to be sorted. It is important that #<map><var># be of equal or greather length than vec.>ENDIF#
!>
!> @param #<data><var>#
!> The vector to be sorted.
!>
#IF<map><true>#!> @param #<map><var>#
#IF<map><true>#!> Vector of intergers which is permuted in the same way that #<data><var># is permuted. If #<map><name># is a vector of values 1 to N, where N is the length of #<data><var>#, than, #<map><name># is the map from #<data><var>#'s original value to #<data><var>#'s final value after being sort.
#IF<map><true>#!>
!> @param ascending
!> Indicates if the ordering is ascending or decending.
!>
!> | Value | Meaning         |
!> |:-----:| :-------------- |
!> | true  | Ascending order |
!> | false | Decending order |
!>
#IF<map><true>#!> @warning
#IF<map><true>#!> The length of #<map><var># is not checked within this subroutine. If the length of #<data><var># is larger than the length of #<map><name>#, than this subroutine may attempt to access elements of #<map><name># that do not exist.
#IF<map><true>#!>
!> @author Arin R. Bratt
!> @date 11/07/2014

pure subroutine #<sort><name>##map##data#(#<data><var>##<map><arg>#,ascending)
  ! Arguments
  #<data><full>#, intent(inout) :: #<data><var>#(:)
  #IF<map><true>#integer, intent(inout) :: #<map><var>#(:)
  logical, intent(in) :: ascending

  ! Main
  if (ascending) then
     call #<sort><name>#_asc#map##data#(#<data><var>##<map><arg>#)
  else
     call #<sort><name>#_dec#map##data#(#<data><var>##<map><arg>#)
  end if
end subroutine #<sort><name>##map##data#

max_stars_repo_name: arinrb/bsort

max_issues_repo_name: arinrb/bsort

max_forks_repo_name: arinrb/bsort


text:
The West Nile Ramblers play western garage & country noir, a crunchy combination of western swing, gypsy folk, rock & roll, and jazz that inspires people to cut loose and have a good time. The Pepper Peddler Alex Roth and http://www.andylentz.com Andy Lentz bring their decadelong history of musical collaboration to bear in fronting the band and captivating the crowd. Their infectious energy brought them to the stage opening for Paul Brady as winners of the 2012 Ballyshannon Folk Festival band contest in Ireland and has kept them returning to Davis Music Fest the last three years as one of the most popular acts of the festival.

Their latest release, “Western Garage & Country Noir,” is an homage to the vinyl era that features an Aside of electrified studio recordings and a Bside of bluegrassstyle live tracks taped at a rustic cabin in Lake Tahoe. Strains of country, traditional, folk and funk are expertly handled in the songwriting and musicianship of Roth and Lentz, and by their stellar cast of accompanying players including Bart van der Zeeuw (formerly of Ks Choice), Katie Weed of Old Man Markley, and Cooper McBean of The Devil Makes Three.

Also, they have been known to bicycle to their gigs with all their gear.  Gotta give them some credit for that.

20120303 18:24:00 nbsp They were at the Davis Flea Market.  I liked them. Users/Mbrowne

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      FUNCTION GAMMP(A,X)
      IF(X.LT.0..OR.A.LE.0.)PAUSE
      IF(X.LT.A+1.)THEN
        CALL GSER(GAMSER,A,X,GLN)
        GAMMP=GAMSER
      ELSE
        CALL GCF(GAMMCF,A,X,GLN)
        GAMMP=1.-GAMMCF
      ENDIF
      RETURN
      END

ext: for


text:
!      Ken Klare, LANL P-4      (c)1990,1991
!***********************************************************************
!********************* 3.2 Cubic Spline Evaluation *********************
!-----------------------------------------------------------------------
      Real Function CSVAL(x,nintv,break,c)
      Real      x      !input, the point to evaluate
      Integer nintv      !input, the number of polynomial pieces
      Real      break(*)!input, [nintv+1] breakpoints
      Real      c(4,*)      !input, [4,nintv+1] polynomial coefficients
!+Cubic spline evaluation.
      Real      PPDER
      CSVAL = PPDER(0,x,4,nintv,break,c)
      End
!-----------------------------------------------------------------------
      Real Function CSDER(ideriv,x,nintv,break,c)
      Integer ideriv      !input, the order of deriavtive to be evaluated
      Real      x      !input, the point to evaluate
      Integer nintv      !input, the number of polynomial pieces
      Real      break(*)!input, [nintv+1] breakpoints
      Real      c(4,*)      !input, [4,nintv+1] polynomial coefficients
!+Cubic spline derivative evaluation.
      Real      PPDER
      CSDER = PPDER(ideriv,x,4,nintv,break,c)
      End
!-----------------------------------------------------------------------
      Real Function CSITG(a,b,nintv,break,c)
      Real      a,b      !input, the range to evaluate
      Integer nintv      !input, the number of polynomial pieces
      Real      break(*)!input, [nintv+1] breakpoints
      Real      c(4,*)      !input, [4,nintv+1] polynomial coefficients
!+Cubic spline evaluation.
      Real      PPITG
      CSITG = PPITG(a,b,4,nintv,break,c)
      End
!***********************************************************************
!********************* Vector extensions *******************************
!-----------------------------------------------------------------------
      Subroutine CSVALV(nx,x,nintv,break,c,y)
!     ENTRY CSVALV_(nx,x,nintv,break,c,y)
      Integer      nx      !input, the number of points
      Real      x(*)      !input, the [nx] points to evaluate
      Integer nintv      !input, the number of polynomial pieces
      Real      break(*)!input, [nintv+1] breakpoints
      Real      c(4,*)      !input, [4,nintv+1] polynomial coefficients
      Real      y(*)      !output, the evaluated points
!+Cubic spline evaluation.
      Call PPDERV(0,nx,x,4,nintv,break,c,y)
      End
!-----------------------------------------------------------------------
      Subroutine CSDERV(ideriv,nx,x,nintv,break,c,y)
!     ENTRY CSDERV_(ideriv,nx,x,nintv,break,c,y)
      Integer ideriv      !input, the order of deriavtive to be evaluated
      Integer      nx      !input, the number of points
      Real      x(*)      !input, the [nx] points to evaluate
      Integer nintv      !input, the number of polynomial pieces
      Real      break(*)!input, [nintv+1] breakpoints
      Real      c(4,*)      !input, [4,nintv+1] polynomial coefficients
      Real      y(*)      !output, the evaluated points
!+Cubic spline derivative evaluation.
      Call PPDERV(ideriv,nx,x,4,nintv,break,c,y)
      End
!-----------------------------------------------------------------------
      Subroutine CSITGV(na,a,nb,b,nintv,break,c,y)
!     ENTRY CSITGV_(na,a,nb,b,nintv,break,c,y)
      Integer      na      !input, the number of points
      Real      a(*)      !input, the [na] points to evaluate
      Integer      nb      !input, the number of points
      Real      b(*)      !input, the [nb] points to evaluate
      Integer nintv      !input, the number of polynomial pieces
      Real      break(*)!input, [nintv+1] breakpoints
      Real      c(4,*)      !input, [4,nintv+1] polynomial coefficients
      Real      y(*)      !output, the evaluated points
!+Cubic spline evaluation.
      Call PPITGV(na,a,nb,b,4,nintv,break,c,y)
      End

max_stars_repo_name: fluffynukeit/mdsplus

max_issues_repo_name: fluffynukeit/mdsplus

max_forks_repo_name: fluffynukeit/mdsplus


text:
C
      SUBROUTINE SNEWBILCAL (DGLNG,DGLAT,RADI,DYEAR,DALNG,
     $                       DALAT,DBABS,DBN,DBE,DBD,DXL,ISTAT)
c  ------------------------------------------------------------------------
c  Part of "poa_calfos" (STECF/IPMG)
c  Part of "NEWBILCAL" program to produce IGRF geomagnetic model (GSFC)
c
c  History:
c  --------
c  Version   Date        Author          Description
c      1.0       Jul 00  M. Rosa         POA calibration code
c      1.2       Jun 01  A. Alexov       Added IDS's common block for file
c                                        header access, ISTAT
c  ------------------------------------------------------------------------
C
C *************************************************************
C ** IGRF MAGNETIC FIELD MODEL +  SHELLG L-VALUE CALCULATION **
C *************************************************************
C adapted by Michael Rosa STECF 22 DEC 1996 from BILCAL, VERSION 3.0,
C AUGUST 1995, includes also AACGM code
C
C  UMR     =   ATAN(1.0)*4./180.   <DEGREE>*UMR=<RADIANT>
C  ERA         EARTH RADIUS FOR NORMALIZATION OF CARTESIAN COORDINATES
C  EREQU       MAJOR HALF AXIS FOR EARTH ELLIPSOID (6378.160 KM)
C  ERPOL       MINOR HALF AXIS FOR EARTH ELLIPSOID (6356.775 KM)
C  AQUAD       SQUARE OF MAJOR HALF AXIS FOR EARTH ELLIPSOID
C  BQUAD       SQUARE OF MINOR HALF AXIS FOR EARTH ELLIPSOID
C  ERA, EREQU and ERPOL as recommended by the IAU
C---------------------------------------------------------------------
      IMPLICIT NONE
      INTEGER*4   ICODE
      REAL*8 DGLNG,DGLAT,DALNG,DALAT,DBABS,DXL,RADI
      REAL*8 DBN,DBE,DBD,DIMO,DIP,DEC,DYEAR
c      REAL*8 DGLNG,DGLAT,RADI,DIMO,DIP,DEC

      REAL*4      GLNG,GLAT,ALAT,ALNG,BABS,BAB1,XL
      REAL*4      UMR,ERAD,AQUAD,BQUAD
      REAL*4      HEIGHT,YEAR,BN,BE,BD

      COMMON/GENER/       UMR,ERAD,AQUAD,BQUAD
C
C local variables
C
      INTEGER ISTAT
      CHARACTER*80 CONTXT
      INTEGER STDOUT
      PARAMETER (STDOUT = 1)
      INTEGER STDERR
      PARAMETER (STDERR = 2)
C
C Common block containing input/output file descriptions
C
C       IDS - file id numbers
C       GCOUNT - group counts
C       NAXIS - naxis
C       NAXIS1 - first dimensions
C       NAXIS2 - second dimensions
C       FILL - Fill values
C
        INTEGER IDS(30),NAXIS(30),NAXIS1(30),NAXIS2(30),GCOUNT(30)
        REAL FILL(30)
        COMMON /IOCOM/IDS,GCOUNT,NAXIS,NAXIS1,NAXIS2,FILL

C     initialize common /GENER/
      if (ERAD.LT.6300.0) CALL INITIZE
C
      HEIGHT = SNGL(RADI/1000.D0)-ERAD
      GLNG = SNGL(DGLNG)
      GLAT = SNGL(DGLAT)
      YEAR = SNGL(DYEAR)

C Get field coefficients for actual decimal year
      CALL FELDCOF (YEAR,DIMO,ISTAT)
      IF (ISTAT .NE. 0) THEN
             GO TO 999
      END IF
C Calculate field
      CALL FELDG (GLAT,GLNG,HEIGHT,BN,BE,BD,BABS)
C Calculate shell parameters
      CALL SHELLG (GLAT,GLNG,HEIGHT,DIMO,XL,ICODE,BAB1)

      DIP = ASIN(BD/BABS)/UMR
      DEC = ASIN(BE/SQRT(BE*BE+BN*BN))/UMR

      DBABS = DBLE(BABS)
      DBN = BN
      DBE = BE
      DBD = BD
      DXL = DBLE(XL)

      CALL CONV_GEO_CRD (GLAT,GLNG,HEIGHT,ALAT,ALNG,1,ISTAT)
C check the output error
      IF (ISTAT .NE. 0) THEN
             GO TO 999
      END IF

      DALNG = DBLE(ALNG)
      DALAT = DBLE(ALAT)

      IF (ISTAT .NE. 0) THEN
             GO TO 999
      END IF

      ISTAT=0
      GO TO 1000
999   CALL YMSPUT(CONTXT,STDOUT+STDERR,0,ISTAT)
      ISTAT=1
1000  RETURN
      END

max_stars_repo_name: iraf-community/stsdas


text:
!-------------------------------------------------------------------------
!         NASA/GSFC, Data Assimilation Office, Code 910.3, GEOS/DAS      !
!-------------------------------------------------------------------------
!BOP
!
! !MODULE: m_odsmeta --- Defines the ODS metadata conventions
!
! !INTERFACE:

      MODULE m_odsmeta

! !USES:

      Implicit none

! !DESCRIPTION:
! \label{MODS:Meta}
!  Defines ODS metadata conventions.
!
! !TO DO:
!
!  1. Define history and exclusion marks for superobbing
!  2. Implement bitwise encoding of history mark
!  3. Add kx/kt parameters
!
! !REVISION HISTORY:
!
!  01Feb99 (Dee)      See v1.1 revision history for previous revisions
!  13Dec99 (RT)       Added preprocessing flags defined by Meta and Yelena
!  20Dec99 (Dee/RT)   Removed restriction on preprocessing and online
!                      flag values; made it into a module.
!  23dec99  Todling/
!           da Silva  Created from qc_flag.h, added GEOS-4 specific flags.
!  28Dec99  Todling   Adding kt info (use to live in kt_max.h).
!  09Feb00  Dee       Added X_NOSTATS, X_PSAS
!  09Feb00  Todling   Defined kxmax here.
!  06Jun00  Todling   Updated kt list.
!  04Oct00  Todling   Increased kxmax to 512.
!  13May02  Dee       Modified some SQC flags and removed obsolete ones
!  03Jun02  Dee       Updated kt list; added X_BADLAYER
!  23Apr03  Sienkiewicz   Add history & exclusion marks for preprocessing
!  04Jun04  Todling   Added kxmod (per Jing Guo's m_Sndx)
!
!EOP
!-------------------------------------------------------------------------

      include 'odsmeta.h'

      end MODULE m_odsmeta

max_stars_repo_name: GEOS-ESM/AeroApps

max_forks_repo_name: GEOS-ESM/AeroApps


!----------------------------------------------------------------------------------------------
! MODULE: read_txt_mod
!> @author Marin Sapunar, Ruđer Bošković Institute
!> @date August, 2018
!
!> @brief Contains subroutines for reading arrays from simple ASCII files.
!----------------------------------------------------------------------------------------------
module read_txt_mod
    use global_defs
    implicit none

    private
    public :: read_txt

    !----------------------------------------------------------------------------------------------
    ! SUBROUTINE: read_txt
    !
    !> @brief Read a real(dp) array from a text file.
    !> @details
    !! The first line of the file should contain the dimensions of the array.
    !! The subsequent lines loop over dimensions of the array from the outermost to the second
    !! dimension. The first dimension is written as a vector in each line.
    !! If the array dimensions are given as input (have values different than zero on calling the
    !! subroutine), an error is thrown if the dimensions given in the file don't match the input.
    !----------------------------------------------------------------------------------------------
    interface read_txt
        module procedure read_txt_dim2
        module procedure read_txt_dim3
    end interface read_txt

contains

    !----------------------------------------------------------------------------------------------
    ! SUBROUTINE: read_txt_dim2
    !> @brief Read a 2-dimensional array from a text file. See read_txt.
    !----------------------------------------------------------------------------------------------
    subroutine read_txt_dim2(fname, n1, n2, array)
        character(len=*), intent(in) :: fname
        integer, intent(inout) :: n1
        integer, intent(inout) :: n2
        real(dp), allocatable, intent(out) :: array(:, :)
        integer :: d1, d2, iunit, iotest

        call check_file_exists(fname)
        open(newunit=iunit, file=fname)
        read(iunit, *, iostat=iotest) d1, d2
        if (iotest /= 0) call read_error_dim(fname, 2, iotest)
        if ((n1 /= 0) .and. (n1 /= d1)) call dim_error(fname, 1, n1, d1)
        if ((n2 /= 0) .and. (n2 /= d2)) call dim_error(fname, 2, n2, d2)
        n1 = d1
        n2 = d2
        if (allocated(array)) deallocate(array)
        allocate(array(n1, n2))
        read(iunit, *, iostat=iotest) array
        if (iotest /= 0) call read_error_val(fname, iotest)
        close(iunit)
    end subroutine read_txt_dim2

    !----------------------------------------------------------------------------------------------
    ! SUBROUTINE: read_txt_dim3
    !> @brief Read a 3-dimensional array from a text file. See read_txt.
    !----------------------------------------------------------------------------------------------
    subroutine read_txt_dim3(fname, n1, n2, n3, array)
        character(len=*), intent(in) :: fname
        integer, intent(inout) :: n1
        integer, intent(inout) :: n2
        integer, intent(inout) :: n3
        real(dp), allocatable, intent(out) :: array(:, :, :)
        integer :: d1, d2, d3, iunit, iotest

        call check_file_exists(fname)
        open(newunit=iunit, file=fname)
        read(iunit, *, iostat=iotest) d1, d2, d3
        if (iotest /= 0) call read_error_dim(fname, 3, iotest)
        if ((n1 /= 0) .and. (n1 /= d1)) call dim_error(fname, 1, n1, d1)
        if ((n2 /= 0) .and. (n2 /= d2)) call dim_error(fname, 2, n2, d2)
        if ((n3 /= 0) .and. (n3 /= d3)) call dim_error(fname, 3, n3, d3)
        n1 = d1
        n2 = d2
        n3 = d3
        if (allocated(array)) deallocate(array)
        allocate(array(n1, n2, n3))
        read(iunit, *, iostat=iotest) array
        if (iotest /= 0) call read_error_val(fname, iotest)
        close(iunit)
    end subroutine read_txt_dim3

    !----------------------------------------------------------------------------------------------
    ! SUBROUTINE: check_file_exists
    !> @brief Print error message and abort program if file is not found.
    !----------------------------------------------------------------------------------------------
    subroutine check_file_exists(file_name, error_message)
        character(len=*), intent(in) :: file_name
        character(len=*), intent(in), optional :: error_message
        logical :: chk

        inquire(file=file_name, exist=chk)
        if (.not. chk) then
            write(stderr, *) ' Failed to find required file: '//file_name
            if (present(error_message)) write(stderr, *) error_message
            call abort()
        end if
    end subroutine check_file_exists

    !----------------------------------------------------------------------------------------------
    ! SUBROUTINE: dim_error
    !> @brief Exit the program with an error if an array dimension doesn't match expected value.
    !----------------------------------------------------------------------------------------------
    subroutine dim_error(fname, cdim, n, d)
        character(len=*), intent(in) :: fname
        integer, intent(in) :: cdim
        integer, intent(in) :: n
        integer, intent(in) :: d

        write(stderr, '(3a)') 'Error reading array from file ', trim(fname), '.'
        write(stderr, '(1x, a, i0)') 'For array dimension: ', cdim
        write(stderr, '(3x, a, i0, a, i0, a)') 'Expected size ', n, ' instead of ', d, '.'
        call abort()
    end subroutine dim_error

    !----------------------------------------------------------------------------------------------
    ! SUBROUTINE: read_error_dim
    !> @brief Exit the program with an error if array dimensions can't be read from the first line.
    !----------------------------------------------------------------------------------------------
    subroutine read_error_dim(fname, cdim, io)
        character(len=*), intent(in) :: fname
        integer, intent(in) :: cdim
        integer, intent(in) :: io

        write(stderr, '(3a)') 'Error reading array dimensions from file ', trim(fname), '.'
        write(stderr, '(1x, a, i0, a)') 'First line should contain ', cdim, &
        &                               ' integers with array dimensions.'
        write(stderr, '(1x,a,i0)') 'iostat: ', io
        call abort()
    end subroutine read_error_dim

    !----------------------------------------------------------------------------------------------
    ! SUBROUTINE: read_error_val
    !> @brief Exit the program with an error if array dimensions can't be read from the first line.
    !----------------------------------------------------------------------------------------------
    subroutine read_error_val(fname, io)
        character(len=*), intent(in) :: fname
        integer, intent(in) :: io

        write(stderr, '(3a)') 'Error reading array values from file ', trim(fname), '.'
        write(stderr, '(1x,a,i0)') 'iostat: ', io
        call abort()
    end subroutine read_error_val

end module read_txt_mod


text:
      integer function pt_handle_rules(tgt_info,cur_target,wlist,mode)

      implicit none

      include 'opdim.h'
      include 'stdunit.h'
      include 'mdef_target_info.h'
      include 'def_word_list.h'
      include 'par_parse_targets.h'
      include 'ifc_targets.h'

      character(len=15), parameter ::
     &    i_am = 'pt_handle_rules'

      type(target_info), intent(inout) ::
     &    tgt_info
      type(word_list), intent(inout) ::
     &    wlist
      integer, intent(inout) ::
     &    mode
      character(len=*), intent(inout) ::
     &    cur_target

      integer ::
     &    error, arg_type, idx_arg, idx_cmd, ndim, idim, nargs
      character ::
     &    sep
      character(len=256) ::
     &    command, label, word

      type(action), pointer ::
     &    cmd_pnt
      type(action_arg), pointer ::
     &    arg_pnt

      integer, pointer ::
     &    list_int(:)
      real(8), pointer ::
     &    list_rl8(:)
      logical, pointer ::
     &    list_log(:), need_arg(:)
      character(len=len_target_name), pointer ::
     &    list_label(:)

      logical, external ::
     &    advance_word_list_entry
      integer, external ::
     &    idx_action, idx_command_arg

      error = no_error

      ! get current entry again
      call get_word_list_entry(command,sep,wlist)

      ! compare with all the defined rules
      call uppcas(command)
      idx_cmd = idx_action(command,tgt_info)

      if (idx_cmd.le.0) then
        error = error_unknown_comkey
      else if (mode.ne.scan_for_rules) then
        error = error_here_no_rule
      else if (sep.ne.'(') then
        error = error_exp_bob
      end if

      pt_handle_rules = error
      if (error.ne.no_error) return

      ! if these tests are passed, we can register the rule
      call set_rule2(cur_target,command,tgt_info)

      cmd_pnt => tgt_info%act_array(idx_cmd)%act

      ! make an array to identify required arguments
      nargs = cmd_pnt%n_arguments
      allocate(need_arg(max(1,nargs)))
      do idx_arg = 1, nargs
        arg_pnt => cmd_pnt%arg(idx_arg)
        need_arg(idx_arg) = arg_pnt%required
      end do

      ! now process the list of arguments
      if (.not.advance_word_list_entry(wlist,'d')) then
        call quit(1,i_am,'unexpected absence of d-branch in wlist')
      end if

      scan_arg: do
        call get_word_list_entry(label,sep,wlist)
        ! accept a line-break
        if (len_trim(label).eq.0.and.sep.eq.'E') then
          if (.not.advance_word_list_entry(wlist,' ')) then
            pt_handle_rules = error_unexp_eob
            return
          end if
          call get_word_list_entry(label,sep,wlist)
        end if

        call uppcas(label)
        ! check whether the argument label exists
        ! and get the type of expected argument
        idx_arg = idx_command_arg(label,cmd_pnt)
        if (idx_arg.le.0) then
          pt_handle_rules = error_unknown_arglab
          return
        end if
        arg_pnt => cmd_pnt%arg(idx_arg)
        arg_type = arg_pnt%type
        ! remember that argument is set
        need_arg(idx_arg) = .false.

        if (sep.eq.')'.and.trim(label).eq.'') exit scan_arg

        if (sep.ne.'=') then
          pt_handle_rules = error_exp_eq
          return
        end if

        ! get next entry
        if (.not.advance_word_list_entry(wlist,' ')) then
          error = error_exp_arg
        end if

        call get_word_list_entry(word,sep,wlist)
        ! an array upcoming?
        if (trim(word).eq.'') then
          if (sep.ne.'(') then
            pt_handle_rules = error_exp_bob
            return
          end if

          ! move down and count entries
          if (.not.advance_word_list_entry(wlist,'d')) then
           call quit(1,i_am,'unexpected error when descending in wlist')
          end if

          ndim = 0
          cnt_lp: do
            call get_word_list_entry(word,sep,wlist)
            if (sep.ne.','.and.sep.ne.')') then
              pt_handle_rules = error_exp_copa
              return
            end if
            ndim = ndim+1
            if (.not.advance_word_list_entry(wlist,' ')) exit cnt_lp
          end do cnt_lp

          ! return to upper level
          if (.not.advance_word_list_entry(wlist,'u')) then
            call quit(1,i_am,'unexpected error when ascending in wlist')
          end if
          ! and go back
          if (.not.advance_word_list_entry(wlist,'d')) then
           call quit(1,i_am,'unexpected error when descending in wlist')
          end if
          ! (we now are back where we started!)

        else ! single argument
          ndim = 1
        end if

        ! allocate space for arguments (if needed)
        select case(arg_type)
        case(aatype_label)
          allocate(list_label(ndim))
        case(aatype_log)
          allocate(list_log(ndim))
        case(aatype_int)
          allocate(list_int(ndim))
        case(aatype_rl8)
          allocate(list_rl8(ndim))
        case(aatype_str)
          if (ndim.gt.1) then
            pt_handle_rules = error_no_str_arr
            return
          end if
        case default
          call quit(1,i_am,
     &              'unsupported argument type for '//trim(label))
        end select

        ! read in argument lists
        do idim = 1, ndim
          call get_word_list_entry(word,sep,wlist)
          select case(arg_type)
          case(aatype_label)
            list_label(idim) = word
          case(aatype_log)
            read(word,*,err=101) list_log(idim)
          case(aatype_int)
            read(word,*,err=102) list_int(idim)
          case(aatype_rl8)
            read(word,*,err=103) list_rl8(idim)
          case(aatype_str)
            ! use word directly
          end select

          if (idim.lt.ndim) then
            if (.not.advance_word_list_entry(wlist,' '))
     &        call quit(1,i_am,
     &                  'unexpected end of list when re-scanning')
          end if

        end do

        ! register
        select case(arg_type)
        case(aatype_label)
          call set_arg(cur_target,command,label,ndim,tgt_info,
     &         val_label=list_label)
          deallocate(list_label)
        case(aatype_log)
          call set_arg(cur_target,command,label,ndim,tgt_info,
     &         val_log=list_log)
          deallocate(list_log)
        case(aatype_int)
          call set_arg(cur_target,command,label,ndim,tgt_info,
     &         val_int=list_int)
          deallocate(list_int)
        case(aatype_rl8)
          call set_arg(cur_target,command,label,ndim,tgt_info,
     &         val_rl8=list_rl8)
          deallocate(list_rl8)
        case(aatype_str)
          call set_arg(cur_target,command,label,1,tgt_info,
     &         val_str=word)
        end select

        ! re-position the pointer on wlist
        if (ndim.gt.1) then
          ! move up again
          if (.not.advance_word_list_entry(wlist,'u')) then
            call quit(1,i_am,'unexpected error while ascending')
          end if
          ! and proceed by one
          if (.not.advance_word_list_entry(wlist,' ')) then
            call quit(1,i_am,'unexpected error after ascending')
          end if
          call get_word_list_entry(label,sep,wlist)
        end if

        if (sep.ne.','.and.sep.ne.')') then
          pt_handle_rules = error_exp_copa
          return
        end if
        if (.not.advance_word_list_entry(wlist,' ')) exit scan_arg

        ! skip possible line-break
        call get_word_list_entry(label,sep,wlist)
        if (len_trim(label).eq.0.and.sep.eq.'E') then
          if (.not.advance_word_list_entry(wlist,' ')) then
            pt_handle_rules = error_unexp_eob
            return
          end if
          call get_word_list_entry(label,sep,wlist)
          ! now at the end?
          if (len_trim(label).eq.0.and.sep.eq.')') exit scan_arg
        end if

      end do scan_arg

      ! check whether all required targets were set
      nargs = cmd_pnt%n_arguments
      do idx_arg = 1, nargs
        if (need_arg(idx_arg)) then
          arg_pnt => cmd_pnt%arg(idx_arg)
          write(lulog,*)'argument "',trim(arg_pnt%arg_label),'" not set'
          error = error_arg_not_set
        end if
      end do

      deallocate(need_arg)

      pt_handle_rules = error
      if (error.ne.no_error) return

      if (.not.advance_word_list_entry(wlist,'u')) then
        call quit(1,i_am,'unexpected absence of u-branch in wlist')
      end if

      if (.not.advance_word_list_entry(wlist,' ')) then
        pt_handle_rules = error_unexp_eob
        return
      end if

      pt_handle_rules = no_error

      return

  101 pt_handle_rules = error_read_log
      return
  102 pt_handle_rules = error_read_int
      return
  103 pt_handle_rules = error_read_rl8
      return

      end function

max_stars_repo_name: ak-ustutt/GeCCo-public

max_issues_repo_name: ak-ustutt/GeCCo-public

max_forks_repo_name: ak-ustutt/GeCCo-public


text:
      subroutine numlin(size,n)
      common/gflmn/sz,ks,kp,met
      sz=size
      ks=n
      kp=0
      met=666
      return
      end

max_stars_repo_name: sergev/grafor

max_issues_repo_name: sergev/grafor

max_forks_repo_name: sergev/grafor


PROGRAM MAINF90
  USE TEST_MODULE
  PRINT *,'Sum is',TEST_MODULE_FUNCTION(1., 2.)
END PROGRAM MAINF90


text:
*DECK STOUT
      SUBROUTINE STOUT (N, NELT, IA, JA, A, ISYM, SOLN, RHS, IUNIT, JOB)
C***BEGIN PROLOGUE  STOUT
C***PURPOSE  Write out SLAP Triad Format Linear System.
C            Routine to write out a SLAP Triad format matrix and right
C            hand side and solution to the system, if known.
C***LIBRARY   SLATEC (SLAP)
C***CATEGORY  N1
C***TYPE      SINGLE PRECISION (STOUT-S, DTOUT-D)
C***KEYWORDS  DIAGNOSTICS, LINEAR SYSTEM, SLAP SPARSE
C***AUTHOR  Seager, Mark K., (LLNL)
C             Lawrence Livermore National Laboratory
C             PO BOX 808, L-60
C             Livermore, CA 94550 (510) 423-3141
C             seager@llnl.gov
C***DESCRIPTION
C
C *Usage:
C     INTEGER N, NELT, IA(NELT), JA(NELT), ISYM, IUNIT, JOB
C     REAL    A(NELT), SOLN(N), RHS(N)
C
C     CALL STOUT( N, NELT, IA, JA, A, ISYM, SOLN, RHS, IUNIT, JOB )
C
C *Arguments:
C N      :IN       Integer
C         Order of the Matrix.
C NELT   :IN       Integer.
C         Number of non-zeros stored in A.
C IA     :IN       Integer IA(NELT).
C JA     :IN       Integer JA(NELT).
C A      :IN       Real A(NELT).
C         These arrays should hold the matrix A in the SLAP
C         Triad format.  See "Description", below.
C ISYM   :IN       Integer.
C         Flag to indicate symmetric storage format.
C         If ISYM=0, all non-zero entries of the matrix are stored.
C         If ISYM=1, the matrix is symmetric, and only the lower
C         triangle of the matrix is stored.
C SOLN   :IN       Real SOLN(N).
C         The solution to the linear system, if known.  This array
C         is accessed if and only if JOB is set to print it out,
C         see below.
C RHS    :IN       Real RHS(N).
C         The right hand side vector.  This array is accessed if and
C         only if JOB is set to print it out, see below.
C IUNIT  :IN       Integer.
C         Fortran logical I/O device unit number to write the matrix
C         to.  This unit must be connected in a system dependent fashion
C         to a file or the console or you will get a nasty message
C         from the Fortran I/O libraries.
C JOB    :IN       Integer.
C         Flag indicating what I/O operations to perform.
C         JOB = 0 => Print only the matrix.
C             = 1 => Print matrix and RHS.
C             = 2 => Print matrix and SOLN.
C             = 3 => Print matrix, RHS and SOLN.
C
C *Description:
C       The format for the output is as follows.  On  the first line
C       are counters and flags: N, NELT, ISYM, IRHS, ISOLN.  N, NELT
C       and ISYM are described above.  IRHS is  a flag indicating if
C       the RHS was  written out (1 is  yes, 0 is  no).  ISOLN  is a
C       flag indicating if the SOLN was written out  (1 is yes, 0 is
C       no).  The format for the fist line is: 5i10.  Then comes the
C       NELT Triad's IA(I), JA(I) and A(I), I = 1, NELT.  The format
C       for  these lines is   :  1X,I5,1X,I5,1X,E16.7.   Then  comes
C       RHS(I), I = 1, N, if IRHS = 1.  Then  comes SOLN(I), I  = 1,
C       N, if ISOLN = 1.  The format for these lines is: 1X,E16.7.
C
C       =================== S L A P Triad format ===================
C       This routine requires that the  matrix A be   stored in  the
C       SLAP  Triad format.  In  this format only the non-zeros  are
C       stored.  They may appear in  *ANY* order.  The user supplies
C       three arrays of  length NELT, where  NELT is  the number  of
C       non-zeros in the matrix: (IA(NELT), JA(NELT), A(NELT)).  For
C       each non-zero the user puts the row and column index of that
C       matrix element  in the IA and  JA arrays.  The  value of the
C       non-zero  matrix  element is  placed   in  the corresponding
C       location of the A array.   This is  an  extremely  easy data
C       structure to generate.  On  the  other hand it   is  not too
C       efficient on vector computers for  the iterative solution of
C       linear systems.  Hence,   SLAP changes   this  input    data
C       structure to the SLAP Column format  for  the iteration (but
C       does not change it back).
C
C       Here is an example of the  SLAP Triad   storage format for a
C       5x5 Matrix.  Recall that the entries may appear in any order.
C
C           5x5 Matrix      SLAP Triad format for 5x5 matrix on left.
C                              1  2  3  4  5  6  7  8  9 10 11
C       |11 12  0  0 15|   A: 51 12 11 33 15 53 55 22 35 44 21
C       |21 22  0  0  0|  IA:  5  1  1  3  1  5  5  2  3  4  2
C       | 0  0 33  0 35|  JA:  1  2  1  3  5  3  5  2  5  4  1
C       | 0  0  0 44  0|
C       |51  0 53  0 55|
C
C *Cautions:
C     This routine will attempt to write to the Fortran logical output
C     unit IUNIT, if IUNIT .ne. 0.  Thus, the user must make sure that
C     this logical unit is attached to a file or terminal before calling
C     this routine with a non-zero value for IUNIT.  This routine does
C     not check for the validity of a non-zero IUNIT unit number.
C***REFERENCES  (NONE)
C***ROUTINES CALLED  (NONE)
C***REVISION HISTORY  (YYMMDD)
C   871119  DATE WRITTEN
C   881213  Previous REVISION DATE
C   890915  Made changes requested at July 1989 CML Meeting.  (MKS)
C   890922  Numerous changes to prologue to make closer to SLATEC
C           standard.  (FNF)
C   890929  Numerous changes to reduce SP/DP differences.  (FNF)
C   910411  Prologue converted to Version 4.0 format.  (BAB)
C   920511  Added complete declaration section.  (WRB)
C   930701  Updated CATEGORY section.  (FNF, WRB)
C***END PROLOGUE  STOUT
C     .. Scalar Arguments ..
      INTEGER ISYM, IUNIT, JOB, N, NELT
C     .. Array Arguments ..
      REAL A(NELT), RHS(N), SOLN(N)
      INTEGER IA(NELT), JA(NELT)
C     .. Local Scalars ..
      INTEGER I, IRHS, ISOLN
C***FIRST EXECUTABLE STATEMENT  STOUT
C
C         If RHS and SOLN are to be printed also.
C         Write out the information heading.
C
      IRHS = 0
      ISOLN = 0
      IF( JOB.EQ.1 .OR. JOB.EQ.3 ) IRHS = 1
      IF( JOB.GT.1 ) ISOLN = 1
      WRITE(IUNIT,1000) N, NELT, ISYM, IRHS, ISOLN
C
C         Write out the matrix non-zeros in Triad format.
      DO 10 I = 1, NELT
         WRITE(IUNIT,1010) IA(I), JA(I), A(I)
 10   CONTINUE
C
C         If requested, write out the rhs.
      IF( IRHS.EQ.1 ) THEN
         WRITE(IUNIT,1020) (RHS(I),I=1,N)
      ENDIF
C
C         If requested, write out the solution.
      IF( ISOLN.EQ.1 ) THEN
         WRITE(IUNIT,1020) (SOLN(I),I=1,N)
      ENDIF
      RETURN
 1000 FORMAT(5I10)
 1010 FORMAT(1X,I5,1X,I5,1X,E16.7)
 1020 FORMAT(1X,E16.7)
C------------- LAST LINE OF STOUT FOLLOWS ----------------------------
      END


text:
      subroutine token(x,y,mrk,fz,n,h,th)
      call tmf(x,y,xf,yf)
      call move(xf,yf,0)
      call marker(mrk)
      call number(xf,yf,h,fz,n,th)
      return
      end


text:
!> Random Number generator from Numerical Recipes
!! To be used for reproducibility of the results
module f_random
  use f_precisions
  implicit none

  private

  integer, parameter :: ntab=32

  logical :: start = .true.
  integer :: iy = 0
  integer :: idum = 0
  integer, dimension(NTAB), save :: iv

!  public :: builtin_rand
  interface f_random_number
     module procedure f_random_number_i0,f_random_number_d0,f_random_number_r0,f_random_number_d3
     module procedure f_random_number_i1,f_random_number_d1,f_random_number_i1_i1
     !module procedure f_random_number_d4
  end interface f_random_number

  public :: f_random_number,f_random_seed

  contains

    subroutine f_random_seed(seed)
      implicit none
      integer, intent(in) :: seed
      idum=seed
    end subroutine f_random_seed

    subroutine f_random_number_d0(harvest,seed,reset)
      implicit none
      real(f_double), intent(out) :: harvest
      integer, intent(in), optional :: seed
      logical, intent(in), optional :: reset

      if (present(seed)) idum=seed
      harvest = real(builtin_rand(idum,reset=reset),f_double)

    end subroutine f_random_number_d0

    subroutine f_random_number_r0(harvest,seed,reset)
      implicit none
      real(f_simple), intent(out) :: harvest
      integer, intent(in), optional :: seed
      logical, intent(in), optional :: reset

      if (present(seed)) idum=seed
      harvest = builtin_rand(idum,reset=reset)

    end subroutine f_random_number_r0

    subroutine f_random_number_d1(harvest,seed,reset)
      implicit none
      real(f_double), dimension(:), intent(out) :: harvest
      integer, intent(in), optional :: seed
      logical, intent(in), optional :: reset
      !local variables
      integer :: i

      if (present(seed)) idum=seed

      do i=1,size(harvest)
         harvest(i) = real(builtin_rand(idum,reset=reset),f_double)
      end do
    end subroutine f_random_number_d1

    subroutine f_random_number_d3(harvest,seed,reset)
      implicit none
      real(f_double), dimension(:,:,:), intent(out) :: harvest
      integer, intent(in), optional :: seed
      logical, intent(in), optional :: reset
      !local variables
      integer :: i1,i2,i3

      if (present(seed)) idum=seed

      do i3=lbound(harvest,3),ubound(harvest,3)
         do i2=lbound(harvest,2),ubound(harvest,2)
            do i1=lbound(harvest,1),ubound(harvest,1)
               harvest(i1,i2,i3) = real(builtin_rand(idum,reset=reset),f_double)
            end do
         end do
      end do
    end subroutine f_random_number_d3

    subroutine f_random_number_i0(harvest,range,seed,reset)
      implicit none
      integer, intent(in) :: range
      integer, intent(out) :: harvest
      integer, intent(in), optional :: seed
      logical, intent(in), optional :: reset

      if (present(seed)) idum=seed
      harvest=nint(builtin_rand(idum,reset=reset)*real(range,f_simple))

    end subroutine f_random_number_i0

    subroutine f_random_number_i1(harvest,range,seed,reset)
      implicit none
      integer, intent(in) :: range
      integer, dimension(:), intent(out) :: harvest
      integer, intent(in), optional :: seed
      logical, intent(in), optional :: reset
      !local variables
      integer :: i

      if (present(seed)) idum=seed

      do i=1,size(harvest)
         harvest(i) =nint(builtin_rand(idum,reset=reset)*real(range,f_simple))
      end do
    end subroutine f_random_number_i1

    subroutine f_random_number_i1_i1(harvest,ranges,seed,reset)
      implicit none
      integer, dimension(:), intent(in) :: ranges
      integer, dimension(:), intent(out) :: harvest
      integer, intent(in), optional :: seed
      logical, intent(in), optional :: reset
      !local variables
      integer :: i

      if (present(seed)) idum=seed

      do i=1,size(harvest)
         harvest(i) =nint(builtin_rand(idum,reset=reset)*real(ranges(i),f_simple))
      end do
    end subroutine f_random_number_i1_i1

    function builtin_rand(idum, reset)
      !use randomData, only : ntab, iy, iv, start

      implicit none

      integer, intent(inout) :: idum
      real(f_simple) :: builtin_rand
      logical,intent(in),optional :: reset
      !local variables
      integer, parameter :: ia=16807,im=2147483647,iq=127773,ir=2836,ndiv=1+(im-1)/ntab
      real(f_simple), parameter :: am=1.e0/real(im,kind=f_simple),eps=1.2e-7,rnmx=1.0-eps
      integer :: j,k
      logical :: reset_

      reset_ = .false.
      if (present(reset)) reset_ = reset

      if (reset_) start = .true.

      if (start) then
         iv(:) = 0
         start = .false.
      end if
      if (idum <= 0.or. iy == 0) then
         idum=max(-idum,1)
         do j=ntab+8,1,-1
            k=idum/iq
            idum=ia*(idum-k*iq)-ir*k
            if (idum < 0) idum=idum+im
            if (j <= ntab) iv(j)=idum
         end do
         iy=iv(1)
      endif
      k=idum/iq
      idum=ia*(idum-k*iq)-ir*k
      if (idum <= 0) idum=idum+im
      j=1+iy/ndiv
      iy=iv(j)
      iv(j)=idum
      builtin_rand=min(am*iy,rnmx)
    END FUNCTION builtin_rand

end module f_random

max_stars_repo_name: adegomme/aiida-bigdft-plugin


text:
	SUBROUTINE AF_PTOP  ( istlyr, ietlyr, iptr, htop, iret )
C************************************************************************
C* AF_PTOP								*
C*									*
C* This subroutine decodes the top height of an element within a PIREP  *
C* report.								*
C*									*
C* AF_PTOP ( ISTLYR, IETLYR, IPTR, HTOP, IRET ) 	                *
C*									*
C* Input parameters:							*
C*	ISTLYR		INTEGER		Index of "like-type" group which*
C*					contains start of layer		*
C*	IETLYR		INTEGER		Index of "like-type" group which*
C*					contains end of layer		*
C*									*
C* Output parameters:							*
C*	IPTR		INTEGER		Pointer to next element         *
C*	HTOP		REAL		Top of element in feet    	*
C*	IRET		INTEGER		Return code 			*
C*					  0 = normal return 		*
C*									*
C**									*
C* Log:									*
C* D. Kidwell/NCEP	 8/99                                           *
C* J. Wu/GSC             7/00   Moved INCLUDE 'ERMISS.FNC' before the   *
C*                              DATA statement                          *
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
	INCLUDE		'afcmn.cmn'
C*
	PARAMETER	( NTIDS = 2 )
	CHARACTER	tids ( NTIDS )*3
C*
	LOGICAL		tophgt
C*
	INCLUDE		'ERMISS.FNC'
	INCLUDE		'affnc.fnc'
C*
	DATA		tids
     +			/ 'TOP', 'TPS' /
C*
C-----------------------------------------------------------------------
C
C*	Initialize variables.
C
	iret  = 0
	ifptr = istlyr
	htop  = RMISSD
C
C*	Check for a top height.
C
	tophgt = .false.
	IF ( ifptr .le. ietlyr ) THEN
	    IF ( ( lensf ( ifptr ) .eq. 1 ) .and.
     +	         ( fields ( ifptr ) (1:1) .eq. '-' ) ) THEN
	        ifptr  = ifptr + 1
		tophgt = .true.
	    END IF
C
C*	    Check at most the next four "like-type" groups for a TOP
C*	    keyword.
C
	    ii    = ifptr
	    maxii = IEDX ( ii, 3, ietlyr )
	    it    = 0
	    DO WHILE ( ( ii .le. maxii ) .and. ( it .eq. 0 ) )
		CALL ST_FIND ( fields (ii) (1:3), tids, NTIDS, it, ier )
		ii = ii + 1
		IF ( it .gt. 0 ) THEN
		    ifptr  = ii
		    tophgt = .true.
C
C*		    A TOP keyword was found.  Skip the next group if
C*		    it is not numeric.
C
		    IF ( itypsf ( ifptr ) .ne. NMR ) ifptr = ifptr + 1
		END IF
	    END DO
C
	    IF ( tophgt .and. ( ifptr .le. ietlyr ) ) THEN
		IF ( itypsf ( ifptr ) .eq. NMR ) THEN
C
C*		    Get the top height.
C
		    CALL AF_HHFM ( fields ( ifptr ) ( :lensf (ifptr) ),
     +				   htop, ier )
C
C*		    Check to see if this group is followed by a hyphen
C*		    and a second top height.
C
		    IF ( ( ifptr + 2 ) .le. ietlyr ) THEN
		        IF ( ( fields ( ifptr + 1 ) ( 1:1 ) .eq. '-' )
     +			  .and. ( itypsf ( ifptr + 2 ) .eq. NMR ) ) THEN
			    ifptr = ifptr + 2
			    CALL AF_HHFM ( fields ( ifptr )
     +			         ( :lensf ( ifptr ) ), htop2, ier )
			    IF ( .not. ERMISS ( htop2 ) ) htop = htop2
			END IF
		    END IF
		    ifptr = ifptr + 1
		END IF
	      ELSE IF ( tophgt ) THEN
C
C*		The TOP keyword or hyphen ended the layer.  Do not
C*		mistake a second base height for a top height.
C
		htop = RMISSD
	    END IF
	END IF
	iptr = ifptr
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
*deck dsmit
      subroutine dsmit(n1,c,s,v,y)
c
c  return an n1 by n1 triangular matrix c such that
c  c(t)sc=i, where i is an n1 by n1 identity matrix.
c
*mdc*if real16 vax ibm cray
*      implicit real*16 (a-h,o-z)
*mdc*else
      implicit real*8 (a-h,o-z)
*mdc*endif
*mdc*if real16q ibm vax
*c *** this mdc block is not necessary since the constants will be ***
*c *** correctly converted at compile time.  should it be removed? ***
*      parameter (a0=0.0q0, a1=1.0q0)
*mdc*else
      parameter (a0=0.0d0, a1=1.0d0)
*mdc*endif
c     common /ffm/ x,nn,n1,k1
      dimension c(n1,*), s(n1,*), v(*), y(*)
c
      do 20 i=1,n1
        do 10 j=1,i
          c(i,j)=a0
   10   continue
   20 continue
      do 80 j=1,n1
        kmax=j-1
        fac=s(j,j)
        do 30 k=1,kmax
          v(k)=a0
          y(k)=s(k,j)
   30   continue
        do 60 k=1,kmax
          dot=a0
          do 40 i=1,k
            dot=c(i,k)*y(i)+dot
   40     continue
          do 50 i=1,k
            v(i)=v(i)-dot*c(i,k)
   50     continue
          fac = fac - dot * dot
   60   continue
        fac=a1/sqrt(fac)
        c(j,j)=fac
        do 70 k=1,kmax
          c(k,j)=fac*v(k)
   70   continue
   80 continue
      return
      end

max_stars_repo_name: MOLFDIR/MOLFDIR

max_issues_repo_name: yingxingcheng/MOLFDIR

max_forks_repo_name: yingxingcheng/MOLFDIR