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

Last change on this file since 6140 was 6140, checked in by timgraham, 9 years ago

Merge of branches/2015/dev_merge_2015 back into trunk. Merge excludes NEMOGCM/TOOLS/OBSTOOLS/ for now due to issues with the change of file type. Will sort these manually with further commits.

Branch merged as follows:
In the working copy of branch ran:
svn merge ​svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk@HEAD
Small conflicts due to bug fixes applied to trunk since the dev_merge_2015 was copied. Bug fixes were applied to the branch as well so these were easy to resolve.
Branch committed at this stage

In working copy run:
svn switch ​svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk
to switch working copy

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