MODULE traadv_ubs !!============================================================================== !! *** MODULE traadv_ubs *** !! Ocean active tracers: horizontal & vertical advective trend !!============================================================================== !! History : 9.0 ! 06-08 (L. Debreu, R. Benshila) Original code !!---------------------------------------------------------------------- !!---------------------------------------------------------------------- !! 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 USE trdmod_oce 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 prtctl IMPLICIT NONE PRIVATE PUBLIC tra_adv_ubs ! routine called by traadv module REAL(wp), DIMENSION(jpi,jpj) :: e1e2tr ! = 1/(e1t * e2t) !! * Substitutions # include "domzgr_substitute.h90" # include "vectopt_loop_substitute.h90" !!---------------------------------------------------------------------- !! OPA 9.0 , LOCEAN-IPSL (2006) !! $Id: traadv_ubs.F90 1528 2009-07-23 14:38:47Z rblod $ !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) !!---------------------------------------------------------------------- CONTAINS SUBROUTINE tra_adv_ubs( kt, pun, pvn, pwn ) !!---------------------------------------------------------------------- !! *** 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 (ta,sa) 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 oce, ONLY : zwx => ua ! use ua as workspace USE oce, ONLY : zwy => va ! use va as workspace !! INTEGER , INTENT(in) :: kt ! ocean time-step index REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) :: pun ! effective ocean velocity, u_component REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) :: pvn ! effective ocean velocity, v_component REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) :: pwn ! effective ocean velocity, w_component !! INTEGER :: ji, jj, jk ! dummy loop indices REAL(wp) :: zta, zsa, zbtr, zcoef ! temporary scalars REAL(wp) :: zfui, zfp_ui, zfm_ui, zcenut, zcenus ! " " REAL(wp) :: zfvj, zfp_vj, zfm_vj, zcenvt, zcenvs ! " " REAL(wp) :: z_hdivn_x, z_hdivn_y, z_hdivn ! " " REAL(wp), DIMENSION(jpi,jpj) :: zeeu, zeev ! temporary 2D workspace REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwz , zww ! temporary 3D workspace REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztu , ztv , zltu , zltv, ztrdt ! " " REAL(wp), DIMENSION(jpi,jpj,jpk) :: zsu , zsv , zlsu , zlsv, ztrds ! " " !!---------------------------------------------------------------------- zltu(:,:,:) = 0.e0 zltv(:,:,:) = 0.e0 zlsu(:,:,:) = 0.e0 zlsv(:,:,:) = 0.e0 IF( kt == nit000 ) THEN IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) 'tra_adv_ubs : horizontal UBS advection scheme' IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' ! e1e2tr(:,:) = 1. / ( e1t(:,:) * e2t(:,:) ) ENDIF ! Save ta and sa trends ztrdt(:,:,:) = ta(:,:,:) ztrds(:,:,:) = sa(:,:,:) zcoef = 1./6. ! ! =============== DO jk = 1, jpkm1 ! Horizontal slab ! ! =============== ! Initialization of metric arrays (for z- or s-coordinates) DO jj = 1, jpjm1 DO ji = 1, fs_jpim1 ! vector opt. #if defined key_zco ! z-coordinates, no vertical scale factors zeeu(ji,jj) = e2u(ji,jj) / e1u(ji,jj) * umask(ji,jj,jk) zeev(ji,jj) = e1v(ji,jj) / e2v(ji,jj) * vmask(ji,jj,jk) #else ! s-coordinates, vertical scale factor are used zeeu(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj) * umask(ji,jj,jk) zeev(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj) * vmask(ji,jj,jk) #endif END DO END DO ! Laplacian ! First derivative (gradient) DO jj = 1, jpjm1 DO ji = 1, fs_jpim1 ! vector opt. ztu(ji,jj,jk) = zeeu(ji,jj) * ( tb(ji+1,jj ,jk) - tb(ji,jj,jk) ) zsu(ji,jj,jk) = zeeu(ji,jj) * ( sb(ji+1,jj ,jk) - sb(ji,jj,jk) ) ztv(ji,jj,jk) = zeev(ji,jj) * ( tb(ji ,jj+1,jk) - tb(ji,jj,jk) ) zsv(ji,jj,jk) = zeev(ji,jj) * ( sb(ji ,jj+1,jk) - sb(ji,jj,jk) ) END DO END DO ! Second derivative (divergence) DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. #if ! defined key_zco zcoef = 1. / ( 6. * fse3t(ji,jj,jk) ) #endif zltu(ji,jj,jk) = ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk) ) * zcoef zlsu(ji,jj,jk) = ( zsu(ji,jj,jk) - zsu(ji-1,jj,jk) ) * zcoef zltv(ji,jj,jk) = ( ztv(ji,jj,jk) - ztv(ji,jj-1,jk) ) * zcoef zlsv(ji,jj,jk) = ( zsv(ji,jj,jk) - zsv(ji,jj-1,jk) ) * zcoef END DO END DO ! ! ================= END DO ! End of slab ! ! ================= ! Lateral boundary conditions on the laplacian (zlt,zls) (unchanged sgn) CALL lbc_lnk( zltu, 'T', 1. ) ; CALL lbc_lnk( zlsu, 'T', 1. ) CALL lbc_lnk( zltv, 'T', 1. ) ; CALL lbc_lnk( zlsv, 'T', 1. ) ! ! =============== DO jk = 1, jpkm1 ! Horizontal slab ! ! =============== ! Horizontal advective fluxes DO jj = 1, jpjm1 DO ji = 1, fs_jpim1 ! vector opt. ! volume fluxes * 1/2 #if defined key_zco zfui = 0.5 * e2u(ji,jj) * pun(ji,jj,jk) zfvj = 0.5 * e1v(ji,jj) * pvn(ji,jj,jk) #else zfui = 0.5 * e2u(ji,jj) * fse3u(ji,jj,jk) * pun(ji,jj,jk) zfvj = 0.5 * e1v(ji,jj) * fse3v(ji,jj,jk) * pvn(ji,jj,jk) #endif ! upstream scheme zfp_ui = zfui + ABS( zfui ) zfp_vj = zfvj + ABS( zfvj ) zfm_ui = zfui - ABS( zfui ) zfm_vj = zfvj - ABS( zfvj ) ! centered scheme zcenut = zfui * ( tn(ji,jj,jk) + tn(ji+1,jj ,jk) ) zcenvt = zfvj * ( tn(ji,jj,jk) + tn(ji ,jj+1,jk) ) zcenus = zfui * ( sn(ji,jj,jk) + sn(ji+1,jj ,jk) ) zcenvs = zfvj * ( sn(ji,jj,jk) + sn(ji ,jj+1,jk) ) ! mixed centered / upstream 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) zww(ji,jj,jk) = zcenus - zfp_ui * zlsu(ji,jj,jk) -zfm_ui * zlsu(ji+1,jj,jk) zwz(ji,jj,jk) = zcenvs - zfp_vj * zlsv(ji,jj,jk) -zfm_vj * zlsv(ji,jj+1,jk) END DO END DO ! Tracer flux divergence at t-point added to the general trend DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. ! horizontal advective trends #if defined key_zco zbtr = e1e2tr(ji,jj) #else zbtr = e1e2tr(ji,jj) / fse3t(ji,jj,jk) #endif zta = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk) & & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk) ) zsa = - zbtr * ( zww(ji,jj,jk) - zww(ji-1,jj ,jk) & & + zwz(ji,jj,jk) - zwz(ji ,jj-1,jk) ) ! add it to the general tracer trends ta(ji,jj,jk) = ta(ji,jj,jk) + zta sa(ji,jj,jk) = sa(ji,jj,jk) + zsa END DO END DO ! ! =============== END DO ! End of slab ! ! =============== ! Horizontal trend used in tra_adv_ztvd subroutine zltu(:,:,:) = ta(:,:,:) - ztrdt(:,:,:) zlsu(:,:,:) = sa(:,:,:) - ztrds(:,:,:) ! 3. Save the horizontal advective trends for diagnostic ! ------------------------------------------------------ IF( l_trdtra ) THEN ! Recompute the hoizontal advection zta & zsa trends computed ! at the step 2. above in making the difference between the new ! trends and the previous one ta()/sa - ztrdt()/ztrds() and add ! the term tn()/sn()*hdivn() to recover the Uh gradh(T/S) trends ztrdt(:,:,:) = 0.e0 ; ztrds(:,:,:) = 0.e0 ! ! T/S ZONAL advection trends DO jk = 1, jpkm1 DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. !-- Compute zonal divergence by splitting hdivn (see divcur.F90) #if defined key_zco zbtr = e1e2tr(ji,jj) z_hdivn_x = ( e2u(ji ,jj) * pun(ji ,jj,jk) & & - e2u(ji-1,jj) * pun(ji-1,jj,jk) ) * zbtr #else zbtr = e1e2tr(ji,jj) / fse3t(ji,jj,jk) z_hdivn_x = ( e2u(ji ,jj) * fse3u(ji ,jj,jk) * pun(ji ,jj,jk) & & - e2u(ji-1,jj) * fse3u(ji-1,jj,jk) * pun(ji-1,jj,jk) ) * zbtr #endif ztrdt(ji,jj,jk) = - ( zwx(ji,jj,jk) - zwx(ji-1,jj,jk) ) * zbtr + tn(ji,jj,jk) * z_hdivn_x ztrds(ji,jj,jk) = - ( zww(ji,jj,jk) - zww(ji-1,jj,jk) ) * zbtr + sn(ji,jj,jk) * z_hdivn_x END DO END DO END DO CALL trd_mod(ztrdt, ztrds, jptra_trd_xad, 'TRA', kt) ! save the trends ! ! T/S MERIDIONAL advection trends DO jk = 1, jpkm1 DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. !-- Compute merid. divergence by splitting hdivn (see divcur.F90) #if defined key_zco zbtr = e1e2tr(ji,jj) z_hdivn_y = ( e1v(ji, jj) * pvn(ji,jj ,jk) & & - e1v(ji,jj-1) * pvn(ji,jj-1,jk) ) * zbtr #else zbtr = e1e2tr(ji,jj) / fse3t(ji,jj,jk) z_hdivn_y = ( e1v(ji, jj) * fse3v(ji,jj ,jk) * pvn(ji,jj ,jk) & & - e1v(ji,jj-1) * fse3v(ji,jj-1,jk) * pvn(ji,jj-1,jk) ) * zbtr #endif ztrdt(ji,jj,jk) = - ( zwy(ji,jj,jk) - zwy(ji,jj-1,jk) ) * zbtr + tn(ji,jj,jk) * z_hdivn_y ztrds(ji,jj,jk) = - ( zwz(ji,jj,jk) - zwz(ji,jj-1,jk) ) * zbtr + sn(ji,jj,jk) * z_hdivn_y END DO END DO END DO CALL trd_mod(ztrdt, ztrds, jptra_trd_yad, 'TRA', kt) ! save the trends ! ENDIF IF(ln_ctl) CALL prt_ctl( tab3d_1=ta, clinfo1=' ubs had - Ta: ', mask1=tmask, & & tab3d_2=sa, clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' ) ! "zonal" mean advective heat and salt transport IF( ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN IF( lk_zco ) THEN DO jk = 1, jpkm1 DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. zwy(ji,jj,jk) = zwy(ji,jj,jk) * fse3v(ji,jj,jk) zwz(ji,jj,jk) = zwz(ji,jj,jk) * fse3v(ji,jj,jk) END DO END DO END DO ENDIF pht_adv(:) = ptr_vj( zwy(:,:,:) ) pst_adv(:) = ptr_vj( zwz(:,:,:) ) ENDIF ! II. Vertical advection ! ---------------------- IF( l_trdtra ) THEN ! Save ta and sa trends ztrdt(:,:,:) = ta(:,:,:) ztrds(:,:,:) = sa(:,:,:) ENDIF ! TVD scheme the vertical direction CALL tra_adv_ztvd(kt, pwn, zltu, zlsu) IF( l_trdtra ) THEN ! Save the final vertical advective trends DO jk = 1, jpkm1 DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. #if defined key_zco zbtr = e1e2tr(ji,jj) z_hdivn_x = e2u(ji,jj)*pun(ji,jj,jk) - e2u(ji-1,jj)*pun(ji-1,jj,jk) z_hdivn_y = e1v(ji,jj)*pvn(ji,jj,jk) - e1v(ji,jj-1)*pvn(ji,jj-1,jk) #else zbtr = e1e2tr(ji,jj) / fse3t(ji,jj,jk) z_hdivn_x = e2u(ji,jj)*fse3u(ji,jj,jk)*pun(ji,jj,jk) - e2u(ji-1,jj)*fse3u(ji-1,jj,jk)*pun(ji-1,jj,jk) z_hdivn_y = e1v(ji,jj)*fse3v(ji,jj,jk)*pvn(ji,jj,jk) - e1v(ji,jj-1)*fse3v(ji,jj-1,jk)*pvn(ji,jj-1,jk) #endif z_hdivn = (z_hdivn_x + z_hdivn_y) * zbtr zbtr = e1e2tr(ji,jj) / fse3t(ji,jj,jk) ztrdt(ji,jj,jk) = ta(ji,jj,jk) - ztrdt(ji,jj,jk) - tn(ji,jj,jk) * z_hdivn ztrds(ji,jj,jk) = sa(ji,jj,jk) - ztrds(ji,jj,jk) - sn(ji,jj,jk) * z_hdivn END DO END DO END DO CALL trd_mod(ztrdt, ztrds, jptra_trd_zad, 'TRA', kt) ! <<< ADD TO PREVIOUSLY COMPUTED ! ENDIF IF(ln_ctl) CALL prt_ctl( tab3d_1=ta, clinfo1=' ubs zad - Ta: ', mask1=tmask, & & tab3d_2=sa, clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra') ! END SUBROUTINE tra_adv_ubs SUBROUTINE tra_adv_ztvd( kt, pwn, zttrd, zstrd ) !!---------------------------------------------------------------------- !! *** ROUTINE tra_adv_ztvd *** !! !! ** Purpose : Compute the now trend due to total advection of !! tracers and add it to the general trend of tracer equations !! !! ** Method : TVD scheme, i.e. 2nd order centered scheme with !! corrected flux (monotonic correction) !! note: - this advection scheme needs a leap-frog time scheme !! !! ** Action : - update (ta,sa) with the now advective tracer trends !! - save the trends in (ztrdt,ztrds) ('key_trdtra') !!---------------------------------------------------------------------- INTEGER , INTENT(in) :: kt ! ocean time-step REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) :: pwn ! verical effective velocity REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) :: zttrd, zstrd ! lateral advective trends on T & S !! INTEGER :: ji, jj, jk ! dummy loop indices REAL(wp) :: z2dtt, zbtr, zew, z2 ! temporary scalar REAL(wp) :: ztak, zfp_wk ! " " REAL(wp) :: zsak, zfm_wk ! " " REAL(wp), DIMENSION (jpi,jpj,jpk) :: zti, ztw ! temporary 3D workspace REAL(wp), DIMENSION (jpi,jpj,jpk) :: zsi, zsw ! " " !!---------------------------------------------------------------------- IF( kt == nit000 .AND. lwp ) THEN WRITE(numout,*) WRITE(numout,*) 'tra_adv_ztvd : vertical TVD advection scheme' WRITE(numout,*) '~~~~~~~~~~~~' ENDIF IF( neuler == 0 .AND. kt == nit000 ) THEN ; z2 = 1. ELSE ; z2 = 2. ENDIF ! Bottom value : flux set to zero ! -------------- ztw(:,:,jpk) = 0.e0 ; zsw(:,:,jpk) = 0.e0 zti (:,:,:) = 0.e0 ; zsi (:,:,:) = 0.e0 ! upstream advection with initial mass fluxes & intermediate update ! ------------------------------------------------------------------- ! Surface value IF( lk_vvl ) THEN ! variable volume : flux set to zero ztw(:,:,1) = 0.e0 zsw(:,:,1) = 0.e0 ELSE ! free surface-constant volume DO jj = 1, jpj DO ji = 1, jpi zew = e1t(ji,jj) * e2t(ji,jj) * pwn(ji,jj,1) ztw(ji,jj,1) = zew * tb(ji,jj,1) zsw(ji,jj,1) = zew * sb(ji,jj,1) END DO END DO ENDIF ! Interior value DO jk = 2, jpkm1 DO jj = 1, jpj DO ji = 1, jpi zew = 0.5 * e1t(ji,jj) * e2t(ji,jj) * pwn(ji,jj,jk) zfp_wk = zew + ABS( zew ) zfm_wk = zew - ABS( zew ) ztw(ji,jj,jk) = zfp_wk * tb(ji,jj,jk) + zfm_wk * tb(ji,jj,jk-1) zsw(ji,jj,jk) = zfp_wk * sb(ji,jj,jk) + zfm_wk * sb(ji,jj,jk-1) END DO END DO END DO ! update and guess with monotonic sheme DO jk = 1, jpkm1 z2dtt = z2 * rdttra(jk) DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. zbtr = 1./ ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) ztak = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * zbtr zsak = - ( zsw(ji,jj,jk) - zsw(ji,jj,jk+1) ) * zbtr ta(ji,jj,jk) = ta(ji,jj,jk) + ztak sa(ji,jj,jk) = sa(ji,jj,jk) + zsak zti (ji,jj,jk) = ( tb(ji,jj,jk) + z2dtt * ( ztak + zttrd(ji,jj,jk) ) ) * tmask(ji,jj,jk) zsi (ji,jj,jk) = ( sb(ji,jj,jk) + z2dtt * ( zsak + zstrd(ji,jj,jk) ) ) * tmask(ji,jj,jk) END DO END DO END DO ! Lateral boundary conditions on zti, zsi (unchanged sign) CALL lbc_lnk( zti, 'T', 1. ) CALL lbc_lnk( zsi, 'T', 1. ) ! antidiffusive flux : high order minus low order ! ------------------------------------------------- ! Surface value ztw(:,:,1) = 0.e0 ; zsw(:,:,1) = 0.e0 ! Interior value DO jk = 2, jpkm1 DO jj = 1, jpj DO ji = 1, jpi zew = 0.5 * e1t(ji,jj) * e2t(ji,jj) * pwn(ji,jj,jk) ztw(ji,jj,jk) = zew * ( tn(ji,jj,jk) + tn(ji,jj,jk-1) ) - ztw(ji,jj,jk) zsw(ji,jj,jk) = zew * ( sn(ji,jj,jk) + sn(ji,jj,jk-1) ) - zsw(ji,jj,jk) END DO END DO END DO ! monotonicity algorithm ! ------------------------ CALL nonosc_z( tb, ztw, zti, z2 ) CALL nonosc_z( sb, zsw, zsi, z2 ) ! final trend with corrected fluxes ! ----------------------------------- DO jk = 1, jpkm1 DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) ! k- vertical advective trends ztak = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * zbtr zsak = - ( zsw(ji,jj,jk) - zsw(ji,jj,jk+1) ) * zbtr ! add them to the general tracer trends ta(ji,jj,jk) = ta(ji,jj,jk) + ztak sa(ji,jj,jk) = sa(ji,jj,jk) + zsak END DO END DO END DO ! END SUBROUTINE tra_adv_ztvd SUBROUTINE nonosc_z( pbef, pcc, paft, prdt ) !!--------------------------------------------------------------------- !! *** 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 !!---------------------------------------------------------------------- REAL(wp), INTENT(in ) :: prdt ! ??? REAL(wp), INTENT(inout), 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 !! INTEGER :: ji, jj, jk ! dummy loop indices INTEGER :: ikm1 REAL(wp) :: zpos, zneg, zbt, za, zb, zc, zbig, zrtrn, z2dtt REAL(wp), DIMENSION (jpi,jpj,jpk) :: zbetup, zbetdo !!---------------------------------------------------------------------- zbig = 1.e+40 zrtrn = 1.e-15 zbetup(:,:,:) = 0.e0 ; zbetdo(:,:,:) = 0.e0 ! 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 DO jk = 1, jpkm1 ikm1 = MAX(jk-1,1) DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. 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 DO jk = 1, jpkm1 ikm1 = MAX(jk-1,1) DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. 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 ! ------------------------------------------------------ DO jk = 1, jpkm1 z2dtt = prdt * rdttra(jk) DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. ! 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 ! ------------------------------------------- DO jk = 2, jpkm1 DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. 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 ! END SUBROUTINE nonosc_z !!====================================================================== END MODULE traadv_ubs