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

Last change on this file since 1152 was 1152, checked in by rblod, 16 years ago

Convert cvs header to svn Id, step II
  • Property svn:executable set to *
  • Property svn:keywords set to Id
File size: 19.2 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 trp_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   !! $Id$
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#if defined key_lim3 || defined key_lim2
140      REAL(wp) ::                           &
141         ztfreez                               ! freezing point
142#endif
143      CHARACTER (len=22) :: charout
144      !!----------------------------------------------------------------------
145
146      IF( kt == nittrc000 ) THEN
147         IF(lwp) WRITE(numout,*)
148         IF(lwp) WRITE(numout,*) 'trc_adv_cen2 : 2nd order centered advection scheme'
149         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~   Vector optimization case'
150         IF(lwp) WRITE(numout,*)
151 
152         upsmsk(:,:) = 0.e0                              ! not upstream by default
153         IF( cp_cfg == "orca" )   CALL ups_orca_set      ! set mixed Upstream/centered scheme near some straits
154         !                                               ! and in closed seas (orca2 and orca4 only)
155         zbtr2(:,:) = 1. / ( e1t(:,:) * e2t(:,:) )
156      ENDIF
157
158#if defined key_trcbbl_adv
159
160      ! Advective bottom boundary layer
161      ! -------------------------------
162      zun(:,:,:) = un(:,:,:) - u_trc_bbl(:,:,:)
163      zvn(:,:,:) = vn(:,:,:) - v_trc_bbl(:,:,:)
164      zwn(:,:,:) = wn(:,:,:) + w_trc_bbl(:,:,:)
165#endif
166
167      ! Upstream / centered scheme indicator
168      ! ------------------------------------
169      DO jk = 1, jpk
170         DO jj = 1, jpj
171            DO ji = 1, jpi
172#if defined key_lim3 || defined key_lim2
173               ztfreez = ( - 0.0575 + 1.710523e-3 * SQRT( sn(ji,jj,1) )   &
174                 &                  - 2.154996e-4 *       sn(ji,jj,1)   ) * sn(ji,jj,1)
175
176               zind(ji,jj,jk) = MAX (   &
177                  rnfmsk(ji,jj) * rnfmsk_z(jk),      &  ! near runoff mouths (& closed sea outflows)
178                  upsmsk(ji,jj)                      &  ! some of some straits
179                  !                                     ! below ice covered area (if tn < "freezing"+0.1 )
180                  , MAX(  0., SIGN( 1., ztfreez + 0.1 - tn(ji,jj,jk) )  ) * tmask(ji,jj,jk)   &
181                  &                  )
182
183#else
184               zind(ji,jj,jk) = MAX (   &
185                  rnfmsk(ji,jj) * rnfmsk_z(jk),      &  ! near runoff mouths (& closed sea outflows)
186                  upsmsk(ji,jj)                      &  ! some of some straits
187                  &                  )
188#endif
189            END DO
190         END DO
191      END DO
192
193
194
195      DO jn = 1, jptra
196         ! I. Horizontal advective fluxes
197         ! ------------------------------
198         
199         ! Second order centered tracer flux at u and v-points
200         
201         !                                                ! ===============
202         DO jk = 1, jpkm1                                 ! Horizontal slab
203            !                                             ! ===============
204            DO jj = 1, jpjm1
205               DO ji = 1, fs_jpim1   ! vector opt.
206               ! upstream indicator
207                  zcofi = MAX( zind(ji+1,jj,jk), zind(ji,jj,jk) )
208                  zcofj = MAX( zind(ji,jj+1,jk), zind(ji,jj,jk) )
209                  ! volume fluxes * 1/2
210#if ! defined key_zco
211                  zfui = 0.5 * e2u(ji,jj) * fse3u(ji,jj,jk) * zun(ji,jj,jk)
212                  zfvj = 0.5 * e1v(ji,jj) * fse3v(ji,jj,jk) * zvn(ji,jj,jk)
213#else
214                  zfui = 0.5 * e2u(ji,jj) * zun(ji,jj,jk)
215                  zfvj = 0.5 * e1v(ji,jj) * zvn(ji,jj,jk)
216#endif
217               ! upstream scheme
218                  zfp_ui = zfui + ABS( zfui )
219                  zfp_vj = zfvj + ABS( zfvj )
220                  zfm_ui = zfui - ABS( zfui )
221                  zfm_vj = zfvj - ABS( zfvj )
222                  zupsut = zfp_ui * trb(ji,jj,jk,jn) + zfm_ui * trb(ji+1,jj  ,jk,jn)
223                  zupsvt = zfp_vj * trb(ji,jj,jk,jn) + zfm_vj * trb(ji  ,jj+1,jk,jn)
224                  ! centered scheme
225                  zcenut = zfui * ( trn(ji,jj,jk,jn) + trn(ji+1,jj  ,jk,jn) )
226                  zcenvt = zfvj * ( trn(ji,jj,jk,jn) + trn(ji  ,jj+1,jk,jn) )
227                  ! mixed centered / upstream scheme
228                  zwx(ji,jj,jk) = zcofi * zupsut + (1.-zcofi) * zcenut
229                  zwy(ji,jj,jk) = zcofj * zupsvt + (1.-zcofj) * zcenvt               
230               END DO
231            END DO
232
233
234            ! 2. Tracer flux divergence at t-point added to the general trend
235            ! -------------------------
236
237            DO jj = 2, jpjm1
238               DO ji = fs_2, fs_jpim1   ! vector opt.
239#if ! defined key_zco
240                  zbtr = zbtr2(ji,jj) / fse3t(ji,jj,jk)
241#else
242                  zbtr = zbtr2(ji,jj) 
243#endif
244                  ! horizontal advective trends
245                  ztra = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk)   &
246                     &             + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk)  )
247
248                  ! add it to the general tracer trends
249                  tra(ji,jj,jk,jn) = tra(ji,jj,jk,jn) + ztra
250
251#if defined key_trc_diatrd 
252                  ! recompute the trends in i- and j-direction as Uh gradh(T)
253#if ! defined key_zco
254                  zfui = 0.5 * e2u(ji  ,jj) * fse3u(ji,  jj,jk) * zun(ji,  jj,jk)
255                  zfui1= 0.5 * e2u(ji-1,jj) * fse3u(ji-1,jj,jk) * zun(ji-1,jj,jk)
256                  zfvj = 0.5 * e1v(ji,jj  ) * fse3v(ji,jj  ,jk) * zvn(ji,jj  ,jk)
257                  zfvj1= 0.5 * e1v(ji,jj-1) * fse3v(ji,jj-1,jk) * zvn(ji,jj-1,jk)
258# else
259                  zfui = 0.5 * e2u(ji  ,jj) * zun(ji,  jj,jk)
260                  zfui1= 0.5 * e2u(ji-1,jj) * zun(ji-1,jj,jk)
261                  zfvj = 0.5 * e1v(ji,jj  ) * zvn(ji,jj  ,jk)
262                  zfvj1= 0.5 * e1v(ji,jj-1) * zvn(ji,jj-1,jk)
263# endif
264                  ztai = - zbtr * ( zfui  * ( trn(ji+1,jj  ,jk,jn) - trn(ji,  jj,jk,jn) )   &
265                     &                + zfui1 * ( trn(ji,  jj,  jk,jn) - trn(ji-1,jj,jk,jn) ) )
266                  ztaj = - zbtr * ( zfvj  * ( trn(ji  ,jj+1,jk,jn) - trn(ji,jj  ,jk,jn) )    &
267                     &                + zfvj1 * ( trn(ji  ,jj  ,jk,jn) - trn(ji,jj-1,jk,jn) ) )
268                  ! save i- and j- advective trends computed as Uh gradh(T)
269                  IF (luttrd(jn)) trtrd(ji,jj,jk,ikeep(jn),1) = ztai
270                  IF (luttrd(jn)) trtrd(ji,jj,jk,ikeep(jn),2) = ztaj
271#endif
272               END DO
273            END DO
274            !                                             ! ===============
275         END DO                                           !   End of slab
276         !                                                ! ===============
277      ENDDO
278
279      IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
280         WRITE(charout, FMT="('centered2 - had')")
281         CALL prt_ctl_trc_info(charout)
282         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm,clinfo2='trd')
283      ENDIF
284     
285      ! II. Vertical advection
286      ! ----------------------
287      DO jn = 1, jptra
288
289         ! Bottom value : flux set to zero
290         zwx(:,:,jpk) = 0.e0 
291
292         ! Surface value
293         IF ( lk_dynspg_rl ) THEN       ! rigid lid : flux set to zero
294            zwx(:,:, 1 ) = 0.e0 
295         ELSE                           ! free surface-constant volume
296            zwx(:,:, 1 ) = zwn(:,:,1) * trn(:,:,1,jn)
297         ENDIF
298
299         ! 1. Vertical advective fluxes
300         ! ----------------------------
301
302         ! Second order centered tracer flux at w-point
303
304         DO jk = 2, jpk
305            DO jj = 2, jpjm1
306               DO ji = fs_2, fs_jpim1   ! vector opt.
307                  ! upstream indicator
308                  zcofk = MAX( zind(ji,jj,jk-1), zind(ji,jj,jk) )
309                  ! velocity * 1/2
310                  zhw = 0.5 * zwn(ji,jj,jk)
311                  ! upstream scheme
312                  zfp_w = zhw + ABS( zhw )
313                  zfm_w = zhw - ABS( zhw )
314                  zupst = zfp_w * trb(ji,jj,jk,jn) + zfm_w * trb(ji,jj,jk-1,jn)
315                  ! centered scheme
316                  zcent = zhw * ( trn(ji,jj,jk,jn) + trn(ji,jj,jk-1,jn) )
317                  ! centered scheme
318                  zwx(ji,jj,jk) = zcofk * zupst + (1.-zcofk) * zcent
319               END DO
320            END DO
321         END DO
322
323
324         ! 2. Tracer flux divergence at t-point added to the general trend
325         ! -------------------------
326
327         DO jk = 1, jpkm1
328            DO jj = 2, jpjm1
329               DO ji = fs_2, fs_jpim1   ! vector opt.
330                  ze3tr = 1. / fse3t(ji,jj,jk)
331                  ! vertical advective trends
332                  ztra = - ze3tr * ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) )
333                  ! add it to the general tracer trends
334                  tra(ji,jj,jk,jn) =  tra(ji,jj,jk,jn) + ztra
335#if defined key_trc_diatrd 
336                  ! save the vertical advective trends computed as w gradz(T)
337                  IF (luttrd(jn)) trtrd(ji,jj,jk,ikeep(jn),3) = ztra - trn(ji,jj,jk,jn) * hdivn(ji,jj,jk)
338#endif
339               END DO
340            END DO
341         END DO
342
343      END DO
344
345      IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
346         WRITE(charout, FMT="('centered - zad')")
347         CALL prt_ctl_trc_info(charout)
348         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm,clinfo2='trd')
349      ENDIF
350
351   END SUBROUTINE trc_adv_cen2
352
353   SUBROUTINE ups_orca_set
354      !!----------------------------------------------------------------------
355      !!                  ***  ROUTINE ups_orca_set  ***
356      !!
357      !! ** Purpose :   add a portion of upstream scheme in area where the
358      !!                centered scheme generates too strong overshoot
359      !!
360      !! ** Method  :   orca (R4 and R2) confiiguration setting. Set upsmsk
361      !!                array to nozero value in some straith.
362      !!
363      !! ** Action : - upsmsk set to 1 at some strait, 0 elsewhere for orca
364      !!----------------------------------------------------------------------
365      INTEGER  ::   ii0, ii1, ij0, ij1      ! temporary integers
366      !!----------------------------------------------------------------------
367
368      ! mixed upstream/centered scheme near river mouths
369      ! ------------------------------------------------
370      SELECT CASE ( jp_cfg )
371      !                                        ! =======================
372      CASE ( 4 )                               !  ORCA_R4 configuration
373         !                                     ! =======================
374         !                                          ! Gibraltar Strait
375         ii0 =  70   ;   ii1 =  71
376         ij0 =  52   ;   ij1 =  53   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50
377         !
378         !                                     ! =======================
379      CASE ( 2 )                               !  ORCA_R2 configuration
380         !                                     ! =======================
381         !                                          ! Gibraltar Strait
382         ij0 = 102   ;   ij1 = 102
383         ii0 = 138   ;   ii1 = 138   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.20
384         ii0 = 139   ;   ii1 = 139   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.40
385         ii0 = 140   ;   ii1 = 140   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50
386         ij0 = 101   ;   ij1 = 102
387         ii0 = 141   ;   ii1 = 141   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50
388         !                                          ! Bab el Mandeb Strait
389         ij0 =  87   ;   ij1 =  88
390         ii0 = 164   ;   ii1 = 164   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.10
391         ij0 =  88   ;   ij1 =  88
392         ii0 = 163   ;   ii1 = 163   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.25
393         ii0 = 162   ;   ii1 = 162   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.40
394         ii0 = 160   ;   ii1 = 161   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50
395         ij0 =  89   ;   ij1 =  89
396         ii0 = 158   ;   ii1 = 160   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.25
397         ij0 =  90   ;   ij1 =  90
398         ii0 = 160   ;   ii1 = 160   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.25
399         !                                          ! Sound Strait
400         ij0 = 116   ;   ij1 = 116
401         ii0 = 145   ;   ii1 = 147   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50
402         !
403      END SELECT
404
405   END SUBROUTINE ups_orca_set
406
407#else
408
409   !!----------------------------------------------------------------------
410   !!   Default option                                         Empty module
411   !!----------------------------------------------------------------------
412CONTAINS
413   SUBROUTINE trc_adv_cen2( kt ) 
414      INTEGER, INTENT(in) :: kt
415      WRITE(*,*) 'trc_adv_cen2: You should not have seen this print! error?', kt
416   END SUBROUTINE trc_adv_cen2
417#endif
418   !!======================================================================
419END MODULE trcadv_cen2
Note: See TracBrowser for help on using the repository browser.