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

source: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90 @ 7508

Last change on this file since 7508 was 7508, checked in by mocavero, 8 years ago

changes on code duplication and workshare construct

  • Property svn:keywords set to Id
File size: 18.3 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) 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!$OMP PARALLEL
149!$OMP DO schedule(static) private(jk, jj, ji)
150      DO jk = 1, jpk
151         DO jj = 1, jpj
152            DO ji = 1, jpi
153               tmask(ji,jj,jk) = 0._wp
154               IF( REAL( mbathy(ji,jj) - jk, wp ) + 0.1_wp >= 0._wp )   tmask(ji,jj,jk) = 1._wp
155            END DO 
156         END DO 
157      END DO 
158     
159      ! (ISF) define barotropic mask and mask the ice shelf point
160!$OMP DO schedule(static) private(jj, ji)
161      DO jj = 1, jpj
162         DO ji = 1, jpi
163            ssmask(ji,jj)=tmask(ji,jj,1) ! at this stage ice shelf is not masked
164         END DO 
165      END DO 
166!$OMP DO schedule(static) private(jk, jj, ji)
167      DO jk = 1, jpk
168         DO jj = 1, jpj
169            DO ji = 1, jpi
170               IF( REAL( misfdep(ji,jj) - jk, wp ) - 0.1_wp >= 0._wp )   THEN
171                  tmask(ji,jj,jk) = 0._wp
172               END IF
173            END DO 
174         END DO 
175      END DO 
176
177      ! Interior domain mask (used for global sum)
178      ! --------------------
179!$OMP DO schedule(static) private(jj, ji)
180      DO jj = 1, jpj
181         DO ji = 1, jpi
182            tmask_i(ji,jj) = ssmask(ji,jj)            ! (ISH) tmask_i = 1 even on the ice shelf
183            tmask_h(ji,jj) = 1._wp                 ! 0 on the halo and 1 elsewhere
184         END DO 
185      END DO 
186!$OMP END PARALLEL
187      iif = jpreci                         ! ???
188      iil = nlci - jpreci + 1
189      ijf = jprecj                         ! ???
190      ijl = nlcj - jprecj + 1
191
192      tmask_h( 1 :iif,   :   ) = 0._wp      ! first columns
193      tmask_h(iil:jpi,   :   ) = 0._wp      ! last  columns (including mpp extra columns)
194      tmask_h(   :   , 1 :ijf) = 0._wp      ! first rows
195      tmask_h(   :   ,ijl:jpj) = 0._wp      ! last  rows (including mpp extra rows)
196
197      ! north fold mask
198      ! ---------------
199      tpol(1:jpiglo) = 1._wp 
200      fpol(1:jpiglo) = 1._wp
201      IF( jperio == 3 .OR. jperio == 4 ) THEN      ! T-point pivot
202         tpol(jpiglo/2+1:jpiglo) = 0._wp
203         fpol(     1    :jpiglo) = 0._wp
204         IF( mjg(nlej) == jpjglo ) THEN                  ! only half of the nlcj-1 row
205            DO ji = iif+1, iil-1
206               tmask_h(ji,nlej-1) = tmask_h(ji,nlej-1) * tpol(mig(ji))
207            END DO
208         ENDIF
209      ENDIF
210     
211!$OMP PARALLEL DO schedule(static) private(jj, ji)
212      DO jj = 1, jpj
213         DO ji = 1, jpi
214            tmask_i(ji,jj) = tmask_i(ji,jj) * tmask_h(ji,jj)
215         END DO
216      END DO
217
218      IF( jperio == 5 .OR. jperio == 6 ) THEN      ! F-point pivot
219         tpol(     1    :jpiglo) = 0._wp
220         fpol(jpiglo/2+1:jpiglo) = 0._wp
221      ENDIF
222
223      ! 2. Ocean/land mask at u-,  v-, and z-points (computed from tmask)
224      ! -------------------------------------------
225!$OMP PARALLEL
226!$OMP DO schedule(static) private(jk, jj, ji)
227      DO jk = 1, jpk
228         DO jj = 1, jpjm1
229            DO ji = 1, fs_jpim1   ! vector loop
230               umask(ji,jj,jk) = tmask(ji,jj  ,jk) * tmask(ji+1,jj  ,jk)
231               vmask(ji,jj,jk) = tmask(ji,jj  ,jk) * tmask(ji  ,jj+1,jk)
232            END DO
233            DO ji = 1, jpim1      ! NO vector opt.
234               fmask(ji,jj,jk) = tmask(ji,jj  ,jk) * tmask(ji+1,jj  ,jk)   &
235                  &            * tmask(ji,jj+1,jk) * tmask(ji+1,jj+1,jk)
236            END DO
237         END DO
238      END DO
239      ! (ISF) MIN(1,SUM(umask)) is here to check if you have effectively at least 1 wet cell at u point
240!$OMP DO schedule(static) private(jj, ji)
241      DO jj = 1, jpjm1
242         DO ji = 1, fs_jpim1   ! vector loop
243            ssumask(ji,jj)  = ssmask(ji,jj) * ssmask(ji+1,jj  )  * MIN(1._wp,SUM(umask(ji,jj,:)))
244            ssvmask(ji,jj)  = ssmask(ji,jj) * ssmask(ji  ,jj+1)  * MIN(1._wp,SUM(vmask(ji,jj,:)))
245         END DO
246         DO ji = 1, jpim1      ! NO vector opt.
247            ssfmask(ji,jj) =  ssmask(ji,jj  ) * ssmask(ji+1,jj  )   &
248               &            * ssmask(ji,jj+1) * ssmask(ji+1,jj+1) * MIN(1._wp,SUM(fmask(ji,jj,:)))
249         END DO
250      END DO
251!$OMP END DO NOWAIT
252!$OMP END PARALLEL
253      CALL lbc_lnk( umask  , 'U', 1._wp )      ! Lateral boundary conditions
254      CALL lbc_lnk( vmask  , 'V', 1._wp )
255      CALL lbc_lnk( fmask  , 'F', 1._wp )
256      CALL lbc_lnk( ssumask, 'U', 1._wp )      ! Lateral boundary conditions
257      CALL lbc_lnk( ssvmask, 'V', 1._wp )
258      CALL lbc_lnk( ssfmask, 'F', 1._wp )
259
260      ! 3. Ocean/land mask at wu-, wv- and w points
261      !----------------------------------------------
262!$OMP PARALLEL
263!$OMP DO schedule(static) private(jj, ji)
264      DO jj = 1, jpj
265         DO ji = 1, jpi
266            wmask (ji,jj,1) = tmask(ji,jj,1)     ! surface
267            wumask(ji,jj,1) = umask(ji,jj,1)
268            wvmask(ji,jj,1) = vmask(ji,jj,1)
269         END DO 
270      END DO 
271!$OMP DO schedule(static) private(jk)
272      DO jk = 2, jpk                   ! interior values
273         wmask (:,:,jk) = tmask(:,:,jk) * tmask(:,:,jk-1)
274         wumask(:,:,jk) = umask(:,:,jk) * umask(:,:,jk-1)   
275         wvmask(:,:,jk) = vmask(:,:,jk) * vmask(:,:,jk-1)
276      END DO
277!$OMP END DO NOWAIT
278!$OMP END PARALLEL
279
280      ! Lateral boundary conditions on velocity (modify fmask)
281      ! ---------------------------------------     
282      DO jk = 1, jpk
283         zwf(:,:) = fmask(:,:,jk)         
284         DO jj = 2, jpjm1
285            DO ji = fs_2, fs_jpim1   ! vector opt.
286               IF( fmask(ji,jj,jk) == 0._wp ) THEN
287                  fmask(ji,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,jj), zwf(ji,jj+1),   &
288                     &                                           zwf(ji-1,jj), zwf(ji,jj-1)  )  )
289               ENDIF
290            END DO
291         END DO
292         DO jj = 2, jpjm1
293            IF( fmask(1,jj,jk) == 0._wp ) THEN
294               fmask(1  ,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1) ) )
295            ENDIF
296            IF( fmask(jpi,jj,jk) == 0._wp ) THEN
297               fmask(jpi,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) )
298            ENDIF
299         END DO         
300         DO ji = 2, jpim1
301            IF( fmask(ji,1,jk) == 0._wp ) THEN
302               fmask(ji, 1 ,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1) ) )
303            ENDIF
304            IF( fmask(ji,jpj,jk) == 0._wp ) THEN
305               fmask(ji,jpj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpjm1) ) )
306            ENDIF
307         END DO
308      END DO
309      !
310      IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN   ! ORCA_R2 configuration
311         !                                                 ! Increased lateral friction near of some straits
312         !                                ! Gibraltar strait  : partial slip (fmask=0.5)
313         ij0 = 101   ;   ij1 = 101
314         ii0 = 139   ;   ii1 = 140   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  0.5_wp
315         ij0 = 102   ;   ij1 = 102
316         ii0 = 139   ;   ii1 = 140   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  0.5_wp
317         !
318         !                                ! Bab el Mandeb : partial slip (fmask=1)
319         ij0 =  87   ;   ij1 =  88
320         ii0 = 160   ;   ii1 = 160   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  1._wp
321         ij0 =  88   ;   ij1 =  88
322         ii0 = 159   ;   ii1 = 159   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  1._wp
323         !
324         !                                ! Danish straits  : strong slip (fmask > 2)
325! We keep this as an example but it is instable in this case
326!         ij0 = 115   ;   ij1 = 115
327!         ii0 = 145   ;   ii1 = 146   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 4._wp
328!         ij0 = 116   ;   ij1 = 116
329!         ii0 = 145   ;   ii1 = 146   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 4._wp
330         !
331      ENDIF
332      !
333      IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN   ! ORCA R1 configuration
334         !                                                 ! Increased lateral friction near of some straits
335         ! This dirty section will be suppressed by simplification process:
336         ! all this will come back in input files
337         ! Currently these hard-wired indices relate to configuration with
338         ! extend grid (jpjglo=332)
339         !
340         isrow = 332 - jpjglo
341         !
342         IF(lwp) WRITE(numout,*)
343         IF(lwp) WRITE(numout,*) '   orca_r1: increase friction near the following straits : '
344         IF(lwp) WRITE(numout,*) '      Gibraltar '
345         ii0 = 282           ;   ii1 = 283        ! Gibraltar Strait
346         ij0 = 241 - isrow   ;   ij1 = 241 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp 
347
348         IF(lwp) WRITE(numout,*) '      Bhosporus '
349         ii0 = 314           ;   ii1 = 315        ! Bhosporus Strait
350         ij0 = 248 - isrow   ;   ij1 = 248 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp 
351
352         IF(lwp) WRITE(numout,*) '      Makassar (Top) '
353         ii0 =  48           ;   ii1 =  48        ! Makassar Strait (Top)
354         ij0 = 189 - isrow   ;   ij1 = 190 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp 
355
356         IF(lwp) WRITE(numout,*) '      Lombok '
357         ii0 =  44           ;   ii1 =  44        ! Lombok Strait
358         ij0 = 164 - isrow   ;   ij1 = 165 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp 
359
360         IF(lwp) WRITE(numout,*) '      Ombai '
361         ii0 =  53           ;   ii1 =  53        ! Ombai Strait
362         ij0 = 164 - isrow   ;   ij1 = 165 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp 
363
364         IF(lwp) WRITE(numout,*) '      Timor Passage '
365         ii0 =  56           ;   ii1 =  56        ! Timor Passage
366         ij0 = 164 - isrow   ;   ij1 = 165 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp 
367
368         IF(lwp) WRITE(numout,*) '      West Halmahera '
369         ii0 =  58           ;   ii1 =  58        ! West Halmahera Strait
370         ij0 = 181 - isrow   ;   ij1 = 182 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp 
371
372         IF(lwp) WRITE(numout,*) '      East Halmahera '
373         ii0 =  55           ;   ii1 =  55        ! East Halmahera Strait
374         ij0 = 181 - isrow   ;   ij1 = 182 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp 
375         !
376      ENDIF
377      !
378      CALL lbc_lnk( fmask, 'F', 1._wp )      ! Lateral boundary conditions on fmask
379      !
380      ! CAUTION : The fmask may be further modified in dyn_vor_init ( dynvor.F90 )
381      !
382      CALL wrk_dealloc( jpi, jpj, imsk )
383      CALL wrk_dealloc( jpi, jpj, zwf  )
384      !
385      IF( nn_timing == 1 )  CALL timing_stop('dom_msk')
386      !
387   END SUBROUTINE dom_msk
388   
389   !!======================================================================
390END MODULE dommsk
Note: See TracBrowser for help on using the repository browser.