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

Last change on this file since 18 was 3, checked in by opalod, 20 years ago

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