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.
dynldf_lap_blp.F90 in NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/DYN – NEMO

source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/DYN/dynldf_lap_blp.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: 7.9 KB
Line 
1MODULE dynldf_lap_blp
2   !!======================================================================
3   !!                   ***  MODULE  dynldf_lap_blp  ***
4   !! Ocean dynamics:  lateral viscosity trend (laplacian and bilaplacian)
5   !!======================================================================
6   !! History : 3.7  ! 2014-01  (G. Madec, S. Masson)  Original code, re-entrant laplacian
7   !!----------------------------------------------------------------------
8
9   !!----------------------------------------------------------------------
10   !!   dyn_ldf_lap   : update the momentum trend with the lateral viscosity using an iso-level   laplacian operator
11   !!   dyn_ldf_blp   : update the momentum trend with the lateral viscosity using an iso-level bilaplacian operator
12   !!----------------------------------------------------------------------
13   USE oce            ! ocean dynamics and tracers
14   USE dom_oce        ! ocean space and time domain
15   USE ldfdyn         ! lateral diffusion: eddy viscosity coef.
16   USE ldfslp         ! iso-neutral slopes
17   USE zdf_oce        ! ocean vertical physics
18   !
19   USE in_out_manager ! I/O manager
20   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
21
22   IMPLICIT NONE
23   PRIVATE
24
25   PUBLIC dyn_ldf_lap  ! called by dynldf.F90
26   PUBLIC dyn_ldf_blp  ! called by dynldf.F90
27
28   !! * Substitutions
29#  include "vectopt_loop_substitute.h90"
30#  include "do_loop_substitute.h90"
31   !!----------------------------------------------------------------------
32   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
33   !! $Id$
34   !! Software governed by the CeCILL license (see ./LICENSE)
35   !!----------------------------------------------------------------------
36CONTAINS
37
38   SUBROUTINE dyn_ldf_lap( kt, Kbb, Kmm, pu, pv, pu_rhs, pv_rhs, kpass )
39      !!----------------------------------------------------------------------
40      !!                     ***  ROUTINE dyn_ldf_lap  ***
41      !!                       
42      !! ** Purpose :   Compute the before horizontal momentum diffusive
43      !!      trend and add it to the general trend of momentum equation.
44      !!
45      !! ** Method  :   The Laplacian operator apply on horizontal velocity is
46      !!      writen as :   grad_h( ahmt div_h(U )) - curl_h( ahmf curl_z(U) )
47      !!
48      !! ** Action : - pu_rhs, pv_rhs increased by the harmonic operator applied on pu, pv.
49      !!----------------------------------------------------------------------
50      INTEGER                         , INTENT(in   ) ::   kt         ! ocean time-step index
51      INTEGER                         , INTENT(in   ) ::   Kbb, Kmm   ! ocean time level indices
52      INTEGER                         , INTENT(in   ) ::   kpass      ! =1/2 first or second passage
53      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pu, pv     ! before velocity  [m/s]
54      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_rhs   ! velocity trend   [m/s2]
55      !
56      INTEGER  ::   ji, jj, jk   ! dummy loop indices
57      REAL(wp) ::   zsign        ! local scalars
58      REAL(wp) ::   zua, zva     ! local scalars
59      REAL(wp), DIMENSION(jpi,jpj) ::   zcur, zdiv
60      !!----------------------------------------------------------------------
61      !
62      IF( kt == nit000 .AND. lwp ) THEN
63         WRITE(numout,*)
64         WRITE(numout,*) 'dyn_ldf : iso-level harmonic (laplacian) operator, pass=', kpass
65         WRITE(numout,*) '~~~~~~~ '
66      ENDIF
67      !
68      IF( kpass == 1 ) THEN   ;   zsign =  1._wp      ! bilaplacian operator require a minus sign
69      ELSE                    ;   zsign = -1._wp      !  (eddy viscosity coef. >0)
70      ENDIF
71      !
72      !                                                ! ===============
73      DO jk = 1, jpkm1                                 ! Horizontal slab
74         !                                             ! ===============
75         DO_2D_01_01
76            !                                      ! ahm * e3 * curl  (computed from 1 to jpim1/jpjm1)
77!!gm open question here : e3f  at before or now ?    probably now...
78!!gm note that ahmf has already been multiplied by fmask
79            zcur(ji-1,jj-1) = ahmf(ji-1,jj-1,jk) * e3f(ji-1,jj-1,jk) * r1_e1e2f(ji-1,jj-1)       &
80               &     * (  e2v(ji  ,jj-1) * pv(ji  ,jj-1,jk) - e2v(ji-1,jj-1) * pv(ji-1,jj-1,jk)  &
81               &        - e1u(ji-1,jj  ) * pu(ji-1,jj  ,jk) + e1u(ji-1,jj-1) * pu(ji-1,jj-1,jk)  )
82            !                                      ! ahm * div        (computed from 2 to jpi/jpj)
83!!gm note that ahmt has already been multiplied by tmask
84            zdiv(ji,jj)     = ahmt(ji,jj,jk) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kbb)                                         &
85               &     * (  e2u(ji,jj)*e3u(ji,jj,jk,Kbb) * pu(ji,jj,jk) - e2u(ji-1,jj)*e3u(ji-1,jj,jk,Kbb) * pu(ji-1,jj,jk)  &
86               &        + e1v(ji,jj)*e3v(ji,jj,jk,Kbb) * pv(ji,jj,jk) - e1v(ji,jj-1)*e3v(ji,jj-1,jk,Kbb) * pv(ji,jj-1,jk)  )
87         END_2D
88         !
89         DO_2D_00_00
90            pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zsign * (                                                 &
91               &              - ( zcur(ji  ,jj) - zcur(ji,jj-1) ) * r1_e2u(ji,jj) / e3u(ji,jj,jk,Kmm)   &
92               &              + ( zdiv(ji+1,jj) - zdiv(ji,jj  ) ) * r1_e1u(ji,jj)                     )
93               !
94            pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zsign * (                                                 &
95               &                ( zcur(ji,jj  ) - zcur(ji-1,jj) ) * r1_e1v(ji,jj) / e3v(ji,jj,jk,Kmm)   &
96               &              + ( zdiv(ji,jj+1) - zdiv(ji  ,jj) ) * r1_e2v(ji,jj)                     )
97         END_2D
98         !                                             ! ===============
99      END DO                                           !   End of slab
100      !                                                ! ===============
101      !
102   END SUBROUTINE dyn_ldf_lap
103
104
105   SUBROUTINE dyn_ldf_blp( kt, Kbb, Kmm, pu, pv, pu_rhs, pv_rhs )
106      !!----------------------------------------------------------------------
107      !!                 ***  ROUTINE dyn_ldf_blp  ***
108      !!                   
109      !! ** Purpose :   Compute the before lateral momentum viscous trend
110      !!              and add it to the general trend of momentum equation.
111      !!
112      !! ** Method  :   The lateral viscous trends is provided by a bilaplacian
113      !!      operator applied to before field (forward in time).
114      !!      It is computed by two successive calls to dyn_ldf_lap routine
115      !!
116      !! ** Action :   pt(:,:,:,:,Krhs)   updated with the before rotated bilaplacian diffusion
117      !!----------------------------------------------------------------------
118      INTEGER                         , INTENT(in   ) ::   kt         ! ocean time-step index
119      INTEGER                         , INTENT(in   ) ::   Kbb, Kmm   ! ocean time level indices
120      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pu, pv     ! before velocity fields
121      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_rhs   ! momentum trend
122      !
123      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zulap, zvlap   ! laplacian at u- and v-point
124      !!----------------------------------------------------------------------
125      !
126      IF( kt == nit000 )  THEN
127         IF(lwp) WRITE(numout,*)
128         IF(lwp) WRITE(numout,*) 'dyn_ldf_blp : bilaplacian operator momentum '
129         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
130      ENDIF
131      !
132      zulap(:,:,:) = 0._wp
133      zvlap(:,:,:) = 0._wp
134      !
135      CALL dyn_ldf_lap( kt, Kbb, Kmm, pu, pv, zulap, zvlap, 1 )   ! rotated laplacian applied to pt (output in zlap,Kbb)
136      !
137      CALL lbc_lnk_multi( 'dynldf_lap_blp', zulap, 'U', -1., zvlap, 'V', -1. )             ! Lateral boundary conditions
138      !
139      CALL dyn_ldf_lap( kt, Kbb, Kmm, zulap, zvlap, pu_rhs, pv_rhs, 2 )   ! rotated laplacian applied to zlap (output in pt(:,:,:,:,Krhs))
140      !
141   END SUBROUTINE dyn_ldf_blp
142
143   !!======================================================================
144END MODULE dynldf_lap_blp
Note: See TracBrowser for help on using the repository browser.