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 @ 9513

Last change on this file since 9513 was 9513, checked in by mathiot, 6 years ago

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