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

Last change on this file since 881 was 833, checked in by rblod, 16 years ago

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