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/2019/dev_r11943_MERGE_2019/src/OCE/DOM – NEMO

source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/DOM/dommsk.F90 @ 12340

Last change on this file since 12340 was 12340, checked in by acc, 4 years ago

Branch 2019/dev_r11943_MERGE_2019. This commit introduces basic do loop macro
substitution to the 2019 option 1, merge branch. These changes have been SETTE
tested. The only addition is the do_loop_substitute.h90 file in the OCE directory but
the macros defined therein are used throughout the code to replace identifiable, 2D-
and 3D- nested loop opening and closing statements with single-line alternatives. Code
indents are also adjusted accordingly.

The following explanation is taken from comments in the new header file:

This header file contains preprocessor definitions and macros used in the do-loop
substitutions introduced between version 4.0 and 4.2. The primary aim of these macros
is to assist in future applications of tiling to improve performance. This is expected
to be achieved by alternative versions of these macros in selected locations. The
initial introduction of these macros simply replaces all identifiable nested 2D- and
3D-loops with single line statements (and adjusts indenting accordingly). Do loops
are identifiable if they comform to either:

DO jk = ....

DO jj = .... DO jj = ...

DO ji = .... DO ji = ...
. OR .
. .

END DO END DO

END DO END DO

END DO

and white-space variants thereof.

Additionally, only loops with recognised jj and ji loops limits are treated; these are:
Lower limits of 1, 2 or fs_2
Upper limits of jpi, jpim1 or fs_jpim1 (for ji) or jpj, jpjm1 or fs_jpjm1 (for jj)

The macro naming convention takes the form: DO_2D_BT_LR where:

B is the Bottom offset from the PE's inner domain;
T is the Top offset from the PE's inner domain;
L is the Left offset from the PE's inner domain;
R is the Right offset from the PE's inner domain

So, given an inner domain of 2,jpim1 and 2,jpjm1, a typical example would replace:

DO jj = 2, jpj

DO ji = 1, jpim1
.
.

END DO

END DO

with:

DO_2D_01_10
.
.
END_2D

similar conventions apply to the 3D loops macros. jk loop limits are retained
through macro arguments and are not restricted. This includes the possibility of
strides for which an extra set of DO_3DS macros are defined.

In the example definition below the inner PE domain is defined by start indices of
(kIs, kJs) and end indices of (kIe, KJe)

#define DO_2D_00_00 DO jj = kJs, kJe ; DO ji = kIs, kIe
#define END_2D END DO ; END DO

TO DO:

Only conventional nested loops have been identified and replaced by this step. There are constructs such as:

DO jk = 2, jpkm1

z2d(:,:) = z2d(:,:) + e3w(:,:,jk,Kmm) * z3d(:,:,jk) * wmask(:,:,jk)

END DO

which may need to be considered.
  • Property svn:keywords set to Id
File size: 13.9 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        ! open boundary
29   !
30   USE in_out_manager ! I/O manager
31   USE iom            ! IOM library
32   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
33   USE lib_mpp        ! Massively Parallel Processing library
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#  include "do_loop_substitute.h90"
48   !!----------------------------------------------------------------------
49   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
50   !! $Id$
51   !! Software governed by the CeCILL license (see ./LICENSE)
52   !!----------------------------------------------------------------------
53CONTAINS
54
55   SUBROUTINE dom_msk( k_top, k_bot )
56      !!---------------------------------------------------------------------
57      !!                 ***  ROUTINE dom_msk  ***
58      !!
59      !! ** Purpose :   Compute land/ocean mask arrays at tracer points, hori-
60      !!      zontal velocity points (u & v), vorticity points (f) points.
61      !!
62      !! ** Method  :   The ocean/land mask  at t-point is deduced from ko_top
63      !!      and ko_bot, the indices of the fist and last ocean t-levels which
64      !!      are either defined in usrdef_zgr or read in zgr_read.
65      !!                The velocity masks (umask, vmask, wmask, wumask, wvmask)
66      !!      are deduced from a product of the two neighboring tmask.
67      !!                The vorticity mask (fmask) is deduced from tmask taking
68      !!      into account the choice of lateral boundary condition (rn_shlat) :
69      !!         rn_shlat = 0, free slip  (no shear along the coast)
70      !!         rn_shlat = 2, no slip  (specified zero velocity at the coast)
71      !!         0 < rn_shlat < 2, partial slip   | non-linear velocity profile
72      !!         2 < rn_shlat, strong slip        | in the lateral boundary layer
73      !!
74      !!      tmask_i : interior ocean mask at t-point, i.e. excluding duplicated
75      !!                rows/lines due to cyclic or North Fold boundaries as well
76      !!                as MPP halos.
77      !!      tmask_h : halo mask at t-point, i.e. excluding duplicated rows/lines
78      !!                due to cyclic or North Fold boundaries as well as MPP halos.
79      !!
80      !! ** Action :   tmask, umask, vmask, wmask, wumask, wvmask : land/ocean mask
81      !!                         at t-, u-, v- w, wu-, and wv-points (=0. or 1.)
82      !!               fmask   : land/ocean mask at f-point (=0., or =1., or
83      !!                         =rn_shlat along lateral boundaries)
84      !!               tmask_i : interior ocean mask
85      !!               tmask_h : halo mask
86      !!               ssmask , ssumask, ssvmask, ssfmask : 2D ocean mask
87      !!----------------------------------------------------------------------
88      INTEGER, DIMENSION(:,:), INTENT(in) ::   k_top, k_bot   ! first and last ocean level
89      !
90      INTEGER  ::   ji, jj, jk     ! dummy loop indices
91      INTEGER  ::   iif, iil       ! local integers
92      INTEGER  ::   ijf, ijl       !   -       -
93      INTEGER  ::   iktop, ikbot   !   -       -
94      INTEGER  ::   ios, inum
95      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zwf   ! 2D workspace
96      !!
97      NAMELIST/namlbc/ rn_shlat, ln_vorlat
98      NAMELIST/nambdy/ ln_bdy ,nb_bdy, ln_coords_file, cn_coords_file,         &
99         &             ln_mask_file, cn_mask_file, cn_dyn2d, nn_dyn2d_dta,     &
100         &             cn_dyn3d, nn_dyn3d_dta, cn_tra, nn_tra_dta,             &
101         &             ln_tra_dmp, ln_dyn3d_dmp, rn_time_dmp, rn_time_dmp_out, &
102         &             cn_ice, nn_ice_dta,                                     &
103         &             ln_vol, nn_volctl, nn_rimwidth
104      !!---------------------------------------------------------------------
105      !
106      READ  ( numnam_ref, namlbc, IOSTAT = ios, ERR = 901 )
107901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namlbc in reference namelist' )
108      READ  ( numnam_cfg, namlbc, IOSTAT = ios, ERR = 902 )
109902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namlbc in configuration namelist' )
110      IF(lwm) WRITE ( numond, namlbc )
111     
112      IF(lwp) THEN                  ! control print
113         WRITE(numout,*)
114         WRITE(numout,*) 'dommsk : ocean mask '
115         WRITE(numout,*) '~~~~~~'
116         WRITE(numout,*) '   Namelist namlbc'
117         WRITE(numout,*) '      lateral momentum boundary cond.    rn_shlat  = ',rn_shlat
118         WRITE(numout,*) '      consistency with analytical form   ln_vorlat = ',ln_vorlat 
119      ENDIF
120      !
121      IF(lwp) WRITE(numout,*)
122      IF     (      rn_shlat == 0.               ) THEN   ;   IF(lwp) WRITE(numout,*) '   ==>>>   ocean lateral  free-slip'
123      ELSEIF (      rn_shlat == 2.               ) THEN   ;   IF(lwp) WRITE(numout,*) '   ==>>>   ocean lateral  no-slip'
124      ELSEIF ( 0. < rn_shlat .AND. rn_shlat < 2. ) THEN   ;   IF(lwp) WRITE(numout,*) '   ==>>>   ocean lateral  partial-slip'
125      ELSEIF ( 2. < rn_shlat                     ) THEN   ;   IF(lwp) WRITE(numout,*) '   ==>>>   ocean lateral  strong-slip'
126      ELSE
127         CALL ctl_stop( 'dom_msk: wrong value for rn_shlat (i.e. a negalive value). We stop.' )
128      ENDIF
129
130      !  Ocean/land mask at t-point  (computed from ko_top and ko_bot)
131      ! ----------------------------
132      !
133      tmask(:,:,:) = 0._wp
134      DO_2D_11_11
135         iktop = k_top(ji,jj)
136         ikbot = k_bot(ji,jj)
137         IF( iktop /= 0 ) THEN       ! water in the column
138            tmask(ji,jj,iktop:ikbot  ) = 1._wp
139         ENDIF
140      END_2D
141      !
142      ! the following call is mandatory
143      ! it masks boundaries (bathy=0) where needed depending on the configuration (closed, periodic...) 
144      CALL lbc_lnk( 'dommsk', tmask  , 'T', 1._wp )      ! Lateral boundary conditions
145
146     ! Mask corrections for bdy (read in mppini2)
147      READ  ( numnam_ref, nambdy, IOSTAT = ios, ERR = 903)
148903   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambdy in reference namelist' )
149      READ  ( numnam_cfg, nambdy, IOSTAT = ios, ERR = 904 )
150904   IF( ios >  0 )   CALL ctl_nam ( ios , 'nambdy in configuration namelist' )
151      ! ------------------------
152      IF ( ln_bdy .AND. ln_mask_file ) THEN
153         CALL iom_open( cn_mask_file, inum )
154         CALL iom_get ( inum, jpdom_data, 'bdy_msk', bdytmask(:,:) )
155         CALL iom_close( inum )
156         DO_3D_11_11( 1, jpkm1 )
157            tmask(ji,jj,jk) = tmask(ji,jj,jk) * bdytmask(ji,jj)
158         END_3D
159      ENDIF
160         
161      !  Ocean/land mask at u-, v-, and f-points   (computed from tmask)
162      ! ----------------------------------------
163      ! NB: at this point, fmask is designed for free slip lateral boundary condition
164      DO jk = 1, jpk
165         DO jj = 1, jpjm1
166            DO ji = 1, fs_jpim1   ! vector loop
167               umask(ji,jj,jk) = tmask(ji,jj  ,jk) * tmask(ji+1,jj  ,jk)
168               vmask(ji,jj,jk) = tmask(ji,jj  ,jk) * tmask(ji  ,jj+1,jk)
169            END DO
170            DO ji = 1, jpim1      ! NO vector opt.
171               fmask(ji,jj,jk) = tmask(ji,jj  ,jk) * tmask(ji+1,jj  ,jk)   &
172                  &            * tmask(ji,jj+1,jk) * tmask(ji+1,jj+1,jk)
173            END DO
174         END DO
175      END DO
176      CALL lbc_lnk_multi( 'dommsk', umask, 'U', 1., vmask, 'V', 1., fmask, 'F', 1. )      ! Lateral boundary conditions
177 
178      ! Ocean/land mask at wu-, wv- and w points    (computed from tmask)
179      !-----------------------------------------
180      wmask (:,:,1) = tmask(:,:,1)     ! surface
181      wumask(:,:,1) = umask(:,:,1)
182      wvmask(:,:,1) = vmask(:,:,1)
183      DO jk = 2, jpk                   ! interior values
184         wmask (:,:,jk) = tmask(:,:,jk) * tmask(:,:,jk-1)
185         wumask(:,:,jk) = umask(:,:,jk) * umask(:,:,jk-1)   
186         wvmask(:,:,jk) = vmask(:,:,jk) * vmask(:,:,jk-1)
187      END DO
188
189
190      ! Ocean/land column mask at t-, u-, and v-points   (i.e. at least 1 wet cell in the vertical)
191      ! ----------------------------------------------
192      ssmask (:,:) = MAXVAL( tmask(:,:,:), DIM=3 )
193      ssumask(:,:) = MAXVAL( umask(:,:,:), DIM=3 )
194      ssvmask(:,:) = MAXVAL( vmask(:,:,:), DIM=3 )
195
196
197      ! Interior domain mask  (used for global sum)
198      ! --------------------
199      !
200      iif = nn_hls   ;   iil = nlci - nn_hls + 1
201      ijf = nn_hls   ;   ijl = nlcj - nn_hls + 1
202      !
203      !                          ! halo mask : 0 on the halo and 1 elsewhere
204      tmask_h(:,:) = 1._wp                 
205      tmask_h( 1 :iif,   :   ) = 0._wp      ! first columns
206      tmask_h(iil:jpi,   :   ) = 0._wp      ! last  columns (including mpp extra columns)
207      tmask_h(   :   , 1 :ijf) = 0._wp      ! first rows
208      tmask_h(   :   ,ijl:jpj) = 0._wp      ! last  rows (including mpp extra rows)
209      !
210      !                          ! north fold mask
211      tpol(1:jpiglo) = 1._wp 
212      fpol(1:jpiglo) = 1._wp
213      IF( jperio == 3 .OR. jperio == 4 ) THEN      ! T-point pivot
214         tpol(jpiglo/2+1:jpiglo) = 0._wp
215         fpol(     1    :jpiglo) = 0._wp
216         IF( mjg(nlej) == jpjglo ) THEN                  ! only half of the nlcj-1 row for tmask_h
217            DO ji = iif+1, iil-1
218               tmask_h(ji,nlej-1) = tmask_h(ji,nlej-1) * tpol(mig(ji))
219            END DO
220         ENDIF
221      ENDIF
222      !
223      IF( jperio == 5 .OR. jperio == 6 ) THEN      ! F-point pivot
224         tpol(     1    :jpiglo) = 0._wp
225         fpol(jpiglo/2+1:jpiglo) = 0._wp
226      ENDIF
227      !
228      !                          ! interior mask : 2D ocean mask x halo mask
229      tmask_i(:,:) = ssmask(:,:) * tmask_h(:,:)
230
231
232      ! Lateral boundary conditions on velocity (modify fmask)
233      ! --------------------------------------- 
234      IF( rn_shlat /= 0 ) THEN      ! Not free-slip lateral boundary condition
235         !
236         ALLOCATE( zwf(jpi,jpj) )
237         !
238         DO jk = 1, jpk
239            zwf(:,:) = fmask(:,:,jk)         
240            DO_2D_00_00
241               IF( fmask(ji,jj,jk) == 0._wp ) THEN
242                  fmask(ji,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,jj), zwf(ji,jj+1),   &
243                     &                                           zwf(ji-1,jj), zwf(ji,jj-1)  )  )
244               ENDIF
245            END_2D
246            DO jj = 2, jpjm1
247               IF( fmask(1,jj,jk) == 0._wp ) THEN
248                  fmask(1  ,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1) ) )
249               ENDIF
250               IF( fmask(jpi,jj,jk) == 0._wp ) THEN
251                  fmask(jpi,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) )
252               ENDIF
253            END DO         
254            DO ji = 2, jpim1
255               IF( fmask(ji,1,jk) == 0._wp ) THEN
256                  fmask(ji, 1 ,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1) ) )
257               ENDIF
258               IF( fmask(ji,jpj,jk) == 0._wp ) THEN
259                  fmask(ji,jpj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpjm1) ) )
260               ENDIF
261            END DO
262#if defined key_agrif 
263            IF( .NOT. AGRIF_Root() ) THEN
264               IF ((nbondi ==  1).OR.(nbondi == 2)) fmask(nlci-1 , :     ,jk) = 0.e0      ! east
265               IF ((nbondi == -1).OR.(nbondi == 2)) fmask(1      , :     ,jk) = 0.e0      ! west
266               IF ((nbondj ==  1).OR.(nbondj == 2)) fmask(:      ,nlcj-1 ,jk) = 0.e0      ! north
267               IF ((nbondj == -1).OR.(nbondj == 2)) fmask(:      ,1      ,jk) = 0.e0      ! south
268            ENDIF 
269#endif
270         END DO
271         !
272         DEALLOCATE( zwf )
273         !
274         CALL lbc_lnk( 'dommsk', fmask, 'F', 1._wp )      ! Lateral boundary conditions on fmask
275         !
276         ! CAUTION : The fmask may be further modified in dyn_vor_init ( dynvor.F90 ) depending on ln_vorlat
277         !
278      ENDIF
279     
280      ! User defined alteration of fmask (use to reduce ocean transport in specified straits)
281      ! --------------------------------
282      !
283      CALL usr_def_fmask( cn_cfg, nn_cfg, fmask )
284      !
285   END SUBROUTINE dom_msk
286   
287   !!======================================================================
288END MODULE dommsk
Note: See TracBrowser for help on using the repository browser.