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

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

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