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/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

source: branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90 @ 5802

Last change on this file since 5802 was 5802, checked in by mathiot, 8 years ago

ice sheet coupling: reproducibility fixed in MISOMIP configuration + contrain bathy to be constant during all the run

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