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

source: branches/UKMO/dev_r5021_nn_etau_revision/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90 @ 5446

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

Update UKMO/dev_r5021_nn_etau_revision branch to revision 5442 of the trunk.
File size: 31.2 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: all this will come back in input files
405         ! Currently these hard-wired indices relate to the original (pre-v3.6) configuration
406         ! which had a grid-size of 362x292.
407         ! This grid has been extended southwards for use with the under ice-shelf options (isf) introduced in v3.6.
408         ! The original domain can still be used optionally if the isf code is not activated.
409         ! An adjustment (isrow) is made to the hard-wired indices if the extended domain (362x332) is being used.
410         !
411         IF    ( jpjglo == 292 ) THEN  ;  isrow = 0  ! Using pre-v3.6 files or adjusted start row from isf-extended grid
412         ELSEIF( jpjglo == 332 ) THEN  ;  isrow = 40 ! Using full isf­extended domain.
413         ENDIF     
414
415         IF(lwp) WRITE(numout,*)
416         IF(lwp) WRITE(numout,*) '   orca_r1: increase friction near the following straits : '
417         IF(lwp) WRITE(numout,*) '      Gibraltar '
418         ii0 = 282           ;   ii1 = 283        ! Gibraltar Strait
419         ij0 = 201 + isrow   ;   ij1 = 201 + isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp 
420
421         IF(lwp) WRITE(numout,*) '      Bhosporus '
422         ii0 = 314           ;   ii1 = 315        ! Bhosporus Strait
423         ij0 = 208 + isrow   ;   ij1 = 208 + isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp 
424
425         IF(lwp) WRITE(numout,*) '      Makassar (Top) '
426         ii0 =  48           ;   ii1 =  48        ! Makassar Strait (Top)
427         ij0 = 149 + isrow   ;   ij1 = 150 + isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp 
428
429         IF(lwp) WRITE(numout,*) '      Lombok '
430         ii0 =  44           ;   ii1 =  44        ! Lombok Strait
431         ij0 = 124 + isrow   ;   ij1 = 125 + isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp 
432
433         IF(lwp) WRITE(numout,*) '      Ombai '
434         ii0 =  53           ;   ii1 =  53        ! Ombai Strait
435         ij0 = 124 + isrow   ;   ij1 = 125 + isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp 
436
437         IF(lwp) WRITE(numout,*) '      Timor Passage '
438         ii0 =  56           ;   ii1 =  56        ! Timor Passage
439         ij0 = 124 + isrow   ;   ij1 = 125 + isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp 
440
441         IF(lwp) WRITE(numout,*) '      West Halmahera '
442         ii0 =  58           ;   ii1 =  58        ! West Halmahera Strait
443         ij0 = 141 + isrow   ;   ij1 = 142 + isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp 
444
445         IF(lwp) WRITE(numout,*) '      East Halmahera '
446         ii0 =  55           ;   ii1 =  55        ! East Halmahera Strait
447         ij0 = 141 + isrow   ;   ij1 = 142 + isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp 
448         !
449      ENDIF
450      !
451      CALL lbc_lnk( fmask, 'F', 1._wp )      ! Lateral boundary conditions on fmask
452
453      ! CAUTION : The fmask may be further modified in dyn_vor_init ( dynvor.F90 )
454           
455      IF( nprint == 1 .AND. lwp ) THEN      ! Control print
456         imsk(:,:) = INT( tmask_i(:,:) )
457         WRITE(numout,*) ' tmask_i : '
458         CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1,   &
459               &                           1, jpj, 1, 1, numout)
460         WRITE (numout,*)
461         WRITE (numout,*) ' dommsk: tmask for each level'
462         WRITE (numout,*) ' ----------------------------'
463         DO jk = 1, jpk
464            imsk(:,:) = INT( tmask(:,:,jk) )
465
466            WRITE(numout,*)
467            WRITE(numout,*) ' level = ',jk
468            CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1,   &
469               &                              1, jpj, 1, 1, numout)
470         END DO
471         WRITE(numout,*)
472         WRITE(numout,*) ' dom_msk: vmask for each level'
473         WRITE(numout,*) ' -----------------------------'
474         DO jk = 1, jpk
475            imsk(:,:) = INT( vmask(:,:,jk) )
476            WRITE(numout,*)
477            WRITE(numout,*) ' level = ',jk
478            CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1,   &
479               &                              1, jpj, 1, 1, numout)
480         END DO
481         WRITE(numout,*)
482         WRITE(numout,*) ' dom_msk: fmask for each level'
483         WRITE(numout,*) ' -----------------------------'
484         DO jk = 1, jpk
485            imsk(:,:) = INT( fmask(:,:,jk) )
486            WRITE(numout,*)
487            WRITE(numout,*) ' level = ',jk
488            CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1,   &
489               &                              1, jpj, 1, 1, numout )
490         END DO
491         WRITE(numout,*)
492         WRITE(numout,*) ' dom_msk: bmask '
493         WRITE(numout,*) ' ---------------'
494         WRITE(numout,*)
495         imsk(:,:) = INT( bmask(:,:) )
496         CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1,   &
497            &                              1, jpj, 1, 1, numout )
498      ENDIF
499      !
500      CALL wrk_dealloc( jpi, jpj, imsk )
501      CALL wrk_dealloc( jpi, jpj, zwf  )
502      !
503      IF( nn_timing == 1 )  CALL timing_stop('dom_msk')
504      !
505   END SUBROUTINE dom_msk
506
507#if defined key_noslip_accurate
508   !!----------------------------------------------------------------------
509   !!   'key_noslip_accurate' :         accurate no-slip boundary condition
510   !!----------------------------------------------------------------------
511   
512   SUBROUTINE dom_msk_nsa
513      !!---------------------------------------------------------------------
514      !!                 ***  ROUTINE dom_msk_nsa  ***
515      !!
516      !! ** Purpose :
517      !!
518      !! ** Method  :
519      !!
520      !! ** Action :
521      !!----------------------------------------------------------------------
522      INTEGER  ::   ji, jj, jk, jl      ! dummy loop indices
523      INTEGER  ::   ine, inw, ins, inn, itest, ierror, iind, ijnd
524      REAL(wp) ::   zaa
525      !!---------------------------------------------------------------------
526      !
527      IF( nn_timing == 1 )  CALL timing_start('dom_msk_nsa')
528      !
529      IF(lwp) WRITE(numout,*)
530      IF(lwp) WRITE(numout,*) 'dom_msk_nsa : noslip accurate boundary condition'
531      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   using Schchepetkin and O Brian scheme'
532      IF( lk_mpp )   CALL ctl_stop( ' mpp version is not yet implemented' )
533
534      ! mask for second order calculation of vorticity
535      ! ----------------------------------------------
536      ! noslip boundary condition: fmask=1  at convex corner, store
537      ! index of straight coast meshes ( 'west', refering to a coast,
538      ! means west of the ocean, aso)
539     
540      DO jk = 1, jpk
541         DO jl = 1, 4
542            npcoa(jl,jk) = 0
543            DO ji = 1, 2*(jpi+jpj)
544               nicoa(ji,jl,jk) = 0
545               njcoa(ji,jl,jk) = 0
546            END DO
547         END DO
548      END DO
549     
550      IF( jperio == 2 ) THEN
551         WRITE(numout,*) ' '
552         WRITE(numout,*) ' symetric boundary conditions need special'
553         WRITE(numout,*) ' treatment not implemented. we stop.'
554         STOP
555      ENDIF
556     
557      ! convex corners
558     
559      DO jk = 1, jpkm1
560         DO jj = 1, jpjm1
561            DO ji = 1, jpim1
562               zaa = tmask(ji  ,jj,jk) + tmask(ji  ,jj+1,jk)   &
563                  &+ tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk)
564               IF( ABS(zaa-3._wp) <= 0.1_wp )   fmask(ji,jj,jk) = 1._wp
565            END DO
566         END DO
567      END DO
568
569      ! north-south straight coast
570
571      DO jk = 1, jpkm1
572         inw = 0
573         ine = 0
574         DO jj = 2, jpjm1
575            DO ji = 2, jpim1
576               zaa = tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk)
577               IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN
578                  inw = inw + 1
579                  nicoa(inw,1,jk) = ji
580                  njcoa(inw,1,jk) = jj
581                  IF( nprint == 1 ) WRITE(numout,*) ' west  : ', jk, inw, ji, jj
582               ENDIF
583               zaa = tmask(ji,jj,jk) + tmask(ji,jj+1,jk)
584               IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN
585                  ine = ine + 1
586                  nicoa(ine,2,jk) = ji
587                  njcoa(ine,2,jk) = jj
588                  IF( nprint == 1 ) WRITE(numout,*) ' east  : ', jk, ine, ji, jj
589               ENDIF
590            END DO
591         END DO
592         npcoa(1,jk) = inw
593         npcoa(2,jk) = ine
594      END DO
595
596      ! west-east straight coast
597
598      DO jk = 1, jpkm1
599         ins = 0
600         inn = 0
601         DO jj = 2, jpjm1
602            DO ji =2, jpim1
603               zaa = tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk)
604               IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN
605                  ins = ins + 1
606                  nicoa(ins,3,jk) = ji
607                  njcoa(ins,3,jk) = jj
608                  IF( nprint == 1 ) WRITE(numout,*) ' south : ', jk, ins, ji, jj
609               ENDIF
610               zaa = tmask(ji+1,jj,jk) + tmask(ji,jj,jk)
611               IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN
612                  inn = inn + 1
613                  nicoa(inn,4,jk) = ji
614                  njcoa(inn,4,jk) = jj
615                  IF( nprint == 1 ) WRITE(numout,*) ' north : ', jk, inn, ji, jj
616               ENDIF
617            END DO
618         END DO
619         npcoa(3,jk) = ins
620         npcoa(4,jk) = inn
621      END DO
622
623      itest = 2 * ( jpi + jpj )
624      DO jk = 1, jpk
625         IF( npcoa(1,jk) > itest .OR. npcoa(2,jk) > itest .OR.   &
626             npcoa(3,jk) > itest .OR. npcoa(4,jk) > itest ) THEN
627           
628            WRITE(ctmp1,*) ' level jk = ',jk
629            WRITE(ctmp2,*) ' straight coast index arraies are too small.:'
630            WRITE(ctmp3,*) ' npe, npw, nps, npn = ', npcoa(1,jk), npcoa(2,jk),   &
631                &                                     npcoa(3,jk), npcoa(4,jk)
632            WRITE(ctmp4,*) ' 2*(jpi+jpj) = ',itest,'. we stop.'
633            CALL ctl_stop( ctmp1, ctmp2, ctmp3, ctmp4 )
634        ENDIF
635      END DO
636
637      ierror = 0
638      iind = 0
639      ijnd = 0
640      IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 )   iind = 2
641      IF( nperio == 3 .OR. nperio == 4 .OR. nperio == 5 .OR. nperio == 6 )   ijnd = 2
642      DO jk = 1, jpk
643         DO jl = 1, npcoa(1,jk)
644            IF( nicoa(jl,1,jk)+3 > jpi+iind ) THEN
645               ierror = ierror+1
646               icoord(ierror,1) = nicoa(jl,1,jk)
647               icoord(ierror,2) = njcoa(jl,1,jk)
648               icoord(ierror,3) = jk
649            ENDIF
650         END DO
651         DO jl = 1, npcoa(2,jk)
652            IF(nicoa(jl,2,jk)-2 < 1-iind ) THEN
653               ierror = ierror + 1
654               icoord(ierror,1) = nicoa(jl,2,jk)
655               icoord(ierror,2) = njcoa(jl,2,jk)
656               icoord(ierror,3) = jk
657            ENDIF
658         END DO
659         DO jl = 1, npcoa(3,jk)
660            IF( njcoa(jl,3,jk)+3 > jpj+ijnd ) THEN
661               ierror = ierror + 1
662               icoord(ierror,1) = nicoa(jl,3,jk)
663               icoord(ierror,2) = njcoa(jl,3,jk)
664               icoord(ierror,3) = jk
665            ENDIF
666         END DO
667         DO jl = 1, npcoa(4,jk)
668            IF( njcoa(jl,4,jk)-2 < 1) THEN
669               ierror=ierror + 1
670               icoord(ierror,1) = nicoa(jl,4,jk)
671               icoord(ierror,2) = njcoa(jl,4,jk)
672               icoord(ierror,3) = jk
673            ENDIF
674         END DO
675      END DO
676     
677      IF( ierror > 0 ) THEN
678         IF(lwp) WRITE(numout,*)
679         IF(lwp) WRITE(numout,*) '              Problem on lateral conditions'
680         IF(lwp) WRITE(numout,*) '                 Bad marking off at points:'
681         DO jl = 1, ierror
682            IF(lwp) WRITE(numout,*) 'Level:',icoord(jl,3),   &
683               &                  '  Point(',icoord(jl,1),',',icoord(jl,2),')'
684         END DO
685         CALL ctl_stop( 'We stop...' )
686      ENDIF
687      !
688      IF( nn_timing == 1 )  CALL timing_stop('dom_msk_nsa')
689      !
690   END SUBROUTINE dom_msk_nsa
691
692#else
693   !!----------------------------------------------------------------------
694   !!   Default option :                                      Empty routine
695   !!----------------------------------------------------------------------
696   SUBROUTINE dom_msk_nsa       
697   END SUBROUTINE dom_msk_nsa
698#endif
699   
700   !!======================================================================
701END MODULE dommsk
Note: See TracBrowser for help on using the repository browser.