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

Last change on this file since 467 was 457, checked in by opalod, 18 years ago

nemo_v1_update_049:RB: reorganization of tracers part, remove traadv_cen2_atsk.h90 traldf_iso_zps.F90 trazdf_iso.F90 trazdf_iso_vopt.F90, change atsk routines to jki
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 20.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_oce      ! choice/control of key cpp for surface pressure gradient
18   USE trabbl          ! tracers: bottom boundary layer
19   USE lib_mpp
20   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
21   USE diaptr          ! poleward transport diagnostics
22   USE prtctl          ! Print control
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, pun, pvn, pwn )
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      !! * Arguments
65      INTEGER , INTENT( in ) ::   kt          ! ocean time-step index
66      REAL(wp), INTENT( in ), DIMENSION(jpi,jpj,jpk) ::   &
67         pun, pvn, pwn                        ! now ocean velocity fields
68
69
70      !! * Local declarations
71      INTEGER ::   ji, jj, jk               ! dummy loop indices
72      REAL(wp) ::   &
73         zu, zv, zw, zeu, zev,           & 
74         zew, zbtr, zstep,               &
75         z0u, z0v, z0w,                  &
76         zzt1, zzt2, zalpha,             &
77         zzs1, zzs2, z2,                 &
78         zta, zsa
79      REAL(wp), DIMENSION (jpi,jpj,jpk) ::   &
80         zt1, zt2, ztp1, ztp2,   &
81         zs1, zs2, zsp1, zsp2,   &
82         ztdta, ztdsa
83      !!----------------------------------------------------------------------
84
85      IF( kt == nit000 .AND. lwp ) THEN
86         WRITE(numout,*)
87         WRITE(numout,*) 'tra_adv_muscl2 : MUSCL2 advection scheme'
88         WRITE(numout,*) '~~~~~~~~~~~~~~~'
89      ENDIF
90
91      IF( neuler == 0 .AND. kt == nit000 ) THEN
92          z2=1.
93      ELSE
94          z2=2.
95      ENDIF
96
97      ! Save ta and sa trends
98      IF( l_trdtra )   THEN
99         ztdta(:,:,:) = ta(:,:,:) 
100         ztdsa(:,:,:) = sa(:,:,:) 
101         l_adv = 'mu2'
102      ENDIF
103
104
105      ! I. Horizontal advective fluxes
106      ! ------------------------------
107
108      ! first guess of the slopes
109      ! interior values
110      DO jk = 1, jpkm1
111         DO jj = 1, jpjm1     
112            DO ji = 1, fs_jpim1   ! vector opt.
113               zt1(ji,jj,jk) = umask(ji,jj,jk) * ( tb(ji+1,jj,jk) - tb(ji,jj,jk) )
114               zs1(ji,jj,jk) = umask(ji,jj,jk) * ( sb(ji+1,jj,jk) - sb(ji,jj,jk) )
115               zt2(ji,jj,jk) = vmask(ji,jj,jk) * ( tb(ji,jj+1,jk) - tb(ji,jj,jk) )
116               zs2(ji,jj,jk) = vmask(ji,jj,jk) * ( sb(ji,jj+1,jk) - sb(ji,jj,jk) )
117            END DO
118         END DO
119      END DO
120      ! bottom values
121      zt1(:,:,jpk) = 0.e0    ;    zt2(:,:,jpk) = 0.e0
122      zs1(:,:,jpk) = 0.e0    ;    zs2(:,:,jpk) = 0.e0
123
124      ! lateral boundary conditions on zt1, zt2 ; zs1, zs2   (changed sign)
125      CALL lbc_lnk( zt1, 'U', -1. )   ;   CALL lbc_lnk( zs1, 'U', -1. )
126      CALL lbc_lnk( zt2, 'V', -1. )   ;   CALL lbc_lnk( zs2, 'V', -1. )
127
128      ! Slopes
129      ! interior values
130      DO jk = 1, jpkm1
131         DO jj = 2, jpj
132            DO ji = fs_2, jpi   ! vector opt.
133               ztp1(ji,jj,jk) =                    ( zt1(ji,jj,jk) + zt1(ji-1,jj  ,jk) )   &
134                  &           * ( 0.25 + SIGN( 0.25, zt1(ji,jj,jk) * zt1(ji-1,jj  ,jk) ) )
135               zsp1(ji,jj,jk) =                    ( zs1(ji,jj,jk) + zs1(ji-1,jj  ,jk) )   &
136                  &           * ( 0.25 + SIGN( 0.25, zs1(ji,jj,jk) * zs1(ji-1,jj  ,jk) ) )
137               ztp2(ji,jj,jk) =                    ( zt2(ji,jj,jk) + zt2(ji  ,jj-1,jk) )   &
138                  &           * ( 0.25 + SIGN( 0.25, zt2(ji,jj,jk) * zt2(ji  ,jj-1,jk) ) )
139               zsp2(ji,jj,jk) =                    ( zs2(ji,jj,jk) + zs2(ji  ,jj-1,jk) )   &
140                  &           * ( 0.25 + SIGN( 0.25, zs2(ji,jj,jk) * zs2(ji  ,jj-1,jk) ) )
141            END DO
142         END DO
143      END DO
144      ! bottom values
145      ztp1(:,:,jpk) = 0.e0    ;    ztp2(:,:,jpk) = 0.e0
146      zsp1(:,:,jpk) = 0.e0    ;    zsp2(:,:,jpk) = 0.e0
147
148      ! Slopes limitation
149      DO jk = 1, jpkm1
150         DO jj = 2, jpj
151            DO ji = fs_2, jpi   ! vector opt.
152               ztp1(ji,jj,jk) = SIGN( 1., ztp1(ji,jj,jk) )   &
153                  &           * MIN(    ABS( ztp1(ji  ,jj,jk) ),   &
154                  &                  2.*ABS( zt1 (ji-1,jj,jk) ),   &
155                  &                  2.*ABS( zt1 (ji  ,jj,jk) ) )
156               zsp1(ji,jj,jk) = SIGN( 1., zsp1(ji,jj,jk) )   &
157                  &           * MIN(    ABS( zsp1(ji  ,jj,jk) ),   &
158                  &                  2.*ABS( zs1 (ji-1,jj,jk) ),   &
159                  &                  2.*ABS( zs1 (ji  ,jj,jk) ) )
160               ztp2(ji,jj,jk) = SIGN( 1., ztp2(ji,jj,jk) )   &
161                  &           * MIN(    ABS( ztp2(ji,jj  ,jk) ),   &
162                  &                  2.*ABS( zt2 (ji,jj-1,jk) ),   &
163                  &                  2.*ABS( zt2 (ji,jj  ,jk) ) )
164               zsp2(ji,jj,jk) = SIGN( 1., zsp2(ji,jj,jk) )   &
165                  &           * MIN(    ABS( zsp2(ji,jj  ,jk) ),   &
166                  &                  2.*ABS( zs2 (ji,jj-1,jk) ),   &
167                  &                  2.*ABS( zs2 (ji,jj  ,jk) ) )
168            END DO
169         END DO
170      END DO       
171
172      ! Advection terms
173      ! interior values
174      DO jk = 1, jpkm1
175         zstep  = z2 * rdttra(jk)
176         DO jj = 2, jpjm1     
177            DO ji = fs_2, fs_jpim1   ! vector opt.
178               ! volume fluxes
179#if defined key_zco
180               zeu = e2u(ji,jj)                   * pun(ji,jj,jk)
181               zev = e1v(ji,jj)                   * pvn(ji,jj,jk)
182#else
183               zeu = e2u(ji,jj) * fse3u(ji,jj,jk) * pun(ji,jj,jk)
184               zev = e1v(ji,jj) * fse3v(ji,jj,jk) * pvn(ji,jj,jk)
185#endif
186               ! MUSCL fluxes
187               z0u = SIGN( 0.5, pun(ji,jj,jk) )           
188               zalpha = 0.5 - z0u
189               zu  = z0u - 0.5 * pun(ji,jj,jk) * zstep / e1u(ji,jj)
190               zzt1 = tb(ji+1,jj,jk) + zu*ztp1(ji+1,jj,jk)
191               zzt2 = tb(ji  ,jj,jk) + zu*ztp1(ji  ,jj,jk)
192               zzs1 = sb(ji+1,jj,jk) + zu*zsp1(ji+1,jj,jk)
193               zzs2 = sb(ji  ,jj,jk) + zu*zsp1(ji  ,jj,jk)
194               zt1(ji,jj,jk) = zeu * ( zalpha * zzt1 + (1.-zalpha) * zzt2 )
195               zs1(ji,jj,jk) = zeu * ( zalpha * zzs1 + (1.-zalpha) * zzs2 )
196
197               z0v = SIGN( 0.5, pvn(ji,jj,jk) )           
198               zalpha = 0.5 - z0v
199               zv  = z0v - 0.5 * pvn(ji,jj,jk) * zstep / e2v(ji,jj)
200               zzt1 = tb(ji,jj+1,jk) + zv*ztp2(ji,jj+1,jk)
201               zzt2 = tb(ji,jj  ,jk) + zv*ztp2(ji,jj  ,jk)
202               zzs1 = sb(ji,jj+1,jk) + zv*zsp2(ji,jj+1,jk)
203               zzs2 = sb(ji,jj  ,jk) + zv*zsp2(ji,jj  ,jk)
204               zt2(ji,jj,jk) = zev * ( zalpha * zzt1 + (1.-zalpha) * zzt2 )
205               zs2(ji,jj,jk) = zev * ( zalpha * zzs1 + (1.-zalpha) * zzs2 )
206            END DO
207         END DO
208      END DO
209
210      !!!!  centered scheme at lateral b.C. if off-shore velocity
211      DO jk = 1, jpkm1
212        DO jj = 2, jpjm1
213            DO ji = fs_2, fs_jpim1   ! vector opt.
214#if defined key_zco
215               IF( umask(ji,jj,jk) == 0. ) THEN
216                  IF( pun(ji+1,jj,jk) > 0. .AND. ji /= jpi ) THEN
217                     zt1(ji+1,jj,jk) = e2u(ji+1,jj) * pun(ji+1,jj,jk) * ( tb(ji+1,jj,jk) + tb(ji+2,jj,jk) ) * 0.5
218                     zs1(ji+1,jj,jk) = e2u(ji+1,jj) * pun(ji+1,jj,jk) * ( sb(ji+1,jj,jk) + sb(ji+2,jj,jk) ) * 0.5
219                  ENDIF
220                  IF( pun(ji-1,jj,jk) < 0. ) THEN
221                     zt1(ji-1,jj,jk) = e2u(ji-1,jj) * pun(ji-1,jj,jk) * ( tb(ji-1,jj,jk) + tb(ji  ,jj,jk) ) * 0.5
222                     zs1(ji-1,jj,jk) = e2u(ji-1,jj) * pun(ji-1,jj,jk) * ( sb(ji-1,jj,jk) + sb(ji  ,jj,jk) ) * 0.5
223                  ENDIF
224               ENDIF
225               IF( vmask(ji,jj,jk) == 0. ) THEN
226                  IF( pvn(ji,jj+1,jk) > 0. .AND. jj /= jpj ) THEN
227                     zt2(ji,jj+1,jk) = e1v(ji,jj+1) * pvn(ji,jj+1,jk) * ( tb(ji,jj+1,jk) + tb(ji,jj+2,jk) ) * 0.5
228                     zs2(ji,jj+1,jk) = e1v(ji,jj+1) * pvn(ji,jj+1,jk) * ( sb(ji,jj+1,jk) + sb(ji,jj+2,jk) ) * 0.5
229                  ENDIF
230                  IF( pvn(ji,jj-1,jk) < 0. ) THEN
231                     zt2(ji,jj-1,jk) = e1v(ji,jj-1) * pvn(ji,jj-1,jk) * ( tb(ji,jj-1,jk) + tb(ji  ,jj,jk) ) * 0.5
232                     zs2(ji,jj-1,jk) = e1v(ji,jj-1) * pvn(ji,jj-1,jk) * ( sb(ji,jj-1,jk) + sb(ji  ,jj,jk) ) * 0.5
233                  ENDIF
234               ENDIF
235#else
236               IF( umask(ji,jj,jk) == 0. ) THEN
237                  IF( pun(ji+1,jj,jk) > 0. .AND. ji /= jpi ) THEN
238                     zt1(ji+1,jj,jk) = e2u(ji+1,jj)* fse3u(ji+1,jj,jk)   &
239                        &            * pun(ji+1,jj,jk) * ( tb(ji+1,jj,jk) + tb(ji+2,jj,jk) ) * 0.5
240                     zs1(ji+1,jj,jk) = e2u(ji+1,jj)* fse3u(ji+1,jj,jk)   &
241                        &            * pun(ji+1,jj,jk) * ( sb(ji+1,jj,jk) + sb(ji+2,jj,jk) ) * 0.5
242                  ENDIF
243                  IF( pun(ji-1,jj,jk) < 0. ) THEN
244                     zt1(ji-1,jj,jk) = e2u(ji-1,jj)* fse3u(ji-1,jj,jk)   &
245                        &            * pun(ji-1,jj,jk) * ( tb(ji-1,jj,jk) + tb(ji  ,jj,jk) ) * 0.5
246                     zs1(ji-1,jj,jk) = e2u(ji-1,jj)* fse3u(ji-1,jj,jk)   &
247                        &            * pun(ji-1,jj,jk) * ( sb(ji-1,jj,jk) + sb(ji  ,jj,jk) ) * 0.5
248                  ENDIF
249               ENDIF
250               IF( vmask(ji,jj,jk) == 0. ) THEN
251                  IF( pvn(ji,jj+1,jk) > 0. .AND. jj /= jpj ) THEN
252                     zt2(ji,jj+1,jk) = e1v(ji,jj+1) * fse3v(ji,jj+1,jk)   &
253                        &            * pvn(ji,jj+1,jk) * ( tb(ji,jj+1,jk) + tb(ji,jj+2,jk) ) * 0.5
254                     zs2(ji,jj+1,jk) = e1v(ji,jj+1) * fse3v(ji,jj+1,jk)   &
255                        &            * pvn(ji,jj+1,jk) * ( sb(ji,jj+1,jk) + sb(ji,jj+2,jk) ) * 0.5
256                  ENDIF
257                  IF( pvn(ji,jj-1,jk) < 0. ) THEN
258                     zt2(ji,jj-1,jk) = e1v(ji,jj-1)* fse3v(ji,jj-1,jk)   &
259                        &            * pvn(ji,jj-1,jk) * ( tb(ji,jj-1,jk) + tb(ji  ,jj,jk) ) * 0.5
260                     zs2(ji,jj-1,jk) = e1v(ji,jj-1)* fse3v(ji,jj-1,jk)   &
261                        &            * pvn(ji,jj-1,jk) * ( sb(ji,jj-1,jk) + sb(ji  ,jj,jk) ) * 0.5
262                  ENDIF
263               ENDIF
264#endif
265            END DO
266         END DO
267      END DO
268
269      ! lateral boundary conditions on zt1, zt2 ; zs1, zs2   (changed sign)
270      CALL lbc_lnk( zt1, 'U', -1. )   ;   CALL lbc_lnk( zs1, 'U', -1. ) 
271      CALL lbc_lnk( zt2, 'V', -1. )   ;   CALL lbc_lnk( zs2, 'V', -1. )
272
273      ! Save MUSCL fluxes to compute i- & j- horizontal
274      ! advection trends in the MLD
275      IF( l_trdtra )   THEN
276         ! save i- terms
277         tladi(:,:,:) = zt1(:,:,:) 
278         sladi(:,:,:) = zs1(:,:,:) 
279         ! save j- terms
280         tladj(:,:,:) = zt2(:,:,:) 
281         sladj(:,:,:) = zs2(:,:,:) 
282      ENDIF
283
284      ! Compute & add the horizontal advective trend
285
286      DO jk = 1, jpkm1
287         DO jj = 2, jpjm1     
288            DO ji = fs_2, fs_jpim1   ! vector opt.
289#if defined key_zco
290               zbtr = 1. / ( e1t(ji,jj)*e2t(ji,jj) )
291#else
292               zbtr = 1. / ( e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk) )
293#endif
294               ! horizontal advective trends
295               zta = - zbtr * ( zt1(ji,jj,jk) - zt1(ji-1,jj  ,jk  )   &
296                  &           + zt2(ji,jj,jk) - zt2(ji  ,jj-1,jk  ) )
297               zsa = - zbtr * ( zs1(ji,jj,jk) - zs1(ji-1,jj  ,jk  )   &
298                  &           + zs2(ji,jj,jk) - zs2(ji  ,jj-1,jk  ) ) 
299               ! add it to the general tracer trends
300               ta(ji,jj,jk) = ta(ji,jj,jk) + zta
301               sa(ji,jj,jk) = sa(ji,jj,jk) + zsa
302            END DO
303        END DO
304      END DO       
305
306      ! Save the horizontal advective trends for diagnostic
307
308      IF( l_trdtra )   THEN
309         ! Recompute the horizontal advection zta & zsa trends computed
310         ! at the step 2. above in making the difference between the new
311         ! trends and the previous one ta()/sa - ztdta()/ztdsa() and add
312         ! the term tn()/sn()*hdivn() to recover the Uh gradh(T/S) trends
313         ztdta(:,:,:) = ta(:,:,:) - ztdta(:,:,:) + tn(:,:,:) * hdivn(:,:,:) 
314         ztdsa(:,:,:) = sa(:,:,:) - ztdsa(:,:,:) + sn(:,:,:) * hdivn(:,:,:)
315
316         CALL trd_mod(ztdta, ztdsa, jpttdlad, 'TRA', kt)
317
318         ! Save the new ta and sa trends
319         ztdta(:,:,:) = ta(:,:,:) 
320         ztdsa(:,:,:) = sa(:,:,:) 
321
322      ENDIF
323
324      IF(ln_ctl) THEN
325         CALL prt_ctl(tab3d_1=ta, clinfo1=' muscl2 had  - Ta: ', mask1=tmask, &
326            &         tab3d_2=sa, clinfo2=' Sa: ', mask2=tmask, clinfo3='tra')
327      ENDIF
328
329      ! "zonal" mean advective heat and salt transport
330      IF( ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN
331         IF( lk_zco ) THEN
332            DO jk = 1, jpkm1
333               DO jj = 2, jpjm1
334                  DO ji = fs_2, fs_jpim1   ! vector opt.
335                    zt2(ji,jj,jk) = zt2(ji,jj,jk) * fse3v(ji,jj,jk)
336                    zs2(ji,jj,jk) = zs2(ji,jj,jk) * fse3v(ji,jj,jk)
337                  END DO
338               END DO
339            END DO
340         ENDIF
341         pht_adv(:) = ptr_vj( zt2(:,:,:) )
342         pst_adv(:) = ptr_vj( zs2(:,:,:) )
343      ENDIF
344
345      ! II. Vertical advective fluxes
346      ! -----------------------------
347     
348      ! First guess of the slope
349      ! interior values
350      DO jk = 2, jpkm1
351         zt1(:,:,jk) = tmask(:,:,jk) * ( tb(:,:,jk-1) - tb(:,:,jk) )
352         zs1(:,:,jk) = tmask(:,:,jk) * ( sb(:,:,jk-1) - sb(:,:,jk) )
353      END DO
354      ! surface & bottom boundary conditions
355      zt1 (:,:, 1 ) = 0.e0    ;    zt1 (:,:,jpk) = 0.e0
356      zs1 (:,:, 1 ) = 0.e0    ;    zs1 (:,:,jpk) = 0.e0
357
358      ! Slopes
359      DO jk = 2, jpkm1
360         DO jj = 1, jpj
361            DO ji = 1, jpi
362               ztp1(ji,jj,jk) =                    ( zt1(ji,jj,jk) + zt1(ji,jj,jk+1) )   &
363                  &           * ( 0.25 + SIGN( 0.25, zt1(ji,jj,jk) * zt1(ji,jj,jk+1) ) )
364               zsp1(ji,jj,jk) =                    ( zs1(ji,jj,jk) + zs1(ji,jj,jk+1) )   &
365                  &           * ( 0.25 + SIGN( 0.25, zs1(ji,jj,jk) * zs1(ji,jj,jk+1) ) )
366            END DO
367         END DO
368      END DO
369
370      ! Slopes limitation
371      ! interior values
372      DO jk = 2, jpkm1
373         DO jj = 1, jpj
374            DO ji = 1, jpi
375               ztp1(ji,jj,jk) = SIGN( 1., ztp1(ji,jj,jk) )   &
376                  &           * MIN(    ABS( ztp1(ji,jj,jk  ) ),   &
377                  &                  2.*ABS( zt1 (ji,jj,jk+1) ),   &
378                  &                  2.*ABS( zt1 (ji,jj,jk  ) ) )
379               zsp1(ji,jj,jk) = SIGN( 1., zsp1(ji,jj,jk) )   &
380                  &           * MIN(    ABS( zsp1(ji,jj,jk  ) ),   &
381                  &                  2.*ABS( zs1 (ji,jj,jk+1) ),   &
382                  &                  2.*ABS( zs1 (ji,jj,jk  ) ) )
383            END DO
384         END DO
385      END DO
386      ! surface values
387      ztp1(:,:,1) = 0.e0
388      zsp1(:,:,1) = 0.e0
389
390      ! vertical advective flux
391      ! interior values
392      DO jk = 1, jpkm1
393         zstep  = z2 * rdttra(jk)
394         DO jj = 2, jpjm1     
395            DO ji = fs_2, fs_jpim1   ! vector opt.
396               zew = pwn(ji,jj,jk+1)
397               z0w = SIGN( 0.5, pwn(ji,jj,jk+1) )
398               zalpha = 0.5 + z0w
399               zw  = z0w - 0.5 * pwn(ji,jj,jk+1)*zstep / fse3w(ji,jj,jk+1)
400               zzt1 = tb(ji,jj,jk+1) + zw*ztp1(ji,jj,jk+1)
401               zzt2 = tb(ji,jj,jk  ) + zw*ztp1(ji,jj,jk  )
402               zzs1 = sb(ji,jj,jk+1) + zw*zsp1(ji,jj,jk+1)
403               zzs2 = sb(ji,jj,jk  ) + zw*zsp1(ji,jj,jk  )
404               zt1(ji,jj,jk+1) = zew * ( zalpha * zzt1 + (1.-zalpha)*zzt2 )
405               zs1(ji,jj,jk+1) = zew * ( zalpha * zzs1 + (1.-zalpha)*zzs2 )
406            END DO
407         END DO
408      END DO
409      DO jk = 2, jpkm1
410        DO jj = 2, jpjm1
411            DO ji = fs_2, fs_jpim1   ! vector opt.
412               IF( tmask(ji,jj,jk+1) == 0. ) THEN
413                  IF( pwn(ji,jj,jk) > 0. ) THEN
414                     zt1(ji,jj,jk) = pwn(ji,jj,jk) * ( tb(ji,jj,jk-1) + tb(ji,jj,jk) ) * 0.5
415                     zs1(ji,jj,jk) = pwn(ji,jj,jk) * ( sb(ji,jj,jk-1) + sb(ji,jj,jk) ) * 0.5
416                  ENDIF
417               ENDIF
418            END DO
419         END DO
420      END DO
421
422      ! surface values
423      IF( lk_dynspg_rl ) THEN                           ! rigid lid : flux set to zero
424         zt1(:,:, 1 ) = 0.e0
425         zs1(:,:, 1 ) = 0.e0
426      ELSE                                              ! free surface
427         zt1(:,:, 1 ) = pwn(:,:,1) * tb(:,:,1)
428         zs1(:,:, 1 ) = pwn(:,:,1) * sb(:,:,1)
429      ENDIF
430
431      ! bottom values
432      zt1(:,:,jpk) = 0.e0
433      zs1(:,:,jpk) = 0.e0
434
435
436      ! Compute & add the vertical advective trend
437
438      DO jk = 1, jpkm1
439         DO jj = 2, jpjm1     
440            DO ji = fs_2, fs_jpim1   ! vector opt.
441               zbtr = 1. / fse3t(ji,jj,jk)
442               ! horizontal advective trends
443               zta = - zbtr * ( zt1(ji,jj,jk) - zt1(ji,jj,jk+1) )
444               zsa = - zbtr * ( zs1(ji,jj,jk) - zs1(ji,jj,jk+1) )
445               ! add it to the general tracer trends
446               ta(ji,jj,jk) =  ta(ji,jj,jk) + zta
447               sa(ji,jj,jk) =  sa(ji,jj,jk) + zsa
448            END DO
449         END DO
450      END DO
451
452      ! Save the vertical advective trends for diagnostic
453
454      IF( l_trdtra )   THEN
455         ! Recompute the vertical advection zta & zsa trends computed
456         ! at the step 2. above in making the difference between the new
457         ! trends and the previous one: ta()/sa - ztdta()/ztdsa() and substract
458         ! the term tn()/sn()*hdivn() to recover the W gradz(T/S) trends
459         ztdta(:,:,:) = ta(:,:,:) - ztdta(:,:,:) - tn(:,:,:) * hdivn(:,:,:) 
460         ztdsa(:,:,:) = sa(:,:,:) - ztdsa(:,:,:) - sn(:,:,:) * hdivn(:,:,:)
461
462         CALL trd_mod(ztdta, ztdsa, jpttdzad, 'TRA', kt)
463      ENDIF
464
465      IF(ln_ctl) THEN
466         CALL prt_ctl(tab3d_1=ta, clinfo1=' muscl2 zad  - Ta: ', mask1=tmask, &
467            &         tab3d_2=sa, clinfo2=' Sa: ', mask2=tmask, clinfo3='tra')
468
469      ENDIF
470
471   END SUBROUTINE tra_adv_muscl2
472
473   !!======================================================================
474END MODULE traadv_muscl2
Note: See TracBrowser for help on using the repository browser.