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/UKMO/dev_r8864_restart_date/NEMOGCM/TOOLS/DOMAINcfg/src – NEMO

source: branches/UKMO/dev_r8864_restart_date/NEMOGCM/TOOLS/DOMAINcfg/src/dommsk.f90 @ 9235

Last change on this file since 9235 was 9235, checked in by davestorkey, 6 years ago

UKMO/dev_r8864_restart_date : clear SVN keywords.
File size: 18.1 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   !!----------------------------------------------------------------------
45   !!                   ***  vectopt_loop_substitute  ***
46   !!----------------------------------------------------------------------
47   !! ** purpose :   substitute the inner loop start/end indices with CPP macro
48   !!                allow unrolling of do-loop (useful with vector processors)
49   !!----------------------------------------------------------------------
50   !!----------------------------------------------------------------------
51   !! NEMO/OPA 3.7 , NEMO Consortium (2014)
52   !! $Id$
53   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
54   !!----------------------------------------------------------------------
55   !!----------------------------------------------------------------------
56   !! NEMO/OPA 3.2 , LODYC-IPSL  (2009)
57   !! $Id$
58   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
59   !!----------------------------------------------------------------------
60CONTAINS
61
62   SUBROUTINE dom_msk
63      !!---------------------------------------------------------------------
64      !!                 ***  ROUTINE dom_msk  ***
65      !!
66      !! ** Purpose :   Compute land/ocean mask arrays at tracer points, hori-
67      !!      zontal velocity points (u & v), vorticity points (f) points.
68      !!
69      !! ** Method  :   The ocean/land mask is computed from the basin bathy-
70      !!      metry in level (mbathy) which is defined or read in dommba.
71      !!      mbathy equals 0 over continental T-point
72      !!      and the number of ocean level over the ocean.
73      !!
74      !!      At a given position (ji,jj,jk) the ocean/land mask is given by:
75      !!      t-point : 0. IF mbathy( ji ,jj) =< 0
76      !!                1. IF mbathy( ji ,jj) >= jk
77      !!      u-point : 0. IF mbathy( ji ,jj)  or mbathy(ji+1, jj ) =< 0
78      !!                1. IF mbathy( ji ,jj) and mbathy(ji+1, jj ) >= jk.
79      !!      v-point : 0. IF mbathy( ji ,jj)  or mbathy( ji ,jj+1) =< 0
80      !!                1. IF mbathy( ji ,jj) and mbathy( ji ,jj+1) >= jk.
81      !!      f-point : 0. IF mbathy( ji ,jj)  or mbathy( ji ,jj+1)
82      !!                   or mbathy(ji+1,jj)  or mbathy(ji+1,jj+1) =< 0
83      !!                1. IF mbathy( ji ,jj) and mbathy( ji ,jj+1)
84      !!                  and mbathy(ji+1,jj) and mbathy(ji+1,jj+1) >= jk.
85      !!      tmask_i : interior ocean mask at t-point, i.e. excluding duplicated
86      !!                rows/lines due to cyclic or North Fold boundaries as well
87      !!                as MPP halos.
88      !!
89      !!        The lateral friction is set through the value of fmask along
90      !!      the coast and topography. This value is defined by rn_shlat, a
91      !!      namelist parameter:
92      !!         rn_shlat = 0, free slip  (no shear along the coast)
93      !!         rn_shlat = 2, no slip  (specified zero velocity at the coast)
94      !!         0 < rn_shlat < 2, partial slip   | non-linear velocity profile
95      !!         2 < rn_shlat, strong slip        | in the lateral boundary layer
96      !!
97      !!      N.B. If nperio not equal to 0, the land/ocean mask arrays
98      !!      are defined with the proper value at lateral domain boundaries.
99      !!
100      !!      In case of open boundaries (lk_bdy=T):
101      !!        - tmask is set to 1 on the points to be computed bay the open
102      !!          boundaries routines.
103      !!
104      !! ** Action :   tmask    : land/ocean mask at t-point (=0. or 1.)
105      !!               umask    : land/ocean mask at u-point (=0. or 1.)
106      !!               vmask    : land/ocean mask at v-point (=0. or 1.)
107      !!               fmask    : land/ocean mask at f-point (=0. or 1.)
108      !!                          =rn_shlat along lateral boundaries
109      !!               tmask_i  : interior ocean mask
110      !!----------------------------------------------------------------------
111      INTEGER  ::   ji, jj, jk               ! dummy loop indices
112      INTEGER  ::   iif, iil, ii0, ii1, ii   ! local integers
113      INTEGER  ::   ijf, ijl, ij0, ij1       !   -       -
114      INTEGER  ::   ios
115      INTEGER  ::   isrow                    ! index for ORCA1 starting row
116      INTEGER , POINTER, DIMENSION(:,:) ::  imsk
117      REAL(wp), POINTER, DIMENSION(:,:) ::  zwf
118      !!
119      NAMELIST/namlbc/ rn_shlat, ln_vorlat
120      !!---------------------------------------------------------------------
121      !
122      IF( nn_timing == 1 )  CALL timing_start('dom_msk')
123      !
124      CALL wrk_alloc( jpi, jpj, imsk )
125      CALL wrk_alloc( jpi, jpj, zwf  )
126      !
127      REWIND( numnam_ref )              ! Namelist namlbc in reference namelist : Lateral momentum boundary condition
128      READ  ( numnam_ref, namlbc, IOSTAT = ios, ERR = 901 )
129901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namlbc in reference namelist', lwp )
130
131      REWIND( numnam_cfg )              ! Namelist namlbc in configuration namelist : Lateral momentum boundary condition
132      READ  ( numnam_cfg, namlbc, IOSTAT = ios, ERR = 902 )
133902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namlbc in configuration namelist', lwp )
134      IF(lwm) WRITE ( numond, namlbc )
135     
136      IF(lwp) THEN                  ! control print
137         WRITE(numout,*)
138         WRITE(numout,*) 'dommsk : ocean mask '
139         WRITE(numout,*) '~~~~~~'
140         WRITE(numout,*) '   Namelist namlbc'
141         WRITE(numout,*) '      lateral momentum boundary cond.    rn_shlat  = ',rn_shlat
142         WRITE(numout,*) '      consistency with analytical form   ln_vorlat = ',ln_vorlat 
143      ENDIF
144
145      IF     (      rn_shlat == 0.               ) THEN   ;   IF(lwp) WRITE(numout,*) '   ocean lateral  free-slip '
146      ELSEIF (      rn_shlat == 2.               ) THEN   ;   IF(lwp) WRITE(numout,*) '   ocean lateral  no-slip '
147      ELSEIF ( 0. < rn_shlat .AND. rn_shlat < 2. ) THEN   ;   IF(lwp) WRITE(numout,*) '   ocean lateral  partial-slip '
148      ELSEIF ( 2. < rn_shlat                     ) THEN   ;   IF(lwp) WRITE(numout,*) '   ocean lateral  strong-slip '
149      ELSE
150         WRITE(ctmp1,*) ' rn_shlat is negative = ', rn_shlat
151         CALL ctl_stop( ctmp1 )
152      ENDIF
153
154      ! 1. Ocean/land mask at t-point (computed from mbathy)
155      ! -----------------------------
156      ! N.B. tmask has already the right boundary conditions since mbathy is ok
157      !
158      tmask(:,:,:) = 0._wp
159      DO jk = 1, jpk
160         DO jj = 1, jpj
161            DO ji = 1, jpi
162               IF( REAL( mbathy(ji,jj) - jk, wp ) + 0.1_wp >= 0._wp )   tmask(ji,jj,jk) = 1._wp
163            END DO 
164         END DO 
165      END DO 
166     
167      ! (ISF) define barotropic mask and mask the ice shelf point
168      ssmask(:,:)=tmask(:,:,1) ! at this stage ice shelf is not masked
169     
170      DO jk = 1, jpk
171         DO jj = 1, jpj
172            DO ji = 1, jpi
173               IF( REAL( misfdep(ji,jj) - jk, wp ) - 0.1_wp >= 0._wp )   THEN
174                  tmask(ji,jj,jk) = 0._wp
175               END IF
176            END DO 
177         END DO 
178      END DO 
179
180      ! Interior domain mask (used for global sum)
181      ! --------------------
182      tmask_i(:,:) = ssmask(:,:)            ! (ISH) tmask_i = 1 even on the ice shelf
183
184      tmask_h(:,:) = 1._wp                 ! 0 on the halo and 1 elsewhere
185      iif = jpreci                         ! ???
186      iil = nlci - jpreci + 1
187      ijf = jprecj                         ! ???
188      ijl = nlcj - jprecj + 1
189
190      tmask_h( 1 :iif,   :   ) = 0._wp      ! first columns
191      tmask_h(iil:jpi,   :   ) = 0._wp      ! last  columns (including mpp extra columns)
192      tmask_h(   :   , 1 :ijf) = 0._wp      ! first rows
193      tmask_h(   :   ,ijl:jpj) = 0._wp      ! last  rows (including mpp extra rows)
194
195      ! north fold mask
196      ! ---------------
197      tpol(1:jpiglo) = 1._wp 
198      fpol(1:jpiglo) = 1._wp
199      IF( jperio == 3 .OR. jperio == 4 ) THEN      ! T-point pivot
200         tpol(jpiglo/2+1:jpiglo) = 0._wp
201         fpol(     1    :jpiglo) = 0._wp
202         IF( mjg(nlej) == jpjglo ) THEN                  ! only half of the nlcj-1 row
203            DO ji = iif+1, iil-1
204               tmask_h(ji,nlej-1) = tmask_h(ji,nlej-1) * tpol(mig(ji))
205            END DO
206         ENDIF
207      ENDIF
208     
209      tmask_i(:,:) = tmask_i(:,:) * tmask_h(:,:)
210
211      IF( jperio == 5 .OR. jperio == 6 ) THEN      ! F-point pivot
212         tpol(     1    :jpiglo) = 0._wp
213         fpol(jpiglo/2+1:jpiglo) = 0._wp
214      ENDIF
215
216      ! 2. Ocean/land mask at u-,  v-, and z-points (computed from tmask)
217      ! -------------------------------------------
218      DO jk = 1, jpk
219         DO jj = 1, jpjm1
220            DO ji = 1, jpim1   ! vector loop
221               umask(ji,jj,jk) = tmask(ji,jj  ,jk) * tmask(ji+1,jj  ,jk)
222               vmask(ji,jj,jk) = tmask(ji,jj  ,jk) * tmask(ji  ,jj+1,jk)
223            END DO
224            DO ji = 1, jpim1      ! NO vector opt.
225               fmask(ji,jj,jk) = tmask(ji,jj  ,jk) * tmask(ji+1,jj  ,jk)   &
226                  &            * tmask(ji,jj+1,jk) * tmask(ji+1,jj+1,jk)
227            END DO
228         END DO
229      END DO
230      ! (ISF) MIN(1,SUM(umask)) is here to check if you have effectively at least 1 wet cell at u point
231      DO jj = 1, jpjm1
232         DO ji = 1, 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      CALL lbc_lnk( umask  , 'U', 1._wp )      ! Lateral boundary conditions
242      CALL lbc_lnk( vmask  , 'V', 1._wp )
243      CALL lbc_lnk( fmask  , 'F', 1._wp )
244      CALL lbc_lnk( ssumask, 'U', 1._wp )      ! Lateral boundary conditions
245      CALL lbc_lnk( ssvmask, 'V', 1._wp )
246      CALL lbc_lnk( ssfmask, 'F', 1._wp )
247
248      ! 3. Ocean/land mask at wu-, wv- and w points
249      !----------------------------------------------
250      wmask (:,:,1) = tmask(:,:,1)     ! surface
251      wumask(:,:,1) = umask(:,:,1)
252      wvmask(:,:,1) = vmask(:,:,1)
253      DO jk = 2, jpk                   ! interior values
254         wmask (:,:,jk) = tmask(:,:,jk) * tmask(:,:,jk-1)
255         wumask(:,:,jk) = umask(:,:,jk) * umask(:,:,jk-1)   
256         wvmask(:,:,jk) = vmask(:,:,jk) * vmask(:,:,jk-1)
257      END DO
258
259      ! Lateral boundary conditions on velocity (modify fmask)
260      ! ---------------------------------------     
261      DO jk = 1, jpk
262         zwf(:,:) = fmask(:,:,jk)         
263         DO jj = 2, jpjm1
264            DO ji = 2, jpim1   ! vector opt.
265               IF( fmask(ji,jj,jk) == 0._wp ) THEN
266                  fmask(ji,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,jj), zwf(ji,jj+1),   &
267                     &                                           zwf(ji-1,jj), zwf(ji,jj-1)  )  )
268               ENDIF
269            END DO
270         END DO
271         DO jj = 2, jpjm1
272            IF( fmask(1,jj,jk) == 0._wp ) THEN
273               fmask(1  ,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1) ) )
274            ENDIF
275            IF( fmask(jpi,jj,jk) == 0._wp ) THEN
276               fmask(jpi,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) )
277            ENDIF
278         END DO         
279         DO ji = 2, jpim1
280            IF( fmask(ji,1,jk) == 0._wp ) THEN
281               fmask(ji, 1 ,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1) ) )
282            ENDIF
283            IF( fmask(ji,jpj,jk) == 0._wp ) THEN
284               fmask(ji,jpj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpjm1) ) )
285            ENDIF
286         END DO
287      END DO
288      !
289      IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN   ! ORCA_R2 configuration
290         !                                                 ! Increased lateral friction near of some straits
291         !                                ! Gibraltar strait  : partial slip (fmask=0.5)
292         ij0 = 101   ;   ij1 = 101
293         ii0 = 139   ;   ii1 = 140   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  0.5_wp
294         ij0 = 102   ;   ij1 = 102
295         ii0 = 139   ;   ii1 = 140   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  0.5_wp
296         !
297         !                                ! Bab el Mandeb : partial slip (fmask=1)
298         ij0 =  87   ;   ij1 =  88
299         ii0 = 160   ;   ii1 = 160   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  1._wp
300         ij0 =  88   ;   ij1 =  88
301         ii0 = 159   ;   ii1 = 159   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  1._wp
302         !
303         !                                ! Danish straits  : strong slip (fmask > 2)
304! We keep this as an example but it is instable in this case
305!         ij0 = 115   ;   ij1 = 115
306!         ii0 = 145   ;   ii1 = 146   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 4._wp
307!         ij0 = 116   ;   ij1 = 116
308!         ii0 = 145   ;   ii1 = 146   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 4._wp
309         !
310      ENDIF
311      !
312      IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN   ! ORCA R1 configuration
313         !                                                 ! Increased lateral friction near of some straits
314         ! This dirty section will be suppressed by simplification process:
315         ! all this will come back in input files
316         ! Currently these hard-wired indices relate to configuration with
317         ! extend grid (jpjglo=332)
318         !
319         isrow = 332 - jpjglo
320         !
321         IF(lwp) WRITE(numout,*)
322         IF(lwp) WRITE(numout,*) '   orca_r1: increase friction near the following straits : '
323         IF(lwp) WRITE(numout,*) '      Gibraltar '
324         ii0 = 282           ;   ii1 = 283        ! Gibraltar Strait
325         ij0 = 241 - isrow   ;   ij1 = 241 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp 
326
327         IF(lwp) WRITE(numout,*) '      Bhosporus '
328         ii0 = 314           ;   ii1 = 315        ! Bhosporus Strait
329         ij0 = 248 - isrow   ;   ij1 = 248 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp 
330
331         IF(lwp) WRITE(numout,*) '      Makassar (Top) '
332         ii0 =  48           ;   ii1 =  48        ! Makassar Strait (Top)
333         ij0 = 189 - isrow   ;   ij1 = 190 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp 
334
335         IF(lwp) WRITE(numout,*) '      Lombok '
336         ii0 =  44           ;   ii1 =  44        ! Lombok Strait
337         ij0 = 164 - isrow   ;   ij1 = 165 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp 
338
339         IF(lwp) WRITE(numout,*) '      Ombai '
340         ii0 =  53           ;   ii1 =  53        ! Ombai Strait
341         ij0 = 164 - isrow   ;   ij1 = 165 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp 
342
343         IF(lwp) WRITE(numout,*) '      Timor Passage '
344         ii0 =  56           ;   ii1 =  56        ! Timor Passage
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,*) '      West Halmahera '
348         ii0 =  58           ;   ii1 =  58        ! West Halmahera Strait
349         ij0 = 181 - isrow   ;   ij1 = 182 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp 
350
351         IF(lwp) WRITE(numout,*) '      East Halmahera '
352         ii0 =  55           ;   ii1 =  55        ! East Halmahera Strait
353         ij0 = 181 - isrow   ;   ij1 = 182 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp 
354         !
355      ENDIF
356      !
357      CALL lbc_lnk( fmask, 'F', 1._wp )      ! Lateral boundary conditions on fmask
358      !
359      ! CAUTION : The fmask may be further modified in dyn_vor_init ( dynvor.F90 )
360      !
361      CALL wrk_dealloc( jpi, jpj, imsk )
362      CALL wrk_dealloc( jpi, jpj, zwf  )
363      !
364      IF( nn_timing == 1 )  CALL timing_stop('dom_msk')
365      !
366   END SUBROUTINE dom_msk
367   
368   !!======================================================================
369END MODULE dommsk
Note: See TracBrowser for help on using the repository browser.