source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90 @ 2392

Last change on this file since 2392 was 2392, checked in by gm, 10 years ago

v3.3beta: Cross Land Advection (ticket #127) full rewriting + MPP bug corrections

  • Property svn:keywords set to Id
File size: 26.8 KB
Line 
1MODULE dommsk
2   !!======================================================================
3   !!                       ***  MODULE dommsk   ***
4   !! Ocean initialization : domain land/sea mask
5   !!======================================================================
6   !! History :  OPA  ! 1987-07  (G. Madec)  Original code
7   !!            6.0  ! 1993-03  (M. Guyon)  symetrical conditions (M. Guyon)
8   !!            7.0  ! 1996-01  (G. Madec)  suppression of common work arrays
9   !!             -   ! 1996-05  (G. Madec)  mask computed from tmask and sup-
10   !!                 !                      pression of the double computation of bmask
11   !!            8.0  ! 1997-02  (G. Madec)  mesh information put in domhgr.F
12   !!            8.1  ! 1997-07  (G. Madec)  modification of mbathy and fmask
13   !!             -   ! 1998-05  (G. Roullet)  free surface
14   !!            8.2  ! 2000-03  (G. Madec)  no slip accurate
15   !!             -   ! 2001-09  (J.-M. Molines)  Open boundaries
16   !!   NEMO     1.0  ! 2002-08  (G. Madec)  F90: Free form and module
17   !!             -   ! 2005-11  (V. Garnier) Surface pressure gradient organization
18   !!            3.2  ! 2009-07  (R. Benshila) Suppression of rigid-lid option
19   !!----------------------------------------------------------------------
20
21   !!----------------------------------------------------------------------
22   !!   dom_msk        : compute land/ocean mask
23   !!   dom_msk_nsa    : update land/ocean mask when no-slip accurate option is used.
24   !!----------------------------------------------------------------------
25   USE oce             ! ocean dynamics and tracers
26   USE dom_oce         ! ocean space and time domain
27   USE obc_oce         ! ocean open boundary conditions
28   USE in_out_manager  ! I/O manager
29   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
30   USE lib_mpp
31   USE dynspg_oce      ! choice/control of key cpp for surface pressure gradient
32
33   IMPLICIT NONE
34   PRIVATE
35
36   PUBLIC   dom_msk    ! routine called by inidom.F90
37
38   !                            !!* Namelist namlbc : lateral boundary condition *
39   REAL(wp) ::   rn_shlat = 2.   ! type of lateral boundary condition on velocity
40   
41   !! * Substitutions
42#  include "vectopt_loop_substitute.h90"
43   !!----------------------------------------------------------------------
44   !! NEMO/OPA 3.2 , LODYC-IPSL  (2009)
45   !! $Id$
46   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
47   !!----------------------------------------------------------------------
48CONTAINS
49   
50   SUBROUTINE dom_msk
51      !!---------------------------------------------------------------------
52      !!                 ***  ROUTINE dom_msk  ***
53      !!
54      !! ** Purpose :   Compute land/ocean mask arrays at tracer points, hori-
55      !!      zontal velocity points (u & v), vorticity points (f) and baro-
56      !!      tropic stream function  points (b).
57      !!        Set mbathy to the number of non-zero w-levels of a water column
58      !!
59      !! ** Method  :   The ocean/land mask is computed from the basin bathy-
60      !!      metry in level (mbathy) which is defined or read in dommba.
61      !!      mbathy equals 0 over continental T-point
62      !!      and the number of ocean level over the ocean.
63      !!
64      !!      At a given position (ji,jj,jk) the ocean/land mask is given by:
65      !!      t-point : 0. IF mbathy( ji ,jj) =< 0
66      !!                1. IF mbathy( ji ,jj) >= jk
67      !!      u-point : 0. IF mbathy( ji ,jj)  or mbathy(ji+1, jj ) =< 0
68      !!                1. IF mbathy( ji ,jj) and mbathy(ji+1, jj ) >= jk.
69      !!      v-point : 0. IF mbathy( ji ,jj)  or mbathy( ji ,jj+1) =< 0
70      !!                1. IF mbathy( ji ,jj) and mbathy( ji ,jj+1) >= jk.
71      !!      f-point : 0. IF mbathy( ji ,jj)  or mbathy( ji ,jj+1)
72      !!                   or mbathy(ji+1,jj)  or mbathy(ji+1,jj+1) =< 0
73      !!                1. IF mbathy( ji ,jj) and mbathy( ji ,jj+1)
74      !!                and mbathy(ji+1,jj) and mbathy(ji+1,jj+1) >= jk.
75      !!      b-point : the same definition as for f-point of the first ocean
76      !!                level (surface level) but with 0 along coastlines.
77      !!
78      !!        The lateral friction is set through the value of fmask along
79      !!      the coast and topography. This value is defined by rn_shlat, a
80      !!      namelist parameter:
81      !!         rn_shlat = 0, free slip  (no shear along the coast)
82      !!         rn_shlat = 2, no slip  (specified zero velocity at the coast)
83      !!         0 < rn_shlat < 2, partial slip   | non-linear velocity profile
84      !!         2 < rn_shlat, strong slip        | in the lateral boundary layer
85      !!
86      !!      N.B. If nperio not equal to 0, the land/ocean mask arrays
87      !!      are defined with the proper value at lateral domain boundaries,
88      !!      but bmask. indeed, bmask defined the domain over which the
89      !!      barotropic stream function is computed. this domain cannot
90      !!      contain identical columns because the matrix associated with
91      !!      the barotropic stream function equation is then no more inverti-
92      !!      ble. therefore bmask is set to 0 along lateral domain boundaries
93      !!      even IF nperio is not zero.
94      !!
95      !!      In case of open boundaries (lk_obc=T):
96      !!        - tmask is set to 1 on the points to be computed bay the open
97      !!          boundaries routines.
98      !!        - bmask is  set to 0 on the open boundaries.
99      !!
100      !!      Set mbathy to the number of non-zero w-levels of a water column
101      !!                  mbathy = min( mbathy, 1 ) + 1
102      !!      (note that the minimum value of mbathy is 2).
103      !!
104      !! ** Action :   tmask    : land/ocean mask at t-point (=0. or 1.)
105      !!               umask    : land/ocean mask at u-point (=0. or 1.)
106      !!               vmask    : land/ocean mask at v-point (=0. or 1.)
107      !!               fmask    : land/ocean mask at f-point (=0. or 1.)
108      !!                          =rn_shlat along lateral boundaries
109      !!               bmask    : land/ocean mask at barotropic stream
110      !!                          function point (=0. or 1.) and set to 0 along lateral boundaries
111      !!               mbathy   : number of non-zero w-levels
112      !!----------------------------------------------------------------------
113      INTEGER  ::   ji, jj, jk      ! dummy loop indices
114      INTEGER  ::   iif, iil, ii0, ii1, ii
115      INTEGER  ::   ijf, ijl, ij0, ij1
116      INTEGER , DIMENSION(jpi,jpj) ::  imsk
117      REAL(wp), DIMENSION(jpi,jpj) ::   zwf
118      !!
119      NAMELIST/namlbc/ rn_shlat
120      !!---------------------------------------------------------------------
121     
122      REWIND( numnam )              ! Namelist namlbc : lateral momentum boundary condition
123      READ  ( numnam, namlbc )
124     
125      IF(lwp) THEN                  ! control print
126         WRITE(numout,*)
127         WRITE(numout,*) 'dommsk : ocean mask '
128         WRITE(numout,*) '~~~~~~'
129         WRITE(numout,*) '   Namelist namlbc'
130         WRITE(numout,*) '      lateral momentum boundary cond.    rn_shlat = ',rn_shlat
131      ENDIF
132
133      IF     (      rn_shlat == 0.               ) THEN   ;   IF(lwp) WRITE(numout,*) '   ocean lateral  free-slip '
134      ELSEIF (      rn_shlat == 2.               ) THEN   ;   IF(lwp) WRITE(numout,*) '   ocean lateral  no-slip '
135      ELSEIF ( 0. < rn_shlat .AND. rn_shlat < 2. ) THEN   ;   IF(lwp) WRITE(numout,*) '   ocean lateral  partial-slip '
136      ELSEIF ( 2. < rn_shlat                     ) THEN   ;   IF(lwp) WRITE(numout,*) '   ocean lateral  strong-slip '
137      ELSE
138         WRITE(ctmp1,*) ' rn_shlat is negative = ', rn_shlat
139         CALL ctl_stop( ctmp1 )
140      ENDIF
141
142      ! 1. Ocean/land mask at t-point (computed from mbathy)
143      ! -----------------------------
144      ! N.B. tmask has already the right boundary conditions since mbathy is ok
145      !
146      tmask(:,:,:) = 0.e0
147      DO jk = 1, jpk
148         DO jj = 1, jpj
149            DO ji = 1, jpi
150               IF( REAL( mbathy(ji,jj) - jk ) +.1 >= 0.e0 )   tmask(ji,jj,jk) = 1.e0
151            END DO 
152         END DO 
153      END DO 
154
155!!gm  ????
156#if defined key_zdfkpp
157      IF( cp_cfg == 'orca' ) THEN
158         IF( jp_cfg == 2 )   THEN       ! land point on Bab el Mandeb zonal section
159            ij0 =  87   ;   ij1 =  88
160            ii0 = 160   ;   ii1 = 161
161            tmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 0.e0
162         ELSE
163            IF(lwp) WRITE(numout,*)
164            IF(lwp) WRITE(numout,cform_war)
165            IF(lwp) WRITE(numout,*)
166            IF(lwp) WRITE(numout,*)'          A mask must be applied on Bab el Mandeb strait'
167            IF(lwp) WRITE(numout,*)'          in case of ORCAs configurations'
168            IF(lwp) WRITE(numout,*)'          This is a problem which is not yet solved'
169            IF(lwp) WRITE(numout,*)
170         ENDIF
171      ENDIF
172#endif
173!!gm end
174
175      ! Interior domain mask (used for global sum)
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      ! ---------------
190      tpol(1:jpiglo) = 1.e0 
191      fpol(1:jpiglo) = 1.e0
192      IF( jperio == 3 .OR. jperio == 4 ) THEN      ! T-point pivot
193         tpol(jpiglo/2+1:jpiglo) = 0.e0
194         fpol(     1    :jpiglo) = 0.e0
195         IF( mjg(nlej) == jpjglo ) THEN                  ! only half of the nlcj-1 row
196            DO ji = iif+1, iil-1
197               tmask_i(ji,nlej-1) = tmask_i(ji,nlej-1) * tpol(mig(ji))
198            END DO
199         ENDIF
200      ENDIF
201      IF( jperio == 5 .OR. jperio == 6 ) THEN      ! F-point pivot
202         tpol(     1    :jpiglo) = 0.e0
203         fpol(jpiglo/2+1:jpiglo) = 0.e0
204      ENDIF
205
206      ! 2. Ocean/land mask at u-,  v-, and z-points (computed from tmask)
207      ! -------------------------------------------
208      DO jk = 1, jpk
209         DO jj = 1, jpjm1
210            DO ji = 1, fs_jpim1   ! vector loop
211               umask(ji,jj,jk) = tmask(ji,jj  ,jk) * tmask(ji+1,jj  ,jk)
212               vmask(ji,jj,jk) = tmask(ji,jj  ,jk) * tmask(ji  ,jj+1,jk)
213            END DO
214            DO ji = 1, jpim1      ! NO vector opt.
215               fmask(ji,jj,jk) = tmask(ji,jj  ,jk) * tmask(ji+1,jj  ,jk)   &
216                  &            * tmask(ji,jj+1,jk) * tmask(ji+1,jj+1,jk)
217            END DO
218         END DO
219      END DO
220      CALL lbc_lnk( umask, 'U', 1. )      ! Lateral boundary conditions
221      CALL lbc_lnk( vmask, 'V', 1. )
222      CALL lbc_lnk( fmask, 'F', 1. )
223
224
225      ! 4. ocean/land mask for the elliptic equation
226      ! --------------------------------------------
227      bmask(:,:) = tmask(:,:,1)       ! elliptic equation is written at t-point
228      !
229      !                               ! Boundary conditions
230      !                                    ! cyclic east-west : bmask must be set to 0. on rows 1 and jpi
231      IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN
232         bmask( 1 ,:) = 0.e0
233         bmask(jpi,:) = 0.e0
234      ENDIF
235      IF( nperio == 2 ) THEN               ! south symmetric :  bmask must be set to 0. on row 1
236         bmask(:, 1 ) = 0.e0
237      ENDIF
238      !                                    ! north fold :
239      IF( nperio == 3 .OR. nperio == 4 ) THEN   ! T-pt pivot : bmask set to 0. on row jpj and on half jpjglo-1 row
240         DO ji = 1, jpi                     
241            ii = ji + nimpp - 1
242            bmask(ji,jpj-1) = bmask(ji,jpj-1) * tpol(ii)
243            bmask(ji,jpj  ) = 0.e0
244         END DO
245      ENDIF
246      IF( nperio == 5 .OR. nperio == 6 ) THEN   ! F-pt pivot and T-pt elliptic eq. : bmask set to 0. on row jpj
247         bmask(:,jpj) = 0.e0
248      ENDIF
249      !
250      IF( lk_mpp ) THEN                    ! mpp specificities
251         !                                      ! bmask is set to zero on the overlap region
252         IF( nbondi /= -1 .AND. nbondi /= 2 )   bmask(  1 :jpreci,:) = 0.e0
253         IF( nbondi /=  1 .AND. nbondi /= 2 )   bmask(nlci:jpi   ,:) = 0.e0
254         IF( nbondj /= -1 .AND. nbondj /= 2 )   bmask(:,  1 :jprecj) = 0.e0
255         IF( nbondj /=  1 .AND. nbondj /= 2 )   bmask(:,nlcj:jpj   ) = 0.e0
256         !
257         IF( npolj == 3 .OR. npolj == 4 ) THEN  ! north fold : bmask must be set to 0. on rows jpj-1 and jpj
258            DO ji = 1, nlci
259               ii = ji + nimpp - 1
260               bmask(ji,nlcj-1) = bmask(ji,nlcj-1) * tpol(ii)
261               bmask(ji,nlcj  ) = 0.e0
262            END DO
263         ENDIF
264         IF( npolj == 5 .OR. npolj == 6 ) THEN  ! F-pt pivot and T-pt elliptic eq. : bmask set to 0. on row jpj
265            DO ji = 1, nlci
266               bmask(ji,nlcj  ) = 0.e0
267            END DO
268         ENDIF
269      ENDIF
270
271
272      ! mask for second order calculation of vorticity
273      ! ----------------------------------------------
274      CALL dom_msk_nsa
275
276     
277      ! Lateral boundary conditions on velocity (modify fmask)
278      ! ---------------------------------------     
279      DO jk = 1, jpk
280         zwf(:,:) = fmask(:,:,jk)         
281         DO jj = 2, jpjm1
282            DO ji = fs_2, fs_jpim1   ! vector opt.
283               IF( fmask(ji,jj,jk) == 0. ) THEN
284                  fmask(ji,jj,jk) = rn_shlat * MIN( 1., MAX( zwf(ji+1,jj), zwf(ji,jj+1),   &
285                     &                                       zwf(ji-1,jj), zwf(ji,jj-1)  )  )
286               ENDIF
287            END DO
288         END DO
289         DO jj = 2, jpjm1
290            IF( fmask(1,jj,jk) == 0. ) THEN
291               fmask(1  ,jj,jk) = rn_shlat * MIN( 1., MAX( zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1) ) )
292            ENDIF
293            IF( fmask(jpi,jj,jk) == 0. ) THEN
294               fmask(jpi,jj,jk) = rn_shlat * MIN( 1., MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) )
295            ENDIF
296         END DO         
297         DO ji = 2, jpim1
298            IF( fmask(ji,1,jk) == 0. ) THEN
299               fmask(ji, 1 ,jk) = rn_shlat * MIN( 1., MAX( zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1) ) )
300            ENDIF
301            IF( fmask(ji,jpj,jk) == 0. ) THEN
302               fmask(ji,jpj,jk) = rn_shlat * MIN( 1., MAX( zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpjm1) ) )
303            ENDIF
304         END DO
305      END DO
306      !
307      IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN   ! ORCA_R2 configuration
308         !                                                 ! Increased lateral friction near of some straits
309         IF( nn_cla == 0 ) THEN
310            !                                ! Gibraltar strait  : partial slip (fmask=0.5)
311            ij0 = 101   ;   ij1 = 101
312            ii0 = 139   ;   ii1 = 140   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  0.5e0
313            ij0 = 102   ;   ij1 = 102
314            ii0 = 139   ;   ii1 = 140   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  0.5e0
315            !
316            !                                ! Bab el Mandeb : partial slip (fmask=1)
317            ij0 =  87   ;   ij1 =  88
318            ii0 = 160   ;   ii1 = 160   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  1.e0
319            ij0 =  88   ;   ij1 =  88
320            ii0 = 159   ;   ii1 = 159   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  1.e0
321            !
322         ENDIF
323         !                                ! Danish straits  : strong slip (fmask > 2)
324! We keep this as an example but it is instable in this case
325!         ij0 = 115   ;   ij1 = 115
326!         ii0 = 145   ;   ii1 = 146   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 4.0e0
327!         ij0 = 116   ;   ij1 = 116
328!         ii0 = 145   ;   ii1 = 146   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 4.0e0
329         !
330      ENDIF
331      !
332      IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN   ! ORCA R1 configuration
333         !                                                 ! Increased lateral friction near of some straits
334         IF(lwp) WRITE(numout,*)
335         IF(lwp) WRITE(numout,*) '   orca_r1: increase friction near the following straits : '
336         IF(lwp) WRITE(numout,*) '      Gibraltar '
337         ii0 = 283   ;   ii1 = 284        ! Gibraltar Strait
338         ij0 = 200   ;   ij1 = 200   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) =  2.0 
339
340         IF(lwp) WRITE(numout,*) '      Bhosporus '
341         ii0 = 314   ;   ii1 = 315        ! Bhosporus Strait
342         ij0 = 208   ;   ij1 = 208   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) =  2.0 
343
344         IF(lwp) WRITE(numout,*) '      Makassar (Top) '
345         ii0 =  48   ;   ii1 =  48        ! Makassar Strait (Top)
346         ij0 = 149   ;   ij1 = 150   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) =  3.0 
347
348         IF(lwp) WRITE(numout,*) '      Lombok '
349         ii0 =  44   ;   ii1 =  44        ! Lombok Strait
350         ij0 = 124   ;   ij1 = 125   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) =  2.0 
351
352         IF(lwp) WRITE(numout,*) '      Ombai '
353         ii0 =  53   ;   ii1 =  53        ! Ombai Strait
354         ij0 = 124   ;   ij1 = 125   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) = 2.0 
355
356         IF(lwp) WRITE(numout,*) '      Timor Passage '
357         ii0 =  56   ;   ii1 =  56        ! Timor Passage
358         ij0 = 124   ;   ij1 = 125   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) = 2.0 
359
360         IF(lwp) WRITE(numout,*) '      West Halmahera '
361         ii0 =  58   ;   ii1 =  58        ! West Halmahera Strait
362         ij0 = 141   ;   ij1 = 142   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) = 3.0 
363
364         IF(lwp) WRITE(numout,*) '      East Halmahera '
365         ii0 =  55   ;   ii1 =  55        ! East Halmahera Strait
366         ij0 = 141   ;   ij1 = 142   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) = 3.0 
367         !
368      ENDIF
369      !
370      CALL lbc_lnk( fmask, 'F', 1. )      ! Lateral boundary conditions on fmask
371
372     
373      ! Mbathy set to the number of w-level (minimum value 2)
374      ! -----------------------------------
375      DO jj = 1, jpj
376         DO ji = 1, jpi
377            mbathy(ji,jj) = MAX( 1, mbathy(ji,jj) ) + 1
378         END DO
379      END DO
380     
381      IF( nprint == 1 .AND. lwp ) THEN      ! Control print
382         imsk(:,:) = INT( tmask_i(:,:) )
383         WRITE(numout,*) ' tmask_i : '
384         CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1,   &
385               &                           1, jpj, 1, 1, numout)
386         WRITE (numout,*)
387         WRITE (numout,*) ' dommsk: tmask for each level'
388         WRITE (numout,*) ' ----------------------------'
389         DO jk = 1, jpk
390            imsk(:,:) = INT( tmask(:,:,jk) )
391
392            WRITE(numout,*)
393            WRITE(numout,*) ' level = ',jk
394            CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1,   &
395               &                              1, jpj, 1, 1, numout)
396         END DO
397         WRITE(numout,*)
398         WRITE(numout,*) ' dom_msk: vmask for each level'
399         WRITE(numout,*) ' -----------------------------'
400         DO jk = 1, jpk
401            imsk(:,:) = INT( vmask(:,:,jk) )
402            WRITE(numout,*)
403            WRITE(numout,*) ' level = ',jk
404            CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1,   &
405               &                              1, jpj, 1, 1, numout)
406         END DO
407         WRITE(numout,*)
408         WRITE(numout,*) ' dom_msk: fmask for each level'
409         WRITE(numout,*) ' -----------------------------'
410         DO jk = 1, jpk
411            imsk(:,:) = INT( fmask(:,:,jk) )
412            WRITE(numout,*)
413            WRITE(numout,*) ' level = ',jk
414            CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1,   &
415               &                              1, jpj, 1, 1, numout )
416         END DO
417         WRITE(numout,*)
418         WRITE(numout,*) ' dom_msk: bmask '
419         WRITE(numout,*) ' ---------------'
420         WRITE(numout,*)
421         imsk(:,:) = INT( bmask(:,:) )
422         CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1,   &
423               &                           1, jpj, 1, 1, numout )
424      ENDIF
425      !
426   END SUBROUTINE dom_msk
427
428#if defined key_noslip_accurate
429   !!----------------------------------------------------------------------
430   !!   'key_noslip_accurate' :         accurate no-slip boundary condition
431   !!----------------------------------------------------------------------
432   
433   SUBROUTINE dom_msk_nsa
434      !!---------------------------------------------------------------------
435      !!                 ***  ROUTINE dom_msk_nsa  ***
436      !!
437      !! ** Purpose :
438      !!
439      !! ** Method  :
440      !!
441      !! ** Action :
442      !!
443      !!----------------------------------------------------------------------
444      INTEGER  :: ji, jj, jk, jl      ! dummy loop indices
445      INTEGER  ::   ine, inw, ins, inn, itest, ierror, iind, ijnd
446      REAL(wp) ::   zaa
447      INTEGER, DIMENSION(jpi*jpj*jpk,3) ::  icoord
448      !!---------------------------------------------------------------------
449     
450
451      IF(lwp)WRITE(numout,*)
452      IF(lwp)WRITE(numout,*) 'dom_msk_nsa : noslip accurate boundary condition'
453      IF(lwp)WRITE(numout,*) '~~~~~~~~~~~   using Schchepetkin and O Brian scheme'
454      IF( lk_mpp ) CALL ctl_stop( ' mpp version is not yet implemented' )
455
456      ! mask for second order calculation of vorticity
457      ! ----------------------------------------------
458      ! noslip boundary condition: fmask=1  at convex corner, store
459      ! index of straight coast meshes ( 'west', refering to a coast,
460      ! means west of the ocean, aso)
461     
462      DO jk = 1, jpk
463         DO jl = 1, 4
464            npcoa(jl,jk) = 0
465            DO ji = 1, 2*(jpi+jpj)
466               nicoa(ji,jl,jk) = 0
467               njcoa(ji,jl,jk) = 0
468            END DO
469         END DO
470      END DO
471     
472      IF( jperio == 2 ) THEN
473         WRITE(numout,*) ' '
474         WRITE(numout,*) ' symetric boundary conditions need special'
475         WRITE(numout,*) ' treatment not implemented. we stop.'
476         STOP
477      ENDIF
478     
479      ! convex corners
480     
481      DO jk = 1, jpkm1
482         DO jj = 1, jpjm1
483            DO ji = 1, jpim1
484               zaa = tmask(ji  ,jj,jk) + tmask(ji  ,jj+1,jk)   &
485                  &+ tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk)
486               IF( ABS(zaa-3.) <= 0.1 )   fmask(ji,jj,jk) = 1.
487            END DO
488         END DO
489      END DO
490
491      ! north-south straight coast
492
493      DO jk = 1, jpkm1
494         inw = 0
495         ine = 0
496         DO jj = 2, jpjm1
497            DO ji = 2, jpim1
498               zaa = tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk)
499               IF( ABS(zaa-2.) <= 0.1 .AND. fmask(ji,jj,jk) == 0 ) THEN
500                  inw = inw + 1
501                  nicoa(inw,1,jk) = ji
502                  njcoa(inw,1,jk) = jj
503                  IF( nprint == 1 ) WRITE(numout,*) ' west  : ', jk, inw, ji, jj
504               ENDIF
505               zaa = tmask(ji,jj,jk) + tmask(ji,jj+1,jk)
506               IF( ABS(zaa-2.) <= 0.1 .AND. fmask(ji,jj,jk) == 0 ) THEN
507                  ine = ine + 1
508                  nicoa(ine,2,jk) = ji
509                  njcoa(ine,2,jk) = jj
510                  IF( nprint == 1 ) WRITE(numout,*) ' east  : ', jk, ine, ji, jj
511               ENDIF
512            END DO
513         END DO
514         npcoa(1,jk) = inw
515         npcoa(2,jk) = ine
516      END DO
517
518      ! west-east straight coast
519
520      DO jk = 1, jpkm1
521         ins = 0
522         inn = 0
523         DO jj = 2, jpjm1
524            DO ji =2, jpim1
525               zaa = tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk)
526               IF( ABS(zaa-2.) <= 0.1 .AND. fmask(ji,jj,jk) == 0 ) THEN
527                  ins = ins + 1
528                  nicoa(ins,3,jk) = ji
529                  njcoa(ins,3,jk) = jj
530                  IF( nprint == 1 ) WRITE(numout,*) ' south : ', jk, ins, ji, jj
531               ENDIF
532               zaa = tmask(ji+1,jj,jk) + tmask(ji,jj,jk)
533               IF( ABS(zaa-2.) <= 0.1 .AND. fmask(ji,jj,jk) == 0 ) THEN
534                  inn = inn + 1
535                  nicoa(inn,4,jk) = ji
536                  njcoa(inn,4,jk) = jj
537                  IF( nprint == 1 ) WRITE(numout,*) ' north : ', jk, inn, ji, jj
538               ENDIF
539            END DO
540         END DO
541         npcoa(3,jk) = ins
542         npcoa(4,jk) = inn
543      END DO
544
545      itest = 2 * ( jpi + jpj )
546      DO jk = 1, jpk
547         IF( npcoa(1,jk) > itest .OR. npcoa(2,jk) > itest .OR.   &
548             npcoa(3,jk) > itest .OR. npcoa(4,jk) > itest ) THEN
549           
550            WRITE(ctmp1,*) ' level jk = ',jk
551            WRITE(ctmp2,*) ' straight coast index arraies are too small.:'
552            WRITE(ctmp3,*) ' npe, npw, nps, npn = ', npcoa(1,jk), npcoa(2,jk),   &
553                &                                     npcoa(3,jk), npcoa(4,jk)
554            WRITE(ctmp4,*) ' 2*(jpi+jpj) = ',itest,'. we stop.'
555            CALL ctl_stop( ctmp1, ctmp2, ctmp3, ctmp4 )
556        ENDIF
557      END DO
558
559      ierror = 0
560      iind = 0
561      ijnd = 0
562      IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) iind = 2
563      IF( nperio == 3 .OR. nperio == 4 .OR. nperio == 5 .OR. nperio == 6 ) ijnd = 2
564      DO jk = 1, jpk
565         DO jl = 1, npcoa(1,jk)
566            IF( nicoa(jl,1,jk)+3 > jpi+iind ) THEN
567               ierror = ierror+1
568               icoord(ierror,1) = nicoa(jl,1,jk)
569               icoord(ierror,2) = njcoa(jl,1,jk)
570               icoord(ierror,3) = jk
571            ENDIF
572         END DO
573         DO jl = 1, npcoa(2,jk)
574            IF(nicoa(jl,2,jk)-2 < 1-iind ) THEN
575               ierror = ierror + 1
576               icoord(ierror,1) = nicoa(jl,2,jk)
577               icoord(ierror,2) = njcoa(jl,2,jk)
578               icoord(ierror,3) = jk
579            ENDIF
580         END DO
581         DO jl = 1, npcoa(3,jk)
582            IF( njcoa(jl,3,jk)+3 > jpj+ijnd ) THEN
583               ierror = ierror + 1
584               icoord(ierror,1) = nicoa(jl,3,jk)
585               icoord(ierror,2) = njcoa(jl,3,jk)
586               icoord(ierror,3) = jk
587            ENDIF
588         END DO
589         DO jl=1,npcoa(4,jk)
590            IF( njcoa(jl,4,jk)-2 < 1) THEN
591               ierror=ierror+1
592               icoord(ierror,1)=nicoa(jl,4,jk)
593               icoord(ierror,2)=njcoa(jl,4,jk)
594               icoord(ierror,3)=jk
595            ENDIF
596         END DO
597      END DO
598     
599      IF( ierror > 0 ) THEN
600         IF(lwp) WRITE(numout,*)
601         IF(lwp) WRITE(numout,*) '              Problem on lateral conditions'
602         IF(lwp) WRITE(numout,*) '                 Bad marking off at points:'
603         DO jl = 1, ierror
604            IF(lwp) WRITE(numout,*) 'Level:',icoord(jl,3),   &
605               &                  '  Point(',icoord(jl,1),',',icoord(jl,2),')'
606         END DO
607         CALL ctl_stop( 'We stop...' )
608      ENDIF
609
610   END SUBROUTINE dom_msk_nsa
611
612#else
613   !!----------------------------------------------------------------------
614   !!   Default option :                                      Empty routine
615   !!----------------------------------------------------------------------
616   SUBROUTINE dom_msk_nsa       
617   END SUBROUTINE dom_msk_nsa
618#endif
619   
620   !!======================================================================
621END MODULE dommsk
Note: See TracBrowser for help on using the repository browser.