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

Last change on this file since 247 was 247, checked in by opalod, 19 years ago

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