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