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

source: branches/2013/dev_r3853_CNRS9_ConfSetting/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90 @ 3901

Last change on this file since 3901 was 3901, checked in by clevy, 11 years ago

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