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 tags/nemo_dev_x1/NEMO/OPA_SRC/TRA – NEMO

source: tags/nemo_dev_x1/NEMO/OPA_SRC/TRA/traadv_cen2.F90 @ 1512

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

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