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

Last change on this file since 467 was 454, checked in by opalod, 18 years ago

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