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_atsk.h90 in trunk/NEMO/OPA_SRC/TRA – NEMO

source: trunk/NEMO/OPA_SRC/TRA/traadv_cen2_atsk.h90 @ 216

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

CT : UPDATE151 : New trends organization

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