source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_tvd.F90 @ 7771

Last change on this file since 7771 was 7771, checked in by frrh, 4 years ago

Apply optimisations to various areas of code replacing the use of
allocated pointers with straightforward direct ALLOCATE and DEALLOCATE
operations.

These optimisations largely have an impact in models featuring MEDUSA,
i.e. those with significant numbers of tracers, although they are
expected to have a small impact in all configurations.

Code developed and tested in NEMO branch branches/UKMO/dev_r5518_optim_GO6_alloc
Tested in stand-alone GO6-GSI8, GO6-GSI8-MEDUSA and UKESM coupled models.
NEMO ticket #1821 documents this change further.

File size: 35.0 KB
Line 
1MODULE traadv_tvd
2   !!==============================================================================
3   !!                       ***  MODULE  traadv_tvd  ***
4   !! Ocean  tracers:  horizontal & vertical advective trend
5   !!==============================================================================
6   !! History :  OPA  !  1995-12  (L. Mortier)  Original code
7   !!                 !  2000-01  (H. Loukos)  adapted to ORCA
8   !!                 !  2000-10  (MA Foujols E.Kestenare)  include file not routine
9   !!                 !  2000-12  (E. Kestenare M. Levy)  fix bug in trtrd indexes
10   !!                 !  2001-07  (E. Durand G. Madec)  adaptation to ORCA config
11   !!            8.5  !  2002-06  (G. Madec)  F90: Free form and module
12   !!    NEMO    1.0  !  2004-01  (A. de Miranda, G. Madec, J.M. Molines ): advective bbl
13   !!            2.0  !  2008-04  (S. Cravatte) add the i-, j- & k- trends computation
14   !!             -   !  2009-11  (V. Garnier) Surface pressure gradient organization
15   !!            3.3  !  2010-05  (C. Ethe, G. Madec)  merge TRC-TRA + switch from velocity to transport
16   !!----------------------------------------------------------------------
17
18   !!----------------------------------------------------------------------
19   !!   tra_adv_tvd   : update the tracer trend with the 3D advection trends using a TVD scheme
20   !!   nonosc        : compute monotonic tracer fluxes by a non-oscillatory algorithm
21   !!----------------------------------------------------------------------
22   USE oce            ! ocean dynamics and active tracers
23   USE dom_oce        ! ocean space and time domain
24   USE trc_oce        ! share passive tracers/Ocean variables
25   USE trd_oce        ! trends: ocean variables
26   USE trdtra         ! tracers trends
27   USE dynspg_oce     ! choice/control of key cpp for surface pressure gradient
28   USE diaptr         ! poleward transport diagnostics
29   !
30   USE lib_mpp        ! MPP library
31   USE lbclnk         ! ocean lateral boundary condition (or mpp link)
32   USE in_out_manager ! I/O manager
33   USE wrk_nemo       ! Memory Allocation
34   USE timing         ! Timing
35   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
36   USE iom
37
38   IMPLICIT NONE
39   PRIVATE
40
41   PUBLIC   tra_adv_tvd        ! routine called by traadv.F90
42   PUBLIC   tra_adv_tvd_zts    ! routine called by traadv.F90
43
44   LOGICAL ::   l_trd   ! flag to compute trends
45   LOGICAL ::   l_trans   ! flag to output vertically integrated transports
46
47   !! * Substitutions
48#  include "domzgr_substitute.h90"
49#  include "vectopt_loop_substitute.h90"
50   !!----------------------------------------------------------------------
51   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
52   !! $Id$
53   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
54   !!----------------------------------------------------------------------
55CONTAINS
56
57   SUBROUTINE tra_adv_tvd ( kt, kit000, cdtype, p2dt, pun, pvn, pwn,      &
58      &                                       ptb, ptn, pta, kjpt )
59      !!----------------------------------------------------------------------
60      !!                  ***  ROUTINE tra_adv_tvd  ***
61      !!
62      !! **  Purpose :   Compute the now trend due to total advection of
63      !!       tracers and add it to the general trend of tracer equations
64      !!
65      !! **  Method  :   TVD scheme, i.e. 2nd order centered scheme with
66      !!       corrected flux (monotonic correction)
67      !!       note: - this advection scheme needs a leap-frog time scheme
68      !!
69      !! ** Action : - update (pta) with the now advective tracer trends
70      !!             - save the trends
71      !!----------------------------------------------------------------------
72      USE oce     , ONLY:   zwx => ua        , zwy => va          ! (ua,va) used as workspace
73      !
74      INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index
75      INTEGER                              , INTENT(in   ) ::   kit000          ! first time step index
76      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator)
77      INTEGER                              , INTENT(in   ) ::   kjpt            ! number of tracers
78      REAL(wp), DIMENSION(        jpk     ), INTENT(in   ) ::   p2dt            ! vertical profile of tracer time-step
79      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pun, pvn, pwn   ! 3 ocean velocity components
80      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb, ptn        ! before and now tracer fields
81      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta             ! tracer trend
82      !
83      INTEGER  ::   ji, jj, jk, jn           ! dummy loop indices
84      INTEGER  ::   ik 
85      REAL(wp) ::   z2dtt, zbtr, ztra        ! local scalar
86      REAL(wp) ::   zfp_ui, zfp_vj, zfp_wk   !   -      -
87      REAL(wp) ::   zfm_ui, zfm_vj, zfm_wk   !   -      -
88      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zwi, zwz
89      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrdx, ztrdy, ztrdz, zptry
90      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   :: z2d
91      !!----------------------------------------------------------------------
92      !
93      IF( nn_timing == 1 )  CALL timing_start('tra_adv_tvd')
94      !
95      ALLOCATE(zwi(1:jpi, 1:jpj, 1:jpk))
96      ALLOCATE(zwz(1:jpi, 1:jpj, 1:jpk))
97
98      !
99      IF( kt == kit000 )  THEN
100         IF(lwp) WRITE(numout,*)
101         IF(lwp) WRITE(numout,*) 'tra_adv_tvd : TVD advection scheme on ', cdtype
102         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
103         !
104      ENDIF
105      l_trd = .FALSE.
106      l_trans = .FALSE.
107      IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.
108      IF( cdtype == 'TRA' .AND. (iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") ) ) l_trans = .TRUE.
109      !
110      IF( l_trd .OR. l_trans )  THEN
111         ALLOCATE(ztrdx(1:jpi, 1:jpj, 1:jpk))
112         ALLOCATE(ztrdy(1:jpi, 1:jpj, 1:jpk))
113         ALLOCATE(ztrdz(1:jpi, 1:jpj, 1:jpk))
114         ztrdx(:,:,:) = 0.e0   ;    ztrdy(:,:,:) = 0.e0   ;   ztrdz(:,:,:) = 0.e0
115         ALLOCATE(z2d(1:jpi, 1:jpj))
116      ENDIF
117      !
118      IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN 
119         ALLOCATE(zptry(1:jpi, 1:jpj, 1:jpk))
120         zptry(:,:,:) = 0._wp
121      ENDIF
122      !
123      zwi(:,:,:) = 0.e0 ; 
124      !
125      !                                                          ! ===========
126      DO jn = 1, kjpt                                            ! tracer loop
127         !                                                       ! ===========
128         ! 1. Bottom and k=1 value : flux set to zero
129         ! ----------------------------------
130         zwx(:,:,jpk) = 0.e0    ;    zwz(:,:,jpk) = 0.e0
131         zwy(:,:,jpk) = 0.e0    ;    zwi(:,:,jpk) = 0.e0
132         
133         zwz(:,:,1  ) = 0._wp
134         ! 2. upstream advection with initial mass fluxes & intermediate update
135         ! --------------------------------------------------------------------
136         ! upstream tracer flux in the i and j direction
137         DO jk = 1, jpkm1
138            DO jj = 1, jpjm1
139               DO ji = 1, fs_jpim1   ! vector opt.
140                  ! upstream scheme
141                  zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) )
142                  zfm_ui = pun(ji,jj,jk) - ABS( pun(ji,jj,jk) )
143                  zfp_vj = pvn(ji,jj,jk) + ABS( pvn(ji,jj,jk) )
144                  zfm_vj = pvn(ji,jj,jk) - ABS( pvn(ji,jj,jk) )
145                  zwx(ji,jj,jk) = 0.5 * ( zfp_ui * ptb(ji,jj,jk,jn) + zfm_ui * ptb(ji+1,jj  ,jk,jn) )
146                  zwy(ji,jj,jk) = 0.5 * ( zfp_vj * ptb(ji,jj,jk,jn) + zfm_vj * ptb(ji  ,jj+1,jk,jn) )
147               END DO
148            END DO
149         END DO
150
151         ! upstream tracer flux in the k direction
152         ! Interior value
153         DO jk = 2, jpkm1
154            DO jj = 1, jpj
155               DO ji = 1, jpi
156                  zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) )
157                  zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) )
158                  zwz(ji,jj,jk) = 0.5 * ( zfp_wk * ptb(ji,jj,jk,jn) + zfm_wk * ptb(ji,jj,jk-1,jn) ) * wmask(ji,jj,jk)
159               END DO
160            END DO
161         END DO
162         ! Surface value
163         IF( lk_vvl ) THEN   
164            IF ( ln_isfcav ) THEN
165               DO jj = 1, jpj
166                  DO ji = 1, jpi
167                     zwz(ji,jj, mikt(ji,jj) ) = 0.e0          ! volume variable
168                  END DO
169               END DO
170            ELSE
171               zwz(:,:,1) = 0.e0          ! volume variable
172            END IF
173         ELSE               
174            IF ( ln_isfcav ) THEN
175               DO jj = 1, jpj
176                  DO ji = 1, jpi
177                     zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn)   ! linear free surface
178                  END DO
179               END DO   
180            ELSE
181               zwz(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn)   ! linear free surface
182            END IF
183         ENDIF
184
185         ! total advective trend
186         DO jk = 1, jpkm1
187            z2dtt = p2dt(jk)
188            DO jj = 2, jpjm1
189               DO ji = fs_2, fs_jpim1   ! vector opt.
190                  ! total intermediate advective trends
191                  ztra = - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
192                     &      + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   &
193                     &      + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) / e1e2t(ji,jj)
194                  ! update and guess with monotonic sheme
195                  pta(ji,jj,jk,jn) =                       pta(ji,jj,jk,jn) +         ztra   / fse3t_n(ji,jj,jk) * tmask(ji,jj,jk)
196                  zwi(ji,jj,jk)    = ( fse3t_b(ji,jj,jk) * ptb(ji,jj,jk,jn) + z2dtt * ztra ) / fse3t_a(ji,jj,jk) * tmask(ji,jj,jk)
197               END DO
198            END DO
199         END DO
200         !                             ! Lateral boundary conditions on zwi  (unchanged sign)
201         CALL lbc_lnk( zwi, 'T', 1. ) 
202
203         !                                 ! trend diagnostics (contribution of upstream fluxes)
204         IF( l_trd .OR. l_trans )  THEN 
205            ! store intermediate advective trends
206            ztrdx(:,:,:) = zwx(:,:,:)   ;    ztrdy(:,:,:) = zwy(:,:,:)  ;   ztrdz(:,:,:) = zwz(:,:,:)
207         END IF
208         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
209         IF( cdtype == 'TRA' .AND. ln_diaptr )    zptry(:,:,:) = zwy(:,:,:) 
210
211         ! 3. antidiffusive flux : high order minus low order
212         ! --------------------------------------------------
213         ! antidiffusive flux on i and j
214         DO jk = 1, jpkm1
215            DO jj = 1, jpjm1
216               DO ji = 1, fs_jpim1   ! vector opt.
217                  zwx(ji,jj,jk) = 0.5 * pun(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj,jk,jn) ) - zwx(ji,jj,jk)
218                  zwy(ji,jj,jk) = 0.5 * pvn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj+1,jk,jn) ) - zwy(ji,jj,jk)
219               END DO
220            END DO
221         END DO
222     
223         ! antidiffusive flux on k
224         ! Interior value
225         DO jk = 2, jpkm1                   
226            DO jj = 1, jpj
227               DO ji = 1, jpi
228                  zwz(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj,jk-1,jn) ) - zwz(ji,jj,jk)
229               END DO
230            END DO
231         END DO
232         ! surface value
233         IF ( ln_isfcav ) THEN
234            DO jj = 1, jpj
235               DO ji = 1, jpi
236                  zwz(ji,jj,mikt(ji,jj)) = 0.e0
237               END DO
238            END DO
239         ELSE
240            zwz(:,:,1) = 0.e0
241         END IF
242         CALL lbc_lnk( zwx, 'U', -1. )   ;   CALL lbc_lnk( zwy, 'V', -1. )         ! Lateral bondary conditions
243         CALL lbc_lnk( zwz, 'W',  1. )
244
245         ! 4. monotonicity algorithm
246         ! -------------------------
247         CALL nonosc( ptb(:,:,:,jn), zwx, zwy, zwz, zwi, p2dt )
248
249
250         ! 5. final trend with corrected fluxes
251         ! ------------------------------------
252         DO jk = 1, jpkm1
253            DO jj = 2, jpjm1
254               DO ji = fs_2, fs_jpim1   ! vector opt. 
255                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
256                  ! total advective trends
257                  ztra = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
258                     &             + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   &
259                     &             + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) )
260                  ! add them to the general tracer trends
261                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra * tmask(ji,jj,jk)
262               END DO
263            END DO
264         END DO
265
266         !                                 ! trend diagnostics (contribution of upstream fluxes)
267         IF( l_trd .OR. l_trans )  THEN
268            ztrdx(:,:,:) = ztrdx(:,:,:) + zwx(:,:,:)  ! <<< Add to previously computed
269            ztrdy(:,:,:) = ztrdy(:,:,:) + zwy(:,:,:)  ! <<< Add to previously computed
270            ztrdz(:,:,:) = ztrdz(:,:,:) + zwz(:,:,:)  ! <<< Add to previously computed
271         ENDIF
272         
273         IF( l_trd ) THEN
274            CALL trd_tra( kt, cdtype, jn, jptra_xad, ztrdx, pun, ptn(:,:,:,jn) )
275            CALL trd_tra( kt, cdtype, jn, jptra_yad, ztrdy, pvn, ptn(:,:,:,jn) )
276            CALL trd_tra( kt, cdtype, jn, jptra_zad, ztrdz, pwn, ptn(:,:,:,jn) )
277         END IF
278
279         IF( l_trans .AND. jn==jp_tem ) THEN
280            z2d(:,:) = 0._wp 
281            DO jk = 1, jpkm1
282               DO jj = 2, jpjm1
283                  DO ji = fs_2, fs_jpim1   ! vector opt.
284                     z2d(ji,jj) = z2d(ji,jj) + ztrdx(ji,jj,jk) 
285                  END DO
286               END DO
287            END DO
288            CALL lbc_lnk( z2d, 'U', -1. )
289            CALL iom_put( "uadv_heattr", rau0_rcp * z2d )       ! heat transport in i-direction
290              !
291            z2d(:,:) = 0._wp 
292            DO jk = 1, jpkm1
293               DO jj = 2, jpjm1
294                  DO ji = fs_2, fs_jpim1   ! vector opt.
295                     z2d(ji,jj) = z2d(ji,jj) + ztrdy(ji,jj,jk) 
296                  END DO
297               END DO
298            END DO
299            CALL lbc_lnk( z2d, 'V', -1. )
300            CALL iom_put( "vadv_heattr", rau0_rcp * z2d )       ! heat transport in j-direction
301         ENDIF
302         ! "Poleward" heat and salt transports (contribution of upstream fluxes)
303         IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN 
304            zptry(:,:,:) = zptry(:,:,:) + zwy(:,:,:)  ! <<< Add to previously computed
305            CALL dia_ptr_ohst_components( jn, 'adv', zptry(:,:,:) )
306         ENDIF
307         !
308      END DO
309      !
310      DEALLOCATE( zwi )
311      DEALLOCATE( zwz )
312      IF( l_trd .OR. l_trans )  THEN
313         DEALLOCATE( ztrdx )
314         DEALLOCATE( ztrdy )
315         DEALLOCATE( ztrdz )
316         DEALLOCATE( z2d )
317      ENDIF
318      IF( cdtype == 'TRA' .AND. ln_diaptr ) DEALLOCATE( zptry )
319      !
320      IF( nn_timing == 1 )  CALL timing_stop('tra_adv_tvd')
321      !
322   END SUBROUTINE tra_adv_tvd
323
324   SUBROUTINE tra_adv_tvd_zts ( kt, kit000, cdtype, p2dt, pun, pvn, pwn,      &
325      &                                       ptb, ptn, pta, kjpt )
326      !!----------------------------------------------------------------------
327      !!                  ***  ROUTINE tra_adv_tvd_zts  ***
328      !!
329      !! **  Purpose :   Compute the now trend due to total advection of
330      !!       tracers and add it to the general trend of tracer equations
331      !!
332      !! **  Method  :   TVD ZTS scheme, i.e. 2nd order centered scheme with
333      !!       corrected flux (monotonic correction). This version use sub-
334      !!       timestepping for the vertical advection which increases stability
335      !!       when vertical metrics are small.
336      !!       note: - this advection scheme needs a leap-frog time scheme
337      !!
338      !! ** Action : - update (pta) with the now advective tracer trends
339      !!             - save the trends
340      !!----------------------------------------------------------------------
341      USE oce     , ONLY:   zwx => ua        , zwy => va          ! (ua,va) used as workspace
342      !
343      INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index
344      INTEGER                              , INTENT(in   ) ::   kit000          ! first time step index
345      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator)
346      INTEGER                              , INTENT(in   ) ::   kjpt            ! number of tracers
347      REAL(wp), DIMENSION(        jpk     ), INTENT(in   ) ::   p2dt            ! vertical profile of tracer time-step
348      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pun, pvn, pwn   ! 3 ocean velocity components
349      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb, ptn        ! before and now tracer fields
350      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta             ! tracer trend
351      !
352      REAL(wp), DIMENSION( jpk )                           ::   zts             ! length of sub-timestep for vertical advection
353      REAL(wp), DIMENSION( jpk )                           ::   zr_p2dt         ! reciprocal of tracer timestep
354      INTEGER  ::   ji, jj, jk, jl, jn       ! dummy loop indices 
355      INTEGER  ::   jnzts = 5       ! number of sub-timesteps for vertical advection
356      INTEGER  ::   jtb, jtn, jta   ! sub timestep pointers for leap-frog/euler forward steps
357      INTEGER  ::   jtaken          ! toggle for collecting appropriate fluxes from sub timesteps
358      REAL(wp) ::   z_rzts          ! Fractional length of Euler forward sub-timestep for vertical advection
359      REAL(wp) ::   z2dtt, zbtr, ztra        ! local scalar
360      REAL(wp) ::   zfp_ui, zfp_vj, zfp_wk   !   -      -
361      REAL(wp) ::   zfm_ui, zfm_vj, zfm_wk   !   -      -
362      REAL(wp), ALLOCATABLE, DIMENSION(:,:  ) :: zwx_sav , zwy_sav
363      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zwi, zwz, zhdiv, zwz_sav, zwzts
364      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrdx, ztrdy, ztrdz
365      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zptry
366      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) :: ztrs
367      !!----------------------------------------------------------------------
368      !
369      IF( nn_timing == 1 )  CALL timing_start('tra_adv_tvd_zts')
370      !
371      ALLOCATE(zwx_sav(1:jpi, 1:jpj))
372      ALLOCATE(zwy_sav(1:jpi, 1:jpj))
373      ALLOCATE(zwi(1:jpi, 1:jpj, 1:jpk))
374      ALLOCATE(zwz(1:jpi, 1:jpj, 1:jpk))       
375      ALLOCATE(zhdiv(1:jpi, 1:jpj, 1:jpk))       
376      ALLOCATE(zwz_sav(1:jpi, 1:jpj, 1:jpk))       
377      ALLOCATE(zwzts(1:jpi, 1:jpj, 1:jpk)) 
378      ALLOCATE(ztrs(1:jpi, 1:jpj, 1:jpk, 1:kjpt+1))
379      !
380      IF( kt == kit000 )  THEN
381         IF(lwp) WRITE(numout,*)
382         IF(lwp) WRITE(numout,*) 'tra_adv_tvd_zts : TVD ZTS advection scheme on ', cdtype
383         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
384      ENDIF
385      !
386      l_trd = .FALSE.
387      IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.
388      !
389      IF( l_trd )  THEN
390         ALLOCATE(ztrdx(1:jpi, 1:jpj, 1:jpk))       
391         ALLOCATE(ztrdy(1:jpi, 1:jpj, 1:jpk))       
392         ALLOCATE(ztrdz(1:jpi, 1:jpj, 1:jpk))       
393         ztrdx(:,:,:) = 0._wp  ;    ztrdy(:,:,:) = 0._wp  ;   ztrdz(:,:,:) = 0._wp
394      ENDIF
395      !
396      IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN 
397         ALLOCATE(zptry(1:jpi, 1:jpj, 1:jpk))       
398         zptry(:,:,:) = 0._wp
399      ENDIF
400      !
401      zwi(:,:,:) = 0._wp
402      z_rzts = 1._wp / REAL( jnzts, wp )
403      zr_p2dt(:) = 1._wp / p2dt(:)
404      !
405      !                                                          ! ===========
406      DO jn = 1, kjpt                                            ! tracer loop
407         !                                                       ! ===========
408         ! 1. Bottom value : flux set to zero
409         ! ----------------------------------
410         zwx(:,:,jpk) = 0._wp   ;    zwz(:,:,jpk) = 0._wp
411         zwy(:,:,jpk) = 0._wp   ;    zwi(:,:,jpk) = 0._wp
412
413         ! 2. upstream advection with initial mass fluxes & intermediate update
414         ! --------------------------------------------------------------------
415         ! upstream tracer flux in the i and j direction
416         DO jk = 1, jpkm1
417            DO jj = 1, jpjm1
418               DO ji = 1, fs_jpim1   ! vector opt.
419                  ! upstream scheme
420                  zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) )
421                  zfm_ui = pun(ji,jj,jk) - ABS( pun(ji,jj,jk) )
422                  zfp_vj = pvn(ji,jj,jk) + ABS( pvn(ji,jj,jk) )
423                  zfm_vj = pvn(ji,jj,jk) - ABS( pvn(ji,jj,jk) )
424                  zwx(ji,jj,jk) = 0.5_wp * ( zfp_ui * ptb(ji,jj,jk,jn) + zfm_ui * ptb(ji+1,jj  ,jk,jn) )
425                  zwy(ji,jj,jk) = 0.5_wp * ( zfp_vj * ptb(ji,jj,jk,jn) + zfm_vj * ptb(ji  ,jj+1,jk,jn) )
426               END DO
427            END DO
428         END DO
429
430         ! upstream tracer flux in the k direction
431         ! Interior value
432         DO jk = 2, jpkm1
433            DO jj = 1, jpj
434               DO ji = 1, jpi
435                  zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) )
436                  zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) )
437                  zwz(ji,jj,jk) = 0.5_wp * ( zfp_wk * ptb(ji,jj,jk,jn) + zfm_wk * ptb(ji,jj,jk-1,jn) )
438               END DO
439            END DO
440         END DO
441         ! Surface value
442         IF( lk_vvl ) THEN
443            IF ( ln_isfcav ) THEN
444               DO jj = 1, jpj
445                  DO ji = 1, jpi
446                     zwz(ji,jj, mikt(ji,jj) ) = 0.e0          ! volume variable +    isf
447                  END DO
448               END DO
449            ELSE
450               zwz(:,:,1) = 0.e0                              ! volume variable + no isf
451            END IF
452         ELSE
453            IF ( ln_isfcav ) THEN
454               DO jj = 1, jpj
455                  DO ji = 1, jpi
456                     zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn)   ! linear free surface +    isf
457                  END DO
458               END DO
459            ELSE
460               zwz(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn)                                               ! linear free surface + no isf
461            END IF
462         ENDIF
463
464         ! total advective trend
465         DO jk = 1, jpkm1
466            z2dtt = p2dt(jk)
467            DO jj = 2, jpjm1
468               DO ji = fs_2, fs_jpim1   ! vector opt.
469                  ! total intermediate advective trends
470                  ztra = - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
471                     &      + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   &
472                     &      + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) / e1e2t(ji,jj)
473                  ! update and guess with monotonic sheme
474                  pta(ji,jj,jk,jn) =                       pta(ji,jj,jk,jn) +         ztra   / fse3t_n(ji,jj,jk) * tmask(ji,jj,jk)
475                  zwi(ji,jj,jk)    = ( fse3t_b(ji,jj,jk) * ptb(ji,jj,jk,jn) + z2dtt * ztra ) / fse3t_a(ji,jj,jk) * tmask(ji,jj,jk)
476               END DO
477            END DO
478         END DO
479         !                             ! Lateral boundary conditions on zwi  (unchanged sign)
480         CALL lbc_lnk( zwi, 'T', 1. ) 
481
482         !                                 ! trend diagnostics (contribution of upstream fluxes)
483         IF( l_trd )  THEN 
484            ! store intermediate advective trends
485            ztrdx(:,:,:) = zwx(:,:,:)   ;    ztrdy(:,:,:) = zwy(:,:,:)  ;   ztrdz(:,:,:) = zwz(:,:,:)
486         END IF
487         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
488         IF( cdtype == 'TRA' .AND. ln_diaptr )  zptry(:,:,:) = zwy(:,:,:)
489
490         ! 3. antidiffusive flux : high order minus low order
491         ! --------------------------------------------------
492         ! antidiffusive flux on i and j
493         !
494         DO jk = 1, jpkm1
495            !
496            DO jj = 1, jpjm1
497               DO ji = 1, fs_jpim1   ! vector opt.
498                  zwx_sav(ji,jj) = zwx(ji,jj,jk)
499                  zwy_sav(ji,jj) = zwy(ji,jj,jk)
500
501                  zwx(ji,jj,jk) = 0.5_wp * pun(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj,jk,jn) )
502                  zwy(ji,jj,jk) = 0.5_wp * pvn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj+1,jk,jn) )
503               END DO
504            END DO
505
506            DO jj = 2, jpjm1         ! partial horizontal divergence
507               DO ji = fs_2, fs_jpim1
508                  zhdiv(ji,jj,jk) = (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk)   &
509                     &               + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk)  )
510               END DO
511            END DO
512
513            DO jj = 1, jpjm1
514               DO ji = 1, fs_jpim1   ! vector opt.
515                  zwx(ji,jj,jk) = zwx(ji,jj,jk)  - zwx_sav(ji,jj)
516                  zwy(ji,jj,jk) = zwy(ji,jj,jk)  - zwy_sav(ji,jj)
517               END DO
518            END DO
519         END DO
520     
521         ! antidiffusive flux on k
522         zwz(:,:,1) = 0._wp        ! Surface value
523         zwz_sav(:,:,:) = zwz(:,:,:)
524         !
525         ztrs(:,:,:,1) = ptb(:,:,:,jn)
526         ztrs(:,:,1,2) = ptb(:,:,1,jn)
527         ztrs(:,:,1,3) = ptb(:,:,1,jn)
528         zwzts(:,:,:) = 0._wp
529
530         DO jl = 1, jnzts                   ! Start of sub timestepping loop
531
532            IF( jl == 1 ) THEN              ! Euler forward to kick things off
533              jtb = 1   ;   jtn = 1   ;   jta = 2
534              zts(:) = p2dt(:) * z_rzts
535              jtaken = MOD( jnzts + 1 , 2)  ! Toggle to collect every second flux
536                                            ! starting at jl =1 if jnzts is odd;
537                                            ! starting at jl =2 otherwise
538            ELSEIF( jl == 2 ) THEN          ! First leapfrog step
539              jtb = 1   ;   jtn = 2   ;   jta = 3
540              zts(:) = 2._wp * p2dt(:) * z_rzts
541            ELSE                            ! Shuffle pointers for subsequent leapfrog steps
542              jtb = MOD(jtb,3) + 1
543              jtn = MOD(jtn,3) + 1
544              jta = MOD(jta,3) + 1
545            ENDIF
546            DO jk = 2, jpkm1          ! Interior value
547               DO jj = 2, jpjm1
548                  DO ji = fs_2, fs_jpim1
549                     zwz(ji,jj,jk) = 0.5_wp * pwn(ji,jj,jk) * ( ztrs(ji,jj,jk,jtn) + ztrs(ji,jj,jk-1,jtn) )
550                     IF( jtaken == 0 ) zwzts(ji,jj,jk) = zwzts(ji,jj,jk) + zwz(ji,jj,jk)*zts(jk)           ! Accumulate time-weighted vertcal flux
551                  END DO
552               END DO
553            END DO
554
555            jtaken = MOD( jtaken + 1 , 2 )
556
557            DO jk = 2, jpkm1          ! Interior value
558               DO jj = 2, jpjm1
559                  DO ji = fs_2, fs_jpim1
560                     zbtr = 1._wp / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
561                     ! total advective trends
562                     ztra = - zbtr * (  zhdiv(ji,jj,jk) + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) )
563                     ztrs(ji,jj,jk,jta) = ztrs(ji,jj,jk,jtb) + zts(jk) * ztra
564                  END DO
565               END DO
566            END DO
567
568         END DO
569
570         DO jk = 2, jpkm1          ! Anti-diffusive vertical flux using average flux from the sub-timestepping
571            DO jj = 2, jpjm1
572               DO ji = fs_2, fs_jpim1
573                  zwz(ji,jj,jk) = zwzts(ji,jj,jk) * zr_p2dt(jk) - zwz_sav(ji,jj,jk)
574               END DO
575            END DO
576         END DO
577         CALL lbc_lnk( zwx, 'U', -1. )   ;   CALL lbc_lnk( zwy, 'V', -1. )         ! Lateral bondary conditions
578         CALL lbc_lnk( zwz, 'W',  1. )
579
580         ! 4. monotonicity algorithm
581         ! -------------------------
582         CALL nonosc( ptb(:,:,:,jn), zwx, zwy, zwz, zwi, p2dt )
583
584
585         ! 5. final trend with corrected fluxes
586         ! ------------------------------------
587         DO jk = 1, jpkm1
588            DO jj = 2, jpjm1
589               DO ji = fs_2, fs_jpim1   ! vector opt. 
590                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
591                  ! total advective trends
592                  ztra = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
593                     &             + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   &
594                     &             + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) )
595                  ! add them to the general tracer trends
596                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
597               END DO
598            END DO
599         END DO
600
601         !                                 ! trend diagnostics (contribution of upstream fluxes)
602         IF( l_trd )  THEN
603            ztrdx(:,:,:) = ztrdx(:,:,:) + zwx(:,:,:)  ! <<< Add to previously computed
604            ztrdy(:,:,:) = ztrdy(:,:,:) + zwy(:,:,:)  ! <<< Add to previously computed
605            ztrdz(:,:,:) = ztrdz(:,:,:) + zwz(:,:,:)  ! <<< Add to previously computed
606           
607            CALL trd_tra( kt, cdtype, jn, jptra_xad, ztrdx, pun, ptn(:,:,:,jn) )   
608            CALL trd_tra( kt, cdtype, jn, jptra_yad, ztrdy, pvn, ptn(:,:,:,jn) ) 
609            CALL trd_tra( kt, cdtype, jn, jptra_zad, ztrdz, pwn, ptn(:,:,:,jn) ) 
610         END IF
611         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
612         IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN 
613            zptry(:,:,:) = zptry(:,:,:) + zwy(:,:,:) 
614            CALL dia_ptr_ohst_components( jn, 'adv', zptry(:,:,:) )
615         ENDIF
616         !
617      END DO
618      !
619      DEALLOCATE(zwi) 
620      DEALLOCATE(zwz) 
621      DEALLOCATE(zhdiv) 
622      DEALLOCATE(zwz_sav) 
623      DEALLOCATE(zwzts)
624      DEALLOCATE(ztrs )
625      DEALLOCATE(zwx_sav) 
626      DEALLOCATE(zwy_sav )
627
628      IF( l_trd )  THEN
629          DEALLOCATE(ztrdx) 
630          DEALLOCATE(ztrdy) 
631          DEALLOCATE(ztrdz)
632      END IF
633
634      IF( cdtype == 'TRA' .AND. ln_diaptr ) DEALLOCATE(zptry )
635      !
636      IF( nn_timing == 1 )  CALL timing_stop('tra_adv_tvd_zts')
637      !
638   END SUBROUTINE tra_adv_tvd_zts
639
640
641   SUBROUTINE nonosc( pbef, paa, pbb, pcc, paft, p2dt )
642      !!---------------------------------------------------------------------
643      !!                    ***  ROUTINE nonosc  ***
644      !!     
645      !! **  Purpose :   compute monotonic tracer fluxes from the upstream
646      !!       scheme and the before field by a nonoscillatory algorithm
647      !!
648      !! **  Method  :   ... ???
649      !!       warning : pbef and paft must be masked, but the boundaries
650      !!       conditions on the fluxes are not necessary zalezak (1979)
651      !!       drange (1995) multi-dimensional forward-in-time and upstream-
652      !!       in-space based differencing for fluid
653      !!----------------------------------------------------------------------
654      REAL(wp), DIMENSION(jpk)         , INTENT(in   ) ::   p2dt            ! vertical profile of tracer time-step
655      REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(in   ) ::   pbef, paft      ! before & after field
656      REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(inout) ::   paa, pbb, pcc   ! monotonic fluxes in the 3 directions
657      !
658      INTEGER  ::   ji, jj, jk   ! dummy loop indices
659      INTEGER  ::   ikm1         ! local integer
660      REAL(wp) ::   zpos, zneg, zbt, za, zb, zc, zbig, zrtrn, z2dtt   ! local scalars
661      REAL(wp) ::   zau, zbu, zcu, zav, zbv, zcv, zup, zdo            !   -      -
662      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zbetup, zbetdo, zbup, zbdo
663      !!----------------------------------------------------------------------
664      !
665      IF( nn_timing == 1 )  CALL timing_start('nonosc')
666      !
667      ALLOCATE(zbetup(1:jpi, 1:jpj, 1:jpk))
668      ALLOCATE(zbetdo(1:jpi, 1:jpj, 1:jpk))
669      ALLOCATE(zbup(1:jpi, 1:jpj, 1:jpk))
670      ALLOCATE(zbdo(1:jpi, 1:jpj, 1:jpk))
671      !
672      zbig  = 1.e+40_wp
673      zrtrn = 1.e-15_wp
674      zbetup(:,:,:) = 0._wp   ;   zbetdo(:,:,:) = 0._wp
675
676      ! Search local extrema
677      ! --------------------
678      ! max/min of pbef & paft with large negative/positive value (-/+zbig) inside land
679      zbup = MAX( pbef * tmask - zbig * ( 1._wp - tmask ),   &
680         &        paft * tmask - zbig * ( 1._wp - tmask )  )
681      zbdo = MIN( pbef * tmask + zbig * ( 1._wp - tmask ),   &
682         &        paft * tmask + zbig * ( 1._wp - tmask )  )
683
684      DO jk = 1, jpkm1
685         ikm1 = MAX(jk-1,1)
686         z2dtt = p2dt(jk)
687         DO jj = 2, jpjm1
688            DO ji = fs_2, fs_jpim1   ! vector opt.
689
690               ! search maximum in neighbourhood
691               zup = MAX(  zbup(ji  ,jj  ,jk  ),   &
692                  &        zbup(ji-1,jj  ,jk  ), zbup(ji+1,jj  ,jk  ),   &
693                  &        zbup(ji  ,jj-1,jk  ), zbup(ji  ,jj+1,jk  ),   &
694                  &        zbup(ji  ,jj  ,ikm1), zbup(ji  ,jj  ,jk+1)  )
695
696               ! search minimum in neighbourhood
697               zdo = MIN(  zbdo(ji  ,jj  ,jk  ),   &
698                  &        zbdo(ji-1,jj  ,jk  ), zbdo(ji+1,jj  ,jk  ),   &
699                  &        zbdo(ji  ,jj-1,jk  ), zbdo(ji  ,jj+1,jk  ),   &
700                  &        zbdo(ji  ,jj  ,ikm1), zbdo(ji  ,jj  ,jk+1)  )
701
702               ! positive part of the flux
703               zpos = MAX( 0., paa(ji-1,jj  ,jk  ) ) - MIN( 0., paa(ji  ,jj  ,jk  ) )   &
704                  & + MAX( 0., pbb(ji  ,jj-1,jk  ) ) - MIN( 0., pbb(ji  ,jj  ,jk  ) )   &
705                  & + MAX( 0., pcc(ji  ,jj  ,jk+1) ) - MIN( 0., pcc(ji  ,jj  ,jk  ) )
706
707               ! negative part of the flux
708               zneg = MAX( 0., paa(ji  ,jj  ,jk  ) ) - MIN( 0., paa(ji-1,jj  ,jk  ) )   &
709                  & + MAX( 0., pbb(ji  ,jj  ,jk  ) ) - MIN( 0., pbb(ji  ,jj-1,jk  ) )   &
710                  & + MAX( 0., pcc(ji  ,jj  ,jk  ) ) - MIN( 0., pcc(ji  ,jj  ,jk+1) )
711
712               ! up & down beta terms
713               zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) / z2dtt
714               zbetup(ji,jj,jk) = ( zup            - paft(ji,jj,jk) ) / ( zpos + zrtrn ) * zbt
715               zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zdo            ) / ( zneg + zrtrn ) * zbt
716            END DO
717         END DO
718      END DO
719      CALL lbc_lnk( zbetup, 'T', 1. )   ;   CALL lbc_lnk( zbetdo, 'T', 1. )   ! lateral boundary cond. (unchanged sign)
720
721      ! 3. monotonic flux in the i & j direction (paa & pbb)
722      ! ----------------------------------------
723      DO jk = 1, jpkm1
724         DO jj = 2, jpjm1
725            DO ji = fs_2, fs_jpim1   ! vector opt.
726               zau = MIN( 1._wp, zbetdo(ji,jj,jk), zbetup(ji+1,jj,jk) )
727               zbu = MIN( 1._wp, zbetup(ji,jj,jk), zbetdo(ji+1,jj,jk) )
728               zcu =       ( 0.5  + SIGN( 0.5 , paa(ji,jj,jk) ) )
729               paa(ji,jj,jk) = paa(ji,jj,jk) * ( zcu * zau + ( 1._wp - zcu) * zbu )
730
731               zav = MIN( 1._wp, zbetdo(ji,jj,jk), zbetup(ji,jj+1,jk) )
732               zbv = MIN( 1._wp, zbetup(ji,jj,jk), zbetdo(ji,jj+1,jk) )
733               zcv =       ( 0.5  + SIGN( 0.5 , pbb(ji,jj,jk) ) )
734               pbb(ji,jj,jk) = pbb(ji,jj,jk) * ( zcv * zav + ( 1._wp - zcv) * zbv )
735
736      ! monotonic flux in the k direction, i.e. pcc
737      ! -------------------------------------------
738               za = MIN( 1., zbetdo(ji,jj,jk+1), zbetup(ji,jj,jk) )
739               zb = MIN( 1., zbetup(ji,jj,jk+1), zbetdo(ji,jj,jk) )
740               zc =       ( 0.5  + SIGN( 0.5 , pcc(ji,jj,jk+1) ) )
741               pcc(ji,jj,jk+1) = pcc(ji,jj,jk+1) * ( zc * za + ( 1._wp - zc) * zb )
742            END DO
743         END DO
744      END DO
745      CALL lbc_lnk( paa, 'U', -1. )   ;   CALL lbc_lnk( pbb, 'V', -1. )   ! lateral boundary condition (changed sign)
746      !
747      DEALLOCATE(zbetup)
748      DEALLOCATE(zbetdo) 
749      DEALLOCATE(zbup)
750      DEALLOCATE(zbdo)
751      !
752      IF( nn_timing == 1 )  CALL timing_stop('nonosc')
753      !
754   END SUBROUTINE nonosc
755
756   !!======================================================================
757END MODULE traadv_tvd
Note: See TracBrowser for help on using the repository browser.