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

Last change on this file since 236 was 236, checked in by opalod, 19 years ago

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