MODULE dynadv_ubs !!====================================================================== !! *** MODULE dynadv_ubs *** !! Ocean dynamics: Update the momentum trend with the flux form advection !! trend using a 3rd order upstream biased scheme !!====================================================================== !! History : 9.0 ! 06-08 (R. Benshila, L. Debreu) Original code !!---------------------------------------------------------------------- !!---------------------------------------------------------------------- !! dyn_adv_ubs : flux form momentum advection using (ln_dynadv=T) !! an 3rd order Upstream Biased Scheme or Quick scheme !! combined with 2nd or 4th order finite differences !!---------------------------------------------------------------------- USE oce ! ocean dynamics and tracers USE dom_oce ! ocean space and time domain USE dynspg_oce ! surface pressure gradient USE in_out_manager ! I/O manager USE dynspg_rl ! I/O manager USE lbclnk ! ocean lateral boundary conditions (or mpp link) IMPLICIT NONE PRIVATE REAL(wp), PARAMETER :: gamma1 = 1._wp/4._wp ! =1/4 quick ; =1/3 3rd order UBS REAL(wp), PARAMETER :: gamma2 = 1._wp/8._wp ! =0 2nd order ; =1/8 4th order centred !! * Routine accessibility PUBLIC dyn_adv_ubs ! routine called by step.F90 !! * Substitutions # include "domzgr_substitute.h90" # include "vectopt_loop_substitute.h90" !!---------------------------------------------------------------------- !! OPA 9.0 , LODYC-IPSL (2006) !! $Header$ !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) !!---------------------------------------------------------------------- CONTAINS SUBROUTINE dyn_adv_ubs( kt ) !!---------------------------------------------------------------------- !! *** ROUTINE dyn_adv_ubs *** !! !! ** Purpose : Compute the now momentum advection trend in flux form !! and the general trend of the momentum equation. !! !! ** Method : The scheme is the one implemeted in ROMS. It depends !! on two parameter gamma1 and gamma2. The former control the !! upstream baised part of the scheme and the later the centred !! part: gamma1 = 0 pure centered (no diffusive part) !! = 1/4 Quick scheme !! = 1/3 3rd order Upstream biased scheme !! gamma2 = 0 2nd order finite differencing !! = 1/8 4th order finite differencing !! For stability reasons, the first term of the fluxes which cor- !! responds 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). !! Default value (hard coded in the begining of the module) are !! gamma1=1/4 and gamma2=1/8. !! !! ** Action : - update (ua,va) with the 3D advective momentum trends !! !! Reference : Shchepetkin & McWilliams, 2005, Ocean Modelling. !!---------------------------------------------------------------------- USE oce, ONLY: zfu => ta, & ! use ta as 3D workspace zfv => sa ! use sa as 3D workspace INTEGER, INTENT(in) :: kt ! ocean time-step index INTEGER :: ji, jj, jk ! dummy loop indices REAL(wp) :: zua, zva, zbu, zbv ! temporary scalars REAL(wp) :: zui, zvj, zfuj, zfvi, zl_u, zl_v! temporary scalars REAL(wp), DIMENSION(jpi,jpj,jpk) :: zfu_t, zfu_f ! temporary workspace REAL(wp), DIMENSION(jpi,jpj,jpk) :: zfv_t, zfv_f ! " " REAL(wp), DIMENSION(jpi,jpj,jpk) :: zfw, zfu_uw, zfv_vw REAL(wp), DIMENSION(jpi,jpj,jpk,2) :: zlu_uu, zlu_uv ! temporary workspace REAL(wp), DIMENSION(jpi,jpj,jpk,2) :: zlv_vv, zlv_vu ! temporary workspace !!---------------------------------------------------------------------- IF( kt == nit000 ) THEN IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) 'dyn_adv_ubs : UBS flux form momentum advection' IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' ENDIF zfu_t(:,:,:) = 0.e0 zfv_t(:,:,:) = 0.e0 zfu_f(:,:,:) = 0.e0 zfv_f(:,:,:) = 0.e0 zlu_uu(:,:,:,:) = 0.e0 zlv_vv(:,:,:,:) = 0.e0 zlu_uv(:,:,:,:) = 0.e0 zlv_vu(:,:,:,:) = 0.e0 ! ! =============== DO jk = 1, jpkm1 ! Horizontal slab ! ! =============== ! Laplacian ! --------- zfu(:,:,jk) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk) zfv(:,:,jk) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk) DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. zlu_uu(ji,jj,jk,1) = ( ub (ji+1,jj,jk)-2.*ub (ji,jj,jk)+ub (ji-1,jj,jk) ) * umask(ji,jj,jk) zlv_vv(ji,jj,jk,1) = ( vb (ji,jj+1,jk)-2.*vb (ji,jj,jk)+vb (ji,jj-1,jk) ) * vmask(ji,jj,jk) zlu_uv(ji,jj,jk,1) = ( ub (ji,jj+1,jk)-2.*ub (ji,jj,jk)+ub (ji,jj-1,jk) ) * umask(ji,jj,jk) zlv_vu(ji,jj,jk,1) = ( vb (ji+1,jj,jk)-2.*vb (ji,jj,jk)+vb (ji-1,jj,jk) ) * vmask(ji,jj,jk) zlu_uu(ji,jj,jk,2) = ( zfu(ji+1,jj,jk)-2.*zfu(ji,jj,jk)+zfu(ji-1,jj,jk) ) * umask(ji,jj,jk) zlv_vv(ji,jj,jk,2) = ( zfv(ji,jj+1,jk)-2.*zfv(ji,jj,jk)+zfv(ji,jj-1,jk) ) * vmask(ji,jj,jk) zlu_uv(ji,jj,jk,2) = ( zfu(ji,jj+1,jk)-2.*zfu(ji,jj,jk)+zfu(ji,jj-1,jk) ) * umask(ji,jj,jk) zlv_vu(ji,jj,jk,2) = ( zfv(ji+1,jj,jk)-2.*zfv(ji,jj,jk)+zfv(ji-1,jj,jk) ) * vmask(ji,jj,jk) END DO END DO ! ! =============== END DO ! End of slab ! ! =============== CALL lbc_lnk( zlu_uu(:,:,:,1), 'U', -1.) CALL lbc_lnk( zlu_uu(:,:,:,2), 'U', -1.) CALL lbc_lnk( zlv_vv(:,:,:,1), 'V', -1.) CALL lbc_lnk( zlv_vv(:,:,:,2), 'V', -1.) CALL lbc_lnk( zlu_uv(:,:,:,1), 'U', -1.) CALL lbc_lnk( zlu_uv(:,:,:,2), 'U', -1.) CALL lbc_lnk( zlv_vu(:,:,:,1), 'V', -1.) CALL lbc_lnk( zlv_vu(:,:,:,2), 'V', -1.) ! I. Horizontal advection ! ----------------------- ! ! =============== DO jk = 1, jpkm1 ! Horizontal slab ! ! =============== ! horizontal volume fluxes zfu(:,:,jk) = 0.25 * e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk) zfv(:,:,jk) = 0.25 * e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk) ! horizontal momentum fluxes at T- and F-point DO jj = 1, jpjm1 DO ji = 1, fs_jpim1 ! vector opt. zui = ( un(ji,jj,jk) + un(ji+1,jj ,jk) ) zvj = ( vn(ji,jj,jk) + vn(ji ,jj+1,jk) ) IF (zui > 0) THEN ; zl_u = zlu_uu(ji ,jj,jk,1) ELSE ; zl_u = zlu_uu(ji+1,jj,jk,1) ENDIF IF (zvj > 0) THEN ; zl_v = zlv_vv(ji,jj ,jk,1) ELSE ; zl_v = zlv_vv(ji,jj+1,jk,1) ENDIF zfu_t(ji+1,jj ,jk) = ( zfu(ji,jj,jk) + zfu(ji+1,jj ,jk) & & - gamma2 * ( zlu_uu(ji,jj,jk,2) + zlu_uu(ji+1,jj ,jk,2) ) ) & & * ( zui - gamma1 * zl_u) zfv_t(ji ,jj+1,jk) = ( zfv(ji,jj,jk) + zfv(ji ,jj+1,jk) & & - gamma2 * ( zlv_vv(ji,jj,jk,2) + zlv_vv(ji ,jj+1,jk,2) ) ) & & * ( zvj - gamma1 * zl_v) zfuj = ( zfu(ji,jj,jk) + zfu(ji ,jj+1,jk) ) zfvi = ( zfv(ji,jj,jk) + zfv(ji+1,jj ,jk) ) IF (zfuj > 0) THEN ; zl_v = zlv_vu( ji ,jj ,jk,1) ELSE ; zl_v = zlv_vu( ji+1,jj,jk,1) ENDIF IF (zfvi > 0) THEN ; zl_u = zlu_uv( ji,jj ,jk,1) ELSE ; zl_u = zlu_uv( ji,jj+1,jk,1) ENDIF zfv_f(ji ,jj ,jk) = ( zfvi - gamma2 * ( zlv_vu(ji,jj,jk,2) + zlv_vu(ji+1,jj ,jk,2) ) ) & & * ( un(ji,jj,jk) + un(ji ,jj+1,jk) - gamma1 * zl_u ) zfu_f(ji ,jj ,jk) = ( zfuj - gamma2 * ( zlu_uv(ji,jj,jk,2) + zlu_uv(ji ,jj+1,jk,2) ) ) & & * ( vn(ji,jj,jk) + vn(ji+1,jj ,jk) - gamma1 * zl_v ) END DO END DO ! divergence of horizontal momentum fluxes DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. zbu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) zbv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) ! horizontal advective trends zua = - ( zfu_t(ji+1,jj ,jk) - zfu_t(ji ,jj ,jk) & & + zfv_f(ji ,jj ,jk) - zfv_f(ji ,jj-1,jk) ) / zbu zva = - ( zfu_f(ji ,jj ,jk) - zfu_f(ji-1,jj ,jk) & & + zfv_t(ji ,jj+1,jk) - zfv_t(ji ,jj ,jk) ) / zbv ! add it to the general tracer trends ua(ji,jj,jk) = ua(ji,jj,jk) + zua va(ji,jj,jk) = va(ji,jj,jk) + zva #if defined key_trddyn utrd(ji,jj,jk,1) = zua ! save the horizontal advective trend of momentum vtrd(ji,jj,jk,1) = zva #endif END DO END DO ! ! =============== END DO ! End of slab ! ! =============== ! II. Vertical advection ! ---------------------- ! Second order centered tracer flux at w-point DO jk = 1, jpkm1 ! Vertical volume fluxes zfw(:,:,jk) = 0.25 * e1t(:,:) * e2t(:,:) * wn(:,:,jk) ! Vertical advective fluxes IF( jk == 1 ) THEN zfu_uw(:,:,jpk) = 0.e0 ! Bottom value : flux set to zero zfv_vw(:,:,jpk) = 0.e0 ! ! Surface value IF( lk_dynspg_rl ) THEN ! rigid lid : flux set to zero zfu_uw(:,:, 1 ) = 0.e0 zfv_vw(:,:, 1 ) = 0.e0 ELSE ! free surface-constant volume DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 zfu_uw(ji,jj, 1 ) = 2.e0 * ( zfw(ji,jj,1) + zfw(ji+1,jj ,1) ) * un(ji,jj,1) zfv_vw(ji,jj, 1 ) = 2.e0 * ( zfw(ji,jj,1) + zfw(ji ,jj+1,1) ) * vn(ji,jj,1) END DO END DO ENDIF ELSE ! ! interior fluxes DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. zfu_uw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji+1,jj ,jk) ) * ( un(ji,jj,jk) + un(ji,jj,jk-1) ) zfv_vw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji ,jj+1,jk) ) * ( vn(ji,jj,jk) + vn(ji,jj,jk-1) ) END DO END DO ENDIF END DO ! momentum flux divergence at u-, v-points added to the general trend DO jk = 1, jpkm1 DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. zua = - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) ) & & / ( e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) ) zva = - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) ) & & / ( e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) ) ! add it to the general tracer trends ua(ji,jj,jk) = ua(ji,jj,jk) + zua va(ji,jj,jk) = va(ji,jj,jk) + zva END DO END DO END DO END SUBROUTINE dyn_adv_ubs !!============================================================================== END MODULE dynadv_ubs