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

Last change on this file since 719 was 719, checked in by ctlod, 14 years ago

get back to the nemo_v2_3 version for trunk

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 26.7 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   !! $Header$
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               fmask(ji,jj,jk) = tmask(ji,jj  ,jk) * tmask(ji+1,jj  ,jk)   &
232                  &            * tmask(ji,jj+1,jk) * tmask(ji+1,jj+1,jk)
233            END DO
234         END DO
235      END DO
236
237      IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN
238         !                                           ! =======================
239         ! modified vmask value in                   !  ORCA_R2 configuration
240         ! the vicinity of some straits              ! =======================
241
242         IF( n_cla == 1 ) THEN 
243            !                                ! vmask = 0. on Gibraltar zonal section
244            ij0 = 101   ;   ij1 = 101
245            ii0 = 138   ;   ii1 = 139   ;   vmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 19:jpk ) = 0.e0
246            !                                ! vmask = 0. on Bab el Mandeb zonal section
247            ij0 =  87   ;   ij1 =  87
248            ii0 = 161   ;   ii1 = 163   ;   vmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 18:jpk ) = 0.e0
249         ENDIF
250
251      ENDIF
252     
253      ! Lateral boundary conditions
254      CALL lbc_lnk( umask, 'U', 1. )
255      CALL lbc_lnk( vmask, 'V', 1. )
256      CALL lbc_lnk( fmask, 'F', 1. )
257
258
259      ! 4. ocean/land mask for the elliptic equation
260      ! --------------------------------------------
261     
262      ! Computation
263      IF( lk_dynspg_rl ) THEN
264         bmask(:,:) = fmask(:,:,1)       ! elliptic equation is written at f-point
265      ELSE
266         bmask(:,:) = tmask(:,:,1)       ! elliptic equation is written at t-point
267      ENDIF
268     
269      ! Boundary conditions
270      !   cyclic east-west : bmask must be set to 0. on rows 1 and jpi
271      IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN
272         bmask( 1 ,:) = 0.e0
273         bmask(jpi,:) = 0.e0
274      ENDIF
275     
276      !   south symmetric :  bmask must be set to 0. on row 1
277      IF( nperio == 2 ) THEN
278         bmask(:, 1 ) = 0.e0
279      ENDIF
280     
281      !   north fold :
282      IF( nperio == 3 .OR. nperio == 4 ) THEN
283         IF( lk_dynspg_rl ) THEN
284            ! T-pt pivot and F-pt elliptic eq. : bmask set to 0. on rows jpj-1 and jpj
285            bmask(:,jpj-1) = 0.e0
286            bmask(:,jpj  ) = 0.e0
287         ELSE
288            ! T-pt pivot and T-pt elliptic eq. : bmask set to 0. on row jpj and on half jpjglo-1 row
289            DO ji = 1, jpi
290               ii = ji + nimpp - 1
291               bmask(ji,jpj-1) = bmask(ji,jpj-1) * tpol(ii)
292               bmask(ji,jpj  ) = 0.e0
293            END DO
294         ENDIF
295      ENDIF
296      IF( nperio == 5 .OR. nperio == 6 ) THEN
297         IF( lk_dynspg_rl ) THEN
298            ! F-pt pivot and F-pt elliptic eq. : bmask set to 0. on row jpj and on half jpjglo-1 row
299            DO ji = 1, jpi
300               ii = ji + nimpp - 1
301               bmask(ji,jpj-1) = bmask(ji,jpj-1) * fpol(ii)
302               bmask(ji,jpj  ) = 0.e0
303            END DO
304         ELSE
305            ! F-pt pivot and T-pt elliptic eq. : bmask set to 0. on row jpj
306            bmask(:,jpj) = 0.e0
307         ENDIF
308      ENDIF
309
310      ! Mpp boundary conditions: bmask is set to zero on the overlap
311      ! region for all elliptic solvers
312
313      IF( lk_mpp ) THEN
314         IF( nbondi /= -1 .AND. nbondi /= 2 )   bmask(  1 :jpreci,:) = 0.e0
315         IF( nbondi /=  1 .AND. nbondi /= 2 )   bmask(nlci:jpi   ,:) = 0.e0
316         IF( nbondj /= -1 .AND. nbondj /= 2 )   bmask(:,  1 :jprecj) = 0.e0
317         IF( nbondj /=  1 .AND. nbondj /= 2 )   bmask(:,nlcj:jpj   ) = 0.e0
318     
319         ! north fold : bmask must be set to 0. on rows jpj-1 and jpj
320         IF( npolj == 3 .OR. npolj == 4 ) THEN
321            IF( lk_dynspg_rl ) THEN
322               DO ji = 1, nlci
323                  bmask(ji,nlcj-1) = 0.e0
324                  bmask(ji,nlcj  ) = 0.e0
325               END DO
326            ELSE
327               DO ji = 1, nlci
328                  ii = ji + nimpp - 1
329                  bmask(ji,nlcj-1) = bmask(ji,nlcj-1) * tpol(ii)
330                  bmask(ji,nlcj  ) = 0.e0
331               END DO
332            ENDIF
333         ENDIF
334         IF( npolj == 5 .OR. npolj == 6 ) THEN
335            IF( lk_dynspg_rl ) THEN
336               DO ji = 1, nlci
337                  ii = ji + nimpp - 1
338                  bmask(ji,nlcj-1) = bmask(ji,nlcj-1) * fpol(ii)
339                  bmask(ji,nlcj  ) = 0.e0
340               END DO
341            ELSE
342               DO ji = 1, nlci
343                  bmask(ji,nlcj  ) = 0.e0
344               END DO
345            ENDIF
346         ENDIF
347      ENDIF
348
349
350      ! mask for second order calculation of vorticity
351      ! ----------------------------------------------
352     
353      CALL dom_msk_nsa
354
355     
356      ! Lateral boundary conditions on velocity (modify fmask)
357      ! ---------------------------------------
358     
359      DO jk = 1, jpk
360
361         zwf(:,:) = fmask(:,:,jk)
362         
363         DO jj = 2, jpjm1
364            DO ji = fs_2, fs_jpim1   ! vector opt.
365               IF( fmask(ji,jj,jk) == 0. ) THEN
366                  fmask(ji,jj,jk) = shlat * MIN( 1., MAX( zwf(ji+1,jj), zwf(ji,jj+1),   &
367                     &                                    zwf(ji-1,jj), zwf(ji,jj-1)  )  )
368               ENDIF
369            END DO
370         END DO
371         
372         DO jj = 2, jpjm1
373            IF( fmask(1,jj,jk) == 0. ) THEN
374               fmask(1  ,jj,jk) = shlat * MIN( 1., MAX( zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1) ) )
375            ENDIF
376            IF( fmask(jpi,jj,jk) == 0. ) THEN
377               fmask(jpi,jj,jk) = shlat * MIN( 1., MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) )
378            ENDIF
379         END DO
380         
381         DO ji = 2, jpim1
382            IF( fmask(ji,1,jk) == 0. ) THEN
383               fmask(ji, 1 ,jk) = shlat * MIN( 1., MAX( zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1) ) )
384            ENDIF
385            IF( fmask(ji,jpj,jk) == 0. ) THEN
386               fmask(ji,jpj,jk) = shlat * MIN( 1., MAX( zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpjm1) ) )
387            ENDIF
388         END DO
389      END DO
390     
391
392      IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN
393         !                                           ! =======================
394         ! Increased lateral friction in             !  ORCA_R2 configuration
395         ! the vicinity of some straits              ! =======================
396         !
397         IF( n_cla == 0 ) THEN
398            !                                ! Sound  strait
399            ij0 = 116   ;   ij1 = 117
400            ii0 = 147   ;   ii1 = 148   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 10.0e0
401         ELSE
402            !                                ! Gibraltar strait and Gulf of Cadiz
403            ij0 = 102   ;   ij1 = 102
404            ii0 = 137   ;   ii1 = 139   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  0.0e0
405            ij0 = 101   ;   ij1 = 101
406            ii0 = 139   ;   ii1 = 139   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  0.0e0
407            ij0 = 100   ;   ij1 = 100
408            ii0 = 137   ;   ii1 = 139   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  0.0e0
409            !                                ! Sound  strait
410            ij0 = 116   ;   ij1 = 117
411            ii0 = 147   ;   ii1 = 148   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 10.0e0
412         ENDIF
413         !
414      ENDIF
415     
416      ! Lateral boundary conditions on fmask
417      CALL lbc_lnk( fmask, 'F', 1. )
418     
419      ! Mbathy set to the number of w-level (minimum value 2)
420      ! -----------------------------------
421      IF( lk_isl ) THEN
422         ! this is done at the end of solver_init routine
423      ELSE
424         DO jj = 1, jpj
425            DO ji = 1, jpi
426               mbathy(ji,jj) = MAX( 1, mbathy(ji,jj) ) + 1
427            END DO
428         END DO
429      ENDIF
430     
431      ! Control print
432      ! -------------
433      IF( nprint == 1 .AND. lwp ) THEN
434         imsk(:,:) = INT( tmask_i(:,:) )
435         WRITE(numout,*) ' tmask_i : '
436         CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1,   &
437               &                           1, jpj, 1, 1, numout)
438         WRITE (numout,*)
439         WRITE (numout,*) ' dommsk: tmask for each level'
440         WRITE (numout,*) ' ----------------------------'
441         DO jk = 1, jpk
442            imsk(:,:) = INT( tmask(:,:,jk) )
443
444            WRITE(numout,*)
445            WRITE(numout,*) ' level = ',jk
446            CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1,   &
447               &                              1, jpj, 1, 1, numout)
448         END DO
449         WRITE(numout,*)
450         WRITE(numout,*) ' dom_msk: vmask for each level'
451         WRITE(numout,*) ' -----------------------------'
452         DO jk = 1, jpk
453            imsk(:,:) = INT( vmask(:,:,jk) )
454            WRITE(numout,*)
455            WRITE(numout,*) ' level = ',jk
456            CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1,   &
457               &                              1, jpj, 1, 1, numout)
458         END DO
459         WRITE(numout,*)
460         WRITE(numout,*) ' dom_msk: fmask for each level'
461         WRITE(numout,*) ' -----------------------------'
462         DO jk = 1, jpk
463            imsk(:,:) = INT( fmask(:,:,jk) )
464            WRITE(numout,*)
465            WRITE(numout,*) ' level = ',jk
466            CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1,   &
467               &                              1, jpj, 1, 1, numout )
468         END DO
469         WRITE(numout,*)
470         WRITE(numout,*) ' dom_msk: bmask '
471         WRITE(numout,*) ' ---------------'
472         WRITE(numout,*)
473         imsk(:,:) = INT( bmask(:,:) )
474         CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1,   &
475               &                           1, jpj, 1, 1, numout )
476      ENDIF
477
478   END SUBROUTINE dom_msk
479
480#if defined key_noslip_accurate
481   !!----------------------------------------------------------------------
482   !!   'key_noslip_accurate' :         accurate no-slip boundary condition
483   !!----------------------------------------------------------------------
484   
485   SUBROUTINE dom_msk_nsa
486      !!---------------------------------------------------------------------
487      !!                 ***  ROUTINE dom_msk_nsa  ***
488      !!
489      !! ** Purpose :
490      !!
491      !! ** Method  :
492      !!
493      !! ** Action :
494      !!
495      !! History :
496      !!        !  00-03  (G. Madec)  no slip accurate
497      !!----------------------------------------------------------------------
498      INTEGER  :: ji, jj, jk, jl      ! dummy loop indices
499      INTEGER ::   ine, inw, ins, inn, itest, ierror, iind, ijnd, ii
500      INTEGER, DIMENSION(jpi*jpj*jpk,3) ::  icoord
501      REAL(wp) ::   zaa
502      !!---------------------------------------------------------------------
503     
504
505      IF(lwp)WRITE(numout,*)
506      IF(lwp)WRITE(numout,*) 'dom_msk_nsa : noslip accurate boundary condition'
507      IF(lwp)WRITE(numout,*) '~~~~~~~~~~~   using Schchepetkin and O Brian scheme'
508      IF( lk_mpp ) CALL ctl_stop( ' mpp version is not yet implemented' )
509
510      ! mask for second order calculation of vorticity
511      ! ----------------------------------------------
512      ! noslip boundary condition: fmask=1  at convex corner, store
513      ! index of straight coast meshes ( 'west', refering to a coast,
514      ! means west of the ocean, aso)
515     
516      DO jk = 1, jpk
517         DO jl = 1, 4
518            npcoa(jl,jk) = 0
519            DO ji = 1, 2*(jpi+jpj)
520               nicoa(ji,jl,jk) = 0
521               njcoa(ji,jl,jk) = 0
522            END DO
523         END DO
524      END DO
525     
526      IF( jperio == 2 ) THEN
527         WRITE(numout,*) ' '
528         WRITE(numout,*) ' symetric boundary conditions need special'
529         WRITE(numout,*) ' treatment not implemented. we stop.'
530         STOP
531      ENDIF
532     
533      ! convex corners
534     
535      DO jk = 1, jpkm1
536         DO jj = 1, jpjm1
537            DO ji = 1, jpim1
538               zaa = tmask(ji  ,jj,jk) + tmask(ji  ,jj+1,jk)   &
539                  &+ tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk)
540               IF( ABS(zaa-3.) <= 0.1 )   fmask(ji,jj,jk) = 1.
541            END DO
542         END DO
543      END DO
544
545      ! north-south straight coast
546
547      DO jk = 1, jpkm1
548         inw = 0
549         ine = 0
550         DO jj = 2, jpjm1
551            DO ji = 2, jpim1
552               zaa = tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk)
553               IF( ABS(zaa-2.) <= 0.1 .AND. fmask(ji,jj,jk) == 0 ) THEN
554                  inw = inw + 1
555                  nicoa(inw,1,jk) = ji
556                  njcoa(inw,1,jk) = jj
557                  IF( nprint == 1 ) WRITE(numout,*) ' west  : ', jk, inw, ji, jj
558               ENDIF
559               zaa = tmask(ji,jj,jk) + tmask(ji,jj+1,jk)
560               IF( ABS(zaa-2.) <= 0.1 .AND. fmask(ji,jj,jk) == 0 ) THEN
561                  ine = ine + 1
562                  nicoa(ine,2,jk) = ji
563                  njcoa(ine,2,jk) = jj
564                  IF( nprint == 1 ) WRITE(numout,*) ' east  : ', jk, ine, ji, jj
565               ENDIF
566            END DO
567         END DO
568         npcoa(1,jk) = inw
569         npcoa(2,jk) = ine
570      END DO
571
572      ! west-east straight coast
573
574      DO jk = 1, jpkm1
575         ins = 0
576         inn = 0
577         DO jj = 2, jpjm1
578            DO ji =2, jpim1
579               zaa = tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk)
580               IF( ABS(zaa-2.) <= 0.1 .AND. fmask(ji,jj,jk) == 0 ) THEN
581                  ins = ins + 1
582                  nicoa(ins,3,jk) = ji
583                  njcoa(ins,3,jk) = jj
584                  IF( nprint == 1 ) WRITE(numout,*) ' south : ', jk, ins, ji, jj
585               ENDIF
586               zaa = tmask(ji+1,jj,jk) + tmask(ji,jj,jk)
587               IF( ABS(zaa-2.) <= 0.1 .AND. fmask(ji,jj,jk) == 0 ) THEN
588                  inn = inn + 1
589                  nicoa(inn,4,jk) = ji
590                  njcoa(inn,4,jk) = jj
591                  IF( nprint == 1 ) WRITE(numout,*) ' north : ', jk, inn, ji, jj
592               ENDIF
593            END DO
594         END DO
595         npcoa(3,jk) = ins
596         npcoa(4,jk) = inn
597      END DO
598
599      itest = 2 * ( jpi + jpj )
600      DO jk = 1, jpk
601         IF( npcoa(1,jk) > itest .OR. npcoa(2,jk) > itest .OR.   &
602             npcoa(3,jk) > itest .OR. npcoa(4,jk) > itest ) THEN
603           
604            WRITE(ctmp1,*) ' level jk = ',jk
605            WRITE(ctmp2,*) ' straight coast index arraies are too small.:'
606            WRITE(ctmp3,*) ' npe, npw, nps, npn = ', npcoa(1,jk), npcoa(2,jk),   &
607                &                                     npcoa(3,jk), npcoa(4,jk)
608            WRITE(ctmp4,*) ' 2*(jpi+jpj) = ',itest,'. we stop.'
609            CALL ctl_stop( ctmp1, ctmp2, ctmp3, ctmp4 )
610        ENDIF
611      END DO
612
613      ierror = 0
614      iind = 0
615      ijnd = 0
616      IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) iind = 2
617      IF( nperio == 3 .OR. nperio == 4 .OR. nperio == 5 .OR. nperio == 6 ) ijnd = 2
618      DO jk = 1, jpk
619         DO jl = 1, npcoa(1,jk)
620            IF( nicoa(jl,1,jk)+3 > jpi+iind ) THEN
621               ierror = ierror+1
622               icoord(ierror,1) = nicoa(jl,1,jk)
623               icoord(ierror,2) = njcoa(jl,1,jk)
624               icoord(ierror,3) = jk
625            ENDIF
626         END DO
627         DO jl = 1, npcoa(2,jk)
628            IF(nicoa(jl,2,jk)-2 < 1-iind ) THEN
629               ierror = ierror + 1
630               icoord(ierror,1) = nicoa(jl,2,jk)
631               icoord(ierror,2) = njcoa(jl,2,jk)
632               icoord(ierror,3) = jk
633            ENDIF
634         END DO
635         DO jl = 1, npcoa(3,jk)
636            IF( njcoa(jl,3,jk)+3 > jpj+ijnd ) THEN
637               ierror = ierror + 1
638               icoord(ierror,1) = nicoa(jl,3,jk)
639               icoord(ierror,2) = njcoa(jl,3,jk)
640               icoord(ierror,3) = jk
641            ENDIF
642         END DO
643         DO jl=1,npcoa(4,jk)
644            IF( njcoa(jl,4,jk)-2 < 1) THEN
645               ierror=ierror+1
646               icoord(ierror,1)=nicoa(jl,4,jk)
647               icoord(ierror,2)=njcoa(jl,4,jk)
648               icoord(ierror,3)=jk
649            ENDIF
650         END DO
651      END DO
652     
653      IF( ierror > 0 ) THEN
654         IF(lwp) WRITE(numout,*)
655         IF(lwp) WRITE(numout,*) '              Problem on lateral conditions'
656         IF(lwp) WRITE(numout,*) '                 Bad marking off at points:'
657         DO jl = 1, ierror
658            IF(lwp) WRITE(numout,*) 'Level:',icoord(jl,3),   &
659               &                  '  Point(',icoord(jl,1),',',icoord(jl,2),')'
660         END DO
661         CALL ctl_stop( 'We stop...' )
662      ENDIF
663
664   END SUBROUTINE dom_msk_nsa
665
666#else
667   !!----------------------------------------------------------------------
668   !!   Default option :                                      Empty routine
669   !!----------------------------------------------------------------------
670   SUBROUTINE dom_msk_nsa       
671   END SUBROUTINE dom_msk_nsa
672#endif
673   
674   !!======================================================================
675END MODULE dommsk
Note: See TracBrowser for help on using the repository browser.