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

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

#2050 first version. Compiled OK in moci test suite

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