Column model pr

[ntlewis]

This pull request adds single-column model functionality to the Isca 'framework' for modelling planetary atmospheres. When the column model is used, the dynamical core is bypassed but Isca's 'physics' routines can still be called.

This functionality is useful for testing new convection / radiation parametrizations as the column model runs a gazillion times faster (especially if you're only simulating one column) than the full GCM. Due to the way the code has been rewritten, the user can in principle simulate many columns (in lat and lon) at the same time, however currently only using one core. Attempting to use >1 core will cause the model to crash.

This development should very much be considered to be in the BETA phase, as it has not yet been sufficiently tested, and because it may in time become clear that some things implemented by this pull request should be done in a different way.

The column model is currently initiated as a bit of a hack. The line

from isca import ColumnCodeBase

in a python experiment configuration script sets a compiler flag -DCOLUMN_MODEL that tells the model to use the following files:

atmos_column/column.F90
atmos_column/column_grid.F90
atmos_column/column_init_cond.F90
atmos_column/column_initialize_fields.F90

to initialize the model (including constructing the model grid), do the model timestepping
(using a leapfrog scheme as before), and handle input/output.

Some known issues at present:

1. As mentioned previously, Isca running in column model 'mode' can only use one core.
2. When configured as a single column, interpolator_mod will not interpolate input netcdf files as intended, but instead return an array full of zeros.
3. Spatially varying topography cannot be specified, however the user can impose a spatially homogeneous non-zero geopotential height if they really want to.

Issue 1) likely arises because the code is trying to make use of the spec_mpp routines, in spite of the fact that when the column model is configured a spectral space is not defined. This issue in principle could be resolved by implementing the mpp rountines from an FMS gridpoint GCM, however I have yet to investigate this possibility sufficiently.

Issue 2) can probably be resolved by creating an interpolator_1D interface, however even in this scenario it is not clear to me what the desired functionality here would be (e.g. averaging over an input file with lat and lon? or just selecting the input data that corresponds to the single column's location in space?).

A hacky work around for Issue 2 is implemented in rrtm_radiation.f90 to enable the user to supply a single vertical profile for ozone via the rrtm_radiation namelist. This code works in an identical way to that in vert_coordinate.F90 where the user is able to supply pressure levels via the namelist.

Example test cases column_test.py and column_test_rrtm.py are provided to show how the model can be configured. Additionally supplied in support of column_test_rrtm.py is scm_interp_routine.py which can be used to interpolate model input files (e.g. ozone) to the model pressure levels, and select a single vertical profile (or create a globally averaged profile of the standard model input) before running the model (i.e. preprocessing). The example configured by column_test_rrtm.py makes use of this functionality to create a globally averaged vertical profile for ozone which is supplied to rrtm via the namelist.

It is worth mentioning the following:

The wind is prescribed (it needs to be non-zero at the surface to allow for latent and
sensible surface heat fluxes). Currently the user can set a namelist variable 'surface_wind'
that sets u_surf and v_surf = surface_wind / sqrt(2), so that wind_surf = sqrt(u_surf^2 +
v_surf^2) = surface_wind. u and v at all other altitudes are set to zero (hardcoded).

At the moment the model needs to use the vertical turbulent diffusion parameterization in order
for the mixed layer code to work. This is not very consistent as the u and v wind are prescribed
and so the u,v tendenency from the diffusion is thrown away. Hence an implicit assumption when
using the column model is that 'the dynamics' would restore the surface winds to their prescribed
speed, so that du/dt total is zero. I have checked relatively extensively what using this parameterization is doing and it seems to be negligible / zero, so this is not really a concern. For example, when I tested the column by configuring it to solve for pure grey-gas radiative equilibirum, there was no appreciable difference between the result of the column model simulation and the analytical solution for grey gas radiative equilibrium.

[ntlewis]

I've added Brett's modifications to the column model. With these included, I'm happy for the SCM to go into the master. Just re-running the trip tests to make sure nothing is broken.

[ntlewis]

All trip tests pass. Nb: I had to modify the ape_aquaplanet test case python script to get it to compile and run with gfortran. There were a couple of namelist options set as floats when they should have been integers. When running the trip tests, I have compared this branch with my master (which includes two new commits to allow the socrates related tests to compile and run). @sit23 can you approve this so I can merge the SCM into the master?

[daw538]

I've made use of the column model in my private branches of Isca so I am familiar with the scheme and am satisfied that it works appropriately. There are a couple of additional changes compared to my version relating to diffusivity but I am happy with those.

At some point it would be great to have some docs added but I'm not going to hold back the acceptance of this into the main branch. Wondering if it might be useful to add to the header of column.F90 a text link to this PR just so users are aware of the 'health warnings' when using the model (which are listed at the top of the PR).

Otherwise, I am happy to approve this!