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

source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90 @ 2380

Last change on this file since 2380 was 2380, checked in by acc, 13 years ago

nemo_v3_3beta. ORCA_R1 settings (step 2, see ticket #758). Introduces key_orca_r1 (46 level default, 75 level if key_orca_r1=75)
  • Property svn:keywords set to Id
File size: 27.2 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   !!             -   ! 1993-03  (M. Guyon)  symetrical conditions (M. Guyon)
8   !!             -   ! 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   !!             -   ! 1997-02  (G. Madec)  mesh information put in domhgr.F
12   !!             -   ! 1997-07  (G. Madec)  modification of mbathy and fmask
13   !!             -   ! 1998-05  (G. Roullet)  free surface
14   !!             -   ! 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
38   !                            !!* Namelist namlbc : lateral boundary condition *
39   REAL(wp) ::   rn_shlat = 2.   ! type of lateral boundary condition on velocity
40   
41   !! * Substitutions
42#  include "vectopt_loop_substitute.h90"
43   !!----------------------------------------------------------------------
44   !! NEMO/OPA 3.2 , LODYC-IPSL  (2009)
45   !! $Id$
46   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
47   !!----------------------------------------------------------------------
48
49CONTAINS
50   
51   SUBROUTINE dom_msk
52      !!---------------------------------------------------------------------
53      !!                 ***  ROUTINE dom_msk  ***
54      !!
55      !! ** Purpose :   Compute land/ocean mask arrays at tracer points, hori-
56      !!      zontal velocity points (u & v), vorticity points (f) and baro-
57      !!      tropic stream function  points (b).
58      !!        Set mbathy to the number of non-zero w-levels of a water column
59      !!
60      !! ** Method  :   The ocean/land mask is computed from the basin bathy-
61      !!      metry in level (mbathy) which is defined or read in dommba.
62      !!      mbathy equals 0 over continental T-point
63      !!      and the number of ocean level over the ocean.
64      !!
65      !!      At a given position (ji,jj,jk) the ocean/land mask is given by:
66      !!      t-point : 0. IF mbathy( ji ,jj) =< 0
67      !!                1. IF mbathy( ji ,jj) >= jk
68      !!      u-point : 0. IF mbathy( ji ,jj)  or mbathy(ji+1, jj ) =< 0
69      !!                1. IF mbathy( ji ,jj) and mbathy(ji+1, jj ) >= jk.
70      !!      v-point : 0. IF mbathy( ji ,jj)  or mbathy( ji ,jj+1) =< 0
71      !!                1. IF mbathy( ji ,jj) and mbathy( ji ,jj+1) >= jk.
72      !!      f-point : 0. IF mbathy( ji ,jj)  or mbathy( ji ,jj+1)
73      !!                   or mbathy(ji+1,jj)  or mbathy(ji+1,jj+1) =< 0
74      !!                1. IF mbathy( ji ,jj) and mbathy( ji ,jj+1)
75      !!                and mbathy(ji+1,jj) and mbathy(ji+1,jj+1) >= jk.
76      !!      b-point : the same definition as for f-point of the first ocean
77      !!                level (surface level) but with 0 along coastlines.
78      !!
79      !!        The lateral friction is set through the value of fmask along
80      !!      the coast and topography. This value is defined by rn_shlat, a
81      !!      namelist parameter:
82      !!         rn_shlat = 0, free slip  (no shear along the coast)
83      !!         rn_shlat = 2, no slip  (specified zero velocity at the coast)
84      !!         0 < rn_shlat < 2, partial slip   | non-linear velocity profile
85      !!         2 < rn_shlat, strong slip        | in the lateral boundary layer
86      !!
87      !!      N.B. If nperio not equal to 0, the land/ocean mask arrays
88      !!      are defined with the proper value at lateral domain boundaries,
89      !!      but bmask. indeed, bmask defined the domain over which the
90      !!      barotropic stream function is computed. this domain cannot
91      !!      contain identical columns because the matrix associated with
92      !!      the barotropic stream function equation is then no more inverti-
93      !!      ble. therefore bmask is set to 0 along lateral domain boundaries
94      !!      even IF nperio is not zero.
95      !!
96      !!      In case of open boundaries (lk_obc=T):
97      !!        - tmask is set to 1 on the points to be computed bay the open
98      !!          boundaries routines.
99      !!        - bmask is  set to 0 on the open boundaries.
100      !!
101      !!      Set mbathy to the number of non-zero w-levels of a water column
102      !!                  mbathy = min( mbathy, 1 ) + 1
103      !!      (note that the minimum value of mbathy is 2).
104      !!
105      !! ** Action :   tmask    : land/ocean mask at t-point (=0. or 1.)
106      !!               umask    : land/ocean mask at u-point (=0. or 1.)
107      !!               vmask    : land/ocean mask at v-point (=0. or 1.)
108      !!               fmask    : land/ocean mask at f-point (=0. or 1.)
109      !!                          =rn_shlat along lateral boundaries
110      !!               bmask    : land/ocean mask at barotropic stream
111      !!                          function point (=0. or 1.) and set to 0 along lateral boundaries
112      !!               mbathy   : number of non-zero w-levels
113      !!----------------------------------------------------------------------
114      INTEGER  ::   ji, jj, jk      ! dummy loop indices
115      INTEGER  ::   iif, iil, ii0, ii1, ii
116      INTEGER  ::   ijf, ijl, ij0, ij1
117      INTEGER , DIMENSION(jpi,jpj) ::  imsk
118      REAL(wp), DIMENSION(jpi,jpj) ::   zwf
119      !!
120      NAMELIST/namlbc/ rn_shlat
121      !!---------------------------------------------------------------------
122     
123      REWIND( numnam )              ! Namelist namlbc : lateral momentum boundary condition
124      READ  ( numnam, namlbc )
125     
126      IF(lwp) THEN                  ! control print
127         WRITE(numout,*)
128         WRITE(numout,*) 'dommsk : ocean mask '
129         WRITE(numout,*) '~~~~~~'
130         WRITE(numout,*) '   Namelist namlbc'
131         WRITE(numout,*) '      lateral momentum boundary cond.    rn_shlat = ',rn_shlat
132      ENDIF
133
134      IF     (      rn_shlat == 0.               ) THEN   ;   IF(lwp) WRITE(numout,*) '   ocean lateral free-slip '
135      ELSEIF (      rn_shlat == 2.               ) THEN   ;   IF(lwp) WRITE(numout,*) '   ocean lateral  no-slip '
136      ELSEIF ( 0. < rn_shlat .AND. rn_shlat < 2. ) THEN   ;   IF(lwp) WRITE(numout,*) '   ocean lateral  partial-slip '
137      ELSEIF ( 2. < rn_shlat                     ) THEN   ;   IF(lwp) WRITE(numout,*) '   ocean lateral  strong-slip '
138      ELSE
139         WRITE(ctmp1,*) ' rn_shlat is negative = ', rn_shlat
140         CALL ctl_stop( ctmp1 )
141      ENDIF
142
143      ! 1. Ocean/land mask at t-point (computed from mbathy)
144      ! -----------------------------
145      ! N.B. tmask has already the right boundary conditions since mbathy is ok
146      !
147      tmask(:,:,:) = 0.e0
148      DO jk = 1, jpk
149         DO jj = 1, jpj
150            DO ji = 1, jpi
151               IF( REAL( mbathy(ji,jj) - jk ) +.1 >= 0.e0 )   tmask(ji,jj,jk) = 1.e0
152            END DO 
153         END DO 
154      END DO 
155
156!!gm  ????
157#if defined key_zdfkpp
158      IF( cp_cfg == 'orca' ) THEN
159         IF( jp_cfg == 2 )   THEN       ! land point on Bab el Mandeb zonal section
160            ij0 =  87   ;   ij1 =  88
161            ii0 = 160   ;   ii1 = 161
162            tmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 0.e0
163         ELSE
164            IF(lwp) WRITE(numout,*)
165            IF(lwp) WRITE(numout,cform_war)
166            IF(lwp) WRITE(numout,*)
167            IF(lwp) WRITE(numout,*)'          A mask must be applied on Bab el Mandeb strait'
168            IF(lwp) WRITE(numout,*)'          in case of ORCAs configurations'
169            IF(lwp) WRITE(numout,*)'          This is a problem which is not yet solved'
170            IF(lwp) WRITE(numout,*)
171         ENDIF
172      ENDIF
173#endif
174!!gm end
175
176      ! Interior domain mask (used for global sum)
177      ! --------------------
178      tmask_i(:,:) = tmask(:,:,1)
179      iif = jpreci                         ! ???
180      iil = nlci - jpreci + 1
181      ijf = jprecj                         ! ???
182      ijl = nlcj - jprecj + 1
183
184      tmask_i( 1 :iif,   :   ) = 0.e0      ! first columns
185      tmask_i(iil:jpi,   :   ) = 0.e0      ! last  columns (including mpp extra columns)
186      tmask_i(   :   , 1 :ijf) = 0.e0      ! first rows
187      tmask_i(   :   ,ijl:jpj) = 0.e0      ! last  rows (including mpp extra rows)
188
189      ! north fold mask
190      ! ---------------
191      tpol(1:jpiglo) = 1.e0 
192      fpol(1:jpiglo) = 1.e0
193      IF( jperio == 3 .OR. jperio == 4 ) THEN      ! T-point pivot
194         tpol(jpiglo/2+1:jpiglo) = 0.e0
195         fpol(     1    :jpiglo) = 0.e0
196         IF( mjg(nlej) == jpjglo ) THEN                  ! only half of the nlcj-1 row
197            DO ji = iif+1, iil-1
198               tmask_i(ji,nlej-1) = tmask_i(ji,nlej-1) * tpol(mig(ji))
199            END DO
200         ENDIF
201      ENDIF
202      IF( jperio == 5 .OR. jperio == 6 ) THEN      ! F-point pivot
203         tpol(     1    :jpiglo) = 0.e0
204         fpol(jpiglo/2+1:jpiglo) = 0.e0
205      ENDIF
206
207      ! 2. Ocean/land mask at u-,  v-, and z-points (computed from tmask)
208      ! -------------------------------------------
209      DO jk = 1, jpk
210         DO jj = 1, jpjm1
211            DO ji = 1, fs_jpim1   ! vector loop
212               umask(ji,jj,jk) = tmask(ji,jj  ,jk) * tmask(ji+1,jj  ,jk)
213               vmask(ji,jj,jk) = tmask(ji,jj  ,jk) * tmask(ji  ,jj+1,jk)
214            END DO
215            DO ji = 1, jpim1      ! NO vector opt.
216               fmask(ji,jj,jk) = tmask(ji,jj  ,jk) * tmask(ji+1,jj  ,jk)   &
217                  &            * tmask(ji,jj+1,jk) * tmask(ji+1,jj+1,jk)
218            END DO
219         END DO
220      END DO
221      CALL lbc_lnk( umask, 'U', 1. )      ! Lateral boundary conditions
222      CALL lbc_lnk( vmask, 'V', 1. )
223      CALL lbc_lnk( fmask, 'F', 1. )
224
225
226      ! 4. ocean/land mask for the elliptic equation
227      ! --------------------------------------------
228      bmask(:,:) = tmask(:,:,1)       ! elliptic equation is written at t-point
229      !
230      !                               ! Boundary conditions
231      !                                    ! cyclic east-west : bmask must be set to 0. on rows 1 and jpi
232      IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN
233         bmask( 1 ,:) = 0.e0
234         bmask(jpi,:) = 0.e0
235      ENDIF
236      IF( nperio == 2 ) THEN               ! south symmetric :  bmask must be set to 0. on row 1
237         bmask(:, 1 ) = 0.e0
238      ENDIF
239      !                                    ! north fold :
240      IF( nperio == 3 .OR. nperio == 4 ) THEN   ! T-pt pivot : bmask set to 0. on row jpj and on half jpjglo-1 row
241         DO ji = 1, jpi                     
242            ii = ji + nimpp - 1
243            bmask(ji,jpj-1) = bmask(ji,jpj-1) * tpol(ii)
244            bmask(ji,jpj  ) = 0.e0
245         END DO
246      ENDIF
247      IF( nperio == 5 .OR. nperio == 6 ) THEN   ! F-pt pivot and T-pt elliptic eq. : bmask set to 0. on row jpj
248         bmask(:,jpj) = 0.e0
249      ENDIF
250      !
251      IF( lk_mpp ) THEN                    ! mpp specificities
252         !                                      ! bmask is set to zero on the overlap region
253         IF( nbondi /= -1 .AND. nbondi /= 2 )   bmask(  1 :jpreci,:) = 0.e0
254         IF( nbondi /=  1 .AND. nbondi /= 2 )   bmask(nlci:jpi   ,:) = 0.e0
255         IF( nbondj /= -1 .AND. nbondj /= 2 )   bmask(:,  1 :jprecj) = 0.e0
256         IF( nbondj /=  1 .AND. nbondj /= 2 )   bmask(:,nlcj:jpj   ) = 0.e0
257         !
258         IF( npolj == 3 .OR. npolj == 4 ) THEN  ! north fold : bmask must be set to 0. on rows jpj-1 and jpj
259            DO ji = 1, nlci
260               ii = ji + nimpp - 1
261               bmask(ji,nlcj-1) = bmask(ji,nlcj-1) * tpol(ii)
262               bmask(ji,nlcj  ) = 0.e0
263            END DO
264         ENDIF
265         IF( npolj == 5 .OR. npolj == 6 ) THEN  ! F-pt pivot and T-pt elliptic eq. : bmask set to 0. on row jpj
266            DO ji = 1, nlci
267               bmask(ji,nlcj  ) = 0.e0
268            END DO
269         ENDIF
270      ENDIF
271
272
273      ! mask for second order calculation of vorticity
274      ! ----------------------------------------------
275      CALL dom_msk_nsa
276
277     
278      ! Lateral boundary conditions on velocity (modify fmask)
279      ! ---------------------------------------     
280      DO jk = 1, jpk
281         zwf(:,:) = fmask(:,:,jk)         
282         DO jj = 2, jpjm1
283            DO ji = fs_2, fs_jpim1   ! vector opt.
284               IF( fmask(ji,jj,jk) == 0. ) THEN
285                  fmask(ji,jj,jk) = rn_shlat * MIN( 1., MAX( zwf(ji+1,jj), zwf(ji,jj+1),   &
286                     &                                       zwf(ji-1,jj), zwf(ji,jj-1)  )  )
287               ENDIF
288            END DO
289         END DO
290         DO jj = 2, jpjm1
291            IF( fmask(1,jj,jk) == 0. ) THEN
292               fmask(1  ,jj,jk) = rn_shlat * MIN( 1., MAX( zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1) ) )
293            ENDIF
294            IF( fmask(jpi,jj,jk) == 0. ) THEN
295               fmask(jpi,jj,jk) = rn_shlat * MIN( 1., MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) )
296            ENDIF
297         END DO         
298         DO ji = 2, jpim1
299            IF( fmask(ji,1,jk) == 0. ) THEN
300               fmask(ji, 1 ,jk) = rn_shlat * MIN( 1., MAX( zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1) ) )
301            ENDIF
302            IF( fmask(ji,jpj,jk) == 0. ) THEN
303               fmask(ji,jpj,jk) = rn_shlat * MIN( 1., MAX( zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpjm1) ) )
304            ENDIF
305         END DO
306      END DO
307      !
308      IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN   ! ORCA_R2 configuration
309         !                                                 ! Increased lateral friction near of some straits
310         IF( n_cla == 0 ) THEN
311            !                                ! Gibraltar strait  : partial slip (fmask=0.5)
312            ij0 = 101   ;   ij1 = 101
313            ii0 = 139   ;   ii1 = 140   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  0.5e0
314            ij0 = 102   ;   ij1 = 102
315            ii0 = 139   ;   ii1 = 140   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  0.5e0
316            !
317            !                                ! Bab el Mandeb : partial slip (fmask=1)
318            ij0 =  87   ;   ij1 =  88
319            ii0 = 160   ;   ii1 = 160   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  1.e0
320            ij0 =  88   ;   ij1 =  88
321            ii0 = 159   ;   ii1 = 159   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  1.e0
322            !
323         ENDIF
324
325         !                                ! Danish straits  : strong slip (fmask > 2)
326! We keep this as an example but it is instable in this case
327!         ij0 = 115   ;   ij1 = 115
328!         ii0 = 145   ;   ii1 = 146   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 4.0e0
329!         ij0 = 116   ;   ij1 = 116
330!         ii0 = 145   ;   ii1 = 146   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 4.0e0
331         !
332      ENDIF
333            !                                               ! =====================
334         IF( cp_cfg == "orca" .AND. jp_cfg .eq. 1 ) THEN    ! ORCA R1 configuration
335            !                                               ! =====================
336
337            ii0 = 283   ;   ii1 = 284        ! Gibraltar Strait
338            ij0 = 200   ;   ij1 = 200   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) =  2.0 
339            IF(lwp) WRITE(numout,*)
340            IF(lwp) WRITE(numout,*) '             orca_r1: Gibraltar : '
341
342            ii0 = 314   ;   ii1 = 315        ! Bhosporus Strait
343            ij0 = 208   ;   ij1 = 208   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) =  2.0 
344            IF(lwp) WRITE(numout,*)
345            IF(lwp) WRITE(numout,*) '             orca_r1: Bhosporus : '
346
347            ii0 =  48   ;   ii1 =  48        ! Makassar Strait (Top)
348            ij0 = 149   ;   ij1 = 150   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) =  3.0 
349            IF(lwp) WRITE(numout,*)
350            IF(lwp) WRITE(numout,*) '             orca_r1: Makassar (Top) : '
351
352            ii0 =  44   ;   ii1 =  44        ! Lombok Strait
353            ij0 = 124   ;   ij1 = 125   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) =  2.0 
354            IF(lwp) WRITE(numout,*)
355            IF(lwp) WRITE(numout,*) '             orca_r1: Lombok : '
356
357            ii0 =  53   ;   ii1 =  53        ! Ombai Strait
358            ij0 = 124   ;   ij1 = 125   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) = 2.0 
359            IF(lwp) WRITE(numout,*)
360            IF(lwp) WRITE(numout,*) '             orca_r1: Ombai : '
361
362            ii0 =  56   ;   ii1 =  56        ! Timor Passage
363            ij0 = 124   ;   ij1 = 125   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) = 2.0 
364            IF(lwp) WRITE(numout,*)
365            IF(lwp) WRITE(numout,*) '             orca_r1: '
366
367            ii0 =  58   ;   ii1 =  58        ! West Halmahera Strait
368            ij0 = 141   ;   ij1 = 142   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) = 3.0 
369            IF(lwp) WRITE(numout,*)
370            IF(lwp) WRITE(numout,*) '             orca_r1: West Halmahera : '
371
372            ii0 =  55   ;   ii1 =  55        ! East Halmahera Strait
373            ij0 = 141   ;   ij1 = 142   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) = 3.0 
374            IF(lwp) WRITE(numout,*)
375            IF(lwp) WRITE(numout,*) '             orca_r1: East Halmahera : '
376
377            !
378            !
379         ENDIF
380      !
381      CALL lbc_lnk( fmask, 'F', 1. )      ! Lateral boundary conditions on fmask
382
383     
384      ! Mbathy set to the number of w-level (minimum value 2)
385      ! -----------------------------------
386      DO jj = 1, jpj
387         DO ji = 1, jpi
388            mbathy(ji,jj) = MAX( 1, mbathy(ji,jj) ) + 1
389         END DO
390      END DO
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   END SUBROUTINE dom_msk
438
439#if defined key_noslip_accurate
440   !!----------------------------------------------------------------------
441   !!   'key_noslip_accurate' :         accurate no-slip boundary condition
442   !!----------------------------------------------------------------------
443   
444   SUBROUTINE dom_msk_nsa
445      !!---------------------------------------------------------------------
446      !!                 ***  ROUTINE dom_msk_nsa  ***
447      !!
448      !! ** Purpose :
449      !!
450      !! ** Method  :
451      !!
452      !! ** Action :
453      !!
454      !!----------------------------------------------------------------------
455      INTEGER  :: ji, jj, jk, jl      ! dummy loop indices
456      INTEGER  ::   ine, inw, ins, inn, itest, ierror, iind, ijnd
457      REAL(wp) ::   zaa
458      INTEGER, DIMENSION(jpi*jpj*jpk,3) ::  icoord
459      !!---------------------------------------------------------------------
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.) <= 0.1 )   fmask(ji,jj,jk) = 1.
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.) <= 0.1 .AND. fmask(ji,jj,jk) == 0 ) 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.) <= 0.1 .AND. fmask(ji,jj,jk) == 0 ) 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.) <= 0.1 .AND. fmask(ji,jj,jk) == 0 ) 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.) <= 0.1 .AND. fmask(ji,jj,jk) == 0 ) 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.