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

Last change on this file since 1528 was 1528, checked in by rblod, 15 years ago

Suppress rigid-lid option, see ticket #486
  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 24.6 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 dynspg_oce      ! choice/control of key cpp for surface pressure gradient
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 , LOCEAN-IPSL (2005)
35   !! $Id$
36   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
37   !!---------------------------------------------------------------------------------
38
39CONTAINS
40   
41   SUBROUTINE dom_msk
42      !!---------------------------------------------------------------------
43      !!                 ***  ROUTINE dom_msk  ***
44      !!
45      !! ** Purpose :   Compute land/ocean mask arrays at tracer points, hori-
46      !!      zontal velocity points (u & v), vorticity points (f) and baro-
47      !!      tropic stream function  points (b).
48      !!        Set mbathy to the number of non-zero w-levels of a water column
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
53      !!      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 (lk_obc=T):
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      !!   9.0  !  05-11  (V. Garnier) Surface pressure gradient organization
121      !!----------------------------------------------------------------------
122      INTEGER  ::   ji, jj, jk      ! dummy loop indices
123      INTEGER  ::   iif, iil, ii0, ii1, ii
124      INTEGER  ::   ijf, ijl, ij0, ij1
125      INTEGER, DIMENSION(jpi,jpj) ::  imsk
126      REAL(wp), DIMENSION(jpi,jpj) ::   zwf
127
128      NAMELIST/namlbc/ shlat
129      !!---------------------------------------------------------------------
130     
131      ! Namelist namlbc : lateral momentum boundary condition
132      REWIND( numnam )
133      READ  ( numnam, namlbc )
134      IF(lwp) THEN
135         WRITE(numout,*)
136         WRITE(numout,*) 'dommsk : ocean mask '
137         WRITE(numout,*) '~~~~~~'
138         WRITE(numout,*) '         Namelist namlbc'
139         WRITE(numout,*) '            lateral momentum boundary cond. shlat = ',shlat
140      ENDIF
141
142      IF ( shlat == 0. ) THEN
143          IF(lwp) WRITE(numout,*) '         ocean lateral free-slip '
144        ELSEIF ( shlat  ==  2. ) THEN
145          IF(lwp) WRITE(numout,*) '         ocean lateral  no-slip '
146        ELSEIF ( 0. < shlat .AND. shlat < 2. ) THEN
147          IF(lwp) WRITE(numout,*) '         ocean lateral  partial-slip '
148        ELSEIF ( 2. < shlat ) THEN
149          IF(lwp) WRITE(numout,*) '         ocean lateral  strong-slip '
150        ELSE
151          WRITE(ctmp1,*) ' shlat is negative = ', shlat
152          CALL ctl_stop( ctmp1 )
153      ENDIF
154
155      ! 1. Ocean/land mask at t-point (computed from mbathy)
156      ! -----------------------------
157      ! Tmask has already the right boundary conditions since mbathy is ok
158
159      tmask(:,:,:) = 0.e0
160      DO jk = 1, jpk
161         DO jj = 1, jpj
162            DO ji = 1, jpi
163               IF( FLOAT( mbathy(ji,jj)-jk )+.1 >= 0.e0 ) tmask(ji,jj,jk) = 1.e0
164            END DO 
165         END DO 
166      END DO 
167
168#if defined key_zdfkpp
169      IF( cp_cfg == 'orca' )   THEN
170         IF( jp_cfg == 2 )   THEN
171            ! land point on Bab el Mandeb zonal section
172            ij0 =  87   ;   ij1 =  88
173            ii0 = 160   ;   ii1 = 161
174            tmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 0.e0
175         ELSE
176            IF(lwp) WRITE(numout,*)
177            IF(lwp) WRITE(numout,cform_war)
178            IF(lwp) WRITE(numout,*)
179            IF(lwp) WRITE(numout,*)'          A mask must be applied on Bab el Mandeb strait'
180            IF(lwp) WRITE(numout,*)'          in case of ORCAs configurations'
181            IF(lwp) WRITE(numout,*)'          This is a problem which is not yet solved'
182            IF(lwp) WRITE(numout,*)
183         ENDIF
184      ENDIF
185#endif
186
187      ! Interior domain mask (used for global sum)
188      ! --------------------
189
190      tmask_i(:,:) = tmask(:,:,1)
191      iif = jpreci                         ! ???
192      iil = nlci - jpreci + 1
193      ijf = jprecj                         ! ???
194      ijl = nlcj - jprecj + 1
195
196      tmask_i( 1 :iif,   :   ) = 0.e0      ! first columns
197      tmask_i(iil:jpi,   :   ) = 0.e0      ! last  columns (including mpp extra columns)
198      tmask_i(   :   , 1 :ijf) = 0.e0      ! first rows
199      tmask_i(   :   ,ijl:jpj) = 0.e0      ! last  rows (including mpp extra rows)
200
201      ! north fold mask
202      tpol(1:jpiglo) = 1.e0 
203      fpol(1:jpiglo) = 1.e0
204      IF( jperio == 3 .OR. jperio == 4 ) THEN      ! T-point pivot
205         tpol(jpiglo/2+1:jpiglo) = 0.e0
206         fpol(     1    :jpiglo) = 0.e0
207         ! T-point pivot: only half of the nlcj-1 row
208         IF( mjg(nlej) == jpjglo )   THEN
209            DO ji = iif+1, iil-1
210               tmask_i(ji,nlej-1) = tmask_i(ji,nlej-1) * tpol(mig(ji))
211            END DO
212         ENDIF
213      ENDIF
214      IF( jperio == 5 .OR. jperio == 6 ) THEN      ! F-point pivot
215         tpol(     1    :jpiglo) = 0.e0
216         fpol(jpiglo/2+1:jpiglo) = 0.e0
217      ENDIF
218
219      ! 2. Ocean/land mask at u-,  v-, and z-points (computed from tmask)
220      ! -------------------------------------------
221     
222      ! Computation
223      DO jk = 1, jpk
224         DO jj = 1, jpjm1
225            DO ji = 1, fs_jpim1   ! vector loop
226               umask(ji,jj,jk) = tmask(ji,jj  ,jk) * tmask(ji+1,jj  ,jk)
227               vmask(ji,jj,jk) = tmask(ji,jj  ,jk) * tmask(ji  ,jj+1,jk)
228            END DO
229            DO ji = 1, jpim1 
230               fmask(ji,jj,jk) = tmask(ji,jj  ,jk) * tmask(ji+1,jj  ,jk)   &
231                  &            * tmask(ji,jj+1,jk) * tmask(ji+1,jj+1,jk)
232            END DO
233         END DO
234      END DO
235
236      ! Lateral boundary conditions
237      CALL lbc_lnk( umask, 'U', 1. )
238      CALL lbc_lnk( vmask, 'V', 1. )
239      CALL lbc_lnk( fmask, 'F', 1. )
240
241
242      ! 4. ocean/land mask for the elliptic equation
243      ! --------------------------------------------
244     
245      ! Computation
246      bmask(:,:) = tmask(:,:,1)       ! elliptic equation is written at t-point
247     
248      ! Boundary conditions
249      !   cyclic east-west : bmask must be set to 0. on rows 1 and jpi
250      IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN
251         bmask( 1 ,:) = 0.e0
252         bmask(jpi,:) = 0.e0
253      ENDIF
254     
255      !   south symmetric :  bmask must be set to 0. on row 1
256      IF( nperio == 2 ) THEN
257         bmask(:, 1 ) = 0.e0
258      ENDIF
259     
260      !   north fold :
261      IF( nperio == 3 .OR. nperio == 4 ) THEN
262         ! T-pt pivot and T-pt elliptic eq. : bmask set to 0. on row jpj and on half jpjglo-1 row
263         DO ji = 1, jpi
264            ii = ji + nimpp - 1
265            bmask(ji,jpj-1) = bmask(ji,jpj-1) * tpol(ii)
266            bmask(ji,jpj  ) = 0.e0
267         END DO
268      ENDIF
269      IF( nperio == 5 .OR. nperio == 6 ) THEN
270         ! F-pt pivot and T-pt elliptic eq. : bmask set to 0. on row jpj
271         bmask(:,jpj) = 0.e0
272      ENDIF
273
274      ! Mpp boundary conditions: bmask is set to zero on the overlap
275      ! region for all elliptic solvers
276
277      IF( lk_mpp ) THEN
278         IF( nbondi /= -1 .AND. nbondi /= 2 )   bmask(  1 :jpreci,:) = 0.e0
279         IF( nbondi /=  1 .AND. nbondi /= 2 )   bmask(nlci:jpi   ,:) = 0.e0
280         IF( nbondj /= -1 .AND. nbondj /= 2 )   bmask(:,  1 :jprecj) = 0.e0
281         IF( nbondj /=  1 .AND. nbondj /= 2 )   bmask(:,nlcj:jpj   ) = 0.e0
282     
283         ! north fold : bmask must be set to 0. on rows jpj-1 and jpj
284         IF( npolj == 3 .OR. npolj == 4 ) THEN
285            DO ji = 1, nlci
286               ii = ji + nimpp - 1
287               bmask(ji,nlcj-1) = bmask(ji,nlcj-1) * tpol(ii)
288               bmask(ji,nlcj  ) = 0.e0
289            END DO
290         ENDIF
291         IF( npolj == 5 .OR. npolj == 6 ) THEN
292            DO ji = 1, nlci
293               bmask(ji,nlcj  ) = 0.e0
294            END DO
295         ENDIF
296      ENDIF
297
298
299      ! mask for second order calculation of vorticity
300      ! ----------------------------------------------
301     
302      CALL dom_msk_nsa
303
304     
305      ! Lateral boundary conditions on velocity (modify fmask)
306      ! ---------------------------------------
307     
308      DO jk = 1, jpk
309
310         zwf(:,:) = fmask(:,:,jk)
311         
312         DO jj = 2, jpjm1
313            DO ji = fs_2, fs_jpim1   ! vector opt.
314               IF( fmask(ji,jj,jk) == 0. ) THEN
315                  fmask(ji,jj,jk) = shlat * MIN( 1., MAX( zwf(ji+1,jj), zwf(ji,jj+1),   &
316                     &                                    zwf(ji-1,jj), zwf(ji,jj-1)  )  )
317               ENDIF
318            END DO
319         END DO
320         
321         DO jj = 2, jpjm1
322            IF( fmask(1,jj,jk) == 0. ) THEN
323               fmask(1  ,jj,jk) = shlat * MIN( 1., MAX( zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1) ) )
324            ENDIF
325            IF( fmask(jpi,jj,jk) == 0. ) THEN
326               fmask(jpi,jj,jk) = shlat * MIN( 1., MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) )
327            ENDIF
328         END DO
329         
330         DO ji = 2, jpim1
331            IF( fmask(ji,1,jk) == 0. ) THEN
332               fmask(ji, 1 ,jk) = shlat * MIN( 1., MAX( zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1) ) )
333            ENDIF
334            IF( fmask(ji,jpj,jk) == 0. ) THEN
335               fmask(ji,jpj,jk) = shlat * MIN( 1., MAX( zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpjm1) ) )
336            ENDIF
337         END DO
338      END DO
339     
340
341      IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN
342         !                                           ! =======================
343         ! Increased lateral friction in             !  ORCA_R2 configuration
344         ! the vicinity of some straits              ! =======================
345         !
346         IF( n_cla == 0 ) THEN
347            !                                ! Gibraltar strait  : partial slip (fmask=0.5)
348            ij0 = 101   ;   ij1 = 101
349            ii0 = 139   ;   ii1 = 140   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  0.5e0
350            ij0 = 102   ;   ij1 = 102
351            ii0 = 139   ;   ii1 = 140   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  0.5e0
352            !
353            !                                ! Bab el Mandeb : partial slip (fmask=1)
354            ij0 =  87   ;   ij1 =  88
355            ii0 = 160   ;   ii1 = 160   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  1.e0
356            ij0 =  88   ;   ij1 =  88
357            ii0 = 159   ;   ii1 = 159   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  1.e0
358            !
359         ENDIF
360         !                                ! Sound  strait  : strong slip (fmask > 2)
361         ij0 = 115   ;   ij1 = 115
362         ii0 = 145   ;   ii1 = 146   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 4.0e0
363         ij0 = 116   ;   ij1 = 116
364         ii0 = 145   ;   ii1 = 146   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 4.0e0
365         !
366      ENDIF
367     
368      ! Lateral boundary conditions on fmask
369      CALL lbc_lnk( fmask, 'F', 1. )
370     
371      ! Mbathy set to the number of w-level (minimum value 2)
372      ! -----------------------------------
373      DO jj = 1, jpj
374         DO ji = 1, jpi
375            mbathy(ji,jj) = MAX( 1, mbathy(ji,jj) ) + 1
376         END DO
377      END DO
378     
379      ! Control print
380      ! -------------
381      IF( nprint == 1 .AND. lwp ) THEN
382         imsk(:,:) = INT( tmask_i(:,:) )
383         WRITE(numout,*) ' tmask_i : '
384         CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1,   &
385               &                           1, jpj, 1, 1, numout)
386         WRITE (numout,*)
387         WRITE (numout,*) ' dommsk: tmask for each level'
388         WRITE (numout,*) ' ----------------------------'
389         DO jk = 1, jpk
390            imsk(:,:) = INT( tmask(:,:,jk) )
391
392            WRITE(numout,*)
393            WRITE(numout,*) ' level = ',jk
394            CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1,   &
395               &                              1, jpj, 1, 1, numout)
396         END DO
397         WRITE(numout,*)
398         WRITE(numout,*) ' dom_msk: vmask for each level'
399         WRITE(numout,*) ' -----------------------------'
400         DO jk = 1, jpk
401            imsk(:,:) = INT( vmask(:,:,jk) )
402            WRITE(numout,*)
403            WRITE(numout,*) ' level = ',jk
404            CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1,   &
405               &                              1, jpj, 1, 1, numout)
406         END DO
407         WRITE(numout,*)
408         WRITE(numout,*) ' dom_msk: fmask for each level'
409         WRITE(numout,*) ' -----------------------------'
410         DO jk = 1, jpk
411            imsk(:,:) = INT( fmask(:,:,jk) )
412            WRITE(numout,*)
413            WRITE(numout,*) ' level = ',jk
414            CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1,   &
415               &                              1, jpj, 1, 1, numout )
416         END DO
417         WRITE(numout,*)
418         WRITE(numout,*) ' dom_msk: bmask '
419         WRITE(numout,*) ' ---------------'
420         WRITE(numout,*)
421         imsk(:,:) = INT( bmask(:,:) )
422         CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1,   &
423               &                           1, jpj, 1, 1, numout )
424      ENDIF
425
426   END SUBROUTINE dom_msk
427
428#if defined key_noslip_accurate
429   !!----------------------------------------------------------------------
430   !!   'key_noslip_accurate' :         accurate no-slip boundary condition
431   !!----------------------------------------------------------------------
432   
433   SUBROUTINE dom_msk_nsa
434      !!---------------------------------------------------------------------
435      !!                 ***  ROUTINE dom_msk_nsa  ***
436      !!
437      !! ** Purpose :
438      !!
439      !! ** Method  :
440      !!
441      !! ** Action :
442      !!
443      !! History :
444      !!        !  00-03  (G. Madec)  no slip accurate
445      !!----------------------------------------------------------------------
446      INTEGER  :: ji, jj, jk, jl      ! dummy loop indices
447      INTEGER ::   ine, inw, ins, inn, itest, ierror, iind, ijnd
448      INTEGER, DIMENSION(jpi*jpj*jpk,3) ::  icoord
449      REAL(wp) ::   zaa
450      !!---------------------------------------------------------------------
451     
452
453      IF(lwp)WRITE(numout,*)
454      IF(lwp)WRITE(numout,*) 'dom_msk_nsa : noslip accurate boundary condition'
455      IF(lwp)WRITE(numout,*) '~~~~~~~~~~~   using Schchepetkin and O Brian scheme'
456      IF( lk_mpp ) CALL ctl_stop( ' mpp version is not yet implemented' )
457
458      ! mask for second order calculation of vorticity
459      ! ----------------------------------------------
460      ! noslip boundary condition: fmask=1  at convex corner, store
461      ! index of straight coast meshes ( 'west', refering to a coast,
462      ! means west of the ocean, aso)
463     
464      DO jk = 1, jpk
465         DO jl = 1, 4
466            npcoa(jl,jk) = 0
467            DO ji = 1, 2*(jpi+jpj)
468               nicoa(ji,jl,jk) = 0
469               njcoa(ji,jl,jk) = 0
470            END DO
471         END DO
472      END DO
473     
474      IF( jperio == 2 ) THEN
475         WRITE(numout,*) ' '
476         WRITE(numout,*) ' symetric boundary conditions need special'
477         WRITE(numout,*) ' treatment not implemented. we stop.'
478         STOP
479      ENDIF
480     
481      ! convex corners
482     
483      DO jk = 1, jpkm1
484         DO jj = 1, jpjm1
485            DO ji = 1, jpim1
486               zaa = tmask(ji  ,jj,jk) + tmask(ji  ,jj+1,jk)   &
487                  &+ tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk)
488               IF( ABS(zaa-3.) <= 0.1 )   fmask(ji,jj,jk) = 1.
489            END DO
490         END DO
491      END DO
492
493      ! north-south straight coast
494
495      DO jk = 1, jpkm1
496         inw = 0
497         ine = 0
498         DO jj = 2, jpjm1
499            DO ji = 2, jpim1
500               zaa = tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk)
501               IF( ABS(zaa-2.) <= 0.1 .AND. fmask(ji,jj,jk) == 0 ) THEN
502                  inw = inw + 1
503                  nicoa(inw,1,jk) = ji
504                  njcoa(inw,1,jk) = jj
505                  IF( nprint == 1 ) WRITE(numout,*) ' west  : ', jk, inw, ji, jj
506               ENDIF
507               zaa = tmask(ji,jj,jk) + tmask(ji,jj+1,jk)
508               IF( ABS(zaa-2.) <= 0.1 .AND. fmask(ji,jj,jk) == 0 ) THEN
509                  ine = ine + 1
510                  nicoa(ine,2,jk) = ji
511                  njcoa(ine,2,jk) = jj
512                  IF( nprint == 1 ) WRITE(numout,*) ' east  : ', jk, ine, ji, jj
513               ENDIF
514            END DO
515         END DO
516         npcoa(1,jk) = inw
517         npcoa(2,jk) = ine
518      END DO
519
520      ! west-east straight coast
521
522      DO jk = 1, jpkm1
523         ins = 0
524         inn = 0
525         DO jj = 2, jpjm1
526            DO ji =2, jpim1
527               zaa = tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk)
528               IF( ABS(zaa-2.) <= 0.1 .AND. fmask(ji,jj,jk) == 0 ) THEN
529                  ins = ins + 1
530                  nicoa(ins,3,jk) = ji
531                  njcoa(ins,3,jk) = jj
532                  IF( nprint == 1 ) WRITE(numout,*) ' south : ', jk, ins, ji, jj
533               ENDIF
534               zaa = tmask(ji+1,jj,jk) + tmask(ji,jj,jk)
535               IF( ABS(zaa-2.) <= 0.1 .AND. fmask(ji,jj,jk) == 0 ) THEN
536                  inn = inn + 1
537                  nicoa(inn,4,jk) = ji
538                  njcoa(inn,4,jk) = jj
539                  IF( nprint == 1 ) WRITE(numout,*) ' north : ', jk, inn, ji, jj
540               ENDIF
541            END DO
542         END DO
543         npcoa(3,jk) = ins
544         npcoa(4,jk) = inn
545      END DO
546
547      itest = 2 * ( jpi + jpj )
548      DO jk = 1, jpk
549         IF( npcoa(1,jk) > itest .OR. npcoa(2,jk) > itest .OR.   &
550             npcoa(3,jk) > itest .OR. npcoa(4,jk) > itest ) THEN
551           
552            WRITE(ctmp1,*) ' level jk = ',jk
553            WRITE(ctmp2,*) ' straight coast index arraies are too small.:'
554            WRITE(ctmp3,*) ' npe, npw, nps, npn = ', npcoa(1,jk), npcoa(2,jk),   &
555                &                                     npcoa(3,jk), npcoa(4,jk)
556            WRITE(ctmp4,*) ' 2*(jpi+jpj) = ',itest,'. we stop.'
557            CALL ctl_stop( ctmp1, ctmp2, ctmp3, ctmp4 )
558        ENDIF
559      END DO
560
561      ierror = 0
562      iind = 0
563      ijnd = 0
564      IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) iind = 2
565      IF( nperio == 3 .OR. nperio == 4 .OR. nperio == 5 .OR. nperio == 6 ) ijnd = 2
566      DO jk = 1, jpk
567         DO jl = 1, npcoa(1,jk)
568            IF( nicoa(jl,1,jk)+3 > jpi+iind ) THEN
569               ierror = ierror+1
570               icoord(ierror,1) = nicoa(jl,1,jk)
571               icoord(ierror,2) = njcoa(jl,1,jk)
572               icoord(ierror,3) = jk
573            ENDIF
574         END DO
575         DO jl = 1, npcoa(2,jk)
576            IF(nicoa(jl,2,jk)-2 < 1-iind ) THEN
577               ierror = ierror + 1
578               icoord(ierror,1) = nicoa(jl,2,jk)
579               icoord(ierror,2) = njcoa(jl,2,jk)
580               icoord(ierror,3) = jk
581            ENDIF
582         END DO
583         DO jl = 1, npcoa(3,jk)
584            IF( njcoa(jl,3,jk)+3 > jpj+ijnd ) THEN
585               ierror = ierror + 1
586               icoord(ierror,1) = nicoa(jl,3,jk)
587               icoord(ierror,2) = njcoa(jl,3,jk)
588               icoord(ierror,3) = jk
589            ENDIF
590         END DO
591         DO jl=1,npcoa(4,jk)
592            IF( njcoa(jl,4,jk)-2 < 1) THEN
593               ierror=ierror+1
594               icoord(ierror,1)=nicoa(jl,4,jk)
595               icoord(ierror,2)=njcoa(jl,4,jk)
596               icoord(ierror,3)=jk
597            ENDIF
598         END DO
599      END DO
600     
601      IF( ierror > 0 ) THEN
602         IF(lwp) WRITE(numout,*)
603         IF(lwp) WRITE(numout,*) '              Problem on lateral conditions'
604         IF(lwp) WRITE(numout,*) '                 Bad marking off at points:'
605         DO jl = 1, ierror
606            IF(lwp) WRITE(numout,*) 'Level:',icoord(jl,3),   &
607               &                  '  Point(',icoord(jl,1),',',icoord(jl,2),')'
608         END DO
609         CALL ctl_stop( 'We stop...' )
610      ENDIF
611
612   END SUBROUTINE dom_msk_nsa
613
614#else
615   !!----------------------------------------------------------------------
616   !!   Default option :                                      Empty routine
617   !!----------------------------------------------------------------------
618   SUBROUTINE dom_msk_nsa       
619   END SUBROUTINE dom_msk_nsa
620#endif
621   
622   !!======================================================================
623END MODULE dommsk
Note: See TracBrowser for help on using the repository browser.