anchor.diff:  !>     layers and arbitrary bathymetry.
   !>
   !>
  -!>            
  +!>
   !>    @section Grid
   !>
   !>    /\ ------------
  
positive.code:         dx,dy,hfacW,hfacE,nx,ny,layers,rho0,&
          spongeVTimeScale,spongeV,RedGrav,botDrag)
  
  !    Use dh/dt, du/dt and dv/dt to step h, u and v forward in time with
  !    Adams-Bashforth third order linear multistep method
  
  
      unew=u+dt*(23d0*dudt - 16d0*dudtold + 5d0*dudtveryold)/12d0
  
      vnew=v+dt*(23d0*dvdt - 16d0*dvdtold + 5d0*dvdtveryold)/12d0
  
      hnew=h+dt*(23d0*dhdt - 16d0*dhdtold + 5d0*dhdtveryold)/12d0
  
  !    Apply the boundary conditions
  !    Enforce no normal flow boundary condition
  !    and no flow in dry cells.
  !    no/free-slip is done inside dudt and dvdt subroutines.
  !    hfacW and hfacS are zero where the transition between
  !    wet and dry cells occurs. wetmask is 1 in wet cells,
  !    and zero in dry cells.
      do k = 1,layers
          unew(:,:,k) = unew(:,:,k)*hfacW*wetmask(:,:)
          vnew(:,:,k) = vnew(:,:,k)*hfacS*wetmask(:,:)
      enddo
  
  
  !   Do the isopycnal layer physics
      if (.not. RedGrav) then
          ! calculate the barotropic velocities
          call calc_baro_u(ub,unew,hnew,eta,freesurfFac,nx,ny,layers)
          call calc_baro_v(vb,vnew,hnew,eta,freesurfFac,nx,ny,layers)
  
          ! calculate divergence of ub and vb, and solve for the pressure field that removes it
          call calc_eta_star(ub,vb,eta,etastar,freesurfFac,nx,ny,dx,dy,dt)
          !print *, maxval(abs(etastar))
  
          ! prevent barotropic signals from bouncing around outside the wet region of the model.
          !etastar = etastar*wetmask
  
          call SOR_solver(a,etanew,etastar,freesurfFac,nx,ny, &
              dt,rjac,eps,maxits,n)
          !print *, maxval(abs(etanew))
  
          etanew = etanew*wetmask
          call wrap_fields_2D(etanew,nx,ny)
  
          ! now update the velocities using the barotropic tendency due to the pressure gradient
          dudt_bt = 0d0
          dvdt_bt = 0d0
  
          do i = 1,nx
              do j = 0,ny
                  do k = 1,layers
                      dudt_bt(i,j) = -g_vec(1)*(etanew(i,j) - etanew(i-1,j))/(dx)
                      unew(i,j,k) = unew(i,j,k) + dt*dudt_bt(i,j) !23d0*dt*dudt_bt(i,j)/12d0
                  end do
              end do
          end do
  
          do i = 0,nx
              do j = 1,ny
                  do k = 1,layers
                      dvdt_bt(i,j) = -g_vec(1)*(etanew(i,j) - etanew(i,j-1))/(dy)
                      vnew(i,j,k) = vnew(i,j,k) + dt*dvdt_bt(i,j)! 23d0*dt*dvdt_bt(i,j)/12d0
                  end do
              end do
          end do
  
  
          ! We now have correct velocities at the next time level, but the layer thicknesses were updated using the velocities without the barotropic pressure contribution. force consistency between layer thicknesses and ocean depth by scaling thicknesses to agree with free surface.
  
          h_norming = (freesurfFac*etanew + depth)/sum(hnew,3)
          do k = 1,layers
              hnew(:,:,k) = hnew(:,:,k)*h_norming
          end do
  
          if (maxval(abs(h_norming - 1d0)).gt. 1e-2) then
              print *, 'substantial inconsistency between h and eta (in %).', maxval(abs(h_norming - 1d0))*100
          end if
  
      !    Enforce no normal flow boundary condition
      !    and no flow in dry cells.
      !    no/free-slip is done inside dudt and dvdt subroutines.
      !    hfacW and hfacS are zero where the transition between
      !    wet and dry cells occurs. wetmask is 1 in wet cells,
      !    and zero in dry cells.
          do k = 1,layers
              unew(:,:,k) = unew(:,:,k)*hfacW*wetmask(:,:)
              vnew(:,:,k) = vnew(:,:,k)*hfacS*wetmask(:,:)
          enddo
      endif
      
  
      
  ! code to stop layers getting too thin
      counter = 0
      
      do k = 1,layers
        do j=1,ny
          do i=1,nx   
            if (hnew(i,j,k) .lt. hmin) then
          hnew(i,j,k) = hmin
          counter = counter + 1
          if (counter .eq. 1) then
            ! write a file saying that the layer thickness value dropped below hmin and this line has been used
            OPEN(UNIT=10, FILE='layer thickness dropped below hmin.txt', ACTION="write", STATUS="unknown", &
            FORM="formatted", POSITION = "append")
            write(10,1111) n
  1111              format( "layer thickness dropped below hmin at time step ", 1i10.10)  
            CLOSE(UNIT=10)
          endif
            endif
          end do
        end do
      end do
  
  ! wrap fields around for periodic simulations
      call wrap_fields_3D(unew,nx,ny,layers)
      call wrap_fields_3D(vnew,nx,ny,layers)
      call wrap_fields_3D(hnew,nx,ny,layers)
      call wrap_fields_2D(etanew,nx,ny)
  
  
  !    Accumulate average fields
      if (avwrite .ne. 0) then
          hav = hav + hnew
          uav = uav + unew
          vav = vav + vnew
          if (.not. RedGrav) then
              etaav = eta + etanew
          endif
      end if
  
  !     shuffle arrays: old -> very old,  present -> old, new -> present
  !    height and velocity fields
      h = hnew
      u = unew
      v = vnew
      if (.not. RedGrav) then
          eta = etanew
      endif
  
  !    tendencies (old -> very old)
      dhdtveryold = dhdtold
      dudtveryold = dudtold
      dvdtveryold = dvdtold
  
  !    tendencies (current -> old)
      dudtold = dudt
      dvdtold = dvdt
      dhdtold = dhdt
  
  
  !     now have new fields in main arrays and old fields in very old arrays
  
  !--------------------------- Code for generating output -------------------------------------
  !     write data to file? 
        if (mod(n-1,nwrite).eq.0) then 
  
          !     
          !----------------------- Snapshot Output -----------------
          write(num,'(i10.10)') n
  
          !output the data to a file
  
          open(unit = 10, status='replace',file='output/snap.h.'//num,form='unformatted')  
          write(10) h(1:nx,1:ny,:)
          close(10) 
          open(unit = 10, status='replace',file='output/snap.u.'//num,form='unformatted')  
          write(10) u(1:nx+1,1:ny,:)
          close(10) 
          open(unit = 10, status='replace',file='output/snap.v.'//num,form='unformatted')  
          write(10) v(1:nx,1:ny+1,:)
          close(10) 
          if (.not. RedGrav) then
              open(unit = 10, status='replace',file='output/snap.eta.'//num,form='unformatted')  
              write(10) eta(1:nx,1:ny)
              close(10) 
          endif
  
          if (DumpWind .eqv. .TRUE.) then
              open(unit = 10, status='replace',file='output/wind_x.'//num,form='unformatted')  
              write(10) wind_x(1:nx+1,1:ny)
              close(10) 
              open(unit = 10, status='replace',file='output/wind_y.'//num,form='unformatted')  
              write(10) wind_y(1:nx,1:ny+1)
              close(10)
          endif 
  
          ! Check if there are NaNs in the data
          call break_if_NaN(h,nx,ny,layers,n)
          !     call break_if_NaN(u,nx,ny,layers,n)
          !     call break_if_NaN(v,nx,ny,layers,n)     
      
      endif
  
          !----------------------- Average Output -----------------
        if (avwrite .eq. 0) then 
          go to 120
        else if (mod(n-1,avwrite).eq.0) then 
  
          hav=hav/real(avwrite)
          uav=uav/real(avwrite)
          vav=vav/real(avwrite)
  
          write(num,'(i10.10)') n
  
          !output the data to a file
  
          open(unit = 10, status='replace',file='output/av.h.'//num,form='unformatted')  
          write(10) hav(1:nx,1:ny,:)
          close(10) 
          open(unit = 10, status='replace',file='output/av.u.'//num,form='unformatted')  
          write(10) uav(1:nx+1,1:ny,:)
          close(10) 
          open(unit = 10, status='replace',file='output/av.v.'//num,form='unformatted')  
          write(10) vav(1:nx,1:ny+1,:)
          close(10) 
          if (.not. RedGrav) then
              open(unit = 10, status='replace',file='output/av.eta.'//num,form='unformatted')  
              write(10) etaav(1:nx,1:ny)
              close(10) 
          endif
  
          ! Check if there are NaNs in the data
          call break_if_NaN(h,nx,ny,layers,n)
          !     call break_if_NaN(u,nx,ny,layers,n)
          !     call break_if_NaN(v,nx,ny,layers,n)     
  
          !    Reset average quantities
          hav = 0.0
          uav = 0.0
          vav = 0.0
          if (.not. RedGrav) then
              etaav = 0.0
          endif
          !       h2av=0.0
      
  120 endif
  
  9999  continue
  
      OPEN(UNIT=10, FILE='run_finished.txt', ACTION="write", STATUS="unknown", &
      FORM="formatted", POSITION = "append")
      write(10,1112) n
      1112      format( "run finished at time step ", 1i10.10)  
      CLOSE(UNIT=10)
  
      print *, 'Execution ended normally' 
  
      stop 0
  
      END PROGRAM MIM
  
  
  
  ! ----------------------------- Beginning of the subroutines --------------------
  
negative.code: !> Minimalist Isopycnal Model (MIM) with n layers
  !!
  !
  !>
  !>     @mainpage Documentation for MIM.f90
  !>
  !>     @section Overview
  