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

source: trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_cen2.F90 @ 3294

Last change on this file since 3294 was 3294, checked in by rblod, 12 years ago

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