This is the mail archive of the
gcc-bugs@gcc.gnu.org
mailing list for the GCC project.
[Bug tree-optimization/37310] gfortran errors in compilation and the making for upgraded compilers
- From: "petermorgan at grapevine dot net dot au" <gcc-bugzilla at gcc dot gnu dot org>
- To: gcc-bugs at gcc dot gnu dot org
- Date: 2 Sep 2008 12:11:23 -0000
- Subject: [Bug tree-optimization/37310] gfortran errors in compilation and the making for upgraded compilers
- References: <bug-37310-16653@http.gcc.gnu.org/bugzilla/>
- Reply-to: gcc-bugzilla at gcc dot gnu dot org
------- 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.
[peterm@currawong 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