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

source: branches/UKMO/test_moci_test_suite_namelist_read/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90 @ 9383

Last change on this file since 9383 was 9383, checked in by andmirek, 6 years ago

#2050 fixes and changes

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