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

Last change on this file since 204 was 204, checked in by opalod, 19 years ago

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