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

source: trunk/NEMO/OPA_SRC/TRA/traadv_cen2_jki.F90 @ 473

Last change on this file since 473 was 458, checked in by opalod, 18 years ago

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