source: branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_cen2.F90 @ 3876

Last change on this file since 3876 was 3876, checked in by gm, 8 years ago

dev_r3858_CNRS3_Ediag: #927 phasing with 2011/dev_r3309_LOCEAN12_Ediag branche + mxl diag update

  • Property svn:keywords set to Id
File size: 18.6 KB
Line 
1MODULE traadv_cen2
2   !!======================================================================
3   !!                     ***  MODULE  traadv_cen2  ***
4   !! Ocean  tracers:  horizontal & vertical advective trend
5   !!======================================================================
6   !! History :  OPA  ! 2001-08  (G. Madec, E. Durand) v8.2 trahad+trazad=traadv
7   !!   NEMO     1.0  ! 2002-06  (G. Madec)  F90: Free form and module
8   !!             -   ! 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 trd_oce         ! trends: ocean variables
24   USE trdtra          ! trends manager: tracers
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 trc_oce         ! share passive tracers/Ocean variables
32   USE lib_mpp         ! MPP library
33   USE wrk_nemo        ! Memory Allocation
34   USE timing          ! Timing
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, ALLOCATABLE, SAVE, DIMENSION(:,:) :: 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 , NEMO Consortium (2010)
51   !! $Id$
52   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
53   !!----------------------------------------------------------------------
54CONTAINS
55
56   SUBROUTINE tra_adv_cen2( kt, kit000, 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 (l_trdtra=T or l_trctra=T): the trend is
88      !!      saved for diagnostics. The trends saved is expressed as
89      !!      Uh.gradh(T), i.e.  save trend = ztra + ptn divn
90      !!
91      !!         Part II : vertical advection
92      !!      For temperature (idem for salinity) the advective trend is com-
93      !!      puted as follows :
94      !!            ztra = 1/e3t dk+1[ zwz ]
95      !!      where the vertical advective flux, zwz, is given by :
96      !!            zwz = zcofk * zupst + (1-zcofk) * zcent
97      !!      with
98      !!        zupsv = upstream flux = wn * (ptb(k) or ptb(k-1) ) [wn>0 or <0]
99      !!        zcenu = centered flux = wn * mk(tn)
100      !!         The surface boundary condition is :
101      !!      variable volume (lk_vvl = T) : zero advective flux
102      !!      lin. free-surf  (lk_vvl = F) : wn(:,:,1) * ptn(:,:,1)
103      !!         Add this trend now to the general trend of tracer (ta,sa):
104      !!             pta = pta + ztra
105      !!         Trend diagnostic (l_trdtra=T or l_trctra=T): the trend is
106      !!      saved for diagnostics. The trends saved is expressed as :
107      !!             save trend =  w.gradz(T) = ztra - ptn divn.
108      !!
109      !! ** Action :  - update pta  with the now advective tracer trends
110      !!              - save trends if needed
111      !!----------------------------------------------------------------------
112      USE oce     , ONLY:   zwx => ua        , zwy  => va          ! (ua,va) used as 3D workspace
113      !
114      INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index
115      INTEGER                              , INTENT(in   ) ::   kit000          ! first time step index
116      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator)
117      INTEGER                              , INTENT(in   ) ::   kjpt            ! number of tracers
118      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pun, pvn, pwn   ! 3 ocean velocity components
119      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb, ptn        ! before and now tracer fields
120      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta             ! tracer trend
121      !
122      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
123      INTEGER  ::   ierr             ! local integer
124      REAL(wp) ::   zbtr, ztra                            ! local scalars
125      REAL(wp) ::   zfp_ui, zfp_vj, zfp_w, zcofi          !   -      -
126      REAL(wp) ::   zfm_ui, zfm_vj, zfm_w, zcofj, zcofk   !   -      -
127      REAL(wp) ::   zupsut, zcenut, zupst                 !   -      -
128      REAL(wp) ::   zupsvt, zcenvt, zcent, zice           !   -      -
129      REAL(wp), POINTER, DIMENSION(:,:  ) :: ztfreez 
130      REAL(wp), POINTER, DIMENSION(:,:,:) :: zwz, zind
131      !!----------------------------------------------------------------------
132      !
133      IF( nn_timing == 1 )  CALL timing_start('tra_adv_cen2')
134      !
135      CALL wrk_alloc( jpi, jpj, ztfreez )
136      CALL wrk_alloc( jpi, jpj, jpk, zwz, zind )
137      !
138
139      IF( kt == kit000 )  THEN
140         IF(lwp) WRITE(numout,*)
141         IF(lwp) WRITE(numout,*) 'tra_adv_cen2 : 2nd order centered advection scheme on ', cdtype
142         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~ '
143         IF(lwp) WRITE(numout,*)
144         !
145         IF( .NOT. ALLOCATED( upsmsk ) )  THEN
146             ALLOCATE( upsmsk(jpi,jpj), STAT=ierr )
147             IF( ierr /= 0 )   CALL ctl_stop('STOP', 'tra_adv_cen2: unable to allocate array')
148         ENDIF
149
150         !
151         upsmsk(:,:) = 0._wp                             ! not upstream by default
152         !
153         IF( cp_cfg == "orca" )   CALL ups_orca_set      ! set mixed Upstream/centered scheme near some straits
154         !                                               ! and in closed seas (orca2 and orca4 only)
155         IF( jp_cfg == 2 .AND. .NOT. ln_rstart ) THEN    ! Increase the background in the surface layers
156            avmb(1) = 10.  * avmb(1)      ;      avtb(1) = 10.  * avtb(1)
157            avmb(2) = 10.  * avmb(2)      ;      avtb(2) = 10.  * avtb(2)
158            avmb(3) =  5.  * avmb(3)      ;      avtb(3) =  5.  * avtb(3)
159            avmb(4) =  2.5 * avmb(4)      ;      avtb(4) =  2.5 * avtb(4)
160         ENDIF
161         !
162         l_trd = .FALSE.
163         IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )   l_trd = .TRUE.
164      ENDIF
165      !
166      ! Upstream / centered scheme indicator
167      ! ------------------------------------
168!!gm  not strickly exact : the freezing point should be computed at each ocean levels...
169!!gm  not a big deal since cen2 is no more used in global ice-ocean simulations
170      ztfreez(:,:) = tfreez( tsn(:,:,1,jp_sal) )
171      DO jk = 1, jpk
172         DO jj = 1, jpj
173            DO ji = 1, jpi
174               !                                        ! below ice covered area (if tn < "freezing"+0.1 )
175               IF( tsn(ji,jj,jk,jp_tem) <= ztfreez(ji,jj) + 0.1 ) THEN   ;   zice = 1.e0
176               ELSE                                                      ;   zice = 0.e0
177               ENDIF
178               zind(ji,jj,jk) = MAX (   &
179                  rnfmsk(ji,jj) * rnfmsk_z(jk),      &  ! near runoff mouths (& closed sea outflows)
180                  upsmsk(ji,jj)               ,      &  ! some of some straits
181                  zice                               &  ! below ice covered area (if tn < "freezing"+0.1 )
182                  &                  ) * tmask(ji,jj,jk)
183            END DO
184         END DO
185      END DO
186
187      DO jn = 1, kjpt
188         !
189         ! I. Horizontal advection
190         !    ====================
191         !
192         DO jk = 1, jpkm1
193            !                        ! Second order centered tracer flux at u- and v-points
194            DO jj = 1, jpjm1
195               !
196               DO ji = 1, fs_jpim1   ! vector opt.
197                  ! upstream indicator
198                  zcofi = MAX( zind(ji+1,jj,jk), zind(ji,jj,jk) )
199                  zcofj = MAX( zind(ji,jj+1,jk), zind(ji,jj,jk) )
200                  !
201                  ! upstream scheme
202                  zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) )
203                  zfm_ui = pun(ji,jj,jk) - ABS( pun(ji,jj,jk) )
204                  zfp_vj = pvn(ji,jj,jk) + ABS( pvn(ji,jj,jk) )
205                  zfm_vj = pvn(ji,jj,jk) - ABS( pvn(ji,jj,jk) )
206                  zupsut = zfp_ui * ptb(ji,jj,jk,jn) + zfm_ui * ptb(ji+1,jj  ,jk,jn)
207                  zupsvt = zfp_vj * ptb(ji,jj,jk,jn) + zfm_vj * ptb(ji  ,jj+1,jk,jn)
208                  ! centered scheme
209                  zcenut = pun(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj  ,jk,jn) )
210                  zcenvt = pvn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji  ,jj+1,jk,jn) )
211                  ! mixed centered / upstream scheme
212                  zwx(ji,jj,jk) = 0.5 * ( zcofi * zupsut + (1.-zcofi) * zcenut )
213                  zwy(ji,jj,jk) = 0.5 * ( zcofj * zupsvt + (1.-zcofj) * zcenvt )
214               END DO
215            END DO
216         END DO
217
218         ! II. Vertical advection
219         !     ==================
220         !
221         !                                                ! Vertical advective fluxes
222         zwz(:,:,jpk) = 0.e0                                   ! Bottom  value : flux set to zero
223         !                                                     ! Surface value :
224         IF( lk_vvl ) THEN   ;   zwz(:,:, 1 ) = 0.e0                         ! volume variable
225         ELSE                ;   zwz(:,:, 1 ) = pwn(:,:,1) * ptn(:,:,1,jn)   ! linear free surface
226         ENDIF
227         !
228         DO jk = 2, jpk              ! Second order centered tracer flux at w-point
229            DO jj = 2, jpjm1
230               DO ji = fs_2, fs_jpim1   ! vector opt.
231                  ! upstream indicator
232                  zcofk = MAX( zind(ji,jj,jk-1), zind(ji,jj,jk) ) 
233                  ! mixed centered / upstream scheme
234                  zfp_w = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) )
235                  zfm_w = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) )
236                  zupst = zfp_w * ptb(ji,jj,jk,jn) + zfm_w * ptb(ji,jj,jk-1,jn)
237                  ! centered scheme
238                  zcent = pwn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj,jk-1,jn) )
239                  ! mixed centered / upstream scheme
240                  zwz(ji,jj,jk) = 0.5 * ( zcofk * zupst + (1.-zcofk) * zcent )
241               END DO
242            END DO
243         END DO
244
245         ! II. Divergence of advective fluxes
246         ! ----------------------------------
247         DO jk = 1, jpkm1
248            DO jj = 2, jpjm1
249               DO ji = fs_2, fs_jpim1   ! vector opt.
250                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) *  fse3t(ji,jj,jk) )
251                  ! advective trends
252                  ztra = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
253                  &                + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   &
254                  &                + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1)  )
255                  ! advective trends added to the general tracer trends
256                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
257               END DO
258            END DO
259         END DO
260
261         !                                 ! trend diagnostics
262         IF( l_trd ) THEN
263            CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptn(:,:,:,jn) )
264            CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptn(:,:,:,jn) )
265            CALL trd_tra( kt, cdtype, jn, jptra_zad, zwz, pwn, ptn(:,:,:,jn) )
266         END IF
267         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
268         IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN 
269           IF( jn == jp_tem )  htr_adv(:) = ptr_vj( zwy(:,:,:) )
270           IF( jn == jp_sal )  str_adv(:) = ptr_vj( zwy(:,:,:) )
271         ENDIF
272         !
273      END DO
274
275      ! ---------------------------  required in restart file to ensure restartability)
276      ! avmb, avtb will be read in zdfini in restart case as they are used in zdftke, kpp etc...
277      IF( lrst_oce .AND. cdtype == 'TRA' ) THEN
278         CALL iom_rstput( kt, nitrst, numrow, 'avmb', avmb )
279         CALL iom_rstput( kt, nitrst, numrow, 'avtb', avtb )
280      ENDIF
281      !
282      CALL wrk_dealloc( jpi, jpj, ztfreez )
283      CALL wrk_dealloc( jpi, jpj, jpk, zwz, zind )
284      !
285      IF( nn_timing == 1 )  CALL timing_stop('tra_adv_cen2')
286      !
287   END SUBROUTINE tra_adv_cen2
288   
289   
290   SUBROUTINE ups_orca_set
291      !!----------------------------------------------------------------------
292      !!                  ***  ROUTINE ups_orca_set  ***
293      !!       
294      !! ** Purpose :   add a portion of upstream scheme in area where the
295      !!                centered scheme generates too strong overshoot
296      !!
297      !! ** Method  :   orca (R4 and R2) confiiguration setting. Set upsmsk
298      !!                array to nozero value in some straith.
299      !!
300      !! ** Action : - upsmsk set to 1 at some strait, 0 elsewhere for orca
301      !!----------------------------------------------------------------------
302      INTEGER  ::   ii0, ii1, ij0, ij1      ! temporary integers
303      !!----------------------------------------------------------------------
304     
305      !
306      IF( nn_timing == 1 )  CALL timing_start('ups_orca_set')
307      !
308      ! mixed upstream/centered scheme near river mouths
309      ! ------------------------------------------------
310      SELECT CASE ( jp_cfg )
311      !                                        ! =======================
312      CASE ( 4 )                               !  ORCA_R4 configuration
313         !                                     ! =======================
314         !                                          ! Gibraltar Strait
315         ii0 =  70   ;   ii1 =  71
316         ij0 =  52   ;   ij1 =  53   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50
317         !
318         !                                     ! =======================
319      CASE ( 2 )                               !  ORCA_R2 configuration
320         !                                     ! =======================
321         !                                          ! Gibraltar Strait
322         ij0 = 102   ;   ij1 = 102
323         ii0 = 138   ;   ii1 = 138   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.20
324         ii0 = 139   ;   ii1 = 139   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.40
325         ii0 = 140   ;   ii1 = 140   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50
326         ij0 = 101   ;   ij1 = 102
327         ii0 = 141   ;   ii1 = 141   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50
328         !                                          ! Bab el Mandeb Strait
329         ij0 =  87   ;   ij1 =  88
330         ii0 = 164   ;   ii1 = 164   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.10
331         ij0 =  88   ;   ij1 =  88
332         ii0 = 163   ;   ii1 = 163   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.25
333         ii0 = 162   ;   ii1 = 162   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.40
334         ii0 = 160   ;   ii1 = 161   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50
335         ij0 =  89   ;   ij1 =  89
336         ii0 = 158   ;   ii1 = 160   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.25
337         ij0 =  90   ;   ij1 =  90
338         ii0 = 160   ;   ii1 = 160   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.25
339         !                                          ! Sound Strait
340         ij0 = 116   ;   ij1 = 116
341         ii0 = 144   ;   ii1 = 144   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.25
342         ii0 = 145   ;   ii1 = 147   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50
343         ii0 = 148   ;   ii1 = 148   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.25
344         !
345      END SELECT 
346     
347      ! mixed upstream/centered scheme over closed seas
348      ! -----------------------------------------------
349      CALL clo_ups( upsmsk(:,:) )
350      !
351      IF( nn_timing == 1 )  CALL timing_stop('ups_orca_set')
352      !
353   END SUBROUTINE ups_orca_set
354
355   !!======================================================================
356END MODULE traadv_cen2
Note: See TracBrowser for help on using the repository browser.