New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
dommsk.F90 in branches/UKMO/dev_isf_flx_UKESM_r9321/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

source: branches/UKMO/dev_isf_flx_UKESM_r9321/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90 @ 10026

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

Compute isf heat content flx using interface S in of far field S + add isfmask to simplify isf masking

File size: 32.3 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
221      ! in case of iceshelf definition of iceshelf surface mask (1 under
222      ! iceshelf and 0 elsewhere
223      ssisfmask = (1._wp - tmask(:,:,1)) * ssmask(:,:)
224
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('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         CALL ctl_stop('STOP', 'NEMO abort from dom_msk_nsa')
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   !!======================================================================
726END MODULE dommsk
Note: See TracBrowser for help on using the repository browser.