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

Last change on this file since 11738 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
RevLine 
[3]1MODULE traadv_muscl
[503]2   !!======================================================================
[3]3   !!                       ***  MODULE  traadv_muscl  ***
[2528]4   !! Ocean  tracers:  horizontal & vertical advective trend
[503]5   !!======================================================================
[2528]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
[3680]10   !!            3.4  !  2012-06  (P. Oddo, M. Vichi) include the upstream where needed
[503]11   !!----------------------------------------------------------------------
[3]12
13   !!----------------------------------------------------------------------
14   !!   tra_adv_muscl : update the tracer trend with the horizontal
15   !!                   and vertical advection trends using MUSCL scheme
16   !!----------------------------------------------------------------------
[3625]17   USE oce            ! ocean dynamics and active tracers
[4990]18   USE trc_oce        ! share passive tracers/Ocean variables
[3625]19   USE dom_oce        ! ocean space and time domain
[4990]20   USE trd_oce        ! trends: ocean variables
21   USE trdtra         ! tracers trends manager
[3625]22   USE dynspg_oce     ! choice/control of key cpp for surface pressure gradient
[5147]23   USE sbcrnf         ! river runoffs
[3625]24   USE diaptr         ! poleward transport diagnostics
[4990]25   !
[3625]26   USE wrk_nemo       ! Memory Allocation
27   USE timing         ! Timing
28   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
[4990]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)
[3]32
33   IMPLICIT NONE
34   PRIVATE
35
[4990]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   
[3]42   !! * Substitutions
43#  include "domzgr_substitute.h90"
44#  include "vectopt_loop_substitute.h90"
45   !!----------------------------------------------------------------------
[2528]46   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
[6486]47   !! $Id$
[2528]48   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
[3]49   !!----------------------------------------------------------------------
50CONTAINS
51
[4990]52   SUBROUTINE tra_adv_muscl( kt, kit000, cdtype, p2dt, pun, pvn, pwn,   &
53      &                                                ptb, pta, kjpt, ld_msc_ups )
[3]54      !!----------------------------------------------------------------------
55      !!                    ***  ROUTINE tra_adv_muscl  ***
[216]56      !!
[3]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      !!
[216]61      !! ** Method  : MUSCL scheme plus centered scheme at ocean boundaries
[3]62      !!
63      !! ** Action  : - update (ta,sa) with the now advective tracer trends
[2528]64      !!              - save trends
[3]65      !!
[503]66      !! References : Estubier, A., and M. Levy, Notes Techn. Pole de Modelisation
67      !!              IPSL, Sept. 2000 (http://www.lodyc.jussieu.fr/opa)
68      !!----------------------------------------------------------------------
[2528]69      INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index
[3294]70      INTEGER                              , INTENT(in   ) ::   kit000          ! first time step index
[2528]71      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator)
72      INTEGER                              , INTENT(in   ) ::   kjpt            ! number of tracers
[3718]73      LOGICAL                              , INTENT(in   ) ::   ld_msc_ups      ! use upstream scheme within muscl
[2528]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
[2715]78      !
[4990]79      INTEGER  ::   ji, jj, jk, jn            ! dummy loop indices
80      INTEGER  ::   ierr                      ! local integer
[2715]81      REAL(wp) ::   zu, z0u, zzwx, zw         ! local scalars
82      REAL(wp) ::   zv, z0v, zzwy, z0w        !   -      -
83      REAL(wp) ::   ztra, zbtr, zdt, zalpha   !   -      -
[7771]84      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zslpx, zslpy   ! 3D workspace
85      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zwx  , zwy     ! -      -
[3]86      !!----------------------------------------------------------------------
[3294]87      !
88      IF( nn_timing == 1 )  CALL timing_start('tra_adv_muscl')
89      !
[7771]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) )
[3294]94      !
95      IF( kt == kit000 )  THEN
[2528]96         IF(lwp) WRITE(numout,*)
97         IF(lwp) WRITE(numout,*) 'tra_adv : MUSCL advection scheme on ', cdtype
[3718]98         IF(lwp) WRITE(numout,*) '        : mixed up-stream           ', ld_msc_ups
[2528]99         IF(lwp) WRITE(numout,*) '~~~~~~~'
[3680]100         IF(lwp) WRITE(numout,*)
[11101]101         IF(lwp .AND. lflush) CALL flush(numout)
[2528]102         !
[3680]103         !
[3718]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
[3680]110         ENDIF
111
[3718]112         IF( .NOT. ALLOCATED( xind ) ) THEN
113             ALLOCATE( xind(jpi,jpj,jpk), STAT=ierr )
[3680]114             IF( ierr /= 0 )   CALL ctl_stop('STOP', 'tra_adv_muscl: unable to allocate zind array')
115         ENDIF
116         !
[3]117
[3718]118         !
[4990]119         ! Upstream / MUSCL scheme indicator
[3718]120         ! ------------------------------------
[4990]121!!gm  useless
[3718]122         xind(:,:,:) = 1._wp                             ! set equal to 1 where up-stream is not needed
[4990]123!!gm
[3718]124         !
125         IF( ld_msc_ups )  THEN
[4990]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
[3718]130            END DO
[3680]131         ENDIF 
[3718]132         !
133      ENDIF 
[4990]134      !     
[2528]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
[3]150         END DO
[2528]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
[3]164            END DO
165         END DO
[503]166         !
[2528]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) ) )
[503]176               END DO
[2528]177           END DO
178         END DO             ! interior values
[216]179
[2528]180         !                                             !-- MUSCL horizontal advective fluxes
181         DO jk = 1, jpkm1                                     ! interior values
182            zdt  = p2dt(jk)
[503]183            DO jj = 2, jpjm1
184               DO ji = fs_2, fs_jpim1   ! vector opt.
[2528]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) )
[4990]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)
[2528]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) )
[4990]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)
[2528]198                  zwy(ji,jj,jk) = pvn(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy )
[503]199               END DO
200            END DO
201         END DO
[2528]202         !                                                    ! lateral boundary conditions on zwx, zwy   (changed sign)
203         CALL lbc_lnk( zwx, 'U', -1. )   ;   CALL lbc_lnk( zwy, 'V', -1. )
[503]204         !
[2528]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
[3]215               END DO
[2528]216           END DO
217         END DO       
218         !                                 ! trend diagnostics (contribution of upstream fluxes)
[4990]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) )
[2528]223         END IF
224         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
[7179]225         IF( cdtype == 'TRA' .AND. ln_diaptr )  CALL dia_ptr_ohst_components( jn, 'adv', zwy(:,:,:)  )
[3]226
[2528]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) )
[3]233         END DO
234
[2528]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
[3]243            END DO
244         END DO
[2528]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
[3]253            END DO
254         END DO
[2528]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 
[4990]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  )
[2528]271                  zwx(ji,jj,jk+1) = pwn(ji,jj,jk+1) * ( zalpha * zzwx + (1.-zalpha) * zzwy )
272               END DO
[3]273            END DO
274         END DO
275
[4990]276         DO jk = 1, jpkm1                    ! Compute & add the vertical advective trend
[2528]277            DO jj = 2, jpjm1     
[503]278               DO ji = fs_2, fs_jpim1   ! vector opt.
[4990]279                  zbtr = 1. / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) )
[2528]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
[503]284               END DO
285            END DO
286         END DO
[2528]287         !                                 ! Save the vertical advective trends for diagnostic
[4990]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) )
[503]291         !
[4990]292      END DO
[503]293      !
[7771]294      DEALLOCATE( zslpx )
295      DEALLOCATE( zslpy )
296      DEALLOCATE( zwx   )
297      DEALLOCATE( zwy   )
[2715]298      !
[3294]299      IF( nn_timing == 1 )  CALL timing_stop('tra_adv_muscl')
300      !
[3]301   END SUBROUTINE tra_adv_muscl
302
303   !!======================================================================
304END MODULE traadv_muscl
Note: See TracBrowser for help on using the repository browser.