mksurfdata_esmf tool: missing percent clay and sand

[aelkouk]

Hi

I followed this tutorial to make high-resolution CTSM surface dataset: Setting up (high-res sparse) regional-grid CTSM simulations · ESCOMP/CTSM · Discussion #1919
The resulting fields for soil texture (percent sand and clay) have no value (same value as ocean) outside the US, e.g., Mexico.
I am not sure if mksurfdata_esmf tool is ready for external use? or maybe am I using the wrong tag?

Thank you!

 

[slevis]

@aelkouk
If you have access to NCAR's cheyenne or casper, pls see sample fsurdat files generated with the new mksurfdata_esmf tool here:
/glade/p/cesmdata/cseg/inputdata/lnd/clm2/surfdata_esmf/ctsm5.2.0
I don't have a good explanation for the issue in your data, given the information that you provided. What tag are you using? What region are you trying to simulate? Are you pointing to the soil data provided by the tool?

The new tool is not released, yet, but it should be ready, so we are interested in confirming that you did not discover a problem.

 

[aelkouk]

@slevis thank you for your reply.
- I used tag: alpha-ctsm5.2.mksrf.13_ctsm5.1.dev111, I also tried an earlier tag: alpha-ctsm5.2.mksrf.09_ctsm5.1.dev090, but it has the same problem
- I am trying to simulate a regional domain of Southwest North America, I have the mesh file on cheyenne: /glade/scratch/elkoukah/inputs/conus404/mesh/lnd_mesh.nc
- Yes the screenshot is from the surface data provided by the tool, please find the tools output here: /glade/work/elkoukah/dev_ctsm/tools/mksurfdata_esmf/

 

[slevis]

Latest tag is alpha-ctsm5.2.mksrf.15_ctsm5.1.dev123, though this shouldn't matter.

I don't have a good explanation for what you're seeing, yet. It may take me a while to troubleshoot this among other priorities, so pls feel free to continue your own troubleshooting in the meantime. If you happen to resolve before I get to it, pls post what you find.

 

[slevis]

First, can you try (or have you already tried) making a global fsurdat file? It would be interesting to know whether you see the same problem.

 

[aelkouk]

I don't have that problem making a global fsurdat file (T42 resolution). Maybe this issue only happens when making a regional fsurdat file?

 

[slevis]

@aelkouk I was discussing your post with @erik and he had a suggestion:
In /tools/site_and_regional, you can run subset_data to subset global fsurdat files to regional. So, if you made a global mesh at the resolution of your regional one, then you could generate a new global fsurdat file with mksurfdata_esmf, and you would subset that using the subset_data tool.

 

[aelkouk]

Thanks a lot! that is a clever work around.

 

[amans]

I am facing the same issue when making regional surface dataset for mid-US from ctsm5.2.mksurfdat/tools/mksurfdata_esmf. The resulting surface file has clipped data for variables PCT_CLAY, PCT_SAND, CFRAG, ORGANIC, ORGC, mapunits, etc. (see attached screenshot). I find the same issue in different resolutions (0.05deg and 0.125deg). The surface, mesh, and logfiles are in /glade/work/amans/data/MRB-CONUS_0.125 and /glade/work/amans/data/MRB-CONUS_0.05.

Checking the logfile, I see:

Attempting to make %sand, %clay, orgc, cfrag, bulk, phaq .....
Input mapunit file is /glade/campaign/cesm/cesmdata/inputdata/lnd/clm2/rawdata/mksrf_soil_mapunits_5x5min_WISE.c220330.nc
Input lookup table file is /glade/campaign/cesm/cesmdata/inputdata/lnd/clm2/rawdata/mksrf_soil_lookup.10level.WISE.c220330.nc
Input mesh/grid file is /glade/campaign/cesm/cesmdata/inputdata/lnd/clm2/mappingdata/grids/UNSTRUCTgrid_5x5min_nomask_cdf5_c200129.nc
reading mesh_i directly in mksoiltex
Reading in mapunit data in mksoiltex
Setting mask in mesh where mapunit data is 0 mksoiltex
before route handle creation mksoiltex
after route handle creation mksoiltex
before call to dynamic mask set creation mksoiltex
after call to dynamic mask set creation mksoiltex
mksoiltex writing out mapunits
WARNING: assigning sand_o = -4 to 99%
WARNING: assigning other sand_o < 0 to 43%
WARNING: assigning clay_o = -4 to 1%
WARNING: assigning other clay_o < 0 to 18%
WARNING: assigning orgc_o = -4 to 1
WARNING: assigning other orgc_o < 0 to 0
WARNING: same warnings for cfrag_o as for orgc_o
WARNING: assigning bulk_o < 0 to 1.5
WARNING: assigning phaq_o < 0 to 7
mksoiltex writing out soil percent sand
mksoiltex writing out soil percent clay
mksoiltex writing out soil organic matter
mksoiltex writing out soil organic carbon content
mksoiltex writing out coarse fragments in soil
mksoiltex writing out soil bulk density
mksoiltex writing out soil pH measured in H2O
Successfully made %sand, %clay, orgc, cfrag, bulk, phaq .....

It looks like the issue is starting from the mapunit file. I looked at the input mapunit file but could't figure out what the issue is.

 

[oleson]

I looked over your setup and I don't see anything wrong with it at first glance, so it seems like it could be a real issue. We'll have to try to repeat the process ourselves and troubleshoot it.

 

[oleson]

I've been able to reproduce this. I've filed an issue here:

PCT_CLAY and PCT_SAND are truncated in the creation of regional surface datasets · Issue #2430 · ESCOMP/CTSM

Brief summary of bug A couple of users have reported that PCT_CLAY and PCT_SAND fields are truncated in the creation of regional surface datasets. See below for what I mean by truncated. The releva...

github.com

 

[oleson]

As suggested above,
an alternative to this is to create a global scrip file at the desired resolution using /glade/work/oleson/release-cesm2.1.3/components/clm/tools/mkmapgrids/mkscripgrid.ncl. Then create a mesh file and a surface dataset, then subset it.

setenv S_LAT -90

setenv N_LAT 90

setenv E_LON 360

setenv W_LON 0

setenv NX 7200

setenv NY 3600

setenv GRIDFILE SCRIPgrid_7200x3600pt_Global_nomask_c240321.nc

ncl mkscripgrid.ncl

Then follow the steps to create a global surface dataset, i.e.,

qcmd -- /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 /glade/work/oleson/release-cesm2.1.3/components/clm/tools/mkmapgrids/SCRIPgrid_7200x3600pt_Global_nomask_c240321.nc lnd_mesh.nc 0

ncap2 -s 'elementMask(:)=0' lnd_mesh.nc mask_mesh.nc

ncdump -v grid_dims glade/work/oleson/release-cesm2.1.3/components/clm/tools/mkmapgrids/SCRIPgrid_7200x3600pt_Global_nomask_c240321.nc (to get nx and ny)

conda activate ctsm_pylib

./gen_mksurfdata_namelist --start-year 2000 --end-year 2000 --model-mesh-nx 7200 --model-mesh-ny 3600 --model-mesh lnd_mesh.nc --res 0.05-conus-nomask

./gen_mksurfdata_jobscript_single --number-of-nodes 24 --tasks-per-node 12 --namelist-file surfdata_0.05-conus-nomask_hist_2000_78pfts_c240321.namelist

This was the resulting global surface dataset for PCT_CLAY:



Then subset the datasets using subset_data in /glade/work/oleson/ctsm5.2.mksurfdata_032124/tools/site_and_regional

 

[amans]

Thank you @oleson.
About the workaround, I have successfully created the global surface dataset. However, I have some issues while subsetting. It seems to take a long time to subset the global surface dataset. For my case (0.05deg/mid US), I first submitted ./subset_data region with qsub -l walltime=03:00:00 -l select=4:ncpus=128qsub which did not look like it would complete in 3 hours. So, I then tried qsub -l walltime=12:00:00 -l select=24:ncpus=128. This job is currently on its 11th run hour. I will update when it reaches the end.

The weird thing is that I can see the subsetted output surface file. The variables seem to be complete; I don't have the issue of truncated %clay, %sand, and others anymore. For instance, output for PCT_CLAY variable:

But I can't view monthly_LAI and SAI. The output file is in /glade/work/amans/data/MRB-CONUS_0.05/surfdata_0.05-conus_hist_2000_78pfts_MRB_0.05_c240322.nc.

Can I expect to get complete output if I submit it for a longer time, say 24 hours?
Also, is there a way to track the progress? I would really appreciate it.

 

[amans]

I then tried qsub -l walltime=12:00:00 -l select=24:ncpus=128. This job is currently on its 11th run hour. I will update when it reaches the end.

PBS: job killed: walltime 43293 exceeded limit 43200

ncview still doesn't display the map for various variables like monthly_LAI and SAI. ncview terminates with this error:
calculating min and maxes for MONTHLY_LAI.netcdf_fi_get_data: error on nc_get_vara_float call
cdfid=65536 variable=MONTHLY_LAI
start, count:
[0]: 0 1
[1]: 0 79
[2]: 0 458
[3]: 0 767
NetCDF: Numeric conversion not representable

Another question I have is, is there a way to resume subset_data from the point it last stopped, i.e., at the end of the 12 hours in the case above?

 

[slevis]

Unfortunately there is no way to continue subset_data from where it stopped. I hope that subset_data will work for you if you ask for 24 hours instead of 12...

 

[slevis]

Did you try plain ./subset_data without the qsub command? If so, did that job run out of memory?

 

[amans]

Update: I recently submitted the subset workaround for 24 hours with qsub -l select=1:ncpus=128. It successfully completed but it consumed the whole 24 hours, equaling around 3000 core-hours. The subsetted surface has no truncated variables and is running for my case.
I did not try plain ./subset_data.

 

[amans]

Did you try plain ./subset_data without the qsub command? If so, did that job run out of memory?

Plain ./subset_data did not work. The job was killed after ~10 minutes and I got a high usage email probably due to huge memory usage.

 

[slevis]

Thank you for the update.