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/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_cen2.F90 @ 2399

Last change on this file since 2399 was 2399, checked in by gm, 13 years ago

v3.3beta: diaptr (poleward heat & salt transports) #759 : rewriting including dynamical allocation + DOCTOR names

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