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.
trcadv_muscl.F90 in trunk/NEMO/TOP_SRC/TRP – NEMO

source: trunk/NEMO/TOP_SRC/TRP/trcadv_muscl.F90 @ 264

Last change on this file since 264 was 261, checked in by opalod, 19 years ago

nemo_v1_update_05:RB+OA:Update and rewritting of (part of) the TOP component
  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 13.3 KB
Line 
1MODULE trcadv_muscl
2   !!==============================================================================
3   !!                       ***  MODULE  trcadv_muscl  ***
4   !! Ocean passive tracers:  horizontal & vertical advective trend
5   !!==============================================================================
6#if defined key_passivetrc
7   !!----------------------------------------------------------------------
8   !!   trc_adv_muscl : update the tracer trend with the horizontal
9   !!                   and vertical advection trends using MUSCL scheme
10   !!----------------------------------------------------------------------
11   !! * Modules used
12   USE oce_trc         ! ocean dynamics and active tracers variables
13   USE trc             ! ocean passive tracers variables
14   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
15   USE trcbbl          ! advective passive tracers in the BBL
16   USE lib_mpp
17
18   IMPLICIT NONE
19   PRIVATE
20
21   !! * Accessibility
22   PUBLIC trc_adv_muscl  ! routine called by trcstp.F90
23
24   !! * Substitutions
25#  include "passivetrc_substitute.h90"
26   !!----------------------------------------------------------------------
27   !!   OPA 9.0 , LODYC-IPSL (2003)
28   !!----------------------------------------------------------------------
29
30CONTAINS
31
32   SUBROUTINE trc_adv_muscl( kt )
33      !!----------------------------------------------------------------------
34      !!                    ***  ROUTINE trc_adv_muscl  ***
35      !!         
36      !! ** Purpose :   Compute the now trend due to total advection of any pas-
37      !!      sive tracer using a MUSCL scheme (Monotone Upstream-centered Scheme
38      !!      for Conservation Laws) and add it to the general tracer trend.
39      !!
40      !! ** Method  :
41      !!
42      !! ** Action  : - update tra with the now advective tracer trends
43      !!              - save trends in trtrd ('key_trc_diatrd')
44      !!
45      !! References :               
46      !!      Estubier, A., and M. Levy, Notes Techn. Pole de Modelisation
47      !! IPSL, Sept. 2000 (http://www.lodyc.jussieu.fr/opa)
48      !!
49      !! History :
50      !!        !  06-00  (A.Estublier)  for passive tracers
51      !!   9.0  !  03-04  (C. Ethe, G. Madec)  F90: Free form and module
52      !!----------------------------------------------------------------------
53      !! * modules used
54#if defined key_trcbbl_adv
55      USE oce_trc            , zun => ua,  &  ! use ua as workspace
56         &                     zvn => va      ! use va as workspace
57      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwn
58#else
59      USE oce_trc            , zun => un,  &  ! When no bbl, zun == un
60                               zvn => vn,  &  !              zvn == vn
61                               zwn => wn      !              zwn == wn
62#endif
63
64      !! * Arguments
65      INTEGER, INTENT( in ) ::   kt         ! ocean time-step
66
67      !! * Local declarations
68      INTEGER ::   ji, jj, jk,jn            ! dummy loop indices
69      REAL(wp), DIMENSION (jpi,jpj,jpk) ::   &
70         zt1, zt2, ztp1, ztp2
71
72      REAL(wp) ::   zu, zv, zw, zeu, zev, zew, zbtr, ztra
73      REAL(wp) ::   z0u, z0v, z0w
74      REAL(wp) ::   zzt1, zzt2, zalpha, z2dtt
75#if defined key_trc_diatrd
76      REAL(wp) ::   ztai, ztaj
77      REAL(wp) ::   zfui, zfvj
78#endif
79
80      !!----------------------------------------------------------------------
81
82
83      IF( kt == nittrc000 .AND. lwp ) THEN
84         WRITE(numout,*)
85         WRITE(numout,*) 'trc_adv : MUSCL advection scheme'
86         WRITE(numout,*) '~~~~~~~'
87      ENDIF
88
89 
90
91#if defined key_trcbbl_adv
92      ! Advective bottom boundary layer
93      ! -------------------------------
94      zun(:,:,:) = un (:,:,:) - u_trc_bbl(:,:,:)
95      zvn(:,:,:) = vn (:,:,:) - v_trc_bbl(:,:,:)
96      zwn(:,:,:) = wn (:,:,:) + w_trc_bbl(:,:,:)
97#endif
98
99 
100
101      DO jn = 1, jptra
102
103         ! I. Horizontal advective fluxes
104         ! ------------------------------
105
106         ! first guess of the slopes
107         ! interior values
108         DO jk = 1, jpkm1
109            DO jj = 1, jpjm1     
110               DO ji = 1, fs_jpim1   ! vector opt.
111                  zt1(ji,jj,jk) = umask(ji,jj,jk) * ( trb(ji+1,jj,jk,jn) - trb(ji,jj,jk,jn) )
112                  zt2(ji,jj,jk) = vmask(ji,jj,jk) * ( trb(ji,jj+1,jk,jn) - trb(ji,jj,jk,jn) )
113               END DO
114            END DO
115         END DO
116         ! bottom values
117         zt1(:,:,jpk) = 0.e0
118         zt2(:,:,jpk) = 0.e0
119
120         ! lateral boundary conditions on zt1, zt2
121         CALL lbc_lnk( zt1, 'U', -1. )   
122         CALL lbc_lnk( zt2, 'V', -1. ) 
123
124
125         ! Slopes
126         ! interior values
127         DO jk = 1, jpkm1
128            DO jj = 2, jpj
129               DO ji = fs_2, jpi   ! vector opt.
130                  ztp1(ji,jj,jk) =                    ( zt1(ji,jj,jk) + zt1(ji-1,jj  ,jk) )   &
131                     &           * ( 0.25 + SIGN( 0.25, zt1(ji,jj,jk) * zt1(ji-1,jj  ,jk) ) )
132                  ztp2(ji,jj,jk) =                    ( zt2(ji,jj,jk) + zt2(ji  ,jj-1,jk) )   &
133                     &           * ( 0.25 + SIGN( 0.25, zt2(ji,jj,jk) * zt2(ji  ,jj-1,jk) ) )
134               END DO
135            END DO
136         END DO
137         ! bottom values
138         ztp1(:,:,jpk) = 0.e0
139         ztp2(:,:,jpk) = 0.e0
140
141         ! Slopes limitation
142         DO jk = 1, jpkm1
143            DO jj = 2, jpj
144               DO ji = fs_2, jpi   ! vector opt.
145                  ztp1(ji,jj,jk) = SIGN( 1., ztp1(ji,jj,jk) )   &
146                     &           * MIN(    ABS( ztp1(ji  ,jj,jk) ),   &
147                     &                  2.*ABS( zt1 (ji-1,jj,jk) ),   &
148                     &                  2.*ABS( zt1 (ji  ,jj,jk) ) )
149
150                  ztp2(ji,jj,jk) = SIGN( 1., ztp2(ji,jj,jk) )   &
151                     &           * MIN(    ABS( ztp2(ji,jj  ,jk) ),   &
152                     &                  2.*ABS( zt2 (ji,jj-1,jk) ),   &
153                     &                  2.*ABS( zt2 (ji,jj  ,jk) ) )
154
155               END DO
156            END DO
157         END DO
158
159         ! Advection terms
160         ! interior values
161         DO jk = 1, jpkm1
162            DO jj = 2, jpjm1     
163               DO ji = fs_2, fs_jpim1   ! vector opt.
164                  ! volume fluxes
165#if defined key_s_coord || defined key_partial_steps
166                  zeu = e2u(ji,jj) * fse3u(ji,jj,jk) * zun(ji,jj,jk)
167                  zev = e1v(ji,jj) * fse3v(ji,jj,jk) * zvn(ji,jj,jk)
168#else
169                  zeu = e2u(ji,jj) * zun(ji,jj,jk)
170                  zev = e1v(ji,jj) * zvn(ji,jj,jk)
171#endif
172                  ! MUSCL fluxes
173                  z2dtt = rdttra(jk) * FLOAT(ndttrc)
174                  z0u = SIGN( 0.5, zun(ji,jj,jk) )           
175                  zalpha = 0.5 - z0u
176                  zu  = z0u - 0.5 * zun(ji,jj,jk) * z2dtt / e1u(ji,jj)
177                  zzt1 = trb(ji+1,jj,jk,jn) + zu*ztp1(ji+1,jj,jk)
178                  zzt2 = trb(ji  ,jj,jk,jn) + zu*ztp1(ji  ,jj,jk)
179                  zt1(ji,jj,jk) = zeu * ( zalpha * zzt1 + (1.-zalpha) * zzt2 )
180                  z0v = SIGN( 0.5, zvn(ji,jj,jk) )           
181                  zalpha = 0.5 - z0v
182                  zv  = z0v - 0.5 * zvn(ji,jj,jk) * z2dtt / e2v(ji,jj)
183                  zzt1 = trb(ji,jj+1,jk,jn) + zv*ztp2(ji,jj+1,jk)
184                  zzt2 = trb(ji,jj  ,jk,jn) + zv*ztp2(ji,jj  ,jk)
185                  zt2(ji,jj,jk) = zev * ( zalpha * zzt1 + (1.-zalpha) * zzt2 )
186               END DO
187            END DO
188         END DO
189
190         ! lateral boundary conditions on zt1, zt2 (changed sign)
191         CALL lbc_lnk( zt1, 'U', -1. ) 
192         CALL lbc_lnk( zt2, 'V', -1. ) 
193
194         ! Compute and add the horizontal advective trend
195
196         DO jk = 1, jpkm1
197            DO jj = 2, jpjm1     
198               DO ji = fs_2, fs_jpim1   ! vector opt.
199#if defined key_s_coord || defined key_partial_steps
200                  zbtr = 1. / ( e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk) )
201#else
202                  zbtr = 1. / ( e1t(ji,jj)*e2t(ji,jj) )
203#endif
204                  ! horizontal advective trends
205                  ztra = - zbtr * ( zt1(ji,jj,jk) - zt1(ji-1,jj  ,jk  )   &
206                     &            + zt2(ji,jj,jk) - zt2(ji  ,jj-1,jk  ) )
207                  ! add it to the general tracer trends
208                  tra(ji,jj,jk,jn) = tra(ji,jj,jk,jn) + ztra
209#if defined key_trc_diatrd
210                  ! recompute the trends in i- and j-direction as Uh gradh(T)
211#   if defined key_s_coord || defined key_partial_steps
212                  zfui =  e2u(ji  ,jj) * fse3u(ji,  jj,jk) * un(ji,  jj,jk)   &
213                     & -  e2u(ji-1,jj) * fse3u(ji-1,jj,jk) * un(ji-1,jj,jk)
214                  zfvj =  e1v(ji,jj  ) * fse3v(ji,jj  ,jk) * vn(ji,jj  ,jk)   &
215                     & -  e1v(ji,jj-1) * fse3v(ji,jj-1,jk) * vn(ji,jj-1,jk)
216#   else
217                  zfui = e2u(ji  ,jj) * un(ji,  jj,jk)   &
218                     & - e2u(ji-1,jj) * un(ji-1,jj,jk)
219                  zfvj = e1v(ji,jj  ) * vn(ji,jj  ,jk)   &
220                     & - e1v(ji,jj-1) * vn(ji,jj-1,jk)
221#   endif
222                  ztai =-zbtr * (  zt1(ji,jj,jk) - zt1(ji-1,jj  ,jk) - trn(ji,jj,jk,jn) * zfui  )
223                  ztaj =-zbtr * (  zt2(ji,jj,jk) - zt2(ji  ,jj-1,jk) - trn(ji,jj,jk,jn) * zfvj  )
224                  ! save i- and j- advective trends computed as Uh gradh(T)
225                  trtrd(ji,jj,jk,jn,1) = ztai
226                  trtrd(ji,jj,jk,jn,2) = ztaj
227#endif
228               END DO
229            END DO
230         END DO
231
232         IF(ln_ctl) THEN         ! print mean trends (used for debugging)
233            ztra = SUM( tra(2:nictle,2:njctle,1:jpkm1,jn) * tmask(2:nictle,2:njctle,1:jpkm1) )
234            WRITE(numout,*) ' trc/had  - ',ctrcnm(jn),' : ', ztra-tra_ctl(jn), ' muscl' 
235            tra_ctl(jn) = ztra 
236         ENDIF
237
238         ! II. Vertical advective fluxes
239         ! -----------------------------
240
241         ! First guess of the slope
242         ! interior values
243         DO jk = 2, jpkm1
244            zt1(:,:,jk) = tmask(:,:,jk) * ( trb(:,:,jk-1,jn) - trb(:,:,jk,jn) )
245         END DO
246         ! surface and bottom boundary conditions
247         zt1 (:,:, 1 ) = 0.e0 
248         zt1 (:,:,jpk) = 0.e0
249         ! Slopes
250         DO jk = 2, jpkm1
251            DO jj = 1, jpj
252               DO ji = 1, jpi
253                  ztp1(ji,jj,jk) =                    ( zt1(ji,jj,jk) + zt1(ji,jj,jk+1) )   &
254                     &           * ( 0.25 + SIGN( 0.25, zt1(ji,jj,jk) * zt1(ji,jj,jk+1) ) )
255               END DO
256            END DO
257         END DO
258
259         ! Slopes limitation
260         ! interior values
261         DO jk = 2, jpkm1
262            DO jj = 1, jpj
263               DO ji = 1, jpi
264                  ztp1(ji,jj,jk) = SIGN( 1., ztp1(ji,jj,jk) )   &
265                     &           * MIN(    ABS( ztp1(ji,jj,jk  ) ),   &
266                     &                  2.*ABS( zt1 (ji,jj,jk+1) ),   &
267                     &                  2.*ABS( zt1 (ji,jj,jk  ) ) )
268               END DO
269            END DO
270         END DO
271         ! surface values
272         ztp1(:,:,1) = 0. 
273         ! vertical advective flux
274         ! interior values
275         DO jk = 1, jpkm1
276            DO jj = 2, jpjm1     
277               DO ji = fs_2, fs_jpim1   ! vector opt.
278                  z2dtt = rdttra(jk) * FLOAT(ndttrc)
279                  zew = zwn(ji,jj,jk+1)
280                  z0w = SIGN( 0.5, zwn(ji,jj,jk+1) )
281                  zalpha = 0.5 + z0w
282                  zw  = z0w - 0.5 * zwn(ji,jj,jk+1)*z2dtt / fse3w(ji,jj,jk+1)
283                  zzt1 = trb(ji,jj,jk+1,jn) + zw*ztp1(ji,jj,jk+1)
284                  zzt2 = trb(ji,jj,jk  ,jn) + zw*ztp1(ji,jj,jk  )
285                  zt1(ji,jj,jk+1) = zew * ( zalpha * zzt1 + (1.-zalpha)*zzt2 )
286               END DO
287            END DO
288         END DO
289         ! surface values
290         IF( lk_dynspg_fsc .OR. lk_dynspg_fsc_tsk ) THEN        ! free surface-constant volume
291            zt1(:,:, 1 ) = zwn(:,:,1) * trb(:,:,1,jn)
292         ELSE                                                   ! rigid lid : flux set to zero
293            zt1(:,:, 1 ) = 0.e0
294         ENDIF
295
296         ! bottom values
297         zt1(:,:,jpk) = 0.e0
298
299         ! Compute & add the vertical advective trend
300
301         DO jk = 1, jpkm1
302            DO jj = 2, jpjm1     
303               DO ji = fs_2, fs_jpim1   ! vector opt.
304                  zbtr = 1. / fse3t(ji,jj,jk)
305                  ! horizontal advective trends
306                  ztra = - zbtr * ( zt1(ji,jj,jk) - zt1(ji,jj,jk+1) )
307                  ! add it to the general tracer trends
308                  tra(ji,jj,jk,jn) =  tra(ji,jj,jk,jn) + ztra
309#if defined key_trc_diatrd
310                  ! save the vertical advective trends computed as w gradz(T)
311                  trtrd(ji,jj,jk,jn,3) = ztra - trn(ji,jj,jk,jn) * hdivn(ji,jj,jk)
312#endif
313               END DO
314            END DO
315         END DO
316
317         IF(ln_ctl) THEN         ! print mean trends (used for debugging)
318            ztra = SUM( tra(2:nictle,2:njctle,1:jpkm1,jn) * tmask(2:nictle,2:njctle,1:jpkm1) )
319            WRITE(numout,*) ' trc/zad  - ',ctrcnm(jn),' : ', ztra-tra_ctl(jn), ' muscl' 
320            tra_ctl(jn) = ztra 
321         ENDIF
322
323      END DO
324
325END SUBROUTINE trc_adv_muscl
326
327#else
328   !!----------------------------------------------------------------------
329   !!   Default option                                         Empty module
330   !!----------------------------------------------------------------------
331CONTAINS
332   SUBROUTINE trc_adv_muscl( kt ) 
333      INTEGER, INTENT(in) :: kt
334      WRITE(*,*) 'trc_adv_muscl: You should not have seen this print! error?', kt
335   END SUBROUTINE trc_adv_muscl
336#endif
337
338   !!======================================================================
339END MODULE trcadv_muscl
Note: See TracBrowser for help on using the repository browser.