source: branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90 @ 4666

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

#1331 : add ISOMIP config files + ice shelf code

  • Property svn:keywords set to Id
File size: 29.8 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      lmask(:,:)=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( micedep(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(:,:) = lmask(:,:)            ! (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)  = lmask(ji,jj) * lmask(ji+1,jj  )  * MIN(1._wp,SUM(umask(ji,jj,:)))
269            vmask_i(ji,jj)  = lmask(ji,jj) * lmask(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) =  lmask(ji,jj  ) * lmask(ji+1,jj  )   &
273               &            * lmask(ji,jj+1) * lmask(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
284      ! 4. ocean/land mask for the elliptic equation
285      ! --------------------------------------------
286      bmask(:,:) = lmask(:,:)       ! elliptic equation is written at t-point
287      !
288      !                               ! Boundary conditions
289      !                                    ! cyclic east-west : bmask must be set to 0. on rows 1 and jpi
290      IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN
291         bmask( 1 ,:) = 0._wp
292         bmask(jpi,:) = 0._wp
293      ENDIF
294      IF( nperio == 2 ) THEN               ! south symmetric :  bmask must be set to 0. on row 1
295         bmask(:, 1 ) = 0._wp
296      ENDIF
297      !                                    ! north fold :
298      IF( nperio == 3 .OR. nperio == 4 ) THEN   ! T-pt pivot : bmask set to 0. on row jpj and on half jpjglo-1 row
299         DO ji = 1, jpi                     
300            ii = ji + nimpp - 1
301            bmask(ji,jpj-1) = bmask(ji,jpj-1) * tpol(ii)
302            bmask(ji,jpj  ) = 0._wp
303         END DO
304      ENDIF
305      IF( nperio == 5 .OR. nperio == 6 ) THEN   ! F-pt pivot and T-pt elliptic eq. : bmask set to 0. on row jpj
306         bmask(:,jpj) = 0._wp
307      ENDIF
308      !
309      IF( lk_mpp ) THEN                    ! mpp specificities
310         !                                      ! bmask is set to zero on the overlap region
311         IF( nbondi /= -1 .AND. nbondi /= 2 )   bmask(  1 :jpreci,:) = 0._wp
312         IF( nbondi /=  1 .AND. nbondi /= 2 )   bmask(nlci:jpi   ,:) = 0._wp
313         IF( nbondj /= -1 .AND. nbondj /= 2 )   bmask(:,  1 :jprecj) = 0._wp
314         IF( nbondj /=  1 .AND. nbondj /= 2 )   bmask(:,nlcj:jpj   ) = 0._wp
315         !
316         IF( npolj == 3 .OR. npolj == 4 ) THEN  ! north fold : bmask must be set to 0. on rows jpj-1 and jpj
317            DO ji = 1, nlci
318               ii = ji + nimpp - 1
319               bmask(ji,nlcj-1) = bmask(ji,nlcj-1) * tpol(ii)
320               bmask(ji,nlcj  ) = 0._wp
321            END DO
322         ENDIF
323         IF( npolj == 5 .OR. npolj == 6 ) THEN  ! F-pt pivot and T-pt elliptic eq. : bmask set to 0. on row jpj
324            DO ji = 1, nlci
325               bmask(ji,nlcj  ) = 0._wp
326            END DO
327         ENDIF
328      ENDIF
329
330
331      ! mask for second order calculation of vorticity
332      ! ----------------------------------------------
333      CALL dom_msk_nsa
334
335     
336      ! Lateral boundary conditions on velocity (modify fmask)
337      ! ---------------------------------------     
338      DO jk = 1, jpk
339         zwf(:,:) = fmask(:,:,jk)         
340         DO jj = 2, jpjm1
341            DO ji = fs_2, fs_jpim1   ! vector opt.
342               IF( fmask(ji,jj,jk) == 0._wp ) THEN
343                  fmask(ji,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,jj), zwf(ji,jj+1),   &
344                     &                                           zwf(ji-1,jj), zwf(ji,jj-1)  )  )
345               ENDIF
346            END DO
347         END DO
348         DO jj = 2, jpjm1
349            IF( fmask(1,jj,jk) == 0._wp ) THEN
350               fmask(1  ,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1) ) )
351            ENDIF
352            IF( fmask(jpi,jj,jk) == 0._wp ) THEN
353               fmask(jpi,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) )
354            ENDIF
355         END DO         
356         DO ji = 2, jpim1
357            IF( fmask(ji,1,jk) == 0._wp ) THEN
358               fmask(ji, 1 ,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1) ) )
359            ENDIF
360            IF( fmask(ji,jpj,jk) == 0._wp ) THEN
361               fmask(ji,jpj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpjm1) ) )
362            ENDIF
363         END DO
364      END DO
365      !
366      IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN   ! ORCA_R2 configuration
367         !                                                 ! Increased lateral friction near of some straits
368         IF( nn_cla == 0 ) THEN
369            !                                ! Gibraltar strait  : partial slip (fmask=0.5)
370            ij0 = 101   ;   ij1 = 101
371            ii0 = 139   ;   ii1 = 140   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  0.5_wp
372            ij0 = 102   ;   ij1 = 102
373            ii0 = 139   ;   ii1 = 140   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  0.5_wp
374            !
375            !                                ! Bab el Mandeb : partial slip (fmask=1)
376            ij0 =  87   ;   ij1 =  88
377            ii0 = 160   ;   ii1 = 160   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  1._wp
378            ij0 =  88   ;   ij1 =  88
379            ii0 = 159   ;   ii1 = 159   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  1._wp
380            !
381         ENDIF
382         !                                ! Danish straits  : strong slip (fmask > 2)
383! We keep this as an example but it is instable in this case
384!         ij0 = 115   ;   ij1 = 115
385!         ii0 = 145   ;   ii1 = 146   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 4._wp
386!         ij0 = 116   ;   ij1 = 116
387!         ii0 = 145   ;   ii1 = 146   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 4._wp
388         !
389      ENDIF
390      !
391      IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN   ! ORCA R1 configuration
392         !                                                 ! Increased lateral friction near of some straits
393         IF(lwp) WRITE(numout,*)
394         IF(lwp) WRITE(numout,*) '   orca_r1: increase friction near the following straits : '
395         IF(lwp) WRITE(numout,*) '      Gibraltar '
396         ii0 = 283   ;   ii1 = 284        ! Gibraltar Strait
397         ij0 = 200   ;   ij1 = 200   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) =  2._wp 
398
399         IF(lwp) WRITE(numout,*) '      Bhosporus '
400         ii0 = 314   ;   ii1 = 315        ! Bhosporus Strait
401         ij0 = 208   ;   ij1 = 208   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) =  2._wp 
402
403         IF(lwp) WRITE(numout,*) '      Makassar (Top) '
404         ii0 =  48   ;   ii1 =  48        ! Makassar Strait (Top)
405         ij0 = 149   ;   ij1 = 150   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) =  3._wp 
406
407         IF(lwp) WRITE(numout,*) '      Lombok '
408         ii0 =  44   ;   ii1 =  44        ! Lombok Strait
409         ij0 = 124   ;   ij1 = 125   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) =  2._wp 
410
411         IF(lwp) WRITE(numout,*) '      Ombai '
412         ii0 =  53   ;   ii1 =  53        ! Ombai Strait
413         ij0 = 124   ;   ij1 = 125   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) = 2._wp 
414
415         IF(lwp) WRITE(numout,*) '      Timor Passage '
416         ii0 =  56   ;   ii1 =  56        ! Timor Passage
417         ij0 = 124   ;   ij1 = 125   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) = 2._wp 
418
419         IF(lwp) WRITE(numout,*) '      West Halmahera '
420         ii0 =  58   ;   ii1 =  58        ! West Halmahera Strait
421         ij0 = 141   ;   ij1 = 142   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) = 3._wp 
422
423         IF(lwp) WRITE(numout,*) '      East Halmahera '
424         ii0 =  55   ;   ii1 =  55        ! East Halmahera Strait
425         ij0 = 141   ;   ij1 = 142   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) = 3._wp 
426         !
427      ENDIF
428      !
429      CALL lbc_lnk( fmask, 'F', 1._wp )      ! Lateral boundary conditions on fmask
430
431      ! CAUTION : The fmask may be further modified in dyn_vor_init ( dynvor.F90 )
432           
433      IF( nprint == 1 .AND. lwp ) THEN      ! Control print
434         imsk(:,:) = INT( tmask_i(:,:) )
435         WRITE(numout,*) ' tmask_i : '
436         CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1,   &
437               &                           1, jpj, 1, 1, numout)
438         WRITE (numout,*)
439         WRITE (numout,*) ' dommsk: tmask for each level'
440         WRITE (numout,*) ' ----------------------------'
441         DO jk = 1, jpk
442            imsk(:,:) = INT( tmask(:,:,jk) )
443
444            WRITE(numout,*)
445            WRITE(numout,*) ' level = ',jk
446            CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1,   &
447               &                              1, jpj, 1, 1, numout)
448         END DO
449         WRITE(numout,*)
450         WRITE(numout,*) ' dom_msk: vmask for each level'
451         WRITE(numout,*) ' -----------------------------'
452         DO jk = 1, jpk
453            imsk(:,:) = INT( vmask(:,:,jk) )
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: fmask for each level'
461         WRITE(numout,*) ' -----------------------------'
462         DO jk = 1, jpk
463            imsk(:,:) = INT( fmask(:,:,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: bmask '
471         WRITE(numout,*) ' ---------------'
472         WRITE(numout,*)
473         imsk(:,:) = INT( bmask(:,:) )
474         CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1,   &
475            &                              1, jpj, 1, 1, numout )
476      ENDIF
477      !
478      CALL wrk_dealloc( jpi, jpj, imsk )
479      CALL wrk_dealloc( jpi, jpj, zwf  )
480      !
481      IF( nn_timing == 1 )  CALL timing_stop('dom_msk')
482      !
483   END SUBROUTINE dom_msk
484
485#if defined key_noslip_accurate
486   !!----------------------------------------------------------------------
487   !!   'key_noslip_accurate' :         accurate no-slip boundary condition
488   !!----------------------------------------------------------------------
489   
490   SUBROUTINE dom_msk_nsa
491      !!---------------------------------------------------------------------
492      !!                 ***  ROUTINE dom_msk_nsa  ***
493      !!
494      !! ** Purpose :
495      !!
496      !! ** Method  :
497      !!
498      !! ** Action :
499      !!----------------------------------------------------------------------
500      INTEGER  ::   ji, jj, jk, jl      ! dummy loop indices
501      INTEGER  ::   ine, inw, ins, inn, itest, ierror, iind, ijnd
502      REAL(wp) ::   zaa
503      !!---------------------------------------------------------------------
504      !
505      IF( nn_timing == 1 )  CALL timing_start('dom_msk_nsa')
506      !
507      IF(lwp) WRITE(numout,*)
508      IF(lwp) WRITE(numout,*) 'dom_msk_nsa : noslip accurate boundary condition'
509      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   using Schchepetkin and O Brian scheme'
510      IF( lk_mpp )   CALL ctl_stop( ' mpp version is not yet implemented' )
511
512      ! mask for second order calculation of vorticity
513      ! ----------------------------------------------
514      ! noslip boundary condition: fmask=1  at convex corner, store
515      ! index of straight coast meshes ( 'west', refering to a coast,
516      ! means west of the ocean, aso)
517     
518      DO jk = 1, jpk
519         DO jl = 1, 4
520            npcoa(jl,jk) = 0
521            DO ji = 1, 2*(jpi+jpj)
522               nicoa(ji,jl,jk) = 0
523               njcoa(ji,jl,jk) = 0
524            END DO
525         END DO
526      END DO
527     
528      IF( jperio == 2 ) THEN
529         WRITE(numout,*) ' '
530         WRITE(numout,*) ' symetric boundary conditions need special'
531         WRITE(numout,*) ' treatment not implemented. we stop.'
532         STOP
533      ENDIF
534     
535      ! convex corners
536     
537      DO jk = 1, jpkm1
538         DO jj = 1, jpjm1
539            DO ji = 1, jpim1
540               zaa = tmask(ji  ,jj,jk) + tmask(ji  ,jj+1,jk)   &
541                  &+ tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk)
542               IF( ABS(zaa-3._wp) <= 0.1_wp )   fmask(ji,jj,jk) = 1._wp
543            END DO
544         END DO
545      END DO
546
547      ! north-south straight coast
548
549      DO jk = 1, jpkm1
550         inw = 0
551         ine = 0
552         DO jj = 2, jpjm1
553            DO ji = 2, jpim1
554               zaa = tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk)
555               IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN
556                  inw = inw + 1
557                  nicoa(inw,1,jk) = ji
558                  njcoa(inw,1,jk) = jj
559                  IF( nprint == 1 ) WRITE(numout,*) ' west  : ', jk, inw, ji, jj
560               ENDIF
561               zaa = tmask(ji,jj,jk) + tmask(ji,jj+1,jk)
562               IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN
563                  ine = ine + 1
564                  nicoa(ine,2,jk) = ji
565                  njcoa(ine,2,jk) = jj
566                  IF( nprint == 1 ) WRITE(numout,*) ' east  : ', jk, ine, ji, jj
567               ENDIF
568            END DO
569         END DO
570         npcoa(1,jk) = inw
571         npcoa(2,jk) = ine
572      END DO
573
574      ! west-east straight coast
575
576      DO jk = 1, jpkm1
577         ins = 0
578         inn = 0
579         DO jj = 2, jpjm1
580            DO ji =2, jpim1
581               zaa = tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk)
582               IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN
583                  ins = ins + 1
584                  nicoa(ins,3,jk) = ji
585                  njcoa(ins,3,jk) = jj
586                  IF( nprint == 1 ) WRITE(numout,*) ' south : ', jk, ins, ji, jj
587               ENDIF
588               zaa = tmask(ji+1,jj,jk) + tmask(ji,jj,jk)
589               IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN
590                  inn = inn + 1
591                  nicoa(inn,4,jk) = ji
592                  njcoa(inn,4,jk) = jj
593                  IF( nprint == 1 ) WRITE(numout,*) ' north : ', jk, inn, ji, jj
594               ENDIF
595            END DO
596         END DO
597         npcoa(3,jk) = ins
598         npcoa(4,jk) = inn
599      END DO
600
601      itest = 2 * ( jpi + jpj )
602      DO jk = 1, jpk
603         IF( npcoa(1,jk) > itest .OR. npcoa(2,jk) > itest .OR.   &
604             npcoa(3,jk) > itest .OR. npcoa(4,jk) > itest ) THEN
605           
606            WRITE(ctmp1,*) ' level jk = ',jk
607            WRITE(ctmp2,*) ' straight coast index arraies are too small.:'
608            WRITE(ctmp3,*) ' npe, npw, nps, npn = ', npcoa(1,jk), npcoa(2,jk),   &
609                &                                     npcoa(3,jk), npcoa(4,jk)
610            WRITE(ctmp4,*) ' 2*(jpi+jpj) = ',itest,'. we stop.'
611            CALL ctl_stop( ctmp1, ctmp2, ctmp3, ctmp4 )
612        ENDIF
613      END DO
614
615      ierror = 0
616      iind = 0
617      ijnd = 0
618      IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 )   iind = 2
619      IF( nperio == 3 .OR. nperio == 4 .OR. nperio == 5 .OR. nperio == 6 )   ijnd = 2
620      DO jk = 1, jpk
621         DO jl = 1, npcoa(1,jk)
622            IF( nicoa(jl,1,jk)+3 > jpi+iind ) THEN
623               ierror = ierror+1
624               icoord(ierror,1) = nicoa(jl,1,jk)
625               icoord(ierror,2) = njcoa(jl,1,jk)
626               icoord(ierror,3) = jk
627            ENDIF
628         END DO
629         DO jl = 1, npcoa(2,jk)
630            IF(nicoa(jl,2,jk)-2 < 1-iind ) THEN
631               ierror = ierror + 1
632               icoord(ierror,1) = nicoa(jl,2,jk)
633               icoord(ierror,2) = njcoa(jl,2,jk)
634               icoord(ierror,3) = jk
635            ENDIF
636         END DO
637         DO jl = 1, npcoa(3,jk)
638            IF( njcoa(jl,3,jk)+3 > jpj+ijnd ) THEN
639               ierror = ierror + 1
640               icoord(ierror,1) = nicoa(jl,3,jk)
641               icoord(ierror,2) = njcoa(jl,3,jk)
642               icoord(ierror,3) = jk
643            ENDIF
644         END DO
645         DO jl = 1, npcoa(4,jk)
646            IF( njcoa(jl,4,jk)-2 < 1) THEN
647               ierror=ierror + 1
648               icoord(ierror,1) = nicoa(jl,4,jk)
649               icoord(ierror,2) = njcoa(jl,4,jk)
650               icoord(ierror,3) = jk
651            ENDIF
652         END DO
653      END DO
654     
655      IF( ierror > 0 ) THEN
656         IF(lwp) WRITE(numout,*)
657         IF(lwp) WRITE(numout,*) '              Problem on lateral conditions'
658         IF(lwp) WRITE(numout,*) '                 Bad marking off at points:'
659         DO jl = 1, ierror
660            IF(lwp) WRITE(numout,*) 'Level:',icoord(jl,3),   &
661               &                  '  Point(',icoord(jl,1),',',icoord(jl,2),')'
662         END DO
663         CALL ctl_stop( 'We stop...' )
664      ENDIF
665      !
666      IF( nn_timing == 1 )  CALL timing_stop('dom_msk_nsa')
667      !
668   END SUBROUTINE dom_msk_nsa
669
670#else
671   !!----------------------------------------------------------------------
672   !!   Default option :                                      Empty routine
673   !!----------------------------------------------------------------------
674   SUBROUTINE dom_msk_nsa       
675   END SUBROUTINE dom_msk_nsa
676#endif
677   
678   !!======================================================================
679END MODULE dommsk
Note: See TracBrowser for help on using the repository browser.