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

source: branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90 @ 2797

Last change on this file since 2797 was 2797, checked in by davestorkey, 13 years ago

Delete BDY module and first implementation of new OBC module.

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