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

Last change on this file since 3294 was 3294, checked in by rblod, 9 years ago

Merge of 3.4beta into the trunk

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