Running FATES with soil thickness

[jiamenglai]

Describe the changes

I tried to run FATES with soil thickness read from the surface data, however I met the error:

if ( use_var_soil_thick ) then
call endrun(msg=' ERROR: use_var_soil_thick and use_fates cannot both be set to true.'//&
errMsg(FILE, LINE))
end if

Then I tried the method proposed here #7319 by changing the ZSOI. It worked for site simulation. However, the variable ZSOI in the surfdata only has one dimension (nlevsoi), which means it cannot vary with spatial locations. I think it might be because that the soil depth in ELM code is defined as spatially constant. If I want to run regional simulations, how can we represent the spatially varying soil depth?

Is your request related to a problem?

No response

Alternative options

No response

Additional context

No response

[mpaiao]

@rgknox @glemieux Do you remember why ELM does not accept variable soil thickness together? It seems CTSM does not have such limitation. I wonder if this is a legacy code, or something we kept until someone decided to test it.

@jiamenglai as an option, maybe give it a try by commenting out this if block and see what happens. It may either work or explain why this barrier is in place...

[wwieder]

CLM does have this issue, but @glemieux has a PR that fixes the issue.

[glemieux]

@rgknox @glemieux Do you remember why ELM does not accept variable soil thickness together? It seems CTSM does not have such limitation. I wonder if this is a legacy code, or something we kept until someone decided to test it.

I'm not sure of the reason here. That particular logic is from about 5 years ago and I think is associated with a much older issue: E3SM-Project/E3SM#2159.

[jiamenglai]

@rgknox @glemieux Do you remember why ELM does not accept variable soil thickness together? It seems CTSM does not have such limitation. I wonder if this is a legacy code, or something we kept until someone decided to test it.

@jiamenglai as an option, maybe give it a try by commenting out this if block and see what happens. It may either work or explain why this barrier is in place...

Thanks for the suggestion. I tried to comment out this block and run the simulation. It gave an error in 'init_soil_depths' around the following lines:

if(abs(this%fates(nc)%bc_in(s)%dz_decomp_sisl(j)-this%fates(nc)%bc_in(s)%dz_sisl(j))>1.e-10_r8)then
write(iulog,) 'when vertical soil decomp dynamics are on'
write(iulog,) 'fates assumes that the decomposition depths equal the soil depths'
write(iulog,) 'layer: ',j
write(iulog,) 'dz_decomp_sisl(j): ',this%fates(nc)%bc_in(s)%dz_decomp_sisl(j)
write(iulog,*) 'dz_sisl(j): ',this%fates(nc)%bc_in(s)%dz_sisl(j)
call endrun(msg=errMsg(sourcefile, LINE))
end if

However, I checked the code in initVerticalMod.F90, the use of 'zbedrock' seems to only affect the value of nlevbed but would not alter the calculation of decomposition and soil depth?

[glemieux]

Per the SE meeting today, @glemieux will review this later this week.

[mpaiao]

I am not familiar with the decomposition code in ELM or FATES, but I checked the code @jiamenglai noted above (elm/src/main/elmfates_interfaceMod.F90) and wondered if there is a bug in there.

The code assigns dz_sisl and dz_decomp_sisl as follows:
https://github.com/E3SM-Project/E3SM/blob/5f825dbd95d57bea04b240b702eb4e69b00982be/components/elm/src/main/elmfates_interfaceMod.F90#L3514-L3518

So dz_decomp_sisl is set for layers 1 through nlevdecomp, whilst dz_sisl is set for layers 1 through nlevsoil. But then the code checks whether or not dz_sisl matches dz_decomp_sisl for layers 1 through nlevsoil.

https://github.com/E3SM-Project/E3SM/blob/5f825dbd95d57bea04b240b702eb4e69b00982be/components/elm/src/main/elmfates_interfaceMod.F90#L3521-L3533

Shouldn't it only check for layers 1 through nlevdecomp? Perhaps nlevdecomp is identical to nlevsoil for standard runs, but not when zbedrock is defined.