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

source: branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90 @ 5098

Last change on this file since 5098 was 5098, checked in by mathiot, 9 years ago

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