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.
traadv_cen2.F90 in branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA – NEMO

source: branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/traadv_cen2.F90 @ 2024

Last change on this file since 2024 was 2024, checked in by cetlod, 14 years ago

Merge of active and passive tracer advection/diffusion modules, see ticket:693
  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 18.1 KB
Line 
1MODULE traadv_cen2
2   !!======================================================================
3   !!                     ***  MODULE  traadv_cen2  ***
4   !! Ocean  tracers:  horizontal & vertical advective trend
5   !!======================================================================
6   !! History :  8.2  ! 2001-08  (G. Madec, E. Durand)  trahad+trazad=traadv
7   !!            1.0  ! 2002-06  (G. Madec)  F90: Free form and module
8   !!            9.0  ! 2004-08  (C. Talandier) New trends organization
9   !!             -   ! 2005-11  (V. Garnier) Surface pressure gradient organization
10   !!            2.0  ! 2006-04  (R. Benshila, G. Madec) Step reorganization
11   !!             -   ! 2006-07  (G. madec)  add ups_orca_set routine
12   !!            3.2  ! 2009-07  (G. Madec) add avmb, avtb in restart for cen2 advection
13   !!            3.3  ! 2010-05  (C. Ethe, G. Madec)  merge TRC-TRA + switch from velocity to transport
14   !!----------------------------------------------------------------------
15
16   !!----------------------------------------------------------------------
17   !!   tra_adv_cen2 : update the tracer trend with the horizontal and
18   !!                  vertical advection trends using a seconder order
19   !!   ups_orca_set : allow mixed upstream/centered scheme in specific
20   !!                  area (set for orca 2 and 4 only)
21   !!----------------------------------------------------------------------
22   USE oce, ONLY: tsn  ! now ocean temperature and salinity
23   USE dom_oce         ! ocean space and time domain
24   USE eosbn2          ! equation of state
25   USE trdmod_oce      ! tracers trends
26   USE trdtra          ! tracers trends
27   USE closea          ! closed sea
28   USE sbcrnf          ! river runoffs
29   USE in_out_manager  ! I/O manager
30   USE iom             ! IOM library
31   USE diaptr          ! poleward transport diagnostics
32   USE zdf_oce         ! ocean vertical physics
33   USE restart         ! ocean restart
34
35   IMPLICIT NONE
36   PRIVATE
37
38   PUBLIC   tra_adv_cen2    ! routine called by step.F90
39   PUBLIC   ups_orca_set    ! routine used by traadv_cen2_jki.F90
40
41   LOGICAL  :: l_trd       ! flag to compute trends
42
43   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: upsmsk    !: mixed upstream/centered scheme near some straits
44   !                                                   !  and in closed seas (orca 2 and 4 configurations)
45   !! * Substitutions
46#  include "domzgr_substitute.h90"
47#  include "vectopt_loop_substitute.h90"
48   !!----------------------------------------------------------------------
49   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009)
50   !! $Id$
51   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
52   !!----------------------------------------------------------------------
53
54CONTAINS
55
56   SUBROUTINE tra_adv_cen2( kt   , cdtype, pun  , pvn, pwn, &
57      &                     ptrab, ptran , ptraa, kjpt   ) 
58      !!----------------------------------------------------------------------
59      !!                  ***  ROUTINE tra_adv_cen2  ***
60      !!                 
61      !! ** Purpose :   Compute the now trend due to the advection of tracers
62      !!      and add it to the general trend of passive tracer equations.
63      !!
64      !! ** Method  :   The advection is evaluated by a second order centered
65      !!      scheme using now fields (leap-frog scheme). In specific areas
66      !!      (vicinity of major river mouths, some straits, or where tn is
67      !!      approaching the freezing point) it is mixed with an upstream
68      !!      scheme for stability reasons.
69      !!         Part 0 : compute the upstream / centered flag
70      !!                  (3D array, zind, defined at T-point (0<zind<1))
71      !!         Part I : horizontal advection
72      !!       * centered flux:
73      !!               zcenu = e2u*e3u  un  mi(ptran)
74      !!               zcenv = e1v*e3v  vn  mj(ptran)
75      !!       * upstream flux:
76      !!               zupsu = e2u*e3u  un  (ptrab(i) or ptrab(i-1) ) [un>0 or <0]
77      !!               zupsv = e1v*e3v  vn  (ptrab(j) or ptrab(j-1) ) [vn>0 or <0]
78      !!       * mixed upstream / centered horizontal advection scheme
79      !!               zcofi = max(zind(i+1), zind(i))
80      !!               zcofj = max(zind(j+1), zind(j))
81      !!               zwx = zcofi * zupsu + (1-zcofi) * zcenu
82      !!               zwy = zcofj * zupsv + (1-zcofj) * zcenv
83      !!       * horizontal advective trend (divergence of the fluxes)
84      !!               ztra = 1/(e1t*e2t*e3t) { di-1[zwx] + dj-1[zwy] }
85      !!       * Add this trend now to the general trend of tracer (ta,sa):
86      !!               ptraa = ptraa + ztra
87      !!       * trend diagnostic ('key_trdtra' defined): the trend is
88      !!      saved for diagnostics. The trends saved is expressed as
89      !!      Uh.gradh(T), i.e.
90      !!                     save trend = ztra + ptran divn
91      !!
92      !!         Part II : vertical advection
93      !!      For temperature (idem for salinity) the advective trend is com-
94      !!      puted as follows :
95      !!            ztra = 1/e3t dk+1[ zwz ]
96      !!      where the vertical advective flux, zwz, is given by :
97      !!            zwz = zcofk * zupst + (1-zcofk) * zcent
98      !!      with
99      !!        zupsv = upstream flux = wn * (ptrab(k) or ptrab(k-1) ) [wn>0 or <0]
100      !!        zcenu = centered flux = wn * mk(tn)
101      !!         The surface boundary condition is :
102      !!      variable volume (lk_vvl = T) : zero advective flux
103      !!      lin. free-surf  (lk_vvl = F) : wn(:,:,1) * ptran(:,:,1)
104      !!         Add this trend now to the general trend of tracer (ta,sa):
105      !!             ptraa = ptraa + ztra
106      !!         Trend diagnostic ('key_trdtra' defined): the trend is
107      !!      saved for diagnostics. The trends saved is expressed as :
108      !!             save trend =  w.gradz(T) = ztra - ptran divn.
109      !!
110      !! ** Action :  - update ptraa  with the now advective tracer trends
111      !!              - save trends if needed
112      !!----------------------------------------------------------------------
113      !!* Module used
114      USE oce         , zwx => ua   ! use ua as workspace
115      USE oce         , zwy => va   ! use va as workspace
116      !!* Arguments
117      INTEGER         , INTENT(in   )                               ::   kt              ! ocean time-step index
118      CHARACTER(len=3), INTENT(in   )                               ::   cdtype          ! =TRA or TRC (tracer indicator)
119      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk)       ::   pun, pvn, pwn   ! 3 ocean velocity components
120      INTEGER         , INTENT(in   )                               ::   kjpt            ! number of tracers
121      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk,kjpt)  ::   ptrab, ptran        ! before and now tracer fields
122      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::   ptraa           ! tracer trend
123      !!* Local declarations
124      INTEGER  ::   ji, jj, jk, jn                   ! dummy loop indices
125      REAL(wp) ::   zbtr, ztra                       ! temporary scalars
126      REAL(wp) ::   zfp_ui, zfp_vj, zfp_w            !    -         -
127      REAL(wp) ::   zfm_ui, zfm_vj, zfm_w            !    -         -
128      REAL(wp) ::   zcofi , zcofj , zcofk            !    -         -
129      REAL(wp) ::   zupsut, zcenut                   !    -         -
130      REAL(wp) ::   zupsvt, zcenvt                   !    -         -
131      REAL(wp) ::   zupst , zcent                    !    -         -
132      REAL(wp) ::   zice                             !    -         -
133      REAL(wp), DIMENSION(jpi,jpj)     ::   ztfreez            ! 2D workspace
134      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwz, zind   ! 3D workspace
135      !!----------------------------------------------------------------------
136
137      IF( kt == nit000 ) THEN
138         IF(lwp) WRITE(numout,*)
139         IF(lwp) WRITE(numout,*) 'tra_adv_cen2 : 2nd order centered advection scheme'
140         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~   Vector optimization case'
141         IF(lwp) WRITE(numout,*)
142         !
143         upsmsk(:,:) = 0.e0                              ! not upstream by default
144         !
145         IF( cp_cfg == "orca" )   CALL ups_orca_set      ! set mixed Upstream/centered scheme near some straits
146         !                                               ! and in closed seas (orca2 and orca4 only)
147         IF( jp_cfg == 2 .AND. .NOT. ln_rstart ) THEN    ! Increase the background in the surface layers
148            avmb(1) = 10.  * avmb(1)      ;      avtb(1) = 10.  * avtb(1)
149            avmb(2) = 10.  * avmb(2)      ;      avtb(2) = 10.  * avtb(2)
150            avmb(3) =  5.  * avmb(3)      ;      avtb(3) =  5.  * avtb(3)
151            avmb(4) =  2.5 * avmb(4)      ;      avtb(4) =  2.5 * avtb(4)
152         ENDIF
153         !
154         l_trd = .FALSE.
155         IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.
156      ENDIF
157      !
158      ! Upstream / centered scheme indicator
159      ! ------------------------------------
160!!gm  not strickly exact : the freezing point should be computed at each ocean levels...
161!!gm  not a big deal since cen2 is no more used in global ice-ocean simulations
162      ztfreez(:,:) = tfreez( tsn(:,:,1,jp_sal) )
163      DO jk = 1, jpk
164         DO jj = 1, jpj
165            DO ji = 1, jpi
166               !                                        ! below ice covered area (if tn < "freezing"+0.1 )
167               IF( tsn(ji,jj,jk,jp_tem) <= ztfreez(ji,jj) + 0.1 ) THEN   ;   zice = 1.e0
168               ELSE                                              ;   zice = 0.e0
169               ENDIF
170               zind(ji,jj,jk) = MAX (   &
171                  rnfmsk(ji,jj) * rnfmsk_z(jk),      &  ! near runoff mouths (& closed sea outflows)
172                  upsmsk(ji,jj)               ,      &  ! some of some straits
173                  zice                               &  ! below ice covered area (if tn < "freezing"+0.1 )
174                  &                  ) * tmask(ji,jj,jk)
175            END DO
176         END DO
177      END DO
178
179      DO jn = 1, kjpt
180         !
181         ! I. Horizontal advection
182         !    ====================
183         !
184         DO jk = 1, jpkm1
185            !                        ! Second order centered tracer flux at u- and v-points
186            DO jj = 1, jpjm1
187               !
188               DO ji = 1, fs_jpim1   ! vector opt.
189                  ! upstream indicator
190                  zcofi = MAX( zind(ji+1,jj,jk), zind(ji,jj,jk) )
191                  zcofj = MAX( zind(ji,jj+1,jk), zind(ji,jj,jk) )
192                  !
193                  ! upstream scheme
194                  zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) )
195                  zfm_ui = pun(ji,jj,jk) - ABS( pun(ji,jj,jk) )
196                  zfp_vj = pvn(ji,jj,jk) + ABS( pvn(ji,jj,jk) )
197                  zfm_vj = pvn(ji,jj,jk) - ABS( pvn(ji,jj,jk) )
198                  zupsut = zfp_ui * ptrab(ji,jj,jk,jn) + zfm_ui * ptrab(ji+1,jj  ,jk,jn)
199                  zupsvt = zfp_vj * ptrab(ji,jj,jk,jn) + zfm_vj * ptrab(ji  ,jj+1,jk,jn)
200                  ! centered scheme
201                  zcenut = pun(ji,jj,jk) * ( ptran(ji,jj,jk,jn) + ptran(ji+1,jj  ,jk,jn) )
202                  zcenvt = pvn(ji,jj,jk) * ( ptran(ji,jj,jk,jn) + ptran(ji  ,jj+1,jk,jn) )
203                  ! mixed centered / upstream scheme
204                  zwx(ji,jj,jk) = 0.5 * ( zcofi * zupsut + (1.-zcofi) * zcenut )
205                  zwy(ji,jj,jk) = 0.5 * ( zcofj * zupsvt + (1.-zcofj) * zcenvt )
206               END DO
207            END DO
208         END DO
209
210         ! II. Vertical advection
211         !     ==================
212         !
213         !                                                ! Vertical advective fluxes
214         zwz(:,:,jpk) = 0.e0                                   ! Bottom  value : flux set to zero
215         !                                                     ! Surface value :
216         IF( lk_vvl ) THEN   ;   zwz(:,:, 1 ) = 0.e0                         ! volume variable
217         ELSE                ;   zwz(:,:, 1 ) = pwn(:,:,1) * ptran(:,:,1,jn)   ! linear free surface
218         ENDIF
219         !
220         DO jk = 2, jpk              ! Second order centered tracer flux at w-point
221            DO jj = 2, jpjm1
222               DO ji = fs_2, fs_jpim1   ! vector opt.
223                  ! upstream indicator
224                  zcofk = MAX( zind(ji,jj,jk-1), zind(ji,jj,jk) ) 
225                  ! mixed centered / upstream scheme
226                  zfp_w = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) )
227                  zfm_w = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) )
228                  zupst = zfp_w * ptrab(ji,jj,jk,jn) + zfm_w * ptrab(ji,jj,jk-1,jn)
229                  ! centered scheme
230                  zcent = pwn(ji,jj,jk) * ( ptran(ji,jj,jk,jn) + ptran(ji,jj,jk-1,jn) )
231                  ! mixed centered / upstream scheme
232                  zwz(ji,jj,jk) = 0.5 * ( zcofk * zupst + (1.-zcofk) * zcent )
233               END DO
234            END DO
235         END DO
236
237         ! II. Divergence of advective fluxes
238         ! ----------------------------------
239         DO jk = 1, jpkm1
240            DO jj = 2, jpjm1
241               DO ji = fs_2, fs_jpim1   ! vector opt.
242                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) *  fse3t(ji,jj,jk) )
243                  ! advective trends
244                  ztra = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
245                  &                + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   &
246                  &                + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1)  )
247                  ! advective trends added to the general tracer trends
248                  ptraa(ji,jj,jk,jn) = ptraa(ji,jj,jk,jn) + ztra
249               END DO
250            END DO
251         END DO
252
253         !                                 ! trend diagnostics (contribution of upstream fluxes)
254         IF( l_trd ) THEN
255            CALL trd_tra( kt, cdtype, jn, jptra_trd_xad, zwx, pun, ptran(:,:,:,jn) )
256            CALL trd_tra( kt, cdtype, jn, jptra_trd_yad, zwy, pvn, ptran(:,:,:,jn) )
257            CALL trd_tra( kt, cdtype, jn, jptra_trd_zad, zwz, pwn, ptran(:,:,:,jn) )
258         END IF
259         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
260         IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN 
261           IF( jn == jp_tem )  pht_adv(:) = ptr_vj( zwy(:,:,:) )
262           IF( jn == jp_sal )  pst_adv(:) = ptr_vj( zwy(:,:,:) )
263         ENDIF
264         !
265      ENDDO
266
267      ! ---------------------------  required in restart file to ensure restartability)
268      ! avmb, avtb will be read in zdfini in restart case as they are used in zdftke, kpp etc...
269      IF( lrst_oce .AND. cdtype == 'TRA' ) THEN
270         CALL iom_rstput( kt, nitrst, numrow, 'avmb', avmb )
271         CALL iom_rstput( kt, nitrst, numrow, 'avtb', avtb )
272      ENDIF
273      !
274   END SUBROUTINE tra_adv_cen2
275   
276   
277   SUBROUTINE ups_orca_set
278      !!----------------------------------------------------------------------
279      !!                  ***  ROUTINE ups_orca_set  ***
280      !!       
281      !! ** Purpose :   add a portion of upstream scheme in area where the
282      !!                centered scheme generates too strong overshoot
283      !!
284      !! ** Method  :   orca (R4 and R2) confiiguration setting. Set upsmsk
285      !!                array to nozero value in some straith.
286      !!
287      !! ** Action : - upsmsk set to 1 at some strait, 0 elsewhere for orca
288      !!----------------------------------------------------------------------
289      INTEGER  ::   ii0, ii1, ij0, ij1      ! temporary integers
290      !!----------------------------------------------------------------------
291     
292      ! mixed upstream/centered scheme near river mouths
293      ! ------------------------------------------------
294      SELECT CASE ( jp_cfg )
295      !                                        ! =======================
296      CASE ( 4 )                               !  ORCA_R4 configuration
297         !                                     ! =======================
298         !                                          ! Gibraltar Strait
299         ii0 =  70   ;   ii1 =  71
300         ij0 =  52   ;   ij1 =  53   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50
301         !
302         !                                     ! =======================
303      CASE ( 2 )                               !  ORCA_R2 configuration
304         !                                     ! =======================
305         !                                          ! Gibraltar Strait
306         ij0 = 102   ;   ij1 = 102
307         ii0 = 138   ;   ii1 = 138   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.20
308         ii0 = 139   ;   ii1 = 139   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.40
309         ii0 = 140   ;   ii1 = 140   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50
310         ij0 = 101   ;   ij1 = 102
311         ii0 = 141   ;   ii1 = 141   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50
312         !                                          ! Bab el Mandeb Strait
313         ij0 =  87   ;   ij1 =  88
314         ii0 = 164   ;   ii1 = 164   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.10
315         ij0 =  88   ;   ij1 =  88
316         ii0 = 163   ;   ii1 = 163   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.25
317         ii0 = 162   ;   ii1 = 162   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.40
318         ii0 = 160   ;   ii1 = 161   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50
319         ij0 =  89   ;   ij1 =  89
320         ii0 = 158   ;   ii1 = 160   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.25
321         ij0 =  90   ;   ij1 =  90
322         ii0 = 160   ;   ii1 = 160   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.25
323         !                                          ! Sound Strait
324         ij0 = 116   ;   ij1 = 116
325         ii0 = 144   ;   ii1 = 144   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.25
326         ii0 = 145   ;   ii1 = 147   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50
327         ii0 = 148   ;   ii1 = 148   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.25
328         !
329      END SELECT 
330     
331      ! mixed upstream/centered scheme over closed seas
332      ! -----------------------------------------------
333      CALL clo_ups( upsmsk(:,:) )
334      !
335   END SUBROUTINE ups_orca_set
336
337   !!======================================================================
338END MODULE traadv_cen2
Note: See TracBrowser for help on using the repository browser.