anchor.diff:          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
  -
  +    ! Use dh/dt, du/dt and dv/dt to step h, u and v forward in time with
  +    ! the 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.
  +    ! Apply the boundary conditions
  +    ! - Enforce no normal flow boundary condition
  +    !   and no flow in dry cells.
  +    ! - no/free-slip is done inside the 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 step, 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
  +
  +      ! Apply the boundary conditions
  +      ! - Enforce no normal flow boundary condition
  +      !   and no flow in dry cells.
  +      ! - no/free-slip is done inside the 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
  +      enddo
       endif
  -    
  -
  -    
  -! code to stop layers getting too thin
  +
  +    ! Code to stop layers from getting too thin
       counter = 0
  -    
  +
       do k = 1,layers
         do j=1,ny
  -        do i=1,nx   
  +        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
  +            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
  +    ! 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
  +    ! 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
  +      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
  +    ! 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
  +      eta = etanew
       endif
   
  -!    tendencies (old -> very old)
  +    ! Tendencies (old -> very old)
       dhdtveryold = dhdtold
       dudtveryold = dudtold
       dvdtveryold = dvdtold
   
  -!    tendencies (current -> old)
  +    ! 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)     
  -    
  +    ! 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
  -    
  +    ! ----------------------- 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
  +9999 continue
   
       OPEN(UNIT=10, FILE='run_finished.txt', ACTION="write", STATUS="unknown", &
  -    FORM="formatted", POSITION = "append")
  +        FORM="formatted", POSITION = "append")
       write(10,1112) n
  -    1112      format( "run finished at time step ", 1i10.10)  
  +1112 format( "run finished at time step ", 1i10.10)
       CLOSE(UNIT=10)
   
  -    print *, 'Execution ended normally' 
  +    print *, 'Execution ended normally'
   
       stop 0
   
  -    END PROGRAM MIM
  -
  +  END PROGRAM MIM
   
   
   ! ----------------------------- Beginning of the subroutines --------------------
  
positive.code: !> Minimalist Isopycnal Model (MIM) with n layers
  !!
  !
  !>       
  !>     @mainpage Documentation for MIM.f90
  !>
  !>     @section Overview
  
negative.code: !>     layers and arbitrary bathymetry.
  !>
  !>
  !>
  !>    @section Grid
  !>
  !>    /\ ------------
  