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

source: branches/2012/dev_LOCEAN_UKMO_2012/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_cen2.F90 @ 3653

Last change on this file since 3653 was 3653, checked in by cetlod, 11 years ago

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