[Do not merge!] Clean handling of domain boundary conditions

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}/FastScape_types.f90
   ${FASTSCAPELIB_SRC_DIR}/LocalMinima.f90
   ${FASTSCAPELIB_SRC_DIR}/Marine.f90
   ${FASTSCAPELIB_SRC_DIR}/Strati.f90

src/Diffusion.f90

@@ -14,7 +14,7 @@ subroutine Diffusion ()
 
   !print*,'Diffusion'
 
-  write (cbc,'(i4)') ibc
+  write (cbc,'(i4)') bounds%ibc
 
   dx=xl/(nx-1)
   dy=yl/(ny-1)

src/FastScape_ctx.f90

@@ -4,10 +4,13 @@ module FastScapeContext
   ! should not be accessed or changed
   ! see API for name of routines and externally accessible variables
 
+  use FastScapeTypes
+
   implicit none
 
   integer :: nx, ny, nn, nstack
-  integer :: step, ibc
+  type(boundaries), target :: bounds
+  integer :: step
   integer :: nGSStreamPowerLaw, nGSMarine
   logical :: setup_has_been_run
   double precision, target, dimension(:), allocatable :: h,u,vx,vy,length,a,erate,etot,catch,catch0,b,precip,kf,kd
@@ -63,6 +66,7 @@ subroutine SetUp()
 
     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))
     allocate (p_mfd_exp(nn))
     allocate (length(nn),a(nn),erate(nn),etot(nn),b(nn),Sedflux(nn),Fmix(nn),kf(nn),kd(nn))
     allocate (lake_depth(nn),hwater(nn),mrec(8,nn),mnrec(nn),mwrec(8,nn),mlrec(8,nn),mstack(nn))
@@ -142,7 +146,7 @@ subroutine Destroy()
     if (allocated(mstack)) deallocate(mstack)
     if (allocated(g)) deallocate(g)
     if (allocated(p_mfd_exp)) deallocate(p_mfd_exp)
-
+    if (allocated(bounds%bc)) deallocate(bounds%bc)
 
     return
 
@@ -379,7 +383,7 @@ subroutine View()
     write (*,*) 'xl,yl',xl,yl
     write (*,*) 'dt',dt
     write (*,*) 'Kf,Kfsed,,m,n,Kd,Kdsed,G1,G2',sum(kf)/nn,kfsed,m,n,sum(kd)/nn,kdsed,g1,g2
-    write (*,*) 'ibc',ibc
+    write (*,*) 'ibc',bounds%ibc
     write (*,*) 'h',minval(h),sum(h)/nn,maxval(h)
     write (*,*) 'u',minval(u),sum(u)/nn,maxval(u)
 
@@ -509,7 +513,7 @@ subroutine Debug ()
 
     implicit none
 
-    integer i,j,counter,ij,i1,i2,j1,j2
+    integer i,j,ij,counter
     character*4 cbc
 
     write (*,*) '--------------------------------------------------------'
@@ -534,17 +538,8 @@ subroutine Debug ()
     write (*,*) 'Total number of self donors',counter
 
     counter=0
-    write (cbc,'(i4)') ibc
-    i1=1
-    i2=nx
-    j1=1
-    j2=ny
-    if (cbc(4:4).eq.'1') i1=2
-    if (cbc(2:2).eq.'1') i2=nx-1
-    if (cbc(1:1).eq.'1') j1=2
-    if (cbc(3:3).eq.'1') j2=ny-1
-    do j=j1,j2
-      do i=i1,i2
+    do j=bounds%j1,bounds%j2
+      do i=bounds%i1,bounds%i2
         ij=(j-1)*nx+i
         if (rec(ij)==ij) counter=counter+1
       enddo
@@ -569,11 +564,31 @@ end subroutine Debug
 
   subroutine SetBC (jbc)
 
-    integer, intent(in) :: jbc
-
-    ibc = jbc
+    implicit none
 
-    return
+    integer, intent(in) :: jbc
+    character*4 :: cbc
+
+    bounds%ibc = jbc
+
+    write (cbc,'(i4)') jbc
+    bounds%bc=.FALSE.
+    bounds%i1=1
+    bounds%i2=nx
+    bounds%j1=1
+    bounds%j2=ny
+    if (cbc(4:4).eq.'1') bounds%i1=2
+    if (cbc(2:2).eq.'1') bounds%i2=nx-1
+    if (cbc(1:1).eq.'1') bounds%j1=2
+    if (cbc(3:3).eq.'1') bounds%j2=ny-1
+    if (cbc(4:4).eq.'1') bounds%bc(1:nn:nx)=.TRUE.
+    if (cbc(2:2).eq.'1') bounds%bc(nx:nn:nx)=.TRUE.
+    if (cbc(1:1).eq.'1') bounds%bc(1:nx)=.TRUE.
+    if (cbc(3:3).eq.'1') bounds%bc(nx*(ny-1)+1:nn)=.TRUE.
+    bounds%xcyclic=.FALSE.
+    bounds%ycyclic=.FALSE.
+    if (cbc(4:4).ne.'1'.and.cbc(2:2).ne.'1') bounds%xcyclic=.TRUE.
+    if (cbc(1:1).ne.'1'.and.cbc(3:3).ne.'1') bounds%ycyclic=.TRUE.
 
   end subroutine SetBC
 
@@ -831,9 +846,7 @@ subroutine compute_fluxes (tectonic_flux, erosion_flux, boundary_flux)
 
       double precision, intent(out) :: tectonic_flux, erosion_flux, boundary_flux
       double precision :: surf
-      logical, dimension(:), allocatable :: bc
       double precision, dimension(:), allocatable :: flux
-      character*4 :: cbc
       integer ij,ijk,k
 
       surf = xl*yl/(nx - 1)/(ny - 1)
@@ -843,9 +856,9 @@ subroutine compute_fluxes (tectonic_flux, erosion_flux, boundary_flux)
 
       ! computes receiver and stack information for multi-direction flow
       !allocate (mrec(8,nn), mnrec(nn), mwrec(8,nn), mlrec(8,nn), mstack(nn), hwater(nn)
-      !call find_mult_rec (h, rec, stack, hwater, mrec, mnrec, mwrec, mlrec, mstack, nx, ny, xl/(nx-1), yl/(ny-1), p, ibc, p_mfd_exp)
+      !call find_mult_rec (h, rec, stack, hwater, mrec, mnrec, mwrec, mlrec, mstack, nx, ny, xl/(nx-1), yl/(ny-1), p, bounds, p_mfd_exp)
       ! computes sediment flux
-      allocate (flux(nn), bc(nn))
+      allocate (flux(nn))
 
       flux = erate
 
@@ -866,15 +879,9 @@ subroutine compute_fluxes (tectonic_flux, erosion_flux, boundary_flux)
       endif
 
       ! compute boundary flux
-      write (cbc,'(i4)') ibc
-      bc=.FALSE.
-      if (cbc(4:4).eq.'1') bc(1:nn:nx) = .TRUE.
-      if (cbc(2:2).eq.'1') bc(nx:nn:nx) = .TRUE.
-      if (cbc(1:1).eq.'1') bc(1:nx) = .TRUE.
-      if (cbc(3:3).eq.'1') bc(nx*(ny - 1) + 1:nn) = .TRUE.
-      boundary_flux = sum(flux,bc)*surf
-
-      deallocate (flux, bc)
+      boundary_flux = sum(flux,bounds%bc)*surf
+
+      deallocate (flux)
 
     end subroutine compute_fluxes
 

src/FastScape_types.f90

@@ -0,0 +1,14 @@
+module FastScapeTypes
+
+  ! Defines derived types used in FastScape
+
+  implicit none
+
+  type, public :: boundaries
+    integer :: ibc
+    integer :: i1, i2, j1, j2
+    logical :: xcyclic, ycyclic
+    logical, dimension(:), allocatable :: bc
+  end type boundaries
+
+end module FastScapeTypes