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

source: branches/UKMO/dev_r5518_GO6_package_text_diagnostics/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90 @ 10774

Last change on this file since 10774 was 10774, checked in by andmirek, 5 years ago

GMED 450 add flush after prints

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