adding a new feature to CTSM: excess ground ice physics

[lca041]

Purpose with this issue:
We propose to add excess ground ice physics into the ctsm5.1.dev014 for better solving soil physics and hydrological impacts associated with the thawing process in the permafrost regions. This will imply changes in the following 8 files:

• bld/namelist_files/namelist_definition_ctsm.xml
• src/biogeophys/SoilHydrologyMod.F90
• src/biogeophys/SoilTemperatureMod.F90
• src/biogeophys/TemperatureType.F90
• src/biogeophys/TotalWaterAndHeatMod.F90
• src/biogeophys/WaterStateType.F90
• src/main/clm_varctl.F90
• src/main/controlMod.F90

Brief science information:
The current CTSM does not include excess ground ice for the permafrost region, which inevitably leads to biases in calculating the soil physics and hydrological change along with the permafrost degradataion process (Burn and Nelson, 2006). We propose to add a feature of excess ice physics, which was initially developed based on CLM4.5 (Lee et al., 2014). Compared to the original development, we suggest to add a namelist entry and switches in the code so ctsm users could to turn on and off the excess ice physics easily.

Pseudocode for all file changes:

• bld/namelist_files/namelist_definition_ctsm.xml: adding an namelist entry in order to turn on and off excess ice physics upon users need.
• src/biogeophys/SoilHydrologyMod.F90: involving the added excess ice into the calculation of total volume of soil ice.
• src/biogeophys/SoilTemperatureMod.F90: calculate the melt of excess ice and its impact to soil temperature and extra soil water.
• src/biogeophys/TemperatureType.F90: resetting the initial soil temperature to a freezing temperature.
• src/biogeophys/TotalWaterAndHeatMod.F90: involving the added excess ice into the calculation of total mass of soil ice.
• src/biogeophys/WaterStateType.F90: setting the initial condition of excess ice by reading variables from surface data (surface data file is revised to include excess ice initial condition variables).
• src/main/clm_varctl.F90: adding a switch to turn on and off excess ice, the default value of which is “false”.
• src/main/controlMod.F90: adding a flag to read the excess ice physics namelist entry.

Long term goal:
To add the suggested feature to the master branch of CTSM.

People involved:
@lca041, Hanna Lee (hale@norceresearch.no), @kjetilaas, @sunnivin, @dlawrenncar, @swensosc. ​

[mvdebolskiy]

Hi @billsacks, @dlawrenncar, @swensosc.
I am working on developing code for this issue.
A "merged" implementation of excess ice can be found at @lca041's fork https://github.com/lca041/ctsm/tree/clm5.0.dev92_exice.
The place where I am developing this implementation into the current CTSM/master can be found here: https://github.com/mvdebolskiy/CTSM/tree/excess_ice.
In addition, I am maintaining a hackmd doc, where I keep a detailed summary of this code (for future documentation and easier understanding: https://hackmd.io/@mvdebolskiy/r1EqT998O.
Right now there few problems that need to be decided on about the code, so it properly fits in the clm:

1. Variable initialization for the cold start
2. Controlling excess ice execution and interference with the existing code in /src/biogeophys/SoilTemperatureMod.F90

Cold Start

Right now the public variables that are associated with the excess_ice execution are appeneded to the /src/biogeophys/WaterStateType.F90:

real(r8), pointer :: excess_ice_col     (:,:)   ! ammount of excess ice per soil layer in a column (kg/m2) 
real(r8), pointer :: init_exice             (:,:)   ! initial value (begining of a timestep) of col excess ice per soil layer (kg/m2)
real(r8), pointer :: exice_melt_lev     (:,:)   ! excess ice melting per soil layer (m)
real(r8), pointer :: exice_melt           (:)    ! excess ice melting per soil column (m) 

The method WaterStateType%InitCold has been modified to read a 2d (lat/lon) variable exice_init from a surface dataset (fsurdat). The single value per column is then assigned to soil layers (excess_ice_col variable) below a hardcoded depth (4 m). The rest of the soil column is assigned 0 excess ice amount. Then init_exice is assigned a value of excess_ice_col.
All the other additional variables are assigned to 0.

This is not consistent with the rest of the model.
I had a problem with initializing the model in @lca041's fork (given fixing some typos to make it to compile), however, the problem might have been due to the use_excess_ice not being broadcasted to all the processors (there is no call mpi_bcast for this variable in his version).

In any case, there are couple of questions:

• The "initial" values are in fact initial though since excess ice is not refreezing in current implementation, they can be viewed as maximum values. Should this remain in fsurdat?
• If the initial values should be in fsurdat shouldn't they be read the same way soil parameters are read in (/src/clm_instMod.F90) clm_instInit()? F.e. after initialization of soilstate_inst a separate subroutine is called SoilStateInitTimeConst() that updates variables in the soilstate_inst. Excess ice can be initialized to NaNs or spval within the WaterStateType%InitCold() and then the values are read from fsurdat and updated later when SoilHydrologyInitTimeConst() is called.

SoilTemperatureMod.F90

The flow is controlled with use_excess_ice flag and associated namelist entry.
The allocation and initialization of the additional variables in WaterStateType are already controlled by it in @lca041's fork. However, the code in /src/biogeophys/SoilTemperatureMod.F90 doesn't use it.
What would be the proper way to control the execution of excess ice code in the /src/biogeophys/SoilTemperatureMod.F90 so it doesn't interfere when the use_excess_ice is false?

1. To have a flagged checked every time excess ice code is executed within existing subroutines?
2. Add a single flag check and and interface to override the usage of SoilTemperatureMod to a separate module that contains existing excess ice implementation? (a lot of duplication, separate file)
3. Add a few flag checks in the SoilTemperature() subroutine where modified private subroutines are used? (no additional files, some duplication)
4. Other way?

Right now SoilTemperature() subroutine has local variables that are used in calculations of excess ice effects (dz2, zi2, z2) that have to be passed to private subroutines (like SoilThermProp()). This is done by passing them as arguments.

Summary

Decisions have to be made:

1. What should be the way excess ice variables are initialized during the cold start?
2. What is the preferred structure for the bulk of the excess ice implementation in /src/biogeophys/SoilTemperatureMod.F90?

[swensosc]

For 1., I think the best way to proceed is to use the input streams functionality. The streams functionality will read in input data from a separate file, and the spatial interpolation will be handled by the model (in this case, there will not be time varying fields). That way, the spatial resolution of the input data does not have to be modified for every surface data file. The 'use_excess_ice' flag can be used to determine whether to read from an excess ice stream file. If it is false, the excess ice variable can be set to zero. To see some examples of how the streams files are used, look for the prescribed soil moisture stream, or the dynamic lai stream. Erik should also be able to give advice on how to set up the stream for the initial condition.

For 2., I think the code can/should be written such that when excess ice is zero, the model behavior reverts to the original code. Specifically, variables like z2 (btw, a more descriptive variable name should be used) and the variables related to the soil thermal properties should revert to z when excess ice is zero. That way, the flag itself is not needed in those routines. Additionally, that is the behavior that should occur for an excess ice simulation after the excess ice has melted, so the case with no initial excess ice, and all excess ice melted should be identical (aside from transient differences resulting from the time when excess ice was nonzero).

[mvdebolskiy]

For 1., I think the best way to proceed is to use the input streams functionality. The streams functionality will read in input data from a separate file, and the spatial interpolation will be handled by the model (in this case, there will not be time varying fields). That way, the spatial resolution of the input data does not have to be modified for every surface data file. The 'use_excess_ice' flag can be used to determine whether to read from an excess ice stream file. If it is false, the excess ice variable can be set to zero. To see some examples of how the streams files are used, look for the prescribed soil moisture stream, or the dynamic lai stream. Erik should also be able to give advice on how to set up the stream for the initial condition.

Should it be within the existing SoilMoistureStreamsMod (I assume the version at /src/cpl/nuopc/ is the most current one).
Or should it be a separate module, since procedures PrescribedSoilMoistureAdvance and PrescribedSoilMoistureInterp are not needed for excess_ice?

For 2., I think the code can/should be written such that when excess ice is zero, the model behavior reverts to the original code. Specifically, variables like z2 (btw, a more descriptive variable name should be used) and the variables related to the soil thermal properties should revert to z when excess ice is zero. That way, the flag itself is not needed in those routines. Additionally, that is the behavior that should occur for an excess ice simulation after the excess ice has melted, so the case with no initial excess ice, and all excess ice melted should be identical (aside from transient differences resulting from the time when excess ice was nonzero).

Would not that be too taxing computationally for use cases where excess_ice is switched off? Should I try to implement this version and a separate-module version and test the performance difference?

@ekluzek what's the best source on making data stream dataset?

[swensosc]

It should be a separate module for excess ice. If you think that it might add significant computational cost, then testing the performance would be a good idea.

[ekluzek]

Should it be within the existing SoilMoistureStreamsMod (I assume the version at /src/cpl/nuopc/ is the most current one).
Or should it be a separate module, since procedures PrescribedSoilMoistureAdvance and PrescribedSoilMoistureInterp are not needed for excess_ice?

I would create a new module for excess ice streams, named something like excessSoilIceStreams.F90. This should go in the src/cpl/nuopc directory. Another module should handle applying the excess ice streams to the model state. You should also create a stub for MCT in src/cpl/mct that doesn't do anything, but just provides a stub interface so you can build with MCT. For clarity you can add an abort if excess ice is tried to be used for MCT and not NUOPC.

For 2., I think the code can/should be written such that when excess ice is zero, the model behavior reverts to the original code. Specifically, variables like z2 (btw, a more descriptive variable name should be used) and the variables related to the soil thermal properties should revert to z when excess ice is zero. That way, the flag itself is not needed in those routines. Additionally, that is the behavior that should occur for an excess ice simulation after the excess ice has melted, so the case with no initial excess ice, and all excess ice melted should be identical (aside from transient differences resulting from the time when excess ice was nonzero).

Would not that be too taxing computationally for use cases where excess_ice is switched off? Should I try to implement this version and a separate-module version and test the performance difference?

I doubt it will matter, as long as the switch is at a high level. If the switch is imbedded deep within a patch loop or something -- then yes maybe it would make a small difference. I would just look at the performance difference with the new code compared to the baseline code version.

@ekluzek what's the best source on making data stream dataset?

I would look at the examples of the existing stream datasets. You might especially look at the ch4 finundated stream dataset, since it's similar in only having a single time. The soil moisture stream dataset will also be useful as your dataset willl be in 3D. The MCT implementation for 3D fields requires the stream dataset to be on the same resolution as the model is running at -- but NUOPC doesn't require that. @mvertens is implementing more streams in NUOPC so when that work is done, it should be helpful for you to see how to do this.

[mvertens]

@erik @mvdebolskiy - I have an upcoming PR that I am finalizing that translates all of the mct stream functionality to calling cdeps share code.
If you want to look at the code before I issue a PR - its the branch feature/stream_refactor in https://github.com/mvertens/CTSM

[mvertens]

@erik @mvdebolskiy - just to clarify - all of the MCT stream functionality is sitll in place. Its just if you use CMEPS, you will use the new CDEPS stream functionality.

[mvdebolskiy]

@ekluzek the datastream will be 2d and the redistribution of the data within the column will be handled by a couple of parameters. I've already made .nc files with 0.9x1.25 resolution (since there was no example of interpolation in /src/cpl/nuopc/SoilMoistureStreamsMod.F90 since the last time I've checked. But I am trying to learn more about Meshes in ESMF/CDEPS and how to make the original dataset (0.1x0.1 degrees) to get interpolated in runtime. If I understand correctly, I will just need to create a meshfile specific to my resolution.

@mvertens does that mean that I can write with the old mct instructions? Also, what would be the proper version of ESMF for me to get on sigma2 machine? If I understood correctly 8.1 should be enough. Is it ok to build it with intel-2019b compiler?

[mvertens]

@mvdebolskiy - I would encourage you to use the new CDEPS framework since that is what we are moving to with all the new science being added to the model. My comment is to clarify that the older MCT streams will be kept around for backwards compatibility (they will now reside in the src/cpl/mct) but no new development will be done with them. We are currently using intel19.1.1. @jedwards4b - can you please respond regarding intel-2019b?

[jedwards4b]

I'm not sure what intel-2019b is but it's probably okay as long as you are consistent with the cesm build.

[mvdebolskiy]

Hi, sorry for the delay, I managed to run latest master with nuopc and cdeps/cmeps on the Norwegian machines with intel-2021b.
I've also managed to get initialization with data streams in the same way it is done with ch4finundated, however, this means that NLFilename has to be passed into the WaterStateType%InitCold, I did it as an optional argument for the subroutines that are calling it from WaterType, you can look at the latest commits in my fork. Right now everything works only on 0.9x1.25 resolution and I will need some help to make mesh file and data file for the original 12.5 by 12.5 km dataset

[ekluzek]

@mvdebolskiy two issues that talk about mesh files are here: #1397and #1513. One of the tools we use is an ESMF tool that converts files from SCRIP grid format to ESMF mesh format ESMF_Scrip2Unstruct. It's an executable built with the ESMF library under the "bin" directory. The ESMF documentation gives the specifications for contents of the SCRIP grid files and the ESMF mesh files. The SCRIP grid format is simpler as it just needs the gridcell center latitude/longitudes and the latitude/longitudes of the four vertices.