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

source: branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90 @ 2633

Last change on this file since 2633 was 2633, checked in by trackstand2, 13 years ago

Renamed wrk_use => wrk_in_use and wrk_release => wrk_not_released

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