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/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA – NEMO

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv_mus.F90 @ 10946

Last change on this file since 10946 was 10946, checked in by acc, 5 years ago

2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps : Convert STO, TRD and USR modules and all knock on effects of these conversions. Note change to USR module may have implications for the TEST CASES (not tested yet). Standard SETTE tested only
  • Property svn:keywords set to Id
File size: 15.1 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 "vectopt_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      !!             - htr_adv, str_adv : 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. ln_diaptr )                                               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 jk = 1, jpkm1                       ! interior values
134            DO jj = 1, jpjm1     
135               DO ji = 1, fs_jpim1   ! vector opt.
136                  zwx(ji,jj,jk) = umask(ji,jj,jk) * ( pt(ji+1,jj,jk,jn,Kbb) - pt(ji,jj,jk,jn,Kbb) )
137                  zwy(ji,jj,jk) = vmask(ji,jj,jk) * ( pt(ji,jj+1,jk,jn,Kbb) - pt(ji,jj,jk,jn,Kbb) )
138               END DO
139           END DO
140         END DO
141         ! lateral boundary conditions   (changed sign)
142         CALL lbc_lnk_multi( 'traadv_mus', zwx, 'U', -1. , zwy, 'V', -1. )
143         !                                !-- Slopes of tracer
144         zslpx(:,:,jpk) = 0._wp                 ! bottom values
145         zslpy(:,:,jpk) = 0._wp
146         DO jk = 1, jpkm1                       ! interior values
147            DO jj = 2, jpj
148               DO ji = fs_2, jpi   ! vector opt.
149                  zslpx(ji,jj,jk) =                    ( zwx(ji,jj,jk) + zwx(ji-1,jj  ,jk) )   &
150                     &            * ( 0.25 + SIGN( 0.25, zwx(ji,jj,jk) * zwx(ji-1,jj  ,jk) ) )
151                  zslpy(ji,jj,jk) =                    ( zwy(ji,jj,jk) + zwy(ji  ,jj-1,jk) )   &
152                     &            * ( 0.25 + SIGN( 0.25, zwy(ji,jj,jk) * zwy(ji  ,jj-1,jk) ) )
153               END DO
154            END DO
155         END DO
156         !
157         DO jk = 1, jpkm1                 !-- Slopes limitation
158            DO jj = 2, jpj
159               DO ji = fs_2, jpi   ! vector opt.
160                  zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) ) * MIN(    ABS( zslpx(ji  ,jj,jk) ),   &
161                     &                                                 2.*ABS( zwx  (ji-1,jj,jk) ),   &
162                     &                                                 2.*ABS( zwx  (ji  ,jj,jk) ) )
163                  zslpy(ji,jj,jk) = SIGN( 1., zslpy(ji,jj,jk) ) * MIN(    ABS( zslpy(ji,jj  ,jk) ),   &
164                     &                                                 2.*ABS( zwy  (ji,jj-1,jk) ),   &
165                     &                                                 2.*ABS( zwy  (ji,jj  ,jk) ) )
166               END DO
167           END DO
168         END DO
169         !
170         DO jk = 1, jpkm1                 !-- MUSCL horizontal advective fluxes
171            DO jj = 2, jpjm1
172               DO ji = fs_2, fs_jpim1   ! vector opt.
173                  ! MUSCL fluxes
174                  z0u = SIGN( 0.5, pU(ji,jj,jk) )
175                  zalpha = 0.5 - z0u
176                  zu  = z0u - 0.5 * pU(ji,jj,jk) * p2dt * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm)
177                  zzwx = pt(ji+1,jj,jk,jn,Kbb) + xind(ji,jj,jk) * zu * zslpx(ji+1,jj,jk)
178                  zzwy = pt(ji  ,jj,jk,jn,Kbb) + xind(ji,jj,jk) * zu * zslpx(ji  ,jj,jk)
179                  zwx(ji,jj,jk) = pU(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy )
180                  !
181                  z0v = SIGN( 0.5, pV(ji,jj,jk) )
182                  zalpha = 0.5 - z0v
183                  zv  = z0v - 0.5 * pV(ji,jj,jk) * p2dt * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm)
184                  zzwx = pt(ji,jj+1,jk,jn,Kbb) + xind(ji,jj,jk) * zv * zslpy(ji,jj+1,jk)
185                  zzwy = pt(ji,jj  ,jk,jn,Kbb) + xind(ji,jj,jk) * zv * zslpy(ji,jj  ,jk)
186                  zwy(ji,jj,jk) = pV(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy )
187               END DO
188            END DO
189         END DO
190         CALL lbc_lnk_multi( 'traadv_mus', zwx, 'U', -1. , zwy, 'V', -1. )   ! lateral boundary conditions   (changed sign)
191         !
192         DO jk = 1, jpkm1                 !-- Tracer advective trend
193            DO jj = 2, jpjm1     
194               DO ji = fs_2, fs_jpim1   ! vector opt.
195                  pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )       &
196                  &                                     + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  ) )     &
197                  &                                   * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
198               END DO
199           END DO
200         END DO       
201         !                                ! trend diagnostics
202         IF( l_trd )  THEN
203            CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_xad, zwx, pU, pt(:,:,:,jn,Kbb) )
204            CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_yad, zwy, pV, pt(:,:,:,jn,Kbb) )
205         END IF
206         !                                 ! "Poleward" heat and salt transports
207         IF( l_ptr )  CALL dia_ptr_hst( jn, 'adv', zwy(:,:,:) )
208         !                                 !  heat transport
209         IF( l_hst )  CALL dia_ar5_hst( jn, 'adv', zwx(:,:,:), zwy(:,:,:) )
210         !
211         !                          !* Vertical advective fluxes
212         !
213         !                                !-- first guess of the slopes
214         zwx(:,:, 1 ) = 0._wp                   ! surface & bottom boundary conditions
215         zwx(:,:,jpk) = 0._wp
216         DO jk = 2, jpkm1                       ! interior values
217            zwx(:,:,jk) = tmask(:,:,jk) * ( pt(:,:,jk-1,jn,Kbb) - pt(:,:,jk,jn,Kbb) )
218         END DO
219         !                                !-- Slopes of tracer
220         zslpx(:,:,1) = 0._wp                   ! surface values
221         DO jk = 2, jpkm1                       ! interior value
222            DO jj = 1, jpj
223               DO ji = 1, jpi
224                  zslpx(ji,jj,jk) =                     ( zwx(ji,jj,jk) + zwx(ji,jj,jk+1) )  &
225                     &            * (  0.25 + SIGN( 0.25, zwx(ji,jj,jk) * zwx(ji,jj,jk+1) )  )
226               END DO
227            END DO
228         END DO
229         DO jk = 2, jpkm1                 !-- Slopes limitation
230            DO jj = 1, jpj                      ! interior values
231               DO ji = 1, jpi
232                  zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) ) * MIN(    ABS( zslpx(ji,jj,jk  ) ),   &
233                     &                                                 2.*ABS( zwx  (ji,jj,jk+1) ),   &
234                     &                                                 2.*ABS( zwx  (ji,jj,jk  ) )  )
235               END DO
236            END DO
237         END DO
238         DO jk = 1, jpk-2                 !-- vertical advective flux
239            DO jj = 2, jpjm1     
240               DO ji = fs_2, fs_jpim1   ! vector opt.
241                  z0w = SIGN( 0.5, pW(ji,jj,jk+1) )
242                  zalpha = 0.5 + z0w
243                  zw  = z0w - 0.5 * pW(ji,jj,jk+1) * p2dt * r1_e1e2t(ji,jj) / e3w(ji,jj,jk+1,Kmm)
244                  zzwx = pt(ji,jj,jk+1,jn,Kbb) + xind(ji,jj,jk) * zw * zslpx(ji,jj,jk+1)
245                  zzwy = pt(ji,jj,jk  ,jn,Kbb) + xind(ji,jj,jk) * zw * zslpx(ji,jj,jk  )
246                  zwx(ji,jj,jk+1) = pW(ji,jj,jk+1) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) * wmask(ji,jj,jk)
247               END DO
248            END DO
249         END DO
250         IF( ln_linssh ) THEN                   ! top values, linear free surface only
251            IF( ln_isfcav ) THEN                      ! ice-shelf cavities (top of the ocean)
252               DO jj = 1, jpj
253                  DO ji = 1, jpi
254                     zwx(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kbb)
255                  END DO
256               END DO   
257            ELSE                                      ! no cavities: only at the ocean surface
258               zwx(:,:,1) = pW(:,:,1) * pt(:,:,1,jn,Kbb)
259            ENDIF
260         ENDIF
261         !
262         DO jk = 1, jpkm1                 !-- vertical advective trend
263            DO jj = 2, jpjm1     
264               DO ji = fs_2, fs_jpim1   ! vector opt.
265                  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)
266               END DO
267            END DO
268         END DO
269         !                                ! send trends for diagnostic
270         IF( l_trd )  CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_zad, zwx, pW, pt(:,:,:,jn,Kbb) )
271         !
272      END DO                     ! end of tracer loop
273      !
274   END SUBROUTINE tra_adv_mus
275
276   !!======================================================================
277END MODULE traadv_mus
Note: See TracBrowser for help on using the repository browser.