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

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

nemo_v1_update_022 : CE + RB + CT : add print control possibility
  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 13.3 KB
Line 
1MODULE trcadv_muscl
2   !!==============================================================================
3   !!                       ***  MODULE  trcadv_muscl  ***
4   !! Ocean passive tracers:  horizontal & vertical advective trend
5   !!==============================================================================
6#if defined key_passivetrc
7   !!----------------------------------------------------------------------
8   !!   trc_adv_muscl : update the tracer trend with the horizontal
9   !!                   and vertical advection trends using MUSCL 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   USE lib_mpp
17   USE prtctl_trc      ! Print control for debbuging
18
19   IMPLICIT NONE
20   PRIVATE
21
22   !! * Accessibility
23   PUBLIC trc_adv_muscl  ! routine called by trcstp.F90
24
25   !! * Substitutions
26#  include "passivetrc_substitute.h90"
27   !!----------------------------------------------------------------------
28   !!   OPA 9.0 , LODYC-IPSL (2003)
29   !!----------------------------------------------------------------------
30
31CONTAINS
32
33   SUBROUTINE trc_adv_muscl( kt )
34      !!----------------------------------------------------------------------
35      !!                    ***  ROUTINE trc_adv_muscl  ***
36      !!         
37      !! ** Purpose :   Compute the now trend due to total advection of any pas-
38      !!      sive tracer using a MUSCL scheme (Monotone Upstream-centered Scheme
39      !!      for Conservation Laws) and add it to the general tracer trend.
40      !!
41      !! ** Method  :
42      !!
43      !! ** Action  : - update tra with the now advective tracer trends
44      !!              - save trends in trtrd ('key_trc_diatrd')
45      !!
46      !! References :               
47      !!      Estubier, A., and M. Levy, Notes Techn. Pole de Modelisation
48      !! IPSL, Sept. 2000 (http://www.lodyc.jussieu.fr/opa)
49      !!
50      !! History :
51      !!        !  06-00  (A.Estublier)  for passive tracers
52      !!   9.0  !  03-04  (C. Ethe, G. Madec)  F90: Free form and module
53      !!----------------------------------------------------------------------
54      !! * modules used
55#if defined key_trcbbl_adv
56      USE oce_trc            , zun => ua,  &  ! use ua as workspace
57         &                     zvn => va      ! use va as workspace
58      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwn
59#else
60      USE oce_trc            , zun => un,  &  ! When no bbl, zun == un
61                               zvn => vn,  &  !              zvn == vn
62                               zwn => wn      !              zwn == wn
63#endif
64
65      !! * Arguments
66      INTEGER, INTENT( in ) ::   kt         ! ocean time-step
67
68      !! * Local declarations
69      INTEGER ::   ji, jj, jk,jn            ! dummy loop indices
70      REAL(wp), DIMENSION (jpi,jpj,jpk) ::   &
71         zt1, zt2, ztp1, ztp2
72
73      REAL(wp) ::   zu, zv, zw, zeu, zev, zew, zbtr, ztra
74      REAL(wp) ::   z0u, z0v, z0w
75      REAL(wp) ::   zzt1, zzt2, zalpha, z2dtt
76#if defined key_trc_diatrd
77      REAL(wp) ::   ztai, ztaj
78      REAL(wp) ::   zfui, zfvj
79#endif
80      CHARACTER (len=22) :: charout
81      !!----------------------------------------------------------------------
82
83
84      IF( kt == nittrc000 .AND. lwp ) THEN
85         WRITE(numout,*)
86         WRITE(numout,*) 'trc_adv : MUSCL advection scheme'
87         WRITE(numout,*) '~~~~~~~'
88      ENDIF
89
90 
91
92#if defined key_trcbbl_adv
93      ! Advective bottom boundary layer
94      ! -------------------------------
95      zun(:,:,:) = un (:,:,:) - u_trc_bbl(:,:,:)
96      zvn(:,:,:) = vn (:,:,:) - v_trc_bbl(:,:,:)
97      zwn(:,:,:) = wn (:,:,:) + w_trc_bbl(:,:,:)
98#endif
99
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
122         CALL lbc_lnk( zt1, 'U', -1. )   
123         CALL lbc_lnk( zt2, 'V', -1. ) 
124
125
126         ! Slopes
127         ! interior values
128         DO jk = 1, jpkm1
129            DO jj = 2, jpj
130               DO ji = fs_2, jpi   ! vector opt.
131                  ztp1(ji,jj,jk) =                    ( zt1(ji,jj,jk) + zt1(ji-1,jj  ,jk) )   &
132                     &           * ( 0.25 + SIGN( 0.25, zt1(ji,jj,jk) * zt1(ji-1,jj  ,jk) ) )
133                  ztp2(ji,jj,jk) =                    ( zt2(ji,jj,jk) + zt2(ji  ,jj-1,jk) )   &
134                     &           * ( 0.25 + SIGN( 0.25, zt2(ji,jj,jk) * zt2(ji  ,jj-1,jk) ) )
135               END DO
136            END DO
137         END DO
138         ! bottom values
139         ztp1(:,:,jpk) = 0.e0
140         ztp2(:,:,jpk) = 0.e0
141
142         ! Slopes limitation
143         DO jk = 1, jpkm1
144            DO jj = 2, jpj
145               DO ji = fs_2, jpi   ! vector opt.
146                  ztp1(ji,jj,jk) = SIGN( 1., ztp1(ji,jj,jk) )   &
147                     &           * MIN(    ABS( ztp1(ji  ,jj,jk) ),   &
148                     &                  2.*ABS( zt1 (ji-1,jj,jk) ),   &
149                     &                  2.*ABS( zt1 (ji  ,jj,jk) ) )
150
151                  ztp2(ji,jj,jk) = SIGN( 1., ztp2(ji,jj,jk) )   &
152                     &           * MIN(    ABS( ztp2(ji,jj  ,jk) ),   &
153                     &                  2.*ABS( zt2 (ji,jj-1,jk) ),   &
154                     &                  2.*ABS( zt2 (ji,jj  ,jk) ) )
155
156               END DO
157            END DO
158         END DO
159
160         ! Advection terms
161         ! interior values
162         DO jk = 1, jpkm1
163            DO jj = 2, jpjm1     
164               DO ji = fs_2, fs_jpim1   ! vector opt.
165                  ! volume fluxes
166#if defined key_s_coord || defined key_partial_steps
167                  zeu = e2u(ji,jj) * fse3u(ji,jj,jk) * zun(ji,jj,jk)
168                  zev = e1v(ji,jj) * fse3v(ji,jj,jk) * zvn(ji,jj,jk)
169#else
170                  zeu = e2u(ji,jj) * zun(ji,jj,jk)
171                  zev = e1v(ji,jj) * zvn(ji,jj,jk)
172#endif
173                  ! MUSCL fluxes
174                  z2dtt = rdttra(jk) * FLOAT(ndttrc)
175                  z0u = SIGN( 0.5, zun(ji,jj,jk) )           
176                  zalpha = 0.5 - z0u
177                  zu  = z0u - 0.5 * zun(ji,jj,jk) * z2dtt / e1u(ji,jj)
178                  zzt1 = trb(ji+1,jj,jk,jn) + zu*ztp1(ji+1,jj,jk)
179                  zzt2 = trb(ji  ,jj,jk,jn) + zu*ztp1(ji  ,jj,jk)
180                  zt1(ji,jj,jk) = zeu * ( zalpha * zzt1 + (1.-zalpha) * zzt2 )
181                  z0v = SIGN( 0.5, zvn(ji,jj,jk) )           
182                  zalpha = 0.5 - z0v
183                  zv  = z0v - 0.5 * zvn(ji,jj,jk) * z2dtt / e2v(ji,jj)
184                  zzt1 = trb(ji,jj+1,jk,jn) + zv*ztp2(ji,jj+1,jk)
185                  zzt2 = trb(ji,jj  ,jk,jn) + zv*ztp2(ji,jj  ,jk)
186                  zt2(ji,jj,jk) = zev * ( zalpha * zzt1 + (1.-zalpha) * zzt2 )
187               END DO
188            END DO
189         END DO
190
191         ! lateral boundary conditions on zt1, zt2 (changed sign)
192         CALL lbc_lnk( zt1, 'U', -1. ) 
193         CALL lbc_lnk( zt2, 'V', -1. ) 
194
195         ! Compute and add the horizontal advective trend
196
197         DO jk = 1, jpkm1
198            DO jj = 2, jpjm1     
199               DO ji = fs_2, fs_jpim1   ! vector opt.
200#if defined key_s_coord || defined key_partial_steps
201                  zbtr = 1. / ( e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk) )
202#else
203                  zbtr = 1. / ( e1t(ji,jj)*e2t(ji,jj) )
204#endif
205                  ! horizontal advective trends
206                  ztra = - zbtr * ( zt1(ji,jj,jk) - zt1(ji-1,jj  ,jk  )   &
207                     &            + zt2(ji,jj,jk) - zt2(ji  ,jj-1,jk  ) )
208                  ! add it to the general tracer trends
209                  tra(ji,jj,jk,jn) = tra(ji,jj,jk,jn) + ztra
210#if defined key_trc_diatrd
211                  ! recompute the trends in i- and j-direction as Uh gradh(T)
212#   if defined key_s_coord || defined key_partial_steps
213                  zfui =  e2u(ji  ,jj) * fse3u(ji,  jj,jk) * un(ji,  jj,jk)   &
214                     & -  e2u(ji-1,jj) * fse3u(ji-1,jj,jk) * un(ji-1,jj,jk)
215                  zfvj =  e1v(ji,jj  ) * fse3v(ji,jj  ,jk) * vn(ji,jj  ,jk)   &
216                     & -  e1v(ji,jj-1) * fse3v(ji,jj-1,jk) * vn(ji,jj-1,jk)
217#   else
218                  zfui = e2u(ji  ,jj) * un(ji,  jj,jk)   &
219                     & - e2u(ji-1,jj) * un(ji-1,jj,jk)
220                  zfvj = e1v(ji,jj  ) * vn(ji,jj  ,jk)   &
221                     & - e1v(ji,jj-1) * vn(ji,jj-1,jk)
222#   endif
223                  ztai =-zbtr * (  zt1(ji,jj,jk) - zt1(ji-1,jj  ,jk) - trn(ji,jj,jk,jn) * zfui  )
224                  ztaj =-zbtr * (  zt2(ji,jj,jk) - zt2(ji  ,jj-1,jk) - trn(ji,jj,jk,jn) * zfvj  )
225                  ! save i- and j- advective trends computed as Uh gradh(T)
226                  trtrd(ji,jj,jk,jn,1) = ztai
227                  trtrd(ji,jj,jk,jn,2) = ztaj
228#endif
229               END DO
230            END DO
231         END DO
232      ENDDO
233
234      IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
235         WRITE(charout, FMT="('muscl - had')")
236         CALL prt_ctl_trc_info(charout)
237         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm,clinfo2='trd')
238      ENDIF
239
240         ! II. Vertical advective fluxes
241         ! -----------------------------
242
243      DO jn = 1, jptra
244         ! First guess of the slope
245         ! interior values
246         DO jk = 2, jpkm1
247            zt1(:,:,jk) = tmask(:,:,jk) * ( trb(:,:,jk-1,jn) - trb(:,:,jk,jn) )
248         END DO
249         ! surface and bottom boundary conditions
250         zt1 (:,:, 1 ) = 0.e0 
251         zt1 (:,:,jpk) = 0.e0
252         ! Slopes
253         DO jk = 2, jpkm1
254            DO jj = 1, jpj
255               DO ji = 1, jpi
256                  ztp1(ji,jj,jk) =                    ( zt1(ji,jj,jk) + zt1(ji,jj,jk+1) )   &
257                     &           * ( 0.25 + SIGN( 0.25, zt1(ji,jj,jk) * zt1(ji,jj,jk+1) ) )
258               END DO
259            END DO
260         END DO
261
262         ! Slopes limitation
263         ! interior values
264         DO jk = 2, jpkm1
265            DO jj = 1, jpj
266               DO ji = 1, jpi
267                  ztp1(ji,jj,jk) = SIGN( 1., ztp1(ji,jj,jk) )   &
268                     &           * MIN(    ABS( ztp1(ji,jj,jk  ) ),   &
269                     &                  2.*ABS( zt1 (ji,jj,jk+1) ),   &
270                     &                  2.*ABS( zt1 (ji,jj,jk  ) ) )
271               END DO
272            END DO
273         END DO
274         ! surface values
275         ztp1(:,:,1) = 0. 
276         ! vertical advective flux
277         ! interior values
278         DO jk = 1, jpkm1
279            DO jj = 2, jpjm1     
280               DO ji = fs_2, fs_jpim1   ! vector opt.
281                  z2dtt = rdttra(jk) * FLOAT(ndttrc)
282                  zew = zwn(ji,jj,jk+1)
283                  z0w = SIGN( 0.5, zwn(ji,jj,jk+1) )
284                  zalpha = 0.5 + z0w
285                  zw  = z0w - 0.5 * zwn(ji,jj,jk+1)*z2dtt / fse3w(ji,jj,jk+1)
286                  zzt1 = trb(ji,jj,jk+1,jn) + zw*ztp1(ji,jj,jk+1)
287                  zzt2 = trb(ji,jj,jk  ,jn) + zw*ztp1(ji,jj,jk  )
288                  zt1(ji,jj,jk+1) = zew * ( zalpha * zzt1 + (1.-zalpha)*zzt2 )
289               END DO
290            END DO
291         END DO
292         ! surface values
293         IF( lk_dynspg_fsc .OR. lk_dynspg_fsc_tsk ) THEN        ! free surface-constant volume
294            zt1(:,:, 1 ) = zwn(:,:,1) * trb(:,:,1,jn)
295         ELSE                                                   ! rigid lid : flux set to zero
296            zt1(:,:, 1 ) = 0.e0
297         ENDIF
298
299         ! bottom values
300         zt1(:,:,jpk) = 0.e0
301
302         ! Compute & add the vertical advective trend
303
304         DO jk = 1, jpkm1
305            DO jj = 2, jpjm1     
306               DO ji = fs_2, fs_jpim1   ! vector opt.
307                  zbtr = 1. / fse3t(ji,jj,jk)
308                  ! horizontal advective trends
309                  ztra = - zbtr * ( zt1(ji,jj,jk) - zt1(ji,jj,jk+1) )
310                  ! add it to the general tracer trends
311                  tra(ji,jj,jk,jn) =  tra(ji,jj,jk,jn) + ztra
312#if defined key_trc_diatrd
313                  ! save the vertical advective trends computed as w gradz(T)
314                  trtrd(ji,jj,jk,jn,3) = ztra - trn(ji,jj,jk,jn) * hdivn(ji,jj,jk)
315#endif
316               END DO
317            END DO
318         END DO
319
320      END DO
321
322      IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
323         WRITE(charout, FMT="('muscl - zad')")
324         CALL prt_ctl_trc_info(charout)
325         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm,clinfo2='trd')
326      ENDIF
327
328END SUBROUTINE trc_adv_muscl
329
330#else
331   !!----------------------------------------------------------------------
332   !!   Default option                                         Empty module
333   !!----------------------------------------------------------------------
334CONTAINS
335   SUBROUTINE trc_adv_muscl( kt ) 
336      INTEGER, INTENT(in) :: kt
337      WRITE(*,*) 'trc_adv_muscl: You should not have seen this print! error?', kt
338   END SUBROUTINE trc_adv_muscl
339#endif
340
341   !!======================================================================
342END MODULE trcadv_muscl
Note: See TracBrowser for help on using the repository browser.