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

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

Use existing dom_uniq routine to mask duplicate points.
This routine does essentailly what I already do and
although it's slightly excessive for what we need,
it allows us to reduce the number of lines
of code we need by calling an existing routine
(which is tried and tested!).
File size: 33.4 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      CALL lbc_lnk( vmask_i_diag, 'V', 1. )
319
320      ! Now mask out any duplicate points
321      CALL dom_uniq( uvt, 'V' )
322      DO jk = 1, jpk
323         vmask_i_diag(:,:,jk) = vmask_i_diag(:,:,jk) * uvt(:,:)
324      END DO
325
326
327
328      ! 3. Ocean/land mask at wu-, wv- and w points
329      !----------------------------------------------
330      wmask (:,:,1) = tmask(:,:,1) ! ????????
331      wumask(:,:,1) = umask(:,:,1) ! ????????
332      wvmask(:,:,1) = vmask(:,:,1) ! ????????
333      DO jk=2,jpk
334         wmask (:,:,jk)=tmask(:,:,jk) * tmask(:,:,jk-1)
335         wumask(:,:,jk)=umask(:,:,jk) * umask(:,:,jk-1)   
336         wvmask(:,:,jk)=vmask(:,:,jk) * vmask(:,:,jk-1)
337      END DO
338
339      ! 4. ocean/land mask for the elliptic equation
340      ! --------------------------------------------
341      bmask(:,:) = ssmask(:,:)       ! elliptic equation is written at t-point
342      !
343      !                               ! Boundary conditions
344      !                                    ! cyclic east-west : bmask must be set to 0. on rows 1 and jpi
345      IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN
346         bmask( 1 ,:) = 0._wp
347         bmask(jpi,:) = 0._wp
348      ENDIF
349      IF( nperio == 2 ) THEN               ! south symmetric :  bmask must be set to 0. on row 1
350         bmask(:, 1 ) = 0._wp
351      ENDIF
352      !                                    ! north fold :
353      IF( nperio == 3 .OR. nperio == 4 ) THEN   ! T-pt pivot : bmask set to 0. on row jpj and on half jpjglo-1 row
354         DO ji = 1, jpi                     
355            ii = ji + nimpp - 1
356            bmask(ji,jpj-1) = bmask(ji,jpj-1) * tpol(ii)
357            bmask(ji,jpj  ) = 0._wp
358         END DO
359      ENDIF
360      IF( nperio == 5 .OR. nperio == 6 ) THEN   ! F-pt pivot and T-pt elliptic eq. : bmask set to 0. on row jpj
361         bmask(:,jpj) = 0._wp
362      ENDIF
363      !
364      IF( lk_mpp ) THEN                    ! mpp specificities
365         !                                      ! bmask is set to zero on the overlap region
366         IF( nbondi /= -1 .AND. nbondi /= 2 )   bmask(  1 :jpreci,:) = 0._wp
367         IF( nbondi /=  1 .AND. nbondi /= 2 )   bmask(nlci:jpi   ,:) = 0._wp
368         IF( nbondj /= -1 .AND. nbondj /= 2 )   bmask(:,  1 :jprecj) = 0._wp
369         IF( nbondj /=  1 .AND. nbondj /= 2 )   bmask(:,nlcj:jpj   ) = 0._wp
370         !
371         IF( npolj == 3 .OR. npolj == 4 ) THEN  ! north fold : bmask must be set to 0. on rows jpj-1 and jpj
372            DO ji = 1, nlci
373               ii = ji + nimpp - 1
374               bmask(ji,nlcj-1) = bmask(ji,nlcj-1) * tpol(ii)
375               bmask(ji,nlcj  ) = 0._wp
376            END DO
377         ENDIF
378         IF( npolj == 5 .OR. npolj == 6 ) THEN  ! F-pt pivot and T-pt elliptic eq. : bmask set to 0. on row jpj
379            DO ji = 1, nlci
380               bmask(ji,nlcj  ) = 0._wp
381            END DO
382         ENDIF
383      ENDIF
384
385
386      ! mask for second order calculation of vorticity
387      ! ----------------------------------------------
388      CALL dom_msk_nsa
389
390     
391      ! Lateral boundary conditions on velocity (modify fmask)
392      ! ---------------------------------------     
393      DO jk = 1, jpk
394         zwf(:,:) = fmask(:,:,jk)         
395         DO jj = 2, jpjm1
396            DO ji = fs_2, fs_jpim1   ! vector opt.
397               IF( fmask(ji,jj,jk) == 0._wp ) THEN
398                  fmask(ji,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,jj), zwf(ji,jj+1),   &
399                     &                                           zwf(ji-1,jj), zwf(ji,jj-1)  )  )
400               ENDIF
401            END DO
402         END DO
403         DO jj = 2, jpjm1
404            IF( fmask(1,jj,jk) == 0._wp ) THEN
405               fmask(1  ,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1) ) )
406            ENDIF
407            IF( fmask(jpi,jj,jk) == 0._wp ) THEN
408               fmask(jpi,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) )
409            ENDIF
410         END DO         
411         DO ji = 2, jpim1
412            IF( fmask(ji,1,jk) == 0._wp ) THEN
413               fmask(ji, 1 ,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1) ) )
414            ENDIF
415            IF( fmask(ji,jpj,jk) == 0._wp ) THEN
416               fmask(ji,jpj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpjm1) ) )
417            ENDIF
418         END DO
419      END DO
420      !
421      IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN   ! ORCA_R2 configuration
422         !                                                 ! Increased lateral friction near of some straits
423         IF( nn_cla == 0 ) THEN
424            !                                ! Gibraltar strait  : partial slip (fmask=0.5)
425            ij0 = 101   ;   ij1 = 101
426            ii0 = 139   ;   ii1 = 140   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  0.5_wp
427            ij0 = 102   ;   ij1 = 102
428            ii0 = 139   ;   ii1 = 140   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  0.5_wp
429            !
430            !                                ! Bab el Mandeb : partial slip (fmask=1)
431            ij0 =  87   ;   ij1 =  88
432            ii0 = 160   ;   ii1 = 160   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  1._wp
433            ij0 =  88   ;   ij1 =  88
434            ii0 = 159   ;   ii1 = 159   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  1._wp
435            !
436         ENDIF
437         !                                ! Danish straits  : strong slip (fmask > 2)
438! We keep this as an example but it is instable in this case
439!         ij0 = 115   ;   ij1 = 115
440!         ii0 = 145   ;   ii1 = 146   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 4._wp
441!         ij0 = 116   ;   ij1 = 116
442!         ii0 = 145   ;   ii1 = 146   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 4._wp
443         !
444      ENDIF
445      !
446      IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN   ! ORCA R1 configuration
447         !                                                 ! Increased lateral friction near of some straits
448         ! This dirty section will be suppressed by simplification process:
449         ! all this will come back in input files
450         ! Currently these hard-wired indices relate to configuration with
451         ! extend grid (jpjglo=332)
452         !
453         isrow = 332 - jpjglo
454         !
455         IF(lwp) WRITE(numout,*)
456         IF(lwp) WRITE(numout,*) '   orca_r1: increase friction near the following straits : '
457         IF(lwp) WRITE(numout,*) '      Gibraltar '
458         ii0 = 282           ;   ii1 = 283        ! Gibraltar Strait
459         ij0 = 241 - isrow   ;   ij1 = 241 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp 
460
461         IF(lwp) WRITE(numout,*) '      Bhosporus '
462         ii0 = 314           ;   ii1 = 315        ! Bhosporus Strait
463         ij0 = 248 - isrow   ;   ij1 = 248 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp 
464
465         IF(lwp) WRITE(numout,*) '      Makassar (Top) '
466         ii0 =  48           ;   ii1 =  48        ! Makassar Strait (Top)
467         ij0 = 189 - isrow   ;   ij1 = 190 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp 
468
469         IF(lwp) WRITE(numout,*) '      Lombok '
470         ii0 =  44           ;   ii1 =  44        ! Lombok Strait
471         ij0 = 164 - isrow   ;   ij1 = 165 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp 
472
473         IF(lwp) WRITE(numout,*) '      Ombai '
474         ii0 =  53           ;   ii1 =  53        ! Ombai Strait
475         ij0 = 164 - isrow   ;   ij1 = 165 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp 
476
477         IF(lwp) WRITE(numout,*) '      Timor Passage '
478         ii0 =  56           ;   ii1 =  56        ! Timor Passage
479         ij0 = 164 - isrow   ;   ij1 = 165 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp 
480
481         IF(lwp) WRITE(numout,*) '      West Halmahera '
482         ii0 =  58           ;   ii1 =  58        ! West Halmahera Strait
483         ij0 = 181 - isrow   ;   ij1 = 182 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp 
484
485         IF(lwp) WRITE(numout,*) '      East Halmahera '
486         ii0 =  55           ;   ii1 =  55        ! East Halmahera Strait
487         ij0 = 181 - isrow   ;   ij1 = 182 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp 
488         !
489      ENDIF
490      !
491      IF( cp_cfg == "orca" .AND. jp_cfg == 025 .AND. rn_shlat == 0.0 ) THEN   
492         !                                              ! ORCA_R025 configuration
493         !                                              ! Increased lateral friction on parts of Antarctic coastline
494         !                                              ! for increased stability
495         !                                              ! NB. This only works to do this here if we have free slip
496         !                                              ! generally, so fmask is zero at coast points.
497         IF(lwp) WRITE(numout,*)
498         IF(lwp) WRITE(numout,*) '   orca_r025: increase friction in following regions : '
499         IF(lwp) WRITE(numout,*) '      whole Antarctic coastline: partial slip shlat=1 '
500
501         zphi_drake_passage = -58.0_wp
502         zshlat_antarc = 1.0_wp
503         zwf(:,:) = fmask(:,:,1)         
504         DO jj = 2, jpjm1
505            DO ji = fs_2, fs_jpim1   ! vector opt.
506               IF( gphif(ji,jj) .lt. zphi_drake_passage .and. fmask(ji,jj,1) == 0._wp ) THEN
507                  fmask(ji,jj,:) = zshlat_antarc * MIN( 1._wp , MAX( zwf(ji+1,jj), zwf(ji,jj+1),   &
508                     &                                           zwf(ji-1,jj), zwf(ji,jj-1)  )  )
509               ENDIF
510            END DO
511         END DO
512      END IF
513      !
514      CALL lbc_lnk( fmask, 'F', 1._wp )      ! Lateral boundary conditions on fmask
515
516      ! CAUTION : The fmask may be further modified in dyn_vor_init ( dynvor.F90 )
517           
518      IF( nprint == 1 .AND. lwp ) THEN      ! Control print
519         imsk(:,:) = INT( tmask_i(:,:) )
520         WRITE(numout,*) ' tmask_i : '
521         CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1,   &
522               &                           1, jpj, 1, 1, numout)
523         WRITE (numout,*)
524         WRITE (numout,*) ' dommsk: tmask for each level'
525         WRITE (numout,*) ' ----------------------------'
526         DO jk = 1, jpk
527            imsk(:,:) = INT( tmask(:,:,jk) )
528
529            WRITE(numout,*)
530            WRITE(numout,*) ' level = ',jk
531            CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1,   &
532               &                              1, jpj, 1, 1, numout)
533         END DO
534         WRITE(numout,*)
535         WRITE(numout,*) ' dom_msk: vmask for each level'
536         WRITE(numout,*) ' -----------------------------'
537         DO jk = 1, jpk
538            imsk(:,:) = INT( vmask(:,:,jk) )
539            WRITE(numout,*)
540            WRITE(numout,*) ' level = ',jk
541            CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1,   &
542               &                              1, jpj, 1, 1, numout)
543         END DO
544         WRITE(numout,*)
545         WRITE(numout,*) ' dom_msk: fmask for each level'
546         WRITE(numout,*) ' -----------------------------'
547         DO jk = 1, jpk
548            imsk(:,:) = INT( fmask(:,:,jk) )
549            WRITE(numout,*)
550            WRITE(numout,*) ' level = ',jk
551            CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1,   &
552               &                              1, jpj, 1, 1, numout )
553         END DO
554         WRITE(numout,*)
555         WRITE(numout,*) ' dom_msk: bmask '
556         WRITE(numout,*) ' ---------------'
557         WRITE(numout,*)
558         imsk(:,:) = INT( bmask(:,:) )
559         CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1,   &
560            &                              1, jpj, 1, 1, numout )
561      ENDIF
562      !
563      CALL wrk_dealloc( jpi, jpj, imsk )
564      CALL wrk_dealloc( jpi, jpj, zwf  )
565      !
566      IF( nn_timing == 1 )  CALL timing_stop('dom_msk')
567      !
568   END SUBROUTINE dom_msk
569
570#if defined key_noslip_accurate
571   !!----------------------------------------------------------------------
572   !!   'key_noslip_accurate' :         accurate no-slip boundary condition
573   !!----------------------------------------------------------------------
574   
575   SUBROUTINE dom_msk_nsa
576      !!---------------------------------------------------------------------
577      !!                 ***  ROUTINE dom_msk_nsa  ***
578      !!
579      !! ** Purpose :
580      !!
581      !! ** Method  :
582      !!
583      !! ** Action :
584      !!----------------------------------------------------------------------
585      INTEGER  ::   ji, jj, jk, jl      ! dummy loop indices
586      INTEGER  ::   ine, inw, ins, inn, itest, ierror, iind, ijnd
587      REAL(wp) ::   zaa
588      !!---------------------------------------------------------------------
589      !
590      IF( nn_timing == 1 )  CALL timing_start('dom_msk_nsa')
591      !
592      IF(lwp) WRITE(numout,*)
593      IF(lwp) WRITE(numout,*) 'dom_msk_nsa : noslip accurate boundary condition'
594      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   using Schchepetkin and O Brian scheme'
595      IF( lk_mpp )   CALL ctl_stop('STOP', ' mpp version is not yet implemented' )
596
597      ! mask for second order calculation of vorticity
598      ! ----------------------------------------------
599      ! noslip boundary condition: fmask=1  at convex corner, store
600      ! index of straight coast meshes ( 'west', refering to a coast,
601      ! means west of the ocean, aso)
602     
603      DO jk = 1, jpk
604         DO jl = 1, 4
605            npcoa(jl,jk) = 0
606            DO ji = 1, 2*(jpi+jpj)
607               nicoa(ji,jl,jk) = 0
608               njcoa(ji,jl,jk) = 0
609            END DO
610         END DO
611      END DO
612     
613      IF( jperio == 2 ) THEN
614         WRITE(numout,*) ' '
615         WRITE(numout,*) ' symetric boundary conditions need special'
616         WRITE(numout,*) ' treatment not implemented. we stop.'
617         CALL ctl_stop('STOP', 'NEMO abort from dom_msk_nsa')
618      ENDIF
619     
620      ! convex corners
621     
622      DO jk = 1, jpkm1
623         DO jj = 1, jpjm1
624            DO ji = 1, jpim1
625               zaa = tmask(ji  ,jj,jk) + tmask(ji  ,jj+1,jk)   &
626                  &+ tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk)
627               IF( ABS(zaa-3._wp) <= 0.1_wp )   fmask(ji,jj,jk) = 1._wp
628            END DO
629         END DO
630      END DO
631
632      ! north-south straight coast
633
634      DO jk = 1, jpkm1
635         inw = 0
636         ine = 0
637         DO jj = 2, jpjm1
638            DO ji = 2, jpim1
639               zaa = tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk)
640               IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN
641                  inw = inw + 1
642                  nicoa(inw,1,jk) = ji
643                  njcoa(inw,1,jk) = jj
644                  IF( nprint == 1 ) WRITE(numout,*) ' west  : ', jk, inw, ji, jj
645               ENDIF
646               zaa = tmask(ji,jj,jk) + tmask(ji,jj+1,jk)
647               IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN
648                  ine = ine + 1
649                  nicoa(ine,2,jk) = ji
650                  njcoa(ine,2,jk) = jj
651                  IF( nprint == 1 ) WRITE(numout,*) ' east  : ', jk, ine, ji, jj
652               ENDIF
653            END DO
654         END DO
655         npcoa(1,jk) = inw
656         npcoa(2,jk) = ine
657      END DO
658
659      ! west-east straight coast
660
661      DO jk = 1, jpkm1
662         ins = 0
663         inn = 0
664         DO jj = 2, jpjm1
665            DO ji =2, jpim1
666               zaa = tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk)
667               IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN
668                  ins = ins + 1
669                  nicoa(ins,3,jk) = ji
670                  njcoa(ins,3,jk) = jj
671                  IF( nprint == 1 ) WRITE(numout,*) ' south : ', jk, ins, ji, jj
672               ENDIF
673               zaa = tmask(ji+1,jj,jk) + tmask(ji,jj,jk)
674               IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN
675                  inn = inn + 1
676                  nicoa(inn,4,jk) = ji
677                  njcoa(inn,4,jk) = jj
678                  IF( nprint == 1 ) WRITE(numout,*) ' north : ', jk, inn, ji, jj
679               ENDIF
680            END DO
681         END DO
682         npcoa(3,jk) = ins
683         npcoa(4,jk) = inn
684      END DO
685
686      itest = 2 * ( jpi + jpj )
687      DO jk = 1, jpk
688         IF( npcoa(1,jk) > itest .OR. npcoa(2,jk) > itest .OR.   &
689             npcoa(3,jk) > itest .OR. npcoa(4,jk) > itest ) THEN
690           
691            WRITE(ctmp1,*) ' level jk = ',jk
692            WRITE(ctmp2,*) ' straight coast index arraies are too small.:'
693            WRITE(ctmp3,*) ' npe, npw, nps, npn = ', npcoa(1,jk), npcoa(2,jk),   &
694                &                                     npcoa(3,jk), npcoa(4,jk)
695            WRITE(ctmp4,*) ' 2*(jpi+jpj) = ',itest,'. we stop.'
696            CALL ctl_stop( ctmp1, ctmp2, ctmp3, ctmp4 )
697        ENDIF
698      END DO
699
700      ierror = 0
701      iind = 0
702      ijnd = 0
703      IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 )   iind = 2
704      IF( nperio == 3 .OR. nperio == 4 .OR. nperio == 5 .OR. nperio == 6 )   ijnd = 2
705      DO jk = 1, jpk
706         DO jl = 1, npcoa(1,jk)
707            IF( nicoa(jl,1,jk)+3 > jpi+iind ) THEN
708               ierror = ierror+1
709               icoord(ierror,1) = nicoa(jl,1,jk)
710               icoord(ierror,2) = njcoa(jl,1,jk)
711               icoord(ierror,3) = jk
712            ENDIF
713         END DO
714         DO jl = 1, npcoa(2,jk)
715            IF(nicoa(jl,2,jk)-2 < 1-iind ) THEN
716               ierror = ierror + 1
717               icoord(ierror,1) = nicoa(jl,2,jk)
718               icoord(ierror,2) = njcoa(jl,2,jk)
719               icoord(ierror,3) = jk
720            ENDIF
721         END DO
722         DO jl = 1, npcoa(3,jk)
723            IF( njcoa(jl,3,jk)+3 > jpj+ijnd ) THEN
724               ierror = ierror + 1
725               icoord(ierror,1) = nicoa(jl,3,jk)
726               icoord(ierror,2) = njcoa(jl,3,jk)
727               icoord(ierror,3) = jk
728            ENDIF
729         END DO
730         DO jl = 1, npcoa(4,jk)
731            IF( njcoa(jl,4,jk)-2 < 1) THEN
732               ierror=ierror + 1
733               icoord(ierror,1) = nicoa(jl,4,jk)
734               icoord(ierror,2) = njcoa(jl,4,jk)
735               icoord(ierror,3) = jk
736            ENDIF
737         END DO
738      END DO
739     
740      IF( ierror > 0 ) THEN
741         IF(lwp) WRITE(numout,*)
742         IF(lwp) WRITE(numout,*) '              Problem on lateral conditions'
743         IF(lwp) WRITE(numout,*) '                 Bad marking off at points:'
744         DO jl = 1, ierror
745            IF(lwp) WRITE(numout,*) 'Level:',icoord(jl,3),   &
746               &                  '  Point(',icoord(jl,1),',',icoord(jl,2),')'
747         END DO
748         CALL ctl_stop( 'We stop...' )
749      ENDIF
750      !
751      IF( nn_timing == 1 )  CALL timing_stop('dom_msk_nsa')
752      !
753   END SUBROUTINE dom_msk_nsa
754
755#else
756   !!----------------------------------------------------------------------
757   !!   Default option :                                      Empty routine
758   !!----------------------------------------------------------------------
759   SUBROUTINE dom_msk_nsa       
760   END SUBROUTINE dom_msk_nsa
761#endif
762   
763   !!======================================================================
764END MODULE dommsk
Note: See TracBrowser for help on using the repository browser.