Gaea, Analysis and GFDL notes

/archive/Thomas.Flannaghan is where data is delivered after runs finish (automatic). This directory is on the analysis machine (also mounted as read-only on public). Data analysis should be done on the analysis machine. Log in via public using a different RSA code.

Data is delivered to rather cryptic folders. I've written a little tool called archtool.py to automate the organising of runs. This tool displays a list of all runs on archive along with when they were created. At present, the only thing that can be done is automatically extracting the data (and run namelists) to a new sub-directory ("run code" in the run index). The tool automatically recognises and handles restart runs.

• gcp allows movement between gaea, analysis and GFDL but not outside it.
• Normal scp and ssh are available on GFDL workstations (i.e. public).
• archive is mounted on public at /archive/$USER/ (read only) so we can do public> scp /archive/$USER/.... atmos02:/... to move data to sayre.

All settings for compiling experiments and for running jobs can be specified in the XML file. frerun remakes all namelist files so if only runtime options in the XML are modified, we do not need to remake.

The "group letter" is found in the XML inside the project tag. This is "gfdl_w". Any references to c1 and c2 should be removed from this string. The job will silently fail if this is wrong or not set :-(. An XML example is

...
<platform name="ncrc2.default">
  <directory stem="$(FRE_STEM)"/>
  <project>gfdl_w</project>
  ...

An experiment is defined as a unique codebase. We do not need new experiments if we are just changing namelists. They are defined in the XML.

1. Copy an existing experiment in the XML
2. Change the name and path of the experiment (see example below)
3. Run fremake but don't run command given at end of fremake
4. Check to see that CVS operations took place during fremake. If not, the name or path was not changed correctly (see step 2) and an existing experiment was in the path given.

The experiment XML should end up like this:

<experiment name="{experiment name}">
  <description>
    Blah blah blah
  </description>
  <component name="fms_spectral_solo" 
             paths="$(rootDir)/{experiment name}/src" 
             includeDir="$(srcDir)/shared/include">
    <source versionControl="cvs">
      ...

where {experiment name} is the name of the experiment (must be unique in the XML file).

There is an interactive mode which is good for debugging compiler and runtime errors. To get an interactive session, use

msub -I -l walltime=01:00:00,size=16

This gives a terminal that will work for 1 hour with 16 cores (the number of cores must be the same as npes in the XML – see *Horizontal resolution). From this interactive terminal, we can make the model using make (with env script) or using the build script, and we can run the model executable using

aprun -n 16 fms_HSt42_.....x

Running the model produces a lot of netcdf files which need combining and sending to archive. I'm not sure how this is done yet. Again, the n option refers to the number of cores and must be the same as npes in the XML.

All runs (regression and production) can be restarted.

The HSt42 model defined in idealized.xml gives a typical T42 resolution run, but we can modify the xml to give different resolutions. The key spectral_dynamics_nml namelist parameters are

namelist parameter	T21	T42	T85
num_fourier	21	42	85
num_spectral	22	43	86
lat_max	32	64	128
lon_max	64	128	256

In addition, we must change the regression or production run specification to indicate the number of processors. For regression runs, the specification for T42 is something like

<regression name="...">
  <run days="..." npes="16" runTimePerJob="..."/>
</regression>

The npes field is something to do with the number of processors required for the run (this is also the number used when triggering an interactive session). For T85, increase npes to "32". When the run actually is submitted, it typically uses twice the number of cores as specified here.

The spectral_dynamics_nml controls the vertical coordinate. The vert_coord_option allows various different ways to define vertical coordinates. The simplest of these are "even_sigma" and "uneven_sigma".

"even_sigma" simply divides the vertical coordinate evenly in sigma, from 0 to 1. num_levels gives the number of levels.

"uneven_sigma" allows more complicated sigma coordinates. Again, num_levels gives the number of levels. The configuration namelist settings are as follows

option	meaning	Default
surf_res	Confusing – see exact definitions	0.1
exponent	Confusing – see exact definitions	2.5
scale_heights	Top level expressed in scale heights	4
zero_top	Forces top sigma level to be 0	.true.

Reasonable (i.e. approx 1 km) resolution at the TTL/lower stratosphere can be obtained with 60 levels, and scale_heights = 10 (the remaining parameters are set to the defaults). Polvani-Kushner use a set of custom specified levels but these lie closest to

=scale_heights=       =11.0,
=surf_res=            = 0.5,
exponent            = 3,

These values are from Martin Jucker (he fitted these parameters to the PK levels).

The exact definitions (from atmos_spectral/init/vert_coordinate.F90):

p = a * p_ref + b * p_surf

where p_surf is the instantaneous surface pressure and p_ref is some constant. uneven_sigma defines a and b as follows

! a is set to zero
a = 0

! b is defined by:
do k = 1,num_levels
   zeta = 1 - float(k - 1) / float(num_levels)
   z    = =surf_res= * zeta + (1 - surf_res) * zeta**exponent
   b(k) = exp(-z * scale_heights)
enddo

! sort out the model surface level and top level
b(num_levels+1) = 1.0
if (zero_top) b(1) = 0.0

Therefore, when surf_res = 1 we have roughly evenly spaced levels in height (exactly even in log-pressure coordinates). When surf_res = 0 we have exponentially spaced levels. This means that surf_res doesn't really have anything to do with the surface – it is more like saying "how close to linear are we?". It's value affects the resolution everywhere.

Here is a simple calculator which also gives the approximate height (log-pressure height).

# namelist parameters:
num_levels = 10
scale_heights = 11
surf_res = 0.5
exponent = 3

import math
for k in range(1, num_levels + 1):
    zeta = 1 - (k - 1) / float(num_levels)
    z = surf_res * zeta + (1 - surf_res) * zeta ** exponent
    b = math.exp(-z * scale_heights)
    b_km = 7 * z * scale_heights       
    print 'b(%d) =\t%0.3f (approx %0.1f km)' % (k, b, b_km)

There is no sponge by default in the HSt42 model defined in idealized.xml. There is only an optional very simple sponge available, which acts on only the top level (defined in spectral_damping.F90). It is defined by the following parameters in spectral_dynamics_nml:

option	meaning	default
eddy_sponge_coeff	Eddy damping in top layer	0
zmu_sponge_coeff	Zonal mean u damping	0
zmv_sponge_coeff	Zonal mean v damping	0

All coefficients are in units of m²/s, and are applied as damping terms like \(coeff * \Delta^2\).

The Polvani-Kushner hs_forcing.F90 file (given to me by Martin) has a more advanced sponge. I have added it to my hs_forcing_zonalacc.F90 file. The following options control the sponge, and are given in the hs_forcing_nml namelist.

option	meaning	default
sponge_flag	Is sponge used?	.false.
sponge_tau_days	The sponge timescale (days)	1
sponge_pbottom	The lowest level to apply sponge (Pascals)	100

The sponge is a Rayleigh friction with coefficient given by sponge_coeff * (1 - p / sponge_pbottom) ** 2. sponge_coeff is the inverse timescale derived from sponge_tau_days.

option		default
t_zero	Surface temperature at equator	315
t_strat	Stratospheric temperature	200
delh	Equator to pole temperature difference	60
delv	Delta Theta per scale height	10
eps	Hemispheric assymetry	0
sigma_b	Boundary layer sigma	0.7
ka	Newtonian damping in free atmos	-40
ks	Newtonian damping in boundary layer	-4
kf	Rayleigh damping in boundary layer	-1
do_conserve_energy	Converves energy dissipated via	.true.
	Rayleigh damping (and other momentum
	forcings) by heating.

The damping rates k are given in units of days if negative. Not sure about units if positive, but probably days^-1. Temperatures are all in Kelvin.

do_conserve_energy does not really affect the solution to the local_forcing-type problems. Both have been tested with minimal differences (<5% difference).

This is my custom code documentation. I have added this behaviour to hs_forcing_zonalacc.F90, which replaces hs_forcing.F90 (using a symbolic link). The namelist parameters go into the hs_forcing namelist. All options are prefixed with local_forcing and are given here:

local_forcing*	meaning	units	default
_option	Defines the forcing mode		''
_xwidth	Longitudinal half-width	degrees	10
_ywidth	Latitudinal half-width	degrees	10
_ycenter	Central latitude	degrees	0
_zwidth	Half-width in height	scale height	0.5
_zcenter	Central height.	scale height	2.2
_amp	Amplitude of forcing.	m/s/day	1
_start	Start time of forcing.	days	-1
_end	End time of forcing.	days	-1

local_forcing_start and local_forcing_end allow control of when the forcing is turned on and off. A value of -1 (the default) removes that time constraint. For example, switching on from 100 days would give local_forcing_start = 100 and local_forcing_end = -1.

local_forcing_option is a string and can have several values:

local_forcing_option	meaning
''	No forcing. The local_forcing code is not run.
'Cos2Dipole'	A two-signed forcing in the vertical.
	For positive amplitude, negative forcing is
	on the top and positive forcing is on the
	bottom. The forcing has the form
	cos(x)2 cos(y)2 sin(z)
'CosDipole'	Vertical structure as above, but with
	cos(x) cos(y) horizontal structure.
'Cos2Single'	Single signed forcing (cos(z) vertical)
'CosSingle'	– ditto –
'CosDipoleZM'	Zonal mean CosDipole – Zonal mean caveats
'Cos2DipoleZM'	Zonal mean Cos2Dipole – Zonal mean caveats

from numpy import *
Lx = 30 # degrees
Ly = 10 # degrees
A = 0.3 # zonal mean amplitude

x = linspace(-180, 180, 5000)
y = linspace(-10, 10, 100)

f = cos(0.5 * pi * x / Lx) * cos(0.5 * pi * y[:,None] / Ly) \
    * ((abs(x) < Lx) & (abs(y[:,None]) < Ly))
f *= A / f.mean(axis=-1).max()
print 'zonal mean amplitude, A     ', f.mean(axis=-1).max()
print '10N-10S zonal mean amplitude', f.mean()
print 'local amplitude             ', f.max()
print '10N-10S local amplitude     ', f.mean(axis=0).max()
print 'Theoretical local amplitude ', A * pi * 360 / (4. * Lx)  

The base version of hs_forcing.F90 includes Isidoro's local heating. In my version, I've removed this and replaced it my a local heating routine that mirrors the local forcing routine. All namelist variables do the same thing here as they do for local forcing and have the same defaults, but are named local_heating_xxxx in place of local_forcing_xxxx.

The only difference is that local_heating_amp is in units of K/day.

See above for all local_forcing settings.