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

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

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