program reorder character*80 filename integer passmjd(4000,2),idata(400,2),istore(2), > countall,jstop,choice,passno(4000),dndk, > pntcnt(4000,2),shrtpas(3000,2),shrtcnt, > passrec(4000,2),passrem,cnter,denum,pchoice, > spknum,mincheck,spkvar,nowant(4000),outnum, > minchk,cntsome,inrec,innum,outrec,oif,inf,onf, > innumall,pntall(4000,2) real data(400,27),dstore(27),east,west,diffwe,hilat, > lolat,north,south,percent,totlat,upper,lower, > desdata(400,27),intpdata(400,27) double precision aver(4000,2),ra(4000),cross,passavg(4000,2), > savglon(4000,2),crosss common data(400,27),idata(400,2) common /order1/ aver(4000,2),passmjd(4000,2),cross common /order2/ passno(4000),pntcnt(4000,2),pntall(4000,2) common /order3/ savglon(4000,2),crosss common /hsort/ ra(4000) common /shorty/ shrtpas(3000,2),shrtcnt,passrec(4000,2) common /spike/ desdata(400,27),upper,lower,spkvar common /intb/ intpdata(400,27) c c----------------------------------- program description c c this program takes data in the 2-integers and 27-reals format c and reorders the entire dataset into a sorted file according c to the variable that the user chooses. that variable is usually c the average longitude of each individual pass, such that after c reordering, the dataset will have all passes arranged from LOWEST c average longitude (-179.99) to the HIGHEST average longitude c (+179.99). if the dataset crosses the -180.0 180.0 meridian, c then the eastern (negative) longitudes are incremented to a c positive value by adding 360.0 (see write statement with c variable "cross"). the sorting variable can also be the c average elevation or the pass numbers of each individual pass c (sorting by pass numbers = sorting by time). this program takes c a little time (about 15 minutes on the DECstation 3100). the c program requires DIRECT ACCESS or else it just won't happen! c NOTE: as crazy as it may seem, real*8 is necessary for the c averages because i found two dusk longitude averages to be c EXACTLY the same with real*4. if two averages are the same c then passno(4000) will have the same pass twice which messes c up subroutine reorder2. c NOTE: i usually try to keep all file reads and writes in the c main program but reorder2 and reorder3 are a deviation c from the rule c NOTE: the dataset must be despiked and interpolated to c correctly calculate the longitude averages of the c extended passes. however, i prefer to not write out c the despiked or interpolated data because this program c represents the end of the first processing step, after- c which the data is ready for more involved processing c (ie. correlation filtering, bandpassing ...). c therefore, the output from reorder should be the c original data, only reordered. get it?? c NOTE: for direct access on an ibm rs6000, recl=116. c on a dec3100, recl=29 c c program date: 16 apr 91 c write (*,*) 'INPUT 2I-27R FILE:' read (*,9990) filename 9990 format (a80) open (10, file=filename, status='old',form='unformatted', > access='direct',recl=116) write (*,*) '0 IF THIS IS A DUSK DATASET' write (*,*) '1 IF THIS IS A DAWN DATASET' read (*,*) dndk write (*,*) '1 TO AVERAGE AND REORDER ON LONGITUDE' write (*,*) '2 TO AVERAGE AND REORDER ELEVATION' write (*,*) '3 TO REORDER ON PASS NUMBER' read (*,*) choice c if (choice .eq. 1) then write (*,*) 'OUTPUT FILE OF 2I-27R DATASET REORDERED' read (*,9990) filename open (20, file=filename,form='unformatted',access='DIRECT', > recl=116) write (*,*) 'INTERMEDIATE I/0 FILE NOT REORDERED' write (*,*) '------ DO NOT USE THIS FILE -------' read (*,9990) filename open (21, file=filename,form='unformatted',access='direct', > recl=116) endif if (choice .gt. 1) then write (*,*) 'OUTPUT FILE OF 2I-27R DATASET REORDERED' read (*,9990) filename open (20, file=filename,form='unformatted',access='direct', > recl=116) endif write (*,*) 'OUTPUT FILE OF PASS NUMBERS AND AVERAGED', > ' SORTED VARIABLE' read (*,9990) filename open (22, file=filename,form='formatted') c write (*,*) '0 IF YOU WANT ALL PASSES' write (*,*) '1 IF SOME PASSES NEED TO BE REMOVED' read (*,*) pchoice if (pchoice .eq. 1) then write (*,*) 'INPUT FILE OF PASSES YOU DO NOT WANT' read (*,9990) filename open (11, file=filename,form='formatted') endif write (*,*) 'WHAT IS THE MINIMUM NUMBER OF' write (*,*) 'OBSERVATIONS ALLOWABLE FOR EACH PASS (50)' read (*,*) mincheck c write (*,*) '0 FOR NO DESPIKING OF DATA SET' write (*,*) '1 FOR DESPIKING ONCE' write (*,*) '2 FOR DESPIKING TWICE (this is the usual choice)' write (*,*) '3 ... AND SO ON' read (*,*) spknum if (spknum .gt. 0) then write (*,*) 'WHAT IS THE MAXIMUM nT: (1.0)' write (*,*) 'WHAT IS THE MINIMUM nT: (-1.0)' read (*,*) upper,lower write (*,*) 'WHICH VARIABLE TO DESPIKE: (23)' write (*,*) '1=LAT, 2=LONG,...23=totavgmag..25=resavgmag...' write (*,*) 'lat lon rad mlt invlat diplat bs bv x y z' write (*,*) 'bva xa ya za totfld xfld yfld zfld inc dec' write (*,*) 'totmag totavgmag resid resavgmag ringcur sec ' read (*,*) spkvar endif c---------------------------- the following if statement determines if c the study area includes the -180.0 180.0 c longitude line. for further comments see c subroutine reorder1 cross=0.0 crosss=cross minchk=0 if (choice .eq. 1) then write (*,*) 'WESTERN MOST LONGITUDE OF STUDY AREA' write (*,*) 'EASTERN MOST LONGITUDE OF STUDY AREA' write (*,*) '-180.0 to 180.0 NOT 0.0 to 360.0' read (*,*) west,east diffwe = west - east if (diffwe .gt. 0.0) then cross=360.0 crosss=cross write (*,*) ' ' write (*,*) 'the program has determined that this study' write (*,*) 'area crosses the 180.0, -180.0 meridian' write (*,*) ' ' endif write (*,*) 'NORTHERN AND SOUTHERN MOST LATITUDES' write (*,*) '90.0 to -90.0 NOT 0.0 to 180.0' read (*,*) north,south write (*,*) 'PERCENT OF TOTAL LATITUDE LENGTH TO' write (*,*) 'TO BE CONSIDERED TO FIND SHORT PASSES (90.0)' read (*,*) percent c-------------------------- percent is used to calculate the range that c is used in subroutine shorts to determine c if a pass is a short pass or a long pass. c see shorts for more info. also, since no c passes go below or above 83.0 degrees, the c program resets north and south if needed. if (north .gt. 83.0) north = 83.0 if (south .lt. -83.0) south = -83.0 totlat=abs(north-south) minchk=(int(((100.0-percent)/(2*100.0))*totlat)+1)*3 percent=((100.0-percent)/(2*100.0))*totlat lolat=south+percent hilat=north-percent endif c if (minchk .gt. mincheck) then mincheck=minchk write (*,*) 'minimum observation cut-off increased to', > ' =',mincheck endif c if (pchoice .eq. 1) then do 50 kk=1,4000 read (11,*,end=55) nowant(kk) 50 continue 55 continue pchoice=kk-1 endif c c--------------------------- the main program reads the data to find c which 2i-27r lines belong to a specific pass c number (idata(n,1)) and since it reads one line c of the next pass it stores that line in c memory. c inrec=1 outrec=0 shrtcnt=0 countall=0 jstop=0 passrem=0 cntsome=0 c read (10,rec=inrec) (idata(1,i),i=1,2),(data(1,j),j=1,27) 100 n=2 105 inrec=inrec+1 read (10,rec=inrec,err=110) (idata(n,i),i=1,2),(data(n,j),j=1,27) if (idata(n,1) .ne. idata(n-1,1)) go to 120 n=n+1 go to 105 110 continue jstop=1 120 continue do 130 i=1,2 istore(i)=idata(n,i) 130 continue do 140 i=1,27 dstore(i)=data(n,i) 140 continue c countall=countall+1 innum=n-1 innumall=innum c c------------------------------------ if passes are NOT wanted, remove them if (pchoice .eq. 0) go to 145 do 143 ii=1,pchoice if (idata(innum,1) .eq. nowant(ii)) then write (*,*) 'PASS NUMBER REMOVED ',nowant(ii), > idata(innum,1) passrem=passrem+1 go to 400 endif 143 continue 145 continue c if (innum .lt. mincheck) then write (*,9980) idata(innum,1),innum 9980 format ('PASS REMOVED AT READ =',i6,' OBSERV COUNT =',i5) passrem=passrem+1 go to 400 endif c------------------------------------ search for passes that cross from c -180.0 to 180.0 meridian imerpass=idata(1,1) call meridian (innum,imerpass) c------------------------------------ despike the data if user chooses and c despike the number of times chosen if (spknum .eq. 0) go to 190 cnter=0 150 call despike (innum,denum) cnter=cnter+1 innum=denum do 180 k=1,innum do 180 kk=1,27 data(k,kk)=desdata(k,kk) 180 continue if (cnter .lt. spknum) go to 150 c 190 continue if (innum .lt. mincheck) then write (*,9981) idata(innum,1),innum 9981 format ('PASS REMOVED AFTER DESPIKING =',i6, > ' OBSERV COUNT =',i5) passrem=passrem+1 go to 400 endif c c------------------------------------- interpolate the dataset call interp1 (dndk,innum,outnum) innum=outnum do 210 i=1,innum do 210 j=1,27 data(i,j)=intpdata(i,j) 210 continue if (innum .lt. mincheck) then write (*,9982) idata(1,1),innum 9982 format ('PASS REMOVED AFTER INTERPOLATING =',i6, > ' OBSERV COUNT =',i5) passrem=passrem+1 go to 400 endif c if (choice .eq. 1) oif=21 if (choice .gt. 1) go to 260 do 250 i=1,innum outrec=outrec+1 write (oif,rec=outrec) (idata(1,j),j=1,2),(data(i,j),j=1,27) 250 continue c 260 continue cntsome=cntsome+1 pntcnt(cntsome,2)=innum pntcnt(cntsome,1)=idata(1,1) pntall(cntsome,1)=idata(1,1) pntall(cntsome,2)=innumall c c------------------------------------- subroutine finds all short passes if (choice .eq. 1) call shorts (innum,hilat,lolat,dndk) c c------------------------------------- now call reorder1 which will find c the average longitude and elevation c (not radius) and store them in c aver(4000) as well as storing the c pass number and modified julian day for c the current pass. call reorder1 (innum,cntsome) c------------------------------------- ok, now go back and get more passes c to average until done with the file c 400 continue do 410 i=1,2 idata(1,i)=istore(i) 410 continue do 420 i=1,27 data(1,i)=dstore(i) 420 continue if (jstop .eq. 1) go to 500 go to 100 c 500 continue c c------------------------------------ now sort the chosen variable call sort (cntsome,choice) c do 550 i=1,cntsome do 530 j=1,cntsome if (choice .eq. 1) then if (aver(j,choice) .eq. ra(i)) then passno(i)=passmjd(j,1) passavg(i,1)=dble(passmjd(j,1)) passavg(i,2)=aver(j,choice) write (22,*) (passmjd(j,ii),ii=1,2), > (aver(j,ii),ii=1,2) go to 540 endif elseif (choice .eq. 2) then if (aver(j,choice) .eq. ra(i)) then write (22,*) (passmjd(j,ii),ii=1,2), > (aver(j,ii),ii=1,2) passno(i)=passmjd(j,1) go to 540 endif elseif (choice .eq. 3) then if (passmjd(j,1) .eq. int(ra(i))) then write (22,*) (passmjd(j,ii),ii=1,2), > (aver(j,ii),ii=1,2) passno(i)=passmjd(j,1) go to 540 endif endif 530 continue 540 continue 550 continue c c------------------------------ now reread the file in reorder2 and write c it ordered according to the pass numbers c given to reorder2. c if (choice .eq. 1) then inf=21 onf=20 elseif (choice .gt. 1) then inf=10 onf=20 endif call reorder2 (cntsome,inf,onf) c c------------------------------ if sorting by average longitude, then must c extend the shorter passes and calculate a c new average. see subroutine reorder3 for c more information. c if (choice .eq. 1) then c call reorder3 (cntsome,dndk) c do 600 i=1,cntsome do 610 j=1,shrtcnt if (savglon(j,1) .eq. passavg(i,1)) then passavg(i,2)=savglon(j,2) go to 620 endif 610 continue 620 continue ra(i)=passavg(i,2) 600 continue c call sort (cntsome,4) c write (22,*) ' ' write (22,*) 'new reordering as follows' write (22,*) ' ' do 640 i=1,cntsome do 650 j=1,cntsome if (ra(i) .eq. passavg(j,2)) then passno(i)=int(passavg(j,1)) do 655 k=1,cntsome if (passno(i) .eq. passmjd(k,1)) then i22mjd=passmjd(k,2) x22elev=aver(k,2) go to 656 endif 655 continue 656 continue write (22,*) int(passavg(j,1)),i22mjd, > passavg(j,2),x22elev go to 660 endif 650 continue 660 continue 640 continue write (22,*) ' ' write (22,*) shrtcnt,' short passes as follows' write (22,*) ' ' do 680 i=1,shrtcnt write (22,*) (shrtpas(i,j),j=1,2) 680 continue c call reorder2 (cntsome,10,20) c endif c 999 continue write (*,*) 'total passes read =',countall write (*,*) 'total passes written =',cntsome write (*,*) 'total passes removed =',passrem write (*,*) 'total passes considered to be short =',shrtcnt write (*,*) 'total records read for original input file=',inrec-1 close (10) close (20) close (21) close (22) close (25) stop end c c c subroutine meridian (innum,passnum) real data(400,27) integer idata(400,2),innum,passnum common data(400,27),idata(400,2) c do 100 i=1,innum-1 if (data(i,2) .lt. data(i+1,2)) then write (*,*) passnum,' CROSSES -180.0',data(i,2),data(i+1,2) do 150 ii=1,innum if (data(ii,2) .lt. 0.0) data(ii,2)=data(ii,2)+360.0 150 continue go to 200 endif 100 continue c 200 continue return end c c c subroutine reorder1 (nobs,num) real data(400,27) double precision nobss,along(4000),radd(4000),aalong,aavg,selev, > elev(4000),savg,aver(4000,2),cross integer nobs,num,idata(400,2),passmjd(4000,2),mjd, > passnum common /order1/ aver(4000,2),passmjd(4000,2),cross common data(400,27),idata(400,2) c c----------------------------- subroutine description c c reorder1 takes a given set of longitudes and elevations and c finds the average longitude and elevation for the set. since c some longitudes cross FROM -180.0 TO +180.0 (that is, longitudes c always decrease unless crossing 180) it necessary to correct c the average to the more usual 360 method. therefore the c dataset is ordered from westernmost longitude to eastern most. c NOTE: real*8 is necessary since on some rare occasions the c averages can be the same at real*4 precision. c NOTE: when the study area includes the -180.0 180.0 longitude c line (but does not include all other longitudes) it is c necessary to add 360.0 to the negative (or eastern) c longitudes so that eastern longitudes will be located after c the western longitudes. c NOTE: for datasets that are global (ie. polar datasets or the c whole blasted world) then variable 'west' should be input c as -180.0 and variable 'east' should be input as 180.0. c input as such will produce a map centered on 0.0 c longitude. c nobss=dble(nobs) passnum=idata(1,1) mjd=idata(1,2) do 50 n=1,nobs along(n)=dble(data(n,2)) radd(n)=dble(data(n,3)) 50 continue c aalong=0.0 do 110 n=1,nobs-1 if (along(n) .lt. along(n+1)) then write (*,*) passnum,' CROSSES -180.0 to 180.0', > along(n),along(n+1) go to 130 endif 110 continue do 120 n=1,nobs aalong=aalong+along(n) 120 continue aavg=aalong/nobss go to 150 c 130 aalong=0.0 do 135 n=1,nobs if (along(n) .lt. 0.0) then along(n)=along(n)+360.0 endif aalong=aalong+along(n) 135 continue aavg=aalong/nobss c 150 continue selev=0.0 do 170 n=1,nobs elev(n)=radd(n)-6378.140 selev=selev+elev(n) 170 continue savg=selev/nobss c if (aavg .gt. 180.0) aavg=aavg-360.0 passmjd(num,1)=passnum passmjd(num,2)=mjd aver(num,1)=aavg + cross aver(num,2)=savg c 200 continue return end c c c subroutine reorder2 (nlines,inf,onf) integer nrecord(4000),passno(4000),npoints(4000), > idata(2),pntcnt(4000,2),passrec(4000,2), > shrtpas(3000,2),shrtcnt,inrec2,inf,onf, > pntall(4000,2) real data(27) common /order2/ passno(4000),pntcnt(4000,2),pntall(4000,2) common /shorty/ shrtpas(3000,2),shrtcnt,passrec(4000,2) c c------------------------ the basis for this subroutine was provided c quite generously c by: Dr. D.R.H. O'Connell c Dept. of Geomagicians c Ohio State University c c c------------------------ determine which pass this point belongs to. c if (inf .eq. 10) then do 20 i=1,nlines passrec(i,1)=passno(i) do 30 ii=1,nlines if (passno(i) .eq. pntall(ii,1)) then npoints(i)=pntall(ii,2) go to 35 endif 30 continue 35 continue 20 continue elseif (inf .eq. 21) then do 50 i=1,nlines passrec(i,1)=passno(i) do 40 ii=1,nlines if (passno(i) .eq. pntcnt(ii,1)) then npoints(i)=pntcnt(ii,2) go to 45 endif 40 continue 45 continue 50 continue endif c c------------------------ npoints = number of records to allocate c------------------------ for each pass number. nrecord = the output file c------------------------ record positions for each pass nrecord(1)=1 passrec(1,2)=1 do 60 i=2,nlines i1=i-1 nrecord(i)=nrecord(i1)+npoints(i1) passrec(i,2)=nrecord(i) 60 continue c------------------------ read each data point c c rewind (inf) inrec2=0 70 inrec2=inrec2+1 read (inf,rec=inrec2,err=90) (idata(j),j=1,2),(data(k),k=1,27) c------------------------ determine the matching pass number and its c------------------------ output record number by searching all pass numbers do 80 i=1,nlines if (idata(1) .eq. passno(i)) then write (onf,rec=nrecord(i)) * (idata(j),j=1,2),(data(k),k=1,27) c------------------------ increment output record number for this pass number nrecord(i)=nrecord(i)+1 c------------------------ read next data point goto 70 endif 80 continue go to 70 90 continue c write (*,*) inrec2 - 1,' TOTAL RECORDS READ FOR FILE',inf write (*,*) nrecord(nlines)-1,' TOTAL RECORDS WRITTEN', > ' FOR FILE',onf c------------------------ webe jammin c return end c c c c SUBROUTINE SORT(N,choice) double precision ra(4000),aver(4000,2),rra,cross integer passmjd(4000,2),n,choice common /order1/ aver(4000,2),passmjd(4000,2),cross common /hsort/ ra(4000) c------------------------------this subroutine is written by the authors c of: Numerical Recipes (fortran); c The Art of Scientific Computing c Cambridge University Press c 1989, p. 230 c the routine is referred to as "heapsort" c if (choice .le. 2) then do 10 i=1,n ra(i)=aver(i,choice) 10 continue elseif (choice .eq. 3) then do 30 i=1,n ra(i)=dble(passmjd(i,1)) 30 continue endif c L=N/2+1 IR=N 100 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 ENDIF ENDIF I=L J=L+L 200 IF(J.LE.IR)THEN IF(J.LT.IR)THEN IF(RA(J).LT.RA(J+1))J=J+1 ENDIF IF(RRA.LT.RA(J))THEN RA(I)=RA(J) I=J J=J+J ELSE J=IR+1 ENDIF GO TO 200 ENDIF RA(I)=RRA GO TO 100 END c c c subroutine shorts (innum,hilat,lolat,dndk) integer innum,idata(400,2),dndk,shrtcnt, > shrtpas(3000,2),passrec(4000,2) real data(400,27),hilat,lolat common data(400,27),idata(400,2) common /shorty/ shrtpas(3000,2),shrtcnt,passrec(4000,2) c c------------------------------ subroutine description c shorts determines if the pass is short. a short pass is c a pass which does not extend above the northern-most c (hilat) or below the southern-most (lolat). c if (dndk .eq. 0) then if (data(1,1) .gt. lolat .or. > data(innum,1) .lt. hilat) then shrtcnt=shrtcnt+1 shrtpas(shrtcnt,1)=idata(1,1) shrtpas(shrtcnt,2)=innum endif elseif (dndk .eq. 1) then if (data(innum,1) .gt. lolat .or. > data(1,1) .lt. hilat) then shrtcnt=shrtcnt+1 shrtpas(shrtcnt,1)=idata(1,1) shrtpas(shrtcnt,2)=innum endif endif c return end c c c subroutine reorder3 (allcnt,dndk) integer shrtpas(3000,2),shrtcnt,passrec(4000,2), > allcnt,passnum,recnum,row,recnt,rrow,frow, > isdata(400,2),frec,nrec,rrec,pass,ifdata(400,2), > irdata(400,2),numcnt,minrrow,stocount,minfrow, > dndk,fcnt,rcnt real sdata(400,2),fdata(400,2),rdata(400,2),sftdat(400,2), > ftdata(400,2),rtdata(400,2),totfdiff,totrdiff, > fdiff,rdiff,totsft,totsrt,srtdat(400,2) double precision avgslon,avgrdiff,avgfdiff,savglon(4000,2),crosss common /shorty/ shrtpas(3000,2),shrtcnt,passrec(4000,2) common /trunc/ sdata(400,2),fdata(400,2),rdata(400,2), > sftdat(400,2),ftdata(400,2),srtdat(400,2), > rtdata(400,2) common /order3/ savglon(4000,2),crosss c c---------------------------------------------subroutine description c this subroutine extends all short passes by adding or subtracting c the average longitude difference from the closest full length c pass, then calculates the true average longitude for the short c pass, and stores that true average in an array. the short passes c are found by subroutine shorts with a user defined percentage of total c length to be considered as a long pass. as this percentage is c increased the distance to the closest pass also increases such that c there is more chance for error due to a poor fit of the average c longitude. the closest full length pass is either west or east c of the short pass. if west, then the average difference is added c to the full length pass, however, if east then the average c difference is subtracted. these values are then added to the c longitudes of the short pass to create a set of full length c longitudes which are averaged together to get the true average c longitude. the fundamental principle involved here is that passes c are always parallel and do not actually cross over each other. c therfore, the difference between adjacent passes remains almost c -or at least pretty dog-gone close to almost- constant. c c icnt=1 50 continue do 100 i=1,allcnt if (shrtpas(icnt,1) .eq. passrec(i,1)) then recnum=passrec(i,2) passnum=passrec(i,1) row=shrtpas(icnt,2) recnt=i go to 200 endif 100 continue c 200 continue c---------------------------- read in the short pass do 220 i=1,row read (20,rec=recnum) (isdata(i,j),j=1,2),(sdata(i,j),j=1,2) recnum=recnum+1 220 continue recnum=recnum-1 if (passnum .ne. isdata(row,1)) then write (*,*) 'wrong pass number in reorder3' stop endif c--------------------------- add 360.0 to longitudes that cross the c -180.0 to 180.0 meridian so that averages c are correct do 230 n=1,row-1 if (sdata(n,2) .lt. sdata(n+1,2)) then do 235 i=1,row if (sdata(i,2) .lt. 0.0) sdata(i,2)=sdata(i,2)+360.0 235 continue go to 238 endif 230 continue 238 continue c c--------------------------- find the starting record number for the c nearest east pass frec=0 numcnt=recnt+1 if (numcnt .gt. allcnt) go to 270 nrec=passrec(numcnt,2) 240 continue read (20,rec=nrec) pass do 260 i=1,shrtcnt if (pass .eq. shrtpas(i,1)) then numcnt=numcnt+1 if (numcnt .gt. allcnt) go to 270 nrec=passrec(numcnt,2) go to 240 endif 260 continue frec=nrec c c-------------------------------- find the starting record number for the c nearest west pass 270 continue rrec=0 numcnt=recnt-1 if (numcnt .lt. 1) go to 300 nrec=passrec(numcnt,2) 280 continue read (20,rec=nrec) pass do 290 i=1,shrtcnt if (pass .eq. shrtpas(i,1)) then numcnt=numcnt-1 if (numcnt .lt. 1) go to 300 nrec=passrec(numcnt,2) go to 280 endif 290 continue rrec=nrec c c------------------------------------ calculate the average longitude 300 continue avgfdiff=10000.0 if (frec .gt. 0) then i=1 read (20,rec=frec,err=340) (ifdata(i,j),j=1,2), > (fdata(i,j),j=1,2) i=i+1 frec=frec+1 320 read (20,rec=frec,err=340) (ifdata(i,j),j=1,2), > (fdata(i,j),j=1,2) if (ifdata(i-1,1) .ne. ifdata(i,1)) go to 340 frec=frec+1 i=i+1 go to 320 340 continue frow=i-1 do 344 n=1,frow-1 if (fdata(n,2) .lt. fdata(n+1,2)) then do 342 i=1,frow if (fdata(i,2) .lt. 0.0) fdata(i,2)=fdata(i,2)+360.0 342 continue go to 345 endif 344 continue 345 continue c------------------------------ truncate the short and forward passes c to the same length c call truncate (row,frow,dndk,minfrow,stocount,passnum, > ifdata(frow,1),1) totfdiff=0.0 totsft=0.0 do 350 i=1,minfrow fdiff=abs(sftdat(i,2)-ftdata(i,2)) totfdiff=fdiff+totfdiff totsft=totsft+sftdat(i,2) 350 continue c----------------------------- calculate the average longitude difference avgfdiff=dble(totfdiff/real(minfrow)) endif c c----------------------------- repeat the process for the closest west pass avgrdiff=10000.0 if (rrec .gt. 0) then i=1 read (20,rec=rrec) (irdata(i,j),j=1,2),(rdata(i,j),j=1,2) i=i+1 rrec=rrec+1 360 read (20,rec=rrec) (irdata(i,j),j=1,2),(rdata(i,j),j=1,2) if (irdata(i-1,1) .ne. irdata(i,1)) go to 380 rrec=rrec+1 i=i+1 go to 360 380 continue rrow=i-1 do 384 n=1,rrow-1 if (rdata(n,2) .lt. rdata(n+1,2)) then do 382 i=1,rrow if (rdata(i,2) .lt. 0.0) rdata(i,2)=rdata(i,2)+360.0 382 continue go to 385 endif 384 continue 385 continue call truncate (row,rrow,dndk,minrrow,stocount,passnum, > irdata(rrow,1),-1) totrdiff=0.0 totsrt=0.0 do 390 i=1,minrrow rdiff=abs(srtdat(i,2)-rtdata(i,2)) totrdiff=rdiff+totrdiff totsrt=totsrt+srtdat(i,2) 390 continue avgrdiff=dble(totrdiff/real(minrrow)) endif c c------------------------------------ if the east difference is the smallest c then use the east pass to calculate c average longitude of the extended short c pass if (avgfdiff .lt. avgrdiff) then do 400 i=1,frow if (ftdata(1,2) .eq. fdata(i,2)) then fcnt=i go to 410 endif 400 continue 410 continue if (sftdat(1,2) .lt. ftdata(1,2)) then if (fcnt-1 .eq. 0) go to 425 do 420 i=1,fcnt-1 totsft=totsft+(fdata(i,2)-real(avgfdiff)) 420 continue 425 continue if (minfrow+fcnt .eq. frow) go to 435 do 430 i=minfrow+fcnt,frow totsft=totsft+(fdata(i,2)-real(avgfdiff)) 430 continue 435 continue elseif (sftdat(1,2) .gt. ftdata(1,2)) then if (fcnt-1 .eq. 0) go to 445 do 440 i=1,fcnt-1 totsft=totsft+(fdata(i,2)+real(avgfdiff)) 440 continue 445 continue if (minfrow+fcnt .eq. frow) go to 455 do 450 i=minfrow+fcnt,frow totsft=totsft+(fdata(i,2)+real(avgfdiff)) 450 continue 455 continue endif avgslon=dble(totsft)/dble(frow) c c------------------------------------ repeat the process if the west pass is c the closest elseif (avgrdiff .lt. avgfdiff) then do 500 i=1,rrow if (rtdata(1,2) .eq. rdata(i,2)) then rcnt=i go to 510 endif 500 continue 510 continue if (srtdat(1,2) .lt. rtdata(1,2)) then if (rcnt-1 .eq. 0) go to 525 do 520 i=1,rcnt-1 totsrt=totsrt+(rdata(i,2)-real(avgrdiff)) 520 continue 525 continue if (minrrow+rcnt .eq. rrow) go to 535 do 530 i=minrrow+rcnt,rrow totsrt=totsrt+(rdata(i,2)-real(avgrdiff)) 530 continue 535 continue elseif (srtdat(1,2) .gt. rtdata(1,2)) then if (rcnt-1 .eq. 0) go to 545 do 540 i=1,rcnt-1 totsrt=totsrt+(rdata(i,2)+real(avgrdiff)) 540 continue 545 continue if (minrrow+rcnt .eq. rrow) go to 555 do 550 i=minrrow+rcnt,rrow totsrt=totsrt+(rdata(i,2)+real(avgrdiff)) 550 continue 555 continue endif avgslon=dble(totsrt)/dble(rrow) endif c c------------------------store the average in array if (avgslon .gt. 180.0) avgslon = avgslon - 360.0 savglon(icnt,1)=dble(passnum) savglon(icnt,2)=avgslon + crosss c icnt=icnt+1 if (icnt .gt. shrtcnt) return go to 50 c end c c c subroutine truncate (xrow,yrow,dndk,minrow,stocount, > xpassno,ypassno,fr) integer xrow,yrow,stocount,rowii,rowinc,minrow,fr, > dndk,xpassno,ypassno real xdata(400,2),ydata(400,2),sftdat(400,2),ftdata(400,2), > adata,bdata,diffab,abss, > rtdata(400,2),srtdat(400,2) common /trunc/ sdata(400,2),fdata(400,2),rdata(400,2), > sftdat(400,2),ftdata(400,2),srtdat(400,2), > rtdata(400,2) c c------------------------------ subroutine description c truncate compares the input passes and truncates both c passes to the same overlapping length. c do 20 j=1,xrow do 20 jj=1,2 xdata(j,jj)=sdata(j,jj) 20 continue if (fr) 50,50,60 50 do 55 j=1,yrow do 55 jj=1,2 ydata(j,jj)=rdata(j,jj) 55 continue go to 80 60 do 65 j=1,yrow do 65 jj=1,2 ydata(j,jj)=fdata(j,jj) 65 continue c 80 continue stocount=0 jj=1 rowii=xrow rowinc=yrow c------------------------loops from 90 to 200 increment through the c two input passes and truncate the lengths c to the same length 90 continue adata=xdata(jj,1) bdata=ydata(jj,1) diffab=adata-bdata abss=abs(diffab) if (rowii .eq. 0 .or. rowinc .eq. 0) then c------------------- if this happens, then the findgap subroutine from c movetrunc will have to be implemented to remove c the appropriate pass. so far, there hasn't been c any problems. another alternative would be to c use only the east or west pass instead of c comparing both to the short pass. c write (*,*) 'xrows (s) =',rowii,' yrows (',fr,') =',rowinc write (*,*) 'big problem with pass number =',xpassno,ypassno write (*,*) 'xrow =',xrow,' yrow =',yrow stop endif minrow=min(rowii,rowinc) c------------------------if pass a (ii) matches pass b (inc) at c beginning length then write to xdata and c ydata and race to main program if (abss .lt. 0.33) then if (fr .eq. -1) then do 110 ll=1,minrow do 110 kk=1,2 srtdat(ll,kk)=xdata(ll,kk) rtdata(ll,kk)=ydata(ll,kk) 110 continue elseif (fr .eq. 1) then do 115 ll=1,minrow do 115 kk=1,2 sftdat(ll,kk)=xdata(ll,kk) ftdata(ll,kk)=ydata(ll,kk) 115 continue endif return endif c------------------------if pass a no matcha the b data then find new c a or b depending on whether or not ascending c or descending order of independent variable if (abss .ge. 0.33) then stocount=stocount+1 c------------------------if this is a dusk pass then will count from c -90.0 lat degrees toward the equator if (dndk .eq. 0) then if (xdata(jj,1) .gt. ydata(jj,1)) then rowinc=rowinc-1 do 130 mm=1,rowinc do 130 kk=1,2 ydata(mm,kk)=ydata(mm+1,kk) 130 continue elseif (xdata(jj,1) .lt. ydata(jj,1)) then rowii=rowii-1 do 150 nn=1,rowii do 150 kk=1,2 xdata(nn,kk)=xdata(nn+1,kk) 150 continue endif c------------------------if this is a dawn pass then will count from c the equator toward the south pole c that is decreasing independent variable elseif (dndk .eq. 1) then if (xdata(jj,1) .lt. ydata(jj,1)) then rowinc=rowinc-1 do 160 mm=1,rowinc do 160 kk=1,2 ydata(mm,kk)=ydata(mm+1,kk) 160 continue elseif (xdata(jj,1) .gt. ydata(jj,1)) then rowii=rowii-1 do 170 nn=1,rowii do 170 kk=1,2 xdata(nn,kk)=xdata(nn+1,kk) 170 continue endif endif endif c go to 90 c end c c c c subroutine interp1 (dndk,num,ii) real data(400,27),xdata(27),intpdata(400,27) integer num,idata(400,2),dndk common data(400,27),idata(400,2) common /inta/ xdata(27) common /intb/ intpdata(400,27) c c-------------------------------------- subroutine description c c the basic concept of this subroutine was provided by: c Dr. D.N. "Tiku" Ravat c Dept. of Geology c Purdue University c this subroutine linearly interpolates ALL 27-r variables by c basing the interpolation on the latitudes which are interpolated c at every 0.33 degrees of starting latitude. c ii=0 xlat=real(int(data(1,1)*100.0))/100.0 if (dndk .eq. 0) then if (xlat .lt. data(1,1)) xlat=xlat + 0.33 i=1 100 if (xlat.ge.data(i,1) .and. xlat.le.data(i+1,1)) then call interp2 (i,xlat) xdata(2)=real(int(xdata(2)*100.0))/100.0 ii=ii+1 do 150 j=1,27 intpdata(ii,j)=xdata(j) 150 continue xlat=xlat + 0.33 if (xlat .gt. data(num,1)) return go to 100 elseif (xlat .gt. data(i+1,1)) then i=i+1 go to 100 endif c elseif (dndk .eq. 1) then if (xlat .gt. data(1,1)) xlat=xlat - 0.33 i=1 180 if (xlat.le.data(i,1) .and. xlat.ge.data(i+1,1)) then call interp2 (i,xlat) xdata(2)=real(int(xdata(2)*100.0))/100.0 ii=ii+1 do 200 j=1,27 intpdata(ii,j)=xdata(j) 200 continue xlat=xlat - 0.33 if (xlat .lt. data(num,1)) return go to 180 elseif (xlat .lt. data(i+1,1)) then i=i+1 go to 180 endif endif c end c c c subroutine interp2 (inum,xlat) real data(400,27),diffdata(27),xdata(27) integer inum,idata(400,2) common data(400,27),idata(400,2) common /inta/ xdata(27) c c----------------------------- this subroutine is also from Tiku and is c the interpolator (not to confused with the c terminator!) c do 100 i=1,27 diffdata(i)=data(inum,i)-data(inum+1,i) 100 continue do 120 i=1,27 xdata(i)=data(inum,i)+(xlat-data(inum,1))* > (diffdata(i)/diffdata(1)) 120 continue c return end c c c subroutine despike (npts,outnum) real data(400,27),desdata(400,27),upper,lower integer ic(400),outnum,var,idata(400,2) common data(400,27),idata(400,2) common /spike/ desdata(400,27),upper,lower,var c c--------------------------- this subroutine also provided by Tiku c c c PROGRAM DESPIKE C******************************************************************* C PROGRAM DESPIKE REMOVES MOST SPIKES FROM THE INPUT DATA SET. C HOWEVER, FOR BEST RESULTS, IT IS SUGGESTED TO RUN DESPIKE C AT LEAST THREE TIMES --- FOR EXAMPLE: C C INPUT1 ===DESPIKE===> OUTPUT1 C (OUTPUT1 = INPUT2) ===DESPIKE===> OUTPUT2 C (OUTPUT2 = INPUT3) ===DESPIKE===> OUTPUT3. C C STILL, AFTER RUNNING DESPIKE THREE TIMES, IT FAILS TO ELIMINATE C ORBITS WITH DISCONTINUOUS RESID VS LATITUDE PROFILES. C PROGRAM DEGAP ATTEMPTS TO TAKE CARE OF SUCH PASSES. C C C PARAMETERS TO CHECK: "UPPER" AND "LOWER" (IN NANOTESLAS): C C IF PROGRAM DESPIKE HAS DETERMINED OBSERVATION N TO BE C A GOOD POINT, IT THEN SETS OUT TO DETERMINE IF POINT N+1 C IS A GOOD POINT. IT DOES THIS BY CHECKING THE POINTS C N, N+1, AND N+2. OBSERVATION N+1 WILL BE A GOOD POINT C IF THE RESIDUAL DIFFERENCE BETWEEN POINT N+1 AND THE C POINT ABOVE IT (N OR N+2) IS LESS THAN "UPPER" AND C IF THE RESIDUAL DIFFERENCE BETWEEN IT AND THE POINT C BELOW IT (N+2 OR N) IS GREATER THAN "LOWER". C****************************************************************** C DO 2 I=1,400 IC(I)=1 2 CONTINUE C C****************************************************************** C ARE THE FIRST NEW POINTS SPIKES? C NOTE: DATA(U,23) = RESID2(U) C I=1 15 SL1=(DATA(I+1,var)-DATA(I,var)) SL2=(DATA(I+2,var)-DATA(I,var)) SL3=(DATA(I+3,var)-DATA(I,var)) SL4=(DATA(I+4,var)-DATA(I,var)) XSL=ABS(SL1+SL2+SL3+SL4)/4.0 C S2=(DATA(I+2,var)-DATA(I+1,var)) C IF(ABS(SL1).GT.(3.0*XSL).OR.ABS(SL1).GT.(ABS(3.0*S2))) IC(I)=0 IF(IC(I).EQ.0) THEN I=I+1 GO TO 15 ENDIF C C****************************************************************** C ARE THE MID POINTS SPIKES? C DO 20 J=I,NPTS-2 SL2=(DATA(J+1,var)-DATA(J+2,var)) IF(SL1.GT.UPPER.AND.SL2.LT.LOWER) IC(J+1)=0 IF(SL1.LT.LOWER.AND.SL2.GT.UPPER) IC(J+1)=0 SL1=SL2 20 CONTINUE C C****************************************************************** C IS THE LAST POINT A SPIKE? C K=NPTS-2 25 IF(IC(K).EQ.0) THEN K=K-1 GO TO 25 ENDIF C SL1=ABS(DATA(K,var)-DATA(NPTS-1,var)) SL2=ABS(DATA(NPTS-1,var)-DATA(NPTS,var)) SL3=ABS(DATA(K,var)-DATA(NPTS,var)) IF(IC(NPTS-1).EQ.0) THEN IF(SL1.GT.(3.0*SL2)) IC(NPTS)=0 IF(SL3.GT.(3.0*UPPER)) IC(NPTS)=0 ENDIF IF(IC(NPTS-1).EQ.1) THEN IF(SL2.GT.(3.0*SL1)) IC(NPTS)=0 ENDIF C C NOBS=0 C DO 30 I=1,NPTS C IF(IC(I).EQ.1) NOBS=NOBS+1 C 30 CONTINUE C WRITE(6,*) IDATA(1,1), NOBS C outnum=0 DO 35 I=1,NPTS IF(IC(I).EQ.1) THEN outnum=outnum+1 do 32 m=1,27 desdata(outnum,m)=data(i,m) 32 continue ENDIF 35 CONTINUE C return END