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 branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

source: branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_mus.F90 @ 5870

Last change on this file since 5870 was 5870, checked in by acc, 8 years ago

Branch 2015/dev_r5803_NOC_WAD. Merge in trunk changes from 5803 to 5869 in preparation for merge. Also tidied and reorganised some wetting and drying code. Renamed wadlmt.F90 to wetdry.F90. Wetting drying code changes restricted to domzgr.F90, domvvl.F90 nemogcm.F90 sshwzv.F90, dynspg_ts.F90, wetdry.F90 and dynhpg.F90. Code passes full SETTE tests with ln_wd=.false.. Still awaiting test case for checking with ln_wd=.false.
  • Property svn:keywords set to Id
File size: 15.7 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 dynspg_oce     ! choice/control of key cpp for surface pressure gradient
24   USE sbcrnf         ! river runoffs
25   USE diaptr         ! poleward transport diagnostics
26   !
27   USE wrk_nemo       ! Memory Allocation
28   USE timing         ! Timing
29   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
30   USE in_out_manager ! I/O manager
31   USE lib_mpp        ! distribued memory computing
32   USE lbclnk         ! ocean lateral boundary condition (or mpp link)
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 4 configurations)
41   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   xind     !: mixed upstream/centered index
42   
43   !! * Substitutions
44#  include "domzgr_substitute.h90"
45#  include "vectopt_loop_substitute.h90"
46   !!----------------------------------------------------------------------
47   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
48   !! $Id$
49   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
50   !!----------------------------------------------------------------------
51CONTAINS
52
53   SUBROUTINE tra_adv_mus( kt, kit000, cdtype, p2dt, pun, pvn, pwn,             &
54      &                                              ptb, pta, kjpt, ld_msc_ups )
55      !!----------------------------------------------------------------------
56      !!                    ***  ROUTINE tra_adv_mus  ***
57      !!
58      !! ** Purpose :   Compute the now trend due to total advection of tracers
59      !!              using a MUSCL scheme (Monotone Upstream-centered Scheme for
60      !!              Conservation Laws) and add it to the general tracer trend.
61      !!
62      !! ** Method  : MUSCL scheme plus centered scheme at ocean boundaries
63      !!              ld_msc_ups=T :
64      !!
65      !! ** Action  : - update (ta,sa) with the now advective tracer trends
66      !!              - send trends to trdtra module for further diagnostcs
67      !!
68      !! References : Estubier, A., and M. Levy, Notes Techn. Pole de Modelisation
69      !!              IPSL, Sept. 2000 (http://www.lodyc.jussieu.fr/opa)
70      !!----------------------------------------------------------------------
71      INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index
72      INTEGER                              , INTENT(in   ) ::   kit000          ! first time step index
73      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator)
74      INTEGER                              , INTENT(in   ) ::   kjpt            ! number of tracers
75      LOGICAL                              , INTENT(in   ) ::   ld_msc_ups      ! use upstream scheme within muscl
76      REAL(wp), DIMENSION(        jpk     ), INTENT(in   ) ::   p2dt            ! vertical profile of tracer time-step
77      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pun, pvn, pwn   ! 3 ocean velocity components
78      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb             ! before tracer field
79      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta             ! tracer trend
80      !
81      INTEGER  ::   ji, jj, jk, jn       ! dummy loop indices
82      INTEGER  ::   ierr                 ! local integer
83      REAL(wp) ::   zu, z0u, zzwx, zw    ! local scalars
84      REAL(wp) ::   zv, z0v, zzwy, z0w   !   -      -
85      REAL(wp) ::   zdt, zalpha          !   -      -
86      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zslpx, zslpy   ! 3D workspace
87      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zwx  , zwy     ! -      -
88      !!----------------------------------------------------------------------
89      !
90      IF( nn_timing == 1 )  CALL timing_start('tra_adv_mus')
91      !
92      CALL wrk_alloc( jpi,jpj,jpk,   zslpx, zslpy, zwx, zwy )
93      !
94      IF( kt == kit000 )  THEN
95         IF(lwp) WRITE(numout,*)
96         IF(lwp) WRITE(numout,*) 'tra_adv : MUSCL advection scheme on ', cdtype
97         IF(lwp) WRITE(numout,*) '        : mixed up-stream           ', ld_msc_ups
98         IF(lwp) WRITE(numout,*) '~~~~~~~'
99         IF(lwp) WRITE(numout,*)
100         !
101         ! Upstream / MUSCL scheme indicator
102         !
103         ALLOCATE( xind(jpi,jpj,jpk), STAT=ierr )
104         xind(:,:,:) = 1._wp              ! set equal to 1 where up-stream is not needed
105         !
106         IF( ld_msc_ups ) THEN            ! define the upstream indicator (if asked)
107            ALLOCATE( upsmsk(jpi,jpj), STAT=ierr )
108            upsmsk(:,:) = 0._wp                             ! not upstream by default
109            !
110            DO jk = 1, jpkm1
111               xind(:,:,jk) = 1._wp                              &                 ! =>1 where up-stream is not needed
112                  &         - MAX ( rnfmsk(:,:) * rnfmsk_z(jk),  &                 ! =>0 near runoff mouths (& closed sea outflows)
113                  &                 upsmsk(:,:)                ) * tmask(:,:,jk)   ! =>0 in some user defined area
114            END DO
115         ENDIF 
116         !
117      ENDIF 
118      !     
119      !                                                     ! ===========
120      DO jn = 1, kjpt                                       ! tracer loop
121         !                                                  ! ===========
122         ! I. Horizontal advective fluxes
123         ! ------------------------------
124         ! first guess of the slopes
125         zwx(:,:,jpk) = 0.e0   ;   zwy(:,:,jpk) = 0.e0        ! bottom values
126         ! interior values
127         DO jk = 1, jpkm1
128            DO jj = 1, jpjm1     
129               DO ji = 1, fs_jpim1   ! vector opt.
130                  zwx(ji,jj,jk) = umask(ji,jj,jk) * ( ptb(ji+1,jj,jk,jn) - ptb(ji,jj,jk,jn) )
131                  zwy(ji,jj,jk) = vmask(ji,jj,jk) * ( ptb(ji,jj+1,jk,jn) - ptb(ji,jj,jk,jn) )
132               END DO
133           END DO
134         END DO
135         !
136         CALL lbc_lnk( zwx, 'U', -1. )                        ! lateral boundary conditions on zwx, zwy   (changed sign)
137         CALL lbc_lnk( zwy, 'V', -1. )
138         !                                             !-- Slopes of tracer
139         zslpx(:,:,jpk) = 0.e0   ;   zslpy(:,:,jpk) = 0.e0    ! bottom values
140         DO jk = 1, jpkm1                                     ! interior values
141            DO jj = 2, jpj
142               DO ji = fs_2, jpi   ! vector opt.
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 DO
148            END DO
149         END DO
150         !
151         DO jk = 1, jpkm1                                     ! Slopes limitation
152            DO jj = 2, jpj
153               DO ji = fs_2, jpi   ! vector opt.
154                  zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) ) * MIN(    ABS( zslpx(ji  ,jj,jk) ),   &
155                     &                                                 2.*ABS( zwx  (ji-1,jj,jk) ),   &
156                     &                                                 2.*ABS( zwx  (ji  ,jj,jk) ) )
157                  zslpy(ji,jj,jk) = SIGN( 1., zslpy(ji,jj,jk) ) * MIN(    ABS( zslpy(ji,jj  ,jk) ),   &
158                     &                                                 2.*ABS( zwy  (ji,jj-1,jk) ),   &
159                     &                                                 2.*ABS( zwy  (ji,jj  ,jk) ) )
160               END DO
161           END DO
162         END DO
163         !
164         !                                             !-- MUSCL horizontal advective fluxes
165         DO jk = 1, jpkm1                                     ! interior values
166            zdt  = p2dt(jk)
167            DO jj = 2, jpjm1
168               DO ji = fs_2, fs_jpim1   ! vector opt.
169                  ! MUSCL fluxes
170                  z0u = SIGN( 0.5, pun(ji,jj,jk) )
171                  zalpha = 0.5 - z0u
172                  zu  = z0u - 0.5 * pun(ji,jj,jk) * zdt / ( e1e2u(ji,jj) * fse3u(ji,jj,jk) )
173                  zzwx = ptb(ji+1,jj,jk,jn) + xind(ji,jj,jk) * zu * zslpx(ji+1,jj,jk)
174                  zzwy = ptb(ji  ,jj,jk,jn) + xind(ji,jj,jk) * zu * zslpx(ji  ,jj,jk)
175                  zwx(ji,jj,jk) = pun(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy )
176                  !
177                  z0v = SIGN( 0.5, pvn(ji,jj,jk) )
178                  zalpha = 0.5 - z0v
179                  zv  = z0v - 0.5 * pvn(ji,jj,jk) * zdt / ( e1e2v(ji,jj) * fse3v(ji,jj,jk) )
180                  zzwx = ptb(ji,jj+1,jk,jn) + xind(ji,jj,jk) * zv * zslpy(ji,jj+1,jk)
181                  zzwy = ptb(ji,jj  ,jk,jn) + xind(ji,jj,jk) * zv * zslpy(ji,jj  ,jk)
182                  zwy(ji,jj,jk) = pvn(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy )
183               END DO
184            END DO
185         END DO
186         CALL lbc_lnk( zwx, 'U', -1. )   ;   CALL lbc_lnk( zwy, 'V', -1. )   ! lateral boundary conditions   (changed sign)
187         !
188         DO jk = 1, jpkm1         ! Tracer flux divergence at t-point added to the general trend
189            DO jj = 2, jpjm1     
190               DO ji = fs_2, fs_jpim1   ! vector opt.
191                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )       &
192                  &                                     + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  ) )     &
193                  &                                   / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) )
194               END DO
195           END DO
196         END DO       
197         !                                 ! trend diagnostics (contribution of upstream fluxes)
198         IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR.   &
199            &( cdtype == 'TRC' .AND. l_trdtrc )      )  THEN
200            CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptb(:,:,:,jn) )
201            CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptb(:,:,:,jn) )
202         END IF
203         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
204         IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN 
205            IF( jn == jp_tem )  htr_adv(:) = ptr_sj( zwy(:,:,:) )
206            IF( jn == jp_sal )  str_adv(:) = ptr_sj( zwy(:,:,:) )
207         ENDIF
208
209         ! II. Vertical advective fluxes
210         ! -----------------------------
211         !                                !-- first guess of the slopes
212         zwx(:,:, 1 ) = 0._wp                   ! surface & bottom boundary conditions
213         zwx(:,:,jpk) = 0._wp                   ! surface & bottom boundary conditions
214         DO jk = 2, jpkm1                                     ! interior values
215            zwx(:,:,jk) = tmask(:,:,jk) * ( ptb(:,:,jk-1,jn) - ptb(:,:,jk,jn) )
216         END DO
217
218         !                                !-- Slopes of tracer
219         zslpx(:,:,1) = 0._wp                   ! surface values
220         DO jk = 2, jpkm1                       ! interior value
221            DO jj = 1, jpj
222               DO ji = 1, jpi
223                  zslpx(ji,jj,jk) =                    ( zwx(ji,jj,jk) + zwx(ji,jj,jk+1) )   &
224                     &            * ( 0.25 + SIGN( 0.25, zwx(ji,jj,jk) * zwx(ji,jj,jk+1) ) )
225               END DO
226            END DO
227         END DO
228         !                                !-- Slopes limitation
229         DO jk = 2, jpkm1                       ! interior values
230            DO jj = 1, jpj
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         !                                !-- vertical advective flux
239         DO jk = 1, jpkm1                       ! interior values
240            zdt  = p2dt(jk)
241            DO jj = 2, jpjm1     
242               DO ji = fs_2, fs_jpim1   ! vector opt.
243                  z0w = SIGN( 0.5, pwn(ji,jj,jk+1) )
244                  zalpha = 0.5 + z0w
245                  zw  = z0w - 0.5 * pwn(ji,jj,jk+1) * zdt  / ( e1e2t(ji,jj) * fse3w(ji,jj,jk+1) )
246                  zzwx = ptb(ji,jj,jk+1,jn) + xind(ji,jj,jk) * zw * zslpx(ji,jj,jk+1)
247                  zzwy = ptb(ji,jj,jk  ,jn) + xind(ji,jj,jk) * zw * zslpx(ji,jj,jk  )
248                  zwx(ji,jj,jk+1) = pwn(ji,jj,jk+1) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) * wmask(ji,jj,jk)
249               END DO
250            END DO
251         END DO
252         !                                      ! top values  (bottom already set to zero)
253         IF( lk_vvl ) THEN                            !  variable volume
254            zwx(:,:, 1 ) = 0._wp                            ! k=1 only as zwx has been multiplied by wmask
255         ELSE                                         ! linear free surface
256            IF( ln_isfcav ) THEN                            ! ice-shelf cavities (top of the ocean)
257               DO jj = 1, jpj
258                  DO ji = 1, jpi
259                     zwx(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn)
260                  END DO
261               END DO   
262            ELSE                                             ! no cavities: only at the ocean surface
263               zwx(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn)
264            ENDIF
265         ENDIF
266         !
267         DO jk = 1, jpkm1                    ! Compute & add the vertical advective trend
268            DO jj = 2, jpjm1     
269               DO ji = fs_2, fs_jpim1   ! vector opt.
270                  pta(ji,jj,jk,jn) =  pta(ji,jj,jk,jn) - ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) ) / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) )
271               END DO
272            END DO
273         END DO
274         !                                 ! Save the vertical advective trends for diagnostic
275         IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR.     &
276            &( cdtype == 'TRC' .AND. l_trdtrc )      )   &
277            CALL trd_tra( kt, cdtype, jn, jptra_zad, zwx, pwn, ptb(:,:,:,jn) )
278         !
279      END DO
280      !
281      CALL wrk_dealloc( jpi,jpj,jpk,   zslpx, zslpy, zwx, zwy )
282      !
283      IF( nn_timing == 1 )  CALL timing_stop('tra_adv_mus')
284      !
285   END SUBROUTINE tra_adv_mus
286
287   !!======================================================================
288END MODULE traadv_mus
Note: See TracBrowser for help on using the repository browser.