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/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/OCE/TRA – NEMO

source: NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/OCE/TRA/traadv_mus.F90 @ 13900

Last change on this file since 13900 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
Line 
1MODULE traadv_mus
2   !!======================================================================
3   !!                       ***  MODULE  traadv_mus  ***
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   !!            3.7  !  2015-09  (G. Madec) add the ice-shelf cavities boundary condition
12   !!----------------------------------------------------------------------
13
14   !!----------------------------------------------------------------------
15   !!   tra_adv_mus   : update the tracer trend with the horizontal
16   !!                   and vertical advection trends using MUSCL scheme
17   !!----------------------------------------------------------------------
18   USE oce            ! ocean dynamics and active tracers
19   USE trc_oce        ! share passive tracers/Ocean variables
20   USE dom_oce        ! ocean space and time domain
21   USE trd_oce        ! trends: ocean variables
22   USE trdtra         ! tracers trends manager
23   USE sbcrnf         ! river runoffs
24   USE diaptr         ! poleward transport diagnostics
25   USE diaar5         ! AR5 diagnostics
26
27   !
28   USE iom            ! XIOS library
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   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
33
34   IMPLICIT NONE
35   PRIVATE
36
37   PUBLIC   tra_adv_mus   ! routine called by traadv.F90
38   
39   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   upsmsk   !: mixed upstream/centered scheme near some straits
40   !                                                           !  and in closed seas (orca 2 and 1 configurations)
41   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   xind     !: mixed upstream/centered index
42   
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
47   !! * Substitutions
48#  include "do_loop_substitute.h90"
49   !!----------------------------------------------------------------------
50   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
51   !! $Id$
52   !! Software governed by the CeCILL license (see ./LICENSE)
53   !!----------------------------------------------------------------------
54CONTAINS
55
56   SUBROUTINE tra_adv_mus( kt, kit000, cdtype, p2dt, pU, pV, pW,             &
57      &                    Kbb, Kmm, pt, kjpt, Krhs, ld_msc_ups )
58      !!----------------------------------------------------------------------
59      !!                    ***  ROUTINE tra_adv_mus  ***
60      !!
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.
64      !!
65      !! ** Method  : MUSCL scheme plus centered scheme at ocean boundaries
66      !!              ld_msc_ups=T :
67      !!
68      !! ** Action : - update pt(:,:,:,:,Krhs)  with the now advective tracer trends
69      !!             - send trends to trdtra module for further diagnostcs (l_trdtra=T)
70      !!             - poleward advective heat and salt transport (ln_diaptr=T)
71      !!
72      !! References : Estubier, A., and M. Levy, Notes Techn. Pole de Modelisation
73      !!              IPSL, Sept. 2000 (http://www.lodyc.jussieu.fr/opa)
74      !!----------------------------------------------------------------------
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
84      !
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   ! -      -
91      !!----------------------------------------------------------------------
92      !
93      IF( kt == kit000 )  THEN
94         IF(lwp) WRITE(numout,*)
95         IF(lwp) WRITE(numout,*) 'tra_adv : MUSCL advection scheme on ', cdtype
96         IF(lwp) WRITE(numout,*) '        : mixed up-stream           ', ld_msc_ups
97         IF(lwp) WRITE(numout,*) '~~~~~~~'
98         IF(lwp) WRITE(numout,*)
99         !
100         ! Upstream / MUSCL scheme indicator
101         !
102         ALLOCATE( xind(jpi,jpj,jpk), STAT=ierr )
103         xind(:,:,:) = 1._wp              ! set equal to 1 where up-stream is not needed
104         !
105         IF( ld_msc_ups ) THEN            ! define the upstream indicator (if asked)
106            ALLOCATE( upsmsk(jpi,jpj), STAT=ierr )
107            upsmsk(:,:) = 0._wp                             ! not upstream by default
108            !
109            DO jk = 1, jpkm1
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
113            END DO
114         ENDIF 
115         !
116      ENDIF 
117      !     
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.
122      IF(   cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) )  )   l_ptr = .TRUE. 
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      !
126      DO jn = 1, kjpt            !==  loop over the tracers  ==!
127         !
128         !                          !* Horizontal advective fluxes
129         !
130         !                                !-- first guess of the slopes
131         zwx(:,:,jpk) = 0._wp                   ! bottom values
132         zwy(:,:,jpk) = 0._wp 
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
137         ! lateral boundary conditions   (changed sign)
138         CALL lbc_lnk_multi( 'traadv_mus', zwx, 'U', -1. , zwy, 'V', -1. )
139         !                                !-- Slopes of tracer
140         zslpx(:,:,jpk) = 0._wp                 ! bottom values
141         zslpy(:,:,jpk) = 0._wp
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
148         !
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
157         !
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
174         CALL lbc_lnk_multi( 'traadv_mus', zwx, 'U', -1. , zwy, 'V', -1. )   ! lateral boundary conditions   (changed sign)
175         !
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
181         !                                ! trend diagnostics
182         IF( l_trd )  THEN
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) )
185         END IF
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(:,:,:) )
190         !
191         !                          !* Vertical advective fluxes
192         !
193         !                                !-- first guess of the slopes
194         zwx(:,:, 1 ) = 0._wp                   ! surface & bottom boundary conditions
195         zwx(:,:,jpk) = 0._wp
196         DO jk = 2, jpkm1                       ! interior values
197            zwx(:,:,jk) = tmask(:,:,jk) * ( pt(:,:,jk-1,jn,Kbb) - pt(:,:,jk,jn,Kbb) )
198         END DO
199         !                                !-- Slopes of tracer
200         zslpx(:,:,1) = 0._wp                   ! surface values
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
218         IF( ln_linssh ) THEN                   ! top values, linear free surface only
219            IF( ln_isfcav ) THEN                      ! ice-shelf cavities (top of the ocean)
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
223            ELSE                                      ! no cavities: only at the ocean surface
224               zwx(:,:,1) = pW(:,:,1) * pt(:,:,1,jn,Kbb)
225            ENDIF
226         ENDIF
227         !
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
231         !                                ! send trends for diagnostic
232         IF( l_trd )  CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_zad, zwx, pW, pt(:,:,:,jn,Kbb) )
233         !
234      END DO                     ! end of tracer loop
235      !
236   END SUBROUTINE tra_adv_mus
237
238   !!======================================================================
239END MODULE traadv_mus
Note: See TracBrowser for help on using the repository browser.