c
c
c "OSadcp.f"  processes RDI-LOGADCP ascii data
c             (earlier versions were 	test.f
c					9204.f
c					adcp.f
c					bbadcp.f
c
c                 x is latitude north
c                 y is longitude east
c                 t is in minutes
c
c input files:  #7  "xxxx.files", where xxxx is the experiment name, i.e., 9204
c               #9  Navigation files
c               #10 ADCP .raw data files
c    
c output files #3  adcp.csh for sending plotxy plots to printer
c               #4  adcp.log for monitoring the parameters of the run
c
c               #24 transect data (by file)
c               #25 tilt data (by file) no screen, no average
c               #26 ship velocity data (by file)
c               #27 shp velocity data (all)
c               #28 velocity data (by file)
c            #36 calibration data (by file)
c
c plotxy files  "pltvs.*"  ship speed from BT and GPS
c               "pltxy.*"  tilt sensors
c
c 09-14-1992 amue (9204.f)
c 03-04-1993 amue (modified for U of D use, test.f)
c 04-22-1993 amue
c 06-06-1993 amue (rotations after averaging, directional averaging of tilt)
c 11-07-1994 amue (modified for BB-adcp data at Rutgers University
c 12-24-1996 amue (modified code to accomodate 20 deg. angles;
c                                to screen for excessive tilt angles;
c                                to improve "passing" statistics
c 03-17-1997 amue (error in BT velocity computation removed)
c 09-24-1998 amue new time base, jan.1=day 0
c 08-12-2003 amue major overhaul for Ocean Surveyor data
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	common /raw/line1(6),line2(17),line3(5),line4(8),line5(13),line6(5),
     1			line(85,18),id(85),d1(85),d2(85),d3(85),d4(85),
     1			bt1,bt2,bt3,bt4,h1,h2,h3,h4,ncol
	common /ini0/navg,nensemble,ivel,i_old,ntot,ical,nbt,nxy,
     1	            cheadbar,sheadbar,cpbar,spbar,crbar,srbar
	common /ini1/xbar,ybar,btxbar,btybar,btzbar,hbar,tbar,hsdbar,xship,yship
	common /ini2/c30,s30,c,fact,fact1,xsc_lat,pi,d2r,idum
c
	parameter (nz=85)
	parameter (nn=99000)
	integer line1,line2,line3,line4,line5,line6
	integer line,id,mon(12),mon_n(12),mon_l(12),day0,zoff
	integer iq(nn),ibot(4),ibot2(4),nzmax1,ntot,nzmax,iww2
	real ka,ka0,m(4),sc(4),d1,d2,d3,d4,tilt(8),df(4),dcor(4),xcor(4),xint(4)
	real c30,s30,c,fact,fact1,pi,d2r,sig_up,xship,yship
	real h_nav(nn),p_nav(nn),r_nav(nn),hadcp,havg,h(1000)
	character file1*9,file2*12,label*3,path*23,adcp(9)*49,nav*49,all1*93,nav2*49 
	character name*36,input*12,cal_name*15,all*91,name1*9,zname*9,btname*9
	character biname*9,trans*98
	integer nensemble,bt1,bt2,bt3,bt4,h1,h2,h3,h4,nbar(nz),cdate,sec
	real lat,lon,zsg(4),head,headbar1,headbar,hsdev,hsdbar,hdif,hvar,heading(1000)
	real xx(nn),yy(nn),tt(nn),uu(nn),vv(nn),tstick(nn),hsdev_crit
	real gyro(nn),ww(nz),cbar(nz),sbar(nz),cor_min(nz)
	real aspeed,adir,gpx(nn),gpy(nn),ugps(1000),vgps(1000)
	real h_seab(100,200),hh,lon0,lat0,dlon,dlat
	real tsca1,tsca2,tsca3,tsca4,dt
	real sxx(1000),sxvar(nz,1000),syvar(nz,1000),szvar(nz,1000)
	real vsx(1000),vsy(1000),vsz(1000)
	double precision tnav1,tnav2,dt_nav,dx_nav,dy_nav,gpt(nn)
	double precision xsc_lat,vsx1,vsy1,dtt,dx,dy
	double precision sheadbar,cheadbar,spbar,srbar,cpbar
	double precision xbar,ybar,btxbar,btybar,btzbar,bt1bar,bt2bar,bt3bar,bt4bar
	double precision ubar(nz),vbar(nz),wbar(nz),ebar(nz)
c
c	data (mon(i),i=1,11)/31,28,31,30,31,30,31,31,30,31,30/
c leap year, 1/1 00:01 is time=0.0
	data (mon_n(i),i=1,12)/-1,30,58,89,119,150,180,211,242,272,303,333/
	data (mon_l(i),i=1,12)/-1,30,59,90,120,151,181,212,243,273,304,334/
c 
	data z_adcp/15./,ipg_crit/25/
	data ipg_bt_crit/25/,bte_crit/50./,ve_crit/50./,hsdev_crit/8./
	open (unit=4,file='adcp.log')
	open (unit=3,file='adcp.csh')
	write(3,*)'#! /bin/csh'
c
	call attention
c
	pi = 4.*atan(1.)
c
c vertical misalignment angle for beam #1+2 (arg1) and beam #3+4 (arg2) 
c
c	Healy Sept.-2003
	alpha2 =  0.94
	alpha3 = -0.90
c Healy July/August 2003
	alpha2 = 0.76
	alpha3 = -1.10
c
	arg1 = alpha2*pi/180.
	arg2 = alpha3*pi/180.
	blanking = 10.
	e_urick = 0.025
	db_scale = 0.45
c
	ct_crit = cos(pi*15./180.)
	ct_crit = cos(pi*10./180.)
	tsca1 = 0
	tsca2 = 24*60*60
	tsca3 = 60*60
	tsca4 = 60
c
	write(6,*)
	write(6,*)'Screening Variables:'
	write(6,*)
	write(6,*)'Percent good water  pings (%)       : ',ipg_crit
	write(6,*)'Percent good BT pings (%)           : ',ipg_bt_crit
	write(6,*)'BT Error Velocity (cm/s)            : ',bte_crit
	write(6,*)'Water Error Velocity (cm/s)         : ',ve_crit
	write(6,*)
c
	    write(4,*)'Screening Variables:'
	    write(4,*)
	    write(4,*)'Percent good water  pings (%)       : ',ipg_crit
	    write(4,*)'Percent good BT pings (%)           : ',ipg_bt_crit
	    write(4,*)'BT Error Velocity (cm/s)            : ',bte_crit
	    write(4,*)'Water Error Velocity (cm/s)         : ',ve_crit
	    write(4,*)
c
	do 18 i=1,12
		mon(i) = mon_n(i)
18	continue
c
c averaging intervall dt in minutes
c
	    write(6,*)' enter averaging intervall (min)'
	    read(5,*)dt
c		dt = 5.
	    write(4,143)dt
143	format(/'averaging intervall (min) ',f10.1/)
c
c enter filename that contains the list of file to be procesed
c
		input = '0307.files'
c	write(6,*)'Enter input file'
c	read(5,*)input
112	format(a12)
	open (unit=7,file=input)
c
c 	read input parameters common to all files
c
c mon0	 month of time reference
c day0	 day of time reference
c nfile	 # of adcp files to be processed
c path	 directory where the adcp and navigation file are to be found
c iplt	 plotting on (1) or off (0)
c
	    read(7,102)mon0,day0,nfile,path,zoff,iplt,ical0,
     1				cal_a,cal_b,cal_c,cal_d,cal_e,a
	    write(6,102)mon0,day0,nfile,path,zoff,iplt,ical0,
     1				cal_a,cal_b,cal_c,cal_d,cal_e,a
	    write(4,*)
	    write(4,102)mon0,day0,nfile,path,zoff,iplt,ical0,
     1				cal_a,cal_b,cal_c,cal_d,cal_e,a
	    write(4,*)
c
c process one file at a time, total of nfile files
c
	    write(4,*)'=====> start processing ADCP data <======= '
	    write(4,*)
	    do 1 ifile=1,nfile
		izbar = 0
		n76 = 0
	 	call initial(ifile)
		write(6,*)ifile,xbar,ybar
c
c	read input parameter from file
c
	    read(7,103)file1,nens,label,ncol,ictd,izcut,
     1		           c1head,s1head,c2head,s2head,alpha,beta
	    if (ical0.eq.1) then
		    alpha = cal_a
		    c1head = cal_b
		    s1head = cal_c
		    c2head = cal_d
		    s2head = cal_e
		    beta = a
	    else if (ical0.eq.0) then 
		    c1head = 0.
		    s1head = 0.
		    c2head = 0.
		    s2head = 0.
		    alpha = 0.
		    beta = 1.
		    alpha2 = 0.
		    alpha3 = 0.
		    arg1 = 0.
		    arg2 = 0.
	    end if
	    write(4,140)ical0,c1head,s1head,c2head,s2head,alpha,alpha2,alpha3,beta
	    write(4,141)file1,nens,ncol,ictd,izcut,zoff,c
	    write(6,141)file1,nens,ncol,ictd,izcut,zoff,c
140	format(	'Calibration mode (0-none,1-uniform,2-individual ',i12/
     1		'Coefficient for cos(heading)                    ',f12.1/
     1		'Coefficient for sin(heading)                    ',f12.1/
     1		'Coefficient for cos(2*heading)                  ',f12.1/
     1		'Coefficient for sin(2*heading)                  ',f12.1/
     1		'Misalignment (x-y plane) alpha-1                ',f12.2/
     1		'Misalignment (x-z plane) alpha-2                ',f12.2/
     1		'Misalignment (y-z plane) alpha-3                ',f12.2/
     1		'Scaling (x-y plane only) beta                   ',f12.5/)
141	format(	'file1--> adcp data file                         ',a12/
     1		'nens --> #of ensembles in the adcp and nav file ',i12/
     1		'ncol --> # of columns in the ADCP data file     ',i12/
     1		'ictd --> format of the navigation data file     ',i12/
     1		'izcut--> # of adcp bins cut prior to processing ',i12/
     1		'zoff -->  depth of ADCP transducers             ',i12/
     1		'c    -->  speed of sound (m/s)                  ',f12.0)
c
c get full filename (with path) of data files in a character array
c
 	    open(unit=10,file='file.name')
 	    write(10,127)path,file1
127	format(a23,'/ascii/',a9,'.l1')
 	    write(10,128)path,file1
128	format(a23,'/ascii/',a9,'.l2')
 	    write(10,129)path,file1
129	format(a23,'/ascii/',a9,'.l3')
 	    write(10,130)path,file1
130	format(a23,'/ascii/',a9,'.l4')
 	    write(10,131)path,file1
131	format(a23,'/ascii/',a9,'.line1')
 	    write(10,132)path,file1
132	format(a23,'/ascii/',a9,'.line2')
 	    write(10,134)path,file1
134	format(a23,'/ascii/',a9,'.line4')
 	    write(10,135)path,file1
135	format(a23,'/ascii/',a9,'.line5')
 	    write(10,136)path,file1
136	format(a23,'/ascii/',a9,'.line6')
	if (input.eq.'0307.files') then
 	    write(10,144)path,file1
	else
 	    write(10,137)path,file1
	end if
137	format(a23,'/ascii/',a9,'.nav')
144	format(a23,'/ascii/',a9,'.n1r')
 	    write(10,239)path,file1
239	format(a23,'/ascii/',a9,'.n2r')
	    close(10)
	open(unit=10,file='file.name')
	read(10,108)(adcp(i),i=1,9)
	read(10,108)nav
	read(10,108)nav2
	close(10)
c
	open(unit=88,file=adcp(4))
	open(unit=89,file=adcp(3))
	open(unit=91,file=adcp(5))
	open(unit=92,file=adcp(6))
	open(unit=94,file=adcp(7))
	open(unit=95,file=adcp(8))
	open(unit=96,file=adcp(9))
	open(unit=97,file=adcp(1))
c Cor:
	open(unit=98,file=adcp(2))
	open(unit=9,file=nav)
	open(unit=11,file=nav2)
c
c water depth from SeaBeam:
c
	open(unit=87,file='/Users/a2/0309_pro/seabeam/grid/seabeam.avg')
	read(87,*)nseabeam,lon0,lat0,dlon,dlat
	do 87 i=1,nseabeam
		read(87,187)hh,ix,iy
187		format(18x,f8.1,16x,2i4)
		h_seab(ix,iy) = -hh
87	continue
	close (87)
c
c nav file from .ENS file (aft P-code GPS on Healy 2003
c
	do 125 i=1,nens
		read(9,200)iensnav,imon,iday,ihr,imin,sec,msec,h_nav(i),xx(i),yy(i),ispd,dir
		if (abs(xx(i)).lt.15 .or. abs(yy(i)).lt.10) then
			xx(i) = xx(i-1)
			yy(i) = yy(i-1)
		end if
		read(11,201)iensnav2,h_nav(i),p_nav(i),r_nav(i),iq(i)
201		format(i9,40x,3f8.2,22x,i2)
200		format(i9,6i3,f8.2,2f14.8,i6,f8.2)
c
c units are hours
c
		
		gpt(i) = mon(imon)*24.+iday*24.+ihr+float(imin)/60.+float(sec)/3600.+float(msec)/360000.-6096.
125	continue
	told = gpt(1)*60+5
c
c Lanczos raised cosine filter with cut-off at 30 pings and 151 pings filter weights
c
	riww = 71
	tco = 15
	iww2 = (riww-1)/2+6
	call filter(1.,nens,xx,riww,tco)
	call filter(1.,nens,yy,riww,tco)
	call filter(1.,nens,h_nav,riww,tco)
	call filter(1.,nens,p_nav,riww,tco)
	call filter(1.,nens,r_nav,riww,tco)
	do 124 i=3,nens-2
		dt_nav = gpt(i+2) - gpt(i-2)
		arg = pi*yy(i)/180.
c
c units of GPS derived ship's velocities are m/s
c
		gpx(i) = (xx(i+2)-xx(i-2))*1852./60.*cos(arg)/dt_nav
		gpy(i) = (yy(i+2)-yy(i-2))*1852./60./dt_nav
124	continue
c	call filter(1.,nens,gpx,riww,tco)
c	call filter(1.,nens,gpy,riww,tco)
c
	close (9)
	close (11)
	    write(6,*)
	    write(6,*)'files open'
	    write(6,*)
	    open(unit=26,file='vship.dat')
	    open(unit=25,file='tilt.dat')
	    open(unit=28,file='velocity.dat')
	    open(unit=24,file='transect.dat')
c	    open(unit=36,file='calibration.dat')
c
	    nens = nens-2
	    nens_avg = 0
	    call read_os(iens,ncell,nz)
c ADCP time (in seconds)
	    t1 = mon(line1(2))
	    t2 = (t1+line1(3))*tsca2
	    t3 = line1(4)*tsca3
	    t4 = line1(5)*tsca4
	    tt(1) = t1+t2+t3+t4+line1(6)
c
c scale to convert counts to cm/s
c
	call co_sys(c,scale,line2,zsg,sig_up)
		vxscale = scale/2./s30
		vzscale = scale/4./c30
c
	sca1228 = 1
	if (line2(17).eq.307 .or. line2(17).eq.300) then
		write(6,*)
		write(6,*)'Assuming 307 kHz ADCP'
		write(6,*)
		z_crit = 80
	else if (line2(17).eq.153) then
		write(6,*)
		write(6,*)'Assuming 153 kHz ADCP'
		write(6,*)
		z_crit = 80
	else if (line2(17).eq.75) then
		write(6,*)
		write(6,*)'Assuming 75 kHz ADCP'
		write(6,*)
		sca1228 = 0.1
		z_crit = 80
	else if (line2(17).eq.1228) then
		write(6,*)
		write(6,*)'Assuming 1228 kHz ADCP'
		write(6,*)
		z_crit = 3000
		sca1228 = 0.1
	else
		write(6,*)'need z_crit'
		stop
	end if
	write(4,*)'z_crit: ',z_crit
	write(4,*)'sca1228:',sca1228
c		
c time in minutes
c
	    time = tt(1)/60.
	    lat = yy(1)
	    lon = xx(1)
	    time_old = time
c
	    do 13 iz=1,nz
		    ubar(iz) = 0.
		    vbar(iz) = 0.
		    wbar(iz) = 0.
		    ebar(iz) = 0.
		    cbar(iz) = 0.
		    sbar(iz) = 0.
c		    cor_min(iz) = 999.
		    nbar(iz) = 0
13	continue
	ndata = 0
	nprof = 0
c
c =============>START PROCESSING OF DATA, LOOP #2 GOES OVER ALL ENEMBLES <=========
c
	do 23 i=2,iww2
		call read_os(iens,ncell,nz)
c ADCP time in seconds 
c
c month
		    t1 = mon(line1(2))
c day
		    t2 = (t1+line1(3))*tsca2
c hour
		    t3 = line1(4)*tsca3
c min
		    t4 = line1(5)*tsca4
		    tt(i) = t1+t2+t3+t4+line1(6)
		    time = tt(i)/60
23	continue
	time_old = time
	    do 2 i=iww2+1,nens-iww2
		call read_os(iens,ncell,nz)
		if (i.gt.iens) then
			iens = i
		end if
c ADCP time in seconds 
		    t1 = mon(line1(2))
		    t2 = (t1+line1(3))*tsca2
		    t3 = line1(4)*tsca3
		    t4 = line1(5)*tsca4
		    tt(i) = t1+t2+t3+t4+line1(6)
c time in minutes
		time = tt(i)/60
		lat = yy(i)
		lon = xx(i)
c
c convert tilt data to physical units
c correct for magnetic deviation
c
c ADCP gyro:
	    head = line4(4)*0.01 
	    hsdev = line4(7)*0.01
	    pitch = line4(2)*0.01
	    psdev = line4(5)*0.01
	    roll = line4(3)*0.01
	    rsdev = line4(6)*0.01
	    temp = line4(8)*0.01
c Ashtech:
	if (iq(i) .eq. 0) then
	    	head = h_nav(i)
 		pitch = p_nav(i)
 		roll = r_nav(i)
	end if
c
	    if (pitch.gt.20.) pitch = pitch-360.
	    if (roll.gt.20.) roll = roll-360.
c
c compute bin mapping coefficients
c
	    rr = roll*d2r
	    pp = pitch*d2r
c	    ka = sqrt( 1.-sin(pp)*sin(pp)*sin(rr)*sin(rr) )
c	    pp = asin( sin(pp)*cos(rr)/ka )
c
	    cp = cos(pp)
	    sp = sin(pp)
	    cr = cos(rr)
	    sr = sin(rr)
c
c discard data with excessive tilt or roll (typically ct_crit=cos(20 deg.)
c
	if (cp.lt.ct_crit .or. cr.lt.ct_crit) then
		write(6,*)'Bad attitude: ',i,pitch,roll,head
		goto 2
	end if
c
c tilt sensors fixed to ADCP
c
	    m(1) = -sr*cp *s30
	    m(2) = sp     *s30
	    m(3) = cp*cr  *c30
c
	    sc(1) = c30 / (m(3)+zsg(1)*m(1))
	    sc(2) = c30 / (m(3)+zsg(2)*m(1))
	    sc(3) = c30 / (m(3)+zsg(3)*m(2))
	    sc(4) = c30 / (m(3)+zsg(4)*m(2))

c Vs-BT
c
c 3-beam solutions
c
		btflag = 0
	    if (bt1.eq.99999 .and. max(abs(bt2),abs(bt3),abs(bt4)).lt.9999) then
		bt1 = bt3+bt4-bt2 
		h1 = h2
		line5(10) = 255
		btflag = 1
	    else if (bt2.eq.99999 .and. max(abs(bt1),abs(bt3),abs(bt4)).lt.9999) then
		bt2 = bt3+bt4-bt1
		h2 = h1
		line5(11) = 255
		btflag = 1
	    else if (bt3.eq.99999 .and. max(abs(bt1),abs(bt2),abs(bt4)).lt.9999) then
		bt3 = bt1+bt2-bt4
		h3 = h4
		line5(12) = 255
		btflag = 1
	    else if (bt4.eq.99999 .and. max(abs(bt1),abs(bt3),abs(bt3)).lt.9999) then
		bt4 = bt1+bt2-bt3
		h4 = h3
		line5(13) = 255
		btflag = 1
	    end if
c
		btcor = min(line5(10),line5(11),line5(12),line5(13))
c
c flags
	    btmax = max(abs(bt1),abs(bt2),abs(bt3),abs(bt4))
	    btrack = max(h1,h2,h3,h4)*0.01
c
c single ping bottom track data 
c
	if (btrack.lt.1100 .and.
     1	    btmax.lt.99999 .and. btcor.ge.220 .and. iq(i).eq.0) then
	  write(62,162)iens,bt1,bt2,bt3,bt4,gpx(i),gpy(i),head,pitch,roll
	end if
162	format(i6,4i8,2f20.10,3f10.3)
c
	    btx = vxscale * sig_up*(zsg(1)*bt1+zsg(2)*bt2)
	    bty = vxscale * sig_up*(zsg(3)*bt3+zsg(4)*bt4)
c
c The factor of two here accounts for larger factor in vzscale
c
	    btz1 = vzscale * sig_up*(bt1+bt2)*2.
	    btz2 = vzscale * sig_up*(bt3+bt4)*2.
c	    btz  = vzscale * sig_up*(bt1+bt2+bt3+bt4)
c	    bte = vzscale * (bt1+bt2-bt3-bt4)
c
c vertical-plane misalignment
c
		xr  =  btx*cos(arg1)+btz1*sin(arg1)
		xzr = -btx*sin(arg1)+btz1*cos(arg1)
		yr  =  bty*cos(arg2)+btz2*sin(arg2)
		yzr = -bty*sin(arg2)+btz2*cos(arg2)
		btx = xr
		bty = yr
		btz = (xzr+yzr)*0.5
		bte = (xzr-yzr)*0.5
c
		speed = sqrt(btx*btx+bty*bty)*0.01
		dir = atan2(bty,btx)*180./pi
c
c	write(37,238)iens,time,bt1,bt2,bt3,bt4,btx,bty,btz,bte
238	format(i5,f10.2,4i6,4f8.1)
c	write(25,105)iens,time,head,hsdev,pitch,psdev,roll,rsdev,temp,0.01*bte,lon
c     1				,gyro(i)
c
c h is the true vertical distance from the transducer heads to the bottom
c   (Chris Humprey, RDI, feb.1994)
c
	    hens = (h1+h2+h3+h4)*0.25*0.01
	    hprof = hens*100.*c30
	    if (btrack.gt.1300) then
		    hprof = line2(7)*ncell
	    end if
c
c ----->Averaging<------
c
	    if (time.lt.time_old+dt) then
		ntot = ntot+1
		ix = (lon-lon0)/dlon
		iy = (lat-lat0)/dlat
		if (ix.lt.100 .and. iy.lt.200 .and. ix.gt.0 .and. iy.gt.0) then
			bot = h_seab(ix,iy)
			h(ntot) = bot
			havg = havg+bot
		else
			bot = 9999
			h(ntot) = bot
			havg = havg+bot
		end if
c
c	Convert intensities to backscatter
c
		do 154 ib=1,ncell
		  range = float(ib*line2(7)+1000)*0.01
		  do 155 j=1,4
			line(ib,10+j) = db_scale*float(line(ib,10+j))+20*log(range)+2.*e_urick*range
155	  	  continue
154		continue
		goto 994
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
c
c Estimating bottom from echo-intensity maximum (NOT WORKING YET)
c
		xint(1) = -999.
		xint(2) = -999.
		xint(3) = -999.
		xint(4) = -999.
		xxint = -999.
		do 16 j=1,4
			ibot(j) = ncell
16		continue		
		do 14 ib=2,ncell
		do 15 j=1,4
c
c difference in backscatter between bin-ib and bin-1:
c
			df(j) = (line(ib,10+j) - line(1,10+j))
			if (df(j).gt.xint(j)) then
				xint(j) = df(j)
				ibot(j) = ib
			end if
15		continue
14		continue
		bot1 = (ibot(1)+ibot(2)+ibot(3)+ibot(4))*0.25
c
c sorting (Press et al., 1992, p.320)
c
c		do 17 j=2,4
c			xbot = ibot(j)
c			xxint = xint(j)
c			do 117 j2=j-1,1,-1
c				if (ibot(j2).le.xbot) goto 217
c				ibot(j2+1) = ibot(j2)
c				xint(j2+1) = xint(j2)
117			continue
c			j2 = 0
217			xint(j2+1) = xxint
c   			ibot(j2+1) = xbot
17		continue
c
c Nov. 21, 2003:
		nzmin2 = min(ibot(1),ibot(2),ibot(3),ibot(4))
 		nzmax2 = max(ibot(1),ibot(2),ibot(3),ibot(4))
		if (nzmax2-nzmin2.gt.2) then
			nzmax1 = 0
		else
			nzmax1 = nzmin2
		end if
		hadcp = nzmax1*line2(7)*0.01+blanking-line2(7)*0.005
		h(ntot) = bot
		havg = havg+bot
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
994	continue
		nzmax1 = int( 0.85*(bot-8.6)/(float(line2(7))*0.01) )-1
		if (btrack.lt.1100 .and. abs(bte).lt.bte_crit .and.
     1		    btmax.lt.99999 .and. btcor.ge.220 .and. iq(i).eq.0) then
		 	nbt = nbt+1
			btxbar = btxbar+btx
			btybar = btybar+bty
			bt1bar = bt1bar+bt1
			bt2bar = bt2bar+bt2
			bt3bar = bt3bar+bt3
			bt4bar = bt4bar+bt4
			vsx(nbt) = btx
			vsy(nbt) = bty
			vsz(nbt) = btz
		 	btzbar = btzbar+btz
			btebar = btebar+bte
c
c gpxbar and gpybar are used for BT calibrations only and thus are averaged only with BT data
c
			gpxbar = gpxbar+gpx(i)
			gpybar = gpybar+gpy(i)
			hbar = hbar+hens
		 	nzmax = int((hprof*0.85)/line2(7))-1
			iref = 1
		    else
			nzmax = nzmax1
			iref = 2
		    end if
cccc --->
	      	    nzmax = min(ncell,nzmax)
cccc --->
		    tbar = tbar+temp
		    cheadbar = cheadbar+cos(head*d2r)
		    sheadbar = sheadbar+sin(head*d2r)
		    cpbar = cpbar+cos(pitch*d2r)
		    spbar = spbar+sin(pitch*d2r)
		    crbar = crbar+cos(roll*d2r)
		    srbar = srbar+sin(roll*d2r)
		    hsdbar = hsdbar+hsdev
		    if (abs(lat).gt.1) then
			xbar = xbar+lat
			ybar = ybar+lon
			gpxbar1 = gpxbar1+gpx(i)
			gpybar1 = gpybar1+gpy(i)
			headbar = headbar + head
			nxy = nxy+1
			ugps(nxy) = gpx(i)
			vgps(nxy) = gpy(i)
			heading(nxy) = head
		    end if
c
	    do 4 ib=1,nzmax
		    j1 = ib*sc(1)+0.5
		    j2 = ib*sc(2)+0.5
		    j3 = ib*sc(3)+0.5
		    j4 = ib*sc(4)+0.5
		    jmax = max(j1,j2,j3,j4)
c
	    if (id(ib).lt.9999 .and. jmax.le.nzmax) then
c
c downward looking xducers
c
		    vxe = vxscale * sig_up*(zsg(1)*d1(j1)+zsg(2)*d2(j2))
		    vye = vxscale * sig_up*(zsg(3)*d3(j3)+zsg(4)*d4(j4))
		    vz1 = vzscale * sig_up*(d1(j1)+d2(j2))*2.
		    vz2 = vzscale * sig_up*(d3(j3)+d4(j4))*2.
c		    vze = vzscale * sig_up*(d1(j1)+d2(j2)+d3(j3)+d4(j4))
c		    ve  = vzscale *        (d1(j1)+d2(j2)-d3(j3)-d4(j4))
c
c vertical-plane misalignment
c
		xr  =  vxe*cos(arg1)+vz1*sin(arg1)
		xzr = -vxe*sin(arg1)+vz1*cos(arg1)
		yr  =  vye*cos(arg2)+vz2*sin(arg2)
		yzr = -vye*sin(arg2)+vz2*cos(arg2)
		vxe = xr
		vye = yr
		vze = (xzr+yzr)*0.5
		ve  = (xzr-yzr)*0.5
c
		    ivel = ivel+1
		    pg_min = float(min(line(ib,15),line(ib,16),line(ib,17),line(ib,18)))
		    cor = float(min(line(ib,7),line(ib,8),line(ib,9),line(ib,10)))
		    scatter = float(min(line(ib,11),line(ib,12),line(ib,13),line(ib,14)))
		    if (abs(vze).lt.ve_crit .and. (pg_min.gt.50. .or. ncol.eq.10)) then
			    ubar(ib) = ubar(ib)+vxe
			    vbar(ib) = vbar(ib)+vye
			    wbar(ib) = wbar(ib)+vze
			    ebar(ib) = ebar(ib)+ve
			    cbar(ib) = cbar(ib)+cor
			    sbar(ib) = sbar(ib)+scatter
			    nbar(ib) = nbar(ib)+1
		sxvar(ib,nbar(ib)) = vxe
		syvar(ib,nbar(ib)) = vye
		szvar(ib,nbar(ib)) = vze
		    end if
	    end if
4	continue
c
c ---------------> FORMING THE AVERAGE <--------------------------
c
	    else
		nprof = nprof+1
c
c assign a bottom depth hbar if % good bottom track exceeds ipg_bt_crit
c
		    ipg_bt = int(100*float(nbt)/float(ntot))
		    if (ipg_bt.gt.ipg_bt_crit) then
			hbar = hbar/float(nbt)
		    else
			hbar = havg/float(ntot)
			hdev = 0.
			do 52 itot=1,ntot
				dh = h(itot)-hbar
				hdev = hdev+dh*dh
52			continue
			hdev = sqrt(hdev/float(ntot-1))
c			hbar = -1
		    end if
c
		    if (nxy.gt.0) then
			xbar = xbar/float(nxy)
			ybar = ybar/float(nxy)
			gpxbar1 = gpxbar1/float(nxy)
			gpybar1 = gpybar1/float(nxy)
		median = nxy/2
		call hpsort(nxy,ugps)
		svar1 = gpxbar1-ugps(median)
		call hpsort(nxy,vgps)
		svar2 = gpybar1-vgps(median)
		gps_var = sqrt(svar1*svar1+svar2*svar2)*100.
		call hpsort(nxy,heading)
			headbar = headbar/float(nxy)
			hvar = headbar-heading(median)
			hsdev = sqrt(hvar*hvar)
		    end if
		    tbar = tbar/float(ntot)
		    cheadbar = cheadbar/float(ntot)
		    sheadbar = sheadbar/float(ntot)
		    cpbar = cpbar/float(ntot)
		    spbar = spbar/float(ntot)
		    crbar = crbar/float(ntot)
		    srbar = srbar/float(ntot)
		    headbar = atan2(sheadbar,cheadbar)/d2r
c
c compass calibration
c
	arg = headbar*d2r
	headbar = headbar-(alpha+s1head*sin(arg)+s2head*sin(2.*arg)+
     1                    	 c1head*cos(arg)+c2head*cos(2.*arg))
	arg = headbar*d2r
	pbar = atan2(spbar,cpbar)/d2r
	rbar = atan2(srbar,crbar)/d2r
c
c convert data to earth co-ordinates
c
	rr0 = rbar*d2r
	pp0 = pbar*d2r
c	ka0 = sqrt( 1.-sin(pp0)*sin(pp0)*sin(rr0)*sin(rr0) )
c	pp0 = asin( sin(pp0)*cos(rr0)/ka0 )
c
	cp0 = cos(pp0)
	sp0 = sin(pp0)
	cr0 = cos(rr0)
	sr0 = sin(rr0)
	ch0 = cos(arg)
	sh0 = sin(arg)
c
	hsdbar = hsdbar/float(ntot)
	if (nbt.gt.0) then
		xnbt = 1./float(nbt)
		btxbar = btxbar*xnbt
		btybar = btybar*xnbt
		btzbar = btzbar*xnbt
		btebar = btebar*xnbt
		median = nbt/2
		call hpsort(nbt,vsx)
		svar1 = btxbar-vsx(median)
		call hpsort(nbt,vsy)
		svar2 = btybar-vsy(median)
		bt_var = sqrt(svar1*svar1+svar2*svar2)
		ibt1 = int(bt1bar*xnbt)
		ibt2 = int(bt2bar*xnbt)
		ibt3 = int(bt3bar*xnbt)
		ibt4 = int(bt4bar*xnbt)
		gpxbar = gpxbar*xnbt
		gpybar = gpybar*xnbt
c light screen
c		if (bt_var.lt.40 .and. hsdev.lt.80.) then
	  write(72,162)iens,ibt1,ibt2,ibt3,ibt4,gpxbar,gpybar,headbar,pbar,rbar
c		end if
c
		btxx =  btxbar*(ch0*cr0+sh0*sr0*sp0)
     1		       +btybar*sh0*cp0
     1		       +btzbar*(ch0*sr0-sh0*cr0*sp0)
		btyy = -btxbar*(sh0*cr0-ch0*sr0*sp0)
     1		       +btybar*ch0*cp0
     1		       -btzbar*(sh0*sr0+ch0*sp0*cr0)
		btzz = -btxbar*sr0*cp0
     1                 +btybar*sp0
     1                 +btzbar*cp0*cr0
		btee = btebar
c
147		format(i6,3i5,2f9.3,2f6.0,i4,4f7.1,f7.1)
c
		speed = sqrt(btxx*btxx+btyy*btyy)*0.01
		dir = atan2(btyy,btxx)*180./pi
c
		aspeed = sqrt(gpxbar*gpxbar+gpybar*gpybar)
		adir = atan2(gpybar,gpxbar)*180./pi
	 	write(6,147)iens,nxy,ntot,nbt,xbar,ybar,hbar,hdev,ipg_bt,btxx,btyy,btzz,btee,headbar
 	 	write(46,147)iens,nxy,ntot,nbt,xbar,ybar,hbar,hdev,ipg_bt,btxx,btyy,btzz,btee,headbar
c
c calibration
c
c		ical = ical+1
		speed = sqrt(btxx*btxx+btyy*btyy)*beta
c		write(36,105)iens,time,headbar,pbar,rbar,temp,btxx*beta,btyy*beta,temp,xbar,ybar,speed
	else
		vsvar = gps_var
	 	write(6,147)iens,nxy,ntot,nbt,xbar,ybar,hbar,hdev,ipg_bt,gpxbar1,gpybar1,gps_var,btee,headbar
 	 	write(46,147)iens,nxy,ntot,nbt,xbar,ybar,hbar,hdev,ipg_bt,btxx,btyy,btzz,btee,headbar
		btxx = 0.
		btyy = 0.
		btzz = 0.
	end if
c
c estimate velocity of ship:
c
		iref = 0
		shipbar = sqrt(btxx*btxx+btyy*btyy)
		if (ipg_bt.gt.ipg_bt_crit .and. shipbar.lt.1000) then
			xship = btxx*beta
			yship = btyy*beta
			t = time-dt*0.5
			ax = 1.
			bx = 0.
			ay = 1.
			by = 0.
			vsvar = bt_var
			iref = 1
		else
c
c turn GPS velocity from m/s into cm/s
c
			xship = -gpxbar1*100.
			yship = -gpybar1*100.
			vsvar = gps_var
			iref = 2
		end if
c
c Use only GPS for ship speed determination
c
c		xship = -gpxbar1*100.
c		yship = -gpybar1*100.
c		vsvar = gps_var
		if (vsvar.gt.999) vsvar = 999
	vs_bt  = sqrt(xship*xship+yship*yship)
	vh_bt =  atan2(yship,xship)/d2r
	vs_bt = vs_bt*0.01
		shipbar = sqrt(xship*xship+yship*yship)*0.01
		uzbar = 0.
		vzbar = 0.
		wzbar = 0.
		ezbar = 0.
		uzbar2 = 0.
		vzbar2 = 0.
		nzbar = 0
		flag = 0
		blank = float(line2(9))*0.01
		dz = float(line2(7))*0.01
		offset = blank+float(zoff)+dz/2.
		do 7 ib=ncell,1,-1
c		do 7 ib=nzmax,1,-1
			ibin_z   = int(offset+ib*dz)
			if (nbar(ib).ge.2) then
				ubar(ib) = ubar(ib)/float(nbar(ib))
				vbar(ib) = vbar(ib)/float(nbar(ib))
				wbar(ib) = wbar(ib)/float(nbar(ib))
				ebar(ib) = ebar(ib)/float(nbar(ib))
				cbar(ib) = cbar(ib)/float(nbar(ib))
				sbar(ib) = sbar(ib)/float(nbar(ib))
c
c tilt sensors fixed to adcp
c
				vxe =  ubar(ib)*(ch0*cr0+sh0*sr0*sp0)
     1				      +vbar(ib)*         sh0*cp0
     1				      +wbar(ib)*(ch0*sr0-sh0*cr0*sp0)
				vye = -ubar(ib)*(sh0*cr0-ch0*sr0*sp0)
     1            		      +vbar(ib)*       ch0*cp0
     1           		      -wbar(ib)*(sh0*sr0+ch0*sp0*cr0)
				vze = -ubar(ib)*sr0*cp0
     1             		      +vbar(ib)*sp0
     1              		      +wbar(ib)*cp0*cr0
				t = time-dt*0.5
				u = (vxe*beta-xship)
				v = (vye*beta-yship)
				w = vze
				e = ebar(ib)
				ww(ib) = w
				ipg = int(nbar(ib)/float(ntot)*100.)
c
c 				compute depth averages
c
				if (ipg.gt.ipg_crit .and. ib.lt.nzmax) then
					uzbar = uzbar+u
					vzbar = vzbar+v
					wzbar = wzbar+w
					ezbar = ezbar+e
					uzbar2 = uzbar2+u*u
					vzbar2 = vzbar2+v*v
					nzbar = nzbar+1
				end if
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
cccc
cccc 			WRITE AWAY DATA
cccc
		uspd = sqrt(u*u+v*v)
		if (shipbar.lt.8 .and. abs(xbar).gt.0 .and. uspd.lt.200) then
c
c
c ubar(ib), sxvar(ib,ibar), etc. are all in ADCP co-ordinates (cm/s)
c
		do 12 ibar=1,nbar(ib)
			sxx(ibar) = sxvar(ib,ibar)
12		continue
		call hpsort(nbar(ib),sxx)
		median = nbar(ib)/2
		svar1 = ubar(ib)-sxx(median)
		do 3 ibar=1,nbar(ib)
			sxx(ibar) = syvar(ib,ibar)
3		continue
		call hpsort(nbar(ib),sxx)
		svar2 = vbar(ib)-sxx(median)
		svar = sqrt(svar1*svar1+svar2*svar2)
			ndata = ndata+1
			write(26,114)t,ybar,xbar,ib,w,e
     1			,hbar,shipbar,headbar,ipg,ipg_bt,tbar,sbar(ib),svar,vsvar,cbar(ib),iref,ibin_z
			write(24,114)t,ybar,xbar,ib,u,v
     1			,hbar,shipbar,headbar,ipg,ipg_bt,tbar,sbar(ib),hsdev,vsvar,cbar(ib),iref,ibin_z
		else
			write(44,114)t,ybar,xbar,ib,u,v
     1			,hbar,shipbar,headbar,ipg,ipg_bt,tbar,sbar(ib),hsdev,vsvar,cbar(ib),iref,ibin_z
		end if
114		format(f8.1,f10.4,f10.4,i6,2f7.1
     1			,f6.0,f4.1,f6.0,2i4,3f7.1,2f6.1,i2,i4)
115		format(f8.1,f10.4,f10.4,i6,2f7.1
     1			,f6.0,f4.1,f6.0,2i4,3f7.1,2f6.1,i2,f8.1,i4,i5)
cccc
cccc
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
			end if
7		continue
c
		if (nzbar.gt.1 .and. shipbar.lt.10 .and. abs(xbar).gt.0) then
			uzbar = uzbar/float(nzbar)
			vzbar = vzbar/float(nzbar)
			uzsdev = sqrt((uzbar2-uzbar*uzbar/nzbar)/float(nzbar-1))
			vzsdev = sqrt((vzbar2-vzbar*vzbar/nzbar)/float(nzbar-1))
			wzbar = wzbar/float(nzbar)
			ezbar = ezbar/float(nzbar)
			write(76,178)t,nzbar,wzbar,btxx*beta,btyy*beta,btzz,btee,ipg_bt,gpxbar1*100,gpybar1*100,var*100
178			format(f8.1,i4,5f7.1,i4,3f7.1)
			n76 = n76+1
			z = 999.
			izbar = izbar+1

			write(82,118)t,ybar,xbar,nzbar,uzbar,vzbar
     1			,hbar,shipbar,headbar,ipg,ipg_bt,tbar,uzsdev,hsdev,vsvar,vzsdev,iref,ibin_z
     1			,(line1(k),k=2,5)
118			format(f8.1,f10.4,f10.4,i6,2f7.1
     1			,f6.0,f4.1,f6.0,2i4,3f7.1,2f6.1,i2,4i2,i4)
		end if
c		write(79,179)t,u,v,w,e,wzbar,ezbar,btzz*beta,iens
179		format(f8.0,7f7.1,i6)
c
c		write(77,114)t,ybar,xbar,nzbar,uzbar,vzbar
c     1			,hbar,shipbar,headbar,ipg,ipg_bt,tbar,vh_bt,hsdbar,vsvar,cbar(1),iref
c
		t = time+dt*0.5
c
c ---------->NEW AVERAGE (INITIALIZING)<-------------------
c
		iflag2 = 0
		do 5 ib=1,nzmax
			j1 = ib*sc(1)+0.5
			j2 = ib*sc(2)+0.5
			j3 = ib*sc(3)+0.5
			j4 = ib*sc(4)+0.5
			jmax = max(j1,j2,j3,j4)
		if (id(ib).lt.9999 .and. jmax.le.nzmax) then
c
c downward looking xducers
c
			vxe = vxscale * sig_up*(zsg(1)*d1(j1)+zsg(2)*d2(j2))
			vye = vxscale * sig_up*(zsg(3)*d3(j3)+zsg(4)*d4(j4))
c			vze = vzscale * sig_up*(d1(j1)+d2(j2)+d3(j3)+d4(j4))
		        vz1 = vzscale * sig_up*(d1(j1)+d2(j2))*2.
		        vz2 = vzscale * sig_up*(d3(j3)+d4(j4))*2.
c			ve  = vzscale *        (d1(j1)+d2(j2)-d3(j3)-d4(j4))
c
c vertical-plane misalignment
c
			xr  =  vxe*cos(arg1)+vz1*sin(arg1)
			xzr = -vxe*sin(arg1)+vz1*cos(arg1)
			yr  =  vye*cos(arg2)+vz2*sin(arg2)
			yzr = -vye*sin(arg2)+vz2*cos(arg2)
			vxe = xr
			vye = yr
			vze = (xzr+yzr)*0.5
			ve  = (xzr-yzr)*0.5
c
			if (abs(vze).lt.ve_crit) then
				ubar(ib) = vxe
				vbar(ib) = vye
				wbar(ib) = vze
				ebar(ib) = ve
				cbar(ib) = cor
				sbar(ib) = scatter
c				cor_min(ib) = 999.
				nbar(ib) = 1
		sxvar(ib,nbar(ib)) = vxe
		syvar(ib,nbar(ib)) = vye
		szvar(ib,nbar(ib)) = vze
			else
				ubar(ib) = 0.
				vbar(ib) = 0.
				wbar(ib) = 0.
				ebar(ib) = 0.
				cbar(ib) = 0.
				sbar(ib) = 0.
c				cor_min(ib) = 999.
				nbar(ib) = 0
			end if
		else
			ubar(ib) = 0.
			vbar(ib) = 0.
			wbar(ib) = 0.
			ebar(ib) = 0.
			cbar(ib) = 0.
			sbar(ib) = 0.
c			cor_min(ib) = 999.
			nbar(ib) = 0
		end if
5		continue
		do 6 ib=nzmax+1,nz
			ubar(ib) = 0.
			vbar(ib) = 0.
			wbar(ib) = 0.
			ebar(ib) = 0.
			cbar(ib) = 0.
			sbar(ib) = 0.
c			cor_min(ib) = 999.
			nbar(ib) = 0
6		continue
c
		if (abs(lat).gt.1.) then
			xbar = lat
			ybar = lon
			gpxbar1 = gpx(i)
			gpybar1 = gpy(i)
			headbar = head
			nxy = 1
			ugps(nxy) = gpx(i)
			vgps(nxy) = gpy(i)
		else
			xbar = 0.
			ybar = 0.
			gpxbar1 = 0.
			gpybar1 = 0.
			headbar = 0.
			nxy = 0
		end if
c
		ntot = 1
		cheadbar = cos(head*d2r)
		sheadbar = sin(head*d2r)
		cpbar = cos(pitch*d2r)
		spbar = sin(pitch*d2r)
		crbar = cos(roll*d2r)
		srbar = sin(roll*d2r)
		hsdbar = hsdev
		tbar = temp
c
		if (btrack.lt.1100 .and. abs(bte).lt.bte_crit .and.
     1		    btmax.lt.99999 .and. btcor.ge.220 .and. iq(i).eq.0) then
			hbar = hens
			havg = 0.
			hdev = 0.
			btxbar = btx
			btybar = bty
			btzbar = btz
			btebar = bte
			bt1bar = bt1
			bt2bar = bt2
			bt3bar = bt3
			bt4bar = bt4
			gpxbar = gpx(i)
			gpybar = gpy(i)
			nbt = 1
		else
			hbar = 0.
			havg = bot
			btxbar = 0.
			btybar = 0.
			btzbar = 0.
			btebar = 0.
			bt1bar = 0.
			bt2bar = 0.
			bt3bar = 0.
			bt4bar = 0.
			gpxbar = 0.
			gpybar = 0.
			xship = 0.
			yship = 0.
			nbt = 0
		end if
		time_old = time
		i_old = i
	  end if
2	continue
	close (8)
	close (9)
	close (26)
	close (25)
	close (28)
	close (24)
c	close (36)
	close (76)
c	close (82)
	nens = nens-1
c
	write(6,*)'writing to scratch',label
	open (unit=10,file='scratch')
	write(10,121)path,label,label,label,label,label,label
121	format(a23,'/traos.',a3,7x,'calibration.',a3,'trbad.',a3,'depth.',a3,
     1		'botrk.',a3,'bin02.',a3)
	close(10)
	write(6,*)'reading from scratch',label
	open (unit=10,file='scratch')
	read(10,122)name,cal_name,name1,zname,btname,biname
122	format(a36,4x,a15,4a9)
	close(10)
c calibration 
c	open (unit=36,file='calibration.dat')
c	open (unit=37,file=cal_name)
146	format(a91)
c	do 36 i=1,ical
c		read(36,146)all
c		write(37,146)all
36	continue
c	close (36)
c	close (37)
c bottom track
c		open(unit=76,file='fort.76')
c		open(unit=66,file=btname)
c		do 38 i=1,n76
c			read(76,138)all1
c			write(66,139)all1
138			format(a51)
38		continue
c		close (76)
c		close (66)
c depth avg.
c		open(unit=82,file='fort.82')
c		open(unit=85,file=zname)
c		do 37 i=1,izbar
c			read(82,139)all1
c			write(85,139)all1
139			format(a93)
37		continue
c		close (82)
c
c =======>SCREENING VELOCITY DATA<================
c
	open (unit=13,file=name)
c	open (unit=38,file=name1)
	open (unit=38,file='trans.screen')
	open (unit=10,file='transect.dat')
	open (unit=26,file='vship.dat')
	nprofile = 0
	npass = 0
c
	if (ndata.gt.0) then
		read(26,114)t,ybar,xbar,ib,w,e,
     1			hbar,shipbar,headbar,ipg,ipg_bt,tbar,scatter,hdev,vsvar,cor,iref,ibin_z
		read(10,114)t,ybar,xbar,ib,u,v,
     1			hbar,vs_bt,headbar,ipg,ipg_bt,tbar,scatter,hsdev,vsvar,cor,iref,ibin_z
		ib_old = ib
		ib_old1 = ib
		ndata0 = 1
		do 8 i=2,ndata
			read(26,114)t,ybar,xbar,ib,w,e,
     1				hbar,shipbar,headbar,ipg,ipg_bt,tbar,scatter,hdev,vsvar,cor,iref,ibin_z
			read(10,114)t,ybar,xbar,ib1,u,v
     1				,hbar,vs_bt1,headbar1,ipg,ipg_bt,tbar,dummy,hsdbar,vsvar,cor,iref,ibin_z
			if (ib.gt.ib_old) then
				nprofile = nprofile+1	
			end if
			ib_old = ib
c very light sreen
		if (ibin_z.lt.hbar*0.85 .and. ipg.gt.ipg_crit) then
c light screen
c		if (ibin_z.lt.hbar*0.85 .and. ipg.gt.ipg_crit .and. abs(e).lt.15 .and. vsvar.lt.40 .and. hsdbar.lt.80.) then
c heavy screen
c		if (ipg.gt.ipg_crit .and. abs(e).lt.20 .and. vsvar.lt.16 .and. hsdbar.lt.16.) then
			if (ib.gt.ib_old1) then
				npass = npass+1	
			end if
			ib_old1 = ib
			write(13,115)t,ybar,xbar,ib1,u,v
     1			,hbar,vs_bt1,headbar1,ipg,ipg_bt,hsdbar,w,e
     1			,vsvar,cor,iref,hdev,ibin_z,npass+1
			ndata0 = ndata0+1
		else
148			format(f8.1,f10.4,f9.4,i6,2f7.1
     1				,f6.1,f8.1,i4,2f7.1)
			write(38,148)t,ybar,xbar,ib1,u,v,dvs,dvh,ipg,hbar,hsdev1
		end if
8	   continue
	end if
	if (nens.gt.0) then
		write(6,142)float(npass)/float(nprofile)*100.,npass,nprofile,ndata0
		write(4,142)float(npass)/float(nprofile)*100.,npass,nprofile,ndata0
	end if
142	format(/'percent of ensemble passing screening: ',f7.1/
     1		'number of ensemble passing data      : ',i6/
     1		'number of ensemble                   : ',i6/
     1		'number of data passing               : ',i6/)
	  close (13)
	  close (26)
	  close (38)
	  close(10)
1	continue 
c
101	format(2f8.1,i3,f10.6,5f7.1,f5.1,f5.0)
102	format(2i2,i3,1x,a23,2x,i2,3x,i1,3x,i1,5f6.1,f8.5)
103	format(a9,1x,i6,1x,a3,i3,2i2,5f6.1,f7.4)
104	format(f10.6,5i6)
105	format(1x,i5,f15.2,2f7.2,4f7.2,f6.2,f11.5,2f11.5)
106	format(i5,f10.4,4i5,2f10.2,2f7.0)
108	format(a45)
109	format(f7.3,6i5)
110	format(2i5,f15.5,4f10.1)
944 	format(a98)
c
	stop
	end
	subroutine attention
	write(6,*)'#################################################'
	write(6,*)
	write(6,*)'SCREENING WORKS ONLY FOR BINS ORDER SEQUENTIALLY'
	write(6,*)'	IN REVERSE ORDER'
	write(6,*)
	write(6,*)'#################################################'
	return
	end
	function select(k,n,arr)
	integer k,n
	real select,arr(n)
	integer i,ir,j,l,mid
	real a,temp
	l = 1
	ir = n
1	if (ir-l.eq.1) then
	   if (ir-l.eq.1) then
		if (arr(ir).lt.arr(l)) then
		   temp = arr(l)
		   arr(l) = arr(ir)
		   arr(ir) = temp
		end if
	   end if
	   select = arr(k)
	   return
	else
	   mid = (l+ir)/2
	   temp = arr(mid)
	   arr(mid) = arr(l+1)
	   arr(l+1) = temp
	   if (arr(l).gt.arr(ir)) then
		temp = arr(l)
		arr(l) = arr(ir)
		arr(ir) = temp
	   end if
	   if (arr(l+1).gt.arr(ir)) then
		temp = arr(l+1)
		arr(l+1) = arr(ir)
		arr(ir) = temp
	   end if
	   if (arr(l).gt.arr(l+1)) then
		temp = arr(l)
		arr(l) = arr(l+1)
		arr(l+1) = temp
	   end if
	   i = l+1
	   j = ir
	   a = arr(l+1)
3	   continue
		i = i+1
	   if (arr(i).lt.a) goto 3
4	   continue
		j = j-1
	   if (arr(j).gt.a) goto 4
	   if (j.lt.i) goto 5
	   temp = arr(i)
	   arr(i) = arr(j)
	   arr(j) = temp
	   goto 3
5	   arr(l+1) = arr(j)
	   arr(j) = a
	   if (j.ge.k) ir=j-1
	   if (j.le.k) l=i
	end if
	goto 1
	end
ccccccccccccccccccccccccccc
	subroutine hpsort(n,ra)
	integer n
	real ra(n)
	integer i,ir,j,l
	real rra
	if (n.lt.2) return
	l = n/2+1
	ir = n
10	continue
		if (l.gt.1) then
		  l = l-1
		  rra = ra(l)
		else
		  rra = ra(ir)
		  ra(ir) = ra(1)
		  ir = ir-1
		  if (ir.eq.1) then
			ra(1) = rra
			return
		  end if
		end if
		i = l
		j = l+1
20		if (j.le.ir) then
		  if (j.lt.ir) then
			if (ra(j).lt.ra(j+1)) j = j+1
		  end if
		  if (rra.lt.ra(j)) then
			ra(i) = ra(j)
			i = j
			j = j+j
		  else
			j = ir+1
		  end if
	goto 20
	end if
	ra(i) = rra
	goto 10
	end