Context Navigation

routing_native_irrig.f90 @ 8227

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

New version of routing simple => become routing_native ; after validation routing_simple will be supress

more robust routing scheme on native grid
rearchitecturing of routing by spliting routing sub-process into separate file/module => smaller file
works on global as well on regional (orchidee) grid.
works on standard 50km routing grid and highres routing MERIT grid (1-2km), for regional studies
integrate irrigation old scheme et new scheme (pedro arboleda)
more water conservative (1e-16 Vs 1e-7 relative for simple routing)
To be tested : ICOLMDZOR grid and ICOLAMDZ-LAM grid that are expeted to work also.

It is a first guess, not the definitive native routing package that will continue to evolve.
No side effect is expected on other configurations since it is just some files adding.
Native routing is activated using the run.def key :
ROUTING_METHOD=native

File size: 34.0 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, 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) :: reinfiltration(nbpt)
230	REAL(r_std), INTENT(IN) :: irrigated_next(nbpt)
231	REAL(r_std), INTENT(IN) :: irrig_frac(nbpt) !! Irrig. fraction interpolated in routing, and saved to pass to slowproc if irrigated_soiltile = .TRUE.
232	REAL(r_std), INTENT(IN) :: root_deficit(nbpt) !! soil water deficit
233	REAL(r_std), INTENT(IN) :: soiltile(nbpt,nstm) !! Fraction of each soil tile within vegtot (0-1, unitless)
234	REAL(r_std), INTENT(IN) :: fraction_aeirrig_sw(nbpt) !! Fraction of area equipped for irrigation from surface water, of irrig_frac
235
236	REAL(r_std) :: irrig_needs_r(nbpt_r) !! Total irrigation requirement (water requirements by the crop for its optimal growth) (kg)
237	REAL(r_std) :: irrig_actual_r(nbpt_r) !! Possible irrigation according to the water availability in the reservoirs (kg)
238	REAL(r_std) :: irrig_actual(nbpt) !! Possible irrigation according to the water availability in the reservoirs (kg)
239
240	IF (do_irrigation) THEN
241	IF (irrig_map_dynamic_flag ) irrigated=irrigated_next
242
243	IF (old_irrig_scheme) THEN
244	CALL irrigation_compute_requested_old(reinfiltration, irrigated, irrig_needs_r)
245	CALL irrigation_old(irrig_needs_r, irrig_actual_r)
246	ELSE
247	CALL irrigation_compute_requested_new(dt_routing, root_deficit, soiltile, irrig_frac, irrig_needs_r)
248	CALL irrigation_new(fraction_aeirrig_sw, irrig_needs_r, irrig_actual_r)
249	ENDIF
250
251	IF (is_omp_root) THEN
252	CALL xios_send_field("routing_irrigation_r",irrig_actual_r)
253	CALL xios_recv_field("routing_irrigation",irrig_actual)
254	ENDIF
255	CALL scatter_omp(irrig_actual,irrigation_mean)
256	irrigation_mean(:)=irrigation_mean(:)/area(:)
257	CALL irrigation_mean_reset()
258	ELSE
259	irrigation_mean(:) = 0
260	ENDIF
261
262	END SUBROUTINE irrigation_main
263
264
265
266
267
268	SUBROUTINE irrigation_new(fraction_aeirrig_sw, irrig_needs_r, irrig_actual_r)
269	USE constantes
270	USE routing_native_flow_mod, ONLY : routing_mask_r, routing_flow_get, routing_flow_set
271	USE mod_orchidee_para
272	USE xios
273	IMPLICIT NONE
274
275	REAL(r_std),INTENT(IN) :: fraction_aeirrig_sw(nbpt) !! Total irrigation requirement (water requirements by the crop for its optimal growth) (kg)
276	REAL(r_std),INTENT(IN) :: irrig_needs_r(nbpt_r) !! Total irrigation requirement (water requirements by the crop for its optimal growth) (kg)
277	REAL(r_std),INTENT(OUT) :: irrig_actual_r(nbpt_r) !! Possible irrigation according to the water availability in the reservoirs (kg)
278
279	REAL(r_std) :: irrig_deficit_r(nbpt_r) !! Amount of water missing for irrigation (kg)
280
281	LOGICAL :: IsFail_slow !! Logical to ask if slow reserv. does not fit irrigation demand
282	LOGICAL :: IsFail_fast !! Logical to ask if fast reserv. does not fit irrigation demand
283	LOGICAL :: IsFail_stre !! Logical to ask if stream reserv. does not fit irrigation demand
284
285	REAL(r_std) :: Count_failure_slow(nbpt_r) !! Counter times slow reserv. does not fit irrigation demand
286	REAL(r_std) :: Count_failure_fast(nbpt_r) !! Counter times fast reserv. does not fit irrigation demand
287	REAL(r_std) :: Count_failure_stre(nbpt_r) !! Counter times stream reserv. does not fit irrigation demand
288
289	REAL(r_std) :: pot_slow_wdr_dummy, pot_fast_wdr_dummy, pot_stre_wdr_dummy
290	REAL(r_std) :: slow_wdr_dummy, fast_wdr_dummy, stre_wdr_dummy
291
292	REAL(r_std) :: slow_reservoir_r(nbpt_r)
293	REAL(r_std) :: fast_reservoir_r(nbpt_r)
294	REAL(r_std) :: stream_reservoir_r(nbpt_r)
295	REAL(r_std) :: fraction_aeirrig_sw_mpi(nbp_mpi)
296	REAL(r_std) :: fraction_aeirrig_sw_r(nbpt_r)
297	REAL(r_std) :: pcent_vol_irrig
298	INTEGER :: ig
299
300
301	CALL gather_omp(fraction_aeirrig_sw, fraction_aeirrig_sw_mpi)
302
303	IF (is_omp_root) THEN
304
305	CALL xios_send_field("fraction_aeirrig_sw", fraction_aeirrig_sw_mpi)
306	CALL xios_recv_field("fraction_aeirrig_sw_r", fraction_aeirrig_sw_r) ! ==> need conservative quantity interp
307
308
309	Count_failure_slow(:) = zero !
310	Count_failure_fast(:) = zero !
311	Count_failure_stre(:) = zero !
312
313	CALL routing_flow_get(slow_reservoir_r=slow_reservoir_r, fast_reservoir_r=fast_reservoir_r, stream_reservoir_r=stream_reservoir_r)
314
315	DO ig=1,nbpt_r
316
317	IF (routing_mask_r(ig)) THEN
318
319	IF (select_source_irrig) THEN
320
321	! For irrig. scheme, available_reserve gives the amount of water available for irrigation in every reservoir
322	! --> avail_reserve is a vector of dimension=3, BY DEFINITION i=1 for streamflow, i=2 fast, and i=3 slow reservoir
323
324	! The new priorization scheme takes into account irrig. infrastructur according to GMIA map
325	! It also withdraw water according to availability, it means that it wont seek for all the water in then
326	! stream reservoir, even if this one could respond to the demand by itself
327
328	pot_slow_wdr_dummy = ( 1 - fraction_aeirrig_sw_r(ig)) * avail_reserve(3)*slow_reservoir_r(ig)
329	pot_fast_wdr_dummy = fraction_aeirrig_sw_r(ig) * avail_reserve(2) * fast_reservoir_r(ig)
330	pot_stre_wdr_dummy = fraction_aeirrig_sw_r(ig) * avail_reserve(1) * stream_reservoir_r(ig)
331	pcent_vol_irrig = zero
332	IsFail_slow = .FALSE. !
333	IsFail_fast = .FALSE. !
334	IsFail_stre = .FALSE. !
335
336	irrig_actual_r(ig) = MIN(irrig_needs_r(ig), pot_stre_wdr_dummy + pot_fast_wdr_dummy + pot_slow_wdr_dummy)
337
338	!! additional IF to calculate pcent_vol_irrig, in the case the total avail.
339	!! water is zero, I.E. when there is no water in surface and fraction_ae = 0,
340	!! so GW is not taken into account
341	!! Note on pcent_vol_irrig: It correspond to the fraction of available water in surface,
342	!! considering environmental needs and irrigation equipement by source from map. It controls
343	!! how the source of water withdrawl, especially when requirements < available water
344
345	IF ( (pot_stre_wdr_dummy + pot_fast_wdr_dummy + pot_slow_wdr_dummy) .GT. min_sechiba ) THEN
346
347	pcent_vol_irrig = ( pot_stre_wdr_dummy + pot_fast_wdr_dummy ) / &
348	( pot_stre_wdr_dummy + pot_fast_wdr_dummy + pot_slow_wdr_dummy )
349
350	!Irrig_actual set to zero, because there is no available water.
351	!Put to avoid negative values due to problems in the Min function
352
353	irrig_actual_r(ig) = MAX(irrig_actual_r(ig), zero)
354
355	!Already zero because pcent_vol_irrig initialized to zero
356	!Put here to readability but not necessary
357	!ELSE
358	! pcent_vol_irrig = zero
359
360	ENDIF
361
362	!Note for irrig_gw_source(ig): first we add the slow_reservoir volume. Then we substract the updated slow_reservoir. It should be the
363	!Volume used for irrigation that comes from GW
364	! Idem for irrig_fast_source and irrig_str_source
365
366	slow_wdr_dummy = slow_reservoir_r(ig)
367	slow_reservoir_r(ig) = MAX( (un - ( un - fraction_aeirrig_sw_r(ig) ) * avail_reserve(3) ) * &
368	slow_reservoir_r(ig), slow_reservoir_r(ig) + &
369	MIN( - irrig_actual_r(ig) * (un - pcent_vol_irrig ), &
370	avail_reserve(2) * fraction_aeirrig_sw_r(ig) * fast_reservoir_r(ig) + &
371	MIN(zero, avail_reserve(1) * fraction_aeirrig_sw_r(ig) * stream_reservoir_r(ig) - &
372	pcent_vol_irrig * irrig_actual_r(ig) ) ) )
373
374	slow_wdr_dummy = slow_wdr_dummy - slow_reservoir_r(ig)
375	irrig_gw_source_r(ig) = irrig_gw_source_r(ig) + slow_wdr_dummy
376
377	fast_wdr_dummy = fast_reservoir_r(ig)
378	fast_reservoir_r(ig) = MAX( (un - avail_reserve(2) * fraction_aeirrig_sw_r(ig) ) * fast_reservoir_r(ig) , &
379	fast_reservoir_r(ig) + MIN(zero, avail_reserve(1) * fraction_aeirrig_sw_r(ig) * stream_reservoir_r(ig) - &
380	pcent_vol_irrig * irrig_actual_r(ig) ) )
381	fast_wdr_dummy = fast_wdr_dummy - fast_reservoir_r(ig)
382	irrig_fast_source_r(ig) = irrig_fast_source_r(ig) + fast_wdr_dummy
383
384	stre_wdr_dummy = stream_reservoir_r(ig)
385	stream_reservoir_r(ig) = MAX((un - avail_reserve(1)* fraction_aeirrig_sw_r(ig) )*stream_reservoir_r(ig), &
386	stream_reservoir_r(ig) - pcent_vol_irrig * irrig_actual_r(ig) )
387	stre_wdr_dummy = stre_wdr_dummy - stream_reservoir_r(ig)
388	irrig_str_source_r(ig) = irrig_str_source_r(ig) + stre_wdr_dummy
389
390	irrig_deficit_r(ig) = irrig_needs_r(ig)-irrig_actual_r(ig)
391
392	!A reservoir is failing to give water for infiltration if pot. req > pot. withdrawal
393	!We assume that the pot. requirement = Needs * fraction of area equipped for SW/GW
394	!In the case of surface. we also sustract the withdrawal from Fast/Stream, because both are
395	! considered as surface water
396
397	IsFail_slow = ( ( irrig_needs_r(ig)*(un - fraction_aeirrig_sw_r(ig)) ) > pot_slow_wdr_dummy )
398	IsFail_fast = ( ( irrig_needs_r(ig)*fraction_aeirrig_sw_r(ig) - stre_wdr_dummy ) > pot_fast_wdr_dummy )
399	IsFail_stre = ( ( irrig_needs_r(ig)*fraction_aeirrig_sw_r(ig) - fast_wdr_dummy ) > pot_stre_wdr_dummy )
400
401	IF( IsFail_stre ) Count_failure_stre(ig) = un
402	IF( IsFail_fast ) Count_failure_fast(ig) = un
403	IF( IsFail_slow ) Count_failure_slow(ig) = un
404
405	ELSE IF (.NOT. select_source_irrig) THEN
406
407	! For irrig. scheme, available_reserve gives the amount of water available for irrigation in every reservoir
408	! --> avail_reserve is a vector of dimension=3, BY DEFINITION i=1 for streamflow, i=2 fast, and i=3 slow reservoir
409
410	pot_slow_wdr_dummy = avail_reserve(3)*slow_reservoir_r(ig)
411	pot_fast_wdr_dummy = avail_reserve(2)*fast_reservoir_r(ig)
412	pot_stre_wdr_dummy = avail_reserve(1)*stream_reservoir_r(ig)
413	IsFail_slow = .FALSE. !
414	IsFail_fast = .FALSE. !
415	IsFail_stre = .FALSE. !
416
417	irrig_actual_r(ig) = MIN(irrig_needs_r(ig), pot_stre_wdr_dummy + pot_fast_wdr_dummy + pot_slow_wdr_dummy )
418
419	!Note for irrig_gw_source(ig): first we add the slow_reservoir volume. Then we substract the updated slow_reservoir. It should be the
420	!Volume used for irrigation that comes from GW
421	! Idem for irrig_fast_source and irrig_str_source
422
423	slow_wdr_dummy = slow_reservoir_r(ig)
424	slow_reservoir_r(ig) = MAX( (1-avail_reserve(3) )*slow_reservoir_r(ig), slow_reservoir_r(ig) &
425	+ MIN(zero, avail_reserve(2)*fast_reservoir_r(ig) &
426	+ MIN(zero, avail_reserve(1)*stream_reservoir_r(ig)-irrig_actual_r(ig))))
427	slow_wdr_dummy = slow_wdr_dummy - slow_reservoir_r(ig)
428	irrig_gw_source_r(ig) = irrig_gw_source_r(ig) + slow_wdr_dummy
429
430	fast_wdr_dummy = fast_reservoir_r(ig)
431	fast_reservoir_r(ig) = MAX( (1-avail_reserve(2) )*fast_reservoir_r(ig) , fast_reservoir_r(ig) &
432	+ MIN(zero, avail_reserve(1)*stream_reservoir_r(ig)-irrig_actual_r(ig)))
433	fast_wdr_dummy = fast_wdr_dummy - fast_reservoir_r(ig)
434	irrig_fast_source_r(ig) = irrig_fast_source_r(ig) + fast_wdr_dummy
435
436	stre_wdr_dummy = stream_reservoir_r(ig)
437	stream_reservoir_r(ig) = MAX( (1-avail_reserve(1) )*stream_reservoir_r(ig), stream_reservoir_r(ig)-irrig_actual_r(ig) )
438	stre_wdr_dummy = stre_wdr_dummy - stream_reservoir_r(ig)
439	irrig_str_source_r(ig) = irrig_str_source_r(ig) + stre_wdr_dummy
440
441	irrig_deficit_r(ig) = irrig_needs_r(ig)-irrig_actual_r(ig)
442
443	!A reservoir is failing to give water for infiltration if pot. req > pot. withdrawal
444	! Because it follows the old scheme, we do not separate between surface/gw, but consider that
445	! priority is given in this order: River, Fast and Slow reservoir.
446
447	IsFail_slow = ( ( irrig_needs_r(ig) - stre_wdr_dummy - fast_wdr_dummy ) > pot_slow_wdr_dummy )
448	IsFail_fast = ( ( irrig_needs_r(ig) - stre_wdr_dummy ) > pot_fast_wdr_dummy )
449	IsFail_stre = ( irrig_needs_r(ig) > pot_stre_wdr_dummy )
450
451	IF( IsFail_stre ) Count_failure_stre(ig) = un
452	IF( IsFail_fast ) Count_failure_fast(ig) = un
453	IF( IsFail_slow ) Count_failure_slow(ig) = un
454
455	ENDIF
456
457	ENDIF
458	ENDDO
459
460	CALL routing_flow_set(slow_reservoir_r=slow_reservoir_r, fast_reservoir_r=fast_reservoir_r, stream_reservoir_r=stream_reservoir_r)
461
462	ENDIF
463
464	END SUBROUTINE irrigation_new
465
466
467	SUBROUTINE irrigation_old(irrig_needs_r, irrig_actual_r)
468	USE constantes
469	USE mod_orchidee_para
470	USE routing_native_flow_mod, ONLY : routing_mask_r, routing_flow_get, routing_flow_set
471	IMPLICIT NONE
472	REAL(r_std),INTENT(IN) :: irrig_needs_r(nbpt_r) !! Total irrigation requirement (water requirements by the crop for its optimal growth) (kg)
473	REAL(r_std),INTENT(OUT) :: irrig_actual_r(nbpt_r) !! Possible irrigation according to the water availability in the reservoirs (kg)
474
475	REAL(r_std) :: irrig_deficit_r(nbpt_r) !! Amount of water missing for irrigation (kg)
476
477	LOGICAL :: IsFail_slow !! Logical to ask if slow reserv. does not fit irrigation demand
478	LOGICAL :: IsFail_fast !! Logical to ask if fast reserv. does not fit irrigation demand
479	LOGICAL :: IsFail_stre !! Logical to ask if stream reserv. does not fit irrigation demand
480
481	REAL(r_std) :: Count_failure_slow(nbpt_r) !! Counter times slow reserv. does not fit irrigation demand
482	REAL(r_std) :: Count_failure_fast(nbpt_r) !! Counter times fast reserv. does not fit irrigation demand
483	REAL(r_std) :: Count_failure_stre(nbpt_r) !! Counter times stream reserv. does not fit irrigation demand
484
485	REAL(r_std) :: pot_slow_wdr_dummy, pot_fast_wdr_dummy, pot_stre_wdr_dummy
486	REAL(r_std) :: slow_wdr_dummy, fast_wdr_dummy, stre_wdr_dummy
487
488	REAL(r_std) :: slow_reservoir_r(nbpt_r)
489	REAL(r_std) :: fast_reservoir_r(nbpt_r)
490	REAL(r_std) :: stream_reservoir_r(nbpt_r)
491
492	INTEGER :: ig
493
494	IF (is_omp_root) THEN
495
496	Count_failure_slow(:) = zero !
497	Count_failure_fast(:) = zero !
498	Count_failure_stre(:) = zero !
499
500	CALL routing_flow_get(slow_reservoir_r=slow_reservoir_r, fast_reservoir_r=fast_reservoir_r, stream_reservoir_r=stream_reservoir_r)
501
502	DO ig=1,nbpt_r
503
504	IF (routing_mask_r(ig)) THEN
505
506	! Old irrigation scheme as in tag 2.0
507	! Note for irrig_gw_source(ig): first we add the slow_reservoir volume. Then we substract the updated slow_reservoir. It should be the
508	! Volume used for irrigation that comes from GW
509	! Idem for irrig_fast_source and irrig_str_source
510
511	pot_slow_wdr_dummy = slow_reservoir_r(ig)
512	pot_fast_wdr_dummy = fast_reservoir_r(ig)
513	pot_stre_wdr_dummy = stream_reservoir_r(ig)
514	IsFail_slow = .FALSE. !
515	IsFail_fast = .FALSE. !
516	IsFail_stre = .FALSE. !
517
518	irrig_actual_r(ig) = MIN(irrig_needs_r(ig), stream_reservoir_r(ig) + fast_reservoir_r(ig) + slow_reservoir_r(ig) )
519
520	slow_wdr_dummy = slow_reservoir_r(ig)
521	slow_reservoir_r(ig) = MAX(zero, slow_reservoir_r(ig) + MIN(zero, fast_reservoir_r(ig) &
522	+ MIN(zero, stream_reservoir_r(ig)-irrig_actual_r(ig))))
523	slow_wdr_dummy = slow_wdr_dummy - slow_reservoir_r(ig)
524	irrig_gw_source_r(ig) = irrig_gw_source_r(ig) + slow_wdr_dummy
525
526	fast_wdr_dummy = fast_reservoir_r(ig)
527	fast_reservoir_r(ig) = MAX( zero, fast_reservoir_r(ig) + MIN(zero, stream_reservoir_r(ig)-irrig_actual_r(ig)))
528	fast_wdr_dummy = fast_wdr_dummy - fast_reservoir_r(ig)
529	irrig_fast_source_r(ig) = irrig_fast_source_r(ig) + fast_wdr_dummy
530
531	stre_wdr_dummy = stream_reservoir_r(ig)
532	stream_reservoir_r(ig) = MAX(zero, stream_reservoir_r(ig)-irrig_actual_r(ig) )
533	stre_wdr_dummy = stre_wdr_dummy - stream_reservoir_r(ig)
534	irrig_str_source_r(ig) = irrig_str_source_r(ig) + stre_wdr_dummy
535
536	irrig_deficit_r(ig) = irrig_needs_r(ig)-irrig_actual_r(ig)
537
538	! A reservoir is failing to give water for infiltration if pot. req > pot. withdrawal
539	! Because it follows the old scheme, we do not separate between surface/gw, but consider that
540	! priority is given in this order: River, Fast and Slow reservoir.
541
542	IsFail_slow = ( ( irrig_needs_r(ig) - stre_wdr_dummy - fast_wdr_dummy ) > pot_slow_wdr_dummy )
543	IsFail_fast = ( ( irrig_needs_r(ig) - stre_wdr_dummy ) > pot_fast_wdr_dummy )
544	IsFail_stre = ( irrig_needs_r(ig) > pot_stre_wdr_dummy )
545
546	IF( IsFail_stre ) Count_failure_stre(ig) = un
547	IF( IsFail_fast ) Count_failure_fast(ig) = un
548	IF( IsFail_slow ) Count_failure_slow(ig) = un
549	ELSE
550	irrig_actual_r(ig) = 0
551	ENDIF
552	ENDDO
553
554	CALL routing_flow_set(slow_reservoir_r=slow_reservoir_r, fast_reservoir_r=fast_reservoir_r, stream_reservoir_r=stream_reservoir_r)
555
556	ENDIF
557
558	END SUBROUTINE irrigation_old
559
560
561	SUBROUTINE irrigation_compute_requested_new(dt_routing, root_deficit, soiltile, irrig_frac, irrig_netereq_r)
562	USE constantes
563	USE constantes_soil
564	USE xios
565	USE mod_orchidee_para
566	USE grid, ONLY : area
567	IMPLICIT NONE
568	REAL(r_std),INTENT(IN) :: dt_routing
569	REAL(r_std),INTENT(IN) :: root_deficit(nbpt)
570	REAL(r_std),INTENT(IN) :: soiltile(nbpt,nstm) !! Fraction of each soil tile within vegtot (0-1, unitless)
571	REAL(r_std),INTENT(IN) :: irrig_frac(nbpt) !! Irrig. fraction interpolated in routing, and saved to pass to slowproc if irrigated_soiltile = .TRUE.
572	REAL(r_std),INTENT(OUT) :: irrig_netereq_r(nbpt_r)
573
574	REAL(r_std) :: irrig_netereq(nbpt)
575	REAL(r_std) :: irrig_netereq_mpi(nbp_mpi)
576	REAL(r_std) :: area_mpi(nbp_mpi)
577
578	INTEGER :: ig
579
580	irrig_netereq(:)=0
581
582	DO ig=1,nbpt
583	!It enters to the new irrigation module only if there is an irrigated fraction, if not irrig_netereq = zero for that cell
584	IF ( irrig_frac(ig) .GT. min_sechiba ) THEN
585
586	irrig_netereq(ig) = irrig_netereq(ig) + MIN( irrig_dosmax/3600*dt_routing, &
587	root_deficit(ig) ) * soiltile(ig, irrig_st) * vegtot_mean(ig)
588	! By definition, irrig_dosmax is in kg/m2 of soil tile/hour,dividing by 3600(seconds/hour) * DT_ROUTING !
589	! = kg/m2 of soil tile/(routing timestep)
590	! irrig_netereq(kg/m2 of grid cell / routing timstep ) is equal to
591	! root_deficit (kg/m2 of soil tile) * soiltile*vegtot (fraction of soil tile at cell level) = kg/m2 of grid cell
592
593	IF (.NOT. irrigated_soiltile .AND. ( soiltile(ig,irrig_st) .GT. min_sechiba ) .AND. (vegtot_mean(ig) .GT. min_sechiba) ) THEN
594	! Irrigated_soiltile asks if there is an independent soil tile for irrigated crops. If not,
595	! actual volume calculated for irrig_netereq assumed that the whole SOILTILE was irrigated, but in this case
596	! just a fraction of the irrig_st (irrigated soil tile, by default = 3) is actually irrigated,
597	! 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)
598	! Demand(ST3)irrig_frac/(soiltile(3)vegtot) = irrig_netereq_In_ST3, then
599	!irrig_netereq_In_ST3 * (soiltile(3)vegtot) = irrig_netereq at grid scale = Demand(ST3)irrig_frac.
600	irrig_netereq(ig) = irrig_netereq(ig) * irrig_frac(ig) / ( soiltile(ig,irrig_st) * vegtot_mean(ig) )
601	!irrig_netereq = kg/m2 of grid cell
602	ENDIF
603	ENDIF
604	ENDDO
605
606	CALL gather_omp(irrig_netereq, irrig_netereq_mpi)
607	CALL gather_omp(area, area_mpi)
608
609	IF (is_omp_root) THEN
610	CALL xios_send_field("irrig_netereq",irrig_netereq_mpi*area_mpi) !! contfrac ?
611	CALL xios_recv_field("irrig_netereq_r",irrig_netereq_r) ! ==> need conservative quantity interp
612	ENDIF
613
614	END SUBROUTINE irrigation_compute_requested_new
615
616
617	SUBROUTINE irrigation_compute_requested_old(reinfiltration, irrigated, irrig_need_r)
618	USE constantes
619	USE xios
620	USE mod_orchidee_para
621	USE grid, ONLY : area
622	IMPLICIT NONE
623	REAL(r_std),INTENT(IN) :: reinfiltration(nbpt)
624	REAL(r_std),INTENT(IN) :: irrigated(nbpt)
625	REAL(r_std),INTENT(OUT) :: irrig_need_r(nbpt_r)
626
627	REAL(r_std) :: irrig_netereq(nbpt)
628	REAL(r_std) :: irrig_netereq_mpi(nbp_mpi)
629	REAL(r_std) :: area_mpi(nbp_mpi)
630
631	INTEGER :: ig, ir
632
633	irrig_netereq(:)=0
634
635	DO ig=1,nbpt
636	IF ((vegtot_mean(ig) .GT. min_sechiba) .AND. (humrel_mean(ig) .LT. un-min_sechiba) .AND. (runoff_mean(ig) .LT. min_sechiba) ) THEN
637	irrig_netereq(ig) = (irrigated(ig)/area(ig)) * MAX(zero, transpot_mean(ig) - (precip_mean(ig)+reinfiltration(ig))) ! ==> ok kg
638	ELSE
639	irrig_netereq(ig) = 0
640	ENDIF
641	ENDDO
642
643	CALL gather_omp(irrig_netereq, irrig_netereq_mpi)
644	CALL gather_omp(area, area_mpi)
645
646	IF (is_omp_root) THEN
647	CALL xios_send_field("irrig_netereq",irrig_netereq_mpi*area_mpi) ! contfrac ??
648	CALL xios_recv_field("irrig_netereq_r",irrig_need_r) ! ==> need conservative quantity interp
649	ENDIF
650
651	END SUBROUTINE irrigation_compute_requested_old
652
653	SUBROUTINE test
654	USE routing_native_flow_mod, ONLY : routing_mask_r
655
656	END SUBROUTINE TEST
657
658
659
660
661
662
663	!for now SUBROUTINE abduction
664	!for now REAL(r_std), DIMENSION(nbpt) :: irrig_adduct !! Amount of water carried over from other basins for irrigation (kg)
665	!for now
666	!for now DO ig=1,nbpt
667
668	!for now !
669	!for now ! Check if we cannot find the missing water in another basin of the same grid (stream reservoir only).
670	!for now ! If we find that then we create some adduction from that subbasin to the one where we need it for
671	!for now ! irrigation.
672	!for now !
673
674	!for now !> If crops water requirements have not been supplied (irrig_deficit>0), we check if we cannot find the missing water
675	!for now !> in another basin of the same grid. If there is water in the stream reservoir of this subbasin, we create some adduction
676	!for now !> from that subbasin to the one where we need it for irrigation.
677	!for now !>
678	!for now
679	!for now DO ib=1,nbasmax
680	!for now stream_tot = a_stream_adduction * SUM(stream_reservoir(ig,:))
681
682	!for now DO WHILE ( irrig_deficit(ig,ib) > min_sechiba .AND. stream_tot > min_sechiba)
683	!for now fi = MAXLOC(stream_reservoir(ig,:))
684	!for now ib2 = fi(1)
685
686	!for now irrig_adduct(ig,ib) = MIN(irrig_deficit(ig,ib), a_stream_adduction * stream_reservoir(ig,ib2))
687	!for now stream_reservoir(ig,ib2) = stream_reservoir(ig,ib2)-irrig_adduct(ig,ib)
688	!for now irrig_deficit(ig,ib) = irrig_deficit(ig,ib)-irrig_adduct(ig,ib)
689	!for now stream_tot = a_stream_adduction * SUM(stream_reservoir(ig,:))
690
691	!for now ENDDO
692	!for now ENDDO
693	!for now
694	!for now ENDDO
695
696	!for now !
697	!for now ! If we are at higher resolution we might need to look at neighboring grid boxes to find the streams
698	!for now ! which can feed irrigation
699	!for now !
700	!for now !> At higher resolution (grid box smaller than 100x100km), we can import water from neighboring grid boxes
701	!for now !> to the one where we need it for irrigation.
702	!for now !
703	!for now IF (is_root_prc) THEN
704	!for now ALLOCATE(irrig_deficit_glo(nbp_glo, nbasmax), stream_reservoir_glo(nbp_glo, nbasmax), &
705	!for now irrig_adduct_glo(nbp_glo, nbasmax), stat=ier)
706	!for now ELSE
707	!for now ALLOCATE(irrig_deficit_glo(0, 0), stream_reservoir_glo(0, 0), irrig_adduct_glo(0, 0), stat=ier)
708	!for now ENDIF
709	!for now IF (ier /= 0) CALL ipslerr_p(3,'routing_flow','Pb in allocate for irrig_deficit_glo, stream_reservoir_glo,...','','')
710
711	!for now CALL gather(irrig_deficit, irrig_deficit_glo)
712	!for now CALL gather(stream_reservoir, stream_reservoir_glo)
713	!for now CALL gather(irrig_adduct, irrig_adduct_glo)
714
715	!for now IF (is_root_prc) THEN
716	!for now !
717	!for now DO ig=1,nbp_glo
718	!for now ! Only work if the grid box is smaller than 100x100km. Else the piplines we build
719	!for now ! here would be too long to be reasonable.
720	!for now IF ( resolution_g(ig,1) < 100000. .AND. resolution_g(ig,2) < 100000. ) THEN
721	!for now
722	!for now DO ib=1,nbasmax
723	!for now !
724	!for now IF ( irrig_deficit_glo(ig,ib) > min_sechiba ) THEN
725	!for now !
726	!for now streams_around(:,:) = zero
727	!for now !
728	!for now DO in=1,NbNeighb
729	!for now ig2 = neighbours_g(ig,in)
730	!for now IF (ig2 .GT. 0 ) THEN
731	!for now streams_around(in,:) = a_stream_adduction * stream_reservoir_glo(ig2,:)
732	!for now igrd(in) = ig2
733	!for now ENDIF
734	!for now ENDDO
735	!for now !
736	!for now IF ( MAXVAL(streams_around) .GT. zero ) THEN
737	!for now !
738	!for now ff=MAXLOC(streams_around)
739	!for now ig2=igrd(ff(1))
740	!for now ib2=ff(2)
741	!for now !
742	!for now IF ( routing_area_glo(ig2,ib2) .GT. 0 .AND. a_stream_adduction * stream_reservoir_glo(ig2,ib2) > zero ) THEN
743	!for now adduction = MIN(irrig_deficit_glo(ig,ib), a_stream_adduction * stream_reservoir_glo(ig2,ib2))
744	!for now stream_reservoir_glo(ig2,ib2) = stream_reservoir_glo(ig2,ib2) - adduction
745	!for now irrig_deficit_glo(ig,ib) = irrig_deficit_glo(ig,ib) - adduction
746	!for now irrig_adduct_glo(ig,ib) = irrig_adduct_glo(ig,ib) + adduction
747	!for now ENDIF
748	!for now !
749	!for now ENDIF
750	!for now !
751	!for now ENDIF
752	!for now !
753	!for now ENDDO
754	!for now
755	!for now ENDIF
756
757	!for now ENDDO
758	!for now !
759	!for now ENDIF
760
761	!for now CALL scatter(irrig_deficit_glo, irrig_deficit)
762	!for now CALL scatter(stream_reservoir_glo, stream_reservoir)
763	!for now CALL scatter(irrig_adduct_glo, irrig_adduct)
764
765	!for now DEALLOCATE(irrig_deficit_glo, stream_reservoir_glo, irrig_adduct_glo)
766	!for now
767	!for now
768	!for now END SUBROUTINE abduction
769
770
771
772	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 @ 8227

Download in other formats: