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

Last change on this file since 941 was 941, checked in by cetlod, 16 years ago

phasing the passive tracer transport module to the new version of NEMO, see ticket 143
  • Property svn:executable set to *
File size: 16.5 KB
Line 
1MODULE trcadv_muscl2
2   !!==============================================================================
3   !!                       ***  MODULE  trcadv_muscl2  ***
4   !! Ocean passive tracers:  horizontal & vertical advective trend
5   !!==============================================================================
6#if defined key_top
7   !!----------------------------------------------------------------------
8   !!   'key_top'                                                TOP models
9   !!----------------------------------------------------------------------
10   !!   tra_adv_muscl2 : update the tracer trend with the horizontal
11   !!                    and vertical advection trends using MUSCL2 scheme
12   !!----------------------------------------------------------------------
13   !! * Modules used
14   USE oce_trc         ! ocean dynamics and active tracers variables
15   USE trc             ! ocean passive tracers variables
16   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
17   USE trcbbl          ! advective passive tracers in the BBL
18   USE prtctl_trc          ! Print control for debbuging
19
20   IMPLICIT NONE
21   PRIVATE
22
23   !! * Accessibility
24   PUBLIC trc_adv_muscl2        ! routine called by trcstp.F90
25
26   !! * Module variable
27   REAL(wp), DIMENSION(jpk) ::   &
28      rdttrc                     ! vertical profile of tracer time-step
29
30   !! * Substitutions
31#  include "top_substitute.h90"
32   !!----------------------------------------------------------------------
33   !!   TOP 1.0 , LOCEAN-IPSL (2005)
34   !! $Header$
35   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
36   !!----------------------------------------------------------------------
37
38CONTAINS
39
40   SUBROUTINE trc_adv_muscl2( kt )
41      !!----------------------------------------------------------------------
42      !!                   ***  ROUTINE trc_adv_muscl2  ***
43      !!
44      !! ** Purpose :   Compute the now trend due to total advection of passi-
45      !!      ve tracer using a MUSCL scheme (Monotone Upstream-
46      !!      Centered Scheme for Conservation Laws) and add it to the general
47      !!      tracer trend.
48      !!
49      !! ** Method : MUSCL scheme plus centered scheme at ocean boundaries
50      !!
51      !! ** Action : - update tra with the now advective tracer trends
52      !!             - save trends in trtrd ('key_trc_diatrd')
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      !!   9.0  !  03-04  (C. Ethe, G. Madec)  F90: Free form and module
61      !!----------------------------------------------------------------------
62      !! * modules used
63#if defined key_trcbbl_adv
64      USE oce_trc            , zun => ua,  &  ! use ua as workspace
65         &                     zvn => va      ! use va as workspace
66      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwn
67#else
68      USE oce_trc            , zun => un,  &  ! When no bbl, zun == un
69                               zvn => vn,  &  !              zvn == vn
70                               zwn => wn      !              zwn == wn
71#endif
72      !! * Arguments
73      INTEGER, INTENT( in ) ::   kt         ! ocean time-step
74
75      !! * Local declarations
76      INTEGER ::   ji, jj, jk,jn               ! dummy loop indices
77      REAL(wp), DIMENSION (jpi,jpj,jpk) ::   &
78         zt1, zt2, ztp1, ztp2
79      REAL(wp) ::   zu, zv, zw, zeu, zev, zew, zbtr, ztra
80      REAL(wp) ::   z0u, z0v, z0w
81      REAL(wp) ::   zzt1, zzt2, zalpha
82
83#if defined key_trc_diatrd
84      REAL(wp) ::   ztai, ztaj
85      REAL(wp) ::   zfui, zfvj
86#endif
87      CHARACTER (len=22) :: charout
88      !!----------------------------------------------------------------------
89
90      IF( kt == nittrc000 .AND. lwp ) THEN
91         WRITE(numout,*)
92         WRITE(numout,*) 'trc_adv_muscl2 : MUSCL2 advection scheme'
93         WRITE(numout,*) '~~~~~~~~~~~~~~~'
94         rdttrc(:) =  rdttra(:) * FLOAT(ndttrc)
95      ENDIF
96     
97#if defined key_trcbbl_adv       
98      ! Advective bottom boundary layer
99      ! -------------------------------
100      zun(:,:,:) = un (:,:,:) - u_trc_bbl(:,:,:)
101      zvn(:,:,:) = vn (:,:,:) - v_trc_bbl(:,:,:)
102      zwn(:,:,:) = wn (:,:,:) + w_trc_bbl( :,:,:)
103#endif
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  (changed sign)
126         CALL lbc_lnk( zt1, 'U', -1. )
127         CALL lbc_lnk( zt2, 'V', -1. )
128
129         ! Slopes
130         ! interior values
131         DO jk = 1, jpkm1
132            DO jj = 2, jpj
133               DO ji = fs_2, jpi   ! vector opt.
134                  ztp1(ji,jj,jk) =                    ( zt1(ji,jj,jk) + zt1(ji-1,jj  ,jk) )   &
135                     &           * ( 0.25 + SIGN( 0.25, zt1(ji,jj,jk) * zt1(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               END DO
139            END DO
140         END DO
141         ! bottom values
142         ztp1(:,:,jpk) = 0.e0 
143         ztp2(:,:,jpk) = 0.e0
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               END DO
159            END DO
160         END DO
161
162         ! Advection terms
163         ! interior values
164         DO jk = 1, jpkm1
165            DO jj = 2, jpjm1     
166               DO ji = fs_2, fs_jpim1   ! vector opt.
167                  ! volume fluxes
168#if ! defined key_zco
169                  zeu = e2u(ji,jj) * fse3u(ji,jj,jk) * zun(ji,jj,jk)
170                  zev = e1v(ji,jj) * fse3v(ji,jj,jk) * zvn(ji,jj,jk)
171#else
172                  zeu = e2u(ji,jj) * zun(ji,jj,jk)
173                  zev = e1v(ji,jj) * zvn(ji,jj,jk)
174#endif
175                  ! MUSCL fluxes
176                  z0u = SIGN( 0.5, zun(ji,jj,jk) )           
177                  zalpha = 0.5 - z0u
178                  zu  = z0u - 0.5 * zun(ji,jj,jk) * rdttrc(jk) / e1u(ji,jj)
179                  zzt1 = trb(ji+1,jj,jk,jn) + zu*ztp1(ji+1,jj,jk)
180                  zzt2 = trb(ji  ,jj,jk,jn) + zu*ztp1(ji  ,jj,jk)
181                  zt1(ji,jj,jk) = zeu * ( zalpha * zzt1 + (1.-zalpha) * zzt2 )
182
183                  z0v = SIGN( 0.5, zvn(ji,jj,jk) )           
184                  zalpha = 0.5 - z0v
185                  zv  = z0v - 0.5 * zvn(ji,jj,jk) * rdttrc(jk) / e2v(ji,jj)
186                  zzt1 = trb(ji,jj+1,jk,jn) + zv*ztp2(ji,jj+1,jk)
187                  zzt2 = trb(ji,jj  ,jk,jn) + zv*ztp2(ji,jj  ,jk)
188                  zt2(ji,jj,jk) = zev * ( zalpha * zzt1 + (1.-zalpha) * zzt2 )
189
190               END DO
191            END DO
192         END DO
193
194
195         DO jk = 1, jpkm1
196            DO jj = 2, jpjm1
197               DO ji = fs_2, fs_jpim1   ! vector opt.
198#if ! defined key_zco
199                  zev = e1v(ji,jj) * fse3v(ji,jj,jk)
200                  IF( umask(ji,jj,jk) == 0. ) THEN
201                     IF( zun(ji+1,jj,jk) > 0. .AND. ji /= jpi ) THEN
202                        zt1(ji+1,jj,jk) = e2u(ji+1,jj)* fse3u(ji+1,jj,jk)   &
203                           &            * zun(ji+1,jj,jk) * ( trb(ji+1,jj,jk,jn) + trb(ji+2,jj,jk,jn) ) * 0.5
204                     ENDIF
205                     IF( zun(ji-1,jj,jk) < 0. ) THEN
206                        zt1(ji-1,jj,jk) = e2u(ji-1,jj)* fse3u(ji-1,jj,jk)   &
207                           &            * zun(ji-1,jj,jk) * ( trb(ji-1,jj,jk,jn) + trb(ji  ,jj,jk,jn) ) * 0.5
208                     ENDIF
209                  ENDIF
210                  IF( vmask(ji,jj,jk) == 0. ) THEN
211                     IF( zvn(ji,jj+1,jk) > 0. .AND. jj /= jpj ) THEN
212                        zt2(ji,jj+1,jk) = e1v(ji,jj+1) * fse3v(ji,jj+1,jk)   &
213                           &            * zvn(ji,jj+1,jk) * ( trb(ji,jj+1,jk,jn) + trb(ji,jj+2,jk,jn) ) * 0.5
214                     ENDIF
215                     IF( zvn(ji,jj-1,jk) < 0. ) THEN
216                        zt2(ji,jj-1,jk) = e1v(ji,jj-1)* fse3v(ji,jj-1,jk)   &
217                           &            * zvn(ji,jj-1,jk) * ( trb(ji,jj-1,jk,jn) + trb(ji  ,jj,jk,jn) ) * 0.5
218                     ENDIF
219                  ENDIF
220
221#else
222                  IF( umask(ji,jj,jk) == 0. ) THEN
223                     IF( zun(ji+1,jj,jk) > 0. .AND. ji /= jpi ) THEN
224                        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
225                     ENDIF
226                     IF( zun(ji-1,jj,jk) < 0. ) THEN
227                        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
228                     ENDIF
229                  ENDIF
230                  IF( vmask(ji,jj,jk) == 0. ) THEN
231                     IF( zvn(ji,jj+1,jk) > 0. .AND. jj /= jpj ) THEN
232                        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
233                     ENDIF
234                     IF( zvn(ji,jj-1,jk) < 0. ) THEN
235                        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
236                     ENDIF
237                  ENDIF
238#endif
239               END DO
240            END DO
241         END DO
242
243         ! lateral boundary conditions on zt1, zt2   (changed sign)
244         CALL lbc_lnk( zt1, 'U', -1. )
245         CALL lbc_lnk( zt2, 'V', -1. )
246
247         ! Compute and add the horizontal advective trend
248
249         DO jk = 1, jpkm1
250            DO jj = 2, jpjm1     
251               DO ji = fs_2, fs_jpim1   ! vector opt.
252#if ! defined key_zco
253                  zbtr = 1. / ( e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk) )
254#else
255                  zbtr = 1. / ( e1t(ji,jj)*e2t(ji,jj) )
256#endif
257                  ! horizontal advective trends
258                  ztra = - zbtr * ( zt1(ji,jj,jk) - zt1(ji-1,jj  ,jk  )   &
259                     &            + zt2(ji,jj,jk) - zt2(ji  ,jj-1,jk  ) )
260                  ! add it to the general tracer trends
261                  tra(ji,jj,jk,jn) = tra(ji,jj,jk,jn) + ztra
262#if defined key_trc_diatrd
263                  ! recompute the trends in i- and j-direction as Uh gradh(T)
264#if ! defined key_zco
265                  zfui =  e2u(ji  ,jj) * fse3u(ji,  jj,jk) * un(ji,  jj,jk)   &
266                     & -  e2u(ji-1,jj) * fse3u(ji-1,jj,jk) * un(ji-1,jj,jk)
267                  zfvj =  e1v(ji,jj  ) * fse3v(ji,jj  ,jk) * vn(ji,jj  ,jk)   &
268                     & -  e1v(ji,jj-1) * fse3v(ji,jj-1,jk) * vn(ji,jj-1,jk)
269#   else
270                  zfui = e2u(ji  ,jj) * un(ji,  jj,jk)   &
271                     & - e2u(ji-1,jj) * un(ji-1,jj,jk)
272                  zfvj = e1v(ji,jj  ) * vn(ji,jj  ,jk)   &
273                     & - e1v(ji,jj-1) * vn(ji,jj-1,jk)
274#   endif
275                  ztai =-zbtr * (  zt1(ji,jj,jk) - zt1(ji-1,jj  ,jk) - trn(ji,jj,jk,jn) * zfui  )
276                  ztaj =-zbtr * (  zt2(ji,jj,jk) - zt2(ji  ,jj-1,jk) - trn(ji,jj,jk,jn) * zfvj  )
277                  ! save i- and j- advective trends computed as Uh gradh(T)
278                  IF (luttrd(jn)) trtrd(ji,jj,jk,ikeep(jn),1) = ztai
279                  IF (luttrd(jn)) trtrd(ji,jj,jk,ikeep(jn),2) = ztaj
280
281#endif
282               END DO
283            END DO
284         END DO
285      ENDDO
286
287      IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
288         WRITE(charout, FMT="('muscl2 - had')")
289         CALL prt_ctl_trc_info(charout)
290         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm,clinfo2='trd')
291      ENDIF
292
293      ! II. Vertical advective fluxes
294      ! -----------------------------
295
296      DO jn = 1, jptra
297
298         ! First guess of the slope
299         ! interior values
300         DO jk = 2, jpkm1
301            zt1(:,:,jk) = tmask(:,:,jk) * ( trb(:,:,jk-1,jn) - trb(:,:,jk,jn) )
302         END DO
303         ! surface and bottom boundary conditions
304         zt1 (:,:, 1 ) = 0.e0 
305         zt1 (:,:,jpk) = 0.e0
306
307         ! Slopes
308         DO jk = 2, jpkm1
309            DO jj = 1, jpj
310               DO ji = 1, jpi
311                  ztp1(ji,jj,jk) =                    ( zt1(ji,jj,jk) + zt1(ji,jj,jk+1) )   &
312                     &           * ( 0.25 + SIGN( 0.25, zt1(ji,jj,jk) * zt1(ji,jj,jk+1) ) )
313               END DO
314            END DO
315         END DO
316
317         ! Slopes limitation
318         ! interior values
319         DO jk = 2, jpkm1
320            DO jj = 1, jpj
321               DO ji = 1, jpi
322                  ztp1(ji,jj,jk) = SIGN( 1., ztp1(ji,jj,jk) )   &
323                     &           * MIN(    ABS( ztp1(ji,jj,jk  ) ),   &
324                     &                  2.*ABS( zt1 (ji,jj,jk+1) ),   &
325                     &                  2.*ABS( zt1 (ji,jj,jk  ) ) )
326               END DO
327            END DO
328         END DO
329
330         ! surface values
331         ztp1(:,:,1) = 0. 
332
333         ! vertical advective flux
334         ! interior values
335         DO jk = 1, jpkm1
336            DO jj = 2, jpjm1     
337               DO ji = fs_2, fs_jpim1   ! vector opt.
338                  zew = zwn(ji,jj,jk+1)
339                  z0w = SIGN( 0.5, zwn(ji,jj,jk+1) )
340                  zalpha = 0.5 + z0w
341                  zw  = z0w - 0.5 * zwn(ji,jj,jk+1)* rdttrc(jk)/ fse3w(ji,jj,jk+1)
342                  zzt1 = trb(ji,jj,jk+1,jn) + zw*ztp1(ji,jj,jk+1)
343                  zzt2 = trb(ji,jj,jk  ,jn) + zw*ztp1(ji,jj,jk  )
344                  zt1(ji,jj,jk+1) = zew * ( zalpha * zzt1 + (1.-zalpha)*zzt2 )
345               END DO
346            END DO
347         END DO
348         DO jk = 2, jpkm1
349            DO jj = 2, jpjm1
350               DO ji = fs_2, fs_jpim1   ! vector opt.
351                  IF( tmask(ji,jj,jk+1) == 0. ) THEN
352                     IF( zwn(ji,jj,jk) > 0. ) THEN
353                        zt1(ji,jj,jk) = zwn(ji,jj,jk) * ( trb(ji,jj,jk-1,jn) + trb(ji,jj,jk,jn) ) * 0.5
354                     ENDIF
355                  ENDIF
356               END DO
357            END DO
358         END DO
359
360         ! surface values
361         IF( lk_dynspg_rl ) THEN        ! rigid lid : flux set to zero
362            zt1(:,:, 1 ) = 0.e0
363         ELSE                           ! free surface
364            zt1(:,:, 1 ) = zwn(:,:,1) * trb(:,:,1,jn)
365         ENDIF
366
367         ! bottom values
368         zt1(:,:,jpk) = 0.e0
369
370         ! Compute & add the vertical advective trend
371
372         DO jk = 1, jpkm1
373            DO jj = 2, jpjm1     
374               DO ji = fs_2, fs_jpim1   ! vector opt.
375                  zbtr = 1. / fse3t(ji,jj,jk)
376                  ! horizontal advective trends
377                  ztra = - zbtr * ( zt1(ji,jj,jk) - zt1(ji,jj,jk+1) )
378                  ! add it to the general tracer trends
379                  tra(ji,jj,jk,jn) =  tra(ji,jj,jk,jn) + ztra
380#if defined key_trc_diatrd
381                  ! save the vertical advective trends computed as w gradz(T)
382                  IF (luttrd(jn)) trtrd(ji,jj,jk,ikeep(jn),3) = ztra - trn(ji,jj,jk,jn) * hdivn(ji,jj,jk)
383#endif
384               END DO
385            END DO
386         END DO
387
388      END DO
389
390      IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
391         WRITE(charout, FMT="('muscl2 - zad')")
392         CALL prt_ctl_trc_info(charout)
393         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm,clinfo2='trd')
394      ENDIF
395
396   END SUBROUTINE trc_adv_muscl2
397
398#else
399
400   !!----------------------------------------------------------------------
401   !!   Default option                                         Empty module
402   !!----------------------------------------------------------------------
403CONTAINS
404   SUBROUTINE trc_adv_muscl2( kt ) 
405      INTEGER, INTENT(in) :: kt
406      WRITE(*,*) 'trc_adv_muscl2: You should not have seen this print! error?', kt
407   END SUBROUTINE trc_adv_muscl2
408#endif
409
410   !!======================================================================
411END MODULE trcadv_muscl2
Note: See TracBrowser for help on using the repository browser.