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

source: branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90 @ 11559

Last change on this file since 11559 was 11559, checked in by mattmartin, 5 years ago

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