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