source: branches/UKMO/dev_r5518_GO6_under_ice_relax_dr_hook/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90 @ 11738

Last change on this file since 11738 was 11738, checked in by marc, 12 months ago

The Dr Hook changes from my perl code.

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