anchor.diff:      call evaluate_zeta(zeta,uhalf,vhalf,nx,ny,layers,dx,dy)
   
   !    Now calculate d/dt of u,v,h and store as dhdtold, dudtold and dvdtold
  -    call evaluate_dhdt(dhdtold, hhalf,uhalf,vhalf,ah,dx,dy,nx,ny,layers, spongeHTimeScale,spongeH,wetmask)
  +    call evaluate_dhdt(dhdtold, hhalf,uhalf,vhalf,ah,dx,dy,nx,ny,layers, spongeHTimeScale,spongeH,wetmask,RedGrav)
       call evaluate_dudt(dudtold, hhalf,uhalf,vhalf,b,zeta,wind_x,&
           fu, au,ar,slip,dx,dy,hfacN,hfacS,nx,ny,layers,rho0,&
           spongeUTimeScale,spongeU,RedGrav,botDrag)
  
positive.code: !------------------------------------------------------------------------------------------
  !> Calculate the tendency of layer thickness for each of the active layers
  !! dh/dt is in the centre of each grid point.
      subroutine evaluate_dhdt(dhdt, h,u,v,ah,dx,dy,nx,ny,layers, spongeTimeScale,spongeH,wetmask)
  !    dhdt is evaluated at the centre of the grid box
      integer nx,ny,layers
      integer i,j,k
      double precision dhdt(0:nx,0:ny,layers), h(0:nx,0:ny,layers)
      double precision u(nx,0:ny,layers),v(0:nx,ny,layers)
      double precision spongeTimeScale(0:nx,0:ny,layers)
      double precision spongeH(0:nx,0:ny,layers)
      double precision wetmask(0:nx,0:ny)
      double precision dx, dy
      double precision ah(layers)
  
      dhdt = 0d0
  
      do k = 1,layers
        do j=1,ny-1
          do i=1,nx-1
            dhdt(i,j,k)=ah(k)*(h(i+1,j,k)*wetmask(i+1,j) + &
                              (1d0 - wetmask(i+1,j))*h(i,j,k) & ! boundary
                          + h(i-1,j,k)*wetmask(i-1,j) + &
                              (1d0 - wetmask(i-1,j))*h(i,j,k) & ! boundary
                          - 2*h(i,j,k))/ &
          (dx*dx) + & ! x-componenet 
            ah(k)*(h(i,j+1,k)*wetmask(i,j+1) + &
                  (1d0 - wetmask(i,j+1))*h(i,j,k) & ! reflect value around boundary
              + h(i,j-1,k)*wetmask(i,j-1) + &
                  (1d0 - wetmask(i,j-1))*h(i,j,k) & ! reflect value around boundary
              - 2*h(i,j,k))/ &
          (dy*dy) - & !y-component horizontal diffusion
          ((h(i,j,k)+h(i+1,j,k))*u(i+1,j,k) - (h(i-1,j,k)+h(i,j,k))*u(i,j,k))/(dx*2d0) - & !d(hu)/dx
          ((h(i,j,k)+h(i,j+1,k))*v(i,j+1,k) - (h(i,j-1,k)+h(i,j,k))*v(i,j,k))/(dy*2d0) + & !d(hv)/dy
          spongeTimeScale(i,j,k)*(spongeH(i,j,k)-h(i,j,k)) ! forced relaxtion in the sponge regions.
  
negative.code:     call evaluate_zeta(zeta,u,v,nx,ny,layers,dx,dy)
  
  !     Calculate dhdt, dudt, dvdt at current time step
      call evaluate_dhdt(dhdtold, h,u,v,ah,dx,dy,nx,ny,layers, spongeHTimeScale,spongeH,wetmask,RedGrav)
  
      call evaluate_dudt(dudtold, h,u,v,b,zeta,wind_x,fu, & 
          au,ar,slip,dx,dy,hfacN,hfacS,nx,ny,layers,rho0, &
  