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

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

Initial revision
  • 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 , 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 ('key_s_coord' defined) or
29      !!            * z-coordinate with partial steps ('key_partial_steps'),
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 (default key), e3t=e3u=e3v:
34      !!               zcenu = e2u  un  mi(tn)
35      !!               zcenv = e1v  vn  mj(tn)
36      !!       * upstream flux:
37      !!            * s-coordinate ('key_s_coord' defined) or
38      !!            * z-coordinate with partial steps ('key_partial_steps')
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 ('key_s_coord' defined)
51      !!              or z-coordinate with partial steps ('key_partial_steps')
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 (key_dynspg_frd = T) : zero advective flux
80      !!      free-surf (key_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      !!----------------------------------------------------------------------
94      !! * Modules used
95      USE oce              , zwx => ua,  &  ! use ua as workspace
96         &                   zwy => va      ! use va as workspace
97#if defined key_trabbl_adv
98      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  &  ! temporary arrays
99         &         zun, zvn, zwn
100#else
101      USE oce                , zun => un,  &  ! When no bbl, zun == un
102         &                     zvn => vn,  &  ! When no bbl, zvn == vn
103         &                     zwn => wn      ! When no bbl, zwn == wn
104#endif
105      !! * Arguments
106      INTEGER, INTENT( in ) ::   kt         ! ocean time-step index
107
108      !! * Local save
109      REAL(wp), DIMENSION(jpi,jpj), SAVE ::   &
110         zbtr2
111
112      !! * Local declarations
113      INTEGER ::   ji, jj, jk               ! dummy loop indices
114      REAL(wp) ::   zbtr, zta, zsa, zfui, zfvj
115      REAL(wp) ::   zhw, ze3tr
116      REAL(wp) ::   zcofi, zcofj, zupsut, zupsvt, zupsus, zupsvs,   &
117                               zfp_ui, zfp_vj, zfm_ui, zfm_vj
118      REAL(wp) ::   zcofk, zupst, zupss, zcent, zcens, zfp_w, zfm_w
119      REAL(wp) ::   zcenut, zcenvt, zcenus, zcenvs
120      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &
121         zwz, zww, zind                     ! workspace
122#if defined key_trdtra || defined key_trdmld
123      REAL(wp) ::   ztai, ztaj, zsai, zsaj
124      REAL(wp) ::   zfui1, zfvj1
125#endif
126      !!----------------------------------------------------------------------
127
128      IF( kt == nit000 ) THEN
129         IF(lwp) WRITE(numout,*)
130         IF(lwp) WRITE(numout,*) 'tra_adv_cen2 : 2nd order centered advection scheme'
131         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~   Auto-tasking case'
132         IF(lwp) WRITE(numout,*)
133         
134         zbtr2(:,:) = 1. / ( e1t(:,:) * e2t(:,:) )
135      ENDIF
136
137      !                                                ! ===============
138      DO jk = 1, jpkm1                                 ! Horizontal slab
139         !                                             ! ===============
140#if defined key_trabbl_adv
141
142         ! Advective bottom boundary layer
143         ! -------------------------------
144         zun(:,:,jk) = un (:,:,jk) - u_bbl(:,:,jk)
145         zvn(:,:,jk) = vn (:,:,jk) - v_bbl(:,:,jk)
146         zwn(:,:,jk) = wn (:,:,jk) + w_bbl(:,:,jk)
147#endif
148
149
150         ! 0. Upstream / centered scheme indicator
151         ! ---------------------------------------
152         DO jj = 1, jpj
153            DO ji = 1, jpi
154               zind(ji,jj,jk) =  MAX (              &
155                  upsrnfh(ji,jj) * upsrnfz(jk),     &  ! changing advection scheme near runoff
156                  upsadv(ji,jj)                     &  ! in the vicinity of some straits
157#if defined key_ice_lim
158                  , tmask(ji,jj,jk)                 &  ! half upstream tracer fluxes if tn < ("freezing"+0.1 )
159                  * MAX( 0., SIGN( 1., fzptn(ji,jj)+0.1-tn(ji,jj,jk) ) )   &
160#endif
161               )
162            END DO
163         END DO
164
165
166         ! I. Horizontal advective fluxes
167         ! ------------------------------
168         ! Second order centered tracer flux at u and v-points
169         DO jj = 1, jpjm1
170            DO ji = 1, fs_jpim1   ! vector opt.
171               ! upstream indicator
172               zcofi = MAX( zind(ji+1,jj,jk), zind(ji,jj,jk) )
173               zcofj = MAX( zind(ji,jj+1,jk), zind(ji,jj,jk) )
174               ! volume fluxes * 1/2
175#if defined key_s_coord || defined key_partial_steps
176               zfui = 0.5 * e2u(ji,jj) * fse3u(ji,jj,jk) * zun(ji,jj,jk)
177               zfvj = 0.5 * e1v(ji,jj) * fse3v(ji,jj,jk) * zvn(ji,jj,jk)
178#else
179               zfui = 0.5 * e2u(ji,jj) * zun(ji,jj,jk)
180               zfvj = 0.5 * e1v(ji,jj) * zvn(ji,jj,jk)
181#endif
182               ! upstream scheme
183               zfp_ui = zfui + ABS( zfui )
184               zfp_vj = zfvj + ABS( zfvj )
185               zfm_ui = zfui - ABS( zfui )
186               zfm_vj = zfvj - ABS( zfvj )
187               zupsut = zfp_ui * tb(ji,jj,jk) + zfm_ui * tb(ji+1,jj  ,jk)
188               zupsvt = zfp_vj * tb(ji,jj,jk) + zfm_vj * tb(ji  ,jj+1,jk)
189               zupsus = zfp_ui * sb(ji,jj,jk) + zfm_ui * sb(ji+1,jj  ,jk)
190               zupsvs = zfp_vj * sb(ji,jj,jk) + zfm_vj * sb(ji  ,jj+1,jk)
191               ! centered scheme
192               zcenut = zfui * ( tn(ji,jj,jk) + tn(ji+1,jj  ,jk) )
193               zcenvt = zfvj * ( tn(ji,jj,jk) + tn(ji  ,jj+1,jk) )
194               zcenus = zfui * ( sn(ji,jj,jk) + sn(ji+1,jj  ,jk) )
195               zcenvs = zfvj * ( sn(ji,jj,jk) + sn(ji  ,jj+1,jk) )
196               ! mixed centered / upstream scheme
197               zwx(ji,jj,jk) = zcofi * zupsut + (1.-zcofi) * zcenut
198               zwy(ji,jj,jk) = zcofj * zupsvt + (1.-zcofj) * zcenvt
199               zww(ji,jj,jk) = zcofi * zupsus + (1.-zcofi) * zcenus
200               zwz(ji,jj,jk) = zcofj * zupsvs + (1.-zcofj) * zcenvs
201            END DO
202         END DO
203
204         ! 2. Tracer flux divergence at t-point added to the general trend
205         ! -------------------------
206         DO jj = 2, jpjm1
207            DO ji = fs_2, fs_jpim1   ! vector opt.
208#if defined key_s_coord || defined key_partial_steps
209               zbtr = zbtr2(ji,jj) / fse3t(ji,jj,jk)
210#else
211               zbtr = zbtr2(ji,jj)
212#endif
213               ! horizontal advective trends
214               zta = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj,jk) + zwy(ji,jj,jk) - zwy(ji,jj-1,jk)  )
215               zsa = - zbtr * (  zww(ji,jj,jk) - zww(ji-1,jj,jk) + zwz(ji,jj,jk) - zwz(ji,jj-1,jk)  )
216               ! add it to the general tracer trends
217               ta(ji,jj,jk) = ta(ji,jj,jk) + zta
218               sa(ji,jj,jk) = sa(ji,jj,jk) + zsa
219#if defined key_trdtra || defined key_trdmld
220               ! save the horizontal advective trend of tracer
221               ttrd(ji,jj,jk,1) = zta + tn(ji,jj,jk) * hdivn(ji,jj,jk)
222               strd(ji,jj,jk,1) = zsa + sn(ji,jj,jk) * hdivn(ji,jj,jk)
223               ! recompute the trends in i- and j-direction as Uh gradh(T)
224#   if defined key_s_coord || defined key_partial_steps
225               zfui = 0.5 * e2u(ji  ,jj) * fse3u(ji,  jj,jk) * zun(ji,  jj,jk)
226               zfui1= 0.5 * e2u(ji-1,jj) * fse3u(ji-1,jj,jk) * zun(ji-1,jj,jk)
227               zfvj = 0.5 * e1v(ji,jj  ) * fse3v(ji,jj  ,jk) * zvn(ji,jj  ,jk)
228               zfvj1= 0.5 * e1v(ji,jj-1) * fse3v(ji,jj-1,jk) * zvn(ji,jj-1,jk)
229#   else
230               zfui = 0.5 * e2u(ji  ,jj) * zun(ji,  jj,jk)
231               zfui1= 0.5 * e2u(ji-1,jj) * zun(ji-1,jj,jk)
232               zfvj = 0.5 * e1v(ji,jj  ) * zvn(ji,jj  ,jk)
233               zfvj1= 0.5 * e1v(ji,jj-1) * zvn(ji,jj-1,jk)
234#   endif
235               ztai = -zbtr * (  zfui * ( tn(ji+1,jj,jk) - tn(ji,jj,jk) ) + zfui1* ( tn(ji,jj,jk) - tn(ji-1,jj,jk) )  )
236               zsai = -zbtr * (  zfui * ( sn(ji+1,jj,jk) - sn(ji,jj,jk) ) + zfui1* ( sn(ji,jj,jk) - sn(ji-1,jj,jk) )  )
237               ztaj = -zbtr * (  zfvj * ( tn(ji,jj+1,jk) - tn(ji,jj,jk) ) + zfvj1* ( tn(ji,jj,jk) - tn(ji,jj-1,jk) )  )
238               zsaj = -zbtr * (  zfvj * ( sn(ji,jj+1,jk) - sn(ji,jj,jk) ) + zfvj1* ( sn(ji,jj,jk) - sn(ji,jj-1,jk) )  )
239               ! save i- and j- advective trends computed as Uh gradh(T)
240               ttrdh(ji,jj,jk,1) = ztai
241               ttrdh(ji,jj,jk,2) = ztaj
242               strdh(ji,jj,jk,1) = zsai
243               strdh(ji,jj,jk,2) = zsaj
244#endif
245            END DO
246         END DO
247         !                                             ! ===============
248      END DO                                           !   End of slab
249      !                                                ! ===============
250
251      IF( l_ctl .AND. lwp ) THEN
252         zta = SUM( ta(2:jpim1,2:jpjm1,1:jpkm1) * tmask(2:jpim1,2:jpjm1,1:jpkm1) )
253         zsa = SUM( sa(2:jpim1,2:jpjm1,1:jpkm1) * tmask(2:jpim1,2:jpjm1,1:jpkm1) )
254         WRITE(numout,*) ' had  - Ta: ', zta-t_ctl, ' Sa: ', zsa-s_ctl
255         t_ctl = zta   ;   s_ctl = zsa
256      ENDIF
257
258#if defined key_diaptr
259      ! "zonal" mean advective heat and salt transport
260      IF( MOD( kt, nf_ptr ) == 0 ) THEN
261# if defined key_s_coord || defined key_partial_steps
262         pht_adv(:,:) = prt_vj( zwy(:,:,:) )
263         pst_adv(:,:) = prt_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(:,:) = prt_vj( zwy(:,:,:) )
274         pst_adv(:,:) = prt_vj( zwz(:,:,:) )
275# endif
276      ENDIF
277#endif
278
279      ! II. Vertical advection
280      ! ----------------------
281      !                                                ! ===============
282      DO jj = 2, jpjm1                                 !  Vertical slab
283         !                                             ! ===============
284         ! Bottom value : flux and indicator set to zero
285         zwz (:,jj,jpk) = 0.e0     ;    zww(:,jj,jpk) = 0.e0
286         zind(:,jj,jpk) = 0.e0
287
288         ! Surface value
289#if defined key_dynspg_fsc
290         ! free surface-constant volume
291         zwz(:,jj, 1 ) = zwn(:,jj,1) * tn(:,jj,1)
292         zww(:,jj, 1 ) = zwn(:,jj,1) * sn(:,jj,1)
293#endif
294#if defined key_dynspg_rl
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#if defined key_trdtra || defined key_trdmld
334               ! save the vertical advective trends computed as w gradz(T)
335               ttrd(ji,jj,jk,2) = zta - tn(ji,jj,jk) * hdivn(ji,jj,jk)
336               strd(ji,jj,jk,2) = zsa - sn(ji,jj,jk) * hdivn(ji,jj,jk)
337#endif
338            END DO
339         END DO
340         !                                             ! ===============
341      END DO                                           !   End of slab
342      !                                                ! ===============
343
344      IF( l_ctl .AND. lwp ) THEN
345         zta = SUM( ta(2:jpim1,2:jpjm1,1:jpkm1) * tmask(2:jpim1,2:jpjm1,1:jpkm1) )
346         zsa = SUM( sa(2:jpim1,2:jpjm1,1:jpkm1) * tmask(2:jpim1,2:jpjm1,1:jpkm1) )
347         WRITE(numout,*) ' zad  - Ta: ', zta-t_ctl, ' Sa: ', zsa-s_ctl, ' centered2 autotsk'
348         t_ctl = zta   ;   s_ctl = zsa
349      ENDIF
350
351   END SUBROUTINE tra_adv_cen2
Note: See TracBrowser for help on using the repository browser.