SUBROUTINE initatmgrid(rank, nbcore, nbpt, nbseg, polygons_in, centre_in, area_in, contfrac_in, neighbours_in) ! USE grid USE ioipsl ! IMPLICIT NONE ! !! INPUT VARIABLES INTEGER, INTENT(in) :: rank !! rank of processor INTEGER, INTENT(in) :: nbcore INTEGER, INTENT(in) :: nbpt !! Atmospheric Domain size (unitless) INTEGER, INTENT(in) :: nbseg !! Number of segments in the polygone REAL, INTENT(in) :: polygons_in(nbpt,2*nbseg+1,2) REAL, INTENT(in) :: centre_in(nbpt,2) REAL, INTENT(in) :: contfrac_in(nbpt) REAL, INTENT(in) :: area_in(nbpt) INTEGER, INTENT(in) :: neighbours_in(nbpt,nbseg*2) ! !! Local INTEGER :: i INTEGER :: neighbours_loc(nbpt,nbseg*2) CHARACTER(LEN=33) :: outfname ! DO i=1,nbpt ! Move from C to F indexing. neighbours_loc(i,:) = neighbours_in(i,:)+1 ENDDO ! CALL grid_init(nbpt, nbseg, polygons_in, centre_in, area_in, contfrac_in, neighbours_loc) ! numout=100+rank WRITE(outfname,'(A15,I4.4,A6,I4.4,A4)') "Out_RoutingReg_",rank,"_from_",nbcore,".txt" OPEN(unit=numout, file=outfname) ! END SUBROUTINE initatmgrid SUBROUTINE closeinterface() ! USE ioipsl ! CLOSE(unit=numout) ! END SUBROUTINE closeinterface SUBROUTINE gethydrogrid(nbpt, nbvmax_in, nbxmax_in, sub_pts, sub_index, sub_area, max_basins, min_topoind, & & lon_rel, lat_rel, trip, basins, topoindex, fac, hierarchy, nbi, nbj, area_bx, trip_bx, basin_bx, & & topoind_bx, fac_bx, hierarchy_bx, lon_bx, lat_bx, lshead_bx) ! USE ioipsl USE grid USE routing_tools USE routing_reg ! IMPLICIT NONE ! INTEGER, INTENT(in) :: nbvmax_in, nbxmax_in ! INTEGER, INTENT(in) :: nbpt !! Domain size (unitless) INTEGER, INTENT(in) :: sub_pts(nbpt) !! Number of high resolution points on this grid (unitless) INTEGER, INTENT(in) :: sub_index(nbpt,nbvmax_in,2) !! Indices of the points we need on the fine grid (unitless) REAL, INTENT(inout) :: max_basins !! The current maximum of basins REAL, INTENT(in) :: min_topoind !! The current minimum of topographic index (m) REAL, INTENT(in) :: sub_area(nbpt,nbvmax_in) !! Area on the fine grid (m^2) ! REAL, INTENT(in) :: lon_rel(nbpt,nbvmax_in) !! REAL, INTENT(in) :: lat_rel(nbpt,nbvmax_in) !! coordinates of the fine grid ! REAL, INTENT(inout) :: trip(nbpt,nbvmax_in) !! The trip field (unitless) REAL, INTENT(inout) :: basins(nbpt,nbvmax_in) !! data on the fine grid REAL, INTENT(inout) :: topoindex(nbpt,nbvmax_in) !! Topographic index of the residence time (m) REAL, INTENT(inout) :: hierarchy(nbpt,nbvmax_in) !! data on the fine grid REAL, INTENT(inout) :: fac(nbpt,nbvmax_in) !! data on the fine grid ! !! OUTPUT VARIABLES INTEGER, INTENT(out) :: nbi(nbpt), nbj(nbpt) !! Number of point in x and y within the grid (unitless) REAL, INTENT(out) :: area_bx(nbpt,nbxmax_in,nbxmax_in) !! Area of each small box in the grid box (m^2) REAL, INTENT(out) :: hierarchy_bx(nbpt,nbxmax_in,nbxmax_in) !! Level in the basin of the point REAL, INTENT(out) :: fac_bx(nbpt,nbxmax_in,nbxmax_in) !! Flow accumulation REAL, INTENT(out) :: lon_bx(nbpt,nbxmax_in,nbxmax_in) !! REAL, INTENT(out) :: lat_bx(nbpt,nbxmax_in,nbxmax_in) !! REAL, INTENT(out) :: lshead_bx(nbpt,nbxmax_in,nbxmax_in) !! Large scale heading for outflow points. REAL, INTENT(out) :: topoind_bx(nbpt,nbxmax_in,nbxmax_in) !! Topographic index of the residence time for each of the smaller boxes (m) INTEGER, INTENT(out) :: trip_bx(nbpt,nbxmax_in,nbxmax_in) !! The trip field for each of the smaller boxes (unitless) INTEGER, INTENT(out) :: basin_bx(nbpt,nbxmax_in,nbxmax_in) !! ! INTEGER :: ii, ib REAL :: resolution(nbpt,2) ! ! These two parameters are part of the routing_reg module saved variables. They are transfered here. ! nbvmax = nbvmax_in nbxmax = nbxmax_in ! DO ii=1,nbpt resolution(ii,1) = SQRT(area(ii)) resolution(ii,2) = SQRT(area(ii)) ENDDO DO ib=1,nbpt CALL 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, & & nbi(ib), nbj(ib), area_bx(ib,:,:), trip_bx(ib,:,:), basin_bx(ib,:,:), topoind_bx(ib,:,:), & & fac_bx(ib,:,:), hierarchy_bx(ib,:,:), lon_bx(ib,:,:), lat_bx(ib,:,:), lshead_bx(ib,:,:)) ENDDO ! END SUBROUTINE gethydrogrid SUBROUTINE findbasins(nbpt, nbvmax_in, nbxmax_in, nbi, nbj, trip_bx, basin_bx, fac_bx, hierarchy_bx, topoind_bx, lshead_bx, & & nb_basin, basin_inbxid, basin_outlet, basin_outtp, basin_sz, basin_bxout, basin_bbout, basin_pts, & & basin_lshead, coast_pts, lontmp, lattmp) ! USE ioipsl USE grid USE routing_tools USE routing_reg ! ! Arguments ! INTEGER, INTENT(in) :: nbvmax_in, nbxmax_in INTEGER, INTENT(in) :: nbpt !! Domain size (unitless) INTEGER, INTENT(in) :: nbi(nbpt), nbj(nbpt) !! Number of point in x and y within the grid (unitless) INTEGER, INTENT(inout) :: trip_bx(nbpt,nbxmax_in,nbxmax_in) !! The trip field for each of the smaller boxes (unitless) INTEGER, INTENT(inout) :: basin_bx(nbpt,nbxmax_in,nbxmax_in) !! REAL, INTENT(in) :: fac_bx(nbpt,nbxmax_in,nbxmax_in) !! Flow accumulation REAL, INTENT(in) :: hierarchy_bx(nbpt,nbxmax_in,nbxmax_in) !! Level in the basin of the point REAL, INTENT(inout) :: topoind_bx(nbpt,nbxmax_in,nbxmax_in) !! Topographic index of the residence time for each of the smaller boxes (m) REAL, INTENT(in) :: lshead_bx(nbpt,nbxmax_in,nbxmax_in) !! Large scale heading for outflow points. ! ! ! OUTPUT VARIABLES INTEGER, INTENT(out) :: nb_basin(nbpt) !! Number of sub-basins (unitless) INTEGER, INTENT(out) :: basin_inbxid(nbpt,nbvmax_in) !! REAL, INTENT(out) :: basin_outlet(nbpt,nbvmax_in,2) !! REAL, INTENT(out) :: basin_outtp(nbpt,nbvmax_in) !! INTEGER, INTENT(out) :: basin_sz(nbpt,nbvmax_in) !! INTEGER, INTENT(out) :: basin_bxout(nbpt,nbvmax_in) !! INTEGER, INTENT(out) :: basin_bbout(nbpt,nbvmax_in) !! INTEGER, INTENT(out) :: basin_pts(nbpt,nbvmax_in, nbvmax_in, 2) !! REAL, INTENT(out) :: basin_lshead(nbpt,nbvmax_in) !! INTEGER, INTENT(out) :: coast_pts(nbpt,nbvmax_in) !! The coastal flow points (unitless) ! ! For debug and get coordinate of river outlet ! REAL, INTENT(in) :: lontmp(nbpt,nbxmax_in,nbxmax_in) !! Longitude REAL, INTENT(in) :: lattmp(nbpt,nbxmax_in,nbxmax_in) !! Latitude ! ! Local ! INTEGER :: ib ! diaglalo(1,:) = (/ 39.6791, 2.6687 /) ! IF ( nbvmax_in .NE. nbvmax .OR. nbxmax_in .NE. nbxmax ) THEN WRITE(*,*) "nbvmax or nbxmax have changed !!" STOP ENDIF DO ib=1,nbpt CALL routing_reg_findbasins(ib, nbi(ib), nbj(ib), trip_bx(ib,:,:), basin_bx(ib,:,:), fac_bx(ib,:,:), hierarchy_bx(ib,:,:), & & topoind_bx(ib,:,:), lshead_bx(ib,:,:), diaglalo, & & nb_basin(ib), basin_inbxid(ib,:), basin_outlet(ib,:,:), basin_outtp(ib,:), basin_sz(ib,:), basin_bxout(ib,:), & & basin_bbout(ib,:), basin_pts(ib,:,:,:), basin_lshead(ib,:), coast_pts(ib,:), lontmp(ib,:,:), lattmp(ib,:,:)) ENDDO END SUBROUTINE findbasins SUBROUTINE globalize(nbpt, nbvmax_in, nbxmax_in, area_bx, lon_bx, lat_bx, trip_bx, hierarchy_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_fac, basin_topoind, basin_id, basin_coor, basin_type, basin_flowdir, basin_lshead, & & outflow_grid, outflow_basin, nbcoastal, coastal_basin) ! USE ioipsl USE grid USE routing_tools USE routing_reg ! ! Arguments ! !! INPUT VARIABLES INTEGER(i_std), INTENT (in) :: nbpt !! Domain size (unitless) INTEGER(i_std), INTENT(in) :: nbvmax_in, nbxmax_in REAL(r_std), INTENT (in), DIMENSION(nbpt,nbxmax_in,nbxmax_in) :: area_bx !! Area of each small box in the grid box (m^2) REAL(r_std), INTENT (in), DIMENSION(nbpt,nbxmax_in,nbxmax_in) :: lon_bx !! Longitude of each small box in the grid box REAL(r_std), INTENT (in), DIMENSION(nbpt,nbxmax_in,nbxmax_in) :: lat_bx !! Latitude of each small box in the grid box INTEGER(i_std), INTENT (in), DIMENSION(nbpt,nbxmax_in,nbxmax_in) :: trip_bx !! The trip field for each of the smaller boxes (unitless) REAL(r_std), INTENT (in), DIMENSION(nbpt,nbxmax_in,nbxmax_in) :: hierarchy_bx !! Level in the basin of the point REAL(r_std), INTENT (in), DIMENSION(nbpt,nbxmax_in,nbxmax_in) :: fac_bx !! Level in the basin of the point REAL(r_std), INTENT (in), DIMENSION(nbpt,nbxmax_in,nbxmax_in) :: topoind_bx !! Topographic index of the residence time for each of the smaller boxes (m) REAL(r_std), INTENT (in) :: min_topoind !! The current minimum of topographic index (m) INTEGER(i_std), INTENT (in), DIMENSION(nbpt) :: nb_basin !! Number of sub-basins (unitless) INTEGER(i_std), INTENT (in), DIMENSION(nbpt,nbvmax_in) :: basin_inbxid, basin_sz !! ID of basin, number of points in the basin REAL(r_std), INTENT (in), DIMENSION(nbpt,nbvmax_in,2) :: basin_outlet !! Outlet coordinate of each subgrid basin REAL(r_std), INTENT (in), DIMENSION(nbpt,nbvmax_in) :: basin_outtp !! INTEGER(i_std), INTENT (in), DIMENSION(nbpt,nbvmax_in,nbvmax_in,2) :: basin_pts !! Points in each basin INTEGER(i_std), INTENT (in), DIMENSION(nbpt,nbvmax_in) :: basin_bxout !! outflow direction INTEGER(i_std), INTENT (in), DIMENSION(nbpt,nbvmax_in) :: basin_bbout !! outflow sub-basin REAL(r_std), INTENT (in), DIMENSION(nbpt,nbvmax_in) :: lshead !! Large scale heading of outflow. INTEGER(i_std), INTENT (in), DIMENSION(nbpt,nbvmax_in) :: coast_pts !! The coastal flow points (unitless) !! Output VARIABLES INTEGER(i_std), INTENT (in) :: nwbas !! INTEGER(i_std), INTENT (out), DIMENSION(nbpt) :: basin_count !! INTEGER(i_std), INTENT (out), DIMENSION(nbpt) :: basin_notrun !! INTEGER(i_std), INTENT (out), DIMENSION(nbpt,nwbas) :: basin_id !! REAL(r_std), INTENT (out), DIMENSION(nbpt,nwbas,2) :: basin_coor !! REAL(r_std), INTENT (out), DIMENSION(nbpt,nwbas) :: basin_type !! INTEGER(i_std), INTENT (out), DIMENSION(nbpt,nwbas) :: basin_flowdir !! Water flow directions in the basin (unitless) REAL(r_std), INTENT (out), DIMENSION(nbpt,nwbas) :: basin_area !! REAL(r_std), INTENT (out), DIMENSION(nbpt,nwbas,2) :: basin_cg !! Centre of gravity of the HTU REAL(r_std), INTENT (out), DIMENSION(nbpt,nwbas) :: basin_hierarchy !! REAL(r_std), INTENT (out), DIMENSION(nbpt,nwbas) :: basin_fac !! REAL(r_std), INTENT (out), DIMENSION(nbpt,nwbas) :: basin_topoind !! Topographic index of the residence time for a basin (m) REAL(r_std), INTENT (out), DIMENSION(nbpt,nwbas) :: basin_lshead !! Large scale heading out of the grid box. INTEGER(i_std), INTENT (out), DIMENSION(nbpt,nwbas) :: outflow_grid !! Type of outflow on the grid box (unitless) INTEGER(i_std), INTENT (out), DIMENSION(nbpt,nwbas) :: outflow_basin !! INTEGER(i_std), INTENT (out), DIMENSION(nbpt) :: nbcoastal !! INTEGER(i_std), INTENT (out), DIMENSION(nbpt,nwbas) :: coastal_basin !! !! INTEGER(i_std) :: ib !! IF ( nbvmax_in .NE. nbvmax .OR. nbxmax_in .NE. nbxmax ) THEN WRITE(*,*) "GLOBALIZE : nbvmax or nbxmax have changed !!" STOP ENDIF !! DO ib=1,nbpt CALL routing_reg_globalize(nbpt, ib, neighbours, area_bx(ib,:,:), lon_bx(ib,:,:), lat_bx(ib,:,:), trip_bx(ib,:,:), & & hierarchy_bx(ib,:,:), fac_bx(ib,:,:), & & topoind_bx(ib,:,:), min_topoind, nb_basin(ib), basin_inbxid(ib,:), basin_outlet(ib,:,:), basin_outtp(ib,:), & & basin_sz(ib,:), basin_pts(ib,:,:,:), basin_bxout(ib,:), basin_bbout(ib,:), lshead(ib,:), coast_pts(ib,:), nwbas, & & basin_count, basin_notrun, basin_area, basin_cg, basin_hierarchy, basin_fac, basin_topoind, basin_id, basin_coor, & & basin_type, basin_flowdir, basin_lshead, outflow_grid, outflow_basin, nbcoastal, coastal_basin) ENDDO ! END SUBROUTINE globalize SUBROUTINE linkup(nbpt, nbxmax_in, nwbas, basin_count, basin_area, basin_id, basin_flowdir, & & basin_lshead, basin_hierarchy, basin_fac, outflow_grid, outflow_basin, inflow_number, inflow_grid, & & inflow_basin, nbcoastal, coastal_basin, invented_basins) ! USE ioipsl USE grid USE routing_tools USE routing_reg ! ! Arguments ! !! INPUT VARIABLES INTEGER(i_std), INTENT (in) :: nbpt !! Domain size (unitless) INTEGER(i_std), INTENT (in) :: nbxmax_in 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 !! 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_in*8) :: inflow_number !! INTEGER(i_std), INTENT(out), DIMENSION(nbpt,nbxmax_in*8,nbxmax_in*8) :: inflow_basin !! INTEGER(i_std), INTENT(out), DIMENSION(nbpt,nbxmax_in*8,nbxmax_in*8) :: inflow_grid !! ! INTEGER(i_std), INTENT(inout), DIMENSION(nbpt) :: nbcoastal !! INTEGER(i_std), INTENT(inout), DIMENSION(nbpt,nwbas) :: coastal_basin !! ! ! IF ( nbxmax_in .NE. nbxmax ) THEN WRITE(*,*) "LINKUP : nbxmax has changed !!" STOP ENDIF CALL 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) ! END SUBROUTINE linkup SUBROUTINE areanorm(nbpt, nwbas, basin_count, basin_area, outflow_grid, basin_area_norm) ! USE grid ! !! INPUT VARIABLES INTEGER(i_std), INTENT (in) :: nbpt !! Domain size (unitless) INTEGER(i_std), INTENT (in) :: nwbas !! 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(inout) :: outflow_grid !! Type of outflow on the grid box (unitless) REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(out) :: basin_area_norm !! ! !! LOCAL VARIABLES INTEGER(i_std) :: ib, ig REAL(r_std) :: contarea !! REAL(r_std) :: totbasins !! ! ! Normalize the area of all basins ! DO ig=1,nbpt ! totbasins = SUM(basin_area(ig,1:basin_count(ig))) contarea = area(ig)*contfrac(ig) ! DO ib=1,basin_count(ig) basin_area_norm(ig,ib) = basin_area(ig,ib)/totbasins*contarea ! ! Simplify the outflow condition. Large rivers will be reset later in rivclassification. ! We set all outflow points to coastal flow. This will be adjusted later once all procs have ! exchanged their information. So we will have outflow_grid = -2. ! IF ( outflow_grid(ig,ib) .EQ. -1) THEN outflow_grid(ig,ib) = -2 ENDIF ENDDO ! ENDDO ! END SUBROUTINE areanorm SUBROUTINE fetch(nbpt, nwbas, nbcore, corepts, basin_count, basin_area, basin_id, basin_hierarchy, basin_fac, & & outflow_grid, outflow_basin, partial_sum, fetch_basin, outflow_uparea) ! USE ioipsl USE grid USE routing_tools USE routing_reg ! ! Arguments ! !! INPUT VARIABLES INTEGER(i_std), INTENT (in) :: nbpt !! Domain size (unitless) INTEGER(i_std), INTENT (in) :: 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 ! !! IN-OUTPUT VARIABLES REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: fetch_basin !! ! Output REAL(r_std), DIMENSION(nbpt*nwbas), INTENT(out) :: outflow_uparea ! ! !! Local INTEGER(i_std) :: ic INTEGER(i_std), DIMENSION(nbcore) :: fcorepts ! ! DO ic=1,nbcore fcorepts(ic) = corepts(ic)+1 ENDDO ! CALL routing_reg_fetch(nbpt, area, contfrac, nwbas, nbcore, fcorepts, basin_count, basin_area, basin_id, & & basin_hierarchy, basin_fac, outflow_grid, outflow_basin, partial_sum, fetch_basin, outflow_uparea) ! END SUBROUTINE fetch SUBROUTINE rivclassification(nbpt, nwbas, nbcore, corepts, basin_count, outflow_grid, outflow_basin, fetch_basin, & & largest_rivarea, num_largest) ! USE defprec ! !! INPUT VARIABLES INTEGER(i_std), INTENT (in) :: nbpt !! Domain size (unitless) INTEGER(i_std), INTENT (in) :: nwbas !! INTEGER(i_std), INTENT (in) :: nbcore INTEGER(i_std), DIMENSION(nbcore), INTENT(in) :: corepts INTEGER(i_std), DIMENSION(nbpt), INTENT(in) :: basin_count !! INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: outflow_grid INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: outflow_basin !! REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(in) :: fetch_basin !! REAL(r_std), INTENT(in) :: largest_rivarea INTEGER(i_std), INTENT (out) :: num_largest ! !! LOCAL VARIABLES INTEGER(i_std) :: i, ib, ig INTEGER(i_std), DIMENSION(nbcore) :: fcorepts ! ! DO i=1,nbcore fcorepts(i) = corepts(i)+1 ENDDO ! num_largest = 0 ! ! Just a dummy use to avoid f2py warnings outflow_basin(1,1) = outflow_basin(1,1) ! DO i=1,nbcore ig = fcorepts(i) ! DO ib=1,basin_count(ig) ! IF (outflow_grid(ig,ib) .LT. 0 .AND. fetch_basin(ig,ib) .GE. largest_rivarea) THEN num_largest = num_largest + 1 outflow_grid(ig,ib) = -1 ENDIF ! ENDDO ENDDO ! END SUBROUTINE rivclassification SUBROUTINE finish_truncate(nbpt, nbxmax_in, nbasmax, nwbas, num_largest, gridarea, cfrac, basin_count, basin_notrun, basin_area, & & basin_cg, basin_topoind, fetch_basin, basin_id, basin_coor, basin_type, basin_flowdir, outflow_grid, outflow_basin, & & inflow_number, inflow_grid, inflow_basin, routing_area, routing_cg, topo_resid, route_nbbasin, route_togrid, & & route_tobasin, route_nbintobas, global_basinid, route_outlet, route_type, origin_nbintobas, routing_fetch) ! USE ioipsl USE grid USE routing_tools USE routing_reg ! ! Arguments ! !! INPUT VARIABLES INTEGER(i_std), INTENT (in) :: nbpt !! Domain size (unitless) INTEGER(i_std), INTENT (in) :: nbxmax_in, nbasmax INTEGER(i_std), INTENT (in) :: nwbas !! INTEGER(i_std), INTENT (in) :: num_largest REAL(r_std), DIMENSION(nbpt), INTENT(in) :: gridarea REAL(r_std), DIMENSION(nbpt), INTENT(in) :: cfrac ! 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,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_in), INTENT(inout) :: inflow_number !! INTEGER(i_std), DIMENSION(nbpt,nbxmax_in,nbxmax_in), INTENT(inout) :: inflow_basin !! INTEGER(i_std), DIMENSION(nbpt,nbxmax_in,nbxmax_in), INTENT(inout) :: inflow_grid !! ! ! Output variables ! REAL(r_std), DIMENSION(nbpt,nbasmax), INTENT(out) :: routing_area !! Surface of basin (m^2) REAL(r_std), DIMENSION(nbpt,nbasmax), INTENT(out) :: routing_fetch !! Upstream are flowing into HTU (m^2) REAL(r_std), DIMENSION(nbpt,nbasmax,2), INTENT(out) :: routing_cg !! Centre of gravity of HTU (Lat, Lon) REAL(r_std), DIMENSION(nbpt,nbasmax), INTENT(out) :: topo_resid !! Topographic index of the retention time (m) INTEGER(i_std), DIMENSION(nbpt), INTENT(out) :: route_nbbasin !! INTEGER(i_std), DIMENSION(nbpt,nbasmax), INTENT(out) :: route_togrid !! Grid into which the basin flows (unitless) INTEGER(i_std), DIMENSION(nbpt,nbasmax), INTENT(out) :: route_tobasin !! Basin in to which the water goes (unitless) INTEGER(i_std), DIMENSION(nbpt), INTENT(out) :: route_nbintobas !! Number of basin into current one (unitless) INTEGER(i_std), DIMENSION(nbpt), INTENT(out) :: origin_nbintobas !! Number of sub-grid basin into current one before truncation (unitless) INTEGER(i_std), DIMENSION(nbpt,nbasmax), INTENT(out) :: global_basinid !! ID of basin (unitless) REAL(r_std), DIMENSION(nbpt,nbasmax,2), INTENT(out) :: route_outlet !! Coordinate of outlet (-) REAL(r_std), DIMENSION(nbpt,nbasmax), INTENT(out) :: route_type !! Coordinate of outlet (-) ! ! ! CALL routing_reg_end_truncate(nbpt, nbasmax, gridarea, cfrac, nwbas, num_largest, basin_count, & & basin_notrun, basin_area, basin_cg, basin_topoind, fetch_basin, basin_id, basin_coor, basin_type, basin_flowdir, & & outflow_grid, outflow_basin, inflow_number, inflow_grid, inflow_basin) routing_area(:,:) = routing_area_glo(:,:) routing_cg(:,:,:) = routing_cg_glo(:,:,:) topo_resid(:,:) = topo_resid_glo(:,:) route_nbbasin(:) = route_count_glo(:) route_togrid(:,:) = route_togrid_glo(:,:) route_tobasin(:,:) = route_tobasin_glo(:,:) routing_fetch(:,:) = route_fetch_glo(:,:) route_nbintobas(:) = route_nbintobas_glo(:) global_basinid(:,:) = global_basinid_glo(:,:) route_outlet(:,:,:) = route_outlet_glo(:,:,:) route_type(:,:) = route_type_glo(:,:) origin_nbintobas(:) = origin_nbintobas_glo(:) END SUBROUTINE finish_truncate SUBROUTINE killbas(nbpt, nbxmax_in, nbasmax, nwbas,ops, tokill, totakeover, numops, basin_count, basin_area, & & basin_cg, basin_topoind, fetch_basin, basin_id, basin_coor, basin_type, basin_flowdir, outflow_grid, outflow_basin, & & inflow_number, inflow_grid, inflow_basin) ! ! USE ioipsl USE grid USE routing_tools USE routing_reg ! !! INPUT VARIABLES INTEGER(i_std), INTENT (in) :: nbpt !! Domain size (unitless) INTEGER(i_std), INTENT (in) :: nbxmax_in, nbasmax INTEGER(i_std), INTENT (in) :: nwbas !! INTEGER(i_std), INTENT (in) :: ops !! INTEGER(i_std), DIMENSION(nbpt, ops), INTENT (in) :: tokill, totakeover INTEGER(i_std), DIMENSION(nbpt), INTENT (in) :: numops ! INTEGER(i_std), DIMENSION(nbpt), INTENT(inout) :: basin_count !! 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,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_in), INTENT(inout) :: inflow_number !! INTEGER(i_std), DIMENSION(nbpt,nbxmax_in,nbxmax_in), INTENT(inout) :: inflow_basin !! INTEGER(i_std), DIMENSION(nbpt,nbxmax_in,nbxmax_in), INTENT(inout) :: inflow_grid !! !! LOCAL INTEGER(i_std) :: ib, op, tok, totak DO ib=1,nbpt DO op=1,numops(ib) IF (basin_count(ib) > nbasmax) THEN tok = tokill(ib,op) totak = totakeover(ib,op) IF (tok .GT. 0) THEN WRITE(numout,*) "nbpt", ib, "tokill", tok, "totakover", totak CALL routing_reg_killbas(nbpt, ib, tok, totak, nwbas, basin_count, basin_area, basin_cg, basin_topoind,& & fetch_basin, basin_id, basin_coor, basin_type, basin_flowdir, outflow_grid, outflow_basin, inflow_number,& & inflow_grid, inflow_basin) END IF END IF END DO END DO ! REPRENDRE LA FIN DE TRUNCATE END SUBROUTINE killbas SUBROUTINE checkrouting(nbpt, nwbas, outflow_grid, outflow_basin, basin_count) ! USE ioipsl USE grid USE routing_tools USE routing_reg ! !! INPUT VARIABLES INTEGER(i_std), INTENT (in) :: nbpt !! Domain size (unitless) INTEGER(i_std), INTENT (in) :: nwbas !! 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 !! INTEGER(i_std), DIMENSION(nbpt), INTENT(in) :: basin_count !! ! Local INTEGER(i_std), DIMENSION(nbpt,nwbas) :: flag !! INTEGER(i_std) :: ig, ib, it, igrif, ibasf, basnum, test flag(:,:) = 0 WRITE(numout,*) "Checking routing" test = 0 DO ig=1,nbpt basnum = basin_count(ig) DO ib=1,basnum IF (flag(ig,ib) .EQ. 0) THEN flag(ig,ib) = -1 igrif = outflow_grid(ig,ib) ibasf = outflow_basin(ig,ib) ! DO WHILE ((igrif .GT. 0) .AND. (flag(igrif,ibasf) .NE. -1) .AND. (flag(igrif,ibasf) .NE. -99)) flag(igrif, ibasf) = -1 it = outflow_grid(igrif,ibasf) ibasf = outflow_basin(igrif,ibasf) igrif = it END DO IF ((flag(igrif,ibasf) .EQ. -99) .OR. (igrif .LE. 0)) THEN flag(ig,ib) = -99 igrif = outflow_grid(ig,ib) ibasf = outflow_basin(ig,ib) DO WHILE ((igrif .GT. 0) .AND. (flag(igrif,ibasf) .NE. -99)) flag(igrif, ibasf) = -99 it = outflow_grid(igrif,ibasf) ibasf = outflow_basin(igrif,ibasf) igrif = it END DO ELSE IF (flag(igrif,ibasf) .EQ. -1) THEN test = test + 1 flag(ig,ib) = -99 igrif = outflow_grid(ig,ib) ibasf = outflow_basin(ig,ib) DO WHILE (flag(igrif,ibasf) .EQ. -1) flag(igrif, ibasf) = -99 it = outflow_grid(igrif,ibasf) ibasf = outflow_basin(igrif,ibasf) igrif = it END DO END IF END IF END DO END DO WRITE(numout,*) "**** Fetch has",test, "loop errors" END SUBROUTINE checkrouting SUBROUTINE checkfetch(nbpt, nwbas, fetch_basin, outflow_grid, outflow_basin, basin_count) ! USE ioipsl USE grid USE routing_tools USE routing_reg ! !! INPUT VARIABLES INTEGER(i_std), INTENT (in) :: nbpt !! Domain size (unitless) INTEGER(i_std), INTENT (in) :: nwbas !! REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(in) :: fetch_basin !! 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 !! INTEGER(i_std), DIMENSION(nbpt), INTENT(in) :: basin_count !! ! Local INTEGER(i_std) :: ig, ib, it, bt, igrif, ibasf, basnum, test WRITE(numout,*) "Checking Fetch coherence" test = 0 DO ig=1,nbpt basnum = basin_count(ig) DO ib=1,basnum igrif = outflow_grid(ig,ib) ibasf = outflow_basin(ig,ib) IF (igrif .GT. 0) THEN IF (fetch_basin(ig,ib) .GT. fetch_basin(igrif,ibasf)) THEN WRITE(numout,*) "****ERROR Fetch, nbpt:",ig WRITE(numout,*) igrif, ibasf WRITE(numout,*) it, bt test = 1 END IF END IF END DO END DO IF (test .EQ. 0) WRITE(numout,*) "**** Fetch has no error" END SUBROUTINE checkfetch SUBROUTINE finalfetch(nbpt, nbasmax, nbcore, corepts, routing_area, basin_count, route_togrid, route_tobasin, partial_sum, & & fetch_basin, outflow_uparea) ! USE defprec USE ioipsl USE constantes_var ! !! INPUT VARIABLES INTEGER(i_std), INTENT(in) :: nbpt !! Domain size (unitless) INTEGER(i_std), INTENT(in) :: nbasmax !! 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,nbasmax), INTENT(in) :: routing_area !! Surface of basin (m^2) INTEGER(i_std), DIMENSION(nbpt,nbasmax), INTENT(in) :: route_togrid !! Grid into which the basin flows (unitless) INTEGER(i_std), DIMENSION(nbpt,nbasmax), INTENT(in) :: route_tobasin !! Basin in to which the water goes (unitless) REAL(r_std), DIMENSION(nbpt,nbasmax), INTENT(in) :: partial_sum ! !! INPUT/OUTPUT VARIABLES REAL(r_std), DIMENSION(nbpt,nbasmax), INTENT(inout) :: fetch_basin !! ! Out REAL(r_std), DIMENSION(nbpt*nbasmax), INTENT(out) :: outflow_uparea ! ! LOCAL INTEGER(i_std) :: ic, ib, ig, igrif, ibasf, itt, ff(1) INTEGER(i_std), PARAMETER :: maxiterations=10000 INTEGER(i_std), DIMENSION(nbcore) :: fcorepts ! ! DO ic=1,nbcore fcorepts(ic) = corepts(ic)+1 ENDDO ! ! 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 = route_togrid(ig,ib) ibasf = route_tobasin(ig,ib) ! itt = 0 ! DO WHILE ( ibasf .LE. nbasmax ) ! fetch_basin(igrif,ibasf) = MAX(fetch_basin(igrif,ibasf), partial_sum(ig,ib)) ! it = route_togrid(igrif,ibasf) ibasf = route_tobasin(igrif,ibasf) igrif = it itt = itt + 1 IF ( itt .GT. maxiterations) THEN IF ( itt .GT. maxiterations+20) THEN CALL ipslerr_p(3,'routing_interface','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 = fcorepts(ic) ! DO ib=1,basin_count(ig) ! fetch_basin(ig,ib) = fetch_basin(ig,ib) + routing_area(ig,ib) ! igrif = route_togrid(ig,ib) ibasf = route_tobasin(ig,ib) ! itt = 0 ! DO WHILE ( ibasf .LE. nbasmax ) ! fetch_basin(igrif,ibasf) = fetch_basin(igrif,ibasf) + routing_area(ig,ib) ! it = route_togrid(igrif,ibasf) ibasf = route_tobasin(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 ! IF ( itt .GT. maxiterations+20) THEN CALL ipslerr_p(3,'routing_reg_fetch','Problem with computing final fetch','','') ENDIF ENDIF ENDDO ENDDO ENDDO ! ! Get the largest upstream areas in this core ! nboutflow = 0 outflow_uparea(:) = zero ! DO ic=1,nbcore ig = fcorepts(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 ( route_tobasin(ig,ib) .GT. nbasmax) THEN nboutflow = nboutflow + 1 IF ( nboutflow <= nbpt*nbasmax) 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 finalfetch SUBROUTINE finalrivclass(nbpt, nbasmax, nbcore, corepts, route_togrid, route_tobasin, routing_fetch, largest_rivarea, num_largest) ! USE routing_reg ! INTEGER(i_std), INTENT (in) :: nbpt !! Domain size (unitless) INTEGER(i_std), INTENT (in) :: nbasmax !! INTEGER(i_std), INTENT (in) :: nbcore INTEGER(i_std), DIMENSION(nbcore), INTENT(in) :: corepts INTEGER(i_std), DIMENSION(nbpt,nbasmax), INTENT(inout) :: route_togrid !! Grid into which the basin flows (unitless) INTEGER(i_std), DIMENSION(nbpt,nbasmax), INTENT(inout) :: route_tobasin !! Basin in to which the water goes (unitless) REAL(r_std), DIMENSION(nbpt,nbasmax), INTENT(inout) :: routing_fetch !! Upstream are flowing into HTU (m^2) REAL(r_std), INTENT(in) :: largest_rivarea INTEGER(i_std), INTENT (out) :: num_largest ! ! LOCAL INTEGER(i_std) :: i, ib, ig INTEGER(i_std), DIMENSION(nbcore) :: fcorepts ! ! DO i=1,nbcore fcorepts(i) = corepts(i)+1 ENDDO ! ! Just a dummy use to avoid f2py warnings route_togrid(1,1) = route_togrid(1,1) ! ! Redo the the distinction between river outflow and coastal flow. We can not ! take into account the return flow points. ! num_largest = 0 DO i = 1, nbcore ig = fcorepts(i) DO ib = 1, nbasmax IF (route_tobasin_glo(ig,ib) .GT. nbasmax) THEN IF (routing_fetch(ig,ib) .GE. largest_rivarea ) THEN num_largest = num_largest + 1 route_tobasin_glo(ig,ib) = nbasmax + 3 route_tobasin(ig,ib) = nbasmax + 3 ENDIF ENDIF ENDDO ENDDO ! WRITE(numout,*) 'NUMBER OF RIVERS :', COUNT(route_tobasin_glo .GE. nbasmax + 3) ! END SUBROUTINE finalrivclass