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 @ 4400

Last change on this file since 4400 was 3211, checked in by spickles2, 12 years ago

Stephen Pickles, 11 Dec 2011

Commit to bring the rest of the DCSE NEMO development branch
in line with the latest development version. This includes
array index re-ordering of all OPA_SRC/.

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