program velpert c--this program reads mukappa model (Su & Dziewonski 97) directly c--from Weijia Su's files, converts it into perturbations in c--Vs and Vp (which requires reading prem). parameter( kmax=13) parameter( lmax=20) parameter(ndim=200) parameter( llmax=((lmax+2)*(lmax+1)/2-1)) c--set irem to 1 to remove the average from the model at any given depth parameter(irem=1) real mu real kappa real*8 R(ndim),VP(ndim),VS(ndim),Q(ndim),RF(ndim) real*8 qkappa(ndim),rho(ndim) dimension ac(0:llmax),bc(0:llmax) ! crust perturbation dimension a670(0:llmax),b670(0:llmax) ! 670 discontinuity dimension aw(0:llmax,0:kmax) ! model coefficients dimension bw(0:llmax,0:kmax) integer degmask(0:lmax),laymask(0:kmax+2) ! masks dimension ar(0:llmax),br(0:llmax),ylm(0:llmax) ! working arrays dimension cmph(0:lmax),smph(0:lmax) ! working arrays character*132 filename c--READ IN PREM call read_prem(ndim,r,vp,vs,q,rf,qkappa,rho) c--open output file and select depth at which we want to read the model. print*,'please input coordinates of the point' print*,'depth in km (less than 2900)?' read*,dep rad=6371.-dep print*,'longitude' read*,xxlo print*,'latitude' read*,xxla c--determine density, shear and P velocity at that depth according to prem. do ide=1,ndim if((rad.gt.r(ide)).and.(rad.lt.r(ide+1)))then dens=(rho(ide)+rho(ide+1))/2. alpha=(vp(ide)+vp(ide+1))/2. beta=(vs(ide)+vs(ide+1))/2. elseif(rad.eq.r(ide))then dens=rho(ide) alpha=vp(ide) beta=vs(ide) endif enddo c===================================================================== c===================================================================== c--READ MU c--(delta (mu over rho) )/(mu over rho) filename='mk12wm13_mu.dat' call RdHrvModel(degmask,laymask,ac,bc,a670,b670,aw,bw,kk,ll, . irftype,mantle,lmax,kmax,llmax,FILENAME,IREM) lon=xxlo klat=90.-xxla if((lon.gt.360).or.(klat.gt.180).or.(klat.lt.0))then print*,'error, wrong value of latitude' stop endif rr=rad c--use colatitude... theta=1.*klat phi=lon*1. theta=theta*3.14159265/180.0 phi=phi*3.14159265/180.0 c--evalulate at (rr,theta,phi) call CalVelPert(mantle,degmask,cmph,smph,ar,br,ylm,aw,bw,kk,ll, . irftype,theta,phi,rr,llmax,vel) c--The perturbation is in dv/v, you should mutiply it by c--100 to get perturbation in percent. mu=vel*100. c===================================================================== c--READ KAPPA c--(delta (kappa over rho) )/(kappa over rho) filename='mk12wm13_kappa.dat' call RdHrvModel(degmask,laymask,ac,bc,a670,b670,aw,bw,kk,ll, . irftype,mantle,lmax,kmax,llmax,FILENAME,IREM) lon=xxlo klat=90.-xxla rr=rad c--use colatitude... theta=1.*klat phi=lon*1. theta=theta*3.14159265/180.0 phi=phi*3.14159265/180.0 c--evalulate at (rr,theta,phi) call CalVelPert(mantle,degmask,cmph,smph,ar,br,ylm,aw,bw,kk,ll, . irftype,theta,phi,rr,llmax,vel) c--The perturbation is in dv/v, you should mutiply it by c--100 to get perturbation in percent. kappa=vel*100. c===================================================================== c===================================================================== c--determine P and S velocity from mu and kappa theta=1.*klat phi=lon*1. theta=theta*3.14159265/180.0 phi=phi*3.14159265/180.0 c--S velocity perturbation delvs_vs=.5*mu c--Bulk sound velocity perturbation delbu_bu=.5*kappa c--P velocity perturbation (see Su & Dziewonski (PEPI,1997)) sratiops=(beta*beta)/(alpha*alpha) delvp_vp=((4./3.)*sratiops)*delvs_vs delvp_vp=delvp_vp+((1.-((4./3.)*sratiops))*delbu_bu) c--write out the model(at that point) print*,'what follows is percent perturbation in mu over rho, & kappa over rho, bulk sound velocity, S velocity, P velocity, & respectively' write(*,*)mu,kappa,delbu_bu,delvs_vs,delvp_vp stop end c===================================================================== subroutine RdHrvModel(degmask,laymask,ac,bc,a670,b670,aw,bw,kk,ll, . irftype,mantle,lmax,kmax,llmax,FILENAME,IREM) integer degmask(0:lmax),laymask(0:kmax+2) dimension ac(0:llmax),bc(0:llmax),a670(0:llmax),b670(0:llmax) dimension aw(0:llmax,0:kmax),bw(0:llmax,0:kmax) character*132 filename c c initialize model coefficients c do i=0,llmax ac(i)=0.0 bc(i)=0.0 enddo do k=0,kmax do i=0,llmax aw(i,k)=0.0 bw(i,k)=0.0 enddo enddo c print *,"Harvard whole mantle model file:" c read(*,"(a)")filename io=1 open(io,file=filename,status="old") print *,"model file=",filename c c read the first line of the model file c call Read1stLine(io,degmask,laymask,irftype,kk,ll,lmax, &kmax,mantle,IREM) if(laymask(0).eq.1) ! read crust . call ReadLayer(io,ac,bc,ll) if(laymask(1).eq.1) ! read 670 . call ReadLayer(io,a670,b670,ll) do k=0,kk ! read model coefficients call ReadLayer(io,aw(0,k),bw(0,k),ll) enddo close(io) return end c===================================================================== subroutine Read1stLine(io,degmask,laymask, &irftype,kk,ll,lmax,kmax,mantle,IREM) integer degmask(0:lmax),laymask(0:kmax+2) character buff*132,ichar0*1,ichar1*1 include "velpert.h" read(io,"(a132)")buff read(buff,*)ll if(ll.gt.lmax)then print *, " angular order exceeds limit" print *, " l=", ll, " lmax=",lmax call exit(1) endif ichar0=' ' ii=0 do i=1,132 ichar1=buff(i:i) if(ichar0.eq.' '.and.ichar1.ne.' ')then ii = ii + 1 if (ii.eq.2)then read(buff(i:),"(132i1)")(degmask(j),j=0,ll) c--REMOVE THE AVERAGE? IF(IREM.EQ.1)DEGMASK(0)=0. else if(ii.eq.3)then read(buff(i:),*)kk else if(ii.eq.4)then read(buff(i:),"(132i1)")(laymask(j),j=0,kk-1) else if(ii.eq.5)then read(buff(i:),*)irftype goto 1 endif endif ichar0=ichar1 enddo write(*,"('invalid first line in the model file')") stop 1 kstr=kk kk=kk-2-1 write(*,"('spherical harmonic degree = ',i5)")ll write(*,"('radial order k = ',i5)")kk write(*,"('number of structures = ',i5)")kstr write(*,"('degree mask = ',132i1)")(degmask(j),j=0,ll) write(*,"('structure mask = ',132i1)")(laymask(j),j=0,kstr-1) mantle=WHOLEMANTLE ! fixed for WHOLEMANTLE model now if(irftype.eq.LEGENDRE)then print *," radial function = LEGENDRE" else if(irftype.eq.CHEBYSHEV)then print *," radial function = CHEBYSHEV" else if(irftype.eq.SHELL)then print *, " radial function = SHELL" else if(irftype.eq.BSPLINE)then print *, " radial function = BSPLINE" else print *, " unknown radial function=",irftype call exit(1) endif if(kk.gt.kmax)then print *, " order of radial function exceeds limit" print *, " k=", kk, " kmax=",kmax call exit(1) endif return end subroutine ReadLayer(io,ak,bk,ll) dimension ak(0:((ll+2)*(ll+1)/2-1)) dimension bk(0:((ll+2)*(ll+1)/2-1)) lm=0 do i=0,ll read(io,*)ak(lm),(ak(m+lm),bk(m+lm),m=1,i) bk(lm) = 0.0 lm=lm+i+1 enddo return end subroutine CalVelPert(mantle,degmask,cmph,smph,ar,br,ylm,aw,bw, . kk,ll,irftype,theta,phi,rr,llmax,vel) integer degmask(0:ll) dimension aw(0:llmax,0:kk),bw(0:llmax,0:kk) dimension ar(0:llmax),br(0:llmax),ylm(0:llmax) dimension cmph(0:ll),smph(0:ll) call SumRad(rr,ar,br,aw,bw,mantle,irftype,kk,ll,llmax) call CalYlm(ll,theta,ylm) call CSmPhi(phi,cmph,smph,ll) call SumSurCoeff(degmask,ar,br,ll,cmph,smph,ylm,vel,llmax) return end subroutine CalYlm(ll,theta,y) dimension y(0:((ll+2)*(ll+1)/2-1)) data y00,y10/0.28209479177388,0.4886025119029/ x=cos(theta) y0=y00 ! sqrt(1.0/(4*pi)) y(0)=y0 y(1)=y10*x ! sqrt(3.0/(4*pi))*x fac=2.0 sqrx2=sqrt(1.0-x*x) j=2 jd=2 do i=1,ll y0 = -sqrt(1.0+1.0/fac)*sqrx2*y0 y(j) = y0 fac = fac+2.0 if(i.ne.ll)y(j+i+1)= sqrt(fac+1.0)*x*y0 jd=jd+1 j=j+jd enddo tl=2.0 j=0 jd=0 do l=2,ll tl =tl+2.0 tl1 =tl-1.0 tlm1 =sqrt(tl1) tlp1 =sqrt(tl+1.0) tlm3 =sqrt(tl-3.0) do m=0,l-2 xlpm = (l-m)*(l+m) fac = tlp1/sqrt(xlpm) xlm = sqrt(xlpm-tl1)/tlm3 l2=j+m l1=l2+l-1 l0=l1+l y(l0)=(x*tlm1*y(l1)-xlm*y(l2))*fac enddo jd=jd+1 j=j+jd enddo return end subroutine CSmPhi(phi,cmph,smph,ll) dimension cmph(0:ll),smph(0:ll) cp=cos(phi) sp=sin(phi) cmphi=1.0 smphi=0.0 do m=0,ll cmph(m)=cmphi smph(m)=smphi tmphi=cmphi*cp-smphi*sp smphi=smphi*cp+cmphi*sp cmphi=tmphi enddo return end subroutine SumSurCoeff(degmask,az,bz,ll,cmphi,smphi,y,vel,llmax) integer degmask(0:ll) dimension az(0:llmax),bz(0:llmax) dimension y(0:llmax) dimension cmphi(0:ll),smphi(0:ll) vel=0.0 lm=0 do i=0,ll if(degmask(i).ne.0)vel=vel+az(lm)*y(lm) lm=lm+1 do m=1,i if(degmask(i).ne.0) . vel=vel+(az(lm)*cmphi(m)+bz(lm)*smphi(m))*y(lm) lm=lm+1 enddo enddo return end subroutine SumRad(r,ar,br,aw,bw,mantle,irftype,kk,ll,llmax) dimension aw(0:llmax,0:kk),bw(0:llmax,0:kk) dimension ar(0:llmax),br(0:llmax) dimension f(0:100) if(ll.gt.100)then print *," Error in SumRad, max radial order (kk) exceeds limit" print *," Please increase the size of f(0:100)" call exit(1) endif call RadFunc(r,f,mantle,irftype,kk) ! calculate radial function lm=0 do i=0,ll ar(lm)=0.0 br(lm)=0.0 lm=lm+1 do m=1,i ar(lm)=0.0 br(lm)=0.0 lm=lm+1 enddo enddo c c summerize model coefficients at radius r c do k=0,kk lm=0 do i=0,ll ar(lm)=ar(lm)+f(k)*aw(lm,k) lm=lm+1 do m=1,i ar(lm)=ar(lm)+f(k)*aw(lm,k) br(lm)=br(lm)+f(k)*bw(lm,k) lm=lm+1 enddo enddo enddo return end subroutine RadFunc(r,f,mantle,irftype,kk) dimension f(0:kk) include "velpert.h" if (mantle.eq.WHOLEMANTLE)then rtop=RMOHO rbot=RCMB else print*," not yet implemented for options other & than WHOLEMANTLE model" call exit(1) endif u = (r+r-rtop-rbot)/(rtop-rbot) if(irftype.eq.CHEBYSHEV)then call ChebyshevFun(u,kk,f) else print *," not yet implemented for radial & functions other than CHEBYSHEV" call exit(1) endif return end subroutine ChebyshevFun(x,kk,f) dimension f(0:kk) dimension cheby_coef(0:13) data cheby_coef/ . 0.70710678118655, 1.2247448713916, . 1.0350983390135, 1.0145993123918, . 1.00803225754840, 1.0050890913907, . 1.0035149493262, 1.0025740068320, . 1.00196657023780, 1.0015515913133, . 1.0012554932754, 1.0010368069141, . 1.00087070107920, 1.0007415648034/ if(kk.gt.13)then print *," error in function ChebyshevFun, kk>13" call exit(1) endif f(0) = 1.0 f(1) = x twox = 2.0*x do i=2,kk ! recursive relation f(i) = twox*f(i-1)-f(i-2) enddo do i=0,kk ! scale them f(i)=f(i)*cheby_coef(i) enddo return end c========================================================================== c========================================================================== subroutine read_prem(ndim,r,vp,vs,q,rf,qkappa,rho) IMPLICIT DOUBLE PRECISION(A-H,O-Z) DOUBLE PRECISION R,VP,VS,Q,RADIUS,RN,QKAPPA DIMENSION R(ndim),VP(ndim),VS(ndim),Q(ndim),RF(ndim),COEF(5), & QKAPPA(ndim) DIMENSION RHO(ndim) DFLOAT(I)=I iu=1 open(iu,file='foanis05.mod',status='old') READ(IU,1000,END=99) NREG,NIC,NOC,RN KNT=0 KNT1=0 c=============================== c=============================== c--loop over all the prem "regions" DO 10 NN=1,NREG READ(IU,1001) NLAY,IFS,R1,R2 DR=(R2-R1)/DFLOAT(NLAY-1) DO 2 I=1,NLAY RF(I)=R1+DR*DFLOAT(I-1) KNT=KNT+1 IF(KNT.EQ.1) RF(I)=.1D0 IF(NN.NE.NREG.AND.I.EQ.NLAY) RF(I)=RF(I)-.01D0 2 R(KNT)=RF(I) c=============================== c--loop over all the set of coefficients(physical obs: vs, vp, etc.) DO 1 J=1,8 READ(IU,1002) (COEF(I),I=1,5) IF(J.GT.5) GO TO 1 c=============================== c--loop over all the depths in that region DO 3 I=1,NLAY IND=I+KNT1 RAT=RF(I)/RN IF(IFS.GT.0.AND.J.LE.3) GO TO 40 VAL=COEF(1) DO 4 K=2,5 4 VAL=VAL+COEF(K)*(RAT**(K-1)) c--for that set of coefficients and that depth, calculate the model: c IF(J.EQ.1) RHO(I)=VAL IF(J.EQ.1) RHO(ind)=VAL IF(J.EQ.2) VP(IND)=VAL IF(J.EQ.3) VS(IND)=VAL IF(J.EQ.4) Q(IND)=VAL IF(J.EQ.5) QKAPPA(IND)=VAL GO TO 3 40 IF(J.GT.1) GO TO 41 Q(IND)=COEF(1) DO 42 K=2,5 42 Q(IND)=Q(IND)+COEF(K)*(RAT**(K-1)) GO TO 3 41 IF(J.EQ.3) GO TO 43 RHO0=COEF(4) c RHO(I)=RHO0*((1.D0-2.D0*Q(IND))**1.5D0) RHO(ind)=RHO0*((1.D0-2.D0*Q(IND))**1.5D0) 43 ST=Q(IND) E=1.D0-2.D0*ST c VAL=10.D0*DSQRT((E**2.5D0)*(COEF(1)+ST*(COEF(2)+.5D0*ST* c 1COEF(3)))/RHO(I)) VAL=10.D0*DSQRT((E**2.5D0)*(COEF(1)+ST*(COEF(2)+.5D0*ST* 1COEF(3)))/RHO(ind)) IF(J.EQ.2) VP(IND)=VAL IF(J.EQ.3) VS(IND)=VAL 3 CONTINUE c=============================== 1 CONTINUE c=============================== 10 KNT1=KNT1+NLAY c=============================== c=============================== NLEVEL=KNT c--knt is just the total number of layers RADIUS=RN NRAD=NLEVEL NLAYER=NRAD-1 NN=NLEVEL 50 IF(VS(NN).GT.0.D0) GO TO 51 NN=NN-1 GO TO 50 51 IF(NN.EQ.NLEVEL) RETURN NNP1=NN+1 DO 52 I=NNP1,NLEVEL VP(I)=VP(NN) 52 VS(I)=VS(NN) RETURN 99 STOP 1000 FORMAT(3I5,D15.8) 1001 FORMAT(2I5,2D15.8) 1002 FORMAT(5D16.9) END