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

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

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