Context Navigation

source: branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/routing_native_irrig.f90 @ 8366

Last change on this file since 8366 was 8366, checked in by yann.meurdesoif, 6 months ago

Update irrigation for native routing

fix problem due to an uninitialized array that must be set to 0 for Pedro irrigation

YM

File size: 34.2 KB

Line
1	MODULE routing_native_irrig_mod
2	USE constantes
3
4
5	PRIVATE
6
7	INTEGER, SAVE :: nbpt
8	!$OMP THREADPRIVATE(nbpt)
9
10	INTEGER, SAVE :: nbpt_r
11	!$OMP THREADPRIVATE(nbpt_r)
12
13	REAL(r_std), ALLOCATABLE, SAVE :: irrig_gw_source_r(:) ! diag
14	REAL(r_std), ALLOCATABLE, SAVE :: irrig_fast_source_r(:) ! diag
15	REAL(r_std), ALLOCATABLE, SAVE :: irrig_str_source_r(:) ! diag
16
17	REAL(r_std), ALLOCATABLE, SAVE :: vegtot_mean(:)
18	!$OMP THREADPRIVATE(vegtot_mean)
19
20	REAL(r_std), ALLOCATABLE, SAVE :: humrel_mean(:)
21	!$OMP THREADPRIVATE(humrel_mean)
22
23	REAL(r_std), ALLOCATABLE, SAVE :: transpot_mean(:)
24	!$OMP THREADPRIVATE(transpot_mean)
25
26	REAL(r_std), ALLOCATABLE, SAVE :: runoff_mean(:)
27	!$OMP THREADPRIVATE(runoff_mean)
28
29	REAL(r_std), ALLOCATABLE, SAVE :: precip_mean(:)
30	!$OMP THREADPRIVATE(precip_mean)
31
32	REAL(r_std), ALLOCATABLE, SAVE :: irrigation_mean(:)
33	!$OMP THREADPRIVATE(irrigation_mean)
34
35	REAL(r_std), ALLOCATABLE, SAVE :: irrigated(:)
36	!$OMP THREADPRIVATE(irrigated)
37
38	PUBLIC irrigation_initialize, irrigation_main, irrigation_mean_make, irrigation_get, irrigation_finalize
39
40	CONTAINS
41
42	SUBROUTINE irrigation_get(irrigation_mean)
43	IMPLICIT NONE
44	REAL(r_std),OPTIONAL, INTENT(OUT) :: irrigation_mean(:)
45
46	CALL irrigation_get_(irrigation_mean)
47	END SUBROUTINE irrigation_get
48
49	SUBROUTINE irrigation_get_(irrigation_mean_)
50	IMPLICIT NONE
51	REAL(r_std),OPTIONAL, INTENT(OUT) :: irrigation_mean_(:)
52	IF (PRESENT(irrigation_mean_)) irrigation_mean_ = irrigation_mean
53
54	END SUBROUTINE irrigation_get_
55
56
57	SUBROUTINE irrigation_initialize(kjit, rest_id, nbpt_, nbpt_r_, irrigated_next )
58	USE constantes
59	IMPLICIT NONE
60	INTEGER, INTENT(IN) :: kjit
61	INTEGER, INTENT(IN) :: rest_id
62	INTEGER, INTENT(IN) :: nbpt_
63	INTEGER, INTENT(IN) :: nbpt_r_
64	REAL(r_std), INTENT(IN) :: irrigated_next(nbpt_)
65
66	CALL irrigation_local_init(kjit, rest_id, nbpt_, nbpt_r_,irrigated_next)
67	CALL irrigation_mean_init(kjit, rest_id)
68
69	END SUBROUTINE irrigation_initialize
70
71
72	SUBROUTINE irrigation_finalize(kjit, rest_id)
73	IMPLICIT NONE
74	INTEGER, INTENT(IN) :: kjit
75	INTEGER, INTENT(IN) :: rest_id
76
77	CALL irrigation_mean_finalize(kjit, rest_id)
78
79	END SUBROUTINE irrigation_finalize
80
81
82	SUBROUTINE irrigation_local_init(kjit, rest_id, nbpt_, nbpt_r_, irrigated_next)
83	IMPLICIT NONE
84	INTEGER, INTENT(IN) :: kjit
85	INTEGER, INTENT(IN) :: rest_id
86	INTEGER, INTENT(IN) :: nbpt_
87	INTEGER, INTENT(IN) :: nbpt_r_
88	REAL(r_std), INTENT(IN) :: irrigated_next(nbpt_)
89
90	nbpt = nbpt_
91	nbpt_r = nbpt_r_
92	IF (do_irrigation) THEN
93	ALLOCATE(irrigated(nbpt))
94	irrigated(:)=irrigated_next(:)
95	ALLOCATE(irrig_gw_source_r(nbpt_r))
96	ALLOCATE(irrig_fast_source_r(nbpt_r))
97	ALLOCATE(irrig_str_source_r(nbpt_r))
98	ENDIF
99
100	END SUBROUTINE irrigation_local_init
101
102
103
104	SUBROUTINE irrigation_mean_init(kjit, rest_id)
105	USE ioipsl_para
106	USE grid
107	USE sechiba_io_p
108	IMPLICIT NONE
109	INTEGER, INTENT(IN) :: kjit
110	INTEGER, INTENT(IN) :: rest_id
111	CHARACTER(LEN=80) :: var_name !! To store variables names for I/O (unitless)
112	INTEGER(i_std) :: ier
113
114	ALLOCATE(irrigation_mean(nbpt), stat=ier)
115	irrigation_mean(:) = 0
116	IF (do_irrigation) THEN
117	IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for irrigation_mean','','')
118	var_name = 'irrigation'
119	CALL ioconf_setatt_p('UNITS', 'Kg/dt')
120	CALL ioconf_setatt_p('LONG_NAME','Artificial irrigation flux')
121	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., irrigation_mean, "gather", nbp_glo, index_g)
122	CALL setvar_p (irrigation_mean, val_exp, 'NO_KEYWORD', zero)
123
124	ALLOCATE(vegtot_mean(nbpt), stat=ier)
125	ALLOCATE(humrel_mean(nbpt), stat=ier)
126	ALLOCATE(transpot_mean(nbpt), stat=ier)
127	ALLOCATE (runoff_mean(nbpt), stat=ier)
128	ALLOCATE(precip_mean(nbpt), stat=ier)
129	CALL irrigation_mean_reset
130	ENDIF
131
132	END SUBROUTINE irrigation_mean_init
133
134
135	SUBROUTINE irrigation_mean_reset
136	IMPLICIT NONE
137
138	vegtot_mean(:) = 0
139	humrel_mean(:) = 0
140	transpot_mean(:) = 0
141	runoff_mean(:) = 0
142	precip_mean(:) = 0
143
144	END SUBROUTINE irrigation_mean_reset
145
146	SUBROUTINE irrigation_mean_make(dt_routing, veget_max, humrel, transpot, runoff ,precip )
147	USE pft_parameters
148	IMPLICIT NONE
149	REAL(r_std),INTENT(IN) :: dt_routing
150	REAL(r_std),INTENT(IN) :: veget_max(:,:)
151	REAL(r_std),INTENT(IN) :: humrel(:,:)
152	REAL(r_std),INTENT(IN) :: transpot(:,:)
153	REAL(r_std),INTENT(IN) :: runoff(:)
154	REAL(r_std),INTENT(IN) :: precip(:)
155
156	REAL(r_std), DIMENSION(nbpt) :: tot_vegfrac_nowoody !! Total fraction occupied by grass (0-1,unitless)
157	INTEGER :: jv, ig
158
159	IF (do_irrigation) THEN
160	runoff_mean(:) = runoff_mean(:) + runoff
161	precip_mean(:) = precip_mean(:) + precip
162
163	!! Computes the total fraction occupied by the grasses and the crops for each grid cell
164	tot_vegfrac_nowoody(:) = zero
165	DO jv = 1, nvm
166	IF ( (jv /= ibare_sechiba) .AND. .NOT.(is_tree(jv)) ) THEN
167	tot_vegfrac_nowoody(:) = tot_vegfrac_nowoody(:) + veget_max(:,jv)
168	END IF
169	END DO
170
171	DO ig = 1, nbpt
172	IF ( tot_vegfrac_nowoody(ig) .GT. min_sechiba ) THEN
173	DO jv = 1,nvm
174	IF ( (jv /= ibare_sechiba) .AND. .NOT.(is_tree(jv)) ) THEN
175	transpot_mean(ig) = transpot_mean(ig) + transpot(ig,jv) * veget_max(ig,jv)/tot_vegfrac_nowoody(ig)
176	END IF
177	END DO
178	ELSE
179	IF (MAXVAL(veget_max(ig,2:nvm)) .GT. min_sechiba) THEN
180	DO jv = 2, nvm
181	transpot_mean(ig) = transpot_mean(ig) + transpot(ig,jv) * veget_max(ig,jv)/ SUM(veget_max(ig,2:nvm))
182	ENDDO
183	ENDIF
184	ENDIF
185	ENDDO
186
187	! New irrigation scheme uses mean of vegtot with jv 1 to nvm
188	! Old scheme keeps jv 2 to nvm, even if possibly an error
189	IF ( .NOT. old_irrig_scheme ) THEN
190	DO jv=1,nvm
191	DO ig=1,nbpt
192	humrel_mean(ig) = humrel_mean(ig) + humrel(ig,jv)veget_max(ig,jv)dt_sechiba/dt_routing
193	vegtot_mean(ig) = vegtot_mean(ig) + veget_max(ig,jv)*dt_sechiba/dt_routing
194	ENDDO
195	ENDDO
196	ELSE
197	DO jv=2,nvm
198	DO ig=1,nbpt
199	humrel_mean(ig) = humrel_mean(ig) + humrel(ig,jv)veget_max(ig,jv)dt_sechiba/dt_routing
200	vegtot_mean(ig) = vegtot_mean(ig) + veget_max(ig,jv)*dt_sechiba/dt_routing
201	ENDDO
202	ENDDO
203	ENDIF
204	ENDIF !do_irrigation
205	END SUBROUTINE irrigation_mean_make
206
207	SUBROUTINE irrigation_mean_finalize(kjit, rest_id)
208	USE grid
209	USE ioipsl_para
210	IMPLICIT NONE
211	INTEGER, INTENT(IN) :: kjit
212	INTEGER, INTENT(IN) :: rest_id
213
214	IF (do_irrigation) THEN
215	CALL restput_p (rest_id, 'irrigation', nbp_glo, 1, 1, kjit, irrigation_mean, 'scatter', nbp_glo, index_g)
216	ENDIF
217
218	END SUBROUTINE irrigation_mean_finalize
219
220
221	SUBROUTINE irrigation_main(dt_routing, contfrac, reinfiltration, irrigated_next, irrig_frac, root_deficit, soiltile, &
222	fraction_aeirrig_sw )
223	USE constantes
224	USE constantes_soil
225	USE grid, ONLY : area
226	USE xios
227	IMPLICIT NONE
228	REAL(r_std), INTENT(IN) :: dt_routing
229	REAL(r_std), INTENT(IN) :: contfrac(nbpt)
230	REAL(r_std), INTENT(IN) :: reinfiltration(nbpt)
231	REAL(r_std), INTENT(IN) :: irrigated_next(nbpt)
232	REAL(r_std), INTENT(IN) :: irrig_frac(nbpt) !! Irrig. fraction interpolated in routing, and saved to pass to slowproc if irrigated_soiltile = .TRUE.
233	REAL(r_std), INTENT(IN) :: root_deficit(nbpt) !! soil water deficit
234	REAL(r_std), INTENT(IN) :: soiltile(nbpt,nstm) !! Fraction of each soil tile within vegtot (0-1, unitless)
235	REAL(r_std), INTENT(IN) :: fraction_aeirrig_sw(nbpt) !! Fraction of area equipped for irrigation from surface water, of irrig_frac
236
237	REAL(r_std) :: irrig_needs_r(nbpt_r) !! Total irrigation requirement (water requirements by the crop for its optimal growth) (kg)
238	REAL(r_std) :: irrig_actual_r(nbpt_r) !! Possible irrigation according to the water availability in the reservoirs (kg)
239	REAL(r_std) :: irrig_actual(nbpt) !! Possible irrigation according to the water availability in the reservoirs (kg)
240
241	IF (do_irrigation) THEN
242	IF (irrig_map_dynamic_flag ) irrigated=irrigated_next
243
244	IF (old_irrig_scheme) THEN
245	CALL irrigation_compute_requested_old(contfrac, reinfiltration, irrigated, irrig_needs_r)
246	CALL irrigation_old(irrig_needs_r, irrig_actual_r)
247	ELSE
248	CALL irrigation_compute_requested_new(dt_routing, contfrac, root_deficit, soiltile, irrig_frac, irrig_needs_r)
249	CALL irrigation_new(contfrac,fraction_aeirrig_sw, irrig_needs_r, irrig_actual_r)
250	ENDIF
251
252	IF (is_omp_root) THEN
253	CALL xios_send_field("routing_irrigation_r",irrig_actual_r)
254	CALL xios_recv_field("routing_irrigation",irrig_actual)
255	ENDIF
256	CALL scatter_omp(irrig_actual,irrigation_mean)
257	irrigation_mean(:)=irrigation_mean(:)/(area(:)*contfrac(:))
258	CALL irrigation_mean_reset()
259	ELSE
260	irrigation_mean(:) = 0
261	ENDIF
262
263	END SUBROUTINE irrigation_main
264
265
266
267
268
269	SUBROUTINE irrigation_new(contfrac, fraction_aeirrig_sw, irrig_needs_r, irrig_actual_r)
270	USE constantes
271	USE routing_native_flow_mod, ONLY : routing_mask_r, routing_flow_get, routing_flow_set
272	USE mod_orchidee_para
273	USE xios
274	IMPLICIT NONE
275	REAL(r_std),INTENT(IN) :: contfrac(nbpt)
276	REAL(r_std),INTENT(IN) :: fraction_aeirrig_sw(nbpt)
277	REAL(r_std),INTENT(IN) :: irrig_needs_r(nbpt_r) !! Total irrigation requirement (water requirements by the crop for its optimal growth) (kg)
278	REAL(r_std),INTENT(OUT) :: irrig_actual_r(nbpt_r) !! Possible irrigation according to the water availability in the reservoirs (kg)
279
280	REAL(r_std) :: irrig_deficit_r(nbpt_r) !! Amount of water missing for irrigation (kg)
281
282	LOGICAL :: IsFail_slow !! Logical to ask if slow reserv. does not fit irrigation demand
283	LOGICAL :: IsFail_fast !! Logical to ask if fast reserv. does not fit irrigation demand
284	LOGICAL :: IsFail_stre !! Logical to ask if stream reserv. does not fit irrigation demand
285
286	REAL(r_std) :: Count_failure_slow(nbpt_r) !! Counter times slow reserv. does not fit irrigation demand
287	REAL(r_std) :: Count_failure_fast(nbpt_r) !! Counter times fast reserv. does not fit irrigation demand
288	REAL(r_std) :: Count_failure_stre(nbpt_r) !! Counter times stream reserv. does not fit irrigation demand
289
290	REAL(r_std) :: pot_slow_wdr_dummy, pot_fast_wdr_dummy, pot_stre_wdr_dummy
291	REAL(r_std) :: slow_wdr_dummy, fast_wdr_dummy, stre_wdr_dummy
292
293	REAL(r_std) :: slow_reservoir_r(nbpt_r)
294	REAL(r_std) :: fast_reservoir_r(nbpt_r)
295	REAL(r_std) :: stream_reservoir_r(nbpt_r)
296	REAL(r_std) :: fraction_aeirrig_sw_mpi(nbp_mpi)
297	REAL(r_std) :: fraction_aeirrig_sw_r(nbpt_r)
298	REAL(r_std) :: pcent_vol_irrig
299	INTEGER :: ig
300
301
302	CALL gather_omp(fraction_aeirrig_sw, fraction_aeirrig_sw_mpi)
303	!! fraction moyenne => need contfrac ??
304
305	IF (is_omp_root) THEN
306
307	CALL xios_send_field("fraction_aeirrig_sw", fraction_aeirrig_sw_mpi)
308	CALL xios_recv_field("fraction_aeirrig_sw_r", fraction_aeirrig_sw_r) ! ==> need conservative quantity interp => no flux => to check !!!
309
310
311	Count_failure_slow(:) = zero !
312	Count_failure_fast(:) = zero !
313	Count_failure_stre(:) = zero !
314
315	CALL routing_flow_get(slow_reservoir_r=slow_reservoir_r, fast_reservoir_r=fast_reservoir_r, stream_reservoir_r=stream_reservoir_r)
316
317	DO ig=1,nbpt_r
318
319	IF (routing_mask_r(ig)) THEN
320
321	IF (select_source_irrig) THEN
322
323	! For irrig. scheme, available_reserve gives the amount of water available for irrigation in every reservoir
324	! --> avail_reserve is a vector of dimension=3, BY DEFINITION i=1 for streamflow, i=2 fast, and i=3 slow reservoir
325
326	! The new priorization scheme takes into account irrig. infrastructur according to GMIA map
327	! It also withdraw water according to availability, it means that it wont seek for all the water in then
328	! stream reservoir, even if this one could respond to the demand by itself
329
330	pot_slow_wdr_dummy = ( 1 - fraction_aeirrig_sw_r(ig)) * avail_reserve(3)*slow_reservoir_r(ig)
331	pot_fast_wdr_dummy = fraction_aeirrig_sw_r(ig) * avail_reserve(2) * fast_reservoir_r(ig)
332	pot_stre_wdr_dummy = fraction_aeirrig_sw_r(ig) * avail_reserve(1) * stream_reservoir_r(ig)
333	pcent_vol_irrig = zero
334	IsFail_slow = .FALSE. !
335	IsFail_fast = .FALSE. !
336	IsFail_stre = .FALSE. !
337
338	irrig_actual_r(ig) = MIN(irrig_needs_r(ig), pot_stre_wdr_dummy + pot_fast_wdr_dummy + pot_slow_wdr_dummy)
339
340	!! additional IF to calculate pcent_vol_irrig, in the case the total avail.
341	!! water is zero, I.E. when there is no water in surface and fraction_ae = 0,
342	!! so GW is not taken into account
343	!! Note on pcent_vol_irrig: It correspond to the fraction of available water in surface,
344	!! considering environmental needs and irrigation equipement by source from map. It controls
345	!! how the source of water withdrawl, especially when requirements < available water
346
347	IF ( (pot_stre_wdr_dummy + pot_fast_wdr_dummy + pot_slow_wdr_dummy) .GT. min_sechiba ) THEN
348
349	pcent_vol_irrig = ( pot_stre_wdr_dummy + pot_fast_wdr_dummy ) / &
350	( pot_stre_wdr_dummy + pot_fast_wdr_dummy + pot_slow_wdr_dummy )
351
352	!Irrig_actual set to zero, because there is no available water.
353	!Put to avoid negative values due to problems in the Min function
354
355	irrig_actual_r(ig) = MAX(irrig_actual_r(ig), zero)
356
357	!Already zero because pcent_vol_irrig initialized to zero
358	!Put here to readability but not necessary
359	!ELSE
360	! pcent_vol_irrig = zero
361
362	ENDIF
363
364	!Note for irrig_gw_source(ig): first we add the slow_reservoir volume. Then we substract the updated slow_reservoir. It should be the
365	!Volume used for irrigation that comes from GW
366	! Idem for irrig_fast_source and irrig_str_source
367
368	slow_wdr_dummy = slow_reservoir_r(ig)
369	slow_reservoir_r(ig) = MAX( (un - ( un - fraction_aeirrig_sw_r(ig) ) * avail_reserve(3) ) * &
370	slow_reservoir_r(ig), slow_reservoir_r(ig) + &
371	MIN( - irrig_actual_r(ig) * (un - pcent_vol_irrig ), &
372	avail_reserve(2) * fraction_aeirrig_sw_r(ig) * fast_reservoir_r(ig) + &
373	MIN(zero, avail_reserve(1) * fraction_aeirrig_sw_r(ig) * stream_reservoir_r(ig) - &
374	pcent_vol_irrig * irrig_actual_r(ig) ) ) )
375
376	slow_wdr_dummy = slow_wdr_dummy - slow_reservoir_r(ig)
377	irrig_gw_source_r(ig) = irrig_gw_source_r(ig) + slow_wdr_dummy
378
379	fast_wdr_dummy = fast_reservoir_r(ig)
380	fast_reservoir_r(ig) = MAX( (un - avail_reserve(2) * fraction_aeirrig_sw_r(ig) ) * fast_reservoir_r(ig) , &
381	fast_reservoir_r(ig) + MIN(zero, avail_reserve(1) * fraction_aeirrig_sw_r(ig) * stream_reservoir_r(ig) - &
382	pcent_vol_irrig * irrig_actual_r(ig) ) )
383	fast_wdr_dummy = fast_wdr_dummy - fast_reservoir_r(ig)
384	irrig_fast_source_r(ig) = irrig_fast_source_r(ig) + fast_wdr_dummy
385
386	stre_wdr_dummy = stream_reservoir_r(ig)
387	stream_reservoir_r(ig) = MAX((un - avail_reserve(1)* fraction_aeirrig_sw_r(ig) )*stream_reservoir_r(ig), &
388	stream_reservoir_r(ig) - pcent_vol_irrig * irrig_actual_r(ig) )
389	stre_wdr_dummy = stre_wdr_dummy - stream_reservoir_r(ig)
390	irrig_str_source_r(ig) = irrig_str_source_r(ig) + stre_wdr_dummy
391
392	irrig_deficit_r(ig) = irrig_needs_r(ig)-irrig_actual_r(ig)
393
394	!A reservoir is failing to give water for infiltration if pot. req > pot. withdrawal
395	!We assume that the pot. requirement = Needs * fraction of area equipped for SW/GW
396	!In the case of surface. we also sustract the withdrawal from Fast/Stream, because both are
397	! considered as surface water
398
399	IsFail_slow = ( ( irrig_needs_r(ig)*(un - fraction_aeirrig_sw_r(ig)) ) > pot_slow_wdr_dummy )
400	IsFail_fast = ( ( irrig_needs_r(ig)*fraction_aeirrig_sw_r(ig) - stre_wdr_dummy ) > pot_fast_wdr_dummy )
401	IsFail_stre = ( ( irrig_needs_r(ig)*fraction_aeirrig_sw_r(ig) - fast_wdr_dummy ) > pot_stre_wdr_dummy )
402
403	IF( IsFail_stre ) Count_failure_stre(ig) = un
404	IF( IsFail_fast ) Count_failure_fast(ig) = un
405	IF( IsFail_slow ) Count_failure_slow(ig) = un
406
407	ELSE IF (.NOT. select_source_irrig) THEN
408
409	! For irrig. scheme, available_reserve gives the amount of water available for irrigation in every reservoir
410	! --> avail_reserve is a vector of dimension=3, BY DEFINITION i=1 for streamflow, i=2 fast, and i=3 slow reservoir
411
412	pot_slow_wdr_dummy = avail_reserve(3)*slow_reservoir_r(ig)
413	pot_fast_wdr_dummy = avail_reserve(2)*fast_reservoir_r(ig)
414	pot_stre_wdr_dummy = avail_reserve(1)*stream_reservoir_r(ig)
415	IsFail_slow = .FALSE. !
416	IsFail_fast = .FALSE. !
417	IsFail_stre = .FALSE. !
418
419	irrig_actual_r(ig) = MIN(irrig_needs_r(ig), pot_stre_wdr_dummy + pot_fast_wdr_dummy + pot_slow_wdr_dummy )
420
421	!Note for irrig_gw_source(ig): first we add the slow_reservoir volume. Then we substract the updated slow_reservoir. It should be the
422	!Volume used for irrigation that comes from GW
423	! Idem for irrig_fast_source and irrig_str_source
424
425	slow_wdr_dummy = slow_reservoir_r(ig)
426	slow_reservoir_r(ig) = MAX( (1-avail_reserve(3) )*slow_reservoir_r(ig), slow_reservoir_r(ig) &
427	+ MIN(zero, avail_reserve(2)*fast_reservoir_r(ig) &
428	+ MIN(zero, avail_reserve(1)*stream_reservoir_r(ig)-irrig_actual_r(ig))))
429	slow_wdr_dummy = slow_wdr_dummy - slow_reservoir_r(ig)
430	irrig_gw_source_r(ig) = irrig_gw_source_r(ig) + slow_wdr_dummy
431
432	fast_wdr_dummy = fast_reservoir_r(ig)
433	fast_reservoir_r(ig) = MAX( (1-avail_reserve(2) )*fast_reservoir_r(ig) , fast_reservoir_r(ig) &
434	+ MIN(zero, avail_reserve(1)*stream_reservoir_r(ig)-irrig_actual_r(ig)))
435	fast_wdr_dummy = fast_wdr_dummy - fast_reservoir_r(ig)
436	irrig_fast_source_r(ig) = irrig_fast_source_r(ig) + fast_wdr_dummy
437
438	stre_wdr_dummy = stream_reservoir_r(ig)
439	stream_reservoir_r(ig) = MAX( (1-avail_reserve(1) )*stream_reservoir_r(ig), stream_reservoir_r(ig)-irrig_actual_r(ig) )
440	stre_wdr_dummy = stre_wdr_dummy - stream_reservoir_r(ig)
441	irrig_str_source_r(ig) = irrig_str_source_r(ig) + stre_wdr_dummy
442
443	irrig_deficit_r(ig) = irrig_needs_r(ig)-irrig_actual_r(ig)
444
445	!A reservoir is failing to give water for infiltration if pot. req > pot. withdrawal
446	! Because it follows the old scheme, we do not separate between surface/gw, but consider that
447	! priority is given in this order: River, Fast and Slow reservoir.
448
449	IsFail_slow = ( ( irrig_needs_r(ig) - stre_wdr_dummy - fast_wdr_dummy ) > pot_slow_wdr_dummy )
450	IsFail_fast = ( ( irrig_needs_r(ig) - stre_wdr_dummy ) > pot_fast_wdr_dummy )
451	IsFail_stre = ( irrig_needs_r(ig) > pot_stre_wdr_dummy )
452
453	IF( IsFail_stre ) Count_failure_stre(ig) = un
454	IF( IsFail_fast ) Count_failure_fast(ig) = un
455	IF( IsFail_slow ) Count_failure_slow(ig) = un
456
457	ENDIF
458	ELSE
459	irrig_actual_r(ig) = 0
460	ENDIF
461	ENDDO
462
463	CALL routing_flow_set(slow_reservoir_r=slow_reservoir_r, fast_reservoir_r=fast_reservoir_r, stream_reservoir_r=stream_reservoir_r)
464
465	ENDIF
466
467	END SUBROUTINE irrigation_new
468
469
470	SUBROUTINE irrigation_old(irrig_needs_r, irrig_actual_r)
471	USE constantes
472	USE mod_orchidee_para
473	USE routing_native_flow_mod, ONLY : routing_mask_r, routing_flow_get, routing_flow_set
474	IMPLICIT NONE
475	REAL(r_std),INTENT(IN) :: irrig_needs_r(nbpt_r) !! Total irrigation requirement (water requirements by the crop for its optimal growth) (kg)
476	REAL(r_std),INTENT(OUT) :: irrig_actual_r(nbpt_r) !! Possible irrigation according to the water availability in the reservoirs (kg)
477
478	REAL(r_std) :: irrig_deficit_r(nbpt_r) !! Amount of water missing for irrigation (kg)
479
480	LOGICAL :: IsFail_slow !! Logical to ask if slow reserv. does not fit irrigation demand
481	LOGICAL :: IsFail_fast !! Logical to ask if fast reserv. does not fit irrigation demand
482	LOGICAL :: IsFail_stre !! Logical to ask if stream reserv. does not fit irrigation demand
483
484	REAL(r_std) :: Count_failure_slow(nbpt_r) !! Counter times slow reserv. does not fit irrigation demand
485	REAL(r_std) :: Count_failure_fast(nbpt_r) !! Counter times fast reserv. does not fit irrigation demand
486	REAL(r_std) :: Count_failure_stre(nbpt_r) !! Counter times stream reserv. does not fit irrigation demand
487
488	REAL(r_std) :: pot_slow_wdr_dummy, pot_fast_wdr_dummy, pot_stre_wdr_dummy
489	REAL(r_std) :: slow_wdr_dummy, fast_wdr_dummy, stre_wdr_dummy
490
491	REAL(r_std) :: slow_reservoir_r(nbpt_r)
492	REAL(r_std) :: fast_reservoir_r(nbpt_r)
493	REAL(r_std) :: stream_reservoir_r(nbpt_r)
494
495	INTEGER :: ig
496
497	IF (is_omp_root) THEN
498
499	Count_failure_slow(:) = zero !
500	Count_failure_fast(:) = zero !
501	Count_failure_stre(:) = zero !
502
503	CALL routing_flow_get(slow_reservoir_r=slow_reservoir_r, fast_reservoir_r=fast_reservoir_r, stream_reservoir_r=stream_reservoir_r)
504
505	DO ig=1,nbpt_r
506
507	IF (routing_mask_r(ig)) THEN
508
509	! Old irrigation scheme as in tag 2.0
510	! Note for irrig_gw_source(ig): first we add the slow_reservoir volume. Then we substract the updated slow_reservoir. It should be the
511	! Volume used for irrigation that comes from GW
512	! Idem for irrig_fast_source and irrig_str_source
513
514	pot_slow_wdr_dummy = slow_reservoir_r(ig)
515	pot_fast_wdr_dummy = fast_reservoir_r(ig)
516	pot_stre_wdr_dummy = stream_reservoir_r(ig)
517	IsFail_slow = .FALSE. !
518	IsFail_fast = .FALSE. !
519	IsFail_stre = .FALSE. !
520
521	irrig_actual_r(ig) = MIN(irrig_needs_r(ig), stream_reservoir_r(ig) + fast_reservoir_r(ig) + slow_reservoir_r(ig) )
522
523	slow_wdr_dummy = slow_reservoir_r(ig)
524	slow_reservoir_r(ig) = MAX(zero, slow_reservoir_r(ig) + MIN(zero, fast_reservoir_r(ig) &
525	+ MIN(zero, stream_reservoir_r(ig)-irrig_actual_r(ig))))
526	slow_wdr_dummy = slow_wdr_dummy - slow_reservoir_r(ig)
527	irrig_gw_source_r(ig) = irrig_gw_source_r(ig) + slow_wdr_dummy
528
529	fast_wdr_dummy = fast_reservoir_r(ig)
530	fast_reservoir_r(ig) = MAX( zero, fast_reservoir_r(ig) + MIN(zero, stream_reservoir_r(ig)-irrig_actual_r(ig)))
531	fast_wdr_dummy = fast_wdr_dummy - fast_reservoir_r(ig)
532	irrig_fast_source_r(ig) = irrig_fast_source_r(ig) + fast_wdr_dummy
533
534	stre_wdr_dummy = stream_reservoir_r(ig)
535	stream_reservoir_r(ig) = MAX(zero, stream_reservoir_r(ig)-irrig_actual_r(ig) )
536	stre_wdr_dummy = stre_wdr_dummy - stream_reservoir_r(ig)
537	irrig_str_source_r(ig) = irrig_str_source_r(ig) + stre_wdr_dummy
538
539	irrig_deficit_r(ig) = irrig_needs_r(ig)-irrig_actual_r(ig)
540
541	! A reservoir is failing to give water for infiltration if pot. req > pot. withdrawal
542	! Because it follows the old scheme, we do not separate between surface/gw, but consider that
543	! priority is given in this order: River, Fast and Slow reservoir.
544
545	IsFail_slow = ( ( irrig_needs_r(ig) - stre_wdr_dummy - fast_wdr_dummy ) > pot_slow_wdr_dummy )
546	IsFail_fast = ( ( irrig_needs_r(ig) - stre_wdr_dummy ) > pot_fast_wdr_dummy )
547	IsFail_stre = ( irrig_needs_r(ig) > pot_stre_wdr_dummy )
548
549	IF( IsFail_stre ) Count_failure_stre(ig) = un
550	IF( IsFail_fast ) Count_failure_fast(ig) = un
551	IF( IsFail_slow ) Count_failure_slow(ig) = un
552	ELSE
553	irrig_actual_r(ig) = 0
554	ENDIF
555	ENDDO
556
557	CALL routing_flow_set(slow_reservoir_r=slow_reservoir_r, fast_reservoir_r=fast_reservoir_r, stream_reservoir_r=stream_reservoir_r)
558
559	ENDIF
560
561	END SUBROUTINE irrigation_old
562
563
564	SUBROUTINE irrigation_compute_requested_new(dt_routing, contfrac, root_deficit, soiltile, irrig_frac, irrig_netereq_r)
565	USE constantes
566	USE constantes_soil
567	USE xios
568	USE mod_orchidee_para
569	USE grid, ONLY : area
570	IMPLICIT NONE
571	REAL(r_std),INTENT(IN) :: dt_routing
572	REAL(r_std),INTENT(IN) :: contfrac(nbpt)
573	REAL(r_std),INTENT(IN) :: root_deficit(nbpt)
574	REAL(r_std),INTENT(IN) :: soiltile(nbpt,nstm) !! Fraction of each soil tile within vegtot (0-1, unitless)
575	REAL(r_std),INTENT(IN) :: irrig_frac(nbpt) !! Irrig. fraction interpolated in routing, and saved to pass to slowproc if irrigated_soiltile = .TRUE.
576	REAL(r_std),INTENT(OUT) :: irrig_netereq_r(nbpt_r)
577
578	REAL(r_std) :: irrig_netereq(nbpt)
579	REAL(r_std) :: irrig_netereq_mpi(nbp_mpi)
580
581	INTEGER :: ig
582
583	irrig_netereq(:)=0
584
585	DO ig=1,nbpt
586	!It enters to the new irrigation module only if there is an irrigated fraction, if not irrig_netereq = zero for that cell
587	IF ( irrig_frac(ig) .GT. min_sechiba ) THEN
588
589	irrig_netereq(ig) = irrig_netereq(ig) + MIN( irrig_dosmax/3600*dt_routing, &
590	root_deficit(ig) ) * soiltile(ig, irrig_st) * vegtot_mean(ig)
591	! By definition, irrig_dosmax is in kg/m2 of soil tile/hour,dividing by 3600(seconds/hour) * DT_ROUTING !
592	! = kg/m2 of soil tile/(routing timestep)
593	! irrig_netereq(kg/m2 of grid cell / routing timstep ) is equal to
594	! root_deficit (kg/m2 of soil tile) * soiltile*vegtot (fraction of soil tile at cell level) = kg/m2 of grid cell
595
596	IF (.NOT. irrigated_soiltile .AND. ( soiltile(ig,irrig_st) .GT. min_sechiba ) .AND. (vegtot_mean(ig) .GT. min_sechiba) ) THEN
597	! Irrigated_soiltile asks if there is an independent soil tile for irrigated crops. If not,
598	! actual volume calculated for irrig_netereq assumed that the whole SOILTILE was irrigated, but in this case
599	! just a fraction of the irrig_st (irrigated soil tile, by default = 3) is actually irrigated,
600	! and irrig_netereq needs to be reduced by irrig_frac/( soiltile * vegtot ) (note that it is max = 1 thanks to irrig_frac calculation in l. 424)
601	! Demand(ST3)irrig_frac/(soiltile(3)vegtot) = irrig_netereq_In_ST3, then
602	!irrig_netereq_In_ST3 * (soiltile(3)vegtot) = irrig_netereq at grid scale = Demand(ST3)irrig_frac.
603	irrig_netereq(ig) = irrig_netereq(ig) * irrig_frac(ig) / ( soiltile(ig,irrig_st) * vegtot_mean(ig) )
604	!irrig_netereq = kg/m2 of grid cell
605	ENDIF
606	ENDIF
607	ENDDO
608
609	CALL gather_omp(irrig_netereqareacontfrac, irrig_netereq_mpi)
610
611	IF (is_omp_root) THEN
612	CALL xios_send_field("irrig_netereq",irrig_netereq_mpi) !! contfrac => ok ?
613	CALL xios_recv_field("irrig_netereq_r",irrig_netereq_r) ! ==> need conservative quantity interp
614	ENDIF
615
616	END SUBROUTINE irrigation_compute_requested_new
617
618
619	SUBROUTINE irrigation_compute_requested_old(contfrac, reinfiltration, irrigated, irrig_need_r)
620	USE constantes
621	USE xios
622	USE mod_orchidee_para
623	USE grid, ONLY : area
624	IMPLICIT NONE
625	REAL(r_std),INTENT(IN) :: contfrac(nbpt)
626	REAL(r_std),INTENT(IN) :: reinfiltration(nbpt)
627	REAL(r_std),INTENT(IN) :: irrigated(nbpt)
628	REAL(r_std),INTENT(OUT) :: irrig_need_r(nbpt_r)
629
630	REAL(r_std) :: irrig_netereq(nbpt)
631	REAL(r_std) :: irrig_netereq_mpi(nbp_mpi)
632
633	INTEGER :: ig, ir
634
635	irrig_netereq(:)=0
636
637	DO ig=1,nbpt
638	IF ((vegtot_mean(ig) .GT. min_sechiba) .AND. (humrel_mean(ig) .LT. un-min_sechiba) .AND. (runoff_mean(ig) .LT. min_sechiba) ) THEN
639	irrig_netereq(ig) = (irrigated(ig)/area(ig)) * MAX(zero, transpot_mean(ig) - (precip_mean(ig)+reinfiltration(ig))) ! ==> ok kg
640	ELSE
641	irrig_netereq(ig) = 0
642	ENDIF
643	ENDDO
644
645	CALL gather_omp(irrig_netereqareacontfrac, irrig_netereq_mpi)
646
647	IF (is_omp_root) THEN
648	CALL xios_send_field("irrig_netereq",irrig_netereq_mpi) ! contfrac ??
649	CALL xios_recv_field("irrig_netereq_r",irrig_need_r) ! ==> need conservative quantity interp
650	ENDIF
651
652	END SUBROUTINE irrigation_compute_requested_old
653
654	SUBROUTINE test
655	USE routing_native_flow_mod, ONLY : routing_mask_r
656
657	END SUBROUTINE TEST
658
659
660
661
662
663
664	!for now SUBROUTINE abduction
665	!for now REAL(r_std), DIMENSION(nbpt) :: irrig_adduct !! Amount of water carried over from other basins for irrigation (kg)
666	!for now
667	!for now DO ig=1,nbpt
668
669	!for now !
670	!for now ! Check if we cannot find the missing water in another basin of the same grid (stream reservoir only).
671	!for now ! If we find that then we create some adduction from that subbasin to the one where we need it for
672	!for now ! irrigation.
673	!for now !
674
675	!for now !> If crops water requirements have not been supplied (irrig_deficit>0), we check if we cannot find the missing water
676	!for now !> in another basin of the same grid. If there is water in the stream reservoir of this subbasin, we create some adduction
677	!for now !> from that subbasin to the one where we need it for irrigation.
678	!for now !>
679	!for now
680	!for now DO ib=1,nbasmax
681	!for now stream_tot = a_stream_adduction * SUM(stream_reservoir(ig,:))
682
683	!for now DO WHILE ( irrig_deficit(ig,ib) > min_sechiba .AND. stream_tot > min_sechiba)
684	!for now fi = MAXLOC(stream_reservoir(ig,:))
685	!for now ib2 = fi(1)
686
687	!for now irrig_adduct(ig,ib) = MIN(irrig_deficit(ig,ib), a_stream_adduction * stream_reservoir(ig,ib2))
688	!for now stream_reservoir(ig,ib2) = stream_reservoir(ig,ib2)-irrig_adduct(ig,ib)
689	!for now irrig_deficit(ig,ib) = irrig_deficit(ig,ib)-irrig_adduct(ig,ib)
690	!for now stream_tot = a_stream_adduction * SUM(stream_reservoir(ig,:))
691
692	!for now ENDDO
693	!for now ENDDO
694	!for now
695	!for now ENDDO
696
697	!for now !
698	!for now ! If we are at higher resolution we might need to look at neighboring grid boxes to find the streams
699	!for now ! which can feed irrigation
700	!for now !
701	!for now !> At higher resolution (grid box smaller than 100x100km), we can import water from neighboring grid boxes
702	!for now !> to the one where we need it for irrigation.
703	!for now !
704	!for now IF (is_root_prc) THEN
705	!for now ALLOCATE(irrig_deficit_glo(nbp_glo, nbasmax), stream_reservoir_glo(nbp_glo, nbasmax), &
706	!for now irrig_adduct_glo(nbp_glo, nbasmax), stat=ier)
707	!for now ELSE
708	!for now ALLOCATE(irrig_deficit_glo(0, 0), stream_reservoir_glo(0, 0), irrig_adduct_glo(0, 0), stat=ier)
709	!for now ENDIF
710	!for now IF (ier /= 0) CALL ipslerr_p(3,'routing_flow','Pb in allocate for irrig_deficit_glo, stream_reservoir_glo,...','','')
711
712	!for now CALL gather(irrig_deficit, irrig_deficit_glo)
713	!for now CALL gather(stream_reservoir, stream_reservoir_glo)
714	!for now CALL gather(irrig_adduct, irrig_adduct_glo)
715
716	!for now IF (is_root_prc) THEN
717	!for now !
718	!for now DO ig=1,nbp_glo
719	!for now ! Only work if the grid box is smaller than 100x100km. Else the piplines we build
720	!for now ! here would be too long to be reasonable.
721	!for now IF ( resolution_g(ig,1) < 100000. .AND. resolution_g(ig,2) < 100000. ) THEN
722	!for now
723	!for now DO ib=1,nbasmax
724	!for now !
725	!for now IF ( irrig_deficit_glo(ig,ib) > min_sechiba ) THEN
726	!for now !
727	!for now streams_around(:,:) = zero
728	!for now !
729	!for now DO in=1,NbNeighb
730	!for now ig2 = neighbours_g(ig,in)
731	!for now IF (ig2 .GT. 0 ) THEN
732	!for now streams_around(in,:) = a_stream_adduction * stream_reservoir_glo(ig2,:)
733	!for now igrd(in) = ig2
734	!for now ENDIF
735	!for now ENDDO
736	!for now !
737	!for now IF ( MAXVAL(streams_around) .GT. zero ) THEN
738	!for now !
739	!for now ff=MAXLOC(streams_around)
740	!for now ig2=igrd(ff(1))
741	!for now ib2=ff(2)
742	!for now !
743	!for now IF ( routing_area_glo(ig2,ib2) .GT. 0 .AND. a_stream_adduction * stream_reservoir_glo(ig2,ib2) > zero ) THEN
744	!for now adduction = MIN(irrig_deficit_glo(ig,ib), a_stream_adduction * stream_reservoir_glo(ig2,ib2))
745	!for now stream_reservoir_glo(ig2,ib2) = stream_reservoir_glo(ig2,ib2) - adduction
746	!for now irrig_deficit_glo(ig,ib) = irrig_deficit_glo(ig,ib) - adduction
747	!for now irrig_adduct_glo(ig,ib) = irrig_adduct_glo(ig,ib) + adduction
748	!for now ENDIF
749	!for now !
750	!for now ENDIF
751	!for now !
752	!for now ENDIF
753	!for now !
754	!for now ENDDO
755	!for now
756	!for now ENDIF
757
758	!for now ENDDO
759	!for now !
760	!for now ENDIF
761
762	!for now CALL scatter(irrig_deficit_glo, irrig_deficit)
763	!for now CALL scatter(stream_reservoir_glo, stream_reservoir)
764	!for now CALL scatter(irrig_adduct_glo, irrig_adduct)
765
766	!for now DEALLOCATE(irrig_deficit_glo, stream_reservoir_glo, irrig_adduct_glo)
767	!for now
768	!for now
769	!for now END SUBROUTINE abduction
770
771
772
773	END MODULE routing_native_irrig_mod

Note: See TracBrowser for help on using the repository browser.

Download in other formats: