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_muscl.F90 in branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_muscl.F90 @ 4401

Last change on this file since 4401 was 3211, checked in by spickles2, 12 years ago

Stephen Pickles, 11 Dec 2011

Commit to bring the rest of the DCSE NEMO development branch
in line with the latest development version. This includes
array index re-ordering of all OPA_SRC/.
  • Property svn:keywords set to Id
File size: 16.4 KB
Line 
1MODULE traadv_muscl
2   !!======================================================================
3   !!                       ***  MODULE  traadv_muscl  ***
4   !! Ocean  tracers:  horizontal & vertical advective trend
5   !!======================================================================
6   !! History :       !  2000-06  (A.Estublier)  for passive tracers
7   !!                 !  2001-08  (E.Durand, G.Madec)  adapted for T & S
8   !!   NEMO     1.0  !  2002-06  (G. Madec)  F90: Free form and module
9   !!            3.2  !  2010-05  (C. Ethe, G. Madec)  merge TRC-TRA + switch from velocity to transport
10   !!----------------------------------------------------------------------
11
12   !!----------------------------------------------------------------------
13   !!   tra_adv_muscl : update the tracer trend with the horizontal
14   !!                   and vertical advection trends using MUSCL scheme
15   !!----------------------------------------------------------------------
16   USE oce             ! ocean dynamics and active tracers
17   USE dom_oce         ! ocean space and time domain
18   USE trdmod_oce      ! tracers trends
19   USE trdtra      ! tracers trends
20   USE in_out_manager  ! I/O manager
21   USE dynspg_oce      ! choice/control of key cpp for surface pressure gradient
22   USE trabbl          ! tracers: bottom boundary layer
23   USE lib_mpp         ! distribued memory computing
24   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
25   USE diaptr          ! poleward transport diagnostics
26   USE trc_oce         ! share passive tracers/Ocean variables
27
28
29   IMPLICIT NONE
30   PRIVATE
31
32   PUBLIC   tra_adv_muscl  ! routine called by step.F90
33
34   LOGICAL  :: l_trd       ! flag to compute trends
35
36   !! * Control permutation of array indices
37#  include "oce_ftrans.h90"
38#  include "dom_oce_ftrans.h90"
39#  include "trc_oce_ftrans.h90"
40
41   !! * Substitutions
42#  include "domzgr_substitute.h90"
43#  include "vectopt_loop_substitute.h90"
44   !!----------------------------------------------------------------------
45   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
46   !! $Id$
47   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
48   !!----------------------------------------------------------------------
49CONTAINS
50
51   SUBROUTINE tra_adv_muscl( kt, cdtype, p2dt, pun, pvn, pwn, &
52      &                                        ptb, pta, kjpt )
53      !!----------------------------------------------------------------------
54      !!                    ***  ROUTINE tra_adv_muscl  ***
55      !!
56      !! ** Purpose :   Compute the now trend due to total advection of T and
57      !!      S using a MUSCL scheme (Monotone Upstream-centered Scheme for
58      !!      Conservation Laws) and add it to the general tracer trend.
59      !!
60      !! ** Method  : MUSCL scheme plus centered scheme at ocean boundaries
61      !!
62      !! ** Action  : - update (ta,sa) with the now advective tracer trends
63      !!              - save trends
64      !!
65      !! References : Estubier, A., and M. Levy, Notes Techn. Pole de Modelisation
66      !!              IPSL, Sept. 2000 (http://www.lodyc.jussieu.fr/opa)
67      !!----------------------------------------------------------------------
68      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
69      USE oce     , ONLY:   zwx   => ua       , zwy   => va          ! (ua,va) used as workspace
70      USE wrk_nemo, ONLY:   zslpx => wrk_3d_1 , zslpy => wrk_3d_2    ! 3D workspace
71
72      !! DCSE_NEMO: need additional directives for renamed module variables
73!FTRANS zwx zwy zslpx zslpy :I :I :z
74
75      !
76      INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index
77      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator)
78      INTEGER                              , INTENT(in   ) ::   kjpt            ! number of tracers
79      REAL(wp), DIMENSION(        jpk     ), INTENT(in   ) ::   p2dt            ! vertical profile of tracer time-step
80
81      !! DCSE_NEMO: This style defeats ftrans
82!     REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pun, pvn, pwn   ! 3 ocean velocity components
83!     REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb             ! before tracer field
84!     REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta             ! tracer trend
85
86!FTRANS pun pvn pwn :I :I :z
87!FTRANS ptb :I :I :z :
88!FTRANS pta :I :I :z :
89      REAL(wp), INTENT(in   ) ::   pun(jpi,jpj,jpk)        ! ocean velocity component (u)
90      REAL(wp), INTENT(in   ) ::   pvn(jpi,jpj,jpk)        ! ocean velocity component (v)
91      REAL(wp), INTENT(in   ) ::   pwn(jpi,jpj,jpk)        ! ocean velocity component (w)
92      REAL(wp), INTENT(in   ) ::   ptb(jpi,jpj,jpk,kjpt)   ! tracer fields (before)
93      REAL(wp), INTENT(inout) ::   pta(jpi,jpj,jpk,kjpt)   ! tracer trend
94
95      !
96      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
97      REAL(wp) ::   zu, z0u, zzwx, zw         ! local scalars
98      REAL(wp) ::   zv, z0v, zzwy, z0w        !   -      -
99      REAL(wp) ::   ztra, zbtr, zdt, zalpha   !   -      -
100      !!----------------------------------------------------------------------
101
102      IF( wrk_in_use(3, 1,2) ) THEN
103         CALL ctl_stop('tra_adv_muscl: requested workspace arrays unavailable')   ;   RETURN
104      ENDIF
105
106      IF( kt == nit000 )  THEN
107         IF(lwp) WRITE(numout,*)
108         IF(lwp) WRITE(numout,*) 'tra_adv : MUSCL advection scheme on ', cdtype
109         IF(lwp) WRITE(numout,*) '~~~~~~~'
110         !
111         l_trd = .FALSE.
112         IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.
113      ENDIF
114
115      !                                                     ! ===========
116      DO jn = 1, kjpt                                       ! tracer loop
117         !                                                  ! ===========
118         ! I. Horizontal advective fluxes
119         ! ------------------------------
120         ! first guess of the slopes
121         zwx(:,:,jpk) = 0.e0   ;   zwy(:,:,jpk) = 0.e0        ! bottom values
122         ! interior values
123#if defined key_z_first
124         DO jj = 1, jpjm1     
125            DO ji = 1, jpim1
126               DO jk = 1, jpkm1
127#else
128         DO jk = 1, jpkm1
129            DO jj = 1, jpjm1     
130               DO ji = 1, fs_jpim1   ! vector opt.
131#endif
132                  zwx(ji,jj,jk) = umask(ji,jj,jk) * ( ptb(ji+1,jj,jk,jn) - ptb(ji,jj,jk,jn) )
133                  zwy(ji,jj,jk) = vmask(ji,jj,jk) * ( ptb(ji,jj+1,jk,jn) - ptb(ji,jj,jk,jn) )
134               END DO
135           END DO
136         END DO
137         !
138         CALL lbc_lnk( zwx, 'U', -1. )                        ! lateral boundary conditions on zwx, zwy   (changed sign)
139         CALL lbc_lnk( zwy, 'V', -1. )
140         !                                             !-- Slopes of tracer
141         zslpx(:,:,jpk) = 0.e0   ;   zslpy(:,:,jpk) = 0.e0    ! bottom values
142#if defined key_z_first
143         DO jj = 2, jpj                                       ! interior values
144            DO ji = 2, jpi
145               DO jk = 1, jpkm1
146#else
147         DO jk = 1, jpkm1                                     ! interior values
148            DO jj = 2, jpj
149               DO ji = fs_2, jpi   ! vector opt.
150#endif
151                  zslpx(ji,jj,jk) =                    ( zwx(ji,jj,jk) + zwx(ji-1,jj  ,jk) )   &
152                     &            * ( 0.25 + SIGN( 0.25, zwx(ji,jj,jk) * zwx(ji-1,jj  ,jk) ) )
153                  zslpy(ji,jj,jk) =                    ( zwy(ji,jj,jk) + zwy(ji  ,jj-1,jk) )   &
154                     &            * ( 0.25 + SIGN( 0.25, zwy(ji,jj,jk) * zwy(ji  ,jj-1,jk) ) )
155               END DO
156            END DO
157         END DO
158         !
159#if defined key_z_first
160         DO jj = 2, jpj                                       ! Slopes limitation
161            DO ji = 2, jpi
162               DO jk = 1, jpkm1
163#else
164         DO jk = 1, jpkm1                                     ! Slopes limitation
165            DO jj = 2, jpj
166               DO ji = fs_2, jpi   ! vector opt.
167#endif
168                  zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) ) * MIN(    ABS( zslpx(ji  ,jj,jk) ),   &
169                     &                                                 2.*ABS( zwx  (ji-1,jj,jk) ),   &
170                     &                                                 2.*ABS( zwx  (ji  ,jj,jk) ) )
171                  zslpy(ji,jj,jk) = SIGN( 1., zslpy(ji,jj,jk) ) * MIN(    ABS( zslpy(ji,jj  ,jk) ),   &
172                     &                                                 2.*ABS( zwy  (ji,jj-1,jk) ),   &
173                     &                                                 2.*ABS( zwy  (ji,jj  ,jk) ) )
174               END DO
175           END DO
176         END DO             ! interior values
177
178         !                                             !-- MUSCL horizontal advective fluxes
179#if defined key_z_first
180         DO jj = 2, jpjm1                                     ! interior values
181            DO ji = 2, jpim1
182               DO jk = 1, jpkm1
183                  zdt  = p2dt(jk)
184#else
185         DO jk = 1, jpkm1                                     ! interior values
186            zdt  = p2dt(jk)
187            DO jj = 2, jpjm1
188               DO ji = fs_2, fs_jpim1   ! vector opt.
189#endif
190                  ! MUSCL fluxes
191                  z0u = SIGN( 0.5, pun(ji,jj,jk) )
192                  zalpha = 0.5 - z0u
193                  zu  = z0u - 0.5 * pun(ji,jj,jk) * zdt / ( e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) )
194                  zzwx = ptb(ji+1,jj,jk,jn) + zu * zslpx(ji+1,jj,jk)
195                  zzwy = ptb(ji  ,jj,jk,jn) + zu * zslpx(ji  ,jj,jk)
196                  zwx(ji,jj,jk) = pun(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy )
197                  !
198                  z0v = SIGN( 0.5, pvn(ji,jj,jk) )
199                  zalpha = 0.5 - z0v
200                  zv  = z0v - 0.5 * pvn(ji,jj,jk) * zdt / ( e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) )
201                  zzwx = ptb(ji,jj+1,jk,jn) + zv * zslpy(ji,jj+1,jk)
202                  zzwy = ptb(ji,jj  ,jk,jn) + zv * zslpy(ji,jj  ,jk) 
203                  zwy(ji,jj,jk) = pvn(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy )
204               END DO
205            END DO
206         END DO
207         !                                                    ! lateral boundary conditions on zwx, zwy   (changed sign)
208         CALL lbc_lnk( zwx, 'U', -1. )   ;   CALL lbc_lnk( zwy, 'V', -1. )
209         !
210         ! Tracer flux divergence at t-point added to the general trend
211#if defined key_z_first
212         DO jj = 2, jpjm1     
213            DO ji = 2, jpim1
214               DO jk = 1, jpkm1
215#else
216         DO jk = 1, jpkm1
217            DO jj = 2, jpjm1     
218               DO ji = fs_2, fs_jpim1   ! vector opt.
219#endif
220                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
221                  ! horizontal advective trends
222                  ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
223                  &               + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  ) )
224                  ! add it to the general tracer trends
225                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
226               END DO
227           END DO
228         END DO       
229         !                                 ! trend diagnostics (contribution of upstream fluxes)
230         IF( l_trd )  THEN
231            CALL trd_tra( kt, cdtype, jn, jptra_trd_xad, zwx, pun, ptb(:,:,:,jn) )
232            CALL trd_tra( kt, cdtype, jn, jptra_trd_yad, zwy, pvn, ptb(:,:,:,jn) )
233         END IF
234         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
235         IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN 
236            IF( jn == jp_tem )  htr_adv(:) = ptr_vj( zwy(:,:,:) )
237            IF( jn == jp_sal )  str_adv(:) = ptr_vj( zwy(:,:,:) )
238         ENDIF
239
240         ! II. Vertical advective fluxes
241         ! -----------------------------
242         !                                             !-- first guess of the slopes
243#if defined key_z_first
244         DO jj = 1, jpj
245            DO ji = 1, jpi
246               zwx(ji,jj,1) = 0.e0                             ! surface boundary conditions
247               DO jk = 2, jpkm1                                ! interior values
248                  zwx(ji,jj,jk) = tmask(ji,jj,jk) * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) )
249               END DO
250               zwx(ji,jj,jpk) = 0.e0                           ! bottom boundary conditions
251            END DO
252         END DO
253#else
254         zwx (:,:, 1 ) = 0.e0    ;    zwx (:,:,jpk) = 0.e0    ! surface & bottom boundary conditions
255         DO jk = 2, jpkm1                                     ! interior values
256            zwx(:,:,jk) = tmask(:,:,jk) * ( ptb(:,:,jk-1,jn) - ptb(:,:,jk,jn) )
257         END DO
258#endif
259
260         !                                             !-- Slopes of tracer
261#if defined key_z_first
262         DO jj = 1, jpj
263            DO ji = 1, jpi
264               zslpx(ji,jj,1) = 0.e0                          ! surface values
265               DO jk = 2, jpkm1                               ! interior value
266#else
267         zslpx(:,:,1) = 0.e0                                  ! surface values
268         DO jk = 2, jpkm1                                     ! interior value
269            DO jj = 1, jpj
270               DO ji = 1, jpi
271#endif
272                  zslpx(ji,jj,jk) =                    ( zwx(ji,jj,jk) + zwx(ji,jj,jk+1) )   &
273                     &            * ( 0.25 + SIGN( 0.25, zwx(ji,jj,jk) * zwx(ji,jj,jk+1) ) )
274               END DO
275            END DO
276         END DO
277         !                                             !-- Slopes limitation
278#if defined key_z_first
279         DO jj = 1, jpj   
280            DO ji = 1, jpi
281               DO jk = 2, jpkm1                               ! interior values
282#else
283         DO jk = 2, jpkm1                                     ! interior values
284            DO jj = 1, jpj
285               DO ji = 1, jpi
286#endif
287                  zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) ) * MIN(    ABS( zslpx(ji,jj,jk  ) ),   &
288                     &                                                 2.*ABS( zwx  (ji,jj,jk+1) ),   &
289                     &                                                 2.*ABS( zwx  (ji,jj,jk  ) )  )
290               END DO
291            END DO
292         END DO
293         !                                             !-- vertical advective flux
294         !                                                    ! surface values  (bottom already set to zero)
295         IF( lk_vvl ) THEN    ;   zwx(:,:, 1 ) = 0.e0                      !  variable volume
296         ELSE                 ;   zwx(:,:, 1 ) = pwn(:,:,1) * ptb(:,:,1,jn)   ! linear free surface
297         ENDIF 
298         !
299#if defined key_z_first
300         DO jj = 2, jpjm1                                     ! interior values
301            DO ji = 2, jpim1
302               DO jk = 1, jpkm1
303                  zdt  = p2dt(jk)
304#else
305         DO jk = 1, jpkm1                                     ! interior values
306            zdt  = p2dt(jk)
307            DO jj = 2, jpjm1     
308               DO ji = fs_2, fs_jpim1   ! vector opt.
309#endif
310                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3w(ji,jj,jk+1) )
311                  z0w = SIGN( 0.5, pwn(ji,jj,jk+1) )
312                  zalpha = 0.5 + z0w
313                  zw  = z0w - 0.5 * pwn(ji,jj,jk+1) * zdt * zbtr 
314                  zzwx = ptb(ji,jj,jk+1,jn) + zw * zslpx(ji,jj,jk+1)
315                  zzwy = ptb(ji,jj,jk  ,jn) + zw * zslpx(ji,jj,jk  )
316                  zwx(ji,jj,jk+1) = pwn(ji,jj,jk+1) * ( zalpha * zzwx + (1.-zalpha) * zzwy )
317               END DO
318            END DO
319         END DO
320
321         ! Compute & add the vertical advective trend
322#if defined key_z_first
323         DO jj = 2, jpjm1     
324            DO ji = 2, jpim1
325               DO jk = 1, jpkm1
326#else
327         DO jk = 1, jpkm1
328            DO jj = 2, jpjm1     
329               DO ji = fs_2, fs_jpim1   ! vector opt.
330#endif
331                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
332                  ! vertical advective trends
333                  ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) )
334                  ! add it to the general tracer trends
335                  pta(ji,jj,jk,jn) =  pta(ji,jj,jk,jn) + ztra
336               END DO
337            END DO
338         END DO
339         !                                 ! Save the vertical advective trends for diagnostic
340         IF( l_trd ) CALL trd_tra( kt, cdtype, jn, jptra_trd_zad, zwx, pwn, ptb(:,:,:,jn) )
341         !
342      ENDDO
343      !
344      IF( wrk_not_released(3, 1,2) )   CALL ctl_stop('tra_adv_muscl: requested workspace arrays unavailable')
345      !
346   END SUBROUTINE tra_adv_muscl
347
348   !!======================================================================
349END MODULE traadv_muscl
Note: See TracBrowser for help on using the repository browser.