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

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

Create 3D U and V interior masks to mask out wrap columns
and 1st and last rows. Apply to appropriate grids
at dfinition phase.

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