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/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

source: trunk/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90 @ 2528

Last change on this file since 2528 was 2528, checked in by rblod, 13 years ago

Update NEMOGCM from branch nemo_v3_3_beta

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