Context Navigation

routing_native_irrig.f90 @ 8504

Last change on this file since 8504 was 8504, checked in by josefine.ghattas, 7 weeks ago
Integrated correction done in the trunk [8503] for compilation without XIOS.
File size: 34.3 KB

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

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

Context Navigation

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

Download in other formats: