This is a user story for how to create / adapt the GFDL CM2.1 model for novel topographies, for deep time paleoclimate runs. I will write this story in stages… currently a work in progress.
0. Supporting files
See the example repository GitHub - dkhutch/cm2.1_paleo: Sharing tools and input files for creating new boundary conditions using GFDL CM2.1, for paleoclimate simulations for a working example of the inputs needed to create a new land-sea mask. This repo has 3 folders:
a55_c1_grids: A set of input grids for a 55 Ma simulation
tools: A set of executables (compiled on Gadi) for generating new grid
bin: the CM2.1 model executable and collation tool.
See also the input file directory for the Eocene 55 Ma, 1xPI CO2 run:
1. Adapting the land-sea mask
First, look at the a55_c1_grids folder. This folder contains grid files for the ocean, atmosphere and land:
atmos_hgrid.nc
land_hgrid.nc
ocean_hgrid_1x1.5_82S.nc
Note: each of these grids is a ‘supergrid’, containing both the cell centres and cell vertices, i.e. (2n+1) as many x and y points in each direction as the model itself.
To generate a new land-sea mask, you first need to create a new ocean topography file. The working example given is called: topog.nc
The model uses a tripolar grid, with a transition from regular lat-lon south of 65N, and a bipolar grid north of 65N that has poles over Siberia and North America. To interpolate a new topography onto the model grid, use make_topog from the tools folder, using a command like this:
./make_topog --mosaic ocean_mosaic.nc --vgrid ocean_vgrid.nc --topog_type realistic --topog_file input_topog.nc --topog_field topo --scale_factor -1 --deepen_shallow
Here, the tool takes an input topography file input_topog.nc, which has topography variable topo, interpolates it to the model grid, and applies a scale factor of -1 to account for depth being positive downwards. NOTE: input_topog.nc can be any topography input file you like (on a regular lat-lon grid), such as a modern elevation dataset or a paleogeography of your choice.
The output file will be
topog.nc
like in the example repository. You can of course edit topog.nc using your favourite netcdf editing tool, just bear in mind that the Arctic region is complicated by the tripolar grid.
To generate a land-sea mask for this new topography, you then run the command make_coupler_mosaic from the tools folder, like this:
./make_coupler_mosaic --atmos_mosaic atmos_mosaic.nc --ocean_mosaic ocean_mosaic.nc --ocean_topog topog.nc --land_mosaic land_mosaic.nc
Here, the input grids for each model component are specified via the mosaic file such as ocean_mosaic.nc, which points to the input grid and specifies that is has periodic longitudinal boundary conditions. Once run successfully, make_coupler_mosaic will output the exchange grids:
atmos_mosaic_tile1Xland_mosaic_tile1.nc land_mosaic_tile1Xocean_mosaic_tile1.nc
atmos_mosaic_tile1Xocean_mosaic_tile1.nc
and the land-sea mask files:
land_mask.nc ocean_mask.nc
It will also create a “master” mosaic file called mosaic.nc. Before using in the model, you need to rename this to grid_spec.nc
mv mosaic.nc grid_spec.nc
The model will always look for grid_spec.nc for its input grid information.
Now you have a new land-sea mask and exchange grid! The next step is to create a new atmosphere topography.
2. Create a new atmosphere topography
Here we want to recreate an input topography for an arbitrary paleogeography. In the input directory:
There is a restart file called
fv_rst.res.nc
This restart file contains a variable called Surface_geopotential, with dimensions (Time, lat, lon), where Time is a single timestep. The surface geopotential is the topography multiplied by gravitational acceleration grav, which in our model has the exact value:
grav = 9.8
To create a new version of the topography, simply interpolate your desired paleogeography onto the model grid (in this case, 60 x 96 grid cells) and multiply by grav. Also, all negative values (i.e. below sea level) should be set to zero, since they are seen by the atmosphere as ocean. You may also want to consider enforcing that any points which have a land fraction of 0 (from land_mask.nc in Step 1), should also be set to zero, after your interpolation step. However, this is not a strict requirement of the model, just helps for physical consistency.
All other variables in fv_rst.res.nc can be left unchanged, i.e. you can use the same input file in the example but just overwrite the Surface_geopotential.
3 Create input file for mountain drag parameterisation
In the run directory for the a55_c1 simulation,
The main input file input.nml contains a namelist called topography_nml, which specifies an input file topography file:
topog_file = 'INPUT/atmos_topog_Herold.nc'
WARNING: THIS IS NOT THE TOPOGRAPHY INPUT. This is actually the input file for the mountain drag parameterisation. This input file is usually specified at higher resolution than the model resolution, i.e. for the pre-industrial simulation it is specified at 0.16 degree resolution. Ideally, a high resolution allows the mountain drag scheme to more accurately calculate topography gradients, which determines the mountain drag scheme.
For the Eocene simulations given in our example, we only had a 1 degree resolution, hence we used that as our input to the mountain drag scheme. Bear in mind, the latitude and longitude for this input file are specified in radians, and are labelled xdat and ydat, varying from 0 to 2*pi in longitude, and from -pi/2 to pi/2 in latitude. Also, the latitude and longitude values are specified at cell edges, with (n+1) points in each direction. Follow the example of atmos_topog_Herold.nc for a guide on how to structure the variables. The topography itself is simply called zdat (elevation) and all values below sea level have been set to zero.
If you have a higher resolution paleogeography, simply use the high-resolution version as an input, and adapt the field names to be equivalent to xdat, ydat and zdat as in the example.
4 Create new river runoff
The river runoff is specified through the input file river_destination_field. This is an ASCII input file that contains the following:
1st line: (number of x points) (number of y points) (starting longitude),
In the example this is: 96 60 0
2nd line: a fortran format statement specifying how the data should be read in. In our case this is (20i4), i.e. 20 integers occupying 4 columns of text
Remaining lines: A list of indices specifying the x and y indices of river runoff for each grid cell. This is given as an ordered list of x-indices, followed by an ordered list of y-indices.
Since this input file is difficult to comprehend, it is best to retrace the method used to create the indices, which will enable you to visualise the river basins themselves.
The code and examples for this section are found in this folder:
4.1 Run a downslope relocation algorithm from a new topography
In the python script downslope_coarse.py I have written a function to take an input topography (using the input file Herold_drainage.nc), and use a downslope relocation algorithm to generate a river basin map. For each land grid cell, the algorithm searches for the lowest grid cell in the immediate neighbours (firstly in the N, S, E and W directions), and then in the diagonal neighbours if none are found (NW, NE, SE, SW). The algorithm then incrementally steps to the lower point, and searches again for a lower neighbour, until either (a) it reaches the coast, or (b) it reaches a maxcount threshold.
I set the maxcount threshold at 100 steps - suitable for a coarse resolution grid - to avoid getting stuck in an infinite loop. When the coast is found, the algorithm stops and saves the x and y indices for the destination of the river flow (named idest and jdest). If the algorithm gets stuck in a land sink, the maxcount is reached, and the location is saved in the land_sinks variable.
The idest and jdest variables are then used to construct a basin map. Land points are allocated a positive integer, where each positive integer is a unique river basin, river destination points at the coast are allocated a corresponding negative integer where the river runoff ends up. I.e. A river basin with points labelled +n will flow to the destination -n. Ocean points are labelled 0.
This output is then saved in a text output file, in this case labelled Herold_drainage_basins.txt. If you view this file in a text editor with no word wrapping, you can visualise the river basins from the integers themselves. You can also view a plot of the integers in Herold_drainage_basins.png. This illustrates how the basins with positive numbers flow to the negative numbers on the coast.
4.2 Eliminate any land-based sinks (lakes)
When you apply this algorithm to a new topography, you will inevitably generate some land sink (i.e. lakes). The model cannot handle lakes, so you have to get rid of them. The downslope_coarse.py creates an output file called land_sinks.nc which displays where the algorithm got stuck. It also creates an output plot called counts.png. This helps to visualise where the land sinks are, as well as providing a graphical representation of the idest and jdest variables.
When you find a land sink, the simplest solution is to edit the topography to make it go away. An example of how I have done this is in fix_drainage.py. In that file, make a number of manual adjustments to the topography, to provide a pathway for river runoff to ‘break out’ of a land sink. Those choices were arbitrary and guided by eyeballing the topography and choosing a path of least disturbance. These adjustments are then fed into regenerating the Herold_drainage.nc file, to make an updated river runoff map (Note, the Herold_drainage.nc topography is only used to make river runoff, it is not used as an actual topography.)
Making manual adjustments to fix the drainage is admittedly a messy solution, but the output file land_sinks.nc at least guides the user as to where adjustments need to be made. More sophisticated algorithms to enforce drainage exist (e.g. ANUDEM ) but I have taken a bare bones approach.
4.3 Convert the basin map to a river_destination_field
Once you have a basin map (in our case Herold_drainage_basins.txt), you can then convert it to the river_destination_field input that is required for the model. We use a fortran program called make_river_dest.f90, in the folder given above. This takes the Herold_drainage_basins.txt input and converts it to the fortran formatted idest, jdest indices that are read into the model. Here we use fortran indices (using base 1). The output file river_dest_Herold_drainage is the final output. I would then rename this as river_destination_field when feeding back into the input files for the model.
river_destination_field itself is difficult to comprehend for a human being. For this reason, I like to run a double-checking script to re-calculate the river basins from the final output.
4.4 Double check the basin map from the river_destination_field
As a sanity check, you can run the script reverse_rivers.py, which takes the final output (in our case river_dest_Herold_drainage) and calculates a basin map, saved in basins_reversed.txt, with graphical output in basins.png and netcdf output basins.nc. This allows you to inspect an exisiting river_destination_field input file that has been previously made (by anyone) to see how the river basins look in an arbitrary configuration.
5 Vegetation and land surface boundary conditions
Files for section 5 are provided in this Dropbox folder:
5.1 Vegetation file: cover_type_field
Here we convert the Herold boundary condition file: biome_1x1.nc, which presents Eocene vegetation types using the BIOME4 plant types. The first step is to convert the BIOME types to the CM2.1 model types, which we do by running the script:
python biome_to_veg.py
This creates the file cover_eo.nc, which contains vegetation types in the CM2.1 format. The types are converted using the function biome_convert() within that script. Each integer represents a particular plant type. There is an explanation within biome_to_veg.py of what the different plant types mean in the CM2.1 model.
The next step is to run:
python make_eo_cover.py
This converts cover_eo.nc to Herold_cover_type_field. This file is then renamed to cover_type_field and used as an input file, i.e.
cp Herold_cover_type_field cover_type_field
Note: we are leaving the vegetation input file on the 1x1 deg grid as in the Herold boundary conditions. The CM2.1 model interpolates this automatically on the land grid, i.e. we don’t need to match it exactly to our model grid.
Last step is to run the following script:
make_cover_field_nc.py
This is mainly to convert an existing cover_type_field into a netcdf file so it can be easily visualised, using e.g. ncview. Running that script will create a file called Herold_cover_type_field.nc, which has the same data as the cover_type_field. But, I also use the Herold_cover_type_field.nc as an input for creating the soil type (below).
5.2 Soil type: ground_type_field
In this step, I set the soil type to be uniform for all land points everywhere. This is a bit of a quick and dirty solution, as one could ideally set the soil type to be “congruous” with the vegetation type. However, it does make it easy to modify to new land configurations. Simply run:
python uniform_land.py
This takes the land points as found in the Herold_cover_type_field.nc and converts them all to soil type 2, while the ocean points are set to 0. As in the cover_type_field, we leave it at 1x1 deg resolution and let the model do the automatic interpolation.
5.3 Groundwater field:
For the groundwater residence time field, we do a similarly crude homogenisation. Reading in the groundwater_residence_time_field from the pre-industrial boundary conditions (which is generated from realistic values), we take an area-average mean of all land points, to give a single value that we apply to all grid points. In this case we can simply apply it to all points (even ocean points, since they get ignored). To do this, run the script:
python ave_groundwater.py
This creates an input file groundwater_residence_time_field which is the same as the a55_c1 Eocene example.
6 Greenhouse gases
There are 3 greenhouse gas forcing files which you might want to change. They are:
- co2_gblannualdata: CO2 forcing input (ppm). These specify dates and CO2 level. In my example, I have two dates, 1859.5 and 1860.5, and both specify the same CO2 level for the 1xCO2 run (286.0 ppm). If you want to change the CO2, simply change the 286.0 to a new value on both lines. The model reads in CO2 based on a calendar input of January 1860. The dates specified in co2_gblannualdata are chosen simply to “cover” January 1860. (The “1860” doesn’t increment in time, unless you specify that in the input.nml file.)
- ch4_gblannualdata: CH4 forcing input (ppb). As with CO2, if you change the value, change it on both lines of the input file.
- n20_gblannualdata: N20 forcing input (ppb). Same as CO2 and CH4, change the value on both lines.
7 Aerosol forcing
Aerosol forcing can be left unchanged from the existing a55_c1 run. Leaving this section here as a placeholder.
8 Ocean tidal forcing boundary conditions
For this section we use the following files:
8.1 Make tidal amplitude file
The input file tideamp.nc specifies a tidal mixing amplitude to be used by the model. We generate this using a pre-industrial input file tideamp_cm2.1_pi.nc, and homogenise it for all ocean seafloor values. This is performed by running the script:
python homog_tid_Her.py
This creates the output tideamp_homog_Her.nc. Note, there are two output variables in the output, tideamp, which is a homogenised tidal mixing amplitude, and wet, which is an ocean land-sea mask. In the variable wet, 1=ocean, 0=land.
If you modify the land-sea mask, you need to update the variable wet to match your new land-sea mask. You don’t need to modify the tideamp variable because it is the same everywhere.
The output file, in this case tideamp_homog_Her.nc should be renamed to tideamp.nc as an input file to the model.
8.2 Make roughness amplitude
Similar to the tideamp.nc file, we also create a homogenised roughness_amp.nc input file, from a pre-industrial input file. If you run the script
python homog_rough_Her.py
This will create the output roughness_amp_homog_Her.nc. This file contains variables roughness_amp, which is uniform everywhere and can be left unchanged, and wet, which you need to update to your new land-sea mask as with tideamp. Finally, rename the output file roughness_amp.nc, as an input file to the model.
8.3 Make temperature salinity restart file
I have included the temperature salinity restart file ocean_temp_salt.res.a55_c1.nc, which has been spun-up from a 6000-year Eocene simulation at 1xCO2. In this file, grid cells that are not ocean are simply set to zero. To make it applicable to a new topography, i.e. where new ocean grid cells are created, a simple approach is to use nearest neighbour interpolation. I have included a script to do this, you can run:
python extrapolate_temp_salt.py
This simply extrapolates the existing temperature and salinity fields so that they are well-defined everywhere. I.e. You don’t accidentally set temperature or salinity to zero. One can also define a temperature-salinity restart file by very idealised means, e.g. set the ocean to 10 degC everywhere, and salinity to 34.7 psu everywhere. Such an input file can work for a new simulation, if spinning up to a new equilibrium is not an issue.