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 @ 786

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

dev_001_GM - merge TRC-TRA on OPA only, trabbl & zpshde not done and trdmld not OK - compilation OK
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 17.5 KB
Line 
1MODULE traadv_tvd
2   !!==============================================================================
3   !!                       ***  MODULE  traadv_tvd  ***
4   !! Ocean active tracers:  horizontal & vertical advective trend
5   !!==============================================================================
6   !! History :  7.0  !  95-12  (L. Mortier)  Original code
7   !!            8.0  !  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   !!   NEMO     1.0  !  04-01  (A. de Miranda, G. Madec, J.M. Molines ): advective bbl
13   !!             -   !  08-04  (S. Cravatte) add the i-, j- & k- trends computation
14   !!             -   !  05-11  (V. Garnier) Surface pressure gradient organization
15   !!            2.4  !  08-01  (G. Madec) Merge TRA-TRC
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) ::   ztai, ztaj, ztak
74      REAL(wp) ::   z2dtt, zbtr, zeu, zev   ! temporary scalar
75      REAL(wp) ::   zew, z2                          ! temporary scalar
76      REAL(wp) ::   zfp_ui, zfp_vj, zfp_wk           !    "         "
77      REAL(wp) ::   zfm_ui, zfm_vj, zfm_wk            !    "         "
78      REAL(wp), DIMENSION (jpi,jpj,jpk) ::   zti, ztu, ztv, ztw   ! temporary workspace
79      !!----------------------------------------------------------------------
80
81      zti(:,:,:) = 0.e0 
82
83      IF( kt == nit000 .AND. lwp ) THEN
84         WRITE(numout,*)
85         WRITE(numout,*) 'tra_adv_tvd : TVD advection scheme'
86         WRITE(numout,*) '~~~~~~~~~~~'
87      ENDIF
88
89      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;    z2 = 1.
90      ELSE                                        ;    z2 = 2.
91      ENDIF
92
93      ! 1. Bottom value : flux set to zero
94      ! ---------------
95      ztu(:,:,jpk) = 0.e0   ;   ztv(:,:,jpk) = 0.e0
96      ztw(:,:,jpk) = 0.e0   ;   zti(:,:,jpk) = 0.e0
97
98
99      ! 2. upstream advection with initial mass fluxes & intermediate update
100      ! --------------------------------------------------------------------
101      ! upstream tracer flux in the i and j direction
102      DO jk = 1, jpkm1
103         DO jj = 1, jpjm1
104            DO ji = 1, fs_jpim1   ! vector opt.
105               zeu = 0.5 * e2u(ji,jj) * fse3u(ji,jj,jk) * pun(ji,jj,jk)
106               zev = 0.5 * e1v(ji,jj) * fse3v(ji,jj,jk) * pvn(ji,jj,jk)
107               ! upstream scheme
108               zfp_ui = zeu + ABS( zeu )
109               zfm_ui = zeu - ABS( zeu )
110               zfp_vj = zev + ABS( zev )
111               zfm_vj = zev - ABS( zev )
112               ztu(ji,jj,jk) = zfp_ui * ptb(ji,jj,jk) + zfm_ui * ptb(ji+1,jj  ,jk)
113               ztv(ji,jj,jk) = zfp_vj * ptb(ji,jj,jk) + zfm_vj * ptb(ji  ,jj+1,jk)
114            END DO
115         END DO
116      END DO
117
118      ! upstream tracer flux in the k direction
119      ! Surface value
120      IF( lk_dynspg_rl .OR. lk_vvl ) THEN               ! rigid lid or variable volume: flux set to zero
121         ztw(:,:,1) = 0.e0
122      ELSE                                              ! free surface
123         ztw(:,:,1) = e1t(:,:) * e2t(:,:) * pwn(:,:,1) * ptb(:,:,1)
124      ENDIF
125
126      ! Interior value
127      DO jk = 2, jpkm1
128         DO jj = 1, jpj
129            DO ji = 1, jpi
130               zew = 0.5 * e1t(ji,jj) * e2t(ji,jj) * pwn(ji,jj,jk)
131               zfp_wk = zew + ABS( zew )
132               zfm_wk = zew - ABS( zew )
133               ztw(ji,jj,jk) = zfp_wk * ptb(ji,jj,jk) + zfm_wk * ptb(ji,jj,jk-1)
134            END DO
135         END DO
136      END DO
137
138      ! total advective trend
139      DO jk = 1, jpkm1
140         DO jj = 2, jpjm1
141            DO ji = fs_2, fs_jpim1   ! vector opt.
142               zbtr = 1./ ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
143               ! i- j- horizontal & k- vertical advective trends
144               ztai = - ( ztu(ji,jj,jk) - ztu(ji-1,jj  ,jk  ) ) * zbtr
145               ztaj = - ( ztv(ji,jj,jk) - ztv(ji  ,jj-1,jk  ) ) * zbtr
146               ztak = - ( ztw(ji,jj,jk) - ztw(ji  ,jj  ,jk+1) ) * zbtr
147               ! total intermediate advective trends
148               zti(ji,jj,jk) = ztai + ztaj + ztak
149            END DO
150         END DO
151      END DO
152
153      !  Save the horizontal advective trends for diagnostic
154      ! -----------------------------------------------------
155      IF( l_trdtra ) THEN
156         CALL trd_tra_adv( kt, ktra, jpt_trd_xad, cdtype, ztu, pun, ptn )
157         CALL trd_tra_adv( kt, ktra, jpt_trd_yad, cdtype, ztv, pvn, ptn )
158         CALL trd_tra_adv( kt, ktra, jpt_trd_zad, cdtype, ztw, pwn, ptn )
159      ENDIF
160
161      ! update and guess with monotonic sheme
162      DO jk = 1, jpkm1
163         z2dtt = z2 * rdttra(jk)
164         DO jj = 2, jpjm1
165            DO ji = fs_2, fs_jpim1   ! vector opt.
166               pta(ji,jj,jk) =  pta(ji,jj,jk) + zti(ji,jj,jk)
167               zti (ji,jj,jk) = ( ptb(ji,jj,jk) + z2dtt * zti(ji,jj,jk) ) * tmask(ji,jj,jk)
168            END DO
169         END DO
170      END DO
171
172      ! Lateral boundary conditions on zti   (unchanged sign)
173      CALL lbc_lnk( zti, 'T', 1. )
174
175
176      ! 3. antidiffusive flux : high order minus low order
177      ! --------------------------------------------------
178      ! antidiffusive flux on i and j
179      DO jk = 1, jpkm1
180         DO jj = 1, jpjm1
181            DO ji = 1, fs_jpim1   ! vector opt.
182               zeu = 0.5 * e2u(ji,jj) * fse3u(ji,jj,jk) * pun(ji,jj,jk)
183               zev = 0.5 * e1v(ji,jj) * fse3v(ji,jj,jk) * pvn(ji,jj,jk)
184               ztu(ji,jj,jk) = zeu * ( ptn(ji,jj,jk) + ptn(ji+1,jj,jk) ) - ztu(ji,jj,jk)
185               ztv(ji,jj,jk) = zev * ( ptn(ji,jj,jk) + ptn(ji,jj+1,jk) ) - ztv(ji,jj,jk)
186            END DO
187         END DO
188      END DO
189     
190      ! antidiffusive flux on k
191      ! Surface value
192      ztw(:,:,1) = 0.e0
193
194      ! Interior value
195      DO jk = 2, jpkm1
196         DO jj = 1, jpj
197            DO ji = 1, jpi
198               zew = 0.5 * e1t(ji,jj) * e2t(ji,jj) * pwn(ji,jj,jk)
199               ztw(ji,jj,jk) = zew * ( ptn(ji,jj,jk) + ptn(ji,jj,jk-1) ) - ztw(ji,jj,jk)
200            END DO
201         END DO
202      END DO
203
204      ! Lateral bondary conditions
205      CALL lbc_lnk( ztu, 'U', -1. )
206      CALL lbc_lnk( ztv, 'V', -1. )
207      CALL lbc_lnk( ztw, 'W',  1. )
208
209      ! 4. monotonicity algorithm
210      ! -------------------------
211      CALL nonosc( ptb, ztu, ztv, ztw, zti, z2 )
212
213
214      ! 5. final trend with corrected fluxes
215      ! ------------------------------------
216      DO jk = 1, jpkm1
217         DO jj = 2, jpjm1
218            DO ji = fs_2, fs_jpim1   ! vector opt. 
219               zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
220               ! i- j- horizontal & k- vertical advective trends
221               ztai = - ( ztu(ji,jj,jk) - ztu(ji-1,jj  ,jk  )) * zbtr
222               ztaj = - ( ztv(ji,jj,jk) - ztv(ji  ,jj-1,jk  )) * zbtr
223               ztak = - ( ztw(ji,jj,jk) - ztw(ji  ,jj  ,jk+1)) * zbtr
224
225               ! add them to the general tracer trends
226               pta(ji,jj,jk) = pta(ji,jj,jk) + ztai + ztaj + ztak
227            END DO
228         END DO
229      END DO
230
231!!gm  the transport computation is wrong, the upstream part is missing !
232      ! "zonal" mean advective heat and salt transport
233      IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN
234         IF( ktra == jp_tem)   pht_adv(:) = ptr_vj( ztv(:,:,:) )
235         IF( ktra == jp_sal)   pst_adv(:) = ptr_vj( ztv(:,:,:) )
236      ENDIF
237
238      !  Save the horizontal advective trends for diagnostic
239      ! -----------------------------------------------------
240      IF( l_trdtra ) THEN
241         CALL trd_tra_adv( kt, ktra, jpt_trd_xad, cdtype, ztu, pun, ptn, cnbpas='bis' )   ! <<< Add to iad trend
242         CALL trd_tra_adv( kt, ktra, jpt_trd_yad, cdtype, ztv, pvn, ptn, cnbpas='bis' )   ! <<< Add to jad trend
243         CALL trd_tra_adv( kt, ktra, jpt_trd_zad, cdtype, ztw, pwn, ptn, cnbpas='bis' )   ! <<< Add to zad trend
244      ENDIF
245
246      IF(ln_ctl)   CALL prt_ctl( tab3d_1=pta, clinfo1=' tvd - adv: ', mask1=tmask, clinfo3=cdtype )
247      !
248   END SUBROUTINE tra_adv_tvd
249
250
251   SUBROUTINE nonosc( pbef, paa, pbb, pcc, paft, prdt )
252      !!---------------------------------------------------------------------
253      !!                    ***  ROUTINE nonosc  ***
254      !!     
255      !! **  Purpose :   compute monotonic tracer fluxes from the upstream
256      !!       scheme and the before field by a nonoscillatory algorithm
257      !!
258      !! **  Method  :   ... ???
259      !!       warning : pbef and paft must be masked, but the boundaries
260      !!       conditions on the fluxes are not necessary zalezak (1979)
261      !!       drange (1995) multi-dimensional forward-in-time and upstream-
262      !!       in-space based differencing for fluid
263      !!----------------------------------------------------------------------
264      REAL(wp), INTENT(in   ) ::   prdt                               ! ???
265      REAL(wp), INTENT(inout), DIMENSION (jpi,jpj,jpk) ::   pbef, paft       ! before & after field
266      REAL(wp), INTENT(inout), DIMENSION (jpi,jpj,jpk) ::   paa, pbb, pcc    ! monotonic flux in the 3 directions
267      !!
268      INTEGER ::   ji, jj, jk               ! dummy loop indices
269      INTEGER ::   ikm1
270      REAL(wp), DIMENSION (jpi,jpj,jpk) ::   zbetup, zbetdo
271      REAL(wp) ::   zpos, zneg, zbt, za, zb, zc, zbig, zrtrn, z2dtt
272      !!----------------------------------------------------------------------
273
274      zbig = 1.e+40
275      zrtrn = 1.e-15
276      zbetup(:,:,:) = 0.e0   ;   zbetdo(:,:,:) = 0.e0 
277
278      ! Search local extrema
279      ! --------------------
280      ! large negative value (-zbig) inside land
281      pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:) - zbig * ( 1.e0 - tmask(:,:,:) )
282      paft(:,:,:) = paft(:,:,:) * tmask(:,:,:) - zbig * ( 1.e0 - tmask(:,:,:) )
283      ! search maximum in neighbourhood
284      DO jk = 1, jpkm1
285         ikm1 = MAX(jk-1,1)
286         DO jj = 2, jpjm1
287            DO ji = fs_2, fs_jpim1   ! vector opt.
288               zbetup(ji,jj,jk) = MAX(  pbef(ji  ,jj  ,jk  ), paft(ji  ,jj  ,jk  ),   &
289                  &                     pbef(ji-1,jj  ,jk  ), pbef(ji+1,jj  ,jk  ),   &
290                  &                     paft(ji-1,jj  ,jk  ), paft(ji+1,jj  ,jk  ),   &
291                  &                     pbef(ji  ,jj-1,jk  ), pbef(ji  ,jj+1,jk  ),   &
292                  &                     paft(ji  ,jj-1,jk  ), paft(ji  ,jj+1,jk  ),   &
293                  &                     pbef(ji  ,jj  ,ikm1), pbef(ji  ,jj  ,jk+1),   &
294                  &                     paft(ji  ,jj  ,ikm1), paft(ji  ,jj  ,jk+1)  )
295            END DO
296         END DO
297      END DO
298      ! large positive value (+zbig) inside land
299      pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:) + zbig * ( 1.e0 - tmask(:,:,:) )
300      paft(:,:,:) = paft(:,:,:) * tmask(:,:,:) + zbig * ( 1.e0 - tmask(:,:,:) )
301      ! search minimum in neighbourhood
302      DO jk = 1, jpkm1
303         ikm1 = MAX(jk-1,1)
304         DO jj = 2, jpjm1
305            DO ji = fs_2, fs_jpim1   ! vector opt.
306               zbetdo(ji,jj,jk) = MIN(  pbef(ji  ,jj  ,jk  ), paft(ji  ,jj  ,jk  ),   &
307                  &                     pbef(ji-1,jj  ,jk  ), pbef(ji+1,jj  ,jk  ),   &
308                  &                     paft(ji-1,jj  ,jk  ), paft(ji+1,jj  ,jk  ),   &
309                  &                     pbef(ji  ,jj-1,jk  ), pbef(ji  ,jj+1,jk  ),   &
310                  &                     paft(ji  ,jj-1,jk  ), paft(ji  ,jj+1,jk  ),   &
311                  &                     pbef(ji  ,jj  ,ikm1), pbef(ji  ,jj  ,jk+1),   &
312                  &                     paft(ji  ,jj  ,ikm1), paft(ji  ,jj  ,jk+1)  )
313            END DO
314         END DO
315      END DO
316
317      ! restore masked values to zero
318      pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:)
319      paft(:,:,:) = paft(:,:,:) * tmask(:,:,:)
320 
321
322      ! 2. Positive and negative part of fluxes and beta terms
323      ! ------------------------------------------------------
324
325      DO jk = 1, jpkm1
326         z2dtt = prdt * rdttra(jk)
327         DO jj = 2, jpjm1
328            DO ji = fs_2, fs_jpim1   ! vector opt.
329               ! positive & negative part of the flux
330               zpos = MAX( 0., paa(ji-1,jj  ,jk  ) ) - MIN( 0., paa(ji  ,jj  ,jk  ) )   &
331                  & + MAX( 0., pbb(ji  ,jj-1,jk  ) ) - MIN( 0., pbb(ji  ,jj  ,jk  ) )   &
332                  & + MAX( 0., pcc(ji  ,jj  ,jk+1) ) - MIN( 0., pcc(ji  ,jj  ,jk  ) )
333               zneg = MAX( 0., paa(ji  ,jj  ,jk  ) ) - MIN( 0., paa(ji-1,jj  ,jk  ) )   &
334                  & + MAX( 0., pbb(ji  ,jj  ,jk  ) ) - MIN( 0., pbb(ji  ,jj-1,jk  ) )   &
335                  & + MAX( 0., pcc(ji  ,jj  ,jk  ) ) - MIN( 0., pcc(ji  ,jj  ,jk+1) )
336               ! up & down beta terms
337               zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) / z2dtt
338               zbetup(ji,jj,jk) = ( zbetup(ji,jj,jk) - paft(ji,jj,jk) ) / (zpos+zrtrn) * zbt
339               zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zbetdo(ji,jj,jk) ) / (zneg+zrtrn) * zbt
340            END DO
341         END DO
342      END DO
343
344      ! lateral boundary condition on zbetup & zbetdo   (unchanged sign)
345      CALL lbc_lnk( zbetup, 'T', 1. )
346      CALL lbc_lnk( zbetdo, 'T', 1. )
347
348
349      ! 3. monotonic flux in the i & j direction (paa & pbb)
350      ! ----------------------------------------
351      DO jk = 1, jpkm1
352         DO jj = 2, jpjm1
353            DO ji = fs_2, fs_jpim1   ! vector opt.
354               za = MIN( 1.e0, zbetdo(ji,jj,jk), zbetup(ji+1,jj,jk) )
355               zb = MIN( 1.e0, zbetup(ji,jj,jk), zbetdo(ji+1,jj,jk) )
356               zc = 0.5 * ( 1.e0 + SIGN( 1.e0, paa(ji,jj,jk) ) )
357               paa(ji,jj,jk) = paa(ji,jj,jk) * ( zc * za + ( 1.e0 - zc) * zb )
358
359               za = MIN( 1.e0, zbetdo(ji,jj,jk), zbetup(ji,jj+1,jk) )
360               zb = MIN( 1.e0, zbetup(ji,jj,jk), zbetdo(ji,jj+1,jk) )
361               zc = 0.5 * ( 1.e0 + SIGN( 1.e0, pbb(ji,jj,jk) ) )
362               pbb(ji,jj,jk) = pbb(ji,jj,jk) * ( zc * za + ( 1.e0 - zc) * zb )
363            END DO
364         END DO
365      END DO
366
367
368      ! monotonic flux in the k direction, i.e. pcc
369      ! -------------------------------------------
370      DO jk = 2, jpkm1
371         DO jj = 2, jpjm1
372            DO ji = fs_2, fs_jpim1   ! vector opt.
373
374               za = MIN( 1., zbetdo(ji,jj,jk), zbetup(ji,jj,jk-1) )
375               zb = MIN( 1., zbetup(ji,jj,jk), zbetdo(ji,jj,jk-1) )
376               zc = 0.5 * ( 1.e0 + SIGN( 1.e0, pcc(ji,jj,jk) ) )
377               pcc(ji,jj,jk) = pcc(ji,jj,jk) * ( zc * za + ( 1.e0 - zc) * zb )
378            END DO
379         END DO
380      END DO
381
382      ! lateral boundary condition on paa, pbb, pcc
383      CALL lbc_lnk( paa, 'U', -1. )      ! changed sign
384      CALL lbc_lnk( pbb, 'V', -1. )      ! changed sign
385      CALL lbc_lnk( pcc, 'W',  1. )      ! NO changed sign
386      !
387   END SUBROUTINE nonosc
388
389   !!======================================================================
390END MODULE traadv_tvd
Note: See TracBrowser for help on using the repository browser.