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_muscl2.F90 in branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_muscl2.F90 @ 11101

Last change on this file since 11101 was 11101, checked in by frrh, 5 years ago

Merge changes from Met Office GMED ticket 450 to reduce unnecessary
text output from NEMO.
This output, which is typically not switchable, is rarely of interest
in normal (non-debugging) runs and simply redunantley consumes extra
file space.
Further, the presence of this text output has been shown to
significantly degrade performance of models which are run during
Met Office HPC RAID (disk) checks.
The new code introduces switches which are configurable via the
changes made in the associated Met Office MOCI ticket 399.
File size: 15.0 KB
Line 
1MODULE traadv_muscl2
2   !!======================================================================
3   !!                       ***  MODULE  traadv_muscl2  ***
4   !! Ocean  tracers:  horizontal & vertical advective trend
5   !!======================================================================
6   !! History :  1.0  !  2002-06  (G. Madec) from traadv_muscl
7   !!            3.2  !  2010-05  (C. Ethe, G. Madec)  merge TRC-TRA + switch from velocity to transport
8   !!----------------------------------------------------------------------
9
10   !!----------------------------------------------------------------------
11   !!   tra_adv_muscl2 : update the tracer trend with the horizontal
12   !!                    and vertical advection trends using MUSCL2 scheme
13   !!----------------------------------------------------------------------
14   USE oce             ! ocean dynamics and active tracers
15   USE trc_oce         ! share passive tracers/Ocean variables
16   USE dom_oce         ! ocean space and time domain
17   USE trd_oce         ! trends: ocean variables
18   USE trdtra          ! trends manager: tracers
19   USE in_out_manager  ! I/O manager
20   USE dynspg_oce      ! choice/control of key cpp for surface pressure gradient
21   USE diaptr          ! poleward transport diagnostics
22   !
23   USE lib_mpp         ! distribued memory computing
24   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
25   USE wrk_nemo        ! Memory Allocation
26   USE timing          ! Timing
27   USE lib_fortran     ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
28
29   IMPLICIT NONE
30   PRIVATE
31
32   PUBLIC   tra_adv_muscl2        ! routine called by step.F90
33
34   !! * Substitutions
35#  include "domzgr_substitute.h90"
36#  include "vectopt_loop_substitute.h90"
37   !!----------------------------------------------------------------------
38   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
39   !! $Id$
40   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
41   !!----------------------------------------------------------------------
42CONTAINS
43
44   SUBROUTINE tra_adv_muscl2( kt, kit000, cdtype, p2dt, pun, pvn, pwn,      &
45      &                                         ptb, ptn, pta, kjpt )
46      !!----------------------------------------------------------------------
47      !!                   ***  ROUTINE tra_adv_muscl2  ***
48      !!
49      !! ** Purpose :   Compute the now trend due to total advection of T and
50      !!      S using a MUSCL scheme (Monotone Upstream-centered Scheme for
51      !!      Conservation Laws) and add it to the general tracer trend.
52      !!
53      !! ** Method  : MUSCL scheme plus centered scheme at ocean boundaries
54      !!
55      !! ** Action  : - update (pta) with the now advective tracer trends
56      !!              - save trends
57      !!
58      !! References : Estubier, A., and M. Levy, Notes Techn. Pole de Modelisation
59      !!              IPSL, Sept. 2000 (http://www.lodyc.jussieu.fr/opa)
60      !!----------------------------------------------------------------------
61      INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index
62      INTEGER                              , INTENT(in   ) ::   kit000          ! first time step index
63      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator)
64      INTEGER                              , INTENT(in   ) ::   kjpt            ! number of tracers
65      REAL(wp), DIMENSION(        jpk     ), INTENT(in   ) ::   p2dt            ! vertical profile of tracer time-step
66      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pun, pvn, pwn   ! 3 ocean velocity components
67      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb, ptn        ! before & now tracer fields
68      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta             ! tracer trend
69      !!
70      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
71      REAL(wp) ::   zu, z0u, zzwx, zw         ! local scalars
72      REAL(wp) ::   zv, z0v, zzwy, z0w        !   -      -
73      REAL(wp) ::   ztra, zbtr, zdt, zalpha   !   -      -
74      REAL(wp), POINTER, DIMENSION(:,:,:) :: zslpx, zslpy , zwx, zwy
75      !!----------------------------------------------------------------------
76      !
77      IF( nn_timing == 1 )  CALL timing_start('tra_adv_muscl2')
78      !
79      CALL wrk_alloc( jpi, jpj, jpk, zslpx, zslpy, zwx, zwy )
80      !
81
82      IF( kt == kit000 )  THEN
83         IF(lwp) WRITE(numout,*)
84         IF(lwp) WRITE(numout,*) 'tra_adv_muscl2 : MUSCL2 advection scheme on ', cdtype
85         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~'
86         IF(lwp .AND. lflush) CALL flush(numout)
87      ENDIF
88      !
89      !                                                          ! ===========
90      DO jn = 1, kjpt                                            ! tracer loop
91         !                                                       ! ===========
92         ! I. Horizontal advective fluxes
93         ! ------------------------------
94         ! first guess of the slopes
95         zwx(:,:,jpk) = 0.e0   ;   zwy(:,:,jpk) = 0.e0        ! bottom values
96         ! interior values
97         DO jk = 1, jpkm1
98            DO jj = 1, jpjm1     
99               DO ji = 1, fs_jpim1   ! vector opt.
100                  zwx(ji,jj,jk) = umask(ji,jj,jk) * ( ptb(ji+1,jj,jk,jn) - ptb(ji,jj,jk,jn) )
101                  zwy(ji,jj,jk) = vmask(ji,jj,jk) * ( ptb(ji,jj+1,jk,jn) - ptb(ji,jj,jk,jn) )
102               END DO
103           END DO
104         END DO
105         !
106         CALL lbc_lnk( zwx, 'U', -1. )                        ! lateral boundary conditions on zwx, zwy   (changed sign)
107         CALL lbc_lnk( zwy, 'V', -1. )
108         !                                             !-- Slopes of tracer
109         zslpx(:,:,jpk) = 0.e0   ;   zslpy(:,:,jpk) = 0.e0    ! bottom values
110         DO jk = 1, jpkm1                                     ! interior values
111            DO jj = 2, jpj
112               DO ji = fs_2, jpi   ! vector opt.
113                  zslpx(ji,jj,jk) =                    ( zwx(ji,jj,jk) + zwx(ji-1,jj  ,jk) )   &
114                     &            * ( 0.25 + SIGN( 0.25, zwx(ji,jj,jk) * zwx(ji-1,jj  ,jk) ) )
115                  zslpy(ji,jj,jk) =                    ( zwy(ji,jj,jk) + zwy(ji  ,jj-1,jk) )   &
116                     &            * ( 0.25 + SIGN( 0.25, zwy(ji,jj,jk) * zwy(ji  ,jj-1,jk) ) )
117               END DO
118            END DO
119         END DO
120         !
121         DO jk = 1, jpkm1                                     ! Slopes limitation
122            DO jj = 2, jpj
123               DO ji = fs_2, jpi   ! vector opt.
124                  zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) ) * 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) ) * MIN(    ABS( zslpy(ji,jj  ,jk) ),   &
128                     &                                                 2.*ABS( zwy  (ji,jj-1,jk) ),   &
129                     &                                                 2.*ABS( zwy  (ji,jj  ,jk) ) )
130               END DO
131           END DO
132         END DO             ! interior values
133
134        !                                             !-- MUSCL horizontal advective fluxes
135         DO jk = 1, jpkm1                                     ! interior values
136            zdt  = p2dt(jk)
137            DO jj = 2, jpjm1
138               DO ji = fs_2, fs_jpim1   ! vector opt.
139                  ! MUSCL fluxes
140                  z0u = SIGN( 0.5, pun(ji,jj,jk) )
141                  zalpha = 0.5 - z0u
142                  zu  = z0u - 0.5 * pun(ji,jj,jk) * zdt / ( e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) )
143                  zzwx = ptb(ji+1,jj,jk,jn) + zu * zslpx(ji+1,jj,jk)
144                  zzwy = ptb(ji  ,jj,jk,jn) + zu * zslpx(ji  ,jj,jk)
145                  zwx(ji,jj,jk) = pun(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy )
146                  !
147                  z0v = SIGN( 0.5, pvn(ji,jj,jk) )
148                  zalpha = 0.5 - z0v
149                  zv  = z0v - 0.5 * pvn(ji,jj,jk) * zdt / ( e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) )
150                  zzwx = ptb(ji,jj+1,jk,jn) + zv * zslpy(ji,jj+1,jk)
151                  zzwy = ptb(ji,jj  ,jk,jn) + zv * zslpy(ji,jj  ,jk)
152                  zwy(ji,jj,jk) = pvn(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy )
153               END DO
154            END DO
155         END DO
156
157         !!  centered scheme at lateral b.C. if off-shore velocity
158         DO jk = 1, jpkm1
159            DO jj = 2, jpjm1
160               DO ji = fs_2, fs_jpim1   ! vector opt.
161                  IF( umask(ji,jj,jk) == 0. ) THEN
162                     IF( pun(ji+1,jj,jk) > 0. .AND. ji /= jpi ) THEN
163                        zwx(ji+1,jj,jk) = 0.5 * pun(ji+1,jj,jk) * ( ptn(ji+1,jj,jk,jn) + ptn(ji+2,jj,jk,jn) )
164                     ENDIF
165                     IF( pun(ji-1,jj,jk) < 0. ) THEN
166                        zwx(ji-1,jj,jk) = 0.5 * pun(ji-1,jj,jk) * ( ptn(ji-1,jj,jk,jn) + ptn(ji,jj,jk,jn) ) 
167                     ENDIF
168                  ENDIF
169                  IF( vmask(ji,jj,jk) == 0. ) THEN
170                     IF( pvn(ji,jj+1,jk) > 0. .AND. jj /= jpj ) THEN
171                        zwy(ji,jj+1,jk) = 0.5 * pvn(ji,jj+1,jk) * ( ptn(ji,jj+1,jk,jn) + ptn(ji,jj+2,jk,jn) )
172                     ENDIF
173                     IF( pvn(ji,jj-1,jk) < 0. ) THEN
174                        zwy(ji,jj-1,jk) = 0.5 * pvn(ji,jj-1,jk) * ( ptn(ji,jj-1,jk,jn) + ptn(ji,jj,jk,jn) ) 
175                     ENDIF
176                  ENDIF
177               END DO
178            END DO
179         END DO
180         CALL lbc_lnk( zwx, 'U', -1. )   ;   CALL lbc_lnk( zwy, 'V', -1. )   ! lateral boundary condition (changed sign)
181
182         ! Tracer flux divergence at t-point added to the general trend
183         DO jk = 1, jpkm1
184            DO jj = 2, jpjm1
185               DO ji = fs_2, fs_jpim1   ! vector opt.
186                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
187                  ! horizontal advective trends
188                  ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
189                  &               + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  ) )
190                  ! added to the general tracer trends
191                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
192               END DO
193           END DO
194         END DO
195         !                                 ! trend diagnostics (contribution of upstream fluxes)
196         IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR.   &
197            &( cdtype == 'TRC' .AND. l_trdtrc )      ) THEN
198            CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptb(:,:,:,jn) )
199            CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptb(:,:,:,jn) )
200         END IF
201
202         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
203         IF( cdtype == 'TRA' .AND. ln_diaptr ) CALL dia_ptr_ohst_components( jn, 'adv', zwy(:,:,:)  )
204
205         ! II. Vertical advective fluxes
206         ! -----------------------------
207         !                                             !-- first guess of the slopes
208         zwx (:,:, 1 ) = 0.e0    ;    zwx (:,:,jpk) = 0.e0    ! surface & bottom boundary conditions
209         DO jk = 2, jpkm1                                     ! interior values
210            zwx(:,:,jk) = tmask(:,:,jk) * ( ptb(:,:,jk-1,jn) - ptb(:,:,jk,jn) )
211         END DO
212
213         !                                             !-- Slopes of tracer
214         zslpx(:,:,1) = 0.e0                                  ! surface values
215         DO jk = 2, jpkm1                                     ! interior value
216            DO jj = 1, jpj
217               DO ji = 1, jpi
218                  zslpx(ji,jj,jk) =                    ( zwx(ji,jj,jk) + zwx(ji,jj,jk+1) )   &
219                     &            * ( 0.25 + SIGN( 0.25, zwx(ji,jj,jk) * zwx(ji,jj,jk+1) ) )
220               END DO
221            END DO
222         END DO
223         !                                             !-- Slopes limitation
224         DO jk = 2, jpkm1                                     ! interior values
225            DO jj = 1, jpj
226               DO ji = 1, jpi
227                  zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) ) * MIN(    ABS( zslpx(ji,jj,jk  ) ),   &
228                     &                                                 2.*ABS( zwx  (ji,jj,jk+1) ),   &
229                     &                                                 2.*ABS( zwx  (ji,jj,jk  ) )  )
230               END DO
231            END DO
232         END DO
233         !                                             !-- vertical advective flux
234         !                                                    ! surface values  (bottom already set to zero)
235         IF( lk_vvl ) THEN    ;   zwx(:,:, 1 ) = 0.e0                      !  variable volume
236         ELSE                 ;   zwx(:,:, 1 ) = pwn(:,:,1) * ptb(:,:,1,jn)   ! linear free surface
237         ENDIF
238         !
239         DO jk = 1, jpkm1                                     ! interior values
240            zdt  = p2dt(jk)
241            DO jj = 2, jpjm1
242               DO ji = fs_2, fs_jpim1   ! vector opt.
243                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3w(ji,jj,jk+1) )
244                  z0w = SIGN( 0.5, pwn(ji,jj,jk+1) )
245                  zalpha = 0.5 + z0w
246                  zw  = z0w - 0.5 * pwn(ji,jj,jk+1) * zdt * zbtr
247                  zzwx = ptb(ji,jj,jk+1,jn) + zw * zslpx(ji,jj,jk+1)
248                  zzwy = ptb(ji,jj,jk  ,jn) + zw * zslpx(ji,jj,jk  )
249                  zwx(ji,jj,jk+1) = pwn(ji,jj,jk+1) * ( zalpha * zzwx + (1.-zalpha) * zzwy )
250               END DO
251            END DO
252         END DO
253         !
254         DO jk = 2, jpkm1        ! centered near the bottom
255            DO jj = 2, jpjm1
256               DO ji = fs_2, fs_jpim1   ! vector opt.
257                  IF( tmask(ji,jj,jk+1) == 0. ) THEN
258                     IF( pwn(ji,jj,jk) > 0. ) THEN
259                        zwx(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk-1,jn) + ptn(ji,jj,jk,jn) ) 
260                     ENDIF
261                  ENDIF
262               END DO
263            END DO
264         END DO
265         !
266         DO jk = 1, jpkm1        ! Compute & add the vertical advective trend
267            DO jj = 2, jpjm1     
268               DO ji = fs_2, fs_jpim1   ! vector opt.
269                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
270                  ! vertical advective trends
271                  ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) )
272                  ! added to the general tracer trends
273                  pta(ji,jj,jk,jn) =  pta(ji,jj,jk,jn) + ztra
274               END DO
275            END DO
276         END DO
277         !                       ! trend diagnostics (contribution of upstream fluxes)
278         IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR.     &
279            &( cdtype == 'TRC' .AND. l_trdtrc )      )   &
280            CALL trd_tra( kt, cdtype, jn, jptra_zad, zwx, pwn, ptb(:,:,:,jn) )
281         !
282      END DO
283      !
284      CALL wrk_dealloc( jpi, jpj, jpk, zslpx, zslpy, zwx, zwy )
285      !
286      IF( nn_timing == 1 )  CALL timing_stop('tra_adv_muscl2')
287      !
288   END SUBROUTINE tra_adv_muscl2
289
290   !!======================================================================
291END MODULE traadv_muscl2
Note: See TracBrowser for help on using the repository browser.