diff --git a/CMakeLists.txt b/CMakeLists.txt index ad9c5b2..c2ec9af 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -13,6 +13,7 @@ set(FASTSCAPELIB_SRC_FILES ${FASTSCAPELIB_SRC_DIR}/Diffusion.f90 ${FASTSCAPELIB_SRC_DIR}/FastScape_api.f90 ${FASTSCAPELIB_SRC_DIR}/FastScape_ctx.f90 + ${FASTSCAPELIB_SRC_DIR}/FastScapeErrorCodes.f90 ${FASTSCAPELIB_SRC_DIR}/LocalMinima.f90 ${FASTSCAPELIB_SRC_DIR}/Marine.f90 ${FASTSCAPELIB_SRC_DIR}/Strati.f90 @@ -39,14 +40,19 @@ if(WIN32) add_compile_definitions(ON_WINDOWS) endif() +if(SKBUILD) + add_compile_definitions(NO_ERROR_MSG) +endif() + set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -cpp") # Flags # ===== -if(CMAKE_Fortran_COMPILER_ID MATCHES "Flang|GNU") - set(dialect "-ffree-form -std=f2008 -fimplicit-none -fall-intrinsics") +if(CMAKE_Fortran_COMPILER_ID MATCHES "GNU") + set(dialect "-ffree-form -std=f2018 -fimplicit-none -fall-intrinsics -ffree-line-length-none") set(ioflags "-fconvert=big-endian") + set(pythonflags "-ffree-line-length-none") set(f77flags "-std=legacy -ffixed-form") set(bounds "-fbounds-check") set(warnings "-Wall") @@ -58,9 +64,10 @@ if(CMAKE_Fortran_COMPILER_ID MATCHES "Intel") set(bounds "-check bounds") set(warnings "-warn all") endif() -if(CMAKE_Fortran_COMPILER_ID MATCHES "PGI") +if(CMAKE_Fortran_COMPILER_ID MATCHES "Flang|PGI") set(dialect "-Mfreeform -Mdclchk -Mstandard -Mallocatable=03") set(ioflags "") + set(pythonflags "-Mextend") set(f77flags "-Mfixed") set(bounds "-C") set(warnings "-Wall") @@ -69,8 +76,11 @@ endif() set(CMAKE_Fortran_FLAGS_DEBUG "${CMAKE_Fortran_FLAGS_DEBUG} ${bounds} ${warnings}") -# let F2PY (scikit-build) configure flags for F77 and F90 source -if(NOT SKBUILD) +# let F2PY (scikit-build) configure dialect flags for F77 and F90 source +# (some source files are generated by F2PY) +if(SKBUILD) + set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} ${pythonflags}") +else() set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} ${dialect} ${ioflags}") endif() @@ -99,9 +109,15 @@ endif() # Fortran library (static/shared) # =============================== +set(CMAKE_Fortran_MODULE_DIRECTORY ${CMAKE_BINARY_DIR}/mod) + if(BUILD_FASTSCAPELIB_STATIC OR BUILD_FASTSCAPELIB_SHARED) + # build temp library set(FASTSCAPELIB_OBJECTS libfastscape_objects) add_library(${FASTSCAPELIB_OBJECTS} OBJECT ${FASTSCAPELIB_SRC_FILES}) + + # install api Fortran module + install(FILES ${CMAKE_Fortran_MODULE_DIRECTORY}/fastscapeapi.mod DESTINATION include) endif() if(BUILD_FASTSCAPELIB_STATIC) diff --git a/examples/DippingDyke.f90 b/examples/DippingDyke.f90 index 46edbb2..97ab474 100644 --- a/examples/DippingDyke.f90 +++ b/examples/DippingDyke.f90 @@ -5,6 +5,8 @@ program DippingDyke ! and is progressively exhumed by erosion; for this we use the total erosion ! to define the erodibility array kf + use FastScapeAPI + implicit none integer :: nx, ny, istep, nstep, i, j diff --git a/examples/Fan.f90 b/examples/Fan.f90 index 29a47f5..3c09e2f 100644 --- a/examples/Fan.f90 +++ b/examples/Fan.f90 @@ -11,9 +11,11 @@ program Fan ! The evolution of the sedimentary flux out of the system and out of the ! plateau only is stored in the Fluxes.txt file + use FastScapeAPI + implicit none - integer :: nx,ny,istep,nstep,nn,ibc + integer :: nx,ny,istep,nstep,nn,ibc,ierr double precision, dimension(:), allocatable :: h,hp,x,y,kf,kd,b real :: time_in,time_out double precision :: kfsed,m,n,kdsed,g1,g2,expp @@ -30,13 +32,13 @@ program Fan ! initialize FastScape call FastScape_Init () - call FastScape_Set_NX_NY (nx,ny) - call FastScape_Setup () + call FastScape_Set_NX_NY (nx,ny,ierr) + call FastScape_Setup (ierr) ! set model dimensions xl=10.d3 yl=20.d3 - call FastScape_Set_XL_YL (xl,yl) + call FastScape_Set_XL_YL (xl,yl,ierr) ! construct nodal coordinate arrays x and y allocate (x(nx*ny),y(nx*ny)) @@ -45,7 +47,7 @@ program Fan ! set time step dt=2.d3 - call FastScape_Set_DT (dt) + call FastScape_Set_DT (dt,ierr) ! we make the sediment slightly more easily erodible allocate (kf(nn),kd(nn)) @@ -58,26 +60,26 @@ program Fan g1=1.d0 g2=1.d0 expp=1.d0 - call FastScape_Set_Erosional_Parameters (kf,kfsed,m,n,kd,kdsed,g1,g2,expp) + call FastScape_Set_Erosional_Parameters (kf,kfsed,m,n,kd,kdsed,g1,g2,expp,ierr) ! bottom side is fixed only ibc=1000 - call FastScape_Set_BC (ibc) + call FastScape_Set_BC (ibc,ierr) ! initial topography is a 1000 m high plateau allocate (h(nn),b(nn),hp(nn)) call random_number (h) where (y.gt.yl/2.d0) h=h+1000.d0 - call FastScape_Init_H (h) + call FastScape_Init_H (h, ierr) ! set number of time steps nstep = 200 ! echo model setup - call FastScape_View () + call FastScape_View (ierr) ! initializes time step - call FastScape_Get_Step (istep) + call FastScape_Get_Step (istep,ierr) ! set vertical exaggeration vex = 3.d0 @@ -87,23 +89,23 @@ program Fan do while (istep.lt.nstep) ! execute FastScape step - call FastScape_Execute_Step () - call FastScape_Get_Step (istep) + call FastScape_Execute_Step (ierr) + call FastScape_Get_Step (istep,ierr) ! output vtk with sediment thickness information - call FastScape_Copy_H (h) - call FastScape_Copy_Basement (b) - call FastScape_VTK (h-b, vex) + call FastScape_Copy_H (h,ierr) + call FastScape_Copy_Basement (b,ierr) + call FastScape_VTK (h-b, vex, ierr) enddo ! display timing information - call FastScape_Debug() + call FastScape_Debug(ierr) call cpu_time (time_out) print*,'Total run time',time_out-time_in ! exits FastScape - call FastScape_Destroy () + call FastScape_Destroy (ierr) ! deallocate memory deallocate (h,x,y,kf,kd,b,hp) diff --git a/examples/Margin.f90 b/examples/Margin.f90 index afc5c9c..0de1b9f 100644 --- a/examples/Margin.f90 +++ b/examples/Margin.f90 @@ -13,6 +13,8 @@ program Margin ! of the nodes on the grid; these are used to define the uplift function and ! the initial topography + use FastScapeAPI + implicit none integer :: nx, ny, istep, nstep, nfreq, i, j diff --git a/examples/Mountain.f90 b/examples/Mountain.f90 index c0b1b2c..2d0363a 100644 --- a/examples/Mountain.f90 +++ b/examples/Mountain.f90 @@ -8,6 +8,8 @@ program Mountain ! nonlinear erosion law (n=1.5, m=0.6) ! transport coefficient g = 1 + use FastScapeAPI + implicit none integer :: nx, ny, istep, nstep diff --git a/examples/Strati_test.f90 b/examples/Strati_test.f90 index 975bfb6..e770267 100644 --- a/examples/Strati_test.f90 +++ b/examples/Strati_test.f90 @@ -13,6 +13,8 @@ program Strati_test ! of the nodes on the grid; these are used to define the uplift function and ! the initial topography + use FastScapeAPI + implicit none integer :: nx, ny, istep, nstep, nfreq, i, j, nreflector diff --git a/examples/flexure_test.f90 b/examples/flexure_test.f90 index df5ff31..b70e916 100644 --- a/examples/flexure_test.f90 +++ b/examples/flexure_test.f90 @@ -2,6 +2,7 @@ program flexure_test ! test example to demonstrate how to use the flexure routine ! see relevant section below for explanations + use FastScapeAPI implicit none diff --git a/src/Diffusion.f90 b/src/Diffusion.f90 index 6b1c2fa..3c3af64 100644 --- a/src/Diffusion.f90 +++ b/src/Diffusion.f90 @@ -1,4 +1,5 @@ -subroutine Diffusion () +#include "Error.fpp" +subroutine Diffusion (ierr) ! subroutine to solve the diffusion equation by ADI @@ -11,6 +12,7 @@ subroutine Diffusion () integer i,j,ij double precision factxp,factxm,factyp,factym,dx,dy character cbc*4 + integer ierr !print*,'Diffusion' @@ -79,7 +81,7 @@ subroutine Diffusion () inf(nx)=-factxm f(nx)=zintp(nx,j)+factyp*zintp(nx,j+1)-(factyp+factym)*zintp(nx,j)+factym*zintp(nx,j-1) endif - call tridag (inf,diag,sup,f,res,nx) + call tridag (inf,diag,sup,f,res,nx,ierr);FSCAPE_CHKERR(ierr) do i=1,nx zint(i,j)=res(i) enddo @@ -131,7 +133,7 @@ subroutine Diffusion () inf(ny)=-factym f(ny)=zint(i,ny)+factxp*zint(i+1,ny)-(factxp+factxm)*zint(i,ny)+factxm*zint(i-1,ny) endif - call tridag (inf,diag,sup,f,res,ny) + call tridag (inf,diag,sup,f,res,ny,ierr);FSCAPE_CHKERR(ierr) do j=1,ny zintp(i,j)=res(j) enddo @@ -161,8 +163,9 @@ end subroutine Diffusion ! subroutine to solve a tri-diagonal system of equations (from Numerical Recipes) - SUBROUTINE tridag(a,b,c,r,u,n) + SUBROUTINE tridag(a,b,c,r,u,n,ierr) + use FastScapeErrorCodes implicit none INTEGER n @@ -170,10 +173,13 @@ SUBROUTINE tridag(a,b,c,r,u,n) INTEGER j double precision bet double precision,dimension(:),allocatable::gam - + integer,intent(inout):: ierr allocate (gam(n)) - if(b(1).eq.0.d0) stop 'in tridag' + if(b(1).eq.0.d0) then + FSCAPE_RAISE_MESSAGE('Tridag error b(1)=0: cannot solve system',ERR_TridiagNotSolvable,ierr);FSCAPE_CHKERR(ierr) + end if + !if(b(1).eq.0.d0) stop 'in tridag' ! first pass @@ -183,8 +189,9 @@ SUBROUTINE tridag(a,b,c,r,u,n) gam(j)=c(j-1)/bet bet=b(j)-a(j)*gam(j) if(bet.eq.0.) then - print*,'tridag failed' - stop + FSCAPE_RAISE_MESSAGE('Tridag error bet=0: cannot solve system',ERR_TridiagNotSolvable,ierr);FSCAPE_CHKERR(ierr) + ! print*,'tridag failed' + ! stop endif u(j)=(r(j)-a(j)*u(j-1))/bet 11 continue diff --git a/src/Error.fpp b/src/Error.fpp new file mode 100644 index 0000000..62a1850 --- /dev/null +++ b/src/Error.fpp @@ -0,0 +1,99 @@ + +! Error.fpp + +! +! Macros for error handling. +! Fortran preprocessor must be enabled: -fpp. +! Note: For Python bindings no error msg is shown (error codes should be handled in Python) +! + +! +! Initialize ierr argument in API functions +! +#ifdef NO_ERROR_MSG +#define FSCAPE_INITERR(ierr, ierr_, fname) \ + ierr_ = 0; \ + if (present(ierr)) ierr = ierr_; +#else +#define FSCAPE_INITERR(ierr, ierr_, fname) \ + ierr_ = 0; \ + if (present(ierr)) then; \ + ierr = ierr_; \ + else; \ + print '(A)', "*** FastScape depreciation warning *** "; \ + print '(A)', "When calling " // TRIM((fname)); \ + print '(A)', "Calling API functions without 'ierr' argument (integer) is depreciated! Please update your code!"; \ + end if; +#endif + +! +! Raise an exception by pushing the name of the exception + file name, line number to stdout +! Sets ierr to the error type specified by err_type +! +#ifdef NO_ERROR_MSG +#define FSCAPE_RAISE(err_type, ierr) ierr = (err_type); +#else +#define FSCAPE_RAISE(err_type, ierr) \ + ierr = (err_type); \ + print '(A)', "*** FastScape exception *** -> " // TRIM(err_names((ierr))); \ + print '(A,I4,A)', "In file " // __FILE__ // ", line", __LINE__; +#endif + +! +! Raise an exception by pushing custom user provided message to stdout, +! along with the name of the exception + file name, line number to stdout +! Sets ierr to the error type specified by err_type +! +#ifdef NO_ERROR_MSG +#define FSCAPE_RAISE_MESSAGE(msg, err_type, ierr) FSCAPE_RAISE(err_type, ierr); +#else +#define FSCAPE_RAISE_MESSAGE(msg, err_type, ierr) \ + print '(A)', "*** FastScape exception *** "; \ + print '(A)', TRIM((msg)) // "!"; \ + FSCAPE_RAISE(err_type, ierr); +#endif + +! +! Calls return if error code is non-zero. +! +#define FSCAPE_CHKERR(ierr) \ + if (ierr /= 0) then; \ + return; \ + end if; + +! +! Either calls return or stop if error code is non-zero, dependening on whether the +! the ierr argument is present. +! +#ifdef NO_ERROR_MSG +#define FSCAPE_CHKERR_OPT(ierr, ierr_) \ + if (ierr_ /= 0) then; \ + if (present(ierr)) then; \ + ierr = ierr_; \ + return; \ + else; \ + stop; \ + end if; \ + end if; +#else +#define FSCAPE_CHKERR_OPT(ierr, ierr_) \ + if (ierr_ /= 0) then; \ + if (present(ierr)) then; \ + ierr = ierr_; \ + return; \ + else; \ + print '(A)', "*** Execution being halted by calling stop!"; \ + print '(A)', "Call all API routine with 'ierr' argument to avoid this."; \ + stop; \ + end if; \ + end if; +#endif + +! +! Force stop if error code is non-zero. +! This should only be called in driver code - never within library code or the API. +! +#define FSCAPE_CHKERR_ABORT(ierr) \ + if (ierr /= 0) then; \ + stop; \ + end if; diff --git a/src/FastScapeErrorCodes.f90 b/src/FastScapeErrorCodes.f90 new file mode 100644 index 0000000..b99511f --- /dev/null +++ b/src/FastScapeErrorCodes.f90 @@ -0,0 +1,23 @@ + +module FastScapeErrorCodes + implicit none + integer, parameter :: ERR_None = 0, & + ERR_Default = 1, & + ERR_FileNotFound = 2, & + ERR_FileOpenFailed = 3, & + ERR_ParameterInvalid = 4, & + ERR_ParameterOutOfRange = 5, & + ERR_TridiagNotSolvable = 6, & + ERR_SetupNotRun = 7, & + ERR_NotConverged = 8 + + character(len=50), dimension(8) :: err_names = [character(len=50) :: "Default run time error", & + "File error: File not found", & + "File error: File open failed", & + "Parameter error: Input invalid", & + "Parameter error: Out of range" , & + "Solver error: Tridiag not solvable", & + "Order invalid: Must call SetUp() first", & + "Solver error: not converged" ] + +end module FastScapeErrorCodes diff --git a/src/FastScape_api.f90 b/src/FastScape_api.f90 index 8455cea..e6af7e2 100644 --- a/src/FastScape_api.f90 +++ b/src/FastScape_api.f90 @@ -142,633 +142,860 @@ ! writes debugging information to the default output ! ----------------------------------------------------------------------------------------- +#include "Error.fpp" -subroutine FastScape_Init() +module FastScapeAPI - use FastScapeContext + use iso_c_binding - implicit none + contains - call Init() + subroutine FastScape_Init(ierr) bind(C, name='fastscape_init') - return + use FastScapeContext + implicit none -end subroutine FastScape_Init + integer(c_int), optional, intent(out) :: ierr + integer :: ierr_ -!-------------------------------------------------------------------------- + FSCAPE_INITERR(ierr, ierr_, 'FastScape_Init()') -subroutine FastScape_Setup() + call Init() - use FastScapeContext + !return - implicit none + end subroutine FastScape_Init - call SetUp() + !-------------------------------------------------------------------------- - return + subroutine FastScape_Setup(ierr) bind(C, name='fastscape_setup') -end subroutine FastScape_Setup + use FastScapeContext + implicit none -!-------------------------------------------------------------------------- + integer(c_int), optional, intent(out) :: ierr + integer :: ierr_ -subroutine FastScape_Destroy() + FSCAPE_INITERR(ierr, ierr_, 'FastScape_Setup()') - use FastScapeContext + if (nx.eq.0) then + FSCAPE_RAISE_MESSAGE('FastScape_Setup(): nx cannot be zero',ERR_ParameterInvalid,ierr_) + end if + if (nx.le.0) then + FSCAPE_RAISE_MESSAGE('FastScape_Setup(): nx cannot be negative',ERR_ParameterOutOfRange,ierr_) + end if + if (ny.eq.0) then + FSCAPE_RAISE_MESSAGE('FastScape_Setup(): ny cannot be zero',ERR_ParameterInvalid,ierr_) + end if + if (ny.le.0) then + FSCAPE_RAISE_MESSAGE('FastScape_Setup(): ny cannot be negative',ERR_ParameterOutOfRange,ierr_) + end if + FSCAPE_CHKERR_OPT(ierr, ierr_) ! Call FSCAPE_CHKERR() so that all possible exceptions above will be displayed - implicit none + call SetUp() - call Destroy() + return - return + end subroutine FastScape_Setup -end subroutine FastScape_Destroy + !-------------------------------------------------------------------------- -!-------------------------------------------------------------------------- + subroutine FastScape_Destroy(ierr) bind(C, name='fastscape_destroy') -subroutine FastScape_View() + use FastScapeContext - use FastScapeContext + implicit none - implicit none + integer(c_int), optional, intent(out) :: ierr + integer :: ierr_ - call View() + FSCAPE_INITERR(ierr, ierr_, 'FastScape_Destroy()') - return + call Destroy() -end subroutine FastScape_View + return -!-------------------------------------------------------------------------- -subroutine FastScape_Execute_Step() + end subroutine FastScape_Destroy - use FastScapeContext + !-------------------------------------------------------------------------- - implicit none + subroutine FastScape_View(ierr) bind(C, name='fastscape_view') - real :: time_in, time_out + use FastScapeContext - if (runAdvect) then - call cpu_time (time_in) - call Advect () - call cpu_time (time_out) - timeAdvect = timeAdvect + time_out-time_in - endif + implicit none - if (runUplift) then - call cpu_time (time_in) - call Uplift() - call cpu_time (time_out) - timeUplift = timeUplift + time_out-time_in - endif + integer(c_int), optional, intent(out) :: ierr + integer :: ierr_ - if (runSPL) then - call cpu_time (time_in) - if (SingleFlowDirection) then - call FlowRoutingSingleFlowDirection () - call FlowAccumulationSingleFlowDirection () - call StreamPowerLawSingleFlowDirection () - else - call FlowRouting () - call FlowAccumulation () - call StreamPowerLaw () + FSCAPE_INITERR(ierr, ierr_, 'FastScape_View()') + + call View() + + return + + end subroutine FastScape_View + + !-------------------------------------------------------------------------- + subroutine FastScape_Execute_Step(ierr) bind(C, name='fastscape_execute_step') + + use FastScapeContext + + implicit none + + real :: time_in, time_out + integer(c_int), optional, intent(out) :: ierr + integer :: ierr_ + + FSCAPE_INITERR(ierr, ierr_, 'FastScape_Execute_Step()') + + if (runAdvect) then + call cpu_time (time_in) + call Advect () + call cpu_time (time_out) + timeAdvect = timeAdvect + time_out-time_in + endif + + if (runUplift) then + call cpu_time (time_in) + call Uplift() + call cpu_time (time_out) + timeUplift = timeUplift + time_out-time_in + endif + + if (runSPL) then + call cpu_time (time_in) + if (SingleFlowDirection) then + call FlowRoutingSingleFlowDirection (ierr_);FSCAPE_CHKERR_OPT(ierr, ierr_) + call FlowAccumulationSingleFlowDirection () + call StreamPowerLawSingleFlowDirection () + else + call FlowRouting (ierr_);FSCAPE_CHKERR_OPT(ierr, ierr_) + call FlowAccumulation () + call StreamPowerLaw () + endif + call cpu_time (time_out) + timeSPL = timeSPL + time_out-time_in + endif + + if (runDiffusion) then + call cpu_time (time_in) + call Diffusion (ierr_);FSCAPE_CHKERR_OPT(ierr, ierr_) + call cpu_time (time_out) + timeDiffusion = timeDiffusion + time_out-time_in + endif + + if (runMarine) then + call cpu_time (time_in) + call Marine (ierr_);FSCAPE_CHKERR_OPT(ierr, ierr_) + call cpu_time (time_out) + timeMarine = timeMarine + time_out-time_in + endif + + if (runStrati) then + call cpu_time (time_in) + call Run_Strati () + call cpu_time (time_out) + timeStrati = timeStrati + time_out-time_in endif - call cpu_time (time_out) - timeSPL = timeSPL + time_out-time_in - endif - if (runDiffusion) then - call cpu_time (time_in) - call Diffusion () - call cpu_time (time_out) - timeDiffusion = timeDiffusion + time_out-time_in - endif + step=step+1 + + return + + end subroutine FastScape_Execute_Step + + !-------------------------------------------------------------------------- + + subroutine FastScape_Init_H(hp,ierr) bind(C, name='fastscape_init_h') + + use FastScapeContext + + implicit none + + real(c_double), intent(inout), dimension(*) :: hp + integer(c_int), optional, intent(out) :: ierr + integer :: ierr_ + + FSCAPE_INITERR(ierr, ierr_, 'FastScape_Init_H()') + + if (.not.setup_has_been_run) then + FSCAPE_RAISE(ERR_SetupNotRun,ierr_);FSCAPE_CHKERR_OPT(ierr, ierr_) + end if + + call InitH(hp) + + return + + end subroutine FastScape_Init_H + + !-------------------------------------------------------------------------- + + subroutine FastScape_Init_F(Fmixp,ierr) bind(C, name='fastscape_init_f') + + use FastScapeContext + + implicit none + + real(c_double), intent(inout), dimension(*) :: Fmixp + integer(c_int), optional, intent(out) :: ierr + integer :: ierr_ + + FSCAPE_INITERR(ierr, ierr_, 'FastScape_Init_F()') + + if (.not.setup_has_been_run) then + FSCAPE_RAISE(ERR_SetupNotRun,ierr_);FSCAPE_CHKERR_OPT(ierr, ierr_) + end if + + call InitF (Fmixp) + + return + + end subroutine FastScape_Init_F + + !-------------------------------------------------------------------------- + + subroutine FastScape_Copy_H(hp,ierr) bind(C, name='fastscape_copy_h') + + use FastScapeContext + + implicit none + + real(c_double), intent(inout), dimension(*) :: hp + integer(c_int), optional, intent(out) :: ierr + integer :: ierr_ + + FSCAPE_INITERR(ierr, ierr_, 'FastScape_Copy_H()') + + if (.not.setup_has_been_run) then + FSCAPE_RAISE(ERR_SetupNotRun,ierr_);FSCAPE_CHKERR_OPT(ierr, ierr_) + end if + + call CopyH(hp) + + return + + end subroutine FastScape_Copy_H + + !-------------------------------------------------------------------------- + + subroutine FastScape_Copy_Basement(bp,ierr) bind(C, name='fastscape_copy_basement') + + use FastScapeContext + + implicit none + + real(c_double), intent(inout), dimension(*) :: bp + integer(c_int), optional, intent(out) :: ierr + integer :: ierr_ + + FSCAPE_INITERR(ierr, ierr_, 'FastScape_Copy_Basement()') + + if (.not.setup_has_been_run) then + FSCAPE_RAISE(ERR_SetupNotRun,ierr_);FSCAPE_CHKERR_OPT(ierr, ierr_) + end if + + call CopyBasement(bp) + + return + + end subroutine FastScape_Copy_Basement + + !-------------------------------------------------------------------------- + + subroutine FastScape_Copy_Total_Erosion (etotp,ierr) bind(C, name='fastscape_copy_total_erosion') + + use FastScapeContext - if (runMarine) then - call cpu_time (time_in) - call Marine () - call cpu_time (time_out) - timeMarine = timeMarine + time_out-time_in - endif + implicit none - if (runStrati) then - call cpu_time (time_in) - call Run_Strati () - call cpu_time (time_out) - timeStrati = timeStrati + time_out-time_in - endif + real(c_double), intent(inout), dimension(*) :: etotp + integer(c_int), optional, intent(out) :: ierr + integer :: ierr_ - step=step+1 + FSCAPE_INITERR(ierr, ierr_, 'FastScape_Copy_Total_Erosion()') - return + if (.not.setup_has_been_run) then + FSCAPE_RAISE(ERR_SetupNotRun,ierr_);FSCAPE_CHKERR_OPT(ierr, ierr_) + end if -end subroutine FastScape_Execute_Step + call CopyEtot(etotp) -!-------------------------------------------------------------------------- + return -subroutine FastScape_Init_H(hp) + end subroutine FastScape_Copy_Total_Erosion - use FastScapeContext + !-------------------------------------------------------------------------- - implicit none + subroutine FastScape_Copy_Drainage_Area (ap,ierr) bind(C, name='fastscape_copy_drainage_area') - double precision, intent(inout), dimension(*) :: hp + use FastScapeContext - call InitH(hp) + implicit none - return + real(c_double), intent(inout), dimension(*) :: ap + integer(c_int), optional, intent(out) :: ierr + integer :: ierr_ -end subroutine FastScape_Init_H + FSCAPE_INITERR(ierr, ierr_, 'FastScape_Copy_Drainage_Area()') -!-------------------------------------------------------------------------- + if (.not.setup_has_been_run) then + FSCAPE_RAISE(ERR_SetupNotRun,ierr_);FSCAPE_CHKERR_OPT(ierr, ierr_) + end if -subroutine FastScape_Init_F(Fmixp) + call CopyArea(ap) - use FastScapeContext + return - implicit none + end subroutine FastScape_Copy_Drainage_Area - double precision, intent(inout), dimension(*) :: Fmixp + !-------------------------------------------------------------------------- - call InitF (Fmixp) + subroutine FastScape_Copy_Erosion_Rate (eratep,ierr) bind(C, name='fastscape_copy_erosion_rate') - return + use FastScapeContext -end subroutine FastScape_Init_F + implicit none -!-------------------------------------------------------------------------- + real(c_double), intent(inout), dimension(*) :: eratep + integer(c_int), optional, intent(out) :: ierr + integer :: ierr_ -subroutine FastScape_Copy_H(hp) + FSCAPE_INITERR(ierr, ierr_, 'FastScape_Copy_Erosion_Rate()') - use FastScapeContext + if (.not.setup_has_been_run) then + FSCAPE_RAISE(ERR_SetupNotRun,ierr_);FSCAPE_CHKERR_OPT(ierr, ierr_) + end if - implicit none + call CopyERate(eratep) - double precision, intent(inout), dimension(*) :: hp + return - call CopyH(hp) + end subroutine FastScape_Copy_Erosion_Rate - return + !-------------------------------------------------------------------------- -end subroutine FastScape_Copy_H + subroutine FastScape_Copy_Chi (chip,ierr) bind(C, name='fastscape_copy_chi') -!-------------------------------------------------------------------------- + use FastScapeContext -subroutine FastScape_Copy_Basement(bp) + implicit none - use FastScapeContext + real(c_double), intent(inout), dimension(*) :: chip + integer(c_int), optional, intent(out) :: ierr + integer :: ierr_ - implicit none + FSCAPE_INITERR(ierr, ierr_, 'FastScape_Copy_Chi()') - double precision, intent(inout), dimension(*) :: bp + if (.not.setup_has_been_run) then + FSCAPE_RAISE(ERR_SetupNotRun,ierr_);FSCAPE_CHKERR_OPT(ierr, ierr_) + end if - call CopyBasement(bp) + call CopyChi(chip) - return + return -end subroutine FastScape_Copy_Basement + end subroutine FastScape_Copy_Chi -!-------------------------------------------------------------------------- + !-------------------------------------------------------------------------- -subroutine FastScape_Copy_Total_Erosion (etotp) + subroutine FastScape_Copy_Slope (slopep,ierr) bind(C, name='fastscape_copy_slope') - use FastScapeContext + use FastScapeContext - implicit none + implicit none - double precision, intent(inout), dimension(*) :: etotp + real(c_double), intent(inout), dimension(*) :: slopep + integer(c_int), optional, intent(out) :: ierr + integer :: ierr_ - call CopyEtot(etotp) + FSCAPE_INITERR(ierr, ierr_, 'FastScape_Copy_Slope()') - return + if (.not.setup_has_been_run) then + FSCAPE_RAISE(ERR_SetupNotRun,ierr_);FSCAPE_CHKERR_OPT(ierr, ierr_) + end if -end subroutine FastScape_Copy_Total_Erosion + call CopySlope(slopep) -!-------------------------------------------------------------------------- + return -subroutine FastScape_Copy_Drainage_Area (ap) + end subroutine FastScape_Copy_Slope - use FastScapeContext + !-------------------------------------------------------------------------- - implicit none + subroutine FastScape_Copy_Curvature (curvaturep,ierr) bind(C, name='fastscape_copy_curvature') - double precision, intent(inout), dimension(*) :: ap + use FastScapeContext - call CopyArea(ap) + implicit none - return + real(c_double), intent(inout), dimension(*) :: curvaturep + integer(c_int), optional, intent(out) :: ierr + integer :: ierr_ -end subroutine FastScape_Copy_Drainage_Area + FSCAPE_INITERR(ierr, ierr_, 'FastScape_Copy_Curvature()') -!-------------------------------------------------------------------------- + if (.not.setup_has_been_run) then + FSCAPE_RAISE(ERR_SetupNotRun,ierr_);FSCAPE_CHKERR_OPT(ierr, ierr_) + end if -subroutine FastScape_Copy_Erosion_Rate (eratep) + call CopyCurvature(curvaturep) - use FastScapeContext + return - implicit none + end subroutine FastScape_Copy_Curvature - double precision, intent(inout), dimension(*) :: eratep + !-------------------------------------------------------------------------- - call CopyERate(eratep) + subroutine FastScape_Copy_Catchment (catchp,ierr) bind(C, name='fastscape_copy_catchment') - return + use FastScapeContext -end subroutine FastScape_Copy_Erosion_Rate + implicit none -!-------------------------------------------------------------------------- + real(c_double), intent(inout), dimension(*) :: catchp + integer(c_int), optional, intent(out) :: ierr + integer :: ierr_ -subroutine FastScape_Copy_Chi (chip) + FSCAPE_INITERR(ierr, ierr_, 'FastScape_Copy_Catchment()') - use FastScapeContext + if (.not.setup_has_been_run) then + FSCAPE_RAISE(ERR_SetupNotRun,ierr_);FSCAPE_CHKERR_OPT(ierr, ierr_) + end if - implicit none + call CopyCatchment (catchp) - double precision, intent(inout), dimension(*) :: chip + return - call CopyChi(chip) + end subroutine FastScape_Copy_Catchment - return + !-------------------------------------------------------------------------- -end subroutine FastScape_Copy_Chi + subroutine FastScape_Copy_F(Fmixp,ierr) bind(C, name='fastscape_copy_f') -!-------------------------------------------------------------------------- + use FastScapeContext -subroutine FastScape_Copy_Slope (slopep) + implicit none - use FastScapeContext + real(c_double), intent(inout), dimension(*) :: Fmixp + integer(c_int), optional, intent(out) :: ierr + integer :: ierr_ - implicit none + FSCAPE_INITERR(ierr, ierr_, 'FastScape_Copy_F()') - double precision, intent(inout), dimension(*) :: slopep + if (.not.setup_has_been_run) then + FSCAPE_RAISE(ERR_SetupNotRun,ierr_);FSCAPE_CHKERR_OPT(ierr, ierr_) + end if - call CopySlope(slopep) + call CopyF(Fmixp) - return + return -end subroutine FastScape_Copy_Slope + end subroutine FastScape_Copy_F -!-------------------------------------------------------------------------- + !-------------------------------------------------------------------------- -subroutine FastScape_Copy_Curvature (curvaturep) + subroutine FastScape_Copy_Lake_Depth(Lp,ierr) bind(C, name='fastscape_copy_lake_depth') - use FastScapeContext + use FastScapeContext - implicit none + implicit none - double precision, intent(inout), dimension(*) :: curvaturep + real(c_double), intent(inout), dimension(*) :: Lp + integer(c_int), optional, intent(out) :: ierr + integer :: ierr_ - call CopyCurvature(curvaturep) + FSCAPE_INITERR(ierr, ierr_, 'FastScape_Copy_Lake_Depth()') - return + if (.not.setup_has_been_run) then + FSCAPE_RAISE(ERR_SetupNotRun,ierr_);FSCAPE_CHKERR_OPT(ierr, ierr_) + end if -end subroutine FastScape_Copy_Curvature + call CopyLakeDepth(Lp) -!-------------------------------------------------------------------------- + return -subroutine FastScape_Copy_Catchment (catchp) + end subroutine FastScape_Copy_Lake_Depth - use FastScapeContext + !-------------------------------------------------------------------------- - implicit none + subroutine FastScape_Set_NX_NY (nnx,nny,ierr) bind(C, name='fastscape_set_nx_ny') - double precision, intent(inout), dimension(*) :: catchp + use FastScapeContext - call CopyCatchment (catchp) + implicit none - return + integer(c_int), intent(in) :: nnx,nny + integer(c_int), optional, intent(out) :: ierr + integer :: ierr_ -end subroutine FastScape_Copy_Catchment + FSCAPE_INITERR(ierr, ierr_, 'FastScape_Set_NX_NY()') -!-------------------------------------------------------------------------- + call SetNXNY (nnx,nny) -subroutine FastScape_Copy_F(Fmixp) + return - use FastScapeContext + end subroutine FastScape_Set_NX_NY - implicit none + !-------------------------------------------------------------------------- - double precision, intent(inout), dimension(*) :: Fmixp + subroutine FastScape_Set_XL_YL (xxl,yyl,ierr) bind(C, name='fastscape_set_xl_yl') - call CopyF(Fmixp) + use FastScapeContext - return + implicit none -end subroutine FastScape_Copy_F + real(c_double), intent(in) :: xxl,yyl + integer(c_int), optional, intent(out) :: ierr + integer :: ierr_ -!-------------------------------------------------------------------------- + FSCAPE_INITERR(ierr, ierr_, 'FastScape_Set_XL_YL()') -subroutine FastScape_Copy_Lake_Depth(Lp) + call SetXLYL (xxl,yyl) - use FastScapeContext + return - implicit none + end subroutine FastScape_Set_XL_YL - double precision, intent(inout), dimension(*) :: Lp + !-------------------------------------------------------------------------- - call CopyLakeDepth(Lp) + subroutine FastScape_Set_DT (dtt,ierr) bind(C, name='fastscape_set_dt') - return + use FastScapeContext -end subroutine FastScape_Copy_Lake_Depth + implicit none -!-------------------------------------------------------------------------- + real(c_double), intent(in) :: dtt + integer(c_int), optional, intent(out) :: ierr + integer :: ierr_ -subroutine FastScape_Set_NX_NY (nnx,nny) + FSCAPE_INITERR(ierr, ierr_, 'FastScape_Set_DT()') - use FastScapeContext + call SetDT (dtt) - implicit none + return - integer, intent(in) :: nnx,nny + end subroutine FastScape_Set_DT - call SetNXNY (nnx,nny) + !-------------------------------------------------------------------------- - return + subroutine FastScape_Set_Erosional_Parameters (kkf,kkfsed,mm,nnn,kkd,kkdsed,gg1,gg2,pp,ierr) bind(C, name='fastscape_set_erosional_parameters') -end subroutine FastScape_Set_NX_NY + use FastScapeContext -!-------------------------------------------------------------------------- + implicit none -subroutine FastScape_Set_XL_YL (xxl,yyl) + real(c_double), intent(in), dimension(*) :: kkf,kkd + real(c_double), intent(in) :: kkfsed,mm,nnn,kkdsed,gg1,gg2,pp + integer(c_int), optional, intent(out) :: ierr + integer :: ierr_ - use FastScapeContext + FSCAPE_INITERR(ierr, ierr_, 'FastScape_Set_Erosional_Parameters()') - implicit none + call SetErosionalParam (kkf,kkfsed,mm,nnn,kkd,kkdsed,gg1,gg2,pp) - double precision, intent(in) :: xxl,yyl + return - call SetXLYL (xxl,yyl) + end subroutine FastScape_Set_Erosional_Parameters - return + !-------------------------------------------------------------------------- -end subroutine FastScape_Set_XL_YL + subroutine FastScape_Set_Marine_Parameters (sl, p1, p2, z1, z2, r, l, kds1, kds2,ierr) bind(C, name='fastscape_set_marine_parameters') -!-------------------------------------------------------------------------- + use FastScapeContext -subroutine FastScape_Set_DT (dtt) + implicit none - use FastScapeContext + real(c_double), intent(in) :: sl, p1, p2, z1, z2, r, l, kds1, kds2 + integer(c_int), optional, intent(out) :: ierr + integer :: ierr_ - implicit none + FSCAPE_INITERR(ierr, ierr_, 'FastScape_Set_Marine_Parameters()') - double precision, intent(in) :: dtt + call SetMarineParam (sl, p1, p2, z1, z2, r, l, kds1, kds2) - call SetDT (dtt) + return - return + end subroutine FastScape_Set_Marine_Parameters -end subroutine FastScape_Set_DT + !-------------------------------------------------------------------------- -!-------------------------------------------------------------------------- + subroutine FastScape_Get_Sizes (nnx,nny,ierr) bind(C, name='fastscape_get_sizes') -subroutine FastScape_Set_Erosional_Parameters (kkf,kkfsed,mm,nnn,kkd,kkdsed,gg1,gg2,pp) + use FastScapeContext - use FastScapeContext + implicit none - implicit none + integer(c_int), intent(out) :: nnx,nny + integer(c_int), optional, intent(out) :: ierr + integer :: ierr_ - double precision, intent(in), dimension(*) :: kkf,kkd - double precision, intent(in) :: kkfsed,mm,nnn,kkdsed,gg1,gg2,pp + FSCAPE_INITERR(ierr, ierr_, 'FastScape_Get_Sizes()') - call SetErosionalParam (kkf,kkfsed,mm,nnn,kkd,kkdsed,gg1,gg2,pp) + call GetSizes (nnx,nny) - return + return -end subroutine FastScape_Set_Erosional_Parameters + end subroutine FastScape_Get_Sizes -!-------------------------------------------------------------------------- + !-------------------------------------------------------------------------- -subroutine FastScape_Set_Marine_Parameters (sl, p1, p2, z1, z2, r, l, kds1, kds2) + subroutine FastScape_Get_Step (sstep,ierr) bind(C, name='fastscape_get_step') -use FastScapeContext + use FastScapeContext -implicit none + implicit none -double precision, intent(in) :: sl, p1, p2, z1, z2, r, l, kds1, kds2 + integer(c_int), intent(out) :: sstep + integer(c_int), optional, intent(out) :: ierr + integer :: ierr_ -call SetMarineParam (sl, p1, p2, z1, z2, r, l, kds1, kds2) + FSCAPE_INITERR(ierr, ierr_, 'FastScape_Get_Step()') -return + call GetStep (sstep) -end subroutine FastScape_Set_Marine_Parameters + return -!-------------------------------------------------------------------------- + end subroutine FastScape_Get_Step -subroutine FastScape_Get_Sizes (nnx,nny) + !-------------------------------------------------------------------------- - use FastScapeContext + subroutine FastScape_Debug(ierr) bind(C, name='fastscape_debug') - implicit none + use FastScapeContext - integer, intent(out) :: nnx,nny + implicit none - call GetSizes (nnx,nny) + integer(c_int), optional, intent(out) :: ierr + integer :: ierr_ - return + FSCAPE_INITERR(ierr, ierr_, 'FastScape_Debug()') -end subroutine FastScape_Get_Sizes + call Debug() -!-------------------------------------------------------------------------- + return -subroutine FastScape_Get_Step (sstep) + end subroutine FastScape_Debug - use FastScapeContext + !-------------------------------------------------------------------------- - implicit none + subroutine FastScape_Set_BC(jbc,ierr) bind(C, name='fastscape_set_bc') - integer, intent(out) :: sstep + use FastScapeContext - call GetStep (sstep) + implicit none - return + integer(c_int), intent(in) :: jbc + integer(c_int), optional, intent(out) :: ierr + integer :: ierr_ -end subroutine FastScape_Get_Step + FSCAPE_INITERR(ierr, ierr_, 'FastScape_Set_BC()') -!-------------------------------------------------------------------------- + call SetBC (jbc) -subroutine FastScape_Debug() + return - use FastScapeContext + end subroutine FastScape_Set_BC - implicit none + !-------------------------------------------------------------------------- - call Debug() + subroutine FastScape_Set_U (up,ierr) bind(C, name='fastscape_set_u') - return + use FastScapeContext -end subroutine FastScape_Debug + implicit none -!-------------------------------------------------------------------------- + real(c_double), intent(in), dimension(*) :: up + integer(c_int), optional, intent(out) :: ierr + integer :: ierr_ -subroutine FastScape_Set_BC(jbc) + FSCAPE_INITERR(ierr, ierr_, 'FastScape_Set_U()') - use FastScapeContext + call SetU(up) - implicit none + return - integer, intent(in) :: jbc + end subroutine FastScape_Set_U - call SetBC (jbc) + !-------------------------------------------------------------------------- - return + subroutine FastScape_Set_V (ux,uy,ierr) bind(C, name='fastscape_set_v') -end subroutine FastScape_Set_BC + use FastScapeContext -!-------------------------------------------------------------------------- + implicit none -subroutine FastScape_Set_U (up) + real(c_double), intent(in), dimension(*) :: ux,uy + integer(c_int), optional, intent(out) :: ierr + integer :: ierr_ - use FastScapeContext + FSCAPE_INITERR(ierr, ierr_, 'FastScape_Set_V()') - implicit none + call SetV(ux,uy) - double precision, intent(in), dimension(*) :: up + return - call SetU(up) + end subroutine FastScape_Set_V - return + !-------------------------------------------------------------------------- -end subroutine FastScape_Set_U + subroutine FastScape_Reset_Cumulative_Erosion (ierr) bind(C, name='fastscape_reset_cumulative_erosion') -!-------------------------------------------------------------------------- + use FastScapeContext -subroutine FastScape_Set_V (ux,uy) + implicit none - use FastScapeContext + integer(c_int), optional, intent(out) :: ierr + integer :: ierr_ - implicit none + FSCAPE_INITERR(ierr, ierr_, 'FastScape_Reset_Cumulative_Erosion()') - double precision, intent(in), dimension(*) :: ux,uy + call ResetCumulativeErosion () - call SetV(ux,uy) + return - return + end subroutine FastScape_Reset_Cumulative_Erosion -end subroutine FastScape_Set_V + !-------------------------------------------------------------------------- -!-------------------------------------------------------------------------- + subroutine FastScape_Set_H(hp,ierr) bind(C, name='fastscape_set_h') -subroutine FastScape_Reset_Cumulative_Erosion () + use FastScapeContext - use FastScapeContext + implicit none - implicit none + real(c_double), intent(inout), dimension(*) :: hp + integer(c_int), optional, intent(out) :: ierr + integer :: ierr_ - call ResetCumulativeErosion () + FSCAPE_INITERR(ierr, ierr_, 'FastScape_Set_H()') - return + call SetH(hp) -end subroutine FastScape_Reset_Cumulative_Erosion + return -!-------------------------------------------------------------------------- + end subroutine FastScape_Set_H -subroutine FastScape_Set_H(hp) + !-------------------------------------------------------------------------- - use FastScapeContext + subroutine FastScape_Set_All_Layers (dhp,ierr) bind(C, name='fastscape_set_all_layers') - implicit none + use FastScapeContext - double precision, intent(inout), dimension(*) :: hp + implicit none - call SetH(hp) + real(c_double), intent(inout), dimension(*) :: dhp + integer(c_int), optional, intent(out) :: ierr + integer :: ierr_ - return + FSCAPE_INITERR(ierr, ierr_, 'FastScape_Set_All_Layers()') -end subroutine FastScape_Set_H + call SetAllLayers(dhp) -!-------------------------------------------------------------------------- + return -subroutine FastScape_Set_All_Layers (dhp) + end subroutine FastScape_Set_All_Layers - use FastScapeContext + !-------------------------------------------------------------------------- - implicit none + subroutine FastScape_Set_Basement(bp,ierr) bind(C, name='fastscape_set_basement') - double precision, intent(inout), dimension(*) :: dhp + use FastScapeContext - call SetAllLayers(dhp) + implicit none - return + real(c_double), intent(inout), dimension(*) :: bp + integer(c_int), optional, intent(out) :: ierr + integer :: ierr_ -end subroutine FastScape_Set_All_Layers + FSCAPE_INITERR(ierr, ierr_, 'FastScape_Set_Basement()') -!-------------------------------------------------------------------------- + call SetBasement(bp) -subroutine FastScape_Set_Basement(bp) + return - use FastScapeContext + end subroutine FastScape_Set_Basement - implicit none + !-------------------------------------------------------------------------- - double precision, intent(inout), dimension(*) :: bp + subroutine FastScape_Set_Precip (precipp,ierr) bind(C, name='fastscape_set_precip') - call SetBasement(bp) + use FastScapeContext - return + implicit none -end subroutine FastScape_Set_Basement + real(c_double), intent(inout), dimension(*) :: precipp + integer(c_int), optional, intent(out) :: ierr + integer :: ierr_ -!-------------------------------------------------------------------------- + FSCAPE_INITERR(ierr, ierr_, 'FastScape_Set_Precip()') -subroutine FastScape_Set_Precip (precipp) + call SetPrecip (precipp) - use FastScapeContext + return - implicit none + end subroutine FastScape_Set_Precip - double precision, intent(inout), dimension(*) :: precipp + !-------------------------------------------------------------------------- - call SetPrecip (precipp) + subroutine FastScape_VTK (fp, vexp, ierr) bind(C, name='fastscape_vtk') - return + use FastScapeContext -end subroutine FastScape_Set_Precip + implicit none -!-------------------------------------------------------------------------- + real(c_double), intent(in), dimension(*) :: fp + real(c_double), intent(in) :: vexp + integer(c_int), optional, intent(out) :: ierr + integer :: ierr_ -subroutine FastScape_VTK (fp, vexp) + FSCAPE_INITERR(ierr, ierr_, 'FastScape_VTK()') - use FastScapeContext + call Make_VTK (fp, vexp) - implicit none + return - double precision, intent(inout), dimension(*) :: fp - double precision, intent(inout) :: vexp + end subroutine FastScape_VTK - call Make_VTK (fp, vexp) + !-------------------------------------------------------------------------- - return + subroutine FastScape_Strati (nstepp, nreflectorp, nfreqp, vexp, ierr) bind(C, name='fastscape_strati') -end subroutine FastScape_VTK + use FastScapeContext -!-------------------------------------------------------------------------- + implicit none -subroutine FastScape_Strati (nstepp, nreflectorp, nfreqp, vexp) + integer(c_int), intent(inout) :: nstepp, nreflectorp, nfreqp + real(c_double), intent(inout) :: vexp + integer(c_int), optional, intent(out) :: ierr + integer :: ierr_ - use FastScapeContext + FSCAPE_INITERR(ierr, ierr_, 'FastScape_Strati()') - implicit none + call Activate_Strati (nstepp, nreflectorp, nfreqp, vexp) - integer, intent(inout) :: nstepp, nreflectorp, nfreqp - double precision, intent(inout) :: vexp + return - call Activate_Strati (nstepp, nreflectorp, nfreqp, vexp) + end subroutine FastScape_Strati - return + !-------------------------------------------------------------------------- -end subroutine FastScape_Strati + subroutine FastScape_Get_Fluxes (ttectonic_flux, eerosion_flux, bboundary_flux, ierr) bind(C, name='fastscape_get_fluxes') -!-------------------------------------------------------------------------- + use FastScapeContext -subroutine FastScape_Get_Fluxes (ttectonic_flux, eerosion_flux, bboundary_flux) + implicit none - use FastScapeContext + real(c_double), intent(out) :: ttectonic_flux, eerosion_flux, bboundary_flux + integer(c_int), optional, intent(out) :: ierr + integer :: ierr_ - implicit none + FSCAPE_INITERR(ierr, ierr_, 'FastScape_Get_Fluxes()') - double precision, intent(out) :: ttectonic_flux, eerosion_flux, bboundary_flux + call compute_fluxes (ttectonic_flux, eerosion_flux, bboundary_flux) - call compute_fluxes (ttectonic_flux, eerosion_flux, bboundary_flux) + return - return + end subroutine FastScape_Get_Fluxes -end subroutine FastScape_Get_Fluxes +end module FastScapeAPI diff --git a/src/FastScape_ctx.f90 b/src/FastScape_ctx.f90 index 18c63b5..4ee574b 100644 --- a/src/FastScape_ctx.f90 +++ b/src/FastScape_ctx.f90 @@ -1,8 +1,12 @@ + +#include "Error.fpp" + module FastScapeContext ! Context module for FastScape api ! should not be accessed or changed ! see API for name of routines and externally accessible variables + use FastScapeErrorCodes implicit none @@ -35,9 +39,12 @@ module FastScapeContext integer, dimension(:,:), allocatable :: mrec double precision, dimension(:,:), allocatable :: mwrec,mlrec + + contains subroutine Init() + implicit none nx=0 ny=0 @@ -55,16 +62,11 @@ end subroutine Init !--------------------------------------------------------------- subroutine SetUp() - implicit none - if (nx.eq.0) stop 'FastScapeSetup - You need to set nx first' - if (ny.eq.0) stop 'FastScapeSetup - You need to set ny first' - nn=nx*ny call Destroy() - allocate (h(nn),u(nn),vx(nn),vy(nn),stack(nn),ndon(nn),rec(nn),don(8,nn),catch0(nn),catch(nn),precip(nn)) allocate (g(nn)) allocate (bounds_bc(nn)) @@ -159,8 +161,6 @@ subroutine CopyH (hp) double precision, intent(out), dimension(*) :: hp - if (.not.setup_has_been_run) stop 'CopyH - You need to run SetUp first' - hp(1:nn)=h return @@ -173,8 +173,6 @@ subroutine CopyBasement (bp) double precision, intent(out), dimension(*) :: bp - if (.not.setup_has_been_run) stop 'CopyB - You need to run SetUp first' - bp(1:nn)=b return @@ -187,8 +185,6 @@ subroutine CopyEtot (etotp) double precision, intent(inout), dimension(*) :: etotp - if (.not.setup_has_been_run) stop 'CopyEtot - You need to run SetUp first' - etotp(1:nn)=etot return @@ -201,8 +197,6 @@ subroutine CopyArea (ap) double precision, intent(inout), dimension(*) :: ap - if (.not.setup_has_been_run) stop 'CopyArea - You need to run SetUp first' - ap(1:nn)=a return @@ -215,8 +209,6 @@ subroutine CopyErate (eratep) double precision, intent(inout), dimension(*) :: eratep - if (.not.setup_has_been_run) stop 'CopyErate - You need to run SetUp first' - eratep(1:nn)=erate return @@ -232,8 +224,6 @@ subroutine Copychi (chip) integer ij,ijk double precision dx,dy,a0 - if (.not.setup_has_been_run) stop 'CopyChi - You need to run SetUp first' - allocate (chi(nn)) chi=0.d0 dx=xl/(nx-1) @@ -258,8 +248,6 @@ subroutine CopySlope (slopep) double precision, dimension(:), allocatable :: s double precision dx,dy - if (.not.setup_has_been_run) stop 'CopySlope - You need to run SetUp first' - allocate (s(nn)) dx=xl/(nx-1) dy=yl/(ny-1) @@ -279,8 +267,6 @@ subroutine CopyCurvature (curvaturep) double precision, dimension(:), allocatable :: c double precision dx,dy - if (.not.setup_has_been_run) stop 'CopyCurvature - You need to run SetUp first' - allocate (c(nn)) dx=xl/(nx-1) dy=yl/(ny-1) @@ -298,8 +284,6 @@ subroutine CopyCatchment (catchp) double precision, intent(inout), dimension(*) :: catchp - if (.not.setup_has_been_run) stop 'CopyCatchment - You need to run SetUp first' - catchp(1:nn)=catch return @@ -312,8 +296,6 @@ subroutine CopyF (Fmixp) double precision, intent(out), dimension(*) :: Fmixp - if (.not.setup_has_been_run) stop 'CopyF - You need to run SetUp first' - Fmixp(1:nn) = Fmix return @@ -326,8 +308,6 @@ subroutine CopyLakeDepth (Lp) double precision, intent(out), dimension(*) :: Lp - if (.not.setup_has_been_run) stop 'CopyLakeDepth - You need to run SetUp first' - Lp(1:nn) = lake_depth return @@ -340,8 +320,6 @@ subroutine InitH (hp) double precision, intent(in), dimension(*) :: hp - if (.not.setup_has_been_run) stop 'InitH - You need to run SetUp first' - h = hp(1:nn) b = h @@ -355,8 +333,6 @@ subroutine InitF (Fmixp) double precision, intent(in), dimension(*) :: Fmixp - if (.not.setup_has_been_run) stop 'InitF - You need to run SetUp first' - Fmix = Fmixp(1:nn) return diff --git a/src/FlowRouting.f90 b/src/FlowRouting.f90 index 995ae42..79c291a 100644 --- a/src/FlowRouting.f90 +++ b/src/FlowRouting.f90 @@ -1,12 +1,13 @@ +#include "Error.fpp" !-------------------------------------------------------------------------------------------- - -subroutine FlowRouting () +subroutine FlowRouting (ierr) use FastScapeContext implicit none double precision :: dx, dy + integer, intent(inout):: ierr dx = xl/(nx - 1) dy = yl/(ny - 1) @@ -25,11 +26,11 @@ subroutine FlowRouting () call find_stack (rec, don, ndon, nn, catch0, stack, catch) ! removes local minima - call LocalMinima (stack,rec,bounds_bc,ndon,don,h,length,nx,ny,dx,dy) + call LocalMinima (stack,rec,bounds_bc,ndon,don,h,length,nx,ny,dx,dy,ierr);FSCAPE_CHKERR(ierr) ! computes receiver and stack information for mult-direction flow call find_mult_rec (h,rec,stack,hwater,mrec,mnrec,mwrec,mlrec,mstack,nx,ny,dx,dy,p,p_mfd_exp, & - bounds_i1, bounds_i2, bounds_j1, bounds_j2, bounds_xcyclic, bounds_ycyclic) + bounds_i1, bounds_i2, bounds_j1, bounds_j2, bounds_xcyclic, bounds_ycyclic,ierr);FSCAPE_CHKERR(ierr) ! compute lake depth lake_depth = hwater - h @@ -40,7 +41,7 @@ end subroutine FlowRouting !-------------------------------------------------------------------------------------------- -subroutine FlowRoutingSingleFlowDirection () +subroutine FlowRoutingSingleFlowDirection (ierr) use FastScapeContext @@ -48,6 +49,7 @@ subroutine FlowRoutingSingleFlowDirection () integer :: i, ijk, ijr double precision :: dx, dy,deltah + integer, intent(inout):: ierr dx = xl/(nx - 1) dy = yl/(ny - 1) @@ -66,7 +68,7 @@ subroutine FlowRoutingSingleFlowDirection () call find_stack (rec, don, ndon, nn, catch0, stack, catch) ! removes local minima - call LocalMinima (stack,rec,bounds_bc,ndon,don,h,length,nx,ny,dx,dy) + call LocalMinima (stack,rec,bounds_bc,ndon,don,h,length,nx,ny,dx,dy,ierr);FSCAPE_CHKERR(ierr) ! find hwater @@ -145,8 +147,11 @@ end subroutine FlowAccumulationSingleFlowDirection !-------------------------------------------------------------------------------------------- subroutine find_mult_rec (h,rec0,stack0,water,rec,nrec,wrec,lrec,stack,nx,ny,dx,dy,p,p_mfd_exp, & - bounds_i1, bounds_i2, bounds_j1, bounds_j2, bounds_xcyclic, bounds_ycyclic) + bounds_i1, bounds_i2, bounds_j1, bounds_j2, bounds_xcyclic, bounds_ycyclic, ierr) + use FastScapeErrorCodes + + implicit none ! subroutine to find multiple receiver information ! in input: ! h is topography @@ -175,6 +180,7 @@ subroutine find_mult_rec (h,rec0,stack0,water,rec,nrec,wrec,lrec,stack,nx,ny,dx, integer, dimension(:), allocatable :: ndon,vis,parse integer, dimension(:,:), allocatable :: don double precision, dimension(:), allocatable :: h0 + integer, intent(inout):: ierr nn=nx*ny @@ -289,10 +295,13 @@ subroutine find_mult_rec (h,rec0,stack0,water,rec,nrec,wrec,lrec,stack,nx,ny,dx, enddo enddo enddo - if (nstack.ne.nn) stop 'error in stack' deallocate (ndon,don,vis,parse,h0) + if (nstack.ne.nn) then + FSCAPE_RAISE_MESSAGE('Find_mult_rec: error in stack',ERR_Default,ierr);FSCAPE_CHKERR(ierr) + end if + return end subroutine find_mult_rec diff --git a/src/LocalMinima.f90 b/src/LocalMinima.f90 index fce61d4..d588166 100644 --- a/src/LocalMinima.f90 +++ b/src/LocalMinima.f90 @@ -1,4 +1,5 @@ -subroutine LocalMinima (stack,rec,bc,ndon,donor,h,length,nx,ny,dx,dy) +#include "Error.fpp" +subroutine LocalMinima (stack,rec,bc,ndon,donor,h,length,nx,ny,dx,dy,ierr) ! subroutine to compute and remove local inima by recomputing the receiver connectivty ! using Guillaume Cordonnier s algorithm as helped by Benoit Bovy and debuged with Jean Braun @@ -45,6 +46,7 @@ subroutine LocalMinima (stack,rec,bc,ndon,donor,h,length,nx,ny,dx,dy) integer nbasins,basin0,nconn,nconn_max,tree_size integer nn,k,nlocmin + integer, intent(inout):: ierr !continuous_flow = .true. continuous_flow = .false. @@ -91,7 +93,7 @@ subroutine LocalMinima (stack,rec,bc,ndon,donor,h,length,nx,ny,dx,dy) allocate (tree(nbasins-1)) tree_size=0 - call mst_kruskal(conn_weights,conn_basins,nbasins,nconn,tree,tree_size) + call mst_kruskal(conn_weights,conn_basins,nbasins,nconn,tree,tree_size,ierr);FSCAPE_CHKERR(ierr) allocate (sills(nbasins,2)) allocate (basin_stack(nbasins)) @@ -330,7 +332,7 @@ end subroutine connect_basins !---------------------- -subroutine mst_kruskal(conn_weights, conn_basins, nbasins, nconn, mstree, mstree_size) +subroutine mst_kruskal(conn_weights, conn_basins, nbasins, nconn, mstree, mstree_size,ierr) !kruskal algorithm to compute a Minimal Spanning Tree @@ -351,12 +353,13 @@ subroutine mst_kruskal(conn_weights, conn_basins, nbasins, nconn, mstree, mstree integer, dimension(:), allocatable :: sort_id,parent,rank integer mstree_size,eid,eeid,f0,f1,b0,b1 + integer, intent(inout):: ierr allocate (sort_id(nconn)) mstree_size = 0 ! sort edges - call loc_min_3_indexx (nconn,conn_weights,sort_id) + call loc_min_3_indexx (nconn,conn_weights,sort_id,ierr);FSCAPE_CHKERR(ierr) !print*,'weights',conn_weights(sort_id) allocate (parent(nbasins),rank(nbasins)) @@ -944,7 +947,9 @@ end subroutine loc_min_3_find_stack_recurs !----------------------- -subroutine loc_min_3_indexx(n,arr,indx) +subroutine loc_min_3_indexx(n,arr,indx,ierr) + + use FastScapeErrorCodes implicit none @@ -953,6 +958,7 @@ subroutine loc_min_3_indexx(n,arr,indx) parameter (M=7,NSTACK=50) integer i,indxt,ir,itemp,j,jstack,k,l,istack(NSTACK) double precision a + integer, intent(inout):: ierr do j=1,n indx(j)=j @@ -1027,7 +1033,9 @@ subroutine loc_min_3_indexx(n,arr,indx) indx(j)=indxt jstack=jstack+2 - if(jstack.gt.NSTACK) stop 'NSTACK too small in loc_min_3_indexx' + if(jstack.gt.NSTACK) then + FSCAPE_RAISE_MESSAGE('loc_min_3_indexx error: NSTACK too small',ERR_Default,ierr);FSCAPE_CHKERR(ierr) + end if if(ir-i+1.ge.j-l)then istack(jstack)=ir diff --git a/src/Marine.f90 b/src/Marine.f90 index 53c0525..2f1ccc5 100644 --- a/src/Marine.f90 +++ b/src/Marine.f90 @@ -1,4 +1,5 @@ -subroutine Marine() +#include "Error.fpp" +subroutine Marine(ierr) ! Marine transport component ! using silt and sand coupling diffusion solver @@ -15,6 +16,7 @@ subroutine Marine() double precision, dimension(:,:), allocatable :: mmwrec,mmlrec double precision shelfslope,ratio1,ratio2,dx,dy integer ij,ijr,ijk,k + integer, intent(inout):: ierr allocate (flux(nn),shelfdepth(nn),ht(nn),Fs(nn),dh(nn),dh1(nn),dh2(nn),Fmixt(nn),flag(nn)) allocate (dhs(nn),dhs1(nn),F1(nn),F2(nn),zi(nn),zo(nn)) @@ -54,7 +56,7 @@ subroutine Marine() allocate (mmrec(8,nn),mmnrec(nn),mmwrec(8,nn),mmlrec(8,nn),mmstack(nn),mwater(nn)) call find_mult_rec (h,rec,stack,mwater,mmrec,mmnrec,mmwrec,mmlrec,mmstack,nx,ny,dx,dy,0.d0,p_mfd_exp, & - bounds_i1, bounds_i2, bounds_j1, bounds_j2, bounds_xcyclic, bounds_ycyclic) + bounds_i1, bounds_i2, bounds_j1, bounds_j2, bounds_xcyclic, bounds_ycyclic, ierr) !print*,count(flux>0.and.mmnrec==0),count(flux>0),count(mmstack==0) @@ -110,7 +112,7 @@ subroutine Marine() ! silt and sand coupling diffusion in ocean call SiltSandCouplingDiffusion (h,Fmix,flux*Fs,flux*(1.d0-Fs), & - nx,ny,dx,dy,dt,sealevel,layer,kdsea1,kdsea2,nGSMarine,flag,bounds_ibc) + nx,ny,dx,dy,dt,sealevel,layer,kdsea1,kdsea2,nGSMarine,flag,bounds_ibc,ierr);FSCAPE_CHKERR(ierr) ! pure silt and sand during deposition/erosion dh1=((h-ht)*Fmix+layer*(Fmix-Fmixt))*(1.d0-poro1) @@ -169,7 +171,9 @@ end subroutine Marine !---------------------------------------------------------------------------------- subroutine SiltSandCouplingDiffusion (h,f,Q1,Q2,nx,ny,dx,dy,dt, & - sealevel,L,kdsea1,kdsea2,niter,flag,ibc) + sealevel,L,kdsea1,kdsea2,niter,flag,ibc,ierr) + + use FastScapeErrorCodes implicit none @@ -183,6 +187,7 @@ subroutine SiltSandCouplingDiffusion (h,f,Q1,Q2,nx,ny,dx,dy,dt, & double precision dx,dy,dt,sealevel,L,kdsea1,kdsea2 double precision K1,K2,tol,err1,err2 double precision Ap,Bp,Cp,Dp,Ep,Mp,Np + integer, intent(inout):: ierr character cbc*4 @@ -483,8 +488,8 @@ subroutine SiltSandCouplingDiffusion (h,f,Q1,Q2,nx,ny,dx,dy,dt, & !print*,'niter',niter,minval(h-hp),sum(h-hp)/nn,maxval(h-hp),err1 if (niter.gt.1000) then - print*,'Multi-lithology diffusion not convergning; decrease time step' - stop + FSCAPE_RAISE_MESSAGE('Marine error: Multi-lithology diffusion not converging; decrease time step',ERR_NotConverged,ierr) + FSCAPE_CHKERR(ierr) endif ! end of iteration