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 branches/dev_001_GM/NEMO/TOP_SRC/TRP – NEMO

source: branches/dev_001_GM/NEMO/TOP_SRC/TRP/trcadv_muscl2.F90 @ 772

Last change on this file since 772 was 772, checked in by gm, 16 years ago

dev_001_GM - change the name of cpp key to key_top, key_lobster, key_pisces, key_kriest and the corresponding lk_

  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
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 "passivetrc_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.