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

source: branches/UKMO/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90 @ 11134

Last change on this file since 11134 was 11134, checked in by jcastill, 5 years ago

Full set of changes as in the original branch

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