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_mus.F90 in NEMO/trunk/src/OCE/TRA – NEMO

source: NEMO/trunk/src/OCE/TRA/traadv_mus.F90 @ 12808

Last change on this file since 12808 was 12377, checked in by acc, 4 years ago

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge --ignore-ancestry \

​svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
​svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The --ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

2020 developments can be started from any trunk revision later than this one.
  • Property svn:keywords set to Id
File size: 13.4 KB
RevLine 
[5770]1MODULE traadv_mus
[503]2   !!======================================================================
[5770]3   !!                       ***  MODULE  traadv_mus  ***
[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
[5770]11   !!            3.7  !  2015-09  (G. Madec) add the ice-shelf cavities boundary condition
[503]12   !!----------------------------------------------------------------------
[3]13
14   !!----------------------------------------------------------------------
[5770]15   !!   tra_adv_mus   : update the tracer trend with the horizontal
[3]16   !!                   and vertical advection trends using MUSCL scheme
17   !!----------------------------------------------------------------------
[3625]18   USE oce            ! ocean dynamics and active tracers
[4990]19   USE trc_oce        ! share passive tracers/Ocean variables
[3625]20   USE dom_oce        ! ocean space and time domain
[4990]21   USE trd_oce        ! trends: ocean variables
22   USE trdtra         ! tracers trends manager
[5147]23   USE sbcrnf         ! river runoffs
[3625]24   USE diaptr         ! poleward transport diagnostics
[7646]25   USE diaar5         ! AR5 diagnostics
26
[4990]27   !
[9019]28   USE iom            ! XIOS library
[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)
[9124]32   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
[3]33
34   IMPLICIT NONE
35   PRIVATE
36
[5770]37   PUBLIC   tra_adv_mus   ! routine called by traadv.F90
[4990]38   
39   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   upsmsk   !: mixed upstream/centered scheme near some straits
[7646]40   !                                                           !  and in closed seas (orca 2 and 1 configurations)
[4990]41   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   xind     !: mixed upstream/centered index
42   
[7646]43   LOGICAL  ::   l_trd   ! flag to compute trends
44   LOGICAL  ::   l_ptr   ! flag to compute poleward transport
45   LOGICAL  ::   l_hst   ! flag to compute heat/salt transport
46
[3]47   !! * Substitutions
[12377]48#  include "do_loop_substitute.h90"
[3]49   !!----------------------------------------------------------------------
[9598]50   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
[1152]51   !! $Id$
[10068]52   !! Software governed by the CeCILL license (see ./LICENSE)
[3]53   !!----------------------------------------------------------------------
54CONTAINS
55
[12377]56   SUBROUTINE tra_adv_mus( kt, kit000, cdtype, p2dt, pU, pV, pW,             &
57      &                    Kbb, Kmm, pt, kjpt, Krhs, ld_msc_ups )
[3]58      !!----------------------------------------------------------------------
[5770]59      !!                    ***  ROUTINE tra_adv_mus  ***
[216]60      !!
[5770]61      !! ** Purpose :   Compute the now trend due to total advection of tracers
62      !!              using a MUSCL scheme (Monotone Upstream-centered Scheme for
63      !!              Conservation Laws) and add it to the general tracer trend.
[3]64      !!
[216]65      !! ** Method  : MUSCL scheme plus centered scheme at ocean boundaries
[5770]66      !!              ld_msc_ups=T :
[3]67      !!
[12377]68      !! ** Action : - update pt(:,:,:,:,Krhs)  with the now advective tracer trends
[6140]69      !!             - send trends to trdtra module for further diagnostcs (l_trdtra=T)
[12377]70      !!             - poleward advective heat and salt transport (ln_diaptr=T)
[3]71      !!
[503]72      !! References : Estubier, A., and M. Levy, Notes Techn. Pole de Modelisation
73      !!              IPSL, Sept. 2000 (http://www.lodyc.jussieu.fr/opa)
74      !!----------------------------------------------------------------------
[12377]75      INTEGER                                  , INTENT(in   ) ::   kt              ! ocean time-step index
76      INTEGER                                  , INTENT(in   ) ::   Kbb, Kmm, Krhs  ! ocean time level indices
77      INTEGER                                  , INTENT(in   ) ::   kit000          ! first time step index
78      CHARACTER(len=3)                         , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator)
79      INTEGER                                  , INTENT(in   ) ::   kjpt            ! number of tracers
80      LOGICAL                                  , INTENT(in   ) ::   ld_msc_ups      ! use upstream scheme within muscl
81      REAL(wp)                                 , INTENT(in   ) ::   p2dt            ! tracer time-step
82      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pU, pV, pW      ! 3 ocean volume flux components
83      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! tracers and RHS of tracer equation
[2715]84      !
[9019]85      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
86      INTEGER  ::   ierr             ! local integer
87      REAL(wp) ::   zu, z0u, zzwx, zw , zalpha   ! local scalars
88      REAL(wp) ::   zv, z0v, zzwy, z0w           !   -      -
89      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwx, zslpx   ! 3D workspace
90      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwy, zslpy   ! -      -
[3]91      !!----------------------------------------------------------------------
[3294]92      !
93      IF( kt == kit000 )  THEN
[2528]94         IF(lwp) WRITE(numout,*)
95         IF(lwp) WRITE(numout,*) 'tra_adv : MUSCL advection scheme on ', cdtype
[3718]96         IF(lwp) WRITE(numout,*) '        : mixed up-stream           ', ld_msc_ups
[2528]97         IF(lwp) WRITE(numout,*) '~~~~~~~'
[3680]98         IF(lwp) WRITE(numout,*)
[2528]99         !
[5770]100         ! Upstream / MUSCL scheme indicator
[3680]101         !
[5770]102         ALLOCATE( xind(jpi,jpj,jpk), STAT=ierr )
[7753]103         xind(:,:,:) = 1._wp              ! set equal to 1 where up-stream is not needed
[5770]104         !
105         IF( ld_msc_ups ) THEN            ! define the upstream indicator (if asked)
106            ALLOCATE( upsmsk(jpi,jpj), STAT=ierr )
[7753]107            upsmsk(:,:) = 0._wp                             ! not upstream by default
[5770]108            !
[4990]109            DO jk = 1, jpkm1
[7753]110               xind(:,:,jk) = 1._wp                              &                 ! =>1 where up-stream is not needed
111                  &         - MAX ( rnfmsk(:,:) * rnfmsk_z(jk),  &                 ! =>0 near runoff mouths (& closed sea outflows)
112                  &                 upsmsk(:,:)                ) * tmask(:,:,jk)   ! =>0 in some user defined area
[3718]113            END DO
[3680]114         ENDIF 
[3718]115         !
116      ENDIF 
[4990]117      !     
[7646]118      l_trd = .FALSE.
119      l_hst = .FALSE.
120      l_ptr = .FALSE.
121      IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )      l_trd = .TRUE.
[12377]122      IF(   cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) )  )   l_ptr = .TRUE. 
[7646]123      IF(   cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. &
124         &                          iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) ) l_hst = .TRUE.
125      !
[6140]126      DO jn = 1, kjpt            !==  loop over the tracers  ==!
127         !
128         !                          !* Horizontal advective fluxes
129         !
130         !                                !-- first guess of the slopes
[7753]131         zwx(:,:,jpk) = 0._wp                   ! bottom values
132         zwy(:,:,jpk) = 0._wp 
[12377]133         DO_3D_10_10( 1, jpkm1 )
134            zwx(ji,jj,jk) = umask(ji,jj,jk) * ( pt(ji+1,jj,jk,jn,Kbb) - pt(ji,jj,jk,jn,Kbb) )
135            zwy(ji,jj,jk) = vmask(ji,jj,jk) * ( pt(ji,jj+1,jk,jn,Kbb) - pt(ji,jj,jk,jn,Kbb) )
136         END_3D
[9094]137         ! lateral boundary conditions   (changed sign)
[10425]138         CALL lbc_lnk_multi( 'traadv_mus', zwx, 'U', -1. , zwy, 'V', -1. )
[6140]139         !                                !-- Slopes of tracer
[7753]140         zslpx(:,:,jpk) = 0._wp                 ! bottom values
141         zslpy(:,:,jpk) = 0._wp
[12377]142         DO_3D_01_01( 1, jpkm1 )
143            zslpx(ji,jj,jk) =                    ( zwx(ji,jj,jk) + zwx(ji-1,jj  ,jk) )   &
144               &            * ( 0.25 + SIGN( 0.25, zwx(ji,jj,jk) * zwx(ji-1,jj  ,jk) ) )
145            zslpy(ji,jj,jk) =                    ( zwy(ji,jj,jk) + zwy(ji  ,jj-1,jk) )   &
146               &            * ( 0.25 + SIGN( 0.25, zwy(ji,jj,jk) * zwy(ji  ,jj-1,jk) ) )
147         END_3D
[503]148         !
[12377]149         DO_3D_01_01( 1, jpkm1 )
150            zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) ) * MIN(    ABS( zslpx(ji  ,jj,jk) ),   &
151               &                                                 2.*ABS( zwx  (ji-1,jj,jk) ),   &
152               &                                                 2.*ABS( zwx  (ji  ,jj,jk) ) )
153            zslpy(ji,jj,jk) = SIGN( 1., zslpy(ji,jj,jk) ) * MIN(    ABS( zslpy(ji,jj  ,jk) ),   &
154               &                                                 2.*ABS( zwy  (ji,jj-1,jk) ),   &
155               &                                                 2.*ABS( zwy  (ji,jj  ,jk) ) )
156         END_3D
[5770]157         !
[12377]158         DO_3D_00_00( 1, jpkm1 )
159            ! MUSCL fluxes
160            z0u = SIGN( 0.5, pU(ji,jj,jk) )
161            zalpha = 0.5 - z0u
162            zu  = z0u - 0.5 * pU(ji,jj,jk) * p2dt * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm)
163            zzwx = pt(ji+1,jj,jk,jn,Kbb) + xind(ji,jj,jk) * zu * zslpx(ji+1,jj,jk)
164            zzwy = pt(ji  ,jj,jk,jn,Kbb) + xind(ji,jj,jk) * zu * zslpx(ji  ,jj,jk)
165            zwx(ji,jj,jk) = pU(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy )
166            !
167            z0v = SIGN( 0.5, pV(ji,jj,jk) )
168            zalpha = 0.5 - z0v
169            zv  = z0v - 0.5 * pV(ji,jj,jk) * p2dt * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm)
170            zzwx = pt(ji,jj+1,jk,jn,Kbb) + xind(ji,jj,jk) * zv * zslpy(ji,jj+1,jk)
171            zzwy = pt(ji,jj  ,jk,jn,Kbb) + xind(ji,jj,jk) * zv * zslpy(ji,jj  ,jk)
172            zwy(ji,jj,jk) = pV(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy )
173         END_3D
[10425]174         CALL lbc_lnk_multi( 'traadv_mus', zwx, 'U', -1. , zwy, 'V', -1. )   ! lateral boundary conditions   (changed sign)
[503]175         !
[12377]176         DO_3D_00_00( 1, jpkm1 )
177            pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )       &
178            &                                     + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  ) )     &
179            &                                   * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
180         END_3D
[6140]181         !                                ! trend diagnostics
[7646]182         IF( l_trd )  THEN
[12377]183            CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_xad, zwx, pU, pt(:,:,:,jn,Kbb) )
184            CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_yad, zwy, pV, pt(:,:,:,jn,Kbb) )
[2528]185         END IF
[7646]186         !                                 ! "Poleward" heat and salt transports
187         IF( l_ptr )  CALL dia_ptr_hst( jn, 'adv', zwy(:,:,:) )
188         !                                 !  heat transport
189         IF( l_hst )  CALL dia_ar5_hst( jn, 'adv', zwx(:,:,:), zwy(:,:,:) )
[6140]190         !
191         !                          !* Vertical advective fluxes
192         !
[5770]193         !                                !-- first guess of the slopes
[7753]194         zwx(:,:, 1 ) = 0._wp                   ! surface & bottom boundary conditions
195         zwx(:,:,jpk) = 0._wp
[6140]196         DO jk = 2, jpkm1                       ! interior values
[12377]197            zwx(:,:,jk) = tmask(:,:,jk) * ( pt(:,:,jk-1,jn,Kbb) - pt(:,:,jk,jn,Kbb) )
[3]198         END DO
[5770]199         !                                !-- Slopes of tracer
[7753]200         zslpx(:,:,1) = 0._wp                   ! surface values
[12377]201         DO_3D_11_11( 2, jpkm1 )
202            zslpx(ji,jj,jk) =                     ( zwx(ji,jj,jk) + zwx(ji,jj,jk+1) )  &
203               &            * (  0.25 + SIGN( 0.25, zwx(ji,jj,jk) * zwx(ji,jj,jk+1) )  )
204         END_3D
205         DO_3D_11_11( 2, jpkm1 )
206            zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) ) * MIN(    ABS( zslpx(ji,jj,jk  ) ),   &
207               &                                                 2.*ABS( zwx  (ji,jj,jk+1) ),   &
208               &                                                 2.*ABS( zwx  (ji,jj,jk  ) )  )
209         END_3D
210         DO_3D_00_00( 1, jpk-2 )
211            z0w = SIGN( 0.5, pW(ji,jj,jk+1) )
212            zalpha = 0.5 + z0w
213            zw  = z0w - 0.5 * pW(ji,jj,jk+1) * p2dt * r1_e1e2t(ji,jj) / e3w(ji,jj,jk+1,Kmm)
214            zzwx = pt(ji,jj,jk+1,jn,Kbb) + xind(ji,jj,jk) * zw * zslpx(ji,jj,jk+1)
215            zzwy = pt(ji,jj,jk  ,jn,Kbb) + xind(ji,jj,jk) * zw * zslpx(ji,jj,jk  )
216            zwx(ji,jj,jk+1) = pW(ji,jj,jk+1) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) * wmask(ji,jj,jk)
217         END_3D
[6140]218         IF( ln_linssh ) THEN                   ! top values, linear free surface only
219            IF( ln_isfcav ) THEN                      ! ice-shelf cavities (top of the ocean)
[12377]220               DO_2D_11_11
221                  zwx(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kbb)
222               END_2D
[6140]223            ELSE                                      ! no cavities: only at the ocean surface
[12377]224               zwx(:,:,1) = pW(:,:,1) * pt(:,:,1,jn,Kbb)
[5770]225            ENDIF
226         ENDIF
227         !
[12377]228         DO_3D_00_00( 1, jpkm1 )
229            pt(ji,jj,jk,jn,Krhs) =  pt(ji,jj,jk,jn,Krhs) - ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) ) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
230         END_3D
[6140]231         !                                ! send trends for diagnostic
[12377]232         IF( l_trd )  CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_zad, zwx, pW, pt(:,:,:,jn,Kbb) )
[503]233         !
[6140]234      END DO                     ! end of tracer loop
[503]235      !
[5770]236   END SUBROUTINE tra_adv_mus
[3]237
238   !!======================================================================
[5770]239END MODULE traadv_mus
Note: See TracBrowser for help on using the repository browser.