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

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

nemo_v1_update_022 : CE + RB + CT : add print control possibility
  • 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   USE prtctl_trc
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      CHARACTER (len=22) :: charout
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      ! 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      DO jn = 1, jptra
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      ENDDO
258
259      IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
260         WRITE(charout, FMT="('centered2 - had')")
261         CALL prt_ctl_trc_info(charout)
262         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm,clinfo2='trd')
263      ENDIF
264     
265      ! II. Vertical advection
266      ! ----------------------
267      DO jn = 1, jptra
268
269         ! Bottom value : flux set to zero
270         zwx(:,:,jpk) = 0.e0 
271
272         ! Surface value
273         IF( lk_dynspg_fsc ) THEN       ! free surface-constant volume
274            zwx(:,:, 1 ) = zwn(:,:,1) * trn(:,:,1,jn)
275         ELSE                           ! 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      END DO
324
325      IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
326         WRITE(charout, FMT="('centered - zad')")
327         CALL prt_ctl_trc_info(charout)
328         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm,clinfo2='trd')
329      ENDIF
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.