	PROGRAM scEDIT2
C
c	version 2: computes and outputs sigma-theta as well as t and s
c	and works for units with pressure
c
c
C	EDITS SeaCat FILES FROM SCREAD.
C	OUTPUTS LINES NSTART THROUGH NEND AND INTERPOLATES
C	DATA OUTSIDE RANGE TMIN TO TMAX, SMIN to SMAX OR LINES WHICH ARE
C	FLAGGED IN ARRAY NB.
C	ALSO HAS OPTION TO ADD IA LINES BEFORE LINE NUMBER NA
C
C	WRITTEN BY PAUL FREITAG, JUNE 88
C
c	not tested for nbad or nadd > 0
c
	character*24 pdate
	character*40 form
	character*50 head
	character*80 fnam1,fnam2
	DIMENSION HOLD(200,3)
	DIMENSION NB(200),NA(25),IA(25)
	integer*2 itime(5,200)
C
	nrun=0
100	read (5,*,end=999,err=999) FNAM1,fnam2,nstart,nend,nbad,nadd
	call fOPEN(1,FNAM1,'OLD')
	READ (1,102) HEAD
102	format(a)
	nrun=nrun+1
	if(nrun.gt.1) call pej
	WRITE (6,111) FNAM1,HEAD
111	FORMAT(' PROGRAM SCEDIT2'//10X,'INPUT FILE NAME: ',A/
     1  10X,'INPUT FILE HEADER: ',A)
	call fOPEN(2,FNAM2,'NEW')
	read(head(25:30),*) dep
	read(head(17:18),'(i2)') iv
	iprs=0
	if(iv.lt.0) iprs=1
	head(31:40)=' SCEDIT2'
	CALL fDATE(pdate)
	head(40:50)=pdate(5:10)//pdate(20:24)
	WRITE (6,127) FNAM2,HEAD
127	FORMAT(//10X,'OUTPUT FILE NAME: 'A,/
     1  10X,'OUTPUT FILE HEADER: 'A/10x,'Run date: ',A)
	WRITE (2,102) HEAD
	WRITE (6,132) NSTART,NEND,NBAD,NADD
132	FORMAT(//10X,'NSTART  NEND  NBAD  NADD'/10X,4I6)
	write(6,133) dep
133	format(10x,'Depth = ',f8.1)
	IF(NBAD.LE.0) GO TO 135
	IF(NBAD.GT.200) NBAD=200
	read (5,*) (NB(I),I=1,NBAD)
	WRITE (6,134) (NB(I),I=1,NBAD)
134	FORMAT(//10X,'FOLLOWING LINES INTERPOLATED'//
     1  (10X,10I6))
135	IF(NADD.LE.0) GO TO 140
	IF(NADD.GT.25) NADD=25
	read (5,*) (NA(I),IA(I),I=1,NADD)
	WRITE (6,139) (IA(I),NA(I),I=1,NADD)
139	FORMAT(10X,I5,' LINES ADDED BEFORE LINE ',I5)
140	WRITE (6,141) 
141	FORMAT(//10X,'VALUES OUTSIDE LIMITS INTERPOLATED'//
     1  8X,'VARIABLE MINIMUM MAXIMUM')
	if(iprs.eq.0) then
	read (5,*) TMIN,TMAX,SMIN,SMAX
	WRITE (6,151) TMIN,TMAX,SMIN,SMAX
151	FORMAT(6X,'Temperature',2F8.3/8x,'Salinity',2f8.3
     1	//' TIME DDMMYY   Temp.  Salin.'//)
	else
	read (5,*) TMIN,TMAX,SMIN,SMAX,pmin,pmax
	WRITE (6,1151) TMIN,TMAX,SMIN,SMAX,pmin,pmax
1151	FORMAT(6X,'Temperature',2F8.3/8x,'Salinity',2f8.3/
     1	8x,'Pressure',2f8.3//' TIME DDMMYY   Temp.  Salin.   Pres.'//)
	endif
	NDATA=NEND-NSTART+1
	NVAR=2+iprs
	write(form,191) nvar+1
191	format('(1x,2i2,1x,3i2,',i1,'f8.3,i6)')
	DO 190 N=1,NSTART
190	READ(1,form) (ITIME(I,1),I=1,5),(HOLD(1,I),I=1,NVAR)
	NB1=1
	NA1=1
	NR=NSTART
	NINT=0
	NOUT=1
	t0=theta(hold(1,2),hold(1,1),dep,0.0)
	spvol= SVAN(hold(1,2),t0,0.0,SIGMA)
	WRITE (2,form) (ITIME(I,1),I=1,5),(HOLD(1,I),I=1,NVAR),
     1	sigma,NOUT
	WRITE (6,form) (ITIME(I,1),I=1,5),(HOLD(1,I),I=1,NVAR),
     1	sigma,NOUT
	NH=1
	ND9=NDATA-9
	DO 300 NN=2,NDATA
	NH=NH+1
	READ (1,form,END=310) (ITIME(I,NH),I=1,5),
     1  (HOLD(NH,I),I=1,NVAR)
	NR=NR+1
C
C	CHECK AND INTERPOLATE SCANS IF REQUESTED
C
201	IF(NB1.GT.NBAD) GO TO 212
	IF(NR-NB(NB1)) 212,202,211
202	HOLD(NH,1)=TMAX+10.
	HOLD(NH,2)=SMAX+10.
	if(iprs.ne.0) HOLD(NH,3)=pMAX+10.
	NB1=NB1+1
	GO TO 300
211	NB1=NB1+1
	GO TO 201
C
C	CHECK AND ADD SCANS IF REQUESTED
C
212	IF(NA1.GT.NADD) GO TO 220
	NFLAG=0
	IF(NR-NA(NA1)) 220,213,219
213	NFLAG=1
	IADD=IA(NA1)
	ND9=ND9+IADD
	CALL ADDSCN(ITIME,HOLD,IADD,NOUT,NH,DELT,NVAR,FORM)
219	NA1=NA1+1
	IF(NFLAG.EQ.0) GO TO 212
220	T=HOLD(NH,1)
	S=HOLD(NH,2)
	if(iprs.ne.0) p=HOLD(NH,3)
	IF(T.LT.TMIN.OR.T.GT.TMAX) GO TO 300
	IF(S.LT.SMIN.OR.S.GT.SMAX) GO TO 300
	IF(iprs.ne.0.and.(p.LT.pMIN.OR.p.GT.pMAX)) GO TO 300
	IF(NH.EQ.2) GO TO 250
	CALL DESP(HOLD(1,1),NH,TMIN,TMAX)
	CALL DESP(HOLD(1,2),NH,SMIN,SMAX)
	if(iprs.ne.0) CALL DESP(HOLD(1,3),NH,pMIN,pMAX)
250	DO 280 N=2,NH
	t0=theta(hold(n,2),hold(n,1),dep,0.0)
	spvol= SVAN(hold(n,2),t0,0.0,SIGMA)
	NOUT=NOUT+1
	IF(NOUT.LE.10.OR.NOUT.GE.ND9)
     1  WRITE (6,form) (ITIME(I,N),I=1,5),(HOLD(N,I),I=1,NVAR),
     1	sigma,NOUT
	IF(NH.EQ.2) GO TO 280
	if(n.ne.nh) NINT=NINT+1
	IF(N.EQ.2.OR.N.EQ.NH-1)THEN
	WRITE (6,form) (ITIME(I,N),I=1,5),(HOLD(N,I),I=1,NVAR),
     1	sigma,NOUT
	WRITE (6,261) N-1
261	FORMAT('+',34X,'**',I6)
	ENDIF
280	WRITE (2,form) (ITIME(I,N),I=1,5),(HOLD(N,I),I=1,NVAR),
     1	sigma,NOUT
	DO 290 I=1,5
290	ITIME(I,1)=ITIME(I,NH)
	DO 295 I=1,NVAR
295	HOLD(1,I)=HOLD(NH,I)
	NH=1
300	CONTINUE
	GO TO 320
310	WRITE (6,311)
311	FORMAT(//10X,'EOF ON UNIT 1 BEFORE REACHING LINE NEND')
320	CLOSE(UNIT=1)
	CLOSE(UNIT=2)
	PER=FLOAT(NINT)/NOUT*100.
	WRITE (6,323) NINT,PER
323	FORMAT(//10X,I5,' LINES INTERPOLATED - ',F5.1,' PERCENT')
	GO TO 100
999	STOP
	END
	SUBROUTINE ADDSCN(ITIME,HOLD,IADD,NOUT,NH,DELT,NSENS,FORM)
	character*40 form
	DIMENSION TT(2),TX(2)
	DIMENSION ITIME(5,200),HOLD(200,2),ntime(5)
	FADD=IADD+1
	N=NH-1
	DO 214 K=1,NSENS
	TT(K)=HOLD(N,K)
214	TX(K)=(HOLD(NH,K)-TT(K))/FADD
	DO 218 I=1,IADD
	CALL GETDT2(ITIME,1,DELT,I+1,ntime)
	NOUT=NOUT+1
	RR=I
	DO 217 K=1,NSENS
217	TT(K)=TT(K)+TX(K)
	IF(I.EQ.1.OR.I.EQ.IADD) WRITE (6,FORM) ntime,(TT(K),K
     1	=1,NSENS),RR,NOUT
218	WRITE (2,FORM) ntime,(TT(K),K=1,NSENS),RR,NOUT
	RETURN
	END
	SUBROUTINE DESP(X,NX,XMIN,XMAX)
	DIMENSION X(1)
	NB=0
	X1=X(1)
	DO 50 I=2,NX
5	IF(X(I).GE.XMIN.AND.X(I).LE.XMAX) GO TO 20
	NB=NB+1
	GO TO 50
20	IF(NB.EQ.0) GO TO 50
	FNB=NB+1
	XD=(X(I)-X1)/FNB
	DO 30 J=1,NB
	J1=I-NB-1+J
30	X(J1)=X(J1-1)+XD
	NB=0
	X1=X(I)
50	CONTINUE
	RETURN
	END
	REAL FUNCTION SVAN (S,T,PO,SIGMA)
C
C**************************************************
C 
C SPECIFIC VOLUME ANOMALY (STERIC ANOMALY) BASED ON 1980 EQUATION
C OF STATE FOR SEAWATER AND 1978 PRACTICAL SALINITY SCALE.
C
C--REFERENCES:	MILLERO,ET AL (1980) DEEP-SEA RES.,27A,255-264.
C 		MILLERO AND POISSON 1981,DEEP-SEA RES.,28A,625-629.
C  BOTH ABOVE REVERENCES ARE ALSO FOUND IN UNESCO REPORT NO. 38 (1981)
C
C--UNITS:
C	PRESSURE	PO	DECIBARS
C	TEMPERATURE	T	DEG CELSIUS (IPTS-68)
C	SALINITY 	S	(PSS-78)
C	SPEC VOL ANOM	SVAN	1.0E-8M**3/KG
C	DENSITY ANOM	SIGMA	KG/M**3
C
C--CHECK VALUE:
C	SVAN = 981.30210 E-8 M**3/KG FOR S=40 (PSS-78),
C		T = 40 DEG C, PO = 10000 DECIBARS
C	SIGMA = 59.82037 KG/M**3 FOR S = 40 (PSS-78),
C		T = 40 DEG C, PO = 10000 DECIBARS
C
C--FROM FOFONOFF,N. AND R.C.MILLARD,1983,DRAFT COPY. MAY 84,LJM.
C
C***********************************************************
C
	REAL P,T,S,SIG,SR,R1,R2,R3,R4
   	REAL A,B,C,D,E,A1,B1,AW,BW,K,K0,KW,K35
C--EQUIVALENCE
	EQUIVALENCE (E,D,B1),(BW,B,R3),(C,A1,R2)
	EQUIVALENCE (AW,A,R1),(KW,K0,K)
C--DATA
	DATA R3500,R4/1028.1063,4.8314E-4/
	DATA DR350/28.106331/
C--R4 IS REFERRED TO AS C IN MILLERO AND POISSON 1981
C--CONVERT PRESSURE TO BARS AND TAKE SQUARE ROOT OF SALINITY
	P = PO/10.
  	SR = SQRT(ABS(S))
C--PURE WATER DENSITY AT ATMOSPHERIC PRESSURE
C--BIGG P.H.(1967) BR. J. APPLIED PHYSICS 8 PP 521-537.
	R1 = ((((6.536332E-9*T-1.120083E-6)*T+1.001685E-4)*T
     1 - 9.095290E-3)*T+6.793952E-2)*T-28.263737
C--SEAWATER DENSITY ATM PRESS
C--COEFFICIENTS INVOLVING SALINITY
C--R2 = A  IN NOTATION OF MILLERO AND POISSON 1981
	R2 = (((5.3875E-9*T-8.2467E-7)*T+7.6438E-5)*T-4.0899E-3)*T
     1 + 8.24493E-1
C--R3 = B IN NOTATION OF MILLERO AND POISSON 1981
	R3 = (-1.6546E-6*T+1.0227E-4)*T-5.72466E-3
C--INTERNATIONAL ONE-ATMOSPHERE EQUATION OF STATE OF SEAWATER
	SIG = (R4*S + R3*SR + R2)*S + R1
C--SPECIFIC VOLUME AT ATMOSPHERIC PRESSURE
	V350P = 1.0/R3500
	SVA = -SIG*V350P/(R3500+SIG)
	SIGMA = SIG + DR350
C--SCALE SPECIFIC VOL ANOMALY TO NORMALLY REPORTED UNITS
	SVAN = SVA*1.0E+8
	IF (P.EQ.0.0) RETURN
C
C--NEW HIGH PRESSURE EQUATION OF STATE FOR SEAWATER
C
C--MILLERO, ET AL, 1980 DSR 27A, PP 255-264
C--CONSTANT NOTATION FOLLOWS ARTICLE
C
C--COMPUTE COMPRESSION TERMS
	E = (9.1697E-10*T+2.0816E-8)*T-9.9348E-7
   	BW = (5.2787E-8*T-6.12293E-6)*T+3.47718E-5
	B = BW + E*S
C
	D = 1.91075E-4
	C = (-1.6078E-6*T-1.0981E-5)*T+2.2838E-3
	AW = ((-5.77905E-7*T+1.16092E-4)*T + 1.43713E-3)*T
     1 - 0.1194975
	A = (D*SR + C)*S + AW
C
	B1 = (-5.3009E-4*T+1.6483E-2)*T+7.944E-2
	A1 = ((-6.1670E-5*T+1.09987E-2)*T-0.603459)*T+54.6746
	KW = (((-5.155288E-5*T+1.360477E-2)*T-2.327105)*T
     1 + 148.4206)*T - 1930.06
	K0 = (B1*SR + A1)*S + KW
C--EVALUATE PRESSURE POLYNOMIAL WHERE:
C	K EQUALS SECANT BULT MODULUS OF SEAWATER
C	DK=K(S,T,P)-K(35,0,P)
C	K35=K(35,0,P)
C
	DK = (B*P + A)*P + K0
	K35 = (5.03217E-5*P+3.359406)*P + 21582.27
	GAM = P/K35
	PK = 1.0 - GAM
	SVA = SVA*PK + (V350P+SVA)*P*DK/(K35*(K35+DK))
C--SCALE SPECIFIC VOLUME ANOMALY TO NORMALLY REPORTED UNITS
	SVAN = SVA*1.0E+8
	V350P = V350P*PK
C
C--COMPUTE DENSITY ANOMALY WITH RESPECT TO 1000.0 KG/M**3
C	1) DR350: DENSITY ANOMALY AT 35 PSS, 0 DEG C AND 0 DBARS
C	2) DR35P: DENSITY ANOMALY AT 35 PSS, 0 C, PRES. VARIATION
C	3) DVAN:  DENSITY ANOMALY VARIATIONS INVOLVING SPECIFIC
C		  VOLUME ANOMALY
C
C--CHECK VALUE: SIGMA = 59.82037 KG/M**3 FOR S = 40 PSS,
C		T = 40 DEG C, PO = 10000 DB
C
	DR35P = GAM/V350P
	DVAN = SVA/(V350P*(V350P+SVA))
	SIGMA = DR350 + DR35P - DVAN
	RETURN
	END
