Context Navigation

routing_native_flow.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: 74.1 KB

Line
1	MODULE routing_native_flow_mod
2
3	USE ioipsl
4	USE xios_orchidee
5	USE ioipsl_para
6	USE constantes
7	USE constantes_soil
8	USE pft_parameters
9	USE sechiba_io_p
10	USE interpol_help
11	USE grid
12	USE mod_orchidee_para
13
14
15	IMPLICIT NONE
16	PRIVATE
17
18	PUBLIC :: routing_flow_xios_initialize, routing_flow_set, routing_flow_get, routing_flow_main
19	PUBLIC :: routing_flow_initialize, routing_flow_finalize, routing_flow_clear, routing_flow_make_mean
20	PUBLIC :: compute_coastline
21
22	INTEGER,SAVE :: nbpt
23	!$OMP THREADPRIVATE(nbpt)
24
25	!! PARAMETERS
26	REAL(r_std), SAVE :: fast_tcst = 3.0 !! Property of the fast reservoir (day/m)
27	!$OMP THREADPRIVATE(fast_tcst)
28	REAL(r_std), SAVE :: slow_tcst = 25.0 !! Property of the slow reservoir (day/m)
29	!$OMP THREADPRIVATE(slow_tcst)
30	REAL(r_std), SAVE :: stream_tcst = 0.24 !! Property of the stream reservoir (day/m)
31	!$OMP THREADPRIVATE(stream_tcst)
32
33
34	REAL(r_std),SAVE,ALLOCATABLE :: runoff_mean(:)
35	!$OMP THREADPRIVATE(runoff_mean)
36
37	REAL(r_std),SAVE,ALLOCATABLE :: drainage_mean(:)
38	!$OMP THREADPRIVATE(drainage_mean)
39
40	REAL(r_std),SAVE,ALLOCATABLE,PUBLIC,PROTECTED,TARGET :: fast_reservoir_r(:) !! Water amount in the fast reservoir (kg) - (local routing grid)
41	!$OMP THREADPRIVATE(fast_reservoir_r)
42
43	REAL(r_std),SAVE,ALLOCATABLE,PUBLIC,PROTECTED,TARGET:: slow_reservoir_r(:) !! Water amount in the slow reservoir (kg) - (local routing grid)
44	!$OMP THREADPRIVATE(slow_reservoir_r)
45
46	REAL(r_std),SAVE,ALLOCATABLE,PUBLIC,PROTECTED,TARGET :: stream_reservoir_r(:) !! Water amount in the stream reservoir (kg) - (local routing grid)
47	!$OMP THREADPRIVATE(stream_reservoir_r)
48
49	REAL(r_std),SAVE,ALLOCATABLE,PUBLIC,PROTECTED :: riverflow_mean(:) !! Water amount in the stream reservoir (kg) - (local routing grid)
50	!$OMP THREADPRIVATE(riverflow_mean)
51
52	REAL(r_std),SAVE,ALLOCATABLE,PUBLIC,PROTECTED :: coastalflow_mean(:) !! Water amount in the stream reservoir (kg) - (local routing grid)
53	!$OMP THREADPRIVATE(coastalflow_mean)
54
55	REAL(r_std),SAVE,ALLOCATABLE,PUBLIC,PROTECTED :: lakeinflow_mean(:) !! Water amount in the stream reservoir (kg) - (local routing grid)
56	!$OMP THREADPRIVATE(lakeinflow_mean)
57
58	!
59	! Specific variables for simple routing
60	!
61	REAL(r_std),SAVE,ALLOCATABLE :: topoind_r(:) !! Topographic index of the retention time (m) index - (local routing grid)
62	!$OMP THREADPRIVATE(topoind_r)
63	INTEGER,SAVE,ALLOCATABLE :: route_flow_rp1(:) !! flow index from cell to neighboring cell following the trip direction - (local routing grid + halo)
64	!$OMP THREADPRIVATE(route_flow_rp1)
65	LOGICAL,SAVE,ALLOCATABLE :: is_lakeinflow_r(:) !! is lake inflow point - (local routing grid)
66	!$OMP THREADPRIVATE(is_lakeinflow_r)
67	LOGICAL,SAVE,ALLOCATABLE :: is_coastalflow_r(:) !! is coastal flow point - (local routing grid)
68	!$OMP THREADPRIVATE(is_coastalflow_r)
69	LOGICAL,SAVE,ALLOCATABLE :: is_riverflow_r(:) !! is river flow point - (local routing grid)
70	!$OMP THREADPRIVATE(is_riverflow_r)
71	LOGICAL,SAVE,ALLOCATABLE :: is_coastline(:) !! is coastline point - (local native grid)
72	!$OMP THREADPRIVATE(is_coastline)
73	LOGICAL,SAVE,ALLOCATABLE :: is_streamflow_r(:) !! is stream flow point - (local routing grid)
74	!$OMP THREADPRIVATE(is_streamflow_r)
75	LOGICAL,SAVE,ALLOCATABLE,PUBLIC,PROTECTED :: routing_mask_r(:) !! valid routing point - (local routing grid)
76	!$OMP THREADPRIVATE(routing_mask_r)
77	LOGICAL,SAVE,ALLOCATABLE :: coast_mask(:) !! is a coast point - (local native grid)
78	!$OMP THREADPRIVATE(coast_mask)
79	INTEGER,SAVE :: total_coast_points !! global number of coast point - (local native grid)
80	!$OMP THREADPRIVATE(total_coast_points)
81	INTEGER,SAVE :: nbpt_r !! number of point in local routing grid
82	!$OMP THREADPRIVATE(nbpt_r)
83	INTEGER,SAVE :: nbpt_rp1 !! number of point in local routing grid with halo of 1
84	!$OMP THREADPRIVATE(nbpt_rp1)
85	REAL(r_std),SAVE,ALLOCATABLE :: routing_weight(:) !! Weight to transform runoff and drainage flux to water quantity in a conservative way (local native grid -> local routing grid)
86	!$OMP THREADPRIVATE(routing_weight)
87	REAL(r_std),SAVE,ALLOCATABLE :: routing_weight_in(:) !! Weight to transform runoff and drainage flux to water quantity in a conservative way (local native grid -> local routing grid)
88	!$OMP THREADPRIVATE(routing_weight_in)
89	REAL(r_std),SAVE,ALLOCATABLE :: unrouted_weight(:) !! Weight to transform runoff and drainage flux to water quantity in a conservative way (local native grid -> local routing grid)
90	!$OMP THREADPRIVATE(unrouted_weight)
91
92	REAL(r_std),SAVE,ALLOCATABLE :: weight_coast_to_coast_r(:)
93	!$OMP THREADPRIVATE(weight_coast_to_coast_r)
94
95	REAL(r_std),SAVE,ALLOCATABLE :: weight_coast_to_lake_r(:)
96	!$OMP THREADPRIVATE(weight_coast_to_lake_r)
97
98	REAL(r_std),SAVE,ALLOCATABLE :: weight_lake_to_coast_r(:)
99	!$OMP THREADPRIVATE(weight_lake_to_coast_r)
100
101	REAL(r_std),SAVE,ALLOCATABLE :: weight_lake_to_lake_r(:)
102	!$OMP THREADPRIVATE(weight_lake_to_lake_r)
103
104	REAL(r_std),SAVE,ALLOCATABLE :: basins_area_r(:)
105	!$OMP THREADPRIVATE(basins_area_r)
106
107	!! when doing interpolation from local routing grid to local native grid (river flow+coastal flow)
108	REAL(r_std),SAVE,ALLOCATABLE :: basins_extended_r(:) !! basins riverflow id (local routing grid)
109	!$OMP THREADPRIVATE(basins_extended_r)
110	INTEGER(i_std) :: basins_count !! number of basins (local routing grid)
111	!$OMP THREADPRIVATE(basins_count)
112	INTEGER(i_std) :: basins_out !! number of basins to output for diag
113	!$OMP THREADPRIVATE(basins_out)
114
115	INTEGER(i_std) :: split_routing !! time spliting for routing
116	!$OMP THREADPRIVATE(split_routing)
117
118
119	INTEGER :: nb_station
120	CHARACTER(LEN=60),ALLOCATABLE :: station(:)
121	REAL,ALLOCATABLE :: station_lonlat(:,:)
122	INTEGER,ALLOCATABLE :: station_index(:)
123	INTEGER :: station_ts = 0
124
125
126	! INTEGER(i_std), PARAMETER :: nb_stations=14
127	! REAL(r_std),PARAMETER :: station_lon(nb_stations) = &
128	! (/ -162.8830, -90.9058, -55.5110, -49.3242, -133.7447, -63.6000, 28.7167, &
129	! 15.3000, 66.5300, 89.6700, 86.5000, 127.6500, 3.3833, 117.6200 /)
130	! REAL(r_std),PARAMETER :: station_lat(nb_stations) = &
131	! (/ 61.9340, 32.3150, -1.9470, -5.1281, 67.4583, 8.1500, 45.2167, &
132	! -4.3000, 66.5700, 25.1800, 67.4800, 70.7000, 11.8667, 30.7700 /)
133	! CHARACTER(LEN=17),PARAMETER :: station_name(nb_stations) = &
134	! (/ "Pilot station ", "Vicksburg ", "Obidos ", &
135	! "Itupiranga ", "Arctic red river ", "Puente Angostura ", &
136	! "Ceatal Izmail ", "Kinshasa ", "Salekhard ", &
137	! "Bahadurabad ", "Igarka ", "Kusur ", &
138	! "Malanville ", "Datong " /)
139
140	CONTAINS
141
142
143	SUBROUTINE routing_flow_get(slow_reservoir_r, fast_reservoir_r, stream_reservoir_r, riverflow_mean, coastalflow_mean, &
144	lakeinflow_mean)
145	IMPLICIT NONE
146	REAL(r_std),OPTIONAL, INTENT(OUT) :: slow_reservoir_r(:)
147	REAL(r_std),OPTIONAL, INTENT(OUT) :: fast_reservoir_r(:)
148	REAL(r_std),OPTIONAL, INTENT(OUT) :: stream_reservoir_r(:)
149	REAL(r_std),OPTIONAL, INTENT(OUT) :: riverflow_mean(:)
150	REAL(r_std),OPTIONAL, INTENT(OUT) :: coastalflow_mean(:)
151	REAL(r_std),OPTIONAL, INTENT(OUT) :: lakeinflow_mean(:)
152
153	CALL routing_flow_get_(slow_reservoir_r, fast_reservoir_r, stream_reservoir_r, riverflow_mean, coastalflow_mean, lakeinflow_mean)
154	END SUBROUTINE routing_flow_get
155
156	SUBROUTINE routing_flow_get_(slow_reservoir_r_, fast_reservoir_r_, stream_reservoir_r_, riverflow_mean_, &
157	coastalflow_mean_, lakeinflow_mean_)
158	IMPLICIT NONE
159	REAL(r_std),OPTIONAL, INTENT(OUT) :: slow_reservoir_r_(:)
160	REAL(r_std),OPTIONAL, INTENT(OUT) :: fast_reservoir_r_(:)
161	REAL(r_std),OPTIONAL, INTENT(OUT) :: stream_reservoir_r_(:)
162	REAL(r_std),OPTIONAL, INTENT(OUT) :: riverflow_mean_(:)
163	REAL(r_std),OPTIONAL, INTENT(OUT) :: coastalflow_mean_(:)
164	REAL(r_std),OPTIONAL, INTENT(OUT) :: lakeinflow_mean_(:)
165
166	IF (PRESENT(slow_reservoir_r_)) slow_reservoir_r_ = slow_reservoir_r
167	IF (PRESENT(fast_reservoir_r_)) fast_reservoir_r_ = fast_reservoir_r
168	IF (PRESENT(stream_reservoir_r_)) stream_reservoir_r_ = stream_reservoir_r
169	IF (PRESENT(riverflow_mean_)) riverflow_mean_ = riverflow_mean
170	IF (PRESENT(coastalflow_mean_)) coastalflow_mean_ = coastalflow_mean
171	IF (PRESENT(lakeinflow_mean_)) lakeinflow_mean_ = lakeinflow_mean
172
173	END SUBROUTINE routing_flow_get_
174
175	SUBROUTINE routing_flow_set(slow_reservoir_r, fast_reservoir_r, stream_reservoir_r, riverflow_mean, coastalflow_mean, &
176	lakeinflow_mean)
177	IMPLICIT NONE
178	REAL(r_std),OPTIONAL, INTENT(IN) :: slow_reservoir_r(:)
179	REAL(r_std),OPTIONAL, INTENT(IN) :: fast_reservoir_r(:)
180	REAL(r_std),OPTIONAL, INTENT(IN) :: stream_reservoir_r(:)
181	REAL(r_std),OPTIONAL, INTENT(IN) :: riverflow_mean(:)
182	REAL(r_std),OPTIONAL, INTENT(IN) :: coastalflow_mean(:)
183	REAL(r_std),OPTIONAL, INTENT(IN) :: lakeinflow_mean(:)
184
185	CALL routing_flow_set_(slow_reservoir_r, fast_reservoir_r, stream_reservoir_r, riverflow_mean, coastalflow_mean, lakeinflow_mean)
186
187	END SUBROUTINE routing_flow_set
188
189	SUBROUTINE routing_flow_set_(slow_reservoir_r_, fast_reservoir_r_, stream_reservoir_r_, riverflow_mean_, &
190	coastalflow_mean_, lakeinflow_mean_)
191	IMPLICIT NONE
192	REAL(r_std),OPTIONAL, INTENT(IN) :: slow_reservoir_r_(:)
193	REAL(r_std),OPTIONAL, INTENT(IN) :: fast_reservoir_r_(:)
194	REAL(r_std),OPTIONAL, INTENT(IN) :: stream_reservoir_r_(:)
195	REAL(r_std),OPTIONAL, INTENT(IN) :: riverflow_mean_(:)
196	REAL(r_std),OPTIONAL, INTENT(IN) :: coastalflow_mean_(:)
197	REAL(r_std),OPTIONAL, INTENT(IN) :: lakeinflow_mean_(:)
198
199	IF (PRESENT(slow_reservoir_r_)) slow_reservoir_r = slow_reservoir_r_
200	IF (PRESENT(fast_reservoir_r_)) fast_reservoir_r = fast_reservoir_r_
201	IF (PRESENT(stream_reservoir_r_)) stream_reservoir_r = stream_reservoir_r_
202	IF (PRESENT(riverflow_mean_)) riverflow_mean = riverflow_mean_
203	IF (PRESENT(coastalflow_mean_)) coastalflow_mean = coastalflow_mean_
204	IF (PRESENT(lakeinflow_mean_)) lakeinflow_mean = lakeinflow_mean_
205
206	END SUBROUTINE routing_flow_set_
207	!! =============================================================================================================================
208	!! SUBROUTINE: routing_simple_xios_initialize
209	!!
210	!>\BRIEF Initialize xios dependant defintion before closing context defintion
211	!!
212	!! DESCRIPTION: Initialize xios dependant defintion before closing context defintion.
213	!! This subroutine is called before the xios context is closed.
214	!!
215	!! RECENT CHANGE(S): None
216	!!
217	!! REFERENCE(S): None
218	!!
219	!! FLOWCHART: None
220	!! \n
221	!_ ==============================================================================================================================
222
223	SUBROUTINE routing_flow_xios_initialize
224	USE xios
225	USE routing, ONLY : routing_names
226	IMPLICIT NONE
227
228	INTEGER(i_std) ::ib
229	INTEGER :: nbasmax=1
230	!! 0 Variable and parameter description
231	CHARACTER(LEN=60),ALLOCATABLE :: label(:)
232	LOGICAL :: file_exists
233
234	IF (is_omp_root) THEN
235	CALL xios_get_axis_attr("basins", n_glo=basins_out) ! get nb basins to output
236	ALLOCATE(label(basins_out))
237	CALL routing_names(basins_out,label)
238	CALL xios_set_axis_attr("basins", label=label) ! set riverflow basins name
239	INQUIRE(FILE="routing_start.nc", EXIST=file_exists)
240	IF (file_exists) CALL xios_set_file_attr("routing_start", enabled=.TRUE.)
241	ENDIF
242
243	!! Define XIOS axis size needed for the model output
244	! Add axis for homogeneity between all routing schemes, these dimensions are currently not used in this scheme
245	CALL xios_orchidee_addaxis("nbhtu", nbasmax, (/(REAL(ib,r_std),ib=1,nbasmax)/))
246	CALL xios_orchidee_addaxis("nbasmon", 1, (/(REAL(ib,r_std),ib=1,1)/))
247
248	END SUBROUTINE routing_flow_xios_initialize
249
250
251	SUBROUTINE routing_flow_initialize(kjit, rest_id, nbpt_, dt_routing, contfrac, nbpt_r_, riverflow, coastalflow)
252	IMPLICIT NONE
253	INTEGER ,INTENT(IN) :: kjit
254	INTEGER ,INTENT(IN) :: rest_id
255	INTEGER, INTENT(IN) :: nbpt_ !! nb points orchidee grid
256	REAL(r_std), INTENT(IN) :: dt_routing
257	REAL(r_std), INTENT(IN) :: contfrac(nbpt_) !! fraction of land
258	INTEGER, INTENT(OUT) :: nbpt_r_ !! nb points routing native grid
259	REAL(r_std), INTENT(OUT) :: riverflow(nbpt_) !! Outflow of the major rivers. The flux will be located on the continental grid but this should be a coastal point (kg/dt)
260	REAL(r_std), INTENT(OUT) :: coastalflow(nbpt_) !! Outflow on coastal points by small basins. This is the water which flows in a disperse way into the ocean (kg/dt)
261	INTEGER :: ier
262	CHARACTER(LEN=80) :: var_name !! To store variables names for I/O (unitless)
263
264	nbpt = nbpt_
265	CALL routing_flow_init_local(contfrac, nbpt_r_)
266	nbpt_r = nbpt_r_
267	CALL routing_flow_init_mean(kjit, rest_id)
268	CALL initialize_stations(dt_routing)
269	!
270	! Put into the restart file the fluxes so that they can be regenerated at restart.
271	!
272	ALLOCATE (lakeinflow_mean(nbpt), stat=ier)
273	IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for lakeinflow_mean','','')
274	var_name = 'lakeinflow'
275	CALL ioconf_setatt_p('UNITS', 'Kg/dt')
276	CALL ioconf_setatt_p('LONG_NAME','Lake inflow')
277	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., lakeinflow_mean, "gather", nbp_glo, index_g)
278	CALL setvar_p (lakeinflow_mean, val_exp, 'NO_KEYWORD', zero)
279
280	ALLOCATE (riverflow_mean(nbpt), stat=ier)
281	IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for riverflow_mean','','')
282	var_name = 'riverflow'
283	CALL ioconf_setatt_p('UNITS', 'Kg/dt')
284	CALL ioconf_setatt_p('LONG_NAME','River flux into the sea')
285	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., riverflow_mean, "gather", nbp_glo, index_g)
286	CALL setvar_p (riverflow_mean, val_exp, 'NO_KEYWORD', zero)
287
288	ALLOCATE (coastalflow_mean(nbpt), stat=ier)
289	IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for coastalflow_mean','','')
290	var_name = 'coastalflow'
291	CALL ioconf_setatt_p('UNITS', 'Kg/dt')
292	CALL ioconf_setatt_p('LONG_NAME','Diffuse flux into the sea')
293	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., coastalflow_mean, "gather", nbp_glo, index_g)
294	CALL setvar_p (coastalflow_mean, val_exp, 'NO_KEYWORD', zero)
295
296	riverflow(:) = riverflow_mean(:)
297	coastalflow(:) = coastalflow_mean(:)
298
299	END SUBROUTINE routing_flow_initialize
300
301
302	SUBROUTINE compute_coastline(contfrac, coastline)
303	USE mod_orchidee_para
304	USE grid
305	IMPLICIT NONE
306	REAL(r_std),INTENT(IN) :: contfrac(nbpt) ! INPUT : fraction of continent (unitless)
307	LOGICAL,INTENT(OUT) :: coastline(nbpt) ! OUTPUT : coastline mask (true : on coastaline, else false
308
309	REAL(r_std) :: contfrac_glo(nbp_glo)
310	REAL(r_std) :: contfrac2D_glo(iim_g*jjm_g)
311	LOGICAL :: coastline_glo(nbp_glo)
312	LOGICAL :: coastline2D_glo(iim_g*jjm_g)
313	INTEGER(i_std) :: i,j,ij,next_i,next_j,next_ij,m,k
314	REAL(r_std), PARAMETER :: epsilon=1e-5
315	LOGICAL :: is_periodic
316
317	is_periodic = global
318	CALL gather(contfrac,contfrac_glo)
319	contfrac2D_glo(:)=0
320	DO i=1,nbp_glo
321	contfrac2D_glo(index_g(i)) = contfrac_glo(i)
322	ENDDO
323
324	IF (is_mpi_root .AND. is_omp_root) THEN
325
326	coastline2D_glo(:)=.FALSE.
327	DO j=1,jjm_g
328	DO i=1,iim_g
329	ij=(j-1)*iim_g+i
330	IF (contfrac2D_glo(ij)>epsilon .AND. contfrac2D_glo(ij)<1-epsilon) THEN
331	coastline2D_glo(ij)=.TRUE.
332	ELSE IF (contfrac2D_glo(ij)>=1-epsilon .AND. grid_type/=unstructured ) THEN
333	DO k=-1,1
334	DO m=-1,1
335	IF (k==0 .AND. m==0) CYCLE
336	next_i=i+k
337	next_j=j+m
338
339	IF (next_i==0) THEN ! manage periodicity
340	IF (is_periodic) THEN
341	next_i=iim_g
342	ELSE
343	coastline2D_glo(ij)=.TRUE.
344	CYCLE
345	ENDIF
346	ENDIF
347
348	IF (next_i==iim_g+1) THEN! manage periodicity
349	IF (is_periodic) THEN
350	next_i=1 !! periodic but not true for limited area
351	ELSE
352	coastline2D_glo(ij)=.TRUE.
353	CYCLE
354	ENDIF
355	ENDIF
356	IF (next_j==0 .OR. next_j==jjm_g+1) THEN
357	IF (.NOT. is_periodic) coastline2D_glo(ij)=.TRUE.
358	CYCLE
359	ENDIF
360	next_ij = (next_j-1)*iim_g+next_i
361	IF (contfrac2D_glo(next_ij)<=epsilon) coastline2D_glo(ij)=.TRUE.
362	ENDDO
363	ENDDO
364	ENDIF
365	ENDDO
366	ENDDO
367
368	DO i=1,nbp_glo
369	coastline_glo(i) = coastline2D_glo(index_g(i))
370	ENDDO
371
372	ENDIF
373
374	CALL scatter(coastline_glo, coastline)
375
376	END SUBROUTINE compute_coastline
377
378
379
380	SUBROUTINE routing_flow_correct_riverflow(ni, nj, contfrac, coastline, trip_r, trip_extended_r, topoind_r)
381	USE xios
382	USE grid, ONLY : global
383	IMPLICIT NONE
384	INCLUDE "mpif.h"
385	INTEGER, INTENT(IN) :: ni ! INPUT : size i (longitude) of the local routing domain
386	INTEGER, INTENT(IN) :: nj ! INPUT : size i (latitude) of the local routing domain
387	REAL(r_std), INTENT(IN) :: contfrac(nbp_mpi) ! INPUT : continental fraction (unittless)
388	LOGICAL, INTENT(IN) :: coastline(nbp_mpi) ! INPUT/OUTPUT : coastline mask
389	REAL(r_std), INTENT(INOUT) :: trip_r(ni,nj) ! INPUT/OUTPUT : diection of flow, which will be modified by the routine
390	REAL(r_std), INTENT(INOUT) :: trip_extended_r(ni,nj) ! INPUT : direction of flow extended to ocean points
391	REAL(r_std), INTENT(INOUT) :: topoind_r(ni,nj) ! INPUT/OUTPUT : topographic index which will be modified by the routine if extended to ocean
392
393	REAL(r_std) :: contfrac_r(ni,nj)
394	REAL(r_std) :: mask_coast(nbp_mpi)
395	REAL(r_std) :: frac_coast_r(ni,nj)
396	REAL(r_std) :: trip_extended_rp1(0:ni+1,0:nj+1)
397	REAL(r_std) :: trip_rp1(0:ni+1,0:nj+1)
398	REAL(r_std) :: state_r(ni,nj)
399	REAL(r_std) :: state_rp1(0:ni+1,0:nj+1)
400
401	INTEGER, PARAMETER :: is_ter=0
402	INTEGER, PARAMETER :: is_coast=1
403	INTEGER, PARAMETER :: is_oce=2
404	INTEGER :: updated
405	INTEGER :: it, i,j,k,m, nexti,nextj,previ,prevj,ierr
406	INTEGER :: ibegin, jbegin, ni_glo, nj_glo
407	REAL(r_std) :: lonvalue_1d(ni), latvalue_1d(nj)
408	REAL(r_std),PARAMETER :: epsilon=1e-3
409
410	TYPE(xios_duration) :: ts
411	TYPE(xios_domain) :: domain_hdl
412	TYPE(xios_domaingroup) :: domain_def_hdl
413	TYPE(xios_field) :: field_hdl
414	TYPE(xios_fieldgroup) :: field_def_hdl
415	TYPE(xios_file) :: file_hdl
416	TYPE(xios_filegroup) :: file_def_hdl
417	TYPE(xios_expand_domain) :: domain_expand_hdl
418	TYPE(xios_date) :: start_date, time_origin
419	CHARACTER(LEN=20) :: calendar_type
420	LOGICAL :: true=.TRUE.
421
422
423	CALL xios_get_domain_attr("routing_domain", ibegin=ibegin, jbegin=jbegin, ni_glo=ni_glo,nj_glo=nj_glo) ! get routing domain dimension
424	CALL xios_get_domain_attr("routing_domain", lonvalue_1d=lonvalue_1d, latvalue_1d=latvalue_1d) ! get routing domain dimension
425
426	! manage non periodic boundary in case of LAM
427	IF (.NOT. global) THEN
428	IF (ibegin==0) THEN
429	i=1 ;
430	DO j=1,nj
431	IF (trip_r(i,j)>=1 .AND. trip_r(i,j)<=8) trip_r(i,j) = 97
432	IF (trip_extended_r(i,j)>=1 .AND. trip_extended_r(i,j)<=8) trip_extended_r(i,j) = 97
433	ENDDO
434	ENDIF
435
436	IF (ibegin+ni==ni_glo) THEN
437	i=ni ;
438	DO j=1,nj
439	IF (trip_r(i,j)>=1 .AND. trip_r(i,j)<=8) trip_r(i,j) = 97
440	IF (trip_extended_r(i,j)>=1 .AND. trip_extended_r(i,j)<=8) trip_extended_r(i,j) = 97
441	ENDDO
442	ENDIF
443
444	IF (jbegin==0) THEN
445	j=1 ;
446	DO i=1,ni
447	IF (trip_r(i,j)>=1 .AND. trip_r(i,j)<=8) trip_r(i,j) = 97
448	IF (trip_extended_r(i,j)>=1 .AND. trip_extended_r(i,j)<=8) trip_extended_r(i,j) = 97
449	ENDDO
450	ENDIF
451
452	IF (jbegin+nj==nj_glo) THEN
453	j=nj ;
454	DO i=1,ni
455	IF (trip_r(i,j)>=1 .AND. trip_r(i,j)<=8) trip_r(i,j) = 97
456	IF (trip_extended_r(i,j)>=1 .AND. trip_extended_r(i,j)<=8) trip_extended_r(i,j) = 97
457	ENDDO
458	ENDIF
459	ENDIF
460
461
462	CALL xios_send_field("routing_contfrac",contfrac)
463	CALL xios_recv_field("routing_contfrac_r", contfrac_r)
464
465	CALL xios_send_field("trip_ext_r",trip_extended_r)
466	CALL xios_recv_field("trip_ext_rp1",trip_extended_rp1)
467
468	mask_coast=0
469	WHERE(coastline) mask_coast=1
470	CALL xios_send_field("mask_coastline",mask_coast)
471	CALL xios_recv_field("frac_coastline_r",frac_coast_r)
472
473	DO j=1,nj
474	DO i=1,ni
475	IF (contfrac_r(i,j)<epsilon .OR. contfrac_r(i,j)> 2 ) THEN
476	state_r(i,j)=is_oce
477	ELSE IF (contfrac_r(i,j)>1-epsilon) THEN
478	state_r(i,j)=is_ter
479	ELSE
480	state_r(i,j)=is_coast
481	ENDIF
482
483	IF (frac_coast_r(i,j)>epsilon .AND. frac_coast_r(i,j)<=1+epsilon ) state_r(i,j)=is_coast
484	ENDDO
485	ENDDO
486
487
488	DO j=1,nj
489	DO i=1,ni
490	IF (trip_r(i,j)>=100) THEN
491	IF (state_r(i,j) /= is_oce) trip_r(i,j)=trip_extended_r(i,j)
492	ENDIF
493	ENDDO
494	ENDDO
495
496	! create a new XIOS context to solve iterative process of correction
497	! so we can use filter to update halo
498	! CALL xios_get_calendar_type(calendar_type)
499	! CALL xios_get_start_date(start_date)
500	! CALL xios_get_time_origin(time_origin)
501	!
502	! CALL xios_context_initialize("orchidee_routing_correction",MPI_COMM_ORCH)
503	! CALL xios_orchidee_change_context("orchidee_routing_correction")
504	! ts.second=1
505	! calendar_type="gregorian"
506	! CALL xios_define_calendar(type=calendar_type,start_date=start_date,time_origin=time_origin, timestep=ts )
507
508	! CALL xios_get_handle("domain_definition",domain_def_hdl)
509	! CALL xios_add_child(domain_def_hdl, domain_hdl, "routing_domain")
510	! CALL xios_set_attr(domain_hdl, type="rectilinear", ibegin=ibegin, jbegin=jbegin, ni=ni, nj=nj, ni_glo=ni_glo, nj_glo=nj_glo) ! set routing domain dimension
511	! CALL xios_set_attr(domain_hdl, lonvalue_1d=lonvalue_1d, latvalue_1d=latvalue_1d) ! set routing domain dimension
512	! CALL xios_add_child(domain_def_hdl, domain_hdl, "routing_domain_expand")
513	! CALL xios_set_attr(domain_hdl, domain_ref="routing_domain")
514	! CALL xios_add_child(domain_hdl,domain_expand_hdl)
515	! CALL xios_set_attr(domain_expand_hdl, type="edge", i_periodic=.TRUE., j_periodic=.TRUE.)
516	!
517	! CALL xios_get_handle("field_definition",field_def_hdl)
518	! CALL xios_add_child(field_def_hdl,field_hdl,"trip_r_init")
519	! CALL xios_set_attr(field_hdl,domain_ref="routing_domain")
520	! CALL xios_add_child(field_def_hdl,field_hdl,"trip_r_final")
521	! CALL xios_set_attr(field_hdl,domain_ref="routing_domain")
522
523	! CALL xios_add_child(field_def_hdl,field_hdl,"trip_r")
524	! CALL xios_set_attr(field_hdl,domain_ref="routing_domain")
525	! CALL xios_add_child(field_def_hdl,field_hdl,"trip_rp1")
526	! CALL xios_set_attr(field_hdl, field_ref="trip_r", domain_ref="routing_domain_expand", read_access=.TRUE.)
527
528	! CALL xios_add_child(field_def_hdl,field_hdl,"state_r")
529	! CALL xios_set_attr(field_hdl,domain_ref="routing_domain")
530	! CALL xios_add_child(field_def_hdl,field_hdl,"state_rp1")
531	! CALL xios_set_attr(field_hdl, field_ref="state_r", domain_ref="routing_domain_expand", read_access=.TRUE.)
532	!
533	! CALL xios_get_handle("file_definition",file_def_hdl)
534	! CALL xios_add_child(file_def_hdl,file_hdl,"routing_correction")
535	! CALL xios_set_attr(file_hdl,type="one_file", output_freq=ts, sync_freq=ts, enabled=.TRUE.)
536	! CALL xios_add_child(file_hdl, field_hdl)
537	! CALL xios_set_attr(field_hdl, field_ref="trip_r", operation="instant")
538	! CALL xios_add_child(file_hdl, field_hdl)
539	! CALL xios_set_attr(field_hdl, field_ref="state_r", operation="instant")
540	! CALL xios_add_child(file_hdl, field_hdl)
541	! CALL xios_set_attr(field_hdl, field_ref="trip_r_init", operation="once")
542	! CALL xios_add_child(file_hdl, field_hdl)
543	! CALL xios_set_attr(field_hdl, field_ref="state_r_init", operation="once")
544
545
546	CALL xios_context_initialize("orchidee_init_routing",MPI_COMM_ORCH)
547	CALL xios_orchidee_change_context("orchidee_init_routing")
548	CALL xios_set_domain_attr("routing_domain", ibegin=ibegin, jbegin=jbegin, ni=ni, nj=nj, ni_glo=ni_glo, nj_glo=nj_glo) ! get routing domain dimension
549	CALL xios_set_domain_attr("routing_domain", lonvalue_1d=lonvalue_1d, latvalue_1d=latvalue_1d) ! get routing domain dimension
550	CALL xios_close_context_definition()
551	CALL xios_update_calendar(1)
552
553	CALL xios_send_field("trip_r_init",trip_r)
554
555	CALL xios_send_field("trip_r",trip_r)
556	CALL xios_recv_field("trip_rp1",trip_rp1)
557
558	CALL xios_send_field("state_r",state_r)
559	CALL xios_recv_field("state_rp1",state_rp1)
560
561	DO j=1,nj
562	DO i=1,ni
563	IF (trip_rp1(i,j)<=8) THEN
564	CALL next_trip(trip_rp1, ni, nj, i, j, nexti, nextj)
565	IF (trip_rp1(nexti,nextj)>=100) trip_rp1(i,j)=98 ! unrouted point becomes coastal flow
566	ENDIF
567	ENDDO
568	ENDDO
569
570
571	state_r(1:ni,1:nj) = state_rp1(1:ni,1:nj)
572	trip_r(1:ni,1:nj) = trip_rp1(1:ni,1:nj)
573
574	updated=1
575	it=2
576	DO WHILE(updated>0 .AND. it<MAX(ni_glo,nj_glo))
577	updated=0
578
579	CALL xios_update_calendar(it)
580
581	CALL xios_send_field("trip_r", trip_r)
582	CALL xios_recv_field("trip_rp1",trip_rp1)
583
584	CALL xios_send_field("state_r",state_r)
585	CALL xios_recv_field("state_rp1",state_rp1)
586
587	state_r(1:ni,1:nj) = state_rp1(1:ni,1:nj)
588	trip_r(1:ni,1:nj) = trip_rp1(1:ni,1:nj)
589
590
591
592	DO j=0,nj+1
593	DO i=0,ni+1
594
595	IF (trip_rp1(i,j)==97 .OR. trip_rp1(i,j)==98 .OR. trip_rp1(i,j)==99) THEN
596
597	IF (state_rp1(i,j)==is_coast) THEN ! ok
598
599	IF (trip_rp1(i,j)==97) THEN
600	IF (i>=1 .AND. i<=ni .AND. j>=1 .AND. j<=nj) trip_r(i,j)=98 ! => lakinflow will becomes coastalflow
601	ELSE
602	! ok for riverflow
603	ENDIF
604
605	ELSE IF (state_rp1(i,j) == is_ter) THEN ! too short
606
607	IF (trip_rp1(i,j)==97) THEN
608	!ok for lakinflow
609	ELSE ! => for riverflow/coastalflow
610	IF (trip_extended_rp1(i,j)==97) THEN ! definitively a lake that cannot be routed further away
611	IF (i>=1 .AND. i<=ni .AND. j>=1 .AND. j<=nj) THEN
612	trip_r(i,j)=97
613	updated=updated+1
614	ENDIF
615	ELSE
616	CALL next_trip(trip_extended_rp1, ni, nj, i, j, nexti, nextj)
617	IF (nexti/=-1 .AND. nextj/=-1) THEN
618	IF (nexti==i .AND. nextj==j) THEN
619	!! routed on itself in midle of land => lake
620	IF (nexti>=1 .AND. nexti<=ni .AND. nextj>=1 .AND. nextj<=nj) trip_r(i,j)=97
621	ELSE
622
623	IF (state_rp1(nexti,nextj)==is_coast .OR. state_rp1(nexti,nextj)==is_ter) THEN
624
625	IF (i>=1 .AND. i<=ni .AND. j>=1 .AND. j<=nj) THEN
626	trip_r(i,j)=trip_extended_rp1(i,j)
627	state_r(i,j)= is_ter
628	topoind_r(i,j) = 1e-10 ! flow instantanesly
629	updated=updated+1
630	ENDIF
631
632	IF (nexti>=1 .AND. nexti<=ni .AND. nextj>=1 .AND. nextj<=nj) THEN
633	IF (trip_r(nexti,nextj)/=trip_rp1(i,j)) THEN
634	trip_r(nexti,nextj)=trip_rp1(i,j)
635	updated=updated+1
636	ENDIF
637	ENDIF
638
639	ELSE
640	IF (i>=1 .AND. i<=ni .AND. j>=1 .AND. j<=nj) THEN
641	STOP 'not possible'
642	ENDIF
643	ENDIF
644	ENDIF
645	ENDIF
646	ENDIF
647	ENDIF
648
649	ELSE IF (state_rp1(i,j) == is_oce) THEN ! too long
650
651	DO k=-1,1
652	DO m=-1,1
653	IF (k==0 .AND. m==0) CYCLE
654	prevj = j+k ; previ = i+m
655	CALL next_trip(trip_rp1, ni, nj, previ, prevj, nexti, nextj)
656	! IF (nexti==i .AND. nextj==j .AND. previ>=1 .AND. previ<=ni .AND. prevj>=1 .AND. prevj<=nj) trip_r(previ,prevj)=trip_rp1(i,j)
657	IF (nexti==i .AND. nextj==j .AND. previ>=1 .AND. previ<=ni .AND. prevj>=1 .AND. prevj<=nj) trip_r(previ,prevj)=trip_rp1(i,j)
658	ENDDO
659	ENDDO
660
661	IF (i>=1 .AND. i<=ni .AND. j>=1 .AND. j<=nj) THEN
662	trip_r(i,j)=1000
663	updated=updated+1
664	ENDIF
665	ENDIF
666
667	ENDIF
668	ENDDO
669	ENDDO
670
671	CALL MPI_ALLREDUCE(MPI_IN_PLACE , updated, 1, MPI_INT_ORCH, MPI_SUM, MPI_COMM_ORCH, ierr)
672	it = it +1
673	ENDDO
674
675	CALL xios_send_field("trip_r_final",trip_r)
676
677	CALL xios_context_finalize
678	CALL xios_orchidee_change_context("orchidee")
679
680	CONTAINS
681
682	SUBROUTINE next_trip(trip, ni, nj, i, j, nexti, nextj)
683	IMPLICIT NONE
684	REAL, INTENT(IN) :: trip(0:ni+1,0:nj+1)
685	INTEGER, INTENT(IN) :: ni
686	INTEGER, INTENT(IN) :: nj
687	INTEGER, INTENT(IN) :: i
688	INTEGER, INTENT(IN) :: j
689	INTEGER, INTENT(OUT) :: nexti
690	INTEGER, INTENT(OUT) :: nextj
691
692	IF (i<0 .OR. i>ni+1 .OR. j<0 .OR. j>nj+1) THEN
693	nexti=-1
694	nextj=-1
695	RETURN
696	ENDIF
697
698	SELECT CASE(NINT(trip(i,j)))
699	CASE (1)
700	nextj=j-1 ; nexti=i ! north
701	CASE (2)
702	nextj=j-1 ; nexti=i+1 ! north-east
703	CASE (3)
704	nextj=j ; nexti=i+1 ! east
705	CASE (4)
706	nextj=j+1 ; nexti=i+1 ! south-east
707	CASE (5)
708	nextj=j+1 ; nexti=i ! south
709	CASE(6)
710	nextj=j+1 ; nexti=i-1 ! south-west
711	CASE (7)
712	nextj=j ; nexti=i-1 ! west
713	CASE (8)
714	nextj=j-1 ; nexti=i-1 ! north-west
715	CASE DEFAULT
716	nextj=-1 ; nexti=-1 ! undefined
717	END SELECT
718
719	IF (nexti<0 .OR. nexti>ni+1) nexti=-1
720	IF (nextj<0 .OR. nextj>nj+1) nextj=-1
721
722	END SUBROUTINE next_trip
723
724	END SUBROUTINE routing_flow_correct_riverflow
725
726
727
728
729
730
731
732	SUBROUTINE routing_flow_init_mean(kjit, rest_id)
733	USE ioipsl
734	IMPLICIT NONE
735	INTEGER(i_std), INTENT(in) :: kjit
736	INTEGER(i_std), INTENT(in) :: rest_id !! Restart file identifier (unitless)
737
738	INTEGER :: ier
739	CHARACTER(LEN=80) :: var_name !! To store variables names for I/O (unitless)
740
741
742	! Get from the restart the fluxes we accumulated.
743	!
744
745	ALLOCATE (runoff_mean(nbpt), stat=ier)
746	runoff_mean(:) = 0
747	IF (ier /= 0) CALL ipslerr_p(3,'routing_flow_init_mean','Pb in allocate for runoff_mean','','')
748	var_name = 'runoff_route'
749	CALL ioconf_setatt_p('UNITS', 'Kg')
750	CALL ioconf_setatt_p('LONG_NAME','Accumulated runoff for routing')
751	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., runoff_mean, "gather", nbp_glo, index_g)
752	CALL setvar_p (runoff_mean, val_exp, 'NO_KEYWORD', zero)
753
754	ALLOCATE(drainage_mean(nbpt), stat=ier)
755	drainage_mean(:) = 0
756	IF (ier /= 0) CALL ipslerr_p(3,'routing_flow_init_mean','Pb in allocate for drainage_mean','','')
757	var_name = 'drainage_route'
758	CALL ioconf_setatt_p('UNITS', 'Kg')
759	CALL ioconf_setatt_p('LONG_NAME','Accumulated drainage for routing')
760	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., drainage_mean, "gather", nbp_glo, index_g)
761	CALL setvar_p (drainage_mean, val_exp, 'NO_KEYWORD', zero)
762
763	END SUBROUTINE routing_flow_init_mean
764
765	SUBROUTINE routing_flow_make_mean(runoff, drainage)
766	IMPLICIT NONE
767	REAL(r_std),INTENT(IN) :: runoff(:)
768	REAL(r_std),INTENT(IN) :: drainage(:)
769
770	runoff_mean(:) = runoff_mean(:) + runoff
771	drainage_mean(:) = drainage_mean(:) + drainage
772
773	END SUBROUTINE routing_flow_make_mean
774
775
776	SUBROUTINE routing_flow_reset_mean
777	IMPLICIT NONE
778	runoff_mean(:) = 0
779	drainage_mean(:) = 0
780
781	END SUBROUTINE routing_flow_reset_mean
782
783
784	SUBROUTINE routing_flow_finalize_mean(kjit, rest_id)
785	USE ioipsl
786	IMPLICIT NONE
787	INTEGER,INTENT(IN) :: kjit
788	INTEGER,INTENT(IN) :: rest_id
789
790	CALL restput_p (rest_id, 'runoff_route', nbp_glo, 1, 1, kjit, runoff_mean, 'scatter', nbp_glo, index_g)
791	CALL restput_p (rest_id, 'drainage_route', nbp_glo, 1, 1, kjit, drainage_mean, 'scatter', nbp_glo, index_g)
792
793	DEALLOCATE(runoff_mean)
794	DEALLOCATE(drainage_mean)
795
796	END SUBROUTINE routing_flow_finalize_mean
797
798
799	SUBROUTINE routing_flow_init_local(contfrac, nbpt_r)
800	USE xios
801	USE grid, ONLY : area
802	IMPLICIT NONE
803	INCLUDE "mpif.h"
804
805	!! 0 Variable and parameter description
806	!! 0.1 Input variables
807	REAL(r_std),INTENT(IN) :: contfrac(nbpt) !! fraction of land
808	INTEGER,INTENT(OUT) :: nbpt_r !! nb points routing native grid
809
810	!! 0.2 Local variables
811	INTEGER :: ni !! longitude dimension of local routing grid
812	INTEGER :: nj !! latitude dimension of local routing grid
813	REAL(r_std) :: contfrac_mpi(nbp_mpi)
814	REAL(r_std),ALLOCATABLE :: trip_rp1(:,:) !! direction of flow (1-8) or river flow (99) or coastal flow (98) or lake inflow (97) - local routing grid + halo (0:ni+1,0:nj+1)
815	REAL(r_std),ALLOCATABLE :: trip_extended_r(:,:) !! direction of flow (1-8) or river flow (99) or coastal flow (98) or lake inflow (97) - local routing grid (ni,nj)
816	!! routing is artificially computed on sea and endoric basins
817	REAL(r_std),ALLOCATABLE :: diag_r(:) !! fraction of routing cell intesected by coastal cells of native grid
818	REAL(r_std),ALLOCATABLE :: frac_lake_r(:) !! fraction of routing cell intesected by cells of native grid
819	REAL(r_std),ALLOCATABLE :: frac_coast_r(:)
820
821	REAL(r_std) :: diag(nbp_mpi) !! fraction of routing cell intesected by coastal cells of native grid
822
823	LOGICAL :: coastline(nbpt)
824	LOGICAL :: coastline_mpi(nbp_mpi)
825	REAL(r_std) :: area_mpi(nbp_mpi)
826
827	INTEGER :: ij, ij_r, ij_rp1, i,j,jp1,jm1,ip1,im1,jr,ir
828	INTEGER :: basins_count_mpi
829	INTEGER :: nb_coast_points
830	INTEGER :: ierr
831	LOGICAL :: file_exists
832	REAL(r_std) :: epsilon = 1e-5
833
834	!_ ================================================================================================================================
835	!
836	!> A value for property of each reservoir (in day/m) is given to compute a time constant (in day)
837	!> for each reservoir (product of tcst and topo_resid).
838	!> The value of tcst has been calibrated for the three reservoirs over the Senegal river basin only,
839	!> during the 1 degree NCEP Corrected by Cru (NCC) resolution simulations (Ngo-Duc et al., 2005, Ngo-Duc et al., 2006) and
840	!> generalized for all the basins of the world. The "slow reservoir" and the "fast reservoir"
841	!> have the highest value in order to simulate the groundwater.
842	!> The "stream reservoir", which represents all the water of the stream, has the lowest value.
843	!> Those figures are the same for all the basins of the world.
844
845	CALL getin_p('SLOW_TCST', slow_tcst)
846	!
847	!Config Key = FAST_TCST
848	!Config Desc = Time constant for the fast reservoir
849	!Config If = RIVER_ROUTING
850	!Config Def = 3.0
851	!Config Help = This parameters allows the user to fix the
852	!Config time constant (in days) of the fast reservoir
853	!Config in order to get better river flows for
854	!Config particular regions.
855	!Config Units = [days]
856
857	CALL getin_p('FAST_TCST', fast_tcst)
858	!
859	!Config Key = STREAM_TCST
860	!Config Desc = Time constant for the stream reservoir
861	!Config If = RIVER_ROUTING
862	!Config Def = 0.24
863	!Config Help = This parameters allows the user to fix the
864	!Config time constant (in days) of the stream reservoir
865	!Config in order to get better river flows for
866	!Config particular regions.
867	!Config Units = [days]
868
869	CALL getin_p('STREAM_TCST', stream_tcst)
870	!
871	!Config Key = FLOOD_TCST
872	!Config Desc = Time constant for the flood reservoir
873	!Config If = RIVER_ROUTING
874	!Config Def = 4.0
875	!Config Help = This parameters allows the user to fix the
876	!Config time constant (in days) of the flood reservoir
877	!Config in order to get better river flows for
878	!Config particular regions.
879	!Config Units = [days]
880
881	split_routing=1
882	CALL getin_p("SPLIT_ROUTING",split_routing)
883
884	CALL compute_coastline(contfrac, coastline)
885
886	CALL gather_omp(contfrac, contfrac_mpi)
887	CALL gather_omp(area, area_mpi)
888	CALL gather_omp(coastline, coastline_mpi)
889
890	ALLOCATE(is_coastline(nbp_mpi))
891	is_coastline=coastline_mpi
892
893	IF (is_omp_root) THEN
894
895	WHERE(contfrac_mpi<=epsilon) contfrac_mpi=0
896	WHERE(contfrac_mpi>=epsilon) contfrac_mpi=1
897
898	diag=0
899	WHERE(coastline_mpi) diag=1
900	CALL xios_send_field("is_coastline", diag)
901
902	CALL xios_get_domain_attr("routing_domain", ni=ni, nj=nj) ! get routing domain dimension
903
904	nbpt_r= ni*nj
905	nbpt_rp1= (ni+2)*(nj+2)
906
907	! Allocate module variable
908	ALLOCATE(fast_reservoir_r(ni*nj))
909	ALLOCATE(slow_reservoir_r(ni*nj))
910	ALLOCATE(stream_reservoir_r(ni*nj))
911
912	ALLOCATE(is_lakeinflow_r(ni*nj))
913	ALLOCATE(is_coastalflow_r(ni*nj))
914	ALLOCATE(is_riverflow_r(ni*nj))
915	ALLOCATE(is_streamflow_r(ni*nj))
916
917	ALLOCATE(topoind_r(ni*nj))
918	ALLOCATE(routing_mask_r(ni*nj))
919	ALLOCATE(basins_extended_r(nbpt_r))
920	ALLOCATE(route_flow_rp1((ni+2)*(nj+2)))
921
922	ALLOCATE(routing_weight(nbp_mpi))
923	ALLOCATE(routing_weight_in(nbp_mpi))
924	ALLOCATE(unrouted_weight(nbp_mpi))
925
926	ALLOCATE(diag_r(nbpt_r)) ! for diags
927	ALLOCATE(frac_coast_r(nbpt_r))
928	ALLOCATE(frac_lake_r(nbpt_r))
929	ALLOCATE(weight_coast_to_coast_r(nbpt_r))
930	ALLOCATE(weight_coast_to_lake_r(nbpt_r))
931	ALLOCATE(weight_lake_to_coast_r(nbpt_r))
932	ALLOCATE(weight_lake_to_lake_r(nbpt_r))
933
934	ALLOCATE(trip_extended_r(ni,nj))
935
936	is_lakeinflow_r(:) = .FALSE.
937	is_coastalflow_r(:) = .FALSE.
938	is_riverflow_r(:) = .FALSE.
939
940
941	INQUIRE(FILE="routing_start.nc", EXIST=file_exists)
942
943	IF (file_exists) THEN
944	CALL xios_recv_field("fast_reservoir_start",fast_reservoir_r)
945	CALL xios_recv_field("slow_reservoir_start",slow_reservoir_r)
946	CALL xios_recv_field("stream_reservoir_start",stream_reservoir_r)
947	ELSE
948	fast_reservoir_r(:)=0
949	slow_reservoir_r(:)=0
950	stream_reservoir_r(:)=0
951	ENDIF
952
953
954	ALLOCATE(trip_rp1(0:ni+1,0:nj+1))
955
956	trip_rp1(:,:)=1e10
957	CALL xios_recv_field("trip_r",trip_rp1(1:ni,1:nj)) ! recv trip array with halo of 1
958	CALL xios_recv_field("trip_extended_r",trip_extended_r) ! recv extended trip array from file
959	CALL xios_recv_field("topoind_r",topoind_r) ! recv topo index array from file
960	CALL xios_recv_field("basins_extended_r",basins_extended_r) ! recv basins index from file
961
962	CALL routing_flow_correct_riverflow(ni, nj, contfrac_mpi, coastline, trip_rp1(1:ni,1:nj), trip_extended_r, topoind_r)
963
964	CALL xios_send_field("trip_update_r",trip_rp1(1:ni,1:nj)) ! send to xios trip array to update for halo
965	CALL xios_recv_field("trip_rp1",trip_rp1) ! recv trip array with halo of 1
966
967
968
969	!! Compute the routing
970	!! Loop on all point point of the local routing grid + halo
971
972	route_flow_rp1(:)=-1
973
974	DO j=0,nj+1
975	jp1=j+1
976	jm1=j-1
977	DO i=0,ni+1
978	ij_rp1=i+(ni+2)*j+1
979	ij_r=i+ni*(j-1)
980	ip1=i+1
981	im1=i-1
982
983	IF (trip_rp1(i,j) < 100) THEN
984
985	IF (i>=1 .AND. i<=ni .AND. j>=1 .AND. j<=nj) routing_mask_r(ij_r)=.TRUE.
986
987	ir=-1 ; jr=-1 ! -> -1 for 97,98,99
988	SELECT CASE (NINT(trip_rp1(i,j))) ! get the trip value for each points
989	CASE (1)
990	jr=jm1 ; ir=i ! north
991	CASE (2)
992	jr=jm1 ; ir=ip1 ! north-east
993	CASE (3)
994	jr=j ; ir=ip1 ! east
995	CASE (4)
996	jr=jp1 ; ir=ip1 ! south-east
997	CASE (5)
998	jr=jp1 ; ir=i ! south
999	CASE(6)
1000	jr=jp1 ; ir=im1 ! south-west
1001	CASE (7)
1002	jr=j ; ir=im1 ! west
1003	CASE (8)
1004	jr=jm1 ; ir=im1 ! north-west
1005	CASE (97)
1006	IF ( i>0 .AND. i<ni+1 .AND. j>0 .AND. j<nj+1) THEN ! if inside my local domain
1007	is_lakeinflow_r(ij_r)=.TRUE. ! I am a lakeinflow point and route to myself
1008	jr=j ; ir=i
1009	ENDIF
1010	CASE (98)
1011	IF ( i>0 .AND. i<ni+1 .AND. j>0 .AND. j<nj+1) THEN ! if inside my local domain
1012	is_coastalflow_r(ij_r)=.TRUE. ! I am a coastal flow point and route to myself
1013	jr=j ; ir=i
1014	ENDIF
1015	CASE (99)
1016	IF ( i>0 .AND. i<ni+1 .AND. j>0 .AND. j<nj+1) THEN ! if inside my local domain
1017	is_riverflow_r(ij_r)=.TRUE. ! I am a riverflow point and route to myself
1018	jr=j ; ir=i
1019	ENDIF
1020	END SELECT
1021
1022	IF (ir<0 .OR. ir>ni+1 .OR. jr<0 .OR. jr>nj+1) THEN
1023	route_flow_rp1(ij_rp1)=-1 ! if route outside my local domain+halo, no routing (will be done by other process)
1024	ELSE
1025	IF (ir<1 .OR. ir>ni .OR. jr<1 .OR. jr>nj) THEN
1026	route_flow_rp1(ij_rp1)=-1 ! if route outside my local domain, no routing (will be done by other process)
1027	ELSE
1028	route_flow_rp1(ij_rp1)=ir+ni*(jr-1) ! define the cell where to flow
1029	IF (trip_rp1(ir,jr)>99) STOP 'Pb point not routed to outflow'
1030	ENDIF
1031	ENDIF
1032	ELSE
1033	IF (i>=1 .AND. i<=ni .AND. j>=1 .AND. j<=nj) routing_mask_r(ij_r)=.FALSE.
1034	ENDIF
1035	ENDDO
1036	ENDDO
1037	is_streamflow_r(:) = .NOT. (is_lakeinflow_r(:) .OR. is_coastalflow_r(:) .OR. is_riverflow_r(:)) .AND. routing_mask_r(:)
1038
1039	diag_r(:)=0
1040	WHERE(routing_mask_r) diag_r=1
1041	CALL xios_send_field("routing_mask_r", diag_r)
1042	CALL xios_recv_field("frac_routing", diag)
1043	routing_weight_in (:) = 0
1044	unrouted_weight (:) = 0
1045	WHERE (diag > 0) routing_weight_in = contfrac_mpi*area_mpi/diag
1046	WHERE (diag == 0 ) unrouted_weight = contfrac_mpi*area_mpi
1047
1048
1049	! compute the number of coast cells to distribute lakeinflow onto
1050	diag=0
1051	WHERE (coastline_mpi) diag=1
1052	nb_coast_points=SUM(diag)
1053	CALL reduce_sum_mpi(nb_coast_points, total_coast_points)
1054	CALL bcast_mpi(total_coast_points)
1055
1056
1057	DO ij=1,nbp_mpi
1058	routing_weight(ij)=contfrac_mpi(ij)*area_mpi(ij)
1059	ENDDO
1060
1061	! looking for basins
1062	basins_count_mpi=0
1063	DO ij=1,nbpt_r
1064	IF (basins_extended_r(ij) > basins_count_mpi) basins_count_mpi=basins_extended_r(ij)
1065	ENDDO
1066	CALL MPI_ALLREDUCE(basins_count_mpi,basins_count,1,MPI_INT_ORCH, MPI_MAX,MPI_COMM_ORCH,ierr)
1067
1068	diag(:) = 0
1069	WHERE (coastline_mpi) diag=1
1070	CALL xios_send_field("mask_coastal",diag)
1071	CALL xios_recv_field("frac_coastal_r",diag_r)
1072	WHERE (diag_r > 100.) ! missing_value
1073	diag_r=0
1074	ELSE WHERE (diag_r < epsilon)
1075	diag_r=0
1076	ELSEWHERE (diag_r > 1-epsilon)
1077	diag_r=1
1078	END WHERE
1079	frac_coast_r=diag_r ;
1080
1081	diag(:) = 0
1082	WHERE (.NOT. coastline_mpi) diag=1
1083	CALL xios_send_field("mask_lake",diag)
1084	CALL xios_recv_field("frac_lake_r",diag_r)
1085	WHERE (diag_r > 100.) ! missing_value
1086	diag_r=0.
1087	ELSEWHERE (diag_r < epsilon)
1088	diag_r=0.
1089	ELSEWHERE (diag_r > 1-epsilon)
1090	diag_r=1.
1091	END WHERE
1092	frac_lake_r = diag_r
1093
1094	weight_coast_to_coast_r(:)=0
1095	weight_coast_to_lake_r(:)=0
1096	weight_lake_to_coast_r(:)=0
1097	weight_lake_to_lake_r(:)=0
1098
1099	DO ij=1,nbpt_r
1100	IF (is_riverflow_r(ij) .OR. is_coastalflow_r(ij) ) THEN
1101	IF (frac_coast_r(ij)==0 .AND. frac_lake_r(ij)==0) THEN
1102	PRINT*,"riverflow or costalflow point on routing grid can not be interpolated on orchidee grid, this migh not happen"
1103	STOP
1104	ELSE IF (frac_coast_r(ij)==0 .AND. frac_lake_r(ij)>0) THEN
1105	weight_coast_to_lake_r(ij) = 1./frac_lake_r(ij)
1106	weight_coast_to_coast_r(ij) = 0
1107	ELSE IF (frac_coast_r(ij)>0 ) THEN
1108	weight_coast_to_coast_r(ij) = 1./(frac_coast_r(ij) + frac_lake_r(ij))
1109	weight_coast_to_lake_r(ij) = 0
1110	ENDIF
1111	ELSE IF (is_lakeinflow_r(ij)) THEN
1112	IF ( frac_coast_r(ij)==0 .AND. frac_lake_r(ij)==0) THEN
1113	PRINT*,"lakeinflow on routing grid can not be interpolated on orchidee grid, this might not happen"
1114	STOP
1115	ELSE IF (frac_lake_r(ij)==0 .AND. frac_coast_r(ij)>0) THEN
1116	weight_lake_to_coast_r(ij) = 1./frac_coast_r(ij)
1117	weight_lake_to_lake_r(ij) = 0
1118	ELSE IF (frac_lake_r(ij)>0 ) THEN
1119	weight_lake_to_lake_r(ij) = 1./(frac_coast_r(ij) + frac_lake_r(ij))
1120	weight_lake_to_coast_r(ij) = 0
1121	ENDIF
1122	ENDIF
1123	ENDDO
1124
1125	CALL compute_basins_area(nbpt_r,nbpt_rp1)
1126
1127	ELSE
1128
1129	nbpt_r=0
1130
1131	ENDIF !is_omp_root
1132
1133
1134	END SUBROUTINE routing_flow_init_local
1135
1136
1137	SUBROUTINE compute_basins_area(nbpt_r,nbpt_rp1 )
1138	USE xios
1139	USE routing_native_para
1140	IMPLICIT NONE
1141	INCLUDE 'mpif.h'
1142	INTEGER, INTENT(IN) :: nbpt_r
1143	INTEGER, INTENT(IN) :: nbpt_rp1
1144	REAL(r_std),PARAMETER :: PI=atan(1.)*4
1145	REAL(r_std),PARAMETER :: earth_radius=6.371009E6 ! in meter
1146
1147	INTEGER :: ni, nj
1148	REAL(r_std),ALLOCATABLE :: lon(:), lat(:)
1149	REAL(r_std) :: area_rp1(nbpt_rp1)
1150	REAL(r_std) :: send_area_rp1(nbpt_rp1)
1151	REAL(r_std) :: recv_area_r(nbpt_rp1)
1152	REAL :: delta_lon, delta_lat
1153	INTEGER :: i,j,ij,it,ig,ierr
1154	REAL(r_std) :: sum_area
1155
1156	ALLOCATE(basins_area_r(nbpt_r))
1157
1158	CALL xios_get_domain_attr("routing_domain", ni=ni, nj=nj)
1159	ALLOCATE(lon(ni),lat(nj))
1160	CALL xios_get_domain_attr("routing_domain",lonvalue_1d=lon, latvalue_1d=lat)
1161
1162	! compute routing cell area
1163	delta_lon=ABS(lon(2)-lon(1))*Pi/180.
1164	delta_lat=ABS(lat(2)-lat(1))*Pi/180.
1165
1166	DO j=1,nj
1167	DO i=1,ni
1168	ig=(j+1-1)*(ni+2)+i+1
1169	area_rp1(ig) = delta_londelta_latcos(lat(j)Pi/180.)earth_radius*earth_radius ! pretty good approximation
1170	ij = (j-1)*ni+i
1171	basins_area_r(ij) = delta_londelta_latcos(lat(j)Pi/180.)earth_radius*earth_radius ! pretty good approximation
1172	ENDDO
1173	ENDDO
1174
1175	CALL update_halo(area_rp1)
1176
1177
1178	send_area_rp1(:) = area_rp1(:)
1179	sum_area=1.
1180	DO WHILE (sum_area/=0)
1181
1182	recv_area_r(:) = 0
1183	DO ig=1,nbpt_rp1
1184	IF ( route_flow_rp1(ig) > 0 ) THEN
1185	basins_area_r(route_flow_rp1(ig)) = basins_area_r(route_flow_rp1(ig)) + send_area_rp1(ig)
1186	recv_area_r(route_flow_rp1(ig)) = recv_area_r(route_flow_rp1(ig)) + send_area_rp1(ig)
1187	ENDIF
1188	ENDDO
1189
1190	DO j=1,nj
1191	DO i=1,ni
1192	ij=((j-1)*ni)+i
1193	ig=(j+1-1)*(ni+2)+i+1
1194	send_area_rp1(ig) = recv_area_r(ij)
1195	IF (is_riverflow_r(ij) .OR. is_coastalflow_r(ij) .OR. is_lakeinflow_r(ij)) send_area_rp1(ig)=0
1196	ENDDO
1197	ENDDO
1198
1199	CALL update_halo(send_area_rp1)
1200	sum_area = sum(send_area_rp1)
1201	CALL MPI_ALLREDUCE(MPI_IN_PLACE, sum_area, 1, MPI_REAL_ORCH, MPI_SUM, MPI_COMM_ORCH, ierr)
1202
1203	ENDDO
1204
1205	END SUBROUTINE compute_basins_area
1206
1207
1208	SUBROUTINE initialize_stations(dt_routing)
1209	USE xios
1210	IMPLICIT NONE
1211	INCLUDE 'mpif.h'
1212	REAL(r_std), INTENT (in) :: dt_routing !! Routing time step (s)
1213
1214	REAL(r_std),PARAMETER :: PI=atan(1.)*4
1215	REAL(r_std),PARAMETER :: earth_radius=6.371009E6 ! in meter
1216	INTEGER :: ni, nj
1217	REAL(r_std),ALLOCATABLE :: lon(:), lat(:)
1218
1219	INTEGER :: nb_mouth
1220	CHARACTER(LEN=60),ALLOCATABLE :: mouth(:)
1221	INTEGER,ALLOCATABLE :: mouth_id(:)
1222	INTEGER,ALLOCATABLE :: mouth_index(:)
1223
1224	REAL,ALLOCATABLE :: station_lonlat(:,:)
1225	INTEGER,ALLOCATABLE :: station_prec(:)
1226
1227	INTEGER ::i,j,r,k
1228
1229	TYPE(xios_duration) :: ts
1230	TYPE(xios_scalar) :: scalar_hdl
1231	TYPE(xios_scalargroup) :: scalar_def_hdl
1232	TYPE(xios_field) :: field_hdl
1233	TYPE(xios_fieldgroup) :: field_def_hdl, fieldgroup_hdl
1234	TYPE(xios_file) :: file_hdl
1235	TYPE(xios_filegroup) :: file_def_hdl
1236	TYPE(xios_date) :: start_date, time_origin
1237	CHARACTER(LEN=20) :: calendar_type
1238
1239	CHARACTER(LEN=256) :: str_station_ind
1240	REAL :: lon_a,lat_a,lon_b,lat_b,dist,max_area, max_max_area, min_dist
1241	INTEGER :: min_dist_index
1242	INTEGER :: ierr
1243
1244
1245	CALL xios_get_domain_attr("routing_domain", ni=ni, nj=nj)
1246	ALLOCATE(lon(ni),lat(nj))
1247	CALL xios_get_domain_attr("routing_domain",lonvalue_1d=lon, latvalue_1d=lat)
1248
1249	nb_station=0
1250	CALL getin("nb_station",nb_station)
1251	ALLOCATE(station(nb_station))
1252	ALLOCATE(station_lonlat(nb_station,2))
1253	ALLOCATE(station_index(nb_station))
1254	ALLOCATE(station_prec(nb_station))
1255
1256	DO k=1,nb_station
1257	WRITE(str_station_ind,*) k
1258	str_station_ind=ADJUSTL(str_station_ind)
1259	CALL getin("station"//TRIM(str_station_ind)//"_id",station(k))
1260	CALL getin("station"//TRIM(str_station_ind)//"_coor",station_lonlat(k,:))
1261	station_prec(k)=50000
1262	CALL getin("station"//TRIM(str_station_ind)//"_prec",station_prec(k))
1263
1264	lon_a = station_lonlat(k,1)*Pi/180.
1265	lat_a = station_lonlat(k,2)*Pi/180.
1266
1267	max_area=0
1268	min_dist=HUGE(min_dist)
1269	DO j=1,nj
1270	DO i=1,ni
1271	r=((j-1)*ni)+i
1272	IF (routing_mask_r(r)) THEN
1273	lon_b=lon(i)*Pi/180.
1274	lat_b=lat(j)*Pi/180.
1275	dist=earth_radiusacos(sin(lon_a)sin(lon_b)+cos(lon_a)cos(lon_b)cos(lat_b-lat_a))
1276	IF (dist<station_prec(k)) THEN
1277	IF (dist<min_dist) THEN
1278	min_dist=dist
1279	min_dist_index = r
1280	ENDIF
1281
1282	IF (basins_area_r(r) > max_area) THEN
1283	max_area = basins_area_r(r)
1284	station_index(k) = r
1285	ENDIF
1286	ENDIF
1287	ENDIF
1288	ENDDO
1289	ENDDO
1290
1291	CALL MPI_ALLREDUCE(max_area, max_max_area, 1, MPI_REAL_ORCH, MPI_MAX, MPI_COMM_ORCH, ierr)
1292	IF (max_area /= max_max_area) station_index(k)=-1
1293	IF (max_area == 0) station_index(k)=-1
1294
1295	! PRINT*,"Station ",TRIM(station(k))," old coor :",lon(min_dist_index),lat(min_dist_index)," basin_area (km2)", basins_area_r(min_dist_index)/1000/1000
1296	! PRINT*,"Station ",TRIM(station(k)),"new coor ",lon(station_index(k)),lat(station_index(k))," basin_area(km2) ", basins_area_r(station_index(k))/1000/1000
1297	ENDDO
1298
1299	CALL xios_get_calendar_type(calendar_type)
1300	CALL xios_get_start_date(start_date)
1301	CALL xios_get_time_origin(time_origin)
1302
1303	CALL xios_context_initialize("orchidee_routing_out",MPI_COMM_ORCH)
1304	CALL xios_orchidee_change_context("orchidee_routing_out")
1305	ts.second=dt_routing
1306	calendar_type="gregorian"
1307	CALL xios_define_calendar(type=calendar_type,start_date=start_date,time_origin=time_origin, timestep=ts )
1308
1309	CALL xios_get_handle("scalar_definition",scalar_def_hdl)
1310
1311	! CALL xios_get_handle("mouths",field_def_hdl)
1312	! DO k=1,nb_mouth
1313	! CALL xios_add_child(scalar_def_hdl,scalar_hdl,TRIM(mouth(k))//"_scalar")
1314	! CALL xios_set_attr(scalar_hdl,label=TRIM(mouth(k)))
1315	! CALL xios_add_child(field_def_hdl,field_hdl, TRIM(mouth(k)))
1316	! CALL xios_set_attr(field_hdl,scalar_ref=TRIM(mouth(k))//"_scalar")
1317	! CALL xios_set_attr(field_hdl, freq_op = freq_op*xios_second)
1318	! ENDDO
1319
1320	CALL xios_get_handle("field_definition",field_def_hdl)
1321	CALL xios_add_child(field_def_hdl, fieldgroup_hdl, "stations")
1322	DO k=1,nb_station
1323	CALL xios_add_child(scalar_def_hdl,scalar_hdl,TRIM(station(k))//"_scalar")
1324	CALL xios_set_attr(scalar_hdl,label=TRIM(station(k)))
1325	CALL xios_add_child(fieldgroup_hdl,field_hdl, TRIM(station(k)))
1326	CALL xios_set_attr(field_hdl,scalar_ref=TRIM(station(k))//"_scalar")
1327	ENDDO
1328
1329	CALL xios_get_handle("file_definition",file_def_hdl)
1330	CALL xios_add_child(file_def_hdl, file_hdl, "stations")
1331	CALL xios_set_attr(file_hdl, type="one_file", output_freq=ts, sync_freq=ts, enabled=.TRUE.)
1332	CALL xios_add_child(file_hdl, fieldgroup_hdl)
1333	CALL xios_set_attr(fieldgroup_hdl, group_ref="stations", operation="average")
1334	CALL xios_close_context_definition()
1335	CALL xios_orchidee_change_context("orchidee")
1336	END SUBROUTINE initialize_stations
1337
1338	SUBROUTINE routing_flow_main(dt_routing) !
1339
1340	USE xios
1341	USE grid, ONLY : area
1342	USE routing_native_para
1343	IMPLICIT NONE
1344	INCLUDE "mpif.h"
1345
1346	!! 0 Variable and parameter description
1347	!! 0.1 Input variables
1348	REAL(r_std), INTENT (in) :: dt_routing !! Routing time step (s)
1349
1350	!! 0.4 Local variables
1351	REAL(r_std) :: runoff(nbp_mpi) !! Grid-point runoff (kg/dt)
1352	REAL(r_std) :: runoff_in(nbp_mpi) !! Grid-point runoff (kg/dt)
1353	REAL(r_std) :: drainage(nbp_mpi) !! Grid-point drainage (kg/dt)
1354	REAL(r_std) :: drainage_in(nbp_mpi) !! Grid-point drainage (kg/dt)
1355	REAL(r_std) :: riverflow(nbp_mpi)
1356	REAL(r_std) :: coastalflow(nbp_mpi)
1357	REAL(r_std) :: lakeinflow(nbp_mpi)
1358	REAL(r_std) :: fast_diag_mpi(nbp_mpi)
1359	REAL(r_std) :: slow_diag_mpi(nbp_mpi)
1360	REAL(r_std) :: stream_diag_mpi(nbp_mpi)
1361	REAL(r_std) :: area_mpi(nbp_mpi) ! cell area
1362	REAL(r_std) :: flow_coast(nbp_mpi)
1363	REAL(r_std) :: flow_lake(nbp_mpi)
1364
1365	! from input model -> routing_grid
1366	REAL(r_std) :: runoff_r(nbpt_r) !! Grid-point runoff (kg/m^2/dt)
1367	REAL(r_std) :: drainage_r(nbpt_r) !! Grid-point drainage (kg/m^2/dt)
1368
1369	REAL(r_std), DIMENSION(nbpt_r) :: fast_flow_r !! Outflow from the fast reservoir (kg/dt)
1370	REAL(r_std), DIMENSION(nbpt_r) :: slow_flow_r !! Outflow from the slow reservoir (kg/dt)
1371	REAL(r_std), DIMENSION(nbpt_r) :: stream_flow_r !! Outflow from the stream reservoir (kg/dt)
1372	REAL(r_std), DIMENSION(nbpt_r) :: hydrographs_r !! hydrograph (kg/dt)
1373	REAL(r_std), DIMENSION(nbpt_r) :: transport_r !! Water transport between basins (kg/dt)
1374
1375	INTEGER(i_std) :: ig !! Indices (unitless)
1376	INTEGER(i_std) :: isplit
1377
1378	LOGICAL, PARAMETER :: check_reservoir = .TRUE. !! Logical to choose if we write informations when a negative amount of water is occurring in a reservoir (true/false)
1379
1380	REAL(r_std), DIMENSION(nbpt_r) :: lakeinflow_r
1381	REAL(r_std), DIMENSION(nbpt_r) :: coastalflow_r
1382	REAL(r_std), DIMENSION(nbpt_r) :: riverflow_r
1383
1384	REAL(r_std), DIMENSION(nbpt_r) :: flow_r
1385	REAL(r_std), DIMENSION(nbpt_rp1) :: flow_rp1
1386	REAL(r_std) :: flow !! Outflow computation for the reservoirs (kg/dt)
1387	REAL(r_std) :: basins_riverflow_mpi(0:basins_count)
1388	REAL(r_std) :: basins_riverflow(0:basins_count)
1389	REAL(r_std) :: water_balance_before, water_balance_after
1390	REAL(r_std) :: sum_water_before, sum_water_after
1391	REAL(r_std) :: value
1392	INTEGER(i_std) :: k
1393	INTEGER :: ierr
1394
1395	!_ ================================================================================================================================
1396	!> The outflow fluxes from the three reservoirs are computed.
1397	!> The outflow of volume of water Vi into the reservoir i is assumed to be linearly related to its volume.
1398	!> The water travel simulated by the routing scheme is dependent on the water retention index topo_resid
1399	!> given by a 0.5 degree resolution map for each pixel performed from a simplification of Manning's formula
1400	!> (Dingman, 1994; Ducharne et al., 2003).
1401	!> The resulting product of tcst (in day/m) and topo_resid (in m) represents the time constant (day)
1402	!> which is an e-folding time, the time necessary for the water amount
1403	!> in the stream reservoir to decrease by a factor e. Hence, it gives an order of
1404	!> magnitude of the travel time through this reservoir between
1405	!> the sub-basin considered and its downstream neighbor.
1406
1407
1408
1409	CALL gather_omp(runoff_mean,runoff)
1410	CALL gather_omp(drainage_mean, drainage)
1411	CALL gather_omp(area, area_mpi)
1412
1413
1414	IF (is_omp_root) THEN
1415	CALL xios_send_field("routing_basins_area",basins_area_r)
1416
1417	hydrographs_r(:)=0
1418	! water balance before
1419
1420	sum_water_before = sum((runoff(:)+drainage(:))*routing_weight(:))+sum(fast_reservoir_r(:)+slow_reservoir_r(:)+stream_reservoir_r(:))
1421	CALL MPI_ALLREDUCE(sum_water_before, water_balance_before,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr)
1422
1423
1424	runoff_in=runoff*routing_weight_in
1425	CALL xios_send_field("routing_runoff",runoff_in) ! interp conservative model -> routing
1426	CALL xios_recv_field("routing_runoff_r",runoff_r)
1427	WHERE(.NOT. routing_mask_r) runoff_r = 0
1428
1429	drainage_in=drainage*routing_weight_in
1430	CALL xios_send_field("routing_drainage",drainage_in) ! interp conservative model -> routing
1431	CALL xios_recv_field("routing_drainage_r",drainage_r)
1432	WHERE(.NOT. routing_mask_r) drainage_r = 0
1433
1434	CALL MPI_ALLREDUCE(sum(runoff*routing_weight),sum_water_before,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr)
1435	CALL MPI_ALLREDUCE(sum(runoff_r),sum_water_after,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr)
1436	IF (sum_water_after/=0) THEN
1437	IF (is_mpi_root) PRINT *,"runoff fixer : ", sum_water_before/sum_water_after
1438	DO ig=1,nbpt_r
1439	runoff_r(ig) = runoff_r(ig) * sum_water_before/sum_water_after
1440	ENDDO
1441	ENDIF
1442
1443	CALL MPI_ALLREDUCE(sum(drainage*routing_weight),sum_water_before,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr)
1444	CALL MPI_ALLREDUCE(sum(drainage_r),sum_water_after,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr)
1445	IF (sum_water_after/=0) THEN
1446	IF (is_mpi_root) PRINT *,"drainage fixer : ", sum_water_before/sum_water_after
1447	DO ig=1,nbpt_r
1448	drainage_r(ig) = drainage_r(ig) * sum_water_before/sum_water_after
1449	ENDDO
1450	ENDIF
1451
1452
1453	CALL MPI_ALLREDUCE(sum(runoffrouting_weight)+sum(drainagerouting_weight),sum_water_before,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr)
1454	CALL MPI_ALLREDUCE(sum(runoff_r) +sum(drainage_r),sum_water_after,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr)
1455
1456	IF (is_mpi_root) PRINT *,"Loose water by interpolation ; before : ", sum_water_before," ; after : ",sum_water_after, &
1457	" ; delta : ", 100.(sum_water_after-sum_water_before)/(0.5(sum_water_after+sum_water_before)),"%"
1458
1459
1460	runoff_in(:) = runoff(:)*unrouted_weight(:)
1461	drainage_in(:) = drainage(:)*unrouted_weight(:)
1462	runoff_r(:) = runoff_r(:) / split_routing
1463	drainage_r(:) = drainage_r(:) / split_routing
1464	hydrographs_r(:) = 0
1465	transport_r(:)=0
1466
1467	DO isplit=1,split_routing
1468
1469	DO ig=1,nbpt_r
1470	IF ( routing_mask_r(ig) ) THEN
1471
1472	flow = MIN(fast_reservoir_r(ig)/((topoind_r(ig)/1000.)fast_tcstone_day/(dt_routing/split_routing)),&
1473	& fast_reservoir_r(ig)-min_sechiba)
1474	fast_flow_r(ig) = MAX(flow, zero)
1475	!
1476	flow = MIN(slow_reservoir_r(ig)/((topoind_r(ig)/1000.)slow_tcstone_day/(dt_routing/split_routing)),&
1477	& slow_reservoir_r(ig)-min_sechiba)
1478	slow_flow_r(ig) = MAX(flow, zero)
1479	!
1480	flow = MIN(stream_reservoir_r(ig)/((topoind_r(ig)/1000.)stream_tcstone_day/(dt_routing/split_routing)),&
1481	& stream_reservoir_r(ig)-min_sechiba)
1482	stream_flow_r(ig) = MAX(flow, zero)
1483	!
1484	!
1485	ELSE
1486	fast_flow_r(ig) = zero
1487	slow_flow_r(ig) = zero
1488	stream_flow_r(ig) = zero
1489	ENDIF
1490	ENDDO
1491
1492
1493	!-
1494	!- Compute the transport
1495	!-
1496
1497	DO ig=1,nbpt_r
1498	flow_rp1(r_to_rp1(ig))=fast_flow_r(ig) + slow_flow_r(ig) + stream_flow_r(ig)
1499	ENDDO
1500
1501	CALL update_halo(flow_rp1) ! transfert halo
1502
1503	DO ig=1,nbpt_rp1
1504	IF ( route_flow_rp1(ig) > 0 ) THEN
1505	transport_r(route_flow_rp1(ig))=transport_r(route_flow_rp1(ig))+ flow_rp1(ig)
1506	ENDIF
1507	ENDDO
1508
1509
1510	DO ig=1,nbpt_r
1511	IF ( routing_mask_r(ig) ) THEN
1512	fast_reservoir_r(ig) = fast_reservoir_r(ig) + runoff_r(ig) - fast_flow_r(ig)
1513	slow_reservoir_r(ig) = slow_reservoir_r(ig) + drainage_r(ig) - slow_flow_r(ig)
1514	stream_reservoir_r(ig) = stream_reservoir_r(ig) - stream_flow_r(ig)
1515	IF (is_streamflow_r(ig)) THEN
1516	stream_reservoir_r(ig)= stream_reservoir_r(ig) + transport_r(ig)
1517	transport_r(ig) = 0
1518	ENDIF
1519
1520	IF ( stream_reservoir_r(ig) .LT. zero ) THEN
1521	IF ( check_reservoir ) THEN
1522	WRITE(numout,*) "WARNING : negative stream reservoir at :", ig, ". Problem is being corrected."
1523	WRITE(numout,*) "stream_reservoir, transport, stream_flow,: ", &
1524	stream_reservoir_r(ig), transport_r(ig), stream_flow_r(ig)
1525	ENDIF
1526	fast_reservoir_r(ig) = fast_reservoir_r(ig) + stream_reservoir_r(ig)
1527	stream_reservoir_r(ig) = zero
1528	ENDIF
1529	!
1530	IF ( fast_reservoir_r(ig) .LT. zero ) THEN
1531	IF ( check_reservoir ) THEN
1532	WRITE(numout,*) "WARNING : negative fast reservoir at :", ig, ". Problem is being corrected."
1533	WRITE(numout,*) "fast_reservoir, runoff, fast_flow : ", fast_reservoir_r(ig), &
1534	&runoff_r(ig), fast_flow_r(ig)
1535	ENDIF
1536	slow_reservoir_r(ig) = slow_reservoir_r(ig) + fast_reservoir_r(ig)
1537	fast_reservoir_r(ig) = zero
1538	ENDIF
1539
1540	IF ( slow_reservoir_r(ig) .LT. - min_sechiba ) THEN
1541	WRITE(numout,*) 'WARNING : There is a negative reservoir at :', ig
1542	WRITE(numout,*) 'WARNING : slowr, slow_flow, drainage', &
1543	& slow_reservoir_r(ig), slow_flow_r(ig), drainage_r(ig)
1544	CALL ipslerr_p(2, 'routing_simple_flow', 'WARNING negative slow_reservoir.','','')
1545	ENDIF
1546	ENDIF
1547	ENDDO
1548
1549	DO ig=1,nbpt_r
1550	IF ( routing_mask_r(ig) ) THEN
1551	hydrographs_r(ig)=hydrographs_r(ig)+fast_flow_r(ig)+slow_flow_r(ig)+stream_flow_r(ig)
1552	ENDIF
1553	ENDDO
1554
1555	ENDDO ! isplit
1556
1557	lakeinflow_r(:)=0
1558	coastalflow_r(:)=0
1559	riverflow_r(:)=0
1560	basins_riverflow_mpi(:)=0
1561
1562	DO ig=1,nbpt_r
1563	IF ( routing_mask_r(ig) ) THEN
1564	IF (is_lakeinflow_r(ig)) THEN
1565	lakeinflow_r(ig) = transport_r(ig)
1566	basins_riverflow_mpi(basins_extended_r(ig)) = basins_riverflow_mpi(basins_extended_r(ig))+lakeinflow_r(ig)
1567	ENDIF
1568
1569	IF (is_coastalflow_r(ig)) THEN
1570	coastalflow_r(ig) = transport_r(ig)
1571	basins_riverflow_mpi(basins_extended_r(ig)) = &
1572	basins_riverflow_mpi(basins_extended_r(ig))+coastalflow_r(ig)
1573	ENDIF
1574
1575	IF (is_riverflow_r(ig)) THEN
1576	riverflow_r(ig) = transport_r(ig)
1577	basins_riverflow_mpi(basins_extended_r(ig)) = &
1578	basins_riverflow_mpi(basins_extended_r(ig))+riverflow_r(ig)
1579	ENDIF
1580	ENDIF
1581	ENDDO
1582
1583	! now send riverflow, coastalflow and lakeinflow on orchidee grid
1584
1585	coastalflow(:)=0.
1586	riverflow(:)=0.
1587	lakeinflow(:)=0.
1588
1589	CALL xios_send_field("routing_coastalflow_to_coast_r" ,coastalflow_r*weight_coast_to_coast_r)
1590	CALL xios_recv_field("routing_coastalflow_to_coast" ,flow_coast)
1591	WHERE(.NOT.is_coastline) flow_coast=0
1592
1593	CALL xios_send_field("routing_coastalflow_to_lake_r" ,coastalflow_r*weight_coast_to_lake_r)
1594	CALL xios_recv_field("routing_coastalflow_to_lake" ,flow_lake)
1595	WHERE(is_coastline) flow_lake=0.
1596
1597	CALL MPI_ALLREDUCE(sum(coastalflow_r),sum_water_before,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr)
1598	CALL MPI_ALLREDUCE(sum(flow_coast+flow_lake),sum_water_after,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr)
1599	IF (sum_water_after/=0) THEN
1600	IF (is_mpi_root) PRINT *,"coastalflow fixer : ", sum_water_before/sum_water_after
1601	flow_coast(:)=flow_coast(:)*(sum_water_before/sum_water_after)
1602	flow_lake(:)=flow_lake(:)*(sum_water_before/sum_water_after)
1603	ENDIF
1604
1605	coastalflow=coastalflow+flow_coast
1606	lakeinflow=lakeinflow+flow_lake
1607
1608	CALL xios_send_field("routing_riverflow_to_coast_r" ,riverflow_r*weight_coast_to_coast_r)
1609	CALL xios_recv_field("routing_riverflow_to_coast" ,flow_coast)
1610	WHERE(.NOT.is_coastline) flow_coast=0
1611
1612	CALL xios_send_field("routing_riverflow_to_lake_r" ,riverflow_r*weight_coast_to_lake_r)
1613	CALL xios_recv_field("routing_riverflow_to_lake" ,flow_lake)
1614	WHERE(is_coastline) flow_lake=0.
1615
1616	CALL MPI_ALLREDUCE(sum(riverflow_r),sum_water_before,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr)
1617	CALL MPI_ALLREDUCE(sum(flow_coast+flow_lake),sum_water_after,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr)
1618	IF (sum_water_after/=0) THEN
1619	IF (is_mpi_root) PRINT *,"riverflow fixer : ", sum_water_before/sum_water_after
1620	flow_coast(:)=flow_coast(:)*(sum_water_before/sum_water_after)
1621	flow_lake(:)=flow_lake(:)*(sum_water_before/sum_water_after)
1622	ENDIF
1623
1624	riverflow=riverflow+flow_coast
1625	lakeinflow=lakeinflow+flow_lake
1626
1627	CALL xios_send_field("routing_lakeinflow_to_coast_r" ,lakeinflow_r*weight_lake_to_coast_r)
1628	CALL xios_recv_field("routing_lakeinflow_to_coast" ,flow_coast)
1629	WHERE(.NOT.is_coastline) flow_coast=0
1630
1631	CALL xios_send_field("routing_lakeinflow_to_lake_r" ,lakeinflow_r*weight_lake_to_lake_r)
1632	CALL xios_recv_field("routing_lakeinflow_to_lake" ,flow_lake)
1633	WHERE(is_coastline) flow_lake=0.
1634
1635	CALL MPI_ALLREDUCE(sum(lakeinflow_r),sum_water_before,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr)
1636	CALL MPI_ALLREDUCE(sum(flow_coast+flow_lake),sum_water_after,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr)
1637	IF (sum_water_after/=0) THEN
1638	IF (is_mpi_root) PRINT *,"lakeinflow fixer : ", sum_water_before/sum_water_after
1639	flow_coast(:)=flow_coast(:)*(sum_water_before/sum_water_after)
1640	flow_lake(:)=flow_lake(:)*(sum_water_before/sum_water_after)
1641	ENDIF
1642
1643	coastalflow=coastalflow+flow_coast
1644	lakeinflow=lakeinflow+flow_lake
1645
1646
1647	WHERE(is_coastline) coastalflow = coastalflow + runoff_in + drainage_in
1648	WHERE(.NOT.is_coastline) lakeinflow = lakeinflow + runoff_in + drainage_in
1649
1650	sum_water_after = sum(coastalflow(:) + riverflow(:) + lakeinflow(:))+sum(fast_reservoir_r(:)+slow_reservoir_r(:)+stream_reservoir_r(:))
1651	CALL MPI_ALLREDUCE(sum_water_after, water_balance_after,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr)
1652
1653	IF (is_mpi_root) PRINT *,"routing water Balance ; before : ", water_balance_before," ; after : ",water_balance_after, &
1654	" ; delta : ", 100(water_balance_after-water_balance_before)/(0.5(water_balance_after+water_balance_before)),"%"
1655	! diag
1656	CALL xios_send_field("routing_fast_reservoir_r" , fast_reservoir_r)
1657	CALL xios_send_field("routing_slow_reservoir_r" , slow_reservoir_r)
1658	CALL xios_send_field("routing_stream_reservoir_r" , stream_reservoir_r)
1659	CALL xios_send_field("routing_riverflow_r" , riverflow_r)
1660	CALL xios_send_field("routing_coastalflow_r" , coastalflow_r)
1661	CALL xios_send_field("routing_lakeinflow_r" , lakeinflow_r)
1662	CALL xios_send_field("out_flow",lakeinflow+coastalflow+riverflow)
1663	CALL xios_send_field("routing_hydrographs_r", (hydrographs_r+lakeinflow_r+coastalflow_r+riverflow_r)/1000./dt_routing)
1664	CALL xios_send_field("routing_riverflow" , riverflow)
1665	CALL xios_send_field("routing_coastalflow" , coastalflow)
1666	CALL xios_send_field("routing_lakeinflow" , lakeinflow)
1667	!
1668	! stations
1669	CALL xios_orchidee_change_context("orchidee_routing_out")
1670	station_ts = station_ts +1
1671	CALL xios_update_calendar(station_ts)
1672	DO k=1,nb_station
1673	IF (station_index(k) /=-1) THEN
1674	value = hydrographs_r(station_index(k))/1000./dt_routing
1675	ELSE
1676	value = 0
1677	ENDIF
1678	CALL MPI_ALLREDUCE(MPI_IN_PLACE,value,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr)
1679	CALL xios_send_field(TRIM(station(k)),value)
1680	ENDDO
1681	CALL xios_orchidee_change_context("orchidee")
1682	ENDIF ! is_omp_root
1683
1684	CALL scatter_omp(riverflow,riverflow_mean)
1685	CALL scatter_omp(coastalflow,coastalflow_mean)
1686	CALL scatter_omp(lakeinflow,lakeinflow_mean)
1687
1688	CALL routing_flow_reset_mean
1689
1690	END SUBROUTINE routing_flow_main
1691
1692
1693	!! =============================================================================================================================
1694	!! SUBROUTINE: routing_simple_finalize
1695	!!
1696	!>\BRIEF Write to restart file
1697	!!
1698	!! DESCRIPTION: Write module variables to restart file
1699	!!
1700	!! RECENT CHANGE(S)
1701	!!
1702	!! REFERENCE(S)
1703	!!
1704	!! FLOWCHART
1705	!! \n
1706	!_ ==============================================================================================================================
1707
1708	SUBROUTINE routing_flow_finalize(kjit, rest_id)
1709	USE xios
1710	USE ioipsl
1711	IMPLICIT NONE
1712	INTEGER, INTENT(IN) :: kjit
1713	INTEGER, INTENT(IN) :: rest_id
1714	!_ ================================================================================================================================
1715
1716	IF (is_omp_root) THEN
1717	CALL xios_send_field("fast_reservoir_restart",fast_reservoir_r)
1718	CALL xios_send_field("slow_reservoir_restart",slow_reservoir_r)
1719	CALL xios_send_field("stream_reservoir_restart",stream_reservoir_r)
1720	ENDIF
1721
1722	CALL restput_p (rest_id, 'riverflow', nbp_glo, 1, 1, kjit, riverflow_mean, 'scatter', nbp_glo, index_g)
1723	CALL restput_p (rest_id, 'coastalflow', nbp_glo, 1, 1, kjit, coastalflow_mean, 'scatter', nbp_glo, index_g)
1724	CALL restput_p (rest_id, 'lakeinflow', nbp_glo, 1, 1, kjit, lakeinflow_mean, 'scatter', nbp_glo, index_g)
1725	CALL routing_flow_finalize_mean(kjit, rest_id)
1726
1727	END SUBROUTINE routing_flow_finalize
1728
1729
1730	!! ================================================================================================================================
1731	!! SUBROUTINE : routing_simple_clear
1732	!!
1733	!>\BRIEF This subroutine deallocates the block memory previously allocated.
1734	!!
1735	!! DESCRIPTION: This subroutine deallocates the block memory previously allocated.
1736	!!
1737	!! RECENT CHANGE(S): None
1738	!!
1739	!! MAIN OUTPUT VARIABLE(S):
1740	!!
1741	!! REFERENCES : None
1742	!!
1743	!! FLOWCHART :None
1744	!! \n
1745	!_ ================================================================================================================================
1746
1747	SUBROUTINE routing_flow_clear
1748	IMPLICIT NONE
1749
1750	IF (is_omp_root) THEN
1751	IF (ALLOCATED(topoind_r)) DEALLOCATE(topoind_r)
1752	IF (ALLOCATED(route_flow_rp1)) DEALLOCATE(route_flow_rp1)
1753	IF (ALLOCATED(routing_mask_r)) DEALLOCATE(routing_mask_r)
1754	IF (ALLOCATED(fast_reservoir_r)) DEALLOCATE(fast_reservoir_r)
1755	IF (ALLOCATED(slow_reservoir_r)) DEALLOCATE(slow_reservoir_r)
1756	IF (ALLOCATED(is_lakeinflow_r)) DEALLOCATE(is_lakeinflow_r)
1757	IF (ALLOCATED(is_coastalflow_r)) DEALLOCATE(is_coastalflow_r)
1758	IF (ALLOCATED(is_riverflow_r)) DEALLOCATE(is_riverflow_r)
1759	ENDIF
1760
1761	END SUBROUTINE routing_flow_clear
1762
1763
1764	END MODULE routing_native_flow_mod

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

Context Navigation

source: branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/routing_native_flow.f90 @ 8227

Download in other formats: