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_tvd.F90 in branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

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

Last change on this file since 11101 was 11101, checked in by frrh, 5 years ago

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