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

source: trunk/NEMO/TOP_SRC/TRP/trcadv_muscl2.F90 @ 204

Last change on this file since 204 was 204, checked in by opalod, 19 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: 16.4 KB
Line 
1MODULE trcadv_muscl2
2   !!==============================================================================
3   !!                       ***  MODULE  trcadv_muscl2  ***
4   !! Ocean passive tracers:  horizontal & vertical advective trend
5   !!==============================================================================
6#if defined key_passivetrc
7   !!----------------------------------------------------------------------
8   !!   tra_adv_muscl2 : update the tracer trend with the horizontal
9   !!                    and vertical advection trends using MUSCL2 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
17   IMPLICIT NONE
18   PRIVATE
19
20   !! * Accessibility
21   PUBLIC trc_adv_muscl2        ! 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_muscl2( kt )
36      !!----------------------------------------------------------------------
37      !!                   ***  ROUTINE trc_adv_muscl2  ***
38      !!
39      !! ** Purpose :   Compute the now trend due to total advection of passi-
40      !!      ve tracer using a MUSCL scheme (Monotone Upstream-
41      !!      Centered Scheme for Conservation Laws) and add it to the general
42      !!      tracer trend.
43      !!
44      !! ** Method : MUSCL scheme plus centered scheme at ocean boundaries
45      !!
46      !! ** Action : - update tra with the now advective tracer trends
47      !!             - save trends in trtrd ('key_trc_diatrd')
48      !!
49      !! References :               
50      !!      Estubier, A., and M. Levy, Notes Techn. Pole de Modelisation
51      !! IPSL, Sept. 2000 (http://www.lodyc.jussieu.fr/opa)
52      !!
53      !! History :
54      !!        !  06-00  (A.Estublier)  for passive tracers
55      !!   9.0  !  03-04  (C. Ethe, G. Madec)  F90: Free form and module
56      !!----------------------------------------------------------------------
57      !! * modules used
58#if defined key_trcbbl_adv
59      USE oce_trc            , zun => ua,  &  ! use ua as workspace
60         &                     zvn => va      ! use va as workspace
61      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwn
62#else
63      USE oce_trc            , zun => un,  &  ! When no bbl, zun == un
64                               zvn => vn,  &  !              zvn == vn
65                               zwn => wn      !              zwn == wn
66#endif
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      REAL(wp) ::   zu, zv, zw, zeu, zev, zew, zbtr, ztra
75      REAL(wp) ::   z0u, z0v, z0w
76      REAL(wp) ::   zzt1, zzt2, zalpha
77
78#if defined key_trc_diatrd
79      REAL(wp) ::   ztai, ztaj
80      REAL(wp) ::   zfui, zfvj
81#endif
82      !!----------------------------------------------------------------------
83
84      IF( kt == nittrc000 .AND. lwp ) THEN
85         WRITE(numout,*)
86         WRITE(numout,*) 'trc_adv_muscl2 : MUSCL2 advection scheme'
87         WRITE(numout,*) '~~~~~~~~~~~~~~~'
88         rdttrc(:) =  rdttra(:) * FLOAT(ndttrc)
89      ENDIF
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      DO jn = 1, jptra
101
102         ! I. Horizontal advective fluxes
103         ! ------------------------------
104
105         ! first guess of the slopes
106         ! interior values
107         DO jk = 1, jpkm1
108            DO jj = 1, jpjm1     
109               DO ji = 1, fs_jpim1   ! vector opt.
110                  zt1(ji,jj,jk) = umask(ji,jj,jk) * ( trb(ji+1,jj,jk,jn) - trb(ji,jj,jk,jn) )
111                  zt2(ji,jj,jk) = vmask(ji,jj,jk) * ( trb(ji,jj+1,jk,jn) - trb(ji,jj,jk,jn) )
112               END DO
113            END DO
114         END DO
115         ! bottom values
116         zt1(:,:,jpk) = 0.e0
117         zt2(:,:,jpk) = 0.e0
118
119         ! lateral boundary conditions on zt1, zt2  (changed sign)
120         CALL lbc_lnk( zt1, 'U', -1. )
121         CALL lbc_lnk( zt2, 'V', -1. )
122
123         ! Slopes
124         ! interior values
125         DO jk = 1, jpkm1
126            DO jj = 2, jpj
127               DO ji = fs_2, jpi   ! vector opt.
128                  ztp1(ji,jj,jk) =                    ( zt1(ji,jj,jk) + zt1(ji-1,jj  ,jk) )   &
129                     &           * ( 0.25 + SIGN( 0.25, zt1(ji,jj,jk) * zt1(ji-1,jj  ,jk) ) )
130                  ztp2(ji,jj,jk) =                    ( zt2(ji,jj,jk) + zt2(ji  ,jj-1,jk) )   &
131                     &           * ( 0.25 + SIGN( 0.25, zt2(ji,jj,jk) * zt2(ji  ,jj-1,jk) ) )
132               END DO
133            END DO
134         END DO
135         ! bottom values
136         ztp1(:,:,jpk) = 0.e0 
137         ztp2(:,:,jpk) = 0.e0
138
139         ! Slopes limitation
140         DO jk = 1, jpkm1
141            DO jj = 2, jpj
142               DO ji = fs_2, fs_jpim1   ! vector opt.
143                  ztp1(ji,jj,jk) = SIGN( 1., ztp1(ji,jj,jk) )   &
144                     &           * MIN(    ABS( ztp1(ji  ,jj,jk) ),   &
145                     &                  2.*ABS( zt1 (ji-1,jj,jk) ),   &
146                     &                  2.*ABS( zt1 (ji  ,jj,jk) ) )
147
148                  ztp2(ji,jj,jk) = SIGN( 1., ztp2(ji,jj,jk) )   &
149                     &           * MIN(    ABS( ztp2(ji,jj  ,jk) ),   &
150                     &                  2.*ABS( zt2 (ji,jj-1,jk) ),   &
151                     &                  2.*ABS( zt2 (ji,jj  ,jk) ) )
152               END DO
153            END DO
154         END DO
155
156         ! Advection terms
157         ! interior values
158         DO jk = 1, jpkm1
159            DO jj = 2, jpjm1     
160               DO ji = fs_2, fs_jpim1   ! vector opt.
161                  ! volume fluxes
162#if defined key_s_coord || defined key_partial_steps
163                  zeu = e2u(ji,jj) * fse3u(ji,jj,jk) * zun(ji,jj,jk)
164                  zev = e1v(ji,jj) * fse3v(ji,jj,jk) * zvn(ji,jj,jk)
165#else
166                  zeu = e2u(ji,jj) * zun(ji,jj,jk)
167                  zev = e1v(ji,jj) * zvn(ji,jj,jk)
168#endif
169                  ! MUSCL fluxes
170                  z0u = SIGN( 0.5, zun(ji,jj,jk) )           
171                  zalpha = 0.5 - z0u
172                  zu  = z0u - 0.5 * zun(ji,jj,jk) * rdttrc(jk) / e1u(ji,jj)
173                  zzt1 = trb(ji+1,jj,jk,jn) + zu*ztp1(ji+1,jj,jk)
174                  zzt2 = trb(ji  ,jj,jk,jn) + zu*ztp1(ji  ,jj,jk)
175                  zt1(ji,jj,jk) = zeu * ( zalpha * zzt1 + (1.-zalpha) * zzt2 )
176
177                  z0v = SIGN( 0.5, zvn(ji,jj,jk) )           
178                  zalpha = 0.5 - z0v
179                  zv  = z0v - 0.5 * zvn(ji,jj,jk) * rdttrc(jk) / e2v(ji,jj)
180                  zzt1 = trb(ji,jj+1,jk,jn) + zv*ztp2(ji,jj+1,jk)
181                  zzt2 = trb(ji,jj  ,jk,jn) + zv*ztp2(ji,jj  ,jk)
182                  zt2(ji,jj,jk) = zev * ( zalpha * zzt1 + (1.-zalpha) * zzt2 )
183
184               END DO
185            END DO
186         END DO
187
188
189         DO jk = 1, jpkm1
190            DO jj = 2, jpjm1
191               DO ji = fs_2, fs_jpim1   ! vector opt.
192#if defined key_s_coord || defined key_partial_steps
193                  zev = e1v(ji,jj) * fse3v(ji,jj,jk)
194                  IF( umask(ji,jj,jk) == 0. ) THEN
195                     IF( zun(ji+1,jj,jk) > 0. .AND. ji /= jpi ) THEN
196                        zt1(ji+1,jj,jk) = e2u(ji+1,jj)* fse3u(ji+1,jj,jk)   &
197                           &            * zun(ji+1,jj,jk) * ( trb(ji+1,jj,jk,jn) + trb(ji+2,jj,jk,jn) ) * 0.5
198                     ENDIF
199                     IF( zun(ji-1,jj,jk) < 0. ) THEN
200                        zt1(ji-1,jj,jk) = e2u(ji-1,jj)* fse3u(ji-1,jj,jk)   &
201                           &            * zun(ji-1,jj,jk) * ( trb(ji-1,jj,jk,jn) + trb(ji  ,jj,jk,jn) ) * 0.5
202                     ENDIF
203                  ENDIF
204                  IF( vmask(ji,jj,jk) == 0. ) THEN
205                     IF( zvn(ji,jj+1,jk) > 0. .AND. jj /= jpj ) THEN
206                        zt2(ji,jj+1,jk) = e1v(ji,jj+1) * fse3v(ji,jj+1,jk)   &
207                           &            * zvn(ji,jj+1,jk) * ( trb(ji,jj+1,jk,jn) + trb(ji,jj+2,jk,jn) ) * 0.5
208                     ENDIF
209                     IF( zvn(ji,jj-1,jk) < 0. ) THEN
210                        zt2(ji,jj-1,jk) = e1v(ji,jj-1)* fse3v(ji,jj-1,jk)   &
211                           &            * zvn(ji,jj-1,jk) * ( trb(ji,jj-1,jk,jn) + trb(ji  ,jj,jk,jn) ) * 0.5
212                     ENDIF
213                  ENDIF
214
215#else
216                  IF( umask(ji,jj,jk) == 0. ) THEN
217                     IF( zun(ji+1,jj,jk) > 0. .AND. ji /= jpi ) THEN
218                        zt1(ji+1,jj,jk) = e2u(ji+1,jj) * zun(ji+1,jj,jk) * ( trb(ji+1,jj,jk,jn) + trb(ji+2,jj,jk,jn) ) * 0.5
219                     ENDIF
220                     IF( zun(ji-1,jj,jk) < 0. ) THEN
221                        zt1(ji-1,jj,jk) = e2u(ji-1,jj) * zun(ji-1,jj,jk) * ( trb(ji-1,jj,jk,jn) + trb(ji  ,jj,jk,jn) ) * 0.5
222                     ENDIF
223                  ENDIF
224                  IF( vmask(ji,jj,jk) == 0. ) THEN
225                     IF( zvn(ji,jj+1,jk) > 0. .AND. jj /= jpj ) THEN
226                        zt2(ji,jj+1,jk) = e1v(ji,jj+1) * zvn(ji,jj+1,jk) * ( trb(ji,jj+1,jk,jn) + trb(ji,jj+2,jk,jn) ) * 0.5
227                     ENDIF
228                     IF( zvn(ji,jj-1,jk) < 0. ) THEN
229                        zt2(ji,jj-1,jk) = e1v(ji,jj-1) * zvn(ji,jj-1,jk) * ( trb(ji,jj-1,jk,jn) + trb(ji  ,jj,jk,jn) ) * 0.5
230                     ENDIF
231                  ENDIF
232#endif
233               END DO
234            END DO
235         END DO
236
237         ! lateral boundary conditions on zt1, zt2   (changed sign)
238         CALL lbc_lnk( zt1, 'U', -1. )
239         CALL lbc_lnk( zt2, 'V', -1. )
240
241         ! Compute and add the horizontal advective trend
242
243         DO jk = 1, jpkm1
244            DO jj = 2, jpjm1     
245               DO ji = fs_2, fs_jpim1   ! vector opt.
246#if defined key_s_coord || defined key_partial_steps
247                  zbtr = 1. / ( e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk) )
248#else
249                  zbtr = 1. / ( e1t(ji,jj)*e2t(ji,jj) )
250#endif
251                  ! horizontal advective trends
252                  ztra = - zbtr * ( zt1(ji,jj,jk) - zt1(ji-1,jj  ,jk  )   &
253                     &            + zt2(ji,jj,jk) - zt2(ji  ,jj-1,jk  ) )
254                  ! add it to the general tracer trends
255                  tra(ji,jj,jk,jn) = tra(ji,jj,jk,jn) + ztra
256#if defined key_trc_diatrd
257                  ! recompute the trends in i- and j-direction as Uh gradh(T)
258#   if defined key_s_coord || defined key_partial_steps
259                  zfui =  e2u(ji  ,jj) * fse3u(ji,  jj,jk) * un(ji,  jj,jk)   &
260                     & -  e2u(ji-1,jj) * fse3u(ji-1,jj,jk) * un(ji-1,jj,jk)
261                  zfvj =  e1v(ji,jj  ) * fse3v(ji,jj  ,jk) * vn(ji,jj  ,jk)   &
262                     & -  e1v(ji,jj-1) * fse3v(ji,jj-1,jk) * vn(ji,jj-1,jk)
263#   else
264                  zfui = e2u(ji  ,jj) * un(ji,  jj,jk)   &
265                     & - e2u(ji-1,jj) * un(ji-1,jj,jk)
266                  zfvj = e1v(ji,jj  ) * vn(ji,jj  ,jk)   &
267                     & - e1v(ji,jj-1) * vn(ji,jj-1,jk)
268#   endif
269                  ztai =-zbtr * (  zt1(ji,jj,jk) - zt1(ji-1,jj  ,jk) - trn(ji,jj,jk,jn) * zfui  )
270                  ztaj =-zbtr * (  zt2(ji,jj,jk) - zt2(ji  ,jj-1,jk) - trn(ji,jj,jk,jn) * zfvj  )
271                  ! save i- and j- advective trends computed as Uh gradh(T)
272                  trtrd(ji,jj,jk,jn,1) = ztai
273                  trtrd(ji,jj,jk,jn,2) = ztaj
274
275#endif
276               END DO
277            END DO
278         END DO
279
280         IF(l_ctl) THEN         ! print mean trends (used for debugging)
281            ztra = SUM( tra(2:nictl,2:njctl,1:jpkm1,jn) * tmask(2:nictl,2:njctl,1:jpkm1) )
282            WRITE(numout,*) ' trc/had  - ',ctrcnm(jn),' : ', ztra-tra_ctl(jn), ' muscl' 
283            tra_ctl(jn) = ztra 
284         ENDIF
285
286         ! II. Vertical advective fluxes
287         ! -----------------------------
288
289         ! First guess of the slope
290         ! interior values
291         DO jk = 2, jpkm1
292            zt1(:,:,jk) = tmask(:,:,jk) * ( trb(:,:,jk-1,jn) - trb(:,:,jk,jn) )
293         END DO
294         ! surface and bottom boundary conditions
295         zt1 (:,:, 1 ) = 0.e0 
296         zt1 (:,:,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               END DO
305            END DO
306         END DO
307
308         ! Slopes limitation
309         ! interior values
310         DO jk = 2, jpkm1
311            DO jj = 1, jpj
312               DO ji = 1, jpi
313                  ztp1(ji,jj,jk) = SIGN( 1., ztp1(ji,jj,jk) )   &
314                     &           * MIN(    ABS( ztp1(ji,jj,jk  ) ),   &
315                     &                  2.*ABS( zt1 (ji,jj,jk+1) ),   &
316                     &                  2.*ABS( zt1 (ji,jj,jk  ) ) )
317               END DO
318            END DO
319         END DO
320
321         ! surface values
322         ztp1(:,:,1) = 0. 
323
324         ! vertical advective flux
325         ! interior values
326         DO jk = 1, jpkm1
327            DO jj = 2, jpjm1     
328               DO ji = fs_2, fs_jpim1   ! vector opt.
329                  zew = zwn(ji,jj,jk+1)
330                  z0w = SIGN( 0.5, zwn(ji,jj,jk+1) )
331                  zalpha = 0.5 + z0w
332                  zw  = z0w - 0.5 * zwn(ji,jj,jk+1)* rdttrc(jk)/ fse3w(ji,jj,jk+1)
333                  zzt1 = trb(ji,jj,jk+1,jn) + zw*ztp1(ji,jj,jk+1)
334                  zzt2 = trb(ji,jj,jk  ,jn) + zw*ztp1(ji,jj,jk  )
335                  zt1(ji,jj,jk+1) = zew * ( zalpha * zzt1 + (1.-zalpha)*zzt2 )
336               END DO
337            END DO
338         END DO
339         DO jk = 2, jpkm1
340            DO jj = 2, jpjm1
341               DO ji = fs_2, fs_jpim1   ! vector opt.
342                  IF( tmask(ji,jj,jk+1) == 0. ) THEN
343                     IF( zwn(ji,jj,jk) > 0. ) THEN
344                        zt1(ji,jj,jk) = zwn(ji,jj,jk) * ( trb(ji,jj,jk-1,jn) + trb(ji,jj,jk,jn) ) * 0.5
345                     ENDIF
346                  ENDIF
347               END DO
348            END DO
349         END DO
350
351         ! surface values
352         IF( lk_dynspg_fsc .OR. lk_dynspg_fsc_tsk ) THEN        ! free surface-constant volume
353            zt1(:,:, 1 ) = zwn(:,:,1) * trb(:,:,1,jn)
354         ELSE                                                   ! rigid lid : flux set to zero
355            zt1(:,:, 1 ) = 0.e0
356         ENDIF
357
358         ! bottom values
359         zt1(:,:,jpk) = 0.e0
360
361         ! Compute & add the vertical advective trend
362
363         DO jk = 1, jpkm1
364            DO jj = 2, jpjm1     
365               DO ji = fs_2, fs_jpim1   ! vector opt.
366                  zbtr = 1. / fse3t(ji,jj,jk)
367                  ! horizontal advective trends
368                  ztra = - zbtr * ( zt1(ji,jj,jk) - zt1(ji,jj,jk+1) )
369                  ! add it to the general tracer trends
370                  tra(ji,jj,jk,jn) =  tra(ji,jj,jk,jn) + ztra
371#if defined key_trc_diatrd
372                  ! save the vertical advective trends computed as w gradz(T)
373                  trtrd(ji,jj,jk,jn,3) = ztra - trn(ji,jj,jk,jn) * hdivn(ji,jj,jk)
374#endif
375               END DO
376            END DO
377         END DO
378
379         IF(l_ctl) THEN         ! print mean trends (used for debugging)
380            ztra = SUM( tra(2:nictl,2:njctl,1:jpkm1,jn) * tmask(2:nictl,2:njctl,1:jpkm1) )
381            WRITE(numout,*) ' trc/zad  - ',ctrcnm(jn),' : ', ztra-tra_ctl(jn), ' muscl2' 
382            tra_ctl(jn) = ztra 
383         ENDIF
384
385      END DO
386
387   END SUBROUTINE trc_adv_muscl2
388
389#else
390
391   !!----------------------------------------------------------------------
392   !!   Default option                                         Empty module
393   !!----------------------------------------------------------------------
394CONTAINS
395   SUBROUTINE trc_adv_muscl2( kt ) 
396      INTEGER, INTENT(in) :: kt
397      WRITE(*,*) 'trc_adv_muscl2: You should not have seen this print! error?', kt
398   END SUBROUTINE trc_adv_muscl2
399#endif
400
401   !!======================================================================
402END MODULE trcadv_muscl2
Note: See TracBrowser for help on using the repository browser.