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

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