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

Last change on this file since 247 was 247, checked in by opalod, 19 years ago

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