Line data Source code
1 47 : subroutine rdfreq(icase,ids,cs,l,nord,sig,frq,ekin,ierr)
2 : c
3 : c read frequency data from d/s ids. case depends on icase.
4 : c
5 : c icase = 1: grand summary. returned in cs(1-50). frequency from
6 : c variational frequency or, if not available from
7 : c eigenfrequency, possibly corrected for perturbation
8 : c in gravitational potential
9 : c icase = 2: short summary. returned in cs(1-7)
10 : c icase = 3: observed frequencies.
11 : c icase = 4: grand summary. returned in cs(1-50). sig and frq
12 : c obtained from eigenfrequency in cs(20).
13 : c Note that this allows setting Cowling
14 : c approximation frequency.
15 : c icase = 5: grand summary. returned in cs(1-50). frq obtained from
16 : c Richardson extrapolation frequency in cs(37),
17 : c if this is set. Otherwise variational frequency is used.
18 : c icase = 6: grand summary. returned in cs(1-50). sig and frq
19 : c obtained from eigenfrequency in cs(21),
20 : c hence possibly corrected for effect of
21 : c gravitational potential, for Cowling calculation.
22 : c
23 : c if icase .lt. 0, read from single-precision dataset, according
24 : c to value of abs(icase)
25 : c
26 : c returns mode degree, order, dimensionless sigma**2, frequency
27 : c (in microHz) and dimensionless energy in l, nord, sig, frq and ekin.
28 : c
29 : c Note: for icase = 1 returns corrected sigma**2, if frequency
30 : c calculation was in Cowling approximation. This is normally
31 : c also what is set in the short summary.
32 : c
33 : c ierr is returned as 1 or 2, on end of file or error in read.
34 : c
35 : c If data relevant for icase are not found, other data are used in the
36 : c order of preference: variational frequency, eigenfrequency. In that
37 : c case, ierr is returned as -icase_actual, where icase_actual is the
38 : c value of icase corresponding to the frequency used.
39 : c
40 : c original version: 26/9/86.
41 : c
42 : c Modified 1/1/08, to return negative ierr if proper data not
43 : c available. So far implemented only for icase = 5.
44 : c
45 : c ................................
46 : c
47 : c Double precision version
48 : c ++++++++++++++++++++++++
49 : c
50 : c Dated 10/3/90
51 : c
52 : implicit double precision (a-h, o-z)
53 : real css, sss
54 : dimension cs(*),csd(50),icsd(8), ssd(7),issd(2),
55 : * css(50),icss(8), sss(7),isss(2)
56 : c
57 : c common defining standard input and output
58 : c
59 : common/cstdio/ istdin, istdou, istdpr, istder
60 : equivalence(csd(39), icsd(1)),(ssd(6), issd(1))
61 : equivalence(css(39), icss(1)),(sss(6), isss(1))
62 : c
63 : data initrd, idsp /0, -1/
64 : data icerr /0/
65 : c
66 47 : cgrav=6.67232d-8
67 47 : twopi=8.d0*atan(1.d0)
68 : c
69 47 : icasea=iabs(icase)
70 : c
71 47 : ierr = 0
72 : c
73 47 : if(icasea.eq.2) then
74 : c
75 : c short summary
76 : c
77 0 : if(icase.eq.2) then
78 0 : read(ids,end=92,err=95) (cs(i),i=1,7)
79 : else
80 : c
81 : c single precision read
82 : c
83 0 : read(ids,end=92,err=95) (css(i),i=1,7)
84 0 : do 12 i=1,5
85 0 : 12 csd(i)=css(i)
86 0 : do 14 i=1,2
87 0 : 14 icsd(i)=icss(i)
88 0 : do 16 i=1,7
89 0 : 16 cs(i)=csd(i)
90 : end if
91 : c
92 0 : l=nint(cs(1))
93 0 : nord=nint(cs(2))
94 0 : sig=cs(3)
95 0 : frq=1000*cs(5)
96 0 : ekin=cs(4)
97 : c
98 0 : go to 90
99 : c
100 47 : else if(icasea.eq.3) then
101 : c
102 : c observed frequencies
103 : c
104 : c test for first read, file in the format for rdofrq
105 : c
106 0 : if(initrd.eq.0.or.ids.ne.idsp) then
107 0 : call skpcom(ids)
108 0 : read(ids,*,end=92,err=95) xlobsn
109 0 : if(xlobsn.lt.0) then
110 0 : read(ids,*,end=92,err=95) xnnobs
111 : else
112 0 : backspace ids
113 : end if
114 : c
115 0 : idsp=ids
116 0 : initrd=1
117 : end if
118 : c
119 0 : read(ids,*,end=92,err=95) l,nord,frq
120 0 : sig=0
121 0 : ekin=0
122 : c
123 0 : go to 90
124 : c
125 : end if
126 : c
127 : c other options assume grand summary input
128 : c
129 : c read grand summary
130 : c
131 47 : 20 if(icase.gt.0) then
132 2397 : 22 read(ids,end=92,err=95) (cs(i),i=1,50)
133 : c
134 : c test for proper summary (rather than x-record in mode file)
135 : c
136 46 : if(cs(2).le.1) go to 22
137 : else
138 : c
139 : c single precision read
140 : c
141 0 : 24 read(ids,end=92,err=95) (css(i),i=1,50)
142 : c
143 : c test for proper summary (rather than x-record in mode file)
144 : c
145 0 : if(css(2).le.1) go to 24
146 : c
147 0 : do 32 i=1,38
148 0 : 32 csd(i)=css(i)
149 0 : do 34 i=1,8
150 0 : 34 icsd(i)=icss(i)
151 0 : do 36 i=1,50
152 0 : 36 cs(i)=csd(i)
153 : c
154 : end if
155 : c
156 46 : l=nint(cs(18))
157 46 : nord=nint(cs(19))
158 46 : ekin=cs(24)
159 : c
160 : c setting of frequencies depend on value of icasea
161 : c
162 46 : if(icasea.eq.1) then
163 : c
164 : c use variational frequency, if available
165 : c
166 0 : sig=cs(21)
167 0 : if(sig.le.0) sig=cs(20)
168 0 : frq=1000*cs(27)
169 0 : if(frq.le.0) then
170 0 : frq=sqrt(sig/cs(20))*16666.6666667d0/cs(25)
171 0 : frq1=1.d6/twopi*sqrt(cgrav*cs(2)*sig/cs(3)**3)
172 0 : if(abs(frq/frq1-1).gt.1.e-5.and.icerr.le.20) then
173 0 : icerr=icerr+1
174 0 : if(istdpr.gt.0) write(istdpr,110) frq, frq1
175 0 : if(icerr.le.2.and.istdpr.gt.0)
176 0 : * write(istdpr,115) cs(2), cs(3), cs(20), cs(21),
177 0 : * cs(25)
178 : end if
179 : end if
180 : c
181 46 : else if(icasea.eq.4) then
182 : c
183 : c grand summary, frequency from uncorrected eigenfrequency
184 : c (note that this is what is given in cs(25)
185 : c include a test of frequency factor, for the time being,
186 : c for cases where cgrav .ne. 6.67232e-8
187 : c
188 0 : sig=cs(20)
189 0 : frq=16666.6666667d0/cs(25)
190 : c
191 46 : else if(icasea.eq.5) then
192 : c
193 : c grand summary, frequency from cs(37) (Richardson extrapolation)
194 : c if not available, use variational frequency or eigenfrequency
195 : c
196 0 : sig=cs(21)
197 0 : if(sig.le.0) sig=cs(20)
198 0 : if(cs(37).ne.0) then
199 0 : frq=1000*cs(37)
200 : else
201 0 : icerr=icerr+1
202 0 : if(icerr.le.20.and.istdpr.gt.0) write(istdpr,120)
203 0 : if(icerr.le.1) write(istder,120)
204 0 : if(cs(27).gt.0) then
205 0 : frq=1000*cs(27)
206 0 : ierr=-1
207 : else
208 0 : frq=16666.666667d0/cs(25)
209 0 : ierr=-1
210 : end if
211 : end if
212 : c
213 46 : else if(icasea.eq.6) then
214 : c
215 : c use eigenfrequency, possibly corrected for potential perturbation
216 : c
217 46 : sig=cs(21)
218 46 : if(sig.le.0) sig=cs(20)
219 46 : frq=sqrt(sig/cs(20))*16666.6666667d0/cs(25)
220 46 : frq1=1.d6/twopi*sqrt(cgrav*cs(2)*sig/cs(3)**3)
221 46 : if(abs(frq/frq1-1).gt.1.e-5.and.icerr.le.20) then
222 0 : icerr=icerr+1
223 0 : if(istdpr.gt.0) write(istdpr,110) icasea, frq, frq1
224 0 : if(icerr.le.2.and.istdpr.gt.0)
225 0 : * write(istdpr,115) cs(2), cs(3), cs(20), cs(21),
226 0 : * cs(25)
227 : end if
228 : c
229 : else
230 : c
231 0 : if(istdpr.gt.0) write(istdpr,130) icase
232 0 : stop
233 : end if
234 : c
235 : 90 continue
236 : c
237 : c
238 46 : return
239 : c
240 : c EOF or error
241 : c
242 1 : 92 ierr=1
243 1 : return
244 : c
245 0 : 95 ierr=2
246 0 : if(istdpr.gt.0)
247 0 : * write(istdpr,*) 'Error on reading d/s ',ids,' in s/r rdfreq'
248 0 : return
249 : 100 format(2i5,f8.2)
250 : 110 format(' *** warning. in s/r rdfreq, with icase = ',i1,
251 : * ', frq =',1pe13.6,' .ne. frq1 =',e13.6)
252 : 115 format(' cs(2), cs(3) =',1p2e15.6,' cs(20), cs(21), cs(25) =',
253 : * 3e15.6)
254 : 120 format(' *** error in s/r rdfreq with icase = 5. cs(37) not set')
255 : 130 format(//' ***** error. icase =',i5,' is illegal in rdfreq')
256 : end
257 0 : subroutine rdobsf(icase,ids,l,nord,frq,error,ierr)
258 : c
259 : c read observed frequency data from d/s ids. case depends on icase.
260 : c unused.
261 : c
262 : c icase = 3: no error data. each record consists of
263 : c l order frequency
264 : c
265 : c icase = 4: file contains error data. each record consists of
266 : c l order frequency error
267 : c
268 : c returns mode degree, order, frequency
269 : c (in microHz) and possibly error in l, nord, frq and error.
270 : c
271 : c ierr is returned as 1 or 2, on end of file or error in read.
272 : c
273 : c original version: 18/8/89.
274 : c
275 : c ................................
276 : c
277 : c for the moment skip test for rdofrq format. think about later
278 : c
279 : implicit double precision (a-h, o-z)
280 : c
281 : c common defining standard input and output
282 : c
283 : common/cstdio/ istdin, istdou, istdpr, istder
284 : data initrd /0/
285 : c
286 0 : ierr = 0
287 : c
288 : c test for first read, file in the format for rdofrq
289 : c
290 0 : if(initrd.eq.0) then
291 0 : call skpcom(ids)
292 0 : read(ids,*,end=92,err=95) xlobsn
293 0 : if(xlobsn.lt.0) then
294 0 : read(ids,*,end=92,err=95) xnnobs
295 : else
296 0 : backspace ids
297 : end if
298 : c
299 0 : initrd=1
300 : end if
301 : c
302 : c test for case
303 : c
304 0 : if(icase.eq.3) then
305 : c
306 0 : read(ids,*,end=92,err=95) l,nord,frq
307 0 : error=0
308 : else
309 0 : read(ids,*,end=92,err=95) l,nord,frq,error
310 : end if
311 : c
312 : 90 continue
313 : c
314 : c
315 0 : return
316 : c
317 : c EOF or error
318 : c
319 0 : 92 ierr=1
320 0 : return
321 : c
322 0 : 95 ierr=2
323 0 : if(istdpr.gt.0)
324 0 : * write(istdpr,*) 'Error on reading d/s ',ids,' in s/r rdobsf'
325 0 : return
326 : end
|
