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

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

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