text:
	SUBROUTINE SF_OFIL  ( filnam, wrtflg, shrflg, isffln, iflsrc,
     +			      nparm,  parms,  iret )
C************************************************************************
C* SF_OFIL								*
C*									*
C* This subroutine opens an existing surface data file.  It is		*
C* called internally by the user-callable open routines.		*
C*									*
C* SF_OFIL  ( FILNAM, WRTFLG, SHRFLG, ISFFLN, IFLSRC, NPARM, PARMS,	*
C*            IRET )							*
C*									*
C* Input parameters:							*
C*	FILNAM		CHAR*		Surface file name		*
C*	WRTFLG		LOGICAL		Write access flag		*
C*	SHRFLG		LOGICAL		Share flag			*
C*									*
C* Output parameters:							*
C*	ISFFLN		INTEGER		File number			*
C*	IFLSRC		INTEGER		Data source			*
C*	NPARM		INTEGER		Number of parameters		*
C*	PARMS (NPARM)	CHAR*4		Parameter names			*
C*	IRET		INTEGER		Return code			*
C*					  0 = normal return		*
C*					 -2 = file could not be opened	*
C*					 -6 = file not surface file	*
C**									*
C* Log: 								*
C* I. Graffman/RDS	 8/87						*
C* M. desJardins/GSFC	 6/88	Documentation				*
C* K. Brill/NMC		 8/93	Changes for STD2 and SPRI		*
C* D. Keiser/GSC	 3/96	Added choice between SFTX and SFDT	*
C* T. Lee/GSC		10/97	Added KSWFO and KWFO2 in DM_LSSF	*
C*				calling sequences			*
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
	INCLUDE		'GMBDTA.CMN'
	INCLUDE		'sfcmn.cmn'
C*
	CHARACTER*(*)	filnam, parms (*)
	LOGICAL		wrtflg, shrflg
C*
	INTEGER		iscale (MMPARM), iofset (MMPARM), ibit (MMPARM)
C------------------------------------------------------------------------
C*	Open the file.
C
	CALL DM_OPEN  ( filnam, wrtflg, shrflg, isffln, iftyp, iflsrc,
     +			nrow,   ncol,   nprt,   nfhdrs, iret )
	IF  ( iret .ne. 0 )  THEN
	    iret = -2
	    CALL ER_WMSG  ( 'SF', iret, filnam, ier )
	    RETURN
	END IF
C
C*	Set file open flag.
C
	isfcfn (isffln) = isffln
C
C*	Check that this is a surface file.
C
	IF  ( iftyp .ne. MFSF )  THEN
	    iret = -6
	    CALL ER_WMSG  ( 'SF', iret, filnam, ier )
	    CALL SF_CLOS  ( isffln, ier )
	    RETURN
	END IF
C
C*	Get parameter information.
C
	IF ( ( iflsrc .eq. 100 ) .and. ( nprt .eq. 1 ) ) THEN
	    CALL DM_PART  ( isffln, 'SFTX', lenh, it, nparm, parms,
     +			    iscale, iofset, ibit, iret )
	ELSE
	    CALL DM_PART  ( isffln, 'SFDT', lenh, it, nparm, parms,
     +			    iscale, iofset, ibit, iret )
	END IF
	IF  ( iret .ne. 0 )  THEN
	    iret = -6
	    CALL ER_WMSG  ( 'SF', iret, filnam, ier )
	    CALL SF_CLOS  ( isffln, ier )
	    RETURN
	END IF
C
C*	Retrieve date/time keys.
C
	CALL DM_LTIM  ( isffln, dttype (isffln), kdate (isffln),
     +			ktime (isffln), iret )
	IF  ( iret .ne. 0 )  THEN
	    CALL ER_WMSG  ( 'DM', iret, filnam, ier )
	    iret = -6
	    CALL ER_WMSG  ( 'SF', iret, filnam, ier )
	    CALL SF_CLOS  ( isffln, ier )
	    RETURN
	END IF
C
C*	Retrieve station information.
C
	CALL DM_LSSF  ( isffln, sttype (isffln), kstid (isffln),
     +			kstnm (isffln), kslat (isffln),
     +			kslon (isffln), kselv (isffln),
     +			kstat (isffln), kcoun (isffln),
     +			kstd2 (isffln), kspri (isffln),
     +			kswfo (isffln), kwfo2 (isffln),
     +			iret )
	IF  ( iret .ne. 0 )  THEN
	    CALL ER_WMSG  ( 'DM', iret, filnam, ier )
	    iret = -6
	    CALL SF_CLOS  ( isffln, ier )
	    RETURN
	END IF
C
C*	Save number of parameters, reset row and column.
C
	kparm  (isffln) = nparm
	krow   (isffln) = 0
	kcol   (isffln) = 0
	timset (isffln) = .false.
	stnset (isffln) = .false.
	ftmset (isffln) = .false.
	curtim (isffln) = ' '
	curstn (isffln) = ' '
C*
	RETURN
	END

ext: f

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
MODULE VCA_VARS_GLOBAL
  USE VCA_INPUT_VARS
  USE VCA_SPARSE_MATRIX
  !
  USE SF_CONSTANTS
  USE SF_ARRAYS,  only: arange,linspace
  USE SF_IOTOOLS, only: str
#ifdef _MPI
  USE MPI
  USE SF_MPI
#endif
  implicit none

  !-------------------- CUSTOM OBSERVABLE STRUCTURE ----------------------!
  type observable
    complex(8),dimension(:,:,:,:,:,:,:),allocatable :: sij
    character(len=32)                               :: o_name
    real(8)                                         :: o_value
  end type observable

  type custom_observables
     type(observable),dimension(:),allocatable               :: item
     integer                                                 :: N_asked
     integer                                                 :: N_filled
     logical                                                 :: init=.false.
  end type custom_observables

  !-------------------- EFFECTIVE BATH STRUCTURE ----------------------!
  type effective_bath
     real(8),dimension(:,:,:,:),allocatable        :: e     !local energies [1//Nlat][Nspin][1//Norb][Nbath]
     real(8),dimension(:,:,:,:),allocatable        :: v     !spin-keep hyb. [Nlat][Nspin][Norb][Nbath]
     logical                                       :: status=.false.
  end type effective_bath

  !------------------ FULL HAMILTONIAN STRUCTURE ---------------------!
  type full_espace
     real(8),dimension(:),pointer      :: e
     complex(8),dimension(:,:),pointer :: M
  end type full_espace

  !---------------- SECTOR-TO-FOCK SPACE STRUCTURE -------------------!
  type sector_map
     integer,dimension(:),allocatable :: map
     logical                          :: status=.false.
  end type sector_map

  interface map_allocate
     module procedure :: map_allocate_scalar
     module procedure :: map_allocate_vector
  end interface map_allocate

  interface map_deallocate
     module procedure :: map_deallocate_scalar
     module procedure :: map_deallocate_vector
  end interface map_deallocate

  !-------------- GMATRIX FOR FAST EVALUATION OF GF ------------------!
  !note that we use a single Qmatrix here which must be intended as
  !component corresponding to the poles.
  type GFspectrum
     complex(8),dimension(:),allocatable :: weight
     complex(8),dimension(:),allocatable :: poles
  end type GFspectrum
  type GFchannel
     type(GFspectrum),dimension(:),allocatable :: channel !N_channel = 2 (c,cdag), 4 (c,cdag,c pm cdag)
  end type GFchannel
 type GFmatrix
     type(GFchannel),dimension(:),allocatable :: state !state_list%size = # of state in the spectrum
  end type GFmatrix

  interface GFmatrix_allocate
     module procedure :: allocate_GFmatrix_Nstate
     module procedure :: allocate_GFmatrix_Nchan
     module procedure :: allocate_GFmatrix_Nexc
  end interface GFmatrix_allocate

  !------------------ ABTRACT INTERFACES PROCEDURES ------------------!
  !SPARSE MATRIX-VECTOR PRODUCTS USED IN ED_MATVEC
  !dbleMat*dbleVec
  abstract interface
     subroutine cc_sparse_HxV(Nloc,v,Hv)
       integer                    :: Nloc
       complex(8),dimension(Nloc) :: v
       complex(8),dimension(Nloc) :: Hv
     end subroutine cc_sparse_HxV
  end interface

 !-------------------------- ED  VARIABLES --------------------------!

  !SIZE OF THE PROBLEM
  !=========================================================
  integer,save                                       :: Ns       !Number of levels per spin
  integer,save                                       :: Nsectors !Number of sectors
  integer,save                                       :: Ns_orb
  integer,save                                       :: Ns_ud

  !non-interacting cluster Hamiltonian and full system hopping matrix
  !=========================================================
  complex(8),dimension(:,:,:,:,:,:),allocatable     :: impHloc ![Nlat][Nlat][Nspin][Nspin][Norb][Norb]
  complex(8),dimension(:,:,:,:,:,:,:),allocatable   :: impHk   ![Nlat][Nlat][Nspin][Nspin][Norb][Norb][Nktot]
  complex(8),dimension(:,:),allocatable             :: embeddedHloc ![Nlat*Nspin*Norb*(1+Nbath)][Nlat*Nspin*Norb*(1+Nbath)]
  complex(8),dimension(:,:,:),allocatable           :: embeddedHk ![Nlat*Nspin*Norb*(1+Nbath)][Nlat*Nspin*Norb*(1+Nbath)][Nktot]

  !Some maps between sectors and full Hilbert space (pointers)  CHECK!!!!!!!!!!!!!!!!!!!!!
  !=========================================================
  integer,allocatable,dimension(:,:)              :: getsector
  integer,allocatable,dimension(:,:,:)            :: getCsector
  integer,allocatable,dimension(:,:,:)            :: getCDGsector
  integer,allocatable,dimension(:)                :: getDim,getDimUp,getDimDw
  integer,allocatable,dimension(:)                :: getNup,getNdw
  integer,allocatable,dimension(:,:,:)            :: getBathStride
  integer,allocatable,dimension(:,:)              :: impIndex
  logical,allocatable,dimension(:)                :: twin_mask
  logical,allocatable,dimension(:)                :: sectors_mask

  !Effective Bath used in the VCA code (this is opaque to user)
  !=========================================================
  type(effective_bath)                            :: vca_bath

  !Eigenvalues,Eigenvectors FULL DIAGONALIZATION
  !=========================================================
  type(full_espace),dimension(:),allocatable        :: espace

  !Sparse matrix for Lanczos diagonalization.
  !=========================================================
  type(sparse_matrix_csr)                            :: spH0d !diagonal part
  type(sparse_matrix_csr)                            :: spH0nd !non-diagonal part
  type(sparse_matrix_csr),dimension(:),allocatable   :: spH0ups,spH0dws !reduced UP and DW parts
  procedure(cc_sparse_HxV),pointer                   :: spHtimesV_p=>null()

  !Variables for DIAGONALIZATION
  !=========================================================
  integer,allocatable,dimension(:)                  :: neigen_sector
  !--------------- LATTICE WRAP VARIABLES -----------------!
  integer,allocatable,dimension(:,:)                 :: neigen_sectorii
  integer,allocatable,dimension(:)                   :: neigen_totalii
  logical                                            :: trim_state_list=.false.

  !Partition function, Omega potential, SFT potential,original lattice bandwidth
  !=========================================================
  real(8)                                           :: zeta_function
  real(8)                                           :: omega_potential
  real(8)                                           :: sft_potential

  !Cluster Green's functions
  !(Nlat,Nlat,Nspin,Nspin,Norb,Norb,:)
  !=========================================================
  complex(8),allocatable,dimension(:,:,:,:,:,:,:)   :: impGmats ![Nlat][Nlat][Nspin][Nspin][Norb][Norb][L]
  complex(8),allocatable,dimension(:,:,:,:,:,:,:)   :: impGreal ![Nlat][Nlat][Nspin][Nspin][Norb][Norb][L]
  !
  complex(8),allocatable,dimension(:,:,:,:,:,:,:)   :: impG0mats ![Nlat][Nlat][Nspin][Nspin][Norb][Norb][L]
  complex(8),allocatable,dimension(:,:,:,:,:,:,:)   :: impG0real ![Nlat][Nlat][Nspin][Nspin][Norb][Norb][L]
  !
  complex(8),allocatable,dimension(:,:,:,:,:,:,:)   :: impSmats ![Nlat][Nlat][Nspin][Nspin][Norb][Norb][L]
  complex(8),allocatable,dimension(:,:,:,:,:,:,:)   :: impSreal ![Nlat][Nlat][Nspin][Nspin][Norb][Norb][L]
  !
  type(GFmatrix),allocatable,dimension(:,:,:,:,:,:) :: impGmatrix

  !Spin Susceptibilities
  !=========================================================
  real(8),allocatable,dimension(:,:)                 :: spinChi_tau
  complex(8),allocatable,dimension(:,:)              :: spinChi_w
  complex(8),allocatable,dimension(:,:)              :: spinChi_iv

  !Cluster local observables:
  !=========================================================
  real(8),dimension(:,:),allocatable                ::  imp_dens    ![Nlat][Norb]
  real(8),dimension(:,:),allocatable                ::  imp_dens_up ![Nlat][Norb]
  real(8),dimension(:,:),allocatable                ::  imp_dens_dw ![Nlat][Norb]
  real(8),dimension(:,:),allocatable                ::  imp_docc    ![Nlat][Norb]

  !Custom observables:
  !=========================================================
  type(custom_observables)                             ::  custom_o

  !Suffix string attached to the output files.
  !=========================================================
  character(len=64)                                 :: file_suffix=""
  logical                                           :: offdiag_gf_flag=.false.

  !This is the internal Mpi Communicator and variables.
  !=========================================================
#ifdef _MPI
  integer                                            :: MpiComm_Global=MPI_COMM_NULL
  integer                                            :: MpiComm=MPI_COMM_NULL
#endif
  integer                                            :: MpiGroup_Global=MPI_GROUP_NULL
  integer                                            :: MpiGroup=MPI_GROUP_NULL
  logical                                            :: MpiStatus=.false.
  logical                                            :: MpiMaster=.true.
  integer                                            :: MpiRank=0
  integer                                            :: MpiSize=1
  integer,allocatable,dimension(:)                   :: MpiMembers
  integer                                            :: mpiQup=0
  integer                                            :: mpiRup=0
  integer                                            :: mpiQdw=0
  integer                                            :: mpiRdw=0
  integer                                            :: mpiQ=0
  integer                                            :: mpiR=0
  integer                                            :: mpiIstart
  integer                                            :: mpiIend
  integer                                            :: mpiIshift
  logical                                            :: mpiAllThreads=.true.

  !Frequency and time arrays:
  !=========================================================
  !real(8),dimension(:),allocatable                  :: wm,tau,wr,vm

  !flag for finite temperature calculation
  !=========================================================
  logical                                           :: finiteT

contains

  !> Get stride position in the one-particle many-body space
  function index_stride_lso(ilat,ispin,iorb) result(indx)
    integer :: ilat
    integer :: iorb
    integer :: ispin
    integer :: indx
    indx = iorb + (ilat-1)*Norb + (ispin-1)*Norb*Nlat
  end function index_stride_lso

  !>Allocate and Deallocate Hilbert space maps (sector<-->Fock)
  subroutine map_allocate_scalar(H,N)
    type(sector_map) :: H
    integer          :: N
    if(H%status) call map_deallocate_scalar(H)
    allocate(H%map(N))
    H%status=.true.
  end subroutine map_allocate_scalar
  !
  subroutine map_allocate_vector(H,N)
    type(sector_map),dimension(:)       :: H
    integer,dimension(size(H))          :: N
    integer                             :: i
    do i=1,size(H)
       call map_allocate_scalar(H(i),N(i))
    enddo
  end subroutine map_allocate_vector

  subroutine map_deallocate_scalar(H)
    type(sector_map) :: H
    if(.not.H%status)then
       write(*,*) "WARNING map_deallocate_scalar: H is not allocated"
       return
    endif
    if(allocated(H%map))deallocate(H%map)
    H%status=.false.
  end subroutine map_deallocate_scalar
  !
  subroutine map_deallocate_vector(H)
    type(sector_map),dimension(:) :: H
    integer                       :: i
    do i=1,size(H)
       call map_deallocate_scalar(H(i))
    enddo
  end subroutine map_deallocate_vector

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

 subroutine vca_set_MpiComm(comm)
#ifdef _MPI
    integer :: comm,ierr
    MpiComm_Global = comm
    MpiComm        = MpiComm_Global
    MpiStatus      = .true.
    MpiSize        = get_Size_MPI(MpiComm_Global)
    MpiRank        = get_Rank_MPI(MpiComm_Global)
    MpiMaster      = get_Master_MPI(MpiComm_Global)
    call Mpi_Comm_group(MpiComm_Global,MpiGroup_Global,ierr)
#else
    integer,optional :: comm
#endif
  end subroutine vca_set_MpiComm

  subroutine vca_del_MpiComm()
#ifdef _MPI
    MpiComm_Global = MPI_UNDEFINED
    MpiComm        = MPI_UNDEFINED
    MpiStatus      = .false.
    MpiSize        = 1
    MpiRank        = 0
    MpiMaster      = .true.
#endif
  end subroutine vca_del_MpiComm

  !Allocate the channels in GFmatrix structure
subroutine allocate_gfmatrix_Nstate(self,Nstate)
    type(GFmatrix) :: self
    integer        :: Nstate
    if(allocated(self%state))deallocate(self%state)
    allocate(self%state(Nstate))
  end subroutine allocate_gfmatrix_Nstate

  subroutine allocate_gfmatrix_Nchan(self,istate,Nchan)
    type(GFmatrix) :: self
    integer        :: istate,Nchan
    if(allocated(self%state(istate)%channel))deallocate(self%state(istate)%channel)
    allocate(self%state(istate)%channel(Nchan))
  end subroutine allocate_gfmatrix_Nchan

  !Allocate the Excitations spectrum at a given channel
  subroutine allocate_gfmatrix_Nexc(self,istate,ichan,Nexc)
    type(GFmatrix) :: self
    integer        :: istate,ichan
    integer        :: Nexc
    if(allocated(self%state(istate)%channel(ichan)%weight))deallocate(self%state(istate)%channel(ichan)%weight)
    if(allocated(self%state(istate)%channel(ichan)%poles))deallocate(self%state(istate)%channel(ichan)%poles)
    allocate(self%state(istate)%channel(ichan)%weight(Nexc))
    allocate(self%state(istate)%channel(ichan)%poles(Nexc))
  end subroutine allocate_gfmatrix_Nexc

END MODULE VCA_VARS_GLOBAL

max_stars_repo_name: QcmPlab/VCA

max_issues_repo_name: QcmPlab/VCA

max_forks_repo_name: QcmPlab/VCA


!-------------------------------------------------------------------------------
! Copyright (c) 2016 The University of Tokyo
! This software is released under the MIT License, see LICENSE.txt
!-------------------------------------------------------------------------------
!> This module provide a function to ECHO for IFSTR solver
module m_static_echo
contains
  !C
  !C***
  !> ECHO for IFSTR solver
  !C***
  !C
  subroutine FSTR_ECHO (hecMESH)
    use m_fstr
    use m_hecmw2fstr_mesh_conv
    implicit none
    type (hecmwST_local_mesh):: hecMESH

    !** Local variables
    integer(kind=kint) :: i,j,is,iE,nn,ic_type
    integer(kind=kint) :: icel,isect,ig1,iS0,iE0,ik
    integer(kind=kint) :: in,nid, itype
    real(kind=kreal) :: x,y,z
    integer(kind=kint):: nids(20)
    !C====
    !C +-------------------------------+
    !C | NODE                          |
    !C +-------------------------------+
    !C===
    !C** nodal coordinate
    write(ILOG,*) '### Number of nodes',hecMESH%n_node
    write(ILOG,*) 'ID X Y Z'
    do i=1,hecMESH%n_node
      nid = hecMESH%global_node_ID(i)
      x = hecMESH%node(3*i-2)
      y = hecMESH%node(3*i-1)
      z = hecMESH%node(3*i)
      write(ILOG,'(i8,e15.5,e15.5,e15.5)') nid,x,y,z
    enddo
    !C
    !C +-------------------------------+
    !C | ELEMENT                       |
    !C +-------------------------------+
    !C===
    call fstr2hecmw_mesh_conv( hecMESH )
    write(ILOG,*) '### Elements', hecMESH%n_elem
    do itype= 1, hecMESH%n_elem_type
      is= hecMESH%elem_type_index(itype-1) + 1
      iE= hecMESH%elem_type_index(itype  )
      ic_type= hecMESH%elem_type_item(itype)
      !C** Set number of nodes
      nn = hecmw_get_max_node(ic_type)
      !C** element loop
      do icel= is, iE
        !C** node ID
        is= hecMESH%elem_node_index(icel-1)
        do j=1,nn
          if( hecMESH%n_refine > 0 ) then
            nids(j)= hecMESH%elem_node_item (is+j)
          else
            nids(j)= hecMESH%global_node_ID( hecMESH%elem_node_item (is+j))
          endif
        enddo
        !C** section  ID
        isect= hecMESH%section_ID(icel)
        write(ILOG,*) '### Element ID=',ic_type,isect,hecMESH%global_elem_id(icel)
        write(ILOG,*) (nids(j),j=1,nn)
      enddo
    enddo
    call hecmw2fstr_mesh_conv( hecMESH )
    !C +-------------------------------+
    !C | NODE GROUP                    |
    !C +-------------------------------+
    write(ILOG,*) '### Ngroup'
    do ig1= 1, hecMESH%node_group%n_grp
      write(ILOG,*)
      write(ILOG,'(a80)') hecMESH%node_group%grp_name(ig1)
      iS0= hecMESH%node_group%grp_index(ig1-1) + 1
      iE0= hecMESH%node_group%grp_index(ig1  )
      do ik= iS0, iE0
        in   = hecMESH%node_group%grp_item(ik)
        write(ILOG,*) hecMESH%global_node_ID(in)
      enddo
    enddo
    !C +-------------------------------+
    !C | ELEMEN GROUP                  |
    !C +-------------------------------+
    write(ILOG,*) '### Egroup'
    do ig1= 1, hecMESH%elem_group%n_grp
      write(ILOG,*)
      write(ILOG,'(a80)') hecMESH%elem_group%grp_name(ig1)
      iS0= hecMESH%elem_group%grp_index(ig1-1) + 1
      iE0= hecMESH%elem_group%grp_index(ig1  )
      do ik= iS0, iE0
        in   = hecMESH%elem_group%grp_item(ik)
        write(ILOG,*) hecMESH%global_elem_ID(in)
      enddo
    enddo
    write(ILOG,*) '### Reftemp',ref_temp
    !C====
  end subroutine FSTR_ECHO
end module m_static_echo


text:
      SUBROUTINE GETOKE(CBUFO,NCHARS)
c  ****   Watch out - getoke writes a null after the last char,
c         so cbufo must be bigger than the number of characters
C
C   GETOKE RETURNS CONSECUTIVE TOKENS (ITEMS BETWEEN A DELIMITER), ONE PER CALL,
C   FROM CHARACTER STRING in COMMON/sioln1/. AN ALPHA STRING IS RETURNED IN
C   TOKEN WHEN IT STARTS AND ENDS WHEN SINGLE QUOTES.  (THE QUOTES ARE NOT
C   RETURNED.  THE STRING MUST BE TERMINATED WITH A QUOTE AND A BLANK, SO THAT
C   QUOTES MAY BE INCLUDED IN THE STRING SO LONG AS THE QUOTE IS NOT FOLLOWED BY
C   A BLANK).
c       nchars will be a zero if end of line is detected or an comment character
c   is detected.  A comment character is either ! (VMS & Cray), { (Apollo),
c   # (Unix). Any of these comment characters will work on any machine!
C
C  ARGUMENTS:
C    CBUFO  - THE CHARACTER*(*) STRING SET BY GETOKE CONTAINING THE NEXT TOKEN
C             FOUND.  CBUFO MUST BE NCHARS+1 CHARACTERS LONG (SINCE C STRINGS
C             MUST BE TERMINATED WITH A NULL).
C    NCHARS - THE NUMBER OF CHARACTERS IN THE TOKEN RETURNED IN CBUFO. A 0
C             (ZERO) NUMBER OF CHARACTERS MEANS THAT NO TOKEN WAS FOUND AND
C             THAT ANOTHER LINE SHOULD BE READ AND GETOKE CALLED AGAIN WITH
C             jchar=1.
C
C
C
C  WRITTEN AND COPYRIGHTED BY:
C  PAUL HENKART, SCRIPPS INSTITUTION OF OCEANOGRAPHY, 29 FEBRUARY 1984
C  ALL RIGHTS ARE RESERVED BY THE AUTHOR.  PERMISSION TO COPY OR REPRODUCE THIS
C  SUBROUTINE, BY COMPUTER OR OTHER MEANS, MAY BE OBTAINED ONLY FROM THE AUTHOR.
c
c  mod 25 Aug 98 - Make a NULL be a delimiter
c  mod 19 Jul 02 - Add entry getoke1 to go along with rline1, using different
c                  common blocks.
c  mod 2 May 04 - Add entry getokec for getoke with a comma delimiter.
c
      CHARACTER*1 CDELIM,QUOTE,tab, null
      CHARACTER*3000 CBUFIN
      COMMON /SIOLN1/ CBUFIN
      COMMON /sioln2/ jchar,NCBUF
      CHARACTER* (*) CBUFO
      DATA CDELIM/' '/, QUOTE/''''/

      cbufo = ' '
      NCHARS=0                                                          ! COUNT THE NON BLANK CHARACTERS IN THE TOKEN
      IQUOTE=0                                                          ! COUNT THE QUOTES IN THE TOKEN
      tab = CHAR(9)                                                     ! the tab character
      null = CHAR(0)                                                    ! NULL
      IF(jchar.LT.1) jchar=1
   10 CONTINUE
      IF( CBUFIN(jchar:jchar) .NE. CDELIM .AND.
     *   cbufin(jchar:jchar) .NE. null .AND.
     *   cbufin(jchar:jchar) .NE. tab ) GOTO 20                         ! STRIP OFF LEADING BLANKS
      jchar=jchar+1
      IF(jchar.GT.NCBUF) RETURN
      GO TO 10
   20 IF(CBUFIN(jchar:jchar).NE.QUOTE) GO TO 30                         ! IS IT A QUOTE?
      ISTART=jchar+1                                                    ! STRIP OF THE LEADING QUOTE
      IQUOTE=1                                                          ! SIGNAL THAT THE STRING STARTED WITH A QUOTE
      GO TO 40
   30 CONTINUE
      IF( cbufin(jchar:jchar) .EQ. '!' .OR.
     *    cbufin(jchar:jchar) .EQ. '{' .OR.
     *    cbufin(jchar:jchar) .EQ. '#' ) THEN
             nchars = 0
             GOTO 100
      ENDIF
c**** toss out non ASCII characters
      IF( ICHAR(cbufin(jchar:jchar)) .LT. 32 ) THEN
          jchar = jchar + 1
          GOTO 10
      ENDIF
      ISTART=jchar                                                      ! THE FIRST CHARACTER OF THE TOKEN TO BE RETURNED
   40 CONTINUE                                                          ! NOW FIND THE END OF THE TOKEN
      IF(jchar.GT.NCBUF) GO TO 110                                      ! ARE WE AT THE END OF THE BUFFER?
      IF(CBUFIN(jchar:jchar).EQ.CDELIM.AND.IQUOTE.NE.1) GOTO 100        ! WAS IT A BLANK?
      IF(CBUFIN(jchar:jchar).EQ. tab .AND. IQUOTE.NE.1) GOTO 100        ! WAS IT A BLANK?
      IF(CBUFIN(jchar:jchar).EQ. null .AND. IQUOTE.NE.1) GOTO 100       ! WAS IT A NULL?
      jchar=jchar+1
      IF( CBUFIN(jchar:jchar) .EQ. QUOTE .AND. IQUOTE .EQ. 1 ) THEN
          IQUOTE = 2
          GOTO 100
      ENDIF
      NCHARS=NCHARS+1                                                   ! THE CURREN CHARACTER IS NOT A BLANK OR A QUOTE
      GO TO 40                                                          ! GO LOOK AT THE NEXT CHARACTER
  100 CONTINUE
      IF( NCHARS .EQ. 0 ) RETURN                                        ! DON'T TRY TO MOVE ZERO CHARACTERS!!
      IF(IQUOTE.ne.1) GO TO 110                                         ! THE CURRENT CHARACTER IS A BLANK WITHIN QUOTES
      jchar=jchar+1
      nchars=nchars+1
      go to 40
c**** end the returned string with a null character so that c rcognizes
c**** the end of string.
  110 IF( nchars .GT. 0 ) CBUFO(1:NCHARS)=CBUFIN(ISTART:jchar)
      cbufo(nchars+1:nchars+1) = null
      jchar=jchar+1                                                     ! STRIP OFF THE BLANK
      RETURN
      END

      SUBROUTINE GETOKE1(CBUFO,NCHARS)
c**** Use rline1 if rline is being used by someone else.  This allows
c**** rdline and rline to coexist without rline using rdline's common.
c**** Use getoke1 to get the tokens from rline1 lines (common sioln3).
C
C   GETOKE RETURNS CONSECUTIVE TOKENS (ITEMS BETWEEN A DELIMITER), ONE PER CALL,
C   FROM CHARACTER STRING in COMMON/sioln1/. AN ALPHA STRING IS RETURNED IN
C   TOKEN WHEN IT STARTS AND ENDS WHEN SINGLE QUOTES.  (THE QUOTES ARE NOT
C   RETURNED.  THE STRING MUST BE TERMINATED WITH A QUOTE AND A BLANK, SO THAT
C   QUOTES MAY BE INCLUDED IN THE STRING SO LONG AS THE QUOTE IS NOT FOLLOWED BY
C   A BLANK).
c       nchars will be a zero if end of line is detected or an comment character
c   is detected.  A comment character is either ! (VMS & Cray), { (Apollo),
c   # (Unix). Any of these comment characters will work on any machine!
C
C  ARGUMENTS:
C    CBUFO  - THE CHARACTER*(*) STRING SET BY GETOKE CONTAINING THE NEXT TOKEN
C             FOUND.  CBUFO MUST BE NCHARS+1 CHARACTERS LONG (SINCE C STRINGS
C             MUST BE TERMINATED WITH A NULL).
C    NCHARS - THE NUMBER OF CHARACTERS IN THE TOKEN RETURNED IN CBUFO. A 0
C             (ZERO) NUMBER OF CHARACTERS MEANS THAT NO TOKEN WAS FOUND AND
C             THAT ANOTHER LINE SHOULD BE READ AND GETOKE CALLED AGAIN WITH
C             jchar=1.
C
C
C
C  WRITTEN AND COPYRIGHTED BY:
C  PAUL HENKART, SCRIPPS INSTITUTION OF OCEANOGRAPHY, 29 FEBRUARY 1984
C  ALL RIGHTS ARE RESERVED BY THE AUTHOR.  PERMISSION TO COPY OR REPRODUCE THIS
C  SUBROUTINE, BY COMPUTER OR OTHER MEANS, MAY BE OBTAINED ONLY FROM THE AUTHOR.
c
c  mod 25 Aug 98 - Make a NULL be a delimiter
c
      CHARACTER*1 CDELIM,QUOTE,tab, null
      CHARACTER*200 CBUFIN
      COMMON /SIOLN3/ CBUFIN
      COMMON /sioln4/ jchar,NCBUF
      CHARACTER* (*) CBUFO
      DATA CDELIM/' '/, QUOTE/''''/

      cbufo = ' '
      NCHARS=0                                                          ! COUNT THE NON BLANK CHARACTERS IN THE TOKEN
      IQUOTE=0                                                          ! COUNT THE QUOTES IN THE TOKEN
      tab = CHAR(9)                                                     ! the tab character
      null = CHAR(0)                                                    ! NULL
      IF(jchar.LT.1) jchar=1
   10 CONTINUE
      IF( CBUFIN(jchar:jchar) .NE. CDELIM .AND.
     *   cbufin(jchar:jchar) .NE. null .AND.
     *   cbufin(jchar:jchar) .NE. tab ) GOTO 20                         ! STRIP OFF LEADING BLANKS
      jchar=jchar+1
      IF(jchar.GT.NCBUF) RETURN
      GO TO 10
   20 IF(CBUFIN(jchar:jchar).NE.QUOTE) GO TO 30                         ! IS IT A QUOTE?
      ISTART=jchar+1                                                    ! STRIP OF THE LEADING QUOTE
      IQUOTE=1                                                          ! SIGNAL THAT THE STRING STARTED WITH A QUOTE
      GO TO 40
   30 CONTINUE
      IF( cbufin(jchar:jchar) .EQ. '!' .OR.
     *    cbufin(jchar:jchar) .EQ. '{' .OR.
     *    cbufin(jchar:jchar) .EQ. '#' ) THEN
             nchars = 0
             GOTO 100
      ENDIF
c**** toss out non ASCII characters
      IF( ICHAR(cbufin(jchar:jchar)) .LT. 32 ) THEN
          jchar = jchar + 1
          GOTO 10
      ENDIF
      ISTART=jchar                                                      ! THE FIRST CHARACTER OF THE TOKEN TO BE RETURNED
   40 CONTINUE                                                          ! NOW FIND THE END OF THE TOKEN
      IF(jchar.GT.NCBUF) GO TO 110                                      ! ARE WE AT THE END OF THE BUFFER?
      IF(CBUFIN(jchar:jchar).EQ.CDELIM.AND.IQUOTE.NE.1) GOTO 100        ! WAS IT A BLANK?
      IF(CBUFIN(jchar:jchar).EQ. tab .AND. IQUOTE.NE.1) GOTO 100        ! WAS IT A BLANK?
      IF(CBUFIN(jchar:jchar).EQ. null .AND. IQUOTE.NE.1) GOTO 100       ! WAS IT A NULL?
      jchar=jchar+1
      IF( CBUFIN(jchar:jchar) .EQ. QUOTE .AND. IQUOTE .EQ. 1 ) THEN
          IQUOTE = 2
          GOTO 100
      ENDIF
      NCHARS=NCHARS+1                                                   ! THE CURREN CHARACTER IS NOT A BLANK OR A QUOTE
      GO TO 40                                                          ! GO LOOK AT THE NEXT CHARACTER
  100 CONTINUE
      IF( NCHARS .EQ. 0 ) RETURN                                        ! DON'T TRY TO MOVE ZERO CHARACTERS!!
      IF(IQUOTE.ne.1) GO TO 110                                         ! THE CURRENT CHARACTER IS A BLANK WITHIN QUOTES
      jchar=jchar+1
      nchars=nchars+1
      go to 40
c**** end the returned string with a null character so that c rcognizes
c**** the end of string.
  110 IF( nchars .GT. 0 ) CBUFO(1:NCHARS)=CBUFIN(ISTART:jchar)
      cbufo(nchars+1:nchars+1) = null
      jchar=jchar+1                                                     ! STRIP OFF THE BLANK
      RETURN
      END

      SUBROUTINE GETOKEC(CBUFO,NCHARS)
c     This is the same as GETOKE, except COMMA is also a delimiter.
c   ( I don't know what will break if I just add it to geoke, so make
c   it a unique entry point).
c  mod 25 Aug 98 - Make a NULL be a delimiter
c  mod 2 May 98 - Make a COMMA be a delimiter
c               - Don't honor the comment convention.
c mod 22 July 2007 - Look out for ,, - vacuous field
c
      CHARACTER*1 CDELIM,QUOTE,tab, null, comma
      CHARACTER*200 CBUFIN
      COMMON /SIOLN1/ CBUFIN
      COMMON /sioln2/ jchar,NCBUF
      CHARACTER* (*) CBUFO
      DATA CDELIM/' '/, QUOTE/''''/, comma/','/

      cbufo = ' '
      NCHARS=0                                                          ! COUNT THE NON BLANK CHARACTERS IN THE TOKEN
c****   look out for vacuous field ( ,, )
      IF( cbufin(jchar:jchar) .EQ. comma ) GOTO 110
      IQUOTE=0                                                          ! COUNT THE QUOTES IN THE TOKEN
      tab = CHAR(9)                                                     ! the tab character
      null = CHAR(0)                                                    ! NULL
      IF(jchar.LT.1) jchar=1
   10 CONTINUE
      IF( CBUFIN(jchar:jchar) .NE. CDELIM .AND.
     *   cbufin(jchar:jchar) .NE. null .AND.
     *   cbufin(jchar:jchar) .NE. comma .AND.
     *   cbufin(jchar:jchar) .NE. tab ) GOTO 20                         ! STRIP OFF LEADING BLANKS
      jchar=jchar+1
      IF(jchar.GT.NCBUF) RETURN
      GO TO 10
   20 IF(CBUFIN(jchar:jchar).NE.QUOTE) GO TO 30                         ! IS IT A QUOTE?
      ISTART=jchar+1                                                    ! STRIP OF THE LEADING QUOTE
      IQUOTE=1                                                          ! SIGNAL THAT THE STRING STARTED WITH A QUOTE
      GO TO 40
   30 CONTINUE
c**** toss out non ASCII characters
      IF( ICHAR(cbufin(jchar:jchar)) .LT. 32 ) THEN
          jchar = jchar + 1
          GOTO 10
      ENDIF
      ISTART=jchar                                                      ! THE FIRST CHARACTER OF THE TOKEN TO BE RETURNED
   40 CONTINUE                                                          ! NOW FIND THE END OF THE TOKEN
      IF(jchar.GT.NCBUF) GO TO 110                                      ! ARE WE AT THE END OF THE BUFFER?
      IF(CBUFIN(jchar:jchar).EQ.CDELIM.AND.IQUOTE.NE.1) GOTO 100        ! WAS IT A BLANK?
      IF(CBUFIN(jchar:jchar).EQ. tab .AND. IQUOTE.NE.1) GOTO 100        ! WAS IT A BLANK?
      IF(CBUFIN(jchar:jchar).EQ. null .AND. IQUOTE.NE.1) GOTO 100       ! WAS IT A NULL?
      IF(CBUFIN(jchar:jchar).EQ. comma .AND. IQUOTE.NE.1) GOTO 100       ! WAS IT A COMMA?
      jchar=jchar+1
      IF( CBUFIN(jchar:jchar) .EQ. QUOTE .AND. IQUOTE .EQ. 1 ) THEN
          IQUOTE = 2
          GOTO 100
      ENDIF
      NCHARS=NCHARS+1                                                   ! THE CURREN CHARACTER IS NOT A BLANK OR A QUOTE
      GO TO 40                                                          ! GO LOOK AT THE NEXT CHARACTER
  100 CONTINUE
      IF( NCHARS .EQ. 0 ) RETURN                                        ! DON'T TRY TO MOVE ZERO CHARACTERS!!
      IF(IQUOTE.ne.1) GO TO 110                                         ! THE CURRENT CHARACTER IS A BLANK WITHIN QUOTES
      jchar=jchar+1
      nchars=nchars+1
      go to 40
c**** end the returned string with a null character so that c rcognizes
c**** the end of string.
  110 IF( nchars .GT. 0 ) CBUFO(1:NCHARS)=CBUFIN(ISTART:jchar)
      cbufo(nchars+1:nchars+1) = null
      jchar=jchar+1                                                     ! STRIP OFF THE BLANK
      RETURN
      END
c
c
      SUBROUTINE GETOKE1C(CBUFO,NCHARS)
c     This is the same as GETOKE, except COMMA is also a delimiter.
c   ( I don't know what will break if I just add it to geoke, so make
c   it a unique entry point).
c  mod 25 Aug 98 - Make a NULL be a delimiter
c  mod 2 May 98 - Make a COMMA be a delimiter
c               - Don't honor the comment convention.
c mod 22 July 2007 - Look out for ,, - vacuous field
c
      CHARACTER*1 CDELIM,QUOTE,tab, null, comma
      CHARACTER*200 CBUFIN
      COMMON /SIOLN3/ CBUFIN
      COMMON /sioln4/ jchar,NCBUF
      CHARACTER* (*) CBUFO
      DATA CDELIM/' '/, QUOTE/''''/, comma/','/

      cbufo = ' '
      NCHARS=0                                                          ! COUNT THE NON BLANK CHARACTERS IN THE TOKEN
c****   look out for vacuous field ( ,, )
      IF( cbufin(jchar:jchar) .EQ. comma ) GOTO 110
      IQUOTE=0                                                          ! COUNT THE QUOTES IN THE TOKEN
      tab = CHAR(9)                                                     ! the tab character
      null = CHAR(0)                                                    ! NULL
      IF(jchar.LT.1) jchar=1
   10 CONTINUE
      IF( CBUFIN(jchar:jchar) .NE. CDELIM .AND.
     *   cbufin(jchar:jchar) .NE. null .AND.
     *   cbufin(jchar:jchar) .NE. comma .AND.
     *   cbufin(jchar:jchar) .NE. tab ) GOTO 20                         ! STRIP OFF LEADING BLANKS
      jchar=jchar+1
      IF(jchar.GT.NCBUF) RETURN
      GO TO 10
   20 IF(CBUFIN(jchar:jchar).NE.QUOTE) GO TO 30                         ! IS IT A QUOTE?
      ISTART=jchar+1                                                    ! STRIP OF THE LEADING QUOTE
      IQUOTE=1                                                          ! SIGNAL THAT THE STRING STARTED WITH A QUOTE
      GO TO 40
   30 CONTINUE
c**** toss out non ASCII characters
      IF( ICHAR(cbufin(jchar:jchar)) .LT. 32 ) THEN
          jchar = jchar + 1
          GOTO 10
      ENDIF
      ISTART=jchar                                                      ! THE FIRST CHARACTER OF THE TOKEN TO BE RETURNED
   40 CONTINUE                                                          ! NOW FIND THE END OF THE TOKEN
      IF(jchar.GT.NCBUF) GO TO 110                                      ! ARE WE AT THE END OF THE BUFFER?
      IF(CBUFIN(jchar:jchar).EQ.CDELIM.AND.IQUOTE.NE.1) GOTO 100        ! WAS IT A BLANK?
      IF(CBUFIN(jchar:jchar).EQ. tab .AND. IQUOTE.NE.1) GOTO 100        ! WAS IT A BLANK?
      IF(CBUFIN(jchar:jchar).EQ. null .AND. IQUOTE.NE.1) GOTO 100       ! WAS IT A NULL?
      IF(CBUFIN(jchar:jchar).EQ. comma .AND. IQUOTE.NE.1) GOTO 100       ! WAS IT A COMMA?
      jchar=jchar+1
      IF( CBUFIN(jchar:jchar) .EQ. QUOTE .AND. IQUOTE .EQ. 1 ) THEN
          IQUOTE = 2
          GOTO 100
      ENDIF
      NCHARS=NCHARS+1                                                   ! THE CURREN CHARACTER IS NOT A BLANK OR A QUOTE
      GO TO 40                                                          ! GO LOOK AT THE NEXT CHARACTER
  100 CONTINUE
      IF( NCHARS .EQ. 0 ) RETURN                                        ! DON'T TRY TO MOVE ZERO CHARACTERS!!
      IF(IQUOTE.ne.1) GO TO 110                                         ! THE CURRENT CHARACTER IS A BLANK WITHIN QUOTES
      jchar=jchar+1
      nchars=nchars+1
      go to 40
c**** end the returned string with a null character so that c rcognizes
c**** the end of string.
  110 IF( nchars .GT. 0 ) CBUFO(1:NCHARS)=CBUFIN(ISTART:jchar)
      cbufo(nchars+1:nchars+1) = null
      jchar=jchar+1                                                     ! STRIP OFF THE BLANK
      RETURN
      END

max_stars_repo_name: apthorpe/sioseis

max_issues_repo_name: henkart/sioseis

max_forks_repo_name: henkart/sioseis


text:
!< Jiang-Shu and Gerolymos-Senechal-Vallet optimal weights.
module wenoof_optimal_weights_js
!< Jiang-Shu and Gerolymos-Senechal-Vallet optimal weights.
!<
!< @note The provided WENO optimal weights implements the optimal weights defined in *Efficient Implementation of Weighted ENO
!< Schemes*, Guang-Shan Jiang, Chi-Wang Shu, JCP, 1996, vol. 126, pp. 202--228, doi:10.1006/jcph.1996.0130 and
!< *Very-high-order weno schemes*, G. A. Gerolymos, D. Senechal, I. Vallet, JCP, 2009, vol. 228, pp. 8481-8524,
!< doi:10.1016/j.jcp.2009.07.039

use penf, only : I_P, R_P
use wenoof_optimal_weights

implicit none
private
public :: optimal_weights_js
public :: optimal_weights_js_constructor
public :: create_optimal_weights_js_constructor

type, extends(optimal_weights_constructor) :: optimal_weights_js_constructor
  !< Jiang-Shu and Gerolymos-Senechal-Vallet optimal weights object constructor.
endtype optimal_weights_js_constructor

type, extends(optimal_weights):: optimal_weights_js
  !< Jiang-Shu and Gerolymos-Senechal-Vallet optimal weights object.
  !<
  !< @note The provided WENO optimal weights implements the optimal weights defined in *Efficient Implementation of Weighted ENO
  !< Schemes*, Guang-Shan Jiang, Chi-Wang Shu, JCP, 1996, vol. 126, pp. 202--228, doi:10.1006/jcph.1996.0130 and
  !< *Very-high-order weno schemes*, G. A. Gerolymos, D. Senechal, I. Vallet, JCP, 2009, vol. 228, pp. 8481-8524,
  !< doi:10.1016/j.jcp.2009.07.039
  contains
    ! deferred public methods
    procedure, pass(self) :: compute     !< Compute weights.
    procedure, nopass     :: description !< Return weights string-description.
endtype optimal_weights_js

contains
  ! public non TBP
  subroutine create_optimal_weights_js_constructor(S, constructor)
  !< Create optimal weights constructor.
  integer(I_P),                                    intent(in)  :: S           !< Stencils dimension.
  class(optimal_weights_constructor), allocatable, intent(out) :: constructor !< Optimal weights constructor.

  allocate(optimal_weights_js_constructor :: constructor)
  constructor%S = S
  endsubroutine create_optimal_weights_js_constructor

  ! deferred public methods
  pure subroutine compute(self, S)
  !< Compute weights.
  class(optimal_weights_js), intent(inout) :: self !< Optimal weights.
  integer(I_P),              intent(in)    :: S    !< Number of stencils used.

  associate(opt => self%opt)
    select case(S)
      case(2) ! 3rd order
        ! 1 => left interface (i-1/2)
        opt(1, 0) = 2._R_P/3._R_P ! stencil 0
        opt(1, 1) = 1._R_P/3._R_P ! stencil 1
        ! 2 => right interface (i+1/2)
        opt(2, 0) = 1._R_P/3._R_P ! stencil 0
        opt(2, 1) = 2._R_P/3._R_P ! stencil 1
      case(3) ! 5th order
        ! 1 => left interface (i-1/2)
        opt(1, 0) = 0.3_R_P ! stencil 0
        opt(1, 1) = 0.6_R_P ! stencil 1
        opt(1, 2) = 0.1_R_P ! stencil 2
        ! 2 => right interface (i+1/2)
        opt(2, 0) = 0.1_R_P ! stencil 0
        opt(2, 1) = 0.6_R_P ! stencil 1
        opt(2, 2) = 0.3_R_P ! stencil 2
      case(4) ! 7th order
        ! 1 => left interface (i-1/2)
        opt(1, 0) =  4._R_P/35._R_P ! stencil 0
        opt(1, 1) = 18._R_P/35._R_P ! stencil 1
        opt(1, 2) = 12._R_P/35._R_P ! stencil 2
        opt(1, 3) =  1._R_P/35._R_P ! stencil 3
        ! 2 => right interface (i+1/2)
        opt(2, 0) =  1._R_P/35._R_P ! stencil 0
        opt(2, 1) = 12._R_P/35._R_P ! stencil 1
        opt(2, 2) = 18._R_P/35._R_P ! stencil 2
        opt(2, 3) =  4._R_P/35._R_P ! stencil 3
      case(5) ! 9th order
        ! 1 => left interface (i-1/2)
        opt(1, 0) =  5._R_P/126._R_P ! stencil 0
        opt(1, 1) = 20._R_P/63._R_P  ! stencil 1
        opt(1, 2) = 10._R_P/21._R_P  ! stencil 2
        opt(1, 3) = 10._R_P/63._R_P  ! stencil 3
        opt(1, 4) =  1._R_P/126._R_P ! stencil 4
        ! 2 => right interface (i+1/2)
        opt(2, 0) =  1._R_P/126._R_P ! stencil 0
        opt(2, 1) = 10._R_P/63._R_P  ! stencil 1
        opt(2, 2) = 10._R_P/21._R_P  ! stencil 2
        opt(2, 3) = 20._R_P/63._R_P  ! stencil 3
        opt(2, 4) =  5._R_P/126._R_P ! stencil 4
      case(6) ! 11th order
        ! 1 => left interface (i-1/2)
        opt(1, 0) =   1._R_P/77._R_P  ! stencil 0
        opt(1, 1) =  25._R_P/154._R_P ! stencil 1
        opt(1, 2) = 100._R_P/231._R_P ! stencil 2
        opt(1, 3) =  25._R_P/77._R_P  ! stencil 3
        opt(1, 4) =   5._R_P/77._R_P  ! stencil 4
        opt(1, 5) =   1._R_P/462._R_P ! stencil 5
        ! 2 => right interface (i+1/2)
        opt(2, 0) =   1._R_P/462._R_P ! stencil 0
        opt(2, 1) =   5._R_P/77._R_P  ! stencil 1
        opt(2, 2) =  25._R_P/77._R_P  ! stencil 2
        opt(2, 3) = 100._R_P/231._R_P ! stencil 3
        opt(2, 4) =  25._R_P/154._R_P ! stencil 4
        opt(2, 5) =   1._R_P/77._R_P  ! stencil 5
      case(7) ! 13th order
        ! 1 => left interface (i-1/2)
        opt(1, 0) =   7._R_P/1716._R_P ! stencil 0
        opt(1, 1) =  21._R_P/286._R_P  ! stencil 1
        opt(1, 2) = 175._R_P/572._R_P  ! stencil 2
        opt(1, 3) = 175._R_P/429._R_P  ! stencil 3
        opt(1, 4) = 105._R_P/572._R_P  ! stencil 4
        opt(1, 5) =   7._R_P/286._R_P  ! stencil 5
        opt(1, 6) =   1._R_P/1716._R_P ! stencil 6
        ! 2 => right interface (i+1/2)
        opt(2, 0) =   1._R_P/1716._R_P ! stencil 0
        opt(2, 1) =   7._R_P/286._R_P  ! stencil 1
        opt(2, 2) = 105._R_P/572._R_P  ! stencil 2
        opt(2, 3) = 175._R_P/429._R_P  ! stencil 3
        opt(2, 4) = 175._R_P/572._R_P  ! stencil 4
        opt(2, 5) =  21._R_P/286._R_P  ! stencil 5
        opt(2, 6) =   7._R_P/1716._R_P ! stencil 6
      case(8) ! 15th order
        ! 1 => left interface (i-1/2)
        opt(1, 0) =   8._R_P/6435._R_P ! stencil 0
        opt(1, 1) = 196._R_P/6435._R_P ! stencil 1
        opt(1, 2) = 392._R_P/2145._R_P ! stencil 2
        opt(1, 3) = 490._R_P/1287._R_P ! stencil 3
        opt(1, 4) = 392._R_P/1287._R_P ! stencil 4
        opt(1, 5) = 196._R_P/2145._R_P ! stencil 5
        opt(1, 6) =  56._R_P/6435._R_P ! stencil 6
        opt(1, 7) =   1._R_P/6435._R_P ! stencil 7
        ! 2 => right interface (i+1/2)
        opt(2, 0) =   1._R_P/6435._R_P ! stencil 0
        opt(2, 1) =  56._R_P/6435._R_P ! stencil 1
        opt(2, 2) = 196._R_P/2145._R_P ! stencil 2
        opt(2, 3) = 392._R_P/1287._R_P ! stencil 3
        opt(2, 4) = 490._R_P/1287._R_P ! stencil 4
        opt(2, 5) = 392._R_P/2145._R_P ! stencil 5
        opt(2, 6) = 196._R_P/6435._R_P ! stencil 6
        opt(2, 7) =   8._R_P/6435._R_P ! stencil 7
      case(9) ! 17th order
        ! 1 => left interface (i-1/2)
        opt(1, 0) =    9._R_P/24310._R_P ! stencil 0
        opt(1, 1) =  144._R_P/12155._R_P ! stencil 1
        opt(1, 2) = 1176._R_P/12155._R_P ! stencil 2
        opt(1, 3) = 3528._R_P/12155._R_P ! stencil 3
        opt(1, 4) =  882._R_P/2431._R_P  ! stencil 4
        opt(1, 5) = 2352._R_P/12155._R_P ! stencil 5
        opt(1, 6) =  504._R_P/12155._R_P ! stencil 6
        opt(1, 7) =   36._R_P/12155._R_P ! stencil 7
        opt(1, 8) =    1._R_P/24310._R_P ! stencil 8
        ! 2 => right interface (i+1/2)
        opt(2, 0) =    1._R_P/24310._R_P ! stencil 0
        opt(2, 1) =   36._R_P/12155._R_P ! stencil 1
        opt(2, 2) =  504._R_P/12155._R_P ! stencil 2
        opt(2, 3) = 2352._R_P/12155._R_P ! stencil 3
        opt(2, 4) =  882._R_P/2431._R_P  ! stencil 4
        opt(2, 5) = 3528._R_P/12155._R_P ! stencil 5
        opt(2, 6) = 1176._R_P/12155._R_P ! stencil 6
        opt(2, 7) =  144._R_P/12155._R_P ! stencil 7
        opt(2, 8) =    9._R_P/24310._R_P ! stencil 8
    endselect
  endassociate
  endsubroutine compute

  pure function description() result(string)
  !< Return string-description of weights.
  character(len=:), allocatable :: string           !< String-description.
  character(len=1), parameter   :: nl=new_line('a') !< New line character.

  string = 'WENO optimal weights,'//nl
  string = string//'  based on the work by Jiang and Shu "Efficient Implementation of Weighted ENO Schemes", see '// &
           'JCP, 1996, vol. 126, pp. 202--228, doi:10.1006/jcph.1996.0130 and'//nl
  string = string//'  on the work by Gerolymos, Senechal and Vallet "Very-high-order weno schemes", see '// &
           'JCP, 2009, vol. 228, pp. 8481--8524, doi:10.1016/j.jcp.2009.07.039'//nl
  string = string//'    The optimal weights are allocated in a two-dimensional array, in which the first index'//nl
  string = string//'    is the face selected (1 => i-1/2, 2 => i+1/2) and the second index is the number of the stencil '//nl
  string = string//'    (from 0 to S-1)'
  endfunction description
endmodule wenoof_optimal_weights_js

max_stars_repo_name: giacrossi/WenOOF

max_issues_repo_name: giacrossi/WenOOF

max_forks_repo_name: giacrossi/WenOOF


!123456
      real :: E_max, E_min
      real :: fermi
      integer :: nedos,spin, nion, soc
      real, allocatable :: s(:,:), p(:,:), d(:,:)
      real, allocatable :: s_up(:,:), s_dwn(:,:), p_up(:,:), p_dwn(:,:)
      real, allocatable :: d_up(:,:), d_dwn(:,:), s_total(:,:)
      real, allocatable :: p_total(:,:), d_total(:,:), DOS(:)
      real, allocatable :: integrated_DOS(:), DOS_up(:)
      real, allocatable :: DOS_dwn(:), integrated_DOS_up(:)
      real, allocatable :: integrated_DOS_dwn(:)
      real, allocatable :: s_x(:,:), s_y(:,:), s_z(:,:)
      real, allocatable :: p_x(:,:), p_y(:,:), p_z(:,:)
      real, allocatable :: d_x(:,:), d_y(:,:), d_z(:,:), energy(:)
      integer :: atom1, atom2

!123456
      OPEN(1,file='band.dat')
      read(1,*) fermi
      read(1,*) soc
      read(1,*) atom1, atom2

      open(15,file='EIGENVAL')
      read(15,*) c1, c2, c3, spin

      open(2,file='DOSCAR')
      read(2,*) nion, c1, c2, c3
      read(2,*)
      read(2,*)
      read(2,*)
      read(2,*)
      read(2,*) E_max, E_min, nedos, E_fermi, c1

      write(6,*) E_max
       If(spin.eq.1 .and. soc.eq.1) then

      allocate(energy(nedos), DOS(nedos),integrated_DOS(nedos))

       do i = 1,nedos

        read(2,*) energy(i), DOS(i), integrated_DOS(i)

         enddo

      open(4,file='result_total.out')
200   format(f14.10,3x,f14.10,3x,f14.10)

      do i = 1,nedos
      write(4,200) energy(i), DOS(i), integrated_DOS(i)
      enddo

       allocate(s(nion,nedos),p(nion,nedos),d(nion,nedos))

       do i = 1, atom1
       read(2,*)

       do j = 1, nedos
       read(2,*) energy(j), s(i,j), p(i,j), d(i,j)
       enddo

       enddo

       do i = atom1+1, atom2
       read(2,*)

       do j = 1, nedos
       read(2,*) energy(j), s(i,j), p(i,j), d(i,j)
       enddo

       enddo

      open(3,file='result_atom1.out')
100   format(f14.10,3x,f14.10,3x,f14.10,3x,f14.10)

      do i = 1, atom1
      do j = 1,nedos
      write(3,100) energy(j), s(i,j),p(i,j),d(i,j)
      enddo
      write(3,*)
      enddo

      open(5,file='result_atom2.out')

      do i = atom1+1, atom2
      do j = 1, nedos
      write(5,100) energy(j), s(i,j), p(i,j), d(i,j)
      enddo
      write(5,*)
      enddo

      endif

      if(spin.eq.1 .and. soc.eq.2) then

      allocate(energy(nedos), DOS(nedos),integrated_DOS(nedos))

       do i = 1,nedos

        read(2,*) energy(i), DOS(i), integrated_DOS(i)

         enddo

      open(4,file='result_total.out')

      do i = 1,nedos
      write(4,200) energy(i), DOS(i), integrated_DOS(i)
      enddo

       allocate(s_total(nion,nedos),s_x(nion,nedos),s_y(nion,nedos),s_z(nion,nedos))
       allocate(p_total(nion,nedos),p_x(nion,nedos),p_y(nion,nedos),p_z(nion,nedos))
       allocate(d_total(nion,nedos),d_x(nion,nedos),d_y(nion,nedos),d_z(nion,nedos))

       do i = 1, atom1
       read(2,*)

       do j = 1, nedos
       read(2,*) energy(j), s_total(i,j), s_x(i,j),s_y(i,j),s_z(i,j),p_total(i,j),p_x(i,j),p_y(i,j),p_z(i,j),d_total(i,j),d_x(i,j),d_y(i,j),d_z(i,j)
       enddo

       enddo

       do i = atom1+1, atom2
       read(2,*)

       do j = 1, nedos
       read(2,*) energy(j), s_total(i,j),s_x(i,j),s_y(i,j),s_z(i,j),p_total(i,j),p_x(i,j),p_y(i,j),p_z(i,j),d_total(i,j),d_x(i,j),d_y(i,j),d_z(i,j)
       enddo

       enddo

      open(3,file='result_atom1.out')
101   format(f14.10,3x,f14.10,3x,f14.10,3x,f14.10,3x,f14.10,3x,f14.10,3x,f14.10,3x,f14.10,3x,f14.10,3x,f14.10,3x,f14.10,3x,f14.10,3x,f14.10)

      do i = 1, atom1
      do j = 1,nedos
      write(3,101) energy(j), s_total(i,j),s_x(i,j),s_y(i,j),s_z(i,j),p_total(i,j),p_x(i,j),p_y(i,j),p_z(i,j),d_total(i,j),d_x(i,j),d_y(i,j),d_z(i,j)
      enddo
      write(3,*)
      enddo

      open(5,file='result_atom2.out')

      do i = atom1+1, atom2
      do j = 1,nedos
      write(5,101) energy(j),s_total(i,j),s_x(i,j),s_y(i,j),s_z(i,j),p_total(i,j),p_x(i,j),p_y(i,j),p_z(i,j),d_total(i,j),d_x(i,j),d_y(i,j),d_z(i,j)
      enddo
      write(5,*)
      enddo

      endif

       If(spin.eq.2) then

      allocate(energy(nedos), DOS_up(nedos),integrated_DOS_up(nedos))
      allocate(DOS_dwn(nedos),integrated_DOS_dwn(nedos))

       do i = 1,nedos

        read(2,*) energy(i), DOS_up(i),DOS_dwn(i),integrated_DOS_up(i),integrated_DOS_dwn(i)

         enddo

      open(4,file='result_total.out')
201   format(f14.10,3x,f14.10,3x,f14.10,3x,f14.10,3x,f14.10)
      do i = 1,nedos

     write(4,201) energy(i),DOS_up(i),DOS_dwn(i),integrated_DOS_up(i),integrated_DOS_dwn(i)
      enddo

       allocate(s_up(nion,nedos),p_up(nion,nedos),d_up(nion,nedos))
       allocate(s_dwn(nion,nedos),p_dwn(nion,nedos),d_dwn(nion,nedos))

       do i = 1, atom1
       read(2,*)

       do j = 1, nedos
       read(2,*) energy(j), s_up(i,j), s_dwn(i,j),p_up(i,j),p_dwn(i,j), d_up(i,j),d_dwn(i,j)
       enddo

       enddo

       do i = atom1+1,atom2
       read(2,*)

       do j = 1, nedos
       read(2,*) energy(j), s_up(i,j), s_dwn(i,j),p_up(i,j),p_dwn(i,j),d_up(i,j),d_dwn(i,j)
       enddo

       enddo

      open(3,file='result_atom1.out')
102   format(f14.10,3x,f14.10,3x,f14.10,3x,f14.10,3x,f14.10,3x,f14.10,3x,f14.10)

      do i = 1, atom1
      do j = 1,nedos
      write(3,102) energy(j), s_up(i,j),s_dwn(i,j),p_up(i,j),p_dwn(i,j),d_up(i,j),d_dwn(i,j)
      enddo
      write(3,*)
      enddo

      open(5,file='result_atom2.out')

      do i = atom1+1, atom2
      do j = 1,nedos
      write(5,102) energy(j),s_up(i,j),s_dwn(i,j),p_up(i,j),p_dwn(i,j),d_up(i,j),d_dwn(i,j)
      enddo
      write(5,*)
      enddo

      endif

      stop

      end


module read_wrf_mass_guess_mod
use abstract_read_wrf_mass_guess_mod
  type, extends(abstract_read_wrf_mass_guess_class) :: read_wrf_mass_guess_class
  contains
    procedure, pass(this) :: read_wrf_mass_binary_guess => read_wrf_mass_binary_guess_dummy
    procedure, pass(this) :: read_wrf_mass_netcdf_guess => read_wrf_mass_netcdf_guess_dummy
  end type read_wrf_mass_guess_class
contains
  subroutine read_wrf_mass_binary_guess_dummy(this,mype)
  !$$$  subprogram documentation block
  !                .      .    .                                       .
  ! subprogram:    read_wrf_mass_binary_guess
  !   prgmmr:
  !
  ! abstract:
  !
  ! program history log:
  !   2009-12-07  lueken - added subprogram doc block and implicit none
  !
  !   input argument list:
  !
  !   output argument list:
  !
  ! attributes:
  !   language: f90
  !   machine:  ibm RS/6000 SP
  !
  !$$$ end documentation block
    use kinds,only: i_kind
    implicit none
    class(read_wrf_mass_guess_class),intent(inout) :: this
    integer(i_kind),intent(in)::mype
    write(6,*)'READ_WRF_MASS_BINARY_GUESS:  dummy routine, does nothing!'
  end subroutine read_wrf_mass_binary_guess_dummy

  subroutine read_wrf_mass_netcdf_guess_dummy(this,mype)
  !$$$  subprogram documentation block
  !                .      .    .                                       .
  ! subprogram:    read_wrf_mass_netcdf_guess
  !   prgmmr:
  !
  ! abstract:
  !
  ! program history log:
  !   2009-12-07  lueken - added subprogram doc block and implicit none
  !
  !   input argument list:
  !
  !   output argument list:
  !
  ! attributes:
  !   language: f90
  !   machine:  ibm RS/6000 SP
  !
  !$$$ end documentation block
    use kinds,only: i_kind
    implicit none
    class(read_wrf_mass_guess_class),intent(inout) :: this
    integer(i_kind),intent(in)::mype
    write(6,*)'READ_WRF_MASS_NETCDF_GUESS:  dummy routine, does nothing!'
  end subroutine read_wrf_mass_netcdf_guess_dummy
end module read_wrf_mass_guess_mod


text:
      subroutine res (n, t, y, v, r, ires)
c This subroutine computes the residual vector
c   r = g(t,y) - A(t,y)*v .
c It uses nm1 = n - 1 from Common.
c If ires = -1, only g(t,y) is returned in r, since A(t,y) does
c not depend on y.
c
      integer i, ires, n, nm1
      double precision t, y, v, r, r4d, eodsq, one, four, six,
     1   fact1, fact4
      dimension y(*), v(*), r(*)
      dimension n(*)
      integer itemp
      common /test1/ r4d, eodsq, nm1
      data one /1.0d0/, four /4.0d0/, six /6.0d0/
c
      itemp = n(1)
      call gfun (itemp, t, y, r)
      if (ires .eq. -1) return
c
      fact1 = one/six
      fact4 = four/six
      r(1) = r(1) - (fact4*v(1) + fact1*v(2))
      do 10 i = 2, nm1
  10   r(i) = r(i) - (fact1*v(i-1) + fact4*v(i) + fact1*v(i+1))
      r(n(1)) = r(n(1)) - (fact1*v(nm1) + fact4*v(n(1)))
      return
c end of subroutine res for the DLSODI demonstration problem.
      end

max_stars_repo_name: sbwhitecap/clocc-hg

max_issues_repo_name: sbwhitecap/clocc-hg

max_forks_repo_name: sbwhitecap/clocc-hg


SUBROUTINE UpdateGhostLayerYCurvilinear(var,S,NNx,NNy,CompGrid,alpha,beta,GhostGridX,GhostGridY)
! By Allan P. Engsig-Karup.
USE Precision
USE DataTypes
IMPLICIT NONE
TYPE (Level_def) :: CompGrid
INTEGER :: NNx, NNy, rank, alpha, beta, i,j, GhostGridX, GhostGridY, diffa
REAL(KIND=long), DIMENSION(NNx,NNy) :: var, S
REAL(KIND=long), DIMENSION(2*alpha+1) :: Stencile, Stenciln
REAL(KIND=long) :: ex, ey, nx, ny, tmp, diagn
INTEGER, DIMENSION(2*alpha+1) :: idx
IF (alpha/=beta) THEN
	print*,'ERROR: alpha/=beta. (UpdateGhostLayerYCurvilinear)'
END IF
rank = 2*beta+1 !GD change: 2*alpha+1
DO i = -alpha,alpha
	idx(i+alpha+1) = i
END DO
! two boundaries
j=1
DO i = 1+GhostGridX, NNx-GhostGridX
	if (i-1<alpha+1) then
		diffa = (alpha+1)-(i-1)
	else if (i+1>NNx-alpha) then
		diffa = (NNx-alpha)-(i+1) !GD: change (NNx-beta)-(i+1)
	else
		diffa = 0
	end if

	! use geometric factors for wall boundary
!	ex = CompGrid%CurvilinearStuff%ex(i,j+GhostGridY)
	ey = CompGrid%CurvilinearStuff%ey(i,j+GhostGridY)
!	nx = CompGrid%CurvilinearStuff%nx(i,j+GhostGridY)
	ny = CompGrid%CurvilinearStuff%ny(i,j+GhostGridY)

	Stencile = CompGrid%CurvilinearStuff%DiffStencils%StencilG(1:rank,1+alpha-diffa,1)
	Stenciln = CompGrid%CurvilinearStuff%DiffStencils%StencilG(1:rank,1+GhostGridY,1)
	diagn = ny*Stenciln(1)

	tmp = ey * DOT_PRODUCT(Stencile,var(i+idx+diffa,j+GhostGridY) ) + ny * DOT_PRODUCT( Stenciln(2:rank),var(i,2:rank) )

    var(i,j) = (S(i,j+GhostGridY) - tmp) / diagn

END DO
j = NNy
DO i = 1+GhostGridX, NNx-GhostGridX
	if (i-1<alpha+1) then
		diffa = (alpha+1)-(i-1)
	else if (i+1>NNx-alpha) then
		diffa = (NNx-alpha)-(i+1)
	else
		diffa = 0
	end if

	! use geometric factors for wall boundary
!	ex = CompGrid%CurvilinearStuff%ex(i,j-GhostGridY)
	ey = CompGrid%CurvilinearStuff%ey(i,j-GhostGridY)
!	nx = CompGrid%CurvilinearStuff%nx(i,j-GhostGridY)
	ny = CompGrid%CurvilinearStuff%ny(i,j-GhostGridY)

	Stencile = CompGrid%CurvilinearStuff%DiffStencils%StencilG(1:rank,1+alpha-diffa,1)
	Stenciln = CompGrid%CurvilinearStuff%DiffStencils%StencilG(1:rank,rank-GhostGridY,1)

	diagn = ny*Stenciln(rank)

	tmp = ey * DOT_PRODUCT(Stencile,var(i+idx+diffa,j-GhostGridY) ) + ny * DOT_PRODUCT( Stenciln(1:rank-1),var(i,j-rank+1:j-1) )

    var(i,j) = (S(i,j-GhostGridY)- tmp) / diagn
END DO
END SUBROUTINE UpdateGhostLayerYCurvilinear


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

      SUBROUTINE OVRLAP (A, COORD, IDESS, NEESS, NNESS, IPEESS, IPNESS,
     *   LTEESS, LTNESS, FACESS, DISP, NUMNP, NDIM, NUMESS,
     *   TIME, ITMSEL, TITLE, IMAS, ISLV, NUMEL)

      DIMENSION A(*), COORD(NUMNP,*), IDESS(*), NEESS(*),
     *   NNESS(*), IPEESS(*), IPNESS(*), LTEESS(*), LTNESS(*),
     *   FACESS(*), TIME(*), DISP(NUMNP,*)
      LOGICAL ITMSEL(*)
      CHARACTER*80 TITLE, STRA
      include 'nu_mass.blk'
      include 'nu_logs.blk'
      include 'nu_ptim.blk'
      include 'nu_io.blk'
      LOGICAL ERROR

      DIMENSION CPTIME(10)

      DO 10 I=1,10
         CPTIME(I) = 0.0
   10 CONTINUE

      IFLGM = LOCINT (IMAS, NUMESS, IDESS)
      IFLGS = LOCINT (ISLV, NUMESS, IDESS)

      ERROR = .FALSE.
      IF (IFLGM .EQ. 0) THEN
         WRITE (STRA, 30) 'Master', IMAS
         CALL SQZSTR (STRA, LSTR)
         CALL PRTERR ('ERROR', STRA(:LSTR))
   30    FORMAT (1X,A,' Surface Flag ',I5,' not found. ')
         ERROR = .TRUE.
      END IF
      IF (IFLGS .EQ. 0) THEN
         WRITE (STRA, 30) 'Slave', ISLV
         CALL SQZSTR (STRA, LSTR)
         CALL PRTERR ('ERROR', STRA(:LSTR))
         ERROR = .TRUE.
      END IF
      IF (ERROR) RETURN

      WRITE (STRA, 50) IMAS, ISLV
      CALL SQZSTR (STRA, LSTR)
      DO 40 IO=IOMIN, IOMAX
         WRITE (IO, 55) STRA(:LSTR)
   40 CONTINUE
   50 FORMAT ('Checking Master Surface ',I5,' Versus Slave Surface ',I5)
   55 FORMAT (/,1X,A,/)

      NSEGM = NEESS(IFLGM)
      IPTRM = IPNESS(IFLGM)
      IEPTM = IPEESS(IFLGM)

      NSEGS = NEESS(IFLGS)
      IPTRS = IPNESS(IFLGS)

      MULT = 2 * NDIM - 2

C ... PROCESS SLAVE SET TO REMOVE DUPLICATE NODES

      CALL MDRSRV ('MAPSLV', IMPSL, MULT*NSEGS)
      CALL MDRSRV ('ITEMP',  ITMP,  MAX(NUMNP,3*NSEGM))
      CALL UNIQUE (LTNESS(IPTRS), MULT*NSEGS, A(IMPSL), A(ITMP),
     *   NIQS, NUMNP)
      CALL MDRSRV ('NIQSLV', INQS, NIQS)
      CALL TRANIQ (LTNESS(IPTRS), A(IMPSL), A(INQS), MULT*NSEGS, 1)

      CALL MDRSRV ('MINMAX', IMNMX, 2*NDIM*NSEGM)
      CALL MDRSRV ('LFACE',  ILFAC, 2*NDIM*NUMEL)
      CALL MDSTAT (NERRS, NUSED)
      IF (NERRS .GT. 0) THEN
         CALL MEMERR
         STOP
      END IF

C ... Beginning of Time Step Loop

      IF (EXODUS .AND. ISDIS) THEN
         CALL GETDSP (A(IR), DISP, NDIM, NUMNP, TIME, ITMSEL,
     *        'R', ISTAT)
         IF (ISTAT .NE. 0) GO TO 150
      END IF

 60   CONTINUE
      IF (EXODUS) THEN
         CALL GETDSP (A(IR), DISP, NDIM, NUMNP, TIME, ITMSEL,
     *        'A', ISTAT)
         IF (ISTAT .NE. 0) GO TO 150
         DO 160 IO=IOMIN, IOMAX
           WRITE (IO, 56) TREAD
 160      CONTINUE
 56      FORMAT (' At Time = ', 1PE15.8)
         IF (NDIM .EQ. 3) THEN
            CALL OVRMX3 (LTEESS(IEPTM), DISP, A(IX), NSEGM, A(IMNMX),
     *           A(INQS), NIQS, A(ITMP), LTNESS(IPTRM),
     *           NUMIN, NUMFAC, NUMON, NUMEL, A(ILFAC), NUMNP)
         ELSE
            CALL OVRMX2 (LTEESS(IEPTM), DISP, A(IX), NSEGM, A(IMNMX),
     *           A(INQS), NIQS, A(ITMP), LTNESS(IPTRM),
     *           NUMIN, NUMFAC, NUMON, NUMEL, A(ILFAC), NUMNP)
         END IF

      ELSE
C ... Not EXODUS
         IF (NDIM .EQ. 3) THEN
            CALL OVRMX3 (LTEESS(IEPTM), COORD, A(IX), NSEGM, A(IMNMX),
     *           A(INQS), NIQS, A(ITMP), LTNESS(IPTRM),
     *           NUMIN, NUMFAC, NUMON, NUMEL, A(ILFAC), NUMNP)
         ELSE
            CALL OVRMX2 (LTEESS(IEPTM), COORD, A(IX), NSEGM, A(IMNMX),
     *           A(INQS), NIQS, A(ITMP), LTNESS(IPTRM),
     *           NUMIN, NUMFAC, NUMON, NUMEL, A(ILFAC), NUMNP)
         END IF
      END IF

      IF (NUMIN .GT. 0) THEN
         DO 70 IO=IOMIN, IOMAX
            WRITE (IO, 80) NUMIN
   70    CONTINUE
   80    FORMAT (' Warning -',I5,' Nodes penetrate master surface.')
      ELSE
         DO 90 IO=IOMIN, IOMAX
            WRITE (IO, 100)
   90    CONTINUE
  100    FORMAT (' No nodes found penetrating master surface.')
      END IF
      IF (NUMFAC .GT. 0) THEN
         DO 110 IO=IOMIN, IOMAX
            WRITE (IO, 120) NUMFAC
  110    CONTINUE
  120    FORMAT (' ',I5,' Nodes on an element face.')
      END IF
      IF (NUMON .GT. 0) THEN
         DO 130 IO=IOMIN, IOMAX
            WRITE (IO, 140) NUMON
  130    CONTINUE
  140    FORMAT (' ',I5,' Nodes in both side sets.')
      END IF
      IF (EXODUS) GO TO 60
 150  CONTINUE

      CALL MDDEL ('LFACE')
      CALL MDDEL ('MINMAX')
      CALL MDDEL ('NIQSLV')
      CALL MDDEL ('ITEMP' )
      CALL MDDEL ('MAPSLV')
      RETURN
      END

max_stars_repo_name: jschueller/seacas

max_issues_repo_name: jschueller/seacas

max_forks_repo_name: jschueller/seacas


text:
module hip_tap_plugin
    use :: hip_tap_driver_module, only : hip_tap_driver

    implicit none
    public
end module hip_tap_plugin

max_stars_repo_name: cheshyre/ntcl-tensor

max_issues_repo_name: cheshyre/ntcl-tensor

max_forks_repo_name: cheshyre/ntcl-tensor


text:
REAL*8 FUNCTION MakeGridPoint(cilm, lmax, lat, longitude, csphase, norm, dealloc)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!	This function will determine the value at a given latitude and
!	longitude corresponding to the given set of spherical harmonics.
!	Latitude and Longitude are assumed to be in DEGREES!
!
!	Calling Parameters:
!		IN
!			cilm:		Spherical harmonic coefficients, with dimensions
!					(2, lmax+1, lmax+1).
!			lmax:		Maximum degree used in the expansion.
!			lat:		latitude (degrees).
!			long:		longitude(degrees).
!		OPTIONAL (IN)
!			norm		Spherical harmonic normalization:
!						(1) "geodesy" (default)
!						(2) Schmidt
!						(3) unnormalized
!						(4) orthonormalized
!			csphase:		1: Do not include the phase factor of (-1)^m
!						-1: Apply the phase factor of (-1)^m.
!			dealloc		If (1) Deallocate saved memory in Legendre function routines.
!					Default (0) is not to deallocate memory.
!
!
!	Dependencies:		PlmBar, PlBar, PlmSchmidt, PlmON, CSPHASE_DEFAULT
!
!	Notes:
!		1.	If lmax is greater than the the maximum spherical harmonic
!			degree of the input file, then this file will be ZERO PADDED!
!			(i.e., those degrees after lmax are assumed to be zero).
!
!	Written by Mark Wieczorek (June 2004)
!
!	August 8, 2012. Modified to precompute sin and cosine terms using a recusion relationship. Also, the default has been modified to
!			save variables used in the Legendre function routines. To deallocate this space, specify
!			the new optional parameter DEALLOC to be 1.
!
!	Copyright (c) 2005-2012, Mark A. Wieczorek
!	All rights reserved.
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	use SHTOOLS, only: PlmBar, PLegendreA, PlmSchmidt, PlmON

	implicit none

	real*8, intent(in)            :: cilm(2,lmax+1,lmax+1), lat, longitude
	integer, intent(in)           :: lmax, csphase
	integer, intent(in), optional :: norm, dealloc

	real*8                        :: pi, x, expand, lon
	integer                       :: index, l, m, l1, m1, astat(3)
	real*8, allocatable           :: pl(:), cosm(:), sinm(:)

 	allocate(pl((lmax+1)*(lmax+2)/2), stat = astat(1))
 	allocate(cosm(lmax+1), stat = astat(2))
 	allocate(sinm(lmax+1), stat = astat(3))

	if (sum(astat(1:3)) /= 0) then
		print*, "Error --- MakeGridPoint"
		print*, "Cannot allocate memory for arrays PL, MCOS and MSIN", astat(1), astat(2), astat(3)
		stop
	endif

	pi  = acos(-1.0d0)
	x   = sin(lat*pi/180.0d0)
	lon = longitude * pi/180.0d0

	if (present(norm)) then
		if (norm == 1) call PlmBar(pl, lmax, x, csphase = csphase)
		if (norm == 2) call PlmSchmidt(pl, lmax, x, csphase = csphase)
		if (norm == 3) call PLegendreA(pl, lmax, x, csphase = csphase)
		if (norm == 4) call PlmON(pl, lmax, x, csphase = csphase)
	else
		call PlmBar(pl, lmax, x, csphase = csphase)
	endif

	expand = 0.0d0

	! Precompute sines and cosines. Use multiple angle identity to minimize number of calls to SIN and COS.
	sinm(1) = 0.0d0
	cosm(1) = 1.0d0
	if (lmax > 0) then
		sinm(2) = sin(lon)
		cosm(2) = cos(lon)
	endif
	do m=2, lmax, 1
		m1 = m+1
		sinm(m1) = 2 * sinm(m)*cosm(2) - sinm(m-1)
		cosm(m1) = 2 * cosm(m)*cosm(2) - cosm(m-1)
	enddo

	do l = lmax, 0, -1
		l1 = l+1
		index = (l+1)*l/2 + 1
		expand = expand + cilm(1,l1,1) * pl(index)

		do m = 1, l, 1
			m1 = m+1
			index = index + 1
			expand = expand + ( cilm(1,l1,m1) * cosm(m1) + &
				cilm(2,l1,m1) * sinm(m1) ) * pl(index)
		enddo
	enddo

	MakeGridPoint = expand

	! deallocate memory
	if (present(dealloc)) then
		if (dealloc == 1) then
			if (present(norm)) then
				if (norm == 1) call PlmBar(pl, -1, x, csphase = csphase)
				if (norm == 2) call PlmSchmidt(pl, -1, x, csphase = csphase)
				if (norm == 4) call PlmON(pl, -1, x, csphase = csphase)
			else
				call PlmBar(pl, -1, x, csphase = csphase)
			endif
		endif
	endif

	deallocate(pl)
	deallocate(cosm)
	deallocate(sinm)

END FUNCTION MakeGridPoint

max_stars_repo_name: dchandan/PySHTOOLS

max_issues_repo_name: dchandan/PySHTOOLS

max_forks_repo_name: dchandan/PySHTOOLS


text:
FUNCTION ppchi2(p, v, g) RESULT(fn_val)

! N.B. Argument IFAULT has been removed.

! This version by Alan Miller
! amiller @ bigpond.net.au
! Latest revision - 27 October 2000

!  Algorithm AS 91   Appl. Statist. (1975) Vol.24, P.35

!  To evaluate the percentage points of the chi-squared
!  probability distribution function.

!  p must lie in the range 0.000002 to 0.999998,
!  v must be positive,
!  g must be supplied and should be equal to ln(gamma(v/2.0))

!  Incorporates the suggested changes in AS R85 (vol.40(1), pp.233-5, 1991)
!  which should eliminate the need for the limited range for p above,
!  though these limits have not been removed from the routine.

!  If IFAULT = 4 is returned, the result is probably as accurate as
!  the machine will allow.

!  Auxiliary routines required: PPND = AS 111 (or AS 241) and GAMMAD = AS 239.

IMPLICIT NONE
INTEGER, PARAMETER    :: dp = SELECTED_REAL_KIND(12, 60)

REAL (dp), INTENT(IN)  :: p
REAL (dp), INTENT(IN)  :: v
REAL (dp), INTENT(IN)  :: g
REAL (dp)              :: fn_val

INTERFACE
  FUNCTION gammad(x, p) RESULT(fn_val)
    IMPLICIT NONE
    INTEGER, PARAMETER    :: dp = SELECTED_REAL_KIND(12, 60)
    REAL (dp), INTENT(IN) :: x, p
    REAL (dp)             :: fn_val
  END FUNCTION gammad

  SUBROUTINE ppnd16 (p, normal_dev, ifault)
    IMPLICIT NONE
    INTEGER, PARAMETER      :: dp = SELECTED_REAL_KIND(12, 60)
    REAL (dp), INTENT(IN)   :: p
    INTEGER, INTENT(OUT)    :: ifault
    REAL (dp), INTENT(OUT)  :: normal_dev
  END SUBROUTINE ppnd16
END INTERFACE

! Local variables

REAL (dp)  :: a, b, c, p1, p2, q, s1, s2, s3, s4, s5, s6, t, x, xx
INTEGER    :: i, if1

INTEGER, PARAMETER    :: maxit = 20
REAL (dp), PARAMETER  :: aa = 0.6931471806_dp, e = 0.5e-06_dp,         &
                         pmin = 0.000002_dp, pmax = 0.999998_dp,       &
                         zero = 0.0_dp, half = 0.5_dp, one = 1.0_dp,   &
                         two = 2.0_dp, three = 3.0_dp, six = 6.0_dp,   &
                         c1 = 0.01_dp, c2 = 0.222222_dp, c3 = 0.32_dp, &
                         c4 = 0.4_dp, c5 = 1.24_dp, c6 = 2.2_dp,       &
                         c7 = 4.67_dp, c8 = 6.66_dp, c9 = 6.73_dp,     &
                         c10 = 13.32_dp, c11 = 60.0_dp, c12 = 70.0_dp, &
                         c13 = 84.0_dp, c14 = 105.0_dp, c15 = 120.0_dp, &
                         c16 = 127.0_dp, c17 = 140.0_dp, c18 = 175.0_dp, &
                         c19 = 210.0_dp, c20 = 252.0_dp, c21 = 264.0_dp, &
                         c22 = 294.0_dp, c23 = 346.0_dp, c24 = 420.0_dp, &
                         c25 = 462.0_dp, c26 = 606.0_dp, c27 = 672.0_dp, &
                         c28 = 707.0_dp, c29 = 735.0_dp, c30 = 889.0_dp, &
                         c31 = 932.0_dp, c32 = 966.0_dp, c33 = 1141.0_dp, &
                         c34 = 1182.0_dp, c35 = 1278.0_dp, c36 = 1740.0_dp, &
                         c37 = 2520.0_dp, c38 = 5040.0_dp

!       Test arguments and initialise

fn_val = -one
IF (p < pmin .OR. p > pmax) THEN
  WRITE(*, *) 'Error in PPCHI2: p must be between 0.000002 & 0.999998'
  RETURN
END IF
IF (v <= zero) THEN
  WRITE(*, *) 'Error in PPCHI2: Number of deg. of freedom <= 0'
  RETURN
END IF

xx = half * v
c = xx - one

!       Starting approximation for small chi-squared

IF (v < -c5 * LOG(p)) THEN
  fn_val = (p * xx * EXP(g + xx * aa)) ** (one/xx)
  IF (fn_val < e) GO TO 6
  GO TO 4
END IF

!       Starting approximation for v less than or equal to 0.32

IF (v > c3) GO TO 3
fn_val = c4
a = LOG(one-p)

2 q = fn_val
p1 = one + fn_val * (c7+fn_val)
p2 = fn_val * (c9 + fn_val * (c8 + fn_val))
t = -half + (c7 + two * fn_val) / p1 - (c9 + fn_val * (c10 + three * fn_val)) / p2
fn_val = fn_val - (one - EXP(a + g + half * fn_val + c * aa) * p2 / p1) / t
IF (ABS(q / fn_val - one) > c1) GO TO 2
GO TO 4

!       Call to algorithm AS 241 - note that p has been tested above.

3 CALL ppnd16(p, x, if1)

!       Starting approximation using Wilson and Hilferty estimate

p1 = c2 / v
fn_val = v * (x * SQRT(p1) + one - p1) ** 3

!       Starting approximation for p tending to 1

IF (fn_val > c6 * v + six) fn_val = -two * (LOG(one-p) - c * LOG(half * fn_val) + g)

!       Call to algorithm AS 239 and calculation of seven term Taylor series

4 DO i = 1, maxit
  q = fn_val
  p1 = half * fn_val
  p2 = p - gammad(p1, xx)

  t = p2 * EXP(xx * aa + g + p1 - c * LOG(fn_val))
  b = t / fn_val
  a = half * t - b * c
  s1 = (c19 + a * (c17 + a * (c14 + a * (c13 + a * (c12 + c11 * a))))) / c24
  s2 = (c24 + a * (c29 + a * (c32 + a * (c33 + c35 * a)))) / c37
  s3 = (c19 + a * (c25 + a * (c28 + c31 * a))) / c37
  s4 = (c20 + a * (c27 + c34 * a) + c * (c22 + a * (c30 + c36 * a))) / c38
  s5 = (c13 + c21 * a + c * (c18 + c26 * a)) / c37
  s6 = (c15 + c * (c23 + c16 * c)) / c38
  fn_val = fn_val + t * (one + half * t * s1 - b * c * (s1 - b *   &
           (s2 - b * (s3 - b * (s4 - b * (s5 - b * s6))))))
  IF (ABS(q / fn_val - one) > e) RETURN
END DO

WRITE(*, *) 'Error in PPCHI2: Max. number of iterations exceeded'

6 RETURN
END FUNCTION ppchi2

max_stars_repo_name: agforero/FTFramework

max_issues_repo_name: agforero/fortran-testing-framework

max_forks_repo_name: agforero/fortran-testing-framework


text:
  function sinvla (gsmObj, j, a, z, e)

! ======================================================================
!
!     Written in 2001 by K.K. Gudima
!     Edited in 2001 by S.G. Mashnik
!     Modified in October, 2003 by A. J. Sierk
!     Edited by AJS, LANL T-2, December, 2011.
!
! ======================================================================
!
!  Calculates the inelastic cross section (mb) for interaction of the
!  particle "j" of energy e (MeV) with the nucleus (a,z) using NASA
!  systematics (NIM B 117 (1996) 347) for all charged particles and
!  neutrons with energy above the peak in the cross section, and
!  Kalbach systematics (J. Phys. G: Nucl. Phys., 24 (1998) 847) for
!  lower energy neutrons;
!
!  j = 1, 2, 3, 4, 5, 6 means n, p, d, t, He3, and He4
!
! ======================================================================

    use, intrinsic:: iso_fortran_env, only: int32, real64
    use gsm_params, only: zro, one
    use generalizedSpallationData, only: enmax, scalf

    implicit none
    class(GSM),     intent(inout) :: gsmObj
    integer(int32), intent(in   ) :: j
    real(real64),   intent(in   ) :: a
    real(real64),   intent(in   ) :: z
    real(real64),   intent(in   ) :: e
    real(real64)                  :: sinvla

    integer(int32) :: iat
    real(real64)   :: ea, enm, scf

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

    real(real64), dimension( 6 ), parameter :: &
         & ap = [ 1.0_real64, 1.0_real64, 2.0_real64, 3.0_real64, 3.0_real64, 4.0_real64 ], &
         & zp = [ 0.0_real64, 1.0_real64, 1.0_real64, 1.0_real64, 2.0_real64, 2.0_real64 ]

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

    iat = nint(a)
    enm = enmax(iat-1)
    scf = scalf(iat-1)
    if (j > 6) then
       write(gsmObj%io%message, 1100)"42"
       call gsmObj%io%print(3, 3, gsmObj%io%message)
    end if
    if (j == 1 .and. e < enm .and. (iat > 2 .and. iat < 300)) then
!   Kalbach's renormalized parametrization is used only
!   for neutron at e<enm
       sinvla = gsmObj%cskalb (1, zro, one, z, a, e)*scf
    else
!   NASA's parametrization is used for p,d,t,he-3,he-4 and
!   for neutron at e>enm
       ea = e/ap(j)
       sinvla = gsmObj%xabs (zp(j), ap(j), z, a, ea)
    endif
    return

! ======================================================================
1100 format("Exceeded ap() array in 'siginv.f90', line ", A)
! ======================================================================
  end function sinvla

  function xabs (gsmObj, zp, ap, zt, at, e)

! ======================================================================
!
!   Absorption xsec revised version RKT-97/4;
!   neutron data from Barashenkov.
!   This gives absorption xsec for given zp, ap, zt, at, e
!   (Mev/nucleon); it can be used for neutrons also.
!
!> Link to the citation:
!> https://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19980003828.pdf
!> \todo Made a *.bib file and add a \c \cite command to here.
!
!     Modified by A. J. Sierk, LANL T-16, June, 2002.
!     Modified by A. J. Sierk, October, 2003.
!     Edited by AJS, LANL T-2, December, 2011.
!     Edited by LMK, XCP-3, July 2013 (included error protection)
!
! ======================================================================

    use, intrinsic:: iso_fortran_env, only: int32, real64
    use gsm_params, only: zro, one, two, thr, four, pi, twpi, ato3rd

    implicit none
    class(GSM),      intent(inout) :: gsmObj
    real(real64),    intent(in   ) :: zp
    real(real64),    intent(in   ) :: ap
    real(real64),    intent(in   ) :: zt
    real(real64),    intent(in   ) :: at
    real(real64),    intent(in   ) :: e
    real(real64)                   :: xabs

    integer(int32) :: iap, iat, izp, izt
    real(real64)   :: abi, apthrd, atthrd, aux1, aux2, axe, bcm, beta, &
         & bigb, bigr, ce, const, const1, delta, density, dumm1, e20,  &
         & e40, ecm, ecm13, ecmp, ecmt, gcm, plab, rp, rt, sl, t1, temp, &
         & term1, twxsec, vp, vt, x1, xat, xm, xq1, xq2, xq3, xq4, xzt

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

    real(real64), parameter :: c13 = one/thr
    real(real64), parameter :: c43 = four*pi/thr

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

!   Save original particle identities:
    xq1 = zp
    xq2 = ap
    xq3 = zt
    xq4 = at
!   Reverse identities so smaller Z particle is the projectile:
    if (xq3 < xq1) then
       xzt = xq1
       xat = xq2
       xq1 = xq3
       xq2 = xq4
       xq3 = xzt
       xq4 = xat
    endif
    iap = nint(xq2)
    iat = nint(xq4)
    izp = nint(xq1)
    izt = nint(xq3)
!
!  Nucleon-nucleon inelastic xsec not included here:
!
    if (izp*izt == 1 .or. (izp + izt) == 1) then
       xabs = zro
       return
    endif
    rp = gsmObj%radius (iap)
    rt = gsmObj%radius (iat)
    vp =  c43*rp**3
    vt =  c43*rt**3
    if (vp < div0Lim .and. vp > -div0Lim) then
       vp = div0Lim
       write(gsmObj%io%message,1000) "132"
       call gsmObj%io%print(4, 3, gsmObj%io%message)
    end if
    if (vt < div0Lim .and. vt > -div0Lim) then
       vt = div0Lim
       write(gsmObj%io%message,1000) "133"
       call gsmObj%io%print(4, 3, gsmObj%io%message)
    end if
    density = 0.5d0*((xq2/vp) + (xq4/vt))
    const = 1.75d0*density/8.824728d-2
    if (xq4 < div0Lim .and. xq4 > -div0Lim) then
       xq4 = div0Lim
       write(gsmObj%io%message,1000) "139"
       call gsmObj%io%print(4, 3, gsmObj%io%message)
    end if
    abi = one - two*xq3/xq4
    apthrd = ato3rd(iap)
    atthrd = ato3rd(iat)
    if (iap == 1) const = 2.05d0
    if (izp == 2 .and. iap == 4) then
       const1 = 2.77d0 - xq4*8.0d-3 + xq4*xq4*1.8d-5
    endif
    if (izp == 3) const = const/thr
    t1 = 40.d0
    axe = one + e/938.d0
    temp = xq2**2 + xq4**2 + two*xq2*axe*xq4
    if (temp < div0Lim .and. temp > -div0Lim) then
       temp = div0Lim
       write(gsmObj%io%message,1000) "154"
       call gsmObj%io%print(4, 3, gsmObj%io%message)
    end if
    gcm = (xq2*axe + xq4)/sqrt(abs(temp))
    if (gcm < div0Lim .and. gcm > -div0Lim) then
       gcm = div0Lim
       write(gsmObj%io%message,1000) "159"
       call gsmObj%io%print(4, 3, gsmObj%io%message)
    end if
    bcm = sqrt(abs(one - one/gcm**2))
    plab = xq2*sqrt(abs(two*938.d0*e + e*e))
    ecmp = gcm*axe*938.d0*xq2 - bcm*gcm*plab - xq2*938.d0
    ecmt = 938.d0*xq4*(gcm - one)
    ecm = ecmp + ecmt
    ecm13 = ecm**c13
    if (ecm13 < div0Lim .and. ecm13 > -div0Lim) then
       ecm13 = div0Lim
       write(gsmObj%io%message,1000) "169"
       call gsmObj%io%print(4, 3, gsmObj%io%message)
    end if
    bigr = rp + rt + 1.2d0*(apthrd + atthrd)/ecm13
    if (bigr < div0Lim .and. bigr > -div0Lim) then
       bigr = div0Lim
       write(gsmObj%io%message,1000) "174"
       call gsmObj%io%print(4, 3, gsmObj%io%message)
    end if
    bigb = 1.43996517d0*xq1*xq3/bigr
    if (izp == 1 .and. iat > 56) bigb = 0.90d0*bigb
    if (iap == 1) then
!   Proton:  (neutron has bigb = 0.)
       if (iat < 4) then
          bigb = 21.d0*bigb
       elseif (iat == 4) then
          bigb = 27.d0*bigb
       elseif (iat == 12) then
          bigb = 3.5d0*bigb
       elseif (iat <= 16 .and. iat >= 13) then
          bigb = (xq4/7.d0)*bigb
       endif
       if (izt == 12) then
          bigb = 1.8d0*bigb
       elseif (izt == 14) then
          bigb = 1.4d0*bigb
       elseif (izt == 20) then
          bigb = 1.3d0*bigb
       endif
    endif
!  Possible for d or t incident on H
    if (iap < 4 .and. iat == 1) bigb = 21.d0*bigb
    xm = one
    if (izp == 0) then
!  Neutron:
       if (iat >= 11 .and. iat < 40) t1 = 30.d0
       if (izt == 14) t1 = 35.d0
       if (izt == 26) t1 = 30.d0
       if (izt == 0) then
          bigb = zro
       elseif (iat < 200) then
          x1 = 2.83d0 - 3.1d-2*xq4 + 1.7d-4*xq4*xq4
          x1 = max (one, x1)
          if (iat < 12) then
             sl = 0.6d0
          elseif (iat == 12) then
             sl = 1.6d0
          else
             sl = one
          endif
          temp = sl*x1
          if (temp < div0Lim .and. temp > -div0Lim) then
             temp = div0Lim
             write(gsmObj%io%message,1000) "220"
             call gsmObj%io%print(4, 3, gsmObj%io%message)
          end if
          xm = (one - x1*exp(-e/(temp)))
       else
!   7/11/96
          xm = (  one - 0.3d0 * exp( -(e - one)/15.d0 )  ) * &
               & (  one - exp( -(e - 0.9d0) )  )
       endif
    endif
    temp = one + exp((250.d0 - e)/75.d0)
    aux1 = one/(temp)
    if (izp == 2 .and. iap == 4) then
!  Alpha:
       const = const1 - 0.8d0*aux1
    endif
    if (izp == 1 .and. iap == 1) then
!  Proton:
       const = 2.05d0 - 0.05d0*aux1
       if (iat < 4) then
          t1 = 55.d0
          const = 1.7d0
       elseif (iat > 45) then
          t1 = 40.d0 + xq4/thr
       endif
       e20 = exp((20.d0 - e)/10.d0)
       temp = one + e20
       aux2 = e20/(temp)
       if (izt == 12) then
          t1 = 40.d0
          const = 2.05d0 - thr*aux2
       elseif (izt == 14) then
          t1 = 40.d0
          const = 2.05d0 - 1.75d0*aux2
       elseif (izt == 18) then
          t1 = 40.d0
          const = 2.05d0 - two*aux2
       elseif (izt == 20) then
          t1 = 40.d0
          e40 = exp((40.d0 - e)/10.d0)
          temp = one + e40
          const = 2.05d0 - one*e40/(temp)
       endif
    elseif (izp == 0 .and. iap == 1) then
!  Neutron:
       if (density < div0Lim .and. density > -div0Lim) then
          density = div0Lim
          write(gsmObj%io%message,1000) "278"
          call gsmObj%io%print(4, 3, gsmObj%io%message)
       end if
       const = two*0.134457d0/density
       if (iat < 60) then
          const = const - 1.5d0*abi
          ! Adjust cross section for shell corrections (NOTE for large E it's
          ! negligible)
          if (iat <= 40 .and. ((0.01_real64 * e - 1.7_real64) <= 12)) then
             temp = one + exp(-1.70d0 + 0.01d0*e)
             const = const + 0.25d0/(temp)
         end if
       elseif (iat > 140 .and. iat < 200) then
          const = const - 1.5d0*abi
       endif
       temp = xq4 - xq3
       if (temp < div0Lim .and. temp > -div0Lim) then
          temp = div0Lim
          write(gsmObj%io%message,1000) "295"
          call gsmObj%io%print(4, 3, gsmObj%io%message)
       end if
       if (izt > 82) const = const - xq3/(temp)
       if (izt <= 20 .and. izt >= 10) then
          if ((0.1_real64 * e - two) <= 12) then
             temp = one + exp(0.1d0*e - two)
             const = const - one/(temp)
          end if
       elseif (izt >= 82) then
          temp = one + exp(0.05d0*e - one)
          const = const - two/(temp)
       endif
    endif
    if (t1 < div0Lim .and. t1 > -div0Lim) then
       t1 = div0Lim
       write(gsmObj%io%message,1000) "316"
       call gsmObj%io%print(4, 3, gsmObj%io%message)
    end if
    ce = const*(one - exp(-e/t1)) - 0.292d0*exp(-e/792.d0) * &
         & cos(0.229d0*e**0.453d0)
    dumm1 = atthrd + apthrd
    if (dumm1 < div0Lim .and. dumm1 > -div0Lim) then
       dumm1 = div0Lim
       write(gsmObj%io%message,1000) "323"
       call gsmObj%io%print(4, 3, gsmObj%io%message)
    end if
    term1 = atthrd*apthrd/dumm1
    delta = 1.615d0*term1 - 0.873d0*ce
    if (ecm13 < div0Lim .and. ecm13 > -div0Lim) then
       ecm13 = div0Lim
       write(gsmObj%io%message,1000) "329"
       call gsmObj%io%print(4, 3, gsmObj%io%message)
    end if
    delta = delta + 0.140d0*term1/ecm13
    if (xq2 < div0Lim .and. xq2 > -div0Lim) then
       xq2 = div0Lim
       write(gsmObj%io%message,1000) "334"
       call gsmObj%io%print(4, 3, gsmObj%io%message)
    end if
    delta = delta + 0.794d0*abi*xq1/xq2
    delta = -delta
    beta = one
    twxsec = 10.d0*pi*1.26d0*1.26d0*beta*(0.873d0*dumm1 - delta)**2
    if (ecm < div0Lim .and. ecm > -div0Lim) then
       ecm = div0Lim
       write(gsmObj%io%message,1000) "342"
       call gsmObj%io%print(4, 3, gsmObj%io%message)
    end if
    xabs = twxsec*(one - bigb/ecm)*xm
    xabs = max (xabs, zro)

    return

! ======================================================================
1000 format("Divide by zero error prevented in 'siginv.f90', line(s) ", A)
! ======================================================================
  end function xabs

  function radius (gsmObj, ia)

! ======================================================================
!
!     Edited in September, 2003 by A. J. Sierk
!     Edited by AJS, LANL T-2, December, 2011.
!
! ======================================================================

    use, intrinsic:: iso_fortran_env, only: int32, real64
    use gsm_params, only: ato3rd

    implicit none
    class(GSM),     intent(inout) :: gsmObj
    integer(int32), intent(in   ) :: ia
    real(real64)                  :: radius

    integer(int32) :: i

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

    integer(int32), parameter, dimension(23) :: na =             &
         & [ 1,     2,     3,     4,     6,     7,     9,    10, &
         &  11,    12,    13,    14,    15,    16,    17,    18, &
         &  19,    20,    22,    23,    24,    25,    26         ]

    real(real64),   parameter, dimension(23) :: rms =              &
         & [ 0.850d0, 2.095d0, 1.976d0, 1.671d0, 2.570d0, 2.410d0, &
         &   2.519d0, 2.45d0 , 2.420d0, 2.471d0, 2.440d0, 2.580d0, &
         &   2.611d0, 2.730d0, 2.662d0, 2.727d0, 2.900d0, 3.040d0, &
         &   2.969d0, 2.940d0, 3.075d0, 3.110d0, 3.060d0           ]

    real(real64), parameter :: fact = sqrt( 5.0_real64 / 3.0_real64 )

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

    if (ia == 5 .or. ia == 8 .or. ia == 21 .or. ia > 26) then
       radius = fact * (0.84d0*ato3rd(ia)  + 0.55d0)
    else
       do i = 1,23
          if (ia == na(i)) radius = fact*rms(i)
       end do
    endif

    return

! ======================================================================
  end function radius

  function cskalb (gsmObj, kp, zp, ap, zt, at, elab)

! ======================================================================
!
!     Written in 1982; revised 1990
!
!     Calculate optical model reaction cross sections
!     using empirical parameterization
!     of Narasimha Murthy, Chaterjee, and Gupta
!     going over to the geometrical limit at high energy
!
!           Proton cross sections scaled down with signor for a<100
!           (appropriate for Becchetti-Greenlees potential)
!           Neutron cross sections scaled down sith signor for a<40
!           (appropriate for Mani et al potential)
!
!     parameter values set in subroutine sigpar
!
!     Edited in September, 2003 by A. J. Sierk
!     Edited by AJS, LANL T-2, December, 2011.
!     Edited by LMK, XCP-3, July 2013 (included error protection)
!
! ======================================================================

    use, intrinsic:: iso_fortran_env, only: int32, real64
    use gsm_params, only: zro, one, two, ato3rd
    use generalizedSpallationData, only: xl0, xl1, xm0, xm1, xn0, xn1, &
         & xn2, xp0, xp1, xp2

    implicit none
    class(GSM),     intent(inout) :: gsmObj
    integer(int32), intent(in   ) :: kp
    real(real64),   intent(in   ) :: zp
    real(real64),   intent(in   ) :: ap
    real(real64),   intent(in   ) :: zt
    real(real64),   intent(in   ) :: at
    real(real64),   intent(in   ) :: elab
    real(real64)                  :: cskalb

    integer(int32) :: iat, jout
    real(real64)   :: a, ares, athrd, b, c, cut, ec, ecsq, ecut, ecut2, &
         & etest, flow, geom, p, ra, sig, signor, signor2, spill, &
         & temp, w, xlamb, xmu, xnu, xnulam

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

    flow = 1.d-250
    spill = 1.d+250
    jout = nint (ap)
    iat = nint (at)
    ares = at + ap
    athrd = ato3rd(iat + jout)
    signor = one
!    signor reduces p and n result for light targs as per expt.
    if (kp == 1) then
       if (ares < 40.d0) signor = 0.7d0 + ares*0.0075d0
       xlamb = xl0(1)/athrd + xl1(1)
       xmu = xm0(1)*athrd + xm1(1)*athrd*athrd
       xnu = xn0(1)*athrd*ares + xn1(1)*athrd*athrd + xn2(1)
       ec = 0.5d0
       ecsq = 0.25d0
       p = xp0(1)
       xnulam = one
       etest = 32.d0
!     etest is the energy above which the rxn cross section is
!     compared with the geometrical limit and the max taken.
!     xnulam here is a dummy value to be used later.
       ra = zro
    else
       ra = 1.20d0
       if (kp == 2) then
          ra = zro
          if (ares < 60.d0) then
             signor = 0.92d0
          elseif (ares < 100.d0) then
             signor = 0.8d0 + ares*0.002d0
          endif
       endif
       temp = 1.5d0*athrd + ra
       ec = 1.439965173d0*zp*zt/(temp)
       if (ec < div0Lim .and. ec > -div0Lim) then
          ec = div0Lim
          write(gsmObj%io%message,1000) "493"
          call gsmObj%io%print(4, 3, gsmObj%io%message)
       end if
       ecsq = ec * ec
       p = xp0(kp) + xp1(kp)/ec + xp2(kp)/ecsq
       xlamb = xl0(kp)*ares + xl1(kp)
       a = ares**xm1(kp)
       xmu = xm0(kp) * a
       xnu = a*(xn0(kp) + xn1(kp)*ec + xn2(kp)*ecsq)
       if (jout == 2 .or. jout == 3) ra = 0.8d0
!    new values of ra are for calculating the geometrical limit
!    to the cross section.
       if (kp == 2) then
          c = min(3.15d0, ec*0.5d0)
          w = 0.7d0*c/3.15d0
!    c and w are for the global corr'n factor for elab<ec
!    for light targs they are scaled down from global values
       endif
       if (xlamb < div0Lim .and. xlamb > -div0Lim) then
          xlamb = div0Lim
          write(gsmObj%io%message,1000) "511"
          call gsmObj%io%print(4, 3, gsmObj%io%message)
       end if
       xnulam = xnu/xlamb
       if (xnulam > spill) xnulam = zro
       if (xnulam >= flow) then
          if (kp == 2) then
             etest = sqrt(xnulam) + 7.d0
          else
             etest = 1.2d0*sqrt(xnulam)
          endif
! for xnulam > 0, sig reaches a maximum at sqrt(xnulam).
       endif
    endif
    a = -two*p*ec + xlamb - xnu/ecsq
    b = p*ecsq + xmu + two*xnu/ec
    ecut = zro
    cut = a*a - 4.d0*p*b
    if (cut > zro) ecut = sqrt(abs(cut))
    ecut = (ecut - a)/(p + p)
    ecut2 = ecut
    if (cut < zro) ecut2 = ecut - two
    sig = zro
    if (elab <= ec) then
       if (elab > ecut2) then
          sig = (p*elab*elab + a*elab + b) * signor
          if (kp == 2) then
             if (w < div0Lim .and. w > -div0Lim) then
                w = div0Lim
                write(gsmObj%io%message,1000) "539"
                call gsmObj%io%print(4, 3, gsmObj%io%message)
             end if
             signor2 = (ec - elab - c)/w
             signor2 = one + exp(signor2)
             sig = sig/signor2
          endif
!    first signor gives empirical global corr'ns at low elab
!    second signor corrects values near elab=0; light nuclei
       endif
    else
       sig = (xlamb*elab + xmu + xnu/elab) * signor
       geom = zro
       if (xnulam >= flow .and. elab >= etest) then
          geom = sqrt(abs(ap*elab))
          if (geom < div0Lim .and. geom > -div0Lim) then
             geom = div0Lim
             write(gsmObj%io%message,1000) "563"
             call gsmObj%io%print(4, 3, gsmObj%io%message)
          end if
          geom = 1.23d0*athrd + ra + 4.573d0/geom
          geom = 31.416d0*geom*geom
          sig = max(geom, sig)
       endif
    endif

    cskalb = sig
    return
! ======================================================================
1000 format("Divide by zero error prevented in 'siginv.f90', line(s) ", A)
! ======================================================================
  end function cskalb

max_stars_repo_name: lanl/generalized-spallation-model

max_stars_repo_licenses: Intel Unlicense

max_issues_repo_name: lanl/generalized-spallation-model

max_issues_repo_licenses: Intel Unlicense

max_forks_repo_name: lanl/generalized-spallation-model

max_forks_repo_licenses: Intel Unlicense


module precision

  use iso_fortran_env, only: real32, real64, int32

  integer(int32), parameter :: ip = int32
  integer(ip), parameter :: sp = real32
  integer(ip), parameter :: dp = real64

end module precision


! { dg-do run }
    integer (kind = 4) :: a, a2
    integer (kind = 2) :: b, b2
    real :: c
    double precision :: d, d2, c2
    integer, dimension (10) :: e
    e(:) = 5
    e(7) = 9
!$omp atomic write seq_cst
    a = 1
!$omp atomic seq_cst, write
    b = 2
!$omp atomic write, seq_cst
    c = 3
!$omp atomic seq_cst write
    d = 4
!$omp atomic capture seq_cst
    a2 = a
    a = a + 4
!$omp end atomic
!$omp atomic capture, seq_cst
    b = b - 18
    b2 = b
!$omp end atomic
!$omp atomic seq_cst, capture
    c2 = c
    c = 2.0 * c
!$omp end atomic
!$omp atomic seq_cst capture
    d = d / 2.0
    d2 = d
!$omp end atomic
    if (a2 .ne. 1 .or. b2 .ne. -16 .or. c2 .ne. 3 .or. d2 .ne. 2) call abort
!$omp atomic read seq_cst
    a2 = a
!$omp atomic seq_cst, read
    c2 = c
    if (a2 .ne. 5 .or. b2 .ne. -16 .or. c2 .ne. 6 .or. d2 .ne. 2) call abort
    a2 = 10
    if (a2 .ne. 10) call abort
!$omp atomic capture
    a2 = a
    a = e(1) + e(6) + e(7) * 2
!$omp endatomic
    if (a2 .ne. 5) call abort
!$omp atomic read
    a2 = a
!$omp end atomic
    if (a2 .ne. 28) call abort
!$omp atomic capture seq_cst
    b2 = b
    b = e(1) + e(7) + e(5) * 2
!$omp end atomic
    if (b2 .ne. -16) call abort
!$omp atomic seq_cst, read
    b2 = b
!$omp end atomic
    if (b2 .ne. 24) call abort
end


text:
program random_integers
    implicit none
    integer :: nr_values, lower, upper, i

    call get_values(nr_values, lower, upper)
    do i = 1, nr_values
        print '(I0)', random_integer(lower, upper)
    end do

contains

    !> @brief copmute a random number in a given range
    !>
    !> The function returns a random integer value between a given lower and
    !> upper limits, inclusively.
    !>
    !> @param[in] lower lower value for the random value
    !> @param[in] upper upper value for the random value
    !> @returns integer value
    function random_integer(lower, upper) result(rand_int)
        implicit none
        integer, value :: lower, upper
        integer :: rand_int
        real :: rand_real

        call random_number(rand_real)
        rand_int = lower + floor(rand_real*(upper - lower + 1))
    end function random_integer

    subroutine get_values(nr_values, lower, upper)
        use, intrinsic :: iso_fortran_env, only : error_unit
        implicit none
        integer, intent(out) :: nr_values, lower, upper

        if (command_argument_count() /= 3) then
            write (unit=error_unit, fmt='(A)') &
                'error: argument nr_values lower upper expected'
            stop 1
        end if
        nr_values = get_integer(1)
        lower = get_integer(2)
        upper = get_integer(3)
        if (lower > upper) then
            write (unit=error_unit, fmt='(A)') &
                'error: upper should be larger or equal to lower'
            stop 3
        end if
    end subroutine get_values

    function get_integer(i) result(value)
        use, intrinsic :: iso_fortran_env, only : error_unit
        implicit none
        integer, value :: i
        integer :: value
        character(len=32) :: buffer
        integer :: istat
        character(len=1024) :: msg

        call get_command_argument(i, buffer)
        read (buffer, fmt=*, iostat=istat, iomsg=msg) value
        if (istat /= 0) then
            write (unit=error_unit, fmt='(2A)') 'error: ', msg
            stop 2
        end if
    end function get_integer

end program random_integers

max_stars_repo_name: caguerra/Fortran-MOOC

max_issues_repo_name: caguerra/Fortran-MOOC

max_forks_repo_name: caguerra/Fortran-MOOC


text:
use plantfem
implicit none

type(IO_) :: f
type(Math_) :: math
type(Random_) :: random
complex(complex64),allocatable :: t(:), Ft(:),Fw(:),Fw_bar(:),Ft_(:),omega(:)
real(real64) :: domega, n
integer(int32) :: i

! Time
t = linspace([0.0d0, 10.0d0],1024)

! Frequency
omega = linspace([0.0d0, 1024.0d0/10.0d0],1024)

Ft = zeros(1024)
!Ft(:) = 1.0d0/2.0d0/Math%PI*domega*(sin( (2.0d0*n + 1)/2.0d0*domega*t ) )/( sin( 1.0d0/2.0d0*domega*t  ))
do i=1,size(Ft)
    Ft(i) = random%gauss(mu=0.0d0,sigma=0.50d0)
enddo
Ft(:) = Ft(:)+5.0d0*sin(2.0d0*Math%PI*t)

Fw = FFT(Ft)

Fw( 516-400:516+400 )=0.0d0

Ft_ = IFFT(Fw)

Fw_bar = Fw
Fw_bar = Fw*Fw

call f%open("Fw.txt","w")
!call f%write(real(t), real(Fw) )
call f%write(real(omega(1:size(t)/2 ) ), real(Fw(1:size(t)/2) ) )
call f%close()

call f%open("Ft.txt","w")
!call f%write(real(t), real(Fw) )
call f%write(real(t), real(Ft) )
call f%close()

call f%open("Ft_.txt","w")
!call f%write(real(t), real(Fw) )
call f%write(real(t), real(Ft_) )
call f%close()

call f%plot("Ft.txt","w l")
call f%plot("Fw.txt","w l")
call f%plot("Ft_.txt","w l")
end

max_stars_repo_name: kazulagi/plantFEM


text:
      SUBROUTINE GFSCOM(AWY,NUY,KC,IDENT,AC,SCR)
C
C     ROUTINE TO COMPUTE THE FLUID COMPRESSIBILTY MATRIX
C
C     THIS MATRIX CONTAINS THE SPRING FACTOR WHICH COUPLES THE
C     STRUCTURE AND FREE SURFACE TO PREVENT VOLUME CHANGES
C
      DOUBLE PRECISION        DZ(1)    ,DKCOMP   ,VAL
C
      REAL          KCOMP    ,RZ(1)
C
      INTEGER       KC       ,AC       ,Z        ,SYSBUF   ,AWY
     1             ,MCB(7)   ,NAME(2)  ,TI1      ,TO1      ,TO2
     2             ,SCR
C
C
C     MODULE PARAMETERS
C
      COMMON /BLANK/     NOGRAV   ,NOFREE   ,KCOMP   ,COMPTP
     1                       ,FORM     ,LMODES
C
C     OPEN CORE
C
      COMMON / ZZZZZZ /        Z(1)
C
C     SYSTEM COMMON
C
      COMMON / SYSTEM /       SYSBUF
C
C     PACK - UNPACK COMMON BLOCKS
C
      COMMON / PACKX /        TI1      ,TO1      ,I1       ,N1
     1                       ,INCR1
      COMMON / UNPAKX /       TO2      ,I2       ,N2       ,INCR2
      COMMON / ZBLPKX /       A(4)     ,IROW
C
      EQUIVALENCE   ( Z(1) , RZ(1) , DZ(1) )
     1             ,( VAL , A(1) )
C
      DATA NAME / 4HGFSC , 4HOM   /
C
C     ALLOCATE CORE
C
      NZ = KORSZ(Z(1))
      IBUF = NZ - SYSBUF
      NZ = IBUF - 1
      IF(NZ .LT. NUY) GO TO 1008
C
C     FORM A COLUMN VECTOR OF ONES
C
      TI1 = 1
      TO1 = 2
      I1 = 1
      N1 = NUY
      INCR1 = 1
      DO 30 I=1,NUY
   30 RZ(I) = 1.0
      CALL MAKMCB(MCB,IDENT,NUY,2,2)
      CALL GOPEN(IDENT,Z(IBUF),1)
      CALL PACK(RZ(1),IDENT,MCB)
      CALL CLOSE(IDENT,1)
      CALL WRTTRL(MCB)
C
      CALL SSG2B(AWY,IDENT,0,AC,0,2,1,SCR)
C
C     PERFORM MULTIPLY TO GET COMPRESSIBLITY MATRIX
C
C
C     UNPACK ROW OF AC INTO CORE
C
      MCB(1) = AC
      CALL RDTRL(MCB)
      NROW = MCB(3)
      IF(NZ .LT. 2*NROW) GO TO 1008
      TO2 = 2
      I2 = 1
      N2 = NROW
      INCR2 = 1
C
      CALL GOPEN(AC,Z(IBUF),0)
      CALL UNPACK(*40,AC,DZ(1))
      GO TO 60
C
C     AC IS NULL
C
   40 DO 50 I=1,NROW
   50 DZ(I) = 0.0D0
C
C     SET UP TO CREATE KC MATRIX
C
   60 CALL CLOSE(AC,1)
C
      DKCOMP = DBLE(KCOMP)
      CALL GOPEN(KC,Z(IBUF),1)
      CALL MAKMCB(MCB,KC,NROW,1,2)
C
C     LOOP OVER NON-ZERO TERMS OF AC TO CREATE KC
C
      DO 90 I=1,NROW
      CALL BLDPK(2,2,KC,0,0)
      IF(DZ(I) .EQ. 0.0D0) GO TO 80
      DO 70 J=1,NROW
      IROW = J
      VAL = DKCOMP * DZ(J) * DZ(I)
      CALL ZBLPKI
   70 CONTINUE
   80 CALL BLDPKN(KC,0,MCB)
   90 CONTINUE
      CALL CLOSE(KC,1)
C
C     WRITE TRAILER
C
      CALL WRTTRL(MCB)
      RETURN
C
C     ERRORS
C
 1008 CALL MESAGE(-8,0,NAME)
      RETURN
      END