Skip to content

GitLab

Commit da8221a4 authored by POLCHER Jan's avatar POLCHER Jan 🚴🏾
Browse files

Add the diagnostic of the centre of gravity of the HTU to thegraph description.

parent 98f59ec8
Show whitespace changes
Inline Side-by-side
Showing with 124 additions and 62 deletions
DocumentationInterface
View file @ da8221a4
......@@ -97,13 +97,15 @@ basin_pts : rank-4 array('i') with bounds (nbpt,nbvmax_in,nbvmax_in,2)
basin_lshead : rank-2 array('f') with bounds (nbpt,nbvmax_in)
coast_pts : rank-2 array('i') with bounds (nbpt,nbvmax_in)
====================================================================
basin_count,basin_notrun,basin_area,basin_hierarchy,basin_fac,basin_topoind,basin_id,basin_coor,basin_type,basin_flowdir,basin_lshead,outflow_grid,outflow_basin,nbcoastal,coastal_basin = globalize(area_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,[nbpt,nbvmax_in,nbxmax_in])
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 = globalize(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,[nbpt,nbvmax_in,nbxmax_in])
Wrapper for ``globalize``.
Parameters
----------
area_bx : input rank-3 array('f') with bounds (nbpt,nbxmax_in,nbxmax_in)
lon_bx : input rank-3 array('f') with bounds (nbpt,nbxmax_in,nbxmax_in)
lat_bx : input rank-3 array('f') with bounds (nbpt,nbxmax_in,nbxmax_in)
trip_bx : input rank-3 array('i') with bounds (nbpt,nbxmax_in,nbxmax_in)
hierarchy_bx : input rank-3 array('f') with bounds (nbpt,nbxmax_in,nbxmax_in)
fac_bx : input rank-3 array('f') with bounds (nbpt,nbxmax_in,nbxmax_in)
......@@ -135,6 +137,7 @@ Returns
basin_count : rank-1 array('i') with bounds (nbpt)
basin_notrun : rank-1 array('i') with bounds (nbpt)
basin_area : rank-2 array('f') with bounds (nbpt,nwbas)
basin_cg : rank-3 array('f') with bounds (nbpt,nwbas,2)
basin_hierarchy : rank-2 array('f') with bounds (nbpt,nwbas)
basin_fac : rank-2 array('f') with bounds (nbpt,nwbas)
basin_topoind : rank-2 array('f') with bounds (nbpt,nwbas)
......@@ -206,7 +209,7 @@ Returns
-------
fetch_basin : rank-2 array('f') with bounds (nbpt,nwbas)
====================================================================
routing_area,topo_resid,route_togrid,route_tobasin,route_nbintobas,global_basinid,route_outlet,route_type,origin_nbintobas = truncate(nbasmax,basin_count,basin_notrun,basin_area,basin_topoind,fetch_basin,basin_id,basin_coor,basin_type,basin_flowdir,outflow_grid,outflow_basin,inflow_number,inflow_grid,inflow_basin,[nbpt,nbxmax_in,nwbas])
routing_area,routing_cg,topo_resid,route_togrid,route_tobasin,route_nbintobas,global_basinid,route_outlet,route_type,origin_nbintobas = truncate(nbasmax,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,[nbpt,nbxmax_in,nwbas])
Wrapper for ``truncate``.
......@@ -216,6 +219,7 @@ nbasmax : input int
basin_count : in/output rank-1 array('i') with bounds (nbpt)
basin_notrun : in/output rank-1 array('i') with bounds (nbpt)
basin_area : in/output rank-2 array('f') with bounds (nbpt,nwbas)
basin_cg : in/output rank-3 array('f') with bounds (nbpt,nwbas,2)
basin_topoind : in/output rank-2 array('f') with bounds (nbpt,nwbas)
fetch_basin : in/output rank-2 array('f') with bounds (nbpt,nwbas)
basin_id : in/output rank-2 array('i') with bounds (nbpt,nwbas)
......@@ -240,6 +244,7 @@ nwbas : input int, optional
Returns
-------
routing_area : rank-2 array('f') with bounds (nbpt,nbasmax)
routing_cg : rank-3 array('f') with bounds (nbpt,nbasmax,2)
topo_resid : rank-2 array('f') with bounds (nbpt,nbasmax)
route_togrid : rank-2 array('i') with bounds (nbpt,nbasmax)
route_tobasin : rank-2 array('i') with bounds (nbpt,nbasmax)
......
F90subroutines/routing_interface.f90
View file @ da8221a4
......@@ -154,9 +154,9 @@ SUBROUTINE findbasins(nbpt, nbvmax_in, nbxmax_in, nbi, nbj, trip_bx, basin_bx, f
END SUBROUTINE findbasins
SUBROUTINE globalize(nbpt, nbvmax_in, nbxmax_in, area_bx, trip_bx, hierarchy_bx, fac_bx, topoind_bx, &
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_hierarchy, &
& 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)
!
......@@ -171,6 +171,8 @@ SUBROUTINE globalize(nbpt, nbvmax_in, nbxmax_in, area_bx, trip_bx, hierarchy_bx,
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
......@@ -194,6 +196,7 @@ SUBROUTINE globalize(nbpt, nbvmax_in, nbxmax_in, area_bx, trip_bx, hierarchy_bx,
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)
......@@ -211,10 +214,11 @@ SUBROUTINE globalize(nbpt, nbvmax_in, nbxmax_in, area_bx, trip_bx, hierarchy_bx,
ENDIF
!!
DO ib=1,nbpt
CALL routing_reg_globalize(nbpt, ib, neighbours, area_bx(ib,:,:), trip_bx(ib,:,:), hierarchy_bx(ib,:,:), fac_bx(ib,:,:), &
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_hierarchy, basin_fac, basin_topoind, basin_id, basin_coor, &
& 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
!
......@@ -296,9 +300,9 @@ SUBROUTINE fetch(nbpt, nwbas, basin_count, basin_area, basin_id, basin_hierarchy
END SUBROUTINE fetch
SUBROUTINE truncate(nbpt, nbxmax_in, nbasmax, nwbas, basin_count, basin_notrun, basin_area, basin_topoind,&
SUBROUTINE truncate(nbpt, nbxmax_in, nbasmax, nwbas, 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, topo_resid, route_togrid, route_tobasin, route_nbintobas, &
& inflow_grid, inflow_basin, routing_area, routing_cg, topo_resid, route_togrid, route_tobasin, route_nbintobas, &
& global_basinid, route_outlet, route_type, origin_nbintobas)
!
USE ioipsl
......@@ -320,6 +324,7 @@ SUBROUTINE truncate(nbpt, nbxmax_in, nbasmax, nwbas, basin_count, basin_notrun,
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)
......@@ -334,6 +339,7 @@ SUBROUTINE truncate(nbpt, nbxmax_in, nbasmax, nwbas, basin_count, basin_notrun,
! Output variables
!
REAL(r_std), DIMENSION(nbpt,nbasmax), INTENT(out) :: routing_area !! Surface of basin (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,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)
......@@ -348,11 +354,12 @@ SUBROUTINE truncate(nbpt, nbxmax_in, nbasmax, nwbas, basin_count, basin_notrun,
STOP
ENDIF
CALL routing_reg_truncate(nbpt, nbasmax, area, contfrac, nwbas, basin_count, basin_notrun, basin_area, basin_topoind,&
& fetch_basin, basin_id, basin_coor, basin_type, basin_flowdir, outflow_grid, outflow_basin, inflow_number,&
CALL routing_reg_truncate(nbpt, nbasmax, area, contfrac, nwbas, 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_togrid(:,:) = route_togrid_glo(:,:)
route_tobasin(:,:) = route_tobasin_glo(:,:)
......
F90subroutines/routing_reg.f90
View file @ da8221a4
......@@ -16,6 +16,7 @@ MODULE routing_reg
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_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_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)
......@@ -1633,9 +1634,9 @@ CONTAINS
!! \n
!_ ================================================================================================================================
SUBROUTINE routing_reg_globalize(nbpt, ib, neighbours, area_bx, trip_bx, hierarchy_bx, fac_bx, topoind_bx, &
SUBROUTINE routing_reg_globalize(nbpt, ib, neighbours, 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_hierarchy, &
& 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)
!
......@@ -1648,6 +1649,8 @@ SUBROUTINE routing_reg_globalize(nbpt, ib, neighbours, area_bx, trip_bx, hierarc
!! (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) :: fac_bx !! Level in the basin of the point
......@@ -1655,7 +1658,7 @@ SUBROUTINE routing_reg_globalize(nbpt, ib, neighbours, area_bx, trip_bx, hierarc
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
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
......@@ -1667,10 +1670,11 @@ SUBROUTINE routing_reg_globalize(nbpt, ib, neighbours, area_bx, trip_bx, hierarc
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 !!
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_fac !!
REAL(r_std), DIMENSION(nbpt,nwbas) :: basin_topoind !! Topographic index of the residence time for a basin (m)
......@@ -1747,6 +1751,7 @@ SUBROUTINE routing_reg_globalize(nbpt, ib, neighbours, area_bx, trip_bx, hierarc
! 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
!
......@@ -1762,6 +1767,11 @@ SUBROUTINE routing_reg_globalize(nbpt, ib, neighbours, area_bx, trip_bx, hierarc
!
basin_area(ib,ij) = basin_area(ib,ij) + area_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
......@@ -2847,9 +2857,9 @@ SUBROUTINE routing_reg_fetch(nbpt, gridarea, contfrac, nwbas, basin_count, basin
!! \n
!_ ================================================================================================================================
SUBROUTINE routing_reg_truncate(nbpt, nbasmax, gridarea, contfrac, nwbas, basin_count, basin_notrun, basin_area, basin_topoind,&
& fetch_basin, basin_id, basin_coor, basin_type, basin_flowdir, outflow_grid, outflow_basin, inflow_number,&
& inflow_grid, inflow_basin)
SUBROUTINE routing_reg_truncate(nbpt, nbasmax, gridarea, contfrac, nwbas, 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)
!
IMPLICIT NONE
!
......@@ -2872,6 +2882,7 @@ SUBROUTINE routing_reg_truncate(nbpt, nbasmax, gridarea, contfrac, nwbas, basin_
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)
......@@ -3137,7 +3148,7 @@ SUBROUTINE routing_reg_truncate(nbpt, nbasmax, gridarea, contfrac, nwbas, basin_
! Then we call routing_reg_killbas to clean up the basins in this grid
!
IF ( kbas .GT. 0 .AND. sbas .GT. 0 ) THEN
CALL routing_reg_killbas(nbpt, ib, kbas, sbas, nwbas, basin_count, basin_area, basin_topoind,&
CALL routing_reg_killbas(nbpt, ib, kbas, sbas, 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)
ENDIF
......@@ -3189,7 +3200,7 @@ SUBROUTINE routing_reg_truncate(nbpt, nbasmax, gridarea, contfrac, nwbas, basin_
! Then we call routing_reg_killbas to clean up the basins in this grid
!
IF ( kbas .GT. 0 .AND. sbas .GT. 0 ) THEN
CALL routing_reg_killbas(nbpt, ib, kbas, sbas, nwbas, basin_count, basin_area, basin_topoind,&
CALL routing_reg_killbas(nbpt, ib, kbas, sbas, 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)
ENDIF
......@@ -3248,7 +3259,7 @@ SUBROUTINE routing_reg_truncate(nbpt, nbasmax, gridarea, contfrac, nwbas, basin_
ENDIF
!
IF ( kbas .GT. 0 .AND. sbas .GT. 0 ) THEN
CALL routing_reg_killbas(nbpt, ib, kbas, sbas, nwbas, basin_count, basin_area, basin_topoind,&
CALL routing_reg_killbas(nbpt, ib, kbas, sbas, 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)
ENDIF
......@@ -3282,6 +3293,7 @@ SUBROUTINE routing_reg_truncate(nbpt, nbasmax, gridarea, contfrac, nwbas, basin_
DO ij=1,basin_count(ib)
!
routing_area_glo(ib,ij) = basin_area(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)
......@@ -3441,7 +3453,7 @@ SUBROUTINE routing_reg_truncate(nbpt, nbasmax, gridarea, contfrac, nwbas, basin_
!! \n
!_ ================================================================================================================================
SUBROUTINE routing_reg_killbas(nbpt, ib, tokill, totakeover, nwbas, basin_count, basin_area, basin_topoind,&
SUBROUTINE routing_reg_killbas(nbpt, ib, tokill, totakeover, 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)
!
......@@ -3460,6 +3472,7 @@ SUBROUTINE routing_reg_killbas(nbpt, ib, tokill, totakeover, nwbas, basin_count,
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 !!
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)
......@@ -3482,6 +3495,10 @@ SUBROUTINE routing_reg_killbas(nbpt, ib, tokill, totakeover, nwbas, basin_count,
!
!
basin_area(ib, totakeover) = basin_area(ib, totakeover) + basin_area(ib, tokill)
!
basin_cg(ib, totakeover, 1) = (basin_cg(ib, totakeover, 1) + basin_cg(ib, tokill, 1))/2.0
basin_cg(ib, totakeover, 2) = (basin_cg(ib, totakeover, 2) + basin_cg(ib, tokill, 2))/2.0
!
basin_topoind(ib, totakeover) = (basin_topoind(ib, totakeover) + basin_topoind(ib, tokill))/2.0
!
! Add the fetch of the basin will kill to the one which gets the water
......@@ -3563,6 +3580,8 @@ SUBROUTINE routing_reg_killbas(nbpt, ib, tokill, totakeover, nwbas, basin_count,
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))
......@@ -3615,6 +3634,7 @@ SUBROUTINE routing_reg_killbas(nbpt, ib, tokill, totakeover, nwbas, basin_count,
nbasmax_save = nbasmax
ALLOCATE (routing_area_glo(nbpt,nbasmax), stat=ier)
ALLOCATE (routing_cg_glo(nbpt,nbasmax,2), stat=ier)
ALLOCATE (topo_resid_glo(nbpt,nbasmax), stat=ier)
ALLOCATE (route_togrid_glo(nbpt,nbasmax), stat=ier)
ALLOCATE (route_tobasin_glo(nbpt,nbasmax), stat=ier)
......
Interface.py
View file @ da8221a4
......@@ -60,6 +60,7 @@ def initatmgrid(nbpt, modelgrid) :
return
#
#
#
class HydroOverlap :
#
def __init__(self, nbpt, nbvmax, sub_pts, sub_index_in, sub_area_in, sub_lon_in, sub_lat_in, modelgrid, hydrodata) :
......@@ -110,11 +111,11 @@ class HydroOverlap :
# Clean-up these arrays so that they are easy to use in Python.
self.lon_bx[self.lon_bx > 360.]=np.nan
self.lat_bx[self.lat_bx > 90.]=np.nan
# with open("Grids.txt", 'wb') as fp:
# pickle.dump(lon_bx, fp)
# pickle.dump(lat_bx, fp)
#
return
#
#
#
class HydroSuper :
def __init__(self, nbvmax, hydrodata, hydrooverlap) :
#
......@@ -128,13 +129,18 @@ class HydroSuper :
#
# Call Globalize
#
self.basin_count, self.basin_notrun, self.basin_area, self.basin_hierarchy, self.basin_fac, self.basin_topoind, \
self.basin_id, self.basin_coor, self.basin_type, self.basin_flowdir, \
lon_bx_tmp = hydrooverlap.lon_bx
lon_bx_tmp[np.isnan(lon_bx_tmp)] = undef_int
lat_bx_tmp = hydrooverlap.lat_bx
lat_bx_tmp[np.isnan(lat_bx_tmp)] = undef_int
self.basin_count, self.basin_notrun, self.basin_area, self.basin_cg, self.basin_hierarchy, self.basin_fac, self.basin_topoind, \
self.basin_id, self.basin_outcoor, self.basin_type, self.basin_flowdir, \
self.basin_lshead, self.outflow_grid, self.outflow_basin, self.nbcoastal, self.coastal_basin = \
routing_interface.globalize(hydrooverlap.area_bx, hydrooverlap.trip_bx, hydrooverlap.hierarchy_bx,\
hydrooverlap.fac_bx, hydrooverlap.topoind_bx, hydrodata.topoindmin, nb_basin, \
basin_inbxid, basin_outlet, basin_outtp, self.basin_sz, self.basin_pts, basin_bxout, \
routing_interface.globalize(hydrooverlap.area_bx, lon_bx_tmp, lat_bx_tmp, hydrooverlap.trip_bx, \
hydrooverlap.hierarchy_bx, hydrooverlap.fac_bx, hydrooverlap.topoind_bx, hydrodata.topoindmin, \
nb_basin, basin_inbxid, basin_outlet, basin_outtp, self.basin_sz, self.basin_pts, basin_bxout, \
basin_bbout, basin_lshead, coast_pts, hydrooverlap.nwbas)
return
#
def linkup(self, hydrodata) :
#
......@@ -145,22 +151,27 @@ class HydroSuper :
self.basin_flowdir, self.basin_lshead, self.basin_hierarchy, \
self.basin_fac, self.outflow_grid, self.outflow_basin, \
self.nbcoastal, self.coastal_basin, float(hydrodata.basinsmax))
return
#
def fetch(self) :
self.fetch_basin = routing_interface.fetch(self.basin_count, self.basin_area, self.basin_id, self.basin_hierarchy, \
self.basin_fac, self.outflow_grid, self.outflow_basin)
return
#
#
#
class HydroGraph :
def __init__(self, nbasmax, hydrosuper) :
self.nbasmax = nbasmax
self.routing_area, self.topo_resid, self.route_togrid, self.route_tobasin, self.route_nbintobas, self.global_basinid, \
self.routing_area, self.routing_cg, self.topo_resid, self.route_togrid, self.route_tobasin, self.route_nbintobas, self.global_basinid, \
self.route_outlet, self.route_type, self.origin_nbintobas = \
routing_interface.truncate(nbasmax, hydrosuper.basin_count, hydrosuper.basin_notrun, hydrosuper.basin_area, \
hydrosuper.basin_topoind, hydrosuper.fetch_basin, hydrosuper.basin_id, \
hydrosuper.basin_coor, hydrosuper.basin_type, hydrosuper.basin_flowdir, \
hydrosuper.basin_cg, hydrosuper.basin_topoind, hydrosuper.fetch_basin, hydrosuper.basin_id, \
hydrosuper.basin_outcoor, hydrosuper.basin_type, hydrosuper.basin_flowdir, \
hydrosuper.outflow_grid, hydrosuper.outflow_basin, \
hydrosuper.inflow_number,hydrosuper.inflow_grid,hydrosuper.inflow_basin)
return
#
def dumpnetcdf(self, filename, globalgrid, procgrid, part) :
#
NCFillValue=1.0e20
......@@ -232,7 +243,7 @@ class HydroGraph :
# Once gathered on root_proc we transform them into float64. Careful, Integer variables do not have NaN !
#
rarea = procgrid.landscatter(self.routing_area[:,:], order='F')
rarea = rarea.astype(np.float64, copy=False)
rarea = rarea.astype(vtyp, copy=False)
rarea[np.isnan(rarea)] = NCFillValue
if part.rank == 0 :
routingarea = outnf.createVariable("routingarea", vtyp, ('htu','y','x'), fill_value=NCFillValue)
......@@ -243,7 +254,7 @@ class HydroGraph :
routingarea[:,:,:] = part.gather(rarea)
#
rgrid = procgrid.landscatter(self.route_togrid[:,:], order='F')
rgrid = rgrid.astype(np.float64, copy=False)
rgrid = rgrid.astype(vtyp, copy=False)
rgrid[rgrid >= RPP.IntFillValue] = NCFillValue
if part.rank == 0 :
routetogrid = outnf.createVariable("routetogrid", vtyp, ('htu','y','x'), fill_value=NCFillValue)
......@@ -254,7 +265,7 @@ class HydroGraph :
routetogrid[:,:,:] = part.gather(rgrid)
#
rtobasin = procgrid.landscatter(self.route_tobasin[:,:], order='F')
rtobasin = rtobasin.astype(np.float64, copy=False)
rtobasin = rtobasin.astype(vtyp, copy=False)
rtobasin[rtobasin >= RPP.IntFillValue] = NCFillValue
if part.rank == 0 :
routetobasin = outnf.createVariable("routetobasin", vtyp, ('htu','y','x'), fill_value=NCFillValue)
......@@ -265,7 +276,7 @@ class HydroGraph :
routetobasin[:,:,:] = part.gather(rtobasin)
#
rid = procgrid.landscatter(self.global_basinid[:,:], order='F')
rid = rid.astype(np.float64, copy=False)
rid = rid.astype(vtyp, copy=False)
rid[rid >= RPP.IntFillValue] = NCFillValue
if part.rank == 0 :
basinid = outnf.createVariable("basinid", vtyp, ('htu','y','x'), fill_value=NCFillValue)
......@@ -276,7 +287,7 @@ class HydroGraph :
basinid[:,:,:] = part.gather(rid)
#
rintobas = procgrid.landscatter(self.route_nbintobas[:])
rintobas = rintobas.astype(np.float64, copy=False)
rintobas = rintobas.astype(vtyp, copy=False)
rintobas[rintobas >= RPP.IntFillValue] = NCFillValue
if part.rank == 0 :
routenbintobas = outnf.createVariable("routenbintobas", vtyp, ('y','x'), fill_value=NCFillValue)
......@@ -287,7 +298,7 @@ class HydroGraph :
routenbintobas[:,:] = part.gather(rintobas)
#
onbintobas = procgrid.landscatter(self.origin_nbintobas[:])
onbintobas = onbintobas.astype(np.float64, copy=False)
onbintobas = onbintobas.astype(vtyp, copy=False)
onbintobas[onbintobas >= RPP.IntFillValue] = NCFillValue
if part.rank == 0 :
originnbintobas = outnf.createVariable("originnbintobas", vtyp, ('y','x'), fill_value=NCFillValue)
......@@ -298,7 +309,7 @@ class HydroGraph :
originnbintobas[:,:] = part.gather(onbintobas)
#
olat = procgrid.landscatter(self.route_outlet[:,:,0], order='F')
olat = olat.astype(np.float64, copy=False)
olat = olat.astype(vtyp, copy=False)
olat[np.isnan(olat)] = NCFillValue
if part.rank == 0 :
outletlat = outnf.createVariable("outletlat", vtyp, ('htu','y','x'), fill_value=NCFillValue)
......@@ -309,7 +320,7 @@ class HydroGraph :
outletlat[:,:,:] = part.gather(olat)
#
olon = procgrid.landscatter(self.route_outlet[:,:,1], order='F')
olon = olon.astype(np.float64, copy=False)
olon = olon.astype(vtyp, copy=False)
olon[np.isnan(olon)] = NCFillValue
if part.rank == 0 :
outletlon = outnf.createVariable("outletlon", vtyp, ('htu','y','x'), fill_value=NCFillValue)
......@@ -320,7 +331,7 @@ class HydroGraph :
outletlon[:,:,:] = part.gather(olon)
#
otype = procgrid.landscatter(self.route_type[:,:], order='F')
otype = otype.astype(np.float64, copy=False)
otype = otype.astype(vtyp, copy=False)
otype[np.isnan(otype)] = NCFillValue
if part.rank == 0 :
outlettype = outnf.createVariable("outlettype", vtyp, ('htu','y','x'), fill_value=NCFillValue)
......@@ -331,7 +342,7 @@ class HydroGraph :
outlettype[:,:,:] = part.gather(otype)
#
tind = procgrid.landscatter(self.topo_resid[:,:], order='F')
tind = tind.astype(np.float64, copy=False)
tind = tind.astype(vtyp, copy=False)
tind[np.isnan(tind)] = NCFillValue
if part.rank == 0 :
topoindex = outnf.createVariable("topoindex", vtyp, ('htu','y','x'), fill_value=NCFillValue)
......@@ -341,9 +352,28 @@ class HydroGraph :
topoindex = np.zeros((1,1,1))
topoindex[:,:,:] = part.gather(tind)
#
# Save centre of gravity of HTU
#
cg = procgrid.landscatter(self.routing_cg[:,:,:], order='F')
cg = cg.astype(vtyp, copy=False)
cg[np.isnan(cg)] = NCFillValue
if part.rank == 0 :
CG_lon = outnf.createVariable("CG_lon", vtyp, ('htu','y','x'), fill_value=NCFillValue)
CG_lon.title = "Longitude of centre of gravity of HTU"
CG_lon.units = "degrees east"
CG_lat = outnf.createVariable("CG_lat", vtyp, ('htu','y','x'), fill_value=NCFillValue)
CG_lat.title = "Latitude of centre of gravity of HTU"
CG_lat.units = "degrees north"
else :
CG_lon = np.zeros((1,1,1))
CG_lat = np.zeros((1,1,1))
CG_lon[:,:,:] = part.gather(cg[1,:,:,:])
CG_lat[:,:,:] = part.gather(cg[0,:,:,:])
#
if part.rank == 0 :
outnf.close()
#
return
......
TestConfigs/run.def.WRF_local
View file @ da8221a4
......@@ -19,7 +19,7 @@ nbasmax = 35
#
# Output
#
GraphFile = test_graph.nc
GraphFile = MEDCORDEX_test_graph.nc
#
# Diagnostics
# You need to provide an interval in longitude and Latitude.
......
run.def
View file @ da8221a4
......@@ -19,7 +19,7 @@ nbasmax = 35
#
# Output
#
GraphFile = test_graph.nc
GraphFile = MEDCORDEX_test_graph.nc
#
# Diagnostics
# You need to provide an interval in longitude and Latitude.
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment