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

Last change on this file since 197 was 186, checked in by opalod, 20 years ago

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