MODULE traadv_ubs !!============================================================================== !! *** MODULE traadv_ubs *** !! Ocean active tracers: horizontal & vertical advective trend !!============================================================================== !! History : 1.0 ! 2006-08 (L. Debreu, R. Benshila) Original code !! 3.3 ! 2010-05 (C. Ethe, G. Madec) merge TRC-TRA + switch from velocity to transport !!---------------------------------------------------------------------- !!---------------------------------------------------------------------- !! tra_adv_ubs : update the tracer trend with the horizontal !! advection trends using a third order biaised scheme !!---------------------------------------------------------------------- USE oce ! ocean dynamics and active tracers USE dom_oce ! ocean space and time domain USE trdmod_oce ! ocean space and time domain USE trdtra USE lib_mpp USE lbclnk ! ocean lateral boundary condition (or mpp link) USE in_out_manager ! I/O manager USE diaptr ! poleward transport diagnostics USE dynspg_oce ! choice/control of key cpp for surface pressure gradient USE trc_oce ! share passive tracers/Ocean variables IMPLICIT NONE PRIVATE PUBLIC tra_adv_ubs ! routine called by traadv module LOGICAL :: l_trd ! flag to compute trends or not !! * Control permutation of array indices # include "oce_ftrans.h90" # include "dom_oce_ftrans.h90" # include "trc_oce_ftrans.h90" !! * Substitutions # include "domzgr_substitute.h90" # include "vectopt_loop_substitute.h90" !!---------------------------------------------------------------------- !! NEMO/OPA 3.3 , NEMO Consortium (2010) !! $Id$ !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) !!---------------------------------------------------------------------- CONTAINS SUBROUTINE tra_adv_ubs ( kt, cdtype, p2dt, pun, pvn, pwn, & & ptb, ptn, pta, kjpt ) !!---------------------------------------------------------------------- !! *** ROUTINE tra_adv_ubs *** !! !! ** Purpose : Compute the now trend due to the advection of tracers !! and add it to the general trend of passive tracer equations. !! !! ** Method : The upstream biased third (UBS) is order scheme based !! on an upstream-biased parabolic interpolation (Shchepetkin and McWilliams 2005) !! It is only used in the horizontal direction. !! For example the i-component of the advective fluxes are given by : !! ! e1u e3u un ( mi(Tn) - zltu(i ) ) if un(i) >= 0 !! zwx = ! or !! ! e1u e3u un ( mi(Tn) - zltu(i+1) ) if un(i) < 0 !! where zltu is the second derivative of the before temperature field: !! zltu = 1/e3t di[ e2u e3u / e1u di[Tb] ] !! This results in a dissipatively dominant (i.e. hyper-diffusive) !! truncation error. The overall performance of the advection scheme !! is similar to that reported in (Farrow and Stevens, 1995). !! For stability reasons, the first term of the fluxes which corresponds !! to a second order centered scheme is evaluated using the now velocity !! (centered in time) while the second term which is the diffusive part !! of the scheme, is evaluated using the before velocity (forward in time). !! Note that UBS is not positive. Do not use it on passive tracers. !! On the vertical, the advection is evaluated using a TVD scheme, as !! the UBS have been found to be too diffusive. !! !! ** Action : - update (pta) with the now advective tracer trends !! !! Reference : Shchepetkin, A. F., J. C. McWilliams, 2005, Ocean Modelling, 9, 347-404. !! Farrow, D.E., Stevens, D.P., 1995, J. Phys. Ocean. 25, 1731Ð1741. !!---------------------------------------------------------------------- USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released USE oce , ONLY: zwx => ua , zwy => va ! (ua,va) used as workspace USE wrk_nemo, ONLY: ztu => wrk_3d_1 , ztv => wrk_3d_2 ! 3D workspace USE wrk_nemo, ONLY: zltu => wrk_3d_3 , zltv => wrk_3d_4 ! - - USE wrk_nemo, ONLY: zti => wrk_3d_5 , ztw => wrk_3d_6 ! - - !! DCSE_NEMO: need additional directives for renamed module variables !FTRANS zwx zwy ztu ztv zltu zltv zti ztw :I :I :z ! INTEGER , INTENT(in ) :: kt ! ocean time-step index CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) INTEGER , INTENT(in ) :: kjpt ! number of tracers REAL(wp), DIMENSION( jpk ), INTENT(in ) :: p2dt ! vertical profile of tracer time-step ! REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pun, pvn, pwn ! 3 ocean velocity components ! REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb, ptn ! before and now tracer fields ! REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend !FTRANS pun pvn pwn :I :I :z !FTRANS ptb ptn pta :I :I :z : REAL(wp), INTENT(in ) :: pun(jpi,jpj,jpk) ! ocean velocity component (u) REAL(wp), INTENT(in ) :: pvn(jpi,jpj,jpk) ! ocean velocity component (v) REAL(wp), INTENT(in ) :: pwn(jpi,jpj,jpk) ! ocean velocity component (w) !! DCSE_NEMO: Next two arguments made inout to silence the cray compile, !! which rightly complains about the call to nonosc_v (which also has them !! as inout) REAL(wp), INTENT(inout) :: ptb(jpi,jpj,jpk,kjpt) ! tracer fields (before) REAL(wp), INTENT(inout) :: ptn(jpi,jpj,jpk,kjpt) ! tracer fields (now) REAL(wp), INTENT(inout) :: pta(jpi,jpj,jpk,kjpt) ! tracer trend ! INTEGER :: ji, jj, jk, jn ! dummy loop indices REAL(wp) :: ztra, zbtr, zcoef, z2dtt ! local scalars REAL(wp) :: zfp_ui, zfm_ui, zcenut, ztak, zfp_wk, zfm_wk ! - - REAL(wp) :: zfp_vj, zfm_vj, zcenvt, zeeu, zeev, z_hdivn ! - - !!---------------------------------------------------------------------- IF( wrk_in_use(3, 1,2,3,4,5,6) )THEN CALL ctl_stop('tra_adv_ubs: requested workspace arrays unavailable') ; RETURN ENDIF IF( kt == nit000 ) THEN IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) 'tra_adv_ubs : horizontal UBS advection scheme on ', cdtype IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' ! l_trd = .FALSE. IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE. ENDIF ! ! ! =========== DO jn = 1, kjpt ! tracer loop ! ! =========== ! 1. Bottom value : flux set to zero ! ---------------------------------- zltu(:,:,jpk) = 0.e0 ; zltv(:,:,jpk) = 0.e0 ! #if defined key_z_first DO jj = 1, jpjm1 DO ji = 1, jpim1 DO jk = 1, jpkm1 #else DO jk = 1, jpkm1 ! Horizontal slab ! ! Laplacian DO jj = 1, jpjm1 ! First derivative (gradient) DO ji = 1, fs_jpim1 ! vector opt. #endif zeeu = e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj) * umask(ji,jj,jk) zeev = e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj) * vmask(ji,jj,jk) ztu(ji,jj,jk) = zeeu * ( ptb(ji+1,jj ,jk,jn) - ptb(ji,jj,jk,jn) ) ztv(ji,jj,jk) = zeev * ( ptb(ji ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) ) END DO END DO #if defined key_z_first END DO DO jj = 2, jpjm1 ! Second derivative (divergence) DO ji = 2, jpim1 DO jk = 1, jpkm1 #else DO jj = 2, jpjm1 ! Second derivative (divergence) DO ji = fs_2, fs_jpim1 ! vector opt. #endif zcoef = 1. / ( 6. * fse3t(ji,jj,jk) ) zltu(ji,jj,jk) = ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk) ) * zcoef zltv(ji,jj,jk) = ( ztv(ji,jj,jk) - ztv(ji,jj-1,jk) ) * zcoef END DO END DO ! END DO ! End of slab CALL lbc_lnk( zltu, 'T', 1. ) ; CALL lbc_lnk( zltv, 'T', 1. ) ! Lateral boundary cond. (unchanged sgn) ! ! Horizontal advective fluxes #if defined key_z_first DO jj = 1, jpjm1 DO ji = 1, jpim1 DO jk = 1, jpkm1 #else DO jk = 1, jpkm1 ! Horizontal slab DO jj = 1, jpjm1 DO ji = 1, fs_jpim1 ! vector opt. #endif ! upstream transport zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) ) zfm_ui = pun(ji,jj,jk) - ABS( pun(ji,jj,jk) ) zfp_vj = pvn(ji,jj,jk) + ABS( pvn(ji,jj,jk) ) zfm_vj = pvn(ji,jj,jk) - ABS( pvn(ji,jj,jk) ) ! centered scheme zcenut = 0.5 * pun(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj ,jk,jn) ) zcenvt = 0.5 * pvn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji ,jj+1,jk,jn) ) ! UBS scheme zwx(ji,jj,jk) = zcenut - zfp_ui * zltu(ji,jj,jk) - zfm_ui * zltu(ji+1,jj,jk) zwy(ji,jj,jk) = zcenvt - zfp_vj * zltv(ji,jj,jk) - zfm_vj * zltv(ji,jj+1,jk) END DO END DO END DO ! End of slab zltu(:,:,:) = pta(:,:,:,jn) ! store pta trends ! Horizontal advective trends #if defined key_z_first DO jj = 2, jpjm1 DO ji = 2, jpim1 DO jk = 1, jpkm1 #else DO jk = 1, jpkm1 ! Tracer flux divergence at t-point added to the general trend DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. #endif zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) ! horizontal advective ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk) & & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk) ) ! add it to the general tracer trends pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra END DO END DO ! END DO ! End of slab ! Horizontal trend used in tra_adv_ztvd subroutine zltu(:,:,:) = pta(:,:,:,jn) - zltu(:,:,:) ! 3. Save the horizontal advective trends for diagnostic ! ------------------------------------------------------ ! ! trend diagnostics (contribution of upstream fluxes) IF( l_trd ) THEN CALL trd_tra( kt, cdtype, jn, jptra_trd_xad, zwx, pun, ptn(:,:,:,jn) ) CALL trd_tra( kt, cdtype, jn, jptra_trd_yad, zwy, pvn, ptn(:,:,:,jn) ) END IF ! ! "Poleward" heat and salt transports (contribution of upstream fluxes) IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN IF( jn == jp_tem ) htr_adv(:) = ptr_vj( zwy(:,:,:) ) IF( jn == jp_sal ) str_adv(:) = ptr_vj( zwy(:,:,:) ) ENDIF ! TVD scheme for the vertical direction ! ---------------------- IF( l_trd ) zltv(:,:,:) = pta(:,:,:,jn) ! store pta if trend diag. ! Bottom value : flux set to zero ztw(:,:,jpk) = 0.e0 ; zti(:,:,jpk) = 0.e0 ! Surface value IF( lk_vvl ) THEN ; ztw(:,:,1) = 0.e0 ! variable volume : flux set to zero ELSE ; ztw(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn) ! free constant surface ENDIF ! upstream advection with initial mass fluxes & intermediate update ! ------------------------------------------------------------------- ! Interior value #if defined key_z_first DO jj = 1, jpj DO ji = 1, jpi DO jk = 2, jpkm1 #else DO jk = 2, jpkm1 DO jj = 1, jpj DO ji = 1, jpi #endif zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) ) zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) ) ztw(ji,jj,jk) = 0.5 * ( zfp_wk * ptb(ji,jj,jk,jn) + zfm_wk * ptb(ji,jj,jk-1,jn) ) END DO END DO END DO ! update and guess with monotonic sheme #if defined key_z_first DO jj = 2, jpjm1 DO ji = 2, jpim1 DO jk = 1, jpkm1 z2dtt = p2dt(jk) #else DO jk = 1, jpkm1 z2dtt = p2dt(jk) DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. #endif zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) ztak = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * zbtr pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztak zti(ji,jj,jk) = ( ptb(ji,jj,jk,jn) + z2dtt * ( ztak + zltu(ji,jj,jk) ) ) * tmask(ji,jj,jk) END DO END DO END DO ! CALL lbc_lnk( zti, 'T', 1. ) ! Lateral boundary conditions on zti, zsi (unchanged sign) ! antidiffusive flux : high order minus low order #if defined key_z_first DO jj = 1, jpj DO ji = 1, jpi ztw(ji,jj,1) = 0.e0 ! Surface value DO jk = 2, jpkm1 ! Interior value #else ztw(:,:,1) = 0.e0 ! Surface value DO jk = 2, jpkm1 ! Interior value DO jj = 1, jpj DO ji = 1, jpi #endif ztw(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj,jk-1,jn) ) - ztw(ji,jj,jk) END DO END DO END DO ! CALL nonosc_z( ptb(:,:,:,jn), ztw, zti, p2dt ) ! monotonicity algorithm ! final trend with corrected fluxes #if defined key_z_first DO jj = 2, jpjm1 DO ji = 2, jpim1 DO jk = 1, jpkm1 #else DO jk = 1, jpkm1 DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. #endif zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) ! k- vertical advective trends ztra = - zbtr * ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) ! added to the general tracer trends pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra END DO END DO END DO ! Save the final vertical advective trends IF( l_trd ) THEN ! vertical advective trend diagnostics ! (compute -w.dk[ptn]= -dk[w.ptn] + ptn.dk[w]) #if defined key_z_first DO jj = 2, jpjm1 DO ji = 2, jpim1 DO jk = 1, jpkm1 #else DO jk = 1, jpkm1 DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. #endif zbtr = 1.e0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) z_hdivn = ( pwn(ji,jj,jk) - pwn(ji,jj,jk+1) ) * zbtr zltv(ji,jj,jk) = pta(ji,jj,jk,jn) - zltv(ji,jj,jk) + ptn(ji,jj,jk,jn) * z_hdivn END DO END DO END DO CALL trd_tra( kt, cdtype, jn, jptra_trd_zad, zltv ) ENDIF ! ENDDO ! IF( wrk_not_released(3, 1,2,3,4,5,6) ) CALL ctl_stop('tra_adv_ubs: failed to release workspace arrays') ! !! * Reset control of array index permutation !FTRANS CLEAR # include "oce_ftrans.h90" # include "dom_oce_ftrans.h90" # include "trc_oce_ftrans.h90" END SUBROUTINE tra_adv_ubs SUBROUTINE nonosc_z( pbef, pcc, paft, p2dt ) !!--------------------------------------------------------------------- !! *** ROUTINE nonosc_z *** !! !! ** Purpose : compute monotonic tracer fluxes from the upstream !! scheme and the before field by a nonoscillatory algorithm !! !! ** Method : ... ??? !! warning : pbef and paft must be masked, but the boundaries !! conditions on the fluxes are not necessary zalezak (1979) !! drange (1995) multi-dimensional forward-in-time and upstream- !! in-space based differencing for fluid !!---------------------------------------------------------------------- USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released USE wrk_nemo, ONLY: zbetup => wrk_3d_1, zbetdo => wrk_3d_2 ! 3D workspace !! DCSE_NEMO: need additional directives for renamed module variables !FTRANS zbetup zbetdo :I :I :z ! REAL(wp), INTENT(in ), DIMENSION(jpk) :: p2dt ! vertical profile of tracer time-step !! DCSE_NEMO: This style defeats ftrans ! REAL(wp), DIMENSION (jpi,jpj,jpk) :: pbef ! before field ! REAL(wp), INTENT(inout), DIMENSION (jpi,jpj,jpk) :: paft ! after field ! REAL(wp), INTENT(inout), DIMENSION (jpi,jpj,jpk) :: pcc ! monotonic flux in the k direction !FTRANS pbef paft pcc :I :I :z REAL(wp), INTENT(inout) :: pbef(jpi,jpj,jpk) ! before field REAL(wp), INTENT(inout) :: paft(jpi,jpj,jpk) ! after field REAL(wp), INTENT(inout) :: pcc(jpi,jpj,jpk) ! monotonic flux in the k direction ! INTEGER :: ji, jj, jk ! dummy loop indices INTEGER :: ikm1 ! local integer REAL(wp) :: zpos, zneg, zbt, za, zb, zc, zbig, zrtrn, z2dtt ! local scalars !!---------------------------------------------------------------------- IF( wrk_in_use(3, 1,2) ) THEN CALL ctl_stop('nonosc_z: requested workspace arrays unavailable') ; RETURN ENDIF zbig = 1.e+40_wp zrtrn = 1.e-15_wp zbetup(:,:,:) = 0._wp ; zbetdo(:,:,:) = 0._wp ! Search local extrema ! -------------------- ! large negative value (-zbig) inside land pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:) - zbig * ( 1.e0 - tmask(:,:,:) ) paft(:,:,:) = paft(:,:,:) * tmask(:,:,:) - zbig * ( 1.e0 - tmask(:,:,:) ) ! search maximum in neighbourhood #if defined key_z_first DO jj = 2, jpjm1 DO ji = 2, jpim1 DO jk = 1, jpkm1 ikm1 = MAX(jk-1,1) #else DO jk = 1, jpkm1 ikm1 = MAX(jk-1,1) DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. #endif zbetup(ji,jj,jk) = MAX( pbef(ji ,jj ,jk ), paft(ji ,jj ,jk ), & & pbef(ji ,jj ,ikm1), pbef(ji ,jj ,jk+1), & & paft(ji ,jj ,ikm1), paft(ji ,jj ,jk+1) ) END DO END DO END DO ! large positive value (+zbig) inside land pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:) + zbig * ( 1.e0 - tmask(:,:,:) ) paft(:,:,:) = paft(:,:,:) * tmask(:,:,:) + zbig * ( 1.e0 - tmask(:,:,:) ) ! search minimum in neighbourhood #if defined key_z_first DO jj = 2, jpjm1 DO ji = 2, jpim1 DO jk = 1, jpkm1 ikm1 = MAX(jk-1,1) #else DO jk = 1, jpkm1 ikm1 = MAX(jk-1,1) DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. #endif zbetdo(ji,jj,jk) = MIN( pbef(ji ,jj ,jk ), paft(ji ,jj ,jk ), & & pbef(ji ,jj ,ikm1), pbef(ji ,jj ,jk+1), & & paft(ji ,jj ,ikm1), paft(ji ,jj ,jk+1) ) END DO END DO END DO ! restore masked values to zero pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:) paft(:,:,:) = paft(:,:,:) * tmask(:,:,:) ! 2. Positive and negative part of fluxes and beta terms ! ------------------------------------------------------ #if defined key_z_first DO jj = 2, jpjm1 DO ji = 2, jpim1 DO jk = 1, jpkm1 z2dtt = p2dt(jk) #else DO jk = 1, jpkm1 z2dtt = p2dt(jk) DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. #endif ! positive & negative part of the flux zpos = MAX( 0., pcc(ji ,jj ,jk+1) ) - MIN( 0., pcc(ji ,jj ,jk ) ) zneg = MAX( 0., pcc(ji ,jj ,jk ) ) - MIN( 0., pcc(ji ,jj ,jk+1) ) ! up & down beta terms zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) / z2dtt zbetup(ji,jj,jk) = ( zbetup(ji,jj,jk) - paft(ji,jj,jk) ) / (zpos+zrtrn) * zbt zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zbetdo(ji,jj,jk) ) / (zneg+zrtrn) * zbt END DO END DO END DO ! monotonic flux in the k direction, i.e. pcc ! ------------------------------------------- #if defined key_z_first DO jj = 2, jpjm1 DO ji = 2, jpim1 DO jk = 2, jpkm1 #else DO jk = 2, jpkm1 DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. #endif za = MIN( 1., zbetdo(ji,jj,jk), zbetup(ji,jj,jk-1) ) zb = MIN( 1., zbetup(ji,jj,jk), zbetdo(ji,jj,jk-1) ) zc = 0.5 * ( 1.e0 + SIGN( 1.e0, pcc(ji,jj,jk) ) ) pcc(ji,jj,jk) = pcc(ji,jj,jk) * ( zc * za + ( 1.e0 - zc) * zb ) END DO END DO END DO ! IF( wrk_not_released(3, 1,2) ) CALL ctl_stop('nonosc_z: failed to release workspace arrays') ! END SUBROUTINE nonosc_z !!====================================================================== END MODULE traadv_ubs