1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
| module SectorWaterMod
#include "shr_assert.h"
! specify all module/subroutines used
! here I will need to update the list by taking out things I didn't really used
...
!PUBLIC TYPES:
! Define sector water usage module parameters
type, public :: sectorwater_params_type
! Time of day to start domestic and livestock water use, seconds (0 = midnight).
! We start satisfying the demand in the time step FOLLOWING this time,
integer :: dom_and_liv_start_time
! Time of day to start industrial (thermoelectric, manufacturing and mining) water use, seconds (0 = midnight).
! We start satisfying the demand in the time step FOLLOWING this time,
integer :: ind_start_time
! Legth in seconds other which domestic and livestock water demand is satisfied.
! Actual time may differ if this is not a multiple of dtime.
! SectorWater module won't work properly if dtime > secsperday.
integer :: dom_and_liv_length
! Legth in seconds other which industrial water demand is satisfied.
integer :: ind_length
! Threshold for river water volume below which sectorwater usage is
! shut off, if limit_sectorwater is .true. (fraction of available river water).
! A threshold of 0 means all river water can be used; a threshold of 0.1 means 90% of the river volume can be used;
real(r8) :: sectorwater_river_volume_threshold
! Whether sectorwater usage is limited based on river storage.
! This only applies if ROF is enabled
! (i.e., rof_prognostic is .true.) - otherwise we don't limit
! sectorwater usage, regardless of the value of this flag.
logical :: limit_sectorwater_if_rof_enabled
! use groundwater supply for sectorwater usage
! (in addition to surface water)
logical :: use_groundwater_sectorwater
end type sectorwater_params_type
! Define the sectorwater_type to store dynamical quantities related to sectoral
! water usage, as well as the required procedures
type, public :: sectorwater_type
type(sectorwater_params_type) :: params
integer :: dtime ! land model time step (sec)
integer :: dom_and_liv_nsteps_per_day ! number of time steps per day in which we satisfy domestic and livestock demand
integer :: ind_nsteps_per_day ! number of time steps per day in which we satisfy industrial demand (thermoelectric, manufacturing and mining)
! Private data members; time-varying:
! naming: dom = domestic, liv = livestock, elec = thermoelectric,
! mfc = manufacturing, min = mining
! naming: withd = withdrawal, cons = consumption, rf = return flow
real(r8), pointer :: input_mon_dom_withd_grc (:) ! input expected withdrawal for current month
real(r8), pointer :: input_mon_dom_cons_grc (:) ! input expected consumption for current month
real(r8), pointer :: dom_withd_grc (:) ! expected withdrawal flux for the day [mm/s]
real(r8), pointer :: dom_cons_grc (:) ! expected consumption flux for the day [mm/s]
real(r8), pointer :: dom_withd_actual_grc (:) ! actual withdrawal flux for the day [mm/s]
real(r8), pointer :: dom_cons_actual_grc (:) ! actual consumption flux for the day [mm/s]
real(r8), pointer :: dom_rf_actual_grc (:) ! actual return flow flux for the day [mm/s]
! define same quantities for other sectors
...
integer , pointer :: n_dom_and_liv_steps_left_grc (:) ! number of time steps for which we still need to satisfy domestic and livestock demand (if 0, ignore)
integer , pointer :: n_ind_steps_left_grc (:) ! number of time steps for which we still need to satisfy industrial demand (if 0, ignore)
contains
! Public routines
procedure, public :: Init => SectorWaterInit
! procedure, public :: Restart
procedure, public :: ReadSectorWaterData
! procedure, public :: CalcSectorWaterFluxes
procedure, public :: CalcSectorWaterNeeded
procedure, public :: Clean => SectorWaterClean ! deallocate memory
! Private routines
procedure, private :: ReadNamelist
procedure, private :: CheckNamelistValidity ! Check for validity of input parameters
procedure, private :: InitAllocate => SectorWaterInitAllocate
procedure, private :: InitHistory => SectorWaterInitHistory
procedure, private :: InitCold => SectorWaterInitCold
procedure, private :: CalcSectorDemandVolrLimited ! calculate demand limited by river volume for each patch
end type sectorwater_type
interface sectorwater_params_type
module procedure sectorwater_params_constructor
end interface sectorwater_params_type
real(r8), parameter :: m3_over_km2_to_mm = 1.e-3_r8
character(len=*), parameter, private :: sourcefile = &
__FILE__
contains
! ========================================================================
! Constructors
! ========================================================================
function sectorwater_params_constructor(dom_and_liv_start_time, ind_start_time, &
dom_and_liv_length, ind_length, sectorwater_river_volume_threshold, &
limit_sectorwater_if_rof_enabled, use_groundwater_sectorwater) &
result(this)
! !DESCRIPTION:
! Create an sectorwater_params instance
! !ARGUMENTS:
type(sectorwater_params_type) :: this ! function result
integer , intent(in) :: dom_and_liv_start_time
integer , intent(in) :: ind_start_time
integer , intent(in) :: dom_and_liv_length
integer , intent(in) :: ind_length
real(r8), intent(in) :: sectorwater_river_volume_threshold
logical , intent(in) :: limit_sectorwater_if_rof_enabled
logical , intent(in) :: use_groundwater_sectorwater
! !LOCAL VARIABLES:
character(len=*), parameter :: subname = 'sectorwater_params_constructor'
!-----------------------------------------------------------------------
this%dom_and_liv_start_time = dom_and_liv_start_time
this%ind_start_time = ind_start_time
this%dom_and_liv_length = dom_and_liv_length
this%ind_length = ind_length
this%sectorwater_river_volume_threshold = sectorwater_river_volume_threshold
this%limit_sectorwater_if_rof_enabled = limit_sectorwater_if_rof_enabled
this%use_groundwater_sectorwater = use_groundwater_sectorwater
end function sectorwater_params_constructor
! ========================================================================
! Infrastructure routines (initialization, restart, etc.)
! ========================================================================
subroutine SectorWaterInit(this, bounds, NLFilename, use_aquifer_layer)
class(sectorwater_type) , intent(inout) :: this
type(bounds_type) , intent(in) :: bounds
character(len=*) , intent(in) :: NLFilename ! Namelist filename
logical , intent(in) :: use_aquifer_layer
call this%ReadNamelist(NLFilename, use_aquifer_layer)
call this%InitAllocate(bounds) ! whether an aquifer layer is used in this run
call this%InitHistory(bounds)
call this%InitCold(bounds)
end subroutine SectorWaterInit
!-----------------------------------------------------------------------
subroutine ReadNamelist(this, NLFilename, use_aquifer_layer)
! !DESCRIPTION:
! Read the sectorwater namelist
! !USES:
...
! !ARGUMENTS:
class(sectorwater_type) , intent(inout) :: this
character(len=*), intent(in) :: NLFilename ! Namelist filename
logical, intent(in) :: use_aquifer_layer ! whether an aquifer layer is used in this run
! !LOCAL VARIABLES:
! variables from sectorwater_params_type
integer :: dom_and_liv_start_time
integer :: ind_start_time
integer :: dom_and_liv_length
integer :: ind_length
real(r8) :: sectorwater_river_volume_threshold
logical :: limit_sectorwater_if_rof_enabled
logical :: use_groundwater_sectorwater
integer :: ierr ! error code
integer :: unitn ! unit for namelist file
character(len=*), parameter :: nmlname = 'sectorwater_inparm'
character(len=*), parameter :: subname = 'ReadNamelist'
!-----------------------------------------------------------------------
namelist /sectorwater_inparm/ dom_and_liv_start_time, ind_start_time, dom_and_liv_length, &
ind_length, sectorwater_river_volume_threshold, limit_sectorwater_if_rof_enabled, &
use_groundwater_sectorwater
! Initialize options to garbage defaults, forcing all to be specified
! explicitly in order to get reasonable results
dom_and_liv_start_time = 0
ind_start_time = 0
dom_and_liv_length = 0
ind_length = 0
sectorwater_river_volume_threshold = nan
limit_sectorwater_if_rof_enabled = .false.
use_groundwater_sectorwater = .false.
! Read the Namelist and fill the sectorwater_inparm namelist variables
...
! MPI_BCAST all the sectorwater_params_type variables values to
! all processes
...
! Update the sectorwater_type params with the values from namelist file
this%params = sectorwater_params_type( &
dom_and_liv_start_time = dom_and_liv_start_time, &
ind_start_time = ind_start_time, &
dom_and_liv_length = dom_and_liv_length, &
ind_length = ind_length, &
sectorwater_river_volume_threshold = sectorwater_river_volume_threshold, &
limit_sectorwater_if_rof_enabled = limit_sectorwater_if_rof_enabled, &
use_groundwater_sectorwater = use_groundwater_sectorwater)
! Write to the log file the parameters values obtained from namelist file
! Later log files can be checked to confirm that the right values were used
! Here we also call the CheckNamelistValidity(use_aquifer_layer) subroutine
...
end subroutine ReadNamelist
!-----------------------------------------------------------------------
subroutine CheckNamelistValidity(this, use_aquifer_layer)
! !DESCRIPTION:
! Check for validity of input parameters.
! Assumes that the inputs have already been packed into 'this%params'.
! Only needs to be called by the master task, since parameters are the same
! for all tasks.
! !ARGUMENTS:
class(sectorwater_type), intent(in) :: this
logical, intent(in) :: use_aquifer_layer ! whether an aquifer layer is used in this run
! !LOCAL VARIABLES:
character(len=*), parameter :: subname = 'CheckNamelistValidity'
!-----------------------------------------------------------------------
! Use associate to access the this%params% variables with shortcuts
...
! Use if statements to check if params values are in admissible values range
...
end subroutine CheckNamelistValidity
!-----------------------------------------------------------------------
subroutine SectorWaterInitAllocate(this, bounds)
! !DESCRIPTION:
! Initialize sector water data structure
! !USES:
...
! !ARGUMENTS:
class(sectorwater_type) , intent(inout) :: this
type(bounds_type) , intent(in) :: bounds
! !LOCAL VARIABLES:
integer :: begg, endg ! bounds limits
character(len=*), parameter :: subname = 'InitAllocate'
!-----------------------------------------------------------------------
begg = bounds%begg; endg= bounds%endg
! Allocate memory for variables of sectorwater_type in the limits of the
! clumps bounds (+ give nan as default values)
allocate(this%input_mon_dom_withd_grc (begg:endg)) ; this%input_mon_dom_withd_grc (:) = nan
allocate(this%input_mon_dom_cons_grc (begg:endg)) ; this%input_mon_dom_cons_grc (:) = nan
allocate(this%dom_withd_grc (begg:endg)) ; this%dom_withd_grc (:) = nan
allocate(this%dom_cons_grc (begg:endg)) ; this%dom_cons_grc (:) = nan
allocate(this%dom_withd_actual_grc (begg:endg)) ; this%dom_withd_actual_grc (:) = nan
allocate(this%dom_cons_actual_grc (begg:endg)) ; this%dom_cons_actual_grc (:) = nan
allocate(this%dom_rf_actual_grc (begg:endg)) ; this%dom_rf_actual_grc (:) = nan
! Do the same for other sectors
...
end subroutine SectorWaterInitAllocate
!-----------------------------------------------------------------------
subroutine SectorWaterInitHistory(this, bounds)
! !DESCRIPTION:
! Initialize sectoral water use history fields
! !USES:
...
! !ARGUMENTS:
class(sectorWater_type) , intent(inout) :: this
type(bounds_type) , intent(in) :: bounds
! !LOCAL VARIABLES:
integer :: begg, endg
character(len=*), parameter :: subname = 'InitHistory'
!-----------------------------------------------------------------------
begg = bounds%begg; endg= bounds%endg
! Add actual withdrawal as history fields for outputs
! Maybe later we can add other fields too
! But for the first tests, this is enough
this%dom_withd_actual_grc(begg:endg) = spval
call hist_addfld1d (fname='DOM_ACTUAL_WITHD', units='mm/s', &
avgflag='A', long_name='domestic actual withdrawal flux', &
ptr_patch=this%dom_withd_actual_grc, default='inactive')
! Do the same for other sectors
...
end subroutine SectorWaterInitHistory
!-----------------------------------------------------------------------
subroutine SectorWaterInitCold(this, bounds)
! !DESCRIPTION:
! Do cold-start initialization for sector water data structure
! !ARGUMENTS:
class(sectorwater_type) , intent(inout) :: this
type(bounds_type) , intent(in) :: bounds
character(len=*), parameter :: subname = 'InitCold'
!-----------------------------------------------------------------------
! Initialize the sectorwater_type variables
this%dtime = get_step_size()
this%ind_nsteps_per_day = 0._r8
this%dom_and_liv_nsteps_per_day = 0._r8
! actual withdrawl
this%dom_withd_actual_grc(bounds%begg:bounds%endg) = 0._r8
... ! same for other sectors
! actual consumption
this%dom_cons_actual_grc(bounds%begg:bounds%endg) = 0._r8
... ! same for other sectors
! actual return flow
this%dom_rf_actual_grc(bounds%begg:bounds%endg) = 0._r8
... ! same for other sectors
end subroutine SectorWaterInitCold
!-----------------------------------------------------------------------
pure function Calc_dom_and_liv_NstepsPerDay(this, dtime) result(dom_and_liv_nsteps_per_day)
! !DESCRIPTION:
! Given dtime (sec), determine number of steps per day to satisfy demand for domestic and livestock sectors
...
end function Calc_dom_and_liv_NstepsPerDay
!-----------------------------------------------------------------------
pure function Calc_ind_NstepsPerDay(this, dtime) result(ind_nsteps_per_day)
! !DESCRIPTION:
! Given dtime (sec), determine number of steps per day to satisfy demand for industrial sectors
...
end function Calc_ind_NstepsPerDay
!-----------------------------------------------------------------------
subroutine SectorWaterClean(this)
! !DESCRIPTION:
! Deallocate memory for all sectorwater_type variables
...
end subroutine sectorWaterClean
! ========================================================================
! Science routines
! ========================================================================
subroutine ReadSectorWaterData (this, bounds, mon)
! !DESCRIPTION:
! read the input data from fsurfdata (withdrawal and consumption)
! at the moment fsurfdata for sector water module test have the data
! for year 2000 in monthly format for all sectors
! when we will finish development we will not rely on surfdata, but on
! alternative methods more convenient for transient data
! !USES:
...
! !ARGUMENTS:
class(sectorwater_type), intent(inout) :: this
type(bounds_type) , intent(in) :: bounds
integer , intent(in) :: mon ! month (1, ..., 12) for nstep+1
! !LOCAL VARIABLES:
type(file_desc_t) :: ncid ! netcdf id
integer :: ier ! error code
integer :: g ! gridcell index
integer :: ni,nj,ns ! indices
integer :: dimid,varid ! input netCDF id's
integer :: ntim ! number of input data time samples
integer :: nlon_i ! number of input data longitudes
integer :: nlat_i ! number of input data latitudes
logical :: isgrid2d ! true => file is 2d
real(r8), pointer :: mon_dom_withd(:) ! monthly domestic withdrawal read from input files
real(r8), pointer :: mon_dom_cons(:) ! monthly domestic consumption read from input files
... ! Same for other sectors
character(len=256) :: locfn ! local file name
character(len=32) :: subname = 'ReadSectorWaterData'
!-----------------------------------------------------------------------
if (masterproc) then
write (iulog,*) 'Attempting to read annual sectoral water usage data .....'
end if
allocate(&
mon_dom_withd(bounds%begg:bounds%endg), &
mon_dom_cons(bounds%begg:bounds%endg), &
..., & ! same for other sectors
stat=ier)
! Determine necessary indices
call getfil(fsurdat, locfn, 0)
call ncd_pio_openfile (ncid, trim(locfn), 0)
call ncd_inqfdims (ncid, isgrid2d, ni, nj, ns)
! Check if ldomain and input file have matching dimensions (ni, nj and ns)
! if not then abort run
...
! read the data for current month
call ncd_io(ncid=ncid, varname='withd_dom', flag='read', data=mon_dom_withd, &
dim1name=grlnd, nt=mon)
call ncd_io(ncid=ncid, varname='cons_dom', flag='read', data=mon_dom_cons, &
dim1name=grlnd, nt=mon)
... ! do the same for other sectors
call ncd_pio_closefile(ncid) ! close input file
! fill the sectorwater_type input fields with current month
! withdrawal and consumption (loop over each gridcell)
do g = bounds%begg,bounds%endg
this%input_mon_dom_withd_grc(g) = mon_dom_withd(g)
this%input_mon_dom_cons_grc(g) = mon_dom_cons(g)
... ! do the same for other sectors
end do
end subroutine ReadSectorWaterData
!-----------------------------------------------------------------------
subroutine CalcSectorWaterNeeded(this, bounds, volr, rof_prognostic)
! !DESCRIPTION:
! Calculate sector water fluxes (withdrawal, consumption)
! !USES:
use shr_const_mod , only : SHR_CONST_TKFRZ
use clm_time_manager , only : get_curr_date, is_end_curr_month, get_curr_days_per_year
! !ARGUMENTS:
class(sectorwater_type) , intent(inout) :: this
type(bounds_type) , intent(in) :: bounds
! river water volume (m3) (ignored if rof_prognostic is .false.)
real(r8), intent(in) :: volr(bounds%begg:bounds%endg)
! whether we're running with a prognostic ROF component;
! this is needed to determine whether we can limit demand based on river volume.
logical, intent(in) :: rof_prognostic
! !LOCAL VARIABLES:
integer :: g ! gridcell index
integer :: year ! year (0, ...) for nstep+1
integer :: mon ! month (1, ..., 12) for nstep+1
integer :: day ! day of month (1, ..., 31) for nstep+1
integer :: sec ! seconds into current date for nstep+1
real(r8) :: dayspyr ! days per year
real(r8) :: dayspm ! days per month
real(r8) :: secs_per_day ! seconds per day
real(r8) :: dom_and_liv_flux_factor ! factor to transform the demand from mm per month to mm/s for the given day
real(r8) :: ind_flux_factor ! factor to transform the demand from mm per month to mm/s for the given day
real(r8) :: dom_demand(bounds%begg:bounds%endg)
real(r8) :: dom_demand_volr_limited(bounds%begg:bounds%endg)
real(r8) :: dom_consumption(bounds%begg:bounds%endg)
real(r8) :: dom_consumption_volr_limited(bounds%begg:bounds%endg)
... ! Same for other sectors
! Whether we should limit deficits by available volr
logical :: limit_sectorwater
character(len=*), parameter :: subname = 'CalcSectorWaterNeeded'
!-----------------------------------------------------------------------
! Get current date
call get_curr_date(year, mon, day, sec)
! Get number of days in current year
dayspyr = get_curr_days_per_year()
! Compute average number of days per month for the current year
dayspm = dayspyr/12_r8
! Compute the flux factors to transform from mm/month to mm/s
dom_and_liv_flux_factor = ((1_r8/dayspm)/this%params%dom_and_liv_length)
ind_flux_factor = ((1_r8/dayspm)/this%params%ind_length)
! Read input for new month if end of month
! This will update the input fields of withdrawal and consumption
! of sectorwater_type
if (is_end_curr_month()) then
call this%ReadSectorWaterData(bounds, mon)
endif
! Compute demand [mm]
! First initialize demand to 0 everywhere;
dom_demand(bounds%begg:bounds%endg) = 0._r8
dom_consumption(bounds%begg:bounds%endg) = 0._r8
... ! same for other sectors
! Update this day fluxes based on input data
do g = bounds%begg,bounds%endg
dom_demand(g) = this%input_mon_dom_withd_grc(g)
dom_consumption(g) = this%input_mon_dom_cons_grc(g)
... ! same for other sectors
end do ! end loop over gridcels
! Limit sectoral withdrawal based on volr (river volume)
! Note that we cannot do this limiting if running without a prognostic
! river model, since we need river volume for this limiting.
limit_sectorwater = (this%params%limit_sectorwater_if_rof_enabled .and. rof_prognostic)
if (limit_sectorwater) then
! Compute the allowed withdrawal if water limitation is taken into account
call this%CalcSectorDemandVolrLimited( &
bounds = bounds, &
dom_demand = dom_demand(bounds%begg:bounds%endg), &
dom_consumption = dom_consumption(bounds%begg:bounds%endg), &
liv_demand = liv_demand(bounds%begg:bounds%endg), &
liv_consumption = liv_consumption(bounds%begg:bounds%endg), &
elec_demand = elec_demand(bounds%begg:bounds%endg), &
elec_consumption = elec_consumption(bounds%begg:bounds%endg), &
mfc_demand = mfc_demand(bounds%begg:bounds%endg), &
mfc_consumption = mfc_consumption(bounds%begg:bounds%endg), &
min_demand = min_demand(bounds%begg:bounds%endg), &
min_consumption = min_consumption(bounds%begg:bounds%endg), &
volr = volr(bounds%begg:bounds%endg), &
dom_demand_volr_limited = dom_demand_volr_limited(bounds%begg:bounds%endg), &
dom_consumption_volr_limited = dom_consumption_volr_limited(bounds%begg:bounds%endg), &
liv_demand_volr_limited = liv_demand_volr_limited(bounds%begg:bounds%endg), &
liv_consumption_volr_limited = liv_consumption_volr_limited(bounds%begg:bounds%endg), &
elec_demand_volr_limited = elec_demand_volr_limited(bounds%begg:bounds%endg), &
elec_consumption_volr_limited = elec_consumption_volr_limited(bounds%begg:bounds%endg), &
mfc_demand_volr_limited = mfc_demand_volr_limited(bounds%begg:bounds%endg), &
mfc_consumption_volr_limited = mfc_consumption_volr_limited(bounds%begg:bounds%endg), &
min_demand_volr_limited = min_demand_volr_limited(bounds%begg:bounds%endg), &
min_consumption_volr_limited = min_consumption_volr_limited(bounds%begg:bounds%endg))
else
! If water limitation is not considered
! the consider the volr_limited withdrawals to be the same
! as input based withdrawl
dom_demand_volr_limited(bounds%begg:bounds%endg) = dom_demand(bounds%begg:bounds%endg)
dom_consumption_volr_limited(bounds%begg:bounds%endg) = dom_consumption(bounds%begg:bounds%endg)
... ! do the same for the other sectors
end if
! Convert demand to withdrawal rates [mm/s]
do g = bounds%begg,bounds%endg
! compute expected and actual withdrawal and consumption
this%dom_withd_grc(g) = dom_demand(g)*dom_and_liv_flux_factor
this%dom_withd_actual_grc(g) = dom_demand_volr_limited(g)*dom_and_liv_flux_factor
this%dom_cons_grc(g) = dom_consumption(g)*dom_and_liv_flux_factor
this%dom_cons_actual_grc(g) = dom_consumption_volr_limited(g)*dom_and_liv_flux_factor
! computed actual return flow
this%dom_rf_actual_grc(g) = this%dom_withd_actual_grc(g) - this%dom_cons_actual_grc(g)
! Do the same for other sectors
...
end do
end subroutine CalcSectorWaterNeeded
!-----------------------------------------------------------------------
subroutine CalcSectorDemandVolrLimited(this, bounds, dom_demand, dom_consumption, liv_demand, liv_consumption, elec_demand, elec_consumption, &
mfc_demand, mfc_consumption, min_demand, min_consumption, volr, dom_demand_volr_limited, dom_consumption_volr_limited, liv_demand_volr_limited, &
liv_consumption_volr_limited, elec_demand_volr_limited, elec_consumption_volr_limited, mfc_demand_volr_limited, &
mfc_consumption_volr_limited, min_demand_volr_limited, min_consumption_volr_limited)
! !ARGUMENTS:
class(sectorwater_type) , intent(in) :: this
type(bounds_type) , intent(in) :: bounds
real(r8), intent(in) :: dom_demand( bounds%begg:bounds%endg)
real(r8), intent(in) :: dom_consumption( bounds%begg:bounds%endg)
... ! same for other sectors
! river water volume [m3]
real(r8), intent(in) :: volr( bounds%begg:bounds%endg)
real(r8), intent(out) :: dom_demand_volr_limited( bounds%begg:bounds%endg)
real(r8), intent(out) :: dom_consumption_volr_limited( bounds%begg:bounds%endg)
... ! same for other sectors
! !LOCAL VARIABLES:
integer :: g ! gridcell index
real(r8) :: available_volr ! volr available for withdrawal [m3]
real(r8) :: max_demand_supported_by_volr ! [kg/m2] [i.e., mm]
! ratio of demand_volr_limited to demand for each grid cell
real(r8) :: dom_demand_limited_ratio_grc(bounds%begg:bounds%endg)
... ! same for other sectors
character(len=*), parameter :: subname = 'CalcSectorDemandVolrLimited'
!-----------------------------------------------------------------------
do g = bounds%begg, bounds%endg
if (volr(g) > 0._r8) then
available_volr = volr(g) * (1._r8 - this%params%sectorwater_river_volume_threshold)
max_demand_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_demand_supported_by_volr = 0._r8
end if
if (dom_demand(g) > max_demand_supported_by_volr) then
dom_demand_limited_ratio_grc(g) = max_demand_supported_by_volr / dom_demand(g)
liv_demand_limited_ratio_grc(g) = 0._r8
elec_demand_limited_ratio_grc(g) = 0._r8
mfc_demand_limited_ratio_grc(g) = 0._r8
min_demand_limited_ratio_grc(g) = 0._r8
else if (liv_demand(g) > (max_demand_supported_by_volr - dom_demand(g))) then
dom_demand_limited_ratio_grc(g) = 1._r8
liv_demand_limited_ratio_grc(g) = max_demand_supported_by_volr / liv_demand(g)
elec_demand_limited_ratio_grc(g) = 0._r8
mfc_demand_limited_ratio_grc(g) = 0._r8
min_demand_limited_ratio_grc(g) = 0._r8
else if (elec_demand(g) > (max_demand_supported_by_volr - dom_demand(g) - liv_demand(g))) then
dom_demand_limited_ratio_grc(g) = 1._r8
liv_demand_limited_ratio_grc(g) = 1._r8
elec_demand_limited_ratio_grc(g) = max_demand_supported_by_volr/ elec_demand(g)
mfc_demand_limited_ratio_grc(g) = 0._r8
min_demand_limited_ratio_grc(g) = 0._r8
else if (mfc_demand(g) > (max_demand_supported_by_volr - dom_demand(g) - liv_demand(g) - elec_demand(g))) then
dom_demand_limited_ratio_grc(g) = 1._r8
liv_demand_limited_ratio_grc(g) = 1._r8
elec_demand_limited_ratio_grc(g) = 1._r8
mfc_demand_limited_ratio_grc(g) = max_demand_supported_by_volr / mfc_demand(g)
min_demand_limited_ratio_grc(g) = 0._r8
else if (min_demand(g) > (max_demand_supported_by_volr - dom_demand(g) - liv_demand(g) - elec_demand(g) - mfc_demand(g))) then
dom_demand_limited_ratio_grc(g) = 1._r8
liv_demand_limited_ratio_grc(g) = 1._r8
elec_demand_limited_ratio_grc(g) = 1._r8
mfc_demand_limited_ratio_grc(g) = 1._r8
min_demand_limited_ratio_grc(g) = max_demand_supported_by_volr / min_demand(g)
else
dom_demand_limited_ratio_grc(g) = 1._r8
liv_demand_limited_ratio_grc(g) = 1._r8
elec_demand_limited_ratio_grc(g) = 1._r8
mfc_demand_limited_ratio_grc(g) = 1._r8
min_demand_limited_ratio_grc(g) = 1._r8
end if
end do
dom_demand_volr_limited(bounds%begg:bounds%endg) = 0._r8
dom_consumption_volr_limited(bounds%begg:bounds%endg) = 0._r8
... ! same for other sectors
do g = bounds%begg, bounds%endg
dom_demand_volr_limited(g) = dom_demand(g) * dom_demand_limited_ratio_grc(g)
dom_consumption_volr_limited(g) = dom_consumption(g) * dom_demand_limited_ratio_grc(g)
... ! same for other sectors
end do
end subroutine CalcSectorDemandVolrLimited
end module SectorWaterMod
|