diff --git a/CMakeLists.txt b/CMakeLists.txt
index ad9c5b2..1b119ee 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -21,6 +21,7 @@ set(FASTSCAPELIB_SRC_FILES
   ${FASTSCAPELIB_SRC_DIR}/Uplift.f90
   ${FASTSCAPELIB_SRC_DIR}/VTK.f90
   ${FASTSCAPELIB_SRC_DIR}/TerrainDerivatives.f90
+  ${FASTSCAPELIB_SRC_DIR}/Named_VTK.f90
   )
 
 set(FLEXURE_SRC_DIR Flexure2D_v1.0/src)
diff --git a/src/Marine.f90 b/src/Marine.f90
index 53c0525..c919a45 100644
--- a/src/Marine.f90
+++ b/src/Marine.f90
@@ -7,20 +7,22 @@ subroutine Marine()
   use FastScapeContext
 
   implicit none
+ 
+  double precision dhs(nn),dhs1(nn)
 
   double precision, dimension(:), allocatable :: flux,shelfdepth,ht,Fs,dh,dh1,dh2,Fmixt,mwater
-  double precision, dimension(:), allocatable :: dhs, dhs1, F1, F2, zi, zo
-  integer, dimension(:), allocatable :: flag,mmnrec,mmstack
+  double precision, dimension(:), allocatable :: F1, F2, zi, zo
+  integer, dimension(:), allocatable :: COTflag,mmnrec,mmstack
   integer, dimension(:,:), allocatable :: mmrec
   double precision, dimension(:,:), allocatable :: mmwrec,mmlrec
   double precision shelfslope,ratio1,ratio2,dx,dy
-  integer ij,ijr,ijk,k
+  integer ij,ijr,ijk,k,istep,i,j
 
-  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))
+  allocate (flux(nn),shelfdepth(nn),ht(nn),Fs(nn),dh(nn),dh1(nn),dh2(nn),Fmixt(nn),COTflag(nn))
+  allocate (F1(nn),F2(nn),zi(nn),zo(nn))
 
   ! set nodes at transition between ocean and continent
-  flag=0
+  COTflag=0
 
   dx=xl/(nx-1)
   dy=yl/(ny-1)
@@ -40,7 +42,14 @@ subroutine Marine()
   where (h.gt.sealevel) flux=0.d0
 
   ! set nodes at transition between ocean and continent
-  !where (flux.gt.tiny(flux)) flag=1
+  !where (flux.gt.tiny(flux)) COTflag=1
+  ! use single flow direction to set the flag that marks
+  ! the recieving node below/at sealevel of a node above/at
+  ! sealevel as a continent-ocean transition node. 
+  do ij=1,nn
+    ijr=rec(ij)
+    if (h(ij).ge.sealevel.and.h(ijr).le.sealevel) COTflag(ijr)=1
+  enddo
 
   ! decompact volume of pure solid phase (silt and sand) from onshore
   ratio1=ratio/(1.d0-poro1)
@@ -110,7 +119,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,COTflag,bounds_ibc)
 
   ! pure silt and sand during deposition/erosion
   dh1=((h-ht)*Fmix+layer*(Fmix-Fmixt))*(1.d0-poro1)
@@ -169,7 +178,7 @@ 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,COTflag,ibc)
 
   implicit none
 
@@ -178,7 +187,7 @@ subroutine SiltSandCouplingDiffusion (h,f,Q1,Q2,nx,ny,dx,dy,dt, &
   double precision h(nx*ny),f(nx*ny),Q1(nx*ny),Q2(nx*ny)
   double precision, dimension(:), allocatable :: hp,fp,ht,ft,hhalf,fhalf,fhalfp
   double precision, dimension(:), allocatable :: diag,sup,inf,rhs,res,tint
-  integer flag(nx*ny)
+  integer COTflag(nx*ny)
 
   double precision dx,dy,dt,sealevel,L,kdsea1,kdsea2
   double precision K1,K2,tol,err1,err2
@@ -224,7 +233,7 @@ subroutine SiltSandCouplingDiffusion (h,f,Q1,Q2,nx,ny,dx,dy,dt, &
         ijp=(j)*nx+i
         ijm=(j-2)*nx+i
         ! in ocean and not at ocean-continent transition
-        if (ht(ij).le.sealevel.and.flag(ij).eq.0) then
+        if (ht(ij).le.sealevel.and.COTflag(ij).eq.0) then
           if (i.eq.1) then
             if (cbc(4:4).eq.'1') then
               diag(i)=1.d0
@@ -298,7 +307,7 @@ subroutine SiltSandCouplingDiffusion (h,f,Q1,Q2,nx,ny,dx,dy,dt, &
         ijp=(j)*nx+i
         ijm=(j-2)*nx+i
         ! in ocean and not at ocean-continent transition
-        if (ht(ij).le.sealevel.and.flag(ij).eq.0) then
+        if (ht(ij).le.sealevel.and.COTflag(ij).eq.0) then
           ! deposition
           if (hhalf(ij).ge.(1.d0+1.d-6)*ht(ij)) then
             Dp=(hhalf(ij)-ht(ij))/dt
@@ -355,7 +364,7 @@ subroutine SiltSandCouplingDiffusion (h,f,Q1,Q2,nx,ny,dx,dy,dt, &
         ijp=(j)*nx+i
         ijm=(j-2)*nx+i
         ! in ocean and not at ocean-continent transition
-        if (ht(ij).le.sealevel.and.flag(ij).eq.0) then
+        if (ht(ij).le.sealevel.and.COTflag(ij).eq.0) then
           if (j.eq.1) then
             if (cbc(1:1).eq.'1') then
               diag(j)=1.d0
@@ -429,7 +438,7 @@ subroutine SiltSandCouplingDiffusion (h,f,Q1,Q2,nx,ny,dx,dy,dt, &
         ijp=(j)*nx+i
         ijm=(j-2)*nx+i
         ! in ocean and not at ocean-continent transition
-        if (ht(ij).le.sealevel.and.flag(ij).eq.0) then
+        if (ht(ij).le.sealevel.and.COTflag(ij).eq.0) then
           ! deposition
           if (h(ij).ge.(1.d0+1.d-6)*hhalf(ij)) then
             Dp=(h(ij)-hhalf(ij))/dt
diff --git a/src/Named_VTK.f90 b/src/Named_VTK.f90
new file mode 100644
index 0000000..9794137
--- /dev/null
+++ b/src/Named_VTK.f90
@@ -0,0 +1,96 @@
+subroutine Fastscape_Named_VTK (f, vex, istep, foldername, k)
+
+    ! Subroutine modified from VTK.f90 to create a simple
+    ! VTK file inside of a specified folder path.
+    ! This additionally outputs the topography, basement
+    ! erosion rate, total erosion, drainage area, and
+    ! catchment on top of the given HHHHH field.
+    ! This assumes the output folder is created beforehand.
+
+    use FastScapeContext
+    implicit none
+
+    integer, intent(in) :: k, istep
+    double precision, intent(in) :: vex
+    double precision, intent(in), dimension(*) :: f
+    character(len=k), intent(in) :: foldername
+    character cstep*7
+
+    integer nheader,nfooter,npart1,npart2, ftime
+    character header*1024,footer*1024,part1*1024,part2*1024,nxc*6,nyc*6,nnc*12
+    integer i,j
+    double precision dx,dy
+
+    dx = xl/(nx - 1)
+    dy = yl/(ny - 1)
+
+    write (cstep,'(i7)') istep
+    if (istep.lt.10) cstep(1:6)='000000'
+    if (istep.lt.100) cstep(1:5)='00000'
+    if (istep.lt.1000) cstep(1:4)='0000'
+    if (istep.lt.10000) cstep(1:3)='000'
+    if (istep.lt.100000) cstep(1:2)='00'
+    if (istep.lt.1000000) cstep(1:1)='0'
+
+    write (nxc,'(i6)') nx
+    write (nyc,'(i6)') ny
+    write (nnc,'(i12)') nn
+
+    header(1:1024)=''
+    header='# vtk DataFile Version 3.0'//char(10)//'FastScape'//char(10) &
+    //'BINARY'//char(10)//'DATASET STRUCTURED_GRID'//char(10) &
+    //'DIMENSIONS '//nxc//' '//nyc//' 1'//char(10)//'POINTS' &
+    //nnc//' float'//char(10)
+    nheader=len_trim(header)
+    footer(1:1024)=''
+    footer='POINT_DATA'//nnc//char(10)
+    nfooter=len_trim(footer)
+    part1(1:1024)=''
+    part1='SCALARS '
+    npart1=len_trim(part1)+1
+    part2(1:1024)=''
+    part2=' float 1'//char(10)//'LOOKUP_TABLE default'//char(10)
+    npart2=len_trim(part2)
+
+    ! Formatting to output the correct size for each visualization field.
+    open(unit=77,file=(foldername//'/Topography'//cstep//'.vtk'),status='unknown', &
+    form='unformatted',access='direct', recl=nheader+3*4*nn+nfooter+(npart1+10+npart2+4*nn) &
+    +(npart1+5+npart2+4*nn)+(npart1+8+npart2+4*nn)+(npart1+12+npart2+4*nn)+(npart1+13+npart2+4*nn) &
+    +(npart1+13+npart2+4*nn)+(npart1+9+npart2+4*nn))
+    write (77,rec=1) &
+    header(1:nheader), &
+    ((sngl(dx*(i-1)),sngl(dy*(j-1)),sngl(h(i+(j-1)*nx)*abs(vex)),i=1,nx),j=1,ny), &
+    footer(1:nfooter), &
+    part1(1:npart1)//'topography'//part2(1:npart2),sngl(h(1:nn)), &
+    part1(1:npart1)//'HHHHH'//part2(1:npart2),sngl(f(1:nn)), &
+    part1(1:npart1)//'basement'//part2(1:npart2),sngl(b(1:nn)), &
+    part1(1:npart1)//'erosion_rate'//part2(1:npart2),sngl(erate(1:nn)), &
+    part1(1:npart1)//'total_erosion'//part2(1:npart2),sngl(etot(1:nn)), &
+    part1(1:npart1)//'drainage_area'//part2(1:npart2),sngl(a(1:nn)) , &
+    part1(1:npart1)//'catchment'//part2(1:npart2),sngl(catch(1:nn))
+    close(77)
+
+    if (vex.lt.0.d0) then
+      open(unit=77,file=(foldername//'/Basement'//cstep//'.vtk'),status='unknown',form='unformatted',access='direct', &
+      recl=nheader+3*4*nn+nfooter+(npart1+1+npart2+4*nn) &
+      +(npart1+5+npart2+4*nn))
+      write (77,rec=1) &
+      header(1:nheader), &
+      ((sngl(dx*(i-1)),sngl(dy*(j-1)),sngl(b(i+(j-1)*nx)*abs(vex)),i=1,nx),j=1,ny), &
+      footer(1:nfooter), &
+      part1(1:npart1)//'B'//part2(1:npart2),sngl(b(1:nn)), &
+      part1(1:npart1)//'HHHHH'//part2(1:npart2),sngl(f(1:nn))
+      close(77)
+      open(unit=77,file=(foldername//'/SeaLevel'//cstep//'.vtk'),status='unknown',form='unformatted',access='direct', &
+      recl=nheader+3*4*nn+nfooter+(npart1+2+npart2+4*nn))
+      write (77,rec=1) &
+      header(1:nheader), &
+      ((sngl(dx*(i-1)),sngl(dy*(j-1)),sngl(sealevel*abs(vex)),i=1,nx),j=1,ny), &
+      footer(1:nfooter), &
+      part1(1:npart1)//'SL'//part2(1:npart2),(sngl(sealevel),i=1,nn)
+      close(77)
+    endif
+
+
+    return
+  end subroutine Fastscape_Named_VTK