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

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