program remap_subdata_extrapolate_set c modified from remap_subdata_extrapolate_1more to allow input of c number of bins to extrapolate c modified from remap_subdata_extrapolate to extrapolate 1 more bin shallower c made remap_subdata_extrapolate_1more c remap data from upward looking ADCPs c uses historical mean sound velocity profile to compute proper bin widths c and range error c reads pressure file to compute transducer depth c editted 10sep02 to extrapolate two shallowest values to one more shallow c value for interpolating to allow shallower data values c modified from remap_subdata to allow input of depth error evident from c comparison of target strength trans depth and seacat/mtpr c modified again 5 mar 02 to allow input to increase cutoff depth by c input amount c will take difference (example: 5m) and incrementally change the actual depth of the c data over the data profile using transducer depth as total depth c (before 5aug03 used average depth as total depth) c itype as of 8/26/96: c sets data bad due to sidelobe interference near surface if shallower than c itype=1: no cutoff c itype=2: trans_depth*(1-cos(theta))+bin_width +sidelobe adj c where theta is the transducer angle c reads offset between transducer head and pressure c sensor, mean density c write output file with time,transdepth,computed good data cutoff depth, c min_remapped_gooddata_depth, minimum_original_gooddata_depth, c 2nd_original_gooddata_depth, shallower_original_depth c changed 8/8/96 to write actual_depth,u,v,err rather than actual_depth,u,v,w c or for amp actual_depth,beam1,beam2,beam4 c writes output file in rd file format with 4 variables: actual_depth,u,v,err c changed 7/99 to allow input of transducer angle and nominal bw, pulse & blank c changed 11/99 to allow an amp code to add 4.5m to range compared with vel c range c c changed 2000 to input lat and trans depth to compute gravity, which should be c near 9.7806 for equatorial, it is printed in log file c $! c $! nArrays new_ymin new_ydel new_numbins outfile theta nom_bw nom_pulse nom_blnk c vel_or_amp tlat mean_tdepth(z) depth_adj sidelobe_deepening bin_ext c $! depthoffset,density,itype,transdepth_outfilename,actual_depth_outfilename c $! nftype start end ny nv delta_time_in_days !ny,nv=0 here, comp in read_array c $! input_file T T T T F F 1 37 c $! nftype start end ny nv delta_time_in_days !ny,nv=1 for pressure file c $! input_file format depth xbad c $! num_ctd_depths num_ctd_var c $! ctd_file_name nskip ctd_ydel c 1 10. 5. 60 'datka7p/ka07_10_5_adjdepth_ext.vel' 20. 8. 8. 4. 'vel' 0.0 258.6 6.48 2. 3 c 2.1 1.02505 2 'datka7p/ka07_1stdepth_adjdepth_ext.out' 'junk.out' c 2 10 00 20 6 02 23 00 25 12 02 0 0 .04166667 c 'datka7/ka07_tc.vel' T T T T F F 1 44 c 0 10 00 20 6 02 23 00 25 12 02 1 1 .04166667 c 'datka7/ka7hlf.300sc' '(1x,2i2,1x,3i2,16x,f8.3/)' 250 -99.999 c 500 2 c 'pprog/input/mean_snd.170w' 1 1. c $ exit c c c (time-x, depth-y, component-v) 3-D array name is "uv". c array remap is uv after sndvel depth correction and depth remapping. c array iuv is remap converted to integer for writing output file. c c c Dimensions of uv array - nx,ny,nv - are determined as follows: c c for RD arrays: c c nx determined from start/end time and increment c ny determined by bin range input c nv determined by input logical array "lu" c lu is hard coded in this program rather than in command file c c Input parameters: c c RECORD 1 c this changed c read(5,*,end=900) narrays, new_ymin, new_ydel, new_numbins, c outname theta nom_bw nom_pulse nom_blnk 'data_type' c tlat mean_tdepth depth_adj sidelobe_deepening c bin_ext c read(5,*,end=900) del_sensor, density, itype, outname2 outname3 c where c Narrays = number of input arrays (should be 1) c new_ymin = 1st depth of new remapped array c new_ydel = depth increment of new remapped array c new_numbins = numbins of new remapped array c del_sensor = offset depth to be subtracted from pressure sensor data to c get depth of transducer c density = mean density calculated from historical ctd's c itype = code to determine cutoff depth equation near surface from c sidelobe interference c 1 = no cutoff depth c 2 = cutoff depth eq Tdepth*{1-cos theta}+binwidth + sidelobe adj c outname2 = file name of output transdepth and min_good_data_depth c outname3 = file name of rd file with 4 variables: actual_depth,u,v,err c theta = transducer angle (20 or 30) c nom_bw = nominal bin width = 8 or 4 c nom_pulse = nomimal pulse length = 8 or 4 c nom_blnk = nomimal blanking = 4 or 2 c data_type = 'vel' or 'amp' c tlat = latidude for computing gravity c z = average depth of deployment c adjdepth = depth error evident from comparison c of target strength trans depth and seacat/mtpr c sladj is sidelobe adjustment, input to deeper cutoff depth for sidelobe error c bin_ext is number of bins shallower to linearly interpolate c c read rdi data c read(5,*) nftype,ImatDAT,LmatDAT,ny,nv,XDEL !ny and nv are zero, determined c ! from max/min bin# and read_array c where c nftype = input array type c 0: conventional time series (1 file per depth) c 1: RD data set from Proteus c 2: RD data set from memory c 3: any ASCII file with tabulated uv data c 4: atlas climatology file c 5: gpcp satellite rain data c 6: new format composite time series c -1: Proteus type format, but can adjust depths c -2: RD memory type format, but can adjust depths c ImatDAT(1-5) = for nftype= 0-2; start time/date for array. c LmatDAT(1-5) = for nftype= 0-2; end time/date for array c ny = array depth dimension (not used for nftype=1,2) c nv = array component dimension: not for nftype=1,2 c nv=iao = array operator : for nftype=1,2 only c 0, do nothing c 1, add to and store in previous array c 2, subtract from and store in previous array c !!! make sure previous array has same dimensions!!! c c XDEL = array time increment in days c c FOR nftype = 1 or 2 , NEXT RECORD c read(5,*) fnam,lu,nymin,nymax c if nftype = -1 or -2 also read ymin,ydel c where c FNAM = input file name c lu(1-4) = logical array c lu(1)= T; input zonal velocity c lu(2)= T; input meridional velocity c lu(3)= T; input number of hourly ensembles (nftype = 1) c or input vertical velocity (nftype = 2) c lu(4)= T; input error velocity (nftype = 2) c lu(5)= T; calculate speed from u,v ! must be false c lu(6)= T; calculate direction from u,v ! must be false c nymin = minimum bin number to use c nymax = maximum bin number to use c c read pressure data c read(5,*) nftype,ImatDAT,LmatDAT,ny,nv,XDEL !ny and nv = 1, nftype=0 c FOR nftype = 0, NEXT ny RECORDS c read (5,*) FNAM,FORM,YDEP,(xflg(iv),iv=1,nv) c FNAM = input file name c FORM = input file format c YDEP = depth value c xflg = bad data flag, one for each variable; if nftype = 6 and c the input file has a FLAGS header record, then the flags c read from the header will be used as the bad value flags c and xflg will be read but not used. c read snd velocity data: c read(5,*) nyerr,nverr c where nyerr = number of ctd depths == le 500 c nverr = number of ctd variables == le 4 c read(5,*) ctdfnam, nskip, ydel c where ctdfnam = ctd file name c nskip = number of lines to skip in ctd file c ydel = delta depth of ctd data, should be 1m c c parameter (max_matrix=2500000,max_arrays=10, c 1 max_work=60000,max_input_numbins=44) !! 40 bins in original rdi file parameter (max_matrix=5000000,max_matrix2=5000000,max_arrays=2, 1 max_work=60000,max_input_numbins=44) !! 40 bins in original rdi file parameter (max_matrix3=10000) dimension uv(max_matrix),y(200,max_arrays), 1 iuv(max_matrix2) !! iuv is interger output from remap dimension remap(max_matrix2),tr(max_matrix3) real actual_depth(max_input_numbins),transdepth real time_int real new_ymin,new_ydel,range_err(500,5) character*40 form,outname,outname2,outname3 COMMON/NAME/PROJ(2),MOOR,VNAM,ID,DELTin,DEPTH,REST COMMON/INFORM/ FORM,IDATE(5),LDATE(5) common/check/nchk,xbad COMMON/matrix/ ImatDAT(5),LmatDAT(5),XDEL common/peprdi/nrec,NBins,NBout,NBstart,NBend integer nx,ny,nv,narrays,out,ix,iy,iv integer ry,new_numbins,out2,itype,out3 c common/pepfile/out data xbad /-999.99/ integer*4 calc_size,iactual_depth(max_input_numbins) real theta,tlat,sladj real ctd(500,4),get_trans_depth,g,gp,z,density,del_sensor real nom_bw, nom_pulse, nom_blnk, bw,p,blnk integer nyerr,nverr,n,nskip,lx,ly,lv,ibin_ext character*50 ctdfnam character*80 line character*3 data_type out=4 ! output file number out2 = 3 out3 = 2 100 read(5,*,end=900,err=900) narrays,new_ymin,new_ydel,new_numbins, 1 outname,theta,nom_bw,nom_pulse,nom_blnk,data_type, 1 tlat,z,adjdepth,sladj,ibin_ext read(5,*,end=900,err=900) 1 del_sensor,density,itype,outname2,outname3 write(6,107) 'Output File Name: ',outname, 1 'Depth Output File Name: ',outname2, 1 ' new_ymin:',new_ymin,' new_ydel:',new_ydel, 1 ' depth subtracted from pressure data to get trans depth:', 1 del_sensor, 1 ' density:', density, ' #bins extrap:',ibin_ext 107 format(/2x,a,a/2x,a,a/a,f8.2,a,f6.2/a,f6.2/a,f9.5/a,i4) write(6,'(a,a)') ' Actual depth output file:',outname3 c if (theta.eq.20.) then bw = nom_bw * 1.085 p = nom_pulse * 1.085 blnk = nom_blnk * 1.085 elseif (theta.eq.30.) then bw = nom_bw p = nom_pulse blnk = nom_blnk else write(6,'(/a)') ' Error: transducer angle not 20 or 30 degrees' go to 900 endif write(6,'(4(a,f5.2)/)') 1 ' bw:',bw,' pulse:',p,' blnk:',blnk write(6,'(a,f6.2,a,f7.2,a,f6.2,a,f5.2)') 1 ' lat:',tlat,' avg_depth:',z,' depth_adj:',adjdepth, 1 ' sidelobe_adj:',sladj call fopen(out,outname,'new') call fopen(out2,outname2,'new') call fopen(out3,outname3,'new') if (itype.eq.1) then write(6,'(a)') ' no cutoff depth' write(out2,'(a,a)') 1 ' no_cutoff_depth transd calc_gd_d mapp_gd_d', 1 ' act_dep1 act_dep2 nxt_sha_act_d' elseif (itype.eq.2) then if (data_type.eq.'vel') then write(6,'(a,f3.0,a,f3.1)') 1 ' cutoff depth eq Tdepth*{1-cos',theta,'}+ bw + ',sladj write(out2,'(a,f3.0,a,f3.1,a,a)') 1 ' D*{1-cos',theta,'} + bw + ',sladj, 1 '} transd calc_gd_d mapp_gd_d', 1 ' act_dep1 act_dep2 nxt_sha_act_d' elseif (data_type.eq.'amp') then write(6,'(a,f3.0,a,f3.1)') 1 ' cutoff depth eq Tdepth*{1-cos',theta,'}+ .75bw + ',sladj write(out2,'(a,f3.0,a,f3.1,a,a)') 1 ' D*{1-cos',theta,'} + .75bw + ',sladj, 1 '} transd calc_gd_d mapp_gd_d', 1 ' act_dep1 act_dep2 nxt_sha_act_d' else write(6,'(a)') ' Error: must be amp or vel' endif else write(6,'(a)') ' Error* itype must be between 1 or 2' go to 900 endif ry = new_numbins c Now read RDI data (subroutine read_array from .lnk file) read(5,*) nftype,ImatDAT,LmatDAT,ny,nv,XDEL nf = 1 ! set array number = 1 10 call read_array(uv,y(1,nf),nx,ny,nv,nftype) write(6,'(a,i4,a,i4)') ' array #:',nf,' Array type:',nftype write(6,102) ' start: ',ImatDAT,' end: ',LmatDAT,' delta: ',XDEL 102 format(a,5i3,a,5i3,a,f10.5) c c c *** Check that matrix not to big *** c calc_size = nx*ny*nv if(calc_size.gt.max_matrix) then write(6,1011) 1011 format(' Increase uv matrix dimensions - max_matrix') go to 900 endif c c *** Check that matrix not to big *** c calc_size = nx*ry*nv if(calc_size.gt.max_matrix2) then write(6,'(a)') ' Increase remap matrix dimensions - max_matrix' go to 900 endif c Now read pressure data (subroutine read_array from .lnk file) read(5,*) nftype,ImatDAT,LmatDAT,ly,lv,XDEL nf = 2 ! now do array # 2 transducer depth call read_array(tr,y(1,nf),lx,ly,lv,nftype) write(6,'(a,i4,a,i4)') ' array #:',nf,' Array type:',nftype write(6,105) ' start: ',ImatDAT,' end: ',LmatDAT,' delta: ',XDEL 105 format(a,5i3,a,5i3,a,f10.5) c c *** Check that matrix not to big *** c calc_size = lx*ly*lv if(calc_size.gt.max_matrix3) then write(6,'(a)') ' Increase pres matrix dimensions - max_matrix' go to 900 endif c call dmp_array(uv,nx,ny,nv) c call dmp_array(tr,lx,ly,lv) c c read ctd file info and compute range_error and corrected range as a function c of 1m uncorrected vertical range values (from transducer) read(5,*) nyerr,nverr write(6,'(1x,2i4)') nyerr,nverr if (nyerr.gt.500.or.nverr.gt.4) then write(6,'(a)') ' Error, increase matrix dimensions for ctd data' go to 900 endif c ydel should be 1m read(5,*) ctdfnam,nskip,ydel WRITE (6,99) CTDFNAM 99 Format(//10x,'Input Filename: ',a) call fopen(1,CTDFNAM,'old') do n=1,nskip read(1,'(a80)',end=50) line end do c expect ctd file to have scan, pressure, soundvel, scan is not considered a c variable, ie ctd(iy,1) are pressures and ctd(iy,2) are snd vels c this program assumes sound vel is ctd(iy,2) do 33 iy=1,nyerr read (1,'(8x,f9.2,1x,f8.1)') (ctd(iy,iv),iv=1,nverr) 33 end do !iy 50 close(1) c c For date of test days written to log file only c Get time increment (XDEL in days) and convert to time_int (in minutes) c Then recreate date from getdt2 subroutine by incrementing time (ix) c if (XDEL.eq.1) then time_int = 1440. elseif (XDEL.ge..04.and.XDEL.lt..05) then time_int = 60. elseif (XDEL.ge..02.and.XDEL.lt..03) then time_int = 30. elseif (XDEL.ge..01.and.XDEL.lt..02) then time_int = 15. else write(6,'(a)') ' Error: time increment not a day, an hour', 1 ' 30 or 15 minutes' go to 900 endif c Convert pressure sensor to transducer depth: use pres=density*g*depth c 2 c g' = 9.80616(1-.002644 cos 2@ + .000007 cos 2@) where @ is latitude c -6 c g = g' + 2.202 * 10 z c tlat = 0. c z = 140. gp = 9.80616*(1-.002644*(cosd(2.*tlat)) 1 + .000007*(cosd(2.*tlat))**2) c for lat=0 gp=9.7803012 c for z=140 addition to gp below is .00030828 g = gp + 2.202 * z * .001 * .001 write(6,'(a,f10.6)') ' g:',g c comment out following c g = 9.780609 ! gravity at 0 lat, average depth 140m c density = 1.02446 ! mean 0 to 275m at 156E from 9 historical ctds c del_sensor = dist m pres sensor below transducer. imput do ix = 1,nx transdepth=get_trans_depth(ix,tr,lx,ly,lv,density,g,del_sensor) c compute corrected vertical range array using correct sound velocity c assumes ctd file has sound velocities from 1 to 500m call vert_range_error(range_err,ctd,nyerr,nverr,transdepth, 1 data_type) c get actual depths of bins using corrected range array and transducer depth c bin 1 is still the deepest bin call get_depth(ny,range_err,transdepth,actual_depth, 1 bw,p,blnk,adjdepth) c write output file in rdi format with 4 variables actualdepth,u,v,err c actual depth has deepest depth in bin 1, write each ix c if (ix.eq.1) then c write(6,'(2x,a)') outname3 c endif call write_output_actual(uv,nx,ny,nv,iuv,out3, 1 ix,actual_depth,iactual_depth) c call dump_test(uv,nx,ny,nv,ix,transdepth,actual_depth) c call dmp_array(uv,nx,ny,nv) c remap velocities for bins beginning at depth new_ymin with increment new_ydel c bin 1 will now be shallowest, no extrapolation performed c shallow bins affected by side lobes are set bad if shallower than c transdepth*(1-cos(theta)) + bw+sidelobe adj if (ix.eq.1) then write(6,'(a,i5)') ' itype:',itype endif call depth_adj(uv,nx,ny,nv,ix,remap,ry,new_ymin, 1 new_ydel,actual_depth,transdepth,out2,itype,theta,bw, 1 sladj,ibin_ext,data_type) enddo !ix c write output file of remapped vels in original rdi integer format nftype = 2 write(6,'(2x,a)') outname c call dmp_array(remap,nx,ry,nv) call write_output(remap,nx,ry,nv,iuv,out,nftype, 1 new_numbins,new_ymin,new_ydel) close(out) close(out2) close(out3) 900 stop end c c ** function to compute transducer depth from pressure** c use pres=density*g*depth/10, trans_depth = depth-del_sensor c function get_trans_depth(ix,tr,lx,ly,lv,density,g,del_sensor) real get_trans_depth,density,g,del_sensor dimension tr(lx,ly,lv) get_trans_depth = tr(ix,1,1) * 10. /(g * density) - del_sensor c write(6,'(/a,f8.3,a,f8.3,a,f9.6,a,f8.5,a,f5.2)') c 1 ' transdepth:',get_trans_depth,' tr:',tr(ix,1,1),' g:', c 1 g,' dens:',density,' delsens:',del_sensor return end c *** c Subroutine to compute range error versus every 1m of vertical range from c transducer due to sound velocity adjustment only. This adjustment is c needed because RDI assumes a sound velocity of 1475.1 m/s when c converting time to meters. c c Transducer is upward looking at depth transdepth c c Subtract range_err from range in order to compute corrected range c c range_err(iy,5) has 5 variables: range,range_err,corrected_range, c meanSndVel at range, and pressure at range. subroutine vert_range_error(range_err,ctd,nyerr,nverr 1 ,transdepth,data_type) integer nyerr,nverr,iy real ctd(nyerr,nverr) integer itransdepth real last,transdepth,range_err(500,5) character*3 data_type last = 0. c compute versus vertical range from transducer rather than pressure itransdepth = transdepth + .5 c do versus 1 m increments of range rather than bins c write(6,'(/a,i5,a,f8.3/)') ' itrans:',itransdepth, c 1 ' trans:',transdepth do iy = 1,nyerr ky = itransdepth - iy if (ky.lt.1) then ky=1 endif range_err(iy,1) = iy ! range from transduder range_err(iy,2) = 1.-(ctd(ky,2)/1475.1)+last ! accumulated range err if (data_type .eq. 'vel') then range_err(iy,3) = iy - range_err(iy,2) ! corrected range elseif (data_type .eq. 'amp') then c changed feb05 range_err(iy,3) = 1 iy - range_err(iy,2) + (bw/4)*(ctd(ky,2)/1475.1) c range_err(iy,3) = iy - range_err(iy,2) + 4.5 c write(6,'(a)') ' amprangeerror' else write(6,'(a)') ' Error -- data_type must be amp or vel' go to 900 endif range_err(iy,4) = ctd(ky,2) ! sound vel at range range_err(iy,5) = ky ! Pressure at range last = range_err(iy,2) enddo 900 return end c c **** Subroutine to get actual depths of bins using corrected range array and c transducer depth, bin 1 is still deepest bin c subroutine get_depth(ny,range_err,transdepth,actual_depth, 1 bw,p,blnk,adjdepth) integer ny,iy real cor_vert_r,range_err(500,5),actual_depth(ny) real blnk,transdepth real p,bw,uncor_r c common/check/nchk,xbad c logical bad_data(40,4) c c p ~ 8.68 ! pulse length for 20 deg system c bw~ 8.68 ! Bin width for 20 deg system c blnk ~ 4.34 ! blanking for 20 deg system, P = 8.68, BW = 8.68 c n = bin number c c The uncorrected vertical range (R) to the depth cell (N): c c R = Blnk + abs(P-BW)/2 + (N*BW) c then get corrected range for each bin range from range_err array that c was computed incrementally for each meter with ratio c/1475.1 c do iy = 1,ny ! number of depths c uncor_r = blnk + abs(p-bw) + iy*bw c get even integer uncorrectd range values which bracket uncor_R value I1_uncor_R = uncor_r I2_uncor_R = uncor_r + 1. c get interpolated corrected range value cor_vert_r = ((uncor_R - I1_uncor_R)/(I2_uncor_R - I1_uncor_R) 1 * (range_err(I2_uncor_r,3) - range_err(I1_uncor_r,3))) 1 + range_err(I1_uncor_r,3) actual_depth(iy) = transdepth - cor_vert_r c 5aug03 decided to replace z (avg t_depth) with transdepth c if (actual_depth(iy).lt.z) then ! don't want neg value added to act_depth c actual_depth(iy) = actual_depth(iy)+ c 1 ((z-actual_depth(iy))/z)*adjdepth c endif c adjdepth is applied incrementally over whole water column c ie for 200m, every 2m starting from t_depth has additional .01 times adjdepth c thus, an adjustment to a given actual_depth = percentage of watercolumn from c transducer depth times adjdepth c if (actual_depth(iy).lt.transdepth) then ! don't want neg value added for pos adjdepth actual_depth(iy) = actual_depth(iy)+ 1 ((transdepth-actual_depth(iy))/transdepth)*adjdepth endif c write(6,'(1x,i3,i5,f9.3,i5,f11.3,2f9.3,f11.3)') c 1 iy,I1_uncor_r,uncor_r,I2_uncor_r, c 1 range_err(I1_uncor_r,3),cor_vert_r,range_err(I2_uncor_r,3), c 1 actual_depth(iy) enddo ! ny -- depth 990 return end c c **** Subroutine to remap RDI velocities for new depth configuration after c original bin depths corrected by avg sound vel profile and transdepth c c actual_depths still have bin 1 deepest c new depths begin at depth new_ymin with increment new_ydel c bin 1 will now be shallowest, no extrapolation of vels performed c shallow bins affected by side lobes are set bad if shallower than c transdepth*(1-cos(theta))+bw+sidelobe adj c c array is 3 dimensional in time, bin number and variable, c rdi variables are u,v,w,error c c aa is original velocity array, bb is remapped array c subroutine depth_adj(aa,nx,ny,nv,ix,bb,ry,new_ymin, 1 new_ydel,actual_depth,transdepth,out2,itype,theta,bw, 1 sladj,ibin_ext,data_type) integer nx,ny,nv,ry,iv,iy,ix,iry,iiy,itype,out2,ibin_ext dimension aa(nx,ny,nv),bb(nx,ry,nv) real mapped_depth(600),bw real vel_interp real transdepth,theta real time_int real save_depth_iy,save_depth_iyp1 real sh_map_d,save_depth_up real actual_depth(ny),new_ymin,new_ydel common/check/nchk,xbad logical bad_data(44,4),not_found character*40 form COMMON /INFORM/ FORM,IDATE(5),LDATE(5) COMMON/matrix/ ImatDAT(5),LmatDAT(5),XDEL character*3 data_type c sladj is sidelobe adjustment, input to deepen cutoff depth for sidelobe error c logical not_found used when looking for two shallowest depths that bracket c the first bin. May be same if exactly same value as data depth c depth above which there is surface sidelobe interference c = D*(1-cos(theta)+bw+sidelobe adj c where D is depth of transducer c write(6,'(2(a,f10.5))') 'theta:',theta,' bw:',bw c modified 5 mar 02 to allow input to deeper cutoff depth for sidelobe error if (itype.eq.1) then data_good_depth = -500. elseif (itype.eq.2) then if (data_type.eq.'vel') then data_good_depth = (transdepth * (1.-cosd(theta)) ) +bw+sladj elseif (data_type.eq.'amp') then c amp mapped shallower because taken at end of pulse c note using approximate nominal bw here rather than actual bw data_good_depth = (transdepth * (1.-cosd(theta)) ) +.75*bw+sladj else write(6,'(a)') 'Error: must be vel or amp' endif else write(6,'(a)') ' Error* itype must be between 1 or 2' go to 990 endif c if (ix.eq.2.or.ix.eq.100.or.ix.eq.400.or.ix.eq.1000) then c write(6,'(/a,i5,a,f8.2,a,f10.4/)') c 1 ' ix:',ix,' transdepth:',transdepth, c 1 ' data_good_depth:',data_good_depth c endif do iy = 1,ny if (actual_depth(iy).lt.data_good_depth) then do iv = 1,nv aa(ix,iy,iv) = xbad enddo c if (ix.eq.1) then c write(6,'(a,i4,2(a,f8.2))') ' iy:',iy, c 1 ' act_d:',actual_depth(iy), c 1 ' d_gd_d:',data_good_depth c endif endif do iv = 1,nv bad_data(iy,iv) = .false. if (aa(ix,iy,iv).eq.xbad) then bad_data(iy,iv) = .true. endif enddo ! nv enddo ! ny c extrapolate to two depths shallower, shallowest depth largest iy do iy = 3,ny-2 if(actual_depth(iy).lt.data_good_depth.and. 1 actual_depth(iy-1).ge.data_good_depth.and. 1 actual_depth(iy-2).ge.data_good_depth) then c ***changed this 29 jul 04 (iy-2) above not (iy-1) c ***1 actual_depth(iy-1).ge.data_good_depth) then do iv = 1,nv c changed to add line below 4nov04 so don't extrapolate with 5120 if depths c ok but data bad. Should come up with a way to search down for consecutive gd data c but later... if(bad_data(iy-1,iv).or.bad_data(iy-2,iv)) go to 222 aa(ix,iy,iv) = vel_interp(aa(ix,iy-1,iv),aa(ix,iy-2,iv), 1 actual_depth(iy),actual_depth(iy-1),actual_depth(iy-2)) c write(6,'(2x,2i4,3(f8.2,f8.3,3x))') c 1 ix,iy,actual_depth(iy),aa(ix,iy,iv),actual_depth(iy-1), c 1 aa(ix,iy-1,iv),actual_depth(iy-2),aa(ix,iy-2,iv) bad_data(iy,iv) = .false. c next shallower depth if (ibin_ext.lt.2) go to 222 aa(ix,iy+1,iv) = vel_interp(aa(ix,iy,iv),aa(ix,iy-1,iv), 1 actual_depth(iy+1),actual_depth(iy),actual_depth(iy-1)) c write(6,'(2x,2i4,3(f8.2,f8.3,3x))') ix,iy+1, c 1 actual_depth(iy+1),aa(ix,iy+1,iv),actual_depth(iy), c 1 aa(ix,iy,iv),actual_depth(iy-1),aa(ix,iy-1,iv) bad_data(iy+1,iv) = .false. c next shallower depth if (ibin_ext.lt.3) go to 222 aa(ix,iy+2,iv) = vel_interp(aa(ix,iy+1,iv),aa(ix,iy,iv), 1 actual_depth(iy+2),actual_depth(iy+1),actual_depth(iy)) c write(6,'(2x,2i4,3(f8.2,f8.3,3x))') ix,iy+2, c 1 actual_depth(iy+2),aa(ix,iy+2,iv),actual_depth(iy+1), c 1 aa(ix,iy+1,iv),actual_depth(iy),aa(ix,iy,iv) bad_data(iy+2,iv) = .false. c next shallower depth added (not in remap_subdata_extrapolate) if (ibin_ext.lt.4) go to 222 aa(ix,iy+3,iv) = vel_interp(aa(ix,iy+2,iv),aa(ix,iy+1,iv), 1 actual_depth(iy+3),actual_depth(iy+2),actual_depth(iy+1)) if (ibin_ext.gt.4) then write(6,'(a)') ' Error: # bins to extrapolate must be 1 to 4' endif 222 continue enddo !nv endif enddo c Compute corrected velocities not_found = .TRUE. do iry = 1,ry mapped_depth(iry) = new_ymin+(iry-1)*new_ydel !new remapped depth values c write(6,'(a,i4,a,f8.2)') c 1 ' iry:',iry,' mapped_d:',mapped_depth(iry) iy = 2 801 continue if (mapped_depth(iry).GT.actual_depth(iy-1)) then do iv = 1,nv bb(ix,iry,iv) = xbad enddo ! nv elseif (mapped_depth(iry).eq.actual_depth(iy-1)) then do iv = 1,nv if (bad_data(iy-1,iv)) then bb(ix,iry,iv)=xbad else bb(ix,iry,iv)=aa(ix,iy-1,iv) if (not_found) then save_depth_iy=actual_depth(iy-1) save_depth_iyp1=actual_depth(iy-1) sh_map_d=mapped_depth(iry) save_depth_up=actual_depth(iy) endif not_found = .FALSE. endif enddo elseif (mapped_depth(iry).LT.actual_depth(iy-1).and. 1 mapped_depth(iry).GT.actual_depth(iy)) then do iv = 1,nv if ( (bad_data(iy,iv)).or.(bad_data(iy-1,iv)) ) then bb(ix,iry,iv) = xbad else ! good data both sides, interpolate bb(ix,iry,iv) = vel_interp(aa(ix,iy-1,iv),aa(ix,iy,iv), 1 mapped_depth(iry),actual_depth(iy-1),actual_depth(iy)) if (not_found) then save_depth_iy=actual_depth(iy-1) save_depth_iyp1=actual_depth(iy) sh_map_d=mapped_depth(iry) if (iy.eq.ny) then save_depth_up=actual_depth(iy) elseif (iy.lt.ny) then save_depth_up=actual_depth(iy+1) endif endif not_found = .FALSE. endif enddo !nv variables elseif (mapped_depth(iry).eq.actual_depth(iy)) then do iv = 1,nv if (bad_data(iy,iv)) then bb(ix,iry,iv)=xbad else bb(ix,iry,iv)=aa(ix,iy,iv) if (not_found) then save_depth_iy=actual_depth(iy) save_depth_iyp1=actual_depth(iy) sh_map_d=mapped_depth(iry) if (iy.eq.ny) then save_depth_up=actual_depth(iy) elseif (iy.lt.ny) then save_depth_up=actual_depth(iy+1) endif endif not_found = .FALSE. endif enddo elseif (mapped_depth(iry).lt.actual_depth(iy)) then iy = iy+1 if (iy.le.ny) then go to 801 else do iv = 1,nv bb(ix,iry,iv) = xbad enddo !nv endif else write(6,'(a)') ' Error in depth comparisons' go to 899 endif enddo ! iry, remapped depth c c Get time increment (XDEL in days) and convert to time_int (in minutes) c Then recreate date from getdt2 subroutine by incrementing time (ix) c Reconvert back to RDI time format c if (XDEL.eq.1) then time_int = 1440. elseif (XDEL.ge..04.and.XDEL.lt..05) then time_int = 60. elseif (XDEL.ge..02.and.XDEL.lt..03) then time_int = 30. elseif (XDEL.ge..01.and.XDEL.lt..02) then time_int = 15. else write(6,'(a)') ' Error: time increment not a day, an hour', 1 ' 30 or 15 minutes' go to 990 endif call getdt2(ImatDAT,1,time_int,ix,Ldate) if (save_depth_iyp1.lt.0) then save_depth_iyp1 = 0. endif if (save_depth_iy.lt.0) then save_depth_iy = 0. endif if (save_depth_up.lt.0) then save_depth_up = 0. endif c added test numbers to help with plotting straight lines at given depth over c vel plot test1 = 1. test2 = 2. write(out2,'(1x,2i2,1x,3i2,4x,2f8.2,4f10.2,2f3.0)') 1 (ldate(ii),ii=1,5),transdepth,data_good_depth, 1 sh_map_d,save_depth_iyp1,save_depth_iy, 1 save_depth_up,test1,test2 if (ix.eq.1) then write(6,'(/a)') ' Mapped depths ' do iry = 1,ry,8 if ((iry+7).le.ry) then write(6,'(8x,8f8.2)') (mapped_depth(iiy),iiy=iry,iry+7) else write(6,'(8x,8f8.2)') (mapped_depth(iiy),iiy=iry,ry) endif enddo endif 899 continue 990 return end c c *** function to extrapolate or interpolate between vels at two old depths to c a new vel at a new depth. c x1,x2 and y1,y2 are old vels and depths. y is new depth. function vel_interp(x1,x2,y,y1,y2) real x1,x2,y,y1,y2,vel_interp vel_interp = ( (y-y1)/(y2-y1)*(x2-x1) ) + x1 return end c c **** Subroutine to write output file ***** c subroutine write_output(remap,nx,ry,nv,iuv,out,nftype, 1 new_numbins,new_ymin,new_ydel) integer nx,ry,nv,ix,iy,iv,out,nftype,ntype integer ndate(5) common/peprdi/nrec,NBins,NBout,NBstart,NBend character*40 form COMMON/NAME/PROJ(2),MOOR,VNAM,ID,DELTin,DEPTH,REST COMMON /INFORM/ FORM,IDATE(5),LDATE(5) real new_ymin,new_ydel real time_int COMMON/matrix/ ImatDAT(5),LmatDAT(5),XDEL dimension remap(nx,ry,nv) dimension iuv(nx,ry,nv) common/check/nchk,xbad integer new_numbins,nb c call dmp_array(uv,nx,ny,nv) c call dmp_array(remap,nx,ry,nv) c c set ntype to be negative 1 or negative 2 for calling inverse of cnv_date c that is convert from convertional time back to RDI time if (nftype.lt.0) then ntype = nftype else ntype = - nftype endif write(6,'(/)') if (nv.eq.1) then form = '(10(1x,i6))' write(6,'(a,a)') 1 ' The 1 variable RDI output file format = ',form elseif (nv.eq.2) then form = '(5(1x,2i6))' write(6,'(a,a)') 1 ' The 2 variable RDI output file format = ',form elseif (nv.eq.3) then form = '(4(1x,3i6))' write(6,'(a,a)') 1 ' The 3 variable RDI output file format = ',form elseif (nv.eq.4) then form = '(3(1x,4i6))' write(6,'(a,a)') 1 ' The 4 variable RDI output file format = ',form endif nb=index(form,'))')+1 ! assumes standard rdi memory format (*(*x,*i*)) do 335 ix = 1,nx ! time c c Get time increment (XDEL in days) and convert to time_int (in minutes) c Then recreate date from getdt2 subroutine by incrementing time (ix) c Reconvert back to RDI time format c if (XDEL.eq.1) then time_int = 1440. elseif (XDEL.ge..04.and.XDEL.lt..05) then time_int = 60. elseif (XDEL.ge..02.and.XDEL.lt..03) then time_int = 30. elseif (XDEL.ge..01.and.XDEL.lt..02) then time_int = 15. else write(6,'(a)') ' Error: time increment not a day, an hour', 1 ' 30 or 15 minutes' go to 99 endif call getdt2(ImatDAT,1,time_int,ix,Ldate) call cnv_date(Ldate,ndate,ntype) ! ntype negative reverses conversion c c Write out header for each time increment (ix) c c !old format ** write(out,229) ix,ndate,NBins,NBout,NBStart,NBend,form write(out,229) ix,ndate,NBins,new_numbins,NBStart,new_numbins, 1 form(1:nb),' !',nv,new_ymin,new_ydel 229 format(I6,5I3,4I4,2X,A,a,i2,2x,f6.1,f5.1) do 47 iy = 1,ry do 48 iv = 1,nv if(remap(ix,iy,iv).le.xbad) then ! xbad is set to -999 in readarray iuv(ix,iy,iv) = 5120 else if (remap(ix,iy,iv).lt.0.) then ! .05 for integer roundoff iuv(ix,iy,iv) = (remap(ix,iy,iv)-.05) * 10. else iuv(ix,iy,iv) = (remap(ix,iy,iv)+.05) * 10. endif endif 48 enddo !iv 47 enddo !iy write(out,form) ((iuv(ix,iy,iv),iv=1,nv),iy=1,ry) c write out header to log file for first two and last ensembles if (ix.le.2.or.ix.eq.100.or.ix.eq.400.or.ix.eq.nx) then write(6,'(/)') write(6,229) ix,ndate,NBins,new_numbins,NBStart,new_numbins, 1 form(1:nb),' !',nv,new_ymin,new_ydel write(6,form) ((iuv(ix,iy,iv),iv=1,nv),iy=1,ry) endif 335 enddo ! ix 99 return end c c **** Subroutine to dump test values ***** c subroutine dump_test(uv,nx,ny,nv,ix,transdepth,actual_depth) dimension uv(nx,ny,nv) real actual_depth(ny),transdepth write(6,'(a/)') 1 ' bin_no actual_depth u v w err' do iy = 1,ny write(6,'(1x,i6,f14.4,4f8.2)') 1 iy,actual_depth(iy),(uv(ix,iy,iv),iv=1,4) enddo !iy return end c c **** Subroutine to dump array values ***** c subroutine dmp_array(uv,nx,ny,nv) dimension uv(nx,ny,nv) do iv=1,nv write(6,101) iv 101 format(//10x,'Variable ',i1) do iy=1,ny write(6,102) (uv(ix,iy,iv),ix=1,nx) 102 format(1x,15f8.2) end do end do return end c c **** Subroutine to write output file of actual depths and uv ***** c subroutine write_output_actual(uv,nx,ny,nv,iuv,out3, 1 ix,actual_depth,iactual_depth) integer nx,ny,nv,ix,iy,iv,out3,nb integer ndate(5) integer*4 iactual_depth(ny) common/peprdi/nrec,NBins,NBout,NBstart,NBend character*40 form COMMON/NAME/PROJ(2),MOOR,VNAM,ID,DELTin,DEPTH,REST COMMON /INFORM/ FORM,IDATE(5),LDATE(5) real time_int COMMON/matrix/ ImatDAT(5),LmatDAT(5),XDEL dimension uv(nx,ny,nv) real actual_depth(ny) dimension iuv(nx,ny,nv) common/check/nchk,xbad form = '(3(1x,4i6))' c form = '(2(1x,5i6))' if (ix.eq.1) then write(6,'(/)') write(6,'(a/a,a)') 1 ' The 4 variables are depth u v err', 1 ' The 4 variable RDI output file format = ',form endif nb=index(form,'))')+1 ! assumes standard rdi memory format (*(*x,*i*)) c c Get time increment (XDEL in days) and convert to time_int (in minutes) c Then recreate date from getdt2 subroutine by incrementing time (ix) c Reconvert back to RDI time format c if (XDEL.eq.1) then time_int = 1440. elseif (XDEL.ge..04.and.XDEL.lt..05) then time_int = 60. elseif (XDEL.ge..02.and.XDEL.lt..03) then time_int = 30. elseif (XDEL.ge..01.and.XDEL.lt..02) then time_int = 15. else write(6,'(a)') ' Error: time increment not a day, an hour', 1 '30 or 15 minutes' go to 99 endif call getdt2(ImatDAT,1,time_int,ix,Ldate) call cnv_date(Ldate,ndate,-2) ! ntype negative reverses conversion c c Write out header for each time increment (ix) c write(out3,229) ix,ndate,NBins,NBins,NBStart,Nbins, 1 form(1:nb),' !',nv 229 format(I6,5I3,4I4,2X,A,a,i2) c write out header to log file for first two and last ensembles if (ix.le.2.or.ix.eq.100.or.ix.eq.400.or.ix.eq.nx) then write(6,'(/)') write(6,229) ix,ndate,NBins,Nbins,NBStart,Nbins, 1 form(1:nb),' !',nv endif do 22 iy = 1,ny do 32 iv = 1,nv if(uv(ix,iy,iv).le.xbad) then ! xbad is set to -999 in readarray iuv(ix,iy,iv) = 5120 else if (uv(ix,iy,iv).lt.0.) then ! .05 for integer roundoff iuv(ix,iy,iv) = (uv(ix,iy,iv)-.05) * 10. else iuv(ix,iy,iv) = (uv(ix,iy,iv)+.05) * 10. endif endif 32 enddo !iv if (actual_depth(iy).lt.0.) then iactual_depth(iy)=(actual_depth(iy)-.05) * 10. else iactual_depth(iy)=(actual_depth(iy)+.05) * 10. endif 22 enddo !iy c changed 8/8/96 to write actual_depth,u,v,err rather than actual_depth,u,v,w write(out3,form) 1 (iactual_depth(iy),iuv(ix,iy,1),iuv(ix,iy,2), 1 iuv(ix,iy,4),iy=1,ny) c write to log file first two and last ensembles if (ix.le.2.or.ix.eq.100.or.ix.eq.400.or.ix.eq.nx) then write(6,form) 1 (iactual_depth(iy),iuv(ix,iy,1),iuv(ix,iy,2), 1 iuv(ix,iy,4),iy=1,ny) endif 99 return end