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

source: branches/2014/dev_CNRS1_10_TEOS10_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_cen2.F90 @ 4915

Last change on this file since 4915 was 4619, checked in by gm, 10 years ago

#1294 : TEOS-10 and Ediag

  • Property svn:keywords set to Id
File size: 18.4 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 traadv.F90
40
41   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: upsmsk !: mixed upstream/centered scheme near some straits
42   !                                                             !  and in closed seas (orca 2 and 4 configurations)
43   !! * Substitutions
44#  include "domzgr_substitute.h90"
45#  include "vectopt_loop_substitute.h90"
46   !!----------------------------------------------------------------------
47   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
48   !! $Id$
49   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
50   !!----------------------------------------------------------------------
51CONTAINS
52
53   SUBROUTINE tra_adv_cen2( kt, kit000, cdtype, pun, pvn, pwn,     &
54      &                                         ptb, ptn, pta, kjpt   ) 
55      !!----------------------------------------------------------------------
56      !!                  ***  ROUTINE tra_adv_cen2  ***
57      !!                 
58      !! ** Purpose :   Compute the now trend due to the advection of tracers
59      !!      and add it to the general trend of passive tracer equations.
60      !!
61      !! ** Method  :   The advection is evaluated by a second order centered
62      !!      scheme using now fields (leap-frog scheme). In specific areas
63      !!      (vicinity of major river mouths, some straits, or where tn is
64      !!      approaching the freezing point) it is mixed with an upstream
65      !!      scheme for stability reasons.
66      !!         Part 0 : compute the upstream / centered flag
67      !!                  (3D array, zind, defined at T-point (0<zind<1))
68      !!         Part I : horizontal advection
69      !!       * centered flux:
70      !!               zcenu = e2u*e3u  un  mi(ptn)
71      !!               zcenv = e1v*e3v  vn  mj(ptn)
72      !!       * upstream flux:
73      !!               zupsu = e2u*e3u  un  (ptb(i) or ptb(i-1) ) [un>0 or <0]
74      !!               zupsv = e1v*e3v  vn  (ptb(j) or ptb(j-1) ) [vn>0 or <0]
75      !!       * mixed upstream / centered horizontal advection scheme
76      !!               zcofi = max(zind(i+1), zind(i))
77      !!               zcofj = max(zind(j+1), zind(j))
78      !!               zwx = zcofi * zupsu + (1-zcofi) * zcenu
79      !!               zwy = zcofj * zupsv + (1-zcofj) * zcenv
80      !!       * horizontal advective trend (divergence of the fluxes)
81      !!               ztra = 1/(e1t*e2t*e3t) { di-1[zwx] + dj-1[zwy] }
82      !!       * Add this trend now to the general trend of tracer (ta,sa):
83      !!               pta = pta + ztra
84      !!       * trend diagnostic (l_trdtra=T or l_trctra=T): the trend is
85      !!      saved for diagnostics. The trends saved is expressed as
86      !!      Uh.gradh(T), i.e.  save trend = ztra + ptn divn
87      !!
88      !!         Part II : vertical advection
89      !!      For temperature (idem for salinity) the advective trend is com-
90      !!      puted as follows :
91      !!            ztra = 1/e3t dk+1[ zwz ]
92      !!      where the vertical advective flux, zwz, is given by :
93      !!            zwz = zcofk * zupst + (1-zcofk) * zcent
94      !!      with
95      !!        zupsv = upstream flux = wn * (ptb(k) or ptb(k-1) ) [wn>0 or <0]
96      !!        zcenu = centered flux = wn * mk(tn)
97      !!         The surface boundary condition is :
98      !!      variable volume (lk_vvl = T) : zero advective flux
99      !!      lin. free-surf  (lk_vvl = F) : wn(:,:,1) * ptn(:,:,1)
100      !!         Add this trend now to the general trend of tracer (ta,sa):
101      !!             pta = pta + ztra
102      !!         Trend diagnostic (l_trdtra=T or l_trctra=T): the trend is
103      !!      saved for diagnostics. The trends saved is expressed as :
104      !!             save trend =  w.gradz(T) = ztra - ptn divn.
105      !!
106      !! ** Action :  - update pta  with the now advective tracer trends
107      !!              - save trends if needed
108      !!----------------------------------------------------------------------
109      INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index
110      INTEGER                              , INTENT(in   ) ::   kit000          ! first time step index
111      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator)
112      INTEGER                              , INTENT(in   ) ::   kjpt            ! number of tracers
113      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pun, pvn, pwn   ! 3 ocean velocity components
114      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb, ptn        ! before and now tracer fields
115      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta             ! tracer trend
116      !
117      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
118      INTEGER  ::   ierr             ! local integer
119      REAL(wp) ::   zbtr, ztra                            ! local scalars
120      REAL(wp) ::   zfp_ui, zfp_vj, zfp_w, zcofi          !   -      -
121      REAL(wp) ::   zfm_ui, zfm_vj, zfm_w, zcofj, zcofk   !   -      -
122      REAL(wp) ::   zupsut, zcenut, zupst                 !   -      -
123      REAL(wp) ::   zupsvt, zcenvt, zcent, zice           !   -      -
124      REAL(wp), POINTER, DIMENSION(:,:)   :: zfzp         ! 2D workspace
125      REAL(wp), POINTER, DIMENSION(:,:,:) :: zwx, zwy     ! 3D     -
126      REAL(wp), POINTER, DIMENSION(:,:,:) :: zwz, zind    !  -     -
127      !!----------------------------------------------------------------------
128      !
129      IF( nn_timing == 1 )  CALL timing_start('tra_adv_cen2')
130      !
131      CALL wrk_alloc( jpi, jpj, zfzp )
132      CALL wrk_alloc( jpi, jpj, jpk, zwx, zwy, zwz, zind )
133      !
134
135      IF( kt == kit000 )  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         IF( .NOT. ALLOCATED( upsmsk ) )  THEN
142             ALLOCATE( upsmsk(jpi,jpj), STAT=ierr )
143             IF( ierr /= 0 )   CALL ctl_stop('STOP', 'tra_adv_cen2: unable to allocate array')
144         ENDIF
145
146         !
147         upsmsk(:,:) = 0._wp                             ! not upstream by default
148         !
149         IF( cp_cfg == "orca" )   CALL ups_orca_set      ! set mixed Upstream/centered scheme near some straits
150         !                                               ! and in closed seas (orca2 and orca4 only)
151         IF( jp_cfg == 2 .AND. .NOT. ln_rstart ) THEN    ! Increase the background in the surface layers
152            avmb(1) = 10.  * avmb(1)      ;      avtb(1) = 10.  * avtb(1)
153            avmb(2) = 10.  * avmb(2)      ;      avtb(2) = 10.  * avtb(2)
154            avmb(3) =  5.  * avmb(3)      ;      avtb(3) =  5.  * avtb(3)
155            avmb(4) =  2.5 * avmb(4)      ;      avtb(4) =  2.5 * avtb(4)
156         ENDIF
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      zfzp(:,:) = eos_fzp( 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) <= zfzp(ji,jj) + 0.1 ) THEN   ;   zice = 1._wp
169               ELSE                                                   ;   zice = 0._wp
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
255         IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR.    &
256            &( cdtype == 'TRC' .AND. l_trdtrc ) )   THEN
257            CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptn(:,:,:,jn) )
258            CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptn(:,:,:,jn) )
259            CALL trd_tra( kt, cdtype, jn, jptra_zad, zwz, pwn, ptn(:,:,:,jn) )
260         END IF
261         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
262         IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN 
263           IF( jn == jp_tem )   htr_adv(:) = ptr_vj( zwy(:,:,:) )
264           IF( jn == jp_sal )   str_adv(:) = ptr_vj( zwy(:,:,:) )
265         ENDIF
266         !
267      END DO
268
269      ! ---------------------------  required in restart file to ensure restartability)
270      ! avmb, avtb will be read in zdfini in restart case as they are used in zdftke, kpp etc...
271      IF( lrst_oce .AND. cdtype == 'TRA' ) THEN
272         CALL iom_rstput( kt, nitrst, numrow, 'avmb', avmb )
273         CALL iom_rstput( kt, nitrst, numrow, 'avtb', avtb )
274      ENDIF
275      !
276      CALL wrk_dealloc( jpi, jpj, zfzp )
277      CALL wrk_dealloc( jpi, jpj, jpk, zwx, zwy, zwz, zind )
278      !
279      IF( nn_timing == 1 )  CALL timing_stop('tra_adv_cen2')
280      !
281   END SUBROUTINE tra_adv_cen2
282   
283   
284   SUBROUTINE ups_orca_set
285      !!----------------------------------------------------------------------
286      !!                  ***  ROUTINE ups_orca_set  ***
287      !!       
288      !! ** Purpose :   add a portion of upstream scheme in area where the
289      !!                centered scheme generates too strong overshoot
290      !!
291      !! ** Method  :   orca (R4 and R2) confiiguration setting. Set upsmsk
292      !!                array to nozero value in some straith.
293      !!
294      !! ** Action : - upsmsk set to 1 at some strait, 0 elsewhere for orca
295      !!----------------------------------------------------------------------
296      INTEGER  ::   ii0, ii1, ij0, ij1      ! temporary integers
297      !!----------------------------------------------------------------------
298      !
299      IF( nn_timing == 1 )  CALL timing_start('ups_orca_set')
300      !
301      ! mixed upstream/centered scheme near river mouths
302      ! ------------------------------------------------
303      SELECT CASE ( jp_cfg )
304      !                                        ! =======================
305      CASE ( 4 )                               !  ORCA_R4 configuration
306         !                                     ! =======================
307         !                                          ! Gibraltar Strait
308         ii0 =  70   ;   ii1 =  71
309         ij0 =  52   ;   ij1 =  53   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50
310         !
311         !                                     ! =======================
312      CASE ( 2 )                               !  ORCA_R2 configuration
313         !                                     ! =======================
314         !                                          ! Gibraltar Strait
315         ij0 = 102   ;   ij1 = 102
316         ii0 = 138   ;   ii1 = 138   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.20
317         ii0 = 139   ;   ii1 = 139   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.40
318         ii0 = 140   ;   ii1 = 140   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50
319         ij0 = 101   ;   ij1 = 102
320         ii0 = 141   ;   ii1 = 141   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50
321         !                                          ! Bab el Mandeb Strait
322         ij0 =  87   ;   ij1 =  88
323         ii0 = 164   ;   ii1 = 164   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.10
324         ij0 =  88   ;   ij1 =  88
325         ii0 = 163   ;   ii1 = 163   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.25
326         ii0 = 162   ;   ii1 = 162   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.40
327         ii0 = 160   ;   ii1 = 161   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50
328         ij0 =  89   ;   ij1 =  89
329         ii0 = 158   ;   ii1 = 160   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.25
330         ij0 =  90   ;   ij1 =  90
331         ii0 = 160   ;   ii1 = 160   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.25
332         !                                          ! Sound Strait
333         ij0 = 116   ;   ij1 = 116
334         ii0 = 144   ;   ii1 = 144   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.25
335         ii0 = 145   ;   ii1 = 147   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50
336         ii0 = 148   ;   ii1 = 148   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.25
337         !
338      END SELECT 
339     
340      ! mixed upstream/centered scheme over closed seas
341      ! -----------------------------------------------
342      CALL clo_ups( upsmsk(:,:) )
343      !
344      IF( nn_timing == 1 )  CALL timing_stop('ups_orca_set')
345      !
346   END SUBROUTINE ups_orca_set
347
348   !!======================================================================
349END MODULE traadv_cen2
Note: See TracBrowser for help on using the repository browser.