PROGRAM DEFSOLVE C *************************************************************** C C solves for parameters which best reproduce observed deflection C C **************************************************************** implicit real*8 (a-h,o-z) parameter nt= 64 ! dimension of topo array parameter nk= nt/2 + 2 ! dimension of 1-d filter array parameter nfft= 150 + 6*nt ! dimension of FFT work arrays parameter ns= 9 ! number of Laplace domain samples parameter nsm= ns - 1 parameter nl=29 ! number of layers in Earth model parameter nm=1 ! number of parallel Maxwell elements parameter np= (nl)*(nm) + 6 ! number of parameters to solve parameter npp = np + 1 parameter npp2= np + 2 parameter nh=15 ! number of points in lake history parameter ntt=4 ! number of target times parameter nttm= ntt - 1 parameter mlk=4 ! number of lakes real*8 fzks(nk,ns),frks(nk,ns) real*8 kk(nk),tt(ns) real*8 a(nsm,nsm),b(nsm),f(nsm),s(nsm) real*8 fez(nk),bbz(nk,nsm) real*8 fer(nk),bbr(nk,nsm) real*8 kmax,kyr,norm real*8 fzet(nk),fzvt(nk) real*8 fz1d(nk),fz2d(nt,nt) real*8 fret(nk),frzt(nk) real*8 fr1d(nk),fr2d(nt,nt) real*8 z(nl),rho(nl),lam(nl) real*8 mu(nl,nm),eta(nl,nm) real*8 mul(nl) real*8 elev(nh),time(nh),ttime(ntt) real*8 htobs(nttm,400),htsig(nttm,400),dea(nttm,400),dno(nttm,400) real*8 htcomp(nttm,400) real*8 nomin,noref,eamin,earef real*8 z1(nttm,mlk),z2(nttm,mlk),dz1(nttm,mlk),dz2(nttm,mlk) real*8 thr(ntt),dht(nttm) real*8 p(np),pnom(np),delp(np),psig(np) real*8 resid(nttm,400),comp(nttm,400,npp) real*8 gtg(np,np),gtd(np),gg(np),gnorm(np) real*4 topo(nt,nt),ctopo(nt,nt) complex*16 load(nt,nt),ddflct(nt,nt) complex*16 dflct(nt,nt,nh),trdflct(nt,nt,nh) complex*16 ttdflct(nt,nt,ntt) complex*16 cwk(nt) real*8 rwk(nfft) integer*4 iwk(nfft) integer*4 nsh(nttm,mlk),nshtot(nttm),nlk(nttm) character*11 fname common /utmref/ nomin,noref,eamin,earef c weight= 10.0d0 c weight= 4.0d0 c weight= 1.0d0 weight= 0.25d0 zer= 0.0d0 one= 1.0d0 two= 2.0d0 pi= 3.141592653589793d0 cent= 3.1558d9 ! seconds in 100 years arrsz= 1024000.d0 ! topo array size ps= arrsz/dfloat(nt) rt2= dsqrt(two) nomin= 3973000. noref= 4485000. eamin= -182000. earef= 330000. kmax= two*pi/rt2/ps smax= 3.169d-9 ! 10 years smin= 3.169d-13 ! 10^5 years call getemod(nl,nm, z,rho,lam,mu,eta) call gettopo(nt, topo) call getshore(ps, nsh,nshtot,nlk,htobs,htsig,dea,dno) c begin perturbation loop do ip=1,npp2 call gethist(nh,ntt, time,elev,ttime) call perturbem(ip,np,nl,nm,z,rho,eta,nh,time,ntt,ttime, * p,pnom,delp) call fksgen(nk,kmax,ns,smin,smax,nl,nm,z,rho,lam,mu,mul,eta, * kk,tt,frks,fzks) call fks2t(nk,ns,nsm,kk,tt,fzks,a,b,f,s,fez,bbz) c call fks2t(nk,ns,nsm,kk,tt,frks,a,b,f,s,fer,bbr) c iterate deflection calculation at load input epochs itrmx= 4 ! number of iterations do itr=1,itrmx if(itr.ne.1) then call thresh(nh,nt,ntt,nttm,nsh,dea,dno,htobs,ps,dflct, * time,ttime,elev,thr) end if if(ip.gt.npp) then write(6,*) 'Lake Elevation history at iteration', itr write(6,100) elev 100 format(' ',8f10.2) c check deflection at center and threshold if(itr.eq.itrmx) then write(6,*) ' it time Zenda Stockton Lakeside' kyr= 3.1558d10 ish=1 do it=1,nh call interp(nt,nh,nttm,dno,dea,ps,dflct,ish, 1,it, defobs1) call interp(nt,nh,nttm,dno,dea,ps,dflct,ish,30,it, defobs2) call interp(nt,nh,nttm,dno,dea,ps,dflct,ish,49,it, defobs3) t= time(it)/kyr write(6,40) it,t,defobs1,defobs2,defobs3 40 format(' ',i4,f8.2,3f12.3) end do end if end if do it=2,nh-1 ! load specification times if(itr.eq.1) then call loadgen(nt,topo,elev(it), load) else call diffem(nt,topo,dflct(1,1,it), ctopo) call loadgen(nt,ctopo,elev(it), load) end if call fft2d(load,nt,nt,+1,nfft,iwk,rwk,cwk) t0= time(it-1) ! begin load t1= time(it) ! peak load t2= time(it+1) ! end load dt1= t1-t0 ! rise time dt2= t2-t1 ! fall time do jt=it,nh ! accumulate effects at (jt) due to load at (it) t= time(jt) dt0= t-t1 ! time past peak tps= dt0 + dt1 ! time past start if(tps.gt.zer) then call fktgen(nk,ns,nsm,tt,dt0,dt1,dt2,fez,bbz,fzet,fzvt) do ik=1,nk fz1d(ik)= fzet(ik) + fzvt(ik) end do call zfilt(nt,nk,fz1d, fz2d) call aplfilt(nt,fz2d,load, ddflct) call addem(nt,ddflct,trdflct(1,1,jt)) end if end do end do do it=2,nh call fft2d(trdflct(1,1,it),nt,nt,-1,nfft,iwk,rwk,cwk) call movem(nt,trdflct(1,1,it),dflct(1,1,it)) end do end do c compute deflections at target times do it=2,nh-1 call diffem(nt,topo,dflct(1,1,it), ctopo) call loadgen(nt,ctopo,elev(it), load) call fft2d(load,nt,nt,+1,nfft,iwk,rwk,cwk) t0= time(it-1) ! begin load t1= time(it) ! peak load t2= time(it+1) ! end load dt1= t1-t0 ! rise time dt2= t2-t1 ! fall time do jt=1,ntt t= ttime(jt) dt0= t-t1 ! time past peak tps= dt0 + dt1 ! time past start if(tps.gt.zer) then call fktgen(nk,ns,nsm,tt,dt0,dt1,dt2,fez,bbz,fzet,fzvt) do ik=1,nk fz1d(ik)= fzet(ik) + fzvt(ik) end do call zfilt(nt,nk,fz1d, fz2d) call aplfilt(nt,fz2d,load, ddflct) call addem(nt,ddflct,ttdflct(1,1,jt)) end if end do end do do itt=1,ntt call fft2d(ttdflct(1,1,itt),nt,nt,-1,nfft,iwk,rwk,cwk) end do call extract(nt,ntt,nttm,nshtot,dea,dno,ps,ttdflct, htcomp) call stats(ntt,nttm,nsh,nlk,dea,dno,htobs,htsig,htcomp, * z1,z2,dz1,dz2) if(ip.le.np) then write(6,*) write(6,*) '*********************************************' write(6,*) ' perturbed parameter number: ',ip write(6,*) else if(ip.eq.npp) then write(6,*) write(6,*) '*********************************************' write(6,*) ' perturbed parameter number: none' write(6,*) else write(6,*) write(6,*) '*********************************************' write(6,*) ' new solution ' write(6,*) end if call writem(nl,nm,z,rho,lam,mu,eta, nh,ntt,time,elev,ttime, * nttm,nlk,nsh,z1,z2,dz1,dz2,thr) if(ip.le.npp) then call savecomp(ip,npp,nttm,nshtot,htobs,htcomp,resid,comp) end if if(ip.eq.npp) then call solvem(np,npp,p,pnom,delp,psig,nttm,nsh,nlk, * resid,comp,htsig,gtg,gtd,gg,gnorm,weight) end if if(ip.eq.npp2) then call listem(ntt,nttm,nshtot,htobs,htcomp,htsig,thr,dht) call writemod(nl,nm,z,rho,lam,mu,eta) end if end do stop end