------- Comment #10 from petermorgan at grapevine dot net dot au 2008-09-02 12:11 ------- Subject: Re: gfortran errors in compilation and the making for upgraded compilers
Dear Guys Here are the requests that you asked for. The first run is without any mention of f2c in option names. It still fails Running unimake to create Makefile for blsum System name: Linux currawong 2.6.25.14-108.fc9.x86_64 #1 SMP Mon Aug 4 13:46:35 EDT 2008 x86_64 x86_64 x86_64 GNU/Linux System release number translated to 2625 No i86 compiler specification--assuming gfortran (gcc 4.2x) Makefile for blsum remade by unimake gfortran -O -Wuninitialized -fno-automatic -fno-backslash blsum.f ../gen_util/gen_util_lib.a ../../libraries/matrix/kinv_lib.a ../../libraries/comlib/com_lib.a -o blsum rm -f blsum.o gfortran -O -Wuninitialized -fno-automatic -fno-backslash ensum.f ../gen_util/gen_util_lib.a ../../libraries/matrix/kinv_lib.a ../../libraries/comlib/com_lib.a -o ensum rm -f ensum.o gfortran -O -Wuninitialized -fno-automatic -fno-backslash tssum.f ../gen_util/gen_util_lib.a ../../libraries/matrix/kinv_lib.a ../../libraries/comlib/com_lib.a -o tssum tssum.f: In function 'remove_ejmp': tssum.f:712: internal compiler error: in set_uids_in_ptset, at tree-ssa-structalias.c:4789 Please submit a full bug report, with preprocessed source if appropriate. See <http://gcc.gnu.org/bugs.html> for instructions. make: *** [tssum] Error 1 Failure in make_globk -- install_software terminated ----------------------------------------------------------------------- Here is a second run with only -O set It also has the same error in blsum Running unimake to create Makefile for blsum System name: Linux currawong 2.6.25.14-108.fc9.x86_64 #1 SMP Mon Aug 4 13:46:35 EDT 2008 x86_64 x86_64 x86_64 GNU/Linux System release number translated to 2625 No i86 compiler specification--assuming gfortran (gcc 4.2x) Makefile for blsum remade by unimake gfortran -O blsum.f ../gen_util/gen_util_lib.a ../../libraries/matrix/kinv_lib.a ../../libraries/comlib/com_lib.a -o blsum rm -f blsum.o gfortran -O ensum.f ../gen_util/gen_util_lib.a ../../libraries/matrix/kinv_lib.a ../../libraries/comlib/com_lib.a -o ensum rm -f ensum.o gfortran -O tssum.f ../gen_util/gen_util_lib.a ../../libraries/matrix/kinv_lib.a ../../libraries/comlib/com_lib.a -o tssum tssum.f: In function 'remove_ejmp': tssum.f:712: internal compiler error: in set_uids_in_ptset, at tree-ssa-structalias.c:4789 Please submit a full bug report, with preprocessed source if appropriate. See <http://gcc.gnu.org/bugs.html> for instructions. make: *** [tssum] Error 1 Failure in make_globk -- install_software terminated ------------------------------------------------- Now I know that I had version 10.33 running under 4.2.2 on a 32 bit machine so I built a version 4.2.2 language system on my 64 bit machine. [EMAIL PROTECTED] blsum]$ gcc -v Using built-in specs. Target: x86_64-unknown-linux-gnu Configured with: ../gcc-4.2.2.source/configure --prefix=/opt --with-mpfr-lib=/opt/lib --with-mpfr-include=/opt/include --disable-multilib Thread model: posix gcc version 4.2.2 There are no errors in a compilation of GAMIT/GLOBK when I use this package. I have attached the code modules for you. As I said they compile okay under 4.2.2 Please let me know how this helps. cheers Peter kargl at gcc dot gnu dot org wrote: > ------- Comment #9 from kargl at gcc dot gnu dot org 2008-09-02 03:46 ------- > (In reply to comment #8) > >> gfortran -O3 -Wuninitialized -fno-f2c -ffast-math -fno-automatic >> -fno-backslash tssum.f ../gen_util/gen_util_lib.a >> ../../libraries/matrix/kinv_lib.a ../../libraries/comlib/com_lib.a -o tssum >> tssum.f: In function 'remove_ejmp': >> tssum.f:712: internal compiler error: in set_uids_in_ptset, at >> tree-ssa-structalias.c:4789 >> Please submit a full bug report, >> > > What happens if you use a sane set of options? In particular, > DO NOT USE any option with f2c in it name. DO NOT USE --ffast-math, > which is misnamed. What happens if you use -O or -O2. > > Finally, is the source available? > > > -- *************************************** * * * Peter and Carol Morgan * * 20 Goodparla St * * Hawker, ACT, 2614 * * Australia * * * * Home Phone +61 (0)2 6254 0137 * * Peter's Mobile +61 (0)4 1854 0137 * * * * * *************************************** program tssum * Program to generate PBO time series files * * Usages: * tssum <dir> <prod_id> <-R> <list of .org files> * * where <dir> -- directory to put the time series in. This directory * will be checked for existing files and these will be appended to * <prod_id> is product id with form: pbo.final_frame. Characters * 5:9 will used for time series type (normally rapid or final) * <-R> causes exisiting time series files to be ignored * <list of .org files> glred/globk org files with output option PBOP * output option set. * * PROD_ID types valid in PBO: * Format <cen>.<series>_<frame>+<type> * <cen> - Center 3 characters (bsl cwu pbo mit) or special product * code e.g., aug for Augustine Volcano, aku for Akutan volcano * <series> - Orbit series: final or rapid * <frame> - Frame type: loose or frame * <type> - Optional type. If not given series name is used. Additional * type put in the final series is supplemental run (suppl). * * Stsndard PROD_ID * pbo.rapid_frame * pbo.final_frame * pbo.final_frame+suppl include 'tssum.h' integer*4 len_run, nr, ierr, ns, rcpar, trimlen character*64 runstring **** OK, Read the runsting (for the moment just by position) len_run = rcpar(1,tsdir) ln_tsdir = len_run if( len_run.eq.0 ) then call proper_runstring('tssum.hlp','tssum',1) endif len_run = rcpar(2,prod_id) * MOD TAH 051129: Extract type from prod_id: Format pbo.final_frame+suppl ts_ref_type = prod_id(5:9) if( len_run.gt.15 ) then ts_ref_type = prod_id(17:21) prod_id = prod_id(1:15) endif ln_prod_id = trimlen(prod_id) len_run = rcpar(3,runstring) if( runstring(1:2).eq.'-R' ) then ts_refresh = .true. else ts_refresh = .false. endif ***** OK Loop over input files num_ent = 0 num_site = 0 num_code = 0 nr = 3 len_run = 1 do while ( len_run.gt.0 ) nr = nr + 1 len_run = rcpar(nr,in_file) if( len_run.gt.0 ) then call read_in endif end do call systime( date_rel, sec_rel) **** OK, now loop over all the codes do ns = 1, num_code * Generate TS file name ts_file = tsdir(1:ln_tsdir) // '/' // in_code(ns) // '.' // . prod_id(1:ln_prod_id) // '.pos' ***** Try to open file unless we are refreshing if( .not.ts_refresh ) then open(200,file=ts_file,iostat=ierr,status='old') if( ierr.ne.0 ) then new_ts = .true. else new_ts = .false. end if else new_ts = .true. end if ***** See if we should read old series num_ts = 0 if( .not.new_ts ) then call read_ts(200) else call gen_refdata(ns) endif * MOD TAH 080325: Before merging correct input series for any * East jumps due to latitude shifts call remove_ejmp(ns) ***** OK, merge new entries with old call merge_ts(ns) ***** Now close the infile and re-open to write close(200) open(200,file=ts_file,iostat=ierr,status='unknown') call write_ts(200, ns) end do ***** Thats all end CITTLE READ_IN subroutine read_in include 'tssum.h' character*8 gsite_name character*16 gsite_full integer*4 date(5), ierr, jerr, j, ns, cs, ne integer*4 trimlen, indx real*8 gmjd, pos_xyz_fin(3),xyz_std(6), unc_geod(3),llu_std(3), . pos_neu_fin(3), neu_std(6), sec character*512 line **** Open the input file open(100,file=in_file, iostat=ierr, status='old') call report_error('IOSTAT',ierr,'open',in_file,0,'tssum') sec = 0.d0 do while ( ierr.eq.0 ) read(100,'(a)', iostat=ierr) line * MOD TAH 060921: Replaced XPS with X since eq-renames may change the name if( ierr.eq.0 .and. line(1:4).eq.'pbo.' .and. . line(11:11).ne.'X' ) then read(line,100,iostat=jerr) gsite_name, gsite_full, . date, gmjd, . (pos_xyz_fin(j), j=1,3), (xyz_std(j),j=1,6), . (unc_geod(j),j=1,3), (llu_std(j),j=1,3), . (pos_neu_fin(j),j=1,3), (neu_std(j),j=1,6) 100 format(5x,a8,1x,a16,1x,i5,4(1x,i2.2),1x,F10.4, . 1x,3F15.5,3F8.5,3F7.3,3x, . 2F16.10,1x,F10.5,1x,2F8.1,1x,F10.5,3x, . 2F15.5,1x,F10.5,3F8.5,3F7.3) call report_error('IOSTAT',jerr,'read',line,0,'read_in') c100 format('pbo. ',a8,1x,a16,1x,i5,4(1x,i2.2),1x,F10.4, c . 1x,3F15.5,3F8.5,3F7.3,' | ', c . 2F16.10,1x,F10.5,1x,2F8.1,1x,F10.5,' | ', c . 2F15.5,1x,F10.5,3F8.5,3F7.3) call ymdhms_to_mjd(date, sec, gmjd) if( jerr.eq.0 ) then * See if we can match site name indx = 0 call get_cmd(gsite_name,in_site,num_site,ns,indx) if( ns.le.0 ) then num_site = num_site + 1 ns = num_site in_site(ns) = gsite_name end if indx = 0 call get_cmd(gsite_name(1:4),in_code,num_code,cs,indx) if( cs.le.0 ) then num_code = num_code + 1 cs = num_code in_full(cs) = gsite_full in_code(cs) = gsite_name(1:4) end if ***** OK save this entry num_ent = num_ent + 1 ne = num_ent in_ns(ne) = ns in_cs(ne) = cs in_mjd(ne) = gmjd do j = 1,3 in_xyz(j,ne) = pos_xyz_fin(j) in_neu(j,ne) = pos_neu_fin(j) in_llu(j,ne) = unc_geod(j) end do do j = 1,6 in_xyz_std(j,ne) = xyz_std(j) in_neu_std(j,ne) = neu_std(j) end do end if end if end do **** Tell user were we are write(*,110) in_file(1:trimlen(in_file)), num_site, num_code, . num_ent 110 format('File: ',a,' Sites ',i5,' Codes ',i5,' Entries ',i10) return end subroutine read_ts(unit) include 'tssum.h' integer*4 unit, ierr, jerr, trimlen, date(5), i, j, isec integer*4 ts_srt(max_ts) ! indexs to sorted time series integer*4 kn, cnt, in real*8 sec, rot(3,3) real*8 max_mjd, min_mjd ! Max and min m,jd used in sorting logical line_read ! Set true when line has been read checking for ! Field dscriptions. logical time_order_out ! Set true to denote file out of time order logical set ! Set true when ranked value put in ts_srt character*512 line **** Read the header records read(unit,'(a)',iostat=ierr) line * MOD TAH 051129: See if version line next call report_error('IOSTAT',ierr,'read',line,1,'read_ts') read(unit,'(a)',iostat=ierr) line if( index(line,'Version').gt.0 ) then read(line,105,iostat=ierr) ts_ver_read 105 format(16x,a) c105 format('Format Version: ',a) call report_error('IOSTAT',ierr,'read',line,1,'read_ts') read(unit,'(a)',iostat=ierr) line endif read(line,110,iostat=ierr) ts_code 110 format(16x,a) c110 format('4-character ID: ',a4) call report_error('IOSTAT',ierr,'read',line,1,'read_ts') read(unit,120,iostat=ierr) ts_full 120 format(16x,a16) c120 format('Station name : ',a16) call report_error('IOSTAT',ierr,'read','Station',1,'read_ts') read(unit,'(a)',iostat=ierr) line read(line,130,iostat=ierr) date, isec 130 format(16x,i4,i2,i2,1x,i2,i2,i2) sec = isec c130 format('First Epoch : ', i4,i2.2,i2.2,1x,i2.2,i2.2,i2.2) call report_error('IOSTAT',ierr,'read',line,1,'read_ts') call ymdhms_to_mjd(date,sec,ts_first) cFirst Epoch : 20060521 120000 read(unit,140,iostat=ierr) date, isec sec = isec 140 format(16x, i4,i2,i2,1x,i2,i2,i2) call report_error('IOSTAT',ierr,'read','Last Epoch',1,'read_ts') c140 format('Last Epoch : ', i4,i2.2,i2.2,1x,i2.2,i2.2,i2.2) call ymdhms_to_mjd(date,sec,ts_last) C**** Do not read release date, since we will generate a new one here C read(unit,150,iostat=ierr) date_rel, nint(sec_rel) read(unit,'(a)',iostat=ierr) line 150 format(16x, i4,i2,i2,1x,i2,i2,i2) c150 format('Release Data : ', i4,i2.2,i2.2,1x,i2.2,i2.2,i2.2) call report_error('IOSTAT',ierr,'read',line,1,'tssum') read(unit,200,iostat=ierr) ref_xyz 200 format(25x,3F15.5) c200 format('XYZ Reference position : ',3F15.5) call report_error('IOSTAT',ierr,'read',line,1,'tssum') call XYZ_to_GEOD(rot,ref_xyz, ref_llu) call loc_to_geod(ref_llu, ref_neu) C Use compute version C read(unit,220,iostat=ierr) (pi/2-ref_llu(1))*180/pi, C . ref_llu(2)*180/pi,ref_llu(3) read(unit,'(a)',iostat=ierr) line 220 format(25x,2F16.10,1x,F10.5) call report_error('IOSTAT',ierr,'read',line,1,'tssum') c220 format('NEU Reference position : ',2F16.10,1x,F10.5) * MOD TAH 080108: See if field descrription lines are present read(unit,'(a)',iostat=ierr) line line_read = .true. if( index(line,'Start Field Description').gt.0 ) then do while (index(line,'End Field Description').eq.0 ) read(unit,'(a)',iostat=ierr) line call report_error('IOSTAT',ierr,'read',ts_file,0, . 'tssum/Field Description') enddo line_read = .false. endif **** OK, Now read in the time series i = 0 time_order_out = .false. do while ( ierr.eq.0 ) if( .not.line_read ) then read(unit,'(a)',iostat=ierr) line endif line_read = .false. if( ierr.eq.0 ) then i = i + 1 read(line,300,iostat=jerr) date, isec, ts_mjd(i), . (ts_xyz(j,i),j=1,3), (ts_xyz_std(j,i),j=1,6), . (ts_llu(j,i),j=1,3), . (ts_neu(j,i),j=1,3), (ts_neu_std(j,i),j=1,6), . ts_type(i) sec = isec 300 format(1x,i4,i2,i2,1x,i2,i2,i2,1x,F10.4, . 3F15.5,3F9.5,3F7.3,3x,2F16.10,1x,F10.5,3x, . 3(F9.5,1x),1x,3F9.5,3F7.3,1x,a5) if( jerr.ne.0 ) then write(*,310) jerr, line(1:trimlen(line)) 310 format('JERR ',i4,1x,a) endif call ymdhms_to_mjd(date,sec,ts_mjd(i)) do j = 1,3 ts_neu(j,i) = ts_neu(j,i) + ref_neu(j) end do **** Make sure that time series is in time order. if( i.gt.1 .and. ts_mjd(i).lt.ts_mjd(i-1) ) then write(*,320) ts_file(1:trimlen(ts_file)),i 320 format('**WARNING** Time series in ',a, . ' out of time order at record ',i5) time_order_out = .true. endif end if end do * MOD TAH 080104: If time seris out of order, resort now. num_ts = i if( time_order_out ) then * Buildup up list of indices in time order. In this * approach dupliucate times are replaced with later * values. do i = 1,num_ts ts_srt(i) = 0 end do max_mjd = -1e20 cnt = 0 do i = 1, num_ts min_mjd = 1e20 set = .false. do j = 1, num_ts if( ts_mjd(j) .le. min_mjd .and. . ts_mjd(j) .gt. max_mjd ) then ts_srt(i) = j min_mjd = ts_mjd(j) set = .true. endif enddo if( set ) then cnt = cnt + 1 max_mjd = ts_mjd(ts_srt(i)) end if enddo * ts_srt now is sorted list. Sort the list of values do i = 1, num_ts in = ts_srt(i) * if in > 0 then swap the ith entry with the in'th * value call save_ts(i,0) ! cp i'th entries to save values (0) call save_ts(in,i) ! cp in'th entries to ith values call save_ts(0,in) ! cp saved values to in'th entry * Now update ts_srt to show values have moved kn = 0 do j = i+1,num_ts if( ts_srt(j).eq.i ) then kn =j endif enddo if( kn.gt. 0 ) then ts_srt(kn) = in endif end do * Save number of values in time series num_ts = cnt endif write(*,400) ts_file(1:trimlen(ts_file)), num_ts 400 format('TS_file ',a,' sorted entries ',i5) **** Thats all return end CTITLE SAVE_TS subroutine save_ts(in,out) * Routine to save time series values from index in to index out * If out is zero, save to or from save arrrays include 'tssum.h' integer*4 in, out * LOCAL integer*4 i **** See which way we are saving if( in.gt.0 .and. out.gt.0 ) then ts_mjd(out) = ts_mjd(in) ts_type(out) = ts_type(in) do i = 1,3 ts_xyz(i,out) = ts_xyz(i,in) ts_llu(i,out) = ts_llu(i,in) ts_neu(i,out) = ts_neu(i,in) enddo do i = 1,6 ts_xyz_std(i,out) = ts_xyz_std(i,in) ts_neu_std(i,out) = ts_neu_std(i,in) enddo elseif( out .eq. 0 ) then sv_mjd = ts_mjd(in) sv_type = ts_type(in) do i = 1,3 sv_xyz(i) = ts_xyz(i,in) sv_llu(i) = ts_llu(i,in) sv_neu(i) = ts_neu(i,in) enddo do i = 1,6 sv_xyz_std(i) = ts_xyz_std(i,in) sv_neu_std(i) = ts_neu_std(i,in) enddo elseif( in.eq.0 ) then ts_mjd(out) = sv_mjd ts_type(out) = sv_type do i = 1,3 ts_xyz(i,out) = sv_xyz(i) ts_llu(i,out) = sv_llu(i) ts_neu(i,out) = sv_neu(i) enddo do i = 1,6 ts_xyz_std(i,out) = sv_xyz_std(i) ts_neu_std(i,out) = sv_neu_std(i) enddo else call report_error('PROG',-1,'sav','timeseries',0,'save_ts') endif return end CTITLE GEN_REFDATA subroutine gen_refdata(ns) include 'tssum.h' integer*4 ns ! code site number real*8 av_xyz(3) integer*4 na, i, j **** Routine to generate reference XYZ do j = 1,3 av_xyz(j) = 0.0d0 end do na = 0 do i = 1, num_ent if( in_cs(i).eq.ns) then if( in_xyz_std(1,i).lt.0.1d0 .and. . in_xyz_std(2,i).lt.0.1d0 .and. . in_xyz_std(3,i).lt.0.1d0 ) then na = na + 1 do j = 1,3 av_xyz(j) = av_xyz(j) + in_xyz(j,i) end do end if end if end do if( na.gt.0 ) then do j = 1,3 ref_xyz(j) = av_xyz(j)/na end do else * Nothing with sigmas less than 10 cm: So what to do? * Use all available data? do i = 1, num_ent if( in_cs(i).eq.ns) then na = na + 1 do j = 1,3 av_xyz(j) = av_xyz(j) + in_xyz(j,i) end do end if end do if( na.gt.0 ) then do j = 1,3 ref_xyz(j) = av_xyz(j)/na end do end if endif **** Thats all return end subroutine merge_ts(ns) * Routine to merge the new time series values with the old ones * If a time matches with an old value, the older value will be replaced. include 'tssum.h' integer*4 ns ! code site number integer*4 i,j,k, n, l logical replaced ! Logical set true if entry is being replaced with ! a new one. do i = 1, num_ent if( in_cs(i).eq.ns) then **** The site code matches, scan the old time series values * to see we need to replace the value, insert a new one or * append to the old list (later would be most common). * Times may differ by a minute or so; so test with toleranace C if( num_ts.gt.0 .and. in_mjd(i).le.ts_mjd(num_ts) ) then replaced = .false. if( num_ts.gt.0 .and. . in_mjd(i)-ts_mjd(num_ts).lt.1.d-3 ) then * This entry is earlier. Find out where to put it. n = num_ts do while ( n.ge.1 ) C if ( ts_mjd(n).eq.in_mjd(i) ) then if ( abs(ts_mjd(n)-in_mjd(i)).le.1.d-3 ) then * Matches the time, so replace the old values * with the new k = n n = 0 replaced = .true. elseif ( ts_mjd(n).lt.in_mjd(i) .or.n.eq.1) then * Need to move up all the old values and insert * the new entry in at n+1 k = n+1 * Move up the old values do l = num_ts+1, k, -1 ts_mjd(l) = ts_mjd(l-1) ts_type(l) = ts_type(l-1) do j = 1, 3 ts_xyz(j,l) = ts_xyz(j,l-1) ts_neu(j,l) = ts_neu(j,l-1) ts_llu(j,l) = ts_llu(j,l-1) enddo do j = 1, 6 ts_xyz_std(j,l) = ts_xyz_std(j,l-1) ts_neu_std(j,l) = ts_neu_std(j,l-1) enddo end do num_ts = num_ts + 1 n = 0 else n = n -1 end if end do else * Normal case of just adding new results to end num_ts = num_ts + 1 k = num_ts endif * Now add the new values ts_mjd(k) = in_mjd(i) * Get which code we should add for this point. if( replaced .and. ts_ref_type.eq.'suppl' ) then ts_type(k) = 'suppf' elseif( replaced .and. ts_ref_type.eq.'campd' ) then ts_type(k) = 'campf' else ts_type(k) = ts_ref_type endif do j = 1, 3 ts_xyz(j,k) = in_xyz(j,i) ts_neu(j,k) = in_neu(j,i) ts_llu(j,k) = in_llu(j,i) enddo do j = 1, 6 ts_xyz_std(j,k) = in_xyz_std(j,i) ts_neu_std(j,k) = in_neu_std(j,i) enddo end if end do **** Thats all return end subroutine write_ts(unit,ns) include '../includes/const_param.h' include 'tssum.h' integer*4 ns ! code site number integer*4 unit ! unit number integer*4 ierr, i, j, date(5) real*8 sec, rot(3,3) real*8 ts_mag, Cts_mag ! Magnitude of XYZ/NEU and function to ! compute it. ***** Scan the ts list and get start and stop times ts_first = 1.d6 ts_last = 0.d0 do i = 1, num_ts if( ts_mjd(i).lt.ts_first ) ts_first = ts_mjd(i) if( ts_mjd(i).gt.ts_last ) ts_last = ts_mjd(i) end do ***** Write out the header lines write(unit,100,iostat=ierr) 100 format('PBO Station Position Time Series') * MOD TAH 051129: Write out version number write(unit,105,iostat=ierr) ts_ver 105 format('Format Version: ',a) write(unit,110,iostat=ierr) in_code(ns) 110 format('4-character ID: ',a4) write(unit,120,iostat=ierr) in_full(ns) 120 format('Station name : ',a16) call jd_to_ymdhms(ts_first+1d-6, date,sec) write(unit,130,iostat=ierr) date, nint(sec) 130 format('First Epoch : ', i4,i2.2,i2.2,1x,i2.2,i2.2,i2.2) call jd_to_ymdhms(ts_last+1d-6, date,sec) write(unit,140,iostat=ierr) date, nint(sec) 140 format('Last Epoch : ', i4,i2.2,i2.2,1x,i2.2,i2.2,i2.2) write(unit,150,iostat=ierr) date_rel, nint(sec_rel) 150 format('Release Data : ', i4,i2.2,i2.2,1x,i2.2,i2.2,i2.2) write(unit,200,iostat=ierr) ref_xyz 200 format('XYZ Reference position : ',3F15.5) call XYZ_to_GEOD(rot,ref_xyz, ref_llu) call loc_to_geod(ref_llu, ref_neu) write(unit,220,iostat=ierr) (pi/2-ref_llu(1))*180/pi, . ref_llu(2)*180/pi,ref_llu(3) 220 format('NEU Reference position : ',2F16.10,1x,F10.5) **** OK, Now write out the time series do i = 1, num_ts call jd_to_ymdhms(ts_mjd(i)+1.d-6,date,sec) * * Check size ts_mag = Cts_mag(ts_neu(1,i),ref_neu) if( ts_mag.lt.99.999d0 ) then * Use Standard format write(unit,300) date, nint(sec), ts_mjd(i), . (ts_xyz(j,i),j=1,3), (ts_xyz_std(j,i),j=1,6), . (ts_llu(j,i),j=1,3), . (ts_neu(j,i)-ref_neu(j),j=1,3), . (ts_neu_std(j,i),j=1,6), . ts_type(i) 300 format(1x,i4,i2.2,i2.2,1x,i2.2,i2.2,i2.2,1x,F10.4, . 3F15.5,3F9.5,3F7.3,3x,2F16.10,1x,F10.5,3x, . 3(F9.5,1x),1x,3F9.5,3F7.3,1x,a) else write(unit,310) date, nint(sec), ts_mjd(i), . (ts_xyz(j,i),j=1,3), (ts_xyz_std(j,i),j=1,6), . (ts_llu(j,i),j=1,3), . (ts_neu(j,i)-ref_neu(j),j=1,3), . (ts_neu_std(j,i),j=1,6), . ts_type(i) 310 format(1x,i4,i2.2,i2.2,1x,i2.2,i2.2,i2.2,1x,F10.4, . 3F15.5,3F9.5,3F7.3,3x,2F16.10,1x,F10.5,3x, . 3(F9.3,1x),1x,3F9.5,3F7.3,1x,a) endif end do **** Thats all return end CTITLE CTS_MAG real*8 function CTS_mag(NEU, Ref) implicit none * routine to compute magnitude of max element real*8 NEU(3), Ref(3) real*8 dNEU(3) integer*4 i CTS_mag = 0.d0 do i = 1,3 dneu(i) = NEU(i)-REF(i) if( abs(dneu(i)).gt.CTS_mag ) CTS_mag = abs(dneu(i)) end do return end CTITLE REMOVE_EJMP subroutine remove_ejmp(ns) * Subroutine to remove east jumps that can occur when a site * moves too far in latitude. Usually only a problem for ice * sites. Here we compare the reference lat with the latitude * in the individual entries include '../includes/const_param.h' include 'tssum.h' integer*4 ns ! code site number integer*4 i real*8 locref(3), locdat(3) ! Geod Co-lat, long and heigt ! with reference latitude and ! day value of latitude real*8 rot(3,3) real*8 neuref(3), neudat(3) ! NEU values for reference and ! daily latiude. If no latitude ! induced East jump, then should ! be the same real*8 dEast, old_dEast ! East change and previous east change ! (latter used for output) * ***** Get the reference positions for this site. call XYZ_to_GEOD(rot,ref_xyz, ref_llu) * Loop over entries finding the one for this site old_dEast = 0.0d0 do i = 1, num_ent if( in_cs(i).eq.ns) then * Test to see if change. Call loc_to_geod with * actual and refence latitude locdat(1) = pi/2-in_llu(1,i)*pi/180 locref(2) = in_llu(2,i)*pi/180 locdat(2) = in_llu(2,i)*pi/180 locref(3) = in_llu(3,i) locdat(3) = in_llu(3,i) * Now change reference lat in locref locref(1) = ref_llu(1) * Convert to NEU call loc_to_geod(locref, neuref) call loc_to_geod(locdat, neudat) c print *,' Ref ',i, locref, neuref c print *,' Dat ',i, locdat, neudat * Get the East difference dEast = neuref(2)-neudat(2) if( abs(dEast).gt.1.d-4 ) then in_neu(2,i) = in_neu(2,i)+dEast endif if( abs(dEast-old_dEast).gt.1.d-2) then write(*,220) in_code(ns), in_mjd(i), dEast 220 format('East Offset at ',a4,' MJD ',F12.4, . ' dEast ',F15.3,' m') endif old_dEast = dEast endif end do ***** Thats all return end * DECLARATIONS of tssum integer*4 max_ent ! Maximum number of total entries integer*4 max_site ! Maximum number of sites integer*4 max_ts ! Maximum number of entries per site character*(*) ts_ver ! Version of time series files parameter ( max_ent = 5000000 ) parameter ( max_site = 2048 ) parameter ( max_ts = 10000 ) * * MOD TAH 051129: Added version 1.0.1: Version number in file and type * designation (rapid/final) parameter ( ts_ver = '1.0.1' ) * Declaratons integer*4 num_ent, num_site, num_code, num_ts * Data records from .org files integer*4 in_ns(max_ent), in_cs(max_ent) integer*4 ln_tsdir, ln_prod_id integer*4 date_rel(5) real*8 sec_rel real*8 in_mjd(max_ent), in_xyz(3,max_ent), in_xyz_std(6,max_ent), . in_llu(3,max_ent), . in_neu(3,max_ent),in_neu_std(6,max_ent) real*8 ts_mjd(max_ts), ts_xyz(3,max_ts), ts_xyz_std(6,max_ts), . ts_llu(3,max_ts), . ts_neu(3,max_ts),ts_neu_std(6,max_ts) real*8 sv_mjd, sv_xyz(3), sv_xyz_std(6), . sv_llu(3), . sv_neu(3),sv_neu_std(6) real*8 ts_first, ts_last real*8 ref_xyz(3), ref_llu(3), ref_neu(3) character*8 in_site(max_site) character*4 in_code(max_site), ts_code character*128 tsdir, prod_id character*5 ts_type(max_ts) ! Based on prod_id, rapid/final and any new ! type passed in the prod_id. Normally updated from ! the value read in. ., sv_type ., ts_ref_type ! New reference type based on prod_id character*128 in_file, ts_file character*16 in_full(max_site), ts_full character*16 ts_ver_read logical ts_refresh, new_ts common / ts_i4 / num_ent, num_site, num_code, num_ts, . in_ns, in_cs, date_rel, ln_tsdir, ln_prod_id, . ts_refresh, new_ts common / ts_r8 /sec_rel, in_mjd, in_xyz, in_xyz_std, . in_llu, in_neu,in_neu_std, ts_mjd, ts_xyz, ts_xyz_std, . ts_llu, ts_neu,ts_neu_std, ts_first, ts_last, . ref_xyz, ref_llu, ref_neu, . sv_mjd, sv_xyz, sv_xyz_std, sv_llu, sv_neu,sv_neu_std common / ts_ch /in_site,in_code, ts_code, tsdir, . prod_id, in_file, ts_file, in_full, ts_full, ts_ver_read, . ts_type, sv_type, ts_ref_type *--------------------------------------------------------------------- * CONST_PARAM.FTNI * This include file gives all of the constants which are used * in the Kalman filter processing software. * * T.Herring 12:41 PM TUE., 13 JAN., 1987 *--------------------------------------------------------------------- * earth_flat - Earth's flattening * earth_rad - Equatorial radius of the Earth (m) -- Also the * - semi-major-axis of the ellipsoid. * earth_to_moon - Mass ratio of earth and moom * g_earth - Gravitional acceleration at the equator (m/s**2) * GM_moon - GM for moon * GM_sun - GM for sun * GM_earth - GM for Earth * G_univ - Gravitional constant * pi - PI * rad_to_deg - Conversion from radians to degs. * rad_to_mas - Conversion from radians to milliarcsecs. * sec_per_day - Number of seconds in 24 hours * vel_light - speed of light in m/s * DJ2000 - Julian date of J2000 * sec360 - number of seconds in 360 degreees. * solar_to_sidereal - Conversion from solar days to sidereal * - days (at J2000) * fL1, fL2 - GPS frequencies in Hz at L1 and L2 * dfsf, sfdf - Difference of frequency divided by the sum of * - frequencies (used form widelane and * - narrowlane) * lcf1, lcf2 - Multipliers for LC from L1 and L2 frequencies * lgf1, lgf2 - Multipliers for LG from L1 and L2 frequencies * pcf1, pcf2 - Multipliers for PC from P1 and P2 frequencies real*8 earth_flat, earth_rad, earth_to_moon, g_earth, . GM_earth, GM_moon, GM_sun, G_univ, . pi, rad_to_deg, rad_to_mas, sec_per_day, . vel_light, DJ2000, sec360, solar_to_sidereal, fL1, fL2 real*8 dfsf, sfdf, lcf1, lcf2, lgf1, lgf2, pcf1, pcf2 * Values for the parameters C parameter ( earth_flat = 0.003352891869D0 ) * ! m C parameter ( earth_rad = 6378145.D0 ) * WGS-84 parameters for the elliposoid (a and 1/f) parameter ( earth_rad = 6378137.D0 ) parameter ( earth_flat = 1.d0/298.257222101d0 ) parameter ( earth_to_moon = 81.30065918D0 ) * ! m/sec**2 parameter ( g_earth = 9.780318458D0 ) * ! m**3/sec**2 parameter ( GM_moon = 0.49027975D+13 ) * ! m**3/sec**2 parameter ( GM_sun = 0.132712499D+21 ) * ! m**3/sec**2 * IERS Value 3.986004418d+14; solve value 3.9860346D+14 parameter ( GM_earth = 3.986004418d+14 ) * ! m**2/(kg sec**2) parameter ( G_univ = 0.66732D-10 ) parameter ( pi = 3.1415926535897932D0 ) parameter ( sec_per_day = 86400.D0 ) * ! m/s parameter ( sec360 = 1296000.d0 ) parameter ( vel_light = 299792458.D0 ) * ! Julian days parameter ( DJ2000 = 2451545.d0 ) parameter ( solar_to_sidereal = 1.002737909d0 ) parameter ( fL1 = 154*10.23d6 ) parameter ( fL2 = 120*10.23d6 ) * Computed quanities parameter ( rad_to_deg = 180.d0 /pi ) parameter ( rad_to_mas = 648000.d3/pi ) parameter ( dfsf = (fL1-fL2)/(fL1+fL2) ) parameter ( sfdf = (fL1+fL2)/(fL1-fL2) ) parameter ( lcf1 = 1.d0/(1.d0 - (fL2/fL1)**2) ) parameter ( lcf2 = -(fL2/fL1)/(1.d0 - (fL2/fL1)**2) ) parameter ( lgf1 = -fL2/fL1) parameter ( lgf2 = 1.d0 ) parameter ( pcf1 = fL1**2/(fL1**2-fL2**2) ) parameter ( pcf2 = -fL2**2/(fL1**2-fL2**2) ) -- http://gcc.gnu.org/bugzilla/show_bug.cgi?id=37310
