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

Last change on this file since 473 was 457, checked in by opalod, 18 years ago

nemo_v1_update_049:RB: reorganization of tracers part, remove traadv_cen2_atsk.h90 traldf_iso_zps.F90 trazdf_iso.F90 trazdf_iso_vopt.F90, change atsk routines to jki
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 16.3 KB
Line 
1MODULE traadv_cen2
2   !!==============================================================================
3   !!                       ***  MODULE  traadv_cen2  ***
4   !! Ocean active tracers:  horizontal & vertical advective trend
5   !!==============================================================================
6   !! History :
7   !!   8.2  !  01-08  (G. Madec, E. Durand)  trahad+trazad = traadv
8   !!   8.5  !  02-06  (G. Madec)  F90: Free form and module
9   !!   9.0  !  04-08  (C. Talandier) New trends organization
10   !!    "   !  05-11  (V. Garnier) Surface pressure gradient organization
11   !!    "   !  06-04  (R. Benshila, G. Madec) Step reorganization
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   !! * Modules used
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   !! * Accessibility
36   PUBLIC tra_adv_cen2    ! routine called by step.F90
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   !! This software is governed by the CeCILL licence see 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      !!         C A U T I O N : the trend saved is the centered trend only.
89      !!      It doesn't take into account the upstream part of the scheme.
90      !!         NB:in z-coordinate - full step (ln_zco=T) e3u=e3v=e3t, so
91      !!      they vanish from the expression of the flux and divergence.
92      !!
93      !!         Part II : vertical advection
94      !!      For temperature (idem for salinity) the advective trend is com-
95      !!      puted as follows :
96      !!            zta = 1/e3t dk+1[ zwz ]
97      !!      where the vertical advective flux, zwz, is given by :
98      !!            zwz = zcofk * zupst + (1-zcofk) * zcent
99      !!      with
100      !!        zupsv = upstream flux = wn * (tb(k) or tb(k-1) ) [wn>0 or <0]
101      !!        zcenu = centered flux = wn * mk(tn)
102      !!         The surface boundary condition is :
103      !!      rigid-lid (lk_dynspg_frd = T) : zero advective flux
104      !!      free-surf (lk_dynspg_fsc = T) : wn(:,:,1) * tn(:,:,1)
105      !!         Add this trend now to the general trend of tracer (ta,sa):
106      !!            (ta,sa) = (ta,sa) + ( zta , zsa )
107      !!         Trend diagnostic ('key_trdtra' defined): the trend is
108      !!      saved for diagnostics. The trends saved is expressed as :
109      !!             save trend =  w.gradz(T) = zta - tn divn.
110      !!
111      !! ** Action :  - update (ta,sa) with the now advective tracer trends
112      !!              - save trends in (ttrdh,ttrd,strdhi,strd) ('key_trdtra')
113      !!
114      !!----------------------------------------------------------------------
115      !! * Modules used
116      USE oce                , zwx => ua,  &  ! use ua as workspace
117         &                     zwy => va      ! use va as workspace
118
119      !! * Arguments
120      INTEGER , INTENT( in ) ::   kt          ! ocean time-step index
121      REAL(wp), INTENT( in ), DIMENSION(jpi,jpj,jpk) ::   &
122         pun, pvn, pwn                        ! now ocean velocity fields
123 
124      !! * Local save
125      REAL(wp), DIMENSION(jpi,jpj), SAVE ::   &
126         zbtr2
127 
128      !! * Local declarations
129      INTEGER  ::   ji, jj, jk                 ! dummy loop indices
130      REAL(wp) ::                           &
131         zbtr, zta, zsa, zfui, zfvj,        &  ! temporary scalars
132         zhw, ze3tr, zcofi, zcofj,          &  !    "         "
133         zupsut, zupsvt, zupsus, zupsvs,    &  !    "         "
134         zfp_ui, zfp_vj, zfm_ui, zfm_vj,    &  !    "         "
135         zcofk, zupst, zupss, zcent,        &  !    "         "
136         zcens, zfp_w, zfm_w,               &  !    "         "
137         zcenut, zcenvt, zcenus, zcenvs        !    "         "
138      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &
139         zwz, zww, zind,                    &  ! temporary workspace arrays
140         ztdta, ztdsa                          !    "         "
141      !!----------------------------------------------------------------------
142
143      IF( kt == nit000 ) THEN
144         IF(lwp) WRITE(numout,*)
145         IF(lwp) WRITE(numout,*) 'tra_adv_cen2 : 2nd order centered advection scheme'
146         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~   Vector optimization case'
147         IF(lwp) WRITE(numout,*)
148   
149         zbtr2(:,:) = 1. / ( e1t(:,:) * e2t(:,:) )
150      ENDIF
151
152      ! Save ta and sa trends
153      IF( l_trdtra )   THEN
154         ztdta(:,:,:) = ta(:,:,:) 
155         ztdsa(:,:,:) = sa(:,:,:) 
156         l_adv = 'ce2'
157      ENDIF
158
159      ! Upstream / centered scheme indicator
160      ! ------------------------------------
161 
162      DO jk = 1, jpk
163         DO jj = 1, jpj
164            DO ji = 1, jpi
165               zind(ji,jj,jk) =  MAX ( upsrnfh(ji,jj) * upsrnfz(jk),     &  ! changing advection scheme near runoff
166                  &                    upsadv(ji,jj)                     &  ! in the vicinity of some straits
167#if defined key_ice_lim
168                  &                  , tmask(ji,jj,jk)                   &  ! half upstream tracer fluxes
169                  &                  * MAX( 0., SIGN( 1., fzptn(ji,jj)   &  ! if tn < ("freezing"+0.1 )
170                  &                                +0.1-tn(ji,jj,jk) ) ) &
171#endif
172                  &                  )
173            END DO
174         END DO
175      END DO
176
177
178      ! I. Horizontal advective fluxes
179      ! ------------------------------
180
181      ! 1. Second order centered tracer flux at u and v-points
182      ! -------------------------------------------------------
183
184      !                                                ! ===============
185      DO jk = 1, jpkm1                                 ! Horizontal slab
186         !                                             ! ===============
187         DO jj = 1, jpjm1
188            DO ji = 1, fs_jpim1   ! vector opt.
189               ! upstream indicator
190               zcofi = MAX( zind(ji+1,jj,jk), zind(ji,jj,jk) )
191               zcofj = MAX( zind(ji,jj+1,jk), zind(ji,jj,jk) )
192               ! volume fluxes * 1/2
193#if defined key_zco
194               zfui = 0.5 * e2u(ji,jj) * pun(ji,jj,jk)
195               zfvj = 0.5 * e1v(ji,jj) * pvn(ji,jj,jk)
196#else
197               zfui = 0.5 * e2u(ji,jj) * fse3u(ji,jj,jk) * pun(ji,jj,jk)
198               zfvj = 0.5 * e1v(ji,jj) * fse3v(ji,jj,jk) * pvn(ji,jj,jk)
199#endif
200               ! upstream scheme
201               zfp_ui = zfui + ABS( zfui )
202               zfp_vj = zfvj + ABS( zfvj )
203               zfm_ui = zfui - ABS( zfui )
204               zfm_vj = zfvj - ABS( zfvj )
205               zupsut = zfp_ui * tb(ji,jj,jk) + zfm_ui * tb(ji+1,jj  ,jk)
206               zupsvt = zfp_vj * tb(ji,jj,jk) + zfm_vj * tb(ji  ,jj+1,jk)
207               zupsus = zfp_ui * sb(ji,jj,jk) + zfm_ui * sb(ji+1,jj  ,jk)
208               zupsvs = zfp_vj * sb(ji,jj,jk) + zfm_vj * sb(ji  ,jj+1,jk)
209               ! centered scheme
210               zcenut = zfui * ( tn(ji,jj,jk) + tn(ji+1,jj  ,jk) )
211               zcenvt = zfvj * ( tn(ji,jj,jk) + tn(ji  ,jj+1,jk) )
212               zcenus = zfui * ( sn(ji,jj,jk) + sn(ji+1,jj  ,jk) )
213               zcenvs = zfvj * ( sn(ji,jj,jk) + sn(ji  ,jj+1,jk) )
214               ! mixed centered / upstream scheme
215               zwx(ji,jj,jk) = zcofi * zupsut + (1.-zcofi) * zcenut
216               zwy(ji,jj,jk) = zcofj * zupsvt + (1.-zcofj) * zcenvt
217               zww(ji,jj,jk) = zcofi * zupsus + (1.-zcofi) * zcenus
218               zwz(ji,jj,jk) = zcofj * zupsvs + (1.-zcofj) * zcenvs
219            END DO
220         END DO
221
222         ! 2. Tracer flux divergence at t-point added to the general trend
223         ! ---------------------------------------------------------------
224
225         DO jj = 2, jpjm1
226            DO ji = fs_2, fs_jpim1   ! vector opt.
227#if defined key_zco
228               zbtr = zbtr2(ji,jj)
229#else
230               zbtr = zbtr2(ji,jj) / fse3t(ji,jj,jk)
231#endif
232               ! horizontal advective trends
233               zta = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk)   &
234                  &            + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk)  )
235               zsa = - zbtr * (  zww(ji,jj,jk) - zww(ji-1,jj  ,jk)   &
236                  &            + zwz(ji,jj,jk) - zwz(ji  ,jj-1,jk)  )
237               ! add it to the general tracer trends
238               ta(ji,jj,jk) = ta(ji,jj,jk) + zta
239               sa(ji,jj,jk) = sa(ji,jj,jk) + zsa
240            END DO
241         END DO
242         !                                             ! ===============
243      END DO                                           !   End of slab
244      !                                                ! ===============
245
246      ! 3. Save the horizontal advective trends for diagnostic
247      ! ------------------------------------------------------
248      IF( l_trdtra )   THEN
249         ! Recompute the hoizontal advection zta & zsa trends computed
250         ! at the step 2. above in making the difference between the new
251         ! trends and the previous one ta()/sa - ztdta()/ztdsa() and add
252         ! the term tn()/sn()*hdivn() to recover the Uh gradh(T/S) trends
253         ztdta(:,:,:) = ta(:,:,:) - ztdta(:,:,:) + tn(:,:,:) * hdivn(:,:,:) 
254         ztdsa(:,:,:) = sa(:,:,:) - ztdsa(:,:,:) + sn(:,:,:) * hdivn(:,:,:)
255
256         CALL trd_mod(ztdta, ztdsa, jpttdlad, 'TRA', kt)
257
258         ! Save the new ta and sa trends
259         ztdta(:,:,:) = ta(:,:,:) 
260         ztdsa(:,:,:) = sa(:,:,:) 
261
262      ENDIF
263
264      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ta, clinfo1=' cen2 had  - Ta: ', mask1=tmask, &
265         &                       tab3d_2=sa, clinfo2=            ' Sa: ', mask2=tmask, clinfo3='tra' )
266
267      ! "zonal" mean advective heat and salt transport
268      IF( ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN
269         IF( lk_zco ) THEN
270            DO jk = 1, jpkm1
271               DO jj = 2, jpjm1
272                  DO ji = fs_2, fs_jpim1   ! vector opt.
273                    zwy(ji,jj,jk) = zwy(ji,jj,jk) * fse3v(ji,jj,jk)
274                    zwz(ji,jj,jk) = zwz(ji,jj,jk) * fse3v(ji,jj,jk)
275                  END DO
276               END DO
277            END DO
278         ENDIF
279         pht_adv(:) = ptr_vj( zwy(:,:,:) )
280         pst_adv(:) = ptr_vj( zwz(:,:,:) )
281      ENDIF
282
283      ! II. Vertical advection
284      ! ----------------------
285
286      ! Bottom value : flux set to zero
287      zwx(:,:,jpk) = 0.e0     ;    zwy(:,:,jpk) = 0.e0
288
289      ! Surface value
290      IF( lk_dynspg_rl ) THEN
291         ! rigid lid : flux set to zero
292         zwx(:,:, 1 ) = 0.e0    ;    zwy(:,:, 1 ) = 0.e0
293      ELSE
294         ! free surface
295         zwx(:,:, 1 ) = pwn(:,:,1) * tn(:,:,1)
296         zwy(:,:, 1 ) = pwn(:,:,1) * sn(:,:,1)
297      ENDIF
298
299      ! 1. Vertical advective fluxes
300      ! ----------------------------
301
302      ! Second order centered tracer flux at w-point
303
304      DO jk = 2, jpk
305         DO jj = 2, jpjm1
306            DO ji = fs_2, fs_jpim1   ! vector opt.
307               ! upstream indicator
308               zcofk = MAX( zind(ji,jj,jk-1), zind(ji,jj,jk) )
309               ! velocity * 1/2
310               zhw = 0.5 * pwn(ji,jj,jk)
311               ! upstream scheme
312               zfp_w = zhw + ABS( zhw )
313               zfm_w = zhw - ABS( zhw )
314               zupst = zfp_w * tb(ji,jj,jk) + zfm_w * tb(ji,jj,jk-1)
315               zupss = zfp_w * sb(ji,jj,jk) + zfm_w * sb(ji,jj,jk-1)
316               ! centered scheme
317               zcent = zhw * ( tn(ji,jj,jk) + tn(ji,jj,jk-1) )
318               zcens = zhw * ( sn(ji,jj,jk) + sn(ji,jj,jk-1) )
319               ! mixed centered / upstream scheme
320               zwx(ji,jj,jk) = zcofk * zupst + (1.-zcofk) * zcent
321               zwy(ji,jj,jk) = zcofk * zupss + (1.-zcofk) * zcens
322            END DO
323         END DO
324      END DO
325
326
327      ! 2. Tracer flux divergence at t-point added to the general trend
328      ! -------------------------
329 
330      DO jk = 1, jpkm1
331         DO jj = 2, jpjm1
332            DO ji = fs_2, fs_jpim1   ! vector opt.
333               ze3tr = 1. / fse3t(ji,jj,jk)
334               ! vertical advective trends
335               zta = - ze3tr * ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) )
336               zsa = - ze3tr * ( zwy(ji,jj,jk) - zwy(ji,jj,jk+1) )
337               ! add it to the general tracer trends
338               ta(ji,jj,jk) =  ta(ji,jj,jk) + zta
339               sa(ji,jj,jk) =  sa(ji,jj,jk) + zsa
340            END DO
341         END DO
342      END DO
343
344      ! 3. Save the vertical advective trends for diagnostic
345      ! ----------------------------------------------------
346      IF( l_trdtra )   THEN
347         ! Recompute the vertical advection zta & zsa trends computed
348         ! at the step 2. above in making the difference between the new
349         ! trends and the previous one: ta()/sa - ztdta()/ztdsa() and substract
350         ! the term tn()/sn()*hdivn() to recover the W gradz(T/S) trends
351         ztdta(:,:,:) = ta(:,:,:) - ztdta(:,:,:) - tn(:,:,:) * hdivn(:,:,:) 
352         ztdsa(:,:,:) = sa(:,:,:) - ztdsa(:,:,:) - sn(:,:,:) * hdivn(:,:,:)
353
354         CALL trd_mod(ztdta, ztdsa, jpttdzad, 'TRA', kt)
355      ENDIF
356
357      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ta, clinfo1=' cen2 zad  - Ta: ', mask1=tmask, &
358         &                       tab3d_2=sa, clinfo2=            ' Sa: ', mask2=tmask, clinfo3='tra' )
359
360
361   END SUBROUTINE tra_adv_cen2
362
363   !!======================================================================
364END MODULE traadv_cen2
Note: See TracBrowser for help on using the repository browser.