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

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

Correct U mask footprint and remove N-fold specific masking on
T points in out new 3d T mask.

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