-
Notifications
You must be signed in to change notification settings - Fork 318
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Ctsm sectorwater pull request v2 #1975
base: master
Are you sure you want to change the base?
Ctsm sectorwater pull request v2 #1975
Conversation
…and continue pointing at the ESCOMP repos, but then code will work once the pull request is finished for these 3 dependency repos.
…ut data through mksurfdata tool)
…the consumption flux from sectoral water usage do not cause any errors.
…sumption. Operation performed at column level.
…ould be disposed on surface soil. The sectoral consumption is disposed on natural vegetation soil columns to avoid interference with irrigated columns.
…er competition principles. Irrigation is last in priority, therefore max irrigation withdrawal is limited by not only taking into account VOLR, but also how much was already taken by other sectors.
…routing model. In this commit the required support is added to the CTSM coupler part.
…, only one subroutine call is required. The call is done before the irrigation, as other sectors are higher in their priority. The CalcIrrigationNeeded() is also updated to take into account how much water was withdrawn for other sectors in given day.
…le with sectoral water abstraction for year 2000.
Thanks for opening up this PR, Sabin. This is exciting work. Have you been working with or communicating with anyone else in the LMWG on this work that we should bring into the conversation? Regardless, we'll have a look at your PR and start making a plan to move forward. |
After further discussion we'd like to bring in the CTSM side of this PR with namelist changes so that users have to turn on this feature under limited applications. @samsrabin is likely going to have time to work on this after doing some additional crop model developments (maybe by the end of this summer). The code should also have a science review, maybe @swensosc can have a look on this (also not urgent given other priorities)? |
Has any testing been done? Perhaps to accelerate things, Sabin could run
the short test suite. That would catch it if there is a big problem (i.e.,
that when turned off it is not bit-for-bit. In general, it would be good
to make it possible for external collaborators to easily run the simple
test suite, which could help things move along faster. I know we have
talked about this before and some people have done it, I think.
…On Thu, May 25, 2023 at 10:24 AM will wieder ***@***.***> wrote:
After further discussion we'd like to bring in the CTSM side of this PR
with namelist changes so that users have to turn on this feature under
limited applications. @samsrabin <https://github.com/samsrabin> is likely
going to have time to work on this after doing some additional crop model
developments (maybe by the end of this summer). The code should also have a
science review, maybe @swensosc <https://github.com/swensosc> can have a
look on this (also not urgent given other priorities)?
—
Reply to this email directly, view it on GitHub
<#1975 (comment)>, or
unsubscribe
<https://github.com/notifications/unsubscribe-auth/AFABYVEGTBCP3WLO6XZSJLLXH6BUXANCNFSM6AAAAAAWLWRPFI>
.
You are receiving this because you were mentioned.Message ID:
***@***.***>
|
@wwieder, thanks for the update. As mentioned, we can provide any required assistance in the process. @dlawrenncar, I did some tests (a complete list can be found in the description of this pull request #1975). My intention was to compile all the testing and validation into a series of Jupyter notebooks, which will be published with the development paper. Please, let me know if there is some additional testing, which is required. |
I am referring to our systems tests which we use to make sure things like
exact restarts work and that there aren't any divide by zero possible
situations, etc. These are not science tests, but software tests. The
full suite of tests runs into the hundreds of tests, but there is a smaller
one that catches some of the most common errors.
…On Wed, May 31, 2023 at 2:03 AM Sabin ***@***.***> wrote:
@wwieder <https://github.com/wwieder>, thanks for the update. As
mentioned, we can provide any required assistance in the process.
@dlawrenncar <https://github.com/dlawrenncar>, I did some tests (a
complete list can be found in the description of this pull request #1975
<#1975>). My intention was to compile
all the testing and validation into a series of Jupyter notebooks, which
will be published with the development paper. Please, let me know if there
is some additional testing, which is required.
—
Reply to this email directly, view it on GitHub
<#1975 (comment)>, or
unsubscribe
<https://github.com/notifications/unsubscribe-auth/AFABYVDMJG4OPW2SY3JJB5LXI33LHANCNFSM6AAAAAAWLWRPFI>
.
You are receiving this because you were mentioned.Message ID:
***@***.***>
|
Thanks for the extensive testing you've done already, @TaranuDev ! Dave is right that a bit more using the built-in testing tools would be nice to try and identify any hidden bugs. Fortunately, it's pretty simple to do a few quick tests using the Some notes first:
First, run the test suite with the version of CTSM you branched from (or the most recent version you merged in), which seems to be
The tests will be set up and run from the directory given in the first line of the output of that command. Once they've started, you can check their status by going there and doing The tests whose names start with When that's done, you can re-run the tests with your code like so, which will also compare the outputs bit-for-bit with the previous run:
You can find more details on any failed comparisons using the instructions here. Try looking for output variables with large normalized root mean square (RMS) errors, as these are the ones that see the largest relative difference between runs. |
@dlawrenncar, thanks for the clarifications, we didn't run these tests yet, but we will do it. @samsrabin, thank you for the detailed instructions, it is very helpful! I will run these tests in the near future. If I have any issues, as suggested, I will start a discussion in the dedicated section. |
Just an update about the timeline here: The plan at this point is for me to work on integrating this in the late summer. @TaranuDev, if you could get those tests done in the next couple of months that'd be excellent. It's not a big deal for me to do the tests myself, but if they reveal any issues it'd be best for you to get a head start on fixing them. |
Thanks for the update @samsrabin! I planned to do the tests by the end of the next week. Cheers, |
Quick update on the requested testing: I did the run_sys_test. I expected the NLCOMP to fail since I add new namelist variables. At the moment, I did a quick check of the log files for the ERP_D case, and I found many non-zero differences. In principle, it is possible that there is a domino effect, where an error in one variable can propagate to many others.
For now, I am more inclined to believe that it is the option (2), because when I checked why the NLCOMP failed for ERP_D case, there are a lot of differences unrelated to my development (e.g. a different topography file, extra fields related to O3, etc). I am giving this incomplete update because I promised to finish the testing last week, and experienced a delay due to work on a different project. I will be on vacation the following two weeks, but when I will be back I will redo the testing and make sure that everything works perfectly. Sorry for the delays! Best regards, P.S. for the record here is the location of my tests: |
Hi Sabin, Thanks for doing those tests! This might indeed be a result of choosing the wrong baseline. From your branch's commit history, it looks like you branched from commit f72479b, which is actually tag |
I redid the tests with the correct baseline (ctsm5.1.dev119). For the NLCOMP, the test is failing because: "Found extra variable: 'sectorwater'; Found extra namelist: sectorwater_inparm". For the BASELINE fails, I get the message: ctsm5.1.dev119_clmshort_baseline: FIELDLIST field lists differ (otherwise bit-for-bit). I get the same result for each experiment in the test suit and compiler. The tests results can be found here: /glade/scratch/staranu/tests_0712-040306ch
Please let me know if anything else is required. I will try to finish the paper on this development this August, so that everything, including testing and validation, is well documented. Directly after, I will work on extending the model functionality for future scenarios too, as we previously discussed. Best regards, |
That's great, thanks Sabin! There are more tests to do, but I can handle those. Best wishes for the paper-writing phase; I'm happy to look over a complete draft if you'd like, but totally not necessary. |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
You've done a lot of great work here, Sabin. I've added a few questions and suggestions at various points in the code. In addition, some general questions:
- (1) Sector water input data should probably be provided as a single netCDF file instead of providing
path_sectorwater_input_data
, which is a list of paths with one netCDF file for each year. Unless there's a reason I'm not aware of? - (2) You say that "The consumption flux (considered water lost in the process of usage due to various factors) is disposed on the surface of soil columns with natural vegetation." Could you clarify? "Consumption" to me implies everything that is used; i.e., withdrawals minus conveyance losses (e.g., leaks from pipes). But here it sounds like "consumption" is just anything that's withdrawn for non-irrigation sectors and ends up back on the ground. So, e.g., water evaporated in industrial processes wouldn't be included here. If that's the case, where is that water included?
- (3) Dumping all "consumed" water on the natural column will affect plants and soils. That's true for some consumed water, but in reality some goes directly into rivers (or is evaporated, as discussed above). Is there any way around this? Could dump it all into ocean, but even that wouldn't account for evaporated "consumption."
- (4) What is the reason for passing both withdrawal and return flow from the land to the rivers? Could these be combined into one "net withdrawal" value?
- (5) Could all sectors be combined when passing from land to rivers?
Finally, it would be best if we could split this into two PRs: One that's not dependent on externals updates (MOSART, CMEPS, and CPL7) and one that is. However, we can think about how best (and whether!) to do that once we address everything else.
if (limit_sectorwater_if_rof_enabled) then | ||
if (sectorwater_river_volume_threshold < 0._r8 .or. sectorwater_river_volume_threshold > 1._r8) then | ||
write(iulog,*) ' ERROR: sectorwater_river_volume_threshold must be between 0 and 1' | ||
write(iulog,*) ' sectorwater_river_volume_threshold = ', sectorwater_river_volume_threshold | ||
call endrun(msg=' ERROR: sectorwater_river_volume_threshold must be between 0 and 1 ' // & | ||
errMsg(sourcefile, __LINE__)) | ||
end if | ||
end if |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
- This check of
sectorwater_river_volume_threshold
should probably happen inCLMBuildNamelist.pm
so that the model can't even be submitted with an invalid value. This subroutine could then be deleted. (It actually looks like it's never called?)
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I will do this.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I'll do this.
integer :: day ! day of month (1, ..., 31) for nstep+1 | ||
integer :: sec ! seconds into current date for nstep+1 | ||
|
||
integer :: first_read ! variable to do first read, in future I may prefer to make a subroutine is_beg_curr_month to avoid exception for first reading |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
- Should
first_read
be a logical instead? (Also, "in future" part of comment should be removed.)
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Sabin will do this.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Sabin will do this.
! Get current date | ||
call get_curr_date(year, mon, day, sec) | ||
dayspyr = get_curr_days_per_year() | ||
dayspm = dayspyr/12._r8 |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
- Do the monthly input files not account for the length of the month? If they do, we should use the actual number of days in the month (28–31).
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I'll work on this, but it doesn't make a big difference in results.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I'll do this.
end if | ||
|
||
! Convert demand to withdrawal rates [mm/s] | ||
! Here it also seems like I could directly operate with the this%arrays instead of generating new ones (to check) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
- Yes, I'm pretty sure you can. This would simplify this subroutine and its dummy variables significantly; let's give it a try. (Same thing for
CalcSectorDemandVolrLimited()
.)
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I will do this.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I'll do this.
! I think the algorithm is potentially too conservative | ||
! I will need to check river discharge when there is high amount of unsatisfied demand | ||
! The reason why I am saying this is because we compare the expected demand for the entire day with the current volume available in the river... | ||
! It would make sense to use such an algorithm under the condition that volr does not change much during a day for a given gridcell | ||
! But if this is not the case we maybe underestimate the amount of water available for usage. | ||
! If this would be done once a day, then no problem, but we do it at each time step | ||
! This means that volr get |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
- This temporal mismatch could also cause an overestimate of how much water is available, no?
- Is this how it's done for irrigation?
- Incomplete sentence at line 1163.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Conservative because it's looking at the river once and saying "this is how much water we have." But 3 hours later (timescale of coupling) that river will have flowed some more. So it's basically underestimating the available water by a factor of 8.
Seemingly this is how it's done for irrigation, but it's not as much of an issue there because irrigation is only applied over a 4-hour window (as opposed to 24 for other sectors' demand).
Calculating availability limitation at the 3-hour timestep seems tricky. For now, to see how much difference this might make,Sabin will test multiplying available water in the timestep by 8. Conservation errors will be handled, as for irrigation, by taking extra needed water from ocean.
@swensosc, do you have any thoughts on the best way to handle this?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Conservative because it's looking at the river once and saying "this is how much water we have." But 3 hours later (timescale of coupling) that river will have flowed some more. So it's basically underestimating the available water by a factor of 8.
Seemingly this is how it's done for irrigation, but it's not as much of an issue there because irrigation is only applied over a 4-hour window (as opposed to 24 for other sectors' demand).
Calculating availability limitation at the 3-hour timestep seems tricky. For now, to see how much difference this might make,
- Sabin will test multiplying available water in the timestep by 8. (Conservation errors will be handled, as for irrigation, by taking extra needed water from ocean.)
@swensosc, do you have any thoughts on the best way to handle this?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
@swensosc The factor of 8 is because we're looking at how much water is available in one 3-hour timestep, but demand is actually spread out over 24 hours. So 24/3 = 8.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Just chiming into make sure the final solution we use allows for the coupling timestep from the ROF model to vary. For a case this is the ROF_NCPL and NCPL_BASE_PERIOD variables. This isn't normally sent to CTSM, as ROF is expected to do it's own thing and usually CTSM doesn't need to know. If CTSM needs to know the coupling, we might need to pass it into the CTSM namelist.
! This consumption flux will be applied on surface soil of the natural vegetation column | ||
! The user should make sure that the water usage input data have no NaN values (and use 0 instead) | ||
! This way the sum of the sectors do not risk to become NaN, because one of the sectors have no data for given location | ||
! Of course this kind of problem can be overcomed in the code directly, but there is no good reason why it shouldn't be done at the level of input data |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
- NaN values should be checked for and handled (either set to 0 or error) during netCDF import.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Sabin will do this.
! Here I am not sure if it is needed to have loop over all tracers (it seems that the tracers mechanism is not maintained anymore) | ||
! So I could in principle just use associate : w => water_inst%bulk_and_tracers(1) which correspond to bulk water. | ||
! To stay in aggreement with the legacy code, I am looping over the tracers. |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
- I think it's right to loop over the tracers, so this comment can be deleted.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Sabin will do this.
if (sectorwater .AND. volr(g) > sectorwater_total_actual_withd(g)) then | ||
available_volr = (volr(g) - sectorwater_total_actual_withd(g)) * (1._r8 - this%params%irrig_river_volume_threshold) | ||
max_deficit_supported_by_volr = available_volr / grc%area(g) * m3_over_km2_to_mm | ||
elseif (sectorwater .AND. (volr(g) < sectorwater_total_actual_withd(g) .OR. volr(g) <= 0._r8)) then |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
- Is it possible for
volr(g) < sectorwater_total_actual_withd(g)
to not be satisfied whenvolr(g) <= 0._r8
is? If not, can delete the latter. In fact, it might be possible to rearrange the wholeif
block like this (requires thatsectorwater_total_actual_withd
is never negative, and that it's always 0 whensectorwater
is false):
if (volr(g) - sectorwater_total_actual_withd(g) > 0._r8) then
available_volr = (volr(g) - sectorwater_total_actual_withd(g)) * (1._r8 - this%params%irrig_river_volume_threshold)
max_deficit_supported_by_volr = available_volr / grc%area(g) * m3_over_km2_to_mm
else
! Ensure that negative volr is treated the same as 0 volr
max_deficit_supported_by_volr = 0._r8
end if
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Sabin will write a boolean table to see if these are logically equivalent; this simplification might be possible.
@@ -36,6 +36,17 @@ module Waterlnd2atmType | |||
real(r8), pointer :: qflx_rofice_grc (:) ! rof ice forcing, grc level | |||
real(r8), pointer :: qflx_liq_from_ice_col(:) ! liquid runoff from converted ice runoff | |||
real(r8), pointer :: qirrig_grc (:) ! irrigation flux | |||
real(r8), pointer :: qdom_withd_grc (:) ! sfc domestic actual withdrawal as satisfied by CLM |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
- Does
sfc
here mean "surface"?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I'll delete these.
call this%InitAllocate(bounds) | ||
call this%InitHistory(bounds) | ||
call this%InitCold(bounds) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
- To save memory, we should avoid doing all this when
sectorwater
is false.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I'll do this as one of the last things before integration.
I don't understand where the factor of 8 comes from?
…On Fri, Sep 8, 2023 at 9:35 AM Sam Rabin ***@***.***> wrote:
***@***.**** commented on this pull request.
------------------------------
In src/biogeophys/SectorWaterMod.F90
<#1975 (comment)>:
> +
+ liv_demand_volr_limited(bounds%begg:bounds%endg) = liv_demand(bounds%begg:bounds%endg)
+ liv_consumption_volr_limited(bounds%begg:bounds%endg) = liv_consumption(bounds%begg:bounds%endg)
+
+ elec_demand_volr_limited(bounds%begg:bounds%endg) = elec_demand(bounds%begg:bounds%endg)
+ elec_consumption_volr_limited(bounds%begg:bounds%endg) = elec_consumption(bounds%begg:bounds%endg)
+
+ mfc_demand_volr_limited(bounds%begg:bounds%endg) = mfc_demand(bounds%begg:bounds%endg)
+ mfc_consumption_volr_limited(bounds%begg:bounds%endg) = mfc_consumption(bounds%begg:bounds%endg)
+
+ min_demand_volr_limited(bounds%begg:bounds%endg) = min_demand(bounds%begg:bounds%endg)
+ min_consumption_volr_limited(bounds%begg:bounds%endg) = min_consumption(bounds%begg:bounds%endg)
+ end if
+
+ ! Convert demand to withdrawal rates [mm/s]
+ ! Here it also seems like I could directly operate with the this%arrays instead of generating new ones (to check)
I will do this.
------------------------------
In src/biogeophys/SectorWaterMod.F90
<#1975 (comment)>:
> + ! I think the algorithm is potentially too conservative
+ ! I will need to check river discharge when there is high amount of unsatisfied demand
+ ! The reason why I am saying this is because we compare the expected demand for the entire day with the current volume available in the river...
+ ! It would make sense to use such an algorithm under the condition that volr does not change much during a day for a given gridcell
+ ! But if this is not the case we maybe underestimate the amount of water available for usage.
+ ! If this would be done once a day, then no problem, but we do it at each time step
+ ! This means that volr get
Conservative because it's looking at the river once and saying "this is
how much water we have." But 3 hours later (timescale of coupling) that
river will have flowed some more. So it's basically underestimating the
available water by a factor of 8.
Seemingly this is how it's done for irrigation, but it's not as much of an
issue there because irrigation is only applied over a 4-hour window (as
opposed to 24 for other sectors' demand).
Calculating availability limitation at the 3-hour timestep seems tricky.
For now, to see how much difference this might make,Sabin will test
multiplying available water in the timestep by 8. Conservation errors will
be handled, as for irrigation, by taking extra needed water from ocean.
@swensosc <https://github.com/swensosc>, do you have any thoughts on the
best way to handle this?
------------------------------
In src/biogeophys/SectorWaterMod.F90
<#1975 (comment)>:
> + ! I think the algorithm is potentially too conservative
+ ! I will need to check river discharge when there is high amount of unsatisfied demand
+ ! The reason why I am saying this is because we compare the expected demand for the entire day with the current volume available in the river...
+ ! It would make sense to use such an algorithm under the condition that volr does not change much during a day for a given gridcell
+ ! But if this is not the case we maybe underestimate the amount of water available for usage.
+ ! If this would be done once a day, then no problem, but we do it at each time step
+ ! This means that volr get
Conservative because it's looking at the river once and saying "this is
how much water we have." But 3 hours later (timescale of coupling) that
river will have flowed some more. So it's basically underestimating the
available water by a factor of 8.
Seemingly this is how it's done for irrigation, but it's not as much of an
issue there because irrigation is only applied over a 4-hour window (as
opposed to 24 for other sectors' demand).
Calculating availability limitation at the 3-hour timestep seems tricky.
For now, to see how much difference this might make,Sabin will test
multiplying available water in the timestep by 8. Conservation errors will
be handled, as for irrigation, by taking extra needed water from ocean.
@swensosc <https://github.com/swensosc>, do you have any thoughts on the
best way to handle this?
—
Reply to this email directly, view it on GitHub
<#1975 (comment)>, or
unsubscribe
<https://github.com/notifications/unsubscribe-auth/AGRN57A6RT5YVGBWNPY65S3XZM3NLANCNFSM6AAAAAAWLWRPFI>
.
You are receiving this because you were mentioned.Message ID:
***@***.***>
|
Some notes from our very productive meeting this morning:
Sabin modeled this after how
All consumption is added to precipitation. The thought here was to minimize how directly we mess with coupled components. Sending to the ocean would likely cause similar complaints about how we withdraw from the ocean for unmet irrigation demand in the existing model; sending everything to atmosphere might mess with energy fluxes too much. For what it's worth, these fluxes are pretty minimal—on average ≤ 50 mm/year almost everywhere. I think I'd like to have a discussion with some others about whether there's a better way to do this, but for now we'll leave it.
Eventually, we may want to add some sort of pollution tracers to the return flows (e.g., sediment from mining operations), so we'll keep these separated.
This should be possible, yes. Of course, we'll still want to have them separate within CTSM so we can do the demand competition. The plan moving forward is that we can do one PR now without the changes to how the input data files work (or two PRs, first independent of and then dependent on externals updates). Then, when Sabin is done with the input data file rework, we can do another PR to bring that in. I will also:
|
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
@TaranuDev this is great, glad you are working on this. I noticed an easy issue looking at the discussion with @samsrabin, that I point out in my review. I did a quick look over and didn't go into detail. But, the thing that popped out was that 86400 needs to be a parameter. This is important for code readability, and also if any of the constants change, or need to be coordinated between models. And there are people that use CESM to run different planets so although on earth that doesn't change it is different there. Of course I wouldn't expect them to use sector water -- but maybe if you want to model Mar's colonies? Anyway, just something that helps with code readability.
I do like how you handled that for a different units conversion in m3_over_km2_to_mm. That makes the conversion very clear. That actually might be something we should use elsewhere in the code.
Thanks for the contribution.
! Add output variables | ||
! Domestic sector: | ||
this%input_mon_dom_withd_grc(begg:endg) = 0 | ||
call hist_addfld1d (fname='INPUT_MON_DOM_WITHD', units='mm', & |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
All of the history fields are marked as inactive. Which just means they need to be added to the namelist for them to appear on history. I'm just wondering if when sector water is active -- if some of them should be on? Sector water won't be the default at first, so turning on a few makes sense so you don't have to turn them on by hand. At the same time we want to try to pare down the number of active fields, as we already have a ton of output. Any fields that are 3D we need to be especially conservative with, but things that are 2D less so.
dayspm = dayspyr/12._r8 | ||
first_read = 1 | ||
|
||
flux_transform_from_monthly_to_second = (1._r8/dayspm)/86400._r8 |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
The 86400 here and elsewhere in this file should use the secspday parameter from clm_varcon.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Agreed! This is something I have noted in my "minor fixes I'll do" that I didn't bother putting in my review.
! I think the algorithm is potentially too conservative | ||
! I will need to check river discharge when there is high amount of unsatisfied demand | ||
! The reason why I am saying this is because we compare the expected demand for the entire day with the current volume available in the river... | ||
! It would make sense to use such an algorithm under the condition that volr does not change much during a day for a given gridcell | ||
! But if this is not the case we maybe underestimate the amount of water available for usage. | ||
! If this would be done once a day, then no problem, but we do it at each time step | ||
! This means that volr get |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Just chiming into make sure the final solution we use allows for the coupling timestep from the ROF model to vary. For a case this is the ROF_NCPL and NCPL_BASE_PERIOD variables. This isn't normally sent to CTSM, as ROF is expected to do it's own thing and usually CTSM doesn't need to know. If CTSM needs to know the coupling, we might need to pass it into the CTSM namelist.
! But if this is not the case we maybe underestimate the amount of water available for usage. | ||
! If this would be done once a day, then no problem, but we do it at each time step | ||
! This means that volr get | ||
if (dom_demand(g) * 86400.0 > max_demand_supported_by_volr) then |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
More 86400 here and below.
This is a work in progress, so converting to draft. |
Some potential issue raised by one of the reviewers for the corresponding paper: -> In the manuscript we mention Potential problem: Since we distribute the consumption over natural vegetation land unit, if we go very high resolution, and there will be a grid cell with a single land unit (e.g. urban or cropland), what will happen? (we might add to the code an exception for this case to avoid crashes). |
Description of changes
In this model development, we add support for sectoral water abstractions in CESM.
In many ways, this development is similar to the one done for irrigation abstractions already.
The withdrawal, consumption and return flow fluxes are computed daily in the land component based on monthly input data.
The information about the fluxes is communicated through the coupler to the routing model MOSART, from where the requested water is supplied. At 1 day cycle, withdrawal and recycling of water happens at the same time. For example, if 1km3 is abstracted and there is a local 80% recycling rate, then 0.8km3 is returning to the river network in the same day.
The consumption flux (considered water lost in the process of usage due to various factors) is disposed on the surface of soil columns with natural vegetation. From there this water, will participate in the surface water-energy balance through evaporation, runoff, infiltration...
The supported sectors are domestic, livestock, thermoelectric, manufacturing and mining.
Together with irrigation, this cover most of the known human water usage across the globe.
In addition to expected fluxes, if the user activates through the namelist, to limit sectoral abstraction based on river water availability, the model computes actual water withdrawals. In this case, when water is limited, a sectoral competition is happening between the sectors. At the moment, the sectoral competition algorithm is very simplistic. It is based on the principle that each sector is satisfied in the order of their priority. The current priority from highest to lowest is the following: domestic, livestock, thermoelectric, manufacturing, mining, irrigation (we also connected the irrigation module to our development, so irrigation also competes for water now).
Specific notes
Contributors other than yourself, if any:
I am the sole contributor in terms of code development.
But I did received support and feedback during the development and validation process. See more information about people involved in the corresponding issue section.
CTSM Issues Fixed (include github issue #): #1699
Are answers expected to change (and if so in what way)?
In the land component, there will be some small changes in the surface water-energy balance due to the fact that we apply the consumed water for sectoral usage on top of the soil columns covered by natural vegetation.
In the routing model, there will be some changes to the total water storage (VOLR) and final discharge into the oceans, since the sectoral water usage have a consumption element (meaning that the withdrawal is larger than the return flow, as some water is lost in the process and may not return to the river).
Any User Interface Changes (namelist or namelist defaults changes)?
Yes, we introduce 4 namelist variables and multiple output variables. Normally, everything is very well documented in the code, and the addition of namelist variables and output variables was done following existing CESM practices. The namelist variables added are:
sectorwater = .true. or .false. (if .true. we allow sector water abstractions; default is .false.)
limit_sectorwater_if_rof_enabled = .true. or .false. (if .true. limit abstractions based on river water availability)
sectorwater_river_volume_threshold = [0, 1] (default 0.1, which fraction of current river volume storage can be used to satisfy demand; as for irrigation, simply used to avoid negative routing or taking too much)
path_sectorwater_input_data = 'some_path' (path to input data .txt file, default ' ' )
Then we also have multiple outputs which can be specified by the user in the namelist (both for the CLM and MOSART).
Testing performed, if any:
Many tests were performed to make sure the model is performing as intended.
This include:
For all these tests we have got positive results.
Currently, we are beginning the preparation for our first paper on the topic, where we will present the development.
If during the preparation of the paper any bug will be found, we will swiftly communicate this and make the corrections. But we hope this will not happen, taking into account the extensive testing and validation, and the carefull checking of code on a dozen occasions.
For any questions, please contact me at [email protected].
Cheers,
Sabin