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_muscl.F90 in tags/nemo_dev_x2/NEMO/OPA_SRC/TRA – NEMO

source: tags/nemo_dev_x2/NEMO/OPA_SRC/TRA/traadv_muscl.F90 @ 718

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

CT : BUGFIX041 : Correction on salinity and temperature trends for cpp keys key_trdtra or key_trdmld
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 16.9 KB
Line 
1MODULE traadv_muscl
2   !!==============================================================================
3   !!                       ***  MODULE  traadv_muscl  ***
4   !! Ocean active tracers:  horizontal & vertical advective trend
5   !!==============================================================================
6
7   !!----------------------------------------------------------------------
8   !!   tra_adv_muscl : update the tracer trend with the horizontal
9   !!                   and vertical advection trends using MUSCL 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 boundary layer
19   USE lib_mpp
20   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
21   USE ptr             ! poleward transport diagnostics
22
23   IMPLICIT NONE
24   PRIVATE
25
26   !! * Accessibility
27   PUBLIC tra_adv_muscl  ! routine called by step.F90
28
29   !! * Substitutions
30#  include "domzgr_substitute.h90"
31#  include "vectopt_loop_substitute.h90"
32   !!----------------------------------------------------------------------
33   !!   OPA 9.0 , LODYC-IPSL (2003)
34   !!----------------------------------------------------------------------
35
36CONTAINS
37
38   SUBROUTINE tra_adv_muscl( kt )
39      !!----------------------------------------------------------------------
40      !!                    ***  ROUTINE tra_adv_muscl  ***
41      !!         
42      !! ** Purpose :   Compute the now trend due to total advection of T and
43      !!      S using a MUSCL scheme (Monotone Upstream-centered Scheme for
44      !!      Conservation Laws) and add it to the general tracer trend.
45      !!
46      !! ** Method  :
47      !!
48      !! ** Action  : - update (ta,sa) with the now advective tracer trends
49      !!              - save trends in (ttrdh,ttrd,strdh,strd) ('key_trdtra')
50      !!
51      !! References :               
52      !!      Estubier, A., and M. Levy, Notes Techn. Pole de Modelisation
53      !! IPSL, Sept. 2000 (http://www.lodyc.jussieu.fr/opa)
54      !!
55      !! History :
56      !!        !  06-00  (A.Estublier)  for passive tracers
57      !!        !  01-08  (E.Durand G.Madec)  adapted for T & S
58      !!   8.5  !  02-06  (G. Madec)  F90: Free form and module
59      !!----------------------------------------------------------------------
60      !! * modules used
61#if defined key_trabbl_adv
62      USE oce                , zun => ua,  &  ! use ua as workspace
63         &                     zvn => va      ! use va as workspace
64      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwn
65#else
66      USE oce                , zun => un,  &  ! When no bbl, zun == un
67                               zvn => vn,  &  !              zvn == vn
68                               zwn => wn      !              zwn == wn
69#endif
70
71      !! * Arguments
72      INTEGER, INTENT( in ) ::   kt         ! ocean time-step
73
74      !! * Local declarations
75      INTEGER ::   ji, jj, jk               ! dummy loop indices
76      REAL(wp), DIMENSION (jpi,jpj,jpk) ::   &
77         zt1, zt2, ztp1, ztp2,   &
78         zs1, zs2, zsp1, zsp2
79      REAL(wp) ::   zu, zv, zw, zeu, zev, zew, zbtr, zstep, zta, zsa
80      REAL(wp) ::   z0u, z0v, z0w
81      REAL(wp) ::   z1u, z1v, z1w
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 : MUSCL 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, fs_jpim1   ! vector opt.
141               z0u = zt1(ji,jj,jk) * zt1(ji-1,jj,jk)
142               IF( z0u > 0. ) THEN
143                  ztp1(ji,jj,jk) = 0.5 * ( zt1(ji,jj,jk)+zt1(ji-1,jj,jk) )
144               ELSE
145                  ztp1(ji,jj,jk) = 0.e0
146               ENDIF
147               z1u = zs1(ji,jj,jk) * zs1(ji-1,jj,jk)
148               IF( z1u > 0. ) THEN
149                  zsp1(ji,jj,jk) = 0.5 * ( zs1(ji,jj,jk)+zs1(ji-1,jj,jk) )
150               ELSE
151                  zsp1(ji,jj,jk) = 0.e0
152               ENDIF
153
154               z0v = zt2(ji,jj,jk) * zt2(ji,jj-1,jk)
155               IF( z0v > 0. ) THEN
156                  ztp2(ji,jj,jk) = 0.5 * ( zt2(ji,jj,jk)+zt2(ji,jj-1,jk) )
157               ELSE
158                  ztp2(ji,jj,jk) = 0.e0
159               ENDIF
160               z1v = zs2(ji,jj,jk) * zs2(ji,jj-1,jk)
161               IF( z1v > 0. ) THEN
162                  zsp2(ji,jj,jk) = 0.5 * ( zs2(ji,jj,jk)+zs2(ji,jj-1,jk) )
163               ELSE
164                  zsp2(ji,jj,jk) = 0.e0
165               ENDIF
166            END DO
167         END DO
168      END DO
169      ! bottom values
170      ztp1(:,:,jpk) = 0.e0    ;    ztp2(:,:,jpk) = 0.e0
171      zsp1(:,:,jpk) = 0.e0    ;    zsp2(:,:,jpk) = 0.e0
172
173! Slopes limitation
174      DO jk = 1, jpkm1
175         DO jj = 2, jpj
176            DO ji = fs_2, fs_jpim1   ! vector opt.
177               ztp1(ji,jj,jk) = SIGN( 1., ztp1(ji,jj,jk) )   &
178                  &           * MIN(    ABS( ztp1(ji  ,jj,jk) ),   &
179                  &                  2.*ABS( zt1 (ji-1,jj,jk) ),   &
180                  &                  2.*ABS( zt1 (ji  ,jj,jk) ) )
181               zsp1(ji,jj,jk) = SIGN( 1., zsp1(ji,jj,jk) )   &
182                  &           * MIN(    ABS( zsp1(ji  ,jj,jk) ),   &
183                  &                  2.*ABS( zs1 (ji-1,jj,jk) ),   &
184                  &                  2.*ABS( zs1 (ji  ,jj,jk) ) )
185               ztp2(ji,jj,jk) = SIGN( 1., ztp2(ji,jj,jk) )   &
186                  &           * MIN(    ABS( ztp2(ji,jj  ,jk) ),   &
187                  &                  2.*ABS( zt2 (ji,jj-1,jk) ),   &
188                  &                  2.*ABS( zt2 (ji,jj  ,jk) ) )
189               zsp2(ji,jj,jk) = SIGN( 1., zsp2(ji,jj,jk) )   &
190                  &           * MIN(    ABS( zsp2(ji,jj  ,jk) ),   &
191                  &                  2.*ABS( zs2 (ji,jj-1,jk) ),   &
192                  &                  2.*ABS( zs2 (ji,jj  ,jk) ) )
193            END DO
194         END DO
195      END DO       
196
197      ! Advection terms
198      ! interior values
199      DO jk = 1, jpkm1
200         zstep  = z2 * rdttra(jk)
201         DO jj = 2, jpjm1     
202            DO ji = fs_2, fs_jpim1   ! vector opt.
203               ! volume fluxes
204#if defined key_s_coord || defined key_partial_steps
205               zeu = e2u(ji,jj) * fse3u(ji,jj,jk) * zun(ji,jj,jk)
206               zev = e1v(ji,jj) * fse3v(ji,jj,jk) * zvn(ji,jj,jk)
207#else
208               zeu = e2u(ji,jj) * zun(ji,jj,jk)
209               zev = e1v(ji,jj) * zvn(ji,jj,jk)
210#endif
211               ! MUSCL fluxes
212               z0u = SIGN( 0.5, zun(ji,jj,jk) )           
213               zalpha = 0.5 - z0u
214               zu  = z0u - 0.5 * zun(ji,jj,jk) * zstep / e1u(ji,jj)
215               zzt1 = tb(ji+1,jj,jk) + zu*ztp1(ji+1,jj,jk)
216               zzt2 = tb(ji  ,jj,jk) + zu*ztp1(ji  ,jj,jk)
217               zzs1 = sb(ji+1,jj,jk) + zu*zsp1(ji+1,jj,jk)
218               zzs2 = sb(ji  ,jj,jk) + zu*zsp1(ji  ,jj,jk)
219               zt1(ji,jj,jk) = zeu * ( zalpha * zzt1 + (1.-zalpha) * zzt2 )
220               zs1(ji,jj,jk) = zeu * ( zalpha * zzs1 + (1.-zalpha) * zzs2 )
221
222               z0v = SIGN( 0.5, zvn(ji,jj,jk) )           
223               zalpha = 0.5 - z0v
224               zv  = z0v - 0.5 * zvn(ji,jj,jk) * zstep / e2v(ji,jj)
225               zzt1 = tb(ji,jj+1,jk) + zv*ztp2(ji,jj+1,jk)
226               zzt2 = tb(ji,jj  ,jk) + zv*ztp2(ji,jj  ,jk)
227               zzs1 = sb(ji,jj+1,jk) + zv*zsp2(ji,jj+1,jk)
228               zzs2 = sb(ji,jj  ,jk) + zv*zsp2(ji,jj  ,jk)
229               zt2(ji,jj,jk) = zev * ( zalpha * zzt1 + (1.-zalpha) * zzt2 )
230               zs2(ji,jj,jk) = zev * ( zalpha * zzs1 + (1.-zalpha) * zzs2 )
231            END DO
232         END DO
233      END DO
234
235      ! lateral boundary conditions on zt1, zt2 ; zs1, zs2   (changed sign)
236      CALL lbc_lnk( zt1, 'U', -1. )   ;   CALL lbc_lnk( zs1, 'U', -1. ) 
237      CALL lbc_lnk( zt2, 'V', -1. )   ;   CALL lbc_lnk( zs2, 'V', -1. )
238
239      ! Compute & add the horizontal advective trend
240
241      DO jk = 1, jpkm1
242         DO jj = 2, jpjm1     
243            DO ji = fs_2, fs_jpim1   ! vector opt.
244#if defined key_s_coord || defined key_partial_steps
245               zbtr = 1. / ( e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk) )
246#else
247               zbtr = 1. / ( e1t(ji,jj)*e2t(ji,jj) )
248#endif
249               ! horizontal advective trends
250               zta = - zbtr * ( zt1(ji,jj,jk) - zt1(ji-1,jj  ,jk  )   &
251                  &           + zt2(ji,jj,jk) - zt2(ji  ,jj-1,jk  ) )
252               zsa = - zbtr * ( zs1(ji,jj,jk) - zs1(ji-1,jj  ,jk  )   &
253                  &           + zs2(ji,jj,jk) - zs2(ji  ,jj-1,jk  ) ) 
254               ! add it to the general tracer trends
255               ta(ji,jj,jk) = ta(ji,jj,jk) + zta
256               sa(ji,jj,jk) = sa(ji,jj,jk) + zsa
257#if defined key_trdtra || defined key_trdmld
258               ! save the horizontal advective trend of tracer
259               ttrd(ji,jj,jk,1) = zta + tn(ji,jj,jk) * hdivn(ji,jj,jk)
260               strd(ji,jj,jk,1) = zsa + sn(ji,jj,jk) * hdivn(ji,jj,jk)
261               ! recompute the trends in i- and j-direction as Uh gradh(T)
262#   if defined key_s_coord || defined key_partial_steps
263               zfui =  e2u(ji  ,jj) * fse3u(ji,  jj,jk) * un(ji,  jj,jk)   &
264                  & -  e2u(ji-1,jj) * fse3u(ji-1,jj,jk) * un(ji-1,jj,jk)
265               zfvj =  e1v(ji,jj  ) * fse3v(ji,jj  ,jk) * vn(ji,jj  ,jk)   &
266                  & -  e1v(ji,jj-1) * fse3v(ji,jj-1,jk) * vn(ji,jj-1,jk)
267#   else
268               zfui = e2u(ji  ,jj) * un(ji,  jj,jk)   &
269                  & - e2u(ji-1,jj) * un(ji-1,jj,jk)
270               zfvj = e1v(ji,jj  ) * vn(ji,jj  ,jk)   &
271                  & - e1v(ji,jj-1) * vn(ji,jj-1,jk)
272#   endif
273               ztai =-zbtr * (  zt1(ji,jj,jk) - zt1(ji-1,jj  ,jk) - tn(ji,jj,jk) * zfui  )
274               ztaj =-zbtr * (  zt2(ji,jj,jk) - zt2(ji  ,jj-1,jk) - tn(ji,jj,jk) * zfvj  )
275               zsai =-zbtr * (  zs1(ji,jj,jk) - zs1(ji-1,jj  ,jk) - sn(ji,jj,jk) * zfui  )
276               zsaj =-zbtr * (  zs2(ji,jj,jk) - zs2(ji  ,jj-1,jk) - sn(ji,jj,jk) * zfvj  )
277               ! save i- and j- advective trends computed as Uh gradh(T)
278               ttrdh(ji,jj,jk,1) = ztai
279               ttrdh(ji,jj,jk,2) = ztaj
280               strdh(ji,jj,jk,1) = zsai
281               strdh(ji,jj,jk,2) = zsaj
282#endif
283            END DO
284        END DO
285      END DO       
286
287#if defined key_diaptr 
288      ! "zonal" mean advective heat and salt transport
289      IF( MOD( kt, nf_ptr ) == 0 ) THEN
290# if defined key_s_coord || defined key_partial_steps
291         pht_adv(:) = prt_vj( zt2(:,:,:) )
292         pst_adv(:) = prt_vj( zs2(:,:,:) )
293# else
294         DO jk = 1, jpkm1
295            DO jj = 2, jpjm1
296               DO ji = fs_2, fs_jpim1   ! vector opt.
297                 zt2(ji,jj,jk) = zt2(ji,jj,jk) * fse3v(ji,jj,jk)
298                 zs2(ji,jj,jk) = zs2(ji,jj,jk) * fse3v(ji,jj,jk)
299               END DO
300            END DO
301         END DO
302         pht_adv(:) = prt_vj( zt2(:,:,:) )
303         pst_adv(:) = prt_vj( zs2(:,:,:) )
304# endif
305      ENDIF
306#endif
307
308
309      ! II. Vertical advective fluxes
310      ! -----------------------------
311     
312      ! First guess of the slope
313      ! interior values
314      DO jk = 2, jpkm1
315         zt1(:,:,jk) = tmask(:,:,jk) * ( tb(:,:,jk-1) - tb(:,:,jk) )
316         zs1(:,:,jk) = tmask(:,:,jk) * ( sb(:,:,jk-1) - sb(:,:,jk) )
317      END DO
318      ! surface & bottom boundary conditions
319      zt1 (:,:, 1 ) = 0.e0    ;    zt1 (:,:,jpk) = 0.e0
320      zs1 (:,:, 1 ) = 0.e0    ;    zs1 (:,:,jpk) = 0.e0
321
322      ! Slopes
323      DO jk = 2, jpkm1
324         DO jj = 1, jpj
325            DO ji = 1, jpi
326               z0w = zt1(ji,jj,jk) * zt1(ji,jj,jk+1) 
327               IF( z0w > 0. ) THEN
328                  ztp1(ji,jj,jk) = 0.5 * ( zt1(ji,jj,jk) + zt1(ji,jj,jk+1) )
329               ELSE
330                  ztp1(ji,jj,jk) = 0.e0
331               ENDIF
332               z1w = zs1(ji,jj,jk) * zs1(ji,jj,jk+1) 
333               IF( z1w > 0. ) THEN
334                  zsp1(ji,jj,jk) = 0.5 * ( zs1(ji,jj,jk) + zs1(ji,jj,jk+1) )
335               ELSE
336                  zsp1(ji,jj,jk) = 0.e0
337               ENDIF
338            END DO
339         END DO
340      END DO
341
342      ! Slopes limitation
343      ! interior values
344      DO jk = 2, jpkm1
345         DO jj = 1, jpj
346            DO ji = 1, jpi
347               ztp1(ji,jj,jk) = SIGN( 1., ztp1(ji,jj,jk) )   &
348                  &           * MIN(    ABS( ztp1(ji,jj,jk  ) ),   &
349                  &                  2.*ABS( zt1 (ji,jj,jk+1) ),   &
350                  &                  2.*ABS( zt1 (ji,jj,jk  ) ) )
351               zsp1(ji,jj,jk) = SIGN( 1., zsp1(ji,jj,jk) )   &
352                  &           * MIN(    ABS( zsp1(ji,jj,jk  ) ),   &
353                  &                  2.*ABS( zs1 (ji,jj,jk+1) ),   &
354                  &                  2.*ABS( zs1 (ji,jj,jk  ) ) )
355            END DO
356         END DO
357      END DO
358      ! surface values
359      ztp1(:,:,1) = 0. 
360      zsp1(:,:,1) = 0.
361
362      ! vertical advective flux
363      ! interior values
364      DO jk = 1, jpkm1
365         zstep  = z2 * rdttra(jk)
366         DO jj = 2, jpjm1     
367            DO ji = fs_2, fs_jpim1   ! vector opt.
368               zew = zwn(ji,jj,jk+1)
369               z0w = SIGN( 0.5, zwn(ji,jj,jk+1) )
370               zalpha = 0.5 + z0w
371               zw  = z0w - 0.5 * zwn(ji,jj,jk+1)*zstep / fse3w(ji,jj,jk+1)
372               zzt1 = tb(ji,jj,jk+1) + zw*ztp1(ji,jj,jk+1)
373               zzt2 = tb(ji,jj,jk  ) + zw*ztp1(ji,jj,jk  )
374               zzs1 = sb(ji,jj,jk+1) + zw*zsp1(ji,jj,jk+1)
375               zzs2 = sb(ji,jj,jk  ) + zw*zsp1(ji,jj,jk  )
376               zt1(ji,jj,jk+1) = zew * ( zalpha * zzt1 + (1.-zalpha)*zzt2 )
377               zs1(ji,jj,jk+1) = zew * ( zalpha * zzs1 + (1.-zalpha)*zzs2 )
378            END DO
379         END DO
380      END DO
381      ! surface values
382      IF( lk_dynspg_fsc .OR. lk_dynspg_fsc_tsk ) THEN   ! free surface-constant volume
383         zt1(:,:, 1 ) = zwn(:,:,1) * tb(:,:,1)
384         zs1(:,:, 1 ) = zwn(:,:,1) * sb(:,:,1)
385      ELSE                                  ! rigid lid : flux set to zero
386         zt1(:,:, 1 ) = 0.e0
387         zs1(:,:, 1 ) = 0.e0
388      ENDIF
389      ! bottom values
390      zt1(:,:,jpk) = 0.e0
391      zs1(:,:,jpk) = 0.e0
392
393
394      ! Compute & add the vertical advective trend
395
396      DO jk = 1, jpkm1
397         DO jj = 2, jpjm1     
398            DO ji = fs_2, fs_jpim1   ! vector opt.
399               zbtr = 1. / fse3t(ji,jj,jk)
400               ! horizontal advective trends
401               zta = - zbtr * ( zt1(ji,jj,jk) - zt1(ji,jj,jk+1) )
402               zsa = - zbtr * ( zs1(ji,jj,jk) - zs1(ji,jj,jk+1) )
403               ! add it to the general tracer trends
404               ta(ji,jj,jk) =  ta(ji,jj,jk) + zta
405               sa(ji,jj,jk) =  sa(ji,jj,jk) + zsa
406#if defined key_trdtra || defined key_trdmld
407               ! save the vertical advective trends computed as w gradz(T)
408               ttrd(ji,jj,jk,2) = zta - tn(ji,jj,jk) * hdivn(ji,jj,jk)
409               strd(ji,jj,jk,2) = zsa - sn(ji,jj,jk) * hdivn(ji,jj,jk)
410#endif
411            END DO
412         END DO
413      END DO
414
415      IF( l_ctl .AND. lwp ) THEN
416         zta = SUM( ta(2:jpim1,2:jpjm1,1:jpkm1) * tmask(2:jpim1,2:jpjm1,1:jpkm1) )
417         zsa = SUM( sa(2:jpim1,2:jpjm1,1:jpkm1) * tmask(2:jpim1,2:jpjm1,1:jpkm1) )
418         WRITE(numout,*) ' zad  - Ta: ', zta-t_ctl, ' Sa: ', zsa-s_ctl, ' muscl'
419         t_ctl = zta   ;   s_ctl = zsa
420      ENDIF
421
422   END SUBROUTINE tra_adv_muscl
423
424   !!======================================================================
425END MODULE traadv_muscl
Note: See TracBrowser for help on using the repository browser.