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

source: branches/2015/dev_r5056_CMCC4_simplification/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90 @ 5282

Last change on this file since 5282 was 5282, checked in by diovino, 9 years ago

Dev. branch CMCC4_simplification ticket #1456

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