MODULE traadv_cen2 !!====================================================================== !! *** MODULE traadv_cen2 *** !! Ocean active tracers: horizontal & vertical advective trend !!====================================================================== !! History : 8.2 ! 2001-08 (G. Madec, E. Durand) trahad+trazad=traadv !! 1.0 ! 2002-06 (G. Madec) F90: Free form and module !! 9.0 ! 2004-08 (C. Talandier) New trends organization !! - ! 2005-11 (V. Garnier) Surface pressure gradient organization !! 2.0 ! 2006-04 (R. Benshila, G. Madec) Step reorganization !! - ! 2006-07 (G. madec) add ups_orca_set routine !! 3.2 ! 2009-07 (G. Madec) add avmb, avtb in restart for cen2 advection !!---------------------------------------------------------------------- !!---------------------------------------------------------------------- !! tra_adv_cen2 : update the tracer trend with the horizontal and !! vertical advection trends using a seconder order !! ups_orca_set : allow mixed upstream/centered scheme in specific !! area (set for orca 2 and 4 only) !!---------------------------------------------------------------------- USE oce ! ocean dynamics and active tracers USE dom_oce ! ocean space and time domain USE sbc_oce ! surface boundary condition: ocean USE dynspg_oce ! choice/control of key cpp for surface pressure gradient USE trdmod_oce ! ocean variables trends USE eosbn2 ! equation of state USE trdmod ! ocean active tracers trends USE closea ! closed sea USE trabbl ! advective term in the BBL USE sbcmod ! surface Boundary Condition USE sbcrnf ! river runoffs USE in_out_manager ! I/O manager USE iom ! IOM library USE lib_mpp USE lbclnk ! ocean lateral boundary condition (or mpp link) USE diaptr ! poleward transport diagnostics USE prtctl ! Print control USE zdf_oce ! ocean vertical physics USE restart ! ocean restart IMPLICIT NONE PRIVATE PUBLIC tra_adv_cen2 ! routine called by step.F90 PUBLIC ups_orca_set ! routine used by traadv_cen2_jki.F90 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: upsmsk !: mixed upstream/centered scheme near some straits ! ! and in closed seas (orca 2 and 4 configurations) REAL(wp), DIMENSION(jpi,jpj) :: btr2 ! inverse of T-point surface [1/(e1t*e2t)] !! * Substitutions # include "domzgr_substitute.h90" # include "vectopt_loop_substitute.h90" !!---------------------------------------------------------------------- !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009) !! $Id$ !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) !!---------------------------------------------------------------------- CONTAINS SUBROUTINE tra_adv_cen2( kt, pun, pvn, pwn ) !!---------------------------------------------------------------------- !! *** ROUTINE tra_adv_cen2 *** !! !! ** 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 : !! variable volume (lk_vvl = T) : zero advective flux !! lin. free-surf (lk_vvl = F) : 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, zhw, ze3tr ! temporary scalars REAL(wp) :: zfp_ui, zfp_vj, zfp_w , zfui ! - - REAL(wp) :: zfm_ui, zfm_vj, zfm_w , zfvj ! - - REAL(wp) :: zcofi , zcofj , zcofk ! - - REAL(wp) :: zupsut, zupsus, zcenut, zcenus ! - - REAL(wp) :: zupsvt, zupsvs, zcenvt, zcenvs ! - - REAL(wp) :: zupst , zupss , zcent , zcens ! - - REAL(wp) :: z_hdivn_x, z_hdivn_y, z_hdivn ! - - REAL(wp) :: zice ! - - REAL(wp), DIMENSION(jpi,jpj) :: ztfreez ! 2D workspace 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,*) '~~~~~~~~~~~~ Vector optimization case' IF(lwp) WRITE(numout,*) ! upsmsk(:,:) = 0.e0 ! not upstream by default ! IF( cp_cfg == "orca" ) CALL ups_orca_set ! set mixed Upstream/centered scheme near some straits ! ! and in closed seas (orca2 and orca4 only) ! btr2(:,:) = 1. / ( e1t(:,:) * e2t(:,:) ) ! inverse of T-point surface ! IF( jp_cfg == 2 .AND. .NOT. ln_rstart ) THEN ! Increase the background in the surface layers avmb(1) = 10. * avmb(1) ; avtb(1) = 10. * avtb(1) avmb(2) = 10. * avmb(2) ; avtb(2) = 10. * avtb(2) avmb(3) = 5. * avmb(3) ; avtb(3) = 5. * avtb(3) avmb(4) = 2.5 * avmb(4) ; avtb(4) = 2.5 * avtb(4) ENDIF ENDIF ! Upstream / centered scheme indicator ! ------------------------------------ !!gm not strickly exact : the freezing point should be computed at each ocean levels... !!gm not a big deal since cen2 is no more used in global ice-ocean simulations ztfreez(:,:) = tfreez( sn(:,:,1) ) DO jk = 1, jpk DO jj = 1, jpj DO ji = 1, jpi ! ! below ice covered area (if tn < "freezing"+0.1 ) IF( tn(ji,jj,jk) <= ztfreez(ji,jj) + 0.1 ) THEN ; zice = 1.e0 ELSE ; zice = 0.e0 ENDIF zind(ji,jj,jk) = MAX ( & rnfmsk(ji,jj) * rnfmsk_z(jk), & ! near runoff mouths (& closed sea outflows) upsmsk(ji,jj) , & ! some of some straits zice & ! below ice covered area (if tn < "freezing"+0.1 ) & ) * tmask(ji,jj,jk) END DO END DO END DO ! I. Horizontal advection ! ==================== ! DO jk = 1, jpkm1 ! ! 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 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) ! ! 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. zbtr = btr2(ji,jj) / fse3t(ji,jj,jk) ! ta(ji,jj,jk) = ta(ji,jj,jk) - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk) & & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk) ) sa(ji,jj,jk) = sa(ji,jj,jk) - zbtr * ( zww(ji,jj,jk) - zww(ji-1,jj ,jk) & & + zwz(ji,jj,jk) - zwz(ji ,jj-1,jk) ) END DO END DO END DO IF( l_trdtra ) THEN ! Save the i- and j-advective trends for diagnostic (U.gradz(T) 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 with OBC, BDY, cla, eiv, advective bbl 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 ! 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) ! DO jk = 1, jpkm1 ! T/S MERIDIONAL advection trends DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. 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 ! 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) ! ztrdt(:,:,:) = ta(:,:,:) ; ztrds(:,:,:) = sa(:,:,:) ! Save the horizontal up-to-date ta/sa trends ! ENDIF IF( ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN ! "zonal" mean advective heat and salt transport pht_adv(:) = ptr_vj( zwy(:,:,:) ) pst_adv(:) = ptr_vj( zwz(:,:,:) ) 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' ) ! II. Vertical advection ! ================== ! zwx(:,:,jpk) = 0.e0 ; zwy(:,:,jpk) = 0.e0 ! Bottom value : flux set to zero ! IF( lk_vvl ) THEN ! Surface value : zero in variable volume zwx(:,:, 1 ) = 0.e0 ; zwy(:,:, 1 ) = 0.e0 ELSE ! : linear free surface case zwx(:,:, 1 ) = pwn(:,:,1) * tn(:,:,1) zwy(:,:, 1 ) = pwn(:,:,1) * sn(:,:,1) ENDIF ! DO jk = 2, jpk ! Second order centered tracer flux at w-point DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. zcofk = MAX( zind(ji,jj,jk-1), zind(ji,jj,jk) ) ! upstream indicator zhw = 0.5 * pwn(ji,jj,jk) ! velocity * 1/2 ! zfp_w = zhw + ABS( zhw ) ! upstream scheme 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) ! zcent = zhw * ( tn(ji,jj,jk) + tn(ji,jj,jk-1) ) ! centered scheme zcens = zhw * ( sn(ji,jj,jk) + sn(ji,jj,jk-1) ) ! zwx(ji,jj,jk) = zcofk * zupst + (1.-zcofk) * zcent ! mixed centered / upstream scheme zwy(ji,jj,jk) = zcofk * zupss + (1.-zcofk) * zcens END DO END DO END DO ! DO jk = 1, jpkm1 ! divergence of Tracer flux added to the general trend DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. ze3tr = 1. / fse3t(ji,jj,jk) ta(ji,jj,jk) = ta(ji,jj,jk) - ze3tr * ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) ) sa(ji,jj,jk) = sa(ji,jj,jk) - ze3tr * ( zwy(ji,jj,jk) - zwy(ji,jj,jk+1) ) END DO END DO END DO IF( l_trdtra ) THEN ! Save the vertical advective trends for diagnostic (W gradz(T) trends) DO jk = 1, jpkm1 DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. 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) ! 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 ! write avmb, avtb in restart (traadv_cen2 requires a modified avmb, avtb that are ! --------------------------- required in restart file to ensure restartability) ! avmb, avtb will be read in zdfini in restart case as they are used in zdftke, kpp etc... IF( lrst_oce ) THEN CALL iom_rstput( kt, nitrst, numrow, 'avmb', avmb ) CALL iom_rstput( kt, nitrst, numrow, 'avtb', avtb ) 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 SUBROUTINE ups_orca_set !!---------------------------------------------------------------------- !! *** ROUTINE ups_orca_set *** !! !! ** Purpose : add a portion of upstream scheme in area where the !! centered scheme generates too strong overshoot !! !! ** Method : orca (R4 and R2) confiiguration setting. Set upsmsk !! array to nozero value in some straith. !! !! ** Action : - upsmsk set to 1 at some strait, 0 elsewhere for orca !!---------------------------------------------------------------------- INTEGER :: ii0, ii1, ij0, ij1 ! temporary integers !!---------------------------------------------------------------------- ! mixed upstream/centered scheme near river mouths ! ------------------------------------------------ SELECT CASE ( jp_cfg ) ! ! ======================= CASE ( 4 ) ! ORCA_R4 configuration ! ! ======================= ! ! Gibraltar Strait ii0 = 70 ; ii1 = 71 ij0 = 52 ; ij1 = 53 ; upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50 ! ! ! ======================= CASE ( 2 ) ! ORCA_R2 configuration ! ! ======================= ! ! Gibraltar Strait ij0 = 102 ; ij1 = 102 ii0 = 138 ; ii1 = 138 ; upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.20 ii0 = 139 ; ii1 = 139 ; upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.40 ii0 = 140 ; ii1 = 140 ; upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50 ij0 = 101 ; ij1 = 102 ii0 = 141 ; ii1 = 141 ; upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50 ! ! Bab el Mandeb Strait ij0 = 87 ; ij1 = 88 ii0 = 164 ; ii1 = 164 ; upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.10 ij0 = 88 ; ij1 = 88 ii0 = 163 ; ii1 = 163 ; upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.25 ii0 = 162 ; ii1 = 162 ; upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.40 ii0 = 160 ; ii1 = 161 ; upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50 ij0 = 89 ; ij1 = 89 ii0 = 158 ; ii1 = 160 ; upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.25 ij0 = 90 ; ij1 = 90 ii0 = 160 ; ii1 = 160 ; upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.25 ! ! Sound Strait ij0 = 116 ; ij1 = 116 ii0 = 144 ; ii1 = 144 ; upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.25 ii0 = 145 ; ii1 = 147 ; upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50 ii0 = 148 ; ii1 = 148 ; upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.25 ! END SELECT ! mixed upstream/centered scheme over closed seas ! ----------------------------------------------- CALL clo_ups( upsmsk(:,:) ) ! END SUBROUTINE ups_orca_set !!====================================================================== END MODULE traadv_cen2