! ================================================================================================================================ ! MODULE : routing_highres ! ! 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: None !! !! RECENT CHANGE(S): Now works together with the routing pre-processor : https://gitlab.in2p3.fr/ipsl/lmd/intro/routingpp !! !! REFERENCE(S) : !! !! SVN : !! $HeadURL: svn://forge.ipsl.jussieu.fr/orchidee/branches/ORCHIDEE-ROUTING/ORCHIDEE/src_sechiba/routing.f90 $ !! $Date: 2022-03-24 11:25:05 +0100 (Do, 24 Mär 2022) $ !! $Revision: 7545 $ !! \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_highres USE ioipsl USE xios_orchidee USE ioipsl_para USE constantes USE constantes_var USE constantes_soil USE pft_parameters USE sechiba_io_p USE interpol_help USE grid USE mod_orchidee_para USE haversine IMPLICIT NONE PRIVATE PUBLIC :: routing_highres_main, routing_highres_initialize, routing_highres_finalize, routing_highres_clear, routing_highres_xios_initialize INTERFACE routing_hr_landgather MODULE PROCEDURE routing_hr_landgather_i1, routing_hr_landgather_i2, routing_hr_landgather_r END INTERFACE routing_hr_landgather INTEGER(i_std),PARAMETER :: WaterCp=1000.*4.1813 !! water heat capacity in J/Kg/K !! PARAMETERS INTEGER(i_std), SAVE :: nbasmax=-1 !! The maximum number of basins we wish to have per grid box (truncation of the model) (unitless) INTEGER(i_std), SAVE :: nbasmon = 4 !! Number of basins to be monitored INTEGER(i_std), SAVE :: inflows=-1 !! The maximum number of inflows (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 = -1. !! Property of the fast reservoir, (s/km) !$OMP THREADPRIVATE(fast_tcst) REAL(r_std), SAVE :: slow_tcst = -1. !! Property of the slow reservoir, (s/km) !$OMP THREADPRIVATE(slow_tcst) REAL(r_std), SAVE :: stream_tcst = -1. !! Property of the stream reservoir, (s/km) !$OMP THREADPRIVATE(stream_tcst) REAL(r_std), SAVE :: flood_tcst = -1. !! Property of the floodplains reservoir, (s/km) !$OMP THREADPRIVATE(flood_tcst) REAL(r_std), SAVE :: swamp_cst = -1. !! Fraction of the river transport that flows to the swamps (unitless;0-1) !$OMP THREADPRIVATE(swamp_cst) REAL(r_std), SAVE :: lim_floodcri = -1. !! Minimal orog diff between two consecutive floodplains htu (m) !$OMP THREADPRIVATE(lim_floodcri) ! ! Relation between volume and fraction of floodplains ! 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 :: rfloodmax = 0.5 !! Maximal discharge reducer when there are floodplains !$OMP THREADPRIVATE(rfloodmax) REAL(r_std), SAVE :: overflow_tcst = 5 !! Maximal discharge reducer when there are floodplains !$OMP THREADPRIVATE(overflow_tcst) INTEGER(i_std), SAVE :: overflow_repetition = 1 !! Number of repetition of overflow for each routing step !$OMP THREADPRIVATE(overflow_repetition) ! Soil temperature depth to be used to estimate runoff and drainage temperatures ! REAL(r_std), PARAMETER, DIMENSION(2) :: runofftempdepth = (/ 0.0, 0.3 /) !! Layer which will determine the temperature of runoff REAL(r_std), PARAMETER, DIMENSION(2) :: drainagetempdepth = (/ 3.0, 90.0 /) !! Layer which will determine the temperature of runoff ! ! ! 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 :: ntemp_layer = 4 !! Number of layers to be taken to determine the ground water temperature. !$OMP THREADPRIVATE(ntemp_layer) 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 :: dofloodoverflow = .FALSE. !! Logical to choose if floodplains overflow is activated or not (true/false) !$OMP THREADPRIVATE(dofloodoverflow) 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) REAL(r_std), SAVE :: conduct_factor = 1. !! Adjustment factor for floodplains reinfiltration !$OMP THREADPRIVATE(conduct_factor) ! ! The variables describing the basins and their routing, need to be in the restart file. ! INTEGER(i_std), SAVE :: num_largest = 200 !! 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) ! CHARACTER(LEN=80),SAVE :: graphfilename="routing_graph.nc" !$OMP THREADPRIVATE(graphfilename) REAL(r_std), SAVE :: undef_graphfile !$OMP THREADPRIVATE(undef_graphfile) REAL(r_std), SAVE :: graphfile_version = 0.0 !$OMP THREADPRIVATE(graphfile_version) REAL(r_std), SAVE :: maxtimestep = 1800.0 !! A reasonalble maximum time step. Actual value to be read from graphfile. !$OMP THREADPRIVATE(maxtimestep) 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) REAL(r_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: stream_resid_loc !! Topographic index of the retention time (m) !$OMP THREADPRIVATE(stream_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) INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: HTUdiag_loc !! Variable to diagnose the hydrographs !$OMP THREADPRIVATE(HTUdiag_loc) INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: HTUdiag_glo !! Variable to diagnose the hydrographs !$OMP THREADPRIVATE(HTUdiag_glo) LOGICAL, SAVE :: MonitoringinGraph=.FALSE. LOGICAL, SAVE :: ReadGraph=.FALSE. LOGICAL, SAVE :: ReadMonitoring=.FALSE. REAL(r_std), SAVE :: stream_maxresid !$OMP THREADPRIVATE(stream_maxresid) ! ! 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) REAL(r_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: stream_resid_glo !! Topographic index of the retention time (m) !$OMP THREADPRIVATE(stream_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, 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) REAL(r_std), SAVE, POINTER, DIMENSION(:,:) :: stream_resid !! Topographic index of the retention time (m) !$OMP THREADPRIVATE(stream_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, ALLOCATABLE, DIMENSION(:) :: irrigated !! Area equipped for irrigation in each grid box (m^2) !$OMP THREADPRIVATE(irrigated) REAL(r_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: floodplains_glo !! Maximal surface which can be inundated in each grid box (m^2) !$OMP THREADPRIVATE(floodplains_glo) REAL(r_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: floodplains_loc !! Maximal surface which can be inundated in each grid box (m^2) !$OMP THREADPRIVATE(floodplains_loc) REAL(r_std), SAVE, POINTER, DIMENSION(:,:) :: floodplains !! Maximal surface which can be inundated in each grid box (m^2) !$OMP THREADPRIVATE(floodplains) REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: floodmap !! Floodplains Fraction for each grid point. !$OMP THREADPRIVATE(floodmap) !!! REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: tempdiag_mean !! Averaged soil temperatures !$OMP THREADPRIVATE(tempdiag_mean) ! ! FLOOD OVERFLOW REAL(r_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: orog_min_glo !! !$OMP THREADPRIVATE(orog_min_glo)! REAL(r_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: orog_min_loc !! !$OMP THREADPRIVATE(orog_min_loc) REAL(r_std), SAVE, POINTER, DIMENSION(:,:) :: orog_min !! !$OMP THREADPRIVATE(orog_min) ! INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: route_innum_glo !! !$OMP THREADPRIVATE(route_innum_glo) INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: route_innum_loc !! !$OMP THREADPRIVATE(route_innum_loc) INTEGER(i_std), SAVE, POINTER, DIMENSION(:,:) :: route_innum !! !$OMP THREADPRIVATE(route_innum) ! INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:,:) :: route_ingrid_glo !! !$OMP THREADPRIVATE(route_ingrid_glo) INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:,:) :: route_ingrid_loc !! !$OMP THREADPRIVATE(route_ingrid_loc) INTEGER(i_std), SAVE, POINTER, DIMENSION(:,:,:) :: route_ingrid !! !$OMP THREADPRIVATE(route_ingrid) ! INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:,:) :: route_inbasin_glo !! !$OMP THREADPRIVATE(route_inbasin_glo) INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:,:) :: route_inbasin_loc !! !$OMP THREADPRIVATE(route_inbasin_loc) INTEGER(i_std), SAVE, POINTER, DIMENSION(:,:,:) :: route_inbasin !! !$OMP THREADPRIVATE(route_inbasin) ! !!! REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: swamp !! Maximal surface of swamps in each grid box (m^2) !$OMP THREADPRIVATE(swamp) REAL(r_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: fp_beta_glo !! Parameter to fix the shape of the floodplain (>1 for convex edges, <1 for concave edges) (unitless) !$OMP THREADPRIVATE(fp_beta_glo) REAL(r_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: fp_beta_loc !! Parameter to fix the shape of the floodplain (>1 for convex edges, <1 for concave edges) (unitless) !$OMP THREADPRIVATE(fp_beta_loc) REAL(r_std), SAVE, POINTER, DIMENSION(:,:) :: fp_beta !! Parameter to fix the shape of the floodplain (>1 for convex edges, <1 for concave edges) (unitless) !$OMP THREADPRIVATE(fp_beta) REAL(r_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: floodcri_glo !! Potential height for which all the basin is a pond (mm) !$OMP THREADPRIVATE(floodcri_glo) REAL(r_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: floodcri_loc !! Potential height for which all the basin is a pond (mm) !$OMP THREADPRIVATE(floodcri_loc) REAL(r_std), SAVE, POINTER, DIMENSION(:,:) :: floodcri !! Potential height for which all the basin is a pond (mm) !$OMP THREADPRIVATE(floodcri) ! ! The reservoirs, also to be put into 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) ! ! Reservoir temperatures ! REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: fast_temp !! Water temperature in the fast reservoir (K) !$OMP THREADPRIVATE(fast_temp) REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: slow_temp !! Water temperature in the slow reservoir (K) !$OMP THREADPRIVATE(slow_temp) REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: stream_temp !! Water temperature in the stream reservoir (K) !$OMP THREADPRIVATE(stream_temp) REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: streamlimit !! !$OMP THREADPRIVATE(streamlimit) ! ! 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 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) REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: hydrotemp !! Temperature of the largest river (in the HTUdiag sense) in the grid (K) !$OMP THREADPRIVATE(hydrotemp) REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: HTUhgmon !! Hydrographs to be monitored on specific HTUs (kg/dt) !$OMP THREADPRIVATE(HTUhgmon) REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: HTUhgmon_glo !! Hydrographs to be monitored on specific HTUs (kg/dt) !$OMP THREADPRIVATE(HTUhgmon_glo) REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: HTUtempmon !! Temperature to be monitored on specific HTUs (K) !$OMP THREADPRIVATE(HTUtempmon) REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: HTUtempmon_glo !! Temperature to be monitored on specific HTUs (K) !$OMP THREADPRIVATE(HTUtempmon_glo) ! ! 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(:) :: mask_coast !! Mask with coastal gridcells on local grid(1/0) !$OMP THREADPRIVATE(mask_coast) REAL(r_std), SAVE :: max_lake_reservoir !! Maximum limit of water in lake_reservoir [kg/m2] !$OMP THREADPRIVATE(max_lake_reservoir) INTEGER(i_std), SAVE :: nb_coast_gridcells !! Number of gridcells which can receive coastalflow !$OMP THREADPRIVATE(nb_coast_gridcells) CONTAINS !! ============================================================================================================================= !! SUBROUTINE: routing_highres_initialize !! !>\BRIEF Initialize the routing module !! !! DESCRIPTION: Initialize the routing 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_highres_initialize( kjit, nbpt, index, & rest_id, hist_id, hist2_id, lalo, & neighbours, resolution, contfrac, tempdiag, & returnflow, reinfiltration, irrigation, riverflow, & coastalflow, flood_frac, flood_res ) IMPLICIT NONE !! 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,NbNeighb) !! 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) :: tempdiag(nbpt,ngrnd) !! 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 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 REAL(r_std), DIMENSION(nbp_glo) :: mask_coast_glo !! Mask with coastal gridcells on global grid (1/0) 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) INTEGER :: ig, ib, rtg, rtb !! Index REAL(r_std) :: stream_tcst_orig INTEGER :: ier !! Error handeling !_ ================================================================================================================================ ! ! 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_hr_init (kjit, nbpt, index, returnflow, reinfiltration, irrigation, & riverflow, coastalflow, flood_frac, flood_res, tempdiag, rest_id) routing_area => routing_area_loc floodplains => floodplains_loc topo_resid => topo_resid_loc stream_resid => stream_resid_loc route_togrid => route_togrid_loc route_tobasin => route_tobasin_loc global_basinid => global_basinid_loc hydrodiag => hydrodiag_loc fp_beta => fp_beta_loc floodcri => floodcri_loc ! route_innum => route_innum_loc route_ingrid => route_ingrid_loc route_inbasin => route_inbasin_loc orog_min => orog_min_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 ( ReadGraph .OR. ReadMonitoring) THEN CALL routing_hr_basins_p(nbpt, lalo, neighbours, resolution, contfrac) ENDIF ! Keep the information so we can check the time step. stream_tcst_orig = stream_tcst ! IF (stream_tcst .LE. 0 .OR. fast_tcst .LE. 0 .OR. slow_tcst .LE. 0 .OR. flood_tcst .LE. 0 ) THEN CALL ipslerr(3,'routing_highres_initialize',' The time constants of the routing reservoirs were not initialized. ', & 'Please check if they are present in the HTU graph file', ' ') ELSE ! !> The time constants for the various reservoirs should be read from the graph file !> produced by routingpp (https://gitlab.in2p3.fr/ipsl/lmd/intro/routingpp). They are !> also saved in the restart file so that we do not need to read the graph file at each restart. !> But once they are set in the model the user can changed them through the run.def. !> This is a useful option to test values but should not be an operational solution. The !> correct value should be given to the model through the graph file. !> The getin_p operation cannot be done earlier as in routing_hr_init above these constant !> might not yet be known. ! !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 (s/km) of the slow reservoir !Config in order to get better river flows for !Config particular regions. !Config Units = [days] ! 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 (s/km) 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 (s/km) 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 (s/km) 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 = LIM_FLOODCRI !Config Desc = Difference of orography between floodplains HTUs. !Config If = RIVER_ROUTING !Config Def = 0.3 !Config Help = This parameters allows the user to fix the !Config minimal difference of orography between two consecutive !Config floodplains HTU. !Config Units = [meter] CALL getin_p('LIM_FLOODCRI', lim_floodcri) ! ENDIF ! ! Verify that the time step is compatible with the graph file. ! If the user has changed the time constant of the stream reservoir then ! the maximum time step needs to be adjusted. ! IF ( stream_tcst_orig == 0 ) THEN WRITE(*,*) "routing_highres_initialize : Update stream_tcst ", stream_tcst_orig, stream_tcst stream_tcst_orig = stream_tcst ENDIF IF ( dt_routing > maxtimestep/stream_tcst_orig*stream_tcst ) THEN WRITE(*,*) "routing_highres_initialize : Chosen time step : ", dt_routing WRITE(*,*) "routing_highres_initialize : Recommended time step : ", maxtimestep/stream_tcst_orig*stream_tcst CALL ipslerr_p(2,'routing_highres_initialize','The chosen time step is larger than the value recommended','in the graph file.','') ENDIF ! ! ! IF (dofloodoverflow) THEN CALL routing_hr_inflows(nbp_glo, nbasmax, inflows, floodplains_glo,route_innum_glo,route_ingrid_glo,route_inbasin_glo) END IF !! Create a mask containing all possible coastal gridcells and count total number of coastal gridcells IF (is_root_prc) THEN mask_coast_glo(:)=0 DO ib=1,nbasmax DO ig=1,nbp_glo rtg = route_togrid_glo(ig,ib) rtb = route_tobasin_glo(ig,ib) ! Coastal gridcells are stored in nbasmax+2 IF ( rtb == nbasmax+2) THEN mask_coast_glo(rtg) = 1 END IF END DO END DO nb_coast_gridcells=SUM(mask_coast_glo) IF (printlev>=3) WRITE(numout,*) 'Number of coastal gridcells = ', nb_coast_gridcells IF (nb_coast_gridcells == 0)THEN CALL ipslerr(3,'routing_highres_initialize',& 'Number of coastal gridcells is zero for routing. ', & 'If this is a global run, this is an error.',& 'If this is a regional run, please check to make sure your region includes a full basin or turn routing off.') ENDIF ENDIF CALL bcast(nb_coast_gridcells) ALLOCATE(mask_coast(nbpt), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_inititalize','Pb in allocate for mask_coast','','') CALL scatter(mask_coast_glo, mask_coast) ! ! 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_hr_irrigmap(nbpt, index, lalo, neighbours, resolution, & ! contfrac, init_irrig, irrigated, init_flood, floodplains, init_swamp, swamp, hist_id, hist2_id) !ENDIF IF (printlev >= 5) WRITE(numout,*) 'End of routing_highres_initialize' END SUBROUTINE routing_highres_initialize !! ============================================================================================================================= !! SUBROUTINE: routing_highres_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_highres_xios_initialize USE xios_orchidee IMPLICIT NONE INTEGER(i_std) ::ib ! ! If the routing_graph file is available we will extract the information in the dimensions ! and parameters. ! !Config Key = ROUTING_FILE !Config Desc = Name of file which contains the routing information graph on the model grid !Config If = RIVER_ROUTING !Config Def = routing.nc !Config Help = The file provided here should allows to route the water from one HTU !Config to another. The RoutingPP code needs to be used in order to generate !Config the routing graph for the model grid. !Config More details on : https://gitlab.in2p3.fr/ipsl/lmd/intro/routingpp !Config Units = [FILE] ! graphfilename = 'routing_graph.nc' CALL getin('ROUTING_FILE',graphfilename) CALL routing_hr_graphinfo(graphfilename, nbasmax, inflows, nbasmon, undef_graphfile, stream_tcst, fast_tcst, slow_tcst, & & flood_tcst, swamp_cst, lim_floodcri) CALL xios_orchidee_addaxis("nbhtu", nbasmax, (/(REAL(ib,r_std),ib=1,nbasmax)/)) CALL xios_orchidee_addaxis("nbasmon", nbasmon, (/(REAL(ib,r_std),ib=1,nbasmon)/)) END SUBROUTINE routing_highres_xios_initialize !! ================================================================================================================================ !! SUBROUTINE : routing_highres_main !! !>\BRIEF This module routes the water over the continents (runoff and !! drainage produced by the 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_highres_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, & & tempdiag, reinf_slope, returnflow, reinfiltration, irrigation, riverflow, coastalflow, rest_id, hist_id, hist2_id) IMPLICIT NONE !! 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,NbNeighb) !! 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) :: tempdiag(nbpt,ngrnd) !! 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 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 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), DIMENSION(nbpt) :: fast_diag_old !! Reservoir in the beginning of the time step REAL(r_std), DIMENSION(nbpt) :: slow_diag_old !! Reservoir in the beginning of the time step REAL(r_std), DIMENSION(nbpt) :: stream_diag_old !! Reservoir in the beginning of the time step REAL(r_std), DIMENSION(nbpt) :: lake_diag_old !! Reservoir in the beginning of the time step REAL(r_std), DIMENSION(nbpt) :: pond_diag_old !! Reservoir in the beginning of the time step REAL(r_std), DIMENSION(nbpt) :: flood_diag_old !! Reservoir in the beginning of the time step !! For water budget check in the three routing reservoirs (positive if input > output) !! Net fluxes averaged over each grid cell in kg/m^2/dt REAL(r_std), DIMENSION(nbpt) :: netflow_stream_diag !! Input - Output flow to stream reservoir REAL(r_std), DIMENSION(nbpt) :: netflow_fast_diag !! Input - Output flow to fast reservoir REAL(r_std), DIMENSION(nbpt) :: netflow_slow_diag !! Input - Output flow to slow reservoir ! REAL(r_std), DIMENSION(nbpt,nbasmax) :: stemp_total_tend, stemp_advec_tend, stemp_relax_tend ! !_ ================================================================================================================================ ! Save reservoirs in beginning of time step to calculate the water budget fast_diag_old = fast_diag slow_diag_old = slow_diag stream_diag_old = stream_diag lake_diag_old = lake_diag pond_diag_old = pond_diag flood_diag_old = flood_diag ! !! 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(:) = tempdiag(:,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 tempdiag_mean(:,:) = tempdiag_mean(:,:) + tempdiag(:,:)*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 ! !! Computes the transport of water in the various reservoirs ! CALL routing_hr_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, & & tempdiag_mean, reinf_slope, lakeinflow_mean, returnflow_mean, reinfiltration_mean, & & irrigation_mean, riverflow_mean, coastalflow_mean, hydrographs, slowflow_diag, flood_frac, & & flood_res, netflow_stream_diag, netflow_fast_diag, netflow_slow_diag, & & stemp_total_tend, stemp_advec_tend, stemp_relax_tend) ! !! Responsible for storing the water in lakes ! CALL routing_hr_lake(nbpt, dt_routing, lakeinflow_mean, humrel_mean, return_lakes) ! returnflow_mean(:) = returnflow_mean(:) + return_lakes(:) 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 tempdiag_mean(:,:) = zero vegtot_mean(:) = zero ! Change the units of the routing fluxes from kg/dt_routing into kg/dt_sechiba hydrographs(:) = hydrographs(:)/dt_routing*dt_sechiba HTUhgmon(:,:) = HTUhgmon(:,:)/dt_routing*dt_sechiba slowflow_diag(:) = slowflow_diag(:)/dt_routing*dt_sechiba ! Change the units of the routing fluxes from kg/m^2/dt_routing into kg/m^2/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_routing to m^3/dt_sechiba riverflow_mean(:) = riverflow_mean(:)/dt_routing*dt_sechiba/mille coastalflow_mean(:) = coastalflow_mean(:)/dt_routing*dt_sechiba/mille ! Water budget residu of the three routing reservoirs (in kg/m^2/s) ! Note that these diagnostics are done using local variables only calculated ! during the time steps when the routing is calculated CALL xios_orchidee_send_field("wbr_stream",(stream_diag - stream_diag_old - netflow_stream_diag)/dt_routing) CALL xios_orchidee_send_field("wbr_fast", (fast_diag - fast_diag_old - netflow_fast_diag)/dt_routing) CALL xios_orchidee_send_field("wbr_slow", (slow_diag - slow_diag_old - netflow_slow_diag)/dt_routing) CALL xios_orchidee_send_field("wbr_lake", (lake_diag - lake_diag_old - & lakeinflow_mean + return_lakes)/dt_routing) CALL xios_orchidee_send_field("StreamT_TotTend", stemp_total_tend) CALL xios_orchidee_send_field("StreamT_AdvTend", stemp_advec_tend) CALL xios_orchidee_send_field("StreamT_RelTend", stemp_relax_tend) 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("mask_coast",mask_coast) IF ( do_irrigation ) THEN CALL xios_orchidee_send_field("irrigmap",irrigated) ENDIF IF ( do_floodplains ) THEN !! May be improved by performing the operation with XIOS floodmap(:) = 0.0 DO ig=1,nbpt floodmap(ig) = SUM(floodplains(ig,:)) / (area(ig)*contfrac(ig)) END DO CALL xios_orchidee_send_field("floodmap",floodmap) ENDIF IF ( doswamps ) THEN CALL xios_orchidee_send_field("swampmap",swamp) ENDIF ! ! Water storage in reservoirs [kg/m^2] 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("floodr",flood_diag) CALL xios_orchidee_send_field("floodh",flood_height) ! FLOODPLAINS CALL xios_orchidee_send_field("flood_frac",flood_frac) ! Difference between the end and the beginning of the routing time step [kg/m^2] CALL xios_orchidee_send_field("delfastr", fast_diag - fast_diag_old) CALL xios_orchidee_send_field("delslowr", slow_diag - slow_diag_old) CALL xios_orchidee_send_field("delstreamr", stream_diag - stream_diag_old) CALL xios_orchidee_send_field("dellaker", lake_diag - lake_diag_old) CALL xios_orchidee_send_field("delpondr", pond_diag - pond_diag_old) CALL xios_orchidee_send_field("delfloodr", flood_diag - flood_diag_old) ! Water fluxes converted from kg/m^2/dt_sechiba into kg/m^2/s CALL xios_orchidee_send_field("irrigation",irrigation/dt_sechiba) CALL xios_orchidee_send_field("netirrig",irrig_netereq/dt_sechiba) CALL xios_orchidee_send_field("riversret",returnflow/dt_sechiba) CALL xios_orchidee_send_field("reinfiltration",reinfiltration/dt_sechiba) ! Transform from kg/dt_sechiba into m^3/s CALL xios_orchidee_send_field("hydrographs",hydrographs/mille/dt_sechiba) CALL xios_orchidee_send_field("htuhgmon",HTUhgmon/mille/dt_sechiba) CALL xios_orchidee_send_field("htutempmon",HTUtempmon) CALL xios_orchidee_send_field("hydrotemp", hydrotemp) CALL xios_orchidee_send_field("streamlimit", streamlimit) 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. xios_orchidee_ok) THEN 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/mille, 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 CALL histwrite_p(hist_id, 'SurfStor', kjit, flood_diag+pond_diag+lake_diag, nbpt, index) CALL histwrite_p(hist_id, 'Dis', kjit, hydrographs/mille, nbpt, index) ! CALL histwrite_p(hist_id, 'slowr', kjit, slow_diag, nbpt, index) CALL histwrite_p(hist_id, 'fastr', kjit, fast_diag, nbpt, index) CALL histwrite_p(hist_id, 'streamr', kjit, stream_diag, nbpt, index) CALL histwrite_p(hist_id, 'lakevol', kjit, lake_diag, nbpt, index) CALL histwrite_p(hist_id, 'pondr', kjit, pond_diag, 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/mille, 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 ! CALL histwrite_p(hist2_id, 'SurfStor', kjit, flood_diag+pond_diag+lake_diag, nbpt, index) CALL histwrite_p(hist2_id, 'Dis', kjit, hydrographs/mille, nbpt, index) ! ENDIF ENDIF ENDIF ! ! END SUBROUTINE routing_highres_main !! ============================================================================================================================= !! SUBROUTINE: routing_highres_finalize !! !>\BRIEF Write to restart file !! !! DESCRIPTION: Write module variables to restart file !! !! RECENT CHANGE(S) !! !! REFERENCE(S) !! !! FLOWCHART !! \n !_ ============================================================================================================================== SUBROUTINE routing_highres_finalize( kjit, nbpt, rest_id, flood_frac, flood_res ) IMPLICIT NONE !! 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 !_ ================================================================================================================================ ! ! Write restart variables ! CALL restput_p (rest_id, 'routingcounter', kjit, time_counter) CALL restput_p (rest_id, 'streamtcst', kjit, stream_tcst) CALL restput_p (rest_id, 'slowtcst', kjit, slow_tcst) CALL restput_p (rest_id, 'fasttcst', kjit, fast_tcst) CALL restput_p (rest_id, 'floodtcst', kjit, flood_tcst) CALL restput_p (rest_id, 'swampcst', kjit, swamp_cst) CALL restput_p (rest_id, 'lim_floodcri', kjit, lim_floodcri) CALL restput_p (rest_id, 'nbasmax', kjit, nbasmax) CALL restput_p (rest_id, 'nbasmon', kjit, nbasmon) CALL restput_p (rest_id, 'inflows', kjit, inflows) 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, 'routenbintobas', nbp_glo, nbasmax, 1, kjit, REAL(route_nbintobas,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, 'topoindex_stream', nbp_glo, nbasmax, 1, kjit, stream_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, nbasmax, 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, 'fasttemp', nbp_glo, nbasmax, 1, kjit, fast_temp, 'scatter', nbp_glo, index_g) CALL restput_p (rest_id, 'slowtemp', nbp_glo, nbasmax, 1, kjit, slow_temp, 'scatter', nbp_glo, index_g) CALL restput_p (rest_id, 'streamtemp', nbp_glo, nbasmax, 1, kjit, stream_temp, '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, 'htuhgmon', nbp_glo, nbasmon, 1, kjit, HTUhgmon, 'scatter', nbp_glo, index_g) CALL restput_p (rest_id, 'slowflow_diag', nbp_glo, 1, 1, kjit, slowflow_diag, 'scatter', nbp_glo, index_g) CALL restput_p (rest_id, 'hydrotemp', nbp_glo, 1, 1, kjit, hydrotemp, 'scatter', nbp_glo, index_g) CALL restput_p (rest_id, 'htutempmon', nbp_glo, nbasmon, 1, kjit, HTUtempmon, '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) CALL restput_p (rest_id, 'tempdiag_route', nbp_glo, ngrnd, 1, kjit, tempdiag_mean, 'scatter', nbp_glo, index_g) CALL restput_p (rest_id, 'gridrephtu', nbp_glo, 1, 1, kjit, REAL(hydrodiag,r_std), 'scatter', nbp_glo, index_g) CALL restput_p (rest_id, 'htudiag', nbp_glo, nbasmon, 1, kjit, REAL(HTUdiag_loc,r_std), '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, nbasmax, 1, kjit, floodplains, 'scatter', nbp_glo, index_g) CALL restput_p (rest_id, 'floodcri', nbp_glo, nbasmax, 1, kjit, floodcri, 'scatter', nbp_glo, index_g) CALL restput_p (rest_id, 'floodp_beta', nbp_glo, nbasmax, 1, kjit, fp_beta, 'scatter', nbp_glo, index_g) ENDIF IF ( dofloodoverflow ) THEN CALL restput_p (rest_id, 'orog_min', nbp_glo, nbasmax, 1,kjit,orog_min, 'scatter', nbp_glo, index_g) END IF IF ( doswamps ) THEN CALL restput_p (rest_id, 'swamp', nbp_glo, 1, 1, kjit, swamp, 'scatter', nbp_glo, index_g) ENDIF END SUBROUTINE routing_highres_finalize !! ================================================================================================================================ !! SUBROUTINE : routing_hr_init !! !>\BRIEF This subroutine allocates the memory and get the fixed fields from the restart file. !! !! DESCRIPTION (definitions, functional, design, flags) : None !! !! RECENT CHANGE(S): None !! !! MAIN OUTPUT VARIABLE(S): !! !! REFERENCES : None !! !! FLOWCHART :None !! \n !_ ================================================================================================================================ SUBROUTINE routing_hr_init(kjit, nbpt, index, returnflow, reinfiltration, irrigation, & & riverflow, coastalflow, flood_frac, flood_res, tempdiag, rest_id) ! IMPLICIT NONE ! ! interface description ! !! 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,ngrnd),INTENT(in) :: tempdiag !! Temperature profile in soil INTEGER(i_std), INTENT(in) :: rest_id !! Restart file identifier (unitless) ! !! 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 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) ! !! 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), ALLOCATABLE, DIMENSION(:) :: tmp_real ! REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: tmp_real_g2 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, im, ib, ipn(1), nbhtumon !! Indices (unitless) REAL(r_std) :: nbasmon_tmp, nbasmax_tmp, inflows_tmp !_ ================================================================================================================================ ! ! ! 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 = dt_sechiba CALL getin_p('DT_ROUTING', dt_routing) ! ! ! !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. IF ( do_floodplains ) CALL getin_p('DO_FLOODINFILT', dofloodinfilt) ! !Config Key = CONDUCT_FACTOR_FP !Config Desc = Adjustment factor for floodplains reinfiltration !Config If = RIVER_ROUTING !Config Def = n !Config Help = Factor used to reduce the infiltration from the !Config floodplains. For a value of 1, the infiltration is !Config unchanged, for a value of 0 there is no infiltration. !Config Units = - ! conduct_factor = 1.0 IF ( do_floodplains ) CALL getin_p('CONDUCT_FACTOR_FP', conduct_factor) ! ! !Config Key = DO_FLOODOVERFLOW !Config Desc = Should floodplains overflow to upstream HTUs floodplains !Config If = RIVER_ROUTING !Config Def = n !Config Help = This parameters allows the user to ask the model !Config to take into account the overflow of the !Config floodplains. The water can flow to the upstream !Config floodplains reservoir if the current flood height !Config is higher than the upstream one. !Config Units = [FLAG] ! dofloodoverflow = .FALSE. IF ( do_floodplains ) CALL getin_p('DO_FLOODOVERFLOW', dofloodoverflow) ! !Config Key = OVERFLOW_REPETITION !Config Desc = Repetition of overflow at each routing time step !Config If = RIVER_ROUTING !Config Def = n !Config Help = This parameters allows the user to ask the model !Config repeat the overflow a certain amount of time !Config in order to have more stability with lower !Config overflow time step. !Config Units = [FLAG] ! overflow_repetition = 1 IF ( do_floodplains ) CALL getin_p('OVERFLOW_REPETITION', overflow_repetition) ! !Config Key = R_FLOODMAX !Config Desc = Maximal values for R factor !Config If = DO_FLOODPLAINS !Config Def = 0.5 !Config Help = R is the factor of reduction of the stream discharge !Config if there is floodplains. This is the maximal value !Config when the HTU is fully filled. !Config R = 1 -> discharge = 0 !Config R = 0 -> Maximal discharge ! rfloodmax = 0.5 IF ( do_floodplains ) CALL getin_p('R_FLOODMAX', rfloodmax) ! !Config Key = OVERFLOW_TCST !Config Desc = Time Constant for overflow in day !Config If = DO_FLOODPLAINS !Config Def = 1 !Config Help = OVERFLOW_TCST is the time constant !Config For the floodplains overflow ! overflow_tcst = 1 IF ( do_floodplains ) CALL getin_p('OVERFLOW_TCST', overflow_tcst) ! !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 = 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 = [-] ! ANTHONY OLD FLOODPLAINS !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] ! ANTHONY OLD FLOODPLAINS !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 ! ! If the routing_graph file is available we will extract the information in the dimensions ! and parameters. ! !Config Key = ROUTING_FILE !Config Desc = Name of file which contains the routing information graph on the model grid !Config If = RIVER_ROUTING !Config Def = routing.nc !Config Help = The file provided here should allows to route the water from one HTU !Config to another. The RoutingPP code needs to be used in order to generate !Config the routing graph for the model grid. !Config More details on : https://gitlab.in2p3.fr/ipsl/lmd/intro/routingpp !Config Units = [FILE] ! !graphfilename = 'routing_graph.nc' !CALL getin('ROUTING_FILE',graphfilename) !CALL routing_hr_graphinfo(graphfilename, nbasmax, inflows, nbasmon, undef_graphfile, stream_tcst, fast_tcst, slow_tcst, & & flood_tcst, swamp_cst, lim_floodcri) ! At this stage we could have an option to force reading of graph ! ! Constants which can be in the restart file ! var_name ="routingcounter" CALL ioconf_setatt_p('UNITS', 's') CALL ioconf_setatt_p('LONG_NAME','Time counter for the routing scheme') CALL restget_p (rest_id, var_name, kjit, .TRUE., zero, time_counter) CALL setvar_p (time_counter, val_exp, 'NO_KEYWORD', zero) ! Parameters which are in the restart file IF (stream_tcst .LE. 0 ) THEN var_name ="streamtcst" CALL ioconf_setatt_p('UNITS', 's/km') CALL ioconf_setatt_p('LONG_NAME','Time constant for the stream reservoir') CALL restget_p (rest_id, var_name, kjit, .TRUE., zero, stream_tcst) ENDIF IF (slow_tcst .LE. 0 ) THEN var_name ="slowtcst" CALL ioconf_setatt_p('UNITS', 's/km') CALL ioconf_setatt_p('LONG_NAME','Time constant for the slow reservoir') CALL restget_p (rest_id, var_name, kjit, .TRUE., zero, slow_tcst) ENDIF IF (fast_tcst .LE. 0 ) THEN var_name ="fasttcst" CALL ioconf_setatt_p('UNITS', 's/km') CALL ioconf_setatt_p('LONG_NAME','Time constant for the fast reservoir') CALL restget_p (rest_id, var_name, kjit, .TRUE., zero, fast_tcst) ENDIF IF (flood_tcst .LE. 0 ) THEN var_name ="floodtcst" CALL ioconf_setatt_p('UNITS', 's/km') CALL ioconf_setatt_p('LONG_NAME','Time constant for the flood reservoir') CALL restget_p (rest_id, var_name, kjit, .TRUE., zero, flood_tcst) ENDIF IF (swamp_cst .LE. 0 ) THEN var_name ="swampcst" CALL ioconf_setatt_p('UNITS', '-') CALL ioconf_setatt_p('LONG_NAME','Fraction of the river transport that flows to the swamps') CALL restget_p (rest_id, var_name, kjit, .TRUE., zero, swamp_cst) ENDIF IF (lim_floodcri .LE. 0 ) THEN var_name ="lim_floodcri" CALL ioconf_setatt_p('UNITS', 'm') CALL ioconf_setatt_p('LONG_NAME','Minimal difference of orography consecutive floodplains HTUs') CALL restget_p (rest_id, var_name, kjit, .TRUE., zero, lim_floodcri) ENDIF ! ! Number of HTUs ! var_name ="nbasmax" CALL ioconf_setatt_p('UNITS', '-') CALL ioconf_setatt_p('LONG_NAME','Number of HTU per grid box') CALL restget_p (rest_id, var_name, kjit, .TRUE., zero, nbasmax_tmp) CALL routing_hr_restartconsistency(var_name, nbasmax, nbasmax_tmp) ! ! Number of inflows ! var_name ="inflows" CALL ioconf_setatt_p('UNITS', '-') CALL ioconf_setatt_p('LONG_NAME','Maximum number of inflows per HTU') CALL restget_p (rest_id, var_name, kjit, .TRUE., zero, inflows_tmp) CALL routing_hr_restartconsistency(var_name, inflows, inflows_tmp) ! ! Dimension of HTU monitoring variable ! var_name ="nbasmon" CALL ioconf_setatt_p('UNITS', '-') CALL ioconf_setatt_p('LONG_NAME','Number of HTU to be monitored') CALL restget_p (rest_id, var_name, kjit, .TRUE., zero, nbasmon_tmp) CALL routing_hr_restartconsistency(var_name, nbasmon, nbasmon_tmp) ! ! Continuation of extraction from restart file. ! ALLOCATE (routing_area_loc(nbpt,nbasmax), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for routing_area_loc','','') ALLOCATE (routing_area_glo(nbp_glo,nbasmax)) IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for routing_area_glo','','') IF ( .NOT. ReadGraph ) THEN 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 ENDIF CALL scatter(routing_area_glo,routing_area_loc) routing_area=>routing_area_loc IF ( do_floodplains ) THEN !!!!!!!!!!!!!!!!!!!!!!!!!!! !! ANTHONY - BETA ALLOCATE (fp_beta_loc(nbpt,nbasmax), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for fp_beta_loc','','') ALLOCATE (fp_beta_glo(nbp_glo,nbasmax)) IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for fp_beta_glo','','') IF ( .NOT. ReadGraph ) THEN IF (is_root_prc) THEN var_name = 'floodp_beta' CALL ioconf_setatt('UNITS', '-') CALL ioconf_setatt('LONG_NAME','Beta parameter for floodplains') CALL restget (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., fp_beta_glo, "gather", nbp_glo, index_g) ENDIF CALL scatter(fp_beta_glo,fp_beta_loc) fp_beta=>fp_beta_loc END IF !!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!! !! ANTHONY - h0 - floodcri ALLOCATE (floodcri_loc(nbpt,nbasmax), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for floodcri_loc','','') ALLOCATE (floodcri_glo(nbp_glo,nbasmax)) IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for floodcri_glo','','') IF ( .NOT. ReadGraph ) THEN IF (is_root_prc) THEN var_name = 'floodcri' CALL ioconf_setatt('UNITS', 'mm') CALL ioconf_setatt('LONG_NAME','Height of complete flood') CALL restget (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., floodcri_glo, "gather", nbp_glo, index_g) END IF CALL scatter(floodcri_glo,floodcri_loc) floodcri=>floodcri_loc ENDIF END IF ALLOCATE (tmp_real_g(nbp_glo,nbasmax), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for tmp_real_g','','') ALLOCATE (route_togrid_loc(nbpt,nbasmax), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for route_togrid_loc','','') ALLOCATE (route_togrid_glo(nbp_glo,nbasmax), stat=ier) ! used in global in routing_hr_flow IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for route_togrid_glo','','') IF ( .NOT. ReadGraph ) THEN 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_hr_flow CALL scatter(route_togrid_glo,route_togrid_loc) route_togrid=>route_togrid_loc ENDIF ! ALLOCATE (route_tobasin_loc(nbpt,nbasmax), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for route_tobasin_loc','','') ALLOCATE (route_tobasin_glo(nbp_glo,nbasmax), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for route_tobasin_glo','','') IF ( .NOT. ReadGraph ) THEN 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 num_largest = COUNT(route_tobasin_glo .EQ. nbasmax+3) ENDIF CALL scatter(route_tobasin_glo,route_tobasin_loc) CALL bcast(num_largest) route_tobasin=>route_tobasin_loc ENDIF ! ! nbintobasin ! ALLOCATE (route_nbintobas_loc(nbpt,nbasmax), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for route_nbintobas_loc','','') ALLOCATE (route_nbintobas_glo(nbp_glo,nbasmax), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for route_nbintobas_glo','','') IF ( .NOT. ReadGraph ) THEN 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 ENDIF ! ALLOCATE (global_basinid_loc(nbpt,nbasmax), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for global_basinid_loc','','') ALLOCATE (global_basinid_glo(nbp_glo,nbasmax), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for global_basinid_glo','','') IF ( .NOT. ReadGraph ) THEN 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 ENDIF ! ALLOCATE (topo_resid_loc(nbpt,nbasmax), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for topo_resid_loc','','') ALLOCATE (topo_resid_glo(nbp_glo,nbasmax), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for topo_resid_glo','','') IF ( .NOT. ReadGraph ) THEN IF (is_root_prc) THEN var_name = 'topoindex' CALL ioconf_setatt('UNITS', 'km') 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 ENDIF ! ALLOCATE (stream_resid_loc(nbpt,nbasmax), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for stream_resid_loc','','') ALLOCATE (stream_resid_glo(nbp_glo,nbasmax), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for stream_resid_glo','','') IF ( .NOT. ReadGraph ) THEN IF (is_root_prc) THEN var_name = 'topoindex_stream' CALL ioconf_setatt('UNITS', 'km') CALL ioconf_setatt('LONG_NAME','Topographic index of the residence time') CALL restget (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., stream_resid_glo, "gather", nbp_glo, index_g) stream_maxresid=MAXVAL(stream_resid_glo, MASK=stream_resid_glo .LT. undef_graphfile) ENDIF CALL bcast(stream_maxresid) CALL scatter(stream_resid_glo,stream_resid_loc) stream_resid=>stream_resid_loc ENDIF ! ALLOCATE (fast_reservoir(nbpt,nbasmax), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','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_hr_init','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_hr_init','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_hr_init','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_hr_init','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, nbasmax), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','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, nbasmax, 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_hr_init','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_hr_init','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_hr_init','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_hr_init','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 (floodmap(nbpt), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for floodmap','','') ALLOCATE (floodplains_loc(nbpt,nbasmax), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for floodplains_loc','','') ALLOCATE (floodplains_glo(nbp_glo,nbasmax), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for floodplains_glo','','') IF ( .NOT. ReadGraph ) THEN IF (is_root_prc) THEN var_name = 'floodplains' CALL ioconf_setatt_p('UNITS', 'm^2') CALL ioconf_setatt_p('LONG_NAME','Surface which can be flooded') CALL restget (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., floodplains_glo, "gather", nbp_glo, index_g) END IF CALL scatter(floodplains_glo,floodplains_loc) floodplains=>floodplains_loc END IF ENDIF !!! !!! ANTHONY : OVERFLOW !!! IF ( dofloodoverflow) THEN ALLOCATE (orog_min_loc(nbpt,nbasmax), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for orog_min_loc','','') ALLOCATE (orog_min_glo(nbp_glo,nbasmax), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for orog_min_glo','','') !! IF ( .NOT. ReadGraph ) THEN IF (is_root_prc) THEN var_name = 'orog_min' CALL ioconf_setatt('UNITS', 'm') CALL ioconf_setatt('LONG_NAME','HTU minimum orography') CALL restget (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., orog_min_glo, "gather", nbp_glo, index_g) END IF CALL scatter(orog_min_glo,orog_min_loc) orog_min=>orog_min_loc END IF ! ALLOCATE (route_innum_loc(nbpt,nbasmax), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for route_innum_loc','','') ALLOCATE (route_innum_glo(nbp_glo,nbasmax), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for route_innum_glo','','') CALL scatter(route_innum_glo,route_innum_loc) route_innum=>route_innum_loc ! ALLOCATE (route_ingrid_loc(nbpt,nbasmax, inflows), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for route_ingrid_loc','','') ALLOCATE (route_ingrid_glo(nbp_glo,nbasmax,inflows), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for route_ingrid_glo','','') CALL scatter(route_ingrid_glo,route_ingrid_loc) route_ingrid=>route_ingrid_loc ! ALLOCATE (route_inbasin_loc(nbpt,nbasmax, inflows), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for route_inbasin_loc','','') ALLOCATE (route_inbasin_glo(nbp_glo,nbasmax, inflows), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for route_inbasin_glo','','') CALL scatter(route_inbasin_glo,route_inbasin_loc) route_inbasin=>route_inbasin_loc END IF !!! !!! !!! IF ( doswamps ) THEN ALLOCATE (swamp(nbpt), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','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_hr_init','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_hr_init','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_hr_init','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_hr_init','Pb in allocate for irrigation_mean','','') ALLOCATE (irrig_netereq(nbpt), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','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 (fast_temp(nbpt,nbasmax), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for fast_temp','','') var_name = 'fasttemp' CALL ioconf_setatt_p('UNITS', 'K') CALL ioconf_setatt_p('LONG_NAME','Water temperature in the fast reservoir') CALL restget_p (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., fast_temp, "gather", nbp_glo, index_g) ALLOCATE (slow_temp(nbpt,nbasmax), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for slow_temp','','') var_name = 'slowtemp' CALL ioconf_setatt_p('UNITS', 'K') CALL ioconf_setatt_p('LONG_NAME','Water temperature in the slow reservoir') CALL restget_p (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., slow_temp, "gather", nbp_glo, index_g) IF ( COUNT(fast_temp == val_exp) == nbpt*nbasmax ) THEN CALL groundwatertemp(nbpt, nbasmax, ngrnd, tempdiag, znt, dlt, fast_temp, slow_temp) ENDIF ALLOCATE (stream_temp(nbpt,nbasmax), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for stream_temp','','') var_name = 'streamtemp' CALL ioconf_setatt_p('UNITS', 'K') CALL ioconf_setatt_p('LONG_NAME','Water temperature in the stream reservoir') CALL restget_p (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., stream_temp, "gather", nbp_glo, index_g) IF ( COUNT(stream_temp == val_exp) == nbpt*nbasmax ) THEN DO ig=1,nbpt stream_temp(ig,:) = tempdiag(ig,1) ENDDO ENDIF ALLOCATE (riverflow_mean(nbpt), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','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_hr_init','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(zlt-2)) floodtemp_lev = ipn(1) ALLOCATE (floodtemp(nbpt), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for floodtemp','','') floodtemp(:) = tempdiag(:,floodtemp_lev) ALLOCATE(hydrographs(nbpt), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for hydrographs','','') var_name = 'hydrographs' CALL ioconf_setatt_p('UNITS', 'kg/dt_sechiba') 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(hydrotemp(nbpt), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for hydrotemp','','') var_name = 'hydrotemp' CALL ioconf_setatt_p('UNITS', 'K') CALL ioconf_setatt_p('LONG_NAME','Temperature of most significant river of grid') CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., hydrotemp, "gather", nbp_glo, index_g) CALL setvar_p (hydrotemp, val_exp, 'NO_KEYWORD', ZeroCelsius) ALLOCATE(slowflow_diag(nbpt), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for slowflow_diag','','') var_name = 'slowflow_diag' CALL ioconf_setatt_p('UNITS', 'kg/dt_sechiba') 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) ! ! Grid diagnostic at representative HTU ! ALLOCATE(hydrodiag_loc(nbpt),hydrodiag_glo(nbp_glo),stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for hydrodiag_glo','','') IF ( .NOT. ReadGraph ) THEN IF (is_root_prc) THEN ALLOCATE(tmp_real(nbp_glo)) var_name = 'gridrephtu' CALL ioconf_setatt('UNITS', '-') CALL ioconf_setatt('LONG_NAME','Representative HTU for the grid') CALL restget(rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE.,tmp_real, "gather", nbp_glo, index_g) hydrodiag_glo(:) = 1 WHERE ( tmp_real .LT. val_exp ) hydrodiag_glo = NINT(tmp_real) ENDWHERE DEALLOCATE(tmp_real) ENDIF CALL scatter(hydrodiag_glo, hydrodiag_loc) ENDIF ! ! Station diagnostics ! ALLOCATE(HTUdiag_loc(nbpt,nbasmon), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for HTUdiag_loc','','') ALLOCATE(HTUdiag_glo(nbp_glo,nbasmon), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for HTUdiag_glo','','') ALLOCATE(tmp_real_g2(nbp_glo,nbasmon), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for tmp_real_g2','','') ! IF (is_root_prc) THEN var_name = 'htudiag' CALL ioconf_setatt('UNITS', 'index') CALL ioconf_setatt('LONG_NAME','Index of HTU to be monitored') CALL restget(rest_id, var_name, nbp_glo, nbasmon, 1, kjit, .TRUE., tmp_real_g2, "gather", nbp_glo, index_g) HTUdiag_glo(:,:) = -1 WHERE ( tmp_real_g2 .LT. val_exp ) HTUdiag_glo = NINT(tmp_real_g2) ENDWHERE nbhtumon = 0 DO ig=1,nbp_glo DO im=1,nbasmon IF ( HTUdiag_glo(ig,im) > 0 ) THEN nbhtumon = nbhtumon + 1 ENDIF ENDDO ENDDO WRITE(numout,*) "After restget : Found a total of ", nbhtumon, " HTUs to be monitored and written into HTUhgmon" ENDIF CALL scatter(HTUdiag_glo, HTUdiag_loc) CALL bcast(nbhtumon) DEALLOCATE(tmp_real_g2) ! ALLOCATE(HTUhgmon(nbpt,nbasmon), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for HTUhgmon','','') ALLOCATE(HTUhgmon_glo(nbp_glo,nbasmon), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for HTUhgmon_glo','','') ! IF (is_root_prc) THEN var_name = 'htuhgmon' CALL ioconf_setatt('UNITS', 'kg/dt_sechiba') CALL ioconf_setatt('LONG_NAME','Hydrograph at selected HTU of grid') CALL restget(rest_id, var_name, nbp_glo, nbasmon, 1, kjit, .TRUE., HTUhgmon_glo, "gather", nbp_glo, index_g) WHERE ( HTUhgmon_glo .GE. val_exp ) HTUhgmon_glo = zero ENDWHERE ENDIF CALL scatter(HTUhgmon_glo, HTUhgmon) DEALLOCATE(HTUhgmon_glo) ! ! Restart of the temperature monitoring ! ALLOCATE(HTUtempmon(nbpt,nbasmon), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for HTUtempmon','','') ALLOCATE(HTUtempmon_glo(nbp_glo,nbasmon), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for HTUtempmon_glo','','') ! IF (is_root_prc) THEN var_name = 'htutempmon' CALL ioconf_setatt('UNITS', 'K') CALL ioconf_setatt('LONG_NAME','Temperature at selected HTU of grid') CALL restget(rest_id, var_name, nbp_glo, nbasmon, 1, kjit, .TRUE., HTUtempmon_glo, "gather", nbp_glo, index_g) WHERE ( HTUtempmon_glo .GE. val_exp ) HTUtempmon_glo = ZeroCelsius ENDWHERE HTUtempmon_glo(:,:) = ZeroCelsius ENDIF CALL scatter(HTUtempmon_glo, HTUtempmon) DEALLOCATE(HTUtempmon_glo) ! ! 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), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for fast_diag,..','','') fast_diag(:) = zero slow_diag(:) = zero stream_diag(:) = zero flood_diag(:) = zero pond_diag(:) = zero lake_diag(:) = 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 better 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_hr_init','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_hr_init','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_hr_init','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_hr_init','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_hr_init','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_hr_init','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_hr_init','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_hr_init','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_hr_init','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) ! ALLOCATE(tempdiag_mean(nbpt,ngrnd), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_init','Pb in allocate for tempdiag_mean','','') var_name = 'tempdiag_route' CALL ioconf_setatt_p('UNITS', 'K') CALL ioconf_setatt_p('LONG_NAME','Mean temperature profile') CALL restget_p (rest_id, var_name, nbp_glo, ngrnd, 1, kjit, .TRUE., tempdiag_mean, "gather", nbp_glo, index_g) CALL setvar_p (tempdiag_mean, val_exp, 'NO_KEYWORD', Zero) ! DEALLOCATE(tmp_real_g) ! ! Other variables ! ALLOCATE(streamlimit(nbpt), stat=ier) ! END SUBROUTINE routing_hr_init ! !! ================================================================================================================================ !! SUBROUTINE : routing_hr_restartconsistency !! !>\BRIEF : This subroutine will verify that the important dimensions for the routing exist in the restart file !! or can be read from the routing_graph.nc file. This ensures that if the restart and graph file are !! not consistent the latter is read and overwrite whatever was in the restart. Then the user will be !! using the new routing graph from the routing_graph.nc and not whatever is in the restart. !! \n !_ ================================================================================================================================ SUBROUTINE routing_hr_restartconsistency(varname, dimgraph, dimrestart) ! IMPLICIT NONE ! !! INPUT VARIABLES CHARACTER(LEN=80), INTENT(in) :: varname !! Name of dimension INTEGER(i_std), INTENT(inout) :: dimgraph !! Dimension read in the graph file and also final result. REAL(r_std), INTENT(in) :: dimrestart !! Dimension read in the restart file. ! ! Case when the routing_hr_graphinfo could not get any information from the routing_graph.nc IF ( dimgraph < un ) THEN IF ( dimrestart > zero ) THEN ! Only information from the restart. dimgraph = NINT(dimrestart) ELSE WRITE(*,*) "Problem : No information in the routing_graph file and no routing information in restart ", TRIM(varname) CALL ipslerr(3,'routing_hr_restartconsistency',& 'No routing_graph file availble and no information in restart.', & 'Cannot perform routing in ORCHIDEE.', ' ') ENDIF ! Information from the routing_graph.nc file exists ! ELSE IF ( dimgraph .NE. NINT(dimrestart) ) THEN WRITE(*,*) "Problem for ", TRIM(varname)," in restart is not the same as in routing_graph.nc " WRITE(*,*) "Value of ", TRIM(varname), " in restart file : ", dimrestart WRITE(*,*) "Value of ", TRIM(varname), " in routing_graph.nc file : ", dimgraph CALL ipslerr(2,'routing_hr_restartconsistency',& 'The value of dimension provided is not consistant with the one in routing_graph file.', & 'We will read a new graph from the given file.', ' ') ReadGraph = .TRUE. ELSE !! Nothing to do ENDIF ENDIF END SUBROUTINE routing_hr_restartconsistency !! ================================================================================================================================ !! SUBROUTINE : routing_highres_clear !! !>\BRIEF : This subroutine deallocates the block memory previously allocated. !! \n !_ ================================================================================================================================ SUBROUTINE routing_highres_clear() 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(stream_resid_loc)) DEALLOCATE(stream_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(stream_resid_glo)) DEALLOCATE(stream_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(fast_temp)) DEALLOCATE(fast_temp) IF (ALLOCATED(slow_temp)) DEALLOCATE(slow_temp) IF (ALLOCATED(stream_temp)) DEALLOCATE(stream_temp) 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(tempdiag_mean)) DEALLOCATE(tempdiag_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(hydrographs)) DEALLOCATE(hydrographs) IF (ALLOCATED(hydrotemp)) DEALLOCATE(hydrotemp) IF (ALLOCATED(HTUhgmon)) DEALLOCATE(HTUhgmon) IF (ALLOCATED(HTUhgmon_glo)) DEALLOCATE(HTUhgmon_glo) IF (ALLOCATED(HTUtempmon)) DEALLOCATE(HTUtempmon) IF (ALLOCATED(HTUtempmon_glo)) DEALLOCATE(HTUtempmon_glo) IF (ALLOCATED(slowflow_diag)) DEALLOCATE(slowflow_diag) IF (ALLOCATED(irrigation_mean)) DEALLOCATE(irrigation_mean) IF (ALLOCATED(irrigated)) DEALLOCATE(irrigated) IF (ALLOCATED(floodplains_glo)) DEALLOCATE(floodplains_glo) IF (ALLOCATED(floodplains_loc)) DEALLOCATE(floodplains_loc) 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(route_innum_loc)) DEALLOCATE(route_innum_loc) IF (ALLOCATED(route_ingrid_loc)) DEALLOCATE(route_ingrid_loc) IF (ALLOCATED(route_inbasin_loc)) DEALLOCATE(route_inbasin_loc) IF (ALLOCATED(route_innum_glo)) DEALLOCATE(route_innum_glo) IF (ALLOCATED(route_ingrid_glo)) DEALLOCATE(route_ingrid_glo) IF (ALLOCATED(route_inbasin_glo)) DEALLOCATE(route_inbasin_glo) ! IF (ALLOCATED(orog_min_loc)) DEALLOCATE(orog_min_loc) IF (ALLOCATED(orog_min_glo)) DEALLOCATE(orog_min_glo) ! IF (ALLOCATED(floodcri_loc)) DEALLOCATE(floodcri_loc) IF (ALLOCATED(floodcri_glo)) DEALLOCATE(floodcri_glo) IF (ALLOCATED(fp_beta_loc)) DEALLOCATE(fp_beta_loc) IF (ALLOCATED(fp_beta_glo)) DEALLOCATE(fp_beta_glo) END SUBROUTINE routing_highres_clear ! !! ================================================================================================================================ !! SUBROUTINE : routing_hr_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_hr_flow(nbpt, dt_routing, lalo, floodout, runoff, drainage, & & vegtot, totnobio, transpot_mean, precip, humrel, k_litt, floodtemp, tempdiag, & & reinf_slope, lakeinflow, returnflow, reinfiltration, irrigation, riverflow, & & coastalflow, hydrographs, slowflow_diag, flood_frac, flood_res, & & netflow_stream_diag, netflow_fast_diag, netflow_slow_diag, & & stemp_total_tend, stemp_advec_tend, stemp_relax_tend) ! IMPLICIT NONE ! !! 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(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(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) :: tempdiag(nbpt,ngrnd) !! Soil temperature profiles (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) REAL(r_std), INTENT(out) :: lakeinflow(nbpt) !! Water inflow to the lakes (kg/dt) ! !! OUTPUT VARIABLES 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 module to allow re-evaporation (kg/m^2/dt_routing) 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_routing) 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_routing) 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_routing) 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) REAL(r_std), INTENT(out) :: netflow_stream_diag(nbpt) !! Input - Output flow to stream reservoir REAL(r_std), INTENT(out) :: netflow_fast_diag(nbpt) !! Input - Output flow to fast reservoir REAL(r_std), INTENT(out) :: netflow_slow_diag(nbpt) !! Input - Output flow to slow reservoir REAL(r_std), INTENT(out) :: stemp_total_tend(nbpt, nbasmax) !! Total tendency in GJ/s computed for the stream reservoir. REAL(r_std), INTENT(out) :: stemp_advec_tend(nbpt, nbasmax) !! Tendency (GJ/s) produced by advection REAL(r_std), INTENT(out) :: stemp_relax_tend(nbpt, nbasmax) !! Tendency (GJ/s) produced by relaxation ! !! LOCAL VARIABLES REAL(r_std), DIMENSION(nbpt, nbasmax) :: fast_flow !! Outflow from the fast reservoir (kg/dt) REAL(r_std), DIMENSION(nbpt, nbasmax) :: slow_flow !! Outflow from the slow reservoir (kg/dt) REAL(r_std), DIMENSION(nbpt, nbasmax) :: stream_flow !! Outflow from the stream reservoir (kg/dt) REAL(r_std), DIMENSION(nbpt, nbasmax) :: flood_flow !! Outflow from the floodplain reservoir (kg/dt) REAL(r_std), DIMENSION(nbpt, nbasmax) :: pond_inflow !! Inflow to the pond reservoir (kg/dt) REAL(r_std), DIMENSION(nbpt, nbasmax) :: pond_drainage !! Drainage from pond (kg/m^2/dt) REAL(r_std), DIMENSION(nbpt, nbasmax) :: flood_drainage !! Drainage from floodplains (kg/m^2/dt) REAL(r_std), DIMENSION(nbpt, nbasmax) :: return_swamp !! Inflow to the swamp (kg/dt) REAL(r_std), DIMENSION(nbpt, nbasmax) :: source REAL(r_std), DIMENSION(nbpt, nbasmax) :: ewh ! ! Irrigation per basin ! REAL(r_std), DIMENSION(nbpt, nbasmax) :: irrig_needs !! Total irrigation requirement (water requirements by the crop for its optimal growth) (kg) REAL(r_std), DIMENSION(nbpt, nbasmax) :: irrig_actual !! Possible irrigation according to the water availability in the reservoirs (kg) REAL(r_std), DIMENSION(nbpt, nbasmax) :: irrig_deficit !! Amount of water missing for irrigation (kg) REAL(r_std), DIMENSION(nbpt, nbasmax) :: irrig_adduct !! Amount of water carried over from other basins for irrigation (kg) ! ! The transport terms are over a larger indexing space so that outlfows to ocean and lakes do not generate out of bounds issues. ! Non existing HTU have their index set to zero and their memory will end-up in index 0 of transport. ! REAL(r_std), DIMENSION(nbpt, 0:nbasmax+3) :: transport !! Water transport between basins (kg/dt) REAL(r_std), DIMENSION(nbp_glo, 0:nbasmax+3) :: transport_glo !! Water transport between basins (kg/dt) REAL(r_std), DIMENSION(nbpt, 0:nbasmax+3) :: transport_temp !! Temperature transport between grids REAL(r_std), DIMENSION(nbp_glo, 0:nbasmax+3) :: transport_temp_glo !! Temperature transport global for transfers ! REAL(r_std) :: oldtemp REAL(r_std) :: oldstream INTEGER(i_std), SAVE :: nbunpy=0 ! REAL(r_std), DIMENSION(nbpt, nbasmax) :: floods !! Water flow in to the floodplains (kg/dt) REAL(r_std), DIMENSION(nbpt, nbasmax) :: potflood !! 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) 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) :: 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(nbp_glo) :: lake_overflow_g !! Removed water from lake reservoir on global grid (kg/gridcell/dt_routing) REAL(r_std), DIMENSION(nbpt) :: lake_overflow !! Removed water from lake reservoir on local grid (kg/gridcell/dt_routing) REAL(r_std), DIMENSION(nbpt) :: lake_overflow_coast !! lake_overflow distributed on coast gridcells, only diag(kg/gridcell/dt_routing) REAL(r_std) :: total_lake_overflow !! Sum of lake_overflow over full grid (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, im !! Indices (unitless) INTEGER(i_std) :: rtg, rtb, in, ing, inb,inf!! Indices (unitless) !INTEGER(i_std) :: numflood !! INTEGER(i_std) :: ier, negslow !! Error handling INTEGER(i_std), DIMENSION(20) :: negig, negib REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: fast_flow_g !! Outflow from the fast reservoir (kg/dt) REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: slow_flow_g !! Outflow from the slow reservoir (kg/dt) REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: stream_flow_g !! Outflow from the stream reservoir (kg/dt) REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: fast_temp_g !! Temperature of the fast reservoir (K) REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: slow_temp_g !! Temperature of the slow reservoir (K) REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: stream_temp_g !! Temperature of the stream reservoir (K) !REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: flood_height_g !! Floodplains height (m) !REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: flood_frac_bas_g !! Fraction of the HTU flooded REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: irrig_deficit_glo !! Amount of water missing for irrigation (kg) REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: stream_reservoir_glo !! Water amount in the stream reservoir (kg) !REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: flood_reservoir_glo !! Water amount in the stream reservoir (kg) REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: irrig_adduct_glo !! Amount of water carried over from other basins for irrigation (kg) REAL(r_std) :: reduced !! Discharge reduction due to floodplains REAL(r_std) :: htmp, hscale !! Water height scalingfor temperature relaxation REAL(r_std) :: krelax, den !! PARAMETERS LOGICAL, PARAMETER :: check_reservoir = .FALSE. !! Logical to choose if we write informations when a negative amount of water is occurring in a reservoir (true/false) !_ ================================================================================================================================ ! ! hscale = 1. CALL getin_p('ROUTING_HSCALEKH',hscale) ! transport(:,:) = zero transport_glo(:,:) = zero transport_temp(:,:) = zero !tp_00 its a transport, not a temperature !! transport_temp_glo(:,:) = zero !tp_00 irrig_netereq(:) = zero irrig_needs(:,:) = zero irrig_actual(:,:) = zero irrig_deficit(:,:) = zero irrig_adduct(:,:) = zero totarea(:) = zero totflood(:) = zero ! ! Compute all the fluxes ! DO ib=1,nbasmax DO ig=1,nbpt ! totarea(ig) = totarea(ig) + routing_area(ig,ib) totflood(ig) = totflood(ig) + flood_reservoir(ig,ib) ENDDO ENDDO ! !> 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 s/km) and topo_resid (in km) represents the time constant (s) !> 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 groundwatertemp(nbpt, nbasmax, ngrnd, tempdiag, znt, dlt, fast_temp, slow_temp) ! streamlimit(:) = zero ! DO ib=1,nbasmax DO ig=1,nbpt IF ( route_tobasin(ig,ib) .GT. 0 ) THEN ! ! Each of the fluxes is limited by the water in the reservoir and a small margin ! (min_reservoir) to avoid rounding errors. ! flow = MIN(fast_reservoir(ig,ib)/(topo_resid(ig,ib)*fast_tcst/dt_routing),& & fast_reservoir(ig,ib)-min_sechiba) fast_flow(ig,ib) = MAX(flow, zero) flow = MIN(slow_reservoir(ig,ib)/(topo_resid(ig,ib)*slow_tcst/dt_routing),& & slow_reservoir(ig,ib)-min_sechiba) slow_flow(ig,ib) = MAX(flow, zero) ! Need to adjust the reduction of the flow reduced = MAX(1-SQRT(MIN(flood_frac_bas(ig,ib),rfloodmax)), min_sechiba) ! Add the reduction flow parameter flow = stream_reservoir(ig,ib)/(stream_resid(ig,ib)*stream_tcst/dt_routing)*reduced flow = MIN(flow, stream_reservoir(ig,ib)-min_sechiba) stream_flow(ig,ib) = MAX(flow, zero) IF ( stream_flow(ig,ib) .GE. stream_reservoir(ig,ib)-min_sechiba .AND. stream_flow(ig,ib) > zero .AND. & & routing_area(ig,ib) > zero ) THEN streamlimit(ig) = streamlimit(ig)+1.0 ENDIF ! ELSE fast_flow(ig,ib) = zero slow_flow(ig,ib) = zero stream_flow(ig,ib) = zero ENDIF ENDDO ENDDO !- !- Compute the fluxes out of the floodplains and ponds if they exist. !- IF (do_floodplains .OR. doponds) THEN DO ig=1,nbpt IF (flood_frac(ig) .GT. min_sechiba) THEN !!!! ! PONDS : not actualized ! !flow = MIN(floodout(ig)*totarea(ig)*pond_frac(ig)/flood_frac(ig), pond_reservoir(ig)+totflood(ig)) !pondex = MAX(flow - pond_reservoir(ig), zero) !pond_reservoir(ig) = pond_reservoir(ig) - (flow - pondex) ! ! If demand was over reservoir size, we will take it out from floodplains ! !pond_excessflow(:) = zero !DO ib=1,nbasmax ! pond_excessflow(ib) = MIN(pondex*flood_frac_bas(ig,ib)/(flood_frac(ig)-pond_frac(ig)),& ! & flood_reservoir(ig,ib)) ! pondex = pondex - pond_excessflow(ib) !ENDDO ! !IF ( pondex .GT. min_sechiba) THEN ! WRITE(numout,*) "Unable to redistribute the excess pond outflow over the water available in the floodplain." ! WRITE(numout,*) "Pondex = ", pondex ! WRITE(numout,*) "pond_excessflow(:) = ", pond_excessflow(:) !ENDIF ! DO ib=1,nbasmax ! ! when ponds actualized : add pond_excessflow to flow ! This is the flow out of the reservoir due to ET (+ pond excessflow(ig), suppressed here) !flow = floodout(ig)*routing_area(ig,ib)*flood_frac_bas(ig,ib)/flood_frac(ig) flow = floodout(ig)*routing_area(ig,ib)*flood_frac_bas(ig,ib) ! flood_reservoir(ig,ib) = flood_reservoir(ig,ib) - flow ! ! IF (flood_reservoir(ig,ib) .LT. min_sechiba) THEN flood_reservoir(ig,ib) = zero ENDIF IF (pond_reservoir(ig) .LT. min_sechiba) THEN pond_reservoir(ig) = zero ENDIF ENDDO ENDIF ENDDO ENDIF !- !- Computing the drainage and outflow from floodplains !> Drainage from floodplains is depending on a averaged conductivity (k_litt) !> for saturated infiltration in the 'litter' layer. Flood_drainage will be !> a component of the total reinfiltration that leaves the routing scheme. !- IF (do_floodplains) THEN IF (dofloodinfilt) THEN DO ib=1,nbasmax DO ig=1,nbpt flood_drainage(ig,ib) = MAX(zero, MIN(flood_reservoir(ig,ib), & & flood_frac_bas(ig,ib)* routing_area(ig,ib) * k_litt(ig) * & & conduct_factor * dt_routing/one_day)) flood_reservoir(ig,ib) = flood_reservoir(ig,ib) - flood_drainage(ig,ib) ENDDO ENDDO ELSE DO ib=1,nbasmax DO ig=1,nbpt flood_drainage(ig,ib) = zero ENDDO ENDDO ENDIF !> Outflow from floodplains is computed depending a delay. This delay is characterized by a time constant !> function of the surface of the floodplains and the product of topo_resid and flood_tcst. flood_tcst !> has been calibrated through observations in the Niger Inner Delta (D'Orgeval, 2006). ! DO ib=1,nbasmax DO ig=1,nbpt IF ( route_tobasin(ig,ib) .GT. 0 ) THEN IF (flood_reservoir(ig,ib) .GT. min_sechiba) THEN flow = MIN(flood_reservoir(ig,ib)/(stream_resid(ig,ib)*flood_tcst/dt_routing),& & flood_reservoir(ig,ib)-min_sechiba) flow = MAX(flow, zero) ELSE flow = zero ENDIF flood_flow(ig,ib) = flow ELSE flood_flow(ig,ib) = zero ENDIF ENDDO ENDDO ELSE DO ib=1,nbasmax DO ig=1,nbpt flood_drainage(ig,ib) = zero flood_flow(ig,ib) = zero flood_reservoir(ig,ib) = zero ENDDO ENDDO ENDIF !- !- Computing drainage and inflow for ponds !> Drainage from ponds is computed in the same way than for floodplains. !> Reinfiltrated fraction from the runoff (i.e. the outflow from the fast reservoir) !> is the inflow of the pond reservoir. !- IF (doponds) THEN ! If used, the slope coef is not used in hydrol for water2infilt DO ib=1,nbasmax DO ig=1,nbpt pond_inflow(ig,ib) = fast_flow(ig,ib) * reinf_slope(ig) pond_drainage(ig,ib) = MIN(pond_reservoir(ig)*routing_area(ig,ib)/totarea(ig), & & pond_frac(ig)*routing_area(ig,ib)*k_litt(ig)*dt_routing/one_day) fast_flow(ig,ib) = fast_flow(ig,ib) - pond_inflow(ig,ib) ENDDO ENDDO ELSE DO ib=1,nbasmax DO ig=1,nbpt pond_inflow(ig,ib) = zero pond_drainage(ig,ib) = zero pond_reservoir(ig) = zero ENDDO ENDDO ENDIF source(:,:) = fast_flow(:,:) + slow_flow(:,:) + stream_flow(:,:) CALL downstreamsum(nbpt, nbasmax, source, transport) source(:,:) = fast_flow(:,:)*fast_temp(:,:) + slow_flow(:,:)*slow_temp(:,:) + & & stream_flow(:,:)*stream_temp(:,:) CALL downstreamsum(nbpt, nbasmax, source, transport_temp) !- !- Do the floodings - First initialize !- return_swamp(:,:)=zero floods(:,:)=zero !- !> Over swamp areas, a fraction of water (return_swamp) is withdrawn from the river depending on the !> parameter swamp_cst. !> It will be transferred into soil moisture and thus does not return directly to the river. ! !- 1. Swamps: Take out water from the river to put it to the swamps !- ! IF ( doswamps ) THEN tobeflooded(:) = swamp(:) DO ib=1,nbasmax DO ig=1,nbpt potflood(ig,ib) = transport(ig,ib) ! IF ( tobeflooded(ig) > 0. .AND. potflood(ig,ib) > 0. .AND. floodtemp(ig) > tp_00 ) THEN ! IF (routing_area(ig,ib) > tobeflooded(ig)) THEN floodindex = tobeflooded(ig) / routing_area(ig,ib) ELSE floodindex = 1.0 ENDIF return_swamp(ig,ib) = swamp_cst * potflood(ig,ib) * floodindex ! tobeflooded(ig) = tobeflooded(ig) - routing_area(ig,ib) ! ENDIF ENDDO ENDDO ENDIF !- !- 2. Floodplains: Update the reservoir with the flux computed above. !- IF ( do_floodplains ) THEN DO ig=1,nbpt DO ib=1,nbasmax IF (floodplains(ig, ib) .GT. min_sechiba .AND. floodtemp(ig) .GT. tp_00) THEN floods(ig,ib) = transport(ig,ib) - return_swamp(ig,ib) ENDIF ENDDO ENDDO ENDIF ! ! 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. ! negslow = 0 DO ig=1,nbpt DO ib=1,nbasmax ! fast_reservoir(ig,ib) = fast_reservoir(ig,ib) + runoff(ig)*routing_area(ig,ib) - & & fast_flow(ig,ib) - pond_inflow(ig,ib) ! slow_reservoir(ig,ib) = slow_reservoir(ig,ib) + drainage(ig)*routing_area(ig,ib) - & & slow_flow(ig,ib) ! oldstream = stream_reservoir(ig, ib) * stream_temp(ig,ib) ! stream_reservoir(ig,ib) = stream_reservoir(ig,ib) + flood_flow(ig,ib) + transport(ig,ib) - & & stream_flow(ig,ib) - return_swamp(ig,ib) - floods(ig,ib) ! ! Diagnostics of the stream reservoir ! IF ( routing_area(ig,ib) > zero ) THEN ! 1000 to transform kg into m^3 htmp = stream_reservoir(ig,ib)*1000/routing_area(ig,ib) ewh(ig,ib) = 1.0/(1.0+htmp*hscale) ELSE ewh(ig,ib) = 1.0 ENDIF ! !reste du calcul ! krelax = ewh(ig,ib) ! den = 1.0/(1.0+dt_routing*krelax) IF ( stream_reservoir(ig,ib) > 1.e-6 ) THEN oldtemp = stream_temp(ig,ib) stream_temp(ig,ib) = den * dt_routing * krelax * fast_temp(ig,ib) + & & den * oldstream/stream_reservoir(ig,ib) + & & den * transport_temp(ig, ib)/stream_reservoir(ig,ib) - & & den * oldtemp*stream_flow(ig,ib)/stream_reservoir(ig,ib) ! !Stream_temp [K], stream_reservoir [kg], WaterCp [J/g/K] yields tendencies in GJ/s ! stemp_total_tend(ig,ib) = WaterCp*1.e-6*(stream_temp(ig,ib)*stream_reservoir(ig,ib) - oldstream)/dt_routing stemp_advec_tend(ig,ib) = WaterCp*1.e-6*(transport_temp(ig, ib) - oldtemp*stream_flow(ig,ib))/dt_routing stemp_relax_tend(ig,ib) = WaterCp*1.e-6*stream_reservoir(ig,ib)*krelax*(fast_temp(ig,ib)-stream_temp(ig,ib)) ELSE stream_temp(ig,ib) = MAX(fast_temp(ig,ib), ZeroCelsius) stemp_total_tend(ig,ib) = zero stemp_advec_tend(ig,ib) = zero stemp_relax_tend(ig,ib) = zero ENDIF ! flood_reservoir(ig,ib) = flood_reservoir(ig,ib) + floods(ig,ib) - & & flood_flow(ig,ib) ! pond_reservoir(ig) = pond_reservoir(ig) + pond_inflow(ig,ib) - pond_drainage(ig,ib) ! IF ( flood_reservoir(ig,ib) .LT. zero ) THEN IF ( check_reservoir ) THEN WRITE(numout,*) "WARNING : negative flood reservoir at :", ig, ib, ". Problem is being corrected." WRITE(numout,*) "flood_reservoir, floods, flood_flow : ", flood_reservoir(ig,ib), floods(ig,ib), & & flood_flow(ig,ib) ENDIF stream_reservoir(ig,ib) = stream_reservoir(ig,ib) + flood_reservoir(ig,ib) flood_reservoir(ig,ib) = zero ENDIF ! IF ( stream_reservoir(ig,ib) .LT. zero ) THEN IF ( check_reservoir ) THEN WRITE(numout,*) "WARNING : negative stream reservoir at :", ig, ib, ". Problem is being corrected." WRITE(numout,*) "stream_reservoir, flood_flow, transport : ", stream_reservoir(ig,ib), flood_flow(ig,ib), & & transport(ig,ib) WRITE(numout,*) "stream_flow, return_swamp, floods :", stream_flow(ig,ib), return_swamp(ig,ib), floods(ig,ib) ENDIF fast_reservoir(ig,ib) = fast_reservoir(ig,ib) + stream_reservoir(ig,ib) stream_reservoir(ig,ib) = zero ENDIF ! IF ( fast_reservoir(ig,ib) .LT. zero ) THEN IF ( check_reservoir ) THEN WRITE(numout,*) "WARNING : negative fast reservoir at :", ig, ib, ". Problem is being corrected." WRITE(numout,*) "fast_reservoir, runoff, fast_flow, ponf_inflow : ", fast_reservoir(ig,ib), & &runoff(ig), fast_flow(ig,ib), pond_inflow(ig,ib) ENDIF slow_reservoir(ig,ib) = slow_reservoir(ig,ib) + fast_reservoir(ig,ib) fast_reservoir(ig,ib) = zero ENDIF IF ( slow_reservoir(ig,ib) .LT. - min_sechiba ) THEN IF ( negslow < 20 ) THEN negslow = negslow + 1 negig(negslow) = ig negib(negslow) = ib ENDIF ENDIF ENDDO ENDDO IF ( negslow > 0 ) THEN DO ier = 1,negslow ig = negig(ier) ib = negib(ier) WRITE(numout,*) 'WARNING : There is a negative reservoir at :', ig, ib,lalo(ig,:) WRITE(numout,*) 'WARNING : slowr, slow_flow, drainage', & & slow_reservoir(ig,ib), slow_flow(ig,ib), drainage(ig) WRITE(numout,*) 'WARNING : pondr, pond_inflow, pond_drainage', & & pond_reservoir(ig), pond_inflow(ig,ib), pond_drainage(ig,ib) CALL ipslerr_p(2, 'routing_hr_flow', 'WARNING negative slow_reservoir.','','') ENDDO ENDIF totflood(:) = zero DO ig=1,nbpt DO ib=1,nbasmax totflood(ig) = totflood(ig) + flood_reservoir(ig,ib) ENDDO ENDDO ! ! ESTIMATE the flooded fraction ! IF (do_floodplains .OR. doponds) THEN CALL routing_hr_flood(nbpt, flood_frac, totarea, totflood) ELSE flood_frac(:) = zero flood_height(:,:) = zero flood_frac_bas(:,:) = zero ENDIF !! ANTHONY : OVERFLOW !! CALCULATE TRANSFER BETWEEN FLOODPLAINS RESERVOIR IF (do_floodplains .AND. dofloodoverflow) Then ! The overflow is repeated "overflow_repetition" times ! This is in order to have more stability and ! be able to use lower "overflow_tcst". DO ier = 1,overflow_repetition CALL routing_hr_overflow(nbpt, nbasmax) END DO ! Once done we update the floodplains fraction and the floodplains height CALL routing_hr_flood(nbpt, flood_frac, totarea, totflood) END IF !- !- Compute the total reinfiltration and returnflow to the grid box !> A term of returnflow is computed including the water from the swamps that does not return directly to the river !> but will be put into soil moisture (see hydrol module). !> A term of reinfiltration is computed including the water that reinfiltrated from the ponds and floodplains areas. !> It will be put into soil moisture (see hydrol module). !- IF (do_floodplains .OR. doswamps .OR. doponds) THEN returnflow(:) = zero reinfiltration(:) = zero ! DO ib=1,nbasmax DO ig=1,nbpt returnflow(ig) = returnflow(ig) + return_swamp(ig,ib) reinfiltration(ig) = reinfiltration(ig) + pond_drainage(ig,ib) + flood_drainage(ig,ib) ENDDO ENDDO ! DO ig=1,nbpt returnflow(ig) = returnflow(ig)/totarea(ig) reinfiltration(ig) = reinfiltration(ig)/totarea(ig) ENDDO ELSE returnflow(:) = zero reinfiltration(:) = zero ENDIF ! ! Compute the net irrigation requirement from Univ of Kassel ! ! This is a very low priority process and thus only applies if ! there is some water left in the reservoirs after all other things. ! !> The computation of the irrigation is performed here. !> * First step !> In a first time, the water requirements (irrig_netereq) by the crops for their optimal growth are calculated !> over each irrigated fraction (irrigated(ig)/totarea(ig)). It is the difference !> between the maximal water loss by the crops (transpot_mean) and the net water amount kept by the soil !> (precipitation and reinfiltration). Transpot_mean is computed in the routines enerbil and diffuco. It !> is derived from the effective transpiration parametrization under stress-free conditions, called potential transpiration. !> Crop_coef was used by a previous parametrization of irrigation in the code. Here, its value is equal to one. !> The crop coefficient was constant in space and time to represent a mean resistance of the vegetation to the potential evaporation. !> Now, the term crop_coef*Epot is substituted by transpot_mean (see Guimberteau et al., 2011). !> * Second step !> We compute irrigation needs in order to supply Irrig_netereq. Water for irrigation (irrig_actual) is withdrawn !> from the reservoirs. The amount of water is withdrawn in priority from the stream reservoir. !> If the irrigation requirement is higher than the water availability of the reservoir, water is withdrawn !> from the fast reservoir or, in the extreme case, from the slow reservoir. !> * Third step !> We compute a deficit in water for irrigation. If it is positive, irrigation (depending on water availibility in the reservoirs) !> has not supplied the crops requirements. ! IF ( do_irrigation ) THEN DO ig=1,nbpt ! IF ((vegtot(ig) .GT. min_sechiba) .AND. (humrel(ig) .LT. un-min_sechiba) .AND. & & (runoff(ig) .LT. min_sechiba) ) THEN irrig_netereq(ig) = (irrigated(ig) / totarea(ig) ) * MAX(zero, transpot_mean(ig) - & & (precip(ig)+reinfiltration(ig)) ) ENDIF ! DO ib=1,nbasmax IF ( routing_area(ig,ib) .GT. 0 ) THEN irrig_needs(ig,ib) = irrig_netereq(ig) * routing_area(ig,ib) irrig_actual(ig,ib) = MIN(irrig_needs(ig,ib),& & stream_reservoir(ig,ib) + fast_reservoir(ig,ib) + slow_reservoir(ig,ib) ) slow_reservoir(ig,ib) = MAX(zero, slow_reservoir(ig,ib) + & & MIN(zero, fast_reservoir(ig,ib) + MIN(zero, stream_reservoir(ig,ib)-irrig_actual(ig,ib)))) fast_reservoir(ig,ib) = MAX( zero, & & fast_reservoir(ig,ib) + MIN(zero, stream_reservoir(ig,ib)-irrig_actual(ig,ib))) stream_reservoir(ig,ib) = MAX(zero, stream_reservoir(ig,ib)-irrig_actual(ig,ib) ) irrig_deficit(ig,ib) = irrig_needs(ig,ib)-irrig_actual(ig,ib) ENDIF ENDDO ! ! Check if we cannot find the missing water in another basin of the same grid (stream reservoir only). ! If we find that then we create some adduction from that subbasin to the one where we need it for ! irrigation. ! !> If crops water requirements have not been supplied (irrig_deficit>0), we check if we cannot find the missing water !> in another basin of the same grid. If there is water in the stream reservoir of this subbasin, we create some adduction !> from that subbasin to the one where we need it for irrigation. !> DO ib=1,nbasmax stream_tot = SUM(stream_reservoir(ig,:)) DO WHILE ( irrig_deficit(ig,ib) > min_sechiba .AND. stream_tot > min_sechiba) fi = MAXLOC(stream_reservoir(ig,:)) ib2 = fi(1) irrig_adduct(ig,ib) = MIN(irrig_deficit(ig,ib), stream_reservoir(ig,ib2)) stream_reservoir(ig,ib2) = stream_reservoir(ig,ib2)-irrig_adduct(ig,ib) irrig_deficit(ig,ib) = irrig_deficit(ig,ib)-irrig_adduct(ig,ib) stream_tot = SUM(stream_reservoir(ig,:)) ENDDO ENDDO ! ENDDO ! ! If we are at higher resolution we might need to look at neighboring grid boxes to find the streams ! which can feed irrigation ! !> At higher resolution (grid box smaller than 100x100km), we can import water from neighboring grid boxes !> to the one where we need it for irrigation. ! IF (is_root_prc) THEN ALLOCATE(irrig_deficit_glo(nbp_glo, nbasmax), stream_reservoir_glo(nbp_glo, nbasmax), & & irrig_adduct_glo(nbp_glo, nbasmax), stat=ier) ELSE ALLOCATE(irrig_deficit_glo(0, 0), stream_reservoir_glo(0, 0), & & irrig_adduct_glo(0, 0), stat=ier) ENDIF IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_flow','Pb in allocate for irrig_deficit_glo, stream_reservoir_glo,...','','') CALL gather(irrig_deficit, irrig_deficit_glo) CALL gather(stream_reservoir, stream_reservoir_glo) CALL gather(irrig_adduct, irrig_adduct_glo) IF (is_root_prc) THEN ! DO ig=1,nbp_glo ! Only work if the grid box is smaller than 100x100km. Else the piplines we build ! here would be too long to be reasonable. IF ( resolution_g(ig,1) < 100000. .AND. resolution_g(ig,2) < 100000. ) THEN DO ib=1,nbasmax ! IF ( irrig_deficit_glo(ig,ib) > min_sechiba ) THEN ! streams_around(:,:) = zero ! DO in=1,NbNeighb ig2 = neighbours_g(ig,in) IF (ig2 .GT. 0 ) THEN streams_around(in,:) = stream_reservoir_glo(ig2,:) igrd(in) = ig2 ENDIF ENDDO ! IF ( MAXVAL(streams_around) .GT. zero ) THEN ! ff=MAXLOC(streams_around) ig2=igrd(ff(1)) ib2=ff(2) ! IF ( routing_area_glo(ig2,ib2) .GT. 0 .AND. stream_reservoir_glo(ig2,ib2) > zero ) THEN adduction = MIN(irrig_deficit_glo(ig,ib), stream_reservoir_glo(ig2,ib2)) stream_reservoir_glo(ig2,ib2) = stream_reservoir_glo(ig2,ib2) - adduction irrig_deficit_glo(ig,ib) = irrig_deficit_glo(ig,ib) - adduction irrig_adduct_glo(ig,ib) = irrig_adduct_glo(ig,ib) + adduction ENDIF ! ENDIF ! ENDIF ! ENDDO ENDIF ENDDO ! ENDIF ! CALL scatter(irrig_deficit_glo, irrig_deficit) CALL scatter(stream_reservoir_glo, stream_reservoir) CALL scatter(irrig_adduct_glo, irrig_adduct) DEALLOCATE(irrig_deficit_glo, stream_reservoir_glo, irrig_adduct_glo) ENDIF !! Calculate the net water flow to each routing reservoir (in kg/dt) !! to further diagnose the corresponding water budget residu !! in routing_highres_main netflow_fast_diag(:) = zero netflow_slow_diag(:) = zero netflow_stream_diag(:) = zero DO ib=1,nbasmax DO ig=1,nbpt netflow_fast_diag(ig) = netflow_fast_diag(ig) + runoff(ig)*routing_area(ig,ib) & - fast_flow(ig,ib) - pond_inflow(ig,ib) netflow_slow_diag(ig) = netflow_slow_diag(ig) + drainage(ig)*routing_area(ig,ib) & - slow_flow(ig,ib) netflow_stream_diag(ig) = netflow_stream_diag(ig) + flood_flow(ig,ib) + transport(ig,ib) & - stream_flow(ig,ib) - return_swamp(ig,ib) - floods(ig,ib) ENDDO ENDDO !! Grid cell averaging DO ig=1,nbpt netflow_fast_diag(ig) = netflow_fast_diag(ig)/totarea(ig) netflow_slow_diag(ig) = netflow_slow_diag(ig)/totarea(ig) netflow_stream_diag(ig) = netflow_stream_diag(ig)/totarea(ig) ENDDO ! ! ! Compute the fluxes which leave the routing scheme ! ! Lakeinflow is in Kg/dt ! returnflow is in Kg/m^2/dt ! hydrographs(:) = zero hydrotemp(:) = zero HTUhgmon(:,:) = zero HTUtempmon(:,:) = zero slowflow_diag(:) = zero fast_diag(:) = zero slow_diag(:) = zero stream_diag(:) = zero flood_diag(:) = zero pond_diag(:) = zero irrigation(:) = zero ! ! DO ib=1,nbasmax ! DO ig=1,nbpt ! DO im=1,nbasmon IF (HTUdiag_loc(ig,im) > 0 .AND. HTUdiag_loc(ig,im) .EQ. ib ) THEN HTUhgmon(ig,im) = fast_flow(ig,ib) + slow_flow(ig,ib) + stream_flow(ig,ib) HTUtempmon(ig,im) = stream_temp(ig,ib) ENDIF ENDDO ! IF (hydrodiag(ig) == ib) THEN hydrographs(ig) = fast_flow(ig,ib) + slow_flow(ig,ib) + stream_flow(ig,ib) hydrotemp(ig) = stream_temp(ig,ib) slowflow_diag(ig) = slowflow_diag(ig) + slow_flow(ig,ib) ENDIF 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) irrigation (ig) = irrigation (ig) + irrig_actual(ig,ib) + irrig_adduct(ig,ib) ENDDO ENDDO ! DO ig=1,nbpt fast_diag(ig) = fast_diag(ig)/totarea(ig) slow_diag(ig) = slow_diag(ig)/totarea(ig) stream_diag(ig) = stream_diag(ig)/totarea(ig) flood_diag(ig) = flood_diag(ig)/totarea(ig) pond_diag(ig) = pond_reservoir(ig)/totarea(ig) ! irrigation(ig) = irrigation(ig)/totarea(ig) ! ! The three output types for the routing : endoheric basins,, rivers and ! diffuse coastal flow. ! lakeinflow(ig) = transport(ig,nbasmax+1) coastalflow(ig) = transport(ig,nbasmax+2) riverflow(ig) = transport(ig,nbasmax+3) ! ENDDO ! flood_res = flood_diag + pond_diag !! 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 DO ig=1,nbpt lake_overflow(ig) = MAX(0., lake_reservoir(ig) - max_lake_reservoir*totarea(ig)) lake_reservoir(ig) = lake_reservoir(ig) - lake_overflow(ig) END DO ! Transform lake_overflow from kg/grid-cell/dt_routing into kg/m^2/s CALL xios_orchidee_send_field("lake_overflow",lake_overflow(:)/totarea(:)/dt_routing) ! Calculate the sum of the lake_overflow and distribute it uniformly over all gridboxes CALL gather(lake_overflow,lake_overflow_g) IF (is_root_prc) THEN total_lake_overflow=SUM(lake_overflow_g) END IF CALL bcast(total_lake_overflow) ! Distribute the lake_overflow uniformly over all coastal gridcells ! lake_overflow_coast is only calculated to be used as diagnostics if needed DO ig=1,nbpt coastalflow(ig) = coastalflow(ig) + total_lake_overflow/nb_coast_gridcells * mask_coast(ig) lake_overflow_coast(ig) = total_lake_overflow/nb_coast_gridcells * mask_coast(ig) END DO ! Transform from kg/grid-cell/dt_routing into m^3/grid-cell/s to match output unit of coastalflow CALL xios_orchidee_send_field("lake_overflow_coast",lake_overflow_coast/mille/dt_routing) END SUBROUTINE routing_hr_flow ! !! ================================================================================================================================ !! SUBROUTINE : groundwatertemp !! !>\BRIEF : This subroutine computes the temperature of the groundwater leaving the HTU !! !! 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 groundwatertemp(nbpt, nbasmax, nl, tempdiag, lev, dlz, fast_temp, slow_temp) ! INPUT INTEGER(i_std), INTENT(in) :: nbpt, nbasmax, nl REAL(r_std), INTENT(in) :: tempdiag(nbpt,nl) REAL(r_std), INTENT(in) :: lev(nl), dlz(nl) REAL(r_std), INTENT(inout) :: slow_temp(nbpt,nbasmax), fast_temp(nbpt,nbasmax) ! OUTPUT ! LOCAL INTEGER(i_std) :: ig, ib, im REAL(r_std) :: sw REAL(r_std) :: rw(nl), dw(nl) LOGICAL, SAVE :: alltop=.FALSE. LOGICAL, SAVE :: FirstCall=.TRUE. ! IF ( FirstCall ) THEN !Config Key = ROUTING_ALLTOPT !Config Desc = Should drainage have the temperature of the top soil (0.3m) ? !Config Def = False !Config Help = The default behaviour of the scheme is that runoff has the temperature !Config Help of the top 30 cm of soil. Drainage will have the temperature of the lowest !Config Help soil layer (3-17m). If set to True this flag will give drainage the same !Config Help temperature as runoff. !Config Units = Logical alltop=.FALSE. CALL getin_p('ROUTING_ALLTOPT', alltop) ! WRITE(numout,*) "Runoff will have the average soil temperature of layers from ", runofftempdepth(1),& & " to ", runofftempdepth(2), "[m]" ! IF ( alltop ) THEN WRITE(numout,*) "Drainage will have the average soil temperature of layers from ", runofftempdepth(1),& & " to ", runofftempdepth(2), "[m]" ELSE WRITE(numout,*) "Drainage will have the average soil temperature of layers from ", drainagetempdepth(1),& & " to ", MIN(drainagetempdepth(2), SUM(dlz)), "[m]" ENDIF FirstCall=.FALSE. ENDIF ! CALL tempdepthweight(nl, dlz, runofftempdepth(1), runofftempdepth(2), rw) CALL tempdepthweight(nl, dlz, drainagetempdepth(1), MIN(drainagetempdepth(2), SUM(dlz)), dw) ! slow_temp(:,:) = zero fast_temp(:,:) = zero ! Compute for each HTU the temperature of runoff and drainage water. DO im = 1,nl DO ib=1,nbasmax DO ig=1,nbpt fast_temp(ig,ib) = fast_temp(ig,ib) + tempdiag(ig,im)*rw(im) ! The option to have drainage water at the same temperature as runoff IF ( alltop ) THEN slow_temp(ig,ib) = slow_temp(ig,ib) + tempdiag(ig,im)*rw(im) ELSE slow_temp(ig,ib) = slow_temp(ig,ib) + tempdiag(ig,im)*dw(im) ENDIF ENDDO ENDDO ENDDO END SUBROUTINE groundwatertemp SUBROUTINE tempdepthweight(n, dz, top, bot, w) ! Input INTEGER(i_std), INTENT(in) :: n REAL(r_std), INTENT(in) :: dz(n) REAL(r_std), INTENT(in) :: top, bot ! Output REAL(r_std), INTENT(out) :: w(n) ! Local INTEGER(i_std) :: i REAL(r_std) :: sw w(:) = zero sw = zero DO i=1,n w(i) = MAX(zero, MIN(sw+dz(i), bot) - MAX(top, sw)) sw = sw + dz(i) ENDDO w(:) = w(:)/(bot-top) END SUBROUTINE tempdepthweight !! ================================================================================================================================ !! SUBROUTINE : downstreamsum !! !>\BRIEF : This subroutine sums the input variables onto the downstream HTU in the river graph. !! !! DESCRIPTION : We assume that the downstream HTU is defined by route_togrid and route_tobas. As these !! donwstream HTU can be on another processor we do this job on the root processor. So before we need to !! transfer all the data onto that processor and then redistribute the result. !! Keep in mind that if an HTU does not exit then route_tobas = 0. So the result array needs !! to have this index. The end of the rivers are between nbmax+1 and nbmax+3 so this indexing space is also !! needed in the result array. !! !! RECENT CHANGE(S): None !! !! MAIN OUTPUT VARIABLE(S): !! !! REFERENCES : None !! !! FLOWCHART :None !! \n !_ ================================================================================================================================ !- SUBROUTINE downstreamsum(nbpt, nbmax, v, t) ! Input INTEGER(i_std), INTENT(in) :: nbpt, nbmax REAL(r_std), INTENT(in), DIMENSION(nbpt, nbmax) :: v ! Output REAL(r_std), INTENT(out), DIMENSION(nbpt, 0:nbmax+3) :: t ! ! Local ! INTEGER(i_std) :: ig, ib, rtg, rtb INTEGER(i_std) :: ier REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: v_g, t_g ! ! Allocate memory if needed. Should only happen only once in order to reduce computing time. ! IF ( .NOT. ALLOCATED(v_g) ) THEN IF (is_root_prc) THEN ALLOCATE(v_g(nbp_glo,nbmax), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'downstreamsum','Pb in allocate for v_g','','') ELSE ALLOCATE(v_g(1,1)) ENDIF ENDIF IF ( .NOT. ALLOCATED(t_g) ) THEN IF (is_root_prc) THEN ALLOCATE(t_g(nbp_glo,0:nbmax+3), stat=ier) IF (ier /= 0) CALL ipslerr_p(3,'downstreamsum','Pb in allocate for t_g','','') ELSE ALLOCATE(t_g(1,1)) ENDIF ENDIF ! ! Gather the source variable on the root processor. ! CALL gather(v, v_g) ! ! The downstream sum is performed only on the root processor. ! IF (is_root_prc) THEN t_g(:,:) = zero DO ib=1,nbmax DO ig=1,nbp_glo rtg = route_togrid_glo(ig,ib) rtb = route_tobasin_glo(ig,ib) t_g(rtg,rtb) = t_g(rtg,rtb) + v_g(ig,ib) ENDDO ENDDO ENDIF ! ! Redistribute the downstream field to the all processors. ! CALL scatter(t_g, t) ! END SUBROUTINE downstreamsum !! ================================================================================================================================ !! SUBROUTINE : routing_hr_flood !! !>\BRIEF : This subroutine estimate the flood fraction and the flood height for each HTU !! !! 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_hr_flood(nbpt, flood_frac, totarea, totflood) ! IMPLICIT NONE ! !! INPUT VARIABLES INTEGER(i_std), INTENT(in) :: nbpt !! Domain size (unitless) REAL(r_std), INTENT(in), DIMENSION(nbpt) :: totflood !! Total amount of water in the floodplains reservoir (kg) REAL(r_std), INTENT(in), DIMENSION(nbpt) :: totarea !! Total area of basin (m^2) !! Flooded fraction of the grid box (unitless;0-1) ! !! OUTPUT VARIABLES REAL(r_std), INTENT(inout) :: flood_frac(nbpt) ! !! LOCAL VARIABLES INTEGER(i_std) :: ig, ib !! Indices (unitless) REAL(r_std) :: diff, voltemp !! Discharge reduction due to floodplains !_ ================================================================================================================================ ! ! ! Initialize the variables flood_frac(:) = zero flood_height(:,:) = zero flood_frac_bas(:,:) = zero DO ig=1, nbpt IF (totflood(ig) .GT. min_sechiba) THEN DO ib=1,nbasmax IF (floodplains(ig,ib) .GT. min_sechiba) THEN ! We have to convert h0 to m and the flood_reservoir in m^3 flood_frac_bas(ig,ib) = ((fp_beta(ig,ib)+un) * flood_reservoir(ig,ib) / 1000) / ( floodcri(ig,ib) / 1000 * floodplains(ig,ib)) flood_frac_bas(ig,ib) = (flood_frac_bas(ig,ib)) ** (fp_beta(ig,ib)/(fp_beta(ig,ib)+1)) flood_frac_bas(ig,ib) = MIN(flood_frac_bas(ig,ib), floodplains(ig,ib)/ routing_area(ig,ib) ) ! flood_height is in mm ! there is two cases: flood_height < h0, flood_height >= h0 (this corresponds to flood_frac_bas = 1 ) IF ( flood_frac_bas(ig,ib) .EQ. floodplains(ig,ib) / routing_area(ig,ib) ) THEN ! voltemp is on m^3 ! Calculation of volume corresponding to h0 voltemp = floodplains(ig,ib)/(fp_beta(ig,ib)+un) * ( floodcri(ig,ib) / 1000 ) voltemp = flood_reservoir(ig,ib) / 1000 - voltemp ! flood height is in mm flood_height(ig, ib) = voltemp / floodplains(ig,ib) * 1000 + floodcri(ig,ib) ELSE ! flood height is in mm flood_height(ig, ib) = (flood_frac_bas(ig,ib)) ** (1/fp_beta(ig,ib)) * floodcri(ig,ib) END IF ENDIF ENDDO ENDIF DO ib=1,nbasmax flood_frac(ig) = flood_frac(ig) + flood_frac_bas(ig,ib) * routing_area(ig,ib) / totarea(ig) END DO flood_frac(ig) = flood_frac(ig) + pond_frac(ig) ! ENDDO END SUBROUTINE routing_hr_flood ! !! ================================================================================================================================ !! SUBROUTINE : routing_hr_overflow !! !>\BRIEF : This subroutine performs the overflow fluxes !! !! DESCRIPTION (definitions, functional, design, flags): \n !! !! RECENT CHANGE(S): None !! !! MAIN OUTPUT VARIABLE(S): !! !! REFERENCES : None !! !! FLOWCHART :None !! \n !_ ================================================================================================================================ !- SUBROUTINE routing_hr_overflow(nbpt, nbasmax) ! IMPLICIT NONE ! !! INPUT VARIABLES INTEGER(i_std), INTENT(in) :: nbpt,nbasmax !! Domain size (unitless) ! !! LOCAL VARIABLES REAL(r_std), DIMENSION(nbpt,nbasmax) :: transport_overflow !! Water transport between floodplains - flood overflow (kg/dt) REAL(r_std), DIMENSION(nbp_glo,nbasmax) :: transport_overflow_glo !! Water transport between floodplains - flood overflow (kg/dt) REAL(r_std), DIMENSION(nbpt,nbasmax) :: overflow_loss !! Water loss from flood overflow (kg/dt) REAL(r_std), DIMENSION(nbp_glo,nbasmax) :: overflow_loss_glo !! Water loss from flood overflow (kg/dt) ! REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: flood_height_g !! Floodplains height (m) REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: flood_frac_bas_g !! Fraction of the HTU flooded REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: flood_reservoir_glo !! Water amount in the stream reservoir (kg) ! REAL(r_std), ALLOCATABLE, DIMENSION(:) :: DH,DH_temp !! Difference of height - flood overflow (kg/dt) ! INTEGER(i_std) :: numflood !! ! INTEGER(i_std) :: ig, ib, inf,inb,ing !! Indices (unitless) REAL(r_std) :: diff !! Discharge reduction due to floodplains REAL(r_std) :: flow !! Outflow computation for the reservoirs (kg/dt) REAL(r_std) :: dorog !! Discharge reduction due to floodplains INTEGER(i_std) :: ier !! Error handling !_ ================================================================================================================================ ! !! ANTHONY : OVERFLOW !! CALCULATE TRANSFER BETWEEN FLOODPLAINS RESERVOIR IF (is_root_prc) THEN ALLOCATE( flood_height_g(nbp_glo, nbasmax), flood_frac_bas_g(nbp_glo, nbasmax), stat=ier) ALLOCATE( flood_reservoir_glo(nbp_glo, nbasmax), stat=ier) ELSE ALLOCATE( flood_height_g(1,1), flood_frac_bas_g(1,1), stat=ier) ALLOCATE( flood_reservoir_glo(1, 1), stat=ier) ENDIF ! IF (ier /= 0) CALL ipslerr_p(3,'routing_hr_flow','Pb in allocate for flood_height_glo/floog_frac_glo','','') ! CALL gather(flood_height,flood_height_g) CALL gather(flood_frac_bas,flood_frac_bas_g) CALL gather(flood_reservoir,flood_reservoir_glo) ! IF (is_root_prc) THEN transport_overflow_glo(:,:) = 0 overflow_loss_glo(:,:) = 0 DO ib=1,nbasmax DO ig=1,nbp_glo IF ( floodplains_glo(ig,ib)/routing_area_glo(ig,ib) .GT. 0.5) THEN numflood = 0 ! Number of inflows for overflow ALLOCATE(DH(route_innum_glo(ig,ib))) DH(:) = 0 DH_temp(:) = -1 DO inf=1,route_innum_glo(ig,ib) ing = route_ingrid_glo(ig,ib,inf) inb = route_inbasin_glo(ig,ib,inf) IF ( floodplains_glo(ing,inb)/routing_area_glo(ing,inb) .GT. 0 ) THEN ! Minimum of deltaorog is defined at lim_floodcri (0.3 m ! can be used). dorog = MAX(orog_min_glo(ing,inb)- orog_min_glo(ig,ib), lim_floodcri) ! flood_height is in mm and orog min in m diff = (flood_height_g(ig,ib)- flood_height_g(ing,inb))/1000 - dorog DH(inf) = max(diff, 0.) ! ! Flux is estimated via floodplains_glo ! Then factor 1000 is to convert m^3 to kg ! OVERFLOW_TCST is in seconds flow = DH(inf) * (floodplains_glo(ig,ib)* floodplains_glo(ing,inb))/(floodplains_glo(ig,ib)+floodplains_glo(ing,inb))*1000 / overflow_tcst * dt_routing / one_day transport_overflow_glo(ing,inb) = transport_overflow_glo(ing,inb) + flow overflow_loss_glo(ig,ib) = overflow_loss_glo(ig,ib) + flow END IF END DO DEALLOCATE(DH) END IF ENDDO ENDDO END IF ! Send to local variables CALL scatter(transport_overflow_glo, transport_overflow) CALL scatter(overflow_loss_glo, overflow_loss) ! Apply the volume changes DO ig=1,nbpt DO ib=1,nbasmax IF ( floodplains(ig,ib) .GT. 0 ) THEN flood_reservoir(ig,ib) = flood_reservoir(ig,ib) + transport_overflow(ig,ib) - overflow_loss(ig,ib) ! NEED to check if flood reservoir is less than 0, this may be a critical issue ! Solved by an adequate use of an higher overflow time constant ! To obtain the same result as with a lower overflow parameter ! -> repeat a few time the operation with and higher overflow parameter IF ( flood_reservoir(ig,ib) .LT. 0 ) THEN WRITE(*,*) "Issue of flood reservoir < 0 due to overflow at ", ig, ib stream_reservoir(ig,ib) = stream_reservoir(ig,ib) + stream_reservoir(ig,ib) ! + because negative ! flood_reservoir(ig,ib) = 0 END IF END IF END DO END DO DEALLOCATE( flood_height_g, flood_frac_bas_g) END SUBROUTINE routing_hr_overflow ! !! ================================================================================================================================ !! SUBROUTINE : routing_hr_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_hr_lake(nbpt, dt_routing, lakeinflow, humrel, return_lakes) ! IMPLICIT NONE ! !! 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(out) :: lakeinflow(nbpt) !! Water inflow to the lakes (kg/dt) REAL(r_std), INTENT(in) :: humrel(nbpt) !! Soil moisture stress, root extraction potential (unitless) ! !! OUTPUT VARIABLES REAL(r_std), INTENT(out) :: return_lakes(nbpt) !! Water from lakes flowing back into soil moisture (kg/m^2/dt) ! !! 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 = SUM(routing_area(ig,:)) ! lake_reservoir(ig) = lake_reservoir(ig) + lakeinflow(ig) IF ( doswamps ) THEN ! Calculate a return flow that will be extracted from the lake reservoir and reinserted in the soil in hydrol ! 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 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 better to scale it to the size of the lake ! but this information is not yet available. lake_diag(ig) = lake_reservoir(ig)/total_area lakeinflow(ig) = lakeinflow(ig)/total_area ENDDO ! END SUBROUTINE routing_hr_lake ! !! ================================================================================================================================ !! SUBROUTINE : routing_hr_basins_p !! !>\BRIEF This routing read the file created by RoutingPreProc : https://gitlab.in2p3.fr/ipsl/lmd/intro/routingpp !! !! DESCRIPTION (definitions, functional, design, flags) : None !! Once the atmospheric grid is defined and the land/sea mask set, RoutingPreProc has to used to generate the !! HTU graphs for the domain. This can be done either on the basis of the HydroSHEDS, MERIT or the old Vörösmarty map !! of catchments. During this step all the information will be created to allow ORCHIDEE to route the water and !! and monitor the flows at given stations. !! For the moment the ROUTING_FILE (Perhaps to renamed RoutingGraph) is read using IOIPSL but that should evolve toward XIOS. !! !! RECENT CHANGE(S): None !! !! MAIN OUTPUT VARIABLE(S): !! !! REFERENCES : None !! !! FLOWCHART : None !! \n !_ ================================================================================================================================ SUBROUTINE routing_hr_basins_p(nbpt, lalo, neighbours, resolution, contfrac) ! IMPLICIT NONE ! !! INPUT VARIABLES INTEGER(i_std), INTENT(in) :: nbpt !! Domain size (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,NbNeighb) !! Vector of neighbours for each grid point (1=North and then clockwise) (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) ! ! LOCAL ! INTEGER(i_std) :: iml, jml, lml, tml INTEGER(i_std) :: i, j, ni, fid, ib, ig, ic, ign, ibn, og, ob, ier, im REAL(r_std) :: corr REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:) :: tmpvar_glo REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: tmpvar REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: lon, lat, landindex INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:) :: indextab INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:) :: landfileindex INTEGER(i_std), ALLOCATABLE, DIMENSION(:) :: land2land INTEGER(i_std) :: nbhtumon ! !_ ================================================================================================================================ ! ! ! IF (is_root_prc) THEN ! CALL flininfo(graphfilename, iml, jml, lml, tml, fid) ! IF (iml .NE. iim_g .AND. jml .NE. jjm_g ) THEN CALL ipslerr(3,'routing_hr_basins_p',& 'The routing graph file does not have the right dimensions for the model.', & 'Are you sure you are using the right routing graph file ?', ' ') ENDIF ! ! ALLOCATE(tmpvar_glo(iml,jml,nbasmax)) ALLOCATE(tmpvar(iml,jml)) ALLOCATE(lon(iml,jml)) ALLOCATE(lat(iml,jml)) ALLOCATE(landindex(iml,jml)) ALLOCATE(indextab(iml,jml)) ALLOCATE(landfileindex(iml,jml)) ALLOCATE(land2land(iml*jml)) ! CALL flinget(fid, 'lon', iml, jml, 1, tml, 1, 0, lon) CALL flinget(fid, 'lat', iml, jml, 1, tml, 1, 0, lat) CALL flinget(fid, 'nbpt_glo', iml, jml, 1, tml, 1, 0, landindex) ! ! Replace NaN and other undef values ! DO i=1,iml DO j=1,jml IF ( landindex(i,j) /= landindex(i,j) .OR. landindex(i,j) >= undef_graphfile) THEN landindex(i,j) = -1 ENDIF ENDDO ENDDO ! ! Compute land index for file data. Information could be in file ! ! ni=NINT(MAXVAL(landindex)) IF ( ni .NE. nbp_glo) THEN WRITE(numout,*) "Error routing_hr_basins_p : ni, nbp_glo : ", ni, nbp_glo, undef_graphfile CALL ipslerr(3,'routing_hr_basins_p',& 'The routing graph file does not have the same number', & 'of land points as the model.',& ' ') ENDIF ! CALL routing_hr_indexfilegrid(iml, jml, nbp_glo, lon, lat, landindex, indextab, land2land) ! CALL flinget(fid, 'basin_area', iml, jml, nbasmax, tml, 1, 0, tmpvar_glo) CALL routing_hr_landgather(iml, jml, nbasmax, nbp_glo, indextab, tmpvar_glo, routing_area_glo, & & zero) IF ( do_floodplains ) THEN CALL flinget(fid, 'basin_floodp', iml, jml, nbasmax, tml, 1, 0, tmpvar_glo) CALL routing_hr_landgather(iml, jml, nbasmax, nbp_glo, indextab, tmpvar_glo, floodplains_glo, & & zero) ! CALL flinget(fid, 'floodcri', iml, jml, nbasmax, tml, 1, 0, tmpvar_glo) CALL routing_hr_landgather(iml, jml, nbasmax, nbp_glo, indextab, tmpvar_glo, floodcri_glo, & & un) ! CALL flinget(fid, 'basin_beta_fp', iml, jml, nbasmax, tml, 1, 0, tmpvar_glo) CALL routing_hr_landgather(iml, jml, nbasmax, nbp_glo, indextab, tmpvar_glo, fp_beta_glo, & & un) END IF CALL flinget(fid, 'topoindex', iml, jml, nbasmax, tml, 1, 0, tmpvar_glo) CALL routing_hr_landgather(iml, jml, nbasmax, nbp_glo, indextab, tmpvar_glo, topo_resid_glo, & & undef_graphfile) IF ( graphfile_version >= 2.0) THEN CALL flinget(fid, 'topoindex_stream', iml, jml, nbasmax, tml, 1, 0, tmpvar_glo) CALL routing_hr_landgather(iml, jml, nbasmax, nbp_glo, indextab, tmpvar_glo, stream_resid_glo, & & undef_graphfile) CALL ipslerr(1,'routing_hr_basins_p',& 'The topoindex_stream variable was found in routing_graph.nc', & 'It will be used the topographic index of the stream store.',& ' ') ELSE stream_resid_glo(:,:) = topo_resid_glo(:,:) ENDIF stream_maxresid=MAXVAL(stream_resid_glo, MASK=stream_resid_glo .LT. undef_graphfile) CALL flinget(fid, 'basinid', iml, jml, nbasmax, tml, 1, 0, tmpvar_glo) CALL routing_hr_landgather(iml, jml, nbasmax, nbp_glo, indextab, tmpvar_glo, global_basinid_glo, & & undef_int) CALL flinget(fid, 'routetogrid', iml, jml, nbasmax, tml, 1, 0, tmpvar_glo) CALL routing_hr_landgather(iml, jml, nbasmax, nbp_glo, indextab, tmpvar_glo, route_togrid_glo, & & undef_int) CALL routing_hr_convertlandpts(nbp_glo, nbasmax, land2land, route_togrid_glo) CALL flinget(fid, 'routetobasin', iml, jml, nbasmax, tml, 1, 0, tmpvar_glo) CALL routing_hr_landgather(iml, jml, nbasmax, nbp_glo, indextab, tmpvar_glo, route_tobasin_glo, 0) CALL flinget(fid, 'routenbintobas', iml, jml, nbasmax, tml, 1, 0, tmpvar_glo) CALL routing_hr_landgather(iml, jml, nbasmax, nbp_glo, indextab, tmpvar_glo, route_nbintobas_glo, 0) !! IF ( dofloodoverflow ) THEN CALL flinget(fid, 'basin_orog_min', iml, jml, nbasmax, tml, 1, 0, tmpvar_glo) CALL routing_hr_landgather(iml, jml, nbasmax, nbp_glo, indextab, tmpvar_glo, orog_min_glo, un) END IF !! IF ( graphfile_version >= 2.6) THEN CALL flinget(fid, 'gridrephtu', iml, jml, 1, tml, 1, 0, tmpvar) CALL routing_hr_landgather(iml, jml, nbp_glo, indextab, tmpvar, hydrodiag_glo, -1) ELSE hydrodiag_glo(:) = 1 ENDIF !! IF ( MonitoringinGraph ) THEN CALL flinget(fid, 'HTUmonitor', iml, jml, nbasmon, tml, 1, 0, tmpvar_glo) CALL routing_hr_landgather(iml, jml, nbasmon, nbp_glo, indextab, tmpvar_glo, HTUdiag_glo, -1) ELSE HTUdiag_glo(:,:) = -1 ENDIF ! CALL flinclo(fid) DEALLOCATE(indextab) DEALLOCATE(lon) DEALLOCATE(lat) DEALLOCATE(tmpvar_glo) DEALLOCATE(tmpvar) ! ! Convert floodplains fraction into floodplains surface IF ( do_floodplains ) THEN !floodplains_glo(:, :) = 0 DO ig = 1,nbp_glo DO ib = 1,nbasmax floodplains_glo(ig, ib) = routing_area_glo(ig,ib) * floodplains_glo(ig, ib) END DO END DO END IF ! ! Verifications of the routing graph. ! nbhtumon = 0 DO ig = 1,nbp_glo ! Noramlize the areas so that differences in precision of area compution by RoutingPP do not affect the model ! corr = contfrac_g(ig)*area_g(ig)/SUM(routing_area_glo(ig,:)) IF (ABS(1 - corr) > 0.0002 ) THEN WRITE(*,*) "Correcting the HTU area to take into account contfrac", corr IF ( ABS(1 - corr) > 0.1) THEN WRITE(*,*) "Coordinates : ", lalo_g(ig,1), lalo_g(ig,2) WRITE(*,*) "Contfrac and area in model : ", contfrac_g(ig), area_g(ig) WRITE(*,*) "Total grid area in graph file : ", SUM(routing_area_glo(ig,:)) WRITE(*,*) "The new areas are : ", SUM(routing_area_glo(ig,:)), contfrac_g(ig)*area_g(ig) WRITE(*,*) "Correction factor : ", corr CALL ipslerr(3,'routing_hr_basins_p',& 'There is a mismatch in the area of the grid', & 'Either there are issues with the projection of the grid ',' or contfrac mismatches.') ELSE CALL ipslerr(2,'routing_hr_basins_p',& 'The area of the grid had to be adjusted by less than 10% :', & ' ',' ') ENDIF ENDIF DO ib = 1,nbasmax routing_area_glo(ig,ib) = corr*routing_area_glo(ig,ib) ENDDO ! ! DO ib = 1,nbasmax ! IF (topo_resid_glo(ig,ib) <= zero .AND. route_tobasin_glo(ig, ib) .LE. nbasmax+3) THEN ! If the basin has no surface we change silently as it does not matter. IF ( routing_area_glo(ig,ib) > zero ) THEN CALL ipslerr(2,'routing_hr_basins_p',& 'Some zero topo_resid (topoindex) values were encoutered and replaced here :', & ' ',' ') WRITE(*,*) "routing_hr_basins_p : topo_resid_glo : ", topo_resid_glo(ig,ib), routing_area_glo(ig,ib) WRITE(*,*) "routing_hr_basins_p : Coordinates : ", lalo_g(ig,1), lalo_g(ig,2) topo_resid_glo(ig,ib) = 10 stream_resid_glo(ig,ib) = 10 WRITE(*,*) "routing_hr_basins_p : New topo_resid_glo : ", topo_resid_glo(ig,ib) ELSE topo_resid_glo(ig,ib) = 10 stream_resid_glo(ig,ib) = 10 ENDIF ENDIF ! ! IF ( route_togrid_glo(ig, ib) > nbp_glo ) THEN IF ( route_tobasin_glo(ig,ib) <= nbasmax+3 ) THEN WRITE(*,*) "Issues with the global grid : ", ig, ib, route_togrid_glo(ig, ib), route_tobasin_glo(ig,ib) CALL ipslerr(3,'routing_hr_basins_p','route_togrid is not compatible with the model configuration', & ' ',' ') ENDIF ELSE ic = 0 ign = ig ibn = ib ! Locate outflow point DO WHILE (ibn .GT. 0 .AND. ibn .LE. nbasmax .AND. ic .LT. nbasmax*nbp_glo) ic = ic + 1 og = ign ob = ibn ign = route_togrid_glo(og, ob) ibn = route_tobasin_glo(og, ob) ! IF (ibn .GT. nbasmax+3 .OR. ign .GT. nbp_glo) THEN WRITE(*,*) "Reached point ", ign, ibn, " on condition ", nbasmax+3, nbp_glo WRITE(*,*) "Why do we flow into basin :", route_tobasin_glo(og, ob), " at ", og,ob WRITE(*,*) "Coordinates : ", lalo_g(ob,1), lalo_g(ob,2) WRITE(*,*) "neighbours_g : ", MINVAL(neighbours_g(ob,:)) CALL ipslerr(3,'routing_hr_basins_p','The river flows into a place outside of the grid.', & ' ',' ') ENDIF ENDDO IF ( ic .GE. nbasmax*nbp_glo) THEN WRITE(*,*) "Some river did not converge on point ", ig, ib, ic WRITE(*,*) "The start point in the graph was : ", lalo_g(ig,2), lalo_g(ig,1), ib WRITE(*,*) "The last point we passed through was : ", lalo_g(og,2), lalo_g(og,1), ob WRITE(*,*) "The next one would be : ", lalo_g(ign,2), lalo_g(ign,1), ibn og = route_togrid_glo(ign, ibn) ob = route_tobasin_glo(ign, ibn) WRITE(*,*) "The after next HTU would be : ", lalo_g(og,2), lalo_g(og,1), ob WRITE(*,*) "Last information : ", ign, ibn CALL ipslerr(3,'routing_hr_basins_p','The river never flows into an outflow point.', & ' ',' ') ENDIF ENDIF ENDDO ! ! Count stations to be monitored ! DO im=1,nbasmon IF ( HTUdiag_glo(ig,im) > 0 ) THEN nbhtumon = nbhtumon + 1 ENDIF ENDDO ENDDO WRITE(numout,*) "Found a total of ", nbhtumon, " HTUs to be monitored and written into HTUhgmon" ! ! Compute num_largest ! num_largest = COUNT(route_tobasin_glo .EQ. nbasmax+3) WRITE(numout,*) "After _basins_p : Number of largest rivers : ", COUNT(route_tobasin_glo .EQ. nbasmax+3) ENDIF ! CALL bcast(num_largest) CALL bcast(nbasmax) CALL bcast(nbasmon) CALL bcast(inflows) ! CALL scatter(routing_area_glo,routing_area_loc) IF ( do_floodplains ) THEN CALL scatter(floodplains_glo,floodplains_loc) CALL scatter(floodcri_glo, floodcri_loc) CALL scatter(fp_beta_glo, fp_beta_loc) END IF CALL scatter(global_basinid_glo, global_basinid_loc) CALL scatter(topo_resid_glo, topo_resid_loc) CALL scatter(stream_resid_glo, stream_resid_loc) CALL scatter(route_togrid_glo, route_togrid_loc) CALL scatter(route_tobasin_glo, route_tobasin_loc) CALL scatter(route_nbintobas_glo, route_nbintobas_loc) CALL scatter(hydrodiag_glo, hydrodiag_loc) CALL scatter(HTUdiag_glo, HTUdiag_loc) IF ( do_floodplains .AND. dofloodoverflow ) THEN CALL scatter(orog_min_glo, orog_min_loc) END IF ! CALL bcast(stream_tcst) CALL bcast(fast_tcst) CALL bcast(slow_tcst) CALL bcast(flood_tcst) CALL bcast(swamp_cst) CALL bcast(lim_floodcri) CALL bcast(stream_maxresid) ! END SUBROUTINE routing_hr_basins_p ! !! ================================================================================================================================ !! SUBROUTINE : routing_hr_graphinfo !! !>\BRIEF Extract some basic information from the routing graph file which cannot be obtained through IOIPSL. !! !! ================================================================================================================================ SUBROUTINE routing_hr_graphinfo(filename, basmax, infmax, basmon, undef, tstream, tfast, tslow, tflood, cswamp, lfpcri) ! USE netcdf ! IMPLICIT NONE ! !! 0. Variables and parameter declaration !! 0.1 Input variables CHARACTER(LEN=*), INTENT(in) :: filename !! filename: name of the file to open INTEGER(i_std), INTENT(inout) :: basmax !! maximum number of HTUs INTEGER(i_std), INTENT(inout) :: basmon !! Number of HTUs to be monitored by grid box. INTEGER(i_std), INTENT(inout) :: infmax !! Maximum number of inflows. REAL(r_std), INTENT(out) :: undef REAL(r_std), INTENT(out) :: tstream, tfast, tslow, tflood, cswamp !! Time constants to be extracted REAL(r_std), INTENT(out) :: lfpcri !! Constant lim_floodcri to be taken from graph file. ! INTEGER(i_std) :: rcode, nid, dimid, ndims, nvars INTEGER(i_std), DIMENSION(4) :: dimids INTEGER(i_std) :: iv, ndimsvar CHARACTER(LEN=20) :: dname, varname ! ! IF (is_root_prc) THEN ! rcode = nf90_open(TRIM(filename), NF90_NOWRITE, nid) IF (rcode == NF90_NOERR) THEN ! ! Get graph file version ! rcode = nf90_get_att(nid, NF90_GLOBAL, "RoutingPPVersion", graphfile_version) IF (rcode /= NF90_NOERR) THEN graphfile_version = 0.0 ENDIF ! ! Assumes that the number of HTUs is in the dimension z ! rcode = nf90_inq_dimid(nid, "z", dimid) IF (rcode /= NF90_NOERR) CALL ipslerr_p(3, 'routing_hr_graphinfo', 'Error in nf90_inq_dimid for z', & TRIM(nf90_strerror(rcode)),'') rcode = nf90_inquire_dimension(nid, dimid, dname, basmax) IF (rcode /= NF90_NOERR) CALL ipslerr_p(3, 'routing_hr_graphinfo', 'Error in nf90_inquire_dimension basmax', & TRIM(nf90_strerror(rcode)),'') ! rcode = nf90_inq_dimid(nid, "inflow", dimid) IF (rcode /= NF90_NOERR) CALL ipslerr_p(3, 'routing_hr_graphinfo', 'Error in nf90_inq_dimid for inflow', & TRIM(nf90_strerror(rcode)),'') rcode = nf90_inquire_dimension(nid, dimid, dname, infmax) IF (rcode /= NF90_NOERR) CALL ipslerr_p(3, 'routing_hr_graphinfo', 'Error in nf90_inquire_dimension inflows', & TRIM(nf90_strerror(rcode)),'') ! rcode = nf90_inq_dimid(nid, "htumon", dimid) IF (rcode /= NF90_NOERR) THEN MonitoringinGraph = .FALSE. basmon = 1 ELSE rcode = nf90_inquire_dimension(nid, dimid, dname, basmon) IF (rcode /= NF90_NOERR) CALL ipslerr_p(3, 'routing_hr_graphinfo', 'Error in nf90_inquire_dimension for basmon', & TRIM(nf90_strerror(rcode)),'') MonitoringinGraph = .TRUE. ENDIF ! ! rcode = NF90_INQUIRE (nid, nDimensions=ndims, nVariables=nvars) IF (rcode /= NF90_NOERR) CALL ipslerr_p(3, 'routing_hr_graphinfo', 'Error in nf90_inquire', & TRIM(nf90_strerror(rcode)),'') ! DO iv=1,nvars ! rcode = NF90_INQUIRE_VARIABLE(nid, iv, name=varname, ndims=ndimsvar, dimids=dimids) ! SELECT CASE (varname) CASE ("basin_area") rcode = NF90_GET_ATT(nid, iv, "missing_value", undef) IF (rcode /= NF90_NOERR) THEN IF ( rcode == NF90_ENOTATT ) THEN rcode = NF90_GET_ATT(nid, iv, "_FillValue", undef) IF (rcode /= NF90_NOERR) CALL ipslerr_p(3, 'routing_hr_graphinfo', 'Did not get FillValue with nf90_get_att', & TRIM(nf90_strerror(rcode)),'') ELSE IF (rcode /= NF90_NOERR) CALL ipslerr_p(3, 'routing_hr_graphinfo', 'Error in nf90_get_att', & TRIM(nf90_strerror(rcode)),'') ENDIF ENDIF CASE("StreamTimeCst") rcode = NF90_GET_VAR(nid,iv,tstream) IF (rcode /= NF90_NOERR) CALL ipslerr_p(3, 'routing_hr_graphinfo', 'Error in nf90_get_var for variable StreamTimeCst', '','') ! If in an old version convert 10^3d/km into s/km IF (graphfile_version < 1.0) THEN tstream = tstream/1000*one_day ENDIF CASE("FastTimeCst") rcode = NF90_GET_VAR(nid,iv,tfast) IF (rcode /= NF90_NOERR) CALL ipslerr_p(3, 'routing_hr_graphinfo', 'Error in nf90_get_var for variable FastTimeCst', '','') ! If in an old version convert 10^3d/km into s/km IF (graphfile_version < 1.0) THEN tfast = tfast/1000*one_day ENDIF CASE("SlowTimeCst") rcode = NF90_GET_VAR(nid,iv,tslow) IF (rcode /= NF90_NOERR) CALL ipslerr_p(3, 'routing_hr_graphinfo', 'Error in nf90_get_var for variable SlowTimeCst', '','') ! If in an old version convert 10^3d/km into s/km IF (graphfile_version < 1.0) THEN tslow = tslow/1000*one_day ENDIF CASE("FloodTimeCst") rcode = NF90_GET_VAR(nid,iv,tflood) IF (rcode /= NF90_NOERR) CALL ipslerr_p(3, 'routing_hr_graphinfo', 'Error in nf90_get_var for variable FloodTimeCst', '','') ! If in an old version convert 10^3d/km into s/km IF (graphfile_version < 1.0) THEN tflood = tflood/1000*one_day ENDIF CASE("SwampCst") rcode = NF90_GET_VAR (nid,iv,cswamp) IF (rcode /= NF90_NOERR) CALL ipslerr_p(3, 'routing_hr_graphinfo', 'Error in nf90_get_var for variable SwampCst', '','') CASE("MaxTimeStep") rcode = NF90_GET_VAR (nid,iv,maxtimestep) IF (rcode /= NF90_NOERR) CALL ipslerr_p(3, 'routing_hr_graphinfo', 'Error in nf90_get_var for variable MaxTimeStep', '','') CASE("LimFloodcri") rcode = NF90_GET_VAR(nid,iv,lfpcri) IF (rcode /= NF90_NOERR) CALL ipslerr_p(3, 'routing_hr_graphinfo', 'Error in nf90_get_var for variable LimFloodcri', '','') END SELECT ENDDO rcode = NF90_CLOSE(nid) IF (rcode /= NF90_NOERR) CALL ipslerr_p(3, 'routing_hr_graphinfo', 'Error in nf90_close', & TRIM(nf90_strerror(rcode)),'') ! ! Before RoutingGraph version 2.5 the lfpcri parameter was hardcoded in routing.f90 and set to 2m. ! This information is preserved here. IF ( graphfile_version < 2.5 ) THEN lfpcri = 2.0 ENDIF ! ELSE ! Case without Graphfile. So we consider that the information will be found in the restart. CALL ipslerr_p(2, 'routing_hr_graphinfo', 'Could not open the rotung_graph.nc file', & "Expect to find all the information needed in the restart file.",'') ! MonitoringinGraph = .FALSE. basmax = -1 basmon = -1 infmax = -1 undef = undef_sechiba tstream = -1 tfast = -1 tslow = -1 tflood = -1 cswamp = -1 lfpcri = -1 ENDIF ENDIF !! CALL bcast(MonitoringinGraph) CALL bcast(basmax) CALL bcast(basmon) CALL bcast(infmax) CALL bcast(undef) CALL bcast(tstream) CALL bcast(tfast) CALL bcast(tslow) CALL bcast(tflood) CALL bcast(cswamp) CALL bcast(lfpcri) !! !! END SUBROUTINE routing_hr_graphinfo ! !! ================================================================================================================================ !! SUBROUTINE : routing_hr_indexfilegrid !! !>\BRIEF Locates all the points of the routing graph file on the model grid. This ensure that no assumption is made on the !! orientation of the grid in the file. !! !! ================================================================================================================================ SUBROUTINE routing_hr_indexfilegrid(im, jm, nbl, lon, lat, landindex, indextab, il2il) INTEGER(i_std), INTENT(IN) :: im, jm, nbl REAL(r_std), INTENT(IN) :: lon(im,jm), lat(im,jm) REAL(r_std), INTENT(IN) :: landindex(im,jm) INTEGER(i_std), INTENT(INOUT) :: indextab(im,jm) INTEGER(i_std), INTENT(OUT) :: il2il(nbl) ! INTEGER(i_std) :: il,i,j INTEGER(i_std) :: f(2) INTEGER(i_std) :: ih, jh, ir, jr REAL(r_std) :: nd REAL(r_std), DIMENSION(im,jm) :: dist REAL(r_std) :: mindist = 1000. !! Minimum distance in m between two points to be matched. ! indextab(:,:) = -1 ih = INT(im/2.) ir = NINT(im/2.)+1 jh = INT(jm/2.) jr = NINT(jm/2.)+1 ! DO il=1,nbl dist(:,:) = undef_sechiba DO i=MAX(1,ih-ir),MIN(im,ih+ir) DO j=MAX(1,jh-jr),MIN(jm,jh+jr) dist(i,j) = haversine_distance(lon(i,j), lat(i,j), lalo_g(il,2), lalo_g(il,1)) ENDDO ENDDO f=MINLOC(dist) IF ( dist(f(1),f(2)) < mindist ) THEN indextab(f(1),f(2)) = il il2il(NINT(landindex(f(1),f(2)))) = il dist(f(1),f(2)) = undef_sechiba ELSE CALL ipslerr(3,'routing_hr_indexfilegrid',& 'Distance of the closest point in the two grids is too large. ', & 'Are you sure the routing graph file is on the correct grid ?',& ' ') ENDIF ! ! See if the next point is close by ! nd = haversine_distance(lalo_g(il,2), lalo_g(il,1), lalo_g(MIN(il+1,nbl),2), lalo_g(MIN(il+1,nbl),1)) IF ( nd < MINVAL(dist)*3 ) THEN ! The next point is close so zoom in ih = f(1) ir = 4 jh = f(2) jr = 4 ELSE ! Back to starting conditions as the next point is far away ih = INT(im/2.) ir = NINT(im/2.)+1 jh = INT(jm/2.) jr = NINT(jm/2.)+1 ENDIF ENDDO END SUBROUTINE routing_hr_indexfilegrid ! !! ================================================================================================================================ !! SUBROUTINE : routing_hr_convertlandpts !! !>\BRIEF In case the order of land points was different in RoutingPreProc and the model. The route_togrid is corrected. !! !! ================================================================================================================================ ! SUBROUTINE routing_hr_convertlandpts(nbl, nbas, land2land, route_togrid) INTEGER(i_std), INTENT(IN) :: nbl, nbas INTEGER(i_std), INTENT(IN), DIMENSION(nbl) :: land2land INTEGER(i_std), INTENT(INOUT), DIMENSION(nbl,nbas) :: route_togrid ! INTEGER(i_std) :: ip, ib ! DO ip=1,nbl DO ib=1,nbas IF ( route_togrid(ip,ib) < undef_int .AND. route_togrid(ip,ib) > 0 ) THEN route_togrid(ip,ib) = land2land(route_togrid(ip,ib)) ELSE route_togrid(ip,ib) = ip ENDIF ENDDO ENDDO ! END SUBROUTINE routing_hr_convertlandpts ! !! ================================================================================================================================ !! SUBROUTINE : routing_hr_inflows !! !>\BRIEF Calculate the inflows from the outflows information. !! !! ================================================================================================================================ ! SUBROUTINE routing_hr_inflows(nbl, nbas, inf, floodplains_glo,route_innum_glo,route_ingrid_glo,route_inbasin_glo) IMPLICIT None INTEGER(i_std), INTENT(IN) :: nbl, nbas, inf REAL(r_std), INTENT(IN), DIMENSION(nbl,nbas) :: floodplains_glo INTEGER(i_std), INTENT(INOUT), DIMENSION(nbl,nbas) :: route_innum_glo INTEGER(i_std), INTENT(INOUT), DIMENSION(nbl,nbas, inf) :: route_ingrid_glo, route_inbasin_glo ! INTEGER(i_std) :: ig, ib, og, ob ! route_innum_glo(:,:) = 0 route_ingrid_glo(:,:,:) = 0 route_inbasin_glo(:,:,:) = 0 DO ig=1,nbl DO ib=1,nbas IF (floodplains_glo(ig,ib) .GT. 0) THEN og = route_togrid_glo(ig,ib) ob = route_tobasin_glo(ig,ib) IF (ob .LE. nbasmax) THEN IF (floodplains_glo(og,ob) .GT. 0) THEN route_innum_glo(og, ob) = route_innum_glo(og, ob) + 1 route_ingrid_glo(og,ob,route_innum_glo(og, ob)) = ig route_inbasin_glo(og,ob,route_innum_glo(og, ob)) = ib END IF END IF END IF ENDDO ENDDO END SUBROUTINE routing_hr_inflows ! !! !! ================================================================================================================================ !! SUBROUTINE : routing_hr_landgather !! !>\BRIEF Gathers the routing information onto landpoints, i.e. goes from an X/Y grid to a list of land points. !! !! ================================================================================================================================ ! SUBROUTINE routing_hr_landgather_r(im,jm,nbas,nbl,indextab,ijfield,landfield,def) ! INTEGER(i_std), INTENT(IN) :: im,jm,nbas,nbl INTEGER(i_std), INTENT(IN) :: indextab(im,jm) REAL(r_std), INTENT(IN), DIMENSION(im,jm,nbas) :: ijfield REAL(r_std), INTENT(OUT), DIMENSION(nbl,nbas) :: landfield REAL(r_std), INTENT(IN) :: def ! INTEGER(i_std) :: i,j,k ! DO i=1,im DO j=1,jm IF ( indextab(i,j) > 0 ) THEN DO k=1,nbas ! Catch undef or NaN values IF (ijfield(i,j,k) >= undef_graphfile .OR. ijfield(i,j,k) /= ijfield(i,j,k)) THEN landfield(indextab(i,j),k) = def ELSE landfield(indextab(i,j),k) = ijfield(i,j,k) ENDIF ENDDO ENDIF ENDDO ENDDO END SUBROUTINE routing_hr_landgather_r SUBROUTINE routing_hr_landgather_i2(im,jm,nbas,nbl,indextab,ijfield,landfield,def) ! INTEGER(i_std), INTENT(IN) :: im,jm,nbas,nbl INTEGER(i_std), INTENT(IN) :: indextab(im,jm) REAL(r_std), INTENT(IN), DIMENSION(im,jm,nbas) :: ijfield INTEGER(i_std), INTENT(OUT), DIMENSION(nbl,nbas) :: landfield INTEGER(i_std), INTENT(IN) :: def ! INTEGER(i_std) :: i,j,in ! DO i=1,im DO j=1,jm IF ( indextab(i,j) > 0 ) THEN DO in=1,nbas IF (ijfield(i,j,in) .GE. undef_int ) THEN landfield(indextab(i,j),in) = def ELSE landfield(indextab(i,j),in) = ijfield(i,j,in) ENDIF ENDDO ENDIF ENDDO ENDDO END SUBROUTINE routing_hr_landgather_i2 SUBROUTINE routing_hr_landgather_i1(im,jm,nbl,indextab,ijfield,landfield,def) ! INTEGER(i_std), INTENT(IN) :: im,jm,nbl INTEGER(i_std), INTENT(IN) :: indextab(im,jm) REAL(r_std), INTENT(IN), DIMENSION(im,jm) :: ijfield INTEGER(i_std), INTENT(OUT), DIMENSION(nbl) :: landfield INTEGER(i_std), INTENT(IN) :: def ! INTEGER(i_std) :: i,j ! DO i=1,im DO j=1,jm IF ( indextab(i,j) > 0 ) THEN IF (ijfield(i,j) .GE. undef_int ) THEN landfield(indextab(i,j)) = def ELSE landfield(indextab(i,j)) = ijfield(i,j) ENDIF ENDIF ENDDO ENDDO END SUBROUTINE routing_hr_landgather_i1 ! ! !! ================================================================================================================================ !! SUBROUTINE : routing_hr_irrigmap !! !>\BRIEF This subroutine interpolates the 0.5x0.5 degree based map of irrigated areas to the resolution of the model. !! !! DESCRIPTION (definitions, functional, design, flags) : None !! !! RECENT CHANGE(S): None !! !! MAIN OUTPUT VARIABLE(S): !! !! REFERENCES : None !! !! FLOWCHART : None !! \n !_ ================================================================================================================================ SUBROUTINE routing_hr_irrigmap (nbpt, index, lalo, neighbours, resolution, contfrac, & & init_irrig, irrigated, init_flood, init_swamp, swamp, hist_id, hist2_id) ! IMPLICIT NONE ! !! PARAMETERS INTEGER(i_std), PARAMETER :: ilake = 1 !! Number of type of lakes area (unitless) INTEGER(i_std), PARAMETER :: idam = 2 !! Number of type of dams area (unitless) INTEGER(i_std), PARAMETER :: iflood = 3 !! Number of type of floodplains area (unitless) INTEGER(i_std), PARAMETER :: iswamp = 4 !! Number of type of swamps area (unitless) INTEGER(i_std), PARAMETER :: isal = 5 !! Number of type of salines area (unitless) INTEGER(i_std), PARAMETER :: ipond = 6 !! Number of type of ponds area (unitless) INTEGER(i_std), PARAMETER :: ntype = 6 !! Number of types of flooded surfaces (unitless) !! INPUT VARIABLES INTEGER(i_std), INTENT(in) :: nbpt !! Domain size (unitless) INTEGER(i_std), INTENT(in) :: index(nbpt) !! Index on the global map. REAL(r_std), INTENT(in) :: lalo(nbpt,2) !! Vector of latitude and longitudes (beware of the order !) INTEGER(i_std), INTENT(in) :: neighbours(nbpt,NbNeighb)!! Vector of neighbours for each grid point 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) 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) LOGICAL, INTENT(in) :: init_irrig !! Logical to initialize the irrigation (true/false) LOGICAL, INTENT(in) :: init_flood !! Logical to initialize the floodplains (true/false) LOGICAL, INTENT(in) :: init_swamp !! Logical to initialize the swamps (true/false) ! !! OUTPUT VARIABLES REAL(r_std), INTENT(out) :: irrigated(:) !! Irrigated surface in each grid box (m^2) !! REAL(r_std), INTENT(out) :: floodplains(:) !! Surface which can be inundated in each grid box (m^2) REAL(r_std), INTENT(out) :: swamp(:) !! Surface which can be swamp in each grid box (m^2) ! !! LOCAL VARIABLES ! Interpolation variables ! INTEGER(i_std) :: nbpmax, nix, njx, fopt !! CHARACTER(LEN=30) :: callsign !! REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:) :: resol_lu !! Resolution read on the map INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:) :: mask !! Mask to exclude some points (unitless) REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: irrsub_area !! Area on the fine grid (m^2) INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:,:) :: irrsub_index !! Indices of the points we need on the fine grid (unitless) INTEGER :: ALLOC_ERR !! LOGICAL :: ok_interpol = .FALSE. !! Flag for interpolation (true/false) ! CHARACTER(LEN=80) :: filename !! Name of the netcdf file (unitless) INTEGER(i_std) :: iml, jml, lml, tml, fid, ib, ip, jp, itype !! Indices (unitless) REAL(r_std) :: lev(1), date, dt, coslat !! INTEGER(i_std) :: itau(1) !! REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: latrel !! Latitude REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: lonrel !! Longitude REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: irrigated_frac !! Irrigated fraction of the grid box (unitless;0-1) REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:) :: flood_fracmax !! Maximal flooded fraction of the grid box (unitless;0-1) REAL(r_std) :: area_irrig !! Irrigated surface in the grid box (m^2) REAL(r_std) :: area_flood(ntype) !! Flooded surface in the grid box (m^2) !!$ REAL(r_std) :: irrigmap(nbpt) !!$ REAL(r_std) :: swampmap(nbpt) !_ ================================================================================================================================ ! !Config Key = IRRIGATION_FILE !Config Desc = Name of file which contains the map of irrigated areas !Config Def = floodplains.nc !Config If = DO_IRRIGATION OR DO_FLOODPLAINS !Config Help = The name of the file to be opened to read the field !Config with the area in m^2 of the area irrigated within each !Config 0.5 0.5 deg grid box. The map currently used is the one !Config developed by the Center for Environmental Systems Research !Config in Kassel (1995). !Config Units = [FILE] ! filename = 'floodplains.nc' CALL getin_p('IRRIGATION_FILE',filename) ! IF (is_root_prc) THEN CALL flininfo(filename,iml, jml, lml, tml, fid) CALL flinclo(fid) ELSE iml = 0 jml = 0 lml = 0 tml = 0 ENDIF ! CALL bcast(iml) CALL bcast(jml) CALL bcast(lml) CALL bcast(tml) ! ! ! ALLOCATE (latrel(iml,jml), STAT=ALLOC_ERR) IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'routing_hr_irrigmap','Pb in allocate for latrel','','') ALLOCATE (lonrel(iml,jml), STAT=ALLOC_ERR) IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'routing_hr_irrigmap','Pb in allocate for lonrel','','') ALLOCATE (irrigated_frac(iml,jml), STAT=ALLOC_ERR) IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'routing_hr_irrigmap','Pb in allocate for irrigated_frac','','') ALLOCATE (flood_fracmax(iml,jml,ntype), STAT=ALLOC_ERR) IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'routing_hr_irrigmap','Pb in allocate for flood_fracmax','','') IF (is_root_prc) CALL flinopen(filename, .FALSE., iml, jml, lml, lonrel, latrel, lev, tml, itau, date, dt, fid) CALL bcast(lonrel) CALL bcast(latrel) ! IF (is_root_prc) CALL flinget(fid, 'irrig', iml, jml, lml, tml, 0, 0, irrigated_frac) CALL bcast(irrigated_frac) IF (is_root_prc) CALL flinget(fid, 'lake', iml, jml, lml, tml, 0, 0, flood_fracmax(:,:,ilake)) IF (is_root_prc) CALL flinget(fid, 'dam', iml, jml, lml, tml, 0, 0, flood_fracmax(:,:,idam)) IF (is_root_prc) CALL flinget(fid, 'flood', iml, jml, lml, tml, 0, 0, flood_fracmax(:,:,iflood)) IF (is_root_prc) CALL flinget(fid, 'swamp', iml, jml, lml, tml, 0, 0, flood_fracmax(:,:,iswamp)) IF (is_root_prc) CALL flinget(fid, 'saline', iml, jml, lml, tml, 0, 0, flood_fracmax(:,:,isal)) IF (is_root_prc) CALL flinget(fid, 'pond', iml, jml, lml, tml, 0, 0, flood_fracmax(:,:,ipond)) CALL bcast(flood_fracmax) ! IF (is_root_prc) CALL flinclo(fid) ! ! Set to zero all fraction which are less than 0.5% ! DO ip=1,iml DO jp=1,jml ! IF ( irrigated_frac(ip,jp) .LT. undef_sechiba-un) THEN irrigated_frac(ip,jp) = irrigated_frac(ip,jp)/100. IF ( irrigated_frac(ip,jp) < 0.005 ) irrigated_frac(ip,jp) = zero ENDIF ! DO itype=1,ntype IF ( flood_fracmax(ip,jp,itype) .LT. undef_sechiba-1.) THEN flood_fracmax(ip,jp,itype) = flood_fracmax(ip,jp,itype)/100 IF ( flood_fracmax(ip,jp,itype) < 0.005 ) flood_fracmax(ip,jp,itype) = zero ENDIF ENDDO ! ENDDO ENDDO IF (printlev>=2) THEN WRITE(numout,*) 'lonrel : ', MAXVAL(lonrel), MINVAL(lonrel) WRITE(numout,*) 'latrel : ', MAXVAL(latrel), MINVAL(latrel) WRITE(numout,*) 'irrigated_frac : ', MINVAL(irrigated_frac, MASK=irrigated_frac .GT. 0), & MAXVAL(irrigated_frac, MASK=irrigated_frac .LT. undef_sechiba) WRITE(numout,*) 'flood_fracmax : ', MINVAL(flood_fracmax, MASK=flood_fracmax .GT. 0), & MAXVAL(flood_fracmax, MASK=flood_fracmax .LT. undef_sechiba) END IF ! Consider all points a priori ! ALLOCATE(resol_lu(iml,jml,2), STAT=ALLOC_ERR) IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'routing_hr_irrigmap','Pb in allocate for resol_lu','','') ALLOCATE(mask(iml,jml), STAT=ALLOC_ERR) IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'routing_hr_irrigmap','Pb in allocate for mask','','') mask(:,:) = 0 DO ip=1,iml DO jp=1,jml ! ! Exclude the points where we are close to the missing value. ! !MG This condition cannot be applied in floodplains/swamps configuration because ! the same mask would be used for the interpolation of irrigation, floodplains and swamps maps. ! IF ( irrigated_frac(ip,jp) < undef_sechiba ) THEN mask(ip,jp) = 1 ! ENDIF ! ! Resolution in longitude ! coslat = MAX( COS( latrel(ip,jp) * pi/180. ), mincos ) IF ( ip .EQ. 1 ) THEN resol_lu(ip,jp,1) = ABS( lonrel(ip+1,jp) - lonrel(ip,jp) ) * pi/180. * R_Earth * coslat ELSEIF ( ip .EQ. iml ) THEN resol_lu(ip,jp,1) = ABS( lonrel(ip,jp) - lonrel(ip-1,jp) ) * pi/180. * R_Earth * coslat ELSE resol_lu(ip,jp,1) = ABS( lonrel(ip+1,jp) - lonrel(ip-1,jp) )/2. * pi/180. * R_Earth * coslat ENDIF ! ! Resolution in latitude ! IF ( jp .EQ. 1 ) THEN resol_lu(ip,jp,2) = ABS( latrel(ip,jp) - latrel(ip,jp+1) ) * pi/180. * R_Earth ELSEIF ( jp .EQ. jml ) THEN resol_lu(ip,jp,2) = ABS( latrel(ip,jp-1) - latrel(ip,jp) ) * pi/180. * R_Earth ELSE resol_lu(ip,jp,2) = ABS( latrel(ip,jp-1) - latrel(ip,jp+1) )/2. * pi/180. * R_Earth ENDIF ! ENDDO ENDDO ! ! The number of maximum vegetation map points in the GCM grid is estimated. ! Some lmargin is taken. ! callsign = 'Irrigation map' ok_interpol = .FALSE. IF (is_root_prc) THEN nix=INT(MAXVAL(resolution_g(:,1))/MAXVAL(resol_lu(:,:,1)))+2 njx=INT(MAXVAL(resolution_g(:,2))/MAXVAL(resol_lu(:,:,2)))+2 nbpmax = nix*njx*2 IF (printlev>=1) THEN WRITE(numout,*) "Projection arrays for ",callsign," : " WRITE(numout,*) "nbpmax = ",nbpmax, nix, njx END IF ENDIF CALL bcast(nbpmax) ALLOCATE(irrsub_index(nbpt, nbpmax, 2), STAT=ALLOC_ERR) IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'routing_hr_irrigmap','Pb in allocate for irrsub_index','','') irrsub_index(:,:,:)=0 ALLOCATE(irrsub_area(nbpt, nbpmax), STAT=ALLOC_ERR) IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'routing_hr_irrigmap','Pb in allocate for irrsub_area','','') irrsub_area(:,:)=zero CALL aggregate_p(nbpt, lalo, neighbours, resolution, contfrac, & & iml, jml, lonrel, latrel, mask, callsign, & & nbpmax, irrsub_index, irrsub_area, ok_interpol) ! ! WHERE (irrsub_area < 0) irrsub_area=zero ! ! Test here if not all sub_area are larger than 0 if so, then we need to increase nbpmax ! DO ib=1,nbpt ! area_irrig = 0.0 area_flood = 0.0 ! DO fopt=1,COUNT(irrsub_area(ib,:) > zero) ! ip = irrsub_index(ib, fopt, 1) jp = irrsub_index(ib, fopt, 2) ! IF (irrigated_frac(ip,jp) .LT. undef_sechiba-1.) THEN area_irrig = area_irrig + irrsub_area(ib,fopt)*irrigated_frac(ip,jp) ENDIF ! DO itype=1,ntype IF (flood_fracmax(ip,jp,itype) .LT. undef_sechiba-1.) THEN area_flood(itype) = area_flood(itype) + irrsub_area(ib,fopt)*flood_fracmax(ip,jp,itype) ENDIF ENDDO ENDDO ! ! Put the total irrigated and flooded areas in the output variables ! IF ( init_irrig ) THEN irrigated(ib) = MIN(area_irrig, resolution(ib,1)*resolution(ib,2)*contfrac(ib)) IF ( irrigated(ib) < 0 ) THEN WRITE(numout,*) 'We have a problem here : ', irrigated(ib) WRITE(numout,*) 'resolution :', resolution(ib,1), resolution(ib,2) WRITE(numout,*) area_irrig CALL ipslerr_p(3,'routing_hr_irrigmap','Problem with irrigated...','','') ENDIF !!$ ! Compute a diagnostic of the map. !!$ IF(contfrac(ib).GT.zero) THEN !!$ irrigmap (ib) = irrigated(ib) / ( resolution(ib,1)*resolution(ib,2)*contfrac(ib) ) !!$ ELSE !!$ irrigmap (ib) = zero !!$ ENDIF ! ENDIF ! ! ! IF ( init_swamp ) THEN swamp(ib) = MIN(area_flood(iswamp), resolution(ib,1)*resolution(ib,2)*contfrac(ib)) IF ( swamp(ib) < 0 ) THEN WRITE(numout,*) 'We have a problem here : ', swamp(ib) WRITE(numout,*) 'resolution :', resolution(ib,1), resolution(ib,2) WRITE(numout,*) area_flood CALL ipslerr_p(3,'routing_hr_irrigmap','Problem with swamp...','','') ENDIF !!$ ! Compute a diagnostic of the map. !!$ IF(contfrac(ib).GT.zero) THEN !!$ swampmap(ib) = swamp(ib) / ( resolution(ib,1)*resolution(ib,2)*contfrac(ib) ) !!$ ELSE !!$ swampmap(ib) = zero !!$ ENDIF ENDIF ! ! ENDDO ! ! IF (printlev>=1) THEN IF ( init_irrig ) WRITE(numout,*) "Diagnostics irrigated :", MINVAL(irrigated), MAXVAL(irrigated) IF ( init_flood ) WRITE(numout,*) "Diagnostics floodplains :", MINVAL(floodplains), MAXVAL(floodplains) IF ( init_swamp ) WRITE(numout,*) "Diagnostics swamp :", MINVAL(swamp), MAXVAL(swamp) END IF ! No compensation is done for overlapping floodplains, swamp and irrig. At least overlapping will not ! happen between floodplains and swamp alone ! IF ( init_irrig .AND. init_flood ) THEN ! DO ib = 1, nbpt ! surp = (floodplains(ib)+swamp(ib)+irrigated(ib)) / (resolution(ib,1)*resolution(ib,2)*contfrac(ib)) ! IF ( surp .GT. un ) THEN ! floodplains(ib) = floodplains(ib) / surp ! swamp(ib) = swamp(ib) / surp ! irrigated(ib) = irrigated(ib) / surp ! ENDIF ! ENDDO ! ENDIF ! DEALLOCATE (irrsub_area) DEALLOCATE (irrsub_index) ! DEALLOCATE (mask) DEALLOCATE (resol_lu) ! DEALLOCATE (lonrel) DEALLOCATE (latrel) ! END SUBROUTINE routing_hr_irrigmap ! !! ================================================================================================================================ !! SUBROUTINE : routing_hr_waterbal !! !>\BRIEF This subroutine checks the water balance in the routing module. !! !! DESCRIPTION (definitions, functional, design, flags) : None !! !! RECENT CHANGE(S): None !! !! MAIN OUTPUT VARIABLE(S): !! !! REFERENCES : None !! !! FLOWCHART : None !! \n !_ ================================================================================================================================ SUBROUTINE routing_hr_waterbal(nbpt, reinit, floodout, runoff, drainage, returnflow, & & reinfiltration, irrigation, riverflow, coastalflow) ! IMPLICIT NONE ! !! INPUT VARIABLES INTEGER(i_std), INTENT(in) :: nbpt !! Domain size (unitless) LOGICAL, INTENT(in) :: reinit !! Controls behaviour (true/false) 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) :: returnflow(nbpt) !! The water flow from lakes and swamps which returns to the grid box. !! This water will go back into the hydrol module to allow re-evaporation (kg/m^2/dt) REAL(r_std), INTENT(in) :: reinfiltration(nbpt) !! Water flow from ponds and floodplains which returns to the grid box (kg/m^2/dt) REAL(r_std), INTENT(in) :: 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(in) :: 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(in) :: coastalflow(nbpt) !! Outflow on coastal points by small basins. This is the water which flows in a disperse way into the ocean (kg/dt) ! ! We sum-up all the water we have in the warious reservoirs ! REAL(r_std), SAVE :: totw_flood !! Sum of all the water amount in the floodplains reservoirs (kg) !$OMP THREADPRIVATE(totw_flood) REAL(r_std), SAVE :: totw_stream !! Sum of all the water amount in the stream reservoirs (kg) !$OMP THREADPRIVATE(totw_stream) REAL(r_std), SAVE :: totw_fast !! Sum of all the water amount in the fast reservoirs (kg) !$OMP THREADPRIVATE(totw_fast) REAL(r_std), SAVE :: totw_slow !! Sum of all the water amount in the slow reservoirs (kg) !$OMP THREADPRIVATE(totw_slow) REAL(r_std), SAVE :: totw_lake !! Sum of all the water amount in the lake reservoirs (kg) !$OMP THREADPRIVATE(totw_lake) REAL(r_std), SAVE :: totw_pond !! Sum of all the water amount in the pond reservoirs (kg) !$OMP THREADPRIVATE(totw_pond) REAL(r_std), SAVE :: totw_in !! Sum of the water flow in to the routing scheme !$OMP THREADPRIVATE(totw_in) REAL(r_std), SAVE :: totw_out !! Sum of the water flow out to the routing scheme !$OMP THREADPRIVATE(totw_out) REAL(r_std), SAVE :: totw_return !! !$OMP THREADPRIVATE(totw_return) REAL(r_std), SAVE :: totw_irrig !! !$OMP THREADPRIVATE(totw_irrig) REAL(r_std), SAVE :: totw_river !! !$OMP THREADPRIVATE(totw_river) REAL(r_std), SAVE :: totw_coastal !! !$OMP THREADPRIVATE(totw_coastal) REAL(r_std) :: totarea !! Total area of basin (m^2) REAL(r_std) :: area !! Total area of routing (m^2) INTEGER(i_std) :: ig !! ! ! Just to make sure we do not get too large numbers ! ! !! PARAMETERS REAL(r_std), PARAMETER :: scaling = 1.0E+6 !! REAL(r_std), PARAMETER :: allowed_err = 50. !! !_ ================================================================================================================================ ! IF ( reinit ) THEN ! totw_flood = zero totw_stream = zero totw_fast = zero totw_slow = zero totw_lake = zero totw_pond = zero totw_in = zero ! DO ig=1,nbpt ! totarea = SUM(routing_area(ig,:)) ! totw_flood = totw_flood + SUM(flood_reservoir(ig,:)/scaling) totw_stream = totw_stream + SUM(stream_reservoir(ig,:)/scaling) totw_fast = totw_fast + SUM(fast_reservoir(ig,:)/scaling) totw_slow = totw_slow + SUM(slow_reservoir(ig,:)/scaling) totw_lake = totw_lake + lake_reservoir(ig)/scaling totw_pond = totw_pond + pond_reservoir(ig)/scaling ! totw_in = totw_in + (runoff(ig)*totarea + drainage(ig)*totarea - floodout(ig)*totarea)/scaling ! ENDDO ! ELSE ! totw_out = zero totw_return = zero totw_irrig = zero totw_river = zero totw_coastal = zero area = zero ! DO ig=1,nbpt ! totarea = SUM(routing_area(ig,:)) ! totw_flood = totw_flood - SUM(flood_reservoir(ig,:)/scaling) totw_stream = totw_stream - SUM(stream_reservoir(ig,:)/scaling) totw_fast = totw_fast - SUM(fast_reservoir(ig,:)/scaling) totw_slow = totw_slow - SUM(slow_reservoir(ig,:)/scaling) totw_lake = totw_lake - lake_reservoir(ig)/scaling totw_pond = totw_pond - pond_reservoir(ig)/scaling ! totw_return = totw_return + (reinfiltration(ig)+returnflow(ig))*totarea/scaling totw_irrig = totw_irrig + irrigation(ig)*totarea/scaling totw_river = totw_river + riverflow(ig)/scaling totw_coastal = totw_coastal + coastalflow(ig)/scaling ! area = area + totarea ! ENDDO totw_out = totw_return + totw_irrig + totw_river + totw_coastal ! ! Now we have all the information to balance our water ! IF ( ABS((totw_flood + totw_stream + totw_fast + totw_slow + totw_lake + totw_pond) - & & (totw_out - totw_in)) > allowed_err ) THEN WRITE(numout,*) 'WARNING : Water not conserved in routing. Limit at ', allowed_err, ' 10^6 kg' WRITE(numout,*) '--Water-- change : flood stream fast ', totw_flood, totw_stream, totw_fast WRITE(numout,*) '--Water-- change : slow, lake ', totw_slow, totw_lake WRITE(numout,*) '--Water>>> change in the routing res. : ', totw_flood + totw_stream + totw_fast + totw_slow + totw_lake WRITE(numout,*) '--Water input : ', totw_in WRITE(numout,*) '--Water output : ', totw_out WRITE(numout,*) '--Water output : return, irrig ', totw_return, totw_irrig WRITE(numout,*) '--Water output : river, coastal ',totw_river, totw_coastal WRITE(numout,*) '--Water>>> change by fluxes : ', totw_out - totw_in, ' Diff [mm/dt]: ', & & ((totw_flood + totw_stream + totw_fast + totw_slow + totw_lake) - (totw_out - totw_in))/area ! Stop the model CALL ipslerr_p(3, 'routing_hr_waterbal', 'Water is not conserved in routing.','','') ENDIF ! ENDIF ! END SUBROUTINE routing_hr_waterbal ! ! END MODULE routing_highres