MODULE routing_reg USE defprec USE grid USE constantes_var USE routing_tools IMPLICIT NONE INTEGER(i_std), SAVE :: nbvmax=440 !! The maximum number of basins we can handle at any time during the generation of the maps (unitless) INTEGER(i_std), SAVE :: nbxmax=63 !! The maximum number of points in one direction (NS or WE) using in routing_reg_getgrid (unitless) INTEGER(i_std), SAVE :: maxpercent=2 !! The maximum area percentage of a sub-basin in a grib-box (default = 2%) INTEGER(i_std), SAVE :: nbpt_save INTEGER(i_std), SAVE :: nbasmax_save REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: routing_area_glo !! Surface of basin (m^2) REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: routing_orog_glo !! Mean topography (m) REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: routing_floodp_glo !! Surface of floodplains (m^2) REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: routing_cg_glo !! Centre of gravity of HTU (Lat, Lon) REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: topo_resid_glo !! Topographic index of the retention time (m) INTEGER(i_std), SAVE, ALLOCATABLE, DIMENSION(:) :: route_count_glo !! Number of basins per grid INTEGER(i_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: route_togrid_glo !! Grid into which the basin flows (unitless) INTEGER(i_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: route_tobasin_glo !! Basin in to which the water goes (unitless) INTEGER(i_std), SAVE, ALLOCATABLE, DIMENSION(:) :: route_nbintobas_glo !! Number of basin into current one (unitless) INTEGER(i_std), SAVE, ALLOCATABLE, DIMENSION(:) :: origin_nbintobas_glo !! Number of sub-grid basin into current one before truncation (unitless) INTEGER(i_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: global_basinid_glo !! ID of basin (unitless) REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: route_outlet_glo !! Coordinate of outlet (-) REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: route_type_glo !! Coordinate of outlet (-) REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: route_fetch_glo !! Upstream area at each HTU CONTAINS !! ================================================================================================================================ !! SUBROUTINE : routing_reg_getgrid !! !>\BRIEF This subroutine extracts from the global high resolution fields !! the data for the current grid box we are dealing with. !! !! DESCRIPTION (definitions, functional, design, flags) : !! 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).:\n !! !! RECENT CHANGE(S): None !! !! MAIN OUTPUT VARIABLE(S): !! !! REFERENCES : None !! !! FLOWCHART : None !! \n !_ ================================================================================================================================ SUBROUTINE routing_reg_getgrid(nbpt, ib, sub_pts, sub_index, sub_area, max_basins, min_topoind, & & lon_rel, lat_rel, lalo, resolution, contfrac, trip, basins, topoindex, fac, hierarchy, & & orog, floodp,& & nbi, nbj, area_bx, trip_bx, basin_bx, topoind_bx, fac_bx, hierarchy_bx,orog_bx, floodp_bx, lon_bx, lat_bx, lshead_bx) ! IMPLICIT NONE ! !! INPUT VARIABLES INTEGER(i_std), INTENT(in) :: nbpt !! Domain size (unitless) INTEGER(i_std), INTENT(in) :: ib !! Current basin (unitless) INTEGER(i_std), INTENT(in) :: sub_pts(nbpt) !! Number of high resolution points on this grid (unitless) INTEGER(i_std), INTENT(in) :: sub_index(nbpt,nbvmax,2) !! Indices of the points we need on the fine grid (unitless) REAL(r_std), INTENT(inout) :: max_basins !! The current maximum of basins REAL(r_std), INTENT(in) :: min_topoind !! The current minimum of topographic index (m) REAL(r_std), INTENT(in) :: sub_area(nbpt,nbvmax) !! Area on the fine grid (m^2) ! REAL(r_std), INTENT(in) :: lon_rel(nbpt,nbvmax) !! REAL(r_std), INTENT(in) :: lat_rel(nbpt,nbvmax) !! 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 of each grid box in X and Y (m) REAL(r_std), INTENT(in) :: contfrac(nbpt) !! Fraction of land in each grid box (unitless;0-1) ! REAL(r_std), INTENT(inout) :: trip(nbpt,nbvmax) !! The trip field (unitless) REAL(r_std), INTENT(inout) :: basins(nbpt,nbvmax) !! data on the fine grid REAL(r_std), INTENT(inout) :: topoindex(nbpt,nbvmax) !! Topographic index of the residence time (m) REAL(r_std), INTENT(inout) :: hierarchy(nbpt,nbvmax) !! data on the fine grid REAL(r_std), INTENT(inout) :: fac(nbpt,nbvmax) !! data on the fine grid REAL(r_std), INTENT(inout) :: orog(nbpt,nbvmax) !! data on the fine grid REAL(r_std), INTENT(inout) :: floodp(nbpt,nbvmax) !! data on the fine grid ! !! OUTPUT VARIABLES INTEGER(i_std), INTENT(out) :: nbi, nbj !! Number of point in x and y within the grid (unitless) REAL(r_std), INTENT(out) :: area_bx(nbxmax,nbxmax) !! Area of each small box in the grid box (m^2) REAL(r_std), INTENT(out) :: hierarchy_bx(nbxmax,nbxmax) !! Level in the basin of the point REAL(r_std), INTENT(out) :: fac_bx(nbxmax,nbxmax) !! Flow accumulation REAL(r_std), INTENT(out) :: lon_bx(nbxmax,nbxmax) !! REAL(r_std), INTENT(out) :: lat_bx(nbxmax,nbxmax) !! REAL(r_std), INTENT(out) :: lshead_bx(nbxmax,nbxmax) !! Large scale heading for outflow points. REAL(r_std), INTENT(out) :: topoind_bx(nbxmax,nbxmax) !! Topographic index of the residence time for each of the smaller boxes (m) INTEGER(i_std), INTENT(out) :: trip_bx(nbxmax,nbxmax) !! The trip field for each of the smaller boxes (unitless) INTEGER(i_std), INTENT(out) :: basin_bx(nbxmax,nbxmax) !! REAL(i_std), INTENT(out) :: orog_bx(nbxmax,nbxmax) !! REAL(i_std), INTENT(out) :: floodp_bx(nbxmax,nbxmax) !! ! !! LOCAL VARIABLES INTEGER(i_std) :: ip, jp, ll(1), iloc, jloc !! Indices (unitless) INTEGER(i_std) :: ipp, jpp INTEGER(i_std), DIMENSION(8,2) :: inc REAL(r_std) :: cenlon, cenlat, dlon, dlat, deslon, deslat, facti, factj REAL(r_std) :: lonstr(nbxmax*nbxmax) !! REAL(r_std) :: latstr(nbxmax*nbxmax) !! ! !_ ================================================================================================================================ ! ! 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 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 fac_bx(:,:) = undef_sechiba lon_bx(:,:) = undef_sechiba lat_bx(:,:) = undef_sechiba lshead_bx(:,:) = undef_sechiba orog_bx(:,:) = undef_sechiba floodp_bx(:,:) = undef_sechiba cenlon = zero cenlat = zero ! IF ( sub_pts(ib) > 0 ) THEN ! DO ip=1,sub_pts(ib) IF ( ip >nbxmax*nbxmax ) THEN CALL ipslerr_p(3,'routing_reg_getgrid','nbxmax too small when filling lonstr',& & 'Please change method to estimate nbxmax','') ENDIF lonstr(ip) = lon_rel(ib, ip) latstr(ip) = lat_rel(ib, ip) 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) ! ! Verify dimension and allocated space ! IF ( nbi > nbxmax .OR. nbj > nbxmax ) THEN WRITE(numout,*) "size of area : nbi=",nbi,"nbj=",nbj, "nbxmax=", nbxmax CALL ipslerr_p(3,'routing_reg_getgrid','nbxmax too small','Please change method to estimate nbxmax','') ENDIF ! ! 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(ib, ip))) iloc = ll(1) ll = MINLOC(ABS(latstr(1:nbj) - lat_rel(ib, ip))) jloc = ll(1) ! trip_bx(iloc, jloc) = NINT(trip(ib, ip)) basin_bx(iloc, jloc) = NINT(basins(ib, ip)) area_bx(iloc, jloc) = sub_area(ib, ip) topoind_bx(iloc, jloc) = topoindex(ib, ip) hierarchy_bx(iloc, jloc) = hierarchy(ib, ip) fac_bx(iloc, jloc) = fac(ib, ip) orog_bx(iloc, jloc) = orog(ib, ip) floodp_bx(iloc, jloc) = floodp(ib, ip) * sub_area(ib, ip) lon_bx(iloc, jloc) = lon_rel(ib, ip) lat_bx(iloc, jloc) = lat_rel(ib, ip) cenlon = cenlon + lon_rel(ib, ip)/sub_pts(ib) cenlat = cenlat + lat_rel(ib, ip)/sub_pts(ib) 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 ! Be careful here: We deal with basin_hierarchy(ib,ij) = min_topoind later! hierarchy_bx(iloc, jloc) = min_topoind fac_bx(iloc, jloc) = zero 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 DO ip=1,nbi ! ! Compute the destination of the flow on the high resolution grid (if possible !) to see if it still within the points ! belonging to the grid we are working on. ! IF ( trip_bx(ip,jp) < 97 ) THEN ! ipp = ip+inc(trip_bx(ip,jp),1) jpp = jp+inc(trip_bx(ip,jp),2) ! ! Check if the indices are outside of the box, then we have a point on the domain border ! IF ( ipp < 1 .OR. ipp > nbi .OR. jpp < 1 .OR. jpp > nbj ) THEN IF (routing_nextgrid(ib,trip_bx(ip,jp)) < -1 ) THEN trip_bx(ip,jp) = 98 ELSE trip_bx(ip,jp) = trip_bx(ip,jp) + 100 ENDIF ELSE ! ! It can also be a border point if the neighbour is not defined ! IF ( basin_bx(ipp,jpp) > undef_int-1 ) THEN IF (routing_nextgrid(ib,trip_bx(ip,jp)) < -1 ) THEN trip_bx(ip,jp) = 98 ELSE trip_bx(ip,jp) = trip_bx(ip,jp) + 100 ENDIF ENDIF ENDIF ELSE IF ( trip_bx(ip,jp) > 100 .AND. trip_bx(ip,jp) < 109 ) THEN WRITE(numout,*) 'WARNING : Point flows our of routing.nc file ' WRITE(numout,*) 'WARNING : Point : basin = ', basin_bx(ip,jp) WRITE(numout,*) 'WARNING : Point : coord = ', lon_bx(ip, jp), lat_bx(ip, jp) WRITE(numout,*) 'WARNING : Please consider using a larger domaine for the routing.nc file' trip_bx(ip,jp) = 98 ENDIF ! ENDDO ENDDO ! ! Compute the large scale flow direction for each outflow point (trip_bx > 100). This ! is done by walking into the direction indicated by the small scale trip. The distance of this walk ! will depend how far we are from the corner of the polygon. If we are in ! the middle of a segment we walk less far than when we are on either corner. ! dlon=ABS(corners_g(ib,NbSegments,1)-corners_g(ib,1,1)) dlat=ABS(corners_g(ib,NbSegments,2)-corners_g(ib,1,2)) DO ip=1,NbSegments-1 dlon = MAX(dlon,ABS(corners_g(ib,ip+1,1)-corners_g(ib,ip,1))) dlat = MAX(dlat,ABS(corners_g(ib,ip+1,2)-corners_g(ib,ip,2))) ENDDO ! DO jp=1,nbj DO ip=1,nbi ! IF ( trip_bx(ip,jp) < undef_int ) THEN lshead_bx(ip,jp) = zero IF ( trip_bx(ip,jp) > 100 ) THEN IF ( nbi > 1 ) THEN facti = REAL(ABS((nbi-1)/2.0-(ip-1)), r_std)/((nbi-1)/2.0) ELSE facti = un ENDIF deslon = lon_bx(ip,jp) + inc(trip_bx(ip,jp)-100,1)*dlon/2.0*facti IF ( nbj > 1 ) THEN factj = REAL(ABS((nbj-1)/2.0-(jp-1)), r_std)/((nbj-1)/2.0) ELSE factj = un ENDIF deslat = lat_bx(ip,jp) - inc(trip_bx(ip,jp)-100,2)*dlat/2.0*factj lshead_bx(ip,jp) = MOD(haversine_heading(deslon, deslat, cenlon, cenlat)+180.0, 360.0) ENDIF ENDIF ENDDO ENDDO ! END SUBROUTINE routing_reg_getgrid ! !! ================================================================================================================================ !! SUBROUTINE : routing_reg_findbasins !! !>\BRIEF This subroutine finds 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. This is the faster version of the old routing_reg_findbasins !! !! DESCRIPTION (definitions, functional, design, flags) : !! 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 \n !! ------------------- \n !! Inside of the box : \n !! trip = 1 : flow = N \n !! trip = 2 : flow = NE \n !! trip = 3 : flow = E \n !! trip = 4 : flow = SE \n !! trip = 5 : flow = S \n !! trip = 6 : flow = SW \n !! trip = 7 : flow = W \n !! trip = 8 : flow = NW \n !! trip = 97 : return flow into the ground \n !! trip = 98 : coastal flow (diffuse flow into the oceans) These values are created here \n !! trip = 99 : river flow into the oceans \n !! !! Out flow from the grid : \n !! trip = 101 : flow = N out of the coarse grid \n !! trip = 102 : flow = NE out of the coarse grid \n !! trip = 103 : flow = E out of the coarse grid \n !! trip = 104 : flow = SE out of the coarse grid \n !! trip = 105 : flow = S out of the coarse grid \n !! trip = 106 : flow = SW out of the coarse grid \n !! trip = 107 : flow = W out of the coarse grid \n !! trip = 108 : flow = NW out of the coarse grid! \n !! RECENT CHANGE(S): None !! !! MAIN OUTPUT VARIABLE(S): !! !! REFERENCES : None !! !! FLOWCHART : None !! \n !_ ================================================================================================================================ SUBROUTINE routing_reg_findbasins(ib, nbi, nbj, trip, basin, fac, hierarchy, topoind, lshead, diaglalo, & & nb_basin, basin_inbxid, basin_outlet, basin_outtp, basin_sz, basin_bxout, basin_bbout, basin_pts, & & basin_lshead, coast_pts, lontmp, lattmp) ! IMPLICIT NONE ! !! INPUT VARIABLES INTEGER(i_std), INTENT(in) :: ib !! INTEGER(i_std), INTENT(in) :: nbi !! Number of point in x within the grid (unitless) INTEGER(i_std), INTENT(in) :: nbj !! Number of point in y within the grid (unitless) REAL(r_std), INTENT(in) :: hierarchy(:,:) !! REAL(r_std), INTENT(in) :: fac(:,:) !! REAL(r_std), INTENT(in) :: lshead(:,:) REAL(r_std), DIMENSION(:,:), INTENT(in) :: diaglalo !! Point (in Lat/Lon) where diagnostics will be printed. ! ! Modified ! INTEGER(i_std), INTENT(inout) :: trip(:,:) !! The trip field (unitless) INTEGER(i_std), INTENT(inout) :: basin(:,:) !! REAL(r_std), INTENT(inout) :: topoind(:,:) !! Topographic index of the residence time (m) ! ! For debug and get coordinate of river outlet ! REAL(r_std), INTENT(in) :: lontmp(:,:) !! Longitude REAL(r_std), INTENT(in) :: lattmp(:,:) !! Latitude ! !! OUTPUT VARIABLES INTEGER(i_std), INTENT(out) :: nb_basin !! Number of sub-basins (unitless) INTEGER(i_std), INTENT(out) :: basin_inbxid(nbvmax) !! REAL(r_std), INTENT(out) :: basin_outlet(nbvmax,2) !! REAL(r_std), INTENT(out) :: basin_outtp(nbvmax) !! INTEGER(i_std), INTENT(out) :: basin_sz(nbvmax) !! INTEGER(i_std), INTENT(out) :: basin_bxout(nbvmax) !! REAL(r_std), INTENT(out) :: basin_lshead(nbvmax) !! INTEGER(i_std), INTENT(out) :: basin_bbout(nbvmax) !! INTEGER(i_std), INTENT(out) :: basin_pts(nbvmax, nbvmax, 2) !! INTEGER(i_std), INTENT(out) :: coast_pts(nbvmax) !! The coastal flow points (unitless) ! !! LOCAL VARIABLES LOGICAL, PARAMETER :: debug=.TRUE. CHARACTER(LEN=7) :: fmt !! CHARACTER(LEN=9) :: fmtr !! ! LOGICAL :: newpoint !! (true/false) INTEGER(i_std), DIMENSION(8,2) :: inc !! INTEGER(i_std), DIMENSION(nbi,nbj) :: outtmp !! ! INTEGER(i_std) :: nbb, p, cnt, ibas !! To count number of subbasin INTEGER(i_std) :: il, jl, ill, jll, ij, iz!! INTEGER(i_std) :: totsz, checksz !! ! INTEGER(i_std) :: tbname(nbvmax) !! INTEGER(i_std) :: tsz(nbvmax) !! INTEGER(i_std) :: tpts(nbvmax,nbvmax,2) !! INTEGER(i_std) :: toutdir(nbvmax) !! INTEGER(i_std) :: toutbas(nbvmax) !! INTEGER(i_std) :: toutloc(nbvmax,2) !! REAL(r_std) :: tlon(nbvmax) !! REAL(r_std) :: tlat(nbvmax) !! REAL(r_std) :: touttp(nbvmax) !! REAL(r_std) :: toutlshead(nbvmax) !! ! INTEGER(i_std) :: tmpsz(nbvmax) !! INTEGER(i_std) :: ip, jp, jpp(1), ipb !! INTEGER(i_std) :: ie, je !! INTEGER(i_std) :: sortind(nbvmax) !! ! INTEGER(i_std) :: itrans !! INTEGER(i_std) :: trans(nbvmax) !! ! INTEGER(i_std) :: new_nb !! Number of sub-basins (unitless) INTEGER(i_std) :: mp !! Number of dividing sub-basin in mainstream (unitless) INTEGER(i_std) :: new_bname(nbvmax) !! INTEGER(i_std) :: new_outdir(nbvmax) !! REAL(r_std) :: new_heading(nbvmax)!! INTEGER(i_std) :: new_outbas(nbvmax) !! INTEGER(i_std) :: new_outloc(nbvmax,2) !! REAL(r_std) :: new_lon(nbvmax) !! REAL(r_std) :: new_lat(nbvmax) !! INTEGER(i_std) :: new_sz(nbvmax) !! INTEGER(i_std) :: new_pts(nbvmax, nbvmax, 2) !! INTEGER(i_std) :: oldorder, neworder !! ! INTEGER(i_std) :: option !! Option to calculate topoindex !_ ================================================================================================================================ ! IF ( debug .AND. routing_diagbox(ib, diaglalo) ) THEN WRITE(numout,*) "Point: ", ib WRITE(numout,*) "Coord:", lalo(ib,2), lalo(ib,1) WRITE(numout,*) "Per: ", nbi, nbj ! WRITE(fmt,"('(',I3,'I6)')") nbi WRITE(numout,*) '+++++++++++++++++++ TRIP IN FINDBASINS ++++++++++++++++++++' DO je=1,nbj WRITE(numout,fmt) trip(1:nbi,je) ENDDO ENDIF ! ! 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 ! ! 1.0 Find number of outflow point (>= 97) ! So basically, we get nb_basin, basin_bxout, basin_inbxid ! (with nbb, toutdir, tbname) ! nbb = 0 tbname(:) = undef_int toutdir(:) = 0 toutloc(:,:) = 0 outtmp(:,:) = -1 ! toutbas(:) = undef_int touttp(:) = undef_int toutlshead(:) = undef_sechiba ! DO ip=1,nbi DO jp=1,nbj ! ! IF ( basin(ip,jp) .LT. undef_int) THEN ! Because, there is maybe difference between undef_int and undef in routing.nc !IF ( trip(ip,jp) .GE. 97 .AND. trip(ip,jp) .LT. undef_int ) THEN IF ( trip(ip,jp) .GE. 97 .AND. trip(ip,jp) .LT. 109 ) THEN nbb = nbb + 1 IF ( nbb .GT. nbvmax ) CALL ipslerr_p(3,'routing_reg_findbasins','nbvmax too small','first section','') tbname(nbb) = basin(ip,jp) toutdir(nbb) = trip(ip,jp)-100 toutlshead(nbb) = lshead(ip,jp) touttp(nbb) = trip(ip,jp) toutloc(nbb,1) = ip toutloc(nbb,2) = jp tlat(nbb) = lattmp(ip,jp) tlon(nbb) = lontmp(ip,jp) outtmp(ip,jp) = nbb ENDIF ENDIF ENDDO ENDDO ! ! 2.0 Follow the flow of the water to get size and gather point for each ! sub-basin. So we get basin_sz, basin_pts (with tsz, tpts) ! totsz = 0 tsz(:) = 0 tpts(:,:,:) = 0 ! DO ip=1,nbi DO jp=1,nbj IF ( basin(ip,jp) .LT. undef_int) THEN ! Think about basin .LT. undef_int here !IF ( trip(ip,jp) .GT. 0 .AND. trip(ip,jp) .LT. undef_int ) THEN IF ( trip(ip,jp) .GT. 0 .AND. trip(ip,jp) .LT. 109 ) THEN p = trip(ip,jp) cnt = 0 il = ip ; jl = jp newpoint = .TRUE. DO WHILE ( p .GT. 0 .AND. p .LT. 9 .AND. cnt .LT. nbi*nbj .AND. newpoint ) cnt = cnt + 1 ill = il + inc(p,1) jll = jl + inc(p,2) il = ill ; jl = jll p = trip(il,jl) IF ( outtmp(il,jl) .NE. -1 ) newpoint = .FALSE. ENDDO ! ! IF ( cnt .EQ. nbi*nbj) THEN IF ( debug ) THEN WRITE(numout,*) "Error: cnt .EQ. nbi*nbj, ib= ",ib WRITE(numout,*) "Point: ", ip, jp WRITE(numout,*) "Per: ", nbi, nbj WRITE(numout,*) "Number of outlet: ", nbb, cnt WRITE(numout,*) '+++++++++++++++++++ OUTLET IN FINDBASINS ++++++++++++++++++++' WRITE(fmt,"('(',I3,'I6)')") nbi DO je=1,nbj WRITE(numout,fmt) outtmp(1:nbi,je) ENDDO ! WRITE(numout,*) '+++++++++++++++++++ TRIP IN FINDBASINS ++++++++++++++++++++' DO je=1,nbj WRITE(numout,fmt) trip(1:nbi,je) ENDDO ! WRITE(numout,*) '+++++++++++++++++++ BASIN IDs IN FINDBASINS ++++++++++++++++++++' DO je=1,nbj WRITE(numout,fmt) basin(1:nbi,je) ENDDO ! WRITE(numout,*) '+++++++++++++++++++ LONGITUDE IN FINDBASINS ++++++++++++++++++++' WRITE(fmtr,"('(',I3,'F8.1)')") nbi DO je=1,nbj WRITE(numout,fmtr) lontmp(1:nbi,je) ENDDO ! WRITE(numout,*) '+++++++++++++++++++ LATITUDE IN FINDBASINS ++++++++++++++++++++' WRITE(fmtr,"('(',I3,'F8.1)')") nbi DO je=1,nbj WRITE(numout,fmtr) lattmp(1:nbi,je) ENDDO ! WRITE(numout,*) '+++++++++++++++++++ HIERARCHY IN FINDBASINS ++++++++++++++++++++' WRITE(fmtr,"('(',I3,'F8.1)')") nbi DO je=1,nbj WRITE(numout,fmtr) hierarchy(1:nbi,je) ENDDO CALL FLUSH(numout) ENDIF CALL ipslerr_p(3,'routing_reg_findbasins','We could not route point','','') ENDIF ! IF ( outtmp(il,jl) .NE. -1 ) THEN ibas = outtmp(il,jl) tsz(ibas) = tsz(ibas) + 1 IF ( tsz(ibas) .GT. nbvmax ) CALL ipslerr_p(3,'routing_reg_findbasins','nbvmax too small','second section','') tpts(ibas, tsz(ibas), 1) = ip tpts(ibas, tsz(ibas), 2) = jp outtmp(ip,jp) = outtmp(il,jl) ELSE IF ( debug ) THEN WRITE(numout,*) "Error: outtmp(il,jl) .EQ. -1 ", ib WRITE(numout,*) "Point: ", ip, jp WRITE(numout,*) "Per: ", nbi, nbj WRITE(numout,*) "Number of outlet: ", nbb WRITE(numout,*) '+++++++++++++++++++ OUTLET IN FINDBASINS ++++++++++++++++++++' WRITE(fmt,"('(',I3,'I6)')") nbi DO je=1,nbj WRITE(numout,fmt) outtmp(1:nbi,je) ENDDO ! WRITE(numout,*) '+++++++++++++++++++ TRIP IN FINDBASINS ++++++++++++++++++++' DO je=1,nbj WRITE(numout,fmt) trip(1:nbi,je) ENDDO ! WRITE(numout,*) '+++++++++++++++++++ BASIN IDs IN FINDBASINS ++++++++++++++++++++' DO je=1,nbj WRITE(numout,fmt) basin(1:nbi,je) ENDDO ! WRITE(numout,*) '+++++++++++++++++++ LONGITUDE IN FINDBASINS ++++++++++++++++++++' WRITE(fmtr,"('(',I3,'F8.1)')") nbi DO je=1,nbj WRITE(numout,fmtr) lontmp(1:nbi,je) ENDDO ! WRITE(numout,*) '+++++++++++++++++++ LATITUDE IN FINDBASINS ++++++++++++++++++++' WRITE(fmtr,"('(',I3,'F8.1)')") nbi DO je=1,nbj WRITE(numout,fmtr) lattmp(1:nbi,je) ENDDO ! WRITE(numout,*) '+++++++++++++++++++ HIERARCHY IN FINDBASINS ++++++++++++++++++++' WRITE(fmtr,"('(',I3,'F8.1)')") nbi DO je=1,nbj WRITE(numout,fmtr) hierarchy(1:nbi,je) ENDDO ENDIF CALL FLUSH(numout) WRITE (numout,*) ' outtmp(il,jl):', outtmp(il,jl) WRITE (numout,*) ' trip(il,jl):', trip(il,jl) CALL ipslerr_p(3,'routing_reg_findbasins','We could not find right outlet','','') ENDIF ! totsz = totsz + 1 ENDIF ! Testing this IF ENDIF ENDDO ENDDO ! IF ( debug .AND. routing_diagbox(ib, diaglalo) ) THEN WRITE(numout,*) "=== To check grid before dividing sub-basin ===" WRITE(numout,*) "Number of points: ", nbi, nbj WRITE(numout,*) "Number of outlet: ", nbb WRITE(numout,*) '+++++++++++++++++++ OUTLET IN FINDBASINS ++++++++++++++++++++' WRITE(fmt,"('(',I3,'I6)')") nbi DO je=1,nbj WRITE(numout,fmt) outtmp(1:nbi,je) ENDDO ! WRITE(numout,*) '+++++++++++++++++++ TRIP IN FINDBASINS ++++++++++++++++++++' DO je=1,nbj WRITE(numout,fmt) trip(1:nbi,je) ENDDO ENDIF ! ! 4.0 Use the level 1 of the Pfafstetter coding system to divide the large ! sub-grid basin. The large sub-grid basin here is defined by the area ! larger than 9% of grid box area. This threshold should be discussed later. ! ! Make a loop through all sub-grid basins to find if there is big sub-grid basin. ! ! We need at least one subbasin with size greater than 4 (to have at least one ! tributary). If you run ORCHIDEE with high resolution, such as 5 km, you ! don't need to divide the subbasin too small. ! ! Check if there is still subbasin with area larger than [maxpercent] % of total area ! DO WHILE ( COUNT(tsz(:)/REAL(totsz) >= maxpercent/REAL(100) ) > 0 .AND. COUNT(tsz(:) >= 4 ) > 0 ) ! cnt = 0 DO ibas = 1, nbb ! Only take the subbasin is bigger than 4 points IF (( tsz(ibas)/REAL(totsz) .GE. maxpercent/REAL(100) ) .AND. ( tsz(ibas) .GE. 4 )) THEN cnt = cnt + 1 trans(cnt) = ibas ENDIF ENDDO ! ! cnt should be greater or equal 1 here DO ipb = 1, cnt ! ibas = trans(ipb) ! CALL routing_reg_divbas(nbi, nbj, ibas, toutloc(ibas,1), toutloc(ibas,2), tsz(ibas), toutbas(ibas), toutdir(ibas), & & toutlshead(ibas), tpts, trip, basin, fac, lontmp, lattmp, & & new_nb, mp, new_bname, new_outdir, new_heading, new_outbas, new_lat, new_lon, new_outloc, new_sz, new_pts) ! ! If the dividing step is ok ! IF ( new_nb .NE. 0 ) THEN ! ! If we have split the basin then we need to change the old one ! tbname(ibas) = new_bname(1) toutdir(ibas) = new_outdir(1) toutlshead(ibas) = new_heading(1) ! IF ( new_outbas(1) .NE. undef_int ) THEN IF ( mp .NE. 1 ) THEN toutbas(ibas) = new_outbas(1)+nbb-1 ELSE toutbas(ibas) = new_outbas(1) ENDIF ENDIF ! toutloc(ibas,:) = new_outloc(1,:) tlat(ibas) = new_lat(1) tlon(ibas) = new_lon(1) ! touttp(ibas) = trip(new_outloc(1,1),new_outloc(1,2)) ! tsz(ibas) = new_sz(1) tpts(ibas,:,:) = new_pts(1,:,:) DO p = 1, tsz(ibas) outtmp(tpts(ibas,p,1),tpts(ibas,p,2)) = ibas ENDDO ! ! Add the new sub-grid basins to the end of the list ! IF ( nbb+new_nb-1 .LE. nbvmax) THEN DO ip = 2, new_nb tbname(nbb+ip-1) = new_bname(ip) toutdir(nbb+ip-1) = new_outdir(ip) toutlshead(nbb+ip-1) = new_heading(ip) ! ! Trung: Should come back here soon (small bug can come from ! Trung: this point) ! IF ( new_outbas(ip) .EQ. 1 .AND. ip .NE. mp ) THEN toutbas(nbb+ip-1) = ibas ELSE IF ( ip .NE. mp ) THEN toutbas(nbb+ip-1) = new_outbas(ip)+nbb-1 ELSE toutbas(nbb+ip-1) = new_outbas(ip) ENDIF ! toutloc(nbb+ip-1,:) = new_outloc(ip,:) tlat(nbb+ip-1) = new_lat(ip) tlon(nbb+ip-1) = new_lon(ip) ! touttp(nbb+ip-1) = trip(new_outloc(ip,1),new_outloc(ip,2)) ! tsz(nbb+ip-1) = new_sz(ip) tpts(nbb+ip-1,:,:) = new_pts(ip,:,:) DO p = 1, tsz(nbb+ip-1) outtmp(tpts(nbb+ip-1,p,1),tpts(nbb+ip-1,p,2)) = nbb+ip-1 ENDDO ENDDO nbb = nbb+new_nb-1 ELSE WRITE(numout,*) 'Increase nbvmax. It is too small to contain all the basins (routing_reg_findbasins)' CALL ipslerr_p(3,'routing_reg_findbasins','Increase nbvmax.','It is too small to contain all the basins','') ENDIF ! ELSE ! ! Print some information if the division failed ! WRITE(numout,*) 'Can not divide sub-basins (routing_reg_findbasins): ',ibas CALL ipslerr_p(3,'routing_reg_findbasins','new_nb = 0','Something is wrong with routing_reg_divbas','') ! ENDIF ENDDO ! ENDDO ! ! 3.0 Gather coastal points ! itrans = 0 coast_pts(:) = undef_int ! ! Get all the points we can collect (negative value of IDs in routing.nc) ! DO ip=1,nbb IF ( tsz(ip) .EQ. 1 .AND. trip(tpts(ip,1,1),tpts(ip,1,2)) .EQ. 98 .AND. tbname(ip) .LT. 0 ) THEN itrans = itrans + 1 trans(itrans) = ip ENDIF ENDDO ! ! Put everything in the first basin (beware of tbname(trans(1)) can be negative) ! IF ( itrans .GT. 1) THEN ipb = trans(1) coast_pts(tsz(ipb)) = tbname(ipb) tbname(ipb) = -1 DO ip=2,itrans tsz(ipb) = tsz(ipb) + 1 coast_pts(tsz(ipb)) = tbname(trans(ip)) tsz(trans(ip)) = 0 tpts(ipb, tsz(ipb), 1) = tpts(trans(ip), 1, 1) tpts(ipb, tsz(ipb), 2) = tpts(trans(ip), 1, 2) ! outtmp(tpts(ipb, tsz(ipb), 1),tpts(ipb, tsz(ipb), 2)) = ipb ENDDO ENDIF ! ! A.0 To check each grid when testing this subroutine ! IF ( debug .AND. routing_diagbox(ib, diaglalo) ) THEN WRITE(numout,*) "=== To check grid after dividing sub-basin ===" WRITE(numout,*) "Number of points: ", nbi, nbj WRITE(numout,*) "Number of outlet: ", nbb WRITE(numout,*) '+++++++++++++++++++ OUTLET IN FINDBASINS ++++++++++++++++++++' WRITE(fmt,"('(',I3,'I6)')") nbi DO je=1,nbj WRITE(numout,fmt) outtmp(1:nbi,je) ENDDO ! WRITE(numout,*) '+++++++++++++++++++ TRIP IN FINDBASINS ++++++++++++++++++++' DO je=1,nbj WRITE(numout,fmt) trip(1:nbi,je) ENDDO ! WRITE(numout,*) '+++++++++++++++++++ BASIN IDs IN FINDBASINS ++++++++++++++++++++' DO je=1,nbj WRITE(numout,fmt) basin(1:nbi,je) ENDDO ! WRITE(numout,*) '+++++++++++++++++++ LONGITUDE IN FINDBASINS ++++++++++++++++++++' WRITE(fmtr,"('(',I3,'F8.1)')") nbi DO je=1,nbj WRITE(numout,fmtr) lontmp(1:nbi,je) ENDDO ! WRITE(numout,*) '+++++++++++++++++++ LATITUDE IN FINDBASINS ++++++++++++++++++++' WRITE(fmtr,"('(',I3,'F8.1)')") nbi DO je=1,nbj WRITE(numout,fmtr) lattmp(1:nbi,je) ENDDO ! WRITE(numout,*) '+++++++++++++++++++ HIERARCHY IN FINDBASINS ++++++++++++++++++++' WRITE(fmtr,"('(',I3,'F8.1)')") nbi DO je=1,nbj WRITE(numout,fmtr) hierarchy(1:nbi,je) ENDDO WRITE(numout,*) "=== ====================================== ===" ! WRITE(numout,*) '+++++++++++++++++++ FAC IN FINDBASINS ++++++++++++++++++++' WRITE(fmtr,"('(',I3,'F8.1)')") nbi DO je=1,nbj WRITE(numout,fmtr) fac(1:nbi,je) ENDDO WRITE(numout,*) "=== ====================================== ===" ENDIF ! ! 5.0 Sort the output by size of the various basins. ! nb_basin = COUNT(tsz(1:nbb) .GT. 0) tmpsz(:) = -1 tmpsz(1:nbb) = tsz(1:nbb) DO ip=1,nbb jpp = MAXLOC(tmpsz(:)) IF ( tsz(jpp(1)) .GT. 0) THEN sortind(ip) = jpp(1) tmpsz(jpp(1)) = -1 ENDIF ENDDO basin_inbxid(1:nb_basin) = tbname(sortind(1:nb_basin)) basin_outlet(1:nb_basin,1) = tlat(sortind(1:nb_basin)) basin_outlet(1:nb_basin,2) = tlon(sortind(1:nb_basin)) basin_outtp(1:nb_basin) = touttp(sortind(1:nb_basin)) basin_sz(1:nb_basin) = tsz(sortind(1:nb_basin)) basin_pts(1:nb_basin,:,:) = tpts(sortind(1:nb_basin),:,:) basin_bxout(1:nb_basin) = toutdir(sortind(1:nb_basin)) basin_lshead(1:nb_basin) = toutlshead(sortind(1:nb_basin)) ! ! Correct toutbas after sorting ! basin_bbout(1:nb_basin) = undef_int DO ip = 1, nb_basin IF ( toutbas(sortind(ip)) .NE. undef_int ) THEN neworder = 0 oldorder = toutbas(sortind(ip)) DO ie = 1, nb_basin IF ( sortind(ie) .EQ. oldorder ) THEN neworder = ie ENDIF ENDDO IF ( neworder .NE. 0 ) THEN basin_bbout(ip) = neworder ELSE WRITE(numout,*) 'ip, sortind(ip) :', ip, sortind(ip) WRITE(numout,*) 'toutbas(sortind(ip)) :', toutbas(sortind(ip)) WRITE(numout,*) 'Found new order :', neworder ! WRITE(numout,*) 'nb_basin & nbb :', nb_basin, nbb WRITE(numout,*) 'sortind(1:nb_basin) :', sortind(1:nb_basin) WRITE(numout,*) 'tsz(sortind(1:nb_basin)) :', tsz(sortind(1:nb_basin)) WRITE(numout,*) 'tsz(1:nbb) :', tsz(1:nbb) WRITE(numout,*) 'toutbas(sortind(1:nb_basin)) :', toutbas(sortind(1:nb_basin)) WRITE(numout,*) 'toutdir(sortind(1:nb_basin)) :', toutdir(sortind(1:nb_basin)) WRITE(numout,*) 'tbname(sortind(1:nb_basin)) :', tbname(sortind(1:nb_basin)) WRITE(numout,*) 'toutbas(1:nbb) :', toutbas(1:nbb) ! CALL ipslerr_p(3,'routing_reg_findbasins','We got problem when sorting toutbas','','') ENDIF ENDIF ENDDO ! IF ( debug .AND. routing_diagbox(ib, diaglalo) ) THEN WRITE(numout,*) "Grid box: ", ib WRITE(numout,*) "Number of points: ", nbi, nbj WRITE(numout,*) mpi_rank, "=== To check coastal point simplify ===" WRITE(numout,*) "Number of outlet: " WRITE(numout,*) " + Before gather: ", nbb WRITE(numout,*) " + After gather: ", nb_basin WRITE(numout,*) mpi_rank, "=== To check dividing sub-basin process ===" ip = COUNT(basin_bxout(1:nb_basin) .GE. -3 .AND. basin_bxout(1:nb_basin) .LE. -1 ) WRITE(numout,*) "Outlet to lake, ocean: ", ip ip = COUNT(basin_bxout(1:nb_basin) .GE. 1 .AND. basin_bxout(1:nb_basin) .LE. 8 ) WRITE(numout,*) "Outlet to other grid box: ", ip ip = COUNT(basin_bxout(1:nb_basin) .EQ. -4) WRITE(numout,*) "Outlet in the same grid box: ", ip ip = COUNT(basin_bbout(1:nb_basin) .NE. undef_int) WRITE(numout,*) "Should be the same value as above: ", ip ENDIF ! ! 6.0 Make simple verification ! What is the size of the region behind each outflow point ? ! We can only check if we have at least as many outflows as basins ! checksz = 0 DO ip=1,nbb checksz = checksz + tsz(ip) ENDDO IF ( checksz .NE. totsz) THEN WRITE(numout,*) 'Water got lost while I tried to find basins' WRITE(numout,*) checksz, ' /= ', totsz CALL ipslerr_p(3,'routing_reg_findbasins','Water got lost while I tried to find basins','','') ENDIF ! 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. -4) IF ( ip .LT. nb_basin ) THEN WRITE(numout,*) 'We found some coastal point here :', ip WRITE(numout,*) 'nb_basin :', nb_basin WRITE(numout,*) 'Basin sized :', basin_sz(1:nb_basin) ENDIF ! ! 7.0 If we calculate topoindex with option 3 ! we need to accumulate the topoind_bx along the flow ! ! Get option for calculating topoindex ! !Config Key = TOPOINDEX !Config Desc = Options to calculate topoindex !Config Def = 3 !Config Help = This flag allows to calculate topoindex with different options !Config Units = [-] ! option = 3 !!$ CALL getin('TOPOINDEX', option) !!$ ! !!$ ! Checking option: until now (May 2016) there is only 3 acceptable options !!$ ! = 1 : topoindex .EQ. 1000. everywhere !!$ ! = 2 : simple average as the method in the older routing scheme !!$ ! = 3 : accumulate the "topoind" from routing.nc along the flow (as Agnes DUCHARNE's suggestion) !!$ ! !!$ IF ( option .LT. 1 .OR. option .GT. 3 ) THEN !!$ WRITE(numout,*) ' You chose wrong option for calculating topoindex: ', option !!$ WRITE(numout,*) ' It should be: 1, 2 or 3' !!$ STOP 'routing_reg_globalize' !!$ ENDIF ! ! Work only using option number 3 ! IF ( option .EQ. 3 ) THEN ! Loop for all sub-basins and each point belongs to that sub-basin DO ij=1,nb_basin DO iz=1,basin_sz(ij) ! Get the coordinate of this point ip = basin_pts(ij,iz,1) jp = basin_pts(ij,iz,2) ! Start to accumulate "topoind" if this is not a outlet point IF ( trip(ip,jp) .GT. 0 .AND. trip(ip,jp) .LT. 97 ) THEN p = trip(ip,jp) cnt = 0 il = ip ; jl = jp DO WHILE ( p .GT. 0 .AND. p .LT. 9 .AND. cnt .LT. nbi*nbj ) cnt = cnt + 1 ill = il + inc(p,1) jll = jl + inc(p,2) ! Of course, if you are careful, you will want to check if ! this point still belongs to current sub-basin. Here I ! assume that the sub-routine routing_reg_divbas worked well. ! So we don't need to check. topoind(basin_pts(ij,iz,1),basin_pts(ij,iz,2))=topoind(basin_pts(ij,iz,1),basin_pts(ij,iz,2))+topoind(ill,jll) ! il = ill ; jl = jll p = trip(il,jl) ENDDO ENDIF ! ENDDO ENDDO ENDIF ! END SUBROUTINE routing_reg_findbasins ! !! ================================================================================================================================ !! SUBROUTINE : routing_reg_divbas !! !>\BRIEF This subroutine divide the large sub-grid basins to smaller ones. It is !! based on The Pfafstetter Coding System for Watershed Identification. You can !! find information here (for example): !! http://ponce.sdsu.edu/pfafstetter_system_revised_presentation.html !! !! DESCRIPTION (definitions, functional, design, flags) : !! !! RECENT CHANGE(S): None !! !! MAIN OUTPUT VARIABLE(S): !! !! REFERENCES : None !! !! FLOWCHART : None !! \n !_ ================================================================================================================================ SUBROUTINE routing_reg_divbas(nbi, nbj, ibas, iloc, jloc, tsz, tout, toutd, headd, & & tpts, trip, basin, fac, lon, lat, & & new_nb, oic, new_bname, new_outdir, new_head, new_outbas, new_lat, new_lon, new_outloc, new_sz, new_pts) ! IMPLICIT NONE ! !! INPUT VARIABLES INTEGER(i_std), INTENT(in) :: nbi !! Number of point in x within the grid (unitless) INTEGER(i_std), INTENT(in) :: nbj !! Number of point in y within the grid (unitless) INTEGER(i_std), INTENT(in) :: ibas !! Order of the basin will be divided INTEGER(i_std), INTENT(in) :: iloc, jloc !! Outlet location ! INTEGER(i_std), INTENT(in) :: tsz !! Size of sub-basin INTEGER(i_std), INTENT(in) :: tout !! Outlet type INTEGER(i_std), INTENT(in) :: toutd !! Outlet direction REAL(r_std), INTENT(in) :: headd INTEGER(i_std), INTENT(in) :: tpts(:,:,:) !! INTEGER(i_std), INTENT(in) :: basin(:,:) !! REAL(r_std), INTENT(in) :: fac(:,:) !! REAL(r_std), INTENT(in) :: lon(:,:), lat(:,:) !! ! !! MODIFIED VARIABLES INTEGER(i_std), INTENT(inout) :: trip(:,:) !! ! !! OUTPUT VARIABLES INTEGER(i_std), INTENT(out) :: new_nb !! Number of sub-basins (unitless) INTEGER(i_std), INTENT(out) :: oic INTEGER(i_std), INTENT(out) :: new_bname(nbvmax) !! INTEGER(i_std), INTENT(out) :: new_outdir(nbvmax) !! REAL(r_std), INTENT(out) :: new_head(nbvmax) !! INTEGER(i_std), INTENT(out) :: new_outbas(nbvmax) !! INTEGER(i_std), INTENT(out) :: new_outloc(nbvmax,2) !! REAL(r_std), INTENT(out) :: new_lon(nbvmax) !! REAL(r_std), INTENT(out) :: new_lat(nbvmax) !! INTEGER(i_std), INTENT(out) :: new_sz(nbvmax) !! INTEGER(i_std), INTENT(out) :: new_pts(nbvmax, nbvmax, 2) !! ! !! LOCAL VARIABLES REAL(r_std), PARAMETER :: flag=-9999. !! LOGICAL, PARAMETER :: debug=.FALSE. LOGICAL, PARAMETER :: checkgrid=.FALSE. CHARACTER(LEN=7) :: fmt !! CHARACTER(LEN=9) :: fmtr !! CHARACTER(LEN=11) :: afmt !! CHARACTER(LEN=13) :: afmtr !! ! ! REAL(r_std), DIMENSION(nbi,nbj) :: factmp !! INTEGER(i_std), DIMENSION(nbi,nbj) :: triptmp !! INTEGER(i_std), DIMENSION(nbi,nbj) :: triptemp !! INTEGER(i_std), DIMENSION(nbi,nbj) :: basintmp !! REAL(r_std), DIMENSION(nbi,nbj) :: lontmp, lattmp !! ! INTEGER(i_std) :: il, jl, ill, jll !! INTEGER(i_std) :: ie, je !! INTEGER(i_std) :: p, cnt, k, l, ic, ik !! INTEGER(i_std) :: ip, isz, checksz !! INTEGER(i_std) :: ff(1) !! ! LOGICAL :: okpoint=.TRUE. ! ! Number of neighbours on the HydroShed grid (regular Lat/Lon) ! INTEGER(i_std), PARAMETER :: nbne=8 INTEGER(i_std), DIMENSION(nbne,2) :: inc !! ! ! Explanation for dimension of below arrays: ! ! + Firstly, I suggest that you should read about the Pfafstetter Coding ! System, here for example: ! http://ponce.sdsu.edu/pfafstetter_system_revised_presentation.html ! ! * Note: It will be easier if you keep in mind that the original ! information given in grid box of nbi*nbj. But here, to deal with main ! stem and tributaries, I need to store them to cnt*nbne (exactly, ! cnt*[nbne-2]). Where, cnt is the length of main stem and nbne is number ! of neighbour points, and usually 2 neighbour points will belongs to ! main stem. ! ! + You will see steps for coding the sub-basin are: ! ! - From the sub-basin outlet, trace upstream along the main stem of the ! river, and get the tributaries. To store location and flow accumulation ! of points belongs to main stream and tributaries: ! INTEGER(i_std) :: main_loc(nbvmax,2), tri_loc(nbne,nbvmax,2) REAL(r_std) :: main_fac(nbvmax), tri_fac(nbne,nbvmax) REAL(r_std) :: tmptri_fac(nbne) ! Flow accumulation of all neighbour points INTEGER(i_std) :: sortedtrifac(nbne) ! And sort of tmptri_fac(nbne) REAL(r_std) :: alltri_fac(nbvmax) ! Sort flow accumulation of all tributaries INTEGER(i_std) :: alltri_loc1(nbvmax), alltri_loc2(nbvmax) ! And their location INTEGER(i_std) :: sortedallfac(nbvmax) ! Sort alltri_fac(nbvmax) ! ! * Note: again, there are few dimensions should be exactly [nbne-1] ! or [nbne-2]. Because when you follow upstream the mainstream, there are ! always at least ONE neighbour points will not belong to TRIBUTARY. ! But to make your life easier, I put all "nbne" here. ! ! - Identify the 4 greatest tributaries: ! REAL(r_std) :: tmpmain_fac(4) ! Flow accumulations of 4 greatest tributaries INTEGER(i_std) :: tmpmain_loc(4,2), sortedmainfac(4) ! And their locations and sorted values INTEGER(i_std) :: numtri(4), usetri_loc(4,4,2) ! Counting and store location when processing these 4 points ! ! - Inter-basins are the watersheds that contribute flow to the main ! stem: there can be 5 inter-basins. So we need 9 dividing points: ! INTEGER(i_std) :: divloc(9,2) ! Location of dividing points (maximum is 9) ! !_ ================================================================================================================================ ! ! 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 ! ! 0.0 Get the information of all the points belongs to the subbasin ! which need to divide. Flag for flow accumulation is -9999 (the value of ! flow accumulation should be always non negative). ! factmp(:,:) = flag basintmp(:,:) = 0 lontmp(:,:) = 0 lattmp(:,:) = 0 triptmp(:,:) = -1 triptemp(:,:) = -1 ! DO isz = 1, tsz il = tpts(ibas, isz, 1) jl = tpts(ibas, isz, 2) ! factmp(il,jl) = fac(il,jl) triptmp(il,jl) = trip(il,jl) triptemp(il,jl) = trip(il,jl) basintmp(il,jl) = basin(il,jl) lontmp(il,jl) = lon(il,jl) lattmp(il,jl) = lat(il,jl) ENDDO ! ! Print out information of grid box ! IF ( checkgrid ) THEN ! WRITE(numout,*) " Routing_reg_divbas: Grid before dividing " WRITE(numout,*) " Size: ", nbi, nbj WRITE(fmt,"('(',I3,'I6)')") nbi WRITE(fmtr,"('(',I3,'F8.1)')") nbi ! WRITE(numout,*) '+++++++++++++++++++ TRIP IN DIVBAS ++++++++++++++++++++' DO je=1,nbj WRITE(numout,fmt) triptmp(1:nbi,je) ENDDO ! WRITE(numout,*) '+++++++++++++++++++ BASIN IDs IN DIVBAS ++++++++++++++++++++' DO je=1,nbj WRITE(numout,fmt) basintmp(1:nbi,je) ENDDO ! WRITE(numout,*) '+++++++++++++++++++ LONGITUDE IN DIVBAS ++++++++++++++++++++' DO je=1,nbj WRITE(numout,fmtr) lontmp(1:nbi,je) ENDDO ! WRITE(numout,*) '+++++++++++++++++++ LATITUDE IN DIVBAS ++++++++++++++++++++' DO je=1,nbj WRITE(numout,fmtr) lattmp(1:nbi,je) ENDDO ! WRITE(numout,*) '+++++++++++++++++++ FAC IN DIVBAS ++++++++++++++++++++' DO je=1,nbj WRITE(numout,fmtr) factmp(1:nbi,je) ENDDO CALL FLUSH(numout) ENDIF ! ! 1.0 If the first version of this subroutine doesn't work well, ! we need to write the validation step of input arrays here (for ! example, the input sub-basin should only has 1 outlet). ! (October 2016: it still works well :)) ! IF ( triptmp(iloc,jloc) .LT. -4 .OR. triptmp(iloc,jloc) .GE. 109 ) THEN CALL ipslerr_p(3,'routing_reg_divbas','The TRIP of outlet is wrong','','') ENDIF IF ( debug ) THEN WRITE(numout,*) 'Basin number: ', ibas WRITE(numout,*) 'Basin size: ', tsz WRITE(numout,*) 'Basin out: ', tout, toutd ENDIF ! ! 2.0 From the outlet of the subbasin, trace back to find the main stream ! and identify the 4 greatest tributaries. If there are less than 4 ! tributraries, take them all. ! ! Trace back upstream from the outlet of the sub-basin ! and define the mainstream and all tributaries ! main_loc(:,:) = 0 main_fac(:) = flag tri_loc(:,:,:) = 0 tri_fac(:,:) = flag ! ! Start! ! cnt = 0 okpoint = .TRUE. il = iloc ; jl = jloc DO WHILE ( il .GT. 0 .AND. il .LE. nbi .AND. jl .GT. 0 .AND. jl .LE. nbj & & .AND. cnt .LT. nbi*nbj .AND. okpoint ) ! ! Count the length of the mainstream ! cnt = cnt + 1 main_loc(cnt,1) = il main_loc(cnt,2) = jl main_fac(cnt) = factmp(il,jl) ! ! Look for [nbne] neighbour points to find the tributaries ! October 2016: so far, nbne = 8 (for 8 directions of TRIP) ! l = 0 DO k = 1, nbne ill = il + inc(k,1) jll = jl + inc(k,2) ! If this point is still in the grid box IF ( ill .GT. 0 .AND. ill .LE. nbi .AND. jll .GT. 0 .AND. jll .LE. nbj ) THEN ! If this point still belongs to the current sub-basin IF ( triptmp(ill,jll) .GT. 0 .AND. triptmp(ill,jll) .LE. nbne ) THEN ie = ill + inc(INT(triptmp(ill,jll)),1) je = jll + inc(INT(triptmp(ill,jll)),2) ! If this neighbour point contributes water for current point IF ( ie .EQ. il .AND. je .EQ. jl ) THEN ! Store all points here l = l + 1 tri_fac(l,cnt) = factmp(ill,jll) tri_loc(l,cnt,1) = ill tri_loc(l,cnt,2) = jll ! Mark the processed points here triptmp(ill,jll) = -2 ENDIF ENDIF ENDIF ENDDO ! ! Only go further if we have found tributaries (there is still neighbour ! point flows into the current point) ! IF ( l > 0 ) THEN ! ! Sort the FAC of all neighbour points tmptri_fac(:) = flag sortedtrifac(:) = 0 ! tmptri_fac(:) = tri_fac(:,cnt) DO k = 1, l ff = MAXLOC(tmptri_fac) sortedtrifac(k) = ff(1) tmptri_fac(ff(1)) = flag ENDDO ! Continue trace upstream: the neighbour point with highest flow ! accumulation will belongs to main stream. We move to this point ! (change value of il, jl). IF ( tri_fac(sortedtrifac(1),cnt) .NE. flag ) THEN il = tri_loc(sortedtrifac(1),cnt,1) jl = tri_loc(sortedtrifac(1),cnt,2) ! IF ( checkgrid ) THEN WRITE(numout,*) "cnt = ", cnt, " fac = ", tri_fac(sortedtrifac(1),cnt) ENDIF ! tri_fac(sortedtrifac(1),cnt) = flag ELSE okpoint = .FALSE. ENDIF ELSE okpoint = .FALSE. ENDIF ! IF ( cnt .EQ. nbi*nbj) THEN IF ( debug ) THEN WRITE(numout,*) "Error: cnt .EQ. nbi*nbj " WRITE(numout,*) "Point: ", il, jl WRITE(numout,*) "Per: ", nbi, nbj WRITE(fmt,"('(',I3,'I6)')") nbi WRITE(fmtr,"('(',I3,'F8.1)')") nbi ! WRITE(numout,*) '+++++++++++++++++++ TRIP IN DIVBAS ++++++++++++++++++++' DO je=1,nbj WRITE(numout,fmt) triptemp(1:nbi,je) ENDDO ! WRITE(numout,*) '+++++++++++++++++++ BASIN IDs IN DIVBAS ++++++++++++++++++++' DO je=1,nbj WRITE(numout,fmt) basintmp(1:nbi,je) ENDDO ! WRITE(numout,*) '+++++++++++++++++++ LONGITUDE IN DIVBAS ++++++++++++++++++++' DO je=1,nbj WRITE(numout,fmtr) lontmp(1:nbi,je) ENDDO ! WRITE(numout,*) '+++++++++++++++++++ LATITUDE IN DIVBAS ++++++++++++++++++++' DO je=1,nbj WRITE(numout,fmtr) lattmp(1:nbi,je) ENDDO ! WRITE(numout,*) '+++++++++++++++++++ FAC IN DIVBAS ++++++++++++++++++++' DO je=1,nbj WRITE(numout,fmtr) factmp(1:nbi,je) ENDDO CALL FLUSH(numout) ENDIF CALL ipslerr_p(3,'routing_reg_divbas','We could not route point','','') ENDIF ! ENDDO IF ( debug ) WRITE (numout,*) 'Length of mainstream: ', cnt ! ! 3.0 Sort the flow accumulation of the tributaries and find the ! divided location. ! alltri_fac(:) = flag alltri_loc1(:) = 0 alltri_loc2(:) = 0 divloc(:,:) = 0 ! ip = 0 DO k = 1, cnt ! Actually, we should only DO l = 1, nbne-2 ! Because, there are usually 2 neighbour points belong to main stream ! and maximum (nbne-2) neighbour points can be tributaries. ! DO l = 1, nbne-2 DO l = 1, nbne IF ( tri_fac(l,k) .NE. flag ) THEN ip = ip + 1 alltri_fac(ip) = tri_fac(l,k) ! Attention: here we store location of each tributary in main stream ! ( k <= cnt ) and its direction compared with main stream ( l <= nbne ) alltri_loc1(ip) = k alltri_loc2(ip) = l ENDIF ENDDO ENDDO IF ( debug ) WRITE (numout,*) 'Number of tributaries: ', ip ! If we have at least one tributary IF ( ip .GT. 0 ) THEN ! Original output oic = 0 ! sortedallfac(:) = 0 DO k = 1, ip ff = MAXLOC(alltri_fac) sortedallfac(k) = ff(1) alltri_fac(ff(1)) = flag ENDDO ! Get the points in the mainstream tmpmain_fac(:) = flag tmpmain_loc(:,:) = 0 numtri(:) = 0 usetri_loc(:,:,:) = 0 ! Using l to count the actual number of points in main stream need to use ! for dividing (not always equal 4 as ideal case). l = 0 ! According to Pfafstetter coding, we only look for 4 greatest ! tributaries: DO k = 1, 4 ! IF ( sortedallfac(k) .NE. 0 ) THEN ! ik = alltri_loc1(sortedallfac(k)) il = alltri_loc2(sortedallfac(k)) ! Starting store location and flow accumulation of points in main stream: IF ( COUNT ( tmpmain_fac(:) .EQ. main_fac(ik) ) .EQ. 0 ) THEN l = l + 1 tmpmain_loc(l,1) = main_loc(ik,1) tmpmain_loc(l,2) = main_loc(ik,2) tmpmain_fac(l) = main_fac(ik) ! numtri(l) = 1 usetri_loc(numtri(l),l,1) = tri_loc(il,ik,1) usetri_loc(numtri(l),l,2) = tri_loc(il,ik,2) ! ELSE ! If this point in main stream have more than 1 tributary ! belongs to top 4 greatest tributaries numtri(l) = numtri(l) + 1 usetri_loc(numtri(l),l,1) = tri_loc(il,ik,1) usetri_loc(numtri(l),l,2) = tri_loc(il,ik,2) ! ENDIF ENDIF ENDDO ! Sort the mainstream points sortedmainfac(:) = 0 DO k = 1, l ff = MAXLOC(tmpmain_fac) sortedmainfac(k) = ff(1) tmpmain_fac(ff(1)) = flag ENDDO ! ic = 0 DO k = l, 1, -1 IF ( tmpmain_loc(sortedmainfac(k),1) .NE. 0 ) THEN ic = ic + 1 divloc(ic,1) = tmpmain_loc(sortedmainfac(k),1) divloc(ic,2) = tmpmain_loc(sortedmainfac(k),2) ! new_outbas(ic) = ic + 1 ! ENDIF ENDDO ! If the last divide point is not the original outlet IF ( divloc(ic,1) .NE. iloc .OR. divloc(ic,2) .NE. jloc ) THEN ic = ic + 1 divloc(ic,1) = iloc divloc(ic,2) = jloc ENDIF ! new_outbas(ic) = tout ! Number of mainstream sub-basin oic = ic ! DO k = l, 1, -1 ! !DO il = 1, nbne-2 ! Maximum 4 tributaries: DO il = 1, 4 IF ( usetri_loc(il,sortedmainfac(k),1) .NE. 0 ) THEN ic = ic + 1 divloc(ic,1) = usetri_loc(il,sortedmainfac(k),1) divloc(ic,2) = usetri_loc(il,sortedmainfac(k),2) ! !new_outbas(ic) = new_outbas(l-k+1) new_outbas(ic) = l-k+1 ! ENDIF ENDDO ENDDO ! Save number of dividing points l = ic ! ! If we don't have any tributary ELSE oic = 2 l = 2 cnt = cnt/2 + 1 divloc(1,1) = main_loc(cnt,1) divloc(1,2) = main_loc(cnt,2) divloc(2,1) = iloc divloc(2,2) = jloc ! new_outbas(1) = 2 new_outbas(2) = tout ! ENDIF IF ( debug ) WRITE (numout,*) 'Number of new sub-basin: ', l, oic ! ! Now, cut it ! Cut the sub-basin ! Release The Kraken ! ! new_nb = l ! DO ik = 1, l ! A small trick here IF ( oic .EQ. 2 .AND. l .EQ. 2 ) THEN k = ik ELSE IF ( ik .LE. (l - oic) ) THEN k = ik + oic ELSE k = ik - (l - oic) ENDIF ENDIF ! I'm so stupid careful here with this IF IF ( divloc(k,1) .NE. 0 .AND. divloc(k,2) .NE. 0 ) THEN ! new_sz(k) = 1 new_pts(k,new_sz(k),1) = divloc(k,1) new_pts(k,new_sz(k),2) = divloc(k,2) ! new_bname(k) = basin(divloc(k,1),divloc(k,2)) new_outloc(k,1) = divloc(k,1) new_outloc(k,2) = divloc(k,2) new_lon(k) = lon(divloc(k,1),divloc(k,2)) new_lat(k) = lat(divloc(k,1),divloc(k,2)) ! ! Change the TRIP here for routing_reg_globalize (don't change the ! original outlet before dividing). Change new_outdir for ! routing_reg_linkup. ! IF ( k .NE. oic ) THEN trip(divloc(k,1), divloc(k,2)) = -4 new_outdir(k) = -4 new_head(k) = undef_sechiba ELSE new_outdir(k) = toutd new_head(k) = headd ENDIF ! ! triptemp(divloc(k,1), divloc(k,2)) = -1 ! DO isz = 1, tsz il = tpts(ibas, isz, 1) jl = tpts(ibas, isz, 2) IF ( triptemp(il,jl) .GT. 0 .AND. triptemp(il,jl) .LT. 9 ) THEN p = triptemp(il,jl) cnt = 0 ie = il ; je = jl okpoint = .TRUE. DO WHILE ( p .GT. 0 .AND. p .LT. 9 .AND. cnt .LT. nbi*nbj .AND. okpoint ) cnt = cnt + 1 ill = ie + inc(p,1) jll = je + inc(p,2) ie = ill ; je = jll p = trip(ie,je) IF ( ie .EQ. divloc(k,1) .AND. je .EQ. divloc(k,2) ) THEN okpoint = .FALSE. new_sz(k) = new_sz(k) + 1 new_pts(k,new_sz(k),1) = il new_pts(k,new_sz(k),2) = jl triptemp(tpts(ibas, isz, 1), tpts(ibas, isz, 2)) = -1 ENDIF ENDDO ENDIF ENDDO ENDIF ENDDO ! ! 4.0 Make simple verification ! IF ( new_nb .GE. 0 ) THEN checksz = 0 DO ip=1,new_nb checksz = checksz + new_sz(ip) ENDDO ! IF ( checksz .NE. tsz) THEN WRITE(numout,*) ' Water got lost while I tried to divide basins' WRITE(numout,*) ' Number of new sub-basin : ', new_nb ip = COUNT(divloc(:,1) .NE. 0) WRITE(numout,*) ' Number of mainstrem sub-basin : ', oic WRITE(numout,*) ' Number of tributaries : ', ip - oic WRITE(numout,*) checksz, ' /= ', tsz ! WRITE(afmt,"('(A10,', I3,'I6)')") ip WRITE(afmtr,"('(A10,', I3,'F8.1)')") ip ! WRITE(numout,afmt) 'new_sz = ', new_sz(1:ip) WRITE(numout,*) 'new_sz(oic) = ', new_sz(oic) WRITE(numout,afmt) 'new_bname', new_bname(1:ip) WRITE(numout,afmt) 'new_outdir', new_outdir(1:ip) WRITE(numout,afmt) 'new_outbas', new_outbas(1:ip) WRITE(numout,afmt) 'new_outloc1', new_outloc(1:ip,1) WRITE(numout,afmt) 'new_outloc2', new_outloc(1:ip,2) WRITE(numout,afmtr) 'new_lon', new_lon(1:ip) WRITE(numout,afmtr) 'new_lat', new_lat(1:ip) ! ip = COUNT(triptemp(1:nbi,1:nbj) .EQ. -1) WRITE(numout,*) 'Unprocessed points :', nbi*nbj-ip ! WRITE(fmt,"('(',I3,'I6)')") nbi WRITE(fmtr,"('(',I3,'F8.1)')") nbi ! WRITE(numout,*) '+++++++++++++++++++ TRIP IN DIVBAS ++++++++++++++++++++' DO je=1,nbj WRITE(numout,fmt) triptemp(1:nbi,je) ENDDO ! CALL ipslerr_p(3,'routing_reg_divbas','Water got lost while I tried to divide sub-basins','','') ENDIF ENDIF ! ip = COUNT(triptemp(1:nbi,1:nbj) .EQ. -1) IF ( ip .NE. nbi*nbj ) THEN WRITE(numout,*) 'Unprocessed points :', nbi*nbj-ip ENDIF ! END SUBROUTINE routing_reg_divbas ! ! !! ================================================================================================================================ !! SUBROUTINE : routing_reg_globalize !! !>\BRIEF This subroutine puts the basins found for grid box in the global map. !! Connection can only be made later when all information is together. !! !! DESCRIPTION (definitions, functional, design, flags) : None !! !! RECENT CHANGE(S): None !! !! MAIN OUTPUT VARIABLE(S): !! 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 !! !! REFERENCES : None !! !! FLOWCHART : None !! \n !_ ================================================================================================================================ SUBROUTINE routing_reg_globalize(nbpt, ib, neighbours, area_bx, lon_bx, lat_bx, trip_bx, hierarchy_bx, & & orog_bx, floodp_bx, fac_bx, topoind_bx, & & min_topoind, nb_basin, basin_inbxid, basin_outlet, basin_outtp, basin_sz, basin_pts, basin_bxout, & & basin_bbout, lshead, coast_pts, nwbas, basin_count, basin_notrun, basin_area, basin_cg, basin_hierarchy, & & basin_orog, basin_floodp, & & basin_fac, basin_topoind, basin_id, basin_coor, basin_type, basin_flowdir, basin_lshead, & & outflow_grid, outflow_basin, nbcoastal, coastal_basin) ! IMPLICIT NONE ! !! INPUT VARIABLES INTEGER(i_std), INTENT (in) :: nbpt !! Domain size (unitless) INTEGER(i_std), INTENT (in) :: ib !! Current basin (unitless) INTEGER(i_std), INTENT(in) :: neighbours(nbpt,NbNeighb)!! Vector of neighbours for each grid point !! (1=North and then clockwise) !! LOCAL VARIABLES REAL(r_std), DIMENSION(nbxmax,nbxmax) :: area_bx !! Area of each small box in the grid box (m^2) REAL(r_std), DIMENSION(nbxmax,nbxmax) :: lon_bx !! Longitude of each small box in the grid box REAL(r_std), DIMENSION(nbxmax,nbxmax) :: lat_bx !! Latitude of each small box in the grid box INTEGER(i_std), DIMENSION(nbxmax,nbxmax) :: trip_bx !! The trip field for each of the smaller boxes (unitless) REAL(r_std), DIMENSION(nbxmax,nbxmax) :: hierarchy_bx !! Level in the basin of the point REAL(r_std), DIMENSION(nbxmax,nbxmax) :: orog_bx !! Level in the basin of the point REAL(r_std), DIMENSION(nbxmax,nbxmax) :: floodp_bx !! Level in the basin of the point REAL(r_std), DIMENSION(nbxmax,nbxmax) :: fac_bx !! Level in the basin of the point REAL(r_std), DIMENSION(nbxmax,nbxmax) :: topoind_bx !! Topographic index of the residence time for each of the smaller boxes (m) REAL(r_std) :: min_topoind !! The current minimum of topographic index (m) INTEGER(i_std) :: nb_basin !! Number of sub-basins (unitless) INTEGER(i_std), DIMENSION(nbvmax) :: basin_inbxid, basin_sz !! ID of basin, number of points in the basin REAL(r_std), DIMENSION(nbvmax,2) :: basin_outlet !! Outlet coordinate of each subgrid basin (lat,lon) REAL(r_std), DIMENSION(nbvmax) :: basin_outtp !! 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), DIMENSION(nbvmax) :: basin_bbout !! outflow sub-basin REAL(r_std), DIMENSION(nbvmax) :: lshead !! Large scale heading of outflow. INTEGER(i_std) :: coast_pts(nbvmax) !! The coastal flow points (unitless) ! global maps INTEGER(i_std) :: nwbas !! INTEGER(i_std), DIMENSION(nbpt) :: basin_count !! INTEGER(i_std), DIMENSION(nbpt) :: basin_notrun !! INTEGER(i_std), DIMENSION(nbpt,nwbas) :: basin_id !! REAL(r_std), DIMENSION(nbpt,nwbas,2) :: basin_coor !! Coordinates of the outflow point REAL(r_std), DIMENSION(nbpt,nwbas) :: basin_type !! INTEGER(i_std), DIMENSION(nbpt,nwbas) :: basin_flowdir !! Water flow directions in the basin (unitless) REAL(r_std), DIMENSION(nbpt,nwbas) :: basin_area !! REAL(r_std), DIMENSION(nbpt,nwbas,2) :: basin_cg !! Centre of gravity of the HTU in latitude, longitude REAL(r_std), DIMENSION(nbpt,nwbas) :: basin_hierarchy !! REAL(r_std), DIMENSION(nbpt,nwbas) :: basin_orog !! REAL(r_std), DIMENSION(nbpt,nwbas) :: basin_floodp !! REAL(r_std), DIMENSION(nbpt,nwbas) :: basin_fac !! REAL(r_std), DIMENSION(nbpt,nwbas) :: basin_topoind !! Topographic index of the residence time for a basin (m) REAL(r_std), DIMENSION(nbpt,nwbas) :: basin_lshead !! Large scale heading out of the grid box. INTEGER(i_std), DIMENSION(nbpt,nwbas) :: outflow_grid !! Type of outflow on the grid box (unitless) INTEGER(i_std), DIMENSION(nbpt,nwbas) :: outflow_basin !! INTEGER(i_std), DIMENSION(nbpt) :: nbcoastal !! INTEGER(i_std), DIMENSION(nbpt,nwbas) :: coastal_basin !! ! INTEGER(i_std) :: ij, iz !! Indices (unitless) CHARACTER(LEN=4) :: hierar_method = 'OUTP' !! ! ! To calculate topoindex in new method ! INTEGER(i_std) :: option !! Option to calculate topoindex REAL(r_std) :: area_frac !! Area fraction for calculate topoindex ! LOGICAL, PARAMETER :: debug = .FALSE. !! Add new logical variable to debug ! !_ ================================================================================================================================ ! ! Get option for calculating topoindex ! !Config Key = TOPOINDEX !Config Desc = Options to calculate topoindex !Config Def = 3 !Config Help = This flag allows to calculate topoindex with different options !Config Units = [-] ! option = 3 !!$ CALL getin('TOPOINDEX', option) ! ! Checking option: until now (May 2016) there is only 3 acceptable options ! = 1 : topoindex .EQ. 1000. everywhere ! = 2 : simple average as the method in the older routing scheme ! = 3 : accumulate the "topoind" from routing.nc along the flow (as Agnes DUCHARNE's suggestion) ! IF ( option .LT. 1 .OR. option .GT. 3 ) THEN WRITE(numout,*) ' You chose wrong option for calculating topoindex: ', option WRITE(numout,*) ' It should be: 1, 2 or 3' STOP 'routing_reg_globalize' ENDIF ! 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 WRITE(numout,*) 'nwbas=',nwbas WRITE(numout,*) 'basin_count(ib)=',basin_count(ib) ! call ipslerr_p(3,'routing_reg_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) ! basin_coor(ib,basin_count(ib),1) = basin_outlet(ij,1) basin_coor(ib,basin_count(ib),2) = basin_outlet(ij,2) ! basin_type(ib,basin_count(ib)) = basin_outtp(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) = zero basin_cg(ib,ij,:) = zero basin_hierarchy(ib,ij) = zero basin_fac(ib,ij) = zero basin_orog(ib,ij) = zero basin_floodp(ib,ij) = zero ! SELECT CASE (hierar_method) ! CASE("MINI") basin_hierarchy(ib,ij) = undef_sechiba ! END SELECT basin_topoind(ib,ij) = zero ! 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)) ! Numerator for mean orography, then need to be divided by the HTUs area basin_orog(ib,ij) = basin_orog(ib,ij) + orog_bx(basin_pts(ij,iz,1),basin_pts(ij,iz,2)) & & * area_bx(basin_pts(ij,iz,1),basin_pts(ij,iz,2)) basin_floodp(ib,ij) = basin_floodp(ib,ij) + floodp_bx(basin_pts(ij,iz,1),basin_pts(ij,iz,2)) ! ! Compute centre of gravity ! basin_cg(ib,ij,1) = basin_cg(ib,ij,1) + lat_bx(basin_pts(ij,iz,1),basin_pts(ij,iz,2))/basin_sz(ij) basin_cg(ib,ij,2) = basin_cg(ib,ij,2) + lon_bx(basin_pts(ij,iz,1),basin_pts(ij,iz,2))/basin_sz(ij) ! ! 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)) basin_fac(ib,ij) = basin_fac(ib,ij) + & & fac_bx(basin_pts(ij,iz,1),basin_pts(ij,iz,2)) CASE("MINI") ! The smallest 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)) ! Trung: We should take the "fac" value at the same point basin_fac(ib,ij) = fac_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 ) & & .OR. ( trip_bx(basin_pts(ij,iz,1),basin_pts(ij,iz,2)) .EQ. -4 ) ) THEN basin_hierarchy(ib,ij) = hierarchy_bx(basin_pts(ij,iz,1),basin_pts(ij,iz,2)) ENDIF ! Trung: Be careful of statement " .LT. 0" IF ( ( trip_bx(basin_pts(ij,iz,1),basin_pts(ij,iz,2)) .GE. 97 ) & & .OR. ( trip_bx(basin_pts(ij,iz,1),basin_pts(ij,iz,2)) .EQ. -4 ) ) THEN basin_fac(ib,ij) = fac_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' CALL ipslerr_p(3,'routing_reg_globalize','Unknown method for computing the hierarchy of the basin','','') END SELECT ! ENDDO ! ! Finish the calculation of mean orography basin_orog(ib,ij) = basin_orog(ib,ij)/basin_area(ib,ij) ! SELECT CASE (hierar_method) ! CASE("MEAN") basin_hierarchy(ib,ij) = basin_hierarchy(ib,ij)/REAL(basin_sz(ij),r_std) basin_fac(ib,ij) = basin_fac(ib,ij)/REAL(basin_sz(ij),r_std) ! END SELECT ! ! Calculating topoindex depends on options ! basin_topoind(ib,ij) = zero ! Constant value (1000.) IF ( option .EQ. 1 ) THEN basin_topoind(ib,ij) = 1000. ENDIF ! Simple average (as the old method) IF ( option .EQ. 2 ) THEN DO iz=1,basin_sz(ij) basin_topoind(ib,ij) = basin_topoind(ib,ij) + topoind_bx(basin_pts(ij,iz,1),basin_pts(ij,iz,2)) ENDDO basin_topoind(ib,ij) = basin_topoind(ib,ij)/REAL(basin_sz(ij),r_std) ENDIF ! Minus outlet hierarchy IF ( option .EQ. 3 ) THEN DO iz=1,basin_sz(ij) area_frac = area_bx(basin_pts(ij,iz,1),basin_pts(ij,iz,2))/REAL(basin_area(ib,ij)) basin_topoind(ib,ij) = basin_topoind(ib,ij) + topoind_bx(basin_pts(ij,iz,1),basin_pts(ij,iz,2))*area_frac ENDDO ! ! Trung: In routing_reg_globalize: ! Trung: Checking negative value of basin_topoind ! IF ( debug ) THEN IF ( basin_topoind(ib,ij) .LT. 0. ) THEN WRITE (numout,*) 'In routing_reg_globalize: We got negative value in ',ib,ij WRITE (numout,*) 'Start checking ... ' WRITE (numout,*) 'basin_area ', REAL(basin_area(ib,ij)) WRITE (numout,*) ' ==================> ' DO iz=1,basin_sz(ij) area_frac = area_bx(basin_pts(ij,iz,1),basin_pts(ij,iz,2))/REAL(basin_area(ib,ij)) WRITE (numout,*) 'iz = ',iz WRITE (numout,*) 'topoind_bx :', topoind_bx(basin_pts(ij,iz,1),basin_pts(ij,iz,2)) WRITE (numout,*) 'area_bx ', area_bx(basin_pts(ij,iz,1),basin_pts(ij,iz,2)) WRITE (numout,*) 'area_frac ', area_frac WRITE (numout,*) 'topoind_bx * area_frac ', topoind_bx(basin_pts(ij,iz,1),basin_pts(ij,iz,2))*area_frac ENDDO ENDIF ENDIF ENDIF ! ! 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 .AND. basin_bxout(ij) .NE. -4) THEN IF ( option .NE. 1 ) THEN basin_hierarchy(ib,ij) = min_topoind !basin_fac(ib,ij) = min_topoind ELSE basin_hierarchy(ib,ij) = 0.00 !basin_fac(ib,ij) = 0.00 ENDIF 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 basin_lshead(ib,ij) = lshead(ij) outflow_grid(ib,ij) = routing_nextgrid(ib,basin_bxout(ij)) ELSE basin_lshead(ib,ij) = undef_sechiba outflow_grid(ib,ij) = basin_bxout(ij) ENDIF ! outflow_basin(ib,ij) = undef_int IF (basin_bxout(ij) .EQ. -4) THEN outflow_basin(ib,ij) = basin_bbout(ij) ENDIF ! ! ENDDO ! ! Get the number of sub-grid basin before truncation ! basin_notrun = basin_count ! END SUBROUTINE routing_reg_globalize ! !! ================================================================================================================================ !! SUBROUTINE : routing_reg_linkup !! !>\BRIEF This subroutine makes the connections between the basins and ensure global coherence. !! !! DESCRIPTION (definitions, functional, design, flags) : !! The convention for outflow_grid is : !! outflow_grid > : The grid bow in which the basin is supposed to flow. !! outflow_grid = -1 : River flow !! outflow_grid = -2 : Coastal flow !! outflow_grid = -3 : Return flow !! outflow_grid = -4 : Flows into a basin in the same grid !! For outflow_basin we have the following conventions : !! outflow_basin = basin id in the case of the basin not flowing out of the current grid (i.e. outflow_grid=-4). !! Else outflow_basin = undef_int and the objective of this routine is to find what that basin is. This work !! essentially occurs in the routine routing_reg_bestsubbasin. But for that we need to find the grid box where !! we will look for the right basin. This operation is performed here in 5 successive steps. The order correspond !! to the solution we would prefer. !! 1.0 : We will look in the grid provided by outflow_grid for the right basin. This neighbour has been obtained !! from the small scale flow direction given in routing.nc. !! 2.0 : Try the grid box given by the large scale flow direction. !! 3.0 : Have another attempt at the outflow_grid by looking at the neighbour just to the right and the left. !! 4.0 : Here we look at half of the neighbouring grid boxes around the large scale flow direction. !! 5.0 : If all the above failed then we look within the current grid box if we find a suitable outflow basin. !! !! RECENT CHANGE(S): None !! !! MAIN OUTPUT VARIABLE(S): !! !! REFERENCES : None !! !! FLOWCHART : None !! \n !_ ================================================================================================================================ SUBROUTINE routing_reg_linkup(nbpt, neighbours, nwbas, basin_count, basin_area, basin_id, basin_flowdir, & & basin_lshead, basin_hierarchy, basin_fac, diaglalo, outflow_grid, outflow_basin, inflow_number, inflow_grid, & & inflow_basin, nbcoastal, coastal_basin, invented_basins) ! IMPLICIT NONE ! !! INPUT VARIABLES INTEGER(i_std), INTENT (in) :: nbpt !! Domain size (unitless) INTEGER(i_std), INTENT (in), DIMENSION(nbpt,NbNeighb) :: neighbours !! REAL(r_std), INTENT(in) :: invented_basins !! ! INTEGER(i_std), INTENT (in) :: nwbas !! INTEGER(i_std), INTENT (in), DIMENSION(nbpt) :: basin_count !! INTEGER(i_std), INTENT (in), DIMENSION(nbpt,nwbas) :: basin_id !! INTEGER(i_std), INTENT (in), DIMENSION(nbpt,nwbas) :: basin_flowdir !! Water flow directions in the basin (unitless) REAL(r_std), INTENT (in), DIMENSION(nbpt,nwbas) :: basin_lshead !! Large scale flow direction out of the basin. REAL(r_std), INTENT (in), DIMENSION(nbpt,nwbas) :: basin_area !! REAL(r_std), INTENT (inout), DIMENSION(nbpt,nwbas) :: basin_hierarchy !! REAL(r_std), INTENT (in), DIMENSION(nbpt,nwbas) :: basin_fac !! REAL(r_std), INTENT(in), DIMENSION(:,:) :: diaglalo !! Point (in Lat/Lon) where diagnostics will be printed. INTEGER(i_std), INTENT(inout), DIMENSION(nbpt,nwbas) :: outflow_grid !! Type of outflow on the grid box (unitless) INTEGER(i_std), INTENT(inout), DIMENSION(nbpt,nwbas) :: outflow_basin !! ! INTEGER(i_std), INTENT(out), DIMENSION(nbpt,nbxmax*8) :: inflow_number !! INTEGER(i_std), INTENT(out), DIMENSION(nbpt,nbxmax*8,nbxmax*8) :: inflow_basin !! INTEGER(i_std), INTENT(out), DIMENSION(nbpt,nbxmax*8,nbxmax*8) :: inflow_grid !! ! INTEGER(i_std), INTENT(inout), DIMENSION(nbpt) :: nbcoastal !! INTEGER(i_std), INTENT(inout), DIMENSION(nbpt,nwbas) :: coastal_basin !! ! !! LOCAL VARIABLES INTEGER(i_std) :: sp, sb, sbl, inp, outdm1, outdp1 !! Indices (unitless) INTEGER(i_std) :: spa, sba, sbint INTEGER(i_std) :: i, nbocean, it REAL(r_std) :: wocean INTEGER(i_std) :: dp1, dm1, dm1i, dls, bp1, bm1, bls !! Indices (unitless) INTEGER(i_std) :: dop, bop, nbtotest !! INTEGER(i_std) :: fbas(nwbas), nbfbas !! REAL(r_std) :: ang, angp1, angm1, bopqual, blsqual, bp1qual, bm1qual, crit REAL(r_std) :: fbas_hierarchy(nwbas) !! INTEGER(i_std) :: ff(1) !! INTEGER(i_std), DIMENSION(NbNeighb) :: gridstotest, gridbasin REAL(r_std), DIMENSION(NbNeighb) :: diffangle, gridangle INTEGER(i_std), DIMENSION(nbpt,5) :: solved INTEGER(i_std) :: unsolved, testbasinid REAL(r_std), ALLOCATABLE, DIMENSION(:) :: hatoutflow INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:) :: minhloc ! INTEGER(i_std) , DIMENSION(NbNeighb) :: order_ref INTEGER(i_std) :: nb_sym, nb, val, ind_nb, onb INTEGER(i_std) :: found ! !! PARAMETERS LOGICAL, PARAMETER :: debug = .TRUE. !! (true/false) ! !_ ================================================================================================================================ ! ! testbasinid = 195 ! IF ( debug ) WRITE (numout,*) 'SIZE inflow_grid:',SIZE(inflow_grid,1),SIZE(inflow_grid,2),SIZE(inflow_grid,3) IF ( debug ) WRITE (numout,*) 'SIZE inflow_basin:',SIZE(inflow_basin,1),SIZE(inflow_basin,2),SIZE(inflow_basin,3) CALL FLUSH(numout) ! inflow_number(:,:) = 0 ! ! Variable to keep track of how many points were solved at each stage ! solved(:,:) = 0 ! ! Some preparatory work. We need the minimum hierarchy for each basin and its location. ! This comes from the fact that for coastal basins it is often difficult to decide if they ! first flow into another grid or directly into the ocean. So this information will be ! useful. ! ALLOCATE(hatoutflow(INT(invented_basins))) ALLOCATE(minhloc(INT(invented_basins),2)) hatoutflow(:) = undef_sechiba DO sp=1,nbpt DO sb=1,basin_count(sp) IF ( basin_id(sp,sb) < INT(invented_basins) .AND. basin_id(sp,sb) > 0 ) THEN IF ( hatoutflow(basin_id(sp,sb)) > basin_hierarchy(sp,sb) ) THEN hatoutflow(basin_id(sp,sb)) = basin_hierarchy(sp,sb) minhloc(basin_id(sp,sb),1) = sp minhloc(basin_id(sp,sb),2) = sb ENDIF ENDIF ENDDO ENDDO ! ! 1.0 Follow the flow direction as given by the high resolution basin description. This case also treats the ! basins which remain in the current grid box, rivers, coastal and return flows. ! DO sp=1,nbpt ! DO sb=1,basin_count(sp) ! ! We only work on this point if it does not flow into the ocean ! or flow to another sub-basin in the same grid box. ! At this point any of the outflows is designated by a negative values in ! outflow_grid ! WRITE(numout,*) "*****" WRITE(numout,*) "Linkup 0 - sp, sb = ", sp, sb WRITE(numout,*) "Linkup 0 - outflow_grid = ", outflow_grid(sp,sb) IF ( outflow_grid(sp,sb) == 0 ) WRITE(numout,*) "Linkup 1.0 -- Flow out of Halo zone" ! IF ( outflow_grid(sp,sb) < 0 ) THEN IF ( outflow_grid(sp,sb) .EQ. -4 ) THEN ! Flow into a basin of the same grid bop = outflow_basin(sp,sb) CALL routing_updateflow(sp, sb, sp, bop, nbpt, nwbas, nbxmax*8, outflow_grid, outflow_basin, & & inflow_number, inflow_grid, inflow_basin) IF ( outflow_basin(sp,sb) == bop ) THEN solved(sp,1) = solved(sp,1) + 1 ENDIF WRITE(numout,*) sp,sb,"flow in the same grid, basin:",bop ! ELSE IF ( outflow_grid(sp,sb) .EQ. -3 ) THEN ! Return flow ! Nothing to do but just remember it is done. solved(sp,1) = solved(sp,1) + 1 WRITE(numout,*) sp,sb,"is a return flow" ELSE IF ( outflow_grid(sp,sb) .EQ. -2 ) THEN ! Coastal flow ! Nothing to do but just remember it is done. solved(sp,1) = solved(sp,1) + 1 WRITE(numout,*) sp,sb,"is a coastal flow" ELSE IF ( outflow_grid(sp,sb) .EQ. -1 ) THEN ! River flow ! Nothing to do but just remember it is done. solved(sp,1) = solved(sp,1) + 1 WRITE(numout,*) sp,sb,"is a river outflow" ENDIF ELSE IF ( outflow_grid(sp,sb) .GT. 0 ) THEN found = 0 ! ! Deal with the first one as usual ! inp = outflow_grid(sp,sb) ! CALL routing_reg_bestsubbasin(sp, sb, basin_id(sp,sb), basin_hierarchy(sp,sb), basin_fac(sp,sb), & & basin_flowdir(sp,sb), invented_basins, & & nbpt, nwbas, inp, basin_count, basin_id, basin_hierarchy, basin_fac, & & basin_flowdir, nbcoastal, coastal_basin, bop, bopqual) ! IF ( bop .LT. undef_int ) THEN ! CALL routing_updateflow(sp, sb, inp, bop, nbpt, nwbas, nbxmax*8, outflow_grid, outflow_basin, & & inflow_number, inflow_grid, inflow_basin) IF ( outflow_basin(sp,sb) == bop ) THEN solved(sp,1) = solved(sp,1) + 1 found = 1 WRITE(numout,*) "Solution found in the original outflow_grid" ENDIF ! ENDIF ! IF ( found == 0 ) THEN WRITE (numout,*) "Establishing the list of neighbours" ! Organize the location of the neighbours to visit by order of priority ! first the outflow grid then 2 by 2 till the opposite side (by +1/-1 - +2/-2 ...) ! if NbNeighb is odd we have to had the opposite neighbour ! ex : if NbNeighb = 8 and outflow grid is at the bottom right ! ! 8 | 7 | 5 ! ---|---|--- ! 6 | x | 3 ! ---|---|--- ! 4 | 2 | 1 ! order_ref(:) = -1 inp = outflow_grid(sp,sb) ! The first one is the location of the outflow grid order_ref(1) = minloc(abs(neighbours_g(sp,:) - inp), 1) ! Then we have to adjust in function of NbNeighb to add the location 2 by 2 nb_sym = (NbNeighb-1)/2 ind_nb = 2 DO nb=1,nb_sym DO onb=-1,1,2 val = modulo(order_ref(1)+nb*onb,NbNeighb) IF ( val == 0) val = NbNeighb order_ref(ind_nb) = val ind_nb = ind_nb + 1 ENDDO ENDDO ! If NbNeighb is odd IF ( modulo(NbNeighb,2) == 0 ) THEN val = modulo(order_ref(1)+4,NbNeighb) IF ( val == 0) val = NbNeighb order_ref(NbNeighb) = val ENDIF ! ! ! If not found, now test the next points 2 by 2 ! IF ( found == 0) THEN nb = 1 DO WHILE (( found == 0 ) .AND. (nb .LE. nb_sym)) ! ! From the grid point + and - get bp and bpqual ! First point : dpi1 -> bp1 et bp1qual dp1 = neighbours_g(sp,order_ref(nb*2)) ! we know it's wrong but it will serve to make decision later angp1 = routing_anglediff_g(sp, dp1, basin_flowdir(sp,sb)) ! Check if grid point IF (dp1 .GT. 0) THEN CALL routing_reg_bestsubbasin(sp, sb, basin_id(sp,sb), basin_hierarchy(sp,sb), basin_fac(sp,sb), & & basin_flowdir(sp,sb), invented_basins, & & nbpt, nwbas, dp1, basin_count, basin_id, basin_hierarchy, basin_fac, & & basin_flowdir, nbcoastal, coastal_basin, bp1, bp1qual) ELSE bp1 = undef_int bp1qual = 0 ENDIF ! ! Second point ! dm1 -> bm1 et bm1qual dm1 = neighbours_g(sp,order_ref(nb*2+1)) ! we know it's wrong but it will serve to make decision later angm1 = routing_anglediff_g(sp, dm1, basin_flowdir(sp,sb)) IF (dm1 .GT. 0) THEN CALL routing_reg_bestsubbasin(sp, sb, basin_id(sp,sb), basin_hierarchy(sp,sb), basin_fac(sp,sb), & & basin_flowdir(sp,sb), invented_basins, & & nbpt, nwbas, dm1, basin_count, basin_id, basin_hierarchy, basin_fac, & & basin_flowdir, nbcoastal, coastal_basin, bm1, bm1qual) ELSE bm1 = undef_int bm1qual = 0 ENDIF ! ! Decide between dp1 and dm1 which is the better solution ! dop = undef_int IF ( bp1 < undef_int .OR. bm1 < undef_int ) THEN IF ( bp1 < undef_int ) THEN dop = dp1 bop = bp1 ELSE IF ( bm1 < undef_int ) THEN dop = dm1 bop = bm1 ELSE IF ( bp1qual > bm1qual ) THEN dop = dp1 bop = bp1 ELSE IF ( bp1qual > bm1qual ) THEN dop = dm1 bop = bm1 ELSE ! ATTENTION !! ! In the case both point are the same quality ! we need of another factor of decision ! for now we use the angles even if we know it's wrong IF ( angm1 < angp1 ) THEN dop = dm1 bop = bm1 ELSE dop = dp1 bop = bp1 ENDIF ENDIF ENDIF ENDIF IF ( dop > 0 .AND. dop < undef_int .AND. bop < undef_int ) THEN ! CALL routing_updateflow(sp, sb, dop, bop, nbpt, nwbas, nbxmax*8, outflow_grid, outflow_basin, & & inflow_number, inflow_grid, inflow_basin) ! IF ( outflow_basin(sp,sb) == bop ) THEN solved(sp,2) = solved(sp,2) + 1 found = 1 WRITE (numout,*) "Neighbours, output found at:" , nb, "level" ENDIF ! ELSE nb = nb+1 ENDIF ENDDO ENDIF ! ! Should we check the opposite point if odd NeighNb? !IF ( modulo(NbNeighb,2) == 0 ) THEN ! ! Do like first one ! ! loc(NbNeighb) ENDIF ! ! Looking for a solution in the grid -> HTU with a similar hierarchy or lowest hierarchy ! IF ( found == 0 ) THEN WRITE (numout,*) "Looking for a solution in the grid" sbint = undef_int DO sba=1,basin_count(sp) IF ( sba .NE. sb ) THEN IF ( basin_id(sp,sb) == basin_id(sp,sba) ) THEN ! ! Only look for suitable hierrachy if we have a solution for the sub-basin : it has an outflow grid ! and basin or flows into the ocean. ! IF ( (outflow_grid(sp,sba) > 0 .AND. outflow_basin(sp,sba) < undef_int) .OR. & & outflow_grid(sp,sba) .EQ. -1 .OR. outflow_grid(sp,sba) .EQ. -2 ) THEN ! ! calculation of the criterion for considering that the hierarchy are close enough. ! Watch out hierarchy can be zero. 1 is considered the ! minimum here. ! ! Trung: If current sub-basin (sp,sb) is the closest ! point to the outlet of river, the value of ! basin_hierarchy(sp,sb) is always lower than ! basin_hierarchy(sp,sba). It means that crit ! is negative and statement IF (crit < 0.1) is ! not useful. => I changed it a little bit. ! crit = (basin_hierarchy(sp,sb)-basin_hierarchy(sp,sba))/(basin_hierarchy(sp,sb)+1) IF ( crit < 0.01 .AND. crit >= zero ) THEN ! Trung: Put condition of flow accumulation here to ! avoid error of circulation in routing_reg_fetch IF (basin_fac(sp,sba) .GE. basin_fac(sp,sb)) THEN sbint = sba ENDIF ENDIF ENDIF ENDIF ENDIF ENDDO ! ! If this grid contains the lowest hierarchy of this basin ID, we can also merge without second thought. ! IF ( basin_id(sp,sb) < invented_basins ) THEN IF ( sp == minhloc(basin_id(sp,sb),1) ) THEN sbint = minhloc(basin_id(sp,sb),2) ENDIF ENDIF IF (( sbint < undef_int ) .AND. (sbint .GT. 0) .AND. (sbint .NE. sb)) THEN basin_hierarchy(sp,sb) = basin_hierarchy(sp,sbint) CALL routing_updateflow(sp, sb, sp, sbint, nbpt, nwbas, nbxmax*8, outflow_grid, outflow_basin, & & inflow_number, inflow_grid,inflow_basin) found = 1 WRITE (numout,*) "Lowest basin hierarchy in the grid file" ENDIF ENDIF ! ! Looking if the minimun hierarchy is in a neighbour grid point IF (found .EQ. 0) THEN dop = undef_int bop = undef_int IF ( basin_id(sp,sb) < invented_basins ) THEN DO nb=1,NbNeighb IF (( neighbours_g(sp,order_ref(nb)) == minhloc(basin_id(sp,sb),1)) .AND. (found == 0)) THEN dop = minhloc(basin_id(sp,sb),1) bop = minhloc(basin_id(sp,sb),2) IF ((dop < undef_int ) .AND. (dop .GT. 0) .AND. (bop < undef_int) .AND. (bop .GT. 0)) THEN CALL routing_updateflow(sp, sb, dop, bop, nbpt,nwbas, nbxmax*8, outflow_grid, outflow_basin, & & inflow_number, inflow_grid, inflow_basin) ! It is possible that the lowest hierarchy is in two grid cells ! In this case, if we have an error for this routing_updateflow ! We go to next step and make it a coastal flow IF (outflow_basin(sp, sb) < undef_int) THEN found = 1 WRITE(numout,*) "Lowest basin hierarchy in the neighbours grid points" ELSE WRITE(numout,*) "Lowest hierarchy may be in two different grid points" END IF WRITE(numout,*) "hierarch of (sp,sb)", basin_hierarchy(sp,sb) WRITE(numout,*) "Lowest basinid hierarch", basin_hierarchy(dop,bop) END IF END IF END DO END IF ENDIF ! ! Last option : coastal outflow IF (found .EQ. 0) THEN ! ! Trung: Well, if you come to this step, it should be a special ! (or curious, I said) case. For example, when processing ! Sulina branch of river Danube with the forcing data ! from E2OFD (/bdd/ORCHIDEE_Forcing/BC/OOL/OL2/E2OFD/) ! in resolution of 0.25 degree. The grid box [45.125N, 29.375E] ! is a ocean point in the middle of 7 land points ("a ! hole"). But HydroSHEDS shows that Sulina branch flows ! through this grid cell (from [45.125N, 29.125E] to [45.125N, 29.625E]). ! Trung: If the "JUMP" does not work, it means a worse case. ! Or it means a worse case. ! The outlet of this river falls in the ocean point of ! forcing data. For example: with 0.25 degree forcing ! data of E2OFD, the grid point [ 45.625N ; 13.875E ] ! lost the outlet of basin ID 69666. ! We need to use this point (which should be the last point ! of this river with shortest hierarchy and highest flow accumulation) ! to become new outlet. ! ! This is not good solution because there can be few ! outlet for one basin ID. ! But I need the model works now !!! => so, come back ! here later ! ! WRITE (numout,*) "Linkup : Made a NEW OUTLET at sp & sb: ",sp,sb ! Coastal flow or river flow is both ok here outflow_grid(sp,sb) = -2 basin_hierarchy(sp,sb) = 0.00 ENDIF IF ( outflow_basin(sp,sb) < undef_int ) solved(sp,4) = solved(sp,4)+1 ENDIF ENDDO ENDDO ! END SUBROUTINE routing_reg_linkup ! ! !! ================================================================================================================================ !! SUBROUTINE : routing_reg_fetch !! !>\BRIEF This subroutine computes 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. !! !! DESCRIPTION (definitions, functional, design, flags) : None !! !! RECENT CHANGE(S): None !! !! MAIN OUTPUT VARIABLE(S): !! !! REFERENCES : None !! !! FLOWCHART : None !! \n !_ ================================================================================================================================ SUBROUTINE routing_reg_fetch(nbpt, gridarea, contfrac, nwbas, nbcore, corepts, basin_count, basin_area, basin_id, & & basin_hierarchy, basin_fac, outflow_grid, outflow_basin, partial_sum, fetch_basin, outflow_uparea) ! IMPLICIT NONE ! !! INPUT VARIABLES INTEGER(i_std), INTENT(in) :: nbpt !! Domain size (unitless) ! REAL(r_std), DIMENSION(nbpt), INTENT(in) :: gridarea !! The size of each grid box in X and Y (m) REAL(r_std), DIMENSION(nbpt), INTENT(in) :: contfrac !! Fraction of land in each grid box (unitless;0-1) ! INTEGER(i_std) :: nwbas !! INTEGER(i_std), INTENT (in) :: nbcore INTEGER(i_std), DIMENSION(nbcore), INTENT(in) :: corepts INTEGER(i_std), DIMENSION(nbpt), INTENT(in) :: basin_count !! REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(in) :: basin_area !! INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(in) :: basin_id !! REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(in) :: basin_hierarchy REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(in) :: basin_fac INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(in) :: outflow_grid !! Type of outflow on the grid box (unitless) INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(in) :: outflow_basin !! REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(in) :: partial_sum ! !! INPUT/OUTPUT VARIABLES REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: fetch_basin !! REAL(r_std), DIMENSION(nbpt*nwbas), INTENT(out) :: outflow_uparea ! !! LOCAL VARIABLES INTEGER(i_std) :: ig, ib, ic, ff(1), it, itt, igrif, ibasf, nboutflow !! Indices (unitless) REAL(r_std), DIMENSION(nbpt*nbvmax) :: tmp_area !! INTEGER(i_std), DIMENSION(nbpt*nbvmax,2) :: tmpindex !! ! INTEGER(i_std), PARAMETER :: maxiterations=10000 ! !_ ================================================================================================================================ ! ! Propagate first the partial sums from the halo into the domain ! DO ig=1,nbpt ! DO ib=1,basin_count(ig) ! IF (partial_sum(ig,ib) > EPSILON(partial_sum)) THEN ! fetch_basin(ig, ib) = MAX(fetch_basin(ig, ib), partial_sum(ig,ib)) ! igrif = outflow_grid(ig,ib) ibasf = outflow_basin(ig,ib) ! itt = 0 ! DO WHILE (igrif .GT. 0) ! fetch_basin(igrif,ibasf) = MAX(fetch_basin(igrif,ibasf), partial_sum(ig, ib)) ! it = outflow_grid(igrif, ibasf) ibasf = outflow_basin(igrif, ibasf) igrif = it itt = itt + 1 IF ( itt .GT. maxiterations) THEN IF ( itt .GT. maxiterations+20) THEN CALL ipslerr_p(3,'routing_reg_fetch','Problem in propagating partial sum','','') ENDIF ENDIF ENDDO ENDIF ENDDO ENDDO ! ! Add the areas contributed by the core region of the domain. ! DO ic=1,nbcore ig = corepts(ic) ! DO ib=1,basin_count(ig) ! fetch_basin(ig,ib) = fetch_basin(ig,ib) + basin_area(ig,ib) ! igrif = outflow_grid(ig,ib) ibasf = outflow_basin(ig,ib) ! itt = 0 ! DO WHILE (igrif .GT. 0 ) ! fetch_basin(igrif,ibasf) = fetch_basin(igrif,ibasf) + basin_area(ig, ib) ! it = outflow_grid(igrif, ibasf) ibasf = outflow_basin(igrif, ibasf) igrif = it itt = itt + 1 IF ( itt .GT. maxiterations) THEN WRITE(numout,& "('Grid ',I5, ' and basin ',I5, ' did not converge after iteration ',I5)") ig, ib, itt WRITE(numout, "('==> We flow through grid ',I5,' and basin ',I5)") igrif, ibasf WRITE(numout,*) '==> Basin ID :', basin_id(igrif,ibasf), "Hierarchy :", basin_hierarchy(igrif,ibasf) WRITE(numout,*) "==> Flow accumulation :", basin_fac(igrif,ibasf) WRITE(numout,*) "==> Coordinates : ", lalo_g(igrif,2), lalo_g(igrif,1) ! IF ( itt .GT. maxiterations+20) THEN CALL ipslerr_p(3,'routing_reg_fetch','Problem...','','') ENDIF ENDIF ENDDO ! ENDDO ! ENDDO ! WRITE(numout,*) 'The smallest FETCH :', MINVAL(fetch_basin) WRITE(numout,*) 'The largest FETCH :', MAXVAL(fetch_basin) ! ! nboutflow = 0 outflow_uparea(:) = zero ! DO ic=1,nbcore ig = corepts(ic) ! DO ib=1,basin_count(ig) ! ! We do not need any more the river flow flag as we are going to reset it. ! We archive the upstream area of all outflow points so that we can sort them. ! IF ( outflow_grid(ig,ib) .LT. 0) THEN nboutflow = nboutflow + 1 IF ( nboutflow <= nbpt*nwbas) THEN outflow_uparea(nboutflow) = fetch_basin(ig,ib) ELSE ! Not enought place so we replace a smaller one. IF (fetch_basin(ig,ib) > MINVAL(outflow_uparea)) THEN ff = MINLOC(outflow_uparea) outflow_uparea(ff(1)) = fetch_basin(ig,ib) ELSE ! Ignore basin ENDIF ENDIF ENDIF ! ENDDO ENDDO ! END SUBROUTINE routing_reg_fetch ! ! !! ================================================================================================================================ !! SUBROUTINE : routing_reg_truncate !! !>\BRIEF This subroutine reduces the number of basins per grid 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. !! !! DESCRIPTION (definitions, functional, design, flags) : !! 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.:\n !! !! RECENT CHANGE(S): None !! !! MAIN OUTPUT VARIABLE(S): !! !! REFERENCES : None !! !! FLOWCHART : None !! \n !_ ================================================================================================================================ SUBROUTINE routing_reg_end_truncate(nbpt, nbasmax, gridarea, contfrac, nwbas, num_largest, & & basin_count, basin_notrun, basin_area, basin_orog, basin_floodp, & & basin_cg, basin_topoind, fetch_basin, basin_id, basin_coor, basin_type, basin_flowdir, & & outflow_grid, outflow_basin, inflow_number, inflow_grid, inflow_basin) ! IMPLICIT NONE ! !! PARAMETERS ! INTEGER(i_std), PARAMETER :: pickmax = 2000 !! !! INPUT VARIABLES INTEGER(i_std), INTENT(in) :: nbpt !! Domain size (unitless) INTEGER(i_std), INTENT(in) :: nbasmax ! REAL(r_std), DIMENSION(nbpt), INTENT(in) :: gridarea !! The size of each grid box in X and Y (m) REAL(r_std), DIMENSION(nbpt), INTENT(in) :: contfrac !! Fraction of land in each grid box (unitless;0-1) ! INTEGER(i_std), INTENT(in) :: nwbas !! INTEGER(i_std), INTENT (in) :: num_largest INTEGER(i_std), DIMENSION(nbpt), INTENT(inout) :: basin_count !! INTEGER(i_std), DIMENSION(nbpt), INTENT(inout) :: basin_notrun !! INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: basin_id !! REAL(r_std), DIMENSION(nbpt,nwbas,2), INTENT(inout) :: basin_coor !! REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: basin_type !! INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: basin_flowdir !! Water flow directions in the basin (unitless) REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: basin_area !! REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: basin_orog !! REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: basin_floodp !! REAL(r_std), DIMENSION(nbpt,nwbas,2), INTENT(inout) :: basin_cg !! REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: basin_topoind !! Topographic index of the residence time for a basin (m) REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: fetch_basin !! INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: outflow_grid !! Type of outflow on the grid box (unitless) INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: outflow_basin !! ! ! Changed nwbas to nbxmax*4 and *8 to save the memory ! INTEGER(i_std), DIMENSION(nbpt,nbxmax*8), INTENT(inout) :: inflow_number !! INTEGER(i_std), DIMENSION(nbpt,nbxmax*8,nbxmax*8), INTENT(inout) :: inflow_basin !! INTEGER(i_std), DIMENSION(nbpt,nbxmax*8,nbxmax*8), INTENT(inout) :: inflow_grid !! ! !! LOCAL VARIABLES INTEGER(i_std) :: ib, ij, ic, ibf, ijf, igrif, ibasf, cnt, pold, bold, ff(2) !! Indices (unitless) INTEGER(i_std) :: ii, kbas, sbas, ik, iter, ibt, obj !! Indices (unitless) REAL(r_std), DIMENSION(nbpt,nbasmax) :: floflo !! REAL(r_std), DIMENSION(nbpt) :: gridarea_tmp !! REAL(r_std), DIMENSION(nbpt) :: gridbasinarea !! REAL(r_std) :: ratio !! INTEGER(i_std), DIMENSION(pickmax,2) :: largest_basins !! INTEGER(i_std), DIMENSION(pickmax) :: tmp_ids !! INTEGER(i_std) :: multbas !! INTEGER(i_std) :: iml(1) !! X resolution of the high resolution grid 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 !! !$OMP THREADPRIVATE(indextrunc) ! LOGICAL :: trunc_equa = .TRUE. !! Logical to choose if equality truncation is activated or not (true/false) ! !_ ================================================================================================================================ ! ! ALLOCATE memory ! CALL allocategraph(nbpt, nbasmax) ! ! Read option for truncation method ! ! 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_glo(ib,:) = ib route_tobasin_glo(ib,:) = 0 routing_area_glo(ib,:) = zero route_fetch_glo(ib,:) = zero route_count_glo(ib) = basin_count(ib) route_nbintobas_glo(ib) = basin_count(ib) origin_nbintobas_glo(ib) = basin_notrun(ib) ! ENDDO ! ! Transfer the info into the definitive variables ! DO ib=1,nbpt DO ij=1,basin_count(ib) ! routing_area_glo(ib,ij) = basin_area(ib,ij) routing_orog_glo(ib,ij) = basin_orog(ib,ij) routing_floodp_glo(ib,ij) = basin_floodp(ib,ij) routing_cg_glo(ib,ij,:) = basin_cg(ib,ij,:) topo_resid_glo(ib,ij) = basin_topoind(ib,ij) global_basinid_glo(ib,ij) = basin_id(ib,ij) route_outlet_glo(ib,ij,1) = basin_coor(ib,ij,1) route_outlet_glo(ib,ij,2) = basin_coor(ib,ij,2) ! route_fetch_glo(ib,ij) = fetch_basin(ib,ij) ! route_type_glo(ib,ij) = basin_type(ib,ij) ! route_togrid_glo(ib,ij) = outflow_grid(ib,ij) route_tobasin_glo(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) ! Flowing out of the grid IF ( route_togrid_glo(ib,ij) .EQ. 0 ) THEN route_tobasin_glo(ib,ij) = nbasmax + 2 route_togrid_glo(ib,ij) = ib ! River flows ELSE IF ( route_togrid_glo(ib,ij) .EQ. -1 ) THEN route_tobasin_glo(ib,ij) = nbasmax + 2 route_togrid_glo(ib,ij) = ib ! Coastal flows ELSE IF ( route_togrid_glo(ib,ij) .EQ. -2 ) THEN route_tobasin_glo(ib,ij) = nbasmax + 2 route_togrid_glo(ib,ij) = ib ! Return flow ELSE IF ( route_togrid_glo(ib,ij) .EQ. -3 ) THEN route_tobasin_glo(ib,ij) = nbasmax + 1 route_togrid_glo(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_glo(ib,ij) ijf = route_tobasin_glo(ib,ij) IF ( ibf .GT. 0 ) THEN 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) CALL ipslerr_p(3,'routing_reg_truncate','Problem with routing..','','') ENDIF ENDIF ENDDO ENDDO ! ! Verify areas of the continents ! floflo(:,:) = zero gridarea_tmp(:) = contfrac(:)*gridarea(:) DO ib=1,nbpt gridbasinarea(ib) = SUM(routing_area_glo(ib,:)) ENDDO ! DO ib=1,nbpt IF ( gridbasinarea(ib) > zero ) THEN ratio = gridarea_tmp(ib)/gridbasinarea(ib) routing_area_glo(ib,:) = routing_area_glo(ib,:)*ratio ELSE WRITE(numout,*) 'gridbasinarea(ib) <= zero. We should stop here :', ib ENDIF ENDDO ! WRITE(numout,*) 'Sum of area of all outflow areas :',SUM(routing_area_glo) WRITE(numout,*) 'Surface of all continents :', SUM(gridarea_tmp) ! END SUBROUTINE routing_reg_end_truncate ! !! ================================================================================================================================ !! SUBROUTINE : routing_reg_killbas !! !>\BRIEF The aim of this subroutine 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. !! !! DESCRIPTION (definitions, functional, design, flags) : None !! !! RECENT CHANGE(S): None !! !! MAIN OUTPUT VARIABLE(S): !! !! REFERENCES : None !! !! FLOWCHART : None !! \n !_ ================================================================================================================================ SUBROUTINE routing_reg_killbas(nbpt, ib, tokill, totakeover, nwbas, basin_count, basin_area, & & basin_orog, basin_floodp, basin_cg, basin_topoind,& & fetch_basin, basin_id, basin_coor, basin_type, basin_flowdir, outflow_grid, outflow_basin, inflow_number,& & inflow_grid, inflow_basin) ! ! IMPLICIT NONE ! INTEGER(i_std) :: tokill !! INTEGER(i_std) :: totakeover !! INTEGER(i_std) :: nbpt !! Domain size (unitless) INTEGER(i_std) :: ib !! Current basin (unitless) ! INTEGER(i_std) :: nwbas !! INTEGER(i_std), DIMENSION(nbpt) :: basin_count !! INTEGER(i_std), DIMENSION(nbpt,nwbas) :: basin_id !! REAL(r_std), DIMENSION(nbpt,nwbas,2) :: basin_coor !! REAL(r_std), DIMENSION(nbpt,nwbas) :: basin_type !! INTEGER(i_std), DIMENSION(nbpt,nwbas) :: basin_flowdir !! Water flow directions in the basin (unitless) REAL(r_std), DIMENSION(nbpt,nwbas) :: basin_area !! REAL(r_std), DIMENSION(nbpt,nwbas) :: basin_orog !! REAL(r_std), DIMENSION(nbpt,nwbas) :: basin_floodp !! REAL(r_std), DIMENSION(nbpt,nwbas,2) :: basin_cg !! REAL(r_std), DIMENSION(nbpt,nwbas) :: basin_topoind !! Topographic index of the residence time for a basin (m) REAL(r_std), DIMENSION(nbpt,nwbas) :: fetch_basin !! INTEGER(i_std), DIMENSION(nbpt,nwbas) :: outflow_grid !! Type of outflow on the grid box (unitless) INTEGER(i_std), DIMENSION(nbpt,nwbas) :: outflow_basin !! ! ! Saving memory: changed nwbas to nbxmax*4 and *8 ! INTEGER(i_std), DIMENSION(nbpt,nbxmax*8) :: inflow_number !! INTEGER(i_std), DIMENSION(nbpt,nbxmax*8,nbxmax*8) :: inflow_basin !! INTEGER(i_std), DIMENSION(nbpt,nbxmax*8,nbxmax*8) :: inflow_grid !! ! !! LOCAL VARIABLES INTEGER(i_std) :: inf, ibs, ing, inb, ibasf, igrif, it !! Indices (unitless) LOGICAL :: doshift !! (true/false) !_ ================================================================================================================================ ! ! Update the information needed in the basin "totakeover" ! For the moment only area ! basin_cg(ib, totakeover, 1) = (basin_cg(ib, totakeover, 1)*basin_area(ib, totakeover) & & + basin_cg(ib, tokill, 1)*basin_area(ib, tokill))/(basin_area(ib, totakeover) + basin_area(ib, tokill)) ! basin_cg(ib, totakeover, 2) = (basin_cg(ib, totakeover, 2)*basin_area(ib, totakeover)& & + basin_cg(ib, tokill, 2)*basin_area(ib, tokill))/(basin_area(ib, totakeover) + basin_area(ib, tokill)) ! basin_topoind(ib, totakeover) = (basin_topoind(ib, totakeover)*basin_area(ib, totakeover) & & + basin_topoind(ib, tokill)*basin_area(ib, tokill))/(basin_area(ib, totakeover) + basin_area(ib, tokill)) ! basin_area(ib, totakeover) = basin_area(ib, totakeover) + basin_area(ib,tokill) ! basin_floodp(ib, totakeover) = basin_floodp(ib, totakeover) + basin_floodp(ib,tokill) ! basin_orog(ib, totakeover) = (basin_orog(ib, totakeover)*basin_area(ib, totakeover) & & + basin_orog(ib, tokill)*basin_area(ib, tokill))/(basin_area(ib, totakeover) + basin_area(ib, tokill)) ! ! 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' CALL ipslerr_p(3,'routing_reg_killbas','Basin not found','','') 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_coor(ib, tokill:basin_count(ib)-1,1) = basin_coor(ib, tokill+1:basin_count(ib),1) basin_coor(ib, tokill:basin_count(ib)-1,2) = basin_coor(ib, tokill+1:basin_count(ib),2) basin_type(ib, tokill:basin_count(ib)-1) = basin_type(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) = zero basin_cg(ib, tokill:basin_count(ib)-1,:) = basin_cg(ib, tokill+1:basin_count(ib),:) basin_cg(ib, basin_count(ib):nwbas,:) = zero basin_topoind(ib, tokill:basin_count(ib)-1) = basin_topoind(ib, tokill+1:basin_count(ib)) basin_topoind(ib, basin_count(ib):nwbas) = zero fetch_basin(ib, tokill:basin_count(ib)-1) = fetch_basin(ib, tokill+1:basin_count(ib)) fetch_basin(ib, basin_count(ib):nwbas) = zero ! ! 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_reg_killbas ! SUBROUTINE allocategraph(nbpt, nbasmax) INTEGER(i_std), INTENT(in) :: nbpt, nbasmax INTEGER(i_std) :: ier nbpt_save = nbpt nbasmax_save = nbasmax ALLOCATE (routing_area_glo(nbpt,nbasmax), stat=ier) ALLOCATE (routing_orog_glo(nbpt,nbasmax), stat=ier) ALLOCATE (routing_floodp_glo(nbpt,nbasmax), stat=ier) ALLOCATE (route_fetch_glo(nbpt,nbasmax), stat=ier) ALLOCATE (routing_cg_glo(nbpt,nbasmax,2), stat=ier) ALLOCATE (topo_resid_glo(nbpt,nbasmax), stat=ier) ALLOCATE (route_count_glo(nbpt), stat=ier) ALLOCATE (route_togrid_glo(nbpt,nbasmax), stat=ier) ALLOCATE (route_tobasin_glo(nbpt,nbasmax), stat=ier) ALLOCATE (route_nbintobas_glo(nbpt), stat=ier) ALLOCATE (global_basinid_glo(nbpt,nbasmax), stat=ier) ALLOCATE (route_outlet_glo(nbpt,nbasmax,2), stat=ier) ALLOCATE (route_type_glo(nbpt,nbasmax), stat=ier) ALLOCATE (origin_nbintobas_glo(nbpt), stat=ier) END SUBROUTINE allocategraph END MODULE routing_reg