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_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

source: branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90 @ 6004

Last change on this file since 6004 was 6004, checked in by gm, 8 years ago

#1613: vvl by default, step III: Merge with the trunk (free surface simplification) (see wiki)
  • Property svn:keywords set to Id
File size: 17.5 KB
Line 
1MODULE dommsk
2   !!======================================================================
3   !!                       ***  MODULE dommsk   ***
4   !! Ocean initialization : domain land/sea mask
5   !!======================================================================
6   !! History :  OPA  ! 1987-07  (G. Madec)  Original code
7   !!            6.0  ! 1993-03  (M. Guyon)  symetrical conditions (M. Guyon)
8   !!            7.0  ! 1996-01  (G. Madec)  suppression of common work arrays
9   !!             -   ! 1996-05  (G. Madec)  mask computed from tmask
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
12   !!             -   ! 1998-05  (G. Roullet)  free surface
13   !!            8.2  ! 2000-03  (G. Madec)  no slip accurate
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
18   !!            3.6  ! 2015-05  (P. Mathiot) ISF: add wmask,wumask and wvmask
19   !!----------------------------------------------------------------------
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
26   !
27   USE in_out_manager  ! I/O manager
28   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
29   USE lib_mpp         !
30   USE wrk_nemo        ! Memory allocation
31   USE timing          ! Timing
32
33   IMPLICIT NONE
34   PRIVATE
35
36   PUBLIC   dom_msk    ! routine called by inidom.F90
37
38   !                            !!* Namelist namlbc : lateral boundary condition *
39   REAL(wp)        :: rn_shlat   ! type of lateral boundary condition on velocity
40   LOGICAL, PUBLIC :: ln_vorlat  !  consistency of vorticity boundary condition
41   !                                            with analytical eqs.
42
43   !! * Substitutions
44#  include "vectopt_loop_substitute.h90"
45   !!----------------------------------------------------------------------
46   !! NEMO/OPA 3.2 , LODYC-IPSL  (2009)
47   !! $Id$
48   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
49   !!----------------------------------------------------------------------
50CONTAINS
51
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.
62      !!      mbathy equals 0 over continental T-point
63      !!      and the number of ocean level over the ocean.
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)
75      !!                  and mbathy(ji+1,jj) and mbathy(ji+1,jj+1) >= jk.
76      !!      b-point : the same definition as for f-point of the first ocean
77      !!                level (surface level) but with 0 along coastlines.
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.
81      !!
82      !!        The lateral friction is set through the value of fmask along
83      !!      the coast and topography. This value is defined by rn_shlat, a
84      !!      namelist parameter:
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
89      !!
90      !!      N.B. If nperio not equal to 0, the land/ocean mask arrays
91      !!      are defined with the proper value at lateral domain boundaries.
92      !!
93      !!      In case of open boundaries (lk_bdy=T):
94      !!        - tmask is set to 1 on the points to be computed bay the open
95      !!          boundaries routines.
96      !!
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.)
101      !!                          =rn_shlat along lateral boundaries
102      !!               tmask_i  : interior ocean mask
103      !!----------------------------------------------------------------------
104      INTEGER  ::   ji, jj, jk               ! dummy loop indices
105      INTEGER  ::   iif, iil, ii0, ii1, ii   ! local integers
106      INTEGER  ::   ijf, ijl, ij0, ij1       !   -       -
107      INTEGER  ::   ios
108      INTEGER  ::   isrow                    ! index for ORCA1 starting row
109      INTEGER , POINTER, DIMENSION(:,:) ::  imsk
110      REAL(wp), POINTER, DIMENSION(:,:) ::  zwf
111      !!
112      NAMELIST/namlbc/ rn_shlat, ln_vorlat
113      !!---------------------------------------------------------------------
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      !
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 )
127      IF(lwm) WRITE ( numond, namlbc )
128     
129      IF(lwp) THEN                  ! control print
130         WRITE(numout,*)
131         WRITE(numout,*) 'dommsk : ocean mask '
132         WRITE(numout,*) '~~~~~~'
133         WRITE(numout,*) '   Namelist namlbc'
134         WRITE(numout,*) '      lateral momentum boundary cond.    rn_shlat  = ',rn_shlat
135         WRITE(numout,*) '      consistency with analytical form   ln_vorlat = ',ln_vorlat 
136      ENDIF
137
138      IF     (      rn_shlat == 0.               ) THEN   ;   IF(lwp) WRITE(numout,*) '   ocean lateral  free-slip '
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 )
145      ENDIF
146
147      ! 1. Ocean/land mask at t-point (computed from mbathy)
148      ! -----------------------------
149      ! N.B. tmask has already the right boundary conditions since mbathy is ok
150      !
151      tmask(:,:,:) = 0._wp
152      DO jk = 1, jpk
153         DO jj = 1, jpj
154            DO ji = 1, jpi
155               IF( REAL( mbathy(ji,jj) - jk, wp ) + 0.1_wp >= 0._wp )   tmask(ji,jj,jk) = 1._wp
156            END DO 
157         END DO 
158      END DO 
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 
172
173      ! Interior domain mask (used for global sum)
174      ! --------------------
175      tmask_i(:,:) = ssmask(:,:)            ! (ISH) tmask_i = 1 even on the ice shelf
176      iif = jpreci                         ! ???
177      iil = nlci - jpreci + 1
178      ijf = jprecj                         ! ???
179      ijl = nlcj - jprecj + 1
180
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)
185
186      ! north fold mask
187      ! ---------------
188      tpol(1:jpiglo) = 1._wp 
189      fpol(1:jpiglo) = 1._wp
190      IF( jperio == 3 .OR. jperio == 4 ) THEN      ! T-point pivot
191         tpol(jpiglo/2+1:jpiglo) = 0._wp
192         fpol(     1    :jpiglo) = 0._wp
193         IF( mjg(nlej) == jpjglo ) THEN                  ! only half of the nlcj-1 row
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
198      ENDIF
199      IF( jperio == 5 .OR. jperio == 6 ) THEN      ! F-point pivot
200         tpol(     1    :jpiglo) = 0._wp
201         fpol(jpiglo/2+1:jpiglo) = 0._wp
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)
211            END DO
212            DO ji = 1, jpim1      ! NO vector opt.
213               fmask(ji,jj,jk) = tmask(ji,jj  ,jk) * tmask(ji+1,jj  ,jk)   &
214                  &            * tmask(ji,jj+1,jk) * tmask(ji+1,jj+1,jk)
215            END DO
216         END DO
217      END DO
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
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 )
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 )
235
236      ! 3. Ocean/land mask at wu-, wv- and w points
237      !----------------------------------------------
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)
245      END DO
246
247      ! Lateral boundary conditions on velocity (modify fmask)
248      ! ---------------------------------------     
249      DO jk = 1, jpk
250         zwf(:,:) = fmask(:,:,jk)         
251         DO jj = 2, jpjm1
252            DO ji = fs_2, fs_jpim1   ! vector opt.
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)  )  )
256               ENDIF
257            END DO
258         END DO
259         DO jj = 2, jpjm1
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) ) )
262            ENDIF
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) ) )
265            ENDIF
266         END DO         
267         DO ji = 2, jpim1
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) ) )
270            ENDIF
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) ) )
273            ENDIF
274         END DO
275      END DO
276      !
277      IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN   ! ORCA_R2 configuration
278         !                                                 ! Increased lateral friction near of some straits
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         !
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
294!         ii0 = 145   ;   ii1 = 146   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 4._wp
295!         ij0 = 116   ;   ij1 = 116
296!         ii0 = 145   ;   ii1 = 146   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 4._wp
297         !
298      ENDIF
299      !
300      IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN   ! ORCA R1 configuration
301         !                                                 ! Increased lateral friction near of some straits
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         !
309         IF(lwp) WRITE(numout,*)
310         IF(lwp) WRITE(numout,*) '   orca_r1: increase friction near the following straits : '
311         IF(lwp) WRITE(numout,*) '      Gibraltar '
312         ii0 = 282           ;   ii1 = 283        ! Gibraltar Strait
313         ij0 = 241 - isrow   ;   ij1 = 241 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp 
314
315         IF(lwp) WRITE(numout,*) '      Bhosporus '
316         ii0 = 314           ;   ii1 = 315        ! Bhosporus Strait
317         ij0 = 248 - isrow   ;   ij1 = 248 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp 
318
319         IF(lwp) WRITE(numout,*) '      Makassar (Top) '
320         ii0 =  48           ;   ii1 =  48        ! Makassar Strait (Top)
321         ij0 = 189 - isrow   ;   ij1 = 190 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp 
322
323         IF(lwp) WRITE(numout,*) '      Lombok '
324         ii0 =  44           ;   ii1 =  44        ! Lombok Strait
325         ij0 = 164 - isrow   ;   ij1 = 165 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp 
326
327         IF(lwp) WRITE(numout,*) '      Ombai '
328         ii0 =  53           ;   ii1 =  53        ! Ombai Strait
329         ij0 = 164 - isrow   ;   ij1 = 165 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp 
330
331         IF(lwp) WRITE(numout,*) '      Timor Passage '
332         ii0 =  56           ;   ii1 =  56        ! Timor Passage
333         ij0 = 164 - isrow   ;   ij1 = 165 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp 
334
335         IF(lwp) WRITE(numout,*) '      West Halmahera '
336         ii0 =  58           ;   ii1 =  58        ! West Halmahera Strait
337         ij0 = 181 - isrow   ;   ij1 = 182 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp 
338
339         IF(lwp) WRITE(numout,*) '      East Halmahera '
340         ii0 =  55           ;   ii1 =  55        ! East Halmahera Strait
341         ij0 = 181 - isrow   ;   ij1 = 182 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp 
342         !
343      ENDIF
344      !
345      CALL lbc_lnk( fmask, 'F', 1._wp )      ! Lateral boundary conditions on fmask
346      !
347      ! CAUTION : The fmask may be further modified in dyn_vor_init ( dynvor.F90 )
348      !
349      CALL wrk_dealloc( jpi, jpj, imsk )
350      CALL wrk_dealloc( jpi, jpj, zwf  )
351      !
352      IF( nn_timing == 1 )  CALL timing_stop('dom_msk')
353      !
354   END SUBROUTINE dom_msk
355   
356   !!======================================================================
357END MODULE dommsk
Note: See TracBrowser for help on using the repository browser.