Context Navigation

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

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

Some update for routing native :

stations can be now also selected using the basin area information
Add diagnostics on orchidee grid, in the standard sechiba_history file :
- hydrograph
- streamr
- fastr
- slowr
- irrigation

File size: 96.0 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 :: routing_flow_diags
21	PUBLIC :: compute_coastline
22
23	INTEGER,SAVE :: nbpt
24	!$OMP THREADPRIVATE(nbpt)
25
26	!! PARAMETERS
27	REAL(r_std), SAVE :: fast_tcst = 9e-3 !! Property of the fast reservoir (day/km)
28	!$OMP THREADPRIVATE(fast_tcst)
29	REAL(r_std), SAVE :: slow_tcst = 1.2e-2 !! Property of the slow reservoir (day/km)
30	!$OMP THREADPRIVATE(slow_tcst)
31	REAL(r_std), SAVE :: stream_tcst = 3e-5 !! Property of the stream reservoir (day/km)
32	!$OMP THREADPRIVATE(stream_tcst)
33
34
35	REAL(r_std),SAVE,ALLOCATABLE :: runoff_mean(:)
36	!$OMP THREADPRIVATE(runoff_mean)
37
38	REAL(r_std),SAVE,ALLOCATABLE :: drainage_mean(:)
39	!$OMP THREADPRIVATE(drainage_mean)
40
41	REAL(r_std),SAVE,ALLOCATABLE,PUBLIC,PROTECTED,TARGET :: fast_reservoir_r(:) !! Water amount in the fast reservoir (kg) - (local routing grid)
42	!$OMP THREADPRIVATE(fast_reservoir_r)
43
44	REAL(r_std),SAVE,ALLOCATABLE,PUBLIC,PROTECTED,TARGET:: slow_reservoir_r(:) !! Water amount in the slow reservoir (kg) - (local routing grid)
45	!$OMP THREADPRIVATE(slow_reservoir_r)
46
47	REAL(r_std),SAVE,ALLOCATABLE,PUBLIC,PROTECTED,TARGET :: stream_reservoir_r(:) !! Water amount in the stream reservoir (kg) - (local routing grid)
48	!$OMP THREADPRIVATE(stream_reservoir_r)
49
50
51	REAL(r_std),SAVE,ALLOCATABLE,PUBLIC,PROTECTED,TARGET :: fast_diag(:) !! Diag on water amount in the fast reservoir (kg/m^2) - (orchidee grid)
52	!$OMP THREADPRIVATE(fast_diag)
53
54	REAL(r_std),SAVE,ALLOCATABLE,PUBLIC,PROTECTED,TARGET:: slow_diag(:) !! Diag on water amount in the slow reservoir (kg/m^2) - (orchidee grid)
55	!$OMP THREADPRIVATE(slow_diag)
56
57	REAL(r_std),SAVE,ALLOCATABLE,PUBLIC,PROTECTED,TARGET :: stream_diag(:) !! Diag on water amount in the stream reservoir (kg/m^2) - (orchidee grid)
58	!$OMP THREADPRIVATE(stream_diag)
59
60	REAL(r_std),SAVE,ALLOCATABLE,PUBLIC,PROTECTED,TARGET :: hydrographs_diag(:) !! Diag on water amount in the stream reservoir (kg/m^2) - (orchidee grid)
61	!$OMP THREADPRIVATE(hydrographs_diag)
62
63
64	REAL(r_std),SAVE,ALLOCATABLE,PUBLIC,PROTECTED :: riverflow_mean(:) !! Water amount in the stream reservoir (kg) - (local routing grid)
65	!$OMP THREADPRIVATE(riverflow_mean)
66
67	REAL(r_std),SAVE,ALLOCATABLE,PUBLIC,PROTECTED :: coastalflow_mean(:) !! Water amount in the stream reservoir (kg) - (local routing grid)
68	!$OMP THREADPRIVATE(coastalflow_mean)
69
70	REAL(r_std),SAVE,ALLOCATABLE,PUBLIC,PROTECTED :: lakeinflow_mean(:) !! Water amount in the stream reservoir (kg) - (local routing grid)
71	!$OMP THREADPRIVATE(lakeinflow_mean)
72
73	!
74	! Specific variables for simple routing
75	!
76	REAL(r_std),SAVE,ALLOCATABLE :: topoind_r(:) !! Topographic index of the retention time (m) index - (local routing grid)
77	!$OMP THREADPRIVATE(topoind_r)
78	REAL(r_std),SAVE :: topoind_factor !! conversion factor for topographic index (merit : m, standard : km), topo index must be in m
79	!$OMP THREADPRIVATE(topoind_factor)
80	INTEGER,SAVE,ALLOCATABLE :: route_flow_rp1(:) !! flow index from cell to neighboring cell following the trip direction - (local routing grid + halo)
81	!$OMP THREADPRIVATE(route_flow_rp1)
82	LOGICAL,SAVE,ALLOCATABLE :: is_lakeinflow_r(:) !! is lake inflow point - (local routing grid)
83	!$OMP THREADPRIVATE(is_lakeinflow_r)
84	LOGICAL,SAVE,ALLOCATABLE :: is_coastalflow_r(:) !! is coastal flow point - (local routing grid)
85	!$OMP THREADPRIVATE(is_coastalflow_r)
86	LOGICAL,SAVE,ALLOCATABLE :: is_riverflow_r(:) !! is river flow point - (local routing grid)
87	!$OMP THREADPRIVATE(is_riverflow_r)
88	LOGICAL,SAVE,ALLOCATABLE :: is_coastline(:) !! is coastline point - (local native grid)
89	!$OMP THREADPRIVATE(is_coastline)
90	LOGICAL,SAVE,ALLOCATABLE :: is_streamflow_r(:) !! is stream flow point - (local routing grid)
91	!$OMP THREADPRIVATE(is_streamflow_r)
92	LOGICAL,SAVE,ALLOCATABLE,PUBLIC,PROTECTED :: routing_mask_r(:) !! valid routing point - (local routing grid)
93	!$OMP THREADPRIVATE(routing_mask_r)
94	LOGICAL,SAVE,ALLOCATABLE :: coast_mask(:) !! is a coast point - (local native grid)
95	!$OMP THREADPRIVATE(coast_mask)
96	INTEGER,SAVE :: total_coast_points !! global number of coast point - (local native grid)
97	!$OMP THREADPRIVATE(total_coast_points)
98	INTEGER,SAVE :: nbpt_r !! number of point in local routing grid
99	!$OMP THREADPRIVATE(nbpt_r)
100	INTEGER,SAVE :: nbpt_rp1 !! number of point in local routing grid with halo of 1
101	!$OMP THREADPRIVATE(nbpt_rp1)
102	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)
103	!$OMP THREADPRIVATE(routing_weight)
104	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)
105	!$OMP THREADPRIVATE(routing_weight_in)
106	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)
107	!$OMP THREADPRIVATE(unrouted_weight)
108
109	REAL(r_std),SAVE,ALLOCATABLE :: weight_coast_to_coast_r(:)
110	!$OMP THREADPRIVATE(weight_coast_to_coast_r)
111
112	REAL(r_std),SAVE,ALLOCATABLE :: weight_coast_to_lake_r(:)
113	!$OMP THREADPRIVATE(weight_coast_to_lake_r)
114
115	REAL(r_std),SAVE,ALLOCATABLE :: weight_lake_to_coast_r(:)
116	!$OMP THREADPRIVATE(weight_lake_to_coast_r)
117
118	REAL(r_std),SAVE,ALLOCATABLE :: weight_lake_to_lake_r(:)
119	!$OMP THREADPRIVATE(weight_lake_to_lake_r)
120
121	REAL(r_std),SAVE,ALLOCATABLE :: basins_area_r(:)
122	!$OMP THREADPRIVATE(basins_area_r)
123
124	!! when doing interpolation from local routing grid to local native grid (river flow+coastal flow)
125	REAL(r_std),SAVE,ALLOCATABLE :: basins_extended_r(:) !! basins riverflow id (local routing grid)
126	!$OMP THREADPRIVATE(basins_extended_r)
127	INTEGER(i_std) :: basins_count !! number of basins (local routing grid)
128	!$OMP THREADPRIVATE(basins_count)
129	INTEGER(i_std) :: basins_out !! number of basins to output for diag
130	!$OMP THREADPRIVATE(basins_out)
131
132	INTEGER(i_std) :: split_routing !! time spliting for routing
133	!$OMP THREADPRIVATE(split_routing)
134
135
136	INTEGER :: nb_station
137	CHARACTER(LEN=60),ALLOCATABLE :: station(:)
138	REAL,ALLOCATABLE :: station_lonlat(:,:)
139	INTEGER,ALLOCATABLE :: station_index(:)
140	INTEGER :: station_ts = 0
141
142
143	! INTEGER(i_std), PARAMETER :: nb_stations=14
144	! REAL(r_std),PARAMETER :: station_lon(nb_stations) = &
145	! (/ -162.8830, -90.9058, -55.5110, -49.3242, -133.7447, -63.6000, 28.7167, &
146	! 15.3000, 66.5300, 89.6700, 86.5000, 127.6500, 3.3833, 117.6200 /)
147	! REAL(r_std),PARAMETER :: station_lat(nb_stations) = &
148	! (/ 61.9340, 32.3150, -1.9470, -5.1281, 67.4583, 8.1500, 45.2167, &
149	! -4.3000, 66.5700, 25.1800, 67.4800, 70.7000, 11.8667, 30.7700 /)
150	! CHARACTER(LEN=17),PARAMETER :: station_name(nb_stations) = &
151	! (/ "Pilot station ", "Vicksburg ", "Obidos ", &
152	! "Itupiranga ", "Arctic red river ", "Puente Angostura ", &
153	! "Ceatal Izmail ", "Kinshasa ", "Salekhard ", &
154	! "Bahadurabad ", "Igarka ", "Kusur ", &
155	! "Malanville ", "Datong " /)
156
157	CONTAINS
158
159
160	SUBROUTINE routing_flow_get(slow_reservoir_r, fast_reservoir_r, stream_reservoir_r, riverflow_mean, coastalflow_mean, &
161	lakeinflow_mean)
162	IMPLICIT NONE
163	REAL(r_std),OPTIONAL, INTENT(OUT) :: slow_reservoir_r(:)
164	REAL(r_std),OPTIONAL, INTENT(OUT) :: fast_reservoir_r(:)
165	REAL(r_std),OPTIONAL, INTENT(OUT) :: stream_reservoir_r(:)
166	REAL(r_std),OPTIONAL, INTENT(OUT) :: riverflow_mean(:)
167	REAL(r_std),OPTIONAL, INTENT(OUT) :: coastalflow_mean(:)
168	REAL(r_std),OPTIONAL, INTENT(OUT) :: lakeinflow_mean(:)
169
170	CALL routing_flow_get_(slow_reservoir_r, fast_reservoir_r, stream_reservoir_r, riverflow_mean, coastalflow_mean, lakeinflow_mean)
171	END SUBROUTINE routing_flow_get
172
173	SUBROUTINE routing_flow_get_(slow_reservoir_r_, fast_reservoir_r_, stream_reservoir_r_, riverflow_mean_, &
174	coastalflow_mean_, lakeinflow_mean_)
175	IMPLICIT NONE
176	REAL(r_std),OPTIONAL, INTENT(OUT) :: slow_reservoir_r_(:)
177	REAL(r_std),OPTIONAL, INTENT(OUT) :: fast_reservoir_r_(:)
178	REAL(r_std),OPTIONAL, INTENT(OUT) :: stream_reservoir_r_(:)
179	REAL(r_std),OPTIONAL, INTENT(OUT) :: riverflow_mean_(:)
180	REAL(r_std),OPTIONAL, INTENT(OUT) :: coastalflow_mean_(:)
181	REAL(r_std),OPTIONAL, INTENT(OUT) :: lakeinflow_mean_(:)
182
183	IF (PRESENT(slow_reservoir_r_)) slow_reservoir_r_ = slow_reservoir_r
184	IF (PRESENT(fast_reservoir_r_)) fast_reservoir_r_ = fast_reservoir_r
185	IF (PRESENT(stream_reservoir_r_)) stream_reservoir_r_ = stream_reservoir_r
186	IF (PRESENT(riverflow_mean_)) riverflow_mean_ = riverflow_mean
187	IF (PRESENT(coastalflow_mean_)) coastalflow_mean_ = coastalflow_mean
188	IF (PRESENT(lakeinflow_mean_)) lakeinflow_mean_ = lakeinflow_mean
189
190	END SUBROUTINE routing_flow_get_
191
192	SUBROUTINE routing_flow_set(slow_reservoir_r, fast_reservoir_r, stream_reservoir_r, riverflow_mean, coastalflow_mean, &
193	lakeinflow_mean)
194	IMPLICIT NONE
195	REAL(r_std),OPTIONAL, INTENT(IN) :: slow_reservoir_r(:)
196	REAL(r_std),OPTIONAL, INTENT(IN) :: fast_reservoir_r(:)
197	REAL(r_std),OPTIONAL, INTENT(IN) :: stream_reservoir_r(:)
198	REAL(r_std),OPTIONAL, INTENT(IN) :: riverflow_mean(:)
199	REAL(r_std),OPTIONAL, INTENT(IN) :: coastalflow_mean(:)
200	REAL(r_std),OPTIONAL, INTENT(IN) :: lakeinflow_mean(:)
201
202	CALL routing_flow_set_(slow_reservoir_r, fast_reservoir_r, stream_reservoir_r, riverflow_mean, coastalflow_mean, lakeinflow_mean)
203
204	END SUBROUTINE routing_flow_set
205
206	SUBROUTINE routing_flow_set_(slow_reservoir_r_, fast_reservoir_r_, stream_reservoir_r_, riverflow_mean_, &
207	coastalflow_mean_, lakeinflow_mean_)
208	IMPLICIT NONE
209	REAL(r_std),OPTIONAL, INTENT(IN) :: slow_reservoir_r_(:)
210	REAL(r_std),OPTIONAL, INTENT(IN) :: fast_reservoir_r_(:)
211	REAL(r_std),OPTIONAL, INTENT(IN) :: stream_reservoir_r_(:)
212	REAL(r_std),OPTIONAL, INTENT(IN) :: riverflow_mean_(:)
213	REAL(r_std),OPTIONAL, INTENT(IN) :: coastalflow_mean_(:)
214	REAL(r_std),OPTIONAL, INTENT(IN) :: lakeinflow_mean_(:)
215
216	IF (PRESENT(slow_reservoir_r_)) slow_reservoir_r = slow_reservoir_r_
217	IF (PRESENT(fast_reservoir_r_)) fast_reservoir_r = fast_reservoir_r_
218	IF (PRESENT(stream_reservoir_r_)) stream_reservoir_r = stream_reservoir_r_
219	IF (PRESENT(riverflow_mean_)) riverflow_mean = riverflow_mean_
220	IF (PRESENT(coastalflow_mean_)) coastalflow_mean = coastalflow_mean_
221	IF (PRESENT(lakeinflow_mean_)) lakeinflow_mean = lakeinflow_mean_
222
223	END SUBROUTINE routing_flow_set_
224	!! =============================================================================================================================
225	!! SUBROUTINE: routing_simple_xios_initialize
226	!!
227	!>\BRIEF Initialize xios dependant defintion before closing context defintion
228	!!
229	!! DESCRIPTION: Initialize xios dependant defintion before closing context defintion.
230	!! This subroutine is called before the xios context is closed.
231	!!
232	!! RECENT CHANGE(S): None
233	!!
234	!! REFERENCE(S): None
235	!!
236	!! FLOWCHART: None
237	!! \n
238	!_ ==============================================================================================================================
239
240	SUBROUTINE routing_flow_xios_initialize
241	USE xios
242	USE routing, ONLY : routing_names
243	IMPLICIT NONE
244
245	INTEGER(i_std) ::ib
246	INTEGER :: nbasmax=1
247	!! 0 Variable and parameter description
248	CHARACTER(LEN=60),ALLOCATABLE :: label(:)
249	LOGICAL :: file_exists
250
251	IF (is_omp_root) THEN
252	CALL xios_get_axis_attr("basins", n_glo=basins_out) ! get nb basins to output
253	ALLOCATE(label(basins_out))
254	CALL routing_names(basins_out,label)
255	CALL xios_set_axis_attr("basins", label=label) ! set riverflow basins name
256	INQUIRE(FILE="routing_start.nc", EXIST=file_exists)
257	IF (file_exists) CALL xios_set_file_attr("routing_start", enabled=.TRUE.)
258	ENDIF
259
260	!! Define XIOS axis size needed for the model output
261	! Add axis for homogeneity between all routing schemes, these dimensions are currently not used in this scheme
262	CALL xios_orchidee_addaxis("nbhtu", nbasmax, (/(REAL(ib,r_std),ib=1,nbasmax)/))
263	CALL xios_orchidee_addaxis("nbasmon", 1, (/(REAL(ib,r_std),ib=1,1)/))
264
265	END SUBROUTINE routing_flow_xios_initialize
266
267
268	SUBROUTINE routing_flow_initialize(kjit, rest_id, nbpt_, dt_routing, contfrac, nbpt_r_, riverflow, coastalflow)
269	USE grid, ONLY : area
270	USE xios
271	IMPLICIT NONE
272	INTEGER ,INTENT(IN) :: kjit
273	INTEGER ,INTENT(IN) :: rest_id
274	INTEGER, INTENT(IN) :: nbpt_ !! nb points orchidee grid
275	REAL(r_std), INTENT(IN) :: dt_routing
276	REAL(r_std), INTENT(IN) :: contfrac(nbpt_) !! fraction of land
277	INTEGER, INTENT(OUT) :: nbpt_r_ !! nb points routing native grid
278	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)
279	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)
280	INTEGER :: ier
281	CHARACTER(LEN=80) :: var_name !! To store variables names for I/O (unitless)
282
283	nbpt = nbpt_
284	CALL routing_flow_init_local(contfrac, nbpt_r_)
285	nbpt_r = nbpt_r_
286	CALL routing_flow_init_mean(kjit, rest_id)
287	CALL initialize_stations(dt_routing)
288	!
289	! Put into the restart file the fluxes so that they can be regenerated at restart.
290	!
291	ALLOCATE (lakeinflow_mean(nbpt), stat=ier)
292	IF (ier /= 0) CALL ipslerr_p(3,'routing_flow_initialize','Pb in allocate for lakeinflow_mean','','')
293	var_name = 'lakeinflow'
294	CALL ioconf_setatt_p('UNITS', 'Kg/dt')
295	CALL ioconf_setatt_p('LONG_NAME','Lake inflow')
296	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., lakeinflow_mean, "gather", nbp_glo, index_g)
297	CALL setvar_p (lakeinflow_mean, val_exp, 'NO_KEYWORD', zero)
298
299	ALLOCATE (riverflow_mean(nbpt), stat=ier)
300	IF (ier /= 0) CALL ipslerr_p(3,'routing_flow_initialize','Pb in allocate for riverflow_mean','','')
301	var_name = 'riverflow'
302	CALL ioconf_setatt_p('UNITS', 'Kg/dt')
303	CALL ioconf_setatt_p('LONG_NAME','River flux into the sea')
304	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., riverflow_mean, "gather", nbp_glo, index_g)
305	CALL setvar_p (riverflow_mean, val_exp, 'NO_KEYWORD', zero)
306
307	ALLOCATE (coastalflow_mean(nbpt), stat=ier)
308	IF (ier /= 0) CALL ipslerr_p(3,'routing_flow_initialize','Pb in allocate for coastalflow_mean','','')
309	var_name = 'coastalflow'
310	CALL ioconf_setatt_p('UNITS', 'Kg/dt')
311	CALL ioconf_setatt_p('LONG_NAME','Diffuse flux into the sea')
312	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., coastalflow_mean, "gather", nbp_glo, index_g)
313	CALL setvar_p (coastalflow_mean, val_exp, 'NO_KEYWORD', zero)
314
315	riverflow(:) = riverflow_mean(:)
316	coastalflow(:) = coastalflow_mean(:)
317
318	CALL routing_flow_initialize_diag(kjit, rest_id, contfrac)
319
320	END SUBROUTINE routing_flow_initialize
321
322	SUBROUTINE routing_flow_initialize_diag(kjit, rest_id, contfrac)
323	USE xios
324	USE grid, ONLY : area
325	IMPLICIT NONE
326	INTEGER ,INTENT(IN) :: kjit
327	INTEGER ,INTENT(IN) :: rest_id
328	REAL(r_std), INTENT(IN) :: contfrac(nbpt) !! fraction of land
329	INTEGER :: ier
330	CHARACTER(LEN=80) :: var_name !! To store variables names for I/O (unitless)
331
332
333	REAL(r_std) :: contfrac_mpi(nbp_mpi)
334	REAL(r_std) :: area_mpi(nbp_mpi)
335	REAL(r_std) :: fast_diag_mpi(nbp_mpi)
336	REAL(r_std) :: slow_diag_mpi(nbp_mpi)
337	REAL(r_std) :: stream_diag_mpi(nbp_mpi)
338	REAL(r_std) :: hydrographs_diag_mpi(nbp_mpi)
339
340	CALL gather_omp(contfrac, contfrac_mpi)
341	CALL gather_omp(area, area_mpi)
342
343	IF (is_omp_root) THEN
344	ALLOCATE(fast_diag(nbpt), stat=ier)
345	IF (ier /= 0) CALL ipslerr_p(3,'routing_flow_initialize_diag','Pb in allocate for fast_diag','','')
346	ALLOCATE(slow_diag(nbpt), stat=ier)
347	IF (ier /= 0) CALL ipslerr_p(3,'routing_flow_initialize_diag','Pb in allocate for slow_diag','','')
348	ALLOCATE(stream_diag(nbpt), stat=ier)
349	IF (ier /= 0) CALL ipslerr_p(3,'routing_flow_initialize_diag','Pb in allocate for stream_diag','','')
350
351
352	CALL xios_send_field("routing_fast_diag0_r", fast_reservoir_r)
353	CALL xios_recv_field("routing_fast_diag0", fast_diag_mpi)
354	CALL xios_send_field("routing_slow_diag0_r", slow_reservoir_r)
355	CALL xios_recv_field("routing_slow_diag0", slow_diag_mpi)
356	CALL xios_send_field("routing_stream_diag0_r", stream_reservoir_r)
357	CALL xios_recv_field("routing_stream_diag0", stream_diag_mpi)
358
359	fast_diag_mpi=fast_diag_mpi/(area_mpi*contfrac_mpi) !! kg => kg/m^2
360	slow_diag_mpi=slow_diag_mpi/(area_mpi*contfrac_mpi) !! kg => kg/m^2
361	stream_diag_mpi=stream_diag_mpi/(area_mpi*contfrac_mpi) !! kg => kg/m^2
362	ENDIF
363
364	CALL scatter_omp(fast_diag_mpi,fast_diag)
365	CALL scatter_omp(slow_diag_mpi,slow_diag)
366	CALL scatter_omp(stream_diag_mpi,stream_diag)
367
368	ALLOCATE(hydrographs_diag(nbpt), stat=ier)
369	IF (ier /= 0) CALL ipslerr_p(3,'routing_flow_initialize_diag','Pb in allocate for hydrographs_diag','','')
370
371	var_name = 'hydrographs'
372	CALL ioconf_setatt_p('UNITS', 'kg/dt_sechiba')
373	CALL ioconf_setatt_p('LONG_NAME','Hydrograph at outlow of grid')
374	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., hydrographs_diag, "gather", nbp_glo, index_g)
375	CALL setvar_p (hydrographs_diag, val_exp, 'NO_KEYWORD', zero)
376
377	END SUBROUTINE routing_flow_initialize_diag
378
379
380	SUBROUTINE compute_coastline(contfrac, coastline)
381	USE mod_orchidee_para
382	USE grid
383	IMPLICIT NONE
384	REAL(r_std),INTENT(IN) :: contfrac(nbpt) ! INPUT : fraction of continent (unitless)
385	LOGICAL,INTENT(OUT) :: coastline(nbpt) ! OUTPUT : coastline mask (true : on coastaline, else false
386
387	REAL(r_std) :: contfrac_glo(nbp_glo)
388	REAL(r_std) :: contfrac2D_glo(iim_g*jjm_g)
389	LOGICAL :: coastline_glo(nbp_glo)
390	LOGICAL :: coastline2D_glo(iim_g*jjm_g)
391	LOGICAL :: mask1d(nbpt)
392	LOGICAL :: mask1d_glo(nbp_glo)
393	LOGICAL :: mask2d_glo(iim_g*jjm_g)
394	INTEGER(i_std) :: i,j,ij,next_i,next_j,next_ij,m,k
395	REAL(r_std), PARAMETER :: epsilon=1e-5
396	LOGICAL :: is_periodic
397
398	is_periodic = global
399	mask1d=.TRUE. ;
400	CALL gather(mask1d,mask1d_glo)
401	CALL gather(contfrac,contfrac_glo)
402
403
404	IF (is_mpi_root .AND. is_omp_root) THEN
405	contfrac2D_glo(:)=0
406	mask2d_glo=.FALSE.
407	DO i=1,nbp_glo
408	contfrac2D_glo(index_g(i)) = contfrac_glo(i)
409	mask2d_glo(index_g(i))=mask1d_glo(i)
410	ENDDO
411
412	coastline2D_glo(:)=.FALSE.
413	DO j=1,jjm_g
414	DO i=1,iim_g
415	ij=(j-1)*iim_g+i
416	IF (grid_type==unstructured) THEN
417	IF (contfrac2D_glo(ij)>epsilon .AND. contfrac2D_glo(ij)<1-epsilon) THEN
418	coastline2D_glo(ij)=.TRUE.
419	ENDIF
420	ELSE
421	DO k=-1,1
422	DO m=-1,1
423	IF (k==0 .AND. m==0) CYCLE
424	next_i=i+k
425	next_j=j+m
426
427	IF (next_i==0) THEN ! manage periodicity
428	IF (is_periodic) THEN
429	next_i=iim_g
430	ELSE
431	coastline2D_glo(ij)=.TRUE.
432	CYCLE
433	ENDIF
434	ENDIF
435
436	IF (next_i==iim_g+1) THEN! manage periodicity
437	IF (is_periodic) THEN
438	next_i=1 !! periodic but not true for limited area
439	ELSE
440	coastline2D_glo(ij)=.TRUE.
441	CYCLE
442	ENDIF
443	ENDIF
444	IF (next_j==0 .OR. next_j==jjm_g+1) THEN
445	IF (.NOT. is_periodic) coastline2D_glo(ij)=.TRUE.
446	CYCLE
447	ENDIF
448	next_ij = (next_j-1)*iim_g+next_i
449	IF (.NOT. mask2d_glo(next_ij)) coastline2D_glo(ij)=.TRUE.
450	! IF (contfrac2D_glo(next_ij)<=epsilon) coastline2D_glo(ij)=.TRUE.
451	ENDDO
452	ENDDO
453	ENDIF
454	ENDDO
455	ENDDO
456
457	DO i=1,nbp_glo
458	coastline_glo(i) = coastline2D_glo(index_g(i))
459	ENDDO
460
461	ENDIF
462
463	CALL scatter(coastline_glo, coastline)
464
465	END SUBROUTINE compute_coastline
466
467
468	SUBROUTINE new_routing_flow_correct_riverflow(ni, nj, contfrac, coastline, trip_r, trip_extended_r, topoind_r)
469	USE xios
470	! USE grid, ONLY : global
471	IMPLICIT NONE
472	INCLUDE "mpif.h"
473	INTEGER, INTENT(IN) :: ni ! INPUT : size i (longitude) of the local routing domain
474	INTEGER, INTENT(IN) :: nj ! INPUT : size i (latitude) of the local routing domain
475	REAL(r_std), INTENT(IN) :: contfrac(nbp_mpi) ! INPUT : continental fraction (unittless)
476	LOGICAL, INTENT(IN) :: coastline(nbp_mpi) ! INPUT/OUTPUT : coastline mask
477	REAL(r_std), INTENT(INOUT) :: trip_r(ni,nj) ! INPUT/OUTPUT : diection of flow, which will be modified by the routine
478	REAL(r_std), INTENT(INOUT) :: trip_extended_r(ni,nj) ! INPUT : direction of flow extended to ocean points
479	REAL(r_std), INTENT(INOUT) :: topoind_r(ni,nj) ! INPUT/OUTPUT : topographic index which will be modified by the routine if extended to ocean
480
481	REAL(r_std) :: mask(nbp_mpi)
482	REAL(r_std) :: contfrac_r(ni,nj)
483	REAL(r_std) :: mask_coast(nbp_mpi)
484	REAL(r_std) :: frac_coast_r(ni,nj)
485	REAL(r_std) :: trip_extended_rp1(0:ni+1,0:nj+1)
486	REAL(r_std) :: trip_rp1(0:ni+1,0:nj+1)
487	REAL(r_std) :: state_r(ni,nj)
488	REAL(r_std) :: state_rp1(0:ni+1,0:nj+1)
489	REAL(r_std) :: done_rp1(0:ni+1,0:nj+1)
490	REAL(r_std) :: trip_out_r(0:ni+1,0:nj+1)
491
492	INTEGER, PARAMETER :: is_ter=0
493	INTEGER, PARAMETER :: is_coast=1
494	INTEGER, PARAMETER :: is_oce=2
495	INTEGER, PARAMETER :: riverflow=99
496	INTEGER, PARAMETER :: coastalflow=98
497	INTEGER, PARAMETER :: lakeinflow=97
498	INTEGER, PARAMETER :: norouting=100
499
500	INTEGER :: updated
501	INTEGER :: it, i,j,k,m, nexti,nextj,previ,prevj,ierr
502	INTEGER :: ibegin, jbegin, ni_glo, nj_glo
503	REAL(r_std) :: lonvalue_1d(ni), latvalue_1d(nj)
504	REAL(r_std),PARAMETER :: epsilon=1e-5
505
506	TYPE(xios_duration) :: ts
507	TYPE(xios_domain) :: domain_hdl
508	TYPE(xios_domaingroup) :: domain_def_hdl
509	TYPE(xios_field) :: field_hdl
510	TYPE(xios_fieldgroup) :: field_def_hdl
511	TYPE(xios_file) :: file_hdl
512	TYPE(xios_filegroup) :: file_def_hdl
513	TYPE(xios_expand_domain) :: domain_expand_hdl
514	TYPE(xios_date) :: start_date, time_origin
515	CHARACTER(LEN=20) :: calendar_type
516	LOGICAL :: true=.TRUE.
517	LOGICAL :: global=.FALSE.
518	LOGICAL :: ok
519
520	CALL xios_get_domain_attr("routing_domain", ibegin=ibegin, jbegin=jbegin, ni_glo=ni_glo,nj_glo=nj_glo) ! get routing domain dimension
521	CALL xios_get_domain_attr("routing_domain", lonvalue_1d=lonvalue_1d, latvalue_1d=latvalue_1d) ! get routing domain dimension
522
523	! manage non periodic boundary in case of LAM
524	IF (.NOT. global) THEN
525	IF (ibegin==0) THEN
526	i=1 ;
527	DO j=1,nj
528	IF (trip_r(i,j)>=1 .AND. trip_r(i,j)<=8) trip_r(i,j) = coastalflow
529	IF (trip_extended_r(i,j)>=1 .AND. trip_extended_r(i,j)<=8) trip_extended_r(i,j) = coastalflow
530	ENDDO
531	ENDIF
532
533	IF (ibegin+ni==ni_glo) THEN
534	i=ni ;
535	DO j=1,nj
536	IF (trip_r(i,j)>=1 .AND. trip_r(i,j)<=8) trip_r(i,j) = coastalflow
537	IF (trip_extended_r(i,j)>=1 .AND. trip_extended_r(i,j)<=8) trip_extended_r(i,j) = coastalflow
538	ENDDO
539	ENDIF
540
541	IF (jbegin==0) THEN
542	j=1 ;
543	DO i=1,ni
544	IF (trip_r(i,j)>=1 .AND. trip_r(i,j)<=8) trip_r(i,j) = coastalflow
545	IF (trip_extended_r(i,j)>=1 .AND. trip_extended_r(i,j)<=8) trip_extended_r(i,j) = coastalflow
546	ENDDO
547	ENDIF
548
549	IF (jbegin+nj==nj_glo) THEN
550	j=nj ;
551	DO i=1,ni
552	IF (trip_r(i,j)>=1 .AND. trip_r(i,j)<=8) trip_r(i,j) = coastalflow
553	IF (trip_extended_r(i,j)>=1 .AND. trip_extended_r(i,j)<=8) trip_extended_r(i,j) = coastalflow
554	ENDDO
555	ENDIF
556	ENDIF
557
558	mask = 1
559	CALL xios_send_field("routing_contfrac", mask)
560	CALL xios_recv_field("routing_contfrac_r", contfrac_r)
561
562	CALL xios_send_field("trip_ext_r",trip_extended_r)
563	CALL xios_recv_field("trip_ext_rp1",trip_extended_rp1)
564
565	mask_coast=0
566	WHERE(coastline) mask_coast=1
567	CALL xios_send_field("mask_coastline",mask_coast)
568	CALL xios_recv_field("frac_coastline_r",frac_coast_r)
569
570	DO j=1,nj
571	DO i=1,ni
572	IF (frac_coast_r(i,j)>epsilon) THEN
573	state_r(i,j)=is_coast
574	ELSE
575	IF (contfrac_r(i,j)> epsilon) THEN
576	state_r(i,j)=is_ter
577	ELSE
578	state_r(i,j)=is_oce
579	ENDIF
580	ENDIF
581	ENDDO
582	ENDDO
583
584
585	DO j=1,nj
586	DO i=1,ni
587	IF (trip_r(i,j)>=100) THEN
588	trip_r(i,j)=norouting
589	ENDIF
590	ENDDO
591	ENDDO
592
593	CALL xios_context_initialize("orchidee_init_routing",MPI_COMM_ORCH)
594	CALL xios_orchidee_change_context("orchidee_init_routing")
595	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
596	CALL xios_set_domain_attr("routing_domain", lonvalue_1d=lonvalue_1d, latvalue_1d=latvalue_1d) ! get routing domain dimension
597	CALL xios_close_context_definition()
598	CALL xios_update_calendar(1)
599
600	CALL xios_send_field("trip_r_init",trip_r)
601
602	CALL xios_send_field("trip_r",trip_r)
603	CALL xios_recv_field("trip_rp1",trip_rp1)
604
605	CALL xios_send_field("state_r",state_r)
606	CALL xios_recv_field("state_rp1",state_rp1)
607
608	! DO j=1,nj
609	! DO i=1,ni
610	! IF (trip_rp1(i,j)<=8) THEN
611	! CALL next_trip(trip_rp1, ni, nj, i, j, nexti, nextj)
612	! IF (trip_rp1(nexti,nextj)>=100) trip_rp1(i,j)=98 ! unrouted point becomes coastal flow
613	! ENDIF
614	! ENDDO
615	! ENDDO
616
617	state_r(1:ni,1:nj) = state_rp1(1:ni,1:nj)
618	trip_r(1:ni,1:nj) = trip_rp1(1:ni,1:nj)
619
620	updated=1
621	it=2
622	DO WHILE(updated>0)
623	updated=0
624
625	CALL xios_update_calendar(it)
626
627	CALL xios_send_field("trip_r", trip_r)
628	CALL xios_recv_field("trip_rp1",trip_rp1)
629
630	CALL xios_send_field("state_r",state_r)
631	CALL xios_recv_field("state_rp1",state_rp1)
632
633	state_r(1:ni,1:nj) = state_rp1(1:ni,1:nj)
634	trip_r(1:ni,1:nj) = trip_rp1(1:ni,1:nj)
635
636	! first => riverflow and coastalflow routing to ocean must be routed to coast => to long
637	DO j=1,nj
638	DO i=1,ni
639	CALL next_trip(trip_rp1, ni, nj, i, j, nexti, nextj, ok)
640	IF (ok) THEN
641	IF (state_rp1(nexti,nextj)==is_oce .AND. (trip_rp1(nexti,nextj)==riverflow .OR. trip_rp1(nexti,nextj)==coastalflow .OR. trip_rp1(nexti,nextj)==lakeinflow)) THEN
642	trip_r(i,j) = trip_rp1(nexti,nextj)
643	updated=updated + 1
644	ENDIF
645	ENDIF
646	ENDDO
647	ENDDO
648
649	DO j=1,nj
650	DO i=1,ni
651	IF (state_rp1(i,j)==is_oce) THEN
652	IF (trip_rp1(i,j)==riverflow .OR. trip_rp1(i,j)==coastalflow .OR. trip_rp1(i,j)==lakeinflow) THEN
653	trip_r(i,j) = norouting
654	updated=updated+1
655	ENDIF
656	ENDIF
657	ENDDO
658	ENDDO
659
660	CALL MPI_ALLREDUCE(MPI_IN_PLACE , updated, 1, MPI_INT_ORCH, MPI_SUM, MPI_COMM_ORCH, ierr)
661	it = it +1
662	ENDDO
663
664	! second prolongate riverflow and coastalflow to coast if they are in land
665	updated=1
666	DO WHILE(updated>0)
667	updated=0
668
669	CALL xios_update_calendar(it)
670
671	CALL xios_send_field("trip_r", trip_r)
672	CALL xios_recv_field("trip_rp1",trip_rp1)
673
674	CALL xios_send_field("state_r",state_r)
675	CALL xios_recv_field("state_rp1",state_rp1)
676
677	state_r(1:ni,1:nj) = state_rp1(1:ni,1:nj)
678	trip_r(1:ni,1:nj) = trip_rp1(1:ni,1:nj)
679
680	DO j=1,nj
681	DO i=1,ni
682	IF (state_rp1(i,j)==is_ter .AND. (trip_rp1(i,j)==riverflow .OR. trip_rp1(i,j)==coastalflow .OR. trip_rp1(i,j)==lakeinflow)) THEN
683	CALL next_trip(trip_extended_rp1, ni, nj, i, j, nexti, nextj, ok)
684	IF (ok) THEN
685	IF (state_rp1(nexti,nextj)==is_ter .OR. state_rp1(nexti,nextj)==is_coast) THEN
686	trip_r(i,j) = trip_extended_rp1(i,j)
687	topoind_r(i,j) = 1e-10 ! flow instantanesly for now but better to take the topind of the incomming flux
688	updated=updated+1
689	ENDIF
690	ENDIF
691	ENDIF
692	ENDDO
693	ENDDO
694
695	DO j=1,nj
696	DO i=1,ni
697	IF (state_rp1(i,j)==is_ter .OR. state_rp1(i,j)==is_coast) THEN
698	DO k=-1,1
699	DO m=-1,1
700	IF (k==0 .AND. m==0) CYCLE
701	prevj = j+k ; previ = i+m
702	CALL next_trip(trip_extended_rp1, ni, nj, previ, prevj, nexti, nextj, ok)
703	IF (ok) THEN
704	IF (nexti==i .AND. nextj==j) THEN
705	IF (state_rp1(previ,prevj)==is_ter .AND. (trip_rp1(previ,prevj)==riverflow .OR. trip_rp1(previ,prevj)==coastalflow .OR. trip_rp1(previ,prevj)==lakeinflow)) THEN
706	IF (trip_r(i,j)/=riverflow) THEN
707	trip_r(i,j)=trip_rp1(previ,prevj)
708	ELSE
709	trip_r(i,j)=riverflow
710	ENDIF
711	updated=updated+1
712	ENDIF
713	ENDIF
714	ENDIF
715	ENDDO
716	ENDDO
717	ENDIF
718	ENDDO
719	ENDDO
720
721	CALL MPI_ALLREDUCE(MPI_IN_PLACE , updated, 1, MPI_INT_ORCH, MPI_SUM, MPI_COMM_ORCH, ierr)
722	it = it +1
723	ENDDO
724
725	CALL xios_update_calendar(it)
726
727	CALL xios_send_field("trip_r", trip_r)
728	CALL xios_recv_field("trip_rp1",trip_rp1)
729
730	CALL xios_send_field("state_r",state_r)
731	CALL xios_recv_field("state_rp1",state_rp1)
732	DO j=1,nj
733	DO i=1,ni
734	IF (state_rp1(i,j)==is_coast) THEN
735	IF (trip_r(i,j)==norouting) trip_r(i,j) = coastalflow
736	ENDIF
737	ENDDO
738	ENDDO
739
740	! it = it +1
741	! CALL xios_update_calendar(it)
742	!
743	! CALL xios_send_field("trip_r", trip_r)
744	! CALL xios_recv_field("trip_rp1",trip_rp1)
745	!
746	! CALL xios_send_field("state_r",state_r)
747	! CALL xios_recv_field("state_rp1",state_rp1)
748	!
749	! DO j=1,nj
750	! DO i=1,ni
751	! IF (state_rp1(i,j)==is_oce) THEN
752	! trip_r(i,j) = norouting
753	! ENDIF
754	! ENDDO
755	! ENDDO
756
757	it = it +1
758	CALL xios_update_calendar(it)
759
760	CALL xios_send_field("trip_r", trip_r)
761	CALL xios_recv_field("trip_rp1",trip_rp1)
762
763	CALL xios_send_field("state_r",state_r)
764	CALL xios_recv_field("state_rp1",state_rp1)
765	! check
766	updated=0
767	DO j=1,nj
768	DO i=1,ni
769	IF (state_rp1(i,j)==is_ter) THEN
770	IF (trip_r(i,j) == coastalflow .OR. trip_r(i,j) == riverflow .OR. trip_r(i,j) == norouting) THEN
771	updated=updated+1
772	trip_r(i,j) = lakeinflow
773	ENDIF
774	ELSE IF (state_rp1(i,j)==is_coast) THEN
775	IF (trip_r(i,j) == norouting) THEN
776	updated=updated+1
777	trip_r(i,j) = coastalflow
778	ENDIF
779	ENDIF
780	ENDDO
781	ENDDO
782
783
784	CALL MPI_ALLREDUCE(MPI_IN_PLACE , updated, 1, MPI_INT_ORCH, MPI_SUM, MPI_COMM_ORCH, ierr)
785	IF (updated/=0 .AND. is_mpi_root) PRINT*,"WARNING : ",updated," riverflow or coastal flow are not on the coast, ie in middle of land ==> converted into lakeinflow"
786	CALL MPI_BARRIER(MPI_COMM_ORCH, ierr)
787
788	updated=0
789	DO j=1,nj
790	DO i=1,ni
791	CALL next_trip(trip_rp1,ni,nj,i,j,nexti,nextj,ok)
792	IF (ok) THEN
793	IF (trip_rp1(nexti,nextj)==norouting) THEN
794	updated=updated+1
795	PRINT*,"point",ibegin+i,jbegin+j," (trip(i,j)=",trip_rp1(i,j),") is routed to nothing => ",ibegin+nexti,jbegin+nextj, &
796	" (trip(i,j)=",trip_rp1(nexti,nextj),")"
797	ENDIF
798	ENDIF
799	ENDDO
800	ENDDO
801
802	CALL MPI_BARRIER(MPI_COMM_ORCH, ierr)
803	CALL MPI_ALLREDUCE(MPI_IN_PLACE , updated, 1, MPI_INT_ORCH, MPI_SUM, MPI_COMM_ORCH, ierr)
804	IF (updated/=0) THEN
805	IF (is_mpi_root) PRINT*,"ERROR : ",updated," points are not correctly routed (routed to nothing)"
806	CALL xios_context_finalize
807	STOP
808	ENDIF
809
810
811	updated=0
812	DO j=1,nj
813	DO i=1,ni
814	CALL next_trip(trip_rp1,ni,nj,i,j,nexti,nextj,ok)
815	IF (ok) THEN
816	IF (state_rp1(nexti,nextj)==is_oce .AND. (trip_rp1(nexti,nextj)==riverflow .OR. trip_rp1(nexti,nextj)==coastalflow .OR. trip_rp1(nexti,nextj)==lakeinflow )) THEN
817	updated=updated+1
818	PRINT*,"point",ibegin+i,jbegin+j," (trip(i,j)=",trip_rp1(i,j),") is routed to coastalflow/riverflow/lakeinflow point that is not on land grid => ",ibegin+nexti,jbegin+nextj, &
819	" (trip(i,j)=",trip_rp1(nexti,nextj),")"
820	ENDIF
821	ENDIF
822	ENDDO
823	ENDDO
824
825	CALL MPI_BARRIER(MPI_COMM_ORCH, ierr)
826	CALL MPI_ALLREDUCE(MPI_IN_PLACE , updated, 1, MPI_INT_ORCH, MPI_SUM, MPI_COMM_ORCH, ierr)
827	IF (updated/=0) THEN
828	IF (is_mpi_root) PRINT*,"ERROR : ",updated," points are routed to coastalflow/riverflow/lakeinflow points that are not on land grid"
829	CALL xios_context_finalize
830	STOP
831	ENDIF
832
833
834
835	it = it +1
836	CALL xios_update_calendar(it)
837
838	CALL xios_send_field("trip_r", trip_r)
839	CALL xios_recv_field("trip_rp1",trip_rp1)
840
841	CALL xios_send_field("state_r",state_r)
842	CALL xios_recv_field("state_rp1",state_rp1)
843	state_r(1:ni,1:nj) = state_rp1(1:ni,1:nj)
844	trip_r(1:ni,1:nj) = trip_rp1(1:ni,1:nj)
845
846	CALL xios_send_field("trip_r_final",trip_r)
847
848	CALL xios_context_finalize
849	CALL xios_orchidee_change_context("orchidee")
850
851	CONTAINS
852
853	SUBROUTINE next_trip(trip, ni, nj, i, j, nexti, nextj, ok)
854	IMPLICIT NONE
855	REAL, INTENT(IN) :: trip(0:ni+1,0:nj+1)
856	INTEGER, INTENT(IN) :: ni
857	INTEGER, INTENT(IN) :: nj
858	INTEGER, INTENT(IN) :: i
859	INTEGER, INTENT(IN) :: j
860	INTEGER, INTENT(OUT) :: nexti
861	INTEGER, INTENT(OUT) :: nextj
862	LOGICAL, INTENT(OUT) :: ok
863
864	ok=.TRUE.
865
866	SELECT CASE(NINT(trip(i,j)))
867	CASE (1)
868	nextj=j-1 ; nexti=i ! north
869	CASE (2)
870	nextj=j-1 ; nexti=i+1 ! north-east
871	CASE (3)
872	nextj=j ; nexti=i+1 ! east
873	CASE (4)
874	nextj=j+1 ; nexti=i+1 ! south-east
875	CASE (5)
876	nextj=j+1 ; nexti=i ! south
877	CASE(6)
878	nextj=j+1 ; nexti=i-1 ! south-west
879	CASE (7)
880	nextj=j ; nexti=i-1 ! west
881	CASE (8)
882	nextj=j-1 ; nexti=i-1 ! north-west
883	CASE DEFAULT
884	ok=.FALSE.
885	END SELECT
886
887	END SUBROUTINE next_trip
888
889	END SUBROUTINE new_routing_flow_correct_riverflow
890
891
892
893
894	SUBROUTINE routing_flow_correct_riverflow(ni, nj, contfrac, coastline, trip_r, trip_extended_r, topoind_r)
895	USE xios
896	! USE grid, ONLY : global
897	IMPLICIT NONE
898	INCLUDE "mpif.h"
899	INTEGER, INTENT(IN) :: ni ! INPUT : size i (longitude) of the local routing domain
900	INTEGER, INTENT(IN) :: nj ! INPUT : size i (latitude) of the local routing domain
901	REAL(r_std), INTENT(IN) :: contfrac(nbp_mpi) ! INPUT : continental fraction (unittless)
902	LOGICAL, INTENT(IN) :: coastline(nbp_mpi) ! INPUT/OUTPUT : coastline mask
903	REAL(r_std), INTENT(INOUT) :: trip_r(ni,nj) ! INPUT/OUTPUT : diection of flow, which will be modified by the routine
904	REAL(r_std), INTENT(INOUT) :: trip_extended_r(ni,nj) ! INPUT : direction of flow extended to ocean points
905	REAL(r_std), INTENT(INOUT) :: topoind_r(ni,nj) ! INPUT/OUTPUT : topographic index which will be modified by the routine if extended to ocean
906
907	REAL(r_std) :: contfrac_r(ni,nj)
908	REAL(r_std) :: mask_coast(nbp_mpi)
909	REAL(r_std) :: frac_coast_r(ni,nj)
910	REAL(r_std) :: trip_extended_rp1(0:ni+1,0:nj+1)
911	REAL(r_std) :: trip_rp1(0:ni+1,0:nj+1)
912	REAL(r_std) :: state_r(ni,nj)
913	REAL(r_std) :: state_rp1(0:ni+1,0:nj+1)
914
915	INTEGER, PARAMETER :: is_ter=0
916	INTEGER, PARAMETER :: is_coast=1
917	INTEGER, PARAMETER :: is_oce=2
918	INTEGER :: updated
919	INTEGER :: it, i,j,k,m, nexti,nextj,previ,prevj,ierr
920	INTEGER :: ibegin, jbegin, ni_glo, nj_glo
921	REAL(r_std) :: lonvalue_1d(ni), latvalue_1d(nj)
922	REAL(r_std),PARAMETER :: epsilon=1e-5
923
924	TYPE(xios_duration) :: ts
925	TYPE(xios_domain) :: domain_hdl
926	TYPE(xios_domaingroup) :: domain_def_hdl
927	TYPE(xios_field) :: field_hdl
928	TYPE(xios_fieldgroup) :: field_def_hdl
929	TYPE(xios_file) :: file_hdl
930	TYPE(xios_filegroup) :: file_def_hdl
931	TYPE(xios_expand_domain) :: domain_expand_hdl
932	TYPE(xios_date) :: start_date, time_origin
933	CHARACTER(LEN=20) :: calendar_type
934	LOGICAL :: true=.TRUE.
935	LOGICAL :: global=.FALSE.
936
937	CALL xios_get_domain_attr("routing_domain", ibegin=ibegin, jbegin=jbegin, ni_glo=ni_glo,nj_glo=nj_glo) ! get routing domain dimension
938	CALL xios_get_domain_attr("routing_domain", lonvalue_1d=lonvalue_1d, latvalue_1d=latvalue_1d) ! get routing domain dimension
939
940	! manage non periodic boundary in case of LAM
941	IF (.NOT. global) THEN
942	IF (ibegin==0) THEN
943	i=1 ;
944	DO j=1,nj
945	IF (trip_r(i,j)>=1 .AND. trip_r(i,j)<=8) trip_r(i,j) = 97
946	IF (trip_extended_r(i,j)>=1 .AND. trip_extended_r(i,j)<=8) trip_extended_r(i,j) = 97
947	ENDDO
948	ENDIF
949
950	IF (ibegin+ni==ni_glo) THEN
951	i=ni ;
952	DO j=1,nj
953	IF (trip_r(i,j)>=1 .AND. trip_r(i,j)<=8) trip_r(i,j) = 97
954	IF (trip_extended_r(i,j)>=1 .AND. trip_extended_r(i,j)<=8) trip_extended_r(i,j) = 97
955	ENDDO
956	ENDIF
957
958	IF (jbegin==0) THEN
959	j=1 ;
960	DO i=1,ni
961	IF (trip_r(i,j)>=1 .AND. trip_r(i,j)<=8) trip_r(i,j) = 97
962	IF (trip_extended_r(i,j)>=1 .AND. trip_extended_r(i,j)<=8) trip_extended_r(i,j) = 97
963	ENDDO
964	ENDIF
965
966	IF (jbegin+nj==nj_glo) THEN
967	j=nj ;
968	DO i=1,ni
969	IF (trip_r(i,j)>=1 .AND. trip_r(i,j)<=8) trip_r(i,j) = 97
970	IF (trip_extended_r(i,j)>=1 .AND. trip_extended_r(i,j)<=8) trip_extended_r(i,j) = 97
971	ENDDO
972	ENDIF
973	ENDIF
974
975
976	CALL xios_send_field("routing_contfrac",contfrac)
977	CALL xios_recv_field("routing_contfrac_r", contfrac_r)
978
979	CALL xios_send_field("trip_ext_r",trip_extended_r)
980	CALL xios_recv_field("trip_ext_rp1",trip_extended_rp1)
981
982	mask_coast=0
983	WHERE(coastline) mask_coast=1
984	CALL xios_send_field("mask_coastline",mask_coast)
985	CALL xios_recv_field("frac_coastline_r",frac_coast_r)
986
987	DO j=1,nj
988	DO i=1,ni
989	IF (contfrac_r(i,j)<epsilon .OR. contfrac_r(i,j)> 2 ) THEN
990	state_r(i,j)=is_oce
991	ELSE IF (contfrac_r(i,j)>1-epsilon) THEN
992	state_r(i,j)=is_ter
993	ELSE
994	state_r(i,j)=is_coast
995	ENDIF
996
997	IF (frac_coast_r(i,j)>epsilon .AND. frac_coast_r(i,j)<=1+epsilon ) state_r(i,j)=is_coast
998	ENDDO
999	ENDDO
1000
1001
1002	DO j=1,nj
1003	DO i=1,ni
1004	IF (trip_r(i,j)>=100) THEN
1005	IF (state_r(i,j) /= is_oce) trip_r(i,j)=trip_extended_r(i,j)
1006	ENDIF
1007	ENDDO
1008	ENDDO
1009
1010	CALL xios_context_initialize("orchidee_init_routing",MPI_COMM_ORCH)
1011	CALL xios_orchidee_change_context("orchidee_init_routing")
1012	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
1013	CALL xios_set_domain_attr("routing_domain", lonvalue_1d=lonvalue_1d, latvalue_1d=latvalue_1d) ! get routing domain dimension
1014	CALL xios_close_context_definition()
1015	CALL xios_update_calendar(1)
1016
1017	CALL xios_send_field("trip_r_init",trip_r)
1018
1019	CALL xios_send_field("trip_r",trip_r)
1020	CALL xios_recv_field("trip_rp1",trip_rp1)
1021
1022	CALL xios_send_field("state_r",state_r)
1023	CALL xios_recv_field("state_rp1",state_rp1)
1024
1025	DO j=1,nj
1026	DO i=1,ni
1027	IF (trip_rp1(i,j)<=8) THEN
1028	CALL next_trip(trip_rp1, ni, nj, i, j, nexti, nextj)
1029	IF (trip_rp1(nexti,nextj)>=100) trip_rp1(i,j)=98 ! unrouted point becomes coastal flow
1030	ENDIF
1031	ENDDO
1032	ENDDO
1033
1034
1035	state_r(1:ni,1:nj) = state_rp1(1:ni,1:nj)
1036	trip_r(1:ni,1:nj) = trip_rp1(1:ni,1:nj)
1037
1038	updated=1
1039	it=2
1040	DO WHILE(updated>0 .AND. it<MAX(ni_glo,nj_glo))
1041	updated=0
1042
1043	CALL xios_update_calendar(it)
1044
1045	CALL xios_send_field("trip_r", trip_r)
1046	CALL xios_recv_field("trip_rp1",trip_rp1)
1047
1048	CALL xios_send_field("state_r",state_r)
1049	CALL xios_recv_field("state_rp1",state_rp1)
1050
1051	state_r(1:ni,1:nj) = state_rp1(1:ni,1:nj)
1052	trip_r(1:ni,1:nj) = trip_rp1(1:ni,1:nj)
1053
1054
1055
1056	DO j=0,nj+1
1057	DO i=0,ni+1
1058
1059	IF (trip_rp1(i,j)==97 .OR. trip_rp1(i,j)==98 .OR. trip_rp1(i,j)==99) THEN
1060
1061	IF (state_rp1(i,j)==is_coast) THEN ! ok
1062
1063	IF (trip_rp1(i,j)==97) THEN
1064	IF (i>=1 .AND. i<=ni .AND. j>=1 .AND. j<=nj) trip_r(i,j)=98 ! => lakinflow will becomes coastalflow
1065	ELSE
1066	! ok for riverflow
1067	ENDIF
1068
1069	ELSE IF (state_rp1(i,j) == is_ter) THEN ! too short
1070
1071	IF (trip_rp1(i,j)==97) THEN
1072	!ok for lakinflow
1073	ELSE ! => for riverflow/coastalflow
1074	IF (trip_extended_rp1(i,j)==97) THEN ! definitively a lake that cannot be routed further away
1075	IF (i>=1 .AND. i<=ni .AND. j>=1 .AND. j<=nj) THEN
1076	trip_r(i,j)=97
1077	updated=updated+1
1078	ENDIF
1079	ELSE
1080	CALL next_trip(trip_extended_rp1, ni, nj, i, j, nexti, nextj)
1081	IF (nexti/=-1 .AND. nextj/=-1) THEN
1082	IF (nexti==i .AND. nextj==j) THEN
1083	!! routed on itself in midle of land => lake
1084	IF (nexti>=1 .AND. nexti<=ni .AND. nextj>=1 .AND. nextj<=nj) trip_r(i,j)=97
1085	ELSE
1086
1087	IF (state_rp1(nexti,nextj)==is_coast .OR. state_rp1(nexti,nextj)==is_ter) THEN
1088
1089	IF (i>=1 .AND. i<=ni .AND. j>=1 .AND. j<=nj) THEN
1090	trip_r(i,j)=trip_extended_rp1(i,j)
1091	state_r(i,j)= is_ter
1092	topoind_r(i,j) = 1e-10 ! flow instantanesly
1093	updated=updated+1
1094	ENDIF
1095
1096	IF (nexti>=1 .AND. nexti<=ni .AND. nextj>=1 .AND. nextj<=nj) THEN
1097	IF (trip_r(nexti,nextj)/=trip_rp1(i,j)) THEN
1098	trip_r(nexti,nextj)=trip_rp1(i,j)
1099	updated=updated+1
1100	ENDIF
1101	ENDIF
1102
1103	ELSE
1104	IF (i>=1 .AND. i<=ni .AND. j>=1 .AND. j<=nj) THEN
1105	STOP 'not possible'
1106	ENDIF
1107	ENDIF
1108	ENDIF
1109	ENDIF
1110	ENDIF
1111	ENDIF
1112
1113	ELSE IF (state_rp1(i,j) == is_oce) THEN ! too long
1114
1115	DO k=-1,1
1116	DO m=-1,1
1117	IF (k==0 .AND. m==0) CYCLE
1118	prevj = j+k ; previ = i+m
1119	CALL next_trip(trip_rp1, ni, nj, previ, prevj, nexti, nextj)
1120	! 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)
1121	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)
1122	ENDDO
1123	ENDDO
1124
1125	IF (i>=1 .AND. i<=ni .AND. j>=1 .AND. j<=nj) THEN
1126	trip_r(i,j)=1000
1127	updated=updated+1
1128	ENDIF
1129	ENDIF
1130
1131	ENDIF
1132	ENDDO
1133	ENDDO
1134
1135	CALL MPI_ALLREDUCE(MPI_IN_PLACE , updated, 1, MPI_INT_ORCH, MPI_SUM, MPI_COMM_ORCH, ierr)
1136	it = it +1
1137	ENDDO
1138
1139	CALL xios_send_field("trip_r_final",trip_r)
1140
1141	CALL xios_context_finalize
1142	CALL xios_orchidee_change_context("orchidee")
1143
1144	CONTAINS
1145
1146	SUBROUTINE next_trip(trip, ni, nj, i, j, nexti, nextj)
1147	IMPLICIT NONE
1148	REAL, INTENT(IN) :: trip(0:ni+1,0:nj+1)
1149	INTEGER, INTENT(IN) :: ni
1150	INTEGER, INTENT(IN) :: nj
1151	INTEGER, INTENT(IN) :: i
1152	INTEGER, INTENT(IN) :: j
1153	INTEGER, INTENT(OUT) :: nexti
1154	INTEGER, INTENT(OUT) :: nextj
1155
1156	IF (i<0 .OR. i>ni+1 .OR. j<0 .OR. j>nj+1) THEN
1157	nexti=-1
1158	nextj=-1
1159	RETURN
1160	ENDIF
1161
1162	SELECT CASE(NINT(trip(i,j)))
1163	CASE (1)
1164	nextj=j-1 ; nexti=i ! north
1165	CASE (2)
1166	nextj=j-1 ; nexti=i+1 ! north-east
1167	CASE (3)
1168	nextj=j ; nexti=i+1 ! east
1169	CASE (4)
1170	nextj=j+1 ; nexti=i+1 ! south-east
1171	CASE (5)
1172	nextj=j+1 ; nexti=i ! south
1173	CASE(6)
1174	nextj=j+1 ; nexti=i-1 ! south-west
1175	CASE (7)
1176	nextj=j ; nexti=i-1 ! west
1177	CASE (8)
1178	nextj=j-1 ; nexti=i-1 ! north-west
1179	CASE DEFAULT
1180	nextj=-1 ; nexti=-1 ! undefined
1181	END SELECT
1182
1183	IF (nexti<0 .OR. nexti>ni+1) nexti=-1
1184	IF (nextj<0 .OR. nextj>nj+1) nextj=-1
1185
1186	END SUBROUTINE next_trip
1187
1188	END SUBROUTINE routing_flow_correct_riverflow
1189
1190
1191
1192
1193
1194
1195
1196	SUBROUTINE routing_flow_init_mean(kjit, rest_id)
1197	USE ioipsl
1198	IMPLICIT NONE
1199	INTEGER(i_std), INTENT(in) :: kjit
1200	INTEGER(i_std), INTENT(in) :: rest_id !! Restart file identifier (unitless)
1201
1202	INTEGER :: ier
1203	CHARACTER(LEN=80) :: var_name !! To store variables names for I/O (unitless)
1204
1205
1206	! Get from the restart the fluxes we accumulated.
1207	!
1208
1209	ALLOCATE (runoff_mean(nbpt), stat=ier)
1210	runoff_mean(:) = 0
1211	IF (ier /= 0) CALL ipslerr_p(3,'routing_flow_init_mean','Pb in allocate for runoff_mean','','')
1212	var_name = 'runoff_route'
1213	CALL ioconf_setatt_p('UNITS', 'Kg')
1214	CALL ioconf_setatt_p('LONG_NAME','Accumulated runoff for routing')
1215	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., runoff_mean, "gather", nbp_glo, index_g)
1216	CALL setvar_p (runoff_mean, val_exp, 'NO_KEYWORD', zero)
1217
1218	ALLOCATE(drainage_mean(nbpt), stat=ier)
1219	drainage_mean(:) = 0
1220	IF (ier /= 0) CALL ipslerr_p(3,'routing_flow_init_mean','Pb in allocate for drainage_mean','','')
1221	var_name = 'drainage_route'
1222	CALL ioconf_setatt_p('UNITS', 'Kg')
1223	CALL ioconf_setatt_p('LONG_NAME','Accumulated drainage for routing')
1224	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., drainage_mean, "gather", nbp_glo, index_g)
1225	CALL setvar_p (drainage_mean, val_exp, 'NO_KEYWORD', zero)
1226
1227	END SUBROUTINE routing_flow_init_mean
1228
1229	SUBROUTINE routing_flow_make_mean(runoff, drainage)
1230	IMPLICIT NONE
1231	REAL(r_std),INTENT(IN) :: runoff(:)
1232	REAL(r_std),INTENT(IN) :: drainage(:)
1233
1234	runoff_mean(:) = runoff_mean(:) + runoff
1235	drainage_mean(:) = drainage_mean(:) + drainage
1236
1237	END SUBROUTINE routing_flow_make_mean
1238
1239
1240	SUBROUTINE routing_flow_reset_mean
1241	IMPLICIT NONE
1242	runoff_mean(:) = 0
1243	drainage_mean(:) = 0
1244
1245	END SUBROUTINE routing_flow_reset_mean
1246
1247
1248	SUBROUTINE routing_flow_finalize_mean(kjit, rest_id)
1249	USE ioipsl
1250	IMPLICIT NONE
1251	INTEGER,INTENT(IN) :: kjit
1252	INTEGER,INTENT(IN) :: rest_id
1253
1254	CALL restput_p (rest_id, 'runoff_route', nbp_glo, 1, 1, kjit, runoff_mean, 'scatter', nbp_glo, index_g)
1255	CALL restput_p (rest_id, 'drainage_route', nbp_glo, 1, 1, kjit, drainage_mean, 'scatter', nbp_glo, index_g)
1256
1257	DEALLOCATE(runoff_mean)
1258	DEALLOCATE(drainage_mean)
1259
1260	END SUBROUTINE routing_flow_finalize_mean
1261
1262
1263	SUBROUTINE routing_flow_init_local(contfrac, nbpt_r)
1264	USE xios
1265	USE grid, ONLY : area
1266	IMPLICIT NONE
1267	INCLUDE "mpif.h"
1268
1269	!! 0 Variable and parameter description
1270	!! 0.1 Input variables
1271	REAL(r_std),INTENT(IN) :: contfrac(nbpt) !! fraction of land
1272	INTEGER,INTENT(OUT) :: nbpt_r !! nb points routing native grid
1273
1274	!! 0.2 Local variables
1275	INTEGER :: ni !! longitude dimension of local routing grid
1276	INTEGER :: nj !! latitude dimension of local routing grid
1277	REAL(r_std) :: contfrac_mpi(nbp_mpi)
1278	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)
1279	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)
1280	!! routing is artificially computed on sea and endoric basins
1281	REAL(r_std),ALLOCATABLE :: diag_r(:) !! fraction of routing cell intesected by coastal cells of native grid
1282	REAL(r_std),ALLOCATABLE :: frac_lake_r(:) !! fraction of routing cell intesected by cells of native grid
1283	REAL(r_std),ALLOCATABLE :: frac_coast_r(:)
1284
1285	REAL(r_std) :: diag(nbp_mpi) !! fraction of routing cell intesected by coastal cells of native grid
1286
1287	LOGICAL :: coastline(nbpt)
1288	LOGICAL :: coastline_mpi(nbp_mpi)
1289	REAL(r_std) :: area_mpi(nbp_mpi)
1290
1291	INTEGER :: ij, ij_r, ij_rp1, i,j,jp1,jm1,ip1,im1,jr,ir
1292	INTEGER :: basins_count_mpi
1293	INTEGER :: nb_coast_points
1294	INTEGER :: ierr
1295	LOGICAL :: file_exists
1296	REAL(r_std) :: epsilon = 1e-5
1297	CHARACTER(LEN=255) :: routing_file_type
1298	!_ ================================================================================================================================
1299	!
1300	!> A value for property of each reservoir (in day/m) is given to compute a time constant (in day)
1301	!> for each reservoir (product of tcst and topo_resid).
1302	!> The value of tcst has been calibrated for the three reservoirs over the Senegal river basin only,
1303	!> during the 1 degree NCEP Corrected by Cru (NCC) resolution simulations (Ngo-Duc et al., 2005, Ngo-Duc et al., 2006) and
1304	!> generalized for all the basins of the world. The "slow reservoir" and the "fast reservoir"
1305	!> have the highest value in order to simulate the groundwater.
1306	!> The "stream reservoir", which represents all the water of the stream, has the lowest value.
1307	!> Those figures are the same for all the basins of the world.
1308
1309	CALL getin_p('SLOW_TCST', slow_tcst)
1310	!
1311	!Config Key = FAST_TCST
1312	!Config Desc = Time constant for the fast reservoir
1313	!Config If = RIVER_ROUTING
1314	!Config Def = 3.0
1315	!Config Help = This parameters allows the user to fix the
1316	!Config time constant (in days) of the fast reservoir
1317	!Config in order to get better river flows for
1318	!Config particular regions.
1319	!Config Units = [days]
1320
1321	CALL getin_p('FAST_TCST', fast_tcst)
1322	!
1323	!Config Key = STREAM_TCST
1324	!Config Desc = Time constant for the stream reservoir
1325	!Config If = RIVER_ROUTING
1326	!Config Def = 0.24
1327	!Config Help = This parameters allows the user to fix the
1328	!Config time constant (in days) of the stream reservoir
1329	!Config in order to get better river flows for
1330	!Config particular regions.
1331	!Config Units = [days]
1332
1333	CALL getin_p('STREAM_TCST', stream_tcst)
1334	!
1335	!Config Key = FLOOD_TCST
1336	!Config Desc = Time constant for the flood reservoir
1337	!Config If = RIVER_ROUTING
1338	!Config Def = 4.0
1339	!Config Help = This parameters allows the user to fix the
1340	!Config time constant (in days) of the flood reservoir
1341	!Config in order to get better river flows for
1342	!Config particular regions.
1343	!Config Units = [days]
1344
1345	CALL getin_p('SLOW_TCST', slow_tcst)
1346
1347	routing_file_type = "standard"
1348	CALL getin_p("routing_file_type", routing_file_type)
1349
1350	IF (TRIM(routing_file_type)=="standard") THEN
1351	topoind_factor=1000
1352	ELSE IF (TRIM(routing_file_type)=="merit") THEN
1353	topoind_factor=1
1354	ELSE
1355	CALL ipslerr_p(3,'routing_flow_init_local', &
1356	'Error in call routing_flow_init_local','getin routing_file_type : bad value, must be <standard> or <merit>','')
1357	ENDIF
1358
1359	split_routing=1
1360	CALL getin_p("SPLIT_ROUTING",split_routing)
1361
1362	CALL compute_coastline(contfrac, coastline)
1363
1364	CALL gather_omp(contfrac, contfrac_mpi)
1365	CALL gather_omp(area, area_mpi)
1366	CALL gather_omp(coastline, coastline_mpi)
1367
1368	ALLOCATE(is_coastline(nbp_mpi))
1369	is_coastline=coastline_mpi
1370
1371	IF (is_omp_root) THEN
1372	WHERE(contfrac_mpi <= epsilon) contfrac_mpi=0
1373	WHERE(contfrac_mpi >= 1-epsilon) contfrac_mpi=1
1374
1375	CALL xios_get_domain_attr("routing_domain", ni=ni, nj=nj) ! get routing domain dimension
1376
1377	nbpt_r= ni*nj
1378	nbpt_rp1= (ni+2)*(nj+2)
1379
1380	! Allocate module variable
1381	ALLOCATE(fast_reservoir_r(ni*nj))
1382	ALLOCATE(slow_reservoir_r(ni*nj))
1383	ALLOCATE(stream_reservoir_r(ni*nj))
1384
1385	ALLOCATE(is_lakeinflow_r(ni*nj))
1386	ALLOCATE(is_coastalflow_r(ni*nj))
1387	ALLOCATE(is_riverflow_r(ni*nj))
1388	ALLOCATE(is_streamflow_r(ni*nj))
1389
1390	ALLOCATE(topoind_r(ni*nj))
1391	ALLOCATE(routing_mask_r(ni*nj))
1392	ALLOCATE(basins_extended_r(nbpt_r))
1393	ALLOCATE(route_flow_rp1((ni+2)*(nj+2)))
1394
1395	ALLOCATE(routing_weight(nbp_mpi))
1396	ALLOCATE(routing_weight_in(nbp_mpi))
1397	ALLOCATE(unrouted_weight(nbp_mpi))
1398
1399	ALLOCATE(diag_r(nbpt_r)) ! for diags
1400	ALLOCATE(frac_coast_r(nbpt_r))
1401	ALLOCATE(frac_lake_r(nbpt_r))
1402	ALLOCATE(weight_coast_to_coast_r(nbpt_r))
1403	ALLOCATE(weight_coast_to_lake_r(nbpt_r))
1404	ALLOCATE(weight_lake_to_coast_r(nbpt_r))
1405	ALLOCATE(weight_lake_to_lake_r(nbpt_r))
1406
1407	ALLOCATE(trip_extended_r(ni,nj))
1408
1409	! correction on coastline in case of LAM, the cells near the frontier are considered as coastline
1410	! => interpolate model grid to routing grid and reinterpolate to model grid. Fractionnal cells will be considered as coastline too
1411	! diag=1
1412	! CALL xios_send_field("one", diag)
1413	! CALL xios_recv_field("tmp1_coastline", diag_r)
1414	! CALL xios_send_field("tmp2_coastline", diag_r)
1415	! CALL xios_recv_field("lam_coastline", diag)
1416	! WHERE(diag<1-epsilon) is_coastline=.TRUE.
1417
1418	diag=0
1419	WHERE(is_coastline) diag=1
1420	CALL xios_send_field("is_coastline", diag)
1421
1422	is_lakeinflow_r(:) = .FALSE.
1423	is_coastalflow_r(:) = .FALSE.
1424	is_riverflow_r(:) = .FALSE.
1425
1426
1427	INQUIRE(FILE="routing_start.nc", EXIST=file_exists)
1428
1429	IF (file_exists) THEN
1430	CALL xios_recv_field("fast_reservoir_start",fast_reservoir_r)
1431	CALL xios_recv_field("slow_reservoir_start",slow_reservoir_r)
1432	CALL xios_recv_field("stream_reservoir_start",stream_reservoir_r)
1433	ELSE
1434	fast_reservoir_r(:)=0
1435	slow_reservoir_r(:)=0
1436	stream_reservoir_r(:)=0
1437	ENDIF
1438
1439
1440	ALLOCATE(trip_rp1(0:ni+1,0:nj+1))
1441
1442	trip_rp1(:,:)=1e10
1443	CALL xios_recv_field("trip_r",trip_rp1(1:ni,1:nj)) ! recv trip array with halo of 1
1444	CALL xios_recv_field("trip_extended_r",trip_extended_r) ! recv extended trip array from file
1445	CALL xios_recv_field("topoind_r",topoind_r) ! recv topo index array from file
1446	CALL xios_recv_field("basins_extended_r",basins_extended_r) ! recv basins index from file
1447
1448	! CALL routing_flow_correct_riverflow(ni, nj, contfrac_mpi, coastline, trip_rp1(1:ni,1:nj), trip_extended_r, topoind_r)
1449	CALL new_routing_flow_correct_riverflow(ni, nj, contfrac_mpi, is_coastline, trip_rp1(1:ni,1:nj), trip_extended_r, topoind_r)
1450
1451	CALL xios_send_field("trip_update_r",trip_rp1(1:ni,1:nj)) ! send to xios trip array to update for halo
1452	CALL xios_recv_field("trip_rp1",trip_rp1) ! recv trip array with halo of 1
1453
1454
1455
1456	!! Compute the routing
1457	!! Loop on all point point of the local routing grid + halo
1458
1459	route_flow_rp1(:)=-1
1460
1461	DO j=0,nj+1
1462	jp1=j+1
1463	jm1=j-1
1464	DO i=0,ni+1
1465	ij_rp1=i+(ni+2)*j+1
1466	ij_r=i+ni*(j-1)
1467	ip1=i+1
1468	im1=i-1
1469
1470	IF (trip_rp1(i,j) < 100) THEN
1471
1472	IF (i>=1 .AND. i<=ni .AND. j>=1 .AND. j<=nj) routing_mask_r(ij_r)=.TRUE.
1473
1474	ir=-1 ; jr=-1 ! -> -1 for 97,98,99
1475	SELECT CASE (NINT(trip_rp1(i,j))) ! get the trip value for each points
1476	CASE (1)
1477	jr=jm1 ; ir=i ! north
1478	CASE (2)
1479	jr=jm1 ; ir=ip1 ! north-east
1480	CASE (3)
1481	jr=j ; ir=ip1 ! east
1482	CASE (4)
1483	jr=jp1 ; ir=ip1 ! south-east
1484	CASE (5)
1485	jr=jp1 ; ir=i ! south
1486	CASE(6)
1487	jr=jp1 ; ir=im1 ! south-west
1488	CASE (7)
1489	jr=j ; ir=im1 ! west
1490	CASE (8)
1491	jr=jm1 ; ir=im1 ! north-west
1492	CASE (97)
1493	IF ( i>0 .AND. i<ni+1 .AND. j>0 .AND. j<nj+1) THEN ! if inside my local domain
1494	is_lakeinflow_r(ij_r)=.TRUE. ! I am a lakeinflow point and route to myself
1495	jr=j ; ir=i
1496	ENDIF
1497	CASE (98)
1498	IF ( i>0 .AND. i<ni+1 .AND. j>0 .AND. j<nj+1) THEN ! if inside my local domain
1499	is_coastalflow_r(ij_r)=.TRUE. ! I am a coastal flow point and route to myself
1500	jr=j ; ir=i
1501	ENDIF
1502	CASE (99)
1503	IF ( i>0 .AND. i<ni+1 .AND. j>0 .AND. j<nj+1) THEN ! if inside my local domain
1504	is_riverflow_r(ij_r)=.TRUE. ! I am a riverflow point and route to myself
1505	jr=j ; ir=i
1506	ENDIF
1507	END SELECT
1508
1509	IF (ir<0 .OR. ir>ni+1 .OR. jr<0 .OR. jr>nj+1) THEN
1510	route_flow_rp1(ij_rp1)=-1 ! if route outside my local domain+halo, no routing (will be done by other process)
1511	ELSE
1512	IF (ir<1 .OR. ir>ni .OR. jr<1 .OR. jr>nj) THEN
1513	route_flow_rp1(ij_rp1)=-1 ! if route outside my local domain, no routing (will be done by other process)
1514	ELSE
1515	route_flow_rp1(ij_rp1)=ir+ni*(jr-1) ! define the cell where to flow
1516	IF (trip_rp1(ir,jr)>99) STOP 'Pb point not routed to outflow'
1517	ENDIF
1518	ENDIF
1519	ELSE
1520	IF (i>=1 .AND. i<=ni .AND. j>=1 .AND. j<=nj) routing_mask_r(ij_r)=.FALSE.
1521	ENDIF
1522	ENDDO
1523	ENDDO
1524	is_streamflow_r(:) = .NOT. (is_lakeinflow_r(:) .OR. is_coastalflow_r(:) .OR. is_riverflow_r(:)) .AND. routing_mask_r(:)
1525
1526	diag_r(:)=0
1527	WHERE(routing_mask_r) diag_r=1
1528	CALL xios_send_field("routing_mask_r", diag_r)
1529	CALL xios_recv_field("frac_routing", diag)
1530	WHERE (diag < epsilon) diag=0
1531	routing_weight_in (:) = 0
1532	unrouted_weight (:) = 0
1533	WHERE (diag > 0) routing_weight_in = contfrac_mpi*area_mpi/diag
1534	WHERE (diag == 0 ) unrouted_weight = contfrac_mpi*area_mpi
1535
1536
1537	! compute the number of coast cells to distribute lakeinflow onto
1538	diag=0
1539	WHERE (is_coastline) diag=1
1540	nb_coast_points=SUM(diag)
1541	CALL reduce_sum_mpi(nb_coast_points, total_coast_points)
1542	CALL bcast_mpi(total_coast_points)
1543
1544
1545	DO ij=1,nbp_mpi
1546	routing_weight(ij)=contfrac_mpi(ij)*area_mpi(ij)
1547	ENDDO
1548
1549	! looking for basins
1550	basins_count_mpi=0
1551	DO ij=1,nbpt_r
1552	IF (basins_extended_r(ij) > basins_count_mpi) basins_count_mpi=basins_extended_r(ij)
1553	ENDDO
1554	CALL MPI_ALLREDUCE(basins_count_mpi,basins_count,1,MPI_INT_ORCH, MPI_MAX,MPI_COMM_ORCH,ierr)
1555
1556	diag(:) = 0
1557	WHERE (is_coastline) diag=1
1558	CALL xios_send_field("mask_coastal",diag)
1559	CALL xios_recv_field("frac_coastal_r",diag_r)
1560	WHERE (diag_r > 100) ! missing_value
1561	diag_r=0
1562	ELSE WHERE (diag_r < epsilon)
1563	diag_r=0
1564	ELSEWHERE (diag_r > 1-epsilon)
1565	diag_r=1
1566	END WHERE
1567	frac_coast_r=diag_r ;
1568
1569	diag(:) = 0
1570	WHERE (.NOT. is_coastline) diag=1
1571	CALL xios_send_field("mask_lake",diag)
1572	CALL xios_recv_field("frac_lake_r",diag_r)
1573	WHERE (diag_r > 100) ! missing_value
1574	diag_r=0.
1575	ELSEWHERE (diag_r < epsilon)
1576	diag_r=0.
1577	ELSEWHERE (diag_r > 1-epsilon)
1578	diag_r=1.
1579	END WHERE
1580	frac_lake_r = diag_r
1581
1582	weight_coast_to_coast_r(:)=0
1583	weight_coast_to_lake_r(:)=0
1584	weight_lake_to_coast_r(:)=0
1585	weight_lake_to_lake_r(:)=0
1586
1587	DO ij=1,nbpt_r
1588	IF (is_riverflow_r(ij) .OR. is_coastalflow_r(ij) ) THEN
1589	IF (frac_coast_r(ij)==0 .AND. frac_lake_r(ij)==0) THEN
1590	PRINT*,"riverflow or costalflow point on routing grid can not be interpolated on orchidee grid, this migh not happen"
1591	STOP
1592	ELSE IF (frac_coast_r(ij)==0 .AND. frac_lake_r(ij)>0) THEN
1593	weight_coast_to_lake_r(ij) = 1./frac_lake_r(ij)
1594	weight_coast_to_coast_r(ij) = 0
1595	ELSE IF (frac_coast_r(ij)>0) THEN
1596	weight_coast_to_coast_r(ij) = 1./(frac_coast_r(ij))
1597	weight_coast_to_lake_r(ij) = 0
1598	ELSE
1599	ENDIF
1600	ELSE IF (is_lakeinflow_r(ij)) THEN
1601	IF ( frac_coast_r(ij)==0 .AND. frac_lake_r(ij)==0) THEN
1602	PRINT*,"lakeinflow on routing grid can not be interpolated on orchidee grid, this might not happen"
1603	STOP
1604	ELSE IF (frac_lake_r(ij)==0 .AND. frac_coast_r(ij)>0) THEN
1605	weight_lake_to_coast_r(ij) = 1./frac_coast_r(ij)
1606	weight_lake_to_lake_r(ij) = 0
1607	ELSE IF (frac_lake_r(ij)>0 ) THEN
1608	weight_lake_to_lake_r(ij) = 1./(frac_lake_r(ij))
1609	weight_lake_to_coast_r(ij) = 0
1610	ENDIF
1611	ENDIF
1612	ENDDO
1613
1614	CALL compute_basins_area(nbpt_r,nbpt_rp1)
1615
1616	ELSE
1617
1618	nbpt_r=0
1619
1620	ENDIF !is_omp_root
1621
1622
1623	END SUBROUTINE routing_flow_init_local
1624
1625
1626	SUBROUTINE compute_basins_area(nbpt_r,nbpt_rp1 )
1627	USE xios
1628	USE routing_native_para
1629	IMPLICIT NONE
1630	INCLUDE 'mpif.h'
1631	INTEGER, INTENT(IN) :: nbpt_r
1632	INTEGER, INTENT(IN) :: nbpt_rp1
1633	REAL(r_std),PARAMETER :: PI=atan(1.)*4
1634	REAL(r_std),PARAMETER :: earth_radius=6.371009E6 ! in meter
1635
1636	INTEGER :: ni, nj
1637	REAL(r_std),ALLOCATABLE :: lon(:), lat(:)
1638	REAL(r_std) :: area_rp1(nbpt_rp1)
1639	REAL(r_std) :: send_area_rp1(nbpt_rp1)
1640	REAL(r_std) :: recv_area_r(nbpt_rp1)
1641	REAL :: delta_lon, delta_lat
1642	INTEGER :: i,j,ij,it,ig,ierr
1643	REAL(r_std) :: sum_area
1644
1645	ALLOCATE(basins_area_r(nbpt_r))
1646
1647	CALL xios_get_domain_attr("routing_domain", ni=ni, nj=nj)
1648	ALLOCATE(lon(ni),lat(nj))
1649	CALL xios_get_domain_attr("routing_domain",lonvalue_1d=lon, latvalue_1d=lat)
1650
1651	! compute routing cell area
1652	delta_lon=ABS(lon(2)-lon(1))*Pi/180.
1653	delta_lat=ABS(lat(2)-lat(1))*Pi/180.
1654
1655	DO j=1,nj
1656	DO i=1,ni
1657	ig=(j+1-1)*(ni+2)+i+1
1658	area_rp1(ig) = delta_londelta_latcos(lat(j)Pi/180.)earth_radius*earth_radius ! pretty good approximation
1659	ij = (j-1)*ni+i
1660	basins_area_r(ij) = delta_londelta_latcos(lat(j)Pi/180.)earth_radius*earth_radius ! pretty good approximation
1661	ENDDO
1662	ENDDO
1663
1664	CALL update_halo(area_rp1)
1665
1666
1667	send_area_rp1(:) = area_rp1(:)
1668	sum_area=1.
1669	DO WHILE (sum_area/=0)
1670
1671	recv_area_r(:) = 0
1672	DO ig=1,nbpt_rp1
1673	IF ( route_flow_rp1(ig) > 0 ) THEN
1674	basins_area_r(route_flow_rp1(ig)) = basins_area_r(route_flow_rp1(ig)) + send_area_rp1(ig)
1675	recv_area_r(route_flow_rp1(ig)) = recv_area_r(route_flow_rp1(ig)) + send_area_rp1(ig)
1676	ENDIF
1677	ENDDO
1678
1679	DO j=1,nj
1680	DO i=1,ni
1681	ij=((j-1)*ni)+i
1682	ig=(j+1-1)*(ni+2)+i+1
1683	send_area_rp1(ig) = recv_area_r(ij)
1684	IF (is_riverflow_r(ij) .OR. is_coastalflow_r(ij) .OR. is_lakeinflow_r(ij)) send_area_rp1(ig)=0
1685	ENDDO
1686	ENDDO
1687
1688	CALL update_halo(send_area_rp1)
1689	sum_area = sum(send_area_rp1)
1690	CALL MPI_ALLREDUCE(MPI_IN_PLACE, sum_area, 1, MPI_REAL_ORCH, MPI_SUM, MPI_COMM_ORCH, ierr)
1691
1692	ENDDO
1693
1694	END SUBROUTINE compute_basins_area
1695
1696
1697	SUBROUTINE initialize_stations(dt_routing)
1698	USE xios
1699	IMPLICIT NONE
1700	INCLUDE 'mpif.h'
1701	REAL(r_std), INTENT (in) :: dt_routing !! Routing time step (s)
1702
1703	REAL(r_std),PARAMETER :: PI=atan(1.)*4
1704	REAL(r_std),PARAMETER :: earth_radius=6.371009E6 ! in meter
1705	INTEGER :: ni, nj
1706	REAL(r_std),ALLOCATABLE :: lon(:), lat(:)
1707
1708	INTEGER :: nb_mouth
1709	CHARACTER(LEN=60),ALLOCATABLE :: mouth(:)
1710	INTEGER,ALLOCATABLE :: mouth_id(:)
1711	INTEGER,ALLOCATABLE :: mouth_index(:)
1712
1713	REAL,ALLOCATABLE :: station_lonlat(:,:)
1714	REAL,ALLOCATABLE :: station_prec(:)
1715	REAL,ALLOCATABLE :: station_area(:)
1716	INTEGER :: station_index_prec, station_index_prec_i, station_index_prec_j
1717	INTEGER :: station_index_area, station_index_area_i, station_index_area_j
1718	LOGICAL :: has_lonlat
1719	LOGICAL :: has_area
1720	INTEGER ::i,j,r,k
1721
1722	TYPE(xios_duration) :: ts
1723	TYPE(xios_scalar) :: scalar_hdl
1724	TYPE(xios_scalargroup) :: scalar_def_hdl
1725	TYPE(xios_field) :: field_hdl
1726	TYPE(xios_fieldgroup) :: field_def_hdl, fieldgroup_hdl
1727	TYPE(xios_file) :: file_hdl
1728	TYPE(xios_filegroup) :: file_def_hdl
1729	TYPE(xios_date) :: start_date, time_origin
1730	CHARACTER(LEN=20) :: calendar_type
1731
1732	CHARACTER(LEN=256) :: str_station_ind
1733	REAL :: lon_a,lat_a,lon_b,lat_b,dist,max_area, max_max_area, min_dist, min_area, min_min_area
1734	INTEGER :: min_dist_index, min_dist_index_i, min_dist_index_j
1735	INTEGER :: ierr
1736
1737
1738	CALL xios_get_domain_attr("routing_domain", ni=ni, nj=nj)
1739	ALLOCATE(lon(ni),lat(nj))
1740	CALL xios_get_domain_attr("routing_domain",lonvalue_1d=lon, latvalue_1d=lat)
1741
1742	nb_station=0
1743	CALL getin("nb_station",nb_station)
1744	ALLOCATE(station(nb_station))
1745	ALLOCATE(station_lonlat(nb_station,2))
1746	ALLOCATE(station_index(nb_station))
1747	ALLOCATE(station_prec(nb_station))
1748	ALLOCATE(station_area(nb_station))
1749
1750	DO k=1,nb_station
1751	WRITE(str_station_ind,*) k
1752	str_station_ind=ADJUSTL(str_station_ind)
1753	CALL getin("station"//TRIM(str_station_ind)//"_id",station(k))
1754	station_lonlat(k,:) = 1.1e10
1755	has_lonlat=.TRUE.
1756	CALL getin("station"//TRIM(str_station_ind)//"_coor",station_lonlat(k,:))
1757	IF (station_lonlat(k,1)>1e10 .OR. station_lonlat(k,2)>1e10) has_lonlat=.FALSE.
1758	station_prec(k)=50000
1759	CALL getin("station"//TRIM(str_station_ind)//"_prec",station_prec(k))
1760	station_area(k) = -1
1761	has_area=.TRUE.
1762	CALL getin("station"//TRIM(str_station_ind)//"_area",station_area(k))
1763	IF (station_area(k)==-1) has_area=.FALSE.
1764
1765	lon_a = station_lonlat(k,1)*Pi/180.
1766	lat_a = station_lonlat(k,2)*Pi/180.
1767
1768	max_area=0
1769	min_dist=HUGE(min_dist)
1770	min_area=HUGE(min_area)
1771
1772	DO j=1,nj
1773	DO i=1,ni
1774	r=((j-1)*ni)+i
1775	IF (routing_mask_r(r)) THEN
1776	lon_b=lon(i)*Pi/180.
1777	lat_b=lat(j)*Pi/180.
1778	dist=earth_radiusacos(sin(lon_a)sin(lon_b)+cos(lon_a)cos(lon_b)cos(lat_b-lat_a))
1779	IF (dist<station_prec(k)) THEN
1780	IF (dist<min_dist) THEN
1781	min_dist=dist
1782	min_dist_index = r
1783	min_dist_index_i = i
1784	min_dist_index_j = j
1785	ENDIF
1786
1787	IF (basins_area_r(r) > max_area) THEN
1788	max_area = basins_area_r(r)
1789	station_index_prec = r
1790	station_index_prec_i = i
1791	station_index_prec_j = j
1792	ENDIF
1793
1794	IF (has_area) THEN
1795	IF (ABS(basins_area_r(r)-station_area(k)*1e6) < min_area) THEN
1796	min_area = ABS(basins_area_r(r)-station_area(k)*1e6)
1797	station_index_area = r
1798	station_index_area_i = i
1799	station_index_area_j = j
1800	ENDIF
1801	ENDIF
1802
1803	ENDIF
1804	ENDIF
1805	ENDDO
1806	ENDDO
1807
1808	CALL MPI_ALLREDUCE(max_area, max_max_area, 1, MPI_REAL_ORCH, MPI_MAX, MPI_COMM_ORCH, ierr)
1809	IF (max_area /= max_max_area) station_index_prec=-1
1810	IF (max_area == 0) station_index_prec=-1
1811
1812	IF (has_area) THEN
1813	CALL MPI_ALLREDUCE(min_area, min_min_area, 1, MPI_REAL_ORCH, MPI_MIN, MPI_COMM_ORCH, ierr)
1814	IF (min_area /= min_min_area) station_index_area=-1
1815	IF (min_area == HUGE(min_area)) station_index_area=-1
1816	ELSE
1817	station_index_area=-1
1818	ENDIF
1819
1820	IF (station_index_prec/=-1) THEN
1821	PRINT*,"Found station ",TRIM(station(k))," based on coordinates lon-lat",lon(min_dist_index_i),lat(min_dist_index_j), &
1822	" with basin_area (km2) ", basins_area_r(min_dist_index)/1000/1000
1823	PRINT*,"=> at coordinate ",lon(station_index_prec_i),lat(station_index_prec_j)," basin_area(km2) ", &
1824	basins_area_r(station_index_prec)/1000/1000
1825	ENDIF
1826
1827	IF (station_index_area/=-1) THEN
1828	PRINT*,"Found station ",TRIM(station(k))," based on basin area (km2)", station_area(k)
1829	PRINT*,"=> at coordinate ",lon(station_index_area_i),lat(station_index_area_j)," basin_area(km2) ", &
1830	basins_area_r(station_index_area)/1000/1000
1831	ENDIF
1832
1833	IF (station_index_area/=-1 ) THEN
1834	station_index(k) = station_index_area
1835	ELSE
1836	station_index(k) = station_index_prec
1837	ENDIF
1838
1839	ENDDO
1840
1841	CALL xios_get_calendar_type(calendar_type)
1842	CALL xios_get_start_date(start_date)
1843	CALL xios_get_time_origin(time_origin)
1844
1845	CALL xios_context_initialize("orchidee_routing_out",MPI_COMM_ORCH)
1846	CALL xios_orchidee_change_context("orchidee_routing_out")
1847	ts.second=dt_routing
1848	calendar_type="gregorian"
1849	CALL xios_define_calendar(type=calendar_type,start_date=start_date,time_origin=time_origin, timestep=ts )
1850
1851	CALL xios_get_handle("scalar_definition",scalar_def_hdl)
1852
1853	! CALL xios_get_handle("mouths",field_def_hdl)
1854	! DO k=1,nb_mouth
1855	! CALL xios_add_child(scalar_def_hdl,scalar_hdl,TRIM(mouth(k))//"_scalar")
1856	! CALL xios_set_attr(scalar_hdl,label=TRIM(mouth(k)))
1857	! CALL xios_add_child(field_def_hdl,field_hdl, TRIM(mouth(k)))
1858	! CALL xios_set_attr(field_hdl,scalar_ref=TRIM(mouth(k))//"_scalar")
1859	! CALL xios_set_attr(field_hdl, freq_op = freq_op*xios_second)
1860	! ENDDO
1861
1862	CALL xios_get_handle("field_definition",field_def_hdl)
1863	CALL xios_add_child(field_def_hdl, fieldgroup_hdl, "stations")
1864	DO k=1,nb_station
1865	CALL xios_add_child(scalar_def_hdl,scalar_hdl,TRIM(station(k))//"_scalar")
1866	CALL xios_set_attr(scalar_hdl,label=TRIM(station(k)))
1867	CALL xios_add_child(fieldgroup_hdl,field_hdl, TRIM(station(k)))
1868	CALL xios_set_attr(field_hdl,scalar_ref=TRIM(station(k))//"_scalar")
1869	ENDDO
1870
1871	CALL xios_get_handle("file_definition",file_def_hdl)
1872	CALL xios_add_child(file_def_hdl, file_hdl, "stations")
1873	CALL xios_set_attr(file_hdl, type="one_file", output_freq=ts, sync_freq=ts, enabled=.TRUE.)
1874	CALL xios_add_child(file_hdl, fieldgroup_hdl)
1875	CALL xios_set_attr(fieldgroup_hdl, group_ref="stations", operation="average")
1876	CALL xios_close_context_definition()
1877	CALL xios_orchidee_change_context("orchidee")
1878	END SUBROUTINE initialize_stations
1879
1880	SUBROUTINE routing_flow_main(dt_routing, contfrac) !
1881
1882	USE xios
1883	USE grid, ONLY : area
1884	USE routing_native_para
1885	IMPLICIT NONE
1886	INCLUDE "mpif.h"
1887
1888	!! 0 Variable and parameter description
1889	!! 0.1 Input variables
1890	REAL(r_std), INTENT (in) :: dt_routing !! Routing time step (s)
1891	REAL(r_std), INTENT (in) :: contfrac(nbpt) !! Routing time step (s)
1892
1893	!! 0.4 Local variables
1894	REAL(r_std) :: runoff(nbp_mpi) !! Grid-point runoff (kg/dt)
1895	REAL(r_std) :: runoff_in(nbp_mpi) !! Grid-point runoff (kg/dt)
1896	REAL(r_std) :: drainage(nbp_mpi) !! Grid-point drainage (kg/dt)
1897	REAL(r_std) :: drainage_in(nbp_mpi) !! Grid-point drainage (kg/dt)
1898	REAL(r_std) :: riverflow(nbp_mpi)
1899	REAL(r_std) :: coastalflow(nbp_mpi)
1900	REAL(r_std) :: lakeinflow(nbp_mpi)
1901	REAL(r_std) :: area_mpi(nbp_mpi) ! cell area
1902	REAL(r_std) :: contfrac_mpi(nbp_mpi) ! cell area
1903	REAL(r_std) :: flow_coast(nbp_mpi)
1904	REAL(r_std) :: flow_lake(nbp_mpi)
1905
1906	! diags
1907	REAL(r_std) :: fast_diag_mpi(nbp_mpi)
1908	REAL(r_std) :: slow_diag_mpi(nbp_mpi)
1909	REAL(r_std) :: stream_diag_mpi(nbp_mpi)
1910	REAL(r_std) :: hydrographs_diag_mpi(nbp_mpi)
1911
1912	! from input model -> routing_grid
1913	REAL(r_std) :: runoff_r(nbpt_r) !! Grid-point runoff (kg/m^2/dt)
1914	REAL(r_std) :: drainage_r(nbpt_r) !! Grid-point drainage (kg/m^2/dt)
1915
1916	REAL(r_std), DIMENSION(nbpt_r) :: fast_flow_r !! Outflow from the fast reservoir (kg/dt)
1917	REAL(r_std), DIMENSION(nbpt_r) :: slow_flow_r !! Outflow from the slow reservoir (kg/dt)
1918	REAL(r_std), DIMENSION(nbpt_r) :: stream_flow_r !! Outflow from the stream reservoir (kg/dt)
1919	REAL(r_std), DIMENSION(nbpt_r) :: hydrographs_r !! hydrograph (kg/dt)
1920	REAL(r_std), DIMENSION(nbpt_r) :: transport_r !! Water transport between basins (kg/dt)
1921
1922	INTEGER(i_std) :: ig !! Indices (unitless)
1923	INTEGER(i_std) :: isplit
1924
1925	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)
1926
1927	REAL(r_std), DIMENSION(nbpt_r) :: lakeinflow_r
1928	REAL(r_std), DIMENSION(nbpt_r) :: coastalflow_r
1929	REAL(r_std), DIMENSION(nbpt_r) :: riverflow_r
1930
1931	REAL(r_std), DIMENSION(nbpt_r) :: flow_r
1932	REAL(r_std), DIMENSION(nbpt_rp1) :: flow_rp1
1933	REAL(r_std) :: flow !! Outflow computation for the reservoirs (kg/dt)
1934	REAL(r_std) :: basins_riverflow_mpi(0:basins_count)
1935	REAL(r_std) :: basins_riverflow(0:basins_count)
1936	REAL(r_std) :: water_balance_before, water_balance_after
1937	REAL(r_std) :: sum_water_before, sum_water_after
1938	REAL(r_std) :: value
1939	INTEGER(i_std) :: k
1940	INTEGER :: ierr
1941
1942	!_ ================================================================================================================================
1943	!> The outflow fluxes from the three reservoirs are computed.
1944	!> The outflow of volume of water Vi into the reservoir i is assumed to be linearly related to its volume.
1945	!> The water travel simulated by the routing scheme is dependent on the water retention index topo_resid
1946	!> given by a 0.5 degree resolution map for each pixel performed from a simplification of Manning's formula
1947	!> (Dingman, 1994; Ducharne et al., 2003).
1948	!> The resulting product of tcst (in day/m) and topo_resid (in m) represents the time constant (day)
1949	!> which is an e-folding time, the time necessary for the water amount
1950	!> in the stream reservoir to decrease by a factor e. Hence, it gives an order of
1951	!> magnitude of the travel time through this reservoir between
1952	!> the sub-basin considered and its downstream neighbor.
1953
1954
1955
1956	CALL gather_omp(runoff_mean,runoff)
1957	CALL gather_omp(drainage_mean, drainage)
1958	CALL gather_omp(area, area_mpi)
1959	CALL gather_omp(contfrac, contfrac_mpi)
1960
1961	IF (is_omp_root) THEN
1962	CALL xios_send_field("routing_basins_area",basins_area_r)
1963
1964	hydrographs_r(:)=0
1965	! water balance before
1966
1967	sum_water_before = sum((runoff(:)+drainage(:))*routing_weight(:))+sum(fast_reservoir_r(:)+slow_reservoir_r(:)+stream_reservoir_r(:))
1968	CALL MPI_ALLREDUCE(sum_water_before, water_balance_before,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr)
1969
1970
1971	runoff_in=runoff*routing_weight_in
1972	CALL xios_send_field("routing_runoff",runoff_in) ! interp conservative model -> routing
1973	CALL xios_recv_field("routing_runoff_r",runoff_r)
1974	WHERE(.NOT. routing_mask_r) runoff_r = 0
1975
1976	drainage_in=drainage*routing_weight_in
1977	CALL xios_send_field("routing_drainage",drainage_in) ! interp conservative model -> routing
1978	CALL xios_recv_field("routing_drainage_r",drainage_r)
1979	WHERE(.NOT. routing_mask_r) drainage_r = 0
1980
1981	CALL MPI_ALLREDUCE(sum(runoff*(routing_weight-unrouted_weight)),sum_water_before,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr)
1982	CALL MPI_ALLREDUCE(sum(runoff_r),sum_water_after,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr)
1983	IF (sum_water_after/=0) THEN
1984	IF (is_mpi_root) PRINT *,"runoff fixer : ", sum_water_before/sum_water_after
1985	DO ig=1,nbpt_r
1986	runoff_r(ig) = runoff_r(ig) * sum_water_before/sum_water_after
1987	ENDDO
1988	ENDIF
1989
1990	CALL MPI_ALLREDUCE(sum(drainage*(routing_weight-unrouted_weight)),sum_water_before,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr)
1991	CALL MPI_ALLREDUCE(sum(drainage_r),sum_water_after,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr)
1992	IF (sum_water_after/=0) THEN
1993	IF (is_mpi_root) PRINT *,"drainage fixer : ", sum_water_before/sum_water_after
1994	DO ig=1,nbpt_r
1995	drainage_r(ig) = drainage_r(ig) * sum_water_before/sum_water_after
1996	ENDDO
1997	ENDIF
1998
1999
2000	CALL MPI_ALLREDUCE(sum((runoff+drainage)*(routing_weight-unrouted_weight)),sum_water_before,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr)
2001	CALL MPI_ALLREDUCE(sum(runoff_r+drainage_r),sum_water_after,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr)
2002
2003	IF (is_mpi_root) PRINT *,"Loose water by interpolation ; before : ", sum_water_before," ; after : ",sum_water_after, &
2004	" ; delta : ", 100.(sum_water_after-sum_water_before)/(0.5(sum_water_after+sum_water_before)),"%"
2005
2006
2007	runoff_in(:) = runoff(:)*unrouted_weight(:)
2008	drainage_in(:) = drainage(:)*unrouted_weight(:)
2009	runoff_r(:) = runoff_r(:) / split_routing
2010	drainage_r(:) = drainage_r(:) / split_routing
2011	hydrographs_r(:) = 0
2012	transport_r(:)=0
2013
2014	DO isplit=1,split_routing
2015
2016	DO ig=1,nbpt_r
2017	IF ( routing_mask_r(ig) ) THEN
2018
2019	flow = MIN(fast_reservoir_r(ig)/((topoind_r(ig)topoind_factor/1000.)fast_tcst*one_day/(dt_routing/split_routing)),&
2020	& fast_reservoir_r(ig)-min_sechiba)
2021	fast_flow_r(ig) = MAX(flow, zero)
2022	!
2023	flow = MIN(slow_reservoir_r(ig)/((topoind_r(ig)topoind_factor/1000.)slow_tcst*one_day/(dt_routing/split_routing)),&
2024	& slow_reservoir_r(ig)-min_sechiba)
2025	slow_flow_r(ig) = MAX(flow, zero)
2026	!
2027	flow = MIN(stream_reservoir_r(ig)/((topoind_r(ig)topoind_factor/1000.)stream_tcst*one_day/(dt_routing/split_routing)),&
2028	& stream_reservoir_r(ig)-min_sechiba)
2029	stream_flow_r(ig) = MAX(flow, zero)
2030	!
2031	!
2032	ELSE
2033	fast_flow_r(ig) = zero
2034	slow_flow_r(ig) = zero
2035	stream_flow_r(ig) = zero
2036	ENDIF
2037	ENDDO
2038
2039
2040	!-
2041	!- Compute the transport
2042	!-
2043
2044	DO ig=1,nbpt_r
2045	flow_rp1(r_to_rp1(ig))=fast_flow_r(ig) + slow_flow_r(ig) + stream_flow_r(ig)
2046	ENDDO
2047
2048	CALL update_halo(flow_rp1) ! transfert halo
2049
2050	DO ig=1,nbpt_rp1
2051	IF ( route_flow_rp1(ig) > 0 ) THEN
2052	transport_r(route_flow_rp1(ig))=transport_r(route_flow_rp1(ig))+ flow_rp1(ig)
2053	ENDIF
2054	ENDDO
2055
2056
2057	DO ig=1,nbpt_r
2058	IF ( routing_mask_r(ig) ) THEN
2059	fast_reservoir_r(ig) = fast_reservoir_r(ig) + runoff_r(ig) - fast_flow_r(ig)
2060	slow_reservoir_r(ig) = slow_reservoir_r(ig) + drainage_r(ig) - slow_flow_r(ig)
2061	stream_reservoir_r(ig) = stream_reservoir_r(ig) - stream_flow_r(ig)
2062	IF (is_streamflow_r(ig)) THEN
2063	stream_reservoir_r(ig)= stream_reservoir_r(ig) + transport_r(ig)
2064	transport_r(ig) = 0
2065	ENDIF
2066
2067	IF ( stream_reservoir_r(ig) .LT. zero ) THEN
2068	IF ( check_reservoir ) THEN
2069	WRITE(numout,*) "WARNING : negative stream reservoir at :", ig, ". Problem is being corrected."
2070	WRITE(numout,*) "stream_reservoir, transport, stream_flow,: ", &
2071	stream_reservoir_r(ig), transport_r(ig), stream_flow_r(ig)
2072	ENDIF
2073	fast_reservoir_r(ig) = fast_reservoir_r(ig) + stream_reservoir_r(ig)
2074	stream_reservoir_r(ig) = zero
2075	ENDIF
2076	!
2077	IF ( fast_reservoir_r(ig) .LT. zero ) THEN
2078	IF ( check_reservoir ) THEN
2079	WRITE(numout,*) "WARNING : negative fast reservoir at :", ig, ". Problem is being corrected."
2080	WRITE(numout,*) "fast_reservoir, runoff, fast_flow : ", fast_reservoir_r(ig), &
2081	&runoff_r(ig), fast_flow_r(ig)
2082	ENDIF
2083	slow_reservoir_r(ig) = slow_reservoir_r(ig) + fast_reservoir_r(ig)
2084	fast_reservoir_r(ig) = zero
2085	ENDIF
2086
2087	IF ( slow_reservoir_r(ig) .LT. - min_sechiba ) THEN
2088	WRITE(numout,*) 'WARNING : There is a negative reservoir at :', ig
2089	WRITE(numout,*) 'WARNING : slowr, slow_flow, drainage', &
2090	& slow_reservoir_r(ig), slow_flow_r(ig), drainage_r(ig)
2091	CALL ipslerr_p(2, 'routing_simple_flow', 'WARNING negative slow_reservoir.','','')
2092	ENDIF
2093	ENDIF
2094	ENDDO
2095
2096	DO ig=1,nbpt_r
2097	IF ( routing_mask_r(ig) ) THEN
2098	hydrographs_r(ig)=hydrographs_r(ig)+fast_flow_r(ig)+slow_flow_r(ig)+stream_flow_r(ig)
2099	ENDIF
2100	ENDDO
2101
2102	ENDDO ! isplit
2103
2104	lakeinflow_r(:)=0
2105	coastalflow_r(:)=0
2106	riverflow_r(:)=0
2107	basins_riverflow_mpi(:)=0
2108
2109	DO ig=1,nbpt_r
2110	IF ( routing_mask_r(ig) ) THEN
2111	IF (is_lakeinflow_r(ig)) THEN
2112	lakeinflow_r(ig) = transport_r(ig)
2113	basins_riverflow_mpi(basins_extended_r(ig)) = basins_riverflow_mpi(basins_extended_r(ig))+lakeinflow_r(ig)
2114	ENDIF
2115
2116	IF (is_coastalflow_r(ig)) THEN
2117	coastalflow_r(ig) = transport_r(ig)
2118	basins_riverflow_mpi(basins_extended_r(ig)) = &
2119	basins_riverflow_mpi(basins_extended_r(ig))+coastalflow_r(ig)
2120	ENDIF
2121
2122	IF (is_riverflow_r(ig)) THEN
2123	riverflow_r(ig) = transport_r(ig)
2124	basins_riverflow_mpi(basins_extended_r(ig)) = &
2125	basins_riverflow_mpi(basins_extended_r(ig))+riverflow_r(ig)
2126	ENDIF
2127	ENDIF
2128	ENDDO
2129
2130	! now send riverflow, coastalflow and lakeinflow on orchidee grid
2131
2132	coastalflow(:)=0.
2133	riverflow(:)=0.
2134	lakeinflow(:)=0.
2135
2136	CALL xios_send_field("routing_coastalflow_to_coast_r" ,coastalflow_r*weight_coast_to_coast_r)
2137	CALL xios_recv_field("routing_coastalflow_to_coast" ,flow_coast)
2138	WHERE(.NOT.is_coastline) flow_coast=0
2139
2140	CALL xios_send_field("routing_coastalflow_to_lake_r" ,coastalflow_r*weight_coast_to_lake_r)
2141	CALL xios_recv_field("routing_coastalflow_to_lake" ,flow_lake)
2142	WHERE(is_coastline) flow_lake=0.
2143
2144	CALL MPI_ALLREDUCE(sum(coastalflow_r),sum_water_before,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr)
2145	CALL MPI_ALLREDUCE(sum(flow_coast+flow_lake),sum_water_after,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr)
2146	IF (sum_water_after/=0) THEN
2147	IF (is_mpi_root) PRINT *,"coastalflow fixer : ", sum_water_before/sum_water_after
2148	flow_coast(:)=flow_coast(:)*(sum_water_before/sum_water_after)
2149	flow_lake(:)=flow_lake(:)*(sum_water_before/sum_water_after)
2150	ENDIF
2151
2152	coastalflow=coastalflow+flow_coast
2153	lakeinflow=lakeinflow+flow_lake
2154
2155	CALL xios_send_field("routing_riverflow_to_coast_r" ,riverflow_r*weight_coast_to_coast_r)
2156	CALL xios_recv_field("routing_riverflow_to_coast" ,flow_coast)
2157	WHERE(.NOT.is_coastline) flow_coast=0
2158
2159	CALL xios_send_field("routing_riverflow_to_lake_r" ,riverflow_r*weight_coast_to_lake_r)
2160	CALL xios_recv_field("routing_riverflow_to_lake" ,flow_lake)
2161	WHERE(is_coastline) flow_lake=0.
2162
2163	CALL MPI_ALLREDUCE(sum(riverflow_r),sum_water_before,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr)
2164	CALL MPI_ALLREDUCE(sum(flow_coast+flow_lake),sum_water_after,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr)
2165	IF (sum_water_after/=0) THEN
2166	IF (is_mpi_root) PRINT *,"riverflow fixer : ", sum_water_before/sum_water_after
2167	flow_coast(:)=flow_coast(:)*(sum_water_before/sum_water_after)
2168	flow_lake(:)=flow_lake(:)*(sum_water_before/sum_water_after)
2169	ENDIF
2170
2171	riverflow=riverflow+flow_coast
2172	lakeinflow=lakeinflow+flow_lake
2173
2174	CALL xios_send_field("routing_lakeinflow_to_coast_r" ,lakeinflow_r*weight_lake_to_coast_r)
2175	CALL xios_recv_field("routing_lakeinflow_to_coast" ,flow_coast)
2176	WHERE(.NOT.is_coastline) flow_coast=0
2177
2178	CALL xios_send_field("routing_lakeinflow_to_lake_r" ,lakeinflow_r*weight_lake_to_lake_r)
2179	CALL xios_recv_field("routing_lakeinflow_to_lake" ,flow_lake)
2180	WHERE(is_coastline) flow_lake=0.
2181
2182	CALL MPI_ALLREDUCE(sum(lakeinflow_r),sum_water_before,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr)
2183	CALL MPI_ALLREDUCE(sum(flow_coast+flow_lake),sum_water_after,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr)
2184	IF (sum_water_after/=0) THEN
2185	IF (is_mpi_root) PRINT *,"lakeinflow fixer : ", sum_water_before/sum_water_after
2186	flow_coast(:)=flow_coast(:)*(sum_water_before/sum_water_after)
2187	flow_lake(:)=flow_lake(:)*(sum_water_before/sum_water_after)
2188	ENDIF
2189
2190	coastalflow=coastalflow+flow_coast+(runoff+drainage)*unrouted_weight
2191	lakeinflow=lakeinflow+flow_lake
2192
2193
2194	! WHERE(is_coastline) coastalflow = coastalflow + runoff_in + drainage_in
2195	! WHERE(.NOT.is_coastline) lakeinflow = lakeinflow + runoff_in + drainage_in
2196
2197	sum_water_after = sum(coastalflow(:) + riverflow(:) + lakeinflow(:))+sum(fast_reservoir_r(:)+slow_reservoir_r(:)+stream_reservoir_r(:))
2198	CALL MPI_ALLREDUCE(sum_water_after, water_balance_after,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr)
2199
2200	IF (is_mpi_root) PRINT *,"routing water Balance ; before : ", water_balance_before," ; after : ",water_balance_after, &
2201	" ; delta : ", 100(water_balance_after-water_balance_before)/(0.5(water_balance_after+water_balance_before)),"%"
2202	! diag
2203	CALL xios_send_field("routing_fast_reservoir_r" , fast_reservoir_r)
2204	CALL xios_send_field("routing_slow_reservoir_r" , slow_reservoir_r)
2205	CALL xios_send_field("routing_stream_reservoir_r" , stream_reservoir_r)
2206	CALL xios_send_field("routing_riverflow_r" , riverflow_r)
2207	CALL xios_send_field("routing_coastalflow_r" , coastalflow_r)
2208	CALL xios_send_field("routing_lakeinflow_r" , lakeinflow_r)
2209	CALL xios_send_field("out_flow",lakeinflow+coastalflow+riverflow)
2210	CALL xios_send_field("routing_hydrographs_r", (hydrographs_r+lakeinflow_r+coastalflow_r+riverflow_r)/1000./dt_routing)
2211	CALL xios_send_field("routing_riverflow" , riverflow)
2212	CALL xios_send_field("routing_coastalflow" , coastalflow)
2213	CALL xios_send_field("routing_lakeinflow" , lakeinflow)
2214	!
2215	! stations
2216	CALL xios_orchidee_change_context("orchidee_routing_out")
2217	station_ts = station_ts +1
2218	CALL xios_update_calendar(station_ts)
2219	DO k=1,nb_station
2220	IF (station_index(k) /=-1) THEN
2221	value = hydrographs_r(station_index(k))/1000./dt_routing
2222	ELSE
2223	value = 0
2224	ENDIF
2225	CALL MPI_ALLREDUCE(MPI_IN_PLACE,value,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr)
2226	CALL xios_send_field(TRIM(station(k)),value)
2227	ENDDO
2228
2229	CALL xios_orchidee_change_context("orchidee") !! return on orchidee context
2230
2231
2232	!!! for orchidee history diag
2233
2234	CALL xios_send_field("routing_fast_diag_r", fast_reservoir_r)
2235	CALL xios_recv_field("routing_fast_diag", fast_diag_mpi)
2236	CALL xios_send_field("routing_slow_diag_r", slow_reservoir_r)
2237	CALL xios_recv_field("routing_slow_diag", slow_diag_mpi)
2238	CALL xios_send_field("routing_stream_diag_r", stream_reservoir_r)
2239	CALL xios_recv_field("routing_stream_diag", stream_diag_mpi)
2240	CALL xios_send_field("routing_hydrographs_diag_r", hydrographs_r+lakeinflow_r+coastalflow_r+riverflow_r)
2241	CALL xios_recv_field("routing_hydrographs_diag", hydrographs_diag_mpi)
2242
2243	fast_diag_mpi=fast_diag_mpi/(area_mpi*contfrac_mpi) !! kg => kg/m^2
2244	slow_diag_mpi=slow_diag_mpi/(area_mpi*contfrac_mpi) !! kg => kg/m^2
2245	stream_diag_mpi=stream_diag_mpi/(area_mpi*contfrac_mpi) !! kg => kg/m^2
2246	hydrographs_diag_mpi = hydrographs_diag_mpi/1000.
2247
2248	ENDIF ! is_omp_root
2249
2250	CALL scatter_omp(riverflow, riverflow_mean)
2251	CALL scatter_omp(coastalflow, coastalflow_mean)
2252	CALL scatter_omp(lakeinflow, lakeinflow_mean)
2253	CALL scatter_omp(fast_diag_mpi, fast_diag)
2254	CALL scatter_omp(slow_diag_mpi, slow_diag)
2255	CALL scatter_omp(stream_diag_mpi, stream_diag)
2256	CALL scatter_omp(hydrographs_diag_mpi, hydrographs_diag)
2257
2258	CALL routing_flow_reset_mean
2259
2260	END SUBROUTINE routing_flow_main
2261
2262	SUBROUTINE routing_flow_diags(dt_sechiba)
2263	USE xios
2264	IMPLICIT NONE
2265	REAL(r_std) :: dt_sechiba
2266
2267	CALL xios_orchidee_send_field("fastr",fast_diag)
2268	CALL xios_orchidee_send_field("slowr",slow_diag)
2269	CALL xios_orchidee_send_field("streamr",stream_diag)
2270	CALL xios_orchidee_send_field("hydrographs",hydrographs_diag*dt_sechiba)
2271
2272	END SUBROUTINE routing_flow_diags
2273
2274
2275	!! =============================================================================================================================
2276	!! SUBROUTINE: routing_simple_finalize
2277	!!
2278	!>\BRIEF Write to restart file
2279	!!
2280	!! DESCRIPTION: Write module variables to restart file
2281	!!
2282	!! RECENT CHANGE(S)
2283	!!
2284	!! REFERENCE(S)
2285	!!
2286	!! FLOWCHART
2287	!! \n
2288	!_ ==============================================================================================================================
2289
2290	SUBROUTINE routing_flow_finalize(kjit, rest_id)
2291	USE xios
2292	USE ioipsl
2293	IMPLICIT NONE
2294	INTEGER, INTENT(IN) :: kjit
2295	INTEGER, INTENT(IN) :: rest_id
2296	!_ ================================================================================================================================
2297
2298	IF (is_omp_root) THEN
2299	CALL xios_send_field("fast_reservoir_restart",fast_reservoir_r)
2300	CALL xios_send_field("slow_reservoir_restart",slow_reservoir_r)
2301	CALL xios_send_field("stream_reservoir_restart",stream_reservoir_r)
2302	ENDIF
2303
2304	CALL restput_p (rest_id, 'riverflow', nbp_glo, 1, 1, kjit, riverflow_mean, 'scatter', nbp_glo, index_g)
2305	CALL restput_p (rest_id, 'coastalflow', nbp_glo, 1, 1, kjit, coastalflow_mean, 'scatter', nbp_glo, index_g)
2306	CALL restput_p (rest_id, 'lakeinflow', nbp_glo, 1, 1, kjit, lakeinflow_mean, 'scatter', nbp_glo, index_g)
2307	CALL routing_flow_finalize_mean(kjit, rest_id)
2308
2309	CALL routing_flow_finalize_diag(kjit, rest_id)
2310
2311	CALL xios_orchidee_change_context("orchidee_routing_out")
2312	CALL xios_context_finalize()
2313	CALL xios_orchidee_change_context("orchidee")
2314
2315	END SUBROUTINE routing_flow_finalize
2316
2317	SUBROUTINE routing_flow_finalize_diag(kjit, rest_id)
2318	USE ioipsl
2319	IMPLICIT NONE
2320	INTEGER, INTENT(IN) :: kjit
2321	INTEGER, INTENT(IN) :: rest_id
2322
2323	CALL restput_p (rest_id, 'hydrographs', nbp_glo, 1, 1, kjit, hydrographs_diag, 'scatter', nbp_glo, index_g)
2324
2325	END SUBROUTINE routing_flow_finalize_diag
2326
2327	!! ================================================================================================================================
2328	!! SUBROUTINE : routing_simple_clear
2329	!!
2330	!>\BRIEF This subroutine deallocates the block memory previously allocated.
2331	!!
2332	!! DESCRIPTION: This subroutine deallocates the block memory previously allocated.
2333	!!
2334	!! RECENT CHANGE(S): None
2335	!!
2336	!! MAIN OUTPUT VARIABLE(S):
2337	!!
2338	!! REFERENCES : None
2339	!!
2340	!! FLOWCHART :None
2341	!! \n
2342	!_ ================================================================================================================================
2343
2344	SUBROUTINE routing_flow_clear
2345	IMPLICIT NONE
2346
2347	IF (is_omp_root) THEN
2348	IF (ALLOCATED(topoind_r)) DEALLOCATE(topoind_r)
2349	IF (ALLOCATED(route_flow_rp1)) DEALLOCATE(route_flow_rp1)
2350	IF (ALLOCATED(routing_mask_r)) DEALLOCATE(routing_mask_r)
2351	IF (ALLOCATED(fast_reservoir_r)) DEALLOCATE(fast_reservoir_r)
2352	IF (ALLOCATED(slow_reservoir_r)) DEALLOCATE(slow_reservoir_r)
2353	IF (ALLOCATED(is_lakeinflow_r)) DEALLOCATE(is_lakeinflow_r)
2354	IF (ALLOCATED(is_coastalflow_r)) DEALLOCATE(is_coastalflow_r)
2355	IF (ALLOCATED(is_riverflow_r)) DEALLOCATE(is_riverflow_r)
2356	ENDIF
2357
2358	END SUBROUTINE routing_flow_clear
2359
2360
2361	END MODULE routing_native_flow_mod

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

Download in other formats: