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

source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/TRA/traadv_cen.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: 9.7 KB
Line 
1MODULE traadv_cen
2   !!======================================================================
3   !!                     ***  MODULE  traadv_cen  ***
4   !! Ocean  tracers:   advective trend (2nd/4th order centered)
5   !!======================================================================
6   !! History :  3.7  ! 2014-05  (G. Madec)  original code
7   !!----------------------------------------------------------------------
8
9   !!----------------------------------------------------------------------
10   !!   tra_adv_cen   : update the tracer trend with the advection trends using a centered or scheme (2nd or 4th order)
11   !!                   NB: on the vertical it is actually a 4th order COMPACT scheme which is used
12   !!----------------------------------------------------------------------
13   USE dom_oce        ! ocean space and time domain
14   USE eosbn2         ! equation of state
15   USE traadv_fct     ! acces to routine interp_4th_cpt
16   USE trd_oce        ! trends: ocean variables
17   USE trdtra         ! trends manager: tracers
18   USE diaptr         ! poleward transport diagnostics
19   USE diaar5         ! AR5 diagnostics
20   !
21   USE in_out_manager ! I/O manager
22   USE iom            ! IOM library
23   USE trc_oce        ! share passive tracers/Ocean variables
24   USE lib_mpp        ! MPP library
25
26   IMPLICIT NONE
27   PRIVATE
28
29   PUBLIC   tra_adv_cen   ! called by traadv.F90
30   
31   REAL(wp) ::   r1_6 = 1._wp / 6._wp   ! =1/6
32
33   LOGICAL ::   l_trd   ! flag to compute trends
34   LOGICAL ::   l_ptr   ! flag to compute poleward transport
35   LOGICAL ::   l_hst   ! flag to compute heat/salt transport
36
37   !! * Substitutions
38#  include "vectopt_loop_substitute.h90"
39#  include "do_loop_substitute.h90"
40   !!----------------------------------------------------------------------
41   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
42   !! $Id$
43   !! Software governed by the CeCILL license (see ./LICENSE)
44   !!----------------------------------------------------------------------
45CONTAINS
46
47   SUBROUTINE tra_adv_cen( kt, kit000, cdtype, pU, pV, pW,     &
48      &                    Kmm, pt, kjpt, Krhs, kn_cen_h, kn_cen_v ) 
49      !!----------------------------------------------------------------------
50      !!                  ***  ROUTINE tra_adv_cen  ***
51      !!                 
52      !! ** Purpose :   Compute the now trend due to the advection of tracers
53      !!      and add it to the general trend of passive tracer equations.
54      !!
55      !! ** Method  :   The advection is evaluated by a 2nd or 4th order scheme
56      !!               using now fields (leap-frog scheme).
57      !!       kn_cen_h = 2  ==>> 2nd order centered scheme on the horizontal
58      !!                = 4  ==>> 4th order    -        -       -      -
59      !!       kn_cen_v = 2  ==>> 2nd order centered scheme on the vertical
60      !!                = 4  ==>> 4th order COMPACT  scheme     -      -
61      !!
62      !! ** Action : - update pt(:,:,:,:,Krhs)  with the now advective tracer trends
63      !!             - send trends to trdtra module for further diagnostcs (l_trdtra=T)
64      !!             - poleward advective heat and salt transport (l_diaptr=T)
65      !!----------------------------------------------------------------------
66      INTEGER                                  , INTENT(in   ) ::   kt              ! ocean time-step index
67      INTEGER                                  , INTENT(in   ) ::   Kmm, Krhs       ! ocean time level indices
68      INTEGER                                  , INTENT(in   ) ::   kit000          ! first time step index
69      CHARACTER(len=3)                         , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator)
70      INTEGER                                  , INTENT(in   ) ::   kjpt            ! number of tracers
71      INTEGER                                  , INTENT(in   ) ::   kn_cen_h        ! =2/4 (2nd or 4th order scheme)
72      INTEGER                                  , INTENT(in   ) ::   kn_cen_v        ! =2/4 (2nd or 4th order scheme)
73      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pU, pV, pW      ! 3 ocean volume flux components
74      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! tracers and RHS of tracer equation
75      !
76      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
77      INTEGER  ::   ierr             ! local integer
78      REAL(wp) ::   zC2t_u, zC4t_u   ! local scalars
79      REAL(wp) ::   zC2t_v, zC4t_v   !   -      -
80      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwx, zwy, zwz, ztu, ztv, ztw
81      !!----------------------------------------------------------------------
82      !
83      IF( kt == kit000 )  THEN
84         IF(lwp) WRITE(numout,*)
85         IF(lwp) WRITE(numout,*) 'tra_adv_cen : centered advection scheme on ', cdtype, ' order h/v =', kn_cen_h,'/', kn_cen_v
86         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~ '
87      ENDIF
88      !                          ! set local switches
89      l_trd = .FALSE.
90      l_hst = .FALSE.
91      l_ptr = .FALSE.
92      IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )       l_trd = .TRUE.
93      IF(   cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) )  )    l_ptr = .TRUE. 
94      IF(   cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. &
95         &                          iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) )  l_hst = .TRUE.
96      !
97      !                   
98      zwz(:,:, 1 ) = 0._wp       ! surface & bottom vertical flux set to zero for all tracers
99      zwz(:,:,jpk) = 0._wp
100      !
101      DO jn = 1, kjpt            !==  loop over the tracers  ==!
102         !
103         SELECT CASE( kn_cen_h )       !--  Horizontal fluxes  --!
104         !
105         CASE(  2  )                         !* 2nd order centered
106            DO_3D_10_10( 1, jpkm1 )
107               zwx(ji,jj,jk) = 0.5_wp * pU(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj  ,jk,jn,Kmm) )
108               zwy(ji,jj,jk) = 0.5_wp * pV(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji  ,jj+1,jk,jn,Kmm) )
109            END_3D
110            !
111         CASE(  4  )                         !* 4th order centered
112            ztu(:,:,jpk) = 0._wp                   ! Bottom value : flux set to zero
113            ztv(:,:,jpk) = 0._wp
114            DO_3D_00_00( 1, jpkm1 )
115               ztu(ji,jj,jk) = ( pt(ji+1,jj  ,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * umask(ji,jj,jk)
116               ztv(ji,jj,jk) = ( pt(ji  ,jj+1,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * vmask(ji,jj,jk)
117            END_3D
118            CALL lbc_lnk_multi( 'traadv_cen', ztu, 'U', -1. , ztv, 'V', -1. )   ! Lateral boundary cond.
119            !
120            DO_3D_00_10( 1, jpkm1 )
121               zC2t_u = pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj  ,jk,jn,Kmm)   ! C2 interpolation of T at u- & v-points (x2)
122               zC2t_v = pt(ji,jj,jk,jn,Kmm) + pt(ji  ,jj+1,jk,jn,Kmm)
123               !                                                  ! C4 interpolation of T at u- & v-points (x2)
124               zC4t_u =  zC2t_u + r1_6 * ( ztu(ji-1,jj,jk) - ztu(ji+1,jj,jk) )
125               zC4t_v =  zC2t_v + r1_6 * ( ztv(ji,jj-1,jk) - ztv(ji,jj+1,jk) )
126               !                                                  ! C4 fluxes
127               zwx(ji,jj,jk) =  0.5_wp * pU(ji,jj,jk) * zC4t_u
128               zwy(ji,jj,jk) =  0.5_wp * pV(ji,jj,jk) * zC4t_v
129            END_3D
130            !
131         CASE DEFAULT
132            CALL ctl_stop( 'traadv_fct: wrong value for nn_fct' )
133         END SELECT
134         !
135         SELECT CASE( kn_cen_v )       !--  Vertical fluxes  --!   (interior)
136         !
137         CASE(  2  )                         !* 2nd order centered
138            DO_3D_00_00( 2, jpk )
139               zwz(ji,jj,jk) = 0.5 * pW(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji,jj,jk-1,jn,Kmm) ) * wmask(ji,jj,jk)
140            END_3D
141            !
142         CASE(  4  )                         !* 4th order compact
143            CALL interp_4th_cpt( pt(:,:,:,jn,Kmm) , ztw )      ! ztw = interpolated value of T at w-point
144            DO_3D_00_00( 2, jpkm1 )
145               zwz(ji,jj,jk) = pW(ji,jj,jk) * ztw(ji,jj,jk) * wmask(ji,jj,jk)
146            END_3D
147            !
148         END SELECT
149         !
150         IF( ln_linssh ) THEN                !* top value   (linear free surf. only as zwz is multiplied by wmask)
151            IF( ln_isfcav ) THEN                  ! ice-shelf cavities (top of the ocean)
152               DO_2D_11_11
153                  zwz(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kmm) 
154               END_2D
155            ELSE                                   ! no ice-shelf cavities (only ocean surface)
156               zwz(:,:,1) = pW(:,:,1) * pt(:,:,1,jn,Kmm)
157            ENDIF
158         ENDIF
159         !               
160         DO_3D_00_00( 1, jpkm1 )
161            pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs)    &
162               &             - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )    &
163               &                + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )    &
164               &                + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1)  ) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
165         END_3D
166         !                             ! trend diagnostics
167         IF( l_trd ) THEN
168            CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_xad, zwx, pU, pt(:,:,:,jn,Kmm) )
169            CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_yad, zwy, pV, pt(:,:,:,jn,Kmm) )
170            CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_zad, zwz, pW, pt(:,:,:,jn,Kmm) )
171         END IF
172         !                                 ! "Poleward" heat and salt transports
173         IF( l_ptr )   CALL dia_ptr_hst( jn, 'adv', zwy(:,:,:) )
174         !                                 !  heat and salt transport
175         IF( l_hst )   CALL dia_ar5_hst( jn, 'adv', zwx(:,:,:), zwy(:,:,:) )
176         !
177      END DO
178      !
179   END SUBROUTINE tra_adv_cen
180   
181   !!======================================================================
182END MODULE traadv_cen
Note: See TracBrowser for help on using the repository browser.