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

source: branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90 @ 6667

Last change on this file since 6667 was 6667, checked in by gm, 8 years ago

#1692 - branch SIMPLIF_2_usrdef: reduced domain_cfg.nc file: GYRE OK using usrdef or reading file

  • Property svn:keywords set to Id
File size: 17.2 KB
RevLine 
[3]1MODULE dommsk
[1566]2   !!======================================================================
[3]3   !!                       ***  MODULE dommsk   ***
4   !! Ocean initialization : domain land/sea mask
[1566]5   !!======================================================================
6   !! History :  OPA  ! 1987-07  (G. Madec)  Original code
[2528]7   !!            6.0  ! 1993-03  (M. Guyon)  symetrical conditions (M. Guyon)
8   !!            7.0  ! 1996-01  (G. Madec)  suppression of common work arrays
[6140]9   !!             -   ! 1996-05  (G. Madec)  mask computed from tmask
[2528]10   !!            8.0  ! 1997-02  (G. Madec)  mesh information put in domhgr.F
[6667]11   !!            8.1  ! 1997-07  (G. Madec)  modification of kbat and fmask
[1566]12   !!             -   ! 1998-05  (G. Roullet)  free surface
[2528]13   !!            8.2  ! 2000-03  (G. Madec)  no slip accurate
[1566]14   !!             -   ! 2001-09  (J.-M. Molines)  Open boundaries
15   !!   NEMO     1.0  ! 2002-08  (G. Madec)  F90: Free form and module
16   !!             -   ! 2005-11  (V. Garnier) Surface pressure gradient organization
17   !!            3.2  ! 2009-07  (R. Benshila) Suppression of rigid-lid option
[5836]18   !!            3.6  ! 2015-05  (P. Mathiot) ISF: add wmask,wumask and wvmask
[1566]19   !!----------------------------------------------------------------------
[3]20
21   !!----------------------------------------------------------------------
22   !!   dom_msk        : compute land/ocean mask
23   !!----------------------------------------------------------------------
24   USE oce             ! ocean dynamics and tracers
25   USE dom_oce         ! ocean space and time domain
[6140]26   !
[3]27   USE in_out_manager  ! I/O manager
28   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
[6140]29   USE lib_mpp         !
[3294]30   USE wrk_nemo        ! Memory allocation
31   USE timing          ! Timing
[3]32
33   IMPLICIT NONE
34   PRIVATE
35
[6140]36   PUBLIC   dom_msk    ! routine called by inidom.F90
[3]37
[1601]38   !                            !!* Namelist namlbc : lateral boundary condition *
[4147]39   REAL(wp)        :: rn_shlat   ! type of lateral boundary condition on velocity
40   LOGICAL, PUBLIC :: ln_vorlat  !  consistency of vorticity boundary condition
[3294]41   !                                            with analytical eqs.
[2715]42
[3]43   !! * Substitutions
44#  include "vectopt_loop_substitute.h90"
[1566]45   !!----------------------------------------------------------------------
46   !! NEMO/OPA 3.2 , LODYC-IPSL  (2009)
[1152]47   !! $Id$
[2528]48   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
[1566]49   !!----------------------------------------------------------------------
[3]50CONTAINS
[2715]51
[6667]52   SUBROUTINE dom_msk( k_top, k_bot )
[3]53      !!---------------------------------------------------------------------
54      !!                 ***  ROUTINE dom_msk  ***
55      !!
56      !! ** Purpose :   Compute land/ocean mask arrays at tracer points, hori-
[6125]57      !!      zontal velocity points (u & v), vorticity points (f) points.
[3]58      !!
[6667]59      !! ** Method  :   The ocean/land mask  at t-point is deduced from ko_top
60      !!      and ko_bot, the indices of the fist and last ocean t-levels which
61      !!      are either defined in usrdef_zgr or read in zgr_read.
62      !!                The velocity masks (umask, vmask, wmask, wumask, wvmask)
63      !!      are deduced from a product of the two neighboring tmask.
64      !!                The vorticity mask (fmask) is deduced from tmask taking
65      !!      into account the choice of lateral boundary condition (rn_shlat) :
[1601]66      !!         rn_shlat = 0, free slip  (no shear along the coast)
67      !!         rn_shlat = 2, no slip  (specified zero velocity at the coast)
68      !!         0 < rn_shlat < 2, partial slip   | non-linear velocity profile
69      !!         2 < rn_shlat, strong slip        | in the lateral boundary layer
[3]70      !!
[6667]71      !!      tmask_i : interior ocean mask at t-point, i.e. excluding duplicated
72      !!                rows/lines due to cyclic or North Fold boundaries as well
73      !!                as MPP halos.
74      !!      tmask_h : halo mask at t-point, i.e. excluding duplicated rows/lines
75      !!                due to cyclic or North Fold boundaries as well
76      !!                as MPP halos.
[3]77      !!
[4328]78      !!      In case of open boundaries (lk_bdy=T):
[6667]79      !!        - tmask is set to 1 on the points to be computed by the open
[3]80      !!          boundaries routines.
81      !!
[6667]82      !! ** Action :   tmask, umask, vmask, wmask, wumask, wvmask : land/ocean mask
83      !!                         at t-, u-, v- w, wu-, and wv-points (=0. or 1.)
84      !!               fmask   : land/ocean mask at f-point (=0., or =1., or
85      !!                         =rn_shlat along lateral boundaries)
86      !!               tmask_i : interior ocean mask
87      !!               tmask_h : halo mask
88      !!               ssmask , ssumask, ssvmask, ssfmask : 2D ocean mask
[3]89      !!----------------------------------------------------------------------
[6667]90      INTEGER, DIMENSION(:,:), INTENT(in) ::   k_top, k_bot   ! first and last ocean level
91      !
[5836]92      INTEGER  ::   ji, jj, jk               ! dummy loop indices
[2715]93      INTEGER  ::   iif, iil, ii0, ii1, ii   ! local integers
94      INTEGER  ::   ijf, ijl, ij0, ij1       !   -       -
[6667]95      INTEGER  ::   iktop, ikbot             !   -       -
[4147]96      INTEGER  ::   ios
[5385]97      INTEGER  ::   isrow                    ! index for ORCA1 starting row
[3294]98      REAL(wp), POINTER, DIMENSION(:,:) ::  zwf
[1601]99      !!
[3294]100      NAMELIST/namlbc/ rn_shlat, ln_vorlat
[3]101      !!---------------------------------------------------------------------
[3294]102      !
103      IF( nn_timing == 1 )  CALL timing_start('dom_msk')
104      !
[4147]105      REWIND( numnam_ref )              ! Namelist namlbc in reference namelist : Lateral momentum boundary condition
106      READ  ( numnam_ref, namlbc, IOSTAT = ios, ERR = 901 )
107901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namlbc in reference namelist', lwp )
108
109      REWIND( numnam_cfg )              ! Namelist namlbc in configuration namelist : Lateral momentum boundary condition
110      READ  ( numnam_cfg, namlbc, IOSTAT = ios, ERR = 902 )
111902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namlbc in configuration namelist', lwp )
[4624]112      IF(lwm) WRITE ( numond, namlbc )
[1566]113     
114      IF(lwp) THEN                  ! control print
[3]115         WRITE(numout,*)
116         WRITE(numout,*) 'dommsk : ocean mask '
117         WRITE(numout,*) '~~~~~~'
[1566]118         WRITE(numout,*) '   Namelist namlbc'
[3294]119         WRITE(numout,*) '      lateral momentum boundary cond.    rn_shlat  = ',rn_shlat
120         WRITE(numout,*) '      consistency with analytical form   ln_vorlat = ',ln_vorlat 
[3]121      ENDIF
122
[2528]123      IF     (      rn_shlat == 0.               ) THEN   ;   IF(lwp) WRITE(numout,*) '   ocean lateral  free-slip '
[1601]124      ELSEIF (      rn_shlat == 2.               ) THEN   ;   IF(lwp) WRITE(numout,*) '   ocean lateral  no-slip '
125      ELSEIF ( 0. < rn_shlat .AND. rn_shlat < 2. ) THEN   ;   IF(lwp) WRITE(numout,*) '   ocean lateral  partial-slip '
126      ELSEIF ( 2. < rn_shlat                     ) THEN   ;   IF(lwp) WRITE(numout,*) '   ocean lateral  strong-slip '
127      ELSE
128         WRITE(ctmp1,*) ' rn_shlat is negative = ', rn_shlat
129         CALL ctl_stop( ctmp1 )
[3]130      ENDIF
131
[6667]132
133      !  Ocean/land mask at t-point  (computed from ko_top and ko_bot)
134      ! ----------------------------
[1566]135      !
[6667]136      DO jj = 1, jpj
137         DO ji = 1, jpi
138            iktop = k_top(ji,jj)
139            ikbot = k_bot(ji,jj)
140            tmask(ji,jj,      1:iktop-1) = 0._wp      ! mask the iceshelves
141            tmask(ji,jj,iktop  :ikbot  ) = 1._wp
142            tmask(ji,jj,ikbot+1:jpkglo ) = 0._wp      ! mask the ocean topography
[3]143         END DO 
144      END DO 
[6667]145
146      ! 2D ocean mask (=1 if at least one level of the water column is ocean, =0 otherwise)
147      WHERE( k_bot(:,:) > 0 )   ;   ssmask(:,:) = 1._wp
148      ELSEWHERE                 ;   ssmask(:,:) = 0._wp
149      END WHERE
[4990]150     
151     
[6667]152      ! Interior domain mask  (used for global sum)
[3]153      ! --------------------
[6667]154      !
155      iif = jpreci   ;   iil = nlci - jpreci + 1
156      ijf = jprecj   ;   ijl = nlcj - jprecj + 1
157      !
158      !                          ! halo mask : 0 on the halo and 1 elsewhere
159      tmask_h(:,:) = 1._wp                 
[6140]160      tmask_h( 1 :iif,   :   ) = 0._wp      ! first columns
161      tmask_h(iil:jpi,   :   ) = 0._wp      ! last  columns (including mpp extra columns)
162      tmask_h(   :   , 1 :ijf) = 0._wp      ! first rows
163      tmask_h(   :   ,ijl:jpj) = 0._wp      ! last  rows (including mpp extra rows)
[6667]164      !
165      !                          ! north fold mask
[2528]166      tpol(1:jpiglo) = 1._wp 
167      fpol(1:jpiglo) = 1._wp
[3]168      IF( jperio == 3 .OR. jperio == 4 ) THEN      ! T-point pivot
[2528]169         tpol(jpiglo/2+1:jpiglo) = 0._wp
170         fpol(     1    :jpiglo) = 0._wp
[6667]171         IF( mjg(nlej) == jpjglo ) THEN                  ! only half of the nlcj-1 row for tmask_h
[291]172            DO ji = iif+1, iil-1
[6140]173               tmask_h(ji,nlej-1) = tmask_h(ji,nlej-1) * tpol(mig(ji))
[291]174            END DO
175         ENDIF
[3]176      ENDIF
[6667]177      !
[3]178      IF( jperio == 5 .OR. jperio == 6 ) THEN      ! F-point pivot
[2528]179         tpol(     1    :jpiglo) = 0._wp
180         fpol(jpiglo/2+1:jpiglo) = 0._wp
[3]181      ENDIF
[6667]182      !
183      !                          ! interior mask : 2D ocean mask x halo mask
184      tmask_i(:,:) = ssmask(:,:) * tmask_h(:,:)
[3]185
[6667]186
187      !  Ocean/land mask at u-, v-, and z-points (computed from tmask)
188      ! ----------------------------------------
[3]189      DO jk = 1, jpk
190         DO jj = 1, jpjm1
191            DO ji = 1, fs_jpim1   ! vector loop
192               umask(ji,jj,jk) = tmask(ji,jj  ,jk) * tmask(ji+1,jj  ,jk)
193               vmask(ji,jj,jk) = tmask(ji,jj  ,jk) * tmask(ji  ,jj+1,jk)
[1271]194            END DO
[1694]195            DO ji = 1, jpim1      ! NO vector opt.
[3]196               fmask(ji,jj,jk) = tmask(ji,jj  ,jk) * tmask(ji+1,jj  ,jk)   &
[62]197                  &            * tmask(ji,jj+1,jk) * tmask(ji+1,jj+1,jk)
[3]198            END DO
199         END DO
200      END DO
[6140]201      ! (ISF) MIN(1,SUM(umask)) is here to check if you have effectively at least 1 wet cell at u point
[4990]202      DO jj = 1, jpjm1
203         DO ji = 1, fs_jpim1   ! vector loop
[6667]204!!gm  simpler :
205!            ssumask(ji,jj)  = MIN(  1._wp , SUM( umask(ji,jj,:) )  )
206!            ssvmask(ji,jj)  = MIN(  1._wp , SUM( vmask(ji,jj,:) )  )
207!!gm
208!!gm  faster :
209!         ssumask(ji,jj) = ssmask(ji,jj) * tmask(ji+1,jj  )
210!         ssvmask(ji,jj) = ssmask(ji,jj) * tmask(ji  ,jj+1)
211!!gm
[6140]212            ssumask(ji,jj)  = ssmask(ji,jj) * ssmask(ji+1,jj  )  * MIN(1._wp,SUM(umask(ji,jj,:)))
213            ssvmask(ji,jj)  = ssmask(ji,jj) * ssmask(ji  ,jj+1)  * MIN(1._wp,SUM(vmask(ji,jj,:)))
[6667]214!!end
[4990]215         END DO
216         DO ji = 1, jpim1      ! NO vector opt.
[6667]217!!gm faster
218!            ssfmask(ji,jj) =  ssmask(ji,jj  ) * ssmask(ji+1,jj  )   &
219!               &            * ssmask(ji,jj+1) * ssmask(ji+1,jj+1)
220!!gm
[6140]221            ssfmask(ji,jj) =  ssmask(ji,jj  ) * ssmask(ji+1,jj  )   &
[4990]222               &            * ssmask(ji,jj+1) * ssmask(ji+1,jj+1) * MIN(1._wp,SUM(fmask(ji,jj,:)))
[6667]223!!gm
[4990]224         END DO
225      END DO
[6140]226      CALL lbc_lnk( umask  , 'U', 1._wp )      ! Lateral boundary conditions
227      CALL lbc_lnk( vmask  , 'V', 1._wp )
[6667]228!      CALL lbc_lnk( fmask  , 'F', 1._wp )         ! applied after the specification of lateral b.c.
229      CALL lbc_lnk( ssumask, 'U', 1._wp )
[6140]230      CALL lbc_lnk( ssvmask, 'V', 1._wp )
231      CALL lbc_lnk( ssfmask, 'F', 1._wp )
[3]232
[6667]233
234      ! Ocean/land mask at wu-, wv- and w points
[5120]235      !----------------------------------------------
[5836]236      wmask (:,:,1) = tmask(:,:,1)     ! surface
237      wumask(:,:,1) = umask(:,:,1)
238      wvmask(:,:,1) = vmask(:,:,1)
239      DO jk = 2, jpk                   ! interior values
240         wmask (:,:,jk) = tmask(:,:,jk) * tmask(:,:,jk-1)
241         wumask(:,:,jk) = umask(:,:,jk) * umask(:,:,jk-1)   
242         wvmask(:,:,jk) = vmask(:,:,jk) * vmask(:,:,jk-1)
[5120]243      END DO
[3]244
[6667]245
[3]246      ! Lateral boundary conditions on velocity (modify fmask)
[1566]247      ! ---------------------------------------     
[6667]248      CALL wrk_alloc( jpi,jpj,   zwf  )
249      !
[3]250      DO jk = 1, jpk
[1566]251         zwf(:,:) = fmask(:,:,jk)         
[3]252         DO jj = 2, jpjm1
253            DO ji = fs_2, fs_jpim1   ! vector opt.
[2528]254               IF( fmask(ji,jj,jk) == 0._wp ) THEN
255                  fmask(ji,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,jj), zwf(ji,jj+1),   &
256                     &                                           zwf(ji-1,jj), zwf(ji,jj-1)  )  )
[3]257               ENDIF
258            END DO
259         END DO
260         DO jj = 2, jpjm1
[2528]261            IF( fmask(1,jj,jk) == 0._wp ) THEN
262               fmask(1  ,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1) ) )
[3]263            ENDIF
[2528]264            IF( fmask(jpi,jj,jk) == 0._wp ) THEN
265               fmask(jpi,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) )
[3]266            ENDIF
[1566]267         END DO         
[3]268         DO ji = 2, jpim1
[2528]269            IF( fmask(ji,1,jk) == 0._wp ) THEN
270               fmask(ji, 1 ,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1) ) )
[3]271            ENDIF
[2528]272            IF( fmask(ji,jpj,jk) == 0._wp ) THEN
273               fmask(ji,jpj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpjm1) ) )
[3]274            ENDIF
275         END DO
276      END DO
[1566]277      !
[6667]278      CALL wrk_dealloc( jpi,jpj,   zwf  )
279      !
[1566]280      IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN   ! ORCA_R2 configuration
281         !                                                 ! Increased lateral friction near of some straits
[5836]282         !                                ! Gibraltar strait  : partial slip (fmask=0.5)
283         ij0 = 101   ;   ij1 = 101
284         ii0 = 139   ;   ii1 = 140   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  0.5_wp
285         ij0 = 102   ;   ij1 = 102
286         ii0 = 139   ;   ii1 = 140   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  0.5_wp
287         !
288         !                                ! Bab el Mandeb : partial slip (fmask=1)
289         ij0 =  87   ;   ij1 =  88
290         ii0 = 160   ;   ii1 = 160   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  1._wp
291         ij0 =  88   ;   ij1 =  88
292         ii0 = 159   ;   ii1 = 159   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  1._wp
293         !
[1707]294         !                                ! Danish straits  : strong slip (fmask > 2)
295! We keep this as an example but it is instable in this case
296!         ij0 = 115   ;   ij1 = 115
[2528]297!         ii0 = 145   ;   ii1 = 146   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 4._wp
[1707]298!         ij0 = 116   ;   ij1 = 116
[2528]299!         ii0 = 145   ;   ii1 = 146   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 4._wp
[3]300         !
301      ENDIF
[1566]302      !
[2528]303      IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN   ! ORCA R1 configuration
304         !                                                 ! Increased lateral friction near of some straits
[5506]305         ! This dirty section will be suppressed by simplification process:
306         ! all this will come back in input files
307         ! Currently these hard-wired indices relate to configuration with
308         ! extend grid (jpjglo=332)
309         !
310         isrow = 332 - jpjglo
311         !
[2528]312         IF(lwp) WRITE(numout,*)
313         IF(lwp) WRITE(numout,*) '   orca_r1: increase friction near the following straits : '
314         IF(lwp) WRITE(numout,*) '      Gibraltar '
[5385]315         ii0 = 282           ;   ii1 = 283        ! Gibraltar Strait
[5552]316         ij0 = 241 - isrow   ;   ij1 = 241 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp 
[1566]317
[2528]318         IF(lwp) WRITE(numout,*) '      Bhosporus '
[5385]319         ii0 = 314           ;   ii1 = 315        ! Bhosporus Strait
[5552]320         ij0 = 248 - isrow   ;   ij1 = 248 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp 
[2528]321
322         IF(lwp) WRITE(numout,*) '      Makassar (Top) '
[5385]323         ii0 =  48           ;   ii1 =  48        ! Makassar Strait (Top)
[5552]324         ij0 = 189 - isrow   ;   ij1 = 190 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp 
[2528]325
326         IF(lwp) WRITE(numout,*) '      Lombok '
[5385]327         ii0 =  44           ;   ii1 =  44        ! Lombok Strait
[5552]328         ij0 = 164 - isrow   ;   ij1 = 165 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp 
[2528]329
330         IF(lwp) WRITE(numout,*) '      Ombai '
[5385]331         ii0 =  53           ;   ii1 =  53        ! Ombai Strait
[5552]332         ij0 = 164 - isrow   ;   ij1 = 165 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp 
[2528]333
334         IF(lwp) WRITE(numout,*) '      Timor Passage '
[5385]335         ii0 =  56           ;   ii1 =  56        ! Timor Passage
[5552]336         ij0 = 164 - isrow   ;   ij1 = 165 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp 
[2528]337
338         IF(lwp) WRITE(numout,*) '      West Halmahera '
[5385]339         ii0 =  58           ;   ii1 =  58        ! West Halmahera Strait
[5552]340         ij0 = 181 - isrow   ;   ij1 = 182 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp 
[2528]341
342         IF(lwp) WRITE(numout,*) '      East Halmahera '
[5385]343         ii0 =  55           ;   ii1 =  55        ! East Halmahera Strait
[5552]344         ij0 = 181 - isrow   ;   ij1 = 182 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp 
[2528]345         !
346      ENDIF
347      !
348      CALL lbc_lnk( fmask, 'F', 1._wp )      ! Lateral boundary conditions on fmask
[6140]349      !
[3294]350      ! CAUTION : The fmask may be further modified in dyn_vor_init ( dynvor.F90 )
[1566]351      !
[3294]352      IF( nn_timing == 1 )  CALL timing_stop('dom_msk')
353      !
[3]354   END SUBROUTINE dom_msk
355   
356   !!======================================================================
357END MODULE dommsk
Note: See TracBrowser for help on using the repository browser.