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

Last change on this file since 389 was 359, checked in by opalod, 18 years ago

nemo_v1_update_033 : RB + CT : Add new surface pressure gradient algorithms

  • 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 , 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_rl = T) : zero advective flux
82      !!      free-surf                    : 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      !!    "   !  05-11  (V. Garnier) Surface pressure gradient organization
97      !!----------------------------------------------------------------------
98      !! * Modules used
99      USE oce              , zwx => ua,    &  ! use ua as workspace
100         &                   zwy => va        ! use va as workspace
101#if defined key_trabbl_adv
102      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  &  ! temporary arrays
103         &         zun, zvn, zwn
104#else
105      USE oce                , zun => un,  &  ! When no bbl, zun == un
106         &                     zvn => vn,  &  ! When no bbl, zvn == vn
107         &                     zwn => wn      ! When no bbl, zwn == wn
108#endif
109      !! * Arguments
110      INTEGER, INTENT( in ) ::   kt            ! ocean time-step index
111
112      !! * Local save
113      REAL(wp), DIMENSION(jpi,jpj), SAVE ::   &
114         zbtr2
115
116      !! * Local declarations
117      INTEGER  ::   ji, jj, jk                 ! dummy loop indices
118      REAL(wp) ::                           &
119         zbtr, zta, zsa, zfui, zfvj,        &  ! temporary scalars
120         zhw, ze3tr, zcofi, zcofj,          &  !    "         "
121         zupsut, zupsvt, zupsus, zupsvs,    &  !    "         "
122         zfp_ui, zfp_vj, zfm_ui, zfm_vj,    &  !    "         "
123         zcofk, zupst, zupss, zcent,        &  !    "         "
124         zcens, zfp_w, zfm_w,               &  !    "         "
125         zcenut, zcenvt, zcenus, zcenvs,    &  !    "         "
126         zfui1, zfvj1                          !    "         "
127      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &
128         zwz, zww, zind,                    &  ! temporary workspace arrays
129         ztdta, ztdsa                          !    "         "
130      !!----------------------------------------------------------------------
131
132      IF( kt == nit000 ) THEN
133         IF(lwp) WRITE(numout,*)
134         IF(lwp) WRITE(numout,*) 'tra_adv_cen2 : 2nd order centered advection scheme'
135         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~   Auto-tasking case'
136         IF(lwp) WRITE(numout,*)
137         
138         zbtr2(:,:) = 1. / ( e1t(:,:) * e2t(:,:) )
139      ENDIF
140
141      ! Save ta and sa trends
142      IF( l_trdtra )   THEN
143         ztdta(:,:,:) = ta(:,:,:)
144         ztdsa(:,:,:) = sa(:,:,:)
145         l_adv = 'ce2'
146      ENDIF
147
148      !                                                ! ===============
149      DO jk = 1, jpkm1                                 ! Horizontal slab
150         !                                             ! ===============
151#if defined key_trabbl_adv
152         ! Advective bottom boundary layer
153         ! -------------------------------
154         zun(:,:,jk) = un (:,:,jk) - u_bbl(:,:,jk)
155         zvn(:,:,jk) = vn (:,:,jk) - v_bbl(:,:,jk)
156         zwn(:,:,jk) = wn (:,:,jk) + w_bbl(:,:,jk)
157#endif
158
159         ! 0. Upstream / centered scheme indicator
160         ! ---------------------------------------
161         DO jj = 1, jpj
162            DO ji = 1, jpi
163               zind(ji,jj,jk) =  MAX (              &
164                  upsrnfh(ji,jj) * upsrnfz(jk),     &  ! changing advection scheme near runoff
165                  upsadv(ji,jj)                     &  ! in the vicinity of some straits
166#if defined key_ice_lim
167                  , tmask(ji,jj,jk)                 &  ! half upstream tracer fluxes if tn < ("freezing"+0.1 )
168                  * MAX( 0., SIGN( 1., fzptn(ji,jj)+0.1-tn(ji,jj,jk) ) )   &
169#endif
170               )
171            END DO
172         END DO
173
174
175         ! I. Horizontal advective fluxes
176         ! ------------------------------
177         ! Second order centered tracer flux at u and v-points
178         DO jj = 1, jpjm1
179            DO ji = 1, fs_jpim1   ! vector opt.
180               ! upstream indicator
181               zcofi = MAX( zind(ji+1,jj,jk), zind(ji,jj,jk) )
182               zcofj = MAX( zind(ji,jj+1,jk), zind(ji,jj,jk) )
183               ! volume fluxes * 1/2
184#if defined key_s_coord || defined key_partial_steps
185               zfui = 0.5 * e2u(ji,jj) * fse3u(ji,jj,jk) * zun(ji,jj,jk)
186               zfvj = 0.5 * e1v(ji,jj) * fse3v(ji,jj,jk) * zvn(ji,jj,jk)
187#else
188               zfui = 0.5 * e2u(ji,jj) * zun(ji,jj,jk)
189               zfvj = 0.5 * e1v(ji,jj) * zvn(ji,jj,jk)
190#endif
191               ! upstream scheme
192               zfp_ui = zfui + ABS( zfui )
193               zfp_vj = zfvj + ABS( zfvj )
194               zfm_ui = zfui - ABS( zfui )
195               zfm_vj = zfvj - ABS( zfvj )
196               zupsut = zfp_ui * tb(ji,jj,jk) + zfm_ui * tb(ji+1,jj  ,jk)
197               zupsvt = zfp_vj * tb(ji,jj,jk) + zfm_vj * tb(ji  ,jj+1,jk)
198               zupsus = zfp_ui * sb(ji,jj,jk) + zfm_ui * sb(ji+1,jj  ,jk)
199               zupsvs = zfp_vj * sb(ji,jj,jk) + zfm_vj * sb(ji  ,jj+1,jk)
200               ! centered scheme
201               zcenut = zfui * ( tn(ji,jj,jk) + tn(ji+1,jj  ,jk) )
202               zcenvt = zfvj * ( tn(ji,jj,jk) + tn(ji  ,jj+1,jk) )
203               zcenus = zfui * ( sn(ji,jj,jk) + sn(ji+1,jj  ,jk) )
204               zcenvs = zfvj * ( sn(ji,jj,jk) + sn(ji  ,jj+1,jk) )
205               ! mixed centered / upstream scheme
206               zwx(ji,jj,jk) = zcofi * zupsut + (1.-zcofi) * zcenut
207               zwy(ji,jj,jk) = zcofj * zupsvt + (1.-zcofj) * zcenvt
208               zww(ji,jj,jk) = zcofi * zupsus + (1.-zcofi) * zcenus
209               zwz(ji,jj,jk) = zcofj * zupsvs + (1.-zcofj) * zcenvs
210            END DO
211         END DO
212
213         ! 2. Tracer flux divergence at t-point added to the general trend
214         ! ---------------------------------------------------------------
215
216         DO jj = 2, jpjm1
217            DO ji = fs_2, fs_jpim1   ! vector opt.
218#if defined key_s_coord || defined key_partial_steps
219               zbtr = zbtr2(ji,jj) / fse3t(ji,jj,jk)
220#else
221               zbtr = zbtr2(ji,jj)
222#endif
223               ! horizontal advective trends
224               zta = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk)   &
225                  &            + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk)  )
226               zsa = - zbtr * (  zww(ji,jj,jk) - zww(ji-1,jj  ,jk)   &
227                  &            + zwz(ji,jj,jk) - zwz(ji  ,jj-1,jk)  )
228               ! add it to the general tracer trends
229               ta(ji,jj,jk) = ta(ji,jj,jk) + zta
230               sa(ji,jj,jk) = sa(ji,jj,jk) + zsa
231            END DO
232         END DO
233         !                                             ! ===============
234      END DO                                           !   End of slab
235      !                                                ! ===============
236
237      ! 3. Save the horizontal advective trends for diagnostic
238      ! ------------------------------------------------------
239      IF( l_trdtra )   THEN
240         ! Recompute the hoizontal advection zta & zsa trends computed
241         ! at the step 2. above in making the difference between the new
242         ! trends and the previous one ta()/sa - ztdta()/ztdsa() and add
243         ! the term tn()/sn()*hdivn() to recover the Uh gradh(T/S) trends
244         ztdta(:,:,:) = ta(:,:,:) - ztdta(:,:,:) + tn(:,:,:) * hdivn(:,:,:)
245         ztdsa(:,:,:) = sa(:,:,:) - ztdsa(:,:,:) + sn(:,:,:) * hdivn(:,:,:)
246
247         CALL trd_mod(ztdta, ztdsa, jpttdlad, 'TRA', kt)
248
249         ! Save the new ta and sa trends
250         ztdta(:,:,:) = ta(:,:,:)
251         ztdsa(:,:,:) = sa(:,:,:)
252
253      ENDIF
254
255      IF(ln_ctl)   THEN
256          CALL prt_ctl(tab3d_1=ta, clinfo1=' centered2 had  - Ta: ', mask1=tmask, &
257             &         tab3d_2=sa, clinfo2=' Sa: ', mask2=tmask, clinfo3='tra')
258      ENDIF
259
260      ! "zonal" mean advective heat and salt transport
261      IF( ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN
262# if defined key_s_coord || defined key_partial_steps
263         pht_adv(:) = ptr_vj( zwy(:,:,:) )
264         pst_adv(:) = ptr_vj( zwz(:,:,:) )
265# else
266         DO jk = 1, jpkm1
267            DO jj = 2, jpjm1
268               DO ji = fs_2, fs_jpim1   ! vector opt.
269                 zwy(ji,jj,jk) = zwy(ji,jj,jk) * fse3v(ji,jj,jk)
270                 zwz(ji,jj,jk) = zwz(ji,jj,jk) * fse3v(ji,jj,jk)
271               END DO
272            END DO
273         END DO
274         pht_adv(:) = ptr_vj( zwy(:,:,:) )
275         pst_adv(:) = ptr_vj( zwz(:,:,:) )
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( lk_dynspg_rl ) THEN
290            ! rigid lid : flux set to zero
291            zwz(:,jj, 1 ) = 0.e0    ;    zww(:,jj, 1 ) = 0.e0
292         ELSE
293            ! free surface
294            zwz(:,jj, 1 ) = zwn(:,jj,1) * tn(:,jj,1)
295            zww(:,jj, 1 ) = zwn(:,jj,1) * sn(:,jj,1)
296         ENDIF
297         
298         ! 1. Vertical advective fluxes
299         ! ----------------------------
300         ! Second order centered tracer flux at w-point
301         DO jk = 2, jpk
302            DO ji = 2, jpim1
303               ! upstream indicator
304               zcofk = MAX( zind(ji,jj,jk-1), zind(ji,jj,jk) )
305               ! velocity * 1/2
306               zhw = 0.5 * zwn(ji,jj,jk)
307               ! upstream scheme
308               zfp_w = zhw + ABS( zhw )
309               zfm_w = zhw - ABS( zhw )
310               zupst = zfp_w * tb(ji,jj,jk) + zfm_w * tb(ji,jj,jk-1)
311               zupss = zfp_w * sb(ji,jj,jk) + zfm_w * sb(ji,jj,jk-1)
312               ! centered scheme
313               zcent = zhw * ( tn(ji,jj,jk) + tn(ji,jj,jk-1) )
314               zcens = zhw * ( sn(ji,jj,jk) + sn(ji,jj,jk-1) )
315               ! mixed centered / upstream scheme
316               zwz(ji,jj,jk) = zcofk * zupst + (1.-zcofk) * zcent
317               zww(ji,jj,jk) = zcofk * zupss + (1.-zcofk) * zcens
318            END DO
319         END DO
320
321         ! 2. Tracer flux divergence at t-point added to the general trend
322         ! -------------------------
323         DO jk = 1, jpkm1
324            DO ji = 2, jpim1
325               ze3tr = 1. / fse3t(ji,jj,jk)
326               ! vertical advective trends
327               zta = - ze3tr * ( zwz(ji,jj,jk) - zwz(ji,jj,jk+1) )
328               zsa = - ze3tr * ( zww(ji,jj,jk) - zww(ji,jj,jk+1) )
329               ! add it to the general tracer trends
330               ta(ji,jj,jk) =  ta(ji,jj,jk) + zta
331               sa(ji,jj,jk) =  sa(ji,jj,jk) + zsa
332            END DO
333         END DO
334         !                                             ! ===============
335      END DO                                           !   End of slab
336      !                                                ! ===============
337
338      ! 3. Save the vertical advective trends for diagnostic
339      ! ----------------------------------------------------
340      IF( l_trdtra )   THEN
341         ! Recompute the vertical advection zta & zsa trends computed
342         ! at the step 2. above in making the difference between the new
343         ! trends and the previous one: ta()/sa - ztdta()/ztdsa() and substract
344         ! the term tn()/sn()*hdivn() to recover the W gradz(T/S) trends
345         ztdta(:,:,:) = ta(:,:,:) - ztdta(:,:,:) - tn(:,:,:) * hdivn(:,:,:)
346         ztdsa(:,:,:) = sa(:,:,:) - ztdsa(:,:,:) - sn(:,:,:) * hdivn(:,:,:)
347
348         CALL trd_mod(ztdta, ztdsa, jpttdzad, 'TRA', kt)
349      ENDIF
350
351      IF(ln_ctl)   THEN
352          CALL prt_ctl(tab3d_1=ta, clinfo1=' centered2 zad  - Ta: ', mask1=tmask, &
353             &         tab3d_2=sa, clinfo2=' Sa: ', mask2=tmask, clinfo3='tra')
354      ENDIF
355
356   END SUBROUTINE tra_adv_cen2
Note: See TracBrowser for help on using the repository browser.