MODULE traadv_qck !!============================================================================== !! *** MODULE traadv_qck *** !! Ocean active tracers: horizontal & vertical advective trend !!============================================================================== !!---------------------------------------------------------------------- !! tra_adv_qck : update the tracer trend with the horizontal !! advection trends using a 3st order !! finite difference scheme !! The vertical advection scheme is the 2nd centered scheme !!---------------------------------------------------------------------- !! * Modules used USE oce ! ocean dynamics and active tracers USE dom_oce ! ocean space and time domain USE dynspg_oce ! USE trdmod_oce ! ocean variables trends USE trdmod ! ocean active tracers trends USE trabbl ! advective term in the BBL 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 prtctl ! Print control IMPLICIT NONE PRIVATE !! * Accessibility PUBLIC tra_adv_qck ! routine called by step.F90 !! * Module variables REAL(wp), DIMENSION(jpi,jpj), SAVE :: & zbtr2 REAL(wp), DIMENSION(jpi,jpj,jpk), SAVE :: & sl REAL(wp) :: & cst1, cst2, dt, coef1 ! temporary scalars INTEGER :: & ji, jj, jk ! dummy loop indices !!---------------------------------------------------------------------- !! * Substitutions # include "domzgr_substitute.h90" # include "vectopt_loop_substitute.h90" !!---------------------------------------------------------------------- !! OPA 9.0 , LOCEAN-IPSL (2005) !! $Id$ !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt !!---------------------------------------------------------------------- CONTAINS SUBROUTINE tra_adv_qck( kt, pun, pvn, pwn ) !!---------------------------------------------------------------------- !! *** ROUTINE tra_adv_qck *** !! !! ** Purpose : Compute the now trend due to the advection of tracers !! and add it to the general trend of passive tracer equations. !! !! ** Method : The advection is evaluated by a third order scheme !! For a positive velocity u : !! !! !! i-1 i i+1 i+2 !! !! |--FU--|--FC--|--FD--|------| !! u(i)>0 !! !! For a negative velocity u : !! !! |------|--FD--|--FC--|--FU--| !! u(i)<0 !! !! FU is the second upwind point !! FD is the first douwning point !! FC is the central point (or the first upwind point) !! !! Flux(i) = u(i) * {0.5(FC+FD) -0.5C(i)(FD-FC) -((1-C(i)å?)/6)(FU+FD-2FC)} !! with C(i)=|u(i)|dx(i)/dt (Courant number) !! !! dt = 2*rdtra and the scalar values are tb and sb !! !! On the vertical, the simple centered scheme used tn and sn !! !! The fluxes are bounded by the ULTIMATE limiter !! to guarantee the monotonicity of the solution and to !! prevent the appearance of spurious numerical oscillations !! !! ** Action : - update (ta,sa) with the now advective tracer trends !! - save the trends in (ttrdh,strdh) ('key_trdtra') !! !! ** Reference : Leonard (1979, 1991) !! History : !! 9.0 ! 06-09 (G. Reffray) Original code !!---------------------------------------------------------------------- !! * Arguments 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 !! REAL(wp) :: z2 ! temporary scalar REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztrdt, ztrds ! temporary 3D workspace IF( kt == nit000 ) THEN IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) 'tra_adv_qck : 3st order quickest advection scheme' IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~ Vector optimization case' IF(lwp) WRITE(numout,*) zbtr2(:,:) = 1. / ( e1t(:,:) * e2t(:,:) ) cst1 = 1./12. cst2 = 2./3. IF (l_trdtra ) THEN CALL ctl_warn( ' Trends not yet implemented for PPM advection scheme ' ) ENDIF ENDIF IF( neuler == 0 .AND. kt == nit000 ) THEN ; z2 = 1. ELSE ; z2 = 2. ENDIF ! Save ta and sa trends IF( l_trdtra ) THEN ! to be done ztrdt(:,:,:) = ta(:,:,:) ztrds(:,:,:) = sa(:,:,:) l_adv = 'qst' ENDIF ! I. Slope estimation at the T-point for the limiter ULTIMATE ! SL = Sum(1/C_out) with C_out : Courant number for the outflows !--------------------------------------------------------------- sl(:,:,:) = 100. ! =============== DO jk = 1, jpkm1 ! Horizontal slab ! ! =============== dt = z2 * rdttra(jk) DO jj = 2, jpjm1 DO ji = 2, fs_jpim1 ! vector opt. coef1 = 1.e-12 IF (pun(ji-1,jj ,jk ).LT.0.) coef1 = coef1 + ABS(pun(ji-1,jj ,jk ))*dt/e1t(ji,jj) IF (pun(ji ,jj ,jk ).GT.0.) coef1 = coef1 + ABS(pun(ji ,jj ,jk ))*dt/e1t(ji,jj) IF (pvn(ji ,jj-1,jk ).LT.0.) coef1 = coef1 + ABS(pvn(ji ,jj-1,jk ))*dt/e2t(ji,jj) IF (pvn(ji ,jj ,jk ).GT.0.) coef1 = coef1 + ABS(pvn(ji ,jj ,jk ))*dt/e2t(ji,jj) IF (pwn(ji ,jj ,jk+1).LT.0.) coef1 = coef1 + ABS(pwn(ji ,jj ,jk+1))*dt/(fse3t(ji,jj,jk)) IF (pwn(ji ,jj ,jk ).GT.0.) coef1 = coef1 + ABS(pwn(ji ,jj ,jk ))*dt/(fse3t(ji,jj,jk)) sl(ji,jj,jk) = 1./coef1 sl(ji,jj,jk) = MIN(100.,sl(ji,jj,jk)) sl(ji,jj,jk) = MAX(1. ,sl(ji,jj,jk)) ENDDO ENDDO ENDDO !--- Lateral boundary conditions on the limiter slope CALL lbc_lnk( sl(:,:,:), 'T', 1. ) ! II. The horizontal fluxes are computed with the QUICKEST + ULTIMATE scheme !--------------------------------------------------------------------------- CALL tra_adv_qck_hor( kt , pun, pvn, tb , ta , pht_adv , z2) CALL tra_adv_qck_hor( kt , pun, pvn, sb , sa , pst_adv , z2) ! Save the horizontal advective trends for diagnostic ! --------------------------------------------------- ! IF( l_trdtra ) THEN ! to be done ! ! T/S ZONAL advection trends ! ENDIF IF(ln_ctl) THEN CALL prt_ctl(tab3d_1=ta, clinfo1=' centered2 had - Ta: ', mask1=tmask, & & tab3d_2=sa, clinfo2=' Sa: ', mask2=tmask, clinfo3='tra') ENDIF ! III. The vertical fluxes are computed with the 2nd order centered scheme !------------------------------------------------------------------------- CALL tra_adv_qck_ver( pwn, tn , ta, z2 ) CALL tra_adv_qck_ver( pwn, sn , sa, z2 ) ! Save the vertical advective trends for diagnostic ! ------------------------------------------------- ! IF( l_trdtra ) THEN ! to be done ! 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 - ztdta()/ztdsa() and substract ! the term tn()/sn()*hdivn() to recover the W gradz(T/S) trends ! ENDIF IF(ln_ctl) THEN CALL prt_ctl(tab3d_1=ta, clinfo1=' centered2 zad - Ta: ', mask1=tmask, & & tab3d_2=sa, clinfo2=' Sa: ', mask2=tmask, clinfo3='tra') ENDIF END SUBROUTINE tra_adv_qck SUBROUTINE tra_adv_qck_hor ( kt , pun, pvn, tra , traa , phtra_adv ,z2 ) !!---------------------------------------------------------------------- !! !!---------------------------------------------------------------------- !! * Arguments INTEGER, INTENT( in ) :: kt ! ocean time-step index REAL, INTENT( in ) :: z2 REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) :: pun, pvn ! horizontal effective velocity REAL(wp), INTENT ( out ), DIMENSION(jpj) :: & phtra_adv REAL(wp), INTENT ( inout ), DIMENSION(jpi,jpj,jpk) :: & tra, traa REAL(wp) :: & za, zbtr, e1, e2, c, dir, fu, fc, fd, & ! temporary scalars coef2, coef3, fho, mask, dx REAL(wp), DIMENSION(jpi,jpj) :: & zee REAL(wp), DIMENSION(jpi,jpj,jpk) :: & zmask, zlap, dwst, lim !---------------------------------------------------------------------- ! 0. Initialization (should ot be needed on the whole array ???) !---------------------------------------------------------------------- zmask = 0.0 zlap = 0.0 dwst = 0.0 lim = 0.0 !---------------------------------------------------------------------- ! I. Part 1 : x-direction !---------------------------------------------------------------------- ! =============== 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 zee(ji,jj) = e2u(ji,jj) / e1u(ji,jj) * umask(ji,jj,jk) #else ! vertical scale factor are used zee(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj) * umask(ji,jj,jk) #endif END DO END DO ! Laplacian of tracers (at before time step) ! ------------------------------------------ !--- First derivative (gradient) DO jj = 1, jpjm1 DO ji = 1, fs_jpim1 ! vector opt. zmask(ji,jj,jk) = zee(ji,jj) * ( tra(ji+1,jj ,jk) - tra(ji,jj,jk) ) END DO END DO DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. #if defined key_zco zee(ji,jj) = e1t(ji,jj) / e2t(ji,jj) #else zee(ji,jj) = e1t(ji,jj) / (e2t(ji,jj) * fse3t(ji,jj,jk)) #endif zlap(ji,jj,jk) = zee(ji,jj) * ( zmask(ji,jj,jk) - zmask(ji-1,jj,jk) ) END DO END DO !--- Function lim=FU+SL*(FC-FU) used by the limiter !--- Computation of the ustream and downstream lim at the T-points DO jj = 2, jpjm1 DO ji = 2, fs_jpim1 ! vector opt. ! Upstream in the x-direction for the tracer zmask(ji,jj,jk)=tra(ji-1,jj,jk)+sl(ji,jj,jk)*(tra(ji,jj,jk)-tra(ji-1,jj,jk)) ! Downstream in the x-direction for the tracer dwst (ji,jj,jk)=tra(ji+1,jj,jk)+sl(ji,jj,jk)*(tra(ji,jj,jk)-tra(ji+1,jj,jk)) ENDDO ENDDO END DO !--- Lateral boundary conditions on the laplacian (unchanged sgn) CALL lbc_lnk( zlap(:,:,:), 'T', 1. ) !--- Lateral boundary conditions for the lim function CALL lbc_lnk( zmask(:,:,:), 'T', 1. ) ; CALL lbc_lnk( dwst(:,:,:), 'T', 1. ) ! =============== DO jk = 1, jpkm1 ! Horizontal slab ! ! =============== ! --- lim at the U-points in function of the direction of the flow ! ---------------------------------------------------------------- DO jj = 1, jpjm1 DO ji = 1, fs_jpim1 ! vector opt. dir = 0.5 + sign(0.5,pun(ji,jj,jk)) ! if pun>0 : diru = 1 otherwise diru = 0 lim(ji,jj,jk)=dir*zmask(ji,jj,jk)+(1-dir)*dwst(ji+1,jj,jk) ! Mask at the T-points in the x-direction (mask=0 or mask=1) zmask(ji,jj,jk)=tmask(ji-1,jj,jk)+tmask(ji,jj,jk)+tmask(ji+1,jj,jk)-2 END DO END DO END DO !--- Lateral boundary conditions for the mask CALL lbc_lnk( zmask(:,:,:), 'T', 1. ) ! Horizontal advective fluxes ! --------------------------- ! =============== DO jk = 1, jpkm1 ! Horizontal slab ! =============== dt = z2 * rdttra(jk) !--- tracer flux at u and v-points DO jj = 1, jpjm1 DO ji = 1, fs_jpim1 ! vector opt. #if defined key_zco e2 = e2u(ji,jj) #else e2 = e2u(ji,jj) * fse3u(ji,jj,jk) #endif dir = 0.5 + sign(0.5,pun(ji,jj,jk)) ! if pun>0 : diru = 1 otherwise diru = 0 dx = dir * e1t(ji,jj) + (1-dir)* e1t(ji+1,jj) c = ABS(pun(ji,jj,jk))*dt/dx ! (00 : dirv = 1 otherwise dirv = 0 lim(ji,jj,jk)=dir*zmask(ji,jj,jk)+(1-dir)*dwst(ji,jj+1,jk) ! Mask at the T-points in the y-direction (mask=0 or mask=1) zmask(ji,jj,jk)=tmask(ji,jj-1,jk)+tmask(ji,jj,jk)+tmask(ji,jj+1,jk)-2 END DO END DO END DO !--- Lateral boundary conditions for the mask CALL lbc_lnk( zmask(:,:,:), 'T', 1. ) ! Horizontal advective fluxes ! ------------------------------- ! =============== DO jk = 1, jpkm1 ! Horizontal slab ! =============== dt = z2 * rdttra(jk) !--- tracer flux at u and v-points DO jj = 1, jpjm1 DO ji = 1, fs_jpim1 ! vector opt. #if defined key_zco e1 = e1v(ji,jj) #else e1 = e1v(ji,jj) * fse3v(ji,jj,jk) #endif dir = 0.5 + sign(0.5,pvn(ji,jj,jk)) ! if pvn>0 : dirv = 1 otherwise dirv = 0 dx = dir * e2t(ji,jj) + (1-dir)* e2t(ji,jj+1) c = ABS(pvn(ji,jj,jk))*dt/dx ! (00 : dirw = 1 otherwise dirw = 0 fc = dir*tra(ji,jj,jk )*fse3t(ji,jj,jk-1)+(1-dir)*tra(ji,jj,jk-1)*fse3t(ji,jj,jk ) ! FC in the z-direction for T fd = dir*tra(ji,jj,jk-1)*fse3t(ji,jj,jk )+(1-dir)*tra(ji,jj,jk )*fse3t(ji,jj,jk-1) ! FD in the z-direction for T !--- Second order centered scheme sl(ji,jj,jk)=pwn(ji,jj,jk)*(fc+fd)/(fse3t(ji,jj,jk-1)+fse3t(ji,jj,jk)) END DO END DO END DO ! 2. 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. ze3tr = 1. / fse3t(ji,jj,jk) ! vertical advective trends za = - ze3tr * ( sl(ji,jj,jk) - sl(ji,jj,jk+1) ) ! add it to the general tracer trends traa(ji,jj,jk) = traa(ji,jj,jk) + za END DO END DO END DO END SUBROUTINE tra_adv_qck_ver REAL FUNCTION bound(fu,fd,fc,fho) real :: fu,fd,fc,fho,fref1,fref2 fref1 = fu fref2 = MAX(MIN(fc,fd),MIN(MAX(fc,fd),fref1)) bound = MAX(MIN(fho,fc),MIN(MAX(fho,fc),fref2)) END FUNCTION !!====================================================================== END MODULE traadv_qck