text:
!  cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!  L. Zhang
!  Northwestern University
!  This subroutine solves for the residual for all degrees of freedom
!____________________________________________________________________
!  L. Zhang, 06/24/2004
!  Tulane University
!  Revised the subroutine to array
!  cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine blockgmresnew(xloc, dloc, doloc, p, hk, ien, f_fluids,ne_local,ien_local,node_local,nn_local, &
			fden,fvis,I_fluid,rngface)
  use global_constants
  use run_variables
  use fluid_variables
  use solid_variables, only: nn_solid
  use r_common, only: density_solid, vis_solid
  implicit none

  integer ien(nen,ne)
  real* 8 xloc(nsd,nn)
  real* 8 dloc(ndf,nn),doloc(ndf,nn)
  real* 8 p(ndf,nn),hk(ne)

  real* 8 x(nsd,nen)
  real* 8 d(ndf,nen),d_old(ndf,nen)

  real* 8 eft0,det,effd,effm,effc
  real* 8 sh(0:nsd,nen),ph(0:nsd,nen)
  real* 8 xr(nsd,nsd),cf(nsd,nsd),sx(nsd,nsd)

  real* 8 drt(ndf),drs(ndf)
  real* 8 dr(nsd,ndf)
  real* 8 u,v,w,pp,ug
  real* 8 tau(nsd,nsd)
  real* 8 hg,taum,tauc,vel,ree, taul
  real* 8 res_c,res_a(nsd),res_t(nsd)
  real* 8 prs_c,prs_t(nsd),p_vec(3),prs_cc(nsd)
  real* 8 mu,nu,ro,g(nsd)
  real* 8 tempc(ndf),temp
  real* 8 dtinv,oma,ama
  integer inl, ie, isd, iq, node,jsd
  integer ieface,irng !,inface

  real* 8 f_fluids(nsd,nn)
  real* 8 fnode(nsd,nen),fq(nsd)
  integer nn_local
  integer node_local(nn_local)
!-------------------------------------------
  real(8) fden(nn)
  real(8) local_den(nen)
  real(8) fvis(nn)
  real(8) local_vis(nen)
  real(8) I_fluid(nn)
  real(8) kappa
  integer rngface(neface,ne)
!---------------------------------------------
!============================
! MPI varibalbes
  integer ne_local ! # of element on each processor
  integer ien_local(ne_local) ! subregion-->wholeregion element index
  integer ie_local ! loop parameter
  integer icount
!---------------------------------------------

  dtinv = 1.0/dt
  if(steady) dtinv = 0.0
  oma   = 1.0 - alpha
  ama   = 1.0 - oma
  kappa = 1.0e4
 !=================================================
!f_fluids(:,:)=f_fluids(:,:)/(0.0625/6.0)
!dloc(ndf,:)=(1.0 - I_fluid(:)) * dloc(ndf,:)
do icount=1, nn_local
        node=node_local(icount)
p(1:nsd,node)=p(1:nsd,node)+f_fluids(1:nsd,node)
end do

!=================================================
!===================================================
! f_fluids is actually the FSI force at fluid nodes,
! then we will just subscrib it from p(!:nsd) which is the
! residuals for the momentum equations.
! 1 let fnod=0 then fndoe and fq 's effects will be canceled out
! 2 do the subscribition after the elements loop
! Xingshi 09/15/2008
!===================================================
  do ie_local=1,ne_local		! loop over elements
      ie=ien_local(ie_local)
      do inl=1,nen
	     x(1:nsd,inl) = xloc(1:nsd,ien(inl,ie))
!============================================================================
!		 fnode(1:nsd,inl) = f_fluids(1:nsd,ien(inl,ie))
                 fnode(:,inl)=0.0
!============================================================================
		 d(1:ndf,inl) =  dloc(1:ndf,ien(inl,ie))
		 d_old(1:ndf,inl) = doloc(1:ndf,ien(inl,ie))
!-----------------------------------------------------
               local_den(inl)=fden(ien(inl,ie))
		local_vis(inl)=fvis(ien(inl,ie))
	 enddo

	 hg = hk(ie)

	 do iq=1,nquad  ! loop over the quadrature points in each element
!...  calculate the shape function and the weight at quad point
		if (nsd==2) then
		    if (nen.eq.3) then !calculate shape function at quad point
			   include "sh2d3n.h"
			elseif (nen.eq.4) then
				include "sh2d4n.h"
			endif
		elseif (nsd==3) then
		    if (nen.eq.4) then !calculate shape function at quad point
			   include "sh3d4n.h"
			elseif (nen.eq.8) then
				include "sh3d8n.h"
			endif
		endif
!		stop
!	write(*,*) 'shape functions=',sh(0,1),sh(0,2),sh(0,3),sh(0,4)
	    eft0 = abs(det) * wq(iq) ! calculate the weight at each quad pt
!...  initialize d, dd/dx, dd/dy, dd/dz, dd/dt
		drs(:) = 0.0
		drt(:) = 0.0
		dr(1:nsd,1:ndf)=0.0
		fq(:)=0.0

!... calculate vi, dvi/dxj
        do inl=1,nen
		   tempc(1:nsd) = ama*d(1:nsd,inl)+oma*d_old(1:nsd,inl)
		   drs(1:nsd) = drs(1:nsd)+sh(0,inl)*tempc(1:nsd)
		   do isd=1,nsd
			 dr(isd,1:nsd) = dr(isd,1:nsd)+sh(isd,inl)*tempc(1:nsd)
		   enddo
		   fq(:) = fq(:) + sh(0,inl)*fnode(:,inl)
	    enddo

        ro=0.0
	mu=0.0
!... calculate dvi/dt, p, dp/dxi
        do inl=1,nen
		   drt(1:nsd)=drt(1:nsd)+sh(0,inl)*(d(1:nsd,inl)-d_old(1:nsd,inl))*dtinv
		   drs(pdf)=drs(pdf)+sh(0,inl)*d(pdf,inl)
		   dr(1:nsd,pdf)=dr(1:nsd,pdf)+sh(1:nsd,inl)*d(pdf,inl)
!----------------------------------------------------------------------------------------
           ro=ro+sh(0,inl)*local_den(inl)
		   mu=mu+sh(0,inl)*local_vis(inl)
	    enddo

!... define u=v1, v=v2, w=v3, pp=p
		if (nsd==2) then
		    u = drs(udf)
			v = drs(vdf)
			w = 0.0
		elseif (nsd==3) then
			u = drs(udf)
		    v = drs(vdf)
		    w = drs(wdf)
		endif
	    pp= drs(pdf)

	    if(stokes) then ! if stokes flow
		   u = 0.0
		   v = 0.0
		   w = 0.0
	    endif

!....  calculate liquid constant and gravity
		g(1:nsd)  = gravity(1:nsd)  ! gravatitional force

	! believe nu is calculated only for turbulent model
		if (nsd==2) then
	    nu = delta(4)*turb_kappa**2*hg**2 * sqrt(2*dr(1,1)**2+(dr(2,1)+dr(1,2))**2 &
	                                            +2*dr(2,2)**2)
		elseif (nsd==3) then
	    nu = delta(4)*turb_kappa**2*hg**2 * sqrt(2*dr(1,1)**2+(dr(2,1)+dr(1,2))**2 &
	                                            +2*dr(2,2)**2+(dr(3,1)+dr(1,3))**2 &
	                                            +2*dr(3,3)**2+(dr(3,2)+dr(2,3))**2)
		endif
	    mu = mu + nu*ro

!....  calculate each term in the residual equation
	    res_c = 0.0
		if (nsd==2) then
		  do inl=1,nen
		   	 res_c = res_c+sh(xsd,inl)*d(udf,inl) &
	                    +sh(ysd,inl)*d(vdf,inl)
		  enddo
		elseif (nsd==3) then
		  do inl=1,nen
		     res_c = res_c+sh(xsd,inl)*d(udf,inl) &
	                    +sh(ysd,inl)*d(vdf,inl) &
	                    +sh(zsd,inl)*d(wdf,inl)
		  enddo
		endif

        do inl=1,nen
		   node=ien(inl,ie)
		   res_c=res_c+sh(0,inl)*(d(ndf,inl)-d_old(ndf,inl))*dtinv* &
		   (1.0/kappa*I_fluid(node))
		enddo  ! add dp/dt term for artificial fluid

	    do isd = 1, nsd
			if (nsd==2) then
			   res_a(isd)=ro*(drt(isd)+u*dr(1,isd)+v*dr(2,isd)-g(isd))-fq(isd)

			elseif (nsd==3) then
			   res_a(isd)=ro*(drt(isd)+u*dr(1,isd)+v*dr(2,isd)+w*dr(3,isd)-g(isd))-fq(isd)
			endif
	    enddo

		res_t(1:nsd) = dr(1:nsd,pdf) + res_a(1:nsd)

		!  ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        ! Mickael 02/01/2005
        ! TAUm, TAUc and TAUl (l for  lsic), See Tezduyar and Sathe, 2003

        vel  = sqrt(u*u+v*v+w*w)  !magnitude of the velocity

        ree  = vel*hg*ro/(2.0*mu)  !Reynolds number
        if(steady.or.(.not.taudt)) then !stablization, taum
			taum = 1.0/sqrt((2.0*vel/hg)**2+(4.0*mu/hg/hg)**2)
		else
			taum = 1.0/sqrt((2.0/dt)**2+(2.0*vel/hg)**2+(4.0*mu/hg/hg)**2)
        endif
        tauc = taum/ro

        if (ree.le.3.0) then
			taul = hg*vel*ree/6.0
        else
            taul = hg*vel/2.0
        endif
        !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!.....   Density optimization

		ph(0:nsd,1:nen) = sh(0:nsd,1:nen)*eft0

!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!c.....   Galerkin Terms (Look at notes)
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!c....   Calculate stress tau(ij) and pressure
		do isd = 1,nsd
			do jsd = 1,nsd
				tau(isd,jsd) = mu*(dr(isd,jsd) + dr(jsd,isd))
			enddo
		enddo

		prs_t(1:nsd) = res_t(1:nsd)*taum

		!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        ! Mickael 02/01/2005
        ! TAUl (l for  lsic), See Tezduyar and Sathe, 2003

        prs_c = ro*res_c*taul
        !
        !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

	    prs_cc(1:nsd) = res_t(1:nsd)*tauc
!.... calculate the residual at each degree of freedom
	    do inl=1,nen ! loop over number of nodes in an element

		   node = ien(inl,ie)
			if (nsd==2) then
			   temp = ro*(u*ph(xsd,inl)+v*ph(ysd,inl))
			elseif (nsd==3) then
			   temp = ro*(u*ph(xsd,inl)+v*ph(ysd,inl)+w*ph(zsd,inl))
			endif
		! Continuty Equation
		   p(pdf,node) = p(pdf,node)-ph(0,inl)*res_c

		! Momentum Equation (Euler Residual)
!===========================================================
! why originally it is minus?
		   p(1:nsd,node) = p(1:nsd,node)-ph(0,inl)*res_a(1:nsd)
!=============================================================
		! Momentum Equation (C: Diffusion terms)
		    if (nsd==2) then
			  do isd=1,nsd
				p(isd,node)=p(isd,node) + ph(isd,inl)*pp -   &
										  ph(1,inl)*tau(1,isd) -  &
										  ph(2,inl)*tau(2,isd)
				p(isd,node)=p(isd,node)+mu*ph(isd,inl)*(dr(1,1)+dr(2,2))*2.0/3.0
			  enddo
			elseif (nsd==3) then
			  do isd=1,nsd
				p(isd,node)=p(isd,node) + ph(isd,inl)*pp -   &
										  ph(1,inl)*tau(1,isd) -  &
										  ph(2,inl)*tau(2,isd) -  &
										  ph(3,inl)*tau(3,isd)
				p(isd,node)=p(isd,node)+mu*ph(isd,inl)*(dr(1,1)+dr(2,2)+dr(3,3))*2.0/3.0
			  enddo
			endif

		! Stablization with Tau_moment
		   if (nsd==2) then
		   	 p(pdf,node) = p(pdf,node) - ph(xsd,inl)*prs_cc(udf)  &
	                                   - ph(ysd,inl)*prs_cc(vdf)
	       elseif (nsd==3) then
		     p(pdf,node) = p(pdf,node) - ph(xsd,inl)*prs_cc(udf)  &
	                                   - ph(ysd,inl)*prs_cc(vdf)  &
	                                   - ph(zsd,inl)*prs_cc(wdf)
		   endif		! Stablization with Tau_cont
		   p(1:nsd,node) = p(1:nsd,node) - prs_t(1:nsd)*temp - ph(1:nsd,inl)*prs_c
	 enddo

	 enddo ! end of qudrature pts loop
  enddo ! end of element loop
 continue
continue
  return
end subroutine blockgmresnew

max_stars_repo_name: biofluids/IFEM-archive

max_issues_repo_name: biofluids/IFEM-archive

max_forks_repo_name: biofluids/IFEM-archive


! { dg-do compile }
! { dg-options "-O3 -ffast-math -fdump-tree-optimized" }

! Make sure that FRE does not replace c with b in d = c - 5

function test(a)
  real a, b, c, d
  b = a + 5.
  c = (a + 5.)
  d = c - 5.
  call foo(b)
  test = d
end

! { dg-final { scan-tree-dump "- 5" "optimized" } }


text:
subroutine foo_ivybridge
end subroutine

max_stars_repo_name: dvuckovic/ecbuild

max_issues_repo_name: dvuckovic/ecbuild

max_forks_repo_name: dvuckovic/ecbuild


text:
MODULE module_sf_noahutl

  REAL, PARAMETER      :: CP = 1004.5, RD = 287.04, SIGMA = 5.67E-8,    &
                          CPH2O = 4.218E+3,CPICE = 2.106E+3,            &
                          LSUBF = 3.335E+5

CONTAINS

  SUBROUTINE CALTMP(T1, SFCTMP, SFCPRS, ZLVL, Q2, TH2, T1V, TH2V, RHO )

    IMPLICIT NONE

    ! Input:
    REAL, INTENT(IN)       :: T1     ! Skin temperature (K)
    REAL, INTENT(IN)       :: SFCTMP ! Air temperature (K) at level ZLVL
    REAL, INTENT(IN)       :: Q2     ! Specific Humidity (kg/kg) at level ZLVL
    REAL, INTENT(IN)       :: SFCPRS ! Atmospheric pressure (Pa) at level ZLVL
    REAL, INTENT(IN)       :: ZLVL   ! Height (m AGL) where atmospheric fields are valid

    ! Output:
    REAL, INTENT(OUT)      :: TH2    ! Potential temperature (considering the reference pressure to be at the surface)
    REAL, INTENT(OUT)      :: T1V    ! Virtual skin temperature (K)
    REAL, INTENT(OUT)      :: TH2V   ! Virtual potential temperature at ZLVL
    REAL, INTENT(OUT)      :: RHO    ! Density

    ! Local:
    REAL                   :: T2V

    TH2 = SFCTMP + ( 0.0098 * ZLVL)
    T1V= T1 * (1.0+ 0.61 * Q2)
    TH2V = TH2 * (1.0+ 0.61 * Q2)
    T2V = SFCTMP * ( 1.0 + 0.61 * Q2 )
    RHO = SFCPRS/(RD * T2V)

  END SUBROUTINE CALTMP

  SUBROUTINE CALHUM(SFCTMP, SFCPRS, Q2SAT, DQSDT2)

    IMPLICIT NONE

    ! Input:
    REAL, INTENT(IN)       :: SFCTMP
    REAL, INTENT(IN)       :: SFCPRS

    ! Output:
    REAL, INTENT(OUT)      :: Q2SAT   ! Saturated specific humidity
    REAL, INTENT(OUT)      :: DQSDT2

    ! Local
    REAL, PARAMETER        :: A2=17.67,A3=273.15,A4=29.65, ELWV=2.501E6,        &
                              A23M4=A2*(A3-A4), E0=611.0, RV=461.0,             &
                              EPSILON=0.622
    REAL                   :: ES

    ! ES:  e.g. Dutton chapter 8, eq 11
        ES = E0 * EXP ( ELWV/RV*(1./A3 - 1./SFCTMP) )

    ! Q2SAT:
        Q2SAT = EPSILON * ES / (SFCPRS-(1-EPSILON)*ES)

! DQSDT2 is calculated assuming Q2SAT is a specific humidity
        !KWM DQSDT2=(Q2SAT/(1+Q2SAT))*A23M4/(SFCTMP-A4)**2
        DQSDT2=Q2SAT*A23M4/(SFCTMP-A4)**2

        END SUBROUTINE CALHUM

END MODULE module_sf_noahutl

subroutine wrf_error_fatal(string)
  implicit none
  character(len=*), intent(in) :: string
  print*, string
  stop
end subroutine wrf_error_fatal

subroutine wrf_message(msg)
  implicit none
  character(len=*), intent(in) :: msg
  write(*,'(A)') msg
end subroutine wrf_message

logical function wrf_dm_on_monitor() result (return_value)
  implicit none
  return_value = .TRUE.
end function wrf_dm_on_monitor

subroutine wrf_dm_bcast_real(rval, ival)
  implicit none
  real, intent(in) :: rval
  integer, intent(in) :: ival
end subroutine wrf_dm_bcast_real

subroutine wrf_dm_bcast_integer(ival1, ival2)
  implicit none
  real, intent(in) :: ival1
  integer, intent(in) :: ival2
end subroutine wrf_dm_bcast_integer

subroutine wrf_dm_bcast_string(sval, ival)
  implicit none
  character(len=*), intent(in) :: sval
  integer, intent(in) :: ival
end subroutine wrf_dm_bcast_string

max_stars_repo_name: KathyNie/LISF

max_issues_repo_name: dmocko/LISF

max_forks_repo_name: dmocko/LISF


text:
C
C
*+ do_2p_lookup

       subroutine do_2p_lookup(status)
C      -------------------------------
C
C Construct a two-point spectral index vs X table for the current spectrum
C
C Returned:
C     status
       integer      status
C
C A table of spectral-index as a function of X (X = sqrt(f1/f_break))
C is constructed for the supplied frequencies f1, f2.  This may be
C plotted.  The lookup table is available for future fitting of spectral
C index results.
*-
       include '../include/anmap_sys_pars.inc'
       include '/mrao/include/chrlib_functions.inc'

C local variables
       real*4      f1, f2, x1, x2, x, y
       real*4      synch_2p_val
       integer     i
       integer     number_line_segments, iunit
       parameter  (number_line_segments = 100)
       save        f1, f2

C check status on entry
       if (status.ne.0) return

C prompt user for frequencies
       if (f1.le.0.0) f1 = 1.4
       if (f2.le.0.0) f2 = 5.0
       call io_getr('Lower-frequency (GHz) : ','*',f1,status)
       call io_getr('Upper-frequency (GHz) : ','*',f2,status)

C make lookup table
       call synch_2p_make( f1, f2, status )
       if (status.ne.0) goto 999

C open results file
       call io_opefil(iunit,general_results_file,'write',0,status)
       write(iunit,*)'%ndata ',number_line_segments
       write(iunit,*)'%ncols   2'

C plot function
       x = x1
       y = synch_2p_val(x)
       do i = 1,number_line_segments
         x = x1 + i*(x2-x1)/float(number_line_segments)
         y = synch_2p_val(x)
         write(iunit,*) x, y
       end do

C close results file
       close (iunit)

C error report
999    call iocmd_err(status,'MAKE-ALPHA-X-TABLE',' ')

       end

ext: f

max_stars_repo_name: CavendishAstrophysics/anmap

max_issues_repo_name: CavendishAstrophysics/anmap

max_forks_repo_name: CavendishAstrophysics/anmap


text:
MODULE test_mod
    implicit none
    DOUBLE PRECISION, POINTER :: concen(:,:)
    INTEGER :: nsize

END MODULE test_mod

SUBROUTINE set_ptr(fc_concen,n)
    USE iso_c_binding
    USE test_mod

    DOUBLE PRECISION, POINTER :: fc_ptr(:)
    INTEGER :: n
    TYPE(c_ptr)  :: fc_concen

    nsize = n

    CALL c_f_POINTER(fc_concen,fc_ptr,[1])
    concen(0:n-1,0:n-1)  => fc_ptr

END SUBROUTINE

SUBROUTINE assign_array()
    USE test_mod, ONLY: nsize, concen
    IMPLICIT NONE

    INTEGER :: i,j

    WRITE(6,*) 'Shape : ', SHAPE(concen)

    DO j = 0,nsize-1
        DO i = 0,nsize-1
            concen(i,j) = nsize*i + j
        end do
    ENDDO

END SUBROUTINE

SUBROUTINE print_array()
    USE test_mod,  ONLY  : nsize, concen
    IMPLICIT none

    INTEGER :: i,j

    DO i = 0,nsize-1
        DO j = 0,nsize-1
            WRITE(6,'(2F24.3)') concen(i,j)
        END DO
    END DO

END SUBROUTINE print_array

max_stars_repo_name: ECLAIRWaveS/ForestClaw

max_issues_repo_name: ECLAIRWaveS/ForestClaw

max_forks_repo_name: ECLAIRWaveS/ForestClaw


text:
! SELF_Quadrature.f90
!
! Copyright 2020 Fluid Numerics LLC
! Author : Joseph Schoonover (joe@fluidnumerics.com)
! Support : self@higherordermethods.org
!
! //////////////////////////////////////////////////////////////////////////////////////////////// !
!
! Contains routines from D.A. Kopriva, 2009, "Implementing Spectral Methods for Partial
! Differential Equations: Algorithms for Scientists and Engineers", Springer.
!
! Routines are defined for computing Legendre and Chebyshev Gauss and Gauss-Lobatto
! quadrature nodes and weights.

MODULE SELF_Quadrature

  USE ISO_FORTRAN_ENV
  USE SELF_Constants

  IMPLICIT NONE

  PUBLIC  :: ChebyshevQuadrature,LegendreQuadrature
  PRIVATE :: ChebyshevGauss,ChebyshevGaussLobatto, &
             LegendreGauss,LegendreGaussLobatto, &
             LegendreQandL

CONTAINS

! =============================================================================================== !
! LegendreQuadrature
!   Returns the specified Legendre quadrature nodes and integration weights.
!
!   Given a polynomial degree, and quadrature type (Gauss or Gauss Lobatto), this subroutine manages
!   the calls to underlying private routines to generate the desired Legendre quadrature.
!
!   Usage :
!
!     INTEGER    :: N, quadType
!     REAL(real64) :: nodes(0:N), weights(0:N)
!
!       CALL LegendreQuadrature( N, quadType, nodes, weights )
!
!   Parameters :
!
!     N (in)
!       Degree of the quadrature
!
!     quadType (in)
!       Flag specifying the quadrature type. Can be set to GAUSS or GAUSS_LOBATTO
!
!     nodes(0:N) (out)
!       Array of quadrature nodes
!
!     weights(0:N) (out)
!       Array of quadrature weights
!
! =============================================================================================== !

  SUBROUTINE LegendreQuadrature(N,nodes,weights,QuadType)
    IMPLICIT NONE
    INTEGER,INTENT(in)     :: N
    REAL(prec),INTENT(out) :: nodes(0:N)
    REAL(prec),INTENT(out) :: weights(0:N)
    INTEGER,INTENT(in)     :: QuadType
    ! Local
    REAL(real64) :: nodesLocal(0:N)
    REAL(real64) :: weightsLocal(0:N)
    INTEGER :: i

    IF (QuadType == GAUSS_LOBATTO) THEN

      CALL LegendreGaussLobatto(N,nodesLocal,weightsLocal)

    ELSEIF (QuadType == GAUSS) THEN

      CALL LegendreGauss(N,nodesLocal,weightsLocal)

    END IF

    DO i = 0,N
      nodes(i) = REAL(nodesLocal(i),prec)
      weights(i) = REAL(weightsLocal(i),prec)
    END DO

  END SUBROUTINE LegendreQuadrature

! =============================================================================================== !
! ChebyshevQuadrature
!
!   Returns the specified Chebyshev quadrature nodes and integration weights.
!
!   Given a polynomial degree, and quadrature type (Gauss or Gauss Lobatto), this subroutine manages
!   the calls to underlying private routines to generate the desired Chebyshev quadrature.
!
!   Usage :
!
!     INTEGER    :: N, quadType
!     REAL(real64) :: nodes(0:N), weights(0:N)
!
!       CALL ChebyshevQuadrature( N, quadType, nodes, weights )
!
!   Input/Output :
!
!     N (in)
!       Degree of the quadrature
!
!     quadType (in)
!       Flag specifying the quadrature type. Can be set to GAUSS or GAUSS_LOBATTO
!
!     nodes(0:N) (out)
!       Array of quadrature nodes
!
!     weights(0:N) (out)
!       Array of quadrature weights
!
! ================================================================================================ !

  SUBROUTINE ChebyshevQuadrature(N,nodes,weights,quadType)
    IMPLICIT NONE
    INTEGER,INTENT(in)     :: N
    REAL(prec),INTENT(out) :: nodes(0:N)
    REAL(prec),INTENT(out) :: weights(0:N)
    INTEGER,INTENT(in)     :: QuadType
    ! Local
    REAL(real64) :: nodesLocal(0:N)
    REAL(real64) :: weightsLocal(0:N)
    INTEGER :: i

    IF (QuadType == CHEBYSHEV_GAUSS_LOBATTO) then

      CALL ChebyshevGaussLobatto(N,nodesLocal,weightsLocal)

    ELSEIF (QuadType == CHEBYSHEV_GAUSS) then

      CALL ChebyshevGauss(N,nodesLocal,weightsLocal)

    END IF

    DO i = 0,N
      nodes(i) = REAL(nodesLocal(i),prec)
      weights(i) = REAL(weightsLocal(i),prec)
    END DO

  END SUBROUTINE ChebyshevQuadrature

! =============================================================================================== !
! S/R ChebyshevGauss
!   From Kopriva (2009) "Implementing Spectral Methods for Partial Differential Equations", p. 67
!   Algorithm 26
! =============================================================================================== !

  SUBROUTINE ChebyshevGauss(N,nodes,weights)
    IMPLICIT NONE
    INTEGER       :: N
    REAL(real64)    :: nodes(0:N)
    REAL(real64)    :: weights(0:N)
    ! Local
    INTEGER    :: j

    DO j = 0,N

      weights(j) = pi/(REAL(N,real64) + 1.0_real64)
      nodes(j) = -cos(pi*(2.0_real64*REAL(j,real64) + 1.0_real64)/(2.0_real64*REAL(N,real64)+2.0_real64))

    END DO

  END SUBROUTINE ChebyshevGauss

! =============================================================================================== !
! S/R ChebyshevGaussLobatto
!   From Kopriva (2009) "Implementing Spectral Methods for Partial Differential Equations", p. 68
!   Algorithm 27
! =============================================================================================== !

  SUBROUTINE ChebyshevGaussLobatto(N,nodes,weights)
    IMPLICIT NONE
    INTEGER       :: N
    REAL(real64)    :: nodes(0:N)
    REAL(real64)    :: weights(0:N)
    ! LOCAL
    INTEGER    :: j

    DO j = 0,N

      weights(j) = pi/REAL(N,real64)
      nodes(j) = -cos(pi*REAL(j,real64)/REAL(N,real64))

    END DO

    weights(0) = weights(0)*0.5_real64
    weights(N) = weights(N)*0.5_real64

  END SUBROUTINE ChebyshevGaussLobatto

! =============================================================================================== !
! S/R LegendreGauss
!   From Kopriva (2009) "Implementing Spectral Methods for Partial Differential Equations", p. 64
!   Algorithm 23
! =============================================================================================== !

  SUBROUTINE LegendreGauss(N,nodes,weights)
    IMPLICIT NONE
    INTEGER    :: N
    REAL(real64) :: nodes(0:N)
    REAL(real64) :: weights(0:N)
    ! Local
    REAL(real64) :: nodes_local(0:N)
    REAL(real64) :: weights_local(0:N)
    REAL(real64) :: lN1,dlN1
    REAL(real64) :: delta
    INTEGER  :: j,kIt

    IF (N == 0) then

      nodes_local(0) = 0.0_real64
      weights_local(0) = 2.0_real64

    ELSEIF (N == 1) then

      nodes_local(0) = -SQRT(1.0_real64/3.0_real64)
      weights_local(0) = 1.0_real64
      nodes_local(1) = -nodes(0)
      weights_local(1) = weights(0)

    ELSE

      DO j = 0, ((N + 1)/2)

        nodes_local(j) = -cos((2.0_real64*REAL(j,real64) + 1.0_real64)*pi/(2.0_real64*REAL(N,real64) + 1.0_real64))

        DO kIt = 1,newtonMax

          CALL LegendrePolynomial(N + 1,nodes_local(j),lN1,dlN1)
          delta = -lN1/dlN1
          nodes_local(j) = nodes_local(j) + delta
          IF (abs(delta) <= TOL*nodes_local(j)) EXIT

        END DO

        CALL LegendrePolynomial(N + 1,nodes_local(j),lN1,dlN1)
        weights_local(j) = 2.0_real64/((1.0_real64 - nodes_local(j)*nodes_local(j))*dlN1*dlN1)
        weights_local(N - j) = weights_local(j)
        nodes_local(N - j) = -nodes_local(j)

      END DO

    END IF

    IF (MOD(REAL(N,real64),2.0_real64) == 0.0_real64) then

      CALL LegendrePolynomial(N + 1,0.0_real64,lN1,dlN1)
      nodes_local(N/2) = 0.0_real64
      weights_local(N/2) = 2.0/(dlN1*dlN1)

    END IF

    DO j = 0,N
      nodes(j) = REAL(nodes_local(j),real64)
      weights(j) = REAL(weights_local(j),real64)
    END DO

  END SUBROUTINE LegendreGauss

  ! =============================================================================================== !
  ! S/R LegendreGaussLobatto
  !   From Kopriva (2009) "Implementing Spectral Methods for Partial Differential Equations", p. 66
  !   Algorithm 25
  ! =============================================================================================== !

  SUBROUTINE LegendreGaussLobatto(N,nodes,weights)
    IMPLICIT NONE
    INTEGER    :: N
    REAL(real64) :: nodes(0:N)
    REAL(real64) :: weights(0:N)
    ! Local
    REAL(real64) :: nodes_local(0:N)
    REAL(real64) :: weights_local(0:N)
    REAL(real64) :: delta,q,qprime,lN
    INTEGER  :: j,kIt

    IF (N == 1) then

      nodes_local(0) = -1.0_real64
      weights_local(0) = 1.0_real64
      nodes_local(1) = 1.0_real64
      weights_local(1) = 1.0_real64

    ELSE

      nodes_local(0) = -1.0_real64
      weights_local(0) = 2.0_real64/(REAL(N,real64)*(REAL(N,real64) + 1.0_real64))
      nodes_local(N) = 1.0_real64
      weights_local(N) = weights_local(0)

      DO j = 1, ((N + 1)/2 - 1)

        nodes_local(j) = -COS((REAL(j,real64) + 0.25_real64)*pi/REAL(N,real64) - &
                              3.0_real64/(8.0_real64*REAL(N,real64)*pi*(REAL(j,real64) + 0.25_real64)))

        DO kIt = 1,newtonMax

          CALL LegendreQandL(N,nodes_local(j),q,qprime,lN)

          delta = -q/qprime
          nodes_local(j) = nodes_local(j) + delta
          IF (ABS(delta) <= TOL*nodes_local(j)) EXIT

        END DO

        CALL LegendreQandL(N,nodes_local(j),q,qprime,lN)

        weights_local(j) = 2.0_real64/(REAL(N,real64)*(REAL(N,real64) + 1.0_real64)*lN*lN)
        weights_local(N - j) = weights_local(j)
        nodes_local(N - j) = -nodes_local(j)

      END DO

    END IF

    IF (MOD(REAL(N,real64),2.0_real64) == 0.0_real64) THEN

      CALL LegendreQandL(N,0.0_real64,q,qprime,lN)

      nodes_local(N/2) = 0.0_real64
      weights_local(N/2) = 2.0_real64/(REAL(N,real64)*(REAL(N,real64) + 1.0_real64)*lN*lN)

    END IF

    DO j = 0,N
      nodes(j) = REAL(nodes_local(j),real64)
      weights(j) = REAL(weights_local(j),real64)
    END DO

  END SUBROUTINE LegendreGaussLobatto

! =============================================================================================== !
! S/R LegendrePolynomial
!   From Kopriva (2009) "Implementing Spectral Methods for Partial Differential Equations", p. 63
!   Algorithm 22
! =============================================================================================== !

  SUBROUTINE LegendrePolynomial(N,x,lAtX,dLdxAtX)
    IMPLICIT NONE
    INTEGER     :: N
    REAL(real64)    :: x
    REAL(real64)    :: lAtX,dLdxAtX
    ! Local
    REAL(real64) :: lNm1,lNm2,dlNm1,dlNm2
    INTEGER  :: i

    IF (N == 0) then

      lAtX = 1.0_real64
      dLdxAtX = 0.0_real64

    ELSEIF (N == 1) then

      lAtX = x
      dLdxAtX = 1.0_real64

    ELSE

      lnM2 = 1.0_real64
      lnM1 = x
      dlnM2 = 0.0_real64
      dlnM1 = 1.0_real64

      DO i = 2,N

        lAtX = ((2.0_real64*REAL(i,real64) - 1.0_real64)*x*lnM1 - &
                (REAL(i,real64) - 1.0_real64)*lnM2)/(REAL(i,real64))

        dldxAtX = dlnM2 + (2.0_real64*REAL(i,real64) - 1.0_real64)*lnM1
        lnM2 = lnM1
        lnM1 = lAtX
        dlnM2 = dlnM1
        dlnM1 = dldxAtX

      END DO

    END IF

  END SUBROUTINE LegendrePolynomial

! =============================================================================================== !
! S/R LegendreQandL
!   From Kopriva (2009) "Implementing Spectral Methods for Partial Differential Equations", p. 65
!   Algorithm 24
! =============================================================================================== !

  SUBROUTINE LegendreQandL(N,x,q,qprime,lN)
    IMPLICIT NONE
    INTEGER    :: N
    REAL(real64) :: x
    REAL(real64) :: lN,q,qprime
    ! Local
    REAL(real64) :: lNm1,lNm2,dlNm1,dlNm2,dlN,lN1,dlN1
    INTEGER    :: i

    lNm2 = 1.0_real64
    lNm1 = x
    dlNm2 = 0.0_real64
    dlNm1 = 1.0_real64

    DO i = 2,N

      lN = (2.0_real64*i - 1.0_real64)/(REAL(i,real64))*x*lNm1 - (REAL(i,real64) - 1.0_real64)/(REAL(i,real64))*lNm2
      dlN = dlNm2 + (2.0_real64*REAL(i,real64) - 1.0_real64)*lNm1
      lNm2 = lNm1
      lNm1 = lN
      dlNm2 = dlNm1
      dlNm1 = dlN

    END DO

    i = N + 1
    lN1 = (2.0_real64*i - 1.0_real64)/(REAL(i,real64))*x*lN - (REAL(i,real64) - 1.0_real64)/(REAL(i,real64))*lNm2
    dlN1 = dlNm2 + (2.0_real64*REAL(i,real64) - 1.0_real64)*lNm1
    q = lN1 - lNm2
    qprime = dlN1 - dlNm2

  END SUBROUTINE LegendreQandL

END MODULE SELF_Quadrature

max_stars_repo_name: FluidNumerics/self-fluids

max_issues_repo_name: FluidNumerics/self-fluids

max_forks_repo_name: FluidNumerics/SELF-Fluids


text:
program farewell_test
    use farewell_m, only: make_farewell
    use iso_fortran_env, only: error_unit, output_unit

    implicit none

    character(len=:), allocatable :: farewell

    allocate(character(len=0) :: farewell)
    farewell = make_farewell("World")

    if (farewell == "Goodbye, World!") then
        write(output_unit, *) "Passed"
    else
        write(error_unit, *) "Failed"
        call exit(1)
    end if
end program farewell_test

max_stars_repo_name: LKedward/fpm-haskell

max_issues_repo_name: LKedward/fpm-haskell

max_forks_repo_name: LKedward/fpm-haskell


text:
	SUBROUTINE INTERPOLATE_SSP(NAME)

!	Interpolates BC03 model SSP at metallicity Z/Zo indicated in name (after character ',')

!	Variables
	parameter (jb=10)
	parameter (mmod=16)
	include 'csp.dec'
	integer bcopt
	real w(imw),h(imw)
        logical stelib
	character*(*) name,files(mmod)*128
	real zmod(mmod),zx(mmod)
!	data zmod/0.0001, 0.0004,  0.004, 0.008,  0.02,  0.05,  0.10/   						! BC03
!	data zmod/0.0001, 0.0002, 0.0005, 0.001, 0.002, 0.004, 0.006, 0.008, 0.010, 0.014, 0.017, 0.02, 0.03, 0.04/	! CB13

!	Check if interpolation has been requested
	i = index(name,',')
	call s500('c','',0.)
	if (i <= 0) then
!		write (501,*) 0.
!		write (501,*) 1.
!		write (501,*) name(:largo(name))
!		write (501,*) 0.
!		write (501,*) 'empty'
		call s500('f','',0.)
		call s500('f','',1.)
		call s500('a',name(:largo(name)),0.)
		call s500('f','',0.)
		call s500('a','empty',0.)
		return
	endif
	name = name(i+1:)
	ngu = nargu(name)
	if (ngu == 1) then
                read (name,*) z
		idth=0
	elseif (ngu == 2) then
                read (name,*) z,idth
	endif
	call s500('f','',z)

!	Check models in use. Define file names
	bcopt = 0
	call system('/bin/ls *ssp.ised > fort.202')
	read (202,'(a)') files(1)
	read (files(1),'(2x,i4)') bcopt
	close (202)
	call system('/bin/rm fort.202')
	if (bcopt == 2003) then
!		Use BC 2003 models
		nmod      = 6
		zsun      = 0.02
		zmod( 1)  = 0.0001
		zmod( 2)  = 0.0004
		zmod( 3)  = 0.004
		zmod( 4)  = 0.008
		zmod( 5)  = 0.02
		zmod( 6)  = 0.05
!		zmod( 7)  = 0.10
		if (index(files(1),'_stelib_') > 1) then
			files( 1) = 'bc2003_z0001_chab_hr_stelib_ssp.ised'	! = bc2003_hr_stelib_m22_chab_ssp.ised = bc2003_hr_m22_chab_ssp.ised
			files( 2) = 'bc2003_z0004_chab_hr_stelib_ssp.ised'	! = bc2003_hr_stelib_m32_chab_ssp.ised = bc2003_hr_m32_chab_ssp.ised
			files( 3) = 'bc2003_z004_chab_hr_stelib_ssp.ised'	! = bc2003_hr_stelib_m42_chab_ssp.ised = bc2003_hr_m42_chab_ssp.ised
			files( 4) = 'bc2003_z008_chab_hr_stelib_ssp.ised'	! = bc2003_hr_stelib_m52_chab_ssp.ised = bc2003_hr_m52_chab_ssp.ised
			files( 5) = 'bc2003_z020_chab_hr_stelib_ssp.ised'	! = bc2003_hr_stelib_m62_chab_ssp.ised = bc2003_hr_m62_chab_ssp.ised
			files( 6) = 'bc2003_z050_chab_hr_stelib_ssp.ised'	! = bc2003_hr_stelib_m72_chab_ssp.ised = bc2003_hr_m72_chab_ssp.ised
!			files( 7) = 'bc2003_z100_chab_hr_stelib_ssp.ised'	! = bc2003_hr_stelib_m82_chab_ssp.ised = bc2003_hr_m82_chab_ssp.ised
			bcopt     = 1
		elseif (index(files(1),'_xmiless_') > 1) then
			files( 1) = 'bc2003_z0001_chab_hr_xmiless_ssp.ised'	! = bc2003_hr_xmiless_m22_chab_ssp.ised
			files( 2) = 'bc2003_z0004_chab_hr_xmiless_ssp.ised'	! = bc2003_hr_xmiless_m32_chab_ssp.ised
			files( 3) = 'bc2003_z004_chab_hr_xmiless_ssp.ised'	! = bc2003_hr_xmiless_m42_chab_ssp.ised
			files( 4) = 'bc2003_z008_chab_hr_xmiless_ssp.ised'	! = bc2003_hr_xmiless_m52_chab_ssp.ised
			files( 5) = 'bc2003_z020_chab_hr_xmiless_ssp.ised'	! = bc2003_hr_xmiless_m62_chab_ssp.ised
			files( 6) = 'bc2003_z050_chab_hr_xmiless_ssp.ised'	! = bc2003_hr_xmiless_m72_chab_ssp.ised
!			files( 7) = 'bc2003_z100_chab_hr_xmiless_ssp.ised'	! = bc2003_hr_xmiless_m82_chab_ssp.ised
			bcopt     = 2
		else
			write (6,*) 'Unknown series of models: ',files(1)(:largo(files(1)))
			stop
		endif
	elseif (bcopt == 2013 .or. bcopt == 2016 .or. bcopt == 2017) then
!		Use CB 2013 models
		nmod      = 16
		zsun      = 0.017
		zmod( 1)  = 0.0000
		zmod( 2)  = 0.0001
		zmod( 3)  = 0.0002
		zmod( 4)  = 0.0005
		zmod( 5)  = 0.001
		zmod( 6)  = 0.002
		zmod( 7)  = 0.004
		zmod( 8)  = 0.006
		zmod( 9)  = 0.008
		zmod(10)  = 0.010
		zmod(11)  = 0.014
		zmod(12)  = 0.017
		zmod(13)  = 0.020
		zmod(14)  = 0.030
		zmod(15)  = 0.040
		zmod(16)  = 0.060
		if (index(files(1),'_stelib_') > 1) then
			files( 1) = 'cb2016_z0000_chab_hr_stelib_ssp.ised'	! = cb2013_s2_z0000y230n_chab_hr_stelib_ssp.ised
			files( 2) = 'cb2016_z0001_chab_hr_stelib_ssp.ised'	! = cb2013_s2_z0001y249n_chab_hr_stelib_ssp.ised
			files( 3) = 'cb2016_z0002_chab_hr_stelib_ssp.ised'	! = cb2013_s2_z0002y249n_chab_hr_stelib_ssp.ised
			files( 4) = 'cb2016_z0005_chab_hr_stelib_ssp.ised'	! = cb2013_s2_z0005y249n_chab_hr_stelib_ssp.ised
			files( 5) = 'cb2016_z001_chab_hr_stelib_ssp.ised' 	! = cb2013_s2_z001y250n_chab_hr_stelib_ssp.ised
			files( 6) = 'cb2016_z002_chab_hr_stelib_ssp.ised' 	! = cb2013_s2_z002y252n_chab_hr_stelib_ssp.ised
			files( 7) = 'cb2016_z004_chab_hr_stelib_ssp.ised' 	! = cb2013_s2_z004y256n_chab_hr_stelib_ssp.ised
			files( 8) = 'cb2016_z006_chab_hr_stelib_ssp.ised' 	! = cb2013_s2_z006y259n_chab_hr_stelib_ssp.ised
			files( 9) = 'cb2016_z008_chab_hr_stelib_ssp.ised' 	! = cb2013_s2_z008y263n_chab_hr_stelib_ssp.ised
			files(10) = 'cb2016_z010_chab_hr_stelib_ssp.ised' 	! = cb2013_s2_z010y267n_chab_hr_stelib_ssp.ised
			files(11) = 'cb2016_z014_chab_hr_stelib_ssp.ised' 	! = cb2013_s2_z014y273n_chab_hr_stelib_ssp.ised
			files(12) = 'cb2016_z017_chab_hr_stelib_ssp.ised' 	! = cb2013_s2_z017y279n_chab_hr_stelib_ssp.ised
			files(13) = 'cb2016_z020_chab_hr_stelib_ssp.ised' 	! = cb2013_s2_z020y284n_chab_hr_stelib_ssp.ised
			files(14) = 'cb2016_z030_chab_hr_stelib_ssp.ised' 	! = cb2013_s2_z030y302n_chab_hr_stelib_ssp.ised
			files(15) = 'cb2016_z040_chab_hr_stelib_ssp.ised' 	! = cb2013_s2_z040y321n_chab_hr_stelib_ssp.ised
			files(16) = 'cb2016_z060_chab_hr_stelib_ssp.ised' 	! = cb2013_s2_z060y356n_chab_hr_stelib_ssp.ised
			bcopt     = 3
		elseif (index(files(1),'_xmilesi_') > 1) then
			files( 1) = 'cb2017_z0000_chab_hr_xmilesi_ssp.ised'	! = cb2013_s2_z0000y230n_chab_hr_xmilesi_ssp.ised
			files( 2) = 'cb2017_z0001_chab_hr_xmilesi_ssp.ised'	! = cb2013_s2_z0001y249n_chab_hr_xmilesi_ssp.ised
			files( 3) = 'cb2017_z0002_chab_hr_xmilesi_ssp.ised'	! = cb2013_s2_z0002y249n_chab_hr_xmilesi_ssp.ised
			files( 4) = 'cb2017_z0005_chab_hr_xmilesi_ssp.ised'	! = cb2013_s2_z0005y249n_chab_hr_xmilesi_ssp.ised
			files( 5) = 'cb2017_z001_chab_hr_xmilesi_ssp.ised' 	! = cb2013_s2_z001y250n_chab_hr_xmilesi_ssp.ised
			files( 6) = 'cb2017_z002_chab_hr_xmilesi_ssp.ised' 	! = cb2013_s2_z002y252n_chab_hr_xmilesi_ssp.ised
			files( 7) = 'cb2017_z004_chab_hr_xmilesi_ssp.ised' 	! = cb2013_s2_z004y256n_chab_hr_xmilesi_ssp.ised
			files( 8) = 'cb2017_z006_chab_hr_xmilesi_ssp.ised' 	! = cb2013_s2_z006y259n_chab_hr_xmilesi_ssp.ised
			files( 9) = 'cb2017_z008_chab_hr_xmilesi_ssp.ised' 	! = cb2013_s2_z008y263n_chab_hr_xmilesi_ssp.ised
			files(10) = 'cb2017_z010_chab_hr_xmilesi_ssp.ised' 	! = cb2013_s2_z010y267n_chab_hr_xmilesi_ssp.ised
			files(11) = 'cb2017_z014_chab_hr_xmilesi_ssp.ised' 	! = cb2013_s2_z014y273n_chab_hr_xmilesi_ssp.ised
			files(12) = 'cb2017_z017_chab_hr_xmilesi_ssp.ised' 	! = cb2013_s2_z017y279n_chab_hr_xmilesi_ssp.ised
			files(13) = 'cb2017_z020_chab_hr_xmilesi_ssp.ised' 	! = cb2013_s2_z020y284n_chab_hr_xmilesi_ssp.ised
			files(14) = 'cb2017_z030_chab_hr_xmilesi_ssp.ised' 	! = cb2013_s2_z030y302n_chab_hr_xmilesi_ssp.ised
			files(15) = 'cb2017_z040_chab_hr_xmilesi_ssp.ised' 	! = cb2013_s2_z040y321n_chab_hr_xmilesi_ssp.ised
			files(16) = 'cb2017_z060_chab_hr_xmilesi_ssp.ised' 	! = cb2013_s2_z060y356n_chab_hr_xmilesi_ssp.ised
			bcopt     = 4
		else
			write (6,*) 'Unknown series of models: ',files(1)(:largo(files(1)))
			stop
		endif
	else
		write (6,*) 'Unknown model set:',bcopt
		stop
	endif

!	Express Z in solar units
	do i=1,nmod
	zx(i)=zmod(i)/zsun
	if (z == zx(i)) then
		name = files(i)
!		write (501,*) 1.
!		write (501,*) name(:largo(name))
!		write (501,*) 0.
!		write (501,*) 'empty'
		call s500('f','',1.)
		call s500('a',name(:largo(name)),0.)
		call s500('f','',0.)
		call s500('a','empty',0.)
		return
	endif
	enddo

!	Check extreme values
	if ( z <= zx(1) ) then
		name = files(1)
!		write (501,*) 1.
!		write (501,*) name(:largo(name))
!		write (501,*) 0.
!		write (501,*) 'empty'
		call s500('f','',1.)
		call s500('a',name(:largo(name)),0.)
		call s500('f','',0.)
		call s500('a','empty',0.)
		return
	elseif ( z >= zx(nmod) ) then
		name = files(nmod)
!		write (501,*) 1.
!		write (501,*) name(:largo(name))
!		write (501,*) 0.
!		write (501,*) 'empty'
		call s500('f','',1.)
		call s500('a',name(:largo(name)),0.)
		call s500('f','',0.)
		call s500('a','empty',0.)
		return
	endif

!	Find bracketing Z's
	do i=1,nmod-1
	if ( z > zx(i) .and. z <= zx(i+1) ) then
		write (6,*) 'Using z:',i,z,zx(i),zx(i+1)
		z1 = alog10(zx(i))
		z2 = alog10(zx(i+1))
		zo = alog10(z)
		w1 = (z2-zo)/(z2-z1)
		w2 = 1. - w1
		i1 = i
		i2 = i+1
!		write (501,*) w1
!		write (501,*) files(i1)(:largo(files(i1)))
!		write (501,*) w2
!		write (501,*) files(i2)(:largo(files(i2)))
		call s500('f','',w1)
		call s500('a',files(i1)(:largo(files(i1))),0.)
		call s500('f','',w2)
		call s500('a',files(i2)(:largo(files(i2))),0.)
	endif
	enddo
!	Build filename of temporary .ised file
	if (bcopt == 1) then
		name = 'bc2003_hr_stelib_mxx_chab_ssp.ised'
	elseif (bcopt == 2) then
		name = 'bc2003_hr_xmiless_mxx_chab_ssp.ised'
	elseif (bcopt == 3) then
!		name = 'cb2013_s2_mxx_chab_hr_stelib_ssp.ised'
		name = 'cb2016_mxx_chab_hr_stelib_ssp.ised'
	elseif (bcopt == 4) then
!		name = 'cb2013_s2_mxx_chab_hr_xmiless_ssp.ised'
!		name = 'cb2016_mxx_chab_hr_xmiless_ssp.ised'
		name = 'cb2017_mxx_chab_hr_xmilesi_ssp.ised'
	else
		write (6,*) 'Stopping on bcopt =',bcopt
	endif
	call gen_file_name(name,z,idth)
	open  (3,file=name,form='unformatted',status='unknown')

!	Perform interpolation
	write (6,*) 'Reading file: ',files(i1)
	open  (1,file=files(i1),form='unformatted',status='old')
!	Read basic parameters from SSP file
	read  (1) nsteps,(tb(i),i=0,nsteps-1),ml,mu,iseg,(xx(i),lm(i),um(i),baux(i),cn(i),cc(i),i=1,iseg),totm,totn,avs,jo,tauo,id,tcut,ttt,ttt,ttt,id,id,igw,stelib
!	Modify id record
	i=index(id,'X=')
	write (id(i:),'(a,f6.3)') 'SSP interpolated to Z/Zo = ',z
	write (3) nsteps,(tb(i),i=0,nsteps-1),ml,mu,iseg,(xx(i),lm(i),um(i),baux(i),cn(i),cc(i),i=1,iseg),totm,totn,avs,jo,tauo,id,tcut,ttt,ttt,ttt,id,id,igw,stelib
	read  (1) inl,(w(i),i=1,inl)
	write (3) inl,(w(i),i=1,inl)
	write (6,*) 'Reading file: ',files(i2)
	write (6,*) 'Weights =',w1,w2
	write (6,*) 'Writing file: ',name(:largo(name))
	open  (2,file=files(i2),form='unformatted',status='old')
!	Read basic parameters from SSP file
	read  (2)
	read  (2)
	do n=0,nsteps-1
	read  (1) inl,(w(i),i=1,inl)
	read  (2) inl,(h(i),i=1,inl)
	write (3) inl,(w1*w(i)+w2*h(i),i=1,inl)
	enddo
	read  (1,end=1) nsteps,(w(i),i=1,nsteps)
	read  (2,end=1) nsteps,(h(i),i=1,nsteps)
	write (3)       nsteps,(w1*w(i)+w2*h(i),i=1,nsteps)
	read  (1,end=1) nsteps,(w(i),i=1,nsteps)
	read  (2,end=1) nsteps,(h(i),i=1,nsteps)
	write (3)       nsteps,(w1*w(i)+w2*h(i),i=1,nsteps)
	read  (1,end=1) nsteps,(w(i),i=1,nsteps)
	read  (2,end=1) nsteps,(h(i),i=1,nsteps)
	write (3)       nsteps,(w1*w(i)+w2*h(i),i=1,nsteps)
	read  (1,end=1) nsteps,(w(i),i=1,nsteps)
	read  (2,end=1) nsteps,(h(i),i=1,nsteps)
	write (3)       nsteps,(w1*w(i)+w2*h(i),i=1,nsteps)
	read  (1,end=1) nsteps,(w(i),i=1,nsteps)
	read  (2,end=1) nsteps,(h(i),i=1,nsteps)
	write (3)       nsteps,(w1*w(i)+w2*h(i),i=1,nsteps)
	read  (1,end=1) nsteps,(w(i),i=1,nsteps)
	read  (2,end=1) nsteps,(h(i),i=1,nsteps)
	write (3)       nsteps,(w1*w(i)+w2*h(i),i=1,nsteps)
	read  (1,end=1) nsteps,(w(i),i=1,nsteps)
	read  (2,end=1) nsteps,(h(i),i=1,nsteps)
	write (3)       nsteps,(w1*w(i)+w2*h(i),i=1,nsteps)
	read  (1,end=1) nsteps,(w(i),i=1,nsteps)
	read  (2,end=1) nsteps,(h(i),i=1,nsteps)
	write (3)       nsteps,(w1*w(i)+w2*h(i),i=1,nsteps)
	read  (1,end=1) nsteps,(w(i),i=1,nsteps)
	read  (2,end=1) nsteps,(h(i),i=1,nsteps)
	write (3)       nsteps,(w1*w(i)+w2*h(i),i=1,nsteps)
	read  (1,end=1) nsteps,(w(i),i=1,nsteps)
	read  (2,end=1) nsteps,(h(i),i=1,nsteps)
	write (3)       nsteps,(w1*w(i)+w2*h(i),i=1,nsteps)
	read  (1,end=1) nsteps,(w(i),i=1,nsteps)
	read  (2,end=1) nsteps,(h(i),i=1,nsteps)
	write (3)       nsteps,(w1*w(i)+w2*h(i),i=1,nsteps)
	read  (1,end=1) nsteps,(w(i),i=1,nsteps)
	read  (2,end=1) nsteps,(h(i),i=1,nsteps)
	write (3)       nsteps,(w1*w(i)+w2*h(i),i=1,nsteps)
	read  (1,end=1) nsteps,(w(i),i=1,nsteps)
	read  (2,end=1) nsteps,(h(i),i=1,nsteps)
	write (3)       nsteps,(w1*w(i)+w2*h(i),i=1,nsteps)
	read  (1,end=1) nsteps,(w(i),i=1,nsteps)
	read  (2,end=1) nsteps,(h(i),i=1,nsteps)
	write (3)       nsteps,(w1*w(i)+w2*h(i),i=1,nsteps)
	read  (1,end=1) nsteps,(w(i),i=1,nsteps)
	read  (2,end=1) nsteps,(h(i),i=1,nsteps)
	write (3)       nsteps,(w1*w(i)+w2*h(i),i=1,nsteps)
	read  (1,end=1) nsteps,(w(i),i=1,nsteps)
	read  (2,end=1) nsteps,(h(i),i=1,nsteps)
	write (3)       nsteps,(w1*w(i)+w2*h(i),i=1,nsteps)
1	close (1)
	close (2)
	close (3)
	return
	end

	subroutine gen_file_name(name,z,idth)

!	Generates unique and new output file name

	character*(*) name,aux*6,bux*128

	i = index(name,'mxx')+3
	write (bux,'(a,i2,a)') '_t',idth,'_' // name(i+1:largo(name))
	if (bux(3:3) == ' ') bux(3:3) = '0'
	do j=0,10000
	write (aux,'(a,i5)') 'z',int(z*10**4)+j
	do k=1,6
	if (aux(k:k) == ' ') aux(k:k) ='0'
	enddo
	name = name(:i) // aux // bux
	open(99,file=name,status='new',iostat=iopenstatus)
	if (iopenstatus == 0) then
		close (99)
		return
	endif
	enddo
	return
	end

max_stars_repo_name: jchavesmontero/galaxpy

max_issues_repo_name: jchavesmontero/galaxpy

max_forks_repo_name: jchavesmontero/galaxpy


text:
-makelib xcelium_lib/xil_defaultlib \
  "../../../pipline_cpu/pipline_cpu.srcs/sources_1/new/pip_cpu.v" \
-endlib
-makelib xcelium_lib/xil_defaultlib \
  glbl.v
-endlib


text:
SUBROUTINE conini(actio)

  !*** contact data READ

  USE c_input
  USE ctrl_db, ONLY: ctype, mconts, ndime, npoin, nconp, numct, ndimc, cactive, memo
  USE npo_db,  ONLY: fcont
  USE outp_db, ONLY: nreqc, nprqc, iwrit
  IMPLICIT NONE

  CHARACTER(len=*),INTENT(IN):: actio

  INTEGER (kind=4) i, ncont

  INTERFACE
    INCLUDE 'coninp.h'
  END INTERFACE

  CALL listen('CONINI')
  IF (.NOT.exists('CONTAC')) THEN
    IF ( numct == 0 ) THEN
      IF( .NOT.ASSOCIATED(fcont)) ALLOCATE( fcont(1,1) )
      IF (nreqc > 0) THEN
        DEALLOCATE (nprqc)
        ALLOCATE (nprqc(1))
      END IF
    END IF
    backs = .TRUE.

  ELSE

    WRITE (lures,"(/,'  C O N T A C T   P A R A M E T E R S ',/)",ERR=9999)
    DO
      CALL listen('CONINI')
      IF (exists('ENDCON')) EXIT
      IF (exists('NCONT')) THEN
        ncont = getint('NCONT ',ndime,' Contact Type  ....................')
        IF(ncont < 0 .OR. ncont > mconts) CALL runend('INPDAT:NCONT MUST BE 1-4')
      ELSE IF (exists('DRAWBE')) THEN
        ncont = 4
        backs = .TRUE.
      ELSE
        ncont = ndime
        backs = .TRUE.
      END IF
      CALL coninp (actio, ncont, iwrit, maxve=memo)
      cactive(ncont) = .TRUE.   !once activated it will remain
    END DO
  END IF

  numct = 0     !initializes
  ctype = 0

  DO i=1,mconts                 !for each contact algorithm
    IF (nconp(i) > 0) THEN      !if there are contact pairs using this alg.
      numct = numct + 1         !increase the number of cont algors used
      ctype (numct) = i         !keep for this set the element type
    END IF
  END DO
  ndimc = ndime
  IF( nconp(1) > 0) ndimc = 4   !include axial rotation
  RETURN
 9999 CALL runen2('')
END SUBROUTINE conini

max_stars_repo_name: jerebenitez/IFE-simpact-openfoam

max_issues_repo_name: jerebenitez/IFE-simpact-openfoam

max_forks_repo_name: jerebenitez/IFE-simpact-openfoam


text:
C
C	$Id: intmsk.f,v 1.4 2008-07-27 00:23:04 haley Exp $
C
C                Copyright (C)  2000
C        University Corporation for Atmospheric Research
C                All Rights Reserved
C
C The use of this Software is governed by a License Agreement.
C
      SUBROUTINE INTMSK(MASK,MASK3)
C
C     THIS SUBROUTINE INITIALIZES THE MASKS USED BY THE MODEL PICTURE
C
      DIMENSION MASK(*),MASK3(*)
      SAVE
C
C     GET NUMBER OF BITS PER WORD
C
      NBPW = I1MACH(5)
C
C     FILL MASKS
C
      IP1 = 1
      IP2 = 7
      MASK(1)  = 1
      MASK3(1) = 3
C
      DO  5 I=2,64
      MASK(I)  = 0
      MASK3(I) = 0
    5 CONTINUE
C
      DO 10 I=2,NBPW-1
      MASK(I)=ISHIFT(IP1,I-1)
      MASK3(I)=ISHIFT(IP2,I-2)
   10 CONTINUE
      MASK(NBPW)=ISHIFT(IP1,NBPW-1)
      MASK3(NBPW)=ISHIFT(MASK3(1),NBPW-2)
      RETURN
      END

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


text:
module clfind_commons
  use amr_commons, ONLY: qdp,dp

  integer::nparts,nparts_tot,npeaks,npeaks_tot
  real(dp)::tot_mass
  real(dp)::relevance_threshold=1.5
  real(dp)::density_threshold=0.
  real(dp)::mass_threshold=0.
  logical::merge_unbound=.false.,clinfo=.false.

  ! Big array for saddlepoint values
  real(dp),allocatable,dimension(:,:)::saddle_dens,saddle_dens_tot

  ! Peak patch properties
  real(dp),allocatable,dimension(:,:)::clump_size_tot,center_of_mass_tot,clump_momentum_tot,clump_force_tot
  real(dp),allocatable,dimension(:,:,:)::second_moments,second_moments_tot,Icl_d_3by3_tot,Icl_3by3_tot
  real(dp),allocatable,dimension(:)::min_dens_tot,av_dens_tot,phi_min_tot
  real(dp),allocatable,dimension(:)::max_dens_tot,e_kin_int_tot,e_bind_tot,e_thermal_tot
!  real(dp),allocatable,dimension(:)::e_kin_int_tot4,e_bind_tot4,e_thermal_tot4
  real(dp),allocatable,dimension(:)::clump_mass_tot,clump_vol_tot,clump_mass_tot4
  real(dp),allocatable,dimension(:,:)::peak_pos_tot,bulk_momentum_tot
  real(dp),allocatable,dimension(:)::saddle_max_tot
  real(dp),allocatable,dimension(:)::relevance_tot
  real(dp),allocatable,dimension(:)::phi_ref, phi_ref_tot
  real(dp),allocatable,dimension(:)::Psurf,Psurf_tot,v_therm_tot,v_rms_tot,m4_tot
  real(dp),allocatable,dimension(:)::e_bind_iso_tot,e_therm_iso_tot,e_kin_iso_tot,grav_term_tot,clump_virial_tot
  real(dp),allocatable,dimension(:)::peak_check,ball4_check,isodens_check,clump_check
  real(dp),allocatable,dimension(:)::Icl_tot,Icl_d_tot,Icl_dd_tot
  logical,allocatable,dimension(:)::contracting

  ! Test particles properties
  real(dp),allocatable,dimension(:)::denp ! Density of the cell containing a test particle. Davide: used by the clump finder.
  integer,allocatable,dimension(:)::iglobalp,icellp,levp,testp_sort ! Used to sort test particles by density
  integer,allocatable,dimension(:)::n_cells_tot,minmatch_tot,new_peak
  integer,allocatable,dimension(:)::occupied,occupied_all ! Tells whether there is already a sink in a clump.
  integer,allocatable,dimension(:)::form,form_all ! Tells whether a sink has to be formed within a clump.

end module clfind_commons

max_stars_repo_name: rknop/amuse

max_issues_repo_name: rknop/amuse

max_forks_repo_name: rknop/amuse


text:
C    @(#)kmprt2.f	20.3 2/13/96
      function kmprt2(i,j)

      include 'ipfinc/cgrate.inc'

      kmprt2 = lrat(1,i) - lrat(1,j)

      if (kmprt2 .eq. 0) kmprt2=lrat(2,i)-lrat(2,j)
      if (kmprt2 .eq. 0) kmprt2=lrat(3,i)-lrat(3,j)

      return
      end

max_stars_repo_name: mbheinen/bpa-ipf-tsp

max_issues_repo_name: cuihantao/bpa-ipf-tsp

max_forks_repo_name: mbheinen/bpa-ipf-tsp


text:
      FUNCTION GOLDEN(AX,BX,CX,F,TOL,XMIN)
      PARAMETER (R=.61803399,C=.38196602)
      X0=AX
      X3=CX
      IF(ABS(CX-BX).GT.ABS(BX-AX))THEN
        X1=BX
        X2=BX+C*(CX-BX)
      ELSE
        X2=BX
        X1=BX-C*(BX-AX)
      ENDIF
      F1=F(X1)
      F2=F(X2)
1     IF(ABS(X3-X0).GT.TOL*(ABS(X1)+ABS(X2)))THEN
        IF(F2.LT.F1)THEN
          X0=X1
          X1=X2
          X2=R*X1+C*X3
          F0=F1
          F1=F2
          F2=F(X2)
        ELSE
          X3=X2
          X2=X1
          X1=R*X2+C*X0
          F3=F2
          F2=F1
          F1=F(X1)
        ENDIF
      GOTO 1
      ENDIF
      IF(F1.LT.F2)THEN
        GOLDEN=F1
        XMIN=X1
      ELSE
        GOLDEN=F2
        XMIN=X2
      ENDIF
      RETURN
      END

ext: for


text:
!>  Module file pv_coupling.mod
!
!>  An interface to pressure-velocity coupling methods (SIMPLE, etc)

module pv_coupling

    use kinds, only : ccs_int
    use types, only: field, ccs_mesh
    use parallel_types, only: parallel_environment

    implicit none

    private

    public :: solve_nonlinear

    interface

    module subroutine solve_nonlinear(par_env, mesh, cps, it_start, it_end, u, v, p, p_prime, mf)
        class(parallel_environment), allocatable, intent(in) :: par_env
        type(ccs_mesh), intent(in) :: mesh
        integer(ccs_int), intent(in) :: cps, it_start, it_end
        class(field), intent(inout) :: u, v, p, p_prime, mf
    end subroutine solve_nonlinear

    end interface

end module pv_coupling

max_stars_repo_name: asimovpp/asimov-ccs

max_issues_repo_name: asimovpp/asimov-ccs

max_forks_repo_name: asimovpp/asimov-ccs


text:
      subroutine rdcmt(iunt, cmt)
      integer iunt, i1
c      character(300) line
      character(4) cmt
      character tmpcmt(4), ch
      logical cmtlin

      cmtlin = .true.
      do i1 = 1, 4
         tmpcmt(i1) = cmt(i1:i1)
      end do
 5    continue
      read(iunt,*,end=10) ch
      do i1 = 1, 4
         if(ch.eq.tmpcmt(i1)) goto 5
      end do

 10   continue
      backspace(iunt)

      return
      end


text:
c
c
c     ###################################################
c     ##  COPYRIGHT (C)  1990  by  Jay William Ponder  ##
c     ##              All Rights Reserved              ##
c     ###################################################
c
c     ############################################################
c     ##                                                        ##
c     ##  subroutine prtseq  --  output of biopolymer sequence  ##
c     ##                                                        ##
c     ############################################################
c
c
c     "prtseq" writes out a biopolymer sequence to an external
c     disk file with 15 residues per line and distinct chains
c     separated by blank lines
c
c
      subroutine prtseq (iseq)
      implicit none
      include 'sizes.i'
      include 'files.i'
      include 'sequen.i'
      integer i,k,iseq
      integer smax,smin
      integer size,start,stop
      logical opened
      character*1 letter
      character*23 fstr
      character*120 seqfile
c
c
c     open the output unit if not already done
c
      inquire (unit=iseq,opened=opened)
      if (.not. opened) then
         seqfile = filename(1:leng)//'.seq'
         call version (seqfile,'new')
         open (unit=iseq,file=seqfile,status='new')
      end if
c
c     write out a three-letter code sequence file
c
      do i = 1, nchain
         letter = chnnam(i)
         start = ichain(1,i)
         stop = ichain(2,i)
         size = stop - start + 1
         smax = 0
         do while (smax .lt. size)
            smin = smax + 1
            smax = smax + 15
            smax = min(smax,size)
            if (i.ne.1 .and. smin.eq.1)  write (iseq,'()')
            fstr = '(3x,a1,i6,1x,15(1x,a3))'
            write (iseq,fstr)  letter,smin,(seq(k+start-1),k=smin,smax)
         end do
      end do
c
c     close the output unit if opened by this routine
c
      if (.not. opened)  close (unit=iseq)
      return
      end

max_stars_repo_name: htna/HCsbLib

max_issues_repo_name: htna/HCsbLib

max_forks_repo_name: htna/HCsbLib


text:
!nesthd1
      program nesthd1
!----- GPL ---------------------------------------------------------------------
!
!  Copyright (C)  Stichting Deltares, 2011-2016.
!
!  This program is free software: you can redistribute it and/or modify
!  it under the terms of the GNU General Public License as published by
!  the Free Software Foundation version 3.
!
!  This program is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!  GNU General Public License for more details.
!
!  You should have received a copy of the GNU General Public License
!  along with this program.  If not, see <http://www.gnu.org/licenses/>.
!
!  contact: delft3d.support@deltares.nl
!  Stichting Deltares
!  P.O. Box 177
!  2600 MH Delft, The Netherlands
!
!  All indications and logos of, and references to, "Delft3D" and "Deltares"
!  are registered trademarks of Stichting Deltares, and remain the property of
!  Stichting Deltares. All rights reserved.
!
!-------------------------------------------------------------------------------
!  $Id: nesthd1.f90 5717 2016-01-12 11:35:24Z mourits $
!  $HeadURL: https://svn.oss.deltares.nl/repos/delft3d/tags/6686/src/tools_gpl/nesthd1/packages/nesthd1/src/nesthd1.f90 $
!****&******************************************************************
! Deltares                         marine and coastal management
!
! program            : nesthd1
! version            : v1.0
! date               : June 1997
! programmer         : Lamber Hulsen/Theo van der Kaaij
! version            : 1.50.01
!
! function           :    determine nest stations, weight factors and
!                         angles in an overall model to generate
!                         boundary conditions (with NESTHD2) at boundary
!                         support points in a nested model
! notes              : 1. the nest characteristics are determined both
!                         for a nesting at water level points as well as
!                         a nesting at velocity points. In this
!                         way you are free to choose the kind of
!                         boundary forcing when running NESTHD2 without
!                         re-running NESTHD1.
!                      2. NESTHD1 requires that the overall and nested
!                         models are defined in the same co-ordinate
!                         system.
!                      3. This version of NESTHD1 will NOT compose a
!                         grid for spherical or rectilinear models.
!                         For these models a user-defined grid is
!                         required.
!****&******************************************************************

      use precision
      use nesthd1_version_module
!
      integer      , dimension(:)    , allocatable :: ipx
      integer      , dimension(:)    , allocatable :: ipy
      integer      , dimension(:)    , allocatable :: itotpx
      integer      , dimension(:)    , allocatable :: itotpy
!
      integer      , dimension(:,:)  , allocatable :: icom1
      integer      , dimension(:,:)  , allocatable :: icom2
      integer      , dimension(:,:)  , allocatable :: mcbsp
      integer      , dimension(:,:)  , allocatable :: ncbsp
!
      integer      , dimension(:,:,:), allocatable :: mcnes
      integer      , dimension(:,:,:), allocatable :: ncnes
!
      real         , dimension(:)    , allocatable :: angle
!
      real         , dimension(:,:)  , allocatable :: x1
      real         , dimension(:,:)  , allocatable :: x2
      real         , dimension(:,:)  , allocatable :: y1
      real         , dimension(:,:)  , allocatable :: y2
      real         , dimension(:,:)  , allocatable :: xbnd
      real         , dimension(:,:)  , allocatable :: ybnd
!
      real         , dimension(:,:,:), allocatable :: weight
!
      character(1) , dimension(:)    , allocatable :: typbnd
      character(20), dimension(:)    , allocatable :: nambnd
!
      integer       lun   (  7   )
      character*80  filnam (7)
      character*80  CIDENT
      logical       spher_crs, spher_det
!
!     Version string
!
    cident = ' '
    call getfullversionstring_nesthd1(cident)
!-----------------------------------------------------------------------
!---- 0. Initialisation
!-----------------------------------------------------------------------
      lunscr = 6

!-----------------------------------------------------------------------
!---- 1. Open all files
!-----------------------------------------------------------------------
      write (*     ,'(/, 2a)') ' ',trim(cident)
      call opnfl1(lun   , filnam, CIDENT)

!-----------------------------------------------------------------------
!---- 2. Dynamic memory allocation (FMM)
!-----------------------------------------------------------------------
!
!     Get dimensions for allocation
      call getdim(lun   , mmax1 , nmax1  , &
     &            mmax2 , nmax2 , maxnrp , maxbnd )
!
      allocate(ipx   (maxnrp), STAT = istat)
      allocate(ipy   (maxnrp), STAT = istat)
      allocate(itotpx(maxnrp), STAT = istat)
      allocate(itotpy(maxnrp), STAT = istat)
!
      allocate(icom1 ( mmax1 , nmax1 ), STAT = istat)
      allocate(icom2 ( mmax2 , nmax2 ), STAT = istat)
      allocate(mcbsp ( maxbnd, 2     ), STAT = istat)
      allocate(ncbsp ( maxbnd, 2     ), STAT = istat)
!
      allocate(mcnes ( maxbnd, 2, 4     ), STAT = istat)
      allocate(ncnes ( maxbnd, 2, 4     ), STAT = istat)
!
      allocate(angle ( maxbnd), STAT = istat)
!
      allocate(x1    ( mmax1 , nmax1 ), STAT = istat)
      allocate(x2    ( mmax2 , nmax2 ), STAT = istat)
      allocate(y1    ( mmax1 , nmax1 ), STAT = istat)
      allocate(y2    ( mmax2 , nmax2 ), STAT = istat)
      allocate(xbnd  ( maxbnd, 2     ), STAT = istat)
      allocate(ybnd  ( maxbnd, 2     ), STAT = istat)
!
      allocate(weight( maxbnd, 2, 4     ), STAT = istat)
!
      allocate(typbnd (maxbnd), STAT = istat)
      allocate(nambnd (maxbnd), STAT = istat)

!-----------------------------------------------------------------------
!---- 3. Zeroise arrays
!-----------------------------------------------------------------------

      mcbsp = 0
      ncbsp = 0
      icom1 = 0
      icom2 = 0

      mcnes = 0
      ncnes = 0

      angle = 0.0

      x1    = 0.0
      x2    = 0.0
      y1    = 0.0
      y2    = 0.0

      weight= 0.0

!-----------------------------------------------------------------------
!---- 4. Read grid and boundary data models
!-----------------------------------------------------------------------

      call reargf(lun(1)    ,x1       ,y1       ,mmax1  ,nmax1  ,mc1   ,nc1, &
                & spher_crs )

      ipx = 0
      ipy = 0
      itotpx = 0
      itotpy = 0
      call inigrd(lun(2),ifout ,nmax1 ,mmax1 ,icom1, ipx,ipy, &
     &            itotpx,itotpy,idum                          )

      if (ifout .ne. 0) stop 'Error in grid enclosure overall model'

      call reargf(lun(3)   ,x2  ,y2 ,mmax2  ,nmax2  ,mc2   ,nc2   , &
               &  spher_det)
      if (spher_crs .neqv. spher_det) stop 'Coarse and detailled model should be in the same coordinate system'

      ipx = 0
      ipy = 0
      itotpx = 0
      itotpy = 0
      call inigrd(lun(4),ifout ,nmax2 ,mmax2 ,icom2, ipx,ipy, &
     &            itotpx,itotpy,idum                          )
      if (ifout .ne. 0) stop 'Error in grid enclosure detailed model'

      call dimbnd(lun(5), nobnd )
      if (nobnd .gt. maxbnd) stop 'Enlarge maxbnd'
      call reabnd(lun(5),mcbsp, ncbsp, typbnd, nambnd, &
     &            nobnd)

!-----------------------------------------------------------------------
!---- 5. Determine weight factors water level boundary
!-----------------------------------------------------------------------

      call detxy (mcbsp , ncbsp , x2     , y2    , icom2 , &
     &            xbnd  , ybnd  , mmax2  , nmax2 , maxbnd, &
     &            nobnd , 'WL'                          )

      call detnst(x1      , y1      , icom1   , xbnd    , ybnd     , &
     &            mcbsp   , ncbsp   , mcnes   , ncnes   , weight   , &
     &            mmax1   , nmax1   , mc1     , nc1     , maxbnd   , &
     &            nobnd   , spher_crs                   )

!-----------------------------------------------------------------------
!---- 6. Write weight factors water level boundary to administration
!        file and write stations to trisula input file
!-----------------------------------------------------------------------

      call wrinst(lun(6)   , mcbsp    , ncbsp    , mcnes    , ncnes   , &
     &            weight   , angle    , maxbnd   , nobnd    , 'WL'    )

      call wrista(lun(7), mcnes    , ncnes    , maxbnd)

!-----------------------------------------------------------------------
!---- 7. Determine weight factors velocity boundary but first
!        zeroise arrays
!-----------------------------------------------------------------------

      mcnes = 0
      ncnes = 0

      weight = 0.0

      call detxy  (mcbsp    ,ncbsp    ,x2    ,y2    ,icom2    , &
     &             xbnd     ,ybnd     ,mmax2 ,nmax2 ,maxbnd   , &
     &             nobnd    , 'UV'                          )

      call detnst (x1      ,y1      ,icom1    ,xbnd    ,ybnd     , &
     &             mcbsp   ,ncbsp   ,mcnes    ,ncnes   ,weight   , &
     &             mmax1   , nmax1  ,mc1      ,nc1     ,maxbnd   , &
     &             nobnd   , spher_crs                           )

!-----------------------------------------------------------------------
!---- 8. Determine angles velocity boundary
!-----------------------------------------------------------------------

      call detang (xbnd     ,ybnd   ,angle    ,mcbsp    ,ncbsp    , &
     &             icom2    ,mmax2  ,nmax2    ,maxbnd   ,nobnd )

!-----------------------------------------------------------------------
!---- 9. Write weight factors/angles velocity boundary
!-----------------------------------------------------------------------

      call wrinst(lun(6)   ,mcbsp    ,ncbsp    ,mcnes    ,ncnes    , &
     &            weight   ,angle    ,maxbnd   ,nobnd    ,'UV'        )

      call wrista(lun(7), mcnes    , ncnes    , maxbnd)

!-----------------------------------------------------------------------
!---- end program
!-----------------------------------------------------------------------

      call clsfil(lun   ,7     )

 9000 continue
      endprogram nesthd1

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


text:
MODULE class_terrain
    IMPLICIT NONE
    PRIVATE

    TYPE , PUBLIC :: terrain
        CHARACTER(LEN=80) :: filename
        INTEGER , POINTER :: map(:,:) => null()
        INTEGER :: n_cols = 0
        INTEGER :: n_rows = 0
        INTEGER :: cell_size = 0
        INTEGER :: no_value = 0

        CONTAINS
            PROCEDURE :: load => load_data
            PROCEDURE :: dump => dump_data
    END TYPE terrain

    CONTAINS
        SUBROUTINE load_data(this)
            CLASS(terrain) :: this
            INTEGER , DIMENSION(4) :: vars
            INTEGER , PARAMETER :: u_number = 10
            CHARACTER(LEN=512),DIMENSION(4) :: buffer
            INTEGER :: b_start, b_end
            INTEGER :: i, j, k, l, res
            INTEGER , EXTERNAL :: a2i

            OPEN(UNIT=u_number,FILE=this%filename,FORM=’FORMATTED’,IOSTAT=res)
            buffer = ’ ’
            DO i = 1, 4
                READ(UNIT=u_number,FMT=’(A)’,IOSTAT=res)
                buffer(i)
                IF(res /= 0) THEN
                    PRINT*, ’Error in reading line 1, status ’, res
                    CLOSE(UNIT=u_number)
                    RETURN
                END IF
            END DO

            CALL extract_header_lines(buffer,vars)
            this%n_cols = vars(1)
            this%n_rows = vars(2)
            this%cell_size = vars(3)
            this%no_value = vars(4)
            ALLOCATE(this%map(this%n_rows,this%n_cols),STAT=res)
            IF(res /= 0) THEN
                PRINT *, ’Allocationfailure, status ’, res
                CLOSE(UNIT=u_number)
            END IF
        END SUBROUTINE load_data

        SUBROUTINE dump_data(this)
            CLASS(terrain) :: this
            INTEGER , PARAMETER :: u_number = 11
            INTEGER :: i, j, res

            OPEN(UNIT=u_number,FILE=this%filename,FORM=’FORMATTED’,IOSTAT=res)
            WRITE(UNIT=u_number,FMT=’(I5)’,IOSTAT=res) this%n_cols
            WRITE(UNIT=u_number,FMT=’(I5)’,IOSTAT=res) this%n_rows
            WRITE(UNIT=u_number,FMT=’(I5)’,IOSTAT=res) this%cell_size
            WRITE(UNIT=u_number,FMT=’(I5)’,IOSTAT=res) this%no_value
            DO i = 1, this%n_rows
                DO j = 1, this%n_cols
                    WRITE(UNIT=u_number,FMT=’(I4,A2)’,ADVANCE=’NO’,IOSTAT=res) &
                    this%map(i,j), ’, ’
                END DO
                WRITE(UNIT=u_number,FMT=’(A)’,IOSTAT=res) ’ ’
            END DO
            CLOSE(UNIT=u_number)
        END SUBROUTINE dump_data

END MODULE class_terrain

PROGRAM ctest
    USE class_terrain
    IMPLICITNONE
    TYPE(terrain) :: td1

    td1 = terrain(’demfile.dat’)
    CALL td1%fload()
END PROGRAM ctest

max_stars_repo_name: jfitz/code-stat

max_issues_repo_name: jfitz/code-stat

max_forks_repo_name: jfitz/code-stat


text:
! Copyright (c) 2012, NVIDIA CORPORATION.  All rights reserved.
!
! Licensed under the Apache License, Version 2.0 (the "License");
! you may not use this file except in compliance with the License.
! You may obtain a copy of the License at
!
!     http://www.apache.org/licenses/LICENSE-2.0
!
! Unless required by applicable law or agreed to in writing, software
! distributed under the License is distributed on an "AS IS" BASIS,
! WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
! See the License for the specific language governing permissions and
! limitations under the License.
!

module mod
logical rslt(4),expect(4)
type :: stuff
integer :: kk1=selected_int_kind(2)
integer :: ll1 = 10
character(len=10) :: cc
integer dd
integer:: pp(10)
end type

type :: objects(l1,k1)
integer, kind:: k1=selected_int_kind(2)
integer, len :: l1 = 10
character(len=l1) :: c
integer :: p(l1)
integer d
end type

end module

program ppp
use mod
integer y
type(objects)::z
!type(objects(l1=10)):: q

rslt(1) = z%l1 .eq. 10
do i=1,z%l1
  z%p(i)=i
enddo

rslt(2) = .true.
do i=1,z%l1
  if (z%p(i) .ne. i) rslt(2) = .false.
enddo

z%c = 'abcdefghij'
rslt(3) = z%c .eq. 'abcdefghij'

rslt(4) = size(z%p) .eq. z%l1

expect = .true.
call check(rslt,expect,4)

end

max_stars_repo_name: kammerdienerb/flang

max_issues_repo_name: kammerdienerb/flang

max_forks_repo_name: kammerdienerb/flang


text:
! MIT License
! Copyright (c) 2019 Sam Miller
! Permission is hereby granted, free of charge, to any person obtaining a copy
! of this software and associated documentation files (the "Software"), to deal
! in the Software without restriction, including without limitation the rights
! to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the Software is
! furnished to do so, subject to the following conditions:
!
! The above copyright notice and this permission notice shall be included in all
! copies or substantial portions of the Software.
!
! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
! SOFTWARE.

module mod_outlet_bc
  use, intrinsic :: iso_fortran_env, only: ik => int32, rk => real64
  use mod_error, only: error_msg
  use mod_globals, only: enable_debug_print, debug_print
  use mod_field, only: field_2d_t
  use mod_grid_block_2d, only: grid_block_2d_t
  use mod_boundary_conditions, only: boundary_condition_t
  use mod_input, only: input_t
  use mod_eos, only: eos
  use mod_nondimensionalization

  implicit none

  private
  public :: outlet_bc_t, outlet_bc_constructor

  type, extends(boundary_condition_t) :: outlet_bc_t
    real(rk) :: outlet_pressure = 0.0_rk
  contains
    procedure, public :: apply => apply_outlet_primitive_var_bc
    final :: finalize
  endtype
contains

  function outlet_bc_constructor(location, input, grid) result(bc)
    type(outlet_bc_t), pointer :: bc
    character(len=2), intent(in) :: location !< Location (+x, -x, +y, or -y)
    class(input_t), intent(in) :: input
    class(grid_block_2d_t), intent(in) :: grid

    allocate(bc)
    bc%name = 'outlet'
    bc%location = location
    bc%priority = 0
    bc%outlet_pressure = input%constant_bc_pressure_value * press_to_nondim

    if(bc%outlet_pressure <= 0.0_rk) then
      call error_msg(module_name='mod_outlet_bc', class_name='outlet_bc_t', &
                     procedure_name='outlet_bc_constructor', &
                     message="Outlet ambient pressure is <= 0!", &
                     file_name=__FILE__, line_number=__LINE__)
    endif

    call bc%set_indices(grid)
  endfunction outlet_bc_constructor

  subroutine apply_outlet_primitive_var_bc(self, rho, u, v, p)

    class(outlet_bc_t), intent(inout) :: self
    class(field_2d_t), intent(inout) :: rho  !< density
    class(field_2d_t), intent(inout) :: u    !< x-velocity
    class(field_2d_t), intent(inout) :: v    !< y-velocity
    class(field_2d_t), intent(inout) :: p    !< pressure

    real(rk) :: cs, gamma, vel, nx, ny, mach

    integer(ik) :: i, j, ihi, ihi_halo

    gamma = eos%get_gamma()

    select case(self%location)
    case('+x')

      if(rho%on_ihi_bc) then
        if(enable_debug_print) then
          call debug_print('Running outlet_bc_t%apply_outlet_primitive_var_bc() +x', &
                           __FILE__, __LINE__)
        endif

        ihi = rho%ihi
        ihi_halo = rho%ihi_halo

        do i = ihi + 1, ihi_halo
          do j = rho%jlo_halo, rho%jhi_halo

            associate(p_d => p%data(ihi, j), p_b => self%outlet_pressure, rho_0 => rho%data(ihi, j))

              ! Sometimes at startup there are zero rho/p for some reason, but it doesnt break the code
              if(p_b > tiny(1.0_rk) .and. rho_0 > tiny(1.0_rk)) then
                cs = sqrt(gamma * abs(p_d / rho_0))
                vel = sqrt(u%data(ihi, j)**2 + v%data(ihi, j)**2)
                mach = vel / cs

                nx = u%data(ihi, j) / (vel + 1e-30_rk)
                ny = v%data(ihi, j) / (vel + 1e-30_rk)

                if(mach > 1.0_rk) then
                  ! Supersonic extrapolation
                  rho%data(i, j) = rho%data(ihi, j)
                  u%data(i, j) = u%data(ihi, j)
                  v%data(i, j) = v%data(ihi, j)
                  p%data(i, j) = p%data(ihi, j)

                else
                  ! Subsonic extrapolation
                  rho%data(i, j) = rho%data(ihi, j) + (p_b - p_d) / cs**2
                  u%data(i, j) = u%data(ihi, j) + nx * (p_d - p_b) / (rho_0 * cs)
                  v%data(i, j) = v%data(ihi, j) + ny * (p_d - p_b) / (rho_0 * cs)
                  p%data(i, j) = self%outlet_pressure
                endif
              endif
            endassociate
          enddo

        enddo
      endif
    case('-x')
      if(rho%on_ilo_bc) then
        if(enable_debug_print) then
          call debug_print('Running outlet_bc_t%apply_outlet_primitive_var_bc() -x', &
                           __FILE__, __LINE__)
        endif
        do i = 1, self%n_ghost_layers
          rho%data(self%ilo_ghost(i), :) = rho%data(self%ilo, :)
          u%data(self%ilo_ghost(i), :) = u%data(self%ilo, :)
          v%data(self%ilo_ghost(i), :) = v%data(self%ilo, :)
          p%data(self%ilo_ghost(i), :) = p%data(self%ilo, :)
        enddo
      endif
    case('+y')
      if(rho%on_jhi_bc) then
        if(enable_debug_print) then
          call debug_print('Running outlet_bc_t%apply_outlet_primitive_var_bc() +y', &
                           __FILE__, __LINE__)
        endif
        do i = 1, self%n_ghost_layers
          rho%data(:, self%jhi_ghost(i)) = rho%data(:, self%jhi)
          u%data(:, self%jhi_ghost(i)) = u%data(:, self%jhi)
          v%data(:, self%jhi_ghost(i)) = v%data(:, self%jhi)
          p%data(:, self%jhi_ghost(i)) = p%data(:, self%jhi)
        enddo
      endif
    case('-y')
      if(rho%on_jlo_bc) then
        if(enable_debug_print) then
          call debug_print('Running outlet_bc_t%apply_outlet_primitive_var_bc() -y', &
                           __FILE__, __LINE__)
        endif
        do i = 1, self%n_ghost_layers
          rho%data(:, self%jlo_ghost(i)) = rho%data(:, self%jlo)
          u%data(:, self%jlo_ghost(i)) = u%data(:, self%jlo)
          v%data(:, self%jlo_ghost(i)) = v%data(:, self%jlo)
          p%data(:, self%jlo_ghost(i)) = p%data(:, self%jlo)
        enddo
      endif
    case default
      error stop "Unsupported location to apply the bc at in outlet_bc_t%apply_outlet_cell_gradient_bc()"
    endselect

  endsubroutine apply_outlet_primitive_var_bc

  subroutine finalize(self)
    type(outlet_bc_t), intent(inout) :: self
    if(enable_debug_print) call debug_print('Running outlet_bc_t%finalize()', __FILE__, __LINE__)
    if(allocated(self%ilo_ghost)) deallocate(self%ilo_ghost)
    if(allocated(self%ihi_ghost)) deallocate(self%ihi_ghost)
    if(allocated(self%jlo_ghost)) deallocate(self%jlo_ghost)
    if(allocated(self%jhi_ghost)) deallocate(self%jhi_ghost)
  endsubroutine finalize

endmodule mod_outlet_bc

max_issues_repo_name: smillerc/cato

max_forks_repo_name: smillerc/cato


text:
Eric Price

yes that is his name.

He is a student of music at UC Davis. Hes active in many of the musical efforts that this town puts forth.

His instrument is the bass.

Currently Eric is working with the University Symphony Orchestra UC Davis Symphony and serving as Music Manager. Also he plays with local band The Keystones.

Earlier this year Eric performed in the production of the Rocky Horror Show. He played bass in the pit band and wore a pink wig and makeup. The title of pit band is sort of misplaced because the band was actually set into the stage and were an intrigal part of the show.

Please say something. say hi. talk to me. meeting people rocks! and talking to old friends is awesome as well.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      SUBROUTINE SETUPG
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INCLUDE 'SIZES'
      COMMON /EXPONT/ ZS(107),ZP(107),ZD(107)
      COMMON /NATYPE/ NZTYPE(107),MTYPE(30),LTYPE
      COMMON /TEMP/  CC(60,6),ZZ(60,6)
      DIMENSION ALLC(6,5,2), ALLZ(6,5,2)
C     SET-UP THE STEWART'S STO-6G EXPANSIONS
C                                            1S
      ALLZ(1,1,1) =2.310303149D01
      ALLZ(2,1,1) =4.235915534D00
      ALLZ(3,1,1) =1.185056519D00
      ALLZ(4,1,1) =4.070988982D-01
      ALLZ(5,1,1) =1.580884151D-01
      ALLZ(6,1,1) =6.510953954D-02
C
      ALLC(1,1,1) =9.163596280D-03
      ALLC(2,1,1) =4.936149294D-02
      ALLC(3,1,1) =1.685383049D-01
      ALLC(4,1,1) =3.705627997D-01
      ALLC(5,1,1) =4.164915298D-01
      ALLC(6,1,1) =1.303340841D-01
C                                      2S
      ALLZ(1,2,1) =2.768496241D01
      ALLZ(2,2,1) =5.077140627D00
      ALLZ(3,2,1) =1.426786050D00
      ALLZ(4,2,1) =2.040335729D-01
      ALLZ(5,2,1) =9.260298399D-02
      ALLZ(6,2,1) =4.416183978D-02
C
      ALLC(1,2,1) =-4.151277819D-03
      ALLC(2,2,1) =-2.067024148D-02
      ALLC(3,2,1) =-5.150303337D-02
      ALLC(4,2,1) =3.346271174D-01
      ALLC(5,2,1) =5.621061301D-01
      ALLC(6,2,1) =1.712994697D-01
C                                     2P
      ALLZ(1,2,2) =5.868285913D00
      ALLZ(2,2,2) =1.530329631D00
      ALLZ(3,2,2) =5.475665231D-01
      ALLZ(4,2,2) =2.288932733D-01
      ALLZ(5,2,2) =1.046655969D-01
      ALLZ(6,2,2) =4.948220127D-02
C
      ALLC(1,2,2) =7.924233646D-03
      ALLC(2,2,2) =5.144104825D-02
      ALLC(3,2,2) =1.898400060D-01
      ALLC(4,2,2) =4.049863191D-01
      ALLC(5,2,2) =4.012362861D-01
      ALLC(6,2,2) =1.051855189D-01
C                                      3S
      ALLZ(1,3,1) =3.273031938D00
      ALLZ(2,3,1) =9.200611311D-01
      ALLZ(3,3,1) =3.593349765D-01
      ALLZ(4,3,1) =8.636686991D-02
      ALLZ(5,3,1) =4.797373812D-02
      ALLZ(6,3,1) =2.724741144D-02
      ALLC(1,3,1) =-6.775596947D-03
      ALLC(2,3,1) =-5.639325779D-02
      ALLC(3,3,1) =-1.587856086D-01
      ALLC(4,3,1) =5.534527651D-01
      ALLC(5,3,1) =5.015351020D-01
      ALLC(6,3,1) =7.223633674D-02
C                                     3P
      ALLZ(1,3,2) =5.077973607D00
      ALLZ(2,3,2) =1.340786940D00
      ALLZ(3,3,2) =2.248434849D-01
      ALLZ(4,3,2) =1.131741848D-01
      ALLZ(5,3,2) =6.076408893D-02
      ALLZ(6,3,2) =3.315424265D-02
      ALLC(1,3,2) =-3.329929840D-03
      ALLC(2,3,2) =-1.419488340D-02
      ALLC(3,3,2) =1.639395770D-01
      ALLC(4,3,2) =4.485358256D-01
      ALLC(5,3,2) =3.908813050D-01
      ALLC(6,3,2) =7.411456232D-02
C                                     4S
      ALLZ(1,4,1) = 1.365346 D+00
      ALLZ(2,4,1) = 4.393213 D-01
      ALLZ(3,4,1) = 1.877069 D-01
      ALLZ(4,4,1) = 9.360270 D-02
      ALLZ(5,4,1) = 5.052263 D-02
      ALLZ(6,4,1) = 2.809354 D-02
      ALLC(1,4,1) = 3.775056 D-03
      ALLC(2,4,1) =-5.585965 D-02
      ALLC(3,4,1) =-3.192946 D-01
      ALLC(4,4,1) =-2.764780 D-02
      ALLC(5,4,1) = 9.049199 D-01
      ALLC(6,4,1) = 3.406258 D-01
C                                   4P
      ALLC(1,4,2) =-7.052075 D-03
      ALLC(2,4,2) =-5.259505 D-02
      ALLC(3,4,2) =-3.773450 D-02
      ALLC(4,4,2) = 3.874773 D-01
      ALLC(5,4,2) = 5.791672 D-01
      ALLC(6,4,2) = 1.221817 D-01
      ALLZ(1,4,2) = 1.365346 D+00
      ALLZ(2,4,2) = 4.393213 D-01
      ALLZ(3,4,2) = 1.877069 D-01
      ALLZ(4,4,2) = 9.360270 D-02
      ALLZ(5,4,2) = 5.052263 D-02
      ALLZ(6,4,2) = 2.809354 D-02
C                                     5S
      ALLZ(1,5,1) = 7.701420258D-01
      ALLZ(2,5,1) = 2.756268915D-01
      ALLZ(3,5,1) = 1.301847480D-01
      ALLZ(4,5,1) = 6.953441940D-02
      ALLZ(5,5,1) = 4.002545502D-02
      ALLZ(6,5,1) = 2.348388309D-02
      ALLC(1,5,1) = 1.267447151D-02
      ALLC(2,5,1) = 3.266734789D-03
      ALLC(3,5,1) =-4.307553999D-01
      ALLC(4,5,1) =-3.231998963D-01
      ALLC(5,5,1) = 1.104322879D+00
      ALLC(6,5,1) = 4.368498703D-01
C                                      5P
      ALLZ(1,5,2) = 7.701420258D-01
      ALLZ(2,5,2) = 2.756268915D-01
      ALLZ(3,5,2) = 1.301847480D-01
      ALLZ(4,5,2) = 6.953441940D-02
      ALLZ(5,5,2) = 4.002545502D-02
      ALLZ(6,5,2) = 2.348388309D-02
      ALLC(1,5,2) =-1.105673292D-03
      ALLC(2,5,2) =-6.243132446D-02
      ALLC(3,5,2) =-1.628476766D-01
      ALLC(4,5,2) = 3.210328714D-01
      ALLC(5,5,2) = 6.964579592D-01
      ALLC(6,5,2) = 1.493146125D-01
      DO 30 I=1,10
         IF(MTYPE(I).EQ.0)GOTO 30
         NI=MTYPE(I)
         XI=ZS(NI)
         IA=I*4-3
         IB=IA+3
         IF(NI.LT.2) THEN
            NQN=1
         ELSEIF(NI.LT.10)THEN
            NQN=2
         ELSEIF(NI.LT.18)THEN
            NQN=3
         ELSEIF(NI.LT.36)THEN
            NQN=4
         ELSEIF(NI.LT.54)THEN
            NQN=5
         ELSE
            WRITE(6,*)' NO GAUSSIANS AVAILABLE'
            STOP
         ENDIF
         DO 20 K=IA,IB
            L=1
            IF(K.GT.IA) L=2
            IF(K.GT.IA) XI=ZP(NI)
            DO 10 J=1,6
               CC(K,J)=ALLC(J,NQN,L)
   10       ZZ(K,J)=ALLZ(J,NQN,L)*XI**2
   20    CONTINUE
   30 CONTINUE
      RETURN
      END

max_stars_repo_name: openmopac/MOPAC-archive

max_issues_repo_name: openmopac/MOPAC-archive

max_forks_repo_name: openmopac/MOPAC-archive


text:
      SUBROUTINE PHO_EVEINI(Imode,P1,P2,Ip1,Ip2)
C********************************************************************
C
C     prepare /POEVT1/ for new event
C
C     first subroutine called for each event
C
C     input:   P1(4)  particle 1
C              P2(4)  particle 2
C              IMODE  0    general initialization
C                     1    initialization of particles and kinematics
C                     2    initialization after internal rejection
C
C     output:  IP1,IP2  index of interacting particles
C
C********************************************************************
      IMPLICIT NONE
      DOUBLE PRECISION DEPS , EPS , onem , P1 , P2
      INTEGER i , im , Imode , ip , Ip1 , Ip2 , ipcomb , ipos , itmp
      SAVE

      DIMENSION P1(4) , P2(4)

      PARAMETER (EPS=1.D-5,DEPS=1.D-15)

C  input/output channels
      INCLUDE 'inc/poinou'
C  event debugging information
      INCLUDE 'inc/podebg'
C  model switches and parameters
      INCLUDE 'inc/pomdls'
C  general process information
      INCLUDE 'inc/poprcs'
C  gamma-lepton or gamma-hadron vertex information
      INCLUDE 'inc/pofsrc'
C  global event kinematics and particle IDs
      INCLUDE 'inc/pogcms'
C**anfe COMMON not needed here
C  energy-interpolation table
C       include 'inc/potabl'
C  cross sections
      INCLUDE 'inc/pocsec'
C  color string configurations including collapsed strings and hadrons
      INCLUDE 'inc/postrg'
C  standard particle data interface
      INCLUDE 'inc/poevt1'
C  extension to standard particle data interface (PHOJET specific)
      INCLUDE 'inc/poevt2'

C  table of particle indices for recursive PHOJET calls
      INCLUDE 'inc/porecu'
C  event weights and generated cross section
      INCLUDE 'inc/powght'

      DIMENSION im(2)

C  reset debug variables
      KSPom = 0
      KHPom = 0
      KSReg = 0
      KHDir = 0
      KSTrg = 0
      KHTrg = 0
      KSLoo = 0
      KHLoo = 0
      KSDpo = 0
      KSOft = 0
      KHArd = 0
C
      IDNodf = 0
      IDIfr1 = 0
      IDIfr2 = 0
      IDDpom = 0
      ISTr = 0
      IPOix1 = 0
      IF ( ISWmdl(14).GT.0 ) IPOix1 = 1
      IPOix2 = 0
      IPOix3 = 0
C  reset /POEVT1/ and /POEVT2/
      CALL PHO_REGPAR(0,0,0,0,0,0.D0,0.D0,0.D0,0.D0,0,0,0,0,ipos,0)
      CALL PHO_SELCOL(0,0,0,0,0,0,0)
      DO i = 0 , 10
         IPOwgc(i) = 0
      END DO

C  initialization of particle kinematics
C  Initialize elastic scattering with initial kinematics
      IF ( ISWmdl(13).EQ.0 ) THEN
         ipcomb = 0
         IF ( ABS(IFPap(1)).EQ.2212 .AND. IFPap(1).EQ.IFPap(2) ) THEN
            ipcomb = 1
         ELSE IF ( ABS(IFPap(1)).EQ.2212 .AND. IFPap(1).EQ.-IFPap(2) )
     &             THEN
            ipcomb = 2
         ELSE
            IF ( LPRi.GT.4 ) WRITE (LO,*)
     &            "PHO_EVEINI: Unsupported particle" ,
     &           "combination for JLL elastic scattering model."
            CALL PHO_ABORT
         END IF
C* FLUKA has a separate version of this model.
C         CALL JLL_SET(ECM**2, IPCOMB, 1)
      END IF
C  lepton-photon/hadron-photon vertex and initial particles
      im(1) = 0
      im(2) = 0
      IF ( (IPAmdl(11).GT.0) .AND. (IDPsrc(1).NE.0) ) THEN
         CALL PHO_REGPAR(1,IDPsrc(1),IDBsrc(1),0,0,PINi(1,1),PINi(2,1),
     &                   PINi(3,1),PINi(4,1),0,0,0,0,im(1),1)
      ELSE
         CALL PHO_REGPAR(1,IFPap(1),IFPab(1),im(1),0,P1(1),P1(2),P1(3),
     &                   P1(4),0,0,0,0,Ip1,1)
      END IF
      IF ( (IPAmdl(12).GT.0) .AND. (IDPsrc(2).NE.0) ) THEN
         CALL PHO_REGPAR(1,IDPsrc(2),IDBsrc(2),0,0,PINi(1,2),PINi(2,2),
     &                   PINi(3,2),PINi(4,2),0,0,0,0,im(2),1)
      ELSE
         CALL PHO_REGPAR(1,IFPap(2),IFPab(2),im(2),0,P2(1),P2(2),P2(3),
     &                   P2(4),0,0,0,0,Ip2,1)
      END IF
      IF ( (IPAmdl(11).GT.0) .AND. (IDPsrc(1).NE.0) ) THEN
         CALL PHO_REGPAR(1,IDPsrc(1),IDBsrc(1),im(1),0,PFIn(1,1),
     &                   PFIn(2,1),PFIn(3,1),PFIn(4,1),0,10,0,0,ipos,1)
         CALL PHO_REGPAR(1,IFPap(1),IFPab(1),im(1),0,P1(1),P1(2),P1(3),
     &                   P1(4),0,0,0,0,Ip1,1)
      END IF
      IF ( (IPAmdl(12).GT.0) .AND. (IDPsrc(2).NE.0) ) THEN
         CALL PHO_REGPAR(1,IDPsrc(2),IDBsrc(2),im(2),0,PFIn(1,2),
     &                   PFIn(2,2),PFIn(3,2),PFIn(4,2),0,10,0,0,ipos,1)
         CALL PHO_REGPAR(1,IFPap(2),IFPab(2),im(2),0,P2(1),P2(2),P2(3),
     &                   P2(4),0,0,0,0,Ip2,1)
      END IF
      NEVhep = KACcep

      IF ( Imode.LE.1 ) THEN
C  CMS energy
         ECM = SQRT((P1(4)+P2(4))**2-(P1(1)+P2(1))**2-(P1(2)+P2(2))
     &         **2-(P1(3)+P2(3))**2)
C       CALL PHO_PECMS(1,PMASS(1),PMASS(2),ECM,PCM,EE)
         PMAss(1) = PHEp(5,Ip1)
         PVIrt(1) = 0.D0
         IF ( IFPap(1).EQ.22 ) PVIrt(1) = PMAss(1)**2
         PMAss(2) = PHEp(5,Ip2)
         PVIrt(2) = 0.D0
         IF ( IFPap(2).EQ.22 ) PVIrt(2) = PMAss(2)**2
      END IF

C  cross section calculations

      IF ( Imode.NE.1 ) THEN
         ip = 1
         CALL PHO_CSINT(ip,IFPap(1),IFPap(2),IGHel(1),IGHel(2),ECM,
     &                  PVIrt(1),PVIrt(2))
      END IF

      IF ( Imode.LE.0 ) THEN
C  effective cross section
         SIGgen(3) = 0.D0
         IF ( ISWmdl(2).GE.1 ) THEN
            IF ( IPRon(1,1).EQ.1 ) SIGgen(3) = SIGtot - SIGela -
     &           SIGvm(0,0) - SIGcdf(0) - SIGlsd(1) - SIGhsd(1)
     &           - SIGlsd(2) - SIGhsd(2) - SIGldd - SIGhdd - SIGdir
            IF ( IPRon(2,1).EQ.1 ) SIGgen(3) = SIGgen(3) + SIGela
            IF ( IPRon(3,1).EQ.1 ) SIGgen(3) = SIGgen(3) + SIGvm(0,0)
            IF ( IPRon(4,1).EQ.1 ) SIGgen(3) = SIGgen(3) + SIGcdf(0)
            IF ( IPRon(5,1).EQ.1 ) SIGgen(3) = SIGgen(3) + SIGlsd(1)
     &           + SIGhsd(1)
            IF ( IPRon(6,1).EQ.1 ) SIGgen(3) = SIGgen(3) + SIGlsd(2)
     &           + SIGhsd(2)
            IF ( IPRon(7,1).EQ.1 ) SIGgen(3) = SIGgen(3) + SIGldd +
     &           SIGhdd
            IF ( IPRon(8,1).EQ.1 ) SIGgen(3) = SIGgen(3) + SIGdir
C  simulate only hard scatterings
         ELSE
            IF ( IPRon(1,1).EQ.1 ) SIGgen(3) = SIGhar
            IF ( IPRon(8,1).EQ.1 ) SIGgen(3) = SIGgen(3) + SIGdir
         END IF

      END IF

C  reset of mother/daughter relations only (IMODE = 2)

C  debug output
      IF ( IDEb(63).GE.15 ) THEN
         IF ( LPRi.GT.4 ) WRITE (LO,'(/1X,2A,I12,I3)') 'PHO_EVEINI: ' ,
     &        '/POEVT1/ initialized (event/mode)' , KEVent , Imode
         IF ( Imode.LE.0 ) THEN
            IF ( LPRi.GT.4 ) WRITE (LO,'(1X,2A,/,5X,1P,6E12.4)')
     &            'PHO_EVEINI: ' ,
     &        'current suppression factors total-1/2 hard-1/2 diff-1/2:'
     &        , FSUp , FSUh , FSUd
            onem = -1.D0
            itmp = IDEb(57)
            IDEb(57) = MAX(5,itmp)
            CALL PHO_XSECT(1,0,onem)
            IDEb(57) = itmp
         END IF
         CALL PHO_PREVNT(0)
      END IF

      END SUBROUTINE

max_stars_repo_name: pzhristov/DPMJET

max_issues_repo_name: pzhristov/DPMJET

max_forks_repo_name: pzhristov/DPMJET


text:
	SUBROUTINE GSTURB ( szturb, ituwid, iret )
C************************************************************************
C* GSTURB								*
C*									*
C* This subroutine sets the turbulence symbol attributes.		*
C*									*
C* GSTURB ( SZTURB, ITUWID, IRET )					*
C*									*
C* Input parameters:							*
C*	SZTURB		REAL		Turbulence symbol size		*
C*	ITUWID		INTEGER		Turbulence symbol line width	*
C*									*
C* Output parameters:							*
C*	IRET		INTEGER		Return code			*
C**									*
C* Log:									*
C* M. Linda/GSC		 8/96	Based on GSWTHR				*
C************************************************************************
	INCLUDE		'ERROR.PRM'
	INCLUDE		'DEVCHR.CMN'
	INCLUDE		'DEVREQ.CMN'
	INCLUDE		'DEVSET.CMN'
C-----------------------------------------------------------------------
C*
	iret = NORMAL
C
C*	If these are the current requested characteristics, do nothing.
C
	IF  ( szturb .eq. rtursz .and. ituwid .eq. ktuwid ) RETURN
C
C*	Set requested parameters.
C
	IF  ( szturb .gt. 0 ) rtursz = szturb
	IF  ( ituwid .ge. 1 ) ktuwid = ituwid
C
C*	Make changes in device only if the device has been set.
C
	IF  ( ddev .ne. ' ' ) THEN
C
C*	    Set only if requested is different than what is already set.
C
	    IF  ( rtursz .ne. stursz .or. ktuwid .ne. ltuwid)
     +		CALL DSTURB ( rtursz, ktuwid, stursz, ltuwid, iret )
	END IF
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
program testEightByteGet
  use typeSizes
  use netcdf
  implicit none

  ! netcdf related variables
  integer :: ncFileID, dimId, varId

  ! Local variables
  integer, parameter :: numLats = 10
  character (len = *), parameter :: fileName = "eightByteIntExample.nc"
  integer            :: counter
  !integer(kind = FourByteInt), dimension(numLats) :: lats
  integer(kind = EightByteInt), dimension(numLats) :: lats

! -----
  call check(nf90_create(path = trim(fileName), cmode = nf90_clobber, ncid = ncFileID))

  ! Define the dimensions
  call check(nf90_def_dim(ncid = ncFileID, name = "lat", len = numLats, dimid = dimID))
  ! Define the variable
  call check(nf90_def_var(ncFileID, "lat", nf90_int, dimids = dimID, varID = varID) )
  ! Leave define mode
  call check(nf90_enddef(ncfileID))
  ! Write the variable
  lats(:) = (/ (counter, counter = 1, numLats) /)
  print *, lats
  call check(nf90_put_var(ncFileID, varId, lats ))
  call check(nf90_close(ncFileID))

  ! Now read it back in...
  lats(:) = 0
  call check(nf90_open(path = trim(fileName), mode = nf90_nowrite, ncid = ncFileID))
  call check(nf90_inq_varid(ncFileId, "lat", varId))
  call check(nf90_get_var(ncFileID, varId, lats))
  print *, lats
  call check(nf90_close(ncFileID))

contains
  ! Internal subroutine - checks error status after each netcdf, prints out text message each time
  !   an error code is returned.
  subroutine check(status)
    integer, intent ( in) :: status

    if(status /= nf90_noerr) then
      print *, trim(nf90_strerror(status))
    end if
  end subroutine check
end program testEightByteGet

max_stars_repo_licenses: FTL zlib-acknowledgement RSA-MD

max_issues_repo_licenses: FTL zlib-acknowledgement RSA-MD

max_forks_repo_licenses: FTL zlib-acknowledgement RSA-MD


text:
      SUBROUTINE ASFINT (VIS, VISI, VICI, NUMLIN, NGPRIM)

C  ASFINT draws a set of fill areas and fill area sets to illustrate
C  the effect of various interior aspect source flags.

C  Input parameters:
C    VIS    : interior style to be displayed/suppressed
C    VISI   : interior style index to be displayed/suppressed
C    VICI   : interior color index to be displayed/suppressed
C    NUMLIN : number of primitives to display
C  Output parameters:
C    NGPRIM : which primitive is drawn differently

      COMMON /GLOBNU/ CTLHND, ERRSIG, ERRFIL, IERRCT, UNERR,
     1        TESTCT, IFLERR, PASSSW, ERRSW, MAXLIN,
     2        CONID, MEMUN, WKID, WTYPE, GLBLUN, INDLUN,
     3        DUMINT, DUMRL
      INTEGER         CTLHND, ERRSIG, ERRFIL, IERRCT, UNERR,
     1        TESTCT, IFLERR, PASSSW, ERRSW, MAXLIN,
     2        CONID, MEMUN, WKID, WTYPE, GLBLUN, INDLUN,
     3        DUMINT(20), ERRIND
      REAL    DUMRL(20)

C aspect identifier
      INTEGER    PTXCI,    PIS,      PISI,     PICI,     PEDFG
      PARAMETER (PTXCI=10, PIS  =11, PISI =12, PICI =13, PEDFG=14)

C aspect source
C                bundled     individual
      INTEGER    PBUNDL,     PINDIV
      PARAMETER (PBUNDL = 0, PINDIV = 1)

C interior style
      INTEGER    PHOLLO,   PSOLID,   PPATTR,   PHATCH,   PISEMP
      PARAMETER (PHOLLO=0, PSOLID=1, PPATTR=2, PHATCH=3, PISEMP=4)
C off/on switch for edge flag and error handling mode
      INTEGER    POFF,     PON
      PARAMETER (POFF = 0, PON = 1)

C aspect visibility
      INTEGER    DISPLY,     SUPPRS
      PARAMETER (DISPLY = 0, SUPPRS = 1)

      INTEGER    VIS  (DISPLY:SUPPRS), NPTS(1), IBSW, IPRIM, NEXT
      INTEGER    VISI (DISPLY:SUPPRS), VICI (DISPLY:SUPPRS)
      INTEGER    IX, NUMLIN,NGPRIM, RNDINT, STYSRC,IIXSRC,COLSRC
      INTEGER    P1(10),P2(10),P3(10)

      REAL       XA(10),YA(10),ZA(10), YTOP,YLOC,YINCR

C  use IBSW to switch individual/bundled
      IBSW = PINDIV+PBUNDL
C  primitive to be drawn incorrectly
      NGPRIM = RNDINT (1,NUMLIN)
      YINCR = 1.0 / (NUMLIN+1)
      YTOP  = 1-YINCR
      YLOC = YTOP
      CALL NUMLAB (NUMLIN, 0.15, YTOP, YINCR)
      IPRIM = 1
      XA(1) = 0.2
      XA(2) = 0.2
      XA(3) = 0.8
      XA(4) = 0.8
      ZA(1) = 0.5
      ZA(2) = 0.5
      ZA(3) = 0.5
      ZA(4) = 0.5
C  set edge flag = OFF
      CALL PSIASF (PEDFG, PINDIV)
      CALL PSEDFG (POFF)
      NPTS(1) = 4
C  ensure both individual and bundle are picked for each aspect
      CALL RNPERM (NUMLIN, P1)
      CALL RNPERM (NUMLIN, P2)
      CALL RNPERM (NUMLIN, P3)

C  loop to draw primitives
      DO 100 IX = 1,NUMLIN
C  randomly pick source for aspects:
         STYSRC = MOD (P1(IX), 2)
         IIXSRC = MOD (P2(IX), 2)
         COLSRC = MOD (P3(IX), 2)

C  set bundle #ix to displayed aspects, for BUNDLED source
C               and suppressed aspects, for INDIVIDUAL source
         CALL PSIR (WKID, IX, VIS(STYSRC), VISI(IIXSRC), VICI(COLSRC))
         CALL PSII (IX)
C  set individual attributes
C                    to displayed aspects, for INDIVIDUAL source
C                  and suppressed aspects, for BUNDLED source
         CALL PSIS  (VIS (IBSW-STYSRC))
         CALL PSISI (VISI(IBSW-IIXSRC))
         CALL PSICI (VICI(IBSW-COLSRC))
C  set ASF flags in accordance with: stysrc, iixsrc, colsrc
         CALL PSIASF (PIS,  STYSRC)
         CALL PSIASF (PISI, IIXSRC)
         CALL PSIASF (PICI, COLSRC)

         IF (IX .EQ. NGPRIM) THEN
C  set ASF flag NOT in accordance with typsrc
            CALL PSIASF (PIS,  IBSW-STYSRC)
            CALL PSIASF (PISI, IBSW-IIXSRC)
            CALL PSIASF (PICI, IBSW-COLSRC)
            NEXT = IPRIM
         ELSE
            NEXT = IPRIM+1
         ENDIF

         YA(1) = YLOC + 0.3*YINCR
         YA(2) = YLOC - 0.3*YINCR
         YA(3) = YLOC - 0.3*YINCR
         YA(4) = YLOC + 0.3*YINCR

         IF (IPRIM .EQ. 1) THEN
            CALL PFA (4, XA,YA)
         ELSEIF (IPRIM .EQ. 2) THEN
            CALL PFA3 (4, XA,YA,ZA)
         ELSEIF (IPRIM .EQ. 3) THEN
            CALL PFAS (1, NPTS, XA,YA)
         ELSEIF (IPRIM .EQ. 4) THEN
            CALL PFAS3 (1, NPTS, XA,YA,ZA)
         ELSE
            CALL PFA (4, XA,YA)
         ENDIF
         YLOC = YLOC-YINCR
         IPRIM = NEXT

C  next ix
100   CONTINUE
      END


program main

!*****************************************************************************80
!
!! MAIN is the main program for SOBOL_PRB.
!
!  Discussion:
!
!    SOBOL_PRB tests the SOBOL library.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    16 January 2010
!
!  Author:
!
!    John Burkardt

implicit real (kind = 8)(a-h, o-z)

  integer ( kind = 8 ), parameter :: ndim = 4

  integer ( kind = 8 ) dim_num
  integer ( kind = 8 ), dimension ( 3 ) :: dim_num_test = (/ 100, 500, 1000 /)
  integer ( kind = 8 ) i
  real ( kind = 8 ), allocatable, dimension ( : ) :: r
  integer ( kind = 8 ), dimension(ndim) ::  seed
  integer ( kind = 8 ) seed_in
  integer ( kind = 8 ) seed_out
  integer ( kind = 8 ) test
      real*8 gasdev
      real*8, dimension(ndim) :: q,x

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST13'
  write ( *, '(a)' ) '  I8_SOBOL computes the next element of a Sobol sequence.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  In this test, we get a few samples at high dimension.'
  write ( *, '(a)' ) '  We only print the first and last 2 entries of each'
  write ( *, '(a)' ) '  sample.'

    dim_num = 2

      open(100, file = 'random.dat')

    write ( *, '(a)' ) ' '
    write ( *, '(a,i8)' ) '  Using dimension    ', ndim

    seed = (/1,2,3,4/)

      seed_in = 0

      pi = 4.0*atan(1.0)

!      ns = 100000 ! sampling points

      print *, 'how many sampling points you want?'
      read(5,*) ns

      z0 = 0d0
      q = (/0d0,0d0,0d0,0d0/)

      anm = (1d0/2d0/pi)**2

    do i = 0, ns

!      do j=1,ndim
!        x(j) = gasdev(seed_in) + q(j)
!      enddo

      call i8_sobol(ndim, seed_in, x)
      x = 12.0*(x-0.5d0)

      z0 = z0 + x(1)**2*x(2)**2*exp(-x(1)**2/2.0-x(2)**2/2d0-x(3)**2/2d0-x(4)**2/2d0)
    enddo

    write ( 6, * )  'integration  = ', z0*anm/dble(ns)*12.0**4

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Normal end of execution.'
  write ( *, '(a)' ) ' '
  call timestamp ( )

1000 format(2(f14.6,1x))
  stop
end


text:
c  median smooth a spectrum with a window of diameter = iwidth pixels
c  the original spectrum is preserved.

c  11.08.02 - fixed an error, it was setting sout(j) rather than sout(i)
c     this increased flattopping of emission lines, etc, and was worse
c     with larger smoothing

      subroutine medsmooth(n,spin,iwidth,sout)

      real spin(n),sout(n)
      integer iwidth
      real temparr(n)
      external select

      include 'pcredshift.h'

c  if smoothing is 1 or 0 pixels, array is unchanged
      if (iwidth .le. 1) then
         do i=1,n
            sout(i) = spin(i)
         end do
         return
      end if

c compute radii for smoothing window
c if iwidth is odd we get a symmetric window, if it's even,
c truncate low side down, to
c make the high side the long side (no particular reason)
      irlow = int((iwidth-1)/2)
      irhigh = iwidth-irlow-1

      do i=1,n
         i1=max(1,i-irlow)
         i2=min(n,i+irhigh)
         np=0
         do j=i1,i2
            if(spin(j) .gt. bad) then
               np = np+1
               temparr(np) = spin(j)
            end if
         end do
         if (np .eq. 0) then
            sout(i) = badset
         else if (np .eq. 1) then
            sout(i) = temparr(1)
         else if (np .eq. 2) then
            sout(i) = (temparr(1) + temparr(2))/2.0
         else
c do the median right for once
            itmp = np/2
            if (mod(np,2) .ne. 0) then
               sout(i) = select(itmp+1,np,temparr)
            else
               sout(i) = (select(itmp,np,temparr) +
     $              select(itmp+1,np,temparr)) / 2.0
            end if
         end if
      end do

      return
      end