text:
module intrinsic_mm_driver_module
    use, intrinsic :: iso_fortran_env, only : &
            real32, &
            real64, &
            int64

    use :: data_api, only : stream

    use :: tensor_api, only : tensor_fortran_converter

    use :: fortran_data_mm_driver_module, only : fortran_data_mm_driver

    implicit none
    private

    public :: intrinsic_mm_driver

    type, extends(fortran_data_mm_driver) :: intrinsic_mm_driver
    contains
        procedure :: perform_mm_r32 => perform_mm_r32
        procedure :: perform_mm_r64 => perform_mm_r64
        procedure :: perform_mm_c64 => perform_mm_c64
        procedure :: perform_mm_c128 => perform_mm_c128
        procedure :: synchronize => synchronize
        procedure :: cleanup => cleanup
    end type intrinsic_mm_driver

    interface intrinsic_mm_driver
        module procedure constructor
    end interface intrinsic_mm_driver
contains
    function constructor() result(this)
        type(intrinsic_mm_driver) :: this

        this%converter = tensor_fortran_converter()
    end function constructor

    subroutine perform_mm_r32(this, dst, left, right, m, n, k, &
            alpha, beta, transposed_left, transposed_right, astream)
        class(intrinsic_mm_driver), intent(in) :: this
        real(real32), dimension(:,:), contiguous, intent(inout) :: dst
        real(real32), dimension(:,:), contiguous, intent(in) :: left, right
        integer(int64), intent(in) :: m, n, k
        real(real32), intent(in) :: alpha, beta
        logical, intent(in) :: transposed_left, transposed_right
        type(stream), intent(in), optional :: astream

        if ( all(shape(dst) /= [m, n]) ) &
                error stop "intrinsic_mm_driver::perform_mm_r32:Strided mm not implemented."

        if ( transposed_left .and. transposed_right ) then
            if ( all(shape(left) /= [k, m]) .or. all(shape(right) /= [n, k]) ) &
                    error stop "intrinsic_mm_driver::perform_mm_r32:Strided mm not implemented."

            if ( beta == 0.0 ) then
                dst = alpha*matmul(transpose(left), transpose(right))
            else
                dst = beta*dst + alpha*matmul(transpose(left), transpose(right))
            end if
        else if (transposed_left) then
            if ( all(shape(left) /= [k, m]) .or. all(shape(right) /= [k, n]) ) &
                    error stop "intrinsic_mm_driver::perform_mm_r32:Strided mm not implemented."

            if ( beta == 0.0 ) then
                dst = alpha*matmul(transpose(left), right)
            else
                dst = beta*dst + alpha*matmul(transpose(left), right)
            end if
        else if ( transposed_right ) then
            if ( all(shape(left) /= [m, k]) .or. all(shape(right) /= [n, k]) ) &
                    error stop "intrinsic_mm_driver::perform_mm_r32:Strided mm not implemented."

            if ( beta == 0.0 ) then
                dst = alpha*matmul(left, transpose(right))
            else
                dst = beta*dst + alpha*matmul(left, transpose(right))
            end if
        else
            if ( all(shape(left) /= [m, k]) .or. all(shape(right) /= [k, n]) ) &
                    error stop "intrinsic_mm_driver::perform_mm_r32:Strided mm not implemented."

            if ( beta == 0.0 ) then
                dst = alpha*matmul(left, right)
            else
                dst = beta*dst + alpha*matmul(left, right)
            end if
        end if
    end subroutine perform_mm_r32

    subroutine perform_mm_r64(this, dst, left, right, m, n, k, &
            alpha, beta, transposed_left, transposed_right, astream)
        class(intrinsic_mm_driver), intent(in) :: this
        real(real64), dimension(:,:), contiguous, intent(inout) :: dst
        real(real64), dimension(:,:), contiguous, intent(in) :: left, right
        integer(int64), intent(in) :: m, n, k
        real(real64), intent(in) :: alpha, beta
        logical, intent(in) :: transposed_left, transposed_right
        type(stream), intent(in), optional :: astream

        if ( all(shape(dst) /= [m, n]) ) &
                error stop "intrinsic_mm_driver::perform_mm_r64:Strided mm not implemented."

        if ( transposed_left .and. transposed_right ) then
            if ( all(shape(left) /= [k, m]) .or. all(shape(right) /= [n, k]) ) &
                    error stop "intrinsic_mm_driver::perform_mm_r64:Strided mm not implemented."

            if ( beta == 0.0d0 ) then
                dst = alpha*matmul(transpose(left), transpose(right))
            else
                dst = beta*dst + alpha*matmul(transpose(left), transpose(right))
            end if
        else if (transposed_left) then
            if ( all(shape(left) /= [k, m]) .or. all(shape(right) /= [k, n]) ) &
                    error stop "intrinsic_mm_driver::perform_mm_r64:Strided mm not implemented."

            if ( beta == 0.0d0 ) then
                dst = alpha*matmul(transpose(left), right)
            else
                dst = beta*dst + alpha*matmul(transpose(left), right)
            end if
        else if ( transposed_right ) then
            if ( all(shape(left) /= [m, k]) .or. all(shape(right) /= [n, k]) ) &
                    error stop "intrinsic_mm_driver::perform_mm_r64:Strided mm not implemented."

            if ( beta == 0.0d0 ) then
                dst = alpha*matmul(left, transpose(right))
            else
                dst = beta*dst + alpha*matmul(left, transpose(right))
            end if
        else
            if ( all(shape(left) /= [m, k]) .or. all(shape(right) /= [k, n]) ) &
                    error stop "intrinsic_mm_driver::perform_mm_r64:Strided mm not implemented."

            if ( beta == 0.0d0 ) then
                dst = alpha*matmul(left, right)
            else
                dst = beta*dst + alpha*matmul(left, right)
            end if
        end if
    end subroutine perform_mm_r64

    subroutine perform_mm_c64(this, dst, left, right, m, n, k, &
            alpha, beta, transposed_left, transposed_right, astream)
        class(intrinsic_mm_driver), intent(in) :: this
        complex(real32), dimension(:,:), contiguous, intent(inout) :: dst
        complex(real32), dimension(:,:), contiguous, intent(in) :: left, right
        integer(int64), intent(in) :: m, n, k
        complex(real32), intent(in) :: alpha, beta
        logical, intent(in) :: transposed_left, transposed_right
        type(stream), intent(in), optional :: astream

        if ( all(shape(dst) /= [m, n]) ) &
                error stop "intrinsic_mm_driver::perform_mm_c64:Strided mm not implemented."

        if ( transposed_left .and. transposed_right ) then
            if ( all(shape(left) /= [k, m]) .or. all(shape(right) /= [n, k]) ) &
                    error stop "intrinsic_mm_driver::perform_mm_c64:Strided mm not implemented."

            if ( beta == 0.0 ) then
                dst = alpha*matmul(transpose(left), transpose(right))
            else
                dst = beta*dst + alpha*matmul(transpose(left), transpose(right))
            end if
        else if (transposed_left) then
            if ( all(shape(left) /= [k, m]) .or. all(shape(right) /= [k, n]) ) &
                    error stop "intrinsic_mm_driver::perform_mm_c64:Strided mm not implemented."

            if ( beta == 0.0 ) then
                dst = alpha*matmul(transpose(left), right)
            else
                dst = beta*dst + alpha*matmul(transpose(left), right)
            end if
        else if ( transposed_right ) then
            if ( all(shape(left) /= [m, k]) .or. all(shape(right) /= [n, k]) ) &
                    error stop "intrinsic_mm_driver::perform_mm_c64:Strided mm not implemented."

            if ( beta == 0.0 ) then
                dst = alpha*matmul(left, transpose(right))
            else
                dst = beta*dst + alpha*matmul(left, transpose(right))
            end if
        else
            if ( all(shape(left) /= [m, k]) .or. all(shape(right) /= [k, n]) ) &
                    error stop "intrinsic_mm_driver::perform_mm_c64:Strided mm not implemented."

            if ( beta == 0.0 ) then
                dst = alpha*matmul(left, right)
            else
                dst = beta*dst + alpha*matmul(left, right)
            end if
        end if
    end subroutine perform_mm_c64

    subroutine perform_mm_c128(this, dst, left, right, m, n, k, &
            alpha, beta, transposed_left, transposed_right, astream)
        class(intrinsic_mm_driver), intent(in) :: this
        complex(real64), dimension(:,:), contiguous, intent(inout) :: dst
        complex(real64), dimension(:,:), contiguous, intent(in) :: left, right
        integer(int64), intent(in) :: m, n, k
        complex(real64), intent(in) :: alpha, beta
        logical, intent(in) :: transposed_left, transposed_right
        type(stream), intent(in), optional :: astream

        if ( all(shape(dst) /= [m, n]) ) &
                error stop "intrinsic_mm_driver::perform_mm_c128:Strided mm not implemented."

        if ( transposed_left .and. transposed_right ) then
            if ( all(shape(left) /= [k, m]) .or. all(shape(right) /= [n, k]) ) &
                    error stop "intrinsic_mm_driver::perform_mm_c128:Strided mm not implemented."

            if ( beta == 0.0d0 ) then
                dst = alpha*matmul(transpose(left), transpose(right))
            else
                dst = beta*dst + alpha*matmul(transpose(left), transpose(right))
            end if
        else if (transposed_left) then
            if ( all(shape(left) /= [k, m]) .or. all(shape(right) /= [k, n]) ) &
                    error stop "intrinsic_mm_driver::perform_mm_c128:Strided mm not implemented."

            if ( beta == 0.0d0 ) then
                dst = alpha*matmul(transpose(left), right)
            else
                dst = beta*dst + alpha*matmul(transpose(left), right)
            end if
        else if ( transposed_right ) then
            if ( all(shape(left) /= [m, k]) .or. all(shape(right) /= [n, k]) ) &
                    error stop "intrinsic_mm_driver::perform_mm_c128:Strided mm not implemented."

            if ( beta == 0.0d0 ) then
                dst = alpha*matmul(left, transpose(right))
            else
                dst = beta*dst + alpha*matmul(left, transpose(right))
            end if
        else
            if ( all(shape(left) /= [m, k]) .or. all(shape(right) /= [k, n]) ) &
                    error stop "intrinsic_mm_driver::perform_mm_c128:Strided mm not implemented."

            if ( beta == 0.0d0 ) then
                dst = alpha*matmul(left, right)
            else
                dst = beta*dst + alpha*matmul(left, right)
            end if
        end if
    end subroutine perform_mm_c128

    subroutine cleanup(this)
        class(intrinsic_mm_driver), intent(inout) :: this

        continue
    end subroutine cleanup

    subroutine synchronize(this, astream)
        class(intrinsic_mm_driver), intent(in) :: this
        type(stream), intent(in), optional :: astream

        continue
    end subroutine synchronize
end module intrinsic_mm_driver_module

max_stars_repo_name: cheshyre/ntcl-algorithms

max_issues_repo_name: cheshyre/ntcl-algorithms

max_forks_repo_name: cheshyre/ntcl-algorithms


text:
      PROGRAM WIND_STAT_RANGE_NC
      USE netcdf   ! NetCDF fortran 90 interface
      IMPLICIT NONE
C
      CHARACTER*40     CTITLE
      INTEGER          IWI,JWI,N,NREC,IOS
      REAL             XFIN,YFIN,DXIN,DYIN
C
      CHARACTER*240    CFILE
      INTEGER          ncFID,ncDID,ncVID,nV,NF
      INTEGER          KREC
      REAL             JDAY,YEAR
      DOUBLE PRECISION TIME(9000),TIME_NEXT
      REAL             WDAY(9001),SCALE_F,ADD_OFF,WMIN,WMAX
C
      INTEGER*2,        ALLOCATABLE :: W(:,:)
      DOUBLE PRECISION, ALLOCATABLE :: WLON(:),WLAT(:)
C
C**********
C*
C 1)  PRINT MODEL WIND FILE STATISTICS, FOR EACH FIELD.
C
C 2)  NRL-style NetCDF WIND FILE FROM ENVIRONEMENT VARIABLE CDF055.
C
C 3)  ALAN J. WALLCRAFT,  JULY 2010.
C*
C**********
C
C     OPEN THE FILE.
C
      CFILE = ' '
      CALL GETENV('CDF055',CFILE)
      IF     (CFILE.EQ.' ') THEN
        WRITE(0,*) 'wind_stat_range_nc: no CDF055 environment variable'
        CALL EXIT(1)
        STOP
      ENDIF
      call nchek('nf90_open',
     &            nf90_open(trim(CFILE), nf90_nowrite, ncFID))
C
C     READ THE NUMBER OF VARIABLES - INFER NUMBER OF FIELDS
C
      call nchek("nf90_inquire-nV",
     &            nf90_inquire(ncFID,nVariables=nV))
      NF = nV - 4  !skip: MT,Date,Latitude,Longitude
      IF     (NF.LE.0) THEN
        WRITE(0,*) 'wind_stat_range_nc: too few variables'
        WRITE(0,*) '  need: MT,Date,Latitude,Longitude,field1,...'
        CALL EXIT(1)
        STOP
      ENDIF
C
C     READ THE TITLE.
C
      call nchek("nf90_get_att-title",
     &            nf90_get_att(ncFID,nf90_global,
     &                               "title",
     &                               CTITLE))
C
C     READ THE AXES.
C
      call nchek('nf90_inq_dimid-MT',
     &            nf90_inq_dimid(        ncFID, 'MT',ncDID))
      call nchek('nf90_inquire_dimension-MT',
     &            nf90_inquire_dimension(ncFID,      ncDID,
     &                                               len=NREC))
      call nchek('nf90_inq_varid-MT',
     &            nf90_inq_varid(ncFID,'MT',ncVID))
      call nchek('nf90_get_var-MT',
     &            nf90_get_var(  ncFID,     ncVID,TIME(1:NREC)))
      call nchek("nf90_get_att-next_MT",
     &            nf90_get_att(ncFID,       ncVID,
     &                                      "next_MT",
     &                                      TIME_NEXT))
C
      call nchek('nf90_inq_dimid-Latitude',
     &            nf90_inq_dimid(        ncFID, 'Latitude',ncDID))
      call nchek('nf90_inquire_dimension-Latitude',
     &            nf90_inquire_dimension(ncFID,            ncDID,
     &                                                 len=JWI))
      ALLOCATE( WLAT(JWI) )
      call nchek('nf90_inq_varid-Latitude',
     &            nf90_inq_varid(ncFID,'Latitude',ncVID))
      call nchek('nf90_get_var-Latitude',
     &            nf90_get_var(  ncFID,           ncVID,WLAT(:)))
C
      call nchek('nf90_inq_dimid-Longitude',
     &            nf90_inq_dimid(        ncFID, 'Longitude',ncDID))
      call nchek('nf90_inquire_dimension-Longitude',
     &            nf90_inquire_dimension(ncFID,             ncDID,
     &                                                  len=IWI))
      ALLOCATE( WLON(IWI) )
      call nchek('nf90_inq_varid-Longitude',
     &            nf90_inq_varid(ncFID,'Longitude',ncVID))
      call nchek('nf90_get_var-Longitude',
     &            nf90_get_var(  ncFID,            ncVID,WLON(:)))
C
      XFIN = WLON(1)
      DXIN = (WLON(IWI) - WLON(1))/(IWI-1.D0)  !assume uniform spaced
      YFIN = WLAT(1)
      IF     (ABS((WLAT(JWI)-WLAT(1))-
     &            (WLAT(  2)-WLAT(1))*(JWI-1.D0)).LT.1.D-2) THEN
*       write(6,*)     (WLAT(JWI)-WLAT(1))
*       write(6,*)     (WLAT(  2)-WLAT(1))*(JWI-1.D0)
*       write(6,*) ABS((WLAT(JWI)-WLAT(1))-
*    &                 (WLAT(  2)-WLAT(1))*(JWI-1.D0))
        DYIN = (WLAT(JWI) - WLAT(1))/(JWI-1.D0)  !assume uniform spaced
      ELSE
        DYIN = 0.0  ! assume a gaussian grid
      ENDIF
      DO KREC= 1,NREC
        WDAY(KREC) = TIME(KREC)
      ENDDO
      WDAY(NREC+1) = TIME_NEXT
C
C     STATISTICS.
C
      WRITE(6,6000) CTITLE
      WRITE(6,6100) IWI,JWI,XFIN,YFIN,DXIN,DYIN,
     +              NREC
*     CALL FLUSH(6)
C
C     READ NREC WIND RECORDS, AND PRINT RANGES.
C
      ALLOCATE( W(IWI,JWI), STAT=IOS )
      IF     (IOS.NE.0) THEN
        WRITE(6,*) 'Error in wind_stat_range: could not allocate ',
     +             IWI*JWI,' 2-byte words'
        CALL EXIT(2)
      ENDIF
C
      W(:,:) = 0.0
      DO KREC= 1,NREC
        DO N= 1,NF
          ncVID = nV-NF+N
          call nchek('nf90_get_var-ncVID',
     &                nf90_get_var(ncFID, ncVID,
     &                                    W(:,:),
     &                                    (/ 1,1,KREC /) ))
          call nchek("nf90_get_att-scale_factor",
     &                nf90_get_att(ncFID, ncVID,
     &                                    "scale_factor",
     &                                    SCALE_F))
          call nchek("nf90_get_att-add_offset",
     &                nf90_get_att(ncFID, ncVID,
     &                                    "add_offset",
     &                                    ADD_OFF))

          WMIN = SCALE_F*MINVAL(W(:,:)) + ADD_OFF
          WMAX = SCALE_F*MAXVAL(W(:,:)) + ADD_OFF
          WRITE(6,'(f10.3,3x,1p2e16.8)') WDAY(KREC),WMIN,WMAX
        ENDDO !n
      ENDDO !krec
C
C     SUMMARY.
C
      CALL WNDAY(WDAY(1), YEAR,JDAY)
      IF     (YEAR.LT.1904.5) THEN
        WRITE(6,6200) NREC,JDAY,NINT(YEAR),WDAY(NREC+1)-WDAY(1)
      ELSE
        WRITE(6,6250) NREC,JDAY,NINT(YEAR),WDAY(NREC+1)-WDAY(1)
      ENDIF
      CALL EXIT(0)
      STOP
C
 6000 FORMAT(A40)
 6100 FORMAT(
     +      'IWI,JWI =',I4,',',I4,
     +   3X,'XFIN,YFIN =',F9.3,',',F9.3,
     +   3X,'DXIN,DYIN =',F6.3,',',F6.3 /
     +      'NREC =',I5,5X,'WDAY =' / (8F10.3) )
 6200 FORMAT(I5,' RECORD CLIMATOLOGY STARTING ON',F7.2,'/',I4,
     +   ' COVERING',F9.2,' DAYS')
 6250 FORMAT(I5,' WIND RECORDS STARTING ON',F7.2,'/',I4,
     +   ' COVERING',F9.2,' DAYS')
C     END OF WDSTAT.
      END
      SUBROUTINE WNDAY(WDAY, YEAR,DAY)
      IMPLICIT NONE
      REAL WDAY,YEAR,DAY
C
C**********
C*
C  1) CONVERT 'WIND DAY' INTO JULIAN DAY AND YEAR.
C
C  2) THE 'WIND DAY' IS THE NUMBER OF DAYS SINCE 001/1901 (WHICH IS
C      WIND DAY 1.0).
C     FOR EXAMPLE:
C      A) YEAR=1901.0 AND DAY=1.0, REPRESENTS 0000Z HRS ON 001/1901
C         SO WDAY WOULD BE 1.0.
C      B) YEAR=1901.0 AND DAY=2.5, REPRESENTS 1200Z HRS ON 002/1901
C         SO WDAY WOULD BE 2.5.
C     YEAR MUST BE NO LESS THAN 1901.0, AND NO GREATER THAN 2099.0.
C     NOTE THAT YEAR 2000 IS A LEAP YEAR (BUT 1900 AND 2100 ARE NOT).
C
C  3) ALAN J. WALLCRAFT, PLANNING SYSTEMS INC., FEBRUARY 1993.
C*
C**********
C
      INTEGER IYR,NLEAP
      REAL    WDAY1
C
C     FIND THE RIGHT YEAR.
C
      IYR   = (WDAY-1.0)/365.25
      NLEAP = IYR/4
      WDAY1 = 365.0*IYR + NLEAP + 1.0
      DAY   = WDAY - WDAY1 + 1.0
      IF     (WDAY1.GT.WDAY) THEN
        IYR   = IYR - 1
      ELSEIF (DAY.GE.367.0) THEN
        IYR   = IYR + 1
      ELSEIF (DAY.GE.366.0 .AND. MOD(IYR,4).NE.3) THEN
        IYR   = IYR + 1
      ENDIF
      NLEAP = IYR/4
      WDAY1 = 365.0*IYR + NLEAP + 1.0
C
C     RETURN YEAR AND JULIAN DAY.
C
      YEAR = 1901 + IYR
      DAY  = WDAY - WDAY1 + 1.0
      RETURN
C     END OF WNDAY.
      END
      subroutine nchek(cnf90,status)
      use netcdf   ! NetCDF fortran 90 interface
      implicit none
c
      character*(*), intent(in) :: cnf90
      integer,       intent(in) :: status
c
c     subroutine to handle NetCDF errors
c
*     if     (.TRUE. ) then !debug
      if     (.FALSE.) then !nodebug
        write(6,'(a)') trim(cnf90)
      endif

      if (status /= nf90_noerr) then
        write(6,'(/a)')   'error in profout - from NetCDF library'
        write(6,'(a/)')   trim(cnf90)
        write(6,'(a/)')   trim(nf90_strerror(status))
        stop
      end if
      end subroutine nchek

ext: f

max_stars_repo_name: TillRasmussen/HYCOM-tools

max_issues_repo_name: TillRasmussen/HYCOM-tools

max_forks_repo_name: TillRasmussen/HYCOM-tools


text:
Picnic Day 2007 is this years Picnic Day celebration with the motto Making OUR MarkFootnote(The Logo is used on this entry by Fair Use Doctrine, and is not Creative Commons).  As usual, there will be the 150+ exhibits and special events all over campus, the parade, Picnic Day Battle of the Bands, a sledding hill with Ski Demonstrations, Davis Dance Revolution, and so much more! This year, we are also hosting a nonprofit faire at Olson Lawn and we are hoping to showcase the charities and nonprofits in the Yolo/Sacramento area.

 Picnic Day Preweek
  Saturday, April 7, 2007:  Featured Creatures will be held at Davis Farmers Market at 10am! It will consist of our second pet costume contest, and our first ever Hamster Ball Races.  This event is cosponsored by Petco, who will provide the hamster ball racetrack. Come kick off Picnic Day Preweek with this fun event!
  Monday, April 9, 2007: Entertainment Preview, noon, Quad
  Tuesday, April 10, 2007: 1st Year Students Talent Show, 8pm, Regan Main
  Wednesday, April 11, 2007: Annual Cow Milking Contest, noon, Quad
  Friday, April 13, 2007: Directors Hay Ride, noon, all around campus

Post your photos below!  Just click Edit on this page and then click on the Upload Files button.

 Photos

20070410 20:42:08 nbsp So, at the moment, weather.com is predicting showers for this Saturday... Users/PaulAmnuaypayoat

20070412 20:39:49 nbsp Weather.com is now predicting rain all day, with a 90% chance of precipitation.  Looks like we might get a rainout like happened a few years ago. Users/AlexPomeranz

  Weather.com is quite unreliable for this area, in my experience.  It has Woodland right now at 40% chance of showers for this afternoon (full of blazing sunshine) and 60% chance of showers tomorrow.

20070413 08:35:34 nbsp Better buy your umbrellas now, if you dont have one. When it poured in 2003, I wasnt prepared and wandered around campus and downtown for hours in the rain. Most stores sold out of their umbrellas, so I had to pay a lot of money for one at Gottschalks Users/ElleWeber

20070413 09:42:52 nbsp have faith young padawans. The sprit of picnic day will clear the skies and the sun will come out yet! Or well get wet. Users/TomAbeles

20070414 11:12:04 nbsp Its raining now... not worth it to go out as Ive missed the parade already. Users/TusharRawat

20070414 16:55:29 nbsp Are there any official stats of what percentage of picnic days get rained out? 3 out of the past 5 years doesnt seem very good... Users/AaGg

    Lets make a list!  2004, 2003, 2003?  It didnt rain in 2005 or 2006..

        So, 2003, 2004, 2007, thats not 3/5...?
              Wouldnt say 04 was rained out completely, though it did rain and cancel some things partially. And it didnt rain in 02.  So thats 2(+.5)  / 6 I guess.

20070414 17:04:13 nbsp I didnt take photos from the float this year due to the rain, but I saw a lot of cameras in the crowd. There was a delay in the parade, does anyone know what caused it? Users/JasonAller

20070414 17:05:56 nbsp 2003 was completely rained out, pretty much like today.  04 rained, but it wasnt that bad, just a drizzle on and off.  They didnt cancel everything outdoor in 2004, a lot of bands still played outside or in those covered trailor stages, iirc.  Users/EdWins

20070414 19:50:24 nbsp The delay toward the beginning of the parade was just because the walking groups were getting shuffled in from a location separate from where the floats started. Or at least that was my understanding from the limited vantage point inside the KDVS blimp gondola. Users/ElisaHough

20070414 22:50:03 nbsp Nice work with the photos and everything you guys.  I was cooped up in the bar all day so know I feel like I was really there  Users/JarrettNoble

20070415 11:54:03 nbsp so rainy, all that was left to do was drink indoors Users/mynameisbernice1

20070415 17:22:17 nbsp I was driving the double deck in the parade for Unitrans and after I had the hook and ladder fire truck pulled in front of me from the walking groups location we still got stuck way down on California forever. I was just glad I was in a bus. The campus tours people(?) who were behind me got soaked for the most part, but they did get to hang out on my bus while we were delayed. Users/TracyPerkins

20070415 20:56:08 nbsp Ive always wondered why they have Picnic Day in April (April showers) which is prime rainy season around here.  Seems like May would be a much better choice  after a lot of the rain, but before it get baking hot. Users/AlexPomeranz

     They might want to not have to compete with the Whole Earth Festival.

20070415 22:10:37 nbsp I heard the KDVS blimp escaped and the FAA was not pleased. Confirm/Deny? Users/ArlenAbraham
 The KDVS blimp was cut loose and Drake (GM) called the FAA, but they didnt seem particularly concerned. Users/ArlenAbraham

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
C
C	$Id: gflas2.f,v 1.5 2008-07-27 00:17:13 haley Exp $
C
C                Copyright (C)  2000
C        University Corporation for Atmospheric Research
C                All Rights Reserved
C
C The use of this Software is governed by a License Agreement.
C
      SUBROUTINE GFLAS2
C
C_BEGIN PROLOGUE GFLAS2
C
C_PURPOSE
C    A call to this subroutine flags termination of output of
C    plotting instructions to Workstation Independent Segment
C    Storage as initiated by a previous call to GFLAS1.
C
C_DESCRIPTION
C    A call to this subroutine flags termination of output of
C    plotting instructions to Workstation Independent Segment
C    Storage (WISS) as initiated by a previous call to GFLAS1.
C
C_ARGUMENTS
C
C (None)
C
C_I/O
C    GFLAS2 terminates output to WISS.  If using the NCAR GKS
C    package, the output file which was opened with GFLAS1 is
C    closed.
C
C_ROUTINES CALLED
C    PLOTIT, SETER
C
C_COMMON BLOCKS
C    GFLASH
C
C_MAXIMUM GKS LEVEL
C    2A
C
C_LANGUAGE
C    FORTRAN
C
C_REVISION DATE
C    890117 (YYMMDD)
C
C_HISTORY
C    This subroutine was coded at NCAR by Fred Clare.  Its purpose
C    is to duplicate the functionality of the previous NCAR System
C    Plot Package entry FLASH2 in a GKS environment.
C
C_REFERENCES
C    "The System Plot Package", NCAR Technical Note number
C    NCAR-TN/162+IA, May, 1983.
C
C    "NCAR Graphics User's Guide", Version 2.00, NCAR Technical Note
C    number NCAR/TN-283+IA, August, 1987.
C
C_LONG DESCRIPTION
C    GFLAS2 deactivates WISS and re-activates all workstations
C    which were active before the call to GFLAS1 which initiated
C    the current GFLAS1/GFLAS2 sequence.  This results in all
C    subsequent plotting instructions being sent to all workstations
C    which were active prior to the current GFLAS1/GFLAS2 sequence.
C    If one is using the NCAR GKS package, the output file GFBNnn
C    is closed (where "nn" is the value of the argument to GFLAS1
C    in the current GFLAS1/GFLAS2 sequence.)
C
C    It is not legal to call FRAME (or clear workstations) between a
C    GFLAS1 call and a GFLAS2 call.
C
C_EXAMPLES
C    Example 1 --  Save a text string in a flash buffer (see GFLAS1
C                  for usage of that subroutine.)
C          SUBROUTINE SVTXT
C          CALL GFLAS1(9)
C          CALL GTX(.5,.5,'STRING')
C          CALL GFLAS2
C          RETURN
C          END
C
C_END PROLOGUE GFLAS2
C
C-------------------------------------------------------------------------
C
C  Common for FLASH package.
C
      COMMON /GFLASH/MODEF,IOPWKS(100),IOACT(100),NUMOP,IWISSI
C     CHARACTER*80 IDR(1),ODR(1)
C
      SAVE
C
C  Check on proper mode value.
C
      IF (MODEF .NE. 1) GO TO 101
      MODEF = 2
C
C  Flush the PLOTIT buffer.
C
      CALL PLOTIT(0,0,0)
C
      CALL GCLSG
C
C  Restore workstation states to what they were before call to
C  GFLAS1.  If WISS was not active, deactivate it since GFLAS1
C  made it active; reactivate all other previously active
C  workstations.
C
C
C  Determine the index of WISS.
C
      DO 20 I=1,NUMOP
      IF (IOPWKS(I) .EQ. IWISSI) THEN
        NWISS = I
        GO TO 30
      ENDIF
   20 CONTINUE
      GO TO 102
   30 CONTINUE
C
C  Deactive WISS if it was not active before call to GFLAS1.
C
      IF (IOACT(NWISS) .EQ. 0) CALL GDAWK(IWISSI)
C
C  Handle non-WISS workstations.
C
      DO 10 I=1,NUMOP
      IF (IOPWKS(I).NE.IWISSI .AND. IOACT(I).EQ.1) CALL GACWK(IOPWKS(I))
   10 CONTINUE
C
C  Inform NCAR GKS that we are finished in GFLAS2
C
C     IDR(1) = '0'
C     IFID = -1394
C     CALL PLOTIT(0,0,2)
C     CALL GESC(IFID,1,IDR,1,IDUM,ODR)
C
      RETURN
C
  101 CONTINUE
      CALL GECLKS
      CALL SETER(' GFLAS2 -- GFLAS2 CALLED WITHOUT CALLING GFLAS1',6,2)
      STOP
  102 CONTINUE
      CALL GECLKS
      CALL SETER(' GFLAS2 -- WISS NOT OPEN AT GFLAS2 TIME',8,2)
      STOP
C
      END

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


text:
      SUBROUTINE MG5_1_MP_HELAS_CALLS_AMPB_1(P,NHEL,H,IC)
C
      USE MG5_1_POLYNOMIAL_CONSTANTS
      IMPLICIT NONE
C
C     CONSTANTS
C
      INTEGER    NEXTERNAL
      PARAMETER (NEXTERNAL=4)
      INTEGER    NCOMB
      PARAMETER (NCOMB=16)

      INTEGER NBORNAMPS
      PARAMETER (NBORNAMPS=2)
      INTEGER    NLOOPS, NLOOPGROUPS, NCTAMPS
      PARAMETER (NLOOPS=22, NLOOPGROUPS=13, NCTAMPS=54)
      INTEGER    NLOOPAMPS
      PARAMETER (NLOOPAMPS=76)
      INTEGER    NWAVEFUNCS,NLOOPWAVEFUNCS
      PARAMETER (NWAVEFUNCS=8,NLOOPWAVEFUNCS=50)
      REAL*16     ZERO
      PARAMETER (ZERO=0.0E0_16)
      COMPLEX*32     IZERO
      PARAMETER (IZERO=CMPLX(0.0E0_16,0.0E0_16,KIND=16))
C     These are constants related to the split orders
      INTEGER    NSO, NSQUAREDSO, NAMPSO
      PARAMETER (NSO=1, NSQUAREDSO=1, NAMPSO=2)
C
C     ARGUMENTS
C
      REAL*16 P(0:3,NEXTERNAL)
      INTEGER NHEL(NEXTERNAL), IC(NEXTERNAL)
      INTEGER H
C
C     LOCAL VARIABLES
C
      INTEGER I,J,K
      COMPLEX*32 COEFS(MAXLWFSIZE,0:VERTEXMAXCOEFS-1,MAXLWFSIZE)
C
C     GLOBAL VARIABLES
C
      INCLUDE 'mp_coupl_same_name.inc'

      INTEGER GOODHEL(NCOMB)
      LOGICAL GOODAMP(NSQUAREDSO,NLOOPGROUPS)
      COMMON/MG5_1_FILTERS/GOODAMP,GOODHEL

      INTEGER SQSO_TARGET
      COMMON/MG5_1_SOCHOICE/SQSO_TARGET

      LOGICAL UVCT_REQ_SO_DONE,MP_UVCT_REQ_SO_DONE,CT_REQ_SO_DONE
     $ ,MP_CT_REQ_SO_DONE,LOOP_REQ_SO_DONE,MP_LOOP_REQ_SO_DONE
     $ ,CTCALL_REQ_SO_DONE,FILTER_SO
      COMMON/MG5_1_SO_REQS/UVCT_REQ_SO_DONE,MP_UVCT_REQ_SO_DONE
     $ ,CT_REQ_SO_DONE,MP_CT_REQ_SO_DONE,LOOP_REQ_SO_DONE
     $ ,MP_LOOP_REQ_SO_DONE,CTCALL_REQ_SO_DONE,FILTER_SO

      COMPLEX*32 AMP(NBORNAMPS)
      COMMON/MG5_1_MP_AMPS/AMP
      COMPLEX*32 W(20,NWAVEFUNCS)
      COMMON/MG5_1_MP_W/W

      COMPLEX*32 WL(MAXLWFSIZE,0:LOOPMAXCOEFS-1,MAXLWFSIZE,
     $ -1:NLOOPWAVEFUNCS)
      COMPLEX*32 PL(0:3,-1:NLOOPWAVEFUNCS)
      COMMON/MG5_1_MP_WL/WL,PL

      COMPLEX*32 AMPL(3,NCTAMPS)
      COMMON/MG5_1_MP_AMPL/AMPL

C
C     ----------
C     BEGIN CODE
C     ----------

C     The target squared split order contribution is already reached
C      if true.
      IF (FILTER_SO.AND.MP_CT_REQ_SO_DONE) THEN
        GOTO 1001
      ENDIF

      CALL MP_IXXXXX(P(0,1),ZERO,NHEL(1),+1*IC(1),W(1,1))
      CALL MP_OXXXXX(P(0,2),ZERO,NHEL(2),-1*IC(2),W(1,2))
      CALL MP_OXXXXX(P(0,3),ZERO,NHEL(3),+1*IC(3),W(1,3))
      CALL MP_IXXXXX(P(0,4),ZERO,NHEL(4),-1*IC(4),W(1,4))
      CALL MP_FFV1P0_3(W(1,1),W(1,2),GC_5,ZERO,ZERO,W(1,5))
C     Amplitude(s) for born diagram with ID 1
      CALL MP_FFV1_0(W(1,4),W(1,3),W(1,5),GC_5,AMP(1))
      CALL MP_FFV1P0_3(W(1,1),W(1,3),GC_5,ZERO,ZERO,W(1,6))
C     Amplitude(s) for born diagram with ID 2
      CALL MP_FFV1_0(W(1,4),W(1,2),W(1,6),GC_5,AMP(2))
      CALL MP_FFV1P0_3(W(1,4),W(1,3),GC_5,ZERO,ZERO,W(1,7))
C     Counter-term amplitude(s) for loop diagram number 3
      CALL MP_R2_GG_1_0(W(1,5),W(1,7),R2_GGQ,AMPL(1,1))
      CALL MP_R2_GG_1_0(W(1,5),W(1,7),R2_GGQ,AMPL(1,2))
      CALL MP_R2_GG_1_0(W(1,5),W(1,7),R2_GGQ,AMPL(1,3))
      CALL MP_R2_GG_1_0(W(1,5),W(1,7),R2_GGQ,AMPL(1,4))
      CALL MP_FFV1P0_3(W(1,4),W(1,2),GC_5,ZERO,ZERO,W(1,8))
C     Counter-term amplitude(s) for loop diagram number 4
      CALL MP_R2_GG_1_0(W(1,6),W(1,8),R2_GGQ,AMPL(1,5))
      CALL MP_R2_GG_1_0(W(1,6),W(1,8),R2_GGQ,AMPL(1,6))
      CALL MP_R2_GG_1_0(W(1,6),W(1,8),R2_GGQ,AMPL(1,7))
      CALL MP_R2_GG_1_0(W(1,6),W(1,8),R2_GGQ,AMPL(1,8))
C     Counter-term amplitude(s) for loop diagram number 5
      CALL MP_FFV1_0(W(1,4),W(1,3),W(1,5),UV_GQQB_1EPS,AMPL(2,9))
      CALL MP_FFV1_0(W(1,4),W(1,3),W(1,5),UV_GQQB_1EPS,AMPL(2,10))
      CALL MP_FFV1_0(W(1,4),W(1,3),W(1,5),UV_GQQB_1EPS,AMPL(2,11))
      CALL MP_FFV1_0(W(1,4),W(1,3),W(1,5),UV_GQQB_1EPS,AMPL(2,12))
      CALL MP_FFV1_0(W(1,4),W(1,3),W(1,5),UV_GQQB,AMPL(1,13))
      CALL MP_FFV1_0(W(1,4),W(1,3),W(1,5),UV_GQQB_1EPS,AMPL(2,14))
      CALL MP_FFV1_0(W(1,4),W(1,3),W(1,5),UV_GQQT,AMPL(1,15))
      CALL MP_FFV1_0(W(1,4),W(1,3),W(1,5),UV_GQQB_1EPS,AMPL(2,16))
      CALL MP_FFV1_0(W(1,4),W(1,3),W(1,5),UV_GQQG_1EPS,AMPL(2,17))
      CALL MP_FFV1_0(W(1,4),W(1,3),W(1,5),R2_GQQ,AMPL(1,18))
C     Counter-term amplitude(s) for loop diagram number 7
      CALL MP_FFV1_0(W(1,4),W(1,2),W(1,6),UV_GQQB_1EPS,AMPL(2,19))
      CALL MP_FFV1_0(W(1,4),W(1,2),W(1,6),UV_GQQB_1EPS,AMPL(2,20))
      CALL MP_FFV1_0(W(1,4),W(1,2),W(1,6),UV_GQQB_1EPS,AMPL(2,21))
      CALL MP_FFV1_0(W(1,4),W(1,2),W(1,6),UV_GQQB_1EPS,AMPL(2,22))
      CALL MP_FFV1_0(W(1,4),W(1,2),W(1,6),UV_GQQB,AMPL(1,23))
      CALL MP_FFV1_0(W(1,4),W(1,2),W(1,6),UV_GQQB_1EPS,AMPL(2,24))
      CALL MP_FFV1_0(W(1,4),W(1,2),W(1,6),UV_GQQT,AMPL(1,25))
      CALL MP_FFV1_0(W(1,4),W(1,2),W(1,6),UV_GQQB_1EPS,AMPL(2,26))
      CALL MP_FFV1_0(W(1,4),W(1,2),W(1,6),UV_GQQG_1EPS,AMPL(2,27))
      CALL MP_FFV1_0(W(1,4),W(1,2),W(1,6),R2_GQQ,AMPL(1,28))
C     Counter-term amplitude(s) for loop diagram number 9
      CALL MP_FFV1_0(W(1,1),W(1,3),W(1,8),UV_GQQB_1EPS,AMPL(2,29))
      CALL MP_FFV1_0(W(1,1),W(1,3),W(1,8),UV_GQQB_1EPS,AMPL(2,30))
      CALL MP_FFV1_0(W(1,1),W(1,3),W(1,8),UV_GQQB_1EPS,AMPL(2,31))
      CALL MP_FFV1_0(W(1,1),W(1,3),W(1,8),UV_GQQB_1EPS,AMPL(2,32))
      CALL MP_FFV1_0(W(1,1),W(1,3),W(1,8),UV_GQQB,AMPL(1,33))
      CALL MP_FFV1_0(W(1,1),W(1,3),W(1,8),UV_GQQB_1EPS,AMPL(2,34))
      CALL MP_FFV1_0(W(1,1),W(1,3),W(1,8),UV_GQQT,AMPL(1,35))
      CALL MP_FFV1_0(W(1,1),W(1,3),W(1,8),UV_GQQB_1EPS,AMPL(2,36))
      CALL MP_FFV1_0(W(1,1),W(1,3),W(1,8),UV_GQQG_1EPS,AMPL(2,37))
      CALL MP_FFV1_0(W(1,1),W(1,3),W(1,8),R2_GQQ,AMPL(1,38))
C     Counter-term amplitude(s) for loop diagram number 13
      CALL MP_FFV1_0(W(1,1),W(1,2),W(1,7),UV_GQQB_1EPS,AMPL(2,39))
      CALL MP_FFV1_0(W(1,1),W(1,2),W(1,7),UV_GQQB_1EPS,AMPL(2,40))
      CALL MP_FFV1_0(W(1,1),W(1,2),W(1,7),UV_GQQB_1EPS,AMPL(2,41))
      CALL MP_FFV1_0(W(1,1),W(1,2),W(1,7),UV_GQQB_1EPS,AMPL(2,42))
      CALL MP_FFV1_0(W(1,1),W(1,2),W(1,7),UV_GQQB,AMPL(1,43))
      CALL MP_FFV1_0(W(1,1),W(1,2),W(1,7),UV_GQQB_1EPS,AMPL(2,44))
      CALL MP_FFV1_0(W(1,1),W(1,2),W(1,7),UV_GQQT,AMPL(1,45))
      CALL MP_FFV1_0(W(1,1),W(1,2),W(1,7),UV_GQQB_1EPS,AMPL(2,46))
      CALL MP_FFV1_0(W(1,1),W(1,2),W(1,7),UV_GQQG_1EPS,AMPL(2,47))
      CALL MP_FFV1_0(W(1,1),W(1,2),W(1,7),R2_GQQ,AMPL(1,48))
C     Counter-term amplitude(s) for loop diagram number 17
      CALL MP_R2_GG_1_R2_GG_3_0(W(1,5),W(1,7),R2_GGQ,R2_GGB,AMPL(1,49))
C     Counter-term amplitude(s) for loop diagram number 18
      CALL MP_R2_GG_1_R2_GG_3_0(W(1,6),W(1,8),R2_GGQ,R2_GGB,AMPL(1,50))
C     Counter-term amplitude(s) for loop diagram number 19
      CALL MP_R2_GG_1_R2_GG_3_0(W(1,5),W(1,7),R2_GGQ,R2_GGT,AMPL(1,51))
C     Counter-term amplitude(s) for loop diagram number 20
      CALL MP_R2_GG_1_R2_GG_3_0(W(1,6),W(1,8),R2_GGQ,R2_GGT,AMPL(1,52))
C     Counter-term amplitude(s) for loop diagram number 21
      CALL MP_R2_GG_1_R2_GG_2_0(W(1,5),W(1,7),R2_GGG_1,R2_GGG_2,AMPL(1
     $ ,53))
C     Counter-term amplitude(s) for loop diagram number 22
      CALL MP_R2_GG_1_R2_GG_2_0(W(1,6),W(1,8),R2_GGG_1,R2_GGG_2,AMPL(1
     $ ,54))
C     At this point, all CT amps needed for (QCD=6), i.e. of split
C      order ID=1, are computed.
      IF(FILTER_SO.AND.SQSO_TARGET.EQ.1) GOTO 2000

      GOTO 1001
 2000 CONTINUE
      MP_CT_REQ_SO_DONE=.TRUE.
 1001 CONTINUE
      END


text:
C --- Driver for computing sensitivity coefficients w.r.t. all initial values
      PROGRAM driver

      INCLUDE 'KPP_ROOT_Parameters.h'
      INCLUDE 'KPP_ROOT_Global.h'

C ---  Number of sensitivity coefficients to compute
C ---  Note: this value is set for sensitivities w.r.t. all initial values
C ---       it may have to be changed for other applications
      INTEGER NSENSIT
      PARAMETER (NSENSIT = NVAR)

      KPP_REAL  DVAL(NSPEC)
      KPP_REAL  Y(NVAR*(NSENSIT+1))
C ---  Note: Y contains: (1:NVAR) concentrations, followed by
C ---                   (1:NVAR) sensitivities w.r.t. first parameter, followed by
C ---                   etc.,  followed by
C ---                   (1:NVAR) sensitivities w.r.t. NSENSIT's parameter

      INTEGER i

C --- The type of sensitivity coefficients to compute
C --- DDMTYPE = 0 : sensitivities w.r.t. initial values
C --- DDMTYPE = 1 : sensitivities w.r.t. parameters
      DDMTYPE = 0

C ---- TIME VARIABLES ------------------

      TSTART = 0
      TEND = TSTART + 600
      DT = 60.
      TEMP = 298

      STEPMIN = 0.01
      STEPMAX = 900

      PRINT*,'Please provide: RTOL = , ATOL = '
      READ*,RTOLS, ATOLS
      do i=1,NVAR
        RTOL(i) = RTOLS
        ATOL(i) = ATOLS
      end do

C ********** TIME LOOP *************************

      CALL Initialize()

C -- Initialize Concentrations and Sensitivities
      DO i=1,NVAR
         Y(i) = VAR(i)
      END DO

C ---  Note: the initial values below are for sensitivities w.r.t. initial values;
C ---       they may have to be changed for other applications
      DO j=1,NSENSIT
        DO i=1,NVAR
	   Y(i+NVAR*j) = 0.0D0
	END DO
	Y(j+NVAR*j) = 1.0D0
      END DO

      CALL InitSaveData()

      WRITE(6,990) (SPC_NAMES(MONITOR(i)), i=1,NMONITOR),
     *             (SMASS(i), i=1,NMASS )
990   FORMAT('done[%] Time[h] ',20(4X,A6))

      TIME = TSTART
      DO WHILE (TIME .lt. TEND)

        CALL GetMass( C, DVAL )
        WRITE(6,991) (TIME-TSTART)/(TEND-TSTART)*100, TIME/3600.,
     *               (C(MONITOR(i))/CFACTOR, i=1,NMONITOR),
     *               (DVAL(i)/CFACTOR, i=1,NMASS)
991   FORMAT(F6.1,'% ',F7.2,3X,20(E10.4,2X))

        CALL SaveData()

        CALL Update_SUN()
        CALL Update_RCONST()

        CALL INTEGRATE( NSENSIT, Y, TIME, TIME+DT )

        DO i=1,NVAR
           VAR(i) = Y(i)
        END DO

      END DO

      CALL GetMass( C, DVAL )
      WRITE(6,991) (TIME-TSTART)/(TEND-TSTART)*100, TIME/3600.,
     *               (C(MONITOR(i))/CFACTOR, i=1,NMONITOR),
     *               (DVAL(i)/CFACTOR, i=1,NMASS)

C      DO i=1,NSENSIT
C        WRITE(6,992) i, ( Y(NVAR*i+j), j=1,NVAR )
C      END DO
C 992  FORMAT('SEN(',I3,') = ',1000(E10.4,2X))

      CALL SaveData()

C *********** END TIME LOOP ********
      OPEN(20, FILE='KPP_ROOT_results.m')
      WRITE(6,*) '**************************************************'
      WRITE(6,*) '* Concentrations and Sensitivities at final time *'
      WRITE(6,*) '*  were written in the file KPP_ROOT_results.m  *'
      WRITE(6,*) '**************************************************'
      DO i=0,NSENSIT
        WRITE(20,993) ( Y(NVAR*i+j), j=1,NVAR )
      END DO
 993  FORMAT(1000(E24.16,2X))

      CALL CloseSaveData()

      STOP
      END

max_stars_repo_name: matzegoebel/WRF-fluxavg

max_issues_repo_name: teb-model/wrf-teb

max_forks_repo_name: teb-model/wrf-teb


text:
      program curden
c***********************************************************************
c
c     daresbury laboratory dl_poly program for the calculation
c     of the fourier transform of the particle current density
c     in space and time
c
c     original curden program written by w.smith march 1982
c     adapted for dl_poly by w.smith august 1994
c     parallel version
c
c     copyright daresbury laboratory 1994
c     author w smith 1994
c
c***********************************************************************
c
      implicit real*8(a-h,o-z)
      parameter (mcore=116000,msp=1080,iread=5,irite=6)
      parameter (pi=3.141592653589793d0,tpi=6.283185307179586d0)
      logical lden,lcor,ltim,lchg
      dimension xxx(msp),yyy(msp),zzz(msp)
      dimension vxx(msp),vyy(msp),vzz(msp)
      dimension man(10),space(mcore)
      data man/10*0/
      call timchk(0,wtime)

c     determine parallel machine characteristics

      open(iread,file='curden_input')
      open(irite,file='curden_output')
      call machine(idnode,mxnode)
      if(idnode.eq.0)then

         write(irite,'(a)')
     x   '# DL_POLY Current Density Correlation Program'

      endif

c     read in control variables

      call start
     x   (idnode,mxnode,natm,ncon,kmax,ntime,ngap,lden,
     x   lcor,ltim,lchg,tstep)

      klim=((2*kmax+1)**3-1)/2

c     calculate fourier transform of current density

      if(lden)then
        iman=(natm+mxnode-1)/mxnode
        man(1)=1
        man(2)=man(1)+2*iman
        man(3)=man(2)+2*iman
        man(4)=man(3)+2*iman
        man(5)=man(4)+4*klim
        man(6)=man(5)+8*klim
        man(7)=man(6)+8*klim
        man(8)=man(7)+2*iman*(kmax+1)
        man(9)=man(8)+2*iman*(kmax+1)
        man(10)=man(9)+2*iman*(kmax+1)
        call fourdc
     x     (idnode,mxnode,natm,ncon,kmax,lchg,xxx,yyy,zzz,vxx,vyy,vzz,
     x      space(man(1)),space(man(2)),space(man(3)),space(man(4)),
     x      space(man(5)),space(man(6)),space(man(7)),space(man(8)),
     X      space(man(9)))
      endif

c     calculate current density correlation function

      if(lcor)then
        nkt=(klim+mxnode-1)/mxnode
        man(1)=1
        man(2)=man(1)+ntime
        man(3)=man(2)+ntime
        man(4)=man(3)+3*klim
        man(5)=man(4)+8*nkt*ntime
        man(6)=man(5)+8*nkt*ntime
        man(7)=man(6)+8*klim
        call correl
     x     (idnode,mxnode,ncon,kmax,ntime,ngap,tstep,space(man(1)),
     x      space(man(2)),space(man(3)),space(man(4)),space(man(5)),
     x      space(man(6)))
      endif

c     calculate correlation function fourier transform

      if(ltim)then
        nkt=(klim+mxnode-1)/mxnode
        man(1)=1
        man(2)=man(1)+ntime
        man(3)=man(2)+3*klim
        man(4)=man(3)+ntime
        man(5)=man(4)+4*nkt*ntime
        man(6)=man(5)+4*ntime
        man(7)=man(6)+4*ntime
        call curfft
     x     (idnode,mxnode,kmax,ntime,ngap,tstep,space(man(1)),
     x      space(man(2)),space(man(3)),space(man(4)),space(man(5)),
     x      space(man(6)))
      endif
      call exit()
      end
      subroutine start
     x   (idnode,mxnode,natm,ncon,kmax,ntime,ngap,lden,lcor,
     x   ltim,lchg,tstep)

c***********************************************************************
c
c     read in control variables and check parameters
c
c     copyright daresbury laboratory 1994
c     author w smith 1994
c
c***********************************************************************
      implicit real*8(a-h,o-z)
      parameter (mcore=116000,msp=1080)
      parameter (iread=5,irite=6)
      logical lden,ltim,lcor,lchg,kill
      character*80 title
      kill=.false.
      read(iread,'(a80)')title
      read(iread,*)natm
      read(iread,*)ncon
      read(iread,*)kmax
      read(iread,*)lden
      read(iread,*)lcor
      read(iread,*)ltim
      read(iread,*)lchg
      read(iread,*)ntime
      read(iread,*)ngap
      read(iread,*)tstep
      if(idnode.eq.0)then
         write(irite,'(a,i10)')
     x     '# maximum size of dynamic core area         = ',mcore
         write(irite,'(a,i10)')
     x     '# number of atoms in simulation box         = ',natm
         write(irite,'(a,i10)')
     x     '# number of configurations                  = ',ncon
         write(irite,'(a,i10)')
     x     '# largest k vector index                    = ',kmax
         write(irite,'(a,i10)')
     x     '# length of correlation arrays              = ',ntime
         write(irite,'(a,i10)')
     x     '# interval between time origins             = ',ngap
         write(irite,'(a,1p,e12.4)')
     x     '# magnitude of time interval                = ',tstep
      endif
      if(idnode.eq.0)then

         if(lchg)then
            write(irite,'(a)')'# charge current option selected'
         else
            write(irite,'(a)')'# density current option selected'
         endif

      endif

c     check on control parameters

      if(natm.gt.msp)then

         if(idnode.eq.0)write(irite,'(a)')
     x        'error - number of atoms too large'
         kill=.true.

      endif
      if(tstep.le.0.d0) then

         if(idnode.eq.0)write(irite,'(a)')
     x        'error - time step assigned zero value'
         kill=.true.

      endif
      iman=(natm+mxnode-1)/mxnode
      klim=((2*kmax+1)**3-1)/2
      nkt=(klim+mxnode-1)/mxnode
      m1=0
      m2=0
      m3=0
      if(lden)m1=6*iman*(kmax+2)+24*klim
      if(lcor)m2=16*nkt*ntime+7*klim+2*ntime
      if(ltim)m3=ntime*(8*nkt+10)+3*klim
      kcore=max(m1,m2,m3)
      if(mcore.lt.kcore)then

         if(idnode.eq.0)write(irite,'(a,i8,a)')
     x        'error - insufficient core allocated. ',kcore,
     x        ' words required'
         kill=.true.

      endif

      if(nbits(ntime).ne.1)then
         if(idnode.eq.0)write(irite,'(a)')
     x        'error - parameter ntime not of the form 2**n'
         kill=.true.

      endif
      if(nbits(ngap).ne.1)then

         if(idnode.eq.0)write(irite,'(a)')
     x        'error - parameter ngap not of the form 2**m'
         kill=.true.

      endif

      if(kill)call exit()

      return
      end
      subroutine fourdc
     x     (idnode,mxnode,natm,ncon,kmax,lchg,xxx,yyy,zzz,
     x     vxx,vyy,vzz,elc,emc,enc,sxyz,vxyz,buff,ekx,eky,ekz)

c***********************************************************************
c
c     calculate spatial fourier transform of current density
c
c     copyright daresbury laboratory 1994
c     author w smith 1994
c
c***********************************************************************
      implicit real*8(a-h,o-z)
      parameter (irite=6,iconf=7,iwork=8)
      parameter (pi=3.141592653589793d0,tpi=6.283185307179586d0)
      character*8 atname
      character*80 header
      logical lchg
      complex*16 ekx,eky,ekz,elc,emc,enc,vxyz,crdc
      complex*16 ckx,cky,ckz,buff
      dimension xxx(*),yyy(*),zzz(*)
      dimension vxx(*),vyy(*),vzz(*)
      dimension ekx(*),eky(*),ekz(*),buff(*)
      dimension elc(*),emc(*),enc(*),vxyz(*),sxyz(*)
      dimension cell(9),rcell(9)
      data cl/1.d0/
      rcl=tpi/cl
      rnatm=sqrt(1.d0/dble(natm))
      klim=((2*kmax+1)**3-1)/2

c     zero mean square density components

      do i=1,4*klim

         sxyz(i)=0.d0

      enddo

c     open history and density files

      open(iconf,file='HISTORY',status='old',err=200)
      open(iwork,file='CUR_SPC',form='unformatted')

      read(iconf,'(a80)')header
      if(idnode.eq.0)then
         write(irite,'(a)')'# HISTORY file header:'
         write(irite,'(a1,a80)')'#',header
      endif
      read(iconf,'(2i10)')levtrj,imcon

c     check file particulars

      if(imcon.eq.0)then
        if(idnode.eq.0)
     x       write(irite,'(a)')
     x       'error - no periodic boundary in HISTORY file'
        call exit()
      endif
      if(levtrj.lt.1)then
        if(idnode.eq.0)
     x       write(irite,'(a)')
     x       'error - no velocity data in HISTORY file'
        call exit()
      endif

      do nstp=1,ncon

        read(iconf,*,end=100)

c     read in configuration data

        read(iconf,*)cell(1),cell(2),cell(3)
        read(iconf,*)cell(4),cell(5),cell(6)
        read(iconf,*)cell(7),cell(8),cell(9)

c     calculate reciprocal lattice basis

        call invert(cell,rcell,det)

c     read atomic data

        do i=1,natm

           if(lchg)then
              read(iconf,*)atname,idx,mass,chg
           else
              read(iconf,*)atname
           endif
          read(iconf,*)xx,yy,zz
          read(iconf,*)vxx(i),vyy(i),vzz(i)
          if(levtrj.eq.2)read(iconf,*)fxx,fyy,fzz
          xxx(i)=xx*rcell(1)+yy*rcell(4)+zz*rcell(7)
          yyy(i)=xx*rcell(2)+yy*rcell(5)+zz*rcell(8)
          zzz(i)=xx*rcell(3)+yy*rcell(6)+zz*rcell(9)
          if(lchg)then
             vxx(i)=chg*vxx(i)
             vyy(i)=chg*vyy(i)
             vzz(i)=chg*vzz(i)
          endif

        enddo

c     calculate fourier exponential terms

        m=0
        do i=idnode+1,natm,mxnode
          m=m+1
          ekx(m)=(1.d0,0.d0)
          eky(m)=(1.d0,0.d0)
          ekz(m)=(1.d0,0.d0)
          elc(m)=cmplx(cos(rcl*xxx(i)),sin(rcl*xxx(i)))
          emc(m)=cmplx(cos(rcl*yyy(i)),sin(rcl*yyy(i)))
          enc(m)=cmplx(cos(rcl*zzz(i)),sin(rcl*zzz(i)))
        enddo
        matm=m
        do l=1,kmax
          do i=1,matm
            ekx(i+l*matm)=ekx(i+matm*(l-1))*elc(i)
            eky(i+l*matm)=eky(i+matm*(l-1))*emc(i)
            ekz(i+l*matm)=ekz(i+matm*(l-1))*enc(i)
          enddo
        enddo

c     start loop over k vectors

        kkk=0
        mmin=0
        lmin=1
        do n=0,kmax
          rn=dble(n)
          do i=1,matm
            enc(i)= ekz(i+n*matm)
          enddo
          do m=mmin,kmax
            rm=dble(m)
            if(m.ge.0)then
              do i=1,matm
                emc(i)= eky(i+m*matm)*enc(i)
              enddo
            else
              do i=1,matm
                emc(i)= conjg(eky(i-m*matm))*enc(i)
              enddo
            endif
            do l=lmin,kmax
              rl=dble(l)
              if(l.ge.0)then
                do i=1,matm
                  elc(i)= ekx(i+l*matm)*emc(i)
                enddo
              else
                do i=1,matm
                  elc(i)= conjg(ekx(i-l*matm))*emc(i)
                enddo
              endif

              rkx=rcell(1)*rl+rcell(2)*rm+rcell(3)*rn
              rky=rcell(4)*rl+rcell(5)*rm+rcell(6)*rn
              rkz=rcell(7)*rl+rcell(8)*rm+rcell(9)*rn
              cnm=rnatm/sqrt(rkx*rkx+rky*rky+rkz*rkz)
              ckx= cnm*conjg(crdc(matm,vxx(idnode+1),mxnode,elc,1))
              cky= cnm*conjg(crdc(matm,vyy(idnode+1),mxnode,elc,1))
              ckz= cnm*conjg(crdc(matm,vzz(idnode+1),mxnode,elc,1))
              vxyz(1+4*kkk)=rkx*ckx+rky*cky+rkz*ckz
              vxyz(2+4*kkk)=rky*ckz-rkz*cky
              vxyz(3+4*kkk)=rkz*ckx-rkx*ckz
              vxyz(4+4*kkk)=rkx*cky-rky*ckx
              kkk=kkk+1

            enddo
            lmin=-kmax
          enddo
          mmin=-kmax
        enddo

c     global sum of fourier transform of current density

        call gdsum(vxyz,8*klim,buff)

c     accumulate mean square values of current density components

        do i=1,4*klim

           sxyz(i)=sxyz(i)+real(vxyz(i)*conjg(vxyz(i)))

        enddo

c     store fourier transform of current density

        if(idnode.eq.0)write(iwork)(vxyz(j),j=1,4*klim)

      enddo

  100 continue

      nstp=nstp-1

      if(idnode.eq.0)then
         write(irite,'(a,i10)')
     x     '# number of timesteps in simulation data    = ',nstp
         close (iwork)
         write(irite,'(a)')
     x     '# CUR_SPC file written and closed'
         write(irite,'(a,a,a)')'#    L    M    N',
     x     '   |Cl(k)|^2','   |Ct(k)|^2'
         k=0
         mmin=0
         lmin=1
         do n=0,kmax
            do m=mmin,kmax
               do l=lmin,kmax
                  slong=sxyz(k+1)/dble(nstp)
                  stran=(sxyz(k+2)+sxyz(k+3)+sxyz(k+4))/dble(nstp)
                  write(irite,'(a1,3i5,1p,2e12.4)')"#",l,m,n,
     x               slong,stran
                  k=k+4
               enddo
               lmin=-kmax
            enddo
            mmin=-kmax
         enddo
         close (iwork)
         write(irite,'(a)')
     x     '# CUR_DEN file written and closed'
      endif

      close (iconf)
      call timchk(1,wtime)

      return

  200 continue

      if(idnode.eq.0)
     x   write(irite,'(a)')'error - HISTORY file not found'
      call exit()

      end
      subroutine correl
     x  (idnode,mxnode,ncon,kmax,ntime,ngap,tstep,ind,num,kvec,cfkt,
     x   ckr0,ckr)

c***********************************************************************
c
c     calculate current density correlation function
c
c     copyright daresbury laboratory 1994
c     author w smith 1994
c
c***********************************************************************

      implicit real*8(a-h,o-z)

      parameter (irite=6,iwork=8,isave1=19)
      complex*16 ckr,ckr0,cfkt
      dimension cfkt(*),ckr0(*),ckr(*),ind(*),num(*),kvec(3,*)

c     open files isave1 and iwork

      open(iwork,file='CUR_SPC',form='unformatted',status='old',err=200)

      if(idnode.eq.0)then

         open(isave1,file='CUR_COR',form='unformatted')

      endif

c     set control parameters

      norg=ntime/ngap
      klim=((2*kmax+1)**3-1)/2
      nblk=mod(klim,mxnode)
      if(idnode.lt.nblk)then
         iblk=4*(klim/mxnode+1)
         ibgn=idnode*iblk
      else
         iblk=4*(klim/mxnode)
         ibgn=idnode*iblk+4*nblk
      endif

      lor=0
      mor=0

c     initialise cfkt array

      do l=1,iblk*ntime
        cfkt(l)=(0.d0,0.d0)
      enddo

      do l=1,ntime
        num(l)=0
      enddo

c     start of loop over time steps

      do n=0,ncon

        read(iwork,end=100)(ckr(j),j=1,4*klim)

        if(mod(n,ngap).le.0)then

          lor=min(lor+1,norg)
          mor=mod(mor,norg)+1
          ind(mor)=1

          do k=1,iblk
            ckr0(mor+ntime*(k-1))=conjg(ckr(k+ibgn))
          enddo

        endif

        do l=1,lor

          m=ind(l)
          ind(l)=m+1
          num(m)=num(m)+1
          do k=1,iblk
            cfkt(m+ntime*(k-1))=cfkt(m+ntime*(k-1))+
     x                          ckr(k+ibgn)*ckr0(l+ntime*(k-1))
          enddo

        enddo

      enddo
  100 continue

      last=min(n-1,ntime)

c     normalise correlation functions, combine transverse components

      i=0
      do k=1,iblk,4

        rnorm1=real(cfkt(1+(k-1)*ntime))
        if(rnorm1.gt.0.d0)rnorm1=dble(num(1))/rnorm1
        rnorm2=real(cfkt(1+k*ntime)+cfkt(1+(k+1)*ntime)+
     x    cfkt(1+(k+2)*ntime))
        if(rnorm2.gt.0.d0)rnorm2=dble(num(1))/rnorm2

        do l=1,last
            cfkt(l+i*ntime)=rnorm1*cfkt(l+(k-1)*ntime)/dble(num(l))
            cfkt(l+(i+1)*ntime)=rnorm2*(cfkt(l+k*ntime)+
     x        cfkt(l+(k+1)*ntime)+cfkt(l+(k+2)*ntime))/dble(num(l))
        enddo
        i=i+2

      enddo

c     calculate k vector indices for printing

      k=0
      mmin=0
      lmin=1
      do n=0,kmax
         do m=mmin,kmax
            do l=lmin,kmax
               k=k+1
               kvec(1,k)=l
               kvec(2,k)=m
               kvec(3,k)=n
            enddo
            lmin=-kmax
         enddo
         mmin=-kmax
      enddo

c     store correlation functions in disc file

      iblk=iblk/2

      if(idnode.eq.0)then

         write(irite,'(a)')
     x   '# Correlation Function of FT of Current Density'
         write(isave1)(cfkt(j),j=1,ntime*iblk)
         call tabulate(idnode,mxnode,kmax,ntime,tstep,kvec,cfkt)

         do kdnode=1,mxnode-1

            call csend(1744549,kdnode,4,kdnode,0)
            call crecv(1744550,cfkt,16*ntime*iblk)
            write(isave1)(cfkt(j),j=1,ntime*iblk)
            call tabulate(kdnode,mxnode,kmax,ntime,tstep,kvec,cfkt)

         enddo

      else

         call crecv(1744549,kdnode,4)
         call csend(1744550,cfkt,16*ntime*iblk,0,0)

      endif

c     close files

      if(idnode.eq.0)then
         close (isave1)
         write(irite,'(a)')
     x      '# CUR_COR file written and closed'
      endif

      close (iwork)
      call timchk(1,wtime)

      return

  200 continue

      if(idnode.eq.0)
     x   write(irite,'(a)')'error - CUR_SPC file not found'
      call exit

      end
      subroutine curfft
     x (idnode,mxnode,kmax,ntime,ngap,tstep,key,kvec,wind,cfkt,work,fta)
c
c***********************************************************************
c
c    calculate fourier transform of current density correlation function
c
c     copyright daresbury laboratory 1994
c     author w smith 1994
c
c***********************************************************************
c
      implicit real*8(a-h,o-z)
      parameter (irite=6,isave1=19,isave2=20)
      parameter (pi=3.141592653589793d0,tpi=6.283185307179586d0)
      complex*16 cfkt(*),work(*),fta(*)
      dimension key(*),wind(*),kvec(3,*)
      data a0,a1,a2/0.42d0,0.50d0,0.08d0/

c     set control parameters

      ntime2=2*ntime
      norg=ntime/ngap
      klim=((2*kmax+1)**3-1)/2
      nblk=mod(klim,mxnode)
      if(idnode.lt.nblk)then
         iblk=2*(klim/mxnode+1)
         ibgn=idnode*iblk
      else
         iblk=2*(klim/mxnode)
         ibgn=idnode*iblk+2*nblk
      endif
      omega=tpi/(dble(ntime2)*tstep)

c     open files isave1 and isave2

      open(isave1,file='CUR_COR',form='unformatted',status='old',
     x    err=200)

      if(idnode.eq.0)then

         open(isave2,file='CUR_FFT',form='unformatted')

      endif

c     initialise complex fast fourier transform routine

      ind=1
      isw=-1
      call fft(ind,isw,ntime2,key,cfkt,work,cfkt)
      ind=0

c     set up window function (blackman function)

      arg=tpi/dble(ntime2)
      do i=1,ntime

         ccc=cos(arg*dble(i+ntime-1))
         wind(i)=a0-a1*ccc+a2*(2.d0*ccc**2-1.d0)

      enddo

c     read the correlation functions from disc

      do i=0,mxnode-1

         if(i.eq.idnode)then

            read(isave1)(cfkt(j),j=1,ntime*iblk)

         else

            read(isave1)

         endif

      enddo

c     loop over correlation functions (longitudinal and transverse)

      do k=1,iblk

c     apply window function

         do j=1,ntime2

            if(j.le.ntime)then

               fta(j)=cfkt(j+ntime*(k-1))*wind(j)

            else

               fta(j)=(0.d0,0.d0)

            endif

         enddo

         fta(1)=fta(1)/2.d0

c     apply complex fourier transform

         call fft(ind,isw,ntime2,key,fta,work,fta)

c     store fourier coefficients

         do j=1,ntime

            cfkt(j+ntime*(k-1))=fta(j)

         enddo

      enddo

c     calculate k vector indices for printing

      k=0
      mmin=0
      lmin=1
      do n=0,kmax
         do m=mmin,kmax
            do l=lmin,kmax
               k=k+1
               kvec(1,k)=l
               kvec(2,k)=m
               kvec(3,k)=n
            enddo
            lmin=-kmax
         enddo
         mmin=-kmax
      enddo

c     save fourier coefficients

      if(idnode.eq.0)then

         write(irite,'(a)')
     x   '# Fourier Transform of Current Density Correlation Functions'
         write(isave2)(cfkt(j),j=1,ntime*iblk)
         call tabulate(idnode,mxnode,kmax,ntime,omega,kvec,cfkt)

         do kdnode=1,mxnode-1

            call csend(1534549,kdnode,4,kdnode,0)
            call crecv(1534550,cfkt,16*ntime*iblk)
            write(isave2)(cfkt(j),j=1,ntime*iblk)
            call tabulate(kdnode,mxnode,kmax,ntime,omega,kvec,cfkt)

         enddo

      else

         call crecv(1534549,kdnode,4)
         call csend(1534550,cfkt,16*ntime*iblk,0,0)

      endif

c     close files

      close (isave1)
      close (isave2)
      if(idnode.eq.0)write(irite,'(a)')
     x   '# CUR_FFT file written and closed'
      call timchk(1,wtime)

      return

  200 continue

      if(idnode.eq.0)
     x   write(irite,'(a)')'error - CUR_COR file not found'
      call exit()

      end
      complex*16 function crdc(n,a,i,b,j)
c
c***********************************************************************
c
c     dl_poly utility routine for dot product of real array A and
c     complex array B
c
c     copyright daresbury laboratory 1994
c     author w smith
c
c***********************************************************************
c
      implicit real*8(a-h,o-z)
      dimension a(n),b(2*n)
      u=sdot(n,a(1),i,b(1),2*j)
      v=sdot(n,a(1),i,b(2),2*j)
      crdc=cmplx(u,v)
      return
      end
      double precision function sdot(n,a,i,b,j)
c
c***********************************************************************
c
c     dl_poly utility routine for dot product of real array A and
c     real array B
c
c     copyright daresbury laboratory 1994
c     author w smith
c
c***********************************************************************
c
      implicit real*8(a-h,o-z)

      dimension a(*),b(*)

      l=1
      m=1
      sdot=0.d0
      do k=1,n

         sdot=sdot+a(l)*b(m)

         l=l+i
         m=m+j

      enddo

      return
      end
      function nbits(n)
c
c**********************************************************************
c
c     dl_poly utility for counting number of set bits in an integer N
c
c     copyright daresbury laboratory 1994
c     author w smith
c
c**********************************************************************
c
      m=n
      nbits=0

      do i=1,128

         if(m.eq.0)return

         nbits=nbits+m-2*(m/2)
         m=m/2

      enddo

      return
      end
      subroutine invert(a,b,d)
c
c***********************************************************************
c
c     dl_poly subroutine to invert a 3 * 3 matrix using cofactors
c
c     copyright - daresbury laboratory 1992
c     author    - w. smith       april 1992
c
c***********************************************************************
c

      real*8 a,b,d,r

      dimension a(9),b(9)
c
c     calculate adjoint matrix
      b(1)=a(5)*a(9)-a(6)*a(8)
      b(2)=a(3)*a(8)-a(2)*a(9)
      b(3)=a(2)*a(6)-a(3)*a(5)
      b(4)=a(6)*a(7)-a(4)*a(9)
      b(5)=a(1)*a(9)-a(3)*a(7)
      b(6)=a(3)*a(4)-a(1)*a(6)
      b(7)=a(4)*a(8)-a(5)*a(7)
      b(8)=a(2)*a(7)-a(1)*a(8)
      b(9)=a(1)*a(5)-a(2)*a(4)
c
c     calculate determinant
      d=a(1)*b(1)+a(4)*b(2)+a(7)*b(3)
      r=0.d0
      if(abs(d).gt.0.d0)r=1.d0/d
c
c     complete inverse matrix
      b(1)=r*b(1)
      b(2)=r*b(2)
      b(3)=r*b(3)
      b(4)=r*b(4)
      b(5)=r*b(5)
      b(6)=r*b(6)
      b(7)=r*b(7)
      b(8)=r*b(8)
      b(9)=r*b(9)
      return
      end
      subroutine fft(ind,isw,ndiv,key,aaa,wfft,bbb)
c***********************************************************************
c
c     fast fourier transform routine
c
c     copyright daresbury laboratory 1994
c     author w smith
c
c***********************************************************************

      implicit real*8(a-h,o-z)

      logical check
      complex*16 aaa(ndiv),bbb(ndiv),wfft(ndiv),ttt
      dimension key(ndiv)
      data tpi/6.2831853072d0/
   10 format(1h0,'error - number of points not a power of two')

c
c     check that array is of suitable length
      nt=1
      check=.true.
      do i=1,20
         nt=2*nt
         if(nt.eq.ndiv)then
            check=.false.
            nu=i
         endif
      enddo
      if(check)then
         write(*,10)
         call exit()
      endif

      if(ind.gt.0)then
c
c     set reverse bit address array

         do kkk=1,ndiv
            iii=0
            jjj=kkk-1
            do j=1,nu
               jj2=jjj/2
               iii=2*(iii-jj2)+jjj
               jjj=jj2
            enddo
            key(kkk)=iii+1
         enddo
c
c     initialise complex exponential factors

         tpn=tpi/dble(ndiv)
         arg=0.d0
         np1=ndiv+1
         np2=ndiv/2
         wfft(1)=(1.d0,0.d0)
         do i=1,np2
            arg=tpn*dble(i)
            wfft(i+1)=cmplx(cos(arg),sin(arg))
            wfft(np1-i)=conjg(wfft(i+1))
         enddo

         return
      endif

c
c     take conjugate of exponentials if required

      if(isw.lt.0)then

         do i=1,ndiv
            wfft(i)=conjg(wfft(i))
         enddo

      endif

c
c     take copy input array

      do i=1,ndiv
         bbb(i)=aaa(i)
      enddo
c
c     perform fourier transform

      kkk=0
      nu1=nu-1
      np2=ndiv/2
      do l=1,nu

  100    do i=1,np2
            iii=key(kkk/2**nu1+1)
            kk1=kkk+1
            k12=kk1+np2
            ttt=bbb(k12)*wfft(iii)
            bbb(k12)=bbb(kk1)-ttt
            bbb(kk1)=bbb(kk1)+ttt
            kkk=kkk+1
         enddo
         kkk=kkk+np2
         if(kkk.lt.ndiv)go to 100
         kkk=0
         nu1=nu1-1
         np2=np2/2

      enddo
c
c     unscramble the fft using bit address array

      do kkk=1,ndiv
         iii=key(kkk)
         if(iii.gt.kkk)then
            ttt=bbb(kkk)
            bbb(kkk)=bbb(iii)
            bbb(iii)=ttt
         endif
      enddo
c
c     restore exponentials to unconjugated values if necessary

      if(isw.lt.0)then

         do i=1,ndiv
            wfft(i)=conjg(wfft(i))
         enddo

      endif

      return
      end
      subroutine tabulate(idnode,mxnode,kmax,ntime,tstep,kvec,cfkt)
c
c*********************************************************************
c
c     routine to write out blocked arrays of complex numbers
c
c     copyright daresbury laboratory
c     author w smith 1994
c
c*********************************************************************
c
      implicit real*8(a-h,o-z)
      parameter (irite=6)
      complex*16 cfkt(*)
      dimension kvec(3,*)

      klim=((2*kmax+1)**3-1)/2
      nblk=mod(klim,mxnode)
      if(idnode.lt.nblk)then
         iblk=2*(klim/mxnode+1)
         ibgn=idnode*iblk
      else
         iblk=2*(klim/mxnode)
         ibgn=idnode*iblk+2*nblk
      endif

      do k=1,iblk

         kk=(ibgn+k-1)/2+1
         if(mod(k,2).eq.1)then
            write(irite,'(a,3i5)')
     x     '# Longitudinal Component ',kvec(1,kk),kvec(2,kk),kvec(3,kk)
         else
            write(irite,'(a,3i5)')
     x     '# Transverse  Component  ',kvec(1,kk),kvec(2,kk),kvec(3,kk)
         endif
         do j=1,ntime

            ttt=tstep*dble(j-1)
            write(irite,'(1p,3e12.4)')ttt,cfkt(j+(k-1)*ntime)

         enddo

         write(irite,'(a)')'&'

      enddo

      return
      end
      subroutine machine(idnode,mxnode)
c
c*********************************************************************
c
c     dl_poly subroutine for obtaining charcteristics of
c     the computer on which the program is being run
c
c     this version is for the intel hypercube
c
c     copyright daresbury laboratory 1992
c     author - w.smith july 1992
c
c*********************************************************************
c

      implicit real*8(a-h,o-z)

c
c     number of nodes available to program

      mxnode=numnodes()

c
c     identity of executing node

      idnode=mynode()

      return
      end
      subroutine timchk(ktim,time)
c
c***********************************************************************
c
c     timing routine (time elapsed in seconds)
c
c     copyright daresbury laboratory 1994
c     author w. smith
c
c***********************************************************************
c
      parameter (irite=6)
      real*8 time
      save init
      data init/0/
      idnode=mynode()
c
c     synchronise time
      call gsync()
      if(init.eq.0)init=mclock()
      itime=mclock()-init
      time=dble(itime)/1000.d0
      if(ktim.gt.0.and.idnode.eq.0)write(irite,
     x  "('# time elapsed since job start = ',f15.8,' seconds',/)")
     x  time
      return
      end


text:
(*  This version Willem Ouwerkerk - April 2021 - bitbang SPI driver
    With use of some of the data from Mecrisp, but much more
    the GD32VF103 user manual V1.2
    Please note that 27 MHz is the maximum SPI clock frequency!

        40010800 =  PORTA
        40010C00 =  PORTB
        40011000 =  PORTC
        40011400 =  PORTD
        40011800 =  PORTE

        40021000 =  RCU_CTL
        40021018 =  RCU_APB2EN

        40013000 =  SPI0_CTL0
        40013008 =  SPI0_STAT
        4001300C =  SPI0_DATA

\ These bits are connected to flash memory on SEEED dev. board
PA5     = SCK     - 1       40010800
PA6     = MISO    - 8
PA7     = MOSI    - 1
PC0     = CE      - 1       40011000

*)

hex
v: fresh
: SPI-SETUP         ( +n -- )
    drop                    \ Remove compatability dummy
\ AF output mode push-pull 10MHz
    FFF00000 40010800 **bic \ Port_A CRL  Clear pins PA5,6,7
    18100000 40010800 **bis \ Port_A CRL  Set pins PA5,6,7
    40 4001080C **bis       \ Port_A OCTL PA.6 Activate pull-up
\ Output mode push-pull 10MHz
    000000FF 40011000 **bic \ Port_C CRL  Clear pin PC0 & PC1
    00000011 40011000 **bis \ Port_C CRL  Set pin PC0 & PC1
    1 4001100C **bis        \ Port_C out  CS high
    20 4001080C **bic ;     \ Port_A OCTL PA.5  CLK = 0

false [if]  \ Clock high level version = 300 kHz, low level version = 1 MHz

: CLOCK-HI     ( -- )      20 4001080C **bis ;      \ PA5_OCTL  SPI clock
: CLOCK-LOW    ( -- )      20 4001080C **bic ;      \ PA5_OCTL  SPI clock
: CLOCK        ( -- )      clock-hi  clock-low ;
: {SPI         ( -- )      1 4001100C **bic ;       \ PC0_OCTL  SPI on, CS=low
: SPI}         ( -- )      1 4001100C **bis ;       \ PCO_OCTL  SPI off, CS=high

\ Write a bit to the SPI-bus, read abit from the SPI-bus
: WRITE-BIT    ( b -- )
    80 and if 80 4001080C **bis else 80 4001080C **bic then ; \ PA7_OCTL

: READ-BIT     ( -- 0|1 )  40 40010808 bit** 0<> 1 and ; \ PA0_STAT

[else]

code CLOCK-HI     ( -- )
    sun 4001080C li     \ PORTA_ODR  Port-A output address
    day 20 li           \ Bit-5 mask
    w sun ) .mov        \ Read Port-A
    w day .or           \ Set bit-0
    sun ) w .mov        \ Write Port-A
    next
end-code

code CLOCK-LOW    ( -- )
    sun 4001080C li     \ PORTA_ODR  Port-A output address
    day -21 li          \ Bit-5 inverted mask
    w sun ) .mov        \ Read Port-A
    w day .and          \ Clear bit-5
    sun ) w .mov        \ Write Port-A
    next
end-code

code CLOCK      ( -- )
    sun 4001080C li     \ PORTA_ODR  Port-A output address
    day 20 li           \ Bit-5 mask
    w sun ) .mov        \ Read Port-A
    w day .or           \ Set bit-0
    sun ) w .mov        \ Write Port-A
    sun ) w .mov        \ Delay a few cycles
    day -21 li          \ Bit-5 inverted mask
    w sun ) .mov        \ Read Port-A
    w day .and          \ Clear bit-5
    sun ) w .mov        \ Write Port-A
    next
end-code

code  {SPI      ( -- )
    sun 4001100C li     \ PORTC_ODR  Port-C output address
    w sun ) .mov        \ Read Port-C
    w -2 .andi          \ Clear bit-0
    sun ) w .mov        \ Write Port-C
    next
end-code

code SPI}        ( -- )
    sun 4001100C li     \ PORTC_ODR  Port-C output address
    day 1 .li           \ Bit-0 mask
    w sun ) .mov        \ Read Port-C
    w day .or           \ Set bit-0
    sun ) w .mov        \ Write Port-C
    next
end-code

code WRITE-BIT  ( b -- )
    sun 40010800 li     \ Port_A base address
    moon C sun x) .mov  \ Read Port_A OCTL register
    day -81 li          \ Prepare inverted bit mask for bit-7
    moon day .and       \ Remove bit-7 from  MOON
    day 80 li           \ Prepare bit mask for bit-7
    tos day .and        \ Leave only bit-7 from TOS
    moon tos .or        \ Add bit-7 to MOON
    C sun x) moon .mov  \ Write MOON back
    tos sp )+ .mov      \ Pop stack
    next
end-code

code READ-BIT   ( -- 0|1 )
    sp -) tos .mov
    moon 40 li          \ Bit-6 mask to MOON
    sun 40010800 li     \ Port_A base address
    tos 8 sun x) .mov   \ Read Port_A STAT register
    tos moon .and       \ Leave bit-6 only
    tos 6 .srli         \ Bit-6 to bit-0
    next
end-code

[then]

: SPI-I/O       ( b0 -- b1 )
    8 for
        dup write-bit  2*  clock-hi
        read-bit or  clock-lo
    next  FF and ;

: SPI-OUT       ( b -- )
    8 for  dup write-bit  2*  clock  next  drop ;

: SPI-IN        ( -- b )
    0  8 for  2*  clock-hi  read-bit or  clock-lo  next ;

: SPI-ON        ( -- )      1 spi-setup ;

spi-on
v: fresh
shield spi\  freeze

\ End ;;;

max_stars_repo_name: embeddingforth/embeddinForth

max_issues_repo_name: embeddingforth/embeddinForth

max_forks_repo_name: embeddingforth/embeddinForth


text:
        subroutine iddp_aid0(eps,m,n,a,krank,list,proj,rnorms)
c
c       uses routine iddp_id to ID a without modifying its entries
c       (in contrast to the usual behavior of iddp_id).
c
c       input:
c       eps -- precision of the decomposition to be constructed
c       m -- first dimension of a
c       n -- second dimension of a
c
c       output:
c       krank -- numerical rank of the ID
c       list -- indices of the columns in the ID
c       proj -- matrix of coefficients needed to interpolate
c               from the selected columns to the other columns in a;
c               proj doubles as a work array in the present routine, so
c               must be at least m*n real*8 elements long
c
c       work:
c       rnorms -- must be at least n real*8 elements long
c
c       _N.B._: proj must be at least m*n real*8 elements long
c
        implicit none
        integer m,n,krank,list(n),j,k
        real*8 eps,a(m,n),proj(m,n),rnorms(n)
c
c
c       Copy a into proj.
c
        do k = 1,n
          do j = 1,m
            proj(j,k) = a(j,k)
          enddo ! j
        enddo ! k
c
c
c       ID proj.
c
        call iddp_id(eps,m,n,proj,krank,list,rnorms)
c
c
        return
        end
c
c
c
c

max_issues_repo_name: birichie/catboost

max_forks_repo_name: birichie/catboost


text:
      PROGRAM FGKE04
C
C Define the error file, the Fortran unit number, the workstation type,
C and the workstation ID to be used in calls to GKS routines.
C
C     PARAMETER (IERRF=6, LUNIT=2, IWTYPE=1,  IWKID=1)   ! NCGM
C     PARAMETER (IERRF=6, LUNIT=2, IWTYPE=8,  IWKID=1)   ! X Windows
C     PARAMETER (IERRF=6, LUNIT=2, IWTYPE=11, IWKID=1)   ! PDF
C     PARAMETER (IERRF=6, LUNIT=2, IWTYPE=20, IWKID=1)   ! PostScript
C
      PARAMETER (IERRF=6, LUNIT=2, IWTYPE=1,  IWKID=1)

      PARAMETER (IWKID1=5, IWKID2=7)
C
C  Illustrate workstation control functions by intermixing plotting to
C  a metafile with plotting to X window workstations.
C
C
C  Open GKS, open and activate a workstation.
C
      CALL GOPKS (IERRF, ISZDM)
      CALL GOPWK (IWKID, LUNIT, IWTYPE)
C
C  Define color index 6 as red on that workstation.
C
      CALL GSCR(IWKID,6,1.,0.,0.)
C
C  Set character height and set the text alignment to (center, half).
C
      CALL GSCHH(.04)
      CALL GSTXAL(2,3)
C
C  Activate the CGM workstation.
C
      CALL GACWK(IWKID)
C
C  Create picture 1 in the metafile and call FRAME to terminate
C  the first picture.
C
      CALL GSTXCI(6)
      CALL GTX(.5,.5,'Picture 1')
      CALL FRAME
C
C  Put another picture in the metafile.
C
      CALL GTX(.5,.5,'Picture 2')
      CALL FRAME
C
C  Open and activate an X11 workstation and define color index 6
C  to be cyan on that workstation (this workstation will be referred
C  to as "The first window".
C
      CALL GOPWK(IWKID1,0,8)
      CALL GACWK(IWKID1)
      CALL GSCR(IWKID1,6,0.,1.,1.)
C
C  Create picture 3.  This will be plotted to the CGM workstation as
C  well as to the X11 window since they are both active.
C
      CALL GTX(.5,.5,'Picture 3')
C
C  Terminate the metafile picture.
C
      CALL NGPICT(IWKID,1)
C
C  Pause in the X window with a "<READY>" prompt and wait for a
C  mouse click.  The window will be cleared after the mouse click.
C
      CALL NGPICT(IWKID1,4)
C
C  Open and activate another X workstation (to be reffered to as "The
C  second window") and define color index 6 to be green on that
C  workstation.
C
      CALL GOPWK(IWKID2,0,8)
      CALL GACWK(IWKID2)
      CALL GSCR(IWKID2,6,0.,1.,0.)
C
C  Plot picture 4.  This will be sent to the CGM workstation and the
C  two X window workstations since they are all active.
C
      CALL GTX(.5,.5,'Picture 4')
C
C  Terminate picture 4 in the metafile.
C
      CALL NGPICT(IWKID,1)
C
C  Make the second window current.
C
      CALL NGPICT(IWKID2,0)
C
C  Pause in the first window with a "<READY>" prompt and wait for a
C  mouse click.  The window will be cleared after the click.
C
      CALL NGPICT(IWKID1,4)
C
C  Clear the second window.
C
      CALL NGPICT(IWKID2,1)
C
C  Deactivate the metafile (but do not close it) and draw picture 5.
C  This will go to the two active X11 workstations, but not the CGM.
C
      CALL GDAWK(IWKID)
      CALL GTX(.5,.5,'Picture 5')
C
C  Make the second window current
C
      CALL NGPICT(IWKID2,0)
C
C  Pause in the first window waiting or a mouse click
C
      CALL NGPICT(IWKID1,4)
C
C  Re-activate the metafile and deactivate the second window.
C
      CALL GACWK(IWKID)
      CALL GDAWK(IWKID2)
C
C  Plot picture 6.  This will go to the first window and to the
C  metafile.
C
      CALL GTX(.5,.5,'Picture 6')
C
C  Terminate the picture in the CGM.
C
      CALL NGPICT(IWKID,1)
C
C  Pause in the first window waiting for a mouse click.
C
      CALL NGPICT(IWKID1,4)
C
C  Reactivate the second window and clear it.
C
      CALL GACWK(IWKID2)
      CALL NGPICT(IWKID2,1)
C
C  Put out picture 7.  This will go to the all active workstations.
C
      CALL GTX(.5,.5,'Picture 7')
C
C  Terminate the picture in the CGM.
C
      CALL NGPICT(IWKID,1)
C
C  Make the first window current.
C
      CALL NGPICT(IWKID1,0)
C
C  Pause in the second window waiting for a mouse click.
C
      CALL NGPICT(IWKID2,4)
C
C  Deactivate and close the first window; deactivate and close the CGM.
C
      CALL GDAWK(IWKID1)
      CALL GDAWK(IWKID)
      CALL GCLWK(IWKID1)
      CALL GCLWK(IWKID)
C
C  Put out picture 8.  This will go to the second window, the only
C  active workstation.
C
      CALL GTX(.5,.5,'Picture 8')
C
C  Pause in the second window with a "<READY>" prompt.
C
      CALL NGPICT(IWKID2,4)
C
C  Deactivate and close the second window.
C
      CALL GDAWK(IWKID2)
      CALL GCLWK(IWKID2)
C
C  Close GKS.
C
      CALL GCLKS
C
      STOP
      END


text:
	PROGRAM SINCAL
C
C	FITS A SINE WAVE TO CALIBRATE TDS A/D CONVERTER
C
	include		'wind_examples:wind_tm_routine_def.for'
	include		'wind_examples:wind_return_code_def.for'
	integer*4	ok
	parameter	event='TDSF'
	integer*4	i,j,k,n,itemp
	integer*4	get_stream_name
	integer*4	wait_for_events			! an entry point
	integer*4	err_count
	integer*4	ichpa(6)
	character*80	stream
	character*4	pa(6,4)
	parameter	size=2048
	integer*4	return_size
	integer*4	tds_channel,ifilf,ifils,ifsps,issps
	integer*4	temp_waves
	character*32	s
	character*32	s_scet
	integer*4	scet(2)
	integer*4	major, minor
	character*80	file
	character*80	title
	character*32	item
	integer*4	ios
	integer*4	NREM,NHRMN,IHRMN
	real ffilter,sfilter,fsps,ssps
!
	common /nrblk/ nrem,NHRMN
	common /headblk/ major,minor,tds_channel,s_scet

C
	CHARACTER*120 JUNK(10)
	CHARACTER*18 PTITLE(10)
	INTEGER*4 TDSCH,hkch
	COMMON /CALBLK/ VOLTS(2048),NDATA(2060)
	COMMON /HEADBL/ JUNK
	COMMON /FITBLK/ NPT,VOLTIN(2048),TM(2048)
	COMMON /ZBLK/ NZER,ZCROSS(128)
	COMMON /PLOTDT/ ITERM,PTITLE
	COMMON /HNTBLK/ NHUNT(25)
	EXTERNAL FITDATA,FITZERO
	DIMENSION IZCROSS(4,128),IZCNTS(128)
	DIMENSION X(25),DX(25),Y(25)
	DIMENSION FFILTER(4),SFILTER(4),FSPS(4),SSPS(4)
C	DATA ITERM /3/
	DATA ITERM /-6/
	DATA TWOPI /6.2831853/
	DATA FFILTER /50000.,12500.,3125.,781./
	DATA SFILTER /3125.,781.,195.,49./
	DATA FSPS /120.,30.,7.5,1.8/
	DATA SSPS /7500.,1875.,469.,117./
	DATA PA /'EXAC','EXAC','EXDC',' BX ',' BY ',' BZ ',
     1           'EXDC','EYAC','EYDC','EXDC','EYDC','EZDC',
     2           '    ','EZAC','EZDC','    ','    ','    ',
     3           '    ','EZAC','EZDC','    ','    ','    '/
	DATA LOLIM /20/
C	INSERT, FREQ,AMPL,SPS,P/A,MODEL,CHANNEL,FILTER
C	DATA FREQ /1000./
	DATA FREQ /95./
C	DATA AMPL /4.0/				! AMPLITUDE IN V p-p
C	DATA AMPL /.4/				! AMPLITUDE IN V p-p
	DATA AMPL /2./				! AMPLITUDE IN V p-p
C	DATA SPS /1875./			! SAMPLE RATE, PER SEC.
	DATA SPS /7500./			! SAMPLE RATE, PER SEC.
C
	PTITLE(3) = 'P/A EXDC'
C	PTITLE(1) = 'ENG. MODEL'
	PTITLE(1) = 'FLIGHT MODEL'
C	PTITLE(2) = 'TDS CH 3'
	WRITE(PTITLE(5),1005) FREQ
 1005	FORMAT(F5.0,' HZ')
	WRITE(PTITLE(6),1006) AMPL
 1006	FORMAT(F5.2,' V P-P')
	PRINT*,'LOLIM=',LOLIM
	PRINT*,' '
C
 1001	FORMAT(I4,1X,20I5)
 1002	FORMAT(A)
C
C	GET STARTED
C
!
	ok = get_stream_name(stream)
	if (.not. ok) stop 'no file supplied.'

	ok = wind_tm_open_channel(tdsch,stream)
	if (.not. ok) stop 'Cannot open tds channel'
	if( ok ) print*,'opened tds channel',tdsch

c	ok = wind_tm_open_channel(hkch,stream)           !pjk
c	if (.not. ok) stop 'Cannot open hk channel'      !pjk

	ok = wind_tm_get_filename(tdsch,file)
	type*, 'file',file

	ok = wind_tm_get_mfmf(tdsch,major,minor)
	if (.not. ok) stop 'cannot get stream position'
	if( ok ) print*,'frame',major,minor
c
	get_tm_stream = 1
c
C           ok = wind_tm_get_next_event(tdsch,major,minor,'HK')
C	   item = 'TEMP_WAVES2'
C	   ok = wind_tm_get_item(tdsch, item, temp_waves, 1, return_size)
C
C	IF OFFLINE, FIND DESIRED TIME
C
        if(stream.eq.'offline') then
 120           ok = wind_tm_get_next_event(tdsch,major,minor,'TDSF')
	   item = 'EVENT_SCET'
	   TYPE*,'AT EVENT_SCET'
	   ok = wind_tm_get_item(tdsch, item, scet, 2, return_size)
	   item = 'EVENT_NUMBER'
	   ok = wind_tm_get_item(tdsch, item, itemp, 1, return_size)
	   if (.not. ok) type *, 'cannot get item ', item, ', ok=', ok
	TYPE*,'SCET',SCET,' event no', itemp
	   ihrmn = scet(2)/100
	   if(ihrmn.lt.nhrmn) then
	     ok = wind_tm_increment_mfmf(major,minor)
	     if ( .not. ok) type *,' cannot find time'
c	     type*,' s/c time' ,ihrmn
c	     if( ok ) go to 110
	     go to 120
	   endif
c		write(26,*) 'data collection at ',scet
	endif
c
	! this is the main program loop

	type*,'going to get next tds event'
C
C	GET NEXT EVENT
C
c 110    if( wind_tm_eof(tdsch,major,minor)) stop
 110    continue

	ok = wind_tm_get_next_event(tdsch,major,minor,event)
	   if (.not. ok) then
	      type *, char(7), '******** missing packet in event ********'
	      type *, 'Cannot get event at MF.mf: ', major, minor
C	      if (wind_tm_realtime(ch)) then
C	         ok = wind_tm_get_latest_mfmf(ch,major,minor)
C	         type *, '-reset to latest position: ', major, minor
C	      else
C	         call wind_tm_increment_packet(major,minor)
C	         type *, '-incrementing packet...'
C	      end if
	   end if

	   item = 'CHANNEL'
	   ok = wind_tm_get_item(tdsch, item, tds_channel, 1, return_size)
	TYPE*,'CHANNEL',tds_channel
	   item = 'EVENT_NUMBER'
	   ok = wind_tm_get_item(tdsch, item, itemp, 1, return_size)
	   type*,'event number',itemp
	   item = 'SLOW_RX_SPEED'
	   ok = wind_tm_get_item(tdsch, item, issps, 1, return_size)
	   type*,'slow rx speed',issps
	   item = 'FAST_RX_SPEED'
	   ok = wind_tm_get_item(tdsch, item, ifsps, 1, return_size)
	   type*,'fast rx speed',ifsps
	   item = 'SLOW_RX_FILTER'
	   ok = wind_tm_get_item(tdsch, item, ifils, 1, return_size)
	   type*,'slow rx filter',ifils
	   item = 'FAST_RX_FILTER'
	   ok = wind_tm_get_item(tdsch, item, ifilf, 1, return_size)
	   type*,'fast rx filter',ifilf
	   item = 'SOURCE_CHAN_1'
	   ok = wind_tm_get_item(tdsch, item, ichpa(1), 1, return_size)
	   item = 'SOURCE_CHAN_2'
	   ok = wind_tm_get_item(tdsch, item, ichpa(2), 1, return_size)
	   item = 'SOURCE_CHAN_3'
	   ok = wind_tm_get_item(tdsch, item, ichpa(3), 1, return_size)
	   item = 'SOURCE_CHAN_4'
	   ok = wind_tm_get_item(tdsch, item, ichpa(4), 1, return_size)
	   item = 'SOURCE_CHAN_5'
	   ok = wind_tm_get_item(tdsch, item, ichpa(5), 1, return_size)
	   item = 'SOURCE_CHAN_6'
	   ok = wind_tm_get_item(tdsch, item, ichpa(6), 1, return_size)
C
	   WRITE(PTITLE(2),1009)  TDS_CHANNEL
	   ipa = ichpa(tds_channel)
	   WRITE(PTITLE(3),1007) pa(tds_channel,ipa+1)
 1007	   FORMAT('P/A ',A4)
	   IF(TDS_CHANNEL.LE.2) THEN
	      WRITE(PTITLE(4),1004) FSPS(IFSPS+1)
	      SPS = 1000.*FSPS(IFSPS+1)
	      WRITE(PTITLE(7),1008) FFILTER(IFILF+1)
	   ELSE
	      WRITE(PTITLE(4),1004) SSPS(ISSPS+1)
	      SPS = SSPS(ISSPS+1)
	      WRITE(PTITLE(7),1008) SFILTER(IFILS+1)
	   ENDIF
 1004	   FORMAT('SPS',F8.1)
 1008	   FORMAT('FILTER',F7.0,' HZ')
 1009	   FORMAT('TDS CHANNEL',I4)

C	   if (ok) type *, 'waves2 temperature', temp_waves2
CW		write(16,*) 'temp in t/m counts',temp_waves
C		write(26,*) 'temp in t/m counts',temp_waves
c	   ok = wind_tm_xlate_item(ch, 'HK', item, temp_waves, title)
c		type*,'temp in degrees?',title
c	   if (.not. ok) type *, 'cannot get item ', item, ', ok=', ok
c
C
	ITEM = 'DATA'
	OK = WIND_TM_GET_ITEM(TDSCH,ITEM,NDATA,SIZE,RETURN_SIZE)
	IF(.NOT.OK) TYPE*,'CANNOT GET ITEM = DATA'
c	  PRINT*,'DATA'
c	  PRINT 222, (NDATA(J),J=1,2048)
 222	FORMAT(10I6)
	DO IK = 1,2048
	  NDATA(IK) = NDATA(IK)-128
	ENDDO
C
C	CALCULATE NOMINAL WAVE PERIOD IN SAMPLES
C
	WPER = SPS/FREQ
	PRINT*,'WPER',WPER
C
C	FIND FIRST ZERO CROSSING
C
	DO IL = 1,2047
	  IZ = IL
	  IF(NDATA(IL).EQ.0) PRINT*,'ZERO DATA',IL,NDATA(IL-1),
     1   NDATA(IL),NDATA(IL+1)
	  IF(NDATA(IL)*NDATA(IL+1).LE.0) GO TO 20
	ENDDO
 20	CONTINUE
C	PRINT*,'FIRST ZERO CROSSING AT IZ=',IZ
C
C
C	  CALCULATE ZERO CROSSINGS IN GROUPS FOR UP TO 4 GROUPS
C
	  IZCROSS(1,1) = IZ
	  IZCNT = 1
	  IL1 = IZ+1
	  IL2 = IL1 + WPER/8.
	  IL2 = MIN0(IL2,2047)
	  DO IGP = 1,4
 	    DO IL = IL1,IL2
	      IF(NDATA(IL)*NDATA(IL+1).LE.0) THEN
	        IZCNT = MIN0(IZCNT+1,128)
	        IZCROSS(IGP,IZCNT) = IL
	      ENDIF
	    ENDDO
	    IL1 = IZCROSS(IGP,1) + 3.*WPER/8.
	    IL2 = IL1 + WPER/4.
	    IL2 = MIN0(IL2,2047)
	    IZCNTS(IGP) = IZCNT
	    IZCNT = 0
  	  ENDDO
C
C	  CALCULATE AVERAGES FOR MULTIPLE ZERO CROSSINGS
C
	  DO IGP = 1,4
	    IZCNT = IZCNTS(IGP)
	    PRINT*,'ZEROS',IGP,(IZCROSS(IGP,I),I=1,IZCNT)
	    ZCROSS(IGP) = 0.
	    DO I = 1,IZCNT
	      ZCROSS(IGP) = ZCROSS(IGP) + IZCROSS(IGP,I)
	    ENDDO
	    IF(IZCNT.NE.0) ZCROSS(IGP) = ZCROSS(IGP)/IZCNT + .5
	  ENDDO
	  PRINT*,'Z CROSSINGS',(ZCROSS(J),J=1,4)
	  WPER1 = .5*(ZCROSS(3) - ZCROSS(1) + ZCROSS(4) - ZCROSS(2))
	  PRINT*,'NEW WPER',WPER1
	  PHI0 = .25*(ZCROSS(1)+ZCROSS(2)+ZCROSS(3)+ZCROSS(4)-2.*WPER1)
	  PRINT*,'PHI0',PHI0
	  NPHI0 = PHI0 + .5
C
C	IF FREQUENCY IS HIGH, THERE WILL BE A LOT OF ZERO CROSSINGS AND
C	ALL WILL BE SINGLE, SO FIND SOME MORE
C
	IF(WPER.LT.820.) THEN
	    PRINT*,'FIND SOME MORE ZEROS'
	  DO JT = 1,127
	    IL1 = ZCROSS(IZCNT) + 3.*WPER/8.
	    IL2 = IL1 + WPER/4.
	    IL2 = MIN0(IL2,2047)
	    IF(IL1.GT.2047) GO TO 25
C	    PRINT*,'IL1,IL2',IL1,IL2
 	      DO IL = IL1,IL2
	        IF(NDATA(IL)*NDATA(IL+1).LE.0) THEN
		  IZCNT = IZCNT+1
	          ZCROSS(IZCNT) = IL+.5
	          IL1 = IL
		  GO TO 22
	        ENDIF
	      ENDDO
 22	      CONTINUE
  	  ENDDO
C
 25	  CONTINUE
C
	  PRINT*,IZCNT,' ZEROS FOUND'
	  PRINT*,(ZCROSS(J),J=1,IZCNT)
	  PRINT*,'INTERVALS'
	  PRINT*,((ZCROSS(J+1)-ZCROSS(J)),J=1,IZCNT-1)
	  WPER1 = 2.*(ZCROSS(IZCNT)-ZCROSS(1))/(IZCNT-1)
	  PRINT*,'NEW WPER',WPER1
	  X(1) = ZCROSS(1)
	  X(2) = 2./WPER1
	  DX(1) = .5
	  DX(2) = .05*X(2)
	  CALL HUNTMN(2,X,DX,Y,FITZERO,SUMSQ)
	  WPER1 = 2./X(2)
	  PHI0 = X(1) + .25*WPER1
	  PRINT*,'NEW WPER,PHI0=',WPER1,PHI0
	  NPHI0 = PHI0 + .5
C
C	  END OF FURTHER ZERO CROSSINGS FOR HIGH FREQUENCIES
C
	ENDIF
C
C
C	ALLOW FOR POSSIBILITY OF NEGATIVE PEAK
C
 	AMPT = AMPL
	IF(NDATA(NPHI0).LT.0) AMPT = -AMPL
	DO J = 1,2048
	  PHI = (J-PHI0)*TWOPI/WPER1
	  VOLTS(J) = .5*AMPT*COS(PHI)
	ENDDO
c	WRITE(27,1003) (J,VOLTS(J), NDATA(J),J=1,2048)
 1003	FORMAT(3(I9,E11.3,I5))
C
C	NOW A SINE WAVE HAS BEEN FITTED THROUGH THE ZERO CROSSINGS AND
C	THE KNOWN AMPLITUDE, SO VOLTAGE INPUT FOR EACH DATA POINT HAS
C	BEEN CALCULATED
C
C	FIT UPPER PART OF CALIBRATION CURVE, WITH ZERO OFFSET
C
	NPT = 0
	DO J = 1,2048
	  IF(IABS(NDATA(J)).GT.71.AND.NDATA(J).GT.LOLIM) THEN
	    NPT = NPT+1
	    VOLTIN(NPT) = ABS(VOLTS(J))
	    TM(NPT) = NDATA(J)
	  ENDIF
	ENDDO
C
	NHUNT(2) = 0
C	X(1) =20./.7238
	X(1) = 31.
	X(2) = 0.
	X(3) = 2.8E-3
	X(3) = X(1)*ALOG10(X(3))
C
	DX(1) = .05*X(1)
	DX(2) = 5.E-3
	DX(3) = .1*X(3)
	IF(NPT.LE.3) GO TO 30
	CALL FITDATA(X,SUMSQS)
	PRINT*,' '
	PRINT*,' '
	PRINT*,' FIRST TRY AT HIGH FIT, +DATA,',NPT,'  POINTS, SUMSQ='
     1  ,SUMSQS
	DO IT = 1,20
	  CALL HUNTMN(3,X,DX,Y,FITDATA,SUMSQ)
	  PRINT*,'HUNTMN OUTPUT',X(1),X(2),X(3),SUMSQ
	  IF(SUMSQ.GT..9995*SUMSQS) GO TO 30
	  SUMSQS = SUMSQ
	ENDDO
	PRINT*,'NOT ENOUGH ITERATIONS, POOR HIGH FIT'
 30	CONTINUE
	IF(X(1).NE.0.) VZERO = 10.**(X(3)/X(1))
	IF(NPT.GT.3) THEN
	  RMS = SQRT(SUMSQ/(NPT-2))
	  PRINT*,' VZERO = OLD X(3) =',VZERO,'  RMS ERR',RMS,'  TM UNITS'
	  DO IT = 1,5
	    IF(DX(3).EQ.0.) DX(3) = 1.E-5
	    CALL FINMIN(3,X,DX,Y,FITDATA,SUMSQF)
	    CALL FINMIN(3,Y,DX,X,FITDATA,SUMSQF)
	    PRINT*,' FINMIN GIVES',(X(N),N=1,3),SUMSQF
	  ENDDO
	  V71 = 10.**((71.+X(3))/X(1)) - X(2)
	  PRINT*,'VOLTS TO GIVE TM = 71', V71
	ELSE
	  PRINT*,'NOT ENOUGH DATA FOR HIGH, +DATA FIT'
	ENDIF
C
	NPT = 0
	DO J = 1,2048
	  IF(IABS(NDATA(J)).GT.71.AND.NDATA(J).LT.-LOLIM) THEN
	    NPT = NPT+1
	    VOLTIN(NPT) = ABS(VOLTS(J))
	    TM(NPT) = -NDATA(J)
	  ENDIF
	ENDDO
C
	NHUNT(2) = 0
	X(2) = 0.
C
	DX(1) = .05*X(1)
	DX(2) = 5.E-3
	DX(3) = .1*X(3)
	IF(NPT.LE.3) GO TO 35
	CALL FITDATA(X,SUMSQS)
	PRINT*,' '
	PRINT*,' '
	PRINT*,'/FIRST TRY AT HIGH FIT, -DATA,',NPT,'  POINTS, SUMSQ='
     1  ,SUMSQS
	DO IT = 1,20
	  CALL HUNTMN(3,X,DX,Y,FITDATA,SUMSQ)
	  PRINT*,'HUNTMN OUTPUT',X(1),X(2),X(3),SUMSQ
	  IF(SUMSQ.GT..9995*SUMSQS) GO TO 35
	  SUMSQS = SUMSQ
	ENDDO
	PRINT*,'NOT ENOUGH ITERATIONS, POOR HIGH FIT'
 35	CONTINUE
	VZERO = 10.**(X(3)/X(1))
	IF(NPT.GT.3) THEN
	  RMS = SQRT(SUMSQ/(NPT-2))
	  PRINT*,' VZERO = OLD X(3) =',VZERO,'  RMS ERR',RMS,'  TM UNITS'
	  DO IT = 1,5
	    IF(DX(3).EQ.0.) DX(3) = 1.E-5
	    CALL FINMIN(3,X,DX,Y,FITDATA,SUMSQF)
	    CALL FINMIN(3,Y,DX,X,FITDATA,SUMSQF)
	    PRINT*,' FINMIN GIVES',(X(N),N=1,3),SUMSQF
	  ENDDO
	  V71 = 10.**((71.+X(3))/X(1)) - X(2)
	  PRINT*,'VOLTS TO GIVE TM = 71', V71
	ELSE
	  PRINT*,'NOT ENOUGH DATA FOR HIGH, -DATA FIT'
	ENDIF
C
C
C	FIT LOWER PART OF CALIBRATION CURVE, WITH NONZERO OFFSET
C
	NPT = 0
	DO J = 1,2048
	  IF(IABS(NDATA(J)).LE.71.AND.NDATA(J).GT.LOLIM) THEN
	    NPT = NPT+1
	    VOLTIN(NPT) = ABS(VOLTS(J))
	    TM(NPT) = NDATA(J)
	  ENDIF
	ENDDO
	NHUNT(2) = 1
	X(2) = 2.3E-2
	DX(1) = .05*X(1)
	DX(2) = .1*X(2)
	DX(3) = 1.E-5
	CALL FITDATA(X,SUMSQS)
	PRINT*,' '
	PRINT*,' '
	PRINT*,' FIRST TRY AT FULL FIT, +DATA,',NPT,'  POINTS, SUMSQ='
     1  ,SUMSQS
	DO IT = 1,20
	  CALL HUNTMN(3,X,DX,Y,FITDATA,SUMSQ)
	  PRINT*,'HUNTMN OUTPUT',X(1),X(2),X(3),SUMSQ
	  IF(SUMSQ.GT..9995*SUMSQS) GO TO 40
	  SUMSQS = SUMSQ
	ENDDO
	PRINT*,'NOT ENOUGH ITERATIONS, POOR FULL FIT'
 40	CONTINUE
	VZERO = 10.**(X(3)/X(1))
	IF(NPT.GE.4) RMS = SQRT(SUMSQ/(NPT-3))
	PRINT*,' VZERO = OLD X(3) =',VZERO,'  RMS ERR',RMS,'  TM UNITS'
	PRINT*,'HUNTMN OUTPUT',X(1),X(2),X(3),SUMSQ
	SUMSQS = SUMSQ
	  DO IT = 1,10
	    CALL FINMIN(3,X,DX,Y,FITDATA,SUMSQF)
	    CALL FINMIN(3,Y,DX,X,FITDATA,SUMSQF)
	    PRINT*,' FINMIN GIVES',(X(N),N=1,3),SUMSQF
	    IF(SUMSQF.GT..9999*SUMSQS) GO TO 47
	    SUMSQS = SUMSQF
	  ENDDO
  47	  V71 = 10.**((71.+X(3))/X(1)) - X(2)
	  PRINT*,' '
	  PRINT*,' FINMIN GIVES',(X(N),N=1,3),SUMSQF
	  PRINT*,'VOLTS TO GIVE TM = 71', V71
C
	NPT = 0
	DO J = 1,2048
	  IF(IABS(NDATA(J)).LE.71.AND.NDATA(J).LT.-LOLIM) THEN
	    NPT = NPT+1
	    VOLTIN(NPT) = ABS(VOLTS(J))
	    TM(NPT) = -NDATA(J)
	  ENDIF
	ENDDO
	NHUNT(2) = 1
	DX(1) = .05*X(1)
	DX(2) = .1*X(2)
	DX(3) = 1.E-5
	CALL FITDATA(X,SUMSQS)
	PRINT*,' '
	PRINT*,' '
	PRINT*,' FIRST TRY AT FULL FIT, -DATA,',NPT,'  POINTS, SUMSQ='
     1  ,SUMSQS
	DO IT = 1,20
	  CALL HUNTMN(3,X,DX,Y,FITDATA,SUMSQ)
	  PRINT*,'HUNTMN OUTPUT',X(1),X(2),X(3),SUMSQ
	  IF(SUMSQ.GT..9995*SUMSQS) GO TO 45
	  SUMSQS = SUMSQ
	ENDDO
	PRINT*,'NOT ENOUGH ITERATIONS, POOR FULL FIT'
 45	CONTINUE
	VZERO = 10.**(X(3)/X(1))
	IF(NPT.GE.4) RMS = SQRT(SUMSQ/(NPT-3))
	PRINT*,' VZERO = OLD X(3) =',VZERO,'  RMS ERR',RMS,'  TM UNITS'
	PRINT*,'HUNTMN OUTPUT',X(1),X(2),X(3),SUMSQ
	SUMSQS = SUMSQ
	  DO IT = 1,10
	    CALL FINMIN(3,X,DX,Y,FITDATA,SUMSQF)
	    CALL FINMIN(3,Y,DX,X,FITDATA,SUMSQF)
	    PRINT*,' FINMIN GIVES',(X(N),N=1,3),SUMSQF
	    IF(SUMSQF.GT..9999*SUMSQS) GO TO 48
	    SUMSQS = SUMSQF
	  ENDDO
  48	  PRINT*,' '
	  PRINT*,' FINMIN GIVES',(X(N),N=1,3),SUMSQF
	  V71 = 10.**((71.+X(3))/X(1)) - X(2)
	  PRINT*,'VOLTS TO GIVE TM = 71', V71
	  PRINT*,' '
	CALL CALPLOT
	STOP
	END
	SUBROUTINE CALPLOT
C
	CHARACTER*18 TITLE(10)
	CHARACTER*120 JUNK(10)
	CHARACTER*120 STR
	COMMON /HEADBL/ JUNK
	COMMON /CALBLK/ VOLTS(2048),NDATA(2060)
	COMMON /PLOTDT/ ITERM,TITLE
C
	COMMON /MONGOPAR/
     1  X1,X2,Y1,Y2,GX1,GX2,GY1,GY2,LX1,LX2,LY1,LY2,
     1  GX,GY,CX,CY,
     1  EXPAND,ANGLE,LTYPE,LWEIGHT,
     1  CHEIGHT,CWIDTH,CXDEF,CYDEF,PSDEF,PYDEF,COFF,
     1  TERMOUT,XYSWAPPED,NUMDEV,
     1  PI,USERVAR(10),AUTODOT
	INTEGER*4 LX1,LX2,LY1,LY2,LTYPE,LWEIGHT,NUMDEV
C
	DIMENSION YY(500),PP(500)
C
	  CALL MGOINIT
	  CALL MGOSETUP(ITERM)
	  CALL MGOERASE
	IF(ITERM.LT.0) THEN
	  CALL MGOSETLOC(400.,330.,3000.,2300.)
	ENDIF
C
C	  WRITE(STR,703) NYR,NDOY,NHRMIN,NFCLKS
C 703	  FORMAT('\\t',I5,' DAY',I4,1X,I4.4,'  S/C CLK=',I12)
	  CALL MGOSETEXPAND(.85)
	  IF(ITERM.GT.0) THEN
	    CALL MGOGRELOCATE(10.,0.)                      ! maxch, crt
	  ELSE
	    CALL MGOGRELOCATE(400.,50.)                      ! hardcopy
	  ENDIF
C	  CALL MGOPUTLABEL(53,STR,9)
c	  WRITE(STR,704) SAA()
c 704	  FORMAT('\\tSOLAR ASPECT',F6.1,' DEG.')
c	  CALL MGOPUTLABEL(55,STR,9)
	  CALL MGOSETEXPAND(1.)
C
	CALL MGOTICKSIZE(-1.,0.,2.,10.)
	YLOLIM = -9.
	CALL MGOSETLIM(-5.,YLOLIM-1.,1.,130.)
	CALL MGOGRID(0)
	CALL MGOSETLTYPE(1)
	CALL MGOGRID(1)
	CALL MGOSETLTYPE(0)
	CALL MGOSETEXPAND(.6)
	  DO I = 1,2048
	    IF(VOLTS(I).GT.0.) THEN
	      VOLTL = ALOG10(VOLTS(I))
	      DATPL = NDATA(I)
	      CALL MGORELOCATE(VOLTL,DATPL)
	      CALL MGOSETANGLE(45.)
	      CALL MGOPOINT(4,1)
	      CALL MGOSETANGLE(0.)
	    ELSE
	      VOLTL = ALOG10(-VOLTS(I))
	      DATPL = -NDATA(I)
	      DATPL = AMAX1(DATPL,YLOLIM)
	      CALL MGORELOCATE(VOLTL,DATPL)
	      CALL MGOPOINT(4,1)
	    ENDIF
	  ENDDO
	CALL MGOSETEXPAND(1.)
	  CALL MGOBOX(1,2)
C	  WRITE(STR,704) JHP,PFRFREQ(JHP+1)
C 704	  FORMAT('\\tFREQ',I3,F6.2,'\\t kHz')
	CALL MGOSETEXPAND(.8)
	  CALL MGOXLABEL(5,'VOLTS')
	  CALL MGOYLABEL(10,'T/M NUMBER')
	  CALL MGOSETEXPAND(.55)
	  CALL MGORELOCATE(-4.9,129.)
	  CALL MGOLABEL(120,JUNK(2))
	  CALL MGORELOCATE(-4.9,119.)
	  CALL MGOLABEL(120,JUNK(1))
	  CALL MGORELOCATE(-4.5,113.)
	  CALL MGOLABEL(18,TITLE(2))
	  CALL MGORELOCATE(-4.5,107.)
	  CALL MGOLABEL(18,TITLE(3))
	  CALL MGORELOCATE(-4.5,101.)
	  CALL MGOLABEL(18,TITLE(4))
	  CALL MGORELOCATE(-4.5,95.)
	  CALL MGOLABEL(18,TITLE(5))
	  CALL MGORELOCATE(-4.5, 89.)
	  CALL MGOLABEL(18,TITLE(6))
	  CALL MGORELOCATE(-4.5, 81.)
	  CALL MGOLABEL(18,TITLE(7))
	  CALL MGORELOCATE(-4.6, 75.)
	  CALL MGOLABEL(12,'INPUT POS. ')
	  CALL MGOSETANGLE(45.)
	  CALL MGOPOINT(4,1)
	  CALL MGOSETANGLE(0.)
	  CALL MGORELOCATE(-4.6, 71.)
	  CALL MGOLABEL(12,'INPUT NEG. ')
	  CALL MGOPOINT(4,1)
	CALL MGOSETEXPAND(1.)
C
C
	CALL MGOSETEXPAND(.8)
	CALL MGOPLOTID(TITLE(8),'SINCAL')
	CALL MGOSETEXPAND(1.)
	IF(ITERM.LT.0) THEN
C	  CALL MGOTIDLE
	  CALL MGOPRNTPLOT(NVEC)
	  PRINT*,' NO. VECTORS PLOTTED',NVEC
	  NPLOTS = NPLOTS + 1
	ELSE
	  CALL MGOTCLOSE
	ENDIF
C
	RETURN
C
	END
	SUBROUTINE FITDATA(X,SUMSQ)
C
C	FITS T/M NO. = X(1)*ALOG10( VOLTS + VOFFS)/X(3))
C		WITH VOFFS = X(2)
C
	COMMON /FITBLK/ NPT,VOLTS(2048),TM(2048)
	DIMENSION X(25),WT(2048)
	DATA WT /2048*1./
C
	SUMSQ = 0.
c	X(3) = AMAX1(X(3),1.E-6)
C	X(2) = AMAX1(X(2),0.)
	DO N = 1,NPT
	  TMA = ABS(TM(N))
	  VT = ABS(VOLTS(N))
	  IF((VT+X(2)).GT.0.) THEN
C	    SUMSQ = SUMSQ + WT(N)*(TMA - X(1)*ALOG10((VT + X(2))/X(3)))**2
	    SUMSQ = SUMSQ + WT(N)*(TMA - X(1)*ALOG10(VT + X(2)) + X(3))**2
	  ELSE
	    SUMSQ = SUMSQ*(1.+2./NPT)
	  ENDIF
	ENDDO
	RETURN
	END
	SUBROUTINE FITZERO(X,SUMSQ)
C
	COMMON /ZBLK/ NPT,ZCROSS(128)
	DIMENSION X(25),WT(2048)
	DATA WT /2048*1./
C
	SUMSQ = 0.
	DO N = 1,NPT
	  SUMSQ = SUMSQ + (ZCROSS(N) - X(1) - X(2)*(N-1))**2
	ENDDO
	RETURN
	END
	options/extend_source
!------------------------------------------------------------------------------
	integer*4	function	get_stream_name(stream)
! This routine gets the user's TM stream type specification.
!
	implicit	none
	character*(*)	stream
	common /nrblk/ nrem,NHRMN
	integer*4	iq,NREM,NHRMN

  6	format(1x,'Enter TM stream type [O=offline(default), R=realtime ]: ',$)
  5	format(q,a)
  4	format(1x,'enter number of events to find and process')
  3	format(q,i10)

 10	write(6,6)
	read(5,5,err=10,end=20) iq, stream

	if (iq .lt. 1) then
	   stream = 'offline'
	else if (stream(1:1) .eq. 'o' .or. stream(1:1) .eq. 'O') then
	   stream = 'offline'
	else if (stream(1:1) .eq. 'r' .or. stream(1:1) .eq. 'R') then
	   stream = 'realtime'
	else
	   ! assume the user entered the name of an offline file
	end if

	get_stream_name = 1

        if(stream.eq.'offline') then
	  write(6,*)  'type hr,min to start, e.g. 0412'
	  read(5,3,err=10,end=20) iq, NHRMN
	  type*,NHRMN
C	  write(6,4)
C	  read(5,3,err=10,end=20) iq, NREM
C	  type*,nrem
	endif

 20	return
	end
	SUBROUTINE OLDATA
	COMMON /CALBLK/ VOLTS(2048),NDATA(2060)
	CHARACTER*120 JUNK(10)
	DO IL = 1,6
	  READ(19,1002) JUNK(IL)
	  PRINT *,IL,JUNK(IL)
	ENDDO
C
	DO IL = 1,103
	  IK2 = 20*IL
	  IK1 = IK2-19
	  READ(19,1001,END=200) LINNO,(NDATA(IK),IK=IK1,IK2)
C	  PRINT*,'DATA',LINNO,NDATA(IK1),NDATA(IK2)
	ENDDO
C
 1001	FORMAT(I4,1X,20I5)
 1002	FORMAT(A)
  200	RETURN
	END

ext: for