! ================================================================================================================================= ! MODULE : routing_simple ! ! CONTACT : orchidee-help _at_ listes.ipsl.fr ! ! LICENCE : IPSL (2006) ! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC ! !>\BRIEF This module routes the water over the continents into the oceans and computes the water !! stored in floodplains or taken for irrigation. !! !!\n DESCRIPTION: The subroutines in this subroutine is only called when ROUTING_METHOD=simple is set in run.def. !! The method can be used for regular latitude-longitude grid or for unstructured grid. !! !! RECENT CHANGE(S): None !! !! REFERENCE(S) : !! !! SVN : !! $HeadURL$ !! $Date$ !! $Revision$ !! \n !_ ================================================================================================================================ ! Histoire Salee !--------------- ! La douce riviere ! Sortant de son lit ! S'est jetee ma chere ! dans les bras mais oui ! du beau fleuve ! ! L'eau coule sous les ponts ! Et puis les flots s'emeuvent ! - N'etes vous pas au courant ? ! Il parait que la riviere ! Va devenir mer ! Roland Bacri ! MODULE routing_simple USE ioipsl USE xios_orchidee USE ioipsl_para USE constantes USE constantes_soil USE pft_parameters USE sechiba_io_p USE interpol_help USE grid USE mod_orchidee_para IMPLICIT NONE PRIVATE PUBLIC :: routing_simple_main, routing_simple_xios_initialize PUBLIC :: routing_simple_initialize, routing_simple_finalize, routing_simple_clear !! PARAMETERS INTEGER(i_std), PARAMETER :: nbasmax=5 !! The maximum number of basins we wish to have per grid box (truncation of the model) (unitless) INTEGER(i_std), SAVE :: nbvmax !! The maximum number of basins we can handle at any time during the generation of the maps (unitless) !$OMP THREADPRIVATE(nbvmax) REAL(r_std), SAVE :: fast_tcst = 3.0 !! Property of the fast reservoir (day/m) !$OMP THREADPRIVATE(fast_tcst) REAL(r_std), SAVE :: slow_tcst = 25.0 !! Property of the slow reservoir (day/m) !$OMP THREADPRIVATE(slow_tcst) REAL(r_std), SAVE :: stream_tcst = 0.24 !! Property of the stream reservoir (day/m) !$OMP THREADPRIVATE(stream_tcst) REAL(r_std), SAVE :: flood_tcst = 4.0 !! Property of the floodplains reservoir (day/m) !$OMP THREADPRIVATE(flood_tcst) REAL(r_std), SAVE :: swamp_cst = 0.2 !! Fraction of the river transport that flows to the swamps (unitless;0-1) !$OMP THREADPRIVATE(swamp_cst) ! ! Relation between volume and fraction of floodplains ! REAL(r_std), SAVE :: beta = 2.0 !! Parameter to fix the shape of the floodplain (>1 for convex edges, <1 for concave edges) (unitless) !$OMP THREADPRIVATE(beta) REAL(r_std), SAVE :: betap = 0.5 !! Ratio of the basin surface intercepted by ponds and the maximum surface of ponds (unitless;0-1) !$OMP THREADPRIVATE(betap) REAL(r_std), SAVE :: floodcri = 2000.0 !! Potential height for which all the basin is flooded (mm) !$OMP THREADPRIVATE(floodcri) ! ! Relation between maximum surface of ponds and basin surface, and drainage (mm/j) to the slow_res ! REAL(r_std), PARAMETER :: pond_bas = 50.0 !! [DISPENSABLE] - not used REAL(r_std), SAVE :: pondcri = 2000.0 !! Potential height for which all the basin is a pond (mm) !$OMP THREADPRIVATE(pondcri) REAL(r_std), PARAMETER :: maxevap_lake = 7.5/86400. !! Maximum evaporation rate from lakes (kg/m^2/s) REAL(r_std),SAVE :: dt_routing !! Routing time step (s) !$OMP THREADPRIVATE(dt_routing) INTEGER(i_std), SAVE :: diagunit = 87 !! Diagnostic file unit (unitless) !$OMP THREADPRIVATE(diagunit) ! Logicals to control model configuration ! LOGICAL, SAVE :: dofloodinfilt = .FALSE. !! Logical to choose if floodplains infiltration is activated or not (true/false) !$OMP THREADPRIVATE(dofloodinfilt) LOGICAL, SAVE :: doswamps = .FALSE. !! Logical to choose if swamps are activated or not (true/false) !$OMP THREADPRIVATE(doswamps) LOGICAL, SAVE :: doponds = .FALSE. !! Logical to choose if ponds are activated or not (true/false) !$OMP THREADPRIVATE(doponds) INTEGER(i_std), SAVE :: num_largest !! Number of largest river basins which should be treated as independently as rivers !! (not flow into ocean as diffusion coastal flow) (unitless) !$OMP THREADPRIVATE(num_largest) REAL(r_std), SAVE :: time_counter !! Time counter (s) !$OMP THREADPRIVATE(time_counter) REAL(r_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: routing_area_loc !! Surface of basin (m^2) !$OMP THREADPRIVATE(routing_area_loc) REAL(r_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: topo_resid_loc !! Topographic index of the retention time (m) !$OMP THREADPRIVATE(topo_resid_loc) INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: route_togrid_loc !! Grid into which the basin flows (unitless) !$OMP THREADPRIVATE(route_togrid_loc) INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: route_tobasin_loc !! Basin in to which the water goes (unitless) !$OMP THREADPRIVATE(route_tobasin_loc) INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: route_nbintobas_loc !! Number of basin into current one (unitless) !$OMP THREADPRIVATE(route_nbintobas_loc) INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: global_basinid_loc !! ID of basin (unitless) !$OMP THREADPRIVATE(global_basinid_loc) INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: hydrodiag_loc !! Variable to diagnose the hydrographs !$OMP THREADPRIVATE(hydrodiag_loc) REAL(r_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:) :: hydroupbasin_loc !! The area upstream of the gauging station (m^2) !$OMP THREADPRIVATE(hydroupbasin_loc) ! parallelism REAL(r_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: routing_area_glo !! Surface of basin (m^2) !$OMP THREADPRIVATE(routing_area_glo) REAL(r_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: topo_resid_glo !! Topographic index of the retention time (m) !$OMP THREADPRIVATE(topo_resid_glo) INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: route_togrid_glo !! Grid into which the basin flows (unitless) !$OMP THREADPRIVATE(route_togrid_glo) INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: route_tobasin_glo !! Basin in to which the water goes (unitless) !$OMP THREADPRIVATE(route_tobasin_glo) INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: route_nbintobas_glo !! Number of basin into current one (unitless) !$OMP THREADPRIVATE(route_nbintobas_glo) INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: global_basinid_glo !! ID of basin (unitless) !$OMP THREADPRIVATE(global_basinid_glo) INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: hydrodiag_glo !! Variable to diagnose the hydrographs !$OMP THREADPRIVATE(hydrodiag_glo) REAL(r_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:) :: hydroupbasin_glo !! The area upstream of the gauging station (m^2) !$OMP THREADPRIVATE(hydroupbasin_glo) REAL(r_std), SAVE, POINTER, DIMENSION(:,:) :: routing_area !! Surface of basin (m^2) !$OMP THREADPRIVATE(routing_area) REAL(r_std), SAVE, POINTER, DIMENSION(:,:) :: topo_resid !! Topographic index of the retention time (m) !$OMP THREADPRIVATE(topo_resid) INTEGER(i_std), SAVE, POINTER, DIMENSION(:,:) :: route_togrid !! Grid into which the basin flows (unitless) !$OMP THREADPRIVATE(route_togrid) INTEGER(i_std), SAVE, POINTER, DIMENSION(:,:) :: route_tobasin !! Basin in to which the water goes (unitless) !$OMP THREADPRIVATE(route_tobasin) INTEGER(i_std), SAVE, POINTER, DIMENSION(:,:) :: route_nbintobas !! Number of basin into current one (unitless) !$OMP THREADPRIVATE(route_nbintobas) INTEGER(i_std), SAVE, POINTER, DIMENSION(:,:) :: global_basinid !! ID of basin (unitless) !$OMP THREADPRIVATE(global_basinid) INTEGER(i_std), SAVE, POINTER, DIMENSION(:,:) :: hydrodiag !! Variable to diagnose the hydrographs !$OMP THREADPRIVATE(hydrodiag) REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: slowflow_diag !! Diagnostic slow flow hydrographs (kg/dt) !$OMP THREADPRIVATE(slowflow_diag) REAL(r_std), SAVE, POINTER, DIMENSION(:) :: hydroupbasin !! The area upstream of the gauging station (m^2) !$OMP THREADPRIVATE(hydroupbasin) REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: irrigated !! Area equipped for irrigation in each grid box (m^2) !$OMP THREADPRIVATE(irrigated) REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: floodplains !! Maximal surface which can be inundated in each grid box (m^2) !$OMP THREADPRIVATE(floodplains) REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: swamp !! Maximal surface of swamps in each grid box (m^2) !$OMP THREADPRIVATE(swamp) ! ! The reservoirs, variables kept in the restart file. ! REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: fast_reservoir !! Water amount in the fast reservoir (kg) !$OMP THREADPRIVATE(fast_reservoir) REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: slow_reservoir !! Water amount in the slow reservoir (kg) !$OMP THREADPRIVATE(slow_reservoir) REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: stream_reservoir !! Water amount in the stream reservoir (kg) !$OMP THREADPRIVATE(stream_reservoir) REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: flood_reservoir !! Water amount in the floodplains reservoir (kg) !$OMP THREADPRIVATE(flood_reservoir) REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: lake_reservoir !! Water amount in the lake reservoir (kg) !$OMP THREADPRIVATE(lake_reservoir) REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: pond_reservoir !! Water amount in the pond reservoir (kg) !$OMP THREADPRIVATE(pond_reservoir) REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: flood_frac_bas !! Flooded fraction per basin (unitless;0-1) !$OMP THREADPRIVATE(flood_frac_bas) REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: pond_frac !! Pond fraction per grid box (unitless;0-1) !$OMP THREADPRIVATE(pond_frac) REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: flood_height !! Floodplain height (mm) !$OMP THREADPRIVATE(flood_height) ! ! The accumulated fluxes. ! REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: floodout_mean !! Accumulated flow out of floodplains (kg/m^2/dt) !$OMP THREADPRIVATE(floodout_mean) REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: runoff_mean !! Accumulated runoff (kg/m^2/dt) !$OMP THREADPRIVATE(runoff_mean) REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: drainage_mean !! Accumulated drainage (kg/m^2/dt) !$OMP THREADPRIVATE(drainage_mean) REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: transpot_mean !! Mean potential transpiration from the plants (kg/m^2/dt) !$OMP THREADPRIVATE(transpot_mean) REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: precip_mean !! Accumulated precipitation (kg/m^2/dt) !$OMP THREADPRIVATE(precip_mean) REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: humrel_mean !! Mean soil moisture stress, mean root extraction potential (unitless) !$OMP THREADPRIVATE(humrel_mean) REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: totnobio_mean !! Mean last total fraction of no bio (unitless;0-1) !$OMP THREADPRIVATE(totnobio_mean) REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: vegtot_mean !! Mean potentially vegetated fraction (unitless;0-1) !$OMP THREADPRIVATE(vegtot_mean) REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: k_litt_mean !! Mean averaged conductivity for saturated infiltration in the 'litter' layer (kg/m^2/dt) !$OMP THREADPRIVATE(k_litt_mean) ! ! The averaged outflow fluxes. ! REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: lakeinflow_mean !! Mean lake inflow (kg/m^2/dt) !$OMP THREADPRIVATE(lakeinflow_mean) REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: returnflow_mean !! Mean water flow from lakes and swamps which returns to the grid box. !! This water will go back into the hydrol or hydrolc module to allow re-evaporation (kg/m^2/dt) !$OMP THREADPRIVATE(returnflow_mean) REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: reinfiltration_mean !! Mean water flow which returns to the grid box (kg/m^2/dt) !$OMP THREADPRIVATE(reinfiltration_mean) REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: irrigation_mean !! Mean irrigation flux. !! This is the water taken from the reservoirs and beeing put into the upper layers of the soil (kg/m^2/dt) !$OMP THREADPRIVATE(irrigation_mean) REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: riverflow_mean !! Mean Outflow of the major rivers. !! The flux will be located on the continental grid but this should be a coastal point (kg/dt) !$OMP THREADPRIVATE(riverflow_mean) REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: coastalflow_mean !! Mean outflow on coastal points by small basins. !! This is the water which flows in a disperse way into the ocean (kg/dt) !$OMP THREADPRIVATE(coastalflow_mean) REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: floodtemp !! Temperature to decide if floodplains work (K) !$OMP THREADPRIVATE(floodtemp) INTEGER(i_std), SAVE :: floodtemp_lev !! Temperature level to decide if floodplains work (K) !$OMP THREADPRIVATE(floodtemp_lev) ! ! Diagnostic variables ... well sort of ! ! REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: irrig_netereq !! Irrigation requirement (water requirements by the crop for its optimal growth (kg/m^2/dt) !$OMP THREADPRIVATE(irrig_netereq) REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: hydrographs !! Hydrographs at the outflow of the grid box for major basins (kg/dt) !$OMP THREADPRIVATE(hydrographs) ! ! Diagnostics for the various reservoirs we use (Kg/m^2) ! REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: fast_diag !! Diagnostic for the fast reservoir (kg/m^2) !$OMP THREADPRIVATE(fast_diag) REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: slow_diag !! Diagnostic for the slow reservoir (kg/m^2) !$OMP THREADPRIVATE(slow_diag) REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: stream_diag !! Diagnostic for the stream reservoir (kg/m^2) !$OMP THREADPRIVATE(stream_diag) REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: flood_diag !! Diagnostic for the floodplain reservoir (kg/m^2) !$OMP THREADPRIVATE(flood_diag) REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: pond_diag !! Diagnostic for the pond reservoir (kg/m^2) !$OMP THREADPRIVATE(pond_diag) REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: lake_diag !! Diagnostic for the lake reservoir (kg/m^2) !$OMP THREADPRIVATE(lake_diag) REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: delsurfstor !! Diagnostic of the change in surface water storage (flood, pond and lake reservoirs) (kg/m^2) !$OMP THREADPRIVATE(delsurfstor) ! ! Specific variables for simple routing ! REAL(r_std),SAVE,ALLOCATABLE :: topoind_r(:) !! Topographic index of the retention time (m) index - (local routing grid) !$OMP THREADPRIVATE(topoind_r) INTEGER,SAVE,ALLOCATABLE :: route_flow_rp1(:) !! flow index from cell to neighboring cell following the trip direction - (local routing grid + halo) !$OMP THREADPRIVATE(route_flow_rp1) REAL(r_std),SAVE,ALLOCATABLE :: fast_reservoir_r(:) !! Water amount in the fast reservoir (kg) - (local routing grid) !$OMP THREADPRIVATE(fast_reservoir_r) REAL(r_std),SAVE,ALLOCATABLE :: slow_reservoir_r(:) !! Water amount in the slow reservoir (kg) - (local routing grid) !$OMP THREADPRIVATE(slow_reservoir_r) REAL(r_std),SAVE,ALLOCATABLE :: stream_reservoir_r(:) !! Water amount in the stream reservoir (kg) - (local routing grid) !$OMP THREADPRIVATE(stream_reservoir_r) LOGICAL,SAVE,ALLOCATABLE :: is_lakeinflow_r(:) !! is lake inflow point - (local routing grid) !$OMP THREADPRIVATE(is_lakeinflow_r) LOGICAL,SAVE,ALLOCATABLE :: is_coastalflow_r(:) !! is coastal flow point - (local routing grid) !$OMP THREADPRIVATE(is_coastalflow_r) LOGICAL,SAVE,ALLOCATABLE :: is_riverflow_r(:) !! is river flow point - (local routing grid) !$OMP THREADPRIVATE(is_riverflow_r) LOGICAL,SAVE,ALLOCATABLE :: is_streamflow_r(:) !! is stream flow point - (local routing grid) !$OMP THREADPRIVATE(is_streamflow_r) LOGICAL,SAVE,ALLOCATABLE :: routing_mask_r(:) !! valid routing point - (local routing grid) !$OMP THREADPRIVATE(routing_mask_r) LOGICAL,SAVE,ALLOCATABLE :: coast_mask(:) !! is a coast point - (local native grid) !$OMP THREADPRIVATE(coast_mask) INTEGER,SAVE :: total_coast_points !! global number of coast point - (local native grid) !$OMP THREADPRIVATE(total_coast_points) INTEGER,SAVE :: nbpt_r !! number of point in local routing grid !$OMP THREADPRIVATE(nbpt_r) INTEGER,SAVE :: nbpt_rp1 !! number of point in local routing grid with halo of 1 !$OMP THREADPRIVATE(nbpt_rp1) 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) !$OMP THREADPRIVATE(routing_weight) 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) !$OMP THREADPRIVATE(routing_weight_r) REAL(r_std),SAVE,ALLOCATABLE :: routing_weight_coast_r(:) !! Weight to transfer lost water into intersecting cell on coast !$OMP THREADPRIVATE(routing_weight_coast_r) !! when doing interpolation from local routing grid to local native grid (river flow+coastal flow) REAL(r_std),SAVE,ALLOCATABLE :: basins_extended_r(:) !! basins riverflow id (local routing grid) !$OMP THREADPRIVATE(basins_extended_r) INTEGER(i_std) :: basins_count !! number of basins (local routing grid) !$OMP THREADPRIVATE(basins_count) INTEGER(i_std) :: basins_out !! number of basins to output for diag !$OMP THREADPRIVATE(basins_out) INTEGER(i_std) :: split_routing !! time spliting for routing !$OMP THREADPRIVATE(split_routing) REAL(r_std), SAVE :: max_lake_reservoir !! Maximum limit of water in lake_reservoir [kg/m2] !$OMP THREADPRIVATE(max_lake_reservoir) INTEGER(i_std), PARAMETER :: nb_stations=14 REAL(r_std),PARAMETER :: station_lon(nb_stations) = & (/ -162.8830, -90.9058, -55.5110, -49.3242, -133.7447, -63.6000, 28.7167, & 15.3000, 66.5300, 89.6700, 86.5000, 127.6500, 3.3833, 117.6200 /) REAL(r_std),PARAMETER :: station_lat(nb_stations) = & (/ 61.9340, 32.3150, -1.9470, -5.1281, 67.4583, 8.1500, 45.2167, & -4.3000, 66.5700, 25.1800, 67.4800, 70.7000, 11.8667, 30.7700 /) CHARACTER(LEN=17),PARAMETER :: station_name(nb_stations) = & (/ "Pilot station ", "Vicksburg ", "Obidos ", & "Itupiranga ", "Arctic red river ", "Puente Angostura ", & "Ceatal Izmail ", "Kinshasa ", "Salekhard ", & "Bahadurabad ", "Igarka ", "Kusur ", & "Malanville ", "Datong " /) CONTAINS !! ============================================================================================================================= !! SUBROUTINE: routing_simple_xios_initialize !! !>\BRIEF Initialize xios dependant defintion before closing context defintion !! !! DESCRIPTION: Initialize xios dependant defintion before closing context defintion. !! This subroutine is called before the xios context is closed. !! !! RECENT CHANGE(S): None !! !! REFERENCE(S): None !! !! FLOWCHART: None !! \n !_ ============================================================================================================================== SUBROUTINE routing_simple_xios_initialize USE xios USE routing, ONLY : routing_names IMPLICIT NONE INTEGER(i_std) ::ib !! 0 Variable and parameter description CHARACTER(LEN=60),ALLOCATABLE :: label(:) LOGICAL :: file_exists IF (is_omp_root) THEN CALL xios_get_axis_attr("basins", n_glo=basins_out) ! get nb basins to output ALLOCATE(label(basins_out)) CALL routing_names(basins_out,label) CALL xios_set_axis_attr("basins", label=label) ! set riverflow basins name INQUIRE(FILE="routing_start.nc", EXIST=file_exists) IF (file_exists) CALL xios_set_file_attr("routing_start", enabled=.TRUE.) ENDIF !! Define XIOS axis size needed for the model output ! Add axis for homogeneity between all routing schemes, these dimensions are currently not used in this scheme CALL xios_orchidee_addaxis("nbhtu", nbasmax, (/(REAL(ib,r_std),ib=1,nbasmax)/)) CALL xios_orchidee_addaxis("nbasmon", 1, (/(REAL(ib,r_std),ib=1,1)/)) END SUBROUTINE routing_simple_xios_initialize !! ============================================================================================================================= !! SUBROUTINE: routing_simple_init_2 !! !>\BRIEF Privat subroutine !! !! DESCRIPTION: Privat subroutine to the module routing_simple. This subroutine is called in the end of routing_simple_initialize !! !! RECENT CHANGE(S): None !! !! REFERENCE(S): None !! !! FLOWCHART: None !! \n !_ ============================================================================================================================== SUBROUTINE routing_simple_init_2(nbpt, contfrac) USE xios USE grid, ONLY : area IMPLICIT NONE INCLUDE "mpif.h" !! 0 Variable and parameter description !! 0.1 Input variables INTEGER,INTENT(in) :: nbpt !! nb points native grid REAL,INTENT(in) :: contfrac(nbpt) !! fraction of land !! 0.2 Local variables INTEGER :: ni !! longitude dimension of local routing grid INTEGER :: nj !! latitude dimension of local routing grid 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) 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) !! routing is artificially computed on sea and endoric basins LOGICAL,ALLOCATABLE :: routing_mask_rp1(:,:) !! valid routing point - local routing grid+halo (0:ni+1,0:nj+1) REAL(r_std),ALLOCATABLE :: mask_native(:) !! some mask on native grid (nbpt) REAL(r_std),ALLOCATABLE :: frac_routing_r(:) !! fraction of routing cell intesected by cells of native grid REAL(r_std),ALLOCATABLE :: frac_routing_coast_r(:) !! fraction of routing cell intesected by coastal cells of native grid INTEGER :: ij, ij_r, ij_rp1, i,j,jp1,jm1,ip1,im1,jr,ir REAL(r_std) :: sum1,sum2,sum_frac_routing_r REAL(r_std) :: contfrac_mpi(nbp_mpi) REAL(r_std) :: area_mpi(nbp_mpi) INTEGER :: basins_count_mpi INTEGER :: nb_coast_points INTEGER :: ierr LOGICAL :: file_exists !_ ================================================================================================================================ split_routing=1 CALL getin_p("SPLIT_ROUTING",split_routing) CALL gather_omp(contfrac,contfrac_mpi) CALL gather_omp(area, area_mpi) IF (is_omp_root) THEN CALL xios_get_domain_attr("routing_domain", ni=ni, nj=nj) ! get routing domain dimension nbpt_r= ni*nj nbpt_rp1= (ni+2)*(nj+2) ALLOCATE(trip_extended_r(ni,nj)) ALLOCATE(topoind_r(ni*nj)) ALLOCATE(is_lakeinflow_r(ni*nj)) ALLOCATE(is_coastalflow_r(ni*nj)) ALLOCATE(is_riverflow_r(ni*nj)) ALLOCATE(is_streamflow_r(ni*nj)) ALLOCATE(routing_mask_r(ni*nj)) ALLOCATE(fast_reservoir_r(ni*nj)) ALLOCATE(slow_reservoir_r(ni*nj)) ALLOCATE(stream_reservoir_r(ni*nj)) ALLOCATE(mask_native(nbp_mpi)) ALLOCATE(coast_mask(nbp_mpi)) ALLOCATE(frac_routing_r(nbpt_r)) ALLOCATE(frac_routing_coast_r(nbpt_r)) ALLOCATE(routing_weight(nbp_mpi)) ALLOCATE(routing_weight_r(nbpt_r)) ALLOCATE(routing_weight_coast_r(nbpt_r)) ALLOCATE(basins_extended_r(nbpt_r)) coast_mask=.FALSE. WHERE( contfrac_mpi(:)< 1.-1.e-5) coast_mask(:)=.TRUE. ! create mask for coastal cells on native grid INQUIRE(FILE="routing_start.nc", EXIST=file_exists) IF (file_exists) THEN CALL xios_recv_field("fast_reservoir_start",fast_reservoir_r) CALL xios_recv_field("slow_reservoir_start",slow_reservoir_r) CALL xios_recv_field("stream_reservoir_start",stream_reservoir_r) ELSE fast_reservoir_r(:)=0 slow_reservoir_r(:)=0 stream_reservoir_r(:)=0 ENDIF ALLOCATE(trip_rp1(0:ni+1,0:nj+1)) trip_rp1(:,:)=1e10 ALLOCATE(routing_mask_rp1(0:ni+1,0:nj+1)) ALLOCATE(route_flow_rp1((ni+2)*(nj+2))) CALL xios_send_field("routing_contfrac",contfrac_mpi) ! diag CALL xios_recv_field("trip_r",trip_rp1(1:ni,1:nj)) ! recv trip array with halo of 1 CALL xios_recv_field("trip_extended_r",trip_extended_r) ! recv extended trip array from file CALL xios_recv_field("topoind_r",topoind_r) ! recv topo index array from file CALL xios_recv_field("basins_extended_r",basins_extended_r) ! recv basins index from file is_lakeinflow_r(:) = .FALSE. is_coastalflow_r(:) = .FALSE. is_riverflow_r(:) = .FALSE. route_flow_rp1(:)=-1 routing_mask_rp1=.TRUE. DO j=1,nj DO i=1,ni IF (trip_rp1(i,j)>99) routing_mask_rp1(i,j)=.FALSE. ! if ocean or endoric basins => masked point ENDDO ENDDO mask_native(:)=0 frac_routing_r(:)=0 WHERE( .NOT. coast_mask(:)) mask_native(:)=1 nb_coast_points=SUM(1-mask_native) CALL reduce_sum_mpi(nb_coast_points, total_coast_points) CALL bcast_mpi(total_coast_points) CALL xios_send_field("mask_native_lake",mask_native) ! send full land point to XIOS (native grid) CALL xios_recv_field("frac_routing_lake_r",frac_routing_r) ! receive fraction of intersected cell by full land, on routing grid DO j=1,nj DO i=1,ni ij=i+ni*(j-1) IF (frac_routing_r(ij)>1-1e-5) THEN ! if full land point (no part of coast) 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 routing_mask_rp1(i,j)=.TRUE. ! transform to non masked trip_rp1(i,j)=trip_extended_r(i,j) ! replace by extended trip value ENDIF ENDIF ENDDO ENDDO mask_native(:)=0 frac_routing_coast_r(:)=0 WHERE( coast_mask(:)) mask_native(:)=1 CALL xios_send_field("mask_native_coast",mask_native) ! send only coastal point to XIOS (native grid) CALL xios_recv_field("frac_routing_coast_r",frac_routing_coast_r) ! receive fraction of intersected cell by coastal point, on routing grid DO j=1,nj DO i=1,ni ij=i+ni*(j-1) IF (frac_routing_coast_r(ij) > 0) THEN ! fraction of coastal point ? IF (.NOT. routing_mask_rp1(i,j)) THEN ! if masked, transform to coastal flow to conserve the water balance routing_mask_rp1(i,j)=.TRUE. trip_rp1(i,j)=98 ! transform to coastal flow topoind_r(ij)=1e-10 ! no residence time ENDIF ENDIF ENDDO ENDDO mask_native(:)=1 frac_routing_r(:)=0 CALL xios_send_field("mask_native",mask_native) ! send only all point to XIOS (native grid) CALL xios_recv_field("frac_routing_r",frac_routing_r) ! receive fraction of intersected cell, (routing grid) DO j=1,nj DO i=1,ni ij=i+ni*(j-1) IF (frac_routing_r(ij)>0) THEN IF (.NOT. routing_mask_rp1(i,j)) THEN routing_mask_rp1(i,j)=.TRUE. ! may never happen ENDIF ELSE routing_mask_rp1(i,j)=.FALSE. ! if no intersection from native cells, mask the point (no routage) trip_rp1(i,j)=1e10 ENDIF ENDDO ENDDO CALL xios_send_field("trip_update_r",trip_rp1(1:ni,1:nj)) ! send to xios trip array to update for halo CALL xios_recv_field("trip_rp1",trip_rp1) ! recv trip array with halo of 1 routing_mask_rp1=.TRUE. DO j=0,nj+1 DO i=0,ni+1 IF (trip_rp1(i,j)>99) routing_mask_rp1(i,j)=.FALSE. ! update mask with halo of 1 ENDDO ENDDO !! Compute the routing !! Loop on all point point of the local routing grid + halo DO j=0,nj+1 jp1=j+1 jm1=j-1 DO i=0,ni+1 ij_rp1=i+(ni+2)*j+1 ij_r=i+ni*(j-1) ip1=i+1 im1=i-1 IF (trip_rp1(i,j) < 100) THEN ir=-1 ; jr=-1 ! -> -1 for 97,98,99 SELECT CASE (NINT(trip_rp1(i,j))) ! get the trip value for each points CASE (1) jr=jm1 ; ir=i ! north CASE (2) jr=jm1 ; ir=ip1 ! north-east CASE (3) jr=j ; ir=ip1 ! east CASE (4) jr=jp1 ; ir=ip1 ! south-east CASE (5) jr=jp1 ; ir=i ! south CASE(6) jr=jp1 ; ir=im1 ! south-west CASE (7) jr=j ; ir=im1 ! west CASE (8) jr=jm1 ; ir=im1 ! north-west CASE (97) IF ( i>0 .AND. i0 .AND. j0 .AND. i0 .AND. j0 .AND. i0 .AND. jni+1 .OR. jr<0 .OR. jr>nj+1) THEN route_flow_rp1(ij_rp1)=-1 ! if route outside my local domain+halo, no routing (will be done by other process) ELSE IF ( .NOT. routing_mask_rp1(ir,jr)) THEN ! if routing to a masked point, cut the flow jr=j ; ir=i ; ! and route to myself IF (i>0 .AND. i0 .AND. j transform to riverflow ELSE is_coastalflow_r(ij_r)=.TRUE. ! else transform to coastalflow ENDIF ENDIF ENDIF IF (ir<1 .OR. ir>ni .OR. jr<1 .OR. jr>nj) THEN route_flow_rp1(ij_rp1)=-1 ! if route outside my local domain, no routing (will be done by other process) ELSE route_flow_rp1(ij_rp1)=ir+ni*(jr-1) ! define the cell where to flow ENDIF ENDIF ENDIF ENDDO ENDDO routing_mask_r(:)=reshape(routing_mask_rp1(1:ni,1:nj),(/ni*nj/)) is_streamflow_r(:)=.NOT. (is_lakeinflow_r(:) .OR. is_coastalflow_r(:) .OR. is_riverflow_r(:)) .AND. routing_mask_r(:) DO ij=1,nbp_mpi routing_weight(ij)=contfrac_mpi(ij)*area_mpi(ij) ENDDO routing_weight_r(:)=0 DO ij=1,nbpt_r IF (frac_routing_r(ij)>0) routing_weight_r(ij)=1./frac_routing_r(ij) ENDDO routing_weight_coast_r(:)=0. DO ij=1,nbpt_r IF (frac_routing_coast_r(ij)>0) routing_weight_coast_r(ij)=1./frac_routing_coast_r(ij) ENDDO CALL xios_send_field("routing_weight_coast_r",routing_weight_coast_r) ! looking for basins basins_count_mpi=0 DO ij=1,nbpt_r IF (basins_extended_r(ij) > basins_count_mpi) basins_count_mpi=basins_extended_r(ij) ENDDO CALL MPI_ALLREDUCE(basins_count_mpi,basins_count,1,MPI_INT_ORCH, MPI_MAX,MPI_COMM_ORCH,ierr) ELSE nbpt_r=0 ENDIF END SUBROUTINE routing_simple_init_2 !! ================================================================================================================================ !! SUBROUTINE : routing_simple_flow !! !>\BRIEF This subroutine computes the transport of water in the various reservoirs !! (including ponds and floodplains) and the water withdrawals from the reservoirs for irrigation. !! !! DESCRIPTION (definitions, functional, design, flags) : !! This will first compute the amount of water which flows out of each of the 3 reservoirs using the assumption of an !! exponential decrease of water in the reservoir (see Hagemann S and Dumenil L. (1998)). Then we compute the fluxes !! for floodplains and ponds. All this will then be used in order to update each of the basins : taking water out of !! the up-stream basin and adding it to the down-stream one. !! As this step happens globaly we have to stop the parallel processing in order to exchange the information. Once !! all reservoirs are updated we deal with irrigation. The final step is to compute diagnostic fluxes. Among them !! the hydrographs of the largest rivers we have chosen to monitor. !! !! RECENT CHANGE(S): None !! !! MAIN OUTPUT VARIABLE(S): lakeinflow, returnflow, reinfiltration, irrigation, riverflow, coastalflow, hydrographs, flood_frac, flood_res !! !! REFERENCES : !! - Ngo-Duc, T., K. Laval, G. Ramillien, J. Polcher, and A. Cazenave (2007) !! Validation of the land water storage simulated by Organising Carbon and Hydrology in Dynamic Ecosystems (ORCHIDEE) with Gravity Recovery and Climate Experiment (GRACE) data. !! Water Resour. Res., 43, W04427, doi:10.1029/2006WR004941. !! * Irrigation: !! - de Rosnay, P., J. Polcher, K. Laval, and M. Sabre (2003) !! Integrated parameterization of irrigation in the land surface model ORCHIDEE. Validation over Indian Peninsula. !! Geophys. Res. Lett., 30(19), 1986, doi:10.1029/2003GL018024. !! - A.C. Vivant (2003) !! Les plaines d'inondations et l'irrigation dans ORCHIDEE, impacts de leur prise en compte. !! , , 51pp. !! - N. Culson (2004) !! Impact de l'irrigation sur le cycle de l'eau !! Master thesis, Paris VI University, 55pp. !! - X.-T. Nguyen-Vinh (2005) !! Analyse de l'impact de l'irrigation en Amerique du Nord - plaine du Mississippi - sur la climatologie regionale !! Master thesis, Paris VI University, 33pp. !! - M. Guimberteau (2006) !! Analyse et modifications proposees de la modelisation de l'irrigation dans un modele de surface. !! Master thesis, Paris VI University, 46pp. !! - Guimberteau M. (2010) !! Modelisation de l'hydrologie continentale et influences de l'irrigation sur le cycle de l'eau. !! Ph.D. thesis, Paris VI University, 195pp. !! - Guimberteau M., Laval K., Perrier A. and Polcher J. (2011). !! Global effect of irrigation and its impact on the onset of the Indian summer monsoon. !! In press, Climate Dynamics, doi: 10.1007/s00382-011-1252-5. !! * Floodplains: !! - A.C. Vivant (2002) !! L'ecoulement lateral de l'eau sur les surfaces continentales. Prise en compte des plaines d'inondations dans ORCHIDEE. !! Master thesis, Paris VI University, 46pp. !! - A.C. Vivant (2003) !! Les plaines d'inondations et l'irrigation dans ORCHIDEE, impacts de leur prise en compte. !! , , 51pp. !! - T. d'Orgeval (2006) !! Impact du changement climatique sur le cycle de l'eau en Afrique de l'Ouest: modelisation et incertitudes. !! Ph.D. thesis, Paris VI University, 188pp. !! - T. d'Orgeval, J. Polcher, and P. de Rosnay (2008) !! Sensitivity of the West African hydrological cycle in ORCHIDEE to infiltration processes. !! Hydrol. Earth Syst. Sci., 12, 1387-1401 !! - M. Guimberteau, G. Drapeau, J. Ronchail, B. Sultan, J. Polcher, J.-M. Martinez, C. Prigent, J.-L. Guyot, G. Cochonneau, !! J. C. Espinoza, N. Filizola, P. Fraizy, W. Lavado, E. De Oliveira, R. Pombosa, L. Noriega, and P. Vauchel (2011) !! Discharge simulation in the sub-basins of the Amazon using ORCHIDEE forced by new datasets. !! Hydrol. Earth Syst. Sci. Discuss., 8, 11171-11232, doi:10.5194/hessd-8-11171-2011 !! !! FLOWCHART :None !! \n !_ ================================================================================================================================ SUBROUTINE routing_simple_flow(nbpt, dt_routing, lalo, floodout, runoff_omp, drainage_omp, & vegtot, totnobio, transpot_mean, precip, humrel, k_litt, floodtemp, reinf_slope, & lakeinflow_omp, returnflow, reinfiltration, irrigation, riverflow_omp, & coastalflow_omp, hydrographs, slowflow_diag, flood_frac, flood_res) USE xios USE grid, ONLY : area IMPLICIT NONE INCLUDE "mpif.h" !! 0 Variable and parameter description !! 0.1 Input variables INTEGER(i_std), INTENT(in) :: nbpt !! Domain size (unitless) REAL(r_std), INTENT (in) :: dt_routing !! Routing time step (s) REAL(r_std), INTENT(in) :: lalo(nbpt,2) !! Vector of latitude and longitudes REAL(r_std), INTENT(in) :: runoff_omp(nbpt) !! Grid-point runoff (kg/m^2/dt) REAL(r_std), INTENT(in) :: floodout(nbpt) !! Grid-point flow out of floodplains (kg/m^2/dt) REAL(r_std), INTENT(in) :: drainage_omp(nbpt) !! Grid-point drainage (kg/m^2/dt) REAL(r_std), INTENT(in) :: vegtot(nbpt) !! Potentially vegetated fraction (unitless;0-1) REAL(r_std), INTENT(in) :: totnobio(nbpt) !! Other areas which can not have vegetation REAL(r_std), INTENT(in) :: transpot_mean(nbpt) !! Mean potential transpiration of the vegetation (kg/m^2/dt) REAL(r_std), INTENT(in) :: precip(nbpt) !! Rainfall (kg/m^2/dt) REAL(r_std), INTENT(in) :: humrel(nbpt) !! Soil moisture stress, root extraction potential (unitless) REAL(r_std), INTENT(in) :: k_litt(nbpt) !! Averaged conductivity for saturated infiltration in the 'litter' layer (kg/m^2/dt) REAL(r_std), INTENT(in) :: floodtemp(nbpt) !! Temperature to decide if floodplains work (K) 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) !! 0.2 Output variables REAL(r_std), INTENT(out) :: lakeinflow_omp(nbpt) !! Water inflow to the lakes (kg/dt) REAL(r_std), INTENT(out) :: returnflow(nbpt) !! The water flow from lakes and swamps which returns into the grid box. !! This water will go back into the hydrol or hydrolc module to allow re-evaporation (kg/m^2/dt) REAL(r_std), INTENT(out) :: reinfiltration(nbpt) !! Water flow from ponds and floodplains which returns to the grid box (kg/m^2/dt) 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) 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) 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) REAL(r_std), INTENT(out) :: hydrographs(nbpt) !! Hydrographs at the outflow of the grid box for major basins (kg/dt) REAL(r_std), INTENT(out) :: slowflow_diag(nbpt) !! Hydrographs of slow_flow = routed slow_flow for major basins (kg/dt) REAL(r_std), INTENT(out) :: flood_frac(nbpt) !! Flooded fraction of the grid box (unitless;0-1) REAL(r_std), INTENT(out) :: flood_res(nbpt) !! Diagnostic of water amount in the floodplains reservoir (kg) !! 0.4 Local variables REAL(r_std) :: runoff(nbp_mpi) !! Grid-point runoff (kg/dt) REAL(r_std) :: drainage(nbp_mpi) !! Grid-point drainage (kg/dt) REAL(r_std) :: riverflow(nbp_mpi) REAL(r_std) :: coastalflow(nbp_mpi) REAL(r_std) :: lakeinflow(nbp_mpi) REAL(r_std) :: fast_diag_mpi(nbp_mpi) REAL(r_std) :: slow_diag_mpi(nbp_mpi) REAL(r_std) :: stream_diag_mpi(nbp_mpi) REAL(r_std) :: area_mpi(nbp_mpi) ! cell area REAL(r_std) :: lake_reservoir_mpi(nbp_mpi) ! cell area ! from input model -> routing_grid REAL(r_std) :: runoff_r(nbpt_r) !! Grid-point runoff (kg/m^2/dt) REAL(r_std) :: floodout_r(nbpt_r) !! Grid-point flow out of floodplains (kg/m^2/dt) REAL(r_std) :: drainage_r(nbpt_r) !! Grid-point drainage (kg/m^2/dt) REAL(r_std) :: vegtot_r(nbpt_r) !! Potentially vegetated fraction (unitless;0-1) REAL(r_std) :: totnobio_r(nbpt_r) !! Other areas which can not have vegetation REAL(r_std) :: transpot_mean_r(nbpt_r) !! Mean potential transpiration of the vegetation (kg/m^2/dt) REAL(r_std) :: precip_r(nbpt_r) !! Rainfall (kg/m^2/dt) REAL(r_std) :: humrel_r(nbpt_r) !! Soil moisture stress, root extraction potential (unitless) REAL(r_std) :: k_litt_r(nbpt_r) !! Averaged conductivity for saturated infiltration in the 'litter' layer (kg/m^2/dt) REAL(r_std) :: floodtemp_r(nbpt_r) !! Temperature to decide if floodplains work (K) 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) REAL(r_std), DIMENSION(nbpt_r) :: fast_flow_r !! Outflow from the fast reservoir (kg/dt) REAL(r_std), DIMENSION(nbpt_r) :: slow_flow_r !! Outflow from the slow reservoir (kg/dt) REAL(r_std), DIMENSION(nbpt_r) :: stream_flow_r !! Outflow from the stream reservoir (kg/dt) REAL(r_std), DIMENSION(nbpt_r) :: hydrographs_r !! hydrograph (kg/dt) REAL(r_std), DIMENSION(nbpt_r) :: flood_flow_r !! Outflow from the floodplain reservoir (kg/dt) REAL(r_std), DIMENSION(nbpt_r) :: pond_inflow_r !! Inflow to the pond reservoir (kg/dt) REAL(r_std), DIMENSION(nbpt_r) :: pond_drainage_r !! Drainage from pond (kg/m^2/dt) REAL(r_std), DIMENSION(nbpt_r) :: flood_drainage_r !! Drainage from floodplains (kg/m^2/dt) REAL(r_std), DIMENSION(nbpt_r) :: return_swamp_r !! Inflow to the swamp (kg/dt) ! ! Irrigation per basin ! REAL(r_std), DIMENSION(nbpt_r) :: irrig_needs_r !! Total irrigation requirement (water requirements by the crop for its optimal growth) (kg) REAL(r_std), DIMENSION(nbpt_r) :: irrig_actual_r !! Possible irrigation according to the water availability in the reservoirs (kg) REAL(r_std), DIMENSION(nbpt_r) :: irrig_deficit_r !! Amount of water missing for irrigation (kg) REAL(r_std), DIMENSION(nbpt_r) :: irrig_adduct_r !! Amount of water carried over from other basins for irrigation (kg) ! REAL(r_std), DIMENSION(nbpt_r) :: transport_r !! Water transport between basins (kg/dt) REAL(r_std), DIMENSION(nbpt_r) :: floods_r !! Water flow in to the floodplains (kg/dt) REAL(r_std), DIMENSION(nbpt_r) :: potflood_r !! Potential inflow to the swamps (kg/dt) REAL(r_std), DIMENSION(nbpt) :: tobeflooded !! Maximal surface which can be inundated in each grid box (m^2) REAL(r_std), DIMENSION(nbpt) :: totarea !! Total area of basin (m^2) REAL(r_std), DIMENSION(nbpt) :: totflood !! Total amount of water in the floodplains reservoir (kg) (not used) REAL(r_std), DIMENSION(nbasmax) :: pond_excessflow !! REAL(r_std) :: flow !! Outflow computation for the reservoirs (kg/dt) REAL(r_std) :: floodindex !! Fraction of grid box area inundated (unitless;0-1) REAL(r_std) :: pondex !! REAL(r_std) :: flood_frac_pot !! Total fraction of the grid box which is flooded at optimum repartition (unitless;0-1) REAL(r_std) :: stream_tot !! Total water amount in the stream reservoirs (kg) REAL(r_std) :: adduction !! Importation of water from a stream reservoir of a neighboring grid box (kg) REAL(r_std), DIMENSION(8,nbasmax) :: streams_around !! Stream reservoirs of the neighboring grid boxes (kg) INTEGER(i_std), DIMENSION(8) :: igrd !! INTEGER(i_std), DIMENSION(2) :: ff !! INTEGER(i_std), DIMENSION(1) :: fi !! INTEGER(i_std) :: ig, ib, ib2, ig2 !! Indices (unitless) INTEGER(i_std) :: rtg, rtb, in !! Indices (unitless) INTEGER(i_std) :: ier !! Error handling INTEGER(i_std) :: isplit 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) REAL(r_std), DIMENSION(nbpt_r) :: flood_frac_bas_r REAL(r_std), DIMENSION(nbpt_r) :: flood_reservoir_r REAL(r_std), DIMENSION(nbpt_r) :: pond_reservoir_r REAL(r_std), DIMENSION(nbpt_r) :: lakeinflow_r REAL(r_std), DIMENSION(nbpt_r) :: coastalflow_r REAL(r_std), DIMENSION(nbpt_r) :: riverflow_r REAL(r_std), DIMENSION(nbpt_r) :: flow_r REAL(r_std), DIMENSION(nbpt_rp1) :: flow_rp1 REAL(r_std) :: water_balance_before, water_balance_after REAL(r_std) :: basins_riverflow_mpi(0:basins_count) REAL(r_std) :: basins_riverflow(0:basins_count) REAL(r_std) :: lake_overflow,sum_lake_overflow, total_lake_overflow INTEGER :: ierr !_ ================================================================================================================================ ! irrig_netereq(:) = zero totarea(:) = zero totflood(:) = zero flood_frac_bas_r(:)=zero ! !> The outflow fluxes from the three reservoirs are computed. !> The outflow of volume of water Vi into the reservoir i is assumed to be linearly related to its volume. !> The water travel simulated by the routing scheme is dependent on the water retention index topo_resid !> given by a 0.5 degree resolution map for each pixel performed from a simplification of Manning's formula !> (Dingman, 1994; Ducharne et al., 2003). !> The resulting product of tcst (in day/m) and topo_resid (in m) represents the time constant (day) !> which is an e-folding time, the time necessary for the water amount !> in the stream reservoir to decrease by a factor e. Hence, it gives an order of !> magnitude of the travel time through this reservoir between !> the sub-basin considered and its downstream neighbor. CALL gather_omp(runoff_omp,runoff) CALL gather_omp(drainage_omp, drainage) CALL gather_omp(area, area_mpi) CALL gather_omp(lake_reservoir, lake_reservoir_mpi) IF (is_omp_root) THEN hydrographs_r(:)=0 DO isplit=1,split_routing irrig_needs_r(:) = zero irrig_actual_r(:) = zero irrig_deficit_r(:) = zero irrig_adduct_r(:) = zero DO ig=1,nbpt_r IF ( routing_mask_r(ig) ) THEN ! ! Each of the fluxes is limited by the water in the reservoir and a small margin ! (min_reservoir) to avoid rounding errors. ! ! IF (fast_reservoir_r(ig)>1e-30) THEN ! PRINT*,"fast_reservoir > 0 ",ig, fast_reservoir_r(ig), topoind_r(ig) ! ENDIF flow = MIN(fast_reservoir_r(ig)/((topoind_r(ig)/1000.)*fast_tcst*one_day/(dt_routing/split_routing)),& & fast_reservoir_r(ig)-min_sechiba) fast_flow_r(ig) = MAX(flow, zero) ! flow = MIN(slow_reservoir_r(ig)/((topoind_r(ig)/1000.)*slow_tcst*one_day/(dt_routing/split_routing)),& & slow_reservoir_r(ig)-min_sechiba) slow_flow_r(ig) = MAX(flow, zero) ! flow = MIN(stream_reservoir_r(ig)/((topoind_r(ig)/1000.)*stream_tcst* & & MAX(un-SQRT(flood_frac_bas_r(ig)),min_sechiba)*one_day/(dt_routing/split_routing)),& & stream_reservoir_r(ig)-min_sechiba) stream_flow_r(ig) = MAX(flow, zero) ! ! ELSE fast_flow_r(ig) = zero slow_flow_r(ig) = zero stream_flow_r(ig) = zero ENDIF ENDDO DO ig=1,nbpt_r flood_drainage_r(ig) = zero flood_flow_r(ig) = zero flood_reservoir_r(ig) = zero ENDDO DO ig=1,nbpt_r pond_inflow_r(ig) = zero pond_drainage_r(ig) = zero pond_reservoir_r(ig) = zero ENDDO !- !- Compute the transport !- flow_r(:)=fast_flow_r(:) + slow_flow_r(:) + stream_flow_r(:) CALL xios_send_field("flow_r",flow_r) ! transfer halo CALL xios_recv_field("flow_rp1",flow_rp1) transport_r(:)=0 DO ig=1,nbpt_rp1 IF ( route_flow_rp1(ig) > 0 ) THEN transport_r(route_flow_rp1(ig))=transport_r(route_flow_rp1(ig))+ flow_rp1(ig) ENDIF ENDDO !- !- Do the floodings - First initialize !- return_swamp_r(:)=zero floods_r(:)=zero ! ! Update all reservoirs !> The slow and deep reservoir (slow_reservoir) collect the deep drainage whereas the !> fast_reservoir collects the computed surface runoff. Both discharge into a third reservoir !> (stream_reservoir) of the next sub-basin downstream. !> Water from the floodplains reservoir (flood_reservoir) flows also into the stream_reservoir of the next sub-basin downstream. !> Water that flows into the pond_reservoir is withdrawn from the fast_reservoir. ! runoff=runoff*routing_weight CALL xios_send_field("routing_runoff",runoff) ! interp conservative model -> routing CALL xios_recv_field("routing_runoff_r",runoff_r) drainage=drainage*routing_weight CALL xios_send_field("routing_drainage",drainage) ! interp conservative model -> routing CALL xios_recv_field("routing_drainage_r",drainage_r) CALL xios_send_field("fast_reservoir_r",fast_reservoir_r) CALL xios_send_field("slow_reservoir_r",slow_reservoir_r) CALL xios_send_field("stream_reservoir_r",stream_reservoir_r) WHERE (.NOT. routing_mask_r) runoff_r(:)=0 WHERE (.NOT. routing_mask_r) drainage_r(:)=0 CALL MPI_ALLREDUCE(sum(runoff+drainage),water_balance_before,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr) CALL MPI_ALLREDUCE(sum(runoff_r+drainage_r),water_balance_after,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr) PRINT *,"routing water Balance ; before : ", water_balance_before," ; after : ",water_balance_after, & " ; delta : ", 100.*(water_balance_after-water_balance_before)/(0.5*(water_balance_after+water_balance_before)),"%" CALL xios_recv_field("water_balance_before",water_balance_before) DO ig=1,nbpt_r IF ( routing_mask_r(ig) ) THEN ! fast_reservoir_r(ig) = fast_reservoir_r(ig) + runoff_r(ig) - & & fast_flow_r(ig) - pond_inflow_r(ig) ! slow_reservoir_r(ig) = slow_reservoir_r(ig) + drainage_r(ig) - & & slow_flow_r(ig) ! stream_reservoir_r(ig) = stream_reservoir_r(ig) + flood_flow_r(ig) - & & stream_flow_r(ig) - return_swamp_r(ig) - floods_r(ig) IF (is_streamflow_r(ig)) stream_reservoir_r(ig)= stream_reservoir_r(ig) + transport_r(ig) ! flood_reservoir_r(ig) = flood_reservoir_r(ig) + floods_r(ig) - & & flood_flow_r(ig) ! pond_reservoir_r(ig) = pond_reservoir_r(ig) + pond_inflow_r(ig) - pond_drainage_r(ig) ! IF ( flood_reservoir_r(ig) .LT. zero ) THEN IF ( check_reservoir ) THEN WRITE(numout,*) "WARNING : negative flood reservoir at :", ig, ". Problem is being corrected." WRITE(numout,*) "flood_reservoir, floods, flood_flow : ", flood_reservoir_r(ig), floods_r(ig), & & flood_flow_r(ig) ENDIF stream_reservoir_r(ig) = stream_reservoir_r(ib) + flood_reservoir_r(ib) flood_reservoir_r(ib) = zero ENDIF ! IF ( stream_reservoir_r(ig) .LT. zero ) THEN IF ( check_reservoir ) THEN WRITE(numout,*) "WARNING : negative stream reservoir at :", ig, ". Problem is being corrected." WRITE(numout,*) "stream_reservoir, flood_flow, transport : ", stream_reservoir_r(ig), flood_flow_r(ig), & & transport_r(ig) WRITE(numout,*) "stream_flow, return_swamp, floods :", stream_flow_r(ig), return_swamp_r(ig), floods_r(ig) ENDIF fast_reservoir_r(ig) = fast_reservoir_r(ig) + stream_reservoir_r(ig) stream_reservoir_r(ig) = zero ENDIF ! IF ( fast_reservoir_r(ig) .LT. zero ) THEN IF ( check_reservoir ) THEN WRITE(numout,*) "WARNING : negative fast reservoir at :", ig, ". Problem is being corrected." WRITE(numout,*) "fast_reservoir, runoff, fast_flow, ponf_inflow : ", fast_reservoir_r(ig), & &runoff_r(ig), fast_flow_r(ig), pond_inflow_r(ig) ENDIF slow_reservoir_r(ig) = slow_reservoir_r(ig) + fast_reservoir_r(ig) fast_reservoir_r(ig) = zero ENDIF IF ( slow_reservoir_r(ig) .LT. - min_sechiba ) THEN WRITE(numout,*) 'WARNING : There is a negative reservoir at :', ig WRITE(numout,*) 'WARNING : slowr, slow_flow, drainage', & & slow_reservoir_r(ig), slow_flow_r(ig), drainage_r(ig) WRITE(numout,*) 'WARNING : pondr, pond_inflow, pond_drainage', & & pond_reservoir_r(ig), pond_inflow_r(ig), pond_drainage_r(ig) CALL ipslerr_p(2, 'routing_simple_flow', 'WARNING negative slow_reservoir.','','') ENDIF ENDIF ENDDO DO ig=1,nbpt_r IF ( routing_mask_r(ig) ) THEN hydrographs_r(ig)=hydrographs_r(ig)+fast_flow_r(ig)+slow_flow_r(ig)+stream_flow_r(ig) ENDIF ENDDO ENDDO ! isplit ! ! ! ! The three output types for the routing : endoheric basins,, rivers and ! diffuse coastal flow. ! lakeinflow_r(:)=0 coastalflow_r(:)=0 riverflow_r(:)=0 basins_riverflow_mpi(:)=0 DO ig=1,nbpt_r IF ( routing_mask_r(ig) ) THEN IF (is_lakeinflow_r(ig)) THEN lakeinflow_r(ig) = transport_r(ig) basins_riverflow_mpi(basins_extended_r(ig)) = & basins_riverflow_mpi(basins_extended_r(ig))+lakeinflow_r(ig) ENDIF IF (is_coastalflow_r(ig)) THEN coastalflow_r(ig) = transport_r(ig) basins_riverflow_mpi(basins_extended_r(ig)) = & basins_riverflow_mpi(basins_extended_r(ig))+coastalflow_r(ig) ENDIF IF (is_riverflow_r(ig)) THEN riverflow_r(ig) = transport_r(ig) basins_riverflow_mpi(basins_extended_r(ig)) = & basins_riverflow_mpi(basins_extended_r(ig))+riverflow_r(ig) ENDIF ENDIF ENDDO CALL MPI_ALLREDUCE(basins_riverflow_mpi,basins_riverflow,basins_count+1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr) CALL xios_send_field("basins_riverflow",basins_riverflow(1:basins_out)/1000./dt_routing) ! flood_res = flood_diag + pond_diag ! CALL xios_send_field("routing_lakeinflow_r" ,lakeinflow_r*routing_weight_r) CALL xios_recv_field("routing_lakeinflow" ,lakeinflow) CALL xios_send_field("routing_coastalflow_r" ,coastalflow_r*routing_weight_coast_r) CALL xios_recv_field("routing_coastalflow_temp" ,coastalflow) WHERE(.NOT. coast_mask(:)) coastalflow(:)=0 CALL xios_send_field("routing_coastalflow" ,coastalflow) CALL xios_send_field("routing_riverflow_r" ,riverflow_r*routing_weight_coast_r) CALL xios_recv_field("routing_riverflow_temp" ,riverflow) WHERE(.NOT. coast_mask(:)) riverflow(:)=0 CALL xios_send_field("routing_riverflow" ,riverflow) CALL xios_send_field("out_flow",lakeinflow+coastalflow+riverflow) CALL xios_send_field("routing_hydrographs_r", hydrographs_r/1000./dt_routing) ! diag CALL xios_send_field("routing_fast_reservoir_r" , fast_reservoir_r) CALL xios_recv_field("routing_fast_reservoir" , fast_diag_mpi) CALL xios_send_field("routing_slow_reservoir_r" , slow_reservoir_r) CALL xios_recv_field("routing_slow_reservoir" , slow_diag_mpi) CALL xios_send_field("routing_stream_reservoir_r" , stream_reservoir_r) CALL xios_recv_field("routing_stream_reservoir" , stream_diag_mpi) CALL xios_recv_field("water_balance_after" ,water_balance_after) PRINT *,"routing water Balance ; before : ", water_balance_before," ; after : ",water_balance_after, & " ; delta : ", 100*(water_balance_after-water_balance_before)/(0.5*(water_balance_after+water_balance_before)),"%" !! Remove water from lake reservoir if it exceeds the maximum limit and distribute it !! uniformly over all possible the coastflow gridcells ! Calculate lake_overflow and remove it from lake_reservoir sum_lake_overflow=0 DO ig=1,nbp_mpi lake_overflow = MAX(0., lake_reservoir_mpi(ig) - max_lake_reservoir*area_mpi(ig)) lake_reservoir_mpi(ig) = lake_reservoir_mpi(ig) - lake_overflow sum_lake_overflow = sum_lake_overflow+lake_overflow END DO ! Calculate the sum of the lake_overflow and distribute it uniformly over all gridboxes CALL reduce_sum_mpi(sum_lake_overflow,total_lake_overflow) CALL bcast_mpi(total_lake_overflow) WHERE(coast_mask) coastalflow = coastalflow + total_lake_overflow/total_coast_points ENDIF ! is_omp_root CALL scatter_omp(riverflow,riverflow_omp) CALL scatter_omp(coastalflow,coastalflow_omp) CALL scatter_omp(lakeinflow,lakeinflow_omp) CALL scatter_omp(lake_reservoir_mpi,lake_reservoir) CALL scatter_omp(fast_diag_mpi,fast_diag) CALL scatter_omp(slow_diag_mpi,slow_diag) CALL scatter_omp(stream_diag_mpi,stream_diag) DO ig=1,nbpt fast_diag(ig)=fast_diag(ig)/area(ig) slow_diag(ig)=slow_diag(ig)/area(ig) stream_diag(ig)=stream_diag(ig)/area(ig) ENDDO flood_frac(:) = zero flood_height(:) = zero flood_frac_bas_r(:) = zero returnflow(:) = zero reinfiltration(:) = zero ! ! ! Compute the fluxes which leave the routing scheme ! ! Lakeinflow is in Kg/dt ! returnflow is in Kg/m^2/dt ! delsurfstor(:) = -flood_diag(:)-pond_diag(:)-lake_diag(:) hydrographs(:) = zero slowflow_diag(:) = zero flood_diag(:) = zero pond_diag(:) = zero irrigation(:) = zero END SUBROUTINE routing_simple_flow !! ============================================================================================================================= !! SUBROUTINE: routing_simple_initialize !! !>\BRIEF Initialize the routing_simple module !! !! DESCRIPTION: Initialize the routing_simple module. Read from restart file or read the routing.nc file to initialize the !! routing scheme. !! !! RECENT CHANGE(S) !! !! REFERENCE(S) !! !! FLOWCHART !! \n !_ ============================================================================================================================== SUBROUTINE routing_simple_initialize( kjit, nbpt, index, & rest_id, hist_id, hist2_id, lalo, & neighbours, resolution, contfrac, stempdiag, & returnflow, reinfiltration, irrigation, riverflow, & coastalflow, flood_frac, flood_res ) IMPLICIT NONE !! 0 Variable and parameter description !! 0.1 Input variables INTEGER(i_std), INTENT(in) :: kjit !! Time step number (unitless) INTEGER(i_std), INTENT(in) :: nbpt !! Domain size (unitless) INTEGER(i_std), INTENT(in) :: index(nbpt) !! Indices of the points on the map (unitless) INTEGER(i_std),INTENT(in) :: rest_id !! Restart file identifier (unitless) INTEGER(i_std),INTENT(in) :: hist_id !! Access to history file (unitless) INTEGER(i_std),INTENT(in) :: hist2_id !! Access to history file 2 (unitless) REAL(r_std), INTENT(in) :: lalo(nbpt,2) !! Vector of latitude and longitudes (beware of the order !) 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) REAL(r_std), INTENT(in) :: resolution(nbpt,2) !! The size of each grid box in X and Y (m) REAL(r_std), INTENT(in) :: contfrac(nbpt) !! Fraction of land in each grid box (unitless;0-1) REAL(r_std), INTENT(in) :: stempdiag(nbpt,nslm) !! Diagnostic soil temperature profile !! 0.2 Output variables REAL(r_std), INTENT(out) :: returnflow(nbpt) !! The water flow from lakes and swamps which returns to the grid box. !! This water will go back into the hydrol or hydrolc module to allow re-evaporation (kg/m^2/dt) REAL(r_std), INTENT(out) :: reinfiltration(nbpt) !! Water flow from ponds and floodplains which returns to the grid box (kg/m^2/dt) 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) 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) 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) REAL(r_std), INTENT(out) :: flood_frac(nbpt) !! Flooded fraction of the grid box (unitless;0-1) REAL(r_std), INTENT(out) :: flood_res(nbpt) !! Diagnostic of water amount in the floodplains reservoir (kg) !! 0.3 Local variables LOGICAL :: init_irrig !! Logical to initialize the irrigation (true/false) LOGICAL :: init_flood !! Logical to initialize the floodplains (true/false) LOGICAL :: init_swamp !! Logical to initialize the swamps (true/false) !_ ================================================================================================================================ ! ! do initialisation ! nbvmax = 440 ! Here we will allocate the memory and get the fixed fields from the restart file. ! If the info is not found then we will compute the routing map. ! CALL routing_simple_init_1 (kjit, nbpt, index, returnflow, reinfiltration, irrigation, & riverflow, coastalflow, flood_frac, flood_res, stempdiag, rest_id) routing_area => routing_area_loc topo_resid => topo_resid_loc route_togrid => route_togrid_loc route_tobasin => route_tobasin_loc global_basinid => global_basinid_loc hydrodiag => hydrodiag_loc ! This routine computes the routing map if the route_togrid_glo is undefined. This means that the ! map has not been initialized during the restart process.. ! !! Reads in the map of the basins and flow directions to construct the catchments of each grid box ! IF ( COUNT(route_togrid_glo .GE. undef_int) .GT. 0 ) THEN ! CALL routing_basins_p(nbpt, lalo, neighbours, resolution, contfrac) ENDIF ! ! Do we have what we need if we want to do irrigation !! Initialisation of flags for irrigated land, flood plains and swamps ! init_irrig = .FALSE. IF ( do_irrigation ) THEN IF (COUNT(irrigated .GE. undef_sechiba-1) > 0) init_irrig = .TRUE. END IF init_flood = .FALSE. IF ( do_floodplains ) THEN IF (COUNT(floodplains .GE. undef_sechiba-1) > 0) init_flood = .TRUE. END IF init_swamp = .FALSE. IF ( doswamps ) THEN IF (COUNT(swamp .GE. undef_sechiba-1) > 0 ) init_swamp = .TRUE. END IF !! If we have irrigated land, flood plains or swamps then we need to interpolate the 0.5 degree !! base data set to the resolution of the model. IF ( init_irrig .OR. init_flood .OR. init_swamp ) THEN ! CALL routing_irrigmap(nbpt, index, lalo, neighbours, resolution, & ! contfrac, init_irrig, irrigated, init_flood, floodplains, init_swamp, swamp, hist_id, hist2_id) ENDIF IF ( do_irrigation ) THEN CALL xios_orchidee_send_field("irrigmap",irrigated) WRITE(numout,*) 'Verification : range of irrigated : ', MINVAL(irrigated), MAXVAL(irrigated) IF ( .NOT. almaoutput ) THEN CALL histwrite_p(hist_id, 'irrigmap', 1, irrigated, nbpt, index) ELSE CALL histwrite_p(hist_id, 'IrrigationMap', 1, irrigated, nbpt, index) ENDIF IF ( hist2_id > 0 ) THEN IF ( .NOT. almaoutput ) THEN CALL histwrite_p(hist2_id, 'irrigmap', 1, irrigated, nbpt, index) ELSE CALL histwrite_p(hist2_id, 'IrrigationMap', 1, irrigated, nbpt, index) ENDIF ENDIF ENDIF IF ( do_floodplains ) THEN CALL xios_orchidee_send_field("floodmap",floodplains) WRITE(numout,*) 'Verification : range of floodplains : ', MINVAL(floodplains), MAXVAL(floodplains) IF ( .NOT. almaoutput ) THEN CALL histwrite_p(hist_id, 'floodmap', 1, floodplains, nbpt, index) ELSE CALL histwrite_p(hist_id, 'FloodplainsMap', 1, floodplains, nbpt, index) ENDIF IF ( hist2_id > 0 ) THEN IF ( .NOT. almaoutput ) THEN CALL histwrite_p(hist2_id, 'floodmap', 1, floodplains, nbpt, index) ELSE CALL histwrite_p(hist2_id, 'FloodplainsMap', 1, floodplains, nbpt, index) ENDIF ENDIF ENDIF IF ( doswamps ) THEN CALL xios_orchidee_send_field("swampmap",swamp) WRITE(numout,*) 'Verification : range of swamp : ', MINVAL(swamp), MAXVAL(swamp) IF ( .NOT. almaoutput ) THEN CALL histwrite_p(hist_id, 'swampmap', 1, swamp, nbpt, index) ELSE CALL histwrite_p(hist_id, 'SwampMap', 1, swamp, nbpt, index) ENDIF IF ( hist2_id > 0 ) THEN IF ( .NOT. almaoutput ) THEN CALL histwrite_p(hist2_id, 'swampmap', 1, swamp, nbpt, index) ELSE CALL histwrite_p(hist2_id, 'SwampMap', 1, swamp, nbpt, index) ENDIF ENDIF ENDIF !! This routine gives a diagnostic of the basins used. ! CALL routing_diagnostic_p(nbpt, index, lalo, resolution, contfrac, hist_id, hist2_id) CALL routing_simple_init_2(nbpt, contfrac) END SUBROUTINE routing_simple_initialize !! ================================================================================================================================ !! SUBROUTINE : routing_simple_main !! !>\BRIEF This module routes the water over the continents (runoff and !! drainage produced by the hydrolc or hydrol module) into the oceans. !! !! DESCRIPTION (definitions, functional, design, flags): !! The routing scheme (Polcher, 2003) carries the water from the runoff and drainage simulated by SECHIBA !! to the ocean through reservoirs, with some delay. The routing scheme is based on !! a parametrization of the water flow on a global scale (Miller et al., 1994; Hagemann !! and Dumenil, 1998). Given the global map of the main watersheds (Oki et al., 1999; !! Fekete et al., 1999; Vorosmarty et al., 2000) which delineates the boundaries of subbasins !! and gives the eight possible directions of water flow within the pixel, the surface !! runoff and the deep drainage are routed to the ocean. The time-step of the routing is one day. !! The scheme also diagnoses how much water is retained in the foodplains and thus return to soil !! moisture or is taken out of the rivers for irrigation. \n !! !! RECENT CHANGE(S): None !! !! MAIN OUTPUT VARIABLE(S): !! The result of the routing are 3 fluxes : !! - riverflow : The water which flows out from the major rivers. The flux will be located !! on the continental grid but this should be a coastal point. !! - coastalflow : This is the water which flows in a disperse way into the ocean. Essentially these !! are the outflows from all of the small rivers. !! - returnflow : This is the water which flows into a land-point - typically rivers which end in !! the desert. This water will go back into the hydrol module to allow re-evaporation. !! - irrigation : This is water taken from the reservoir and is being put into the upper !! layers of the soil. !! The two first fluxes are in kg/dt and the last two fluxes are in kg/(m^2dt).\n !! !! REFERENCE(S) : !! - Miller JR, Russell GL, Caliri G (1994) !! Continental-scale river flow in climate models. !! J. Clim., 7:914-928 !! - Hagemann S and Dumenil L. (1998) !! A parametrization of the lateral waterflow for the global scale. !! Clim. Dyn., 14:17-31 !! - Oki, T., T. Nishimura, and P. Dirmeyer (1999) !! Assessment of annual runoff from land surface models using total runoff integrating pathways (TRIP) !! J. Meteorol. Soc. Jpn., 77, 235-255 !! - Fekete BM, Charles V, Grabs W (2000) !! Global, composite runoff fields based on observed river discharge and simulated water balances. !! Technical report, UNH/GRDC, Global Runoff Data Centre, Koblenz !! - Vorosmarty, C., B. Fekete, B. Meybeck, and R. Lammers (2000) !! Global system of rivers: Its role in organizing continental land mass and defining land-to-ocean linkages !! Global Biogeochem. Cycles, 14, 599-621 !! - Vivant, A-C. (?? 2002) !! Développement du schéma de routage et des plaines d'inondation, MSc Thesis, Paris VI University !! - J. Polcher (2003) !! Les processus de surface a l'echelle globale et leurs interactions avec l'atmosphere !! Habilitation a diriger les recherches, Paris VI University, 67pp. !! !! FLOWCHART : !! \latexonly !! \includegraphics[scale=0.75]{routing_main_flowchart.png} !! \endlatexonly !! \n !_ ================================================================================================================================ SUBROUTINE routing_simple_main(kjit, nbpt, index, & lalo, neighbours, resolution, contfrac, totfrac_nobio, veget_max, floodout, runoff, & drainage, transpot, precip_rain, humrel, k_litt, flood_frac, flood_res, & stempdiag, reinf_slope, returnflow, reinfiltration, irrigation, riverflow, coastalflow, & rest_id, hist_id, hist2_id) IMPLICIT NONE !! 0 Variable and parameter description !! 0.1 Input variables INTEGER(i_std), INTENT(in) :: kjit !! Time step number (unitless) INTEGER(i_std), INTENT(in) :: nbpt !! Domain size (unitless) INTEGER(i_std),INTENT(in) :: rest_id !! Restart file identifier (unitless) INTEGER(i_std),INTENT(in) :: hist_id !! Access to history file (unitless) INTEGER(i_std),INTENT(in) :: hist2_id !! Access to history file 2 (unitless) INTEGER(i_std), INTENT(in) :: index(nbpt) !! Indices of the points on the map (unitless) REAL(r_std), INTENT(in) :: lalo(nbpt,2) !! Vector of latitude and longitudes (beware of the order !) 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) REAL(r_std), INTENT(in) :: resolution(nbpt,2) !! The size of each grid box in X and Y (m) REAL(r_std), INTENT(in) :: contfrac(nbpt) !! Fraction of land in each grid box (unitless;0-1) REAL(r_std), INTENT(in) :: totfrac_nobio(nbpt) !! Total fraction of no-vegetation (continental ice, lakes ...) (unitless;0-1) REAL(r_std), INTENT(in) :: veget_max(nbpt,nvm) !! Maximal fraction of vegetation (unitless;0-1) REAL(r_std), INTENT(in) :: floodout(nbpt) !! Grid-point flow out of floodplains (kg/m^2/dt) REAL(r_std), INTENT(in) :: runoff(nbpt) !! Grid-point runoff (kg/m^2/dt) REAL(r_std), INTENT(in) :: drainage(nbpt) !! Grid-point drainage (kg/m^2/dt) REAL(r_std), INTENT(in) :: transpot(nbpt,nvm) !! Potential transpiration of the vegetation (kg/m^2/dt) REAL(r_std), INTENT(in) :: precip_rain(nbpt) !! Rainfall (kg/m^2/dt) REAL(r_std), INTENT(in) :: k_litt(nbpt) !! Averaged conductivity for saturated infiltration in the 'litter' layer (kg/m^2/dt) REAL(r_std), INTENT(in) :: humrel(nbpt,nvm) !! Soil moisture stress, root extraction potential (unitless) REAL(r_std), INTENT(in) :: stempdiag(nbpt,nslm) !! Diagnostic soil temperature profile 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) !! 0.2 Output variables REAL(r_std), INTENT(out) :: returnflow(nbpt) !! The water flow from lakes and swamps which returns to the grid box. !! This water will go back into the hydrol or hydrolc module to allow re-evaporation (kg/m^2/dt) REAL(r_std), INTENT(out) :: reinfiltration(nbpt) !! Water flow from ponds and floodplains which returns to the grid box (kg/m^2/dt) 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) 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) 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) REAL(r_std), INTENT(out) :: flood_frac(nbpt) !! Flooded fraction of the grid box (unitless;0-1) REAL(r_std), INTENT(out) :: flood_res(nbpt) !! Diagnostic of water amount in the floodplains reservoir (kg) !! 0.3 Local variables CHARACTER(LEN=30) :: var_name !! To store variables names for I/O (unitless) REAL(r_std), DIMENSION(1) :: tmp_day !! REAL(r_std), DIMENSION(nbpt) :: return_lakes !! Water from lakes flowing back into soil moisture (kg/m^2/dt) INTEGER(i_std) :: ig, jv !! Indices (unitless) REAL(r_std), DIMENSION(nbpt) :: tot_vegfrac_nowoody !! Total fraction occupied by grass (0-1,unitless) REAL(r_std),PARAMETER :: delta_lon=(360./144.)/2. REAL(r_std),PARAMETER :: delta_lat=(180./144.)/2. !_ ================================================================================================================================ ! !! Computes the variables averaged between routing time steps and which will be used in subsequent calculations ! floodout_mean(:) = floodout_mean(:) + floodout(:) runoff_mean(:) = runoff_mean(:) + runoff(:) drainage_mean(:) = drainage_mean(:) + drainage(:) floodtemp(:) = stempdiag(:,floodtemp_lev) precip_mean(:) = precip_mean(:) + precip_rain(:) ! !! Computes the total fraction occupied by the grasses and the crops for each grid cell tot_vegfrac_nowoody(:) = zero DO jv = 1, nvm IF ( (jv /= ibare_sechiba) .AND. .NOT.(is_tree(jv)) ) THEN tot_vegfrac_nowoody(:) = tot_vegfrac_nowoody(:) + veget_max(:,jv) END IF END DO DO ig = 1, nbpt IF ( tot_vegfrac_nowoody(ig) .GT. min_sechiba ) THEN DO jv = 1,nvm IF ( (jv /= ibare_sechiba) .AND. .NOT.(is_tree(jv)) ) THEN transpot_mean(ig) = transpot_mean(ig) + transpot(ig,jv) * veget_max(ig,jv)/tot_vegfrac_nowoody(ig) END IF END DO ELSE IF (MAXVAL(veget_max(ig,2:nvm)) .GT. min_sechiba) THEN DO jv = 2, nvm transpot_mean(ig) = transpot_mean(ig) + transpot(ig,jv) * veget_max(ig,jv)/ SUM(veget_max(ig,2:nvm)) ENDDO ENDIF ENDIF ENDDO ! ! Averaged variables (i.e. *dt_sechiba/dt_routing). This accounts for the difference between the shorter ! timestep dt_sechiba of other parts of the model and the long dt_routing timestep (set to one day at present) ! totnobio_mean(:) = totnobio_mean(:) + totfrac_nobio(:)*dt_sechiba/dt_routing k_litt_mean(:) = k_litt_mean(:) + k_litt(:)*dt_sechiba/dt_routing ! ! Only potentially vegetated surfaces are taken into account. At the start of ! the growing seasons we will give more weight to these areas. ! DO jv=2,nvm DO ig=1,nbpt humrel_mean(ig) = humrel_mean(ig) + humrel(ig,jv)*veget_max(ig,jv)*dt_sechiba/dt_routing vegtot_mean(ig) = vegtot_mean(ig) + veget_max(ig,jv)*dt_sechiba/dt_routing ENDDO ENDDO ! time_counter = time_counter + dt_sechiba ! ! If the time has come we do the routing. ! IF ( NINT(time_counter) .GE. NINT(dt_routing) ) THEN ! ! Check the water balance if needed ! ! IF ( check_waterbal ) THEN ! CALL routing_waterbal(nbpt, .TRUE., floodout_mean, runoff_mean, drainage_mean, returnflow_mean, & ! & reinfiltration_mean, irrigation_mean, riverflow_mean, coastalflow_mean) ! ENDIF ! ! Make sure we do not flood north of 49N as there freezing processes start to play a role and they ! are not yet well treated in ORCHIDEE. ! DO ig=1,nbpt IF ( lalo(ig,1) > 49.0 ) THEN floodtemp(ig) = tp_00 - un ENDIF ENDDO ! !! Computes the transport of water in the various reservoirs ! !ym runoff_mean=0 !ym drainage_mean=0 !ym DO ig=1,nbpt !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 !ym ENDDO CALL routing_simple_flow(nbpt, dt_routing, lalo, floodout_mean, runoff_mean, drainage_mean, & & vegtot_mean, totnobio_mean, transpot_mean, precip_mean, humrel_mean, k_litt_mean, floodtemp, reinf_slope, & & lakeinflow_mean, returnflow_mean, reinfiltration_mean, irrigation_mean, riverflow_mean, & & coastalflow_mean, hydrographs, slowflow_diag, flood_frac, flood_res) ! !! Responsible for storing the water in lakes ! CALL routing_simple_lake(nbpt, dt_routing, lakeinflow_mean, humrel_mean, contfrac, return_lakes) ! returnflow_mean(:) = returnflow_mean(:) + return_lakes(:) ! !! Check the water balance in the routing scheme ! ! IF ( check_waterbal ) THEN ! CALL routing_waterbal(nbpt, .FALSE., floodout_mean, runoff_mean, drainage_mean, returnflow_mean, & ! & reinfiltration_mean, irrigation_mean, riverflow_mean, coastalflow_mean) ! ENDIF ! time_counter = zero ! floodout_mean(:) = zero runoff_mean(:) = zero drainage_mean(:) = zero transpot_mean(:) = zero precip_mean(:) = zero ! humrel_mean(:) = zero totnobio_mean(:) = zero k_litt_mean(:) = zero vegtot_mean(:) = zero ! ! Change the units of the routing fluxes from kg/dt_routing into kg/dt_sechiba ! and from m^3/dt_routing into m^3/dt_sechiba ! returnflow_mean(:) = returnflow_mean(:)/dt_routing*dt_sechiba reinfiltration_mean(:) = reinfiltration_mean(:)/dt_routing*dt_sechiba irrigation_mean(:) = irrigation_mean(:)/dt_routing*dt_sechiba irrig_netereq(:) = irrig_netereq(:)/dt_routing*dt_sechiba ! ! Change units as above but at the same time transform the kg/dt_sechiba to m^3/dt_sechiba ! riverflow_mean(:) = riverflow_mean(:)/dt_routing*dt_sechiba/mille coastalflow_mean(:) = coastalflow_mean(:)/dt_routing*dt_sechiba/mille hydrographs(:) = hydrographs(:)/dt_routing*dt_sechiba/mille slowflow_diag(:) = slowflow_diag(:)/dt_routing*dt_sechiba/mille ENDIF ! ! Return the fraction of routed water for this time step. ! returnflow(:) = returnflow_mean(:) reinfiltration(:) = reinfiltration_mean(:) irrigation(:) = irrigation_mean(:) riverflow(:) = riverflow_mean(:) coastalflow(:) = coastalflow_mean(:) ! ! Write diagnostics ! CALL xios_orchidee_send_field("riversret",returnflow*one_day/dt_sechiba) CALL xios_orchidee_send_field("hydrographs",hydrographs/dt_sechiba) CALL xios_orchidee_send_field("slowflow",slowflow_diag/dt_sechiba) ! Qb in m3/s CALL xios_orchidee_send_field("reinfiltration",reinfiltration) CALL xios_orchidee_send_field("fastr",fast_diag) CALL xios_orchidee_send_field("slowr",slow_diag) CALL xios_orchidee_send_field("streamr",stream_diag) CALL xios_orchidee_send_field("laker",lake_diag) CALL xios_orchidee_send_field("pondr",pond_diag) CALL xios_orchidee_send_field("irrigation",irrigation*one_day/dt_sechiba) CALL xios_orchidee_send_field("netirrig",irrig_netereq*one_day/dt_sechiba) CALL xios_orchidee_send_field("floodh",flood_height) CALL xios_orchidee_send_field("floodr",flood_diag) delsurfstor = delsurfstor + flood_diag + pond_diag + lake_diag CALL xios_orchidee_send_field("SurfStor",flood_diag+pond_diag+lake_diag) ! Transform from kg/dt_sechiba into m^3/s CALL xios_orchidee_send_field("hydrographs",hydrographs/mille/dt_sechiba) CALL xios_orchidee_send_field("slowflow",slowflow_diag/mille/dt_sechiba) ! previous id name: Qb CALL xios_orchidee_send_field("coastalflow",coastalflow/dt_sechiba) CALL xios_orchidee_send_field("riverflow",riverflow/dt_sechiba) IF ( .NOT. almaoutput ) THEN ! CALL histwrite_p(hist_id, 'riversret', kjit, returnflow, nbpt, index) IF (do_floodplains .OR. doponds) THEN CALL histwrite_p(hist_id, 'reinfiltration', kjit, reinfiltration, nbpt, index) ENDIF CALL histwrite_p(hist_id, 'hydrographs', kjit, hydrographs, nbpt, index) ! CALL histwrite_p(hist_id, 'fastr', kjit, fast_diag, nbpt, index) CALL histwrite_p(hist_id, 'slowr', kjit, slow_diag, nbpt, index) CALL histwrite_p(hist_id, 'streamr', kjit, stream_diag, nbpt, index) IF ( do_floodplains ) THEN CALL histwrite_p(hist_id, 'floodr', kjit, flood_diag, nbpt, index) CALL histwrite_p(hist_id, 'floodh', kjit, flood_height, nbpt, index) ENDIF CALL histwrite_p(hist_id, 'pondr', kjit, pond_diag, nbpt, index) CALL histwrite_p(hist_id, 'lakevol', kjit, lake_diag, nbpt, index) ! IF ( do_irrigation ) THEN CALL histwrite_p(hist_id, 'irrigation', kjit, irrigation, nbpt, index) CALL histwrite_p(hist_id, 'returnflow', kjit, returnflow, nbpt, index) CALL histwrite_p(hist_id, 'netirrig', kjit, irrig_netereq, nbpt, index) ENDIF ! ELSE ! delsurfstor = delsurfstor + flood_diag + pond_diag + lake_diag CALL histwrite_p(hist_id, 'DelSurfStor', kjit, delsurfstor, nbpt, index) CALL histwrite_p(hist_id, 'SurfStor', kjit, flood_diag+pond_diag+lake_diag, nbpt, index) CALL histwrite_p(hist_id, 'Dis', kjit, hydrographs, nbpt, index) IF ( do_irrigation ) THEN CALL histwrite_p(hist_id, 'Qirrig', kjit, irrigation, nbpt, index) CALL histwrite_p(hist_id, 'Qirrig_req', kjit, irrig_netereq, nbpt, index) ENDIF ! ENDIF IF ( hist2_id > 0 ) THEN IF ( .NOT. almaoutput ) THEN ! CALL histwrite_p(hist2_id, 'riversret', kjit, returnflow, nbpt, index) IF (do_floodplains .OR. doponds) THEN CALL histwrite_p(hist2_id, 'reinfiltration', kjit, reinfiltration, nbpt, index) ENDIF CALL histwrite_p(hist2_id, 'hydrographs', kjit, hydrographs, nbpt, index) ! CALL histwrite_p(hist2_id, 'fastr', kjit, fast_diag, nbpt, index) CALL histwrite_p(hist2_id, 'slowr', kjit, slow_diag, nbpt, index) IF ( do_floodplains ) THEN CALL histwrite_p(hist2_id, 'floodr', kjit, flood_diag, nbpt, index) CALL histwrite_p(hist2_id, 'floodh', kjit, flood_height, nbpt, index) ENDIF CALL histwrite_p(hist2_id, 'pondr', kjit, pond_diag, nbpt, index) CALL histwrite_p(hist2_id, 'streamr', kjit, stream_diag, nbpt, index) CALL histwrite_p(hist2_id, 'lakevol', kjit, lake_diag, nbpt, index) ! IF ( do_irrigation ) THEN CALL histwrite_p(hist2_id, 'irrigation', kjit, irrigation, nbpt, index) CALL histwrite_p(hist2_id, 'returnflow', kjit, returnflow, nbpt, index) CALL histwrite_p(hist2_id, 'netirrig', kjit, irrig_netereq, nbpt, index) ENDIF ! ELSE ! delsurfstor=delsurfstor + flood_diag + pond_diag + lake_diag CALL histwrite_p(hist2_id, 'DelSurfStor', kjit, delsurfstor, nbpt, index) CALL histwrite_p(hist2_id, 'SurfStor', kjit, flood_diag+pond_diag+lake_diag, nbpt, index) CALL histwrite_p(hist2_id, 'Dis', kjit, hydrographs, nbpt, index) ! ENDIF ENDIF ! ! END SUBROUTINE routing_simple_main !! ============================================================================================================================= !! SUBROUTINE: routing_simple_finalize !! !>\BRIEF Write to restart file !! !! DESCRIPTION: Write module variables to restart file !! !! RECENT CHANGE(S) !! !! REFERENCE(S) !! !! FLOWCHART !! \n !_ ============================================================================================================================== SUBROUTINE routing_simple_finalize( kjit, nbpt, rest_id, flood_frac, flood_res ) USE xios IMPLICIT NONE !! 0 Variable and parameter description !! 0.1 Input variables INTEGER(i_std), INTENT(in) :: kjit !! Time step number (unitless) INTEGER(i_std), INTENT(in) :: nbpt !! Domain size (unitless) INTEGER(i_std),INTENT(in) :: rest_id !! Restart file identifier (unitless) REAL(r_std), INTENT(in) :: flood_frac(nbpt) !! Flooded fraction of the grid box (unitless;0-1) REAL(r_std), INTENT(in) :: flood_res(nbpt) !! Diagnostic of water amount in the floodplains reservoir (kg) !! 0.2 Local variables REAL(r_std), DIMENSION(1) :: tmp_day !_ ================================================================================================================================ IF (is_omp_root) THEN CALL xios_send_field("fast_reservoir_restart",fast_reservoir_r) CALL xios_send_field("slow_reservoir_restart",slow_reservoir_r) CALL xios_send_field("stream_reservoir_restart",stream_reservoir_r) ENDIF ! ! Write restart variables ! tmp_day(1) = time_counter IF (is_root_prc) CALL restput (rest_id, 'routingcounter', 1, 1, 1, kjit, tmp_day) CALL restput_p (rest_id, 'routingarea', nbp_glo, nbasmax, 1, kjit, routing_area, 'scatter', nbp_glo, index_g) CALL restput_p (rest_id, 'routetogrid', nbp_glo, nbasmax, 1, kjit, REAL(route_togrid,r_std), 'scatter', & nbp_glo, index_g) CALL restput_p (rest_id, 'routetobasin', nbp_glo, nbasmax, 1, kjit, REAL(route_tobasin,r_std), 'scatter', & nbp_glo, index_g) CALL restput_p (rest_id, 'basinid', nbp_glo, nbasmax, 1, kjit, REAL(global_basinid,r_std), 'scatter', & nbp_glo, index_g) CALL restput_p (rest_id, 'topoindex', nbp_glo, nbasmax, 1, kjit, topo_resid, 'scatter', nbp_glo, index_g) CALL restput_p (rest_id, 'fastres', nbp_glo, nbasmax, 1, kjit, fast_reservoir, 'scatter', nbp_glo, index_g) CALL restput_p (rest_id, 'slowres', nbp_glo, nbasmax, 1, kjit, slow_reservoir, 'scatter', nbp_glo, index_g) CALL restput_p (rest_id, 'streamres', nbp_glo, nbasmax, 1, kjit, stream_reservoir, 'scatter',nbp_glo,index_g) CALL restput_p (rest_id, 'floodres', nbp_glo, nbasmax, 1, kjit, flood_reservoir, 'scatter', nbp_glo, index_g) CALL restput_p (rest_id, 'floodh', nbp_glo, 1, 1, kjit, flood_height, 'scatter', nbp_glo, index_g) CALL restput_p (rest_id, 'flood_frac_bas', nbp_glo, nbasmax, 1, kjit, flood_frac_bas, 'scatter', nbp_glo, index_g) CALL restput_p (rest_id, 'pond_frac', nbp_glo, 1, 1, kjit, pond_frac, 'scatter', nbp_glo, index_g) CALL restput_p (rest_id, 'flood_frac', nbp_glo, 1, 1, kjit, flood_frac, 'scatter', nbp_glo, index_g) CALL restput_p (rest_id, 'flood_res', nbp_glo, 1, 1, kjit, flood_res, 'scatter', nbp_glo, index_g) CALL restput_p (rest_id, 'lakeres', nbp_glo, 1, 1, kjit, lake_reservoir, 'scatter', nbp_glo, index_g) CALL restput_p (rest_id, 'pondres', nbp_glo, 1, 1, kjit, pond_reservoir, 'scatter', nbp_glo, index_g) CALL restput_p (rest_id, 'lakeinflow', nbp_glo, 1, 1, kjit, lakeinflow_mean, 'scatter', nbp_glo, index_g) CALL restput_p (rest_id, 'returnflow', nbp_glo, 1, 1, kjit, returnflow_mean, 'scatter', nbp_glo, index_g) CALL restput_p (rest_id, 'reinfiltration', nbp_glo, 1, 1, kjit, reinfiltration_mean, 'scatter', nbp_glo, index_g) CALL restput_p (rest_id, 'riverflow', nbp_glo, 1, 1, kjit, riverflow_mean, 'scatter', nbp_glo, index_g) CALL restput_p (rest_id, 'coastalflow', nbp_glo, 1, 1, kjit, coastalflow_mean, 'scatter', nbp_glo, index_g) CALL restput_p (rest_id, 'hydrographs', nbp_glo, 1, 1, kjit, hydrographs, 'scatter', nbp_glo, index_g) CALL restput_p (rest_id, 'slowflow_diag', nbp_glo, 1, 1, kjit, slowflow_diag, 'scatter', nbp_glo, index_g) ! ! Keep track of the accumulated variables ! CALL restput_p (rest_id, 'floodout_route', nbp_glo, 1, 1, kjit, floodout_mean, 'scatter', nbp_glo, index_g) CALL restput_p (rest_id, 'runoff_route', nbp_glo, 1, 1, kjit, runoff_mean, 'scatter', nbp_glo, index_g) CALL restput_p (rest_id, 'drainage_route', nbp_glo, 1, 1, kjit, drainage_mean, 'scatter', nbp_glo, index_g) CALL restput_p (rest_id, 'transpot_route', nbp_glo, 1, 1, kjit, transpot_mean, 'scatter', nbp_glo, index_g) CALL restput_p (rest_id, 'precip_route', nbp_glo, 1, 1, kjit, precip_mean, 'scatter', nbp_glo, index_g) CALL restput_p (rest_id, 'humrel_route', nbp_glo, 1, 1, kjit, humrel_mean, 'scatter', nbp_glo, index_g) CALL restput_p (rest_id, 'totnobio_route', nbp_glo, 1, 1, kjit, totnobio_mean, 'scatter', nbp_glo, index_g) CALL restput_p (rest_id, 'k_litt_route', nbp_glo, 1, 1, kjit, k_litt_mean, 'scatter', nbp_glo, index_g) CALL restput_p (rest_id, 'vegtot_route', nbp_glo, 1, 1, kjit, vegtot_mean, 'scatter', nbp_glo, index_g) IF ( do_irrigation ) THEN CALL restput_p (rest_id, 'irrigated', nbp_glo, 1, 1, kjit, irrigated, 'scatter', nbp_glo, index_g) CALL restput_p (rest_id, 'irrigation', nbp_glo, 1, 1, kjit, irrigation_mean, 'scatter', nbp_glo, index_g) ENDIF IF ( do_floodplains ) THEN CALL restput_p (rest_id, 'floodplains', nbp_glo, 1, 1, kjit, floodplains, 'scatter', nbp_glo, index_g) ENDIF IF ( doswamps ) THEN CALL restput_p (rest_id, 'swamp', nbp_glo, 1, 1, kjit, swamp, 'scatter', nbp_glo, index_g) ENDIF END SUBROUTINE routing_simple_finalize !! ================================================================================================================================ !! SUBROUTINE : routing_simple_init_1 !! !>\BRIEF This subroutine allocates the memory and get the fixed fields from the restart file. !! !! DESCRIPTION: Privat subroutine to the module routing_simple. This subroutine is called in the begining !! of routing_simple_initialize !! !! RECENT CHANGE(S): None !! !! MAIN OUTPUT VARIABLE(S): !! !! REFERENCES : None !! !! FLOWCHART :None !! \n !_ ================================================================================================================================ SUBROUTINE routing_simple_init_1(kjit, nbpt, index, returnflow, reinfiltration, irrigation, & riverflow, coastalflow, flood_frac, flood_res, stempdiag, rest_id) IMPLICIT NONE !! 0 Variable and parameter description !! 0.1 Input variables INTEGER(i_std), INTENT(in) :: kjit !! Time step number (unitless) INTEGER(i_std), INTENT(in) :: nbpt !! Domain size (unitless) INTEGER(i_std), DIMENSION (nbpt), INTENT(in) :: index !! Indices of the points on the map (unitless) REAL(r_std), DIMENSION(nbpt,nslm),INTENT(in) :: stempdiag !! Temperature profile in soil INTEGER(i_std), INTENT(in) :: rest_id !! Restart file identifier (unitless) !! 0.2 Output variables REAL(r_std), DIMENSION (nbpt),INTENT(out) :: returnflow !! The water flow from lakes and swamps which returns into the grid box. !! This water will go back into the hydrol or hydrolc module to allow re-evaporation (kg/m^2/dt) REAL(r_std), DIMENSION (nbpt),INTENT(out) :: reinfiltration !! Water flow from ponds and floodplains which returns to the grid box (kg/m^2/dt) 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) 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) 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) REAL(r_std), DIMENSION (nbpt),INTENT(out) :: flood_frac !! Flooded fraction of the grid box (unitless;0-1) REAL(r_std), DIMENSION (nbpt),INTENT(out) :: flood_res !! Diagnostic of water amount in the floodplains reservoir (kg) !! 0.3 Local variables CHARACTER(LEN=80) :: var_name !! To store variables names for I/O (unitless) REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: tmp_real_g !! A temporary real array for the integers REAL(r_std), DIMENSION(1) :: tmp_day !! REAL(r_std) :: ratio !! Diagnostic ratio to check that dt_routing is a multiple of dt_sechiba (unitless) REAL(r_std) :: totarea !! Total area of basin (m^2) INTEGER(i_std) :: ier, ig, ib, ipn(1) !! Indices (unitless) !_ ================================================================================================================================ ! ! ! These variables will require the configuration infrastructure ! !Config Key = DT_ROUTING !Config If = RIVER_ROUTING !Config Desc = Time step of the routing scheme !Config Def = one_day !Config Help = This values gives the time step in seconds of the routing scheme. !Config It should be multiple of the main time step of ORCHIDEE. One day !Config is a good value. !Config Units = [seconds] ! dt_routing = one_day CALL getin_p('DT_ROUTING', dt_routing) ! !Config Key = ROUTING_RIVERS !Config If = RIVER_ROUTING !Config Desc = Number of rivers !Config Def = 50 !Config Help = This parameter chooses the number of largest river basins !Config which should be treated as independently as rivers and not !Config flow into the oceans as diffusion coastal flow. !Config Units = [-] num_largest = 50 CALL getin_p('ROUTING_RIVERS', num_largest) ! !Config Key = DO_FLOODINFILT !Config Desc = Should floodplains reinfiltrate into the soil !Config If = RIVER_ROUTING !Config Def = n !Config Help = This parameters allows the user to ask the model !Config to take into account the flood plains reinfiltration !Config into the soil moisture. It then can go !Config back to the slow and fast reservoirs !Config Units = [FLAG] ! dofloodinfilt = .FALSE. CALL getin_p('DO_FLOODINFILT', dofloodinfilt) ! !Config Key = DO_SWAMPS !Config Desc = Should we include swamp parameterization !Config If = RIVER_ROUTING !Config Def = n !Config Help = This parameters allows the user to ask the model !Config to take into account the swamps and return !Config the water into the bottom of the soil. It then can go !Config back to the atmopshere. This tried to simulate !Config internal deltas of rivers. !Config Units = [FLAG] ! doswamps = .FALSE. CALL getin_p('DO_SWAMPS', doswamps) ! !Config Key = DO_PONDS !Config Desc = Should we include ponds !Config If = RIVER_ROUTING !Config Def = n !Config Help = This parameters allows the user to ask the model !Config to take into account the ponds and return !Config the water into the soil moisture. It then can go !Config back to the atmopshere. This tried to simulate !Config little ponds especially in West Africa. !Config Units = [FLAG] ! doponds = .FALSE. CALL getin_p('DO_PONDS', doponds) !Config Key = SLOW_TCST !Config Desc = Time constant for the slow reservoir !Config If = RIVER_ROUTING !Config Def = 25.0 !Config Help = This parameters allows the user to fix the !Config time constant (in days) of the slow reservoir !Config in order to get better river flows for !Config particular regions. !Config Units = [days] ! !> A value for property of each reservoir (in day/m) is given to compute a time constant (in day) !> for each reservoir (product of tcst and topo_resid). !> The value of tcst has been calibrated for the three reservoirs over the Senegal river basin only, !> during the 1 degree NCEP Corrected by Cru (NCC) resolution simulations (Ngo-Duc et al., 2005, Ngo-Duc et al., 2006) and !> generalized for all the basins of the world. The "slow reservoir" and the "fast reservoir" !> have the highest value in order to simulate the groundwater. !> The "stream reservoir", which represents all the water of the stream, has the lowest value. !> Those figures are the same for all the basins of the world. CALL getin_p('SLOW_TCST', slow_tcst) ! !Config Key = FAST_TCST !Config Desc = Time constant for the fast reservoir !Config If = RIVER_ROUTING !Config Def = 3.0 !Config Help = This parameters allows the user to fix the !Config time constant (in days) of the fast reservoir !Config in order to get better river flows for !Config particular regions. !Config Units = [days] CALL getin_p('FAST_TCST', fast_tcst) ! !Config Key = STREAM_TCST !Config Desc = Time constant for the stream reservoir !Config If = RIVER_ROUTING !Config Def = 0.24 !Config Help = This parameters allows the user to fix the !Config time constant (in days) of the stream reservoir !Config in order to get better river flows for !Config particular regions. !Config Units = [days] CALL getin_p('STREAM_TCST', stream_tcst) ! !Config Key = FLOOD_TCST !Config Desc = Time constant for the flood reservoir !Config If = RIVER_ROUTING !Config Def = 4.0 !Config Help = This parameters allows the user to fix the !Config time constant (in days) of the flood reservoir !Config in order to get better river flows for !Config particular regions. !Config Units = [days] CALL getin_p('FLOOD_TCST', flood_tcst) !Config Key = SWAMP_CST !Config Desc = Fraction of the river that flows back to swamps !Config If = RIVER_ROUTING !Config Def = 0.2 !Config Help = This parameters allows the user to fix the !Config fraction of the river transport !Config that flows to swamps !Config Units = [-] CALL getin_p('SWAMP_CST', swamp_cst) !Config Key = FLOOD_BETA !Config Desc = Parameter to fix the shape of the floodplain !Config If = RIVER_ROUTING !Config Def = 2.0 !Config Help = Parameter to fix the shape of the floodplain !Config (>1 for convex edges, <1 for concave edges) !Config Units = [-] CALL getin_p("FLOOD_BETA", beta) !Config Key = POND_BETAP !Config Desc = Ratio of the basin surface intercepted by ponds and the maximum surface of ponds !Config If = RIVER_ROUTING !Config Def = 0.5 !Config Help = !Config Units = [-] CALL getin_p("POND_BETAP", betap) !Config Key = FLOOD_CRI !Config Desc = Potential height for which all the basin is flooded !Config If = DO_FLOODPLAINS or DO_PONDS !Config Def = 2000. !Config Help = !Config Units = [mm] CALL getin_p("FLOOD_CRI", floodcri) !Config Key = POND_CRI !Config Desc = Potential height for which all the basin is a pond !Config If = DO_FLOODPLAINS or DO_PONDS !Config Def = 2000. !Config Help = !Config Units = [mm] CALL getin_p("POND_CRI", pondcri) !Config Key = MAX_LAKE_RESERVOIR !Config Desc = Maximum limit of water in lake_reservoir !Config If = RIVER_ROUTING !Config Def = 7000 !Config Help = !Config Units = [kg/m2(routing area)] max_lake_reservoir = 7000 CALL getin_p("MAX_LAKE_RESERVOIR", max_lake_reservoir) ! In order to simplify the time cascade check that dt_routing ! is a multiple of dt_sechiba ! ratio = dt_routing/dt_sechiba IF ( ABS(NINT(ratio) - ratio) .GT. 10*EPSILON(ratio)) THEN WRITE(numout,*) 'WARNING -- WARNING -- WARNING -- WARNING' WRITE(numout,*) "The chosen time step for the routing is not a multiple of the" WRITE(numout,*) "main time step of the model. We will change dt_routing so that" WRITE(numout,*) "this condition os fulfilled" dt_routing = NINT(ratio) * dt_sechiba WRITE(numout,*) 'THE NEW DT_ROUTING IS : ', dt_routing ENDIF ! IF ( dt_routing .LT. dt_sechiba) THEN WRITE(numout,*) 'WARNING -- WARNING -- WARNING -- WARNING' WRITE(numout,*) 'The routing timestep can not be smaller than the one' WRITE(numout,*) 'of the model. We reset its value to the model''s timestep.' WRITE(numout,*) 'The old DT_ROUTING is : ', dt_routing dt_routing = dt_sechiba WRITE(numout,*) 'THE NEW DT_ROUTING IS : ', dt_routing ENDIF ! var_name ="routingcounter" IF (is_root_prc) THEN CALL ioconf_setatt('UNITS', 's') CALL ioconf_setatt('LONG_NAME','Time counter for the routing scheme') CALL restget (rest_id, var_name, 1, 1, 1, kjit, .TRUE., tmp_day) IF (tmp_day(1) == val_exp) THEN time_counter = zero ELSE time_counter = tmp_day(1) ENDIF ENDIF CALL bcast(time_counter) ALLOCATE (routing_area_loc(nbpt,nbasmax), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for routing_area_loc','','') ALLOCATE (routing_area_glo(nbp_glo,nbasmax)) IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for routing_area_glo','','') var_name = 'routingarea' IF (is_root_prc) THEN CALL ioconf_setatt('UNITS', 'm^2') CALL ioconf_setatt('LONG_NAME','Area of basin') CALL restget (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., routing_area_glo, "gather", nbp_glo, index_g) ENDIF CALL scatter(routing_area_glo,routing_area_loc) routing_area=>routing_area_loc ALLOCATE (tmp_real_g(nbp_glo,nbasmax), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for tmp_real_g','','') ALLOCATE (route_togrid_loc(nbpt,nbasmax), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for route_togrid_loc','','') ALLOCATE (route_togrid_glo(nbp_glo,nbasmax), stat=ier) ! used in global in routing_simple_flow IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for route_togrid_glo','','') IF (is_root_prc) THEN var_name = 'routetogrid' CALL ioconf_setatt('UNITS', '-') CALL ioconf_setatt('LONG_NAME','Grid into which the basin flows') CALL restget (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., tmp_real_g, "gather", nbp_glo, index_g) route_togrid_glo(:,:) = undef_int WHERE ( tmp_real_g .LT. val_exp ) route_togrid_glo = NINT(tmp_real_g) ENDWHERE ENDIF CALL bcast(route_togrid_glo) ! used in global in routing_simple_flow CALL scatter(route_togrid_glo,route_togrid_loc) route_togrid=>route_togrid_loc ! ALLOCATE (route_tobasin_loc(nbpt,nbasmax), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for route_tobasin_loc','','') ALLOCATE (route_tobasin_glo(nbp_glo,nbasmax), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for route_tobasin_glo','','') IF (is_root_prc) THEN var_name = 'routetobasin' CALL ioconf_setatt('UNITS', '-') CALL ioconf_setatt('LONG_NAME','Basin in to which the water goes') CALL restget (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., tmp_real_g, "gather", nbp_glo, index_g) route_tobasin_glo = undef_int WHERE ( tmp_real_g .LT. val_exp ) route_tobasin_glo = NINT(tmp_real_g) ENDWHERE ENDIF CALL scatter(route_tobasin_glo,route_tobasin_loc) route_tobasin=>route_tobasin_loc ! ! nbintobasin ! ALLOCATE (route_nbintobas_loc(nbpt,nbasmax), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for route_nbintobas_loc','','') ALLOCATE (route_nbintobas_glo(nbp_glo,nbasmax), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for route_nbintobas_glo','','') IF (is_root_prc) THEN var_name = 'routenbintobas' CALL ioconf_setatt('UNITS', '-') CALL ioconf_setatt('LONG_NAME','Number of basin into current one') CALL restget (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., tmp_real_g, "gather", nbp_glo, index_g) route_nbintobas_glo = undef_int WHERE ( tmp_real_g .LT. val_exp ) route_nbintobas_glo = NINT(tmp_real_g) ENDWHERE ENDIF CALL scatter(route_nbintobas_glo,route_nbintobas_loc) route_nbintobas=>route_nbintobas_loc ! ALLOCATE (global_basinid_loc(nbpt,nbasmax), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for global_basinid_loc','','') ALLOCATE (global_basinid_glo(nbp_glo,nbasmax), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for global_basinid_glo','','') IF (is_root_prc) THEN var_name = 'basinid' CALL ioconf_setatt('UNITS', '-') CALL ioconf_setatt('LONG_NAME','ID of basin') CALL restget (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., tmp_real_g, "gather", nbp_glo, index_g) global_basinid_glo = undef_int WHERE ( tmp_real_g .LT. val_exp ) global_basinid_glo = NINT(tmp_real_g) ENDWHERE ENDIF CALL scatter(global_basinid_glo,global_basinid_loc) global_basinid=>global_basinid_loc ! ALLOCATE (topo_resid_loc(nbpt,nbasmax), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for topo_resid_loc','','') ALLOCATE (topo_resid_glo(nbp_glo,nbasmax), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for topo_resid_glo','','') IF (is_root_prc) THEN var_name = 'topoindex' CALL ioconf_setatt('UNITS', 'm') CALL ioconf_setatt('LONG_NAME','Topographic index of the residence time') CALL restget (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., topo_resid_glo, "gather", nbp_glo, index_g) ENDIF CALL scatter(topo_resid_glo,topo_resid_loc) topo_resid=>topo_resid_loc ALLOCATE (fast_reservoir(nbpt,nbasmax), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for fast_reservoir','','') var_name = 'fastres' CALL ioconf_setatt_p('UNITS', 'Kg') CALL ioconf_setatt_p('LONG_NAME','Water in the fast reservoir') CALL restget_p (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., fast_reservoir, "gather", nbp_glo, index_g) CALL setvar_p (fast_reservoir, val_exp, 'NO_KEYWORD', zero) ALLOCATE (slow_reservoir(nbpt,nbasmax), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for slow_reservoir','','') var_name = 'slowres' CALL ioconf_setatt_p('UNITS', 'Kg') CALL ioconf_setatt_p('LONG_NAME','Water in the slow reservoir') CALL restget_p (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., slow_reservoir, "gather", nbp_glo, index_g) CALL setvar_p (slow_reservoir, val_exp, 'NO_KEYWORD', zero) ALLOCATE (stream_reservoir(nbpt,nbasmax), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for stream_reservoir','','') var_name = 'streamres' CALL ioconf_setatt_p('UNITS', 'Kg') CALL ioconf_setatt_p('LONG_NAME','Water in the stream reservoir') CALL restget_p (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., stream_reservoir, "gather", nbp_glo, index_g) CALL setvar_p (stream_reservoir, val_exp, 'NO_KEYWORD', zero) ALLOCATE (flood_reservoir(nbpt,nbasmax), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for flood_reservoir','','') var_name = 'floodres' CALL ioconf_setatt_p('UNITS', 'Kg') CALL ioconf_setatt_p('LONG_NAME','Water in the flood reservoir') CALL restget_p (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., flood_reservoir, "gather", nbp_glo, index_g) CALL setvar_p (flood_reservoir, val_exp, 'NO_KEYWORD', zero) ALLOCATE (flood_frac_bas(nbpt,nbasmax), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for flood_frac_bas','','') var_name = 'flood_frac_bas' CALL ioconf_setatt_p('UNITS', '-') CALL ioconf_setatt_p('LONG_NAME','Flooded fraction per basin') CALL restget_p (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., flood_frac_bas, "gather", nbp_glo, index_g) CALL setvar_p (flood_frac_bas, val_exp, 'NO_KEYWORD', zero) ALLOCATE (flood_height(nbpt), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for flood_height','','') var_name = 'floodh' CALL ioconf_setatt_p('UNITS', '-') CALL ioconf_setatt_p('LONG_NAME','') CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., flood_height, "gather", nbp_glo, index_g) CALL setvar_p (flood_height, val_exp, 'NO_KEYWORD', zero) ALLOCATE (pond_frac(nbpt), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for pond_frac','','') var_name = 'pond_frac' CALL ioconf_setatt_p('UNITS', '-') CALL ioconf_setatt_p('LONG_NAME','Pond fraction per grid box') CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., pond_frac, "gather", nbp_glo, index_g) CALL setvar_p (pond_frac, val_exp, 'NO_KEYWORD', zero) var_name = 'flood_frac' CALL ioconf_setatt_p('UNITS', '-') CALL ioconf_setatt_p('LONG_NAME','Flooded fraction per grid box') CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., flood_frac, "gather", nbp_glo, index_g) CALL setvar_p (flood_frac, val_exp, 'NO_KEYWORD', zero) var_name = 'flood_res' CALL ioconf_setatt_p('UNITS','mm') CALL ioconf_setatt_p('LONG_NAME','Flooded quantity (estimation)') CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., flood_res, "gather", nbp_glo, index_g) CALL setvar_p (flood_res, val_exp, 'NO_KEYWORD', zero) ! flood_res = zero ALLOCATE (lake_reservoir(nbpt), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for lake_reservoir','','') var_name = 'lakeres' CALL ioconf_setatt_p('UNITS', 'Kg') CALL ioconf_setatt_p('LONG_NAME','Water in the lake reservoir') CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., lake_reservoir, "gather", nbp_glo, index_g) CALL setvar_p (lake_reservoir, val_exp, 'NO_KEYWORD', zero) ALLOCATE (pond_reservoir(nbpt), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for pond_reservoir','','') var_name = 'pondres' CALL ioconf_setatt_p('UNITS', 'Kg') CALL ioconf_setatt_p('LONG_NAME','Water in the pond reservoir') CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., pond_reservoir, "gather", nbp_glo, index_g) CALL setvar_p (pond_reservoir, val_exp, 'NO_KEYWORD', zero) ! ! Map of irrigated areas ! IF ( do_irrigation ) THEN ALLOCATE (irrigated(nbpt), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for irrigated','','') var_name = 'irrigated' CALL ioconf_setatt_p('UNITS', 'm^2') CALL ioconf_setatt_p('LONG_NAME','Surface of irrigated area') CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., irrigated, "gather", nbp_glo, index_g) CALL setvar_p (irrigated, val_exp, 'NO_KEYWORD', undef_sechiba) ENDIF IF ( do_floodplains ) THEN ALLOCATE (floodplains(nbpt), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for floodplains','','') var_name = 'floodplains' CALL ioconf_setatt_p('UNITS', 'm^2') CALL ioconf_setatt_p('LONG_NAME','Surface which can be flooded') CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., floodplains, "gather", nbp_glo, index_g) CALL setvar_p (floodplains, val_exp, 'NO_KEYWORD', undef_sechiba) ENDIF IF ( doswamps ) THEN ALLOCATE (swamp(nbpt), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for swamp','','') var_name = 'swamp' CALL ioconf_setatt_p('UNITS', 'm^2') CALL ioconf_setatt_p('LONG_NAME','Surface which can become swamp') CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., swamp, "gather", nbp_glo, index_g) CALL setvar_p (swamp, val_exp, 'NO_KEYWORD', undef_sechiba) ENDIF ! ! Put into the restart file the fluxes so that they can be regenerated at restart. ! ALLOCATE (lakeinflow_mean(nbpt), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for lakeinflow_mean','','') var_name = 'lakeinflow' CALL ioconf_setatt_p('UNITS', 'Kg/dt') CALL ioconf_setatt_p('LONG_NAME','Lake inflow') CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., lakeinflow_mean, "gather", nbp_glo, index_g) CALL setvar_p (lakeinflow_mean, val_exp, 'NO_KEYWORD', zero) ALLOCATE (returnflow_mean(nbpt), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for returnflow_mean','','') var_name = 'returnflow' CALL ioconf_setatt_p('UNITS', 'Kg/m^2/dt') CALL ioconf_setatt_p('LONG_NAME','Deep return flux') CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., returnflow_mean, "gather", nbp_glo, index_g) CALL setvar_p (returnflow_mean, val_exp, 'NO_KEYWORD', zero) returnflow(:) = returnflow_mean(:) ALLOCATE (reinfiltration_mean(nbpt), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for reinfiltration_mean','','') var_name = 'reinfiltration' CALL ioconf_setatt_p('UNITS', 'Kg/m^2/dt') CALL ioconf_setatt_p('LONG_NAME','Top return flux') CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., reinfiltration_mean, "gather", nbp_glo, index_g) CALL setvar_p (reinfiltration_mean, val_exp, 'NO_KEYWORD', zero) reinfiltration(:) = reinfiltration_mean(:) ALLOCATE (irrigation_mean(nbpt), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for irrigation_mean','','') ALLOCATE (irrig_netereq(nbpt), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for irrig_netereq','','') irrig_netereq(:) = zero IF ( do_irrigation ) THEN var_name = 'irrigation' CALL ioconf_setatt_p('UNITS', 'Kg/dt') CALL ioconf_setatt_p('LONG_NAME','Artificial irrigation flux') CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., irrigation_mean, "gather", nbp_glo, index_g) CALL setvar_p (irrigation_mean, val_exp, 'NO_KEYWORD', zero) ELSE irrigation_mean(:) = zero ENDIF irrigation(:) = irrigation_mean(:) ALLOCATE (riverflow_mean(nbpt), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for riverflow_mean','','') var_name = 'riverflow' CALL ioconf_setatt_p('UNITS', 'Kg/dt') CALL ioconf_setatt_p('LONG_NAME','River flux into the sea') CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., riverflow_mean, "gather", nbp_glo, index_g) CALL setvar_p (riverflow_mean, val_exp, 'NO_KEYWORD', zero) riverflow(:) = riverflow_mean(:) ALLOCATE (coastalflow_mean(nbpt), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for coastalflow_mean','','') var_name = 'coastalflow' CALL ioconf_setatt_p('UNITS', 'Kg/dt') CALL ioconf_setatt_p('LONG_NAME','Diffuse flux into the sea') CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., coastalflow_mean, "gather", nbp_glo, index_g) CALL setvar_p (coastalflow_mean, val_exp, 'NO_KEYWORD', zero) coastalflow(:) = coastalflow_mean(:) ! Locate it at the 2m level ipn = MINLOC(ABS(diaglev-2)) floodtemp_lev = ipn(1) ALLOCATE (floodtemp(nbpt), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for floodtemp','','') floodtemp(:) = stempdiag(:,floodtemp_lev) ALLOCATE(hydrographs(nbpt), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for hydrographs','','') var_name = 'hydrographs' CALL ioconf_setatt_p('UNITS', 'm^3/dt') CALL ioconf_setatt_p('LONG_NAME','Hydrograph at outlow of grid') CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., hydrographs, "gather", nbp_glo, index_g) CALL setvar_p (hydrographs, val_exp, 'NO_KEYWORD', zero) ALLOCATE(slowflow_diag(nbpt), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for slowflow_diag','','') var_name = 'slowflow_diag' CALL ioconf_setatt_p('UNITS', 'm^3/dt') CALL ioconf_setatt_p('LONG_NAME','Slowflow hydrograph at outlow of grid') CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE.,slowflow_diag, "gather", nbp_glo, index_g) CALL setvar_p (slowflow_diag, val_exp, 'NO_KEYWORD', zero) ! ! The diagnostic variables, they are initialized from the above restart variables. ! ALLOCATE(fast_diag(nbpt), slow_diag(nbpt), stream_diag(nbpt), flood_diag(nbpt), & & pond_diag(nbpt), lake_diag(nbpt), delsurfstor(nbpt), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for fast_diag,..','','') fast_diag(:) = zero slow_diag(:) = zero stream_diag(:) = zero flood_diag(:) = zero pond_diag(:) = zero lake_diag(:) = zero delsurfstor(:) = zero DO ig=1,nbpt totarea = zero DO ib=1,nbasmax totarea = totarea + routing_area(ig,ib) fast_diag(ig) = fast_diag(ig) + fast_reservoir(ig,ib) slow_diag(ig) = slow_diag(ig) + slow_reservoir(ig,ib) stream_diag(ig) = stream_diag(ig) + stream_reservoir(ig,ib) flood_diag(ig) = flood_diag(ig) + flood_reservoir(ig,ib) ENDDO ! fast_diag(ig) = fast_diag(ig)/totarea slow_diag(ig) = slow_diag(ig)/totarea stream_diag(ig) = stream_diag(ig)/totarea flood_diag(ig) = flood_diag(ig)/totarea ! ! This is the volume of the lake scaled to the entire grid. ! It would be batter to scale it to the size of the lake ! but this information is not yet available. ! lake_diag(ig) = lake_reservoir(ig)/totarea ! ENDDO ! ! Get from the restart the fluxes we accumulated. ! ALLOCATE (floodout_mean(nbpt), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for floodout_mean','','') var_name = 'floodout_route' CALL ioconf_setatt_p('UNITS', 'Kg') CALL ioconf_setatt_p('LONG_NAME','Accumulated flow out of floodplains for routing') CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., floodout_mean, "gather", nbp_glo, index_g) CALL setvar_p (floodout_mean, val_exp, 'NO_KEYWORD', zero) ALLOCATE (runoff_mean(nbpt), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for runoff_mean','','') var_name = 'runoff_route' CALL ioconf_setatt_p('UNITS', 'Kg') CALL ioconf_setatt_p('LONG_NAME','Accumulated runoff for routing') CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., runoff_mean, "gather", nbp_glo, index_g) CALL setvar_p (runoff_mean, val_exp, 'NO_KEYWORD', zero) ALLOCATE(drainage_mean(nbpt), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for drainage_mean','','') var_name = 'drainage_route' CALL ioconf_setatt_p('UNITS', 'Kg') CALL ioconf_setatt_p('LONG_NAME','Accumulated drainage for routing') CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., drainage_mean, "gather", nbp_glo, index_g) CALL setvar_p (drainage_mean, val_exp, 'NO_KEYWORD', zero) ALLOCATE(transpot_mean(nbpt), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for transpot_mean','','') var_name = 'transpot_route' CALL ioconf_setatt_p('UNITS', 'Kg/m^2') CALL ioconf_setatt_p('LONG_NAME','Accumulated potential transpiration for routing/irrigation') CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., transpot_mean, "gather", nbp_glo, index_g) CALL setvar_p (transpot_mean, val_exp, 'NO_KEYWORD', zero) ALLOCATE(precip_mean(nbpt), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for precip_mean','','') var_name = 'precip_route' CALL ioconf_setatt_p('UNITS', 'Kg/m^2') CALL ioconf_setatt_p('LONG_NAME','Accumulated rain precipitation for irrigation') CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., precip_mean, "gather", nbp_glo, index_g) CALL setvar_p (precip_mean, val_exp, 'NO_KEYWORD', zero) ALLOCATE(humrel_mean(nbpt), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for humrel_mean','','') var_name = 'humrel_route' CALL ioconf_setatt_p('UNITS', '-') CALL ioconf_setatt_p('LONG_NAME','Mean humrel for irrigation') CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., humrel_mean, "gather", nbp_glo, index_g) CALL setvar_p (humrel_mean, val_exp, 'NO_KEYWORD', un) ALLOCATE(k_litt_mean(nbpt), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for k_litt_mean','','') var_name = 'k_litt_route' CALL ioconf_setatt_p('UNITS', '-') CALL ioconf_setatt_p('LONG_NAME','Mean cond. for litter') CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., k_litt_mean, "gather", nbp_glo, index_g) CALL setvar_p (k_litt_mean, val_exp, 'NO_KEYWORD', zero) ALLOCATE(totnobio_mean(nbpt), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for totnobio_mean','','') var_name = 'totnobio_route' CALL ioconf_setatt_p('UNITS', '-') CALL ioconf_setatt_p('LONG_NAME','Last Total fraction of no bio for irrigation') CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., totnobio_mean, "gather", nbp_glo, index_g) CALL setvar_p (totnobio_mean, val_exp, 'NO_KEYWORD', zero) ALLOCATE(vegtot_mean(nbpt), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for vegtot_mean','','') var_name = 'vegtot_route' CALL ioconf_setatt_p('UNITS', '-') CALL ioconf_setatt_p('LONG_NAME','Last Total fraction of vegetation') CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., vegtot_mean, "gather", nbp_glo, index_g) CALL setvar_p (vegtot_mean, val_exp, 'NO_KEYWORD', un) ! ! DEALLOCATE(tmp_real_g) ! ! Allocate diagnostic variables ! ALLOCATE(hydrodiag_loc(nbpt,nbasmax),hydrodiag_glo(nbp_glo,nbasmax),stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for hydrodiag_glo','','') hydrodiag=>hydrodiag_loc ALLOCATE(hydroupbasin_loc(nbpt),hydroupbasin_glo(nbp_glo), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for hydroupbasin_glo','','') hydroupbasin=>hydroupbasin_loc END SUBROUTINE routing_simple_init_1 !! ================================================================================================================================ !! SUBROUTINE : routing_simple_clear !! !>\BRIEF This subroutine deallocates the block memory previously allocated. !! !! DESCRIPTION: This subroutine deallocates the block memory previously allocated. !! !! RECENT CHANGE(S): None !! !! MAIN OUTPUT VARIABLE(S): !! !! REFERENCES : None !! !! FLOWCHART :None !! \n !_ ================================================================================================================================ SUBROUTINE routing_simple_clear() IF (is_omp_root) THEN IF (ALLOCATED(topoind_r)) DEALLOCATE(topoind_r) IF (ALLOCATED(route_flow_rp1)) DEALLOCATE(route_flow_rp1) IF (ALLOCATED(fast_reservoir_r)) DEALLOCATE(fast_reservoir_r) IF (ALLOCATED(slow_reservoir_r)) DEALLOCATE(slow_reservoir_r) IF (ALLOCATED(is_lakeinflow_r)) DEALLOCATE(is_lakeinflow_r) IF (ALLOCATED(is_coastalflow_r)) DEALLOCATE(is_coastalflow_r) IF (ALLOCATED(is_riverflow_r)) DEALLOCATE(is_riverflow_r) ENDIF IF (ALLOCATED(routing_area_loc)) DEALLOCATE(routing_area_loc) IF (ALLOCATED(route_togrid_loc)) DEALLOCATE(route_togrid_loc) IF (ALLOCATED(route_tobasin_loc)) DEALLOCATE(route_tobasin_loc) IF (ALLOCATED(route_nbintobas_loc)) DEALLOCATE(route_nbintobas_loc) IF (ALLOCATED(global_basinid_loc)) DEALLOCATE(global_basinid_loc) IF (ALLOCATED(topo_resid_loc)) DEALLOCATE(topo_resid_loc) IF (ALLOCATED(routing_area_glo)) DEALLOCATE(routing_area_glo) IF (ALLOCATED(route_togrid_glo)) DEALLOCATE(route_togrid_glo) IF (ALLOCATED(route_tobasin_glo)) DEALLOCATE(route_tobasin_glo) IF (ALLOCATED(route_nbintobas_glo)) DEALLOCATE(route_nbintobas_glo) IF (ALLOCATED(global_basinid_glo)) DEALLOCATE(global_basinid_glo) IF (ALLOCATED(topo_resid_glo)) DEALLOCATE(topo_resid_glo) IF (ALLOCATED(fast_reservoir)) DEALLOCATE(fast_reservoir) IF (ALLOCATED(slow_reservoir)) DEALLOCATE(slow_reservoir) IF (ALLOCATED(stream_reservoir)) DEALLOCATE(stream_reservoir) IF (ALLOCATED(flood_reservoir)) DEALLOCATE(flood_reservoir) IF (ALLOCATED(flood_frac_bas)) DEALLOCATE(flood_frac_bas) IF (ALLOCATED(flood_height)) DEALLOCATE(flood_height) IF (ALLOCATED(pond_frac)) DEALLOCATE(pond_frac) IF (ALLOCATED(lake_reservoir)) DEALLOCATE(lake_reservoir) IF (ALLOCATED(pond_reservoir)) DEALLOCATE(pond_reservoir) IF (ALLOCATED(returnflow_mean)) DEALLOCATE(returnflow_mean) IF (ALLOCATED(reinfiltration_mean)) DEALLOCATE(reinfiltration_mean) IF (ALLOCATED(riverflow_mean)) DEALLOCATE(riverflow_mean) IF (ALLOCATED(coastalflow_mean)) DEALLOCATE(coastalflow_mean) IF (ALLOCATED(lakeinflow_mean)) DEALLOCATE(lakeinflow_mean) IF (ALLOCATED(runoff_mean)) DEALLOCATE(runoff_mean) IF (ALLOCATED(floodout_mean)) DEALLOCATE(floodout_mean) IF (ALLOCATED(drainage_mean)) DEALLOCATE(drainage_mean) IF (ALLOCATED(transpot_mean)) DEALLOCATE(transpot_mean) IF (ALLOCATED(precip_mean)) DEALLOCATE(precip_mean) IF (ALLOCATED(humrel_mean)) DEALLOCATE(humrel_mean) IF (ALLOCATED(k_litt_mean)) DEALLOCATE(k_litt_mean) IF (ALLOCATED(totnobio_mean)) DEALLOCATE(totnobio_mean) IF (ALLOCATED(vegtot_mean)) DEALLOCATE(vegtot_mean) IF (ALLOCATED(floodtemp)) DEALLOCATE(floodtemp) IF (ALLOCATED(hydrodiag_loc)) DEALLOCATE(hydrodiag_loc) IF (ALLOCATED(hydrodiag_glo)) DEALLOCATE(hydrodiag_glo) IF (ALLOCATED(hydroupbasin_loc)) DEALLOCATE(hydroupbasin_loc) IF (ALLOCATED(hydroupbasin_glo)) DEALLOCATE(hydroupbasin_glo) IF (ALLOCATED(hydrographs)) DEALLOCATE(hydrographs) IF (ALLOCATED(slowflow_diag)) DEALLOCATE(slowflow_diag) IF (ALLOCATED(irrigation_mean)) DEALLOCATE(irrigation_mean) IF (ALLOCATED(irrigated)) DEALLOCATE(irrigated) IF (ALLOCATED(floodplains)) DEALLOCATE(floodplains) IF (ALLOCATED(swamp)) DEALLOCATE(swamp) IF (ALLOCATED(fast_diag)) DEALLOCATE(fast_diag) IF (ALLOCATED(slow_diag)) DEALLOCATE(slow_diag) IF (ALLOCATED(stream_diag)) DEALLOCATE(stream_diag) IF (ALLOCATED(flood_diag)) DEALLOCATE(flood_diag) IF (ALLOCATED(pond_diag)) DEALLOCATE(pond_diag) IF (ALLOCATED(lake_diag)) DEALLOCATE(lake_diag) IF (ALLOCATED(delsurfstor)) DEALLOCATE(delsurfstor) ! END SUBROUTINE routing_simple_clear ! !! ================================================================================================================================ !! SUBROUTINE : routing_simple_lake !! !>\BRIEF : This subroutine stores water in lakes so that it does not cycle through the runoff. !! For the moment it only works for endoheric lakes but I can be extended in the future. !! !! DESCRIPTION (definitions, functional, design, flags): The return flow to the soil moisture reservoir !! is based on a maximum lake evaporation rate (maxevap_lake). \n !! !! RECENT CHANGE(S): None !! !! MAIN OUTPUT VARIABLE(S): !! !! REFERENCES : None !! !! FLOWCHART :None !! \n !_ ================================================================================================================================ SUBROUTINE routing_simple_lake(nbpt, dt_routing, lakeinflow, humrel, contfrac, return_lakes) USE grid, ONLY : area IMPLICIT NONE !! 0 Variable and parameter description !! 0.1 Input variables INTEGER(i_std), INTENT(in) :: nbpt !! Domain size (unitless) REAL(r_std), INTENT (in) :: dt_routing !! Routing time step (s) REAL(r_std), INTENT(in) :: lakeinflow(nbpt) !! Water inflow to the lakes (kg/dt) REAL(r_std), INTENT(in) :: humrel(nbpt) !! Soil moisture stress, root extraction potential (unitless) REAL(r_std), INTENT(in) :: contfrac(nbpt) !! Fraction of land in each grid box (unitless;0-1) !! 0.2 Output variables REAL(r_std), INTENT(out) :: return_lakes(nbpt) !! Water from lakes flowing back into soil moisture (kg/m^2/dt) !! 0.3 Local variables INTEGER(i_std) :: ig !! Indices (unitless) REAL(r_std) :: refill !! REAL(r_std) :: total_area !! Sum of all the surfaces of the basins (m^2) !_ ================================================================================================================================ DO ig=1,nbpt ! total_area = area(ig)*contfrac(ig) ! lake_reservoir(ig) = lake_reservoir(ig) + lakeinflow(ig) !uptake in Kg/dt refill = MAX(zero, maxevap_lake * (un - humrel(ig)) * dt_routing * total_area) return_lakes(ig) = MIN(refill, lake_reservoir(ig)) lake_reservoir(ig) = lake_reservoir(ig) - return_lakes(ig) !Return in Kg/m^2/dt IF ( doswamps ) THEN return_lakes(ig) = return_lakes(ig)/total_area ELSE return_lakes(ig) = zero ENDIF ! ! This is the volume of the lake scaled to the entire grid. ! It would be batter to scale it to the size of the lake ! but this information is not yet available. lake_diag(ig) = lake_reservoir(ig)/total_area ! ENDDO END SUBROUTINE routing_simple_lake END MODULE routing_simple