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

Last change on this file since 216 was 216, checked in by opalod, 19 years ago

CT : UPDATE151 : New trends organization
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 20.9 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 trdmod          ! ocean active tracers trends
15   USE trdmod_oce      ! ocean variables trends
16   USE in_out_manager  ! I/O manager
17   USE dynspg_fsc      ! surface pressure gradient
18   USE dynspg_fsc_atsk ! autotasked 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
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 T and
44      !!      S using a MUSCL scheme (Monotone Upstream-centered Scheme for
45      !!      Conservation Laws) and add it to the general tracer trend.
46      !!
47      !! ** Method  : MUSCL scheme plus centered scheme at ocean boundaries
48      !!
49      !! ** Action  : - update (ta,sa) with the now advective tracer trends
50      !!              - save trends in (ztdta,ztdsa) ('key_trdtra')
51      !!
52      !! References :               
53      !!      Estubier, A., and M. Levy, Notes Techn. Pole de Modelisation
54      !! IPSL, Sept. 2000 (http://www.lodyc.jussieu.fr/opa)
55      !!
56      !! History :
57      !!        !  06-00  (A.Estublier)  for passive tracers
58      !!        !  01-08  (E.Durand G.Madec)  adapted for T & S
59      !!   8.5  !  02-06  (G. Madec)  F90: Free form and module
60      !!   9.0  !  04-08  (C. Talandier) New trends organization
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
73      !! * Arguments
74      INTEGER, INTENT( in ) ::   kt         ! ocean time-step
75
76      !! * Local declarations
77      INTEGER ::   ji, jj, jk               ! dummy loop indices
78      REAL(wp) ::   &
79         zu, zv, zw, zeu, zev,           & 
80         zew, zbtr, zstep,               &
81         z0u, z0v, z0w,                  &
82         zzt1, zzt2, zalpha,             &
83         zzs1, zzs2, z2,                 &
84         zta, zsa
85      REAL(wp), DIMENSION (jpi,jpj,jpk) ::   &
86         zt1, zt2, ztp1, ztp2,   &
87         zs1, zs2, zsp1, zsp2,   &
88         ztdta, ztdsa
89      !!----------------------------------------------------------------------
90
91      IF( kt == nit000 .AND. lwp ) THEN
92         WRITE(numout,*)
93         WRITE(numout,*) 'tra_adv_muscl2 : MUSCL2 advection scheme'
94         WRITE(numout,*) '~~~~~~~~~~~~~~~'
95      ENDIF
96
97      IF( neuler == 0 .AND. kt == nit000 ) THEN
98          z2=1.
99      ELSE
100          z2=2.
101      ENDIF
102
103      ! Save ta and sa trends
104      IF( l_trdtra )   THEN
105         ztdta(:,:,:) = ta(:,:,:) 
106         ztdsa(:,:,:) = sa(:,:,:) 
107         l_adv = 'mu2'
108      ENDIF
109
110#if defined key_trabbl_adv
111      ! Advective bottom boundary layer
112      ! -------------------------------
113      zun(:,:,:) = un (:,:,:) - u_bbl(:,:,:)
114      zvn(:,:,:) = vn (:,:,:) - v_bbl(:,:,:)
115      zwn(:,:,:) = wn (:,:,:) + w_bbl( :,:,:)
116#endif
117
118
119      ! I. Horizontal advective fluxes
120      ! ------------------------------
121
122      ! first guess of the slopes
123      ! interior values
124      DO jk = 1, jpkm1
125         DO jj = 1, jpjm1     
126            DO ji = 1, fs_jpim1   ! vector opt.
127               zt1(ji,jj,jk) = umask(ji,jj,jk) * ( tb(ji+1,jj,jk) - tb(ji,jj,jk) )
128               zs1(ji,jj,jk) = umask(ji,jj,jk) * ( sb(ji+1,jj,jk) - sb(ji,jj,jk) )
129               zt2(ji,jj,jk) = vmask(ji,jj,jk) * ( tb(ji,jj+1,jk) - tb(ji,jj,jk) )
130               zs2(ji,jj,jk) = vmask(ji,jj,jk) * ( sb(ji,jj+1,jk) - sb(ji,jj,jk) )
131            END DO
132         END DO
133      END DO
134      ! bottom values
135      zt1(:,:,jpk) = 0.e0    ;    zt2(:,:,jpk) = 0.e0
136      zs1(:,:,jpk) = 0.e0    ;    zs2(:,:,jpk) = 0.e0
137
138      ! lateral boundary conditions on zt1, zt2 ; zs1, zs2   (changed sign)
139      CALL lbc_lnk( zt1, 'U', -1. )   ;   CALL lbc_lnk( zs1, 'U', -1. )
140      CALL lbc_lnk( zt2, 'V', -1. )   ;   CALL lbc_lnk( zs2, 'V', -1. )
141
142      ! Slopes
143      ! interior values
144      DO jk = 1, jpkm1
145         DO jj = 2, jpj
146            DO ji = fs_2, jpi   ! vector opt.
147               ztp1(ji,jj,jk) =                    ( zt1(ji,jj,jk) + zt1(ji-1,jj  ,jk) )   &
148                  &           * ( 0.25 + SIGN( 0.25, zt1(ji,jj,jk) * zt1(ji-1,jj  ,jk) ) )
149               zsp1(ji,jj,jk) =                    ( zs1(ji,jj,jk) + zs1(ji-1,jj  ,jk) )   &
150                  &           * ( 0.25 + SIGN( 0.25, zs1(ji,jj,jk) * zs1(ji-1,jj  ,jk) ) )
151               ztp2(ji,jj,jk) =                    ( zt2(ji,jj,jk) + zt2(ji  ,jj-1,jk) )   &
152                  &           * ( 0.25 + SIGN( 0.25, zt2(ji,jj,jk) * zt2(ji  ,jj-1,jk) ) )
153               zsp2(ji,jj,jk) =                    ( zs2(ji,jj,jk) + zs2(ji  ,jj-1,jk) )   &
154                  &           * ( 0.25 + SIGN( 0.25, zs2(ji,jj,jk) * zs2(ji  ,jj-1,jk) ) )
155            END DO
156         END DO
157      END DO
158      ! bottom values
159      ztp1(:,:,jpk) = 0.e0    ;    ztp2(:,:,jpk) = 0.e0
160      zsp1(:,:,jpk) = 0.e0    ;    zsp2(:,:,jpk) = 0.e0
161
162      ! Slopes limitation
163      DO jk = 1, jpkm1
164         DO jj = 2, jpj
165            DO ji = fs_2, jpi   ! vector opt.
166               ztp1(ji,jj,jk) = SIGN( 1., ztp1(ji,jj,jk) )   &
167                  &           * MIN(    ABS( ztp1(ji  ,jj,jk) ),   &
168                  &                  2.*ABS( zt1 (ji-1,jj,jk) ),   &
169                  &                  2.*ABS( zt1 (ji  ,jj,jk) ) )
170               zsp1(ji,jj,jk) = SIGN( 1., zsp1(ji,jj,jk) )   &
171                  &           * MIN(    ABS( zsp1(ji  ,jj,jk) ),   &
172                  &                  2.*ABS( zs1 (ji-1,jj,jk) ),   &
173                  &                  2.*ABS( zs1 (ji  ,jj,jk) ) )
174               ztp2(ji,jj,jk) = SIGN( 1., ztp2(ji,jj,jk) )   &
175                  &           * MIN(    ABS( ztp2(ji,jj  ,jk) ),   &
176                  &                  2.*ABS( zt2 (ji,jj-1,jk) ),   &
177                  &                  2.*ABS( zt2 (ji,jj  ,jk) ) )
178               zsp2(ji,jj,jk) = SIGN( 1., zsp2(ji,jj,jk) )   &
179                  &           * MIN(    ABS( zsp2(ji,jj  ,jk) ),   &
180                  &                  2.*ABS( zs2 (ji,jj-1,jk) ),   &
181                  &                  2.*ABS( zs2 (ji,jj  ,jk) ) )
182            END DO
183         END DO
184      END DO       
185
186      ! Advection terms
187      ! interior values
188      DO jk = 1, jpkm1
189         zstep  = z2 * rdttra(jk)
190         DO jj = 2, jpjm1     
191            DO ji = fs_2, fs_jpim1   ! vector opt.
192               ! volume fluxes
193#if defined key_s_coord || defined key_partial_steps
194               zeu = e2u(ji,jj) * fse3u(ji,jj,jk) * zun(ji,jj,jk)
195               zev = e1v(ji,jj) * fse3v(ji,jj,jk) * zvn(ji,jj,jk)
196#else
197               zeu = e2u(ji,jj) * zun(ji,jj,jk)
198               zev = e1v(ji,jj) * zvn(ji,jj,jk)
199#endif
200               ! MUSCL fluxes
201               z0u = SIGN( 0.5, zun(ji,jj,jk) )           
202               zalpha = 0.5 - z0u
203               zu  = z0u - 0.5 * zun(ji,jj,jk) * zstep / e1u(ji,jj)
204               zzt1 = tb(ji+1,jj,jk) + zu*ztp1(ji+1,jj,jk)
205               zzt2 = tb(ji  ,jj,jk) + zu*ztp1(ji  ,jj,jk)
206               zzs1 = sb(ji+1,jj,jk) + zu*zsp1(ji+1,jj,jk)
207               zzs2 = sb(ji  ,jj,jk) + zu*zsp1(ji  ,jj,jk)
208               zt1(ji,jj,jk) = zeu * ( zalpha * zzt1 + (1.-zalpha) * zzt2 )
209               zs1(ji,jj,jk) = zeu * ( zalpha * zzs1 + (1.-zalpha) * zzs2 )
210
211               z0v = SIGN( 0.5, zvn(ji,jj,jk) )           
212               zalpha = 0.5 - z0v
213               zv  = z0v - 0.5 * zvn(ji,jj,jk) * zstep / e2v(ji,jj)
214               zzt1 = tb(ji,jj+1,jk) + zv*ztp2(ji,jj+1,jk)
215               zzt2 = tb(ji,jj  ,jk) + zv*ztp2(ji,jj  ,jk)
216               zzs1 = sb(ji,jj+1,jk) + zv*zsp2(ji,jj+1,jk)
217               zzs2 = sb(ji,jj  ,jk) + zv*zsp2(ji,jj  ,jk)
218               zt2(ji,jj,jk) = zev * ( zalpha * zzt1 + (1.-zalpha) * zzt2 )
219               zs2(ji,jj,jk) = zev * ( zalpha * zzs1 + (1.-zalpha) * zzs2 )
220            END DO
221         END DO
222      END DO
223
224      !!!!  centered scheme at lateral b.C. if off-shore velocity
225      DO jk = 1, jpkm1
226        DO jj = 2, jpjm1
227            DO ji = fs_2, fs_jpim1   ! vector opt.
228#if defined key_s_coord || defined key_partial_steps
229               zev = e1v(ji,jj) * fse3v(ji,jj,jk)
230               IF( umask(ji,jj,jk) == 0. ) THEN
231                  IF( zun(ji+1,jj,jk) > 0. .AND. ji /= jpi ) THEN
232                     zt1(ji+1,jj,jk) = e2u(ji+1,jj)* fse3u(ji+1,jj,jk)   &
233                        &            * zun(ji+1,jj,jk) * ( tb(ji+1,jj,jk) + tb(ji+2,jj,jk) ) * 0.5
234                     zs1(ji+1,jj,jk) = e2u(ji+1,jj)* fse3u(ji+1,jj,jk)   &
235                        &            * zun(ji+1,jj,jk) * ( sb(ji+1,jj,jk) + sb(ji+2,jj,jk) ) * 0.5
236                  ENDIF
237                  IF( zun(ji-1,jj,jk) < 0. ) THEN
238                     zt1(ji-1,jj,jk) = e2u(ji-1,jj)* fse3u(ji-1,jj,jk)   &
239                        &            * zun(ji-1,jj,jk) * ( tb(ji-1,jj,jk) + tb(ji  ,jj,jk) ) * 0.5
240                     zs1(ji-1,jj,jk) = e2u(ji-1,jj)* fse3u(ji-1,jj,jk)   &
241                        &            * zun(ji-1,jj,jk) * ( sb(ji-1,jj,jk) + sb(ji  ,jj,jk) ) * 0.5
242                  ENDIF
243               ENDIF
244               IF( vmask(ji,jj,jk) == 0. ) THEN
245                  IF( zvn(ji,jj+1,jk) > 0. .AND. jj /= jpj ) THEN
246                     zt2(ji,jj+1,jk) = e1v(ji,jj+1) * fse3v(ji,jj+1,jk)   &
247                        &            * zvn(ji,jj+1,jk) * ( tb(ji,jj+1,jk) + tb(ji,jj+2,jk) ) * 0.5
248                     zs2(ji,jj+1,jk) = e1v(ji,jj+1) * fse3v(ji,jj+1,jk)   &
249                        &            * zvn(ji,jj+1,jk) * ( sb(ji,jj+1,jk) + sb(ji,jj+2,jk) ) * 0.5
250                  ENDIF
251                  IF( zvn(ji,jj-1,jk) < 0. ) THEN
252                     zt2(ji,jj-1,jk) = e1v(ji,jj-1)* fse3v(ji,jj-1,jk)   &
253                        &            * zvn(ji,jj-1,jk) * ( tb(ji,jj-1,jk) + tb(ji  ,jj,jk) ) * 0.5
254                     zs2(ji,jj-1,jk) = e1v(ji,jj-1)* fse3v(ji,jj-1,jk)   &
255                        &            * zvn(ji,jj-1,jk) * ( sb(ji,jj-1,jk) + sb(ji  ,jj,jk) ) * 0.5
256                  ENDIF
257               ENDIF
258
259#else
260               IF( umask(ji,jj,jk) == 0. ) THEN
261                  IF( zun(ji+1,jj,jk) > 0. .AND. ji /= jpi ) THEN
262                     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
263                     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
264                  ENDIF
265                  IF( zun(ji-1,jj,jk) < 0. ) THEN
266                     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
267                     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
268                  ENDIF
269               ENDIF
270               IF( vmask(ji,jj,jk) == 0. ) THEN
271                  IF( zvn(ji,jj+1,jk) > 0. .AND. jj /= jpj ) THEN
272                     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
273                     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
274                  ENDIF
275                  IF( zvn(ji,jj-1,jk) < 0. ) THEN
276                     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
277                     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
278                  ENDIF
279               ENDIF
280#endif
281            END DO
282         END DO
283      END DO
284
285      ! lateral boundary conditions on zt1, zt2 ; zs1, zs2   (changed sign)
286      CALL lbc_lnk( zt1, 'U', -1. )   ;   CALL lbc_lnk( zs1, 'U', -1. ) 
287      CALL lbc_lnk( zt2, 'V', -1. )   ;   CALL lbc_lnk( zs2, 'V', -1. )
288
289      ! Save MUSCL fluxes to compute i- & j- horizontal
290      ! advection trends in the MLD
291      IF( l_trdtra )   THEN
292         ! save i- terms
293         tladi(:,:,:) = zt1(:,:,:) 
294         sladi(:,:,:) = zs1(:,:,:) 
295         ! save j- terms
296         tladj(:,:,:) = zt2(:,:,:) 
297         sladj(:,:,:) = zs2(:,:,:) 
298      ENDIF
299
300      ! Compute & add the horizontal advective trend
301
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            END DO
319        END DO
320      END DO       
321
322      ! Save the horizontal advective trends for diagnostic
323
324      IF( l_trdtra )   THEN
325         ! Recompute the hoizontal advection zta & zsa trends computed
326         ! at the step 2. above in making the difference between the new
327         ! trends and the previous one ta()/sa - ztdta()/ztdsa() and add
328         ! the term tn()/sn()*hdivn() to recover the Uh gradh(T/S) trends
329         ztdta(:,:,:) = ta(:,:,:) - ztdta(:,:,:) + tn(:,:,:) * hdivn(:,:,:) 
330         ztdsa(:,:,:) = sa(:,:,:) - ztdsa(:,:,:) + sn(:,:,:) * hdivn(:,:,:)
331
332         CALL trd_mod(ztdta, ztdsa, jpttdlad, 'TRA', kt)
333
334         ! Save the new ta and sa trends
335         ztdta(:,:,:) = ta(:,:,:) 
336         ztdsa(:,:,:) = sa(:,:,:) 
337
338      ENDIF
339
340      IF(l_ctl) THEN
341         zta = SUM( ta(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
342         zsa = SUM( sa(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
343         WRITE(numout,*) ' had  - Ta: ', zta-t_ctl, ' Sa: ', zsa-s_ctl, ' muscl'
344         t_ctl = zta   ;   s_ctl = zsa
345      ENDIF
346
347      ! "zonal" mean advective heat and salt transport
348      IF( ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN
349# if defined key_s_coord || defined key_partial_steps
350         pht_adv(:) = ptr_vj( zt2(:,:,:) )
351         pst_adv(:) = ptr_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(:) = ptr_vj( zt2(:,:,:) )
362         pst_adv(:) = ptr_vj( zs2(:,:,:) )
363# endif
364      ENDIF
365
366      ! II. Vertical advective fluxes
367      ! -----------------------------
368     
369      ! First guess of the slope
370      ! interior values
371      DO jk = 2, jpkm1
372         zt1(:,:,jk) = tmask(:,:,jk) * ( tb(:,:,jk-1) - tb(:,:,jk) )
373         zs1(:,:,jk) = tmask(:,:,jk) * ( sb(:,:,jk-1) - sb(:,:,jk) )
374      END DO
375      ! surface & bottom boundary conditions
376      zt1 (:,:, 1 ) = 0.e0    ;    zt1 (:,:,jpk) = 0.e0
377      zs1 (:,:, 1 ) = 0.e0    ;    zs1 (:,:,jpk) = 0.e0
378
379      ! Slopes
380      DO jk = 2, jpkm1
381         DO jj = 1, jpj
382            DO ji = 1, jpi
383               ztp1(ji,jj,jk) =                    ( zt1(ji,jj,jk) + zt1(ji,jj,jk+1) )   &
384                  &           * ( 0.25 + SIGN( 0.25, zt1(ji,jj,jk) * zt1(ji,jj,jk+1) ) )
385               zsp1(ji,jj,jk) =                    ( zs1(ji,jj,jk) + zs1(ji,jj,jk+1) )   &
386                  &           * ( 0.25 + SIGN( 0.25, zs1(ji,jj,jk) * zs1(ji,jj,jk+1) ) )
387            END DO
388         END DO
389      END DO
390
391      ! Slopes limitation
392      ! interior values
393      DO jk = 2, jpkm1
394         DO jj = 1, jpj
395            DO ji = 1, jpi
396               ztp1(ji,jj,jk) = SIGN( 1., ztp1(ji,jj,jk) )   &
397                  &           * MIN(    ABS( ztp1(ji,jj,jk  ) ),   &
398                  &                  2.*ABS( zt1 (ji,jj,jk+1) ),   &
399                  &                  2.*ABS( zt1 (ji,jj,jk  ) ) )
400               zsp1(ji,jj,jk) = SIGN( 1., zsp1(ji,jj,jk) )   &
401                  &           * MIN(    ABS( zsp1(ji,jj,jk  ) ),   &
402                  &                  2.*ABS( zs1 (ji,jj,jk+1) ),   &
403                  &                  2.*ABS( zs1 (ji,jj,jk  ) ) )
404            END DO
405         END DO
406      END DO
407      ! surface values
408      ztp1(:,:,1) = 0.e0
409      zsp1(:,:,1) = 0.e0
410
411      ! vertical advective flux
412      ! interior values
413      DO jk = 1, jpkm1
414         zstep  = z2 * rdttra(jk)
415         DO jj = 2, jpjm1     
416            DO ji = fs_2, fs_jpim1   ! vector opt.
417               zew = zwn(ji,jj,jk+1)
418               z0w = SIGN( 0.5, zwn(ji,jj,jk+1) )
419               zalpha = 0.5 + z0w
420               zw  = z0w - 0.5 * zwn(ji,jj,jk+1)*zstep / fse3w(ji,jj,jk+1)
421               zzt1 = tb(ji,jj,jk+1) + zw*ztp1(ji,jj,jk+1)
422               zzt2 = tb(ji,jj,jk  ) + zw*ztp1(ji,jj,jk  )
423               zzs1 = sb(ji,jj,jk+1) + zw*zsp1(ji,jj,jk+1)
424               zzs2 = sb(ji,jj,jk  ) + zw*zsp1(ji,jj,jk  )
425               zt1(ji,jj,jk+1) = zew * ( zalpha * zzt1 + (1.-zalpha)*zzt2 )
426               zs1(ji,jj,jk+1) = zew * ( zalpha * zzs1 + (1.-zalpha)*zzs2 )
427            END DO
428         END DO
429      END DO
430      DO jk = 2, jpkm1
431        DO jj = 2, jpjm1
432            DO ji = fs_2, fs_jpim1   ! vector opt.
433               IF( tmask(ji,jj,jk+1) == 0. ) THEN
434                  IF( zwn(ji,jj,jk) > 0. ) THEN
435                     zt1(ji,jj,jk) = zwn(ji,jj,jk) * ( tb(ji,jj,jk-1) + tb(ji,jj,jk) ) * 0.5
436                     zs1(ji,jj,jk) = zwn(ji,jj,jk) * ( sb(ji,jj,jk-1) + sb(ji,jj,jk) ) * 0.5
437                  ENDIF
438               ENDIF
439            END DO
440         END DO
441      END DO
442
443      ! surface values
444      IF( lk_dynspg_fsc .OR. lk_dynspg_fsc_tsk ) THEN   ! free surface-constant volume
445         zt1(:,:, 1 ) = zwn(:,:,1) * tb(:,:,1)
446         zs1(:,:, 1 ) = zwn(:,:,1) * sb(:,:,1)
447      ELSE                                              ! rigid lid : flux set to zero
448         zt1(:,:, 1 ) = 0.e0
449         zs1(:,:, 1 ) = 0.e0
450      ENDIF
451
452      ! bottom values
453      zt1(:,:,jpk) = 0.e0
454      zs1(:,:,jpk) = 0.e0
455
456
457      ! Compute & add the vertical advective trend
458
459      DO jk = 1, jpkm1
460         DO jj = 2, jpjm1     
461            DO ji = fs_2, fs_jpim1   ! vector opt.
462               zbtr = 1. / fse3t(ji,jj,jk)
463               ! horizontal advective trends
464               zta = - zbtr * ( zt1(ji,jj,jk) - zt1(ji,jj,jk+1) )
465               zsa = - zbtr * ( zs1(ji,jj,jk) - zs1(ji,jj,jk+1) )
466               ! add it to the general tracer trends
467               ta(ji,jj,jk) =  ta(ji,jj,jk) + zta
468               sa(ji,jj,jk) =  sa(ji,jj,jk) + zsa
469            END DO
470         END DO
471      END DO
472
473      ! Save the vertical advective trends for diagnostic
474
475      IF( l_trdtra )   THEN
476         ! Recompute the vertical advection zta & zsa trends computed
477         ! at the step 2. above in making the difference between the new
478         ! trends and the previous one: ta()/sa - ztdta()/ztdsa() and substract
479         ! the term tn()/sn()*hdivn() to recover the W gradz(T/S) trends
480         ztdta(:,:,:) = ta(:,:,:) - ztdta(:,:,:) - tn(:,:,:) * hdivn(:,:,:) 
481         ztdsa(:,:,:) = sa(:,:,:) - ztdsa(:,:,:) - sn(:,:,:) * hdivn(:,:,:)
482
483         CALL trd_mod(ztdta, ztdsa, jpttdzad, 'TRA', kt)
484      ENDIF
485
486      IF(l_ctl) THEN
487         zta = SUM( ta(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
488         zsa = SUM( sa(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
489         WRITE(numout,*) ' zad  - Ta: ', zta-t_ctl, ' Sa: ', zsa-s_ctl, ' muscl2'
490         t_ctl = zta   ;   s_ctl = zsa
491      ENDIF
492
493   END SUBROUTINE tra_adv_muscl2
494
495   !!======================================================================
496END MODULE traadv_muscl2
Note: See TracBrowser for help on using the repository browser.