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

source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90 @ 9830

Last change on this file since 9830 was 9830, checked in by frrh, 6 years ago

Merge revisions 9607:9721 of/branches/UKMO/dev_r5518_GO6_diag_bitcomp
into GO6 package branch.

This change ensures most 2D and 3D diagnostics produced by NEMO and MEDUSA
are bit reproducible on different PE decompositions.

Command used:
svn merge -r 9607:9721 ​svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/branches/UKMO/dev_r5518_GO6_diag_bitcomp

Met Office GMED ticket 389 refers.

This change applies a mask to all duplicate grid points
output on diagnistic grids for T, U and V points. i.e. it masks
any wrap columns and duplicated grid points across the N-fold.
Fields affected are all "standard" NEMO diagnostics (scalar and
diaptr diagnostics are not on "normal" grids).

It also introduces some corrections/initialisations to achieve
PE decomposition bit comparison.

Most 2D or 3D fields are now bit comparable on different PE decompositions.
Only diaptr diagnostics can not be guaranteed bit reproducible
(due to their method of computation).

This change does nothing to CICE output.

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