c23456789012345678901234567890123456789012345678901234567890123456789012


      SUBROUTINE flag_vs(opt) 

c     ******************************************************************
c
c     This routine sets a flag for each model variable & each statement 
c
c     ******************************************************************

c     -------------
c     Declarations
c     -------------

c     ''''''''''''''''''''''
      implicit none
      include 'commonvar.f'      
c     ''''''''''''''''''''''

c--   Argument variables

      character	opt*5


c--   Internal variables

      integer	i,j,k




c     --------------
c     Initialization
c     --------------

c---- Flags for bathymetry

      do i=0,itot+1
      do j=0,jtot+1
      do k=0,ktot+1
         flag_b(i,j,k) = .true.
      enddo
      enddo
      enddo


c---- Flags for state variables

      do i=1,itot+1
      do j=1,jtot+1
      do k=1,ktot+1
         flag_r(i,j,k) = .true.
         flag_U(i,j,k) = .true.
         flag_V(i,j,k) = .true.
         flag_W(i,j,k) = .true.
         flag_p(i,j,k) = .true.
      enddo
      enddo
      enddo


c---- Flags for dynamical constraints

      do i=1,itot+1
      do j=1,jtot+1
      do k=1,ktot+1
         flag_mas(i,j,k) = .true.
         flag_twu(i,j,k) = .true.
         flag_twv(i,j,k) = .true.
         flag_vor(i,j,k) = .true.
         flag_den(i,j,k) = .true.
         flag_pos(i,j,k) = .true.
      enddo
      enddo
      enddo




c     -------------------------------------
c     Atlantic domain (5-degree resolution) 
c     -------------------------------------

c---- Test

c     ++++++++++++++++++++++++
      if (opt.eq.'AT_05') then
c     ++++++++++++++++++++++++


c---- Flags for grid boxes 

c     Consider bathymetry

      do i=1,itot
      do j=1,jtot
      do k=1,ktot

         if (z(i,j).le.z_b-k*dz(k)) flag_b(i,j,k) = .false.

      enddo
      enddo
      enddo


c     Remove boxes not removed by bathymetry 
 
      flag_b(1, 5,2) = .false.
      flag_b(2,13,2) = .false.

      flag_b(1, 4,3) = .false.
      flag_b(1, 5,3) = .false.
      flag_b(2,12,3) = .false.
      flag_b(2,13,3) = .false.


c     To facilitate assignation of flag_r(,,) 

      do i=0,itot+1
      do j=0,jtot+1
      do k=0,ktot+1

         if (i.eq.0 .or. i.eq.itot+1 .or.
     >       j.eq.0 .or. j.eq.jtot+1 .or.
     >       k.eq.0 .or. k.eq.ktot+1) flag_b(i,j,k) = .false.

      enddo
      enddo
      enddo


c---- Flags for density 

      do i=1,itot+1
      do j=1,jtot+1
      do k=1,ktot+1

         if (.not. flag_b(i  ,j  ,k  ) .and.
     >	     .not. flag_b(i-1,j  ,k  ) .and.
     >	     .not. flag_b(i-1,j-1,k  ) .and.
     >	     .not. flag_b(i  ,j-1,k  ) .and.
     >	     .not. flag_b(i  ,j  ,k-1) .and.
     >	     .not. flag_b(i-1,j  ,k-1) .and.
     >	     .not. flag_b(i-1,j-1,k-1) .and.
     >	     .not. flag_b(i  ,j-1,k-1)) flag_r(i,j,k) = .false.

      enddo
      enddo
      enddo 


c---- Flags for transport component U 

      do i=1,itot+1
      do j=1,jtot
      do k=1,ktot

         if (.not.flag_b(i-1,j,k) .or.
     >	     .not.flag_b(i  ,j,k)) flag_U(i,j,k) = .false.      

      enddo
      enddo
      enddo


c     Open boundaries 

      flag_U(itot+1, 1,1) = .true.
      flag_U(itot+1, 2,1) = .true.
      flag_U(1     ,15,1) = .true.
      flag_U(1     ,16,1) = .true.
      
      flag_U(itot+1, 1,2) = .true.
      flag_U(itot+1, 2,2) = .true.
      flag_U(itot+1, 3,2) = .true.
      flag_U(1     ,15,2) = .true.
      flag_U(1     ,16,2) = .true.
 
      flag_U(itot+1, 1,3) = .true.
      flag_U(itot+1, 2,3) = .true.
      flag_U(itot+1, 3,3) = .true.
      flag_U(1     ,15,3) = .true.
      flag_U(1     ,16,3) = .true.


c---- Flags for transport component V

      do i=1,itot
      do j=1,jtot+1
      do k=1,ktot

         if (.not.flag_b(i,j-1,k) .or.
     >	     .not.flag_b(i,j  ,k)) flag_V(i,j,k) = .false.    

      enddo
      enddo
      enddo


c     Open boundaries 

      do i=5,10
         flag_V(i,1,1) = .true.
      enddo
      flag_V(  17,1,1) = .true.
      flag_V(itot,1,1) = .true.

      do i=5,11
         flag_V(i,1,2) = .true. 
      enddo
      do i=13,itot
         flag_V(i,1,2) = .true.
      enddo

      do i=5,itot
         flag_V(i,1,3) = .true.
      enddo


c---- Flags for transport component W 

      do i=1,itot
      do j=1,jtot

         flag_W(i,j,1) = .false.
         do k=1+1,ktot+1
            if (.not.flag_b(i,j,k-1)) flag_W(i,j,k) = .false.
         enddo

      enddo
      enddo


c---- Flags for mass conservation

      do i=1,itot
      do j=1,jtot
      do k=1,ktot

         if (.not.flag_b(i,j,k)) flag_mas(i,j,k) = .false.

      enddo
      enddo
      enddo


c---- Flags for thermal wind U
c
c     Note
c     ----
c     1) No thermal wind U at meridional faces within 1 degree of equator

      do i=1,itot+1
      do j=1,jtot
      do k=1,ktot-1

         if (.not.flag_U(i,j,k  ) .or.
     >	     .not.flag_U(i,j,k+1))      flag_twu(i,j,k) = .false.

         if (abs(lat(j  )*deg).le.1.d0) flag_twu(i,j,k) = .false. 
         if (abs(lat(j+1)*deg).le.1.d0) flag_twu(i,j,k) = .false.

      enddo
      enddo
      enddo


c---- Flags for thermal wind V 
c
c     Note
c     ----
c     1) No thermal wind V at zonal faces within 1 degree of equator 

      do i=1,itot
      do j=1,jtot+1
      do k=1,ktot-1

         if (.not.flag_V(i,j,k  ) .or.
     >       .not.flag_V(i,j,k+1))    flag_twv(i,j,k) = .false.

         if (abs(lat(j)*deg).le.1.d0) flag_twv(i,j,k) = .false.


      enddo
      enddo
      enddo


c---- Flags for vorticity 
c
c     Note
c     ----
c     1) No vorticity balance in dry boxes
c     2) No vorticity balance in boxes bounded 
c        by meridional boundary(ies)
c     3) No vorticity balance in boxes touching the equator

      do i=1,itot
      do j=1,jtot
      do k=1,ktot

         if (.not.flag_b(i,j,k))        flag_vor(i,j,k) = .false.

         if (.not.flag_V(i,j  ,k) .or.
     >       .not.flag_V(i,j+1,k))      flag_vor(i,j,k) = .false.

         if (abs(lat(j  )*deg).le.1.d0) flag_vor(i,j,k) = .false.
         if (abs(lat(j+1)*deg).le.1.d0) flag_vor(i,j,k) = .false.

c
         flag_vor(i,j,k) = .false.
c

      enddo
      enddo
      enddo


c---- Flags for density

      do i=1,itot
      do j=1,jtot
      do k=1,ktot

         if (.not.flag_b(i,j,k)) flag_den(i,j,k) = .false.

      enddo
      enddo
      enddo


c---- End of test

c     +++++
      endif
c     +++++




c     -------------------------------------
c     Atlantic domain (3-degree resolution) 
c     -------------------------------------

c---- Test

c     ++++++++++++++++++++++++
      if (opt.eq.'AT_03') then
c     ++++++++++++++++++++++++


c---- Flags for bathymetry  

c     Consider bathymetry

      do i=1,itot
      do j=1,jtot
      do k=1,ktot

         if (z(i,j).le.z_b-k*dz(k)) flag_b(i,j,k) = .false.

      enddo
      enddo
      enddo


c     Remove boxes not removed by bathymetry

      flag_b(1,9,1) = .false.

      do j=5,9
         flag_b(1,j,2) = .false.
      enddo
      flag_b(3,20,2) = .false.
      flag_b(2,21,2) = .false.
      flag_b(3,21,2) = .false.
      flag_b(3,22,2) = .false.  ! To avoid pick-up r(,,) in Carribean Sea (->unrealistic U)

      do j=4,9
         flag_b(1,j,3) = .false.
      enddo
      do i=1,4
         flag_b(i,20,3) = .false.
         flag_b(i,21,3) = .false.
      enddo
      flag_b(3,22,3) = .false.  ! To avoid pick-up r(,,) in Carribean Sea (->unrealistic U)


c     To facilitate assignation of flag_r(,,) 

      do i=0,itot+1
      do j=0,jtot+1
      do k=0,ktot+1

         if (i.eq.0 .or. i.eq.itot+1 .or.
     >       j.eq.0 .or. j.eq.jtot+1 .or.
     >       k.eq.0 .or. k.eq.ktot+1) flag_b(i,j,k) = .false.

      enddo
      enddo
      enddo


c---- Flags for density 

      do i=1,itot+1
      do j=1,jtot+1
      do k=1,ktot+1

         if (.not. flag_b(i  ,j  ,k  ) .and.
     >       .not. flag_b(i-1,j  ,k  ) .and.
     >       .not. flag_b(i-1,j-1,k  ) .and.
     >       .not. flag_b(i  ,j-1,k  ) .and.
     >       .not. flag_b(i  ,j  ,k-1) .and.
     >       .not. flag_b(i-1,j  ,k-1) .and.
     >       .not. flag_b(i-1,j-1,k-1) .and.
     >       .not. flag_b(i  ,j-1,k-1)) flag_r(i,j,k) = .false.

      enddo
      enddo
      enddo


c---- Flags for transport component U 

      do i=1,itot+1
      do j=1,jtot
      do k=1,ktot

         if (.not.flag_b(i-1,j,k) .or.
     >       .not.flag_b(i  ,j,k)) flag_U(i,j,k) = .false.      

      enddo
      enddo
      enddo


c     Open boundaries

      do j=24,26
         flag_U(1,j,1) = .true.
      enddo
      flag_U(itot+1,1,1) = .true.
      flag_U(itot+1,3,1) = .true.
      flag_U(itot+1,4,1) = .true.

      do j=24,27
         flag_U(1,j,2) = .true.
      enddo
      do j=1,4
         flag_U(itot+1,j,2) = .true.
      enddo

      do j=24,27
         flag_U(1,j,3) = .true.
      enddo
      do j=1,5
        flag_U(itot+1,j,3) = .true.
      enddo


c---- Flags for transport component V

      do i=1,itot
      do j=1,jtot+1
      do k=1,ktot

         if (.not.flag_b(i,j-1,k) .or.
     >       .not.flag_b(i,j  ,k)) flag_V(i,j,k) = .false.    

      enddo
      enddo
      enddo


c     Open boundaries

      do i=7,17
         flag_V(i,1,1) = .true.
      enddo
      flag_V(27,1,1) = .true.
      do i=29,itot
         flag_V(i,1,1) = .true.
      enddo

      do i=7,19
         flag_V(i,1,2) = .true.
      enddo
      do i=21,itot
         flag_V(i,1,2) = .true.
      enddo
 
      do i=7,itot
         flag_V(i,1,3) = .true.
      enddo
      flag_V(24,jtot+1,3) = .true.
      flag_V(25,jtot+1,3) = .true.


c---- Flags for transport component W 

      do i=1,itot
      do j=1,jtot

         flag_W(i,j,1) = .false.
         do k=1+1,ktot+1
            if (.not.flag_b(i,j,k-1)) flag_W(i,j,k) = .false.
         enddo

      enddo
      enddo


c---- Flags for mass conservation

      do i=1,itot
      do j=1,jtot
      do k=1,ktot

         if (.not.flag_b(i,j,k)) flag_mas(i,j,k) = .false.

      enddo
      enddo
      enddo


c---- Flags for thermal wind U
c
c     Note
c     ----
c     1) No thermal wind U at meridional faces within 1 degree of equator

      do i=1,itot+1
      do j=1,jtot
      do k=1,ktot-1

         if (.not.flag_U(i,j,k  ) .or.
     >       .not.flag_U(i,j,k+1))      flag_twu(i,j,k) = .false.

         if (abs(lat(j  )*deg).le.1.d0) flag_twu(i,j,k) = .false. 
         if (abs(lat(j+1)*deg).le.1.d0) flag_twu(i,j,k) = .false.

      enddo
      enddo
      enddo


c---- Flags for thermal wind V 
c
c     Note
c     ----
c     1) No thermal wind V at zonal faces within 1 degree of equator 

      do i=1,itot
      do j=1,jtot+1
      do k=1,ktot-1

         if (.not.flag_V(i,j,k  ) .or.
     >       .not.flag_V(i,j,k+1))    flag_twv(i,j,k) = .false.

         if (abs(lat(j)*deg).le.1.d0) flag_twv(i,j,k) = .false.


      enddo
      enddo
      enddo


c---- Flags for vorticity 
c
c     Note
c     ----
c     1) No vorticity balance in dry boxes
c     2) No vorticity balance in boxes bounded 
c        by meridional boundary(ies)
c     3) No vorticity balance in boxes touching the equator

      do i=1,itot
      do j=1,jtot
      do k=1,ktot

         if (.not.flag_b(i,j,k))        flag_vor(i,j,k) = .false.

         if (.not.flag_V(i,j  ,k) .or.
     >       .not.flag_V(i,j+1,k))      flag_vor(i,j,k) = .false.

         if (abs(lat(j  )*deg).le.1.d0) flag_vor(i,j,k) = .false.
         if (abs(lat(j+1)*deg).le.1.d0) flag_vor(i,j,k) = .false.

c
         flag_vor(i,j,k) = .false.
c

      enddo
      enddo
      enddo


c---- Flags for density

      do i=1,itot
      do j=1,jtot
      do k=1,ktot

         if (.not.flag_b(i,j,k)) flag_den(i,j,k) = .false.

      enddo
      enddo
      enddo


c---- End of test

c     +++++
      endif
c     +++++




c     -------------------------------------
c     Atlantic domain (3-degree resolution) 
c     -------------------------------------

c---- Test

c     ++++++++++++++++++++++++
      if (opt.eq.'AT003') then
c     ++++++++++++++++++++++++


c---- Flags for bathymetry  

c     Consider bathymetry

      do i=1,itot
      do j=1,jtot
      do k=1,ktot

         if (z(i,j).le.z_l(k+1)) flag_b(i,j,k) = .false.

      enddo
      enddo
      enddo


c     Remove boxes not removed by bathymetry

      flag_b(1,6,1) = .false.
      flag_b(1,8,1) = .false.
      flag_b(1,9,1) = .false.

      do j=4,9
         flag_b(1,j,2) = .false.
      enddo
      do i=1,3
         flag_b(i,20,2) = .false.
         flag_b(i,21,2) = .false.
      enddo
      do i=1,4
         flag_b(i,22,2) = .false.
      enddo
      flag_b(14,33,2) = .true.	! Charlie-Gibbs
      flag_b(15,33,2) = .true.

      do j=4,9
         flag_b(1,j,3) = .false.
      enddo
      do i=1,5
         flag_b(i,20,3) = .false.
         flag_b(i,21,3) = .false.
      enddo
      do i=1,4
         flag_b(i,22,3) = .false.
      enddo
      flag_b(itot-3,29,3) = .false.
      do i=itot-7,itot
         flag_b(i,jtot,3) = .false.
      enddo


c     Add boxes to allow AABW prescription

      flag_b(13,5,1) = .true.   ! Vema Channel
      flag_b(13,6,1) = .true.
      flag_b(13,7,1) = .true.
    
      flag_b(16,4,1) = .true.   ! Hunter Channel
      flag_b(16,5,1) = .true.

      flag_b(14,15,1) = .true.	! Equatorial passage (36W)
      flag_b(14,16,1) = .true.	
      flag_b(13,16,1) = .true.

      flag_b(12,19,1) = .true.	! Vema Fracture Zone


c     To facilitate assignation of flag_r(,,) 

      do i=0,itot+1
      do j=0,jtot+1
      do k=0,ktot+1

         if (i.eq.0 .or. i.eq.itot+1 .or.
     >       j.eq.0 .or. j.eq.jtot+1 .or.
     >       k.eq.0 .or. k.eq.ktot+1) flag_b(i,j,k) = .false.

      enddo
      enddo
      enddo


c---- Flags for density 

      do i=1,itot+1
      do j=1,jtot+1
      do k=1,ktot+1

         if (.not. flag_b(i  ,j  ,k  ) .and.
     >       .not. flag_b(i-1,j  ,k  ) .and.
     >       .not. flag_b(i-1,j-1,k  ) .and.
     >       .not. flag_b(i  ,j-1,k  ) .and.
     >       .not. flag_b(i  ,j  ,k-1) .and.
     >       .not. flag_b(i-1,j  ,k-1) .and.
     >       .not. flag_b(i-1,j-1,k-1) .and.
     >       .not. flag_b(i  ,j-1,k-1)) flag_r(i,j,k) = .false.

      enddo
      enddo
      enddo


c---- Flags for transport component U 

      do i=1,itot+1
      do j=1,jtot
      do k=1,ktot

         if (.not.flag_b(i-1,j,k) .or.
     >       .not.flag_b(i  ,j,k)) flag_U(i,j,k) = .false.      

      enddo
      enddo
      enddo


c     Open boundaries

      flag_U(itot+1,1,1) = .true.
      flag_U(itot+1,3,1) = .true.
      flag_U(itot+1,4,1) = .true.

      do j=1,5
         flag_U(itot+1,j,2) = .true.
      enddo

      do j=1,5
         flag_U(itot+1,j,3) = .true.
      enddo
      flag_U(23,28    ,3) = .true.	! MOW
      flag_U(21,jtot-1,3) = .true.	! NADW


c---- Flags for transport component V

      do i=1,itot
      do j=1,jtot+1
      do k=1,ktot

         if (.not.flag_b(i,j-1,k) .or.
     >       .not.flag_b(i,j  ,k)) flag_V(i,j,k) = .false.    

      enddo
      enddo
      enddo


c     Open boundaries

      do i=7,17
         flag_V(i,1,1) = .true.
      enddo
      do i=itot-3,itot
         flag_V(i,1,1) = .true.
      enddo

      do i=7,itot
         flag_V(i,1,2) = .true.
      enddo
      flag_V(20,34,2)  = .true.		! NADW

      do i=7,itot
         flag_V(i,1,3) = .true.
      enddo


c---- Flags for transport component W 

      do i=1,itot
      do j=1,jtot

         flag_W(i,j,1) = .false.
         do k=1+1,ktot+1
            if (.not.flag_b(i,j,k-1)) flag_W(i,j,k) = .false.
         enddo

      enddo
      enddo


c---- Flags for tracer 

c++
      if (iflag.ge.1) then
c++


      do i=1,itot+1
      do j=1,jtot+1
      do k=1,ktot+1

         if (.not. flag_b(i  ,j  ,k  ) .and.
     >       .not. flag_b(i-1,j  ,k  ) .and.
     >       .not. flag_b(i-1,j-1,k  ) .and.
     >       .not. flag_b(i  ,j-1,k  ) .and.
     >       .not. flag_b(i  ,j  ,k-1) .and.
     >       .not. flag_b(i-1,j  ,k-1) .and.
     >       .not. flag_b(i-1,j-1,k-1) .and.
     >       .not. flag_b(i  ,j-1,k-1)) flag_p(i,j,k) = .false.

      enddo
      enddo
      enddo

c++
      else 
c++

      do i=1,itot+1
      do j=1,jtot+1
      do k=1,ktot+1
         flag_p(i,j,k) = .false.
      enddo
      enddo
      enddo

c++
      endif
c++


c---- Flags for mass conservation

      do i=1,itot
      do j=1,jtot
      do k=1,ktot

         if (.not.flag_b(i,j,k)) flag_mas(i,j,k) = .false.

      enddo
      enddo
      enddo


c---- Flags for thermal wind U
c
c     Note
c     ----
c     1) No thermal wind U at meridional faces within 1 degree of equator

      do i=1,itot+1
      do j=1,jtot
      do k=1,ktot-1

         if (.not.flag_U(i,j,k  ) .or.
     >       .not.flag_U(i,j,k+1))      flag_twu(i,j,k) = .false.

         if (abs(lat(j  )*deg).le.5.d0) flag_twu(i,j,k) = .false. 
         if (abs(lat(j+1)*deg).le.5.d0) flag_twu(i,j,k) = .false.

      enddo
      enddo
      enddo


c---- Flags for thermal wind V 
c
c     Note
c     ----
c     1) No thermal wind V at zonal faces within 1 degree of equator 

      do i=1,itot
      do j=1,jtot+1
      do k=1,ktot-1

         if (.not.flag_V(i,j,k  ) .or.
     >       .not.flag_V(i,j,k+1))    flag_twv(i,j,k) = .false.

         if (abs(lat(j)*deg).le.5.d0) flag_twv(i,j,k) = .false.


      enddo
      enddo
      enddo


c---- Flags for vorticity equation 
c
c     Note
c     ----
c     1) No vorticity balance in dry boxes
c     2) No vorticity balance in boxes bounded 
c        by meridional boundary(ies)
c     3) No vorticity balance in boxes touching the equator

      do i=1,itot
      do j=1,jtot
      do k=1,ktot

         if (.not.flag_b(i,j,k))        flag_vor(i,j,k) = .false.

         if (.not.flag_V(i,j  ,k) .or.
     >       .not.flag_V(i,j+1,k))      flag_vor(i,j,k) = .false.

         if (abs(lat(j  )*deg).le.5.d0) flag_vor(i,j,k) = .false.
         if (abs(lat(j+1)*deg).le.5.d0) flag_vor(i,j,k) = .false.

c
c         flag_vor(i,j,k) = .false.
c

      enddo
      enddo
      enddo


c---- Flags for density equation

      do i=1,itot
      do j=1,jtot
      do k=1,ktot

         if (.not.flag_b(i,j,k)) flag_den(i,j,k) = .false.

      enddo
      enddo
      enddo


c---- Flags for tracer equation

c++
      if (iflag.ge.1) then
c++

      do i=1,itot
      do j=1,jtot
      do k=1,ktot

         if (.not.flag_b(i,j,k)) flag_pos(i,j,k) = .false.

      enddo
      enddo
      enddo

c++
      else 
c++

      do i=1,itot
      do j=1,jtot
      do k=1,ktot

         flag_pos(i,j,k) = .false.

      enddo
      enddo
      enddo

c++
      endif
c++


c---- End of test

c     +++++
      endif
c     +++++




c     -------------------------------------
c     Atlantic domain (2-degree resolution) 
c     -------------------------------------

c---- Test

c     ++++++++++++++++++++++++
      if (opt.eq.'AT002') then
c     ++++++++++++++++++++++++


c---- Flags for bathymetry  

c     Consider bathymetry

      do i=1,itot
      do j=1,jtot
      do k=1,ktot

         if (z(i,j).le.z_l(k+1)) flag_b(i,j,k) = .false.

      enddo
      enddo
      enddo


c     Remove boxes not removed by bathymetry

c     k = 1

      do i=1,3
      do j=1,14
         flag_b(i,j,1) = .false.
      enddo
      enddo
      
      do i=5,6
      do j=30,31
         flag_b(i,j,1) = .false.
      enddo
      enddo 

c     k = 2

      do i=1,3
      do j=1,15
         flag_b(i,j,2) = .false.
      enddo
      enddo

      do i=1,7
      do j=30,33
         flag_b(i,j,2) = .false.
      enddo
      enddo

      flag_b(24,49,2) = .true.		! Charlie-Gibbs

      flag_b(itot-9,jtot,2) = .false.
      flag_b(itot-8,jtot,2) = .false.

c     k = 3

      do i=1,3
      do j=1,15
         flag_b(i,j,3) = .false.
      enddo
      enddo

      do i=1,8
      do j=29,32
         flag_b(i,j,3) = .false.
      enddo
      enddo

      do i=1,7
         flag_b(i,33,3) = .false.
      enddo

      do i=1,3
         flag_b(i,34,3) = .false.
      enddo

      do i=itot-7,itot
      do j=42,45
         flag_b(i,j,3) = .false.
      enddo
      enddo

      do i=itot-10,itot
      do j=jtot-1 ,jtot
         flag_b(i,j,3) = .false.
      enddo
      enddo


c     Add boxes to allow AABW prescription

      flag_b(20, 8,1) = .true.	! Vema Channel
      flag_b(20, 9,1) = .true.

      flag_b(25, 6,1) = .true.	! Hunter Channel

      flag_b(20,24,1) = .true.	! Equatorial Channel at 36W
      flag_b(21,24,1) = .true.

      flag_b(17,29,1) = .true.	! Vema FZ
      flag_b(18,29,1) = .true.	


c     To facilitate assignation of flag_r(,,) 

      do i=0,itot+1
      do j=0,jtot+1
      do k=0,ktot+1

         if (i.eq.0 .or. i.eq.itot+1 .or.
     >       j.eq.0 .or. j.eq.jtot+1 .or.
     >       k.eq.0 .or. k.eq.ktot+1) flag_b(i,j,k) = .false.

      enddo
      enddo
      enddo


c---- Flags for density 

      do i=1,itot+1
      do j=1,jtot+1
      do k=1,ktot+1

         if (.not. flag_b(i  ,j  ,k  ) .and.
     >       .not. flag_b(i-1,j  ,k  ) .and.
     >       .not. flag_b(i-1,j-1,k  ) .and.
     >       .not. flag_b(i  ,j-1,k  ) .and.
     >       .not. flag_b(i  ,j  ,k-1) .and.
     >       .not. flag_b(i-1,j  ,k-1) .and.
     >       .not. flag_b(i-1,j-1,k-1) .and.
     >       .not. flag_b(i  ,j-1,k-1)) flag_r(i,j,k) = .false.

      enddo
      enddo
      enddo


c---- Flags for transport component U 

      do i=1,itot+1
      do j=1,jtot
      do k=1,ktot

         if (.not.flag_b(i-1,j,k) .or.
     >       .not.flag_b(i  ,j,k)) flag_U(i,j,k) = .false.      

      enddo
      enddo
      enddo


c     Open boundaries

      flag_U(itot+1,1,1) = .true.
      flag_U(itot+1,2,1) = .true.
      flag_U(itot+1,4,1) = .true.
      flag_U(itot+1,5,1) = .true.
      flag_U(itot+1,6,1) = .true.

      do j=1,7
         flag_U(itot+1,j,2) = .true.
      enddo

      do j=1,8
         flag_U(itot+1,j,3) = .true.
      enddo

      flag_U(35,42,3) = .true.	! MOW
      

c---- Flags for transport component V

      do i=1,itot
      do j=1,jtot+1
      do k=1,ktot

         if (.not.flag_b(i,j-1,k) .or.
     >       .not.flag_b(i,j  ,k)) flag_V(i,j,k) = .false.    

      enddo
      enddo
      enddo


c     Open boundaries

      do i=11,27
         flag_V(i,1,1) = .true.
      enddo
      do i=itot-7,itot
         flag_V(i,1,1) = .true.
      enddo

      do i=11,itot
         flag_V(i,1,2) = .true.
      enddo
      flag_V(30,51,2) = .true.		! NADW

      do i=10,itot
         flag_V(i,1,3) = .true.
      enddo

      flag_V(33,jtot-3,3) = .true.	! NADW
      flag_V(31,jtot-1,3) = .true.


c---- Flags for transport component W 

      do i=1,itot
      do j=1,jtot

         flag_W(i,j,1) = .false.
         do k=1+1,ktot+1
            if (.not.flag_b(i,j,k-1)) flag_W(i,j,k) = .false.
         enddo

      enddo
      enddo



c---- Flags for tracer 

c++
      if (iflag.eq.1 .or. iflag.eq.3) then
c++


      do i=1,itot+1
      do j=1,jtot+1
      do k=1,ktot+1

         if (.not. flag_b(i  ,j  ,k  ) .and.
     >       .not. flag_b(i-1,j  ,k  ) .and.
     >       .not. flag_b(i-1,j-1,k  ) .and.
     >       .not. flag_b(i  ,j-1,k  ) .and.
     >       .not. flag_b(i  ,j  ,k-1) .and.
     >       .not. flag_b(i-1,j  ,k-1) .and.
     >       .not. flag_b(i-1,j-1,k-1) .and.
     >       .not. flag_b(i  ,j-1,k-1)) flag_p(i,j,k) = .false.

      enddo
      enddo
      enddo

c++
      else 
c++

      do i=1,itot+1
      do j=1,jtot+1
      do k=1,ktot+1
         flag_p(i,j,k) = .false.
      enddo
      enddo
      enddo

c++
      endif
c++


c---- Flags for mass conservation

      do i=1,itot
      do j=1,jtot
      do k=1,ktot

         if (.not.flag_b(i,j,k)) flag_mas(i,j,k) = .false.

      enddo
      enddo
      enddo


c---- Flags for thermal wind U
c
c     Note
c     ----
c     1) No thermal wind U at meridional faces within 1 degree of equator

      do i=1,itot+1
      do j=1,jtot
      do k=1,ktot-1

         if (.not.flag_U(i,j,k  ) .or.
     >       .not.flag_U(i,j,k+1))      flag_twu(i,j,k) = .false.

         if (abs(lat(j  )*deg).le.5.d0) flag_twu(i,j,k) = .false. 
         if (abs(lat(j+1)*deg).le.5.d0) flag_twu(i,j,k) = .false.

      enddo
      enddo
      enddo


c---- Flags for thermal wind V 
c
c     Note
c     ----
c     1) No thermal wind V at zonal faces within 1 degree of equator 

      do i=1,itot
      do j=1,jtot+1
      do k=1,ktot-1

         if (.not.flag_V(i,j,k  ) .or.
     >       .not.flag_V(i,j,k+1))    flag_twv(i,j,k) = .false.

         if (abs(lat(j)*deg).le.5.d0) flag_twv(i,j,k) = .false.

      enddo
      enddo
      enddo


c---- Flags for vorticity equation 
c
c     Note
c     ----
c     1) No vorticity balance in dry boxes
c     2) No vorticity balance in boxes bounded 
c        by meridional boundary(ies)
c     3) No vorticity balance in boxes touching the equator

      do i=1,itot
      do j=1,jtot
      do k=1,ktot

         if (.not.flag_b(i,j,k))        flag_vor(i,j,k) = .false.

         if (.not.flag_V(i,j  ,k) .or.
     >       .not.flag_V(i,j+1,k))      flag_vor(i,j,k) = .false.

         if (abs(lat(j  )*deg).le.5.d0) flag_vor(i,j,k) = .false.
         if (abs(lat(j+1)*deg).le.5.d0) flag_vor(i,j,k) = .false.

c         flag_vor(i,j,k) = .false.

      enddo
      enddo
      enddo


c---- Flags for density equation

      do i=1,itot
      do j=1,jtot
      do k=1,ktot

         if (.not.flag_b(i,j,k)) flag_den(i,j,k) = .false.

      enddo
      enddo
      enddo


c---- Flags for tracer equation

c++
      if (iflag.eq.1 .or. iflag.eq.3) then
c++

      do i=1,itot
      do j=1,jtot
      do k=1,ktot

         if (.not.flag_b(i,j,k)) flag_pos(i,j,k) = .false.

      enddo
      enddo
      enddo

c++
      else 
c++

      do i=1,itot
      do j=1,jtot
      do k=1,ktot

         flag_pos(i,j,k) = .false.

      enddo
      enddo
      enddo

c++
      endif
c++


c---- End of test

c     +++++
      endif
c     +++++




c     ------------------
c     Flags if iflag = 2 
c     ------------------

      if (iflag.eq.2) then

        do i=1,itot+1
        do j=1,jtot+1
        do k=1,ktot+1
           flag_p  (i,j,k) = .false.
           flag_mas(i,j,k) = .false.
           flag_vor(i,j,k) = .false.
           flag_pos(i,j,k) = .false.
        enddo
        enddo
        enddo

      endif




c     ------------------
c     Flags if iflag = 3 
c     ------------------

      if (iflag.eq.3) then

        do i=1,itot+1
        do j=1,jtot+1
        do k=1,ktot+1
           flag_r  (i,j,k) = .false.
           flag_mas(i,j,k) = .false.
           flag_twu(i,j,k) = .false.
           flag_twv(i,j,k) = .false.
           flag_vor(i,j,k) = .false.
           flag_den(i,j,k) = .false.
        enddo
        enddo
        enddo

      endif




c     ----------------------
c     End of the subroutine
c     ----------------------

      return

      end

c     ******************************************************************
