Skip to content

GitLab

Docker-in-Docker (DinD) capabilities of public runners deactivated. More info

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

Merge branch 'newdivbas' into 'master' 

This is the end of the newdivbas branch

See merge request !2
parents 14d02f54 e7f2480b
Expand all Hide whitespace changes
Inline Side-by-side
Showing with 1181 additions and 241 deletions
F90subroutines/Makefile.JeanZay 0 → 100644
View file @ 57afb669
## Jean Zay compilation
# From environtment
# module purge
# source /gpfslocalsup/spack_soft/environment-modules/current/init/ksh
# module load python/3.7.3
# conda activate MPI
#
#
obj = defprec.o ioipsl.o constantes_var.o haversine.o grid.o routing_reg.o routing_tools.o
#
FC = /gpfslocalsup/pub/anaconda-py3/2019.03/bin/x86_64-conda_cos6-linux-gnu-gfortran
FFLAG = -fPIC
#FDBG = -g -fbounds-check
#
FPY = /linkhome/rech/genlmd01/rron974/.conda/envs/MPI/bin/f2py
#PDBG = --debug
#
EXT_SUFFIX := $(shell python3-config --extension-suffix)
#
all : routing_interface$(EXT_SUFFIX) diagst$(EXT_SUFFIX)
routing_interface : routing_interface$(EXT_SUFFIX)
diagst : diagst$(EXT_SUFFIX)
clean :
rm -f compilation.log *.o *.so *.mod routing_interface.so diagst.so *~
routing_interface$(EXT_SUFFIX) : $(obj) routing_interface.f90
rm -f compilation.log
$(FPY) -c routing_interface.f90 $(PDBG) --fcompiler=gfortran -m routing_interface $(obj) > compilation.log
diagst$(EXT_SUFFIX) : diagst.f90
$(FPY) -c diagst.f90 $(PDBG) --fcompiler=gfortran -m diagst > compilation.log
defprec.o : defprec.f90
$(FC) $(FFLAG) $(FDBG) -c defprec.f90
ioipsl.o : ioipsl.f90
$(FC) $(FFLAG) $(FDBG) -c ioipsl.f90
constantes_var.o : constantes_var.f90
$(FC) $(FFLAG) $(FDBG) -c constantes_var.f90
haversine.o : haversine.f90
$(FC) $(FFLAG) $(FDBG) -c haversine.f90
grid.o : defprec.f90 haversine.o constantes_var.o grid.f90
$(FC) $(FFLAG) $(FDBG) -c grid.f90
routing_tools.o : defprec.o grid.o constantes_var.o routing_tools.f90
$(FC) $(FFLAG) $(FDBG) -c routing_tools.f90
routing_reg.o : defprec.f90 grid.o constantes_var.o routing_tools.o routing_reg.f90
$(FC) $(FFLAG) $(FDBG) -c routing_reg.f90
F90subroutines/Makefile.JeanZay_ifort 0 → 100644
View file @ 57afb669
## Jean Zay compilation
# From environtment
# module purge
#module purge
# module load intel-all/19.0.4
# eval "$(/gpfslocalsup/pub/anaconda-py3/2021.05/bin/conda shell.bash hook)"
# conda create -y -n fastmpi python=3.8
# conda activate fastmpi
#
# export MPICC=$(which mpicc)
# echo $MPICC
# pip install mpi4py --user --no-cache-dir
# conda install numba astropy netcdf4 matplotlib basemap ConfigParser cartopy
## conda install -c anaconda gfortran_linux-64
# git clone https://github.com/spacetelescope/spherical_geometry.git
# cd spherical_geometry
# python setup.py install --user
#
#
obj = defprec.o ioipsl.o constantes_var.o haversine.o grid.o routing_reg.o routing_tools.o
#
FC = /gpfslocalsys/intel/parallel_studio_xe_2019_update4_cluster_edition/compilers_and_libraries_2019.4.243/linux/mpi/intel64/bin/mpiifort
FFLAG = -fPIC
#FDBG = -g -fbounds-check
#
FPY = /linkhome/rech/genlmd01/rron974/.conda/envs/fastmpi/bin/f2py
#PDBG = --debug
#
EXT_SUFFIX := $(shell python3-config --extension-suffix)
#
all : routing_interface$(EXT_SUFFIX) diagst$(EXT_SUFFIX)
routing_interface : routing_interface$(EXT_SUFFIX)
diagst : diagst$(EXT_SUFFIX)
clean :
rm -f compilation.log *.o *.so *.mod routing_interface.so diagst.so *~
routing_interface$(EXT_SUFFIX) : $(obj) routing_interface.f90
rm -f compilation.log
$(FPY) -c routing_interface.f90 $(PDBG) --fcompiler=intelem -m routing_interface $(obj) > compilation.log
diagst$(EXT_SUFFIX) : diagst.f90
$(FPY) -c diagst.f90 $(PDBG) --fcompiler=intelem -m diagst > compilation.log
defprec.o : defprec.f90
$(FC) $(FFLAG) $(FDBG) -c defprec.f90
ioipsl.o : ioipsl.f90
$(FC) $(FFLAG) $(FDBG) -c ioipsl.f90
constantes_var.o : constantes_var.f90
$(FC) $(FFLAG) $(FDBG) -c constantes_var.f90
haversine.o : haversine.f90
$(FC) $(FFLAG) $(FDBG) -c haversine.f90
grid.o : defprec.f90 haversine.o constantes_var.o grid.f90
$(FC) $(FFLAG) $(FDBG) -c grid.f90
routing_tools.o : defprec.o grid.o constantes_var.o routing_tools.f90
$(FC) $(FFLAG) $(FDBG) -c routing_tools.f90
routing_reg.o : defprec.f90 grid.o constantes_var.o routing_tools.o routing_reg.f90
$(FC) $(FFLAG) $(FDBG) -c routing_reg.f90
F90subroutines/Makefile.gcc_sys 0 → 100644
View file @ 57afb669
#
#
obj = defprec.o ioipsl.o constantes_var.o haversine.o grid.o routing_reg.o routing_tools.o
#
FC = gfortran
FFLAG = -fPIC
#FDBG = -g -fbounds-check
#
FPY = f2py
#PDBG = --debug
#
EXT_SUFFIX := $(shell python3-config --extension-suffix)
#
all : routing_interface$(EXT_SUFFIX) diagst$(EXT_SUFFIX)
routing_interface : routing_interface$(EXT_SUFFIX)
diagst : diagst$(EXT_SUFFIX)
clean :
rm -f compilation.log *.o *.so *.mod routing_interface.so diagst.so *~
routing_interface$(EXT_SUFFIX) : $(obj) routing_interface.f90
rm -f compilation.log
$(FPY) -c routing_interface.f90 $(PDBG) --fcompiler=gfortran -m routing_interface $(obj) > compilation.log
diagst$(EXT_SUFFIX) : diagst.f90
$(FPY) -c diagst.f90 $(PDBG) --fcompiler=gfortran -m diagst > compilation.log
defprec.o : defprec.f90
$(FC) $(FFLAG) $(FDBG) -c defprec.f90
ioipsl.o : ioipsl.f90
$(FC) $(FFLAG) $(FDBG) -c ioipsl.f90
constantes_var.o : constantes_var.f90
$(FC) $(FFLAG) $(FDBG) -c constantes_var.f90
haversine.o : haversine.f90
$(FC) $(FFLAG) $(FDBG) -c haversine.f90
grid.o : defprec.f90 haversine.o constantes_var.o grid.f90
$(FC) $(FFLAG) $(FDBG) -c grid.f90
routing_tools.o : defprec.o grid.o constantes_var.o routing_tools.f90
$(FC) $(FFLAG) $(FDBG) -c routing_tools.f90
routing_reg.o : defprec.f90 grid.o constantes_var.o routing_tools.o routing_reg.f90
$(FC) $(FFLAG) $(FDBG) -c routing_reg.f90
F90subroutines/haversine.f90
View file @ 57afb669
......@@ -57,7 +57,8 @@ MODULE haversine
PRIVATE
PUBLIC :: haversine_polyheadings, haversine_polysort, &
& haversine_polyseglen, haversine_laloseglen, haversine_clockwise, haversine_polyarea, &
& haversine_laloarea, haversine_xyarea, haversine_diffangle, haversine_heading
& haversine_laloarea, haversine_xyarea, haversine_diffangle, haversine_heading, &
& haversine_distance
CONTAINS
!!
......
F90subroutines/routing_interface.f90
View file @ 57afb669
SUBROUTINE initatmgrid(rankmpi, nbcore, nbpt, nbseg, polygons_in, centre_in, area_in, contfrac_in, neighbours_in)
SUBROUTINE initatmgrid(rankmpi, outfile, nbpt, nbseg, polygons_in, centre_in, area_in, contfrac_in, neighbours_in)
!
USE grid
USE ioipsl
......@@ -8,7 +8,7 @@ SUBROUTINE initatmgrid(rankmpi, nbcore, nbpt, nbseg, polygons_in, centre_in, are
!
!! INPUT VARIABLES
INTEGER, INTENT(in) :: rankmpi !! rank of processor
INTEGER, INTENT(in) :: nbcore
CHARACTER(LEN=*) :: outfile
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)
......@@ -30,7 +30,7 @@ SUBROUTINE initatmgrid(rankmpi, nbcore, nbpt, nbseg, polygons_in, centre_in, are
CALL grid_init(nbpt, nbseg, polygons_in, centre_in, area_in, contfrac_in, neighbours_loc)
!
numout=100+rankmpi
WRITE(outfname,'(A15,I4.4,A6,I4.4,A4)') "Out_RoutingReg_",INT(rankmpi),"_from_",INT(nbcore),".txt"
outfname = outfile
OPEN(unit=numout, file=outfname)
!
END SUBROUTINE initatmgrid
......@@ -44,9 +44,9 @@ SUBROUTINE closeinterface()
END SUBROUTINE closeinterface
SUBROUTINE gethydrogrid(nbpt, nbvmax_in, ijdimmax, maxpercent_in, sub_pts, sub_index, sub_area, max_basins, &
& min_topoind, meanrlen, meanrdz, minrdz, lon_rel, lat_rel, trip, basins, topoindex, rlen, rdz, fac, hierarchy, &
& orog, floodp, nbi, nbj, area_bx, trip_bx, basin_bx, topoind_bx, rlen_bx, rdz_bx, rweight_bx, fac_bx, &
& hierarchy_bx, orog_bx, floodp_bx, lon_bx, lat_bx, lshead_bx)
& min_topoind, meanrlen, meanrdz, minrdz, dzprecision, lon_rel, lat_rel, trip, basins, topoindex, rlen, &
& rdz, fac, hierarchy, orog, floodp, nbi, nbj, area_bx, trip_bx, basin_bx, topoind_bx, rlen_bx, rdz_bx, &
& rweight_bx, fac_bx, hierarchy_bx, orog_bx, floodp_bx, lon_bx, lat_bx, lshead_bx)
!
USE ioipsl
USE grid
......@@ -63,6 +63,7 @@ SUBROUTINE gethydrogrid(nbpt, nbvmax_in, ijdimmax, maxpercent_in, sub_pts, sub_i
REAL, INTENT(inout) :: max_basins !! The current maximum of basins
REAL, INTENT(in) :: min_topoind !! The current minimum of topographic index (km)
REAL, INTENT(in) :: meanrlen, meanrdz, minrdz !! Mean length and altitude change in the original grid (m)
REAL, INTENT(in) :: dzprecision !! The precision of the orography on the hydrological data
REAL, INTENT(in) :: sub_area(nbpt,nbvmax_in) !! Area on the fine grid (m^2)
!
REAL, INTENT(in) :: lon_rel(nbpt,nbvmax_in) !!
......@@ -114,11 +115,11 @@ SUBROUTINE gethydrogrid(nbpt, nbvmax_in, ijdimmax, maxpercent_in, sub_pts, sub_i
DO ib=1,nbpt
CALL routing_reg_getgrid(nbpt, ib, ijdimmax, sub_pts, sub_index, sub_area, max_basins, min_topoind, &
& meanrlen, meanrdz, minrdz, lon_rel, lat_rel, lalo, resolution, contfrac, trip, basins, topoindex, &
& rlen, rdz, fac, hierarchy, orog, floodp, nbi(ib), nbj(ib), area_bx(ib,:,:), trip_bx(ib,:,:), &
& basin_bx(ib,:,:), topoind_bx(ib,:,:), rlen_bx(ib,:,:), rdz_bx(ib,:,:), fac_bx(ib,:,:), &
& hierarchy_bx(ib,:,:), orog_bx(ib,:,:), floodp_bx(ib,:,:), lon_bx(ib,:,:), lat_bx(ib,:,:), &
& lshead_bx(ib,:,:))
& meanrlen, meanrdz, minrdz, dzprecision, lon_rel, lat_rel, lalo, resolution, contfrac, trip, &
& basins, topoindex, rlen, rdz, fac, hierarchy, orog, floodp, nbi(ib), nbj(ib), area_bx(ib,:,:), &
& trip_bx(ib,:,:), basin_bx(ib,:,:), topoind_bx(ib,:,:), rlen_bx(ib,:,:), rdz_bx(ib,:,:), &
& fac_bx(ib,:,:), hierarchy_bx(ib,:,:), orog_bx(ib,:,:), floodp_bx(ib,:,:), lon_bx(ib,:,:), &
& lat_bx(ib,:,:), lshead_bx(ib,:,:))
ENDDO
!
rweight_bx(:,:,:) = zero
......@@ -193,8 +194,6 @@ SUBROUTINE findbasins(nbpt, nb_htu, nbv, ijdimmax, nbi, nbj, nbasmax, trip_bx, b
& orog_bx(ib,:,:), area_bx(ib,:,:),basin_topoindex_stream(ib,:), basin_rlen_stream(ib,:), basin_dz_stream(ib,:))
ENDDO
WRITE(numout,*) "XXXX range of nb_basin after findbasin : ", MINVAL(nb_basin), sum(nb_basin)/nbpt, MAXVAL(nb_basin)
END SUBROUTINE findbasins
SUBROUTINE globalize(nbpt, nb_htu, nbv, ijdimmax, area_bx, lon_bx, lat_bx, trip_bx, hierarchy_bx, orog_bx, floodp_bx, &
......@@ -241,10 +240,10 @@ SUBROUTINE globalize(nbpt, nb_htu, nbv, ijdimmax, area_bx, lon_bx, lat_bx, trip_
REAL(r_std), INTENT (in), DIMENSION(nbpt,nb_htu) :: basin_rlen_stream, basin_dz_stream
!! 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 !!
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 !! Coordinates of the outflow point
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 !!
......@@ -286,7 +285,8 @@ SUBROUTINE globalize(nbpt, nb_htu, nbv, ijdimmax, area_bx, lon_bx, lat_bx, trip_
!
END SUBROUTINE globalize
SUBROUTINE linkup(nbpt, ijdimmax, nwbas, inflowmax, basin_count, basin_area, basin_id, basin_flowdir, &
SUBROUTINE linkup(nbpt, ijdimmax, nwbas, inflowmax, basin_count, basin_area, basin_id, basin_outcoor, &
& basin_cg, basin_orog_min, basin_flowdir, &
& basin_lshead, basin_hierarchy, basin_fac, outflow_grid, outflow_basin, inflow_number, inflow_grid, &
& inflow_basin, nbcoastal, coastal_basin, invented_basins)
!
......@@ -303,13 +303,17 @@ SUBROUTINE linkup(nbpt, ijdimmax, nwbas, inflowmax, basin_count, basin_area, bas
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 (in), DIMENSION(nbpt) :: basin_count !!
INTEGER(i_std), INTENT (in), DIMENSION(nbpt,nwbas) :: basin_id !!
REAL(r_std), INTENT (in), DIMENSION(nbpt,nwbas,2) :: basin_outcoor !! Coordinates of the outflow point
REAL(r_std), INTENT (in), DIMENSION(nbpt,nwbas,2) :: basin_cg !! Coordinates of the outflow point
REAL(r_std), INTENT (in), DIMENSION(nbpt,nwbas) :: basin_orog_min !!
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 !!
!
......@@ -325,9 +329,10 @@ SUBROUTINE linkup(nbpt, ijdimmax, nwbas, inflowmax, basin_count, basin_area, bas
!
IF ( debug ) WRITE(numout,*) "Memory Mgt Linkup : nbvmax, ijdimmax, nwbas, inflowmax = ", nbvmax, ijdimmax, nwbas, inflowmax
CALL routing_reg_linkup(nbpt, neighbours, nwbas, ijdimmax, inflowmax, 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)
CALL routing_reg_linkup(nbpt, neighbours, nwbas, ijdimmax, inflowmax, basin_count, basin_area, basin_id, basin_outcoor, &
& basin_cg, basin_orog_min, &
& 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
......@@ -373,7 +378,8 @@ SUBROUTINE areanorm(nbpt, nwbas, basin_count, basin_area, outflow_grid, basin_ar
!
END SUBROUTINE areanorm
SUBROUTINE fetch(nbpt, nwbas, nbcore, nbhalo, fhalopts, corepts, basin_count, basin_area, basin_id, basin_hierarchy, basin_fac, &
SUBROUTINE fetch(nbpt, nwbas, nbcore, nbhalo, fhalopts, corepts, basin_count, basin_area, basin_id,&
& basin_outcoor, basin_hierarchy, basin_fac, &
& outflow_grid, outflow_basin, partial_sum, fetch_basin, outflow_uparea)
!
USE ioipsl
......@@ -393,6 +399,7 @@ SUBROUTINE fetch(nbpt, nwbas, nbcore, nbhalo, fhalopts, corepts, basin_count, ba
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,2), INTENT(in) :: basin_outcoor
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)
......@@ -415,7 +422,7 @@ SUBROUTINE fetch(nbpt, nwbas, nbcore, nbhalo, fhalopts, corepts, basin_count, ba
ENDDO
!
CALL routing_reg_fetch(nbpt, area, contfrac, nwbas, nbcore,nbhalo, fhalopts, &
& fcorepts, basin_count, basin_area, basin_id, &
& fcorepts, basin_count, basin_area, basin_id, basin_outcoor, &
& basin_hierarchy, basin_fac, outflow_grid, outflow_basin, partial_sum, fetch_basin, outflow_uparea)
!
......@@ -919,7 +926,7 @@ END SUBROUTINE correct_inflows
SUBROUTINE get_floodcri(nbpt, nwbas,inflowmax, inflow_number, inflow_basin, inflow_grid, &
& basin_count, basin_floodp, basin_orog_min,floodcri)
& basin_count, basin_floodp, basin_orog_min, lim_floodcri, floodcri)
!
USE ioipsl
USE grid
......@@ -937,12 +944,14 @@ SUBROUTINE get_floodcri(nbpt, nwbas,inflowmax, inflow_number, inflow_basin, infl
INTEGER(i_std), DIMENSION(nbpt,nwbas,inflowmax), INTENT(in) :: inflow_basin
INTEGER(i_std), DIMENSION(nbpt,nwbas,inflowmax), INTENT(in) :: inflow_grid
INTEGER(i_std), INTENT(in) :: lim_floodcri
REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(out) :: floodcri
INTEGER :: ig, ib, ug, ub, i_inf, nbas
REAL :: diff, orog, d
floodcri(:,:) = 2
floodcri(:,:) = lim_floodcri
DO ig=1,nbpt
nbas = basin_count(ig)
......@@ -957,7 +966,7 @@ SUBROUTINE get_floodcri(nbpt, nwbas,inflowmax, inflow_number, inflow_basin, infl
d = basin_orog_min(ug,ub) - orog
IF (d<diff) diff = d
END DO
IF (diff < 2) diff = 2
IF (diff < lim_floodcri) diff = lim_floodcri
floodcri(ig, ib) = diff
END IF
END IF
......@@ -1072,7 +1081,7 @@ SUBROUTINE finalfetch(nbpt, nbasmax, nbcore, nbhalo, in_core, halopts, corepts,
END IF
itt = itt + 1
IF ( itt .GT. maxiterations) THEN
CALL ipslerr_p(3,'routing_reg_fetch','Problem in propagating partial sum','','')
CALL ipslerr_p(3,'routing_finalfetch','Problem in propagating partial sum','','')
ENDIF
ENDDO
END IF
......
F90subroutines/routing_reg.f90
View file @ 57afb669
This diff is collapsed.
F90subroutines/routing_tools.f90
View file @ 57afb669
......@@ -3,6 +3,7 @@ MODULE routing_tools
USE ioipsl
USE constantes_var
USE grid
USE haversine
!
IMPLICIT NONE
!
......@@ -203,18 +204,20 @@ CONTAINS
!! \n
!_ ================================================================================================================================
!
INTEGER(i_std) FUNCTION routing_nextgrid(pnt, trip, inc)
INTEGER(i_std) FUNCTION routing_nextgrid(pnt, trip, ilon, ilat, inc)
!
! Input
!
INTEGER(i_std) :: pnt, trip
INTEGER(i_std), OPTIONAL :: inc
REAL(r_std) :: ilon, ilat
!
! Local
!
INTEGER(i_std) :: i
REAL(r_std), DIMENSION(NbNeighb) :: diffangle
INTEGER(i_std), DIMENSION(1) :: imin
REAL(r_std) :: angle, dangle1, dangle2
INTEGER(i_std), DIMENSION(1) :: imin, imin1, imin2, imin3
INTEGER(i_std) :: trip_loc, nbind
!
! Make sure that the flow direction coming from routing.f90 is between 1 and 8. This allows to call the
......@@ -225,18 +228,60 @@ CONTAINS
! Compute the angle between the direction of the routing.nc flow and the headings perpendicular of the segments
! of the polygon.
!
! lalo is (lat, lon)
angle = MOD(haversine_heading(ilon, ilat, lalo_g(pnt,2), lalo_g(pnt,1))+180.0, 360.0)
DO i=1,NbNeighb
diffangle(i) = ABS(haversine_diffangle(headings(pnt,i), rose(trip_loc)))
diffangle(i) = ABS(haversine_diffangle(headings_g(pnt,i), angle))
ENDDO
imin1 = MINLOC(diffangle)
dangle1 = diffangle(imin1(1))
diffangle(imin1(1)) = 999
imin2 = MINLOC(diffangle)
dangle2 = diffangle(imin2(1))
diffangle(imin2(1)) = 999
imin3 = MINLOC(diffangle)
! See if the closest heading is a corner (even number)
IF (MOD(imin1(1),2) .EQ. 0) THEN
! We may consider this corner only if the angle is inferior to 1
IF (dangle1 .LT. 1) THEN
! If we are in the corner (for example the pixel is between W and SW)
IF (trip_loc .EQ. imin1(1)) THEN
! it may follow the diagonal (SW)
nbind = trip_loc
ELSE IF (trip_loc .EQ. imin3(1)) THEN
! The pixel may also go to the other direction (S)
! We may miss the case where the direction is SE and it also go to
! the S direction, this should be improved later
nbind = imin3(1)
ELSE
! the other case is that it go to the grid point of the second
! closest angle (here W)
nbind = imin2(1)
END IF
ELSE
nbind = imin2(1)
END IF
ELSE
nbind = imin1(1)
END IF
! Previous version of the code
! Is also present in the routing_nextgrid_g
!DO i=1,NbNeighb
! diffangle(i) = ABS(haversine_diffangle(headings(pnt,i), rose(trip_loc)))
!ENDDO
!
imin = MINLOC(diffangle)
!imin = MINLOC(diffangle)
!
IF ( PRESENT(inc) ) THEN
! Ensure that the increment the users asks is within the NbNeighb the polygone has.
nbind = MODULO(imin(1)+inc-1, NbNeighb)+1
ELSE
nbind = imin(1)
ENDIF
!IF ( PRESENT(inc) ) THEN
! ! Ensure that the increment the users asks is within the NbNeighb the polygone has.
! nbind = MODULO(imin(1)+inc-1, NbNeighb)+1
!ELSE
! nbind = imin(1)
!ENDIF
!
routing_nextgrid = neighbours(pnt,nbind)
!
......@@ -307,14 +352,17 @@ CONTAINS
!! \n
!_ ================================================================================================================================
!
SUBROUTINE routing_reg_bestsubbasin(src_grid, src_basin, src_id, src_hierarchy, src_fac, src_flowdir, invented_basins, &
& ptmax, basmax, pt, nbsubbas, ids, hierarchy, fac, flowdir, nbcoastals, coastals, &
& outbas, quality)
SUBROUTINE routing_reg_bestsubbasin(src_grid, src_basin, src_id, src_hierarchy, src_fac, src_flowdir, src_outcoord, src_orog,&
& invented_basins, ptmax, basmax, pt, nbsubbas, ids, hierarchy, fac, flowdir, outcoord, &
& out_cg, out_orog_min, &
& nbcoastals, coastals, outbas, quality)
!
! Input
!
INTEGER(i_std), INTENT(in) :: src_grid, src_id, src_flowdir, src_basin
REAL(r_std), INTENT(in) :: src_hierarchy, src_fac ! src_fac should be integer but nevermind
REAL(r_std), INTENT(in) :: src_outcoord(2)
REAL(r_std), INTENT(in) :: src_orog
REAL(r_std), INTENT(in) :: invented_basins
!
INTEGER(i_std), INTENT(in) :: ptmax, basmax, pt
......@@ -322,6 +370,9 @@ CONTAINS
INTEGER(i_std), INTENT(in) :: ids(ptmax,basmax), flowdir(ptmax,basmax)
REAL(r_std), INTENT(in) :: hierarchy(ptmax,basmax)
REAL(r_std), INTENT(in) :: fac(ptmax,basmax) ! fac should be integer but real is still ok
REAL(r_std), INTENT(in) :: outcoord(ptmax,basmax,2) ! Coordinates of outflow point
REAL(r_std), INTENT(in) :: out_cg(ptmax,basmax,2)
REAL(r_std), INTENT(in) :: out_orog_min(ptmax,basmax)
INTEGER(i_std), INTENT(in) :: nbcoastals(ptmax)
INTEGER(i_std), INTENT(in) :: coastals(ptmax,basmax)
!
......@@ -333,17 +384,24 @@ CONTAINS
! Local
!
INTEGER(i_std) :: i, j, sbl, nbfbas, nbopt
INTEGER(i_std), DIMENSION(1) :: ff
INTEGER(i_std) :: ih, id
INTEGER(i_std), DIMENSION(1) :: ff, gg
INTEGER(i_std), DIMENSION(basmax) :: fbas, tbas, sboptions, angle, subbas
INTEGER(i_std), DIMENSION(basmax) :: fbas_flowdir, tbas_flowdir
REAL(r_std), DIMENSION(basmax) :: fbas_hierarchy, tbas_hierarchy
REAL(r_std), DIMENSION(basmax) :: fbas_options
REAL(r_std), DIMENSION(basmax) :: fbas_fac, tbas_fac
REAL(r_std), DIMENSION(basmax) :: fbas_dist, tbas_dist
REAL(r_std), DIMENSION(basmax) :: fbas_dist2, tbas_dist2
REAL(r_std), DIMENSION(basmax) :: fbas_dz, tbas_dz
!
! Debug
!
LOGICAL, PARAMETER :: debug = .FALSE. !! Logical variable to debug
LOGICAL :: follow
LOGICAL, PARAMETER :: verifchoice = .FALSE.
INTEGER(i_std) :: testbasinid
REAL(r_std) :: maxhiediff
REAL(r_std) :: maxhiediff, dist
!
! Set some default values
!
......@@ -370,6 +428,14 @@ CONTAINS
tbas_hierarchy(nbfbas) = hierarchy(pt,sbl)
tbas_fac(nbfbas) = fac(pt,sbl)
tbas_flowdir(nbfbas) = flowdir(pt,sbl)
! Distance between outflow and outflow
tbas_dist(nbfbas) = haversine_distance(src_outcoord(2), src_outcoord(1), &
& outcoord(pt,sbl,2), outcoord(pt,sbl,1))
! Distance between outflow and center of gravity
tbas_dist2(nbfbas) = haversine_distance(src_outcoord(2), src_outcoord(1), &
& out_cg(pt,sbl,2), out_cg(pt,sbl,1))
! Difference of dz between outflow
tbas_dz(nbfbas) = src_orog - out_orog_min(pt,sbl)
ENDIF
ELSE
IF ( COUNT(coastals(pt,1:nbcoastals(pt)) .EQ. src_id) .GT. 0 ) THEN
......@@ -377,6 +443,14 @@ CONTAINS
tbas(nbfbas) = sbl
tbas_hierarchy(nbfbas) = hierarchy(pt,sbl)
tbas_fac(nbfbas) = fac(pt,sbl)
! Distance between outflow and outflow
tbas_dist(nbfbas) = haversine_distance(src_outcoord(2), src_outcoord(1), &
& outcoord(pt,sbl,2), outcoord(pt,sbl,1))
! Distance between outflow and center of gravity
tbas_dist2(nbfbas) = haversine_distance(src_outcoord(2), src_outcoord(1), &
& out_cg(pt,sbl,2), out_cg(pt,sbl,1))
! Difference of dz between outflow
tbas_dz(nbfbas) = src_orog - out_orog_min(pt,sbl)
ENDIF