anchor.diff:                + h(i,j-1,layers)*wetmask(i,j-1)           &
                 + (1d0 - wetmask(i,j-1))*h(i,j,layers)     & ! reflect value around boundary
                 - 2*h(i,j,layers))/(dy*dy) ! y-component horizontal diffusion
  -      enddo
  -    enddo
  +      end do
  +    end do
     else if (.not. RedGrav) then ! using n-layer physics
       ! Calculate bottom layer thickness tendency to balance layers above.
       ! In the flat-bottomed case this will give the same answer.
       dhdt_GM(:,:,layers) = -sum(dhdt_GM(:,:,:layers-1),3)
  -  endif
  +  end if
   
     ! Now add this to the thickness tendency due to the flow field and
     ! sponge regions
  
positive.code:             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
  
negative.code:     vav = 0.0
      if (.not. RedGrav) then
        etaav = 0.0
      end if
    end if
  
    ! For now enforce wetmask to have zeros around the edge.
    wetmask(0,:) = 0d0
  