source: trunk/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90 @ 7698

Last change on this file since 7698 was 7698, checked in by mocavero, 4 years ago

update trunk with OpenMP parallelization

  • Property svn:keywords set to Id
File size: 15.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 oce            ! ocean dynamics and tracers
26   USE dom_oce        ! ocean space and time domain
27   USE usrdef_fmask   ! user defined fmask
28   USE bdy_oce     
29   USE in_out_manager ! I/O manager
30   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
31   USE lib_mpp        ! Massively Parallel Processing library
32   USE wrk_nemo       ! Memory allocation
33   USE timing         ! Timing
34
35   IMPLICIT NONE
36   PRIVATE
37
38   PUBLIC   dom_msk    ! routine called by inidom.F90
39
40   !                            !!* Namelist namlbc : lateral boundary condition *
41   REAL(wp)        :: rn_shlat   ! type of lateral boundary condition on velocity
42   LOGICAL, PUBLIC :: ln_vorlat  !  consistency of vorticity boundary condition
43   !                                            with analytical eqs.
44
45   !! * Substitutions
46#  include "vectopt_loop_substitute.h90"
47   !!----------------------------------------------------------------------
48   !! NEMO/OPA 3.2 , LODYC-IPSL  (2009)
49   !! $Id$
50   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
51   !!----------------------------------------------------------------------
52CONTAINS
53
54   SUBROUTINE dom_msk( k_top, k_bot )
55      !!---------------------------------------------------------------------
56      !!                 ***  ROUTINE dom_msk  ***
57      !!
58      !! ** Purpose :   Compute land/ocean mask arrays at tracer points, hori-
59      !!      zontal velocity points (u & v), vorticity points (f) points.
60      !!
61      !! ** Method  :   The ocean/land mask  at t-point is deduced from ko_top
62      !!      and ko_bot, the indices of the fist and last ocean t-levels which
63      !!      are either defined in usrdef_zgr or read in zgr_read.
64      !!                The velocity masks (umask, vmask, wmask, wumask, wvmask)
65      !!      are deduced from a product of the two neighboring tmask.
66      !!                The vorticity mask (fmask) is deduced from tmask taking
67      !!      into account the choice of lateral boundary condition (rn_shlat) :
68      !!         rn_shlat = 0, free slip  (no shear along the coast)
69      !!         rn_shlat = 2, no slip  (specified zero velocity at the coast)
70      !!         0 < rn_shlat < 2, partial slip   | non-linear velocity profile
71      !!         2 < rn_shlat, strong slip        | in the lateral boundary layer
72      !!
73      !!      tmask_i : interior ocean mask at t-point, i.e. excluding duplicated
74      !!                rows/lines due to cyclic or North Fold boundaries as well
75      !!                as MPP halos.
76      !!      tmask_h : halo mask at t-point, i.e. excluding duplicated rows/lines
77      !!                due to cyclic or North Fold boundaries as well as MPP halos.
78      !!
79      !! ** Action :   tmask, umask, vmask, wmask, wumask, wvmask : land/ocean mask
80      !!                         at t-, u-, v- w, wu-, and wv-points (=0. or 1.)
81      !!               fmask   : land/ocean mask at f-point (=0., or =1., or
82      !!                         =rn_shlat along lateral boundaries)
83      !!               tmask_i : interior ocean mask
84      !!               tmask_h : halo mask
85      !!               ssmask , ssumask, ssvmask, ssfmask : 2D ocean mask
86      !!----------------------------------------------------------------------
87      INTEGER, DIMENSION(:,:), INTENT(in) ::   k_top, k_bot   ! first and last ocean level
88      !
89      INTEGER  ::   ji, jj, jk     ! dummy loop indices
90      INTEGER  ::   iif, iil       ! local integers
91      INTEGER  ::   ijf, ijl       !   -       -
92      INTEGER  ::   iktop, ikbot   !   -       -
93      INTEGER  ::   ios, inum
94      REAL(wp), POINTER, DIMENSION(:,:) ::   zwf   ! 2D workspace
95      !!
96      NAMELIST/namlbc/ rn_shlat, ln_vorlat
97      NAMELIST/nambdy/ ln_bdy ,nb_bdy, ln_coords_file, cn_coords_file,         &
98         &             ln_mask_file, cn_mask_file, cn_dyn2d, nn_dyn2d_dta,     &
99         &             cn_dyn3d, nn_dyn3d_dta, cn_tra, nn_tra_dta,             &
100         &             ln_tra_dmp, ln_dyn3d_dmp, rn_time_dmp, rn_time_dmp_out, &
101         &             cn_ice_lim, nn_ice_lim_dta,                           &
102         &             rn_ice_tem, rn_ice_sal, rn_ice_age,                 &
103         &             ln_vol, nn_volctl, nn_rimwidth, nb_jpk_bdy
104      !!---------------------------------------------------------------------
105      !
106      IF( nn_timing == 1 )  CALL timing_start('dom_msk')
107      !
108      REWIND( numnam_ref )              ! Namelist namlbc in reference namelist : Lateral momentum boundary condition
109      READ  ( numnam_ref, namlbc, IOSTAT = ios, ERR = 901 )
110901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namlbc in reference namelist', lwp )
111
112      REWIND( numnam_cfg )              ! Namelist namlbc in configuration namelist : Lateral momentum boundary condition
113      READ  ( numnam_cfg, namlbc, IOSTAT = ios, ERR = 902 )
114902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namlbc in configuration namelist', lwp )
115      IF(lwm) WRITE ( numond, namlbc )
116     
117      IF(lwp) THEN                  ! control print
118         WRITE(numout,*)
119         WRITE(numout,*) 'dommsk : ocean mask '
120         WRITE(numout,*) '~~~~~~'
121         WRITE(numout,*) '   Namelist namlbc'
122         WRITE(numout,*) '      lateral momentum boundary cond.    rn_shlat  = ',rn_shlat
123         WRITE(numout,*) '      consistency with analytical form   ln_vorlat = ',ln_vorlat 
124      ENDIF
125
126      IF     (      rn_shlat == 0.               ) THEN   ;   IF(lwp) WRITE(numout,*) '   ocean lateral  free-slip '
127      ELSEIF (      rn_shlat == 2.               ) THEN   ;   IF(lwp) WRITE(numout,*) '   ocean lateral  no-slip '
128      ELSEIF ( 0. < rn_shlat .AND. rn_shlat < 2. ) THEN   ;   IF(lwp) WRITE(numout,*) '   ocean lateral  partial-slip '
129      ELSEIF ( 2. < rn_shlat                     ) THEN   ;   IF(lwp) WRITE(numout,*) '   ocean lateral  strong-slip '
130      ELSE
131         WRITE(ctmp1,*) ' rn_shlat is negative = ', rn_shlat
132         CALL ctl_stop( ctmp1 )
133      ENDIF
134
135
136      !  Ocean/land mask at t-point  (computed from ko_top and ko_bot)
137      ! ----------------------------
138      !
139!$OMP PARALLEL
140!$OMP DO schedule(static) private(jk, jj, ji)
141      DO jk = 1, jpk
142         DO jj = 1, jpj
143            DO ji = 1, jpi
144               tmask(ji,jj,jk) = 0._wp
145            END DO
146         END DO
147      END DO
148!$OMP DO schedule(static) private(jj, ji, iktop, ikbot)
149      DO jj = 1, jpj
150         DO ji = 1, jpi
151            iktop = k_top(ji,jj)
152            ikbot = k_bot(ji,jj)
153            IF( iktop /= 0 ) THEN       ! water in the column
154               tmask(ji,jj,iktop:ikbot  ) = 1._wp
155            ENDIF
156         END DO 
157      END DO 
158!$OMP END PARALLEL
159!SF  add here lbc_lnk: bug not still understood : cause now domain configuration is read !
160!!gm I don't understand why... 
161      CALL lbc_lnk( tmask  , 'T', 1._wp )      ! Lateral boundary conditions
162
163     ! Mask corrections for bdy (read in mppini2)
164      REWIND( numnam_ref )              ! Namelist nambdy in reference namelist :Unstructured open boundaries
165      READ  ( numnam_ref, nambdy, IOSTAT = ios, ERR = 903)
166903   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambdy in reference namelist', lwp )
167
168      REWIND( numnam_cfg )              ! Namelist nambdy in configuration namelist :Unstructured open boundaries
169      READ  ( numnam_cfg, nambdy, IOSTAT = ios, ERR = 904 )
170904   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambdy in configuration namelist', lwp )
171      ! ------------------------
172      IF ( ln_bdy .AND. ln_mask_file ) THEN
173!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
174         DO jk = 1, jpkm1
175            DO jj = 1, jpj
176               DO ji = 1, jpi
177                  tmask(ji,jj,jk) = tmask(ji,jj,jk) * bdytmask(ji,jj)
178               END DO
179            END DO
180         END DO
181      ENDIF
182         
183      !  Ocean/land mask at u-, v-, and f-points   (computed from tmask)
184      ! ----------------------------------------
185      ! NB: at this point, fmask is designed for free slip lateral boundary condition
186!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
187      DO jk = 1, jpk
188         DO jj = 1, jpjm1
189            DO ji = 1, fs_jpim1   ! vector loop
190               umask(ji,jj,jk) = tmask(ji,jj  ,jk) * tmask(ji+1,jj  ,jk)
191               vmask(ji,jj,jk) = tmask(ji,jj  ,jk) * tmask(ji  ,jj+1,jk)
192            END DO
193            DO ji = 1, jpim1      ! NO vector opt.
194               fmask(ji,jj,jk) = tmask(ji,jj  ,jk) * tmask(ji+1,jj  ,jk)   &
195                  &            * tmask(ji,jj+1,jk) * tmask(ji+1,jj+1,jk)
196            END DO
197         END DO
198      END DO
199      CALL lbc_lnk( umask  , 'U', 1._wp )      ! Lateral boundary conditions
200      CALL lbc_lnk( vmask  , 'V', 1._wp )
201      CALL lbc_lnk( fmask  , 'F', 1._wp )
202
203 
204      ! Ocean/land mask at wu-, wv- and w points    (computed from tmask)
205      !-----------------------------------------
206!$OMP PARALLEL
207!$OMP DO schedule(static) private(jj, ji)
208      DO jj = 1, jpj
209         DO ji = 1, jpi
210            wmask (ji,jj,1) = tmask(ji,jj,1)     ! surface
211            wumask(ji,jj,1) = umask(ji,jj,1)
212            wvmask(ji,jj,1) = vmask(ji,jj,1)
213         END DO
214      END DO
215!$OMP DO schedule(static) private(jk,jj,ji)
216      DO jk = 2, jpk                   ! interior values
217         DO jj = 1, jpj
218            DO ji = 1, jpi
219               wmask (ji,jj,jk) = tmask(ji,jj,jk) * tmask(ji,jj,jk-1)
220               wumask(ji,jj,jk) = umask(ji,jj,jk) * umask(ji,jj,jk-1)   
221               wvmask(ji,jj,jk) = vmask(ji,jj,jk) * vmask(ji,jj,jk-1)
222            END DO
223         END DO
224      END DO
225!$OMP END PARALLEL
226
227
228      ! Ocean/land column mask at t-, u-, and v-points   (i.e. at least 1 wet cell in the vertical)
229      ! ----------------------------------------------
230      ssmask (:,:) = MAXVAL( tmask(:,:,:), DIM=3 )
231      ssumask(:,:) = MAXVAL( umask(:,:,:), DIM=3 )
232      ssvmask(:,:) = MAXVAL( vmask(:,:,:), DIM=3 )
233
234
235      ! Interior domain mask  (used for global sum)
236      ! --------------------
237      !
238      iif = jpreci   ;   iil = nlci - jpreci + 1
239      ijf = jprecj   ;   ijl = nlcj - jprecj + 1
240      !
241      !                          ! halo mask : 0 on the halo and 1 elsewhere
242!$OMP PARALLEL DO schedule(static) private(jj, ji)
243      DO jj = 1, jpj
244         DO ji = 1, jpi
245            tmask_h(ji,jj) = 1._wp                 
246         END DO
247      END DO
248      tmask_h( 1 :iif,   :   ) = 0._wp      ! first columns
249      tmask_h(iil:jpi,   :   ) = 0._wp      ! last  columns (including mpp extra columns)
250      tmask_h(   :   , 1 :ijf) = 0._wp      ! first rows
251      tmask_h(   :   ,ijl:jpj) = 0._wp      ! last  rows (including mpp extra rows)
252      !
253      !                          ! north fold mask
254      tpol(1:jpiglo) = 1._wp 
255      fpol(1:jpiglo) = 1._wp
256      IF( jperio == 3 .OR. jperio == 4 ) THEN      ! T-point pivot
257         tpol(jpiglo/2+1:jpiglo) = 0._wp
258         fpol(     1    :jpiglo) = 0._wp
259         IF( mjg(nlej) == jpjglo ) THEN                  ! only half of the nlcj-1 row for tmask_h
260            DO ji = iif+1, iil-1
261               tmask_h(ji,nlej-1) = tmask_h(ji,nlej-1) * tpol(mig(ji))
262            END DO
263         ENDIF
264      ENDIF
265      !
266      IF( jperio == 5 .OR. jperio == 6 ) THEN      ! F-point pivot
267         tpol(     1    :jpiglo) = 0._wp
268         fpol(jpiglo/2+1:jpiglo) = 0._wp
269      ENDIF
270      !
271      !                          ! interior mask : 2D ocean mask x halo mask
272!$OMP PARALLEL DO schedule(static) private(jj, ji)
273      DO jj = 1, jpj
274         DO ji = 1, jpi
275            tmask_i(ji,jj) = ssmask(ji,jj) * tmask_h(ji,jj)
276         END DO
277      END DO
278
279
280      ! Lateral boundary conditions on velocity (modify fmask)
281      ! --------------------------------------- 
282      IF( rn_shlat /= 0 ) THEN      ! Not free-slip lateral boundary condition
283         !
284         CALL wrk_alloc( jpi,jpj,   zwf )
285         !
286!$OMP PARALLEL
287         DO jk = 1, jpk
288!$OMP DO schedule(static) private(jj, ji)
289            DO jj = 1, jpj
290               DO ji = 1, jpi
291                  zwf(ji,jj) = fmask(ji,jj,jk)         
292               END DO
293            END DO
294!$OMP DO schedule(static) private(jj, ji)
295            DO jj = 2, jpjm1
296               DO ji = fs_2, fs_jpim1   ! vector opt.
297                  IF( fmask(ji,jj,jk) == 0._wp ) THEN
298                     fmask(ji,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,jj), zwf(ji,jj+1),   &
299                        &                                           zwf(ji-1,jj), zwf(ji,jj-1)  )  )
300                  ENDIF
301               END DO
302            END DO
303!$OMP DO schedule(static) private(jj)
304            DO jj = 2, jpjm1
305               IF( fmask(1,jj,jk) == 0._wp ) THEN
306                  fmask(1  ,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1) ) )
307               ENDIF
308               IF( fmask(jpi,jj,jk) == 0._wp ) THEN
309                  fmask(jpi,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) )
310               ENDIF
311            END DO         
312!$OMP DO schedule(static) private(ji)
313            DO ji = 2, jpim1
314               IF( fmask(ji,1,jk) == 0._wp ) THEN
315                  fmask(ji, 1 ,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1) ) )
316               ENDIF
317               IF( fmask(ji,jpj,jk) == 0._wp ) THEN
318                  fmask(ji,jpj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpjm1) ) )
319               ENDIF
320            END DO
321         END DO
322!$OMP END PARALLEL
323         !
324         CALL wrk_dealloc( jpi,jpj,   zwf )
325         !
326         CALL lbc_lnk( fmask, 'F', 1._wp )      ! Lateral boundary conditions on fmask
327         !
328         ! CAUTION : The fmask may be further modified in dyn_vor_init ( dynvor.F90 ) depending on ln_vorlat
329         !
330      ENDIF
331     
332      ! User defined alteration of fmask (use to reduce ocean transport in specified straits)
333      ! --------------------------------
334      !
335      CALL usr_def_fmask( cn_cfg, nn_cfg, fmask )
336      !
337      !
338      IF( nn_timing == 1 )  CALL timing_stop('dom_msk')
339      !
340   END SUBROUTINE dom_msk
341   
342   !!======================================================================
343END MODULE dommsk
Note: See TracBrowser for help on using the repository browser.