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 branches/dev_001_GM/NEMO/TOP_SRC/TRP – NEMO

source: branches/dev_001_GM/NEMO/TOP_SRC/TRP/trcadv_cen2.F90 @ 776

Last change on this file since 776 was 776, checked in by gm, 16 years ago

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