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 @ 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: 13.9 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 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
21   IMPLICIT NONE
22   PRIVATE
23
24   !! * Accessibility
25   PUBLIC trc_adv_muscl  ! routine called by trcstp.F90
26
27   !! * Substitutions
28#  include "top_substitute.h90"
29   !!----------------------------------------------------------------------
30   !!   TOP 1.0 , LOCEAN-IPSL (2005)
31   !! $Header: /home/opalod/NEMOCVSROOT/NEMO/TOP_SRC/TRP/trcadv_muscl.F90,v 1.13 2007/10/12 09:26:30 opalod Exp $
32   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
33   !!----------------------------------------------------------------------
34
35CONTAINS
36
37   SUBROUTINE trc_adv_muscl( kt )
38      !!----------------------------------------------------------------------
39      !!                    ***  ROUTINE trc_adv_muscl  ***
40      !!         
41      !! ** Purpose :   Compute the now trend due to total advection of any pas-
42      !!      sive tracer using a MUSCL scheme (Monotone Upstream-centered Scheme
43      !!      for Conservation Laws) and add it to the general tracer trend.
44      !!
45      !! ** Method  :
46      !!
47      !! ** Action  : - update tra with the now advective tracer trends
48      !!              - save trends in trtrd ('key_trc_diatrd')
49      !!
50      !! References :               
51      !!      Estubier, A., and M. Levy, Notes Techn. Pole de Modelisation
52      !! IPSL, Sept. 2000 (http://www.lodyc.jussieu.fr/opa)
53      !!
54      !! History :
55      !!        !  06-00  (A.Estublier)  for passive tracers
56      !!   9.0  !  03-04  (C. Ethe, G. Madec)  F90: Free form and module
57      !!----------------------------------------------------------------------
58      !! * modules used
59#if defined key_trcbbl_adv
60      USE oce_trc            , zun => ua,  &  ! use ua as workspace
61         &                     zvn => va      ! use va as workspace
62      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwn
63#else
64      USE oce_trc            , zun => un,  &  ! When no bbl, zun == un
65                               zvn => vn,  &  !              zvn == vn
66                               zwn => wn      !              zwn == wn
67#endif
68
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
77      REAL(wp) ::   zu, zv, zw, zeu, zev, zew, zbtr, ztra
78      REAL(wp) ::   z0u, z0v, z0w
79      REAL(wp) ::   zzt1, zzt2, zalpha, z2dtt
80#if defined key_trc_diatrd
81      REAL(wp) ::   ztai, ztaj
82      REAL(wp) ::   zfui, zfvj
83#endif
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 
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#if defined key_trc_diatrd
108        DO jk = 1,jpk
109           DO jj = 1,jpj
110              DO ji = 1,jpi
111                 IF (luttrd(jn)) trtrd(ji,jj,jk,ikeep(jn),1) = 0.
112                 IF (luttrd(jn)) trtrd(ji,jj,jk,ikeep(jn),2) = 0.
113                 IF (luttrd(jn)) trtrd(ji,jj,jk,ikeep(jn),3) = 0.
114              END DO
115            END DO
116          END DO
117#endif
118         ! I. Horizontal advective fluxes
119         ! ------------------------------
120
121         ! first guess of the slopes
122         ! interior values
123         DO jk = 1, jpkm1
124            DO jj = 1, jpjm1     
125               DO ji = 1, fs_jpim1   ! vector opt.
126                  zt1(ji,jj,jk) = umask(ji,jj,jk) * ( trb(ji+1,jj,jk,jn) - trb(ji,jj,jk,jn) )
127                  zt2(ji,jj,jk) = vmask(ji,jj,jk) * ( trb(ji,jj+1,jk,jn) - trb(ji,jj,jk,jn) )
128               END DO
129            END DO
130         END DO
131         ! bottom values
132         zt1(:,:,jpk) = 0.e0
133         zt2(:,:,jpk) = 0.e0
134
135         ! lateral boundary conditions on zt1, zt2
136         CALL lbc_lnk( zt1, 'U', -1. )   
137         CALL lbc_lnk( zt2, 'V', -1. ) 
138
139
140         ! Slopes
141         ! interior values
142         DO jk = 1, jpkm1
143            DO jj = 2, jpj
144               DO ji = fs_2, jpi   ! vector opt.
145                  ztp1(ji,jj,jk) =                    ( zt1(ji,jj,jk) + zt1(ji-1,jj  ,jk) )   &
146                     &           * ( 0.25 + SIGN( 0.25, zt1(ji,jj,jk) * zt1(ji-1,jj  ,jk) ) )
147                  ztp2(ji,jj,jk) =                    ( zt2(ji,jj,jk) + zt2(ji  ,jj-1,jk) )   &
148                     &           * ( 0.25 + SIGN( 0.25, zt2(ji,jj,jk) * zt2(ji  ,jj-1,jk) ) )
149               END DO
150            END DO
151         END DO
152         ! bottom values
153         ztp1(:,:,jpk) = 0.e0
154         ztp2(:,:,jpk) = 0.e0
155
156         ! Slopes limitation
157         DO jk = 1, jpkm1
158            DO jj = 2, jpj
159               DO ji = fs_2, jpi   ! vector opt.
160                  ztp1(ji,jj,jk) = SIGN( 1., ztp1(ji,jj,jk) )   &
161                     &           * MIN(    ABS( ztp1(ji  ,jj,jk) ),   &
162                     &                  2.*ABS( zt1 (ji-1,jj,jk) ),   &
163                     &                  2.*ABS( zt1 (ji  ,jj,jk) ) )
164
165                  ztp2(ji,jj,jk) = SIGN( 1., ztp2(ji,jj,jk) )   &
166                     &           * MIN(    ABS( ztp2(ji,jj  ,jk) ),   &
167                     &                  2.*ABS( zt2 (ji,jj-1,jk) ),   &
168                     &                  2.*ABS( zt2 (ji,jj  ,jk) ) )
169
170               END DO
171            END DO
172         END DO
173
174         ! Advection terms
175         ! interior values
176         DO jk = 1, jpkm1
177            DO jj = 2, jpjm1     
178               DO ji = fs_2, fs_jpim1   ! vector opt.
179                  ! volume fluxes
180#if ! defined key_zco
181                  zeu = e2u(ji,jj) * fse3u(ji,jj,jk) * zun(ji,jj,jk)
182                  zev = e1v(ji,jj) * fse3v(ji,jj,jk) * zvn(ji,jj,jk)
183#else
184                  zeu = e2u(ji,jj) * zun(ji,jj,jk)
185                  zev = e1v(ji,jj) * zvn(ji,jj,jk)
186#endif
187                  ! MUSCL fluxes
188                  z2dtt = rdttra(jk) * FLOAT(ndttrc)
189                  z0u = SIGN( 0.5, zun(ji,jj,jk) )           
190                  zalpha = 0.5 - z0u
191                  zu  = z0u - 0.5 * zun(ji,jj,jk) * z2dtt / e1u(ji,jj)
192                  zzt1 = trb(ji+1,jj,jk,jn) + zu*ztp1(ji+1,jj,jk)
193                  zzt2 = trb(ji  ,jj,jk,jn) + zu*ztp1(ji  ,jj,jk)
194                  zt1(ji,jj,jk) = zeu * ( zalpha * zzt1 + (1.-zalpha) * zzt2 )
195                  z0v = SIGN( 0.5, zvn(ji,jj,jk) )           
196                  zalpha = 0.5 - z0v
197                  zv  = z0v - 0.5 * zvn(ji,jj,jk) * z2dtt / e2v(ji,jj)
198                  zzt1 = trb(ji,jj+1,jk,jn) + zv*ztp2(ji,jj+1,jk)
199                  zzt2 = trb(ji,jj  ,jk,jn) + zv*ztp2(ji,jj  ,jk)
200                  zt2(ji,jj,jk) = zev * ( zalpha * zzt1 + (1.-zalpha) * zzt2 )
201               END DO
202            END DO
203         END DO
204
205         ! lateral boundary conditions on zt1, zt2 (changed sign)
206         CALL lbc_lnk( zt1, 'U', -1. ) 
207         CALL lbc_lnk( zt2, 'V', -1. ) 
208
209         ! Compute and add the horizontal advective trend
210
211         DO jk = 1, jpkm1
212            DO jj = 2, jpjm1     
213               DO ji = fs_2, fs_jpim1   ! vector opt.
214#if ! defined key_zco
215                  zbtr = 1. / ( e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk) )
216#else
217                  zbtr = 1. / ( e1t(ji,jj)*e2t(ji,jj) )
218#endif
219                  ! horizontal advective trends
220                  ztra = - zbtr * ( zt1(ji,jj,jk) - zt1(ji-1,jj  ,jk  )   &
221                     &            + zt2(ji,jj,jk) - zt2(ji  ,jj-1,jk  ) )
222                  ! add it to the general tracer trends
223                  tra(ji,jj,jk,jn) = tra(ji,jj,jk,jn) + ztra
224#if defined key_trc_diatrd
225                  ! recompute the trends in i- and j-direction as Uh gradh(T)
226#if ! defined key_zco
227                  zfui =  e2u(ji  ,jj) * fse3u(ji,  jj,jk) * un(ji,  jj,jk)   &
228                     & -  e2u(ji-1,jj) * fse3u(ji-1,jj,jk) * un(ji-1,jj,jk)
229                  zfvj =  e1v(ji,jj  ) * fse3v(ji,jj  ,jk) * vn(ji,jj  ,jk)   &
230                     & -  e1v(ji,jj-1) * fse3v(ji,jj-1,jk) * vn(ji,jj-1,jk)
231#   else
232                  zfui = e2u(ji  ,jj) * un(ji,  jj,jk)   &
233                     & - e2u(ji-1,jj) * un(ji-1,jj,jk)
234                  zfvj = e1v(ji,jj  ) * vn(ji,jj  ,jk)   &
235                     & - e1v(ji,jj-1) * vn(ji,jj-1,jk)
236#   endif
237                  ztai =-zbtr * (  zt1(ji,jj,jk) - zt1(ji-1,jj  ,jk) - trn(ji,jj,jk,jn) * zfui  )
238                  ztaj =-zbtr * (  zt2(ji,jj,jk) - zt2(ji  ,jj-1,jk) - trn(ji,jj,jk,jn) * zfvj  )
239                  ! save i- and j- advective trends computed as Uh gradh(T)
240                  IF (luttrd(jn)) trtrd(ji,jj,jk,ikeep(jn),1) = ztai
241                  IF (luttrd(jn)) trtrd(ji,jj,jk,ikeep(jn),2) = ztaj
242#endif
243               END DO
244            END DO
245         END DO
246      ENDDO
247
248      IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
249         WRITE(charout, FMT="('muscl - had')")
250         CALL prt_ctl_trc_info(charout)
251         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm,clinfo2='trd')
252      ENDIF
253
254         ! II. Vertical advective fluxes
255         ! -----------------------------
256
257      DO jn = 1, jptra
258         ! First guess of the slope
259         ! interior values
260         DO jk = 2, jpkm1
261            zt1(:,:,jk) = tmask(:,:,jk) * ( trb(:,:,jk-1,jn) - trb(:,:,jk,jn) )
262         END DO
263         ! surface and bottom boundary conditions
264         zt1 (:,:, 1 ) = 0.e0 
265         zt1 (:,:,jpk) = 0.e0
266         ! Slopes
267         DO jk = 2, jpkm1
268            DO jj = 1, jpj
269               DO ji = 1, jpi
270                  ztp1(ji,jj,jk) =                    ( zt1(ji,jj,jk) + zt1(ji,jj,jk+1) )   &
271                     &           * ( 0.25 + SIGN( 0.25, zt1(ji,jj,jk) * zt1(ji,jj,jk+1) ) )
272               END DO
273            END DO
274         END DO
275
276         ! Slopes limitation
277         ! interior values
278         DO jk = 2, jpkm1
279            DO jj = 1, jpj
280               DO ji = 1, jpi
281                  ztp1(ji,jj,jk) = SIGN( 1., ztp1(ji,jj,jk) )   &
282                     &           * MIN(    ABS( ztp1(ji,jj,jk  ) ),   &
283                     &                  2.*ABS( zt1 (ji,jj,jk+1) ),   &
284                     &                  2.*ABS( zt1 (ji,jj,jk  ) ) )
285               END DO
286            END DO
287         END DO
288         ! surface values
289         ztp1(:,:,1) = 0. 
290         ! vertical advective flux
291         ! interior values
292         DO jk = 1, jpkm1
293            DO jj = 2, jpjm1     
294               DO ji = fs_2, fs_jpim1   ! vector opt.
295                  z2dtt = rdttra(jk) * FLOAT(ndttrc)
296                  zew = zwn(ji,jj,jk+1)
297                  z0w = SIGN( 0.5, zwn(ji,jj,jk+1) )
298                  zalpha = 0.5 + z0w
299                  zw  = z0w - 0.5 * zwn(ji,jj,jk+1)*z2dtt / fse3w(ji,jj,jk+1)
300                  zzt1 = trb(ji,jj,jk+1,jn) + zw*ztp1(ji,jj,jk+1)
301                  zzt2 = trb(ji,jj,jk  ,jn) + zw*ztp1(ji,jj,jk  )
302                  zt1(ji,jj,jk+1) = zew * ( zalpha * zzt1 + (1.-zalpha)*zzt2 )
303               END DO
304            END DO
305         END DO
306         ! surface values
307         IF( lk_dynspg_rl ) THEN        ! rigid lid : flux set to zero
308            zt1(:,:, 1 ) = 0.e0
309         ELSE                           ! free surface
310            zt1(:,:, 1 ) = zwn(:,:,1) * trb(:,:,1,jn)
311         ENDIF
312
313         ! bottom values
314         zt1(:,:,jpk) = 0.e0
315
316         ! Compute & add the vertical advective trend
317
318         DO jk = 1, jpkm1
319            DO jj = 2, jpjm1     
320               DO ji = fs_2, fs_jpim1   ! vector opt.
321                  zbtr = 1. / fse3t(ji,jj,jk)
322                  ! horizontal advective trends
323                  ztra = - zbtr * ( zt1(ji,jj,jk) - zt1(ji,jj,jk+1) )
324                  ! add it to the general tracer trends
325                  tra(ji,jj,jk,jn) =  tra(ji,jj,jk,jn) + ztra
326#if defined key_trc_diatrd
327                  ! save the vertical advective trends computed as w gradz(T)
328                  IF (luttrd(jn)) trtrd(ji,jj,jk,ikeep(jn),3) = ztra - trn(ji,jj,jk,jn) * hdivn(ji,jj,jk)
329#endif
330               END DO
331            END DO
332         END DO
333
334      END DO
335
336      IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
337         WRITE(charout, FMT="('muscl - zad')")
338         CALL prt_ctl_trc_info(charout)
339         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm,clinfo2='trd')
340      ENDIF
341
342END SUBROUTINE trc_adv_muscl
343
344#else
345   !!----------------------------------------------------------------------
346   !!   Default option                                         Empty module
347   !!----------------------------------------------------------------------
348CONTAINS
349   SUBROUTINE trc_adv_muscl( kt ) 
350      INTEGER, INTENT(in) :: kt
351      WRITE(*,*) 'trc_adv_muscl: You should not have seen this print! error?', kt
352   END SUBROUTINE trc_adv_muscl
353#endif
354
355   !!======================================================================
356END MODULE trcadv_muscl
Note: See TracBrowser for help on using the repository browser.