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

source: trunk/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90 @ 6125

Last change on this file since 6125 was 6125, checked in by jchanut, 8 years ago

suppress bmask, #1648

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