C=====================================================
C This subroutine initializes fields and masks arrays
C=====================================================
      SUBROUTINE INI_FIELDS
      implicit none
#include "SIZE.h"
#include "SWE.h"
#include "MASKS.h"

      integer i,j
      integer len_of_rec

      len_of_rec = nx*8

C     First set everything to 0. This will take care of default
C     no normal flow BC
      do j=0,ny+1
         do i=0,nx+1
            u(i,j)=0.d0
            v(i,j)=0.d0
            eta(i,j)=0.d0
            H(i,j)=0.d0
         enddo
      enddo

C     Read initial u, v, and eta from files
      open(unit=22,file=uiniFile,status='old',access='direct',
     &   recl=len_of_rec)
      open(unit=23,file=viniFile,status='old',access='direct',
     &   recl=len_of_rec)
      open(unit=24,file=etainiFile,status='old',access='direct',
     &   recl=len_of_rec)

      write(*,*)'Initial u, v, and eta read from:'
      write(*,*)uiniFile,viniFile,etainiFile

C     Read bathymetry (water depth)
      open(unit=25,file=BathyFile,status='old',access='direct',
     &   recl=len_of_rec)

      write(*,*)'Bathymetry (H) read from'
      write(*,*)BathyFile

      do j=1,ny
         call read_r8seg(nx,u(1,j),22,j)
         call read_r8seg(nx,v(1,j),23,j)
         call read_r8seg(nx,eta(1,j),24,j)
         call read_r8seg(nx,H(1,j),25,j)
      enddo
      close(22)
      close(23)
      close(24)
      close(25)

      if (PeriodicInX) then
         write(*,*) 'Setting periodic in X...'
         do j=0,ny+1
            H(0,j)=H(nx,j)
            H(nx+1,j)=H(1,j)
            eta(0,j)=eta(nx,j)
            u(nx+1,j)=u(1,j)  !  actual periodic BC
            v(0,j)=v(nx,j)  !  for coriolis
         enddo
      endif

      if (PeriodicInY) then
         write(*,*) 'Setting periodic in Y...'
         do i=0,nx+1
            H(i,0)=H(i,ny)
            H(i,ny+1)=H(i,1)
            eta(i,0)=eta(i,ny)
            u(i,0)=u(i,ny)  !  for coriolis
            v(i,ny+1)=v(i,1)  !  actual periodic BC
         enddo
      endif

C     Precompute H on u and v grid: (i,j) below is associated with 
C     eta(i,j). Hup1(i,j) <-> u(i+1,j), Hum1(i,j) <-> u(i,j), 
C     Hvp1(i,j) <-> v(i,j+1), Hvm1(i,j) <-> v(i,j)
      do j=1,ny
         do i=1,nx
            Hup1(i,j) = min(H(i,j),H(i+1,j))
            Hum1(i,j) = min(H(i,j),H(i-1,j))
            Hvp1(i,j) = min(H(i,j),H(i,j+1))
            Hvm1(i,j) = min(H(i,j),H(i,j-1))
         enddo
      enddo

      end