Context Navigation

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

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

Update of routing native

More robust algorithm for correcting routing grid to adjust to coast line model, in global and LAM configuration
More robust water conservation
Reservoir constant are now given in km/day, independently of the routing grid topoindex unit. So a addtionnal parameter is added to the run.def to specify the type of routing file : routing_file_type = standard => 50km legacy routing file routing_file_type = merit => 2km Merit routing file

YM

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

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

Download in other formats: