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_ubs.F90 in tags/nemo_v3_2/nemo_v3_2/NEMO/OPA_SRC/TRA – NEMO

source: tags/nemo_v3_2/nemo_v3_2/NEMO/OPA_SRC/TRA/traadv_ubs.F90 @ 1878

Last change on this file since 1878 was 1878, checked in by flavoni, 14 years ago

initial test for nemogcm
File size: 25.5 KB
RevLine 
[1878]1MODULE traadv_ubs
2   !!==============================================================================
3   !!                       ***  MODULE  traadv_ubs  ***
4   !! Ocean active tracers:  horizontal & vertical advective trend
5   !!==============================================================================
6   !! History :  9.0  !  06-08  (L. Debreu, R. Benshila)  Original code
7   !!----------------------------------------------------------------------
8
9   !!----------------------------------------------------------------------
10   !!   tra_adv_ubs : update the tracer trend with the horizontal
11   !!                 advection trends using a third order biaised scheme 
12   !!----------------------------------------------------------------------
13   USE oce             ! ocean dynamics and active tracers
14   USE dom_oce         ! ocean space and time domain
15   USE trdmod
16   USE trdmod_oce
17   USE lib_mpp
18   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
19   USE in_out_manager  ! I/O manager
20   USE diaptr          ! poleward transport diagnostics
21   USE dynspg_oce      ! choice/control of key cpp for surface pressure gradient
22   USE prtctl
23
24   IMPLICIT NONE
25   PRIVATE
26
27   PUBLIC   tra_adv_ubs   ! routine called by traadv module
28
29   REAL(wp), DIMENSION(jpi,jpj) ::   e1e2tr   ! = 1/(e1t * e2t)
30
31   !! * Substitutions
32#  include "domzgr_substitute.h90"
33#  include "vectopt_loop_substitute.h90"
34   !!----------------------------------------------------------------------
35   !!   OPA 9.0 , LOCEAN-IPSL (2006) 
36   !! $Id: traadv_ubs.F90 1528 2009-07-23 14:38:47Z rblod $
37   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
38   !!----------------------------------------------------------------------
39
40CONTAINS
41
42   SUBROUTINE tra_adv_ubs( kt, pun, pvn, pwn )
43      !!----------------------------------------------------------------------
44      !!                  ***  ROUTINE tra_adv_ubs  ***
45      !!                 
46      !! ** Purpose :   Compute the now trend due to the advection of tracers
47      !!      and add it to the general trend of passive tracer equations.
48      !!
49      !! ** Method  :   The upstream biased third (UBS) is order scheme based
50      !!      on an upstream-biased parabolic interpolation (Shchepetkin and McWilliams 2005)
51      !!      It is only used in the horizontal direction.
52      !!      For example the i-component of the advective fluxes are given by :
53      !!                !  e1u e3u un ( mi(Tn) - zltu(i  ) )   if un(i) >= 0
54      !!          zwx = !  or
55      !!                !  e1u e3u un ( mi(Tn) - zltu(i+1) )   if un(i) < 0
56      !!      where zltu is the second derivative of the before temperature field:
57      !!          zltu = 1/e3t di[ e2u e3u / e1u di[Tb] ]
58      !!      This results in a dissipatively dominant (i.e. hyper-diffusive)
59      !!      truncation error. The overall performance of the advection scheme
60      !!      is similar to that reported in (Farrow and Stevens, 1995).
61      !!      For stability reasons, the first term of the fluxes which corresponds
62      !!      to a second order centered scheme is evaluated using the now velocity
63      !!      (centered in time) while the second term which is the diffusive part
64      !!      of the scheme, is evaluated using the before velocity (forward in time).
65      !!      Note that UBS is not positive. Do not use it on passive tracers.
66      !!      On the vertical, the advection is evaluated using a TVD scheme, as
67      !!      the UBS have been found to be too diffusive.
68      !!
69      !! ** Action : - update (ta,sa) with the now advective tracer trends
70      !!
71      !! Reference : Shchepetkin, A. F., J. C. McWilliams, 2005, Ocean Modelling, 9, 347-404.
72      !!             Farrow, D.E., Stevens, D.P., 1995, J. Phys. Ocean. 25, 1731Ð1741.
73      !!----------------------------------------------------------------------
74      USE oce, ONLY :   zwx => ua   ! use ua as workspace
75      USE oce, ONLY :   zwy => va   ! use va as workspace
76      !!
77      INTEGER , INTENT(in)                         ::   kt             ! ocean time-step index
78      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::  pun   ! effective ocean velocity, u_component
79      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::  pvn   ! effective ocean velocity, v_component
80      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::  pwn   ! effective ocean velocity, w_component
81      !!
82      INTEGER  ::   ji, jj, jk                 ! dummy loop indices
83      REAL(wp) ::   zta, zsa, zbtr, zcoef                  ! temporary scalars
84      REAL(wp) ::   zfui, zfp_ui, zfm_ui, zcenut, zcenus   !    "         "
85      REAL(wp) ::   zfvj, zfp_vj, zfm_vj, zcenvt, zcenvs   !    "         "
86      REAL(wp) ::   z_hdivn_x, z_hdivn_y, z_hdivn          !    "         "
87      REAL(wp), DIMENSION(jpi,jpj)     ::   zeeu, zeev     ! temporary 2D workspace
88      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwz , zww                        ! temporary 3D workspace
89      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ztu , ztv , zltu , zltv, ztrdt   !    "              "
90      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zsu , zsv , zlsu , zlsv, ztrds   !    "              "
91      !!----------------------------------------------------------------------
92
93      zltu(:,:,:) = 0.e0
94      zltv(:,:,:) = 0.e0
95      zlsu(:,:,:) = 0.e0
96      zlsv(:,:,:) = 0.e0
97
98      IF( kt == nit000 ) THEN
99         IF(lwp) WRITE(numout,*)
100         IF(lwp) WRITE(numout,*) 'tra_adv_ubs :  horizontal UBS advection scheme'
101         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
102         !
103         e1e2tr(:,:) = 1. / ( e1t(:,:) * e2t(:,:) )
104      ENDIF
105
106      ! Save ta and sa trends
107      ztrdt(:,:,:) = ta(:,:,:)
108      ztrds(:,:,:) = sa(:,:,:)
109
110      zcoef = 1./6.
111      !                                                ! ===============
112      DO jk = 1, jpkm1                                 ! Horizontal slab
113         !                                             ! ===============
114
115         !  Initialization of metric arrays (for z- or s-coordinates)
116         DO jj = 1, jpjm1
117            DO ji = 1, fs_jpim1   ! vector opt.
118#if defined key_zco
119               ! z-coordinates, no vertical scale factors
120               zeeu(ji,jj) = e2u(ji,jj) / e1u(ji,jj) * umask(ji,jj,jk)
121               zeev(ji,jj) = e1v(ji,jj) / e2v(ji,jj) * vmask(ji,jj,jk)
122#else
123               ! s-coordinates, vertical scale factor are used
124               zeeu(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj) * umask(ji,jj,jk)
125               zeev(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj) * vmask(ji,jj,jk)
126#endif
127            END DO
128         END DO
129
130         !  Laplacian
131         ! First derivative (gradient)
132         DO jj = 1, jpjm1
133            DO ji = 1, fs_jpim1   ! vector opt.
134               ztu(ji,jj,jk) = zeeu(ji,jj) * ( tb(ji+1,jj  ,jk) - tb(ji,jj,jk) )
135               zsu(ji,jj,jk) = zeeu(ji,jj) * ( sb(ji+1,jj  ,jk) - sb(ji,jj,jk) )
136               ztv(ji,jj,jk) = zeev(ji,jj) * ( tb(ji  ,jj+1,jk) - tb(ji,jj,jk) )
137               zsv(ji,jj,jk) = zeev(ji,jj) * ( sb(ji  ,jj+1,jk) - sb(ji,jj,jk) )
138            END DO
139         END DO
140         ! Second derivative (divergence)
141         DO jj = 2, jpjm1
142            DO ji = fs_2, fs_jpim1   ! vector opt.
143#if ! defined key_zco
144               zcoef = 1. / ( 6. * fse3t(ji,jj,jk) )
145#endif         
146               zltu(ji,jj,jk) = (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk)  ) * zcoef
147               zlsu(ji,jj,jk) = (  zsu(ji,jj,jk) - zsu(ji-1,jj,jk)  ) * zcoef
148               zltv(ji,jj,jk) = (  ztv(ji,jj,jk) - ztv(ji,jj-1,jk)  ) * zcoef
149               zlsv(ji,jj,jk) = (  zsv(ji,jj,jk) - zsv(ji,jj-1,jk)  ) * zcoef
150            END DO
151         END DO
152         !                                             ! =================
153      END DO                                           !    End of slab
154      !                                                ! =================
155
156      ! Lateral boundary conditions on the laplacian (zlt,zls)   (unchanged sgn)
157      CALL lbc_lnk( zltu, 'T', 1. )   ;    CALL lbc_lnk( zlsu, 'T', 1. )
158      CALL lbc_lnk( zltv, 'T', 1. )   ;    CALL lbc_lnk( zlsv, 'T', 1. )
159
160      !                                                ! ===============
161      DO jk = 1, jpkm1                                 ! Horizontal slab
162         !                                             ! ===============
163         !  Horizontal advective fluxes
164         DO jj = 1, jpjm1
165            DO ji = 1, fs_jpim1   ! vector opt.
166               ! volume fluxes * 1/2
167#if defined key_zco
168               zfui = 0.5 * e2u(ji,jj) * pun(ji,jj,jk)
169               zfvj = 0.5 * e1v(ji,jj) * pvn(ji,jj,jk)
170#else
171               zfui = 0.5 * e2u(ji,jj) * fse3u(ji,jj,jk) * pun(ji,jj,jk)
172               zfvj = 0.5 * e1v(ji,jj) * fse3v(ji,jj,jk) * pvn(ji,jj,jk)
173#endif
174               ! upstream scheme
175               zfp_ui = zfui + ABS( zfui )
176               zfp_vj = zfvj + ABS( zfvj )
177               zfm_ui = zfui - ABS( zfui )
178               zfm_vj = zfvj - ABS( zfvj )
179               ! centered scheme
180               zcenut = zfui * ( tn(ji,jj,jk) + tn(ji+1,jj  ,jk) )
181               zcenvt = zfvj * ( tn(ji,jj,jk) + tn(ji  ,jj+1,jk) )
182               zcenus = zfui * ( sn(ji,jj,jk) + sn(ji+1,jj  ,jk) )
183               zcenvs = zfvj * ( sn(ji,jj,jk) + sn(ji  ,jj+1,jk) )
184               ! mixed centered / upstream scheme
185               zwx(ji,jj,jk) = zcenut - zfp_ui * zltu(ji,jj,jk) -zfm_ui * zltu(ji+1,jj,jk)
186               zwy(ji,jj,jk) = zcenvt - zfp_vj * zltv(ji,jj,jk) -zfm_vj * zltv(ji,jj+1,jk)
187               zww(ji,jj,jk) = zcenus - zfp_ui * zlsu(ji,jj,jk) -zfm_ui * zlsu(ji+1,jj,jk)
188               zwz(ji,jj,jk) = zcenvs - zfp_vj * zlsv(ji,jj,jk) -zfm_vj * zlsv(ji,jj+1,jk)
189            END DO
190         END DO
191
192         !  Tracer flux divergence at t-point added to the general trend
193         DO jj = 2, jpjm1
194            DO ji = fs_2, fs_jpim1   ! vector opt.
195               ! horizontal advective trends
196#if defined key_zco
197               zbtr = e1e2tr(ji,jj)
198#else
199               zbtr = e1e2tr(ji,jj) / fse3t(ji,jj,jk)
200#endif
201               zta = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk)   &
202                  &            + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk)  )
203               zsa = - zbtr * (  zww(ji,jj,jk) - zww(ji-1,jj  ,jk)   &
204                  &            + zwz(ji,jj,jk) - zwz(ji  ,jj-1,jk)  )
205               ! add it to the general tracer trends
206               ta(ji,jj,jk) = ta(ji,jj,jk) + zta
207               sa(ji,jj,jk) = sa(ji,jj,jk) + zsa
208            END DO
209         END DO
210         !                                             ! ===============
211      END DO                                           !   End of slab
212      !                                                ! ===============
213
214      ! Horizontal trend used in tra_adv_ztvd subroutine
215      zltu(:,:,:) = ta(:,:,:) - ztrdt(:,:,:) 
216      zlsu(:,:,:) = sa(:,:,:) - ztrds(:,:,:) 
217
218      ! 3. Save the horizontal advective trends for diagnostic
219      ! ------------------------------------------------------
220      IF( l_trdtra )   THEN
221         ! Recompute the hoizontal advection zta & zsa trends computed
222         ! at the step 2. above in making the difference between the new
223         ! trends and the previous one ta()/sa - ztrdt()/ztrds() and add
224         ! the term tn()/sn()*hdivn() to recover the Uh gradh(T/S) trends
225         ztrdt(:,:,:) = 0.e0   ;   ztrds(:,:,:) = 0.e0
226         !
227         ! T/S ZONAL advection trends
228         DO jk = 1, jpkm1
229            DO jj = 2, jpjm1
230               DO ji = fs_2, fs_jpim1   ! vector opt.
231                  !-- Compute zonal divergence by splitting hdivn (see divcur.F90)
232#if defined key_zco
233                  zbtr = e1e2tr(ji,jj)
234                  z_hdivn_x = (  e2u(ji  ,jj) * pun(ji  ,jj,jk)          &
235                     &         - e2u(ji-1,jj) * pun(ji-1,jj,jk) ) * zbtr
236#else
237                  zbtr = e1e2tr(ji,jj) / fse3t(ji,jj,jk)
238                  z_hdivn_x = (  e2u(ji  ,jj) * fse3u(ji  ,jj,jk) * pun(ji  ,jj,jk)          &
239                     &         - e2u(ji-1,jj) * fse3u(ji-1,jj,jk) * pun(ji-1,jj,jk) ) * zbtr
240#endif
241                  ztrdt(ji,jj,jk) = - ( zwx(ji,jj,jk) - zwx(ji-1,jj,jk) ) * zbtr + tn(ji,jj,jk) * z_hdivn_x
242                  ztrds(ji,jj,jk) = - ( zww(ji,jj,jk) - zww(ji-1,jj,jk) ) * zbtr + sn(ji,jj,jk) * z_hdivn_x
243               END DO
244            END DO
245         END DO
246         CALL trd_mod(ztrdt, ztrds, jptra_trd_xad, 'TRA', kt)    ! save the trends
247         !
248         ! T/S MERIDIONAL advection trends
249         DO jk = 1, jpkm1
250            DO jj = 2, jpjm1
251               DO ji = fs_2, fs_jpim1   ! vector opt.
252                  !-- Compute merid. divergence by splitting hdivn (see divcur.F90)
253#if defined key_zco
254                  zbtr      = e1e2tr(ji,jj)
255                  z_hdivn_y = (  e1v(ji,  jj) * pvn(ji,jj  ,jk)          &
256                     &         - e1v(ji,jj-1) * pvn(ji,jj-1,jk) ) * zbtr
257#else
258                  zbtr      = e1e2tr(ji,jj) / fse3t(ji,jj,jk)
259                  z_hdivn_y = (  e1v(ji,  jj) * fse3v(ji,jj  ,jk) * pvn(ji,jj  ,jk)          &
260                     &         - e1v(ji,jj-1) * fse3v(ji,jj-1,jk) * pvn(ji,jj-1,jk) ) * zbtr
261#endif
262                  ztrdt(ji,jj,jk) = - ( zwy(ji,jj,jk) - zwy(ji,jj-1,jk) ) * zbtr + tn(ji,jj,jk) * z_hdivn_y         
263                  ztrds(ji,jj,jk) = - ( zwz(ji,jj,jk) - zwz(ji,jj-1,jk) ) * zbtr + sn(ji,jj,jk) * z_hdivn_y
264               END DO
265            END DO
266         END DO
267         CALL trd_mod(ztrdt, ztrds, jptra_trd_yad, 'TRA', kt)     ! save the trends
268         !
269      ENDIF
270
271      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ta, clinfo1=' ubs had  - Ta: ', mask1=tmask,   &
272         &                       tab3d_2=sa, clinfo2=           ' Sa: ', mask2=tmask, clinfo3='tra' )
273
274      ! "zonal" mean advective heat and salt transport
275      IF( ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN
276         IF( lk_zco ) THEN
277            DO jk = 1, jpkm1
278               DO jj = 2, jpjm1
279                  DO ji = fs_2, fs_jpim1   ! vector opt.
280                     zwy(ji,jj,jk) = zwy(ji,jj,jk) * fse3v(ji,jj,jk)
281                     zwz(ji,jj,jk) = zwz(ji,jj,jk) * fse3v(ji,jj,jk)
282                  END DO
283               END DO
284            END DO
285         ENDIF
286         pht_adv(:) = ptr_vj( zwy(:,:,:) )
287         pst_adv(:) = ptr_vj( zwz(:,:,:) )
288      ENDIF
289
290      ! II. Vertical advection
291      ! ----------------------
292      IF( l_trdtra ) THEN          ! Save ta and sa trends
293         ztrdt(:,:,:) = ta(:,:,:)
294         ztrds(:,:,:) = sa(:,:,:)
295      ENDIF
296   
297      ! TVD scheme the vertical direction 
298      CALL tra_adv_ztvd(kt, pwn, zltu, zlsu)
299
300      IF( l_trdtra )   THEN         !  Save the final vertical advective trends
301         DO jk = 1, jpkm1
302            DO jj = 2, jpjm1
303               DO ji = fs_2, fs_jpim1   ! vector opt.
304#if defined key_zco
305                  zbtr      = e1e2tr(ji,jj)
306                  z_hdivn_x = e2u(ji,jj)*pun(ji,jj,jk) - e2u(ji-1,jj)*pun(ji-1,jj,jk)
307                  z_hdivn_y = e1v(ji,jj)*pvn(ji,jj,jk) - e1v(ji,jj-1)*pvn(ji,jj-1,jk)
308#else
309                  zbtr      = e1e2tr(ji,jj) / fse3t(ji,jj,jk)
310                  z_hdivn_x = e2u(ji,jj)*fse3u(ji,jj,jk)*pun(ji,jj,jk) - e2u(ji-1,jj)*fse3u(ji-1,jj,jk)*pun(ji-1,jj,jk)
311                  z_hdivn_y = e1v(ji,jj)*fse3v(ji,jj,jk)*pvn(ji,jj,jk) - e1v(ji,jj-1)*fse3v(ji,jj-1,jk)*pvn(ji,jj-1,jk)
312#endif
313                  z_hdivn   = (z_hdivn_x + z_hdivn_y) * zbtr
314                  zbtr      = e1e2tr(ji,jj) / fse3t(ji,jj,jk)
315                  ztrdt(ji,jj,jk) = ta(ji,jj,jk) - ztrdt(ji,jj,jk) - tn(ji,jj,jk) * z_hdivn
316                  ztrds(ji,jj,jk) = sa(ji,jj,jk) - ztrds(ji,jj,jk) - sn(ji,jj,jk) * z_hdivn
317               END DO
318            END DO
319         END DO
320         CALL trd_mod(ztrdt, ztrds, jptra_trd_zad, 'TRA', kt)   ! <<< ADD TO PREVIOUSLY COMPUTED
321         !
322      ENDIF
323
324      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ta, clinfo1=' ubs zad  - Ta: ', mask1=tmask,   &
325         &                       tab3d_2=sa, clinfo2=           ' Sa: ', mask2=tmask, clinfo3='tra')
326      !
327   END SUBROUTINE tra_adv_ubs
328
329
330   SUBROUTINE tra_adv_ztvd( kt, pwn, zttrd, zstrd )
331      !!----------------------------------------------------------------------
332      !!                  ***  ROUTINE tra_adv_ztvd  ***
333      !!
334      !! **  Purpose :   Compute the now trend due to total advection of
335      !!       tracers and add it to the general trend of tracer equations
336      !!
337      !! **  Method  :   TVD scheme, i.e. 2nd order centered scheme with
338      !!       corrected flux (monotonic correction)
339      !!       note: - this advection scheme needs a leap-frog time scheme
340      !!
341      !! ** Action : - update (ta,sa) with the now advective tracer trends
342      !!             - save the trends in (ztrdt,ztrds) ('key_trdtra')
343      !!----------------------------------------------------------------------
344      INTEGER , INTENT(in)                         ::   kt             ! ocean time-step
345      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::   pwn            ! verical effective velocity
346      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::   zttrd, zstrd   ! lateral advective trends on T & S
347      !!
348      INTEGER  ::   ji, jj, jk              ! dummy loop indices
349      REAL(wp) ::   z2dtt, zbtr, zew, z2    ! temporary scalar 
350      REAL(wp) ::   ztak, zfp_wk            !    "         "
351      REAL(wp) ::   zsak, zfm_wk            !    "         "
352      REAL(wp), DIMENSION (jpi,jpj,jpk) ::   zti, ztw   ! temporary 3D workspace
353      REAL(wp), DIMENSION (jpi,jpj,jpk) ::   zsi, zsw   !    "              "
354      !!----------------------------------------------------------------------
355
356      IF( kt == nit000 .AND. lwp ) THEN
357         WRITE(numout,*)
358         WRITE(numout,*) 'tra_adv_ztvd : vertical TVD advection scheme'
359         WRITE(numout,*) '~~~~~~~~~~~~'
360      ENDIF
361
362      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;    z2 = 1.
363      ELSE                                        ;    z2 = 2.
364      ENDIF
365
366      !  Bottom value : flux set to zero
367      ! --------------
368      ztw(:,:,jpk) = 0.e0   ;   zsw(:,:,jpk) = 0.e0
369      zti  (:,:,:) = 0.e0   ;   zsi  (:,:,:) = 0.e0
370
371
372      !  upstream advection with initial mass fluxes & intermediate update
373      ! -------------------------------------------------------------------
374      ! Surface value
375      IF( lk_vvl ) THEN
376         ! variable volume : flux set to zero
377         ztw(:,:,1) = 0.e0
378         zsw(:,:,1) = 0.e0
379      ELSE
380         ! free surface-constant volume
381         DO jj = 1, jpj
382            DO ji = 1, jpi
383               zew = e1t(ji,jj) * e2t(ji,jj) * pwn(ji,jj,1)
384               ztw(ji,jj,1) = zew * tb(ji,jj,1)
385               zsw(ji,jj,1) = zew * sb(ji,jj,1)
386            END DO
387         END DO
388      ENDIF
389
390      ! Interior value
391      DO jk = 2, jpkm1
392         DO jj = 1, jpj
393            DO ji = 1, jpi
394               zew = 0.5 * e1t(ji,jj) * e2t(ji,jj) * pwn(ji,jj,jk)
395               zfp_wk = zew + ABS( zew )
396               zfm_wk = zew - ABS( zew )
397               ztw(ji,jj,jk) = zfp_wk * tb(ji,jj,jk) + zfm_wk * tb(ji,jj,jk-1)
398               zsw(ji,jj,jk) = zfp_wk * sb(ji,jj,jk) + zfm_wk * sb(ji,jj,jk-1)
399            END DO
400         END DO
401      END DO
402
403      ! update and guess with monotonic sheme
404      DO jk = 1, jpkm1
405         z2dtt = z2 * rdttra(jk)
406         DO jj = 2, jpjm1
407            DO ji = fs_2, fs_jpim1   ! vector opt.
408               zbtr = 1./ ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
409               ztak = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * zbtr
410               zsak = - ( zsw(ji,jj,jk) - zsw(ji,jj,jk+1) ) * zbtr
411               ta(ji,jj,jk) =  ta(ji,jj,jk) + ztak
412               sa(ji,jj,jk) =  sa(ji,jj,jk) + zsak 
413               zti (ji,jj,jk) = ( tb(ji,jj,jk) + z2dtt * ( ztak + zttrd(ji,jj,jk) ) ) * tmask(ji,jj,jk)
414               zsi (ji,jj,jk) = ( sb(ji,jj,jk) + z2dtt * ( zsak + zstrd(ji,jj,jk) ) ) * tmask(ji,jj,jk)
415            END DO
416         END DO
417      END DO
418
419      ! Lateral boundary conditions on zti, zsi   (unchanged sign)
420      CALL lbc_lnk( zti, 'T', 1. )
421      CALL lbc_lnk( zsi, 'T', 1. )
422
423
424      !  antidiffusive flux : high order minus low order
425      ! -------------------------------------------------     
426      ! Surface value
427      ztw(:,:,1) = 0.e0   ;   zsw(:,:,1) = 0.e0
428
429      ! Interior value
430      DO jk = 2, jpkm1
431         DO jj = 1, jpj
432            DO ji = 1, jpi
433               zew = 0.5 * e1t(ji,jj) * e2t(ji,jj) * pwn(ji,jj,jk)
434               ztw(ji,jj,jk) = zew * ( tn(ji,jj,jk) + tn(ji,jj,jk-1) ) - ztw(ji,jj,jk)
435               zsw(ji,jj,jk) = zew * ( sn(ji,jj,jk) + sn(ji,jj,jk-1) ) - zsw(ji,jj,jk)
436            END DO
437         END DO
438      END DO
439
440      !  monotonicity algorithm
441      ! ------------------------
442      CALL nonosc_z( tb, ztw, zti, z2 )
443      CALL nonosc_z( sb, zsw, zsi, z2 )
444
445
446      !  final trend with corrected fluxes
447      ! -----------------------------------
448      DO jk = 1, jpkm1
449         DO jj = 2, jpjm1
450            DO ji = fs_2, fs_jpim1   ! vector opt. 
451               zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
452               ! k- vertical advective trends
453               ztak = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * zbtr
454               zsak = - ( zsw(ji,jj,jk) - zsw(ji,jj,jk+1) ) * zbtr
455               ! add them to the general tracer trends
456               ta(ji,jj,jk) = ta(ji,jj,jk) + ztak
457               sa(ji,jj,jk) = sa(ji,jj,jk) + zsak
458            END DO
459         END DO
460      END DO
461      !
462   END SUBROUTINE tra_adv_ztvd
463
464
465   SUBROUTINE nonosc_z( pbef, pcc, paft, prdt )
466      !!---------------------------------------------------------------------
467      !!                    ***  ROUTINE nonosc_z  ***
468      !!     
469      !! **  Purpose :   compute monotonic tracer fluxes from the upstream
470      !!       scheme and the before field by a nonoscillatory algorithm
471      !!
472      !! **  Method  :   ... ???
473      !!       warning : pbef and paft must be masked, but the boundaries
474      !!       conditions on the fluxes are not necessary zalezak (1979)
475      !!       drange (1995) multi-dimensional forward-in-time and upstream-
476      !!       in-space based differencing for fluid
477      !!----------------------------------------------------------------------
478      REAL(wp), INTENT(in   )                          ::   prdt   ! ???
479      REAL(wp), INTENT(inout), DIMENSION (jpi,jpj,jpk) ::   pbef   ! before field
480      REAL(wp), INTENT(inout), DIMENSION (jpi,jpj,jpk) ::   paft   ! after field
481      REAL(wp), INTENT(inout), DIMENSION (jpi,jpj,jpk) ::   pcc    ! monotonic flux in the k direction
482      !!
483      INTEGER  ::   ji, jj, jk               ! dummy loop indices
484      INTEGER  ::   ikm1
485      REAL(wp) ::   zpos, zneg, zbt, za, zb, zc, zbig, zrtrn, z2dtt
486      REAL(wp), DIMENSION (jpi,jpj,jpk) ::   zbetup, zbetdo
487      !!----------------------------------------------------------------------
488
489      zbig = 1.e+40
490      zrtrn = 1.e-15
491      zbetup(:,:,:) = 0.e0   ;   zbetdo(:,:,:) = 0.e0
492
493      ! Search local extrema
494      ! --------------------
495      ! large negative value (-zbig) inside land
496      pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:) - zbig * ( 1.e0 - tmask(:,:,:) )
497      paft(:,:,:) = paft(:,:,:) * tmask(:,:,:) - zbig * ( 1.e0 - tmask(:,:,:) )
498      ! search maximum in neighbourhood
499      DO jk = 1, jpkm1
500         ikm1 = MAX(jk-1,1)
501         DO jj = 2, jpjm1
502            DO ji = fs_2, fs_jpim1   ! vector opt.
503               zbetup(ji,jj,jk) = MAX(  pbef(ji  ,jj  ,jk  ), paft(ji  ,jj  ,jk  ),   &
504                  &                     pbef(ji  ,jj  ,ikm1), pbef(ji  ,jj  ,jk+1),   &
505                  &                     paft(ji  ,jj  ,ikm1), paft(ji  ,jj  ,jk+1)  )
506            END DO
507         END DO
508      END DO
509      ! large positive value (+zbig) inside land
510      pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:) + zbig * ( 1.e0 - tmask(:,:,:) )
511      paft(:,:,:) = paft(:,:,:) * tmask(:,:,:) + zbig * ( 1.e0 - tmask(:,:,:) )
512      ! search minimum in neighbourhood
513      DO jk = 1, jpkm1
514         ikm1 = MAX(jk-1,1)
515         DO jj = 2, jpjm1
516            DO ji = fs_2, fs_jpim1   ! vector opt.
517               zbetdo(ji,jj,jk) = MIN(  pbef(ji  ,jj  ,jk  ), paft(ji  ,jj  ,jk  ),   &
518                  &                     pbef(ji  ,jj  ,ikm1), pbef(ji  ,jj  ,jk+1),   &
519                  &                     paft(ji  ,jj  ,ikm1), paft(ji  ,jj  ,jk+1)  )
520            END DO
521         END DO
522      END DO
523
524      ! restore masked values to zero
525      pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:)
526      paft(:,:,:) = paft(:,:,:) * tmask(:,:,:)
527 
528
529      ! 2. Positive and negative part of fluxes and beta terms
530      ! ------------------------------------------------------
531
532      DO jk = 1, jpkm1
533         z2dtt = prdt * rdttra(jk)
534         DO jj = 2, jpjm1
535            DO ji = fs_2, fs_jpim1   ! vector opt.
536               ! positive & negative part of the flux
537               zpos = MAX( 0., pcc(ji  ,jj  ,jk+1) ) - MIN( 0., pcc(ji  ,jj  ,jk  ) )
538               zneg = MAX( 0., pcc(ji  ,jj  ,jk  ) ) - MIN( 0., pcc(ji  ,jj  ,jk+1) )
539               ! up & down beta terms
540               zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) / z2dtt
541               zbetup(ji,jj,jk) = ( zbetup(ji,jj,jk) - paft(ji,jj,jk) ) / (zpos+zrtrn) * zbt
542               zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zbetdo(ji,jj,jk) ) / (zneg+zrtrn) * zbt
543            END DO
544         END DO
545      END DO
546
547      ! monotonic flux in the k direction, i.e. pcc
548      ! -------------------------------------------
549      DO jk = 2, jpkm1
550         DO jj = 2, jpjm1
551            DO ji = fs_2, fs_jpim1   ! vector opt.
552
553               za = MIN( 1., zbetdo(ji,jj,jk), zbetup(ji,jj,jk-1) )
554               zb = MIN( 1., zbetup(ji,jj,jk), zbetdo(ji,jj,jk-1) )
555               zc = 0.5 * ( 1.e0 + SIGN( 1.e0, pcc(ji,jj,jk) ) )
556               pcc(ji,jj,jk) = pcc(ji,jj,jk) * ( zc * za + ( 1.e0 - zc) * zb )
557            END DO
558         END DO
559      END DO
560      !
561   END SUBROUTINE nonosc_z
562
563   !!======================================================================
564END MODULE traadv_ubs
Note: See TracBrowser for help on using the repository browser.