source: NEMO/branches/2019/ENHANCE-03_domcfg/src/dommsk.F90 @ 11129

Last change on this file since 11129 was 11129, checked in by mathiot, 15 months ago

simplification of domcfg (rm all var_n and var_b as it is not needed) (ticket #2143)

File size: 14.4 KB
Line 
1MODULE dommsk
2   !!======================================================================
3   !!                       ***  MODULE dommsk   ***
4   !! Ocean initialization : domain land/sea mask
5   !!======================================================================
6   !! History :  OPA  ! 1987-07  (G. Madec)  Original code
7   !!            6.0  ! 1993-03  (M. Guyon)  symetrical conditions (M. Guyon)
8   !!            7.0  ! 1996-01  (G. Madec)  suppression of common work arrays
9   !!             -   ! 1996-05  (G. Madec)  mask computed from tmask
10   !!            8.0  ! 1997-02  (G. Madec)  mesh information put in domhgr.F
11   !!            8.1  ! 1997-07  (G. Madec)  modification of kbat and fmask
12   !!             -   ! 1998-05  (G. Roullet)  free surface
13   !!            8.2  ! 2000-03  (G. Madec)  no slip accurate
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
18   !!            3.6  ! 2015-05  (P. Mathiot) ISF: add wmask,wumask and wvmask
19   !!            4.0  ! 2016-06  (G. Madec, S. Flavoni)  domain configuration / user defined interface
20   !!----------------------------------------------------------------------
21
22   !!----------------------------------------------------------------------
23   !!   dom_msk       : compute land/ocean mask
24   !!----------------------------------------------------------------------
25   USE dom_oce        ! ocean space and time domain
26   USE bdy_oce        ! open boundary
27   !
28   USE in_out_manager ! I/O manager
29   USE iom            ! IOM library
30   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
31   USE lib_mpp        ! Massively Parallel Processing library
32
33   IMPLICIT NONE
34   PRIVATE
35
36   PUBLIC   dom_msk    ! routine called by inidom.F90
37
38   !                            !!* Namelist namlbc : lateral boundary condition *
39   REAL(wp)        :: rn_shlat   ! type of lateral boundary condition on velocity
40   LOGICAL, PUBLIC :: ln_vorlat  !  consistency of vorticity boundary condition
41   !                                            with analytical eqs.
42
43   !! * Substitutions
44#  include "vectopt_loop_substitute.h90"
45   !!----------------------------------------------------------------------
46   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
47   !! $Id: dommsk.F90 10425 2018-12-19 21:54:16Z smasson $
48   !! Software governed by the CeCILL license (see ./LICENSE)
49   !!----------------------------------------------------------------------
50CONTAINS
51
52   SUBROUTINE dom_msk
53      !!---------------------------------------------------------------------
54      !!                 ***  ROUTINE dom_msk  ***
55      !!
56      !! ** Purpose :   Compute land/ocean mask arrays at tracer points, hori-
57      !!      zontal velocity points (u & v), vorticity points (f) points.
58      !!
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) :
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
70      !!
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 as MPP halos.
76      !!
77      !! ** Action :   tmask, umask, vmask, wmask, wumask, wvmask : land/ocean mask
78      !!                         at t-, u-, v- w, wu-, and wv-points (=0. or 1.)
79      !!               fmask   : land/ocean mask at f-point (=0., or =1., or
80      !!                         =rn_shlat along lateral boundaries)
81      !!               tmask_i : interior ocean mask
82      !!               tmask_h : halo mask
83      !!               ssmask , ssumask, ssvmask, ssfmask : 2D ocean mask
84      !!----------------------------------------------------------------------
85      !
86      INTEGER  ::   ji, jj, jk     ! dummy loop indices
87      INTEGER  ::   iif, iil       ! local integers
88      INTEGER  ::   ijf, ijl       !   -       -
89      INTEGER  ::   iktop, ikbot   !   -       -
90      INTEGER  ::   ios, inum
91      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zwf   ! 2D workspace
92      !!
93      NAMELIST/namlbc/ rn_shlat, ln_vorlat
94      NAMELIST/nambdy/ ln_bdy ,nb_bdy, ln_coords_file, cn_coords_file,         &
95         &             ln_mask_file, cn_mask_file, cn_dyn2d, nn_dyn2d_dta,     &
96         &             cn_dyn3d, nn_dyn3d_dta, cn_tra, nn_tra_dta,             &
97         &             ln_tra_dmp, ln_dyn3d_dmp, rn_time_dmp, rn_time_dmp_out, &
98         &             cn_ice, nn_ice_dta,                                     &
99         &             rn_ice_tem, rn_ice_sal, rn_ice_age,                     &
100         &             ln_vol, nn_volctl, nn_rimwidth, nb_jpk_bdy
101      !!---------------------------------------------------------------------
102      !
103      REWIND( numnam_ref )              ! Namelist namlbc in reference namelist : Lateral momentum boundary condition
104      READ  ( numnam_ref, namlbc, IOSTAT = ios, ERR = 901 )
105901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namlbc in reference namelist', lwp )
106      REWIND( numnam_cfg )              ! Namelist namlbc in configuration namelist : Lateral momentum boundary condition
107      READ  ( numnam_cfg, namlbc, IOSTAT = ios, ERR = 902 )
108902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namlbc in configuration namelist', lwp )
109      IF(lwm) WRITE ( numond, namlbc )
110     
111      IF(lwp) THEN                  ! control print
112         WRITE(numout,*)
113         WRITE(numout,*) 'dommsk : ocean mask '
114         WRITE(numout,*) '~~~~~~'
115         WRITE(numout,*) '   Namelist namlbc'
116         WRITE(numout,*) '      lateral momentum boundary cond.    rn_shlat  = ',rn_shlat
117         WRITE(numout,*) '      consistency with analytical form   ln_vorlat = ',ln_vorlat 
118      ENDIF
119      !
120      IF(lwp) WRITE(numout,*)
121      IF     (      rn_shlat == 0.               ) THEN   ;   IF(lwp) WRITE(numout,*) '   ==>>>   ocean lateral  free-slip'
122      ELSEIF (      rn_shlat == 2.               ) THEN   ;   IF(lwp) WRITE(numout,*) '   ==>>>   ocean lateral  no-slip'
123      ELSEIF ( 0. < rn_shlat .AND. rn_shlat < 2. ) THEN   ;   IF(lwp) WRITE(numout,*) '   ==>>>   ocean lateral  partial-slip'
124      ELSEIF ( 2. < rn_shlat                     ) THEN   ;   IF(lwp) WRITE(numout,*) '   ==>>>   ocean lateral  strong-slip'
125      ELSE
126         CALL ctl_stop( 'dom_msk: wrong value for rn_shlat (i.e. a negalive value). We stop.' )
127      ENDIF
128
129     ! 1. Ocean/land mask at t-point (computed from mbathy)
130     ! -----------------------------
131     ! N.B. tmask has already the right boundary conditions since mbathy is ok
132     !
133      tmask(:,:,:) = 0._wp
134      DO jk = 1, jpk
135         DO jj = 1, jpj
136            DO ji = 1, jpi
137               IF(      ( REAL( mbathy (ji,jj) - jk, wp ) + 0.1_wp >= 0._wp )         &
138               &  .AND. ( REAL( misfdep(ji,jj) - jk, wp ) - 0.1_wp <= 0._wp ) ) THEN
139                  tmask(ji,jj,jk) = 1._wp
140               END IF 
141            END DO
142         END DO
143      END DO
144     
145!SF  add here lbc_lnk: bug not still understood : cause now domain configuration is read !
146!!gm I don't understand why... 
147      CALL lbc_lnk( 'dommsk', tmask  , 'T', 1._wp )      ! Lateral boundary conditions
148
149     ! Mask corrections for bdy (read in mppini2)
150      REWIND( numnam_ref )              ! Namelist nambdy in reference namelist :Unstructured open boundaries
151      READ  ( numnam_ref, nambdy, IOSTAT = ios, ERR = 903)
152903   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambdy in reference namelist', lwp )
153      REWIND( numnam_cfg )              ! Namelist nambdy in configuration namelist :Unstructured open boundaries
154      READ  ( numnam_cfg, nambdy, IOSTAT = ios, ERR = 904 )
155904   IF( ios >  0 )   CALL ctl_nam ( ios , 'nambdy in configuration namelist', lwp )
156      ! ------------------------
157      IF ( ln_bdy .AND. ln_mask_file ) THEN
158         CALL iom_open( cn_mask_file, inum )
159         CALL iom_get ( inum, jpdom_data, 'bdy_msk', bdytmask(:,:) )
160         CALL iom_close( inum )
161         DO jk = 1, jpkm1
162            DO jj = 1, jpj
163               DO ji = 1, jpi
164                  tmask(ji,jj,jk) = tmask(ji,jj,jk) * bdytmask(ji,jj)
165               END DO
166            END DO
167         END DO
168      ENDIF
169         
170      !  Ocean/land mask at u-, v-, and f-points   (computed from tmask)
171      ! ----------------------------------------
172      ! NB: at this point, fmask is designed for free slip lateral boundary condition
173      DO jk = 1, jpk
174         DO jj = 1, jpjm1
175            DO ji = 1, fs_jpim1   ! vector loop
176               umask(ji,jj,jk) = tmask(ji,jj  ,jk) * tmask(ji+1,jj  ,jk)
177               vmask(ji,jj,jk) = tmask(ji,jj  ,jk) * tmask(ji  ,jj+1,jk)
178            END DO
179            DO ji = 1, jpim1      ! NO vector opt.
180               fmask(ji,jj,jk) = tmask(ji,jj  ,jk) * tmask(ji+1,jj  ,jk)   &
181                  &            * tmask(ji,jj+1,jk) * tmask(ji+1,jj+1,jk)
182            END DO
183         END DO
184      END DO
185      CALL lbc_lnk_multi( 'dommsk', umask, 'U', 1., vmask, 'V', 1., fmask, 'F', 1. )      ! Lateral boundary conditions
186 
187      ! Ocean/land mask at wu-, wv- and w points    (computed from tmask)
188      !-----------------------------------------
189      wmask (:,:,1) = tmask(:,:,1)     ! surface
190      wumask(:,:,1) = umask(:,:,1)
191      wvmask(:,:,1) = vmask(:,:,1)
192      DO jk = 2, jpk                   ! interior values
193         wmask (:,:,jk) = tmask(:,:,jk) * tmask(:,:,jk-1)
194         wumask(:,:,jk) = umask(:,:,jk) * umask(:,:,jk-1)   
195         wvmask(:,:,jk) = vmask(:,:,jk) * vmask(:,:,jk-1)
196      END DO
197
198
199      ! Ocean/land column mask at t-, u-, and v-points   (i.e. at least 1 wet cell in the vertical)
200      ! ----------------------------------------------
201      ssmask (:,:) = MAXVAL( tmask(:,:,:), DIM=3 )
202      ssumask(:,:) = MAXVAL( umask(:,:,:), DIM=3 )
203      ssvmask(:,:) = MAXVAL( vmask(:,:,:), DIM=3 )
204
205
206      ! Interior domain mask  (used for global sum)
207      ! --------------------
208      !
209      iif = nn_hls   ;   iil = nlci - nn_hls + 1
210      ijf = nn_hls   ;   ijl = nlcj - nn_hls + 1
211      !
212      !                          ! halo mask : 0 on the halo and 1 elsewhere
213      tmask_h(:,:) = 1._wp                 
214      tmask_h( 1 :iif,   :   ) = 0._wp      ! first columns
215      tmask_h(iil:jpi,   :   ) = 0._wp      ! last  columns (including mpp extra columns)
216      tmask_h(   :   , 1 :ijf) = 0._wp      ! first rows
217      tmask_h(   :   ,ijl:jpj) = 0._wp      ! last  rows (including mpp extra rows)
218      !
219      !                          ! north fold mask
220      tpol(1:jpiglo) = 1._wp 
221      fpol(1:jpiglo) = 1._wp
222      IF( jperio == 3 .OR. jperio == 4 ) THEN      ! T-point pivot
223         tpol(jpiglo/2+1:jpiglo) = 0._wp
224         fpol(     1    :jpiglo) = 0._wp
225         IF( mjg(nlej) == jpjglo ) THEN                  ! only half of the nlcj-1 row for tmask_h
226            DO ji = iif+1, iil-1
227               tmask_h(ji,nlej-1) = tmask_h(ji,nlej-1) * tpol(mig(ji))
228            END DO
229         ENDIF
230      ENDIF
231      !
232      IF( jperio == 5 .OR. jperio == 6 ) THEN      ! F-point pivot
233         tpol(     1    :jpiglo) = 0._wp
234         fpol(jpiglo/2+1:jpiglo) = 0._wp
235      ENDIF
236      !
237      !                          ! interior mask : 2D ocean mask x halo mask
238      tmask_i(:,:) = ssmask(:,:) * tmask_h(:,:)
239
240
241      ! Lateral boundary conditions on velocity (modify fmask)
242      ! --------------------------------------- 
243      IF( rn_shlat /= 0 ) THEN      ! Not free-slip lateral boundary condition
244         !
245         ALLOCATE( zwf(jpi,jpj) )
246         !
247         DO jk = 1, jpk
248            zwf(:,:) = fmask(:,:,jk)         
249            DO jj = 2, jpjm1
250               DO ji = fs_2, fs_jpim1   ! vector opt.
251                  IF( fmask(ji,jj,jk) == 0._wp ) THEN
252                     fmask(ji,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,jj), zwf(ji,jj+1),   &
253                        &                                           zwf(ji-1,jj), zwf(ji,jj-1)  )  )
254                  ENDIF
255               END DO
256            END DO
257            DO jj = 2, jpjm1
258               IF( fmask(1,jj,jk) == 0._wp ) THEN
259                  fmask(1  ,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1) ) )
260               ENDIF
261               IF( fmask(jpi,jj,jk) == 0._wp ) THEN
262                  fmask(jpi,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) )
263               ENDIF
264            END DO         
265            DO ji = 2, jpim1
266               IF( fmask(ji,1,jk) == 0._wp ) THEN
267                  fmask(ji, 1 ,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1) ) )
268               ENDIF
269               IF( fmask(ji,jpj,jk) == 0._wp ) THEN
270                  fmask(ji,jpj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpjm1) ) )
271               ENDIF
272            END DO
273#if defined key_agrif 
274            IF( .NOT. AGRIF_Root() ) THEN
275               IF ((nbondi ==  1).OR.(nbondi == 2)) fmask(nlci-1 , :     ,jk) = 0.e0      ! east
276               IF ((nbondi == -1).OR.(nbondi == 2)) fmask(1      , :     ,jk) = 0.e0      ! west
277               IF ((nbondj ==  1).OR.(nbondj == 2)) fmask(:      ,nlcj-1 ,jk) = 0.e0      ! north
278               IF ((nbondj == -1).OR.(nbondj == 2)) fmask(:      ,1      ,jk) = 0.e0      ! south
279            ENDIF 
280#endif
281         END DO
282         !
283         DEALLOCATE( zwf )
284         !
285         CALL lbc_lnk( 'dommsk', fmask, 'F', 1._wp )      ! Lateral boundary conditions on fmask
286         !
287         ! CAUTION : The fmask may be further modified in dyn_vor_init ( dynvor.F90 ) depending on ln_vorlat
288         !
289      ENDIF
290      !
291   END SUBROUTINE dom_msk
292   
293   !!======================================================================
294END MODULE dommsk
Note: See TracBrowser for help on using the repository browser.