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

Last change on this file since 62 was 62, checked in by opalod, 20 years ago

CT : BUGFIX036 : # Loop indice correction for the North fold in rigid lid and mpp case PLUS

# Use logical key "lk_isl" instead of "l_isl"
# Addition of the logical key lk_dynspg_fsc_atsk in parallel to the lk_dynspg_fsc one
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 26.6 KB
Line 
1MODULE dommsk
2   !!==============================================================================
3   !!                       ***  MODULE dommsk   ***
4   !! Ocean initialization : domain land/sea mask
5   !!==============================================================================
6
7   !!----------------------------------------------------------------------
8   !!   dom_msk        : compute land/ocean mask
9   !!   dom_msk_nsa    : update land/ocean mask when no-slip accurate
10   !!                    option is used.
11   !!----------------------------------------------------------------------
12   !! * Modules used
13   USE oce             ! ocean dynamics and tracers
14   USE dom_oce         ! ocean space and time domain
15   USE obc_oce         ! ocean open boundary conditions
16   USE in_out_manager  ! I/O manager
17   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
18   USE lib_mpp
19   USE 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            !                                ! Gibraltar strait and Gulf of Cadiz
384            ij0 = 101   ;   ij1 = 102
385            ii0 = 137   ;   ii1 = 140   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 14.7e0
386            !                                ! Bab el Mandeb strait
387            ij0 =  87   ;   ij1 =  88
388            ii0 = 162   ;   ii1 = 163   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 20.0e0
389            !                                ! Sound  strait
390            ij0 = 116   ;   ij1 = 117
391            ii0 = 147   ;   ii1 = 148   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 10.0e0
392         ELSE
393            !                                ! Gibraltar strait and Gulf of Cadiz
394            ij0 = 102   ;   ij1 = 102
395            ii0 = 137   ;   ii1 = 139   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  0.0e0
396            ij0 = 101   ;   ij1 = 101
397            ii0 = 139   ;   ii1 = 139   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  0.0e0
398            ij0 = 100   ;   ij1 = 100
399            ii0 = 137   ;   ii1 = 139   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  0.0e0
400            !                                ! Sound  strait
401            ij0 = 116   ;   ij1 = 117
402            ii0 = 147   ;   ii1 = 148   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 10.0e0
403         ENDIF
404         !
405      ENDIF
406     
407      ! Lateral boundary conditions on fmask
408      CALL lbc_lnk( fmask, 'F', 1. )
409     
410      ! Mbathy set to the number of w-level (minimum value 2)
411      ! -----------------------------------
412      IF( lk_isl ) THEN
413         ! this is done at the end of solver_init routine
414      ELSE
415         DO jj = 1, jpj
416            DO ji = 1, jpi
417               mbathy(ji,jj) = MAX( 1, mbathy(ji,jj) ) + 1
418            END DO
419         END DO
420      ENDIF
421     
422      ! Control print
423      ! -------------
424      IF( nprint == 1 .AND. lwp ) THEN
425         imsk(:,:) = INT( tmask_i(:,:) )
426         WRITE(numout,*) ' tmask_i : '
427         CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1,   &
428               &                           1, jpj, 1, 1, numout)
429         WRITE (numout,*)
430         WRITE (numout,*) ' dommsk: tmask for each level'
431         WRITE (numout,*) ' ----------------------------'
432         DO jk = 1, jpk
433            imsk(:,:) = INT( tmask(:,:,jk) )
434
435            WRITE(numout,*)
436            WRITE(numout,*) ' level = ',jk
437            CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1,   &
438               &                              1, jpj, 1, 1, numout)
439         END DO
440         WRITE(numout,*)
441         WRITE(numout,*) ' dom_msk: vmask for each level'
442         WRITE(numout,*) ' -----------------------------'
443         DO jk = 1, jpk
444            imsk(:,:) = INT( vmask(:,:,jk) )
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: fmask for each level'
452         WRITE(numout,*) ' -----------------------------'
453         DO jk = 1, jpk
454            imsk(:,:) = INT( fmask(:,:,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: bmask '
462         WRITE(numout,*) ' ---------------'
463         WRITE(numout,*)
464         imsk(:,:) = INT( bmask(:,:) )
465         CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1,   &
466               &                           1, jpj, 1, 1, numout )
467      ENDIF
468
469   END SUBROUTINE dom_msk
470
471#if defined key_noslip_accurate
472   !!----------------------------------------------------------------------
473   !!   'key_noslip_accurate' :         accurate no-slip boundary condition
474   !!----------------------------------------------------------------------
475   
476   SUBROUTINE dom_msk_nsa
477      !!---------------------------------------------------------------------
478      !!                 ***  ROUTINE dom_msk_nsa  ***
479      !!
480      !! ** Purpose :
481      !!
482      !! ** Method  :
483      !!
484      !! ** Action :
485      !!
486      !! History :
487      !!        !  00-03  (G. Madec)  no slip accurate
488      !!----------------------------------------------------------------------
489      !! *Local declarations
490      INTEGER  :: ji, jj, jk, ii     ! dummy loop indices
491      INTEGER ::   ine, inw, ins, inn, itest, ierror, iind, ijnd
492      INTEGER, DIMENSION(jpi*jpj*jpk,3) ::  icoord
493      REAL(wp) ::   zaa
494      !!---------------------------------------------------------------------
495      !!  OPA 9.0, LODYC-IPSL (2003)
496      !!---------------------------------------------------------------------
497     
498
499      IF(lwp)WRITE(numout,*)
500      IF(lwp)WRITE(numout,*) 'dom_msk_nsa : noslip accurate boundary condition'
501      IF(lwp)WRITE(numout,*) '~~~~~~~~~~~   using Schchepetkin and O Brian scheme'
502      IF( lk_mpp ) THEN
503         IF(lwp)WRITE(numout,cform_err)
504         IF(lwp)WRITE(numout,*) ' mpp version is not yet implemented'
505         nstop = nstop + 1
506      ENDIF
507
508      ! mask for second order calculation of vorticity
509      ! ----------------------------------------------
510      ! noslip boundary condition: fmask=1  at convex corner, store
511      ! index of straight coast meshes ( 'west', refering to a coast,
512      ! means west of the ocean, aso)
513     
514      DO jk = 1, jpk
515         DO jl = 1, 4
516            npcoa(jl,jk) = 0
517            DO ji = 1, 2*(jpi+jpj)
518               nicoa(ji,jl,jk) = 0
519               njcoa(ji,jl,jk) = 0
520            END DO
521         END DO
522      END DO
523     
524      IF( jperio == 2 ) THEN
525         WRITE(numout,*) ' '
526         WRITE(numout,*) ' symetric boundary conditions need special'
527         WRITE(numout,*) ' treatment not implemented. we stop.'
528         STOP
529      ENDIF
530     
531      ! convex corners
532     
533      DO jk = 1, jpkm1
534         DO jj = 1, jpjm1
535            DO ji = 1, jpim1
536               zaa = tmask(ji  ,jj,jk) + tmask(ji  ,jj+1,jk)   &
537                  &+ tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk)
538               IF( ABS(zaa-3.) <= 0.1 )   fmask(ji,jj,jk) = 1.
539            END DO
540         END DO
541      END DO
542
543      ! north-south straight coast
544
545      DO jk = 1, jpkm1
546         inw = 0
547         ine = 0
548         DO jj = 2, jpjm1
549            DO ji = 2, jpim1
550               zaa = tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk)
551               IF( ABS(zaa-2.) <= 0.1 .AND. fmask(ji,jj,jk) == 0 ) THEN
552                  inw = inw + 1
553                  nicoa(inw,1,jk) = ji
554                  njcoa(inw,1,jk) = jj
555                  IF( nprint == 1 ) WRITE(numout,*) ' west  : ', jk, inw, ji, jj
556               ENDIF
557               zaa = tmask(ji,jj,jk) + tmask(ji,jj+1,jk)
558               IF( ABS(zaa-2.) <= 0.1 .AND. fmask(ji,jj,jk) == 0 ) THEN
559                  ine = ine + 1
560                  nicoa(ine,2,jk) = ji
561                  njcoa(ine,2,jk) = jj
562                  IF( nprint == 1 ) WRITE(numout,*) ' east  : ', jk, ine, ji, jj
563               ENDIF
564            END DO
565         END DO
566         npcoa(1,jk) = inw
567         npcoa(2,jk) = ine
568      END DO
569
570      ! west-east straight coast
571
572      DO jk = 1, jpkm1
573         ins = 0
574         inn = 0
575         DO jj = 2, jpjm1
576            DO ji =2, jpim1
577               zaa = tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk)
578               IF( ABS(zaa-2.) <= 0.1 .AND. fmask(ji,jj,jk) == 0 ) THEN
579                  ins = ins + 1
580                  nicoa(ins,3,jk) = ji
581                  njcoa(ins,3,jk) = jj
582                  IF( nprint == 1 ) WRITE(numout,*) ' south : ', jk, ins, ji, jj
583               ENDIF
584               zaa = tmask(ji+1,jj,jk) + tmask(ji,jj,jk)
585               IF( ABS(zaa-2.) <= 0.1 .AND. fmask(ji,jj,jk) == 0 ) THEN
586                  inn = inn + 1
587                  nicoa(inn,4,jk) = ji
588                  njcoa(inn,4,jk) = jj
589                  IF( nprint == 1 ) WRITE(numout,*) ' north : ', jk, inn, ji, jj
590               ENDIF
591            END DO
592         END DO
593         npcoa(3,jk) = ins
594         npcoa(4,jk) = inn
595      END DO
596
597      itest = 2 * ( jpi + jpj )
598      DO jk = 1, jpk
599         IF( npcoa(1,jk) > itest .OR. npcoa(2,jk) > itest .OR.   &
600             npcoa(3,jk) > itest .OR. npcoa(4,jk) > itest ) THEN
601            WRITE(numout,*)
602            WRITE(numout,*) ' level jk = ',jk
603            WRITE(numout,*) ' straight coast index arraies are too small.:'
604            WRITE(numout,*) ' npe, npw, nps, npn = ', npcoa(1,jk), npcoa(2,jk),   &
605                &                                     npcoa(3,jk), npcoa(4,jk)
606            WRITE(numout,*) ' 2*(jpi+jpj) = ',itest,'. we stop.'
607            STOP   !!bug nstop to be used
608        ENDIF
609      END DO
610
611      ierror = 0
612      iind = 0
613      ijnd = 0
614      IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) iind = 2
615      IF( nperio == 3 .OR. nperio == 4 .OR. nperio == 5 .OR. nperio == 6 ) ijnd = 2
616      DO jk = 1, jpk
617         DO jl = 1, npcoa(1,jk)
618            IF( nicoa(jl,1,jk)+3 > jpi+iind ) THEN
619               ierror = ierror+1
620               icoord(ierror,1) = nicoa(jl,1,jk)
621               icoord(ierror,2) = njcoa(jl,1,jk)
622               icoord(ierror,3) = jk
623            ENDIF
624         END DO
625         DO jl = 1, npcoa(2,jk)
626            IF(nicoa(jl,2,jk)-2 < 1-iind ) THEN
627               ierror = ierror + 1
628               icoord(ierror,1) = nicoa(jl,2,jk)
629               icoord(ierror,2) = njcoa(jl,2,jk)
630               icoord(ierror,3) = jk
631            ENDIF
632         END DO
633         DO jl = 1, npcoa(3,jk)
634            IF( njcoa(jl,3,jk)+3 > jpj+ijnd ) THEN
635               ierror = ierror + 1
636               icoord(ierror,1) = nicoa(jl,3,jk)
637               icoord(ierror,2) = njcoa(jl,3,jk)
638               icoord(ierror,3) = jk
639            ENDIF
640         END DO
641         DO jl=1,npcoa(4,jk)
642            IF( njcoa(jl,4,jk)-2 < 1) THEN
643               ierror=ierror+1
644               icoord(ierror,1)=nicoa(jl,4,jk)
645               icoord(ierror,2)=njcoa(jl,4,jk)
646               icoord(ierror,3)=jk
647            ENDIF
648         END DO
649      END DO
650     
651      IF( ierror > 0 ) THEN
652         IF(lwp) WRITE(numout,*)
653         IF(lwp) WRITE(numout,*) '              Problem on lateral conditions'
654         IF(lwp) WRITE(numout,*) '                 Bad marking off at points:'
655         DO jl = 1, ierror
656            IF(lwp) WRITE(numout,*) 'Level:',icoord(jl,3),   &
657               &                  '  Point(',icoord(jl,1),',',icoord(jl,2),')'
658         END DO
659         IF(lwp) WRITE(numout,*) 'We stop...'   !!cr print format to be used
660         nstop = nstop + 1
661      ENDIF
662
663   END SUBROUTINE dom_msk_nsa
664
665#else
666   !!----------------------------------------------------------------------
667   !!   Default option :                                      Empty routine
668   !!----------------------------------------------------------------------
669   SUBROUTINE dom_msk_nsa       
670   END SUBROUTINE dom_msk_nsa
671#endif
672   
673   !!======================================================================
674END MODULE dommsk
Note: See TracBrowser for help on using the repository browser.