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 @ 7037

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

ORCA2_LIM_PISCES hybrid version update

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