anchor.diff:  
       call evaluate_dudt(dudt, h,u,v,b,zeta,wind_x,fu, au,ar,slip,&
           dx,dy,hfacN,hfacS,nx,ny,layers,rho0,&
  -        spongeUTimeScale,spongeU)
  +        spongeUTimeScale,spongeU,RedGrav,botDrag)
   
       call evaluate_dvdt(dvdt, h,u,v,b,zeta,wind_y,fv, au,ar,slip, &
  -        dx,dy,hfacW,hfacE,nx,ny,layers,rho0,spongeVTimeScale,spongeV)
  +        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
  
positive.code: 
      call evaluate_dudt(dudtveryold, h,u,v,b,zeta,wind_x,fu, au,ar,&
          slip,dx,dy,hfacN,hfacS,nx,ny,layers,rho0,&
          spongeUTimeScale,spongeU)
  
      call evaluate_dvdt(dvdtveryold, h,u,v,b,zeta,wind_y,fv, &
          au,ar,slip,dx,dy,hfacW,hfacE,nx,ny,layers,rho0,spongeVTimeScale,spongeV)
  
  !    Calculate the values at half the time interval with Forward Euler
      hhalf = h+0.5d0*dt*dhdtveryold
  
negative.code: 
      call evaluate_dudt(dudtveryold, hhalf,uhalf,vhalf,b,zeta, &
          wind_x,fu, au,ar,slip,dx,dy,hfacN,hfacS,nx,ny,layers,rho0, &
          spongeUTimeScale,spongeU,RedGrav,botDrag)
  
      call evaluate_dvdt(dvdtveryold, hhalf,uhalf,vhalf,b,zeta, &
          wind_y,fv, au,ar,slip ,dx,dy,hfacW,hfacE,nx,ny,layers,rho0,& 
          spongeVTimeScale,spongeV,RedGrav,botDrag)
  !    These are the values to be stored in the 'veryold' variables ready to start 
  !    the proper model run.
  
  