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

Last change on this file since 106 was 106, checked in by opalod, 20 years ago

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