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

source: branches/UKMO/dev_r5518_GO6_diag_bitcomp/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90 @ 9707

Last change on this file since 9707 was 9707, checked in by frrh, 7 years ago

Work out N-fold mask points for U and V grids based on
value flipping during N-fold lbc_lnk operations.
It's much easier this way.
Tidy up remaining code.

File size: 34.3 KB
Line 
1
2MODULE dommsk
3   !!======================================================================
4   !!                       ***  MODULE dommsk   ***
5   !! Ocean initialization : domain land/sea mask
6   !!======================================================================
7   !! History :  OPA  ! 1987-07  (G. Madec)  Original code
8   !!            6.0  ! 1993-03  (M. Guyon)  symetrical conditions (M. Guyon)
9   !!            7.0  ! 1996-01  (G. Madec)  suppression of common work arrays
10   !!             -   ! 1996-05  (G. Madec)  mask computed from tmask and sup-
11   !!                 !                      pression of the double computation of bmask
12   !!            8.0  ! 1997-02  (G. Madec)  mesh information put in domhgr.F
13   !!            8.1  ! 1997-07  (G. Madec)  modification of mbathy and fmask
14   !!             -   ! 1998-05  (G. Roullet)  free surface
15   !!            8.2  ! 2000-03  (G. Madec)  no slip accurate
16   !!             -   ! 2001-09  (J.-M. Molines)  Open boundaries
17   !!   NEMO     1.0  ! 2002-08  (G. Madec)  F90: Free form and module
18   !!             -   ! 2005-11  (V. Garnier) Surface pressure gradient organization
19   !!            3.2  ! 2009-07  (R. Benshila) Suppression of rigid-lid option
20   !!----------------------------------------------------------------------
21
22   !!----------------------------------------------------------------------
23   !!   dom_msk        : compute land/ocean mask
24   !!   dom_msk_nsa    : update land/ocean mask when no-slip accurate option is used.
25   !!----------------------------------------------------------------------
26   USE oce             ! ocean dynamics and tracers
27   USE dom_oce         ! ocean space and time domain
28   USE in_out_manager  ! I/O manager
29   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
30   USE lib_mpp
31   USE dynspg_oce      ! choice/control of key cpp for surface pressure gradient
32   USE wrk_nemo        ! Memory allocation
33   USE timing          ! Timing
34
35   IMPLICIT NONE
36   PRIVATE
37
38   PUBLIC   dom_msk         ! routine called by inidom.F90
39   PUBLIC   dom_msk_alloc   ! routine called by nemogcm.F90
40
41   !                            !!* Namelist namlbc : lateral boundary condition *
42   REAL(wp)        :: rn_shlat   ! type of lateral boundary condition on velocity
43   LOGICAL, PUBLIC :: ln_vorlat  !  consistency of vorticity boundary condition
44   !                                            with analytical eqs.
45
46
47   INTEGER, ALLOCATABLE, SAVE, DIMENSION(:,:) ::  icoord ! Workspace for dom_msk_nsa()
48
49   !! * Substitutions
50#  include "vectopt_loop_substitute.h90"
51   !!----------------------------------------------------------------------
52   !! NEMO/OPA 3.2 , LODYC-IPSL  (2009)
53   !! $Id$
54   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
55   !!----------------------------------------------------------------------
56CONTAINS
57   
58   INTEGER FUNCTION dom_msk_alloc()
59      !!---------------------------------------------------------------------
60      !!                 ***  FUNCTION dom_msk_alloc  ***
61      !!---------------------------------------------------------------------
62      dom_msk_alloc = 0
63#if defined key_noslip_accurate
64      ALLOCATE(icoord(jpi*jpj*jpk,3), STAT=dom_msk_alloc)
65#endif
66      IF( dom_msk_alloc /= 0 )   CALL ctl_warn('dom_msk_alloc: failed to allocate icoord array')
67      !
68   END FUNCTION dom_msk_alloc
69
70
71   SUBROUTINE dom_msk
72      !!---------------------------------------------------------------------
73      !!                 ***  ROUTINE dom_msk  ***
74      !!
75      !! ** Purpose :   Compute land/ocean mask arrays at tracer points, hori-
76      !!      zontal velocity points (u & v), vorticity points (f) and baro-
77      !!      tropic stream function  points (b).
78      !!
79      !! ** Method  :   The ocean/land mask is computed from the basin bathy-
80      !!      metry in level (mbathy) which is defined or read in dommba.
81      !!      mbathy equals 0 over continental T-point
82      !!      and the number of ocean level over the ocean.
83      !!
84      !!      At a given position (ji,jj,jk) the ocean/land mask is given by:
85      !!      t-point : 0. IF mbathy( ji ,jj) =< 0
86      !!                1. IF mbathy( ji ,jj) >= jk
87      !!      u-point : 0. IF mbathy( ji ,jj)  or mbathy(ji+1, jj ) =< 0
88      !!                1. IF mbathy( ji ,jj) and mbathy(ji+1, jj ) >= jk.
89      !!      v-point : 0. IF mbathy( ji ,jj)  or mbathy( ji ,jj+1) =< 0
90      !!                1. IF mbathy( ji ,jj) and mbathy( ji ,jj+1) >= jk.
91      !!      f-point : 0. IF mbathy( ji ,jj)  or mbathy( ji ,jj+1)
92      !!                   or mbathy(ji+1,jj)  or mbathy(ji+1,jj+1) =< 0
93      !!                1. IF mbathy( ji ,jj) and mbathy( ji ,jj+1)
94      !!                  and mbathy(ji+1,jj) and mbathy(ji+1,jj+1) >= jk.
95      !!      b-point : the same definition as for f-point of the first ocean
96      !!                level (surface level) but with 0 along coastlines.
97      !!      tmask_i : interior ocean mask at t-point, i.e. excluding duplicated
98      !!                rows/lines due to cyclic or North Fold boundaries as well
99      !!                as MPP halos.
100      !!
101      !!        The lateral friction is set through the value of fmask along
102      !!      the coast and topography. This value is defined by rn_shlat, a
103      !!      namelist parameter:
104      !!         rn_shlat = 0, free slip  (no shear along the coast)
105      !!         rn_shlat = 2, no slip  (specified zero velocity at the coast)
106      !!         0 < rn_shlat < 2, partial slip   | non-linear velocity profile
107      !!         2 < rn_shlat, strong slip        | in the lateral boundary layer
108      !!
109      !!      N.B. If nperio not equal to 0, the land/ocean mask arrays
110      !!      are defined with the proper value at lateral domain boundaries,
111      !!      but bmask. indeed, bmask defined the domain over which the
112      !!      barotropic stream function is computed. this domain cannot
113      !!      contain identical columns because the matrix associated with
114      !!      the barotropic stream function equation is then no more inverti-
115      !!      ble. therefore bmask is set to 0 along lateral domain boundaries
116      !!      even IF nperio is not zero.
117      !!
118      !!      In case of open boundaries (lk_bdy=T):
119      !!        - tmask is set to 1 on the points to be computed bay the open
120      !!          boundaries routines.
121      !!        - bmask is  set to 0 on the open boundaries.
122      !!
123      !! ** Action :   tmask    : land/ocean mask at t-point (=0. or 1.)
124      !!               umask    : land/ocean mask at u-point (=0. or 1.)
125      !!               vmask    : land/ocean mask at v-point (=0. or 1.)
126      !!               fmask    : land/ocean mask at f-point (=0. or 1.)
127      !!                          =rn_shlat along lateral boundaries
128      !!               bmask    : land/ocean mask at barotropic stream
129      !!                          function point (=0. or 1.) and set to 0 along lateral boundaries
130      !!               tmask_i  : interior ocean mask
131      !!----------------------------------------------------------------------
132      !
133      INTEGER  ::   ji, jj, jk      ! dummy loop indices
134      INTEGER  ::   iif, iil, ii0, ii1, ii   ! local integers
135      INTEGER  ::   ijf, ijl, ij0, ij1       !   -       -
136      INTEGER  ::   ios
137      INTEGER  ::   isrow                    ! index for ORCA1 starting row
138      INTEGER , POINTER, DIMENSION(:,:) ::  imsk
139      REAL(wp) ::  zphi_drake_passage, zshlat_antarc
140      REAL(wp), POINTER, DIMENSION(:,:) ::  zwf
141      REAL(wp) :: uv(jpi,jpj)
142      !!
143      NAMELIST/namlbc/ rn_shlat, ln_vorlat
144      !!---------------------------------------------------------------------
145      !
146      IF( nn_timing == 1 )  CALL timing_start('dom_msk')
147      !
148      CALL wrk_alloc( jpi, jpj, imsk )
149      CALL wrk_alloc( jpi, jpj, zwf  )
150      !
151      REWIND( numnam_ref )              ! Namelist namlbc in reference namelist : Lateral momentum boundary condition
152      READ  ( numnam_ref, namlbc, IOSTAT = ios, ERR = 901 )
153901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namlbc in reference namelist', lwp )
154
155      REWIND( numnam_cfg )              ! Namelist namlbc in configuration namelist : Lateral momentum boundary condition
156      READ  ( numnam_cfg, namlbc, IOSTAT = ios, ERR = 902 )
157902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namlbc in configuration namelist', lwp )
158      IF(lwm) WRITE ( numond, namlbc )
159     
160      IF(lwp) THEN                  ! control print
161         WRITE(numout,*)
162         WRITE(numout,*) 'dommsk : ocean mask '
163         WRITE(numout,*) '~~~~~~'
164         WRITE(numout,*) '   Namelist namlbc'
165         WRITE(numout,*) '      lateral momentum boundary cond.    rn_shlat  = ',rn_shlat
166         WRITE(numout,*) '      consistency with analytical form   ln_vorlat = ',ln_vorlat 
167      ENDIF
168
169      IF     (      rn_shlat == 0.               ) THEN   ;   IF(lwp) WRITE(numout,*) '   ocean lateral  free-slip '
170      ELSEIF (      rn_shlat == 2.               ) THEN   ;   IF(lwp) WRITE(numout,*) '   ocean lateral  no-slip '
171      ELSEIF ( 0. < rn_shlat .AND. rn_shlat < 2. ) THEN   ;   IF(lwp) WRITE(numout,*) '   ocean lateral  partial-slip '
172      ELSEIF ( 2. < rn_shlat                     ) THEN   ;   IF(lwp) WRITE(numout,*) '   ocean lateral  strong-slip '
173      ELSE
174         WRITE(ctmp1,*) ' rn_shlat is negative = ', rn_shlat
175         CALL ctl_stop( ctmp1 )
176      ENDIF
177
178      ! 1. Ocean/land mask at t-point (computed from mbathy)
179      ! -----------------------------
180      ! N.B. tmask has already the right boundary conditions since mbathy is ok
181      !
182      tmask(:,:,:) = 0._wp
183      DO jk = 1, jpk
184         DO jj = 1, jpj
185            DO ji = 1, jpi
186               IF( REAL( mbathy(ji,jj) - jk, wp ) + 0.1_wp >= 0._wp )   tmask(ji,jj,jk) = 1._wp
187            END DO 
188         END DO 
189      END DO 
190     
191      ! (ISF) define barotropic mask and mask the ice shelf point
192      ssmask(:,:)=tmask(:,:,1) ! at this stage ice shelf is not masked
193     
194      DO jk = 1, jpk
195         DO jj = 1, jpj
196            DO ji = 1, jpi
197               IF( REAL( misfdep(ji,jj) - jk, wp ) - 0.1_wp >= 0._wp )   THEN
198                  tmask(ji,jj,jk) = 0._wp
199               END IF
200            END DO 
201         END DO 
202      END DO 
203
204!!gm  ????
205#if defined key_zdfkpp
206      IF( cp_cfg == 'orca' ) THEN
207         IF( jp_cfg == 2 )   THEN       ! land point on Bab el Mandeb zonal section
208            ij0 =  87   ;   ij1 =  88
209            ii0 = 160   ;   ii1 = 161
210            tmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 0._wp
211         ELSE
212            IF(lwp) WRITE(numout,*)
213            IF(lwp) WRITE(numout,cform_war)
214            IF(lwp) WRITE(numout,*)
215            IF(lwp) WRITE(numout,*)'          A mask must be applied on Bab el Mandeb strait'
216            IF(lwp) WRITE(numout,*)'          in case of ORCAs configurations'
217            IF(lwp) WRITE(numout,*)'          This is a problem which is not yet solved'
218            IF(lwp) WRITE(numout,*)
219         ENDIF
220      ENDIF
221#endif
222!!gm end
223
224      ! Interior domain mask (used for global sum)
225      ! --------------------
226      tmask_i(:,:) = ssmask(:,:)            ! (ISH) tmask_i = 1 even on the ice shelf
227
228      iif = jpreci                         ! ???
229      iil = nlci - jpreci + 1
230      ijf = jprecj                         ! ???
231      ijl = nlcj - jprecj + 1
232
233      tmask_i( 1 :iif,   :   ) = 0._wp      ! first columns
234      tmask_i(iil:jpi,   :   ) = 0._wp      ! last  columns (including mpp extra columns)
235      tmask_i(   :   , 1 :ijf) = 0._wp      ! first rows
236      tmask_i(   :   ,ijl:jpj) = 0._wp      ! last  rows (including mpp extra rows)
237
238      ! north fold mask
239      ! ---------------
240      tpol(1:jpiglo) = 1._wp 
241      fpol(1:jpiglo) = 1._wp
242      IF( jperio == 3 .OR. jperio == 4 ) THEN      ! T-point pivot
243         tpol(jpiglo/2+1:jpiglo) = 0._wp
244         fpol(     1    :jpiglo) = 0._wp
245         IF( mjg(nlej) == jpjglo ) THEN                  ! only half of the nlcj-1 row
246            DO ji = iif+1, iil-1
247               tmask_i(ji,nlej-1) = tmask_i(ji,nlej-1) * tpol(mig(ji))
248            END DO
249         ENDIF
250      ENDIF
251
252
253      IF( jperio == 5 .OR. jperio == 6 ) THEN      ! F-point pivot
254         tpol(     1    :jpiglo) = 0._wp
255         fpol(jpiglo/2+1:jpiglo) = 0._wp
256      ENDIF
257
258      ! 2. Ocean/land mask at u-,  v-, and z-points (computed from tmask)
259      ! -------------------------------------------
260      DO jk = 1, jpk
261         DO jj = 1, jpjm1
262            DO ji = 1, fs_jpim1   ! vector loop
263               umask(ji,jj,jk) = tmask(ji,jj  ,jk) * tmask(ji+1,jj  ,jk)
264               vmask(ji,jj,jk) = tmask(ji,jj  ,jk) * tmask(ji  ,jj+1,jk)
265            END DO
266            DO ji = 1, jpim1      ! NO vector opt.
267               fmask(ji,jj,jk) = tmask(ji,jj  ,jk) * tmask(ji+1,jj  ,jk)   &
268                  &            * tmask(ji,jj+1,jk) * tmask(ji+1,jj+1,jk)
269            END DO
270         END DO
271      END DO
272      ! (ISF) MIN(1,SUM(umask)) is here to check if you have effectively at least 1 wet u point
273      DO jj = 1, jpjm1
274         DO ji = 1, fs_jpim1   ! vector loop
275            umask_i(ji,jj)  = ssmask(ji,jj) * ssmask(ji+1,jj  )  * MIN(1._wp,SUM(umask(ji,jj,:)))
276            vmask_i(ji,jj)  = ssmask(ji,jj) * ssmask(ji  ,jj+1)  * MIN(1._wp,SUM(vmask(ji,jj,:)))
277         END DO
278         DO ji = 1, jpim1      ! NO vector opt.
279            fmask_i(ji,jj) =  ssmask(ji,jj  ) * ssmask(ji+1,jj  )   &
280               &            * ssmask(ji,jj+1) * ssmask(ji+1,jj+1) * MIN(1._wp,SUM(fmask(ji,jj,:)))
281         END DO
282      END DO
283      CALL lbc_lnk( umask, 'U', 1._wp )      ! Lateral boundary conditions
284      CALL lbc_lnk( vmask, 'V', 1._wp )
285      CALL lbc_lnk( fmask, 'F', 1._wp )
286      CALL lbc_lnk( umask_i, 'U', 1._wp )      ! Lateral boundary conditions
287      CALL lbc_lnk( vmask_i, 'V', 1._wp )
288      CALL lbc_lnk( fmask_i, 'F', 1._wp )
289
290
291      ! Set up mask for diagnostics on T points, to exclude duplicate
292      ! data points in wrap and N-fold regions.
293      DO jk = 1, jpk
294         tmask_i_diag(:,:,jk) = tmask(:,:,jk) * tmask_i(:,:)
295      END DO
296
297      ! Set up mask for diagnostics on U points, to exclude duplicate
298      ! data points in wrap and N-fold regions.
299      umask_i_diag(:,:,:) = 1.0
300      umask_i_diag(2:jpim1,:,:) = tmask(2:jpim1,:,:) + tmask(3:jpi,:,:)
301      CALL lbc_lnk( umask_i_diag, 'U', 1. )
302      ! Now mask out any wrap columns
303      umask_i_diag( 1 :iif,:,:) = 0._wp       ! first columns
304      umask_i_diag(iil:jpi,:,:) = 0._wp       ! last  columns (including mpp extra columns)
305      ! Now mask out any extra bottom rows
306      umask_i_diag(:,1:ijf,:) = 0._wp         ! first rows
307
308      ! Now find which points change sign during a U lbc_lnk to find
309      ! out which points to mask in the N fold.   
310      uv(:,:) = 1.0
311      CALL lbc_lnk( uv, 'U', -1. )
312
313      ! Now find out which points have changed sign and mask them
314      DO jj = 1, jpj
315         DO ji = 1, jpi
316            IF (uv(ji,jj) < 0.0) THEN
317               umask_i_diag(ji,jj,:) = 0.0
318            END IF
319         END DO
320      END DO
321
322      ! Set up mask for diagnostics on V points, to exclude duplicate
323      ! data points in wrap and N-fold regions.
324      vmask_i_diag(:,:,:) = 1.0
325      vmask_i_diag(:,2:jpjm1,:) = tmask(:,2:jpjm1,:) + tmask(:,3:jpj,:)
326      CALL lbc_lnk( vmask_i_diag, 'V', 1. )
327
328      ! Now mask out any wrap columns
329      vmask_i_diag( 1 :iif,:,:) = 0._wp       ! first columns
330      vmask_i_diag(iil:jpi,:,:) = 0._wp       ! last  columns (including mpp extra columns)
331      ! Now mask out any extra rows
332      vmask_i_diag(:,1:ijf,:) = 0._wp         ! first rows
333
334      ! Now find which points change sign during a V lbc_lnk to find
335      ! out which points to mask in the N fold.   
336      uv(:,:) = 1.0
337      CALL lbc_lnk( uv, 'V', -1. )
338
339      ! Now find out which points have changed sign and mask them
340      DO jj = 1, jpj
341         DO ji = 1, jpi
342            IF (uv(ji,jj) < 0.0) THEN
343               vmask_i_diag(ji,jj,:) = 0.0
344            END IF
345         END DO
346      END DO
347
348
349      ! 3. Ocean/land mask at wu-, wv- and w points
350      !----------------------------------------------
351      wmask (:,:,1) = tmask(:,:,1) ! ????????
352      wumask(:,:,1) = umask(:,:,1) ! ????????
353      wvmask(:,:,1) = vmask(:,:,1) ! ????????
354      DO jk=2,jpk
355         wmask (:,:,jk)=tmask(:,:,jk) * tmask(:,:,jk-1)
356         wumask(:,:,jk)=umask(:,:,jk) * umask(:,:,jk-1)   
357         wvmask(:,:,jk)=vmask(:,:,jk) * vmask(:,:,jk-1)
358      END DO
359
360      ! 4. ocean/land mask for the elliptic equation
361      ! --------------------------------------------
362      bmask(:,:) = ssmask(:,:)       ! elliptic equation is written at t-point
363      !
364      !                               ! Boundary conditions
365      !                                    ! cyclic east-west : bmask must be set to 0. on rows 1 and jpi
366      IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN
367         bmask( 1 ,:) = 0._wp
368         bmask(jpi,:) = 0._wp
369      ENDIF
370      IF( nperio == 2 ) THEN               ! south symmetric :  bmask must be set to 0. on row 1
371         bmask(:, 1 ) = 0._wp
372      ENDIF
373      !                                    ! north fold :
374      IF( nperio == 3 .OR. nperio == 4 ) THEN   ! T-pt pivot : bmask set to 0. on row jpj and on half jpjglo-1 row
375         DO ji = 1, jpi                     
376            ii = ji + nimpp - 1
377            bmask(ji,jpj-1) = bmask(ji,jpj-1) * tpol(ii)
378            bmask(ji,jpj  ) = 0._wp
379         END DO
380      ENDIF
381      IF( nperio == 5 .OR. nperio == 6 ) THEN   ! F-pt pivot and T-pt elliptic eq. : bmask set to 0. on row jpj
382         bmask(:,jpj) = 0._wp
383      ENDIF
384      !
385      IF( lk_mpp ) THEN                    ! mpp specificities
386         !                                      ! bmask is set to zero on the overlap region
387         IF( nbondi /= -1 .AND. nbondi /= 2 )   bmask(  1 :jpreci,:) = 0._wp
388         IF( nbondi /=  1 .AND. nbondi /= 2 )   bmask(nlci:jpi   ,:) = 0._wp
389         IF( nbondj /= -1 .AND. nbondj /= 2 )   bmask(:,  1 :jprecj) = 0._wp
390         IF( nbondj /=  1 .AND. nbondj /= 2 )   bmask(:,nlcj:jpj   ) = 0._wp
391         !
392         IF( npolj == 3 .OR. npolj == 4 ) THEN  ! north fold : bmask must be set to 0. on rows jpj-1 and jpj
393            DO ji = 1, nlci
394               ii = ji + nimpp - 1
395               bmask(ji,nlcj-1) = bmask(ji,nlcj-1) * tpol(ii)
396               bmask(ji,nlcj  ) = 0._wp
397            END DO
398         ENDIF
399         IF( npolj == 5 .OR. npolj == 6 ) THEN  ! F-pt pivot and T-pt elliptic eq. : bmask set to 0. on row jpj
400            DO ji = 1, nlci
401               bmask(ji,nlcj  ) = 0._wp
402            END DO
403         ENDIF
404      ENDIF
405
406
407      ! mask for second order calculation of vorticity
408      ! ----------------------------------------------
409      CALL dom_msk_nsa
410
411     
412      ! Lateral boundary conditions on velocity (modify fmask)
413      ! ---------------------------------------     
414      DO jk = 1, jpk
415         zwf(:,:) = fmask(:,:,jk)         
416         DO jj = 2, jpjm1
417            DO ji = fs_2, fs_jpim1   ! vector opt.
418               IF( fmask(ji,jj,jk) == 0._wp ) THEN
419                  fmask(ji,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,jj), zwf(ji,jj+1),   &
420                     &                                           zwf(ji-1,jj), zwf(ji,jj-1)  )  )
421               ENDIF
422            END DO
423         END DO
424         DO jj = 2, jpjm1
425            IF( fmask(1,jj,jk) == 0._wp ) THEN
426               fmask(1  ,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1) ) )
427            ENDIF
428            IF( fmask(jpi,jj,jk) == 0._wp ) THEN
429               fmask(jpi,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) )
430            ENDIF
431         END DO         
432         DO ji = 2, jpim1
433            IF( fmask(ji,1,jk) == 0._wp ) THEN
434               fmask(ji, 1 ,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1) ) )
435            ENDIF
436            IF( fmask(ji,jpj,jk) == 0._wp ) THEN
437               fmask(ji,jpj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpjm1) ) )
438            ENDIF
439         END DO
440      END DO
441      !
442      IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN   ! ORCA_R2 configuration
443         !                                                 ! Increased lateral friction near of some straits
444         IF( nn_cla == 0 ) THEN
445            !                                ! Gibraltar strait  : partial slip (fmask=0.5)
446            ij0 = 101   ;   ij1 = 101
447            ii0 = 139   ;   ii1 = 140   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  0.5_wp
448            ij0 = 102   ;   ij1 = 102
449            ii0 = 139   ;   ii1 = 140   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  0.5_wp
450            !
451            !                                ! Bab el Mandeb : partial slip (fmask=1)
452            ij0 =  87   ;   ij1 =  88
453            ii0 = 160   ;   ii1 = 160   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  1._wp
454            ij0 =  88   ;   ij1 =  88
455            ii0 = 159   ;   ii1 = 159   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  1._wp
456            !
457         ENDIF
458         !                                ! Danish straits  : strong slip (fmask > 2)
459! We keep this as an example but it is instable in this case
460!         ij0 = 115   ;   ij1 = 115
461!         ii0 = 145   ;   ii1 = 146   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 4._wp
462!         ij0 = 116   ;   ij1 = 116
463!         ii0 = 145   ;   ii1 = 146   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 4._wp
464         !
465      ENDIF
466      !
467      IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN   ! ORCA R1 configuration
468         !                                                 ! Increased lateral friction near of some straits
469         ! This dirty section will be suppressed by simplification process:
470         ! all this will come back in input files
471         ! Currently these hard-wired indices relate to configuration with
472         ! extend grid (jpjglo=332)
473         !
474         isrow = 332 - jpjglo
475         !
476         IF(lwp) WRITE(numout,*)
477         IF(lwp) WRITE(numout,*) '   orca_r1: increase friction near the following straits : '
478         IF(lwp) WRITE(numout,*) '      Gibraltar '
479         ii0 = 282           ;   ii1 = 283        ! Gibraltar Strait
480         ij0 = 241 - isrow   ;   ij1 = 241 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp 
481
482         IF(lwp) WRITE(numout,*) '      Bhosporus '
483         ii0 = 314           ;   ii1 = 315        ! Bhosporus Strait
484         ij0 = 248 - isrow   ;   ij1 = 248 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp 
485
486         IF(lwp) WRITE(numout,*) '      Makassar (Top) '
487         ii0 =  48           ;   ii1 =  48        ! Makassar Strait (Top)
488         ij0 = 189 - isrow   ;   ij1 = 190 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp 
489
490         IF(lwp) WRITE(numout,*) '      Lombok '
491         ii0 =  44           ;   ii1 =  44        ! Lombok Strait
492         ij0 = 164 - isrow   ;   ij1 = 165 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp 
493
494         IF(lwp) WRITE(numout,*) '      Ombai '
495         ii0 =  53           ;   ii1 =  53        ! Ombai Strait
496         ij0 = 164 - isrow   ;   ij1 = 165 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp 
497
498         IF(lwp) WRITE(numout,*) '      Timor Passage '
499         ii0 =  56           ;   ii1 =  56        ! Timor Passage
500         ij0 = 164 - isrow   ;   ij1 = 165 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp 
501
502         IF(lwp) WRITE(numout,*) '      West Halmahera '
503         ii0 =  58           ;   ii1 =  58        ! West Halmahera Strait
504         ij0 = 181 - isrow   ;   ij1 = 182 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp 
505
506         IF(lwp) WRITE(numout,*) '      East Halmahera '
507         ii0 =  55           ;   ii1 =  55        ! East Halmahera Strait
508         ij0 = 181 - isrow   ;   ij1 = 182 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp 
509         !
510      ENDIF
511      !
512      IF( cp_cfg == "orca" .AND. jp_cfg == 025 .AND. rn_shlat == 0.0 ) THEN   
513         !                                              ! ORCA_R025 configuration
514         !                                              ! Increased lateral friction on parts of Antarctic coastline
515         !                                              ! for increased stability
516         !                                              ! NB. This only works to do this here if we have free slip
517         !                                              ! generally, so fmask is zero at coast points.
518         IF(lwp) WRITE(numout,*)
519         IF(lwp) WRITE(numout,*) '   orca_r025: increase friction in following regions : '
520         IF(lwp) WRITE(numout,*) '      whole Antarctic coastline: partial slip shlat=1 '
521
522         zphi_drake_passage = -58.0_wp
523         zshlat_antarc = 1.0_wp
524         zwf(:,:) = fmask(:,:,1)         
525         DO jj = 2, jpjm1
526            DO ji = fs_2, fs_jpim1   ! vector opt.
527               IF( gphif(ji,jj) .lt. zphi_drake_passage .and. fmask(ji,jj,1) == 0._wp ) THEN
528                  fmask(ji,jj,:) = zshlat_antarc * MIN( 1._wp , MAX( zwf(ji+1,jj), zwf(ji,jj+1),   &
529                     &                                           zwf(ji-1,jj), zwf(ji,jj-1)  )  )
530               ENDIF
531            END DO
532         END DO
533      END IF
534      !
535      CALL lbc_lnk( fmask, 'F', 1._wp )      ! Lateral boundary conditions on fmask
536
537      ! CAUTION : The fmask may be further modified in dyn_vor_init ( dynvor.F90 )
538           
539      IF( nprint == 1 .AND. lwp ) THEN      ! Control print
540         imsk(:,:) = INT( tmask_i(:,:) )
541         WRITE(numout,*) ' tmask_i : '
542         CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1,   &
543               &                           1, jpj, 1, 1, numout)
544         WRITE (numout,*)
545         WRITE (numout,*) ' dommsk: tmask for each level'
546         WRITE (numout,*) ' ----------------------------'
547         DO jk = 1, jpk
548            imsk(:,:) = INT( tmask(:,:,jk) )
549
550            WRITE(numout,*)
551            WRITE(numout,*) ' level = ',jk
552            CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1,   &
553               &                              1, jpj, 1, 1, numout)
554         END DO
555         WRITE(numout,*)
556         WRITE(numout,*) ' dom_msk: vmask for each level'
557         WRITE(numout,*) ' -----------------------------'
558         DO jk = 1, jpk
559            imsk(:,:) = INT( vmask(:,:,jk) )
560            WRITE(numout,*)
561            WRITE(numout,*) ' level = ',jk
562            CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1,   &
563               &                              1, jpj, 1, 1, numout)
564         END DO
565         WRITE(numout,*)
566         WRITE(numout,*) ' dom_msk: fmask for each level'
567         WRITE(numout,*) ' -----------------------------'
568         DO jk = 1, jpk
569            imsk(:,:) = INT( fmask(:,:,jk) )
570            WRITE(numout,*)
571            WRITE(numout,*) ' level = ',jk
572            CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1,   &
573               &                              1, jpj, 1, 1, numout )
574         END DO
575         WRITE(numout,*)
576         WRITE(numout,*) ' dom_msk: bmask '
577         WRITE(numout,*) ' ---------------'
578         WRITE(numout,*)
579         imsk(:,:) = INT( bmask(:,:) )
580         CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1,   &
581            &                              1, jpj, 1, 1, numout )
582      ENDIF
583      !
584      CALL wrk_dealloc( jpi, jpj, imsk )
585      CALL wrk_dealloc( jpi, jpj, zwf  )
586      !
587      IF( nn_timing == 1 )  CALL timing_stop('dom_msk')
588      !
589   END SUBROUTINE dom_msk
590
591#if defined key_noslip_accurate
592   !!----------------------------------------------------------------------
593   !!   'key_noslip_accurate' :         accurate no-slip boundary condition
594   !!----------------------------------------------------------------------
595   
596   SUBROUTINE dom_msk_nsa
597      !!---------------------------------------------------------------------
598      !!                 ***  ROUTINE dom_msk_nsa  ***
599      !!
600      !! ** Purpose :
601      !!
602      !! ** Method  :
603      !!
604      !! ** Action :
605      !!----------------------------------------------------------------------
606      INTEGER  ::   ji, jj, jk, jl      ! dummy loop indices
607      INTEGER  ::   ine, inw, ins, inn, itest, ierror, iind, ijnd
608      REAL(wp) ::   zaa
609      !!---------------------------------------------------------------------
610      !
611      IF( nn_timing == 1 )  CALL timing_start('dom_msk_nsa')
612      !
613      IF(lwp) WRITE(numout,*)
614      IF(lwp) WRITE(numout,*) 'dom_msk_nsa : noslip accurate boundary condition'
615      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   using Schchepetkin and O Brian scheme'
616      IF( lk_mpp )   CALL ctl_stop('STOP', ' mpp version is not yet implemented' )
617
618      ! mask for second order calculation of vorticity
619      ! ----------------------------------------------
620      ! noslip boundary condition: fmask=1  at convex corner, store
621      ! index of straight coast meshes ( 'west', refering to a coast,
622      ! means west of the ocean, aso)
623     
624      DO jk = 1, jpk
625         DO jl = 1, 4
626            npcoa(jl,jk) = 0
627            DO ji = 1, 2*(jpi+jpj)
628               nicoa(ji,jl,jk) = 0
629               njcoa(ji,jl,jk) = 0
630            END DO
631         END DO
632      END DO
633     
634      IF( jperio == 2 ) THEN
635         WRITE(numout,*) ' '
636         WRITE(numout,*) ' symetric boundary conditions need special'
637         WRITE(numout,*) ' treatment not implemented. we stop.'
638         CALL ctl_stop('STOP', 'NEMO abort from dom_msk_nsa')
639      ENDIF
640     
641      ! convex corners
642     
643      DO jk = 1, jpkm1
644         DO jj = 1, jpjm1
645            DO ji = 1, jpim1
646               zaa = tmask(ji  ,jj,jk) + tmask(ji  ,jj+1,jk)   &
647                  &+ tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk)
648               IF( ABS(zaa-3._wp) <= 0.1_wp )   fmask(ji,jj,jk) = 1._wp
649            END DO
650         END DO
651      END DO
652
653      ! north-south straight coast
654
655      DO jk = 1, jpkm1
656         inw = 0
657         ine = 0
658         DO jj = 2, jpjm1
659            DO ji = 2, jpim1
660               zaa = tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk)
661               IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN
662                  inw = inw + 1
663                  nicoa(inw,1,jk) = ji
664                  njcoa(inw,1,jk) = jj
665                  IF( nprint == 1 ) WRITE(numout,*) ' west  : ', jk, inw, ji, jj
666               ENDIF
667               zaa = tmask(ji,jj,jk) + tmask(ji,jj+1,jk)
668               IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN
669                  ine = ine + 1
670                  nicoa(ine,2,jk) = ji
671                  njcoa(ine,2,jk) = jj
672                  IF( nprint == 1 ) WRITE(numout,*) ' east  : ', jk, ine, ji, jj
673               ENDIF
674            END DO
675         END DO
676         npcoa(1,jk) = inw
677         npcoa(2,jk) = ine
678      END DO
679
680      ! west-east straight coast
681
682      DO jk = 1, jpkm1
683         ins = 0
684         inn = 0
685         DO jj = 2, jpjm1
686            DO ji =2, jpim1
687               zaa = tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk)
688               IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN
689                  ins = ins + 1
690                  nicoa(ins,3,jk) = ji
691                  njcoa(ins,3,jk) = jj
692                  IF( nprint == 1 ) WRITE(numout,*) ' south : ', jk, ins, ji, jj
693               ENDIF
694               zaa = tmask(ji+1,jj,jk) + tmask(ji,jj,jk)
695               IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN
696                  inn = inn + 1
697                  nicoa(inn,4,jk) = ji
698                  njcoa(inn,4,jk) = jj
699                  IF( nprint == 1 ) WRITE(numout,*) ' north : ', jk, inn, ji, jj
700               ENDIF
701            END DO
702         END DO
703         npcoa(3,jk) = ins
704         npcoa(4,jk) = inn
705      END DO
706
707      itest = 2 * ( jpi + jpj )
708      DO jk = 1, jpk
709         IF( npcoa(1,jk) > itest .OR. npcoa(2,jk) > itest .OR.   &
710             npcoa(3,jk) > itest .OR. npcoa(4,jk) > itest ) THEN
711           
712            WRITE(ctmp1,*) ' level jk = ',jk
713            WRITE(ctmp2,*) ' straight coast index arraies are too small.:'
714            WRITE(ctmp3,*) ' npe, npw, nps, npn = ', npcoa(1,jk), npcoa(2,jk),   &
715                &                                     npcoa(3,jk), npcoa(4,jk)
716            WRITE(ctmp4,*) ' 2*(jpi+jpj) = ',itest,'. we stop.'
717            CALL ctl_stop( ctmp1, ctmp2, ctmp3, ctmp4 )
718        ENDIF
719      END DO
720
721      ierror = 0
722      iind = 0
723      ijnd = 0
724      IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 )   iind = 2
725      IF( nperio == 3 .OR. nperio == 4 .OR. nperio == 5 .OR. nperio == 6 )   ijnd = 2
726      DO jk = 1, jpk
727         DO jl = 1, npcoa(1,jk)
728            IF( nicoa(jl,1,jk)+3 > jpi+iind ) THEN
729               ierror = ierror+1
730               icoord(ierror,1) = nicoa(jl,1,jk)
731               icoord(ierror,2) = njcoa(jl,1,jk)
732               icoord(ierror,3) = jk
733            ENDIF
734         END DO
735         DO jl = 1, npcoa(2,jk)
736            IF(nicoa(jl,2,jk)-2 < 1-iind ) THEN
737               ierror = ierror + 1
738               icoord(ierror,1) = nicoa(jl,2,jk)
739               icoord(ierror,2) = njcoa(jl,2,jk)
740               icoord(ierror,3) = jk
741            ENDIF
742         END DO
743         DO jl = 1, npcoa(3,jk)
744            IF( njcoa(jl,3,jk)+3 > jpj+ijnd ) THEN
745               ierror = ierror + 1
746               icoord(ierror,1) = nicoa(jl,3,jk)
747               icoord(ierror,2) = njcoa(jl,3,jk)
748               icoord(ierror,3) = jk
749            ENDIF
750         END DO
751         DO jl = 1, npcoa(4,jk)
752            IF( njcoa(jl,4,jk)-2 < 1) THEN
753               ierror=ierror + 1
754               icoord(ierror,1) = nicoa(jl,4,jk)
755               icoord(ierror,2) = njcoa(jl,4,jk)
756               icoord(ierror,3) = jk
757            ENDIF
758         END DO
759      END DO
760     
761      IF( ierror > 0 ) THEN
762         IF(lwp) WRITE(numout,*)
763         IF(lwp) WRITE(numout,*) '              Problem on lateral conditions'
764         IF(lwp) WRITE(numout,*) '                 Bad marking off at points:'
765         DO jl = 1, ierror
766            IF(lwp) WRITE(numout,*) 'Level:',icoord(jl,3),   &
767               &                  '  Point(',icoord(jl,1),',',icoord(jl,2),')'
768         END DO
769         CALL ctl_stop( 'We stop...' )
770      ENDIF
771      !
772      IF( nn_timing == 1 )  CALL timing_stop('dom_msk_nsa')
773      !
774   END SUBROUTINE dom_msk_nsa
775
776#else
777   !!----------------------------------------------------------------------
778   !!   Default option :                                      Empty routine
779   !!----------------------------------------------------------------------
780   SUBROUTINE dom_msk_nsa       
781   END SUBROUTINE dom_msk_nsa
782#endif
783   
784   !!======================================================================
785END MODULE dommsk
Note: See TracBrowser for help on using the repository browser.