# Adapting GFDL CM2.1 coupled model for deep time paleoclimate experiments **URL:** \ **Category:** Earth System\ **Tags:** experiments, user-story\ **Created:** [1 June 2023 05:02 UTC](https://forum.access-hive.org.au/t/adapting-gfdl-cm2-1-coupled-model-for-deep-time-paleoclimate-experiments/826 "2023-06-01T05:02:19Z")\ **Posts on this page:** 1\ **Page:** 1
**Author:** ![dkhutch](https://sea2.discourse-cdn.com/flex020/user_avatar/forum.access-hive.org.au/dkhutch/32/510_2.png) [@dkhutch](https://forum.access-hive.org.au/u/dkhutch)\ **Post date:** [1 June 2023 05:02 UTC](https://forum.access-hive.org.au/t/adapting-gfdl-cm2-1-coupled-model-for-deep-time-paleoclimate-experiments/826/1 "2023-06-01T05:02:19Z")
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](https://github.com/dkhutch/cm2.1_paleo) 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: > **[Dropbox](https://www.dropbox.com/sh/5hdzsrrsehydo6c/AABoEgQPBqm_Do4CgNW_t9-Fa?dl=0)** # 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: ```auto 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: ```auto ./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: ```auto ./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: ```auto atmos_mosaic_tile1Xland_mosaic_tile1.nc land_mosaic_tile1Xocean_mosaic_tile1.nc atmos_mosaic_tile1Xocean_mosaic_tile1.nc ``` and the land-sea mask files: ```auto 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 ```auto 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: > **[Dropbox](https://www.dropbox.com/sh/5hdzsrrsehydo6c/AABoEgQPBqm_Do4CgNW_t9-Fa?dl=0)** There is a restart file called ```auto 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: ```auto 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, > **[cm2.1\_paleo/rundirs/a55\_c1 at main · dkhutch/cm2.1\_paleo](https://github.com/dkhutch/cm2.1_paleo/tree/main/rundirs/a55_c1)** > > Sharing tools and input files for creating new boundary conditions using GFDL CM2.1, for paleoclimate simulations - dkhutch/cm2.1\_paleo The main input file **input.nml** contains a namelist called **topography\_nml** , which specifies an input file topography file: ```auto 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: > **[Dropbox](https://www.dropbox.com/sh/k0woecs9pt22cs2/AACK05gZZwkTku9cNcX1c4Tza?dl=0)** ## 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](https://fennerschool.anu.edu.au/research/products/anudem#:~:text=for%20more%20information.-,Drainage%20enforcement%20algorithm,in%20input%20sink%20data%20files.) ) 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: > **[Dropbox](https://www.dropbox.com/sh/dwrzcr0qvkjxeov/AAB8_qAmGQBDEKNNgLkp3dola?dl=0)** ## 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: ```auto 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: ```auto 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. ```auto 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: ```auto 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: ```auto 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: ```auto 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: > **[Dropbox](https://www.dropbox.com/sh/ek6kvev2q321flu/AAC0yMEcB2T-kcgYiXZYfKJCa?dl=0)** ## 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: ```auto 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 ```auto 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: ```auto 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.