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

source: branches/UKMO/dev_r7651_GO6pck_shlat2d/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90 @ 7678

Last change on this file since 7678 was 7678, checked in by mathiot, 7 years ago

add option to read shlat2d

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