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.
trcatf.F90 in NEMO/branches/2019/dev_r11943_MERGE_2019/src/TOP/TRP – NEMO

source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/TOP/TRP/trcatf.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.0 KB
Line 
1MODULE trcatf
2   !!======================================================================
3   !!                       ***  MODULE  trcatf  ***
4   !! Ocean passive tracers:  time stepping on passives tracers
5   !!======================================================================
6   !! History :  7.0  !  1991-11  (G. Madec)  Original code
7   !!                 !  1993-03  (M. Guyon)  symetrical conditions
8   !!                 !  1995-02  (M. Levy)   passive tracers
9   !!                 !  1996-02  (G. Madec & M. Imbard)  opa release 8.0
10   !!            8.0  !  1996-04  (A. Weaver)  Euler forward step
11   !!            8.2  !  1999-02  (G. Madec, N. Grima)  semi-implicit pressure grad.
12   !!  NEMO      1.0  !  2002-08  (G. Madec)  F90: Free form and module
13   !!                 !  2002-08  (G. Madec)  F90: Free form and module
14   !!                 !  2002-11  (C. Talandier, A-M Treguier) Open boundaries
15   !!                 !  2004-03  (C. Ethe) passive tracers
16   !!                 !  2007-02  (C. Deltel) Diagnose ML trends for passive tracers
17   !!            2.0  !  2006-02  (L. Debreu, C. Mazauric) Agrif implementation
18   !!            3.0  !  2008-06  (G. Madec)  time stepping always done in trazdf
19   !!            3.1  !  2009-02  (G. Madec, R. Benshila)  re-introduce the vvl option
20   !!            3.3  !  2010-06  (C. Ethe, G. Madec) Merge TRA-TRC
21   !!            4.1  !  2019-08  (A. Coward, D. Storkey) rename trcnxt.F90 -> trcatf.F90. Now only does time filtering.
22   !!----------------------------------------------------------------------
23#if defined key_top
24   !!----------------------------------------------------------------------
25   !!   'key_top'                                                TOP models
26   !!----------------------------------------------------------------------
27   !!   trc_atf     : time stepping on passive tracers
28   !!----------------------------------------------------------------------
29   USE oce_trc         ! ocean dynamics and tracers variables
30   USE trc             ! ocean passive tracers variables
31   USE trd_oce
32   USE trdtra
33   USE traatf
34   USE bdy_oce   , ONLY: ln_bdy
35   USE trcbdy          ! BDY open boundaries
36# if defined key_agrif
37   USE agrif_top_interp
38# endif
39   !
40   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
41   USE prtctl_trc      ! Print control for debbuging
42
43   IMPLICIT NONE
44   PRIVATE
45
46   PUBLIC   trc_atf   ! routine called by step.F90
47
48   REAL(wp) ::   rfact1, rfact2
49
50   !! * Substitutions
51#  include "do_loop_substitute.h90"
52   !!----------------------------------------------------------------------
53   !! NEMO/TOP 4.0 , NEMO Consortium (2018)
54   !! $Id$
55   !! Software governed by the CeCILL license (see ./LICENSE)
56   !!----------------------------------------------------------------------
57CONTAINS
58
59   SUBROUTINE trc_atf( kt, Kbb, Kmm, Kaa, ptr )
60      !!----------------------------------------------------------------------
61      !!                   ***  ROUTINE trcatf  ***
62      !!
63      !! ** Purpose :   Apply the boundary condition on the after passive tracers fields and
64      !!      apply Asselin time filter to the now passive tracer fields if using leapfrog timestep
65      !!
66      !! ** Method  :   Apply lateral boundary conditions on (uu(Kaa),vv(Kaa)) through
67      !!      call to lbc_lnk routine
68      !!
69      !!   For Arakawa or TVD Scheme :
70      !!      A Asselin time filter applied on now tracers tr(Kmm) to avoid
71      !!      the divergence of two consecutive time-steps and tr arrays
72      !!      to prepare the next time_step:
73      !!         (tr(Kmm)) = (tr(Kmm)) + atfp [ (tr(Kbb)) + (tr(Kaa)) - 2 (tr(Kmm)) ]
74      !!
75      !!
76      !! ** Action  : - update tr(Kmm), tr(Kaa)
77      !!----------------------------------------------------------------------
78      INTEGER                                   , INTENT( in )  :: kt             ! ocean time-step index
79      INTEGER                                   , INTENT( in )  :: Kbb, Kmm, Kaa ! time level indices
80      REAL(wp), DIMENSION(jpi,jpj,jpk,jptra,jpt), INTENT(inout) :: ptr            ! passive tracers
81      !
82      INTEGER  ::   jk, jn   ! dummy loop indices
83      REAL(wp) ::   zfact            ! temporary scalar
84      CHARACTER (len=22) :: charout
85      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::   ztrdt    ! 4D workspace
86      !!----------------------------------------------------------------------
87      !
88      IF( ln_timing )   CALL timing_start('trc_atf')
89      !
90      IF( kt == nittrc000 .AND. lwp ) THEN
91         WRITE(numout,*)
92         WRITE(numout,*) 'trc_atf : Asselin time filtering on passive tracers'
93      ENDIF
94      !
95#if defined key_agrif
96      CALL Agrif_trc                   ! AGRIF zoom boundaries
97#endif
98      ! Update after tracer on domain lateral boundaries
99      CALL lbc_lnk( 'trcatf', ptr(:,:,:,:,Kaa), 'T', 1. )   
100
101      IF( ln_bdy )  CALL trc_bdy( kt, Kbb, Kmm, Kaa )
102
103      IF( l_trdtrc )  THEN             ! trends: store now fields before the Asselin filter application
104         ALLOCATE( ztrdt(jpi,jpj,jpk,jptra) )
105         ztrdt(:,:,:,:)  = 0._wp
106         IF( ln_traldf_iso ) THEN                       ! diagnose the "pure" Kz diffusive trend
107            DO jn = 1, jptra
108               CALL trd_tra( kt, Kmm, Kaa, 'TRC', jn, jptra_zdfp, ztrdt(:,:,:,jn) )
109            ENDDO
110         ENDIF
111
112         ! total trend for the non-time-filtered variables.
113         zfact = 1.0 / rdttrc
114         ! G Nurser 23 Mar 2017. Recalculate trend as Delta(e3t*T)/e3tn; e3tn cancel from ts(Kmm) terms
115         IF( ln_linssh ) THEN       ! linear sea surface height only
116            DO jn = 1, jptra
117               DO jk = 1, jpkm1
118                  ztrdt(:,:,jk,jn) = ( ptr(:,:,jk,jn,Kaa)*e3t(:,:,jk,Kaa) / e3t(:,:,jk,Kmm) - ptr(:,:,jk,jn,Kmm)) * zfact
119               END DO
120            END DO
121         ELSE
122            DO jn = 1, jptra
123               DO jk = 1, jpkm1
124                  ztrdt(:,:,jk,jn) = ( ptr(:,:,jk,jn,Kaa) - ptr(:,:,jk,jn,Kmm) ) * zfact
125               END DO
126            END DO
127         ENDIF
128         !
129         DO jn = 1, jptra
130            CALL trd_tra( kt, Kmm, Kaa, 'TRC', jn, jptra_tot, ztrdt(:,:,:,jn) )
131         ENDDO
132         !
133         IF( ln_linssh ) THEN       ! linear sea surface height only
134            ! Store now fields before applying the Asselin filter
135            ! in order to calculate Asselin filter trend later.
136            ztrdt(:,:,:,:) = ptr(:,:,:,:,Kmm) 
137         ENDIF
138
139      ENDIF
140      !                                ! Leap-Frog + Asselin filter time stepping
141      IF( (neuler == 0 .AND. kt == nittrc000) .OR. ln_top_euler ) THEN    ! Euler time-stepping
142         !
143         IF (l_trdtrc .AND. .NOT. ln_linssh ) THEN   ! Zero Asselin filter contribution must be explicitly written out since for vvl
144            !                                        ! Asselin filter is output by tra_nxt_vvl that is not called on this time step
145            ztrdt(:,:,:,:) = 0._wp           
146            DO jn = 1, jptra
147               CALL trd_tra( kt, Kmm, Kaa, 'TRC', jn, jptra_atf, ztrdt(:,:,:,jn) )
148            ENDDO
149         END IF
150         !
151      ELSE     
152         IF( .NOT. l_offline ) THEN ! Leap-Frog + Asselin filter time stepping
153            IF( ln_linssh ) THEN   ;   CALL tra_atf_fix( kt, Kbb, Kmm, Kaa, nittrc000,         'TRC', ptr, jptra )                     !     linear ssh
154            ELSE                   ;   CALL tra_atf_vvl( kt, Kbb, Kmm, Kaa, nittrc000, rdttrc, 'TRC', ptr, sbc_trc, sbc_trc_b, jptra ) ! non-linear ssh
155            ENDIF
156         ELSE
157                                       CALL trc_atf_off( kt, Kbb, Kmm, Kaa, ptr )       ! offline
158         ENDIF
159         !
160         CALL lbc_lnk_multi( 'trcatf', ptr(:,:,:,:,Kmm), 'T', 1._wp, ptr(:,:,:,:,Kaa), 'T', 1._wp, ptr(:,:,:,:,Kaa), 'T', 1._wp )
161      ENDIF
162      !
163      IF( l_trdtrc .AND. ln_linssh ) THEN      ! trend of the Asselin filter (tb filtered - tb)/dt )
164         DO jn = 1, jptra
165            DO jk = 1, jpkm1
166               zfact = 1._wp / r2dttrc 
167               ztrdt(:,:,jk,jn) = ( ptr(:,:,jk,jn,Kbb) - ztrdt(:,:,jk,jn) ) * zfact 
168            END DO
169            CALL trd_tra( kt, Kmm, Kaa, 'TRC', jn, jptra_atf, ztrdt(:,:,:,jn) )
170         END DO
171      END IF
172      IF( l_trdtrc ) DEALLOCATE( ztrdt ) 
173      !
174      IF(sn_cfctl%l_prttrc)   THEN  ! print mean trends (used for debugging)
175         WRITE(charout, FMT="('nxt')")
176         CALL prt_ctl_trc_info(charout)
177         CALL prt_ctl_trc(tab4d=ptr(:,:,:,:,Kmm), mask=tmask, clinfo=ctrcnm)
178      ENDIF
179      !
180      IF( ln_timing )   CALL timing_stop('trc_atf')
181      !
182   END SUBROUTINE trc_atf
183
184
185   SUBROUTINE trc_atf_off( kt, Kbb, Kmm, Kaa, ptr )
186      !!----------------------------------------------------------------------
187      !!                   ***  ROUTINE tra_atf_off  ***
188      !!
189      !!          !!!!!!!!!!!!!!!!! REWRITE HEADER COMMENTS !!!!!!!!!!!!!!
190      !!
191      !! ** Purpose :   Time varying volume: apply the Asselin time filter 
192      !!
193      !! ** Method  : - Apply a thickness weighted Asselin time filter on now fields.
194      !!              - save in (ta,sa) a thickness weighted average over the three
195      !!             time levels which will be used to compute rdn and thus the semi-
196      !!             implicit hydrostatic pressure gradient (ln_dynhpg_imp = T)
197      !!              - swap tracer fields to prepare the next time_step.
198      !!                This can be summurized for tempearture as:
199      !!             ztm = ( e3t_n*tn + rbcp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )   ln_dynhpg_imp = T
200      !!                  /( e3t(:,:,:,Kmm)    + rbcp*[ e3t(:,:,:,Kbb)    - 2 e3t(:,:,:,Kmm)    + e3t(:,:,:,Kaa)    ] )   
201      !!             ztm = 0                                                       otherwise
202      !!             tb  = ( e3t_n*tn + atfp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )
203      !!                  /( e3t(:,:,:,Kmm)    + atfp*[ e3t(:,:,:,Kbb)    - 2 e3t(:,:,:,Kmm)    + e3t(:,:,:,Kaa)    ] )
204      !!             tn  = ta
205      !!             ta  = zt        (NB: reset to 0 after eos_bn2 call)
206      !!
207      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step
208      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T)
209      !!----------------------------------------------------------------------
210      INTEGER                                   , INTENT(in   ) ::  kt            ! ocean time-step index
211      INTEGER                                   , INTENT(in   ) ::  Kbb, Kmm, Kaa ! time level indices
212      REAL(wp), DIMENSION(jpi,jpj,jpk,jptra,jpt), INTENT(inout) ::  ptr           ! passive tracers
213      !!     
214      INTEGER  ::   ji, jj, jk, jn              ! dummy loop indices
215      REAL(wp) ::   ztc_a , ztc_n , ztc_b , ztc_f , ztc_d    ! local scalar
216      REAL(wp) ::   ze3t_b, ze3t_n, ze3t_a, ze3t_f, ze3t_d   !   -      -
217      !!----------------------------------------------------------------------
218      !
219      IF( kt == nittrc000 )  THEN
220         IF(lwp) WRITE(numout,*)
221         IF(lwp) WRITE(numout,*) 'trc_atf_off : Asselin time filtering'
222         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
223         IF( .NOT. ln_linssh ) THEN
224            rfact1 = atfp * rdttrc
225            rfact2 = rfact1 / rau0
226         ENDIF
227        ! 
228      ENDIF
229      !
230      DO jn = 1, jptra     
231         DO_3D_11_11( 1, jpkm1 )
232            ze3t_b = e3t(ji,jj,jk,Kbb)
233            ze3t_n = e3t(ji,jj,jk,Kmm)
234            ze3t_a = e3t(ji,jj,jk,Kaa)
235            !                                         ! tracer content at Before, now and after
236            ztc_b  = ptr(ji,jj,jk,jn,Kbb)  * ze3t_b
237            ztc_n  = ptr(ji,jj,jk,jn,Kmm)  * ze3t_n
238            ztc_a  = ptr(ji,jj,jk,jn,Kaa) * ze3t_a
239            !
240            ze3t_d = ze3t_a - 2. * ze3t_n + ze3t_b
241            ztc_d  = ztc_a  - 2. * ztc_n  + ztc_b
242            !
243            ze3t_f = ze3t_n + atfp * ze3t_d
244            ztc_f  = ztc_n  + atfp * ztc_d
245            !
246            IF( .NOT. ln_linssh .AND. jk == mikt(ji,jj) ) THEN           ! first level
247               ze3t_f = ze3t_f - rfact2 * ( emp_b(ji,jj)      - emp(ji,jj)   ) 
248               ztc_f  = ztc_f  - rfact1 * ( sbc_trc(ji,jj,jn) - sbc_trc_b(ji,jj,jn) )
249            ENDIF
250
251            ze3t_f = 1.e0 / ze3t_f
252            ptr(ji,jj,jk,jn,Kmm) = ztc_f * ze3t_f     ! time filtered "now" field
253            !
254         END_3D
255         !
256      END DO
257      !
258   END SUBROUTINE trc_atf_off
259
260#else
261   !!----------------------------------------------------------------------
262   !!   Default option                                         Empty module
263   !!----------------------------------------------------------------------
264   USE par_oce
265   USE par_trc
266CONTAINS
267   SUBROUTINE trc_atf( kt, Kbb, Kmm, Kaa, ptr ) 
268      INTEGER                                   , INTENT(in)    :: kt
269      INTEGER,                                    INTENT(in   ) :: Kbb, Kmm, Kaa ! time level indices
270      REAL(wp), DIMENSION(jpi,jpj,jpk,jptra,jpt), INTENT(inout) :: ptr           ! passive tracers and RHS of tracer equation
271      WRITE(*,*) 'trc_atf: You should not have seen this print! error?', kt
272   END SUBROUTINE trc_atf
273#endif
274   !!======================================================================
275END MODULE trcatf
Note: See TracBrowser for help on using the repository browser.