MODULE traadv_cen2_jki !!============================================================================== !! *** MODULE traadv_cen2_jki *** !! Ocean active tracers: horizontal & vertical advective trend !!============================================================================== !! History : 8.2 ! 01-08 (G. Madec, E. Durand) trahad+trazad = traadv !! 8.5 ! 02-06 (G. Madec) F90: Free form and module !! 9.0 ! 05-11 (V. Garnier) Surface pressure gradient organization !! " " ! 06-04 (R. Benshila, G. Madec) Step reorganization !!---------------------------------------------------------------------- !!---------------------------------------------------------------------- !! tra_adv_cen2_jki : update the tracer trend with the horizontal and !! vertical advection trends using a seconder order !! centered scheme. Auto-tasking case, k-slab for !! hor. adv., j-slab for vert. adv. (j-k-i loops) !!---------------------------------------------------------------------- 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 flxrnf ! USE trabbl ! advective term in the BBL USE ocfzpt ! 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 ! Print control IMPLICIT NONE PRIVATE PUBLIC tra_adv_cen2_jki ! routine called by step.F90 REAL(wp), DIMENSION(jpi,jpj) :: btr2 ! inverse of T-point surface [1/(e1t*e2t)] !! * Substitutions # include "domzgr_substitute.h90" # include "vectopt_loop_substitute.h90" !!---------------------------------------------------------------------- !! OPA 9.0 , LOCEAN-IPSL (2005) !! $Id$ !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) !!---------------------------------------------------------------------- CONTAINS SUBROUTINE tra_adv_cen2_jki( kt, pun, pvn, pwn ) !!---------------------------------------------------------------------- !! *** ROUTINE tra_adv_cen2_jki *** !! !! ** 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 second order centered !! scheme using now fields (leap-frog scheme). In specific areas !! (vicinity of major river mouths, some straits, or where tn is !! approaching the freezing point) it is mixed with an upstream !! scheme for stability reasons. !! Part 0 : compute the upstream / centered flag !! (3D array, zind, defined at T-point (00 or <0] !! zupsv = e1v*e3v vn (tb(j) or tb(j-1) ) [vn>0 or <0] !! * mixed upstream / centered horizontal advection scheme !! zcofi = max(zind(i+1), zind(i)) !! zcofj = max(zind(j+1), zind(j)) !! zwx = zcofi * zupsu + (1-zcofi) * zcenu !! zwy = zcofj * zupsv + (1-zcofj) * zcenv !! * horizontal advective trend (divergence of the fluxes) !! zta = 1/(e1t*e2t*e3t) { di-1[zwx] + dj-1[zwy] } !! * Add this trend now to the general trend of tracer (ta,sa): !! (ta,sa) = (ta,sa) + ( zta , zsa ) !! * trend diagnostic ('key_trdtra' defined): the trend is !! saved for diagnostics. The trends saved is expressed as !! Uh.gradh(T), i.e. !! save trend = zta + tn divn !! In addition, the advective trend in the two horizontal direc- !! tion is also re-computed as Uh gradh(T). Indeed hadt+tn divn is !! equal to (in s-coordinates, and similarly in z-coord.): !! zta+tn*divn=1/(e1t*e2t*e3t) { mi-1( e2u*e3u un di[tn] ) !! +mj-1( e1v*e3v vn mj[tn] ) } !! NB:in z-coordinate - full step (ln_zco=T) e3u=e3v=e3t, so !! they vanish from the expression of the flux and divergence. !! !! Part II : vertical advection !! For temperature (idem for salinity) the advective trend is com- !! puted as follows : !! zta = 1/e3t dk+1[ zwz ] !! where the vertical advective flux, zwz, is given by : !! zwz = zcofk * zupst + (1-zcofk) * zcent !! with !! zupsv = upstream flux = wn * (tb(k) or tb(k-1) ) [wn>0 or <0] !! zcenu = centered flux = wn * mk(tn) !! The surface boundary condition is : !! rigid-lid (lk_dynspg_frd = T) : zero advective flux !! free-surf (lk_dynspg_fsc = T) : wn(:,:,1) * tn(:,:,1) !! Add this trend now to the general trend of tracer (ta,sa): !! (ta,sa) = (ta,sa) + ( zta , zsa ) !! Trend diagnostic ('key_trdtra' defined): the trend is !! saved for diagnostics. The trends saved is expressed as : !! save trend = w.gradz(T) = zta - tn divn. !! !! ** Action : - update (ta,sa) with the now advective tracer trends !! - save trends in (ztrdt,ztrds) ('key_trdtra') !!---------------------------------------------------------------------- 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 ! 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) :: & zbtr, zta, zsa, zfui, zfvj, & ! temporary scalars zhw, ze3tr, zcofi, zcofj, & ! " " zupsut, zupsvt, zupsus, zupsvs, & ! " " zfp_ui, zfp_vj, zfm_ui, zfm_vj, & ! " " zcofk, zupst, zupss, zcent, & ! " " zcens, zfp_w, zfm_w, & ! " " zcenut, zcenvt, zcenus, zcenvs, & ! " " z_hdivn_x, z_hdivn_y, z_hdivn REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwz, ztrdt, zind ! 3D workspace REAL(wp), DIMENSION(jpi,jpj,jpk) :: zww, ztrds ! " " !!---------------------------------------------------------------------- IF( kt == nit000 ) THEN IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) 'tra_adv_cen2 : 2nd order centered advection scheme' IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~ Auto-tasking case' IF(lwp) WRITE(numout,*) ! btr2(:,:) = 1. / ( e1t(:,:) * e2t(:,:) ) ENDIF !CDIR PARALLEL DO !$OMP PARALLEL DO ! ! =============== DO jk = 1, jpkm1 ! Horizontal slab ! ! =============== ! Upstream / centered scheme indicator ! -------------------------------------- DO jj = 1, jpj DO ji = 1, jpi zind(ji,jj,jk) = MAX ( & upsrnfh(ji,jj) * upsrnfz(jk), & ! changing advection scheme near runoff upsadv(ji,jj) & ! in the vicinity of some straits #if defined key_ice_lim , tmask(ji,jj,jk) & ! half upstream tracer fluxes if tn < ("freezing"+0.1 ) * MAX( 0., SIGN( 1., fzptn(ji,jj)+0.1-tn(ji,jj,jk) ) ) & #endif ) END DO END DO ! Horizontal advective fluxes ! ----------------------------- ! Second order centered tracer flux at u and v-points DO jj = 1, jpjm1 DO ji = 1, fs_jpim1 ! vector opt. ! upstream indicator zcofi = MAX( zind(ji+1,jj,jk), zind(ji,jj,jk) ) zcofj = MAX( zind(ji,jj+1,jk), zind(ji,jj,jk) ) ! 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 ) zupsut = zfp_ui * tb(ji,jj,jk) + zfm_ui * tb(ji+1,jj ,jk) zupsvt = zfp_vj * tb(ji,jj,jk) + zfm_vj * tb(ji ,jj+1,jk) zupsus = zfp_ui * sb(ji,jj,jk) + zfm_ui * sb(ji+1,jj ,jk) zupsvs = zfp_vj * sb(ji,jj,jk) + zfm_vj * sb(ji ,jj+1,jk) ! 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) = zcofi * zupsut + (1.-zcofi) * zcenut zwy(ji,jj,jk) = zcofj * zupsvt + (1.-zcofj) * zcenvt zww(ji,jj,jk) = zcofi * zupsus + (1.-zcofi) * zcenus zwz(ji,jj,jk) = zcofj * zupsvs + (1.-zcofj) * zcenvs 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. #if defined key_zco zbtr = btr2(ji,jj) #else zbtr = btr2(ji,jj) / fse3t(ji,jj,jk) #endif ! horizontal advective trends 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 ! ! =============== IF(ln_ctl) CALL prt_ctl( tab3d_1=ta, clinfo1=' cen2 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 = btr2(ji,jj) z_hdivn_x = ( e2u(ji ,jj) * pun(ji ,jj,jk) & & - e2u(ji-1,jj) * pun(ji-1,jj,jk) ) * zbtr #else zbtr = btr2(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 * ( zwx(ji,jj,jk) - zwx(ji-1,jj,jk) ) + tn(ji,jj,jk) * z_hdivn_x ztrds(ji,jj,jk) = - zbtr * ( zww(ji,jj,jk) - zww(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 = btr2(ji,jj) z_hdivn_y = ( e1v(ji,jj ) * pvn(ji,jj ,jk) & & - e1v(ji,jj-1) * pvn(ji,jj-1,jk) ) * zbtr #else zbtr = btr2(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 * ( zwy(ji,jj,jk) - zwy(ji,jj-1,jk) ) + tn(ji,jj,jk) * z_hdivn_y ztrds(ji,jj,jk) = - zbtr * ( zwz(ji,jj,jk) - zwz(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 IF(ln_ctl) CALL prt_ctl( tab3d_1=ta, clinfo1=' cen2 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 ! Vertical advection ! -------------------- !CDIR PARALLEL DO !$OMP PARALLEL DO ! ! =============== DO jj = 2, jpjm1 ! Vertical slab ! ! =============== ! Bottom value : flux and indicator set to zero zwz (:,jj,jpk) = 0.e0 ; zww(:,jj,jpk) = 0.e0 zind(:,jj,jpk) = 0.e0 ! Surface value IF( lk_dynspg_rl ) THEN ! rigid lid : flux set to zero zwz(:,jj, 1 ) = 0.e0 ; zww(:,jj, 1 ) = 0.e0 ELSE ! free surface-constant volume zwz(:,jj, 1 ) = pwn(:,jj,1) * tn(:,jj,1) zww(:,jj, 1 ) = pwn(:,jj,1) * sn(:,jj,1) ENDIF ! Vertical advective fluxes at w-point DO jk = 2, jpk DO ji = 2, jpim1 ! upstream indicator zcofk = MAX( zind(ji,jj,jk-1), zind(ji,jj,jk) ) ! velocity * 1/2 zhw = 0.5 * pwn(ji,jj,jk) ! upstream scheme zfp_w = zhw + ABS( zhw ) zfm_w = zhw - ABS( zhw ) zupst = zfp_w * tb(ji,jj,jk) + zfm_w * tb(ji,jj,jk-1) zupss = zfp_w * sb(ji,jj,jk) + zfm_w * sb(ji,jj,jk-1) ! centered scheme zcent = zhw * ( tn(ji,jj,jk) + tn(ji,jj,jk-1) ) zcens = zhw * ( sn(ji,jj,jk) + sn(ji,jj,jk-1) ) ! mixed centered / upstream scheme zwz(ji,jj,jk) = zcofk * zupst + (1.-zcofk) * zcent zww(ji,jj,jk) = zcofk * zupss + (1.-zcofk) * zcens END DO END DO ! Tracer flux divergence at t-point added to the general trend DO jk = 1, jpkm1 DO ji = 2, jpim1 ze3tr = 1. / fse3t(ji,jj,jk) ! vertical advective trends zta = - ze3tr * ( zwz(ji,jj,jk) - zwz(ji,jj,jk+1) ) zsa = - ze3tr * ( zww(ji,jj,jk) - zww(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 ! End of slab ! ! =============== ! 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 = btr2(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 = btr2(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=' cen2 zad - Ta: ', mask1=tmask, & & tab3d_2=sa, clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' ) ! END SUBROUTINE tra_adv_cen2_jki !!====================================================================== END MODULE traadv_cen2_jki