SP mode: failed with grass allometry when check tree_lai (time for grass allometry to be separated from trees?)

[XiulinGao]

So I did what @huitang-earth suggested (https://github.com/huitang-earth/fates/blob/367c8fa29b99a35e1f467c4f87ead30f1b96aff7/main/EDPftvarcon.F90#L1717) to temporarily pass the pft mapping check in order get the run continue. Then it crashed again with error messages as below:

in sub-canopy layer in leafc_from_treelai
this is not set up to work for lower canopy layers.
problem in leafc_from_treelai 79518032 1
ENDRUN:ERROR in /home/xiulingao/var-test-runs/CTSM/src/fates/biogeochem/FatesAllometryMod.F90 at line 811 $
ERROR: Unknown error submitted to shr_abort_abort.

The grass allomtery parameters I used are from @xuchongang , which are as below:

"fates_allom_agb1 1 0.000964"
"fates_allom_agb2 1 1.9492"
"fates_allom_agb3 1 0"
"fates_allom_agb4 1 0"
"fates_allom_agb_frac 1 1"
"fates_allom_d2bl1 1 0.000964"
"fates_allom_d2bl2 1 1.9492"
"fates_allom_d2bl3 1 0"
"fates_allom_d2ca_coefficient_max 1 0.03"
"fates_allom_d2ca_coefficient_min 1 0.01"
"fates_allom_d2h1 1 1"
"fates_allom_d2h2 1 1"
"fates_allom_dbh_maxheight 1 1"
"fates_allom_hmode 1 3"
"fates_allom_la_per_sa_int 1 1000"
"fates_allom_sai_scaler 1 0.0012"

According to Chonggang's explanation these are based off tillers of grasses not the whole bunch, which works for his grass species with runners(so it's hard to define individuals as one single parent plant can grow into a large clone). This allometry setting apparently caused problem in tree_lai checking as I ran the model by taking the allometry parameters out (so ran with default allometry parameters) and it works.

[rosiealice]

Hmm. I wonder if there is something up with the conversion from bleaf to treelai using those allometry options?

[XiulinGao]

I don't know if this has something to do with using tiller diameter instead of using stem(or trunk) diameter for grass here. @ckoven might know?

[XiulinGao]

update: I'm now not sure if this has something to do with the allometry parameters. I ran the simulation again with exact same settings(default allom params that I did run successefully last time), and this time it failed with the same error. I think under SP mode there's something going on with this chunk of code here:

fates/biogeochem/FatesAllometryMod.F90

Lines 800 to 812 in 70a3e13

	if( treelai < 0._r8.or. pft == 0 ) then
	write(fates_log(),*) 'negative tree lai in leafc_from_treelai?'
	write(fates_log(),*) 'or.. pft was zero?'
	write(fates_log(),*) 'problem in leafc_from_treelai',treelai,pft
	call endrun(msg=errMsg(sourcefile, __LINE__))
	endif
	if(cl>1)then
	write(fates_log(),*) 'in sub-canopy layer in leafc_from_treelai'
	write(fates_log(),*) 'this is not set up to work for lower canopy layers.'
	write(fates_log(),*) 'problem in leafc_from_treelai',cl,pft
	call endrun(msg=errMsg(sourcefile, __LINE__))
	endif

Is this specific to grass PFT or to all simulations with one PFT? need to figure out.

[rosiealice]

Did it fail with one of those errors in particular? (And if so, did it write the outputs??)

[XiulinGao]

It did fail with the line 807-811, that cl >1. the error message:
in sub-canopy layer in leafc_from_treelai
this is not set up to work for lower canopy layers.
problem in leafc_from_treelai 2 1
ENDRUN:ERROR in /home/xiulingao/var-test-runs/CTSM/src/fates/biogeochem/FatesAllometryMod.F90 at line 811 $
ERROR: Unknown error submitted to shr_abort_abort.

I also noticed that in the successed simulation, there's no understory at all so I assume cl = 1. What is this about SP mode? CL must be one? There are some huge numbers of CL in my previous failed runs.

[rosiealice]

Yes. In principle, all being well, SP mode should create a cohort that exactly fills one canopy layer. The trees will be the height dictated by the surface dataset, and a dbh and canopy area consistent with that. The idea is that this routine calculates, for those trees, the leaf carbon that is conistant with the input TLAI. If CL>1, then something has gone wrong in the calculation of canopy area.

[rosiealice]

Did you try an 'out of the box' run with SP mode? At what point does it start to fail?

[XiulinGao]

Yes, I'm running a single point simulation with only C3 grass. I also substitude MONETHLY_LAI in the surface data with MODIS observation for the site. This error occurs after case.submit
and I also noticed that sometimes it can run but most time it does not.

[rosiealice]

Does it work with the default PFT file? (to rule out other gremlins)
What about a single PFT that is a tree? ( to check if it's a grass-specific or single-PFT specific problem)
sorry aobut this!!!

[rosiealice]

Great. Let me know if you manage to isolate the cause of the failure in shifting to the grasses.

[XiulinGao]

No worries! Thank you for helping me out here! So I ran with single tree PFT using default param file, and ran it twice (so in case the success is just out of luck) and both runs successed! Then I ran with single grass PFT with default param file, ran twice, both simulations ran w/o problem. So I assume there's something wrong with my other param setting. I'll try to figure out which.

Hi Rosie, I'm sorry. I took back what I said. I accidently ran it 4 times with the fates_params_default.cdl, so those 4 runs do not count (guess at least we know it won't fail with default param file with 12 PFTs then). It failed with both single tree and grass PFT (no params changed). So it is a single PFT issue.

Same error BTW.

[rosiealice]

Righto. What was the fail message? The >1 CL ?

[XiulinGao]

yes, CL >1. I want to ask if my fates_hlm_pft_map is [0,0,0,0,0,0,1,1,0,0,0,0,0,0] but in surface data PCT_NAT_PFT is
[0.000000 2.696105 0.000000 0.000000 0.000000 0.000000 0.000000 6.051406 0.000000 0.000000 0.000000 0.000000 0.000000 91.252488 0.000000], there's still only one tree PFT as #9 in PCT_NAT_PFT is 0 right? if this is true then single tree can run. But not matter what there's the CL >1 error with single grass run. I think I have to make sure that the tree PFT I used to run the simulation actually has a percent cover on the site (refer to the PCT_NAT_PFT in surface data), otherwise the model will fail?

[XiulinGao]

update: I changed hlm_pft_map all values to 1 for single grass PFT trying to run SP mode as suggested by Rosie, it failed at the first try. Then I changed DEBUG to TRUE and it ran successfully. In case it was just that one time it did work(as it did before!) I changed DEBUG to FALSE and then TRUE again to run it, model failed when DEBUG off and did work when DEBUG on(and I did try DEBUG on for several times and each time it did run). I'm getting the same canopy layer >1 error message whenever model crashed. I don't understand why DEBUG on can make it run? I haven't looked at the output though.

[rosiealice]

OK I think I made some progress. What seems to be happening is an order of operations problem in init_patches.

Specifically, there is a check that makes sure that the patch areas add up to one exactly. This is occasionally not the case owing to (I think) precision errors in the reading in of the surface dataset.  At the moment there is a check for this that happens after the patches have all been created. The code then adjusts the area of the oldest patch, if there is sufficient area, (and the youngest patch if not) to account for the difference between the total area of the patches and of the column ("AREA" which is perhaps a variable that now deserves a more verbose name...) 
That happens here:

fates/main/EDInitMod.F90

Line 583 in 5d449e4

sites(s)%oldest_patch%area = sites(s)%oldest_patch%area - (tota-area)

The problem with this, as far as I can see, is that the init_cohorts command, which sets the c_area of the single SP cohort, is called BEFORE the change to the patch area in the line linked above. The routine that sets c_area (assign_sp_cohort_properties) is also called in run time, and thus the error is ultimately corrected, but I think we have to go through the checks that are currently failing at the end of initialisation too. 

To fix this, I made another instance of the code that adjusts the patches areas for precision errors, but this time it acts on the incoming PFT area vector (sites(s)%area_pft(nocomp_pft)), and not on the area of the resulting patches. The patch adjusted is arbitrarily the first one, which I guess introduces weighting in favour of earlier PFTs/bare ground, but we are talking 10^-12 order errors here so I'm not too worried about that. 

The code is on a branch here:

https://github.com/rosiealice/fates/commits/SP_onepft_bugfix

This code runs OK with a single PFT in my 4x5 global test when all of the hlm_pft_map values are set to 1.0, and using a single PFT file that Xiulin provided (thanks @XiulinGao!)

I haven't tested it with the original values of hlm_pft_map. Could you have a go with this @XiulinGao and see if it triggers the errors you were seeing?

Noting that

1. I didn't see any of the CL>1 error statements and
2. I don't know why this problem is triggered in the special case of 1 PFT and not with multiple PFTs. It might thus be that this problem masks another. We will see...

Cross fingers!