program momcor c c this program reads the ascii moments from c ship profiler and spits out moments that c are corrected for mean ship motion. c THE STEPS FOR CORRECTION: c 1) rotate oblique radials to a N/E ref. frame c 2) rotate the ship speed vector into oblique c and then a N/E ref. frame c 3) rotate vertical radial into olbique and then c N/E ref. frame c 4) subtract off vertical correction that c will be applied later to correct c wind speed for vertical motion c implicit real*8 (a-h,o-z) real*8 rad(5000,110),snr(5000,110),pow(5000,110), x wid(5000,110),pnois(5000,110),rn(110), x re(110),spd(5000),hed(5000),crs(5000) integer ih(5000,7),mh(5000,8),ibm(5000) character st*3,yr*2,jday*3,rdlin*55,mode*1 pi = 3.141592654 c open(unit=1,file='mom.inp') read(1,2) st,yr 2 format(a3,1x,a2) read(1,*) ibeg read(1,*) read(1,*) azicor,ibma1,ibmp1,ibmv1,ibma2,ibmp2,ibmv2 close(unit=1) c c read the profiler moments c write(jday(1:3),7) ibeg 7 format(i3.3) open(unit=1,file=st//yr//jday//'.raw') read(1,8) rdlin 8 format(a55) do 199 l = 1,5000 read(1,*,end=200) ibm(l),(ih(l,m),m=1,7) read(1,*,end=200) (mh(l,m),m=1,8) do 97 n = 1,mh(l,5) read(1,*,end=200) rad(l,n),snr(l,n), x pow(l,n),wid(l,n),pnois(l,n) 97 continue 199 continue 200 lct = l-1 close(unit=1) c close(unit=1,status='delete') c c output file c open(unit=3,file='.\corrected\'//st//yr//jday//'.raw') write(3,8) rdlin c c read the ship motion data c open(unit=1,file='shipcr'//jday//'.dat') do 299 l = 1,lct read(1,*) nbm,ihr,imin,isec,spd(l),hed(l), x crs(l) 299 continue c close(unit=1,status='delete') close(unit=1) c c correct moments: make sure to start with beam c index 0 or 3 and ensure a moment triplet c (use index 2 and 4 for AOE2001 cruise) c luse = 1 c c loop thru moment triplets c 333 l = luse c c determine whether low mode or high mode and assign beams c if(ibm(l).eq.2) then ibma = ibma1 ibmp = ibmp1 ibmv = ibmv1 mode = 'l' endif if(ibm(l).eq.4) then ibma = ibma2 ibmp = ibmp2 ibmv = ibmv2 mode = 'h' endif if(ibm(l).ne.2.and.ibm(l).ne.4) then luse = luse+1 go to 333 endif c print *, l,ibm(l),ibma,ibmp,ibmv c c check for end of file if(l+2.gt.lct) go to 600 c check to see that beams are sequential c if not, then reset luse c if(ibm(l+1).ne.(ibm(l)+1)) then if(mode.eq.'l'.and.(ibm(l+1).ne.0.and.ibm(l+2).ne.1)) then luse = l+1 go to 333 endif c if(ibm(l+2).ne.(ibm(l)+2)) then if(mode.eq.'h'.and.(ibm(l+1).ne.3.and.ibm(l+2).ne.5)) then luse = l+2 go to 333 endif c average speed,heading,course over triplet if(spd(l).lt.-98.or.(spd(l+1).lt.-98.or.spd(l+2).lt.-98)) x go to 222 if(hed(l).lt.-98.or.(hed(l+1).lt.-98.or.hed(l+2).lt.-98)) x go to 222 if(crs(l).lt.-98.or.(crs(l+1).lt.-98.or.crs(l+2).lt.-98)) x go to 222 cn1 = spd(l)*dcos(crs(l)*pi/180.) cn2 = spd(l+1)*dcos(crs(l+1)*pi/180.) cn3 = spd(l+2)*dcos(crs(l+2)*pi/180.) ce1 = spd(l)*dsin(crs(l)*pi/180.) ce2 = spd(l+1)*dsin(crs(l+1)*pi/180.) ce3 = spd(l+2)*dsin(crs(l+2)*pi/180.) acn = (cn1+cn2+cn3)/3. ace = (ce1+ce2+ce3)/3. hn1 = 1.*dcos(hed(l)*pi/180.) hn2 = 1.*dcos(hed(l+1)*pi/180.) hn3 = 1.*dcos(hed(l+2)*pi/180.) he1 = 1.*dsin(hed(l)*pi/180.) he2 = 1.*dsin(hed(l+1)*pi/180.) he3 = 1.*dsin(hed(l+2)*pi/180.) ahn = (hn1+hn2+hn3)/3. ahe = (he1+he2+he3)/3. us = (acn*acn+ace*ace)**0.5 c = datan2(ace,acn) h = datan2(ahe,ahn) c define beam angles used in correction formulae if(mode.eq.'l') then tha = (float(ih(l+2,2))+azicor)*pi/180. thp = (float(ih(l+1,2))+azicor)*pi/180. z = float(ih(l+1,1))*pi/180. endif if(mode.eq.'h') then tha = (float(ih(l,2))+azicor)*pi/180. thp = (float(ih(l+1+ibmp,2))+azicor)*pi/180. z = float(ih(l+1+ibma,1))*pi/180. endif c gate loop c do 499 n = 1,mh(l,5) c if(mode.eq.'l') then ra = rad(l+2,n) rp = rad(l+1,n) rv = rad(l,n) endif if(mode.eq.'h') then ra = rad(l,n) rp = rad(l+1,n) rv = rad(l+2,n) endif c north radial rn(n) = ra*dcos(tha+h)-rp*dsin(tha+h) x -us*dcos(c)*dcos(z) x -(rv*dsin(z)*dcos(tha+h) x -rv*dsin(z)*dsin(tha+h)-rv*dsin(z)) c east radial re(n) = ra*dsin(tha+h)+rp*dcos(tha+h) x -us*dsin(c)*dcos(z) x -(rv*dsin(z)*dsin(tha+h) x +rv*dsin(z)*dcos(tha+h)-rv*dsin(z)) 499 continue c c write corrected moments to file c for nahru99 only: vert, port, aft, aft, port, vert c then must use special version of qcave.f c use original vertical beam data c if(mode.eq.'l') then write(3,412) ibm(l),(ih(l,m),m=1,7) write(3,413) (mh(l,m),m=1,8) do 439 n = 1,mh(l,5) write(3,414) rad(l,n),snr(l,n), x pow(l,n),wid(l,n),pnois(l,n) 439 continue write(3,411) ibm(l+1),ih(l+1,1),(ih(l+1,m),m=3,7) write(3,413) (mh(l+1,m),m=1,8) do 449 n = 1,mh(l+1,5) write(3,414) re(n),snr(l+1,n),pow(l+1,n), x wid(l+1,n),pnois(l+1,n) 449 continue write(3,410) ibm(l+2),ih(l+2,1),(ih(l+2,m),m=3,7) write(3,413) (mh(l+2,m),m=1,8) do 459 n = 1,mh(l+2,5) write(3,414) rn(n),snr(l+2,n),pow(l+2,n), x wid(l+2,n),pnois(l+2,n) 459 continue endif c if(mode.eq.'h') then write(3,410) ibm(l),ih(l,1),(ih(l,m),m=3,7) write(3,413) (mh(l,m),m=1,8) do 469 n = 1,mh(l,5) write(3,414) rn(n),snr(l,n),pow(l,n), x wid(l,n),pnois(l,n) 469 continue write(3,411) ibm(l+1),ih(l+1,1),(ih(l+1,m),m=3,7) write(3,413) (mh(l+1,m),m=1,8) do 479 n = 1,mh(l+1,5) write(3,414) re(n),snr(l+1,n),pow(l+1,n), x wid(l+1,n),pnois(l+1,n) 479 continue write(3,412) ibm(l+2),(ih(l+2,m),m=1,7) write(3,413) (mh(l+2,m),m=1,8) do 489 n = 1,mh(l+2,5) write(3,414) rad(l+2,n),snr(l+2,n), x pow(l+2,n),wid(l+2,n),pnois(l+2,n) 489 continue endif c 410 format(i2,1x,i2,' 0 ',i2,1x,i3,3(1x,i2)) 411 format(i2,1x,i2,' 90 ',i2,1x,i3,3(1x,i2)) 412 format(i2,1x,i2,1x,i3,1x,i2,1x,i3,3(1x,i2)) 413 format(i3,1x,i4,4(1x,i3),1x,i2,1x,i4) 414 format(5(1x,f7.3),1x,'0') c c if no valid ship data 222 continue c c c update luse and continue thru moments c luse = luse+1 go to 333 c c end loop 600 continue c close(unit=3) c stop end