text:
Module mem_sfcanl

type sfcanl_fields
real, pointer, dimension(:,:)    ::  u,v,pr,t,td,vv,rh,pm,tad  &
                                    ,th,the,hi,ps,vor,qm,qcon,div,thad  &
                                    ,qad,spd,css,vis,fwx,tgd
end type

type(sfcanl_fields) :: sfcanl

! Pointers for renaming arrays in lapsvanl
real, pointer, dimension(:,:) :: &
         u_a,v_a,p_a,t,td,vv,rh,hi,mslp,tadv,theta,thetae,psfc  &
        ,vort,q,qcon,div,thadv,qadv,spd,cssi,vis,fire,tgd_k

Contains

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine alloc_sfcanl_arrays (nxl, nyl)

implicit none
integer :: nxl, nyl

   allocate(sfcanl%u    (nxl,nyl))
   allocate(sfcanl%v    (nxl,nyl))
   allocate(sfcanl%pr   (nxl,nyl))
   allocate(sfcanl%t    (nxl,nyl))
   allocate(sfcanl%td   (nxl,nyl))
   allocate(sfcanl%vv   (nxl,nyl))
   allocate(sfcanl%rh   (nxl,nyl))
   allocate(sfcanl%pm   (nxl,nyl))
   allocate(sfcanl%tad  (nxl,nyl))
   allocate(sfcanl%th   (nxl,nyl))
   allocate(sfcanl%the  (nxl,nyl))
   allocate(sfcanl%hi   (nxl,nyl))
   allocate(sfcanl%ps   (nxl,nyl))
   allocate(sfcanl%vor  (nxl,nyl))
   allocate(sfcanl%qm   (nxl,nyl))
   allocate(sfcanl%qcon (nxl,nyl))
   allocate(sfcanl%div  (nxl,nyl))
   allocate(sfcanl%thad (nxl,nyl))
   allocate(sfcanl%qad  (nxl,nyl))
   allocate(sfcanl%spd  (nxl,nyl))
   allocate(sfcanl%css  (nxl,nyl))
   allocate(sfcanl%vis  (nxl,nyl))
   allocate(sfcanl%fwx  (nxl,nyl))
   allocate(sfcanl%tgd  (nxl,nyl))

return
end subroutine alloc_sfcanl_arrays

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine point_sfcanl_arrays ()

implicit none

u_a   =>  sfcanl%u
v_a   =>  sfcanl%v
p_a   =>  sfcanl%pr
t     =>  sfcanl%t
td    =>  sfcanl%td
vv    =>  sfcanl%vv
rh    =>  sfcanl%rh
hi    =>  sfcanl%pm
mslp  =>  sfcanl%tad
tadv  =>  sfcanl%th
theta =>  sfcanl%the
thetae=>  sfcanl%hi
psfc  =>  sfcanl%ps
vort  =>  sfcanl%vor
q     =>  sfcanl%qm
qcon  =>  sfcanl%qcon
div   =>  sfcanl%div
thadv =>  sfcanl%thad
qadv  =>  sfcanl%qad
spd   =>  sfcanl%spd
cssi  =>  sfcanl%css
vis   =>  sfcanl%vis
fire  =>  sfcanl%fwx
tgd_k =>  sfcanl%tgd

return
end subroutine point_sfcanl_arrays

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine deallocate_sfcanl_arrays()

implicit none

   if (associated(sfcanl%u    )) deallocate(sfcanl%u    )
   if (associated(sfcanl%v    )) deallocate(sfcanl%v    )
   if (associated(sfcanl%pr   )) deallocate(sfcanl%pr   )
   if (associated(sfcanl%t    )) deallocate(sfcanl%t    )
   if (associated(sfcanl%td   )) deallocate(sfcanl%td   )
   if (associated(sfcanl%vv   )) deallocate(sfcanl%vv   )
   if (associated(sfcanl%rh   )) deallocate(sfcanl%rh   )
   if (associated(sfcanl%pm   )) deallocate(sfcanl%pm   )
   if (associated(sfcanl%tad  )) deallocate(sfcanl%tad  )
   if (associated(sfcanl%th   )) deallocate(sfcanl%th   )
   if (associated(sfcanl%the  )) deallocate(sfcanl%the  )
   if (associated(sfcanl%hi   )) deallocate(sfcanl%hi   )
   if (associated(sfcanl%ps   )) deallocate(sfcanl%ps   )
   if (associated(sfcanl%vor  )) deallocate(sfcanl%vor  )
   if (associated(sfcanl%qm   )) deallocate(sfcanl%qm   )
   if (associated(sfcanl%qcon )) deallocate(sfcanl%qcon )
   if (associated(sfcanl%div  )) deallocate(sfcanl%div  )
   if (associated(sfcanl%thad )) deallocate(sfcanl%thad )
   if (associated(sfcanl%qad  )) deallocate(sfcanl%qad  )
   if (associated(sfcanl%spd  )) deallocate(sfcanl%spd  )
   if (associated(sfcanl%css  )) deallocate(sfcanl%css  )
   if (associated(sfcanl%vis  )) deallocate(sfcanl%vis  )
   if (associated(sfcanl%fwx  )) deallocate(sfcanl%fwx  )
   if (associated(sfcanl%tgd  )) deallocate(sfcanl%tgd  )

return
end subroutine deallocate_sfcanl_arrays

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine nullify_sfcanl_arrays ()

implicit none

   if (associated(sfcanl%u    )) nullify(sfcanl%u    )
   if (associated(sfcanl%v    )) nullify(sfcanl%v    )
   if (associated(sfcanl%pr   )) nullify(sfcanl%pr   )
   if (associated(sfcanl%t    )) nullify(sfcanl%t    )
   if (associated(sfcanl%td   )) nullify(sfcanl%td   )
   if (associated(sfcanl%vv   )) nullify(sfcanl%vv   )
   if (associated(sfcanl%rh   )) nullify(sfcanl%rh   )
   if (associated(sfcanl%pm   )) nullify(sfcanl%pm   )
   if (associated(sfcanl%tad  )) nullify(sfcanl%tad  )
   if (associated(sfcanl%th   )) nullify(sfcanl%th   )
   if (associated(sfcanl%the  )) nullify(sfcanl%the  )
   if (associated(sfcanl%hi   )) nullify(sfcanl%hi   )
   if (associated(sfcanl%ps   )) nullify(sfcanl%ps   )
   if (associated(sfcanl%vor  )) nullify(sfcanl%vor  )
   if (associated(sfcanl%qm   )) nullify(sfcanl%qm   )
   if (associated(sfcanl%qcon )) nullify(sfcanl%qcon )
   if (associated(sfcanl%div  )) nullify(sfcanl%div  )
   if (associated(sfcanl%thad )) nullify(sfcanl%thad )
   if (associated(sfcanl%qad  )) nullify(sfcanl%qad  )
   if (associated(sfcanl%spd  )) nullify(sfcanl%spd  )
   if (associated(sfcanl%css  )) nullify(sfcanl%css  )
   if (associated(sfcanl%vis  )) nullify(sfcanl%vis  )
   if (associated(sfcanl%fwx  )) nullify(sfcanl%fwx  )
   if (associated(sfcanl%tgd  )) nullify(sfcanl%tgd  )

return
end subroutine nullify_sfcanl_arrays

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

end Module

max_stars_repo_name: maxinye/laps-mirror

max_stars_repo_licenses: Intel Unlicense OLDAP-2.2.1 NetCDF

max_issues_repo_name: longwosion/laps-mirror

max_issues_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense

max_forks_repo_name: longwosion/laps-mirror

max_forks_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense


text:
c     Subroutine to find CDM qc code
c     AJ_Kettle, 11Dec2019

      SUBROUTINE find_cdm_qc_code(l_rgh,l_lines,
     +    s_vec_airt_qc_flag,s_vec_dewp_qc_flag,
     +    s_vec_stnp_qc_flag,s_vec_slpr_qc_flag,
     +    s_vec_wdir_qc_flag,s_vec_wspd_qc_flag,

     +    s_vec_airt_cdmqc,s_vec_dewp_cdmqc,
     +    s_vec_stnp_cdmqc,s_vec_slpr_cdmqc,
     +    s_vec_wdir_cdmqc,s_vec_wspd_cdmqc)

      IMPLICIT NONE
c************************************************************************
c     Variables passed into subroutine

      INTEGER             :: l_rgh
      INTEGER             :: l_lines

      CHARACTER(LEN=8)    :: s_vec_airt_qc_flag(l_rgh)
      CHARACTER(LEN=8)    :: s_vec_dewp_qc_flag(l_rgh)
      CHARACTER(LEN=8)    :: s_vec_stnp_qc_flag(l_rgh)
      CHARACTER(LEN=8)    :: s_vec_slpr_qc_flag(l_rgh)
      CHARACTER(LEN=8)    :: s_vec_wdir_qc_flag(l_rgh)
      CHARACTER(LEN=8)    :: s_vec_wspd_qc_flag(l_rgh)

c     CDM qc flag
      CHARACTER(LEN=1)    :: s_vec_airt_cdmqc(l_rgh)
      CHARACTER(LEN=1)    :: s_vec_dewp_cdmqc(l_rgh)
      CHARACTER(LEN=1)    :: s_vec_stnp_cdmqc(l_rgh)
      CHARACTER(LEN=1)    :: s_vec_slpr_cdmqc(l_rgh)
      CHARACTER(LEN=1)    :: s_vec_wdir_cdmqc(l_rgh)
      CHARACTER(LEN=1)    :: s_vec_wspd_cdmqc(l_rgh)
c*****
c     Variables used within subroutine

      INTEGER             :: i,j,k,ii,jj,kk

      INTEGER             :: i_cnt_fail_airt
      INTEGER             :: i_cnt_fail_dewp
      INTEGER             :: i_cnt_fail_stnp
      INTEGER             :: i_cnt_fail_slpr
      INTEGER             :: i_cnt_fail_wdir
      INTEGER             :: i_cnt_fail_wspd

c************************************************************************
c      print*,'just entering find_cdm_qc_code'

      i_cnt_fail_airt=0
      i_cnt_fail_dewp=0
      i_cnt_fail_stnp=0
      i_cnt_fail_slpr=0
      i_cnt_fail_wdir=0
      i_cnt_fail_wspd=0

      DO i=1,l_lines

c      Initialize flags to missing
       s_vec_airt_cdmqc(i)='3'
       s_vec_dewp_cdmqc(i)='3'
       s_vec_stnp_cdmqc(i)='3'
       s_vec_slpr_cdmqc(i)='3'
       s_vec_wdir_cdmqc(i)='3'
       s_vec_wspd_cdmqc(i)='3'

c      AIRT
       IF (LEN_TRIM(s_vec_airt_qc_flag(i)).EQ.0) THEN
        s_vec_airt_cdmqc(i)='0'   !if no flag then pass=0
       ENDIF
       IF (LEN_TRIM(s_vec_airt_qc_flag(i)).GT.0) THEN
        s_vec_airt_cdmqc(i)='1'   !if any flag then fail=1
        i_cnt_fail_airt=i_cnt_fail_airt+1
       ENDIF

c      DEWP
       IF (LEN_TRIM(s_vec_dewp_qc_flag(i)).EQ.0) THEN
        s_vec_dewp_cdmqc(i)='0'   !if no flag then pass=0
       ENDIF
       IF (LEN_TRIM(s_vec_dewp_qc_flag(i)).GT.0) THEN
        s_vec_dewp_cdmqc(i)='1'   !if any flag then fail=1
        i_cnt_fail_dewp=i_cnt_fail_dewp+1
       ENDIF

c      STNP
       IF (LEN_TRIM(s_vec_stnp_qc_flag(i)).EQ.0) THEN
        s_vec_stnp_cdmqc(i)='0'   !if no flag then pass=0
       ENDIF
       IF (LEN_TRIM(s_vec_stnp_qc_flag(i)).GT.0) THEN
        s_vec_stnp_cdmqc(i)='1'   !if any flag then fail=1
        i_cnt_fail_stnp=i_cnt_fail_stnp+1
       ENDIF

c      SLPR
       IF (LEN_TRIM(s_vec_slpr_qc_flag(i)).EQ.0) THEN
        s_vec_slpr_cdmqc(i)='0'   !if no flag then pass=0
       ENDIF
       IF (LEN_TRIM(s_vec_slpr_qc_flag(i)).GT.0) THEN
        s_vec_slpr_cdmqc(i)='1'   !if any flag then fail=1
        i_cnt_fail_slpr=i_cnt_fail_slpr+1
       ENDIF

c      WDIR
       IF (LEN_TRIM(s_vec_wdir_qc_flag(i)).EQ.0) THEN
        s_vec_wdir_cdmqc(i)='0'   !if no flag then pass=0
       ENDIF
       IF (LEN_TRIM(s_vec_wdir_qc_flag(i)).GT.0) THEN
        s_vec_wdir_cdmqc(i)='1'   !if any flag then fail=1
        i_cnt_fail_wdir=i_cnt_fail_wdir+1
       ENDIF

c      WSPD
       IF (LEN_TRIM(s_vec_wspd_qc_flag(i)).EQ.0) THEN
        s_vec_wspd_cdmqc(i)='0'   !if no flag then pass=0
       ENDIF
       IF (LEN_TRIM(s_vec_wspd_qc_flag(i)).GT.0) THEN
        s_vec_wspd_cdmqc(i)='1'   !if any flag then fail=1
        i_cnt_fail_wspd=i_cnt_fail_wspd+1
       ENDIF

      ENDDO

c      print*,'i_cnt_fail_airt',i_cnt_fail_airt
c      print*,'i_cnt_fail_dewp',i_cnt_fail_dewp
c      print*,'i_cnt_fail_stnp',i_cnt_fail_stnp
c      print*,'i_cnt_fail_slpr',i_cnt_fail_slpr
c      print*,'i_cnt_fail_wdir',i_cnt_fail_wdir
c      print*,'i_cnt_fail_wspd',i_cnt_fail_wspd

c      print*,'just leaving find_cdm_qc_code'

c      STOP 'find_cdm_qc_code'

      RETURN
      END

ext: f

max_stars_repo_name: glamod/glamod-nuim

max_issues_repo_name: glamod/glamod-nuim

max_forks_repo_name: glamod/glamod-nuim


text:
      module mo_ver_opts

      private
      public :: ver_opts

      contains

      subroutine ver_opts( options, model, machine, march, arch_type, &
                           wrk_dir, cpp_dir, cpp_opts, subfile, diagprnt, &
                           tavgprnt, cpucnt, vec_ftns )
!-----------------------------------------------------------------------
!   	... Set the simulation options
!-----------------------------------------------------------------------

      use io

      implicit none

!-----------------------------------------------------------------------
!   	... Dummy args
!-----------------------------------------------------------------------
      integer, intent(inout) ::      cpucnt
      character(len=8), intent(inout)  ::  machine
      character(len=8), intent(inout)  ::  model
      character(len=8), intent(inout)  ::  march
      character(len=8), intent(inout)  ::  arch_type
      character(len=64), intent(out)   ::  wrk_dir
      character(len=64), intent(inout) ::  cpp_dir, cpp_opts
      character(len=64), intent(out)   ::  subfile
      logical, intent(inout)           ::  diagprnt
      logical, intent(inout)           ::  tavgprnt
      logical, intent(inout)           ::  vec_ftns
      logical, intent(inout)           ::  options(:)

!-----------------------------------------------------------------------
!   	... Local variables
!-----------------------------------------------------------------------
      integer, parameter :: maxparm = 26

      integer      :: kpar, nchar, k
      integer      :: err
      logical      :: entered(maxparm)

      character(len=20) :: keywrd
      character(len=20) :: parkey(maxparm)

      logical      :: found

      integer  ::  lenof

      parkey(1) = 'MACHINE'
      parkey(2) = 'DIFFUSION'
      parkey(3) = 'CONVECTION'
      parkey(4) = 'NORMS'
      parkey(5) = 'CONSERVATION'
      parkey(6) = 'SOURCECODE'
      parkey(7) = 'CPUS'
      parkey(8) = 'MULTITASK'
      parkey(9) = 'FIXER'
      parkey(10) = 'DIAGPRNT'
      parkey(11) = 'RXTNLOOKUP'
      parkey(12) = 'RELHUM'
      parkey(13) = 'F90'
      parkey(14) = 'GEOHEIGHT'
      parkey(15) = 'USERHOOK'
      parkey(16) = 'MODULES'
      parkey(17) = 'WORKDIR'
      parkey(18) = 'NAMEMOD'
      parkey(19) = 'SUBFILE'
      parkey(20) = 'TAVGPRNT'
      parkey(21) = 'VEC_FTNS'
      parkey(22) = 'ARCHITECTURE'
      parkey(23) = 'CPP_DIR'
      parkey(24) = 'CPP_OPTS'
      parkey(25) = 'MODEL'
      parkey(26) = 'MODEL_ARCHITECTURE'

      entered = .false.

!-----------------------------------------------------------------------
!   	... Scan for valid option keyword
!-----------------------------------------------------------------------
      do
         call cardin( lin, buff, nchar )
	 buffh = buff
	 call upcase ( buffh )
         if( buffh == 'ENDVERSIONOPTIONS' ) then
	    if( .not. options(1) ) then       ! not a cray target
	       if( cpucnt > 1 ) then
	          options(10) = .true.        ! "distributed" processing
	       else
	          options(10) = .false.       ! no dist processing
	       end if
	    end if
	    if( .not. options(13) ) then      ! if not fortran90 then no modules
	       options(16) = .false.
	       options(17) = .false.
	    end if
	    if( machine /= 'IBM' ) then
	       vec_ftns = .false.
	    end if
	    if( machine == 'NEC' ) then
               march = 'VECTOR'
	    else if( machine == 'INTEL' ) then
               march = 'SCALAR'
	    end if
	    if( model == 'WRF' ) then
               march = 'SCALAR'
	    end if
            exit
	 end if
	 k = index( buffh(:nchar), '=' )
         if( k /= 0 ) then
	    keywrd = buffh(:k-1)
	    found  = .false.
            do kpar = 1,maxparm
               if( keywrd == parkey(kpar) ) then
		  found = .true.
	          exit
	       end if
	    end do
	 else
            call errmes ( ' option specification has no = operator@', lout, buff, 1, buff )
         end if
	 if( .not. found ) then
!-----------------------------------------------------------------------
!  	... Invalid parameter keyword; terminate the program
!-----------------------------------------------------------------------
            call errmes ( ' # is an invalid options' &
                       // ' parameter keyword@', lout, keywrd, &
                          LENOF(20,keywrd), buff )
         end if

!-----------------------------------------------------------------------
!     	... Valid parameter keyword; now check for duplicate keyword
!-----------------------------------------------------------------------
         if( entered(kpar) ) then
            call errmes( '0 *** # has already been specified@', &
                          lout, parkey(kpar), k, ' ' )
         end if

!-----------------------------------------------------------------------
!     	... Set individual options
!-----------------------------------------------------------------------
         if( kpar == 1 ) then
	    machine = buffh(k+1:nchar)
	    if( machine /= 'CRAYYMP' .and. machine /= 'CRAY2' .and. machine /= 'CRAY3' &
	                             .and. machine /= 'J90' .and. machine /= 'C90' ) then
	       options(1) = .false.
	    end if
         else if( kpar == 6 ) then
	    if( buffh(k+1:nchar) /= 'FULL' ) then
	       options(6) = .false.
	    end if
         else if( kpar == 7 ) then
	    call intcon( buffh(k+1:), &
                         nchar - k, &
                         cpucnt, &
                         err )
	    if( err /= 0 ) then
	       call errmes( ' # is not a valid number@', &
                            lout, &
                            buffh(k+1:), &
                            nchar - k, &
                            buff )
	    end if
         else if( kpar == 21 ) then
	    if( buffh(k+1:nchar) == 'ON' .or. buffh(k+1:nchar) == 'YES' ) then
	       vec_ftns = .true.
	    end if
         else if( kpar == 22 ) then
	    arch_type = buffh(k+1:nchar)
	    if( arch_type /= 'OMP' .and. arch_type /= 'MPI' .and. arch_type /= 'HYBRID' ) then
               call errmes( '0 *** # is not a valid architecture type@', lout, arch_type, 8, ' ' )
	    end if
         else if( kpar == 25 ) then
	    model = buffh(k+1:nchar)
	    if( model /= 'MOZART' .and. model /= 'CAM' .and. model /= 'WRF' ) then
               call errmes( '# is not a valid model type@', lout, buff(k+1:), nchar-k , ' ' )
	    end if
         else if( kpar == 26 ) then
	    march = buffh(k+1:nchar)
	    if( march /= 'SCALAR' .and. march /= 'CACHE' .and. march /= 'VECTOR' ) then
               call errmes( '# is not a valid model architecture type@', lout, buff(k+1:), nchar-k , ' ' )
	    end if
	 else
	    if( buffh(k+1:nchar) == 'ON' .or. buffh(k+1:nchar) == 'YES' ) then
	       if( kpar == 10) then
		  diagprnt = .true.
	       else if( kpar == 20) then
		  tavgprnt = .true.
	       else if( kpar == 18 ) then
	          options(17) = .true.
	       else if( kpar /= 8 .and. kpar /= 9 ) then
	          options(kpar) = .true.
	       else if( kpar == 8 ) then
	          options(10) = .true.
	       else if( kpar == 9 ) then
	          options(9) = .true.
	       end if
	    else
	       if( kpar == 17) then
		  wrk_dir = buff(k+1:nchar)
	       else if( kpar == 23) then
		  cpp_dir = buff(k+1:nchar)
	       else if( kpar == 24) then
		  cpp_opts = buff(k+1:nchar)
	       else if( kpar == 19) then
		  subfile = buff(k+1:nchar)
	       else if( kpar == 10) then
		  diagprnt = .false.
	       else if( kpar == 20) then
		  tavgprnt = .false.
	       else if( kpar /= 8 .and. kpar /= 9 ) then
	          options(kpar) = .false.
	       else if( kpar == 8 ) then
	          options(10) = .false.
	       else if( kpar == 9 ) then
	          options(9) = .false.
	       end if
	    end if
	 end if
	 entered(kpar) = .true.
      end do

      end subroutine ver_opts

      end module mo_ver_opts

max_stars_repo_licenses: NCSA Unlicense

max_issues_repo_licenses: NCSA Unlicense

max_forks_repo_licenses: NCSA Unlicense


text:
module opacity_co
   ! -----------------------------------------------------------------
   ! This module contains all variables relating to the CO opacity
   ! tables. Adapted from Eldridge implemented by SdM
   ! -----------------------------------------------------------------
   ! CO_MT:          opacity table size in log10 T direction
   ! CO_MR:          "" in log R = rho^3/T6 direction
   ! CO_MH:          "" in Hyrdrogen abundance direction
   ! CO_CO:          "" in Carbon and oxygen direction
   ! SPLINE_OPAC_CO: Array for spline coefficients
   ! CBASE, OBASE:   C and O abundance because of metallicity only
   ! opT, opR:       log T and log R on which opacity table is defined
   ! opH, COcompos:  Compositions on which opacity table is defined
   ! -----------------------------------------------------------------

   use real_kind

   implicit none
   integer, parameter :: co_mh = 5
   integer, parameter :: co_mt = 141
   integer, parameter :: co_mr = 31
   integer, parameter :: co_co = 8

   real(double) :: cbase,obase

   real(double) :: opt(co_mt), opr(co_mr)
   real(double), parameter :: oph(co_mh) = &
        (/0., 0.03, 0.1, 0.35, 0.7/)
   real(double), parameter :: cocompos(co_co) =  &
        (/0.0d0, 0.01d0, 0.03d0, 0.1d0, 0.2d0, 0.4d0, 0.6d0, 1.0d0/)

   real(double), save, allocatable :: spline_opac_co(:,:,:,:,:)

contains

! As OPSPLN for the CO-enhanced tables SdM
! Calculate a bicubic spline interpolation fit for the temperature
! and density opacity fit.  This is based on routines from Princeton
! Splines, D. McCune found by L. Dray.
subroutine opspln_co
   use real_kind
   use settings
   use opacity_spline

   implicit none
   real(double) ::  mat(4,co_mt)
   integer :: i, j, k, ic, ic1, ic2

   ! define T and R intervals in opac tables
   ! NB: R = rho/T_6^4 and not rho as in the ld tables
   do i=1,co_mt
      opt(i) = 3.0d0 + 5.0d-2*dble(i-1)
   end do
   do j=1,co_mr
      opr(j) = -8.0d0 + 5.0d-1*dble(j-1)
   end do
   ! Setup CO spline tables
   do k=1,305
      do j=1,co_mr
         do i=1,co_mt
            mat(1,i) = spline_opac_co(1,1,i,j,k)
         end do
         ! Construct splines in the T direction
         call spline(co_mt, opt, mat)
         do i=1,co_mt-1
            do ic1=2,4
               spline_opac_co(ic1,1,i,j,k) = mat(ic1,i)
            end do
         end do
      end do

      ! Construct splines in the R (= rho/T_6^3) direction
      do i=1,co_mt-1
         do ic = 1,4
            do j=1,co_mr
               mat(1,j) = spline_opac_co(ic,1,i,j,k)
            end do
            mat(2,1)  = 0.0d0
            mat(3,1)  = 0.0d0
            mat(2,co_mr) = 0.0d0
            mat(3,co_mr) = 0.0d0

            call spline(co_mr, opr, mat)
            do j=1,co_mr-1
               do ic2 = 2,4
                  spline_opac_co(ic,ic2,i,j,k) = mat(ic2,j)
               end do         !IC2
            end do            !J
         end do               !IC
      end do                  !I
   end do                     !K
end subroutine opspln_co

! Read CO enhanced opacity tables
subroutine load_opacity_co(co_file)
   use real_kind
   use settings
   use constants

   implicit none
   integer :: i,j,k
   real(double) :: temp
   integer, intent(in) :: co_file

   if (.not. allocated(spline_opac_co)) allocate(spline_opac_co(4,4,co_mt,co_mr,305))

   read (co_file,'(e10.2)') czs

   do k=1,305
      read (co_file,*)
      do i=1,co_mt
         read (co_file,'(f5.2, 31f7.3)') temp, (spline_opac_co(1,1,i,j,k), j=1,co_mr)
      end do
   end do
   rewind (co_file)

   ! Calculate initial H and basic CO abundance as scaled from solar
   ! Don't override CH value from init.dat (CH>0)!
   if (ch<0.0d0) ch = 0.76d0 - 3.0d0*czs
   cbase = czs*cc
   obase = czs*co

   ! log (Z/Z_sun), used now and then, eg. Vink mass loss rate
   clogz = log10(max(czs/czsn, 1.0d-40))

   ! Set up coefficients for cubic spline interpolation in opacity
   call opspln_co

end subroutine load_opacity_co

!> -----------------------------------------------------------------------
!! Computes the opacity for a given composition by interpolating in
!! Temperature and density/R recently changed by SdM
!! -----------------------------------------------------------------------
!! INPUT variables - JX: index referring to a certain composition
!!                 - FT: log10 T
!!                 - FR: log10 R = log (T6/rho^3)
!!
!! OUTPUT Variables
!!   - FKL: log opacity for nearest lower composition value in table
!!   - FKH: log opacity for nearest higher composition value in table
!! The actual opacity is found by doing a linear interpolation between
!! FKL and FKH
!!
!! Uses:
!!  + MODULE OPACITY_CO
!!    - SPLINE_OPAC_CO: spline coefficients for CO enhanced opacity tables
!!    - opT(1:CO_MT): log10 T at which CO opacities are defined
!!    - opR(1:CO_MR): log10 R at which CO opacities are defined
!!  + Former common block STAT3
!!
!!  \todo FIXME: CB no longer exists, variable names changed
!!
!!    - F: spline coefficients for Pols opacity tables
!!    - FTM(127): log10 T at which Pols opacities are defined
!!    - FRM(90): log10 Rho at which Pols opacities are defined
!! -----------------------------------------------------------------------
!<
elemental subroutine opacty_co ( jjx, ft, fr, fkl, fkh )
   use real_kind

   implicit none
   integer, intent(in) :: jjx
   real(double), intent(in) :: ft, fr
   real(double), intent(out) :: fkl, fkh

   real(double) :: xft, xfr, dt, dr
   integer :: jx, next_jx, i1, i2

   ! Calculate a bicubic spline interpolation fit for the temperature
   ! and density opacity fit.  Do not stop if the input lies outside the
   ! array but rather use the value at the nearest edge-point.
   xft = max(min(ft, opt(co_mt)), opt(1))
   xfr = max(min(fr, opr(co_mr)), opr(1))
   ! Find interval in which target point lies. Index I1 is index of
   ! nearest lower(!) entry in opT/opR

   i1 = 1 + int((co_mt - 1)*(xft - opt(1))/(opt(co_mt) - opt(1)))
   i2 = 1 + int((co_mr - 1)*(xfr - opr(1))/(opr(co_mr) - opr(1)))

   dt = xft - opt(i1)
   dr = xfr - opr(i2)

   ! Evaluate the splines.
   next_jx = 61
   jx = jjx
   if (jx+next_jx>305) then
      jx = jx-next_jx
   end if

   fkl = spline_opac_co(1, 1, i1, i2, jx) + dr*(spline_opac_co(1, 2, i1, i2, jx)&
        + dr*(spline_opac_co(1, 3, i1, i2, jx) + dr* spline_opac_co(1, 4, i1, i2, jx)))&
        + dt*(spline_opac_co(2, 1, i1, i2, jx) + dr*(spline_opac_co(2, 2, i1, i2, jx)&
        + dr*(spline_opac_co(2, 3, i1, i2, jx) + dr* spline_opac_co(2, 4, i1, i2, jx)))&
        + dt*(spline_opac_co(3, 1, i1, i2, jx) + dr*(spline_opac_co(3, 2, i1, i2, jx)&
        + dr*(spline_opac_co(3, 3, i1, i2, jx) + dr* spline_opac_co(3, 4, i1, i2, jx)))&
        + dt*(spline_opac_co(4, 1, i1, i2, jx) + dr*(spline_opac_co(4, 2, i1, i2, jx)&
        + dr*(spline_opac_co(4, 3, i1, i2, jx) + dr* spline_opac_co(4, 4, i1, i2, jx))))))

   fkh = spline_opac_co(1, 1, i1, i2, jx + next_jx) + dr*(spline_opac_co(1, 2, i1, i2, jx + next_jx)&
        + dr*(spline_opac_co(1, 3, i1, i2, jx + next_jx) + dr* spline_opac_co(1, 4, i1, i2, jx + next_jx)))&
        + dt*(spline_opac_co(2, 1, i1, i2, jx + next_jx) + dr*(spline_opac_co(2, 2, i1, i2, jx + next_jx)&
        + dr*(spline_opac_co(2, 3, i1, i2, jx + next_jx) + dr* spline_opac_co(2, 4, i1, i2, jx + next_jx)))&
        + dt*(spline_opac_co(3, 1, i1, i2, jx + next_jx) + dr*(spline_opac_co(3, 2, i1, i2, jx + next_jx)&
        + dr*(spline_opac_co(3, 3, i1, i2, jx + next_jx) + dr* spline_opac_co(3, 4, i1, i2, jx + next_jx)))&
        + dt*(spline_opac_co(4, 1, i1, i2, jx + next_jx) + dr*(spline_opac_co(4, 2, i1, i2, jx + next_jx)&
        + dr*(spline_opac_co(4, 3, i1, i2, jx + next_jx) + dr* spline_opac_co(4, 4, i1, i2, jx + next_jx))))))

   return      ! We normally don't want any output here...
   !if ( (xft /= ft) .or. (xfr /= fr) ) write (10,100) ft, fr

100 format ('opacity out of range',' ft fr ', 2f9.4)
end subroutine opacty_co

end module opacity_co

max_stars_repo_name: rknop/amuse

max_issues_repo_name: rknop/amuse

max_forks_repo_name: rieder/amuse


text:
      program complexe02
      complex z

c     This is a complex constant
      print *, (5, 6)

c     This is diagnosed as an invalid implicit DO
c     print *, (5, 6, 7)

      print *, (5, 6, I = 1, 2)

      end


text:
!==============================================================================!
  subroutine User_Mod_Initialize(grid)
!------------------------------------------------------------------------------!
!   This is embrio of a function which will be used to initialize variables    !
!   by the user.                                                               !
!------------------------------------------------------------------------------!
!----------------------------------[Modules]-----------------------------------!
  use Grid_Mod
!------------------------------------------------------------------------------!
  implicit none
!---------------------------------[Arguments]----------------------------------!
  type(Grid_Type)    :: grid
!------------------------------------------------------------------------------!

  end subroutine

max_stars_repo_name: MassimoBorrelliPhysicist/T-Flows

max_issues_repo_name: MassimoBorrelliPhysicist/T-Flows

max_forks_repo_name: MassimoBorrelliPhysicist/T-Flows


text:
Subroutine pystfu(kf, x, q2, xpq, jbt)
  Common /hparnt/hipr1(100), ihpr2(50), hint1(100), ihnt2(50)
  Save /hparnt/
  Common /hjcrdn/yp(3, 300), yt(3, 300)
  Save /hjcrdn/
  Common /ludat1/mstu(200), paru(200), mstj(200), parj(200)
  Save /ludat1/
  Common /ludat2/kchg(500, 3), pmas(500, 4), parf(2000), vckm(4, 4)
  Save /ludat2/
  Common /pypars/mstp(200), parp(200), msti(200), pari(200)
  Save /pypars/
  Common /pyint1/mint(400), vint(400)
  Save /pyint1/
  Dimension xpq(-6:6), xq(6), tx(6), tt(6), ts(6), nehlq(8, 2), cehlq(6, 6, 2, 8, 2), cdo(3, 6, 5, 2), cow(3, 5, 4, 2)
  Data nehlq/3, 4, 7, 5, 7, 7, 7, 7, 3, 4, 7, 6, 7, 7, 7, 7/
  Data (((cehlq(ix,it,nx,1,1),ix=1,6),it=1,6), nx=1, 2)/7.677E-01, -2.087E-01, -3.303E-01, -2.517E-02, -1.570E-02, -1.000E-04, -5.326E-01, -2.661E-01, 3.201E-01, 1.192E-01, 2.434E-02, 7.620E-03, 2.162E-01, 1.881E-01, -8.375E-02, -6.515E-02, -1.743E-02, -5.040E-03, -9.211E-02, -9.952E-02, 1.373E-02, 2.506E-02, 8.770E-03, 2.550E-03, 3.670E-02, 4.409E-02, 9.600E-04, -7.960E-03, -3.420E-03, -1.050E-03, -1.549E-02, -2.026E-02, -3.060E-03, 2.220E-03, 1.240E-03, 4.100E-04, 2.395E-01, 2.905E-01, 9.778E-02, 2.149E-02, 3.440E-03, 5.000E-04, 1.751E-02, -6.090E-03, -2.687E-02, -1.916E-02, -7.970E-03, -2.750E-03, -5.760E-03, -5.040E-03, 1.080E-03, 2.490E-03, 1.530E-03, 7.500E-04, 1.740E-03, 1.960E-03, 3.000E-04, -3.400E-04, -2.900E-04, -1.800E-04, -5.300E-04, -6.400E-04, -1.700E-04, 4.000E-05, 6.000E-05, 4.000E-05, 1.700E-04, 2.200E-04, 8.000E-05, 1.000E-05, -1.000E-05, -1.000E-05/
  Data (((cehlq(ix,it,nx,1,2),ix=1,6),it=1,6), nx=1, 2)/7.237E-01, -2.189E-01, -2.995E-01, -1.909E-02, -1.477E-02, 2.500E-04, -5.314E-01, -2.425E-01, 3.283E-01, 1.119E-01, 2.223E-02, 7.070E-03, 2.289E-01, 1.890E-01, -9.859E-02, -6.900E-02, -1.747E-02, -5.080E-03, -1.041E-01, -1.084E-01, 2.108E-02, 2.975E-02, 9.830E-03, 2.830E-03, 4.394E-02, 5.116E-02, -1.410E-03, -1.055E-02, -4.230E-03, -1.270E-03, -1.991E-02, -2.539E-02, -2.780E-03, 3.430E-03, 1.720E-03, 5.500E-04, 2.410E-01, 2.884E-01, 9.369E-02, 1.900E-02, 2.530E-03, 2.400E-04, 1.765E-02, -9.220E-03, -3.037E-02, -2.085E-02, -8.440E-03, -2.810E-03, -6.450E-03, -5.260E-03, 1.720E-03, 3.110E-03, 1.830E-03, 8.700E-04, 2.120E-03, 2.320E-03, 2.600E-04, -4.900E-04, -3.900E-04, -2.300E-04, -6.900E-04, -8.200E-04, -2.000E-04, 7.000E-05, 9.000E-05, 6.000E-05, 2.400E-04, 3.100E-04, 1.100E-04, 0.000E+00, -2.000E-05, -2.000E-05/
  Data (((cehlq(ix,it,nx,2,1),ix=1,6),it=1,6), nx=1, 2)/3.813E-01, -8.090E-02, -1.634E-01, -2.185E-02, -8.430E-03, -6.200E-04, -2.948E-01, -1.435E-01, 1.665E-01, 6.638E-02, 1.473E-02, 4.080E-03, 1.252E-01, 1.042E-01, -4.722E-02, -3.683E-02, -1.038E-02, -2.860E-03, -5.478E-02, -5.678E-02, 8.900E-03, 1.484E-02, 5.340E-03, 1.520E-03, 2.220E-02, 2.567E-02, -3.000E-05, -4.970E-03, -2.160E-03, -6.500E-04, -9.530E-03, -1.204E-02, -1.510E-03, 1.510E-03, 8.300E-04, 2.700E-04, 1.261E-01, 1.354E-01, 3.958E-02, 8.240E-03, 1.660E-03, 4.500E-04, 3.890E-03, -1.159E-02, -1.625E-02, -9.610E-03, -3.710E-03, -1.260E-03, -1.910E-03, -5.600E-04, 1.590E-03, 1.590E-03, 8.400E-04, 3.900E-04, 6.400E-04, 4.900E-04, -1.500E-04, -2.900E-04, -1.800E-04, -1.000E-04, -2.000E-04, -1.900E-04, 0.000E+00, 6.000E-05, 4.000E-05, 3.000E-05, 7.000E-05, 8.000E-05, 2.000E-05, -1.000E-05, -1.000E-05, -1.000E-05/
  Data (((cehlq(ix,it,nx,2,2),ix=1,6),it=1,6), nx=1, 2)/3.578E-01, -8.622E-02, -1.480E-01, -1.840E-02, -7.820E-03, -4.500E-04, -2.925E-01, -1.304E-01, 1.696E-01, 6.243E-02, 1.353E-02, 3.750E-03, 1.318E-01, 1.041E-01, -5.486E-02, -3.872E-02, -1.038E-02, -2.850E-03, -6.162E-02, -6.143E-02, 1.303E-02, 1.740E-02, 5.940E-03, 1.670E-03, 2.643E-02, 2.957E-02, -1.490E-03, -6.450E-03, -2.630E-03, -7.700E-04, -1.218E-02, -1.497E-02, -1.260E-03, 2.240E-03, 1.120E-03, 3.500E-04, 1.263E-01, 1.334E-01, 3.732E-02, 7.070E-03, 1.260E-03, 3.400E-04, 3.660E-03, -1.357E-02, -1.795E-02, -1.031E-02, -3.880E-03, -1.280E-03, -2.100E-03, -3.600E-04, 2.050E-03, 1.920E-03, 9.800E-04, 4.400E-04, 7.700E-04, 5.400E-04, -2.400E-04, -3.900E-04, -2.400E-04, -1.300E-04, -2.600E-04, -2.300E-04, 2.000E-05, 9.000E-05, 6.000E-05, 4.000E-05, 9.000E-05, 1.000E-04, 2.000E-05, -2.000E-05, -2.000E-05, -1.000E-05/
  Data (((cehlq(ix,it,nx,3,1),ix=1,6),it=1,6), nx=1, 2)/6.870E-02, -6.861E-02, 2.973E-02, -5.400E-03, 3.780E-03, -9.700E-04, -1.802E-02, 1.400E-04, 6.490E-03, -8.540E-03, 1.220E-03, -1.750E-03, -4.650E-03, 1.480E-03, -5.930E-03, 6.000E-04, -1.030E-03, -8.000E-05, 6.440E-03, 2.570E-03, 2.830E-03, 1.150E-03, 7.100E-04, 3.300E-04, -3.930E-03, -2.540E-03, -1.160E-03, -7.700E-04, -3.600E-04, -1.900E-04, 2.340E-03, 1.930E-03, 5.300E-04, 3.700E-04, 1.600E-04, 9.000E-05, 1.014E+00, -1.106E+00, 3.374E-01, -7.444E-02, 8.850E-03, -8.700E-04, 9.233E-01, -1.285E+00, 4.475E-01, -9.786E-02, 1.419E-02, -1.120E-03, 4.888E-02, -1.271E-01, 8.606E-02, -2.608E-02, 4.780E-03, -6.000E-04, -2.691E-02, 4.887E-02, -1.771E-02, 1.620E-03, 2.500E-04, -6.000E-05, 7.040E-03, -1.113E-02, 1.590E-03, 7.000E-04, -2.000E-04, 0.000E+00, -1.710E-03, 2.290E-03, 3.800E-04, -3.500E-04, 4.000E-05, 1.000E-05/
  Data (((cehlq(ix,it,nx,3,2),ix=1,6),it=1,6), nx=1, 2)/1.008E-01, -7.100E-02, 1.973E-02, -5.710E-03, 2.930E-03, -9.900E-04, -5.271E-02, -1.823E-02, 1.792E-02, -6.580E-03, 1.750E-03, -1.550E-03, 1.220E-02, 1.763E-02, -8.690E-03, -8.800E-04, -1.160E-03, -2.100E-04, -1.190E-03, -7.180E-03, 2.360E-03, 1.890E-03, 7.700E-04, 4.100E-04, -9.100E-04, 2.040E-03, -3.100E-04, -1.050E-03, -4.000E-04, -2.400E-04, 1.190E-03, -1.700E-04, -2.000E-04, 4.200E-04, 1.700E-04, 1.000E-04, 1.081E+00, -1.189E+00, 3.868E-01, -8.617E-02, 1.115E-02, -1.180E-03, 9.917E-01, -1.396E+00, 4.998E-01, -1.159E-01, 1.674E-02, -1.720E-03, 5.099E-02, -1.338E-01, 9.173E-02, -2.885E-02, 5.890E-03, -6.500E-04, -3.178E-02, 5.703E-02, -2.070E-02, 2.440E-03, 1.100E-04, -9.000E-05, 8.970E-03, -1.392E-02, 2.050E-03, 6.500E-04, -2.300E-04, 2.000E-05, -2.340E-03, 3.010E-03, 5.000E-04, -3.900E-04, 6.000E-05, 1.000E-05/
  Data (((cehlq(ix,it,nx,4,1),ix=1,6),it=1,6), nx=1, 2)/9.482E-01, -9.578E-01, 1.009E-01, -1.051E-01, 3.456E-02, -3.054E-02, -9.627E-01, 5.379E-01, 3.368E-01, -9.525E-02, 1.488E-02, -2.051E-02, 4.300E-01, -8.306E-02, -3.372E-01, 4.902E-02, -9.160E-03, 1.041E-02, -1.925E-01, -1.790E-02, 2.183E-01, 7.490E-03, 4.140E-03, -1.860E-03, 8.183E-02, 1.926E-02, -1.072E-01, -1.944E-02, -2.770E-03, -5.200E-04, -3.884E-02, -1.234E-02, 5.410E-02, 1.879E-02, 3.350E-03, 1.040E-03, 2.948E+01, -3.902E+01, 1.464E+01, -3.335E+00, 5.054E-01, -5.915E-02, 2.559E+01, -3.955E+01, 1.661E+01, -4.299E+00, 6.904E-01, -8.243E-02, -1.663E+00, 1.176E+00, 1.118E+00, -7.099E-01, 1.948E-01, -2.404E-02, -2.168E-01, 8.170E-01, -7.169E-01, 1.851E-01, -1.924E-02, -3.250E-03, 2.088E-01, -4.355E-01, 2.239E-01, -2.446E-02, -3.620E-03, 1.910E-03, -9.097E-02, 1.601E-01, -5.681E-02, -2.500E-03, 2.580E-03, -4.700E-04/
  Data (((cehlq(ix,it,nx,4,2),ix=1,6),it=1,6), nx=1, 2)/2.367E+00, 4.453E-01, 3.660E-01, 9.467E-02, 1.341E-01, 1.661E-02, -3.170E+00, -1.795E+00, 3.313E-02, -2.874E-01, -9.827E-02, -7.119E-02, 1.823E+00, 1.457E+00, -2.465E-01, 3.739E-02, 6.090E-03, 1.814E-02, -1.033E+00, -9.827E-01, 2.136E-01, 1.169E-01, 5.001E-02, 1.684E-02, 5.133E-01, 5.259E-01, -1.173E-01, -1.139E-01, -4.988E-02, -2.021E-02, -2.881E-01, -3.145E-01, 5.667E-02, 9.161E-02, 4.568E-02, 1.951E-02, 3.036E+01, -4.062E+01, 1.578E+01, -3.699E+00, 6.020E-01, -7.031E-02, 2.700E+01, -4.167E+01, 1.770E+01, -4.804E+00, 7.862E-01, -1.060E-01, -1.909E+00, 1.357E+00, 1.127E+00, -7.181E-01, 2.232E-01, -2.481E-02, -2.488E-01, 9.781E-01, -8.127E-01, 2.094E-01, -2.997E-02, -4.710E-03, 2.506E-01, -5.427E-01, 2.672E-01, -3.103E-02, -1.800E-03, 2.870E-03, -1.128E-01, 2.087E-01, -6.972E-02, -2.480E-03, 2.630E-03, -8.400E-04/
  Data (((cehlq(ix,it,nx,5,1),ix=1,6),it=1,6), nx=1, 2)/4.968E-02, -4.173E-02, 2.102E-02, -3.270E-03, 3.240E-03, -6.700E-04, -6.150E-03, -1.294E-02, 6.740E-03, -6.890E-03, 9.000E-04, -1.510E-03, -8.580E-03, 5.050E-03, -4.900E-03, -1.600E-04, -9.400E-04, -1.500E-04, 7.840E-03, 1.510E-03, 2.220E-03, 1.400E-03, 7.000E-04, 3.500E-04, -4.410E-03, -2.220E-03, -8.900E-04, -8.500E-04, -3.600E-04, -2.000E-04, 2.520E-03, 1.840E-03, 4.100E-04, 3.900E-04, 1.600E-04, 9.000E-05, 9.235E-01, -1.085E+00, 3.464E-01, -7.210E-02, 9.140E-03, -9.100E-04, 9.315E-01, -1.274E+00, 4.512E-01, -9.775E-02, 1.380E-02, -1.310E-03, 4.739E-02, -1.296E-01, 8.482E-02, -2.642E-02, 4.760E-03, -5.700E-04, -2.653E-02, 4.953E-02, -1.735E-02, 1.750E-03, 2.800E-04, -6.000E-05, 6.940E-03, -1.132E-02, 1.480E-03, 6.500E-04, -2.100E-04, 0.000E+00, -1.680E-03, 2.340E-03, 4.200E-04, -3.400E-04, 5.000E-05, 1.000E-05/
  Data (((cehlq(ix,it,nx,5,2),ix=1,6),it=1,6), nx=1, 2)/6.478E-02, -4.537E-02, 1.643E-02, -3.490E-03, 2.710E-03, -6.700E-04, -2.223E-02, -2.126E-02, 1.247E-02, -6.290E-03, 1.120E-03, -1.440E-03, -1.340E-03, 1.362E-02, -6.130E-03, -7.900E-04, -9.000E-04, -2.000E-04, 5.080E-03, -3.610E-03, 1.700E-03, 1.830E-03, 6.800E-04, 4.000E-04, -3.580E-03, 6.000E-05, -2.600E-04, -1.050E-03, -3.800E-04, -2.300E-04, 2.420E-03, 9.300E-04, -1.000E-04, 4.500E-04, 1.700E-04, 1.100E-04, 9.868E-01, -1.171E+00, 3.940E-01, -8.459E-02, 1.124E-02, -1.250E-03, 1.001E+00, -1.383E+00, 5.044E-01, -1.152E-01, 1.658E-02, -1.830E-03, 4.928E-02, -1.368E-01, 9.021E-02, -2.935E-02, 5.800E-03, -6.600E-04, -3.133E-02, 5.785E-02, -2.023E-02, 2.630E-03, 1.600E-04, -8.000E-05, 8.840E-03, -1.416E-02, 1.900E-03, 5.800E-04, -2.500E-04, 1.000E-05, -2.300E-03, 3.080E-03, 5.500E-04, -3.700E-04, 7.000E-05, 1.000E-05/
  Data (((cehlq(ix,it,nx,6,1),ix=1,6),it=1,6), nx=1, 2)/9.270E-03, -1.817E-02, 9.590E-03, -6.390E-03, 1.690E-03, -1.540E-03, 5.710E-03, -1.188E-02, 6.090E-03, -4.650E-03, 1.240E-03, -1.310E-03, -3.960E-03, 7.100E-03, -3.590E-03, 1.840E-03, -3.900E-04, 3.400E-04, 1.120E-03, -1.960E-03, 1.120E-03, -4.800E-04, 1.000E-04, -4.000E-05, 4.000E-05, -3.000E-05, -1.800E-04, 9.000E-05, -5.000E-05, -2.000E-05, -4.200E-04, 7.300E-04, -1.600E-04, 5.000E-05, 5.000E-05, 5.000E-05, 8.098E-01, -1.042E+00, 3.398E-01, -6.824E-02, 8.760E-03, -9.000E-04, 8.961E-01, -1.217E+00, 4.339E-01, -9.287E-02, 1.304E-02, -1.290E-03, 3.058E-02, -1.040E-01, 7.604E-02, -2.415E-02, 4.600E-03, -5.000E-04, -2.451E-02, 4.432E-02, -1.651E-02, 1.430E-03, 1.200E-04, -1.000E-04, 1.122E-02, -1.457E-02, 2.680E-03, 5.800E-04, -1.200E-04, 3.000E-05, -7.730E-03, 7.330E-03, -7.600E-04, -2.400E-04, 1.000E-05, 0.000E+00/
  Data (((cehlq(ix,it,nx,6,2),ix=1,6),it=1,6), nx=1, 2)/9.980E-03, -1.945E-02, 1.055E-02, -6.870E-03, 1.860E-03, -1.560E-03, 5.700E-03, -1.203E-02, 6.250E-03, -4.860E-03, 1.310E-03, -1.370E-03, -4.490E-03, 7.990E-03, -4.170E-03, 2.050E-03, -4.400E-04, 3.300E-04, 1.470E-03, -2.480E-03, 1.460E-03, -5.700E-04, 1.200E-04, -1.000E-05, -9.000E-05, 1.500E-04, -3.200E-04, 1.200E-04, -6.000E-05, -4.000E-05, -4.200E-04, 7.600E-04, -1.400E-04, 4.000E-05, 7.000E-05, 5.000E-05, 8.698E-01, -1.131E+00, 3.836E-01, -8.111E-02, 1.048E-02, -1.300E-03, 9.626E-01, -1.321E+00, 4.854E-01, -1.091E-01, 1.583E-02, -1.700E-03, 3.057E-02, -1.088E-01, 8.022E-02, -2.676E-02, 5.590E-03, -5.600E-04, -2.845E-02, 5.164E-02, -1.918E-02, 2.210E-03, -4.000E-05, -1.500E-04, 1.311E-02, -1.751E-02, 3.310E-03, 5.100E-04, -1.200E-04, 5.000E-05, -8.590E-03, 8.380E-03, -9.200E-04, -2.600E-04, 1.000E-05, -1.000E-05/
  Data (((cehlq(ix,it,nx,7,1),ix=1,6),it=1,6), nx=1, 2)/9.010E-03, -1.401E-02, 7.150E-03, -4.130E-03, 1.260E-03, -1.040E-03, 6.280E-03, -9.320E-03, 4.780E-03, -2.890E-03, 9.100E-04, -8.200E-04, -2.930E-03, 4.090E-03, -1.890E-03, 7.600E-04, -2.300E-04, 1.400E-04, 3.900E-04, -1.200E-03, 4.400E-04, -2.500E-04, 2.000E-05, -2.000E-05, 2.600E-04, 1.400E-04, -8.000E-05, 1.000E-04, 1.000E-05, 1.000E-05, -2.600E-04, 3.200E-04, 1.000E-05, -1.000E-05, 1.000E-05, -1.000E-05, 8.029E-01, -1.075E+00, 3.792E-01, -7.843E-02, 1.007E-02, -1.090E-03, 7.903E-01, -1.099E+00, 4.153E-01, -9.301E-02, 1.317E-02, -1.410E-03, -1.704E-02, -1.130E-02, 2.882E-02, -1.341E-02, 3.040E-03, -3.600E-04, -7.200E-04, 7.230E-03, -5.160E-03, 1.080E-03, -5.000E-05, -4.000E-05, 3.050E-03, -4.610E-03, 1.660E-03, -1.300E-04, -1.000E-05, 1.000E-05, -4.360E-03, 5.230E-03, -1.610E-03, 2.000E-04, -2.000E-05, 0.000E+00/
  Data (((cehlq(ix,it,nx,7,2),ix=1,6),it=1,6), nx=1, 2)/8.980E-03, -1.459E-02, 7.510E-03, -4.410E-03, 1.310E-03, -1.070E-03, 5.970E-03, -9.440E-03, 4.800E-03, -3.020E-03, 9.100E-04, -8.500E-04, -3.050E-03, 4.440E-03, -2.100E-03, 8.500E-04, -2.400E-04, 1.400E-04, 5.300E-04, -1.300E-03, 5.600E-04, -2.700E-04, 3.000E-05, -2.000E-05, 2.000E-04, 1.400E-04, -1.100E-04, 1.000E-04, 0.000E+00, 0.000E+00, -2.600E-04, 3.200E-04, 0.000E+00, -3.000E-05, 1.000E-05, -1.000E-05, 8.672E-01, -1.174E+00, 4.265E-01, -9.252E-02, 1.244E-02, -1.460E-03, 8.500E-01, -1.194E+00, 4.630E-01, -1.083E-01, 1.614E-02, -1.830E-03, -2.241E-02, -5.630E-03, 2.815E-02, -1.425E-02, 3.520E-03, -4.300E-04, -7.300E-04, 8.030E-03, -5.780E-03, 1.380E-03, -1.300E-04, -4.000E-05, 3.460E-03, -5.380E-03, 1.960E-03, -2.100E-04, 1.000E-05, 1.000E-05, -4.850E-03, 5.950E-03, -1.890E-03, 2.600E-04, -3.000E-05, 0.000E+00/
  Data (((cehlq(ix,it,nx,8,1),ix=1,6),it=1,6), nx=1, 2)/4.410E-03, -7.480E-03, 3.770E-03, -2.580E-03, 7.300E-04, -7.100E-04, 3.840E-03, -6.050E-03, 3.030E-03, -2.030E-03, 5.800E-04, -5.900E-04, -8.800E-04, 1.660E-03, -7.500E-04, 4.700E-04, -1.000E-04, 1.000E-04, -8.000E-05, -1.500E-04, 1.200E-04, -9.000E-05, 3.000E-05, 0.000E+00, 1.300E-04, -2.200E-04, -2.000E-05, -2.000E-05, -2.000E-05, -2.000E-05, -7.000E-05, 1.900E-04, -4.000E-05, 2.000E-05, 0.000E+00, 0.000E+00, 6.623E-01, -9.248E-01, 3.519E-01, -7.930E-02, 1.110E-02, -1.180E-03, 6.380E-01, -9.062E-01, 3.582E-01, -8.479E-02, 1.265E-02, -1.390E-03, -2.581E-02, 2.125E-02, 4.190E-03, -4.980E-03, 1.490E-03, -2.100E-04, 7.100E-04, 5.300E-04, -1.270E-03, 3.900E-04, -5.000E-05, -1.000E-05, 3.850E-03, -5.060E-03, 1.860E-03, -3.500E-04, 4.000E-05, 0.000E+00, -3.530E-03, 4.460E-03, -1.500E-03, 2.700E-04, -3.000E-05, 0.000E+00/
  Data (((cehlq(ix,it,nx,8,2),ix=1,6),it=1,6), nx=1, 2)/4.260E-03, -7.530E-03, 3.830E-03, -2.680E-03, 7.600E-04, -7.300E-04, 3.640E-03, -6.050E-03, 3.030E-03, -2.090E-03, 5.900E-04, -6.000E-04, -9.200E-04, 1.710E-03, -8.200E-04, 5.000E-04, -1.200E-04, 1.000E-04, -5.000E-05, -1.600E-04, 1.300E-04, -9.000E-05, 3.000E-05, 0.000E+00, 1.300E-04, -2.100E-04, -1.000E-05, -2.000E-05, -2.000E-05, -1.000E-05, -8.000E-05, 1.800E-04, -5.000E-05, 2.000E-05, 0.000E+00, 0.000E+00, 7.146E-01, -1.007E+00, 3.932E-01, -9.246E-02, 1.366E-02, -1.540E-03, 6.856E-01, -9.828E-01, 3.977E-01, -9.795E-02, 1.540E-02, -1.790E-03, -3.053E-02, 2.758E-02, 2.150E-03, -4.880E-03, 1.640E-03, -2.500E-04, 9.200E-04, 4.200E-04, -1.340E-03, 4.600E-04, -8.000E-05, -1.000E-05, 4.230E-03, -5.660E-03, 2.140E-03, -4.300E-04, 6.000E-05, 0.000E+00, -3.890E-03, 5.000E-03, -1.740E-03, 3.300E-04, -4.000E-05, 0.000E+00/
  Data ((cdo(ip,is,1,1),is=1,6), ip=1, 3)/4.190E-01, 3.460E+00, 4.400E+00, 0.000E+00, 0.000E+00, 0.000E+00, 4.000E-03, 7.240E-01, -4.860E+00, 0.000E+00, 0.000E+00, 0.000E+00, -7.000E-03, -6.600E-02, 1.330E+00, 0.000E+00, 0.000E+00, 0.000E+00/
  Data ((cdo(ip,is,1,2),is=1,6), ip=1, 3)/3.740E-01, 3.330E+00, 6.030E+00, 0.000E+00, 0.000E+00, 0.000E+00, 1.400E-02, 7.530E-01, -6.220E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, -7.600E-02, 1.560E+00, 0.000E+00, 0.000E+00, 0.000E+00/
  Data ((cdo(ip,is,2,1),is=1,6), ip=1, 3)/7.630E-01, 4.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, -2.370E-01, 6.270E-01, -4.210E-01, 0.000E+00, 0.000E+00, 0.000E+00, 2.600E-02, -1.900E-02, 3.300E-02, 0.000E+00, 0.000E+00, 0.000E+00/
  Data ((cdo(ip,is,2,2),is=1,6), ip=1, 3)/7.610E-01, 3.830E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, -2.320E-01, 6.270E-01, -4.180E-01, 0.000E+00, 0.000E+00, 0.000E+00, 2.300E-02, -1.900E-02, 3.600E-02, 0.000E+00, 0.000E+00, 0.000E+00/
  Data ((cdo(ip,is,3,1),is=1,6), ip=1, 3)/1.265E+00, 0.000E+00, 8.050E+00, 0.000E+00, 0.000E+00, 0.000E+00, -1.132E+00, -3.720E-01, 1.590E+00, 6.310E+00, -1.050E+01, 1.470E+01, 2.930E-01, -2.900E-02, -1.530E-01, -2.730E-01, -3.170E+00, 9.800E+00/
  Data ((cdo(ip,is,3,2),is=1,6), ip=1, 3)/1.670E+00, 0.000E+00, 9.150E+00, 0.000E+00, 0.000E+00, 0.000E+00, -1.920E+00, -2.730E-01, 5.300E-01, 1.570E+01, -1.010E+02, 2.230E+02, 5.820E-01, -1.640E-01, -7.630E-01, -2.830E+00, 4.470E+01, -1.170E+02/
  Data ((cdo(ip,is,4,1),is=1,6), ip=1, 3)/0.000E+00, -3.600E-02, 6.350E+00, 0.000E+00, 0.000E+00, 0.000E+00, 1.350E-01, -2.220E-01, 3.260E+00, -3.030E+00, 1.740E+01, -1.790E+01, -7.500E-02, -5.800E-02, -9.090E-01, 1.500E+00, -1.130E+01, 1.560E+01/
  Data ((cdo(ip,is,4,2),is=1,6), ip=1, 3)/0.000E+00, -1.200E-01, 3.510E+00, 0.000E+00, 0.000E+00, 0.000E+00, 6.700E-02, -2.330E-01, 3.660E+00, -4.740E-01, 9.500E+00, -1.660E+01, -3.100E-02, -2.300E-02, -4.530E-01, 3.580E-01, -5.430E+00, 1.550E+01/
  Data ((cdo(ip,is,5,1),is=1,6), ip=1, 3)/1.560E+00, 0.000E+00, 6.000E+00, 9.000E+00, 0.000E+00, 0.000E+00, -1.710E+00, -9.490E-01, 1.440E+00, -7.190E+00, -1.650E+01, 1.530E+01, 6.380E-01, 3.250E-01, -1.050E+00, 2.550E-01, 1.090E+01, -1.010E+01/
  Data ((cdo(ip,is,5,2),is=1,6), ip=1, 3)/8.790E-01, 0.000E+00, 4.000E+00, 9.000E+00, 0.000E+00, 0.000E+00, -9.710E-01, -1.160E+00, 1.230E+00, -5.640E+00, -7.540E+00, -5.960E-01, 4.340E-01, 4.760E-01, -2.540E-01, -8.170E-01, 5.500E+00, 1.260E-01/
  Data ((cow(ip,is,1,1),is=1,5), ip=1, 3)/4.0000E-01, 7.0000E-01, 0.0000E+00, 0.0000E+00, 0.0000E+00, -6.2120E-02, 6.4780E-01, 0.0000E+00, 0.0000E+00, 0.0000E+00, -7.1090E-03, 1.3350E-02, 0.0000E+00, 0.0000E+00, 0.0000E+00/
  Data ((cow(ip,is,1,2),is=1,5), ip=1, 3)/4.0000E-01, 6.2800E-01, 0.0000E+00, 0.0000E+00, 0.0000E+00, -5.9090E-02, 6.4360E-01, 0.0000E+00, 0.0000E+00, 0.0000E+00, -6.5240E-03, 1.4510E-02, 0.0000E+00, 0.0000E+00, 0.0000E+00/
  Data ((cow(ip,is,2,1),is=1,5), ip=1, 3)/8.8800E-01, 0.0000E+00, 3.1100E+00, 6.0000E+00, 0.0000E+00, -1.8020E+00, -1.5760E+00, -1.3170E-01, 2.8010E+00, -1.7280E+01, 1.8120E+00, 1.2000E+00, 5.0680E-01, -1.2160E+01, 2.0490E+01/
  Data ((cow(ip,is,2,2),is=1,5), ip=1, 3)/7.9400E-01, 0.0000E+00, 2.8900E+00, 6.0000E+00, 0.0000E+00, -9.1440E-01, -1.2370E+00, 5.9660E-01, -3.6710E+00, -8.1910E+00, 5.9660E-01, 6.5820E-01, -2.5500E-01, -2.3040E+00, 7.7580E+00/
  Data ((cow(ip,is,3,1),is=1,5), ip=1, 3)/9.0000E-01, 0.0000E+00, 5.0000E+00, 0.0000E+00, 0.0000E+00, -2.4280E-01, -2.1200E-01, 8.6730E-01, 1.2660E+00, 2.3820E+00, 1.3860E-01, 3.6710E-03, 4.7470E-02, -2.2150E+00, 3.4820E-01/
  Data ((cow(ip,is,3,2),is=1,5), ip=1, 3)/9.0000E-01, 0.0000E+00, 5.0000E+00, 0.0000E+00, 0.0000E+00, -1.4170E-01, -1.6970E-01, -2.4740E+00, -2.5340E+00, 5.6210E-01, -1.7400E-01, -9.6230E-02, 1.5750E+00, 1.3780E+00, -2.7010E-01/
  Data ((cow(ip,is,4,1),is=1,5), ip=1, 3)/0.0000E+00, -2.2120E-02, 2.8940E+00, 0.0000E+00, 0.0000E+00, 7.9280E-02, -3.7850E-01, 9.4330E+00, 5.2480E+00, 8.3880E+00, -6.1340E-02, -1.0880E-01, -1.0852E+01, -7.1870E+00, -1.1610E+01/
  Data ((cow(ip,is,4,2),is=1,5), ip=1, 3)/0.0000E+00, -8.8200E-02, 1.9240E+00, 0.0000E+00, 0.0000E+00, 6.2290E-02, -2.8920E-01, 2.4240E-01, -4.4630E+00, -8.3670E-01, -4.0990E-02, -1.0820E-01, 2.0360E+00, 5.2090E+00, -4.8400E-02/
  alam = 0.
  Do kfl = -6, 6
    xpq(kfl) = 0.
  End Do
  If (x<0. .Or. x>1.) Then
    Write (mstu(11), 1000) x
    Return
  End If
  kfa = iabs(kf)
  If (kfa/=211 .And. kfa/=2212 .And. kfa/=2112) Then
    Write (mstu(11), 1100) kf
    Return
  End If
  If (mstp(51)==0 .Or. mstp(52)>=2) Then
    kfe = kfa
    If (kfa==2112) kfe = 2212
    Call pystfe(kfe, x, q2, xpq)
    Goto 230
  End If
  If (kfa==211) Goto 200
  If (mstp(51)==1 .Or. mstp(51)==2) Then
    nset = mstp(51)
    If (nset==1) alam = 0.2
    If (nset==2) alam = 0.29
    tmin = log(5./alam**2)
    tmax = log(1E8/alam**2)
    If (mstp(52)==0) Then
      t = tmin
    Else
      t = log(q2/alam**2)
    End If
    vt = max(-1., min(1.,(2.*t-tmax-tmin)/(tmax-tmin)))
    nx = 1
    If (x<=0.1) nx = 2
    If (nx==1) vx = (2.*x-1.1)/0.9
    If (nx==2) vx = max(-1., (2.*log(x)+11.51293)/6.90776)
    cxs = 1.
    If (x<1E-4 .And. abs(parp(51)-1.)>0.01) cxs = (1E-4/x)**(parp(51)-1.)
    tx(1) = 1.
    tx(2) = vx
    tx(3) = 2.*vx**2 - 1.
    tx(4) = 4.*vx**3 - 3.*vx
    tx(5) = 8.*vx**4 - 8.*vx**2 + 1.
    tx(6) = 16.*vx**5 - 20.*vx**3 + 5.*vx
    tt(1) = 1.
    tt(2) = vt
    tt(3) = 2.*vt**2 - 1.
    tt(4) = 4.*vt**3 - 3.*vt
    tt(5) = 8.*vt**4 - 8.*vt**2 + 1.
    tt(6) = 16.*vt**5 - 20.*vt**3 + 5.*vt
    Do kfl = 1, 6
      xqsum = 0.
      Do it = 1, 6
        Do ix = 1, 6
          xqsum = xqsum + cehlq(ix, it, nx, kfl, nset)*tx(ix)*tt(it)
        End Do
      End Do
      xq(kfl) = xqsum*(1.-x)**nehlq(kfl, nset)*cxs
    End Do
    xpq(0) = xq(4)
    xpq(1) = xq(2) + xq(3)
    xpq(2) = xq(1) + xq(3)
    xpq(3) = xq(5)
    xpq(4) = xq(6)
    xpq(-1) = xq(3)
    xpq(-2) = xq(3)
    xpq(-3) = xq(5)
    xpq(-4) = xq(6)
    If (mstp(54)>=5) Then
      If (nset==1) tmin = 8.1905
      If (nset==2) tmin = 7.4474
      If (t<=tmin) Goto 140
      vt = max(-1., min(1.,(2.*t-tmax-tmin)/(tmax-tmin)))
      tt(1) = 1.
      tt(2) = vt
      tt(3) = 2.*vt**2 - 1.
      tt(4) = 4.*vt**3 - 3.*vt
      tt(5) = 8.*vt**4 - 8.*vt**2 + 1.
      tt(6) = 16.*vt**5 - 20.*vt**3 + 5.*vt
      xqsum = 0.
      Do it = 1, 6
        Do ix = 1, 6
          xqsum = xqsum + cehlq(ix, it, nx, 7, nset)*tx(ix)*tt(it)
        End Do
      End Do
      xpq(5) = xqsum*(1.-x)**nehlq(7, nset)
      xpq(-5) = xpq(5)
      140 Continue
    End If
    If (mstp(54)>=6) Then
      If (nset==1) tmin = 11.5528
      If (nset==2) tmin = 10.8097
      tmin = tmin + 2.*log(pmas(6,1)/30.)
      tmax = tmax + 2.*log(pmas(6,1)/30.)
      If (t<=tmin) Goto 160
      vt = max(-1., min(1.,(2.*t-tmax-tmin)/(tmax-tmin)))
      tt(1) = 1.
      tt(2) = vt
      tt(3) = 2.*vt**2 - 1.
      tt(4) = 4.*vt**3 - 3.*vt
      tt(5) = 8.*vt**4 - 8.*vt**2 + 1.
      tt(6) = 16.*vt**5 - 20.*vt**3 + 5.*vt
      xqsum = 0.
      Do it = 1, 6
        Do ix = 1, 6
          xqsum = xqsum + cehlq(ix, it, nx, 8, nset)*tx(ix)*tt(it)
        End Do
      End Do
      xpq(6) = xqsum*(1.-x)**nehlq(8, nset)
      xpq(-6) = xpq(6)
      160 Continue
    End If
  Else If (mstp(51)==3 .Or. mstp(51)==4) Then
    nset = mstp(51) - 2
    If (nset==1) alam = 0.2
    If (nset==2) alam = 0.4
    If (mstp(52)<=0) Then
      sd = 0.
    Else
      sd = log(log(max(q2,4.)/alam**2)/log(4./alam**2))
    End If
    Do kfl = 1, 5
      Do is = 1, 6
        ts(is) = cdo(1, is, kfl, nset) + cdo(2, is, kfl, nset)*sd + cdo(3, is, kfl, nset)*sd**2
      End Do
      If (kfl<=2) Then
        eulbt1 = pygamm(ts(1))*pygamm(ts(2)+1.)/pygamm(ts(1)+ts(2)+1.)
        xq(kfl) = x**ts(1)*(1.-x)**ts(2)*(1.+ts(3)*x)/(eulbt1*(1.+ts(3)*ts(1)/(ts(1)+ts(2)+1.)))
      Else
        xq(kfl) = ts(1)*x**ts(2)*(1.-x)**ts(3)*(1.+ts(4)*x+ts(5)*x**2+ts(6)*x**3)
      End If
    End Do
    xpq(0) = xq(5)
    xpq(1) = xq(2) + xq(3)/6.
    xpq(2) = 3.*xq(1) - xq(2) + xq(3)/6.
    xpq(3) = xq(3)/6.
    xpq(4) = xq(4)
    xpq(-1) = xq(3)/6.
    xpq(-2) = xq(3)/6.
    xpq(-3) = xq(3)/6.
    xpq(-4) = xq(4)
  Else If (mstp(51)>=11 .And. mstp(51)<=13) Then
    Call pystfe(2212, x, q2, xpq)
  Else
    Write (mstu(11), 1200) mstp(51)
  End If
  Goto 230
  200 If ((mstp(51)>=1 .And. mstp(51)<=4) .Or. (mstp(51)>=11 .And. mstp(51)<=13)) Then
    nset = 1
    If (mstp(51)==2 .Or. mstp(51)==4 .Or. mstp(51)==13) nset = 2
    If (nset==1) alam = 0.2
    If (nset==2) alam = 0.4
    If (mstp(52)<=0) Then
      sd = 0.
    Else
      sd = log(log(max(q2,4.)/alam**2)/log(4./alam**2))
    End If
    Do kfl = 1, 4
      Do is = 1, 5
        ts(is) = cow(1, is, kfl, nset) + cow(2, is, kfl, nset)*sd + cow(3, is, kfl, nset)*sd**2
      End Do
      If (kfl==1) Then
        eulbt2 = pygamm(ts(1))*pygamm(ts(2)+1.)/pygamm(ts(1)+ts(2)+1.)
        xq(kfl) = x**ts(1)*(1.-x)**ts(2)/eulbt2
      Else
        xq(kfl) = ts(1)*x**ts(2)*(1.-x)**ts(3)*(1.+ts(4)*x+ts(5)*x**2)
      End If
    End Do
    xpq(0) = xq(2)
    xpq(1) = xq(3)/6.
    xpq(2) = xq(1) + xq(3)/6.
    xpq(3) = xq(3)/6.
    xpq(4) = xq(4)
    xpq(-1) = xq(1) + xq(3)/6.
    xpq(-2) = xq(3)/6.
    xpq(-3) = xq(3)/6.
    xpq(-4) = xq(4)
  Else
    Write (mstu(11), 1200) mstp(51)
  End If
  230 If (kfa==2112) Then
    xps = xpq(1)
    xpq(1) = xpq(2)
    xpq(2) = xps
    xps = xpq(-1)
    xpq(-1) = xpq(-2)
    xpq(-2) = xps
  End If
  If (kf<0) Then
    Do kfl = 1, 4
      xps = xpq(kfl)
      xpq(kfl) = xpq(-kfl)
      xpq(-kfl) = xps
    End Do
  End If
  Do kfl = -6, 6
    xpq(kfl) = max(0., xpq(kfl))
    If (iabs(kfl)>mstp(54)) xpq(kfl) = 0.
  End Do
  If ((jbt/=1 .And. jbt/=2) .Or. ihpr2(6)==0 .Or. ihnt2(16)==1) Goto 400
  atnm = ihnt2(2*jbt-1)
  If (atnm<=1.0) Goto 400
  If (jbt==1) Then
    bbr2 = (yp(1,ihnt2(11))**2+yp(2,ihnt2(11))**2)/1.44/atnm**0.66666
  Else If (jbt==2) Then
    bbr2 = (yt(1,ihnt2(12))**2+yt(2,ihnt2(12))**2)/1.44/atnm**0.66666
  End If
  bbr2 = min(1.0, bbr2)
  abx = (atnm**0.33333333-1.0)
  apx = hipr1(6)*4.0/3.0*abx*sqrt(1.0-bbr2)
  aax = 1.192*alog(atnm)**0.1666666
  rrx = aax*(x**3-1.2*x**2+0.21*x) + 1.0 - (apx-1.079*abx*sqrt(x)/alog(atnm+1.0))*exp(-x**2.0/0.01)
  Do kfl = -6, 6
    xpq(kfl) = xpq(kfl)*rrx
  End Do
  400 Continue
  Return
  1000 Format (' Error: x value outside physical range, x =', 1P, E12.3)
  1100 Format (' Error: illegal particle code for structure function,', ' KF =', I5)
  1200 Format (' Error: bad value of parameter MSTP(51) in PYSTFU,', ' MSTP(51) =', I5)
End Subroutine pystfu

max_stars_repo_name: xiaohaijin/AMPT

max_issues_repo_name: xiaohaijin/AMPT

max_forks_repo_name: xiaohaijin/AMPT


text:
Davis water wells are the source of the tap water in Davis. These are scattered about town and generally the water from the faucet comes from the nearest well(s). Also note that the water towers are used in the distribution scheme to maintain pressure and distribution across the system.

Currently, all the water supplied by the city comes from an aquifer.  In general, the city sits over an aquifer system with two distinct layers.  The upper layer, known as the intermediate depth aquifer, extends approximately 700 feet below the ground surface.  A second layer is known as the deep depth aquifer, is separated from the intermediate depth aquifer by a less permeable clay layer.  Originally, the city supplied all its water from the intermediate depth aquifer.  In recent years however, in order to achieve greater water quality, the city has begun shifting its water supply over to the deep depth aquifer.  Older wells in Davis typically extend 300 to 600 feet below the ground surface, while newer deep wells in Davis may extend 1500 to 1800 feet below the ground surface.  A depth chart of all the well in davis can be seen http://cityofdavis.org/pw/water/pdfs/Well_Profiles_Color.pdf here.

More information on the citys water system can be found in its http://cityofdavis.org/pw/water/pdfs/2005UrbanWaterManagementPlanUpdate.pdf urban water management plan.

Currently Operational Wells
  Davis Well 1
  Davis Well 7
  Davis Well 11
  Davis Well 12
  Davis Well 13
  Davis Well 14
  Davis Well 15  May be abandoned
  Davis Well 18
  Davis Well 19
  Davis Well 20
  Davis Well 21
  Davis Well 22
  Davis Well 23
  Davis Well 24
  Davis Well 25
  Davis Well 26
  Davis Well 27
  Davis Deep Well 28
  Davis Well EM2
  Davis Well EM3
  Davis Deep Well 29
  Davis Deep Well 30
  Davis Deep Well 31
Abandoned Wells
  Davis Well 16

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
HANSARD REVISE * NUMERO 121a
Le jeudi 18 juin 1998
HANSARD REVISE * NUMERO 121a
Le jeudi 18 juin 1998
La seance est ouverte a 21 h 50.
L'huissier du baton noir apporte le message suivant:
En consequence, le President et les deputes se rendent au Senat.
(La seance est levee a 22 h 15.)


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

C=======================================================================
      SUBROUTINE PLTCUR(X,Y,NUM)
      REAL DEVCAP(23)
      REAL DEFOUT(7)
      COMMON /STATUS/DEVCAP,DEFOUT
      REAL DEVP(5)
      COMMON /DEVICE/DEVP
      REAL COLP(3)
      REAL PALETT(3,16)
      COMMON /COLOR/COLP,PALETT
      REAL TEXTP(40)
      COMMON /TEXT/TEXTP
      REAL VECTP(5)
      REAL XCUR
      REAL YCUR
      COMMON /VECTRC/VECTP,XCUR,YCUR
      INTEGER IDEX(200,2)
      INTEGER NVECT(200,2)
      REAL XSIZE(200,2)
      REAL YSIZE(200,2)
      REAL X0(2300,2)
      REAL Y0(2300,2)
      REAL X1(2300,2)
      REAL Y1(2300,2)
      COMMON /FONT/IDEX,NVECT,XSIZE,YSIZE,X0,Y0,X1,Y1
      REAL GRAPHP(100)
      COMMON /GRAPH/GRAPHP
      COMMON /MAPPAR/MAPP(11)
      REAL MAPP
      COMMON /STORAG/MEMORY(1000)
      LOGICAL CPUIFC
      DIMENSION X(*),Y(*),X0T(32),Y0T(32),X1T(32),Y1T(32)
      INTEGER MASKS(1)
      include 'izbit.inc'
      DATA SYMRAT/120./

      IF (NUM.LE.0) THEN
         RETURN

      END IF

      CALL PLTSVV
      CALL PLTSVD
      CALL PLTSVT
      CALL PLTSTV(1,GRAPHP(5))
      CALL PLTSTD(1,GRAPHP(38))
      CALL PLTSTV(2,GRAPHP(63))
      IF (GRAPHP(6).EQ.1.) THEN
         CALL PLTSTT(3,0.)
         CALL PLTSTT(4,0.)
         SYMSIZ = (GRAPHP(3)+GRAPHP(4))*GRAPHP(46)/ (10.*SYMRAT)
         CALL PLTSTT(2,SYMSIZ)
         CALL PLTSTT(11,GRAPHP(66))
      END IF

      J = 0
      NP = 0
      IF (NUM.EQ.1) THEN
         IF (GRAPHP(6).EQ.1. .AND. SYMSIZ.GT.0.0) THEN
            IF (GRAPHP(21).EQ.1.) THEN
               X0T(1) = X(1)
               Y0T(1) = Y(1)

            ELSE IF (GRAPHP(21).EQ.2.) THEN
               IF (X(1).LE.0.) THEN
                  RETURN

               END IF

               X0T(1) = LOG10(X(1))
               Y0T(1) = Y(1)

            ELSE IF (GRAPHP(21).EQ.3.) THEN
               X0T(1) = X(1)
               IF (Y(1).LE.0.) THEN
                  RETURN

               END IF

               Y0T(1) = LOG10(Y(1))

            ELSE IF (GRAPHP(21).EQ.4.) THEN
               IF (X(1).LE.0.) THEN
                  RETURN

               END IF

               IF (Y(1).LE.0.) THEN
                  RETURN

               END IF

               X0T(1) = LOG10(X(1))
               Y0T(1) = LOG10(Y(1))
            END IF

            MASKS(1) = -1
            CALL PLTMP2(GRAPHP(7),1,MASKS,X0T,Y0T,X1T,Y1T)
            JB = IZBIT(1)
            IF (IAND(MASKS(1),JB).NE.0) THEN
               CALL PLTSTD(1,GRAPHP(75))
               CALL PLTXTS(X1T(1),Y1T(1),CHAR(INT(GRAPHP(47))))
            END IF

         END IF

      ELSE
         IF (GRAPHP(21).EQ.1.) THEN
            XSAV = X(1)
            YSAV = Y(1)

         ELSE IF (GRAPHP(21).EQ.2.) THEN
            DO 2000 I = 1,NUM
               IF (X(I).LE.0.) THEN
                  GO TO 2000

               END IF

               XSAV = X(I)
               YSAV = Y(I)
               GO TO 2010

 2000       CONTINUE
 2010       CONTINUE

         ELSE IF (GRAPHP(21).EQ.3.) THEN
            DO 2020 I = 1,NUM
               IF (Y(I).LE.0.) THEN
                  GO TO 2020

               END IF

               XSAV = X(I)
               YSAV = Y(I)
               GO TO 2030

 2020       CONTINUE
 2030       CONTINUE

         ELSE IF (GRAPHP(21).EQ.4.) THEN
            DO 2040 I = 1,NUM
               IF (X(I).LE.0. .OR. Y(I).LE.0.) THEN
                  GO TO 2040

               END IF

               XSAV = X(I)
               YSAV = Y(I)
               GO TO 2050

 2040       CONTINUE
 2050       CONTINUE
         END IF

 2060    IF (.NOT. (J.LT.NUM-1)) GO TO 2070
         K = MIN(NUM-J,32)
         IF (GRAPHP(21).EQ.1.) THEN
            NV = 0
            DO 2080 I = 1,K - 1
               NV = NV + 1
               X0T(I) = XSAV
               Y0T(I) = YSAV
               X1T(I) = X(I+J+1)
               Y1T(I) = Y(I+J+1)
               XSAV = X(I+J+1)
               YSAV = Y(I+J+1)
 2080       CONTINUE

         ELSE IF (GRAPHP(21).EQ.2.) THEN
            NV = 0
            DO 2100 I = 1,K - 1
               IF (X(J+I+1).LE.0.) THEN
                  GO TO 2100

               END IF

               NV = NV + 1
               X0T(NV) = LOG10(XSAV)
               Y0T(NV) = YSAV
               X1T(NV) = LOG10(X(J+I+1))
               Y1T(NV) = Y(J+I+1)
               XSAV = X(I+J+1)
               YSAV = Y(I+J+1)
 2100       CONTINUE

         ELSE IF (GRAPHP(21).EQ.3.) THEN
            NV = 0
            DO 2120 I = 1,K - 1
               IF (Y(J+I+1).LE.0.) THEN
                  GO TO 2120

               END IF

               NV = NV + 1
               X0T(NV) = XSAV
               Y0T(NV) = LOG10(YSAV)
               X1T(NV) = X(J+I+1)
               Y1T(NV) = LOG10(Y(J+I+1))
               XSAV = X(I+J+1)
               YSAV = Y(I+J+1)
 2120       CONTINUE

         ELSE IF (GRAPHP(21).EQ.4.) THEN
            NV = 0
            DO 2140 I = 1,K - 1
               IF (X(J+I+1).LE.0.) THEN
                  GO TO 2140

               END IF

               IF (Y(J+I+1).LE.0.) THEN
                  GO TO 2140

               END IF

               NV = NV + 1
               X0T(NV) = LOG10(XSAV)
               Y0T(NV) = LOG10(YSAV)
               X1T(NV) = LOG10(X(J+I+1))
               Y1T(NV) = LOG10(Y(J+I+1))
               XSAV = X(I+J+1)
               YSAV = Y(I+J+1)
 2140       CONTINUE
         END IF

         CALL PLTDV2(GRAPHP(7),NV,X0T,Y0T,X1T,Y1T)
         IF (GRAPHP(6).EQ.1. .AND. SYMSIZ.GT.0.0) THEN
            CALL PLTSTD(1,GRAPHP(75))
            MASKS(1) = -1
            XT = X1T(NV)
            YT = Y1T(NV)
            CALL PLTMP2(GRAPHP(7),NV,MASKS,X0T,Y0T,X1T,Y1T)
            DO 2160 L = 1,NV
               JB = IZBIT(L)
               IF (IAND(MASKS(1),JB).NE.0 .AND.
     *             MOD(L+NP+INT(GRAPHP(23))-1,INT(GRAPHP(23))).EQ.
     *             0) THEN
                  CALL PLTXTS(X1T(L),Y1T(L),CHAR(INT(GRAPHP(47))))
               END IF

 2160       CONTINUE
            CALL PLTSTD(1,GRAPHP(38))
         END IF

         NP = NP + NV
         J = J + K - 1
         IF (J+1.GE.NUM .AND. (GRAPHP(6).EQ.1..AND.SYMSIZ.GT.0.0)) THEN
            X0T(1) = XT
            Y0T(1) = YT
            CALL PLTSTD(1,GRAPHP(75))
            MASKS(1) = -1
            CALL PLTMP2(GRAPHP(7),1,MASKS,X0T,Y0T,X1T,Y1T)
            JB = IZBIT(1)
            IF (IAND(MASKS(1),JB).NE.0 .AND.
     *          MOD(NP+INT(GRAPHP(23)),INT(GRAPHP(23))).EQ.0) THEN
               CALL PLTXTS(X1T(1),Y1T(1),CHAR(INT(GRAPHP(47))))
            END IF

         END IF

         IF (CPUIFC(.FALSE.)) THEN
            GO TO 2070

         END IF

         GO TO 2060

 2070    CONTINUE
      END IF

      CALL PLTRET
      CALL PLTRED
      CALL PLTREV
      RETURN

      END

max_stars_repo_name: jschueller/seacas

max_issues_repo_name: jschueller/seacas

max_forks_repo_name: jschueller/seacas


text:
!>
!! @copyright Copyright (c) 2016-, Issam SAID <said.issam@gmail.com>
!! All rights reserved.
!!
!! Redistribution and use in source and binary forms, with or without
!! modification, are permitted provided that the following conditions
!! are met:
!!
!! 1. Redistributions of source code must retain the above copyright
!!    notice, this list of conditions and the following disclaimer.
!! 2. Redistributions in binary form must reproduce the above copyright
!!    notice, this list of conditions and the following disclaimer inthe
!!    documentation and/or other materials provided with the distribution.
!! 3. Neither the name of the copyright holder nor the names of its contributors
!!    may be used to endorse or promote products derived from this software
!!    without specific prior written permission.
!!
!! THIS SOFTWARE IS PROVIDED "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,
!! INCLUDING, BUT NOT LIMITED TO, WARRANTIES OF MERCHANTABILITY AND FITNESS
!! FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
!! HOLDER OR ITS CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
!! SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
!! PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
!! PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
!! LIABILITY, WETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
!! NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
!! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
!!
!! @file m_mem_test.f90
!! @author Issam SAID
!! @brief Unit testing file for the ezcu Fortran interface for CUDA
!! memory objects manipulation routines.
!<
module m_mem_test
    use, intrinsic :: iso_c_binding
    use            :: m_ezcu_flags
    use            :: m_ezcu_types
    use            :: m_ezcu_base
    use            :: m_ezcu_dev
    use            :: m_ezcu_mem
    use            :: m_ezcu_util
    use            :: m_handler

    implicit none

    private

    public :: mem_test

    integer, parameter :: N = 64
    type(ezcu_dev_t), pointer :: d

contains

    logical function wrap_buffer() result(status)
        real :: h(N)
        call ezcu_mem_wrap(h, d, N, FLOAT)
        status = associated(d)
    end function wrap_buffer

    logical function wrap_allocatable() result(status)
        real, allocatable :: h(:)
        allocate(h(N))
        call ezcu_mem_wrap(h, d, FLOAT)
        status = associated(d)
        deallocate(h)
    end function wrap_allocatable

    logical function wrap_pointer() result(status)
        real, pointer, dimension(:) :: h
        allocate(h(N))
        call ezcu_mem_wrap_ptr(h, d, FLOAT)
        status = associated(d)
        deallocate(h)
    end function wrap_pointer

    logical function read_access_buffer() result(status)
        real, allocatable :: h(:)
        allocate(h(N))
        call ezcu_mem_wrap(h, d, READ_WRITE)
        call ezcu_mem_update(h, READ_ONLY)
        status = associated(d)
        deallocate(h)
    end function read_access_buffer

    logical function write_access_buffer() result(status)
        real :: h(-3:N-4)
        h(-3:N-4) = 0.0
        h(-3)     = 1.2
        call ezcu_mem_wrap(h, d, N, READ_WRITE)
        call ezcu_mem_update(h, N, READ_ONLY)
        status = h(-3).eq.1.2

    end function write_access_buffer

    logical function write_access_allocatable() result(status)
        real, allocatable :: h(:)
        allocate(h(-3:N-4))
        h(-3:N-4) = 0.0
        h(-3)     = 1.2
        call ezcu_mem_wrap(h, d, READ_WRITE)
        call ezcu_mem_update(h, READ_ONLY)
        status = h(-3).eq.1.2

        deallocate(h)
    end function write_access_allocatable

    logical function write_access_pointer() result(status)
        real, pointer :: h(:)
        allocate(h(-3:N-4))
        h(-3:N-4) = 0.0
        h(-3)     = 1.2
        h(N-10)   = -1.2
        call ezcu_mem_wrap_ptr(h, d, READ_WRITE)
        call ezcu_mem_update_ptr(h, READ_ONLY)
        status = h(-3).eq.1.2.and. h(N-10).eq.-1.2

        deallocate(h)
    end function write_access_pointer

    logical function access_buffer_2d() result(status)
        real, allocatable :: h2d(:,:)
        integer :: i, j, H, W
        H = 16
        W = 4
        status=.true.
        allocate(h2d(H,W))
        call ezcu_mem_wrap(h2d, d, READ_WRITE)
        call ezcu_mem_update(h2d, WRITE_ONLY)
        do j=1,W
            do i=1,H
                h2d(i,j) = (j-1)*H+i
            enddo
        enddo
        call ezcu_mem_update(h2d, READ_ONLY)
        do j=1,W
            do i=1,H
                status = status .and. (abs(h2d(i,j) - ((j-1)*H+i)) < 1.e-4)
            enddo
        enddo
        deallocate(h2d)
    end function access_buffer_2d

    logical function access_buffer_3d() result(status)
        real, allocatable :: h3d(:,:,:)
        integer :: TH, TW, TD
        TH = 16
        TW = 16
        TD = 16

        allocate(h3d(-3:TH+4, -3:TW+4, -3:TD+4))
        h3d = 0.
        h3d(1:TH, 1:TW, 1:TD) = 1.2
        call ezcu_mem_wrap(h3d, d, READ_WRITE)
        call ezcu_mem_update(h3d, READ_ONLY)

        status = .true.
        deallocate(h3d)
    end function access_buffer_3d

    logical function access_buffer_4d() result(status)
        real, allocatable, target         :: h4d(:,:,:,:)
        real, pointer, dimension(:,:,:,:) :: hptr
        integer, parameter                :: TH=16, TW=16, TD=16, TK=32
        real                              :: h(TH,TW,TD,TK)

        allocate(hptr(4,4,4,4))
        allocate(h4d(-3:TH+4, -3:TW+4, -3:TD+4, -3:TK+4))
        h4d = 0.
        h4d(1:TH, 1:TW, 1:TD, 1:TK) = 1.2
        call ezcu_mem_wrap(h4d, d, READ_WRITE)
        call ezcu_mem_wrap_ptr(hptr, d, READ_WRITE)
        call ezcu_mem_wrap(h, d, TH, TW, TD, TK, READ_WRITE)
        call ezcu_mem_update(h4d, READ_ONLY)
        call ezcu_mem_update_ptr(hptr, READ_ONLY)
        call ezcu_mem_update(h, TH, TW, TD, TK, READ_ONLY)

        status = .true.
        deallocate(h4d)
        deallocate(hptr)
    end function access_buffer_4d

    subroutine setup()
        call ezcu_init()
        call ezcu_dev_find(0, d)
    end subroutine setup

    subroutine teardown()
        call ezcu_release()
    end subroutine teardown

    subroutine mem_test()
        call run(setup, teardown, wrap_buffer, &
              "mem_test.wrap_buffer")
        call run(setup, teardown, wrap_allocatable, &
              "mem_test.wrap_allocatable")
        call run(setup, teardown, wrap_pointer, &
              "mem_test.wrap_pointer")
        call run(setup, teardown, read_access_buffer, &
               "mem_test.read_access_buffer")
        call run(setup, teardown, write_access_buffer, &
              "mem_test.write_access_buffer")
        call run(setup, teardown, write_access_allocatable, &
              "mem_test.write_access_allocatable")
        call run(setup, teardown, write_access_pointer, &
              "mem_test.write_access_pointer")
        call run(setup, teardown, access_buffer_2d, &
              "mem_test.access_buffer_2d")
        call run(setup, teardown, access_buffer_3d, &
              "mem_test.access_buffer_3d")
        call run(setup, teardown, access_buffer_4d, &
              "mem_test.access_buffer_4d")
    end subroutine mem_test

end module m_mem_test

max_stars_repo_name: issamsaid/ezCU

max_issues_repo_name: issamsaid/ezCU

max_forks_repo_name: issamsaid/ezCU


text:
integer grid,gridz

real*4,allocatable ::   rho(:,:,:)

double precision ss,ss2,totval

character*80  filename

print*, ' Input densfield'

read (*,'(A)') filename

open(1,file=filename, form='unformatted',status='old')

read (1) grid,grid, gridz

allocate(rho(grid,grid,gridz))

 do kk=1,gridz
    read(1)  ((RHO(II,JJ,KK), ii=1,GRID),jj=1,GRID)
 enddo

close(1)

ss=0.

!$OMP PARALLEL DO DEFAULT (SHARED) PRIVATE(I,J,K) REDUCTION(+:ss)
do K=1,grid
   do J=1,grid
      do I=1,grid
         ss = ss + rho(i,j,k)
      enddo
   enddo
enddo

totval = float(grid)**3

ss = ss / totval

print*, ' Mean valuel=', ss

ss2=0.

!$OMP PARALLEL DO DEFAULT (SHARED) PRIVATE(I,J,K) REDUCTION(+:ss2)
do K=1,grid
   do J=1,grid
      do I=1,grid
         ss2 = (ss-rho(i,j,k))**2 + ss2
      enddo
   enddo
enddo

print*, ' Min-Max values ', Minval(rho), Maxval(rho)

ss2 = sqrt(ss2) / totval

print*, ' RMS=', ss2

stop
end

max_stars_repo_name: JohanComparat/nbody-npt-functions

max_issues_repo_name: JohanComparat/nbody-npt-functions

max_forks_repo_name: JohanComparat/nbody-npt-functions


text:
module unit_test_mod
  implicit none
contains

  function test_array_io() result(nerr) bind(c)
    use array_io_mod
    use iso_c_binding

    integer(kind=c_int) :: nerr
    integer ::  i, j
    real, target :: a(10,3), b(10,3)
    logical :: ok

    character(kind=c_char, len=128), parameter :: &
         filename = c_char_"unit_test_f90_array_io.dat"//C_NULL_CHAR

    do j = 1,3
       do i = 1,10
          a(i,j) = 100*j + i
       end do
    end do
    nerr = 0
    ok = array_io_write(filename, c_loc(a), size(a))
    if (.not. ok) nerr = nerr + 1
    ok = array_io_file_exists(filename)
    if (.not. ok) nerr = nerr + 1
    ok = array_io_read(filename, c_loc(b), size(b))
    if (.not. ok) nerr = nerr + 1
    do j = 1,3
       do i = 1,10
          if (a(i,j) .ne. b(i,j)) nerr = nerr + 1
       end do
    end do
  end function test_array_io

end module unit_test_mod

max_stars_repo_name: AaronDonahue/scream

max_stars_repo_licenses: zlib-acknowledgement FTL RSA-MD

max_issues_repo_name: AaronDonahue/scream

max_issues_repo_licenses: zlib-acknowledgement FTL RSA-MD

max_forks_repo_name: AaronDonahue/scream

max_forks_repo_licenses: zlib-acknowledgement FTL RSA-MD


text:
subroutine rpn2(ixy,maxm,meqn,mwaves,maux,mbc,mx,ql,qr,auxl,auxr,fwave,s,amdq,apdq)

    use amr_module, only: mcapa

    use geoclaw_module, only: g => grav, dry_tolerance
    use geoclaw_module, only: coordinate_system, earth_radius, deg2rad

    implicit none

    ! Input
    integer, intent(in) :: ixy, maxm, mx, meqn, mwaves, mbc, maux
    real(kind=8), intent(in) :: ql(meqn, 1-mbc:maxm+mbc)
    real(kind=8), intent(in) :: qr(meqn, 1-mbc:maxm+mbc)
    real(kind=8), intent(in) :: auxl(maux, 1-mbc:maxm+mbc)
    real(kind=8), intent(in) :: auxr(maux, 1-mbc:maxm+mbc)

    ! Output
    real(kind=8), intent(in out) :: s(mwaves, 1-mbc:maxm+mbc)
    real(kind=8), intent(in out) :: fwave(meqn, mwaves, 1-mbc:maxm+mbc)
    real(kind=8), intent(in out) :: amdq(meqn,1-mbc:maxm+mbc)
    real(kind=8), intent(in out) :: apdq(meqn,1-mbc:maxm+mbc)

    ! Locals
    integer :: normal_dir, tangential_dir, i, mw
    real(kind=8) :: h_l, h_r, hu_l, hu_r, hv_l, hv_r, u_l, u_r, v_l, v_r
    real(kind=8) :: b_l, b_r, phi_l, phi_r

    real(kind=8) :: h_hat, u_hat, v_hat, c_hat

    real(kind=8) :: delta(3), beta(3)
    real(kind=8) :: dxdc

    ! Normal vs. Tangential directions
    if (ixy == 1) then
        normal_dir = 2
        tangential_dir = 3
    else
        normal_dir = 3
        tangential_dir = 2
    endif

    ! Main loop
    do i = 2 - mbc, mx + mbc

        ! Skip this cell interface if both cells are dry
        if (qr(1, i - 1) < dry_tolerance .and. ql(1, i) < dry_tolerance) then
            cycle
        end if

        ! Extract states
        h_l = qr(1, i - 1)
        h_r = ql(1, i)
        hu_l = qr(normal_dir, i - 1)
        hu_r = ql(normal_dir, i)
        hv_l = qr(tangential_dir, i - 1)
        hv_r = ql(tangential_dir, i)
        b_l = auxr(1, i - 1)
        b_r = auxl(1, i)

        ! Determine velocities
        if (h_l < dry_tolerance) then
            h_l = 0.d0
            u_l = 0.d0
            v_l = 0.d0
        else
            u_l = hu_l / h_l
            v_l = hv_l / h_l
        end if
        if (h_r < dry_tolerance) then
            h_r = 0.d0
            u_r = 0.d0
            v_r = 0.d0
        else
            u_r = hu_r / h_r
            v_r = hv_r / h_r
        end if

        ! Momentum flux
        phi_l = h_l * u_l**2 + 0.5d0 * g * h_l**2
        phi_r = h_r * u_r**2 + 0.5d0 * g * h_r**2

        ! Compute Roe averages
        h_hat = 0.5d0 * (h_r + h_l)
        u_hat = (sqrt(h_r) * u_r + sqrt(h_l) * u_l) / (sqrt(h_r) + sqrt(h_l))
        v_hat = (sqrt(h_r) * v_r + sqrt(h_l) * v_l) / (sqrt(h_r) + sqrt(h_l))
        c_hat = sqrt(g * h_hat)

        ! Flux differences
        delta(1) = hu_r - hu_l
        delta(normal_dir) = phi_r - phi_l + 0.5d0 * g * (h_r + h_l) * (b_r - b_l)
        delta(tangential_dir) = h_r * u_r * v_r - h_l * u_l * v_l

        ! Wave speeds - Einfeldt speeds
        s(1, i) = min(u_hat - c_hat, u_l - sqrt(g * h_l))
        s(3, i) = max(u_hat + c_hat, u_r + sqrt(g * h_r))
        s(2, i) = 0.5d0 * (s(1, i) + s(3, i))

        ! Wave strengths
        beta(1) = (s(3, i) * delta(1) - delta(2)) / (s(3, i) - s(1, i))
        beta(3) = (delta(2) - s(1, i) * delta(1)) / (s(3, i) - s(1, i))
        beta(2) = delta(3) - v_l * beta(1) - v_r * beta(3)

        ! F-Waves
        fwave(:, 1, i) = [1.d0, s(1, i), v_l] * beta(1)
        fwave(:, 2, i) = [0.d0, 0.d0, 1.d0] * beta(2)
        fwave(:, 3, i) = [1.d0, s(3, i), v_r] * beta(3)

    end do
    ! End of main loop

    ! Capacity mapping for lat-long grids
    if (mcapa > 0) then
        do i = 2 - mbc, mx + mbc
            if(ixy == 1) then
                dxdc = earth_radius * deg2rad
            else
                dxdc = earth_radius * cos(auxl(3, i)) * deg2rad
            end if

            s(:, i) = dxdc * s(:, i)
            fwave(:, :, i) = dxdc * fwave(:, :, i)
        end do
    end if

    ! Accumulate fluctuations
    do i = 2 - mbc, mx + mbc
        do mw = 1, mwaves
            if (s(mw, i) < 0.d0) then
                amdq(:, i) = amdq(:, i) + fwave(:, mw, i)
            else if (s(mw, i) > 0.d0) then
                apdq(:, i) = apdq(:, i) + fwave(:, mw, i)
            else
                amdq(:, i) = amdq(:, i) + 0.5d0 * fwave(:, mw, i)
                apdq(:, i) = apdq(:, i) + 0.5d0 * fwave(:, mw, i)
            end if
        end do
    end do

    print *, s
    stop

end subroutine rpn2

max_stars_repo_name: mandli/fwave-swe

max_issues_repo_name: mandli/fwave-swe

max_forks_repo_name: mandli/fwave-swe


text: http://visit.ucdavis.edu/TourReg UC Davis Barn Tours are available for groups of preschoolers between 20 and 30 in size. Call (530) 7528111 for more information. If you have to entertain and educate a large group of kids this could be a Good Places to Take Kids good place to take them.


text:
!------------------------------------------------------
! Pijl et al, 4.1.2. Zalesak's rotating disk
!------------------------------------------------------

  FUNCTION InitZalesak( Model, n, t ) RESULT(f)
    USE Types
    USE SolverUtils
    IMPLICIT NONE

    TYPE(Model_t) :: Model
    LOGICAL :: Visited = .FALSE., GotIt
    INTEGER :: n,i,j
    REAL(KIND=dp) :: t,x,y,r,Lx,Ly,x0,y0,r0,w0,dist(4),f,a1,a2,b1,b2,f1,f2

    SAVE Visited, x0, y0, r0, w0, a1, a2, b1, b2

    IF(.NOT. Visited) THEN
      Lx = ListGetConstReal(Model % Simulation,'Zalesak Lx',GotIt)
      IF(.NOT. GotIt) Lx = MAXVAL( Model % Nodes % x )
      Ly = ListGetConstReal(Model % Simulation,'Zalesak Ly',GotIt)
      IF(.NOT. GotIt) Ly = MAXVAL( Model % Nodes % y )

      x0 = 0.5 * Lx
      y0 = 0.75 * Ly
      r0 = (3.0/20.0) * Lx
      w0 = r0 / 3.0d0

      a1 = x0 - w0/2
      a2 = x0 + w0/2
      b1 = -1.0
      b2 = y0 + r0 - w0
      Visited = .TRUE.
    END IF

    x = Model % Nodes % x(n)
    y = Model % Nodes % y(n)

    ! distance from horizontal lines
    IF( x < a1 ) THEN
      dist(1) = SQRT( (x-a1)**2.0 + (y-b1)**2.0)
      dist(2) = SQRT( (x-a1)**2.0 + (y-b2)**2.0)
    ELSE IF( x > a2 ) THEN
      dist(1) = SQRT( (x-a2)**2.0 + (y-b1)**2.0)
      dist(2) = SQRT( (x-a2)**2.0 + (y-b2)**2.0)
    ELSE
      dist(1) = ABS(y - b1)
      dist(2) = ABS(y - b2)
    END IF

    ! distance from vertical lines
    IF( y < b1 ) THEN
      dist(3) = SQRT( (x-a1)**2.0 + (y-b1)**2.0)
      dist(4) = SQRT( (x-a2)**2.0 + (y-b1)**2.0)
    ELSE IF( y > b2 ) THEN
      dist(3) = SQRT( (x-a1)**2.0 + (y-b2)**2.0)
      dist(4) = SQRT( (x-a2)**2.0 + (y-b2)**2.0)
    ELSE
      dist(3) = ABS(x - a1)
      dist(4) = ABS(x - a2)
    END IF

    ! minimum distance to the sides of the rectangle
    f1 = MINVAL( ABS( dist ) )

    ! set negative sign inside the rectangle
    IF( x > a1 .AND. x < a2 .AND. y > b1 .AND. y < b2) THEN
      f1 = -f1
    END IF

    ! inside of a circle
    f2 = r0 - SQRT( (x-x0)**2.0 + (y-y0)**2.0 )

    ! take a union of the rectangle and circle
    IF( f2 > 0.0 ) THEN
      IF( f1 < 0.0 ) THEN
        f = f1
      ELSE
        f = MIN(f1,f2)
      END IF
    ELSE
      IF( f1 > 0.0 ) THEN
        f = f2
      ELSE
        f = MIN(f1,f2)
      END IF
    END IF

  END FUNCTION InitZalesak

  FUNCTION RotateVeloX( Model, n, t ) RESULT(f)
    USE Types
    USE SolverUtils
    IMPLICIT NONE

    TYPE(Model_t) :: Model
    LOGICAL :: Visited = .FALSE., GotIt
    INTEGER :: n,i,j
    REAL(KIND=dp) :: t,x,y,r,Lx,Ly,x0,y0,f

    SAVE Visited, x0, y0

    IF(.NOT. Visited) THEN
      Lx = ListGetConstReal(Model % Simulation,'Zalesak Lx',GotIt)
      IF(.NOT. GotIt) Lx = MAXVAL( Model % Nodes % x )
      Ly = ListGetConstReal(Model % Simulation,'Zalesak Ly',GotIt)
      IF(.NOT. GotIt) Ly = MAXVAL( Model % Nodes % y )

      x0 = Lx / 2.0
      y0 = Ly / 2.0
      Visited = .TRUE.
    END IF

    x = Model % Nodes % x(n)
    y = Model % Nodes % y(n)

    f = -(y-y0)

  END FUNCTION RotateVeloX

  FUNCTION RotateVeloY( Model, n, t ) RESULT(f)
    USE Types
    USE SolverUtils
    IMPLICIT NONE

    TYPE(Model_t) :: Model
    LOGICAL :: Visited = .FALSE., GotIt
    INTEGER :: n,i,j
    REAL(KIND=dp) :: t,x,y,r,Lx,Ly,x0,y0,f

    SAVE Visited, x0, y0

    IF(.NOT. Visited) THEN
      Lx = ListGetConstReal(Model % Simulation,'Zalesak Lx',GotIt)
      IF(.NOT. GotIt) Lx = MAXVAL( Model % Nodes % x )
      Ly = ListGetConstReal(Model % Simulation,'Zalesak Ly',GotIt)
      IF(.NOT. GotIt) Ly = MAXVAL( Model % Nodes % y )

      x0 = Lx / 2.0
      y0 = Ly / 2.0
      Visited = .TRUE.
    END IF

    x = Model % Nodes % x(n)
    y = Model % Nodes % y(n)

    f = x - x0

  END FUNCTION RotateVeloY

max_stars_repo_name: readex-eu/readex-apps

max_issues_repo_name: jcmcmurry/pipelining

max_forks_repo_name: jcmcmurry/pipelining