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

source: trunk/NEMO/OPA_SRC/TRA/traadv_cen2_jki.F90 @ 699

Last change on this file since 699 was 699, checked in by smasson, 17 years ago

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