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

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

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