*******************************************************************************
*         			README				              *
*******************************************************************************

1. General Description
2. Files and Usage
2.1. Program Files
2.2. Input and Output
2.3. Options for Different Usage
2.4. Additional Notes
3. Test Run
3.1. Incompressible Run
3.2. Compressible Run

*******************************************************************************
1. General Description
The code is designed to solve for the love numbers for compressible Earth, and
to obtain the present day geoid/uplift rate by combining the love numbers with
an ice load model. For an incompressible Earth model, the code works as long as
the p-wave velocity is set to be large enough (check 3.2 for details).

2. Files and Usage
2.1 Program Files
The code contains 2 program files:
- test.f       : solve for the love numbers including deg-1 and deg-(2,1) terms
- test.stokes.f: combine the love numbers with ice history

test.f can be run separately, while test.stokes.f needs the
input from test.f.

<test.f> generates the love numbers that can be used to obtain the present day
geoid/uplift rate. For most cases, this routine needs to be run only once, and
its output contains elastic love numbers and viscous residues for all the
degrees.

<test.stokes.f> uses the love numbers to find the geoid/uplift rate. This can be
done with/without sea level equation iteration, with/without the polar wander
effect.


2.2 Input and Output

Input files for test.f:

- Main input file: in.dahlen.viscoelastic 
First line: the fluid love number (the user needs to specify the fluid love number 
            for the particular Earth model they used) 
Second line: the spherical harmonic degrees to consider (minimum degree, maxmum 
             degree, interval) 
Third line: the collocation parameters (minimum of period, maxmum of period, interval). 

#Note: the collocation parameters for the lower degree terms (e.g., degree 1 to 
degrees 11) can also be tuned in the code (see the line 363 to 368 in the code)

- Earth model: model.vm5a or (model.prem+viscosity.model)
First line: name of the model, lithosphere thickness, the number of the data lines
Data section: contains radius of different layers (first column), density profile 
              (second column), P-wave velocity (third column), S-wave velocity 
              (forth column), gravity (fifth column) and viscosity (sixth column).

#Note: The discontinuity boundaries are denoted by repeated radius (e.g., line 161 and 162)

The last second line: the averaged density
The last line: the total number of the discontinuity boundaries followed by the line
                numbers for each boundary.

#Note: If model.prem is used (instead of vm5a), viscosity.model must be included.


Output for test.f
- screen output: surface g value; for each degree, the result of Mark's test is
printed out. If numbers close to 1 are printed out in the first column, the code
is running all right. See more details about Mark's test method in A et al. (2013).

- Output files:  out.dahlen.love, out.dahlen.coll, out.btide, out.rot
These files contain all the love number information(including rotation love numbers), 
both elastic and viscous. Those love numbers are used as input files for in.stokes.ice5 
for the convolution calculation.

Input for test.stokes.f:

- Main input file: in.stokes.ice5  
The first entry of in.stokes.ice5 is always set to 1 here.
Second line: the switch for sea level equation (0 means off; 1 means on), the number of iteration
             for the sea level equation calculation, the switch for the polar wander 
             effects (0 means off; 1 means on))  
Third line: the file name for the ice models (+ static ocean load)
Forth line: the maximum number of degrees in ice models to be included in the convolution
Fifth line: the file name for the output
Sixth line: the fluid love number (should be the same number with that in in.dahlen.viscoelastic)


- Ice models(+static ocean load): ice6g+ocn.lmax100.txt
The spherical harmonics expansion of the ice model + static ocean model for each loading epoch.
Please see the definition of the static ocean variations in A et al. (2013).

- Ice loading epochs: years.reverse (in ka)

#Note: The total number of the loading epochs needs to be specified in the code (e.g., in test.stokes.f search 
"nepochs=" and put the total number of the loading epochs)


- Ocean functions: 
ocn1.0_xyt.txt: land/ocean mask (1 means ocean and 0 means land) in spatial domains (lat/lon)
ocn1.0_lm.txt: land/ocean mask in spectrul domains (spherical harmonics expansion)

Output for test.stokes.f:
- stokes.ice6g: the name can be changed in <in.stokes.ice5>. It contains
the present day geoid (the 3rd and 4th column - Clm and Slm), vertical (the 5th and 6th column) and horizontal(the 7th and 8th column) displacement rate. Those numbers are Stokes coefficients (1/yr).


2.3 Options for Different Usage
This (in)compressible code can run with/without polar wander
effects, with/without sea level equation iterations. Users can always generate
different results by changing options in <in.stokes.ice5>


#Note: if you use the same Earth model for different loading history, no change is needed for test.f, 
it generates all the love numbers from degree 1 to degree 100, or higher if you wish, and can be used 
for different loading history.

2.4 Additional Notes
1.The g (Earth surface gravity) value is computed in <test.f>. But in <test.stokes.f>, they need to be 
specified. g values need to be changed manually in <test.stokes.f>. And g is used when the code 
generates stokes' coefficients. Change has to be made at TWO places. Search for 'change g', and change 
the g value according to either the screen output from <test.f> or any convention you use for stokes' 
coefficients.

2.Change the collocation parameters with care. The current one in the input file
should work for vm5a and PREM Earth model. If a different model
is used, and the screen output of <test.f> might produce numbers that are
significantly different from 1 in the 1st column. In this case, the collocation parameters
 need to be adjusted.

3.This code can also output the surface geoid, vertical displacement and horizontal displacement at each
loading epochs. To do this, user need to modify the code as following:
(1) run test.stokes.f with the setup in the previous section
(2) go to test.stokes.f, uncoment the line 60 and comment out the line 61
(3) run test.stokes.f again to output the results for each epochs
(4) check the results (e.g., PW.ice6g.***, the same format with that for stokes.ice6g)

3. Test Run
3.1 Incompressible Run
The code is initially set up for an incompressible run, using vm5a model, with ICE6G ice loading
history (e.g., Peltier el al. [2015]).

Open <in.stokes.ice5>, change the 2nd line to 1 3 1
Compile <test.f>
Run the exe file
Compile <test.stokes.f>
Run the exe file

#Note: For the incompressible case, the p-wave velocity in the Earth model is set to be large enough.
(e.g., in the line 89 of test.f, the vp = vp*1.e3)
 
#Note: You may find NaN in the screen output from <test.f>, if vp is too large.

3.2 Compressible Run
Users can always run a compressible model as long as the p-wave velocity in
the Earth model is taken as some realistic values (e.g., take the PREM model by commenting out line
86 in test.f).
Compile and run <test.f>. The compressible love numbers will be produced.
 
Reference:
A, G.R., Wahr, J. & Zhong, S.J., 2013. Computations of the viscoelastic response of a 3-d compressible earth to surface loading: an application to glacial isostatic adjustment in Antarctica and Canada, Geophys. J. Int., 192, 557-572.

Peltier, W.R., Argus, D.F. & Drummond, R., 2015. Space geodesy constrains ice age terminal deglaciation: The global ICE-6G_C (VM5a) model [dataset], J. geophys. Res., 120, 450-487.



*******************************************************************************
*                                     END                                     *
*******************************************************************************









