zonchn.f -- code to couple free oscillations for a zonally symmetric earth model

click title to see Google Drive directory JParkCodes

zonchn.f is a code to couple together free oscillations of Earth for a set of highly restricted aspherical earth models. The models are zonally symmetric, dependent only on co-latitude. The code has subroutines that integrate the modes against two types of models (1) a model that depends on a single wavenumber k, and (2) a model that is an equatorial band or a symmetric spherical cap. The first model type was used to examine the wavenumber behavior of surface-wave coupling, and the second type simulates the interaction of surface waves with a sharp boundary. The paper that describes the code's workings is

Park, J., The sensitivity of seismic free oscillations to upper mantle anisotropy I. Zonal symmetry, J. Geophys. Res., v98, 19933-19949, 1993.

Zonal symmetry in the model allows the coupling to vanish for all free-oscillation singlets with differing azimuthal numbers m, m'. In this manner, if there are two modes nSl and n'Sl' with a total of 2l'+1+2l+1 = 2(l+l'+1) singlets, the coupling problem splits apart into (for l'>l) 2l+1 2x2 coupling matrices, and l'-l uncoupled singlets. As a result of this symmetry, zonchn can readily accept a large number of modes in a narrow frequency band without breaking the computer memory or the CPU. The general scheme is to march up in frequency in steps of df and Df: All modes in a half-bandwidth Df about some fiducial frequency fo are coupled, but only those hybrid modes in a narrower bandwidth df about fo are saved. The code increases fo by 2*df and iterates up to some cutoff frequency. The user specifies how many overtones n are included. The uploaded version of zonchn.f has been tested with 0.2.le.f.le.20 mHz and df=0.1 mHz, Df=0.5 mHz and n-max=5.

The code depends on a set of subroutines in the file zonsubs.f, which are posted next in this blog.

The code reads from an input file "phzonchn" that contains

print *,'open(15,file=phzonchn,form=formatted)'
open(15,file='phzonchn',form='formatted')
call yread
print *,'output chain file '
read(15,102)namchn
print *,namchn
102 format(80a)
print *,'initial open of output file in append mode (0=no) '
read(15,*) iappend
print *,'iappend= ',iappend
print *,'enter inner freq halfwidth, outer freq halfwidth (mHz)'
c outer band is modes that couple, inner band contains modes that are saved
read(15,*) hfin,hfout
print *,hfin,hfout
print *,'enter freq band for hybrid modes (mHz)'
read(15,*) fb1,fb2
print *,fb1,fb2
fb1=dmax1(fb1,0.25d0)
print *,'enter max overtone branch to consider'
read(15,*) nbrmax

The option of writing the output file in "append" mode is an artifact of slow computers of the 1990s, in which one might need to stop and restart a long computation. I doubt that the user will need this option in 2009, though my test calculation required the greater part of an hour to complete.

example for phzonchn:

ymodel10_B0.05 --> aspherical model file
ymodel10_B0.05.zon20 --> output file of hybrid modes (linear comb of sph-earth modes, new frequencies)
0 --> dont open output file in append mode (usual case)
0.1 0.5 --> inner band contains modes that are saved, outer band is modes that couple (mHz)
0.2 20 --> frequency band to compute (mHz)
5 --> max number of overtones
/Users/jjpark/Mineos/Models/premiss --> earth model (isotropic PREM in this case)
/Users/jjpark/Mineos/sphnl --> N,L file of spheroidal mode parameters
/Users/jjpark/Mineos/sphe --> direct-access Fortran file of spheroidal modal eigenfunctions
/Users/jjpark/Mineos/tornl --> N,L file of toroidal mode parameters
/Users/jjpark/Mineos/tore --> direct-access Fortran file of toroidal modal eigenfunctions


The model file is based on the A,B,C,D,E Backus parameterization of hexagonally-symmetric anisotropy. A and D are isotropic P and S wavespeed perturbations (peak2peak variation: A=0.05 implies 5%). B is the 2\phi azimuthal Vp term, C is the 4\phi azimuthal Vp term, and E is the 2\phi azimuthal Vs term. The input, as explained in the post for zonsubs.f iare the values of the ABCDE and the wavenumber k of a corrugation in colatitude. The code does not allow the user to choose the equatorial band model on input. You must change a subroutine call and recompile the code. I may fix this soon, as it is cumbersome. Similarly, there is a slight alteration in the coupling formulas that switches the cases for an anisotropic symmetry axis w-hat aligned with theta-hat, or with phi-hat, that is, with colatitude and longitude, respectively. Due to symmetry restrictions, only w-hat parallel to theta-hat or phi-hat is allowed by this code. No tilted axes of symmetry, no general azimuthal angle.