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 tags/nemo_v3_2/nemo_v3_2/NEMO/TOP_SRC/TRP – NEMO

source: tags/nemo_v3_2/nemo_v3_2/NEMO/TOP_SRC/TRP/trcadv_muscl2.F90 @ 1878

Last change on this file since 1878 was 1878, checked in by flavoni, 14 years ago

initial test for nemogcm

File size: 20.0 KB
Line 
1MODULE trcadv_muscl2
2   !!==============================================================================
3   !!                       ***  MODULE  trcadv_muscl2  ***
4   !! Ocean passive tracers:  horizontal & vertical advective trend
5   !!==============================================================================
6#if defined key_top
7   !!----------------------------------------------------------------------
8   !!   'key_top'                                                TOP models
9   !!----------------------------------------------------------------------
10   !!   tra_adv_muscl2 : update the tracer trend with the horizontal
11   !!                    and vertical advection trends using MUSCL2 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 prtctl_trc          ! Print control for debbuging
19   USE trdmld_trc
20   USE trdmld_trc_oce          ! ocean variables trends
21
22   IMPLICIT NONE
23   PRIVATE
24
25   !! * Accessibility
26   PUBLIC trc_adv_muscl2        ! routine called by trcstp.F90
27
28   !! * Module variable
29   REAL(wp), DIMENSION(jpk) ::   &
30      rdttrc                     ! vertical profile of tracer time-step
31
32   !! * Substitutions
33#  include "top_substitute.h90"
34   !!----------------------------------------------------------------------
35   !!   TOP 1.0 , LOCEAN-IPSL (2005)
36   !! $Id: trcadv_muscl2.F90 1528 2009-07-23 14:38:47Z rblod $
37   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
38   !!----------------------------------------------------------------------
39
40CONTAINS
41
42   SUBROUTINE trc_adv_muscl2( kt )
43      !!----------------------------------------------------------------------
44      !!                   ***  ROUTINE trc_adv_muscl2  ***
45      !!
46      !! ** Purpose :   Compute the now trend due to total advection of passi-
47      !!      ve tracer using a MUSCL scheme (Monotone Upstream-
48      !!      Centered Scheme for Conservation Laws) and add it to the general
49      !!      tracer trend.
50      !!
51      !! ** Method : MUSCL scheme plus centered scheme at ocean boundaries
52      !!
53      !! ** Action : - update tra with the now advective tracer trends
54      !!             - save trends ('key_trdmld_trc')
55      !!
56      !! References :               
57      !!      Estubier, A., and M. Levy, Notes Techn. Pole de Modelisation
58      !! IPSL, Sept. 2000 (http://www.lodyc.jussieu.fr/opa)
59      !!
60      !! History :
61      !!        !  06-00  (A.Estublier)  for passive tracers
62      !!   9.0  !  03-04  (C. Ethe, G. Madec)  F90: Free form and module
63      !!----------------------------------------------------------------------
64      !! * modules used
65#if defined key_trcbbl_adv
66      USE oce_trc            , zun => ua,  &  ! use ua as workspace
67         &                     zvn => va      ! use va as workspace
68      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwn
69#else
70      USE oce_trc            , zun => un,  &  ! When no bbl, zun == un
71                               zvn => vn,  &  !              zvn == vn
72                               zwn => wn      !              zwn == wn
73#endif
74      !! * Arguments
75      INTEGER, INTENT( in ) ::   kt         ! ocean time-step
76
77      !! * Local declarations
78      INTEGER ::   ji, jj, jk,jn               ! dummy loop indices
79      REAL(wp), DIMENSION (jpi,jpj,jpk) ::   &
80         zt1, zt2, ztp1, ztp2
81      REAL(wp) ::   zu, zv, zw, zeu, zev, zew, zbtr, ztra
82      REAL(wp) ::   z0u, z0v, z0w
83      REAL(wp) ::   zzt1, zzt2, zalpha
84      REAL(wp) ::   zfui, zfvj
85      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrtrd
86#if defined key_trc_diatrd
87      REAL(wp) ::   ztai, ztaj
88#endif
89      CHARACTER (len=22) :: charout
90      !!----------------------------------------------------------------------
91
92      IF( kt == nittrc000 .AND. lwp ) THEN
93         WRITE(numout,*)
94         WRITE(numout,*) 'trc_adv_muscl2 : MUSCL2 advection scheme'
95         WRITE(numout,*) '~~~~~~~~~~~~~~~'
96         rdttrc(:) =  rdttra(:) * FLOAT(ndttrc)
97      ENDIF
98     
99      IF( l_trdtrc ) ALLOCATE( ztrtrd(jpi,jpj,jpk) )
100
101#if defined key_trcbbl_adv       
102      ! Advective bottom boundary layer
103      ! -------------------------------
104      zun(:,:,:) = un (:,:,:) - u_trc_bbl(:,:,:)
105      zvn(:,:,:) = vn (:,:,:) - v_trc_bbl(:,:,:)
106      zwn(:,:,:) = wn (:,:,:) + w_trc_bbl( :,:,:)
107#endif
108
109
110      DO jn = 1, jptra
111
112         ! I. Horizontal advective fluxes
113         ! ------------------------------
114
115         ! first guess of the slopes
116         ! interior values
117         DO jk = 1, jpkm1
118            DO jj = 1, jpjm1     
119               DO ji = 1, fs_jpim1   ! vector opt.
120                  zt1(ji,jj,jk) = umask(ji,jj,jk) * ( trb(ji+1,jj,jk,jn) - trb(ji,jj,jk,jn) )
121                  zt2(ji,jj,jk) = vmask(ji,jj,jk) * ( trb(ji,jj+1,jk,jn) - trb(ji,jj,jk,jn) )
122               END DO
123            END DO
124         END DO
125         ! bottom values
126         zt1(:,:,jpk) = 0.e0
127         zt2(:,:,jpk) = 0.e0
128
129         ! lateral boundary conditions on zt1, zt2  (changed sign)
130         CALL lbc_lnk( zt1, 'U', -1. )
131         CALL lbc_lnk( zt2, 'V', -1. )
132
133         ! Slopes
134         ! interior values
135         DO jk = 1, jpkm1
136            DO jj = 2, jpj
137               DO ji = fs_2, jpi   ! vector opt.
138                  ztp1(ji,jj,jk) =                    ( zt1(ji,jj,jk) + zt1(ji-1,jj  ,jk) )   &
139                     &           * ( 0.25 + SIGN( 0.25, zt1(ji,jj,jk) * zt1(ji-1,jj  ,jk) ) )
140                  ztp2(ji,jj,jk) =                    ( zt2(ji,jj,jk) + zt2(ji  ,jj-1,jk) )   &
141                     &           * ( 0.25 + SIGN( 0.25, zt2(ji,jj,jk) * zt2(ji  ,jj-1,jk) ) )
142               END DO
143            END DO
144         END DO
145         ! bottom values
146         ztp1(:,:,jpk) = 0.e0 
147         ztp2(:,:,jpk) = 0.e0
148
149         ! Slopes limitation
150         DO jk = 1, jpkm1
151            DO jj = 2, jpj
152               DO ji = fs_2, fs_jpim1   ! vector opt.
153                  ztp1(ji,jj,jk) = SIGN( 1., ztp1(ji,jj,jk) )   &
154                     &           * MIN(    ABS( ztp1(ji  ,jj,jk) ),   &
155                     &                  2.*ABS( zt1 (ji-1,jj,jk) ),   &
156                     &                  2.*ABS( zt1 (ji  ,jj,jk) ) )
157
158                  ztp2(ji,jj,jk) = SIGN( 1., ztp2(ji,jj,jk) )   &
159                     &           * MIN(    ABS( ztp2(ji,jj  ,jk) ),   &
160                     &                  2.*ABS( zt2 (ji,jj-1,jk) ),   &
161                     &                  2.*ABS( zt2 (ji,jj  ,jk) ) )
162               END DO
163            END DO
164         END DO
165
166         ! Advection terms
167         ! interior values
168         DO jk = 1, jpkm1
169            DO jj = 2, jpjm1     
170               DO ji = fs_2, fs_jpim1   ! vector opt.
171                  ! volume fluxes
172#if ! defined key_zco
173                  zeu = e2u(ji,jj) * fse3u(ji,jj,jk) * zun(ji,jj,jk)
174                  zev = e1v(ji,jj) * fse3v(ji,jj,jk) * zvn(ji,jj,jk)
175#else
176                  zeu = e2u(ji,jj) * zun(ji,jj,jk)
177                  zev = e1v(ji,jj) * zvn(ji,jj,jk)
178#endif
179                  ! MUSCL fluxes
180                  z0u = SIGN( 0.5, zun(ji,jj,jk) )           
181                  zalpha = 0.5 - z0u
182                  zu  = z0u - 0.5 * zun(ji,jj,jk) * rdttrc(jk) / e1u(ji,jj)
183                  zzt1 = trb(ji+1,jj,jk,jn) + zu*ztp1(ji+1,jj,jk)
184                  zzt2 = trb(ji  ,jj,jk,jn) + zu*ztp1(ji  ,jj,jk)
185                  zt1(ji,jj,jk) = zeu * ( zalpha * zzt1 + (1.-zalpha) * zzt2 )
186
187                  z0v = SIGN( 0.5, zvn(ji,jj,jk) )           
188                  zalpha = 0.5 - z0v
189                  zv  = z0v - 0.5 * zvn(ji,jj,jk) * rdttrc(jk) / e2v(ji,jj)
190                  zzt1 = trb(ji,jj+1,jk,jn) + zv*ztp2(ji,jj+1,jk)
191                  zzt2 = trb(ji,jj  ,jk,jn) + zv*ztp2(ji,jj  ,jk)
192                  zt2(ji,jj,jk) = zev * ( zalpha * zzt1 + (1.-zalpha) * zzt2 )
193
194               END DO
195            END DO
196         END DO
197
198
199         DO jk = 1, jpkm1
200            DO jj = 2, jpjm1
201               DO ji = fs_2, fs_jpim1   ! vector opt.
202#if ! defined key_zco
203                  zev = e1v(ji,jj) * fse3v(ji,jj,jk)
204                  IF( umask(ji,jj,jk) == 0. ) THEN
205                     IF( zun(ji+1,jj,jk) > 0. .AND. ji /= jpi ) THEN
206                        zt1(ji+1,jj,jk) = e2u(ji+1,jj)* fse3u(ji+1,jj,jk)   &
207                           &            * zun(ji+1,jj,jk) * ( trb(ji+1,jj,jk,jn) + trb(ji+2,jj,jk,jn) ) * 0.5
208                     ENDIF
209                     IF( zun(ji-1,jj,jk) < 0. ) THEN
210                        zt1(ji-1,jj,jk) = e2u(ji-1,jj)* fse3u(ji-1,jj,jk)   &
211                           &            * zun(ji-1,jj,jk) * ( trb(ji-1,jj,jk,jn) + trb(ji  ,jj,jk,jn) ) * 0.5
212                     ENDIF
213                  ENDIF
214                  IF( vmask(ji,jj,jk) == 0. ) THEN
215                     IF( zvn(ji,jj+1,jk) > 0. .AND. jj /= jpj ) THEN
216                        zt2(ji,jj+1,jk) = e1v(ji,jj+1) * fse3v(ji,jj+1,jk)   &
217                           &            * zvn(ji,jj+1,jk) * ( trb(ji,jj+1,jk,jn) + trb(ji,jj+2,jk,jn) ) * 0.5
218                     ENDIF
219                     IF( zvn(ji,jj-1,jk) < 0. ) THEN
220                        zt2(ji,jj-1,jk) = e1v(ji,jj-1)* fse3v(ji,jj-1,jk)   &
221                           &            * zvn(ji,jj-1,jk) * ( trb(ji,jj-1,jk,jn) + trb(ji  ,jj,jk,jn) ) * 0.5
222                     ENDIF
223                  ENDIF
224
225#else
226                  IF( umask(ji,jj,jk) == 0. ) THEN
227                     IF( zun(ji+1,jj,jk) > 0. .AND. ji /= jpi ) THEN
228                        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
229                     ENDIF
230                     IF( zun(ji-1,jj,jk) < 0. ) THEN
231                        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
232                     ENDIF
233                  ENDIF
234                  IF( vmask(ji,jj,jk) == 0. ) THEN
235                     IF( zvn(ji,jj+1,jk) > 0. .AND. jj /= jpj ) THEN
236                        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
237                     ENDIF
238                     IF( zvn(ji,jj-1,jk) < 0. ) THEN
239                        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
240                     ENDIF
241                  ENDIF
242#endif
243               END DO
244            END DO
245         END DO
246
247         ! lateral boundary conditions on zt1, zt2   (changed sign)
248         CALL lbc_lnk( zt1, 'U', -1. )
249         CALL lbc_lnk( zt2, 'V', -1. )
250
251         ! Compute and add the horizontal advective trend
252
253         DO jk = 1, jpkm1
254            DO jj = 2, jpjm1     
255               DO ji = fs_2, fs_jpim1   ! vector opt.
256#if ! defined key_zco
257                  zbtr = 1. / ( e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk) )
258#else
259                  zbtr = 1. / ( e1t(ji,jj)*e2t(ji,jj) )
260#endif
261                  ! horizontal advective trends
262                  ztra = - zbtr * ( zt1(ji,jj,jk) - zt1(ji-1,jj  ,jk  )   &
263                     &            + zt2(ji,jj,jk) - zt2(ji  ,jj-1,jk  ) )
264                  ! add it to the general tracer trends
265                  tra(ji,jj,jk,jn) = tra(ji,jj,jk,jn) + ztra
266#if defined key_trc_diatrd
267                  ! recompute the trends in i- and j-direction as Uh gradh(T)
268#   if ! defined key_zco
269                  zfui =  e2u(ji  ,jj) * fse3u(ji,  jj,jk) * un(ji,  jj,jk)   &
270                     & -  e2u(ji-1,jj) * fse3u(ji-1,jj,jk) * un(ji-1,jj,jk)
271                  zfvj =  e1v(ji,jj  ) * fse3v(ji,jj  ,jk) * vn(ji,jj  ,jk)   &
272                     & -  e1v(ji,jj-1) * fse3v(ji,jj-1,jk) * vn(ji,jj-1,jk)
273#   else
274                  zfui = e2u(ji  ,jj) * un(ji,  jj,jk)   &
275                     & - e2u(ji-1,jj) * un(ji-1,jj,jk)
276                  zfvj = e1v(ji,jj  ) * vn(ji,jj  ,jk)   &
277                     & - e1v(ji,jj-1) * vn(ji,jj-1,jk)
278#   endif
279                  ztai =-zbtr * (  zt1(ji,jj,jk) - zt1(ji-1,jj  ,jk) - trn(ji,jj,jk,jn) * zfui  )
280                  ztaj =-zbtr * (  zt2(ji,jj,jk) - zt2(ji  ,jj-1,jk) - trn(ji,jj,jk,jn) * zfvj  )
281                  ! save i- and j- advective trends computed as Uh gradh(T)
282                  IF (luttrd(jn)) trtrd(ji,jj,jk,ikeep(jn),1) = ztai
283                  IF (luttrd(jn)) trtrd(ji,jj,jk,ikeep(jn),2) = ztaj
284
285#endif
286
287               END DO
288            END DO
289         END DO
290
291         ! 3. Save the horizontal advective trends for diagnostics
292         ! -------------------------------------------------------
293
294         TRDTRC_XY : IF( l_trdtrc ) THEN
295
296            ! 3.1) Passive tracer ZONAL advection trends
297            DO jk = 1, jpkm1
298               DO jj = 2, jpjm1
299                  DO ji = fs_2, fs_jpim1   ! vector opt. 
300#if ! defined key_zco
301                     zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
302                     zfui = e2u(ji  ,jj) * fse3u(ji,  jj,jk) * un(ji,  jj,jk)   &
303                        & - e2u(ji-1,jj) * fse3u(ji-1,jj,jk) * un(ji-1,jj,jk)
304#else
305                     zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) )
306                     zfui = e2u(ji  ,jj) * un(ji,  jj,jk)   &
307                        & - e2u(ji-1,jj) * un(ji-1,jj,jk)
308#endif
309                     ! recompute the trends in i- direction as Uh gradh(T)
310                     ztrtrd(ji,jj,jk) = - zbtr*( zt1(ji,jj,jk) - zt1(ji-1,jj,jk) - trn(ji,jj,jk,jn)*zfui )
311                  END DO
312               END DO
313            END DO
314
315            IF (luttrd(jn)) CALL trd_mod_trc( ztrtrd, jn, jptrc_trd_xad, kt )
316
317            ! 3.2)  Passive tracer MERIDIONAL advection trends
318            DO jk = 1, jpkm1
319               DO jj = 2, jpjm1
320                  DO ji = fs_2, fs_jpim1   ! vector opt.
321                     ! recompute the trends in i- and j-direction as Uh gradh(T)
322#if ! defined key_zco
323                     zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
324                     zfvj = e1v(ji,jj  ) * fse3v(ji,jj  ,jk) * vn(ji,jj  ,jk)   &
325                        & - e1v(ji,jj-1) * fse3v(ji,jj-1,jk) * vn(ji,jj-1,jk)
326#else
327                     zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) )
328                     zfvj = e1v(ji,jj  ) * vn(ji,jj  ,jk)   &
329                        & - e1v(ji,jj-1) * vn(ji,jj-1,jk)
330#endif
331                     ztrtrd(ji,jj,jk) = - zbtr*( zt2(ji,jj,jk) - zt2(ji,jj-1,jk) - trn(ji,jj,jk,jn)*zfvj )
332                  END DO
333               END DO
334            END DO
335
336            IF (luttrd(jn)) CALL trd_mod_trc( ztrtrd, jn, jptrc_trd_yad, kt )
337
338         ENDIF TRDTRC_XY
339
340         !                           !=============
341      END DO                         ! tracer loop
342      !                              !=============
343
344      IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
345         WRITE(charout, FMT="('muscl2 - had')")
346         CALL prt_ctl_trc_info(charout)
347         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm,clinfo2='trd')
348      ENDIF
349
350      ! II. Vertical advective fluxes
351      ! -----------------------------
352
353      DO jn = 1, jptra
354
355         ! First guess of the slope
356         ! interior values
357         DO jk = 2, jpkm1
358            zt1(:,:,jk) = tmask(:,:,jk) * ( trb(:,:,jk-1,jn) - trb(:,:,jk,jn) )
359         END DO
360         ! surface and bottom boundary conditions
361         zt1 (:,:, 1 ) = 0.e0 
362         zt1 (:,:,jpk) = 0.e0
363
364         ! Slopes
365         DO jk = 2, jpkm1
366            DO jj = 1, jpj
367               DO ji = 1, jpi
368                  ztp1(ji,jj,jk) =                    ( zt1(ji,jj,jk) + zt1(ji,jj,jk+1) )   &
369                     &           * ( 0.25 + SIGN( 0.25, zt1(ji,jj,jk) * zt1(ji,jj,jk+1) ) )
370               END DO
371            END DO
372         END DO
373
374         ! Slopes limitation
375         ! interior values
376         DO jk = 2, jpkm1
377            DO jj = 1, jpj
378               DO ji = 1, jpi
379                  ztp1(ji,jj,jk) = SIGN( 1., ztp1(ji,jj,jk) )   &
380                     &           * MIN(    ABS( ztp1(ji,jj,jk  ) ),   &
381                     &                  2.*ABS( zt1 (ji,jj,jk+1) ),   &
382                     &                  2.*ABS( zt1 (ji,jj,jk  ) ) )
383               END DO
384            END DO
385         END DO
386
387         ! surface values
388         ztp1(:,:,1) = 0. 
389
390         ! vertical advective flux
391         ! interior values
392         DO jk = 1, jpkm1
393            DO jj = 2, jpjm1     
394               DO ji = fs_2, fs_jpim1   ! vector opt.
395                  zew = zwn(ji,jj,jk+1)
396                  z0w = SIGN( 0.5, zwn(ji,jj,jk+1) )
397                  zalpha = 0.5 + z0w
398                  zw  = z0w - 0.5 * zwn(ji,jj,jk+1)* rdttrc(jk)/ fse3w(ji,jj,jk+1)
399                  zzt1 = trb(ji,jj,jk+1,jn) + zw*ztp1(ji,jj,jk+1)
400                  zzt2 = trb(ji,jj,jk  ,jn) + zw*ztp1(ji,jj,jk  )
401                  zt1(ji,jj,jk+1) = zew * ( zalpha * zzt1 + (1.-zalpha)*zzt2 )
402               END DO
403            END DO
404         END DO
405         DO jk = 2, jpkm1
406            DO jj = 2, jpjm1
407               DO ji = fs_2, fs_jpim1   ! vector opt.
408                  IF( tmask(ji,jj,jk+1) == 0. ) THEN
409                     IF( zwn(ji,jj,jk) > 0. ) THEN
410                        zt1(ji,jj,jk) = zwn(ji,jj,jk) * ( trb(ji,jj,jk-1,jn) + trb(ji,jj,jk,jn) ) * 0.5
411                     ENDIF
412                  ENDIF
413               END DO
414            END DO
415         END DO
416
417         ! surface values
418         IF( lk_vvl ) THEN
419            ! variable volume: flux set to zero
420            zt1(:,:, 1 ) = 0.e0
421         ELSE
422            ! free surface-constant volume
423            zt1(:,:, 1 ) = zwn(:,:,1) * trb(:,:,1,jn)
424         ENDIF
425
426         ! bottom values
427         zt1(:,:,jpk) = 0.e0
428
429         ! Compute & add the vertical advective trend
430
431         DO jk = 1, jpkm1
432            DO jj = 2, jpjm1     
433               DO ji = fs_2, fs_jpim1   ! vector opt.
434                  zbtr = 1. / fse3t(ji,jj,jk)
435                  ! horizontal advective trends
436                  ztra = - zbtr * ( zt1(ji,jj,jk) - zt1(ji,jj,jk+1) )
437                  ! add it to the general tracer trends
438                  tra(ji,jj,jk,jn) =  tra(ji,jj,jk,jn) + ztra
439#if defined key_trc_diatrd
440                  ! save the vertical advective trends computed as w gradz(T)
441                  IF (luttrd(jn)) trtrd(ji,jj,jk,ikeep(jn),3) = ztra - trn(ji,jj,jk,jn) * hdivn(ji,jj,jk)
442#endif
443
444               END DO
445            END DO
446         END DO
447
448
449         ! 3. Save the vertical advective trends for diagnostic
450         ! ----------------------------------------------------
451
452         TRDTRC_Z : IF( l_trdtrc )THEN
453
454            ! Compute T/S vertical advection trends
455            DO jk = 1, jpkm1
456               DO jj = 2, jpjm1
457                  DO ji = fs_2, fs_jpim1   ! vector opt.
458                     zbtr = 1. / fse3t(ji,jj,jk)
459                     ! horizontal advective trends
460                     ztra = - zbtr * ( zt1(ji,jj,jk) - zt1(ji,jj,jk+1) )
461                     ! save the vertical advective trends computed as w gradz(T)
462                     ztrtrd(ji,jj,jk) = ztra - trn(ji,jj,jk,jn) * hdivn(ji,jj,jk)
463                  END DO
464               END DO
465            END DO
466
467            IF (luttrd(jn)) CALL trd_mod_trc(ztrtrd, jn, jptrc_trd_zad, kt)
468
469         END IF TRDTRC_Z
470
471         !                            !=============
472      END DO                          ! tracer loop
473      !                               !=============
474
475
476      IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
477         WRITE(charout, FMT="('muscl2 - zad')")
478         CALL prt_ctl_trc_info(charout)
479         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm,clinfo2='trd')
480      ENDIF
481
482      IF( l_trdtrc ) DEALLOCATE( ztrtrd )
483
484   END SUBROUTINE trc_adv_muscl2
485
486#else
487
488   !!----------------------------------------------------------------------
489   !!   Default option                                         Empty module
490   !!----------------------------------------------------------------------
491CONTAINS
492   SUBROUTINE trc_adv_muscl2( kt ) 
493      INTEGER, INTENT(in) :: kt
494      WRITE(*,*) 'trc_adv_muscl2: You should not have seen this print! error?', kt
495   END SUBROUTINE trc_adv_muscl2
496#endif
497
498   !!======================================================================
499END MODULE trcadv_muscl2
Note: See TracBrowser for help on using the repository browser.