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

source: branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_cen2.F90 @ 2830

Last change on this file since 2830 was 2830, checked in by kpedwards, 13 years ago

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