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

source: branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90 @ 5770

Last change on this file since 5770 was 5770, checked in by gm, 9 years ago

#1593: LDF-ADV, step II.2: phasing the improvements/simplifications of advective tracer trend (see wiki page)

  • Property svn:keywords set to Id
File size: 30.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   !!            3.6  ! 2015-05  (P. Mathiot) ISF: add wmask,wumask and wvmask
20   !!----------------------------------------------------------------------
21
22   !!----------------------------------------------------------------------
23   !!   dom_msk        : compute land/ocean mask
24   !!   dom_msk_nsa    : update land/ocean mask when no-slip accurate option is used.
25   !!----------------------------------------------------------------------
26   USE oce             ! ocean dynamics and tracers
27   USE dom_oce         ! ocean space and time domain
28   USE in_out_manager  ! I/O manager
29   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
30   USE lib_mpp
31   USE dynspg_oce      ! choice/control of key cpp for surface pressure gradient
32   USE wrk_nemo        ! Memory allocation
33   USE timing          ! Timing
34
35   IMPLICIT NONE
36   PRIVATE
37
38   PUBLIC   dom_msk         ! routine called by inidom.F90
39   PUBLIC   dom_msk_alloc   ! routine called by nemogcm.F90
40
41   !                            !!* Namelist namlbc : lateral boundary condition *
42   REAL(wp)        :: rn_shlat   ! type of lateral boundary condition on velocity
43   LOGICAL, PUBLIC :: ln_vorlat  !  consistency of vorticity boundary condition
44   !                                            with analytical eqs.
45
46
47   INTEGER, ALLOCATABLE, SAVE, DIMENSION(:,:) ::  icoord ! Workspace for dom_msk_nsa()
48
49   !! * Substitutions
50#  include "vectopt_loop_substitute.h90"
51   !!----------------------------------------------------------------------
52   !! NEMO/OPA 3.2 , LODYC-IPSL  (2009)
53   !! $Id$
54   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
55   !!----------------------------------------------------------------------
56CONTAINS
57   
58   INTEGER FUNCTION dom_msk_alloc()
59      !!---------------------------------------------------------------------
60      !!                 ***  FUNCTION dom_msk_alloc  ***
61      !!---------------------------------------------------------------------
62      dom_msk_alloc = 0
63#if defined key_noslip_accurate
64      ALLOCATE(icoord(jpi*jpj*jpk,3), STAT=dom_msk_alloc)
65#endif
66      IF( dom_msk_alloc /= 0 )   CALL ctl_warn('dom_msk_alloc: failed to allocate icoord array')
67      !
68   END FUNCTION dom_msk_alloc
69
70
71   SUBROUTINE dom_msk
72      !!---------------------------------------------------------------------
73      !!                 ***  ROUTINE dom_msk  ***
74      !!
75      !! ** Purpose :   Compute land/ocean mask arrays at tracer points, hori-
76      !!      zontal velocity points (u & v), vorticity points (f) and baro-
77      !!      tropic stream function  points (b).
78      !!
79      !! ** Method  :   The ocean/land mask is computed from the basin bathy-
80      !!      metry in level (mbathy) which is defined or read in dommba.
81      !!      mbathy equals 0 over continental T-point
82      !!      and the number of ocean level over the ocean.
83      !!
84      !!      At a given position (ji,jj,jk) the ocean/land mask is given by:
85      !!      t-point : 0. IF mbathy( ji ,jj) =< 0
86      !!                1. IF mbathy( ji ,jj) >= jk
87      !!      u-point : 0. IF mbathy( ji ,jj)  or mbathy(ji+1, jj ) =< 0
88      !!                1. IF mbathy( ji ,jj) and mbathy(ji+1, jj ) >= jk.
89      !!      v-point : 0. IF mbathy( ji ,jj)  or mbathy( ji ,jj+1) =< 0
90      !!                1. IF mbathy( ji ,jj) and mbathy( ji ,jj+1) >= jk.
91      !!      f-point : 0. IF mbathy( ji ,jj)  or mbathy( ji ,jj+1)
92      !!                   or mbathy(ji+1,jj)  or mbathy(ji+1,jj+1) =< 0
93      !!                1. IF mbathy( ji ,jj) and mbathy( ji ,jj+1)
94      !!                  and mbathy(ji+1,jj) and mbathy(ji+1,jj+1) >= jk.
95      !!      b-point : the same definition as for f-point of the first ocean
96      !!                level (surface level) but with 0 along coastlines.
97      !!      tmask_i : interior ocean mask at t-point, i.e. excluding duplicated
98      !!                rows/lines due to cyclic or North Fold boundaries as well
99      !!                as MPP halos.
100      !!
101      !!        The lateral friction is set through the value of fmask along
102      !!      the coast and topography. This value is defined by rn_shlat, a
103      !!      namelist parameter:
104      !!         rn_shlat = 0, free slip  (no shear along the coast)
105      !!         rn_shlat = 2, no slip  (specified zero velocity at the coast)
106      !!         0 < rn_shlat < 2, partial slip   | non-linear velocity profile
107      !!         2 < rn_shlat, strong slip        | in the lateral boundary layer
108      !!
109      !!      N.B. If nperio not equal to 0, the land/ocean mask arrays
110      !!      are defined with the proper value at lateral domain boundaries,
111      !!      but bmask. indeed, bmask defined the domain over which the
112      !!      barotropic stream function is computed. this domain cannot
113      !!      contain identical columns because the matrix associated with
114      !!      the barotropic stream function equation is then no more inverti-
115      !!      ble. therefore bmask is set to 0 along lateral domain boundaries
116      !!      even IF nperio is not zero.
117      !!
118      !!      In case of open boundaries (lk_bdy=T):
119      !!        - tmask is set to 1 on the points to be computed bay the open
120      !!          boundaries routines.
121      !!        - bmask is  set to 0 on the open boundaries.
122      !!
123      !! ** Action :   tmask    : land/ocean mask at t-point (=0. or 1.)
124      !!               umask    : land/ocean mask at u-point (=0. or 1.)
125      !!               vmask    : land/ocean mask at v-point (=0. or 1.)
126      !!               fmask    : land/ocean mask at f-point (=0. or 1.)
127      !!                          =rn_shlat along lateral boundaries
128      !!               bmask    : land/ocean mask at barotropic stream
129      !!                          function point (=0. or 1.) and set to 0 along lateral boundaries
130      !!               tmask_i  : interior ocean mask
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)     ! surface
287      wumask(:,:,1) = umask(:,:,1)
288      wvmask(:,:,1) = vmask(:,:,1)
289      DO jk = 2, jpk                   ! interior values
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      END DO
376      !
377      IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN   ! ORCA_R2 configuration
378         !                                                 ! Increased lateral friction near of some straits
379         !                                ! Gibraltar strait  : partial slip (fmask=0.5)
380         ij0 = 101   ;   ij1 = 101
381         ii0 = 139   ;   ii1 = 140   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  0.5_wp
382         ij0 = 102   ;   ij1 = 102
383         ii0 = 139   ;   ii1 = 140   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  0.5_wp
384         !
385         !                                ! Bab el Mandeb : partial slip (fmask=1)
386         ij0 =  87   ;   ij1 =  88
387         ii0 = 160   ;   ii1 = 160   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  1._wp
388         ij0 =  88   ;   ij1 =  88
389         ii0 = 159   ;   ii1 = 159   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  1._wp
390         !
391         !                                ! Danish straits  : strong slip (fmask > 2)
392! We keep this as an example but it is instable in this case
393!         ij0 = 115   ;   ij1 = 115
394!         ii0 = 145   ;   ii1 = 146   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 4._wp
395!         ij0 = 116   ;   ij1 = 116
396!         ii0 = 145   ;   ii1 = 146   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 4._wp
397         !
398      ENDIF
399      !
400      IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN   ! ORCA R1 configuration
401         !                                                 ! Increased lateral friction near of some straits
402         ! This dirty section will be suppressed by simplification process:
403         ! all this will come back in input files
404         ! Currently these hard-wired indices relate to configuration with
405         ! extend grid (jpjglo=332)
406         !
407         isrow = 332 - jpjglo
408         !
409         IF(lwp) WRITE(numout,*)
410         IF(lwp) WRITE(numout,*) '   orca_r1: increase friction near the following straits : '
411         IF(lwp) WRITE(numout,*) '      Gibraltar '
412         ii0 = 282           ;   ii1 = 283        ! Gibraltar Strait
413         ij0 = 241 - isrow   ;   ij1 = 241 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp 
414
415         IF(lwp) WRITE(numout,*) '      Bhosporus '
416         ii0 = 314           ;   ii1 = 315        ! Bhosporus Strait
417         ij0 = 248 - isrow   ;   ij1 = 248 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp 
418
419         IF(lwp) WRITE(numout,*) '      Makassar (Top) '
420         ii0 =  48           ;   ii1 =  48        ! Makassar Strait (Top)
421         ij0 = 189 - isrow   ;   ij1 = 190 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp 
422
423         IF(lwp) WRITE(numout,*) '      Lombok '
424         ii0 =  44           ;   ii1 =  44        ! Lombok Strait
425         ij0 = 164 - isrow   ;   ij1 = 165 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp 
426
427         IF(lwp) WRITE(numout,*) '      Ombai '
428         ii0 =  53           ;   ii1 =  53        ! Ombai Strait
429         ij0 = 164 - isrow   ;   ij1 = 165 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp 
430
431         IF(lwp) WRITE(numout,*) '      Timor Passage '
432         ii0 =  56           ;   ii1 =  56        ! Timor Passage
433         ij0 = 164 - isrow   ;   ij1 = 165 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp 
434
435         IF(lwp) WRITE(numout,*) '      West Halmahera '
436         ii0 =  58           ;   ii1 =  58        ! West Halmahera Strait
437         ij0 = 181 - isrow   ;   ij1 = 182 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp 
438
439         IF(lwp) WRITE(numout,*) '      East Halmahera '
440         ii0 =  55           ;   ii1 =  55        ! East Halmahera Strait
441         ij0 = 181 - isrow   ;   ij1 = 182 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp 
442         !
443      ENDIF
444      !
445      CALL lbc_lnk( fmask, 'F', 1._wp )      ! Lateral boundary conditions on fmask
446
447      ! CAUTION : The fmask may be further modified in dyn_vor_init ( dynvor.F90 )
448           
449      IF( nprint == 1 .AND. lwp ) THEN      ! Control print
450         imsk(:,:) = INT( tmask_i(:,:) )
451         WRITE(numout,*) ' tmask_i : '
452         CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1,   &
453               &                           1, jpj, 1, 1, numout)
454         WRITE (numout,*)
455         WRITE (numout,*) ' dommsk: tmask for each level'
456         WRITE (numout,*) ' ----------------------------'
457         DO jk = 1, jpk
458            imsk(:,:) = INT( tmask(:,:,jk) )
459
460            WRITE(numout,*)
461            WRITE(numout,*) ' level = ',jk
462            CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1,   &
463               &                              1, jpj, 1, 1, numout)
464         END DO
465         WRITE(numout,*)
466         WRITE(numout,*) ' dom_msk: vmask for each level'
467         WRITE(numout,*) ' -----------------------------'
468         DO jk = 1, jpk
469            imsk(:,:) = INT( vmask(:,:,jk) )
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: fmask for each level'
477         WRITE(numout,*) ' -----------------------------'
478         DO jk = 1, jpk
479            imsk(:,:) = INT( fmask(:,:,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: bmask '
487         WRITE(numout,*) ' ---------------'
488         WRITE(numout,*)
489         imsk(:,:) = INT( bmask(:,:) )
490         CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1,   &
491            &                              1, jpj, 1, 1, numout )
492      ENDIF
493      !
494      CALL wrk_dealloc( jpi, jpj, imsk )
495      CALL wrk_dealloc( jpi, jpj, zwf  )
496      !
497      IF( nn_timing == 1 )  CALL timing_stop('dom_msk')
498      !
499   END SUBROUTINE dom_msk
500
501#if defined key_noslip_accurate
502   !!----------------------------------------------------------------------
503   !!   'key_noslip_accurate' :         accurate no-slip boundary condition
504   !!----------------------------------------------------------------------
505   
506   SUBROUTINE dom_msk_nsa
507      !!---------------------------------------------------------------------
508      !!                 ***  ROUTINE dom_msk_nsa  ***
509      !!
510      !! ** Purpose :
511      !!
512      !! ** Method  :
513      !!
514      !! ** Action :
515      !!----------------------------------------------------------------------
516      INTEGER  ::   ji, jj, jk, jl      ! dummy loop indices
517      INTEGER  ::   ine, inw, ins, inn, itest, ierror, iind, ijnd
518      REAL(wp) ::   zaa
519      !!---------------------------------------------------------------------
520      !
521      IF( nn_timing == 1 )  CALL timing_start('dom_msk_nsa')
522      !
523      IF(lwp) WRITE(numout,*)
524      IF(lwp) WRITE(numout,*) 'dom_msk_nsa : noslip accurate boundary condition'
525      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   using Schchepetkin and O Brian scheme'
526      IF( lk_mpp )   CALL ctl_stop( ' mpp version is not yet implemented' )
527
528      ! mask for second order calculation of vorticity
529      ! ----------------------------------------------
530      ! noslip boundary condition: fmask=1  at convex corner, store
531      ! index of straight coast meshes ( 'west', refering to a coast,
532      ! means west of the ocean, aso)
533     
534      DO jk = 1, jpk
535         DO jl = 1, 4
536            npcoa(jl,jk) = 0
537            DO ji = 1, 2*(jpi+jpj)
538               nicoa(ji,jl,jk) = 0
539               njcoa(ji,jl,jk) = 0
540            END DO
541         END DO
542      END DO
543     
544      IF( jperio == 2 ) THEN
545         WRITE(numout,*) ' '
546         WRITE(numout,*) ' symetric boundary conditions need special'
547         WRITE(numout,*) ' treatment not implemented. we stop.'
548         STOP
549      ENDIF
550     
551      ! convex corners
552     
553      DO jk = 1, jpkm1
554         DO jj = 1, jpjm1
555            DO ji = 1, jpim1
556               zaa = tmask(ji  ,jj,jk) + tmask(ji  ,jj+1,jk)   &
557                  &+ tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk)
558               IF( ABS(zaa-3._wp) <= 0.1_wp )   fmask(ji,jj,jk) = 1._wp
559            END DO
560         END DO
561      END DO
562
563      ! north-south straight coast
564
565      DO jk = 1, jpkm1
566         inw = 0
567         ine = 0
568         DO jj = 2, jpjm1
569            DO ji = 2, jpim1
570               zaa = tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk)
571               IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN
572                  inw = inw + 1
573                  nicoa(inw,1,jk) = ji
574                  njcoa(inw,1,jk) = jj
575                  IF( nprint == 1 ) WRITE(numout,*) ' west  : ', jk, inw, ji, jj
576               ENDIF
577               zaa = tmask(ji,jj,jk) + tmask(ji,jj+1,jk)
578               IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN
579                  ine = ine + 1
580                  nicoa(ine,2,jk) = ji
581                  njcoa(ine,2,jk) = jj
582                  IF( nprint == 1 ) WRITE(numout,*) ' east  : ', jk, ine, ji, jj
583               ENDIF
584            END DO
585         END DO
586         npcoa(1,jk) = inw
587         npcoa(2,jk) = ine
588      END DO
589
590      ! west-east straight coast
591
592      DO jk = 1, jpkm1
593         ins = 0
594         inn = 0
595         DO jj = 2, jpjm1
596            DO ji =2, jpim1
597               zaa = tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk)
598               IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN
599                  ins = ins + 1
600                  nicoa(ins,3,jk) = ji
601                  njcoa(ins,3,jk) = jj
602                  IF( nprint == 1 ) WRITE(numout,*) ' south : ', jk, ins, ji, jj
603               ENDIF
604               zaa = tmask(ji+1,jj,jk) + tmask(ji,jj,jk)
605               IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN
606                  inn = inn + 1
607                  nicoa(inn,4,jk) = ji
608                  njcoa(inn,4,jk) = jj
609                  IF( nprint == 1 ) WRITE(numout,*) ' north : ', jk, inn, ji, jj
610               ENDIF
611            END DO
612         END DO
613         npcoa(3,jk) = ins
614         npcoa(4,jk) = inn
615      END DO
616
617      itest = 2 * ( jpi + jpj )
618      DO jk = 1, jpk
619         IF( npcoa(1,jk) > itest .OR. npcoa(2,jk) > itest .OR.   &
620             npcoa(3,jk) > itest .OR. npcoa(4,jk) > itest ) THEN
621           
622            WRITE(ctmp1,*) ' level jk = ',jk
623            WRITE(ctmp2,*) ' straight coast index arraies are too small.:'
624            WRITE(ctmp3,*) ' npe, npw, nps, npn = ', npcoa(1,jk), npcoa(2,jk),   &
625                &                                     npcoa(3,jk), npcoa(4,jk)
626            WRITE(ctmp4,*) ' 2*(jpi+jpj) = ',itest,'. we stop.'
627            CALL ctl_stop( ctmp1, ctmp2, ctmp3, ctmp4 )
628        ENDIF
629      END DO
630
631      ierror = 0
632      iind = 0
633      ijnd = 0
634      IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 )   iind = 2
635      IF( nperio == 3 .OR. nperio == 4 .OR. nperio == 5 .OR. nperio == 6 )   ijnd = 2
636      DO jk = 1, jpk
637         DO jl = 1, npcoa(1,jk)
638            IF( nicoa(jl,1,jk)+3 > jpi+iind ) THEN
639               ierror = ierror+1
640               icoord(ierror,1) = nicoa(jl,1,jk)
641               icoord(ierror,2) = njcoa(jl,1,jk)
642               icoord(ierror,3) = jk
643            ENDIF
644         END DO
645         DO jl = 1, npcoa(2,jk)
646            IF(nicoa(jl,2,jk)-2 < 1-iind ) THEN
647               ierror = ierror + 1
648               icoord(ierror,1) = nicoa(jl,2,jk)
649               icoord(ierror,2) = njcoa(jl,2,jk)
650               icoord(ierror,3) = jk
651            ENDIF
652         END DO
653         DO jl = 1, npcoa(3,jk)
654            IF( njcoa(jl,3,jk)+3 > jpj+ijnd ) THEN
655               ierror = ierror + 1
656               icoord(ierror,1) = nicoa(jl,3,jk)
657               icoord(ierror,2) = njcoa(jl,3,jk)
658               icoord(ierror,3) = jk
659            ENDIF
660         END DO
661         DO jl = 1, npcoa(4,jk)
662            IF( njcoa(jl,4,jk)-2 < 1) THEN
663               ierror=ierror + 1
664               icoord(ierror,1) = nicoa(jl,4,jk)
665               icoord(ierror,2) = njcoa(jl,4,jk)
666               icoord(ierror,3) = jk
667            ENDIF
668         END DO
669      END DO
670     
671      IF( ierror > 0 ) THEN
672         IF(lwp) WRITE(numout,*)
673         IF(lwp) WRITE(numout,*) '              Problem on lateral conditions'
674         IF(lwp) WRITE(numout,*) '                 Bad marking off at points:'
675         DO jl = 1, ierror
676            IF(lwp) WRITE(numout,*) 'Level:',icoord(jl,3),   &
677               &                  '  Point(',icoord(jl,1),',',icoord(jl,2),')'
678         END DO
679         CALL ctl_stop( 'We stop...' )
680      ENDIF
681      !
682      IF( nn_timing == 1 )  CALL timing_stop('dom_msk_nsa')
683      !
684   END SUBROUTINE dom_msk_nsa
685
686#else
687   !!----------------------------------------------------------------------
688   !!   Default option :                                      Empty routine
689   !!----------------------------------------------------------------------
690   SUBROUTINE dom_msk_nsa       
691   END SUBROUTINE dom_msk_nsa
692#endif
693   
694   !!======================================================================
695END MODULE dommsk
Note: See TracBrowser for help on using the repository browser.