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

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

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