New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
dommsk.F90 in branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

source: branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90 @ 6060

Last change on this file since 6060 was 6060, checked in by timgraham, 9 years ago

Merged dev_r5836_noc2_VVL_BY_DEFAULT into branch

  • Property svn:keywords set to Id
File size: 17.5 KB
RevLine 
[3]1MODULE dommsk
[1566]2   !!======================================================================
[3]3   !!                       ***  MODULE dommsk   ***
4   !! Ocean initialization : domain land/sea mask
[1566]5   !!======================================================================
6   !! History :  OPA  ! 1987-07  (G. Madec)  Original code
[2528]7   !!            6.0  ! 1993-03  (M. Guyon)  symetrical conditions (M. Guyon)
8   !!            7.0  ! 1996-01  (G. Madec)  suppression of common work arrays
[6060]9   !!             -   ! 1996-05  (G. Madec)  mask computed from tmask
[2528]10   !!            8.0  ! 1997-02  (G. Madec)  mesh information put in domhgr.F
11   !!            8.1  ! 1997-07  (G. Madec)  modification of mbathy and fmask
[1566]12   !!             -   ! 1998-05  (G. Roullet)  free surface
[2528]13   !!            8.2  ! 2000-03  (G. Madec)  no slip accurate
[1566]14   !!             -   ! 2001-09  (J.-M. Molines)  Open boundaries
15   !!   NEMO     1.0  ! 2002-08  (G. Madec)  F90: Free form and module
16   !!             -   ! 2005-11  (V. Garnier) Surface pressure gradient organization
17   !!            3.2  ! 2009-07  (R. Benshila) Suppression of rigid-lid option
[5836]18   !!            3.6  ! 2015-05  (P. Mathiot) ISF: add wmask,wumask and wvmask
[1566]19   !!----------------------------------------------------------------------
[3]20
21   !!----------------------------------------------------------------------
22   !!   dom_msk        : compute land/ocean mask
23   !!----------------------------------------------------------------------
24   USE oce             ! ocean dynamics and tracers
25   USE dom_oce         ! ocean space and time domain
[6060]26   !
[3]27   USE in_out_manager  ! I/O manager
28   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
[6060]29   USE lib_mpp         !
[3294]30   USE wrk_nemo        ! Memory allocation
31   USE timing          ! Timing
[3]32
33   IMPLICIT NONE
34   PRIVATE
35
[6060]36   PUBLIC   dom_msk    ! routine called by inidom.F90
[3]37
[1601]38   !                            !!* Namelist namlbc : lateral boundary condition *
[4147]39   REAL(wp)        :: rn_shlat   ! type of lateral boundary condition on velocity
40   LOGICAL, PUBLIC :: ln_vorlat  !  consistency of vorticity boundary condition
[3294]41   !                                            with analytical eqs.
[2715]42
[3]43   !! * Substitutions
44#  include "vectopt_loop_substitute.h90"
[1566]45   !!----------------------------------------------------------------------
46   !! NEMO/OPA 3.2 , LODYC-IPSL  (2009)
[1152]47   !! $Id$
[2528]48   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
[1566]49   !!----------------------------------------------------------------------
[3]50CONTAINS
[2715]51
[3]52   SUBROUTINE dom_msk
53      !!---------------------------------------------------------------------
54      !!                 ***  ROUTINE dom_msk  ***
55      !!
56      !! ** Purpose :   Compute land/ocean mask arrays at tracer points, hori-
57      !!      zontal velocity points (u & v), vorticity points (f) and baro-
58      !!      tropic stream function  points (b).
59      !!
60      !! ** Method  :   The ocean/land mask is computed from the basin bathy-
61      !!      metry in level (mbathy) which is defined or read in dommba.
[1528]62      !!      mbathy equals 0 over continental T-point
63      !!      and the number of ocean level over the ocean.
[3]64      !!
65      !!      At a given position (ji,jj,jk) the ocean/land mask is given by:
66      !!      t-point : 0. IF mbathy( ji ,jj) =< 0
67      !!                1. IF mbathy( ji ,jj) >= jk
68      !!      u-point : 0. IF mbathy( ji ,jj)  or mbathy(ji+1, jj ) =< 0
69      !!                1. IF mbathy( ji ,jj) and mbathy(ji+1, jj ) >= jk.
70      !!      v-point : 0. IF mbathy( ji ,jj)  or mbathy( ji ,jj+1) =< 0
71      !!                1. IF mbathy( ji ,jj) and mbathy( ji ,jj+1) >= jk.
72      !!      f-point : 0. IF mbathy( ji ,jj)  or mbathy( ji ,jj+1)
73      !!                   or mbathy(ji+1,jj)  or mbathy(ji+1,jj+1) =< 0
74      !!                1. IF mbathy( ji ,jj) and mbathy( ji ,jj+1)
[2528]75      !!                  and mbathy(ji+1,jj) and mbathy(ji+1,jj+1) >= jk.
[3]76      !!      b-point : the same definition as for f-point of the first ocean
77      !!                level (surface level) but with 0 along coastlines.
[2528]78      !!      tmask_i : interior ocean mask at t-point, i.e. excluding duplicated
79      !!                rows/lines due to cyclic or North Fold boundaries as well
80      !!                as MPP halos.
[3]81      !!
82      !!        The lateral friction is set through the value of fmask along
[1601]83      !!      the coast and topography. This value is defined by rn_shlat, a
[3]84      !!      namelist parameter:
[1601]85      !!         rn_shlat = 0, free slip  (no shear along the coast)
86      !!         rn_shlat = 2, no slip  (specified zero velocity at the coast)
87      !!         0 < rn_shlat < 2, partial slip   | non-linear velocity profile
88      !!         2 < rn_shlat, strong slip        | in the lateral boundary layer
[3]89      !!
90      !!      N.B. If nperio not equal to 0, the land/ocean mask arrays
[6060]91      !!      are defined with the proper value at lateral domain boundaries.
[3]92      !!
[4328]93      !!      In case of open boundaries (lk_bdy=T):
[3]94      !!        - tmask is set to 1 on the points to be computed bay the open
95      !!          boundaries routines.
96      !!
[1566]97      !! ** Action :   tmask    : land/ocean mask at t-point (=0. or 1.)
98      !!               umask    : land/ocean mask at u-point (=0. or 1.)
99      !!               vmask    : land/ocean mask at v-point (=0. or 1.)
100      !!               fmask    : land/ocean mask at f-point (=0. or 1.)
[1601]101      !!                          =rn_shlat along lateral boundaries
[2528]102      !!               tmask_i  : interior ocean mask
[3]103      !!----------------------------------------------------------------------
[5836]104      INTEGER  ::   ji, jj, jk               ! dummy loop indices
[2715]105      INTEGER  ::   iif, iil, ii0, ii1, ii   ! local integers
106      INTEGER  ::   ijf, ijl, ij0, ij1       !   -       -
[4147]107      INTEGER  ::   ios
[5385]108      INTEGER  ::   isrow                    ! index for ORCA1 starting row
[3294]109      INTEGER , POINTER, DIMENSION(:,:) ::  imsk
110      REAL(wp), POINTER, DIMENSION(:,:) ::  zwf
[1601]111      !!
[3294]112      NAMELIST/namlbc/ rn_shlat, ln_vorlat
[3]113      !!---------------------------------------------------------------------
[3294]114      !
115      IF( nn_timing == 1 )  CALL timing_start('dom_msk')
116      !
117      CALL wrk_alloc( jpi, jpj, imsk )
118      CALL wrk_alloc( jpi, jpj, zwf  )
119      !
[4147]120      REWIND( numnam_ref )              ! Namelist namlbc in reference namelist : Lateral momentum boundary condition
121      READ  ( numnam_ref, namlbc, IOSTAT = ios, ERR = 901 )
122901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namlbc in reference namelist', lwp )
123
124      REWIND( numnam_cfg )              ! Namelist namlbc in configuration namelist : Lateral momentum boundary condition
125      READ  ( numnam_cfg, namlbc, IOSTAT = ios, ERR = 902 )
126902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namlbc in configuration namelist', lwp )
[4624]127      IF(lwm) WRITE ( numond, namlbc )
[1566]128     
129      IF(lwp) THEN                  ! control print
[3]130         WRITE(numout,*)
131         WRITE(numout,*) 'dommsk : ocean mask '
132         WRITE(numout,*) '~~~~~~'
[1566]133         WRITE(numout,*) '   Namelist namlbc'
[3294]134         WRITE(numout,*) '      lateral momentum boundary cond.    rn_shlat  = ',rn_shlat
135         WRITE(numout,*) '      consistency with analytical form   ln_vorlat = ',ln_vorlat 
[3]136      ENDIF
137
[2528]138      IF     (      rn_shlat == 0.               ) THEN   ;   IF(lwp) WRITE(numout,*) '   ocean lateral  free-slip '
[1601]139      ELSEIF (      rn_shlat == 2.               ) THEN   ;   IF(lwp) WRITE(numout,*) '   ocean lateral  no-slip '
140      ELSEIF ( 0. < rn_shlat .AND. rn_shlat < 2. ) THEN   ;   IF(lwp) WRITE(numout,*) '   ocean lateral  partial-slip '
141      ELSEIF ( 2. < rn_shlat                     ) THEN   ;   IF(lwp) WRITE(numout,*) '   ocean lateral  strong-slip '
142      ELSE
143         WRITE(ctmp1,*) ' rn_shlat is negative = ', rn_shlat
144         CALL ctl_stop( ctmp1 )
[3]145      ENDIF
146
147      ! 1. Ocean/land mask at t-point (computed from mbathy)
148      ! -----------------------------
[1566]149      ! N.B. tmask has already the right boundary conditions since mbathy is ok
150      !
[2528]151      tmask(:,:,:) = 0._wp
[3]152      DO jk = 1, jpk
153         DO jj = 1, jpj
154            DO ji = 1, jpi
[2528]155               IF( REAL( mbathy(ji,jj) - jk, wp ) + 0.1_wp >= 0._wp )   tmask(ji,jj,jk) = 1._wp
[3]156            END DO 
157         END DO 
158      END DO 
[4990]159     
160      ! (ISF) define barotropic mask and mask the ice shelf point
161      ssmask(:,:)=tmask(:,:,1) ! at this stage ice shelf is not masked
162     
163      DO jk = 1, jpk
164         DO jj = 1, jpj
165            DO ji = 1, jpi
166               IF( REAL( misfdep(ji,jj) - jk, wp ) - 0.1_wp >= 0._wp )   THEN
167                  tmask(ji,jj,jk) = 0._wp
168               END IF
169            END DO 
170         END DO 
171      END DO 
[3]172
173      ! Interior domain mask (used for global sum)
174      ! --------------------
[4990]175      tmask_i(:,:) = ssmask(:,:)            ! (ISH) tmask_i = 1 even on the ice shelf
[3]176      iif = jpreci                         ! ???
177      iil = nlci - jpreci + 1
178      ijf = jprecj                         ! ???
179      ijl = nlcj - jprecj + 1
180
[2528]181      tmask_i( 1 :iif,   :   ) = 0._wp      ! first columns
182      tmask_i(iil:jpi,   :   ) = 0._wp      ! last  columns (including mpp extra columns)
183      tmask_i(   :   , 1 :ijf) = 0._wp      ! first rows
184      tmask_i(   :   ,ijl:jpj) = 0._wp      ! last  rows (including mpp extra rows)
[3]185
186      ! north fold mask
[1566]187      ! ---------------
[2528]188      tpol(1:jpiglo) = 1._wp 
189      fpol(1:jpiglo) = 1._wp
[3]190      IF( jperio == 3 .OR. jperio == 4 ) THEN      ! T-point pivot
[2528]191         tpol(jpiglo/2+1:jpiglo) = 0._wp
192         fpol(     1    :jpiglo) = 0._wp
[1566]193         IF( mjg(nlej) == jpjglo ) THEN                  ! only half of the nlcj-1 row
[291]194            DO ji = iif+1, iil-1
195               tmask_i(ji,nlej-1) = tmask_i(ji,nlej-1) * tpol(mig(ji))
196            END DO
197         ENDIF
[3]198      ENDIF
199      IF( jperio == 5 .OR. jperio == 6 ) THEN      ! F-point pivot
[2528]200         tpol(     1    :jpiglo) = 0._wp
201         fpol(jpiglo/2+1:jpiglo) = 0._wp
[3]202      ENDIF
203
204      ! 2. Ocean/land mask at u-,  v-, and z-points (computed from tmask)
205      ! -------------------------------------------
206      DO jk = 1, jpk
207         DO jj = 1, jpjm1
208            DO ji = 1, fs_jpim1   ! vector loop
209               umask(ji,jj,jk) = tmask(ji,jj  ,jk) * tmask(ji+1,jj  ,jk)
210               vmask(ji,jj,jk) = tmask(ji,jj  ,jk) * tmask(ji  ,jj+1,jk)
[1271]211            END DO
[1694]212            DO ji = 1, jpim1      ! NO vector opt.
[3]213               fmask(ji,jj,jk) = tmask(ji,jj  ,jk) * tmask(ji+1,jj  ,jk)   &
[62]214                  &            * tmask(ji,jj+1,jk) * tmask(ji+1,jj+1,jk)
[3]215            END DO
216         END DO
217      END DO
[4990]218      ! (ISF) MIN(1,SUM(umask)) is here to check if you have effectively at least 1 wet u point
219      DO jj = 1, jpjm1
220         DO ji = 1, fs_jpim1   ! vector loop
221            umask_i(ji,jj)  = ssmask(ji,jj) * ssmask(ji+1,jj  )  * MIN(1._wp,SUM(umask(ji,jj,:)))
222            vmask_i(ji,jj)  = ssmask(ji,jj) * ssmask(ji  ,jj+1)  * MIN(1._wp,SUM(vmask(ji,jj,:)))
223         END DO
224         DO ji = 1, jpim1      ! NO vector opt.
225            fmask_i(ji,jj) =  ssmask(ji,jj  ) * ssmask(ji+1,jj  )   &
226               &            * ssmask(ji,jj+1) * ssmask(ji+1,jj+1) * MIN(1._wp,SUM(fmask(ji,jj,:)))
227         END DO
228      END DO
[2528]229      CALL lbc_lnk( umask, 'U', 1._wp )      ! Lateral boundary conditions
230      CALL lbc_lnk( vmask, 'V', 1._wp )
231      CALL lbc_lnk( fmask, 'F', 1._wp )
[4990]232      CALL lbc_lnk( umask_i, 'U', 1._wp )      ! Lateral boundary conditions
233      CALL lbc_lnk( vmask_i, 'V', 1._wp )
234      CALL lbc_lnk( fmask_i, 'F', 1._wp )
[3]235
[5120]236      ! 3. Ocean/land mask at wu-, wv- and w points
237      !----------------------------------------------
[5836]238      wmask (:,:,1) = tmask(:,:,1)     ! surface
239      wumask(:,:,1) = umask(:,:,1)
240      wvmask(:,:,1) = vmask(:,:,1)
241      DO jk = 2, jpk                   ! interior values
242         wmask (:,:,jk) = tmask(:,:,jk) * tmask(:,:,jk-1)
243         wumask(:,:,jk) = umask(:,:,jk) * umask(:,:,jk-1)   
244         wvmask(:,:,jk) = vmask(:,:,jk) * vmask(:,:,jk-1)
[5120]245      END DO
[3]246
247      ! Lateral boundary conditions on velocity (modify fmask)
[1566]248      ! ---------------------------------------     
[3]249      DO jk = 1, jpk
[1566]250         zwf(:,:) = fmask(:,:,jk)         
[3]251         DO jj = 2, jpjm1
252            DO ji = fs_2, fs_jpim1   ! vector opt.
[2528]253               IF( fmask(ji,jj,jk) == 0._wp ) THEN
254                  fmask(ji,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,jj), zwf(ji,jj+1),   &
255                     &                                           zwf(ji-1,jj), zwf(ji,jj-1)  )  )
[3]256               ENDIF
257            END DO
258         END DO
259         DO jj = 2, jpjm1
[2528]260            IF( fmask(1,jj,jk) == 0._wp ) THEN
261               fmask(1  ,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1) ) )
[3]262            ENDIF
[2528]263            IF( fmask(jpi,jj,jk) == 0._wp ) THEN
264               fmask(jpi,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) )
[3]265            ENDIF
[1566]266         END DO         
[3]267         DO ji = 2, jpim1
[2528]268            IF( fmask(ji,1,jk) == 0._wp ) THEN
269               fmask(ji, 1 ,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1) ) )
[3]270            ENDIF
[2528]271            IF( fmask(ji,jpj,jk) == 0._wp ) THEN
272               fmask(ji,jpj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpjm1) ) )
[3]273            ENDIF
274         END DO
275      END DO
[1566]276      !
277      IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN   ! ORCA_R2 configuration
278         !                                                 ! Increased lateral friction near of some straits
[5836]279         !                                ! Gibraltar strait  : partial slip (fmask=0.5)
280         ij0 = 101   ;   ij1 = 101
281         ii0 = 139   ;   ii1 = 140   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  0.5_wp
282         ij0 = 102   ;   ij1 = 102
283         ii0 = 139   ;   ii1 = 140   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  0.5_wp
284         !
285         !                                ! Bab el Mandeb : partial slip (fmask=1)
286         ij0 =  87   ;   ij1 =  88
287         ii0 = 160   ;   ii1 = 160   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  1._wp
288         ij0 =  88   ;   ij1 =  88
289         ii0 = 159   ;   ii1 = 159   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  1._wp
290         !
[1707]291         !                                ! Danish straits  : strong slip (fmask > 2)
292! We keep this as an example but it is instable in this case
293!         ij0 = 115   ;   ij1 = 115
[2528]294!         ii0 = 145   ;   ii1 = 146   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 4._wp
[1707]295!         ij0 = 116   ;   ij1 = 116
[2528]296!         ii0 = 145   ;   ii1 = 146   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 4._wp
[3]297         !
298      ENDIF
[1566]299      !
[2528]300      IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN   ! ORCA R1 configuration
301         !                                                 ! Increased lateral friction near of some straits
[5506]302         ! This dirty section will be suppressed by simplification process:
303         ! all this will come back in input files
304         ! Currently these hard-wired indices relate to configuration with
305         ! extend grid (jpjglo=332)
306         !
307         isrow = 332 - jpjglo
308         !
[2528]309         IF(lwp) WRITE(numout,*)
310         IF(lwp) WRITE(numout,*) '   orca_r1: increase friction near the following straits : '
311         IF(lwp) WRITE(numout,*) '      Gibraltar '
[5385]312         ii0 = 282           ;   ii1 = 283        ! Gibraltar Strait
[5552]313         ij0 = 241 - isrow   ;   ij1 = 241 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp 
[1566]314
[2528]315         IF(lwp) WRITE(numout,*) '      Bhosporus '
[5385]316         ii0 = 314           ;   ii1 = 315        ! Bhosporus Strait
[5552]317         ij0 = 248 - isrow   ;   ij1 = 248 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp 
[2528]318
319         IF(lwp) WRITE(numout,*) '      Makassar (Top) '
[5385]320         ii0 =  48           ;   ii1 =  48        ! Makassar Strait (Top)
[5552]321         ij0 = 189 - isrow   ;   ij1 = 190 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp 
[2528]322
323         IF(lwp) WRITE(numout,*) '      Lombok '
[5385]324         ii0 =  44           ;   ii1 =  44        ! Lombok Strait
[5552]325         ij0 = 164 - isrow   ;   ij1 = 165 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp 
[2528]326
327         IF(lwp) WRITE(numout,*) '      Ombai '
[5385]328         ii0 =  53           ;   ii1 =  53        ! Ombai Strait
[5552]329         ij0 = 164 - isrow   ;   ij1 = 165 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp 
[2528]330
331         IF(lwp) WRITE(numout,*) '      Timor Passage '
[5385]332         ii0 =  56           ;   ii1 =  56        ! Timor Passage
[5552]333         ij0 = 164 - isrow   ;   ij1 = 165 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp 
[2528]334
335         IF(lwp) WRITE(numout,*) '      West Halmahera '
[5385]336         ii0 =  58           ;   ii1 =  58        ! West Halmahera Strait
[5552]337         ij0 = 181 - isrow   ;   ij1 = 182 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp 
[2528]338
339         IF(lwp) WRITE(numout,*) '      East Halmahera '
[5385]340         ii0 =  55           ;   ii1 =  55        ! East Halmahera Strait
[5552]341         ij0 = 181 - isrow   ;   ij1 = 182 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp 
[2528]342         !
343      ENDIF
344      !
345      CALL lbc_lnk( fmask, 'F', 1._wp )      ! Lateral boundary conditions on fmask
[6060]346      !
[3294]347      ! CAUTION : The fmask may be further modified in dyn_vor_init ( dynvor.F90 )
[1566]348      !
[3294]349      CALL wrk_dealloc( jpi, jpj, imsk )
350      CALL wrk_dealloc( jpi, jpj, zwf  )
[2715]351      !
[3294]352      IF( nn_timing == 1 )  CALL timing_stop('dom_msk')
353      !
[3]354   END SUBROUTINE dom_msk
355   
356   !!======================================================================
357END MODULE dommsk
Note: See TracBrowser for help on using the repository browser.