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

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

Improve the merge TRA-TRC, see ticket:701
  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 18.0 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: tn, sn  ! 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.3 , LOCEAN-IPSL (2010)
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      &                                 ptb, ptn, pta, 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(ptn)
74      !!               zcenv = e1v*e3v  vn  mj(ptn)
75      !!       * upstream flux:
76      !!               zupsu = e2u*e3u  un  (ptb(i) or ptb(i-1) ) [un>0 or <0]
77      !!               zupsv = e1v*e3v  vn  (ptb(j) or ptb(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      !!               pta = pta + 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 + ptn 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 * (ptb(k) or ptb(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) * ptn(:,:,1)
104      !!         Add this trend now to the general trend of tracer (ta,sa):
105      !!             pta = pta + 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 - ptn divn.
109      !!
110      !! ** Action :  - update pta  with the now advective tracer trends
111      !!              - save trends if needed
112      !!----------------------------------------------------------------------
113      USE oce         , zwx => ua   ! use ua as workspace
114      USE oce         , zwy => va   ! use va as workspace
115      !!
116      INTEGER         , INTENT(in   )                               ::   kt              ! ocean time-step index
117      CHARACTER(len=3), INTENT(in   )                               ::   cdtype          ! =TRA or TRC (tracer indicator)
118      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk)       ::   pun, pvn, pwn   ! 3 ocean velocity components
119      INTEGER         , INTENT(in   )                               ::   kjpt            ! number of tracers
120      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk,kjpt)  ::   ptb, ptn        ! before and now tracer fields
121      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::   pta           ! tracer trend
122      !!
123      INTEGER  ::   ji, jj, jk, jn                   ! dummy loop indices
124      REAL(wp) ::   zbtr, ztra                       ! temporary scalars
125      REAL(wp) ::   zfp_ui, zfp_vj, zfp_w            !    -         -
126      REAL(wp) ::   zfm_ui, zfm_vj, zfm_w            !    -         -
127      REAL(wp) ::   zcofi , zcofj , zcofk            !    -         -
128      REAL(wp) ::   zupsut, zcenut                   !    -         -
129      REAL(wp) ::   zupsvt, zcenvt                   !    -         -
130      REAL(wp) ::   zupst , zcent                    !    -         -
131      REAL(wp) ::   zice                             !    -         -
132      REAL(wp), DIMENSION(jpi,jpj)     ::   ztfreez            ! 2D workspace
133      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwz, zind   ! 3D workspace
134      !!----------------------------------------------------------------------
135
136      IF( kt == nit000 ) THEN
137         IF(lwp) WRITE(numout,*)
138         IF(lwp) WRITE(numout,*) 'tra_adv_cen2 : 2nd order centered advection scheme'
139         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~   Vector optimization case'
140         IF(lwp) WRITE(numout,*)
141         !
142         upsmsk(:,:) = 0.e0                              ! not upstream by default
143         !
144         IF( cp_cfg == "orca" )   CALL ups_orca_set      ! set mixed Upstream/centered scheme near some straits
145         !                                               ! and in closed seas (orca2 and orca4 only)
146         IF( jp_cfg == 2 .AND. .NOT. ln_rstart ) THEN    ! Increase the background in the surface layers
147            avmb(1) = 10.  * avmb(1)      ;      avtb(1) = 10.  * avtb(1)
148            avmb(2) = 10.  * avmb(2)      ;      avtb(2) = 10.  * avtb(2)
149            avmb(3) =  5.  * avmb(3)      ;      avtb(3) =  5.  * avtb(3)
150            avmb(4) =  2.5 * avmb(4)      ;      avtb(4) =  2.5 * avtb(4)
151         ENDIF
152         !
153         l_trd = .FALSE.
154         IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.
155      ENDIF
156      !
157      ! Upstream / centered scheme indicator
158      ! ------------------------------------
159!!gm  not strickly exact : the freezing point should be computed at each ocean levels...
160!!gm  not a big deal since cen2 is no more used in global ice-ocean simulations
161      ztfreez(:,:) = tfreez( sn(:,:,1) )
162      DO jk = 1, jpk
163         DO jj = 1, jpj
164            DO ji = 1, jpi
165               !                                        ! below ice covered area (if tn < "freezing"+0.1 )
166               IF( tn(ji,jj,jk) <= ztfreez(ji,jj) + 0.1 ) THEN   ;   zice = 1.e0
167               ELSE                                              ;   zice = 0.e0
168               ENDIF
169               zind(ji,jj,jk) = MAX (   &
170                  rnfmsk(ji,jj) * rnfmsk_z(jk),      &  ! near runoff mouths (& closed sea outflows)
171                  upsmsk(ji,jj)               ,      &  ! some of some straits
172                  zice                               &  ! below ice covered area (if tn < "freezing"+0.1 )
173                  &                  ) * tmask(ji,jj,jk)
174            END DO
175         END DO
176      END DO
177
178      DO jn = 1, kjpt
179         !
180         ! I. Horizontal advection
181         !    ====================
182         !
183         DO jk = 1, jpkm1
184            !                        ! Second order centered tracer flux at u- and v-points
185            DO jj = 1, jpjm1
186               !
187               DO ji = 1, fs_jpim1   ! vector opt.
188                  ! upstream indicator
189                  zcofi = MAX( zind(ji+1,jj,jk), zind(ji,jj,jk) )
190                  zcofj = MAX( zind(ji,jj+1,jk), zind(ji,jj,jk) )
191                  !
192                  ! upstream scheme
193                  zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) )
194                  zfm_ui = pun(ji,jj,jk) - ABS( pun(ji,jj,jk) )
195                  zfp_vj = pvn(ji,jj,jk) + ABS( pvn(ji,jj,jk) )
196                  zfm_vj = pvn(ji,jj,jk) - ABS( pvn(ji,jj,jk) )
197                  zupsut = zfp_ui * ptb(ji,jj,jk,jn) + zfm_ui * ptb(ji+1,jj  ,jk,jn)
198                  zupsvt = zfp_vj * ptb(ji,jj,jk,jn) + zfm_vj * ptb(ji  ,jj+1,jk,jn)
199                  ! centered scheme
200                  zcenut = pun(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj  ,jk,jn) )
201                  zcenvt = pvn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji  ,jj+1,jk,jn) )
202                  ! mixed centered / upstream scheme
203                  zwx(ji,jj,jk) = 0.5 * ( zcofi * zupsut + (1.-zcofi) * zcenut )
204                  zwy(ji,jj,jk) = 0.5 * ( zcofj * zupsvt + (1.-zcofj) * zcenvt )
205               END DO
206            END DO
207         END DO
208
209         ! II. Vertical advection
210         !     ==================
211         !
212         !                                                ! Vertical advective fluxes
213         zwz(:,:,jpk) = 0.e0                                   ! Bottom  value : flux set to zero
214         !                                                     ! Surface value :
215         IF( lk_vvl ) THEN   ;   zwz(:,:, 1 ) = 0.e0                         ! volume variable
216         ELSE                ;   zwz(:,:, 1 ) = pwn(:,:,1) * ptn(:,:,1,jn)   ! linear free surface
217         ENDIF
218         !
219         DO jk = 2, jpk              ! Second order centered tracer flux at w-point
220            DO jj = 2, jpjm1
221               DO ji = fs_2, fs_jpim1   ! vector opt.
222                  ! upstream indicator
223                  zcofk = MAX( zind(ji,jj,jk-1), zind(ji,jj,jk) ) 
224                  ! mixed centered / upstream scheme
225                  zfp_w = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) )
226                  zfm_w = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) )
227                  zupst = zfp_w * ptb(ji,jj,jk,jn) + zfm_w * ptb(ji,jj,jk-1,jn)
228                  ! centered scheme
229                  zcent = pwn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj,jk-1,jn) )
230                  ! mixed centered / upstream scheme
231                  zwz(ji,jj,jk) = 0.5 * ( zcofk * zupst + (1.-zcofk) * zcent )
232               END DO
233            END DO
234         END DO
235
236         ! II. Divergence of advective fluxes
237         ! ----------------------------------
238         DO jk = 1, jpkm1
239            DO jj = 2, jpjm1
240               DO ji = fs_2, fs_jpim1   ! vector opt.
241                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) *  fse3t(ji,jj,jk) )
242                  ! advective trends
243                  ztra = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
244                  &                + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   &
245                  &                + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1)  )
246                  ! advective trends added to the general tracer trends
247                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
248               END DO
249            END DO
250         END DO
251
252         !                                 ! trend diagnostics (contribution of upstream fluxes)
253         IF( l_trd ) THEN
254            CALL trd_tra( kt, cdtype, jn, jptra_trd_xad, zwx, pun, ptn(:,:,:,jn) )
255            CALL trd_tra( kt, cdtype, jn, jptra_trd_yad, zwy, pvn, ptn(:,:,:,jn) )
256            CALL trd_tra( kt, cdtype, jn, jptra_trd_zad, zwz, pwn, ptn(:,:,:,jn) )
257         END IF
258         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
259         IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN 
260           IF( jn == jp_tem )  pht_adv(:) = ptr_vj( zwy(:,:,:) )
261           IF( jn == jp_sal )  pst_adv(:) = ptr_vj( zwy(:,:,:) )
262         ENDIF
263         !
264      ENDDO
265
266      ! ---------------------------  required in restart file to ensure restartability)
267      ! avmb, avtb will be read in zdfini in restart case as they are used in zdftke, kpp etc...
268      IF( lrst_oce .AND. cdtype == 'TRA' ) THEN
269         CALL iom_rstput( kt, nitrst, numrow, 'avmb', avmb )
270         CALL iom_rstput( kt, nitrst, numrow, 'avtb', avtb )
271      ENDIF
272      !
273   END SUBROUTINE tra_adv_cen2
274   
275   
276   SUBROUTINE ups_orca_set
277      !!----------------------------------------------------------------------
278      !!                  ***  ROUTINE ups_orca_set  ***
279      !!       
280      !! ** Purpose :   add a portion of upstream scheme in area where the
281      !!                centered scheme generates too strong overshoot
282      !!
283      !! ** Method  :   orca (R4 and R2) confiiguration setting. Set upsmsk
284      !!                array to nozero value in some straith.
285      !!
286      !! ** Action : - upsmsk set to 1 at some strait, 0 elsewhere for orca
287      !!----------------------------------------------------------------------
288      INTEGER  ::   ii0, ii1, ij0, ij1      ! temporary integers
289      !!----------------------------------------------------------------------
290     
291      ! mixed upstream/centered scheme near river mouths
292      ! ------------------------------------------------
293      SELECT CASE ( jp_cfg )
294      !                                        ! =======================
295      CASE ( 4 )                               !  ORCA_R4 configuration
296         !                                     ! =======================
297         !                                          ! Gibraltar Strait
298         ii0 =  70   ;   ii1 =  71
299         ij0 =  52   ;   ij1 =  53   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50
300         !
301         !                                     ! =======================
302      CASE ( 2 )                               !  ORCA_R2 configuration
303         !                                     ! =======================
304         !                                          ! Gibraltar Strait
305         ij0 = 102   ;   ij1 = 102
306         ii0 = 138   ;   ii1 = 138   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.20
307         ii0 = 139   ;   ii1 = 139   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.40
308         ii0 = 140   ;   ii1 = 140   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50
309         ij0 = 101   ;   ij1 = 102
310         ii0 = 141   ;   ii1 = 141   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50
311         !                                          ! Bab el Mandeb Strait
312         ij0 =  87   ;   ij1 =  88
313         ii0 = 164   ;   ii1 = 164   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.10
314         ij0 =  88   ;   ij1 =  88
315         ii0 = 163   ;   ii1 = 163   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.25
316         ii0 = 162   ;   ii1 = 162   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.40
317         ii0 = 160   ;   ii1 = 161   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50
318         ij0 =  89   ;   ij1 =  89
319         ii0 = 158   ;   ii1 = 160   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.25
320         ij0 =  90   ;   ij1 =  90
321         ii0 = 160   ;   ii1 = 160   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.25
322         !                                          ! Sound Strait
323         ij0 = 116   ;   ij1 = 116
324         ii0 = 144   ;   ii1 = 144   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.25
325         ii0 = 145   ;   ii1 = 147   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50
326         ii0 = 148   ;   ii1 = 148   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.25
327         !
328      END SELECT 
329     
330      ! mixed upstream/centered scheme over closed seas
331      ! -----------------------------------------------
332      CALL clo_ups( upsmsk(:,:) )
333      !
334   END SUBROUTINE ups_orca_set
335
336   !!======================================================================
337END MODULE traadv_cen2
Note: See TracBrowser for help on using the repository browser.