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/dev_001_GM/NEMO/OPA_SRC/TRA – NEMO

source: branches/dev_001_GM/NEMO/OPA_SRC/TRA/traadv_tvd.F90 @ 791

Last change on this file since 791 was 791, checked in by gm, 16 years ago

dev_001_GM - TRA/traadv : switch from velocity to transport + optimised traadv_eiv2 introduced - compilation OK

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 17.6 KB
Line 
1MODULE traadv_tvd
2   !!==============================================================================
3   !!                       ***  MODULE  traadv_tvd  ***
4   !! Ocean active tracers:  horizontal & vertical advective trend
5   !!==============================================================================
6   !! History :  7.0  !  1995-12  (L. Mortier)  Original code
7   !!            8.0  !  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   !!             -   !  2008-04  (S. Cravatte) add the i-, j- & k- trends computation
14   !!             -   !  2005-11  (V. Garnier) Surface pressure gradient organization
15   !!            2.4  !  2008-01  (G. Madec)  merge TRC-TRA + switch from velocity to transport
16   !!----------------------------------------------------------------------
17
18   !!----------------------------------------------------------------------
19   !!   tra_adv_tvd  : update the tracer trend with the horizontal and
20   !!                  vertical advection trends using a TVD scheme
21   !!   nonosc       : compute monotonic tracer fluxes by a nonoscillatory algorithm
22   !!----------------------------------------------------------------------
23   USE dom_oce         ! ocean space and time domain
24   USE trdmod          ! ocean active tracers trends
25   USE trdmod_oce      ! ocean variables trends
26   USE in_out_manager  ! I/O manager
27   USE dynspg_oce      ! choice/control of key cpp for surface pressure gradient
28   USE trabbl          ! Advective term of BBL
29   USE lib_mpp         !
30   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
31   USE diaptr          ! poleward transport diagnostics
32   USE prtctl          ! Print control
33
34   IMPLICIT NONE
35   PRIVATE
36
37   PUBLIC   tra_adv_tvd    ! routine called by traadv.F90
38
39   !! * Substitutions
40#  include "domzgr_substitute.h90"
41#  include "vectopt_loop_substitute.h90"
42   !!----------------------------------------------------------------------
43   !! NEMO/OPA 2.4 , LOCEAN-IPSL (2008)
44   !! $Id$
45   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
46   !!----------------------------------------------------------------------
47
48CONTAINS
49
50   SUBROUTINE tra_adv_tvd( kt, cdtype, ktra, pun, pvn, pwn,   &
51      &                                      ptb, ptn, pta )
52      !!----------------------------------------------------------------------
53      !!                  ***  ROUTINE tra_adv_tvd  ***
54      !!
55      !! **  Purpose :   Compute the now trend due to total advection of
56      !!       tracers and add it to the general trend of tracer equations
57      !!
58      !! **  Method  :   TVD scheme, i.e. 2nd order centered scheme with
59      !!       corrected flux (monotonic correction)
60      !!       note: - this advection scheme needs a leap-frog time scheme
61      !!
62      !! ** Action : - update pta with the now advective tracer trends
63      !!             - save the trends in (ztrdt,ztrds) ('key_trdtra')
64      !!----------------------------------------------------------------------
65      INTEGER         , INTENT(in   )                         ::   kt              ! ocean time-step index
66      CHARACTER(len=3), INTENT(in   )                         ::   cdtype          ! =TRA or TRC (tracer indicator)
67      INTEGER         , INTENT(in   )                         ::   ktra            ! tracer index
68      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   pun, pvn, pwn   ! 3 ocean velocity components
69      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   ptb, ptn        ! before and now tracer fields
70      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pta             ! tracer trend
71      !!
72      INTEGER  ::   ji, jj, jk               ! dummy loop indices
73      REAL(wp) ::   z2dtt, zbtr, z2, zzti    ! temporary scalar
74      REAL(wp) ::   zfp_ui, zfp_vj, zfp_wk   !    "         "
75      REAL(wp) ::   zfm_ui, zfm_vj, zfm_wk   !    "         "
76      REAL(wp), DIMENSION (jpi,jpj,jpk) ::   zti, ztu, ztv, ztw   ! 3D workspace
77      !!----------------------------------------------------------------------
78
79      zti(:,:,:) = 0.e0 
80
81      IF( kt == nit000 .AND. lwp ) THEN
82         WRITE(numout,*)
83         WRITE(numout,*) 'tra_adv_tvd : TVD advection scheme'
84         WRITE(numout,*) '~~~~~~~~~~~'
85      ENDIF
86
87      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;    z2 = 1.      ! euler   time-stepping
88      ELSE                                        ;    z2 = 2.      ! leap-frog time-stepping
89      ENDIF
90
91      ! 1. Bottom value : flux set to zero
92      ! ---------------
93      ztu(:,:,jpk) = 0.e0   ;   ztv(:,:,jpk) = 0.e0
94      ztw(:,:,jpk) = 0.e0   ;   zti(:,:,jpk) = 0.e0
95
96
97      ! 2. upstream advection with initial mass fluxes & intermediate update
98      ! --------------------------------------------------------------------
99      ! upstream tracer flux in the i and j direction
100      DO jk = 1, jpkm1
101         DO jj = 1, jpjm1
102            DO ji = 1, fs_jpim1   ! vector opt.
103               zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) )
104               zfm_ui = pun(ji,jj,jk) - ABS( pun(ji,jj,jk) )
105               zfp_vj = pvn(ji,jj,jk) + ABS( pvn(ji,jj,jk) )
106               zfm_vj = pvn(ji,jj,jk) - ABS( pvn(ji,jj,jk) )
107               ztu(ji,jj,jk) = 0.5 * ( zfp_ui * ptb(ji,jj,jk) + zfm_ui * ptb(ji+1,jj  ,jk) )
108               ztv(ji,jj,jk) = 0.5 * ( zfp_vj * ptb(ji,jj,jk) + zfm_vj * ptb(ji  ,jj+1,jk) )
109            END DO
110         END DO
111      END DO
112      !
113      ! upstream tracer flux in the k direction
114      !                                                                           ! Surface value
115      IF( lk_dynspg_rl .OR. lk_vvl ) THEN   ;   ztw(:,:,1) = 0.e0                      ! rigid lid or non-linear fs
116      ELSE                                  ;   ztw(:,:,1) = pwn(:,:,1) * ptb(:,:,1)   ! free surface
117      ENDIF
118      !
119      DO jk = 2, jpkm1                                                            ! Interior value
120         DO jj = 1, jpj
121            DO ji = 1, jpi
122               zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) )
123               zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) )
124               ztw(ji,jj,jk) = 0.5 * ( zfp_wk * ptb(ji,jj,jk) + zfm_wk * ptb(ji,jj,jk-1) )
125            END DO
126         END DO
127      END DO
128      !
129      ! upstream advective trend
130      DO jk = 1, jpkm1
131         z2dtt = z2 * rdttra(jk)
132         DO jj = 2, jpjm1
133            DO ji = fs_2, fs_jpim1   ! vector opt.
134               zbtr = 1./ ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
135               ! total intermediate advective trends
136               zzti = - ( ztu(ji,jj,jk) - ztu(ji-1,jj  ,jk  )            &
137                  &     + ztv(ji,jj,jk) - ztv(ji  ,jj-1,jk  )            &
138                  &     + ztw(ji,jj,jk) - ztw(ji  ,jj  ,jk+1) ) * zbtr
139               ! update and guess with monotonic sheme
140               pta(ji,jj,jk) =   pta(ji,jj,jk) +         zzti
141               zti(ji,jj,jk) = ( ptb(ji,jj,jk) + z2dtt * zzti ) * tmask(ji,jj,jk)
142            END DO
143         END DO
144      END DO
145      !
146      CALL lbc_lnk( zti, 'T', 1. )      ! Lateral boundary conditions on zti   (unchanged sign)
147
148
149      !                                 ! trend diagnostics (contribution of upstream fluxes)
150      IF( l_trdtra ) THEN
151         CALL trd_tra_adv( kt, ktra, jpt_trd_xad, cdtype, ztu, pun, ptn )
152         CALL trd_tra_adv( kt, ktra, jpt_trd_yad, cdtype, ztv, pvn, ptn )
153         CALL trd_tra_adv( kt, ktra, jpt_trd_zad, cdtype, ztw, pwn, ptn )
154      ENDIF
155      !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
156      IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN
157         IF( ktra == jp_tem)   pht_adv(:) = ptr_vj( ztv(:,:,:) ) 
158         IF( ktra == jp_sal)   pst_adv(:) = ptr_vj( ztv(:,:,:) )
159      ENDIF
160
161
162      ! 3. antidiffusive flux : high order minus low order
163      ! --------------------------------------------------
164      !                                 ! anti-diffusive flux on i and j
165      DO jk = 1, jpkm1
166         DO jj = 1, jpjm1
167            DO ji = 1, fs_jpim1   ! vector opt.
168               ztu(ji,jj,jk) = 0.5 * pun(ji,jj,jk) * ( ptn(ji,jj,jk) + ptn(ji+1,jj,jk) ) - ztu(ji,jj,jk)
169               ztv(ji,jj,jk) = 0.5 * pvn(ji,jj,jk) * ( ptn(ji,jj,jk) + ptn(ji,jj+1,jk) ) - ztv(ji,jj,jk)
170            END DO
171         END DO
172      END DO
173      !                                 ! antidiffusive flux on k
174      ztw(:,:,1) = 0.e0                      ! Surface value
175      !
176      DO jk = 2, jpkm1                       ! Interior value
177         DO jj = 1, jpj
178            DO ji = 1, jpi
179               ztw(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk) + ptn(ji,jj,jk-1) ) - ztw(ji,jj,jk)
180            END DO
181         END DO
182      END DO
183      !
184      CALL lbc_lnk( ztu, 'U', -1. )     ! Lateral bondary conditions on upstream fluxes
185      CALL lbc_lnk( ztv, 'V', -1. )
186      CALL lbc_lnk( ztw, 'W',  1. )
187
188      !                                 ! monotonicity algorithm
189      CALL nonosc( ptb, ztu, ztv, ztw, zti, z2 )
190
191
192      ! 4. final trend with anti-diffusive fluxes
193      ! -----------------------------------------
194      DO jk = 1, jpkm1
195         DO jj = 2, jpjm1
196            DO ji = fs_2, fs_jpim1   ! vector opt. 
197               zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
198               ! anti-diffusive trends added to the general tracer trends
199               pta(ji,jj,jk) = pta(ji,jj,jk) - ( ztu(ji,jj,jk) - ztu(ji-1,jj  ,jk  )             &
200                  &                            + ztv(ji,jj,jk) - ztv(ji  ,jj-1,jk  )             &
201                  &                            + ztw(ji,jj,jk) - ztw(ji  ,jj  ,jk+1)  ) * zbtr
202            END DO
203         END DO
204      END DO
205
206      IF( l_trdtra ) THEN               ! trend diagnostic (contribution of anti-diffusive fluxes)
207         CALL trd_tra_adv( kt, ktra, jpt_trd_xad, cdtype, ztu, pun, ptn, cnbpas='bis' )   ! <<< Add to iad trend
208         CALL trd_tra_adv( kt, ktra, jpt_trd_yad, cdtype, ztv, pvn, ptn, cnbpas='bis' )   ! <<< Add to jad trend
209         CALL trd_tra_adv( kt, ktra, jpt_trd_zad, cdtype, ztw, pwn, ptn, cnbpas='bis' )   ! <<< Add to zad trend
210      ENDIF
211      !                                 ! "Poleward" heat and salt transports (contribution of anti-diffusive fluxes)
212      IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN
213         IF( ktra == jp_tem)   pht_adv(:) = pht_adv(:) + ptr_vj( ztv(:,:,:) )
214         IF( ktra == jp_sal)   pst_adv(:) = pst_adv(:) + ptr_vj( ztv(:,:,:) )
215      ENDIF
216      !                                 ! control print
217      IF(ln_ctl)   CALL prt_ctl( tab3d_1=pta, clinfo1=' tvd - adv: ', mask1=tmask, clinfo3=cdtype )
218      !
219   END SUBROUTINE tra_adv_tvd
220
221
222   SUBROUTINE nonosc( pbef, paa, pbb, pcc, paft, prdt )
223      !!---------------------------------------------------------------------
224      !!                    ***  ROUTINE nonosc  ***
225      !!     
226      !! **  Purpose :   compute monotonic tracer fluxes from the upstream
227      !!       scheme and the before field by a nonoscillatory algorithm
228      !!
229      !! **  Method  :   ... ???
230      !!       warning : pbef and paft must be masked, but the boundaries
231      !!       conditions on the fluxes are not necessary zalezak (1979)
232      !!       drange (1995) multi-dimensional forward-in-time and upstream-
233      !!       in-space based differencing for fluid
234      !!----------------------------------------------------------------------
235      REAL(wp), INTENT(in   ) ::   prdt                               ! ???
236      REAL(wp), INTENT(inout), DIMENSION (jpi,jpj,jpk) ::   pbef, paft       ! before & after field
237      REAL(wp), INTENT(inout), DIMENSION (jpi,jpj,jpk) ::   paa, pbb, pcc    ! monotonic flux in the 3 directions
238      !!
239      INTEGER ::   ji, jj, jk               ! dummy loop indices
240      INTEGER ::   ikm1
241      REAL(wp), DIMENSION (jpi,jpj,jpk) ::   zbetup, zbetdo
242      REAL(wp) ::   zpos, zneg, zbt, za, zb, zc, zbig, zrtrn, z2dtt
243      !!----------------------------------------------------------------------
244
245      zbig = 1.e+40
246      zrtrn = 1.e-15
247      zbetup(:,:,:) = 0.e0   ;   zbetdo(:,:,:) = 0.e0 
248
249!!gm   optimisation :  add the optimal version I wrote 1 year ago
250      ! Search local extrema
251      ! --------------------
252      ! large negative value (-zbig) inside land
253      pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:) - zbig * ( 1.e0 - tmask(:,:,:) )
254      paft(:,:,:) = paft(:,:,:) * tmask(:,:,:) - zbig * ( 1.e0 - tmask(:,:,:) )
255      ! search maximum in neighbourhood
256      DO jk = 1, jpkm1
257         ikm1 = MAX(jk-1,1)
258         DO jj = 2, jpjm1
259            DO ji = fs_2, fs_jpim1   ! vector opt.
260               zbetup(ji,jj,jk) = MAX(  pbef(ji  ,jj  ,jk  ), paft(ji  ,jj  ,jk  ),   &
261                  &                     pbef(ji-1,jj  ,jk  ), pbef(ji+1,jj  ,jk  ),   &
262                  &                     paft(ji-1,jj  ,jk  ), paft(ji+1,jj  ,jk  ),   &
263                  &                     pbef(ji  ,jj-1,jk  ), pbef(ji  ,jj+1,jk  ),   &
264                  &                     paft(ji  ,jj-1,jk  ), paft(ji  ,jj+1,jk  ),   &
265                  &                     pbef(ji  ,jj  ,ikm1), pbef(ji  ,jj  ,jk+1),   &
266                  &                     paft(ji  ,jj  ,ikm1), paft(ji  ,jj  ,jk+1)  )
267            END DO
268         END DO
269      END DO
270      ! large positive value (+zbig) inside land
271      pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:) + zbig * ( 1.e0 - tmask(:,:,:) )
272      paft(:,:,:) = paft(:,:,:) * tmask(:,:,:) + zbig * ( 1.e0 - tmask(:,:,:) )
273      ! search minimum in neighbourhood
274      DO jk = 1, jpkm1
275         ikm1 = MAX(jk-1,1)
276         DO jj = 2, jpjm1
277            DO ji = fs_2, fs_jpim1   ! vector opt.
278               zbetdo(ji,jj,jk) = MIN(  pbef(ji  ,jj  ,jk  ), paft(ji  ,jj  ,jk  ),   &
279                  &                     pbef(ji-1,jj  ,jk  ), pbef(ji+1,jj  ,jk  ),   &
280                  &                     paft(ji-1,jj  ,jk  ), paft(ji+1,jj  ,jk  ),   &
281                  &                     pbef(ji  ,jj-1,jk  ), pbef(ji  ,jj+1,jk  ),   &
282                  &                     paft(ji  ,jj-1,jk  ), paft(ji  ,jj+1,jk  ),   &
283                  &                     pbef(ji  ,jj  ,ikm1), pbef(ji  ,jj  ,jk+1),   &
284                  &                     paft(ji  ,jj  ,ikm1), paft(ji  ,jj  ,jk+1)  )
285            END DO
286         END DO
287      END DO
288
289      ! restore masked values to zero
290      pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:)
291      paft(:,:,:) = paft(:,:,:) * tmask(:,:,:)
292 
293
294      ! 2. Positive and negative part of fluxes and beta terms
295      ! ------------------------------------------------------
296
297      DO jk = 1, jpkm1
298         z2dtt = prdt * rdttra(jk)
299         DO jj = 2, jpjm1
300            DO ji = fs_2, fs_jpim1   ! vector opt.
301               ! positive & negative part of the flux
302               zpos = MAX( 0., paa(ji-1,jj  ,jk  ) ) - MIN( 0., paa(ji  ,jj  ,jk  ) )   &
303                  & + MAX( 0., pbb(ji  ,jj-1,jk  ) ) - MIN( 0., pbb(ji  ,jj  ,jk  ) )   &
304                  & + MAX( 0., pcc(ji  ,jj  ,jk+1) ) - MIN( 0., pcc(ji  ,jj  ,jk  ) )
305               zneg = MAX( 0., paa(ji  ,jj  ,jk  ) ) - MIN( 0., paa(ji-1,jj  ,jk  ) )   &
306                  & + MAX( 0., pbb(ji  ,jj  ,jk  ) ) - MIN( 0., pbb(ji  ,jj-1,jk  ) )   &
307                  & + MAX( 0., pcc(ji  ,jj  ,jk  ) ) - MIN( 0., pcc(ji  ,jj  ,jk+1) )
308               ! up & down beta terms
309               zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) / z2dtt
310               zbetup(ji,jj,jk) = ( zbetup(ji,jj,jk) - paft(ji,jj,jk) ) / (zpos+zrtrn) * zbt
311               zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zbetdo(ji,jj,jk) ) / (zneg+zrtrn) * zbt
312            END DO
313         END DO
314      END DO
315
316      ! lateral boundary condition on zbetup & zbetdo   (unchanged sign)
317      CALL lbc_lnk( zbetup, 'T', 1. )
318      CALL lbc_lnk( zbetdo, 'T', 1. )
319
320
321      ! 3. monotonic flux in the i & j direction (paa & pbb)
322      ! ----------------------------------------
323      DO jk = 1, jpkm1
324         DO jj = 2, jpjm1
325            DO ji = fs_2, fs_jpim1   ! vector opt.
326               za = MIN( 1.e0, zbetdo(ji,jj,jk), zbetup(ji+1,jj,jk) )
327               zb = MIN( 1.e0, zbetup(ji,jj,jk), zbetdo(ji+1,jj,jk) )
328               zc = 0.5 * ( 1.e0 + SIGN( 1.e0, paa(ji,jj,jk) ) )
329               paa(ji,jj,jk) = paa(ji,jj,jk) * ( zc * za + ( 1.e0 - zc) * zb )
330
331               za = MIN( 1.e0, zbetdo(ji,jj,jk), zbetup(ji,jj+1,jk) )
332               zb = MIN( 1.e0, zbetup(ji,jj,jk), zbetdo(ji,jj+1,jk) )
333               zc = 0.5 * ( 1.e0 + SIGN( 1.e0, pbb(ji,jj,jk) ) )
334               pbb(ji,jj,jk) = pbb(ji,jj,jk) * ( zc * za + ( 1.e0 - zc) * zb )
335            END DO
336         END DO
337      END DO
338
339
340      ! monotonic flux in the k direction, i.e. pcc
341      ! -------------------------------------------
342      DO jk = 2, jpkm1
343         DO jj = 2, jpjm1
344            DO ji = fs_2, fs_jpim1   ! vector opt.
345
346               za = MIN( 1., zbetdo(ji,jj,jk), zbetup(ji,jj,jk-1) )
347               zb = MIN( 1., zbetup(ji,jj,jk), zbetdo(ji,jj,jk-1) )
348               zc = 0.5 * ( 1.e0 + SIGN( 1.e0, pcc(ji,jj,jk) ) )
349               pcc(ji,jj,jk) = pcc(ji,jj,jk) * ( zc * za + ( 1.e0 - zc) * zb )
350            END DO
351         END DO
352      END DO
353
354      ! lateral boundary condition on paa, pbb, pcc
355      CALL lbc_lnk( paa, 'U', -1. )      ! changed sign
356      CALL lbc_lnk( pbb, 'V', -1. )      ! changed sign
357      CALL lbc_lnk( pcc, 'W',  1. )      ! NO changed sign
358      !
359   END SUBROUTINE nonosc
360
361   !!======================================================================
362END MODULE traadv_tvd
Note: See TracBrowser for help on using the repository browser.