| 1 |
program inv2d |
| 2 |
c 070701 Changed output format to accommodate l=1000 |
| 3 |
|
| 4 |
implicit double precision (a-h,o-z) |
| 5 |
include 'include1.f' |
| 6 |
include 'include2.f' |
| 7 |
|
| 8 |
parameter (maxpoints=maxr1*maxt1,nblock1=10) |
| 9 |
|
| 10 |
double precision split(maxnlm),spliti(maxl) |
| 11 |
double precision sigma(maxnlm),sigmai(maxl) |
| 12 |
double precision sigma1(maxnlm),ss1(maxnlm),ss2(maxnlm) |
| 13 |
double precision rmesh1(maxr1),tmesh1(maxt1) |
| 14 |
double precision rweights(maxr1),tweights(maxt1) |
| 15 |
double precision fr(maxr1),ft(maxr1),fwr(maxr1),fwt(maxt1) |
| 16 |
character*80 cefunc,outfile,splitfile |
| 17 |
real dummy(2),dtime,etime,time1 |
| 18 |
real mur,mut |
| 19 |
real f1(maxmodes,maxr1),f2(maxmodes,maxr1) |
| 20 |
real g1(maxlm,maxt1),g2(maxlm,maxt1) |
| 21 |
double precision f1x(maxmodes),f2x(maxmodes) |
| 22 |
double precision g1x(maxlm),g2x(maxlm) |
| 23 |
double precision f1i(maxr1),f2i(maxr1),g1i(maxt1),g2i(maxt1) |
| 24 |
double precision f1o(maxt1),f2o(maxt1) |
| 25 |
double precision a(0:maxs1-1,maxpoints) |
| 26 |
double precision ax(0:maxs1-1),bx(nblock,nblock) |
| 27 |
integer iax(maxpoints) |
| 28 |
double precision b(maxpoints,maxpoints),right(maxpoints) |
| 29 |
double precision c1(maxr1,maxr1),c2(maxt1,maxt1) |
| 30 |
double precision d(maxpoints,maxpoints),help(3*maxpoints) |
| 31 |
integer iwork(maxpoints) |
| 32 |
c double precision d1(maxpoints,maxpoints),cov(maxpoints,maxpoints) |
| 33 |
double precision var(maxpoints),dib(maxpoints) |
| 34 |
double precision dibx(maxpoints,nblock1) |
| 35 |
c double precision di(maxpoints) |
| 36 |
double precision coeff(maxnlm) |
| 37 |
double precision omega(maxr1,maxt1) |
| 38 |
double precision fro(maxr1,maxt1),fto(maxr1,maxt1) |
| 39 |
integer jav(maxpoints),kav(maxpoints) |
| 40 |
real avkern(maxtav,maxrav) |
| 41 |
integer nnlm1(maxnlm),llist(maxnlm) |
| 42 |
real fnlm(maxnlm) |
| 43 |
logical found |
| 44 |
integer ijk(maxr1,maxt1) |
| 45 |
logical ibig(maxpoints) |
| 46 |
c integer lsplit(maxnlm) |
| 47 |
c integer msplit(maxnlm),isplit(maxnlm),lmsplit(maxnlm) |
| 48 |
integer mindex(maxnlm) |
| 49 |
real x1(maxtav),x2(maxrav),yy2(maxtav,maxrav),mjta(4) |
| 50 |
real xr(10) |
| 51 |
integer ixl(maxmodes),ixh(maxmodes) |
| 52 |
|
| 53 |
|
| 54 |
c nmode,lmode: identification of modes |
| 55 |
c freq: frequency of mode (muHz) |
| 56 |
c radmesh: mesh in radius for kernels and eigenfunctions |
| 57 |
c nexp: number of expansion coefficients |
| 58 |
c cefunc: control file for reading of eigenfunctions |
| 59 |
c iomega: 0 to read omega from input |
| 60 |
c 1 to read omega from file |
| 61 |
c omegafile: file from which omega is read |
| 62 |
c nrin, ntin: number of points where omega is given |
| 63 |
c rin,tin,omegain: points where omega is given and the values there |
| 64 |
c omega is assumed constant outside the given interval |
| 65 |
|
| 66 |
open (20,file='status.log',form='formatted') |
| 67 |
write (20,*) 1 |
| 68 |
close (20) |
| 69 |
c read parameters etc. |
| 70 |
call clnfile(5,10) |
| 71 |
write (6,*) 'ierr, err' |
| 72 |
read (10,*) ierr,err |
| 73 |
write (6,*) 'name of file to read splittings from' |
| 74 |
read (10,'(A80)') splitfile |
| 75 |
write (6,*) 'control file for reading of eigenfunctions' |
| 76 |
read (10,'(A80)') cefunc |
| 77 |
write (6,*) 'number of points to skip in r for setting up f' |
| 78 |
read (10,*) nskipr |
| 79 |
write (6,*) 'number of points in x for setting up the plm''s' |
| 80 |
read (10,*) nsetp |
| 81 |
write (6,*) 'number of points to skip in x for setting up g' |
| 82 |
read (10,*) nskipt |
| 83 |
write (6,*) 'irdf, iwrf, irdg, iwrg, irdb, iwrb' |
| 84 |
read (10,*) irdf,iwrf,irdg,iwrg,irdb,iwrb |
| 85 |
write (6,*) '0 for multiple centerpoints, 1 for 1' |
| 86 |
read (10,*) icenter |
| 87 |
write (6,*) 'mur, mut' |
| 88 |
read (10,*) mur,mut |
| 89 |
nreg=2 |
| 90 |
write (6,*) 'icovar, icoeff, iavkern, nsr, nst, nsr2, nt2' |
| 91 |
read (10,*) icovar, icoeff, iavkern, nsr, nst, nsr2, nt2 |
| 92 |
if (nsr.eq.0) then |
| 93 |
nav=nst |
| 94 |
do 105,i=1,nav |
| 95 |
read (10,*) jav(i),kav(i) |
| 96 |
105 continue |
| 97 |
end if |
| 98 |
c write (6,*) 'name of output file' |
| 99 |
c read (10,'(A80)') outfile |
| 100 |
close(10) |
| 101 |
|
| 102 |
pi=4.0D0*atan(1.0D0) |
| 103 |
pi2=pi/2.0D0 |
| 104 |
write (6,*) 'start reading splittings' |
| 105 |
time1=dtime(dummy) |
| 106 |
open (11,file=splitfile,status='old') |
| 107 |
i=1 |
| 108 |
110 continue |
| 109 |
if (ierr.eq.0) then |
| 110 |
read (11,*,end=150) lnlm(i),nnlm1(i),fnlm(i), |
| 111 |
c mnlm(i),itnlm(i),it1nlm(i),split(i) |
| 112 |
sigma(i)=err |
| 113 |
else |
| 114 |
read (11,*,end=150) lnlm(i),nnlm1(i),fnlm(i), |
| 115 |
c mnlm(i),itnlm(i),it1nlm(i), |
| 116 |
c split(i),sigma(i) |
| 117 |
end if |
| 118 |
sigma1(i)=1.0/sigma(i) |
| 119 |
ss1(i)=split(i)*sigma1(i) |
| 120 |
ss2(i)=split(i)*sigma1(i)**2 |
| 121 |
i=i+1 |
| 122 |
goto 110 |
| 123 |
150 continue |
| 124 |
close(11) |
| 125 |
write (6,*) 'splittings have been read, time=', |
| 126 |
c dtime(dummy),dummy(1),dummy(2) |
| 127 |
nnlm=i-1 |
| 128 |
call setilm |
| 129 |
write (6,*) dtime(dummy),dummy(1),dummy(2) |
| 130 |
lmax1=-1 |
| 131 |
do 10,i=1,nnlm |
| 132 |
if (lnlm(i).gt.lmax1) lmax1=lnlm(i) |
| 133 |
10 continue |
| 134 |
i1=0 |
| 135 |
do 100,l=0,lmax1 |
| 136 |
nl=0 |
| 137 |
do 20,i=1,nnlm |
| 138 |
if (lnlm(i).eq.l) then |
| 139 |
nl=nl+1 |
| 140 |
llist(nl)=i |
| 141 |
end if |
| 142 |
20 continue |
| 143 |
i2=i1+1 |
| 144 |
do 60,i=1,nl |
| 145 |
il=llist(i) |
| 146 |
found=.false. |
| 147 |
do 30,ii=i2,i1 |
| 148 |
if (nnlm1(il).eq.nmode(ii)) then |
| 149 |
found=.true. |
| 150 |
ifound=ii |
| 151 |
end if |
| 152 |
30 continue |
| 153 |
if (found) then |
| 154 |
iln(il)=ifound |
| 155 |
else |
| 156 |
i1=i1+1 |
| 157 |
lmode(i1)=l |
| 158 |
nmode(i1)=nnlm1(il) |
| 159 |
freq(i1)=fnlm(il) |
| 160 |
iln(il)=i1 |
| 161 |
end if |
| 162 |
60 continue |
| 163 |
100 continue |
| 164 |
nmodes=i1 |
| 165 |
write (6,*) dtime(dummy),dummy(1),dummy(2) |
| 166 |
write (6,*) '#nl=',nmodes |
| 167 |
write (6,*) '#nlm=',nnlm |
| 168 |
c open (11,file=splitfile,form='unformatted',status='old') |
| 169 |
c i=1 |
| 170 |
c i1=0 |
| 171 |
c110 read (11,end=150) lmode(i),nmode(i),freq(i), |
| 172 |
c c (spliti(m),m=1,lmode(i)) |
| 173 |
c m0index(i)=i1 |
| 174 |
c do 130,m=1,lmode(i) |
| 175 |
c i1=i1+1 |
| 176 |
c split(i1)=spliti(m) |
| 177 |
c sigma(i1)=1.0E0 |
| 178 |
c sigma1(i1)=1.0D0/sigma(i1) |
| 179 |
c ss1(i1)=split(i1)*sigma1(i1) |
| 180 |
c lsplit(i1)=lmode(i) |
| 181 |
c msplit(i1)=m |
| 182 |
c isplit(i1)=i |
| 183 |
c130 continue |
| 184 |
c i=i+1 |
| 185 |
c goto 110 |
| 186 |
c150 nmodes=i-1 |
| 187 |
c nnlm=i1 |
| 188 |
c close (11) |
| 189 |
write (6,*) 'start reading eigenfunctions' |
| 190 |
time1=dtime(dummy) |
| 191 |
call readefun(cefunc) |
| 192 |
c do i=1,nrad |
| 193 |
c write (6,*) dradmesh(i) |
| 194 |
c end do |
| 195 |
write (6,*) 'eigenfunctions have been read, time=', |
| 196 |
c dtime(dummy),dummy(1),dummy(2) |
| 197 |
if (mod(nrad-1,nskipr).ne.0) then |
| 198 |
write (6,*) 'skip to get inversion mesh in radius does not' |
| 199 |
write (6,*) 'divide into number of kernel meshpoints-1' |
| 200 |
stop |
| 201 |
end if |
| 202 |
if (mod(nrad-1,nsr2).ne.0) then |
| 203 |
write (6,*) 'skip to get av. kernel mesh in radius does not' |
| 204 |
write (6,*) 'divide into number of kernel meshpoints-1' |
| 205 |
stop |
| 206 |
end if |
| 207 |
nrad1=(nrad-1)/nskipr+1 |
| 208 |
nr2=(nrad-1)/nsr2+1 |
| 209 |
do 160,i=1,nrad1 |
| 210 |
rmesh1(i)=radmesh(nskipr*(i-1)+1) |
| 211 |
160 continue |
| 212 |
call dsetint(1,nrad1,rmesh1,rweights) |
| 213 |
if (irdf.eq.0) then |
| 214 |
write (6,*) 'start setting up f matrix' |
| 215 |
time1=dtime(dummy) |
| 216 |
call setf(nskipr,f1,f2) |
| 217 |
write (6,*) 'f matrix has been set up, time=', |
| 218 |
c dtime(dummy),dummy(1),dummy(2) |
| 219 |
else |
| 220 |
write (6,*) 'start reading f matrix' |
| 221 |
time1=dtime(dummy) |
| 222 |
open (21,file=tmpdir//'savef',form='unformatted', |
| 223 |
c status='old') |
| 224 |
read (21) ((f1(i,j),i=1,nmodes),j=1,nrad1) |
| 225 |
read (21) ((f2(i,j),i=1,nmodes),j=1,nrad1) |
| 226 |
close(21) |
| 227 |
write (6,*) 'f matrix has been read, time=', |
| 228 |
c dtime(dummy),dummy(1),dummy(2) |
| 229 |
end if |
| 230 |
if (iwrf.ne.0) then |
| 231 |
write (6,*) 'start writing f matrix' |
| 232 |
time1=dtime(dummy) |
| 233 |
open (21,file=tmpdir//'savef',form='unformatted', |
| 234 |
c status='unknown') |
| 235 |
write (21) ((f1(i,j),i=1,nmodes),j=1,nrad1) |
| 236 |
write (21) ((f2(i,j),i=1,nmodes),j=1,nrad1) |
| 237 |
close(21) |
| 238 |
write (6,*) 'f matrix has been written, time=', |
| 239 |
c dtime(dummy),dummy(1),dummy(2) |
| 240 |
end if |
| 241 |
if (mod(nsetp-1,nskipt).ne.0) then |
| 242 |
write (6,*) 'skip to get av. kernel mesh in latitude does not' |
| 243 |
write (6,*) 'divide into original number of meshpoints-1' |
| 244 |
stop |
| 245 |
end if |
| 246 |
nsetp1=(nsetp-1)/nskipt+1 |
| 247 |
if (irdg.eq.0) then |
| 248 |
write (6,*) 'start setting up g matrix' |
| 249 |
time1=dtime(dummy) |
| 250 |
call setg(0,nsetp,nsetp1,g1,g2,tmesh1) |
| 251 |
write (6,*) 'g matrix has been set up, time=', |
| 252 |
c dtime(dummy),dummy(1),dummy(2) |
| 253 |
else |
| 254 |
write (6,*) 'start reading g matrix' |
| 255 |
time1=dtime(dummy) |
| 256 |
open (21,file=tmpdir//'saveg',form='unformatted', |
| 257 |
c status='old') |
| 258 |
read (21) (tmesh1(i),i=1,nsetp1) |
| 259 |
c read (21) nlm,(lplm(i),mplm(i),i=1,nlm) |
| 260 |
cc read (21) lmax,((iplm(i,j),i=0,lmax),j=0,lmax) |
| 261 |
read (21) ((g1(i,j),i=1,nlm),j=1,nsetp1) |
| 262 |
read (21) ((g2(i,j),i=1,nlm),j=1,nsetp1) |
| 263 |
close(21) |
| 264 |
write (6,*) 'g matrix has been read, time=', |
| 265 |
c dtime(dummy),dummy(1),dummy(2) |
| 266 |
end if |
| 267 |
if (iwrg.ne.0) then |
| 268 |
write (6,*) 'start writing g matrix' |
| 269 |
time1=dtime(dummy) |
| 270 |
open (21,file=tmpdir//'saveg',form='unformatted', |
| 271 |
c status='unknown') |
| 272 |
write (21) (tmesh1(i),i=1,nsetp1) |
| 273 |
c write (21) nlm,(lplm(i),mplm(i),i=1,nlm) |
| 274 |
cc write (21) lmax,((iplm(i,j),i=0,lmax),j=0,lmax) |
| 275 |
write (21) ((g1(i,j),i=1,nlm),j=1,nsetp1) |
| 276 |
write (21) ((g2(i,j),i=1,nlm),j=1,nsetp1) |
| 277 |
close(21) |
| 278 |
write (6,*) 'g matrix has been written, time=', |
| 279 |
c dtime(dummy),dummy(1),dummy(2) |
| 280 |
end if |
| 281 |
call dsetint(1,nsetp1,tmesh1,tweights) |
| 282 |
npoints=nrad1*nsetp1 |
| 283 |
nident=icenter*(nsetp1-1) |
| 284 |
npoi1=npoints-nident |
| 285 |
do 500,i2=1,npoints |
| 286 |
k2=mod(i2-1,nsetp1)+1 |
| 287 |
j2=(i2-1)/nsetp1+1 |
| 288 |
ijk(j2,k2)=max(1,i2-nident) |
| 289 |
500 continue |
| 290 |
do 600,i=1,nnlm |
| 291 |
c lmsplit(i)=iplm(msplit(i),lsplit(i)) |
| 292 |
600 continue |
| 293 |
open (24,file='rmesh.orig',form='formatted') |
| 294 |
do 610,i=1,nrad |
| 295 |
write (24,*) radmesh(i) |
| 296 |
610 continue |
| 297 |
close (24) |
| 298 |
if (iavkern.ne.0) then |
| 299 |
open (24,file='avkern.log',form='formatted') |
| 300 |
end if |
| 301 |
open (25,file='rot.2d',form='formatted') |
| 302 |
if (icovar.eq.1) then |
| 303 |
open (26,file='err.2d',form='formatted') |
| 304 |
end if |
| 305 |
open (27,file='avkern.2d',form='unformatted') |
| 306 |
if (icoeff.eq.1) then |
| 307 |
open (28,file='coeff.2d',form='unformatted') |
| 308 |
end if |
| 309 |
xr(1)=-9.5 |
| 310 |
xr(2)=-8.5 |
| 311 |
xr(3)=-3.5 |
| 312 |
xr(4)=-3.0 |
| 313 |
xr(5)=-2.5 |
| 314 |
xr(6)=-2.0 |
| 315 |
c do 5000,imr=-8,-4,1 |
| 316 |
c mur=10.**imr |
| 317 |
c do 5000,imr=1,6 |
| 318 |
c mur=10.**xr(imr) |
| 319 |
c do 5000,imt=-7,-1,1 |
| 320 |
c mut=10.**imt |
| 321 |
c write (6,*) mur,mut |
| 322 |
if (irdb.eq.0) then |
| 323 |
write (6,*) 'start setting up b matrix' |
| 324 |
time1=dtime(dummy) |
| 325 |
do 950,i=1,npoi1 |
| 326 |
do 940,j=1,npoi1 |
| 327 |
b(j,i)=0.0D0 |
| 328 |
940 continue |
| 329 |
950 continue |
| 330 |
nb1=0 |
| 331 |
nb2=0 |
| 332 |
nb2a=0 |
| 333 |
t2=0 |
| 334 |
t2a=0 |
| 335 |
do 1100,is1=1,nnlm,maxs1 |
| 336 |
write (6,*) is1 |
| 337 |
nbig=0 |
| 338 |
nbig2=0 |
| 339 |
is0=min(maxs1-1,nnlm-is1) |
| 340 |
ix=0 |
| 341 |
do 1020,i2=1,npoints |
| 342 |
ibig(i2)=.false. |
| 343 |
k2=mod(i2-1,nsetp1)+1 |
| 344 |
j2=(i2-1)/nsetp1+1 |
| 345 |
nbig1=0 |
| 346 |
do 1010,isx=0,is0 |
| 347 |
is=is1+isx |
| 348 |
i=iln(is) |
| 349 |
ilmi=ilm(is) |
| 350 |
ax(isx)=sigma1(is)* |
| 351 |
c (f1(i,j2)*g1(ilmi,k2)+f2(i,j2)*g2(ilmi,k2)) |
| 352 |
if (abs(ax(isx)).gt.1e-15) then |
| 353 |
nbig1=nbig1+1 |
| 354 |
end if |
| 355 |
1010 continue |
| 356 |
if (nbig1.ne.0) then |
| 357 |
ibig(i2)=.true. |
| 358 |
nbig=nbig+nbig1 |
| 359 |
nbig2=nbig2+1 |
| 360 |
ix=ix+1 |
| 361 |
iax(ix)=i2 |
| 362 |
do 1015,isx=0,is0 |
| 363 |
a(isx,ix)=ax(isx) |
| 364 |
1015 continue |
| 365 |
end if |
| 366 |
1020 continue |
| 367 |
nx=ix |
| 368 |
nb1=nb1+nbig |
| 369 |
nb2=nb2+nbig2 |
| 370 |
nb2a=nb2a+npoints |
| 371 |
t2=t2+nbig2**2*(is0+1.) |
| 372 |
t2a=t2a+npoints**2*(is0+1.) |
| 373 |
do 1080,i1x=1,nx,nblock |
| 374 |
n1x=min(i1x+nblock-1,nx)-i1x+1 |
| 375 |
do 1070,i2x=i1x,nx,nblock |
| 376 |
n2x=min(i2x+nblock-1,nx)-i2x+1 |
| 377 |
call dgemm('t','n',n2x,n1x,is0+1,1.0D0,a(0,i2x),maxs1, |
| 378 |
c a(0,i1x),maxs1,0.0D0,bx,nblock) |
| 379 |
do 1060,i1=i1x,min(i1x+nblock-1,nx) |
| 380 |
i1a=max(1,iax(i1)-nident) |
| 381 |
do 1050,i2=i2x,min(i2x+nblock-1,nx) |
| 382 |
i2a=max(1,iax(i2)-nident) |
| 383 |
if (i2a.ge.i1a) then |
| 384 |
c b(i2a,i1a)=b(i2a,i1a)+ddot(is0+1,a(0,i2),1,a(0,i1),1) |
| 385 |
b(i2a,i1a)=b(i2a,i1a)+bx(i2-i2x+1,i1-i1x+1) |
| 386 |
end if |
| 387 |
1050 continue |
| 388 |
1060 continue |
| 389 |
1070 continue |
| 390 |
1080 continue |
| 391 |
1100 continue |
| 392 |
write (6,*) (nb1+0.0)/(nnlm*npoints),(nb2+0.0)/nb2a,t2/t2a |
| 393 |
do 1290,i2=1,npoi1 |
| 394 |
do 1280,i1=1,i2-1 |
| 395 |
b(i1,i2)=b(i2,i1) |
| 396 |
1280 continue |
| 397 |
1290 continue |
| 398 |
write (6,*) 'b matrix has been set up, time=', |
| 399 |
c dtime(dummy),dummy(1),dummy(2) |
| 400 |
else |
| 401 |
write (6,*) 'start reading b matrix' |
| 402 |
time1=dtime(dummy) |
| 403 |
open (21,file=tmpdir//'saveb',form='unformatted', |
| 404 |
c status='old') |
| 405 |
read (21) ((b(j,i),j=1,npoi1),i=1,npoi1) |
| 406 |
close(21) |
| 407 |
c do i=1,npoi1 |
| 408 |
c do j=1,npoi1 |
| 409 |
c b(j,i)=real(b(j,i)) |
| 410 |
c end do |
| 411 |
c end do |
| 412 |
write (6,*) 'b matrix has been read, time=', |
| 413 |
c dtime(dummy),dummy(1),dummy(2) |
| 414 |
end if |
| 415 |
c write (6,*) (b(i,i),i=1,npoi1) |
| 416 |
if (iwrb.ne.0) then |
| 417 |
write (6,*) 'start writing b matrix' |
| 418 |
time1=dtime(dummy) |
| 419 |
open (21,file=tmpdir//'saveb',form='unformatted', |
| 420 |
c status='unknown') |
| 421 |
write (21) ((b(j,i),j=1,npoi1),i=1,npoi1) |
| 422 |
close(21) |
| 423 |
write (6,*) 'b matrix has been written, time=', |
| 424 |
c dtime(dummy),dummy(1),dummy(2) |
| 425 |
c to prevent recalculating b if doing trade-off curves |
| 426 |
irdb=1 |
| 427 |
end if |
| 428 |
write (6,*) 'start adding regularization' |
| 429 |
time1=dtime(dummy) |
| 430 |
do 1310,i=1,npoi1 |
| 431 |
do 1300,j=1,npoi1 |
| 432 |
d(j,i)=b(j,i) |
| 433 |
1300 continue |
| 434 |
1310 continue |
| 435 |
c c1 is is the regularization in the 1-d radial case |
| 436 |
c times the radius |
| 437 |
if (nreg.eq.2) then |
| 438 |
c set up weighting functions |
| 439 |
do 1400,i=1,nrad1 |
| 440 |
c fr(i)=1.0 |
| 441 |
fr(i)=rmesh1(i) |
| 442 |
if (rmesh1(i).eq.0) then |
| 443 |
ft(i)=0.0 |
| 444 |
else |
| 445 |
ft(i)=1.0D0/rmesh1(i)**4 |
| 446 |
end if |
| 447 |
c just a try |
| 448 |
c ft(i)=1.0 |
| 449 |
1400 continue |
| 450 |
do 1405,i=1,nrad1 |
| 451 |
fwr(i)=mur*rweights(i)*fr(i) |
| 452 |
c fwr(i)=mur*rweights(i)*rmesh1(i) |
| 453 |
1405 continue |
| 454 |
do 1430,j=1,nrad1 |
| 455 |
do 1420,k=j,min(nrad1,j+2) |
| 456 |
ckj=0.0D0 |
| 457 |
do 1410,i=max(2,j-1,k-1),min(nrad1-1,j+1,k+1) |
| 458 |
xm1=rmesh1(i-1)-rmesh1(i) |
| 459 |
xp1=rmesh1(i+1)-rmesh1(i) |
| 460 |
den=xm1*xp1*(xm1-xp1) |
| 461 |
if (i.eq.j-1) fji=-xm1 |
| 462 |
if (i.eq.j) fji=xm1-xp1 |
| 463 |
if (i.eq.j+1) fji=xp1 |
| 464 |
if (i.eq.k-1) fki=-xm1 |
| 465 |
if (i.eq.k) fki=xm1-xp1 |
| 466 |
if (i.eq.k+1) fki=xp1 |
| 467 |
ckj=ckj+fji*fki*fwr(i)/den**2 |
| 468 |
1410 continue |
| 469 |
c factor of 4=2^2 since second derivative is 2*a where y=a x^2+b x +c |
| 470 |
c1(k,j)=4.*ckj |
| 471 |
c1(j,k)=4.*ckj |
| 472 |
1420 continue |
| 473 |
1430 continue |
| 474 |
do 1500,i=1,nsetp1 |
| 475 |
fwt(i)=mut*tweights(i) |
| 476 |
1500 continue |
| 477 |
do 1530,j=1,nsetp1 |
| 478 |
do 1520,k=j,min(nsetp1,j+2) |
| 479 |
ckj=0.0D0 |
| 480 |
do 1510,i=max(1,j-1,k-1),min(nsetp1,j+1,k+1) |
| 481 |
t=tmesh1(i) |
| 482 |
if (i.gt.1) tm1=tmesh1(i-1) |
| 483 |
if (i.eq.1) tm1=2*tmesh1(1)-tmesh1(2) |
| 484 |
if (i.lt.nsetp1) tp1=tmesh1(i+1) |
| 485 |
if (i.eq.nsetp1) tp1=2*tmesh1(nsetp1)-tmesh1(nsetp1-1) |
| 486 |
xm1=tm1-t |
| 487 |
xp1=tp1-t |
| 488 |
den=xm1*xp1*(xm1-xp1) |
| 489 |
if (i.eq.j-1) fji=-xm1 |
| 490 |
if (i.eq.j) fji=xm1-xp1 |
| 491 |
if (i.eq.j+1) fji=xp1 |
| 492 |
if (i.eq.1.and.j.eq.2) fji=xp1-xm1 |
| 493 |
if (i.eq.nsetp1.and.j.eq.(nsetp1-1)) fji=xp1-xm1 |
| 494 |
if (i.eq.k-1) fki=-xm1 |
| 495 |
if (i.eq.k) fki=xm1-xp1 |
| 496 |
if (i.eq.k+1) fki=xp1 |
| 497 |
if (i.eq.1.and.k.eq.2) fki=xp1-xm1 |
| 498 |
if (i.eq.nsetp1.and.k.eq.(nsetp1-1)) fki=xp1-xm1 |
| 499 |
ckj=ckj+fji*fki*fwt(i)/den**2 |
| 500 |
c write (6,*) i,j,k,ckj |
| 501 |
1510 continue |
| 502 |
c factor of 4=2^2 since second derivative is 2*a where y=a x^2+b x +c |
| 503 |
c2(k,j)=4.*ckj |
| 504 |
c2(j,k)=4.*ckj |
| 505 |
1520 continue |
| 506 |
1530 continue |
| 507 |
end if |
| 508 |
write (6,*) dtime(dummy),dummy(1),dummy(2) |
| 509 |
c first the regularization in theta is applied |
| 510 |
do 1650,j=2,nrad1 |
| 511 |
wr3i=rweights(j)*ft(j) |
| 512 |
c wr3i=rweights(j)/rmesh1(j)**4 |
| 513 |
do 1640,k1=1,nsetp1 |
| 514 |
i1=ijk(j,k1) |
| 515 |
do 1630,k2=max(1,k1-2),min(nsetp1,k1+2) |
| 516 |
i2=ijk(j,k2) |
| 517 |
d(i2,i1)=d(i2,i1)+c2(k2,k1)*wr3i |
| 518 |
1630 continue |
| 519 |
1640 continue |
| 520 |
1650 continue |
| 521 |
write (6,*) dtime(dummy),dummy(1),dummy(2) |
| 522 |
c then the regularization in radius is added |
| 523 |
do 1750,k=1,nsetp1 |
| 524 |
w=tweights(k) |
| 525 |
do 1740,j1=1,nrad1 |
| 526 |
i1=ijk(j1,k) |
| 527 |
do 1730,j2=1,nrad1 |
| 528 |
i2=ijk(j2,k) |
| 529 |
d(i2,i1)=d(i2,i1)+c1(j2,j1)*w |
| 530 |
1730 continue |
| 531 |
1740 continue |
| 532 |
1750 continue |
| 533 |
write (6,*) 'regularization has been added, time=', |
| 534 |
c dtime(dummy),dummy(1),dummy(2) |
| 535 |
write (6,*) 'start setting up right hand side' |
| 536 |
time1=dtime(dummy) |
| 537 |
do 1810,i=1,npoi1 |
| 538 |
right(i)=0.0D0 |
| 539 |
1810 continue |
| 540 |
do 1850,i1=1,npoints |
| 541 |
c i1a=max(1,i1-nident) |
| 542 |
k1=mod(i1-1,nsetp1)+1 |
| 543 |
j1=(i1-1)/nsetp1+1 |
| 544 |
do 1820,i=1,nmodes |
| 545 |
f1x(i)=f1(i,j1) |
| 546 |
f2x(i)=f2(i,j1) |
| 547 |
1820 continue |
| 548 |
do 1830,i=1,nlm |
| 549 |
g1x(i)=g1(i,k1) |
| 550 |
g2x(i)=g2(i,k1) |
| 551 |
1830 continue |
| 552 |
r=0.0D0 |
| 553 |
do 1840,is=1,nnlm |
| 554 |
i=iln(is) |
| 555 |
ilmi=ilm(is) |
| 556 |
r=r+ss2(is)*(f1x(i)*g1x(ilmi)+f2x(i)*g2x(ilmi)) |
| 557 |
1840 continue |
| 558 |
c right(i1a)=right(i1a)+r |
| 559 |
help(i1)=r |
| 560 |
1850 continue |
| 561 |
do 1860,i1=1,npoints |
| 562 |
i1a=max(1,i1-nident) |
| 563 |
right(i1a)=right(i1a)+help(i1) |
| 564 |
1860 continue |
| 565 |
write (6,*) 'right hand side has been set up, time=', |
| 566 |
c dtime(dummy),dummy(1),dummy(2) |
| 567 |
c open (21,file=tmpdir//'saved',form='unformatted', |
| 568 |
c c status='unknown') |
| 569 |
c do i=1,npoi1 |
| 570 |
c write (21) (d(j,i),j=1,npoi1) |
| 571 |
c end do |
| 572 |
c write (21) (right(i),i=1,npoi1) |
| 573 |
c close(21) |
| 574 |
c stop |
| 575 |
write (6,*) 'start decomposing d matrix' |
| 576 |
time1=dtime(dummy) |
| 577 |
c call dpoco(d,maxpoints,npoi1,rcond,help,info) |
| 578 |
anorm=dlansy('1','U',npoi1,d,maxpoints,help) |
| 579 |
call dpotrf('U',npoi1,d,maxpoints,info) |
| 580 |
write (6,*) 'd matrix has been decomposed, time=', |
| 581 |
c dtime(dummy),dummy(1),dummy(2),info |
| 582 |
write (6,*) 'start finding condition number' |
| 583 |
call dpocon('U',npoi1,d,maxpoints,anorm,rcond,help,iwork,info1) |
| 584 |
write (6,*) 'condition number has been found, time=', |
| 585 |
c dtime(dummy),dummy(1),dummy(2) |
| 586 |
write (6,*) rcond,info,info1,anorm |
| 587 |
write (6,*) 'start finding solution' |
| 588 |
time1=dtime(dummy) |
| 589 |
c call dposl(d,maxpoints,npoi1,right) |
| 590 |
call dpotrs('U',npoi1,1,d,maxpoints,right,maxpoints,info) |
| 591 |
write (6,*) 'solution has been found, time=', |
| 592 |
c dtime(dummy),dummy(1),dummy(2) |
| 593 |
c write (21) (right(i),i=1,npoi1) |
| 594 |
c close(21) |
| 595 |
c stop |
| 596 |
write (6,*) 'rcond,info',rcond,info |
| 597 |
do 1875,j=1,nrad1 |
| 598 |
do 1870,k=1,nsetp1 |
| 599 |
omega(j,k)=right(ijk(j,k)) |
| 600 |
1870 continue |
| 601 |
1875 continue |
| 602 |
do 1877,j=1,nrad1 |
| 603 |
write (25,'(101f12.4)') (right(ijk(j,k)),k=1,nsetp1) |
| 604 |
1877 continue |
| 605 |
write (6,*) 'start calculating splittings' |
| 606 |
time1=dtime(dummy) |
| 607 |
c set ixl to the first mode with a given (l,n) |
| 608 |
c set ixh to the last mode with a given (l,n) |
| 609 |
ilast=-1 |
| 610 |
ix=0 |
| 611 |
do 1880,i1=1,nnlm |
| 612 |
i=iln(i1) |
| 613 |
if (i.ne.ilast) then |
| 614 |
ix=ix+1 |
| 615 |
ixl(ix)=i1 |
| 616 |
ilast=i |
| 617 |
end if |
| 618 |
ixh(ix)=i1 |
| 619 |
1880 continue |
| 620 |
nx=ix |
| 621 |
chisq=0.0 |
| 622 |
c*$* assert do (serial) |
| 623 |
open(10,file='splittings.out',form='formatted') |
| 624 |
do 1897,ix=1,nx |
| 625 |
c loop first over each ln |
| 626 |
i=iln(ixl(ix)) |
| 627 |
do 1891,j=1,nrad1 |
| 628 |
f1i(j)=f1(i,j) |
| 629 |
f2i(j)=f2(i,j) |
| 630 |
1891 continue |
| 631 |
call dgemv('t',nrad1,nsetp1,1.0D0,omega,maxr1, |
| 632 |
c f1i,1,0.0D0,f1o,1) |
| 633 |
call dgemv('t',nrad1,nsetp1,1.0D0,omega,maxr1, |
| 634 |
c f2i,1,0.0D0,f2o,1) |
| 635 |
c now do the individual nlm's |
| 636 |
do 1895,i1=ixl(ix),ixh(ix) |
| 637 |
index=ilm(i1) |
| 638 |
c do 1893,k=1,nsetp1 |
| 639 |
c g1i(k)=g1(index,k) |
| 640 |
c g2i(k)=g2(index,k) |
| 641 |
c1893 continue |
| 642 |
c s=ddot(nsetp1,f1o,1,g1i,1)+ddot(nsetp1,f2o,1,g2i,1) |
| 643 |
s=0.0D0 |
| 644 |
do 1893,k=1,nsetp1 |
| 645 |
s=s+g1(index,k)*f1o(k)+g2(index,k)*f2o(k) |
| 646 |
1893 continue |
| 647 |
l=lmode(i) |
| 648 |
n=nmode(i) |
| 649 |
write (10,'(i4,i4,f8.2,i4,2i3,3f12.4)') |
| 650 |
c write (6,'(i4,i4,f8.2,i4,2i3,3f12.4)') |
| 651 |
c l,n,freq(i),mnlm(i1),itnlm(i1),it1nlm(i1), |
| 652 |
c s,split(i1),sigma(i1) |
| 653 |
chisq=chisq+((s-split(i1))*sigma1(i1))**2 |
| 654 |
1895 continue |
| 655 |
1897 continue |
| 656 |
write (6,*) 'splittings have been calculated, time=', |
| 657 |
c dtime(dummy),dummy(1),dummy(2) |
| 658 |
write (6,*) 'start calculating wiggle measures' |
| 659 |
time1=dtime(dummy) |
| 660 |
do 1907,k=1,nsetp1 |
| 661 |
do 1905,j=1,nrad1 |
| 662 |
fro(j,k)=0.0 |
| 663 |
fto(j,k)=0.0 |
| 664 |
c omega(j,k)=sin(tmesh1(k))**2*rmesh1(j)**2 |
| 665 |
c omega(j,k)=rmesh1(j)**3 |
| 666 |
1905 continue |
| 667 |
1907 continue |
| 668 |
do 1930,j=1,nrad1 |
| 669 |
do 1920,i=max(2,j-1),min(nrad1-1,j+1) |
| 670 |
xm1=rmesh1(i-1)-rmesh1(i) |
| 671 |
xp1=rmesh1(i+1)-rmesh1(i) |
| 672 |
den=xm1*xp1*(xm1-xp1) |
| 673 |
if (i.eq.j-1) fji=-xm1 |
| 674 |
if (i.eq.j) fji=xm1-xp1 |
| 675 |
if (i.eq.j+1) fji=xp1 |
| 676 |
do 1910,k=1,nsetp1 |
| 677 |
fro(i,k)=fro(i,k)+omega(j,k)*2*fji/den |
| 678 |
1910 continue |
| 679 |
1920 continue |
| 680 |
1930 continue |
| 681 |
do 1960,k=1,nsetp1 |
| 682 |
do 1950,i=max(1,k-1),min(nsetp1,k+1) |
| 683 |
t=tmesh1(i) |
| 684 |
if (i.gt.1) tm1=tmesh1(i-1) |
| 685 |
if (i.eq.1) tm1=2*tmesh1(1)-tmesh1(2) |
| 686 |
if (i.lt.nsetp1) tp1=tmesh1(i+1) |
| 687 |
if (i.eq.nsetp1) tp1=2*tmesh1(nsetp1)-tmesh1(nsetp1-1) |
| 688 |
xm1=tm1-t |
| 689 |
xp1=tp1-t |
| 690 |
den=xm1*xp1*(xm1-xp1) |
| 691 |
if (i.eq.k-1) fki=-xm1 |
| 692 |
if (i.eq.k) fki=xm1-xp1 |
| 693 |
if (i.eq.k+1) fki=xp1 |
| 694 |
if (i.eq.1.and.k.eq.2) fki=xp1-xm1 |
| 695 |
if (i.eq.nsetp1.and.k.eq.(nsetp1-1)) fki=xp1-xm1 |
| 696 |
do 1940,j=1,nrad1 |
| 697 |
fto(j,i)=fto(j,i)+omega(j,k)*2*fki/den |
| 698 |
1940 continue |
| 699 |
1950 continue |
| 700 |
1960 continue |
| 701 |
c calculate wigr and wigt |
| 702 |
wigr=0.0 |
| 703 |
wigt=0.0 |
| 704 |
do 1980,k=1,nsetp1 |
| 705 |
wt=tweights(k) |
| 706 |
do 1970,j=1,nrad1 |
| 707 |
wigr=wigr+wt*rweights(j)*fr(j)*fro(j,k)**2 |
| 708 |
wigt=wigt+wt*rweights(j)*ft(j)*fto(j,k)**2 |
| 709 |
1970 continue |
| 710 |
1980 continue |
| 711 |
write (6,*) 'wiggle measures have been computed, time=', |
| 712 |
c dtime(dummy),dummy(1),dummy(2) |
| 713 |
write (6,*) 'chisq=',chisq |
| 714 |
write (6,*) 'rms=',sqrt(chisq/nnlm) |
| 715 |
write (6,*) 'wigr, wigt=',wigr,wigt |
| 716 |
write (6,*) '%',mur,mut,chisq,sqrt(chisq/nnlm),wigr,wigt, |
| 717 |
c info,rcond |
| 718 |
c end do |
| 719 |
c end do |
| 720 |
c stop |
| 721 |
if ((icovar.eq.1).or.(iavkern.ge.1)) then |
| 722 |
write (6,*) 'start inverting d matrix' |
| 723 |
time1=dtime(dummy) |
| 724 |
c call dpodi(d,maxpoints,npoi1,det,1) |
| 725 |
call dpotri('U',npoi1,d,maxpoints,info) |
| 726 |
c open (21,file=tmpdir//'saved1',form='unformatted', |
| 727 |
c c status='unknown',access='direct',recl=8*npoi1) |
| 728 |
do 2095,i=1,npoi1 |
| 729 |
c do 2092,j=1,npoi1 |
| 730 |
c di(j)=0.0D0 |
| 731 |
c2092 continue |
| 732 |
c di(i)=1.0D0 |
| 733 |
c call dposl(d,maxpoints,npoi1,di) |
| 734 |
c write (21,rec=i) (di(j),j=1,npoi1) |
| 735 |
c write (21,rec=i) (d(j,i),j=1,i),(d(i,j),j=i+1,npoi1) |
| 736 |
do 2092,j=i+1,npoi1 |
| 737 |
d(j,i)=d(i,j) |
| 738 |
2092 continue |
| 739 |
2095 continue |
| 740 |
c close(21) |
| 741 |
write (6,*) 'd matrix has been inverted, time=', |
| 742 |
c dtime(dummy),dummy(1),dummy(2) |
| 743 |
end if |
| 744 |
if (icovar.ge.1) then |
| 745 |
write (6,*) 'start finding variances' |
| 746 |
c the following loop could be blocked |
| 747 |
c do 2150,i=1,npoi1 |
| 748 |
c call dgemv('n',npoi1,npoi1,1.0D0,b,maxpoints,d(1,i),1, |
| 749 |
c c 0.0D0,dib,1) |
| 750 |
c var(i)=ddot(npoi1,d(1,i),1,dib,1) |
| 751 |
c2150 continue |
| 752 |
do 2150,i=1,npoi1,nblock1 |
| 753 |
ix=min(npoi1,i+nblock1-1) |
| 754 |
nx=ix+1-i |
| 755 |
call dgemm('n','n',npoi1,nx,npoi1,1.0D0,b,maxpoints, |
| 756 |
c d(1,i),maxpoints,0.0D0,dibx,maxpoints) |
| 757 |
do 2140,j=i,ix |
| 758 |
var(j)=ddot1(npoi1,d(1,j),dibx(1,j-i+1)) |
| 759 |
2140 continue |
| 760 |
2150 continue |
| 761 |
do 2190,j=1,nrad1 |
| 762 |
write (26,'(100f11.5)') (sqrt(var(ijk(j,k))),k=1,nsetp1) |
| 763 |
2190 continue |
| 764 |
c close(21) |
| 765 |
write (6,*) 'variances have been found, time=', |
| 766 |
c dtime(dummy),dummy(1),dummy(2) |
| 767 |
c else |
| 768 |
c write (6,*) 'start finding covariance matrix' |
| 769 |
c time1=dtime(dummy) |
| 770 |
c do 2230,i=1,npoi1 |
| 771 |
c do 2220,j=1,npoi1 |
| 772 |
c d(j,i)=ddot1(npoi1,b(1,j),d1(1,i)) |
| 773 |
c2220 continue |
| 774 |
c2230 continue |
| 775 |
c do 2250,i=1,npoi1 |
| 776 |
c do 2240,j=i,npoi1 |
| 777 |
c cov(j,i)=ddot1(npoi1,d1(1,j),d(1,i)) |
| 778 |
c2240 continue |
| 779 |
c2250 continue |
| 780 |
c do 2270,i=1,npoi1 |
| 781 |
c do 2260,j=1,i-1 |
| 782 |
c cov(j,i)=cov(i,j) |
| 783 |
c2260 continue |
| 784 |
c2270 continue |
| 785 |
c do 2280,i=1,nrad1 |
| 786 |
c write (6,*) rmesh1(i) |
| 787 |
c2280 continue |
| 788 |
c do 2290,j=1,nrad1 |
| 789 |
c write (6,'(100f9.5)') (sqrt(cov(ijk(j,k),ijk(j,k))), |
| 790 |
c c k=1,nsetp1) |
| 791 |
c2290 continue |
| 792 |
c write (6,*) 'covariance matrix has been found, time=', |
| 793 |
c c dtime(dummy),dummy(1),dummy(2) |
| 794 |
end if |
| 795 |
if (iavkern.eq.0) then |
| 796 |
write (6,*) 'skipping calculation of averaging kernels' |
| 797 |
else |
| 798 |
write (6,*) 'start finding averaging kernels' |
| 799 |
time1=dtime(dummy) |
| 800 |
write (6,*) nt2,nr2 |
| 801 |
c open (21,file=tmpdir//'saved1',form='unformatted', |
| 802 |
c c status='old',access='direct',recl=8*npoi1) |
| 803 |
c moved up |
| 804 |
c open (27,file='avkern.2d',form='unformatted') |
| 805 |
c open (28,file='coeff',form='unformatted') |
| 806 |
do 2300,j1=1,nr2 |
| 807 |
x2(j1)=radmesh(nsr2*(j1-1)+1) |
| 808 |
2300 continue |
| 809 |
pi=4.0D0*atan(1.0D0) |
| 810 |
pi2=pi/2.0D0 |
| 811 |
do 2301,k1=1,nt2 |
| 812 |
x1(k1)=pi2*(k1-1)/(nt2-1.) |
| 813 |
2301 end do |
| 814 |
write (6,*) nrad1,nsr,nsetp1,nst |
| 815 |
if (nsr.ne.0) then |
| 816 |
i=0 |
| 817 |
do 2303,j=1,nrad1,nsr |
| 818 |
do 2302,k=1,nsetp1,nst |
| 819 |
i=i+1 |
| 820 |
jav(i)=j |
| 821 |
kav(i)=k |
| 822 |
2302 continue |
| 823 |
2303 continue |
| 824 |
nav=i |
| 825 |
end if |
| 826 |
c do 2380,j=1,nrad1,nsr |
| 827 |
c do 2370,k=1,nsetp1,nst |
| 828 |
do 2380,iav=1,nav |
| 829 |
j=jav(iav) |
| 830 |
k=kav(iav) |
| 831 |
ipoi=ijk(j,k) |
| 832 |
c read (21,rec=ipoi) (di(j1),j1=1,npoi1) |
| 833 |
do 2305,is=1,nnlm |
| 834 |
coeff(is)=0.0D0 |
| 835 |
2305 continue |
| 836 |
do 2330,i1=1,npoints |
| 837 |
i1a=max(1,i1-nident) |
| 838 |
k1=mod(i1-1,nsetp1)+1 |
| 839 |
j1=(i1-1)/nsetp1+1 |
| 840 |
do 2310,i=1,nmodes |
| 841 |
f1x(i)=f1(i,j1) |
| 842 |
f2x(i)=f2(i,j1) |
| 843 |
2310 continue |
| 844 |
do 2315,i=1,nlm |
| 845 |
g1x(i)=g1(i,k1) |
| 846 |
g2x(i)=g2(i,k1) |
| 847 |
2315 continue |
| 848 |
c d1i=d1(i1a,ipoi) |
| 849 |
c d1i=di(i1a) |
| 850 |
d1i=d(i1a,ipoi) |
| 851 |
do 2320,is=1,nnlm |
| 852 |
i=iln(is) |
| 853 |
ilmi=ilm(is) |
| 854 |
coeff(is)=coeff(is)+d1i*sigma1(is)* |
| 855 |
c (f1x(i)*g1x(ilmi)+f2x(i)*g2x(ilmi)) |
| 856 |
2320 continue |
| 857 |
2330 continue |
| 858 |
c s1=0.0D0 |
| 859 |
c s2=0.0D0 |
| 860 |
c s3=0.0D0 |
| 861 |
vari=0.0 |
| 862 |
do 2350,is=1,nnlm |
| 863 |
coeff(is)=coeff(is)*sigma1(is) |
| 864 |
c s1=s1+coeff(is)*split(is) |
| 865 |
c s2=s2+mnlm(is)*coeff(is) |
| 866 |
c s3=s3+coeff(is) |
| 867 |
vari=vari+(coeff(is)*sigma(is))**2 |
| 868 |
2350 continue |
| 869 |
c write (6,*) '#',j,k |
| 870 |
c write (6,*) s1,right(ipoi),s1-right(ipoi),s2,s3 |
| 871 |
c call set2davk(nnlm,isplit,msplit,coeff,avkern,nsr2,nt2) |
| 872 |
call set2davk(coeff,avkern,nsr2,nt2) |
| 873 |
if (iavkern.eq.2) then |
| 874 |
write (27) ((avkern(k1,j1),k1=1,nt2),j1=1,nr2) |
| 875 |
end if |
| 876 |
c call fitmn(x1,x2,nt2,nr2,avkern,yy2,mjta,4, |
| 877 |
c c real(tmesh1(k)),real(rmesh1(j))) |
| 878 |
write (24,*) mur,mut,real(rmesh1(j)),real(tmesh1(k)), |
| 879 |
c right(ipoi),sqrt(vari), |
| 880 |
c mjta(3),mjta(1),mjta(4),mjta(2) |
| 881 |
if (icoeff.eq.1) then |
| 882 |
write (28) (real(coeff(is)),is=1,nnlm) |
| 883 |
end if |
| 884 |
c end if |
| 885 |
c2370 continue |
| 886 |
2380 continue |
| 887 |
c close(27) |
| 888 |
c close(28) |
| 889 |
c close(21) |
| 890 |
write (6,*) 'averaging kernels set up, time=', |
| 891 |
c dtime(dummy),dummy(1),dummy(2) |
| 892 |
end if |
| 893 |
5000 continue |
| 894 |
if (iavkern.ne.0) then |
| 895 |
close (24) |
| 896 |
end if |
| 897 |
close (25) |
| 898 |
if (icovar.eq.26) then |
| 899 |
close (26) |
| 900 |
end if |
| 901 |
close (27) |
| 902 |
if (icoeff.eq.28) then |
| 903 |
close (28) |
| 904 |
end if |
| 905 |
open (20,file='status.log',form='formatted') |
| 906 |
write (20,*) 0 |
| 907 |
close (20) |
| 908 |
|
| 909 |
write (6,*) 'total time used:',etime(dummy) |
| 910 |
write (6,*) 'user time: ',dummy(1) |
| 911 |
write (6,*) 'system time: ',dummy(2) |
| 912 |
end |
| 913 |
|
