!> This subroutine is used for calculating energy
!> dispersion with wannier functions
  subroutine ek_slab

     use mpi
     use para
     implicit none

! loop index
     integer :: i
     integer :: j
     integer :: l
     integer :: i1
     integer :: i2
     integer :: NN, nlines, knv3

     integer :: lwork

! wave vector
     real(dp) :: k1(2)
     real(dp) :: k2(2)
     real(Dp) :: k(2), kstart(2), kend(2)
     real(dp) :: kp(16,2), ke(16,2), kpath_stop(16)
     character(4) :: kpath_name(17)

     real(dp) :: cell_volume
     real(dp) :: t1, temp
     real(dp) :: emin, emax

! parameters for zheev
     integer :: info,ierr
     real(Dp), allocatable ::  rwork(:)
     complex(Dp), allocatable :: work(:)

! eigenvalue
     real(Dp) :: eigenvalue(nslab*Num_wann)

! energy dispersion
     real(Dp),allocatable :: ekslab(:,:)
     real(Dp),allocatable :: ekslab_mpi(:,:)

     real(dp), allocatable :: kpoint(:,:)
     real(dp), allocatable :: k_len(:)

     !> color for plot, surface state weight
     real(dp), allocatable :: surf_weight(:, :)
     real(dp), allocatable :: surf_weight_mpi(:, :)

! hamiltonian slab
     complex(Dp),allocatable ::CHamk(:,:)

     lwork= 16*Nslab*Num_wann
     ierr = 0

     kpath_name= ' '
     kp(1,:)=(/0.5d0, 0.0d0/)  ; kpath_name(1)= 'X'
     ke(1,:)=(/0.0d0, 0.0d0/)
     kp(2,:)=(/0.0d0, 0.0d0/)  ; kpath_name(2)= 'G'
     ke(2,:)=(/0.0d0, 0.50d0/)  ! K
     kp(3,:)=(/0.0d0, 0.50d0/) ; kpath_name(3)= 'Y'
     ke(3,:)=(/0.5d0, 0.5d0/)  ! K
     kp(4,:)=(/0.5d0, 0.5d0/)  ; kpath_name(4)= 'M'
     ke(4,:)=(/0.0d0, 0.0d0/)  ; kpath_name(5)= 'G'

     kp(5,:)=(/0.0d0, 0.0d0/)  ! K
     ke(5,:)=(/0.5d0, 0.5d0/)  ! K
     kp(6,:)=(/0.0d0, 0.0d0/)  ! K
     ke(6,:)=(/1.0d0, 1.0d0/)  ! K

     nlines= 4
     NN=Nk
     knv3=NN*nlines
     allocate( kpoint(knv3, 3))
     allocate( k_len (knv3))
     allocate( surf_weight (Nslab* Num_wann, knv3))
     allocate( surf_weight_mpi (Nslab* Num_wann, knv3))
     allocate(ekslab(Nslab*Num_wann,knv3))
     allocate(ekslab_mpi(Nslab*Num_wann,knv3))
     allocate(CHamk(nslab*Num_wann,nslab*Num_wann))
     allocate(work(lwork))
     allocate(rwork(lwork))

     kpoint= 0d0
     surf_weight= 0d0
     surf_weight_mpi= 0d0

     t1=0d0
     k_len=0d0
     kpath_stop= 0d0
     do j=1, nlines
        do i=1, NN
           k = kp(j,:)
           kstart=k
           k1= kstart(1)*Ka2+ kstart(2)*Kb2

           k = ke(j,:)
           kend=k
           k2= kend(1)*Ka2+ kend(2)*Kb2

           kpoint(i+(j-1)*NN,:)= kstart+ (kend-kstart)*dble(i-1)/dble(NN-1)

           temp= dsqrt((k2(1)- k1(1))**2 &
                 +(k2(2)- k1(2))**2)/dble(NN-1)

           if (i.gt.1) then
              t1=t1+temp
           endif
           k_len(i+(j-1)*NN)= t1
        enddo
        kpath_stop(j+1)= t1
     enddo

     ! sweep k
     ekslab=0.0d0
     do i=1+cpuid,knv3,num_cpu
        if (cpuid==0) print *, 'ik ',  i, knv3
        k= kpoint(i, :)
        chamk=0.0d0

        !> no magnetgic field
        if (abs(Bx)<eps9.and. abs(By)<eps9.and. abs(Bz)<eps9)then
           call ham_slab(k,Chamk)
        !> in-plane magnetic field
        elseif (abs(Bx)>eps9 .or. abs(By)>eps9)then
           call ham_slab_parallel_B(k,Chamk)
        !> vertical magnetic field
        else
           print *, 'Error: we only support in-plane magnetic field at present'
           stop 'please set Bz= 0'
        endif

        eigenvalue=0.0d0

        ! diagonal Chamk
        call eigensystem_c('V', 'U', Num_wann*Nslab, CHamk, eigenvalue)

        ekslab(:,i)=eigenvalue

        do j=1, Nslab* Num_wann
           do l=1, Num_wann
              surf_weight(j, i)= surf_weight(j, i) &
                 + abs(CHamk(l, j))**2 & ! first slab
                 + abs(CHamk(Num_wann*Nslab- l+ 1, j))**2 !& ! last slab
                !+ abs(CHamk(Num_wann+ l, j))**2 & ! the second slab
                !+ abs(CHamk(Num_wann*(Nslab-1)- l, j))**2 ! last second slab
           enddo ! l
           surf_weight(j, i)= sqrt(surf_weight(j, i))
        enddo ! j
        if (cpuid==0) write(stdout, *)'SlabEk,k', i, knv3
     enddo ! i
     call mpi_allreduce(ekslab,ekslab_mpi,size(ekslab),&
                       mpi_dp,mpi_sum,mpi_cmw,ierr)
     call mpi_allreduce(surf_weight, surf_weight_mpi,size(surf_weight),&
                       mpi_dp,mpi_sum,mpi_cmw,ierr)

     ekslab=ekslab_mpi
     surf_weight= surf_weight_mpi/ maxval(surf_weight_mpi)

     if(cpuid==0)then
        open(unit=100, file='slabek.dat')
        do j=1, Num_wann*Nslab
           do i=1, knv3
             !write(100,'(3f15.7, i8)')k_len(i), ekslab(j,i), &
             !   (surf_weight(j, i))
              write(100,'(2f15.7, i8)')k_len(i), ekslab(j,i), &
                 int(255-surf_weight(j, i)*255d0)
           enddo
           write(100 , *)''
        enddo
        close(100)
        write(stdout,*) 'calculate energy band  done'
     endif

     emin= minval(ekslab)-0.5d0
     emax= maxval(ekslab)+0.5d0
     !> write script for gnuplot
     if (cpuid==0) then
        open(unit=101, file='slabek.gnu')
        write(101, '(a)')'#set terminal  postscript enhanced color'
        write(101, '(a)')"#set output 'slabek.eps'"
        write(101, '(3a)')'set terminal  pngcairo truecolor enhanced', &
           '  font ",60" size 1920, 1680'
        write(101, '(a)')"set output 'slabek.png'"
        write(101,'(2a)') 'set palette defined ( 0  "green", ', &
           '5 "yellow", 10 "red" )'
        write(101, '(a)')'set style data linespoints'
        write(101, '(a)')'unset ztics'
        write(101, '(a)')'unset key'
        write(101, '(a)')'set pointsize 0.8'
        write(101, '(a)')'set border lw 3 '
        write(101, '(a)')'set view 0,0'
        write(101, '(a)')'#set xtics font ",36"'
        write(101, '(a)')'#set ytics font ",36"'
        write(101, '(a)')'#set ylabel font ",36"'
        write(101, '(a)')'#set xtics offset 0, -1'
        write(101, '(a)')'set ylabel offset -1, 0 '
        write(101, '(a)')'set ylabel "Energy (eV)"'
        write(101, '(a, f10.5, a)')'set xrange [0: ', maxval(k_len), ']'
        write(101, '(a, f10.5, a, f10.5, a)')'set yrange [', emin, ':', emax, ']'
        write(101, 202, advance="no") (trim(kpath_name(i)), kpath_stop(i), i=1, nlines)
        write(101, 203)trim(kpath_name(nlines+1)), kpath_stop(nlines+1)

        do i=1, nlines-1
           write(101, 204)kpath_stop(i+1), emin, kpath_stop(i+1), emax
        enddo
        write(101, '(a)')'rgb(r,g,b) = int(r)*65536 + int(g)*256 + int(b)'
        write(101, '(2a)')"plot 'slabek.dat' u 1:2:(rgb(255,$3, 3)) ",  &
            "w lp lw 2 pt 7  ps 1 lc rgb variable"

        !write(101, '(2a)')"splot 'slabek.dat' u 1:2:3 ",  &
        !   "w lp lw 2 pt 13 palette"
     endif

     202 format('set xtics (',:20('"',A3,'" ',F10.5,','))
     203 format(A3,'" ',F10.5,')')
     204 format('set arrow from ',F10.5,',',F10.5, &
        ' to ',F10.5,',',F10.5, ' nohead')

  return
  end subroutine ek_slab

  !> This subroutine is used for calculating energy dispersion
  !> with wannier functions in plane magnetic field
  !>
  subroutine ek_slab_b

     use mpi
     use para
     implicit none

! loop index
     integer :: i
     integer :: j
     integer :: l
     integer :: i1
     integer :: i2
     integer :: NN, nlines, knv3

     integer :: lwork

! wave vector
     real(dp) :: k1(2)
     real(dp) :: k2(2)
     real(Dp) :: k(2), kstart(2), kend(2)
     real(dp) :: kp(16,2), ke(16,2), kpath_stop(16)
     character(4) :: kpath_name(17)

     real(dp) :: cell_volume
     real(dp) :: t1, temp
     real(dp) :: emin, emax

! parameters for zheev
     integer :: info,ierr
     real(Dp), allocatable ::  rwork(:)
     complex(Dp), allocatable :: work(:)

! eigenvalue
     real(Dp) :: eigenvalue(nslab*Num_wann)

! energy dispersion
     real(Dp),allocatable :: ekslab(:,:)
     real(Dp),allocatable :: ekslab_mpi(:,:)

     real(dp), allocatable :: kpoint(:,:)
     real(dp), allocatable :: k_len(:)

     !> color for plot, surface state weight
     real(dp), allocatable :: surf_weight(:, :)
     real(dp), allocatable :: surf_weight_mpi(:, :)

! hamiltonian slab
     complex(Dp),allocatable ::CHamk(:,:)

     lwork= 16*Nslab*Num_wann
     ierr = 0

     kpath_name= ' '
    !kp(1,:)=(/0.7d0, 0.0d0/)  ; kpath_name(1)= 'X2'
    !ke(1,:)=(/0.5d0, 0.0d0/)
    !kp(2,:)=(/0.5d0, 0.0d0/)  ; kpath_name(2)= 'X'
    !ke(2,:)=(/0.5d0, 0.20d0/)  ! K
     kp(1,:)=(/0.0d0, 0.0d0/)  ; kpath_name(1)= 'G'
     ke(1,:)=(/0.5d0, 0.0d0/)
     kp(2,:)=(/0.5d0, 0.0d0/)  ; kpath_name(2)= 'X'
     ke(2,:)=(/0.5d0, 0.50d0/)  ! K
     kp(3,:)=(/0.5d0, 0.50d0/) ; kpath_name(3)= 'M'
     ke(3,:)=(/0.0d0, 0.0d0/)  ! K
     kp(4,:)=(/0.0d0, 0.0d0/)  ; kpath_name(4)= 'G'
     ke(4,:)=(/0.0d0, 0.5d0/)  ; kpath_name(5)= 'Y'

     kp(5,:)=(/0.0d0, 0.0d0/)  ! K
     ke(5,:)=(/0.5d0, 0.5d0/)  ! K
     kp(6,:)=(/0.0d0, 0.0d0/)  ! K
     ke(6,:)=(/1.0d0, 1.0d0/)  ! K

     nlines= 4
     NN=Nk
     knv3=NN*nlines
     allocate( kpoint(knv3, 3))
     allocate( k_len (knv3))
     allocate( surf_weight (Nslab* Num_wann, knv3))
     allocate( surf_weight_mpi (Nslab* Num_wann, knv3))
     allocate(ekslab(Nslab*Num_wann,knv3))
     allocate(ekslab_mpi(Nslab*Num_wann,knv3))
     allocate(CHamk(nslab*Num_wann,nslab*Num_wann))
     allocate(work(lwork))
     allocate(rwork(lwork))

     kpoint= 0d0
     surf_weight= 0d0
     surf_weight_mpi= 0d0

     t1=0d0
     k_len=0d0
     kpath_stop= 0d0
     do j=1, nlines
        do i=1, NN
           k = kp(j,:)
           kstart=k
           k1= kstart(1)*Ka2+ kstart(2)*Kb2

           k = ke(j,:)
           kend=k
           k2= kend(1)*Ka2+ kend(2)*Kb2

           kpoint(i+(j-1)*NN,:)= kstart+ (kend-kstart)*dble(i-1)/dble(NN-1)

           temp= dsqrt((k2(1)- k1(1))**2 &
                 +(k2(2)- k1(2))**2)/dble(NN-1)

           if (i.gt.1) then
              t1=t1+temp
           endif
           k_len(i+(j-1)*NN)= t1
        enddo
        kpath_stop(j+1)= t1
     enddo

     ! sweep k
     ekslab=0.0d0
     do i=1+cpuid,knv3,num_cpu
        if (cpuid==0) print *, 'ik ',  i, knv3
        k= kpoint(i, :)
        chamk=0.0d0
        eigenvalue=0.0d0
        ! diagonal Chamk
        call eigensystem_c('V', 'U', Num_wann*Nslab, CHamk, eigenvalue)

        ekslab(:,i)=eigenvalue

        do j=1, Nslab* Num_wann
           do l=1, Num_wann
              surf_weight(j, i)= surf_weight(j, i) &
                 + abs(CHamk(l, j))**2 & ! first slab
                 + abs(CHamk(Num_wann*Nslab- l+ 1, j))**2 !& ! last slab
                !+ abs(CHamk(Num_wann+ l, j))**2 & ! the second slab
                !+ abs(CHamk(Num_wann*(Nslab-1)- l+ 1, j))**2 ! last second slab
           enddo ! l
           surf_weight(j, i)= sqrt(surf_weight(j, i))
        enddo ! j
        if (cpuid==0) write(stdout, *)'SlabEk,k', i, knv3
     enddo ! i
     call mpi_allreduce(ekslab,ekslab_mpi,size(ekslab),&
                       mpi_dp,mpi_sum,mpi_cmw,ierr)
     call mpi_allreduce(surf_weight, surf_weight_mpi,size(surf_weight),&
                       mpi_dp,mpi_sum,mpi_cmw,ierr)

     ekslab=ekslab_mpi
     surf_weight= surf_weight_mpi/ maxval(surf_weight_mpi)

     if(cpuid==0)then
        open(unit=100, file='slabek.dat')
        do j=1, Num_wann*Nslab
           do i=1, knv3
             !write(100,'(3f15.7, i8)')k_len(i), ekslab(j,i), &
             !   (surf_weight(j, i))
              write(100,'(2f15.7, i8)')k_len(i), ekslab(j,i), &
                 int(255-surf_weight(j, i)*255d0)
           enddo
           write(100 , *)''
        enddo
        close(100)
        write(stdout,*) 'calculate energy band  done'
     endif

     emin= minval(ekslab)-0.5d0
     emax= maxval(ekslab)+0.5d0
     !> write script for gnuplot
     if (cpuid==0) then
        open(unit=101, file='slabek.gnu')
        write(101, '(a)')'#set terminal  postscript enhanced color'
        write(101, '(a)')"#set output 'slabek.eps'"
        write(101, '(3a)')'set terminal  pngcairo truecolor enhanced', &
           '  font ",60" size 1920, 1680'
        write(101, '(a)')"set output 'slabek.png'"
        write(101,'(2a)') 'set palette defined ( 0  "green", ', &
           '5 "yellow", 10 "red" )'
        write(101, '(a)')'set style data linespoints'
        write(101, '(a)')'unset ztics'
        write(101, '(a)')'unset key'
        write(101, '(a)')'set pointsize 0.8'
        write(101, '(a)')'set border lw 3 '
        write(101, '(a)')'set view 0,0'
        write(101, '(a)')'#set xtics font ",36"'
        write(101, '(a)')'#set ytics font ",36"'
        write(101, '(a)')'#set ylabel font ",36"'
        write(101, '(a)')'#set xtics offset 0, -1'
        write(101, '(a)')'set ylabel offset -1, 0 '
        write(101, '(a)')'set ylabel "Energy (eV)"'
        write(101, '(a, f10.5, a)')'set xrange [0: ', maxval(k_len), ']'
        write(101, '(a, f10.5, a, f10.5, a)')'set yrange [', emin, ':', emax, ']'
        write(101, 202, advance="no") (trim(kpath_name(i)), kpath_stop(i), i=1, nlines)
        write(101, 203)trim(kpath_name(nlines+1)), kpath_stop(nlines+1)

        do i=1, nlines-1
           write(101, 204)kpath_stop(i+1), emin, kpath_stop(i+1), emax
        enddo
        write(101, '(a)')'rgb(r,g,b) = int(r)*65536 + int(g)*256 + int(b)'
        write(101, '(2a)')"plot 'slabek.dat' u 1:2:(rgb(255,$3, 3)) ",  &
            "w lp lw 2 pt 7  ps 1 lc rgb variable"

        !write(101, '(2a)')"splot 'slabek.dat' u 1:2:3 ",  &
        !   "w lp lw 2 pt 13 palette"
     endif

     202 format('set xtics (',:20('"',A3,'" ',F10.5,','))
     203 format(A3,'" ',F10.5,')')
     204 format('set arrow from ',F10.5,',',F10.5,' to ',F10.5,',',F10.5, ' nohead')

  return
  end subroutine ek_slab_b


text:
      PROGRAM CORRIDOR

!***********************************************************************
!* Calculates a Signal Corridor from Y-Series Files
!* and Outputs Average and Cumulative Variance.
!*
!* Language: Fortran
!*
!* Author: Dinesh Sharma
!*         Andrew Orndorff
!*         Shashi Kuppa
!*
!* Date: 2003/07/26
!***********************************************************************

      ! Headers
      INCLUDE 'uds_str.fi'
      INCLUDE 'uds_fxn.fi'
      INCLUDE 'file.fi'

      ! Variables ______________________________________________________

      RECORD /UDS_Y/  S, U
      RECORD /UDS_CONTROL/  C
      RECORD /FILE/  HIT, REPORT, OUTPUT_LIST

      LOGICAL ALL_BATCH, NO_BATCH, NO_INCR, ONE_PASS, NEW_HIT, EOF,
     & BATCH, DONE, DEF_MINIMIZE_RANGE,
     & MINIMIZE_RANGE, WARN_LOW_NUM_POINTS

      INTEGER NUM_ARGS, MAX_ARGS, IOS, IOS_W, I, IE, NUM_FILES,
     & MIN_NFP, MAX_NFP, MIN_NLP, MAX_NLP, NFP_LIMIT, NLP_LIMIT,
     & ANUM_POINTS(NFP_L:NLP_U)

      REAL V_SQ(NFP_L:NLP_U), V(NFP_L:NLP_U), SD(NFP_L:NLP_U),
     & MN(NFP_L:NLP_U), AVERAGE_VARIANCE, CUMULATIVE_VARIANCE

      PARAMETER ( MAX_ARGS = 7 )

      CHARACTER CL_ARG(0:MAX_ARGS)*255, PROG_NAME*8,
     & TEMPLATE_STR*255, DIR_OUT*255,
     & FILE_NAME*255, OUT_FILE*255, DEF_FILE*255,
     & MPS_FILE*255, MMS_FILE*255, OMPS_FILE*255, OMMS_FILE*255,
     & PROMPT_1*79, PROMPT_2*79,
     & XUNITS_1*20, YUNITS_1*20, RANGE_C, MESSAGE*132

      ! Initializations
      PROG_NAME = 'Corridor'
      CALL FILE_INIT( HIT )
      MINIMIZE_RANGE = .FALSE.
      CALL GET_CL( CL_ARG, NUM_ARGS, MAX_ARGS )
      CALL CORRIDOR_CL( CL_ARG, NUM_ARGS, PROG_NAME,
     & DEF_MINIMIZE_RANGE, MINIMIZE_RANGE )
      CALL PROC_CL( CL_ARG, NUM_ARGS, PROG_NAME, TEMPLATE_STR,
     & ALL_BATCH, NO_BATCH, NO_INCR, REPORT, OUTPUT_LIST, DIR_OUT,
     & ONE_PASS, .TRUE., .TRUE. )
      PROMPT_1 = 'First UDS file name/template/hitlist?'
      PROMPT_1(74:) = '[done]'
      PROMPT_2 = 'Next UDS file name/template/hitlist?'
      PROMPT_2(71:) = '[no more]'
      CALL UDS_CONTROL_INIT( C )
      C.CHANFORM = 'Y' ! Only accept Y-series file input

      WRITE( *, '(//23X,A//)' ) '**  Calculate a Signal Corridor **'

      ! Read the first file name/template/hitlist
  100 CALL FILE_PROMPT( TEMPLATE_STR, PROMPT_1, DIR_OUT,
     & ONE_PASS, HIT, NEW_HIT, FILE_NAME, EOF )
      IF ( ( FILE_NAME .EQ. ' ' ) .OR. ( EOF ) ) GO TO 199 ! Done

      ! Batch setup
      CALL BAT_SET( ALL_BATCH, NO_BATCH, NEW_HIT, HIT, BATCH )

      ! Obtain other program settings
  101 IF ( DEF_MINIMIZE_RANGE ) THEN
        MINIMIZE_RANGE = .TRUE.
      ELSE
        WRITE( *, '(/A,25X,A/'' >> '',$)' )
     &     ' Use only data points common to all signals?   (Y/N)',
     &     '[N]'
        READ( *, '(A)', IOSTAT=IOS ) RANGE_C
        IF ( IOS .NE. 0 ) THEN
          CALL MSG_WRITE( '*** Operation skipped', REPORT )
          CALL FILE_DEL( HIT )
          GO TO 190
        END IF
        IF ( ANY_CHARS( RANGE_C, 'Yy' ) ) THEN
          MINIMIZE_RANGE = .TRUE.
        ELSE IF ( ANY_CHARS( RANGE_C, 'Nn ' ) ) THEN
          MINIMIZE_RANGE = .FALSE.
        ELSE
          WRITE( *, * ) '*** Invalid option'
          GO TO 101
        END IF
      END IF

      ! Read first UDS file
      CALL REP_COLS( REPORT )
      C.DIMSYS = ' '
      CALL UDS_GET( FILE_NAME, S, C, REPORT, IOS )
      IF ( IOS .NE. 0 ) THEN ! UDS read error
        IF ( HIT.OPEN ) CALL FILE_DEL( HIT )
        GO TO 190
      END IF

      ! Per-pass initializations
      WARN_LOW_NUM_POINTS = .TRUE.

      AVERAGE_VARIANCE = 0.0
      CUMULATIVE_VARIANCE = 0.0

      ! Process first UDS file
      NFP_LIMIT = 0
      NLP_LIMIT = 0

      MIN_NFP = S.NFP
      MAX_NFP = S.NFP
      MIN_NLP = S.NLP
      MAX_NLP = S.NLP

      XUNITS_1 = S.XUNITS
      YUNITS_1 = S.YUNITS
      C.DIMSYS = S.DIMSYS
      NUM_FILES = 1

      DO I = S.NFP, S.NLP
        ANUM_POINTS(I) = ANUM_POINTS(I) + 1
        V_SQ(I) = S.Y(I)**2
        V(I) = S.Y(I)
      END DO

      ! Read/process the other UDS files
      DONE = .FALSE.
      DO WHILE ( .NOT. DONE )

        ! Read another file name/template/hitlist
  110   CALL FILE_PROMPT( TEMPLATE_STR, PROMPT_2, DIR_OUT,
     &     ( ONE_PASS .OR. HIT.OPEN ), HIT, NEW_HIT, FILE_NAME, EOF )

        IF ( EOF ) THEN
          GO TO 190
        ELSE IF ( FILE_NAME .EQ. ' ' ) THEN ! No more files
          DONE = .TRUE.
        ELSE ! Read/process a file

          ! Read next UDS file
          CALL UDS_GET( FILE_NAME, U, C, REPORT, IOS )
          IF ( IOS .NE. 0 ) THEN ! UDS read error
            IF ( ( HIT.OPEN ) .OR. ( ONE_PASS ) ) THEN ! Skip to next
              CALL FILE_DEL( HIT )
              GO TO 190
            ELSE
              GO TO 110
            END IF
          END IF

          ! Check UDS fields for consistency
          IF ( .NOT. EQUAL_SP( U.DEL, S.DEL, 2.E-7 ) ) THEN
            CALL MSG_WRITE( '*** ERROR - Inconsistent DEL step',
     &         REPORT )
            IF ( ( HIT.OPEN ) .OR. ( ONE_PASS ) ) THEN ! Skip to next
              CALL FILE_DEL( HIT )
              GO TO 190
            ELSE
              GO TO 110
            END IF
          END IF
          IF ( U.YUNITS .NE. YUNITS_1 ) CALL MSG_WRITE(
     &       '*** WARNING - Inconsistent Y-units', REPORT )
          IF ( U.XUNITS .NE. XUNITS_1 ) CALL MSG_WRITE(
     &       '*** WARNING - Inconsistent X-units', REPORT )

          ! Calculate new minimum / maximum NFP and NLP values
          MIN_NFP = MIN(MIN_NFP,U.NFP)
          MAX_NFP = MAX(MAX_NFP,U.NFP)
          MIN_NLP = MIN(MIN_NLP,U.NLP)
          MAX_NLP = MAX(MAX_NLP,U.NLP)

          ! Merge UDS fields
          CALL UDS_MERGE( U, S )

          ! Begin adding signals
          DO I = S.NFP, S.NLP

            ANUM_POINTS(I) = ANUM_POINTS(I) + 1
            IF ( WARN_LOW_NUM_POINTS ) THEN
              IF ( ANUM_POINTS(I) .LT. NUM_FILES ) THEN
                WARN_LOW_NUM_POINTS = .FALSE.
                CALL MSG_WRITE( '*** WARNING - '//
     &             'Number of data points '//
     &             'less than number of curves',
     &             REPORT )
              END IF
            END IF

            S.Y(I) = S.Y(I) + U.Y(I)
            V_SQ(I) = V_SQ(I) + U.Y(I)**2
            V(I) = V(I) + U.Y(I)

          END DO

          ! Update for this file
          NUM_FILES = NUM_FILES + 1
        END IF

      END DO

      ! Compute mean and standard deviation
      IF ( MINIMIZE_RANGE ) THEN
        NFP_LIMIT = MAX_NFP
        NLP_LIMIT = MIN_NLP
      ELSE
        NFP_LIMIT = MIN_NFP
        NLP_LIMIT = MAX_NLP
      END IF

      DO I = NFP_LIMIT, NLP_LIMIT
        S.Y(I) = S.Y(I) / ANUM_POINTS(I)
        MN(I) = S.Y(I)
        SD(I) = SQRT(
     &     ( V_SQ(I) - ( 2 * V(I) * MN(I) ) +
     &     ( ANUM_POINTS(I) * ( MN(I)**2 ) ) ) /
     &     ( ANUM_POINTS(I) - 1 ) )
        CUMULATIVE_VARIANCE = CUMULATIVE_VARIANCE + SD(I)**2
      END DO

      AVERAGE_VARIANCE =
     & CUMULATIVE_VARIANCE / ( NLP_LIMIT - NFP_LIMIT + 1 )

      WRITE( MESSAGE, '(10X,A,F9.5)' )
     &  '*** Average Variance = ', AVERAGE_VARIANCE
      CALL MSG_WRITE( MESSAGE, REPORT )

      WRITE( MESSAGE, '(10X,A,F14.5)' )
     &  '*** Cumulative Variance = ', CUMULATIVE_VARIANCE
      CALL MSG_WRITE( MESSAGE, REPORT )

      WRITE( MESSAGE, '(10X,A,I6)' )
     &  '*** Number of Data Points = ', ( NLP_LIMIT - NFP_LIMIT + 1 )
      CALL MSG_WRITE( MESSAGE, REPORT )

      ! Prepare output fields
      S.NFP = NFP_LIMIT
      S.NLP = NLP_LIMIT
      S.CURNAM = 'AVG'
      S.STATUS = 'COMPUTED'
      S.HIC = 0.
      S.T1 = 0.
      S.T2 = 0.
      S.HICDTUP = 0.
      S.CLIP3M = 0.
      S.CSI = 0.
      S.CD2 = 'MEAN SIGNAL(     )'
      WRITE( S.CD2(13:17), '(I5)', IOSTAT=IOS_W ) NUM_FILES

      ! Set default output file names
      CALL FN_DEF( S, OUT_FILE )
      CALL FN_COMBO( S, OUT_FILE )
      CALL ADD_PATH( DIR_OUT, OUT_FILE )
      IF ( .NOT. NO_INCR ) CALL FN_INCR( OUT_FILE )
      CALL FN_NAME_EXTRACT( OUT_FILE )
      DEF_FILE = OUT_FILE
      IE = FT_POSN_S( DEF_FILE, '.uds' )
      IF ( IE .EQ. 0 ) IE = LEN_TRIM( DEF_FILE ) + 1
      OUT_FILE = DEF_FILE(:IE-1)//'.aps'
      CALL ADD_PATH( DIR_OUT, OUT_FILE )
      IF ( .NOT. NO_INCR ) CALL FN_INCR( OUT_FILE )
      CALL FN_NAME_EXTRACT( OUT_FILE )
      MPS_FILE = OUT_FILE
      OUT_FILE = DEF_FILE(:IE-1)//'.ams'
      CALL ADD_PATH( DIR_OUT, OUT_FILE )
      IF ( .NOT. NO_INCR ) CALL FN_INCR( OUT_FILE )
      CALL FN_NAME_EXTRACT( OUT_FILE )
      MMS_FILE = OUT_FILE

      ! Set the output UDS file names
      IF ( .NOT. BATCH ) THEN
  120   WRITE( *, '(/A/3A/'' >> '',$)' )
     &     ' Output UDS mean file name?',
     &     ' [', DEF_FILE(:L_TRIM(DEF_FILE)), ']'
        READ( *, '(A)', END=190 ) OUT_FILE
        IF ( OUT_FILE .EQ. ' ' ) THEN ! Use default file name
          OUT_FILE = DEF_FILE
        ELSE ! Check file name validity
          CALL ADD_PATH( DIR_OUT, OUT_FILE )
          IF ( .NOT. FN_VALID( OUT_FILE ) ) THEN
            WRITE( *, * ) '*** Invalid file name'
            GO TO 120
          END IF
        END IF
  130   WRITE( *, '(/A/3A/'' >> '',$)' )
     &     ' Output UDS mean+sigma file name?',
     &     ' [', MPS_FILE(:L_TRIM(MPS_FILE)), ']'
        READ( *, '(A)', END=190 ) OMPS_FILE
        IF ( OMPS_FILE .EQ. ' ' ) THEN ! Use default file name
          OMPS_FILE = MPS_FILE
        ELSE ! Check file name validity
          CALL ADD_PATH( DIR_OUT, OMPS_FILE )
          IF ( .NOT. FN_VALID( OMPS_FILE ) ) THEN
            WRITE( *, * ) '*** Invalid file name'
            GO TO 130
          END IF
        END IF
  140   WRITE( *, '(/A/3A/'' >> '',$)' )
     &     ' Output UDS mean-sigma file name?',
     &     ' [', MMS_FILE(:L_TRIM(MMS_FILE)), ']'
        READ( *, '(A)', END=190 ) OMMS_FILE
        IF ( OMMS_FILE .EQ. ' ' ) THEN ! Use default file name
          OMMS_FILE = MMS_FILE
        ELSE ! Check file name validity
          CALL ADD_PATH( DIR_OUT, OMMS_FILE )
          IF ( .NOT. FN_VALID( OMMS_FILE ) ) THEN
            WRITE( *, * ) '*** Invalid file name'
            GO TO 140
          END IF
        END IF
      ELSE ! Batch
        OUT_FILE = DEF_FILE
        OMPS_FILE = MPS_FILE
        OMMS_FILE = MMS_FILE
      END IF

      ! Write the output UDS mean file
      CALL UDS_PUT( OUT_FILE, S, C, REPORT, OUTPUT_LIST,
     & DIR_OUT, .TRUE., IOS )
      IF ( IOS .NE. 0 ) GO TO 190

      ! Compute the mean + sigma signal
      DO I = NFP_LIMIT, NLP_LIMIT
        S.Y(I) = MN(I) + SD(I)
      END DO
      S.CURNAM = 'MPS'
      S.CD2 = 'MEAN + SIGMA(     )'
      WRITE( S.CD2(14:18), '(I5)', IOSTAT=IOS_W ) NUM_FILES

      ! Write the output UDS mean + sigma file
      CALL UDS_PUT( OMPS_FILE, S, C, REPORT, OUTPUT_LIST,
     & DIR_OUT, .TRUE., IOS )
      IF ( IOS .NE. 0 ) GO TO 190

      ! Compute the mean - sigma signal
      DO I = NFP_LIMIT, NLP_LIMIT
        S.Y(I) = MN(I) - SD(I)
      END DO
      S.CURNAM = 'MMS'
      S.CD2 = 'MEAN - SIGMA(     )'
      WRITE( S.CD2(14:18), '(I5)', IOSTAT=IOS_W ) NUM_FILES

      ! Write the output UDS mean - sigma file
      CALL UDS_PUT( OMMS_FILE, S, C, REPORT, OUTPUT_LIST,
     & DIR_OUT, .TRUE., IOS )
      IF ( IOS .NE. 0 ) GO TO 190

      WRITE( *, '(/A,I6,A)', IOSTAT=IOS )
     & ' *** ', NUM_FILES, ' files processed'

      ! Loop for next operation
  190 WRITE( *, '(80(''_''))' )
      GO TO 100

      ! Clean up and stop
  199 CALL CLEAN_UP( HIT, REPORT, OUTPUT_LIST, DIR_OUT )

      END

      SUBROUTINE CORRIDOR_CL( CL_ARG, NUM_ARGS, PROG_NAME,
     & DEF_MINIMIZE_RANGE, MINIMIZE_RANGE )

!***********************************************************************
!* Processes the Program-Specific Command Line Arguments
!***********************************************************************

      ! Arguments ______________________________________________________

      CHARACTER*(*) CL_ARG(0:*) ! Command line argument array

      INTEGER NUM_ARGS ! Number of arguments

      CHARACTER*(*) PROG_NAME ! Program name

      LOGICAL DEF_MINIMIZE_RANGE !

      LOGICAL MINIMIZE_RANGE ! Indicates to use range of common points

      ! Variables ______________________________________________________

      INTEGER I, NS

      CHARACTER ARG*255, C

      ! Process command line arguments
      DEF_MINIMIZE_RANGE = .FALSE.
      MINIMIZE_RANGE = .FALSE.
      DO I = 1, NUM_ARGS
        CALL STR_UPCASE( ARG, CL_ARG(I) )
        C = ARG(1:1)
        IF ( ( C .EQ. '-' ) .OR. ( C .EQ. '/' ) ) THEN
          NS = 2
        ELSE
          NS = 1
        END IF
        IF ( ( ARG(NS:) .EQ. '?' ) .OR.
     &     ( ARG(NS:) .EQ. '??' ) .OR.
     &     ( ARG(NS:) .EQ. 'HELP' ) ) THEN ! Show command line syntax
          CALL SYNTAX_CL( PROG_NAME )
          WRITE( *, '(10X,A)' ) '[MINIMIZE_RANGE]'
          IF ( ( ARG(NS:) .EQ. '??' ) .OR.
     &       ( ARG(NS:) .EQ. 'HELP' ) ) THEN ! Show command line help
            CALL HELP_CL( PROG_NAME )
            WRITE( *, '(10X,A)' )
     &         'MINIMIZE_RANGE : '//
     &         'Use only those points in common for calculations'
          END IF
          STOP ' '
        ELSE IF ( ARG(NS:) .EQ. 'MINIMIZE_RANGE' ) THEN ! Default
          DEF_MINIMIZE_RANGE = .TRUE.
          MINIMIZE_RANGE = .TRUE.
          CL_ARG(I) = ' '
        END IF
      END DO

      RETURN
      END

ext: for

max_stars_repo_name: DeadParrot/NHTSA-Tools

max_issues_repo_name: DeadParrot/NHTSA-Tools

max_forks_repo_name: DeadParrot/NHTSA-Tools


text:
c
c --------------------------------------------------------------------
c
       subroutine icall(val,aux,nrow,ncol,nfil,nvar,naux,
     .                  ilo,ihi,jlo,jhi,klo,khi,level,
     .                  iputst,jputst,kputst)

      use amr_module
      implicit double precision (a-h, o-z)

      dimension val(nvar,nrow,ncol,nfil)
      dimension aux(naux,nrow,ncol,nfil)

       iadd   (ivar,i,j,k) = loc    +    (ivar-1)
     &                              +    (i-1)*nvar
     &                              +    (j-1)*nvar*mitot
     &                              +    (k-1)*nvar*mitot*mjtot
       iaddaux(ivar,i,j,k) = locaux +    (ivar-1)
     &                              +    (i-1)*naux
     &                              +    (j-1)*naux*mitot
     &                              +    (k-1)*naux*mitot*mjtot

c ::::::::::::::::::::::::::: ICALL :::::::::::::::::::::::::::::::
c
c    find intersecting grids at the same level. copy data from
c    intersecting grids to both val and aux arrays.
c
c    use larger definition of grids here - boundary data already in.
c    aux arrays also enlarged size.
c
c    iputst, jputst, kputst: where to copy values into. may not be in
c                            location corresponding to ilo,ihi,etc. if
c                            the patch has been periodically wrapped.

c
c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

       mptr = lstart(level)

 10    if (mptr .eq. 0) go to 99
         iglo = node(ndilo,mptr)
         ighi = node(ndihi,mptr)
         jglo = node(ndjlo,mptr)
         jghi = node(ndjhi,mptr)
         kglo = node(ndklo,mptr)
         kghi = node(ndkhi,mptr)

c         # does it intersect?
!--         ixlo = max(iglo-nghost,ilo)
!--         ixhi = min(ighi+nghost,ihi)
!--         jxlo = max(jglo-nghost,jlo)
!--         jxhi = min(jghi+nghost,jhi)
!--         kxlo = max(kglo-nghost,klo)
!--         kxhi = min(kghi+nghost,khi)
c  how did ghost cells get in the allowable region? They are not filled
c  (since we may be interpolating from newly filled grids, not just grids
c  that have been primed with bcs to be advanced.
         ixlo = max(iglo,ilo)
         ixhi = min(ighi,ihi)
         jxlo = max(jglo,jlo)
         jxhi = min(jghi,jhi)
         kxlo = max(kglo,klo)
         kxhi = min(kghi,khi)

         if ((ixlo .le. ixhi .and. jxlo .le. jxhi) .and.
     &       (                     kxlo .le. kxhi)) then
           loc  = node(store1,mptr)
           locaux = node(storeaux,mptr)
           nx   = ighi - iglo + 1
           ny   = jghi - jglo + 1
           nz   = kghi - kglo + 1
           mitot = nx + 2*nghost
           mjtot = ny + 2*nghost
           mktot = nz + 2*nghost
           do 30 k    = kxlo, kxhi
           do 30 j    = jxlo, jxhi
           do 30 i    = ixlo, ixhi
           do 20 ivar = 1, nvar
              ialloc  =  iadd(ivar,i-iglo+nghost+1,j-jglo+nghost+1,
     &                                             k-kglo+nghost+1)
              val(ivar,i-ilo+iputst,j-jlo+jputst,
     &                              k-klo+kputst) = alloc(ialloc)
 20     continue
           do 25 iaux = 1, naux
             ialloc = iaddaux(iaux,i-iglo+nghost+1,j-jglo+nghost+1,
     &                                             k-kglo+nghost+1)
             aux(iaux,i-ilo+iputst,j-jlo+jputst,
     &                             k-klo+kputst) = alloc(ialloc)
 25        continue
 30       continue
       endif
       mptr = node(levelptr, mptr)
       go to 10

 99   return
      end

ext: f

max_stars_repo_name: navravi/amrclaw

max_issues_repo_name: navravi/amrclaw

max_forks_repo_name: BrisaDavis/amrclaw


text: #redirect Users/LourdesEBorrayo

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
program stack
   use constants
   use io
   use riff
   implicit none

   integer :: i, n, t
   integer(i2) :: c
   real(dp), allocatable :: sound(:, :)
   type(audio), allocatable :: p(:)
   type(audio) :: s

   s%channels  = 1
   s%points    = 0
   s%rate      = 1.0_dp
   s%amplitude = 0.0_dp

   n = command_argument_count()

   allocate(p(n - 1))

   do i = 1, n - 1
      call read_riff(command_argument(i, '/dev/stdin'), p(i))

      s%channels = max(s%channels, p(i)%channels)
      s%points   = max(s%points,   p(i)%points)
      s%rate     = max(s%rate,     p(i)%rate)
   end do

   allocate(s%sound(s%channels, s%points))
   allocate(  sound(s%channels, s%points))

   sound(:, :) = 0.0_dp

   do i = 1, n - 1
      do c = 1, s%channels
         do t = 1, s%points
            sound(c, t) = sound(c, t) + p(i)%amplitude * p(i)%sound( &
               1_i2 + modulo(c - 1_i2, p(i)%channels), &
               1_i2 + modulo(t - 1_i2, p(i)%points))
         end do
      end do
   end do

   sound = sound / i2max

   if (s%points .gt. 0) then
      s%amplitude = maxval(abs(sound))
      s%sound = nint(i2max / s%amplitude * sound, i2)
   end if

   call write_riff(command_argument(-1, '/dev/stdout'), s)
end program stack

max_stars_repo_name: janberges/Tonbandfetzen

max_issues_repo_name: janberges/Tonbandfetzen

max_forks_repo_name: janberges/Tonbandfetzen


text:
*DECK FAC
      FUNCTION FAC (N)
C***BEGIN PROLOGUE  FAC
C***PURPOSE  Compute the factorial function.
C***LIBRARY   SLATEC (FNLIB)
C***CATEGORY  C1
C***TYPE      SINGLE PRECISION (FAC-S, DFAC-D)
C***KEYWORDS  FACTORIAL, FNLIB, SPECIAL FUNCTIONS
C***AUTHOR  Fullerton, W., (LANL)
C***DESCRIPTION
C
C FAC(N) evaluates the factorial function of N.  FAC is single
C precision.  N must be an integer between 0 and 25 inclusive.
C
C***REFERENCES  (NONE)
C***ROUTINES CALLED  GAMLIM, R9LGMC, XERMSG
C***REVISION HISTORY  (YYMMDD)
C   770601  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890531  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
C***END PROLOGUE  FAC
      DIMENSION FACN(26)
      SAVE FACN, SQ2PIL, NMAX
      DATA FACN( 1) / 1.0E0 /
      DATA FACN( 2) / 1.0E0 /
      DATA FACN( 3) / 2.0E0 /
      DATA FACN( 4) / 6.0E0 /
      DATA FACN( 5) / 24.0E0 /
      DATA FACN( 6) / 120.0E0 /
      DATA FACN( 7) / 720.0E0 /
      DATA FACN( 8) / 5040.0E0 /
      DATA FACN( 9) / 40320.0E0 /
      DATA FACN(10) / 362880.0E0 /
      DATA FACN(11) / 3628800.0E0 /
      DATA FACN(12) / 39916800.0E0 /
      DATA FACN(13) / 479001600.0E0 /
      DATA FACN(14) / 6227020800.0E0 /
      DATA FACN(15) / 87178291200.0E0 /
      DATA FACN(16) / 1307674368000.0E0 /
      DATA FACN(17) / 20922789888000.0E0 /
      DATA FACN(18) / 355687428096000.0E0 /
      DATA FACN(19) / 6402373705728000.0E0 /
      DATA FACN(20) /  .12164510040883200E18 /
      DATA FACN(21) /  .24329020081766400E19 /
      DATA FACN(22) /  .51090942171709440E20 /
      DATA FACN(23) /  .11240007277776077E22 /
      DATA FACN(24) /  .25852016738884977E23 /
      DATA FACN(25) /  .62044840173323944E24 /
      DATA FACN(26) /  .15511210043330986E26 /
      DATA SQ2PIL / 0.9189385332 0467274E0/
      DATA NMAX / 0 /
C***FIRST EXECUTABLE STATEMENT  FAC
      IF (NMAX.NE.0) GO TO 10
      CALL GAMLIM (XMIN, XMAX)
      NMAX = XMAX - 1.
C
 10   IF (N .LT. 0) CALL XERMSG ('SLATEC', 'FAC',
     +   'FACTORIAL OF NEGATIVE INTEGER UNDEFINED', 1, 2)
C
      IF (N.LE.25) FAC = FACN(N+1)
      IF (N.LE.25) RETURN
C
      IF (N .GT. NMAX) CALL XERMSG ('SLATEC', 'FAC',
     +   'N SO BIG FACTORIAL(N) OVERFLOWS', 2, 2)
C
      X = N + 1
      FAC = EXP ( (X-0.5)*LOG(X) - X + SQ2PIL + R9LGMC(X) )
C
      RETURN
      END


text:
*DECK PYTHAG
      REAL FUNCTION PYTHAG (A, B)
C***BEGIN PROLOGUE  PYTHAG
C***SUBSIDIARY
C***PURPOSE  Compute the complex square root of a complex number without
C            destructive overflow or underflow.
C***LIBRARY   SLATEC
C***TYPE      SINGLE PRECISION (PYTHAG-S)
C***AUTHOR  (UNKNOWN)
C***DESCRIPTION
C
C     Finds sqrt(A**2+B**2) without overflow or destructive underflow
C
C***SEE ALSO  EISDOC
C***ROUTINES CALLED  (NONE)
C***REVISION HISTORY  (YYMMDD)
C   811101  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900402  Added TYPE section.  (WRB)
C***END PROLOGUE  PYTHAG
      REAL A,B
C
      REAL P,Q,R,S,T
C***FIRST EXECUTABLE STATEMENT  PYTHAG
      P = MAX(ABS(A),ABS(B))
      Q = MIN(ABS(A),ABS(B))
      IF (Q .EQ. 0.0E0) GO TO 20
   10 CONTINUE
         R = (Q/P)**2
         T = 4.0E0 + R
         IF (T .EQ. 4.0E0) GO TO 20
         S = R/T
         P = P + 2.0E0*P*S
         Q = Q*S
      GO TO 10
   20 PYTHAG = P
      RETURN
      END

max_stars_repo_name: homiak/jabaws

max_issues_repo_name: homiak/jabaws

max_forks_repo_name: homiak/jabaws


module correlation

! JLR 30/6/2014
!
! uneven sampled correlation using Gaussian kernel using method of
! Rehfeld et al, Nonlinear Processes in Geophysics, 2011

implicit none

contains

function correlate_gaussian(x,y,tx,ty)
real :: correlate_gaussian

real, intent(in) :: x(:),y(:),tx(:),ty(:)
integer :: n,m,i,j
real :: b,h,d,xmean,ymean,txmean,tymean,delta_t,pi
real :: sdx(size(x)),sdy(size(x)),num(size(x))
real :: sdxg,sdyg,numg,deng

n=size(x,1)
m=size(y,1)

pi=abs(atan2(0.0,-1.0))

! calc mean for x series
xmean=0.0
do i=1,n
  xmean=xmean+x(i)
enddo
xmean=xmean/n
txmean=(tx(n)-tx(1))/(n-1)

! calc mean for y series
ymean=0.0
do j=1,m
  ymean=ymean+y(j)
enddo
ymean=ymean/m
tymean=(ty(m)-ty(1))/(m-1)

sdx=0.0
sdy=0.0
num=0.0
delta_t=max(txmean,tymean)
h=delta_t/4

!$acc kernels loop gang, vector(64)
do i=1,n
  do j=1,m
    d=ty(j)-tx(i)
    b=exp(-d**2/(2*h**2))/sqrt(2*pi*h)
    num(i)=num(i)+(x(i)-xmean)*(y(j)-ymean)*b
    sdx(i)=sdx(i)+b*(x(i)-xmean)**2
    sdy(i)=sdy(i)+b*(y(j)-ymean)**2
  enddo
enddo
!$acc end kernels

sdxg=0.0
sdyg=0.0
numg=0.0

do i=1,n
  sdxg=sdxg+sdx(i)
  sdyg=sdyg+sdy(i)
  numg=numg+num(i)
enddo

correlate_gaussian=numg/sqrt(sdxg*sdyg)

end function

end module


text:
Although the UC system is considered the pinnacle of Californias postsecondary education system, some feel the practical education available to undergraduate students undergraduates is often not as complete as that found at http://www.calstate.edu/ California State Universities, or even many community colleges.

These individuals believe there are a number of reasons for this, but assert that the basic reason is that the UC system focuses very strongly on research and graduate education.  Then again, maybe the fact that the UCs emphasize research and graduate education is a strength rather than a weakness  read below for pros and cons.

In the UCs, there are a number of opportunities for Student Influence.

I originally created this page out of frustration with the way I felt many undergrads were being shortchanged in favor of research and other priorities (disclosure: when I wrote this, I was NOT an undergrad, although I had been in the past, and I was a TA working with undergrads in the past as well). I think however, that at this point a more balanced page would be more helpful (along with a name change). Any objections before I rewrite?  Users/EricKlein

 The benefit of this focus is that the UC system has many of the best researchers in their fields and provides some of the best graduate education to be found anywhere on Earth. The downside is that a lot of choices that the UC makes to optimize research benefit wind up damaging undergraduate education. The way I think of it, the UC system is basically a prep school for grad school, while the CSU system is a prep school for the business/real world.

 Some examples of this phenomenon:
   Sometimes a professor who may be one of the very best researchers in their field may also be a very poor teacher (I had a few of these as an undergrad at UC Santa Barbara).
   The undergrad classes in the UC system usually seem to focus much more on theory than on practice. The upside of this is that many undergrads at UC schools have very good theoretical underpinnings when they graduate, which can serve them well, especially if they plan to attend grad school. The downside is that they are abysmal at the practice of their trade. Im most familiar with how this relates to computer science. I was fortunate in that I worked my way through school by programming, however many of my classmates with less practical experience during college had horrible coding skills, and were not all that in demand by local employers as a result. In contrast students attending Cal Poly SLO,  or even from the local Santa Barbara City College often had better coding skills due to the better teaching at those schools.
    It sounds like you were expecting a software engineering degree instead of a computer science one.

 I dont mean to speak ill of the UC system, because it really is an excellent system (especially for grad school). I feel like they could best serve the undergrads though by doing one of two things:
  1. Being upfront about their priorities, tell students that if they want a practical education, their best bet would be a CSU. Also let the business world know this as well, so that there is no prestige loss with having a CSU degree as opposed to a UC degree. Let everyone know that the UC is all about research, and not making any apologies about it.
  1. Hire more lecturers who are excellent teachers. Let the professors who want to just do research, just do research. Make sure all TAs have not only taken training, but also have passed some basic tests of teaching ability (and communications skills). Make sure those same TAs are given as much support as possible to be good teachers, and to improve. Also, add more internship/coop opportunities with industry, so that students can get more practical/handson experience in their field. Make sure students not only know about these opportunities, but ideally require undergrads to take advantage of them (or very, very strongly encourage it). Emphasize the hiring of experienced people from industry as lecturers (or even professors). It is much easier to learn practical job skills from someone who actually has them.

Ive noticed that few davis students actually go onto to graduate school so for the majority of students its not really a prep for grad school. Most students attend UCs instead of CSUs because they have better reputations. Users/BryanBell

Do you have proof/statistics/studies that prove the claims made on this page?  Users/SummerSong ss
  Its a hard thing to quantify. Who would pay for such a study? Undergrads?  Users/ArlenAbraham arlen
  I absolutely do NOT have hard proof. The statements above were based purely on anecdotal evidence from my own experiences, and the experiences other have told me about. I have no problem with somebody rewriting this to show that it is opinion and not fact. I dont know of any formal studies showing this, nor do I expect any to be funded (although it would be quite interesting). Users/EricKlein ek

I have generally heard the most effective use (financially and educationally) of Californias higher education is 2 years at a JC, Bachelors from a State, and Masters/PhD from a UC.  Davis is a great research university, but with the publishing requirements, a fair portion of the faculty are more concerned about their research, and less about teaching, I have have several professors who were obviously brilliant, but lousy teachers. Users/RogerClark
   I dont necessarily disagree with the sentiments expressed on this page, but I have a hard time accepting conjecture as fact.  Maybe the page should be reworded and presented as an opinion?  Users/SummerSong ss
   I dont think the burdens on a UC professor make them weaker professors.  In fact, UC professors generally have a far lighter teaching load than Community College professors.  And they often teach fewer classes than Cal State professors.  So on the one hand, they could very easily have a great deal more ability to really focus on the classes they teach.  On the other hand, if youre teaching the same course over and over again, its not going to require much preparation anyway.  So it really comes down to who the people are, in my view.  I think the problem (to the extent there is one) is that too much responsibility falls on TAs and also the criteria for which the great intellectuals who are hired at UC Davis more emphasizes research accomplishments, while teaching ability is seen more as an added bonus.  But honestly, hiring the most brilliant people you can isnt a bad way to go.  Some of them may be terrible teachers (and great researchers), but I think the vast majority of them outshine most human beings in both teaching and acamdeic research.  Users/JaimeRaba

20050923 22:39:54 nbsp Hm, interesting to see this discussion, as a professor.  The problem, from our end, is class size.  Honestly.  Its incredibly hard to teach a large class unless its solely factbased.  No matter how good a teacher you wish you were, and how hard you prepare, you know you are losing the top 5% and the bottom 20% of students.  You want to know them, individualize the learning a bit, and you end up just trying to figure out how to get papers back without wasting class time.    You end up relying on memorization/regurgitation way more than you want to, because other teaching techniques are so difficult to implement with lots of students.  Then UC students are inculcated into the lecture system such that other methods can really frustrate them.  On top of all that, research is indeed what gets you promoted, so you cant skimp on that either.  Not offering excuses for bad teaching, here, just want you to know some of our obstacles.   Users/BethFreeman
  Thanks for offering the other side of things in a calm and respectful way. Users/KenjiYamada

20051215 15:10:25 nbsp Ive missed something  is UC supposed to be a highend vocational school?  I thought UC was supposed to provide education, not job training.  A lot of the money for UC is provided by the public, and the public should get something in return.  Traditionally, the argument has been that the public gets in return a body of citizens able to think more clearly and responsibly, especially about public issues.  This is not to say that UC is doing well at providing education even in the sense I have described.  But if you are looking only for job skills which will benefit no one but you and your future employer, then it is unrealistic to want the state (or, in the case of private universities  generous donors) to pay for so much of it. Users/AlexanderWoo
  So... what do you mean by education if not job training? Isnt job training in the public good? Doesnt the public have an interest in training people to fill needed positions in society? I realize that there are other aspects to education aside from job training (such as creating a more well rounded person), but I dont see why that should require that the job training aspects should be neglected. Im not saying that the UC should become a glorified vocational school, but I do think that a college has an obligation to its students to see that they are employable after getting a degree.  Users/EricKlein
   Well, heres the thing. UCs are highereducation. They are! Theres no two ways about it. Traditionally in many countries, schools are set up the same way. Theres the university, which is devoted to education. And then there are vocational schools devoted to more realworld practical experience. I dont think that the UC system is lying about its intentions by offering us an education based on learning about the subjects we ask to learn about. Now, mind you, I think the concern raised here is a tiny bit valid in that I get a sense after all my classes of being pushed toward academic research as a career; which in my case doesnt seem that practical. But I am not leaving this school; I consider it more my fault for not thinking more clearly about the fact that this is a research university. Beyond that I dont think that this school doesnt offer real world experience. More than 2/3 of my friends keep part time jobs and go to school, and just going to Davis can be a responsibility by itself! You learn a lot about timemanagement, improvement of performance, how much performance matters (grades on your record), how to keep a cool head in many situations ... the list goes on and on for me. The difference in the real world (what little Ive experienced it) is that the deadlines are more concrete and there is a bottom line. Unfortunately, UCs dont attempt to teach you that. However, to say that were not employable after a degree is rather bullcrap. The job market has consistently advocated a university degree over many other kinds for the best jobs offered. SS
   As an example of what I mean by education, I refer to an issue tangentially raised by the http://daviswiki.org/IKEA/Talk?actionrecall&version46 IKEA page debate.  I would want anyone with a degree from UC Davis to either (a) know that it was fairly common for German teenagers to join the Hitler Youth, or (b) be responsible and clearheaded enough, and have enough general knowledge about human society, to consider the possibility that being in the Hitler Youth was common and fairly benign and check to see if this was the case before taking any serious action (such as advocating boycotting a store) based purely on someone having been in the Hitler Youth. Users/AlexanderWoo
    Im so remarkably happy for you. My intention of writing that was not to have people boycott IKEA. Ask yourself a question: how many people, oh wait, dare I say, university students, actually knew about Kamprads involvement at all? If you think the number is considerably large in any way, youre delusional. My point is that people who are educated should more readily be educating themselves (I dont think they do as much as would appease me). Thats the main point in addressing something like that about Kamprad. Whether or not people boycott IKEA based on that information is their own choice. Maybe that information isnt enough for them. Maybe reading that, theyll be prompted to find out more. That would be great! SS

20051216 04:17:16 nbsp I, for one, am disgusted by the number of courses that I have taken as an undergrad in which the instructor had such a thick foreign accent (or equivalently bad vocabulary) that his/her words became gnarled to the point of unintelligibility.  I wouldnt say that these instructors have comprised the majority of the instructors in all the courses Ive taken, but they have comprised the majority in some of the most vitally important courses of my education, which is a complete shame.

The UCs intensive research focus aside, they should hire instructors that have a full mastery of the English language, or, permit only people that have mastery of those skills to teach, so long as they have some minimum qualification of understanding the material.

That is, if the goal is education (it should be).

Regarding the UCs focus more on research rather than undergrad care, I have no problem with that, so long as the university clearly sends that message to applicants (they didnt, to me). Users/JohnNapier

20051216 15:18:00 nbsp I had one guy tell us that he wasnt being paid enough to teach us and didnt want to be there.  He said we should drop his class and sign up for the other guys section because he gets paid more and should have to deal with us instead of him. Users/MarieHuynh

20070503 20:00:49 nbsp I cant really respond to the above, but I can fill in things from my own perspective.

The University of California (UC) system is designed to be research oriented; the California State University (CSU) is designed to put teaching before research.  Does that mean youll get a better education at a CSU?  Maybe, maybe not.  In part, it depends on what youre looking for, but also, experiences will vary wildly from school to school, from department to department, and from professor to professor.  Does it matter to you whether your professors are top in their field, actively doing research and staying current?  Then perhaps the UC is a better place for you  but plenty of CSU professors are good researchers, too.  Then again, it might be a mistake to think that CSU professors have tons of time to devote to students  they have very heavy teaching loads, and are expected to do some research, regardless.  Bottom line: the truth is that most faculty become professors because they want to research  grad school is a long, expensive haul for those who want to teach.  Listen to the recommendations of your friends, and choose accordingly.  (That is, choose the professors who obviously care about their teaching  if you choose classes only because they are easy or because the schedule is convenient, then you get no sympathy from me).

Personally, I bucked the norm and went to grad school so that I could teach, because I wanted to teach and couldnt see myself teaching high school.  By the time that the ink was dry, I found that I loved both teaching and research.  But it can be hard to balance both.  Sometimes I think that students dont have a complete understanding of all the demands on a professors time.  (Sometimes professors dont have a complete understanding of all the demands on a students time, either, but thats a story for the students to tell).  So, heres a partial list:

  Research responsibilities:
   Experiments/literature review
   Writing papers and books
   Presenting papers at conferences and colloquia (and the preparation that entails)
   Refereeing papers for journals (academia runs on peer review, but someone has to do it)
   Writing book reviews
   Communicating with colleagues in the field, giving feedback on their papers
   Staying current in the field (reading books and journals, attending colloquia and conferences)

  Teaching responsibilities:
   Doing the assigned reading (I dont have it memorized!)
   Preparing lecture notes
   Lecturing
   Holding office hours
   Holding virtual office hours (email, etc.)  what used to be confined to weekdays has now spread to evenings and weekends
   Grading exams and papers
   Coordinating with TAs
   Supervising graduate dissertations, serving on diss. committees

  Committee responsibilities:
   Department committees
   Collegelevel committees
   Universitylevel committees
   Other positions such as department chair (good luck getting research done), graduate advisor, etc.

Being a professor is a very fulfilling job.  In fact, it feels funny to call it a job, because it ends up permeating every aspect of your life.  But one basically has to accept that there is an infinite amount of work to do and no one to set any boundaries but oneself.  Professors work insane numbers of hours with the benefits being that they can choose those hours and what they work on (for the most part), not to mention the excitement of seeing young minds think and explore.  So, yes, some professors are neglectful of undergraduates, some of them are at UC Davis, and thats a problem.  However, also understand that professors have a lot of obligations to fulfill, and occasionally, we try to have lives, too.

Users/CovertProfessor

  Very good points. Back when I originally wrote this page, I was quite undiplomatic in my naming of the page. The result is something that is really unfair to most of the professors at the UC, and somewhat insulting. My beef was never really with the professors themselves, but with the way the UC system is presented to the public as a whole, and more specifically the way it is presented to prospective undergrads and to a lesser degree, to potential employers. I think the UC system is extremely valuable, and would not want to see it diminished. I very much value my education both at UCSB and here at UCD. My only real point in creating this page was that I felt that quite a few students came to the UC expecting a practical education, the result of which would be real applicable skills that they could use to get a job, and were left with something other than what they expected. Similarly, the UC is promoted as being the pinnacle of Californias education system (and in many ways it is). Many employers therefore view a degree from a UC as somehow more valuable than one from a CSU. In both of these cases, the way the UC is put forth and promoted does a disservice to the public (and to some extent to the UC itself). I dont view an education from a UC as inferior to one from a CSU, but I do view it as substantially different, and in general not as useful for learning skills that can be applied to a career outside of academia. Im not even suggesting that should change, only that the UC should reconsider the way it markets itself. If it doesnt want to remarket itself, it would be nice to see it try to cover the gap between students (and employers) expectations, and reality. Again, I dont believe that this is necessarily something for individual professors themselves to be involved in. More likely this would be something done at the administrative level: either to change marketing, or to initiate programs to ensure better practical skills for those who want them (more accessible coop programs w/ industry, more nutsandbolts sorts of classes offered, etc.). Users/EricKlein ek

   Well, I wasnt insulted by it; I certainly do think that some professors, at the UC, the CSU, and elsewhere, neglect their undergraduates.  I would be surprised to learn that it was any worse at UCD, but I havent been here long enough to know that.  And I cant comment on how UCD markets itself, not having looked at much of the materials.  But your comment raises two questions: 1) Is the job training the primary purpose of a college education?  and 2) To the extent that job training is a goal of college education (even if it isnt the primary goal), is this goal best achieved through skills training?  Is it best for employers or for students?  With regard to question #1, we could obviously get into a long discussion about the purpose of higher education, so I will just claim without argument that the following are at least as important, if not more important: becoming a person who is broadly educated in arts, letters, and sciences, knowledgeable about the world, its history, and its cultures; becoming a person who is in the habit of thinking and questioning, a person who can argue logically and coherently as well evaluate the arguments of others; becoming a person who is excited about learning, so that college is only the beginning of the learning process, and not the end; becoming a person who is prepared to engage in our democracy and our society as an informed and thoughtful citizen.  I do believe that the CSU also has these goals, although how much they are stressed may differ from campus to campus and within campuses as well.  With regard to #2, with economies and technologies changing so fast, it might not actually be in students best interest to graduate trained in a specific skill.  They might have an easier time getting a job upon graduation, but if they have not learned to think they will have a tough time once those particular skills are out of date.  It is true that some employers look for specific skills, rather than smart people who are broadly educated and who can think and learn.  I happen to think that those employers are shortsighted, and they would do much better to hire the aforementioned smart people.  So, if there is anything that needs changing, in my view it is the employers and not the university, and those of us who are connected to academia in any way would do well to make that case whenever and wherever we can.  Users/CovertProfessor

20100422 12:08:36 nbsp I was a community college student, from my experience the majority Professors Graduate from UC have more difficult skills in teaching vs. Professors Graduate from CSU plain and simple. I learn well instructors from CSU! UC may be good with research, but if CSU allow to do research/medicine, give out more doctorate, I see CSU probably come out superior. Users/news1001


 !***********************************************************************************
    SUBROUTINE SortTemplate (imult)
    !c----------------------------------------------------------------------
    !c
    !c----------------------------------------------------------------------
    use varmod
    implicit none

    ! Declare dummy arguments
    INTEGER, INTENT(IN) :: imult

    ! Declare local variables
    INTEGER :: ind,n
    REAL :: xx,yy,zz,hsqd,cont
    REAL, DIMENSION(MAXNOD) :: tmp
    real, DIMENSION(MAXNOD) :: e,f,g,h

    CALL SetRotMat(imult)

    DO ind=1,nltemplate
        xx=ixtemplate(ind)
        yy=iytemplate(ind)
        zz=iztemplate(ind)
        hsqd = 0.0
        DO n=1,3
            cont=rotmat(n,1)*xx+rotmat(n,2)*yy+rotmat(n,3)*zz
            hsqd = hsqd + cont*cont
        END DO
        tmp(ind) = hsqd
    END DO

    call sortem(1,nltemplate,tmp,3,ixtemplate,iytemplate,iztemplate,e,f,g,h)

    end subroutine SortTemplate
    !***********************************************************************************


subroutine compute_tensor(ifile)
!***********************************************************************
! Subroutine handles the computation of both tensors
!***********************************************************************

use tachedata
implicit none

integer :: ifile
real :: percentcount
real,dimension(3,3) :: tensor_element

print*, '-------------'
! If this is an SPH file, require neighbour lists to compute derivatives
if(filetype=='SPH') then

   call get_SPH_neighbours(ifile)

   allocate(tensor(3,3,nelement))

   if (tensorchoice=='tidal') then
      print*, "Calculating Tidal Tensor"
   else if(tensorchoice=='velocity') then
      print*, "Calculating Velocity Shear Tensor"
   endif
   percentcount=0.0

   ! Loop over particles
   do ielement = 1,nelement

      call element_percent_complete(ielement,nelement,percentcount,10.0)

      tensor_element(:,:) =0.0

      ! Choose which tensor to compute
      if(tensorchoice=='velocity') then
         call calc_velocityshear_tensor(ielement,tensor_element)
      else if(tensorchoice=='tidal') then
         call calc_tidal_tensor(ielement,tensor_element)
      endif

      tensor(:,:,ielement) = tensor_element(:,:)

   enddo
   ! End of loop over particles

endif
print*, "Tensors calculated"
print*,'-------------------'

end subroutine compute_tensor


text:
SUBROUTINE clawpack5_setaux(mbc,mx,my, &
     xlower,ylower,dx,dy,maux,aux)
  IMPLICIT NONE

  INTEGER mbc,mx,my,maux
  DOUBLE PRECISION xlower,ylower,dx,dy

  DOUBLE PRECISION aux(maux,1-mbc:mx+mbc,1-mbc:my+mbc)

  DOUBLE PRECISION rhol, cl, rhor, cr
  COMMON /comaux/ rhol,cl,rhor,cr

  INTEGER i,j
  DOUBLE PRECISION xcell

  DO j=1-mbc,my+mbc
     DO i=1-mbc,mx+mbc
        xcell = xlower + (i-0.5d0)*dx
        IF (xcell .LT. 0.5d0) THEN
           aux(1,i,j) = rhol
           aux(2,i,j) = cl
        ELSE
           aux(1,i,j) = rhor
           aux(2,i,j) = cr
        ENDIF
     ENDDO
  ENDDO

  RETURN
END SUBROUTINE clawpack5_setaux

max_stars_repo_name: ECLAIRWaveS/ForestClaw

max_issues_repo_name: ECLAIRWaveS/ForestClaw

max_forks_repo_name: ECLAIRWaveS/ForestClaw


PROGRAM NFLUC

!  I follow the instruction of James' code! More detail you can see:
!  XDR Fortran Interface XTC Example Program with Wrapper
!  2014 (c) James W. Barnett <jbarnet4@tulane.edu>
!  https://github.com/wesbarnett/

! Summary: Concentration (or density) Fluctuation (dela_N) in slabs

! Use the xdr interface
  USE xtc_mod, only: xtcfile

  IMPLICIT NONE
! Declare a variable of type xtcfile
  type(xtcfile) :: xtc

! command
  INTEGER :: NUM_OF_ARGS, I_ARG, N_SKIP
  CHARACTER(LEN=20) :: FILE_NAME,INPUT_NAME,OUTPUT_NAME1,OUTPUT_NAME2,MOL_NAME
  LOGICAL :: LookForInput=.FALSE.,LookForOutput1=.FALSE.,LookForOutput2=.FALSE.
  LOGICAL :: LookForMOL=.FALSE.,LookForDZ=.FALSE.,LookForBLOCK=.FALSE.
  LOGICAL :: FileExist

  REAL :: D_Z
  INTEGER :: SIZE_BLOCK

! Read
  REAL, DIMENSION(:,:,:), ALLOCATABLE :: AR_DATA, TEMP
  INTEGER :: MAX_DATA = 5000, ICR = 2000, N_DATA, N_ATOM, IERR
  REAL, DIMENSION(:,:), ALLOCATABLE :: BOX, TEMP_BOX

! Calculate density fluctuation
  DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: MASS, NPAR_AVG, NPAR_STD2
  REAL, DIMENSION(:), ALLOCATABLE :: T_SLICE
  DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: NPAR
  DOUBLE PRECISION, DIMENSION(:,:,:), ALLOCATABLE :: COM
  DOUBLE PRECISION :: PI = 3.141592653589793
  INTEGER :: N_INTRA, N_SLICE, I_COM
  INTEGER :: I, J, K, L, INDEX
  REAL :: T_SLICE_AVG, T_SLICE_STD

! Set default
  INPUT_NAME = "trajout.xtc"
  OUTPUT_NAME1 = "N-slab.dat"
  OUTPUT_NAME2 = "N-fluc.dat"
  D_Z = 0.12D0
  SIZE_BLOCK = 5
  MOL_NAME = "bf4"

! Command
  NUM_OF_ARGS = command_argument_count()

  IF (NUM_OF_ARGS == 0) THEN
    WRITE(*,*) "USAGE : ./N-fluc.x -i [xtc file] -mol [bf4/pf6] -o1 [file] -o2 [file] -dz [real] -block [integer]"
    WRITE(*,*) "-i trajout.xtc  xtc input file"
    WRITE(*,*) "-o1 N-slab.dat  # of particles in the slab, N(i,t), output file"
    WRITE(*,*) "-o2 N-fluc.dat  Concentation Fluctuation output file"
    WRITE(*,*) "-mol bf4        Which molecules position you saved in the xtc input file (bf4/pf6)"
    WRITE(*,*) "-dz 0.12        Slab thickness (usually cut-off electrostatic interaction is 1.5)"
    WRITE(*,*) "-block 5        size of block (5 = every 5 steps (time) average)"
    WRITE(*,*) "It calcualtes concentration (or denstiy) average and fluctuation in slabs"
    WRITE(*,*) "xtc file should be created by g_traj -jump option"
    WRITE(*,*) "EXAMPLE: make_ndx_mpi -f conf.gro (select only BF4);"
    WRITE(*,*) "         g_traj_mpi -jump -oxt trajout.xtc; ./N-fluc.x -i trajout.xtc"
    STOP
    ELSE IF(NUM_OF_ARGS > 0 ) Then
    !loop across options
      DO I_ARG=1,NUM_OF_ARGS
        CALL get_command_argument(I_ARG,FILE_NAME)
        FILE_NAME = adjustl(FILE_NAME)
        SELECT CASE(FILE_NAME)
          CASE("-i")
            LookForInput=.TRUE.
            CYCLE
          CASE("-o1")
            LookForOutput1=.TRUE.
            CYCLE
          CASE("-o2")
            LookForOutput2=.TRUE.
            CYCLE
          CASE("-mol")
            LookForMOL=.TRUE.
            CYCLE
          CASE("-dz")
            LookForDZ=.TRUE.
            CYCLE
          CASE("-block")
            LookForBLOCK=.TRUE.
            CYCLE
          CASE default
            IF(LookForInput) THEN
              inquire(file=TRIM(FILE_NAME),exist=FileExist)
              IF(.NOT.FileExist) THEN
                WRITE(*,*) 'Input File ',FILE_NAME,' is not found'
                STOP
              EndIF
              INPUT_NAME = FILE_NAME
              LookForInput = .FALSE.
              CYCLE
              ELSE IF(LookForOutput1) THEN
                inquire(file=TRIM(FILE_NAME),exist=FileExist)
                IF(FileExist) THEN
                  WRITE(*,*) 'Output File1 ',FILE_NAME,' already exists'
                EndIF
                OUTPUT_NAME1 = FILE_NAME
                LookForOutput1 = .FALSE.
                CYCLE
              ELSE IF(LookForOutput2) THEN
                inquire(file=TRIM(FILE_NAME),exist=FileExist)
                IF(FileExist) THEN
                  WRITE(*,*) 'Output File2 ',FILE_NAME,' already exists'
                EndIF
                OUTPUT_NAME2 = FILE_NAME
                LookForOutput2 = .FALSE.
                CYCLE
              ELSE IF(LookForMOL) THEN
                READ(FILE_NAME,*) MOL_NAME
                LookForMOL = .FALSE.
                CYCLE
              ELSE IF(LookForDZ) THEN
                READ(FILE_NAME,*) D_Z
                LookForDZ = .FALSE.
                CYCLE
              ELSE IF(LookForBLOCK) THEN
                READ(FILE_NAME,*) SIZE_BLOCK
                LookForBLOCK = .FALSE.
                CYCLE
              ELSE
                PRINT*, "Option ", FILE_NAME, "is unknown"
                STOP
            ENDIF
         END SELECT
      ENDDO
  ENDIF
  PRINT*, "Your command line: N-fluc.x -i ",TRIM(INPUT_NAME),&
   " -o1 ",TRIM(OUTPUT_NAME1)," -o2 ",TRIM(OUTPUT_NAME2),&
   " -mol ",TRIM(MOL_NAME)," -dz ",D_Z,"-block",SIZE_BLOCK

! Read Input file
  PRINT*, "========================"
  PRINT*, "Reading Input File, ",TRIM(INPUT_NAME)

! Initialize it with the name of xtc file you want to read in.
  call xtc%init(TRIM(INPUT_NAME))
  I = 0
  call xtc%read
  IF ( xtc%STAT == 0 ) THEN
    ! check precision of xtc file
    IF( xtc%prec /= 1000.00D0 ) THEN
      PRINT*, "xtc file precision is not 1000.0 (real)!, but ",xtc%prec
      PRINT*, "this program does only support real precision."
      STOP
    ENDIF
!   check x and y of unit cell is same! (I only consider semiisotropic NPZAT ensemble)
    IF(xtc%box(1,1) /= xtc%box(2,2)) THEN
      PRINT*, "pbc box might not for NPZAT. This program does not support this system."
      STOP
    ENDIF
    ! save original position
    N_ATOM = xtc%NATOMS
    ALLOCATE(AR_DATA(1:3,1:N_ATOM,0:MAX_DATA),stat=IERR)
    IF (IERR /= 0) STOP "DATA: problem in attempt to allocate memory"
    AR_DATA = 0.0D0
    AR_DATA(:,:,0) = xtc%pos(:,:)
    ALLOCATE(BOX(1:3,0:MAX_DATA),stat=IERR)
    DO J=1,3
      BOX(J,I) = xtc%box(J,J)
    EndDO
    PRINT*, "semiisotropic initial BOX size (",(BOX(J,I),J=1,3),")"
    ELSE
      PRINT*, TRIM(INPUT_NAME)," has no trajectory. Check again!"
      STOP
  ENDIF
! save coordinates
  DO
    IF(MOD(I+1,500) == 0) PRINT*,"Reading ",I+1,"th time"
    ! read data
    call xtc%read
    !checking
    IF ( xtc%STAT /= 0 ) THEN
      PRINT*, "Reading COMPLETE"
      EXIT
    ENDIF
    IF((xtc%box(1,1) /= BOX(1,0)) .or. (xtc%box(2,2) /= BOX(2,0))) THEN
      PRINT*, "pbc box might be changed. this program does only support NPZAT system."
      STOP
    ENDIF
    I = I + 1
    ! check array memory
    IF(I == MAX_DATA) THEN
      PRINT*, "Memory increases from size of 3 x ",N_ATOM," x ",SIZE(AR_DATA)/(3*N_ATOM)
      ALLOCATE(TEMP(1:3,1:N_ATOM,0:MAX_DATA+ICR),stat=IERR)
      ALLOCATE(TEMP_BOX(1:3,0:MAX_DATA+ICR),stat=IERR)
      IF (IERR /= 0) STOP "reallocate_temp: problem in attempt to allocate memory"
      TEMP = 0.0D0
      TEMP_BOX = 0.0D0
      DO J=1,3
        DO K=1,N_ATOM
          DO L=0,MAX_DATA
            TEMP(J,K,L) = AR_DATA(J,K,L)
          ENDDO
        ENDDO
      ENDDO
      DO J=1,3
        DO L=0,MAX_DATA
          TEMP_BOX(J,L) = BOX(J,L)
        ENDDO
      ENDDO
      DEALLOCATE(AR_DATA)
      DEALLOCATE(BOX)
      ALLOCATE(AR_DATA(1:3,1:N_ATOM,0:MAX_DATA+ICR),stat=IERR)
      ALLOCATE(BOX(1:3,0:MAX_DATA+ICR),stat=IERR)
      IF (IERR /= 0) STOP "reallocate_array: problem in attempt to allocate memory"
      DO J=1,3
        DO K=1,N_ATOM
          DO L=0,MAX_DATA+ICR
            AR_DATA(J,K,L) = TEMP(J,K,L)
          ENDDO
        ENDDO
      ENDDO
      DO J=1,3
        DO L=0,MAX_DATA+ICR
          BOX(J,L) = TEMP_BOX(J,L)
        ENDDO
      ENDDO
      DEALLOCATE(TEMP)
      DEALLOCATE(TEMP_BOX)
      MAX_DATA = MAX_DATA+ICR
    ENDIF
    AR_DATA(:,:,I) = xtc%pos(:,:)
    DO J=1,3
      BOX(J,I) = xtc%box(J,J)
    EndDO
  ENDDO
  N_DATA = I
  PRINT*, "Total frame is ",N_DATA

! Close the file
  call xtc % close

! Just an example to show what was read in
! write(*,'(a,f12.6,a,i0)') " Time (ps): ", xtc % time, "  Step: ", xtc % STEP
! write(*,'(a,f12.6,a,i0)') " Precision: ", xtc % prec, "  No. Atoms: ", xtc % NATOMS
! write(*,'(3f9.3)') xtc % pos
! This is the same order as found in the GRO format fyi
! write(*,'(11f9.5)') xtc % box(1,1), xtc % box(2,2), xtc % box(3,3), &
!                     xtc % box(1,2), xtc % box(1,3), &
!                     xtc % box(2,1), xtc % box(2,3), &
!                     xtc % box(3,1), xtc % box(3,2)

! Calculate COM of the molecule
  PRINT*, "========================"
  PRINT*, "Calculating COM of molecules"

  SELECT CASE(MOL_NAME)
    CASE("bf4") !the molecules are BF4
      N_INTRA = 5
      IF(MOD(N_ATOM,N_INTRA) /= 0) THEN
        PRINT*, "Molecules are not ",TRIM(MOL_NAME)
        STOP
      ENDIF
      ALLOCATE(MASS(1:N_INTRA+1),stat=IERR)
      MASS=0.0D0
      ! mass information of BF4
      MASS(1)=10.811 ! mass of B
      DO I=2,N_INTRA
        MASS(I)=18.9984
      ENDDO
      ! sum of mass of a BF4 molecule
      DO I=1,N_INTRA
        MASS(N_INTRA+1)=MASS(N_INTRA+1)+MASS(I)
      ENDDO
    CASE("pf6") !the molecules are PF6
      N_INTRA = 7
      IF(MOD(N_ATOM,N_INTRA) /= 0) THEN
        PRINT*, "Molecules are not ",TRIM(MOL_NAME)
      ENDIF
      ALLOCATE(MASS(1:N_INTRA+1),stat=IERR)
      MASS=0.0D0
      ! mass information of BF4
      MASS(1)=30.97376 ! mass of P
      DO I=2,N_INTRA
        MASS(I)=18.9984
      ENDDO
      ! sum of mass of a BF4 molecule
      DO I=1,N_INTRA
        MASS(N_INTRA+1)=MASS(N_INTRA+1)+MASS(I)
      ENDDO
    CASE default
      PRINT*, "Cannot match the name of molecules you set!"
      STOP
  END SELECT

! COM
  ALLOCATE(COM(1:3,1:N_ATOM/N_INTRA,0:N_DATA),stat=IERR)
  IF (IERR /= 0) STOP "COM: problem in attempt to allocate memory"
  COM = 0.0D0
  DO I=0,N_DATA
    DO J=1,N_ATOM
      ! find centeral atom
      IF(MOD(J-1,N_INTRA) == 0) THEN
        AR_DATA(:,J,I)=AR_DATA(:,J,I)-BOX(:,I)*FLOOR(AR_DATA(:,J,I)/BOX(:,I))
        ! translate long distance atom
        ELSE
          AR_DATA(:,J,I)=AR_DATA(:,J,I)-BOX(:,I)*NINT((AR_DATA(:,J,I)-AR_DATA(:,J-MOD(J-1,N_INTRA),I))/BOX(:,I))
      ENDIF
      I_COM = (J-1)/N_INTRA+1
      COM(:,I_COM,I)=COM(:,I_COM,I)+AR_DATA(:,J,I)*MASS(MOD(J-1,N_INTRA)+1)
    EndDO
  EndDO
  COM=COM/MASS(N_INTRA+1)
  PRINT*, "calculating COMPLETE"

! Slicing slabs
  PRINT*, "========================"
  PRINT*, "Slicing the simulation box and counting molecules"

  OPEN(UNIT=11,FILE=TRIM(OUTPUT_NAME1),STATUS='unknown')
! # of Slices from initial box_z length
  IF(D_Z <= 0.1D0) THEN
    PRINT*, "thickness of slices, ",D_Z," is too small, < 1 A"
  ENDIF
  N_SLICE = CEILING(BOX(3,0)/D_Z)
  PRINT*, N_SLICE,"# of slices"
  ALLOCATE(NPAR(1:N_SLICE,0:N_DATA),stat=IERR)
  IF (IERR /= 0) STOP "NPAR: problem in attempt to allocate memory"
  NPAR = 0
  ALLOCATE(T_SLICE(0:N_DATA),stat=IERR)
  IF (IERR /= 0) STOP "T_SLICE: problem in attempt to allocate memory"
  T_SLICE = 0
  DO I=0,N_DATA
    T_SLICE(I) = BOX(3,I)/REAL(N_SLICE)
    DO J=1,N_ATOM/N_INTRA
      COM(3,J,I)=COM(3,J,I)-BOX(3,I)*FLOOR(COM(3,J,I)/BOX(3,I))
      INDEX = CEILING(COM(3,J,I)/T_SLICE(I))
!      PRINT*, INDEX, COM(3,J,I)/T_SLICE, CEILING(COM(3,J,I)/T_SLICE), COM(3,J,I)/BOX(3,I),FLOOR(COM(3,J,I)/BOX(3,I))
      IF(INDEX <= 0) THEN
        PRINT*, COM(:,J,I),I,J
        STOP "index of slabs is less than 0 or zero"
      ENDIF
      IF(INDEX > N_SLICE) STOP "index of slabs is greater than MAX"
      NPAR(INDEX,I) = NPAR(INDEX,I) + 1.0D0
!      PRINT*, INDEX,I,NPAR(INDEX,I),INT(NPAR(INDEX,I))
    ENDDO
    DO K=1,N_SLICE
      WRITE(11,'(I8,1X,I3,1X,I5)') I,K,INT(NPAR(K,I))
    ENDDO
    WRITE(11,*)
    WRITe(11,*)
  ENDDO
  CLOSE(11)
  PRINT*, "Writing ",TRIM(OUTPUT_NAME1)," COMPLETE"
  PRINT*, "print NPAR: Time, i-th slab, N_particles in the i-th slab"

  !Block average (time average)
  PRINT*, "========================"
  PRINT*, "Block averaging"
  DO K=1,N_SLICE
    DO I=MOD(N_DATA+1,SIZE_BLOCK),N_DATA-SIZE_BLOCK+1,SIZE_BLOCK ! neglect frist steps in each block
      DO J=I+SIZE_BLOCK-1,I+1,-1 ! sum from last number to first number within size of block
        NPAR(K,J-1) = NPAR(K,J-1)+NPAR(K,J)
      ENDDO
      NPAR(K,I)=NPAR(K,I)/DBLE(SIZE_BLOCK)
!      PRINT*,K,I,NPAR(K,I)
    ENDDO
  ENDDO

  !Calculating avg and std
  PRINT*, "========================"
  PRINT*, "Calculating average and std"
  ALLOCATE(NPAR_AVG(1:N_SLICE),stat=IERR)
  ALLOCATE(NPAR_STD2(1:N_SLICE),stat=IERR)
  NPAR_AVG = 0.0D0
  NPAR_STD2 = 0.0D0
  DO K=1,N_SLICE
    DO I=MOD(N_DATA+1,SIZE_BLOCK),N_DATA-SIZE_BLOCK+1,SIZE_BLOCK
!      PRINT*, "I",I
      NPAR_AVG(K) = NPAR_AVG(K)+NPAR(K,I)
      NPAR_STD2(K) = NPAR_STD2(K)+NPAR(K,I)**2
    ENDDO
    NPAR_AVG(K) = NPAR_AVG(K)/DBLE((N_DATA+1)/SIZE_BLOCK)
    NPAR_STD2(K) = NPAR_STD2(K)/DBLE((N_DATA+1)/SIZE_BLOCK)-NPAR_AVG(K)**2
  ENDDO

  OPEN(UNIT=12,FILE=TRIM(OUTPUT_NAME2),STATUS='unknown')
  DO K=1,N_SLICE
    WRITE(12,'(I5,3(1X,F10.5))') K, NPAR_AVG(K), NPAR_STD2(K), NPAR_STD2(K)/(NPAR_AVG(K)**2)
  ENDDO
  PRINT*, "Writing ",TRIM(OUTPUT_NAME2)," COMPLETE"
  PRINT*, "print NFLUC: i-th slab, avg, std**2, std**2/avg**2"
  CLOSE(12)

  DO I=0,N_DATA
    T_SLICE_AVG = T_SLICE_AVG+T_SLICE(I)
    T_SLICE_STD = T_SLICE_STD+T_SLICE(I)**2
  EndDO
  T_SLICE_AVG = T_SLICE_AVG/(N_DATA+1)
  T_SLICE_STD = SQRT(T_SLICE_STD/(N_DATA+1)-T_SLICE_AVG**2)

  PRINT*, "The following information is not considering block average"
  PRINT*, "avg thickness of slice (nm) =",SUM(T_SLICE)/(N_DATA+1)
  PRINT*, "STD thickness of slice (nm) =",T_SLICE_STD,"(",T_SLICE_STD*100/(SUM(T_SLICE)/(N_DATA+1)),"%)"
  PRINT*, "EXIT N-FLUC PROGRAM"

END PROGRAM NFLUC


text:
      subroutine ldmas(r,u,ui,udt,udtt,id,inod,xmass,yz,fval,
     1        jnod,xdamp,iter,roller,lm,sf,iprec)
c     implicit double precision (a-h,o-z)                                    dp
c
c     concentrated nodal dampers, masses, and roller boundary conditions
c
      common/bk02/ioofc,iphase,imass,lpar(9)
      common/bk03/numdc,imassn,idampn,irller,penstf
      common/bk08/kprint,nstep,ite,ilimit,newstf
      common/bk12/ntlen
      common/bk17/dn1,dn2,nwebuf,ntime,numnp,neq,ibar,mthsol
      common/bk26/dt,a0,a1,a2,a3,a4,a5,a6,a7,a8,a9,a10
      common/bk27/nlcur,nptst,nthpy,nthpz,nthps,xmy,xmz,xms,nload,nptm
      common/fissn1/melemt,nnns,ntpe2,n2g,llls
      common/main_block/ a(1)
c
      dimension r(*),u(*),ui(*),udt(*),udtt(*),id(2,*),inod(*),
     1        xmass(2,*),yz(2,*),fval(*),jnod(*),xdamp(2,*),roller(*),
     2        lm(5*iprec,*),sf(5,*)
c
      equivalence (nnns,nd),(llls,nms)
c
      nnns=2
      llls=5
c
      if (imassn.eq.0) go to 110
c
      lcount=0
      do 100 i=1,imassn
      node=inod(i)
      sf(3,i)=0.
      sf(4,i)=0.
      sf(5,i)=0.
      lm(1,i)=id(1,node)
      lm(2,i)=id(2,node)
      ly=lm(1,i)
      lz=lm(2,i)
      xmassy=xmass(1,i)
      xmassz=xmass(2,i)
c
      if (imass.ne.1) go to 40
c
      if (ly.eq.0) go to 30
      wv=-a2*udt(ly)-a3*udtt(ly)
      if (iter.eq.1) wv=wv+a0*ui(ly)
      r(ly)=r(ly)-wv*xmassy
      sf(3,i)=sf(3,i)+a0*xmassy
   30 if (lz.eq.0) go to 40
      wv=-a2*udt(lz)-a3*udtt(lz)
      if (iter.eq.1) wv=wv+a0*ui(lz)
      r(lz)=r(lz)-wv*xmassz
      sf(5,i)=sf(5,i)+a0*xmassz
c
c     body force loads
c
   40 if (nthpy+nthpz+abs(nthps).eq.0) go to 90
c
      if (nthpy.eq.0.or.ly.eq.0) go to 50
      r(ly)=r(ly)-xmassy*fval(nthpy)*xmy
   50 if (nthpz.eq.0.or.lz.eq.0) go to 60
      r(lz)=r(lz)-xmassz*fval(nthpz)*xmz
   60 if (nthps.eq.0) go to 90
      ispny=0
      if (nthps.lt.0) ispny=1
      nthps=abs(nthps)
      facs=(fval(nthps)*xms)**2
      if (ly.eq.0) go to 70
      radius=yz(1,i)
      radius=radius+u(ly)
      r(ly)=r(ly)-xmassy*radius*facs
   70 if (ispny.eq.0.or.lz.eq.0) go to 80
      radius=yz(2,i)
      radius=radius+u(lz)
      r(lz)=r(lz)-xmassz*radius*facs
   80 if (ispny.eq.1) nthps=-nthps
   90 lcount=lcount+1
      if (lcount.lt.64) go to 100
      if (iphase.eq.3) go to 100
      if (newstf.ne.0) go to 100
      melemt=64
      lcount=0
      if (ioofc.eq.1) call fissln (-5,a(ntlen),a(ntlen),lm)
      if (ioofc.eq.0) call wdiskf (lm,melemt)
  100 continue
      if (lcount.eq.0) go to 110
      if (iphase.eq.3) go to 110
      if (newstf.ne.0) go to 110
      melemt=lcount
      if (ioofc.eq.1) call fissln (-5,a(ntlen),a(ntlen),lm)
      if (ioofc.eq.0) call wdiskf (lm,melemt)
c
  110 if (idampn.eq.0) go to 150
c
      lcount=0
      do 140 i=1,idampn
      node=jnod(i)
      sf(3,i)=0.
      sf(4,i)=0.
      sf(5,i)=0.
      lm(1,i)=id(1,node)
      lm(2,i)=id(2,node)
      ly=lm(1,i)
      lz=lm(2,i)
      xdampy=xdamp(1,i)
      xdampz=xdamp(2,i)
      if (ly.eq.0) go to 120
      wv=-a4*udt(ly)-a5*udtt(ly)
      if (iter.eq.1) wv=wv+a1*ui(ly)
      r(ly)=r(ly)-wv*xdampy
      sf(3,i)=sf(3,i)+a1*xdampy
  120 if (lz.eq.0) go to 130
      wv=-a4*udt(lz)-a5*udtt(lz)
      if (iter.eq.1) wv=wv+a1*ui(lz)
      r(lz)=r(lz)-wv*xdampz
      sf(5,i)=sf(5,i)+a1*xdampz
  130 lcount=lcount+1
      if (lcount.lt.64) go to 140
      if (iphase.eq.3) go to 140
      if (newstf.ne.0) go to 140
      melemt=64
      lcount=0
      if (ioofc.eq.1) call fissln (-5,a(ntlen),a(ntlen),lm)
      if (ioofc.eq.0) call wdiskf (lm,melemt)
  140 continue
      if (lcount.eq.0) go to 150
      if (iphase.eq.3) go to 150
      if (newstf.ne.0) go to 150
      melemt=lcount
      if (ioofc.eq.1) call fissln (-5,a(ntlen),a(ntlen),lm)
      if (ioofc.eq.0) call wdiskf (lm,melemt)
c
  150 if(irller.eq.0) return
c
      lcount=0
      do 160 i=1,numnp
      if(roller(i).eq.0.) go to 160
      angle=roller(i)/57.29577951
      sinthe=sin(angle)
      costhe=cos(angle)
      lcount=lcount+1
      lm(1,lcount)=id(1,i)
      lm(2,lcount)=id(2,i)
      sf(3,lcount)= penstf*sinthe*sinthe
      sf(4,lcount)=-penstf*sinthe*costhe
      sf(5,lcount)= penstf*costhe*costhe
      rtan=r(id(1,i))*costhe+r(id(2,i))*sinthe
      r(id(1,i))=rtan*costhe
      r(id(2,i))=rtan*sinthe
      if (lcount.lt.64) go to 160
      if (iphase.eq.3) go to 160
      if (newstf.ne.0) go to 160
      melemt=64
      lcount=0
      if (ioofc.eq.1) call fissln (-5,a(ntlen),a(ntlen),lm)
      if (ioofc.eq.0) call wdiskf (lm,melemt)
  160 continue
      if (lcount.eq.0) return
      if (iphase.eq.3) return
      if (newstf.ne.0) return
      melemt=lcount
      if (ioofc.eq.1) call fissln (-5,a(ntlen),a(ntlen),lm)
      if (ioofc.eq.0) call wdiskf (lm,melemt)
c
      return
c
      end

max_stars_repo_name: imohame/LabCode

max_issues_repo_name: imohame/LabCode

max_forks_repo_name: imohame/LabCode


text:
C
C	$Id$
C
      SUBROUTINE CSTRFIND (NST,P,N,X,Y,Z,LIST,LPTR,LEND, B1,
     .                   B2,B3,I1,I2,I3)
      INTEGER NST, N, LIST(*), LPTR(*), LEND(N), I1, I2, I3
      DOUBLE PRECISION P(3), X(N), Y(N), Z(N), B1, B2, B3
C
C***********************************************************
C
C                                              From STRIPACK
C                                            Robert J. Renka
C                                  Dept. of Computer Science
C                                       Univ. of North Texas
C                                           renka@cs.unt.edu
C                                                   11/30/99
C
C   This subroutine locates a point P relative to a triangu-
C lation created by Subroutine CSTRMESH.  If P is contained in
C a triangle, the three vertex indexes and barycentric coor-
C dinates are returned.  Otherwise, the indexes of the
C visible boundary nodes are returned.
C
C
C On input:
C
C       NST = Index of a node at which CSTRFIND begins its
C             search.  Search time depends on the proximity
C             of this node to P.
C
C       P = Array of length 3 containing the x, y, and z
C           coordinates (in that order) of the point P to be
C           located.
C
C       N = Number of nodes in the triangulation.  N .GE. 3.
C
C       X,Y,Z = Arrays of length N containing the Cartesian
C               coordinates of the triangulation nodes (unit
C               vectors).  (X(I),Y(I),Z(I)) defines node I
C               for I = 1 to N.
C
C       LIST,LPTR,LEND = Data structure defining the trian-
C                        gulation.  Refer to Subroutine
C                        CSTRMESH.
C
C Input parameters are not altered by this routine.
C
C On output:
C
C       B1,B2,B3 = Unnormalized barycentric coordinates of
C                  the central projection of P onto the un-
C                  derlying planar triangle if P is in the
C                  convex hull of the nodes.  These parame-
C                  ters are not altered if I1 = 0.
C
C       I1,I2,I3 = Counterclockwise-ordered vertex indexes
C                  of a triangle containing P if P is con-
C                  tained in a triangle.  If P is not in the
C                  convex hull of the nodes, I1 and I2 are
C                  the rightmost and leftmost (boundary)
C                  nodes that are visible from P, and
C                  I3 = 0.  (If all boundary nodes are vis-
C                  ible from P, then I1 and I2 coincide.)
C                  I1 = I2 = I3 = 0 if P and all of the
C                  nodes are coplanar (lie on a common great
C                  circle.
C
C Modules required by CSTRFIND:  CSJRAND, CSLSTPTR, CSSTORE
C
C Intrinsic function called by CSTRFIND:  ABS
C
C***********************************************************
C
      INTEGER CSJRAND, CSLSTPTR
      INTEGER IX, IY, IZ, LP, N0, N1, N1S, N2, N2S, N3, N4,
     .        NEXT, NF, NL
      DOUBLE PRECISION CSSTORE
      DOUBLE PRECISION DET, EPS, PTN1, PTN2, Q(3), S12, TOL,
     .                 XP, YP, ZP
      DOUBLE PRECISION X0, X1, X2, Y0, Y1, Y2, Z0, Z1, Z2
C
      SAVE    IX, IY, IZ
      DATA    IX/1/, IY/2/, IZ/3/
C
C Local parameters:
C
C EPS =      Machine precision
C IX,IY,IZ = Integer seeds for CSJRAND
C LP =       LIST pointer
C N0,N1,N2 = Nodes in counterclockwise order defining a
C              cone (with vertex N0) containing P, or end-
C              points of a boundary edge such that P Right
C              N1->N2
C N1S,N2S =  Initially-determined values of N1 and N2
C N3,N4 =    Nodes opposite N1->N2 and N2->N1, respectively
C NEXT =     Candidate for I1 or I2 when P is exterior
C NF,NL =    First and last neighbors of N0, or first
C              (rightmost) and last (leftmost) nodes
C              visible from P when P is exterior to the
C              triangulation
C PTN1 =     Scalar product <P,N1>
C PTN2 =     Scalar product <P,N2>
C Q =        (N2 X N1) X N2  or  N1 X (N2 X N1) -- used in
C              the boundary traversal when P is exterior
C S12 =      Scalar product <N1,N2>
C TOL =      Tolerance (multiple of EPS) defining an upper
C              bound on the magnitude of a negative bary-
C              centric coordinate (B1 or B2) for P in a
C              triangle -- used to avoid an infinite number
C              of restarts with 0 <= B3 < EPS and B1 < 0 or
C              B2 < 0 but small in magnitude
C XP,YP,ZP = Local variables containing P(1), P(2), and P(3)
C X0,Y0,Z0 = Dummy arguments for DET
C X1,Y1,Z1 = Dummy arguments for DET
C X2,Y2,Z2 = Dummy arguments for DET
C
C Statement function:
C
C DET(X1,...,Z0) .GE. 0 if and only if (X0,Y0,Z0) is in the
C                       (closed) left hemisphere defined by
C                       the plane containing (0,0,0),
C                       (X1,Y1,Z1), and (X2,Y2,Z2), where
C                       left is defined relative to an ob-
C                       server at (X1,Y1,Z1) facing
C                       (X2,Y2,Z2).
C
      DET (X1,Y1,Z1,X2,Y2,Z2,X0,Y0,Z0) = X0*(Y1*Z2-Y2*Z1)
     .     - Y0*(X1*Z2-X2*Z1) + Z0*(X1*Y2-X2*Y1)
C
C Initialize variables.
C
      XP = P(1)
      YP = P(2)
      ZP = P(3)
      N0 = NST
      IF (N0 .LT. 1  .OR.  N0 .GT. N)
     .  N0 = CSJRAND(N, IX,IY,IZ )
C
C Compute the relative machine precision EPS and TOL.
C
      EPS = 1.D0
    1 EPS = EPS/2.D0
        IF (CSSTORE(EPS+1.D0) .GT. 1.D0) GO TO 1
      EPS = 2.D0*EPS
      TOL = 4.D0*EPS
C
C Set NF and NL to the first and last neighbors of N0, and
C   initialize N1 = NF.
C
    2 LP = LEND(N0)
      NL = LIST(LP)
      LP = LPTR(LP)
      NF = LIST(LP)
      N1 = NF
C
C Find a pair of adjacent neighbors N1,N2 of N0 that define
C   a wedge containing P:  P CSLEFT N0->N1 and P RIGHT N0->N2.
C
      IF (NL .GT. 0) THEN
C
C   N0 is an interior node.  Find N1.
C
    3   IF ( DET(X(N0),Y(N0),Z(N0),X(N1),Y(N1),Z(N1),
     .           XP,YP,ZP) .LT. 0.D0 ) THEN
          LP = LPTR(LP)
          N1 = LIST(LP)
          IF (N1 .EQ. NL) GO TO 6
          GO TO 3
        ENDIF
      ELSE
C
C   N0 is a boundary node.  Test for P exterior.
C
        NL = -NL
        IF ( DET(X(N0),Y(N0),Z(N0),X(NF),Y(NF),Z(NF),
     .           XP,YP,ZP) .LT. 0.D0 ) THEN
C
C   P is to the right of the boundary edge N0->NF.
C
          N1 = N0
          N2 = NF
          GO TO 9
        ENDIF
        IF ( DET(X(NL),Y(NL),Z(NL),X(N0),Y(N0),Z(N0),
     .           XP,YP,ZP) .LT. 0.D0 ) THEN
C
C   P is to the right of the boundary edge NL->N0.
C
          N1 = NL
          N2 = N0
          GO TO 9
        ENDIF
      ENDIF
C
C P is to the left of arcs N0->N1 and NL->N0.  Set N2 to the
C   next neighbor of N0 (following N1).
C
    4 LP = LPTR(LP)
        N2 = ABS(LIST(LP))
        IF ( DET(X(N0),Y(N0),Z(N0),X(N2),Y(N2),Z(N2),
     .           XP,YP,ZP) .LT. 0.D0 ) GO TO 7
        N1 = N2
        IF (N1 .NE. NL) GO TO 4
      IF ( DET(X(N0),Y(N0),Z(N0),X(NF),Y(NF),Z(NF),
     .         XP,YP,ZP) .LT. 0.D0 ) GO TO 6
C
C P is left of or on arcs N0->NB for all neighbors NB
C   of N0.  Test for P = +/-N0.
C
      IF (CSSTORE(ABS(X(N0)*XP + Y(N0)*YP + Z(N0)*ZP))
     .   .LT. 1.D0-4.D0*EPS) THEN
C
C   All points are collinear iff P Left NB->N0 for all
C     neighbors NB of N0.  Search the neighbors of N0.
C     Note:  N1 = NL and LP points to NL.
C
    5   IF ( DET(X(N1),Y(N1),Z(N1),X(N0),Y(N0),Z(N0),
     .           XP,YP,ZP) .GE. 0.D0 ) THEN
          LP = LPTR(LP)
          N1 = ABS(LIST(LP))
          IF (N1 .EQ. NL) GO TO 14
          GO TO 5
        ENDIF
      ENDIF
C
C P is to the right of N1->N0, or P = +/-N0.  Set N0 to N1
C   and start over.
C
      N0 = N1
      GO TO 2
C
C P is between arcs N0->N1 and N0->NF.
C
    6 N2 = NF
C
C P is contained in a wedge defined by geodesics N0-N1 and
C   N0-N2, where N1 is adjacent to N2.  Save N1 and N2 to
C   test for cycling.
C
    7 N3 = N0
      N1S = N1
      N2S = N2
C
C Top of edge-hopping loop:
C
    8 B3 = DET(X(N1),Y(N1),Z(N1),X(N2),Y(N2),Z(N2),XP,YP,ZP)
      IF (B3 .LT. 0.D0) THEN
C
C   Set N4 to the first neighbor of N2 following N1 (the
C     node opposite N2->N1) unless N1->N2 is a boundary arc.
C
        LP = CSLSTPTR(LEND(N2),N1,LIST,LPTR)
        IF (LIST(LP) .LT. 0) GO TO 9
        LP = LPTR(LP)
        N4 = ABS(LIST(LP))
C
C   Define a new arc N1->N2 which intersects the geodesic
C     N0-P.
C
        IF ( DET(X(N0),Y(N0),Z(N0),X(N4),Y(N4),Z(N4),
     .           XP,YP,ZP) .LT. 0.D0 ) THEN
          N3 = N2
          N2 = N4
          N1S = N1
          IF (N2 .NE. N2S  .AND.  N2 .NE. N0) GO TO 8
        ELSE
          N3 = N1
          N1 = N4
          N2S = N2
          IF (N1 .NE. N1S  .AND.  N1 .NE. N0) GO TO 8
        ENDIF
C
C   The starting node N0 or edge N1-N2 was encountered
C     again, implying a cycle (infinite loop).  Restart
C     with N0 randomly selected.
C
        N0 = CSJRAND(N, IX,IY,IZ )
        GO TO 2
      ENDIF
C
C P is in (N1,N2,N3) unless N0, N1, N2, and P are collinear
C   or P is close to -N0.
C
      IF (B3 .GE. EPS) THEN
C
C   B3 .NE. 0.
C
        B1 = DET(X(N2),Y(N2),Z(N2),X(N3),Y(N3),Z(N3),
     .           XP,YP,ZP)
        B2 = DET(X(N3),Y(N3),Z(N3),X(N1),Y(N1),Z(N1),
     .           XP,YP,ZP)
        IF (B1 .LT. -TOL  .OR.  B2 .LT. -TOL) THEN
C
C   Restart with N0 randomly selected.
C
          N0 = CSJRAND(N, IX,IY,IZ )
          GO TO 2
        ENDIF
      ELSE
C
C   B3 = 0 and thus P lies on N1->N2. Compute
C     B1 = Det(P,N2 X N1,N2) and B2 = Det(P,N1,N2 X N1).
C
        B3 = 0.D0
        S12 = X(N1)*X(N2) + Y(N1)*Y(N2) + Z(N1)*Z(N2)
        PTN1 = XP*X(N1) + YP*Y(N1) + ZP*Z(N1)
        PTN2 = XP*X(N2) + YP*Y(N2) + ZP*Z(N2)
        B1 = PTN1 - S12*PTN2
        B2 = PTN2 - S12*PTN1
        IF (B1 .LT. -TOL  .OR.  B2 .LT. -TOL) THEN
C
C   Restart with N0 randomly selected.
C
          N0 = CSJRAND(N, IX,IY,IZ )
          GO TO 2
        ENDIF
      ENDIF
C
C P is in (N1,N2,N3).
C
      I1 = N1
      I2 = N2
      I3 = N3
      IF (B1 .LT. 0.0) B1 = 0.0
      IF (B2 .LT. 0.0) B2 = 0.0
      RETURN
C
C P Right N1->N2, where N1->N2 is a boundary edge.
C   Save N1 and N2, and set NL = 0 to indicate that
C   NL has not yet been found.
C
    9 N1S = N1
      N2S = N2
      NL = 0
C
C           Counterclockwise Boundary Traversal:
C
   10 LP = LEND(N2)
      LP = LPTR(LP)
      NEXT = LIST(LP)
      IF ( DET(X(N2),Y(N2),Z(N2),X(NEXT),Y(NEXT),Z(NEXT),
     .         XP,YP,ZP) .GE. 0.D0 ) THEN
C
C   N2 is the rightmost visible node if P Forward N2->N1
C     or NEXT Forward N2->N1.  Set Q to (N2 X N1) X N2.
C
        S12 = X(N1)*X(N2) + Y(N1)*Y(N2) + Z(N1)*Z(N2)
        Q(1) = X(N1) - S12*X(N2)
        Q(2) = Y(N1) - S12*Y(N2)
        Q(3) = Z(N1) - S12*Z(N2)
        IF (XP*Q(1) + YP*Q(2) + ZP*Q(3) .GE. 0.D0) GO TO 11
        IF (X(NEXT)*Q(1) + Y(NEXT)*Q(2) + Z(NEXT)*Q(3)
     .      .GE. 0.D0) GO TO 11
C
C   N1, N2, NEXT, and P are nearly collinear, and N2 is
C     the leftmost visible node.
C
        NL = N2
      ENDIF
C
C Bottom of counterclockwise loop:
C
      N1 = N2
      N2 = NEXT
      IF (N2 .NE. N1S) GO TO 10
C
C All boundary nodes are visible from P.
C
      I1 = N1S
      I2 = N1S
      I3 = 0
      RETURN
C
C N2 is the rightmost visible node.
C
   11 NF = N2
      IF (NL .EQ. 0) THEN
C
C Restore initial values of N1 and N2, and begin the search
C   for the leftmost visible node.
C
        N2 = N2S
        N1 = N1S
C
C           Clockwise Boundary Traversal:
C
   12   LP = LEND(N1)
        NEXT = -LIST(LP)
        IF ( DET(X(NEXT),Y(NEXT),Z(NEXT),X(N1),Y(N1),Z(N1),
     .           XP,YP,ZP) .GE. 0.D0 ) THEN
C
C   N1 is the leftmost visible node if P or NEXT is
C     forward of N1->N2.  Compute Q = N1 X (N2 X N1).
C
          S12 = X(N1)*X(N2) + Y(N1)*Y(N2) + Z(N1)*Z(N2)
          Q(1) = X(N2) - S12*X(N1)
          Q(2) = Y(N2) - S12*Y(N1)
          Q(3) = Z(N2) - S12*Z(N1)
          IF (XP*Q(1) + YP*Q(2) + ZP*Q(3) .GE. 0.D0)
     .      GO TO 13
          IF (X(NEXT)*Q(1) + Y(NEXT)*Q(2) + Z(NEXT)*Q(3)
     .        .GE. 0.D0) GO TO 13
C
C   P, NEXT, N1, and N2 are nearly collinear and N1 is the
C     rightmost visible node.
C
          NF = N1
        ENDIF
C
C Bottom of clockwise loop:
C
        N2 = N1
        N1 = NEXT
        IF (N1 .NE. N1S) GO TO 12
C
C All boundary nodes are visible from P.
C
        I1 = N1
        I2 = N1
        I3 = 0
        RETURN
C
C N1 is the leftmost visible node.
C
   13   NL = N1
      ENDIF
C
C NF and NL have been found.
C
      I1 = NF
      I2 = NL
      I3 = 0
      RETURN
C
C All points are collinear (coplanar).
C
   14 I1 = 0
      I2 = 0
      I3 = 0
      RETURN
      END

max_stars_repo_name: xylar/cdat

max_issues_repo_name: xylar/cdat

max_forks_repo_name: CDAT/uvcdat


text:
	SUBROUTINE DMESG  ( messag, iret )
C************************************************************************
C* DMESG								*
C* 									*
C* DMESG  ( MESSAG, IRET )						*
C*									*
C* Input parameters:							*
C* 	MESSAG		CHAR*		Message
C* 									*
C* Output parameters:							*
C* 	IRET		INTEGER		Return code			*
C**									*
C* Log:									*
C* S. Jacobs/EAI	11/92						*
C************************************************************************
	INCLUDE		'ERROR.PRM'
C*
	CHARACTER*(*)	messag
C------------------------------------------------------------------------
	iret = NORMAL
C*
	CALL HMESG ( messag, ier )
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
{ A quine. Prints its own source code.
  Written by Mathew Hendry <m.hendry@dial.pipex.com>. Taken from <http://rosettacode.org/wiki/Quine#FALSE>. }

["'[,34,$!34,'],!"]'[,34,$!34,'],!

max_stars_repo_name: tadeassoucek/FalseCS

max_issues_repo_name: tadeassoucek/FalseCS

max_forks_repo_name: tadeassoucek/FalseCS