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

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90 @ 4427

Last change on this file since 4427 was 4427, checked in by trackstand2, 10 years ago

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