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

Last change on this file since 8485 was 7771, checked in by frrh, 7 years ago

Apply optimisations to various areas of code replacing the use of
allocated pointers with straightforward direct ALLOCATE and DEALLOCATE
operations.

These optimisations largely have an impact in models featuring MEDUSA,
i.e. those with significant numbers of tracers, although they are
expected to have a small impact in all configurations.

Code developed and tested in NEMO branch branches/UKMO/dev_r5518_optim_GO6_alloc
Tested in stand-alone GO6-GSI8, GO6-GSI8-MEDUSA and UKESM coupled models.
NEMO ticket #1821 documents this change further.
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,*)
[2528]101         !
[3680]102         !
[3718]103         IF( ld_msc_ups ) THEN
104            IF( .NOT. ALLOCATED( upsmsk ) )  THEN
105                ALLOCATE( upsmsk(jpi,jpj), STAT=ierr )
106                IF( ierr /= 0 )   CALL ctl_stop('STOP', 'tra_adv_muscl: unable to allocate upsmsk array')
107            ENDIF
108            upsmsk(:,:) = 0._wp                             ! not upstream by default
[3680]109         ENDIF
110
[3718]111         IF( .NOT. ALLOCATED( xind ) ) THEN
112             ALLOCATE( xind(jpi,jpj,jpk), STAT=ierr )
[3680]113             IF( ierr /= 0 )   CALL ctl_stop('STOP', 'tra_adv_muscl: unable to allocate zind array')
114         ENDIF
115         !
[3]116
[3718]117         !
[4990]118         ! Upstream / MUSCL scheme indicator
[3718]119         ! ------------------------------------
[4990]120!!gm  useless
[3718]121         xind(:,:,:) = 1._wp                             ! set equal to 1 where up-stream is not needed
[4990]122!!gm
[3718]123         !
124         IF( ld_msc_ups )  THEN
[4990]125            DO jk = 1, jpkm1
126               xind(:,:,jk) = 1._wp                              &                 ! =>1 where up-stream is not needed
127                  &         - MAX ( rnfmsk(:,:) * rnfmsk_z(jk),  &                 ! =>0 near runoff mouths (& closed sea outflows)
128                  &                 upsmsk(:,:)                ) * tmask(:,:,jk)   ! =>0 near some straits
[3718]129            END DO
[3680]130         ENDIF 
[3718]131         !
132      ENDIF 
[4990]133      !     
[2528]134      !                                                     ! ===========
135      DO jn = 1, kjpt                                       ! tracer loop
136         !                                                  ! ===========
137         ! I. Horizontal advective fluxes
138         ! ------------------------------
139         ! first guess of the slopes
140         zwx(:,:,jpk) = 0.e0   ;   zwy(:,:,jpk) = 0.e0        ! bottom values
141         ! interior values
142         DO jk = 1, jpkm1
143            DO jj = 1, jpjm1     
144               DO ji = 1, fs_jpim1   ! vector opt.
145                  zwx(ji,jj,jk) = umask(ji,jj,jk) * ( ptb(ji+1,jj,jk,jn) - ptb(ji,jj,jk,jn) )
146                  zwy(ji,jj,jk) = vmask(ji,jj,jk) * ( ptb(ji,jj+1,jk,jn) - ptb(ji,jj,jk,jn) )
147               END DO
148           END DO
[3]149         END DO
[2528]150         !
151         CALL lbc_lnk( zwx, 'U', -1. )                        ! lateral boundary conditions on zwx, zwy   (changed sign)
152         CALL lbc_lnk( zwy, 'V', -1. )
153         !                                             !-- Slopes of tracer
154         zslpx(:,:,jpk) = 0.e0   ;   zslpy(:,:,jpk) = 0.e0    ! bottom values
155         DO jk = 1, jpkm1                                     ! interior values
156            DO jj = 2, jpj
157               DO ji = fs_2, jpi   ! vector opt.
158                  zslpx(ji,jj,jk) =                    ( zwx(ji,jj,jk) + zwx(ji-1,jj  ,jk) )   &
159                     &            * ( 0.25 + SIGN( 0.25, zwx(ji,jj,jk) * zwx(ji-1,jj  ,jk) ) )
160                  zslpy(ji,jj,jk) =                    ( zwy(ji,jj,jk) + zwy(ji  ,jj-1,jk) )   &
161                     &            * ( 0.25 + SIGN( 0.25, zwy(ji,jj,jk) * zwy(ji  ,jj-1,jk) ) )
162               END DO
[3]163            END DO
164         END DO
[503]165         !
[2528]166         DO jk = 1, jpkm1                                     ! Slopes limitation
167            DO jj = 2, jpj
168               DO ji = fs_2, jpi   ! vector opt.
169                  zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) ) * MIN(    ABS( zslpx(ji  ,jj,jk) ),   &
170                     &                                                 2.*ABS( zwx  (ji-1,jj,jk) ),   &
171                     &                                                 2.*ABS( zwx  (ji  ,jj,jk) ) )
172                  zslpy(ji,jj,jk) = SIGN( 1., zslpy(ji,jj,jk) ) * MIN(    ABS( zslpy(ji,jj  ,jk) ),   &
173                     &                                                 2.*ABS( zwy  (ji,jj-1,jk) ),   &
174                     &                                                 2.*ABS( zwy  (ji,jj  ,jk) ) )
[503]175               END DO
[2528]176           END DO
177         END DO             ! interior values
[216]178
[2528]179         !                                             !-- MUSCL horizontal advective fluxes
180         DO jk = 1, jpkm1                                     ! interior values
181            zdt  = p2dt(jk)
[503]182            DO jj = 2, jpjm1
183               DO ji = fs_2, fs_jpim1   ! vector opt.
[2528]184                  ! MUSCL fluxes
185                  z0u = SIGN( 0.5, pun(ji,jj,jk) )
186                  zalpha = 0.5 - z0u
187                  zu  = z0u - 0.5 * pun(ji,jj,jk) * zdt / ( e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) )
[4990]188                  zzwx = ptb(ji+1,jj,jk,jn) + xind(ji,jj,jk) * zu * zslpx(ji+1,jj,jk)
189                  zzwy = ptb(ji  ,jj,jk,jn) + xind(ji,jj,jk) * zu * zslpx(ji  ,jj,jk)
[2528]190                  zwx(ji,jj,jk) = pun(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy )
191                  !
192                  z0v = SIGN( 0.5, pvn(ji,jj,jk) )
193                  zalpha = 0.5 - z0v
194                  zv  = z0v - 0.5 * pvn(ji,jj,jk) * zdt / ( e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) )
[4990]195                  zzwx = ptb(ji,jj+1,jk,jn) + xind(ji,jj,jk) * zv * zslpy(ji,jj+1,jk)
196                  zzwy = ptb(ji,jj  ,jk,jn) + xind(ji,jj,jk) * zv * zslpy(ji,jj  ,jk)
[2528]197                  zwy(ji,jj,jk) = pvn(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy )
[503]198               END DO
199            END DO
200         END DO
[2528]201         !                                                    ! lateral boundary conditions on zwx, zwy   (changed sign)
202         CALL lbc_lnk( zwx, 'U', -1. )   ;   CALL lbc_lnk( zwy, 'V', -1. )
[503]203         !
[2528]204         ! Tracer flux divergence at t-point added to the general trend
205         DO jk = 1, jpkm1
206            DO jj = 2, jpjm1     
207               DO ji = fs_2, fs_jpim1   ! vector opt.
208                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
209                  ! horizontal advective trends
210                  ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
211                  &               + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  ) )
212                  ! add it to the general tracer trends
213                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
[3]214               END DO
[2528]215           END DO
216         END DO       
217         !                                 ! trend diagnostics (contribution of upstream fluxes)
[4990]218         IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR.   &
219            &( cdtype == 'TRC' .AND. l_trdtrc )      )  THEN
220            CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptb(:,:,:,jn) )
221            CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptb(:,:,:,jn) )
[2528]222         END IF
223         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
[7179]224         IF( cdtype == 'TRA' .AND. ln_diaptr )  CALL dia_ptr_ohst_components( jn, 'adv', zwy(:,:,:)  )
[3]225
[2528]226         ! II. Vertical advective fluxes
227         ! -----------------------------
228         !                                             !-- first guess of the slopes
229         zwx (:,:, 1 ) = 0.e0    ;    zwx (:,:,jpk) = 0.e0    ! surface & bottom boundary conditions
230         DO jk = 2, jpkm1                                     ! interior values
231            zwx(:,:,jk) = tmask(:,:,jk) * ( ptb(:,:,jk-1,jn) - ptb(:,:,jk,jn) )
[3]232         END DO
233
[2528]234         !                                             !-- Slopes of tracer
235         zslpx(:,:,1) = 0.e0                                  ! surface values
236         DO jk = 2, jpkm1                                     ! interior value
237            DO jj = 1, jpj
238               DO ji = 1, jpi
239                  zslpx(ji,jj,jk) =                    ( zwx(ji,jj,jk) + zwx(ji,jj,jk+1) )   &
240                     &            * ( 0.25 + SIGN( 0.25, zwx(ji,jj,jk) * zwx(ji,jj,jk+1) ) )
241               END DO
[3]242            END DO
243         END DO
[2528]244         !                                             !-- Slopes limitation
245         DO jk = 2, jpkm1                                     ! interior values
246            DO jj = 1, jpj
247               DO ji = 1, jpi
248                  zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) ) * MIN(    ABS( zslpx(ji,jj,jk  ) ),   &
249                     &                                                 2.*ABS( zwx  (ji,jj,jk+1) ),   &
250                     &                                                 2.*ABS( zwx  (ji,jj,jk  ) )  )
251               END DO
[3]252            END DO
253         END DO
[2528]254         !                                             !-- vertical advective flux
255         !                                                    ! surface values  (bottom already set to zero)
256         IF( lk_vvl ) THEN    ;   zwx(:,:, 1 ) = 0.e0                      !  variable volume
257         ELSE                 ;   zwx(:,:, 1 ) = pwn(:,:,1) * ptb(:,:,1,jn)   ! linear free surface
258         ENDIF 
259         !
260         DO jk = 1, jpkm1                                     ! interior values
261            zdt  = p2dt(jk)
262            DO jj = 2, jpjm1     
263               DO ji = fs_2, fs_jpim1   ! vector opt.
264                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3w(ji,jj,jk+1) )
265                  z0w = SIGN( 0.5, pwn(ji,jj,jk+1) )
266                  zalpha = 0.5 + z0w
267                  zw  = z0w - 0.5 * pwn(ji,jj,jk+1) * zdt * zbtr 
[4990]268                  zzwx = ptb(ji,jj,jk+1,jn) + xind(ji,jj,jk) * zw * zslpx(ji,jj,jk+1)
269                  zzwy = ptb(ji,jj,jk  ,jn) + xind(ji,jj,jk) * zw * zslpx(ji,jj,jk  )
[2528]270                  zwx(ji,jj,jk+1) = pwn(ji,jj,jk+1) * ( zalpha * zzwx + (1.-zalpha) * zzwy )
271               END DO
[3]272            END DO
273         END DO
274
[4990]275         DO jk = 1, jpkm1                    ! Compute & add the vertical advective trend
[2528]276            DO jj = 2, jpjm1     
[503]277               DO ji = fs_2, fs_jpim1   ! vector opt.
[4990]278                  zbtr = 1. / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) )
[2528]279                  ! vertical advective trends
280                  ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) )
281                  ! add it to the general tracer trends
282                  pta(ji,jj,jk,jn) =  pta(ji,jj,jk,jn) + ztra
[503]283               END DO
284            END DO
285         END DO
[2528]286         !                                 ! Save the vertical advective trends for diagnostic
[4990]287         IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR.     &
288            &( cdtype == 'TRC' .AND. l_trdtrc )      )   &
289            CALL trd_tra( kt, cdtype, jn, jptra_zad, zwx, pwn, ptb(:,:,:,jn) )
[503]290         !
[4990]291      END DO
[503]292      !
[7771]293      DEALLOCATE( zslpx )
294      DEALLOCATE( zslpy )
295      DEALLOCATE( zwx   )
296      DEALLOCATE( zwy   )
[2715]297      !
[3294]298      IF( nn_timing == 1 )  CALL timing_stop('tra_adv_muscl')
299      !
[3]300   END SUBROUTINE tra_adv_muscl
301
302   !!======================================================================
303END MODULE traadv_muscl
Note: See TracBrowser for help on using the repository browser.