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

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

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