560 new calculation for accounting of root distribution on the transpiration #591
Conversation
ccarouge
force-pushed
the
560-New-calculation-for-accounting-of-root-distribution-on-the-transpiration
branch
from
6dda783 to
4d0093f
Compare
ccarouge
requested review from
har917,
a team and
SeanBryan51
and removed request for
a team
|
@ccarouge I'll look through in detail later - nothing major jumped out so far except
Nevertheless 'we' do need to resolve why the intel build and markdown tests have failed before this can go further. One question - how encompassing do you want this issue to go? We could easily extend it to include rationalising the evaluation of |
|
@ccarouge Could you add some plots of the water balance ( |
There was a problem hiding this comment.
@ccarouge Looks mostly okay from a science component:
- I do wonder whether we should extend this to cover consistency changes in
fwsoil - We do need to think how to deal with
diff(k=ms)/=0.0at the end ofremove_transshould be handled (quite likely given the inconsistencies infwsoil)
Technically - there's the broken markdown and build - I also suspect that MPI is broken since the KIND for ssnow%evapfbl has changed and so I would have expected that the r1len's need to go to r2len
| @@ -7222,13 +7208,6 @@ SUBROUTINE master_restart_types (comm, canopy, air, bgc) | |||
| ! & types(bidx), ierr) | |||
| ! blocks(bidx) = 1 | |||
|
|
|||
| bidx = bidx + 1 | |||
There was a problem hiding this comment.
This change here is consistent with the intent of the PR but it highlights something missing (I think)
Since you've changed the TYPE of ssnow%evapfbl you should also need to change the length of the associated MPI declarations (should be r2len not r1len) - I don't see these changes in the changeset.
There was a problem hiding this comment.
Just a check - does the blen(bidx) = r1len at line 2513 also need changing?
There was a problem hiding this comment.
I believe the line 2513 you are referring to has been removed. It was for canopy%evapfbl that I removed. Line 2513 now has the MPI code for canopy%epot.
There was a problem hiding this comment.
I can't find any reference to ssnow%evapfbl in mpimaster or mpiworker that is using r1len anymore.
There was a problem hiding this comment.
This call to
MPI_Type_Create_hvectorseems to be creating something (types(bidx)?) based on bidx being an r1len variable at this point in the sequence - but that's no longer valid.
I think this must have been fine since the kind was not specified for canopy%evapfbl and defaulted to r1.
There was a problem hiding this comment.
@SeanBryan51 But would it be fine with the PR's changes given that bidx is pointing (not in the sense of a POINTER) to ssnow%evapfbl which is now r_2?
There was a problem hiding this comment.
@har917 , bidx is only used as an array index. types is an array and we add the information for evapfbl in vector form to that array (my guess is that it stores some kind of total length in bytes). The index used for types does not change depending on the type of the variable (r_2 or r_1). types contains INTEGERs, so I'd assume the integer for a r_2 array is bigger than for a r_1 array but it's still only 1 number that only needs to be stored in one place in the types array. So bidx (the index counter for the types array) is still incrementing by one. And we only need one element of the types array types(bidx) to store the information.
There was a problem hiding this comment.
My concern isn't really bidx or types but blen(bidx) - this seems (e.g. line 5659) to be used when setting up how the workers-masters will communicate. So setting this up to have the length 1 or r1len when %evapfbl is an r2len double precision REAL seems ... odd.
However - some of the variables are clearly being dealt with differently - see line 1983 for veg%froot which also does the call to MPI_Type_Create_hvector then blen(bidx)=1 - so there is some consistency (even if I don't follow the why).
Perhaps what I'm really getting at here is that this probably needs a short MPI test of just the change in KIND on ssnow%evapfbl separated out from the rest of the science change (? - more work though).
There was a problem hiding this comment.
@SeanBryan51 But would it be fine with the PR's changes given that bidx is pointing (not in the sense of a POINTER) to
ssnow%evapfblwhich is now r_2?
Yes I think so as long as worker_cable_params is consistent with how the blocks have been organised in master_cable_params (i.e. blen and bidx) there should be no issues, and this looks to be the case.
| SUM(veg%froot * MAX(0.0,MIN(1.0, REAL(ssnow%wb) - & | ||
| SPREAD(soil%swilt, 2, ms))),2) /(soil%sfc-soil%swilt)) | ||
| SUM(veg%froot * MAX(0.0,MIN(1.0, REAL((ssnow%wbliq) - & | ||
| SPREAD(soil%swilt, 2, ms)))),2) /(soil%sfc-soil%swilt)) |
There was a problem hiding this comment.
Not sure I like the use of SPREAD(soil%swilt,2,ms) here - or %sfc- %swilt in the denominator. Seems inconsistent with the push to use either normal or _vec forms consistently ...
... but okay scientifically.
There was a problem hiding this comment.
For fwsoil calculations, I've only changed wb to wbliq and didn't change the use of the parameters. So what you see here for swilt and sfc is what is already in main.
I feel we will need to decide on one form of the soil parameters (2D or 1D) and change everything to that. Because I was changing remove_trans() code, I introduced the _vec forms to avoid code duplications. But I would generally hold off in expanding this unless it greatly simplifies the code, until we decide which form we keep.
There was a problem hiding this comment.
Sorry, forget what I said above @har917 , I didn't check on the context of the code.
The reason I didn't use the _vec form here is I don't know how the non-linear calculation is supposed to work with a vertical profile in soil properties.
I've noticed the groundwater code does not solve the problem.
It could be we need to calculate a rwater term for each soil layer and sum at the end?
Anyway, this is beyond the current PR, I believe.
There was a problem hiding this comment.
In the first instance I would be going with
rwater = MAX(1.0e-9, &
SUM(veg%froot * MAX(0.0,MIN(1.0, (REAL(ssnow%wbliq) - &
soil%swilt_vec)/(soil%sfc_vec-soil%swilt_vec) )) , 2)
which I think is the direct equivalent once you substitute in the _vec forms filled with the single layer values.
I was 'simply' trying to avoid this need to carry both single and _vec forms of the soil parameters (and we really shouldn't need the REAL() here either - or perhaps we need it but with an explicit change to a KIND(r1))
There was a problem hiding this comment.
My problem is the rest of the algorithm:
CABLE/src/science/canopy/cbl_fwsoil.F90
Lines 58 to 80 in 7e83ba7
The way it is written relies on swilt and sfc being scalars. If they vary in depth, then my guess is we need to add a loop on the soil layers and sum at the end. But I don't know if the ti values should be kept as 0, 0.9 and 1 in that case.
I certainly don't want to introduce the _vec form for only half of the function and the scalar form for the rest.
Co-authored-by: Sean Bryan <39685865+SeanBryan51@users.noreply.github.com>
There was a problem hiding this comment.
@ccarouge I'm happy with the science - and @SeanBryan51 is okay with the rest of the technical changes (barring coupled model stuff).
|
Thanks for approving. Unfortunately, I believe this has to wait until we get the library build ticked off. At least, I need to talk to Lachlan before merging in code changes that modify the results! |
CABLE
Thank you for submitting a pull request to the CABLE Project.
Description
Fixes #560
A few changes come with this pull request:
cbl_dryLeaf.F90andcbl_remove_trans.F90.cbl_dryLeaf.F90and stored inssnow%evapfblto update the soil water content incbl_remove_trans.F90instead of recalculating it. It's simplifies the code and has a very minor impact.canopy%evapfbl, only keepssnow%evapfbland give it a double precision. The precision change is in the groundwater updates that has started this work.wbliqinstead ofwbso the calculation only takes into account the liquid water amount. EDIT: Changed the calculations offwsoilto usewbliqconsistently.zse_vecandswilt_vecto make the same calculation work for standard and groundwater. It introduces the use of the *_vec variables in the main part of the code. The alternative is to create copies of 2 subroutines and some interfaces to handle the different dimensions and precisions. It seems overly complex for something that small. I'm also thinking we will want to carry only one version of the soil properties at some point. However, that's open for discussion.zse_vecwasn't properly initialised in CABLE and I have corrected that. I have introduced the initialisation to ESM1.6 and it is already correct in AM3.swilt_vecis initialised correctly in all current applications. If we go with this version of the code, I'll test it works in the coupled applications as well.I also looked at
getrex_1dused in the Haverd2013 case. There are 2 versions in CABLE (cbl_dryLeaf.F90 and cable_sli_roots.F90). Unfortunately, they have slightly diverged from each other. For the moment, I have simply added a warning to the code documentation with a list of the differences I have found. But I didn't want to make this PR bigger than required. I'll leave that part to another time.Testing
modelevaluation results
In addition to the benchcab testing, I've also tested the effect of various steps.
Changing to use
wbliqhas the largest effect for points with icy fraction, differences at Hyytalia for all soil layers:EDIT: Figure updated with the fwsoil calculation changes.
Harmonising the calculations between the canopy and soil has the second largest effect, differences at Tumbarumba for the top soil layer:
Please add a reviewer when ready for review.
📚 Documentation preview 📚: https://cable--591.org.readthedocs.build/en/591/