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

source: branches/UKMO/dev_isf_remapping_UKESM_GO6package_r9314/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90

Last change on this file was 11248, checked in by mathiot, 5 years ago

add condition to fill isolated grid point in the bathymetry in zgr_isf to ensure compatibility with ice shelf draft and properly mask bathy,risfdep,mbathy,misfdep after filling subglacial lakes

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