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

Last change on this file since 941 was 941, checked in by cetlod, 16 years ago

phasing the passive tracer transport module to the new version of NEMO, see ticket 143
  • Property svn:executable set to *
File size: 18.6 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   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   upsmsk    !: mixed upstream/centered scheme near some straits
25   !                                                   !  and in closed seas (orca 2 and 4 configurations)
26
27   !! * Substitutions
28#  include "top_substitute.h90"
29   !!----------------------------------------------------------------------
30   !!   TOP 1.0 , LOCEAN-IPSL (2005)
31   !! $Header$
32   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
33   !!----------------------------------------------------------------------
34
35CONTAINS
36
37   !!----------------------------------------------------------------------
38   !!   Default option :             2nd order centered scheme (k-j-i loop)
39   !!----------------------------------------------------------------------
40
41   SUBROUTINE trc_adv_cen2( kt )
42      !!----------------------------------------------------------------------
43      !!                  ***  ROUTINE trc_adv_cen2  ***
44      !!                 
45      !! ** Purpose :   Compute the now trend due to the advection of tracers
46      !!      and add it to the general trend of passive tracer equations.
47      !!
48      !! ** Method  :   The advection is evaluated by a second order centered
49      !!      scheme using now fields (leap-frog scheme). In specific areas
50      !!      (vicinity of major river mouths, some straits, or where tn is
51      !!      is approaching the freezing point) it is mixed with an upstream
52      !!      scheme for stability reasons.
53      !!        Part 0 : compute the upstream / centered flag
54      !!                 (3D array, zind, defined at T-point (0<zind<1))
55      !!        Part I : horizontal advection
56      !!      * centered flux:
57      !!         * s-coordinate (ln_sco=T) or
58      !!         * z-coordinate with partial steps (ln_zps=T),
59      !!        the vertical scale factors e3. are inside the derivatives:
60      !!               zcenu = e2u*e3u  un  mi(tn)
61      !!               zcenv = e1v*e3v  vn  mj(tn)
62      !!         * z-coordinate (default key), e3t=e3u=e3v:
63      !!               zcenu = e2u  un  mi(tn)
64      !!               zcenv = e1v  vn  mj(tn)
65      !!      * horizontal advective trend (divergence of the fluxes)
66      !!         * s-coordinate (ln_sco=T) or
67      !!         * z-coordinate with partial steps (ln_zps=T)
68      !!               ztra = 1/(e1t*e2t*e3t) { di-1[zwx] + dj-1[zwy] }
69      !!         * z-coordinate (default key), e3t=e3u=e3v:
70      !!               ztra = 1/(e1t*e2t) { di-1[zwx] + dj-1[zwy] }
71      !!      * Add this trend now to the general trend of tracer tra:
72      !!              tra = tra + ztra
73      !!      * trend diagnostic ('key_trc_diatrd'): the trend is saved
74      !!      for diagnostics. The trends saved is expressed as
75      !!      Uh.gradh(T)
76      !!
77      !!         Part II : vertical advection
78      !!      For any tracer  the advective trend is computed as follows :
79      !!            ztra = 1/e3t dk+1[ zwz ]
80      !!      where the vertical advective flux, zwz, is given by :
81      !!            zwz = zcofk * zupst + (1-zcofk) * zcent
82      !!      with
83      !!        zupsv = upstream flux = wn * (trb(k) or trb(k-1) ) [wn>0 or <0]
84      !!        zcenu = centered flux = wn * mk(trn)
85      !!         The surface boundary condition is :
86      !!      rigid-lid (default option) : zero advective flux
87      !!      free-surf ("key_fresurf_cstvol") : wn(:,:,1) * trn(:,:,1)
88      !!         Add this trend now to the general trend of tracer tra :
89      !!            tra = tra + ztra
90      !!         Trend diagnostic ('key_trc_diatrd'): the trend is saved for
91      !!      diagnostics. The trends saved is expressed as :
92      !!             save trend =  w.gradz(T) = ztra - trn divn.
93      !!
94      !! ** Action : - update tra with the now advective tracer trends
95      !!             - save the trends in trtrd ('key_trc_diatrd')
96      !!
97      !! History :
98      !!   8.2  !  01-08  (M-A Filiberti, and M.Levy)  trahad+trazad = traadv
99      !!   8.5  !  02-06  (G. Madec, C. Ethe)  F90: Free form and module
100      !!----------------------------------------------------------------------
101      !! * Modules used
102      USE oce_trc            , zwx => ua,  &  ! use ua as workspace
103         &                     zwy => va      ! use va as workspace
104#if defined key_trcbbl_adv
105      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  &  ! temporary arrays
106         &         zun, zvn, zwn
107#else
108      USE oce_trc            , zun => un,  &  ! When no bbl, zun == un
109         &                     zvn => vn,  &  ! When no bbl, zvn == vn
110         &                     zwn => wn      ! When no bbl, zwn == wn
111#endif
112
113 
114      !! * Arguments
115      INTEGER, INTENT( in ) ::   kt           ! ocean time-step index
116 
117      !! * Local save
118      REAL(wp), DIMENSION(jpi,jpj), SAVE ::   &
119         zbtr2
120 
121      !! * Local declarations
122      INTEGER  ::   ji, jj, jk, jn             ! dummy loop indices
123      REAL(wp) ::                           &
124         zbtr, ztra, zfui, zfvj,            &  ! temporary scalars
125         zhw, ze3tr, zcofi, zcofj,          &  !    "         "
126         zupsut, zupsvt,                    &  !    "         "
127         zfp_ui, zfp_vj, zfm_ui, zfm_vj,    &  !    "         "
128         zcofk, zupst, zcent,               &  !    "         "
129         zfp_w, zfm_w,                      &  !    "         "
130         zcenut, zcenvt                        !
131
132      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &
133         zind                              ! temporary workspace arrays
134#if defined key_trc_diatrd
135      REAL(wp) ::                           &
136         ztai, ztaj,                        &  ! temporary scalars
137         zfui1, zfvj1                          !    "         "
138#endif
139      CHARACTER (len=22) :: charout
140      !!----------------------------------------------------------------------
141
142      IF( kt == nittrc000 ) THEN
143         IF(lwp) WRITE(numout,*)
144         IF(lwp) WRITE(numout,*) 'trc_adv_cen2 : 2nd order centered advection scheme'
145         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~   Vector optimization case'
146         IF(lwp) WRITE(numout,*)
147 
148         upsmsk(:,:) = 0.e0                              ! not upstream by default
149         IF( cp_cfg == "orca" )   CALL ups_orca_set      ! set mixed Upstream/centered scheme near some straits
150         !                                               ! and in closed seas (orca2 and orca4 only)
151         zbtr2(:,:) = 1. / ( e1t(:,:) * e2t(:,:) )
152      ENDIF
153
154#if defined key_trcbbl_adv
155
156      ! Advective bottom boundary layer
157      ! -------------------------------
158      zun(:,:,:) = un(:,:,:) - u_trc_bbl(:,:,:)
159      zvn(:,:,:) = vn(:,:,:) - v_trc_bbl(:,:,:)
160      zwn(:,:,:) = wn(:,:,:) + w_trc_bbl(:,:,:)
161#endif
162
163      ! Upstream / centered scheme indicator
164      ! ------------------------------------
165      DO jk = 1, jpk
166         DO jj = 1, jpj
167            DO ji = 1, jpi
168               zind(ji,jj,jk) = MAX (   &
169                  rnfmsk(ji,jj) * rnfmsk_z(jk),      &  ! near runoff mouths (& closed sea outflows)
170                  upsmsk(ji,jj)                      &  ! some of some straits
171#if defined key_lim3 || defined key_lim2
172                  !                                     ! below ice covered area (if tn < "freezing"+0.1 )
173                  , MAX(  0., SIGN( 1., fzptn(ji,jj) + 0.1 - tn(ji,jj,jk) )  ) * tmask(ji,jj,jk)   &
174#endif
175                  &                  )
176            END DO
177         END DO
178      END DO
179
180
181
182      DO jn = 1, jptra
183         ! I. Horizontal advective fluxes
184         ! ------------------------------
185         
186         ! Second order centered tracer flux at u and v-points
187         
188         !                                                ! ===============
189         DO jk = 1, jpkm1                                 ! Horizontal slab
190            !                                             ! ===============
191            DO jj = 1, jpjm1
192               DO ji = 1, fs_jpim1   ! vector opt.
193               ! upstream indicator
194                  zcofi = MAX( zind(ji+1,jj,jk), zind(ji,jj,jk) )
195                  zcofj = MAX( zind(ji,jj+1,jk), zind(ji,jj,jk) )
196                  ! volume fluxes * 1/2
197#if ! defined key_zco
198                  zfui = 0.5 * e2u(ji,jj) * fse3u(ji,jj,jk) * zun(ji,jj,jk)
199                  zfvj = 0.5 * e1v(ji,jj) * fse3v(ji,jj,jk) * zvn(ji,jj,jk)
200#else
201                  zfui = 0.5 * e2u(ji,jj) * zun(ji,jj,jk)
202                  zfvj = 0.5 * e1v(ji,jj) * zvn(ji,jj,jk)
203#endif
204               ! upstream scheme
205                  zfp_ui = zfui + ABS( zfui )
206                  zfp_vj = zfvj + ABS( zfvj )
207                  zfm_ui = zfui - ABS( zfui )
208                  zfm_vj = zfvj - ABS( zfvj )
209                  zupsut = zfp_ui * trb(ji,jj,jk,jn) + zfm_ui * trb(ji+1,jj  ,jk,jn)
210                  zupsvt = zfp_vj * trb(ji,jj,jk,jn) + zfm_vj * trb(ji  ,jj+1,jk,jn)
211                  ! centered scheme
212                  zcenut = zfui * ( trn(ji,jj,jk,jn) + trn(ji+1,jj  ,jk,jn) )
213                  zcenvt = zfvj * ( trn(ji,jj,jk,jn) + trn(ji  ,jj+1,jk,jn) )
214                  ! mixed centered / upstream scheme
215                  zwx(ji,jj,jk) = zcofi * zupsut + (1.-zcofi) * zcenut
216                  zwy(ji,jj,jk) = zcofj * zupsvt + (1.-zcofj) * zcenvt               
217               END DO
218            END DO
219
220
221            ! 2. Tracer flux divergence at t-point added to the general trend
222            ! -------------------------
223
224            DO jj = 2, jpjm1
225               DO ji = fs_2, fs_jpim1   ! vector opt.
226#if ! defined key_zco
227                  zbtr = zbtr2(ji,jj) / fse3t(ji,jj,jk)
228#else
229                  zbtr = zbtr2(ji,jj) 
230#endif
231                  ! horizontal advective trends
232                  ztra = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk)   &
233                     &             + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk)  )
234
235                  ! add it to the general tracer trends
236                  tra(ji,jj,jk,jn) = tra(ji,jj,jk,jn) + ztra
237
238#if defined key_trc_diatrd 
239                  ! recompute the trends in i- and j-direction as Uh gradh(T)
240#if ! defined key_zco
241                  zfui = 0.5 * e2u(ji  ,jj) * fse3u(ji,  jj,jk) * zun(ji,  jj,jk)
242                  zfui1= 0.5 * e2u(ji-1,jj) * fse3u(ji-1,jj,jk) * zun(ji-1,jj,jk)
243                  zfvj = 0.5 * e1v(ji,jj  ) * fse3v(ji,jj  ,jk) * zvn(ji,jj  ,jk)
244                  zfvj1= 0.5 * e1v(ji,jj-1) * fse3v(ji,jj-1,jk) * zvn(ji,jj-1,jk)
245# else
246                  zfui = 0.5 * e2u(ji  ,jj) * zun(ji,  jj,jk)
247                  zfui1= 0.5 * e2u(ji-1,jj) * zun(ji-1,jj,jk)
248                  zfvj = 0.5 * e1v(ji,jj  ) * zvn(ji,jj  ,jk)
249                  zfvj1= 0.5 * e1v(ji,jj-1) * zvn(ji,jj-1,jk)
250# endif
251                  ztai = - zbtr * ( zfui  * ( trn(ji+1,jj  ,jk,jn) - trn(ji,  jj,jk,jn) )   &
252                     &                + zfui1 * ( trn(ji,  jj,  jk,jn) - trn(ji-1,jj,jk,jn) ) )
253                  ztaj = - zbtr * ( zfvj  * ( trn(ji  ,jj+1,jk,jn) - trn(ji,jj  ,jk,jn) )    &
254                     &                + zfvj1 * ( trn(ji  ,jj  ,jk,jn) - trn(ji,jj-1,jk,jn) ) )
255                  ! save i- and j- advective trends computed as Uh gradh(T)
256                  IF (luttrd(jn)) trtrd(ji,jj,jk,ikeep(jn),1) = ztai
257                  IF (luttrd(jn)) trtrd(ji,jj,jk,ikeep(jn),2) = ztaj
258#endif
259               END DO
260            END DO
261            !                                             ! ===============
262         END DO                                           !   End of slab
263         !                                                ! ===============
264      ENDDO
265
266      IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
267         WRITE(charout, FMT="('centered2 - had')")
268         CALL prt_ctl_trc_info(charout)
269         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm,clinfo2='trd')
270      ENDIF
271     
272      ! II. Vertical advection
273      ! ----------------------
274      DO jn = 1, jptra
275
276         ! Bottom value : flux set to zero
277         zwx(:,:,jpk) = 0.e0 
278
279         ! Surface value
280         IF ( lk_dynspg_rl ) THEN       ! rigid lid : flux set to zero
281            zwx(:,:, 1 ) = 0.e0 
282         ELSE                           ! free surface-constant volume
283            zwx(:,:, 1 ) = zwn(:,:,1) * trn(:,:,1,jn)
284         ENDIF
285
286         ! 1. Vertical advective fluxes
287         ! ----------------------------
288
289         ! Second order centered tracer flux at w-point
290
291         DO jk = 2, jpk
292            DO jj = 2, jpjm1
293               DO ji = fs_2, fs_jpim1   ! vector opt.
294                  ! upstream indicator
295                  zcofk = MAX( zind(ji,jj,jk-1), zind(ji,jj,jk) )
296                  ! velocity * 1/2
297                  zhw = 0.5 * zwn(ji,jj,jk)
298                  ! upstream scheme
299                  zfp_w = zhw + ABS( zhw )
300                  zfm_w = zhw - ABS( zhw )
301                  zupst = zfp_w * trb(ji,jj,jk,jn) + zfm_w * trb(ji,jj,jk-1,jn)
302                  ! centered scheme
303                  zcent = zhw * ( trn(ji,jj,jk,jn) + trn(ji,jj,jk-1,jn) )
304                  ! centered scheme
305                  zwx(ji,jj,jk) = zcofk * zupst + (1.-zcofk) * zcent
306               END DO
307            END DO
308         END DO
309
310
311         ! 2. Tracer flux divergence at t-point added to the general trend
312         ! -------------------------
313
314         DO jk = 1, jpkm1
315            DO jj = 2, jpjm1
316               DO ji = fs_2, fs_jpim1   ! vector opt.
317                  ze3tr = 1. / fse3t(ji,jj,jk)
318                  ! vertical advective trends
319                  ztra = - ze3tr * ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) )
320                  ! add it to the general tracer trends
321                  tra(ji,jj,jk,jn) =  tra(ji,jj,jk,jn) + ztra
322#if defined key_trc_diatrd 
323                  ! save the vertical advective trends computed as w gradz(T)
324                  IF (luttrd(jn)) trtrd(ji,jj,jk,ikeep(jn),3) = ztra - trn(ji,jj,jk,jn) * hdivn(ji,jj,jk)
325#endif
326               END DO
327            END DO
328         END DO
329
330      END DO
331
332      IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
333         WRITE(charout, FMT="('centered - zad')")
334         CALL prt_ctl_trc_info(charout)
335         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm,clinfo2='trd')
336      ENDIF
337
338   END SUBROUTINE trc_adv_cen2
339
340   SUBROUTINE ups_orca_set
341      !!----------------------------------------------------------------------
342      !!                  ***  ROUTINE ups_orca_set  ***
343      !!
344      !! ** Purpose :   add a portion of upstream scheme in area where the
345      !!                centered scheme generates too strong overshoot
346      !!
347      !! ** Method  :   orca (R4 and R2) confiiguration setting. Set upsmsk
348      !!                array to nozero value in some straith.
349      !!
350      !! ** Action : - upsmsk set to 1 at some strait, 0 elsewhere for orca
351      !!----------------------------------------------------------------------
352      INTEGER  ::   ii0, ii1, ij0, ij1      ! temporary integers
353      !!----------------------------------------------------------------------
354
355      ! mixed upstream/centered scheme near river mouths
356      ! ------------------------------------------------
357      SELECT CASE ( jp_cfg )
358      !                                        ! =======================
359      CASE ( 4 )                               !  ORCA_R4 configuration
360         !                                     ! =======================
361         !                                          ! Gibraltar Strait
362         ii0 =  70   ;   ii1 =  71
363         ij0 =  52   ;   ij1 =  53   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50
364         !
365         !                                     ! =======================
366      CASE ( 2 )                               !  ORCA_R2 configuration
367         !                                     ! =======================
368         !                                          ! Gibraltar Strait
369         ij0 = 102   ;   ij1 = 102
370         ii0 = 138   ;   ii1 = 138   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.20
371         ii0 = 139   ;   ii1 = 139   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.40
372         ii0 = 140   ;   ii1 = 140   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50
373         ij0 = 101   ;   ij1 = 102
374         ii0 = 141   ;   ii1 = 141   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50
375         !                                          ! Bab el Mandeb Strait
376         ij0 =  87   ;   ij1 =  88
377         ii0 = 164   ;   ii1 = 164   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.10
378         ij0 =  88   ;   ij1 =  88
379         ii0 = 163   ;   ii1 = 163   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.25
380         ii0 = 162   ;   ii1 = 162   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.40
381         ii0 = 160   ;   ii1 = 161   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50
382         ij0 =  89   ;   ij1 =  89
383         ii0 = 158   ;   ii1 = 160   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.25
384         ij0 =  90   ;   ij1 =  90
385         ii0 = 160   ;   ii1 = 160   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.25
386         !                                          ! Sound Strait
387         ij0 = 116   ;   ij1 = 116
388         ii0 = 145   ;   ii1 = 147   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50
389         !
390      END SELECT
391
392   END SUBROUTINE ups_orca_set
393
394#else
395
396   !!----------------------------------------------------------------------
397   !!   Default option                                         Empty module
398   !!----------------------------------------------------------------------
399CONTAINS
400   SUBROUTINE trc_adv_cen2( kt ) 
401      INTEGER, INTENT(in) :: kt
402      WRITE(*,*) 'trc_adv_cen2: You should not have seen this print! error?', kt
403   END SUBROUTINE trc_adv_cen2
404#endif
405   !!======================================================================
406END MODULE trcadv_cen2
Note: See TracBrowser for help on using the repository browser.