source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90 @ 8280

Last change on this file since 8280 was 8280, checked in by timgraham, 3 years ago

331: Merge of MEDUSA stable branch and HadGEM3 coupling branches into GO6 package branch.

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