issue ed-hydro with specific allometry parameters

[sdeherto]

Hi
I have been facing an error when running ED2 with plant hydraulics on that occurs in plant_hydro.f90 which I would like some input on from people who know this part of the code (@xiangtaoxu )?

I have recently created a new IALLOM setting that uses the dbh-height relationship from IALLOM 2 and all other allometric relationships from IALLOM 3 (see my code here:https://github.com/sdeherto/ED2 ). I also provide some parameter values related to the height allometryt for my site via the config files for three custom made PFTs for my site. My runs crash during the first timestep when entering the plant_hydraulics with the message

===========================================================================
Invalid leaf_psi or wood_psi detected.

Note that this implementation works for ED2 when plant hydrualics are not turned on.
So it is a rather specific issue as I noticed that it only occurs for certain values of my parameters, e.g. my pft 4 has values for which it crashes (b1Ht: 0.0575, b2Ht=0.644, hgt_ref=55), if I set these parameter values equal to the other PFTs (i.e. b1Ht: 0.0178, b2Ht=1.27, hgt_ref=36.8) it runs without issue.

Any thoughts on this would be much appreciated! I have been looking through this with Marcos Longo a few times and there doesn't appear to be immediate issues in the iallom implementation I did. I should also note that when printing these values out at the start of the plant_hydro code it already shows positive wood_psi values while when I print the same variables out at the start of the initialisation code all values still look fine. Does anyone have an idea where wood_psi could get changed do that it becomes positive?

Cheers
Steven

[mpaiao]

@sdeherto I was looking into the code, and found one bug in subroutine tw2rwc (file ED/src/dynamics/plant_hydro.f90), likely introduced by me through hastily copying and pasting code... I am not sure this would cause the crash because it is an edge case, but this may have the potential to change the value of wood_rwc.

ED2/ED/src/dynamics/plant_hydro.f90

Lines 1395 to 1401 in bae4504

	if (tot_water_sat > tiny_num) then
	wood_rwc = wood_water_int / tot_water_sat
	elseif (is_small) then
	wood_rwc = op_buff * small_rwc_min(ipft)
	elseif (is_small) then
	wood_rwc = op_buff * wood_rwc_min(ipft)
	end if

Based on the logic for leaf, I believe this should be:

      if (tot_water_sat > tiny_num) then
         wood_rwc = wood_water_int / tot_water_sat
      elseif (is_small) then
         wood_rwc = op_buff * small_rwc_min(ipft)
      else
         wood_rwc = op_buff * wood_rwc_min(ipft)
      end if

If this resolves the crash, then the next question is figuring out why your tot_water_sat is so low...

[xiangtaoxu]

@sdeherto , can you post/upload the full error message. Knowing the error flags would be helpful. And if it crashes in the first time-step, which means the error might be induced by initialization (as Marcos suggested), it could be helpful to run the model in debug mode (KIND=B in compilation)

[sdeherto]

Hi @xiangtaoxu and @mpaiao

Thanks for the quick responses! I quickly tried the edits you suggested Marcos but it didn't resolve the error. The error message is attached to this post (note that a lot of values are printed this is from checking some values such as wood_psi in the initialisation script and the plant hydraulics scripts).

Thanks for the suggestion on using debug compilers, I will look into this!

test.txt

Cheers
Steven

[xiangtaoxu]

@sdeherto , the error message suggests wood_psi becomes positive without any water flux. This is likely an error happening due to conversion between rwc, tw (total water), and psi. I don't have a good sense right now where the problem is but I did have (bad) experience when playing with iallom, which has various cascading effect in ed_params...

Since the water potential is already anomalous before calling plant_hydro, my bet would be to check conversion functions such as rwc2psi, tw2rwc and their counter parts.

Can you try print out wood_psi and wood_rwc for the problematic cohort before rwc2psi call in this block of codes in ed_type_init?

https://github.com/sdeherto/ED2/blob/75c566bf8eadc3e75f99202dd2ef640cb291b75a/ED/src/init/ed_type_init.f90#L239-L249