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

source: tags/nemo_dev_x6/NEMO/OPA_SRC/TRA/traadv_muscl2.F90 @ 6934

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

CT : UPDATE082 : Finalization of the poleward transport diagnostics
  • 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          ! ocean lateral boundary condition (or mpp link)
21   USE diaptr          ! 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) ::   zzt1, zzt2, zalpha
83      REAL(wp) ::   zzs1, zzs2
84      REAL(wp) ::   z2
85#if defined key_trdtra || defined key_trdmld
86      REAL(wp) ::   ztai, ztaj, zsai, zsaj
87      REAL(wp) ::   zfui, zfvj
88#endif
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#if defined key_trabbl_adv
104
105      ! Advective bottom boundary layer
106      ! -------------------------------
107      zun(:,:,:) = un (:,:,:) - u_bbl(:,:,:)
108      zvn(:,:,:) = vn (:,:,:) - v_bbl(:,:,:)
109      zwn(:,:,:) = wn (:,:,:) + w_bbl( :,:,:)
110#endif
111
112
113      ! I. Horizontal advective fluxes
114      ! ------------------------------
115
116      ! first guess of the slopes
117      ! interior values
118      DO jk = 1, jpkm1
119         DO jj = 1, jpjm1     
120            DO ji = 1, fs_jpim1   ! vector opt.
121               zt1(ji,jj,jk) = umask(ji,jj,jk) * ( tb(ji+1,jj,jk) - tb(ji,jj,jk) )
122               zs1(ji,jj,jk) = umask(ji,jj,jk) * ( sb(ji+1,jj,jk) - sb(ji,jj,jk) )
123               zt2(ji,jj,jk) = vmask(ji,jj,jk) * ( tb(ji,jj+1,jk) - tb(ji,jj,jk) )
124               zs2(ji,jj,jk) = vmask(ji,jj,jk) * ( sb(ji,jj+1,jk) - sb(ji,jj,jk) )
125            END DO
126         END DO
127      END DO
128      ! bottom values
129      zt1(:,:,jpk) = 0.e0    ;    zt2(:,:,jpk) = 0.e0
130      zs1(:,:,jpk) = 0.e0    ;    zs2(:,:,jpk) = 0.e0
131
132      ! lateral boundary conditions on zt1, zt2 ; zs1, zs2   (changed sign)
133      CALL lbc_lnk( zt1, 'U', -1. )   ;   CALL lbc_lnk( zs1, 'U', -1. )
134      CALL lbc_lnk( zt2, 'V', -1. )   ;   CALL lbc_lnk( zs2, 'V', -1. )
135
136      ! Slopes
137      ! interior values
138      DO jk = 1, jpkm1
139         DO jj = 2, jpj
140            DO ji = fs_2, jpi   ! vector opt.
141               ztp1(ji,jj,jk) =                    ( zt1(ji,jj,jk) + zt1(ji-1,jj  ,jk) )   &
142                  &           * ( 0.25 + SIGN( 0.25, zt1(ji,jj,jk) * zt1(ji-1,jj  ,jk) ) )
143               zsp1(ji,jj,jk) =                    ( zs1(ji,jj,jk) + zs1(ji-1,jj  ,jk) )   &
144                  &           * ( 0.25 + SIGN( 0.25, zs1(ji,jj,jk) * zs1(ji-1,jj  ,jk) ) )
145               ztp2(ji,jj,jk) =                    ( zt2(ji,jj,jk) + zt2(ji  ,jj-1,jk) )   &
146                  &           * ( 0.25 + SIGN( 0.25, zt2(ji,jj,jk) * zt2(ji  ,jj-1,jk) ) )
147               zsp2(ji,jj,jk) =                    ( zs2(ji,jj,jk) + zs2(ji  ,jj-1,jk) )   &
148                  &           * ( 0.25 + SIGN( 0.25, zs2(ji,jj,jk) * zs2(ji  ,jj-1,jk) ) )
149            END DO
150         END DO
151      END DO
152      ! bottom values
153      ztp1(:,:,jpk) = 0.e0    ;    ztp2(:,:,jpk) = 0.e0
154      zsp1(:,:,jpk) = 0.e0    ;    zsp2(:,:,jpk) = 0.e0
155
156      ! Slopes limitation
157      DO jk = 1, jpkm1
158         DO jj = 2, jpj
159            DO ji = fs_2, jpi   ! vector opt.
160               ztp1(ji,jj,jk) = SIGN( 1., ztp1(ji,jj,jk) )   &
161                  &           * MIN(    ABS( ztp1(ji  ,jj,jk) ),   &
162                  &                  2.*ABS( zt1 (ji-1,jj,jk) ),   &
163                  &                  2.*ABS( zt1 (ji  ,jj,jk) ) )
164               zsp1(ji,jj,jk) = SIGN( 1., zsp1(ji,jj,jk) )   &
165                  &           * MIN(    ABS( zsp1(ji  ,jj,jk) ),   &
166                  &                  2.*ABS( zs1 (ji-1,jj,jk) ),   &
167                  &                  2.*ABS( zs1 (ji  ,jj,jk) ) )
168               ztp2(ji,jj,jk) = SIGN( 1., ztp2(ji,jj,jk) )   &
169                  &           * MIN(    ABS( ztp2(ji,jj  ,jk) ),   &
170                  &                  2.*ABS( zt2 (ji,jj-1,jk) ),   &
171                  &                  2.*ABS( zt2 (ji,jj  ,jk) ) )
172               zsp2(ji,jj,jk) = SIGN( 1., zsp2(ji,jj,jk) )   &
173                  &           * MIN(    ABS( zsp2(ji,jj  ,jk) ),   &
174                  &                  2.*ABS( zs2 (ji,jj-1,jk) ),   &
175                  &                  2.*ABS( zs2 (ji,jj  ,jk) ) )
176            END DO
177         END DO
178      END DO       
179
180      ! Advection terms
181      ! interior values
182      DO jk = 1, jpkm1
183         zstep  = z2 * rdttra(jk)
184         DO jj = 2, jpjm1     
185            DO ji = fs_2, fs_jpim1   ! vector opt.
186               ! volume fluxes
187#if defined key_s_coord || defined key_partial_steps
188               zeu = e2u(ji,jj) * fse3u(ji,jj,jk) * zun(ji,jj,jk)
189               zev = e1v(ji,jj) * fse3v(ji,jj,jk) * zvn(ji,jj,jk)
190#else
191               zeu = e2u(ji,jj) * zun(ji,jj,jk)
192               zev = e1v(ji,jj) * zvn(ji,jj,jk)
193#endif
194               ! MUSCL fluxes
195               z0u = SIGN( 0.5, zun(ji,jj,jk) )           
196               zalpha = 0.5 - z0u
197               zu  = z0u - 0.5 * zun(ji,jj,jk) * zstep / e1u(ji,jj)
198               zzt1 = tb(ji+1,jj,jk) + zu*ztp1(ji+1,jj,jk)
199               zzt2 = tb(ji  ,jj,jk) + zu*ztp1(ji  ,jj,jk)
200               zzs1 = sb(ji+1,jj,jk) + zu*zsp1(ji+1,jj,jk)
201               zzs2 = sb(ji  ,jj,jk) + zu*zsp1(ji  ,jj,jk)
202               zt1(ji,jj,jk) = zeu * ( zalpha * zzt1 + (1.-zalpha) * zzt2 )
203               zs1(ji,jj,jk) = zeu * ( zalpha * zzs1 + (1.-zalpha) * zzs2 )
204
205               z0v = SIGN( 0.5, zvn(ji,jj,jk) )           
206               zalpha = 0.5 - z0v
207               zv  = z0v - 0.5 * zvn(ji,jj,jk) * zstep / e2v(ji,jj)
208               zzt1 = tb(ji,jj+1,jk) + zv*ztp2(ji,jj+1,jk)
209               zzt2 = tb(ji,jj  ,jk) + zv*ztp2(ji,jj  ,jk)
210               zzs1 = sb(ji,jj+1,jk) + zv*zsp2(ji,jj+1,jk)
211               zzs2 = sb(ji,jj  ,jk) + zv*zsp2(ji,jj  ,jk)
212               zt2(ji,jj,jk) = zev * ( zalpha * zzt1 + (1.-zalpha) * zzt2 )
213               zs2(ji,jj,jk) = zev * ( zalpha * zzs1 + (1.-zalpha) * zzs2 )
214            END DO
215         END DO
216      END DO
217
218      !!!!  centered scheme at lateral b.C. if off-shore velocity
219      DO jk = 1, jpkm1
220        DO jj = 2, jpjm1
221            DO ji = fs_2, fs_jpim1   ! vector opt.
222#if defined key_s_coord || defined key_partial_steps
223               zev = e1v(ji,jj) * fse3v(ji,jj,jk)
224               IF( umask(ji,jj,jk) == 0. ) THEN
225                  IF( zun(ji+1,jj,jk) > 0. .AND. ji /= jpi ) THEN
226                     zt1(ji+1,jj,jk) = e2u(ji+1,jj)* fse3u(ji+1,jj,jk)   &
227                        &            * zun(ji+1,jj,jk) * ( tb(ji+1,jj,jk) + tb(ji+2,jj,jk) ) * 0.5
228                     zs1(ji+1,jj,jk) = e2u(ji+1,jj)* fse3u(ji+1,jj,jk)   &
229                        &            * zun(ji+1,jj,jk) * ( sb(ji+1,jj,jk) + sb(ji+2,jj,jk) ) * 0.5
230                  ENDIF
231                  IF( zun(ji-1,jj,jk) < 0. ) 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  ,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  ,jj,jk) ) * 0.5
236                  ENDIF
237               ENDIF
238               IF( vmask(ji,jj,jk) == 0. ) THEN
239                  IF( zvn(ji,jj+1,jk) > 0. .AND. jj /= jpj ) THEN
240                     zt2(ji,jj+1,jk) = e1v(ji,jj+1) * fse3v(ji,jj+1,jk)   &
241                        &            * zvn(ji,jj+1,jk) * ( tb(ji,jj+1,jk) + tb(ji,jj+2,jk) ) * 0.5
242                     zs2(ji,jj+1,jk) = e1v(ji,jj+1) * fse3v(ji,jj+1,jk)   &
243                        &            * zvn(ji,jj+1,jk) * ( sb(ji,jj+1,jk) + sb(ji,jj+2,jk) ) * 0.5
244                  ENDIF
245                  IF( zvn(ji,jj-1,jk) < 0. ) 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,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,jk) ) * 0.5
250                  ENDIF
251               ENDIF
252
253#else
254               IF( umask(ji,jj,jk) == 0. ) THEN
255                  IF( zun(ji+1,jj,jk) > 0. .AND. ji /= jpi ) THEN
256                     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
257                     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
258                  ENDIF
259                  IF( zun(ji-1,jj,jk) < 0. ) THEN
260                     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
261                     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
262                  ENDIF
263               ENDIF
264               IF( vmask(ji,jj,jk) == 0. ) THEN
265                  IF( zvn(ji,jj+1,jk) > 0. .AND. jj /= jpj ) THEN
266                     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
267                     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
268                  ENDIF
269                  IF( zvn(ji,jj-1,jk) < 0. ) THEN
270                     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
271                     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
272                  ENDIF
273               ENDIF
274#endif
275            END DO
276         END DO
277      END DO
278
279      ! lateral boundary conditions on zt1, zt2 ; zs1, zs2   (changed sign)
280      CALL lbc_lnk( zt1, 'U', -1. )   ;   CALL lbc_lnk( zs1, 'U', -1. ) 
281      CALL lbc_lnk( zt2, 'V', -1. )   ;   CALL lbc_lnk( zs2, 'V', -1. )
282
283      ! Compute & add the horizontal advective trend
284      DO jk = 1, jpkm1
285         DO jj = 2, jpjm1     
286            DO ji = fs_2, fs_jpim1   ! vector opt.
287#if defined key_s_coord || defined key_partial_steps
288               zbtr = 1. / ( e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk) )
289#else
290               zbtr = 1. / ( e1t(ji,jj)*e2t(ji,jj) )
291#endif
292               ! horizontal advective trends
293               zta = - zbtr * ( zt1(ji,jj,jk) - zt1(ji-1,jj  ,jk  )   &
294                  &           + zt2(ji,jj,jk) - zt2(ji  ,jj-1,jk  ) )
295               zsa = - zbtr * ( zs1(ji,jj,jk) - zs1(ji-1,jj  ,jk  )   &
296                  &           + zs2(ji,jj,jk) - zs2(ji  ,jj-1,jk  ) ) 
297               ! add it to the general tracer trends
298               ta(ji,jj,jk) = ta(ji,jj,jk) + zta
299               sa(ji,jj,jk) = sa(ji,jj,jk) + zsa
300#if defined key_trdtra || defined key_trdmld
301               ! save the horizontal advective trend of tracer
302               ttrd(ji,jj,jk,1) = zta + tn(ji,jj,jk) * hdivn(ji,jj,jk)
303               strd(ji,jj,jk,1) = zsa + sn(ji,jj,jk) * hdivn(ji,jj,jk)
304               ! recompute the trends in i- and j-direction as Uh gradh(T)
305#   if defined key_s_coord || defined key_partial_steps
306               zfui =  e2u(ji  ,jj) * fse3u(ji,  jj,jk) * un(ji,  jj,jk)   &
307                  & -  e2u(ji-1,jj) * fse3u(ji-1,jj,jk) * un(ji-1,jj,jk)
308               zfvj =  e1v(ji,jj  ) * fse3v(ji,jj  ,jk) * vn(ji,jj  ,jk)   &
309                  & -  e1v(ji,jj-1) * fse3v(ji,jj-1,jk) * vn(ji,jj-1,jk)
310#   else
311               zfui = e2u(ji  ,jj) * un(ji,  jj,jk)   &
312                  & - e2u(ji-1,jj) * un(ji-1,jj,jk)
313               zfvj = e1v(ji,jj  ) * vn(ji,jj  ,jk)   &
314                  & - e1v(ji,jj-1) * vn(ji,jj-1,jk)
315#   endif
316               ztai =-zbtr * (  zt1(ji,jj,jk) - zt1(ji-1,jj  ,jk) - tn(ji,jj,jk) * zfui  )
317               ztaj =-zbtr * (  zt2(ji,jj,jk) - zt2(ji  ,jj-1,jk) - tn(ji,jj,jk) * zfvj  )
318               zsai =-zbtr * (  zs1(ji,jj,jk) - zs1(ji-1,jj  ,jk) - sn(ji,jj,jk) * zfui  )
319               zsaj =-zbtr * (  zs2(ji,jj,jk) - zs2(ji  ,jj-1,jk) - sn(ji,jj,jk) * zfvj  )
320               ! save i- and j- advective trends computed as Uh gradh(T)
321               ttrdh(ji,jj,jk,1) = ztai
322               ttrdh(ji,jj,jk,2) = ztaj
323               strdh(ji,jj,jk,1) = zsai
324               strdh(ji,jj,jk,2) = zsaj
325#endif
326            END DO
327        END DO
328      END DO       
329
330      IF(l_ctl) THEN
331         zta = SUM( ta(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
332         zsa = SUM( sa(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
333         WRITE(numout,*) ' had  - Ta: ', zta-t_ctl, ' Sa: ', zsa-s_ctl, ' muscl'
334         t_ctl = zta   ;   s_ctl = zsa
335      ENDIF
336
337      ! "zonal" mean advective heat and salt transport
338      IF( ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN
339# if defined key_s_coord || defined key_partial_steps
340         pht_adv(:) = ptr_vj( zt2(:,:,:) )
341         pst_adv(:) = ptr_vj( zs2(:,:,:) )
342# else
343         DO jk = 1, jpkm1
344            DO jj = 2, jpjm1
345               DO ji = fs_2, fs_jpim1   ! vector opt.
346                 zt2(ji,jj,jk) = zt2(ji,jj,jk) * fse3v(ji,jj,jk)
347                 zs2(ji,jj,jk) = zs2(ji,jj,jk) * fse3v(ji,jj,jk)
348               END DO
349            END DO
350         END DO
351         pht_adv(:) = ptr_vj( zt2(:,:,:) )
352         pst_adv(:) = ptr_vj( zs2(:,:,:) )
353# endif
354      ENDIF
355
356      ! II. Vertical advective fluxes
357      ! -----------------------------
358     
359      ! First guess of the slope
360      ! interior values
361      DO jk = 2, jpkm1
362         zt1(:,:,jk) = tmask(:,:,jk) * ( tb(:,:,jk-1) - tb(:,:,jk) )
363         zs1(:,:,jk) = tmask(:,:,jk) * ( sb(:,:,jk-1) - sb(:,:,jk) )
364      END DO
365      ! surface & bottom boundary conditions
366      zt1 (:,:, 1 ) = 0.e0    ;    zt1 (:,:,jpk) = 0.e0
367      zs1 (:,:, 1 ) = 0.e0    ;    zs1 (:,:,jpk) = 0.e0
368
369      ! Slopes
370      DO jk = 2, jpkm1
371         DO jj = 1, jpj
372            DO ji = 1, jpi
373               ztp1(ji,jj,jk) =                    ( zt1(ji,jj,jk) + zt1(ji,jj,jk+1) )   &
374                  &           * ( 0.25 + SIGN( 0.25, zt1(ji,jj,jk) * zt1(ji,jj,jk+1) ) )
375               zsp1(ji,jj,jk) =                    ( zs1(ji,jj,jk) + zs1(ji,jj,jk+1) )   &
376                  &           * ( 0.25 + SIGN( 0.25, zs1(ji,jj,jk) * zs1(ji,jj,jk+1) ) )
377            END DO
378         END DO
379      END DO
380
381      ! Slopes limitation
382      ! interior values
383      DO jk = 2, jpkm1
384         DO jj = 1, jpj
385            DO ji = 1, jpi
386               ztp1(ji,jj,jk) = SIGN( 1., ztp1(ji,jj,jk) )   &
387                  &           * MIN(    ABS( ztp1(ji,jj,jk  ) ),   &
388                  &                  2.*ABS( zt1 (ji,jj,jk+1) ),   &
389                  &                  2.*ABS( zt1 (ji,jj,jk  ) ) )
390               zsp1(ji,jj,jk) = SIGN( 1., zsp1(ji,jj,jk) )   &
391                  &           * MIN(    ABS( zsp1(ji,jj,jk  ) ),   &
392                  &                  2.*ABS( zs1 (ji,jj,jk+1) ),   &
393                  &                  2.*ABS( zs1 (ji,jj,jk  ) ) )
394            END DO
395         END DO
396      END DO
397      ! surface values
398      ztp1(:,:,1) = 0.e0
399      zsp1(:,:,1) = 0.e0
400
401      ! vertical advective flux
402      ! interior values
403      DO jk = 1, jpkm1
404         zstep  = z2 * rdttra(jk)
405         DO jj = 2, jpjm1     
406            DO ji = fs_2, fs_jpim1   ! vector opt.
407               zew = zwn(ji,jj,jk+1)
408               z0w = SIGN( 0.5, zwn(ji,jj,jk+1) )
409               zalpha = 0.5 + z0w
410               zw  = z0w - 0.5 * zwn(ji,jj,jk+1)*zstep / fse3w(ji,jj,jk+1)
411               zzt1 = tb(ji,jj,jk+1) + zw*ztp1(ji,jj,jk+1)
412               zzt2 = tb(ji,jj,jk  ) + zw*ztp1(ji,jj,jk  )
413               zzs1 = sb(ji,jj,jk+1) + zw*zsp1(ji,jj,jk+1)
414               zzs2 = sb(ji,jj,jk  ) + zw*zsp1(ji,jj,jk  )
415               zt1(ji,jj,jk+1) = zew * ( zalpha * zzt1 + (1.-zalpha)*zzt2 )
416               zs1(ji,jj,jk+1) = zew * ( zalpha * zzs1 + (1.-zalpha)*zzs2 )
417            END DO
418         END DO
419      END DO
420      DO jk = 2, jpkm1
421        DO jj = 2, jpjm1
422            DO ji = fs_2, fs_jpim1   ! vector opt.
423               IF( tmask(ji,jj,jk+1) == 0. ) THEN
424                  IF( zwn(ji,jj,jk) > 0. ) THEN
425                     zt1(ji,jj,jk) = zwn(ji,jj,jk) * ( tb(ji,jj,jk-1) + tb(ji,jj,jk) ) * 0.5
426                     zs1(ji,jj,jk) = zwn(ji,jj,jk) * ( sb(ji,jj,jk-1) + sb(ji,jj,jk) ) * 0.5
427                  ENDIF
428               ENDIF
429            END DO
430         END DO
431      END DO
432
433      ! surface values
434      IF( lk_dynspg_fsc .OR. lk_dynspg_fsc_tsk ) THEN   ! free surface-constant volume
435         zt1(:,:, 1 ) = zwn(:,:,1) * tb(:,:,1)
436         zs1(:,:, 1 ) = zwn(:,:,1) * sb(:,:,1)
437      ELSE                                      ! rigid lid : flux set to zero
438         zt1(:,:, 1 ) = 0.e0
439         zs1(:,:, 1 ) = 0.e0
440      ENDIF
441
442      ! bottom values
443      zt1(:,:,jpk) = 0.e0
444      zs1(:,:,jpk) = 0.e0
445
446
447      ! Compute & add the vertical advective trend
448
449      DO jk = 1, jpkm1
450         DO jj = 2, jpjm1     
451            DO ji = fs_2, fs_jpim1   ! vector opt.
452               zbtr = 1. / fse3t(ji,jj,jk)
453               ! horizontal advective trends
454               zta = - zbtr * ( zt1(ji,jj,jk) - zt1(ji,jj,jk+1) )
455               zsa = - zbtr * ( zs1(ji,jj,jk) - zs1(ji,jj,jk+1) )
456               ! add it to the general tracer trends
457               ta(ji,jj,jk) =  ta(ji,jj,jk) + zta
458               sa(ji,jj,jk) =  sa(ji,jj,jk) + zsa
459#if defined key_trdtra || defined key_trdmld
460               ! save the vertical advective trends computed as w gradz(T)
461               ttrd(ji,jj,jk,2) = zta - tn(ji,jj,jk) * hdivn(ji,jj,jk)
462               strd(ji,jj,jk,2) = zsa - sn(ji,jj,jk) * hdivn(ji,jj,jk)
463#endif
464            END DO
465         END DO
466      END DO
467
468      IF(l_ctl) THEN
469         zta = SUM( ta(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
470         zsa = SUM( sa(2:jpim1,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
471         WRITE(numout,*) ' zad  - Ta: ', zta-t_ctl, ' Sa: ', zsa-s_ctl, ' muscl2'
472         t_ctl = zta   ;   s_ctl = zsa
473      ENDIF
474
475   END SUBROUTINE tra_adv_muscl2
476
477   !!======================================================================
478END MODULE traadv_muscl2
Note: See TracBrowser for help on using the repository browser.