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.
traadv_muscl.F90 in branches/dev_001_GM/NEMO/OPA_SRC/TRA – NEMO

source: branches/dev_001_GM/NEMO/OPA_SRC/TRA/traadv_muscl.F90 @ 786

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

dev_001_GM - merge TRC-TRA on OPA only, trabbl & zpshde not done and trdmld not OK - compilation OK
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 12.4 KB
Line 
1MODULE traadv_muscl
2   !!======================================================================
3   !!                       ***  MODULE  traadv_muscl  ***
4   !! Ocean active tracers:  horizontal & vertical advective trend
5   !!======================================================================
6   !! History :       !  06-00  (A.Estublier)  for passive tracers
7   !!                 !  01-08  (E.Durand, G.Madec)  adapted for T & S
8   !!   NEMO     1.0  !  02-06  (G. Madec)  F90: Free form and module
9   !!            2.4  !  08-01  (G. Madec) Merge TRA-TRC
10   !!----------------------------------------------------------------------
11
12   !!----------------------------------------------------------------------
13   !!   tra_adv_muscl : update the tracer trend with the horizontal
14   !!                   and vertical advection trends using MUSCL scheme
15   !!----------------------------------------------------------------------
16   USE dom_oce         ! ocean space and time domain
17   USE trdmod          ! ocean active tracers trends
18   USE trdmod_oce      ! ocean variables trends
19   USE in_out_manager  ! I/O manager
20   USE dynspg_oce      ! choice/control of key cpp for surface pressure gradient
21   USE trabbl          ! tracers: bottom boundary layer
22   USE lib_mpp         ! distribued memory computing
23   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
24   USE diaptr          ! poleward transport diagnostics
25   USE prtctl          ! Print control
26
27   IMPLICIT NONE
28   PRIVATE
29
30   PUBLIC   tra_adv_muscl  ! routine called by step.F90
31
32   !! * Substitutions
33#  include "domzgr_substitute.h90"
34#  include "vectopt_loop_substitute.h90"
35   !!----------------------------------------------------------------------
36   !! NEMO/OPA 2.4 , LOCEAN-IPSL (2008)
37   !! $Id:$
38   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
39   !!----------------------------------------------------------------------
40
41CONTAINS
42
43   SUBROUTINE tra_adv_muscl( kt, cdtype, ktra, pun, pvn, pwn,   &
44      &                                        ptb     , pta )
45      !!----------------------------------------------------------------------
46      !!                    ***  ROUTINE tra_adv_muscl  ***
47      !!
48      !! ** Purpose :   Compute the now trend due to total advection of T and
49      !!      S using a MUSCL scheme (Monotone Upstream-centered Scheme for
50      !!      Conservation Laws) and add it to the general tracer trend.
51      !!
52      !! ** Method  : MUSCL scheme plus centered scheme at ocean boundaries
53      !!
54      !! ** Action  : - update (pta,sa) with the now advective tracer trends
55      !!              - save trends in (ztrdt,ztrds) ('key_trdtra')
56      !!
57      !! References : Estubier, A., and M. Levy, Notes Techn. Pole de Modelisation
58      !!              IPSL, Sept. 2000 (http://www.lodyc.jussieu.fr/opa)
59      !!----------------------------------------------------------------------
60      INTEGER         , INTENT(in   )                         ::   kt              ! ocean time-step index
61      CHARACTER(len=3), INTENT(in   )                         ::   cdtype          ! =TRA or TRC (tracer indicator)
62      INTEGER         , INTENT(in   )                         ::   ktra            ! tracer index
63      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   pun, pvn, pwn   ! 3 ocean velocity components
64      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   ptb             ! before tracer fields
65      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pta             ! tracer trend
66      !!
67      INTEGER  ::   ji, jj, jk   ! dummy loop indices
68      REAL(wp) ::   zu, zv, zw, zeu, zev 
69      REAL(wp) ::   zew, zbtr, z2, zstep 
70      REAL(wp) ::   z0u, z0v, z0w 
71      REAL(wp) ::   zzwx, zzwy, zalpha 
72      REAL(wp) ::   zta
73      REAL(wp), DIMENSION (jpi,jpj,jpk) ::   zwx, zwy, zslpx, zslpy   ! 3D workspace
74      !!----------------------------------------------------------------------
75
76      IF( kt == nit000 .AND. lwp ) THEN
77         WRITE(numout,*)
78         WRITE(numout,*) 'tra_adv : MUSCL advection scheme'
79         WRITE(numout,*) '~~~~~~~'
80      ENDIF
81
82      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;    z2 = 1.
83      ELSE                                        ;    z2 = 2.
84      ENDIF
85
86      ! I. Horizontal advective fluxes
87      ! ------------------------------
88      ! first guess of the slopes
89      ! interior values
90      DO jk = 1, jpkm1
91         DO jj = 1, jpjm1     
92            DO ji = 1, fs_jpim1   ! vector opt.
93               zwx(ji,jj,jk) = umask(ji,jj,jk) * ( ptb(ji+1,jj,jk) - ptb(ji,jj,jk) )
94               zwy(ji,jj,jk) = vmask(ji,jj,jk) * ( ptb(ji,jj+1,jk) - ptb(ji,jj,jk) )
95            END DO
96         END DO
97      END DO
98      ! bottom values
99      zwx(:,:,jpk) = 0.e0    ;    zwy(:,:,jpk) = 0.e0
100
101      ! lateral boundary conditions on zwx, zwy   (changed sign)
102      CALL lbc_lnk( zwx, 'U', -1. )   ;   CALL lbc_lnk( zwy, 'V', -1. )
103
104      ! Slopes
105      ! interior values
106      DO jk = 1, jpkm1
107         DO jj = 2, jpj
108            DO ji = fs_2, jpi   ! vector opt.
109               zslpx(ji,jj,jk) =                    ( zwx(ji,jj,jk) + zwx(ji-1,jj  ,jk) )   &
110                  &           * ( 0.25 + SIGN( 0.25, zwx(ji,jj,jk) * zwx(ji-1,jj  ,jk) ) )
111               zslpy(ji,jj,jk) =                    ( zwy(ji,jj,jk) + zwy(ji  ,jj-1,jk) )   &
112                  &           * ( 0.25 + SIGN( 0.25, zwy(ji,jj,jk) * zwy(ji  ,jj-1,jk) ) )
113            END DO
114         END DO
115      END DO
116      ! bottom values
117      zslpx(:,:,jpk) = 0.e0    ;    zslpy(:,:,jpk) = 0.e0
118
119      ! Slopes limitation
120      DO jk = 1, jpkm1
121         DO jj = 2, jpj
122            DO ji = fs_2, jpi   ! vector opt.
123               zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) )   &
124                  &            * MIN(    ABS( zslpx(ji  ,jj,jk) ),   &
125                  &                   2.*ABS( zwx (ji-1,jj,jk) ),   &
126                  &                   2.*ABS( zwx (ji  ,jj,jk) ) )
127               zslpy(ji,jj,jk) = SIGN( 1., zslpy(ji,jj,jk) )   &
128                  &            * MIN(    ABS( zslpy(ji,jj  ,jk) ),   &
129                  &                   2.*ABS( zwy (ji,jj-1,jk) ),   &
130                  &                   2.*ABS( zwy (ji,jj  ,jk) ) )
131            END DO
132         END DO
133      END DO       
134
135      ! Advection terms
136      ! interior values
137      DO jk = 1, jpkm1
138         zstep  = z2 * rdttra(jk)
139         DO jj = 2, jpjm1     
140            DO ji = fs_2, fs_jpim1   ! vector opt.
141               ! volume fluxes
142#if defined key_zco
143               zeu = e2u(ji,jj)                   * pun(ji,jj,jk)
144               zev = e1v(ji,jj)                   * pvn(ji,jj,jk)
145#else
146               zeu = e2u(ji,jj) * fse3u(ji,jj,jk) * pun(ji,jj,jk)
147               zev = e1v(ji,jj) * fse3v(ji,jj,jk) * pvn(ji,jj,jk)
148#endif
149               ! MUSCL fluxes
150               z0u = SIGN( 0.5, pun(ji,jj,jk) )           
151               zalpha = 0.5 - z0u
152               zu  = z0u - 0.5 * pun(ji,jj,jk) * zstep / e1u(ji,jj)
153               zzwx = ptb(ji+1,jj,jk) + zu * zslpx(ji+1,jj,jk)
154               zzwy = ptb(ji  ,jj,jk) + zu * zslpx(ji  ,jj,jk)
155               zwx(ji,jj,jk) = zeu * ( zalpha * zzwx + (1.-zalpha) * zzwy )
156               !
157               z0v = SIGN( 0.5, pvn(ji,jj,jk) )           
158               zalpha = 0.5 - z0v
159               zv  = z0v - 0.5 * pvn(ji,jj,jk) * zstep / e2v(ji,jj)
160               zzwx = ptb(ji,jj+1,jk) + zv * zslpy(ji,jj+1,jk)
161               zzwy = ptb(ji,jj  ,jk) + zv * zslpy(ji,jj  ,jk)
162               zwy(ji,jj,jk) = zev * ( zalpha * zzwx + (1.-zalpha) * zzwy )
163            END DO
164         END DO
165      END DO
166
167!!gm bug?   there is too many lbc: this have to be checked
168      ! lateral boundary conditions on zwx, zwy   (changed sign)
169      CALL lbc_lnk( zwx, 'U', -1. )   ;   CALL lbc_lnk( zwy, 'V', -1. )
170
171      ! Tracer flux divergence at t-point added to the general trend
172      DO jk = 1, jpkm1
173         DO jj = 2, jpjm1     
174            DO ji = fs_2, fs_jpim1   ! vector opt.
175#if defined key_zco
176               zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) )
177#else
178               zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
179#endif
180               ! horizontal advective trends
181               zta = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
182                  &           + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  ) ) 
183               ! add it to the general tracer trends
184               pta(ji,jj,jk) = pta(ji,jj,jk) + zta
185            END DO
186        END DO
187      END DO       
188
189      IF(ln_ctl)   CALL prt_ctl( tab3d_1=pta, clinfo1=' muscl - had: ', mask1=tmask, clinfo3=cdtype )
190
191      ! Save the horizontal advective trends for diagnostics
192      IF( l_trdtra ) THEN
193         CALL trd_tra_adv( kt, ktra, jpt_trd_xad, cdtype, zwx, pun, ptb )
194         CALL trd_tra_adv( kt, ktra, jpt_trd_yad, cdtype, zwy, pvn, ptb )
195      ENDIF
196
197      IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN
198         IF( lk_zco ) THEN
199            DO jk = 1, jpkm1
200               DO jj = 2, jpjm1
201                  DO ji = fs_2, fs_jpim1   ! vector opt.
202                    zwy(ji,jj,jk) = zwy(ji,jj,jk) * fse3v(ji,jj,jk)
203                  END DO
204               END DO
205            END DO
206         ENDIF
207         IF( ktra == jp_tem)   pht_adv(:) = ptr_vj( zwy(:,:,:) )
208         IF( ktra == jp_sal)   pst_adv(:) = ptr_vj( zwy(:,:,:) )
209      ENDIF
210
211
212      ! II. Vertical advective fluxes
213      ! -----------------------------
214     
215      ! First guess of the slope
216      ! interior values
217      DO jk = 2, jpkm1
218         zwx(:,:,jk) = tmask(:,:,jk) * ( ptb(:,:,jk-1) - ptb(:,:,jk) )
219      END DO
220      ! surface & bottom boundary conditions
221      zwx (:,:, 1 ) = 0.e0    ;    zwx (:,:,jpk) = 0.e0
222
223      ! Slopes
224      DO jk = 2, jpkm1
225         DO jj = 1, jpj
226            DO ji = 1, jpi
227               zslpx(ji,jj,jk) =                    ( zwx(ji,jj,jk) + zwx(ji,jj,jk+1) )   &
228                  &            * ( 0.25 + SIGN( 0.25, zwx(ji,jj,jk) * zwx(ji,jj,jk+1) ) )
229            END DO
230         END DO
231      END DO
232
233      ! Slopes limitation
234      DO jk = 2, jpkm1        ! interior values
235         DO jj = 1, jpj
236            DO ji = 1, jpi
237               zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) )   &
238                  &            * MIN(    ABS( zslpx(ji,jj,jk  ) ),   &
239                  &                   2.*ABS( zwx (ji,jj,jk+1) ),   &
240                  &                   2.*ABS( zwx (ji,jj,jk  ) ) )
241            END DO
242         END DO
243      END DO
244      zslpx(:,:,1) = 0.e0      ! surface values
245
246      ! vertical advective flux
247      DO jk = 1, jpkm1        ! interior values
248         zstep  = z2 * rdttra(jk)
249         DO jj = 2, jpjm1     
250            DO ji = fs_2, fs_jpim1   ! vector opt.
251               zew = pwn(ji,jj,jk+1)
252               z0w = SIGN( 0.5, pwn(ji,jj,jk+1) )
253               zalpha = 0.5 + z0w
254               zw  = z0w - 0.5 * pwn(ji,jj,jk+1) * zstep / fse3w(ji,jj,jk+1)
255               zzwx = ptb(ji,jj,jk+1) + zw * zslpx(ji,jj,jk+1)
256               zzwy = ptb(ji,jj,jk  ) + zw * zslpx(ji,jj,jk  )
257               zwx(ji,jj,jk+1) = zew * ( zalpha * zzwx + (1.-zalpha) * zzwy )
258            END DO
259         END DO
260      END DO
261      !                       ! surface values
262      IF( lk_dynspg_rl .OR. lk_vvl) THEN                ! rigid lid or variable volume: flux set to zero
263         zwx(:,:, 1 ) = 0.e0                                 ! surface
264         zwx(:,:,jpk) = 0.e0                                 ! bottom
265      ELSE                                              ! free surface
266         zwx(:,:, 1 ) = pwn(:,:,1) * ptb(:,:,1)              ! Surface
267         zwx(:,:,jpk) = 0.e0                                 ! bottom
268
269      ENDIF
270
271      ! Compute & add the vertical advective trend
272      DO jk = 1, jpkm1
273         DO jj = 2, jpjm1     
274            DO ji = fs_2, fs_jpim1   ! vector opt.
275               zbtr = 1. / fse3t(ji,jj,jk)
276               ! horizontal advective trends
277               zta = - zbtr * ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) )
278               ! add it to the general tracer trends
279               pta(ji,jj,jk) =  pta(ji,jj,jk) + zta
280            END DO
281         END DO
282      END DO
283
284      ! Save the vertical advective trends for diagnostic
285      ! -------------------------------------------------
286      IF( l_trdtra )   CALL trd_tra_adv( kt, ktra, jpt_trd_zad, cdtype, zwx, pwn, ptb )
287
288      IF(ln_ctl)   CALL prt_ctl( tab3d_1=pta, clinfo1=' muscl - zad: ', mask1=tmask, clinfo3=cdtype )
289      !
290   END SUBROUTINE tra_adv_muscl
291
292   !!======================================================================
293END MODULE traadv_muscl
Note: See TracBrowser for help on using the repository browser.