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 branches/dev_001_GM/NEMO/OPA_SRC/TRA – NEMO

source: branches/dev_001_GM/NEMO/OPA_SRC/TRA/traadv_cen2.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: 13.9 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   !!   NEMO     1.0  !  05-11  (V. Garnier) Surface pressure gradient organization
9   !!             -   !  05-11  (V. Garnier) Surface pressure gradient organization
10   !!             -   !  06-04  (R. Benshila, G. Madec) Step reorganization
11   !!            2.4  !  08-01  (G. Madec) Merge TRA-TRC
12   !!----------------------------------------------------------------------
13
14   !!----------------------------------------------------------------------
15   !!   tra_adv_cen2 : update the tracer trend with the horizontal and vertical
16   !!                  advection trends using a 2nd order centered scheme
17   !!----------------------------------------------------------------------
18   USE oce, ONLY: tn   ! now ocean temperature
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 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   !! NEMO/OPA 2.4 , LOCEAN-IPSL (2008)
43   !! $Id:$
44   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
45   !!----------------------------------------------------------------------
46
47CONTAINS
48
49   SUBROUTINE tra_adv_cen2( kt, cdtype, ktra, pun, pvn, pwn,   &
50      &                                       ptb, ptn, pta )
51      !!----------------------------------------------------------------------
52      !!                  ***  ROUTINE tra_adv_cen2  ***
53      !!                 
54      !! ** Purpose :   Compute the now trend due to the advection of a tracer
55      !!      and add it to the corresponding general trend of 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(ptn)
67      !!               zcenv = e1v*e3v  vn  mj(ptn)
68      !!       * upstream flux:
69      !!               zupsu = e2u*e3u  un  (ptb(i) or ptb(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 (pta):
79      !!               pta = pta + zta
80      !!       * trend diagnostic (lk_trdtra=T): the trend is
81      !!      saved for diagnostics. The trends saved is expressed as
82      !!      Uh.gradh(T), i.e.
83      !!                     save trend = zta + ptn divn
84      !!         In addition, the advective trend in the two horizontal direc-
85      !!      tion is also re-computed as Uh gradh(T). Indeed hadt+ptn divn is
86      !!      equal to (in s-coordinates, and similarly in z-coord.):
87      !!         zta+ptn*divn=1/(e1t*e2t*e3t) { mi-1( e2u*e3u  un  di[ptn] )
88      !!                                       +mj-1( e1v*e3v  vn  mj[ptn] )  }
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      !!      the advective trend is computed as follows :
94      !!            zta = 1/e3t dk+1[ zwx ]
95      !!      where the vertical advective flux, zwx, is given by :
96      !!            zwx = 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(ptn)
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) * ptn(:,:,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 - ptn divn.
108      !!
109      !! ** Action :  - update (ta,sa) with the now advective tracer trends
110      !!              - trend diagnostics (lk_trdtra=T)
111      !!----------------------------------------------------------------------
112      INTEGER         , INTENT(in   )                         ::   kt              ! ocean time-step index
113      CHARACTER(len=3), INTENT(in   )                         ::   cdtype          ! =TRA or TRC (tracer indicator)
114      INTEGER         , INTENT(in   )                         ::   ktra            ! tracer index
115      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   pun, pvn, pwn   ! 3 ocean velocity components
116      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   ptb, ptn        ! before and now tracer fields
117      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pta             ! tracer trend
118      !!
119      INTEGER  ::   ji, jj, jk                                    ! dummy loop indices
120      REAL(wp) ::   zbtr, zta, zhw, ze3tr                         ! temporary scalars
121      REAL(wp) ::   zcofi, zfui, zcenut, zupsut, zfp_ui, zfm_ui   !    "         "
122      REAL(wp) ::   zcofj, zfvj, zcenvt, zupsvt, zfp_vj, zfm_vj   !    "         "
123      REAL(wp) ::   zcofk,       zcent , zupst , zfp_w , zfm_w    !    "         "
124      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwx, zwy, zind   ! 3D workspace
125      !!----------------------------------------------------------------------
126
127      IF( kt == nit000 ) THEN
128         IF(lwp) WRITE(numout,*)
129         IF(lwp) WRITE(numout,*) 'tra_adv_cen2 : 2nd order centered advection scheme'
130         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
131         !
132         btr2(:,:) = 1. / ( e1t(:,:) * e2t(:,:) )
133      ENDIF
134
135      ! Upstream / centered scheme indicator
136      ! ------------------------------------
137      DO jk = 1, jpk
138         DO jj = 1, jpj
139            DO ji = 1, jpi
140               zind(ji,jj,jk) =  MAX ( upsrnfh(ji,jj) * upsrnfz(jk),     &  ! changing advection scheme near runoff
141                  &                    upsadv(ji,jj)                     &  ! in the vicinity of some straits
142#if defined key_ice_lim
143                  &                  , tmask(ji,jj,jk)                   &  ! half upstream tracer fluxes
144                  &                  * MAX( 0., SIGN( 1., fzptn(ji,jj)   &  ! if tn < ("freezing"+0.1 )
145                  &                                +0.1-tn(ji,jj,jk) ) ) &
146#endif
147                  &                  )
148            END DO
149         END DO
150      END DO
151
152      ! I. Horizontal advection
153      ! -----------------------
154
155      !                                                ! ===============
156      DO jk = 1, jpkm1                                 ! Horizontal slab
157         !                                             ! ===============
158         !
159         DO jj = 1, jpjm1                              !  Horizontal advective fluxes
160            DO ji = 1, fs_jpim1   ! vector opt.
161               ! upstream indicator
162               zcofi = MAX( zind(ji+1,jj,jk), zind(ji,jj,jk) )
163               zcofj = MAX( zind(ji,jj+1,jk), zind(ji,jj,jk) )
164               ! volume fluxes * 1/2
165#if defined key_zco
166               zfui = 0.5 * e2u(ji,jj) * pun(ji,jj,jk)
167               zfvj = 0.5 * e1v(ji,jj) * pvn(ji,jj,jk)
168#else
169               zfui = 0.5 * e2u(ji,jj) * fse3u(ji,jj,jk) * pun(ji,jj,jk)
170               zfvj = 0.5 * e1v(ji,jj) * fse3v(ji,jj,jk) * pvn(ji,jj,jk)
171#endif
172               ! upstream scheme
173               zfp_ui = zfui + ABS( zfui )
174               zfp_vj = zfvj + ABS( zfvj )
175               zfm_ui = zfui - ABS( zfui )
176               zfm_vj = zfvj - ABS( zfvj )
177               zupsut = zfp_ui * ptb(ji,jj,jk) + zfm_ui * ptb(ji+1,jj  ,jk)
178               zupsvt = zfp_vj * ptb(ji,jj,jk) + zfm_vj * ptb(ji  ,jj+1,jk)
179               ! centered scheme
180               zcenut = zfui * ( ptn(ji,jj,jk) + ptn(ji+1,jj  ,jk) )
181               zcenvt = zfvj * ( ptn(ji,jj,jk) + ptn(ji  ,jj+1,jk) )
182               ! mixed centered / upstream scheme
183               zwx(ji,jj,jk) = zcofi * zupsut + (1.-zcofi) * zcenut
184               zwy(ji,jj,jk) = zcofj * zupsvt + (1.-zcofj) * zcenvt
185            END DO
186         END DO
187
188         DO jj = 2, jpjm1                              !  horizontal tracer flux divergence added to the general trend
189            DO ji = fs_2, fs_jpim1   ! vector opt.
190#if defined key_zco
191               zbtr = btr2(ji,jj)
192#else
193               zbtr = btr2(ji,jj) / fse3t(ji,jj,jk)
194#endif
195               ! horizontal advective trends
196               zta = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk)   &
197                  &            + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk)  )
198               ! add it to the general tracer trends
199               pta(ji,jj,jk) = pta(ji,jj,jk) + zta
200            END DO
201         END DO
202         !                                             ! ===============
203      END DO                                           !   End of slab
204      !                                                ! ===============
205
206      !  Save the horizontal advective trends for diagnostic
207      ! -----------------------------------------------------
208      IF( l_trdtra ) THEN
209         CALL trd_tra_adv( kt, ktra, jpt_trd_xad, cdtype, zwx, pun, ptn )
210         CALL trd_tra_adv( kt, ktra, jpt_trd_yad, cdtype, zwy, pvn, ptn )
211      ENDIF
212
213      IF(ln_ctl)   CALL prt_ctl( tab3d_1=pta, clinfo1=' cen2 - had: ', mask1=tmask, clinfo3=cdtype )
214
215
216      ! "Poleward" heat and salt transport
217      ! ----------------------------------
218      IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN
219         IF( lk_zco ) THEN
220            DO jk = 1, jpkm1
221               DO jj = 2, jpjm1
222                  DO ji = fs_2, fs_jpim1   ! vector opt.
223                    zwy(ji,jj,jk) = zwy(ji,jj,jk) * fse3v(ji,jj,jk)
224                  END DO
225               END DO
226            END DO
227         ENDIF
228         IF( ktra == jp_tem)   pht_adv(:) = ptr_vj( zwy(:,:,:) )
229         IF( ktra == jp_sal)   pst_adv(:) = ptr_vj( zwy(:,:,:) )
230      ENDIF
231
232
233      ! II. Vertical advection
234      ! ----------------------
235
236      IF( lk_dynspg_rl .OR. lk_vvl ) THEN         ! rigid lid or non-linear free surface
237         zwx(:,:, 1 ) = 0.e0                           ! Surface value : zero flux
238         zwx(:,:,jpk) = 0.e0                           ! Bottom  value : flux set to zero
239      ELSE                                        ! linear free surface
240         zwx(:,:, 1 ) = pwn(:,:,1) * ptn(:,:,1)        ! Surface :  : advection through z=0
241         zwx(:,:,jpk) = 0.e0                           ! Bottom  : flux set to zero
242      ENDIF
243
244      DO jk = 2, jpk                                   ! Vertical advective fluxes (at w-point)
245         DO jj = 2, jpjm1
246            DO ji = fs_2, fs_jpim1   ! vector opt.
247               ! upstream indicator
248               zcofk = MAX( zind(ji,jj,jk-1), zind(ji,jj,jk) )
249               ! velocity * 1/2
250               zhw = 0.5 * pwn(ji,jj,jk)
251               ! upstream scheme
252               zfp_w = zhw + ABS( zhw )
253               zfm_w = zhw - ABS( zhw )
254               zupst = zfp_w * ptb(ji,jj,jk) + zfm_w * ptb(ji,jj,jk-1)
255               ! centered scheme
256               zcent = zhw * ( ptn(ji,jj,jk) + ptn(ji,jj,jk-1) )
257               ! mixed centered / upstream scheme
258               zwx(ji,jj,jk) = zcofk * zupst + (1.-zcofk) * zcent
259            END DO
260         END DO
261      END DO
262
263      DO jk = 1, jpkm1                                 ! Tracer flux divergence at t-point added to the general trend
264         DO jj = 2, jpjm1
265            DO ji = fs_2, fs_jpim1   ! vector opt.
266               ze3tr = 1. / fse3t(ji,jj,jk)
267               ! vertical advective trends
268               zta = - ze3tr * ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) )
269               ! add it to the general tracer trends
270               pta(ji,jj,jk) =  pta(ji,jj,jk) + zta
271            END DO
272         END DO
273      END DO
274
275      ! 3. Save the vertical advective trends for diagnostic
276      ! ----------------------------------------------------
277      IF( l_trdtra )   CALL trd_tra_adv( kt, ktra, jpt_trd_zad, cdtype, zwx, pwn, ptn )
278
279      IF(ln_ctl)   CALL prt_ctl( tab3d_1=pta, clinfo1=' cen2 - zad : ', mask1=tmask, clinfo3=cdtype )
280      !
281   END SUBROUTINE tra_adv_cen2
282
283   !!======================================================================
284END MODULE traadv_cen2
Note: See TracBrowser for help on using the repository browser.