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 NEMO/branches/UKMO/dev_r10037_shlat2d/src/OCE/DOM – NEMO

source: NEMO/branches/UKMO/dev_r10037_shlat2d/src/OCE/DOM/dommsk.F90 @ 10642

Last change on this file since 10642 was 10642, checked in by davestorkey, 5 years ago

UKMO dev_r10037_shlat2d branch: commit code changes.

File size: 16.0 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
[7646]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
[7646]19   !!            4.0  ! 2016-06  (G. Madec, S. Flavoni)  domain configuration / user defined interface
[1566]20   !!----------------------------------------------------------------------
[3]21
22   !!----------------------------------------------------------------------
[7646]23   !!   dom_msk       : compute land/ocean mask
[3]24   !!----------------------------------------------------------------------
[7646]25   USE oce            ! ocean dynamics and tracers
26   USE dom_oce        ! ocean space and time domain
27   USE usrdef_fmask   ! user defined fmask
[9124]28   USE bdy_oce        ! open boundary
29   !
[7646]30   USE in_out_manager ! I/O manager
[9600]31   USE iom            ! IOM library
[7646]32   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
33   USE lib_mpp        ! Massively Parallel Processing library
[10642]34   USE iom             ! For shlat2d
35   USE fldread         ! for sn_shlat2d
[3]36
37   IMPLICIT NONE
38   PRIVATE
39
[6140]40   PUBLIC   dom_msk    ! routine called by inidom.F90
[3]41
[1601]42   !                            !!* Namelist namlbc : lateral boundary condition *
[4147]43   REAL(wp)        :: rn_shlat   ! type of lateral boundary condition on velocity
44   LOGICAL, PUBLIC :: ln_vorlat  !  consistency of vorticity boundary condition
[3294]45   !                                            with analytical eqs.
[2715]46
[3]47   !! * Substitutions
48#  include "vectopt_loop_substitute.h90"
[1566]49   !!----------------------------------------------------------------------
[9598]50   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
[9950]51   !! $Id$
[9598]52   !! Software governed by the CeCILL licence     (./LICENSE)
[1566]53   !!----------------------------------------------------------------------
[3]54CONTAINS
[2715]55
[7646]56   SUBROUTINE dom_msk( k_top, k_bot )
[3]57      !!---------------------------------------------------------------------
58      !!                 ***  ROUTINE dom_msk  ***
59      !!
60      !! ** Purpose :   Compute land/ocean mask arrays at tracer points, hori-
[6125]61      !!      zontal velocity points (u & v), vorticity points (f) points.
[3]62      !!
[7646]63      !! ** Method  :   The ocean/land mask  at t-point is deduced from ko_top
64      !!      and ko_bot, the indices of the fist and last ocean t-levels which
65      !!      are either defined in usrdef_zgr or read in zgr_read.
66      !!                The velocity masks (umask, vmask, wmask, wumask, wvmask)
67      !!      are deduced from a product of the two neighboring tmask.
68      !!                The vorticity mask (fmask) is deduced from tmask taking
69      !!      into account the choice of lateral boundary condition (rn_shlat) :
[1601]70      !!         rn_shlat = 0, free slip  (no shear along the coast)
71      !!         rn_shlat = 2, no slip  (specified zero velocity at the coast)
72      !!         0 < rn_shlat < 2, partial slip   | non-linear velocity profile
73      !!         2 < rn_shlat, strong slip        | in the lateral boundary layer
[3]74      !!
[7646]75      !!      tmask_i : interior ocean mask at t-point, i.e. excluding duplicated
76      !!                rows/lines due to cyclic or North Fold boundaries as well
77      !!                as MPP halos.
78      !!      tmask_h : halo mask at t-point, i.e. excluding duplicated rows/lines
79      !!                due to cyclic or North Fold boundaries as well as MPP halos.
[3]80      !!
[7646]81      !! ** Action :   tmask, umask, vmask, wmask, wumask, wvmask : land/ocean mask
82      !!                         at t-, u-, v- w, wu-, and wv-points (=0. or 1.)
83      !!               fmask   : land/ocean mask at f-point (=0., or =1., or
84      !!                         =rn_shlat along lateral boundaries)
85      !!               tmask_i : interior ocean mask
86      !!               tmask_h : halo mask
87      !!               ssmask , ssumask, ssvmask, ssfmask : 2D ocean mask
[3]88      !!----------------------------------------------------------------------
[7646]89      INTEGER, DIMENSION(:,:), INTENT(in) ::   k_top, k_bot   ! first and last ocean level
90      !
91      INTEGER  ::   ji, jj, jk     ! dummy loop indices
92      INTEGER  ::   iif, iil       ! local integers
93      INTEGER  ::   ijf, ijl       !   -       -
94      INTEGER  ::   iktop, ikbot   !   -       -
95      INTEGER  ::   ios, inum
[9019]96      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zwf   ! 2D workspace
[1601]97      !!
[10642]98      INTEGER  :: inum             !  logical unit for shlat2d
99      REAL(wp) :: zshlat           !: locally modified shlat for some strait
100      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zshlat2d
101      LOGICAL                         :: ln_shlat2d
102      CHARACTER(len = 256)            :: cn_shlat2d_file, cn_shlat2d_var 
103      !!
104      NAMELIST/namlbc/ rn_shlat, ln_vorlat, ln_shlat2d, cn_shlat2d_file, cn_shlat2d_var
[7646]105      NAMELIST/nambdy/ ln_bdy ,nb_bdy, ln_coords_file, cn_coords_file,         &
106         &             ln_mask_file, cn_mask_file, cn_dyn2d, nn_dyn2d_dta,     &
107         &             cn_dyn3d, nn_dyn3d_dta, cn_tra, nn_tra_dta,             &
108         &             ln_tra_dmp, ln_dyn3d_dmp, rn_time_dmp, rn_time_dmp_out, &
[9657]109         &             cn_ice, nn_ice_dta,                                     &
110         &             rn_ice_tem, rn_ice_sal, rn_ice_age,                     &
[7646]111         &             ln_vol, nn_volctl, nn_rimwidth, nb_jpk_bdy
[3]112      !!---------------------------------------------------------------------
[3294]113      !
[4147]114      REWIND( numnam_ref )              ! Namelist namlbc in reference namelist : Lateral momentum boundary condition
115      READ  ( numnam_ref, namlbc, IOSTAT = ios, ERR = 901 )
[9168]116901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namlbc in reference namelist', lwp )
[4147]117      REWIND( numnam_cfg )              ! Namelist namlbc in configuration namelist : Lateral momentum boundary condition
118      READ  ( numnam_cfg, namlbc, IOSTAT = ios, ERR = 902 )
[9168]119902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namlbc in configuration namelist', lwp )
[4624]120      IF(lwm) WRITE ( numond, namlbc )
[1566]121     
122      IF(lwp) THEN                  ! control print
[3]123         WRITE(numout,*)
124         WRITE(numout,*) 'dommsk : ocean mask '
125         WRITE(numout,*) '~~~~~~'
[1566]126         WRITE(numout,*) '   Namelist namlbc'
[3294]127         WRITE(numout,*) '      lateral momentum boundary cond.    rn_shlat  = ',rn_shlat
128         WRITE(numout,*) '      consistency with analytical form   ln_vorlat = ',ln_vorlat 
[3]129      ENDIF
[9169]130      !
131      IF(lwp) WRITE(numout,*)
132      IF     (      rn_shlat == 0.               ) THEN   ;   IF(lwp) WRITE(numout,*) '   ==>>>   ocean lateral  free-slip'
133      ELSEIF (      rn_shlat == 2.               ) THEN   ;   IF(lwp) WRITE(numout,*) '   ==>>>   ocean lateral  no-slip'
134      ELSEIF ( 0. < rn_shlat .AND. rn_shlat < 2. ) THEN   ;   IF(lwp) WRITE(numout,*) '   ==>>>   ocean lateral  partial-slip'
135      ELSEIF ( 2. < rn_shlat                     ) THEN   ;   IF(lwp) WRITE(numout,*) '   ==>>>   ocean lateral  strong-slip'
[1601]136      ELSE
[9527]137         CALL ctl_stop( 'dom_msk: wrong value for rn_shlat (i.e. a negalive value). We stop.' )
[3]138      ENDIF
139
[10642]140      IF ( ln_shlat2d ) THEN
141         IF(lwp) WRITE(numout,*) '         READ shlat as a 2D coefficient in a file '
142         ALLOCATE( zshlat2d(jpi,jpj) )
143         CALL iom_open(TRIM(cn_shlat2d_file), inum)
144         CALL iom_get (inum, jpdom_data, TRIM(cn_shlat2d_var), zshlat2d, 1) !
145         CALL iom_close(inum)
146      ENDIF
147
[7646]148      !  Ocean/land mask at t-point  (computed from ko_top and ko_bot)
149      ! ----------------------------
[1566]150      !
[7753]151      tmask(:,:,:) = 0._wp
[7646]152      DO jj = 1, jpj
153         DO ji = 1, jpi
154            iktop = k_top(ji,jj)
155            ikbot = k_bot(ji,jj)
156            IF( iktop /= 0 ) THEN       ! water in the column
157               tmask(ji,jj,iktop:ikbot  ) = 1._wp
158            ENDIF
[3]159         END DO 
160      END DO 
[7646]161!SF  add here lbc_lnk: bug not still understood : cause now domain configuration is read !
162!!gm I don't understand why... 
163      CALL lbc_lnk( tmask  , 'T', 1._wp )      ! Lateral boundary conditions
[3]164
[7646]165     ! Mask corrections for bdy (read in mppini2)
166      REWIND( numnam_ref )              ! Namelist nambdy in reference namelist :Unstructured open boundaries
167      READ  ( numnam_ref, nambdy, IOSTAT = ios, ERR = 903)
[9168]168903   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambdy in reference namelist', lwp )
[7646]169      REWIND( numnam_cfg )              ! Namelist nambdy in configuration namelist :Unstructured open boundaries
170      READ  ( numnam_cfg, nambdy, IOSTAT = ios, ERR = 904 )
[9168]171904   IF( ios >  0 )   CALL ctl_nam ( ios , 'nambdy in configuration namelist', lwp )
[7646]172      ! ------------------------
173      IF ( ln_bdy .AND. ln_mask_file ) THEN
[9600]174         CALL iom_open( cn_mask_file, inum )
175         CALL iom_get ( inum, jpdom_data, 'bdy_msk', bdytmask(:,:) )
176         CALL iom_close( inum )
[7646]177         DO jk = 1, jpkm1
178            DO jj = 1, jpj
179               DO ji = 1, jpi
180                  tmask(ji,jj,jk) = tmask(ji,jj,jk) * bdytmask(ji,jj)
181               END DO
[291]182            END DO
[7646]183         END DO
[3]184      ENDIF
[7646]185         
186      !  Ocean/land mask at u-, v-, and f-points   (computed from tmask)
187      ! ----------------------------------------
188      ! NB: at this point, fmask is designed for free slip lateral boundary condition
[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
[9090]201      CALL lbc_lnk_multi( umask, 'U', 1., vmask, 'V', 1., fmask, 'F', 1. )      ! Lateral boundary conditions
[7646]202 
203      ! Ocean/land mask at wu-, wv- and w points    (computed from tmask)
204      !-----------------------------------------
[7753]205      wmask (:,:,1) = tmask(:,:,1)     ! surface
206      wumask(:,:,1) = umask(:,:,1)
207      wvmask(:,:,1) = vmask(:,:,1)
[5836]208      DO jk = 2, jpk                   ! interior values
[7753]209         wmask (:,:,jk) = tmask(:,:,jk) * tmask(:,:,jk-1)
210         wumask(:,:,jk) = umask(:,:,jk) * umask(:,:,jk-1)   
211         wvmask(:,:,jk) = vmask(:,:,jk) * vmask(:,:,jk-1)
[5120]212      END DO
[3]213
[7646]214
215      ! Ocean/land column mask at t-, u-, and v-points   (i.e. at least 1 wet cell in the vertical)
216      ! ----------------------------------------------
217      ssmask (:,:) = MAXVAL( tmask(:,:,:), DIM=3 )
218      ssumask(:,:) = MAXVAL( umask(:,:,:), DIM=3 )
219      ssvmask(:,:) = MAXVAL( vmask(:,:,:), DIM=3 )
220
221
222      ! Interior domain mask  (used for global sum)
223      ! --------------------
224      !
[9019]225      iif = nn_hls   ;   iil = nlci - nn_hls + 1
226      ijf = nn_hls   ;   ijl = nlcj - nn_hls + 1
[7646]227      !
228      !                          ! halo mask : 0 on the halo and 1 elsewhere
[7753]229      tmask_h(:,:) = 1._wp                 
[7646]230      tmask_h( 1 :iif,   :   ) = 0._wp      ! first columns
231      tmask_h(iil:jpi,   :   ) = 0._wp      ! last  columns (including mpp extra columns)
232      tmask_h(   :   , 1 :ijf) = 0._wp      ! first rows
233      tmask_h(   :   ,ijl:jpj) = 0._wp      ! last  rows (including mpp extra rows)
234      !
235      !                          ! north fold mask
236      tpol(1:jpiglo) = 1._wp 
237      fpol(1:jpiglo) = 1._wp
238      IF( jperio == 3 .OR. jperio == 4 ) THEN      ! T-point pivot
239         tpol(jpiglo/2+1:jpiglo) = 0._wp
240         fpol(     1    :jpiglo) = 0._wp
241         IF( mjg(nlej) == jpjglo ) THEN                  ! only half of the nlcj-1 row for tmask_h
242            DO ji = iif+1, iil-1
243               tmask_h(ji,nlej-1) = tmask_h(ji,nlej-1) * tpol(mig(ji))
244            END DO
245         ENDIF
246      ENDIF
247      !
248      IF( jperio == 5 .OR. jperio == 6 ) THEN      ! F-point pivot
249         tpol(     1    :jpiglo) = 0._wp
250         fpol(jpiglo/2+1:jpiglo) = 0._wp
251      ENDIF
252      !
253      !                          ! interior mask : 2D ocean mask x halo mask
[7753]254      tmask_i(:,:) = ssmask(:,:) * tmask_h(:,:)
[7646]255
256
[3]257      ! Lateral boundary conditions on velocity (modify fmask)
[7646]258      ! --------------------------------------- 
[10642]259      IF( rn_shlat /= 0 .or. ln_shlat2d ) THEN      ! Not free-slip lateral boundary condition everywhere
[7646]260         !
[9019]261         ALLOCATE( zwf(jpi,jpj) )
[7646]262         !
263         DO jk = 1, jpk
[7753]264            zwf(:,:) = fmask(:,:,jk)         
[10642]265            IF (  ln_shlat2d ) THEN
266               DO jj = 2, jpjm1
267                  DO ji = fs_2, fs_jpim1   ! vector opt.
268                     IF( fmask(ji,jj,jk) == 0._wp ) THEN
269                        fmask(ji,jj,jk) = zshlat2d(ji,jj) * MIN( 1._wp , MAX( zwf(ji+1,jj), zwf(ji,jj+1),   &
270                           &                                           zwf(ji-1,jj), zwf(ji,jj-1)  )  )
271                     ENDIF
272                  END DO
[7646]273               END DO
[10642]274            ELSE
275               DO jj = 2, jpjm1
276                  DO ji = fs_2, fs_jpim1   ! vector opt.
277                     IF( fmask(ji,jj,jk) == 0._wp ) THEN
278                        fmask(ji,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,jj), zwf(ji,jj+1),   &
279                           &                                           zwf(ji-1,jj), zwf(ji,jj-1)  )  )
280                     ENDIF
281                  END DO
282               END DO
283            ENDIF
[7646]284            DO jj = 2, jpjm1
285               IF( fmask(1,jj,jk) == 0._wp ) THEN
286                  fmask(1  ,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1) ) )
[3]287               ENDIF
[7646]288               IF( fmask(jpi,jj,jk) == 0._wp ) THEN
289                  fmask(jpi,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) )
290               ENDIF
291            END DO         
292            DO ji = 2, jpim1
293               IF( fmask(ji,1,jk) == 0._wp ) THEN
294                  fmask(ji, 1 ,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1) ) )
295               ENDIF
296               IF( fmask(ji,jpj,jk) == 0._wp ) THEN
297                  fmask(ji,jpj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpjm1) ) )
298               ENDIF
[3]299            END DO
[9023]300#if defined key_agrif 
[9090]301            IF( .NOT. AGRIF_Root() ) THEN
302               IF ((nbondi ==  1).OR.(nbondi == 2)) fmask(nlci-1 , :     ,jk) = 0.e0      ! east
303               IF ((nbondi == -1).OR.(nbondi == 2)) fmask(1      , :     ,jk) = 0.e0      ! west
304               IF ((nbondj ==  1).OR.(nbondj == 2)) fmask(:      ,nlcj-1 ,jk) = 0.e0      ! north
305               IF ((nbondj == -1).OR.(nbondj == 2)) fmask(:      ,1      ,jk) = 0.e0      ! south
306            ENDIF 
[9023]307#endif
[3]308         END DO
[5836]309         !
[9019]310         DEALLOCATE( zwf )
[10642]311         IF( ln_shlat2d ) DEALLOCATE( zshlat2d )
[5836]312         !
[7646]313         CALL lbc_lnk( fmask, 'F', 1._wp )      ! Lateral boundary conditions on fmask
[3]314         !
[7646]315         ! CAUTION : The fmask may be further modified in dyn_vor_init ( dynvor.F90 ) depending on ln_vorlat
[5506]316         !
[2528]317      ENDIF
[7646]318     
319      ! User defined alteration of fmask (use to reduce ocean transport in specified straits)
320      ! --------------------------------
[2528]321      !
[7646]322      CALL usr_def_fmask( cn_cfg, nn_cfg, fmask )
[6140]323      !
[3]324   END SUBROUTINE dom_msk
325   
326   !!======================================================================
327END MODULE dommsk
Note: See TracBrowser for help on using the repository browser.