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

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

hybrid version update

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