Setting up (high-res sparse) regional-grid CTSM simulations

[slevis-lmwg]

Context
This started as issue #1773 from which I linked a google doc (now out of date and obsolete) with step-by-step instructions intended initially for a case that @XiulinGao needed to run. To the google doc I had added step-by-step instructions for a simpler case that @jkshuman needed to run. I am copying the step-by-step instructions for the two cases here, and periodically including updates that are not going on the google doc.

Another use-case is the PPE sparse grid that @ekluzek is working on in #1731.
Documentation for non-sparse regional cases appears here.

These workflows should likely ultimately appear in the CTSM User's Guide.

Note
If subsetting from a global fsurdat file while not changing land properties or the land mask, follow the instructions for running regional FATES at regular 1-degree resolution (see Table of Contents).

Contents
Instructions for running CTSM-FATES on a sparse regional grid with WRF datm driver files
Instructions for running regional CTSM-FATES at regular 1-degree resolution (alternate workflow appears here)
Addendum for running CTSM-FATES single-point with WRF datm driver files

[slevis-lmwg]

Instructions for running CTSM-FATES on a sparse regional grid with WRF datm driver files

Slevis placing files from this work in /glade/work/slevis/data_UCB/ca_fates_input_data/wrf_CA_sparse_grid
though Xiulin’s datm data is in
/glade/work/slevis/data_UCB/ca_fates_input_data/datm_from_wrf
and fsurdat is here
/glade/work/slevis/data_UCB/ca_fates_input_data/fsurdat/wrf_CA

Slevis: From github issue #1773 I compiled the following series of steps to set up simulations on cheyenne. Let us know if you encounter issues due to the switch to derecho, and we will attempt to update the instructions accordingly.

A) Generate datm files from wrf output

SKIP IF using one of the CTSM datm datasets.

To run this high resolution regional case, we used downscaled 9 km WRF data for the western US (from Stefan Rahimi at UCLA). WRF output files differ significantly from CLM forcing files given the data format. In order to run CLM-FATES with WRF forcing, one has to reprocess original WRF files first. Xiulin Gao wrote a R script to do this(https://github.com/XiulinGao/FATES-CA-grassland-exp/blob/main/scripts/WRF-CLM-forcing.Rmd). If you use Python instead and want to write your own script, a few things to notice before you start:

1. WRF grids are curvilinear grids with 2D longitude and latitude. Actual LONG and LATI cannot be found in the forcing file, one needs a WRF output metadata to obtain the 2D XLONG and XLAT variables.
2. Each WRF file is an hourly file. In order to convert them into monthly files, one needs to concatenate those files. The suggestion is to make an annual file for each year. Columns in the annual file should include XLONG, XLAT, time, and all other climate variables you need to make your forcing.
3. With the annual file from step 2, one can calculate precipitation rate and wind speed. WRF writes out precipitation as cumulative precipitation, so D(RAINC + RAINNC)/Dt should give you the precipitation rate, and then convert to mm/s. Wind speed is calculated as sqrt(U*U + V*V). For wind speed, no grid rotation is required.
4. When writing out a monthly netcdf file, make sure the time stamp is in the middle of each observation window. And rename climate variables to the correct names required by CLM. See detailed requirements for preparing your own forcing to run CLM in the user guide.

B) Generate landmask.nc

SKIP IF subsetting from global to regional while not changing land properties or the land mask.

In the specific example we wish to limit the regional simulation to the grid cells specified by the landmask variable in landmask.nc (in this example, wrf_CA_ann_grass_bilinear.nc which was generated by Xiulin). The file must contain

• Dimensions lon and lat (naming is flexible)
• Variables lsmlon and lsmlat (naming is flexible), landmask and mod_lnd_props
• mod_lnd_props = landmask

This file’s grid matches the grid found in the datm files discussed previously. Circularity in some choices in the current process may make this a requirement. In particular, generating the fsurdat file ← from the mesh file ← from a datm file while modifying the mesh file from the landmask file. Theoretically, though the datm grid may span a larger region (even global) than the CTSM will run on.

If you already have your fsurdat file, generate landmask.nc by copying your fsurdat to landmask.nc, then renaming the mask variable to landmask, then copying the fsurdat mask variable to landmask.nc again, and renaming that to mod_lnd_props. These steps can be accomplished with netcdf operators:

• cp fsurdat.nc landmask.nc
• ncrename -v PFTDATA_MASK,landmask landmask.nc
• ncks -A -v PFTDATA_MASK fsurdat.nc landmask.nc
• ncrename -v PFTDATA_MASK,mod_lnd_props landmask.nc

C) git clone CTSM

git clone https://github.com/ESCOMP/CTSM.git
cd CTSM
./manage_externals/checkout_externals
module unload python
module load conda
./py_env_create
conda activate ctsm_py

D) Generate mesh file from a 2d lat/lon file containing your regional grid & update mask…

NB1. Alternate instructions exist here.
NB2. A CTSM fsurdat (or domain?) file (if you have one or can create one with the subset_data tool) could serve as the 2d lat/lon file with which to generate the mesh file in steps (2-3 below).

1. SKIP (1) IF you have a .nc file with lat/lon variables in 1d or 2d (dimensions).
   If lat/lon have more than 2d, execute a command like this to reduce to 2d variables:
   ncwa -a <unnecessary_dimension> <file_with_lat_lon_3d>.nc <file_with_lat_lon_2d>.nc
2. Generate scrip.nc from <file_with_lat_lon_2d>.nc or from 1981-01.nc if skipped step (1)
   ncks --rgr infer --rgr scrip=scrip.nc 1981-01.nc foo.nc
   where foo.nc contains metadata and will not be used.
3. Generate mesh file from scrip.nc

module load esmf
ESMF_Scrip2Unstruct scrip.nc lnd_mesh.nc 0

NB. This mesh file’s mask = 1 everywhere.
4) Update the mesh file’s mask according to the landmask variable
4a) EITHER DO THIS (shortcut that allows you to skip 4b)
If not changing the land mask, execute only this line:
ncap2 -s 'elementMask(:)=0' lnd_mesh.nc mask_mesh.nc
Now you can skip step (4b).
4b) OR DO THIS

cd tools/modify_mesh_mask
./mesh_mask_modifier modify_mesh_wrf_CA.cfg

Before executing the line above, we set options as follows in the .cfg file:

mesh_mask_in = lnd_mesh.nc
mesh_mask_out = mask_mesh.nc
landmask_file = wrf_CA_ann_grass_bilinear.nc
lat_dimname = lat
lon_dimname = lon
lat_varname = lsmlat
lon_varname = lsmlon

The mask in mask_mesh.nc represents the inverse of the land mask, so
land points = 0
non-land points = 1

E) git checkout ctsm5.2.mksurfdata

SKIP IF you don’t still need to generate an fsurdat file.

Use of the mksurfdata_esmf tool requires the ctsm5.2.mksurfdata branch. One could but should not have to clone the ctsm a second time to check out this branch. When merged with the main trunk of the CTSM (not as of the last edit of this post), users will not need to execute the extra step of checking out this branch:

git checkout ctsm5.2.mksurfdata
./manage_externals/checkout_externals
module unload python
module load conda
./py_env_create
conda activate ctsm_py
cd tools/mksurfdata_esmf

F) Make CTSM surface dataset (fsurdat file) using the mksurfdata_esmf tool

SKIP IF you don’t still need to generate an fsurdat file.

cd tools/mksurfdata_esmf
./gen_mksurfdata_build.sh
./gen_mksurfdata_namelist.py --start-year 2000 --end-year 2000 --model-mesh-nx 147 --model-mesh-ny 151 --nocrop --model-mesh lnd_mesh.nc --res <name_your_grid>
./gen_mksurfdata_jobscript_single.py --number-of-nodes 24 --tasks-per-node 12 --namelist-file <file_generated_in_prev_step>
qsub mksurfdata_jobscript_single

G) Run the CTSM
If you generated a new fsurdat file in (F), then keep working in the ctsm5.2.mksurfdata branch. Otherwise you may keep working in
the CTSM directory that you obtained in step (C).

1. ./create_newcase --case <case_name> --compset I2000Clm50Fates --res CLM_USRDAT --run-unsupported
2. For faster turnaround, consider changing (see slow regional runs with nuopc vs mct #1907 for a detailed investigation)

• env_mach_pes.xml to run on more nodes and processors:
  NTASKS → -4 (Change for all components except ESP)
  NTASKS_PER_INST → 144 (Change for all except ESP)
• env_workflow.xml to run on the regular queue:
  JOB_WALLCLOCK_TIME → 12:00:00 (from 6:00:00 in line 32)
  JOB_QUEUE → regular (from share in line 37)
  PIO_TYPENAME to pnetcdf in env_run.xml

3. ./case.setup
4. Add this line to user_nl_clm:
   fsurdat = '<path and name of fsurdat file discussed above>'
5. Update lines in env_run.xml:

./xmlchange LND_DOMAIN_MESH='/path/to/your/lnd_mesh.nc'
./xmlchange ATM_DOMAIN_MESH='/path/to/your/lnd_mesh.nc'
./xmlchange MASK_MESH='/path/to/your/mask_mesh.nc'
./xmlchange DATM_YR_START=<first year of datm data>
./xmlchange DATM_YR_END=<last year of datm data>

6. Using the regional datm dataset
   SKIP IF using one of the CTSM datm datasets.
   Update user_nl_datm_streams to include the next three lines three times. In the second set of three, replace Solar with Precip and in the third set replace Solar with TPQW. In the datafiles lines include all the years and months in order (only one year shown here as an example):

CLMGSWP3v1.Solar:mapalgo = nn
CLMGSWP3v1.Solar:meshfile = lnd_mesh.nc
CLMGSWP3v1.Solar:datafiles = /glade/work/slevis/data_UCB/ca_fates_input_data/datm_from_wrf/1981-01.nc,/glade/work/slevis/data_UCB/ca_fates_input_data/datm_from_wrf/1981-02.nc,/glade/work/slevis/data_UCB/ca_fates_input_data/datm_from_wrf/1981-03.nc,/glade/work/slevis/data_UCB/ca_fates_input_data/datm_from_wrf/1981-04.nc,/glade/work/slevis/data_UCB/ca_fates_input_data/datm_from_wrf/1981-05.nc,/glade/work/slevis/data_UCB/ca_fates_input_data/datm_from_wrf/1981-06.nc,/glade/work/slevis/data_UCB/ca_fates_input_data/datm_from_wrf/1981-07.nc,/glade/work/slevis/data_UCB/ca_fates_input_data/datm_from_wrf/1981-08.nc,/glade/work/slevis/data_UCB/ca_fates_input_data/datm_from_wrf/1981-09.nc,/glade/work/slevis/data_UCB/ca_fates_input_data/datm_from_wrf/1981-10.nc,/glade/work/slevis/data_UCB/ca_fates_input_data/datm_from_wrf/1981-11.nc,/glade/work/slevis/data_UCB/ca_fates_input_data/datm_from_wrf/1981-12.nc

7. This step addresses restart failures (i.e. CONTINUE = 'TRUE') due to the data atmosphere (datm) model reading nx ny values equal to 1 from env_build.xml. One way around this is to set

ATM_NX and LND_NX to 147
ATM_NY and LND_NY to 151

in env_build.xml so that datm knows the correct size of the regional grid.

A more robust solution allows all new cases to include the correct nx ny values as posted here.

8. ./case.build

9. ./case.submit

[slevis-lmwg]

Instructions for running regional CTSM-FATES at regular 1-degree resolution

slevis: I wrote the following instructions while testing on izumi with this ctsm5.2 tag:
alpha-ctsm5.2.mksrf.10_ctsm5.1.dev111-8-g858539a61
The same instructions should also work on cheyenne and with ctsm5.1.dev111. Let us know if you encounter issues due to the switch to derecho, and we will attempt to update the instructions accordingly.

A) Create case
From your ctsm/cime/scripts directory:
./create_newcase --case <case_name> --res CLM_USRDAT --compset COMPSET=2000_DATM%GSWP3v1_CLM50%FATES_SICE_SOCN_RTM_SGLC_SWAV --mpilib mpi-serial --run-unsupported

NOTE: You need "--mpilib mpi-serial" for simulations running on one processor, which is most common for single-point simulations. For regional simulations, omit this and select a pe-layout that works for your region.

B) Generate fsurdat file
From your ctsm/tools/site_and_regional:
./subset_data region --lat1 -38 --lat2 14 --lon1 275 --lon2 330 --reg south_america --create-surface

NB. The --create-mesh option of this tool did not work for a while. The same functionality is accomplished here externally with step (C).

C) Generate mesh files
In ctsm/tools/site_and_regional/subset_data_regional:

1. module load nco # on derecho OR module load tool/nco/4.7.5 …on izumi
2. Generate scrip.nc: ncks --rgr infer --rgr scrip=scrip.nc <fsurdat_generated_in_B> foo.nc
3. Generate lnd_mesh.nc:

module load esmf
ESMF_Scrip2Unstruct scrip.nc lnd_mesh.nc 0`   # on derecho or casper

OR /project/esmf/PROGS/esmf/8.2.0/mpiuni//gfortran/9.3.0/bin/binO/Linux.gfortran.64.mpiuni.default/ESMF_Scrip2Unstruct scrip.nc lnd_mesh.nc 0 # on izumi
This assumes all grid cells are land in this case so mask is 1. If your mask contains both land and ocean, then look at instructions in the previous example for generating a landmask.nc file and use that file in step 2 to generate the scrip.nc and continue from there.

D) Run the CTSM
In your case directory:

For faster turnaround, consider changing (see #1907 for a detailed investigation)
env_mach_pes.xml to run on more nodes and processors:
NTASKS → -4 (Change for all components except ESP)
NTASKS_PER_INST → 144 (Change for all except ESP)
env_workflow.xml to run on the regular queue:
JOB_WALLCLOCK_TIME → 12:00:00 (from 6:00:00 in line 32)
JOB_QUEUE → regular (from share in line 37)
PIO_TYPENAME to pnetcdf in env_run.xml

./case.setup
Add line to user_nl_clm:
fsurdat = ‘<fsurdat_generated_in_B>’
Update lines in env_run.xml:

./xmlchange LND_DOMAIN_MESH='/path/to/your/lnd_mesh.nc'
./xmlchange ATM_DOMAIN_MESH='/path/to/your/lnd_mesh.nc'
./xmlchange MASK_MESH='/path/to/your/lnd_mesh.nc'

Before you build the case and submit the run, a note about restarts:
If you find that your restart fails, you are likely to address the failure with the equivalent of step 7 from the previous post. The nx, ny for regional 1-degree will be different than in the previous post. You can obtain the values from the fsurdat file generated in step B.

./case.build
./case.submit

[slevis-lmwg]

@olyson thank you for catching that. I changed the instructions to be more generic.

[jcyang0]

Hi @slevis-lmwg, just last week I was able to successfully run step C.3 above on derecho, but now when I try to access /glade/u/apps/ch/opt/esmf-netcdf/8.0.0/intel/19.0.5/bin/bing/Linux.intel.64.mpiuni.default/ESMF_Scrip2Unstruct it says "Permission denied". Has something changed quite recently?

[ekluzek]

Hmmm. I think the ESMF modules are setup a little differently now on Derecho. And actually the above might have been setup on the older Cheyenne space and just became broken as they phase those disk spaces out...

I played around with modules and found a version here for example...

/glade/u/apps/derecho/23.09/spack/opt/spack/esmf/8.6.0/cray-mpich/8.1.27/oneapi/2023.2.1/7haa/bin/ESMF_Scrip2Unstruct

This is with loading the following...

module list

Currently Loaded Modules:

1. ncarenv/23.09 (S) 3) nco/5.1.9 5) ncview/2.1.9 7) intel/2023.2.1 9) cray-mpich/8.1.27 11) netcdf-mpi/4.9.2
2. craype/2.7.20 4) conda/latest 6) cdo/2.3.0 8) ncarcompilers/1.0.0 10) hdf5-mpi/1.12.2 12) esmf/8.6.0

[ekluzek]

@slevis-lmwg we should update the above instructions for Derecho and for this new path.

[slevis-lmwg]

Thank you @ekluzek and @jcyang0.
I updated the path in the above instructions.

[jkshuman]

@slevisconsulting the motivation for putting these instructions into a discussion was to provide a target for users looking for a workflow for regional subsetting. There seem to be at least two cases with different complexity: a basic regional subsetting and the more complex high-resolution sparse grid subsetting. As these are different and distinct workflows, it would make sense to split these into two separate discussion topics. They can certainly have a note to cross-reference to the other workflow, but for sharing with users in the community it would be easier to separate them as two discussion workflows.

[slevis-lmwg]

@jkshuman from follow-up that I have encountered on this topic, I now see that it's easier for me to keep all the information in this one discussion, rather than dispersed. This does not preclude the opening of a separate discussion with the generalized instructions that you had in mind. I hope that's a reasonable compromise.

[jkshuman]

@slevisconsulting adding a link to the generalized instructions from @adrifoster if it useful. #1833

[slevis-lmwg]

Thank you @jkshuman

[niuhanlin]

@slevisconsulting in the fourth step of mesh_generated_in_C I had a WARNING:
ncap2: WARNING assign(): Var being read and written in ASSIGN elementMask.
Is this normal? Or what do I need to do?

[slevis-lmwg]

@niuhanlin I couldn't remember whether I saw this warning when I originally documented this process, so I repeated step (C4) "Generate mask_mesh.nc: ncap2 -s 'elementMask(:)=0' lnd_mesh.nc mask_mesh.nc" on izumi (cheyenne is down for maintenance). Yes, I got the same warning. I checked whether the output file mask_mesh.nc looked the same as before, and I am pleased to say that it is identical. This tells me that you do not need to do anything, and you can use the generated mesh file.

[niuhanlin]

Thank you for your test. I have used the simulation you recommended, but there are errors in between, so if you have suggestions I can keep trying.
Attached is a link to the content of my error report.
#975

[slevis-lmwg]

Similar to the first post in this discussion, a user attempted a single-point simulation using WRF output on a curvilinear grid as their datm driver. The difference from the first post is that this user wants single-point rather than a region. I found that the following changes were necessary for this to work:

Do not specify mesh files in env_run.xml and user_nl_datm_streams. In user_nl_datm_streams, set
CLMGSWP3v1.Solar:meshfile = none
and repeat these lines for Precip and TPQW.

(The user already had a single-point fsurdat file for their case.)

[XiulinGao]

Similar to the first post in this discussion, a user attempted a single-point simulation using WRF output on a curvilinear grid as their datm driver. The difference from the first post is that this user wants single-point rather than a region. I found that the following changes were necessary for this to work:

Do not specify mesh files in env_run.xml and user_nl_datm_streams. In user_nl_datm_streams, set CLMGSWP3v1.Solar:meshfile = none and repeat these lines for Precip and TPQW.

(The user already had a single-point fsurdat file for their case.)

Adam Brown found that if running single-point simulation using WRF forcing and let the model subset the forcing data automatically, the model was not able to grab the correct forcing data given the lon and lat user defined. He found that the forcing being used is somewhere in the desert region, very dry and warm climate, which is similar to what I saw when didn't set up the right mapping algorithm during the regional run, it keeps giving the same forcing values (warm and dry) to all grid cells.

running multiple point-simulations with different lons and lats to check whether the model is using the exact same forcing for all runs will tell us if this is the correct way to use WRF forcing for single-point simulaitons.

[XiulinGao]

Recently running a case using transient historical CO2 with recycled WRF-reanalysis climate driver, ran into issues at the beginning. But after trying and changing model settings many times, I was able to get the model to run (but now facing a mass conservation check error, which should not be related to how I set up the model). As requested by @slevis-lmwg , I am documenting the changes I made to get the case going here:

1. ./xmlchange CLM_CO2_TYPE="diagnostic"
2. ./xmlchange CCSM_BGC="CO2A"
3. ./xmlchange DATM_CO2_TSERIES="20tr"
4. ./xmlchange DATM_PRESAERO="hist"
5. ./xmlchange DATM_PRESNDEP="hist"
6. ./xmlchange DATM_TOPO="observed"

I also used this compset "20TR_DATM%GSWP3v1_CLM50%FATES_SICE_SOCN_SROF_SGLC_SWAV" and seems 4-6 changes are set as above by the compset as default, so you can even skip 4-6.

As the climate driver is not lined up with the CO2, and I recycle the climate driver, so in user_nl_datm_stream file, I also specified "co2tseries.20tr:year_first = yyyy" ,"co2tseries.20tr:year_last = yyyy","co2tseries.20tr:year_align = yyyy","co2tseries.20tr:taxmode = extend". just replace yyyy with whatever year you aim to use.