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

source: branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90 @ 8664

Last change on this file since 8664 was 8664, checked in by jchanut, 6 years ago

Set fmask to 0 beyond AGRIF open boundaries

  • Property svn:keywords set to Id
File size: 31.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   !!            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  ::   isrow                    ! index for ORCA1 starting row
137      INTEGER , POINTER, DIMENSION(:,:) ::  imsk
138      REAL(wp), POINTER, DIMENSION(:,:) ::  zwf
139      !!
140      NAMELIST/namlbc/ rn_shlat, ln_vorlat
141      !!---------------------------------------------------------------------
142      !
143      IF( nn_timing == 1 )  CALL timing_start('dom_msk')
144      !
145      CALL wrk_alloc( jpi, jpj, imsk )
146      CALL wrk_alloc( jpi, jpj, zwf  )
147      !
148      REWIND( numnam_ref )              ! Namelist namlbc in reference namelist : Lateral momentum boundary condition
149      READ  ( numnam_ref, namlbc, IOSTAT = ios, ERR = 901 )
150901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namlbc in reference namelist', lwp )
151
152      REWIND( numnam_cfg )              ! Namelist namlbc in configuration namelist : Lateral momentum boundary condition
153      READ  ( numnam_cfg, namlbc, IOSTAT = ios, ERR = 902 )
154902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namlbc in configuration namelist', lwp )
155      IF(lwm) WRITE ( numond, namlbc )
156     
157      IF(lwp) THEN                  ! control print
158         WRITE(numout,*)
159         WRITE(numout,*) 'dommsk : ocean mask '
160         WRITE(numout,*) '~~~~~~'
161         WRITE(numout,*) '   Namelist namlbc'
162         WRITE(numout,*) '      lateral momentum boundary cond.    rn_shlat  = ',rn_shlat
163         WRITE(numout,*) '      consistency with analytical form   ln_vorlat = ',ln_vorlat 
164      ENDIF
165
166      IF     (      rn_shlat == 0.               ) THEN   ;   IF(lwp) WRITE(numout,*) '   ocean lateral  free-slip '
167      ELSEIF (      rn_shlat == 2.               ) THEN   ;   IF(lwp) WRITE(numout,*) '   ocean lateral  no-slip '
168      ELSEIF ( 0. < rn_shlat .AND. rn_shlat < 2. ) THEN   ;   IF(lwp) WRITE(numout,*) '   ocean lateral  partial-slip '
169      ELSEIF ( 2. < rn_shlat                     ) THEN   ;   IF(lwp) WRITE(numout,*) '   ocean lateral  strong-slip '
170      ELSE
171         WRITE(ctmp1,*) ' rn_shlat is negative = ', rn_shlat
172         CALL ctl_stop( ctmp1 )
173      ENDIF
174
175      ! 1. Ocean/land mask at t-point (computed from mbathy)
176      ! -----------------------------
177      ! N.B. tmask has already the right boundary conditions since mbathy is ok
178      !
179      tmask(:,:,:) = 0._wp
180      DO jk = 1, jpk
181         DO jj = 1, jpj
182            DO ji = 1, jpi
183               IF( REAL( mbathy(ji,jj) - jk, wp ) + 0.1_wp >= 0._wp )   tmask(ji,jj,jk) = 1._wp
184            END DO 
185         END DO 
186      END DO 
187     
188      ! (ISF) define barotropic mask and mask the ice shelf point
189      ssmask(:,:)=tmask(:,:,1) ! at this stage ice shelf is not masked
190     
191      DO jk = 1, jpk
192         DO jj = 1, jpj
193            DO ji = 1, jpi
194               IF( REAL( misfdep(ji,jj) - jk, wp ) - 0.1_wp >= 0._wp )   THEN
195                  tmask(ji,jj,jk) = 0._wp
196               END IF
197            END DO 
198         END DO 
199      END DO 
200
201!!gm  ????
202#if defined key_zdfkpp
203      IF( cp_cfg == 'orca' ) THEN
204         IF( jp_cfg == 2 )   THEN       ! land point on Bab el Mandeb zonal section
205            ij0 =  87   ;   ij1 =  88
206            ii0 = 160   ;   ii1 = 161
207            tmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 0._wp
208         ELSE
209            IF(lwp) WRITE(numout,*)
210            IF(lwp) WRITE(numout,cform_war)
211            IF(lwp) WRITE(numout,*)
212            IF(lwp) WRITE(numout,*)'          A mask must be applied on Bab el Mandeb strait'
213            IF(lwp) WRITE(numout,*)'          in case of ORCAs configurations'
214            IF(lwp) WRITE(numout,*)'          This is a problem which is not yet solved'
215            IF(lwp) WRITE(numout,*)
216         ENDIF
217      ENDIF
218#endif
219!!gm end
220
221      ! Interior domain mask (used for global sum)
222      ! --------------------
223      tmask_i(:,:) = ssmask(:,:)            ! (ISH) tmask_i = 1 even on the ice shelf
224      iif = jpreci                         ! ???
225      iil = nlci - jpreci + 1
226      ijf = jprecj                         ! ???
227      ijl = nlcj - jprecj + 1
228
229      tmask_i( 1 :iif,   :   ) = 0._wp      ! first columns
230      tmask_i(iil:jpi,   :   ) = 0._wp      ! last  columns (including mpp extra columns)
231      tmask_i(   :   , 1 :ijf) = 0._wp      ! first rows
232      tmask_i(   :   ,ijl:jpj) = 0._wp      ! last  rows (including mpp extra rows)
233
234      ! north fold mask
235      ! ---------------
236      tpol(1:jpiglo) = 1._wp 
237      fpol(1:jpiglo) = 1._wp
238      IF( jperio == 3 .OR. jperio == 4 ) THEN      ! T-point pivot
239         tpol(jpiglo/2+1:jpiglo) = 0._wp
240         fpol(     1    :jpiglo) = 0._wp
241         IF( mjg(nlej) == jpjglo ) THEN                  ! only half of the nlcj-1 row
242            DO ji = iif+1, iil-1
243               tmask_i(ji,nlej-1) = tmask_i(ji,nlej-1) * tpol(mig(ji))
244            END DO
245         ENDIF
246      ENDIF
247      IF( jperio == 5 .OR. jperio == 6 ) THEN      ! F-point pivot
248         tpol(     1    :jpiglo) = 0._wp
249         fpol(jpiglo/2+1:jpiglo) = 0._wp
250      ENDIF
251
252      ! 2. Ocean/land mask at u-,  v-, and z-points (computed from tmask)
253      ! -------------------------------------------
254      DO jk = 1, jpk
255         DO jj = 1, jpjm1
256            DO ji = 1, fs_jpim1   ! vector loop
257               umask(ji,jj,jk) = tmask(ji,jj  ,jk) * tmask(ji+1,jj  ,jk)
258               vmask(ji,jj,jk) = tmask(ji,jj  ,jk) * tmask(ji  ,jj+1,jk)
259            END DO
260            DO ji = 1, jpim1      ! NO vector opt.
261               fmask(ji,jj,jk) = tmask(ji,jj  ,jk) * tmask(ji+1,jj  ,jk)   &
262                  &            * tmask(ji,jj+1,jk) * tmask(ji+1,jj+1,jk)
263            END DO
264         END DO
265      END DO
266      ! (ISF) MIN(1,SUM(umask)) is here to check if you have effectively at least 1 wet u point
267      DO jj = 1, jpjm1
268         DO ji = 1, fs_jpim1   ! vector loop
269            umask_i(ji,jj)  = ssmask(ji,jj) * ssmask(ji+1,jj  )  * MIN(1._wp,SUM(umask(ji,jj,:)))
270            vmask_i(ji,jj)  = ssmask(ji,jj) * ssmask(ji  ,jj+1)  * MIN(1._wp,SUM(vmask(ji,jj,:)))
271         END DO
272         DO ji = 1, jpim1      ! NO vector opt.
273            fmask_i(ji,jj) =  ssmask(ji,jj  ) * ssmask(ji+1,jj  )   &
274               &            * ssmask(ji,jj+1) * ssmask(ji+1,jj+1) * MIN(1._wp,SUM(fmask(ji,jj,:)))
275         END DO
276      END DO
277      CALL lbc_lnk( umask, 'U', 1._wp )      ! Lateral boundary conditions
278      CALL lbc_lnk( vmask, 'V', 1._wp )
279      CALL lbc_lnk( fmask, 'F', 1._wp )
280      CALL lbc_lnk( umask_i, 'U', 1._wp )      ! Lateral boundary conditions
281      CALL lbc_lnk( vmask_i, 'V', 1._wp )
282      CALL lbc_lnk( fmask_i, 'F', 1._wp )
283
284      ! 3. Ocean/land mask at wu-, wv- and w points
285      !----------------------------------------------
286      wmask (:,:,1) = tmask(:,:,1) ! ????????
287      wumask(:,:,1) = umask(:,:,1) ! ????????
288      wvmask(:,:,1) = vmask(:,:,1) ! ????????
289      DO jk=2,jpk
290         wmask (:,:,jk)=tmask(:,:,jk) * tmask(:,:,jk-1)
291         wumask(:,:,jk)=umask(:,:,jk) * umask(:,:,jk-1)   
292         wvmask(:,:,jk)=vmask(:,:,jk) * vmask(:,:,jk-1)
293      END DO
294
295      ! 4. ocean/land mask for the elliptic equation
296      ! --------------------------------------------
297      bmask(:,:) = ssmask(:,:)       ! elliptic equation is written at t-point
298      !
299      !                               ! Boundary conditions
300      !                                    ! cyclic east-west : bmask must be set to 0. on rows 1 and jpi
301      IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN
302         bmask( 1 ,:) = 0._wp
303         bmask(jpi,:) = 0._wp
304      ENDIF
305      IF( nperio == 2 ) THEN               ! south symmetric :  bmask must be set to 0. on row 1
306         bmask(:, 1 ) = 0._wp
307      ENDIF
308      !                                    ! north fold :
309      IF( nperio == 3 .OR. nperio == 4 ) THEN   ! T-pt pivot : bmask set to 0. on row jpj and on half jpjglo-1 row
310         DO ji = 1, jpi                     
311            ii = ji + nimpp - 1
312            bmask(ji,jpj-1) = bmask(ji,jpj-1) * tpol(ii)
313            bmask(ji,jpj  ) = 0._wp
314         END DO
315      ENDIF
316      IF( nperio == 5 .OR. nperio == 6 ) THEN   ! F-pt pivot and T-pt elliptic eq. : bmask set to 0. on row jpj
317         bmask(:,jpj) = 0._wp
318      ENDIF
319      !
320      IF( lk_mpp ) THEN                    ! mpp specificities
321         !                                      ! bmask is set to zero on the overlap region
322         IF( nbondi /= -1 .AND. nbondi /= 2 )   bmask(  1 :jpreci,:) = 0._wp
323         IF( nbondi /=  1 .AND. nbondi /= 2 )   bmask(nlci:jpi   ,:) = 0._wp
324         IF( nbondj /= -1 .AND. nbondj /= 2 )   bmask(:,  1 :jprecj) = 0._wp
325         IF( nbondj /=  1 .AND. nbondj /= 2 )   bmask(:,nlcj:jpj   ) = 0._wp
326         !
327         IF( npolj == 3 .OR. npolj == 4 ) THEN  ! north fold : bmask must be set to 0. on rows jpj-1 and jpj
328            DO ji = 1, nlci
329               ii = ji + nimpp - 1
330               bmask(ji,nlcj-1) = bmask(ji,nlcj-1) * tpol(ii)
331               bmask(ji,nlcj  ) = 0._wp
332            END DO
333         ENDIF
334         IF( npolj == 5 .OR. npolj == 6 ) THEN  ! F-pt pivot and T-pt elliptic eq. : bmask set to 0. on row jpj
335            DO ji = 1, nlci
336               bmask(ji,nlcj  ) = 0._wp
337            END DO
338         ENDIF
339      ENDIF
340
341
342      ! mask for second order calculation of vorticity
343      ! ----------------------------------------------
344      CALL dom_msk_nsa
345
346     
347      ! Lateral boundary conditions on velocity (modify fmask)
348      ! ---------------------------------------     
349      DO jk = 1, jpk
350         zwf(:,:) = fmask(:,:,jk)         
351         DO jj = 2, jpjm1
352            DO ji = fs_2, fs_jpim1   ! vector opt.
353               IF( fmask(ji,jj,jk) == 0._wp ) THEN
354                  fmask(ji,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,jj), zwf(ji,jj+1),   &
355                     &                                           zwf(ji-1,jj), zwf(ji,jj-1)  )  )
356               ENDIF
357            END DO
358         END DO
359         DO jj = 2, jpjm1
360            IF( fmask(1,jj,jk) == 0._wp ) THEN
361               fmask(1  ,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1) ) )
362            ENDIF
363            IF( fmask(jpi,jj,jk) == 0._wp ) THEN
364               fmask(jpi,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) )
365            ENDIF
366         END DO         
367         DO ji = 2, jpim1
368            IF( fmask(ji,1,jk) == 0._wp ) THEN
369               fmask(ji, 1 ,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1) ) )
370            ENDIF
371            IF( fmask(ji,jpj,jk) == 0._wp ) THEN
372               fmask(ji,jpj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpjm1) ) )
373            ENDIF
374         END DO
375#if defined key_agrif
376         IF( .NOT. AGRIF_Root() ) THEN
377            IF ((nbondi ==  1).OR.(nbondi == 2)) fmask(nlci-1 , :     ,jk) = 0.e0      ! east
378            IF ((nbondi == -1).OR.(nbondi == 2)) fmask(1      , :     ,jk) = 0.e0      ! west
379            IF ((nbondj ==  1).OR.(nbondj == 2)) fmask(:      ,nlcj-1 ,jk) = 0.e0      ! north
380            IF ((nbondj == -1).OR.(nbondj == 2)) fmask(:      ,1      ,jk) = 0.e0      ! south
381         ENDIF
382#endif
383      END DO
384      !
385      IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN   ! ORCA_R2 configuration
386         !                                                 ! Increased lateral friction near of some straits
387         IF( nn_cla == 0 ) THEN
388            !                                ! Gibraltar strait  : partial slip (fmask=0.5)
389            ij0 = 101   ;   ij1 = 101
390            ii0 = 139   ;   ii1 = 140   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  0.5_wp
391            ij0 = 102   ;   ij1 = 102
392            ii0 = 139   ;   ii1 = 140   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  0.5_wp
393            !
394            !                                ! Bab el Mandeb : partial slip (fmask=1)
395            ij0 =  87   ;   ij1 =  88
396            ii0 = 160   ;   ii1 = 160   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  1._wp
397            ij0 =  88   ;   ij1 =  88
398            ii0 = 159   ;   ii1 = 159   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  1._wp
399            !
400         ENDIF
401         !                                ! Danish straits  : strong slip (fmask > 2)
402! We keep this as an example but it is instable in this case
403!         ij0 = 115   ;   ij1 = 115
404!         ii0 = 145   ;   ii1 = 146   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 4._wp
405!         ij0 = 116   ;   ij1 = 116
406!         ii0 = 145   ;   ii1 = 146   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 4._wp
407         !
408      ENDIF
409      !
410      IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN   ! ORCA R1 configuration
411         !                                                 ! Increased lateral friction near of some straits
412         ! This dirty section will be suppressed by simplification process:
413         ! all this will come back in input files
414         ! Currently these hard-wired indices relate to configuration with
415         ! extend grid (jpjglo=332)
416         !
417         isrow = 332 - jpjglo
418         !
419         IF(lwp) WRITE(numout,*)
420         IF(lwp) WRITE(numout,*) '   orca_r1: increase friction near the following straits : '
421         IF(lwp) WRITE(numout,*) '      Gibraltar '
422         ii0 = 282           ;   ii1 = 283        ! Gibraltar Strait
423         ij0 = 241 - isrow   ;   ij1 = 241 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp 
424
425         IF(lwp) WRITE(numout,*) '      Bhosporus '
426         ii0 = 314           ;   ii1 = 315        ! Bhosporus Strait
427         ij0 = 248 - isrow   ;   ij1 = 248 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp 
428
429         IF(lwp) WRITE(numout,*) '      Makassar (Top) '
430         ii0 =  48           ;   ii1 =  48        ! Makassar Strait (Top)
431         ij0 = 189 - isrow   ;   ij1 = 190 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp 
432
433         IF(lwp) WRITE(numout,*) '      Lombok '
434         ii0 =  44           ;   ii1 =  44        ! Lombok Strait
435         ij0 = 164 - isrow   ;   ij1 = 165 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp 
436
437         IF(lwp) WRITE(numout,*) '      Ombai '
438         ii0 =  53           ;   ii1 =  53        ! Ombai Strait
439         ij0 = 164 - isrow   ;   ij1 = 165 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp 
440
441         IF(lwp) WRITE(numout,*) '      Timor Passage '
442         ii0 =  56           ;   ii1 =  56        ! Timor Passage
443         ij0 = 164 - isrow   ;   ij1 = 165 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp 
444
445         IF(lwp) WRITE(numout,*) '      West Halmahera '
446         ii0 =  58           ;   ii1 =  58        ! West Halmahera Strait
447         ij0 = 181 - isrow   ;   ij1 = 182 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp 
448
449         IF(lwp) WRITE(numout,*) '      East Halmahera '
450         ii0 =  55           ;   ii1 =  55        ! East Halmahera Strait
451         ij0 = 181 - isrow   ;   ij1 = 182 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp 
452         !
453      ENDIF
454      !
455      CALL lbc_lnk( fmask, 'F', 1._wp )      ! Lateral boundary conditions on fmask
456
457      ! CAUTION : The fmask may be further modified in dyn_vor_init ( dynvor.F90 )
458           
459      IF( nprint == 1 .AND. lwp ) THEN      ! Control print
460         imsk(:,:) = INT( tmask_i(:,:) )
461         WRITE(numout,*) ' tmask_i : '
462         CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1,   &
463               &                           1, jpj, 1, 1, numout)
464         WRITE (numout,*)
465         WRITE (numout,*) ' dommsk: tmask for each level'
466         WRITE (numout,*) ' ----------------------------'
467         DO jk = 1, jpk
468            imsk(:,:) = INT( tmask(:,:,jk) )
469
470            WRITE(numout,*)
471            WRITE(numout,*) ' level = ',jk
472            CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1,   &
473               &                              1, jpj, 1, 1, numout)
474         END DO
475         WRITE(numout,*)
476         WRITE(numout,*) ' dom_msk: vmask for each level'
477         WRITE(numout,*) ' -----------------------------'
478         DO jk = 1, jpk
479            imsk(:,:) = INT( vmask(:,:,jk) )
480            WRITE(numout,*)
481            WRITE(numout,*) ' level = ',jk
482            CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1,   &
483               &                              1, jpj, 1, 1, numout)
484         END DO
485         WRITE(numout,*)
486         WRITE(numout,*) ' dom_msk: fmask for each level'
487         WRITE(numout,*) ' -----------------------------'
488         DO jk = 1, jpk
489            imsk(:,:) = INT( fmask(:,:,jk) )
490            WRITE(numout,*)
491            WRITE(numout,*) ' level = ',jk
492            CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1,   &
493               &                              1, jpj, 1, 1, numout )
494         END DO
495         WRITE(numout,*)
496         WRITE(numout,*) ' dom_msk: bmask '
497         WRITE(numout,*) ' ---------------'
498         WRITE(numout,*)
499         imsk(:,:) = INT( bmask(:,:) )
500         CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1,   &
501            &                              1, jpj, 1, 1, numout )
502      ENDIF
503      !
504      CALL wrk_dealloc( jpi, jpj, imsk )
505      CALL wrk_dealloc( jpi, jpj, zwf  )
506      !
507      IF( nn_timing == 1 )  CALL timing_stop('dom_msk')
508      !
509   END SUBROUTINE dom_msk
510
511#if defined key_noslip_accurate
512   !!----------------------------------------------------------------------
513   !!   'key_noslip_accurate' :         accurate no-slip boundary condition
514   !!----------------------------------------------------------------------
515   
516   SUBROUTINE dom_msk_nsa
517      !!---------------------------------------------------------------------
518      !!                 ***  ROUTINE dom_msk_nsa  ***
519      !!
520      !! ** Purpose :
521      !!
522      !! ** Method  :
523      !!
524      !! ** Action :
525      !!----------------------------------------------------------------------
526      INTEGER  ::   ji, jj, jk, jl      ! dummy loop indices
527      INTEGER  ::   ine, inw, ins, inn, itest, ierror, iind, ijnd
528      REAL(wp) ::   zaa
529      !!---------------------------------------------------------------------
530      !
531      IF( nn_timing == 1 )  CALL timing_start('dom_msk_nsa')
532      !
533      IF(lwp) WRITE(numout,*)
534      IF(lwp) WRITE(numout,*) 'dom_msk_nsa : noslip accurate boundary condition'
535      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   using Schchepetkin and O Brian scheme'
536      IF( lk_mpp )   CALL ctl_stop( ' mpp version is not yet implemented' )
537
538      ! mask for second order calculation of vorticity
539      ! ----------------------------------------------
540      ! noslip boundary condition: fmask=1  at convex corner, store
541      ! index of straight coast meshes ( 'west', refering to a coast,
542      ! means west of the ocean, aso)
543     
544      DO jk = 1, jpk
545         DO jl = 1, 4
546            npcoa(jl,jk) = 0
547            DO ji = 1, 2*(jpi+jpj)
548               nicoa(ji,jl,jk) = 0
549               njcoa(ji,jl,jk) = 0
550            END DO
551         END DO
552      END DO
553     
554      IF( jperio == 2 ) THEN
555         WRITE(numout,*) ' '
556         WRITE(numout,*) ' symetric boundary conditions need special'
557         WRITE(numout,*) ' treatment not implemented. we stop.'
558         STOP
559      ENDIF
560     
561      ! convex corners
562     
563      DO jk = 1, jpkm1
564         DO jj = 1, jpjm1
565            DO ji = 1, jpim1
566               zaa = tmask(ji  ,jj,jk) + tmask(ji  ,jj+1,jk)   &
567                  &+ tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk)
568               IF( ABS(zaa-3._wp) <= 0.1_wp )   fmask(ji,jj,jk) = 1._wp
569            END DO
570         END DO
571      END DO
572
573      ! north-south straight coast
574
575      DO jk = 1, jpkm1
576         inw = 0
577         ine = 0
578         DO jj = 2, jpjm1
579            DO ji = 2, jpim1
580               zaa = tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk)
581               IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN
582                  inw = inw + 1
583                  nicoa(inw,1,jk) = ji
584                  njcoa(inw,1,jk) = jj
585                  IF( nprint == 1 ) WRITE(numout,*) ' west  : ', jk, inw, ji, jj
586               ENDIF
587               zaa = tmask(ji,jj,jk) + tmask(ji,jj+1,jk)
588               IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN
589                  ine = ine + 1
590                  nicoa(ine,2,jk) = ji
591                  njcoa(ine,2,jk) = jj
592                  IF( nprint == 1 ) WRITE(numout,*) ' east  : ', jk, ine, ji, jj
593               ENDIF
594            END DO
595         END DO
596         npcoa(1,jk) = inw
597         npcoa(2,jk) = ine
598      END DO
599
600      ! west-east straight coast
601
602      DO jk = 1, jpkm1
603         ins = 0
604         inn = 0
605         DO jj = 2, jpjm1
606            DO ji =2, jpim1
607               zaa = tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk)
608               IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN
609                  ins = ins + 1
610                  nicoa(ins,3,jk) = ji
611                  njcoa(ins,3,jk) = jj
612                  IF( nprint == 1 ) WRITE(numout,*) ' south : ', jk, ins, ji, jj
613               ENDIF
614               zaa = tmask(ji+1,jj,jk) + tmask(ji,jj,jk)
615               IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN
616                  inn = inn + 1
617                  nicoa(inn,4,jk) = ji
618                  njcoa(inn,4,jk) = jj
619                  IF( nprint == 1 ) WRITE(numout,*) ' north : ', jk, inn, ji, jj
620               ENDIF
621            END DO
622         END DO
623         npcoa(3,jk) = ins
624         npcoa(4,jk) = inn
625      END DO
626
627      itest = 2 * ( jpi + jpj )
628      DO jk = 1, jpk
629         IF( npcoa(1,jk) > itest .OR. npcoa(2,jk) > itest .OR.   &
630             npcoa(3,jk) > itest .OR. npcoa(4,jk) > itest ) THEN
631           
632            WRITE(ctmp1,*) ' level jk = ',jk
633            WRITE(ctmp2,*) ' straight coast index arraies are too small.:'
634            WRITE(ctmp3,*) ' npe, npw, nps, npn = ', npcoa(1,jk), npcoa(2,jk),   &
635                &                                     npcoa(3,jk), npcoa(4,jk)
636            WRITE(ctmp4,*) ' 2*(jpi+jpj) = ',itest,'. we stop.'
637            CALL ctl_stop( ctmp1, ctmp2, ctmp3, ctmp4 )
638        ENDIF
639      END DO
640
641      ierror = 0
642      iind = 0
643      ijnd = 0
644      IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 )   iind = 2
645      IF( nperio == 3 .OR. nperio == 4 .OR. nperio == 5 .OR. nperio == 6 )   ijnd = 2
646      DO jk = 1, jpk
647         DO jl = 1, npcoa(1,jk)
648            IF( nicoa(jl,1,jk)+3 > jpi+iind ) THEN
649               ierror = ierror+1
650               icoord(ierror,1) = nicoa(jl,1,jk)
651               icoord(ierror,2) = njcoa(jl,1,jk)
652               icoord(ierror,3) = jk
653            ENDIF
654         END DO
655         DO jl = 1, npcoa(2,jk)
656            IF(nicoa(jl,2,jk)-2 < 1-iind ) THEN
657               ierror = ierror + 1
658               icoord(ierror,1) = nicoa(jl,2,jk)
659               icoord(ierror,2) = njcoa(jl,2,jk)
660               icoord(ierror,3) = jk
661            ENDIF
662         END DO
663         DO jl = 1, npcoa(3,jk)
664            IF( njcoa(jl,3,jk)+3 > jpj+ijnd ) THEN
665               ierror = ierror + 1
666               icoord(ierror,1) = nicoa(jl,3,jk)
667               icoord(ierror,2) = njcoa(jl,3,jk)
668               icoord(ierror,3) = jk
669            ENDIF
670         END DO
671         DO jl = 1, npcoa(4,jk)
672            IF( njcoa(jl,4,jk)-2 < 1) THEN
673               ierror=ierror + 1
674               icoord(ierror,1) = nicoa(jl,4,jk)
675               icoord(ierror,2) = njcoa(jl,4,jk)
676               icoord(ierror,3) = jk
677            ENDIF
678         END DO
679      END DO
680     
681      IF( ierror > 0 ) THEN
682         IF(lwp) WRITE(numout,*)
683         IF(lwp) WRITE(numout,*) '              Problem on lateral conditions'
684         IF(lwp) WRITE(numout,*) '                 Bad marking off at points:'
685         DO jl = 1, ierror
686            IF(lwp) WRITE(numout,*) 'Level:',icoord(jl,3),   &
687               &                  '  Point(',icoord(jl,1),',',icoord(jl,2),')'
688         END DO
689         CALL ctl_stop( 'We stop...' )
690      ENDIF
691      !
692      IF( nn_timing == 1 )  CALL timing_stop('dom_msk_nsa')
693      !
694   END SUBROUTINE dom_msk_nsa
695
696#else
697   !!----------------------------------------------------------------------
698   !!   Default option :                                      Empty routine
699   !!----------------------------------------------------------------------
700   SUBROUTINE dom_msk_nsa       
701   END SUBROUTINE dom_msk_nsa
702#endif
703   
704   !!======================================================================
705END MODULE dommsk
Note: See TracBrowser for help on using the repository browser.