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

source: branches/UKMO/dev_r5518_GO6_under_ice_relax_dr_hook/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_muscl.F90 @ 11738

Last change on this file since 11738 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
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         !
102         !
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
109         ENDIF
110
111         IF( .NOT. ALLOCATED( xind ) ) THEN
112             ALLOCATE( xind(jpi,jpj,jpk), STAT=ierr )
113             IF( ierr /= 0 )   CALL ctl_stop('STOP', 'tra_adv_muscl: unable to allocate zind array')
114         ENDIF
115         !
116
117         !
118         ! Upstream / MUSCL scheme indicator
119         ! ------------------------------------
120!!gm  useless
121         xind(:,:,:) = 1._wp                             ! set equal to 1 where up-stream is not needed
122!!gm
123         !
124         IF( ld_msc_ups )  THEN
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
129            END DO
130         ENDIF 
131         !
132      ENDIF 
133      !     
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
149         END DO
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
163            END DO
164         END DO
165         !
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) ) )
175               END DO
176           END DO
177         END DO             ! interior values
178
179         !                                             !-- MUSCL horizontal advective fluxes
180         DO jk = 1, jpkm1                                     ! interior values
181            zdt  = p2dt(jk)
182            DO jj = 2, jpjm1
183               DO ji = fs_2, fs_jpim1   ! vector opt.
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) )
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)
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) )
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)
197                  zwy(ji,jj,jk) = pvn(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy )
198               END DO
199            END DO
200         END DO
201         !                                                    ! lateral boundary conditions on zwx, zwy   (changed sign)
202         CALL lbc_lnk( zwx, 'U', -1. )   ;   CALL lbc_lnk( zwy, 'V', -1. )
203         !
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
214               END DO
215           END DO
216         END DO       
217         !                                 ! trend diagnostics (contribution of upstream fluxes)
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) )
222         END IF
223         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
224         IF( cdtype == 'TRA' .AND. ln_diaptr )  CALL dia_ptr_ohst_components( jn, 'adv', zwy(:,:,:)  )
225
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) )
232         END DO
233
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
242            END DO
243         END DO
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
252            END DO
253         END DO
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 
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  )
270                  zwx(ji,jj,jk+1) = pwn(ji,jj,jk+1) * ( zalpha * zzwx + (1.-zalpha) * zzwy )
271               END DO
272            END DO
273         END DO
274
275         DO jk = 1, jpkm1                    ! Compute & add the vertical advective trend
276            DO jj = 2, jpjm1     
277               DO ji = fs_2, fs_jpim1   ! vector opt.
278                  zbtr = 1. / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) )
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
283               END DO
284            END DO
285         END DO
286         !                                 ! Save the vertical advective trends for diagnostic
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) )
290         !
291      END DO
292      !
293      DEALLOCATE( zslpx )
294      DEALLOCATE( zslpy )
295      DEALLOCATE( zwx   )
296      DEALLOCATE( zwy   )
297      !
298      IF( nn_timing == 1 )  CALL timing_stop('tra_adv_muscl')
299      !
300   END SUBROUTINE tra_adv_muscl
301
302   !!======================================================================
303END MODULE traadv_muscl
Note: See TracBrowser for help on using the repository browser.