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/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_muscl.F90

Last change on this file 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: 16.0 KB
Line 
1MODULE traadv_muscl
2   !!======================================================================
3   !!                       ***  MODULE  traadv_muscl  ***
4   !! Ocean  tracers:  horizontal & vertical advective trend
5   !!======================================================================
6   !! History :       !  2000-06  (A.Estublier)  for passive tracers
7   !!                 !  2001-08  (E.Durand, G.Madec)  adapted for T & S
8   !!   NEMO     1.0  !  2002-06  (G. Madec)  F90: Free form and module
9   !!            3.2  !  2010-05  (C. Ethe, G. Madec)  merge TRC-TRA + switch from velocity to transport
10   !!            3.4  !  2012-06  (P. Oddo, M. Vichi) include the upstream where needed
11   !!----------------------------------------------------------------------
12
13   !!----------------------------------------------------------------------
14   !!   tra_adv_muscl : update the tracer trend with the horizontal
15   !!                   and vertical advection trends using MUSCL scheme
16   !!----------------------------------------------------------------------
17   USE oce            ! ocean dynamics and active tracers
18   USE trc_oce        ! share passive tracers/Ocean variables
19   USE dom_oce        ! ocean space and time domain
20   USE trd_oce        ! trends: ocean variables
21   USE trdtra         ! tracers trends manager
22   USE dynspg_oce     ! choice/control of key cpp for surface pressure gradient
23   USE sbcrnf         ! river runoffs
24   USE diaptr         ! poleward transport diagnostics
25   !
26   USE wrk_nemo       ! Memory Allocation
27   USE timing         ! Timing
28   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
29   USE in_out_manager ! I/O manager
30   USE lib_mpp        ! distribued memory computing
31   USE lbclnk         ! ocean lateral boundary condition (or mpp link)
32
33   IMPLICIT NONE
34   PRIVATE
35
36   PUBLIC   tra_adv_muscl   ! routine called by traadv.F90
37   
38   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   upsmsk   !: mixed upstream/centered scheme near some straits
39   !                                                           !  and in closed seas (orca 2 and 4 configurations)
40   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   xind     !: mixed upstream/centered index
41   
42   !! * Substitutions
43#  include "domzgr_substitute.h90"
44#  include "vectopt_loop_substitute.h90"
45   !!----------------------------------------------------------------------
46   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
47   !! $Id$
48   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
49   !!----------------------------------------------------------------------
50CONTAINS
51
52   SUBROUTINE tra_adv_muscl( kt, kit000, cdtype, p2dt, pun, pvn, pwn,   &
53      &                                                ptb, pta, kjpt, ld_msc_ups )
54      !!----------------------------------------------------------------------
55      !!                    ***  ROUTINE tra_adv_muscl  ***
56      !!
57      !! ** Purpose :   Compute the now trend due to total advection of T and
58      !!      S using a MUSCL scheme (Monotone Upstream-centered Scheme for
59      !!      Conservation Laws) and add it to the general tracer trend.
60      !!
61      !! ** Method  : MUSCL scheme plus centered scheme at ocean boundaries
62      !!
63      !! ** Action  : - update (ta,sa) with the now advective tracer trends
64      !!              - save trends
65      !!
66      !! References : Estubier, A., and M. Levy, Notes Techn. Pole de Modelisation
67      !!              IPSL, Sept. 2000 (http://www.lodyc.jussieu.fr/opa)
68      !!----------------------------------------------------------------------
69      INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index
70      INTEGER                              , INTENT(in   ) ::   kit000          ! first time step index
71      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator)
72      INTEGER                              , INTENT(in   ) ::   kjpt            ! number of tracers
73      LOGICAL                              , INTENT(in   ) ::   ld_msc_ups      ! use upstream scheme within muscl
74      REAL(wp), DIMENSION(        jpk     ), INTENT(in   ) ::   p2dt            ! vertical profile of tracer time-step
75      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pun, pvn, pwn   ! 3 ocean velocity components
76      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb             ! before tracer field
77      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta             ! tracer trend
78      !
79      INTEGER  ::   ji, jj, jk, jn            ! dummy loop indices
80      INTEGER  ::   ierr                      ! local integer
81      REAL(wp) ::   zu, z0u, zzwx, zw         ! local scalars
82      REAL(wp) ::   zv, z0v, zzwy, z0w        !   -      -
83      REAL(wp) ::   ztra, zbtr, zdt, zalpha   !   -      -
84      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zslpx, zslpy   ! 3D workspace
85      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zwx  , zwy     ! -      -
86      !!----------------------------------------------------------------------
87      !
88      IF( nn_timing == 1 )  CALL timing_start('tra_adv_muscl')
89      !
90      ALLOCATE( zslpx(1:jpi, 1:jpj, 1:jpk) )
91      ALLOCATE( zslpy(1:jpi, 1:jpj, 1:jpk) )
92      ALLOCATE( zwx  (1:jpi, 1:jpj, 1:jpk) )
93      ALLOCATE( zwy  (1:jpi, 1:jpj, 1:jpk) )
94      !
95      IF( kt == kit000 )  THEN
96         IF(lwp) WRITE(numout,*)
97         IF(lwp) WRITE(numout,*) 'tra_adv : MUSCL advection scheme on ', cdtype
98         IF(lwp) WRITE(numout,*) '        : mixed up-stream           ', ld_msc_ups
99         IF(lwp) WRITE(numout,*) '~~~~~~~'
100         IF(lwp) WRITE(numout,*)
101         IF(lwp .AND. lflush) CALL flush(numout)
102         !
103         !
104         IF( ld_msc_ups ) THEN
105            IF( .NOT. ALLOCATED( upsmsk ) )  THEN
106                ALLOCATE( upsmsk(jpi,jpj), STAT=ierr )
107                IF( ierr /= 0 )   CALL ctl_stop('STOP', 'tra_adv_muscl: unable to allocate upsmsk array')
108            ENDIF
109            upsmsk(:,:) = 0._wp                             ! not upstream by default
110         ENDIF
111
112         IF( .NOT. ALLOCATED( xind ) ) THEN
113             ALLOCATE( xind(jpi,jpj,jpk), STAT=ierr )
114             IF( ierr /= 0 )   CALL ctl_stop('STOP', 'tra_adv_muscl: unable to allocate zind array')
115         ENDIF
116         !
117
118         !
119         ! Upstream / MUSCL scheme indicator
120         ! ------------------------------------
121!!gm  useless
122         xind(:,:,:) = 1._wp                             ! set equal to 1 where up-stream is not needed
123!!gm
124         !
125         IF( ld_msc_ups )  THEN
126            DO jk = 1, jpkm1
127               xind(:,:,jk) = 1._wp                              &                 ! =>1 where up-stream is not needed
128                  &         - MAX ( rnfmsk(:,:) * rnfmsk_z(jk),  &                 ! =>0 near runoff mouths (& closed sea outflows)
129                  &                 upsmsk(:,:)                ) * tmask(:,:,jk)   ! =>0 near some straits
130            END DO
131         ENDIF 
132         !
133      ENDIF 
134      !     
135      !                                                     ! ===========
136      DO jn = 1, kjpt                                       ! tracer loop
137         !                                                  ! ===========
138         ! I. Horizontal advective fluxes
139         ! ------------------------------
140         ! first guess of the slopes
141         zwx(:,:,jpk) = 0.e0   ;   zwy(:,:,jpk) = 0.e0        ! bottom values
142         ! interior values
143         DO jk = 1, jpkm1
144            DO jj = 1, jpjm1     
145               DO ji = 1, fs_jpim1   ! vector opt.
146                  zwx(ji,jj,jk) = umask(ji,jj,jk) * ( ptb(ji+1,jj,jk,jn) - ptb(ji,jj,jk,jn) )
147                  zwy(ji,jj,jk) = vmask(ji,jj,jk) * ( ptb(ji,jj+1,jk,jn) - ptb(ji,jj,jk,jn) )
148               END DO
149           END DO
150         END DO
151         !
152         CALL lbc_lnk( zwx, 'U', -1. )                        ! lateral boundary conditions on zwx, zwy   (changed sign)
153         CALL lbc_lnk( zwy, 'V', -1. )
154         !                                             !-- Slopes of tracer
155         zslpx(:,:,jpk) = 0.e0   ;   zslpy(:,:,jpk) = 0.e0    ! bottom values
156         DO jk = 1, jpkm1                                     ! interior values
157            DO jj = 2, jpj
158               DO ji = fs_2, jpi   ! vector opt.
159                  zslpx(ji,jj,jk) =                    ( zwx(ji,jj,jk) + zwx(ji-1,jj  ,jk) )   &
160                     &            * ( 0.25 + SIGN( 0.25, zwx(ji,jj,jk) * zwx(ji-1,jj  ,jk) ) )
161                  zslpy(ji,jj,jk) =                    ( zwy(ji,jj,jk) + zwy(ji  ,jj-1,jk) )   &
162                     &            * ( 0.25 + SIGN( 0.25, zwy(ji,jj,jk) * zwy(ji  ,jj-1,jk) ) )
163               END DO
164            END DO
165         END DO
166         !
167         DO jk = 1, jpkm1                                     ! Slopes limitation
168            DO jj = 2, jpj
169               DO ji = fs_2, jpi   ! vector opt.
170                  zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) ) * MIN(    ABS( zslpx(ji  ,jj,jk) ),   &
171                     &                                                 2.*ABS( zwx  (ji-1,jj,jk) ),   &
172                     &                                                 2.*ABS( zwx  (ji  ,jj,jk) ) )
173                  zslpy(ji,jj,jk) = SIGN( 1., zslpy(ji,jj,jk) ) * MIN(    ABS( zslpy(ji,jj  ,jk) ),   &
174                     &                                                 2.*ABS( zwy  (ji,jj-1,jk) ),   &
175                     &                                                 2.*ABS( zwy  (ji,jj  ,jk) ) )
176               END DO
177           END DO
178         END DO             ! interior values
179
180         !                                             !-- MUSCL horizontal advective fluxes
181         DO jk = 1, jpkm1                                     ! interior values
182            zdt  = p2dt(jk)
183            DO jj = 2, jpjm1
184               DO ji = fs_2, fs_jpim1   ! vector opt.
185                  ! MUSCL fluxes
186                  z0u = SIGN( 0.5, pun(ji,jj,jk) )
187                  zalpha = 0.5 - z0u
188                  zu  = z0u - 0.5 * pun(ji,jj,jk) * zdt / ( e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) )
189                  zzwx = ptb(ji+1,jj,jk,jn) + xind(ji,jj,jk) * zu * zslpx(ji+1,jj,jk)
190                  zzwy = ptb(ji  ,jj,jk,jn) + xind(ji,jj,jk) * zu * zslpx(ji  ,jj,jk)
191                  zwx(ji,jj,jk) = pun(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy )
192                  !
193                  z0v = SIGN( 0.5, pvn(ji,jj,jk) )
194                  zalpha = 0.5 - z0v
195                  zv  = z0v - 0.5 * pvn(ji,jj,jk) * zdt / ( e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) )
196                  zzwx = ptb(ji,jj+1,jk,jn) + xind(ji,jj,jk) * zv * zslpy(ji,jj+1,jk)
197                  zzwy = ptb(ji,jj  ,jk,jn) + xind(ji,jj,jk) * zv * zslpy(ji,jj  ,jk)
198                  zwy(ji,jj,jk) = pvn(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy )
199               END DO
200            END DO
201         END DO
202         !                                                    ! lateral boundary conditions on zwx, zwy   (changed sign)
203         CALL lbc_lnk( zwx, 'U', -1. )   ;   CALL lbc_lnk( zwy, 'V', -1. )
204         !
205         ! Tracer flux divergence at t-point added to the general trend
206         DO jk = 1, jpkm1
207            DO jj = 2, jpjm1     
208               DO ji = fs_2, fs_jpim1   ! vector opt.
209                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
210                  ! horizontal advective trends
211                  ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
212                  &               + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  ) )
213                  ! add it to the general tracer trends
214                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
215               END DO
216           END DO
217         END DO       
218         !                                 ! trend diagnostics (contribution of upstream fluxes)
219         IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR.   &
220            &( cdtype == 'TRC' .AND. l_trdtrc )      )  THEN
221            CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptb(:,:,:,jn) )
222            CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptb(:,:,:,jn) )
223         END IF
224         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
225         IF( cdtype == 'TRA' .AND. ln_diaptr )  CALL dia_ptr_ohst_components( jn, 'adv', zwy(:,:,:)  )
226
227         ! II. Vertical advective fluxes
228         ! -----------------------------
229         !                                             !-- first guess of the slopes
230         zwx (:,:, 1 ) = 0.e0    ;    zwx (:,:,jpk) = 0.e0    ! surface & bottom boundary conditions
231         DO jk = 2, jpkm1                                     ! interior values
232            zwx(:,:,jk) = tmask(:,:,jk) * ( ptb(:,:,jk-1,jn) - ptb(:,:,jk,jn) )
233         END DO
234
235         !                                             !-- Slopes of tracer
236         zslpx(:,:,1) = 0.e0                                  ! surface values
237         DO jk = 2, jpkm1                                     ! interior value
238            DO jj = 1, jpj
239               DO ji = 1, jpi
240                  zslpx(ji,jj,jk) =                    ( zwx(ji,jj,jk) + zwx(ji,jj,jk+1) )   &
241                     &            * ( 0.25 + SIGN( 0.25, zwx(ji,jj,jk) * zwx(ji,jj,jk+1) ) )
242               END DO
243            END DO
244         END DO
245         !                                             !-- Slopes limitation
246         DO jk = 2, jpkm1                                     ! interior values
247            DO jj = 1, jpj
248               DO ji = 1, jpi
249                  zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) ) * MIN(    ABS( zslpx(ji,jj,jk  ) ),   &
250                     &                                                 2.*ABS( zwx  (ji,jj,jk+1) ),   &
251                     &                                                 2.*ABS( zwx  (ji,jj,jk  ) )  )
252               END DO
253            END DO
254         END DO
255         !                                             !-- vertical advective flux
256         !                                                    ! surface values  (bottom already set to zero)
257         IF( lk_vvl ) THEN    ;   zwx(:,:, 1 ) = 0.e0                      !  variable volume
258         ELSE                 ;   zwx(:,:, 1 ) = pwn(:,:,1) * ptb(:,:,1,jn)   ! linear free surface
259         ENDIF 
260         !
261         DO jk = 1, jpkm1                                     ! interior values
262            zdt  = p2dt(jk)
263            DO jj = 2, jpjm1     
264               DO ji = fs_2, fs_jpim1   ! vector opt.
265                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3w(ji,jj,jk+1) )
266                  z0w = SIGN( 0.5, pwn(ji,jj,jk+1) )
267                  zalpha = 0.5 + z0w
268                  zw  = z0w - 0.5 * pwn(ji,jj,jk+1) * zdt * zbtr 
269                  zzwx = ptb(ji,jj,jk+1,jn) + xind(ji,jj,jk) * zw * zslpx(ji,jj,jk+1)
270                  zzwy = ptb(ji,jj,jk  ,jn) + xind(ji,jj,jk) * zw * zslpx(ji,jj,jk  )
271                  zwx(ji,jj,jk+1) = pwn(ji,jj,jk+1) * ( zalpha * zzwx + (1.-zalpha) * zzwy )
272               END DO
273            END DO
274         END DO
275
276         DO jk = 1, jpkm1                    ! Compute & add the vertical advective trend
277            DO jj = 2, jpjm1     
278               DO ji = fs_2, fs_jpim1   ! vector opt.
279                  zbtr = 1. / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) )
280                  ! vertical advective trends
281                  ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) )
282                  ! add it to the general tracer trends
283                  pta(ji,jj,jk,jn) =  pta(ji,jj,jk,jn) + ztra
284               END DO
285            END DO
286         END DO
287         !                                 ! Save the vertical advective trends for diagnostic
288         IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR.     &
289            &( cdtype == 'TRC' .AND. l_trdtrc )      )   &
290            CALL trd_tra( kt, cdtype, jn, jptra_zad, zwx, pwn, ptb(:,:,:,jn) )
291         !
292      END DO
293      !
294      DEALLOCATE( zslpx )
295      DEALLOCATE( zslpy )
296      DEALLOCATE( zwx   )
297      DEALLOCATE( zwy   )
298      !
299      IF( nn_timing == 1 )  CALL timing_stop('tra_adv_muscl')
300      !
301   END SUBROUTINE tra_adv_muscl
302
303   !!======================================================================
304END MODULE traadv_muscl
Note: See TracBrowser for help on using the repository browser.