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

source: branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90 @ 6069

Last change on this file since 6069 was 6069, checked in by timgraham, 8 years ago

Merge of dev_MetOffice_merge_2015 into branch (only NEMO directory for now).

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