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.
trcadv_cen2.F90 in trunk/NEMO/TOP_SRC/TRP – NEMO

source: trunk/NEMO/TOP_SRC/TRP/trcadv_cen2.F90 @ 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:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 15.1 KB
Line 
1MODULE trcadv_cen2
2   !!==============================================================================
3   !!                       ***  MODULE  trcadv_cen2  ***
4   !! Ocean passive tracers:  horizontal & vertical advective tracer trend
5   !!==============================================================================
6#if defined key_passivetrc
7   !!----------------------------------------------------------------------
8   !!   trc_adv_cen2 : update the tracer trend with the horizontal
9   !!                  and vertical advection trends using a 2nd order
10   !!                  centered finite difference scheme
11   !!----------------------------------------------------------------------
12   !! * Modules used
13   USE oce_trc             ! ocean dynamics and active tracers variables
14   USE trc                 ! ocean passive tracers variables
15   USE trcbbl              ! advective passive tracers in the BBL
16
17   IMPLICIT NONE
18   PRIVATE
19
20   !! * Accessibility
21   PUBLIC trc_adv_cen2    ! routine called by trcstp.F90
22
23   !! * Substitutions
24#  include "passivetrc_substitute.h90"
25   !!----------------------------------------------------------------------
26   !!   TOP 1.0 , LOCEAN-IPSL (2005)
27   !! $Header$
28   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
29   !!----------------------------------------------------------------------
30
31CONTAINS
32
33   !!----------------------------------------------------------------------
34   !!   Default option :             2nd order centered scheme (k-j-i loop)
35   !!----------------------------------------------------------------------
36
37   SUBROUTINE trc_adv_cen2( kt )
38      !!----------------------------------------------------------------------
39      !!                  ***  ROUTINE trc_adv_cen2  ***
40      !!                 
41      !! ** Purpose :   Compute the now trend due to the advection of tracers
42      !!      and add it to the general trend of passive tracer equations.
43      !!
44      !! ** Method  :   The advection is evaluated by a second order centered
45      !!      scheme using now fields (leap-frog scheme). In specific areas
46      !!      (vicinity of major river mouths, some straits, or where tn is
47      !!      is approaching the freezing point) it is mixed with an upstream
48      !!      scheme for stability reasons.
49      !!        Part 0 : compute the upstream / centered flag
50      !!                 (3D array, zind, defined at T-point (0<zind<1))
51      !!        Part I : horizontal advection
52      !!      * centered flux:
53      !!         * s-coordinate (lk_sco=T) or
54      !!         * z-coordinate with partial steps (lk_zps=T),
55      !!        the vertical scale factors e3. are inside the derivatives:
56      !!               zcenu = e2u*e3u  un  mi(tn)
57      !!               zcenv = e1v*e3v  vn  mj(tn)
58      !!         * z-coordinate (default key), e3t=e3u=e3v:
59      !!               zcenu = e2u  un  mi(tn)
60      !!               zcenv = e1v  vn  mj(tn)
61      !!      * horizontal advective trend (divergence of the fluxes)
62      !!         * s-coordinate (lk_sco=T) or
63      !!         * z-coordinate with partial steps (lk_zps=T)
64      !!               ztra = 1/(e1t*e2t*e3t) { di-1[zwx] + dj-1[zwy] }
65      !!         * z-coordinate (default key), e3t=e3u=e3v:
66      !!               ztra = 1/(e1t*e2t) { di-1[zwx] + dj-1[zwy] }
67      !!      * Add this trend now to the general trend of tracer tra:
68      !!              tra = tra + ztra
69      !!      * trend diagnostic ('key_trc_diatrd'): the trend is saved
70      !!      for diagnostics. The trends saved is expressed as
71      !!      Uh.gradh(T)
72      !!
73      !!         Part II : vertical advection
74      !!      For any tracer  the advective trend is computed as follows :
75      !!            ztra = 1/e3t dk+1[ zwz ]
76      !!      where the vertical advective flux, zwz, is given by :
77      !!            zwz = zcofk * zupst + (1-zcofk) * zcent
78      !!      with
79      !!        zupsv = upstream flux = wn * (trb(k) or trb(k-1) ) [wn>0 or <0]
80      !!        zcenu = centered flux = wn * mk(trn)
81      !!         The surface boundary condition is :
82      !!      rigid-lid (default option) : zero advective flux
83      !!      free-surf ("key_fresurf_cstvol") : wn(:,:,1) * trn(:,:,1)
84      !!         Add this trend now to the general trend of tracer tra :
85      !!            tra = tra + ztra
86      !!         Trend diagnostic ('key_trc_diatrd'): the trend is saved for
87      !!      diagnostics. The trends saved is expressed as :
88      !!             save trend =  w.gradz(T) = ztra - trn divn.
89      !!
90      !! ** Action : - update tra with the now advective tracer trends
91      !!             - save the trends in trtrd ('key_trc_diatrd')
92      !!
93      !! History :
94      !!   8.2  !  01-08  (M-A Filiberti, and M.Levy)  trahad+trazad = traadv
95      !!   8.5  !  02-06  (G. Madec, C. Ethe)  F90: Free form and module
96      !!----------------------------------------------------------------------
97      !! * Modules used
98      USE oce_trc            , zwx => ua,  &  ! use ua as workspace
99         &                     zwy => va      ! use va as workspace
100#if defined key_trcbbl_adv
101      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  &  ! temporary arrays
102         &         zun, zvn, zwn
103#else
104      USE oce_trc            , 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
109 
110      !! * Arguments
111      INTEGER, INTENT( in ) ::   kt           ! ocean time-step index
112 
113      !! * Local save
114      REAL(wp), DIMENSION(jpi,jpj), SAVE ::   &
115         zbtr2
116 
117      !! * Local declarations
118      INTEGER  ::   ji, jj, jk, jn             ! dummy loop indices
119      REAL(wp) ::                           &
120         zbtr, ztra, zfui, zfvj,            &  ! temporary scalars
121         zhw, ze3tr, zcofi, zcofj,          &  !    "         "
122         zupsut, zupsvt,                    &  !    "         "
123         zfp_ui, zfp_vj, zfm_ui, zfm_vj,    &  !    "         "
124         zcofk, zupst, zcent,               &  !    "         "
125         zfp_w, zfm_w,                      &  !    "         "
126         zcenut, zcenvt                        !
127
128      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &
129         zind                              ! temporary workspace arrays
130#if defined key_trc_diatrd
131      REAL(wp) ::                           &
132         ztai, ztaj,                        &  ! temporary scalars
133         zfui1, zfvj1                          !    "         "
134#endif
135      !!----------------------------------------------------------------------
136
137      IF( kt == nittrc000 ) THEN
138         IF(lwp) WRITE(numout,*)
139         IF(lwp) WRITE(numout,*) 'trc_adv_cen2 : 2nd order centered advection scheme'
140         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~   Vector optimization case'
141         IF(lwp) WRITE(numout,*)
142   
143         zbtr2(:,:) = 1. / ( e1t(:,:) * e2t(:,:) )
144      ENDIF
145
146#if defined key_trcbbl_adv
147
148      ! Advective bottom boundary layer
149      ! -------------------------------
150      zun(:,:,:) = un(:,:,:) - u_trc_bbl(:,:,:)
151      zvn(:,:,:) = vn(:,:,:) - v_trc_bbl(:,:,:)
152      zwn(:,:,:) = wn(:,:,:) + w_trc_bbl(:,:,:)
153#endif
154
155
156      DO jn = 1, jptra
157
158      ! Upstream / centered scheme indicator
159      ! ------------------------------------
160 
161         DO jk = 1, jpk
162            DO jj = 1, jpj
163               DO ji = 1, jpi
164                  zind(ji,jj,jk) =  MAX ( 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
168                     &                  * MAX( 0., SIGN( 1., fzptn(ji,jj)   &  ! if tn < ("freezing"+0.1 )
169                     &                                +0.1-tn(ji,jj,jk) ) ) &
170#endif
171                     &                  )
172               END DO
173            END DO
174         END DO
175
176
177         ! I. Horizontal advective fluxes
178         ! ------------------------------
179         
180         ! Second order centered tracer flux at u and v-points
181         
182         !                                                ! ===============
183         DO jk = 1, jpkm1                                 ! Horizontal slab
184            !                                             ! ===============
185            DO jj = 1, jpjm1
186               DO ji = 1, fs_jpim1   ! vector opt.
187               ! upstream indicator
188                  zcofi = MAX( zind(ji+1,jj,jk), zind(ji,jj,jk) )
189                  zcofj = MAX( zind(ji,jj+1,jk), zind(ji,jj,jk) )
190                  ! volume fluxes * 1/2
191#if defined key_s_coord || defined key_partial_steps
192                  zfui = 0.5 * e2u(ji,jj) * fse3u(ji,jj,jk) * zun(ji,jj,jk)
193                  zfvj = 0.5 * e1v(ji,jj) * fse3v(ji,jj,jk) * zvn(ji,jj,jk)
194#else
195                  zfui = 0.5 * e2u(ji,jj) * zun(ji,jj,jk)
196                  zfvj = 0.5 * e1v(ji,jj) * zvn(ji,jj,jk)
197#endif
198               ! upstream scheme
199                  zfp_ui = zfui + ABS( zfui )
200                  zfp_vj = zfvj + ABS( zfvj )
201                  zfm_ui = zfui - ABS( zfui )
202                  zfm_vj = zfvj - ABS( zfvj )
203                  zupsut = zfp_ui * trb(ji,jj,jk,jn) + zfm_ui * trb(ji+1,jj  ,jk,jn)
204                  zupsvt = zfp_vj * trb(ji,jj,jk,jn) + zfm_vj * trb(ji  ,jj+1,jk,jn)
205                  ! centered scheme
206                  zcenut = zfui * ( trn(ji,jj,jk,jn) + trn(ji+1,jj  ,jk,jn) )
207                  zcenvt = zfvj * ( trn(ji,jj,jk,jn) + trn(ji  ,jj+1,jk,jn) )
208                  ! mixed centered / upstream scheme
209                  zwx(ji,jj,jk) = zcofi * zupsut + (1.-zcofi) * zcenut
210                  zwy(ji,jj,jk) = zcofj * zupsvt + (1.-zcofj) * zcenvt               
211               END DO
212            END DO
213
214
215            ! 2. Tracer flux divergence at t-point added to the general trend
216            ! -------------------------
217
218            DO jj = 2, jpjm1
219               DO ji = fs_2, fs_jpim1   ! vector opt.
220#if defined key_s_coord || defined key_partial_steps
221                  zbtr = zbtr2(ji,jj) / fse3t(ji,jj,jk)
222#else
223                  zbtr = zbtr2(ji,jj) 
224#endif
225                  ! horizontal advective trends
226                  ztra = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk)   &
227                     &             + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk)  )
228
229                  ! add it to the general tracer trends
230                  tra(ji,jj,jk,jn) = tra(ji,jj,jk,jn) + ztra
231
232#if defined key_trc_diatrd 
233                  ! recompute the trends in i- and j-direction as Uh gradh(T)
234# if defined key_s_coord || defined key_partial_steps
235                  zfui = 0.5 * e2u(ji  ,jj) * fse3u(ji,  jj,jk) * zun(ji,  jj,jk)
236                  zfui1= 0.5 * e2u(ji-1,jj) * fse3u(ji-1,jj,jk) * zun(ji-1,jj,jk)
237                  zfvj = 0.5 * e1v(ji,jj  ) * fse3v(ji,jj  ,jk) * zvn(ji,jj  ,jk)
238                  zfvj1= 0.5 * e1v(ji,jj-1) * fse3v(ji,jj-1,jk) * zvn(ji,jj-1,jk)
239# else
240                  zfui = 0.5 * e2u(ji  ,jj) * zun(ji,  jj,jk)
241                  zfui1= 0.5 * e2u(ji-1,jj) * zun(ji-1,jj,jk)
242                  zfvj = 0.5 * e1v(ji,jj  ) * zvn(ji,jj  ,jk)
243                  zfvj1= 0.5 * e1v(ji,jj-1) * zvn(ji,jj-1,jk)
244# endif
245                  ztai = - zbtr * ( zfui  * ( trn(ji+1,jj  ,jk,jn) - trn(ji,  jj,jk,jn) )   &
246                     &                + zfui1 * ( trn(ji,  jj,  jk,jn) - trn(ji-1,jj,jk,jn) ) )
247                  ztaj = - zbtr * ( zfvj  * ( trn(ji  ,jj+1,jk,jn) - trn(ji,jj  ,jk,jn) )    &
248                     &                + zfvj1 * ( trn(ji  ,jj  ,jk,jn) - trn(ji,jj-1,jk,jn) ) )
249                  ! save i- and j- advective trends computed as Uh gradh(T)
250                  trtrd(ji,jj,jk,jn,1) = ztai
251                  trtrd(ji,jj,jk,jn,2) = ztaj
252#endif
253               END DO
254            END DO
255            !                                             ! ===============
256         END DO                                           !   End of slab
257         !                                                ! ===============
258
259         IF(l_ctl) THEN         ! print mean trends (used for debugging)
260            ztra = SUM( tra(2:nictl,2:njctl,1:jpkm1,jn) * tmask(2:nictl,2:njctl,1:jpkm1) )
261            WRITE(numout,*) ' trc/had - ',ctrcnm(jn),' : ', ztra-tra_ctl(jn) 
262            tra_ctl(jn) = ztra 
263         ENDIF
264
265         ! II. Vertical advection
266         ! ----------------------
267
268         ! Bottom value : flux set to zero
269         zwx(:,:,jpk) = 0.e0 
270
271         ! Surface value
272         IF( lk_dynspg_fsc ) THEN       ! free surface-constant volume
273            zwx(:,:, 1 ) = zwn(:,:,1) * trn(:,:,1,jn)
274         ELSE                           ! rigid lid : flux set to zero
275            zwx(:,:, 1 ) = 0.e0 
276         ENDIF
277
278         ! 1. Vertical advective fluxes
279         ! ----------------------------
280
281         ! Second order centered tracer flux at w-point
282
283         DO jk = 2, jpk
284            DO jj = 2, jpjm1
285               DO ji = fs_2, fs_jpim1   ! vector opt.
286                  ! upstream indicator
287                  zcofk = MAX( zind(ji,jj,jk-1), zind(ji,jj,jk) )
288                  ! velocity * 1/2
289                  zhw = 0.5 * zwn(ji,jj,jk)
290                  ! upstream scheme
291                  zfp_w = zhw + ABS( zhw )
292                  zfm_w = zhw - ABS( zhw )
293                  zupst = zfp_w * trb(ji,jj,jk,jn) + zfm_w * trb(ji,jj,jk-1,jn)
294                  ! centered scheme
295                  zcent = zhw * ( trn(ji,jj,jk,jn) + trn(ji,jj,jk-1,jn) )
296                  ! centered scheme
297                  zwx(ji,jj,jk) = zcofk * zupst + (1.-zcofk) * zcent
298               END DO
299            END DO
300         END DO
301
302
303         ! 2. Tracer flux divergence at t-point added to the general trend
304         ! -------------------------
305
306         DO jk = 1, jpkm1
307            DO jj = 2, jpjm1
308               DO ji = fs_2, fs_jpim1   ! vector opt.
309                  ze3tr = 1. / fse3t(ji,jj,jk)
310                  ! vertical advective trends
311                  ztra = - ze3tr * ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) )
312                  ! add it to the general tracer trends
313                  tra(ji,jj,jk,jn) =  tra(ji,jj,jk,jn) + ztra
314#if defined key_trc_diatrd 
315                  ! save the vertical advective trends computed as w gradz(T)
316                  trtrd(ji,jj,jk,jn,3) = ztra - trn(ji,jj,jk,jn) * hdivn(ji,jj,jk)
317#endif
318               END DO
319            END DO
320         END DO
321
322         IF(l_ctl) THEN         ! print mean trends (used for debugging)
323            ztra = SUM( tra(2:nictl,2:njctl,1:jpkm1,jn) * tmask(2:nictl,2:njctl,1:jpkm1) )
324            WRITE(numout,*) ' trc/zad - ',ctrcnm(jn),' : ', ztra-tra_ctl(jn), ' centered2' 
325            tra_ctl(jn) = ztra 
326         ENDIF
327
328      END DO
329
330
331   END SUBROUTINE trc_adv_cen2
332#else
333
334   !!----------------------------------------------------------------------
335   !!   Default option                                         Empty module
336   !!----------------------------------------------------------------------
337CONTAINS
338   SUBROUTINE trc_adv_cen2( kt ) 
339      INTEGER, INTENT(in) :: kt
340      WRITE(*,*) 'trc_adv_cen2: You should not have seen this print! error?', kt
341   END SUBROUTINE trc_adv_cen2
342#endif
343   !!======================================================================
344END MODULE trcadv_cen2
Note: See TracBrowser for help on using the repository browser.