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

source: trunk/NEMO/OPA_SRC/TRA/traadv_muscl.F90 @ 64

Last change on this file since 64 was 3, checked in by opalod, 20 years ago

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