!! !! !! This module routes the water over the continents into the oceans. !! !! 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 !! !! @author Jan Polcher !! @Version : $Revision: 1.41 $, $Date: 2009/01/07 13:39:45 $ !! !! $Header: /home/ssipsl/CVSREP/ORCHIDEE/src_sechiba/routing.f90,v 1.41 2009/01/07 13:39:45 ssipsl Exp $ !! IPSL (2006) !! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC !! MODULE routing ! ! ! routines called : restput, restget ! USE ioipsl ! USE constantes USE constantes_veg USE sechiba_io USE grid USE parallel ! IMPLICIT NONE ! ! public routines : ! PRIVATE PUBLIC :: routing_main, routing_clear ! ! variables used inside hydrol module : declaration and initialisation ! LOGICAL, SAVE :: l_first_routing=.TRUE. !! Initialisation has to be done one time LOGICAL, SAVE :: check_waterbal=.FALSE. !! The check the water balance ! ! The maximum number of basins we wish to have per grid-box (truncation of the model) INTEGER(i_std), PARAMETER :: nbasmax=5 ! The maximum number of bassins we can handle at any time during the generation of the maps. INTEGER(i_std), PARAMETER :: nbvmax = 220 ! ! The time constants are in days. ! REAL(r_std), PARAMETER :: fast_tcst = 3.0, slow_tcst = 25.0, stream_tcst = 0.24 REAL(r_std), PARAMETER :: evap_cst = 0.18, wdelay_cst = 0.7 ! ! Maximum evaporation rate from lakes 7.5 kg/m^2/d transformed into kg/m^2/sec ! REAL(r_std), PARAMETER :: maxevap_lake = 7.5/86400. ! ! Parameter for the Kassel irrigation parametrization linked to the crops ! REAL(r_std), PARAMETER :: crop_coef = 1.5 ! INTEGER(i_std), PARAMETER :: undef_int = 999999999 ! REAL(r_std),SAVE :: dt_routing ! ! Logicals to control model configuration ! LOGICAL, SAVE :: doirrigation = .FALSE. LOGICAL, SAVE :: dofloodplains = .FALSE. ! ! ! The variables describing the basins and their routing, need to be in the restart file. ! INTEGER(i_std), SAVE :: num_largest REAL(r_std), SAVE :: time_counter REAL(r_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: routing_area_loc REAL(r_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: topo_resid_loc INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: route_togrid_loc INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: route_tobasin_loc INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: global_basinid_loc INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: hydrodiag_loc !! Variable to diagnose the hydrographs ! ! parallelism REAL(r_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: routing_area_glo REAL(r_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: topo_resid_glo INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: route_togrid_glo INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: route_tobasin_glo INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: global_basinid_glo INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:) :: hydrodiag_glo !! Variable to diagnose the hydrographs REAL(r_std), SAVE, POINTER, DIMENSION(:,:) :: routing_area REAL(r_std), SAVE, POINTER, DIMENSION(:,:) :: topo_resid INTEGER(i_std), SAVE, POINTER, DIMENSION(:,:) :: route_togrid INTEGER(i_std), SAVE, POINTER, DIMENSION(:,:) :: route_tobasin INTEGER(i_std), SAVE, POINTER, DIMENSION(:,:) :: global_basinid INTEGER(i_std), SAVE, POINTER, DIMENSION(:,:) :: hydrodiag ! Map of irrigated areas and floodplains ! REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: irrigated ! Surface which can be irrigate in each grid-box REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: floodplains ! Surface which can be inondated in each grid-box REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: previous_outflow ! ! The reservoirs, also to be put into the restart file. ! REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: fast_reservoir, slow_reservoir, stream_reservoir REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: lake_reservoir ! ! The accumulated fluxes. ! REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: runoff_mean, drainage_mean, evapot_mean REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: precip_mean, humrel_mean, totnobio_mean, vegtot_mean ! ! The averaged outflow fluxes. ! REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: lakeinflow_mean REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: returnflow_mean, irrigation_mean REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: riverflow_mean, coastalflow_mean REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: floodtemp INTEGER(i_std), SAVE :: floodtemp_lev ! ! Diagnostic variables ... well sort of ! ! REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: irrig_netereq ! ! Hydrographs at the outflow of the grid box for major basins REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: hydrographs ! Diagnostics for the various reservoirs we use (Kg/m^2) REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: fast_diag, slow_diag, stream_diag, lake_diag ! CONTAINS ! !--------------------------------------------------------------------------------- ! SUBROUTINE routing_main(kjit, nbpt, dtradia, ldrestart_read, ldrestart_write, index, & & lalo, neighbours, resolution, contfrac, totfrac_nobio, veget_max, runoff, & & drainage, evapot, precip_rain, humrel, & & stempdiag, returnflow, irrigation, riverflow, coastalflow, rest_id, hist_id, hist2_id) ! IMPLICIT NONE ! ! This module will route the runoff and drainage produced by the hydrol module. These two ! fluxes are provided in kg/m^2dt. 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 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 river reservoir and beeing put into the upper ! layers of the soil. ! The two first fluxes are in kg/dt and the last one is on kg/m^2dt. ! INTEGER(i_std), INTENT(in) :: kjit !! Time step number INTEGER(i_std), INTENT(in) :: nbpt !! Domain size INTEGER(i_std),INTENT(in) :: rest_id,hist_id !! _Restart_ file and _history_ file identifier INTEGER(i_std),INTENT(in) :: hist2_id !! _history_ file 2 identifier REAL(r_std), INTENT(in) :: dtradia !! Time step in seconds LOGICAL, INTENT(in) :: ldrestart_read !! Logical for _restart_ file to read LOGICAL, INTENT(in) :: ldrestart_write !! Logical for _restart_ file to write INTEGER(i_std), INTENT(in) :: index(nbpt) ! Indeces of the points on the 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,8) ! Vector of neighbours for each grid point (1=N, 2=E, 3=S, 4=W) REAL(r_std), INTENT(in) :: resolution(nbpt,2) ! The size in km of each grid-box in X and Y REAL(r_std), INTENT(in) :: contfrac(nbpt) ! Fraction of land in each grid box REAL(r_std), INTENT(in) :: totfrac_nobio(nbpt) ! Total fraction of continental ice+lakes+... REAL(r_std), INTENT(in) :: veget_max(nbpt,nvm) ! Maximum vegetation fraction. We want to have the ! part of the grid which can have some vegetation. REAL(r_std), INTENT(in) :: runoff(nbpt) ! grid-point runoff REAL(r_std), INTENT(in) :: drainage(nbpt) ! grid-point drainage REAL(r_std), INTENT(in) :: evapot(nbpt) ! Potential evaporation REAL(r_std), INTENT(in) :: precip_rain(nbpt) ! Rainfall needed for the irrigation formula REAL(r_std), INTENT(in) :: humrel(nbpt,nvm) ! Soil moisture stress REAL(r_std), INTENT(in) :: stempdiag(nbpt,nbdl)! Temperature profile in soil ! REAL(r_std), INTENT(out) :: returnflow(nbpt) ! The water flow which returns to the grid box (kg/m^2 per dt) REAL(r_std), INTENT(out) :: irrigation(nbpt) ! Irrigation flux (kg/m^2 per dt) REAL(r_std), INTENT(out) :: riverflow(nbpt) ! Outflow of the major rivers REAL(r_std), INTENT(out) :: coastalflow(nbpt) ! Outflow on coastal points by small basins ! ! LOCAL ! CHARACTER(LEN=30) :: var_name REAL(r_std), DIMENSION(1) :: tmp_day REAL(r_std), DIMENSION(nbpt) :: return_lakes INTEGER(i_std) :: ig, jv LOGICAL, SAVE :: init_irrig=.FALSE., init_flood=.FALSE. ! ! do initialisation ! IF (l_first_routing) THEN ! ! 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_init (kjit, nbpt, index, dtradia, returnflow, irrigation, & & riverflow, coastalflow, stempdiag, rest_id) routing_area => routing_area_loc topo_resid => topo_resid_loc route_togrid => route_togrid_loc route_tobasin => route_tobasin_loc global_basinid => global_basinid_loc hydrodiag => hydrodiag_loc ! ! This routine computes the routing map if needed. ! IF ( COUNT(route_togrid_glo .GE. undef_int) .GT. 0 ) THEN CALL routing_basins_p(nbpt, lalo, neighbours, resolution, contfrac) ENDIF ! ! Do we have what we need if we want to do irrigation ! IF ( doirrigation .OR. dofloodplains ) THEN IF ( doirrigation ) THEN IF (COUNT(irrigated .GE. undef_sechiba-1) > 0) THEN init_irrig = .TRUE. ENDIF ENDIF IF ( dofloodplains ) THEN IF (COUNT(floodplains .GE. undef_sechiba-1) > 0) THEN init_flood = .TRUE. ENDIF ENDIF IF ( init_irrig .OR. init_flood ) THEN CALL routing_irrigmap(nbpt, index, lalo, neighbours, resolution, & contfrac, init_irrig, irrigated, init_flood, floodplains, hist_id, hist2_id) ENDIF ENDIF ! ! This routine gives a diag nostic of the basins used. ! CALL routing_diagnostic_p(nbpt, index, resolution, contfrac, hist_id, hist2_id) ! l_first_routing = .FALSE. ! RETURN ! ENDIF ! ! Accumulate the incoming water fluxes ! runoff_mean(:) = runoff_mean(:) + runoff(:) drainage_mean(:) = drainage_mean(:) + drainage(:) evapot_mean(:) = evapot_mean(:) + evapot(:) floodtemp(:) = stempdiag(:,floodtemp_lev) precip_mean(:) = precip_mean(:) + precip_rain(:) ! ! Averaged variables (i.e. *dtradia/dt_routing) ! totnobio_mean(:) = totnobio_mean(:) + totfrac_nobio(:)*dtradia/dt_routing ! ! Only potentially vegetated surfaces are taken into account. At the start of ! the growing seasons (when veget and veget_max are still different) we will ! 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)*dtradia/dt_routing vegtot_mean(ig) = vegtot_mean(ig) + veget_max(ig,jv)*dtradia/dt_routing ENDDO ENDDO ! time_counter = time_counter + dtradia ! ! If the time has come we do the routing. ! IF ( NINT(time_counter) .GE. NINT(dt_routing) ) THEN ! ! Check the water balance if needed ! IF ( check_waterbal ) THEN CALL routing_waterbal(nbpt, .TRUE., runoff_mean, drainage_mean, returnflow_mean, & & irrigation_mean, riverflow_mean, coastalflow_mean) ENDIF ! ! Make sure we do not flood north of 49N ! DO ig=1,nbpt IF ( lalo(ig,1) > 49.0 ) THEN floodtemp(ig) = tp_00 -1. ENDIF ENDDO ! CALL routing_flow(nbpt, dt_routing, runoff_mean, drainage_mean, & & vegtot_mean, totnobio_mean, evapot_mean, precip_mean, humrel_mean, floodtemp, & & lakeinflow_mean, returnflow_mean, irrigation_mean, riverflow_mean, & & coastalflow_mean, hydrographs) ! CALL routing_lake(nbpt, dt_routing, lakeinflow_mean, humrel_mean, return_lakes) ! returnflow_mean(:) = returnflow_mean(:) + return_lakes(:) ! ! Close the water balance if asked for ! IF ( check_waterbal ) THEN CALL routing_waterbal(nbpt, .FALSE., runoff_mean, drainage_mean, returnflow_mean, & & irrigation_mean, riverflow_mean, coastalflow_mean) ENDIF ! time_counter = zero ! runoff_mean(:) = zero drainage_mean(:) = zero evapot_mean(:) = zero precip_mean(:) = zero ! humrel_mean(:) = zero totnobio_mean(:) = zero vegtot_mean(:) = zero ! ! Change the units of the routing fluxes from kg/dt_routing into kg/dtradia ! and from m^3/dt_routing into m^3/dtradia ! returnflow_mean(:) = returnflow_mean(:)/dt_routing*dtradia irrigation_mean(:) = irrigation_mean(:)/dt_routing*dtradia irrig_netereq(:) = irrig_netereq(:)/dt_routing*dtradia ! ! riverflow_mean(:) = riverflow_mean(:)/dt_routing*dtradia coastalflow_mean(:) = coastalflow_mean(:)/dt_routing*dtradia hydrographs(:) = hydrographs(:)/dt_routing*dtradia ! ! Convert from kg/dtradia to m^3/dtradia ! hydrographs(:) = hydrographs(:)/1000. ! ENDIF ! ! Return the fraction of routed water for this time step. ! returnflow(:) = returnflow_mean(:) irrigation(:) = irrigation_mean(:) riverflow(:) = riverflow_mean(:) coastalflow(:) = coastalflow_mean(:) ! ! Write restart ! IF (ldrestart_write) THEN ! var_name ="routingcounter" tmp_day(1) = time_counter IF (is_root_prc) CALL restput (rest_id, var_name, 1, 1, 1, kjit, tmp_day) ! var_name = 'routingarea' CALL restput_p (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, routing_area, 'scatter', nbp_glo, index_g) var_name = 'routetogrid' CALL restput_p (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, REAL(route_togrid,r_std), 'scatter', & & nbp_glo, index_g) var_name = 'routetobasin' CALL restput_p (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, REAL(route_tobasin,r_std), 'scatter', & & nbp_glo, index_g) var_name = 'basinid' CALL restput_p (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, REAL(global_basinid,r_std), 'scatter', & & nbp_glo, index_g) var_name = 'topoindex' CALL restput_p (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, topo_resid, 'scatter', nbp_glo, index_g) var_name = 'fastres' CALL restput_p (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, fast_reservoir, 'scatter', nbp_glo, index_g) var_name = 'slowres' CALL restput_p (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, slow_reservoir, 'scatter', nbp_glo, index_g) var_name = 'streamres' CALL restput_p (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, stream_reservoir, 'scatter',nbp_glo,index_g) ! var_name = 'lakeres' CALL restput_p (rest_id, var_name, nbp_glo, 1, 1, kjit, lake_reservoir, 'scatter', nbp_glo, index_g) ! var_name = 'lakeinflow' CALL restput_p (rest_id, var_name, nbp_glo, 1, 1, kjit, lakeinflow_mean, 'scatter', nbp_glo, index_g) var_name = 'returnflow' CALL restput_p (rest_id, var_name, nbp_glo, 1, 1, kjit, returnflow_mean, 'scatter', nbp_glo, index_g) var_name = 'riverflow' CALL restput_p (rest_id, var_name, nbp_glo, 1, 1, kjit, riverflow_mean, 'scatter', nbp_glo, index_g) var_name = 'coastalflow' CALL restput_p (rest_id, var_name, nbp_glo, 1, 1, kjit, coastalflow_mean, 'scatter', nbp_glo, index_g) var_name = 'hydrographs' CALL restput_p (rest_id, var_name, nbp_glo, 1, 1, kjit, hydrographs, 'scatter', nbp_glo, index_g) ! ! Keep track of the accumulated variables ! var_name = 'runoff_route' CALL restput_p (rest_id, var_name, nbp_glo, 1, 1, kjit, runoff_mean, 'scatter', nbp_glo, index_g) var_name = 'drainage_route' CALL restput_p (rest_id, var_name, nbp_glo, 1, 1, kjit, drainage_mean, 'scatter', nbp_glo, index_g) var_name = 'evapot_route' CALL restput_p (rest_id, var_name, nbp_glo, 1, 1, kjit, evapot_mean, 'scatter', nbp_glo, index_g) var_name = 'precip_route' CALL restput_p (rest_id, var_name, nbp_glo, 1, 1, kjit, precip_mean, 'scatter', nbp_glo, index_g) var_name = 'humrel_route' CALL restput_p (rest_id, var_name, nbp_glo, 1, 1, kjit, humrel_mean, 'scatter', nbp_glo, index_g) var_name = 'totnobio_route' CALL restput_p (rest_id, var_name, nbp_glo, 1, 1, kjit, totnobio_mean, 'scatter', nbp_glo, index_g) var_name = 'vegtot_route' CALL restput_p (rest_id, var_name, nbp_glo, 1, 1, kjit, vegtot_mean, 'scatter', nbp_glo, index_g) ! var_name = 'previous_outflow' CALL restput_p (rest_id, var_name, nbp_glo, nbasmax+3, 1, kjit,previous_outflow,'scatter',nbp_glo,index_g) ! IF ( doirrigation ) THEN var_name = 'irrigated' CALL restput_p (rest_id, var_name, nbp_glo, 1, 1, kjit, irrigated, 'scatter', nbp_glo, index_g) var_name = 'irrigation' CALL restput_p (rest_id, var_name, nbp_glo, 1, 1, kjit, irrigation_mean, 'scatter', nbp_glo, index_g) ENDIF ! IF ( dofloodplains ) THEN var_name = 'floodplains' CALL restput_p (rest_id, var_name, nbp_glo, 1, 1, kjit, floodplains, 'scatter', nbp_glo, index_g) ENDIF ! RETURN ! ENDIF ! ! Write diagnostics ! IF ( .NOT. almaoutput ) THEN ! CALL histwrite(hist_id, 'riversret', kjit, returnflow, nbpt, index) CALL histwrite(hist_id, 'hydrographs', kjit, hydrographs, nbpt, index) ! CALL histwrite(hist_id, 'fastr', kjit, fast_diag, nbpt, index) CALL histwrite(hist_id, 'slowr', kjit, slow_diag, nbpt, index) CALL histwrite(hist_id, 'streamr', kjit, stream_diag, nbpt, index) CALL histwrite(hist_id, 'lakevol', kjit, lake_diag, nbpt, index) ! CALL histwrite(hist_id, 'irrigation', kjit, irrigation, nbpt, index) ! CALL histwrite(hist_id, 'netirrig', kjit, irrig_netereq, nbpt, index) ! ELSE ! CALL histwrite(hist_id, 'dis', kjit, hydrographs, nbpt, index) ! ENDIF IF ( hist2_id > 0 ) THEN IF ( .NOT. almaoutput ) THEN ! CALL histwrite(hist2_id, 'riversret', kjit, returnflow, nbpt, index) CALL histwrite(hist2_id, 'hydrographs', kjit, hydrographs, nbpt, index) ! CALL histwrite(hist2_id, 'fastr', kjit, fast_diag, nbpt, index) CALL histwrite(hist2_id, 'slowr', kjit, slow_diag, nbpt, index) CALL histwrite(hist2_id, 'streamr', kjit, stream_diag, nbpt, index) CALL histwrite(hist2_id, 'lakevol', kjit, lake_diag, nbpt, index) ! CALL histwrite(hist2_id, 'irrigation', kjit, irrigation, nbpt, index) ! CALL histwrite(hist2_id, 'netirrig', kjit, irrig_netereq, nbpt, index) ! ELSE ! CALL histwrite(hist2_id, 'dis', kjit, hydrographs, nbpt, index) ! ENDIF ENDIF ! ! END SUBROUTINE routing_main ! !--------------------------------------------------------------------------------- ! SUBROUTINE routing_init(kjit, nbpt, index, dtradia, returnflow, irrigation, & & riverflow, coastalflow, stempdiag, rest_id) ! IMPLICIT NONE ! ! interface description ! input ! INTEGER(i_std), INTENT(in) :: kjit !! Time step number INTEGER(i_std), INTENT(in) :: nbpt !! Domain size INTEGER(i_std), DIMENSION (nbpt), INTENT(in) :: index !! Indeces of the points on the map REAL(r_std), INTENT(in) :: dtradia !! timestep REAL(r_std), DIMENSION (nbpt),INTENT(out) :: returnflow !! The water flow which returns to the grid box (kg/m^2 per dt) REAL(r_std), DIMENSION (nbpt),INTENT(out) :: irrigation !! Irrigation flow (kg/m^2 per dt) REAL(r_std), DIMENSION (nbpt),INTENT(out) :: riverflow !! Outflow of the major rivers REAL(r_std), DIMENSION (nbpt),INTENT(out) :: coastalflow !! Outflow on coastal points by small basins REAL(r_std), DIMENSION(nbpt,nbdl),INTENT(in) :: stempdiag !! Temperature profile in soil INTEGER(i_std), INTENT(in) :: rest_id !! Restart file identifier ! ! local declaration ! CHARACTER(LEN=80) :: var_name !! To store variables names for I/O REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: tmp_real_g !! A temporary real array for the integers REAL(r_std), DIMENSION(1) :: tmp_day REAL(r_std) :: ratio, totarea INTEGER(i_std) :: ier, ig, ib, ipn(1) ! ! These variables will require the configuration infrastructure ! !Config Key = ROUTING_TIMESTEP !Config If = RIVER_ROUTING !Config Desc = Time step of th routing scheme !Config Def = 86400 !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. ! dt_routing = one_day CALL getin_p('ROUTING_TIMESTEP', dt_routing) ! !Config Key = ROUTING_RIVERS !Config If = RIVER_ROUTING !Config Desc = Number of rivers !Config Def = 50 !Config Help = This parameter chooses the number of largest river basins !Config which should be treated as independently as rivers and not !Config flow into the oceans as diffusion coastal flow. num_largest = 50 CALL getin_p('ROUTING_RIVERS', num_largest) ! !Config Key = DO_IRRIGATION !Config Desc = Should we compute an irrigation flux !Config Def = FALSE !Config Help = This parameters allows the user to ask the model !Config to compute an irigation flux. This performed for the !Config on very simple hypothesis. The idea is to have a good !Config map of irrigated areas and a simple function which estimates !Config the need to irrigate. ! doirrigation = .FALSE. CALL getin_p('DO_IRRIGATION', doirrigation) ! !Config Key = DO_FLOODPLAINS !Config Desc = Should we include floodplains !Config Def = FALSE !Config Help = This parameters allows the user to ask the model !Config to take into account the flood plains and return !Config the water into the soil moisture. It then can go !Config back to the atmopshere. This tried to simulate !Config internal deltas of rivers. ! dofloodplains = .FALSE. CALL getin_p('DO_FLOODPLAINS', dofloodplains) ! ! ! In order to simplify the time cascade check that dt_routing ! is a multiple of dtradia ! ratio = dt_routing/dtradia 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) * dtradia WRITE(numout,*) 'THE NEW DT_ROUTING IS : ', dt_routing ENDIF ! IF ( dt_routing .LT. dtradia) 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.' dt_routing = dtradia WRITE(numout,*) 'THE NEW DT_ROUTING IS : ', dt_routing ENDIF ! var_name ="routingcounter" IF (is_root_prc) THEN CALL ioconf_setatt('UNITS', 's') CALL ioconf_setatt('LONG_NAME','Time counter for the routing scheme') CALL restget (rest_id, var_name, 1, 1, 1, kjit, .TRUE., tmp_day) time_counter = tmp_day(1) CALL setvar (time_counter, val_exp, 'NO_KEYWORD', 0.0_r_std) ENDIF CALL bcast(time_counter) !!$ CALL setvar_p (time_counter, val_exp, 'NO_KEYWORD', 0.0_r_std) ! ALLOCATE (routing_area_loc(nbpt,nbasmax)) ALLOCATE (routing_area_glo(nbp_glo,nbasmax)) var_name = 'routingarea' IF (is_root_prc) THEN CALL ioconf_setatt('UNITS', 'm^2') CALL ioconf_setatt('LONG_NAME','Area of basin') CALL restget (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., routing_area_glo, "gather", nbp_glo, index_g) ENDIF CALL scatter(routing_area_glo,routing_area_loc) routing_area=>routing_area_loc ! ALLOCATE (tmp_real_g(nbp_glo,nbasmax)) ! ALLOCATE (route_togrid_loc(nbpt,nbasmax)) ALLOCATE (route_togrid_glo(nbp_glo,nbasmax)) ! used in global in routing_flow 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_flow CALL scatter(route_togrid_glo,route_togrid_loc) route_togrid=>route_togrid_loc ! ALLOCATE (route_tobasin_loc(nbpt,nbasmax)) ALLOCATE (route_tobasin_glo(nbp_glo,nbasmax)) IF (is_root_prc) THEN var_name = 'routetobasin' CALL ioconf_setatt('UNITS', '-') CALL ioconf_setatt('LONG_NAME','Basin in to which the water goes') CALL restget (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., tmp_real_g, "gather", nbp_glo, index_g) route_tobasin_glo = undef_int WHERE ( tmp_real_g .LT. val_exp ) route_tobasin_glo = NINT(tmp_real_g) ENDWHERE ENDIF CALL scatter(route_tobasin_glo,route_tobasin_loc) route_tobasin=>route_tobasin_loc ! ALLOCATE (global_basinid_loc(nbpt,nbasmax)) ALLOCATE (global_basinid_glo(nbp_glo,nbasmax)) IF (is_root_prc) THEN var_name = 'basinid' CALL ioconf_setatt('UNITS', '-') CALL ioconf_setatt('LONG_NAME','ID of basin') CALL restget (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., tmp_real_g, "gather", nbp_glo, index_g) global_basinid_glo = undef_int WHERE ( tmp_real_g .LT. val_exp ) global_basinid_glo = NINT(tmp_real_g) ENDWHERE ENDIF CALL scatter(global_basinid_glo,global_basinid_loc) global_basinid=>global_basinid_loc ! ALLOCATE (topo_resid_loc(nbpt,nbasmax)) ALLOCATE (topo_resid_glo(nbp_glo,nbasmax)) IF (is_root_prc) THEN var_name = 'topoindex' CALL ioconf_setatt('UNITS', 'm') CALL ioconf_setatt('LONG_NAME','Topographic index of the residence time') CALL restget (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., topo_resid_glo, "gather", nbp_glo, index_g) ENDIF CALL scatter(topo_resid_glo,topo_resid_loc) topo_resid=>topo_resid_loc ! ALLOCATE (fast_reservoir(nbpt,nbasmax)) var_name = 'fastres' CALL ioconf_setatt('UNITS', 'Kg') CALL ioconf_setatt('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', 0.0_r_std) ! ALLOCATE (slow_reservoir(nbpt,nbasmax)) var_name = 'slowres' CALL ioconf_setatt('UNITS', 'Kg') CALL ioconf_setatt('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', 0.0_r_std) ! ALLOCATE (stream_reservoir(nbpt,nbasmax)) var_name = 'streamres' CALL ioconf_setatt('UNITS', 'Kg') CALL ioconf_setatt('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', 0.0_r_std) ! ALLOCATE (lake_reservoir(nbpt)) var_name = 'lakeres' CALL ioconf_setatt('UNITS', 'Kg') CALL ioconf_setatt('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 (lake_reservoir, val_exp, 'NO_KEYWORD', 0.0_r_std) ! ! Map of irrigated areas ! IF ( doirrigation ) THEN ALLOCATE (irrigated(nbpt)) var_name = 'irrigated' CALL ioconf_setatt('UNITS', 'm^2') CALL ioconf_setatt('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 ! ALLOCATE (previous_outflow(nbpt, nbasmax+3)) var_name = 'previous_outflow' CALL ioconf_setatt('UNITS', 'm^3/dt') CALL ioconf_setatt('LONG_NAME','Previous outflow from this basin') CALL restget_p (rest_id, var_name, nbp_glo, nbasmax+3, 1, kjit, .TRUE., previous_outflow, "gather", nbp_glo, index_g) CALL setvar_p (previous_outflow, val_exp, 'NO_KEYWORD', 0.0_r_std) ! IF ( dofloodplains ) THEN ALLOCATE (floodplains(nbpt)) var_name = 'floodplains' CALL ioconf_setatt('UNITS', 'm^2') CALL ioconf_setatt('LONG_NAME','Surface which can be flooded') CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., floodplains, "gather", nbp_glo, index_g) CALL setvar_p (floodplains, val_exp, 'NO_KEYWORD', undef_sechiba) ENDIF ! ! Put into the restart file the fluxes so that they can be regenerated at restart. ! ALLOCATE (lakeinflow_mean(nbpt)) var_name = 'lakeinflow' CALL ioconf_setatt('UNITS', 'Kg/dt') CALL ioconf_setatt('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', 0.0_r_std) ! ALLOCATE (returnflow_mean(nbpt)) var_name = 'returnflow' CALL ioconf_setatt('UNITS', 'Kg/m^2/dt') CALL ioconf_setatt('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', 0.0_r_std) returnflow(:) = returnflow_mean(:) ! ALLOCATE (irrigation_mean(nbpt)) ALLOCATE (irrig_netereq(nbpt)) irrig_netereq(:) = zero ! IF ( doirrigation ) THEN var_name = 'irrigation' CALL ioconf_setatt('UNITS', 'Kg/dt') CALL ioconf_setatt('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', 0.0_r_std) irrigation(:) = irrigation_mean(:) ELSE irrigation_mean(:) = zero ENDIF ! ALLOCATE (riverflow_mean(nbpt)) var_name = 'riverflow' CALL ioconf_setatt('UNITS', 'Kg/dt') CALL ioconf_setatt('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', 0.0_r_std) riverflow(:) = riverflow_mean(:) ! ALLOCATE (coastalflow_mean(nbpt)) var_name = 'coastalflow' CALL ioconf_setatt('UNITS', 'Kg/dt') CALL ioconf_setatt('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', 0.0_r_std) coastalflow(:) = coastalflow_mean(:) ! ! Locate it at the 2m level ipn = MINLOC(ABS(diaglev-2)) floodtemp_lev = ipn(1) ALLOCATE (floodtemp(nbpt)) floodtemp(:) = stempdiag(:,floodtemp_lev) ! ALLOCATE(hydrographs(nbpt)) var_name = 'hydrographs' CALL ioconf_setatt('UNITS', 'm^3/dt') CALL ioconf_setatt('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', 0.0_r_std) ! ! The diagnostic variables, they are initialized from the above restart variables. ! ALLOCATE(fast_diag(nbpt), slow_diag(nbpt), stream_diag(nbpt), lake_diag(nbpt), stat=ier) ! fast_diag(:) = zero slow_diag(:) = zero stream_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) ENDDO ! fast_diag(ig) = fast_diag(ig)/totarea slow_diag(ig) = slow_diag(ig)/totarea stream_diag(ig) = stream_diag(ig)/totarea ! ! This is the volume of the lake scaled to the entire grid. ! It would be batter to scale it to the size of the lake ! but this information is not yet available. ! lake_diag(ig) = lake_reservoir(ig)/totarea ! ENDDO ! ! ! Get from the restart the fluxes we accumulated. ! ALLOCATE (runoff_mean(nbpt)) var_name = 'runoff_route' CALL ioconf_setatt('UNITS', 'Kg') CALL ioconf_setatt('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', 0.0_r_std) ! ALLOCATE(drainage_mean(nbpt)) var_name = 'drainage_route' CALL ioconf_setatt('UNITS', 'Kg') CALL ioconf_setatt('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', 0.0_r_std) ! ALLOCATE(evapot_mean(nbpt)) var_name = 'evapot_route' CALL ioconf_setatt('UNITS', 'Kg/m^2') CALL ioconf_setatt('LONG_NAME','Accumulated potential evaporation for routing') CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., evapot_mean, "gather", nbp_glo, index_g) CALL setvar_p (evapot_mean, val_exp, 'NO_KEYWORD', 0.0_r_std) ! ALLOCATE(precip_mean(nbpt)) var_name = 'precip_route' CALL ioconf_setatt('UNITS', 'Kg/m^2') CALL ioconf_setatt('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', 0.0_r_std) ! ALLOCATE(humrel_mean(nbpt)) var_name = 'humrel_route' CALL ioconf_setatt('UNITS', '-') CALL ioconf_setatt('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', 1.0_r_std) ! ALLOCATE(totnobio_mean(nbpt)) var_name = 'totnobio_route' CALL ioconf_setatt('UNITS', '-') CALL ioconf_setatt('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', 0.0_r_std) ! ALLOCATE(vegtot_mean(nbpt)) var_name = 'vegtot_route' CALL ioconf_setatt('UNITS', '-') CALL ioconf_setatt('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', 1.0_r_std) ! ! DEALLOCATE(tmp_real_g) ! ! Allocate diagnostic variables ! ALLOCATE(hydrodiag_loc(nbpt,nbasmax),hydrodiag_glo(nbp_glo,nbasmax),stat=ier) IF (ier.NE.0) THEN WRITE (numout,*) ' error in hydrodiag allocation. We stop. We need kjpindex words = ', nbpt*nbasmax STOP 'routing_init' END IF hydrodiag=>hydrodiag_loc ! ! END SUBROUTINE routing_init ! !--------------------------------------------------------------------------------- ! SUBROUTINE routing_clear() ! l_first_routing=.TRUE. ! 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(global_basinid_loc)) DEALLOCATE(global_basinid_loc) IF (ALLOCATED(topo_resid_loc)) DEALLOCATE(topo_resid_loc) IF (ALLOCATED(routing_area_glo)) DEALLOCATE(routing_area_glo) IF (ALLOCATED(route_togrid_glo)) DEALLOCATE(route_togrid_glo) IF (ALLOCATED(route_tobasin_glo)) DEALLOCATE(route_tobasin_glo) IF (ALLOCATED(global_basinid_glo)) DEALLOCATE(global_basinid_glo) IF (ALLOCATED(topo_resid_glo)) DEALLOCATE(topo_resid_glo) IF (ALLOCATED(fast_reservoir)) DEALLOCATE(fast_reservoir) IF (ALLOCATED(slow_reservoir)) DEALLOCATE(slow_reservoir) IF (ALLOCATED(stream_reservoir)) DEALLOCATE(stream_reservoir) IF (ALLOCATED(lake_reservoir)) DEALLOCATE(lake_reservoir) IF (ALLOCATED(returnflow_mean)) DEALLOCATE(returnflow_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(drainage_mean)) DEALLOCATE(drainage_mean) IF (ALLOCATED(evapot_mean)) DEALLOCATE(evapot_mean) IF (ALLOCATED(precip_mean)) DEALLOCATE(precip_mean) IF (ALLOCATED(humrel_mean)) DEALLOCATE(humrel_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(irrigation_mean)) DEALLOCATE(irrigation_mean) IF (ALLOCATED(irrigated)) DEALLOCATE(irrigated) IF (ALLOCATED(floodplains)) DEALLOCATE(floodplains) ! END SUBROUTINE routing_clear ! !--------------------------------------------------------------------------------- ! SUBROUTINE routing_flow(nbpt, dt_routing, runoff, drainage, & & vegtot, totnobio, evapot, precip, humrel, floodtemp, & & lakeinflow, returnflow, irrigation, riverflow, & & coastalflow, hydrographs) ! IMPLICIT NONE ! ! ! INPUT ! INTEGER(i_std), INTENT(in) :: nbpt !! Domain size REAL(r_std), INTENT (in) :: dt_routing !! Time step in seconds REAL(r_std), INTENT(in) :: runoff(nbpt) !! grid-point runoff REAL(r_std), INTENT(in) :: drainage(nbpt) !! grid-point drainage REAL(r_std), INTENT(in) :: vegtot(nbpt) !! Potentially vegetated area REAL(r_std), INTENT(in) :: totnobio(nbpt) !! Other areas whichcan not have vegetation REAL(r_std), INTENT(in) :: evapot(nbpt) !! grid-point potential evaporation REAL(r_std), INTENT(in) :: precip(nbpt) !! Precipiation (rainfall here) REAL(r_std), INTENT(in) :: humrel(nbpt) !! Soil moisture stress REAL(r_std), INTENT(in) :: floodtemp(nbpt) !! Temperature to decide if floodplains work REAL(r_std), INTENT(out) :: lakeinflow(nbpt) !! The water flow which flows into lakes (kg/dt) REAL(r_std), INTENT(out) :: returnflow(nbpt) !! Water flowing back into soil moisture (kg/m^2/dt) REAL(r_std), INTENT(out) :: irrigation(nbpt) !! The artificial irrigation (kg/m^2 per dt) REAL(r_std), INTENT(out) :: riverflow(nbpt) !! Outflow of the major rivers (kg/dt) REAL(r_std), INTENT(out) :: coastalflow(nbpt) !! Outflow on coastal points by small basins (kg/dt) REAL(r_std), INTENT(out) :: hydrographs(nbpt) !! Hydrographs at the outflow of the gird box for major basins ! ! LOCAL ! REAL(r_std), DIMENSION(nbpt, nbasmax) :: fast_flow REAL(r_std), DIMENSION(nbpt, nbasmax) :: slow_flow REAL(r_std), DIMENSION(nbpt, nbasmax) :: stream_flow REAL(r_std), DIMENSION(nbpt, nbasmax) :: reinfiltration REAL(r_std), DIMENSION(nbpt, nbasmax) :: baseirrig !! Irrigation uptake from each basin reservoir. REAL(r_std), DIMENSION(nbpt, 0:nbasmax+3) :: transport REAL(r_std), DIMENSION(nbp_glo, 0:nbasmax+3) :: transport_glo REAL(r_std), DIMENSION(nbp_glo, 0:nbasmax+3) :: transport_sum REAL(r_std), DIMENSION(nbpt, nbasmax) :: floods, wdelay, previous_outflow, potflood, inflow REAL(r_std), DIMENSION(nbpt) :: tobeflooded REAL(r_std), DIMENSION(nbpt) :: totarea REAL(r_std) :: flow, floodindex INTEGER(i_std) :: ig, ib, rtg, rtb, ierr REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: fast_flow_g,slow_flow_g,stream_flow_g,floods_g,wdelay_g ! transport(:,:) = zero transport_glo(:,:) = zero previous_outflow(:,:) = zero irrig_netereq(:) = zero baseirrig(:,:) = zero totarea(:) = zero ! ! Compute all the fluxes ! DO ig=1,nbpt DO ib=1,nbasmax ! totarea(ig) = totarea(ig) + routing_area(ig,ib) ! IF ( route_tobasin(ig,ib) .GT. 0 ) THEN ! flow = MIN(fast_reservoir(ig,ib)/((topo_resid(ig,ib)/1000.)*fast_tcst*one_day/dt_routing),& & fast_reservoir(ig,ib)) fast_flow(ig,ib) = flow ! flow = MIN(slow_reservoir(ig,ib)/((topo_resid(ig,ib)/1000.)*slow_tcst*one_day/dt_routing),& & slow_reservoir(ig,ib)) slow_flow(ig,ib) = flow ! flow = MIN(stream_reservoir(ig,ib)/((topo_resid(ig,ib)/1000.)*stream_tcst*one_day/dt_routing),& & stream_reservoir(ig,ib)) stream_flow(ig,ib) = flow ! ! ELSE fast_flow(ig,ib) = 0.0 slow_flow(ig,ib) = 0.0 stream_flow(ig,ib) = 0.0 ENDIF inflow(ig,ib) = fast_flow(ig,ib) + slow_flow(ig,ib) + stream_flow(ig,ib) ENDDO ENDDO ! ! Do the floodings ! IF ( dofloodplains ) THEN DO ig=1,nbpt tobeflooded(ig) = floodplains(ig) DO ib=1,nbasmax potflood(ig,ib) = inflow(ig,ib) - previous_outflow(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 ! floods(ig,ib) = floodindex * evap_cst * potflood(ig,ib) ! wdelay(ig,ib) = floodindex * wdelay_cst * potflood(ig,ib) ! ! tobeflooded(ig) = tobeflooded(ig) - routing_area(ig,ib) ! ELSE floods(ig,ib) = zero wdelay(ig,ib) = zero ENDIF ENDDO ENDDO ELSE DO ig=1,nbpt DO ib=1,nbasmax floods(ig,ib) = zero wdelay(ig,ib) = zero ENDDO ENDDO ENDIF ! DO ig=1,nbpt DO ib=1,nbasmax reinfiltration(ig,ib) = floods(ig,ib) ENDDO ENDDO !ym cette methode conserve les erreurs d'arrondie !ym mais n'est pas la plus efficace IF (is_root_prc) & ALLOCATE( fast_flow_g(nbp_glo, nbasmax), slow_flow_g(nbp_glo, nbasmax), & stream_flow_g(nbp_glo, nbasmax), floods_g(nbp_glo, nbasmax), wdelay_g(nbp_glo, nbasmax) ) CALL gather(fast_flow,fast_flow_g) CALL gather(slow_flow,slow_flow_g) CALL gather(stream_flow,stream_flow_g) CALL gather(floods,floods_g) CALL gather(wdelay,wdelay_g) IF (is_root_prc) THEN DO ig=1,nbp_glo DO ib=1,nbasmax ! rtg = route_togrid_glo(ig,ib) rtb = route_tobasin_glo(ig,ib) transport_glo(rtg,rtb) = transport_glo(rtg,rtb) + fast_flow_g(ig,ib) + slow_flow_g(ig,ib) + & & stream_flow_g(ig,ib) - floods_g(ig,ib) - wdelay_g(ig,ib) ! ENDDO ENDDO ENDIF IF (is_root_prc) & DEALLOCATE( fast_flow_g, slow_flow_g, stream_flow_g, floods_g, wdelay_g ) CALL scatter(transport_glo,transport) !ym DO ig=1,nbpt !ym DO ib=1,nbasmax !ym ! !ym rtg = route_togrid(ig,ib) !ym rtb = route_tobasin(ig,ib) !ym transport_glo(rtg,rtb) = transport_glo(rtg,rtb) + fast_flow(ig,ib) + slow_flow(ig,ib) + & !ym & stream_flow(ig,ib) - floods(ig,ib) - wdelay(ig,ib) !ym ! !ym ENDDO !ym ENDDO !ym !ym CALL reduce_sum(transport_glo,transport_sum) !ym CALL scatter(transport_sum,transport) ! ! Update all reservoirs ! DO ig=1,nbpt DO ib=1,nbasmax fast_reservoir(ig,ib) = fast_reservoir(ig,ib) + runoff(ig)*routing_area(ig,ib) - & & reinfiltration(ig,ib) - fast_flow(ig,ib) + floods(ig,ib) + wdelay(ig,ib) slow_reservoir(ig,ib) = slow_reservoir(ig,ib) + drainage(ig)*routing_area(ig,ib) - & & slow_flow(ig,ib) stream_reservoir(ig,ib) = stream_reservoir(ig,ib) + transport(ig,ib) - & & stream_flow(ig,ib) ! previous_outflow(ig,ib) = fast_flow(ig,ib) + slow_flow(ig,ib) + stream_flow(ig,ib) - & & wdelay(ig,ib) - floods(ig,ib) ! ENDDO ENDDO ! ! Compute the total reinfiltration in the grid box ! returnflow(:) = zero DO ig=1,nbpt ! DO ib=1,nbasmax returnflow(ig) = returnflow(ig) + reinfiltration(ig,ib) ENDDO ! returnflow(ig) = returnflow(ig)/totarea(ig) ! ENDDO ! ! 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. ! IF ( doirrigation ) THEN DO ig=1,nbpt IF ((vegtot(ig) .GT. 0.0) .AND. (humrel(ig) .LT. 0.99)) THEN irrig_netereq(ig) = (irrigated(ig) / totarea(ig) ) * MAX(0.0, & & crop_coef * evapot(ig) - & & MAX(precip(ig)+returnflow(ig)-runoff(ig)-drainage(ig), zero) ) irrig_netereq(ig) = 1 * irrig_netereq(ig) IF(irrig_netereq(ig).LT.0.0) THEN WRITE(numout,*) 'there is a probleme for irrig_netereq',ig,irrig_netereq(ig) ENDIF ENDIF DO ib=1,nbasmax IF ( route_tobasin(ig,ib) .GT. 0 ) THEN baseirrig(ig,ib) = MIN( irrig_netereq(ig) * routing_area(ig,ib),& & stream_reservoir(ig,ib) + fast_reservoir(ig,ib) + slow_reservoir(ig,ib) ) slow_reservoir(ig,ib) = MAX(0.0, slow_reservoir(ig,ib) + & & MIN(0.0, fast_reservoir(ig,ib) + MIN(0.0, stream_reservoir(ig,ib)-baseirrig(ig,ib)))) fast_reservoir(ig,ib) = MAX( 0.0, & & fast_reservoir(ig,ib) + MIN(0.0, stream_reservoir(ig,ib)-baseirrig(ig,ib))) stream_reservoir(ig,ib) = MAX(0.0, stream_reservoir(ig,ib)-baseirrig(ig,ib) ) IF(baseirrig(ig,ib) .LT. 0.0 .OR. slow_reservoir(ig,ib) .LT. 0.0 .OR. & & fast_reservoir(ig,ib) .LT. 0.0 .OR. stream_reservoir(ig,ib) .LT. 0.0) THEN WRITE(numout,*) 'There is negative values related to irrigation', ig,ib,baseirrig(ig,ib), & & slow_reservoir(ig,ib),fast_reservoir(ig,ib),stream_reservoir(ig,ib) ENDIF ENDIF previous_outflow(ig,ib) = fast_flow(ig,ib) + slow_flow(ig,ib) + stream_flow(ig,ib) - & & wdelay(ig,ib) - floods(ig,ib) ENDDO ENDDO ! ENDIF ! ! ! Compute the fluxes which leave the routing scheme ! ! Lakeinflow is in Kg/dt ! returnflow is in Kg/m^2/dt ! hydrographs(:) = zero fast_diag(:) = zero slow_diag(:) = zero stream_diag(:) = zero irrigation(:) = zero ! ! DO ig=1,nbpt ! DO ib=1,nbasmax IF (hydrodiag(ig,ib) > 0 ) THEN hydrographs(ig) = hydrographs(ig) + fast_flow(ig,ib) + slow_flow(ig,ib) + & & stream_flow(ig,ib) - floods(ig,ib) - wdelay(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) irrigation (ig) = irrigation (ig) + baseirrig(ig,ib) ENDDO ! fast_diag(ig) = fast_diag(ig)/totarea(ig) slow_diag(ig) = slow_diag(ig)/totarea(ig) stream_diag(ig) = stream_diag(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) ! IF ( irrigation(ig) .GE. irrig_netereq(ig)+1e4 ) THEN WRITE(numout,*) 'There is a problem here with irrigation',ig,irrigation(ig),irrig_netereq(ig) WRITE(numout,*) irrigated(ig),totarea(ig), evapot(ig), precip_mean(ig),runoff(ig),drainage(ig) STOP ENDIF ! ENDDO ! END SUBROUTINE routing_flow ! !--------------------------------------------------------------------------------- ! SUBROUTINE routing_lake(nbpt, dt_routing, lakeinflow, humrel, returnflow) ! IMPLICIT NONE ! ! This routine 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. ! The return flow to the soil moisture reservoir is based on a maximum ! lake evaporation rate (maxevap_lake). ! INTEGER(i_std), INTENT(in) :: nbpt !! Domain size REAL(r_std), INTENT (in) :: dt_routing !! Time step in seconds REAL(r_std), INTENT(in) :: lakeinflow(nbpt) !! Flow into the lake (kg/dt) REAL(r_std), INTENT(in) :: humrel(nbpt) !! Soil moisture deficit around the lake (Hum !) REAL(r_std), INTENT(out) :: returnflow(nbpt) !! Water flowing back into soil moisture (kg/m^2/dt) ! ! LOCAL ! INTEGER(i_std) :: ig REAL(r_std) :: refill, total_area ! ! ! DO ig=1,nbpt ! total_area = SUM(routing_area(ig,:)) ! lake_reservoir(ig) = lake_reservoir(ig) + lakeinflow(ig) !uptake in Kg/dt refill = MAX(zero, maxevap_lake * (un - humrel(ig)) * dt_routing * total_area) returnflow(ig) = MIN(refill, lake_reservoir(ig)) lake_reservoir(ig) = lake_reservoir(ig) - returnflow(ig) !Return in Kg/m^2/dt returnflow(ig) = returnflow(ig)/total_area ! ! This is the volume of the lake scaled to the entire grid. ! It would be batter to scale it to the size of the lake ! but this information is not yet available. lake_diag(ig) = lake_reservoir(ig)/total_area ! ENDDO ! END SUBROUTINE routing_lake ! !--------------------------------------------------------------------------------- ! SUBROUTINE routing_diagnostic_p(nbpt,index, resolution, contfrac, hist_id, hist2_id) ! IMPLICIT NONE INTEGER(i_std), INTENT(in) :: nbpt !! Domain size INTEGER(i_std), INTENT(in) :: index(nbpt) !! Indeces of the points on the map REAL(r_std), INTENT(in) :: resolution(nbpt,2) !! The size in km of each grid-box in X and Y REAL(r_std), INTENT(in) :: contfrac(nbpt) !! Fraction of land in each grid box. INTEGER(i_std),INTENT (in) :: hist_id !! _history_ file identifier INTEGER(i_std),INTENT (in) :: hist2_id !! _history_ file 2 identifier REAL(r_std), DIMENSION(nbpt) :: nbrivers !! Number of rivers in the grid REAL(r_std), DIMENSION(nbpt) :: basinmap !! Map of basins REAL(r_std), DIMENSION(nbp_glo) :: nbrivers_g !! Number of rivers in the grid REAL(r_std), DIMENSION(nbp_glo) :: basinmap_g !! Map of basins routing_area => routing_area_glo topo_resid => topo_resid_glo route_togrid => route_togrid_glo route_tobasin => route_tobasin_glo global_basinid => global_basinid_glo hydrodiag=>hydrodiag_glo IF (is_root_prc) CALL routing_diagnostic(nbp_glo,index_g, resolution_g, contfrac_g, nbrivers_g,basinmap_g) routing_area => routing_area_loc topo_resid => topo_resid_loc route_togrid => route_togrid_loc route_tobasin => route_tobasin_loc global_basinid => global_basinid_loc hydrodiag=>hydrodiag_loc CALL scatter(nbrivers_g,nbrivers) CALL scatter(basinmap_g,basinmap) CALL scatter(hydrodiag_glo,hydrodiag_loc) IF ( .NOT. almaoutput ) THEN CALL histwrite(hist_id, 'basinmap', 1, basinmap, nbpt, index) CALL histwrite(hist_id, 'nbrivers', 1, nbrivers, nbpt, index) ELSE ENDIF IF ( hist2_id > 0 ) THEN IF ( .NOT. almaoutput ) THEN CALL histwrite(hist2_id, 'basinmap', 1, basinmap, nbpt, index) CALL histwrite(hist2_id, 'nbrivers', 1, nbrivers, nbpt, index) ELSE ENDIF ENDIF END SUBROUTINE routing_diagnostic_p SUBROUTINE routing_diagnostic(nbpt,index, resolution, contfrac,nbrivers,basinmap) ! ! This subroutine will set up a map of the major basins ! ! INPUTS ! INTEGER(i_std), INTENT(in) :: nbpt !! Domain size INTEGER(i_std), INTENT(in) :: index(nbpt) !! Indeces of the points on the map REAL(r_std), INTENT(in) :: resolution(nbpt,2) !! The size in km of each grid-box in X and Y REAL(r_std), INTENT(in) :: contfrac(nbpt) !! Fraction of land in each grid box. ! ! OUTPUTS ! REAL(r_std), DIMENSION(nbpt), INTENT(out) :: nbrivers !! Number of rivers in the grid REAL(r_std), DIMENSION(nbpt), INTENT(out) :: basinmap !! Map of basins ! ! LOCAL ! INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:) :: pts !! list the points belonging to the basin INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:) :: ptbas !! list the basin number for this point INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:) :: outpt !! Outflow point for each basin INTEGER(i_std), ALLOCATABLE, DIMENSION(:) :: nb_pts !! Number of points in the basin REAL(r_std), ALLOCATABLE, DIMENSION(:) :: totarea !! Total area of basin INTEGER(i_std), ALLOCATABLE, DIMENSION(:) :: topids !! The IDs of the first num_largest basins. CHARACTER(LEN=25), ALLOCATABLE, DIMENSION(:) :: basin_names !! Names of the rivers CHARACTER(LEN=25) :: name_str ! INTEGER(i_std) :: ig, ib, og, ob, ign, ibn, ff(1), ic, icc, nb_small ! ALLOCATE(pts(num_largest, nbpt)) ALLOCATE(ptbas(num_largest, nbpt)) ALLOCATE(outpt(num_largest, 2)) ALLOCATE(nb_pts(num_largest)) ALLOCATE(totarea(num_largest)) ALLOCATE(topids(num_largest)) ! ! ! First we get the list of all river outflow points ! We work under the assumption that we only have num_largest basins finishing with ! nbasmax+3. This is checked in routing_truncate. ! nb_small = 1 outpt(:,:) = -1 ic = 0 DO ig=1,nbpt DO ib=1,nbasmax ign = route_togrid(ig, ib) ibn = route_tobasin(ig, ib) IF ( ibn .EQ. nbasmax+3) THEN ic = ic + 1 outpt(ic,1) = ig outpt(ic,2) = ib ! ! Get the largest id of the basins we call a river. This is ! to extract the names of all rivers. ! IF ( global_basinid(ig,ib) > nb_small ) THEN nb_small = global_basinid(ig,ib) ENDIF ENDIF ENDDO ENDDO ! nb_small = MIN(nb_small, 349) ! ALLOCATE(basin_names(nb_small)) ! CALL routing_names(nb_small, basin_names) ! ! Go through all points and basins to see if they outflow as a river and store the ! information needed in the various arrays. ! nb_pts(:) = 0 totarea(:) = 0.0 hydrodiag(:,:) = 0 DO ig=1,nbpt DO ib=1,nbasmax ic = 0 ign = ig ibn = ib ! Locate outflow point DO WHILE (ibn .GT. 0 .AND. ibn .LE. nbasmax .AND. ic .LT. nbasmax*nbpt) ic = ic + 1 og = ign ob = ibn ign = route_togrid(og, ob) ibn = route_tobasin(og, ob) ENDDO ! ! Now that we have an outflow check if it is one of the num_largest rivers. ! In this case we keeps the location so we diagnose it. ! IF ( ibn .EQ. nbasmax + 3) THEN DO icc = 1,num_largest IF ( outpt(icc,1) .EQ. og .AND. outpt(icc,2) .EQ. ob ) THEN ! ! We only keep this point for our map if it is large enough. ! IF ( routing_area(ig,ib) .GT. 0.25*resolution(ig,1)*resolution(ig,2)*contfrac(ig) ) THEN nb_pts(icc) = nb_pts(icc) + 1 pts(icc, nb_pts(icc)) = ig ptbas(icc, nb_pts(icc)) = ib ENDIF totarea(icc) = totarea(icc) + routing_area(ig,ib) ! ID of the river is taken from the last point before the outflow. topids(icc) = global_basinid(og,ob) ! ! On this gridbox and basin we will diagnose the hydrograph ! hydrodiag(ig, ib) = 1 ! ENDIF ENDDO ENDIF ENDDO ENDDO ! ! Construct the map of the largest basins. We take the num_largest basins ! if they have more than 4 points. After that it is of no use. ! ! basinmap(:) = 0.0 DO icc = 1, num_largest ff = MAXLOC(totarea) IF ( nb_pts(ff(1)) .GT. 2 ) THEN DO ig = 1, nb_pts(ff(1)) basinmap(pts(ff(1),ig)) = REAL(icc,r_std) ENDDO ! IF ( topids(ff(1)) .GT. nb_small ) THEN WRITE(name_str, '("NN, Nb : ",I4)') topids(ff(1)) ELSE name_str = basin_names(topids(ff(1))) ENDIF ! WRITE(numout,& '("Basin ID ", I5," ", A15, " Area [km^2] : ", F13.4, " Nb points : ", I4)')& & topids(ff(1)), name_str(1:15), totarea(ff(1))/1.e6, nb_pts(ff(1)) ENDIF totarea(ff(1)) = 0.0 ENDDO ! ! nbrivers(:) = zero DO ig=1,nbpt nbrivers(ig) = COUNT(route_tobasin(ig,1:nbasmax) == nbasmax+3) ENDDO DO ig=1,nbpt IF ( nbrivers(ig) > 1 ) THEN WRITE(numout,*) 'Grid box ', ig, ' has ', NINT(nbrivers(ig)), ' outflow points.' WRITE(numout,*) 'The rivers which flow into the ocean at this point are :' DO icc=1,nbasmax IF ( route_tobasin(ig,icc) == nbasmax+3) THEN IF ( global_basinid(ig,icc) <= nb_small ) THEN WRITE(numout,*) 'ID = ',global_basinid(ig,icc), ' Name = ', basin_names(global_basinid(ig,icc)) ELSE WRITE(numout,*) 'ID = ',global_basinid(ig,icc), ' Problem ===== ID is larger than possible' ENDIF ENDIF ENDDO ENDIF ENDDO ! WRITE(numout,*) 'Maximum topographic index :', MAXVAL(topo_resid) ic = COUNT(topo_resid .GT. 0.) WRITE(numout,*) 'Mean topographic index :', SUM(topo_resid, MASK=topo_resid .GT. 0.)/ic WRITE(numout,*) 'Minimum topographic index :', MINVAL(topo_resid, MASK=topo_resid .GT. 0.) ! DEALLOCATE(pts) DEALLOCATE(outpt) DEALLOCATE(nb_pts) DEALLOCATE(totarea) ! END SUBROUTINE routing_diagnostic ! !--------------------------------------------------------------------------------- ! SUBROUTINE routing_basins_p(nbpt, lalo, neighbours, resolution, contfrac) ! IMPLICIT NONE ! INTEGER(i_std), INTENT(in) :: nbpt ! Number of points for which the data needs to be interpolated REAL(r_std), INTENT(in) :: lalo(nbpt,2) ! Vector of latitude and longitudes (beware of the order !) INTEGER(i_std), INTENT(in) :: neighbours(nbpt,8) ! Vector of neighbours for each grid point (1=N, 2=E, 3=S, 4=W) REAL(r_std), INTENT(in) :: resolution(nbpt,2) ! The size in km of each grid-box in X and Y REAL(r_std), INTENT(in) :: contfrac(nbpt) ! Fraction of land in each grid box. ! INTEGER(i_std) :: neighbours_tmp(nbpt,8) INTEGER(i_std) :: i,j ! DO i=1,nbp_loc ! DO j=1,8 ! IF (neighbours(i,j)==-1) THEN ! neighbours_tmp(i,j)=neighbours(i,j) ! ELSE ! neighbours_tmp(i,j)=neighbours(i,j)+nbp_para_begin(mpi_rank)-1 ! ENDIF ! ENDDO ! ENDDO routing_area => routing_area_glo topo_resid => topo_resid_glo route_togrid => route_togrid_glo route_tobasin => route_tobasin_glo global_basinid => global_basinid_glo IF (is_root_prc) CALL routing_basins(nbp_glo,lalo_g, neighbours_g, resolution_g, contfrac_g) routing_area => routing_area_loc topo_resid => topo_resid_loc route_togrid => route_togrid_loc route_tobasin => route_tobasin_loc global_basinid => global_basinid_loc CALL scatter(routing_area_glo,routing_area_loc) CALL scatter(topo_resid_glo,topo_resid_loc) CALL scatter(route_togrid_glo,route_togrid_loc) CALL scatter(route_tobasin_glo,route_tobasin_loc) CALL scatter(global_basinid_glo,global_basinid_loc) END SUBROUTINE routing_basins_p SUBROUTINE routing_basins(nbpt, lalo, neighbours, resolution, contfrac) ! IMPLICIT NONE ! ! ! This subroutine reads in the map of basins and flow direction to construct the ! the catchments of each grid box. ! ! The work is done in a number of steps which are performed locally on the ! GCM grid: ! 1) First we find the grid-points of the high resolution routing grid which are ! within the coarser grid of GCM. ! 2) When we have these grid points we decompose them into basins in the routine ! routing_findbasins. A number of simplifications are done if needed. ! 3) In the routine routing_globalize we put the basin information of this grid ! into global fields. ! Then we work on the global grid to perform the following tasks : ! 1) We linkup the basins of the various and check the global consistence. ! 2) The area of each outflow point is computed. ! 3) The final step is to reduce the number of basins in order to fit into the truncation. ! ! ! 0.1 INPUT ! INTEGER(i_std), INTENT(in) :: nbpt ! Number of points for which the data needs to be interpolated REAL(r_std), INTENT(in) :: lalo(nbpt,2) ! Vector of latitude and longitudes (beware of the order !) INTEGER(i_std), INTENT(in) :: neighbours(nbpt,8) ! Vector of neighbours for each grid point (1=N, 2=E, 3=S, 4=W) REAL(r_std), INTENT(in) :: resolution(nbpt,2) ! The size in km of each grid-box in X and Y REAL(r_std), INTENT(in) :: contfrac(nbpt) ! Fraction of land in each grid box. ! ! ! 0.3 LOCAL ! REAL(r_std), PARAMETER :: R_Earth = 6378000. ! CHARACTER(LEN=80) :: filename INTEGER(i_std) :: iml, jml, lml, tml, fid, ib, ip, jp, fopt, lastjp, nbexp REAL(r_std) :: lev(1), date, dt, coslat, pi INTEGER(i_std) :: itau(1), sgn REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: trip, basins, topoindex, hierarchy REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: lat_rel, lon_rel, lat_ful, lon_ful REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: loup_rel, lolow_rel, laup_rel, lalow_rel REAL(r_std) :: lon_up, lon_low, lat_up, lat_low ! INTEGER(i_std) :: nbi, nbj REAL(r_std) :: ax, ay, min_topoind, max_basins, invented_basins REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: sub_area INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:,:) :: sub_index INTEGER(i_std), ALLOCATABLE, DIMENSION(:) :: sub_pts REAL(r_std), DIMENSION(nbvmax,nbvmax) :: area_bx, hierarchy_bx, lon_bx, lat_bx, topoind_bx INTEGER(i_std), DIMENSION(nbvmax,nbvmax) :: trip_bx, basin_bx ! INTEGER(i_std) :: coast_pts(nbvmax) REAL(r_std) :: lonrel, louprel, lolowrel ! INTEGER(i_std), ALLOCATABLE, DIMENSION(:) :: basin_count INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:) :: basin_id REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: basin_area, basin_hierarchy, basin_topoind REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: fetch_basin INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:) :: basin_flowdir INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:) :: outflow_grid INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:) :: outflow_basin INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:) :: inflow_number INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:,:) :: inflow_basin, inflow_grid INTEGER(i_std), ALLOCATABLE, DIMENSION(:) :: nbcoastal INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:) :: coastal_basin ! INTEGER(i_std) :: nb_basin, nwbas INTEGER(i_std) :: basin_inbxid(nbvmax), basin_sz(nbvmax), basin_pts(nbvmax,nbvmax,2), basin_bxout(nbvmax) CHARACTER(LEN=7) :: fmt LOGICAL :: debug = .FALSE. INTEGER(i_std), DIMENSION(2) :: diagbox = (/ 1147, 1148 /) ! Test on diagbox and nbpt IF (debug) THEN IF (ANY(diagbox .GT. nbpt)) THEN WRITE(*,*) "Debug diganostics : nbpt, diagbox", nbpt, diagbox call ipslerr(3,'routing_basin', & & 'Problem with diagbox in debug mode.', & & 'diagbox values can''t be greater than land points number.', & & '(decrease diagbox wrong value)') ENDIF ENDIF ! pi = 4. * ATAN(1.) ! ! Needs to be a configurable variable ! ! !Config Key = ROUTING_FILE !Config Desc = Name of file which contains the routing information !Config Def = routing.nc !Config Help = The file provided here should alow the routing module to !Config read the high resolution grid of basins and the flow direction !Config from one mesh to the other. ! filename = 'routing.nc' CALL getin('ROUTING_FILE',filename) ! CALL flininfo(filename,iml, jml, lml, tml, fid) ! ! ALLOCATE (lat_rel(iml,jml)) ALLOCATE (lon_rel(iml,jml)) ALLOCATE (laup_rel(iml,jml)) ALLOCATE (loup_rel(iml,jml)) ALLOCATE (lalow_rel(iml,jml)) ALLOCATE (lolow_rel(iml,jml)) ALLOCATE (lat_ful(iml+2,jml+2)) ALLOCATE (lon_ful(iml+2,jml+2)) ALLOCATE (trip(iml,jml)) ALLOCATE (basins(iml,jml)) ALLOCATE (topoindex(iml,jml)) ALLOCATE (hierarchy(iml,jml)) ! ALLOCATE (sub_area(nbpt,nbvmax)) ALLOCATE (sub_index(nbpt,nbvmax,2)) ALLOCATE (sub_pts(nbpt)) ! ! CALL flinopen(filename, .FALSE., iml, jml, lml, lon_rel, lat_rel, lev, tml, itau, date, dt, fid) ! ! The trip field follows the following convention for the flow of the water : ! trip = 1 : flow = N ! trip = 2 : flow = NE ! trip = 3 : flow = E ! trip = 4 : flow = SE ! trip = 5 : flow = S ! trip = 6 : flow = SW ! trip = 7 : flow = W ! trip = 8 : flow = NW ! trip = 97 : return flow into the ground ! trip = 98 : coastal flow (diffuse flow into the oceans) ! trip = 99 : river flow into the oceans ! ! CALL flinget(fid, 'trip', iml, jml, lml, tml, 1, 1, trip) ! CALL flinget(fid, 'basins', iml, jml, lml, tml, 1, 1, basins) ! CALL flinget(fid, 'topoind', iml, jml, lml, tml, 1, 1, topoindex) ! CALL flinclo(fid) ! nbexp = 0 ! min_topoind = MINVAL(topoindex, MASK=topoindex .LT. undef_sechiba-1.) ! DO ip=1,iml DO jp=1,jml IF ( trip(ip,jp) .LT. 1.e10 .AND. topoindex(ip,jp) .GT. 1.e10) THEN WRITE(numout,*) 'trip exists but not topoind :' WRITE(numout,*) 'ip, jp :', ip, jp WRITE(numout,*) 'trip, topoind : ', trip(ip,jp), topoindex(ip,jp) STOP ENDIF ENDDO ENDDO ! ! Duplicate the border assuming we have a global grid going from west to east ! lon_ful(2:iml+1,2:jml+1) = lon_rel(1:iml,1:jml) lat_ful(2:iml+1,2:jml+1) = lat_rel(1:iml,1:jml) ! IF ( lon_rel(iml,1) .LT. lon_ful(2,2)) THEN lon_ful(1,2:jml+1) = lon_rel(iml,1:jml) lat_ful(1,2:jml+1) = lat_rel(iml,1:jml) ELSE lon_ful(1,2:jml+1) = lon_rel(iml,1:jml)-360 lat_ful(1,2:jml+1) = lat_rel(iml,1:jml) ENDIF IF ( lon_rel(1,1) .GT. lon_ful(iml+1,2)) THEN lon_ful(iml+2,2:jml+1) = lon_rel(1,1:jml) lat_ful(iml+2,2:jml+1) = lat_rel(1,1:jml) ELSE lon_ful(iml+2,2:jml+1) = lon_rel(1,1:jml)+360 lat_ful(iml+2,2:jml+1) = lat_rel(1,1:jml) ENDIF ! sgn = INT(lat_rel(1,1)/ABS(lat_rel(1,1))) lat_ful(2:iml+1,1) = sgn*180 - lat_rel(1:iml,1) sgn = INT(lat_rel(1,jml)/ABS(lat_rel(1,jml))) lat_ful(2:iml+1,jml+2) = sgn*180 - lat_rel(1:iml,jml) lat_ful(1,1) = lat_ful(iml+1,1) lat_ful(iml+2,1) = lat_ful(2,1) lat_ful(1,jml+2) = lat_ful(iml+1,jml+2) lat_ful(iml+2,jml+2) = lat_ful(2,jml+2) ! ! Add the longitude lines to the top and bottom ! lon_ful(:,1) = lon_ful(:,2) lon_ful(:,jml+2) = lon_ful(:,jml+1) ! ! Get the upper and lower limits of each grid box ! DO ip=1,iml DO jp=1,jml ! loup_rel(ip,jp) =MAX(0.5*(lon_ful(ip,jp+1)+lon_ful(ip+1,jp+1)),& & 0.5*(lon_ful(ip+1,jp+1)+lon_ful(ip+2,jp+1))) lolow_rel(ip,jp) =MIN(0.5*(lon_ful(ip,jp+1)+lon_ful(ip+1,jp+1)),& & 0.5*(lon_ful(ip+1,jp+1)+lon_ful(ip+2,jp+1))) laup_rel(ip,jp) =MAX(0.5*(lat_ful(ip+1,jp)+lat_ful(ip+1,jp+1)),& & 0.5*(lat_ful(ip+1,jp+1)+lat_ful(ip+1,jp+2))) lalow_rel(ip,jp) =MIN(0.5*(lat_ful(ip+1,jp)+lat_ful(ip+1,jp+1)),& & 0.5*(lat_ful(ip+1,jp+1)+lat_ful(ip+1,jp+2))) ! ENDDO ENDDO ! ! ! ! Now we take each grid point and find out which values from the forcing we need to average ! DO ib =1, nbpt ! ! We find the 4 limits of the grid-box. As we transform the resolution of the model ! into longitudes and latitudes we do not have the problem of periodicity. ! coslat is a help variable here ! ! coslat = MAX(COS(lalo(ib,1) * pi/180. ), 0.001 )*pi/180. * R_Earth ! lon_up = lalo(ib,2) + resolution(ib,1)/(2.0*coslat) lon_low =lalo(ib,2) - resolution(ib,1)/(2.0*coslat) ! coslat = pi/180. * R_Earth ! lat_up =lalo(ib,1) + resolution(ib,2)/(2.0*coslat) lat_low =lalo(ib,1) - resolution(ib,2)/(2.0*coslat) ! ! ! Find the grid boxes from the data that go into the model's boxes ! We still work as if we had a regular grid ! Well it needs to be localy regular so ! so that the longitude at the latitude of the last found point is close to the one of the next point. ! fopt = 0 lastjp = 1 DO ip=1,iml ! ! Either the center of the data grid point is in the interval of the model grid or ! the East and West limits of the data grid point are on either sides of the border of ! the data grid. ! ! To do that correctly we have to check if the grid box sits on the date-line. ! IF ( lon_low < -180.0 ) THEN lonrel = MOD( lon_rel(ip,lastjp) - 360.0, 360.0) lolowrel = MOD( lolow_rel(ip,lastjp) - 360.0, 360.0) louprel = MOD( loup_rel(ip,lastjp) - 360.0, 360.0) ! ELSE IF ( lon_up > 180.0 ) THEN lonrel = MOD( 360. - lon_rel(ip,lastjp), 360.0) lolowrel = MOD( 360. - lolow_rel(ip,lastjp), 360.0) louprel = MOD( 360. - loup_rel(ip,lastjp), 360.0) ELSE lonrel = lon_rel(ip,lastjp) lolowrel = lolow_rel(ip,lastjp) louprel = loup_rel(ip,lastjp) ENDIF ! ! ! IF ( lonrel > lon_low .AND. lonrel < lon_up .OR. & & lolowrel < lon_low .AND. louprel > lon_low .OR. & & lolowrel < lon_up .AND. louprel > lon_up ) THEN ! DO jp = 1, jml ! ! Now that we have the longitude let us find the latitude ! IF ( lat_rel(ip,jp) > lat_low .AND. lat_rel(ip,jp) < lat_up .OR. & & lalow_rel(ip,jp) < lat_low .AND. laup_rel(ip,jp) > lat_low .OR.& & lalow_rel(ip,jp) < lat_up .AND. laup_rel(ip,jp) > lat_up) THEN ! lastjp = jp ! ! Is it a land point ? ! IF (trip(ip,jp) .LT. 1.e10) THEN ! fopt = fopt + 1 IF ( fopt .GT. nbvmax) THEN WRITE(numout,*) 'Please increase nbvmax in subroutine routing_basins', ib STOP ELSE ! ! If we sit on the date line we need to do the same transformations as above. ! IF ( lon_low < -180.0 ) THEN lolowrel = MOD( lolow_rel(ip,jp) - 360.0, 360.0) louprel = MOD( loup_rel(ip,jp) - 360.0, 360.0) ! ELSE IF ( lon_up > 180.0 ) THEN lolowrel = MOD( 360. - lolow_rel(ip,jp), 360.0) louprel = MOD( 360. - loup_rel(ip,jp), 360.0) ELSE lolowrel = lolow_rel(ip,jp) louprel = loup_rel(ip,jp) ENDIF ! ! Get the area of the fine grid in the model grid ! coslat = MAX( COS( lat_rel(ip,jp) * pi/180. ), 0.001 ) ax = (MIN(lon_up,louprel)-MAX(lon_low, lolowrel))*pi/180. * R_Earth * coslat ay = (MIN(lat_up, laup_rel(ip,jp))-MAX(lat_low,lalow_rel(ip,jp)))*pi/180. * R_Earth sub_area(ib, fopt) = ax*ay sub_index(ib, fopt, 1) = ip sub_index(ib, fopt, 2) = jp ENDIF ENDIF ! sub_pts(ib) = fopt ! ENDIF ! ENDDO ! ENDIF ! ENDDO ! ! ENDDO ! ! Do some memory management. ! DEALLOCATE (laup_rel) DEALLOCATE (loup_rel) DEALLOCATE (lalow_rel) DEALLOCATE (lolow_rel) DEALLOCATE (lat_ful) DEALLOCATE (lon_ful) ! nwbas = MAXVAL(sub_pts) ! ALLOCATE (basin_count(nbpt)) ALLOCATE (basin_area(nbpt,nwbas), basin_hierarchy(nbpt,nwbas), basin_topoind(nbpt,nwbas)) ALLOCATE (fetch_basin(nbpt,nwbas)) ALLOCATE (basin_id(nbpt,nwbas), basin_flowdir(nbpt,nwbas)) ALLOCATE (outflow_grid(nbpt,nwbas),outflow_basin(nbpt,nwbas)) ALLOCATE (inflow_number(nbpt,nwbas)) ALLOCATE (inflow_basin(nbpt,nwbas,nbvmax), inflow_grid(nbpt,nwbas,nbvmax)) ALLOCATE (nbcoastal(nbpt), coastal_basin(nbpt,nwbas)) ! ! Order all sub points in each grid_box and find the sub basins ! ! before we start we set the maps to empty ! basin_id(:,:) = undef_int basin_count(:) = 0 hierarchy(:,:) = undef_sechiba max_basins = MAXVAL(basins, MASK=basins .LT. 1.e10) invented_basins = max_basins nbcoastal(:) = 0 ! CALL routing_hierarchy(iml, jml, trip, topoindex, hierarchy) ! ! DO ib =1, nbpt ! ! ! Set everything to undef to locate easily empty points ! trip_bx(:,:) = undef_int basin_bx(:,:) = undef_int topoind_bx(:,:) = undef_sechiba area_bx(:,:) = undef_sechiba hierarchy_bx(:,:) = undef_sechiba ! ! extract the information for this grid box ! CALL routing_getgrid(nbpt, iml, jml, ib, sub_pts, sub_index, sub_area, max_basins, min_topoind, & & lon_rel, lat_rel, lalo, resolution, contfrac, trip, basins, topoindex, hierarchy, & & nbi, nbj, area_bx, trip_bx, basin_bx, topoind_bx, hierarchy_bx, lon_bx, lat_bx) ! CALL routing_findbasins(nbi, nbj, trip_bx, basin_bx, hierarchy_bx, topoind_bx,& & nb_basin, basin_inbxid, basin_sz, basin_bxout, basin_pts, coast_pts) ! ! Deal with the case where nb_basin=0 for this grid box. In this case all goes into coatal flow ! IF ( debug .AND. (COUNT(diagbox .EQ. ib) .GT. 0) ) THEN WRITE(numout,*) '===================== IB = :', ib WRITE(numout,*) "sub_pts(ib) :", sub_pts(ib) WRITE(numout,*) 'LON LAT of GCM :', lalo(ib,2), lalo(ib,1) WRITE(numout,*) 'Neighbor options :', neighbours(ib,1:8) WRITE(numout,*) 'Resolution :', resolution(ib,1:2) WRITE(fmt,"('(',I3,'I6)')") nbi WRITE(numout,*) '-------------> trip ', trip_bx(1,1) DO jp=1,nbj WRITE(numout,fmt) trip_bx(1:nbi,jp) ENDDO WRITE(numout,*) '-------------> basin ',basin_bx(1,1) DO jp=1,nbj WRITE(numout,fmt) basin_bx(1:nbi,jp) ENDDO WRITE(numout,*) '-------------> hierarchy ',hierarchy_bx(1,1) DO jp=1,nbj WRITE(numout,fmt) INT(hierarchy_bx(1:nbi,jp)/1000.) ENDDO WRITE(numout,*) '-------------> topoindex ',topoind_bx(1,1) DO jp=1,nbj WRITE(numout,fmt) INT(topoind_bx(1:nbi,jp)/1000.) ENDDO ! WRITE(numout,*) '------------> The basins we retain' DO jp=1,nb_basin WRITE(numout,*) 'index, size, bxout, coast :', basin_inbxid(jp), basin_sz(jp),& & basin_bxout(jp), coast_pts(jp) ENDDO ! ENDIF ! CALL routing_globalize(nbpt, ib, neighbours, area_bx, trip_bx, hierarchy_bx, topoind_bx, min_topoind,& & nb_basin, basin_inbxid, basin_sz, basin_pts, basin_bxout, coast_pts, nwbas, basin_count,& & basin_area, basin_hierarchy, basin_topoind, basin_id, basin_flowdir, outflow_grid,& & nbcoastal, coastal_basin) ! IF ( debug .AND. (COUNT(diagbox .EQ. ib) .GT. 0) ) THEN WRITE(numout,*) 'GLOBAL information after routing_globalize for grid ', ib DO jp=1,basin_count(ib) WRITE(numout,*) 'Basin ID : ', basin_id(ib, jp) WRITE(numout,*) 'Basin flowdir :', basin_flowdir(ib, jp) WRITE(numout,*) 'Basin hierarchy :', basin_hierarchy(ib, jp) WRITE(numout,*) 'Basin topoindex :', basin_topoind(ib, jp) WRITE(numout,*) 'Basin outflow grid :', outflow_grid(ib,jp) ENDDO ENDIF ! ENDDO ! CALL routing_linkup(nbpt, neighbours, nwbas, basin_count, basin_area, basin_id, basin_flowdir, & & basin_hierarchy, outflow_grid, outflow_basin, inflow_number, inflow_grid, inflow_basin, & & nbcoastal, coastal_basin, invented_basins) ! WRITE(numout,*) 'The maximum number of basins in any grid :', MAXVAL(basin_count) ! IF ( debug ) THEN DO ib=1,SIZE(diagbox) IF ( diagbox(ib) .GT. 0 ) THEN WRITE(numout,*) 'After routing_linkup information for grid ', diagbox(ib) DO jp=1,basin_count(diagbox(ib)) WRITE(numout,*) 'Basin ID : ', basin_id(diagbox(ib), jp) WRITE(numout,*) 'Basin outflow_grid :', outflow_grid(diagbox(ib), jp) WRITE(numout,*) 'Basin outflow_basin:', outflow_basin(diagbox(ib), jp) WRITE(numout,*) 'Basin hierarchy :', basin_hierarchy(diagbox(ib), jp) ENDDO ENDIF ENDDO ENDIF ! CALL routing_fetch(nbpt, resolution, contfrac, nwbas, basin_count, basin_area, basin_id, outflow_grid, & & outflow_basin, fetch_basin) ! WRITE(numout,*) "Start reducing the number of basins per grid to meet the required truncation." ! CALL routing_truncate(nbpt, resolution, contfrac, nwbas, basin_count, basin_area, basin_topoind,& & fetch_basin, basin_id, basin_flowdir, outflow_grid, outflow_basin, inflow_number,& & inflow_grid, inflow_basin) ! DEALLOCATE (lat_rel) DEALLOCATE (lon_rel) ! DEALLOCATE (trip) DEALLOCATE (basins) DEALLOCATE (topoindex) DEALLOCATE (hierarchy) ! DEALLOCATE (sub_area) DEALLOCATE (sub_index) DEALLOCATE (sub_pts) ! DEALLOCATE (basin_count) DEALLOCATE (basin_area, basin_hierarchy, basin_topoind, fetch_basin) DEALLOCATE (basin_id, basin_flowdir) DEALLOCATE (outflow_grid,outflow_basin) DEALLOCATE (inflow_number) DEALLOCATE (inflow_basin, inflow_grid) DEALLOCATE (nbcoastal, coastal_basin) ! RETURN ! END SUBROUTINE routing_basins ! !----------------------------------------------------------------------- ! SUBROUTINE routing_getgrid(nbpt, iml, jml, ib, sub_pts, sub_index, sub_area, max_basins, min_topoind, & & lon_rel, lat_rel, lalo, resolution, contfrac, trip, basins, topoindex, hierarchy, & & nbi, nbj, area_bx, trip_bx, basin_bx, topoind_bx, hierarchy_bx, lon_bx, lat_bx) ! IMPLICIT NONE ! ! Extracts from the global high resolution fileds the data for the current grid box ! we are dealing with. ! ! Convention for trip on the input : ! The trip field follows the following convention for the flow of the water : ! trip = 1 : flow = N ! trip = 2 : flow = NE ! trip = 3 : flow = E ! trip = 4 : flow = SE ! trip = 5 : flow = S ! trip = 6 : flow = SW ! trip = 7 : flow = W ! trip = 8 : flow = NW ! trip = 97 : return flow into the ground ! trip = 98 : coastal flow (diffuse flow into the oceans) These values are created here ! trip = 99 : river flow into the oceans ! ! On output, the grid boxes of the basin map which flow out of the GCM grid are identified ! by numbers larger than 100 : ! trip = 101 : flow = N out of the coarse grid ! trip = 102 : flow = NE out of the coarse grid ! trip = 103 : flow = E out of the coarse grid ! trip = 104 : flow = SE out of the coarse grid ! trip = 105 : flow = S out of the coarse grid ! trip = 106 : flow = SW out of the coarse grid ! trip = 107 : flow = W out of the coarse grid ! trip = 108 : flow = NW out of the coarse grid ! Inside the grid the convention remains the same as above (ie between 1 and 99). ! ! INPUT ! INTEGER(i_std), INTENT(in) :: nbpt ! Number of points in the global grid INTEGER(i_std), INTENT(in) :: iml, jml ! Resolution of the high resolution grid INTEGER(i_std), INTENT(in) :: ib ! point we are currently dealing with INTEGER(i_std), INTENT(in) :: sub_pts(nbpt) ! Number of high resiolution points on this grid INTEGER(i_std), INTENT(in) :: sub_index(nbpt, nbvmax,2) ! indeces of the points we need on the fine grid REAL(r_std), INTENT(inout) :: max_basins ! The current maximum of basins REAL(r_std), INTENT(in) :: min_topoind ! The current maximum of topographic index REAL(r_std), INTENT(in) :: sub_area(nbpt, nbvmax) ! area on the fine grid REAL(r_std), INTENT(in) :: lon_rel(iml, jml), lat_rel(iml, jml) ! coordinates of the fine grid REAL(r_std), INTENT(in) :: lalo(nbpt,2) ! Vector of latitude and longitudes (beware of the order !) REAL(r_std), INTENT(in) :: resolution(nbpt,2) ! The size in km of each grid-box in X and Y REAL(r_std), INTENT(in) :: contfrac(nbpt) ! Fraction of land in each grid box. REAL(r_std), INTENT(inout) :: trip(iml, jml), basins(iml, jml) ! data on the fine grid REAL(r_std), INTENT(inout) :: topoindex(iml, jml), hierarchy(iml, jml) ! data on the fine grid INTEGER(i_std), INTENT(out):: nbi, nbj ! Number of point ion x and y within the grid REAL(r_std), INTENT(out) :: area_bx(nbvmax,nbvmax), hierarchy_bx(nbvmax,nbvmax) REAL(r_std), INTENT(out) :: lon_bx(nbvmax,nbvmax), lat_bx(nbvmax,nbvmax), topoind_bx(nbvmax,nbvmax) INTEGER(i_std), INTENT(out):: trip_bx(nbvmax,nbvmax), basin_bx(nbvmax,nbvmax) ! ! LOCAL ! INTEGER(i_std) :: ip, jp, ll(1), iloc, jloc REAL(r_std) :: lonstr(nbvmax*nbvmax), latstr(nbvmax*nbvmax) ! IF ( sub_pts(ib) > 0 ) THEN ! DO ip=1,sub_pts(ib) lonstr(ip) = lon_rel(sub_index(ib, ip, 1), sub_index(ib, ip, 2)) latstr(ip) = lat_rel(sub_index(ib, ip, 1), sub_index(ib, ip, 2)) ENDDO ! ! Get the size of the area and order the coordinates to go from North to South and West to East ! CALL routing_sortcoord(sub_pts(ib), lonstr, 'WE', nbi) CALL routing_sortcoord(sub_pts(ib), latstr, 'NS', nbj) ! ! Transfer the data in such a way that (1,1) is the North Western corner and ! (nbi, nbj) the South Eastern. ! DO ip=1,sub_pts(ib) ll = MINLOC(ABS(lonstr(1:nbi) - lon_rel(sub_index(ib, ip, 1), sub_index(ib, ip, 2)))) iloc = ll(1) ll = MINLOC(ABS(latstr(1:nbj) - lat_rel(sub_index(ib, ip, 1), sub_index(ib, ip, 2)))) jloc = ll(1) trip_bx(iloc, jloc) = NINT(trip(sub_index(ib, ip, 1), sub_index(ib, ip, 2))) basin_bx(iloc, jloc) = NINT(basins(sub_index(ib, ip, 1), sub_index(ib, ip, 2))) area_bx(iloc, jloc) = sub_area(ib, ip) topoind_bx(iloc, jloc) = topoindex(sub_index(ib, ip, 1), sub_index(ib, ip, 2)) hierarchy_bx(iloc, jloc) = hierarchy(sub_index(ib, ip, 1), sub_index(ib, ip, 2)) lon_bx(iloc, jloc) = lon_rel(sub_index(ib, ip, 1), sub_index(ib, ip, 2)) lat_bx(iloc, jloc) = lat_rel(sub_index(ib, ip, 1), sub_index(ib, ip, 2)) ENDDO ELSE ! ! This is the case where the model invented a continental point ! nbi = 1 nbj = 1 iloc = 1 jloc = 1 trip_bx(iloc, jloc) = 98 basin_bx(iloc, jloc) = NINT(max_basins + 1) max_basins = max_basins + 1 area_bx(iloc, jloc) = resolution(ib,1)*resolution(ib,2)*contfrac(ib) topoind_bx(iloc, jloc) = min_topoind hierarchy_bx(iloc, jloc) = min_topoind lon_bx(iloc, jloc) = lalo(ib,2) lat_bx(iloc, jloc) = lalo(ib,1) ! ENDIF ! ! Tag in trip all the outflow conditions. The table is thus : ! trip = 100+n : Outflow into another grid-box ! trip = 99 : River outflow into the ocean ! trip = 98 : This will be coastal flow (not organized as a basin) ! trip = 97 : return flow into the soil (local) ! DO jp=1,nbj IF ( trip_bx(1,jp) .EQ. 8 .OR. trip_bx(1,jp) .EQ. 7 .OR. trip_bx(1,jp) .EQ. 6) THEN trip_bx(1,jp) = trip_bx(1,jp) + 100 ENDIF IF ( trip_bx(nbi,jp) .EQ. 2 .OR. trip_bx(nbi,jp) .EQ. 3 .OR. trip_bx(nbi,jp) .EQ. 4) THEN trip_bx(nbi,jp) = trip_bx(nbi,jp) + 100 ENDIF ENDDO DO ip=1,nbi IF ( trip_bx(ip,1) .EQ. 8 .OR. trip_bx(ip,1) .EQ. 1 .OR. trip_bx(ip,1) .EQ. 2) THEN trip_bx(ip,1) = trip_bx(ip,1) + 100 ENDIF IF ( trip_bx(ip,nbj) .EQ. 6 .OR. trip_bx(ip,nbj) .EQ. 5 .OR. trip_bx(ip,nbj) .EQ. 4) THEN trip_bx(ip,nbj) = trip_bx(ip,nbj) + 100 ENDIF ENDDO ! ! ! We simplify the outflow. We only need the direction normal to the ! box boundary and the 4 corners. ! ! Northern border IF ( trip_bx(1,1) .EQ. 102 ) trip_bx(1,1) = 101 IF ( trip_bx(nbi,1) .EQ. 108 ) trip_bx(nbi,1) = 101 DO ip=2,nbi-1 IF ( trip_bx(ip,1) .EQ. 108 .OR. trip_bx(ip,1) .EQ. 102 ) trip_bx(ip,1) = 101 ENDDO ! Southern border IF ( trip_bx(1,nbj) .EQ. 104 ) trip_bx(1,nbj) = 105 IF ( trip_bx(nbi,nbj) .EQ. 106 ) trip_bx(nbi,nbj) = 105 DO ip=2,nbi-1 IF ( trip_bx(ip,nbj) .EQ. 104 .OR. trip_bx(ip,nbj) .EQ. 106 ) trip_bx(ip,nbj) = 105 ENDDO ! Eastern border IF ( trip_bx(nbi,1) .EQ. 104) trip_bx(nbi,1) = 103 IF ( trip_bx(nbi,nbj) .EQ. 102) trip_bx(nbi,nbj) = 103 DO jp=2,nbj-1 IF ( trip_bx(nbi,jp) .EQ. 104 .OR. trip_bx(nbi,jp) .EQ. 102 ) trip_bx(nbi,jp) = 103 ENDDO ! Western border IF ( trip_bx(1,1) .EQ. 106) trip_bx(1,1) = 107 IF ( trip_bx(1,nbj) .EQ. 108) trip_bx(1,nbj) = 107 DO jp=2,nbj-1 IF ( trip_bx(1,jp) .EQ. 106 .OR. trip_bx(1,jp) .EQ. 108 ) trip_bx(1,jp) = 107 ENDDO ! ! END SUBROUTINE routing_getgrid ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! SUBROUTINE routing_sortcoord(nb_in, coords, direction, nb_out) ! IMPLICIT NONE ! ! Input/Output ! INTEGER(i_std), INTENT(in) :: nb_in REAL(r_std), INTENT(inout) :: coords(nb_in) CHARACTER(LEN=2) :: direction INTEGER(i_std), INTENT(out) :: nb_out ! ! Local ! INTEGER(i_std) :: ipos REAL(r_std) :: coords_tmp(nb_in) INTEGER(i_std), DIMENSION(1) :: ll INTEGER(i_std) :: ind(nb_in) ! ipos = 1 nb_out = nb_in ! ! Compress the coordinates array ! DO WHILE ( ipos < nb_in ) IF ( coords(ipos+1) /= undef_sechiba) THEN IF ( COUNT(coords(ipos:nb_out) == coords(ipos)) > 1 ) THEN coords(ipos:nb_out-1) = coords(ipos+1:nb_out) coords(nb_out:nb_in) = undef_sechiba nb_out = nb_out - 1 ELSE ipos = ipos + 1 ENDIF ELSE EXIT ENDIF ENDDO ! ! Sort it now ! ! First we get ready and adjust for the periodicity in longitude ! coords_tmp(:) = undef_sechiba IF ( INDEX(direction, 'WE') == 1 .OR. INDEX(direction, 'EW') == 1) THEN IF ( MAXVAL(ABS(coords(1:nb_out))) .GT. 160 ) THEN coords_tmp(1:nb_out) = MOD(coords(1:nb_out) + 360.0, 360.0) ELSE coords_tmp(1:nb_out) = coords(1:nb_out) ENDIF ELSE IF ( INDEX(direction, 'NS') == 1 .OR. INDEX(direction, 'SN') == 1) THEN coords_tmp(1:nb_out) = coords(1:nb_out) ELSE WRITE(numout,*) 'The chosen direction (', direction,') is not recognized' STOP 'routing_sortcoord' ENDIF ! ! Get it sorted out now ! ipos = 1 ! IF ( INDEX(direction, 'WE') == 1 .OR. INDEX(direction, 'SN') == 1) THEN DO WHILE (COUNT(ABS(coords_tmp(:)-undef_sechiba) > EPSILON(undef_sechiba)*10.) >= 1) ll = MINLOC(coords_tmp(:), coords_tmp /= undef_sechiba) ind(ipos) = ll(1) coords_tmp(ll(1)) = undef_sechiba ipos = ipos + 1 ENDDO ELSE IF ( INDEX(direction, 'EW') == 1 .OR. INDEX(direction, 'NS') == 1) THEN DO WHILE (COUNT(ABS(coords_tmp(:)-undef_sechiba) > EPSILON(undef_sechiba)*10.) >= 1) ll = MAXLOC(coords_tmp(:), coords_tmp /= undef_sechiba) ind(ipos) = ll(1) coords_tmp(ll(1)) = undef_sechiba ipos = ipos + 1 ENDDO ELSE WRITE(numout,*) 'The chosen direction (', direction,') is not recognized (second)' STOP 'routing_sortcoord' ENDIF ! coords(1:nb_out) = coords(ind(1:nb_out)) IF (nb_out < nb_in) THEN coords(nb_out+1:nb_in) = zero ENDIF ! END SUBROUTINE routing_sortcoord ! !------------------------------------------------------------------------------------------------- ! SUBROUTINE routing_findbasins(nbi, nbj, trip, basin, hierarchy, topoind, nb_basin, basin_inbxid, basin_sz,& & basin_bxout, basin_pts, coast_pts) ! IMPLICIT NONE ! ! This subroutine find the basins and does some clean up. The aim is to return the list off all ! points which are within the same basin of the grid_box. ! We will also collect all points which directly flow into the ocean in one basin ! Make sure that we do not have a basin with two outflows and other exceptions. ! ! At this stage no effort is made to come down to the truncation of the model. ! ! Convention for trip ! ------------------- ! Inside of the box : ! trip = 1 : flow = N ! trip = 2 : flow = NE ! trip = 3 : flow = E ! trip = 4 : flow = SE ! trip = 5 : flow = S ! trip = 6 : flow = SW ! trip = 7 : flow = W ! trip = 8 : flow = NW ! trip = 97 : return flow into the ground ! trip = 98 : coastal flow (diffuse flow into the oceans) These values are created here ! trip = 99 : river flow into the oceans ! ! Out flow from the gird : ! trip = 101 : flow = N out of the coarse grid ! trip = 102 : flow = NE out of the coarse grid ! trip = 103 : flow = E out of the coarse grid ! trip = 104 : flow = SE out of the coarse grid ! trip = 105 : flow = S out of the coarse grid ! trip = 106 : flow = SW out of the coarse grid ! trip = 107 : flow = W out of the coarse grid ! trip = 108 : flow = NW out of the coarse grid ! ! Inputs ! INTEGER(i_std) :: nbi, nbj INTEGER(i_std) :: trip(:,:), basin(:,:) REAL(r_std) :: hierarchy(:,:), topoind(:,:) ! ! Outputs ! INTEGER(i_std) :: nb_basin INTEGER(i_std) :: basin_inbxid(nbvmax), basin_sz(nbvmax), basin_bxout(nbvmax) INTEGER(i_std) :: basin_pts(nbvmax, nbvmax, 2) INTEGER(i_std) :: coast_pts(nbvmax) ! ! Local ! INTEGER(i_std) :: ibas, ilf, nbb, nb_in INTEGER(i_std) :: bname(nbvmax), sz(nbvmax), pts(nbvmax,nbvmax,2), nbout(nbvmax) INTEGER(i_std) :: new_nb, new_bname(nbvmax), new_sz(nbvmax), new_pts(nbvmax,nbvmax,2) INTEGER(i_std) :: itrans, trans(nbvmax), outdir(nbvmax) INTEGER(i_std) :: tmpsz(nbvmax) INTEGER(i_std) :: ip, jp, jpp(1), ipb INTEGER(i_std) :: sortind(nbvmax) CHARACTER(LEN=7) :: fmt ! nbb = 0 bname(:) = undef_int sz(:) = 0 nbout(:) = 0 new_pts(:,:,:) = 0 ! ! 1.0 Find all basins within this grid-box ! Sort the variables per basin so that we can more easily ! access data from the same basin (The variables are : ! bname, sz, pts, nbout) ! DO ip=1,nbi DO jp=1,nbj IF ( basin(ip,jp) .LT. undef_int) THEN IF ( COUNT(basin(ip,jp) .EQ. bname(:)) .EQ. 0 ) THEN nbb = nbb + 1 IF ( nbb .GT. nbvmax ) STOP 'nbvmax too small' bname(nbb) = basin(ip,jp) sz(nbb) = 0 ENDIF ! DO ilf=1,nbb IF ( basin(ip,jp) .EQ. bname(ilf) ) THEN ibas = ilf ENDIF ENDDO ! sz(ibas) = sz(ibas) + 1 IF ( sz(ibas) .GT. nbvmax ) STOP 'nbvmax too small' pts(ibas, sz(ibas), 1) = ip pts(ibas, sz(ibas), 2) = jp ! We deal only with outflow and leave flow back into the grid-box for later. IF ( trip(ip,jp) .GE. 97 ) THEN nbout(ibas) = nbout(ibas) + 1 ENDIF ! ENDIF ! ENDDO ENDDO ! ! 2.0 All basins which have size 1 and flow to the ocean are put together. ! itrans = 0 coast_pts(:) = undef_int ! Get all the points we can collect DO ip=1,nbb IF ( sz(ip) .EQ. 1 .AND. trip(pts(ip,1,1),pts(ip,1,2)) .EQ. 99) THEN itrans = itrans + 1 trans(itrans) = ip trip(pts(ip,1,1),pts(ip,1,2)) = 98 ENDIF ENDDO ! put everything in the first basin IF ( itrans .GT. 1) THEN ipb = trans(1) coast_pts(sz(ipb)) = bname(ipb) bname(ipb) = -1 DO ip=2,itrans sz(ipb) = sz(ipb) + 1 coast_pts(sz(ipb)) = bname(trans(ip)) sz(trans(ip)) = 0 pts(ipb, sz(ipb), 1) = pts(trans(ip), 1, 1) pts(ipb, sz(ipb), 2) = pts(trans(ip), 1, 2) ENDDO ENDIF ! ! 3.0 Make sure that we have only one outflow point in each basin ! ! nbb is the number of basins on this grid box. new_nb = 0 DO ip=1,nbb ! We only do this for grid-points which have more than one outflow IF ( sz(ip) .GT. 1 .AND. nbout(ip) .GT. 1) THEN ! ! Pick up all points needed and store them in trans ! itrans = 0 DO jp=1,sz(ip) IF ( trip(pts(ip,jp,1),pts(ip,jp,2)) .GE. 97) THEN itrans = itrans + 1 trans(itrans) = trip(pts(ip,jp,1),pts(ip,jp,2)) ENDIF ENDDO ! ! First issue : We have more than one point of the basin which flows into ! the ocean. In this case we put everything into coastal flow. It will go into ! a separate basin in the routing_globalize routine. ! IF ( (COUNT(trans(1:itrans) .EQ. 99) + COUNT(trans(1:itrans) .EQ. 98)) .GT. 1) THEN DO jp=1,sz(ip) IF ( trip(pts(ip,jp,1),pts(ip,jp,2)) .EQ. 99 ) THEN trip(pts(ip,jp,1),pts(ip,jp,2)) = 98 trans(itrans) = trip(pts(ip,jp,1),pts(ip,jp,2)) ENDIF ENDDO ENDIF ! ! Second issue : We have redundant outflows at the boundaries. That is two small grid ! boxes flowing into the same GCM grid box. ! IF ( COUNT(trans(1:itrans) .GT. 100) .GE. 1) THEN CALL routing_simplify(nbi, nbj, trip, basin, hierarchy, bname(ip)) itrans = 0 DO jp=1,sz(ip) IF ( trip(pts(ip,jp,1),pts(ip,jp,2)) .GE. 9) THEN itrans = itrans + 1 trans(itrans) = trip(pts(ip,jp,1),pts(ip,jp,2)) ENDIF ENDDO ENDIF ! ! Third issue : we have more than one outflow from the boxes. This could be ! - flow into 2 or more neighboring GCM grids ! - flow into a neighboring GCM grids and into the ocean or be a return flow (=97. =98, =99) ! - flow into a neighboring GCM grids or ocean and back into the same GCM grid box ! The only solution is to cut the basin up in as many parts. ! IF ( COUNT(trans(1:itrans) .GE. 97) .GT. 1) THEN ! nb_in = new_nb CALL routing_cutbasin(nbi, nbj, nbb, trip, basin, bname(ip), new_nb, new_bname, new_sz, new_pts) ! ! If we have split the basin then we need to cancel the old one ! IF ( nb_in .NE. new_nb) THEN sz(ip) = 0 ENDIF ! ENDIF ! ENDIF ENDDO ! ! Add the new basins to the end of the list ! If ( nbb+new_nb .LE. nbvmax) THEN DO ip=1,new_nb bname(nbb+ip) = new_bname(ip) sz(nbb+ip) = new_sz(ip) pts(nbb+ip,:,:) = new_pts(ip,:,:) ENDDO nbb = nbb+new_nb ELSE WRITE(numout,*) 'Increase nbvmax. It is too small to contain all the basins (routing_findbasins)' STOP ENDIF ! ! Keep the output direction ! DO ip=1,nbb IF ( sz(ip) .GT. 0 ) THEN trans(:) = 0 DO jp=1,sz(ip) trans(jp) = trip(pts(ip,jp,1),pts(ip,jp,2)) ENDDO outdir(ip) = MAXVAL(trans(1:sz(ip))) IF ( outdir(ip) .GE. 97 ) THEN outdir(ip) = outdir(ip) - 100 ELSE WRITE(numout,*) 'Why are we here and can not find a trip larger than 96' WRITE(numout,*) 'Does this mean that the basin does not have any outflow ', ip, bname(ip) WRITE(fmt,"('(',I3,'I9)')") nbi WRITE(numout,*) '-----------------------> trip' DO jp=1,nbj WRITE(numout,fmt) trip(1:nbi,jp) ENDDO WRITE(numout,*) '-----------------------> basin' DO jp=1,nbj WRITE(numout,fmt) basin(1:nbi,jp) ENDDO STOP 'SUBROUTINE : routing_findbasins' ENDIF ENDIF ENDDO ! ! ! Sort the output by size of the various basins. ! nb_basin = COUNT(sz(1:nbb) .GT. 0) tmpsz(:) = -1 tmpsz(1:nbb) = sz(1:nbb) DO ip=1,nbb jpp = MAXLOC(tmpsz(:)) IF ( sz(jpp(1)) .GT. 0) THEN sortind(ip) = jpp(1) tmpsz(jpp(1)) = -1 ENDIF ENDDO basin_inbxid(1:nb_basin) = bname(sortind(1:nb_basin)) basin_sz(1:nb_basin) = sz(sortind(1:nb_basin)) basin_pts(1:nb_basin,:,:) = pts(sortind(1:nb_basin),:,:) basin_bxout(1:nb_basin) = outdir(sortind(1:nb_basin)) ! ! We can only check if we have at least as many outflows as basins ! ip = COUNT(trip(1:nbi,1:nbj) .GE. 97 .AND. trip(1:nbi,1:nbj) .LT. undef_int) !! ip = ip + COUNT(trip(1:nbi,1:nbj) .EQ. 97) !! IF ( COUNT(trip(1:nbi,1:nbj) .EQ. 98) .GT. 0) ip = ip + 1 IF ( ip .LT. nb_basin ) THEN WRITE(numout,*) 'We have less outflow points than basins :', ip WRITE(fmt,"('(',I3,'I9)')") nbi WRITE(numout,*) '-----------------------> trip' DO jp=1,nbj WRITE(numout,fmt) trip(1:nbi,jp) ENDDO WRITE(numout,*) '-----------------------> basin' DO jp=1,nbj WRITE(numout,fmt) basin(1:nbi,jp) ENDDO WRITE(numout,*) 'nb_basin :', nb_basin WRITE(numout,*) 'Basin sized :', basin_sz(1:nb_basin) STOP 'in routing_findbasins' ENDIF ! END SUBROUTINE routing_findbasins ! ! ------------------------------------------------------------------------------------------ ! SUBROUTINE routing_simplify(nbi, nbj, trip, basin, hierarchy, basin_inbxid) ! IMPLICIT NONE ! ! This subroutine symplifies the routing out of each basin by taking ! out redundancies at the borders of the GCM box. The aim is to have ! only one outflow point per basin. But here we will not change the ! the direction of the outflow. ! ! Inputs INTEGER(i_std) :: nbi, nbj INTEGER(i_std) :: trip(:,:), basin(:,:) REAL(r_std) :: hierarchy(:,:) INTEGER(i_std) :: basin_inbxid ! ! Local ! INTEGER(i_std) :: ip, jp, nbout, basin_sz, iborder INTEGER(i_std), DIMENSION(nbvmax,nbvmax) :: trip_tmp INTEGER(i_std), DIMENSION(nbvmax,nbvmax,2) :: trip_flow INTEGER(i_std), DIMENSION(nbvmax,2) :: outflow INTEGER(i_std), DIMENSION(nbvmax) :: outsz CHARACTER(LEN=7) :: fmt ! INTEGER(i_std), DIMENSION(8,2) :: inc INTEGER(i_std) :: itodo, ill(1), icc, ismall, ibas, iip, jjp, ib, id INTEGER(i_std), DIMENSION(nbvmax) :: todopt, todosz REAL(r_std), DIMENSION(nbvmax) :: todohi LOGICAL :: not_found, debug = .FALSE. ! ! ! The routing code (i=1, j=2) ! inc(1,1) = 0 inc(1,2) = -1 inc(2,1) = 1 inc(2,2) = -1 inc(3,1) = 1 inc(3,2) = 0 inc(4,1) = 1 inc(4,2) = 1 inc(5,1) = 0 inc(5,2) = 1 inc(6,1) = -1 inc(6,2) = 1 inc(7,1) = -1 inc(7,2) = 0 inc(8,1) = -1 inc(8,2) = -1 ! ! ! Symplify the outflow conditions first. We are only interested in the ! outflows which go to different GCM grid-boxes. ! IF ( debug ) THEN WRITE(numout,*) '+++++++++++++++++++ BEFORE ANYTHING ++++++++++++++++++++' WRITE(fmt,"('(',I3,'I6)')") nbi DO jp=1,nbj WRITE(numout,fmt) trip_tmp(1:nbi,jp) ENDDO ENDIF ! ! transfer the trips into an array which only contains the basin we are interested in ! trip_tmp(:,:) = -1 basin_sz = 0 DO ip=1,nbi DO jp=1,nbj IF ( basin(ip,jp) .EQ. basin_inbxid) THEN trip_tmp(ip,jp) = trip(ip,jp) basin_sz = basin_sz + 1 ENDIF ENDDO ENDDO ! ! Determine for each point where it flows to ! CALL routing_findrout(nbi, nbj, trip_tmp, basin_sz, basin_inbxid, nbout, outflow, trip_flow, outsz) ! ! ! ! ! Over the width of a GCM grid box we can have many outflows but we are interested ! in only one for each basin. Thus we wish to collect them all to form only one outflow ! to the neighboring grid-box. ! DO iborder = 101,107,2 ! ! If we have more than one of these outflows then we need to merge the sub-basins ! icc = COUNT(trip_tmp .EQ. iborder)-1 DO WHILE ( icc .GT. 0) ! Pick out all the points we will have to do itodo = 0 DO ip=1,nbout IF (trip_tmp(outflow(ip,1),outflow(ip,2)) .EQ. iborder) THEN itodo = itodo + 1 todopt(itodo) = ip todosz(itodo) = outsz(ip) ! We take the hierarchy of the outflow point as we will try to ! minimize if for the outflow of the entire basin. todohi(itodo) = hierarchy(outflow(ip,1),outflow(ip,2)) ENDIF ENDDO ! ! We change the direction of the smalest basin. ! ill=MAXLOC(todohi(1:itodo)) ismall = todopt(ill(1)) ! DO ip=1,nbi DO jp=1,nbj IF ( trip_flow(ip,jp,1) .EQ. outflow(ismall,1) .AND.& & trip_flow(ip,jp,2) .EQ. outflow(ismall,2) ) THEN ! Now that we have found a point of the smallest sub-basin we ! look around for another sub-basin ib = 1 not_found = .TRUE. DO WHILE ( not_found .AND. ib .LE. itodo ) IF ( ib .NE. ill(1) ) THEN ibas = todopt(ib) DO id=1,8 iip = ip + inc(id,1) jjp = jp + inc(id,2) ! Can we look at this points or is there any need to ? IF ( iip .GE. 1 .AND. iip .LE. nbi .AND. & & jjp .GE. 1 .AND. jjp .LE. nbj .AND. not_found) THEN ! Is this point the one we look for ? IF ( trip_flow(iip,jjp,1) .EQ. outflow(ibas,1) .AND. & & trip_flow(iip,jjp,2) .EQ. outflow(ibas,2)) THEN trip_flow(ip,jp,1) = outflow(ibas,1) trip_flow(ip,jp,2) = outflow(ibas,2) trip_tmp(ip,jp) = id ! This last line ensures that we do not come back to this point ! and that in the end the outer while will stop not_found = .FALSE. ENDIF ENDIF ENDDO ENDIF ib = ib + 1 ENDDO ENDIF ENDDO ENDDO ! icc = icc - 1 ENDDO ! ! ENDDO ! IF ( debug ) THEN WRITE(numout,*) '+++++++++++++++++++ AFTER +++++++++++++++++++++++++++++' WRITE(fmt,"('(',I3,'I6)')") nbi DO jp=1,nbj WRITE(numout,fmt) trip_tmp(1:nbi,jp) ENDDO ENDIF ! ! Put trip_tmp back into trip ! DO ip=1,nbi DO jp=1,nbj IF ( trip_tmp(ip,jp) .GT. 0) THEN trip(ip,jp) = trip_tmp(ip,jp) ENDIF ENDDO ENDDO ! END SUBROUTINE routing_simplify ! !------------------------------------- ! SUBROUTINE routing_cutbasin (nbi, nbj, nbbasins, trip, basin, basin_inbxid, nb, bname, sz, pts) ! IMPLICIT NONE ! ! This subroutine cuts the original basin which has more than one outflow into as ! many subbasins as outflow directions. ! ! Inputs INTEGER(i_std) :: nbi, nbj INTEGER(i_std) :: nbbasins INTEGER(i_std) :: trip(:,:), basin(:,:) INTEGER(i_std) :: basin_inbxid ! ! Outputs ! INTEGER(i_std) :: nb, bname(nbvmax), sz(nbvmax), pts(nbvmax,nbvmax,2) ! ! Local ! INTEGER(i_std) :: ip, jp, iip, jjp, ib, ibb, id, nbout INTEGER(i_std) :: basin_sz, nb_in INTEGER(i_std), DIMENSION(nbvmax,nbvmax) :: trip_tmp INTEGER(i_std), DIMENSION(nbvmax,nbvmax,2) :: trip_flow INTEGER(i_std), DIMENSION(nbvmax,2) :: outflow INTEGER(i_std), DIMENSION(nbvmax) :: outsz CHARACTER(LEN=7) :: fmt LOGICAL :: not_found, debug=.FALSE. ! INTEGER(i_std), DIMENSION(8,2) :: inc ! ! ! The routing code (i=1, j=2) ! inc(1,1) = 0 inc(1,2) = -1 inc(2,1) = 1 inc(2,2) = -1 inc(3,1) = 1 inc(3,2) = 0 inc(4,1) = 1 inc(4,2) = 1 inc(5,1) = 0 inc(5,2) = 1 inc(6,1) = -1 inc(6,2) = 1 inc(7,1) = -1 inc(7,2) = 0 inc(8,1) = -1 inc(8,2) = -1 ! ! Set up a temporary trip field which only contains the values ! for the basin on which we currently work. ! trip_tmp(:,:) = -1 basin_sz = 0 DO ip=1,nbi DO jp=1,nbj IF ( basin(ip,jp) .EQ. basin_inbxid) THEN trip_tmp(ip,jp) = trip(ip,jp) basin_sz = basin_sz + 1 ENDIF ENDDO ENDDO ! CALL routing_findrout(nbi, nbj, trip_tmp, basin_sz, basin_inbxid, nbout, outflow, trip_flow, outsz) ! !! IF ( debug ) THEN !! DO ib = nb_in+1,nb !! DO ip=1,sz(ib) !! trip_tmp(pts(ib, ip, 1),pts(ib, ip, 2)) = ib*(-1)-900 !! ENDDO !! ENDDO !! WRITE(fmt,"('(',I3,'I6)')") nbi !! WRITE(numout,*) 'BEFORE ------------> New basins ' !! WRITE(numout,*) nb, ' sz :', sz(1:nb) !! DO jp=1,nbj !! WRITE(numout,fmt) trip_tmp(1:nbi,jp) !! ENDDO !! ENDIF ! ! Take out the small sub-basins. That is those which have only one grid box ! This is only done if we need to save space in the number of basins. Else we ! can take it easy and keep diverging sub-basins for the moment. ! IF ( nbbasins .GE. nbasmax ) THEN DO ib=1,nbout ! If the sub-basin is of size one and its larger neighbor is flowing into another ! direction then we put them together. IF ( outsz(ib) .EQ. 1 .AND. trip(outflow(ib,1), outflow(ib,2)) .GT. 99 ) THEN ! not_found = .TRUE. DO id=1,8 ip = outflow(ib,1) jp = outflow(ib,2) iip = ip + inc(id,1) jjp = jp + inc(id,2) ! Can we look at this points ? IF ( iip .GE. 1 .AND. iip .LE. nbi .AND. & & jjp .GE. 1 .AND. jjp .LE. nbj .AND. not_found) THEN ! Did we find a direct neighbor which is an outflow point ? IF ( trip_tmp(iip,jjp) .GT. 100 ) THEN ! IF so direct the flow towards it and update the tables. not_found = .FALSE. trip(ip,jp) = id trip_tmp(ip,jp) = id outsz(ib) = 0 ! update the table of this basin DO ibb=1,nbout IF ( iip .EQ. outflow(ibb,1) .AND. jjp .EQ. outflow(ibb,2) ) THEN outsz(ibb) = outsz(ibb)+1 trip_flow(ip,jp,1) = outflow(ibb,1) trip_flow(ip,jp,2) = outflow(ibb,2) ENDIF ENDDO ENDIF ENDIF ENDDO ENDIF ENDDO ENDIF ! ! ! Cut the basin if we have more than 1 left. ! ! IF ( COUNT(outsz(1:nbout) .GT. 0) .GT. 1 ) THEN ! nb_in = nb ! DO ib = 1,nbout IF ( outsz(ib) .GT. 0) THEN nb = nb+1 IF ( nb .GT. nbvmax) THEN WRITE(numout,*) 'nbvmax too small, increase it (routing_cutbasin)' ENDIF bname(nb) = basin_inbxid sz(nb) = 0 DO ip=1,nbi DO jp=1,nbj IF ( (trip_flow(ip,jp,1) + trip_flow(ip,jp,1)) .GT. 0 .AND. & & trip_flow(ip,jp,1) .EQ. outflow(ib,1) .AND. & & trip_flow(ip,jp,2) .EQ. outflow(ib,2) ) THEN sz(nb) = sz(nb) + 1 pts(nb, sz(nb), 1) = ip pts(nb, sz(nb), 2) = jp ENDIF ENDDO ENDDO ENDIF ENDDO ! A short verification IF ( SUM(sz(nb_in+1:nb)) .NE. basin_sz) THEN WRITE(numout,*) 'Lost some points while spliting the basin' WRITE(numout,*) 'nbout :', nbout DO ib = nb_in+1,nb WRITE(numout,*) 'ib, SZ :', ib, sz(ib) ENDDO WRITE(fmt,"('(',I3,'I6)')") nbi WRITE(numout,*) '-------------> trip ' DO jp=1,nbj WRITE(numout,fmt) trip_tmp(1:nbi,jp) ENDDO STOP ENDIF ! IF ( debug ) THEN DO ib = nb_in+1,nb DO ip=1,sz(ib) trip_tmp(pts(ib, ip, 1),pts(ib, ip, 2)) = ib*(-1)-900 ENDDO ENDDO WRITE(fmt,"('(',I3,'I6)')") nbi WRITE(numout,*) 'AFTER-------------> New basins ' WRITE(numout,*) nb, ' sz :', sz(1:nb) DO jp=1,nbj WRITE(numout,fmt) trip_tmp(1:nbi,jp) ENDDO IF ( MAXVAl(trip_tmp(1:nbi,1:nbj)) .GT. 0) THEN STOP ENDIF ENDIF ENDIF ! END SUBROUTINE routing_cutbasin ! !------------------------------------------------------------------------ ! SUBROUTINE routing_hierarchy(iml, jml, trip, topoindex, hierarchy) ! IMPLICIT NONE ! ! This subroutine will find for each point the distance to the outflow point ! along the flowlines of the basin. ! INTEGER(i_std) :: iml, jml REAL(r_std), DIMENSION(iml,jml) :: trip, hierarchy, topoindex ! INTEGER(i_std), DIMENSION(8,2) :: inc INTEGER(i_std) :: ip, jp, ib, ntripi, ntripj, cnt, trp REAL(r_std) :: topohier, topohier_old CHARACTER(LEN=7) :: fmt ! ! The routing code (i=1, j=2) ! inc(1,1) = 0 inc(1,2) = -1 inc(2,1) = 1 inc(2,2) = -1 inc(3,1) = 1 inc(3,2) = 0 inc(4,1) = 1 inc(4,2) = 1 inc(5,1) = 0 inc(5,2) = 1 inc(6,1) = -1 inc(6,2) = 1 inc(7,1) = -1 inc(7,2) = 0 inc(8,1) = -1 inc(8,2) = -1 ! DO ip=1,iml DO jp=1,jml IF ( trip(ip,jp) .LT. undef_sechiba ) THEN ntripi = ip ntripj = jp trp = trip(ip,jp) cnt = 1 ! Warn for extreme numbers IF ( topoindex(ip,jp) .GT. 1.e10 ) THEN WRITE(numout,*) 'We have a very large topographic index for point ', ip, jp WRITE(numout,*) 'This can not be right :', topoindex(ip,jp) STOP ELSE topohier = topoindex(ip,jp) ENDIF ! DO WHILE ( trp .GT. 0 .AND. trp .LT. 9 .AND. cnt .LT. iml*jml) cnt = cnt + 1 ntripi = ntripi + inc(trp,1) IF ( ntripi .LT. 1) ntripi = iml IF ( ntripi .GT. iml) ntripi = 1 ntripj = ntripj + inc(trp,2) topohier_old = topohier topohier = topohier + topoindex(ntripi, ntripj) IF ( topohier_old .GT. topohier) THEN WRITE(numout,*) 'Big Problem, how comes we climb up a hill ?' WRITE(numout,*) 'The old value of topographicaly weighted hierarchy was : ', topohier_old WRITE(numout,*) 'The new one is :', topohier STOP 'routing_hierarchy' ENDIF trp = trip(ntripi, ntripj) ENDDO ! IF ( cnt .EQ. iml*jml) THEN WRITE(numout,*) 'We could not route point (routing_findrout) :', ip, jp WRITE(numout,*) '-------------> trip ' WRITE(fmt,"('(',I3,'I6)')") iml DO ib=1,jml WRITE(numout,fmt) trip(1:iml,ib) ENDDO STOP ENDIF ! hierarchy(ip,jp) = topohier ! ENDIF ENDDO ENDDO ! ! END SUBROUTINE routing_hierarchy ! !------------------------------------------------------------------------ ! SUBROUTINE routing_findrout(nbi, nbj, trip, basin_sz, basinid, nbout, outflow, trip_flow, outsz) ! IMPLICIT NONE ! ! This subroutine simply computes the route to each outflow point and returns ! the outflow point for each point in the basin. ! ! INPUT ! INTEGER(i_std) :: nbi, nbj INTEGER(i_std), DIMENSION(nbvmax,nbvmax) :: trip INTEGER(i_std) :: basin_sz, basinid, nbout INTEGER(i_std), DIMENSION(nbvmax,nbvmax,2) :: trip_flow INTEGER(i_std), DIMENSION(nbvmax,2) :: outflow INTEGER(i_std), DIMENSION(nbvmax) :: outsz ! ! LOCAL ! INTEGER(i_std), DIMENSION(8,2) :: inc INTEGER(i_std) :: ip, jp, ib, cnt, trp, totsz CHARACTER(LEN=7) :: fmt ! ! ! The routing code (i=1, j=2) ! inc(1,1) = 0 inc(1,2) = -1 inc(2,1) = 1 inc(2,2) = -1 inc(3,1) = 1 inc(3,2) = 0 inc(4,1) = 1 inc(4,2) = 1 inc(5,1) = 0 inc(5,2) = 1 inc(6,1) = -1 inc(6,2) = 1 inc(7,1) = -1 inc(7,2) = 0 inc(8,1) = -1 inc(8,2) = -1 ! ! ! Get the outflows and determine for each point to which outflow point it belong ! nbout = 0 trip_flow(:,:,:) = 0 DO ip=1,nbi DO jp=1,nbj IF ( trip(ip,jp) .GT. 9) THEN nbout = nbout + 1 outflow(nbout,1) = ip outflow(nbout,2) = jp ENDIF IF ( trip(ip,jp) .GT. 0) THEN trip_flow(ip,jp,1) = ip trip_flow(ip,jp,2) = jp ENDIF ENDDO ENDDO ! ! Follow the flow of the water ! DO ip=1,nbi DO jp=1,nbj IF ( (trip_flow(ip,jp,1) + trip_flow(ip,jp,2)) .GT. 0) THEN trp = trip(trip_flow(ip,jp,1), trip_flow(ip,jp,2)) cnt = 0 DO WHILE ( trp .GT. 0 .AND. trp .LT. 9 .AND. cnt .LT. nbi*nbj) cnt = cnt + 1 trip_flow(ip,jp,1) = trip_flow(ip,jp,1) + inc(trp,1) trip_flow(ip,jp,2) = trip_flow(ip,jp,2) + inc(trp,2) trp = trip(trip_flow(ip,jp,1), trip_flow(ip,jp,2)) ENDDO IF ( cnt .EQ. nbi*nbj) THEN WRITE(numout,*) 'We could not route point (routing_findrout) :', ip, jp WRITE(numout,*) '-------------> trip ' WRITE(fmt,"('(',I3,'I6)')") nbi DO ib=1,nbj WRITE(numout,fmt) trip(1:nbi,ib) ENDDO STOP ENDIF ENDIF ENDDO ENDDO ! ! What is the size of the region behind each outflow point ? ! totsz = 0 DO ip=1,nbout outsz(ip) = COUNT(trip_flow(:,:,1) .EQ. outflow(ip,1) .AND. trip_flow(:,:,2) .EQ. outflow(ip,2)) totsz = totsz + outsz(ip) ENDDO IF ( basin_sz .NE. totsz) THEN WRITE(numout,*) 'Water got lost while I tried to follow it ' WRITE(numout,*) basin_sz, totsz WRITE(numout,*) 'Basin id :', basinid DO ip=1,nbout WRITE(numout,*) 'ip :', ip, ' outsz :', outsz(ip), ' outflow :', outflow(ip,1), outflow(ip,2) ENDDO WRITE(numout,*) '-------------> trip ' WRITE(fmt,"('(',I3,'I6)')") nbi DO jp=1,nbj WRITE(numout,fmt) trip(1:nbi,jp) ENDDO STOP ENDIF ! END SUBROUTINE routing_findrout ! !---------------------------------------------------------------------------------------------- ! SUBROUTINE routing_globalize(nbpt, ib, neighbours, area_bx, trip_bx, hierarchy_bx, topoind_bx, min_topoind,& & nb_basin, basin_inbxid, basin_sz, basin_pts, basin_bxout, coast_pts, nwbas, basin_count,& & basin_area, basin_hierarchy, basin_topoind, basin_id, basin_flowdir, outflow_grid,& & nbcoastal, coastal_basin) ! IMPLICIT NONE ! ! This subroutine will put the basins found for grid-box in in the global map. Connection can ! only be made later when all information is together. ! ! ! One of the outputs is basin_flowdir. Its convention is 1-8 for the directions from North to North ! West going through South. The negative values will be -3 for return flow, -2 for coastal flow ! and -1 river flow. ! ! LOCAL ! INTEGER(i_std), INTENT (in) :: nbpt, ib ! Current grid-box INTEGER(i_std), INTENT(in) :: neighbours(nbpt,8) ! Vector of neighbours for each grid point ! (1=N, 2=E, 3=S, 4=W) REAL(r_std), DIMENSION(nbvmax,nbvmax) :: area_bx ! Area of each small box in the grid-box INTEGER(i_std), DIMENSION(nbvmax,nbvmax) :: trip_bx ! The trip field for each of the smaler boxes REAL(r_std), DIMENSION(nbvmax,nbvmax) :: hierarchy_bx ! level in the basin of the point REAL(r_std), DIMENSION(nbvmax,nbvmax) :: topoind_bx ! Topographic index REAL(r_std) :: min_topoind ! The default topographic index INTEGER(i_std) :: nb_basin ! number of sub-basins INTEGER(i_std), DIMENSION(nbvmax) :: basin_inbxid, basin_sz ! ID of basin, number of points in the basin INTEGER(i_std), DIMENSION(nbvmax,nbvmax,2):: basin_pts ! Points in each basin INTEGER(i_std), DIMENSION(nbvmax) :: basin_bxout ! outflow direction INTEGER(i_std) :: coast_pts(nbvmax) ! The coastal flow points ! global maps INTEGER(i_std) :: nwbas INTEGER(i_std), DIMENSION(nbpt) :: basin_count INTEGER(i_std), DIMENSION(nbpt,nwbas) :: basin_id, basin_flowdir REAL(r_std), DIMENSION(nbpt,nwbas) :: basin_area, basin_hierarchy, basin_topoind INTEGER(i_std), DIMENSION(nbpt,nwbas) :: outflow_grid INTEGER(i_std), DIMENSION(nbpt) :: nbcoastal INTEGER(i_std), DIMENSION(nbpt,nwbas) :: coastal_basin ! ! LOCAL ! INTEGER(i_std) :: ij, iz CHARACTER(LEN=4) :: hierar_method = 'OUTP' ! ! DO ij=1, nb_basin ! ! Count the basins and keep their ID ! basin_count(ib) = basin_count(ib)+1 if (basin_count(ib) > nwbas) then WRITE(numout,*) 'ib=',ib call ipslerr(3,'routing_globalize', & & 'Problem with basin_count : ', & & 'It is greater than number of allocated basin nwbas.', & & '(stop to count basins)') endif basin_id(ib,basin_count(ib)) = basin_inbxid(ij) ! ! Transfer the list of basins which flow into the ocean as coastal flow. ! IF ( basin_id(ib,basin_count(ib)) .LT. 0) THEN nbcoastal(ib) = basin_sz(ij) coastal_basin(ib,1:nbcoastal(ib)) = coast_pts(1:nbcoastal(ib)) ENDIF ! ! ! Compute the area of the basin ! basin_area(ib,ij) = 0.0 basin_hierarchy(ib,ij) = 0.0 ! SELECT CASE (hierar_method) ! CASE("MINI") basin_hierarchy(ib,ij) = undef_sechiba ! END SELECT basin_topoind(ib,ij) = 0.0 ! DO iz=1,basin_sz(ij) basin_area(ib,ij) = basin_area(ib,ij) + area_bx(basin_pts(ij,iz,1),basin_pts(ij,iz,2)) basin_topoind(ib,ij) = basin_topoind(ib,ij) + topoind_bx(basin_pts(ij,iz,1),basin_pts(ij,iz,2)) ! ! There are a number of ways to determine the hierarchy of the entire basin. ! We allow for three here : ! - Take the mean value ! - Take the minimum value within the basin ! - Take the value at the outflow point ! Probably taking the value of the outflow point is the best solution. ! SELECT CASE (hierar_method) ! CASE("MEAN") ! Mean hierarchy of the basin basin_hierarchy(ib,ij) = basin_hierarchy(ib,ij) + & & hierarchy_bx(basin_pts(ij,iz,1),basin_pts(ij,iz,2)) CASE("MINI") ! The smalest value of the basin IF ( hierarchy_bx(basin_pts(ij,iz,1),basin_pts(ij,iz,2)) .LT. basin_hierarchy(ib,ij)) THEN basin_hierarchy(ib,ij) = hierarchy_bx(basin_pts(ij,iz,1),basin_pts(ij,iz,2)) ENDIF CASE("OUTP") ! Value at the outflow point IF ( trip_bx(basin_pts(ij,iz,1),basin_pts(ij,iz,2)) .GT. 100 ) THEN basin_hierarchy(ib,ij) = hierarchy_bx(basin_pts(ij,iz,1),basin_pts(ij,iz,2)) ENDIF CASE DEFAULT WRITE(numout,*) 'Unknown method for computing the hierarchy of the basin' STOP 'routing_globalize' END SELECT ! ENDDO ! basin_topoind(ib,ij) = basin_topoind(ib,ij)/REAL(basin_sz(ij),r_std) ! SELECT CASE (hierar_method) ! CASE("MEAN") basin_hierarchy(ib,ij) = basin_hierarchy(ib,ij)/REAL(basin_sz(ij),r_std) ! END SELECT ! ! To make sure that it has the lowest number if this is an outflow point we reset basin_hierarchy ! IF (basin_bxout(ij) .LT. 0) THEN basin_hierarchy(ib,ij) = min_topoind basin_topoind(ib,ij) = min_topoind ENDIF ! ! ! Keep the outflow boxes and basin ! basin_flowdir(ib,ij) = basin_bxout(ij) IF (basin_bxout(ij) .GT. 0) THEN outflow_grid(ib,ij) = neighbours(ib,basin_bxout(ij)) ELSE outflow_grid(ib,ij) = basin_bxout(ij) ENDIF ! ! ENDDO ! ! END SUBROUTINE routing_globalize ! !----------------------------------------------------------------------------- ! SUBROUTINE routing_linkup(nbpt, neighbours, nwbas, basin_count, basin_area, basin_id, basin_flowdir, & & basin_hierarchy, outflow_grid, outflow_basin, inflow_number, inflow_grid, inflow_basin, nbcoastal,& & coastal_basin, invented_basins) ! IMPLICIT NONE ! ! This subroutine will make the connections between the basins and ensure global coherence. ! ! The convention for outflow_grid is : ! outflow_grid = -1 : River flow ! outflow_grid = -2 : Coastal flow ! outflow_grid = -3 : Return flow ! ! ! INPUT ! INTEGER(i_std), INTENT (in) :: nbpt INTEGER(i_std), DIMENSION(nbpt,8), INTENT (in) :: neighbours ! INTEGER(i_std) :: nwbas INTEGER(i_std), DIMENSION(nbpt) :: basin_count INTEGER(i_std), DIMENSION(nbpt,nwbas) :: basin_id INTEGER(i_std), DIMENSION(nbpt,nwbas) :: basin_flowdir REAL(r_std), DIMENSION(nbpt,nwbas) :: basin_area, basin_hierarchy INTEGER(i_std), DIMENSION(nbpt,nwbas) :: outflow_grid INTEGER(i_std), DIMENSION(nbpt,nwbas) :: outflow_basin INTEGER(i_std), DIMENSION(nbpt,nwbas) :: inflow_number INTEGER(i_std), DIMENSION(nbpt,nwbas,nbvmax) :: inflow_basin INTEGER(i_std), DIMENSION(nbpt,nwbas,nbvmax) :: inflow_grid INTEGER(i_std), DIMENSION(nbpt) :: nbcoastal INTEGER(i_std), DIMENSION(nbpt,nwbas) :: coastal_basin REAL(r_std), INTENT(in) :: invented_basins ! ! LOCAL ! INTEGER(i_std) :: sp, sb, sbl, inp, bid, outdm1, outdp1 INTEGER(i_std) :: dp1, dm1, dm1i, dp1i, bp1, bm1 INTEGER(i_std) :: dop, bop INTEGER(i_std) :: fbas(nwbas), nbfbas REAL(r_std) :: fbas_hierarchy(nwbas) INTEGER(i_std) :: ff(1) ! outflow_basin(:,:) = undef_int inflow_number(:,:) = 0 ! DO sp=1,nbpt DO sb=1,basin_count(sp) ! inp = outflow_grid(sp,sb) bid = basin_id(sp,sb) ! ! We only work on this point if it does not flow into the ocean ! At this point any of the outflows is designated by a negative values in ! outflow_grid ! IF ( inp .GT. 0 ) THEN ! ! Now find the basin in the onflow point (inp) ! nbfbas = 0 ! ! DO sbl=1,basin_count(inp) ! ! Either it is a standard basin or one aggregated from ocean flow points. ! If we flow into a another grid box we have to make sure that its hierarchy in the ! basin is lower. ! ! IF ( basin_id(inp,sbl) .GT. 0 ) THEN IF ( basin_id(inp,sbl) .EQ. bid .OR. basin_id(inp,sbl) .GT. invented_basins) THEN nbfbas =nbfbas + 1 fbas(nbfbas) = sbl fbas_hierarchy(nbfbas) = basin_hierarchy(inp,sbl) ENDIF ELSE IF ( COUNT(coastal_basin(inp,1:nbcoastal(inp)) .EQ. bid) .GT. 0 ) THEN nbfbas =nbfbas + 1 fbas(nbfbas) = sbl fbas_hierarchy(nbfbas) = basin_hierarchy(inp,sbl) ENDIF ENDIF ! ENDDO ! ! If we have more than one basin we will take the one which is lowest ! in the hierarchy. ! IF (nbfbas .GE. 1) THEN ff = MINLOC(fbas_hierarchy(1:nbfbas)) sbl = fbas(ff(1)) ! bop = undef_int IF ( basin_hierarchy(inp,sbl) .LE. basin_hierarchy(sp,sb) ) THEN IF ( basin_hierarchy(inp,sbl) .LE. basin_hierarchy(sp,sb) ) THEN bop = sbl ELSE ! The same hierarchy is allowed if both grids flow in about ! the same direction : IF ( ( MOD(basin_flowdir(inp,sbl)+1-1, 8)+1 .EQ. basin_flowdir(sp,sb)).OR. & & ( basin_flowdir(inp,sbl) .EQ. basin_flowdir(sp,sb)).OR. & & ( MOD(basin_flowdir(inp,sbl)+7-1, 8)+1 .EQ. basin_flowdir(sp,sb)) ) THEN bop = sbl ENDIF ENDIF ENDIF ! ! If the basin is suitable (bop < undef_int) then take it ! IF ( bop .LT. undef_int ) THEN outflow_basin(sp,sb) = bop inflow_number(inp,bop) = inflow_number(inp,bop) + 1 IF ( inflow_number(inp,bop) .LE. nbvmax ) THEN inflow_grid(inp, bop, inflow_number(inp,bop)) = sp inflow_basin(inp, bop, inflow_number(inp,bop)) = sb ELSE WRITE(numout,*) 'Increase nbvmax' STOP 'routing_linkup' ENDIF ENDIF ENDIF ! ! ENDIF ! ! ! ! Did we find it ? ! ! In case the outflow point was ocean or we did not find the correct basin we start to look ! around. We find two options for the outflow direction (dp1 & dm1) and the corresponding ! basin index (bp1 & bm1). ! ! IF ( outflow_basin(sp,sb) .EQ. undef_int & & .AND. basin_flowdir(sp,sb) .GT. 0) THEN ! dp1i = MOD(basin_flowdir(sp,sb)+1-1, 8)+1 dp1 = neighbours(sp,dp1i) dm1i = MOD(basin_flowdir(sp,sb)+7-1, 8)+1 IF ( dm1i .LT. 1 ) dm1i = 8 dm1 = neighbours(sp,dm1i) ! ! bp1 = -1 IF ( dp1 .GT. 0 ) THEN DO sbl=1,basin_count(dp1) IF (basin_id(dp1,sbl) .EQ. bid .AND.& & basin_hierarchy(sp,sb) .GE. basin_hierarchy(dp1,sbl) .AND. & & bp1 .LT. 0) THEN IF ( basin_hierarchy(sp,sb) .GT. basin_hierarchy(dp1,sbl) ) THEN bp1 = sbl ELSE ! The same hierarchy is allowed if both grids flow in about ! the same direction : IF ( ( MOD(basin_flowdir(dp1,sbl)+1-1, 8)+1 .EQ. basin_flowdir(sp,sb)).OR. & & ( basin_flowdir(dp1,sbl) .EQ. basin_flowdir(sp,sb)).OR. & & ( MOD(basin_flowdir(dp1,sbl)+7-1, 8)+1 .EQ. basin_flowdir(sp,sb)) ) THEN bp1 = sbl ENDIF ENDIF ENDIF ENDDO ENDIF ! bm1 = -1 IF ( dm1 .GT. 0 ) THEN DO sbl=1,basin_count(dm1) IF (basin_id(dm1,sbl) .EQ. bid .AND.& & basin_hierarchy(sp,sb) .GE. basin_hierarchy(dm1,sbl) .AND. & & bm1 .LT. 0) THEN IF ( basin_hierarchy(sp,sb) .GT. basin_hierarchy(dm1,sbl) ) THEN bm1 = sbl ELSE ! The same hierarchy is allowed if both grids flow in about ! the same direction : IF ( ( MOD(basin_flowdir(dm1,sbl)+1-1, 8)+1 .EQ. basin_flowdir(sp,sb)) .OR. & & ( basin_flowdir(dm1,sbl) .EQ. basin_flowdir(sp,sb)) .OR. & & ( MOD(basin_flowdir(dm1,sbl)+7-1, 8)+1 .EQ. basin_flowdir(sp,sb)) ) THEN bm1 = sbl ENDIF ENDIF ENDIF ENDDO ENDIF ! ! ! First deal with the case on land. ! ! For that we need to check if the water will be able to flow out of the grid dp1 or dm1 ! and not return to our current grid. If it is the current grid ! then we can not do anything with that neighbour. Thus we set the ! value of outdm1 and outdp1 back to -1 ! outdp1 = undef_int IF ( dp1 .GT. 0 .AND. bp1 .GT. 0 ) THEN ! if the outflow is into the ocean then we put something less than undef_int in outdp1! IF (basin_flowdir(dp1,bp1) .GT. 0) THEN outdp1 = neighbours(dp1,basin_flowdir(dp1,bp1)) IF ( outdp1 .EQ. sp ) outdp1 = undef_int ELSE outdp1 = nbpt + 1 ENDIF ENDIF outdm1 = undef_int IF ( dm1 .GT. 0 .AND. bm1 .GT. 0 ) THEN IF (basin_flowdir(dm1,bm1) .GT. 0) THEN outdm1 = neighbours(dm1,basin_flowdir(dm1,bm1)) IF ( outdm1 .EQ. sp ) outdm1 = undef_int ELSE outdm1 = nbpt + 1 ENDIF ENDIF ! ! Now that we know our options we need go through them. ! dop = undef_int bop = undef_int IF ( outdp1 .LT. undef_int .AND. outdm1 .LT. undef_int) THEN ! ! In this case we let the current basin flow into the smaller one ! IF ( basin_area(dp1,bp1) .LT. basin_area(dm1,bm1) ) THEN dop = dp1 bop = bp1 ELSE dop = dm1 bop = bm1 ENDIF ! ! ELSE IF ( outdp1 .LT. undef_int ) THEN ! If only the first one is possible dop = dp1 bop = bp1 ELSE IF ( outdm1 .LT. undef_int ) THEN ! If only the second one is possible dop = dm1 bop = bm1 ELSE ! ! Now we are at the point where none of the neighboring points is suitable ! or we have a coastal point. ! ! If there is an option to put the water into the ocean go for it. ! IF ( outflow_grid(sp,sb) .LT. 0 .OR. dm1 .LT. 0 .OR. dp1 .LT. 0 ) THEN dop = -1 ELSE ! ! If we are on a land point with only land neighbors but no one suitable to let the ! water flow into we have to look for a solution in the current grid box. ! ! IF ( bp1 .LT. 0 .AND. bm1 .LT. 0 ) THEN ! ! Do we have more than one basin with the same ID ? ! IF ( COUNT(basin_id(sp,1:basin_count(sp)) .EQ. bid) .GE. 2) THEN ! ! Now we can try the option of flowing into the basin of the same grid box. ! DO sbl=1,basin_count(sp) IF (sbl .NE. sb .AND. basin_id(sp,sbl) .EQ. bid) THEN ! In case this basin has a lower hierarchy or flows into a totaly ! different direction we go for it. IF ( (basin_hierarchy(sp,sb) .GE. basin_hierarchy(sp,sbl)) .OR. & & (basin_flowdir(sp,sbl) .LT. dm1i .AND.& & basin_flowdir(sp,sbl) .GT. dp1i) ) THEN dop = sp bop = sbl IF (basin_hierarchy(sp,sb) .LT. basin_hierarchy(sp,sbl)) THEN WRITE(numout,*) '>>>>>>> POINT CORRECTED against hierarchy :',& & sp, sb, 'into', sbl ENDIF ENDIF ! ENDIF ENDDO ! ENDIF ENDIF ENDIF ! IF ( dop .EQ. undef_int .AND. bop .EQ. undef_int ) THEN WRITE(numout,*) 'Why are we here with point ', sp, sb WRITE(numout,*) 'Coodinates : (lon,lat) = ', lalo(sp,2), lalo(sp,1) WRITE(numout,*) 'Contfrac : = ', contfrac(sp) WRITE(numout,*) 'Local Basin ID :', basin_id(sp,1:basin_count(sp)) WRITE(numout,*) 'Local hierarchies :', basin_hierarchy(sp,1:basin_count(sp)) WRITE(numout,*) 'Local basin_flowdir :', basin_flowdir(sp,1:basin_count(sp)) WRITE(numout,*) 'Local outflowgrid :', outflow_grid(sp,1:basin_count(sp)) WRITE(numout,*) 'outflow_grid :', inp WRITE(numout,*) 'Coodinates outflow : (lon,lat) = ', lalo(inp,2), lalo(inp,1) WRITE(numout,*) 'Contfrac : = ', contfrac(inp) WRITE(numout,*) 'Outflow Basin ID :', basin_id(inp,1:basin_count(inp)) WRITE(numout,*) 'Outflow hierarchies :', basin_hierarchy(inp,1:basin_count(inp)) WRITE(numout,*) 'Outflow basin_flowdir :', basin_flowdir(inp,1:basin_count(inp)) WRITE(numout,*) 'Explored options +1 :', dp1, bp1, outdp1 WRITE(numout,*) 'Explored +1 Basin ID :', basin_id(dp1,1:basin_count(dp1)) WRITE(numout,*) 'Explored +1 hierarchies :', basin_hierarchy(dp1,1:basin_count(dp1)) WRITE(numout,*) 'Explored +1 basin_flowdir :', basin_flowdir(dp1,1:basin_count(dp1)) WRITE(numout,*) 'Explored options -1 :', dm1, bm1, outdm1 WRITE(numout,*) 'Explored -1 Basin ID :', basin_id(dm1,1:basin_count(dm1)) WRITE(numout,*) 'Explored -1 hierarchies :', basin_hierarchy(dm1,1:basin_count(dm1)) WRITE(numout,*) 'Explored -1 basin_flowdir :', basin_flowdir(dm1,1:basin_count(dm1)) WRITE(numout,*) '****************************' IF ( contfrac(sp) > 0.01 ) THEN CALL ipslerr(3,'routing_linkup', & & 'In the routine which make connections between the basins and ensure global coherence,', & & 'there is a problem with outflow linkup without any valid direction.', & & '(Perhaps there is a problem with the grid.)') ENDIF ENDIF ! ENDIF ! ! Now that we know where we want the water to flow to we write the ! the information in the right fields. ! IF ( dop .GT. 0 .AND. dop .NE. undef_int ) THEN outflow_grid(sp,sb) = dop outflow_basin(sp,sb) = bop inflow_number(dop,bop) = inflow_number(dop,bop) + 1 IF ( inflow_number(dop,bop) .LE. nbvmax ) THEN inflow_grid(dop, bop, inflow_number(dop,bop)) = sp inflow_basin(dop, bop, inflow_number(dop,bop)) = sb ELSE WRITE(numout,*) 'Increase nbvmax' STOP 'routing_linkup' ENDIF ! ELSE outflow_grid(sp,sb) = -2 outflow_basin(sp,sb) = undef_int ENDIF ! ENDIF ! ! ! If we still have not found anything then we have to check that there is not a basin ! within the same grid box which has a lower hierarchy. ! ! IF ( outflow_grid(sp,sb) .GT. 0 .AND. outflow_basin(sp,sb) .EQ. undef_int & & .AND. basin_flowdir(sp,sb) .GT. 0) THEN ! WRITE(numout,*) 'There is no reason to here, this part of the code should be dead :', sp,sb ! DO sbl=1,basin_count(sp) ! ! Three conditons are needed to let the water flow into another basin of the ! same grid : ! - another basin than the current one ! - same ID ! - of lower hierarchy. ! IF ( (sbl .NE. sb) .AND. (basin_id(sp,sbl) .EQ. bid)& & .AND. (basin_hierarchy(sp,sb) .GT. basin_hierarchy(sp,sbl)) ) THEN outflow_basin(sp,sb) = sbl inflow_number(sp,sbl) = inflow_number(sp,sbl) + 1 IF ( inflow_number(sp,sbl) .LE. nbvmax ) THEN IF ( sp .EQ. 42 .AND. sbl .EQ. 1) THEN WRITE(numout,*) 'ADD INFLOW (3):', sp, sb ENDIF inflow_grid(sp, sbl, inflow_number(sp,sbl)) = sp inflow_basin(sp, sbl, inflow_number(sp,sbl)) = sb ELSE WRITE(numout,*) 'Increase nbvmax' STOP 'routing_linkup' ENDIF ENDIF ENDDO ENDIF ! ! Ok that is it, we give up :-) ! IF ( outflow_grid(sp,sb) .GT. 0 .AND. outflow_basin(sp,sb) .EQ. undef_int & & .AND. basin_flowdir(sp,sb) .GT. 0) THEN ! WRITE(numout,*) 'We could not find the basin into which we need to flow' WRITE(numout,*) 'Grid point ', sp, ' and basin ', sb WRITE(numout,*) 'Explored neighbours :', dm1, dp1 WRITE(numout,*) 'Outflow direction :', basin_flowdir(sp,sb) WRITE(numout,*) 'Outlfow grid :', outflow_grid(sp,sb) WRITE(numout,*) 'basin ID:',basin_id(sp,sb) WRITE(numout,*) 'Hierarchy :', basin_hierarchy(sp,sb) STOP 'routing_linkup' ENDIF ENDDO ! ENDDO ! ! Check for each outflow basin that it exists ! DO sp=1,nbpt DO sb=1,basin_count(sp) ! inp = outflow_grid(sp,sb) sbl = outflow_basin(sp,sb) IF ( inp .GE. 0 ) THEN IF ( basin_count(inp) .LT. sbl ) THEN WRITE(numout,*) 'point :', sp, ' basin :', sb WRITE(numout,*) 'Flows into point :', inp, ' basin :', sbl WRITE(numout,*) 'But it flows into now here as basin count = ', basin_count(inp) STOP ENDIF ENDIF ENDDO ENDDO ! END SUBROUTINE routing_linkup ! !--------------------------------------------------------------------------- ! SUBROUTINE routing_fetch(nbpt, resolution, contfrac, nwbas, basin_count, basin_area, basin_id,& & outflow_grid, outflow_basin, fetch_basin) ! IMPLICIT NONE ! ! This subroutine will compute the fetch of each basin. This means that for each basin we ! will know how much area is upstream. It will help decide how to procede with the ! the truncation later and allow to set correctly in outflow_grid the distinction ! between coastal and river flow ! ! ! INPUT ! INTEGER(i_std), INTENT(in) :: nbpt ! REAL(r_std), DIMENSION(nbpt,2), INTENT(in) :: resolution REAL(r_std), DIMENSION(nbpt), INTENT(in) :: contfrac ! Fraction of land in each grid box ! INTEGER(i_std) :: nwbas INTEGER(i_std), DIMENSION(nbpt), INTENT(in) :: basin_count REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: basin_area INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(in) :: basin_id INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: outflow_grid INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(in) :: outflow_basin REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(out) :: fetch_basin ! ! LOCAL ! INTEGER(i_std) :: ib, ij, ff(1), it, itt, igrif, ibasf, nboutflow REAL(r_std) :: contarea, totbasins REAL(r_std), DIMENSION(nbpt*nbvmax) :: tmp_area INTEGER(i_std), DIMENSION(nbpt*nbvmax,2) :: tmpindex ! ! ! Normalize the area of all basins ! DO ib=1,nbpt ! totbasins = SUM(basin_area(ib,1:basin_count(ib))) contarea = resolution(ib,1)*resolution(ib,2)*contfrac(ib) ! DO ij=1,basin_count(ib) basin_area(ib,ij) = basin_area(ib,ij)/totbasins*contarea ENDDO ! ENDDO WRITE(numout,*) 'Normalization done' ! ! Compute the area upstream of each basin ! fetch_basin(:,:) = 0.0 ! ! DO ib=1,nbpt ! DO ij=1,basin_count(ib) ! fetch_basin(ib, ij) = fetch_basin(ib, ij) + basin_area(ib,ij) ! igrif = outflow_grid(ib,ij) ibasf = outflow_basin(ib,ij) ! itt = 0 DO WHILE (igrif .GT. 0) fetch_basin(igrif,ibasf) = fetch_basin(igrif,ibasf) + basin_area(ib, ij) it = outflow_grid(igrif, ibasf) ibasf = outflow_basin(igrif, ibasf) igrif = it itt = itt + 1 IF ( itt .GT. 500) THEN WRITE(numout,& "('Grid ',I5, ' and basin ',I5, 'did not converge after iteration ',I5)") ib, ij, itt WRITE(numout,*) 'Basin ID :', basin_id(igrif,ibasf) WRITE(numout,& "('We are stuck with the flow into grid ',I5,' and basin ',I5)") igrif, ibasf IF ( itt .GT. 510) THEN STOP ' routing_fetch' ENDIF ENDIF ENDDO ! ENDDO ! ENDDO ! WRITE(numout,*) 'The smalest FETCH :', MINVAL(fetch_basin) WRITE(numout,*) 'The largest FETCH :', MAXVAL(fetch_basin) ! ! Now we set for the 'num_largest' largest basins the outflow condition as stream flow ! (i.e. outflow_grid = -1) and all other outflow basins are set to coastal flow ! (i.e. outflow_grid = -2). The return flow is not touched. ! nboutflow = 0 ! DO ib=1,nbpt ! DO ij=1,basin_count(ib) ! ! We do not need any more the river flow flag as we are going to reset it. ! IF ( outflow_grid(ib,ij) .EQ. -1) THEN outflow_grid(ib,ij) = -2 ENDIF ! IF ( outflow_grid(ib,ij) .EQ. -2) THEN ! nboutflow = nboutflow + 1 tmp_area(nboutflow) = fetch_basin(ib,ij) tmpindex(nboutflow,1) = ib tmpindex(nboutflow,2) = ij ! ENDIF ! ENDDO ENDDO ! DO ib=1, num_largest ff = MAXLOC(tmp_area(1:nboutflow)) outflow_grid(tmpindex(ff(1),1), tmpindex(ff(1),2)) = -1 tmp_area(ff(1)) = 0.0 ENDDO ! END SUBROUTINE routing_fetch ! !--------------------------------------------------------------------------- ! SUBROUTINE routing_truncate(nbpt, resolution, contfrac, nwbas, basin_count, basin_area, basin_topoind,& & fetch_basin, basin_id, basin_flowdir, outflow_grid, outflow_basin, inflow_number,& & inflow_grid, inflow_basin) ! IMPLICIT NONE ! ! ! This subroutine reduces the number of basins per gird to the value chosen by the user. ! it also computes the final field which will be used to route the water at the ! requested truncation. ! ! INPUT ! INTEGER(i_std) :: nbpt ! REAL(r_std), DIMENSION(nbpt,2) :: resolution REAL(r_std), DIMENSION(nbpt), INTENT(in) :: contfrac ! Fraction of land in each grid box ! INTEGER(i_std) :: nwbas INTEGER(i_std), DIMENSION(nbpt) :: basin_count INTEGER(i_std), DIMENSION(nbpt,nwbas) :: basin_id INTEGER(i_std), DIMENSION(nbpt,nwbas) :: basin_flowdir REAL(r_std), DIMENSION(nbpt,nwbas) :: basin_area REAL(r_std), DIMENSION(nbpt,nwbas) :: basin_topoind REAL(r_std), DIMENSION(nbpt,nwbas) :: fetch_basin INTEGER(i_std), DIMENSION(nbpt,nwbas) :: outflow_grid INTEGER(i_std), DIMENSION(nbpt,nwbas) :: outflow_basin INTEGER(i_std), DIMENSION(nbpt,nwbas) :: inflow_number INTEGER(i_std), DIMENSION(nbpt,nwbas,nwbas) :: inflow_basin INTEGER(i_std), DIMENSION(nbpt,nwbas,nwbas) :: inflow_grid ! ! LOCAL ! INTEGER(i_std), PARAMETER :: pickmax = 200 INTEGER(i_std) :: ib, ij, ibf, ijf, igrif, ibasf, cnt, pold, bold, ff(2) INTEGER(i_std) :: ii, kbas, sbas, ik, iter, ibt, obj REAL(r_std), DIMENSION(nbpt,nbasmax) :: floflo REAL(r_std), DIMENSION(nbpt) :: gridarea, gridbasinarea REAL(r_std) :: ratio INTEGER(i_std), DIMENSION(pickmax,2) :: largest_basins INTEGER(i_std), DIMENSION(pickmax) :: tmp_ids INTEGER(i_std) :: multbas, iml(1) INTEGER(i_std), DIMENSION(pickmax) :: multbas_sz REAL(r_std), DIMENSION(pickmax) :: tmp_area INTEGER(i_std), DIMENSION(pickmax,pickmax) :: multbas_list ! INTEGER(i_std) :: nbtruncate INTEGER(i_std), SAVE, ALLOCATABLE, DIMENSION(:) :: indextrunc ! IF ( .NOT. ALLOCATED(indextrunc)) THEN ALLOCATE(indextrunc(nbpt)) ENDIF ! ! Truncate if needed and find the path closest to the high resolution data. ! ! The algorithm : ! ! We only go through this procedure only as many times as there are basins to take out at most. ! This is important as it allows the simplifications to spread from one grid to the other. ! The for each step of the iteration and at each grid point we check the following options for ! simplifying the pathways of water : ! 1) If the basin of a grid flows into another basin of the same grid. Kill the one which only ! served as a transition ! 2) If in one grid box we have a number of basins which flow into the ocean as coastal flow. ! We kill the smallest one and put into the largest basin. There is no need to manage many ! basins going into the ocean as coastal flows. ! 3) If we have streams run in parallel from one gird box to the others (that is these are ! different basins) we will put the smaller one in the larger one. This may hapen at any ! level of the flow but in theory it should propagate downstream. ! 4) If we have two basins with the same ID but flow into different grid boxes we sacrifice ! the smallest one and route it through the largest. ! ! Obviously if any of the options find something then we skip the rest and take out the basin. ! ! ! We have to go through the grid as least as often as we have to reduce the number of basins ! For good measure we add 3 more passages. ! ! DO iter = 1, MAXVAL(basin_count) - nbasmax +3 ! ! Get the points over which we wish to truncate ! nbtruncate = 0 DO ib = 1, nbpt IF ( basin_count(ib) .GT. nbasmax ) THEN nbtruncate = nbtruncate + 1 indextrunc(nbtruncate) = ib ENDIF ENDDO ! ! Go through the basins which need to be truncated. ! DO ibt=1,nbtruncate ! ib = indextrunc(ibt) ! ! Check if we have basin which flows into a basin in the same grid ! kbas = basin we will have to kill ! sbas = basin which takes over kbas ! kbas = 0 sbas = 0 ! ! 1) Can we find a basin which flows into a basin of the same grid ? ! DO ij=1,basin_count(ib) DO ii=1,basin_count(ib) IF ( outflow_grid(ib,ii) .EQ. ib .AND. outflow_basin(ib, ii) .EQ. ij .AND. kbas*sbas .NE. 0) THEN kbas = ii sbas = ij ENDIF ENDDO ENDDO ! ! 2) Merge two basins which flow into the ocean as coastal or return flow ! (outflow_grid = -2 or -3). Well obviously only if we have more than 1 and ! have not found anything yet! ! IF ( (COUNT(outflow_grid(ib,1:basin_count(ib)) .EQ. -2) .GT. 1 .OR.& & COUNT(outflow_grid(ib,1:basin_count(ib)) .EQ. -3) .GT. 1) .AND.& & kbas*sbas .EQ. 0) THEN ! multbas = 0 multbas_sz(:) = 0 ! IF ( COUNT(outflow_grid(ib,1:basin_count(ib)) .EQ. -2) .GT. 1 ) THEN obj = -2 ELSE obj = -3 ENDIF ! ! First we get the list of all basins which go out as coastal or return flow (obj) ! DO ij=1,basin_count(ib) IF ( outflow_grid(ib,ij) .EQ. obj ) THEN multbas = multbas + 1 multbas_sz(multbas) = ij tmp_area(multbas) = fetch_basin(ib,ij) ENDIF ENDDO ! ! Now the take the smalest to be transfered to the largest ! iml = MAXLOC(tmp_area(1:multbas), MASK = tmp_area(1:multbas) .GT. 0.) sbas = multbas_sz(iml(1)) iml = MINLOC(tmp_area(1:multbas), MASK = tmp_area(1:multbas) .GT. 0.) kbas = multbas_sz(iml(1)) ! ENDIF ! ! 3) If we have basins flowing into the same grid but different basins then we put them ! together. Obviously we first work with the grid which has most streams runing into it ! and puting the smalest in the largests catchments. ! IF ( kbas*sbas .EQ. 0) THEN ! tmp_ids(1:basin_count(ib)) = outflow_grid(ib,1:basin_count(ib)) multbas = 0 multbas_sz(:) = 0 ! ! First obtain the list of basins which flow into the same basin ! DO ij=1,basin_count(ib) IF ( outflow_grid(ib,ij) .GT. 0 .AND.& & COUNT(tmp_ids(1:basin_count(ib)) .EQ. outflow_grid(ib,ij)) .GT. 1) THEN multbas = multbas + 1 DO ii=1,basin_count(ib) IF ( tmp_ids(ii) .EQ. outflow_grid(ib,ij)) THEN multbas_sz(multbas) = multbas_sz(multbas) + 1 multbas_list(multbas,multbas_sz(multbas)) = ii tmp_ids(ii) = -99 ENDIF ENDDO ELSE tmp_ids(ij) = -99 ENDIF ENDDO ! ! Did we come up with any basins to deal with this way ? ! IF ( multbas .GT. 0 ) THEN ! iml = MAXLOC(multbas_sz(1:multbas)) ik = iml(1) ! ! Take the smalest and largest of these basins ! ! DO ii=1,multbas_sz(ik) tmp_area(ii) = fetch_basin(ib, multbas_list(ik,ii)) ENDDO iml = MAXLOC(tmp_area(1:multbas_sz(ik)), MASK = tmp_area(1:multbas_sz(ik)) .GT. 0.) sbas = multbas_list(ik,iml(1)) iml = MINLOC(tmp_area(1:multbas_sz(ik)), MASK = tmp_area(1:multbas_sz(ik)) .GT. 0.) kbas = multbas_list(ik,iml(1)) ! ENDIF ! ENDIF ! ! 4) If we have twice the same basin we put them together even if they flow into different ! directions. If one of them goes to the ocean it takes the advantage. ! IF ( kbas*sbas .EQ. 0) THEN ! tmp_ids(1:basin_count(ib)) = basin_id(ib,1:basin_count(ib)) multbas = 0 multbas_sz(:) = 0 ! ! First obtain the list of basins which have sub-basins in this grid box. ! (these are identified by their IDs) ! DO ij=1,basin_count(ib) IF ( COUNT(tmp_ids(1:basin_count(ib)) .EQ. basin_id(ib,ij)) .GT. 1) THEN multbas = multbas + 1 DO ii=1,basin_count(ib) IF ( tmp_ids(ii) .EQ. basin_id(ib,ij)) THEN multbas_sz(multbas) = multbas_sz(multbas) + 1 multbas_list(multbas,multbas_sz(multbas)) = ii tmp_ids(ii) = -99 ENDIF ENDDO ELSE tmp_ids(ij) = -99 ENDIF ENDDO ! ! We are going to work on the basin with the largest number of sub-basins. ! (IF we have a basin which has subbasins !) ! IF ( multbas .GT. 0 ) THEN ! iml = MAXLOC(multbas_sz(1:multbas)) ik = iml(1) ! ! If one of the basins goes to the ocean then it is going to have the priority ! tmp_area(:) = 0. IF ( COUNT(outflow_grid(ib,multbas_list(ik,1:multbas_sz(ik))) .LT. 0) .GT. 0) THEN DO ii=1,multbas_sz(ik) IF ( outflow_grid(ib,multbas_list(ik,ii)) .LT. 0 .AND. sbas .EQ. 0 ) THEN sbas = multbas_list(ik,ii) ELSE tmp_area(ii) = fetch_basin(ib, multbas_list(ik,ii)) ENDIF ENDDO ! take the smalest of the subbasins iml = MINLOC(tmp_area(1:multbas_sz(ik)), MASK = tmp_area(1:multbas_sz(ik)) .GT. 0.) kbas = multbas_list(ik,iml(1)) ELSE ! ! Else we take simply the largest and smalest ! DO ii=1,multbas_sz(ik) tmp_area(ii) = fetch_basin(ib, multbas_list(ik,ii)) ENDDO iml = MAXLOC(tmp_area(1:multbas_sz(ik)), MASK = tmp_area(1:multbas_sz(ik)) .GT. 0.) sbas = multbas_list(ik,iml(1)) iml = MINLOC(tmp_area(1:multbas_sz(ik)), MASK = tmp_area(1:multbas_sz(ik)) .GT. 0.) kbas = multbas_list(ik,iml(1)) ! ENDIF ! ! ENDIF ENDIF ! ! ! ! Then we call routing_killbas to clean up the basins in this grid ! IF ( kbas .GT. 0 .AND. sbas .GT. 0 ) THEN CALL routing_killbas(nbpt, ib, kbas, sbas, nwbas, basin_count, basin_area, basin_topoind,& & fetch_basin, basin_id, basin_flowdir, outflow_grid, outflow_basin, inflow_number,& & inflow_grid, inflow_basin) ENDIF ! ENDDO ! ! ENDDO ! ! If there are any grids left with too many basins we need to take out the big hammer ! ! We will only do it if this represents less than 5% of all points. ! IF ( COUNT(basin_count .GT. nbasmax) .GT. 0 ) THEN ! ! IF ( COUNT(basin_count .GT. nbasmax)/nbpt*100 .GT. 5 ) THEN WRITE(numout,*) 'We have ', COUNT(basin_count .GT. nbasmax)/nbpt*100, '% of all points which do not yet' WRITE(numout,*) 'have the right trunctaction. That is too much to apply a brutal method' DO ib = 1, nbpt IF ( basin_count(ib) .GT. nbasmax ) THEN ! WRITE(numout,*) 'We did not find a basin which could be supressed. We will' WRITE(numout,*) 'not be able to reduce the truncation in grid ', ib DO ij=1,basin_count(ib) WRITE(numout,*) 'grid, basin nb and id :', ib, ij, basin_id(ib,ij) WRITE(numout,*) 'Outflow grid and basin ->', outflow_grid(ib,ij), outflow_basin(ib, ij) ENDDO ENDIF ENDDO STOP 'routing_truncate' ! ELSE ! ! DO ib = 1,nbpt DO WHILE ( basin_count(ib) .GT. nbasmax ) ! WRITE(numout,*) 'HAMMER, ib, basin_count :', ib, basin_count(ib) ! ! Here we simply put the smallest basins into the largest ones. It is really a brute force ! method but it will only be applied if everything has failed. ! DO ii = 1,basin_count(ib) tmp_area(ii) = fetch_basin(ib, ii) ENDDO ! iml = MAXLOC(tmp_area(1:basin_count(ib)), MASK = tmp_area(1:basin_count(ib)) .GT. 0.) sbas =iml(1) iml = MINLOC(tmp_area(1:basin_count(ib)), MASK = tmp_area(1:basin_count(ib)) .GT. 0.) kbas = iml(1) ! IF ( kbas .GT. 0 .AND. sbas .GT. 0 ) THEN CALL routing_killbas(nbpt, ib, kbas, sbas, nwbas, basin_count, basin_area, basin_topoind,& & fetch_basin, basin_id, basin_flowdir, outflow_grid, outflow_basin, inflow_number,& & inflow_grid, inflow_basin) ENDIF ENDDO ENDDO ! ENDIF ! ! ENDIF ! ! Now that we have reached the right truncation (resolution) we will start ! to produce the variables we will use to route the water. ! DO ib=1,nbpt ! ! For non existing basins the route_tobasin variable is put to zero. This will allow us ! to pick up the number of basin afterwards. ! route_togrid(ib,:) = ib route_tobasin(ib,:) = 0 routing_area(ib,:) = 0.0 ! ENDDO ! ! Transfer the info into the definitive variables ! DO ib=1,nbpt DO ij=1,basin_count(ib) routing_area(ib,ij) = basin_area(ib,ij) topo_resid(ib,ij) = basin_topoind(ib,ij) global_basinid(ib,ij) = basin_id(ib,ij) route_togrid(ib,ij) = outflow_grid(ib,ij) route_tobasin(ib,ij) = outflow_basin(ib,ij) ENDDO ENDDO ! ! ! Set the new convention for the outflow conditions ! Now it is based in the outflow basin and the outflow grid will ! be the same as the current. ! returnflow to the grid : nbasmax + 1 ! coastal flow : nbasmax + 2 ! river outflow : nbasmax + 3 ! ! Here we put everything here in coastal flow. It is later where we will ! put the largest basins into river outflow. ! DO ib=1,nbpt DO ij=1,basin_count(ib) ! River flows IF ( route_togrid(ib,ij) .EQ. -1 ) THEN route_tobasin(ib,ij) = nbasmax + 2 route_togrid(ib,ij) = ib ! Coastal flows ELSE IF ( route_togrid(ib,ij) .EQ. -2 ) THEN route_tobasin(ib,ij) = nbasmax + 2 route_togrid(ib,ij) = ib ! Return flow ELSE IF ( route_togrid(ib,ij) .EQ. -3 ) THEN route_tobasin(ib,ij) = nbasmax + 1 route_togrid(ib,ij) = ib ENDIF ENDDO ENDDO ! ! A second check on the data. Just make sure that each basin flows somewhere. ! DO ib=1,nbpt DO ij=1,basin_count(ib) ibf = route_togrid(ib,ij) ijf = route_tobasin(ib,ij) IF ( ijf .GT. basin_count(ibf) .AND. ijf .LE. nbasmax) THEN WRITE(numout,*) 'Second check' WRITE(numout,*) 'point :', ib, ' basin :', ij WRITE(numout,*) 'Flows into point :', ibf, ' basin :', ijf WRITE(numout,*) 'But it flows into now here as basin count = ', basin_count(ibf) STOP ENDIF ENDDO ENDDO ! ! Verify areas of the contienents ! floflo(:,:) = 0.0 gridarea(:) = contfrac(:)*resolution(:,1)*resolution(:,2) DO ib=1,nbpt gridbasinarea(ib) = SUM(routing_area(ib,:)) ENDDO ! DO ib=1,nbpt DO ij=1,basin_count(ib) cnt = 0 igrif = ib ibasf = ij DO WHILE (ibasf .LE. nbasmax .AND. cnt .LT. nbasmax*nbpt) cnt = cnt + 1 pold = igrif bold = ibasf igrif = route_togrid(pold, bold) ibasf = route_tobasin(pold, bold) IF ( ibasf .GT. basin_count(igrif) .AND. ibasf .LE. nbasmax) THEN WRITE(numout,*) 'We should not be here as the basin flows into the pampa' WRITE(numout,*) 'Last correct point :', pold, bold WRITE(numout,*) 'It pointed to in the new variables :', route_togrid(pold, bold),route_tobasin(pold, bold) WRITE(numout,*) 'The old variables gave :', outflow_grid(pold, bold), outflow_basin(pold, bold) WRITE(numout,*) 'Where we ended up :', igrif,ibasf STOP ENDIF ENDDO ! IF ( ibasf .GT. nbasmax ) THEN floflo(igrif,bold) = floflo(igrif,bold) + routing_area(ib,ij) ELSE WRITE(numout,*) 'The flow did not end up in the ocean or in the grid cell.' WRITE(numout,*) 'For grid ', ib, ' and basin ', ij WRITE(numout,*) 'The last grid was ', igrif, ' and basin ', ibasf STOP 'routing_truncate' ENDIF ENDDO ENDDO ! DO ib=1,nbpt IF ( gridbasinarea(ib) > 0. ) THEN ratio = gridarea(ib)/gridbasinarea(ib) routing_area(ib,:) = routing_area(ib,:)*ratio ENDIF ENDDO ! WRITE(numout,*) 'Sum of area of all outflow areas :',SUM(routing_area) WRITE(numout,*) 'Surface of all continents :', SUM(gridarea) ! ! Redo the the distinction between river outflow and coastal flow. We can not ! take into account the return flow points. ! ibf = 0 DO ib=1, pickmax ff = MAXLOC(floflo) IF ( route_tobasin(ff(1), ff(2)) .EQ. nbasmax + 2) THEN ibf = ibf + 1 largest_basins(ibf,:) = ff(:) ENDIF floflo(ff(1), ff(2)) = 0.0 ENDDO ! ! Put the largest basins into river flows. ! IF ( ibf .LT. num_largest) THEN WRITE(numout,*) 'Not enough basins to choose the ', num_largest, 'largest' STOP 'routing_truncate' ENDIF ! ! ! DO ib=1, num_largest route_tobasin(largest_basins(ib,1),largest_basins(ib,2)) = nbasmax + 3 ENDDO ! WRITE(numout,*) 'NUMBER OF RIVERS :', COUNT(route_tobasin .GE. nbasmax + 3) ! END SUBROUTINE routing_truncate ! !------------------------------------------------------------------------------------- ! SUBROUTINE routing_killbas(nbpt, ib, tokill, totakeover, nwbas, basin_count, basin_area, basin_topoind,& & fetch_basin, basin_id, basin_flowdir, outflow_grid, outflow_basin, inflow_number,& & inflow_grid, inflow_basin) ! ! IMPLICIT NONE ! ! The aim of this routine is to kill a basin (that is put into another larger one). When ! we do this we need to be careful and change all associated variables. ! INTEGER(i_std) :: tokill, totakeover INTEGER(i_std) :: nbpt, ib ! INTEGER(i_std) :: nwbas INTEGER(i_std), DIMENSION(nbpt) :: basin_count INTEGER(i_std), DIMENSION(nbpt,nwbas) :: basin_id INTEGER(i_std), DIMENSION(nbpt,nwbas) :: basin_flowdir REAL(r_std), DIMENSION(nbpt,nwbas) :: basin_area REAL(r_std), DIMENSION(nbpt,nwbas) :: basin_topoind REAL(r_std), DIMENSION(nbpt,nwbas) :: fetch_basin INTEGER(i_std), DIMENSION(nbpt,nwbas) :: outflow_grid INTEGER(i_std), DIMENSION(nbpt,nwbas) :: outflow_basin INTEGER(i_std), DIMENSION(nbpt,nwbas) :: inflow_number INTEGER(i_std), DIMENSION(nbpt,nwbas,nwbas) :: inflow_basin INTEGER(i_std), DIMENSION(nbpt,nwbas,nwbas) :: inflow_grid ! ! LOCAL ! INTEGER(i_std) :: inf, ibs, ing, inb, ibasf, igrif, it LOGICAL :: doshift ! ! Update the information needed in the basin "totakeover" ! For the moment only area !!$ ! !!$ WRITE(numout,*) 'KILL BASIN :', ib, tokill, totakeover, basin_id(ib,tokill), basin_id(ib,totakeover) !!$ ! ! basin_area(ib, totakeover) = basin_area(ib, totakeover) + basin_area(ib, tokill) basin_topoind(ib, totakeover) = (basin_topoind(ib, totakeover) + basin_topoind(ib, tokill))/2.0 ! ! Add the fetch of the basin will kill to the one which gets the water ! fetch_basin(ib, totakeover) = fetch_basin(ib, totakeover) + fetch_basin(ib, tokill) igrif = outflow_grid(ib,totakeover) ibasf = outflow_basin(ib,totakeover) ! inf = 0 DO WHILE (igrif .GT. 0) fetch_basin(igrif,ibasf) = fetch_basin(igrif,ibasf) + fetch_basin(ib, tokill) it = outflow_grid(igrif, ibasf) ibasf = outflow_basin(igrif, ibasf) igrif = it inf = inf + 1 ENDDO ! ! Take out the basin we have just rerouted from the fetch of the basins in which it used to flow. ! igrif = outflow_grid(ib,tokill) ibasf = outflow_basin(ib,tokill) ! DO WHILE (igrif .GT. 0) fetch_basin(igrif,ibasf) = fetch_basin(igrif,ibasf) - fetch_basin(ib, tokill) it = outflow_grid(igrif, ibasf) ibasf = outflow_basin(igrif, ibasf) igrif = it ENDDO ! ! Redirect the flows which went into the basin to be killed before we change everything ! DO inf = 1, inflow_number(ib, tokill) outflow_basin(inflow_grid(ib, tokill, inf), inflow_basin(ib, tokill, inf)) = totakeover inflow_number(ib, totakeover) = inflow_number(ib, totakeover) + 1 inflow_grid(ib, totakeover, inflow_number(ib, totakeover)) = inflow_grid(ib, tokill, inf) inflow_basin(ib, totakeover, inflow_number(ib, totakeover)) = inflow_basin(ib, tokill, inf) ENDDO ! ! Take out the basin to be killed from the list of inflow basins of the downstream basin ! (In case the basin does not flow into an ocean or lake) ! IF ( outflow_grid(ib,tokill) .GT. 0) THEN ! ing = outflow_grid(ib, tokill) inb = outflow_basin(ib, tokill) doshift = .FALSE. ! DO inf = 1, inflow_number(ing, inb) IF ( doshift ) THEN inflow_grid(ing, inb, inf-1) = inflow_grid(ing, inb, inf) inflow_basin(ing, inb, inf-1) = inflow_basin(ing, inb, inf) ENDIF IF ( inflow_grid(ing, inb, inf) .EQ. ib .AND. inflow_basin(ing, inb, inf) .EQ. tokill) THEN doshift = .TRUE. ENDIF ENDDO ! ! This is only to allow for the last check ! inf = inflow_number(ing, inb) IF ( inflow_grid(ing, inb, inf) .EQ. ib .AND. inflow_basin(ing, inb, inf) .EQ. tokill) THEN doshift = .TRUE. ENDIF ! IF ( .NOT. doshift ) THEN WRITE(numout,*) 'Strange we did not find the basin to kill in the downstream basin' STOP 'routing_killbas' ENDIF inflow_number(ing, inb) = inflow_number(ing, inb) - 1 ! ENDIF ! ! Now remove from the arrays the information of basin "tokill" ! basin_id(ib, tokill:basin_count(ib)-1) = basin_id(ib, tokill+1:basin_count(ib)) basin_flowdir(ib, tokill:basin_count(ib)-1) = basin_flowdir(ib, tokill+1:basin_count(ib)) basin_area(ib, tokill:basin_count(ib)-1) = basin_area(ib, tokill+1:basin_count(ib)) basin_area(ib, basin_count(ib):nwbas) = 0.0 basin_topoind(ib, tokill:basin_count(ib)-1) = basin_topoind(ib, tokill+1:basin_count(ib)) basin_topoind(ib, basin_count(ib):nwbas) = 0.0 fetch_basin(ib, tokill:basin_count(ib)-1) = fetch_basin(ib, tokill+1:basin_count(ib)) fetch_basin(ib, basin_count(ib):nwbas) = 0.0 ! ! Before we remove the information from the outflow fields we have to correct the corresponding inflow fields ! of the grids into which the flow goes ! DO ibs = tokill+1,basin_count(ib) ing = outflow_grid(ib, ibs) inb = outflow_basin(ib, ibs) IF ( ing .GT. 0 ) THEN DO inf = 1, inflow_number(ing, inb) IF ( inflow_grid(ing,inb,inf) .EQ. ib .AND. inflow_basin(ing,inb,inf) .EQ. ibs) THEN inflow_basin(ing,inb,inf) = ibs - 1 ENDIF ENDDO ENDIF ENDDO outflow_grid(ib, tokill:basin_count(ib)-1) = outflow_grid(ib, tokill+1:basin_count(ib)) outflow_basin(ib, tokill:basin_count(ib)-1) = outflow_basin(ib, tokill+1:basin_count(ib)) ! ! Basins which moved down also need to redirect their incoming flows. ! DO ibs=tokill+1, basin_count(ib) DO inf = 1, inflow_number(ib, ibs) outflow_basin(inflow_grid(ib, ibs, inf), inflow_basin(ib, ibs, inf)) = ibs-1 ENDDO ENDDO ! ! Shift the inflow basins ! DO it = tokill+1,basin_count(ib) inflow_grid(ib, it-1, 1:inflow_number(ib,it)) = inflow_grid(ib, it, 1:inflow_number(ib,it)) inflow_basin(ib, it-1, 1:inflow_number(ib,it)) = inflow_basin(ib, it, 1:inflow_number(ib,it)) inflow_number(ib,it-1) = inflow_number(ib,it) ENDDO ! basin_count(ib) = basin_count(ib) - 1 ! END SUBROUTINE routing_killbas ! !--------------------------------------------------------------------------------- ! SUBROUTINE routing_names(numlar, basin_names) ! IMPLICIT NONE ! ! Arguments ! INTEGER(i_std), INTENT(in) :: numlar CHARACTER(LEN=*), INTENT(inout) :: basin_names(numlar) ! ! Local ! INTEGER(i_std), PARAMETER :: listleng=349 INTEGER(i_std) :: lenstr, i CHARACTER(LEN=60), DIMENSION(listleng) :: list_names CHARACTER(LEN=60) :: tmp_str ! lenstr = LEN(basin_names(1)) ! list_names(1) = "Amazon" list_names(2) = "Nile" list_names(3) = "Zaire" list_names(4) = "Mississippi" list_names(5) = "Amur" list_names(6) = "Parana" list_names(7) = "Yenisei" list_names(8) = "Ob" list_names(9) = "Lena" list_names(10) = "Niger" list_names(11) = "Zambezi" list_names(12) = "Yangtze" list_names(13) = "Chang Jiang" list_names(14) = "Mackenzie" list_names(15) = "Ganges" list_names(16) = "Chari" list_names(17) = "Volga" list_names(18) = "St. Lawrence" list_names(19) = "Indus" list_names(20) = "Syr-Darya" list_names(21) = "Nelson" list_names(22) = "Orinoco" list_names(23) = "Murray" list_names(24) = "Great Artesian Basin" list_names(25) = "Shatt el Arab" list_names(26) = "Orange" list_names(27) = "Huang He" list_names(28) = "Yukon" list_names(29) = "Senegal" list_names(30) = "Chott Jerid" list_names(31) = "Jubba" list_names(32) = "Colorado (Ari)" list_names(33) = "Rio Grande (US)" list_names(34) = "Danube" list_names(35) = "Mekong" list_names(36) = "Tocantins" list_names(37) = "Wadi al Farigh" list_names(38) = "Tarim" list_names(39) = "Columbia" list_names(40) = "Noname (GHAASBasin49)" list_names(41) = "Kolyma" list_names(42) = "Sao Francisco" list_names(43) = "Amu-Darya" list_names(44) = "Noname (GHAASBasin51)" list_names(45) = "Dnepr" list_names(46) = "Noname (GHAASBasin61)" list_names(47) = "Don" list_names(48) = "Colorado (Arg)" list_names(49) = "Limpopo" list_names(50) = "GHAASBasin50" list_names(51) = "Zhujiang" list_names(52) = "Irrawaddy" list_names(53) = "Volta" list_names(54) = "GHAASBasin54" list_names(55) = "Farah" list_names(56) = "Khatanga" list_names(57) = "Dvina" list_names(58) = "Urugay" list_names(59) = "Qarqan" list_names(60) = "Noname (GHAASBasin75)" list_names(61) = "Parnaiba" list_names(62) = "Noname (GHAASBasin73)" list_names(63) = "Indigirka" list_names(64) = "Churchill (Hud)" list_names(65) = "Godavari" list_names(66) = "Pur - Taz" list_names(67) = "Pechora" list_names(68) = "Baker" list_names(69) = "Ural" list_names(70) = "Neva" list_names(71) = "Liao" list_names(72) = "Salween" list_names(73) = "GHAASBasin73" list_names(74) = "Jordan" list_names(75) = "Noname (GHAASBasin78)" list_names(76) = "Magdalena" list_names(77) = "Krishna" list_names(78) = "Salado" list_names(79) = "Fraser" list_names(80) = "Hai Ho" list_names(81) = "Huai" list_names(82) = "Yana" list_names(83) = "Noname (GHAASBasin95)" list_names(84) = "Noname (GHAASBasin105)" list_names(85) = "Kura" list_names(86) = "Olenek" list_names(87) = "Ogooue" list_names(88) = "Taymyr" list_names(89) = "Negro Arg" list_names(90) = "Chubut" list_names(91) = "GHAASBasin91" list_names(92) = "Noname (GHAASBasin122)" list_names(93) = "Noname (GHAASBasin120)" list_names(94) = "Sacramento" list_names(95) = "Fitzroy West" list_names(96) = "Grande de Santiago" list_names(97) = "Rufiji" list_names(98) = "Wisla" list_names(99) = "Noname (GHAASBasin47)" list_names(100) = "Noname (GHAASBasin127)" list_names(101) = "Hong" list_names(102) = "Noname (GHAASBasin97)" list_names(103) = "Swan-Avon" list_names(104) = "Rhine" list_names(105) = "Cuanza" list_names(106) = "GHAASBasin106" list_names(107) = "Noname (GHAASBasin142)" list_names(108) = "Roviuna" list_names(109) = "Essequibo" list_names(110) = "Elbe" list_names(111) = "Koksoak" list_names(112) = "Chao Phraya" list_names(113) = "Brahmani" list_names(114) = "Noname (GHAASBasin165)" list_names(115) = "Pyasina" list_names(116) = "Fitzroy East" list_names(117) = "Noname (GHAASBasin173)" list_names(118) = "Albany" list_names(119) = "Sanaga" list_names(120) = "GHAASBasin120" list_names(121) = "Noname (GHAASBasin178)" list_names(122) = "Noname (GHAASBasin148)" list_names(123) = "Brazos (Tex)" list_names(124) = "GHAASBasin124" list_names(125) = "Alabama" list_names(126) = "Noname (GHAASBasin174)" list_names(127) = "Noname (GHAASBasin179)" list_names(128) = "Balsas" list_names(129) = "Noname (GHAASBasin172)" list_names(130) = "Burdekin" list_names(131) = "Colorado (Texas)" list_names(132) = "Noname (GHAASBasin150)" list_names(133) = "Odra" list_names(134) = "Loire" list_names(135) = "Noname (GHAASBasin98)" list_names(136) = "Galana" list_names(137) = "Kuskowin" list_names(138) = "Moose" list_names(139) = "Narmada" list_names(140) = "GHAASBasin140" list_names(141) = "GHAASBasin141" list_names(142) = "Flinders" list_names(143) = "Kizil Irmak" list_names(144) = "GHAASBasin144" list_names(145) = "Save" list_names(146) = "Roper" list_names(147) = "Churchill (Atlantic)" list_names(148) = "GHAASBasin148" list_names(149) = "Victoria" list_names(150) = "Back" list_names(151) = "Bandama" list_names(152) = "Severn (Can)" list_names(153) = "Po" list_names(154) = "GHAASBasin154" list_names(155) = "GHAASBasin155" list_names(156) = "GHAASBasin156" list_names(157) = "Rhone" list_names(158) = "Tana (Ken)" list_names(159) = "La Grande" list_names(160) = "GHAASBasin160" list_names(161) = "Cunene" list_names(162) = "Douro" list_names(163) = "GHAASBasin163" list_names(164) = "Nemanus" list_names(165) = "GHAASBasin165" list_names(166) = "Anabar" list_names(167) = "Hayes" list_names(168) = "Mearim" list_names(169) = "GHAASBasin169" list_names(170) = "Panuco" list_names(171) = "GHAASBasin171" list_names(172) = "Doce" list_names(173) = "Gasgoyne" list_names(174) = "GHAASBasin174" list_names(175) = "GHAASBasin175" list_names(176) = "Ashburton" list_names(177) = "GHAASBasin177" list_names(178) = "Peel" list_names(179) = "Daugava" list_names(180) = "GHAASBasin180" list_names(181) = "Ebro" list_names(182) = "Comoe" list_names(183) = "Jacui" list_names(184) = "GHAASBasin184" list_names(185) = "Kapuas" list_names(186) = "GHAASBasin186" list_names(187) = "Penzhina" list_names(188) = "Cauweri" list_names(189) = "GHAASBasin189" list_names(190) = "Mamberamo" list_names(191) = "Sepik" list_names(192) = "GHAASBasin192" list_names(193) = "Sassandra" list_names(194) = "GHAASBasin194" list_names(195) = "GHAASBasin195" list_names(196) = "Nottaway" list_names(197) = "Barito" list_names(198) = "GHAASBasin198" list_names(199) = "Seine" list_names(200) = "Tejo" list_names(201) = "GHAASBasin201" list_names(202) = "Gambia" list_names(203) = "Susquehanna" list_names(204) = "Dnestr" list_names(205) = "Murchinson" list_names(206) = "Deseado" list_names(207) = "Mitchell" list_names(208) = "Mahakam" list_names(209) = "GHAASBasin209" list_names(210) = "Pangani" list_names(211) = "GHAASBasin211" list_names(212) = "GHAASBasin212" list_names(213) = "GHAASBasin213" list_names(214) = "GHAASBasin214" list_names(215) = "GHAASBasin215" list_names(216) = "Bug" list_names(217) = "GHAASBasin217" list_names(218) = "Usumacinta" list_names(219) = "Jequitinhonha" list_names(220) = "GHAASBasin220" list_names(221) = "Corantijn" list_names(222) = "Fuchun Jiang" list_names(223) = "Copper" list_names(224) = "Tapti" list_names(225) = "Menjiang" list_names(226) = "Karun" list_names(227) = "Mezen" list_names(228) = "Guadiana" list_names(229) = "Maroni" list_names(230) = "GHAASBasin230" list_names(231) = "Uda" list_names(232) = "GHAASBasin232" list_names(233) = "Kuban" list_names(234) = "Colville" list_names(235) = "Thaane" list_names(236) = "Alazeya" list_names(237) = "Paraiba do Sul" list_names(238) = "GHAASBasin238" list_names(239) = "Fortesque" list_names(240) = "GHAASBasin240" list_names(241) = "GHAASBasin241" list_names(242) = "Winisk" list_names(243) = "GHAASBasin243" list_names(244) = "GHAASBasin244" list_names(245) = "Ikopa" list_names(246) = "Gilbert" list_names(247) = "Kouilou" list_names(248) = "Fly" list_names(249) = "GHAASBasin249" list_names(250) = "GHAASBasin250" list_names(251) = "GHAASBasin251" list_names(252) = "Mangoky" list_names(253) = "Damodar" list_names(254) = "Onega" list_names(255) = "Moulouya" list_names(256) = "GHAASBasin256" list_names(257) = "Ord" list_names(258) = "GHAASBasin258" list_names(259) = "GHAASBasin259" list_names(260) = "GHAASBasin260" list_names(261) = "GHAASBasin261" list_names(262) = "Narva" list_names(263) = "GHAASBasin263" list_names(264) = "Seal" list_names(265) = "Cheliff" list_names(266) = "Garonne" list_names(267) = "Rupert" list_names(268) = "GHAASBasin268" list_names(269) = "Brahmani" list_names(270) = "Sakarya" list_names(271) = "Gourits" list_names(272) = "Sittang" list_names(273) = "Rajang" list_names(274) = "Evros" list_names(275) = "Appalachicola" list_names(276) = "Attawapiskat" list_names(277) = "Lurio" list_names(278) = "Daly" list_names(279) = "Penner" list_names(280) = "GHAASBasin280" list_names(281) = "GHAASBasin281" list_names(282) = "Guadalquivir" list_names(283) = "Nadym" list_names(284) = "GHAASBasin284" list_names(285) = "Saint John" list_names(286) = "GHAASBasin286" list_names(287) = "Cross" list_names(288) = "Omoloy" list_names(289) = "Oueme" list_names(290) = "GHAASBasin290" list_names(291) = "Gota" list_names(292) = "Nueces" list_names(293) = "Stikine" list_names(294) = "Yalu" list_names(295) = "Arnaud" list_names(296) = "GHAASBasin296" list_names(297) = "Jequitinhonha" list_names(298) = "Kamchatka" list_names(299) = "GHAASBasin299" list_names(300) = "Grijalva" list_names(301) = "GHAASBasin301" list_names(302) = "Kemijoki" list_names(303) = "Olifants" list_names(304) = "GHAASBasin304" list_names(305) = "Tsiribihina" list_names(306) = "Coppermine" list_names(307) = "GHAASBasin307" list_names(308) = "GHAASBasin308" list_names(309) = "Kovda" list_names(310) = "Trinity" list_names(311) = "Glama" list_names(312) = "GHAASBasin312" list_names(313) = "Luan" list_names(314) = "Leichhardt" list_names(315) = "GHAASBasin315" list_names(316) = "Gurupi" list_names(317) = "GR Baleine" list_names(318) = "Aux Feuilles" list_names(319) = "GHAASBasin319" list_names(320) = "Weser" list_names(321) = "GHAASBasin321" list_names(322) = "GHAASBasin322" list_names(323) = "Yesil" list_names(324) = "Incomati" list_names(325) = "GHAASBasin325" list_names(326) = "GHAASBasin326" list_names(327) = "Pungoe" list_names(328) = "GHAASBasin328" list_names(329) = "Meuse" list_names(330) = "Eastmain" list_names(331) = "Araguari" list_names(332) = "Hudson" list_names(333) = "GHAASBasin333" list_names(334) = "GHAASBasin334" list_names(335) = "GHAASBasin335" list_names(336) = "GHAASBasin336" list_names(337) = "Kobuk" list_names(338) = "Altamaha" list_names(339) = "GHAASBasin339" list_names(340) = "Mand" list_names(341) = "Santee" list_names(342) = "GHAASBasin342" list_names(343) = "GHAASBasin343" list_names(344) = "GHAASBasin344" list_names(345) = "Hari" list_names(346) = "GHAASBasin346" list_names(347) = "Wami" list_names(348) = "GHAASBasin348" list_names(349) = "GHAASBasin349" ! basin_names(:) = ' ' ! DO i=1,numlar tmp_str = list_names(i) basin_names(i) = tmp_str(1:MIN(lenstr,LEN_TRIM(tmp_str))) ENDDO ! END SUBROUTINE routing_names ! !--------------------------------------------------------------------------------- ! SUBROUTINE routing_irrigmap (nbpt, index, lalo, neighbours, resolution, contfrac, & & init_irrig, irrigated, init_flood, floodplains, hist_id, hist2_id) ! ! This program will interpoalte the 0.5 x 05 deg based data set to the resolution ! of the model. ! IMPLICIT NONE ! ! INPUT ! INTEGER(i_std), INTENT(in) :: nbpt ! Number of points for which the data needs to be interpolated 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,8) ! Vector of neighbours for each grid point (1=N, 2=E, 3=S, 4=W) REAL(r_std), INTENT(in) :: resolution(nbpt,2) ! The size in m of each grid-box in X and Y REAL(r_std), INTENT(in) :: contfrac(nbpt) ! Fraction of land in each grid box LOGICAL, INTENT(in) :: init_irrig ! Do we need to compute irrigation ? REAL(r_std), INTENT(out) :: irrigated(:) ! Irrigated surface in each grid-box LOGICAL, INTENT(in) :: init_flood ! Do we need to compute floodplains REAL(r_std), INTENT(out) :: floodplains(:) ! Surface which can be inondated in each grid-box INTEGER(i_std), INTENT(in) :: hist_id ! Access to history file INTEGER(i_std), INTENT(in) :: hist2_id ! Access to history file 2 ! ! LOCAL ! REAL(r_std), PARAMETER :: R_Earth = 6378000. ! CHARACTER(LEN=80) :: filename INTEGER(i_std) :: iml, jml, lml, tml, fid, ib, ip, jp, ilf, lastjp, nbexp REAL(r_std) :: lev(1), date, dt, coslat, pi REAL(r_std) :: lon_up, lon_low, lat_up, lat_low INTEGER(i_std) :: itau(1) REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: lat_rel, lon_rel, lat_ful, lon_ful REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: irrigated_frac, flood_frac REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: loup_rel, lolow_rel, laup_rel, lalow_rel REAL(r_std) :: area_irrig, area_flood, ax, ay, sgn, surp REAL(r_std) :: lonrel, louprel, lolowrel REAL(r_std) :: irrigmap(nbpt) ! pi = 4. * ATAN(1.) ! !Config Key = IRRIGATION_FILE !Config Desc = Name of file which contains the map of irrigated areas !Config Def = irrigated.nc !Config If = IRRIGATE !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). ! filename = 'irrigated.nc' CALL getin_p('IRRIGATION_FILE',filename) ! IF (is_root_prc) CALL flininfo(filename,iml, jml, lml, tml, fid) CALL bcast(iml) CALL bcast(jml) CALL bcast(lml) CALL bcast(tml) ! ! ALLOCATE (lat_rel(iml,jml)) ALLOCATE (lon_rel(iml,jml)) ALLOCATE (laup_rel(iml,jml)) ALLOCATE (loup_rel(iml,jml)) ALLOCATE (lalow_rel(iml,jml)) ALLOCATE (lolow_rel(iml,jml)) ALLOCATE (lat_ful(iml+2,jml+2)) ALLOCATE (lon_ful(iml+2,jml+2)) ALLOCATE (irrigated_frac(iml,jml)) ALLOCATE (flood_frac(iml,jml)) ! IF (is_root_prc) CALL flinopen(filename, .FALSE., iml, jml, lml, lon_rel, lat_rel, lev, tml, itau, date, dt, fid) CALL bcast(lon_rel) CALL bcast(lat_rel) CALL bcast(lev) CALL bcast(itau) CALL bcast(date) CALL bcast(dt) ! IF (is_root_prc) CALL flinget(fid, 'irrig', iml, jml, lml, tml, 1, 1, irrigated_frac) CALL bcast(irrigated_frac) IF (is_root_prc) CALL flinget(fid, 'inundated', iml, jml, lml, tml, 1, 1, flood_frac) CALL bcast(flood_frac) ! 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-1.) THEN irrigated_frac(ip,jp) = irrigated_frac(ip,jp)/100. IF ( irrigated_frac(ip,jp) < 0.005 ) irrigated_frac(ip,jp) = 0.0 ENDIF ! IF ( flood_frac(ip,jp) .LT. undef_sechiba-1.) THEN flood_frac(ip,jp) = flood_frac(ip,jp)/100 IF ( flood_frac(ip,jp) < 0.005 ) flood_frac(ip,jp) = 0.0 ENDIF ! ENDDO ENDDO ! WRITE(numout,*) 'lon_rel : ', MAXVAL(lon_rel), MINVAL(lon_rel) WRITE(numout,*) 'lat_rel : ', MAXVAL(lat_rel), MINVAL(lat_rel) WRITE(numout,*) 'irrigated_frac : ', MINVAL(irrigated_frac, MASK=irrigated_frac .GT. 0), & & MAXVAL(irrigated_frac, MASK=irrigated_frac .LT. undef_sechiba) WRITE(numout,*) 'flood_frac : ', MINVAL(flood_frac, MASK=flood_frac .GT. 0), & & MAXVAL(flood_frac, MASK=flood_frac .LT. undef_sechiba) ! nbexp = 0 ! ! Duplicate the border assuming we have a global grid going from west to east ! lon_ful(2:iml+1,2:jml+1) = lon_rel(1:iml,1:jml) lat_ful(2:iml+1,2:jml+1) = lat_rel(1:iml,1:jml) ! IF ( lon_rel(iml,1) .LT. lon_ful(2,2)) THEN lon_ful(1,2:jml+1) = lon_rel(iml,1:jml) lat_ful(1,2:jml+1) = lat_rel(iml,1:jml) ELSE lon_ful(1,2:jml+1) = lon_rel(iml,1:jml)-360 lat_ful(1,2:jml+1) = lat_rel(iml,1:jml) ENDIF IF ( lon_rel(1,1) .GT. lon_ful(iml+1,2)) THEN lon_ful(iml+2,2:jml+1) = lon_rel(1,1:jml) lat_ful(iml+2,2:jml+1) = lat_rel(1,1:jml) ELSE lon_ful(iml+2,2:jml+1) = lon_rel(1,1:jml)+360 lat_ful(iml+2,2:jml+1) = lat_rel(1,1:jml) ENDIF ! sgn = lat_rel(1,1)/ABS(lat_rel(1,1)) lat_ful(2:iml+1,1) = sgn*180 - lat_rel(1:iml,1) sgn = lat_rel(1,jml)/ABS(lat_rel(1,jml)) lat_ful(2:iml+1,jml+2) = sgn*180 - lat_rel(1:iml,jml) lat_ful(1,1) = lat_ful(iml+1,1) lat_ful(iml+2,1) = lat_ful(2,1) lat_ful(1,jml+2) = lat_ful(iml+1,jml+2) lat_ful(iml+2,jml+2) = lat_ful(2,jml+2) ! ! Add the longitude lines to the top and bottom ! lon_ful(:,1) = lon_ful(:,2) lon_ful(:,jml+2) = lon_ful(:,jml+1) ! ! Get the upper and lower limits of each grid box ! DO ip=1,iml DO jp=1,jml loup_rel(ip,jp) =MAX(0.5*(lon_ful(ip,jp+1)+lon_ful(ip+1,jp+1)), 0.5*(lon_ful(ip+1,jp+1)+lon_ful(ip+2,jp+1))) lolow_rel(ip,jp) =MIN(0.5*(lon_ful(ip,jp+1)+lon_ful(ip+1,jp+1)), 0.5*(lon_ful(ip+1,jp+1)+lon_ful(ip+2,jp+1))) laup_rel(ip,jp) =MAX(0.5*(lat_ful(ip+1,jp)+lat_ful(ip+1,jp+1)), 0.5*(lat_ful(ip+1,jp+1)+lat_ful(ip+1,jp+2))) lalow_rel(ip,jp) =MIN(0.5*(lat_ful(ip+1,jp)+lat_ful(ip+1,jp+1)), 0.5*(lat_ful(ip+1,jp+1)+lat_ful(ip+1,jp+2))) ENDDO ENDDO ! ! Now we take each grid point and find out which values from the forcing we need to average ! DO ib =1, nbpt ! ! We find the 4 limits of the grid-box. As we transform the resolution of the model ! into longitudes and latitudes we do not have the problem of periodicity. ! coslat is a help variable here ! ! coslat = MAX(COS(lalo(ib,1) * pi/180. ), 0.001 )*pi/180. * R_Earth ! lon_up = lalo(ib,2)+ resolution(ib,1)/(2.0*coslat) lon_low =lalo(ib,2) - resolution(ib,1)/(2.0*coslat) ! coslat = pi/180. * R_Earth ! lat_up =lalo(ib,1)+resolution(ib,2)/(2.0*coslat) lat_low =lalo(ib,1)-resolution(ib,2)/(2.0*coslat) ! ! ! Find the grid boxes from the data that go into the model's boxes ! We still work as if we had a regular grid ! Well it needs to be localy regular so ! so that the longitude at the latitude of the last found point is close to the one of the next point. ! lastjp = 1 area_irrig = zero area_flood = zero DO ip=1,iml ! ! Either the center of the data grid point is in the interval of the model grid or ! the East and West limits of the data grid point are on either sides of the border of ! the data grid. ! ! ! We find the 4 limits of the grid-box. As we transform the resolution of the model ! into longitudes and latitudes we do not have the problem of periodicity. ! coslat is a help variable here ! ! ! ! To do that correctly we have to check if the grid box sits on the date-line. ! IF ( lon_low < -180.0 ) THEN lonrel = MOD( lon_rel(ip,lastjp) - 360.0, 360.0) lolowrel = MOD( lolow_rel(ip,lastjp) - 360.0, 360.0) louprel = MOD( loup_rel(ip,lastjp) - 360.0, 360.0) ! ELSE IF ( lon_up > 180.0 ) THEN lonrel = MOD( 360. - lon_rel(ip,lastjp), 360.0) lolowrel = MOD( 360. - lolow_rel(ip,lastjp), 360.0) louprel = MOD( 360. - loup_rel(ip,lastjp), 360.0) ELSE lonrel = lon_rel(ip,lastjp) lolowrel = lolow_rel(ip,lastjp) louprel = loup_rel(ip,lastjp) ENDIF ! ! ! IF ( lonrel > lon_low .AND. lonrel < lon_up .OR. & & lolowrel < lon_low .AND. louprel > lon_low .OR. & & lolowrel < lon_up .AND. louprel > lon_up ) THEN ! DO jp = 1, jml ! ! Now that we have the longitude let us find the latitude ! IF ( lat_rel(ip,jp) > lat_low .AND. lat_rel(ip,jp) < lat_up .OR. & & lalow_rel(ip,jp) < lat_low .AND. laup_rel(ip,jp) > lat_low .OR.& & lalow_rel(ip,jp) < lat_up .AND. laup_rel(ip,jp) > lat_up) THEN ! lastjp = jp ! ! Mising values in the file are assumed to be 1e20 ! IF ( lon_low < -180.0 ) THEN lolowrel = MOD( lolow_rel(ip,jp) - 360.0, 360.0) louprel = MOD( loup_rel(ip,jp) - 360.0, 360.0) ! ELSE IF ( lon_up > 180.0 ) THEN lolowrel = MOD( 360. - lolow_rel(ip,jp), 360.0) louprel = MOD( 360. - loup_rel(ip,jp), 360.0) ELSE lolowrel = lolow_rel(ip,jp) louprel = loup_rel(ip,jp) ENDIF ! ! Get the area of the fine grid in the model grid ! coslat = MAX( COS( lat_rel(ip,jp) * pi/180. ), 0.001 ) ax = (MIN(lon_up,louprel)-MAX(lon_low, lolowrel))*pi/180. * R_Earth * coslat ay = (MIN(lat_up, laup_rel(ip,jp))-MAX(lat_low,lalow_rel(ip,jp)))*pi/180. * R_Earth ! IF (irrigated_frac(ip,jp) .LT. undef_sechiba-1.) THEN area_irrig = area_irrig + ax*ay*irrigated_frac(ip,jp) ENDIF ! IF (flood_frac(ip,jp) .LT. undef_sechiba-1.) THEN area_flood = area_flood + ax*ay*flood_frac(ip,jp) ENDIF ! ENDIF ! ENDDO ! ENDIF ! ENDDO ! ! Put the toal 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 STOP ENDIF ! Compute a diagnostic of the map. IF(contfrac(ib).GT.0.0) THEN irrigmap (ib) = irrigated(ib) / ( resolution(ib,1)*resolution(ib,2)*contfrac(ib) ) ELSE irrigmap (ib) = zero ENDIF ! ENDIF ! IF ( init_flood ) THEN floodplains(ib) = MIN(area_flood, resolution(ib,1)*resolution(ib,2)*contfrac(ib)) IF ( floodplains(ib) < 0 ) THEN WRITE(numout,*) 'We have a problem here : ', floodplains(ib) WRITE(numout,*) 'resolution :', resolution(ib,1), resolution(ib,2) WRITE(numout,*) area_flood STOP ENDIF ENDIF ! ! ENDDO ! IF ( init_irrig .AND. init_flood ) THEN DO ib = 1, nbpt IF ( floodplains(ib)+irrigated(ib) > resolution(ib,1)*resolution(ib,2)*contfrac(ib) ) THEN surp = (resolution(ib,1)*resolution(ib,2)*contfrac(ib))/(floodplains(ib)+irrigated(ib)) floodplains(ib) = floodplains(ib) * surp irrigated(ib) = irrigated(ib) * surp ENDIF ENDDO ENDIF ! IF ( init_irrig ) THEN WRITE(numout,*) 'RESULT irrigated : ', MINVAL(irrigated), MAXVAL(irrigated) IF ( .NOT. almaoutput ) THEN CALL histwrite(hist_id, 'irrigmap', 1, irrigmap, nbpt, index) ELSE CALL histwrite(hist_id, 'IrrigationMap', 1, irrigated, nbpt, index) ENDIF IF ( hist2_id > 0 ) THEN IF ( .NOT. almaoutput ) THEN CALL histwrite(hist2_id, 'irrigmap', 1, irrigmap, nbpt, index) ELSE CALL histwrite(hist2_id, 'IrrigationMap', 1, irrigated, nbpt, index) ENDIF ENDIF ENDIF ! IF ( init_flood ) THEN WRITE(numout,*) 'RESULT floodplains : ', MINVAL(floodplains), MAXVAL(floodplains) ENDIF ! END SUBROUTINE routing_irrigmap ! !--------------------------------------------------------------------------------- ! SUBROUTINE routing_waterbal(nbpt, firstcall, runoff, drainage, returnflow, & & irrigation, riverflow, coastalflow) ! IMPLICIT NONE ! ! This subroutine should allow us to check the water balance in the routing module. ! INTEGER(i_std), INTENT(in) :: nbpt ! Domain size LOGICAL, INTENT(in) :: firstcall ! controls behaviour REAL(r_std), INTENT(in) :: runoff(nbpt) ! grid-point runoff REAL(r_std), INTENT(in) :: drainage(nbpt) ! grid-point drainage REAL(r_std), INTENT(in) :: returnflow(nbpt) ! The water flow which returns to the grid box (kg/m^2/dt) REAL(r_std), INTENT(in) :: irrigation(nbpt) ! Irrigation flux (kg/m^2 per dt) REAL(r_std), INTENT(in) :: riverflow(nbpt) ! Outflow of the major rivers (kg/dt) REAL(r_std), INTENT(in) :: coastalflow(nbpt) ! Outflow on coastal points by small basins (kg/dt) ! ! We sum-up all the water we have in the warious reservoirs ! REAL(r_std), SAVE :: totw_stream, totw_fast, totw_slow, totw_lake REAL(r_std), SAVE :: totw_in, totw_out, totw_return, totw_irrig, totw_river, totw_coastal REAL(r_std) :: totarea, area ! ! Just to make sure we do not get too large numbers ! ! REAL(r_std), PARAMETER :: scaling = 1.0E+6 REAL(r_std), PARAMETER :: allowed_err = 50. ! INTEGER(i_std) :: ig ! IF ( firstcall ) THEN ! totw_stream = zero totw_fast = zero totw_slow = zero totw_lake = zero totw_in = zero ! DO ig=1,nbpt ! totarea = SUM(routing_area(ig,:)) ! 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_in = totw_in + (runoff(ig)*totarea + drainage(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_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_return = totw_return + 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_stream + totw_fast + totw_slow + totw_lake) - (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 : stream fast ', totw_stream, totw_fast WRITE(numout,*) '--Water-- change : slow, lake ', totw_slow, totw_lake WRITE(numout,*) '--Water>>> change in the routing res. : ', 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_stream + totw_fast + totw_slow + totw_lake) - (totw_out - totw_in))/area ENDIF ! ENDIF ! END SUBROUTINE routing_waterbal ! ! END MODULE routing