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 @ 67

Last change on this file since 67 was 67, checked in by opalod, 20 years ago

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