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

source: trunk/NEMO/OPA_SRC/TRA/traadv_ubs.F90 @ 507

Last change on this file since 507 was 503, checked in by opalod, 18 years ago

nemo_v1_update_064 : CT : general trends update including the addition of mean windows analysis possibility in the mixed layer
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 23.9 KB
Line 
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   !! $Header$
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  :   ???
50      !!
51      !! ** Action : - update (ta,sa) with the now advective tracer trends
52      !!----------------------------------------------------------------------
53      USE oce, ONLY :   zwx => ua   ! use ua as workspace
54      USE oce, ONLY :   zwy => va   ! use va as workspace
55      !!
56      INTEGER , INTENT(in)                         ::   kt             ! ocean time-step index
57      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::  pun   ! effective ocean velocity, u_component
58      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::  pvn   ! effective ocean velocity, v_component
59      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::  pwn   ! effective ocean velocity, w_component
60      !!
61      INTEGER  ::   ji, jj, jk                 ! dummy loop indices
62      REAL(wp) ::   zta, zsa, zbtr, zcoef                  ! temporary scalars
63      REAL(wp) ::   zfui, zfp_ui, zfm_ui, zcenut, zcenus   !    "         "
64      REAL(wp) ::   zfvj, zfp_vj, zfm_vj, zcenvt, zcenvs   !    "         "
65      REAL(wp) ::   z_hdivn_x, z_hdivn_y, z_hdivn          !    "         "
66      REAL(wp), DIMENSION(jpi,jpj)     ::   zeeu, zeev     ! temporary 2D workspace
67      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwz , zww                        ! temporary 3D workspace
68      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ztu , ztv , zltu , zltv, ztrdt   !    "              "
69      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zsu , zsv , zlsu , zlsv, ztrds   !    "              "
70      !!----------------------------------------------------------------------
71
72      zltu(:,:,:) = 0.e0
73      zltv(:,:,:) = 0.e0
74      zlsu(:,:,:) = 0.e0
75      zlsv(:,:,:) = 0.e0
76
77      IF( kt == nit000 ) THEN
78         IF(lwp) WRITE(numout,*)
79         IF(lwp) WRITE(numout,*) 'tra_adv_ubs :  horizontal UBS advection scheme'
80         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
81         !
82         e1e2tr(:,:) = 1. / ( e1t(:,:) * e2t(:,:) )
83      ENDIF
84
85      ! Save ta and sa trends
86      ztrdt(:,:,:) = ta(:,:,:)
87      ztrds(:,:,:) = sa(:,:,:)
88
89      zcoef = 1./6.
90      !                                                ! ===============
91      DO jk = 1, jpkm1                                 ! Horizontal slab
92         !                                             ! ===============
93
94         !  Initialization of metric arrays (for z- or s-coordinates)
95         DO jj = 1, jpjm1
96            DO ji = 1, fs_jpim1   ! vector opt.
97#if defined key_zco
98               ! z-coordinates, no vertical scale factors
99               zeeu(ji,jj) = e2u(ji,jj) / e1u(ji,jj) * umask(ji,jj,jk)
100               zeev(ji,jj) = e1v(ji,jj) / e2v(ji,jj) * vmask(ji,jj,jk)
101#else
102               ! s-coordinates, vertical scale factor are used
103               zeeu(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj) * umask(ji,jj,jk)
104               zeev(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj) * vmask(ji,jj,jk)
105#endif
106            END DO
107         END DO
108
109         !  Laplacian
110         ! First derivative (gradient)
111         DO jj = 1, jpjm1
112            DO ji = 1, fs_jpim1   ! vector opt.
113               ztu(ji,jj,jk) = zeeu(ji,jj) * ( tb(ji+1,jj  ,jk) - tb(ji,jj,jk) )
114               zsu(ji,jj,jk) = zeeu(ji,jj) * ( sb(ji+1,jj  ,jk) - sb(ji,jj,jk) )
115               ztv(ji,jj,jk) = zeev(ji,jj) * ( tb(ji  ,jj+1,jk) - tb(ji,jj,jk) )
116               zsv(ji,jj,jk) = zeev(ji,jj) * ( sb(ji  ,jj+1,jk) - sb(ji,jj,jk) )
117            END DO
118         END DO
119         ! Second derivative (divergence)
120         DO jj = 2, jpjm1
121            DO ji = fs_2, fs_jpim1   ! vector opt.
122#if ! defined key_zco
123               zcoef = 1. / ( 6. * fse3t(ji,jj,jk) )
124#endif         
125               zltu(ji,jj,jk) = (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk)  ) * zcoef
126               zlsu(ji,jj,jk) = (  zsu(ji,jj,jk) - zsu(ji-1,jj,jk)  ) * zcoef
127               zltv(ji,jj,jk) = (  ztv(ji,jj,jk) - ztv(ji,jj-1,jk)  ) * zcoef
128               zlsv(ji,jj,jk) = (  zsv(ji,jj,jk) - zsv(ji,jj-1,jk)  ) * zcoef
129            END DO
130         END DO
131         !                                             ! =================
132      END DO                                           !    End of slab
133      !                                                ! =================
134
135      ! Lateral boundary conditions on the laplacian (zlt,zls)   (unchanged sgn)
136      CALL lbc_lnk( zltu, 'T', 1. )   ;    CALL lbc_lnk( zlsu, 'T', 1. )
137      CALL lbc_lnk( zltv, 'T', 1. )   ;    CALL lbc_lnk( zlsv, 'T', 1. )
138
139      !                                                ! ===============
140      DO jk = 1, jpkm1                                 ! Horizontal slab
141         !                                             ! ===============
142         !  Horizontal advective fluxes
143         DO jj = 1, jpjm1
144            DO ji = 1, fs_jpim1   ! vector opt.
145               ! volume fluxes * 1/2
146#if defined key_zco
147               zfui = 0.5 * e2u(ji,jj) * pun(ji,jj,jk)
148               zfvj = 0.5 * e1v(ji,jj) * pvn(ji,jj,jk)
149#else
150               zfui = 0.5 * e2u(ji,jj) * fse3u(ji,jj,jk) * pun(ji,jj,jk)
151               zfvj = 0.5 * e1v(ji,jj) * fse3v(ji,jj,jk) * pvn(ji,jj,jk)
152#endif
153               ! upstream scheme
154               zfp_ui = zfui + ABS( zfui )
155               zfp_vj = zfvj + ABS( zfvj )
156               zfm_ui = zfui - ABS( zfui )
157               zfm_vj = zfvj - ABS( zfvj )
158               ! centered scheme
159               zcenut = zfui * ( tn(ji,jj,jk) + tn(ji+1,jj  ,jk) )
160               zcenvt = zfvj * ( tn(ji,jj,jk) + tn(ji  ,jj+1,jk) )
161               zcenus = zfui * ( sn(ji,jj,jk) + sn(ji+1,jj  ,jk) )
162               zcenvs = zfvj * ( sn(ji,jj,jk) + sn(ji  ,jj+1,jk) )
163               ! mixed centered / upstream scheme
164               zwx(ji,jj,jk) = zcenut - zfp_ui * zltu(ji,jj,jk) -zfm_ui * zltu(ji+1,jj,jk)
165               zwy(ji,jj,jk) = zcenvt - zfp_vj * zltv(ji,jj,jk) -zfm_vj * zltv(ji,jj+1,jk)
166               zww(ji,jj,jk) = zcenus - zfp_ui * zlsu(ji,jj,jk) -zfm_ui * zlsu(ji+1,jj,jk)
167               zwz(ji,jj,jk) = zcenvs - zfp_vj * zlsv(ji,jj,jk) -zfm_vj * zlsv(ji,jj+1,jk)
168            END DO
169         END DO
170
171         !  Tracer flux divergence at t-point added to the general trend
172         DO jj = 2, jpjm1
173            DO ji = fs_2, fs_jpim1   ! vector opt.
174               ! horizontal advective trends
175#if defined key_zco
176               zbtr = e1e2tr(ji,jj)
177#else
178               zbtr = e1e2tr(ji,jj) / fse3t(ji,jj,jk)
179#endif
180               zta = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk)   &
181                  &            + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk)  )
182               zsa = - zbtr * (  zww(ji,jj,jk) - zww(ji-1,jj  ,jk)   &
183                  &            + zwz(ji,jj,jk) - zwz(ji  ,jj-1,jk)  )
184               ! add it to the general tracer trends
185               ta(ji,jj,jk) = ta(ji,jj,jk) + zta
186               sa(ji,jj,jk) = sa(ji,jj,jk) + zsa
187            END DO
188         END DO
189         !                                             ! ===============
190      END DO                                           !   End of slab
191      !                                                ! ===============
192
193      ! Horizontal trend used in tra_adv_ztvd subroutine
194      zltu(:,:,:) = ta(:,:,:) - ztrdt(:,:,:) 
195      zlsu(:,:,:) = sa(:,:,:) - ztrds(:,:,:) 
196
197      ! 3. Save the horizontal advective trends for diagnostic
198      ! ------------------------------------------------------
199      IF( l_trdtra )   THEN
200         ! Recompute the hoizontal advection zta & zsa trends computed
201         ! at the step 2. above in making the difference between the new
202         ! trends and the previous one ta()/sa - ztrdt()/ztrds() and add
203         ! the term tn()/sn()*hdivn() to recover the Uh gradh(T/S) trends
204         ztrdt(:,:,:) = 0.e0   ;   ztrds(:,:,:) = 0.e0
205         !
206         ! T/S ZONAL advection trends
207         DO jk = 1, jpkm1
208            DO jj = 2, jpjm1
209               DO ji = fs_2, fs_jpim1   ! vector opt.
210                  !-- Compute zonal divergence by splitting hdivn (see divcur.F90)
211#if defined key_zco
212                  zbtr = e1e2tr(ji,jj)
213                  z_hdivn_x = (  e2u(ji  ,jj) * pun(ji  ,jj,jk)          &
214                     &         - e2u(ji-1,jj) * pun(ji-1,jj,jk) ) * zbtr
215#else
216                  zbtr = e1e2tr(ji,jj) / fse3t(ji,jj,jk)
217                  z_hdivn_x = (  e2u(ji  ,jj) * fse3u(ji  ,jj,jk) * pun(ji  ,jj,jk)          &
218                     &         - e2u(ji-1,jj) * fse3u(ji-1,jj,jk) * pun(ji-1,jj,jk) ) * zbtr
219#endif
220                  ztrdt(ji,jj,jk) = - ( zwx(ji,jj,jk) - zwx(ji-1,jj,jk) ) * zbtr + tn(ji,jj,jk) * z_hdivn_x
221                  ztrds(ji,jj,jk) = - ( zww(ji,jj,jk) - zww(ji-1,jj,jk) ) * zbtr + sn(ji,jj,jk) * z_hdivn_x
222               END DO
223            END DO
224         END DO
225         CALL trd_mod(ztrdt, ztrds, jptra_trd_xad, 'TRA', kt)    ! save the trends
226         !
227         ! T/S MERIDIONAL advection trends
228         DO jk = 1, jpkm1
229            DO jj = 2, jpjm1
230               DO ji = fs_2, fs_jpim1   ! vector opt.
231                  !-- Compute merid. divergence by splitting hdivn (see divcur.F90)
232#if defined key_zco
233                  zbtr      = e1e2tr(ji,jj)
234                  z_hdivn_y = (  e1v(ji,  jj) * pvn(ji,jj  ,jk)          &
235                     &         - e1v(ji,jj-1) * pvn(ji,jj-1,jk) ) * zbtr
236#else
237                  zbtr      = e1e2tr(ji,jj) / fse3t(ji,jj,jk)
238                  z_hdivn_y = (  e1v(ji,  jj) * fse3v(ji,jj  ,jk) * pvn(ji,jj  ,jk)          &
239                     &         - e1v(ji,jj-1) * fse3v(ji,jj-1,jk) * pvn(ji,jj-1,jk) ) * zbtr
240#endif
241                  ztrdt(ji,jj,jk) = - ( zwy(ji,jj,jk) - zwy(ji,jj-1,jk) ) * zbtr + tn(ji,jj,jk) * z_hdivn_y         
242                  ztrds(ji,jj,jk) = - ( zwz(ji,jj,jk) - zwz(ji,jj-1,jk) ) * zbtr + sn(ji,jj,jk) * z_hdivn_y
243               END DO
244            END DO
245         END DO
246         CALL trd_mod(ztrdt, ztrds, jptra_trd_yad, 'TRA', kt)     ! save the trends
247         !
248      ENDIF
249
250      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ta, clinfo1=' ubs had  - Ta: ', mask1=tmask,   &
251         &                       tab3d_2=sa, clinfo2=           ' Sa: ', mask2=tmask, clinfo3='tra' )
252
253      ! "zonal" mean advective heat and salt transport
254      IF( ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN
255         IF( lk_zco ) THEN
256            DO jk = 1, jpkm1
257               DO jj = 2, jpjm1
258                  DO ji = fs_2, fs_jpim1   ! vector opt.
259                     zwy(ji,jj,jk) = zwy(ji,jj,jk) * fse3v(ji,jj,jk)
260                     zwz(ji,jj,jk) = zwz(ji,jj,jk) * fse3v(ji,jj,jk)
261                  END DO
262               END DO
263            END DO
264         ENDIF
265         pht_adv(:) = ptr_vj( zwy(:,:,:) )
266         pst_adv(:) = ptr_vj( zwz(:,:,:) )
267      ENDIF
268
269      ! II. Vertical advection
270      ! ----------------------
271      IF( l_trdtra ) THEN          ! Save ta and sa trends
272         ztrdt(:,:,:) = ta(:,:,:)
273         ztrds(:,:,:) = sa(:,:,:)
274      ENDIF
275   
276      ! TVD scheme the vertical direction 
277      CALL tra_adv_ztvd(kt, pwn, zltu, zlsu)
278
279      IF( l_trdtra )   THEN         !  Save the final vertical advective trends
280         DO jk = 1, jpkm1
281            DO jj = 2, jpjm1
282               DO ji = fs_2, fs_jpim1   ! vector opt.
283#if defined key_zco
284                  zbtr      = e1e2tr(ji,jj)
285                  z_hdivn_x = e2u(ji,jj)*pun(ji,jj,jk) - e2u(ji-1,jj)*pun(ji-1,jj,jk)
286                  z_hdivn_y = e1v(ji,jj)*pvn(ji,jj,jk) - e1v(ji,jj-1)*pvn(ji,jj-1,jk)
287#else
288                  zbtr      = e1e2tr(ji,jj) / fse3t(ji,jj,jk)
289                  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)
290                  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)
291#endif
292                  z_hdivn   = (z_hdivn_x + z_hdivn_y) * zbtr
293                  zbtr      = e1e2tr(ji,jj) / fse3t(ji,jj,jk)
294                  ztrdt(ji,jj,jk) = ta(ji,jj,jk) - ztrdt(ji,jj,jk) - tn(ji,jj,jk) * z_hdivn
295                  ztrds(ji,jj,jk) = sa(ji,jj,jk) - ztrds(ji,jj,jk) - sn(ji,jj,jk) * z_hdivn
296               END DO
297            END DO
298         END DO
299         CALL trd_mod(ztrdt, ztrds, jptra_trd_zad, 'TRA', kt)   ! <<< ADD TO PREVIOUSLY COMPUTED
300         !
301      ENDIF
302
303      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ta, clinfo1=' ubs zad  - Ta: ', mask1=tmask,   &
304         &                       tab3d_2=sa, clinfo2=           ' Sa: ', mask2=tmask, clinfo3='tra')
305      !
306   END SUBROUTINE tra_adv_ubs
307
308
309   SUBROUTINE tra_adv_ztvd( kt, pwn, zttrd, zstrd )
310      !!----------------------------------------------------------------------
311      !!                  ***  ROUTINE tra_adv_ztvd  ***
312      !!
313      !! **  Purpose :   Compute the now trend due to total advection of
314      !!       tracers and add it to the general trend of tracer equations
315      !!
316      !! **  Method  :   TVD scheme, i.e. 2nd order centered scheme with
317      !!       corrected flux (monotonic correction)
318      !!       note: - this advection scheme needs a leap-frog time scheme
319      !!
320      !! ** Action : - update (ta,sa) with the now advective tracer trends
321      !!             - save the trends in (ztrdt,ztrds) ('key_trdtra')
322      !!----------------------------------------------------------------------
323      INTEGER , INTENT(in)                         ::   kt             ! ocean time-step
324      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::   pwn            ! verical effective velocity
325      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::   zttrd, zstrd   ! lateral advective trends on T & S
326      !!
327      INTEGER  ::   ji, jj, jk              ! dummy loop indices
328      REAL(wp) ::   z2dtt, zbtr, zew, z2    ! temporary scalar 
329      REAL(wp) ::   ztak, zfp_wk            !    "         "
330      REAL(wp) ::   zsak, zfm_wk            !    "         "
331      REAL(wp), DIMENSION (jpi,jpj,jpk) ::   zti, ztw   ! temporary 3D workspace
332      REAL(wp), DIMENSION (jpi,jpj,jpk) ::   zsi, zsw   !    "              "
333      !!----------------------------------------------------------------------
334
335      IF( kt == nit000 .AND. lwp ) THEN
336         WRITE(numout,*)
337         WRITE(numout,*) 'tra_adv_ztvd : vertical TVD advection scheme'
338         WRITE(numout,*) '~~~~~~~~~~~~'
339      ENDIF
340
341      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;    z2 = 1.
342      ELSE                                        ;    z2 = 2.
343      ENDIF
344
345      !  Bottom value : flux set to zero
346      ! --------------
347      ztw(:,:,jpk) = 0.e0   ;   zsw(:,:,jpk) = 0.e0
348      zti  (:,:,:) = 0.e0   ;   zsi  (:,:,:) = 0.e0
349
350
351      !  upstream advection with initial mass fluxes & intermediate update
352      ! -------------------------------------------------------------------
353      ! Surface value
354      IF( lk_dynspg_rl ) THEN                           ! rigid lid : flux set to zero
355         ztw(:,:,1) = 0.e0
356         zsw(:,:,1) = 0.e0
357      ELSE                                              ! free surface
358         DO jj = 1, jpj
359            DO ji = 1, jpi
360               zew = e1t(ji,jj) * e2t(ji,jj) * pwn(ji,jj,1)
361               ztw(ji,jj,1) = zew * tb(ji,jj,1)
362               zsw(ji,jj,1) = zew * sb(ji,jj,1)
363            END DO
364         END DO
365      ENDIF
366
367      ! Interior value
368      DO jk = 2, jpkm1
369         DO jj = 1, jpj
370            DO ji = 1, jpi
371               zew = 0.5 * e1t(ji,jj) * e2t(ji,jj) * pwn(ji,jj,jk)
372               zfp_wk = zew + ABS( zew )
373               zfm_wk = zew - ABS( zew )
374               ztw(ji,jj,jk) = zfp_wk * tb(ji,jj,jk) + zfm_wk * tb(ji,jj,jk-1)
375               zsw(ji,jj,jk) = zfp_wk * sb(ji,jj,jk) + zfm_wk * sb(ji,jj,jk-1)
376            END DO
377         END DO
378      END DO
379
380      ! update and guess with monotonic sheme
381      DO jk = 1, jpkm1
382         z2dtt = z2 * rdttra(jk)
383         DO jj = 2, jpjm1
384            DO ji = fs_2, fs_jpim1   ! vector opt.
385               zbtr = 1./ ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
386               ztak = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * zbtr
387               zsak = - ( zsw(ji,jj,jk) - zsw(ji,jj,jk+1) ) * zbtr
388               ta(ji,jj,jk) =  ta(ji,jj,jk) + ztak
389               sa(ji,jj,jk) =  sa(ji,jj,jk) + zsak 
390               zti (ji,jj,jk) = ( tb(ji,jj,jk) + z2dtt * ( ztak + zttrd(ji,jj,jk) ) ) * tmask(ji,jj,jk)
391               zsi (ji,jj,jk) = ( sb(ji,jj,jk) + z2dtt * ( zsak + zstrd(ji,jj,jk) ) ) * tmask(ji,jj,jk)
392            END DO
393         END DO
394      END DO
395
396      ! Lateral boundary conditions on zti, zsi   (unchanged sign)
397      CALL lbc_lnk( zti, 'T', 1. )
398      CALL lbc_lnk( zsi, 'T', 1. )
399
400
401      !  antidiffusive flux : high order minus low order
402      ! -------------------------------------------------     
403      ! Surface value
404      ztw(:,:,1) = 0.e0   ;   zsw(:,:,1) = 0.e0
405
406      ! Interior value
407      DO jk = 2, jpkm1
408         DO jj = 1, jpj
409            DO ji = 1, jpi
410               zew = 0.5 * e1t(ji,jj) * e2t(ji,jj) * pwn(ji,jj,jk)
411               ztw(ji,jj,jk) = zew * ( tn(ji,jj,jk) + tn(ji,jj,jk-1) ) - ztw(ji,jj,jk)
412               zsw(ji,jj,jk) = zew * ( sn(ji,jj,jk) + sn(ji,jj,jk-1) ) - zsw(ji,jj,jk)
413            END DO
414         END DO
415      END DO
416
417      !  monotonicity algorithm
418      ! ------------------------
419      CALL nonosc_z( tb, ztw, zti, z2 )
420      CALL nonosc_z( sb, zsw, zsi, z2 )
421
422
423      !  final trend with corrected fluxes
424      ! -----------------------------------
425      DO jk = 1, jpkm1
426         DO jj = 2, jpjm1
427            DO ji = fs_2, fs_jpim1   ! vector opt. 
428               zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
429               ! k- vertical advective trends
430               ztak = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * zbtr
431               zsak = - ( zsw(ji,jj,jk) - zsw(ji,jj,jk+1) ) * zbtr
432               ! add them to the general tracer trends
433               ta(ji,jj,jk) = ta(ji,jj,jk) + ztak
434               sa(ji,jj,jk) = sa(ji,jj,jk) + zsak
435            END DO
436         END DO
437      END DO
438      !
439   END SUBROUTINE tra_adv_ztvd
440
441
442   SUBROUTINE nonosc_z( pbef, pcc, paft, prdt )
443      !!---------------------------------------------------------------------
444      !!                    ***  ROUTINE nonosc_z  ***
445      !!     
446      !! **  Purpose :   compute monotonic tracer fluxes from the upstream
447      !!       scheme and the before field by a nonoscillatory algorithm
448      !!
449      !! **  Method  :   ... ???
450      !!       warning : pbef and paft must be masked, but the boundaries
451      !!       conditions on the fluxes are not necessary zalezak (1979)
452      !!       drange (1995) multi-dimensional forward-in-time and upstream-
453      !!       in-space based differencing for fluid
454      !!----------------------------------------------------------------------
455      REAL(wp), INTENT(in   )                          ::   prdt   ! ???
456      REAL(wp), INTENT(inout), DIMENSION (jpi,jpj,jpk) ::   pbef   ! before field
457      REAL(wp), INTENT(inout), DIMENSION (jpi,jpj,jpk) ::   paft   ! after field
458      REAL(wp), INTENT(inout), DIMENSION (jpi,jpj,jpk) ::   pcc    ! monotonic flux in the k direction
459      !!
460      INTEGER  ::   ji, jj, jk               ! dummy loop indices
461      INTEGER  ::   ikm1
462      REAL(wp) ::   zpos, zneg, zbt, za, zb, zc, zbig, zrtrn, z2dtt
463      REAL(wp), DIMENSION (jpi,jpj,jpk) ::   zbetup, zbetdo
464      !!----------------------------------------------------------------------
465
466      zbig = 1.e+40
467      zrtrn = 1.e-15
468      zbetup(:,:,:) = 0.e0   ;   zbetdo(:,:,:) = 0.e0
469
470      ! Search local extrema
471      ! --------------------
472      ! large negative value (-zbig) inside land
473      pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:) - zbig * ( 1.e0 - tmask(:,:,:) )
474      paft(:,:,:) = paft(:,:,:) * tmask(:,:,:) - zbig * ( 1.e0 - tmask(:,:,:) )
475      ! search maximum in neighbourhood
476      DO jk = 1, jpkm1
477         ikm1 = MAX(jk-1,1)
478         DO jj = 2, jpjm1
479            DO ji = fs_2, fs_jpim1   ! vector opt.
480               zbetup(ji,jj,jk) = MAX(  pbef(ji  ,jj  ,jk  ), paft(ji  ,jj  ,jk  ),   &
481                  &                     pbef(ji  ,jj  ,ikm1), pbef(ji  ,jj  ,jk+1),   &
482                  &                     paft(ji  ,jj  ,ikm1), paft(ji  ,jj  ,jk+1)  )
483            END DO
484         END DO
485      END DO
486      ! large positive value (+zbig) inside land
487      pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:) + zbig * ( 1.e0 - tmask(:,:,:) )
488      paft(:,:,:) = paft(:,:,:) * tmask(:,:,:) + zbig * ( 1.e0 - tmask(:,:,:) )
489      ! search minimum in neighbourhood
490      DO jk = 1, jpkm1
491         ikm1 = MAX(jk-1,1)
492         DO jj = 2, jpjm1
493            DO ji = fs_2, fs_jpim1   ! vector opt.
494               zbetdo(ji,jj,jk) = MIN(  pbef(ji  ,jj  ,jk  ), paft(ji  ,jj  ,jk  ),   &
495                  &                     pbef(ji  ,jj  ,ikm1), pbef(ji  ,jj  ,jk+1),   &
496                  &                     paft(ji  ,jj  ,ikm1), paft(ji  ,jj  ,jk+1)  )
497            END DO
498         END DO
499      END DO
500
501      ! restore masked values to zero
502      pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:)
503      paft(:,:,:) = paft(:,:,:) * tmask(:,:,:)
504 
505
506      ! 2. Positive and negative part of fluxes and beta terms
507      ! ------------------------------------------------------
508
509      DO jk = 1, jpkm1
510         z2dtt = prdt * rdttra(jk)
511         DO jj = 2, jpjm1
512            DO ji = fs_2, fs_jpim1   ! vector opt.
513               ! positive & negative part of the flux
514               zpos = MAX( 0., pcc(ji  ,jj  ,jk+1) ) - MIN( 0., pcc(ji  ,jj  ,jk  ) )
515               zneg = MAX( 0., pcc(ji  ,jj  ,jk  ) ) - MIN( 0., pcc(ji  ,jj  ,jk+1) )
516               ! up & down beta terms
517               zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) / z2dtt
518               zbetup(ji,jj,jk) = ( zbetup(ji,jj,jk) - paft(ji,jj,jk) ) / (zpos+zrtrn) * zbt
519               zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zbetdo(ji,jj,jk) ) / (zneg+zrtrn) * zbt
520            END DO
521         END DO
522      END DO
523
524      ! monotonic flux in the k direction, i.e. pcc
525      ! -------------------------------------------
526      DO jk = 2, jpkm1
527         DO jj = 2, jpjm1
528            DO ji = fs_2, fs_jpim1   ! vector opt.
529
530               za = MIN( 1., zbetdo(ji,jj,jk), zbetup(ji,jj,jk-1) )
531               zb = MIN( 1., zbetup(ji,jj,jk), zbetdo(ji,jj,jk-1) )
532               zc = 0.5 * ( 1.e0 + SIGN( 1.e0, pcc(ji,jj,jk) ) )
533               pcc(ji,jj,jk) = pcc(ji,jj,jk) * ( zc * za + ( 1.e0 - zc) * zb )
534            END DO
535         END DO
536      END DO
537      !
538   END SUBROUTINE nonosc_z
539
540   !!======================================================================
541END MODULE traadv_ubs
Note: See TracBrowser for help on using the repository browser.