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

source: trunk/NEMO/OPA_SRC/TRA/traadv_tvd.F90 @ 719

Last change on this file since 719 was 719, checked in by ctlod, 17 years ago

get back to the nemo_v2_3 version for trunk
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 24.1 KB
Line 
1MODULE traadv_tvd
2   !!==============================================================================
3   !!                       ***  MODULE  traadv_tvd  ***
4   !! Ocean active tracers:  horizontal & vertical advective trend
5   !!==============================================================================
6   !! History :       !  95-12  (L. Mortier)  Original code
7   !!                 !  00-01  (H. Loukos)  adapted to ORCA
8   !!                 !  00-10  (MA Foujols E.Kestenare)  include file not routine
9   !!                 !  00-12  (E. Kestenare M. Levy)  fix bug in trtrd indexes
10   !!                 !  01-07  (E. Durand G. Madec)  adaptation to ORCA config
11   !!            8.5  !  02-06  (G. Madec)  F90: Free form and module
12   !!            9.0  !  04-01  (A. de Miranda, G. Madec, J.M. Molines ): advective bbl
13   !!            9.0  !  08-04  (S. Cravatte) add the i-, j- & k- trends computation
14   !!            " "  !  05-11  (V. Garnier) Surface pressure gradient organization
15   !!----------------------------------------------------------------------
16
17
18   !!----------------------------------------------------------------------
19   !!   tra_adv_tvd  : update the tracer trend with the horizontal
20   !!                  and vertical advection trends using a TVD scheme
21   !!   nonosc       : compute monotonic tracer fluxes by a nonoscillatory
22   !!                  algorithm
23   !!----------------------------------------------------------------------
24   USE oce             ! ocean dynamics and active tracers
25   USE dom_oce         ! ocean space and time domain
26   USE trdmod          ! ocean active tracers trends
27   USE trdmod_oce      ! ocean variables trends
28   USE in_out_manager  ! I/O manager
29   USE dynspg_oce      ! choice/control of key cpp for surface pressure gradient
30   USE trabbl          ! Advective term of BBL
31   USE lib_mpp
32   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
33   USE diaptr          ! poleward transport diagnostics
34   USE prtctl          ! Print control
35
36
37   IMPLICIT NONE
38   PRIVATE
39
40   PUBLIC   tra_adv_tvd    ! routine called by step.F90
41
42   !! * Substitutions
43#  include "domzgr_substitute.h90"
44#  include "vectopt_loop_substitute.h90"
45   !!----------------------------------------------------------------------
46   !!   OPA 9.0 , LOCEAN-IPSL (2006)
47   !! $Header$
48   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
49   !!----------------------------------------------------------------------
50
51CONTAINS
52
53   SUBROUTINE tra_adv_tvd( kt, pun, pvn, pwn )
54      !!----------------------------------------------------------------------
55      !!                  ***  ROUTINE tra_adv_tvd  ***
56      !!
57      !! **  Purpose :   Compute the now trend due to total advection of
58      !!       tracers and add it to the general trend of tracer equations
59      !!
60      !! **  Method  :   TVD scheme, i.e. 2nd order centered scheme with
61      !!       corrected flux (monotonic correction)
62      !!       note: - this advection scheme needs a leap-frog time scheme
63      !!
64      !! ** Action : - update (ta,sa) with the now advective tracer trends
65      !!             - save the trends in (ztrdt,ztrds) ('key_trdtra')
66      !!----------------------------------------------------------------------
67      USE oce              , ztrdt => ua   ! use ua as workspace
68      USE oce              , ztrds => va   ! use va as workspace
69      !!
70      INTEGER , INTENT(in)                         ::   kt    ! ocean time-step index
71      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::   pun   ! ocean velocity u-component
72      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::   pvn   ! ocean velocity v-component
73      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::   pwn   ! ocean velocity w-component
74      !!
75      INTEGER  ::   ji, jj, jk              ! dummy loop indices
76      REAL(wp) ::                        &  ! temporary scalar
77         ztai, ztaj, ztak,               &  !    "         "   
78         zsai, zsaj, zsak,               &  !    "         "   
79         z_hdivn_x, z_hdivn_y, z_hdivn
80      REAL(wp) ::   &
81         z2dtt, zbtr, zeu, zev,          &  ! temporary scalar
82         zew, z2, zbtr1,                 &  ! temporary scalar
83         zfp_ui, zfp_vj, zfp_wk,         &  !    "         "
84         zfm_ui, zfm_vj, zfm_wk             !    "         "
85      REAL(wp), DIMENSION (jpi,jpj,jpk) ::   zti, ztu, ztv, ztw   ! temporary workspace
86      REAL(wp), DIMENSION (jpi,jpj,jpk) ::   zsi, zsu, zsv, zsw   !    "           "
87      !!----------------------------------------------------------------------
88
89      zti(:,:,:) = 0.e0   ;   zsi(:,:,:) = 0.e0
90
91      IF( kt == nit000 .AND. lwp ) THEN
92         WRITE(numout,*)
93         WRITE(numout,*) 'tra_adv_tvd : TVD advection scheme'
94         WRITE(numout,*) '~~~~~~~~~~~'
95      ENDIF
96
97      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;    z2 = 1.
98      ELSE                                        ;    z2 = 2.
99      ENDIF
100
101      ! 1. Bottom value : flux set to zero
102      ! ---------------
103      ztu(:,:,jpk) = 0.e0   ;   zsu(:,:,jpk) = 0.e0
104      ztv(:,:,jpk) = 0.e0   ;   zsv(:,:,jpk) = 0.e0
105      ztw(:,:,jpk) = 0.e0   ;   zsw(:,:,jpk) = 0.e0
106      zti(:,:,jpk) = 0.e0   ;   zsi(:,:,jpk) = 0.e0
107
108
109      ! 2. upstream advection with initial mass fluxes & intermediate update
110      ! --------------------------------------------------------------------
111      ! upstream tracer flux in the i and j direction
112      DO jk = 1, jpkm1
113         DO jj = 1, jpjm1
114            DO ji = 1, fs_jpim1   ! vector opt.
115               zeu = 0.5 * e2u(ji,jj) * fse3u(ji,jj,jk) * pun(ji,jj,jk)
116               zev = 0.5 * e1v(ji,jj) * fse3v(ji,jj,jk) * pvn(ji,jj,jk)
117               ! upstream scheme
118               zfp_ui = zeu + ABS( zeu )
119               zfm_ui = zeu - ABS( zeu )
120               zfp_vj = zev + ABS( zev )
121               zfm_vj = zev - ABS( zev )
122               ztu(ji,jj,jk) = zfp_ui * tb(ji,jj,jk) + zfm_ui * tb(ji+1,jj  ,jk)
123               ztv(ji,jj,jk) = zfp_vj * tb(ji,jj,jk) + zfm_vj * tb(ji  ,jj+1,jk)
124               zsu(ji,jj,jk) = zfp_ui * sb(ji,jj,jk) + zfm_ui * sb(ji+1,jj  ,jk)
125               zsv(ji,jj,jk) = zfp_vj * sb(ji,jj,jk) + zfm_vj * sb(ji  ,jj+1,jk)
126            END DO
127         END DO
128      END DO
129
130      ! upstream tracer flux in the k direction
131      ! Surface value
132      IF( lk_dynspg_rl .OR. lk_vvl ) THEN               ! rigid lid or variable volume: flux set to zero
133         ztw(:,:,1) = 0.e0
134         zsw(:,:,1) = 0.e0
135      ELSE                                              ! free surface
136         DO jj = 1, jpj
137            DO ji = 1, jpi
138               zew = e1t(ji,jj) * e2t(ji,jj) * pwn(ji,jj,1)
139               ztw(ji,jj,1) = zew * tb(ji,jj,1)
140               zsw(ji,jj,1) = zew * sb(ji,jj,1)
141            END DO
142         END DO
143      ENDIF
144
145      ! Interior value
146      DO jk = 2, jpkm1
147         DO jj = 1, jpj
148            DO ji = 1, jpi
149               zew = 0.5 * e1t(ji,jj) * e2t(ji,jj) * pwn(ji,jj,jk)
150               zfp_wk = zew + ABS( zew )
151               zfm_wk = zew - ABS( zew )
152               ztw(ji,jj,jk) = zfp_wk * tb(ji,jj,jk) + zfm_wk * tb(ji,jj,jk-1)
153               zsw(ji,jj,jk) = zfp_wk * sb(ji,jj,jk) + zfm_wk * sb(ji,jj,jk-1)
154            END DO
155         END DO
156      END DO
157
158      ! total advective trend
159      DO jk = 1, jpkm1
160         DO jj = 2, jpjm1
161            DO ji = fs_2, fs_jpim1   ! vector opt.
162               zbtr = 1./ ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
163               ! i- j- horizontal & k- vertical advective trends
164               ztai = - ( ztu(ji,jj,jk) - ztu(ji-1,jj  ,jk  ) ) * zbtr
165               ztaj = - ( ztv(ji,jj,jk) - ztv(ji  ,jj-1,jk  ) ) * zbtr
166               ztak = - ( ztw(ji,jj,jk) - ztw(ji  ,jj  ,jk+1) ) * zbtr
167               zsai = - ( zsu(ji,jj,jk) - zsu(ji-1,jj  ,jk  ) ) * zbtr
168               zsaj = - ( zsv(ji,jj,jk) - zsv(ji  ,jj-1,jk  ) ) * zbtr
169               zsak = - ( zsw(ji,jj,jk) - zsw(ji  ,jj  ,jk+1) ) * zbtr
170               ! total intermediate advective trends
171               zti(ji,jj,jk) = ztai + ztaj + ztak
172               zsi(ji,jj,jk) = zsai + zsaj + zsak
173            END DO
174         END DO
175      END DO
176
177
178      ! Save the intermediate i / j / k advective trends for diagnostics
179      ! -------------------------------------------------------------------
180      ! Warning : We should use zun instead of un in the computations below, but we
181      ! also use hdivn which is computed with un, vn (check ???). So we use un, vn
182      ! for consistency. Results are therefore approximate with key_trabbl_adv.
183
184      IF( l_trdtra ) THEN
185         ztrdt(:,:,:) = 0.e0   ;   ztrds(:,:,:) = 0.e0
186         !
187         ! T/S ZONAL advection trends
188         DO jk = 1, jpkm1
189            DO jj = 2, jpjm1
190               DO ji = fs_2, fs_jpim1   ! vector opt.
191                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
192                  ztrdt(ji,jj,jk) = - ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk) ) * zbtr
193                  ztrds(ji,jj,jk) = - ( zsu(ji,jj,jk) - zsu(ji-1,jj,jk) ) * zbtr
194               END DO
195            END DO
196         END DO
197         CALL trd_mod(ztrdt, ztrds, jptra_trd_xad, 'TRA', kt)    ! save the trends
198         !
199         ! T/S MERIDIONAL advection trends
200         DO jk = 1, jpkm1
201            DO jj = 2, jpjm1
202               DO ji = fs_2, fs_jpim1   ! vector opt.
203                  zbtr      = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
204                  ztrdt(ji,jj,jk) = - ( ztv(ji,jj,jk) - ztv(ji,jj-1,jk) ) * zbtr
205                  ztrds(ji,jj,jk) = - ( zsv(ji,jj,jk) - zsv(ji,jj-1,jk) ) * zbtr
206               END DO
207            END DO
208         END DO
209         CALL trd_mod(ztrdt, ztrds, jptra_trd_yad, 'TRA', kt)     ! save the trends
210         !
211         ! T/S VERTICAL advection trends
212         DO jk = 1, jpkm1
213            DO jj = 2, jpjm1
214               DO ji = fs_2, fs_jpim1   ! vector opt.         
215                  zbtr      = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
216                  ztrdt(ji,jj,jk) = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * zbtr
217                  ztrds(ji,jj,jk) = - ( zsw(ji,jj,jk) - zsw(ji,jj,jk+1) ) * zbtr
218               END DO
219            END DO
220         END DO
221         CALL trd_mod(ztrdt, ztrds, jptra_trd_zad, 'TRA', kt)     ! save the trends
222         !
223      ENDIF
224
225      ! update and guess with monotonic sheme
226      DO jk = 1, jpkm1
227         z2dtt = z2 * rdttra(jk)
228         DO jj = 2, jpjm1
229            DO ji = fs_2, fs_jpim1   ! vector opt.
230               ta(ji,jj,jk) =  ta(ji,jj,jk) + zti(ji,jj,jk)
231               sa(ji,jj,jk) =  sa(ji,jj,jk) + zsi(ji,jj,jk)
232               zti (ji,jj,jk) = ( tb(ji,jj,jk) + z2dtt * zti(ji,jj,jk) ) * tmask(ji,jj,jk)
233               zsi (ji,jj,jk) = ( sb(ji,jj,jk) + z2dtt * zsi(ji,jj,jk) ) * tmask(ji,jj,jk)
234            END DO
235         END DO
236      END DO
237
238      ! Lateral boundary conditions on zti, zsi   (unchanged sign)
239      CALL lbc_lnk( zti, 'T', 1. )
240      CALL lbc_lnk( zsi, 'T', 1. )
241
242
243      ! 3. antidiffusive flux : high order minus low order
244      ! --------------------------------------------------
245      ! antidiffusive flux on i and j
246      DO jk = 1, jpkm1
247         DO jj = 1, jpjm1
248            DO ji = 1, fs_jpim1   ! vector opt.
249               zeu = 0.5 * e2u(ji,jj) * fse3u(ji,jj,jk) * pun(ji,jj,jk)
250               zev = 0.5 * e1v(ji,jj) * fse3v(ji,jj,jk) * pvn(ji,jj,jk)
251               ztu(ji,jj,jk) = zeu * ( tn(ji,jj,jk) + tn(ji+1,jj,jk) ) - ztu(ji,jj,jk)
252               zsu(ji,jj,jk) = zeu * ( sn(ji,jj,jk) + sn(ji+1,jj,jk) ) - zsu(ji,jj,jk)
253               ztv(ji,jj,jk) = zev * ( tn(ji,jj,jk) + tn(ji,jj+1,jk) ) - ztv(ji,jj,jk)
254               zsv(ji,jj,jk) = zev * ( sn(ji,jj,jk) + sn(ji,jj+1,jk) ) - zsv(ji,jj,jk)
255            END DO
256         END DO
257      END DO
258     
259      ! antidiffusive flux on k
260      ! Surface value
261      ztw(:,:,1) = 0.e0
262      zsw(:,:,1) = 0.e0
263
264      ! Interior value
265      DO jk = 2, jpkm1
266         DO jj = 1, jpj
267            DO ji = 1, jpi
268               zew = 0.5 * e1t(ji,jj) * e2t(ji,jj) * pwn(ji,jj,jk)
269               ztw(ji,jj,jk) = zew * ( tn(ji,jj,jk) + tn(ji,jj,jk-1) ) - ztw(ji,jj,jk)
270               zsw(ji,jj,jk) = zew * ( sn(ji,jj,jk) + sn(ji,jj,jk-1) ) - zsw(ji,jj,jk)
271            END DO
272         END DO
273      END DO
274
275      ! Lateral bondary conditions
276      CALL lbc_lnk( ztu, 'U', -1. )   ;   CALL lbc_lnk( zsu, 'U', -1. )
277      CALL lbc_lnk( ztv, 'V', -1. )   ;   CALL lbc_lnk( zsv, 'V', -1. )
278      CALL lbc_lnk( ztw, 'W',  1. )   ;   CALL lbc_lnk( zsw, 'W',  1. )
279
280      ! 4. monotonicity algorithm
281      ! -------------------------
282      CALL nonosc( tb, ztu, ztv, ztw, zti, z2 )
283      CALL nonosc( sb, zsu, zsv, zsw, zsi, z2 )
284
285
286      ! 5. final trend with corrected fluxes
287      ! ------------------------------------
288      DO jk = 1, jpkm1
289         DO jj = 2, jpjm1
290            DO ji = fs_2, fs_jpim1   ! vector opt. 
291               zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
292               ! i- j- horizontal & k- vertical advective trends
293               ztai = - ( ztu(ji,jj,jk) - ztu(ji-1,jj  ,jk  )) * zbtr
294               ztaj = - ( ztv(ji,jj,jk) - ztv(ji  ,jj-1,jk  )) * zbtr
295               ztak = - ( ztw(ji,jj,jk) - ztw(ji  ,jj  ,jk+1)) * zbtr
296               zsai = - ( zsu(ji,jj,jk) - zsu(ji-1,jj  ,jk  )) * zbtr
297               zsaj = - ( zsv(ji,jj,jk) - zsv(ji  ,jj-1,jk  )) * zbtr
298               zsak = - ( zsw(ji,jj,jk) - zsw(ji  ,jj  ,jk+1)) * zbtr
299
300               ! add them to the general tracer trends
301               ta(ji,jj,jk) = ta(ji,jj,jk) + ztai + ztaj + ztak
302               sa(ji,jj,jk) = sa(ji,jj,jk) + zsai + zsaj + zsak
303            END DO
304         END DO
305      END DO
306
307
308      ! Save the advective trends for diagnostics
309      ! --------------------------------------------
310
311      IF( l_trdtra ) THEN
312         ztrdt(:,:,:) = 0.e0   ;   ztrds(:,:,:) = 0.e0
313         !
314         ! T/S ZONAL advection trends
315         DO jk = 1, jpkm1
316            DO jj = 2, jpjm1
317               DO ji = fs_2, fs_jpim1   ! vector opt.
318                  !-- Compute zonal divergence by splitting hdivn (see divcur.F90)
319                  !   N.B. This computation is not valid along OBCs (if any)
320                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
321                  z_hdivn_x = (  e2u(ji  ,jj) * fse3u(ji  ,jj,jk) * pun(ji  ,jj,jk)          &
322                     &         - e2u(ji-1,jj) * fse3u(ji-1,jj,jk) * pun(ji-1,jj,jk) ) * zbtr
323                  !-- Compute T/S zonal advection trends
324                  ztrdt(ji,jj,jk) = - ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk) ) * zbtr + tn(ji,jj,jk) * z_hdivn_x
325                  ztrds(ji,jj,jk) = - ( zsu(ji,jj,jk) - zsu(ji-1,jj,jk) ) * zbtr + sn(ji,jj,jk) * z_hdivn_x
326               END DO
327            END DO
328         END DO
329         CALL trd_mod(ztrdt, ztrds, jptra_trd_xad, 'TRA', kt, cnbpas='bis')   ! <<< ADD TO PREVIOUSLY COMPUTED
330         !
331         ! T/S MERIDIONAL advection trends
332         DO jk = 1, jpkm1
333            DO jj = 2, jpjm1
334               DO ji = fs_2, fs_jpim1   ! vector opt.
335                  !-- Compute merid. divergence by splitting hdivn (see divcur.F90)
336                  !   N.B. This computation is not valid along OBCs (if any)
337                  zbtr      = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
338                  z_hdivn_y = (  e1v(ji,  jj) * fse3v(ji,jj  ,jk) * pvn(ji,jj  ,jk)          &
339                     &         - e1v(ji,jj-1) * fse3v(ji,jj-1,jk) * pvn(ji,jj-1,jk) ) * zbtr
340                  !-- Compute T/S meridional advection trends
341                  ztrdt(ji,jj,jk) = - ( ztv(ji,jj,jk) - ztv(ji,jj-1,jk) ) * zbtr + tn(ji,jj,jk) * z_hdivn_y         
342                  ztrds(ji,jj,jk) = - ( zsv(ji,jj,jk) - zsv(ji,jj-1,jk) ) * zbtr + sn(ji,jj,jk) * z_hdivn_y         
343               END DO
344            END DO
345         END DO
346         CALL trd_mod(ztrdt, ztrds, jptra_trd_yad, 'TRA', kt, cnbpas='bis')   ! <<< ADD TO PREVIOUSLY COMPUTED
347         !
348         ! T/S VERTICAL advection trends
349         DO jk = 1, jpkm1
350            DO jj = 2, jpjm1
351               DO ji = fs_2, fs_jpim1   ! vector opt.
352                  zbtr1     = 1. / ( e1t(ji,jj) * e2t(ji,jj) )
353#if defined key_zco
354                  zbtr      = zbtr1
355                  z_hdivn_x = e2u(ji,jj)*pun(ji,jj,jk) - e2u(ji-1,jj)*pun(ji-1,jj,jk)
356                  z_hdivn_y = e1v(ji,jj)*pvn(ji,jj,jk) - e1v(ji,jj-1)*pvn(ji,jj-1,jk)
357#else
358                  zbtr      = zbtr1 / fse3t(ji,jj,jk)
359                  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)
360                  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)
361#endif
362                  z_hdivn   = (z_hdivn_x + z_hdivn_y) * zbtr
363                  zbtr      = zbtr1 / fse3t(ji,jj,jk)
364                  ztrdt(ji,jj,jk) = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * zbtr - tn(ji,jj,jk) * z_hdivn
365                  ztrds(ji,jj,jk) = - ( zsw(ji,jj,jk) - zsw(ji,jj,jk+1) ) * zbtr - sn(ji,jj,jk) * z_hdivn
366               END DO
367            END DO
368         END DO
369         CALL trd_mod(ztrdt, ztrds, jptra_trd_zad, 'TRA', kt, cnbpas='bis')   ! <<< ADD TO PREVIOUSLY COMPUTED
370         !
371      ENDIF
372
373      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ta, clinfo1=' tvd adv  - Ta: ', mask1=tmask,   &
374         &                       tab3d_2=sa, clinfo2=           ' Sa: ', mask2=tmask, clinfo3='tra' )
375
376      ! "zonal" mean advective heat and salt transport
377      IF( ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN
378         pht_adv(:) = ptr_vj( ztv(:,:,:) )
379         pst_adv(:) = ptr_vj( zsv(:,:,:) )
380      ENDIF
381      !
382   END SUBROUTINE tra_adv_tvd
383
384
385   SUBROUTINE nonosc( pbef, paa, pbb, pcc, paft, prdt )
386      !!---------------------------------------------------------------------
387      !!                    ***  ROUTINE nonosc  ***
388      !!     
389      !! **  Purpose :   compute monotonic tracer fluxes from the upstream
390      !!       scheme and the before field by a nonoscillatory algorithm
391      !!
392      !! **  Method  :   ... ???
393      !!       warning : pbef and paft must be masked, but the boundaries
394      !!       conditions on the fluxes are not necessary zalezak (1979)
395      !!       drange (1995) multi-dimensional forward-in-time and upstream-
396      !!       in-space based differencing for fluid
397      !!----------------------------------------------------------------------
398      REAL(wp), INTENT( in ) ::   prdt                               ! ???
399      REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT( inout ) ::   &
400         pbef,                            & ! before field
401         paft,                            & ! after field
402         paa,                             & ! monotonic flux in the i direction
403         pbb,                             & ! monotonic flux in the j direction
404         pcc                                ! monotonic flux in the k direction
405      !!
406      INTEGER ::   ji, jj, jk               ! dummy loop indices
407      INTEGER ::   ikm1
408      REAL(wp), DIMENSION (jpi,jpj,jpk) ::   zbetup, zbetdo
409      REAL(wp) ::   zpos, zneg, zbt, za, zb, zc, zbig, zrtrn, z2dtt
410      !!----------------------------------------------------------------------
411
412      zbig = 1.e+40
413      zrtrn = 1.e-15
414      zbetup(:,:,:) = 0.e0   ;   zbetdo(:,:,:) = 0.e0 
415
416      ! Search local extrema
417      ! --------------------
418      ! large negative value (-zbig) inside land
419      pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:) - zbig * ( 1.e0 - tmask(:,:,:) )
420      paft(:,:,:) = paft(:,:,:) * tmask(:,:,:) - zbig * ( 1.e0 - tmask(:,:,:) )
421      ! search maximum in neighbourhood
422      DO jk = 1, jpkm1
423         ikm1 = MAX(jk-1,1)
424         DO jj = 2, jpjm1
425            DO ji = fs_2, fs_jpim1   ! vector opt.
426               zbetup(ji,jj,jk) = MAX(  pbef(ji  ,jj  ,jk  ), paft(ji  ,jj  ,jk  ),   &
427                  &                     pbef(ji-1,jj  ,jk  ), pbef(ji+1,jj  ,jk  ),   &
428                  &                     paft(ji-1,jj  ,jk  ), paft(ji+1,jj  ,jk  ),   &
429                  &                     pbef(ji  ,jj-1,jk  ), pbef(ji  ,jj+1,jk  ),   &
430                  &                     paft(ji  ,jj-1,jk  ), paft(ji  ,jj+1,jk  ),   &
431                  &                     pbef(ji  ,jj  ,ikm1), pbef(ji  ,jj  ,jk+1),   &
432                  &                     paft(ji  ,jj  ,ikm1), paft(ji  ,jj  ,jk+1)  )
433            END DO
434         END DO
435      END DO
436      ! large positive value (+zbig) inside land
437      pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:) + zbig * ( 1.e0 - tmask(:,:,:) )
438      paft(:,:,:) = paft(:,:,:) * tmask(:,:,:) + zbig * ( 1.e0 - tmask(:,:,:) )
439      ! search minimum in neighbourhood
440      DO jk = 1, jpkm1
441         ikm1 = MAX(jk-1,1)
442         DO jj = 2, jpjm1
443            DO ji = fs_2, fs_jpim1   ! vector opt.
444               zbetdo(ji,jj,jk) = MIN(  pbef(ji  ,jj  ,jk  ), paft(ji  ,jj  ,jk  ),   &
445                  &                     pbef(ji-1,jj  ,jk  ), pbef(ji+1,jj  ,jk  ),   &
446                  &                     paft(ji-1,jj  ,jk  ), paft(ji+1,jj  ,jk  ),   &
447                  &                     pbef(ji  ,jj-1,jk  ), pbef(ji  ,jj+1,jk  ),   &
448                  &                     paft(ji  ,jj-1,jk  ), paft(ji  ,jj+1,jk  ),   &
449                  &                     pbef(ji  ,jj  ,ikm1), pbef(ji  ,jj  ,jk+1),   &
450                  &                     paft(ji  ,jj  ,ikm1), paft(ji  ,jj  ,jk+1)  )
451            END DO
452         END DO
453      END DO
454
455      ! restore masked values to zero
456      pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:)
457      paft(:,:,:) = paft(:,:,:) * tmask(:,:,:)
458 
459
460      ! 2. Positive and negative part of fluxes and beta terms
461      ! ------------------------------------------------------
462
463      DO jk = 1, jpkm1
464         z2dtt = prdt * rdttra(jk)
465         DO jj = 2, jpjm1
466            DO ji = fs_2, fs_jpim1   ! vector opt.
467               ! positive & negative part of the flux
468               zpos = MAX( 0., paa(ji-1,jj  ,jk  ) ) - MIN( 0., paa(ji  ,jj  ,jk  ) )   &
469                  & + MAX( 0., pbb(ji  ,jj-1,jk  ) ) - MIN( 0., pbb(ji  ,jj  ,jk  ) )   &
470                  & + MAX( 0., pcc(ji  ,jj  ,jk+1) ) - MIN( 0., pcc(ji  ,jj  ,jk  ) )
471               zneg = MAX( 0., paa(ji  ,jj  ,jk  ) ) - MIN( 0., paa(ji-1,jj  ,jk  ) )   &
472                  & + MAX( 0., pbb(ji  ,jj  ,jk  ) ) - MIN( 0., pbb(ji  ,jj-1,jk  ) )   &
473                  & + MAX( 0., pcc(ji  ,jj  ,jk  ) ) - MIN( 0., pcc(ji  ,jj  ,jk+1) )
474               ! up & down beta terms
475               zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) / z2dtt
476               zbetup(ji,jj,jk) = ( zbetup(ji,jj,jk) - paft(ji,jj,jk) ) / (zpos+zrtrn) * zbt
477               zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zbetdo(ji,jj,jk) ) / (zneg+zrtrn) * zbt
478            END DO
479         END DO
480      END DO
481
482      ! lateral boundary condition on zbetup & zbetdo   (unchanged sign)
483      CALL lbc_lnk( zbetup, 'T', 1. )
484      CALL lbc_lnk( zbetdo, 'T', 1. )
485
486
487      ! 3. monotonic flux in the i & j direction (paa & pbb)
488      ! ----------------------------------------
489      DO jk = 1, jpkm1
490         DO jj = 2, jpjm1
491            DO ji = fs_2, fs_jpim1   ! vector opt.
492               za = MIN( 1.e0, zbetdo(ji,jj,jk), zbetup(ji+1,jj,jk) )
493               zb = MIN( 1.e0, zbetup(ji,jj,jk), zbetdo(ji+1,jj,jk) )
494               zc = 0.5 * ( 1.e0 + SIGN( 1.e0, paa(ji,jj,jk) ) )
495               paa(ji,jj,jk) = paa(ji,jj,jk) * ( zc * za + ( 1.e0 - zc) * zb )
496
497               za = MIN( 1.e0, zbetdo(ji,jj,jk), zbetup(ji,jj+1,jk) )
498               zb = MIN( 1.e0, zbetup(ji,jj,jk), zbetdo(ji,jj+1,jk) )
499               zc = 0.5 * ( 1.e0 + SIGN( 1.e0, pbb(ji,jj,jk) ) )
500               pbb(ji,jj,jk) = pbb(ji,jj,jk) * ( zc * za + ( 1.e0 - zc) * zb )
501            END DO
502         END DO
503      END DO
504
505
506      ! monotonic flux in the k direction, i.e. pcc
507      ! -------------------------------------------
508      DO jk = 2, jpkm1
509         DO jj = 2, jpjm1
510            DO ji = fs_2, fs_jpim1   ! vector opt.
511
512               za = MIN( 1., zbetdo(ji,jj,jk), zbetup(ji,jj,jk-1) )
513               zb = MIN( 1., zbetup(ji,jj,jk), zbetdo(ji,jj,jk-1) )
514               zc = 0.5 * ( 1.e0 + SIGN( 1.e0, pcc(ji,jj,jk) ) )
515               pcc(ji,jj,jk) = pcc(ji,jj,jk) * ( zc * za + ( 1.e0 - zc) * zb )
516            END DO
517         END DO
518      END DO
519
520      ! lateral boundary condition on paa, pbb, pcc
521      CALL lbc_lnk( paa, 'U', -1. )      ! changed sign
522      CALL lbc_lnk( pbb, 'V', -1. )      ! changed sign
523      CALL lbc_lnk( pcc, 'W',  1. )      ! NO changed sign
524      !
525   END SUBROUTINE nonosc
526
527   !!======================================================================
528END MODULE traadv_tvd
Note: See TracBrowser for help on using the repository browser.