1 | ! ================================================================================================================================= |
---|
2 | ! MODULE : routing_simple |
---|
3 | ! |
---|
4 | ! CONTACT : orchidee-help _at_ listes.ipsl.fr |
---|
5 | ! |
---|
6 | ! LICENCE : IPSL (2006) |
---|
7 | ! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC |
---|
8 | ! |
---|
9 | !>\BRIEF This module routes the water over the continents into the oceans and computes the water |
---|
10 | !! stored in floodplains or taken for irrigation. |
---|
11 | !! |
---|
12 | !!\n DESCRIPTION: The subroutines in this subroutine is only called when ROUTING_METHOD=simple is set in run.def. |
---|
13 | !! The method can be used for regular latitude-longitude grid or for unstructured grid. |
---|
14 | !! |
---|
15 | !! RECENT CHANGE(S): None |
---|
16 | !! |
---|
17 | !! REFERENCE(S) : |
---|
18 | !! |
---|
19 | !! SVN : |
---|
20 | !! $HeadURL: svn://forge.ipsl.jussieu.fr/orchidee/branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/routing_simple.f90 $ |
---|
21 | !! $Date: 2019-07-05 14:52:05 +0200 (Fri, 05 Jul 2019) $ |
---|
22 | !! $Revision: 6102 $ |
---|
23 | !! \n |
---|
24 | !_ ================================================================================================================================ |
---|
25 | |
---|
26 | |
---|
27 | |
---|
28 | ! Histoire Salee |
---|
29 | !--------------- |
---|
30 | ! La douce riviere |
---|
31 | ! Sortant de son lit |
---|
32 | ! S'est jetee ma chere |
---|
33 | ! dans les bras mais oui |
---|
34 | ! du beau fleuve |
---|
35 | ! |
---|
36 | ! L'eau coule sous les ponts |
---|
37 | ! Et puis les flots s'emeuvent |
---|
38 | ! - N'etes vous pas au courant ? |
---|
39 | ! Il parait que la riviere |
---|
40 | ! Va devenir mer |
---|
41 | ! Roland Bacri |
---|
42 | ! |
---|
43 | |
---|
44 | |
---|
45 | MODULE routing_simple |
---|
46 | |
---|
47 | USE ioipsl |
---|
48 | USE xios_orchidee |
---|
49 | USE ioipsl_para |
---|
50 | USE constantes |
---|
51 | USE constantes_soil |
---|
52 | USE pft_parameters |
---|
53 | USE sechiba_io_p |
---|
54 | USE interpol_help |
---|
55 | USE grid |
---|
56 | USE mod_orchidee_para |
---|
57 | |
---|
58 | |
---|
59 | IMPLICIT NONE |
---|
60 | PRIVATE |
---|
61 | PUBLIC :: routing_simple_main, routing_simple_xios_initialize |
---|
62 | PUBLIC :: routing_simple_initialize, routing_simple_finalize, routing_simple_clear |
---|
63 | |
---|
64 | !! PARAMETERS |
---|
65 | INTEGER(i_std), PARAMETER :: nbasmax=5 !! The maximum number of basins we wish to have per grid box (truncation of the model) (unitless) |
---|
66 | INTEGER(i_std), SAVE :: nbvmax !! The maximum number of basins we can handle at any time during the generation of the maps (unitless) |
---|
67 | !$OMP THREADPRIVATE(nbvmax) |
---|
68 | REAL(r_std), SAVE :: fast_tcst = 3.0 !! Property of the fast reservoir (day/m) |
---|
69 | !$OMP THREADPRIVATE(fast_tcst) |
---|
70 | REAL(r_std), SAVE :: slow_tcst = 25.0 !! Property of the slow reservoir (day/m) |
---|
71 | !$OMP THREADPRIVATE(slow_tcst) |
---|
72 | REAL(r_std), SAVE :: stream_tcst = 0.24 !! Property of the stream reservoir (day/m) |
---|
73 | !$OMP THREADPRIVATE(stream_tcst) |
---|
74 | REAL(r_std), SAVE :: flood_tcst = 4.0 !! Property of the floodplains reservoir (day/m) |
---|
75 | !$OMP THREADPRIVATE(flood_tcst) |
---|
76 | REAL(r_std), SAVE :: swamp_cst = 0.2 !! Fraction of the river transport that flows to the swamps (unitless;0-1) |
---|
77 | !$OMP THREADPRIVATE(swamp_cst) |
---|
78 | ! |
---|
79 | ! Relation between volume and fraction of floodplains |
---|
80 | ! |
---|
81 | REAL(r_std), SAVE :: beta = 2.0 !! Parameter to fix the shape of the floodplain (>1 for convex edges, <1 for concave edges) (unitless) |
---|
82 | !$OMP THREADPRIVATE(beta) |
---|
83 | REAL(r_std), SAVE :: betap = 0.5 !! Ratio of the basin surface intercepted by ponds and the maximum surface of ponds (unitless;0-1) |
---|
84 | !$OMP THREADPRIVATE(betap) |
---|
85 | REAL(r_std), SAVE :: floodcri = 2000.0 !! Potential height for which all the basin is flooded (mm) |
---|
86 | !$OMP THREADPRIVATE(floodcri) |
---|
87 | ! |
---|
88 | ! Relation between maximum surface of ponds and basin surface, and drainage (mm/j) to the slow_res |
---|
89 | ! |
---|
90 | REAL(r_std), PARAMETER :: pond_bas = 50.0 !! [DISPENSABLE] - not used |
---|
91 | REAL(r_std), SAVE :: pondcri = 2000.0 !! Potential height for which all the basin is a pond (mm) |
---|
92 | !$OMP THREADPRIVATE(pondcri) |
---|
93 | |
---|
94 | REAL(r_std), PARAMETER :: maxevap_lake = 7.5/86400. !! Maximum evaporation rate from lakes (kg/m^2/s) |
---|
95 | REAL(r_std),SAVE :: dt_routing !! Routing time step (s) |
---|
96 | !$OMP THREADPRIVATE(dt_routing) |
---|
97 | INTEGER(i_std), SAVE :: diagunit = 87 !! Diagnostic file unit (unitless) |
---|
98 | !$OMP THREADPRIVATE(diagunit) |
---|
99 | |
---|
100 | ! Logicals to control model configuration |
---|
101 | ! |
---|
102 | LOGICAL, SAVE :: dofloodinfilt = .FALSE. !! Logical to choose if floodplains infiltration is activated or not (true/false) |
---|
103 | !$OMP THREADPRIVATE(dofloodinfilt) |
---|
104 | LOGICAL, SAVE :: doswamps = .FALSE. !! Logical to choose if swamps are activated or not (true/false) |
---|
105 | !$OMP THREADPRIVATE(doswamps) |
---|
106 | LOGICAL, SAVE :: doponds = .FALSE. !! Logical to choose if ponds are activated or not (true/false) |
---|
107 | !$OMP THREADPRIVATE(doponds) |
---|
108 | |
---|
109 | INTEGER(i_std), SAVE :: num_largest !! Number of largest river basins which should be treated as independently as rivers |
---|
110 | !! (not flow into ocean as diffusion coastal flow) (unitless) |
---|
111 | !$OMP THREADPRIVATE(num_largest) |
---|
112 | REAL(r_std), SAVE :: time_counter !! Time counter (s) |
---|
113 | !$OMP THREADPRIVATE(time_counter) |
---|
114 | REAL(r_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: routing_area_loc !! Surface of basin (m^2) |
---|
115 | !$OMP THREADPRIVATE(routing_area_loc) |
---|
116 | REAL(r_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: topo_resid_loc !! Topographic index of the retention time (m) |
---|
117 | !$OMP THREADPRIVATE(topo_resid_loc) |
---|
118 | INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: route_togrid_loc !! Grid into which the basin flows (unitless) |
---|
119 | !$OMP THREADPRIVATE(route_togrid_loc) |
---|
120 | INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: route_tobasin_loc !! Basin in to which the water goes (unitless) |
---|
121 | !$OMP THREADPRIVATE(route_tobasin_loc) |
---|
122 | INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: route_nbintobas_loc !! Number of basin into current one (unitless) |
---|
123 | !$OMP THREADPRIVATE(route_nbintobas_loc) |
---|
124 | INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: global_basinid_loc !! ID of basin (unitless) |
---|
125 | !$OMP THREADPRIVATE(global_basinid_loc) |
---|
126 | INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: hydrodiag_loc !! Variable to diagnose the hydrographs |
---|
127 | !$OMP THREADPRIVATE(hydrodiag_loc) |
---|
128 | REAL(r_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:) :: hydroupbasin_loc !! The area upstream of the gauging station (m^2) |
---|
129 | !$OMP THREADPRIVATE(hydroupbasin_loc) |
---|
130 | |
---|
131 | ! parallelism |
---|
132 | REAL(r_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: routing_area_glo !! Surface of basin (m^2) |
---|
133 | !$OMP THREADPRIVATE(routing_area_glo) |
---|
134 | REAL(r_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: topo_resid_glo !! Topographic index of the retention time (m) |
---|
135 | !$OMP THREADPRIVATE(topo_resid_glo) |
---|
136 | INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: route_togrid_glo !! Grid into which the basin flows (unitless) |
---|
137 | !$OMP THREADPRIVATE(route_togrid_glo) |
---|
138 | INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: route_tobasin_glo !! Basin in to which the water goes (unitless) |
---|
139 | !$OMP THREADPRIVATE(route_tobasin_glo) |
---|
140 | INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: route_nbintobas_glo !! Number of basin into current one (unitless) |
---|
141 | !$OMP THREADPRIVATE(route_nbintobas_glo) |
---|
142 | INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: global_basinid_glo !! ID of basin (unitless) |
---|
143 | !$OMP THREADPRIVATE(global_basinid_glo) |
---|
144 | INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: hydrodiag_glo !! Variable to diagnose the hydrographs |
---|
145 | !$OMP THREADPRIVATE(hydrodiag_glo) |
---|
146 | REAL(r_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:) :: hydroupbasin_glo !! The area upstream of the gauging station (m^2) |
---|
147 | !$OMP THREADPRIVATE(hydroupbasin_glo) |
---|
148 | |
---|
149 | REAL(r_std), SAVE, POINTER, DIMENSION(:,:) :: routing_area !! Surface of basin (m^2) |
---|
150 | !$OMP THREADPRIVATE(routing_area) |
---|
151 | REAL(r_std), SAVE, POINTER, DIMENSION(:,:) :: topo_resid !! Topographic index of the retention time (m) |
---|
152 | !$OMP THREADPRIVATE(topo_resid) |
---|
153 | INTEGER(i_std), SAVE, POINTER, DIMENSION(:,:) :: route_togrid !! Grid into which the basin flows (unitless) |
---|
154 | !$OMP THREADPRIVATE(route_togrid) |
---|
155 | INTEGER(i_std), SAVE, POINTER, DIMENSION(:,:) :: route_tobasin !! Basin in to which the water goes (unitless) |
---|
156 | !$OMP THREADPRIVATE(route_tobasin) |
---|
157 | INTEGER(i_std), SAVE, POINTER, DIMENSION(:,:) :: route_nbintobas !! Number of basin into current one (unitless) |
---|
158 | !$OMP THREADPRIVATE(route_nbintobas) |
---|
159 | INTEGER(i_std), SAVE, POINTER, DIMENSION(:,:) :: global_basinid !! ID of basin (unitless) |
---|
160 | !$OMP THREADPRIVATE(global_basinid) |
---|
161 | INTEGER(i_std), SAVE, POINTER, DIMENSION(:,:) :: hydrodiag !! Variable to diagnose the hydrographs |
---|
162 | !$OMP THREADPRIVATE(hydrodiag) |
---|
163 | REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: slowflow_diag !! Diagnostic slow flow hydrographs (kg/dt) |
---|
164 | !$OMP THREADPRIVATE(slowflow_diag) |
---|
165 | REAL(r_std), SAVE, POINTER, DIMENSION(:) :: hydroupbasin !! The area upstream of the gauging station (m^2) |
---|
166 | !$OMP THREADPRIVATE(hydroupbasin) |
---|
167 | |
---|
168 | REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: irrigated !! Area equipped for irrigation in each grid box (m^2) |
---|
169 | !$OMP THREADPRIVATE(irrigated) |
---|
170 | REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: floodplains !! Maximal surface which can be inundated in each grid box (m^2) |
---|
171 | !$OMP THREADPRIVATE(floodplains) |
---|
172 | REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: swamp !! Maximal surface of swamps in each grid box (m^2) |
---|
173 | !$OMP THREADPRIVATE(swamp) |
---|
174 | ! |
---|
175 | ! The reservoirs, variables kept in the restart file. |
---|
176 | ! |
---|
177 | REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: fast_reservoir !! Water amount in the fast reservoir (kg) |
---|
178 | !$OMP THREADPRIVATE(fast_reservoir) |
---|
179 | REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: slow_reservoir !! Water amount in the slow reservoir (kg) |
---|
180 | !$OMP THREADPRIVATE(slow_reservoir) |
---|
181 | REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: stream_reservoir !! Water amount in the stream reservoir (kg) |
---|
182 | !$OMP THREADPRIVATE(stream_reservoir) |
---|
183 | REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: flood_reservoir !! Water amount in the floodplains reservoir (kg) |
---|
184 | !$OMP THREADPRIVATE(flood_reservoir) |
---|
185 | REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: lake_reservoir !! Water amount in the lake reservoir (kg) |
---|
186 | !$OMP THREADPRIVATE(lake_reservoir) |
---|
187 | REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: pond_reservoir !! Water amount in the pond reservoir (kg) |
---|
188 | !$OMP THREADPRIVATE(pond_reservoir) |
---|
189 | REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: flood_frac_bas !! Flooded fraction per basin (unitless;0-1) |
---|
190 | !$OMP THREADPRIVATE(flood_frac_bas) |
---|
191 | REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: pond_frac !! Pond fraction per grid box (unitless;0-1) |
---|
192 | !$OMP THREADPRIVATE(pond_frac) |
---|
193 | REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: flood_height !! Floodplain height (mm) |
---|
194 | !$OMP THREADPRIVATE(flood_height) |
---|
195 | ! |
---|
196 | ! The accumulated fluxes. |
---|
197 | ! |
---|
198 | REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: floodout_mean !! Accumulated flow out of floodplains (kg/m^2/dt) |
---|
199 | !$OMP THREADPRIVATE(floodout_mean) |
---|
200 | REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: runoff_mean !! Accumulated runoff (kg/m^2/dt) |
---|
201 | !$OMP THREADPRIVATE(runoff_mean) |
---|
202 | REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: drainage_mean !! Accumulated drainage (kg/m^2/dt) |
---|
203 | !$OMP THREADPRIVATE(drainage_mean) |
---|
204 | REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: transpot_mean !! Mean potential transpiration from the plants (kg/m^2/dt) |
---|
205 | !$OMP THREADPRIVATE(transpot_mean) |
---|
206 | REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: precip_mean !! Accumulated precipitation (kg/m^2/dt) |
---|
207 | !$OMP THREADPRIVATE(precip_mean) |
---|
208 | REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: humrel_mean !! Mean soil moisture stress, mean root extraction potential (unitless) |
---|
209 | !$OMP THREADPRIVATE(humrel_mean) |
---|
210 | REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: totnobio_mean !! Mean last total fraction of no bio (unitless;0-1) |
---|
211 | !$OMP THREADPRIVATE(totnobio_mean) |
---|
212 | REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: vegtot_mean !! Mean potentially vegetated fraction (unitless;0-1) |
---|
213 | !$OMP THREADPRIVATE(vegtot_mean) |
---|
214 | REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: k_litt_mean !! Mean averaged conductivity for saturated infiltration in the 'litter' layer (kg/m^2/dt) |
---|
215 | !$OMP THREADPRIVATE(k_litt_mean) |
---|
216 | ! |
---|
217 | ! The averaged outflow fluxes. |
---|
218 | ! |
---|
219 | REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: lakeinflow_mean !! Mean lake inflow (kg/m^2/dt) |
---|
220 | !$OMP THREADPRIVATE(lakeinflow_mean) |
---|
221 | REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: returnflow_mean !! Mean water flow from lakes and swamps which returns to the grid box. |
---|
222 | !! This water will go back into the hydrol or hydrolc module to allow re-evaporation (kg/m^2/dt) |
---|
223 | !$OMP THREADPRIVATE(returnflow_mean) |
---|
224 | REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: reinfiltration_mean !! Mean water flow which returns to the grid box (kg/m^2/dt) |
---|
225 | !$OMP THREADPRIVATE(reinfiltration_mean) |
---|
226 | REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: irrigation_mean !! Mean irrigation flux. |
---|
227 | !! This is the water taken from the reservoirs and beeing put into the upper layers of the soil (kg/m^2/dt) |
---|
228 | !$OMP THREADPRIVATE(irrigation_mean) |
---|
229 | REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: riverflow_mean !! Mean Outflow of the major rivers. |
---|
230 | !! The flux will be located on the continental grid but this should be a coastal point (kg/dt) |
---|
231 | !$OMP THREADPRIVATE(riverflow_mean) |
---|
232 | REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: coastalflow_mean !! Mean outflow on coastal points by small basins. |
---|
233 | !! This is the water which flows in a disperse way into the ocean (kg/dt) |
---|
234 | !$OMP THREADPRIVATE(coastalflow_mean) |
---|
235 | REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: floodtemp !! Temperature to decide if floodplains work (K) |
---|
236 | !$OMP THREADPRIVATE(floodtemp) |
---|
237 | INTEGER(i_std), SAVE :: floodtemp_lev !! Temperature level to decide if floodplains work (K) |
---|
238 | !$OMP THREADPRIVATE(floodtemp_lev) |
---|
239 | ! |
---|
240 | ! Diagnostic variables ... well sort of ! |
---|
241 | ! |
---|
242 | REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: irrig_netereq !! Irrigation requirement (water requirements by the crop for its optimal growth (kg/m^2/dt) |
---|
243 | !$OMP THREADPRIVATE(irrig_netereq) |
---|
244 | REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: hydrographs !! Hydrographs at the outflow of the grid box for major basins (kg/dt) |
---|
245 | !$OMP THREADPRIVATE(hydrographs) |
---|
246 | ! |
---|
247 | ! Diagnostics for the various reservoirs we use (Kg/m^2) |
---|
248 | ! |
---|
249 | REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: fast_diag !! Diagnostic for the fast reservoir (kg/m^2) |
---|
250 | !$OMP THREADPRIVATE(fast_diag) |
---|
251 | REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: slow_diag !! Diagnostic for the slow reservoir (kg/m^2) |
---|
252 | !$OMP THREADPRIVATE(slow_diag) |
---|
253 | REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: stream_diag !! Diagnostic for the stream reservoir (kg/m^2) |
---|
254 | !$OMP THREADPRIVATE(stream_diag) |
---|
255 | REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: flood_diag !! Diagnostic for the floodplain reservoir (kg/m^2) |
---|
256 | !$OMP THREADPRIVATE(flood_diag) |
---|
257 | REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: pond_diag !! Diagnostic for the pond reservoir (kg/m^2) |
---|
258 | !$OMP THREADPRIVATE(pond_diag) |
---|
259 | REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: lake_diag !! Diagnostic for the lake reservoir (kg/m^2) |
---|
260 | !$OMP THREADPRIVATE(lake_diag) |
---|
261 | REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: delsurfstor !! Diagnostic of the change in surface water storage (flood, pond and lake reservoirs) (kg/m^2) |
---|
262 | !$OMP THREADPRIVATE(delsurfstor) |
---|
263 | |
---|
264 | ! |
---|
265 | ! Specific variables for simple routing |
---|
266 | ! |
---|
267 | REAL(r_std),SAVE,ALLOCATABLE :: topoind_r(:) !! Topographic index of the retention time (m) index - (local routing grid) |
---|
268 | !$OMP THREADPRIVATE(topoind_r) |
---|
269 | INTEGER,SAVE,ALLOCATABLE :: route_flow_rp1(:) !! flow index from cell to neighboring cell following the trip direction - (local routing grid + halo) |
---|
270 | !$OMP THREADPRIVATE(route_flow_rp1) |
---|
271 | REAL(r_std),SAVE,ALLOCATABLE :: fast_reservoir_r(:) !! Water amount in the fast reservoir (kg) - (local routing grid) |
---|
272 | !$OMP THREADPRIVATE(fast_reservoir_r) |
---|
273 | REAL(r_std),SAVE,ALLOCATABLE :: slow_reservoir_r(:) !! Water amount in the slow reservoir (kg) - (local routing grid) |
---|
274 | !$OMP THREADPRIVATE(slow_reservoir_r) |
---|
275 | REAL(r_std),SAVE,ALLOCATABLE :: stream_reservoir_r(:) !! Water amount in the stream reservoir (kg) - (local routing grid) |
---|
276 | !$OMP THREADPRIVATE(stream_reservoir_r) |
---|
277 | LOGICAL,SAVE,ALLOCATABLE :: is_lakeinflow_r(:) !! is lake inflow point - (local routing grid) |
---|
278 | !$OMP THREADPRIVATE(is_lakeinflow_r) |
---|
279 | LOGICAL,SAVE,ALLOCATABLE :: is_coastalflow_r(:) !! is coastal flow point - (local routing grid) |
---|
280 | !$OMP THREADPRIVATE(is_coastalflow_r) |
---|
281 | LOGICAL,SAVE,ALLOCATABLE :: is_riverflow_r(:) !! is river flow point - (local routing grid) |
---|
282 | !$OMP THREADPRIVATE(is_riverflow_r) |
---|
283 | LOGICAL,SAVE,ALLOCATABLE :: is_streamflow_r(:) !! is stream flow point - (local routing grid) |
---|
284 | !$OMP THREADPRIVATE(is_streamflow_r) |
---|
285 | LOGICAL,SAVE,ALLOCATABLE :: routing_mask_r(:) !! valid routing point - (local routing grid) |
---|
286 | !$OMP THREADPRIVATE(routing_mask_r) |
---|
287 | LOGICAL,SAVE,ALLOCATABLE :: coast_mask(:) !! is a coast point - (local native grid) |
---|
288 | !$OMP THREADPRIVATE(coast_mask) |
---|
289 | INTEGER,SAVE :: nbpt_r !! number of point in local routing grid |
---|
290 | !$OMP THREADPRIVATE(nbpt_r) |
---|
291 | INTEGER,SAVE :: nbpt_rp1 !! number of point in local routing grid with halo of 1 |
---|
292 | !$OMP THREADPRIVATE(nbpt_rp1) |
---|
293 | 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) |
---|
294 | !$OMP THREADPRIVATE(routing_weight) |
---|
295 | REAL(r_std),SAVE,ALLOCATABLE :: routing_weight_r(:) !! Weight to transfer lost water into intersecting cell when doing interpolation from local routing grid to local native grid (lakeinflow) |
---|
296 | !$OMP THREADPRIVATE(routing_weight_r) |
---|
297 | REAL(r_std),SAVE,ALLOCATABLE :: routing_weight_coast_r(:) !! Weight to transfer lost water into intersecting cell on coast |
---|
298 | !$OMP THREADPRIVATE(routing_weight_coast_r) |
---|
299 | !! when doing interpolation from local routing grid to local native grid (river flow+coastal flow) |
---|
300 | REAL(r_std),SAVE,ALLOCATABLE :: basins_extended_r(:) !! basins riverflow id (local routing grid) |
---|
301 | !$OMP THREADPRIVATE(basins_extended_r) |
---|
302 | INTEGER(i_std) :: basins_count !! number of basins (local routing grid) |
---|
303 | !$OMP THREADPRIVATE(basins_count) |
---|
304 | INTEGER(i_std) :: basins_out !! number of basins to output for diag |
---|
305 | !$OMP THREADPRIVATE(basins_out) |
---|
306 | |
---|
307 | INTEGER(i_std) :: split_routing !! time spliting for routing |
---|
308 | !$OMP THREADPRIVATE(split_routing) |
---|
309 | |
---|
310 | |
---|
311 | INTEGER(i_std), PARAMETER :: nb_stations=14 |
---|
312 | REAL(r_std),PARAMETER :: station_lon(nb_stations) = & |
---|
313 | (/ -162.8830, -90.9058, -55.5110, -49.3242, -133.7447, -63.6000, 28.7167, & |
---|
314 | 15.3000, 66.5300, 89.6700, 86.5000, 127.6500, 3.3833, 117.6200 /) |
---|
315 | REAL(r_std),PARAMETER :: station_lat(nb_stations) = & |
---|
316 | (/ 61.9340, 32.3150, -1.9470, -5.1281, 67.4583, 8.1500, 45.2167, & |
---|
317 | -4.3000, 66.5700, 25.1800, 67.4800, 70.7000, 11.8667, 30.7700 /) |
---|
318 | CHARACTER(LEN=17),PARAMETER :: station_name(nb_stations) = & |
---|
319 | (/ "Pilot station ", "Vicksburg ", "Obidos ", & |
---|
320 | "Itupiranga ", "Arctic red river ", "Puente Angostura ", & |
---|
321 | "Ceatal Izmail ", "Kinshasa ", "Salekhard ", & |
---|
322 | "Bahadurabad ", "Igarka ", "Kusur ", & |
---|
323 | "Malanville ", "Datong " /) |
---|
324 | |
---|
325 | CONTAINS |
---|
326 | |
---|
327 | !! ============================================================================================================================= |
---|
328 | !! SUBROUTINE: routing_simple_xios_initialize |
---|
329 | !! |
---|
330 | !>\BRIEF Initialize xios dependant defintion before closing context defintion |
---|
331 | !! |
---|
332 | !! DESCRIPTION: Initialize xios dependant defintion before closing context defintion. |
---|
333 | !! This subroutine is called before the xios context is closed. |
---|
334 | !! |
---|
335 | !! RECENT CHANGE(S): None |
---|
336 | !! |
---|
337 | !! REFERENCE(S): None |
---|
338 | !! |
---|
339 | !! FLOWCHART: None |
---|
340 | !! \n |
---|
341 | !_ ============================================================================================================================== |
---|
342 | |
---|
343 | SUBROUTINE routing_simple_xios_initialize |
---|
344 | USE xios |
---|
345 | USE routing, ONLY : routing_names |
---|
346 | IMPLICIT NONE |
---|
347 | |
---|
348 | !! 0 Variable and parameter description |
---|
349 | CHARACTER(LEN=60),ALLOCATABLE :: label(:) |
---|
350 | LOGICAL :: file_exists |
---|
351 | |
---|
352 | IF (is_omp_root) THEN |
---|
353 | CALL xios_get_axis_attr("basins", n_glo=basins_out) ! get nb basins to output |
---|
354 | ALLOCATE(label(basins_out)) |
---|
355 | CALL routing_names(basins_out,label) |
---|
356 | CALL xios_set_axis_attr("basins", label=label) ! set riverflow basins name |
---|
357 | INQUIRE(FILE="routing_start.nc", EXIST=file_exists) |
---|
358 | IF (file_exists) CALL xios_set_file_attr("routing_start", enabled=.TRUE.) |
---|
359 | ENDIF |
---|
360 | |
---|
361 | END SUBROUTINE routing_simple_xios_initialize |
---|
362 | |
---|
363 | |
---|
364 | |
---|
365 | !! ============================================================================================================================= |
---|
366 | !! SUBROUTINE: routing_simple_init_2 |
---|
367 | !! |
---|
368 | !>\BRIEF Privat subroutine |
---|
369 | !! |
---|
370 | !! DESCRIPTION: Privat subroutine to the module routing_simple. This subroutine is called in the end of routing_simple_initialize |
---|
371 | !! |
---|
372 | !! RECENT CHANGE(S): None |
---|
373 | !! |
---|
374 | !! REFERENCE(S): None |
---|
375 | !! |
---|
376 | !! FLOWCHART: None |
---|
377 | !! \n |
---|
378 | !_ ============================================================================================================================== |
---|
379 | |
---|
380 | SUBROUTINE routing_simple_init_2(nbpt, contfrac) |
---|
381 | USE xios |
---|
382 | USE grid, ONLY : area |
---|
383 | IMPLICIT NONE |
---|
384 | INCLUDE "mpif.h" |
---|
385 | |
---|
386 | !! 0 Variable and parameter description |
---|
387 | !! 0.1 Input variables |
---|
388 | INTEGER,INTENT(in) :: nbpt !! nb points native grid |
---|
389 | REAL,INTENT(in) :: contfrac(nbpt) !! fraction of land |
---|
390 | |
---|
391 | !! 0.2 Local variables |
---|
392 | INTEGER :: ni !! longitude dimension of local routing grid |
---|
393 | INTEGER :: nj !! latitude dimension of local routing grid |
---|
394 | |
---|
395 | 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) |
---|
396 | 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) |
---|
397 | !! routing is artificially computed on sea and endoric basins |
---|
398 | LOGICAL,ALLOCATABLE :: routing_mask_rp1(:,:) !! valid routing point - local routing grid+halo (0:ni+1,0:nj+1) |
---|
399 | REAL(r_std),ALLOCATABLE :: mask_native(:) !! some mask on native grid (nbpt) |
---|
400 | REAL(r_std),ALLOCATABLE :: frac_routing_r(:) !! fraction of routing cell intesected by cells of native grid |
---|
401 | REAL(r_std),ALLOCATABLE :: frac_routing_coast_r(:) !! fraction of routing cell intesected by coastal cells of native grid |
---|
402 | INTEGER :: ij, ij_r, ij_rp1, i,j,jp1,jm1,ip1,im1,jr,ir |
---|
403 | REAL(r_std) :: sum1,sum2,sum_frac_routing_r |
---|
404 | REAL(r_std) :: contfrac_mpi(nbp_mpi) |
---|
405 | REAL(r_std) :: area_mpi(nbp_mpi) |
---|
406 | INTEGER :: basins_count_mpi |
---|
407 | INTEGER :: ierr |
---|
408 | LOGICAL :: file_exists |
---|
409 | |
---|
410 | !_ ================================================================================================================================ |
---|
411 | |
---|
412 | split_routing=1 |
---|
413 | CALL getin_p("SPLIT_ROUTING",split_routing) |
---|
414 | |
---|
415 | CALL gather_omp(contfrac,contfrac_mpi) |
---|
416 | CALL gather_omp(area, area_mpi) |
---|
417 | |
---|
418 | IF (is_omp_root) THEN |
---|
419 | |
---|
420 | CALL xios_get_domain_attr("routing_domain", ni=ni, nj=nj) ! get routing domain dimension |
---|
421 | |
---|
422 | nbpt_r= ni*nj |
---|
423 | nbpt_rp1= (ni+2)*(nj+2) |
---|
424 | |
---|
425 | ALLOCATE(trip_extended_r(ni,nj)) |
---|
426 | ALLOCATE(topoind_r(ni*nj)) |
---|
427 | ALLOCATE(is_lakeinflow_r(ni*nj)) |
---|
428 | ALLOCATE(is_coastalflow_r(ni*nj)) |
---|
429 | ALLOCATE(is_riverflow_r(ni*nj)) |
---|
430 | ALLOCATE(is_streamflow_r(ni*nj)) |
---|
431 | ALLOCATE(routing_mask_r(ni*nj)) |
---|
432 | ALLOCATE(fast_reservoir_r(ni*nj)) |
---|
433 | ALLOCATE(slow_reservoir_r(ni*nj)) |
---|
434 | ALLOCATE(stream_reservoir_r(ni*nj)) |
---|
435 | |
---|
436 | ALLOCATE(mask_native(nbp_mpi)) |
---|
437 | ALLOCATE(coast_mask(nbp_mpi)) |
---|
438 | ALLOCATE(frac_routing_r(nbpt_r)) |
---|
439 | ALLOCATE(frac_routing_coast_r(nbpt_r)) |
---|
440 | ALLOCATE(routing_weight(nbp_mpi)) |
---|
441 | ALLOCATE(routing_weight_r(nbpt_r)) |
---|
442 | ALLOCATE(routing_weight_coast_r(nbpt_r)) |
---|
443 | ALLOCATE(basins_extended_r(nbpt_r)) |
---|
444 | |
---|
445 | coast_mask=.FALSE. |
---|
446 | WHERE( contfrac_mpi(:)< 1.-1.e-5) coast_mask(:)=.TRUE. ! create mask for coastal cells on native grid |
---|
447 | |
---|
448 | INQUIRE(FILE="routing_start.nc", EXIST=file_exists) |
---|
449 | |
---|
450 | IF (file_exists) THEN |
---|
451 | CALL xios_recv_field("fast_reservoir_start",fast_reservoir_r) |
---|
452 | CALL xios_recv_field("slow_reservoir_start",slow_reservoir_r) |
---|
453 | CALL xios_recv_field("stream_reservoir_start",stream_reservoir_r) |
---|
454 | ELSE |
---|
455 | fast_reservoir_r(:)=0 |
---|
456 | slow_reservoir_r(:)=0 |
---|
457 | stream_reservoir_r(:)=0 |
---|
458 | ENDIF |
---|
459 | |
---|
460 | ALLOCATE(trip_rp1(0:ni+1,0:nj+1)) |
---|
461 | trip_rp1(:,:)=1e10 |
---|
462 | ALLOCATE(routing_mask_rp1(0:ni+1,0:nj+1)) |
---|
463 | ALLOCATE(route_flow_rp1((ni+2)*(nj+2))) |
---|
464 | |
---|
465 | CALL xios_send_field("routing_contfrac",contfrac_mpi) ! diag |
---|
466 | CALL xios_recv_field("trip_r",trip_rp1(1:ni,1:nj)) ! recv trip array with halo of 1 |
---|
467 | CALL xios_recv_field("trip_extended_r",trip_extended_r) ! recv extended trip array from file |
---|
468 | CALL xios_recv_field("topoind_r",topoind_r) ! recv topo index array from file |
---|
469 | CALL xios_recv_field("basins_extended_r",basins_extended_r) ! recv basins index from file |
---|
470 | |
---|
471 | is_lakeinflow_r(:) = .FALSE. |
---|
472 | is_coastalflow_r(:) = .FALSE. |
---|
473 | is_riverflow_r(:) = .FALSE. |
---|
474 | |
---|
475 | route_flow_rp1(:)=-1 |
---|
476 | |
---|
477 | routing_mask_rp1=.TRUE. |
---|
478 | DO j=1,nj |
---|
479 | DO i=1,ni |
---|
480 | IF (trip_rp1(i,j)>99) routing_mask_rp1(i,j)=.FALSE. ! if ocean or endoric basins => masked point |
---|
481 | ENDDO |
---|
482 | ENDDO |
---|
483 | |
---|
484 | mask_native(:)=0 |
---|
485 | frac_routing_r(:)=0 |
---|
486 | WHERE( .NOT. coast_mask(:)) mask_native(:)=1 |
---|
487 | CALL xios_send_field("mask_native_lake",mask_native) ! send full land point to XIOS (native grid) |
---|
488 | CALL xios_recv_field("frac_routing_lake_r",frac_routing_r) ! receive fraction of intersected cell by full land, on routing grid |
---|
489 | |
---|
490 | |
---|
491 | DO j=1,nj |
---|
492 | DO i=1,ni |
---|
493 | ij=i+ni*(j-1) |
---|
494 | IF (frac_routing_r(ij)>1-1e-5) THEN ! if full land point (no part of coast) |
---|
495 | IF ((.NOT. routing_mask_rp1(i,j)) .OR. trip_rp1(i,j)==98 .OR. trip_rp1(i,j)==99) THEN ! if masked point or river flow or coastal flow |
---|
496 | routing_mask_rp1(i,j)=.TRUE. ! transform to non masked |
---|
497 | trip_rp1(i,j)=trip_extended_r(i,j) ! replace by extended trip value |
---|
498 | ENDIF |
---|
499 | ENDIF |
---|
500 | ENDDO |
---|
501 | ENDDO |
---|
502 | |
---|
503 | mask_native(:)=0 |
---|
504 | frac_routing_coast_r(:)=0 |
---|
505 | WHERE( coast_mask(:)) mask_native(:)=1 |
---|
506 | CALL xios_send_field("mask_native_coast",mask_native) ! send only coastal point to XIOS (native grid) |
---|
507 | CALL xios_recv_field("frac_routing_coast_r",frac_routing_coast_r) ! receive fraction of intersected cell by coastal point, on routing grid |
---|
508 | |
---|
509 | DO j=1,nj |
---|
510 | DO i=1,ni |
---|
511 | ij=i+ni*(j-1) |
---|
512 | IF (frac_routing_coast_r(ij) > 0) THEN ! fraction of coastal point ? |
---|
513 | IF (.NOT. routing_mask_rp1(i,j)) THEN ! if masked, transform to coastal flow to conserve the water balance |
---|
514 | routing_mask_rp1(i,j)=.TRUE. |
---|
515 | trip_rp1(i,j)=98 ! transform to coastal flow |
---|
516 | topoind_r(ij)=1e-10 ! no residence time |
---|
517 | ENDIF |
---|
518 | ENDIF |
---|
519 | ENDDO |
---|
520 | ENDDO |
---|
521 | |
---|
522 | |
---|
523 | mask_native(:)=1 |
---|
524 | frac_routing_r(:)=0 |
---|
525 | CALL xios_send_field("mask_native",mask_native) ! send only all point to XIOS (native grid) |
---|
526 | CALL xios_recv_field("frac_routing_r",frac_routing_r) ! receive fraction of intersected cell, (routing grid) |
---|
527 | DO j=1,nj |
---|
528 | DO i=1,ni |
---|
529 | ij=i+ni*(j-1) |
---|
530 | IF (frac_routing_r(ij)>0) THEN |
---|
531 | IF (.NOT. routing_mask_rp1(i,j)) THEN |
---|
532 | routing_mask_rp1(i,j)=.TRUE. ! may never happen |
---|
533 | ENDIF |
---|
534 | ELSE |
---|
535 | routing_mask_rp1(i,j)=.FALSE. ! if no intersection from native cells, mask the point (no routage) |
---|
536 | trip_rp1(i,j)=1e10 |
---|
537 | ENDIF |
---|
538 | ENDDO |
---|
539 | ENDDO |
---|
540 | |
---|
541 | CALL xios_send_field("trip_update_r",trip_rp1(1:ni,1:nj)) ! send to xios trip array to update for halo |
---|
542 | CALL xios_recv_field("trip_rp1",trip_rp1) ! recv trip array with halo of 1 |
---|
543 | |
---|
544 | |
---|
545 | routing_mask_rp1=.TRUE. |
---|
546 | DO j=0,nj+1 |
---|
547 | DO i=0,ni+1 |
---|
548 | IF (trip_rp1(i,j)>99) routing_mask_rp1(i,j)=.FALSE. ! update mask with halo of 1 |
---|
549 | ENDDO |
---|
550 | ENDDO |
---|
551 | |
---|
552 | !! Compute the routing |
---|
553 | !! Loop on all point point of the local routing grid + halo |
---|
554 | DO j=0,nj+1 |
---|
555 | jp1=j+1 |
---|
556 | jm1=j-1 |
---|
557 | DO i=0,ni+1 |
---|
558 | ij_rp1=i+(ni+2)*j+1 |
---|
559 | ij_r=i+ni*(j-1) |
---|
560 | ip1=i+1 |
---|
561 | im1=i-1 |
---|
562 | |
---|
563 | IF (trip_rp1(i,j) < 100) THEN |
---|
564 | |
---|
565 | ir=-1 ; jr=-1 ! -> -1 for 97,98,99 |
---|
566 | SELECT CASE (NINT(trip_rp1(i,j))) ! get the trip value for each points |
---|
567 | CASE (1) |
---|
568 | jr=jm1 ; ir=i ! north |
---|
569 | CASE (2) |
---|
570 | jr=jm1 ; ir=ip1 ! north-east |
---|
571 | CASE (3) |
---|
572 | jr=j ; ir=ip1 ! east |
---|
573 | CASE (4) |
---|
574 | jr=jp1 ; ir=ip1 ! south-east |
---|
575 | CASE (5) |
---|
576 | jr=jp1 ; ir=i ! south |
---|
577 | CASE(6) |
---|
578 | jr=jp1 ; ir=im1 ! south-west |
---|
579 | CASE (7) |
---|
580 | jr=j ; ir=im1 ! west |
---|
581 | CASE (8) |
---|
582 | jr=jm1 ; ir=im1 ! north-west |
---|
583 | CASE (97) |
---|
584 | IF ( i>0 .AND. i<ni+1 .AND. j>0 .AND. j<nj+1) THEN ! if inside my local domain |
---|
585 | is_lakeinflow_r(ij_r)=.TRUE. ! I am a lakeinflow point and route to myself |
---|
586 | jr=j ; ir=i |
---|
587 | ENDIF |
---|
588 | CASE (98) |
---|
589 | IF ( i>0 .AND. i<ni+1 .AND. j>0 .AND. j<nj+1) THEN ! if inside my local domain |
---|
590 | is_coastalflow_r(ij_r)=.TRUE. ! I am a coastal flow point and route to myself |
---|
591 | jr=j ; ir=i |
---|
592 | ENDIF |
---|
593 | CASE (99) |
---|
594 | IF ( i>0 .AND. i<ni+1 .AND. j>0 .AND. j<nj+1) THEN ! if inside my local domain |
---|
595 | is_riverflow_r(ij_r)=.TRUE. ! I am a riverflow point and route to myself |
---|
596 | jr=j ; ir=i |
---|
597 | ENDIF |
---|
598 | END SELECT |
---|
599 | |
---|
600 | IF (ir<0 .OR. ir>ni+1 .OR. jr<0 .OR. jr>nj+1) THEN |
---|
601 | route_flow_rp1(ij_rp1)=-1 ! if route outside my local domain+halo, no routing (will be done by other process) |
---|
602 | ELSE |
---|
603 | |
---|
604 | IF ( .NOT. routing_mask_rp1(ir,jr)) THEN ! if routing to a masked point, cut the flow |
---|
605 | jr=j ; ir=i ; ! and route to myself |
---|
606 | IF (i>0 .AND. i<ni+1 .AND. j>0 .AND. j<nj+1) THEN ! if some flow come from neighbour cell |
---|
607 | IF ( trip_rp1(i-1,j-1)==4 .OR. trip_rp1(i,j-1)==5 .OR. trip_rp1(i+1,j)==6 & |
---|
608 | .OR. trip_rp1(i-1,j)==3 .OR. trip_rp1(i+1,j)==7 & |
---|
609 | .OR. trip_rp1(i-1,j+1)==2 .OR. trip_rp1(i-1,j+1)==1 .OR. trip_rp1(i+1,j+1)==8) THEN |
---|
610 | is_riverflow_r(ij_r)=.TRUE. ! => transform to riverflow |
---|
611 | ELSE |
---|
612 | is_coastalflow_r(ij_r)=.TRUE. ! else transform to coastalflow |
---|
613 | ENDIF |
---|
614 | ENDIF |
---|
615 | ENDIF |
---|
616 | |
---|
617 | IF (ir<1 .OR. ir>ni .OR. jr<1 .OR. jr>nj) THEN |
---|
618 | route_flow_rp1(ij_rp1)=-1 ! if route outside my local domain, no routing (will be done by other process) |
---|
619 | ELSE |
---|
620 | route_flow_rp1(ij_rp1)=ir+ni*(jr-1) ! define the cell where to flow |
---|
621 | ENDIF |
---|
622 | ENDIF |
---|
623 | |
---|
624 | ENDIF |
---|
625 | ENDDO |
---|
626 | ENDDO |
---|
627 | |
---|
628 | routing_mask_r(:)=reshape(routing_mask_rp1(1:ni,1:nj),(/ni*nj/)) |
---|
629 | |
---|
630 | is_streamflow_r(:)=.NOT. (is_lakeinflow_r(:) .OR. is_coastalflow_r(:) .OR. is_riverflow_r(:)) .AND. routing_mask_r(:) |
---|
631 | |
---|
632 | DO ij=1,nbp_mpi |
---|
633 | routing_weight(ij)=contfrac_mpi(ij)*area_mpi(ij) |
---|
634 | ENDDO |
---|
635 | |
---|
636 | routing_weight_r(:)=0 |
---|
637 | DO ij=1,nbpt_r |
---|
638 | IF (frac_routing_r(ij)>0) routing_weight_r(ij)=1./frac_routing_r(ij) |
---|
639 | ENDDO |
---|
640 | |
---|
641 | routing_weight_coast_r(:)=0. |
---|
642 | DO ij=1,nbpt_r |
---|
643 | IF (frac_routing_coast_r(ij)>0) routing_weight_coast_r(ij)=1./frac_routing_coast_r(ij) |
---|
644 | ENDDO |
---|
645 | |
---|
646 | CALL xios_send_field("routing_weight_coast_r",routing_weight_coast_r) |
---|
647 | ! looking for basins |
---|
648 | basins_count_mpi=0 |
---|
649 | DO ij=1,nbpt_r |
---|
650 | IF (basins_extended_r(ij) > basins_count_mpi) basins_count_mpi=basins_extended_r(ij) |
---|
651 | ENDDO |
---|
652 | CALL MPI_ALLREDUCE(basins_count_mpi,basins_count,1,MPI_INT_ORCH, MPI_MAX,MPI_COMM_ORCH,ierr) |
---|
653 | |
---|
654 | ELSE |
---|
655 | |
---|
656 | nbpt_r=0 |
---|
657 | |
---|
658 | ENDIF |
---|
659 | |
---|
660 | END SUBROUTINE routing_simple_init_2 |
---|
661 | |
---|
662 | |
---|
663 | !! ================================================================================================================================ |
---|
664 | !! SUBROUTINE : routing_simple_flow |
---|
665 | !! |
---|
666 | !>\BRIEF This subroutine computes the transport of water in the various reservoirs |
---|
667 | !! (including ponds and floodplains) and the water withdrawals from the reservoirs for irrigation. |
---|
668 | !! |
---|
669 | !! DESCRIPTION (definitions, functional, design, flags) : |
---|
670 | !! This will first compute the amount of water which flows out of each of the 3 reservoirs using the assumption of an |
---|
671 | !! exponential decrease of water in the reservoir (see Hagemann S and Dumenil L. (1998)). Then we compute the fluxes |
---|
672 | !! for floodplains and ponds. All this will then be used in order to update each of the basins : taking water out of |
---|
673 | !! the up-stream basin and adding it to the down-stream one. |
---|
674 | !! As this step happens globaly we have to stop the parallel processing in order to exchange the information. Once |
---|
675 | !! all reservoirs are updated we deal with irrigation. The final step is to compute diagnostic fluxes. Among them |
---|
676 | !! the hydrographs of the largest rivers we have chosen to monitor. |
---|
677 | !! |
---|
678 | !! RECENT CHANGE(S): None |
---|
679 | !! |
---|
680 | !! MAIN OUTPUT VARIABLE(S): lakeinflow, returnflow, reinfiltration, irrigation, riverflow, coastalflow, hydrographs, flood_frac, flood_res |
---|
681 | !! |
---|
682 | !! REFERENCES : |
---|
683 | !! - Ngo-Duc, T., K. Laval, G. Ramillien, J. Polcher, and A. Cazenave (2007) |
---|
684 | !! Validation of the land water storage simulated by Organising Carbon and Hydrology in Dynamic Ecosystems (ORCHIDEE) with Gravity Recovery and Climate Experiment (GRACE) data. |
---|
685 | !! Water Resour. Res., 43, W04427, doi:10.1029/2006WR004941. |
---|
686 | !! * Irrigation: |
---|
687 | !! - de Rosnay, P., J. Polcher, K. Laval, and M. Sabre (2003) |
---|
688 | !! Integrated parameterization of irrigation in the land surface model ORCHIDEE. Validation over Indian Peninsula. |
---|
689 | !! Geophys. Res. Lett., 30(19), 1986, doi:10.1029/2003GL018024. |
---|
690 | !! - A.C. Vivant (2003) |
---|
691 | !! Les plaines d'inondations et l'irrigation dans ORCHIDEE, impacts de leur prise en compte. |
---|
692 | !! , , 51pp. |
---|
693 | !! - N. Culson (2004) |
---|
694 | !! Impact de l'irrigation sur le cycle de l'eau |
---|
695 | !! Master thesis, Paris VI University, 55pp. |
---|
696 | !! - X.-T. Nguyen-Vinh (2005) |
---|
697 | !! Analyse de l'impact de l'irrigation en Amerique du Nord - plaine du Mississippi - sur la climatologie regionale |
---|
698 | !! Master thesis, Paris VI University, 33pp. |
---|
699 | !! - M. Guimberteau (2006) |
---|
700 | !! Analyse et modifications proposees de la modelisation de l'irrigation dans un modele de surface. |
---|
701 | !! Master thesis, Paris VI University, 46pp. |
---|
702 | !! - Guimberteau M. (2010) |
---|
703 | !! Modelisation de l'hydrologie continentale et influences de l'irrigation sur le cycle de l'eau. |
---|
704 | !! Ph.D. thesis, Paris VI University, 195pp. |
---|
705 | !! - Guimberteau M., Laval K., Perrier A. and Polcher J. (2011). |
---|
706 | !! Global effect of irrigation and its impact on the onset of the Indian summer monsoon. |
---|
707 | !! In press, Climate Dynamics, doi: 10.1007/s00382-011-1252-5. |
---|
708 | !! * Floodplains: |
---|
709 | !! - A.C. Vivant (2002) |
---|
710 | !! L'ecoulement lateral de l'eau sur les surfaces continentales. Prise en compte des plaines d'inondations dans ORCHIDEE. |
---|
711 | !! Master thesis, Paris VI University, 46pp. |
---|
712 | !! - A.C. Vivant (2003) |
---|
713 | !! Les plaines d'inondations et l'irrigation dans ORCHIDEE, impacts de leur prise en compte. |
---|
714 | !! , , 51pp. |
---|
715 | !! - T. d'Orgeval (2006) |
---|
716 | !! Impact du changement climatique sur le cycle de l'eau en Afrique de l'Ouest: modelisation et incertitudes. |
---|
717 | !! Ph.D. thesis, Paris VI University, 188pp. |
---|
718 | !! - T. d'Orgeval, J. Polcher, and P. de Rosnay (2008) |
---|
719 | !! Sensitivity of the West African hydrological cycle in ORCHIDEE to infiltration processes. |
---|
720 | !! Hydrol. Earth Syst. Sci., 12, 1387-1401 |
---|
721 | !! - M. Guimberteau, G. Drapeau, J. Ronchail, B. Sultan, J. Polcher, J.-M. Martinez, C. Prigent, J.-L. Guyot, G. Cochonneau, |
---|
722 | !! J. C. Espinoza, N. Filizola, P. Fraizy, W. Lavado, E. De Oliveira, R. Pombosa, L. Noriega, and P. Vauchel (2011) |
---|
723 | !! Discharge simulation in the sub-basins of the Amazon using ORCHIDEE forced by new datasets. |
---|
724 | !! Hydrol. Earth Syst. Sci. Discuss., 8, 11171-11232, doi:10.5194/hessd-8-11171-2011 |
---|
725 | !! |
---|
726 | !! FLOWCHART :None |
---|
727 | !! \n |
---|
728 | !_ ================================================================================================================================ |
---|
729 | |
---|
730 | SUBROUTINE routing_simple_flow(nbpt, dt_routing, lalo, floodout, runoff_omp, drainage_omp, & |
---|
731 | vegtot, totnobio, transpot_mean, precip, humrel, k_litt, floodtemp, reinf_slope, & |
---|
732 | lakeinflow_omp, returnflow, reinfiltration, irrigation, riverflow_omp, & |
---|
733 | coastalflow_omp, hydrographs, slowflow_diag, flood_frac, flood_res) |
---|
734 | |
---|
735 | USE xios |
---|
736 | USE grid, ONLY : area |
---|
737 | IMPLICIT NONE |
---|
738 | INCLUDE "mpif.h" |
---|
739 | |
---|
740 | !! 0 Variable and parameter description |
---|
741 | !! 0.1 Input variables |
---|
742 | |
---|
743 | INTEGER(i_std), INTENT(in) :: nbpt !! Domain size (unitless) |
---|
744 | REAL(r_std), INTENT (in) :: dt_routing !! Routing time step (s) |
---|
745 | REAL(r_std), INTENT(in) :: lalo(nbpt,2) !! Vector of latitude and longitudes |
---|
746 | REAL(r_std), INTENT(in) :: runoff_omp(nbpt) !! Grid-point runoff (kg/m^2/dt) |
---|
747 | REAL(r_std), INTENT(in) :: floodout(nbpt) !! Grid-point flow out of floodplains (kg/m^2/dt) |
---|
748 | REAL(r_std), INTENT(in) :: drainage_omp(nbpt) !! Grid-point drainage (kg/m^2/dt) |
---|
749 | REAL(r_std), INTENT(in) :: vegtot(nbpt) !! Potentially vegetated fraction (unitless;0-1) |
---|
750 | REAL(r_std), INTENT(in) :: totnobio(nbpt) !! Other areas which can not have vegetation |
---|
751 | REAL(r_std), INTENT(in) :: transpot_mean(nbpt) !! Mean potential transpiration of the vegetation (kg/m^2/dt) |
---|
752 | REAL(r_std), INTENT(in) :: precip(nbpt) !! Rainfall (kg/m^2/dt) |
---|
753 | REAL(r_std), INTENT(in) :: humrel(nbpt) !! Soil moisture stress, root extraction potential (unitless) |
---|
754 | REAL(r_std), INTENT(in) :: k_litt(nbpt) !! Averaged conductivity for saturated infiltration in the 'litter' layer (kg/m^2/dt) |
---|
755 | REAL(r_std), INTENT(in) :: floodtemp(nbpt) !! Temperature to decide if floodplains work (K) |
---|
756 | REAL(r_std), INTENT(in) :: reinf_slope(nbpt) !! Coefficient which determines the reinfiltration ratio in the grid box due to flat areas (unitless;0-1) |
---|
757 | |
---|
758 | !! 0.2 Output variables |
---|
759 | |
---|
760 | REAL(r_std), INTENT(out) :: lakeinflow_omp(nbpt) !! Water inflow to the lakes (kg/dt) |
---|
761 | REAL(r_std), INTENT(out) :: returnflow(nbpt) !! The water flow from lakes and swamps which returns into the grid box. |
---|
762 | !! This water will go back into the hydrol or hydrolc module to allow re-evaporation (kg/m^2/dt) |
---|
763 | REAL(r_std), INTENT(out) :: reinfiltration(nbpt) !! Water flow from ponds and floodplains which returns to the grid box (kg/m^2/dt) |
---|
764 | REAL(r_std), INTENT(out) :: irrigation(nbpt) !! Irrigation flux. This is the water taken from the reservoirs and beeing put into the upper layers of the soil (kg/m^2/dt) |
---|
765 | REAL(r_std), INTENT(out) :: riverflow_omp(nbpt) !! Outflow of the major rivers. The flux will be located on the continental grid but this should be a coastal point (kg/dt) |
---|
766 | REAL(r_std), INTENT(out) :: coastalflow_omp(nbpt) !! Outflow on coastal points by small basins. This is the water which flows in a disperse way into the ocean (kg/dt) |
---|
767 | REAL(r_std), INTENT(out) :: hydrographs(nbpt) !! Hydrographs at the outflow of the grid box for major basins (kg/dt) |
---|
768 | REAL(r_std), INTENT(out) :: slowflow_diag(nbpt) !! Hydrographs of slow_flow = routed slow_flow for major basins (kg/dt) |
---|
769 | REAL(r_std), INTENT(out) :: flood_frac(nbpt) !! Flooded fraction of the grid box (unitless;0-1) |
---|
770 | REAL(r_std), INTENT(out) :: flood_res(nbpt) !! Diagnostic of water amount in the floodplains reservoir (kg) |
---|
771 | |
---|
772 | !! 0.4 Local variables |
---|
773 | REAL(r_std) :: runoff(nbp_mpi) !! Grid-point runoff (kg/dt) |
---|
774 | REAL(r_std) :: drainage(nbp_mpi) !! Grid-point drainage (kg/dt) |
---|
775 | REAL(r_std) :: riverflow(nbp_mpi) |
---|
776 | REAL(r_std) :: coastalflow(nbp_mpi) |
---|
777 | REAL(r_std) :: lakeinflow(nbp_mpi) |
---|
778 | REAL(r_std) :: fast_diag_mpi(nbp_mpi) |
---|
779 | REAL(r_std) :: slow_diag_mpi(nbp_mpi) |
---|
780 | REAL(r_std) :: stream_diag_mpi(nbp_mpi) |
---|
781 | |
---|
782 | ! from input model -> routing_grid |
---|
783 | REAL(r_std) :: runoff_r(nbpt_r) !! Grid-point runoff (kg/m^2/dt) |
---|
784 | REAL(r_std) :: floodout_r(nbpt_r) !! Grid-point flow out of floodplains (kg/m^2/dt) |
---|
785 | REAL(r_std) :: drainage_r(nbpt_r) !! Grid-point drainage (kg/m^2/dt) |
---|
786 | REAL(r_std) :: vegtot_r(nbpt_r) !! Potentially vegetated fraction (unitless;0-1) |
---|
787 | REAL(r_std) :: totnobio_r(nbpt_r) !! Other areas which can not have vegetation |
---|
788 | REAL(r_std) :: transpot_mean_r(nbpt_r) !! Mean potential transpiration of the vegetation (kg/m^2/dt) |
---|
789 | REAL(r_std) :: precip_r(nbpt_r) !! Rainfall (kg/m^2/dt) |
---|
790 | REAL(r_std) :: humrel_r(nbpt_r) !! Soil moisture stress, root extraction potential (unitless) |
---|
791 | REAL(r_std) :: k_litt_r(nbpt_r) !! Averaged conductivity for saturated infiltration in the 'litter' layer (kg/m^2/dt) |
---|
792 | REAL(r_std) :: floodtemp_r(nbpt_r) !! Temperature to decide if floodplains work (K) |
---|
793 | REAL(r_std) :: reinf_slope_r(nbpt_r) !! Coefficient which determines the reinfiltration ratio in the grid box due to flat areas (unitless;0-1) |
---|
794 | |
---|
795 | |
---|
796 | |
---|
797 | REAL(r_std), DIMENSION(nbpt_r) :: fast_flow_r !! Outflow from the fast reservoir (kg/dt) |
---|
798 | REAL(r_std), DIMENSION(nbpt_r) :: slow_flow_r !! Outflow from the slow reservoir (kg/dt) |
---|
799 | REAL(r_std), DIMENSION(nbpt_r) :: stream_flow_r !! Outflow from the stream reservoir (kg/dt) |
---|
800 | REAL(r_std), DIMENSION(nbpt_r) :: hydrographs_r !! hydrograph (kg/dt) |
---|
801 | REAL(r_std), DIMENSION(nbpt_r) :: flood_flow_r !! Outflow from the floodplain reservoir (kg/dt) |
---|
802 | REAL(r_std), DIMENSION(nbpt_r) :: pond_inflow_r !! Inflow to the pond reservoir (kg/dt) |
---|
803 | REAL(r_std), DIMENSION(nbpt_r) :: pond_drainage_r !! Drainage from pond (kg/m^2/dt) |
---|
804 | REAL(r_std), DIMENSION(nbpt_r) :: flood_drainage_r !! Drainage from floodplains (kg/m^2/dt) |
---|
805 | REAL(r_std), DIMENSION(nbpt_r) :: return_swamp_r !! Inflow to the swamp (kg/dt) |
---|
806 | ! |
---|
807 | ! Irrigation per basin |
---|
808 | ! |
---|
809 | REAL(r_std), DIMENSION(nbpt_r) :: irrig_needs_r !! Total irrigation requirement (water requirements by the crop for its optimal growth) (kg) |
---|
810 | REAL(r_std), DIMENSION(nbpt_r) :: irrig_actual_r !! Possible irrigation according to the water availability in the reservoirs (kg) |
---|
811 | REAL(r_std), DIMENSION(nbpt_r) :: irrig_deficit_r !! Amount of water missing for irrigation (kg) |
---|
812 | REAL(r_std), DIMENSION(nbpt_r) :: irrig_adduct_r !! Amount of water carried over from other basins for irrigation (kg) |
---|
813 | ! |
---|
814 | REAL(r_std), DIMENSION(nbpt_r) :: transport_r !! Water transport between basins (kg/dt) |
---|
815 | REAL(r_std), DIMENSION(nbpt_r) :: floods_r !! Water flow in to the floodplains (kg/dt) |
---|
816 | REAL(r_std), DIMENSION(nbpt_r) :: potflood_r !! Potential inflow to the swamps (kg/dt) |
---|
817 | REAL(r_std), DIMENSION(nbpt) :: tobeflooded !! Maximal surface which can be inundated in each grid box (m^2) |
---|
818 | REAL(r_std), DIMENSION(nbpt) :: totarea !! Total area of basin (m^2) |
---|
819 | REAL(r_std), DIMENSION(nbpt) :: totflood !! Total amount of water in the floodplains reservoir (kg) (not used) |
---|
820 | REAL(r_std), DIMENSION(nbasmax) :: pond_excessflow !! |
---|
821 | REAL(r_std) :: flow !! Outflow computation for the reservoirs (kg/dt) |
---|
822 | REAL(r_std) :: floodindex !! Fraction of grid box area inundated (unitless;0-1) |
---|
823 | REAL(r_std) :: pondex !! |
---|
824 | REAL(r_std) :: flood_frac_pot !! Total fraction of the grid box which is flooded at optimum repartition (unitless;0-1) |
---|
825 | REAL(r_std) :: stream_tot !! Total water amount in the stream reservoirs (kg) |
---|
826 | REAL(r_std) :: adduction !! Importation of water from a stream reservoir of a neighboring grid box (kg) |
---|
827 | REAL(r_std), DIMENSION(8,nbasmax) :: streams_around !! Stream reservoirs of the neighboring grid boxes (kg) |
---|
828 | INTEGER(i_std), DIMENSION(8) :: igrd !! |
---|
829 | INTEGER(i_std), DIMENSION(2) :: ff !! |
---|
830 | INTEGER(i_std), DIMENSION(1) :: fi !! |
---|
831 | INTEGER(i_std) :: ig, ib, ib2, ig2 !! Indices (unitless) |
---|
832 | INTEGER(i_std) :: rtg, rtb, in !! Indices (unitless) |
---|
833 | INTEGER(i_std) :: ier !! Error handling |
---|
834 | INTEGER(i_std) :: isplit |
---|
835 | |
---|
836 | 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) |
---|
837 | REAL(r_std), DIMENSION(nbpt_r) :: flood_frac_bas_r |
---|
838 | REAL(r_std), DIMENSION(nbpt_r) :: flood_reservoir_r |
---|
839 | REAL(r_std), DIMENSION(nbpt_r) :: pond_reservoir_r |
---|
840 | |
---|
841 | REAL(r_std), DIMENSION(nbpt_r) :: lakeinflow_r |
---|
842 | REAL(r_std), DIMENSION(nbpt_r) :: coastalflow_r |
---|
843 | REAL(r_std), DIMENSION(nbpt_r) :: riverflow_r |
---|
844 | |
---|
845 | REAL(r_std), DIMENSION(nbpt_r) :: flow_r |
---|
846 | REAL(r_std), DIMENSION(nbpt_rp1) :: flow_rp1 |
---|
847 | REAL(r_std) :: water_balance_before, water_balance_after |
---|
848 | REAL(r_std) :: basins_riverflow_mpi(0:basins_count) |
---|
849 | REAL(r_std) :: basins_riverflow(0:basins_count) |
---|
850 | |
---|
851 | INTEGER :: ierr |
---|
852 | |
---|
853 | !_ ================================================================================================================================ |
---|
854 | ! |
---|
855 | |
---|
856 | |
---|
857 | irrig_netereq(:) = zero |
---|
858 | totarea(:) = zero |
---|
859 | totflood(:) = zero |
---|
860 | |
---|
861 | flood_frac_bas_r(:)=zero |
---|
862 | ! |
---|
863 | !> The outflow fluxes from the three reservoirs are computed. |
---|
864 | !> The outflow of volume of water Vi into the reservoir i is assumed to be linearly related to its volume. |
---|
865 | !> The water travel simulated by the routing scheme is dependent on the water retention index topo_resid |
---|
866 | !> given by a 0.5 degree resolution map for each pixel performed from a simplification of Manning's formula |
---|
867 | !> (Dingman, 1994; Ducharne et al., 2003). |
---|
868 | !> The resulting product of tcst (in day/m) and topo_resid (in m) represents the time constant (day) |
---|
869 | !> which is an e-folding time, the time necessary for the water amount |
---|
870 | !> in the stream reservoir to decrease by a factor e. Hence, it gives an order of |
---|
871 | !> magnitude of the travel time through this reservoir between |
---|
872 | !> the sub-basin considered and its downstream neighbor. |
---|
873 | CALL gather_omp(runoff_omp,runoff) |
---|
874 | CALL gather_omp(drainage_omp, drainage) |
---|
875 | |
---|
876 | IF (is_omp_root) THEN |
---|
877 | |
---|
878 | hydrographs_r(:)=0 |
---|
879 | |
---|
880 | DO isplit=1,split_routing |
---|
881 | |
---|
882 | |
---|
883 | irrig_needs_r(:) = zero |
---|
884 | irrig_actual_r(:) = zero |
---|
885 | irrig_deficit_r(:) = zero |
---|
886 | irrig_adduct_r(:) = zero |
---|
887 | |
---|
888 | DO ig=1,nbpt_r |
---|
889 | IF ( routing_mask_r(ig) ) THEN |
---|
890 | ! |
---|
891 | ! Each of the fluxes is limited by the water in the reservoir and a small margin |
---|
892 | ! (min_reservoir) to avoid rounding errors. |
---|
893 | ! |
---|
894 | ! IF (fast_reservoir_r(ig)>1e-30) THEN |
---|
895 | ! PRINT*,"fast_reservoir > 0 ",ig, fast_reservoir_r(ig), topoind_r(ig) |
---|
896 | ! ENDIF |
---|
897 | flow = MIN(fast_reservoir_r(ig)/((topoind_r(ig)/1000.)*fast_tcst*one_day/(dt_routing/split_routing)),& |
---|
898 | & fast_reservoir_r(ig)-min_sechiba) |
---|
899 | fast_flow_r(ig) = MAX(flow, zero) |
---|
900 | ! |
---|
901 | flow = MIN(slow_reservoir_r(ig)/((topoind_r(ig)/1000.)*slow_tcst*one_day/(dt_routing/split_routing)),& |
---|
902 | & slow_reservoir_r(ig)-min_sechiba) |
---|
903 | slow_flow_r(ig) = MAX(flow, zero) |
---|
904 | ! |
---|
905 | flow = MIN(stream_reservoir_r(ig)/((topoind_r(ig)/1000.)*stream_tcst* & |
---|
906 | & MAX(un-SQRT(flood_frac_bas_r(ig)),min_sechiba)*one_day/(dt_routing/split_routing)),& |
---|
907 | & stream_reservoir_r(ig)-min_sechiba) |
---|
908 | stream_flow_r(ig) = MAX(flow, zero) |
---|
909 | ! |
---|
910 | ! |
---|
911 | ELSE |
---|
912 | fast_flow_r(ig) = zero |
---|
913 | slow_flow_r(ig) = zero |
---|
914 | stream_flow_r(ig) = zero |
---|
915 | ENDIF |
---|
916 | ENDDO |
---|
917 | |
---|
918 | DO ig=1,nbpt_r |
---|
919 | flood_drainage_r(ig) = zero |
---|
920 | flood_flow_r(ig) = zero |
---|
921 | flood_reservoir_r(ig) = zero |
---|
922 | ENDDO |
---|
923 | |
---|
924 | DO ig=1,nbpt_r |
---|
925 | pond_inflow_r(ig) = zero |
---|
926 | pond_drainage_r(ig) = zero |
---|
927 | pond_reservoir_r(ig) = zero |
---|
928 | ENDDO |
---|
929 | |
---|
930 | |
---|
931 | !- |
---|
932 | !- Compute the transport |
---|
933 | !- |
---|
934 | |
---|
935 | flow_r(:)=fast_flow_r(:) + slow_flow_r(:) + stream_flow_r(:) |
---|
936 | |
---|
937 | CALL xios_send_field("flow_r",flow_r) ! transfer halo |
---|
938 | CALL xios_recv_field("flow_rp1",flow_rp1) |
---|
939 | |
---|
940 | transport_r(:)=0 |
---|
941 | |
---|
942 | DO ig=1,nbpt_rp1 |
---|
943 | IF ( route_flow_rp1(ig) > 0 ) THEN |
---|
944 | transport_r(route_flow_rp1(ig))=transport_r(route_flow_rp1(ig))+ flow_rp1(ig) |
---|
945 | ENDIF |
---|
946 | ENDDO |
---|
947 | |
---|
948 | |
---|
949 | !- |
---|
950 | !- Do the floodings - First initialize |
---|
951 | !- |
---|
952 | return_swamp_r(:)=zero |
---|
953 | floods_r(:)=zero |
---|
954 | |
---|
955 | ! |
---|
956 | ! Update all reservoirs |
---|
957 | !> The slow and deep reservoir (slow_reservoir) collect the deep drainage whereas the |
---|
958 | !> fast_reservoir collects the computed surface runoff. Both discharge into a third reservoir |
---|
959 | !> (stream_reservoir) of the next sub-basin downstream. |
---|
960 | !> Water from the floodplains reservoir (flood_reservoir) flows also into the stream_reservoir of the next sub-basin downstream. |
---|
961 | !> Water that flows into the pond_reservoir is withdrawn from the fast_reservoir. |
---|
962 | ! |
---|
963 | runoff=runoff*routing_weight |
---|
964 | CALL xios_send_field("routing_runoff",runoff) ! interp conservative model -> routing |
---|
965 | CALL xios_recv_field("routing_runoff_r",runoff_r) |
---|
966 | drainage=drainage*routing_weight |
---|
967 | CALL xios_send_field("routing_drainage",drainage) ! interp conservative model -> routing |
---|
968 | CALL xios_recv_field("routing_drainage_r",drainage_r) |
---|
969 | |
---|
970 | CALL xios_send_field("fast_reservoir_r",fast_reservoir_r) |
---|
971 | CALL xios_send_field("slow_reservoir_r",slow_reservoir_r) |
---|
972 | CALL xios_send_field("stream_reservoir_r",stream_reservoir_r) |
---|
973 | |
---|
974 | WHERE (.NOT. routing_mask_r) runoff_r(:)=0 |
---|
975 | WHERE (.NOT. routing_mask_r) drainage_r(:)=0 |
---|
976 | |
---|
977 | CALL MPI_ALLREDUCE(sum(runoff+drainage),water_balance_before,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr) |
---|
978 | CALL MPI_ALLREDUCE(sum(runoff_r+drainage_r),water_balance_after,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr) |
---|
979 | |
---|
980 | PRINT *,"routing water Balance ; before : ", water_balance_before," ; after : ",water_balance_after, & |
---|
981 | " ; delta : ", 100.*(water_balance_after-water_balance_before)/(0.5*(water_balance_after+water_balance_before)),"%" |
---|
982 | |
---|
983 | CALL xios_recv_field("water_balance_before",water_balance_before) |
---|
984 | |
---|
985 | |
---|
986 | DO ig=1,nbpt_r |
---|
987 | IF ( routing_mask_r(ig) ) THEN |
---|
988 | |
---|
989 | ! |
---|
990 | fast_reservoir_r(ig) = fast_reservoir_r(ig) + runoff_r(ig) - & |
---|
991 | & fast_flow_r(ig) - pond_inflow_r(ig) |
---|
992 | ! |
---|
993 | slow_reservoir_r(ig) = slow_reservoir_r(ig) + drainage_r(ig) - & |
---|
994 | & slow_flow_r(ig) |
---|
995 | ! |
---|
996 | |
---|
997 | stream_reservoir_r(ig) = stream_reservoir_r(ig) + flood_flow_r(ig) - & |
---|
998 | & stream_flow_r(ig) - return_swamp_r(ig) - floods_r(ig) |
---|
999 | |
---|
1000 | IF (is_streamflow_r(ig)) stream_reservoir_r(ig)= stream_reservoir_r(ig) + transport_r(ig) |
---|
1001 | ! |
---|
1002 | flood_reservoir_r(ig) = flood_reservoir_r(ig) + floods_r(ig) - & |
---|
1003 | & flood_flow_r(ig) |
---|
1004 | ! |
---|
1005 | pond_reservoir_r(ig) = pond_reservoir_r(ig) + pond_inflow_r(ig) - pond_drainage_r(ig) |
---|
1006 | ! |
---|
1007 | IF ( flood_reservoir_r(ig) .LT. zero ) THEN |
---|
1008 | IF ( check_reservoir ) THEN |
---|
1009 | WRITE(numout,*) "WARNING : negative flood reservoir at :", ig, ". Problem is being corrected." |
---|
1010 | WRITE(numout,*) "flood_reservoir, floods, flood_flow : ", flood_reservoir_r(ig), floods_r(ig), & |
---|
1011 | & flood_flow_r(ig) |
---|
1012 | ENDIF |
---|
1013 | stream_reservoir_r(ig) = stream_reservoir_r(ib) + flood_reservoir_r(ib) |
---|
1014 | flood_reservoir_r(ib) = zero |
---|
1015 | ENDIF |
---|
1016 | ! |
---|
1017 | IF ( stream_reservoir_r(ig) .LT. zero ) THEN |
---|
1018 | IF ( check_reservoir ) THEN |
---|
1019 | WRITE(numout,*) "WARNING : negative stream reservoir at :", ig, ". Problem is being corrected." |
---|
1020 | WRITE(numout,*) "stream_reservoir, flood_flow, transport : ", stream_reservoir_r(ig), flood_flow_r(ig), & |
---|
1021 | & transport_r(ig) |
---|
1022 | WRITE(numout,*) "stream_flow, return_swamp, floods :", stream_flow_r(ig), return_swamp_r(ig), floods_r(ig) |
---|
1023 | ENDIF |
---|
1024 | fast_reservoir_r(ig) = fast_reservoir_r(ig) + stream_reservoir_r(ig) |
---|
1025 | stream_reservoir_r(ig) = zero |
---|
1026 | ENDIF |
---|
1027 | ! |
---|
1028 | IF ( fast_reservoir_r(ig) .LT. zero ) THEN |
---|
1029 | IF ( check_reservoir ) THEN |
---|
1030 | WRITE(numout,*) "WARNING : negative fast reservoir at :", ig, ". Problem is being corrected." |
---|
1031 | WRITE(numout,*) "fast_reservoir, runoff, fast_flow, ponf_inflow : ", fast_reservoir_r(ig), & |
---|
1032 | &runoff_r(ig), fast_flow_r(ig), pond_inflow_r(ig) |
---|
1033 | ENDIF |
---|
1034 | slow_reservoir_r(ig) = slow_reservoir_r(ig) + fast_reservoir_r(ig) |
---|
1035 | fast_reservoir_r(ig) = zero |
---|
1036 | ENDIF |
---|
1037 | |
---|
1038 | IF ( slow_reservoir_r(ig) .LT. - min_sechiba ) THEN |
---|
1039 | WRITE(numout,*) 'WARNING : There is a negative reservoir at :', ig |
---|
1040 | WRITE(numout,*) 'WARNING : slowr, slow_flow, drainage', & |
---|
1041 | & slow_reservoir_r(ig), slow_flow_r(ig), drainage_r(ig) |
---|
1042 | WRITE(numout,*) 'WARNING : pondr, pond_inflow, pond_drainage', & |
---|
1043 | & pond_reservoir_r(ig), pond_inflow_r(ig), pond_drainage_r(ig) |
---|
1044 | CALL ipslerr_p(2, 'routing_simple_flow', 'WARNING negative slow_reservoir.','','') |
---|
1045 | ENDIF |
---|
1046 | |
---|
1047 | ENDIF |
---|
1048 | ENDDO |
---|
1049 | |
---|
1050 | DO ig=1,nbpt_r |
---|
1051 | IF ( routing_mask_r(ig) ) THEN |
---|
1052 | hydrographs_r(ig)=hydrographs_r(ig)+fast_flow_r(ig)+slow_flow_r(ig)+stream_flow_r(ig) |
---|
1053 | ENDIF |
---|
1054 | ENDDO |
---|
1055 | |
---|
1056 | ENDDO ! isplit |
---|
1057 | |
---|
1058 | ! |
---|
1059 | ! |
---|
1060 | ! |
---|
1061 | ! The three output types for the routing : endoheric basins,, rivers and |
---|
1062 | ! diffuse coastal flow. |
---|
1063 | ! |
---|
1064 | lakeinflow_r(:)=0 |
---|
1065 | coastalflow_r(:)=0 |
---|
1066 | riverflow_r(:)=0 |
---|
1067 | basins_riverflow_mpi(:)=0 |
---|
1068 | |
---|
1069 | DO ig=1,nbpt_r |
---|
1070 | IF ( routing_mask_r(ig) ) THEN |
---|
1071 | |
---|
1072 | IF (is_lakeinflow_r(ig)) THEN |
---|
1073 | lakeinflow_r(ig) = transport_r(ig) |
---|
1074 | basins_riverflow_mpi(basins_extended_r(ig)) = & |
---|
1075 | basins_riverflow_mpi(basins_extended_r(ig))+lakeinflow_r(ig) |
---|
1076 | ENDIF |
---|
1077 | |
---|
1078 | IF (is_coastalflow_r(ig)) THEN |
---|
1079 | coastalflow_r(ig) = transport_r(ig) |
---|
1080 | basins_riverflow_mpi(basins_extended_r(ig)) = & |
---|
1081 | basins_riverflow_mpi(basins_extended_r(ig))+coastalflow_r(ig) |
---|
1082 | ENDIF |
---|
1083 | |
---|
1084 | IF (is_riverflow_r(ig)) THEN |
---|
1085 | riverflow_r(ig) = transport_r(ig) |
---|
1086 | basins_riverflow_mpi(basins_extended_r(ig)) = & |
---|
1087 | basins_riverflow_mpi(basins_extended_r(ig))+riverflow_r(ig) |
---|
1088 | ENDIF |
---|
1089 | |
---|
1090 | ENDIF |
---|
1091 | ENDDO |
---|
1092 | |
---|
1093 | |
---|
1094 | CALL MPI_ALLREDUCE(basins_riverflow_mpi,basins_riverflow,basins_count+1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr) |
---|
1095 | CALL xios_send_field("basins_riverflow",basins_riverflow(1:basins_out)/1000./dt_routing) |
---|
1096 | |
---|
1097 | ! |
---|
1098 | flood_res = flood_diag + pond_diag |
---|
1099 | ! |
---|
1100 | CALL xios_send_field("routing_lakeinflow_r" ,lakeinflow_r*routing_weight_r) |
---|
1101 | CALL xios_recv_field("routing_lakeinflow" ,lakeinflow) |
---|
1102 | CALL xios_send_field("routing_coastalflow_r" ,coastalflow_r*routing_weight_coast_r) |
---|
1103 | CALL xios_recv_field("routing_coastalflow_temp" ,coastalflow) |
---|
1104 | WHERE(.NOT. coast_mask(:)) coastalflow(:)=0 |
---|
1105 | CALL xios_send_field("routing_coastalflow" ,coastalflow) |
---|
1106 | |
---|
1107 | CALL xios_send_field("routing_riverflow_r" ,riverflow_r*routing_weight_coast_r) |
---|
1108 | CALL xios_recv_field("routing_riverflow_temp" ,riverflow) |
---|
1109 | WHERE(.NOT. coast_mask(:)) riverflow(:)=0 |
---|
1110 | CALL xios_send_field("routing_riverflow" ,riverflow) |
---|
1111 | |
---|
1112 | CALL xios_send_field("out_flow",lakeinflow+coastalflow+riverflow) |
---|
1113 | CALL xios_send_field("routing_hydrographs_r", hydrographs_r/1000./dt_routing) |
---|
1114 | |
---|
1115 | ! diag |
---|
1116 | CALL xios_send_field("routing_fast_reservoir_r" , fast_reservoir_r) |
---|
1117 | CALL xios_recv_field("routing_fast_reservoir" , fast_diag_mpi) |
---|
1118 | |
---|
1119 | CALL xios_send_field("routing_slow_reservoir_r" , slow_reservoir_r) |
---|
1120 | CALL xios_recv_field("routing_slow_reservoir" , slow_diag_mpi) |
---|
1121 | |
---|
1122 | CALL xios_send_field("routing_stream_reservoir_r" , stream_reservoir_r) |
---|
1123 | CALL xios_recv_field("routing_stream_reservoir" , stream_diag_mpi) |
---|
1124 | |
---|
1125 | CALL xios_recv_field("water_balance_after" ,water_balance_after) |
---|
1126 | |
---|
1127 | PRINT *,"routing water Balance ; before : ", water_balance_before," ; after : ",water_balance_after, & |
---|
1128 | " ; delta : ", 100*(water_balance_after-water_balance_before)/(0.5*(water_balance_after+water_balance_before)),"%" |
---|
1129 | |
---|
1130 | |
---|
1131 | |
---|
1132 | ENDIF ! is_omp_root |
---|
1133 | |
---|
1134 | CALL scatter_omp(riverflow,riverflow_omp) |
---|
1135 | CALL scatter_omp(coastalflow,coastalflow_omp) |
---|
1136 | CALL scatter_omp(lakeinflow,lakeinflow_omp) |
---|
1137 | CALL scatter_omp(fast_diag_mpi,fast_diag) |
---|
1138 | CALL scatter_omp(slow_diag_mpi,slow_diag) |
---|
1139 | CALL scatter_omp(stream_diag_mpi,stream_diag) |
---|
1140 | |
---|
1141 | DO ig=1,nbpt |
---|
1142 | fast_diag(ig)=fast_diag(ig)/area(ig) |
---|
1143 | slow_diag(ig)=slow_diag(ig)/area(ig) |
---|
1144 | stream_diag(ig)=stream_diag(ig)/area(ig) |
---|
1145 | ENDDO |
---|
1146 | |
---|
1147 | |
---|
1148 | flood_frac(:) = zero |
---|
1149 | flood_height(:) = zero |
---|
1150 | flood_frac_bas_r(:) = zero |
---|
1151 | |
---|
1152 | returnflow(:) = zero |
---|
1153 | reinfiltration(:) = zero |
---|
1154 | |
---|
1155 | ! |
---|
1156 | ! |
---|
1157 | ! Compute the fluxes which leave the routing scheme |
---|
1158 | ! |
---|
1159 | ! Lakeinflow is in Kg/dt |
---|
1160 | ! returnflow is in Kg/m^2/dt |
---|
1161 | ! |
---|
1162 | delsurfstor(:) = -flood_diag(:)-pond_diag(:)-lake_diag(:) |
---|
1163 | hydrographs(:) = zero |
---|
1164 | slowflow_diag(:) = zero |
---|
1165 | flood_diag(:) = zero |
---|
1166 | pond_diag(:) = zero |
---|
1167 | irrigation(:) = zero |
---|
1168 | |
---|
1169 | |
---|
1170 | END SUBROUTINE routing_simple_flow |
---|
1171 | |
---|
1172 | |
---|
1173 | !! ============================================================================================================================= |
---|
1174 | !! SUBROUTINE: routing_simple_initialize |
---|
1175 | !! |
---|
1176 | !>\BRIEF Initialize the routing_simple module |
---|
1177 | !! |
---|
1178 | !! DESCRIPTION: Initialize the routing_simple module. Read from restart file or read the routing.nc file to initialize the |
---|
1179 | !! routing scheme. |
---|
1180 | !! |
---|
1181 | !! RECENT CHANGE(S) |
---|
1182 | !! |
---|
1183 | !! REFERENCE(S) |
---|
1184 | !! |
---|
1185 | !! FLOWCHART |
---|
1186 | !! \n |
---|
1187 | !_ ============================================================================================================================== |
---|
1188 | |
---|
1189 | SUBROUTINE routing_simple_initialize( kjit, nbpt, index, & |
---|
1190 | rest_id, hist_id, hist2_id, lalo, & |
---|
1191 | neighbours, resolution, contfrac, stempdiag, & |
---|
1192 | returnflow, reinfiltration, irrigation, riverflow, & |
---|
1193 | coastalflow, flood_frac, flood_res ) |
---|
1194 | |
---|
1195 | IMPLICIT NONE |
---|
1196 | |
---|
1197 | !! 0 Variable and parameter description |
---|
1198 | !! 0.1 Input variables |
---|
1199 | INTEGER(i_std), INTENT(in) :: kjit !! Time step number (unitless) |
---|
1200 | INTEGER(i_std), INTENT(in) :: nbpt !! Domain size (unitless) |
---|
1201 | INTEGER(i_std), INTENT(in) :: index(nbpt) !! Indices of the points on the map (unitless) |
---|
1202 | INTEGER(i_std),INTENT(in) :: rest_id !! Restart file identifier (unitless) |
---|
1203 | INTEGER(i_std),INTENT(in) :: hist_id !! Access to history file (unitless) |
---|
1204 | INTEGER(i_std),INTENT(in) :: hist2_id !! Access to history file 2 (unitless) |
---|
1205 | REAL(r_std), INTENT(in) :: lalo(nbpt,2) !! Vector of latitude and longitudes (beware of the order !) |
---|
1206 | |
---|
1207 | INTEGER(i_std), INTENT(in) :: neighbours(nbpt,8) !! Vector of neighbours for each grid point |
---|
1208 | !! (1=N, 2=NE, 3=E, 4=SE, 5=S, 6=SW, 7=W, 8=NW) (unitless) |
---|
1209 | REAL(r_std), INTENT(in) :: resolution(nbpt,2) !! The size of each grid box in X and Y (m) |
---|
1210 | REAL(r_std), INTENT(in) :: contfrac(nbpt) !! Fraction of land in each grid box (unitless;0-1) |
---|
1211 | REAL(r_std), INTENT(in) :: stempdiag(nbpt,nslm) !! Diagnostic soil temperature profile |
---|
1212 | |
---|
1213 | !! 0.2 Output variables |
---|
1214 | REAL(r_std), INTENT(out) :: returnflow(nbpt) !! The water flow from lakes and swamps which returns to the grid box. |
---|
1215 | !! This water will go back into the hydrol or hydrolc module to allow re-evaporation (kg/m^2/dt) |
---|
1216 | REAL(r_std), INTENT(out) :: reinfiltration(nbpt) !! Water flow from ponds and floodplains which returns to the grid box (kg/m^2/dt) |
---|
1217 | REAL(r_std), INTENT(out) :: irrigation(nbpt) !! Irrigation flux. This is the water taken from the reservoirs and beeing put into the upper layers of the soil (kg/m^2/dt) |
---|
1218 | 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) |
---|
1219 | |
---|
1220 | 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) |
---|
1221 | REAL(r_std), INTENT(out) :: flood_frac(nbpt) !! Flooded fraction of the grid box (unitless;0-1) |
---|
1222 | REAL(r_std), INTENT(out) :: flood_res(nbpt) !! Diagnostic of water amount in the floodplains reservoir (kg) |
---|
1223 | |
---|
1224 | !! 0.3 Local variables |
---|
1225 | LOGICAL :: init_irrig !! Logical to initialize the irrigation (true/false) |
---|
1226 | LOGICAL :: init_flood !! Logical to initialize the floodplains (true/false) |
---|
1227 | LOGICAL :: init_swamp !! Logical to initialize the swamps (true/false) |
---|
1228 | |
---|
1229 | !_ ================================================================================================================================ |
---|
1230 | |
---|
1231 | ! |
---|
1232 | ! do initialisation |
---|
1233 | ! |
---|
1234 | nbvmax = 440 |
---|
1235 | ! Here we will allocate the memory and get the fixed fields from the restart file. |
---|
1236 | ! If the info is not found then we will compute the routing map. |
---|
1237 | ! |
---|
1238 | CALL routing_simple_init_1 (kjit, nbpt, index, returnflow, reinfiltration, irrigation, & |
---|
1239 | riverflow, coastalflow, flood_frac, flood_res, stempdiag, rest_id) |
---|
1240 | |
---|
1241 | routing_area => routing_area_loc |
---|
1242 | topo_resid => topo_resid_loc |
---|
1243 | route_togrid => route_togrid_loc |
---|
1244 | route_tobasin => route_tobasin_loc |
---|
1245 | global_basinid => global_basinid_loc |
---|
1246 | hydrodiag => hydrodiag_loc |
---|
1247 | |
---|
1248 | ! This routine computes the routing map if the route_togrid_glo is undefined. This means that the |
---|
1249 | ! map has not been initialized during the restart process.. |
---|
1250 | ! |
---|
1251 | !! Reads in the map of the basins and flow directions to construct the catchments of each grid box |
---|
1252 | ! |
---|
1253 | IF ( COUNT(route_togrid_glo .GE. undef_int) .GT. 0 ) THEN |
---|
1254 | ! CALL routing_basins_p(nbpt, lalo, neighbours, resolution, contfrac) |
---|
1255 | ENDIF |
---|
1256 | ! |
---|
1257 | ! Do we have what we need if we want to do irrigation |
---|
1258 | !! Initialisation of flags for irrigated land, flood plains and swamps |
---|
1259 | ! |
---|
1260 | init_irrig = .FALSE. |
---|
1261 | IF ( do_irrigation ) THEN |
---|
1262 | IF (COUNT(irrigated .GE. undef_sechiba-1) > 0) init_irrig = .TRUE. |
---|
1263 | END IF |
---|
1264 | |
---|
1265 | init_flood = .FALSE. |
---|
1266 | IF ( do_floodplains ) THEN |
---|
1267 | IF (COUNT(floodplains .GE. undef_sechiba-1) > 0) init_flood = .TRUE. |
---|
1268 | END IF |
---|
1269 | |
---|
1270 | init_swamp = .FALSE. |
---|
1271 | IF ( doswamps ) THEN |
---|
1272 | IF (COUNT(swamp .GE. undef_sechiba-1) > 0 ) init_swamp = .TRUE. |
---|
1273 | END IF |
---|
1274 | |
---|
1275 | !! If we have irrigated land, flood plains or swamps then we need to interpolate the 0.5 degree |
---|
1276 | !! base data set to the resolution of the model. |
---|
1277 | |
---|
1278 | IF ( init_irrig .OR. init_flood .OR. init_swamp ) THEN |
---|
1279 | ! CALL routing_irrigmap(nbpt, index, lalo, neighbours, resolution, & |
---|
1280 | ! contfrac, init_irrig, irrigated, init_flood, floodplains, init_swamp, swamp, hist_id, hist2_id) |
---|
1281 | ENDIF |
---|
1282 | |
---|
1283 | IF ( do_irrigation ) THEN |
---|
1284 | CALL xios_orchidee_send_field("irrigmap",irrigated) |
---|
1285 | |
---|
1286 | WRITE(numout,*) 'Verification : range of irrigated : ', MINVAL(irrigated), MAXVAL(irrigated) |
---|
1287 | IF ( .NOT. almaoutput ) THEN |
---|
1288 | CALL histwrite_p(hist_id, 'irrigmap', 1, irrigated, nbpt, index) |
---|
1289 | ELSE |
---|
1290 | CALL histwrite_p(hist_id, 'IrrigationMap', 1, irrigated, nbpt, index) |
---|
1291 | ENDIF |
---|
1292 | IF ( hist2_id > 0 ) THEN |
---|
1293 | IF ( .NOT. almaoutput ) THEN |
---|
1294 | CALL histwrite_p(hist2_id, 'irrigmap', 1, irrigated, nbpt, index) |
---|
1295 | ELSE |
---|
1296 | CALL histwrite_p(hist2_id, 'IrrigationMap', 1, irrigated, nbpt, index) |
---|
1297 | ENDIF |
---|
1298 | ENDIF |
---|
1299 | ENDIF |
---|
1300 | |
---|
1301 | IF ( do_floodplains ) THEN |
---|
1302 | CALL xios_orchidee_send_field("floodmap",floodplains) |
---|
1303 | |
---|
1304 | WRITE(numout,*) 'Verification : range of floodplains : ', MINVAL(floodplains), MAXVAL(floodplains) |
---|
1305 | IF ( .NOT. almaoutput ) THEN |
---|
1306 | CALL histwrite_p(hist_id, 'floodmap', 1, floodplains, nbpt, index) |
---|
1307 | ELSE |
---|
1308 | CALL histwrite_p(hist_id, 'FloodplainsMap', 1, floodplains, nbpt, index) |
---|
1309 | ENDIF |
---|
1310 | IF ( hist2_id > 0 ) THEN |
---|
1311 | IF ( .NOT. almaoutput ) THEN |
---|
1312 | CALL histwrite_p(hist2_id, 'floodmap', 1, floodplains, nbpt, index) |
---|
1313 | ELSE |
---|
1314 | CALL histwrite_p(hist2_id, 'FloodplainsMap', 1, floodplains, nbpt, index) |
---|
1315 | ENDIF |
---|
1316 | ENDIF |
---|
1317 | ENDIF |
---|
1318 | |
---|
1319 | IF ( doswamps ) THEN |
---|
1320 | CALL xios_orchidee_send_field("swampmap",swamp) |
---|
1321 | |
---|
1322 | WRITE(numout,*) 'Verification : range of swamp : ', MINVAL(swamp), MAXVAL(swamp) |
---|
1323 | IF ( .NOT. almaoutput ) THEN |
---|
1324 | CALL histwrite_p(hist_id, 'swampmap', 1, swamp, nbpt, index) |
---|
1325 | ELSE |
---|
1326 | CALL histwrite_p(hist_id, 'SwampMap', 1, swamp, nbpt, index) |
---|
1327 | ENDIF |
---|
1328 | IF ( hist2_id > 0 ) THEN |
---|
1329 | IF ( .NOT. almaoutput ) THEN |
---|
1330 | CALL histwrite_p(hist2_id, 'swampmap', 1, swamp, nbpt, index) |
---|
1331 | ELSE |
---|
1332 | CALL histwrite_p(hist2_id, 'SwampMap', 1, swamp, nbpt, index) |
---|
1333 | ENDIF |
---|
1334 | ENDIF |
---|
1335 | ENDIF |
---|
1336 | |
---|
1337 | !! This routine gives a diagnostic of the basins used. |
---|
1338 | ! CALL routing_diagnostic_p(nbpt, index, lalo, resolution, contfrac, hist_id, hist2_id) |
---|
1339 | |
---|
1340 | CALL routing_simple_init_2(nbpt, contfrac) |
---|
1341 | |
---|
1342 | END SUBROUTINE routing_simple_initialize |
---|
1343 | |
---|
1344 | |
---|
1345 | !! ================================================================================================================================ |
---|
1346 | !! SUBROUTINE : routing_simple_main |
---|
1347 | !! |
---|
1348 | !>\BRIEF This module routes the water over the continents (runoff and |
---|
1349 | !! drainage produced by the hydrolc or hydrol module) into the oceans. |
---|
1350 | !! |
---|
1351 | !! DESCRIPTION (definitions, functional, design, flags): |
---|
1352 | !! The routing scheme (Polcher, 2003) carries the water from the runoff and drainage simulated by SECHIBA |
---|
1353 | !! to the ocean through reservoirs, with some delay. The routing scheme is based on |
---|
1354 | !! a parametrization of the water flow on a global scale (Miller et al., 1994; Hagemann |
---|
1355 | !! and Dumenil, 1998). Given the global map of the main watersheds (Oki et al., 1999; |
---|
1356 | !! Fekete et al., 1999; Vorosmarty et al., 2000) which delineates the boundaries of subbasins |
---|
1357 | !! and gives the eight possible directions of water flow within the pixel, the surface |
---|
1358 | !! runoff and the deep drainage are routed to the ocean. The time-step of the routing is one day. |
---|
1359 | !! The scheme also diagnoses how much water is retained in the foodplains and thus return to soil |
---|
1360 | !! moisture or is taken out of the rivers for irrigation. \n |
---|
1361 | !! |
---|
1362 | !! RECENT CHANGE(S): None |
---|
1363 | !! |
---|
1364 | !! MAIN OUTPUT VARIABLE(S): |
---|
1365 | !! The result of the routing are 3 fluxes : |
---|
1366 | !! - riverflow : The water which flows out from the major rivers. The flux will be located |
---|
1367 | !! on the continental grid but this should be a coastal point. |
---|
1368 | !! - coastalflow : This is the water which flows in a disperse way into the ocean. Essentially these |
---|
1369 | !! are the outflows from all of the small rivers. |
---|
1370 | !! - returnflow : This is the water which flows into a land-point - typically rivers which end in |
---|
1371 | !! the desert. This water will go back into the hydrol module to allow re-evaporation. |
---|
1372 | !! - irrigation : This is water taken from the reservoir and is being put into the upper |
---|
1373 | !! layers of the soil. |
---|
1374 | !! The two first fluxes are in kg/dt and the last two fluxes are in kg/(m^2dt).\n |
---|
1375 | !! |
---|
1376 | !! REFERENCE(S) : |
---|
1377 | !! - Miller JR, Russell GL, Caliri G (1994) |
---|
1378 | !! Continental-scale river flow in climate models. |
---|
1379 | !! J. Clim., 7:914-928 |
---|
1380 | !! - Hagemann S and Dumenil L. (1998) |
---|
1381 | !! A parametrization of the lateral waterflow for the global scale. |
---|
1382 | !! Clim. Dyn., 14:17-31 |
---|
1383 | !! - Oki, T., T. Nishimura, and P. Dirmeyer (1999) |
---|
1384 | !! Assessment of annual runoff from land surface models using total runoff integrating pathways (TRIP) |
---|
1385 | !! J. Meteorol. Soc. Jpn., 77, 235-255 |
---|
1386 | !! - Fekete BM, Charles V, Grabs W (2000) |
---|
1387 | !! Global, composite runoff fields based on observed river discharge and simulated water balances. |
---|
1388 | !! Technical report, UNH/GRDC, Global Runoff Data Centre, Koblenz |
---|
1389 | !! - Vorosmarty, C., B. Fekete, B. Meybeck, and R. Lammers (2000) |
---|
1390 | !! Global system of rivers: Its role in organizing continental land mass and defining land-to-ocean linkages |
---|
1391 | !! Global Biogeochem. Cycles, 14, 599-621 |
---|
1392 | !! - Vivant, A-C. (?? 2002) |
---|
1393 | !! Développement du schéma de routage et des plaines d'inondation, MSc Thesis, Paris VI University |
---|
1394 | !! - J. Polcher (2003) |
---|
1395 | !! Les processus de surface a l'echelle globale et leurs interactions avec l'atmosphere |
---|
1396 | !! Habilitation a diriger les recherches, Paris VI University, 67pp. |
---|
1397 | !! |
---|
1398 | !! FLOWCHART : |
---|
1399 | !! \latexonly |
---|
1400 | !! \includegraphics[scale=0.75]{routing_main_flowchart.png} |
---|
1401 | !! \endlatexonly |
---|
1402 | !! \n |
---|
1403 | !_ ================================================================================================================================ |
---|
1404 | |
---|
1405 | |
---|
1406 | SUBROUTINE routing_simple_main(kjit, nbpt, index, & |
---|
1407 | lalo, neighbours, resolution, contfrac, totfrac_nobio, veget_max, floodout, runoff, & |
---|
1408 | drainage, transpot, precip_rain, humrel, k_litt, flood_frac, flood_res, & |
---|
1409 | stempdiag, reinf_slope, returnflow, reinfiltration, irrigation, riverflow, coastalflow, & |
---|
1410 | rest_id, hist_id, hist2_id) |
---|
1411 | |
---|
1412 | IMPLICIT NONE |
---|
1413 | |
---|
1414 | !! 0 Variable and parameter description |
---|
1415 | !! 0.1 Input variables |
---|
1416 | |
---|
1417 | INTEGER(i_std), INTENT(in) :: kjit !! Time step number (unitless) |
---|
1418 | INTEGER(i_std), INTENT(in) :: nbpt !! Domain size (unitless) |
---|
1419 | INTEGER(i_std),INTENT(in) :: rest_id !! Restart file identifier (unitless) |
---|
1420 | INTEGER(i_std),INTENT(in) :: hist_id !! Access to history file (unitless) |
---|
1421 | INTEGER(i_std),INTENT(in) :: hist2_id !! Access to history file 2 (unitless) |
---|
1422 | INTEGER(i_std), INTENT(in) :: index(nbpt) !! Indices of the points on the map (unitless) |
---|
1423 | REAL(r_std), INTENT(in) :: lalo(nbpt,2) !! Vector of latitude and longitudes (beware of the order !) |
---|
1424 | INTEGER(i_std), INTENT(in) :: neighbours(nbpt,8) !! Vector of neighbours for each grid point (1=N, 2=NE, 3=E, 4=SE, 5=S, 6=SW, 7=W, 8=NW) (unitless) |
---|
1425 | REAL(r_std), INTENT(in) :: resolution(nbpt,2) !! The size of each grid box in X and Y (m) |
---|
1426 | REAL(r_std), INTENT(in) :: contfrac(nbpt) !! Fraction of land in each grid box (unitless;0-1) |
---|
1427 | REAL(r_std), INTENT(in) :: totfrac_nobio(nbpt) !! Total fraction of no-vegetation (continental ice, lakes ...) (unitless;0-1) |
---|
1428 | REAL(r_std), INTENT(in) :: veget_max(nbpt,nvm) !! Maximal fraction of vegetation (unitless;0-1) |
---|
1429 | REAL(r_std), INTENT(in) :: floodout(nbpt) !! Grid-point flow out of floodplains (kg/m^2/dt) |
---|
1430 | REAL(r_std), INTENT(in) :: runoff(nbpt) !! Grid-point runoff (kg/m^2/dt) |
---|
1431 | REAL(r_std), INTENT(in) :: drainage(nbpt) !! Grid-point drainage (kg/m^2/dt) |
---|
1432 | REAL(r_std), INTENT(in) :: transpot(nbpt,nvm) !! Potential transpiration of the vegetation (kg/m^2/dt) |
---|
1433 | REAL(r_std), INTENT(in) :: precip_rain(nbpt) !! Rainfall (kg/m^2/dt) |
---|
1434 | REAL(r_std), INTENT(in) :: k_litt(nbpt) !! Averaged conductivity for saturated infiltration in the 'litter' layer (kg/m^2/dt) |
---|
1435 | REAL(r_std), INTENT(in) :: humrel(nbpt,nvm) !! Soil moisture stress, root extraction potential (unitless) |
---|
1436 | REAL(r_std), INTENT(in) :: stempdiag(nbpt,nslm) !! Diagnostic soil temperature profile |
---|
1437 | REAL(r_std), INTENT(in) :: reinf_slope(nbpt) !! Coefficient which determines the reinfiltration ratio in the grid box due to flat areas (unitless;0-1) |
---|
1438 | |
---|
1439 | !! 0.2 Output variables |
---|
1440 | REAL(r_std), INTENT(out) :: returnflow(nbpt) !! The water flow from lakes and swamps which returns to the grid box. |
---|
1441 | !! This water will go back into the hydrol or hydrolc module to allow re-evaporation (kg/m^2/dt) |
---|
1442 | REAL(r_std), INTENT(out) :: reinfiltration(nbpt) !! Water flow from ponds and floodplains which returns to the grid box (kg/m^2/dt) |
---|
1443 | REAL(r_std), INTENT(out) :: irrigation(nbpt) !! Irrigation flux. This is the water taken from the reservoirs and beeing put into the upper layers of the soil (kg/m^2/dt) |
---|
1444 | 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) |
---|
1445 | 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) |
---|
1446 | REAL(r_std), INTENT(out) :: flood_frac(nbpt) !! Flooded fraction of the grid box (unitless;0-1) |
---|
1447 | REAL(r_std), INTENT(out) :: flood_res(nbpt) !! Diagnostic of water amount in the floodplains reservoir (kg) |
---|
1448 | |
---|
1449 | !! 0.3 Local variables |
---|
1450 | CHARACTER(LEN=30) :: var_name !! To store variables names for I/O (unitless) |
---|
1451 | REAL(r_std), DIMENSION(1) :: tmp_day !! |
---|
1452 | REAL(r_std), DIMENSION(nbpt) :: return_lakes !! Water from lakes flowing back into soil moisture (kg/m^2/dt) |
---|
1453 | INTEGER(i_std) :: ig, jv !! Indices (unitless) |
---|
1454 | REAL(r_std), DIMENSION(nbpt) :: tot_vegfrac_nowoody !! Total fraction occupied by grass (0-1,unitless) |
---|
1455 | REAL(r_std),PARAMETER :: delta_lon=(360./144.)/2. |
---|
1456 | REAL(r_std),PARAMETER :: delta_lat=(180./144.)/2. |
---|
1457 | |
---|
1458 | !_ ================================================================================================================================ |
---|
1459 | |
---|
1460 | ! |
---|
1461 | !! Computes the variables averaged between routing time steps and which will be used in subsequent calculations |
---|
1462 | ! |
---|
1463 | floodout_mean(:) = floodout_mean(:) + floodout(:) |
---|
1464 | runoff_mean(:) = runoff_mean(:) + runoff(:) |
---|
1465 | drainage_mean(:) = drainage_mean(:) + drainage(:) |
---|
1466 | floodtemp(:) = stempdiag(:,floodtemp_lev) |
---|
1467 | precip_mean(:) = precip_mean(:) + precip_rain(:) |
---|
1468 | ! |
---|
1469 | !! Computes the total fraction occupied by the grasses and the crops for each grid cell |
---|
1470 | tot_vegfrac_nowoody(:) = zero |
---|
1471 | DO jv = 1, nvm |
---|
1472 | IF ( (jv /= ibare_sechiba) .AND. .NOT.(is_tree(jv)) ) THEN |
---|
1473 | tot_vegfrac_nowoody(:) = tot_vegfrac_nowoody(:) + veget_max(:,jv) |
---|
1474 | END IF |
---|
1475 | END DO |
---|
1476 | |
---|
1477 | DO ig = 1, nbpt |
---|
1478 | IF ( tot_vegfrac_nowoody(ig) .GT. min_sechiba ) THEN |
---|
1479 | DO jv = 1,nvm |
---|
1480 | IF ( (jv /= ibare_sechiba) .AND. .NOT.(is_tree(jv)) ) THEN |
---|
1481 | transpot_mean(ig) = transpot_mean(ig) + transpot(ig,jv) * veget_max(ig,jv)/tot_vegfrac_nowoody(ig) |
---|
1482 | END IF |
---|
1483 | END DO |
---|
1484 | ELSE |
---|
1485 | IF (MAXVAL(veget_max(ig,2:nvm)) .GT. min_sechiba) THEN |
---|
1486 | DO jv = 2, nvm |
---|
1487 | transpot_mean(ig) = transpot_mean(ig) + transpot(ig,jv) * veget_max(ig,jv)/ SUM(veget_max(ig,2:nvm)) |
---|
1488 | ENDDO |
---|
1489 | ENDIF |
---|
1490 | ENDIF |
---|
1491 | ENDDO |
---|
1492 | |
---|
1493 | ! |
---|
1494 | ! Averaged variables (i.e. *dt_sechiba/dt_routing). This accounts for the difference between the shorter |
---|
1495 | ! timestep dt_sechiba of other parts of the model and the long dt_routing timestep (set to one day at present) |
---|
1496 | ! |
---|
1497 | totnobio_mean(:) = totnobio_mean(:) + totfrac_nobio(:)*dt_sechiba/dt_routing |
---|
1498 | k_litt_mean(:) = k_litt_mean(:) + k_litt(:)*dt_sechiba/dt_routing |
---|
1499 | ! |
---|
1500 | ! Only potentially vegetated surfaces are taken into account. At the start of |
---|
1501 | ! the growing seasons we will give more weight to these areas. |
---|
1502 | ! |
---|
1503 | DO jv=2,nvm |
---|
1504 | DO ig=1,nbpt |
---|
1505 | humrel_mean(ig) = humrel_mean(ig) + humrel(ig,jv)*veget_max(ig,jv)*dt_sechiba/dt_routing |
---|
1506 | vegtot_mean(ig) = vegtot_mean(ig) + veget_max(ig,jv)*dt_sechiba/dt_routing |
---|
1507 | ENDDO |
---|
1508 | ENDDO |
---|
1509 | ! |
---|
1510 | time_counter = time_counter + dt_sechiba |
---|
1511 | ! |
---|
1512 | ! If the time has come we do the routing. |
---|
1513 | ! |
---|
1514 | IF ( NINT(time_counter) .GE. NINT(dt_routing) ) THEN |
---|
1515 | ! |
---|
1516 | ! Check the water balance if needed |
---|
1517 | ! |
---|
1518 | ! IF ( check_waterbal ) THEN |
---|
1519 | ! CALL routing_waterbal(nbpt, .TRUE., floodout_mean, runoff_mean, drainage_mean, returnflow_mean, & |
---|
1520 | ! & reinfiltration_mean, irrigation_mean, riverflow_mean, coastalflow_mean) |
---|
1521 | ! ENDIF |
---|
1522 | ! |
---|
1523 | ! Make sure we do not flood north of 49N as there freezing processes start to play a role and they |
---|
1524 | ! are not yet well treated in ORCHIDEE. |
---|
1525 | ! |
---|
1526 | DO ig=1,nbpt |
---|
1527 | IF ( lalo(ig,1) > 49.0 ) THEN |
---|
1528 | floodtemp(ig) = tp_00 - un |
---|
1529 | ENDIF |
---|
1530 | ENDDO |
---|
1531 | ! |
---|
1532 | !! Computes the transport of water in the various reservoirs |
---|
1533 | ! |
---|
1534 | !ym runoff_mean=0 |
---|
1535 | !ym drainage_mean=0 |
---|
1536 | !ym DO ig=1,nbpt |
---|
1537 | !ym IF ( lalo(ig,1)-delta_lat < 0 .AND. lalo(ig,1)+delta_lat > 0 .AND. lalo(ig,2)-delta_lon < 32. .AND. lalo(ig,2)+delta_lon > 32.) runoff_mean(ig)=one_day/dt_routing |
---|
1538 | !ym ENDDO |
---|
1539 | |
---|
1540 | CALL routing_simple_flow(nbpt, dt_routing, lalo, floodout_mean, runoff_mean, drainage_mean, & |
---|
1541 | & vegtot_mean, totnobio_mean, transpot_mean, precip_mean, humrel_mean, k_litt_mean, floodtemp, reinf_slope, & |
---|
1542 | & lakeinflow_mean, returnflow_mean, reinfiltration_mean, irrigation_mean, riverflow_mean, & |
---|
1543 | & coastalflow_mean, hydrographs, slowflow_diag, flood_frac, flood_res) |
---|
1544 | ! |
---|
1545 | !! Responsible for storing the water in lakes |
---|
1546 | ! |
---|
1547 | CALL routing_simple_lake(nbpt, dt_routing, lakeinflow_mean, humrel_mean, contfrac, return_lakes) |
---|
1548 | ! |
---|
1549 | returnflow_mean(:) = returnflow_mean(:) + return_lakes(:) |
---|
1550 | ! |
---|
1551 | !! Check the water balance in the routing scheme |
---|
1552 | ! |
---|
1553 | ! IF ( check_waterbal ) THEN |
---|
1554 | ! CALL routing_waterbal(nbpt, .FALSE., floodout_mean, runoff_mean, drainage_mean, returnflow_mean, & |
---|
1555 | ! & reinfiltration_mean, irrigation_mean, riverflow_mean, coastalflow_mean) |
---|
1556 | ! ENDIF |
---|
1557 | ! |
---|
1558 | time_counter = zero |
---|
1559 | ! |
---|
1560 | floodout_mean(:) = zero |
---|
1561 | runoff_mean(:) = zero |
---|
1562 | drainage_mean(:) = zero |
---|
1563 | transpot_mean(:) = zero |
---|
1564 | precip_mean(:) = zero |
---|
1565 | ! |
---|
1566 | humrel_mean(:) = zero |
---|
1567 | totnobio_mean(:) = zero |
---|
1568 | k_litt_mean(:) = zero |
---|
1569 | vegtot_mean(:) = zero |
---|
1570 | ! |
---|
1571 | ! Change the units of the routing fluxes from kg/dt_routing into kg/dt_sechiba |
---|
1572 | ! and from m^3/dt_routing into m^3/dt_sechiba |
---|
1573 | ! |
---|
1574 | returnflow_mean(:) = returnflow_mean(:)/dt_routing*dt_sechiba |
---|
1575 | reinfiltration_mean(:) = reinfiltration_mean(:)/dt_routing*dt_sechiba |
---|
1576 | irrigation_mean(:) = irrigation_mean(:)/dt_routing*dt_sechiba |
---|
1577 | irrig_netereq(:) = irrig_netereq(:)/dt_routing*dt_sechiba |
---|
1578 | ! |
---|
1579 | ! Change units as above but at the same time transform the kg/dt_sechiba to m^3/dt_sechiba |
---|
1580 | ! |
---|
1581 | riverflow_mean(:) = riverflow_mean(:)/dt_routing*dt_sechiba/mille |
---|
1582 | coastalflow_mean(:) = coastalflow_mean(:)/dt_routing*dt_sechiba/mille |
---|
1583 | hydrographs(:) = hydrographs(:)/dt_routing*dt_sechiba/mille |
---|
1584 | slowflow_diag(:) = slowflow_diag(:)/dt_routing*dt_sechiba/mille |
---|
1585 | ENDIF |
---|
1586 | |
---|
1587 | ! |
---|
1588 | ! Return the fraction of routed water for this time step. |
---|
1589 | ! |
---|
1590 | returnflow(:) = returnflow_mean(:) |
---|
1591 | reinfiltration(:) = reinfiltration_mean(:) |
---|
1592 | irrigation(:) = irrigation_mean(:) |
---|
1593 | riverflow(:) = riverflow_mean(:) |
---|
1594 | coastalflow(:) = coastalflow_mean(:) |
---|
1595 | |
---|
1596 | ! |
---|
1597 | ! Write diagnostics |
---|
1598 | ! |
---|
1599 | |
---|
1600 | CALL xios_orchidee_send_field("riversret",returnflow*one_day/dt_sechiba) |
---|
1601 | CALL xios_orchidee_send_field("hydrographs",hydrographs/dt_sechiba) |
---|
1602 | CALL xios_orchidee_send_field("slowflow",slowflow_diag/dt_sechiba) ! Qb in m3/s |
---|
1603 | CALL xios_orchidee_send_field("reinfiltration",reinfiltration) |
---|
1604 | CALL xios_orchidee_send_field("fastr",fast_diag) |
---|
1605 | CALL xios_orchidee_send_field("slowr",slow_diag) |
---|
1606 | CALL xios_orchidee_send_field("streamr",stream_diag) |
---|
1607 | CALL xios_orchidee_send_field("laker",lake_diag) |
---|
1608 | CALL xios_orchidee_send_field("pondr",pond_diag) |
---|
1609 | CALL xios_orchidee_send_field("irrigation",irrigation*one_day/dt_sechiba) |
---|
1610 | CALL xios_orchidee_send_field("netirrig",irrig_netereq*one_day/dt_sechiba) |
---|
1611 | CALL xios_orchidee_send_field("floodh",flood_height) |
---|
1612 | CALL xios_orchidee_send_field("floodr",flood_diag) |
---|
1613 | delsurfstor = delsurfstor + flood_diag + pond_diag + lake_diag |
---|
1614 | CALL xios_orchidee_send_field("SurfStor",flood_diag+pond_diag+lake_diag) |
---|
1615 | |
---|
1616 | ! Transform from kg/dt_sechiba into m^3/s |
---|
1617 | CALL xios_orchidee_send_field("hydrographs",hydrographs/mille/dt_sechiba) |
---|
1618 | CALL xios_orchidee_send_field("slowflow",slowflow_diag/mille/dt_sechiba) ! previous id name: Qb |
---|
1619 | CALL xios_orchidee_send_field("coastalflow",coastalflow/dt_sechiba) |
---|
1620 | CALL xios_orchidee_send_field("riverflow",riverflow/dt_sechiba) |
---|
1621 | |
---|
1622 | |
---|
1623 | IF ( .NOT. almaoutput ) THEN |
---|
1624 | ! |
---|
1625 | CALL histwrite_p(hist_id, 'riversret', kjit, returnflow, nbpt, index) |
---|
1626 | IF (do_floodplains .OR. doponds) THEN |
---|
1627 | CALL histwrite_p(hist_id, 'reinfiltration', kjit, reinfiltration, nbpt, index) |
---|
1628 | ENDIF |
---|
1629 | CALL histwrite_p(hist_id, 'hydrographs', kjit, hydrographs, nbpt, index) |
---|
1630 | ! |
---|
1631 | CALL histwrite_p(hist_id, 'fastr', kjit, fast_diag, nbpt, index) |
---|
1632 | CALL histwrite_p(hist_id, 'slowr', kjit, slow_diag, nbpt, index) |
---|
1633 | CALL histwrite_p(hist_id, 'streamr', kjit, stream_diag, nbpt, index) |
---|
1634 | IF ( do_floodplains ) THEN |
---|
1635 | CALL histwrite_p(hist_id, 'floodr', kjit, flood_diag, nbpt, index) |
---|
1636 | CALL histwrite_p(hist_id, 'floodh', kjit, flood_height, nbpt, index) |
---|
1637 | ENDIF |
---|
1638 | CALL histwrite_p(hist_id, 'pondr', kjit, pond_diag, nbpt, index) |
---|
1639 | CALL histwrite_p(hist_id, 'lakevol', kjit, lake_diag, nbpt, index) |
---|
1640 | ! |
---|
1641 | IF ( do_irrigation ) THEN |
---|
1642 | CALL histwrite_p(hist_id, 'irrigation', kjit, irrigation, nbpt, index) |
---|
1643 | CALL histwrite_p(hist_id, 'returnflow', kjit, returnflow, nbpt, index) |
---|
1644 | CALL histwrite_p(hist_id, 'netirrig', kjit, irrig_netereq, nbpt, index) |
---|
1645 | ENDIF |
---|
1646 | ! |
---|
1647 | ELSE |
---|
1648 | ! |
---|
1649 | delsurfstor = delsurfstor + flood_diag + pond_diag + lake_diag |
---|
1650 | CALL histwrite_p(hist_id, 'DelSurfStor', kjit, delsurfstor, nbpt, index) |
---|
1651 | CALL histwrite_p(hist_id, 'SurfStor', kjit, flood_diag+pond_diag+lake_diag, nbpt, index) |
---|
1652 | CALL histwrite_p(hist_id, 'Dis', kjit, hydrographs, nbpt, index) |
---|
1653 | IF ( do_irrigation ) THEN |
---|
1654 | CALL histwrite_p(hist_id, 'Qirrig', kjit, irrigation, nbpt, index) |
---|
1655 | CALL histwrite_p(hist_id, 'Qirrig_req', kjit, irrig_netereq, nbpt, index) |
---|
1656 | ENDIF |
---|
1657 | ! |
---|
1658 | ENDIF |
---|
1659 | IF ( hist2_id > 0 ) THEN |
---|
1660 | IF ( .NOT. almaoutput ) THEN |
---|
1661 | ! |
---|
1662 | CALL histwrite_p(hist2_id, 'riversret', kjit, returnflow, nbpt, index) |
---|
1663 | IF (do_floodplains .OR. doponds) THEN |
---|
1664 | CALL histwrite_p(hist2_id, 'reinfiltration', kjit, reinfiltration, nbpt, index) |
---|
1665 | ENDIF |
---|
1666 | CALL histwrite_p(hist2_id, 'hydrographs', kjit, hydrographs, nbpt, index) |
---|
1667 | ! |
---|
1668 | CALL histwrite_p(hist2_id, 'fastr', kjit, fast_diag, nbpt, index) |
---|
1669 | CALL histwrite_p(hist2_id, 'slowr', kjit, slow_diag, nbpt, index) |
---|
1670 | IF ( do_floodplains ) THEN |
---|
1671 | CALL histwrite_p(hist2_id, 'floodr', kjit, flood_diag, nbpt, index) |
---|
1672 | CALL histwrite_p(hist2_id, 'floodh', kjit, flood_height, nbpt, index) |
---|
1673 | ENDIF |
---|
1674 | CALL histwrite_p(hist2_id, 'pondr', kjit, pond_diag, nbpt, index) |
---|
1675 | CALL histwrite_p(hist2_id, 'streamr', kjit, stream_diag, nbpt, index) |
---|
1676 | CALL histwrite_p(hist2_id, 'lakevol', kjit, lake_diag, nbpt, index) |
---|
1677 | ! |
---|
1678 | IF ( do_irrigation ) THEN |
---|
1679 | CALL histwrite_p(hist2_id, 'irrigation', kjit, irrigation, nbpt, index) |
---|
1680 | CALL histwrite_p(hist2_id, 'returnflow', kjit, returnflow, nbpt, index) |
---|
1681 | CALL histwrite_p(hist2_id, 'netirrig', kjit, irrig_netereq, nbpt, index) |
---|
1682 | ENDIF |
---|
1683 | ! |
---|
1684 | ELSE |
---|
1685 | ! |
---|
1686 | delsurfstor=delsurfstor + flood_diag + pond_diag + lake_diag |
---|
1687 | CALL histwrite_p(hist2_id, 'DelSurfStor', kjit, delsurfstor, nbpt, index) |
---|
1688 | CALL histwrite_p(hist2_id, 'SurfStor', kjit, flood_diag+pond_diag+lake_diag, nbpt, index) |
---|
1689 | CALL histwrite_p(hist2_id, 'Dis', kjit, hydrographs, nbpt, index) |
---|
1690 | ! |
---|
1691 | ENDIF |
---|
1692 | ENDIF |
---|
1693 | ! |
---|
1694 | ! |
---|
1695 | END SUBROUTINE routing_simple_main |
---|
1696 | |
---|
1697 | !! ============================================================================================================================= |
---|
1698 | !! SUBROUTINE: routing_simple_finalize |
---|
1699 | !! |
---|
1700 | !>\BRIEF Write to restart file |
---|
1701 | !! |
---|
1702 | !! DESCRIPTION: Write module variables to restart file |
---|
1703 | !! |
---|
1704 | !! RECENT CHANGE(S) |
---|
1705 | !! |
---|
1706 | !! REFERENCE(S) |
---|
1707 | !! |
---|
1708 | !! FLOWCHART |
---|
1709 | !! \n |
---|
1710 | !_ ============================================================================================================================== |
---|
1711 | |
---|
1712 | SUBROUTINE routing_simple_finalize( kjit, nbpt, rest_id, flood_frac, flood_res ) |
---|
1713 | USE xios |
---|
1714 | IMPLICIT NONE |
---|
1715 | |
---|
1716 | !! 0 Variable and parameter description |
---|
1717 | !! 0.1 Input variables |
---|
1718 | INTEGER(i_std), INTENT(in) :: kjit !! Time step number (unitless) |
---|
1719 | INTEGER(i_std), INTENT(in) :: nbpt !! Domain size (unitless) |
---|
1720 | INTEGER(i_std),INTENT(in) :: rest_id !! Restart file identifier (unitless) |
---|
1721 | REAL(r_std), INTENT(in) :: flood_frac(nbpt) !! Flooded fraction of the grid box (unitless;0-1) |
---|
1722 | REAL(r_std), INTENT(in) :: flood_res(nbpt) !! Diagnostic of water amount in the floodplains reservoir (kg) |
---|
1723 | |
---|
1724 | !! 0.2 Local variables |
---|
1725 | REAL(r_std), DIMENSION(1) :: tmp_day |
---|
1726 | |
---|
1727 | !_ ================================================================================================================================ |
---|
1728 | |
---|
1729 | IF (is_omp_root) THEN |
---|
1730 | CALL xios_send_field("fast_reservoir_restart",fast_reservoir_r) |
---|
1731 | CALL xios_send_field("slow_reservoir_restart",slow_reservoir_r) |
---|
1732 | CALL xios_send_field("stream_reservoir_restart",stream_reservoir_r) |
---|
1733 | ENDIF |
---|
1734 | |
---|
1735 | END SUBROUTINE routing_simple_finalize |
---|
1736 | |
---|
1737 | !! ================================================================================================================================ |
---|
1738 | !! SUBROUTINE : routing_simple_init_1 |
---|
1739 | !! |
---|
1740 | !>\BRIEF This subroutine allocates the memory and get the fixed fields from the restart file. |
---|
1741 | !! |
---|
1742 | !! DESCRIPTION: Privat subroutine to the module routing_simple. This subroutine is called in the begining |
---|
1743 | !! of routing_simple_initialize |
---|
1744 | !! |
---|
1745 | !! RECENT CHANGE(S): None |
---|
1746 | !! |
---|
1747 | !! MAIN OUTPUT VARIABLE(S): |
---|
1748 | !! |
---|
1749 | !! REFERENCES : None |
---|
1750 | !! |
---|
1751 | !! FLOWCHART :None |
---|
1752 | !! \n |
---|
1753 | !_ ================================================================================================================================ |
---|
1754 | |
---|
1755 | SUBROUTINE routing_simple_init_1(kjit, nbpt, index, returnflow, reinfiltration, irrigation, & |
---|
1756 | riverflow, coastalflow, flood_frac, flood_res, stempdiag, rest_id) |
---|
1757 | |
---|
1758 | IMPLICIT NONE |
---|
1759 | |
---|
1760 | !! 0 Variable and parameter description |
---|
1761 | !! 0.1 Input variables |
---|
1762 | |
---|
1763 | INTEGER(i_std), INTENT(in) :: kjit !! Time step number (unitless) |
---|
1764 | INTEGER(i_std), INTENT(in) :: nbpt !! Domain size (unitless) |
---|
1765 | INTEGER(i_std), DIMENSION (nbpt), INTENT(in) :: index !! Indices of the points on the map (unitless) |
---|
1766 | REAL(r_std), DIMENSION(nbpt,nslm),INTENT(in) :: stempdiag !! Temperature profile in soil |
---|
1767 | INTEGER(i_std), INTENT(in) :: rest_id !! Restart file identifier (unitless) |
---|
1768 | |
---|
1769 | !! 0.2 Output variables |
---|
1770 | REAL(r_std), DIMENSION (nbpt),INTENT(out) :: returnflow !! The water flow from lakes and swamps which returns into the grid box. |
---|
1771 | !! This water will go back into the hydrol or hydrolc module to allow re-evaporation (kg/m^2/dt) |
---|
1772 | REAL(r_std), DIMENSION (nbpt),INTENT(out) :: reinfiltration !! Water flow from ponds and floodplains which returns to the grid box (kg/m^2/dt) |
---|
1773 | REAL(r_std), DIMENSION (nbpt),INTENT(out) :: irrigation !! Irrigation flux. This is the water taken from the reservoirs and beeing put into the upper layers of the soil.(kg/m^2/dt) |
---|
1774 | REAL(r_std), DIMENSION (nbpt),INTENT(out) :: riverflow !! Outflow of the major rivers. The flux will be located on the continental grid but this should be a coastal point (kg/dt) |
---|
1775 | REAL(r_std), DIMENSION (nbpt),INTENT(out) :: coastalflow !! Outflow on coastal points by small basins. This is the water which flows in a disperse way into the ocean (kg/dt) |
---|
1776 | REAL(r_std), DIMENSION (nbpt),INTENT(out) :: flood_frac !! Flooded fraction of the grid box (unitless;0-1) |
---|
1777 | REAL(r_std), DIMENSION (nbpt),INTENT(out) :: flood_res !! Diagnostic of water amount in the floodplains reservoir (kg) |
---|
1778 | |
---|
1779 | !! 0.3 Local variables |
---|
1780 | CHARACTER(LEN=80) :: var_name !! To store variables names for I/O (unitless) |
---|
1781 | REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: tmp_real_g !! A temporary real array for the integers |
---|
1782 | REAL(r_std), DIMENSION(1) :: tmp_day !! |
---|
1783 | REAL(r_std) :: ratio !! Diagnostic ratio to check that dt_routing is a multiple of dt_sechiba (unitless) |
---|
1784 | REAL(r_std) :: totarea !! Total area of basin (m^2) |
---|
1785 | INTEGER(i_std) :: ier, ig, ib, ipn(1) !! Indices (unitless) |
---|
1786 | |
---|
1787 | !_ ================================================================================================================================ |
---|
1788 | ! |
---|
1789 | ! |
---|
1790 | ! These variables will require the configuration infrastructure |
---|
1791 | ! |
---|
1792 | !Config Key = DT_ROUTING |
---|
1793 | !Config If = RIVER_ROUTING |
---|
1794 | !Config Desc = Time step of the routing scheme |
---|
1795 | !Config Def = one_day |
---|
1796 | !Config Help = This values gives the time step in seconds of the routing scheme. |
---|
1797 | !Config It should be multiple of the main time step of ORCHIDEE. One day |
---|
1798 | !Config is a good value. |
---|
1799 | !Config Units = [seconds] |
---|
1800 | ! |
---|
1801 | dt_routing = one_day |
---|
1802 | CALL getin_p('DT_ROUTING', dt_routing) |
---|
1803 | ! |
---|
1804 | !Config Key = ROUTING_RIVERS |
---|
1805 | !Config If = RIVER_ROUTING |
---|
1806 | !Config Desc = Number of rivers |
---|
1807 | !Config Def = 50 |
---|
1808 | !Config Help = This parameter chooses the number of largest river basins |
---|
1809 | !Config which should be treated as independently as rivers and not |
---|
1810 | !Config flow into the oceans as diffusion coastal flow. |
---|
1811 | !Config Units = [-] |
---|
1812 | num_largest = 50 |
---|
1813 | CALL getin_p('ROUTING_RIVERS', num_largest) |
---|
1814 | ! |
---|
1815 | !Config Key = DO_FLOODINFILT |
---|
1816 | !Config Desc = Should floodplains reinfiltrate into the soil |
---|
1817 | !Config If = RIVER_ROUTING |
---|
1818 | !Config Def = n |
---|
1819 | !Config Help = This parameters allows the user to ask the model |
---|
1820 | !Config to take into account the flood plains reinfiltration |
---|
1821 | !Config into the soil moisture. It then can go |
---|
1822 | !Config back to the slow and fast reservoirs |
---|
1823 | !Config Units = [FLAG] |
---|
1824 | ! |
---|
1825 | dofloodinfilt = .FALSE. |
---|
1826 | CALL getin_p('DO_FLOODINFILT', dofloodinfilt) |
---|
1827 | ! |
---|
1828 | !Config Key = DO_SWAMPS |
---|
1829 | !Config Desc = Should we include swamp parameterization |
---|
1830 | !Config If = RIVER_ROUTING |
---|
1831 | !Config Def = n |
---|
1832 | !Config Help = This parameters allows the user to ask the model |
---|
1833 | !Config to take into account the swamps and return |
---|
1834 | !Config the water into the bottom of the soil. It then can go |
---|
1835 | !Config back to the atmopshere. This tried to simulate |
---|
1836 | !Config internal deltas of rivers. |
---|
1837 | !Config Units = [FLAG] |
---|
1838 | ! |
---|
1839 | doswamps = .FALSE. |
---|
1840 | CALL getin_p('DO_SWAMPS', doswamps) |
---|
1841 | ! |
---|
1842 | !Config Key = DO_PONDS |
---|
1843 | !Config Desc = Should we include ponds |
---|
1844 | !Config If = RIVER_ROUTING |
---|
1845 | !Config Def = n |
---|
1846 | !Config Help = This parameters allows the user to ask the model |
---|
1847 | !Config to take into account the ponds and return |
---|
1848 | !Config the water into the soil moisture. It then can go |
---|
1849 | !Config back to the atmopshere. This tried to simulate |
---|
1850 | !Config little ponds especially in West Africa. |
---|
1851 | !Config Units = [FLAG] |
---|
1852 | ! |
---|
1853 | doponds = .FALSE. |
---|
1854 | CALL getin_p('DO_PONDS', doponds) |
---|
1855 | |
---|
1856 | !Config Key = SLOW_TCST |
---|
1857 | !Config Desc = Time constant for the slow reservoir |
---|
1858 | !Config If = RIVER_ROUTING |
---|
1859 | !Config Def = 25.0 |
---|
1860 | !Config Help = This parameters allows the user to fix the |
---|
1861 | !Config time constant (in days) of the slow reservoir |
---|
1862 | !Config in order to get better river flows for |
---|
1863 | !Config particular regions. |
---|
1864 | !Config Units = [days] |
---|
1865 | ! |
---|
1866 | !> A value for property of each reservoir (in day/m) is given to compute a time constant (in day) |
---|
1867 | !> for each reservoir (product of tcst and topo_resid). |
---|
1868 | !> The value of tcst has been calibrated for the three reservoirs over the Senegal river basin only, |
---|
1869 | !> during the 1 degree NCEP Corrected by Cru (NCC) resolution simulations (Ngo-Duc et al., 2005, Ngo-Duc et al., 2006) and |
---|
1870 | !> generalized for all the basins of the world. The "slow reservoir" and the "fast reservoir" |
---|
1871 | !> have the highest value in order to simulate the groundwater. |
---|
1872 | !> The "stream reservoir", which represents all the water of the stream, has the lowest value. |
---|
1873 | !> Those figures are the same for all the basins of the world. |
---|
1874 | CALL getin_p('SLOW_TCST', slow_tcst) |
---|
1875 | ! |
---|
1876 | !Config Key = FAST_TCST |
---|
1877 | !Config Desc = Time constant for the fast reservoir |
---|
1878 | !Config If = RIVER_ROUTING |
---|
1879 | !Config Def = 3.0 |
---|
1880 | !Config Help = This parameters allows the user to fix the |
---|
1881 | !Config time constant (in days) of the fast reservoir |
---|
1882 | !Config in order to get better river flows for |
---|
1883 | !Config particular regions. |
---|
1884 | !Config Units = [days] |
---|
1885 | CALL getin_p('FAST_TCST', fast_tcst) |
---|
1886 | ! |
---|
1887 | !Config Key = STREAM_TCST |
---|
1888 | !Config Desc = Time constant for the stream reservoir |
---|
1889 | !Config If = RIVER_ROUTING |
---|
1890 | !Config Def = 0.24 |
---|
1891 | !Config Help = This parameters allows the user to fix the |
---|
1892 | !Config time constant (in days) of the stream reservoir |
---|
1893 | !Config in order to get better river flows for |
---|
1894 | !Config particular regions. |
---|
1895 | !Config Units = [days] |
---|
1896 | CALL getin_p('STREAM_TCST', stream_tcst) |
---|
1897 | ! |
---|
1898 | !Config Key = FLOOD_TCST |
---|
1899 | !Config Desc = Time constant for the flood reservoir |
---|
1900 | !Config If = RIVER_ROUTING |
---|
1901 | !Config Def = 4.0 |
---|
1902 | !Config Help = This parameters allows the user to fix the |
---|
1903 | !Config time constant (in days) of the flood reservoir |
---|
1904 | !Config in order to get better river flows for |
---|
1905 | !Config particular regions. |
---|
1906 | !Config Units = [days] |
---|
1907 | CALL getin_p('FLOOD_TCST', flood_tcst) |
---|
1908 | |
---|
1909 | !Config Key = SWAMP_CST |
---|
1910 | !Config Desc = Fraction of the river that flows back to swamps |
---|
1911 | !Config If = RIVER_ROUTING |
---|
1912 | !Config Def = 0.2 |
---|
1913 | !Config Help = This parameters allows the user to fix the |
---|
1914 | !Config fraction of the river transport |
---|
1915 | !Config that flows to swamps |
---|
1916 | !Config Units = [-] |
---|
1917 | CALL getin_p('SWAMP_CST', swamp_cst) |
---|
1918 | |
---|
1919 | !Config Key = FLOOD_BETA |
---|
1920 | !Config Desc = Parameter to fix the shape of the floodplain |
---|
1921 | !Config If = RIVER_ROUTING |
---|
1922 | !Config Def = 2.0 |
---|
1923 | !Config Help = Parameter to fix the shape of the floodplain |
---|
1924 | !Config (>1 for convex edges, <1 for concave edges) |
---|
1925 | !Config Units = [-] |
---|
1926 | CALL getin_p("FLOOD_BETA", beta) |
---|
1927 | |
---|
1928 | !Config Key = POND_BETAP |
---|
1929 | !Config Desc = Ratio of the basin surface intercepted by ponds and the maximum surface of ponds |
---|
1930 | !Config If = RIVER_ROUTING |
---|
1931 | !Config Def = 0.5 |
---|
1932 | !Config Help = |
---|
1933 | !Config Units = [-] |
---|
1934 | CALL getin_p("POND_BETAP", betap) |
---|
1935 | |
---|
1936 | !Config Key = FLOOD_CRI |
---|
1937 | !Config Desc = Potential height for which all the basin is flooded |
---|
1938 | !Config If = DO_FLOODPLAINS or DO_PONDS |
---|
1939 | !Config Def = 2000. |
---|
1940 | !Config Help = |
---|
1941 | !Config Units = [mm] |
---|
1942 | CALL getin_p("FLOOD_CRI", floodcri) |
---|
1943 | |
---|
1944 | !Config Key = POND_CRI |
---|
1945 | !Config Desc = Potential height for which all the basin is a pond |
---|
1946 | !Config If = DO_FLOODPLAINS or DO_PONDS |
---|
1947 | !Config Def = 2000. |
---|
1948 | !Config Help = |
---|
1949 | !Config Units = [mm] |
---|
1950 | CALL getin_p("POND_CRI", pondcri) |
---|
1951 | |
---|
1952 | |
---|
1953 | ! In order to simplify the time cascade check that dt_routing |
---|
1954 | ! is a multiple of dt_sechiba |
---|
1955 | ! |
---|
1956 | ratio = dt_routing/dt_sechiba |
---|
1957 | IF ( ABS(NINT(ratio) - ratio) .GT. 10*EPSILON(ratio)) THEN |
---|
1958 | WRITE(numout,*) 'WARNING -- WARNING -- WARNING -- WARNING' |
---|
1959 | WRITE(numout,*) "The chosen time step for the routing is not a multiple of the" |
---|
1960 | WRITE(numout,*) "main time step of the model. We will change dt_routing so that" |
---|
1961 | WRITE(numout,*) "this condition os fulfilled" |
---|
1962 | dt_routing = NINT(ratio) * dt_sechiba |
---|
1963 | WRITE(numout,*) 'THE NEW DT_ROUTING IS : ', dt_routing |
---|
1964 | ENDIF |
---|
1965 | ! |
---|
1966 | IF ( dt_routing .LT. dt_sechiba) THEN |
---|
1967 | WRITE(numout,*) 'WARNING -- WARNING -- WARNING -- WARNING' |
---|
1968 | WRITE(numout,*) 'The routing timestep can not be smaller than the one' |
---|
1969 | WRITE(numout,*) 'of the model. We reset its value to the model''s timestep.' |
---|
1970 | WRITE(numout,*) 'The old DT_ROUTING is : ', dt_routing |
---|
1971 | dt_routing = dt_sechiba |
---|
1972 | WRITE(numout,*) 'THE NEW DT_ROUTING IS : ', dt_routing |
---|
1973 | ENDIF |
---|
1974 | ! |
---|
1975 | var_name ="routingcounter" |
---|
1976 | IF (is_root_prc) THEN |
---|
1977 | CALL ioconf_setatt('UNITS', 's') |
---|
1978 | CALL ioconf_setatt('LONG_NAME','Time counter for the routing scheme') |
---|
1979 | CALL restget (rest_id, var_name, 1, 1, 1, kjit, .TRUE., tmp_day) |
---|
1980 | IF (tmp_day(1) == val_exp) THEN |
---|
1981 | time_counter = zero |
---|
1982 | ELSE |
---|
1983 | time_counter = tmp_day(1) |
---|
1984 | ENDIF |
---|
1985 | ENDIF |
---|
1986 | CALL bcast(time_counter) |
---|
1987 | |
---|
1988 | |
---|
1989 | ALLOCATE (routing_area_loc(nbpt,nbasmax), stat=ier) |
---|
1990 | IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for routing_area_loc','','') |
---|
1991 | |
---|
1992 | ALLOCATE (routing_area_glo(nbp_glo,nbasmax)) |
---|
1993 | IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for routing_area_glo','','') |
---|
1994 | var_name = 'routingarea' |
---|
1995 | IF (is_root_prc) THEN |
---|
1996 | CALL ioconf_setatt('UNITS', 'm^2') |
---|
1997 | CALL ioconf_setatt('LONG_NAME','Area of basin') |
---|
1998 | CALL restget (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., routing_area_glo, "gather", nbp_glo, index_g) |
---|
1999 | ENDIF |
---|
2000 | CALL scatter(routing_area_glo,routing_area_loc) |
---|
2001 | routing_area=>routing_area_loc |
---|
2002 | |
---|
2003 | ALLOCATE (tmp_real_g(nbp_glo,nbasmax), stat=ier) |
---|
2004 | IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for tmp_real_g','','') |
---|
2005 | |
---|
2006 | ALLOCATE (route_togrid_loc(nbpt,nbasmax), stat=ier) |
---|
2007 | IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for route_togrid_loc','','') |
---|
2008 | ALLOCATE (route_togrid_glo(nbp_glo,nbasmax), stat=ier) ! used in global in routing_simple_flow |
---|
2009 | IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for route_togrid_glo','','') |
---|
2010 | |
---|
2011 | IF (is_root_prc) THEN |
---|
2012 | var_name = 'routetogrid' |
---|
2013 | CALL ioconf_setatt('UNITS', '-') |
---|
2014 | CALL ioconf_setatt('LONG_NAME','Grid into which the basin flows') |
---|
2015 | CALL restget (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., tmp_real_g, "gather", nbp_glo, index_g) |
---|
2016 | route_togrid_glo(:,:) = undef_int |
---|
2017 | WHERE ( tmp_real_g .LT. val_exp ) |
---|
2018 | route_togrid_glo = NINT(tmp_real_g) |
---|
2019 | ENDWHERE |
---|
2020 | ENDIF |
---|
2021 | CALL bcast(route_togrid_glo) ! used in global in routing_simple_flow |
---|
2022 | CALL scatter(route_togrid_glo,route_togrid_loc) |
---|
2023 | route_togrid=>route_togrid_loc |
---|
2024 | ! |
---|
2025 | ALLOCATE (route_tobasin_loc(nbpt,nbasmax), stat=ier) |
---|
2026 | IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for route_tobasin_loc','','') |
---|
2027 | |
---|
2028 | ALLOCATE (route_tobasin_glo(nbp_glo,nbasmax), stat=ier) |
---|
2029 | IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for route_tobasin_glo','','') |
---|
2030 | |
---|
2031 | IF (is_root_prc) THEN |
---|
2032 | var_name = 'routetobasin' |
---|
2033 | CALL ioconf_setatt('UNITS', '-') |
---|
2034 | CALL ioconf_setatt('LONG_NAME','Basin in to which the water goes') |
---|
2035 | CALL restget (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., tmp_real_g, "gather", nbp_glo, index_g) |
---|
2036 | route_tobasin_glo = undef_int |
---|
2037 | WHERE ( tmp_real_g .LT. val_exp ) |
---|
2038 | route_tobasin_glo = NINT(tmp_real_g) |
---|
2039 | ENDWHERE |
---|
2040 | ENDIF |
---|
2041 | CALL scatter(route_tobasin_glo,route_tobasin_loc) |
---|
2042 | route_tobasin=>route_tobasin_loc |
---|
2043 | ! |
---|
2044 | ! nbintobasin |
---|
2045 | ! |
---|
2046 | ALLOCATE (route_nbintobas_loc(nbpt,nbasmax), stat=ier) |
---|
2047 | IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for route_nbintobas_loc','','') |
---|
2048 | ALLOCATE (route_nbintobas_glo(nbp_glo,nbasmax), stat=ier) |
---|
2049 | IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for route_nbintobas_glo','','') |
---|
2050 | |
---|
2051 | IF (is_root_prc) THEN |
---|
2052 | var_name = 'routenbintobas' |
---|
2053 | CALL ioconf_setatt('UNITS', '-') |
---|
2054 | CALL ioconf_setatt('LONG_NAME','Number of basin into current one') |
---|
2055 | CALL restget (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., tmp_real_g, "gather", nbp_glo, index_g) |
---|
2056 | route_nbintobas_glo = undef_int |
---|
2057 | WHERE ( tmp_real_g .LT. val_exp ) |
---|
2058 | route_nbintobas_glo = NINT(tmp_real_g) |
---|
2059 | ENDWHERE |
---|
2060 | ENDIF |
---|
2061 | CALL scatter(route_nbintobas_glo,route_nbintobas_loc) |
---|
2062 | route_nbintobas=>route_nbintobas_loc |
---|
2063 | ! |
---|
2064 | ALLOCATE (global_basinid_loc(nbpt,nbasmax), stat=ier) |
---|
2065 | IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for global_basinid_loc','','') |
---|
2066 | ALLOCATE (global_basinid_glo(nbp_glo,nbasmax), stat=ier) |
---|
2067 | IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for global_basinid_glo','','') |
---|
2068 | |
---|
2069 | IF (is_root_prc) THEN |
---|
2070 | var_name = 'basinid' |
---|
2071 | CALL ioconf_setatt('UNITS', '-') |
---|
2072 | CALL ioconf_setatt('LONG_NAME','ID of basin') |
---|
2073 | CALL restget (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., tmp_real_g, "gather", nbp_glo, index_g) |
---|
2074 | global_basinid_glo = undef_int |
---|
2075 | WHERE ( tmp_real_g .LT. val_exp ) |
---|
2076 | global_basinid_glo = NINT(tmp_real_g) |
---|
2077 | ENDWHERE |
---|
2078 | ENDIF |
---|
2079 | CALL scatter(global_basinid_glo,global_basinid_loc) |
---|
2080 | global_basinid=>global_basinid_loc |
---|
2081 | ! |
---|
2082 | ALLOCATE (topo_resid_loc(nbpt,nbasmax), stat=ier) |
---|
2083 | IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for topo_resid_loc','','') |
---|
2084 | ALLOCATE (topo_resid_glo(nbp_glo,nbasmax), stat=ier) |
---|
2085 | IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for topo_resid_glo','','') |
---|
2086 | |
---|
2087 | IF (is_root_prc) THEN |
---|
2088 | var_name = 'topoindex' |
---|
2089 | CALL ioconf_setatt('UNITS', 'm') |
---|
2090 | CALL ioconf_setatt('LONG_NAME','Topographic index of the residence time') |
---|
2091 | CALL restget (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., topo_resid_glo, "gather", nbp_glo, index_g) |
---|
2092 | ENDIF |
---|
2093 | CALL scatter(topo_resid_glo,topo_resid_loc) |
---|
2094 | topo_resid=>topo_resid_loc |
---|
2095 | |
---|
2096 | ALLOCATE (fast_reservoir(nbpt,nbasmax), stat=ier) |
---|
2097 | IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for fast_reservoir','','') |
---|
2098 | var_name = 'fastres' |
---|
2099 | CALL ioconf_setatt_p('UNITS', 'Kg') |
---|
2100 | CALL ioconf_setatt_p('LONG_NAME','Water in the fast reservoir') |
---|
2101 | CALL restget_p (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., fast_reservoir, "gather", nbp_glo, index_g) |
---|
2102 | CALL setvar_p (fast_reservoir, val_exp, 'NO_KEYWORD', zero) |
---|
2103 | |
---|
2104 | ALLOCATE (slow_reservoir(nbpt,nbasmax), stat=ier) |
---|
2105 | IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for slow_reservoir','','') |
---|
2106 | var_name = 'slowres' |
---|
2107 | CALL ioconf_setatt_p('UNITS', 'Kg') |
---|
2108 | CALL ioconf_setatt_p('LONG_NAME','Water in the slow reservoir') |
---|
2109 | CALL restget_p (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., slow_reservoir, "gather", nbp_glo, index_g) |
---|
2110 | CALL setvar_p (slow_reservoir, val_exp, 'NO_KEYWORD', zero) |
---|
2111 | |
---|
2112 | ALLOCATE (stream_reservoir(nbpt,nbasmax), stat=ier) |
---|
2113 | IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for stream_reservoir','','') |
---|
2114 | var_name = 'streamres' |
---|
2115 | CALL ioconf_setatt_p('UNITS', 'Kg') |
---|
2116 | CALL ioconf_setatt_p('LONG_NAME','Water in the stream reservoir') |
---|
2117 | CALL restget_p (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., stream_reservoir, "gather", nbp_glo, index_g) |
---|
2118 | CALL setvar_p (stream_reservoir, val_exp, 'NO_KEYWORD', zero) |
---|
2119 | |
---|
2120 | ALLOCATE (flood_reservoir(nbpt,nbasmax), stat=ier) |
---|
2121 | IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for flood_reservoir','','') |
---|
2122 | var_name = 'floodres' |
---|
2123 | CALL ioconf_setatt_p('UNITS', 'Kg') |
---|
2124 | CALL ioconf_setatt_p('LONG_NAME','Water in the flood reservoir') |
---|
2125 | CALL restget_p (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., flood_reservoir, "gather", nbp_glo, index_g) |
---|
2126 | CALL setvar_p (flood_reservoir, val_exp, 'NO_KEYWORD', zero) |
---|
2127 | |
---|
2128 | ALLOCATE (flood_frac_bas(nbpt,nbasmax), stat=ier) |
---|
2129 | IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for flood_frac_bas','','') |
---|
2130 | var_name = 'flood_frac_bas' |
---|
2131 | CALL ioconf_setatt_p('UNITS', '-') |
---|
2132 | CALL ioconf_setatt_p('LONG_NAME','Flooded fraction per basin') |
---|
2133 | CALL restget_p (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., flood_frac_bas, "gather", nbp_glo, index_g) |
---|
2134 | CALL setvar_p (flood_frac_bas, val_exp, 'NO_KEYWORD', zero) |
---|
2135 | |
---|
2136 | ALLOCATE (flood_height(nbpt), stat=ier) |
---|
2137 | IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for flood_height','','') |
---|
2138 | var_name = 'floodh' |
---|
2139 | CALL ioconf_setatt_p('UNITS', '-') |
---|
2140 | CALL ioconf_setatt_p('LONG_NAME','') |
---|
2141 | CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., flood_height, "gather", nbp_glo, index_g) |
---|
2142 | CALL setvar_p (flood_height, val_exp, 'NO_KEYWORD', zero) |
---|
2143 | |
---|
2144 | ALLOCATE (pond_frac(nbpt), stat=ier) |
---|
2145 | IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for pond_frac','','') |
---|
2146 | var_name = 'pond_frac' |
---|
2147 | CALL ioconf_setatt_p('UNITS', '-') |
---|
2148 | CALL ioconf_setatt_p('LONG_NAME','Pond fraction per grid box') |
---|
2149 | CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., pond_frac, "gather", nbp_glo, index_g) |
---|
2150 | CALL setvar_p (pond_frac, val_exp, 'NO_KEYWORD', zero) |
---|
2151 | |
---|
2152 | var_name = 'flood_frac' |
---|
2153 | CALL ioconf_setatt_p('UNITS', '-') |
---|
2154 | CALL ioconf_setatt_p('LONG_NAME','Flooded fraction per grid box') |
---|
2155 | CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., flood_frac, "gather", nbp_glo, index_g) |
---|
2156 | CALL setvar_p (flood_frac, val_exp, 'NO_KEYWORD', zero) |
---|
2157 | |
---|
2158 | var_name = 'flood_res' |
---|
2159 | CALL ioconf_setatt_p('UNITS','mm') |
---|
2160 | CALL ioconf_setatt_p('LONG_NAME','Flooded quantity (estimation)') |
---|
2161 | CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., flood_res, "gather", nbp_glo, index_g) |
---|
2162 | CALL setvar_p (flood_res, val_exp, 'NO_KEYWORD', zero) |
---|
2163 | ! flood_res = zero |
---|
2164 | |
---|
2165 | ALLOCATE (lake_reservoir(nbpt), stat=ier) |
---|
2166 | IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for lake_reservoir','','') |
---|
2167 | var_name = 'lakeres' |
---|
2168 | CALL ioconf_setatt_p('UNITS', 'Kg') |
---|
2169 | CALL ioconf_setatt_p('LONG_NAME','Water in the lake reservoir') |
---|
2170 | CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., lake_reservoir, "gather", nbp_glo, index_g) |
---|
2171 | CALL setvar_p (lake_reservoir, val_exp, 'NO_KEYWORD', zero) |
---|
2172 | |
---|
2173 | ALLOCATE (pond_reservoir(nbpt), stat=ier) |
---|
2174 | IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for pond_reservoir','','') |
---|
2175 | var_name = 'pondres' |
---|
2176 | CALL ioconf_setatt_p('UNITS', 'Kg') |
---|
2177 | CALL ioconf_setatt_p('LONG_NAME','Water in the pond reservoir') |
---|
2178 | CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., pond_reservoir, "gather", nbp_glo, index_g) |
---|
2179 | CALL setvar_p (pond_reservoir, val_exp, 'NO_KEYWORD', zero) |
---|
2180 | ! |
---|
2181 | ! Map of irrigated areas |
---|
2182 | ! |
---|
2183 | IF ( do_irrigation ) THEN |
---|
2184 | ALLOCATE (irrigated(nbpt), stat=ier) |
---|
2185 | IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for irrigated','','') |
---|
2186 | var_name = 'irrigated' |
---|
2187 | CALL ioconf_setatt_p('UNITS', 'm^2') |
---|
2188 | CALL ioconf_setatt_p('LONG_NAME','Surface of irrigated area') |
---|
2189 | CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., irrigated, "gather", nbp_glo, index_g) |
---|
2190 | CALL setvar_p (irrigated, val_exp, 'NO_KEYWORD', undef_sechiba) |
---|
2191 | ENDIF |
---|
2192 | |
---|
2193 | IF ( do_floodplains ) THEN |
---|
2194 | ALLOCATE (floodplains(nbpt), stat=ier) |
---|
2195 | IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for floodplains','','') |
---|
2196 | var_name = 'floodplains' |
---|
2197 | CALL ioconf_setatt_p('UNITS', 'm^2') |
---|
2198 | CALL ioconf_setatt_p('LONG_NAME','Surface which can be flooded') |
---|
2199 | CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., floodplains, "gather", nbp_glo, index_g) |
---|
2200 | CALL setvar_p (floodplains, val_exp, 'NO_KEYWORD', undef_sechiba) |
---|
2201 | ENDIF |
---|
2202 | IF ( doswamps ) THEN |
---|
2203 | ALLOCATE (swamp(nbpt), stat=ier) |
---|
2204 | IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for swamp','','') |
---|
2205 | var_name = 'swamp' |
---|
2206 | CALL ioconf_setatt_p('UNITS', 'm^2') |
---|
2207 | CALL ioconf_setatt_p('LONG_NAME','Surface which can become swamp') |
---|
2208 | CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., swamp, "gather", nbp_glo, index_g) |
---|
2209 | CALL setvar_p (swamp, val_exp, 'NO_KEYWORD', undef_sechiba) |
---|
2210 | ENDIF |
---|
2211 | ! |
---|
2212 | ! Put into the restart file the fluxes so that they can be regenerated at restart. |
---|
2213 | ! |
---|
2214 | ALLOCATE (lakeinflow_mean(nbpt), stat=ier) |
---|
2215 | IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for lakeinflow_mean','','') |
---|
2216 | var_name = 'lakeinflow' |
---|
2217 | CALL ioconf_setatt_p('UNITS', 'Kg/dt') |
---|
2218 | CALL ioconf_setatt_p('LONG_NAME','Lake inflow') |
---|
2219 | CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., lakeinflow_mean, "gather", nbp_glo, index_g) |
---|
2220 | CALL setvar_p (lakeinflow_mean, val_exp, 'NO_KEYWORD', zero) |
---|
2221 | |
---|
2222 | ALLOCATE (returnflow_mean(nbpt), stat=ier) |
---|
2223 | IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for returnflow_mean','','') |
---|
2224 | var_name = 'returnflow' |
---|
2225 | CALL ioconf_setatt_p('UNITS', 'Kg/m^2/dt') |
---|
2226 | CALL ioconf_setatt_p('LONG_NAME','Deep return flux') |
---|
2227 | CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., returnflow_mean, "gather", nbp_glo, index_g) |
---|
2228 | CALL setvar_p (returnflow_mean, val_exp, 'NO_KEYWORD', zero) |
---|
2229 | returnflow(:) = returnflow_mean(:) |
---|
2230 | |
---|
2231 | ALLOCATE (reinfiltration_mean(nbpt), stat=ier) |
---|
2232 | IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for reinfiltration_mean','','') |
---|
2233 | var_name = 'reinfiltration' |
---|
2234 | CALL ioconf_setatt_p('UNITS', 'Kg/m^2/dt') |
---|
2235 | CALL ioconf_setatt_p('LONG_NAME','Top return flux') |
---|
2236 | CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., reinfiltration_mean, "gather", nbp_glo, index_g) |
---|
2237 | CALL setvar_p (reinfiltration_mean, val_exp, 'NO_KEYWORD', zero) |
---|
2238 | reinfiltration(:) = reinfiltration_mean(:) |
---|
2239 | |
---|
2240 | ALLOCATE (irrigation_mean(nbpt), stat=ier) |
---|
2241 | IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for irrigation_mean','','') |
---|
2242 | ALLOCATE (irrig_netereq(nbpt), stat=ier) |
---|
2243 | IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for irrig_netereq','','') |
---|
2244 | irrig_netereq(:) = zero |
---|
2245 | |
---|
2246 | IF ( do_irrigation ) THEN |
---|
2247 | var_name = 'irrigation' |
---|
2248 | CALL ioconf_setatt_p('UNITS', 'Kg/dt') |
---|
2249 | CALL ioconf_setatt_p('LONG_NAME','Artificial irrigation flux') |
---|
2250 | CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., irrigation_mean, "gather", nbp_glo, index_g) |
---|
2251 | CALL setvar_p (irrigation_mean, val_exp, 'NO_KEYWORD', zero) |
---|
2252 | ELSE |
---|
2253 | irrigation_mean(:) = zero |
---|
2254 | ENDIF |
---|
2255 | irrigation(:) = irrigation_mean(:) |
---|
2256 | |
---|
2257 | ALLOCATE (riverflow_mean(nbpt), stat=ier) |
---|
2258 | IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for riverflow_mean','','') |
---|
2259 | var_name = 'riverflow' |
---|
2260 | CALL ioconf_setatt_p('UNITS', 'Kg/dt') |
---|
2261 | CALL ioconf_setatt_p('LONG_NAME','River flux into the sea') |
---|
2262 | CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., riverflow_mean, "gather", nbp_glo, index_g) |
---|
2263 | CALL setvar_p (riverflow_mean, val_exp, 'NO_KEYWORD', zero) |
---|
2264 | riverflow(:) = riverflow_mean(:) |
---|
2265 | |
---|
2266 | ALLOCATE (coastalflow_mean(nbpt), stat=ier) |
---|
2267 | IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for coastalflow_mean','','') |
---|
2268 | var_name = 'coastalflow' |
---|
2269 | CALL ioconf_setatt_p('UNITS', 'Kg/dt') |
---|
2270 | CALL ioconf_setatt_p('LONG_NAME','Diffuse flux into the sea') |
---|
2271 | CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., coastalflow_mean, "gather", nbp_glo, index_g) |
---|
2272 | CALL setvar_p (coastalflow_mean, val_exp, 'NO_KEYWORD', zero) |
---|
2273 | coastalflow(:) = coastalflow_mean(:) |
---|
2274 | |
---|
2275 | ! Locate it at the 2m level |
---|
2276 | ipn = MINLOC(ABS(diaglev-2)) |
---|
2277 | floodtemp_lev = ipn(1) |
---|
2278 | ALLOCATE (floodtemp(nbpt), stat=ier) |
---|
2279 | IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for floodtemp','','') |
---|
2280 | floodtemp(:) = stempdiag(:,floodtemp_lev) |
---|
2281 | |
---|
2282 | ALLOCATE(hydrographs(nbpt), stat=ier) |
---|
2283 | IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for hydrographs','','') |
---|
2284 | var_name = 'hydrographs' |
---|
2285 | CALL ioconf_setatt_p('UNITS', 'm^3/dt') |
---|
2286 | CALL ioconf_setatt_p('LONG_NAME','Hydrograph at outlow of grid') |
---|
2287 | CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., hydrographs, "gather", nbp_glo, index_g) |
---|
2288 | CALL setvar_p (hydrographs, val_exp, 'NO_KEYWORD', zero) |
---|
2289 | |
---|
2290 | ALLOCATE(slowflow_diag(nbpt), stat=ier) |
---|
2291 | IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for slowflow_diag','','') |
---|
2292 | var_name = 'slowflow_diag' |
---|
2293 | CALL ioconf_setatt_p('UNITS', 'm^3/dt') |
---|
2294 | CALL ioconf_setatt_p('LONG_NAME','Slowflow hydrograph at outlow of grid') |
---|
2295 | CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE.,slowflow_diag, "gather", nbp_glo, index_g) |
---|
2296 | CALL setvar_p (slowflow_diag, val_exp, 'NO_KEYWORD', zero) |
---|
2297 | |
---|
2298 | ! |
---|
2299 | ! The diagnostic variables, they are initialized from the above restart variables. |
---|
2300 | ! |
---|
2301 | ALLOCATE(fast_diag(nbpt), slow_diag(nbpt), stream_diag(nbpt), flood_diag(nbpt), & |
---|
2302 | & pond_diag(nbpt), lake_diag(nbpt), delsurfstor(nbpt), stat=ier) |
---|
2303 | IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for fast_diag,..','','') |
---|
2304 | |
---|
2305 | fast_diag(:) = zero |
---|
2306 | slow_diag(:) = zero |
---|
2307 | stream_diag(:) = zero |
---|
2308 | flood_diag(:) = zero |
---|
2309 | pond_diag(:) = zero |
---|
2310 | lake_diag(:) = zero |
---|
2311 | delsurfstor(:) = zero |
---|
2312 | |
---|
2313 | DO ig=1,nbpt |
---|
2314 | totarea = zero |
---|
2315 | DO ib=1,nbasmax |
---|
2316 | totarea = totarea + routing_area(ig,ib) |
---|
2317 | fast_diag(ig) = fast_diag(ig) + fast_reservoir(ig,ib) |
---|
2318 | slow_diag(ig) = slow_diag(ig) + slow_reservoir(ig,ib) |
---|
2319 | stream_diag(ig) = stream_diag(ig) + stream_reservoir(ig,ib) |
---|
2320 | flood_diag(ig) = flood_diag(ig) + flood_reservoir(ig,ib) |
---|
2321 | ENDDO |
---|
2322 | ! |
---|
2323 | fast_diag(ig) = fast_diag(ig)/totarea |
---|
2324 | slow_diag(ig) = slow_diag(ig)/totarea |
---|
2325 | stream_diag(ig) = stream_diag(ig)/totarea |
---|
2326 | flood_diag(ig) = flood_diag(ig)/totarea |
---|
2327 | ! |
---|
2328 | ! This is the volume of the lake scaled to the entire grid. |
---|
2329 | ! It would be batter to scale it to the size of the lake |
---|
2330 | ! but this information is not yet available. |
---|
2331 | ! |
---|
2332 | lake_diag(ig) = lake_reservoir(ig)/totarea |
---|
2333 | ! |
---|
2334 | ENDDO |
---|
2335 | ! |
---|
2336 | ! Get from the restart the fluxes we accumulated. |
---|
2337 | ! |
---|
2338 | ALLOCATE (floodout_mean(nbpt), stat=ier) |
---|
2339 | IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for floodout_mean','','') |
---|
2340 | var_name = 'floodout_route' |
---|
2341 | CALL ioconf_setatt_p('UNITS', 'Kg') |
---|
2342 | CALL ioconf_setatt_p('LONG_NAME','Accumulated flow out of floodplains for routing') |
---|
2343 | CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., floodout_mean, "gather", nbp_glo, index_g) |
---|
2344 | CALL setvar_p (floodout_mean, val_exp, 'NO_KEYWORD', zero) |
---|
2345 | |
---|
2346 | ALLOCATE (runoff_mean(nbpt), stat=ier) |
---|
2347 | IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for runoff_mean','','') |
---|
2348 | var_name = 'runoff_route' |
---|
2349 | CALL ioconf_setatt_p('UNITS', 'Kg') |
---|
2350 | CALL ioconf_setatt_p('LONG_NAME','Accumulated runoff for routing') |
---|
2351 | CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., runoff_mean, "gather", nbp_glo, index_g) |
---|
2352 | CALL setvar_p (runoff_mean, val_exp, 'NO_KEYWORD', zero) |
---|
2353 | |
---|
2354 | ALLOCATE(drainage_mean(nbpt), stat=ier) |
---|
2355 | IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for drainage_mean','','') |
---|
2356 | var_name = 'drainage_route' |
---|
2357 | CALL ioconf_setatt_p('UNITS', 'Kg') |
---|
2358 | CALL ioconf_setatt_p('LONG_NAME','Accumulated drainage for routing') |
---|
2359 | CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., drainage_mean, "gather", nbp_glo, index_g) |
---|
2360 | CALL setvar_p (drainage_mean, val_exp, 'NO_KEYWORD', zero) |
---|
2361 | |
---|
2362 | ALLOCATE(transpot_mean(nbpt), stat=ier) |
---|
2363 | IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for transpot_mean','','') |
---|
2364 | var_name = 'transpot_route' |
---|
2365 | CALL ioconf_setatt_p('UNITS', 'Kg/m^2') |
---|
2366 | CALL ioconf_setatt_p('LONG_NAME','Accumulated potential transpiration for routing/irrigation') |
---|
2367 | CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., transpot_mean, "gather", nbp_glo, index_g) |
---|
2368 | CALL setvar_p (transpot_mean, val_exp, 'NO_KEYWORD', zero) |
---|
2369 | |
---|
2370 | ALLOCATE(precip_mean(nbpt), stat=ier) |
---|
2371 | IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for precip_mean','','') |
---|
2372 | var_name = 'precip_route' |
---|
2373 | CALL ioconf_setatt_p('UNITS', 'Kg/m^2') |
---|
2374 | CALL ioconf_setatt_p('LONG_NAME','Accumulated rain precipitation for irrigation') |
---|
2375 | CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., precip_mean, "gather", nbp_glo, index_g) |
---|
2376 | CALL setvar_p (precip_mean, val_exp, 'NO_KEYWORD', zero) |
---|
2377 | |
---|
2378 | ALLOCATE(humrel_mean(nbpt), stat=ier) |
---|
2379 | IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for humrel_mean','','') |
---|
2380 | var_name = 'humrel_route' |
---|
2381 | CALL ioconf_setatt_p('UNITS', '-') |
---|
2382 | CALL ioconf_setatt_p('LONG_NAME','Mean humrel for irrigation') |
---|
2383 | CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., humrel_mean, "gather", nbp_glo, index_g) |
---|
2384 | CALL setvar_p (humrel_mean, val_exp, 'NO_KEYWORD', un) |
---|
2385 | |
---|
2386 | ALLOCATE(k_litt_mean(nbpt), stat=ier) |
---|
2387 | IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for k_litt_mean','','') |
---|
2388 | var_name = 'k_litt_route' |
---|
2389 | CALL ioconf_setatt_p('UNITS', '-') |
---|
2390 | CALL ioconf_setatt_p('LONG_NAME','Mean cond. for litter') |
---|
2391 | CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., k_litt_mean, "gather", nbp_glo, index_g) |
---|
2392 | CALL setvar_p (k_litt_mean, val_exp, 'NO_KEYWORD', zero) |
---|
2393 | |
---|
2394 | ALLOCATE(totnobio_mean(nbpt), stat=ier) |
---|
2395 | IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for totnobio_mean','','') |
---|
2396 | var_name = 'totnobio_route' |
---|
2397 | CALL ioconf_setatt_p('UNITS', '-') |
---|
2398 | CALL ioconf_setatt_p('LONG_NAME','Last Total fraction of no bio for irrigation') |
---|
2399 | CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., totnobio_mean, "gather", nbp_glo, index_g) |
---|
2400 | CALL setvar_p (totnobio_mean, val_exp, 'NO_KEYWORD', zero) |
---|
2401 | |
---|
2402 | ALLOCATE(vegtot_mean(nbpt), stat=ier) |
---|
2403 | IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for vegtot_mean','','') |
---|
2404 | var_name = 'vegtot_route' |
---|
2405 | CALL ioconf_setatt_p('UNITS', '-') |
---|
2406 | CALL ioconf_setatt_p('LONG_NAME','Last Total fraction of vegetation') |
---|
2407 | CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., vegtot_mean, "gather", nbp_glo, index_g) |
---|
2408 | CALL setvar_p (vegtot_mean, val_exp, 'NO_KEYWORD', un) |
---|
2409 | ! |
---|
2410 | ! |
---|
2411 | DEALLOCATE(tmp_real_g) |
---|
2412 | ! |
---|
2413 | ! Allocate diagnostic variables |
---|
2414 | ! |
---|
2415 | ALLOCATE(hydrodiag_loc(nbpt,nbasmax),hydrodiag_glo(nbp_glo,nbasmax),stat=ier) |
---|
2416 | IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for hydrodiag_glo','','') |
---|
2417 | hydrodiag=>hydrodiag_loc |
---|
2418 | |
---|
2419 | ALLOCATE(hydroupbasin_loc(nbpt),hydroupbasin_glo(nbp_glo), stat=ier) |
---|
2420 | IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for hydroupbasin_glo','','') |
---|
2421 | hydroupbasin=>hydroupbasin_loc |
---|
2422 | |
---|
2423 | END SUBROUTINE routing_simple_init_1 |
---|
2424 | |
---|
2425 | |
---|
2426 | !! ================================================================================================================================ |
---|
2427 | !! SUBROUTINE : routing_simple_clear |
---|
2428 | !! |
---|
2429 | !>\BRIEF This subroutine deallocates the block memory previously allocated. |
---|
2430 | !! |
---|
2431 | !! DESCRIPTION: This subroutine deallocates the block memory previously allocated. |
---|
2432 | !! |
---|
2433 | !! RECENT CHANGE(S): None |
---|
2434 | !! |
---|
2435 | !! MAIN OUTPUT VARIABLE(S): |
---|
2436 | !! |
---|
2437 | !! REFERENCES : None |
---|
2438 | !! |
---|
2439 | !! FLOWCHART :None |
---|
2440 | !! \n |
---|
2441 | !_ ================================================================================================================================ |
---|
2442 | |
---|
2443 | SUBROUTINE routing_simple_clear() |
---|
2444 | |
---|
2445 | IF (is_omp_root) THEN |
---|
2446 | IF (ALLOCATED(topoind_r)) DEALLOCATE(topoind_r) |
---|
2447 | IF (ALLOCATED(route_flow_rp1)) DEALLOCATE(route_flow_rp1) |
---|
2448 | IF (ALLOCATED(fast_reservoir_r)) DEALLOCATE(fast_reservoir_r) |
---|
2449 | IF (ALLOCATED(slow_reservoir_r)) DEALLOCATE(slow_reservoir_r) |
---|
2450 | IF (ALLOCATED(is_lakeinflow_r)) DEALLOCATE(is_lakeinflow_r) |
---|
2451 | IF (ALLOCATED(is_coastalflow_r)) DEALLOCATE(is_coastalflow_r) |
---|
2452 | IF (ALLOCATED(is_riverflow_r)) DEALLOCATE(is_riverflow_r) |
---|
2453 | ENDIF |
---|
2454 | |
---|
2455 | |
---|
2456 | IF (ALLOCATED(routing_area_loc)) DEALLOCATE(routing_area_loc) |
---|
2457 | IF (ALLOCATED(route_togrid_loc)) DEALLOCATE(route_togrid_loc) |
---|
2458 | IF (ALLOCATED(route_tobasin_loc)) DEALLOCATE(route_tobasin_loc) |
---|
2459 | IF (ALLOCATED(route_nbintobas_loc)) DEALLOCATE(route_nbintobas_loc) |
---|
2460 | IF (ALLOCATED(global_basinid_loc)) DEALLOCATE(global_basinid_loc) |
---|
2461 | IF (ALLOCATED(topo_resid_loc)) DEALLOCATE(topo_resid_loc) |
---|
2462 | IF (ALLOCATED(routing_area_glo)) DEALLOCATE(routing_area_glo) |
---|
2463 | IF (ALLOCATED(route_togrid_glo)) DEALLOCATE(route_togrid_glo) |
---|
2464 | IF (ALLOCATED(route_tobasin_glo)) DEALLOCATE(route_tobasin_glo) |
---|
2465 | IF (ALLOCATED(route_nbintobas_glo)) DEALLOCATE(route_nbintobas_glo) |
---|
2466 | IF (ALLOCATED(global_basinid_glo)) DEALLOCATE(global_basinid_glo) |
---|
2467 | IF (ALLOCATED(topo_resid_glo)) DEALLOCATE(topo_resid_glo) |
---|
2468 | IF (ALLOCATED(fast_reservoir)) DEALLOCATE(fast_reservoir) |
---|
2469 | IF (ALLOCATED(slow_reservoir)) DEALLOCATE(slow_reservoir) |
---|
2470 | IF (ALLOCATED(stream_reservoir)) DEALLOCATE(stream_reservoir) |
---|
2471 | IF (ALLOCATED(flood_reservoir)) DEALLOCATE(flood_reservoir) |
---|
2472 | IF (ALLOCATED(flood_frac_bas)) DEALLOCATE(flood_frac_bas) |
---|
2473 | IF (ALLOCATED(flood_height)) DEALLOCATE(flood_height) |
---|
2474 | IF (ALLOCATED(pond_frac)) DEALLOCATE(pond_frac) |
---|
2475 | IF (ALLOCATED(lake_reservoir)) DEALLOCATE(lake_reservoir) |
---|
2476 | IF (ALLOCATED(pond_reservoir)) DEALLOCATE(pond_reservoir) |
---|
2477 | IF (ALLOCATED(returnflow_mean)) DEALLOCATE(returnflow_mean) |
---|
2478 | IF (ALLOCATED(reinfiltration_mean)) DEALLOCATE(reinfiltration_mean) |
---|
2479 | IF (ALLOCATED(riverflow_mean)) DEALLOCATE(riverflow_mean) |
---|
2480 | IF (ALLOCATED(coastalflow_mean)) DEALLOCATE(coastalflow_mean) |
---|
2481 | IF (ALLOCATED(lakeinflow_mean)) DEALLOCATE(lakeinflow_mean) |
---|
2482 | IF (ALLOCATED(runoff_mean)) DEALLOCATE(runoff_mean) |
---|
2483 | IF (ALLOCATED(floodout_mean)) DEALLOCATE(floodout_mean) |
---|
2484 | IF (ALLOCATED(drainage_mean)) DEALLOCATE(drainage_mean) |
---|
2485 | IF (ALLOCATED(transpot_mean)) DEALLOCATE(transpot_mean) |
---|
2486 | IF (ALLOCATED(precip_mean)) DEALLOCATE(precip_mean) |
---|
2487 | IF (ALLOCATED(humrel_mean)) DEALLOCATE(humrel_mean) |
---|
2488 | IF (ALLOCATED(k_litt_mean)) DEALLOCATE(k_litt_mean) |
---|
2489 | IF (ALLOCATED(totnobio_mean)) DEALLOCATE(totnobio_mean) |
---|
2490 | IF (ALLOCATED(vegtot_mean)) DEALLOCATE(vegtot_mean) |
---|
2491 | IF (ALLOCATED(floodtemp)) DEALLOCATE(floodtemp) |
---|
2492 | IF (ALLOCATED(hydrodiag_loc)) DEALLOCATE(hydrodiag_loc) |
---|
2493 | IF (ALLOCATED(hydrodiag_glo)) DEALLOCATE(hydrodiag_glo) |
---|
2494 | IF (ALLOCATED(hydroupbasin_loc)) DEALLOCATE(hydroupbasin_loc) |
---|
2495 | IF (ALLOCATED(hydroupbasin_glo)) DEALLOCATE(hydroupbasin_glo) |
---|
2496 | IF (ALLOCATED(hydrographs)) DEALLOCATE(hydrographs) |
---|
2497 | IF (ALLOCATED(slowflow_diag)) DEALLOCATE(slowflow_diag) |
---|
2498 | IF (ALLOCATED(irrigation_mean)) DEALLOCATE(irrigation_mean) |
---|
2499 | IF (ALLOCATED(irrigated)) DEALLOCATE(irrigated) |
---|
2500 | IF (ALLOCATED(floodplains)) DEALLOCATE(floodplains) |
---|
2501 | IF (ALLOCATED(swamp)) DEALLOCATE(swamp) |
---|
2502 | IF (ALLOCATED(fast_diag)) DEALLOCATE(fast_diag) |
---|
2503 | IF (ALLOCATED(slow_diag)) DEALLOCATE(slow_diag) |
---|
2504 | IF (ALLOCATED(stream_diag)) DEALLOCATE(stream_diag) |
---|
2505 | IF (ALLOCATED(flood_diag)) DEALLOCATE(flood_diag) |
---|
2506 | IF (ALLOCATED(pond_diag)) DEALLOCATE(pond_diag) |
---|
2507 | IF (ALLOCATED(lake_diag)) DEALLOCATE(lake_diag) |
---|
2508 | IF (ALLOCATED(delsurfstor)) DEALLOCATE(delsurfstor) |
---|
2509 | ! |
---|
2510 | END SUBROUTINE routing_simple_clear |
---|
2511 | ! |
---|
2512 | |
---|
2513 | !! ================================================================================================================================ |
---|
2514 | !! SUBROUTINE : routing_simple_lake |
---|
2515 | !! |
---|
2516 | !>\BRIEF : This subroutine stores water in lakes so that it does not cycle through the runoff. |
---|
2517 | !! For the moment it only works for endoheric lakes but I can be extended in the future. |
---|
2518 | !! |
---|
2519 | !! DESCRIPTION (definitions, functional, design, flags): The return flow to the soil moisture reservoir |
---|
2520 | !! is based on a maximum lake evaporation rate (maxevap_lake). \n |
---|
2521 | !! |
---|
2522 | !! RECENT CHANGE(S): None |
---|
2523 | !! |
---|
2524 | !! MAIN OUTPUT VARIABLE(S): |
---|
2525 | !! |
---|
2526 | !! REFERENCES : None |
---|
2527 | !! |
---|
2528 | !! FLOWCHART :None |
---|
2529 | !! \n |
---|
2530 | !_ ================================================================================================================================ |
---|
2531 | |
---|
2532 | SUBROUTINE routing_simple_lake(nbpt, dt_routing, lakeinflow, humrel, contfrac, return_lakes) |
---|
2533 | |
---|
2534 | USE grid, ONLY : area |
---|
2535 | |
---|
2536 | IMPLICIT NONE |
---|
2537 | !! 0 Variable and parameter description |
---|
2538 | !! 0.1 Input variables |
---|
2539 | |
---|
2540 | INTEGER(i_std), INTENT(in) :: nbpt !! Domain size (unitless) |
---|
2541 | REAL(r_std), INTENT (in) :: dt_routing !! Routing time step (s) |
---|
2542 | REAL(r_std), INTENT(in) :: lakeinflow(nbpt) !! Water inflow to the lakes (kg/dt) |
---|
2543 | REAL(r_std), INTENT(in) :: humrel(nbpt) !! Soil moisture stress, root extraction potential (unitless) |
---|
2544 | REAL(r_std), INTENT(in) :: contfrac(nbpt) !! Fraction of land in each grid box (unitless;0-1) |
---|
2545 | |
---|
2546 | !! 0.2 Output variables |
---|
2547 | REAL(r_std), INTENT(out) :: return_lakes(nbpt) !! Water from lakes flowing back into soil moisture (kg/m^2/dt) |
---|
2548 | |
---|
2549 | !! 0.3 Local variables |
---|
2550 | INTEGER(i_std) :: ig !! Indices (unitless) |
---|
2551 | REAL(r_std) :: refill !! |
---|
2552 | REAL(r_std) :: total_area !! Sum of all the surfaces of the basins (m^2) |
---|
2553 | |
---|
2554 | !_ ================================================================================================================================ |
---|
2555 | |
---|
2556 | |
---|
2557 | DO ig=1,nbpt |
---|
2558 | ! |
---|
2559 | total_area = area(ig)*contfrac(ig) |
---|
2560 | ! |
---|
2561 | lake_reservoir(ig) = lake_reservoir(ig) + lakeinflow(ig) |
---|
2562 | !uptake in Kg/dt |
---|
2563 | refill = MAX(zero, maxevap_lake * (un - humrel(ig)) * dt_routing * total_area) |
---|
2564 | return_lakes(ig) = MIN(refill, lake_reservoir(ig)) |
---|
2565 | lake_reservoir(ig) = lake_reservoir(ig) - return_lakes(ig) |
---|
2566 | !Return in Kg/m^2/dt |
---|
2567 | IF ( doswamps ) THEN |
---|
2568 | return_lakes(ig) = return_lakes(ig)/total_area |
---|
2569 | ELSE |
---|
2570 | return_lakes(ig) = zero |
---|
2571 | ENDIF |
---|
2572 | ! |
---|
2573 | ! This is the volume of the lake scaled to the entire grid. |
---|
2574 | ! It would be batter to scale it to the size of the lake |
---|
2575 | ! but this information is not yet available. |
---|
2576 | lake_diag(ig) = lake_reservoir(ig)/total_area |
---|
2577 | ! |
---|
2578 | ENDDO |
---|
2579 | |
---|
2580 | END SUBROUTINE routing_simple_lake |
---|
2581 | |
---|
2582 | END MODULE routing_simple |
---|