New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
dommsk.F90 in trunk/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

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

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

First attempt to put dynamic allocation on the trunk

  • 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   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           
393      IF( nprint == 1 .AND. lwp ) THEN      ! Control print
394         imsk(:,:) = INT( tmask_i(:,:) )
395         WRITE(numout,*) ' tmask_i : '
396         CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1,   &
397               &                           1, jpj, 1, 1, numout)
398         WRITE (numout,*)
399         WRITE (numout,*) ' dommsk: tmask for each level'
400         WRITE (numout,*) ' ----------------------------'
401         DO jk = 1, jpk
402            imsk(:,:) = INT( tmask(:,:,jk) )
403
404            WRITE(numout,*)
405            WRITE(numout,*) ' level = ',jk
406            CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1,   &
407               &                              1, jpj, 1, 1, numout)
408         END DO
409         WRITE(numout,*)
410         WRITE(numout,*) ' dom_msk: vmask for each level'
411         WRITE(numout,*) ' -----------------------------'
412         DO jk = 1, jpk
413            imsk(:,:) = INT( vmask(:,:,jk) )
414            WRITE(numout,*)
415            WRITE(numout,*) ' level = ',jk
416            CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1,   &
417               &                              1, jpj, 1, 1, numout)
418         END DO
419         WRITE(numout,*)
420         WRITE(numout,*) ' dom_msk: fmask for each level'
421         WRITE(numout,*) ' -----------------------------'
422         DO jk = 1, jpk
423            imsk(:,:) = INT( fmask(:,:,jk) )
424            WRITE(numout,*)
425            WRITE(numout,*) ' level = ',jk
426            CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1,   &
427               &                              1, jpj, 1, 1, numout )
428         END DO
429         WRITE(numout,*)
430         WRITE(numout,*) ' dom_msk: bmask '
431         WRITE(numout,*) ' ---------------'
432         WRITE(numout,*)
433         imsk(:,:) = INT( bmask(:,:) )
434         CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1,   &
435            &                              1, jpj, 1, 1, numout )
436      ENDIF
437      !
438      IF( wrk_not_released(2, 1)  .OR.   &
439         iwrk_not_released(2, 1)  )   CALL ctl_stop('dom_msk: failed to release workspace arrays')
440      !
441   END SUBROUTINE dom_msk
442
443#if defined key_noslip_accurate
444   !!----------------------------------------------------------------------
445   !!   'key_noslip_accurate' :         accurate no-slip boundary condition
446   !!----------------------------------------------------------------------
447   
448   SUBROUTINE dom_msk_nsa
449      !!---------------------------------------------------------------------
450      !!                 ***  ROUTINE dom_msk_nsa  ***
451      !!
452      !! ** Purpose :
453      !!
454      !! ** Method  :
455      !!
456      !! ** Action :
457      !!----------------------------------------------------------------------
458      INTEGER  ::   ji, jj, jk, jl      ! dummy loop indices
459      INTEGER  ::   ine, inw, ins, inn, itest, ierror, iind, ijnd
460      REAL(wp) ::   zaa
461      !!---------------------------------------------------------------------
462
463      IF(lwp) WRITE(numout,*)
464      IF(lwp) WRITE(numout,*) 'dom_msk_nsa : noslip accurate boundary condition'
465      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   using Schchepetkin and O Brian scheme'
466      IF( lk_mpp )   CALL ctl_stop( ' mpp version is not yet implemented' )
467
468      ! mask for second order calculation of vorticity
469      ! ----------------------------------------------
470      ! noslip boundary condition: fmask=1  at convex corner, store
471      ! index of straight coast meshes ( 'west', refering to a coast,
472      ! means west of the ocean, aso)
473     
474      DO jk = 1, jpk
475         DO jl = 1, 4
476            npcoa(jl,jk) = 0
477            DO ji = 1, 2*(jpi+jpj)
478               nicoa(ji,jl,jk) = 0
479               njcoa(ji,jl,jk) = 0
480            END DO
481         END DO
482      END DO
483     
484      IF( jperio == 2 ) THEN
485         WRITE(numout,*) ' '
486         WRITE(numout,*) ' symetric boundary conditions need special'
487         WRITE(numout,*) ' treatment not implemented. we stop.'
488         STOP
489      ENDIF
490     
491      ! convex corners
492     
493      DO jk = 1, jpkm1
494         DO jj = 1, jpjm1
495            DO ji = 1, jpim1
496               zaa = tmask(ji  ,jj,jk) + tmask(ji  ,jj+1,jk)   &
497                  &+ tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk)
498               IF( ABS(zaa-3._wp) <= 0.1_wp )   fmask(ji,jj,jk) = 1._wp
499            END DO
500         END DO
501      END DO
502
503      ! north-south straight coast
504
505      DO jk = 1, jpkm1
506         inw = 0
507         ine = 0
508         DO jj = 2, jpjm1
509            DO ji = 2, jpim1
510               zaa = tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk)
511               IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN
512                  inw = inw + 1
513                  nicoa(inw,1,jk) = ji
514                  njcoa(inw,1,jk) = jj
515                  IF( nprint == 1 ) WRITE(numout,*) ' west  : ', jk, inw, ji, jj
516               ENDIF
517               zaa = tmask(ji,jj,jk) + tmask(ji,jj+1,jk)
518               IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN
519                  ine = ine + 1
520                  nicoa(ine,2,jk) = ji
521                  njcoa(ine,2,jk) = jj
522                  IF( nprint == 1 ) WRITE(numout,*) ' east  : ', jk, ine, ji, jj
523               ENDIF
524            END DO
525         END DO
526         npcoa(1,jk) = inw
527         npcoa(2,jk) = ine
528      END DO
529
530      ! west-east straight coast
531
532      DO jk = 1, jpkm1
533         ins = 0
534         inn = 0
535         DO jj = 2, jpjm1
536            DO ji =2, jpim1
537               zaa = tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk)
538               IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN
539                  ins = ins + 1
540                  nicoa(ins,3,jk) = ji
541                  njcoa(ins,3,jk) = jj
542                  IF( nprint == 1 ) WRITE(numout,*) ' south : ', jk, ins, ji, jj
543               ENDIF
544               zaa = tmask(ji+1,jj,jk) + tmask(ji,jj,jk)
545               IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN
546                  inn = inn + 1
547                  nicoa(inn,4,jk) = ji
548                  njcoa(inn,4,jk) = jj
549                  IF( nprint == 1 ) WRITE(numout,*) ' north : ', jk, inn, ji, jj
550               ENDIF
551            END DO
552         END DO
553         npcoa(3,jk) = ins
554         npcoa(4,jk) = inn
555      END DO
556
557      itest = 2 * ( jpi + jpj )
558      DO jk = 1, jpk
559         IF( npcoa(1,jk) > itest .OR. npcoa(2,jk) > itest .OR.   &
560             npcoa(3,jk) > itest .OR. npcoa(4,jk) > itest ) THEN
561           
562            WRITE(ctmp1,*) ' level jk = ',jk
563            WRITE(ctmp2,*) ' straight coast index arraies are too small.:'
564            WRITE(ctmp3,*) ' npe, npw, nps, npn = ', npcoa(1,jk), npcoa(2,jk),   &
565                &                                     npcoa(3,jk), npcoa(4,jk)
566            WRITE(ctmp4,*) ' 2*(jpi+jpj) = ',itest,'. we stop.'
567            CALL ctl_stop( ctmp1, ctmp2, ctmp3, ctmp4 )
568        ENDIF
569      END DO
570
571      ierror = 0
572      iind = 0
573      ijnd = 0
574      IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 )   iind = 2
575      IF( nperio == 3 .OR. nperio == 4 .OR. nperio == 5 .OR. nperio == 6 )   ijnd = 2
576      DO jk = 1, jpk
577         DO jl = 1, npcoa(1,jk)
578            IF( nicoa(jl,1,jk)+3 > jpi+iind ) THEN
579               ierror = ierror+1
580               icoord(ierror,1) = nicoa(jl,1,jk)
581               icoord(ierror,2) = njcoa(jl,1,jk)
582               icoord(ierror,3) = jk
583            ENDIF
584         END DO
585         DO jl = 1, npcoa(2,jk)
586            IF(nicoa(jl,2,jk)-2 < 1-iind ) THEN
587               ierror = ierror + 1
588               icoord(ierror,1) = nicoa(jl,2,jk)
589               icoord(ierror,2) = njcoa(jl,2,jk)
590               icoord(ierror,3) = jk
591            ENDIF
592         END DO
593         DO jl = 1, npcoa(3,jk)
594            IF( njcoa(jl,3,jk)+3 > jpj+ijnd ) THEN
595               ierror = ierror + 1
596               icoord(ierror,1) = nicoa(jl,3,jk)
597               icoord(ierror,2) = njcoa(jl,3,jk)
598               icoord(ierror,3) = jk
599            ENDIF
600         END DO
601         DO jl = 1, npcoa(4,jk)
602            IF( njcoa(jl,4,jk)-2 < 1) THEN
603               ierror=ierror + 1
604               icoord(ierror,1) = nicoa(jl,4,jk)
605               icoord(ierror,2) = njcoa(jl,4,jk)
606               icoord(ierror,3) = jk
607            ENDIF
608         END DO
609      END DO
610     
611      IF( ierror > 0 ) THEN
612         IF(lwp) WRITE(numout,*)
613         IF(lwp) WRITE(numout,*) '              Problem on lateral conditions'
614         IF(lwp) WRITE(numout,*) '                 Bad marking off at points:'
615         DO jl = 1, ierror
616            IF(lwp) WRITE(numout,*) 'Level:',icoord(jl,3),   &
617               &                  '  Point(',icoord(jl,1),',',icoord(jl,2),')'
618         END DO
619         CALL ctl_stop( 'We stop...' )
620      ENDIF
621      !
622   END SUBROUTINE dom_msk_nsa
623
624#else
625   !!----------------------------------------------------------------------
626   !!   Default option :                                      Empty routine
627   !!----------------------------------------------------------------------
628   SUBROUTINE dom_msk_nsa       
629   END SUBROUTINE dom_msk_nsa
630#endif
631   
632   !!======================================================================
633END MODULE dommsk
Note: See TracBrowser for help on using the repository browser.