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/UKMO/antarctic_partial_slip/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

source: branches/UKMO/antarctic_partial_slip/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90 @ 5958

Last change on this file since 5958 was 5958, checked in by davestorkey, 9 years ago

UKMO/antarctica_partial_slip branch: remove SVN keywords.
File size: 30.7 KB
Line 
1MODULE dommsk
2   !!======================================================================
3   !!                       ***  MODULE dommsk   ***
4   !! Ocean initialization : domain land/sea mask
5   !!======================================================================
6   !! History :  OPA  ! 1987-07  (G. Madec)  Original code
7   !!            6.0  ! 1993-03  (M. Guyon)  symetrical conditions (M. Guyon)
8   !!            7.0  ! 1996-01  (G. Madec)  suppression of common work arrays
9   !!             -   ! 1996-05  (G. Madec)  mask computed from tmask and sup-
10   !!                 !                      pression of the double computation of bmask
11   !!            8.0  ! 1997-02  (G. Madec)  mesh information put in domhgr.F
12   !!            8.1  ! 1997-07  (G. Madec)  modification of mbathy and fmask
13   !!             -   ! 1998-05  (G. Roullet)  free surface
14   !!            8.2  ! 2000-03  (G. Madec)  no slip accurate
15   !!             -   ! 2001-09  (J.-M. Molines)  Open boundaries
16   !!   NEMO     1.0  ! 2002-08  (G. Madec)  F90: Free form and module
17   !!             -   ! 2005-11  (V. Garnier) Surface pressure gradient organization
18   !!            3.2  ! 2009-07  (R. Benshila) Suppression of rigid-lid option
19   !!----------------------------------------------------------------------
20
21   !!----------------------------------------------------------------------
22   !!   dom_msk        : compute land/ocean mask
23   !!   dom_msk_nsa    : update land/ocean mask when no-slip accurate option is used.
24   !!----------------------------------------------------------------------
25   USE oce             ! ocean dynamics and tracers
26   USE dom_oce         ! ocean space and time domain
27   USE in_out_manager  ! I/O manager
28   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
29   USE lib_mpp
30   USE dynspg_oce      ! choice/control of key cpp for surface pressure gradient
31   USE wrk_nemo        ! Memory allocation
32   USE timing          ! Timing
33
34   IMPLICIT NONE
35   PRIVATE
36
37   PUBLIC   dom_msk         ! routine called by inidom.F90
38   PUBLIC   dom_msk_alloc   ! routine called by nemogcm.F90
39
40   !                            !!* Namelist namlbc : lateral boundary condition *
41   REAL(wp)        :: rn_shlat   ! type of lateral boundary condition on velocity
42   LOGICAL, PUBLIC :: ln_vorlat  !  consistency of vorticity boundary condition
43   !                                            with analytical eqs.
44
45
46   INTEGER, ALLOCATABLE, SAVE, DIMENSION(:,:) ::  icoord ! Workspace for dom_msk_nsa()
47
48   !! * Substitutions
49#  include "vectopt_loop_substitute.h90"
50   !!----------------------------------------------------------------------
51   !! NEMO/OPA 3.2 , LODYC-IPSL  (2009)
52   !! $Id$
53   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
54   !!----------------------------------------------------------------------
55CONTAINS
56   
57   INTEGER FUNCTION dom_msk_alloc()
58      !!---------------------------------------------------------------------
59      !!                 ***  FUNCTION dom_msk_alloc  ***
60      !!---------------------------------------------------------------------
61      dom_msk_alloc = 0
62#if defined key_noslip_accurate
63      ALLOCATE(icoord(jpi*jpj*jpk,3), STAT=dom_msk_alloc)
64#endif
65      IF( dom_msk_alloc /= 0 )   CALL ctl_warn('dom_msk_alloc: failed to allocate icoord array')
66      !
67   END FUNCTION dom_msk_alloc
68
69
70   SUBROUTINE dom_msk
71      !!---------------------------------------------------------------------
72      !!                 ***  ROUTINE dom_msk  ***
73      !!
74      !! ** Purpose :   Compute land/ocean mask arrays at tracer points, hori-
75      !!      zontal velocity points (u & v), vorticity points (f) and baro-
76      !!      tropic stream function  points (b).
77      !!
78      !! ** Method  :   The ocean/land mask is computed from the basin bathy-
79      !!      metry in level (mbathy) which is defined or read in dommba.
80      !!      mbathy equals 0 over continental T-point
81      !!      and the number of ocean level over the ocean.
82      !!
83      !!      At a given position (ji,jj,jk) the ocean/land mask is given by:
84      !!      t-point : 0. IF mbathy( ji ,jj) =< 0
85      !!                1. IF mbathy( ji ,jj) >= jk
86      !!      u-point : 0. IF mbathy( ji ,jj)  or mbathy(ji+1, jj ) =< 0
87      !!                1. IF mbathy( ji ,jj) and mbathy(ji+1, jj ) >= jk.
88      !!      v-point : 0. IF mbathy( ji ,jj)  or mbathy( ji ,jj+1) =< 0
89      !!                1. IF mbathy( ji ,jj) and mbathy( ji ,jj+1) >= jk.
90      !!      f-point : 0. IF mbathy( ji ,jj)  or mbathy( ji ,jj+1)
91      !!                   or mbathy(ji+1,jj)  or mbathy(ji+1,jj+1) =< 0
92      !!                1. IF mbathy( ji ,jj) and mbathy( ji ,jj+1)
93      !!                  and mbathy(ji+1,jj) and mbathy(ji+1,jj+1) >= jk.
94      !!      b-point : the same definition as for f-point of the first ocean
95      !!                level (surface level) but with 0 along coastlines.
96      !!      tmask_i : interior ocean mask at t-point, i.e. excluding duplicated
97      !!                rows/lines due to cyclic or North Fold boundaries as well
98      !!                as MPP halos.
99      !!
100      !!        The lateral friction is set through the value of fmask along
101      !!      the coast and topography. This value is defined by rn_shlat, a
102      !!      namelist parameter:
103      !!         rn_shlat = 0, free slip  (no shear along the coast)
104      !!         rn_shlat = 2, no slip  (specified zero velocity at the coast)
105      !!         0 < rn_shlat < 2, partial slip   | non-linear velocity profile
106      !!         2 < rn_shlat, strong slip        | in the lateral boundary layer
107      !!
108      !!      N.B. If nperio not equal to 0, the land/ocean mask arrays
109      !!      are defined with the proper value at lateral domain boundaries,
110      !!      but bmask. indeed, bmask defined the domain over which the
111      !!      barotropic stream function is computed. this domain cannot
112      !!      contain identical columns because the matrix associated with
113      !!      the barotropic stream function equation is then no more inverti-
114      !!      ble. therefore bmask is set to 0 along lateral domain boundaries
115      !!      even IF nperio is not zero.
116      !!
117      !!      In case of open boundaries (lk_bdy=T):
118      !!        - tmask is set to 1 on the points to be computed bay the open
119      !!          boundaries routines.
120      !!        - bmask is  set to 0 on the open boundaries.
121      !!
122      !! ** Action :   tmask    : land/ocean mask at t-point (=0. or 1.)
123      !!               umask    : land/ocean mask at u-point (=0. or 1.)
124      !!               vmask    : land/ocean mask at v-point (=0. or 1.)
125      !!               fmask    : land/ocean mask at f-point (=0. or 1.)
126      !!                          =rn_shlat along lateral boundaries
127      !!               bmask    : land/ocean mask at barotropic stream
128      !!                          function point (=0. or 1.) and set to 0 along lateral boundaries
129      !!               tmask_i  : interior ocean mask
130      !!----------------------------------------------------------------------
131      !
132      INTEGER  ::   ji, jj, jk      ! dummy loop indices
133      INTEGER  ::   iif, iil, ii0, ii1, ii   ! local integers
134      INTEGER  ::   ijf, ijl, ij0, ij1       !   -       -
135      INTEGER  ::   ios
136      INTEGER  ::   isrow                    ! index for ORCA1 starting row
137      INTEGER , POINTER, DIMENSION(:,:) ::  imsk
138      REAL(wp), POINTER, DIMENSION(:,:) ::  zwf
139      !!
140      NAMELIST/namlbc/ rn_shlat, ln_vorlat
141      !!---------------------------------------------------------------------
142      !
143      IF( nn_timing == 1 )  CALL timing_start('dom_msk')
144      !
145      CALL wrk_alloc( jpi, jpj, imsk )
146      CALL wrk_alloc( jpi, jpj, zwf  )
147      !
148      REWIND( numnam_ref )              ! Namelist namlbc in reference namelist : Lateral momentum boundary condition
149      READ  ( numnam_ref, namlbc, IOSTAT = ios, ERR = 901 )
150901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namlbc in reference namelist', lwp )
151
152      REWIND( numnam_cfg )              ! Namelist namlbc in configuration namelist : Lateral momentum boundary condition
153      READ  ( numnam_cfg, namlbc, IOSTAT = ios, ERR = 902 )
154902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namlbc in configuration namelist', lwp )
155      IF(lwm) WRITE ( numond, namlbc )
156     
157      IF(lwp) THEN                  ! control print
158         WRITE(numout,*)
159         WRITE(numout,*) 'dommsk : ocean mask '
160         WRITE(numout,*) '~~~~~~'
161         WRITE(numout,*) '   Namelist namlbc'
162         WRITE(numout,*) '      lateral momentum boundary cond.    rn_shlat  = ',rn_shlat
163         WRITE(numout,*) '      consistency with analytical form   ln_vorlat = ',ln_vorlat 
164      ENDIF
165
166      IF     (      rn_shlat == 0.               ) THEN   ;   IF(lwp) WRITE(numout,*) '   ocean lateral  free-slip '
167      ELSEIF (      rn_shlat == 2.               ) THEN   ;   IF(lwp) WRITE(numout,*) '   ocean lateral  no-slip '
168      ELSEIF ( 0. < rn_shlat .AND. rn_shlat < 2. ) THEN   ;   IF(lwp) WRITE(numout,*) '   ocean lateral  partial-slip '
169      ELSEIF ( 2. < rn_shlat                     ) THEN   ;   IF(lwp) WRITE(numout,*) '   ocean lateral  strong-slip '
170      ELSE
171         WRITE(ctmp1,*) ' rn_shlat is negative = ', rn_shlat
172         CALL ctl_stop( ctmp1 )
173      ENDIF
174
175      ! 1. Ocean/land mask at t-point (computed from mbathy)
176      ! -----------------------------
177      ! N.B. tmask has already the right boundary conditions since mbathy is ok
178      !
179      tmask(:,:,:) = 0._wp
180      DO jk = 1, jpk
181         DO jj = 1, jpj
182            DO ji = 1, jpi
183               IF( REAL( mbathy(ji,jj) - jk, wp ) + 0.1_wp >= 0._wp )   tmask(ji,jj,jk) = 1._wp
184            END DO 
185         END DO 
186      END DO 
187     
188      ! (ISF) define barotropic mask and mask the ice shelf point
189      ssmask(:,:)=tmask(:,:,1) ! at this stage ice shelf is not masked
190     
191      DO jk = 1, jpk
192         DO jj = 1, jpj
193            DO ji = 1, jpi
194               IF( REAL( misfdep(ji,jj) - jk, wp ) - 0.1_wp >= 0._wp )   THEN
195                  tmask(ji,jj,jk) = 0._wp
196               END IF
197            END DO 
198         END DO 
199      END DO 
200
201!!gm  ????
202#if defined key_zdfkpp
203      IF( cp_cfg == 'orca' ) THEN
204         IF( jp_cfg == 2 )   THEN       ! land point on Bab el Mandeb zonal section
205            ij0 =  87   ;   ij1 =  88
206            ii0 = 160   ;   ii1 = 161
207            tmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 0._wp
208         ELSE
209            IF(lwp) WRITE(numout,*)
210            IF(lwp) WRITE(numout,cform_war)
211            IF(lwp) WRITE(numout,*)
212            IF(lwp) WRITE(numout,*)'          A mask must be applied on Bab el Mandeb strait'
213            IF(lwp) WRITE(numout,*)'          in case of ORCAs configurations'
214            IF(lwp) WRITE(numout,*)'          This is a problem which is not yet solved'
215            IF(lwp) WRITE(numout,*)
216         ENDIF
217      ENDIF
218#endif
219!!gm end
220
221      ! Interior domain mask (used for global sum)
222      ! --------------------
223      tmask_i(:,:) = ssmask(:,:)            ! (ISH) tmask_i = 1 even on the ice shelf
224      iif = jpreci                         ! ???
225      iil = nlci - jpreci + 1
226      ijf = jprecj                         ! ???
227      ijl = nlcj - jprecj + 1
228
229      tmask_i( 1 :iif,   :   ) = 0._wp      ! first columns
230      tmask_i(iil:jpi,   :   ) = 0._wp      ! last  columns (including mpp extra columns)
231      tmask_i(   :   , 1 :ijf) = 0._wp      ! first rows
232      tmask_i(   :   ,ijl:jpj) = 0._wp      ! last  rows (including mpp extra rows)
233
234      ! north fold mask
235      ! ---------------
236      tpol(1:jpiglo) = 1._wp 
237      fpol(1:jpiglo) = 1._wp
238      IF( jperio == 3 .OR. jperio == 4 ) THEN      ! T-point pivot
239         tpol(jpiglo/2+1:jpiglo) = 0._wp
240         fpol(     1    :jpiglo) = 0._wp
241         IF( mjg(nlej) == jpjglo ) THEN                  ! only half of the nlcj-1 row
242            DO ji = iif+1, iil-1
243               tmask_i(ji,nlej-1) = tmask_i(ji,nlej-1) * tpol(mig(ji))
244            END DO
245         ENDIF
246      ENDIF
247      IF( jperio == 5 .OR. jperio == 6 ) THEN      ! F-point pivot
248         tpol(     1    :jpiglo) = 0._wp
249         fpol(jpiglo/2+1:jpiglo) = 0._wp
250      ENDIF
251
252      ! 2. Ocean/land mask at u-,  v-, and z-points (computed from tmask)
253      ! -------------------------------------------
254      DO jk = 1, jpk
255         DO jj = 1, jpjm1
256            DO ji = 1, fs_jpim1   ! vector loop
257               umask(ji,jj,jk) = tmask(ji,jj  ,jk) * tmask(ji+1,jj  ,jk)
258               vmask(ji,jj,jk) = tmask(ji,jj  ,jk) * tmask(ji  ,jj+1,jk)
259            END DO
260            DO ji = 1, jpim1      ! NO vector opt.
261               fmask(ji,jj,jk) = tmask(ji,jj  ,jk) * tmask(ji+1,jj  ,jk)   &
262                  &            * tmask(ji,jj+1,jk) * tmask(ji+1,jj+1,jk)
263            END DO
264         END DO
265      END DO
266      ! (ISF) MIN(1,SUM(umask)) is here to check if you have effectively at least 1 wet u point
267      DO jj = 1, jpjm1
268         DO ji = 1, fs_jpim1   ! vector loop
269            umask_i(ji,jj)  = ssmask(ji,jj) * ssmask(ji+1,jj  )  * MIN(1._wp,SUM(umask(ji,jj,:)))
270            vmask_i(ji,jj)  = ssmask(ji,jj) * ssmask(ji  ,jj+1)  * MIN(1._wp,SUM(vmask(ji,jj,:)))
271         END DO
272         DO ji = 1, jpim1      ! NO vector opt.
273            fmask_i(ji,jj) =  ssmask(ji,jj  ) * ssmask(ji+1,jj  )   &
274               &            * ssmask(ji,jj+1) * ssmask(ji+1,jj+1) * MIN(1._wp,SUM(fmask(ji,jj,:)))
275         END DO
276      END DO
277      CALL lbc_lnk( umask, 'U', 1._wp )      ! Lateral boundary conditions
278      CALL lbc_lnk( vmask, 'V', 1._wp )
279      CALL lbc_lnk( fmask, 'F', 1._wp )
280      CALL lbc_lnk( umask_i, 'U', 1._wp )      ! Lateral boundary conditions
281      CALL lbc_lnk( vmask_i, 'V', 1._wp )
282      CALL lbc_lnk( fmask_i, 'F', 1._wp )
283
284      ! 3. Ocean/land mask at wu-, wv- and w points
285      !----------------------------------------------
286      wmask (:,:,1) = tmask(:,:,1) ! ????????
287      wumask(:,:,1) = umask(:,:,1) ! ????????
288      wvmask(:,:,1) = vmask(:,:,1) ! ????????
289      DO jk=2,jpk
290         wmask (:,:,jk)=tmask(:,:,jk) * tmask(:,:,jk-1)
291         wumask(:,:,jk)=umask(:,:,jk) * umask(:,:,jk-1)   
292         wvmask(:,:,jk)=vmask(:,:,jk) * vmask(:,:,jk-1)
293      END DO
294
295      ! 4. ocean/land mask for the elliptic equation
296      ! --------------------------------------------
297      bmask(:,:) = ssmask(:,:)       ! elliptic equation is written at t-point
298      !
299      !                               ! Boundary conditions
300      !                                    ! cyclic east-west : bmask must be set to 0. on rows 1 and jpi
301      IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN
302         bmask( 1 ,:) = 0._wp
303         bmask(jpi,:) = 0._wp
304      ENDIF
305      IF( nperio == 2 ) THEN               ! south symmetric :  bmask must be set to 0. on row 1
306         bmask(:, 1 ) = 0._wp
307      ENDIF
308      !                                    ! north fold :
309      IF( nperio == 3 .OR. nperio == 4 ) THEN   ! T-pt pivot : bmask set to 0. on row jpj and on half jpjglo-1 row
310         DO ji = 1, jpi                     
311            ii = ji + nimpp - 1
312            bmask(ji,jpj-1) = bmask(ji,jpj-1) * tpol(ii)
313            bmask(ji,jpj  ) = 0._wp
314         END DO
315      ENDIF
316      IF( nperio == 5 .OR. nperio == 6 ) THEN   ! F-pt pivot and T-pt elliptic eq. : bmask set to 0. on row jpj
317         bmask(:,jpj) = 0._wp
318      ENDIF
319      !
320      IF( lk_mpp ) THEN                    ! mpp specificities
321         !                                      ! bmask is set to zero on the overlap region
322         IF( nbondi /= -1 .AND. nbondi /= 2 )   bmask(  1 :jpreci,:) = 0._wp
323         IF( nbondi /=  1 .AND. nbondi /= 2 )   bmask(nlci:jpi   ,:) = 0._wp
324         IF( nbondj /= -1 .AND. nbondj /= 2 )   bmask(:,  1 :jprecj) = 0._wp
325         IF( nbondj /=  1 .AND. nbondj /= 2 )   bmask(:,nlcj:jpj   ) = 0._wp
326         !
327         IF( npolj == 3 .OR. npolj == 4 ) THEN  ! north fold : bmask must be set to 0. on rows jpj-1 and jpj
328            DO ji = 1, nlci
329               ii = ji + nimpp - 1
330               bmask(ji,nlcj-1) = bmask(ji,nlcj-1) * tpol(ii)
331               bmask(ji,nlcj  ) = 0._wp
332            END DO
333         ENDIF
334         IF( npolj == 5 .OR. npolj == 6 ) THEN  ! F-pt pivot and T-pt elliptic eq. : bmask set to 0. on row jpj
335            DO ji = 1, nlci
336               bmask(ji,nlcj  ) = 0._wp
337            END DO
338         ENDIF
339      ENDIF
340
341
342      ! mask for second order calculation of vorticity
343      ! ----------------------------------------------
344      CALL dom_msk_nsa
345
346     
347      ! Lateral boundary conditions on velocity (modify fmask)
348      ! ---------------------------------------     
349      DO jk = 1, jpk
350         zwf(:,:) = fmask(:,:,jk)         
351         DO jj = 2, jpjm1
352            DO ji = fs_2, fs_jpim1   ! vector opt.
353               IF( fmask(ji,jj,jk) == 0._wp ) THEN
354                  fmask(ji,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,jj), zwf(ji,jj+1),   &
355                     &                                           zwf(ji-1,jj), zwf(ji,jj-1)  )  )
356               ENDIF
357            END DO
358         END DO
359         DO jj = 2, jpjm1
360            IF( fmask(1,jj,jk) == 0._wp ) THEN
361               fmask(1  ,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1) ) )
362            ENDIF
363            IF( fmask(jpi,jj,jk) == 0._wp ) THEN
364               fmask(jpi,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) )
365            ENDIF
366         END DO         
367         DO ji = 2, jpim1
368            IF( fmask(ji,1,jk) == 0._wp ) THEN
369               fmask(ji, 1 ,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1) ) )
370            ENDIF
371            IF( fmask(ji,jpj,jk) == 0._wp ) THEN
372               fmask(ji,jpj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpjm1) ) )
373            ENDIF
374         END DO
375      END DO
376      !
377      IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN   ! ORCA_R2 configuration
378         !                                                 ! Increased lateral friction near of some straits
379         IF( nn_cla == 0 ) THEN
380            !                                ! Gibraltar strait  : partial slip (fmask=0.5)
381            ij0 = 101   ;   ij1 = 101
382            ii0 = 139   ;   ii1 = 140   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  0.5_wp
383            ij0 = 102   ;   ij1 = 102
384            ii0 = 139   ;   ii1 = 140   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  0.5_wp
385            !
386            !                                ! Bab el Mandeb : partial slip (fmask=1)
387            ij0 =  87   ;   ij1 =  88
388            ii0 = 160   ;   ii1 = 160   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  1._wp
389            ij0 =  88   ;   ij1 =  88
390            ii0 = 159   ;   ii1 = 159   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  1._wp
391            !
392         ENDIF
393         !                                ! Danish straits  : strong slip (fmask > 2)
394! We keep this as an example but it is instable in this case
395!         ij0 = 115   ;   ij1 = 115
396!         ii0 = 145   ;   ii1 = 146   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 4._wp
397!         ij0 = 116   ;   ij1 = 116
398!         ii0 = 145   ;   ii1 = 146   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 4._wp
399         !
400      ENDIF
401      !
402      IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN   ! ORCA R1 configuration
403         !                                                 ! Increased lateral friction near of some straits
404         ! This dirty section will be suppressed by simplification process:
405         ! all this will come back in input files
406         ! Currently these hard-wired indices relate to configuration with
407         ! extend grid (jpjglo=332)
408         !
409         isrow = 332 - jpjglo
410         !
411         IF(lwp) WRITE(numout,*)
412         IF(lwp) WRITE(numout,*) '   orca_r1: increase friction near the following straits : '
413         IF(lwp) WRITE(numout,*) '      Gibraltar '
414         ii0 = 282           ;   ii1 = 283        ! Gibraltar Strait
415         ij0 = 201 + isrow   ;   ij1 = 241 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp 
416
417         IF(lwp) WRITE(numout,*) '      Bhosporus '
418         ii0 = 314           ;   ii1 = 315        ! Bhosporus Strait
419         ij0 = 208 + isrow   ;   ij1 = 248 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp 
420
421         IF(lwp) WRITE(numout,*) '      Makassar (Top) '
422         ii0 =  48           ;   ii1 =  48        ! Makassar Strait (Top)
423         ij0 = 149 + isrow   ;   ij1 = 190 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp 
424
425         IF(lwp) WRITE(numout,*) '      Lombok '
426         ii0 =  44           ;   ii1 =  44        ! Lombok Strait
427         ij0 = 124 + isrow   ;   ij1 = 165 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp 
428
429         IF(lwp) WRITE(numout,*) '      Ombai '
430         ii0 =  53           ;   ii1 =  53        ! Ombai Strait
431         ij0 = 124 + isrow   ;   ij1 = 165 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp 
432
433         IF(lwp) WRITE(numout,*) '      Timor Passage '
434         ii0 =  56           ;   ii1 =  56        ! Timor Passage
435         ij0 = 124 + isrow   ;   ij1 = 165 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp 
436
437         IF(lwp) WRITE(numout,*) '      West Halmahera '
438         ii0 =  58           ;   ii1 =  58        ! West Halmahera Strait
439         ij0 = 141 + isrow   ;   ij1 = 182 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp 
440
441         IF(lwp) WRITE(numout,*) '      East Halmahera '
442         ii0 =  55           ;   ii1 =  55        ! East Halmahera Strait
443         ij0 = 141 + isrow   ;   ij1 = 182 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp 
444         !
445      ENDIF
446      !
447      CALL lbc_lnk( fmask, 'F', 1._wp )      ! Lateral boundary conditions on fmask
448
449      ! CAUTION : The fmask may be further modified in dyn_vor_init ( dynvor.F90 )
450           
451      IF( nprint == 1 .AND. lwp ) THEN      ! Control print
452         imsk(:,:) = INT( tmask_i(:,:) )
453         WRITE(numout,*) ' tmask_i : '
454         CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1,   &
455               &                           1, jpj, 1, 1, numout)
456         WRITE (numout,*)
457         WRITE (numout,*) ' dommsk: tmask for each level'
458         WRITE (numout,*) ' ----------------------------'
459         DO jk = 1, jpk
460            imsk(:,:) = INT( tmask(:,:,jk) )
461
462            WRITE(numout,*)
463            WRITE(numout,*) ' level = ',jk
464            CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1,   &
465               &                              1, jpj, 1, 1, numout)
466         END DO
467         WRITE(numout,*)
468         WRITE(numout,*) ' dom_msk: vmask for each level'
469         WRITE(numout,*) ' -----------------------------'
470         DO jk = 1, jpk
471            imsk(:,:) = INT( vmask(:,:,jk) )
472            WRITE(numout,*)
473            WRITE(numout,*) ' level = ',jk
474            CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1,   &
475               &                              1, jpj, 1, 1, numout)
476         END DO
477         WRITE(numout,*)
478         WRITE(numout,*) ' dom_msk: fmask for each level'
479         WRITE(numout,*) ' -----------------------------'
480         DO jk = 1, jpk
481            imsk(:,:) = INT( fmask(:,:,jk) )
482            WRITE(numout,*)
483            WRITE(numout,*) ' level = ',jk
484            CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1,   &
485               &                              1, jpj, 1, 1, numout )
486         END DO
487         WRITE(numout,*)
488         WRITE(numout,*) ' dom_msk: bmask '
489         WRITE(numout,*) ' ---------------'
490         WRITE(numout,*)
491         imsk(:,:) = INT( bmask(:,:) )
492         CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1,   &
493            &                              1, jpj, 1, 1, numout )
494      ENDIF
495      !
496      CALL wrk_dealloc( jpi, jpj, imsk )
497      CALL wrk_dealloc( jpi, jpj, zwf  )
498      !
499      IF( nn_timing == 1 )  CALL timing_stop('dom_msk')
500      !
501   END SUBROUTINE dom_msk
502
503#if defined key_noslip_accurate
504   !!----------------------------------------------------------------------
505   !!   'key_noslip_accurate' :         accurate no-slip boundary condition
506   !!----------------------------------------------------------------------
507   
508   SUBROUTINE dom_msk_nsa
509      !!---------------------------------------------------------------------
510      !!                 ***  ROUTINE dom_msk_nsa  ***
511      !!
512      !! ** Purpose :
513      !!
514      !! ** Method  :
515      !!
516      !! ** Action :
517      !!----------------------------------------------------------------------
518      INTEGER  ::   ji, jj, jk, jl      ! dummy loop indices
519      INTEGER  ::   ine, inw, ins, inn, itest, ierror, iind, ijnd
520      REAL(wp) ::   zaa
521      !!---------------------------------------------------------------------
522      !
523      IF( nn_timing == 1 )  CALL timing_start('dom_msk_nsa')
524      !
525      IF(lwp) WRITE(numout,*)
526      IF(lwp) WRITE(numout,*) 'dom_msk_nsa : noslip accurate boundary condition'
527      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   using Schchepetkin and O Brian scheme'
528      IF( lk_mpp )   CALL ctl_stop( ' mpp version is not yet implemented' )
529
530      ! mask for second order calculation of vorticity
531      ! ----------------------------------------------
532      ! noslip boundary condition: fmask=1  at convex corner, store
533      ! index of straight coast meshes ( 'west', refering to a coast,
534      ! means west of the ocean, aso)
535     
536      DO jk = 1, jpk
537         DO jl = 1, 4
538            npcoa(jl,jk) = 0
539            DO ji = 1, 2*(jpi+jpj)
540               nicoa(ji,jl,jk) = 0
541               njcoa(ji,jl,jk) = 0
542            END DO
543         END DO
544      END DO
545     
546      IF( jperio == 2 ) THEN
547         WRITE(numout,*) ' '
548         WRITE(numout,*) ' symetric boundary conditions need special'
549         WRITE(numout,*) ' treatment not implemented. we stop.'
550         STOP
551      ENDIF
552     
553      ! convex corners
554     
555      DO jk = 1, jpkm1
556         DO jj = 1, jpjm1
557            DO ji = 1, jpim1
558               zaa = tmask(ji  ,jj,jk) + tmask(ji  ,jj+1,jk)   &
559                  &+ tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk)
560               IF( ABS(zaa-3._wp) <= 0.1_wp )   fmask(ji,jj,jk) = 1._wp
561            END DO
562         END DO
563      END DO
564
565      ! north-south straight coast
566
567      DO jk = 1, jpkm1
568         inw = 0
569         ine = 0
570         DO jj = 2, jpjm1
571            DO ji = 2, jpim1
572               zaa = tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk)
573               IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN
574                  inw = inw + 1
575                  nicoa(inw,1,jk) = ji
576                  njcoa(inw,1,jk) = jj
577                  IF( nprint == 1 ) WRITE(numout,*) ' west  : ', jk, inw, ji, jj
578               ENDIF
579               zaa = tmask(ji,jj,jk) + tmask(ji,jj+1,jk)
580               IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN
581                  ine = ine + 1
582                  nicoa(ine,2,jk) = ji
583                  njcoa(ine,2,jk) = jj
584                  IF( nprint == 1 ) WRITE(numout,*) ' east  : ', jk, ine, ji, jj
585               ENDIF
586            END DO
587         END DO
588         npcoa(1,jk) = inw
589         npcoa(2,jk) = ine
590      END DO
591
592      ! west-east straight coast
593
594      DO jk = 1, jpkm1
595         ins = 0
596         inn = 0
597         DO jj = 2, jpjm1
598            DO ji =2, jpim1
599               zaa = tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk)
600               IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN
601                  ins = ins + 1
602                  nicoa(ins,3,jk) = ji
603                  njcoa(ins,3,jk) = jj
604                  IF( nprint == 1 ) WRITE(numout,*) ' south : ', jk, ins, ji, jj
605               ENDIF
606               zaa = tmask(ji+1,jj,jk) + tmask(ji,jj,jk)
607               IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN
608                  inn = inn + 1
609                  nicoa(inn,4,jk) = ji
610                  njcoa(inn,4,jk) = jj
611                  IF( nprint == 1 ) WRITE(numout,*) ' north : ', jk, inn, ji, jj
612               ENDIF
613            END DO
614         END DO
615         npcoa(3,jk) = ins
616         npcoa(4,jk) = inn
617      END DO
618
619      itest = 2 * ( jpi + jpj )
620      DO jk = 1, jpk
621         IF( npcoa(1,jk) > itest .OR. npcoa(2,jk) > itest .OR.   &
622             npcoa(3,jk) > itest .OR. npcoa(4,jk) > itest ) THEN
623           
624            WRITE(ctmp1,*) ' level jk = ',jk
625            WRITE(ctmp2,*) ' straight coast index arraies are too small.:'
626            WRITE(ctmp3,*) ' npe, npw, nps, npn = ', npcoa(1,jk), npcoa(2,jk),   &
627                &                                     npcoa(3,jk), npcoa(4,jk)
628            WRITE(ctmp4,*) ' 2*(jpi+jpj) = ',itest,'. we stop.'
629            CALL ctl_stop( ctmp1, ctmp2, ctmp3, ctmp4 )
630        ENDIF
631      END DO
632
633      ierror = 0
634      iind = 0
635      ijnd = 0
636      IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 )   iind = 2
637      IF( nperio == 3 .OR. nperio == 4 .OR. nperio == 5 .OR. nperio == 6 )   ijnd = 2
638      DO jk = 1, jpk
639         DO jl = 1, npcoa(1,jk)
640            IF( nicoa(jl,1,jk)+3 > jpi+iind ) THEN
641               ierror = ierror+1
642               icoord(ierror,1) = nicoa(jl,1,jk)
643               icoord(ierror,2) = njcoa(jl,1,jk)
644               icoord(ierror,3) = jk
645            ENDIF
646         END DO
647         DO jl = 1, npcoa(2,jk)
648            IF(nicoa(jl,2,jk)-2 < 1-iind ) THEN
649               ierror = ierror + 1
650               icoord(ierror,1) = nicoa(jl,2,jk)
651               icoord(ierror,2) = njcoa(jl,2,jk)
652               icoord(ierror,3) = jk
653            ENDIF
654         END DO
655         DO jl = 1, npcoa(3,jk)
656            IF( njcoa(jl,3,jk)+3 > jpj+ijnd ) THEN
657               ierror = ierror + 1
658               icoord(ierror,1) = nicoa(jl,3,jk)
659               icoord(ierror,2) = njcoa(jl,3,jk)
660               icoord(ierror,3) = jk
661            ENDIF
662         END DO
663         DO jl = 1, npcoa(4,jk)
664            IF( njcoa(jl,4,jk)-2 < 1) THEN
665               ierror=ierror + 1
666               icoord(ierror,1) = nicoa(jl,4,jk)
667               icoord(ierror,2) = njcoa(jl,4,jk)
668               icoord(ierror,3) = jk
669            ENDIF
670         END DO
671      END DO
672     
673      IF( ierror > 0 ) THEN
674         IF(lwp) WRITE(numout,*)
675         IF(lwp) WRITE(numout,*) '              Problem on lateral conditions'
676         IF(lwp) WRITE(numout,*) '                 Bad marking off at points:'
677         DO jl = 1, ierror
678            IF(lwp) WRITE(numout,*) 'Level:',icoord(jl,3),   &
679               &                  '  Point(',icoord(jl,1),',',icoord(jl,2),')'
680         END DO
681         CALL ctl_stop( 'We stop...' )
682      ENDIF
683      !
684      IF( nn_timing == 1 )  CALL timing_stop('dom_msk_nsa')
685      !
686   END SUBROUTINE dom_msk_nsa
687
688#else
689   !!----------------------------------------------------------------------
690   !!   Default option :                                      Empty routine
691   !!----------------------------------------------------------------------
692   SUBROUTINE dom_msk_nsa       
693   END SUBROUTINE dom_msk_nsa
694#endif
695   
696   !!======================================================================
697END MODULE dommsk
Note: See TracBrowser for help on using the repository browser.