anchor.diff:    call read_input_fileV(spongeVfile, spongeV, 0.d0, nx, ny, layers)
     call read_input_fileH_2D(wetMaskFile, wetmask, 1.d0, nx, ny)
   
  +  ! For now enforce wetmask to have zeros around the edge.
  +  wetmask(0, :) = 0d0
  +  wetmask(nx+1, :) = 0d0
  +  wetmask(:, 0) = 0d0
  +  wetmask(:, ny+1) = 0d0
  +  ! TODO, this will change one day when the model can do periodic
  +  ! boundary conditions.
  +
  +  call model_run(h, u, v, eta, depth, dx, dy, wetmask, fu, fv, &
  +      dt, au, ar, botDrag, ah, slip, hmin, nTimeSteps, dumpFreq, avFreq, &
  +      maxits, eps, freesurfFac, thickness_error, g_vec, rho0, &
  +      base_wind_x, base_wind_y, wind_mag_time_series, &
  +      spongeHTimeScale, spongeUTimeScale, spongeVTimeScale, &
  +      spongeH, spongeU, spongeV, &
  +      nx, ny, layers, RedGrav, DumpWind)
  +  print *, 'Execution ended normally'
  +  stop 0
  +end program MIM
  +
  +! ------------------------------ Primary routine ----------------------------
  +!> Run the model
  +
  +subroutine model_run(h, u, v, eta, depth, dx, dy, wetmask, fu, fv, &
  +    dt, au, ar, botDrag, ah, slip, hmin, nTimeSteps, dumpFreq, avFreq, &
  +    maxits, eps, freesurfFac, thickness_error, g_vec, rho0, &
  +    base_wind_x, base_wind_y, wind_mag_time_series, &
  +    spongeHTimeScale, spongeUTimeScale, spongeVTimeScale, &
  +    spongeH, spongeU, spongeV, &
  +    nx, ny, layers, RedGrav, DumpWind)
  +  implicit none
  +
  +  double precision, intent(inout) :: h(0:nx+1, 0:ny+1, layers)
  +  double precision, intent(inout) :: u(0:nx+1, 0:ny+1, layers)
  +  double precision, intent(inout) :: v(0:nx+1, 0:ny+1, layers)
  +  double precision, intent(inout) :: eta(0:nx+1, 0:ny+1)
  +  double precision, intent(in) :: depth(0:nx+1, 0:ny+1)
  +  double precision, intent(in) :: dx, dy
  +  double precision, intent(in) :: wetmask(0:nx+1, 0:ny+1)
  +  double precision, intent(in) :: fu(0:nx+1, 0:ny+1)
  +  double precision, intent(in) :: fv(0:nx+1, 0:ny+1)
  +  double precision, intent(in) :: dt, au, ar, botDrag
  +  double precision, intent(in) :: ah(layers)
  +  double precision, intent(in) :: slip, hmin
  +  integer,          intent(in) :: nTimeSteps
  +  double precision, intent(in) :: dumpFreq, avFreq
  +  integer,          intent(in) :: maxits
  +  double precision, intent(in) :: eps, freesurfFac, thickness_error
  +  double precision, intent(in) :: g_vec(layers)
  +  double precision, intent(in) :: rho0
  +  double precision, intent(in) :: base_wind_x(0:nx+1, 0:ny+1)
  +  double precision, intent(in) :: base_wind_y(0:nx+1, 0:ny+1)
  +  double precision, intent(in) :: wind_mag_time_series(nTimeSteps)
  +  double precision, intent(in) :: spongeHTimeScale(0:nx+1, 0:ny+1, layers)
  +  double precision, intent(in) :: spongeUTimeScale(0:nx+1, 0:ny+1, layers)
  +  double precision, intent(in) :: spongeVTimeScale(0:nx+1, 0:ny+1, layers)
  +  double precision, intent(in) :: spongeH(0:nx+1, 0:ny+1, layers)
  +  double precision, intent(in) :: spongeU(0:nx+1, 0:ny+1, layers)
  +  double precision, intent(in) :: spongeV(0:nx+1, 0:ny+1, layers)
  +  integer,          intent(in) :: nx, ny, layers
  +  logical,          intent(in) :: RedGrav, DumpWind
  +
  +  double precision :: dhdt(0:nx+1, 0:ny+1, layers)
  +  double precision :: dhdtold(0:nx+1, 0:ny+1, layers)
  +  double precision :: dhdtveryold(0:nx+1, 0:ny+1, layers)
  +  double precision :: hnew(0:nx+1, 0:ny+1, layers)
  +  ! for initialisation
  +  double precision :: hhalf(0:nx+1, 0:ny+1, layers)
  +  ! for saving average fields
  +  double precision :: hav(0:nx+1, 0:ny+1, layers)
  +
  +  double precision :: dudt(0:nx+1, 0:ny+1, layers)
  +  double precision :: dudtold(0:nx+1, 0:ny+1, layers)
  +  double precision :: dudtveryold(0:nx+1, 0:ny+1, layers)
  +  double precision :: unew(0:nx+1, 0:ny+1, layers)
  +  ! for initialisation
  +  double precision :: uhalf(0:nx+1, 0:ny+1, layers)
  +  ! for saving average fields
  +  double precision :: uav(0:nx+1, 0:ny+1, layers)
  +
  +  double precision :: dvdt(0:nx+1, 0:ny+1, layers)
  +  double precision :: dvdtold(0:nx+1, 0:ny+1, layers)
  +  double precision :: dvdtveryold(0:nx+1, 0:ny+1, layers)
  +  double precision :: vnew(0:nx+1, 0:ny+1, layers)
  +  ! for initialisation
  +  double precision :: vhalf(0:nx+1, 0:ny+1, layers)
  +  ! for saving average fields
  +  double precision :: vav(0:nx+1, 0:ny+1, layers)
  +
  +  double precision :: etanew(0:nx+1, 0:ny+1)
  +  ! for saving average fields
  +  double precision :: etaav(0:nx+1, 0:ny+1)
  +
  +  ! Pressure solver variables
  +  double precision :: a(5, nx, ny)
  +
  +  ! Geometry
  +  double precision :: hfacW(0:nx+1, 0:ny+1)
  +  double precision :: hfacE(0:nx+1, 0:ny+1)
  +  double precision :: hfacN(0:nx+1, 0:ny+1)
  +  double precision :: hfacS(0:nx+1, 0:ny+1)
  +
  +  ! Numerics
  +  double precision :: pi
  +  integer :: nwrite, avwrite
  +  double precision :: rjac
  +
  +  ! Time step loop variable
  +  integer :: n
  +
  +  ! Wind
  +  double precision :: wind_x(0:nx+1, 0:ny+1)
  +  double precision :: wind_y(0:nx+1, 0:ny+1)
  +
  +  nwrite = int(dumpFreq/dt)
  +  avwrite = int(avFreq/dt)
  +
  +  ! Pi, the constant
  +  pi = 3.1415926535897932384
  +
  +  ! Initialize wind fields
  +  wind_x = base_wind_x*wind_mag_time_series(1)
  +  wind_y = base_wind_y*wind_mag_time_series(1)
  +
     ! Initialise the average fields
     if (avwrite .ne. 0) then
       hav = 0.0
  
positive.code:     end if
    end if
  
    ! For now enforce wetmask to have zeros around the edge.
    wetmask(0, :) = 0d0
    wetmask(nx+1, :) = 0d0
    wetmask(:, 0) = 0d0
    wetmask(:, ny+1) = 0d0
    ! TODO, this will change one day when the model can do periodic
    ! boundary conditions.
  
    call calc_boundary_masks(wetmask, hfacW, hfacE, hfacS, hfacN, nx, ny)
  
    call apply_boundary_conditions(u, hfacW, wetmask, nx, ny, layers)
  
negative.code:   character(30) :: wetMaskFile
  
    ! Numerics
    double precision :: dt
    double precision :: au, ar, botDrag
    double precision :: ah(layerwise_input_length)
    double precision :: slip
    double precision :: hmin
    integer nTimeSteps
    double precision :: dumpFreq, avFreq
    double precision, dimension(:), allocatable :: zeros
    integer maxits
    double precision :: eps
    double precision :: freesurfFac
    double precision :: thickness_error
  
  