MODULE traadv_muscl !!====================================================================== !! *** MODULE traadv_muscl *** !! Ocean active tracers: horizontal & vertical advective trend !!====================================================================== !! History : ! 06-00 (A.Estublier) for passive tracers !! ! 01-08 (E.Durand, G.Madec) adapted for T & S !! 8.5 ! 02-06 (G. Madec) F90: Free form and module !!---------------------------------------------------------------------- !!---------------------------------------------------------------------- !! tra_adv_muscl : update the tracer trend with the horizontal !! and vertical advection trends using MUSCL scheme !!---------------------------------------------------------------------- USE oce ! ocean dynamics and active tracers USE dom_oce ! ocean space and time domain USE trdmod ! ocean active tracers trends USE trdmod_oce ! ocean variables trends USE in_out_manager ! I/O manager USE dynspg_oce ! choice/control of key cpp for surface pressure gradient USE trabbl ! tracers: bottom boundary layer USE lib_mpp ! distribued memory computing USE lbclnk ! ocean lateral boundary condition (or mpp link) USE diaptr ! poleward transport diagnostics USE prtctl ! Print control IMPLICIT NONE PRIVATE PUBLIC tra_adv_muscl ! routine called by step.F90 !! * Substitutions # include "domzgr_substitute.h90" # include "vectopt_loop_substitute.h90" !!---------------------------------------------------------------------- !! OPA 9.0 , LOCEAN-IPSL (2006) !! $Header$ !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) !!---------------------------------------------------------------------- CONTAINS SUBROUTINE tra_adv_muscl( kt, pun, pvn, pwn ) !!---------------------------------------------------------------------- !! *** ROUTINE tra_adv_muscl *** !! !! ** Purpose : Compute the now trend due to total advection of T and !! S using a MUSCL scheme (Monotone Upstream-centered Scheme for !! Conservation Laws) and add it to the general tracer trend. !! !! ** Method : MUSCL scheme plus centered scheme at ocean boundaries !! !! ** Action : - update (ta,sa) with the now advective tracer trends !! - save trends in (ztrdt,ztrds) ('key_trdtra') !! !! References : Estubier, A., and M. Levy, Notes Techn. Pole de Modelisation !! IPSL, Sept. 2000 (http://www.lodyc.jussieu.fr/opa) !!---------------------------------------------------------------------- USE oce, ONLY : ztrdt => ua ! use ua as workspace USE oce, ONLY : ztrds => va ! use va as workspace !! INTEGER , INTENT(in) :: kt ! ocean time-step index REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) :: pun ! ocean velocity u-component REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) :: pvn ! ocean velocity v-component REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) :: pwn ! ocean velocity w-component !! INTEGER :: ji, jj, jk ! dummy loop indices REAL(wp) :: & zu, zv, zw, zeu, zev, & zew, zbtr, zstep, & z0u, z0v, z0w, & zzt1, zzt2, zalpha, & zzs1, zzs2, z2, & zta, zsa, & z_hdivn_x, z_hdivn_y, z_hdivn REAL(wp), DIMENSION (jpi,jpj,jpk) :: zt1, zt2, ztp1, ztp2 ! 3D workspace REAL(wp), DIMENSION (jpi,jpj,jpk) :: zs1, zs2, zsp1, zsp2 ! " " !!---------------------------------------------------------------------- IF( kt == nit000 .AND. lwp ) THEN WRITE(numout,*) WRITE(numout,*) 'tra_adv : MUSCL advection scheme' WRITE(numout,*) '~~~~~~~' ENDIF IF( neuler == 0 .AND. kt == nit000 ) THEN ; z2 = 1. ELSE ; z2 = 2. ENDIF ! I. Horizontal advective fluxes ! ------------------------------ ! first guess of the slopes ! interior values DO jk = 1, jpkm1 DO jj = 1, jpjm1 DO ji = 1, fs_jpim1 ! vector opt. zt1(ji,jj,jk) = umask(ji,jj,jk) * ( tb(ji+1,jj,jk) - tb(ji,jj,jk) ) zs1(ji,jj,jk) = umask(ji,jj,jk) * ( sb(ji+1,jj,jk) - sb(ji,jj,jk) ) zt2(ji,jj,jk) = vmask(ji,jj,jk) * ( tb(ji,jj+1,jk) - tb(ji,jj,jk) ) zs2(ji,jj,jk) = vmask(ji,jj,jk) * ( sb(ji,jj+1,jk) - sb(ji,jj,jk) ) END DO END DO END DO ! bottom values zt1(:,:,jpk) = 0.e0 ; zt2(:,:,jpk) = 0.e0 zs1(:,:,jpk) = 0.e0 ; zs2(:,:,jpk) = 0.e0 ! lateral boundary conditions on zt1, zt2 ; zs1, zs2 (changed sign) CALL lbc_lnk( zt1, 'U', -1. ) ; CALL lbc_lnk( zs1, 'U', -1. ) CALL lbc_lnk( zt2, 'V', -1. ) ; CALL lbc_lnk( zs2, 'V', -1. ) ! Slopes ! interior values DO jk = 1, jpkm1 DO jj = 2, jpj DO ji = fs_2, jpi ! vector opt. ztp1(ji,jj,jk) = ( zt1(ji,jj,jk) + zt1(ji-1,jj ,jk) ) & & * ( 0.25 + SIGN( 0.25, zt1(ji,jj,jk) * zt1(ji-1,jj ,jk) ) ) zsp1(ji,jj,jk) = ( zs1(ji,jj,jk) + zs1(ji-1,jj ,jk) ) & & * ( 0.25 + SIGN( 0.25, zs1(ji,jj,jk) * zs1(ji-1,jj ,jk) ) ) ztp2(ji,jj,jk) = ( zt2(ji,jj,jk) + zt2(ji ,jj-1,jk) ) & & * ( 0.25 + SIGN( 0.25, zt2(ji,jj,jk) * zt2(ji ,jj-1,jk) ) ) zsp2(ji,jj,jk) = ( zs2(ji,jj,jk) + zs2(ji ,jj-1,jk) ) & & * ( 0.25 + SIGN( 0.25, zs2(ji,jj,jk) * zs2(ji ,jj-1,jk) ) ) END DO END DO END DO ! bottom values ztp1(:,:,jpk) = 0.e0 ; ztp2(:,:,jpk) = 0.e0 zsp1(:,:,jpk) = 0.e0 ; zsp2(:,:,jpk) = 0.e0 ! Slopes limitation DO jk = 1, jpkm1 DO jj = 2, jpj DO ji = fs_2, jpi ! vector opt. ztp1(ji,jj,jk) = SIGN( 1., ztp1(ji,jj,jk) ) & & * MIN( ABS( ztp1(ji ,jj,jk) ), & & 2.*ABS( zt1 (ji-1,jj,jk) ), & & 2.*ABS( zt1 (ji ,jj,jk) ) ) zsp1(ji,jj,jk) = SIGN( 1., zsp1(ji,jj,jk) ) & & * MIN( ABS( zsp1(ji ,jj,jk) ), & & 2.*ABS( zs1 (ji-1,jj,jk) ), & & 2.*ABS( zs1 (ji ,jj,jk) ) ) ztp2(ji,jj,jk) = SIGN( 1., ztp2(ji,jj,jk) ) & & * MIN( ABS( ztp2(ji,jj ,jk) ), & & 2.*ABS( zt2 (ji,jj-1,jk) ), & & 2.*ABS( zt2 (ji,jj ,jk) ) ) zsp2(ji,jj,jk) = SIGN( 1., zsp2(ji,jj,jk) ) & & * MIN( ABS( zsp2(ji,jj ,jk) ), & & 2.*ABS( zs2 (ji,jj-1,jk) ), & & 2.*ABS( zs2 (ji,jj ,jk) ) ) END DO END DO END DO ! Advection terms ! interior values DO jk = 1, jpkm1 zstep = z2 * rdttra(jk) DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. ! volume fluxes #if defined key_zco zeu = e2u(ji,jj) * pun(ji,jj,jk) zev = e1v(ji,jj) * pvn(ji,jj,jk) #else zeu = e2u(ji,jj) * fse3u(ji,jj,jk) * pun(ji,jj,jk) zev = e1v(ji,jj) * fse3v(ji,jj,jk) * pvn(ji,jj,jk) #endif ! MUSCL fluxes z0u = SIGN( 0.5, pun(ji,jj,jk) ) zalpha = 0.5 - z0u zu = z0u - 0.5 * pun(ji,jj,jk) * zstep / e1u(ji,jj) zzt1 = tb(ji+1,jj,jk) + zu*ztp1(ji+1,jj,jk) zzt2 = tb(ji ,jj,jk) + zu*ztp1(ji ,jj,jk) zzs1 = sb(ji+1,jj,jk) + zu*zsp1(ji+1,jj,jk) zzs2 = sb(ji ,jj,jk) + zu*zsp1(ji ,jj,jk) zt1(ji,jj,jk) = zeu * ( zalpha * zzt1 + (1.-zalpha) * zzt2 ) zs1(ji,jj,jk) = zeu * ( zalpha * zzs1 + (1.-zalpha) * zzs2 ) ! z0v = SIGN( 0.5, pvn(ji,jj,jk) ) zalpha = 0.5 - z0v zv = z0v - 0.5 * pvn(ji,jj,jk) * zstep / e2v(ji,jj) zzt1 = tb(ji,jj+1,jk) + zv*ztp2(ji,jj+1,jk) zzt2 = tb(ji,jj ,jk) + zv*ztp2(ji,jj ,jk) zzs1 = sb(ji,jj+1,jk) + zv*zsp2(ji,jj+1,jk) zzs2 = sb(ji,jj ,jk) + zv*zsp2(ji,jj ,jk) zt2(ji,jj,jk) = zev * ( zalpha * zzt1 + (1.-zalpha) * zzt2 ) zs2(ji,jj,jk) = zev * ( zalpha * zzs1 + (1.-zalpha) * zzs2 ) END DO END DO END DO ! lateral boundary conditions on zt1, zt2 ; zs1, zs2 (changed sign) CALL lbc_lnk( zt1, 'U', -1. ) ; CALL lbc_lnk( zs1, 'U', -1. ) CALL lbc_lnk( zt2, 'V', -1. ) ; CALL lbc_lnk( zs2, 'V', -1. ) ! Tracer flux divergence at t-point added to the general trend DO jk = 1, jpkm1 DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. #if defined key_zco zbtr = 1. / ( e1t(ji,jj)*e2t(ji,jj) ) #else zbtr = 1. / ( e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk) ) #endif ! horizontal advective trends zta = - zbtr * ( zt1(ji,jj,jk) - zt1(ji-1,jj ,jk ) & & + zt2(ji,jj,jk) - zt2(ji ,jj-1,jk ) ) zsa = - zbtr * ( zs1(ji,jj,jk) - zs1(ji-1,jj ,jk ) & & + zs2(ji,jj,jk) - zs2(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 IF(ln_ctl) CALL prt_ctl( tab3d_1=ta, clinfo1=' muscl had - Ta: ', mask1=tmask , & & tab3d_2=sa, clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' ) ! Save the horizontal advective trends for diagnostics IF( l_trdtra ) THEN 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) ! N.B. This computation is not valid along OBCs (if any) #if defined key_zco zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) ) z_hdivn_x = ( e2u(ji ,jj) * pun(ji ,jj,jk) & & - e2u(ji-1,jj) * pun(ji-1,jj,jk) ) * zbtr #else zbtr = 1. / ( e1t(ji,jj) * e2t(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) = - zbtr * ( zt1(ji,jj,jk) - zt1(ji-1,jj,jk) ) + tn(ji,jj,jk) * z_hdivn_x ztrds(ji,jj,jk) = - zbtr * ( zs1(ji,jj,jk) - zs1(ji-1,jj,jk) ) + sn(ji,jj,jk) * z_hdivn_x END DO END DO END DO CALL trd_mod(ztrdt, ztrds, jptra_trd_xad, 'TRA', kt) ! 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) ! N.B. This computation is not valid along OBCs (if any) #if defined key_zco zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) ) z_hdivn_y = ( e1v(ji,jj ) * pvn(ji,jj ,jk) & & - e1v(ji,jj-1) * pvn(ji,jj-1,jk) ) * zbtr #else zbtr = 1. / ( e1t(ji,jj) * e2t(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) = - zbtr * ( zt2(ji,jj,jk) - zt2(ji,jj-1,jk) ) + tn(ji,jj,jk) * z_hdivn_y ztrds(ji,jj,jk) = - zbtr * ( zs2(ji,jj,jk) - zs2(ji,jj-1,jk) ) + 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 up-to-date ta and sa trends ztrdt(:,:,:) = ta(:,:,:) ztrds(:,:,:) = sa(:,:,:) ! ENDIF ! "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. zt2(ji,jj,jk) = zt2(ji,jj,jk) * fse3v(ji,jj,jk) zs2(ji,jj,jk) = zs2(ji,jj,jk) * fse3v(ji,jj,jk) END DO END DO END DO ENDIF pht_adv(:) = ptr_vj( zt2(:,:,:) ) pst_adv(:) = ptr_vj( zs2(:,:,:) ) ENDIF ! II. Vertical advective fluxes ! ----------------------------- ! First guess of the slope ! interior values DO jk = 2, jpkm1 zt1(:,:,jk) = tmask(:,:,jk) * ( tb(:,:,jk-1) - tb(:,:,jk) ) zs1(:,:,jk) = tmask(:,:,jk) * ( sb(:,:,jk-1) - sb(:,:,jk) ) END DO ! surface & bottom boundary conditions zt1 (:,:, 1 ) = 0.e0 ; zt1 (:,:,jpk) = 0.e0 zs1 (:,:, 1 ) = 0.e0 ; zs1 (:,:,jpk) = 0.e0 ! Slopes DO jk = 2, jpkm1 DO jj = 1, jpj DO ji = 1, jpi ztp1(ji,jj,jk) = ( zt1(ji,jj,jk) + zt1(ji,jj,jk+1) ) & & * ( 0.25 + SIGN( 0.25, zt1(ji,jj,jk) * zt1(ji,jj,jk+1) ) ) zsp1(ji,jj,jk) = ( zs1(ji,jj,jk) + zs1(ji,jj,jk+1) ) & & * ( 0.25 + SIGN( 0.25, zs1(ji,jj,jk) * zs1(ji,jj,jk+1) ) ) END DO END DO END DO ! Slopes limitation ! interior values DO jk = 2, jpkm1 DO jj = 1, jpj DO ji = 1, jpi ztp1(ji,jj,jk) = SIGN( 1., ztp1(ji,jj,jk) ) & & * MIN( ABS( ztp1(ji,jj,jk ) ), & & 2.*ABS( zt1 (ji,jj,jk+1) ), & & 2.*ABS( zt1 (ji,jj,jk ) ) ) zsp1(ji,jj,jk) = SIGN( 1., zsp1(ji,jj,jk) ) & & * MIN( ABS( zsp1(ji,jj,jk ) ), & & 2.*ABS( zs1 (ji,jj,jk+1) ), & & 2.*ABS( zs1 (ji,jj,jk ) ) ) END DO END DO END DO ! surface values ztp1(:,:,1) = 0.e0 zsp1(:,:,1) = 0.e0 ! vertical advective flux ! interior values DO jk = 1, jpkm1 zstep = z2 * rdttra(jk) DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. zew = pwn(ji,jj,jk+1) z0w = SIGN( 0.5, pwn(ji,jj,jk+1) ) zalpha = 0.5 + z0w zw = z0w - 0.5 * pwn(ji,jj,jk+1)*zstep / fse3w(ji,jj,jk+1) zzt1 = tb(ji,jj,jk+1) + zw*ztp1(ji,jj,jk+1) zzt2 = tb(ji,jj,jk ) + zw*ztp1(ji,jj,jk ) zzs1 = sb(ji,jj,jk+1) + zw*zsp1(ji,jj,jk+1) zzs2 = sb(ji,jj,jk ) + zw*zsp1(ji,jj,jk ) zt1(ji,jj,jk+1) = zew * ( zalpha * zzt1 + (1.-zalpha)*zzt2 ) zs1(ji,jj,jk+1) = zew * ( zalpha * zzs1 + (1.-zalpha)*zzs2 ) END DO END DO END DO ! surface values IF( lk_dynspg_rl .OR. lk_vvl) THEN ! rigid lid or variable volume: flux set to zero zt1(:,:, 1 ) = 0.e0 zs1(:,:, 1 ) = 0.e0 ELSE ! free surface zt1(:,:, 1 ) = pwn(:,:,1) * tb(:,:,1) zs1(:,:, 1 ) = pwn(:,:,1) * sb(:,:,1) ENDIF ! bottom values zt1(:,:,jpk) = 0.e0 zs1(:,:,jpk) = 0.e0 ! Compute & add the vertical advective trend DO jk = 1, jpkm1 DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. zbtr = 1. / fse3t(ji,jj,jk) ! horizontal advective trends zta = - zbtr * ( zt1(ji,jj,jk) - zt1(ji,jj,jk+1) ) zsa = - zbtr * ( zs1(ji,jj,jk) - zs1(ji,jj,jk+1) ) ! 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 ! Save the vertical advective trends for diagnostic ! ------------------------------------------------- IF( l_trdtra ) THEN ! Recompute the vertical 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 substract ! the term tn()/sn()*hdivn() to recover the W gradz(T/S) trends DO jk = 1, jpkm1 DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. #if defined key_zco zbtr = 1. / ( e1t(ji,jj) * e2t(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 = 1. / ( e1t(ji,jj) * e2t(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 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) ! ENDIF IF(ln_ctl) CALL prt_ctl( tab3d_1=ta, clinfo1=' muscl zad - Ta: ', mask1=tmask , & & tab3d_2=sa, clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' ) ! END SUBROUTINE tra_adv_muscl !!====================================================================== END MODULE traadv_muscl