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

source: trunk/NEMO/OPA_SRC/TRA/traadv_muscl2.F90 @ 1110

Last change on this file since 1110 was 719, checked in by ctlod, 17 years ago

get back to the nemo_v2_3 version for trunk
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 22.5 KB
Line 
1MODULE traadv_muscl2
2   !!==============================================================================
3   !!                       ***  MODULE  traadv_muscl2  ***
4   !! Ocean active tracers:  horizontal & vertical advective trend
5   !!==============================================================================
6   !! History :  9.0  !  02-06  (G. Madec) from traadv_muscl
7   !!----------------------------------------------------------------------
8
9   !!----------------------------------------------------------------------
10   !!   tra_adv_muscl2 : update the tracer trend with the horizontal
11   !!                    and vertical advection trends using MUSCL2 scheme
12   !!----------------------------------------------------------------------
13   USE oce             ! ocean dynamics and active tracers
14   USE dom_oce         ! ocean space and time domain
15   USE trdmod          ! ocean active tracers trends
16   USE trdmod_oce      ! ocean variables trends
17   USE in_out_manager  ! I/O manager
18   USE dynspg_oce      ! choice/control of key cpp for surface pressure gradient
19   USE trabbl          ! tracers: bottom boundary layer
20   USE lib_mpp
21   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
22   USE diaptr          ! poleward transport diagnostics
23   USE prtctl          ! Print control
24
25   IMPLICIT NONE
26   PRIVATE
27
28   !! * Accessibility
29   PUBLIC tra_adv_muscl2        ! routine called by step.F90
30
31   !! * Substitutions
32#  include "domzgr_substitute.h90"
33#  include "vectopt_loop_substitute.h90"
34   !!----------------------------------------------------------------------
35   !!   OPA 9.0 , LOCEAN-IPSL (2006)
36   !! $Header$
37   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
38   !!----------------------------------------------------------------------
39
40CONTAINS
41
42   SUBROUTINE tra_adv_muscl2( kt, pun, pvn, pwn )
43      !!----------------------------------------------------------------------
44      !!                   ***  ROUTINE tra_adv_muscl2  ***
45      !!
46      !! ** Purpose :   Compute the now trend due to total advection of T and
47      !!      S using a MUSCL scheme (Monotone Upstream-centered Scheme for
48      !!      Conservation Laws) and add it to the general tracer trend.
49      !!
50      !! ** Method  : MUSCL scheme plus centered scheme at ocean boundaries
51      !!
52      !! ** Action  : - update (ta,sa) with the now advective tracer trends
53      !!              - save trends in (ztrdt,ztrds) ('key_trdtra')
54      !!
55      !! References : Estubier, A., and M. Levy, Notes Techn. Pole de Modelisation
56      !!              IPSL, Sept. 2000 (http://www.lodyc.jussieu.fr/opa)
57      !!----------------------------------------------------------------------
58      USE oce              , ztrdt => ua   ! use ua as workspace
59      USE oce              , ztrds => va   ! use va as workspace
60      !!
61      INTEGER , INTENT(in)                         ::   kt    ! ocean time-step index
62      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::   pun   ! ocean velocity u-component
63      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::   pvn   ! ocean velocity v-component
64      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::   pwn   ! ocean velocity w-component
65      !!
66      INTEGER ::   ji, jj, jk   ! dummy loop indices
67      REAL(wp) ::   &
68         zu, zv, zw, zeu, zev,           & 
69         zew, zbtr, zstep,               &
70         z0u, z0v, z0w,                  &
71         zzt1, zzt2, zalpha,             &
72         zzs1, zzs2, z2,                 &
73         zta, zsa,                       &
74         z_hdivn_x, z_hdivn_y, z_hdivn
75      REAL(wp), DIMENSION (jpi,jpj,jpk) ::   zt1, zt2, ztp1, ztp2   ! 3D workspace
76      REAL(wp), DIMENSION (jpi,jpj,jpk) ::   zs1, zs2, zsp1, zsp2   !  "      "
77      !!----------------------------------------------------------------------
78
79      IF( kt == nit000 .AND. lwp ) THEN
80         WRITE(numout,*)
81         WRITE(numout,*) 'tra_adv_muscl2 : MUSCL2 advection scheme'
82         WRITE(numout,*) '~~~~~~~~~~~~~~~'
83      ENDIF
84
85      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   z2 = 1.
86      ELSE                                        ;   z2 = 2.
87      ENDIF
88
89      ! I. Horizontal advective fluxes
90      ! ------------------------------
91
92      ! first guess of the slopes
93      ! interior values
94      DO jk = 1, jpkm1
95         DO jj = 1, jpjm1     
96            DO ji = 1, fs_jpim1   ! vector opt.
97               zt1(ji,jj,jk) = umask(ji,jj,jk) * ( tb(ji+1,jj,jk) - tb(ji,jj,jk) )
98               zs1(ji,jj,jk) = umask(ji,jj,jk) * ( sb(ji+1,jj,jk) - sb(ji,jj,jk) )
99               zt2(ji,jj,jk) = vmask(ji,jj,jk) * ( tb(ji,jj+1,jk) - tb(ji,jj,jk) )
100               zs2(ji,jj,jk) = vmask(ji,jj,jk) * ( sb(ji,jj+1,jk) - sb(ji,jj,jk) )
101            END DO
102         END DO
103      END DO
104      ! bottom values
105      zt1(:,:,jpk) = 0.e0    ;    zt2(:,:,jpk) = 0.e0
106      zs1(:,:,jpk) = 0.e0    ;    zs2(:,:,jpk) = 0.e0
107
108      ! lateral boundary conditions on zt1, zt2 ; zs1, zs2   (changed sign)
109      CALL lbc_lnk( zt1, 'U', -1. )   ;   CALL lbc_lnk( zs1, 'U', -1. )
110      CALL lbc_lnk( zt2, 'V', -1. )   ;   CALL lbc_lnk( zs2, 'V', -1. )
111
112      ! Slopes
113      ! interior values
114      DO jk = 1, jpkm1
115         DO jj = 2, jpj
116            DO ji = fs_2, jpi   ! vector opt.
117               ztp1(ji,jj,jk) =                    ( zt1(ji,jj,jk) + zt1(ji-1,jj  ,jk) )   &
118                  &           * ( 0.25 + SIGN( 0.25, zt1(ji,jj,jk) * zt1(ji-1,jj  ,jk) ) )
119               zsp1(ji,jj,jk) =                    ( zs1(ji,jj,jk) + zs1(ji-1,jj  ,jk) )   &
120                  &           * ( 0.25 + SIGN( 0.25, zs1(ji,jj,jk) * zs1(ji-1,jj  ,jk) ) )
121               ztp2(ji,jj,jk) =                    ( zt2(ji,jj,jk) + zt2(ji  ,jj-1,jk) )   &
122                  &           * ( 0.25 + SIGN( 0.25, zt2(ji,jj,jk) * zt2(ji  ,jj-1,jk) ) )
123               zsp2(ji,jj,jk) =                    ( zs2(ji,jj,jk) + zs2(ji  ,jj-1,jk) )   &
124                  &           * ( 0.25 + SIGN( 0.25, zs2(ji,jj,jk) * zs2(ji  ,jj-1,jk) ) )
125            END DO
126         END DO
127      END DO
128      ! bottom values
129      ztp1(:,:,jpk) = 0.e0    ;    ztp2(:,:,jpk) = 0.e0
130      zsp1(:,:,jpk) = 0.e0    ;    zsp2(:,:,jpk) = 0.e0
131
132      ! Slopes limitation
133      DO jk = 1, jpkm1
134         DO jj = 2, jpj
135            DO ji = fs_2, jpi   ! vector opt.
136               ztp1(ji,jj,jk) = SIGN( 1., ztp1(ji,jj,jk) )   &
137                  &           * MIN(    ABS( ztp1(ji  ,jj,jk) ),   &
138                  &                  2.*ABS( zt1 (ji-1,jj,jk) ),   &
139                  &                  2.*ABS( zt1 (ji  ,jj,jk) ) )
140               zsp1(ji,jj,jk) = SIGN( 1., zsp1(ji,jj,jk) )   &
141                  &           * MIN(    ABS( zsp1(ji  ,jj,jk) ),   &
142                  &                  2.*ABS( zs1 (ji-1,jj,jk) ),   &
143                  &                  2.*ABS( zs1 (ji  ,jj,jk) ) )
144               ztp2(ji,jj,jk) = SIGN( 1., ztp2(ji,jj,jk) )   &
145                  &           * MIN(    ABS( ztp2(ji,jj  ,jk) ),   &
146                  &                  2.*ABS( zt2 (ji,jj-1,jk) ),   &
147                  &                  2.*ABS( zt2 (ji,jj  ,jk) ) )
148               zsp2(ji,jj,jk) = SIGN( 1., zsp2(ji,jj,jk) )   &
149                  &           * MIN(    ABS( zsp2(ji,jj  ,jk) ),   &
150                  &                  2.*ABS( zs2 (ji,jj-1,jk) ),   &
151                  &                  2.*ABS( zs2 (ji,jj  ,jk) ) )
152            END DO
153         END DO
154      END DO       
155
156      ! Advection terms
157      ! interior values
158      DO jk = 1, jpkm1
159         zstep  = z2 * rdttra(jk)
160         DO jj = 2, jpjm1     
161            DO ji = fs_2, fs_jpim1   ! vector opt.
162               ! volume fluxes
163#if defined key_zco
164               zeu = e2u(ji,jj)                   * pun(ji,jj,jk)
165               zev = e1v(ji,jj)                   * pvn(ji,jj,jk)
166#else
167               zeu = e2u(ji,jj) * fse3u(ji,jj,jk) * pun(ji,jj,jk)
168               zev = e1v(ji,jj) * fse3v(ji,jj,jk) * pvn(ji,jj,jk)
169#endif
170               ! MUSCL fluxes
171               z0u = SIGN( 0.5, pun(ji,jj,jk) )           
172               zalpha = 0.5 - z0u
173               zu  = z0u - 0.5 * pun(ji,jj,jk) * zstep / e1u(ji,jj)
174               zzt1 = tb(ji+1,jj,jk) + zu*ztp1(ji+1,jj,jk)
175               zzt2 = tb(ji  ,jj,jk) + zu*ztp1(ji  ,jj,jk)
176               zzs1 = sb(ji+1,jj,jk) + zu*zsp1(ji+1,jj,jk)
177               zzs2 = sb(ji  ,jj,jk) + zu*zsp1(ji  ,jj,jk)
178               zt1(ji,jj,jk) = zeu * ( zalpha * zzt1 + (1.-zalpha) * zzt2 )
179               zs1(ji,jj,jk) = zeu * ( zalpha * zzs1 + (1.-zalpha) * zzs2 )
180               !
181               z0v = SIGN( 0.5, pvn(ji,jj,jk) )           
182               zalpha = 0.5 - z0v
183               zv  = z0v - 0.5 * pvn(ji,jj,jk) * zstep / e2v(ji,jj)
184               zzt1 = tb(ji,jj+1,jk) + zv*ztp2(ji,jj+1,jk)
185               zzt2 = tb(ji,jj  ,jk) + zv*ztp2(ji,jj  ,jk)
186               zzs1 = sb(ji,jj+1,jk) + zv*zsp2(ji,jj+1,jk)
187               zzs2 = sb(ji,jj  ,jk) + zv*zsp2(ji,jj  ,jk)
188               zt2(ji,jj,jk) = zev * ( zalpha * zzt1 + (1.-zalpha) * zzt2 )
189               zs2(ji,jj,jk) = zev * ( zalpha * zzs1 + (1.-zalpha) * zzs2 )
190            END DO
191         END DO
192      END DO
193
194      !!!!  centered scheme at lateral b.C. if off-shore velocity
195      DO jk = 1, jpkm1
196        DO jj = 2, jpjm1
197            DO ji = fs_2, fs_jpim1   ! vector opt.
198#if defined key_zco
199               IF( umask(ji,jj,jk) == 0. ) THEN
200                  IF( pun(ji+1,jj,jk) > 0. .AND. ji /= jpi ) THEN
201                     zt1(ji+1,jj,jk) = e2u(ji+1,jj) * pun(ji+1,jj,jk) * ( tb(ji+1,jj,jk) + tb(ji+2,jj,jk) ) * 0.5
202                     zs1(ji+1,jj,jk) = e2u(ji+1,jj) * pun(ji+1,jj,jk) * ( sb(ji+1,jj,jk) + sb(ji+2,jj,jk) ) * 0.5
203                  ENDIF
204                  IF( pun(ji-1,jj,jk) < 0. ) THEN
205                     zt1(ji-1,jj,jk) = e2u(ji-1,jj) * pun(ji-1,jj,jk) * ( tb(ji-1,jj,jk) + tb(ji  ,jj,jk) ) * 0.5
206                     zs1(ji-1,jj,jk) = e2u(ji-1,jj) * pun(ji-1,jj,jk) * ( sb(ji-1,jj,jk) + sb(ji  ,jj,jk) ) * 0.5
207                  ENDIF
208               ENDIF
209               IF( vmask(ji,jj,jk) == 0. ) THEN
210                  IF( pvn(ji,jj+1,jk) > 0. .AND. jj /= jpj ) THEN
211                     zt2(ji,jj+1,jk) = e1v(ji,jj+1) * pvn(ji,jj+1,jk) * ( tb(ji,jj+1,jk) + tb(ji,jj+2,jk) ) * 0.5
212                     zs2(ji,jj+1,jk) = e1v(ji,jj+1) * pvn(ji,jj+1,jk) * ( sb(ji,jj+1,jk) + sb(ji,jj+2,jk) ) * 0.5
213                  ENDIF
214                  IF( pvn(ji,jj-1,jk) < 0. ) THEN
215                     zt2(ji,jj-1,jk) = e1v(ji,jj-1) * pvn(ji,jj-1,jk) * ( tb(ji,jj-1,jk) + tb(ji  ,jj,jk) ) * 0.5
216                     zs2(ji,jj-1,jk) = e1v(ji,jj-1) * pvn(ji,jj-1,jk) * ( sb(ji,jj-1,jk) + sb(ji  ,jj,jk) ) * 0.5
217                  ENDIF
218               ENDIF
219#else
220               IF( umask(ji,jj,jk) == 0. ) THEN
221                  IF( pun(ji+1,jj,jk) > 0. .AND. ji /= jpi ) THEN
222                     zt1(ji+1,jj,jk) = e2u(ji+1,jj)* fse3u(ji+1,jj,jk)   &
223                        &            * pun(ji+1,jj,jk) * ( tb(ji+1,jj,jk) + tb(ji+2,jj,jk) ) * 0.5
224                     zs1(ji+1,jj,jk) = e2u(ji+1,jj)* fse3u(ji+1,jj,jk)   &
225                        &            * pun(ji+1,jj,jk) * ( sb(ji+1,jj,jk) + sb(ji+2,jj,jk) ) * 0.5
226                  ENDIF
227                  IF( pun(ji-1,jj,jk) < 0. ) THEN
228                     zt1(ji-1,jj,jk) = e2u(ji-1,jj)* fse3u(ji-1,jj,jk)   &
229                        &            * pun(ji-1,jj,jk) * ( tb(ji-1,jj,jk) + tb(ji  ,jj,jk) ) * 0.5
230                     zs1(ji-1,jj,jk) = e2u(ji-1,jj)* fse3u(ji-1,jj,jk)   &
231                        &            * pun(ji-1,jj,jk) * ( sb(ji-1,jj,jk) + sb(ji  ,jj,jk) ) * 0.5
232                  ENDIF
233               ENDIF
234               IF( vmask(ji,jj,jk) == 0. ) THEN
235                  IF( pvn(ji,jj+1,jk) > 0. .AND. jj /= jpj ) THEN
236                     zt2(ji,jj+1,jk) = e1v(ji,jj+1) * fse3v(ji,jj+1,jk)   &
237                        &            * pvn(ji,jj+1,jk) * ( tb(ji,jj+1,jk) + tb(ji,jj+2,jk) ) * 0.5
238                     zs2(ji,jj+1,jk) = e1v(ji,jj+1) * fse3v(ji,jj+1,jk)   &
239                        &            * pvn(ji,jj+1,jk) * ( sb(ji,jj+1,jk) + sb(ji,jj+2,jk) ) * 0.5
240                  ENDIF
241                  IF( pvn(ji,jj-1,jk) < 0. ) THEN
242                     zt2(ji,jj-1,jk) = e1v(ji,jj-1)* fse3v(ji,jj-1,jk)   &
243                        &            * pvn(ji,jj-1,jk) * ( tb(ji,jj-1,jk) + tb(ji  ,jj,jk) ) * 0.5
244                     zs2(ji,jj-1,jk) = e1v(ji,jj-1)* fse3v(ji,jj-1,jk)   &
245                        &            * pvn(ji,jj-1,jk) * ( sb(ji,jj-1,jk) + sb(ji  ,jj,jk) ) * 0.5
246                  ENDIF
247               ENDIF
248#endif
249            END DO
250         END DO
251      END DO
252
253      ! lateral boundary conditions on zt1, zt2 ; zs1, zs2   (changed sign)
254      CALL lbc_lnk( zt1, 'U', -1. )   ;   CALL lbc_lnk( zs1, 'U', -1. ) 
255      CALL lbc_lnk( zt2, 'V', -1. )   ;   CALL lbc_lnk( zs2, 'V', -1. )
256
257      ! Compute & add the horizontal advective trend
258
259      DO jk = 1, jpkm1
260         DO jj = 2, jpjm1     
261            DO ji = fs_2, fs_jpim1   ! vector opt.
262#if defined key_zco
263               zbtr = 1. / ( e1t(ji,jj)*e2t(ji,jj) )
264#else
265               zbtr = 1. / ( e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk) )
266#endif
267               ! horizontal advective trends
268               zta = - zbtr * ( zt1(ji,jj,jk) - zt1(ji-1,jj  ,jk  )   &
269                  &           + zt2(ji,jj,jk) - zt2(ji  ,jj-1,jk  ) )
270               zsa = - zbtr * ( zs1(ji,jj,jk) - zs1(ji-1,jj  ,jk  )   &
271                  &           + zs2(ji,jj,jk) - zs2(ji  ,jj-1,jk  ) ) 
272               ! add it to the general tracer trends
273               ta(ji,jj,jk) = ta(ji,jj,jk) + zta
274               sa(ji,jj,jk) = sa(ji,jj,jk) + zsa
275            END DO
276        END DO
277      END DO       
278
279      ! Save the horizontal advective trends for diagnostic
280      IF( l_trdtra ) THEN
281         ztrdt(:,:,:) = 0.e0   ;   ztrds(:,:,:) = 0.e0
282         !
283         ! T/S ZONAL advection trends
284         DO jk = 1, jpkm1
285            DO jj = 2, jpjm1
286               DO ji = fs_2, fs_jpim1   ! vector opt.
287                  !-- Compute zonal divergence by splitting hdivn (see divcur.F90)
288                  !   N.B. This computation is not valid along OBCs (if any)
289#if defined key_zco
290                  zbtr      = 1. / ( e1t(ji,jj) * e2t(ji,jj) )
291                  z_hdivn_x = (  e2u(ji  ,jj) * pun(ji  ,jj,jk)                              &
292                     &         - e2u(ji-1,jj) * pun(ji-1,jj,jk) ) * zbtr
293#else
294                  zbtr      = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
295                  z_hdivn_x = (  e2u(ji  ,jj) * fse3u(ji  ,jj,jk) * pun(ji  ,jj,jk)          &
296                     &         - e2u(ji-1,jj) * fse3u(ji-1,jj,jk) * pun(ji-1,jj,jk) ) * zbtr
297#endif
298                  ztrdt(ji,jj,jk) = - zbtr * ( zt1(ji,jj,jk) - zt1(ji-1,jj,jk) ) + tn(ji,jj,jk) * z_hdivn_x
299                  ztrds(ji,jj,jk) = - zbtr * ( zs1(ji,jj,jk) - zs1(ji-1,jj,jk) ) + sn(ji,jj,jk) * z_hdivn_x
300               END DO
301            END DO
302         END DO
303         CALL trd_mod(ztrdt, ztrds, jptra_trd_xad, 'TRA', kt)
304
305         ! T/S MERIDIONAL advection trends
306         DO jk = 1, jpkm1
307            DO jj = 2, jpjm1
308               DO ji = fs_2, fs_jpim1   ! vector opt.
309                  !-- Compute merid. divergence by splitting hdivn (see divcur.F90)
310                  !   N.B. This computation is not valid along OBCs (if any)
311#if defined key_zco
312                  zbtr      = 1. / ( e1t(ji,jj) * e2t(ji,jj) )
313                  z_hdivn_y = (  e1v(ji,jj  ) * pvn(ji,jj  ,jk)                              &
314                     &         - e1v(ji,jj-1) * pvn(ji,jj-1,jk) ) * zbtr
315#else
316                  zbtr      = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
317                  z_hdivn_y = (  e1v(ji,  jj) * fse3v(ji,jj  ,jk) * pvn(ji,jj  ,jk)          &
318                     &         - e1v(ji,jj-1) * fse3v(ji,jj-1,jk) * pvn(ji,jj-1,jk) ) * zbtr
319#endif
320                  ztrdt(ji,jj,jk) = - zbtr * ( zt2(ji,jj,jk) - zt2(ji,jj-1,jk) ) + tn(ji,jj,jk) * z_hdivn_y         
321                  ztrds(ji,jj,jk) = - zbtr * ( zs2(ji,jj,jk) - zs2(ji,jj-1,jk) ) + sn(ji,jj,jk) * z_hdivn_y
322               END DO
323            END DO
324         END DO
325         CALL trd_mod(ztrdt, ztrds, jptra_trd_yad, 'TRA', kt)
326
327         ! Save the up-to-date ta and sa trends
328         ztrdt(:,:,:) = ta(:,:,:) 
329         ztrds(:,:,:) = sa(:,:,:) 
330         !
331      ENDIF
332
333      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ta, clinfo1=' muscl2 had  - Ta: ', mask1=tmask,   &
334         &                       tab3d_2=sa, clinfo2=              ' Sa: ', mask2=tmask, clinfo3='tra')
335
336      ! "zonal" mean advective heat and salt transport
337      IF( ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN
338         IF( lk_zco ) THEN
339            DO jk = 1, jpkm1
340               DO jj = 2, jpjm1
341                  DO ji = fs_2, fs_jpim1   ! vector opt.
342                    zt2(ji,jj,jk) = zt2(ji,jj,jk) * fse3v(ji,jj,jk)
343                    zs2(ji,jj,jk) = zs2(ji,jj,jk) * fse3v(ji,jj,jk)
344                  END DO
345               END DO
346            END DO
347         ENDIF
348         pht_adv(:) = ptr_vj( zt2(:,:,:) )
349         pst_adv(:) = ptr_vj( zs2(:,:,:) )
350      ENDIF
351
352      ! II. Vertical advective fluxes
353      ! -----------------------------
354     
355      ! First guess of the slope
356      ! interior values
357      DO jk = 2, jpkm1
358         zt1(:,:,jk) = tmask(:,:,jk) * ( tb(:,:,jk-1) - tb(:,:,jk) )
359         zs1(:,:,jk) = tmask(:,:,jk) * ( sb(:,:,jk-1) - sb(:,:,jk) )
360      END DO
361      ! surface & bottom boundary conditions
362      zt1 (:,:, 1 ) = 0.e0    ;    zt1 (:,:,jpk) = 0.e0
363      zs1 (:,:, 1 ) = 0.e0    ;    zs1 (:,:,jpk) = 0.e0
364
365      ! Slopes
366      DO jk = 2, jpkm1
367         DO jj = 1, jpj
368            DO ji = 1, jpi
369               ztp1(ji,jj,jk) =                    ( zt1(ji,jj,jk) + zt1(ji,jj,jk+1) )   &
370                  &           * ( 0.25 + SIGN( 0.25, zt1(ji,jj,jk) * zt1(ji,jj,jk+1) ) )
371               zsp1(ji,jj,jk) =                    ( zs1(ji,jj,jk) + zs1(ji,jj,jk+1) )   &
372                  &           * ( 0.25 + SIGN( 0.25, zs1(ji,jj,jk) * zs1(ji,jj,jk+1) ) )
373            END DO
374         END DO
375      END DO
376
377      ! Slopes limitation
378      ! interior values
379      DO jk = 2, jpkm1
380         DO jj = 1, jpj
381            DO ji = 1, jpi
382               ztp1(ji,jj,jk) = SIGN( 1., ztp1(ji,jj,jk) )   &
383                  &           * MIN(    ABS( ztp1(ji,jj,jk  ) ),   &
384                  &                  2.*ABS( zt1 (ji,jj,jk+1) ),   &
385                  &                  2.*ABS( zt1 (ji,jj,jk  ) ) )
386               zsp1(ji,jj,jk) = SIGN( 1., zsp1(ji,jj,jk) )   &
387                  &           * MIN(    ABS( zsp1(ji,jj,jk  ) ),   &
388                  &                  2.*ABS( zs1 (ji,jj,jk+1) ),   &
389                  &                  2.*ABS( zs1 (ji,jj,jk  ) ) )
390            END DO
391         END DO
392      END DO
393      ! surface values
394      ztp1(:,:,1) = 0.e0
395      zsp1(:,:,1) = 0.e0
396
397      ! vertical advective flux
398      ! interior values
399      DO jk = 1, jpkm1
400         zstep  = z2 * rdttra(jk)
401         DO jj = 2, jpjm1     
402            DO ji = fs_2, fs_jpim1   ! vector opt.
403               zew = pwn(ji,jj,jk+1)
404               z0w = SIGN( 0.5, pwn(ji,jj,jk+1) )
405               zalpha = 0.5 + z0w
406               zw  = z0w - 0.5 * pwn(ji,jj,jk+1)*zstep / fse3w(ji,jj,jk+1)
407               zzt1 = tb(ji,jj,jk+1) + zw*ztp1(ji,jj,jk+1)
408               zzt2 = tb(ji,jj,jk  ) + zw*ztp1(ji,jj,jk  )
409               zzs1 = sb(ji,jj,jk+1) + zw*zsp1(ji,jj,jk+1)
410               zzs2 = sb(ji,jj,jk  ) + zw*zsp1(ji,jj,jk  )
411               zt1(ji,jj,jk+1) = zew * ( zalpha * zzt1 + (1.-zalpha)*zzt2 )
412               zs1(ji,jj,jk+1) = zew * ( zalpha * zzs1 + (1.-zalpha)*zzs2 )
413            END DO
414         END DO
415      END DO
416      DO jk = 2, jpkm1
417        DO jj = 2, jpjm1
418            DO ji = fs_2, fs_jpim1   ! vector opt.
419               IF( tmask(ji,jj,jk+1) == 0. ) THEN
420                  IF( pwn(ji,jj,jk) > 0. ) THEN
421                     zt1(ji,jj,jk) = pwn(ji,jj,jk) * ( tb(ji,jj,jk-1) + tb(ji,jj,jk) ) * 0.5
422                     zs1(ji,jj,jk) = pwn(ji,jj,jk) * ( sb(ji,jj,jk-1) + sb(ji,jj,jk) ) * 0.5
423                  ENDIF
424               ENDIF
425            END DO
426         END DO
427      END DO
428
429      ! surface values
430      IF( lk_dynspg_rl .OR. lk_vvl ) THEN               ! rigid lid or variable volume: flux set to zero
431         zt1(:,:, 1 ) = 0.e0
432         zs1(:,:, 1 ) = 0.e0
433      ELSE                                              ! free surface
434         zt1(:,:, 1 ) = pwn(:,:,1) * tb(:,:,1)
435         zs1(:,:, 1 ) = pwn(:,:,1) * sb(:,:,1)
436      ENDIF
437
438      ! bottom values
439      zt1(:,:,jpk) = 0.e0
440      zs1(:,:,jpk) = 0.e0
441
442
443      ! Compute & add the vertical advective trend
444
445      DO jk = 1, jpkm1
446         DO jj = 2, jpjm1     
447            DO ji = fs_2, fs_jpim1   ! vector opt.
448               zbtr = 1. / fse3t(ji,jj,jk)
449               ! horizontal advective trends
450               zta = - zbtr * ( zt1(ji,jj,jk) - zt1(ji,jj,jk+1) )
451               zsa = - zbtr * ( zs1(ji,jj,jk) - zs1(ji,jj,jk+1) )
452               ! add it to the general tracer trends
453               ta(ji,jj,jk) =  ta(ji,jj,jk) + zta
454               sa(ji,jj,jk) =  sa(ji,jj,jk) + zsa
455            END DO
456         END DO
457      END DO
458
459      ! Save the vertical advective trends for diagnostic
460      IF( l_trdtra )   THEN
461         ! Recompute the vertical advection zta & zsa trends computed
462         ! at the step 2. above in making the difference between the new
463         ! trends and the previous one: ta()/sa - ztrdt()/ztrds() and substract
464         ! the term tn()/sn()*hdivn() to recover the W gradz(T/S) trends
465
466         DO jk = 1, jpkm1
467            DO jj = 2, jpjm1
468               DO ji = fs_2, fs_jpim1   ! vector opt.
469#if defined key_zco
470                  zbtr      = 1. / ( e1t(ji,jj) * e2t(ji,jj) )
471                  z_hdivn_x = e2u(ji,jj)*pun(ji,jj,jk) - e2u(ji-1,jj)*pun(ji-1,jj,jk)
472                  z_hdivn_y = e1v(ji,jj)*pvn(ji,jj,jk) - e1v(ji,jj-1)*pvn(ji,jj-1,jk)
473#else
474                  zbtr      = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
475                  z_hdivn_x = e2u(ji,jj)*fse3u(ji,jj,jk)*pun(ji,jj,jk) - e2u(ji-1,jj)*fse3u(ji-1,jj,jk)*pun(ji-1,jj,jk)
476                  z_hdivn_y = e1v(ji,jj)*fse3v(ji,jj,jk)*pvn(ji,jj,jk) - e1v(ji,jj-1)*fse3v(ji,jj-1,jk)*pvn(ji,jj-1,jk)
477#endif
478                  z_hdivn   = (z_hdivn_x + z_hdivn_y) * zbtr
479                  ztrdt(ji,jj,jk) = ta(ji,jj,jk) - ztrdt(ji,jj,jk) - tn(ji,jj,jk) * z_hdivn 
480                  ztrds(ji,jj,jk) = sa(ji,jj,jk) - ztrds(ji,jj,jk) - sn(ji,jj,jk) * z_hdivn
481               END DO
482            END DO
483         END DO
484         CALL trd_mod(ztrdt, ztrds, jptra_trd_zad, 'TRA', kt)
485         !
486      ENDIF
487
488      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ta, clinfo1=' muscl2 zad  - Ta: ', mask1=tmask,   &
489         &                       tab3d_2=sa, clinfo2=              ' Sa: ', mask2=tmask, clinfo3='tra' )
490      !
491   END SUBROUTINE tra_adv_muscl2
492
493   !!======================================================================
494END MODULE traadv_muscl2
Note: See TracBrowser for help on using the repository browser.