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

source: trunk/NEMO/OPA_SRC/DOM/dommsk.F90 @ 1341

Last change on this file since 1341 was 1273, checked in by ctlod, 15 years ago

update Gibrlatar, Bab El Mandeb and Sound straits in both full & partial steps bathymetry files such as closed seas, see ticket: #305
  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 26.1 KB
Line 
1MODULE dommsk
2   !!==============================================================================
3   !!                       ***  MODULE dommsk   ***
4   !! Ocean initialization : domain land/sea mask
5   !!==============================================================================
6
7   !!----------------------------------------------------------------------
8   !!   dom_msk        : compute land/ocean mask
9   !!   dom_msk_nsa    : update land/ocean mask when no-slip accurate
10   !!                    option is used.
11   !!----------------------------------------------------------------------
12   !! * Modules used
13   USE oce             ! ocean dynamics and tracers
14   USE dom_oce         ! ocean space and time domain
15   USE obc_oce         ! ocean open boundary conditions
16   USE in_out_manager  ! I/O manager
17   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
18   USE lib_mpp
19   USE solisl          ! ???
20   USE dynspg_oce      ! choice/control of key cpp for surface pressure gradient
21
22   IMPLICIT NONE
23   PRIVATE
24
25   !! * Routine accessibility
26   PUBLIC dom_msk        ! routine called by inidom.F90
27
28   !! * Module variables
29   REAL(wp) ::   &
30      shlat = 2.   ! type of lateral boundary condition on velocity (namelist namlbc)
31   
32   !! * Substitutions
33#  include "vectopt_loop_substitute.h90"
34   !!---------------------------------------------------------------------------------
35   !!   OPA 9.0 , LOCEAN-IPSL (2005)
36   !! $Id$
37   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
38   !!---------------------------------------------------------------------------------
39
40CONTAINS
41   
42   SUBROUTINE dom_msk
43      !!---------------------------------------------------------------------
44      !!                 ***  ROUTINE dom_msk  ***
45      !!
46      !! ** Purpose :   Compute land/ocean mask arrays at tracer points, hori-
47      !!      zontal velocity points (u & v), vorticity points (f) and baro-
48      !!      tropic stream function  points (b).
49      !!        Set mbathy to the number of non-zero w-levels of a water column
50      !!      (if island in the domain (lk_isl=T), this is done latter in
51      !!      routine solver_init)
52      !!
53      !! ** Method  :   The ocean/land mask is computed from the basin bathy-
54      !!      metry in level (mbathy) which is defined or read in dommba.
55      !!      mbathy equals 0 over continental T-point, -n over the nth
56      !!      island T-point, and the number of ocean level over the ocean.
57      !!
58      !!      At a given position (ji,jj,jk) the ocean/land mask is given by:
59      !!      t-point : 0. IF mbathy( ji ,jj) =< 0
60      !!                1. IF mbathy( ji ,jj) >= jk
61      !!      u-point : 0. IF mbathy( ji ,jj)  or mbathy(ji+1, jj ) =< 0
62      !!                1. IF mbathy( ji ,jj) and mbathy(ji+1, jj ) >= jk.
63      !!      v-point : 0. IF mbathy( ji ,jj)  or mbathy( ji ,jj+1) =< 0
64      !!                1. IF mbathy( ji ,jj) and mbathy( ji ,jj+1) >= jk.
65      !!      f-point : 0. IF mbathy( ji ,jj)  or mbathy( ji ,jj+1)
66      !!                   or mbathy(ji+1,jj)  or mbathy(ji+1,jj+1) =< 0
67      !!                1. IF mbathy( ji ,jj) and mbathy( ji ,jj+1)
68      !!                and mbathy(ji+1,jj) and mbathy(ji+1,jj+1) >= jk.
69      !!      b-point : the same definition as for f-point of the first ocean
70      !!                level (surface level) but with 0 along coastlines.
71      !!
72      !!        The lateral friction is set through the value of fmask along
73      !!      the coast and topography. This value is defined by shlat, a
74      !!      namelist parameter:
75      !!         shlat = 0, free slip  (no shear along the coast)
76      !!         shlat = 2, no slip  (specified zero velocity at the coast)
77      !!         0 < shlat < 2, partial slip   | non-linear velocity profile
78      !!         2 < shlat, strong slip        | in the lateral boundary layer
79      !!
80      !!      N.B. If nperio not equal to 0, the land/ocean mask arrays
81      !!      are defined with the proper value at lateral domain boundaries,
82      !!      but bmask. indeed, bmask defined the domain over which the
83      !!      barotropic stream function is computed. this domain cannot
84      !!      contain identical columns because the matrix associated with
85      !!      the barotropic stream function equation is then no more inverti-
86      !!      ble. therefore bmask is set to 0 along lateral domain boundaries
87      !!      even IF nperio is not zero.
88      !!
89      !!      In case of open boundaries (lk_obc=T):
90      !!        - tmask is set to 1 on the points to be computed bay the open
91      !!          boundaries routines.
92      !!        - bmask is  set to 0 on the open boundaries.
93      !!
94      !!      Set mbathy to the number of non-zero w-levels of a water column
95      !!                  mbathy = min( mbathy, 1 ) + 1
96      !!      (note that the minimum value of mbathy is 2).
97      !!
98      !! ** Action :
99      !!                     tmask    : land/ocean mask at t-point (=0. or 1.)
100      !!                     umask    : land/ocean mask at u-point (=0. or 1.)
101      !!                     vmask    : land/ocean mask at v-point (=0. or 1.)
102      !!                     fmask    : land/ocean mask at f-point (=0. or 1.)
103      !!                          =shlat along lateral boundaries
104      !!                     bmask    : land/ocean mask at barotropic stream
105      !!                          function point (=0. or 1.) and set to
106      !!                          0 along lateral boundaries
107      !!                   mbathy   : number of non-zero w-levels
108      !!
109      !! History :
110      !!        !  87-07  (G. Madec)  Original code
111      !!        !  91-12  (G. Madec)
112      !!        !  92-06  (M. Imbard)
113      !!        !  93-03  (M. Guyon)  symetrical conditions (M. Guyon)
114      !!        !  96-01  (G. Madec)  suppression of common work arrays
115      !!        !  96-05  (G. Madec)  mask computed from tmask and sup-
116      !!                 pression of the double computation of bmask
117      !!        !  97-02  (G. Madec)  mesh information put in domhgr.F
118      !!        !  97-07  (G. Madec)  modification of mbathy and fmask
119      !!        !  98-05  (G. Roullet)  free surface
120      !!        !  00-03  (G. Madec)  no slip accurate
121      !!        !  01-09  (J.-M. Molines)  Open boundaries
122      !!   8.5  !  02-08  (G. Madec)  F90: Free form and module
123      !!   9.0  !  05-11  (V. Garnier) Surface pressure gradient organization
124      !!----------------------------------------------------------------------
125      INTEGER  ::   ji, jj, jk      ! dummy loop indices
126      INTEGER  ::   iif, iil, ii0, ii1, ii
127      INTEGER  ::   ijf, ijl, ij0, ij1
128      INTEGER, DIMENSION(jpi,jpj) ::  imsk
129      REAL(wp), DIMENSION(jpi,jpj) ::   zwf
130
131      NAMELIST/namlbc/ shlat
132      !!---------------------------------------------------------------------
133     
134      ! Namelist namlbc : lateral momentum boundary condition
135      REWIND( numnam )
136      READ  ( numnam, namlbc )
137      IF(lwp) THEN
138         WRITE(numout,*)
139         WRITE(numout,*) 'dommsk : ocean mask '
140         WRITE(numout,*) '~~~~~~'
141         WRITE(numout,*) '         Namelist namlbc'
142         WRITE(numout,*) '            lateral momentum boundary cond. shlat = ',shlat
143      ENDIF
144
145      IF ( shlat == 0. ) THEN
146          IF(lwp) WRITE(numout,*) '         ocean lateral free-slip '
147        ELSEIF ( shlat  ==  2. ) THEN
148          IF(lwp) WRITE(numout,*) '         ocean lateral  no-slip '
149        ELSEIF ( 0. < shlat .AND. shlat < 2. ) THEN
150          IF(lwp) WRITE(numout,*) '         ocean lateral  partial-slip '
151        ELSEIF ( 2. < shlat ) THEN
152          IF(lwp) WRITE(numout,*) '         ocean lateral  strong-slip '
153        ELSE
154          WRITE(ctmp1,*) ' shlat is negative = ', shlat
155          CALL ctl_stop( ctmp1 )
156      ENDIF
157
158      ! 1. Ocean/land mask at t-point (computed from mbathy)
159      ! -----------------------------
160      ! Tmask has already the right boundary conditions since mbathy is ok
161
162      tmask(:,:,:) = 0.e0
163      DO jk = 1, jpk
164         DO jj = 1, jpj
165            DO ji = 1, jpi
166               IF( FLOAT( mbathy(ji,jj)-jk )+.1 >= 0.e0 ) tmask(ji,jj,jk) = 1.e0
167            END DO 
168         END DO 
169      END DO 
170
171#if defined key_zdfkpp
172      IF( cp_cfg == 'orca' )   THEN
173         IF( jp_cfg == 2 )   THEN
174            ! land point on Bab el Mandeb zonal section
175            ij0 =  87   ;   ij1 =  88
176            ii0 = 160   ;   ii1 = 161
177            tmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 0.e0
178         ELSE
179            IF(lwp) WRITE(numout,*)
180            IF(lwp) WRITE(numout,cform_war)
181            IF(lwp) WRITE(numout,*)
182            IF(lwp) WRITE(numout,*)'          A mask must be applied on Bab el Mandeb strait'
183            IF(lwp) WRITE(numout,*)'          in case of ORCAs configurations'
184            IF(lwp) WRITE(numout,*)'          This is a problem which is not yet solved'
185            IF(lwp) WRITE(numout,*)
186         ENDIF
187      ENDIF
188#endif
189
190      ! Interior domain mask (used for global sum)
191      ! --------------------
192
193      tmask_i(:,:) = tmask(:,:,1)
194      iif = jpreci                         ! ???
195      iil = nlci - jpreci + 1
196      ijf = jprecj                         ! ???
197      ijl = nlcj - jprecj + 1
198
199      tmask_i( 1 :iif,   :   ) = 0.e0      ! first columns
200      tmask_i(iil:jpi,   :   ) = 0.e0      ! last  columns (including mpp extra columns)
201      tmask_i(   :   , 1 :ijf) = 0.e0      ! first rows
202      tmask_i(   :   ,ijl:jpj) = 0.e0      ! last  rows (including mpp extra rows)
203
204      ! north fold mask
205      tpol(1:jpiglo) = 1.e0 
206      fpol(1:jpiglo) = 1.e0
207      IF( jperio == 3 .OR. jperio == 4 ) THEN      ! T-point pivot
208         tpol(jpiglo/2+1:jpiglo) = 0.e0
209         fpol(     1    :jpiglo) = 0.e0
210         ! T-point pivot: only half of the nlcj-1 row
211         IF( mjg(nlej) == jpjglo )   THEN
212            DO ji = iif+1, iil-1
213               tmask_i(ji,nlej-1) = tmask_i(ji,nlej-1) * tpol(mig(ji))
214            END DO
215         ENDIF
216      ENDIF
217      IF( jperio == 5 .OR. jperio == 6 ) THEN      ! F-point pivot
218         tpol(     1    :jpiglo) = 0.e0
219         fpol(jpiglo/2+1:jpiglo) = 0.e0
220      ENDIF
221
222      ! 2. Ocean/land mask at u-,  v-, and z-points (computed from tmask)
223      ! -------------------------------------------
224     
225      ! Computation
226      DO jk = 1, jpk
227         DO jj = 1, jpjm1
228            DO ji = 1, fs_jpim1   ! vector loop
229               umask(ji,jj,jk) = tmask(ji,jj  ,jk) * tmask(ji+1,jj  ,jk)
230               vmask(ji,jj,jk) = tmask(ji,jj  ,jk) * tmask(ji  ,jj+1,jk)
231            END DO
232            DO ji = 1, jpim1 
233               fmask(ji,jj,jk) = tmask(ji,jj  ,jk) * tmask(ji+1,jj  ,jk)   &
234                  &            * tmask(ji,jj+1,jk) * tmask(ji+1,jj+1,jk)
235            END DO
236         END DO
237      END DO
238
239      ! Lateral boundary conditions
240      CALL lbc_lnk( umask, 'U', 1. )
241      CALL lbc_lnk( vmask, 'V', 1. )
242      CALL lbc_lnk( fmask, 'F', 1. )
243
244
245      ! 4. ocean/land mask for the elliptic equation
246      ! --------------------------------------------
247     
248      ! Computation
249      IF( lk_dynspg_rl ) THEN
250         bmask(:,:) = fmask(:,:,1)       ! elliptic equation is written at f-point
251      ELSE
252         bmask(:,:) = tmask(:,:,1)       ! elliptic equation is written at t-point
253      ENDIF
254     
255      ! Boundary conditions
256      !   cyclic east-west : bmask must be set to 0. on rows 1 and jpi
257      IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN
258         bmask( 1 ,:) = 0.e0
259         bmask(jpi,:) = 0.e0
260      ENDIF
261     
262      !   south symmetric :  bmask must be set to 0. on row 1
263      IF( nperio == 2 ) THEN
264         bmask(:, 1 ) = 0.e0
265      ENDIF
266     
267      !   north fold :
268      IF( nperio == 3 .OR. nperio == 4 ) THEN
269         IF( lk_dynspg_rl ) THEN
270            ! T-pt pivot and F-pt elliptic eq. : bmask set to 0. on rows jpj-1 and jpj
271            bmask(:,jpj-1) = 0.e0
272            bmask(:,jpj  ) = 0.e0
273         ELSE
274            ! T-pt pivot and T-pt elliptic eq. : bmask set to 0. on row jpj and on half jpjglo-1 row
275            DO ji = 1, jpi
276               ii = ji + nimpp - 1
277               bmask(ji,jpj-1) = bmask(ji,jpj-1) * tpol(ii)
278               bmask(ji,jpj  ) = 0.e0
279            END DO
280         ENDIF
281      ENDIF
282      IF( nperio == 5 .OR. nperio == 6 ) THEN
283         IF( lk_dynspg_rl ) THEN
284            ! F-pt pivot and F-pt elliptic eq. : bmask set to 0. on row jpj and on half jpjglo-1 row
285            DO ji = 1, jpi
286               ii = ji + nimpp - 1
287               bmask(ji,jpj-1) = bmask(ji,jpj-1) * fpol(ii)
288               bmask(ji,jpj  ) = 0.e0
289            END DO
290         ELSE
291            ! F-pt pivot and T-pt elliptic eq. : bmask set to 0. on row jpj
292            bmask(:,jpj) = 0.e0
293         ENDIF
294      ENDIF
295
296      ! Mpp boundary conditions: bmask is set to zero on the overlap
297      ! region for all elliptic solvers
298
299      IF( lk_mpp ) THEN
300         IF( nbondi /= -1 .AND. nbondi /= 2 )   bmask(  1 :jpreci,:) = 0.e0
301         IF( nbondi /=  1 .AND. nbondi /= 2 )   bmask(nlci:jpi   ,:) = 0.e0
302         IF( nbondj /= -1 .AND. nbondj /= 2 )   bmask(:,  1 :jprecj) = 0.e0
303         IF( nbondj /=  1 .AND. nbondj /= 2 )   bmask(:,nlcj:jpj   ) = 0.e0
304     
305         ! north fold : bmask must be set to 0. on rows jpj-1 and jpj
306         IF( npolj == 3 .OR. npolj == 4 ) THEN
307            IF( lk_dynspg_rl ) THEN
308               DO ji = 1, nlci
309                  bmask(ji,nlcj-1) = 0.e0
310                  bmask(ji,nlcj  ) = 0.e0
311               END DO
312            ELSE
313               DO ji = 1, nlci
314                  ii = ji + nimpp - 1
315                  bmask(ji,nlcj-1) = bmask(ji,nlcj-1) * tpol(ii)
316                  bmask(ji,nlcj  ) = 0.e0
317               END DO
318            ENDIF
319         ENDIF
320         IF( npolj == 5 .OR. npolj == 6 ) THEN
321            IF( lk_dynspg_rl ) THEN
322               DO ji = 1, nlci
323                  ii = ji + nimpp - 1
324                  bmask(ji,nlcj-1) = bmask(ji,nlcj-1) * fpol(ii)
325                  bmask(ji,nlcj  ) = 0.e0
326               END DO
327            ELSE
328               DO ji = 1, nlci
329                  bmask(ji,nlcj  ) = 0.e0
330               END DO
331            ENDIF
332         ENDIF
333      ENDIF
334
335
336      ! mask for second order calculation of vorticity
337      ! ----------------------------------------------
338     
339      CALL dom_msk_nsa
340
341     
342      ! Lateral boundary conditions on velocity (modify fmask)
343      ! ---------------------------------------
344     
345      DO jk = 1, jpk
346
347         zwf(:,:) = fmask(:,:,jk)
348         
349         DO jj = 2, jpjm1
350            DO ji = fs_2, fs_jpim1   ! vector opt.
351               IF( fmask(ji,jj,jk) == 0. ) THEN
352                  fmask(ji,jj,jk) = shlat * MIN( 1., MAX( zwf(ji+1,jj), zwf(ji,jj+1),   &
353                     &                                    zwf(ji-1,jj), zwf(ji,jj-1)  )  )
354               ENDIF
355            END DO
356         END DO
357         
358         DO jj = 2, jpjm1
359            IF( fmask(1,jj,jk) == 0. ) THEN
360               fmask(1  ,jj,jk) = shlat * MIN( 1., MAX( zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1) ) )
361            ENDIF
362            IF( fmask(jpi,jj,jk) == 0. ) THEN
363               fmask(jpi,jj,jk) = shlat * MIN( 1., MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) )
364            ENDIF
365         END DO
366         
367         DO ji = 2, jpim1
368            IF( fmask(ji,1,jk) == 0. ) THEN
369               fmask(ji, 1 ,jk) = shlat * MIN( 1., MAX( zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1) ) )
370            ENDIF
371            IF( fmask(ji,jpj,jk) == 0. ) THEN
372               fmask(ji,jpj,jk) = shlat * MIN( 1., MAX( zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpjm1) ) )
373            ENDIF
374         END DO
375      END DO
376     
377
378      IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN
379         !                                           ! =======================
380         ! Increased lateral friction in             !  ORCA_R2 configuration
381         ! the vicinity of some straits              ! =======================
382         !
383         IF( n_cla == 0 ) THEN
384            !                                ! Gibraltar strait  : partial slip (fmask=0.5)
385            ij0 = 101   ;   ij1 = 101
386            ii0 = 139   ;   ii1 = 140   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  0.5e0
387            ij0 = 102   ;   ij1 = 102
388            ii0 = 139   ;   ii1 = 140   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  0.5e0
389            !
390            !                                ! Bab el Mandeb : partial slip (fmask=1)
391            ij0 =  87   ;   ij1 =  88
392            ii0 = 160   ;   ii1 = 160   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  1.e0
393            ij0 =  88   ;   ij1 =  88
394            ii0 = 159   ;   ii1 = 159   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  1.e0
395            !
396         ENDIF
397         !                                ! Sound  strait  : strong slip (fmask > 2)
398         ij0 = 115   ;   ij1 = 115
399         ii0 = 145   ;   ii1 = 146   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 4.0e0
400         ij0 = 116   ;   ij1 = 116
401         ii0 = 145   ;   ii1 = 146   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 4.0e0
402         !
403      ENDIF
404     
405      ! Lateral boundary conditions on fmask
406      CALL lbc_lnk( fmask, 'F', 1. )
407     
408      ! Mbathy set to the number of w-level (minimum value 2)
409      ! -----------------------------------
410      IF( lk_isl ) THEN
411         ! this is done at the end of solver_init routine
412      ELSE
413         DO jj = 1, jpj
414            DO ji = 1, jpi
415               mbathy(ji,jj) = MAX( 1, mbathy(ji,jj) ) + 1
416            END DO
417         END DO
418      ENDIF
419     
420      ! Control print
421      ! -------------
422      IF( nprint == 1 .AND. lwp ) THEN
423         imsk(:,:) = INT( tmask_i(:,:) )
424         WRITE(numout,*) ' tmask_i : '
425         CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1,   &
426               &                           1, jpj, 1, 1, numout)
427         WRITE (numout,*)
428         WRITE (numout,*) ' dommsk: tmask for each level'
429         WRITE (numout,*) ' ----------------------------'
430         DO jk = 1, jpk
431            imsk(:,:) = INT( tmask(:,:,jk) )
432
433            WRITE(numout,*)
434            WRITE(numout,*) ' level = ',jk
435            CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1,   &
436               &                              1, jpj, 1, 1, numout)
437         END DO
438         WRITE(numout,*)
439         WRITE(numout,*) ' dom_msk: vmask for each level'
440         WRITE(numout,*) ' -----------------------------'
441         DO jk = 1, jpk
442            imsk(:,:) = INT( vmask(:,:,jk) )
443            WRITE(numout,*)
444            WRITE(numout,*) ' level = ',jk
445            CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1,   &
446               &                              1, jpj, 1, 1, numout)
447         END DO
448         WRITE(numout,*)
449         WRITE(numout,*) ' dom_msk: fmask for each level'
450         WRITE(numout,*) ' -----------------------------'
451         DO jk = 1, jpk
452            imsk(:,:) = INT( fmask(:,:,jk) )
453            WRITE(numout,*)
454            WRITE(numout,*) ' level = ',jk
455            CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1,   &
456               &                              1, jpj, 1, 1, numout )
457         END DO
458         WRITE(numout,*)
459         WRITE(numout,*) ' dom_msk: bmask '
460         WRITE(numout,*) ' ---------------'
461         WRITE(numout,*)
462         imsk(:,:) = INT( bmask(:,:) )
463         CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1,   &
464               &                           1, jpj, 1, 1, numout )
465      ENDIF
466
467   END SUBROUTINE dom_msk
468
469#if defined key_noslip_accurate
470   !!----------------------------------------------------------------------
471   !!   'key_noslip_accurate' :         accurate no-slip boundary condition
472   !!----------------------------------------------------------------------
473   
474   SUBROUTINE dom_msk_nsa
475      !!---------------------------------------------------------------------
476      !!                 ***  ROUTINE dom_msk_nsa  ***
477      !!
478      !! ** Purpose :
479      !!
480      !! ** Method  :
481      !!
482      !! ** Action :
483      !!
484      !! History :
485      !!        !  00-03  (G. Madec)  no slip accurate
486      !!----------------------------------------------------------------------
487      INTEGER  :: ji, jj, jk, jl      ! dummy loop indices
488      INTEGER ::   ine, inw, ins, inn, itest, ierror, iind, ijnd
489      INTEGER, DIMENSION(jpi*jpj*jpk,3) ::  icoord
490      REAL(wp) ::   zaa
491      !!---------------------------------------------------------------------
492     
493
494      IF(lwp)WRITE(numout,*)
495      IF(lwp)WRITE(numout,*) 'dom_msk_nsa : noslip accurate boundary condition'
496      IF(lwp)WRITE(numout,*) '~~~~~~~~~~~   using Schchepetkin and O Brian scheme'
497      IF( lk_mpp ) CALL ctl_stop( ' mpp version is not yet implemented' )
498
499      ! mask for second order calculation of vorticity
500      ! ----------------------------------------------
501      ! noslip boundary condition: fmask=1  at convex corner, store
502      ! index of straight coast meshes ( 'west', refering to a coast,
503      ! means west of the ocean, aso)
504     
505      DO jk = 1, jpk
506         DO jl = 1, 4
507            npcoa(jl,jk) = 0
508            DO ji = 1, 2*(jpi+jpj)
509               nicoa(ji,jl,jk) = 0
510               njcoa(ji,jl,jk) = 0
511            END DO
512         END DO
513      END DO
514     
515      IF( jperio == 2 ) THEN
516         WRITE(numout,*) ' '
517         WRITE(numout,*) ' symetric boundary conditions need special'
518         WRITE(numout,*) ' treatment not implemented. we stop.'
519         STOP
520      ENDIF
521     
522      ! convex corners
523     
524      DO jk = 1, jpkm1
525         DO jj = 1, jpjm1
526            DO ji = 1, jpim1
527               zaa = tmask(ji  ,jj,jk) + tmask(ji  ,jj+1,jk)   &
528                  &+ tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk)
529               IF( ABS(zaa-3.) <= 0.1 )   fmask(ji,jj,jk) = 1.
530            END DO
531         END DO
532      END DO
533
534      ! north-south straight coast
535
536      DO jk = 1, jpkm1
537         inw = 0
538         ine = 0
539         DO jj = 2, jpjm1
540            DO ji = 2, jpim1
541               zaa = tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk)
542               IF( ABS(zaa-2.) <= 0.1 .AND. fmask(ji,jj,jk) == 0 ) THEN
543                  inw = inw + 1
544                  nicoa(inw,1,jk) = ji
545                  njcoa(inw,1,jk) = jj
546                  IF( nprint == 1 ) WRITE(numout,*) ' west  : ', jk, inw, ji, jj
547               ENDIF
548               zaa = tmask(ji,jj,jk) + tmask(ji,jj+1,jk)
549               IF( ABS(zaa-2.) <= 0.1 .AND. fmask(ji,jj,jk) == 0 ) THEN
550                  ine = ine + 1
551                  nicoa(ine,2,jk) = ji
552                  njcoa(ine,2,jk) = jj
553                  IF( nprint == 1 ) WRITE(numout,*) ' east  : ', jk, ine, ji, jj
554               ENDIF
555            END DO
556         END DO
557         npcoa(1,jk) = inw
558         npcoa(2,jk) = ine
559      END DO
560
561      ! west-east straight coast
562
563      DO jk = 1, jpkm1
564         ins = 0
565         inn = 0
566         DO jj = 2, jpjm1
567            DO ji =2, jpim1
568               zaa = tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk)
569               IF( ABS(zaa-2.) <= 0.1 .AND. fmask(ji,jj,jk) == 0 ) THEN
570                  ins = ins + 1
571                  nicoa(ins,3,jk) = ji
572                  njcoa(ins,3,jk) = jj
573                  IF( nprint == 1 ) WRITE(numout,*) ' south : ', jk, ins, ji, jj
574               ENDIF
575               zaa = tmask(ji+1,jj,jk) + tmask(ji,jj,jk)
576               IF( ABS(zaa-2.) <= 0.1 .AND. fmask(ji,jj,jk) == 0 ) THEN
577                  inn = inn + 1
578                  nicoa(inn,4,jk) = ji
579                  njcoa(inn,4,jk) = jj
580                  IF( nprint == 1 ) WRITE(numout,*) ' north : ', jk, inn, ji, jj
581               ENDIF
582            END DO
583         END DO
584         npcoa(3,jk) = ins
585         npcoa(4,jk) = inn
586      END DO
587
588      itest = 2 * ( jpi + jpj )
589      DO jk = 1, jpk
590         IF( npcoa(1,jk) > itest .OR. npcoa(2,jk) > itest .OR.   &
591             npcoa(3,jk) > itest .OR. npcoa(4,jk) > itest ) THEN
592           
593            WRITE(ctmp1,*) ' level jk = ',jk
594            WRITE(ctmp2,*) ' straight coast index arraies are too small.:'
595            WRITE(ctmp3,*) ' npe, npw, nps, npn = ', npcoa(1,jk), npcoa(2,jk),   &
596                &                                     npcoa(3,jk), npcoa(4,jk)
597            WRITE(ctmp4,*) ' 2*(jpi+jpj) = ',itest,'. we stop.'
598            CALL ctl_stop( ctmp1, ctmp2, ctmp3, ctmp4 )
599        ENDIF
600      END DO
601
602      ierror = 0
603      iind = 0
604      ijnd = 0
605      IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) iind = 2
606      IF( nperio == 3 .OR. nperio == 4 .OR. nperio == 5 .OR. nperio == 6 ) ijnd = 2
607      DO jk = 1, jpk
608         DO jl = 1, npcoa(1,jk)
609            IF( nicoa(jl,1,jk)+3 > jpi+iind ) THEN
610               ierror = ierror+1
611               icoord(ierror,1) = nicoa(jl,1,jk)
612               icoord(ierror,2) = njcoa(jl,1,jk)
613               icoord(ierror,3) = jk
614            ENDIF
615         END DO
616         DO jl = 1, npcoa(2,jk)
617            IF(nicoa(jl,2,jk)-2 < 1-iind ) THEN
618               ierror = ierror + 1
619               icoord(ierror,1) = nicoa(jl,2,jk)
620               icoord(ierror,2) = njcoa(jl,2,jk)
621               icoord(ierror,3) = jk
622            ENDIF
623         END DO
624         DO jl = 1, npcoa(3,jk)
625            IF( njcoa(jl,3,jk)+3 > jpj+ijnd ) THEN
626               ierror = ierror + 1
627               icoord(ierror,1) = nicoa(jl,3,jk)
628               icoord(ierror,2) = njcoa(jl,3,jk)
629               icoord(ierror,3) = jk
630            ENDIF
631         END DO
632         DO jl=1,npcoa(4,jk)
633            IF( njcoa(jl,4,jk)-2 < 1) THEN
634               ierror=ierror+1
635               icoord(ierror,1)=nicoa(jl,4,jk)
636               icoord(ierror,2)=njcoa(jl,4,jk)
637               icoord(ierror,3)=jk
638            ENDIF
639         END DO
640      END DO
641     
642      IF( ierror > 0 ) THEN
643         IF(lwp) WRITE(numout,*)
644         IF(lwp) WRITE(numout,*) '              Problem on lateral conditions'
645         IF(lwp) WRITE(numout,*) '                 Bad marking off at points:'
646         DO jl = 1, ierror
647            IF(lwp) WRITE(numout,*) 'Level:',icoord(jl,3),   &
648               &                  '  Point(',icoord(jl,1),',',icoord(jl,2),')'
649         END DO
650         CALL ctl_stop( 'We stop...' )
651      ENDIF
652
653   END SUBROUTINE dom_msk_nsa
654
655#else
656   !!----------------------------------------------------------------------
657   !!   Default option :                                      Empty routine
658   !!----------------------------------------------------------------------
659   SUBROUTINE dom_msk_nsa       
660   END SUBROUTINE dom_msk_nsa
661#endif
662   
663   !!======================================================================
664END MODULE dommsk
Note: See TracBrowser for help on using the repository browser.