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

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

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