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_cen2.F90 in trunk/NEMO/OPA_SRC/TRA – NEMO

source: trunk/NEMO/OPA_SRC/TRA/traadv_cen2.F90 @ 216

Last change on this file since 216 was 216, checked in by opalod, 19 years ago

CT : UPDATE151 : New trends organization
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 18.0 KB
Line 
1MODULE traadv_cen2
2   !!==============================================================================
3   !!                       ***  MODULE  traadv_cen2  ***
4   !! Ocean active tracers:  horizontal & vertical advective trend
5   !!==============================================================================
6 
7   !!----------------------------------------------------------------------
8   !!   tra_adv_cen2 : update the tracer trend with the horizontal
9   !!                  and vertical advection trends using a 2nd order
10   !!                  centered finite difference scheme
11   !!----------------------------------------------------------------------
12   !! * Modules used
13   USE oce             ! ocean dynamics and active tracers
14   USE dom_oce         ! ocean space and time domain
15   USE trdmod          ! ocean active tracers trends
16   USE trdmod_oce      ! ocean variables trends
17   USE flxrnf          !
18   USE trabbl          ! advective term in the BBL
19   USE ocfzpt          !
20   USE lib_mpp
21   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
22   USE in_out_manager  ! I/O manager
23   USE diaptr          ! poleward transport diagnostics
24   USE dynspg_fsc      !
25   USE dynspg_fsc_atsk !
26
27   IMPLICIT NONE
28   PRIVATE
29
30   !! * Accessibility
31   PUBLIC tra_adv_cen2    ! routine called by step.F90
32
33   !! * Substitutions
34#  include "domzgr_substitute.h90"
35#  include "vectopt_loop_substitute.h90"
36   !!----------------------------------------------------------------------
37   !!   OPA 9.0 , LODYC-IPSL (2003)
38   !!----------------------------------------------------------------------
39
40CONTAINS
41
42#if defined key_autotasking
43   !!----------------------------------------------------------------------
44   !!   'key_autotasking' :      2nd order centered scheme (k- and j-slabs)
45   !!----------------------------------------------------------------------
46#  include "traadv_cen2_atsk.h90"
47
48#else
49   !!----------------------------------------------------------------------
50   !!   Default option :             2nd order centered scheme (k-j-i loop)
51   !!----------------------------------------------------------------------
52
53   SUBROUTINE tra_adv_cen2( kt )
54      !!----------------------------------------------------------------------
55      !!                  ***  ROUTINE tra_adv_cen2  ***
56      !!                 
57      !! ** Purpose :   Compute the now trend due to the advection of tracers
58      !!      and add it to the general trend of passive tracer equations.
59      !!
60      !! ** Method  :   The advection is evaluated by a second order centered
61      !!      scheme using now fields (leap-frog scheme). In specific areas
62      !!      (vicinity of major river mouths, some straits, or where tn is
63      !!      is approaching the freezing point) it is mixed with an upstream
64      !!      scheme for stability reasons.
65      !!        Part 0 : compute the upstream / centered flag
66      !!                 (3D array, zind, defined at T-point (0<zind<1))
67      !!        Part I : horizontal advection
68      !!      * centered flux:
69      !!         * s-coordinate (lk_sco=T) or
70      !!         * z-coordinate with partial steps (lk_zps=T),
71      !!        the vertical scale factors e3. are inside the derivatives:
72      !!               zcenu = e2u*e3u  un  mi(tn)
73      !!               zcenv = e1v*e3v  vn  mj(tn)
74      !!         * z-coordinate (default key), e3t=e3u=e3v:
75      !!               zcenu = e2u  un  mi(tn)
76      !!               zcenv = e1v  vn  mj(tn)
77      !!      * upstream flux:
78      !!         * s-coordinate (lk_sco=T) or
79      !!         * z-coordinate with partial steps (lk_zps=T)
80      !!               zupsu = e2u*e3u  un  (tb(i) or tb(i-1) ) [un>0 or <0]
81      !!               zupsv = e1v*e3v  vn  (tb(j) or tb(j-1) ) [vn>0 or <0]
82      !!         * z-coordinate (default key)
83      !!               zupsu = e2u*e3u  un  (tb(i) or tb(i-1) ) [un>0 or <0]
84      !!               zupsv = e1v*e3v  vn  (tb(j) or tb(j-1) ) [vn>0 or <0]
85      !!      * mixed upstream / centered horizontal advection scheme
86      !!               zcofi = max(zind(i+1), zind(i))
87      !!               zcofj = max(zind(j+1), zind(j))
88      !!               zwx = zcofi * zupsu + (1-zcofi) * zcenu
89      !!               zwy = zcofj * zupsv + (1-zcofj) * zcenv
90      !!      * horizontal advective trend (divergence of the fluxes)
91      !!         * s-coordinate (lk_sco=T) or
92      !!         * z-coordinate with partial steps (lk_zps=T)
93      !!               zta = 1/(e1t*e2t*e3t) { di-1[zwx] + dj-1[zwy] }
94      !!         * z-coordinate (default key), e3t=e3u=e3v:
95      !!               zta = 1/(e1t*e2t) { di-1[zwx] + dj-1[zwy] }
96      !!      * Add this trend now to the general trend of tracer (ta,sa):
97      !!              (ta,sa) = (ta,sa) + ( zta , zsa )
98      !!      * trend diagnostic ('key_trdtra'): the trend is saved
99      !!      for diagnostics. The trends saved is expressed as
100      !!      Uh.gradh(T), (save trend = zta + tn divn).
101      !!         In addition, the advective trend in the two horizontal direc-
102      !!      tion is also re-computed as Uh gradh(T). Indeed hadt+tn divn is
103      !!      equal to (in s-coordinates, and similarly in z-coord.):
104      !!         zta+tn*divn=1/(e1t*e2t*e3t) { mi-1( e2u*e3u  un  di[tn] )
105      !!                                      +mj-1( e1v*e3v  vn  mj[tn] )  }
106      !!         C A U T I O N : the trend saved is the centered trend only.
107      !!      It doesn't take into account the upstream part of the scheme.
108      !!
109      !!         Part II : vertical advection
110      !!      For temperature (idem for salinity) the advective trend is com-
111      !!      puted as follows :
112      !!            zta = 1/e3t dk+1[ zwz ]
113      !!      where the vertical advective flux, zwz, is given by :
114      !!            zwz = zcofk * zupst + (1-zcofk) * zcent
115      !!      with
116      !!        zupsv = upstream flux = wn * (tb(k) or tb(k-1) ) [wn>0 or <0]
117      !!        zcenu = centered flux = wn * mk(tn)
118      !!         The surface boundary condition is :
119      !!      rigid-lid (default option) : zero advective flux
120      !!      free-surf ("key_fresurf_cstvol") : wn(:,:,1) * tn(:,:,1)
121      !!         Add this trend now to the general trend of tracer (ta,sa):
122      !!            (ta,sa) = (ta,sa) + ( zta , zsa )
123      !!         Trend diagnostic ('key_trdtra'): the trend is saved for
124      !!      diagnostics. The trends saved is expressed as :
125      !!             save trend =  w.gradz(T) = zta - tn divn.
126      !!
127      !! ** Action : - update (ta,sa) with the now advective tracer trends
128      !!             - save the trends in (ttrdh,strdh) ('key_trdtra')
129      !!
130      !! History :
131      !!   8.2  !  01-08  (G. Madec, E. Durand)  trahad+trazad = traadv
132      !!   8.5  !  02-06  (G. Madec)  F90: Free form and module
133      !!   9.0  !  04-08  (C. Talandier) New trends organization
134      !!----------------------------------------------------------------------
135      !! * Modules used
136      USE oce                , zwx => ua,  &  ! use ua as workspace
137         &                     zwy => va      ! use va as workspace
138#if defined key_trabbl_adv
139      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  &  ! temporary arrays
140         &         zun, zvn, zwn
141#else
142      USE oce                , zun => un,  &  ! When no bbl, zun == un
143         &                     zvn => vn,  &  ! When no bbl, zvn == vn
144         &                     zwn => wn      ! When no bbl, zwn == wn
145#endif
146
147 
148      !! * Arguments
149      INTEGER, INTENT( in ) ::   kt           ! ocean time-step index
150 
151      !! * Local save
152      REAL(wp), DIMENSION(jpi,jpj), SAVE ::   &
153         zbtr2
154 
155      !! * Local declarations
156      INTEGER  ::   ji, jj, jk                 ! dummy loop indices
157      REAL(wp) ::                           &
158         zbtr, zta, zsa, zfui, zfvj,        &  ! temporary scalars
159         zhw, ze3tr, zcofi, zcofj,          &  !    "         "
160         zupsut, zupsvt, zupsus, zupsvs,    &  !    "         "
161         zfp_ui, zfp_vj, zfm_ui, zfm_vj,    &  !    "         "
162         zcofk, zupst, zupss, zcent,        &  !    "         "
163         zcens, zfp_w, zfm_w,               &  !    "         "
164         zcenut, zcenvt, zcenus, zcenvs        !    "         "
165      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &
166         zwz, zww, zind,                    &  ! temporary workspace arrays
167         ztdta, ztdsa                          !    "         "
168      !!----------------------------------------------------------------------
169
170      IF( kt == nit000 ) THEN
171         IF(lwp) WRITE(numout,*)
172         IF(lwp) WRITE(numout,*) 'tra_adv_cen2 : 2nd order centered advection scheme'
173         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~   Vector optimization case'
174         IF(lwp) WRITE(numout,*)
175   
176         zbtr2(:,:) = 1. / ( e1t(:,:) * e2t(:,:) )
177      ENDIF
178
179      ! Save ta and sa trends
180      IF( l_trdtra )   THEN
181         ztdta(:,:,:) = ta(:,:,:) 
182         ztdsa(:,:,:) = sa(:,:,:) 
183         l_adv = 'ce2'
184      ENDIF
185
186#if defined key_trabbl_adv
187      ! Advective bottom boundary layer
188      ! -------------------------------
189      zun(:,:,:) = un(:,:,:) - u_bbl(:,:,:)
190      zvn(:,:,:) = vn(:,:,:) - v_bbl(:,:,:)
191      zwn(:,:,:) = wn(:,:,:) + w_bbl(:,:,:)
192#endif
193
194      ! Upstream / centered scheme indicator
195      ! ------------------------------------
196 
197      DO jk = 1, jpk
198         DO jj = 1, jpj
199            DO ji = 1, jpi
200               zind(ji,jj,jk) =  MAX ( upsrnfh(ji,jj) * upsrnfz(jk),     &  ! changing advection scheme near runoff
201                  &                    upsadv(ji,jj)                     &  ! in the vicinity of some straits
202#if defined key_ice_lim
203                  &                  , tmask(ji,jj,jk)                   &  ! half upstream tracer fluxes
204                  &                  * MAX( 0., SIGN( 1., fzptn(ji,jj)   &  ! if tn < ("freezing"+0.1 )
205                  &                                +0.1-tn(ji,jj,jk) ) ) &
206#endif
207                  &                  )
208            END DO
209         END DO
210      END DO
211
212
213      ! I. Horizontal advective fluxes
214      ! ------------------------------
215
216      ! 1. Second order centered tracer flux at u and v-points
217      ! -------------------------------------------------------
218
219      !                                                ! ===============
220      DO jk = 1, jpkm1                                 ! Horizontal slab
221         !                                             ! ===============
222         DO jj = 1, jpjm1
223            DO ji = 1, fs_jpim1   ! vector opt.
224               ! upstream indicator
225               zcofi = MAX( zind(ji+1,jj,jk), zind(ji,jj,jk) )
226               zcofj = MAX( zind(ji,jj+1,jk), zind(ji,jj,jk) )
227               ! volume fluxes * 1/2
228#if defined key_s_coord || defined key_partial_steps
229               zfui = 0.5 * e2u(ji,jj) * fse3u(ji,jj,jk) * zun(ji,jj,jk)
230               zfvj = 0.5 * e1v(ji,jj) * fse3v(ji,jj,jk) * zvn(ji,jj,jk)
231#else
232               zfui = 0.5 * e2u(ji,jj) * zun(ji,jj,jk)
233               zfvj = 0.5 * e1v(ji,jj) * zvn(ji,jj,jk)
234#endif
235               ! upstream scheme
236               zfp_ui = zfui + ABS( zfui )
237               zfp_vj = zfvj + ABS( zfvj )
238               zfm_ui = zfui - ABS( zfui )
239               zfm_vj = zfvj - ABS( zfvj )
240               zupsut = zfp_ui * tb(ji,jj,jk) + zfm_ui * tb(ji+1,jj  ,jk)
241               zupsvt = zfp_vj * tb(ji,jj,jk) + zfm_vj * tb(ji  ,jj+1,jk)
242               zupsus = zfp_ui * sb(ji,jj,jk) + zfm_ui * sb(ji+1,jj  ,jk)
243               zupsvs = zfp_vj * sb(ji,jj,jk) + zfm_vj * sb(ji  ,jj+1,jk)
244               ! centered scheme
245               zcenut = zfui * ( tn(ji,jj,jk) + tn(ji+1,jj  ,jk) )
246               zcenvt = zfvj * ( tn(ji,jj,jk) + tn(ji  ,jj+1,jk) )
247               zcenus = zfui * ( sn(ji,jj,jk) + sn(ji+1,jj  ,jk) )
248               zcenvs = zfvj * ( sn(ji,jj,jk) + sn(ji  ,jj+1,jk) )
249               ! mixed centered / upstream scheme
250               zwx(ji,jj,jk) = zcofi * zupsut + (1.-zcofi) * zcenut
251               zwy(ji,jj,jk) = zcofj * zupsvt + (1.-zcofj) * zcenvt
252               zww(ji,jj,jk) = zcofi * zupsus + (1.-zcofi) * zcenus
253               zwz(ji,jj,jk) = zcofj * zupsvs + (1.-zcofj) * zcenvs
254            END DO
255         END DO
256
257         ! 2. Tracer flux divergence at t-point added to the general trend
258         ! ---------------------------------------------------------------
259
260         DO jj = 2, jpjm1
261            DO ji = fs_2, fs_jpim1   ! vector opt.
262#if defined key_s_coord || defined key_partial_steps
263               zbtr = zbtr2(ji,jj) / fse3t(ji,jj,jk)
264#else
265               zbtr = zbtr2(ji,jj) 
266#endif
267               ! horizontal advective trends
268               zta = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk)   &
269                  &            + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk)  )
270               zsa = - zbtr * (  zww(ji,jj,jk) - zww(ji-1,jj  ,jk)   &
271                  &            + zwz(ji,jj,jk) - zwz(ji  ,jj-1,jk)  )
272               ! add it to the general tracer trends
273               ta(ji,jj,jk) = ta(ji,jj,jk) + zta
274               sa(ji,jj,jk) = sa(ji,jj,jk) + zsa
275            END DO
276         END DO
277         !                                             ! ===============
278      END DO                                           !   End of slab
279      !                                                ! ===============
280
281      ! 3. Save the horizontal advective trends for diagnostic
282      ! ------------------------------------------------------
283      IF( l_trdtra )   THEN
284         ! Recompute the hoizontal advection zta & zsa trends computed
285         ! at the step 2. above in making the difference between the new
286         ! trends and the previous one ta()/sa - ztdta()/ztdsa() and add
287         ! the term tn()/sn()*hdivn() to recover the Uh gradh(T/S) trends
288         ztdta(:,:,:) = ta(:,:,:) - ztdta(:,:,:) + tn(:,:,:) * hdivn(:,:,:) 
289         ztdsa(:,:,:) = sa(:,:,:) - ztdsa(:,:,:) + sn(:,:,:) * hdivn(:,:,:)
290
291         CALL trd_mod(ztdta, ztdsa, jpttdlad, 'TRA', kt)
292
293         ! Save the new ta and sa trends
294         ztdta(:,:,:) = ta(:,:,:) 
295         ztdsa(:,:,:) = sa(:,:,:) 
296
297      ENDIF
298
299      IF(l_ctl) THEN
300         zta = SUM( ta(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
301         zsa = SUM( sa(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
302         WRITE(numout,*) ' had  - Ta: ', zta-t_ctl, ' Sa: ', zsa-s_ctl
303         t_ctl = zta   ;   s_ctl = zsa
304      ENDIF
305
306      ! "zonal" mean advective heat and salt transport
307      IF( ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN
308# if defined key_s_coord || defined key_partial_steps
309         pht_adv(:) = ptr_vj( zwy(:,:,:) )
310         pst_adv(:) = ptr_vj( zwz(:,:,:) )
311# else
312         DO jk = 1, jpkm1
313            DO jj = 2, jpjm1
314               DO ji = fs_2, fs_jpim1   ! vector opt.
315                  zwy(ji,jj,jk) = zwy(ji,jj,jk) * fse3v(ji,jj,jk)
316                  zwz(ji,jj,jk) = zwz(ji,jj,jk) * fse3v(ji,jj,jk)
317               END DO
318            END DO
319         END DO
320         pht_adv(:) = ptr_vj( zwy(:,:,:) )
321         pst_adv(:) = ptr_vj( zwz(:,:,:) )
322# endif
323      ENDIF
324
325
326      ! II. Vertical advection
327      ! ----------------------
328
329      ! Bottom value : flux set to zero
330      zwx(:,:,jpk) = 0.e0     ;    zwy(:,:,jpk) = 0.e0
331
332      ! Surface value
333      IF( lk_dynspg_fsc ) THEN
334         ! free surface-constant volume
335         zwx(:,:, 1 ) = zwn(:,:,1) * tn(:,:,1)
336         zwy(:,:, 1 ) = zwn(:,:,1) * sn(:,:,1)
337      ELSE
338         ! rigid lid : flux set to zero
339         zwx(:,:, 1 ) = 0.e0    ;    zwy(:,:, 1 ) = 0.e0
340      ENDIF
341
342      ! 1. Vertical advective fluxes
343      ! ----------------------------
344
345      ! Second order centered tracer flux at w-point
346
347      DO jk = 2, jpk
348         DO jj = 2, jpjm1
349            DO ji = fs_2, fs_jpim1   ! vector opt.
350               ! upstream indicator
351               zcofk = MAX( zind(ji,jj,jk-1), zind(ji,jj,jk) )
352               ! velocity * 1/2
353               zhw = 0.5 * zwn(ji,jj,jk)
354               ! upstream scheme
355               zfp_w = zhw + ABS( zhw )
356               zfm_w = zhw - ABS( zhw )
357               zupst = zfp_w * tb(ji,jj,jk) + zfm_w * tb(ji,jj,jk-1)
358               zupss = zfp_w * sb(ji,jj,jk) + zfm_w * sb(ji,jj,jk-1)
359               ! centered scheme
360               zcent = zhw * ( tn(ji,jj,jk) + tn(ji,jj,jk-1) )
361               zcens = zhw * ( sn(ji,jj,jk) + sn(ji,jj,jk-1) )
362               ! mixed centered / upstream scheme
363               zwx(ji,jj,jk) = zcofk * zupst + (1.-zcofk) * zcent
364               zwy(ji,jj,jk) = zcofk * zupss + (1.-zcofk) * zcens
365            END DO
366         END DO
367      END DO
368
369
370      ! 2. Tracer flux divergence at t-point added to the general trend
371      ! -------------------------
372 
373      DO jk = 1, jpkm1
374         DO jj = 2, jpjm1
375            DO ji = fs_2, fs_jpim1   ! vector opt.
376               ze3tr = 1. / fse3t(ji,jj,jk)
377               ! vertical advective trends
378               zta = - ze3tr * ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) )
379               zsa = - ze3tr * ( zwy(ji,jj,jk) - zwy(ji,jj,jk+1) )
380               ! add it to the general tracer trends
381               ta(ji,jj,jk) =  ta(ji,jj,jk) + zta
382               sa(ji,jj,jk) =  sa(ji,jj,jk) + zsa
383            END DO
384         END DO
385      END DO
386
387      ! 3. Save the vertical advective trends for diagnostic
388      ! ----------------------------------------------------
389      IF( l_trdtra )   THEN
390         ! Recompute the vertical advection zta & zsa trends computed
391         ! at the step 2. above in making the difference between the new
392         ! trends and the previous one: ta()/sa - ztdta()/ztdsa() and substract
393         ! the term tn()/sn()*hdivn() to recover the W gradz(T/S) trends
394         ztdta(:,:,:) = ta(:,:,:) - ztdta(:,:,:) - tn(:,:,:) * hdivn(:,:,:) 
395         ztdsa(:,:,:) = sa(:,:,:) - ztdsa(:,:,:) - sn(:,:,:) * hdivn(:,:,:)
396
397         CALL trd_mod(ztdta, ztdsa, jpttdzad, 'TRA', kt)
398      ENDIF
399
400      IF(l_ctl) THEN
401         zta = SUM( ta(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
402         zsa = SUM( sa(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
403         WRITE(numout,*) ' zad  - Ta: ', zta-t_ctl, ' Sa: ', zsa-s_ctl, ' centered2'
404         t_ctl = zta   ;   s_ctl = zsa
405      ENDIF
406
407   END SUBROUTINE tra_adv_cen2
408
409#endif
410
411   !!======================================================================
412END MODULE traadv_cen2
Note: See TracBrowser for help on using the repository browser.