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

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, 4 years ago

The Dr Hook changes from my perl code.

File size: 34.8 KB
RevLine 
[10047]1
[3]2MODULE dommsk
[1566]3   !!======================================================================
[3]4   !!                       ***  MODULE dommsk   ***
5   !! Ocean initialization : domain land/sea mask
[1566]6   !!======================================================================
7   !! History :  OPA  ! 1987-07  (G. Madec)  Original code
[2528]8   !!            6.0  ! 1993-03  (M. Guyon)  symetrical conditions (M. Guyon)
9   !!            7.0  ! 1996-01  (G. Madec)  suppression of common work arrays
[1566]10   !!             -   ! 1996-05  (G. Madec)  mask computed from tmask and sup-
11   !!                 !                      pression of the double computation of bmask
[2528]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
[1566]14   !!             -   ! 1998-05  (G. Roullet)  free surface
[2528]15   !!            8.2  ! 2000-03  (G. Madec)  no slip accurate
[1566]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   !!----------------------------------------------------------------------
[3]21
22   !!----------------------------------------------------------------------
23   !!   dom_msk        : compute land/ocean mask
[1566]24   !!   dom_msk_nsa    : update land/ocean mask when no-slip accurate option is used.
[3]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)
[32]30   USE lib_mpp
[367]31   USE dynspg_oce      ! choice/control of key cpp for surface pressure gradient
[3294]32   USE wrk_nemo        ! Memory allocation
[10047]33   USE domwri
[3294]34   USE timing          ! Timing
[3]35
[11738]36   USE yomhook, ONLY: lhook, dr_hook
37   USE parkind1, ONLY: jprb, jpim
38
[3]39   IMPLICIT NONE
40   PRIVATE
41
[2715]42   PUBLIC   dom_msk         ! routine called by inidom.F90
43   PUBLIC   dom_msk_alloc   ! routine called by nemogcm.F90
[3]44
[1601]45   !                            !!* Namelist namlbc : lateral boundary condition *
[4147]46   REAL(wp)        :: rn_shlat   ! type of lateral boundary condition on velocity
47   LOGICAL, PUBLIC :: ln_vorlat  !  consistency of vorticity boundary condition
[3294]48   !                                            with analytical eqs.
[2715]49
[3294]50
[2715]51   INTEGER, ALLOCATABLE, SAVE, DIMENSION(:,:) ::  icoord ! Workspace for dom_msk_nsa()
52
[3]53   !! * Substitutions
54#  include "vectopt_loop_substitute.h90"
[1566]55   !!----------------------------------------------------------------------
56   !! NEMO/OPA 3.2 , LODYC-IPSL  (2009)
[6486]57   !! $Id$
[2528]58   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
[1566]59   !!----------------------------------------------------------------------
[3]60CONTAINS
61   
[2715]62   INTEGER FUNCTION dom_msk_alloc()
[11738]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
[2715]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      !
[11738]80   IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
[2715]81   END FUNCTION dom_msk_alloc
82
83
[3]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.
[1528]94      !!      mbathy equals 0 over continental T-point
95      !!      and the number of ocean level over the ocean.
[3]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)
[2528]107      !!                  and mbathy(ji+1,jj) and mbathy(ji+1,jj+1) >= jk.
[3]108      !!      b-point : the same definition as for f-point of the first ocean
109      !!                level (surface level) but with 0 along coastlines.
[2528]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.
[3]113      !!
114      !!        The lateral friction is set through the value of fmask along
[1601]115      !!      the coast and topography. This value is defined by rn_shlat, a
[3]116      !!      namelist parameter:
[1601]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
[3]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      !!
[4328]131      !!      In case of open boundaries (lk_bdy=T):
[3]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      !!
[1566]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.)
[1601]140      !!                          =rn_shlat along lateral boundaries
[1566]141      !!               bmask    : land/ocean mask at barotropic stream
142      !!                          function point (=0. or 1.) and set to 0 along lateral boundaries
[2528]143      !!               tmask_i  : interior ocean mask
[3]144      !!----------------------------------------------------------------------
[2715]145      !
[454]146      INTEGER  ::   ji, jj, jk      ! dummy loop indices
[2715]147      INTEGER  ::   iif, iil, ii0, ii1, ii   ! local integers
148      INTEGER  ::   ijf, ijl, ij0, ij1       !   -       -
[4147]149      INTEGER  ::   ios
[5385]150      INTEGER  ::   isrow                    ! index for ORCA1 starting row
[3294]151      INTEGER , POINTER, DIMENSION(:,:) ::  imsk
[6491]152      REAL(wp) ::  zphi_drake_passage, zshlat_antarc
[3294]153      REAL(wp), POINTER, DIMENSION(:,:) ::  zwf
[10047]154      REAL(wp) :: uvt(jpi,jpj)   ! dummy array for masking purposes.
[1601]155      !!
[3294]156      NAMELIST/namlbc/ rn_shlat, ln_vorlat
[11738]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
[3]165      !!---------------------------------------------------------------------
[3294]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      !
[4147]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 )
[4624]179      IF(lwm) WRITE ( numond, namlbc )
[1566]180     
181      IF(lwp) THEN                  ! control print
[3]182         WRITE(numout,*)
183         WRITE(numout,*) 'dommsk : ocean mask '
184         WRITE(numout,*) '~~~~~~'
[1566]185         WRITE(numout,*) '   Namelist namlbc'
[3294]186         WRITE(numout,*) '      lateral momentum boundary cond.    rn_shlat  = ',rn_shlat
187         WRITE(numout,*) '      consistency with analytical form   ln_vorlat = ',ln_vorlat 
[3]188      ENDIF
189
[2528]190      IF     (      rn_shlat == 0.               ) THEN   ;   IF(lwp) WRITE(numout,*) '   ocean lateral  free-slip '
[1601]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 )
[3]197      ENDIF
198
199      ! 1. Ocean/land mask at t-point (computed from mbathy)
200      ! -----------------------------
[1566]201      ! N.B. tmask has already the right boundary conditions since mbathy is ok
202      !
[2528]203      tmask(:,:,:) = 0._wp
[3]204      DO jk = 1, jpk
205         DO jj = 1, jpj
206            DO ji = 1, jpi
[2528]207               IF( REAL( mbathy(ji,jj) - jk, wp ) + 0.1_wp >= 0._wp )   tmask(ji,jj,jk) = 1._wp
[3]208            END DO 
209         END DO 
210      END DO 
[4990]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 
[3]224
[1566]225!!gm  ????
[255]226#if defined key_zdfkpp
[1601]227      IF( cp_cfg == 'orca' ) THEN
228         IF( jp_cfg == 2 )   THEN       ! land point on Bab el Mandeb zonal section
[291]229            ij0 =  87   ;   ij1 =  88
230            ii0 = 160   ;   ii1 = 161
[2528]231            tmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 0._wp
[291]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
[255]242#endif
[1566]243!!gm end
[291]244
[3]245      ! Interior domain mask (used for global sum)
246      ! --------------------
[4990]247      tmask_i(:,:) = ssmask(:,:)            ! (ISH) tmask_i = 1 even on the ice shelf
[10047]248
[3]249      iif = jpreci                         ! ???
250      iil = nlci - jpreci + 1
251      ijf = jprecj                         ! ???
252      ijl = nlcj - jprecj + 1
253
[2528]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)
[3]258
259      ! north fold mask
[1566]260      ! ---------------
[2528]261      tpol(1:jpiglo) = 1._wp 
262      fpol(1:jpiglo) = 1._wp
[3]263      IF( jperio == 3 .OR. jperio == 4 ) THEN      ! T-point pivot
[2528]264         tpol(jpiglo/2+1:jpiglo) = 0._wp
265         fpol(     1    :jpiglo) = 0._wp
[1566]266         IF( mjg(nlej) == jpjglo ) THEN                  ! only half of the nlcj-1 row
[291]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
[3]271      ENDIF
[10047]272
273
[3]274      IF( jperio == 5 .OR. jperio == 6 ) THEN      ! F-point pivot
[2528]275         tpol(     1    :jpiglo) = 0._wp
276         fpol(jpiglo/2+1:jpiglo) = 0._wp
[3]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)
[1271]286            END DO
[1694]287            DO ji = 1, jpim1      ! NO vector opt.
[3]288               fmask(ji,jj,jk) = tmask(ji,jj  ,jk) * tmask(ji+1,jj  ,jk)   &
[62]289                  &            * tmask(ji,jj+1,jk) * tmask(ji+1,jj+1,jk)
[3]290            END DO
[10047]291         END DO
[3]292      END DO
[4990]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
[2528]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 )
[4990]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 )
[3]310
[10047]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
[5120]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
[3]356
357      ! 4. ocean/land mask for the elliptic equation
358      ! --------------------------------------------
[4990]359      bmask(:,:) = ssmask(:,:)       ! elliptic equation is written at t-point
[1566]360      !
361      !                               ! Boundary conditions
362      !                                    ! cyclic east-west : bmask must be set to 0. on rows 1 and jpi
[3]363      IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN
[2528]364         bmask( 1 ,:) = 0._wp
365         bmask(jpi,:) = 0._wp
[3]366      ENDIF
[1566]367      IF( nperio == 2 ) THEN               ! south symmetric :  bmask must be set to 0. on row 1
[2528]368         bmask(:, 1 ) = 0._wp
[3]369      ENDIF
[1566]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                     
[1528]373            ii = ji + nimpp - 1
374            bmask(ji,jpj-1) = bmask(ji,jpj-1) * tpol(ii)
[2528]375            bmask(ji,jpj  ) = 0._wp
[1528]376         END DO
[3]377      ENDIF
[1566]378      IF( nperio == 5 .OR. nperio == 6 ) THEN   ! F-pt pivot and T-pt elliptic eq. : bmask set to 0. on row jpj
[2528]379         bmask(:,jpj) = 0._wp
[3]380      ENDIF
[1566]381      !
382      IF( lk_mpp ) THEN                    ! mpp specificities
383         !                                      ! bmask is set to zero on the overlap region
[2528]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
[1566]388         !
389         IF( npolj == 3 .OR. npolj == 4 ) THEN  ! north fold : bmask must be set to 0. on rows jpj-1 and jpj
[1528]390            DO ji = 1, nlci
391               ii = ji + nimpp - 1
392               bmask(ji,nlcj-1) = bmask(ji,nlcj-1) * tpol(ii)
[2528]393               bmask(ji,nlcj  ) = 0._wp
[1528]394            END DO
[32]395         ENDIF
[1566]396         IF( npolj == 5 .OR. npolj == 6 ) THEN  ! F-pt pivot and T-pt elliptic eq. : bmask set to 0. on row jpj
[1528]397            DO ji = 1, nlci
[2528]398               bmask(ji,nlcj  ) = 0._wp
[1528]399            END DO
[32]400         ENDIF
[3]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)
[1566]410      ! ---------------------------------------     
[3]411      DO jk = 1, jpk
[1566]412         zwf(:,:) = fmask(:,:,jk)         
[3]413         DO jj = 2, jpjm1
414            DO ji = fs_2, fs_jpim1   ! vector opt.
[2528]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)  )  )
[3]418               ENDIF
419            END DO
420         END DO
421         DO jj = 2, jpjm1
[2528]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) ) )
[3]424            ENDIF
[2528]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) ) )
[3]427            ENDIF
[1566]428         END DO         
[3]429         DO ji = 2, jpim1
[2528]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) ) )
[3]432            ENDIF
[2528]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) ) )
[3]435            ENDIF
436         END DO
437      END DO
[1566]438      !
439      IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN   ! ORCA_R2 configuration
440         !                                                 ! Increased lateral friction near of some straits
[2528]441         IF( nn_cla == 0 ) THEN
[1273]442            !                                ! Gibraltar strait  : partial slip (fmask=0.5)
443            ij0 = 101   ;   ij1 = 101
[2528]444            ii0 = 139   ;   ii1 = 140   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  0.5_wp
[32]445            ij0 = 102   ;   ij1 = 102
[2528]446            ii0 = 139   ;   ii1 = 140   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  0.5_wp
[1273]447            !
448            !                                ! Bab el Mandeb : partial slip (fmask=1)
449            ij0 =  87   ;   ij1 =  88
[2528]450            ii0 = 160   ;   ii1 = 160   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  1._wp
[1273]451            ij0 =  88   ;   ij1 =  88
[2528]452            ii0 = 159   ;   ii1 = 159   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  1._wp
[1273]453            !
[3]454         ENDIF
[1707]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
[2528]458!         ii0 = 145   ;   ii1 = 146   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 4._wp
[1707]459!         ij0 = 116   ;   ij1 = 116
[2528]460!         ii0 = 145   ;   ii1 = 146   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 4._wp
[3]461         !
462      ENDIF
[1566]463      !
[2528]464      IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN   ! ORCA R1 configuration
465         !                                                 ! Increased lateral friction near of some straits
[5506]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         !
[2528]473         IF(lwp) WRITE(numout,*)
474         IF(lwp) WRITE(numout,*) '   orca_r1: increase friction near the following straits : '
475         IF(lwp) WRITE(numout,*) '      Gibraltar '
[5385]476         ii0 = 282           ;   ii1 = 283        ! Gibraltar Strait
[6487]477         ij0 = 241 - isrow   ;   ij1 = 241 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp 
[1566]478
[2528]479         IF(lwp) WRITE(numout,*) '      Bhosporus '
[5385]480         ii0 = 314           ;   ii1 = 315        ! Bhosporus Strait
[6487]481         ij0 = 248 - isrow   ;   ij1 = 248 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp 
[2528]482
483         IF(lwp) WRITE(numout,*) '      Makassar (Top) '
[5385]484         ii0 =  48           ;   ii1 =  48        ! Makassar Strait (Top)
[6487]485         ij0 = 189 - isrow   ;   ij1 = 190 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp 
[2528]486
487         IF(lwp) WRITE(numout,*) '      Lombok '
[5385]488         ii0 =  44           ;   ii1 =  44        ! Lombok Strait
[6487]489         ij0 = 164 - isrow   ;   ij1 = 165 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp 
[2528]490
491         IF(lwp) WRITE(numout,*) '      Ombai '
[5385]492         ii0 =  53           ;   ii1 =  53        ! Ombai Strait
[6487]493         ij0 = 164 - isrow   ;   ij1 = 165 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp 
[2528]494
495         IF(lwp) WRITE(numout,*) '      Timor Passage '
[5385]496         ii0 =  56           ;   ii1 =  56        ! Timor Passage
[6487]497         ij0 = 164 - isrow   ;   ij1 = 165 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp 
[2528]498
499         IF(lwp) WRITE(numout,*) '      West Halmahera '
[5385]500         ii0 =  58           ;   ii1 =  58        ! West Halmahera Strait
[6487]501         ij0 = 181 - isrow   ;   ij1 = 182 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp 
[2528]502
503         IF(lwp) WRITE(numout,*) '      East Halmahera '
[5385]504         ii0 =  55           ;   ii1 =  55        ! East Halmahera Strait
[6487]505         ij0 = 181 - isrow   ;   ij1 = 182 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp 
[2528]506         !
507      ENDIF
508      !
[6491]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      !
[2528]532      CALL lbc_lnk( fmask, 'F', 1._wp )      ! Lateral boundary conditions on fmask
533
[3294]534      ! CAUTION : The fmask may be further modified in dyn_vor_init ( dynvor.F90 )
[2528]535           
[1566]536      IF( nprint == 1 .AND. lwp ) THEN      ! Control print
[3]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,   &
[2528]578            &                              1, jpj, 1, 1, numout )
[3]579      ENDIF
[1566]580      !
[3294]581      CALL wrk_dealloc( jpi, jpj, imsk )
582      CALL wrk_dealloc( jpi, jpj, zwf  )
[2715]583      !
[3294]584      IF( nn_timing == 1 )  CALL timing_stop('dom_msk')
585      !
[11738]586      IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
[3]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      !!----------------------------------------------------------------------
[2715]604      INTEGER  ::   ji, jj, jk, jl      ! dummy loop indices
[1566]605      INTEGER  ::   ine, inw, ins, inn, itest, ierror, iind, ijnd
606      REAL(wp) ::   zaa
[11738]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
[3]615      !!---------------------------------------------------------------------
[3294]616      !
617      IF( nn_timing == 1 )  CALL timing_start('dom_msk_nsa')
618      !
[2715]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'
[8280]622      IF( lk_mpp )   CALL ctl_stop('STOP', ' mpp version is not yet implemented' )
[3]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.'
[8280]644         CALL ctl_stop('STOP', 'NEMO abort from dom_msk_nsa')
[3]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)   &
[32]653                  &+ tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk)
[2528]654               IF( ABS(zaa-3._wp) <= 0.1_wp )   fmask(ji,jj,jk) = 1._wp
[3]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)
[2528]667               IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN
[3]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)
[2528]674               IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN
[3]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)
[2528]694               IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN
[3]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)
[2528]701               IF( ABS(zaa-2._wp) <= 0.1_wp .AND. fmask(ji,jj,jk) == 0._wp ) THEN
[3]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
[474]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),   &
[32]721                &                                     npcoa(3,jk), npcoa(4,jk)
[474]722            WRITE(ctmp4,*) ' 2*(jpi+jpj) = ',itest,'. we stop.'
723            CALL ctl_stop( ctmp1, ctmp2, ctmp3, ctmp4 )
[3]724        ENDIF
725      END DO
726
727      ierror = 0
728      iind = 0
729      ijnd = 0
[2528]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
[3]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
[2528]757         DO jl = 1, npcoa(4,jk)
[3]758            IF( njcoa(jl,4,jk)-2 < 1) THEN
[2528]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
[3]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),   &
[32]773               &                  '  Point(',icoord(jl,1),',',icoord(jl,2),')'
[3]774         END DO
[474]775         CALL ctl_stop( 'We stop...' )
[3]776      ENDIF
[2715]777      !
[3294]778      IF( nn_timing == 1 )  CALL timing_stop('dom_msk_nsa')
779      !
[11738]780      IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
[3]781   END SUBROUTINE dom_msk_nsa
782
783#else
784   !!----------------------------------------------------------------------
785   !!   Default option :                                      Empty routine
786   !!----------------------------------------------------------------------
787   SUBROUTINE dom_msk_nsa       
[11738]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)
[3]797   END SUBROUTINE dom_msk_nsa
798#endif
799   
800   !!======================================================================
801END MODULE dommsk
Note: See TracBrowser for help on using the repository browser.