ecearth3/sources/nemo-3.6/NEMO/OPA_SRC/DOM/dommsk.F90

MODULE dommsk
   !!======================================================================
   !!                       ***  MODULE dommsk   ***
   !! Ocean initialization : domain land/sea mask 
   !!======================================================================
   !! History :  OPA  ! 1987-07  (G. Madec)  Original code
   !!            6.0  ! 1993-03  (M. Guyon)  symetrical conditions (M. Guyon)
   !!            7.0  ! 1996-01  (G. Madec)  suppression of common work arrays
   !!             -   ! 1996-05  (G. Madec)  mask computed from tmask and sup-
   !!                 !                      pression of the double computation of bmask
   !!            8.0  ! 1997-02  (G. Madec)  mesh information put in domhgr.F
   !!            8.1  ! 1997-07  (G. Madec)  modification of mbathy and fmask
   !!             -   ! 1998-05  (G. Roullet)  free surface
   !!            8.2  ! 2000-03  (G. Madec)  no slip accurate
   !!             -   ! 2001-09  (J.-M. Molines)  Open boundaries
   !!   NEMO     1.0  ! 2002-08  (G. Madec)  F90: Free form and module
   !!             -   ! 2005-11  (V. Garnier) Surface pressure gradient organization
   !!            3.2  ! 2009-07  (R. Benshila) Suppression of rigid-lid option
   !!----------------------------------------------------------------------

   !!----------------------------------------------------------------------
   !!   dom_msk        : compute land/ocean mask
   !!   dom_msk_nsa    : update land/ocean mask when no-slip accurate option is used.
   !!----------------------------------------------------------------------
   USE oce             ! ocean dynamics and tracers
   USE dom_oce         ! ocean space and time domain
   USE in_out_manager  ! I/O manager
   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
   USE lib_mpp
   USE dynspg_oce      ! choice/control of key cpp for surface pressure gradient
   USE wrk_nemo        ! Memory allocation
   USE timing          ! Timing

   IMPLICIT NONE
   PRIVATE

   PUBLIC   dom_msk         ! routine called by inidom.F90
   PUBLIC   dom_msk_alloc   ! routine called by nemogcm.F90

   !                            !!* Namelist namlbc : lateral boundary condition *
   REAL(wp)        :: rn_shlat   ! type of lateral boundary condition on velocity
   LOGICAL, PUBLIC :: ln_vorlat  !  consistency of vorticity boundary condition 
   !                                            with analytical eqs.


   INTEGER, ALLOCATABLE, SAVE, DIMENSION(:,:) ::  icoord ! Workspace for dom_msk_nsa()

   !! * Substitutions
#  include "vectopt_loop_substitute.h90"
   !!----------------------------------------------------------------------
   !! NEMO/OPA 3.2 , LODYC-IPSL  (2009)
   !! $Id: dommsk.F90 5551 2015-07-06 16:38:51Z acc $ 
   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
   !!----------------------------------------------------------------------
CONTAINS
   
   INTEGER FUNCTION dom_msk_alloc()
      !!---------------------------------------------------------------------
      !!                 ***  FUNCTION dom_msk_alloc  ***
      !!---------------------------------------------------------------------
      dom_msk_alloc = 0
#if defined key_noslip_accurate
      ALLOCATE(icoord(jpi*jpj*jpk,3), STAT=dom_msk_alloc)
#endif
      IF( dom_msk_alloc /= 0 )   CALL ctl_warn('dom_msk_alloc: failed to allocate icoord array')
      !
   END FUNCTION dom_msk_alloc


   SUBROUTINE dom_msk
      !!---------------------------------------------------------------------
      !!                 ***  ROUTINE dom_msk  ***
      !!
      !! ** Purpose :   Compute land/ocean mask arrays at tracer points, hori-
      !!      zontal velocity points (u & v), vorticity points (f) and baro-
      !!      tropic stream function  points (b).
      !!
      !! ** Method  :   The ocean/land mask is computed from the basin bathy-
      !!      metry in level (mbathy) which is defined or read in dommba.
      !!      mbathy equals 0 over continental T-point 
      !!      and the number of ocean level over the ocean.
      !!
      !!      At a given position (ji,jj,jk) the ocean/land mask is given by:
      !!      t-point : 0. IF mbathy( ji ,jj) =< 0
      !!                1. IF mbathy( ji ,jj) >= jk
      !!      u-point : 0. IF mbathy( ji ,jj)  or mbathy(ji+1, jj ) =< 0
      !!                1. IF mbathy( ji ,jj) and mbathy(ji+1, jj ) >= jk.
      !!      v-point : 0. IF mbathy( ji ,jj)  or mbathy( ji ,jj+1) =< 0
      !!                1. IF mbathy( ji ,jj) and mbathy( ji ,jj+1) >= jk.
      !!      f-point : 0. IF mbathy( ji ,jj)  or mbathy( ji ,jj+1)
      !!                   or mbathy(ji+1,jj)  or mbathy(ji+1,jj+1) =< 0
      !!                1. IF mbathy( ji ,jj) and mbathy( ji ,jj+1)
      !!                  and mbathy(ji+1,jj) and mbathy(ji+1,jj+1) >= jk.
      !!      b-point : the same definition as for f-point of the first ocean
      !!                level (surface level) but with 0 along coastlines.
      !!      tmask_i : interior ocean mask at t-point, i.e. excluding duplicated
      !!                rows/lines due to cyclic or North Fold boundaries as well
      !!                as MPP halos.
      !!
      !!        The lateral friction is set through the value of fmask along
      !!      the coast and topography. This value is defined by rn_shlat, a
      !!      namelist parameter:
      !!         rn_shlat = 0, free slip  (no shear along the coast)
      !!         rn_shlat = 2, no slip  (specified zero velocity at the coast)
      !!         0 < rn_shlat < 2, partial slip   | non-linear velocity profile
      !!         2 < rn_shlat, strong slip        | in the lateral boundary layer
      !!
      !!      N.B. If nperio not equal to 0, the land/ocean mask arrays
      !!      are defined with the proper value at lateral domain boundaries,
      !!      but bmask. indeed, bmask defined the domain over which the
      !!      barotropic stream function is computed. this domain cannot
      !!      contain identical columns because the matrix associated with
      !!      the barotropic stream function equation is then no more inverti-
      !!      ble. therefore bmask is set to 0 along lateral domain boundaries
      !!      even IF nperio is not zero.
      !!
      !!      In case of open boundaries (lk_bdy=T):
      !!        - tmask is set to 1 on the points to be computed bay the open
      !!          boundaries routines.
      !!        - bmask is  set to 0 on the open boundaries.
      !!
      !! ** Action :   tmask    : land/ocean mask at t-point (=0. or 1.)
      !!               umask    : land/ocean mask at u-point (=0. or 1.)
      !!               vmask    : land/ocean mask at v-point (=0. or 1.)
      !!               fmask    : land/ocean mask at f-point (=0. or 1.)
      !!                          =rn_shlat along lateral boundaries
      !!               bmask    : land/ocean mask at barotropic stream
      !!                          function point (=0. or 1.) and set to 0 along lateral boundaries
      !!               tmask_i  : interior ocean mask
      !!----------------------------------------------------------------------
      !
      INTEGER  ::   ji, jj, jk      ! dummy loop indices
      INTEGER  ::   iif, iil, ii0, ii1, ii   ! local integers
      INTEGER  ::   ijf, ijl, ij0, ij1       !   -       -
      INTEGER  ::   ios
      INTEGER  ::   isrow                    ! index for ORCA1 starting row
      INTEGER , POINTER, DIMENSION(:,:) ::  imsk
      REAL(wp), POINTER, DIMENSION(:,:) ::  zwf
      !!
      NAMELIST/namlbc/ rn_shlat, ln_vorlat
      !!---------------------------------------------------------------------
      !
      IF( nn_timing == 1 )  CALL timing_start('dom_msk')
      !
      CALL wrk_alloc( jpi, jpj, imsk )
      CALL wrk_alloc( jpi, jpj, zwf  )
      !
      REWIND( numnam_ref )              ! Namelist namlbc in reference namelist : Lateral momentum boundary condition
      READ  ( numnam_ref, namlbc, IOSTAT = ios, ERR = 901 )
901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namlbc in reference namelist', lwp )

      REWIND( numnam_cfg )              ! Namelist namlbc in configuration namelist : Lateral momentum boundary condition
      READ  ( numnam_cfg, namlbc, IOSTAT = ios, ERR = 902 )
902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namlbc in configuration namelist', lwp )
      IF(lwm) WRITE ( numond, namlbc )
      
      IF(lwp) THEN                  ! control print
         WRITE(numout,*)
         WRITE(numout,*) 'dommsk : ocean mask '
         WRITE(numout,*) '~~~~~~'
         WRITE(numout,*) '   Namelist namlbc'
         WRITE(numout,*) '      lateral momentum boundary cond.    rn_shlat  = ',rn_shlat
         WRITE(numout,*) '      consistency with analytical form   ln_vorlat = ',ln_vorlat 
      ENDIF

      IF     (      rn_shlat == 0.               ) THEN   ;   IF(lwp) WRITE(numout,*) '   ocean lateral  free-slip '
      ELSEIF (      rn_shlat == 2.               ) THEN   ;   IF(lwp) WRITE(numout,*) '   ocean lateral  no-slip '
      ELSEIF ( 0. < rn_shlat .AND. rn_shlat < 2. ) THEN   ;   IF(lwp) WRITE(numout,*) '   ocean lateral  partial-slip '
      ELSEIF ( 2. < rn_shlat                     ) THEN   ;   IF(lwp) WRITE(numout,*) '   ocean lateral  strong-slip '
      ELSE
         WRITE(ctmp1,*) ' rn_shlat is negative = ', rn_shlat
         CALL ctl_stop( ctmp1 )
      ENDIF

      ! 1. Ocean/land mask at t-point (computed from mbathy)
      ! -----------------------------
      ! N.B. tmask has already the right boundary conditions since mbathy is ok
      !
      tmask(:,:,:) = 0._wp
      DO jk = 1, jpk
         DO jj = 1, jpj
            DO ji = 1, jpi
               IF( REAL( mbathy(ji,jj) - jk, wp ) + 0.1_wp >= 0._wp )   tmask(ji,jj,jk) = 1._wp
            END DO  
         END DO  
      END DO  
      
      ! (ISF) define barotropic mask and mask the ice shelf point
      ssmask(:,:)=tmask(:,:,1) ! at this stage ice shelf is not masked
      
      DO jk = 1, jpk
         DO jj = 1, jpj
            DO ji = 1, jpi
               IF( REAL( misfdep(ji,jj) - jk, wp ) - 0.1_wp >= 0._wp )   THEN
                  tmask(ji,jj,jk) = 0._wp
               END IF
            END DO  
         END DO  
      END DO  

!!gm  ????
#if defined key_zdfkpp
      IF( cp_cfg == 'orca' ) THEN
         IF( jp_cfg == 2 )   THEN       ! land point on Bab el Mandeb zonal section
            ij0 =  87   ;   ij1 =  88
            ii0 = 160   ;   ii1 = 161
            tmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 0._wp
         ELSE
            IF(lwp) WRITE(numout,*)
            IF(lwp) WRITE(numout,cform_war)
            IF(lwp) WRITE(numout,*)
            IF(lwp) WRITE(numout,*)'          A mask must be applied on Bab el Mandeb strait'
            IF(lwp) WRITE(numout,*)'          in case of ORCAs configurations'
            IF(lwp) WRITE(numout,*)'          This is a problem which is not yet solved'
            IF(lwp) WRITE(numout,*)
         ENDIF
      ENDIF
#endif
!!gm end

      ! Interior domain mask (used for global sum)
      ! --------------------
      tmask_i(:,:) = ssmask(:,:)            ! (ISH) tmask_i = 1 even on the ice shelf
      iif = jpreci                         ! ???
      iil = nlci - jpreci + 1
      ijf = jprecj                         ! ???
      ijl = nlcj - jprecj + 1

      tmask_i( 1 :iif,   :   ) = 0._wp      ! first columns
      tmask_i(iil:jpi,   :   ) = 0._wp      ! last  columns (including mpp extra columns)
      tmask_i(   :   , 1 :ijf) = 0._wp      ! first rows
      tmask_i(   :   ,ijl:jpj) = 0._wp      ! last  rows (including mpp extra rows)

      ! north fold mask
      ! ---------------
      tpol(1:jpiglo) = 1._wp 
      fpol(1:jpiglo) = 1._wp
      IF( jperio == 3 .OR. jperio == 4 ) THEN      ! T-point pivot
         tpol(jpiglo/2+1:jpiglo) = 0._wp
         fpol(     1    :jpiglo) = 0._wp
         IF( mjg(nlej) == jpjglo ) THEN                  ! only half of the nlcj-1 row
            DO ji = iif+1, iil-1
               tmask_i(ji,nlej-1) = tmask_i(ji,nlej-1) * tpol(mig(ji))
            END DO
         ENDIF
      ENDIF
      IF( jperio == 5 .OR. jperio == 6 ) THEN      ! F-point pivot
         tpol(     1    :jpiglo) = 0._wp
         fpol(jpiglo/2+1:jpiglo) = 0._wp
      ENDIF

      ! 2. Ocean/land mask at u-,  v-, and z-points (computed from tmask)
      ! -------------------------------------------
      DO jk = 1, jpk
         DO jj = 1, jpjm1
            DO ji = 1, fs_jpim1   ! vector loop
               umask(ji,jj,jk) = tmask(ji,jj  ,jk) * tmask(ji+1,jj  ,jk)
               vmask(ji,jj,jk) = tmask(ji,jj  ,jk) * tmask(ji  ,jj+1,jk)
            END DO
            DO ji = 1, jpim1      ! NO vector opt.
               fmask(ji,jj,jk) = tmask(ji,jj  ,jk) * tmask(ji+1,jj  ,jk)   &
                  &            * tmask(ji,jj+1,jk) * tmask(ji+1,jj+1,jk)
            END DO
         END DO
      END DO
      ! (ISF) MIN(1,SUM(umask)) is here to check if you have effectively at least 1 wet u point
      DO jj = 1, jpjm1
         DO ji = 1, fs_jpim1   ! vector loop
            umask_i(ji,jj)  = ssmask(ji,jj) * ssmask(ji+1,jj  )  * MIN(1._wp,SUM(umask(ji,jj,:)))
            vmask_i(ji,jj)  = ssmask(ji,jj) * ssmask(ji  ,jj+1)  * MIN(1._wp,SUM(vmask(ji,jj,:)))
         END DO
         DO ji = 1, jpim1      ! NO vector opt.
            fmask_i(ji,jj) =  ssmask(ji,jj  ) * ssmask(ji+1,jj  )   &
               &            * ssmask(ji,jj+1) * ssmask(ji+1,jj+1) * MIN(1._wp,SUM(fmask(ji,jj,:)))
         END DO
      END DO
      CALL lbc_lnk( umask, 'U', 1._wp )      ! Lateral boundary conditions
      CALL lbc_lnk( vmask, 'V', 1._wp )
      CALL lbc_lnk( fmask, 'F', 1._wp )
      CALL lbc_lnk( umask_i, 'U', 1._wp )      ! Lateral boundary conditions
      CALL lbc_lnk( vmask_i, 'V', 1._wp )
      CALL lbc_lnk( fmask_i, 'F', 1._wp )

      ! 3. Ocean/land mask at wu-, wv- and w points 
      !----------------------------------------------
      wmask (:,:,1) = tmask(:,:,1) ! ????????
      wumask(:,:,1) = umask(:,:,1) ! ????????
      wvmask(:,:,1) = vmask(:,:,1) ! ????????
      DO jk=2,jpk
         wmask (:,:,jk)=tmask(:,:,jk) * tmask(:,:,jk-1)
         wumask(:,:,jk)=umask(:,:,jk) * umask(:,:,jk-1)   
         wvmask(:,:,jk)=vmask(:,:,jk) * vmask(:,:,jk-1)
      END DO

      ! 4. ocean/land mask for the elliptic equation
      ! --------------------------------------------
      bmask(:,:) = ssmask(:,:)       ! elliptic equation is written at t-point
      !
      !                               ! Boundary conditions
      !                                    ! cyclic east-west : bmask must be set to 0. on rows 1 and jpi
      IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN
         bmask( 1 ,:) = 0._wp
         bmask(jpi,:) = 0._wp
      ENDIF
      IF( nperio == 2 ) THEN               ! south symmetric :  bmask must be set to 0. on row 1
         bmask(:, 1 ) = 0._wp
      ENDIF
      !                                    ! north fold : 
      IF( nperio == 3 .OR. nperio == 4 ) THEN   ! T-pt pivot : bmask set to 0. on row jpj and on half jpjglo-1 row
         DO ji = 1, jpi                      
            ii = ji + nimpp - 1
            bmask(ji,jpj-1) = bmask(ji,jpj-1) * tpol(ii)
            bmask(ji,jpj  ) = 0._wp
         END DO
      ENDIF
      IF( nperio == 5 .OR. nperio == 6 ) THEN   ! F-pt pivot and T-pt elliptic eq. : bmask set to 0. on row jpj
         bmask(:,jpj) = 0._wp
      ENDIF
      !
      IF( lk_mpp ) THEN                    ! mpp specificities
         !                                      ! bmask is set to zero on the overlap region
         IF( nbondi /= -1 .AND. nbondi /= 2 )   bmask(  1 :jpreci,:) = 0._wp
         IF( nbondi /=  1 .AND. nbondi /= 2 )   bmask(nlci:jpi   ,:) = 0._wp
         IF( nbondj /= -1 .AND. nbondj /= 2 )   bmask(:,  1 :jprecj) = 0._wp
         IF( nbondj /=  1 .AND. nbondj /= 2 )   bmask(:,nlcj:jpj   ) = 0._wp
         !
         IF( npolj == 3 .OR. npolj == 4 ) THEN  ! north fold : bmask must be set to 0. on rows jpj-1 and jpj
            DO ji = 1, nlci
               ii = ji + nimpp - 1
               bmask(ji,nlcj-1) = bmask(ji,nlcj-1) * tpol(ii)
               bmask(ji,nlcj  ) = 0._wp
            END DO
         ENDIF
         IF( npolj == 5 .OR. npolj == 6 ) THEN  ! F-pt pivot and T-pt elliptic eq. : bmask set to 0. on row jpj
            DO ji = 1, nlci
               bmask(ji,nlcj  ) = 0._wp
            END DO
         ENDIF
      ENDIF


      ! mask for second order calculation of vorticity
      ! ----------------------------------------------
      CALL dom_msk_nsa

      
      ! Lateral boundary conditions on velocity (modify fmask)
      ! ---------------------------------------     
      DO jk = 1, jpk
         zwf(:,:) = fmask(:,:,jk)         
         DO jj = 2, jpjm1
            DO ji = fs_2, fs_jpim1   ! vector opt.
               IF( fmask(ji,jj,jk) == 0._wp ) THEN
                  fmask(ji,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,jj), zwf(ji,jj+1),   &
                     &                                           zwf(ji-1,jj), zwf(ji,jj-1)  )  )
               ENDIF
            END DO
         END DO
         DO jj = 2, jpjm1
            IF( fmask(1,jj,jk) == 0._wp ) THEN
               fmask(1  ,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1) ) )
            ENDIF
            IF( fmask(jpi,jj,jk) == 0._wp ) THEN
               fmask(jpi,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) )
            ENDIF
         END DO         
         DO ji = 2, jpim1
            IF( fmask(ji,1,jk) == 0._wp ) THEN
               fmask(ji, 1 ,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1) ) )
            ENDIF
            IF( fmask(ji,jpj,jk) == 0._wp ) THEN
               fmask(ji,jpj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpjm1) ) )
            ENDIF
         END DO
#if defined key_agrif
         IF( .NOT. AGRIF_Root() ) THEN
            IF ((nbondi ==  1).OR.(nbondi == 2)) fmask(nlci-1 , :     ,jk) = 0.e0      ! east
            IF ((nbondi == -1).OR.(nbondi == 2)) fmask(1      , :     ,jk) = 0.e0      ! west
            IF ((nbondj ==  1).OR.(nbondj == 2)) fmask(:      ,nlcj-1 ,jk) = 0.e0      ! north
            IF ((nbondj == -1).OR.(nbondj == 2)) fmask(:      ,1      ,jk) = 0.e0      ! south
         ENDIF
#endif
      END DO
      !
      IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN   ! ORCA_R2 configuration
         !                                                 ! Increased lateral friction near of some straits
         IF( nn_cla == 0 ) THEN
            !                                ! Gibraltar strait  : partial slip (fmask=0.5)
            ij0 = 101   ;   ij1 = 101
            ii0 = 139   ;   ii1 = 140   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  0.5_wp
            ij0 = 102   ;   ij1 = 102
            ii0 = 139   ;   ii1 = 140   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  0.5_wp
            !
            !                                ! Bab el Mandeb : partial slip (fmask=1)
            ij0 =  87   ;   ij1 =  88
            ii0 = 160   ;   ii1 = 160   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  1._wp
            ij0 =  88   ;   ij1 =  88
            ii0 = 159   ;   ii1 = 159   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  1._wp
            !
         ENDIF
         !                                ! Danish straits  : strong slip (fmask > 2)
! We keep this as an example but it is instable in this case 
!         ij0 = 115   ;   ij1 = 115
!         ii0 = 145   ;   ii1 = 146   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 4._wp
!         ij0 = 116   ;   ij1 = 116
!         ii0 = 145   ;   ii1 = 146   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 4._wp
         !
      ENDIF
      !
      IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN   ! ORCA R1 configuration
         !                                                 ! Increased lateral friction near of some straits
         ! This dirty section will be suppressed by simplification process:
         ! all this will come back in input files
         ! Currently these hard-wired indices relate to configuration with
         ! extend grid (jpjglo=332)
         !
         isrow = 332 - jpjglo
         !
         IF(lwp) WRITE(numout,*)
         IF(lwp) WRITE(numout,*) '   orca_r1: increase friction near the following straits : '
         IF(lwp) WRITE(numout,*) '      Gibraltar '
         ii0 = 282           ;   ii1 = 283        ! Gibraltar Strait 
         ij0 = 241 - isrow   ;   ij1 = 241 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp  

         IF(lwp) WRITE(numout,*) '      Bhosporus '
         ii0 = 314           ;   ii1 = 315        ! Bhosporus Strait 
         ij0 = 248 - isrow   ;   ij1 = 248 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp  

         IF(lwp) WRITE(numout,*) '      Makassar (Top) '
         ii0 =  48           ;   ii1 =  48        ! Makassar Strait (Top) 
         ij0 = 189 - isrow   ;   ij1 = 190 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp  

         IF(lwp) WRITE(numout,*) '      Lombok '
         ii0 =  44           ;   ii1 =  44        ! Lombok Strait 
         ij0 = 164 - isrow   ;   ij1 = 165 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp  

         IF(lwp) WRITE(numout,*) '      Ombai '
         ii0 =  53           ;   ii1 =  53        ! Ombai Strait 
         ij0 = 164 - isrow   ;   ij1 = 165 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp  

         IF(lwp) WRITE(numout,*) '      Timor Passage '
         ii0 =  56           ;   ii1 =  56        ! Timor Passage 
         ij0 = 164 - isrow   ;   ij1 = 165 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp  

         IF(lwp) WRITE(numout,*) '      West Halmahera '
         ii0 =  58           ;   ii1 =  58        ! West Halmahera Strait 
         ij0 = 181 - isrow   ;   ij1 = 182 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp  

         IF(lwp) WRITE(numout,*) '      East Halmahera '
         ii0 =  55           ;   ii1 =  55        ! East Halmahera Strait 
         ij0 = 181 - isrow   ;   ij1 = 182 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp  
         !
      ENDIF
      !
      CALL lbc_lnk( fmask, 'F', 1._wp )      ! Lateral boundary conditions on fmask

      ! CAUTION : The fmask may be further modified in dyn_vor_init ( dynvor.F90 )
            
      IF( nprint == 1 .AND. lwp ) THEN      ! Control print
         imsk(:,:) = INT( tmask_i(:,:) )
         WRITE(numout,*) ' tmask_i : '
         CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1,   &
               &                           1, jpj, 1, 1, numout)
         WRITE (numout,*)
         WRITE (numout,*) ' dommsk: tmask for each level'
         WRITE (numout,*) ' ----------------------------'
         DO jk = 1, jpk
            imsk(:,:) = INT( tmask(:,:,jk) )

            WRITE(numout,*)
            WRITE(numout,*) ' level = ',jk
            CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1,   &
               &                              1, jpj, 1, 1, numout)
         END DO
         WRITE(numout,*)
         WRITE(numout,*) ' dom_msk: vmask for each level'
         WRITE(numout,*) ' -----------------------------'
         DO jk = 1, jpk
            imsk(:,:) = INT( vmask(:,:,jk) )
            WRITE(numout,*)
            WRITE(numout,*) ' level = ',jk
            CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1,   &
               &                              1, jpj, 1, 1, numout)
         END DO
         WRITE(numout,*)
         WRITE(numout,*) ' dom_msk: fmask for each level'
         WRITE(numout,*) ' -----------------------------'
         DO jk = 1, jpk
            imsk(:,:) = INT( fmask(:,:,jk) )
            WRITE(numout,*)
            WRITE(numout,*) ' level = ',jk
            CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1,   &
               &                              1, jpj, 1, 1, numout )
         END DO
         WRITE(numout,*)
         WRITE(numout,*) ' dom_msk: bmask '
         WRITE(numout,*) ' ---------------'
         WRITE(numout,*)
         imsk(:,:) = INT( bmask(:,:) )
         CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1,   &
            &                              1, jpj, 1, 1, numout )
      ENDIF
      !
      CALL wrk_dealloc( jpi, jpj, imsk )
      CALL wrk_dealloc( jpi, jpj, zwf  )
      !
      IF( nn_timing == 1 )  CALL timing_stop('dom_msk')
      !
   END SUBROUTINE dom_msk

#if defined key_noslip_accurate
   !!----------------------------------------------------------------------
   !!   'key_noslip_accurate' :         accurate no-slip boundary condition
   !!----------------------------------------------------------------------
   
   SUBROUTINE dom_msk_nsa
      !!---------------------------------------------------------------------
      !!                 ***  ROUTINE dom_msk_nsa  ***
      !! 
      !! ** Purpose :
      !!
      !! ** Method  :
      !!
      !! ** Action :
      !!----------------------------------------------------------------------
      INTEGER  ::   ji, jj, jk, jl      ! dummy loop indices
      INTEGER  ::   ine, inw, ins, inn, itest, ierror, iind, ijnd
      REAL(wp) ::   zaa
      !!---------------------------------------------------------------------
      !
      IF( nn_timing == 1 )  CALL timing_start('dom_msk_nsa')
      !
      IF(lwp) WRITE(numout,*)
      IF(lwp) WRITE(numout,*) 'dom_msk_nsa : noslip accurate boundary condition'
      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   using Schchepetkin and O Brian scheme'
      IF( lk_mpp )   CALL ctl_stop( ' mpp version is not yet implemented' )

      ! mask for second order calculation of vorticity
      ! ----------------------------------------------
      ! noslip boundary condition: fmask=1  at convex corner, store
      ! index of straight coast meshes ( 'west', refering to a coast,
      ! means west of the ocean, aso)
      
      DO jk = 1, jpk
         DO jl = 1, 4
            npcoa(jl,jk) = 0
            DO ji = 1, 2*(jpi+jpj)
               nicoa(ji,jl,jk) = 0
               njcoa(ji,jl,jk) = 0
            END DO
         END DO
      END DO
      
      IF( jperio == 2 ) THEN
         WRITE(numout,*) ' '
         WRITE(numout,*) ' symetric boundary conditions need special'
         WRITE(numout,*) ' treatment not implemented. we stop.'
         STOP
      ENDIF
      
      ! convex corners
      
      DO jk = 1, jpkm1
         DO jj = 1, jpjm1
            DO ji = 1, jpim1
               zaa = tmask(ji  ,jj,jk) + tmask(ji  ,jj+1,jk)   &
                  &+ tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk)
               IF( ABS(zaa-3._wp) <= 0.1_wp )   fmask(ji,jj,jk) = 1._wp
            END DO
         END DO
      END DO

      ! north-south straight coast

      DO jk = 1, jpkm1
         inw = 0
         ine = 0
         DO jj = 2, jpjm1
            DO ji = 2, jpim1
               zaa = tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk)
               IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN
                  inw = inw + 1
                  nicoa(inw,1,jk) = ji
                  njcoa(inw,1,jk) = jj
                  IF( nprint == 1 ) WRITE(numout,*) ' west  : ', jk, inw, ji, jj
               ENDIF
               zaa = tmask(ji,jj,jk) + tmask(ji,jj+1,jk)
               IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN
                  ine = ine + 1
                  nicoa(ine,2,jk) = ji
                  njcoa(ine,2,jk) = jj
                  IF( nprint == 1 ) WRITE(numout,*) ' east  : ', jk, ine, ji, jj
               ENDIF
            END DO
         END DO
         npcoa(1,jk) = inw
         npcoa(2,jk) = ine
      END DO

      ! west-east straight coast

      DO jk = 1, jpkm1
         ins = 0
         inn = 0
         DO jj = 2, jpjm1
            DO ji =2, jpim1
               zaa = tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk)
               IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN
                  ins = ins + 1
                  nicoa(ins,3,jk) = ji
                  njcoa(ins,3,jk) = jj
                  IF( nprint == 1 ) WRITE(numout,*) ' south : ', jk, ins, ji, jj
               ENDIF
               zaa = tmask(ji+1,jj,jk) + tmask(ji,jj,jk)
               IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN
                  inn = inn + 1
                  nicoa(inn,4,jk) = ji
                  njcoa(inn,4,jk) = jj
                  IF( nprint == 1 ) WRITE(numout,*) ' north : ', jk, inn, ji, jj
               ENDIF
            END DO
         END DO
         npcoa(3,jk) = ins
         npcoa(4,jk) = inn
      END DO

      itest = 2 * ( jpi + jpj )
      DO jk = 1, jpk
         IF( npcoa(1,jk) > itest .OR. npcoa(2,jk) > itest .OR.   &
             npcoa(3,jk) > itest .OR. npcoa(4,jk) > itest ) THEN
            
            WRITE(ctmp1,*) ' level jk = ',jk
            WRITE(ctmp2,*) ' straight coast index arraies are too small.:'
            WRITE(ctmp3,*) ' npe, npw, nps, npn = ', npcoa(1,jk), npcoa(2,jk),   &
                &                                     npcoa(3,jk), npcoa(4,jk)
            WRITE(ctmp4,*) ' 2*(jpi+jpj) = ',itest,'. we stop.'
            CALL ctl_stop( ctmp1, ctmp2, ctmp3, ctmp4 )
        ENDIF
      END DO

      ierror = 0
      iind = 0
      ijnd = 0
      IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 )   iind = 2
      IF( nperio == 3 .OR. nperio == 4 .OR. nperio == 5 .OR. nperio == 6 )   ijnd = 2
      DO jk = 1, jpk
         DO jl = 1, npcoa(1,jk)
            IF( nicoa(jl,1,jk)+3 > jpi+iind ) THEN
               ierror = ierror+1
               icoord(ierror,1) = nicoa(jl,1,jk)
               icoord(ierror,2) = njcoa(jl,1,jk)
               icoord(ierror,3) = jk
            ENDIF
         END DO
         DO jl = 1, npcoa(2,jk)
            IF(nicoa(jl,2,jk)-2 < 1-iind ) THEN
               ierror = ierror + 1
               icoord(ierror,1) = nicoa(jl,2,jk)
               icoord(ierror,2) = njcoa(jl,2,jk)
               icoord(ierror,3) = jk
            ENDIF
         END DO
         DO jl = 1, npcoa(3,jk)
            IF( njcoa(jl,3,jk)+3 > jpj+ijnd ) THEN
               ierror = ierror + 1
               icoord(ierror,1) = nicoa(jl,3,jk)
               icoord(ierror,2) = njcoa(jl,3,jk)
               icoord(ierror,3) = jk
            ENDIF
         END DO
         DO jl = 1, npcoa(4,jk)
            IF( njcoa(jl,4,jk)-2 < 1) THEN
               ierror=ierror + 1
               icoord(ierror,1) = nicoa(jl,4,jk)
               icoord(ierror,2) = njcoa(jl,4,jk)
               icoord(ierror,3) = jk
            ENDIF
         END DO
      END DO
      
      IF( ierror > 0 ) THEN
         IF(lwp) WRITE(numout,*)
         IF(lwp) WRITE(numout,*) '              Problem on lateral conditions'
         IF(lwp) WRITE(numout,*) '                 Bad marking off at points:'
         DO jl = 1, ierror
            IF(lwp) WRITE(numout,*) 'Level:',icoord(jl,3),   &
               &                  '  Point(',icoord(jl,1),',',icoord(jl,2),')'
         END DO
         CALL ctl_stop( 'We stop...' )
      ENDIF
      !
      IF( nn_timing == 1 )  CALL timing_stop('dom_msk_nsa')
      !
   END SUBROUTINE dom_msk_nsa

#else
   !!----------------------------------------------------------------------
   !!   Default option :                                      Empty routine
   !!----------------------------------------------------------------------
   SUBROUTINE dom_msk_nsa       
   END SUBROUTINE dom_msk_nsa
#endif
   
   !!======================================================================
END MODULE dommsk