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

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

dev_LOCEAN_CMCC_2011:Move one namelist parameter from dynvor to dommsk

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