From da8221a4aa23c0804520824ce0c1c8035046945d Mon Sep 17 00:00:00 2001
From: Jan Polcher <jan.polcher@lmd.jussieu.fr>
Date: Thu, 13 Jun 2019 11:42:18 +0200
Subject: [PATCH] Add the diagnostic of the centre of gravity of the HTU to
thegraph description.
---
DocumentationInterface | 9 ++-
F90subroutines/routing_interface.f90 | 37 +++++++-----
F90subroutines/routing_reg.f90 | 52 +++++++++++------
Interface.py | 84 +++++++++++++++++++---------
TestConfigs/run.def.WRF_local | 2 +-
run.def | 2 +-
6 files changed, 124 insertions(+), 62 deletions(-)
diff --git a/DocumentationInterface b/DocumentationInterface
index e15fe47..8770db9 100644
--- a/DocumentationInterface
+++ b/DocumentationInterface
@@ -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)
diff --git a/F90subroutines/routing_interface.f90 b/F90subroutines/routing_interface.f90
index 03a3c93..e6e4e20 100644
--- a/F90subroutines/routing_interface.f90
+++ b/F90subroutines/routing_interface.f90
@@ -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)
!
@@ -169,8 +169,10 @@ SUBROUTINE globalize(nbpt, nbvmax_in, nbxmax_in, area_bx, trip_bx, hierarchy_bx,
!
!! INPUT VARIABLES
INTEGER(i_std), INTENT (in) :: nbpt !! Domain size (unitless)
- INTEGER(i_std), INTENT(in) :: nbvmax_in, nbxmax_in
+ 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)
@@ -333,26 +338,28 @@ 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), INTENT(out) :: topo_resid !! Topographic index of the retention time (m)
+ 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)
- 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), 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 (-)
+ 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 (-)
!
IF ( nbxmax_in .NE. nbxmax ) THEN
WRITE(*,*) "TRUNCATE : nbxmax has changed !!"
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(:,:)
diff --git a/F90subroutines/routing_reg.f90 b/F90subroutines/routing_reg.f90
index c18bb7f..039c19c 100644
--- a/F90subroutines/routing_reg.f90
+++ b/F90subroutines/routing_reg.f90
@@ -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)
!
@@ -3456,12 +3468,13 @@ SUBROUTINE routing_reg_killbas(nbpt, ib, tokill, totakeover, nwbas, basin_count,
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 !!
+ 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_topoind !! Topographic index of the residence time for a basin (m)
- REAL(r_std), DIMENSION(nbpt,nwbas) :: fetch_basin !!
+ 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)
INTEGER(i_std), DIMENSION(nbpt,nwbas) :: outflow_basin !!
!
@@ -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)
diff --git a/Interface.py b/Interface.py
index 23750b8..6ba1114 100644
--- a/Interface.py
+++ b/Interface.py
@@ -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,14 +129,19 @@ 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) :
#
# Call the linkup routine in routing_reg.
@@ -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
diff --git a/TestConfigs/run.def.WRF_local b/TestConfigs/run.def.WRF_local
index f65aa6e..8c7083b 100644
--- a/TestConfigs/run.def.WRF_local
+++ b/TestConfigs/run.def.WRF_local
@@ -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.
diff --git a/run.def b/run.def
index f65aa6e..8c7083b 100644
--- a/run.def
+++ b/run.def
@@ -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.
--
GitLab