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

source: branches/UKMO/dev_5518_shlat2d/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90 @ 12787

Last change on this file since 12787 was 6568, checked in by mathiot, 8 years ago

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