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