PROGRAM SCREAD4

c same as pauls but change fnam='/home/disk3/mmccarty/calib/seacat/' to
c       fnam='/home/summer2/data/calib/seacat/zeta'
c
c	This version will also work with files that have been
c	dumped with SEASOFT.
c
c	version 4 reads calib. coef. from files
c	[FREITAG.CALIB.SEACAT]:xxx.Tcoef and .Ccoef where xxx
c	is the isn and uses the last entry from the file.
c	Vers. 4 also can process units with pressure
c
c	version 3 differs from ver. s only in ability to process
c	isn>99
C
C	READS RAW SEACAT DATA AND CONVERTS TO TEMPERATURE IN
C	DEGREES AND CONDUCTIVITY/SALINITY IN xxxx/ppt.  
c
C	Scans flagged with low battery (5A5A5A5A) filled with 999.999
C
C	ASSUMES THAT WHEN DATA WAS READ FROM SEACAT THE FOLLOWING
C	COMMANDS WERE USED IN ORDER
C	      DS
C	      DL (ONLY FOR PRE-VERSION 2.0 UNITS)
C	      DH 
C	      DDiii,jjj  (more than one DD okay)
C
C   f77 -o scread4 scread4.f -lhpf
	CHARACTER*70 LINE
	CHARACTER*80 FNAMI,FNAMO
	CHARACTER*50 head
	CHARACTER*8 PROJ,PROG
	CHARACTER*4 MOOR,IFORM
	DIMENSION IDATE(5,10),KDATE(5),DELT(10),NCHG(10)
	dimension mdate(5,6),xdate(6),jdate(5)
C
	DATA PROG /' SCREAD4'/
	data lbat,xbad /'5A5A'X,999.999/
	data jdate /0,0,1,1,87/
        character*1 icr
        icr=char(13)
	LMAX=70
	F1000=1000.
	CS=4.2914
100	READ (5,*,END=900,err=900) FNAMI,FNAMO,PROJ,MOOR,DEPTH,nmark,
     1	((mdate(j,m),j=1,5),m=1,nmark)
	nbi=notbnk(fnami,80)
	nbo=notbnk(fnamo,80)
	WRITE (6,101) PROG,FNAMI(1:nbi),FNAMO(1:nbo)
101	FORMAT(' PROGRAM',A//10X,'INPUT FILE NAME:  ',A/
     1	10X,'OUTPUT FILE NAME: ',A//10X,'INPUT FILE HEADERS:')
	call fOPEN(1,FNAMI,'OLD')
1011	NOUT=0
	nbat=0
	mark=1
110	READ (1,111,END=200) LINE
111	FORMAT(A)
	nb=notbnk(LINE,70)
        if(line(nb:nb).eq.icr) nb=nb-1
	WRITE (6,112) LINE(1:nb)
112	FORMAT(10X,A)
	II=INDEX(LINE(1:20),'SEACAT')
	IF(II.EQ.0) GO TO 110
c
c	decode the software version number and instrument serial number.
c	first see if data has been dumped using SEASOFT (* in col 1)
c	or if software version is 4.0 or greater
c
	isoft=0
	if(line(1:1).eq.'*') isoft=1
	nsmax=0
	II=INDEX(LINE,'.')-1
	read (LINE(II:),119) SVN
119	FORMAT(F3.1)
	if(svn.ge.4.0) isoft=1
	ii=INDEX(LINE,'NO.')+3
	read (line(ii:ii+4),*) isn
	iprs=0
120	READ (1,111,END=200) LINE
	nb=notbnk(LINE,70)
        if(line(nb:nb).eq.icr) nb=nb-1
	WRITE (6,112) LINE(1:nb)
	if(index(line,'pressure').gt.0) iprs=1
	II=INDEX(LINE,'samples =')
	IF(II.EQ.0) GO TO 120
c
c	decode the total number of samples recorded
c
	II=II+10
	IMAX=5
	IM=II+IMAX-1
	ICHR=LNBLNK(LINE(II:IM),IMAX)
	write(IFORM,121) ICHR
121	FORMAT('(I',I1')')
888	FORMAT(I10,A,A)
	read(LINE(II:IM),iform) NSMAX
	NDELT=0
c
c	if SEASOFT is used you need to back up to find start time.
c	(This only works if sample rate is constant)
c
	if(isoft.eq.1) then
	backspace 1
	backspace 1
	endif
130	READ (1,111,END=200) LINE
	nb=notbnk(LINE,70)
        if(line(nb:nb).eq.icr) nb=nb-1
	WRITE (6,112) LINE(1:nb)
	ii=index(line(1:4),'at')
	IF(ii.gt.0) THEN
c
c	decode the start time
c
	NDELT=NDELT+1
	read(line(ii:),131) IDATE(4,NDELT),IDATE(3,NDELT),
     1	IDATE(5,NDELT),IDATE(1,NDELT),IDATE(2,NDELT)
131	FORMAT(2X,3(1X,I2),1X,2(1X,I2))
	II=NOTBNK(LINE,LMAX)
	IL=II-4
	IF(SVN.LT.2.) THEN
	II=INDEX(LINE,'every')
	II=II+6
	ELSE
	II=INDEX(LINE,'val =')
	II=II+6
	ENDIF
	ICHAR=IL-II+1
	write(IFORM,121) ICHAR
	read(LINE(II:),iform) IDELT
	DELT(NDELT)=IDELT
	if(svn.ge.2.3) 	DELT(NDELT)=IDELT/60.
	GO TO 130
	ELSE
	IF(NDELT.EQ.0) GO TO 130
	ENDIF
	DEL0=0.0
	IF(NDELT.EQ.1) DEL0=DELT(1)
	IV=SVN
	if(iprs.ne.0) iv=-iv
	call fOPEN(2,FNAMO,'NEW')
1021	CALL fDATE(LINE)
	WRITE (head,138) PROJ,MOOR,ISN,IV,DEL0,DEPTH,PROG,LINE(4:10),
     1	LINE(20:24)
138	FORMAT(A,A,i4,I2,2F6.1,A,A,a)
	if(isn.le.999) head(13:13)='S'
	write(2,111)head
	WRITE (6,139) head 
139	FORMAT(//10X,'OUTPUT FILE HEAD: ',A)
	call getcal(isn,iprs)
	if(nmark.gt.0) then
	xdel=delt(1)/144.
	write (6,103) ((mdate(j,m),j=1,5),m=1,nmark)
103	format(//10x,'Time marks at ',5i3/(24x,5i3))
	write (6,104)
104	format(//)
	do 105 m=1,nmark
105	xdate(m)=delday(jdate,mdate(1,m))
	endif
140	READ (1,111,END=200) LINE
	nb=notbnk(LINE,70)
        if(line(nb:nb).eq.icr) nb=nb-1
	WRITE (6,112) LINE(1:nb)
	if(line(1:4).eq.'*END') then
	n1=0
	n2=nsmax-1
	else
	IF(LINE(1:4).NE.'S>DD'.and.line(1:4).ne.'S>dd') GO TO 140
	II=INDEX(LINE,',')
	IF(II.GT.0) THEN
	ICHAR=II-5
	write(IFORM,121) ICHAR
	read(LINE(5:II-1),iform) N1
	IL=LNBLNK(LINE(1:20),20)
	ICHAR=IL-II
	write(IFORM,121) ICHAR
	read(LINE(II+1:IL),iform) N2
	nstop=n2
	if(n2.gt.nsmax) n2=nsmax
	ELSE
	N1=0
	N2=NSMAX-1
	nstop=n2
	ENDIF
	endif
	WRITE (6,141) N1,N2
141	FORMAT(//10X,'READING SAMPLES ',I5,' TO ',I5)
	NCHG(1)=0
	NCHG(NDELT+1)=NSMAX
	IF(NDELT.GT.1) THEN
	DO 150 ND=2,NDELT	
	DAYS=DELDAY(IDATE(1,ND-1),IDATE(1,ND))
	NSAMP=DAYS*1440./DELT(ND-1)+.000001
150	NCHG(ND)=NCHG(ND-1)+NSAMP
	ENDIF
	ND=1
	LPR=1
	NTOT=N2-N1+1
	if(iprs.eq.0) then
	WRITE (6,1151)
1151	FORMAT(//5X,'TIME  DATE   TEMP   SALIN',8X,'T FREQ'
     1	8X,'C FREQ',10X,'COND'//)
	else
	WRITE (6,1152)
1152	FORMAT(//5X,'TIME  DATE   TEMP   SALIN   PRESS',8X,'T FREQ'
     1	8X,'C FREQ',10X,'COND   PCOUNT'//)
	endif
	DO 180 N=N1,N2
151	IF(N.GT.NCHG(ND+1)) THEN
	ND=ND+1
	IF(ND.GT.NDELT) GO TO 800
	WRITE (6,162) KDATE,T,S,NOUT
	WRITE (6,152) DELT(ND)
152	FORMAT(/10X,'DELTA-T CHANGES TO ',F6.0)
	xdel=delt(nd)/144.
	LPR=1
	GO TO 151
	ENDIF
	if(iprs.eq.0) then
	READ (1,153,END=200,err=200) ITHEX,ICHEX
153	FORMAT(3Z4)
	else
	READ (1,153,END=200,err=200) ITHEX,ICHEX,iphex
	endif
	if(ithex.eq.lbat.and.ichex.eq.lbat) then
	nbat=nbat+1
	tf=xbad
	t=xbad
	cf=xbad
	c=xbad
	s=xbad
	p=xbad
	lpr=1
	else
	TF=TFREQ(ITHEX)
	T=CATEMP(TF)
	CF=CFREQ(ICHEX)/F1000
	C=CATCON(CF,T)
	CR=C/CS
	S=SAL78(CR,T,DEPTH,0)
	if(iprs.ne.0) p=catprs(iphex)
	endif
	NN=N-NCHG(ND)+1
	CALL GETDT2(IDATE(1,ND),1,DELT(ND),NN,KDATE)
	NOUT=NOUT+1
155	if(mark.le.nmark) then
	yday=delday(jdate,kdate)
	ydel=yday-xdate(mark)
	if(ydel.gt.xdel) then
	mark=mark+1
	write (6,104)
	go to 155
	else
	if(ydel.ge.-xdel) lpr=1
	endif
	endif
	IF(NOUT.LE.10.OR.NOUT.GT.NTOT-10) LPR=1
	if(iprs.eq.0) then
	WRITE (2,161) KDATE,T,S,TF,CF,C,NOUT
161	FORMAT(1X,2I2,1X,3I2,2F8.3,3E14.6,I6)
	IF(LPR.NE.0) WRITE (6,162) KDATE,T,S,TF,CF,C,NOUT
162	FORMAT(5X,2I2,1X,3I2,2F8.3,3E14.6,I6)
	else
	WRITE (2,1161) KDATE,T,S,p,TF,CF,C,iphex,NOUT
1161	FORMAT(1X,2I2,1X,3I2,3F8.3,3E14.6,i8,I6)
	IF(LPR.NE.0) WRITE (6,1162) KDATE,T,S,p,TF,CF,C,iphex,NOUT
1162	FORMAT(5X,2I2,1X,3I2,3F8.3,3E14.6,i8,I6)
	endif
180	LPR=0
c
c	all this is in case they read more data that is recorded
c
	nmore=nstop-n2
	if(nmore.gt.0) then
	do n=1,nmore
	read(1,111,END=200) LINE
	enddo
	endif
	GO TO 140
200	WRITE (6,201) NOUT,nbat
201	FORMAT(//10X,I5,' RECORDS OUTPUT'/10x,i5,' Low Battery')
	CLOSE(UNIT=1)
	CLOSE(UNIT=2)
	call pej
	GO TO 100
800	WRITE (6,801) 
801	FORMAT(//10X,'ERROR IN TIME/DATE')
900	STOP
	END
	FUNCTION LNBLNK(ISTRNG,IMAX)
c
c	this version stops at either ' ' or any control char
c
	CHARACTER*(*) ISTRNG
	I=0
10	I=I+1
	IF(I.GT.IMAX) GO TO 20
	IF(ichar(istrng(i:i)).gt.32) GO TO 10
20	LNBLNK=I-1
	RETURN
	END
	FUNCTION TFREQ(IT)
	DATA A1,A2 /19.,2100./
	TFREQ=IT/A1+A2
	RETURN
	END
	FUNCTION CFREQ(IC)
	DATA A1,A2 /2100.,6.25E6/
	CFREQ=SQRT(IC*A1+A2)
	RETURN
	END
	subroutine getcal(isn,iprs)
	REAL*8 A,B,C,D,F0
	REAL*8 AA,BB,CC,DD
	REAL*4 Ap,Bp,Cp
	REAL*4 M
	character*90 aline
	character*50 fnam
	common /tcoef/ a,b,c,d,f0
	common /ccoef/ aa,bb,cc,dd,m
	common /pcoef/ ap,bp,cp
c changed from
c 	fnam='/home/disk3/mmccarty/calib/seacat/'
	fnam='/home/summer2/data/calib/seacat/zeta/'
	nb=notbnk(fnam,50)
	if(isn.lt.100) then
	write(fnam(nb+1:50),2) isn
2	format(i2,'.tcoef')
	ic=nb+4
	elseif(isn.lt.1000) then
	write(fnam(nb+1:50),3) isn
3	format(i3,'.tcoef')
	ic=nb+5
	else
	write(fnam(nb+1:50),4) isn
4	format(i4,'.tcoef')
	ic=nb+6
	endif
	open(unit=3,file=fnam,type='old',err=100)
10	read (3,11,end=12) aline
11	format(a)
	go to 10
12	close (3)
	read (aline,*) a,b,c,d,f0
	nl=notbnk(aline,90)
	nb=lblnk(aline,nl)
	write(6,13) a,b,c,d,f0,aline(nb:nl)
13	format(/10x,'Temperature calibration coefficients:'/5x,5e15.7/
     1	10x,'Calib. date =',a)
	fnam(ic:ic)='c'
	open(unit=3,file=fnam,type='old',err=200)
20	read (3,11,end=22) aline
	go to 20
22	close (3)
	read (aline,*) aa,bb,cc,dd,m
	nl=notbnk(aline,90)
	nb=lblnk(aline,nl)
	write(6,23) aa,bb,cc,dd,m,aline(nb:nl)
23	format(/10x,'Conductivity calibration coefficients:'/5x,4e15.7,
     1	f6.1/10x,'Calib. date =',a)
	if(iprs.eq.0) return
	fnam(ic:ic)='p'
	open(unit=3,file=fnam,type='old',err=300)
30	read (3,11,end=32) aline
	go to 30
32	close (3)
	read (aline,*) ap,bp,cp
	nl=notbnk(aline,90)
	nb=lblnk(aline,nl)
	write(6,33) ap,bp,cp,aline(nb:nl)
33	format(/10x,'Pressure calibration coefficients:'/5x,3e15.7,
     1	/10x,'Calib. date =',a)
	return
100	WRITE (6,101) ISN,fnam
101	FORMAT(//10X,'NO TEMP COEF. FOR SEACAT S/N ',I3/
     1	10x,'Tried to open file ',a)
  	STOP
200	WRITE (6,201) ISN,fnam
201	FORMAT(//10X,'NO COND. COEF. FOR SEACAT S/N ',I3/
     1	10x,'Tried to open file ',a)
  	STOP
300	WRITE (6,301) ISN,fnam
301	FORMAT(//10X,'NO PRESS. COEF. FOR SEACAT S/N ',I3/
     1	10x,'Tried to open file ',a)
  	STOP
	end
	function catprs(ihex)
c
c--converts pressure counts to psi.  Subtracts 14.7 psi (as per seasoft)
c--then converts psi to dbar
c
	common /pcoef/ ap(3)
	if(ihex.gt.4096) ihex=-1.0*(ihex.and.'FFF'X)
	cnt=ihex
	catprs=(poly(cnt,ap,2)-14.7)*.6894757
	return
	end
	FUNCTION CATEMP(F)
c
c--program to compute SEACAT temperature from frequency
c
	REAL*8 A,B,C,D,F0,FF
	common /tcoef/ a,b,c,d,f0
 	FF=F0/F
	T=1./(A+B*(DLOG(FF))+C*(DLOG(FF))**2 + 
     1  D*(DLOG(FF))**3) - 273.15
	CATEMP=T
	RETURN
	END
c	FUNCTION CATCON(F,T)
c
c        16 jul 97 -- changed to handle
c        ghij coefficients 
c
c--program to compute SEACAT conductivity from frequency
c
c	common /ccoef/ a,b,c,d,m
c	REAL*8 A,B,C,D
c	REAL*4 M
c	CATCON=(A*F**M+B*F**2+C+D*T)/10.
c	RETURN
c 200	STOP
c	END
	FUNCTION CATCON(F,T)
c
c--program to compute SEACAT conductivity from frequency
c--modified 16 Jul 1997 to be able to handle ghij coefficients.
c
	common /ccoef/ a,b,c,d,m
        parameter (ctcor=3.25e-06)
	REAL*8 A,B,C,D
	REAL*4 M
        if (m.gt.0.1)then
  	     CATCON=(A*F**M+B*F**2+C+D*T)/10.
        else
             CATCON = (A+B*F**2+C*F**3+D*F**4)/(10.*(1.0+ctcor*T))
        endif
	RETURN
  200	STOP
	END

      REAL FUNCTION SAL78 (CND,T,P,M)
C
C *************************************************************
C  THE CONDUCTIVITY RATIO (CND) = 1.0000000 FOR SALINITY = 35 PSS-78
C  TEMPERATURE = 15.0 DEG. CELSIUS, AND ATMOSPHERIC PRESSURE.
C
C  WHOI FUNCTION TO CONVERT CONDUCTIVITY RATIO TO SALINITY (M=0)
C  SALINITY TO CONDUCTIVITY RATIO (M=1 AND CND BECOMES INPUT SALINITY)
C
C--REFERENCES:  ALSO LOCATED IN UNESCO REPORT NO. 37 1981
C  PRACTICAL SALINITY SCALE 1978: E.L. LEWIS IEEE OCEAN ENG. JAN. 1980
C
C--UNITS:
C	PRESSURE	P	DECIBARS
C	TEMPERATURE	T	DEG CELSUIS (IPTS-68)
C	CONDUCTIVITY	CND	RATIO 	(M=0)
C	SALINITY	SAL78	(PSS-78) (M=0)
C--CHECKVALUES
C	SAL78=1.888091  : CND=40.0000, T=40 DEG C, P=10000 DB:  M=1
C	SAL78=40.00000	: CND=1.888091,T=40 DEG C, P=10000 DB:  M=0
C
C  SAL78 RATIO :RETURNS ZERO FOR CONDUCTIVITY RATIO: <0.0005
C  SAL78:RETURNS ZERO FOR SALINITY: < 0.02
C
C  PRACTICAL SALINITY SCALES 1978 DEFINITION WITH TEMPERATURE CORRECTION
C  XT=T-15.0 :  XR=SQRT(RT)
C
C  FROM FOFONOFF, N.P. AND R.C. MILLARD, JR., 1983: ALGORITHMS
C    	FOR OCEANOGRAPHIC COMPUTATIONS, DRAFT COPY
C
C *************************************************************
C
      SAL(XR,XT)=((((2.7081*XR-7.0261)*XR+14.0941)*XR+25.3851)*XR
     2-0.1692)*XR + 0.0080 + (XT/(1.0+0.0162*XT))*(((((-0.0144*XR
     3+0.0636)*XR-0.0375)*XR-0.0066)*XR-0.0056)*XR+0.0005)
C--DSAL(XR,XT) FUNCTION FOR DERIVATIVE OF SAL(XR,XT) WITH XR
      DSAL(XR,XT)=((((13.5405*XR-28.1044)*XR+42.2823)*XR+50.7702)*XR
     1-0.1692) + (XT/(1.0+0.0162*XT))*((((-0.0720*XR+0.2544)*XR
     2-0.1125)*XR-0.0132)*XR-0.0056)
C--FUNCTION RT35 : C(35,T,0)/C(35,15,0) VARIATION WITH TEMPERATURE
      RT35(XT)=(((1.0031E-9*XT-6.9698E-7)*XT+1.104259E-4)*XT
     1 + 2.00564E-2)*XT + 0.6766097
C--POLYNOMIALS OF RP: C(S,T,P)/C(S,T,0) VARIATION WITH PRESSURE
C--C(XP) POLYNOMIAL CORRESPONDS TO A1-A3 CONSTANTS: LEWIS 1980
      C(XP) = ((3.989E-15*XP-6.370E-10)*XP + 2.070E-5)*XP
      B(XT) = (4.464E-4*XT+3.426E-2)*XT + 1.0
C--A(XT) POLYNOMIAL CORRESPONDS TO B3 AND B4 CONSTANTS: LEWIS 1980
      A(XT) = -3.107E-3*XT + 0.4215
C--ZERO SALINITY/CONDUCTIVITY TRAP
      SAL78 = 0.0
      IF ((M.EQ.0).AND.(CND.LE.5E-4)) RETURN
      IF ((M.EQ.1).AND.(CND.LE.0.02)) RETURN
C--TEMP ADJUSTMENT
      DT = T - 15.0
C--SELECT BRANCH FOR SALINITY (M=0) OR CONDUCTIVITY (M=1)
      IF (M.EQ.1) GO TO 10
C
C--CONVERT CONDUCTIVITY TO SALINITY
C
      R = CND
      RT = R/(RT35(T)*(1.0 + C(P)/(B(T) + A(T)*R)))
      RT = SQRT (ABS(RT))
      SAL78 = SAL(RT,DT)
      RETURN
C
C--CONVERT SALINITY TO CONDUCTIVITY BY THE 
C--NEWTON-RAPHSON ITERATIVE METHOD
C
C--FIRST APPROXIMATION
   10 RT = SQRT(CND/35.0)
      SI = SAL(RT,DT)
      N = 0
C--ITERATION LOOP BEGINS HERE WITH A MAXIMUM OF 10 CYCLES
   15 RT = RT + (CND - SI)/DSAL(RT,DT)
      SI = SAL(RT,DT)
      N = N + 1
      DELS = ABS(SI-CND)
      IF ((DELS.GT.1.0E-4).AND.(N.LT.10)) GO TO 15
C--COMPUTE CONDUCTIVITY RATIO
      RTT = RT35(T)*RT*RT
      AT = A(T)
      BT = B(T)
      CP = C(P)
      CP = RTT*(CP + BT)
      BT = BT - RTT*AT
C
C--SOLVE QUADRATIC EQUATION FOR R: R = RT35*RT*(1+C/AR+B)
C
      R = SQRT(ABS(BT*BT + 4.0*AT*CP)) - BT
C--CONDUCTIVITY RETURN
      SAL78 = 0.5*R/AT
      RETURN
      END