zonsubs.f -- the subroutines needed to make zonchn.f work

click title to see Google Drive directory JParkCodes

This set of subroutines serves zonchn.f, see the earlier post. The first routine yread reads the parameters of the zonally symmetric model example file here

1 1 ! useless parameters for option I never implemented
0.0 0.05 0.0 0.0 0.0 A,B,C,D,E relative perts to the Backus anisotropy parameters for azim anisotropy with an axis of symmetry.
1 152 181 ! lithosphere in model premiss
10 ! k=10 wavenumber on the sphere

Some relevant passages of code are the following:

c angint is original theta integral for cos{(l+1/2)\theta} structure
c angint2 uses 3-j tricks to speed up the calculation
c angint2X uses 3-j tricks, but for either spherical cap or equatorial band
c angint2XX is s=20 expansion of equatorial band model, no discontinuites
c call angint2(l1,l2,ls(kii),angi)
call angint2X(l1,l2,ls(kii),angi)
c call angint2XX(l1,l2,ls(kii),angi)


c combine isotropic and anisotropic parts, use b,c,e & a,d
c **************************** note that code can incorporate only
c ***************************** one aspherical wavenumber, despite
c ***************************** the do loops
c NOTE: some lines have alternate versions for phi and theta as axis of sym
if(ick.eq.3) then ! tor-tor
do j=1,nom
coup(j)=coup(j)
x +angi(3,j)*(c*ysc(1,1)+e*ysc(2,1)) ! anisotropic
c x -angi(5,j)*e*ysc(3,1)-angi(6,j)*c*ysc(1,1) ! anisotropic-phi
x +angi(5,j)*e*ysc(3,1)-angi(6,j)*c*ysc(1,1) ! anisotropic-theta
x +angi(2,j)*d*ysc(3,1)+angi(3,j)*d*ysc(2,1) ! isotropic
end do

comment or uncomment the phi and theta terms as you wish, but dont comment both or uncomment both.