text:
MODULE Viscid

  USE SIMPLE
  USE SPECIAL
  USE GRID

  IMPLICIT NONE

CONTAINS

  !>
  !! @brief Compute the curvilinear cartesian viscous fluxes in 1 dimension
  !!
  !! retains logic to simplify calculations on uniform or stretched cartesian grids
  !!
  !! See /ref conserve for theory
  !!
  SUBROUTINE StrongFlux1D                                            &
       (numDim, fluxDir, gridSizes, numPoints, opInterval, gridType, &
       gridMetric, tauBuffer, energyBuffer, fluxBuffer)

    IMPLICIT NONE

    INTEGER(KIND=4), INTENT(IN)         :: numDim, fluxDir, gridType
    INTEGER(KIND=8), INTENT(IN)         :: gridSizes(numDim),numPoints
    INTEGER(KIND=8), INTENT(IN)         :: opInterval(2*numDim)
    REAL(KIND=8),    INTENT(IN), TARGET :: gridMetric(numDim*numDim*numPoints)
    REAL(KIND=8),    INTENT(IN), TARGET :: tauBuffer(numPoints*numDim*(numDim+1)/2)
    REAL(KIND=8),    INTENT(IN), TARGET :: energyBuffer(numPoints*numDim)
    REAL(KIND=8),    INTENT(OUT),TARGET :: fluxBuffer(numPoints*(numDim+2))

    INTEGER         :: iDim, tensorIndex, iVel
    INTEGER(KIND=8) :: fluxOffset, tensorOffset, dirOffset, metricOffset
    INTEGER(KIND=8) :: velIndex, fluxIndex, velOffset, heatOffset, dimOffset
    REAL(KIND=8)    :: gridScale
    REAL(KIND=8)    :: minusOne

    REAL(KIND=8),    DIMENSION(:), POINTER :: bufPtr
    REAL(KIND=8),    DIMENSION(:), POINTER :: metricPtr
    REAL(KIND=8),    DIMENSION(:), POINTER :: fluxPtr
    REAL(KIND=8),    DIMENSION(:), POINTER :: tauPtr
    REAL(KIND=8),    DIMENSION(:), POINTER :: energyPtr

    INTEGER, DIMENSION(2,2) :: index2D
    INTEGER, DIMENSION(3,3) :: index3D

    ! Map 0-based index for (2,3)-dimensional symmetric tensor
    index2D = RESHAPE((/ 0, 1, 1, 2 /), SHAPE(index2D))
    index3D = RESHAPE((/ 0, 1, 2, 1, 3, 4, 2, 4, 5 /), SHAPE(index3D))

    dirOffset = (fluxDir-1)*numPoints
    fluxOffset = 0

    IF(gridType < RECTILINEAR) THEN

       gridScale  =  gridMetric(fluxDir)

       ! Momentum terms
       DO iDim = 1, numDim

          IF(numDim == 2) THEN
             tensorIndex = index2D(fluxDir,iDim)
          ELSE
             tensorIndex = index3D(fluxDir,iDim)
          ENDIF

          tensorOffset =  tensorIndex*numPoints
          fluxOffset   =  fluxOffset + numPoints
          fluxPtr      => fluxBuffer(fluxOffset+1:fluxOffset+numPoints)
          tauPtr       => tauBuffer(tensorOffset+1:tensorOffset+numPoints)

          ! Grid metric is scalar and constant, this gets flux = metric * tau_(iDim)(
          CALL YAX(numDim,numPoints,gridSizes,opInterval,gridScale,tauPtr,fluxPtr)

       END DO

       ! Energy
       fluxOffset  =  fluxOffset + numPoints
       fluxPtr     => fluxBuffer(fluxOffset+1:fluxOffset+numPoints)
       energyPtr   => energyBuffer(dirOffset+1:dirOffset+numPoints)

       CALL YAX(numDim,numPoints,gridSizes,opInterval,gridScale,energyPtr,fluxPtr)

    ELSE IF (gridType == RECTILINEAR) THEN

       metricPtr    => gridMetric(dirOffset+1:dirOffset+numPoints)

       ! Momentum terms
       DO iDim = 1, numDim

          IF(numDim == 2) THEN
             tensorIndex = index2D(fluxDir,iDim)
          ELSE
             tensorIndex = index3D(fluxDir,iDim)
          ENDIF

          tensorOffset =  tensorIndex*numPoints
          fluxOffset   =  fluxOffset + numPoints
          fluxPtr      => fluxBuffer(fluxOffset+1:fluxOffset+numPoints)
          tauPtr       => tauBuffer(tensorOffset+1:tensorOffset+numPoints)

          ! Grid metric is scalar, this gets flux = metric * tau_(iDim)(
          CALL ZXY(numDim,numPoints,gridSizes,opInterval,metricPtr,tauPtr,fluxPtr)

       END DO

       ! Energy
       fluxOffset  =  fluxOffset + numPoints
       fluxPtr     => fluxBuffer(fluxOffset+1:fluxOffset+numPoints)
       energyPtr   => energyBuffer(dirOffset+1:dirOffset+numPoints)

       CALL ZXY(numDim,numPoints,gridSizes,opInterval,metricPtr,energyPtr,fluxPtr)

    ELSE ! must be curvilinear

       ! need row 'fluxDir' of the metric tensor
       metricOffset = (fluxDir-1)*numDim*numPoints

       ! For curvilinear
       ! Momentum terms flux_i  = (metric_fluxDir_j * tau_j_i) i=[1:numDim]
       DO iDim = 1, numDim

          fluxOffset   =  fluxOffset + numPoints
          fluxPtr      => fluxBuffer(fluxOffset+1:fluxOffset+numPoints)

          CALL ASSIGNMENTXA(numDim,numPoints,gridSizes,opInterval,0.0_8,fluxPtr)

          DO iVel = 1, numDim

             velOffset = (iVel-1)*numPoints+metricOffset
             metricPtr => gridMetric(velOffset+1:velOffset+numPoints)

             IF(numDim == 2) THEN
                tensorIndex = index2D(iVel,iDim)
             ELSE
                tensorIndex = index3D(iVel,iDim)
             ENDIF

             tensorOffset =  tensorIndex*numPoints
             tauPtr       => tauBuffer(tensorOffset+1:tensorOffset+numPoints)

             ! Get fluxPtr += (metric_fluxDir_iVel * tau_iVel_iDim)
             CALL YWXPY(numDim,numPoints,gridSizes,opInterval,metricPtr,tauPtr,fluxPtr)

          END DO

       END DO

       ! Energy flux is just (metric_fluxDir_i * Q_i), since
       ! both have proper vector storage the DOT operator
       ! can be used
       fluxOffset  =  fluxOffset + numPoints
       fluxPtr     => fluxBuffer(fluxOffset+1:fluxOffset+numPoints)
       metricPtr   => gridMetric(metricOffset+1:metricOffset+numDim*numPoints)

       CALL ZXDOTY(numDim,numPoints,gridSizes,opInterval,numDim,metricPtr,energyBuffer,fluxPtr)

    ENDIF

  END SUBROUTINE StrongFlux1D

  !>
  !! @brief Compute the curvilinear cartesian viscous fluxes in 1 dimension
  !!
  !! retains logic to simplify calculations on uniform or stretched cartesian grids
  !!
  !! See /ref conserve for theory
  !!
  SUBROUTINE ScalarFlux1D                                            &
       (numDim, fluxDir, gridSizes, numPoints, opInterval, gridType, &
       gridMetric, gradScalar, fluxBuffer)

    IMPLICIT NONE

    INTEGER(KIND=4), INTENT(IN)         :: numDim, fluxDir, gridType
    INTEGER(KIND=8), INTENT(IN)         :: gridSizes(numDim),numPoints
    INTEGER(KIND=8), INTENT(IN)         :: opInterval(2*numDim)
    REAL(KIND=8),    INTENT(IN), TARGET :: gridMetric(numDim*numDim*numPoints)
    REAL(KIND=8),    INTENT(IN), TARGET :: gradScalar(numPoints*numDim)
    REAL(KIND=8),    INTENT(OUT)        :: fluxBuffer(numPoints)

    INTEGER(KIND=8) :: dirOffset, metricOffset
    INTEGER(KIND=8) :: velIndex, fluxIndex, velOffset, heatOffset, dimOffset
    REAL(KIND=8)    :: gridScale
    REAL(KIND=8)    :: minusOne

    REAL(KIND=8),    DIMENSION(:), POINTER :: gradPtr
    REAL(KIND=8),    DIMENSION(:), POINTER :: metricPtr

    dirOffset = (fluxDir-1)*numPoints

    IF(gridType < RECTILINEAR) THEN

       gridScale    =  gridMetric(fluxDir)
       gradPtr      => gradScalar(dirOffset+1:dirOffset+numPoints)

       ! Grid metric is scalar and constant, this gets flux = metric * grad(Y)
       CALL YAX(numDim,numPoints,gridSizes,opInterval,gridScale,gradPtr,fluxBuffer)

    ELSE IF (gridType == RECTILINEAR) THEN

       metricPtr    => gridMetric(dirOffset+1:dirOffset+numPoints)
       gradPtr      => gradScalar(dirOffset+1:dirOffset+numPoints)

       ! Grid metric is scalar, this gets flux = metric * tau_(iDim)(
       CALL ZXY(numDim,numPoints,gridSizes,opInterval,metricPtr,gradPtr,fluxBuffer)

    ELSE ! must be curvilinear

       metricOffset = dirOffset*numDim
       metricPtr => gridMetric(metricOffset+1:metricOffset+numDim*numPoints)

       CALL ZXDOTY(numDim,numPoints,gridSizes,opInterval,numDim,metricPtr,gradScalar,fluxBuffer)

    ENDIF

  END SUBROUTINE ScalarFlux1D

END MODULE Viscid

max_stars_repo_name: MTCam/OKLib

max_issues_repo_name: MTCam/OKLib

max_forks_repo_name: MTCam/OKLib


text:
C     The purpose of ths test it to test the compilers ability to differentiate between functions and variables of the same name
C     This is made extra interesting by using a keyword as the variable / function name
      DO(X) = X + 5
      PRINT DO(5)
      X = 50
C     Expected outcome is to print 65
      PRINT X + DO(10)
      END

ext: f

max_stars_repo_name: ElliotAlexander/EDSAC-FORTRAN-Compiler

max_issues_repo_name: ElliotAlexander/EDSAC-FORTRAN-Compiler

max_forks_repo_name: ElliotAlexander/EDSAC-FORTRAN-Compiler


text:
diff --git a/eg.f90 b/eg.f90
index d15ab90..5a99559 100644
--- a/eg.f90
+++ b/eg.f90
@@ -33,7 +33,7 @@ PROGRAM test
   CHARACTER(26) :: name, greet
   INTEGER :: value

-  greet = "Hello"
+  greet = "Hello and welcome"
   PRINT*, greet

   name = "abcde"

max_stars_repo_name: WarwickRSE/PracticalGit

max_issues_repo_name: WarwickRSE/PracticalGit

max_forks_repo_name: WarwickRSE/PracticalGit


text:
	PROGRAM PROCIPSHK
C
C	READS PROCIPSHK.RESULTS AND ADDS SOME INFORMATION
C	IF RTEMP LT 0, IT MAKES A PLOT LIKE MAKEFILE17.  IF RTEMP = ACTUAL
C	VALUE, IT DOES NOTHING EXCEPT REWRITE INPUT LINE TO FOR091.DAT
C	IQUAL = 0 MEANS TIME HAS BEEN DETERMINE FROM LISTS OR NO SHOCK
C	IQUAL = 1 MEANS TIME HAS BEEN DETERMINED FROM SUMMARY PLOTS
C	IQUAL = 2 MEANS TIME HAS BEEN DETERMINED FROM SPLOT
C	IQUAL = 3 MEANS ANGBN HAS BEEN DETERMINED MANUALLY
C	The program also writes a record of TDS events, to for064.dat.
C	afterward, running ORDER64 puts them in chronological order into
C	for065.dat
C
	INTEGER YYMMDD,HHMMSS,YYYY,HH,MM,DDW,OK,OKT,TMCHANNEL,CH
	INTEGER YYYYW,YYMMDDW,HHW,MMW,NDAY
	INTEGER W_CHANNEL_OPEN
	integer*4 OK1,OK2,SCETI4(2),TDS_COUNT
	INTEGER*4 W_CHANNEL_CLOSE,W_EVENT,RET_SIZE,W_MESSAGES_OFF
	INTEGER*4 W_ITEM_I4,W_ITEM_R4,W_ITEM_R8
	REAL BX(50),BY(50),BZ(50),YY3D(50),VX(50),VY(50),VZ(50)
	REAL VTEMP(3),VTEMP1(3),V1(3),V2(3),V3(3),SHNORM(3),DELV(3)
	REAL*8 SCET8,SCETSV,RGSE,SCETTEMP,SCETT(50),SCETINT
	character*32 ITEM
	character*4 event
	character*80	STREAM
	COMMON /FILEDAY/ STREAM
C
	DATA NERR /0/
	DATA INOUT /0/
	DATA RE /6.378E3/
	DATA TWOPI /6.2831853/
	data event /'CDF'/
	DATA INOUT /1/
	DATA IQUAL /0/
C
	OPEN(UNIT=69,FILE='IPSHK.RESULTS',STATUS='OLD',READONLY)
	WRITE(91,*) ' '
	WRITE(91,*) ' '
	WRITE(91,*) ' '
 100	CONTINUE
C
C	READ A NEW LINE FROM THE SHOCK CROSSING LIST
C
	READ(69,*,END=300) SCETDAY,YYMMDD,HHMMSS,XTEMP,YTEMP,ZTEMP,
     1		RTEMP,INOUT
C	IF(SCETDAY.GT.4891.) STOP
	IQUAL = 0
c*********
C	rtemp = -1.
c*********
	print*,' '
	print*,'input time',scetday,yymmdd,hhmmss
	YYYY = YYMMDD/10000.
	MON = (YYMMDD - 10000*YYYY)/100
	NDAY = MOD(YYMMDD,100)
	HH = HHMMSS/10000.
	MM = (HHMMSS - 10000.*HH)/100.
	SS = MOD(HHMMSS,100)
C*******	to correct an error in year 2000 data
	if(yyyy.eq.2000) scetday = scetday - 1.
	SCET8 = DBLE(SCETDAY) + HH/24.D00 + MM/1440.E00 + SS/86400.D00
	print*,'calc scet8,yymmdd ',scet8,yymmdd
C	call w_ur8_to_ymd(scet8,yyyyw,monw,ddw,hhw,mmw,ssw,ms)
C	call w_ur8_from_ymd(scet8,yyyyw,monw,ddw,hhw,mmw,ssw,ms)
C	PRINT*,'return Y M D H M',YYYYW,MONW,DDW,HHW,MMW
C	SCET8 = DBLE(SCETDAY) + HHW/24.D00 + MMW/1440.E00 + SSW/86400.D00
C
C	OPEN A WIND DATA CHANNEL
C
	WRITE(STREAM,30) YYMMDD
 30	FORMAT('WI_LZ_WAV_',I8.8,'_V*.DAT')
	PRINT*,YYMMDD,STREAM
C
C	call w_ur8_to_ydoy(scet8,yyyy,idoy,msec)
C
	print*,'in main, stream=',stream
C
	ok = w_channel_open(tmchannel,stream)
	ch = tmchannel
C
	print*,'open, TMCHANNEL,ch,ok=',TMCHANNEL,ch,ok
	if (ok.ne.1) then
	  print*, 'Cannot open t/m channel'
	endif
c
	scetsv = 0.d00
	print*,'calling channel position 1',scetsv
	ok = w_channel_position(tmchannel,scetsv)
	print*,'file starts at ',scetsv
	okt = w_channel_position(tmchannel,scet8)
C
 	ok = w_event(ch,event)
	ITEM = 'WIND_ORBIT_X(GSE)_R8'
	ok = W_ITEM_R8(ch, item, RGSE, 1, ret_size)
	XGSE = RGSE/RE
	ITEM = 'WIND_ORBIT_Y(GSE)_R8'
	ok = W_ITEM_R8(ch, item, RGSE, 1, ret_size)
	YGSE = RGSE/RE
	ITEM = 'WIND_ORBIT_Z(GSE)_R8'
	ok = W_ITEM_R8(ch, item, RGSE, 1, ret_size)
	ZGSE = RGSE/RE
	RRE = SQRT(XGSE**2 + YGSE**2 + ZGSE**2)
	RTEST = RTEMP/RE
C
	IF(RTEMP.LT.0.) GO TO 200
C
	okt = w_channel_position(tmchannel,scet8)
 	ok = w_event(ch,event)
	print 1107, scetday,yymmdd,hhmmss,xtemp,ytemp,ztemp,rtemp,inout
     1		,iqual,angbn,sumsq,VALF,ALFMACH,VRAT,BRAT,DELB,
     1		TDS_COUNT,TDSM,PWR
	WRITE(91,1107) SCETDAY,YYMMDD,HHMMSS,XTEMP,YTEMP,ZTEMP,RTEMP,INOUT
     1       ,IQUAL,ANGBN,SUMSQ,VALF,ALFMACH,VRAT,BRAT,DELB,
     2	      TDS_COUNT,TDSM,PWR,ISPS
 1107	FORMAT(F8.0,I9,I7,3F8.2,F8.2,2I3,F4.0,F6.3,F6.0,F6.2,3F5.1,
     1	 I4,F6.1,E9.2,I3)
c	at present, 122 characters above
C
c	IF(ABS(RTEST-RRE).LT.1.E-2) GO TO 100
C
	OKT = W_MESSAGES_OFF(ch)
C
 	ok = w_event(ch,event)
C
	ITEM = 'EVENT_SCET'
	ok = W_ITEM_I4(ch, item, SCETI4, 2, ret_size)
C
	ITEM = 'WIND_ORBIT_X(GSE)_R8'
	ok = W_ITEM_R8(ch, item, RGSE, 1, ret_size)
	XGSE = RGSE/RE
	ITEM = 'WIND_ORBIT_Y(GSE)_R8'
	ok = W_ITEM_R8(ch, item, RGSE, 1, ret_size)
	YGSE = RGSE/RE
	ITEM = 'WIND_ORBIT_Z(GSE)_R8'
	ok = W_ITEM_R8(ch, item, RGSE, 1, ret_size)
	ZGSE = RGSE/RE
	ITEM = 'EVENT_SCET_R8'
	OK = W_ITEM_R8(ch, item, SCET8, 1, ret_size)
	SCETINT = DINT(SCET8)
C	print*,'SCET8,SCETINT',SCET8,SCETINT
	RRE = SQRT(XGSE**2 + YGSE**2 + ZGSE**2)
	RTEST = RTEMP/RE
c	go to 200
	ok = w_channel_close(TMCHANNEL)
C	INOUT = MOD(INOUT+1,2)
	GO TO 100
C
 200	CONTINUE
C
	ITEM = 'WIND_ORBIT_X(GSE)_R8'
	ok = W_ITEM_R8(ch, item, RGSE, 1, ret_size)
	XGSE = RGSE/RE
	ITEM = 'WIND_ORBIT_Y(GSE)_R8'
	ok = W_ITEM_R8(ch, item, RGSE, 1, ret_size)
	YGSE = RGSE/RE
	ITEM = 'WIND_ORBIT_Z(GSE)_R8'
	ok = W_ITEM_R8(ch, item, RGSE, 1, ret_size)
	ZGSE = RGSE/RE
	RRE = SQRT(XGSE**2 + YGSE**2 + ZGSE**2)
C
C	AVERAGE B FOR MINUP MINUTES UPSTREAM OF SHOCK
C
	MINUP = 7
	UPN = 2*(MINUP+2)*60./46.
	NPTUPN = UPN
	BXAVR = 0.
	BYAVR = 0.
	BZAVR = 0.
	DENS = 0.
	SAMPL = 1.E-9
	scettemp = scet8
C
C	BACKUP IF THIS IS AN INBOUND PASS, LEAVE TWO MINUTES BLANK
C
C		INTERVAL FOR MFI DATA IS 46 SEC
C	IF(inout.eq.1) then
	  scettemp = scet8 - (minup+2)*1.d00/1440.d00
C	ELSE
C	  scettemp = scet8 + 2.d00/1440.d00
C	ENDIF
	okt = w_channel_position(tmchannel,scettemp)
C
C	GET BOTH DOWNSTREAM AND UPSTREAM DATA
C
 	ok = w_event(ch,event)
 	ITEM = 'WIND_MFI_SCET_R8'
	  ok = W_ITEM_R8(ch, item, SCETT, NPTUPN, ret_size)
	  IF(OK.NE.1) THEN
	    PRINT*,'PROGRAM FAILED TO FIND CDF DATA,OK,YYMMDD=',OK,YYMMDD
	    WRITE(91,1107) SCETDAY,YYMMDD,HHMMSS,XTEMP,YTEMP,ZTEMP,RTEMP,INOUT
     1       ,IQUAL
	    GO TO 100
	  ENDIF
	print*,'after scett item, nptupn,ret_size=', scett(1),scett(RET_SIZE),
     1   nptupn,ret_size
	NPTUPN = RET_SIZE
 	ITEM = 'WIND_MFI_BX(GSE)_R4'
	  ok = W_ITEM_R4(ch, item, BX, NPTUPN, ret_size)
	ITEM = 'WIND_MFI_BY(GSE)_R4'
	  ok = W_ITEM_R4(ch, item, BY, NPTUPN, ret_size)
	ITEM = 'WIND_MFI_BZ(GSE)_R4'
	  ok = W_ITEM_R4(ch, item, BZ, NPTUPN, ret_size)
	ITEM = 'WIND_3DP_ION_DENSITY_R4'
	ok = W_ITEM_R4(ch, item, YY3D, NPTUPN, ret_size)
	ITEM = 'WIND_3DP_ION_VX(GSE)_R4'
	ok = W_ITEM_R4(ch, item, VX, NPTUPN, ret_size)
	ITEM = 'WIND_3DP_ION_VY(GSE)_R4'
	ok = W_ITEM_R4(ch, item, VY, NPTUPN, ret_size)
	ITEM = 'WIND_3DP_ION_VZ(GSE)_R4'
	ok = W_ITEM_R4(ch, item, VZ, NPTUPN, ret_size)
	print*,'B data asked for and got ',nptupn,ret_size
C	print*,'bx',(bx(I),I=1,NPTUPN)
C	print*,'by',(by(I),I=1,NPTUPN)
C	print*,'bz',(bz(I),I=1,NPTUPN)
	print*,'ST',(scett(I),I=1,NPTUPN)
C
C	UPSTREAM AVERAGE
C
	BXAVR = 0.
	BYAVR = 0.
	BZAVR = 0.
	VXAVR = 0.
	VYAVR = 0.
	VZAVR = 0.
	SUMSQ = 0.
	DENS = 0.
	COUNT = 1.E-9
	  IF(INOUT.EQ.0) THEN
C		IT'S OUTBOUND, AVERAGE FROM SHOCK TIME + 2 MIN TO + MINUP+2
	    PPBN1 = DMOD(SCET8,1.d00)*24.D00 + 2.D00/60.D00
	    PPBN2 = PPBN1 + DFLOAT(MINUP)/60.D00
	  ELSE
	    PPBN2 = DMOD(SCET8,1.D00)*24.D00 - 2.D00/60.D00
	    PPBN1 = PPBN2 - DFLOAT(MINUP)/60.D00
	  ENDIF
	SCETINT = DINT(SCET8)
	print*,'ppbn',scet8,SCETINT,ppbn1,ppbn2
	time = DMOD(scet8,1.D00)*24.d00
	N = 0
	DOWHILE(TIME.LT.PPBN1.AND.N.LE.NPTUPN)
	  N = N+1
	  IF(YY3D(N).GT.0.) DENS0 = YY3D(N)
	  time = DMOD(scett(n),1.D00)*24.d00
	ENDDO
	DO N = 1,NPTUPN
	  time = DMOD(scett(n),1.D00)*24.d00
	  IF(TIME.LE.PPBN2.AND.TIME.GE.PPBN1) THEN
C	  print*,'pp check up',time,ppbn2,n,nptupn
	    COUNT = COUNT+1.
	    BXAVR = BXAVR + BX(N)
	    BYAVR = BYAVR + BY(N)
	    BZAVR = BZAVR + BZ(N)
	    SUMSQ = SUMSQ + BX(N)**2 + BY(N)**2 + BZ(N)**2
	    IF(YY3D(N).GT.0.) THEN
	      DENS = DENS + YY3D(N)
	      DENS0 = YY3D(N)
	    ELSE
	      DENS = DENS + DENS0
	    ENDIF
	    VXAVR = VXAVR + VX(N)
	    VYAVR = VYAVR + VY(N)
	    VZAVR = VZAVR + VZ(N)
	  ENDIF
	ENDDO
	BXAVR = BXAVR/COUNT
	BYAVR = BYAVR/COUNT
	BZAVR = BZAVR/COUNT
	VXAVR = VXAVR/COUNT
	VYAVR = VYAVR/COUNT
	VZAVR = VZAVR/COUNT
	DENS = DENS/COUNT
	V1(1) = BXAVR
	V1(2) = BYAVR
	V1(3) = BZAVR
C	IF(DENS.LT.0.) DENS = 7.
	SQAVR = BXAVR**2 + BYAVR**2 + BZAVR**2
	SUMSQ = SUMSQ/COUNT - SQAVR
	print*,'bAVR',bXAVR,BYAVR,BZAVR,COUNT
	IF(SQAVR.NE.0.) SUMSQ = SUMSQ/SQAVR
	BMAG = SQRT(BXAVR**2+BYAVR**2+BZAVR**2)
C	BDOTN = 0.
C	BDOTN = ABS(BDOTN)
C	ANGBN = ACOSD(BDOTN)
C	print*,'bmag,angbn,SUMSQ 1 ',bmag,angbn,SUMSQ
	VALF = 0.
	IF(DENS.GT.0.) VALF = 21.9*BMAG/SQRT(DENS)  !Bmag in nT, dens in cu.cm.
C
	VDOTN = 0.
C	VDOTN = (VXAVR*SNX + VYAVR*SNY + VZAVR*SNZ)/SNORM
C	VDOTN = ABS(VDOTN)
	ALFMACH = 0.
	IF(VALF.NE.0.) ALFMACH = VDOTN/VALF
C
C	DOWNSTREAM AVERAGE
C
	BXAVRD = 0.
	BYAVRD = 0.
	BZAVRD = 0.
	VXAVRD = 0.
	VYAVRD = 0.
	VZAVRD = 0.
	SUMSQD = 0.
	COUNT = 1.E-9
	  IF(INOUT.EQ.0) THEN
C		IT'S OUTBOUND, AVERAGE FROM SHOCK TIME MINUP-2 MIN TO -2
	    PPBN2 = DMOD(SCET8,1.D00)*24.D00 - 2.D00/60.D00
	    PPBN1 = PPBN2 - DFLOAT(MINUP)/60.D00
	  ELSE
	    PPBN1 = DMOD(SCET8,1.d00)*24.D00 + 2.D00/60.D00
	    PPBN2 = PPBN1 + MINUP/60.D00
	  ENDIF
	SCETINT = DINT(SCET8)
	print*,'ppbn',scet8,SCETINT,ppbn1,ppbn2
	DO N = 1,NPTUPN
	  time = DMOD(scett(n),1.D00)*24.d00
	  IF(TIME.LE.PPBN2.AND.TIME.GE.PPBN1) THEN
C	  print*,'pp check down',time,ppbn2,n,nptupn
	    COUNT = COUNT+1.
	    BXAVRD = BXAVRD + BX(N)
	    BYAVRD = BYAVRD + BY(N)
	    BZAVRD = BZAVRD + BZ(N)
	    SUMSQD = SUMSQ + BX(N)**2 + BY(N)**2 + BZ(N)**2
	    VXAVRD = VXAVRD + VX(N)
	    VYAVRD = VYAVRD + VY(N)
	    VZAVRD = VZAVRD + VZ(N)
	  ENDIF
	ENDDO
	BXAVRD = BXAVRD/COUNT
	BYAVRD = BYAVRD/COUNT
	BZAVRD = BZAVRD/COUNT
	VXAVRD = VXAVRD/COUNT
	VYAVRD = VYAVRD/COUNT
	VZAVRD = VZAVRD/COUNT
C	IF(DENS.LT.0.) DENS = 7.
	SQAVRD = BXAVRD**2 + BYAVRD**2 + BZAVRD**2
	SUMSQD = SUMSQD/COUNT - SQAVR
	print*,'bAVRD',bXAVRD,BYAVRD,BZAVRD,COUNT
	IF(SQAVRD.NE.0.) SUMSQD = SUMSQD/SQAVRD
	BMAGD = SQRT(BXAVRD**2+BYAVRD**2+BZAVRD**2)
	DELB = SQRT((BXAVRD-BXAVR)**2+(BYAVRD-BYAVR)**2+(BZAVRD-BZAVR)**2)
C
C	COPLANARITY
C
	SHNORM(1) = BYAVR*BZAVRD - BZAVR*BYAVRD
	SHNORM(2) = BZAVR*BXAVRD - BXAVR*BZAVRD
	SHNORM(3) = BXAVR*BYAVRD - BYAVR*BXAVRD
	SNORM = SQRT(SHNORM(1)**2 + SHNORM(2)**2 + SHNORM(3)**2)
	if(snorm.ne.0.) print*,'shock normal',shnorm(1)/snorm,
     1  shnorm(2)/snorm,shnorm(3)/snorm
C
C	IF MAGNETIC FIELDS NEARLY PARALLEL, USE VELOCITY FOR NORMAL
C
	PRINT*,'SHNORM TEST',BMAG,BMAGD,SNORM
	IF(SNORM.LE..1*SQRT(BMAG*BMAGD)) THEN
	  SHNORM(1) = VXAVR + VXAVRD
	  SHNORM(2) = VYAVR + VYAVRD
	  SHNORM(3) = VZAVR + VZAVRD
	  SNORM = SQRT(SHNORM(1)**2 + SHNORM(2)**2 + SHNORM(3)**2)
	ENDIF
C
C	BDOTN = 0.
	SNX = SHNORM(1)
	SNY = SHNORM(2)
	SNZ = SHNORM(3)
	SNORM = SQRT(SNX**2 + SNY**2 + SNZ**2)
	IF(SNORM*BMAG.NE.0.) print*,'snorm*bmag',snorm*bmag
	IF(SNORM*BMAG.NE.0.)
     1	    BDOTN = (BXAVR*SNX + BYAVR*SNY + BZAVR*SNZ)/(SNORM*BMAG)
C	BDOTN = ABS(BDOTN)
C	ANGBN = ACOSD(BDOTN)
C	print*,'bmag,angbn,SUMSQ',bmag,angbn,SUMSQ
C	if(1) stop
	VALF = 0.
	IF(DENS.GT.0.) VALF = 21.9*BMAG/SQRT(DENS)  !Bmag in nT, dens in cu.cm.
	IF(SNORM.EQ.0.) GO TO 220
C
	VTEMP(1) = BXAVR
	VTEMP(2) = BYAVR
	VTEMP(3) = BZAVR
	SHNORM(1) = SNX/SNORM
	SHNORM(2) = SNY/SNORM
	SHNORM(3) = SNZ/SNORM
	CALL CROSSP(VTEMP,SHNORM,VTEMP1,BTAN)
C	print*,'up',vtemp
C	print*,'up',shnorm
C	print*,'up',v1
C	print*,'up',btan
	VTEMP(1) = BXAVRD
	VTEMP(2) = BYAVRD
	VTEMP(3) = BZAVRD
	CALL CROSSP(VTEMP,SHNORM,VTEMP1,BTAND)
C	print*,'dn',vtemp
C	print*,'dn',shnorm
C	print*,'dn',vtemp1
C	print*,'dn',btand
	BRAT = 0.
	IF(BTAN.NE.0.) BRAT = BTAND/BTAN
	print*,'vt,sh,vt1,bt,btd',vtemp(1),shnorm(1),vtemp1(1),btan,btand
C
C	ASSUME DELV IS NORMAL TO SHOCK
C
	DELV(1) = VXAVR-VXAVRD
	DELV(2) = VYAVR-VYAVRD
	DELV(3) = VZAVR-VZAVRD
c	print*,'avr  ',vxavr,vyavr,vzavr
c	print*,'avrd ',vxavrd,vyavrd,vzavrd
c	print*,'delv ',delv
	ANGBN = -1.
	IF(ABS(DELV(1)).LT.1.E4.AND.ABS(DELV(2)).LT.1.E4.AND.ABS(DELV(3))
     1	  .LT.1.E4) THEN
	  SNORM = SQRT(DELV(1)**2 + DELV(2)**2 + DELV(3)**2)
	  SNX = DELV(1)/SNORM
	  SNY = DELV(2)/SNORM
	  SNZ = DELV(3)/SNORM
	  IF(SNORM*BMAG.NE.0.)
     1	    BDOTN = (BXAVR*SNX + BYAVR*SNY + BZAVR*SNZ)/(SNORM*BMAG)
	  BDOTN = ABS(BDOTN)
	  ANGBN = ACOSD(BDOTN)
	ENDIF
	print*,'bmag,angbn,SUMSQ',bmag,angbn,SUMSQ
C
	VDOTND = (VXAVRD*SNX + VYAVRD*SNY + VZAVRD*SNZ)
	VDOTND = ABS(VDOTND)
	VRAT = 0.
	IF(VDOTND.NE.0.) VRAT = VDOTN/VDOTND
C
C	ALFMACH = 0.
C	IF(VALF.NE.0.) ALFMACH = VDOTN/VALF
C
c	IQUAL=1
C
 220	CONTINUE
c	RTEMP = RRE
	IF(RTEMP.LT.-.5) CALL
     1	 SPLOT(TMCHANNEL,SCET8,ANGBN,MINUP,TDS_COUNT,PWR,TDSM,ISPS,
     1	  SHNORM)
	print 1107, scetday,yymmdd,hhmmss,XGSE,YGSE,ZGSE,RRE,inout
     1	  ,iqual,angbn,sumsq,VALF,ALFMACH,VRAT,BRAT,DELB,
     1		TDS_COUNT,TDSM,PWR,ISPS
C	print*,'delb',delb
	WRITE(91,1107) SCETDAY,YYMMDD,HHMMSS,XGSE,YGSE,ZGSE,RRE,INOUT
     1       ,IQUAL,ANGBN,SUMSQ,VALF,ALFMACH,VRAT,BRAT,DELB,
     2	      TDS_COUNT,TDSM,PWR,ISPS
c	okt = w_channel_position(tmchannel,scet8)
C	if(1) stop
	ok = w_channel_close(TMCHANNEL)
	GO TO 100
 300	CONTINUE
	STOP
	END
	SUBROUTINE SPLOT(CDFCH,SCET,ANGBN,MINUP,TDSC,PWR,TDSM,ISPS,SHNORM)
C
C	PLOT THINGS TO IDENTIFY SHOCK POSITION, NAMELY B, DENSITY,
C		SW VELOCITY
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
	COMMON /FILEDAY/ TODAY
	character*32 ITEM
	character*80 TODAY,STR
	character*4 event
	INTEGER*4 CH,NDAY,ret_size,ok
	INTEGER*4 W_EVENT,W_CHANNEL_OPEN,W_CHANNEL_CLOSE,CDFCH,TDSCH
	INTEGER*4 NDATA(2050),SCETI4(2),NEVT,TDS_CHANNEL,TDS_COUNT,TDSC
	INTEGER*4 BCOUNT,first_channel,last_channel,fft_channel
	REAL YY(600),PP(600),YY3D(600),YY3S(600),YYSD(600),YYSS(600)
	REAL PP1(600),XRE(600),YRE(600),ZRE(600),VDATA(2050),SPECT(1025)
	REAL BDOTN(600),EFFLEN(9),SHNORM(3)
	REAL*8 SCET,SCET1,SCET2,SCSV(600),SCETT,RGSE,DNDAY,SCETTDS1,SCETTDS2
	REAL*8 SCETINT,SCETFFT
	REAL*4 TDSSMIN(100),TDSSMAX(100),TDSSTIM(100)
	CHARACTER*1 TDSSTYPE(100)
	CHARACTER*2 FFTTYPE
	DATA RE /6378./
	DATA AVRF,UPVSD/ 0., 0./
	DATA ISIZE /2048/
	DATA XNOSE,RFLANK/14.6,25.6/
	DATA RNOSE0,RFLANK0 /14.6,25.6/
	DATA EFFLEN /41.1, 3.79, 2.17, 41.1, 3.79, 2.17, 3*1./
	DATA NPLOTS /0/
C
	event = 'CDF'
	ITERM = -1
C	ITERM = -3
C	ITERM = 3
C
C	PLOT 4 HOURS CENTERED ON SCET
C
	NDAY = SCET
	print*,'at call splot, nday,scet of shock crossing=',nday,scet
	SCETINT = DINT(SCET)
	SCET1 = SCET - 2.D00/24.D00
	DNDAY = DBLE(NDAY)
	SCET1 = DMAX1(DNDAY,SCET1)
	SCET2 = SCET1 + 4.D00/24.D00
	SCET2 = DMIN1(SCET2,DNDAY + .999D00)
	PRINT*,'SPLOT CALLED, SCET1,SCET,SCET2=',SCET1,SCET,SCET2
	PRINT*,'TODAY in SPLOT',TODAY
	PRINT*,'NDAY,SCET1,2',NDAY,SCET1,SCET2
C
	CALL MGOINIT
	CALL MGOSETUP(ITERM)
	CALL MGOERASE
	NPLOTS = NPLOTS+1
C
C	PLOT MAGNITUDE OF B
C
C	NUMBER OF MFI POINTS IN FOUR HOURS
	NGET = 240./.75
	YMAX = 5.
C	SET POSITION TO BEGINNING OF DAY
	scett = 0.d00
	call w_channel_position(CDFCH,scett)
	call w_channel_position(CDFCH,scet1)
	NP = 0
C
 	  ok = w_event(CDFCH,EVENT)
	  IF(OK.EQ.1) THEN
	    ITEM = 'WIND_MFI_SCET_R8'
	    ok = W_ITEM_R8(cdfch, item, SCSV, NGET, ret_size)
	    ITEM = 'WIND_MFI_BMAG_R4'
	    ok = W_ITEM_R4(cdfch, item, YY, NGET, ret_size)
C	    print*,'load',n,np,ret_size,scsv(np),yy(np)
	  ELSE
	    IF(OK.EQ.82) GO TO 100
	  ENDIF
C	  DO NP = 1,RET_SIZE
	    NP = 0
	  DOWHILE(SCSV(NP).LT.SCET2)
	    NP = NP+1
	    PP(NP) = (SCSV(NP) - NDAY)*24.D00
	    YMAX = AMAX1(YY(NP),YMAX)
	ENDDO
	NP = MIN0(NP,NGET,RET_SIZE)
	NP = NP-1
C
 100 	  print*,'asked for and got, for MFI',nget,NP
	PRINT*,'1ST,LAST TIME',PP(1),PP(NP)
	PRINT*,'1ST,LAST BMAG',YY(1),YY(NP)
	CALL MGOWINDOW(1,6,1)
C	write(67,*) 'window1',gx1,gx2,gy1,gy2
	XSTRT = PP(1) - .1
	XEND = (SCET2-NDAY)*24.D00 + .1
	CALL MGOSETLIM(XSTRT,0.,XEND,1.1*YMAX)
	IF(NP.GT.0) CALL MGOCONNECT(PP,YY,NP)
	PPSH = (SCET - NDAY)*24.D00
	CALL MGOSETLTYPE(1)
	CALL MGORELOCATE(PPSH,0.)
	CALL MGODRAW(PPSH,1.1*YMAX)
	CALL MGOSETLTYPE(0)
C	CALL MGOSETEXPAND(.8)
	CALL MGOBOX(1,2)
	WRITE(STR,1001) TODAY(11:14),TODAY(15:16),TODAY(17:18)
 1001	FORMAT('HOURS OF ',A4,'/',A2,'/',A2)
	CALL MGOXLABEL(19,STR)
	CALL MGOYLABEL(3,'|B|')
C
C	CALL MGOSETEXPAND(1.)
C	TRANGE = GY2-GY1
C	TINC = .07*TRANGE
C	XTITLE = GX2 +.005*(GX2-GX1)
C	YTITLE = GY2
C	CALL MGOSETEXPAND(.8)
	PRINT*,'TODAY in SPLOT 2',TODAY
C
C	PLOT 3DP,SWE DENSITY,SOLARWIND SPEED
C		SWE IS DOTTED, 3DP IS SOLID
	scett = 0.d00
	call w_channel_position(CDFCH,scett)
	call w_channel_position(CDFCH,scet1)
	print*,'3DP channel set to',cdfch,scet1
C	NUMBER OF 3DP POINTS IN FOUR HOURS
	NGET = 240./.75
 	ok = w_event(CDFCH,EVENT)
        ITEM = 'WIND_SWE_SCET_R8'
	ok = W_ITEM_R8(cdfch, item, SCSV, NGET, ret_size)
 	print*,'asked for SCET8 and got',nget,RET_SIZE
	ITEM = 'WIND_SWE_DENSITY_R4'
	ok = W_ITEM_R4(cdfch, item, YY, NGET, ret_size)
	ITEM = 'WIND_SWE_VX(GSE)_R4'
	ok = W_ITEM_R4(cdfch, item, YYSS, NGET, ret_size)
C	PRINT*,'TODAY in SPLOT 3',TODAY

	YMAX = 5.
C	YSMAX = 1100.
	YSMAX = -10000.
	YSMIN = 10000.
	NP = 0
	N = 1
C	DO N = 1,NGET
	DOWHILE(SCSV(N).LT.SCET2.AND.N.LE.RET_SIZE)
	  IF(YY(N).GT.0.) THEN
	    NP = NP+1
	    PP(NP) = (SCSV(N) - NDAY)*24.D00
	    YY(NP) = YY(N)
	    YYSS(NP) = ABS(YYSS(N))
	    YMAX = AMAX1(YY(NP),YMAX)
	    YSMAX = AMAX1(YYSS(NP),YSMAX)
	    YSMIN = AMIN1(YYSS(NP),YSMIN)
	  ENDIF
	  N = N+1
	ENDDO
 	print*,'asked for SWE_DENS and got',nget,RET_SIZE,NP
	NPS = NP
	IF(NPS.GT.0) THEN
	  PRINT*,'SWE MIN,MAX',YSMIN,YSMAX
	  PRINT*,'1ST,LAST TIME',NPS,PP(1),PP(NPS)
	  PRINT*,'1ST,LAST SWE VX',NPS,YYSS(1),YYSS(NPS)
	ELSE
	  PRINT*,'NO SWE DATA'
	ENDIF
C
        ITEM = 'WIND_3DP_SCET_R8'
	ok = W_ITEM_R8(cdfch, item, SCSV, NGET, ret_size)
	ITEM = 'WIND_3DP_ION_DENSITY_R4'
	ok = W_ITEM_R4(cdfch, item, YY3D, NGET, ret_size)
	ITEM = 'WIND_3DP_ION_VX(GSE)_R4'
	ok = W_ITEM_R4(cdfch, item, YY3S, NGET, ret_size)
C
c***********
c	IF(ITERM.LT.0) THEN
c	  CALL MGOPRNTPLOT(NVEC)
c	  PRINT*,' NO. VECTORS PLOTTED',NVEC
c	ELSE
c	  CALL MGOTCLOSE
c	ENDIF
c	print*,'nget=',nget
c	print*,'return from splot'
c	if(1) return
c***********
	NP = 0
	N = 1
C	DO N = 1,NGET
	DOWHILE(SCSV(N).LT.SCET2.AND.N.LE.RET_SIZE)
	  IF(YY3D(N).GT.0..AND.YY3S(N).GT.-1.E8) THEN
	    NP = NP+1
	    PP1(NP) = (SCSV(N) - NDAY)*24.D00
	    SCSV(NP) = SCSV(N)
	    YY3D(NP) = YY3D(N)
	    YY3S(NP) = ABS(YY3S(N))
	    YMAX  = AMAX1(YY3D(NP),YMAX)
	    YSMAX = AMAX1(YY3S(NP),YSMAX)
	    YSMIN = AMIN1(YY3S(NP),YSMIN)
	  ENDIF
	  N = N+1
	ENDDO
	NP3 = NP-1
	IF(NP3.GT.0) THEN
	  PRINT*,'1ST,LAST TIME',NP3,PP1(1),PP1(NP3)
	  PRINT*,'1ST,LAST 3DP',NP3,YY3D(1),YY3D(NP3)
	ELSE
	  PRINT*,'NO 3DP DATA'
	ENDIF
C
	PRINT*,'TODAY in SPLOT 4',TODAY
	CALL MGOWINDOW(1,6,2)
C	write(67,*) 'window2',gx1,gx2,gy1,gy2
	CALL MGOSETLIM(XSTRT,0.,XEND,YMAX)
	CALL MGOSETLTYPE(1)
	IF(NPS.GT.0) CALL MGOCONNECT(PP,YY,NPS)		! SWE ION DENSITY
	CALL MGOSETLTYPE(0)
	IF(NP3.GT.0) CALL MGOCONNECT(PP1,YY3D,NP3)	! 3DP ION DENSITY
C	CALL MGOSETEXPAND(.8)
	CALL MGOBOX(1,2)
	CALL MGOYLABEL(7,'DENSITY')
C
C	PLOT SOLAR WIND SPEED
C
	CALL MGOWINDOW(1,6,3)
C	write(67,*) 'window3',gx1,gx2,gy1,gy2
	CALL MGOSETLIM(XSTRT,YSMIN,XEND,YSMAX)
	CALL MGOSETLTYPE(1)
	IF(NPS.GT.0) CALL MGOCONNECT(PP,YYSS,NPS)		! SWE VX
	CALL MGOSETLTYPE(0)
	IF(NP3.GT.0) CALL MGOCONNECT(PP1,YY3S,NP3)		! 3DP VX
	CALL MGOBOX(1,2)
	CALL MGOYLABEL(6,'ION VX')
C
	PPSH = (SCET - NDAY)*24.D00
	CALL MGOSETLTYPE(1)
	CALL MGORELOCATE(PPSH,0.)
	CALL MGODRAW(PPSH,1.1*YMAX)
	CALL MGOSETLTYPE(0)
C
C	PLOT B dot N
C
C
	scett = 0.d00
	call w_channel_position(CDFCH,scett)
 	ok = w_event(CDFCH,event)
C	print*,'channel set to 0.',cdfch,scett
	call w_channel_position(CDFCH,scet1)
c 	ok = w_event(CDFCH,event)
	SSMAX = 0.
	SSMIN = 100.
	BDMAX = 90.
	BDMIN = 0.
	print*,'going to 1,np3 loop',np3
	    PPBN1 = (SCET-NDAY)*24.D00 - 2.D00/60.D00
	    PPBN2 = PPBN1 - MINUP/60.D00
	DO N = 1,NP3
  	  ok = w_event(CDFCH,event)
	  SCETT = SCSV(N)
	  ITEM = 'WIND_MFI_BX(GSE)_R4'
	  ok = W_ITEM_R4(cdfch, item, BX, 1, ret_size)
	  ITEM = 'WIND_MFI_BY(GSE)_R4'
	  ok = W_ITEM_R4(cdfch, item, BY, 1, ret_size)
	  ITEM = 'WIND_MFI_BZ(GSE)_R4'
	  ok = W_ITEM_R4(cdfch, item, BZ, 1, ret_size)
	  BNORM = SQRT(BX**2+BY**2+BZ**2)
	  SNORM = SQRT(SHNORM(1)**2 + SHNORM(2)**2 + SHNORM(3)**2)
	  BDOTN(N) = 0.
	  IF(BNORM*SNORM.NE.0.) THEN
	  print*,'bnorm check',snorm,bnorm,bnorm*snorm
     	   BDOTN(N)=(BX*SHNORM(1)+BY*SHNORM(2)+BZ*SHNORM(3))/BNORM/SNORM
	   BDOTN(N) = ACOSD(ABS(BDOTN(N)))
	   BDMAX = AMAX1(BDMAX,BDOTN(N))
	   BDMIN = AMIN1(BDMIN,BDOTN(N))
	  ENDIF
	  time = (scett-nday)*24.d00
	ENDDO
C
	print*,'end np3 loop'
	CALL MGOWINDOW(1,6,4)
C	write(67,*) 'window4',gx1,gx2,gy1,gy2
	CALL MGOSETLIM(XSTRT,BDMIN,XEND,BDMAX)
	IF(NP3.GT.0) CALL MGOCONNECT(PP1,BDOTN,NP3)
C	  DRAW AVERAGE BN ANGLE
	IF(ANGBN.GT.0.) THEN
	  CALL MGOSETLTYPE(1)
	  CALL MGORELOCATE(PPBN1,ANGBN)
	  CALL MGODRAW(PPBN2,ANGBN)
	  CALL MGOSETLTYPE(0)
	ENDIF
	CALL MGOSETLTYPE(1)
	CALL MGORELOCATE(PPSH,0.)
	CALL MGODRAW(PPSH,1.1*YMAX)
	CALL MGOSETLTYPE(0)
	CALL MGOBOX(1,2)
	CALL MGOYLABEL(10,'ANG(B\u.N)')
C
C	PLOT TDSF EVENTS
C
	TDS_COUNT = 0
C	  COUNT EVENTS WITHIN 5 MIN OF CROSSING TIME
	SCETTDS1 = SCET - 5.D00/1440.D00
	SCETTDS2 = SCET + 5.D00/1440.D00
	EVENT = 'TDSF'
	YMAX = 130.
	CALL MGOWINDOW(1,6,5)
C	write(67,*) 'window5',gx1,gx2,gy1,gy2
	CALL MGOSETLIM(XSTRT,0.,XEND,YMAX)
	CALL MGOBOX(1,2)
	CALL MGOYLABEL(4,'TDSF')
C
	tdsch = cdfch
	EVENT = 'TDSF'
	scett = 0.d00
	print*,'going to get TDS events'
	call w_channel_position(TDSCH,scett)
	call w_channel_position(TDSCH,scet1)
	print*,'TDS channel set to',TDSCH,scet1
	PWR = 0.
	TDSM = 0.
	NSLOW = 0
 200	OK = W_EVENT(TDSCH,EVENT)
C	print*,'tds event, event, ok=  ',event,ok
	IF(OK.EQ.82) GO TO 250
	IF(OK.NE.1) GO TO 200
C	ITEM = 'DATA'
C	ok = W_ITEM_I4(TDSCH, item, NDATA, ISIZE, ret_size)
	ITEM = 'EVENT_NUMBER'
	ok = W_ITEM_I4(TDSCH, item, NEVT, 1, ret_size)
C	print*,'tds event, event, ok, no =  ',event,ok,nevt
	IPRC = 4
	CALL TDS_PHYS(TDSCH,IPRC,NDATA,VDATA,SPECT)
	ITEM = 'EVENT_SCET'
	ok = W_ITEM_I4(TDSCH, item, SCETI4, 2, ret_size)
	ITEM = 'FAST_RX_SPEED'
	ok = W_ITEM_I4(cdfch, item, ISPS, 1, ret_size)
C	ITEM = 'EVENT_TRUE_SCET_R8'
C	ITEM = 'EVENT_SCET_R8'
C	ok = W_ITEM_R8(TDSCH, item, SCETT, 1, ret_size)
C
	ITEM = 'CHANNEL'
	ok = W_ITEM_I4(TDSCH, item, TDS_CHANNEL, 1, ret_size)
	EFL = EFFLEN(TDS_CHANNEL)
C
	ITEM = 'EVENT_NUMBER'
	ok = W_ITEM_I4(TDSCH, item, NEVT, 1, ret_size)
C	IF(SCETT.GE.SCETTDS1.AND.SCETT.LE.SCETTDS2) print*,'tds',sceti4,nevt
	ITEM = 'FAST_RX_SPEED_R4'
	ok = W_ITEM_R4(TDSCH, item, SPS, 1, ret_size)
	ITEM = 'FAST_SAMPLER_THRESHOLD'
	ok = W_ITEM_R4(TDSCH, item, VMIN, 1, ret_size)
	NVMAX = MAX0(NDATA(1023),NDATA(1024),NDATA(1025))
	NVMIN = MIN0(NDATA(1023),NDATA(1024),NDATA(1025))
	VMAX = MAX0(NVMAX-128,128-NVMIN)
C	ITEM = 'EVENT_CENTER_SCET_R8'
	ITEM = 'EVENT_SCET_R8'
	ok = W_ITEM_R8(TDSCH, item, SCETT, 1, ret_size)
	IF(SCETT.GE.SCET1.AND.SCETT.LE.SCET2) THEN
	  TDSTIME = (SCETT - NDAY)*24.D00
	  CALL MGORELOCATE(TDSTIME,VMIN)
	  CALL MGODRAW(TDSTIME,VMAX)
C	  CALL TDSDATA(TDSCH,NDATA,SPS,AVRF,UPVSD)
	  IF(TDS_CHANNEL.EQ.1) THEN
	    AVRF1 = AVRF
	    UPVSD1 = UPVSD
	  ELSE
	    AVRF2 = AVRF
	    UPVSD2 = UPVSD
	  ENDIF
C
	  IF(TDS_CHANNEL.EQ.1) WRITE(64,1064) NEVT,SCETI4,'T',TDS_CHANNEL,VMAX
     1		,AVRF1,UPVSD1,UPVSD2
 1064	  FORMAT(I10,I10,I8,A2,I3,F5.0,F7.3,2F6.2)
	ENDIF
C
	IF(SCETT.GE.SCETTDS1.AND.SCETT.LE.SCETTDS2) THEN
	   TDS_COUNT = TDS_COUNT+1
	   TDSC = TDS_COUNT
	  DO N = 1,2048
	    TDSM = AMAX1(TDSM,1000.*ABS(VDATA(N))/EFL)  ! in mV/m
	    PWR = PWR + (VDATA(N)/EFL)**2
	  ENDDO
	ENDIF
C
	GO TO 200
C
 250	EVENT = 'FILL'
	scett = 0.d00
	call w_channel_position(TDSCH,scett)
	call w_channel_position(TDSCH,scet1)
	print*,'FILL channel set to',TDSCH,scet1
 240	CONTINUE
 	OK = W_EVENT(TDSCH,EVENT)
C	print*,'fill event, event, ok=  ',event,ok
	IF(OK.EQ.82) GO TO 400
	IF(OK.NE.1) GO TO 240
C	ITEM = 'DATA'
C	ok = W_ITEM_I4(TDSCH, item, NDATA, ISIZE, ret_size)
	ITEM = 'CHANNEL'
	ok = W_ITEM_I4(TDSCH, item, TDS_CHANNEL, 1, ret_size)
	EFL = EFFLEN(TDS_CHANNEL)
	IF(TDS_CHANNEL.GT.2) GO TO 240
	IPRC = 4
	CALL TDS_PHYS(TDSCH,IPRC,NDATA,VDATA,SPECT)
	ITEM = 'EVENT_SCET'
	ok = W_ITEM_I4(TDSCH, item, SCETI4, 2, ret_size)
C	print*,'item,val,ok ',item,sceti4,ok
	ITEM = 'EVENT_SCET_R8'
	ok = W_ITEM_R8(TDSCH, item, SCETT, 1, ret_size)
C	FILL CHANNELS ARE IN CHRONOLOGICAL ORDER
	IF(SCETT.GT.SCET2+1.D-02) GO TO 400
c	print*,'item,val,ok ',item,scett,ok
	ITEM = 'EVENT_NUMBER'
	ok = W_ITEM_I4(TDSCH, item, NEVT, 1, ret_size)
	ITEM = 'FAST_RX_SPEED_R4'
	IF(TDS_CHANNEL.GT.2) ITEM = 'SLOW_RX_SPEED_R4'
C	ITEM = 'EVENT_TRUE_SCET_R8'
	ok = W_ITEM_R4(TDSCH, item, SPS, 1, ret_size)
	ITEM = 'FAST_SAMPLER_THRESHOLD'
	IF(TDS_CHANNEL.GT.2) ITEM = 'SLOW_SAMPLER_THRESHOLD'
	ok = W_ITEM_R4(TDSCH, item, VMIN, 1, ret_size)
	NVMAX = MAX0(NDATA(1023),NDATA(1024),NDATA(1025))
	NVMIN = MIN0(NDATA(1023),NDATA(1024),NDATA(1025))
	VMAX = MAX0(NVMAX-128,128-NVMIN)
C	ITEM = 'EVENT_CENTER_SCET_R8'
C	print*,'FILL EVENT AT',TDSCH,scett,nevt
	IF(SCETT.GE.SCET1.AND.SCETT.LE.SCET2) THEN
	  TDSTIME = (SCETT - NDAY)*24.D00
	  CALL MGORELOCATE(TDSTIME,VMIN)
	  CALL MGODRAW(TDSTIME,VMAX)
	  CALL TDSDATA(TDSCH,NDATA,SPS,AVRF,UPVSD)
	  IF(TDS_CHANNEL.EQ.1.OR.TDS_CHANNEL.EQ.5) THEN
	    AVRF1 = AVRF
	    UPVSD1 = UPVSD
	  ELSE
	    AVRF2 = AVRF
	    UPVSD2 = UPVSD
	  ENDIF
	  IF(TDS_CHANNEL.GT.2) THEN
	    NSLOW = NSLOW + 1
	    TDSSTIM(NSLOW) = TDSTIME
	    ITEM = 'SOURCE'
	    ok = W_ITEM_I4(TDSCH, item, ISRC, 1, ret_size)
	    IF(ISRC.LE.6) TDSSTYPE(NSLOW) = 'E'
	    IF(ISRC.GT.6) TDSSTYPE(NSLOW) = 'B'
	  ENDIF
C
	IF(SCETT.GE.SCETTDS1.AND.SCETT.LE.SCETTDS2) THEN
	  TDS_COUNT = TDS_COUNT+1
	  TDSC = TDS_COUNT
	  DO N = 1,2048
	    TDSM = AMAX1(TDSM,1000.*ABS(VDATA(N))/EFL)  ! in mV/m
	    PWR = PWR + (VDATA(N)/EFL)**2
	  ENDDO
	ENDIF
C
	  IF(TDS_CHANNEL.EQ.1) WRITE(64,1064) NEVT,SCETI4,'F',TDS_CHANNEL,VMAX
     1		,AVRF,UPVSD1,UPVSD2
	  IF(TDS_CHANNEL.EQ.3) WRITE(64,1064) NEVT,SCETI4,'F',TDS_CHANNEL,VMAX
     1		,AVRF,UPVSD1,UPVSD2
	ENDIF
C
	IF(TDS_CHANNEL.GT.2) GO TO 240
C	  IF(SCETT.GE.SCETTDS1.AND.SCETT.LE.SCETTDS2) TDS_COUNT = TDS_COUNT+1
	GO TO 240			! GET ANOTHER FILL EVENT
C
 400	EVENT = 'TDSS'
	scett = 0.d00
	call w_channel_position(TDSCH,scett)
	call w_channel_position(TDSCH,scet1)
	print*,'FILL channel set to',TDSCH,scet1
C
C	PUT EXTRA STUFF IN WINDOW 6
C
	CALL MGOWINDOW(1,6,6)
	YMAX = 130.
	YMIN = 10.
	CALL MGOSETLIM(XSTRT,0.,XEND,YMAX)
	CALL MGOBOX(1,2)
	CALL MGOYLABEL(4,'TDSS')
	PRINT*,'NSLOW = ',NSLOW
	DO N5 = 1,NSLOW
	  CALL MGORELOCATE(TDSSTIM(N5),YMIN)
	  CALL MGOLABEL(1,TDSSTYPE(N5))
	ENDDO
C
 440	CONTINUE
 	OK = W_EVENT(TDSCH,EVENT)
C	print*,'fill event, event, ok=  ',event,ok
	IF(OK.EQ.82) GO TO 500
	IF(OK.NE.1) GO TO 440
C	ITEM = 'DATA'
C	ok = W_ITEM_I4(TDSCH, item, NDATA, ISIZE, ret_size)
	ITEM = 'CHANNEL'
	ok = W_ITEM_I4(TDSCH, item, TDS_CHANNEL, 1, ret_size)
	EFL = EFFLEN(TDS_CHANNEL)
	IPRC = 4
	CALL TDS_PHYS(TDSCH,IPRC,NDATA,VDATA,SPECT)
	ITEM = 'EVENT_SCET'
	ok = W_ITEM_I4(TDSCH, item, SCETI4, 2, ret_size)
C	print*,'item,val,ok ',item,sceti4,ok
	ITEM = 'EVENT_SCET_R8'
	ok = W_ITEM_R8(TDSCH, item, SCETT, 1, ret_size)
C	FILL CHANNELS ARE IN CHRONOLOGICAL ORDER, BUT TDSS ARE NOT
C	IF(SCETT.GT.SCET2+1.D-02) GO TO 300
c	print*,'item,val,ok ',item,scett,ok
	ITEM = 'EVENT_NUMBER'
	ok = W_ITEM_I4(TDSCH, item, NEVT, 1, ret_size)
C	ITEM = 'FAST_RX_SPEED_R4'
C	IF(TDS_CHANNEL.GT.2) ITEM = 'SLOW_RX_SPEED_R4'
C	ok = W_ITEM_R4(TDSCH, item, SPS, 1, ret_size)
	ITEM = 'FAST_SAMPLER_THRESHOLD'
	IF(TDS_CHANNEL.GT.2) ITEM = 'SLOW_SAMPLER_THRESHOLD'
	ok = W_ITEM_R4(TDSCH, item, VMIN, 1, ret_size)
	NVMAX = MAX0(NDATA(1023),NDATA(1024),NDATA(1025))
	NVMIN = MIN0(NDATA(1023),NDATA(1024),NDATA(1025))
	VMAX = MAX0(NVMAX-128,128-NVMIN)
C	ITEM = 'EVENT_CENTER_SCET_R8'
C	print*,'FILL EVENT AT',TDSCH,scett,nevt
	IF(SCETT.GE.SCET1.AND.SCETT.LE.SCET2) THEN
	  TDSTIME = (SCETT - NDAY)*24.D00
	  CALL MGORELOCATE(TDSTIME,VMIN)
	  CALL MGODRAW(TDSTIME,VMAX)
C
C	IF(SCETT.GE.SCETTDS1.AND.SCETT.LE.SCETTDS2) THEN
C	  TDS_COUNT = TDS_COUNT+1
C	  TDSC = TDS_COUNT
C	  DO N = 1,2048
C	    TDSM = AMAX1(TDSM,1000.*ABS(VDATA(N))/EFL)  ! in mV/m
C	    PWR = PWR + (VDATA(N)/EFL)**2
C	  ENDDO
C	ENDIF
C
C	  IF(TDS_CHANNEL.EQ.1) WRITE(64,1064) NEVT,SCETI4,'F',TDS_CHANNEL,VMAX
C     1		,AVRF,UPVSD1,UPVSD2
C	  IF(TDS_CHANNEL.EQ.3) WRITE(64,1064) NEVT,SCETI4,'F',TDS_CHANNEL,VMAX
C     1		,AVRF,UPVSD1,UPVSD2
	ENDIF
C
C	IF(TDS_CHANNEL.GT.2) GO TO 240
C	  IF(SCETT.GE.SCETTDS1.AND.SCETT.LE.SCETTDS2) TDS_COUNT = TDS_COUNT+1
	GO TO 440			! GET ANOTHER TDSS EVENT
C
 500	CONTINUE
	print*,'going to get FFTs'
	CALL MGOWINDOW(1,6,6)
	YMAX = 130.
	YMIN = 10.
	CALL MGOSETLIM(XSTRT,0.,XEND,YMAX)
	DO NCHGP = 1,3
	  nplot = 0
	  scett = 0.d00
	  call w_channel_position(TDSCH,scett)
	  call w_channel_position(TDSCH,scet1)
	  print*,'FFT channel set to',TDSCH,scet1
	  if(nchgp.eq.1) event = 'FFTH'
	  if(nchgp.eq.2) event = 'FFTM'
	  if(nchgp.eq.3) event = 'FFTL'
	  first_channel = 1
	  if(nchgp.eq.2) first_channel = 3
	  if(nchgp.eq.3) first_channel = 7
	  nw = 2
	  if(nchgp.gt.1) nw = 4
 540	  ok = w_event(tdsch,event)
C
C	print*,'fft event nchgp,ok =',nchgp,ok
	  IF(OK.EQ.82) GO TO 550
	  IF(OK.NE.1) GO TO 540
C
c	  IF(OK.EQ.1) THEN
c	    ITEM = 'WIND_MFI_SCET_R8'
c	    ok = W_ITEM_R8(cdfch, item, SCSV, NGET, ret_size)
c	    ITEM = 'WIND_MFI_BMAG_R4'
c	    ok = W_ITEM_R4(cdfch, item, YY, NGET, ret_size)
C	    print*,'load',n,np,ret_size,scsv(np),yy(np)
c	  ELSE
c	    IF(OK.EQ.82) GO TO 550
c	  ENDIF
c
	  item = 'EVENT_SCET_R8'
	  ok = w_item_r8(tdsch, item, scett, 1, return_size)
	  call w_ur8_to_ydoy(scet8,yyyy,idoy,msec)
C
	  item = 'SOURCE'
	  ok = w_item_i4(tdsch, item, ISRC, 1, return_size)
	  if(isrc.gt.6) Bcount = Bcount + 1
C
	  item = 'channel_number'
	  ok = w_item_i4(tdsch, item, fft_channel, 1, return_size)
	  if(fft_channel.ne.first_channel) then
		last_channel = fft_channel
C		type*,'channel',fft_channel,' is not first, restart'
		go to 540
	  else
		scetfft = scett
	   endif
C
	  EXPSAVE = EXPAND
	  CALL MGOSETEXPAND(.3)
	  FFTTYPE = 'EH'
	  IF(BCOUNT.GE.2.AND.NCHGP.EQ.2) FFTTYPE = 'BM'
	  IF(BCOUNT.LT.2.AND.NCHGP.EQ.2) FFTTYPE = 'EM'
	  IF(BCOUNT.GE.2.AND.NCHGP.EQ.3) FFTTYPE = 'BL'
	  IF(BCOUNT.LT.2.AND.NCHGP.EQ.3) FFTTYPE = 'EL'
C
	  PPFFT = (SCETFFT - NDAY)*24.D00
	  CHGP = .25*YMAX*(4-NCHGP)
	  CALL MGORELOCATE(PPFFT,CHGP)
	  CALL MGOLABEL(2,FFTTYPE)
	  CALL MGOSETEXPAND(EXPSAVE)
C
	  Bcount = 0
	  if(scett.le.scet2) go to 540
 550	  continue
	ENDDO

 900	CONTINUE
C
	CLOSE(UNIT=64)
C
 301	CALL MGOPLOTID('PROCIPSHK','SPLOT')
C
	call w_channel_position(CDFCH,scet1)
 	ok = w_event(CDFCH,event)
	ITEM = 'WIND_ORBIT_X(GSE)_R8'
	ok = W_ITEM_R8(cdfch, item, RGSE, 1, ret_size)
	XGSE1 = RGSE/RE
	ITEM = 'WIND_ORBIT_Y(GSE)_R8'
	ok = W_ITEM_R8(cdfch, item, RGSE, 1, ret_size)
	YGSE1 = RGSE/RE
	ITEM = 'WIND_ORBIT_Z(GSE)_R8'
	ok = W_ITEM_R8(cdfch, item, RGSE, 1, ret_size)
	ZGSE1 = RGSE/RE
	RRE = SQRT(XGSE1**2+YGSE1**2+ZGSE1**2)
	WRITE(STR,1002) XGSE1,YGSE1,ZGSE1,RRE
 1002	FORMAT('X,Y,Z,R=',4F6.1)
	CALL MGORELOCATE(XSTRT,1.2*YMAX)
	CALL MGOPUTLABEL(36,STR,9)
C
	call w_channel_position(CDFCH,scet2)
 	ok = w_event(CDFCH,event)
	ITEM = 'WIND_ORBIT_X(GSE)_R8'
	ok = W_ITEM_R8(cdfch, item, RGSE, 1, ret_size)
	XGSE2 = RGSE/RE
	ITEM = 'WIND_ORBIT_Y(GSE)_R8'
	ok = W_ITEM_R8(cdfch, item, RGSE, 1, ret_size)
	YGSE2 = RGSE/RE
	ITEM = 'WIND_ORBIT_Z(GSE)_R8'
	ok = W_ITEM_R8(cdfch, item, RGSE, 1, ret_size)
	ZGSE2 = RGSE/RE
	RRE = SQRT(XGSE2**2+YGSE2**2+ZGSE2**2)
	WRITE(STR,1002) XGSE2,YGSE2,ZGSE2,RRE
	CALL MGORELOCATE(XEND,1.2*YMAX)
	CALL MGOPUTLABEL(36,STR,7)
	IF(ITERM.LT.0) THEN
	  CALL MGOPRNTPLOT(NVEC)
	  PRINT*,' NO. VECTORS PLOTTED',NVEC
C	ELSE
c	  CALL MGOTCLOSE
	ENDIF
C
	TDSC = TDS_COUNT
	IF(TDS_COUNT.NE.0) PWR =  2.*PWR/TDS_COUNT   ! 2 CHANNELS, EVENTS
C		APPEAR TWICE
	RETURN
	END
	SUBROUTINE TDSDATA(CH,NDATA,SPS,AVRF,UPVSD)
C
C	CALCULATES AVERAGE FREQUENCY AND PERHAPS OTHER THINGS.
C
	INTEGER*4 CH,NDATA(1)
	REAL ZCROSS(1030),DATA(2050),SPECT(1025)
	REAL SPS,AVRF,UPVSDR

C
C	FIND ZERO CROSSINGS, USING ONLY THE CENTER SECTION FROM 512 TO 1536
C
	IZCNT = 0
	IL = 511
	IZ = IL
	IPROCESS = 4
	CALL TDS_PHYS(CH,IPROCESS,NDATA,DATA,SPECT)
C	  IF(NDATA(IL).EQ.0) PRINT*,'ZERO DATA',IL,NDATA(IL),NDATA(IL+1)
	IPROCESS = 0
	DO IL = 512,1536
	  IZ = IL
C	  IF(NDATA(IL).EQ.128) PRINT*,'ZERO DATA',IL,NDATA(IL-1),
C     1   NDATA(IL),NDATA(IL+1)
C
C		COUNT ONLY POS TO NEG
	  IF(NDATA(IL).GT.128.AND.NDATA(IL+1).LE.128) THEN
	        IZCNT = IZCNT+1
		IF(IPROCESS.EQ.0) THEN
		  S1 = NDATA(IL)
		  S2 = NDATA(IL+1)
		ELSE
		  S1 = DATA(IL)
		  S2 = DATA(IL+1)
		ENDIF
	        DELZ = .5
		IF((S1-S2).NE.0.) DELZ = S1/(S1 - S2)
		IF(DELZ.LE.1.AND.DELZ.GE.0.) THEN
			ZCROSS(IZCNT) = IL + DELZ
		ELSE
			ZCROSS(IZCNT) = IL + .5
		ENDIF
	  ENDIF
	ENDDO
C	print*,'ndata',ndata(1023),ndata(1024),ndata(1025),ndata(1026)
C	print*,'data',data(1023),data(1024),data(1025),data(1026)
C	print*,'found zero crossings',izcnt
C
C	CALCULATE UP VS DOWN
C
	VUP = 0.
	VDWN = 0.
	DO N = 2,2047
	  VUP = AMAX1(VUP,DATA(N))
	  VDWN = AMIN1(VDWN,DATA(N))
	ENDDO
	UPVSD = 0.
	IF(VDWN.NE.0.) UPVSD = ABS(VUP)/ABS(VDWN)
	IF(UPVSD.LT.1..AND.UPVSD.GT.0.) UPVSD = 1./UPVSD
C
	SPSKHZ = .001*SPS
C	PRINT*,'IZCNT,SPS',IZCNT,SPS
	AVRPER = 1024.
	IF(IZCNT.GT.1) THEN
	  AVRPER = (ZCROSS(IZCNT)-ZCROSS(1))/(IZCNT-1)
C	  RMS = SQRT(STD/(IZCNT-1) - AVRPER**2)
C	  IF(IZCNT.NE.0) RMS = SPSKHZ*(1./AVRPER - 1./(AVRPER+RMS))
	ENDIF
	AVRF = 1./AVRPER
C
	RETURN
	END

ext: for


text:
module strings
  use iso_c_binding, only: &
       c_loc, &
       c_null_char, &
       c_ptr

  implicit none

  private

  public &
       string, &
       string_pool, &
       string_pool_allocate, &
       string_pool_finalize, &
       string_pool_initialize, &
       string_to_c_string, &
       string_to_c_string_with_pool

  type string
     character(:), allocatable :: value
  end type string

  type string_pointer
     type(string), pointer :: value
  end type string_pointer

  type string_pool
     integer :: offset, capacity
     type(string_pointer), allocatable :: pointers (:)
  end type string_pool

contains

  function string_pool_allocate(pool)
    type(string_pool), intent(inout) :: pool
    type(string), pointer :: string_pool_allocate

    type(string_pointer), allocatable :: buffer (:)

    allocate (string_pool_allocate)

    pool % pointers (pool % offset) % value => string_pool_allocate
    pool % offset = pool % offset + 1

    if (pool % offset > pool % capacity) then
       allocate (buffer (pool % capacity * 2), source=pool % pointers)
       call move_alloc(buffer, pool % pointers)
       pool % capacity = pool % capacity * 2
    end if
  end function string_pool_allocate

  subroutine string_pool_finalize(pool)
    type(string_pool), intent(inout) :: pool

    integer i

    do i = 1, pool % offset - 1
       deallocate (pool % pointers (i) % value)
    end do
  end subroutine string_pool_finalize

  subroutine string_pool_initialize(pool, capacity)
    type(string_pool), intent(inout) :: pool
    integer, intent(in) :: capacity

    pool % offset = 1
    pool % capacity = capacity
    allocate (pool % pointers (capacity))
  end subroutine string_pool_initialize

  function string_to_c_string(fstring, buffer)
    character(*), intent(in) :: fstring
    type(string), intent(inout), target :: buffer
    type(c_ptr) string_to_c_string

    buffer % value = fstring // c_null_char
    string_to_c_string = c_loc(buffer % value)
  end function string_to_c_string

  function string_to_c_string_with_pool(fstring, pool)
    character(*), intent(in) :: fstring
    type(string_pool), intent(inout) :: pool
    type(c_ptr) string_to_c_string_with_pool

    type(string), pointer :: buffer

    buffer => string_pool_allocate(pool)

    buffer % value = fstring // c_null_char
    string_to_c_string_with_pool = c_loc(buffer % value)
  end function string_to_c_string_with_pool

end module strings

max_stars_repo_name: dram/fortran-json-glib

max_issues_repo_name: dram/fortran-json-glib

max_forks_repo_name: dram/fortran-json-glib


text:
less_toxic(y1,u1).
less_toxic(y1,e1).
less_toxic(bb1,jj1).
less_toxic(b1,c1).
less_toxic(h1,c1).
less_toxic(ll1,v1).
less_toxic(k1,c1).
less_toxic(ff1,c1).
less_toxic(ff1,b1).
less_toxic(u1,c1).
less_toxic(kk1,i1).
less_toxic(r1,b1).
less_toxic(ee1,c1).
less_toxic(g1,e1).
less_toxic(bb1,f1).
less_toxic(kk1,c1).
less_toxic(ee1,t1).
less_toxic(r1,d1).
less_toxic(ii1,b1).
less_toxic(k1,u1).
less_toxic(o1,jj1).
less_toxic(w1,dd1).
less_toxic(ee1,cc1).
less_toxic(cc1,dd1).
less_toxic(g1,p1).
less_toxic(s1,jj1).
less_toxic(ll1,t1).
less_toxic(s1,dd1).
less_toxic(l1,c1).
less_toxic(bb1,cc1).
less_toxic(s1,c1).
less_toxic(y1,p1).
less_toxic(r1,aa1).
less_toxic(w1,c1).
less_toxic(l1,f1).
less_toxic(ll1,f1).
less_toxic(dd1,f1).
less_toxic(i1,c1).
less_toxic(j1,aa1).
less_toxic(n1,b1).
less_toxic(b1,f1).
less_toxic(f1,cc1).
less_toxic(v1,y1).
less_toxic(dd1,w1).
less_toxic(cc1,ii1).
less_toxic(aa1,cc1).
less_toxic(p1,u1).
less_toxic(e1,s1).
less_toxic(d1,dd1).
less_toxic(c1,i1).
less_toxic(a1,k1).
less_toxic(l1,ff1).
less_toxic(f1,ee1).

max_stars_repo_name: manoelfranca/cilppp

max_issues_repo_name: manoelfranca/cilppp

max_forks_repo_name: manoelfranca/cilppp


text:
!
! COMMON_r2p1_Module
!
! Module containing procedures for COMMON Record r2p1.
!
!
! CREATION HISTORY:
!       Written by:   Paul van Delst, 24-Dec-2012
!                     paul.vandelst@noaa.gov
!

MODULE COMMON_r2p1_Module

  ! -----------------
  ! Environment setup
  ! -----------------
  ! Module usage
  USE Type_Kinds     , ONLY: fp
  USE File_Utility   , ONLY: File_Open
  USE Message_Handler, ONLY: SUCCESS, FAILURE, INFORMATION, Display_Message
  ! Line-by-line model parameters
  USE LBL_Parameters
  ! Disable implicit typing
  IMPLICIT NONE

  ! ----------
  ! Visibility
  ! ----------
  ! Everything private by default
  PRIVATE
  ! Datatypes
  PUBLIC :: COMMON_r2p1_type
  ! Procedures
  PUBLIC :: COMMON_r2p1_Write

  ! -----------------
  ! Module parameters
  ! -----------------
  CHARACTER(*), PARAMETER :: MODULE_VERSION_ID = &
  ! Message string length
  INTEGER, PARAMETER :: ML = 256
  ! The record I/O format
  CHARACTER(*), PARAMETER :: COMMON_R2P1_FMT = '(1x,i1,i3,i5,f10.2,20x,f8.2,4x,f8.2,5x,f8.3)'

  ! -------------
  ! Derived types
  ! -------------
  TYPE :: COMMON_r2p1_type
    INTEGER  :: iform  = 0       !  Column amount format flag. If == 1, read data es15.7 format
    INTEGER  :: nlayrs = 0       !  Number of layers (maximum of 200)
    INTEGER  :: nmol   = 7       !  Value of highest molecule number used (default = 7; maximum of 35)
    REAL(fp) :: secnt0 = 1.0_fp  !  Column amount scale factor. If +ve, looking up; if -ve looking down
    REAL(fp) :: zh1    = 0.0_fp  !  Observer altitude
    REAL(fp) :: zh2    = 0.0_fp  !  End point altitude
    REAL(fp) :: zangle = 0.0_fp  !  Zenith angle for path calculation
  END TYPE COMMON_r2p1_type

CONTAINS

  FUNCTION COMMON_r2p1_Write(r2p1,fid) RESULT(err_stat)

    ! Arguments
    TYPE(COMMON_r2p1_type), INTENT(IN) :: r2p1
    INTEGER               , INTENT(IN) :: fid
    ! Function result
    INTEGER :: err_stat
    ! Function parameters
    CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'COMMON_r2p1_Write'
    ! Function variables
    CHARACTER(ML) :: msg
    CHARACTER(ML) :: io_msg
    INTEGER :: io_stat

    ! Setup
    err_stat = SUCCESS
    ! ...Check unit is open
    IF ( .NOT. File_Open(fid) ) THEN
      msg = 'File unit is not connected'
      CALL Cleanup(); RETURN
    END IF

    ! Write the record
    WRITE( fid,FMT=COMMON_R2P1_FMT,IOSTAT=io_stat,IOMSG=io_msg) r2p1
    IF ( io_stat /= 0 ) THEN
      msg = 'Error writing record - '//TRIM(io_msg)
      CALL Cleanup(); RETURN
    END IF

  CONTAINS

    SUBROUTINE CleanUp()
      err_stat = FAILURE
      CALL Display_Message( ROUTINE_NAME,msg,err_stat )
    END SUBROUTINE CleanUp

  END FUNCTION COMMON_r2p1_Write

END MODULE COMMON_r2p1_Module

max_stars_repo_name: hsbadr/crtm

max_issues_repo_name: hsbadr/crtm

max_forks_repo_name: hsbadr/crtm


text:
      SUBROUTINE MBMODE(INPUT,OUT,ICOR,NCOR,Z,NI,ND,XD,YD,IS,CR)
C
C     MBMODE BUILDS THE MODE LIKE DATA ON OUT FROM SURFACE SPLINE INTER
C
      DIMENSION Z(1),XD(1),YD(1),NAME(2)
      DATA NAME /4HMBMO,4HDE  /
      NNI = NI*2
      NND = ND*2
      IF(ICOR+NNI+NND.GT.NCOR) CALL MESAGE(-8,0,NAME)
      CALL FREAD(INPUT,Z(ICOR),NNI,0)
      IDP = ICOR + NNI
      L = 0
      DO 10 I=1,ND
      Z(IDP+L) = XD(I)
      Z(IDP+L+1) = YD(I)
      L = L+2
   10 CONTINUE
      ICC = IDP+L
      NCORE = NCOR-ICC
C
C     CALL SSPLIN TO INTERPOLATE
C
      CALL SSPLIN(NI,Z(ICOR),ND,Z(IDP),0,0,1,0,0.0,Z(ICC),NCORE,IS)
      IF(IS.EQ.2) GO TO 1000
C
C     REORDER INTO MACH BOX ORDER
C
      M = IDP+ND
      ICC = ICC-1
      DO 30 I=1,NI
      L = 0
      DO 20 J=1,NND,2
      Z(IDP+L) = Z(ICC+J)
      Z(M+L) = Z(ICC+J+1) * CR
      L = L+1
   20 CONTINUE
      CALL WRITE(OUT,Z(IDP),NND,0)
      ICC = ICC + NND
   30 CONTINUE
 1000 RETURN
      END


text:
This editor can edit this entry and tell us a bit about themselves by clicking the Edit icon.

20090802 16:38:09 nbsp Welcome to the Wiki. Is The Imaginary Zebra a business? If so we can use the business template to include its address and other contact information. Users/JasonAller

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
20041221 20:52:46 nbsp Hi   Users/DanMasiel

20041221 23:50:54 nbsp Thanks for stopping by :)  Users/BryanBell

20050121 23:01:38 nbsp Theres an inversion layer right now that makes low elevation areas super cold and higher elevations hotter. At least thats my understanding of our current weather crisis.  Users/JackHaskel

20060114 22:58:10 nbsp Are you really http://www.math.ucdavis.edu/~bjwbell/ Batman? Users/RoyWright


text:
c     ****************************************************************
c     *                                                              *
c     *                      find_flat_work                          *
c     *                                                              *
c     ****************************************************************
c
      subroutine find_flat_work( combine_stream, combine_text,
     &                           local_count, matl_states_name )
c
      use global_data, only : termout, nodal_stresses_step_list,
     &   nodal_strains_step_list, element_stresses_step_list,
     &   element_strains_step_list, element_states_step_list,
     &   found_node_stresses, found_node_strains, found_ele_stresses,
     &   found_ele_strains, found_ele_states, last_step_no
      implicit none
c
c                parameter declarations
c                ----------------------
c
      integer :: local_count
      logical :: combine_stream, combine_text
      character(len=*) :: matl_states_name
c
c                local variables
c                ---------------
c
      integer :: step
      logical :: l1, l2, l3, l4, l5
      logical, parameter :: local_debug = .false.,
     &                      local_debug2 = .false.
      character(len=80) :: node_stress_fname, node_strain_fname,
     &                     element_stress_fname, element_strain_fname,
     &                     element_states_fname
      character(len=80) :: dummy_file_name
c
c                examine all files in directory to build lists
c                of step numbers to be combined for each type
c                of result file
c
      do step = 1, last_step_no
        nodal_stresses_step_list(step)   = 0
        nodal_strains_step_list(step)    = 0
        element_stresses_step_list(step) = 0
        element_strains_step_list(step)  = 0
        element_states_step_list(step)   = 0
      end do
c
      local_count = 0
      found_node_stresses = .false.
      found_node_strains  = .false.
      found_ele_stresses  = .false.
      found_ele_strains   = .false.
      found_ele_states    = .false.
c
      do step = 1, last_step_no
c
        call build_flat_file_names( step, 0, node_stress_fname,
     &                         node_strain_fname, element_stress_fname,
     &                         element_strain_fname, matl_states_name,
     &                         element_states_fname, termout,
     &                         combine_stream, combine_text,
     &                         dummy_file_name, 1, 1 )
        inquire( file=node_stress_fname,exist=l1 )
        inquire( file=node_strain_fname,exist=l2 )
        inquire( file=element_stress_fname,exist=l3 )
        inquire( file=element_strain_fname,exist=l4 )
        inquire( file=element_states_fname,exist=l5 )
c
        if( local_debug2 .and. step <=  30 ) then
            write(termout,9000) step, node_stress_fname,
     &                         node_strain_fname, element_stress_fname,
     &                         element_strain_fname, matl_states_name,
     &                         element_states_fname
            write(termout,9010) l1, l2, l3, l4, l5
        end if
c
        if( l1 ) then
            nodal_stresses_step_list(step) = 1
            found_node_stresses = .true.
        end if
        if( l2 ) then
            nodal_strains_step_list(step) = 1
            found_node_strains = .true.
        end if
        if( l3 ) then
            element_stresses_step_list(step) = 1
            found_ele_stresses = .true.
        end if
        if( l4 ) then
            element_strains_step_list(step)  = 1
            found_ele_strains = .true.
        end if
        if( l5 ) then
            element_states_step_list(step)  = 1
            found_ele_states = .true.
        end if
        if( l1 .or. l2 .or. l3 .or. l4 .or. l5 ) then
                 local_count = local_count + 1
        end if
c
      end do
c
      if( .not. local_debug ) return
c
      write(termout,*) "... list of flat files to process ..."
      write(termout,*) '      total number of files: ', local_count
      do step = 1, last_step_no
        if( nodal_stresses_step_list(step) == 1 )
     &      write(termout,*) '        step: ',step,' nodal stresses'
        if( nodal_strains_step_list(step) == 1 )
     &      write(termout,*) '        step: ',' nodal strains'
        if( element_stresses_step_list(step) == 1 )
     &      write(termout,*) '        step: ',step,' element stresses'
       if( element_strains_step_list(step) == 1 )
     &      write(termout,*) '        step: ',step,' element strains'
       if( element_states_step_list(step) == 1 )
     &      write(termout,*) '        step: ',step,' element states'
      end do
c
      return
c
 9000 format("... Flat file names for step: ",i7,/,10(10x,a,/))
 9010 format("    Exit flags for step: ",10l2)
c
      end subroutine find_flat_work
c **********************************************************************
c *                                                                    *
c *     support routine:      build_flat_file_names                    *
c *                                                                    *
c **********************************************************************
c
c
      subroutine build_flat_file_names( step, myid, node_stress_fname,
     &                         node_strain_fname, element_stress_fname,
     &                         element_strain_fname, matl_states_name,
     &                         element_states_fname, termout,
     &                         combine_stream, combine_text,
     &                         file_name, result_type, data_type  )
      implicit none
c
c
c                parameter declarations
c                ----------------------
c
      integer :: termout, step, myid, result_type, data_type
      logical :: combine_stream, combine_text
      character(len=*) :: node_stress_fname, node_strain_fname,
     &                    element_stress_fname, element_strain_fname,
     &                    element_states_fname, file_name,
     &                    matl_states_name
c
c                local variables
c                ---------------
c
      integer :: fcount, nchars
      character(len=7) :: step_id
      character(len=4) :: proc_id
      character(len=6) :: format
      logical :: local_debug
c
      local_debug              = .false.
      node_stress_fname(1:)    = ' '
      node_strain_fname(1:)    = ' '
      element_stress_fname(1:) = ' '
      element_strain_fname(1:) = ' '
      element_states_fname(1:) = ' '
c
      write(step_id,9000) step
      write(proc_id,9100) myid
c
      if( combine_stream ) then
         format = 'stream'
         fcount = 6
      elseif( combine_text ) then
         format = 'text'
         fcount = 4
      else
         write(termout,*)
     &       '>> FATAL ERROR: routine build_flat_file_names'
         write(termout,*) '>>   invalid format flag...'
         write(termout,*) '>>   job terminated...'
         stop
      end if
c
      node_stress_fname(1:)    = 'wn' // 's' //
     &                           step_id // '_' //
     &                           format(1:fcount) // "." // proc_id
      node_strain_fname(1:)    = 'wn' // 'e' //
     &                           step_id // '_' //
     &                           format(1:fcount) // "." // proc_id
c
      element_stress_fname(1:) = 'we' // 's' //
     &                           step_id // '_' //
     &                           format(1:fcount) // "." // proc_id
      element_strain_fname(1:) = 'we' // 'e' //
     &                           step_id // '_' //
     &                           format(1:fcount) // "." // proc_id
      nchars = len_trim( matl_states_name )
      element_states_fname(1:) = 'we' // 'm' //
     &                           step_id // '_' // format(1:fcount) //
     &                           '_' // matl_states_name(1:nchars) //
     &                           '.' // proc_id
c
c          data_type = 0   (element results)
c                    = 1   (nodal results)
c
c          result_type  = 0 (stress results)
c                       = 1 (strain results)
c                       = 2 (states results)
c
      if( data_type .eq. 0 ) then
        if( result_type .eq. 1 ) then
           file_name(1:) = element_strain_fname(1:)
        elseif( result_type .eq. 0 ) then
           file_name(1:) = element_stress_fname(1:)
        elseif( result_type .eq. 2 ) then
           file_name(1:) = element_states_fname(1:)
        else
           write(termout,9300)
           stop
        end if
      elseif( data_type .eq. 1 ) then
        if( result_type .eq. 1 ) then
           file_name(1:) = node_strain_fname(1:)
        elseif( result_type .eq. 0 ) then
           file_name(1:) = node_stress_fname(1:)
        else
           write(termout,9300)
           stop
        end if
      else
        write(termout,9400)
        stop
      end if
c
      if( local_debug ) then
        write(termout,*) ' >> local debug in build_flat_file_names...'
        write(termout,9200) node_stress_fname, node_strain_fname,
     &                      element_stress_fname, element_strain_fname,
     &                      element_states_fname
      end if
c
      return
c
 9000 format( i7.7 )
 9100 format( i4.4 )
 9200 format( 5x,' > file names: ', 5(/10x,a) )
 9300 format( '>> FATAL ERROR: routine build_flat_file_names...',
     &  /,    '                invalid results_type',
     &  /,    '                job terminated...' )
 9400 format( '>> FATAL ERROR: routine build_file_names...',
     &  /,    '                invalid data_type',
     &  /,    '                job terminated...' )
c
      end
c **********************************************************************
c *                                                                    *
c *          do_flat_element_step                                      *
c *                                                                    *
c **********************************************************************
c
c
      subroutine do_flat_element_step( step, combine_stream,
     &   combine_text, termout,  results_type, file_no, missing_flg,
     &   matl_state_name, did_flat_stream, did_flat_text,
     &   consistent_num_values )
c
      use global_data, only :  last_proc_id, num_model_elems
      implicit none
c
c
c                parameter declarations
c                ----------------------
c
      integer :: step, termout, results_type, file_no,
     &           consistent_num_values
      logical :: combine_stream, combine_text, missing_flg,
     &           did_flat_stream, did_flat_text
      character(len=*) :: matl_state_name
c
c                local variables
c                ---------------
c
      integer :: proc_id, open_status, num_values,
     &           allocate_status, num_elements,
     &           element_id,miss_count, num_header_lines
      integer, parameter :: max_header_lines = 50
      character(len=80) :: node_stress_fname, node_strain_fname,
     &                     element_stress_fname, element_strain_fname,
     &                     element_states_fname,
     &                     header_lines(max_header_lines)
      character(len=80) :: file_name
      logical :: local_debug, fatal, print_header, print
c
      double precision, allocatable,   dimension(:,:) :: element_values
      double precision, parameter :: zero = 0.0d00
      integer, allocatable, dimension(:)  :: element_list
c
      local_debug = .false.
c
c                 open the root processor file for this load step.
c                 get information about step then close down this
c                 file.
c
      file_name(1:) = " "
      call build_flat_file_names( step, 0, node_stress_fname,
     &                       node_strain_fname, element_stress_fname,
     &                       element_strain_fname, matl_state_name,
     &                       element_states_fname, termout,
     &                       combine_stream, combine_text,
     &                       file_name, results_type, 0 )
      if( local_debug ) write(termout,9700) step,
     &                  file_name(1:len_trim(file_name)), file_no
c
c                open the element result file for processor
c                zero. read the header records to get the
c                number of data values per element and the number
c                of elements. each result file for other processors
c                must have the same values else bad.
c
      call open_old_flat_result_file( file_no, file_name,
     &      combine_stream, combine_text, termout, 1, open_status, 0 )
c
      num_header_lines = 0
      if( combine_text ) call read_flat_header( file_no, header_lines,
     &               num_header_lines, max_header_lines, termout )
c
      if( local_debug )
     &     write(termout,*) '.... num_header_lines: ', num_header_lines
c
c                 allocate space to hold the data values for all
c                 elements in the mode. keep track of which element
c                 values have been read from the files for the step.
c
c                 results_type = 0,1,2 (stresses, strains, states)
c
      if( combine_stream ) read(file_no)   num_values
      if( combine_text )   read(file_no,*) num_values
      close( unit=file_no )
      if( consistent_num_values == -1 )
     &    consistent_num_values = num_values
      if( num_values .ne. consistent_num_values ) then
        write(termout,9730) num_values, consistent_num_values,
     &                      file_name(1:len_trim(file_name))
        stop
      end if
      num_elements = num_model_elems
      if( local_debug ) write(termout,9710) num_values
c
      allocate( element_values(num_values,num_elements),
     &          stat=allocate_status)
      call check_allocate_status( termout, 1, allocate_status )
      allocate( element_list(num_elements), stat=allocate_status )
      call check_allocate_status( termout, 2, allocate_status )
      element_list(1:num_elements) = 0
      element_values(1:num_values,1:num_elements) = zero
      if( local_debug ) write(termout,9720) num_elements
c
c                loop over all the separate processor files of
c                element results for this load step.
c                open the file and read the
c                element results into the single array here. keep
c                track if an element appears other than once... bad.
c
      do proc_id = 0, last_proc_id
c
        call build_flat_file_names( step, proc_id, node_stress_fname,
     &                       node_strain_fname, element_stress_fname,
     &                       element_strain_fname, matl_state_name,
     &                       element_states_fname, termout,
     &                       combine_stream, combine_text,
     &                       file_name, results_type, 0 )
        if( local_debug ) write(termout,9705) proc_id,
     &                    file_name(1:len_trim(file_name)), file_no
c
c                open the element result file for the processor.
c                read all data for a binary or a formmated file. if
c                the open fails, keep looking for processor files
c                until max_procs in case user blew away a
c                processor file.
c
        call open_old_flat_result_file( file_no, file_name,
     &       combine_stream, combine_text, termout, 2, open_status, 1 )
        if( open_status .ne. 0 ) cycle
c
c               skip element results header for this processor.
c
        if( combine_text ) call skip_comment_lines( file_no )
c
c                read data values until end of file. stuff into the
c                single element results table. update count of
c                element appearances.
c
       call read_flat_element_values( file_no, element_values,
     &      num_values, num_elements, element_list, termout,
     &      combine_stream, combine_text )
     &
       close( unit=file_no )
c
c              end loop over processor files for step
      end do
c
c                check that each element has exactly 1 set of
c                values read from all the processor file
c
      fatal = .false.
      print_header = .true.
      miss_count = 0
      print = .true.
c
      do element_id = 1, num_elements
        if( element_list(element_id) .eq. 0 ) missing_flg = .true.
      end do

      if( print ) then
        do element_id = 1, num_elements
         if(  element_list(element_id) .ne. 1 ) then
           if(  print_header ) write(termout,9300)
           print_header = .false.
           write(termout,9400) element_id
           miss_count = miss_count + 1
           if(  miss_count .eq. 5 ) then
              write(termout,9301)
              exit
           end if

         end if
        end do
      end if
c
c               element results data looks ok. create a new,
c               single flat results file for this step.
c               write in all the data.
c
c
      if( combine_text )   did_flat_text   = .true.
      if( combine_stream ) did_flat_stream = .true.
c
      call open_new_flat_result_file( file_no, file_name,
     &     combine_stream, combine_text, termout, 1 )
c
      call write_flat_element_results( file_no, num_header_lines,
     &    header_lines, num_values, num_elements, element_values,
     &    termout, combine_stream, combine_text, element_list )
c
      close(unit=file_no,status='keep')
c
      deallocate( element_values, element_list, stat=allocate_status)
      call check_deallocate_status( termout, 1, allocate_status )
c
      return
c
 9300 format(15x,'> Note: The following elements have an appearance',
     &     /,15x '        count not equal to 1 in Patran result',
     &     /,15x '        files. This could result from any',
     &     /,15x,'        elements which do not emit results.')
 9301 format(15x,'> Note: Further messages of this nature',
     &     /,15x '        are suppressed...')
 9400 format(15x,i7)
c
 9700 format(5x,'> local debug in do_flat_element_step: ',
     & /,10x,'step: ',i7,/,10x,'file_name: ',a,
     & /,10x,'file_no: ',i5 )
 9705 format(5x,'> local debug in do_flat_element_step: ',
     & /,10x,'proc_id: ',i6,/,10x,'file_name: ',a,
     & /,10x,'file_no: ',i5 )
 9710 format(10x,'number of data values: ',i5)
 9720 format(5x,'> data arrays allocated. no. elements: ',i7)
 9730 format( '>> FATAL ERROR: routine do_flat_element_step...',
     &  /,    '                inconsistent values.',
     &  /,    '                num_values, num_consistent: ',2i10,
     &  /,    '                file: ',a,
     &  /,    '                job terminated...' )
c
      end

c **********************************************************************
c *                                                                    *
c *          open_old_flat_result_file                                 *
c *                                                                    *
c **********************************************************************
c
c
      subroutine open_old_flat_result_file( file_no, file_name,
     &             combine_stream, combine_text, termout,
     &             location, open_status, dowhat )
      implicit none
c
c                  parameter declarations
c                  ----------------------
c
      integer :: file_no, termout, location, open_status, dowhat
      character(len=*) :: file_name
      logical :: combine_stream, combine_text
c
c                  local declarations
c                  ------------------
c
c
      if( combine_stream ) then
        open( unit=file_no, file=file_name, status='old',
     &        access='stream', form='unformatted',
     &        iostat= open_status )
      elseif( combine_text ) then
        open( unit=file_no, file=file_name, status='old',
     &        access='sequential', form='formatted',
     &        iostat= open_status )
      else
        write(termout,9200) location
        stop
      end if
c
c                 dowhat = 0 : check open status and stop job if
c                              open failed
c                        = 1 : don't check open status
c
      if( dowhat .eq. 1 ) return
c
      if( open_status .ne. 0 ) then
        write(termout,9100) location
        stop
      end if
c
      return
c
 9100 format( '>> FATAL ERROR: failed flat file open @ ',i2,
     &  /,    '                job terminated...' )
 9200 format( '>> FATAL ERROR: invalid system state @ ',i2,
     &  /,    '                job terminated...' )
c
      end
c **********************************************************************
c *                                                                    *
c *          read_flat_header                                          *
c *                                                                    *
c **********************************************************************
c
c
      subroutine read_flat_header( file_no, header_lines,
     &                             num_header_lines, max_header_lines,
     &                             termout )
      implicit none
c
c                  parameter declarations
c                  ----------------------
c
      integer :: file_no, termout, num_header_lines, max_header_lines
      character(len=80) :: header_lines(max_header_lines)
c
c                  local declarations
c                  ------------------
c
      character :: first_char*1, line*80
c
      header_lines(1:max_header_lines) = " "
      num_header_lines = 0
c
      do
        read(file_no,fmt="(a80)") line
        first_char = line(1:1)
        if( first_char .eq. "#"  .or.
     &      first_char .eq. "!"  .or.
     &      first_char .eq. "c"  .or.
     &      first_char .eq. "C"  ) then
          num_header_lines = num_header_lines + 1
          header_lines(num_header_lines) = line
          cycle
        end if
        backspace(file_no) ! next line has number of values for states
c                            processing
        return
       end do
c
      return
      end

c **********************************************************************
c *                                                                    *
c *          read_flat_element_values                                  *
c *                                                                    *
c **********************************************************************
c
c
      subroutine read_flat_element_values(
     &           file_no, element_values, num_values, num_elements,
     &           element_list, termout, combine_stream,
     &           combine_text )
      implicit none
c
c                  parameter declarations
c                  ----------------------
c
      integer :: file_no, num_values, num_elements, element_list(*),
     &           termout
      logical :: combine_stream, combine_text
      double precision :: element_values(num_values,num_elements)
c
c                  local declarations
c                  ------------------
c
      integer :: read_status, element_id, now_num_values
      double precision, allocatable :: element_data(:)
c
c              next line (value) on both file types should be
c              the number of state variables. This value must
c              be same as value passed in or we have an
c              error somewhere.
c
      if( combine_stream ) read(file_no,iostat=read_status)
     &                       now_num_values
      if( combine_text )   read(file_no,*,iostat=read_status)
     &                       now_num_values
      if( read_status .ne. 0 ) then
          write(termout,9140)
          stop
      end if
      if( now_num_values .ne. num_values ) then
          write(termout,9150) num_values, now_num_values
        stop
      end if
c
      allocate( element_data( num_values) )
c
      do
        if( combine_stream ) then
          read(unit=file_no,iostat=read_status) element_id,
     &        element_data
        elseif( combine_text ) then
          read(unit=file_no,iostat=read_status,fmt=9330)
     &         element_id, element_data
        else
         write(termout,9200)
         stop
        end if
c
        if( read_status .ne. 0 ) return
        if( element_id .le. num_elements ) then
          element_values(1:num_values,element_id) =
     &                    element_data(1:num_values)
          element_list(element_id) = element_list(element_id) + 1
        else
          write(termout,9130) element_id
          stop
        end if
c
      end do
c
      deallocate( element_data )
      return
c
 9130 format( '>> FATAL ERROR: routine read_flat_element_data...',
     &  /,    '                internal error. bad element #: ',i10,
     &  /,    '                job terminated...' )
 9140 format( '>> FATAL ERROR: routine read_flat_element_data...',
     &  /,    '                internal error. ',
     &  /,    '                job terminated...' )
 9150 format( '>> FATAL ERROR: routine read_flat_element_data...',
     &  /,    '                internal error. number of values',
     &  /,    '                inconsistent: ',2i10,
     &  /,    '                job terminated...' )
 9330 format(i8,30e15.6)
 9200 format( '>> FATAL ERROR: invalid system state in ',
     &  /,    '                routine: read_flat_element_data',
     &  /,    '                job terminated...' )
c
      end

c **********************************************************************
c *                                                                    *
c *          open_new_flat_result_file                                 *
c *                                                                    *
c **********************************************************************
c
c
      subroutine open_new_flat_result_file( file_no, file_name,
     &         combine_stream, combine_text, termout, location )
      implicit none
c
c                  parameter declarations
c                  ----------------------
c
      integer :: file_no, termout, location
      character(len=*) :: file_name
      logical :: combine_stream, combine_text
c
c                  local declarations
c                  ------------------
c
      integer :: open_status, dot_pos
      character(len=80) :: local_fname
c
c                  remove the .proc_id from file name
c
      dot_pos = index( file_name, ".")
      if( dot_pos .eq. 0 ) then
        write(termout,9300) file_name(1:len_trim(file_name))
        stop
      end if
      local_fname(1:) = file_name(1:dot_pos-1)
c
      if( combine_stream ) then
        open( unit=file_no, file=local_fname, status='unknown',
     &        access='stream', form='unformatted', iostat= open_status )
      elseif( combine_text ) then
        open( unit=file_no, file=local_fname, status='unknown',
     &        access='sequential', form='formatted',
     &        iostat= open_status )
      else
        write(termout,9200) location
        stop
      end if
c
      if( open_status .ne. 0 ) then
        write(termout,9100) location
        stop
      end if
c
      return
c
 9100 format( '>> FATAL ERROR: failed open_new_flat_result_file @ ',i2,
     &  /,    '                job terminated...' )
 9200 format( '>> FATAL ERROR: invalid system state @ ',i2,
     &  /,    '                file: open_new_flat_result_file',
     &  /,    '                job terminated...' )
 9300 format( '>> FATAL ERROR: invalid system state in routine:',
     &  /,    '                open_new_flat_result_file',
     &  /,    '                bad file name: ',a,
     &  /,    '                job terminated...' )
c
      end
c **********************************************************************
c *                                                                    *
c *          write_flat_element_results                                *
c *                                                                    *
c **********************************************************************
c
c
      subroutine write_flat_element_results(
     &         file_no, num_header_lines, header_lines, num_values,
     &         num_elements, element_values, termout, combine_stream,
     &         combine_text, element_list )
      implicit none
c
c                  parameter declarations
c                  ----------------------
c
      integer :: file_no, termout, num_values, num_elements,
     &           num_header_lines
      logical :: combine_stream, combine_text
      double precision :: element_values(num_values,num_elements)
      integer :: element_list(num_elements)
      character(len=80) :: header_lines(*)
c
c                  local declarations
c                  ------------------
c
      integer ::  element_id, write_status, i
c
      if( combine_stream ) then
c
         do element_id = 1, num_elements
          if( element_list(element_id) .eq. 1 )
     &        write(unit=file_no,iostat=write_status)
     &        (element_values(i,element_id), i = 1, num_values)
          call check_write_status( file_no, write_status, termout, 10 )
         end do
c
      elseif( combine_text ) then
c
         do i = 1, num_header_lines
          write(file_no,9000) trim(header_lines(i))
         end do
         do element_id = 1, num_elements
          if( element_list(element_id) .eq. 1 )
     &       write(unit=file_no,iostat=write_status,fmt=9200)
     &         ( element_values(i,element_id), i = 1, num_values )
          call check_write_status( file_no, write_status, termout, 8 )
         end do
      else
        write(termout,9300)
        stop
c
      end if
c
      return
c
 9000 format(a)
 9300 format( '>> FATAL ERROR: invalid system state in routine:',
     &  /,    '                write_element_results',
     &  /,    '                job terminated...' )
 9200 format(30e15.6)
c
      end

c **********************************************************************
c *                                                                    *
c *          do_flat_nodal_step                                        *
c *                                                                    *
c **********************************************************************
c
c
      subroutine do_flat_nodal_step( step, combine_stream, combine_text,
     &      termout, sig_eps_type, file_no, missing_flg,
     &      did_flat_stream, did_flat_text, consistent_num_values )
c
      use global_data, only :  num_model_nodes, num_model_elems,
     &                         last_proc_id
      implicit none
c
c
c                parameter declarations
c                ----------------------
c
      integer :: step, termout, sig_eps_type, file_no,
     &           consistent_num_values
      logical :: combine_stream, combine_text, missing_flg,
     &           did_flat_stream, did_flat_text
c
c                local variables
c                ---------------
c
      integer:: proc_id, open_status, num_values, allocate_status,
     &          num_nodes, node_id, miss_count, num_header_lines
      integer, allocatable, dimension(:)  :: count_list
      integer, parameter :: max_header_lines = 50
c
      logical ::  local_debug, print_header, print
c
      double precision, allocatable,  dimension(:,:) :: node_values
      double precision, parameter :: zero = 0.0d00
c
      character(len=80) :: node_stress_fname, node_strain_fname,
     &               element_stress_fname, element_strain_fname,
     &               element_states_fname, matl_states_fname,
     &               header_lines(max_header_lines), file_name
c
      local_debug = .false.
c
c                 get name of root processor file for this load step.
c
      matl_states_fname(1:) = "..dummy.."
      call build_flat_file_names( step, 0, node_stress_fname,
     &    node_strain_fname, element_stress_fname,
     &    element_strain_fname, matl_states_fname,
     &    element_states_fname, termout,  combine_stream, combine_text,
     &    file_name, sig_eps_type, 1 )
c
      if( local_debug ) write(termout,9700) step, file_name(1:20),
     &                                      file_no
c
c                open the node result file for processor
c                zero. read/store header lines for text result
c
      call open_old_flat_result_file( file_no, file_name,
     &    combine_stream, combine_text, termout, 1, open_status, 0 )
      num_header_lines = 0
      if( combine_text ) call read_flat_header( file_no, header_lines,
     &                       num_header_lines, max_header_lines,
     &                       termout )
c
      if( local_debug )
     &     write(termout,*) '.... num_header_lines: ', num_header_lines
c
c                 allocate space to hold the data values for all
c                 nodes in the model. keep track of which nodes
c                 have been read from the files for the step.
c
      if( combine_stream ) read(file_no)   num_values
      if( combine_text )   read(file_no,*) num_values
      close( unit=file_no )
      if( consistent_num_values == -1 )
     &    consistent_num_values = num_values
      if( num_values .ne. consistent_num_values ) then
        write(termout,9730) num_values, consistent_num_values,
     &                      file_name(1:len_trim(file_name))
        stop
      end if

      num_nodes = num_model_nodes
      allocate( node_values(num_nodes,num_values), stat=allocate_status)
      call check_allocate_status( termout, 10001, allocate_status )
c
      allocate( count_list(num_nodes), stat=allocate_status )
      call check_allocate_status( termout, 1002, allocate_status )
c
      node_values(1:num_nodes,1:num_values) = zero
      count_list(1:num_nodes)               = 0
      if( local_debug ) write(termout,9720) num_nodes
c
c                loop over all the separate rank files of
c                node results for this load step. open the file and read
c                the node results into the single array here. keep
c                track of how many times a node appears in result files
c                to enable averaging.
c
      do proc_id = 0, last_proc_id ! MPI ranks
c
        call build_flat_file_names( step, proc_id, node_stress_fname,
     &     node_strain_fname, element_stress_fname,
     &     element_strain_fname,  matl_states_fname,
     &     element_states_fname, termout, combine_stream, combine_text,
     &     file_name, sig_eps_type, 1 )
c
c                open the node result file for the processor.
c                read all data for a stream or text file. if
c                the open fails, keep looking for processor files
c                until last_proc_id in case user blew away a
c                processor file.
c
        call open_old_flat_result_file( file_no, file_name,
     &     combine_stream, combine_text, termout, 2, open_status, 1 )
        if( open_status .ne. 0 ) cycle
        if( local_debug ) write(termout,9705) proc_id, file_name(1:20),
     &                                        file_no
c
c               skip header lines for text results files
c
       if( combine_text ) call skip_comment_lines( file_no )
c
c                read node results for this processor. check
c                read data values until end of file. stuff into the
c                single node results table. update count of
c                node appearances.
c
       call read_flat_node_values( file_no, node_values, num_values,
     &    num_nodes, count_list, termout, combine_stream, combine_text )
     &
       close( unit=file_no )
c
      end do  !  over processor files
c
c                issue warning messages for nodes that have no results.
c                disable printing of messages in curent version based
c                on a global message in calling routine.
c
      print_header = .true.
      miss_count   = 0
      print        = .true.
c
      do node_id = 1, num_nodes
         if( count_list(node_id) .eq. 0 ) missing_flg = .true.
      end do
c
      if( print ) then
       do node_id = 1, num_nodes
         if( count_list(node_id) .eq. 0 ) then
           if( print_header ) write(termout,9300)
           print_header = .false.
           write(termout,9400) node_id
           miss_count = miss_count + 1
           if( miss_count .eq. 50 ) then
              write(termout,9301)
              exit
           end if
         end if ! on node_list
       end do
      end if ! on print
c
c               node results data looks ok. create a new,
c               single flat results file for this step.
c               write in all the data.
c
      if( combine_stream ) did_flat_stream = .true.
      if( combine_text )   did_flat_text   = .true.
c
      call open_new_flat_result_file( file_no, file_name,
     &     combine_stream, combine_text, termout, 2 )
c
      call write_flat_node_results( file_no, num_header_lines,
     &      header_lines, num_values, num_nodes, node_values,
     &      count_list, termout, combine_stream, combine_text,
     &      sig_eps_type )
c
      close( unit=file_no )
c
      deallocate( node_values, count_list, stat=allocate_status )
      call check_deallocate_status( termout, 1003, allocate_status )
c
      return
c
 9300 format(10x,'> Note: The following nodes have no results',
     &     /,10x '        in the flat result files...')
 9301 format(10x,'> Note: Further error messages of this nature',
     &     /,10x '        are suppressed...')
 9400 format(15x,i7)
c
 9700 format(5x,'> local debug in do_flat_node_step: ',
     & /,10x,'step: ',i7,/,10x,'file_name: ',a20,
     & /,10x,'file_no: ',i5 )
 9705 format(5x,'> local debug in do_flat_node_step: ',
     & /,10x,'proc_id: ',i6,/,10x,'file_name: ',a20,
     & /,10x,'file_no: ',i5 )
 9720 format(5x,'> data arrays allocated. no. nodes: ',i7)
 9730 format( '>> FATAL ERROR: routine do_flat_nodal_step...',
     &  /,    '                inconsistent values.',
     &  /,    '                num_values, num_consistent: ',2i10,
     &  /,    '                file: ',a,
     &  /,    '                job terminated...' )
c
      end
c **********************************************************************
c *                                                                    *
c *          read_flat_node_values                                     *
c *                                                                    *
c **********************************************************************
c
c
      subroutine read_flat_node_values( file_no, node_values,
     &             num_values, num_nodes, node_counts, termout,
     &             combine_stream, combine_text )
      implicit none
c
c                  parameter declarations
c                  ----------------------
c
      integer :: file_no, num_values, num_nodes, node_counts(*),
     &           termout
      logical :: combine_stream, combine_text
      double precision :: node_values(num_nodes,num_values)
c
c                  local declarations
c                  ------------------
c
      integer :: read_status, node_id, lcount, now_num_values
      double precision :: node_data(100)
c
c              next line (value) on both file types should be
c              the number of state variables. This value must
c              be same as value passed in or we have an
c              error somewhere.
c
      if( combine_stream ) read(file_no,iostat=read_status)
     &                       now_num_values
      if( combine_text )   read(file_no,*,iostat=read_status)
     &                       now_num_values
      if( read_status .ne. 0 ) then
          write(termout,9140)
          stop
      end if
      if( now_num_values .ne. num_values ) then
          write(termout,9150) num_values, now_num_values
        stop
      end if
c
      do
        if( combine_stream ) then
          read(unit=file_no,iostat=read_status) node_id, lcount,
     &          node_data(1:num_values)
        elseif( combine_text ) then
          read(unit=file_no,iostat=read_status,fmt=9330) node_id,
     &          lcount, node_data(1:num_values)
        else
          write(termout,9200)
          stop
        end if
c
        if( read_status .ne. 0 ) return
        if( node_id .le. num_nodes ) then
          node_values(node_id,1:num_values) =
     &                 node_values(node_id,1:num_values) +
     &                 node_data(1:num_values)
          node_counts(node_id) = node_counts(node_id) + lcount
        else
          write(termout,9130)
          stop
        end if
c
      end do
c
      return
c
 9130 format( '>> FATAL ERROR: routine read_flat_node_values...',
     &  /,    '                internal error @ 5',
     &  /,    '                job terminated...' )
 9140 format( '>> FATAL ERROR: routine read_flat_nodal_data...',
     &  /,    '                internal error. ',
     &  /,    '                job terminated...' )
 9150 format( '>> FATAL ERROR: routine read_flat_nodal_data...',
     &  /,    '                internal error. number of values',
     &  /,    '                inconsistent: ',2i10,
     &  /,    '                job terminated...' )
 9330 format(2i9,30d15.6)
 9200 format( '>> FATAL ERROR: invalid system state in ',
     &  /,    '                routine: read_flat_node_values',
     &  /,    '                job terminated...' )
c
      end

c **********************************************************************
c *                                                                    *
c *          write_flat_node_results                                   *
c *                                                                    *
c **********************************************************************
c
c
      subroutine write_flat_node_results( file_no, num_header_lines,
     &         header_lines, num_values, num_nodes, node_values,
     &         node_counts,termout, combine_stream, combine_text,
     &         sig_eps_type )
      implicit none
c
c                  parameter declarations
c                  ----------------------
c
      integer :: file_no, termout, num_values, num_nodes,
     &           sig_eps_type, num_header_lines
      integer :: node_counts(*)
      logical :: combine_stream, combine_text
      double precision :: node_values(num_nodes,num_values)
      character(len=80) :: header_lines(*)
c
c                  local declarations
c                  ------------------
c
      integer :: node_id, write_status, i
      double precision :: scale
      double precision, parameter :: zero = 0.0d0, one = 1.0d0
c
      do node_id = 1, num_nodes
        scale = zero
        if( node_counts(node_id) .gt. 0 )
     &       scale = one /  dble( node_counts(node_id) )
        node_values(node_id,1:num_values) =
     &            node_values(node_id,1:num_values) * scale
      end do
c
      call compute_extra_sig_eps_values( node_values, num_nodes,
     &                                   sig_eps_type )
c
      if( combine_stream ) then
c
         do node_id = 1, num_nodes
           write(unit=file_no,iostat=write_status)
     &          node_values(node_id,1:num_values)
           call check_write_status( file_no, write_status, termout,
     &                              200 )
         end do
c
      elseif( combine_text ) then
c
         do i = 1, num_header_lines
          write(file_no,9000) trim( header_lines(i) )
         end do
c
         do node_id = 1, num_nodes
            write(unit=file_no,iostat=write_status,fmt=9330)
     &            node_values(node_id,1:num_values)
            call check_write_status( file_no, write_status, termout,
     &                               201 )
         end do
c
      else
c
         write(termout,9200)
         stop
c
      end if
c
      return
c
 9000 format(a)
 9200 format( '>> FATAL ERROR: invalid system state in routine:',
     &  /,    '                write_flat_node_results',
     &  /,    '                job terminated...' )
 9330 format(30e15.6)
c
      end