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

source: NEMO/trunk/src/OCE/TRA/traadv_mus_lf.F90 @ 14072

Last change on this file since 14072 was 14072, checked in by laurent, 3 years ago

Merging branch "2020/dev_r13648_ASINTER-04_laurent_bulk_ice", ticket #2369
File size: 14.8 KB
Line 
1MODULE traadv_mus_lf
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 lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
32
33   IMPLICIT NONE
34   PRIVATE
35
36   PUBLIC   tra_adv_mus_lf   ! 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 1 configurations)
40   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   xind     !: mixed upstream/centered index
41   
42   LOGICAL  ::   l_trd   ! flag to compute trends
43   LOGICAL  ::   l_ptr   ! flag to compute poleward transport
44   LOGICAL  ::   l_hst   ! flag to compute heat/salt transport
45
46   !! * Substitutions
47#  include "do_loop_substitute.h90"
48#  include "domzgr_substitute.h90"
49
50#define initial_slop_i(out, ji) \
51        out = umask(ji,jj,jk) * ( pt(ji+1,jj,jk,jn,Kbb) - pt(ji,jj,jk,jn,Kbb) )
52
53#define initial_slop_j(out, jj) \
54        out = vmask(ji,jj,jk) * ( pt(ji,jj+1,jk,jn,Kbb) - pt(ji,jj,jk,jn,Kbb) )
55
56#define initial_slop_k(out, jk) \
57        out = tmask(ji,jj,jk) * ( pt(ji,jj,jk-1,jn,Kbb) - pt(ji,jj,jk,jn,Kbb) )
58
59#define tracer_slop(out, zzw, zzwm1) \
60        out = ( zzw + zzwm1 ) * ( 0.25 + SIGN( 0.25_wp, zzw * zzwm1 ) )
61
62#define limitation_slop(out, zzslp, zzwm1, zzw) \
63        out = SIGN( 1.0_wp, zzslp ) * MIN( ABS( zzslp ), 2.*ABS( zzwm1 ), 2.*ABS( zzw ) )
64
65#define vertical_adv_flux_i(out, jk, slp, slp1) \
66        z0u = SIGN( 0.5_wp, pU(ji,jj,jk) ) ; \
67        zalpha = 0.5 - z0u ; \
68        zu  = z0u - 0.5 * pU(ji,jj,jk) * p2dt * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm) ; \
69        zzwx = pt(ji+1,jj,jk,jn,Kbb) + xind(ji,jj,jk) * zu * slp1 ; \
70        zzwy = pt(ji  ,jj,jk,jn,Kbb) + xind(ji,jj,jk) * zu * slp ; \
71        out = pU(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy )
72
73#define vertical_adv_flux_j(out, jk, slp, slp1) \
74        z0v = SIGN( 0.5_wp, pV(ji,jj,jk) ) ; \
75        zalpha = 0.5 - z0v ; \
76        zv  = z0v - 0.5 * pV(ji,jj,jk) * p2dt * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm) ; \
77        zzwx = pt(ji,jj+1,jk,jn,Kbb) + xind(ji,jj,jk) * zv * slp1 ; \
78        zzwy = pt(ji,jj  ,jk,jn,Kbb) + xind(ji,jj,jk) * zv * slp ; \
79        out = pV(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy )
80
81#define vertical_adv_flux(out, jk, slp, slp1) \
82        z0w = SIGN( 0.5_wp, pW(ji,jj,jk+1) ) ; \
83        zalpha = 0.5 + z0w ; \
84        zw  = z0w - 0.5 * pW(ji,jj,jk+1) * p2dt * r1_e1e2t(ji,jj) / e3w(ji,jj,jk+1,Kmm) ; \
85        zzwx = pt(ji,jj,jk+1,jn,Kbb) + xind(ji,jj,jk) * zw * slp1 ; \
86        zzwy = pt(ji,jj,jk  ,jn,Kbb) + xind(ji,jj,jk) * zw * slp ; \
87        out = pW(ji,jj,jk+1) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) * wmask(ji,jj,jk) 
88
89   !!----------------------------------------------------------------------
90   !!----------------------------------------------------------------------
91   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
92   !! $Id: traadv_mus.F90 13619 2020-10-16 08:41:21Z francesca $
93   !! Software governed by the CeCILL license (see ./LICENSE)
94   !!----------------------------------------------------------------------
95CONTAINS
96
97   SUBROUTINE tra_adv_mus_lf( kt, kit000, cdtype, p2dt, pU, pV, pW,             &
98      &                    Kbb, Kmm, pt, kjpt, Krhs, ld_msc_ups )
99      !!----------------------------------------------------------------------
100      !!                    ***  ROUTINE tra_adv_mus  ***
101      !!
102      !! ** Purpose :   Compute the now trend due to total advection of tracers
103      !!              using a MUSCL scheme (Monotone Upstream-centered Scheme for
104      !!              Conservation Laws) and add it to the general tracer trend.
105      !!
106      !! ** Method  : MUSCL scheme plus centered scheme at ocean boundaries
107      !!              ld_msc_ups=T :
108      !!
109      !! ** Action : - update pt(:,:,:,:,Krhs)  with the now advective tracer trends
110      !!             - send trends to trdtra module for further diagnostcs (l_trdtra=T)
111      !!             - poleward advective heat and salt transport (ln_diaptr=T)
112      !!
113      !! References : Estubier, A., and M. Levy, Notes Techn. Pole de Modelisation
114      !!              IPSL, Sept. 2000 (http://www.lodyc.jussieu.fr/opa)
115      !!----------------------------------------------------------------------
116      INTEGER                                  , INTENT(in   ) ::   kt              ! ocean time-step index
117      INTEGER                                  , INTENT(in   ) ::   Kbb, Kmm, Krhs  ! ocean time level indices
118      INTEGER                                  , INTENT(in   ) ::   kit000          ! first time step index
119      CHARACTER(len=3)                         , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator)
120      INTEGER                                  , INTENT(in   ) ::   kjpt            ! number of tracers
121      LOGICAL                                  , INTENT(in   ) ::   ld_msc_ups      ! use upstream scheme within muscl
122      REAL(wp)                                 , INTENT(in   ) ::   p2dt            ! tracer time-step
123      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pU, pV, pW      ! 3 ocean volume flux components
124      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! tracers and RHS of tracer equation
125      !
126      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
127      INTEGER  ::   ierr             ! local integer
128      REAL(wp) ::   zu, z0u, zw , zalpha   ! local scalars
129      REAL(wp) ::   zv, z0v, z0w           !   -      -
130      REAL(wp) ::   zzwx, zzwxm1, zzwxp1, zzwy, zzwym1, zzwyp1
131      REAL(wp) ::   zzslpx, zzslpxp1, zzslpy, zzslpyp1
132      REAL(wp), TARGET, DIMENSION(jpi,jpj) :: zzwz_buf, zzwzp1_buf, zzwzp2_buf
133      REAL(wp), TARGET, DIMENSION(jpi,jpj) :: zzslpz_buf, zzslpzp1_buf
134      REAL(wp), POINTER, DIMENSION(:,:)    :: tmp, zzwz_ptr, zzwzp1_ptr, zzwzp2_ptr
135      REAL(wp), POINTER, DIMENSION(:,:)    :: zzslpz_ptr, zzslpzp1_ptr 
136      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwz, zwx, zwy   ! 3D workspace
137      !!----------------------------------------------------------------------
138      !
139      IF( kt == kit000 )  THEN
140         IF(lwp) WRITE(numout,*)
141         IF(lwp) WRITE(numout,*) 'tra_adv : MUSCL advection scheme on ', cdtype
142         IF(lwp) WRITE(numout,*) '        : mixed up-stream           ', ld_msc_ups
143         IF(lwp) WRITE(numout,*) '~~~~~~~'
144         IF(lwp) WRITE(numout,*)
145         !
146         ! Upstream / MUSCL scheme indicator
147         !
148         ALLOCATE( xind(jpi,jpj,jpk), STAT=ierr )
149         xind(:,:,:) = 1._wp              ! set equal to 1 where up-stream is not needed
150         !
151         IF( ld_msc_ups ) THEN            ! define the upstream indicator (if asked)
152            ALLOCATE( upsmsk(jpi,jpj), STAT=ierr )
153            upsmsk(:,:) = 0._wp                             ! not upstream by default
154            !
155            DO jk = 1, jpkm1
156               xind(:,:,jk) = 1._wp                              &                 ! =>1 where up-stream is not needed
157                  &         - MAX ( rnfmsk(:,:) * rnfmsk_z(jk),  &                 ! =>0 near runoff mouths (& closed sea outflows)
158                  &                 upsmsk(:,:)                ) * tmask(:,:,jk)   ! =>0 in some user defined area
159            END DO
160         ENDIF 
161         !
162      ENDIF 
163      !     
164      l_trd = .FALSE.
165      l_hst = .FALSE.
166      l_ptr = .FALSE.
167      IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )      l_trd = .TRUE.
168      IF(   cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) )  )   l_ptr = .TRUE. 
169      IF(   cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. &
170         &                          iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) ) l_hst = .TRUE.
171      !
172         zzwz_ptr => zzwz_buf
173         zzwzp1_ptr => zzwzp1_buf
174         zzwzp2_ptr => zzwzp2_buf
175         zzslpz_ptr => zzslpz_buf
176         zzslpzp1_ptr => zzslpzp1_buf
177         !
178      DO jn = 1, kjpt            !==  loop over the tracers  ==!
179         !
180         zwx(:,:,jpk) = 0._wp                   ! bottom values
181         zwy(:,:,jpk) = 0._wp
182         zwz(:,:, 1 ) = 0._wp                   ! surface & bottom boundary conditions
183         zwz(:,:,jpk) = 0._wp
184         !                          !* Horizontal advective fluxes
185         !
186         !!----------------------------------------------------------------------
187         DO_3D( 1, 0, 1, 0, 1, jpkm1 )
188            !-- first guess of the slopes
189            initial_slop_i(zzwxm1, ji-1)
190            initial_slop_i(zzwx, ji)
191            initial_slop_i(zzwxp1, ji+1)
192
193            initial_slop_j(zzwym1, jj-1)
194            initial_slop_j(zzwy, jj)
195            initial_slop_j(zzwyp1, jj+1)
196            !-- Slopes of tracer
197            tracer_slop(zzslpx, zzwx, zzwxm1) 
198            tracer_slop(zzslpxp1, zzwxp1, zzwx) 
199            tracer_slop(zzslpy, zzwy, zzwym1) 
200            tracer_slop(zzslpyp1, zzwyp1, zzwy) 
201            !-- Slopes limitation
202            limitation_slop(zzslpx, zzslpx, zzwxm1, zzwx)
203            limitation_slop(zzslpxp1, zzslpxp1, zzwx, zzwxp1)
204            limitation_slop(zzslpy, zzslpy, zzwym1, zzwy)
205            limitation_slop(zzslpyp1, zzslpyp1, zzwy, zzwyp1)
206            !-- MUSCL horizontal advective fluxes
207            vertical_adv_flux_i(zwx(ji,jj,jk), jk, zzslpx, zzslpxp1) 
208            vertical_adv_flux_j(zwy(ji,jj,jk), jk, zzslpy, zzslpyp1) 
209         END_3D
210         !
211         DO_3D( 0, 0, 0, 0, 1, jpkm1 )    !-- Tracer advective trend
212            pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( zwx(ji,jj,jk) - zwx(ji-1,jj,jk)       &
213            &                                     + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  ) )     &
214            &                                     * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
215         END_3D
216         !                          !* Vertical advective fluxes
217         DO_2D( 0, 0, 0, 0 ) 
218            !-- first guess of the slopes
219            initial_slop_k(zzwz_ptr(ji,jj), 2)
220            initial_slop_k(zzwzp1_ptr(ji,jj), 3)
221            !-- Slopes of tracer
222            tracer_slop(zzslpz_ptr(ji,jj), zzwz_ptr(ji,jj), zzwzp1_ptr(ji,jj))
223            !-- Slopes limitation
224            limitation_slop(zzslpz_ptr(ji,jj), zzslpz_ptr(ji,jj), zzwzp1_ptr(ji,jj), zzwz_ptr(ji,jj))
225            !-- vertical advective flux
226            vertical_adv_flux(zwz(ji,jj,2), 1, 0, zzslpz_ptr(ji,jj))
227         END_2D
228         
229         DO jk = 2, jpk-3   
230            DO_2D( 0, 0, 0, 0 )
231               !-- first guess of the slopes
232               initial_slop_k(zzwzp2_ptr(ji,jj), jk+2)
233               !-- Slopes of tracer
234               tracer_slop(zzslpzp1_ptr(ji,jj), zzwzp1_ptr(ji,jj), zzwzp2_ptr(ji,jj))
235               !-- Slopes limitation
236               limitation_slop(zzslpzp1_ptr(ji,jj), zzslpzp1_ptr(ji,jj), zzwzp2_ptr(ji,jj), zzwzp1_ptr(ji,jj))
237               !-- vertical advective flux
238               vertical_adv_flux(zwz(ji,jj,jk+1), jk, zzslpz_ptr(ji,jj), zzslpzp1_ptr(ji,jj))
239            END_2D
240            tmp => zzwzp1_ptr
241            zzwzp1_ptr => zzwzp2_ptr
242            zzwzp2_ptr => tmp
243
244            tmp => zzslpz_ptr
245            zzslpz_ptr => zzslpzp1_ptr
246            zzslpzp1_ptr => tmp 
247         END DO
248         DO_2D( 0, 0, 0, 0 ) 
249            !-- Slopes of tracer
250            tracer_slop(zzslpzp1_ptr(ji,jj), zzwzp1_ptr(ji,jj), 0)
251            !-- Slopes limitation
252            limitation_slop(zzslpzp1_ptr(ji,jj), zzslpzp1_ptr(ji,jj), 0, zzwzp1_ptr(ji,jj))
253            !-- vertical advective flux
254            vertical_adv_flux(zwz(ji,jj,jpk-1), jpk-2, zzslpz_ptr(ji,jj), zzslpzp1_ptr(ji,jj))
255         END_2D
256
257         IF( ln_linssh ) THEN                   ! top values, linear free surface only
258            IF( ln_isfcav ) THEN                      ! ice-shelf cavities (top of the ocean)
259               DO_2D( 1, 1, 1, 1 )
260                  zwz(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kbb)
261               END_2D
262            ELSE                                      ! no cavities: only at the ocean surface
263               zwz(:,:,1) = pW(:,:,1) * pt(:,:,1,jn,Kbb)
264            ENDIF
265         ENDIF
266         !
267         DO_3D( 0, 0, 0, 0, 1, jpkm1 )     !-- vertical advective trend
268            pt(ji,jj,jk,jn,Krhs) =  pt(ji,jj,jk,jn,Krhs) - ( zwz(ji,jj,jk) - zwz(ji,jj,jk+1) )   &
269               &                                      * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
270         END_3D
271         !                                ! trend horizontal diagnostics
272         IF( l_trd )  THEN
273            CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_xad, zwx, pU, pt(:,:,:,jn,Kbb) )
274            CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_yad, zwy, pV, pt(:,:,:,jn,Kbb) )
275         END IF
276         !                                 ! "Poleward" heat and salt transports
277         IF( l_ptr )  CALL dia_ptr_hst( jn, 'adv', zwy(:,:,:) )
278         !                                 !  heat transport
279         IF( l_hst )  CALL dia_ar5_hst( jn, 'adv', zwx(:,:,:), zwy(:,:,:) )
280         !
281         !                                ! send vertical trends for diagnostic
282         IF( l_trd )  CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_zad, zwz, pW, pt(:,:,:,jn,Kbb) )
283         !
284      END DO                     ! end of tracer loop
285      !
286   END SUBROUTINE tra_adv_mus_lf
287
288   !!======================================================================
289END MODULE traadv_mus_lf
Note: See TracBrowser for help on using the repository browser.