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_mus.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_mus.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.4 KB
Line 
1MODULE traadv_mus
2   !!======================================================================
3   !!                       ***  MODULE  traadv_mus  ***
4   !! Ocean  tracers:  horizontal & vertical advective trend
5   !!======================================================================
6   !! History :       !  2000-06  (A.Estublier)  for passive tracers
7   !!                 !  2001-08  (E.Durand, G.Madec)  adapted for T & S
8   !!   NEMO     1.0  !  2002-06  (G. Madec)  F90: Free form and module
9   !!            3.2  !  2010-05  (C. Ethe, G. Madec)  merge TRC-TRA + switch from velocity to transport
10   !!            3.4  !  2012-06  (P. Oddo, M. Vichi) include the upstream where needed
11   !!            3.7  !  2015-09  (G. Madec) add the ice-shelf cavities boundary condition
12   !!----------------------------------------------------------------------
13
14   !!----------------------------------------------------------------------
15   !!   tra_adv_mus   : update the tracer trend with the horizontal
16   !!                   and vertical advection trends using MUSCL scheme
17   !!----------------------------------------------------------------------
18   USE oce            ! ocean dynamics and active tracers
19   USE trc_oce        ! share passive tracers/Ocean variables
20   USE dom_oce        ! ocean space and time domain
21   USE trd_oce        ! trends: ocean variables
22   USE trdtra         ! tracers trends manager
23   USE sbcrnf         ! river runoffs
24   USE diaptr         ! poleward transport diagnostics
25   USE diaar5         ! AR5 diagnostics
26
27   !
28   USE iom            ! XIOS library
29   USE in_out_manager ! I/O manager
30   USE lib_mpp        ! distribued memory computing
31   USE lbclnk         ! ocean lateral boundary condition (or mpp link)
32   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
33
34   IMPLICIT NONE
35   PRIVATE
36
37   PUBLIC   tra_adv_mus   ! routine called by traadv.F90
38   
39   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   upsmsk   !: mixed upstream/centered scheme near some straits
40   !                                                           !  and in closed seas (orca 2 and 1 configurations)
41   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   xind     !: mixed upstream/centered index
42   
43   LOGICAL  ::   l_trd   ! flag to compute trends
44   LOGICAL  ::   l_ptr   ! flag to compute poleward transport
45   LOGICAL  ::   l_hst   ! flag to compute heat/salt transport
46
47   !! * Substitutions
48#  include "vectopt_loop_substitute.h90"
49#  include "do_loop_substitute.h90"
50   !!----------------------------------------------------------------------
51   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
52   !! $Id$
53   !! Software governed by the CeCILL license (see ./LICENSE)
54   !!----------------------------------------------------------------------
55CONTAINS
56
57   SUBROUTINE tra_adv_mus( kt, kit000, cdtype, p2dt, pU, pV, pW,             &
58      &                    Kbb, Kmm, pt, kjpt, Krhs, ld_msc_ups )
59      !!----------------------------------------------------------------------
60      !!                    ***  ROUTINE tra_adv_mus  ***
61      !!
62      !! ** Purpose :   Compute the now trend due to total advection of tracers
63      !!              using a MUSCL scheme (Monotone Upstream-centered Scheme for
64      !!              Conservation Laws) and add it to the general tracer trend.
65      !!
66      !! ** Method  : MUSCL scheme plus centered scheme at ocean boundaries
67      !!              ld_msc_ups=T :
68      !!
69      !! ** Action : - update pt(:,:,:,:,Krhs)  with the now advective tracer trends
70      !!             - send trends to trdtra module for further diagnostcs (l_trdtra=T)
71      !!             - poleward advective heat and salt transport (ln_diaptr=T)
72      !!
73      !! References : Estubier, A., and M. Levy, Notes Techn. Pole de Modelisation
74      !!              IPSL, Sept. 2000 (http://www.lodyc.jussieu.fr/opa)
75      !!----------------------------------------------------------------------
76      INTEGER                                  , INTENT(in   ) ::   kt              ! ocean time-step index
77      INTEGER                                  , INTENT(in   ) ::   Kbb, Kmm, Krhs  ! ocean time level indices
78      INTEGER                                  , INTENT(in   ) ::   kit000          ! first time step index
79      CHARACTER(len=3)                         , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator)
80      INTEGER                                  , INTENT(in   ) ::   kjpt            ! number of tracers
81      LOGICAL                                  , INTENT(in   ) ::   ld_msc_ups      ! use upstream scheme within muscl
82      REAL(wp)                                 , INTENT(in   ) ::   p2dt            ! tracer time-step
83      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pU, pV, pW      ! 3 ocean volume flux components
84      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! tracers and RHS of tracer equation
85      !
86      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
87      INTEGER  ::   ierr             ! local integer
88      REAL(wp) ::   zu, z0u, zzwx, zw , zalpha   ! local scalars
89      REAL(wp) ::   zv, z0v, zzwy, z0w           !   -      -
90      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwx, zslpx   ! 3D workspace
91      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwy, zslpy   ! -      -
92      !!----------------------------------------------------------------------
93      !
94      IF( kt == kit000 )  THEN
95         IF(lwp) WRITE(numout,*)
96         IF(lwp) WRITE(numout,*) 'tra_adv : MUSCL advection scheme on ', cdtype
97         IF(lwp) WRITE(numout,*) '        : mixed up-stream           ', ld_msc_ups
98         IF(lwp) WRITE(numout,*) '~~~~~~~'
99         IF(lwp) WRITE(numout,*)
100         !
101         ! Upstream / MUSCL scheme indicator
102         !
103         ALLOCATE( xind(jpi,jpj,jpk), STAT=ierr )
104         xind(:,:,:) = 1._wp              ! set equal to 1 where up-stream is not needed
105         !
106         IF( ld_msc_ups ) THEN            ! define the upstream indicator (if asked)
107            ALLOCATE( upsmsk(jpi,jpj), STAT=ierr )
108            upsmsk(:,:) = 0._wp                             ! not upstream by default
109            !
110            DO jk = 1, jpkm1
111               xind(:,:,jk) = 1._wp                              &                 ! =>1 where up-stream is not needed
112                  &         - MAX ( rnfmsk(:,:) * rnfmsk_z(jk),  &                 ! =>0 near runoff mouths (& closed sea outflows)
113                  &                 upsmsk(:,:)                ) * tmask(:,:,jk)   ! =>0 in some user defined area
114            END DO
115         ENDIF 
116         !
117      ENDIF 
118      !     
119      l_trd = .FALSE.
120      l_hst = .FALSE.
121      l_ptr = .FALSE.
122      IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )      l_trd = .TRUE.
123      IF(   cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) )  )   l_ptr = .TRUE. 
124      IF(   cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. &
125         &                          iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) ) l_hst = .TRUE.
126      !
127      DO jn = 1, kjpt            !==  loop over the tracers  ==!
128         !
129         !                          !* Horizontal advective fluxes
130         !
131         !                                !-- first guess of the slopes
132         zwx(:,:,jpk) = 0._wp                   ! bottom values
133         zwy(:,:,jpk) = 0._wp 
134         DO_3D_10_10( 1, jpkm1 )
135            zwx(ji,jj,jk) = umask(ji,jj,jk) * ( pt(ji+1,jj,jk,jn,Kbb) - pt(ji,jj,jk,jn,Kbb) )
136            zwy(ji,jj,jk) = vmask(ji,jj,jk) * ( pt(ji,jj+1,jk,jn,Kbb) - pt(ji,jj,jk,jn,Kbb) )
137         END_3D
138         ! lateral boundary conditions   (changed sign)
139         CALL lbc_lnk_multi( 'traadv_mus', zwx, 'U', -1. , zwy, 'V', -1. )
140         !                                !-- Slopes of tracer
141         zslpx(:,:,jpk) = 0._wp                 ! bottom values
142         zslpy(:,:,jpk) = 0._wp
143         DO_3D_01_01( 1, jpkm1 )
144            zslpx(ji,jj,jk) =                    ( zwx(ji,jj,jk) + zwx(ji-1,jj  ,jk) )   &
145               &            * ( 0.25 + SIGN( 0.25, zwx(ji,jj,jk) * zwx(ji-1,jj  ,jk) ) )
146            zslpy(ji,jj,jk) =                    ( zwy(ji,jj,jk) + zwy(ji  ,jj-1,jk) )   &
147               &            * ( 0.25 + SIGN( 0.25, zwy(ji,jj,jk) * zwy(ji  ,jj-1,jk) ) )
148         END_3D
149         !
150         DO_3D_01_01( 1, jpkm1 )
151            zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) ) * MIN(    ABS( zslpx(ji  ,jj,jk) ),   &
152               &                                                 2.*ABS( zwx  (ji-1,jj,jk) ),   &
153               &                                                 2.*ABS( zwx  (ji  ,jj,jk) ) )
154            zslpy(ji,jj,jk) = SIGN( 1., zslpy(ji,jj,jk) ) * MIN(    ABS( zslpy(ji,jj  ,jk) ),   &
155               &                                                 2.*ABS( zwy  (ji,jj-1,jk) ),   &
156               &                                                 2.*ABS( zwy  (ji,jj  ,jk) ) )
157         END_3D
158         !
159         DO_3D_00_00( 1, jpkm1 )
160            ! MUSCL fluxes
161            z0u = SIGN( 0.5, pU(ji,jj,jk) )
162            zalpha = 0.5 - z0u
163            zu  = z0u - 0.5 * pU(ji,jj,jk) * p2dt * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm)
164            zzwx = pt(ji+1,jj,jk,jn,Kbb) + xind(ji,jj,jk) * zu * zslpx(ji+1,jj,jk)
165            zzwy = pt(ji  ,jj,jk,jn,Kbb) + xind(ji,jj,jk) * zu * zslpx(ji  ,jj,jk)
166            zwx(ji,jj,jk) = pU(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy )
167            !
168            z0v = SIGN( 0.5, pV(ji,jj,jk) )
169            zalpha = 0.5 - z0v
170            zv  = z0v - 0.5 * pV(ji,jj,jk) * p2dt * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm)
171            zzwx = pt(ji,jj+1,jk,jn,Kbb) + xind(ji,jj,jk) * zv * zslpy(ji,jj+1,jk)
172            zzwy = pt(ji,jj  ,jk,jn,Kbb) + xind(ji,jj,jk) * zv * zslpy(ji,jj  ,jk)
173            zwy(ji,jj,jk) = pV(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy )
174         END_3D
175         CALL lbc_lnk_multi( 'traadv_mus', zwx, 'U', -1. , zwy, 'V', -1. )   ! lateral boundary conditions   (changed sign)
176         !
177         DO_3D_00_00( 1, jpkm1 )
178            pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )       &
179            &                                     + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  ) )     &
180            &                                   * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
181         END_3D
182         !                                ! trend diagnostics
183         IF( l_trd )  THEN
184            CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_xad, zwx, pU, pt(:,:,:,jn,Kbb) )
185            CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_yad, zwy, pV, pt(:,:,:,jn,Kbb) )
186         END IF
187         !                                 ! "Poleward" heat and salt transports
188         IF( l_ptr )  CALL dia_ptr_hst( jn, 'adv', zwy(:,:,:) )
189         !                                 !  heat transport
190         IF( l_hst )  CALL dia_ar5_hst( jn, 'adv', zwx(:,:,:), zwy(:,:,:) )
191         !
192         !                          !* Vertical advective fluxes
193         !
194         !                                !-- first guess of the slopes
195         zwx(:,:, 1 ) = 0._wp                   ! surface & bottom boundary conditions
196         zwx(:,:,jpk) = 0._wp
197         DO jk = 2, jpkm1                       ! interior values
198            zwx(:,:,jk) = tmask(:,:,jk) * ( pt(:,:,jk-1,jn,Kbb) - pt(:,:,jk,jn,Kbb) )
199         END DO
200         !                                !-- Slopes of tracer
201         zslpx(:,:,1) = 0._wp                   ! surface values
202         DO_3D_11_11( 2, jpkm1 )
203            zslpx(ji,jj,jk) =                     ( zwx(ji,jj,jk) + zwx(ji,jj,jk+1) )  &
204               &            * (  0.25 + SIGN( 0.25, zwx(ji,jj,jk) * zwx(ji,jj,jk+1) )  )
205         END_3D
206         DO_3D_11_11( 2, jpkm1 )
207            zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) ) * MIN(    ABS( zslpx(ji,jj,jk  ) ),   &
208               &                                                 2.*ABS( zwx  (ji,jj,jk+1) ),   &
209               &                                                 2.*ABS( zwx  (ji,jj,jk  ) )  )
210         END_3D
211         DO_3D_00_00( 1, jpk-2 )
212            z0w = SIGN( 0.5, pW(ji,jj,jk+1) )
213            zalpha = 0.5 + z0w
214            zw  = z0w - 0.5 * pW(ji,jj,jk+1) * p2dt * r1_e1e2t(ji,jj) / e3w(ji,jj,jk+1,Kmm)
215            zzwx = pt(ji,jj,jk+1,jn,Kbb) + xind(ji,jj,jk) * zw * zslpx(ji,jj,jk+1)
216            zzwy = pt(ji,jj,jk  ,jn,Kbb) + xind(ji,jj,jk) * zw * zslpx(ji,jj,jk  )
217            zwx(ji,jj,jk+1) = pW(ji,jj,jk+1) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) * wmask(ji,jj,jk)
218         END_3D
219         IF( ln_linssh ) THEN                   ! top values, linear free surface only
220            IF( ln_isfcav ) THEN                      ! ice-shelf cavities (top of the ocean)
221               DO_2D_11_11
222                  zwx(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kbb)
223               END_2D
224            ELSE                                      ! no cavities: only at the ocean surface
225               zwx(:,:,1) = pW(:,:,1) * pt(:,:,1,jn,Kbb)
226            ENDIF
227         ENDIF
228         !
229         DO_3D_00_00( 1, jpkm1 )
230            pt(ji,jj,jk,jn,Krhs) =  pt(ji,jj,jk,jn,Krhs) - ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) ) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
231         END_3D
232         !                                ! send trends for diagnostic
233         IF( l_trd )  CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_zad, zwx, pW, pt(:,:,:,jn,Kbb) )
234         !
235      END DO                     ! end of tracer loop
236      !
237   END SUBROUTINE tra_adv_mus
238
239   !!======================================================================
240END MODULE traadv_mus
Note: See TracBrowser for help on using the repository browser.