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

source: tags/nemo_dev_x5/NEMO/OPA_SRC/TRA/traadv_muscl.F90 @ 1792

Last change on this file since 1792 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: 16.8 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 diaptr          ! 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) ::   zzt1, zzt2, zalpha
82      REAL(wp) ::   zzs1, zzs2
83      REAL(wp) ::   z2
84#if defined key_trdtra || defined key_trdmld
85      REAL(wp) ::   ztai, ztaj, zsai, zsaj
86      REAL(wp) ::   zfui, zfvj
87#endif
88      !!----------------------------------------------------------------------
89
90      IF( kt == nit000 .AND. lwp ) THEN
91         WRITE(numout,*)
92         WRITE(numout,*) 'tra_adv : MUSCL advection scheme'
93         WRITE(numout,*) '~~~~~~~'
94      ENDIF
95
96      IF( neuler == 0 .AND. kt == nit000 ) THEN
97          z2=1.
98      ELSE
99          z2=2.
100      ENDIF
101
102#if defined key_trabbl_adv
103
104      ! Advective bottom boundary layer
105      ! -------------------------------
106      zun(:,:,:) = un (:,:,:) - u_bbl(:,:,:)
107      zvn(:,:,:) = vn (:,:,:) - v_bbl(:,:,:)
108      zwn(:,:,:) = wn (:,:,:) + w_bbl( :,:,:)
109#endif
110
111
112      ! I. Horizontal advective fluxes
113      ! ------------------------------
114
115      ! first guess of the slopes
116      ! interior values
117      DO jk = 1, jpkm1
118         DO jj = 1, jpjm1     
119            DO ji = 1, fs_jpim1   ! vector opt.
120               zt1(ji,jj,jk) = umask(ji,jj,jk) * ( tb(ji+1,jj,jk) - tb(ji,jj,jk) )
121               zs1(ji,jj,jk) = umask(ji,jj,jk) * ( sb(ji+1,jj,jk) - sb(ji,jj,jk) )
122               zt2(ji,jj,jk) = vmask(ji,jj,jk) * ( tb(ji,jj+1,jk) - tb(ji,jj,jk) )
123               zs2(ji,jj,jk) = vmask(ji,jj,jk) * ( sb(ji,jj+1,jk) - sb(ji,jj,jk) )
124            END DO
125         END DO
126      END DO
127      ! bottom values
128      zt1(:,:,jpk) = 0.e0    ;    zt2(:,:,jpk) = 0.e0
129      zs1(:,:,jpk) = 0.e0    ;    zs2(:,:,jpk) = 0.e0
130
131      ! lateral boundary conditions on zt1, zt2 ; zs1, zs2   (changed sign)
132      CALL lbc_lnk( zt1, 'U', -1. )   ;   CALL lbc_lnk( zs1, 'U', -1. )
133      CALL lbc_lnk( zt2, 'V', -1. )   ;   CALL lbc_lnk( zs2, 'V', -1. )
134
135      ! Slopes
136      ! interior values
137      DO jk = 1, jpkm1
138         DO jj = 2, jpj
139            DO ji = fs_2, jpi   ! vector opt.
140               ztp1(ji,jj,jk) =                    ( zt1(ji,jj,jk) + zt1(ji-1,jj  ,jk) )   &
141                  &           * ( 0.25 + SIGN( 0.25, zt1(ji,jj,jk) * zt1(ji-1,jj  ,jk) ) )
142               zsp1(ji,jj,jk) =                    ( zs1(ji,jj,jk) + zs1(ji-1,jj  ,jk) )   &
143                  &           * ( 0.25 + SIGN( 0.25, zs1(ji,jj,jk) * zs1(ji-1,jj  ,jk) ) )
144               ztp2(ji,jj,jk) =                    ( zt2(ji,jj,jk) + zt2(ji  ,jj-1,jk) )   &
145                  &           * ( 0.25 + SIGN( 0.25, zt2(ji,jj,jk) * zt2(ji  ,jj-1,jk) ) )
146               zsp2(ji,jj,jk) =                    ( zs2(ji,jj,jk) + zs2(ji  ,jj-1,jk) )   &
147                  &           * ( 0.25 + SIGN( 0.25, zs2(ji,jj,jk) * zs2(ji  ,jj-1,jk) ) )
148            END DO
149         END DO
150      END DO
151      ! bottom values
152      ztp1(:,:,jpk) = 0.e0    ;    ztp2(:,:,jpk) = 0.e0
153      zsp1(:,:,jpk) = 0.e0    ;    zsp2(:,:,jpk) = 0.e0
154
155      ! Slopes limitation
156      DO jk = 1, jpkm1
157         DO jj = 2, jpj
158            DO ji = fs_2, jpi   ! vector opt.
159               ztp1(ji,jj,jk) = SIGN( 1., ztp1(ji,jj,jk) )   &
160                  &           * MIN(    ABS( ztp1(ji  ,jj,jk) ),   &
161                  &                  2.*ABS( zt1 (ji-1,jj,jk) ),   &
162                  &                  2.*ABS( zt1 (ji  ,jj,jk) ) )
163               zsp1(ji,jj,jk) = SIGN( 1., zsp1(ji,jj,jk) )   &
164                  &           * MIN(    ABS( zsp1(ji  ,jj,jk) ),   &
165                  &                  2.*ABS( zs1 (ji-1,jj,jk) ),   &
166                  &                  2.*ABS( zs1 (ji  ,jj,jk) ) )
167               ztp2(ji,jj,jk) = SIGN( 1., ztp2(ji,jj,jk) )   &
168                  &           * MIN(    ABS( ztp2(ji,jj  ,jk) ),   &
169                  &                  2.*ABS( zt2 (ji,jj-1,jk) ),   &
170                  &                  2.*ABS( zt2 (ji,jj  ,jk) ) )
171               zsp2(ji,jj,jk) = SIGN( 1., zsp2(ji,jj,jk) )   &
172                  &           * MIN(    ABS( zsp2(ji,jj  ,jk) ),   &
173                  &                  2.*ABS( zs2 (ji,jj-1,jk) ),   &
174                  &                  2.*ABS( zs2 (ji,jj  ,jk) ) )
175            END DO
176         END DO
177      END DO       
178
179      ! Advection terms
180      ! interior values
181      DO jk = 1, jpkm1
182         zstep  = z2 * rdttra(jk)
183         DO jj = 2, jpjm1     
184            DO ji = fs_2, fs_jpim1   ! vector opt.
185               ! volume fluxes
186#if defined key_s_coord || defined key_partial_steps
187               zeu = e2u(ji,jj) * fse3u(ji,jj,jk) * zun(ji,jj,jk)
188               zev = e1v(ji,jj) * fse3v(ji,jj,jk) * zvn(ji,jj,jk)
189#else
190               zeu = e2u(ji,jj) * zun(ji,jj,jk)
191               zev = e1v(ji,jj) * zvn(ji,jj,jk)
192#endif
193               ! MUSCL fluxes
194               z0u = SIGN( 0.5, zun(ji,jj,jk) )           
195               zalpha = 0.5 - z0u
196               zu  = z0u - 0.5 * zun(ji,jj,jk) * zstep / e1u(ji,jj)
197               zzt1 = tb(ji+1,jj,jk) + zu*ztp1(ji+1,jj,jk)
198               zzt2 = tb(ji  ,jj,jk) + zu*ztp1(ji  ,jj,jk)
199               zzs1 = sb(ji+1,jj,jk) + zu*zsp1(ji+1,jj,jk)
200               zzs2 = sb(ji  ,jj,jk) + zu*zsp1(ji  ,jj,jk)
201               zt1(ji,jj,jk) = zeu * ( zalpha * zzt1 + (1.-zalpha) * zzt2 )
202               zs1(ji,jj,jk) = zeu * ( zalpha * zzs1 + (1.-zalpha) * zzs2 )
203
204               z0v = SIGN( 0.5, zvn(ji,jj,jk) )           
205               zalpha = 0.5 - z0v
206               zv  = z0v - 0.5 * zvn(ji,jj,jk) * zstep / e2v(ji,jj)
207               zzt1 = tb(ji,jj+1,jk) + zv*ztp2(ji,jj+1,jk)
208               zzt2 = tb(ji,jj  ,jk) + zv*ztp2(ji,jj  ,jk)
209               zzs1 = sb(ji,jj+1,jk) + zv*zsp2(ji,jj+1,jk)
210               zzs2 = sb(ji,jj  ,jk) + zv*zsp2(ji,jj  ,jk)
211               zt2(ji,jj,jk) = zev * ( zalpha * zzt1 + (1.-zalpha) * zzt2 )
212               zs2(ji,jj,jk) = zev * ( zalpha * zzs1 + (1.-zalpha) * zzs2 )
213            END DO
214         END DO
215      END DO
216
217      ! lateral boundary conditions on zt1, zt2 ; zs1, zs2   (changed sign)
218      CALL lbc_lnk( zt1, 'U', -1. )   ;   CALL lbc_lnk( zs1, 'U', -1. ) 
219      CALL lbc_lnk( zt2, 'V', -1. )   ;   CALL lbc_lnk( zs2, 'V', -1. )
220
221      ! Compute & add the horizontal advective trend
222
223      DO jk = 1, jpkm1
224         DO jj = 2, jpjm1     
225            DO ji = fs_2, fs_jpim1   ! vector opt.
226#if defined key_s_coord || defined key_partial_steps
227               zbtr = 1. / ( e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk) )
228#else
229               zbtr = 1. / ( e1t(ji,jj)*e2t(ji,jj) )
230#endif
231               ! horizontal advective trends
232               zta = - zbtr * ( zt1(ji,jj,jk) - zt1(ji-1,jj  ,jk  )   &
233                  &           + zt2(ji,jj,jk) - zt2(ji  ,jj-1,jk  ) )
234               zsa = - zbtr * ( zs1(ji,jj,jk) - zs1(ji-1,jj  ,jk  )   &
235                  &           + zs2(ji,jj,jk) - zs2(ji  ,jj-1,jk  ) ) 
236               ! add it to the general tracer trends
237               ta(ji,jj,jk) = ta(ji,jj,jk) + zta
238               sa(ji,jj,jk) = sa(ji,jj,jk) + zsa
239#if defined key_trdtra || defined key_trdmld
240               ! save the horizontal advective trend of tracer
241               ttrd(ji,jj,jk,1) = zta + tn(ji,jj,jk) * hdivn(ji,jj,jk)
242               strd(ji,jj,jk,1) = zsa + sn(ji,jj,jk) * hdivn(ji,jj,jk)
243               ! recompute the trends in i- and j-direction as Uh gradh(T)
244#   if defined key_s_coord || defined key_partial_steps
245               zfui =  e2u(ji  ,jj) * fse3u(ji,  jj,jk) * un(ji,  jj,jk)   &
246                  & -  e2u(ji-1,jj) * fse3u(ji-1,jj,jk) * un(ji-1,jj,jk)
247               zfvj =  e1v(ji,jj  ) * fse3v(ji,jj  ,jk) * vn(ji,jj  ,jk)   &
248                  & -  e1v(ji,jj-1) * fse3v(ji,jj-1,jk) * vn(ji,jj-1,jk)
249#   else
250               zfui = e2u(ji  ,jj) * un(ji,  jj,jk)   &
251                  & - e2u(ji-1,jj) * un(ji-1,jj,jk)
252               zfvj = e1v(ji,jj  ) * vn(ji,jj  ,jk)   &
253                  & - e1v(ji,jj-1) * vn(ji,jj-1,jk)
254#   endif
255               ztai =-zbtr * (  zt1(ji,jj,jk) - zt1(ji-1,jj  ,jk) - tn(ji,jj,jk) * zfui  )
256               ztaj =-zbtr * (  zt2(ji,jj,jk) - zt2(ji  ,jj-1,jk) - tn(ji,jj,jk) * zfvj  )
257               zsai =-zbtr * (  zs1(ji,jj,jk) - zs1(ji-1,jj  ,jk) - sn(ji,jj,jk) * zfui  )
258               zsaj =-zbtr * (  zs2(ji,jj,jk) - zs2(ji  ,jj-1,jk) - sn(ji,jj,jk) * zfvj  )
259               ! save i- and j- advective trends computed as Uh gradh(T)
260               ttrdh(ji,jj,jk,1) = ztai
261               ttrdh(ji,jj,jk,2) = ztaj
262               strdh(ji,jj,jk,1) = zsai
263               strdh(ji,jj,jk,2) = zsaj
264#endif
265            END DO
266        END DO
267      END DO       
268
269      IF(l_ctl) THEN
270         zta = SUM( ta(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
271         zsa = SUM( sa(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
272         WRITE(numout,*) ' had  - Ta: ', zta-t_ctl, ' Sa: ', zsa-s_ctl, ' muscl'
273         t_ctl = zta   ;   s_ctl = zsa
274      ENDIF
275
276      ! "zonal" mean advective heat and salt transport
277      IF( ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN
278# if defined key_s_coord || defined key_partial_steps
279         pht_adv(:) = ptr_vj( zt2(:,:,:) )
280         pst_adv(:) = ptr_vj( zs2(:,:,:) )
281# else
282         DO jk = 1, jpkm1
283            DO jj = 2, jpjm1
284               DO ji = fs_2, fs_jpim1   ! vector opt.
285                 zt2(ji,jj,jk) = zt2(ji,jj,jk) * fse3v(ji,jj,jk)
286                 zs2(ji,jj,jk) = zs2(ji,jj,jk) * fse3v(ji,jj,jk)
287               END DO
288            END DO
289         END DO
290         pht_adv(:) = ptr_vj( zt2(:,:,:) )
291         pst_adv(:) = ptr_vj( zs2(:,:,:) )
292# endif
293      ENDIF
294
295      ! II. Vertical advective fluxes
296      ! -----------------------------
297     
298      ! First guess of the slope
299      ! interior values
300      DO jk = 2, jpkm1
301         zt1(:,:,jk) = tmask(:,:,jk) * ( tb(:,:,jk-1) - tb(:,:,jk) )
302         zs1(:,:,jk) = tmask(:,:,jk) * ( sb(:,:,jk-1) - sb(:,:,jk) )
303      END DO
304      ! surface & bottom boundary conditions
305      zt1 (:,:, 1 ) = 0.e0    ;    zt1 (:,:,jpk) = 0.e0
306      zs1 (:,:, 1 ) = 0.e0    ;    zs1 (:,:,jpk) = 0.e0
307
308      ! Slopes
309      DO jk = 2, jpkm1
310         DO jj = 1, jpj
311            DO ji = 1, jpi
312               ztp1(ji,jj,jk) =                    ( zt1(ji,jj,jk) + zt1(ji,jj,jk+1) )   &
313                  &           * ( 0.25 + SIGN( 0.25, zt1(ji,jj,jk) * zt1(ji,jj,jk+1) ) )
314               zsp1(ji,jj,jk) =                    ( zs1(ji,jj,jk) + zs1(ji,jj,jk+1) )   &
315                  &           * ( 0.25 + SIGN( 0.25, zs1(ji,jj,jk) * zs1(ji,jj,jk+1) ) )
316            END DO
317         END DO
318      END DO
319
320      ! Slopes limitation
321      ! interior values
322      DO jk = 2, jpkm1
323         DO jj = 1, jpj
324            DO ji = 1, jpi
325               ztp1(ji,jj,jk) = SIGN( 1., ztp1(ji,jj,jk) )   &
326                  &           * MIN(    ABS( ztp1(ji,jj,jk  ) ),   &
327                  &                  2.*ABS( zt1 (ji,jj,jk+1) ),   &
328                  &                  2.*ABS( zt1 (ji,jj,jk  ) ) )
329               zsp1(ji,jj,jk) = SIGN( 1., zsp1(ji,jj,jk) )   &
330                  &           * MIN(    ABS( zsp1(ji,jj,jk  ) ),   &
331                  &                  2.*ABS( zs1 (ji,jj,jk+1) ),   &
332                  &                  2.*ABS( zs1 (ji,jj,jk  ) ) )
333            END DO
334         END DO
335      END DO
336      ! surface values
337      ztp1(:,:,1) = 0.e0
338      zsp1(:,:,1) = 0.e0
339
340      ! vertical advective flux
341      ! interior values
342      DO jk = 1, jpkm1
343         zstep  = z2 * rdttra(jk)
344         DO jj = 2, jpjm1     
345            DO ji = fs_2, fs_jpim1   ! vector opt.
346               zew = zwn(ji,jj,jk+1)
347               z0w = SIGN( 0.5, zwn(ji,jj,jk+1) )
348               zalpha = 0.5 + z0w
349               zw  = z0w - 0.5 * zwn(ji,jj,jk+1)*zstep / fse3w(ji,jj,jk+1)
350               zzt1 = tb(ji,jj,jk+1) + zw*ztp1(ji,jj,jk+1)
351               zzt2 = tb(ji,jj,jk  ) + zw*ztp1(ji,jj,jk  )
352               zzs1 = sb(ji,jj,jk+1) + zw*zsp1(ji,jj,jk+1)
353               zzs2 = sb(ji,jj,jk  ) + zw*zsp1(ji,jj,jk  )
354               zt1(ji,jj,jk+1) = zew * ( zalpha * zzt1 + (1.-zalpha)*zzt2 )
355               zs1(ji,jj,jk+1) = zew * ( zalpha * zzs1 + (1.-zalpha)*zzs2 )
356            END DO
357         END DO
358      END DO
359      ! surface values
360      IF( lk_dynspg_fsc .OR. lk_dynspg_fsc_tsk ) THEN   ! free surface-constant volume
361         zt1(:,:, 1 ) = zwn(:,:,1) * tb(:,:,1)
362         zs1(:,:, 1 ) = zwn(:,:,1) * sb(:,:,1)
363      ELSE                                  ! rigid lid : flux set to zero
364         zt1(:,:, 1 ) = 0.e0
365         zs1(:,:, 1 ) = 0.e0
366      ENDIF
367      ! bottom values
368      zt1(:,:,jpk) = 0.e0
369      zs1(:,:,jpk) = 0.e0
370
371
372      ! Compute & add the vertical advective trend
373
374      DO jk = 1, jpkm1
375         DO jj = 2, jpjm1     
376            DO ji = fs_2, fs_jpim1   ! vector opt.
377               zbtr = 1. / fse3t(ji,jj,jk)
378               ! horizontal advective trends
379               zta = - zbtr * ( zt1(ji,jj,jk) - zt1(ji,jj,jk+1) )
380               zsa = - zbtr * ( zs1(ji,jj,jk) - zs1(ji,jj,jk+1) )
381               ! add it to the general tracer trends
382               ta(ji,jj,jk) =  ta(ji,jj,jk) + zta
383               sa(ji,jj,jk) =  sa(ji,jj,jk) + zsa
384#if defined key_trdtra || defined key_trdmld
385               ! save the vertical advective trends computed as w gradz(T)
386               ttrd(ji,jj,jk,2) = zta - tn(ji,jj,jk) * hdivn(ji,jj,jk)
387               strd(ji,jj,jk,2) = zsa - sn(ji,jj,jk) * hdivn(ji,jj,jk)
388#endif
389            END DO
390         END DO
391      END DO
392
393      IF(l_ctl) THEN
394         zta = SUM( ta(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
395         zsa = SUM( sa(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
396         WRITE(numout,*) ' zad  - Ta: ', zta-t_ctl, ' Sa: ', zsa-s_ctl, ' muscl'
397         t_ctl = zta   ;   s_ctl = zsa
398      ENDIF
399
400   END SUBROUTINE tra_adv_muscl
401
402   !!======================================================================
403END MODULE traadv_muscl
Note: See TracBrowser for help on using the repository browser.