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

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

update transport modules to take into account new trends organization, see ticket:248
  • Property svn:executable set to *
  • Property svn:keywords set to Id
File size: 16.6 KB
Line 
1MODULE trcadv_muscl
2   !!==============================================================================
3   !!                       ***  MODULE  trcadv_muscl  ***
4   !! Ocean passive tracers:  horizontal & vertical advective trend
5   !!==============================================================================
6#if defined key_top
7   !!----------------------------------------------------------------------
8   !!   'key_top'                                                TOP models
9   !!----------------------------------------------------------------------
10   !!   trc_adv_muscl : update the tracer trend with the horizontal
11   !!                   and vertical advection trends using MUSCL scheme
12   !!----------------------------------------------------------------------
13   !! * Modules used
14   USE oce_trc         ! ocean dynamics and active tracers variables
15   USE trp_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 lib_mpp
19   USE prtctl_trc      ! Print control for debbuging
20   USE trdmld_trc
21   USE trdmld_trc_oce          ! ocean variables trends
22
23   IMPLICIT NONE
24   PRIVATE
25
26   !! * Accessibility
27   PUBLIC trc_adv_muscl  ! routine called by trcstp.F90
28
29   !! * Substitutions
30#  include "top_substitute.h90"
31   !!----------------------------------------------------------------------
32   !!   TOP 1.0 , LOCEAN-IPSL (2005)
33   !! $Header: /home/opalod/NEMOCVSROOT/NEMO/TOP_SRC/TRP/trcadv_muscl.F90,v 1.13 2007/10/12 09:26:30 opalod Exp $
34   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
35   !!----------------------------------------------------------------------
36
37CONTAINS
38
39   SUBROUTINE trc_adv_muscl( kt )
40      !!----------------------------------------------------------------------
41      !!                    ***  ROUTINE trc_adv_muscl  ***
42      !!         
43      !! ** Purpose :   Compute the now trend due to total advection of any pas-
44      !!      sive tracer using a MUSCL scheme (Monotone Upstream-centered Scheme
45      !!      for Conservation Laws) and add it to the general tracer trend.
46      !!
47      !! ** Method  :
48      !!
49      !! ** Action  : - update tra with the now advective tracer trends
50      !!              - save trends ('key_trdmld_trc')
51      !!
52      !! References :               
53      !!      Estubier, A., and M. Levy, Notes Techn. Pole de Modelisation
54      !! IPSL, Sept. 2000 (http://www.lodyc.jussieu.fr/opa)
55      !!
56      !! History :
57      !!        !  06-00  (A.Estublier)  for passive tracers
58      !!   9.0  !  03-04  (C. Ethe, G. Madec)  F90: Free form and module
59      !!----------------------------------------------------------------------
60      !! * modules used
61#if defined key_trcbbl_adv
62      USE oce_trc            , zun => ua,  &  ! use ua as workspace
63         &                     zvn => va      ! use va as workspace
64      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwn
65#else
66      USE oce_trc            , zun => un,  &  ! When no bbl, zun == un
67                               zvn => vn,  &  !              zvn == vn
68                               zwn => wn      !              zwn == wn
69#endif
70
71      !! * Arguments
72      INTEGER, INTENT( in ) ::   kt         ! ocean time-step
73
74      !! * Local declarations
75      INTEGER ::   ji, jj, jk,jn            ! dummy loop indices
76      REAL(wp), DIMENSION (jpi,jpj,jpk) ::   &
77         zt1, zt2, ztp1, ztp2
78
79      REAL(wp) ::   zu, zv, zw, zeu, zev, zew, zbtr, ztra
80      REAL(wp) ::   z0u, z0v, z0w
81      REAL(wp) ::   zzt1, zzt2, zalpha, z2dtt
82      REAL(wp) ::   ztai, ztaj, zfui, zfvj
83      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrtrd
84      CHARACTER (len=22) :: charout
85      !!----------------------------------------------------------------------
86
87
88      IF( kt == nittrc000 .AND. lwp ) THEN
89         WRITE(numout,*)
90         WRITE(numout,*) 'trc_adv : MUSCL advection scheme'
91         WRITE(numout,*) '~~~~~~~'
92      ENDIF
93
94       IF( l_trdtrc ) ALLOCATE( ztrtrd(jpi,jpj,jpk) )
95
96#if defined key_trcbbl_adv
97      ! Advective bottom boundary layer
98      ! -------------------------------
99      zun(:,:,:) = un (:,:,:) - u_trc_bbl(:,:,:)
100      zvn(:,:,:) = vn (:,:,:) - v_trc_bbl(:,:,:)
101      zwn(:,:,:) = wn (:,:,:) + w_trc_bbl(:,:,:)
102#endif
103
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
126         CALL lbc_lnk( zt1, 'U', -1. )   
127         CALL lbc_lnk( zt2, 'V', -1. ) 
128
129
130         ! Slopes
131         ! interior values
132         DO jk = 1, jpkm1
133            DO jj = 2, jpj
134               DO ji = fs_2, jpi   ! vector opt.
135                  ztp1(ji,jj,jk) =                    ( zt1(ji,jj,jk) + zt1(ji-1,jj  ,jk) )   &
136                     &           * ( 0.25 + SIGN( 0.25, zt1(ji,jj,jk) * zt1(ji-1,jj  ,jk) ) )
137                  ztp2(ji,jj,jk) =                    ( zt2(ji,jj,jk) + zt2(ji  ,jj-1,jk) )   &
138                     &           * ( 0.25 + SIGN( 0.25, zt2(ji,jj,jk) * zt2(ji  ,jj-1,jk) ) )
139               END DO
140            END DO
141         END DO
142         ! bottom values
143         ztp1(:,:,jpk) = 0.e0
144         ztp2(:,:,jpk) = 0.e0
145
146         ! Slopes limitation
147         DO jk = 1, jpkm1
148            DO jj = 2, jpj
149               DO ji = fs_2, jpi   ! vector opt.
150                  ztp1(ji,jj,jk) = SIGN( 1., ztp1(ji,jj,jk) )   &
151                     &           * MIN(    ABS( ztp1(ji  ,jj,jk) ),   &
152                     &                  2.*ABS( zt1 (ji-1,jj,jk) ),   &
153                     &                  2.*ABS( zt1 (ji  ,jj,jk) ) )
154
155                  ztp2(ji,jj,jk) = SIGN( 1., ztp2(ji,jj,jk) )   &
156                     &           * MIN(    ABS( ztp2(ji,jj  ,jk) ),   &
157                     &                  2.*ABS( zt2 (ji,jj-1,jk) ),   &
158                     &                  2.*ABS( zt2 (ji,jj  ,jk) ) )
159
160               END DO
161            END DO
162         END DO
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_zco
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                  z2dtt = rdttra(jk) * FLOAT(ndttrc)
179                  z0u = SIGN( 0.5, zun(ji,jj,jk) )           
180                  zalpha = 0.5 - z0u
181                  zu  = z0u - 0.5 * zun(ji,jj,jk) * z2dtt / e1u(ji,jj)
182                  zzt1 = trb(ji+1,jj,jk,jn) + zu*ztp1(ji+1,jj,jk)
183                  zzt2 = trb(ji  ,jj,jk,jn) + zu*ztp1(ji  ,jj,jk)
184                  zt1(ji,jj,jk) = zeu * ( zalpha * zzt1 + (1.-zalpha) * zzt2 )
185                  z0v = SIGN( 0.5, zvn(ji,jj,jk) )           
186                  zalpha = 0.5 - z0v
187                  zv  = z0v - 0.5 * zvn(ji,jj,jk) * z2dtt / 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_zco
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                  IF (luttrd(jn)) trtrd(ji,jj,jk,ikeep(jn),1) = ztai
231                  IF (luttrd(jn)) trtrd(ji,jj,jk,ikeep(jn),2) = ztaj
232#endif
233
234               END DO
235            END DO
236         END DO
237
238         ! 3. Save the horizontal advective trends for diagnostics
239         ! -------------------------------------------------------
240!CDIR BEGIN COLLAPSE
241         TRDTRC_XY : IF( l_trdtrc ) THEN
242
243            ! 3.1) Passive tracer ZONAL advection trends
244            DO jk = 1, jpkm1
245               DO jj = 2, jpjm1
246                  DO ji = fs_2, fs_jpim1   ! vector opt.
247#if ! defined key_zco
248                     zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
249                     zfui = e2u(ji  ,jj) * fse3u(ji,  jj,jk) * zun(ji,  jj,jk)   &
250                        & - e2u(ji-1,jj) * fse3u(ji-1,jj,jk) * zun(ji-1,jj,jk)
251#else
252                     zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) )
253                     zfui = e2u(ji  ,jj) * zun(ji,  jj,jk)  - e2u(ji-1,jj) * zun(ji-1,jj,jk)
254#endif
255                     ! recompute the trends in i- direction as Uh gradh(T)
256                     ztrtrd(ji,jj,jk) = - zbtr*( zt1(ji,jj,jk) - zt1(ji-1,jj,jk) - trn(ji,jj,jk,jn)*zfui )
257                  END DO
258               END DO
259            END DO
260
261            IF (luttrd(jn)) CALL trd_mod_trc( ztrtrd, jn, jptrc_trd_xad, kt )
262
263            ! 3.2)  Passive tracer MERIDIONAL advection trends
264            DO jk = 1, jpkm1
265               DO jj = 2, jpjm1
266                  DO ji = fs_2, fs_jpim1   ! vector opt.
267                     ! recompute the trends in i- and j-direction as Uh gradh(T)
268#if ! defined key_zco
269                     zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
270                     zfvj = e1v(ji,jj  ) * fse3v(ji,jj  ,jk) * zvn(ji,jj  ,jk)   &
271                        & - e1v(ji,jj-1) * fse3v(ji,jj-1,jk) * zvn(ji,jj-1,jk)
272#else
273                     zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) )
274                     zfvj = e1v(ji,jj  ) * zvn(ji,jj  ,jk) - e1v(ji,jj-1) * zvn(ji,jj-1,jk)
275#endif
276                     ztrtrd(ji,jj,jk) = - zbtr*( zt2(ji,jj,jk) - zt2(ji,jj-1,jk) - trn(ji,jj,jk,jn)*zfvj )
277                  END DO
278               END DO
279            END DO
280
281            IF (luttrd(jn)) CALL trd_mod_trc( ztrtrd, jn, jptrc_trd_yad, kt )
282
283         ENDIF TRDTRC_XY
284!CDIR END
285
286      ENDDO
287
288      IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
289         WRITE(charout, FMT="('muscl - had')")
290         CALL prt_ctl_trc_info(charout)
291         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm,clinfo2='trd')
292      ENDIF
293
294         ! II. Vertical advective fluxes
295         ! -----------------------------
296
297      DO jn = 1, jptra
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         ! Slopes
307         DO jk = 2, jpkm1
308            DO jj = 1, jpj
309               DO ji = 1, jpi
310                  ztp1(ji,jj,jk) =                    ( zt1(ji,jj,jk) + zt1(ji,jj,jk+1) )   &
311                     &           * ( 0.25 + SIGN( 0.25, zt1(ji,jj,jk) * zt1(ji,jj,jk+1) ) )
312               END DO
313            END DO
314         END DO
315
316         ! Slopes limitation
317         ! interior values
318         DO jk = 2, jpkm1
319            DO jj = 1, jpj
320               DO ji = 1, jpi
321                  ztp1(ji,jj,jk) = SIGN( 1., ztp1(ji,jj,jk) )   &
322                     &           * MIN(    ABS( ztp1(ji,jj,jk  ) ),   &
323                     &                  2.*ABS( zt1 (ji,jj,jk+1) ),   &
324                     &                  2.*ABS( zt1 (ji,jj,jk  ) ) )
325               END DO
326            END DO
327         END DO
328         ! surface values
329         ztp1(:,:,1) = 0. 
330         ! vertical advective flux
331         ! interior values
332         DO jk = 1, jpkm1
333            DO jj = 2, jpjm1     
334               DO ji = fs_2, fs_jpim1   ! vector opt.
335                  z2dtt = rdttra(jk) * FLOAT(ndttrc)
336                  zew = zwn(ji,jj,jk+1)
337                  z0w = SIGN( 0.5, zwn(ji,jj,jk+1) )
338                  zalpha = 0.5 + z0w
339                  zw  = z0w - 0.5 * zwn(ji,jj,jk+1)*z2dtt / fse3w(ji,jj,jk+1)
340                  zzt1 = trb(ji,jj,jk+1,jn) + zw*ztp1(ji,jj,jk+1)
341                  zzt2 = trb(ji,jj,jk  ,jn) + zw*ztp1(ji,jj,jk  )
342                  zt1(ji,jj,jk+1) = zew * ( zalpha * zzt1 + (1.-zalpha)*zzt2 )
343               END DO
344            END DO
345         END DO
346         ! surface values
347         IF( lk_dynspg_rl ) THEN        ! rigid lid : flux set to zero
348            zt1(:,:, 1 ) = 0.e0
349         ELSE                           ! free surface
350            zt1(:,:, 1 ) = zwn(:,:,1) * trb(:,:,1,jn)
351         ENDIF
352
353         ! bottom values
354         zt1(:,:,jpk) = 0.e0
355
356         ! Compute & add the vertical advective trend
357
358         DO jk = 1, jpkm1
359            DO jj = 2, jpjm1     
360               DO ji = fs_2, fs_jpim1   ! vector opt.
361                  zbtr = 1. / fse3t(ji,jj,jk)
362                  ! horizontal advective trends
363                  ztra = - zbtr * ( zt1(ji,jj,jk) - zt1(ji,jj,jk+1) )
364                  ! add it to the general tracer trends
365                  tra(ji,jj,jk,jn) =  tra(ji,jj,jk,jn) + ztra
366#if defined key_trc_diatrd
367                  ! save the vertical advective trends computed as w gradz(T)
368                  IF (luttrd(jn)) trtrd(ji,jj,jk,ikeep(jn),3) = ztra - trn(ji,jj,jk,jn) * hdivn(ji,jj,jk)
369#endif
370
371               END DO
372            END DO
373         END DO
374
375         ! 3. Save the vertical advective trends for diagnostic
376         ! ----------------------------------------------------
377!CDIR BEGIN COLLAPSE
378         TRDTRC_Z : IF( l_trdtrc )THEN
379
380            ! Compute T/S vertical advection trends
381            DO jk = 1, jpkm1
382               DO jj = 2, jpjm1
383                  DO ji = fs_2, fs_jpim1   ! vector opt.
384                     zbtr = 1. / fse3t(ji,jj,jk)
385                     ! horizontal advective trends
386                     ztra = - zbtr * ( zt1(ji,jj,jk) - zt1(ji,jj,jk+1) )
387                     ! save the vertical advective trends computed as w gradz(T)
388                     ztrtrd(ji,jj,jk) = ztra - trn(ji,jj,jk,jn) * hdivn(ji,jj,jk)
389                  END DO
390               END DO
391            END DO
392
393            IF (luttrd(jn)) CALL trd_mod_trc(ztrtrd, jn, jptrc_trd_zad, kt)
394
395         END IF TRDTRC_Z
396!CDIR END
397      END DO
398
399      IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
400         WRITE(charout, FMT="('muscl - zad')")
401         CALL prt_ctl_trc_info(charout)
402         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm,clinfo2='trd')
403      ENDIF
404
405END SUBROUTINE trc_adv_muscl
406
407#else
408   !!----------------------------------------------------------------------
409   !!   Default option                                         Empty module
410   !!----------------------------------------------------------------------
411CONTAINS
412   SUBROUTINE trc_adv_muscl( kt ) 
413      INTEGER, INTENT(in) :: kt
414      WRITE(*,*) 'trc_adv_muscl: You should not have seen this print! error?', kt
415   END SUBROUTINE trc_adv_muscl
416#endif
417
418   !!======================================================================
419END MODULE trcadv_muscl
Note: See TracBrowser for help on using the repository browser.