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

Last change on this file since 910 was 910, checked in by ctlod, 16 years ago

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