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

source: branches/dev_001_GM/NEMO/TOP_SRC/TRP/trcadv_muscl.F90 @ 776

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

dev_001_GM - passivetrc_substitute.h90 renamed top_substitute.h90 - compilation OK

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