source: branches/2011/dev_LOCEAN_CMCC_2011/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90 @ 3095

Last change on this file since 3095 was 3095, checked in by cetlod, 10 years ago

dev_LOCEAN_CMCC_2011:Move the fmask modification from dommsk to dynvor

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