anchor.diff:  !> Minimalist Isopycnal Model (MIM) with n layers
   !!
   !
  -!>       
  +!>
   !>     @mainpage Documentation for MIM.f90
   !>
   !>     @section Overview
  
positive.code: !>     layers and arbitrary bathymetry.
  !>
  !>
  !>            
  !>    @section Grid
  !>
  !>    /\ ------------
  
negative.code:         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
      ! 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 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
      endif
  
      ! Code to stop layers from getting too thin
      counter = 0
  
      do k = 1,layers
        do j=1,ny
          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
            endif
          end do
        end do
      end do
  
      ! 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
      if (avwrite .ne. 0) then
        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
      h = hnew
      u = unew
      v = vnew
      if (.not. RedGrav) then
        eta = etanew
      endif
  
      ! Tendencies (old -> very old)
      dhdtveryold = dhdtold
      dudtveryold = dudtold
      dvdtveryold = dvdtold
  
      ! 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)
  
      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
  
  120 endif
  
  9999 continue
  
      OPEN(UNIT=10, FILE='run_finished.txt', ACTION="write", STATUS="unknown", &
          FORM="formatted", POSITION = "append")
      write(10,1112) n
  1112 format( "run finished at time step ", 1i10.10)
      CLOSE(UNIT=10)
  
      print *, 'Execution ended normally'
  
      stop 0
  
    END PROGRAM MIM
  
  
  ! ----------------------------- Beginning of the subroutines --------------------
  