!!---------------------------------------------------------------------- !! *** trcbbl_adv.h90 *** !!---------------------------------------------------------------------- !!---------------------------------------------------------------------- !! TOP 1.0 , LOCEAN-IPSL (2005) !! $Header$ !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt !!---------------------------------------------------------------------- SUBROUTINE trc_bbl_adv( kt ) !!---------------------------------------------------------------------- !! *** ROUTINE trc_bbl_adv *** !! !! ** Purpose : Compute the before tracer trend associated !! with the bottom boundary layer and add it to the general trend !! of tracer equations. The bottom boundary layer is supposed to be !! both an advective and diffusive bottom boundary layer. !! !! ** Method : Computes the bottom boundary horizontal and vertical !! advection terms. Add it to the general trend : tra =tra + adv_bbl. !! When the product grad( rho) * grad(h) < 0 (where grad is a !! along bottom slope gradient) an additional lateral 2nd order !! diffusion along the bottom slope is added to the general !! tracer trend, otherwise the additional trend is set to 0. !! Second order operator (laplacian type) with variable coefficient !! computed as follow for temperature (idem on s): !! difft = 1/(e1t*e2t*e3t) { di-1[ ahbt e2u*e3u/e1u di[ztb] ] !! + dj-1[ ahbt e1v*e3v/e2v dj[ztb] ] } !! where ztb is a 2D array: the bottom ocean te;perature and ahtb !! is a time and space varying diffusive coefficient defined by: !! ahbt = zahbp if grad(rho).grad(h) < 0 !! = 0. otherwise. !! Note that grad(.) is the along bottom slope gradient. grad(rho) !! is evaluated using the local density (i.e. referenced at the !! local depth). Typical value of ahbt is 2000 m2/s (equivalent to !! a downslope velocity of 20 cm/s if the condition for slope !! convection is satified) !! Add this before trend to the general trend tra of the !! botton ocean tracer point: !! tra = tra + difft !! !! ** Action : - update tra at the bottom level with the bottom !! boundary layer trend !! !! References : !! Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591. !! !! History : !! 8.5 ! 02-12 (A. de Miranda, G. Madec) Original Code !! 9.0 ! 04-01 (A. de Miranda, G. Madec, J.M. Molines ) !! 9.0 ! 04-03 (C. Ethe) Adaptation for Passive tracers !!---------------------------------------------------------------------- !! * Modules used USE lbclnk ! ocean lateral boundary conditions (or mpp link) !! * Arguments INTEGER, INTENT( in ) :: kt ! ocean time-step !! * Local declarations INTEGER :: ji, jj, jk, jn ! dummy loop indices INTEGER :: ik, iku, ikv ! temporary integers REAL(wp) :: & zsign, zt, zs, zh, zalbet, & ! temporary scalars zgdrho, zbtr, ztra ! " " REAL(wp), DIMENSION(jpi,jpj) :: & zki, zkj, zkw, zkx, zky, zkz, & ! temporary workspace arrays ztnb, zsnb, zdep, ztrb, & ! " " zahu, zahv ! " " REAL(wp), DIMENSION(jpi,jpj) :: & ! temporary workspace arrays zalphax, zwu, zunb, & ! " " zalphay, zwv, zvnb, & ! " " zwx, zwy ! " " REAL(wp), DIMENSION(jpi,jpj,jpk) :: & zhdivn ! temporary workspace arrays REAL(wp) :: & zfui, zfvj, zbt, zsigna ! temporary scalars REAL(wp) :: & fsalbt, pft, pfs, pfh ! statement function CHARACTER (len=22) :: charout !!---------------------------------------------------------------------- ! ratio alpha/beta ! ================ ! fsalbt: ratio of thermal over saline expension coefficients ! pft : potential temperature in degrees celcius ! pfs : salinity anomaly (s-35) in psu ! pfh : depth in meters fsalbt( pft, pfs, pfh ) = & ( ( ( -0.255019e-07 * pft + 0.298357e-05 ) * pft & - 0.203814e-03 ) * pft & + 0.170907e-01 ) * pft & + 0.665157e-01 & +(-0.678662e-05 * pfs - 0.846960e-04 * pft + 0.378110e-02 ) * pfs & + ( ( - 0.302285e-13 * pfh & - 0.251520e-11 * pfs & + 0.512857e-12 * pft * pft ) * pfh & - 0.164759e-06 * pfs & +( 0.791325e-08 * pft - 0.933746e-06 ) * pft & + 0.380374e-04 ) * pfh !!---------------------------------------------------------------------- IF( kt == nittrc000 ) CALL trc_bbl_init ! initialization at first time-step ! 1. 2D fields of bottom temperature and salinity, and bottom slope ! ----------------------------------------------------------------- ! mbathy= number of w-level, minimum value=1 (cf dommsk.F) #if defined key_vectopt_loop && ! defined key_autotasking jj = 1 DO ji = 1, jpij ! vector opt. (forced unrolling) #else DO jj = 1, jpj DO ji = 1, jpi #endif ik = mbkt(ji,jj) ! index of the bottom ocean T-level ztnb(ji,jj) = tn(ji,jj,ik) * tmask(ji,jj,1) ! masked now T at the ocean bottom zsnb(ji,jj) = sn(ji,jj,ik) * tmask(ji,jj,1) ! masked now S at the ocean bottom zdep(ji,jj) = fsdept(ji,jj,ik) ! depth of the ocean bottom T-level #if ! defined key_vectopt_loop || defined key_autotasking END DO #endif END DO #if defined key_vectopt_loop && ! defined key_autotasking jj = 1 DO ji = 1, jpij-jpi ! vector opt. (forced unrolling) zunb(ji,jj) = un(ji,jj,mbku(ji,jj)) * umask(ji,jj,1) zvnb(ji,jj) = vn(ji,jj,mbkv(ji,jj)) * vmask(ji,jj,1) ! retirer le mask en u, v et t ! END DO #else DO jj = 1, jpjm1 DO ji = 1, jpim1 zunb(ji,jj) = un(ji,jj,mbku(ji,jj)) * umask(ji,jj,1) zvnb(ji,jj) = vn(ji,jj,mbkv(ji,jj)) * vmask(ji,jj,1) END DO END DO #endif ! boundary conditions on zunb and zvnb (changed sign) CALL lbc_lnk( zunb, 'U', -1. ) ; CALL lbc_lnk( zvnb, 'V', -1. ) ! Conditional diffusion along the slope in the bottom boundary layer #if defined key_trcbbl_dif # if defined key_vectopt_loop && ! defined key_autotasking jj = 1 DO ji = 1, jpij-jpi ! vector opt. (forced unrolling) # else DO jj = 1, jpjm1 DO ji = 1, jpim1 # endif iku = mbku(ji,jj) ikv = mbkv(ji,jj) zahu(ji,jj) = atrbbl*e2u(ji,jj)*fse3u(ji,jj,iku)/e1u(ji,jj) * umask(ji,jj,1) zahv(ji,jj) = atrbbl*e1v(ji,jj)*fse3v(ji,jj,ikv)/e2v(ji,jj) * vmask(ji,jj,1) # if ! defined key_vectopt_loop || defined key_autotasking END DO # endif END DO #endif ! 2. Criteria of additional bottom diffusivity: grad(rho).grad(h)<0 ! -------------------------------------------- ! Sign of the local density gradient along the i- and j-slopes ! multiplied by the slope of the ocean bottom SELECT CASE ( neos ) CASE ( 0 ) ! Jackett and McDougall (1994) formulation DO jj = 1, jpjm1 DO ji = 1, fs_jpim1 ! vector opt. ! ... temperature, salinity anomalie and depth zt = 0.5 * ( ztnb(ji,jj) + ztnb(ji+1,jj) ) zs = 0.5 * ( zsnb(ji,jj) + zsnb(ji+1,jj) ) - 35.0 zh = 0.5 * ( zdep(ji,jj) + zdep(ji+1,jj) ) ! ... masked ratio alpha/beta zalbet = fsalbt( zt, zs, zh )*umask(ji,jj,1) ! ... local density gradient along i-bathymetric slope zgdrho = zalbet*( ztnb(ji+1,jj) - ztnb(ji,jj) ) & - ( zsnb(ji+1,jj) - zsnb(ji,jj) ) zgdrho = zgdrho * umask(ji,jj,1) ! ... sign of local i-gradient of density multiplied by the i-slope zsign = sign( 0.5, -zgdrho * ( zdep(ji+1,jj) - zdep(ji,jj) ) ) zki(ji,jj) = ( 0.5 - zsign ) * zahu(ji,jj) zsigna= sign(0.5, zunb(ji,jj)*( zdep(ji+1,jj) - zdep(ji,jj) )) zalphax(ji,jj)=(0.5+zsigna)*(0.5-zsign)*umask(ji,jj,1) END DO END DO DO jj = 1, jpjm1 DO ji = 1, fs_jpim1 ! vector opt. ! ... temperature, salinity anomalie and depth zt = 0.5 * ( ztnb(ji,jj+1) + ztnb(ji,jj) ) zs = 0.5 * ( zsnb(ji,jj+1) + zsnb(ji,jj) ) - 35.0 zh = 0.5 * ( zdep(ji,jj+1) + zdep(ji,jj) ) ! ... masked ratio alpha/beta zalbet = fsalbt( zt, zs, zh )*vmask(ji,jj,1) ! ... local density gradient along j-bathymetric slope zgdrho = zalbet*( ztnb(ji,jj+1) - ztnb(ji,jj) ) & - ( zsnb(ji,jj+1) - zsnb(ji,jj) ) zgdrho = zgdrho*vmask(ji,jj,1) ! ... sign of local j-gradient of density multiplied by the j-slope zsign = sign( 0.5, -zgdrho * ( zdep(ji,jj+1) - zdep(ji,jj) ) ) zkj(ji,jj) = ( 0.5 - zsign ) * zahv(ji,jj) zsigna= sign(0.5, zvnb(ji,jj)*(zdep(ji,jj+1) - zdep(ji,jj) ) ) zalphay(ji,jj)=(0.5+zsigna)*(0.5-zsign)*vmask(ji,jj,1) END DO END DO CASE ( 1 ) ! Linear formulation function of temperature only IF(lwp) WRITE(numout,cform_err) IF(lwp) WRITE(numout,*) ' use of linear eos rho(T,S) = rau0 * ( rbeta * S - ralpha * T )' IF(lwp) WRITE(numout,*) ' bbl not implented: easy to do it ' nstop = nstop + 1 CASE ( 2 ) ! Linear formulation function of temperature and salinity IF(lwp) WRITE(numout,cform_err) IF(lwp) WRITE(numout,*) ' use of linear eos rho(T,S) = rau0 * ( rbeta * S - ralpha * T )' IF(lwp) WRITE(numout,*) ' bbl not implented: easy to do it ' nstop = nstop + 1 CASE DEFAULT IF(lwp) WRITE(numout,cform_err) IF(lwp) WRITE(numout,*) ' bad flag value for neos = ', neos nstop = nstop + 1 END SELECT ! lateral boundary conditions on zalphax and zalphay (unchanged sign) CALL lbc_lnk( zalphax, 'U', 1. ) ; CALL lbc_lnk( zalphay, 'V', 1. ) ! 3. Velocities that are exchanged between ajacent bottom boxes. !--------------------------------------------------------------- ! ... is equal to zero but where bbl will work. u_trc_bbl(:,:,:) = 0.e0 v_trc_bbl(:,:,:) = 0.e0 # if defined key_vectopt_loop && ! defined key_autotasking jj = 1 DO ji = 1, jpij-jpi ! vector opt. (forced unrolling) # else DO jj = 1, jpjm1 DO ji = 1, jpim1 # endif iku = mbku(ji,jj) ikv = mbkv(ji,jj) IF( MAX(iku,ikv) > 1 ) THEN u_trc_bbl(ji,jj,iku) = zalphax(ji,jj) * un(ji,jj,iku) * umask(ji,jj,1) v_trc_bbl(ji,jj,ikv) = zalphay(ji,jj) * vn(ji,jj,ikv) * vmask(ji,jj,1) ENDIF # if ! defined key_vectopt_loop || defined key_autotasking END DO # endif END DO ! lateral boundary conditions on u_trc_bbl and v_trc_bbl (changed sign) CALL lbc_lnk( u_trc_bbl, 'U', -1. ) ; CALL lbc_lnk( v_trc_bbl, 'V', -1. ) DO jn = 1, jptra #if defined key_vectopt_loop && ! defined key_autotasking jj = 1 DO ji = 1, jpij ! vector opt. (forced unrolling) #else DO jj = 1, jpj DO ji = 1, jpi #endif ik = mbkt(ji,jj) ! index of the bottom ocean T-level ztrb(ji,jj) = trb(ji,jj,ik,jn) * tmask(ji,jj,1) ! masked now T at the ocean bottom #if ! defined key_vectopt_loop || defined key_autotasking END DO #endif END DO #if defined key_trcbbl_dif ! 4. Additional second order diffusive trends ! ------------------------------------------- ! ... first derivative (gradient) DO jj = 1, jpjm1 DO ji = 1, fs_jpim1 ! vertor opt. zkx(ji,jj) = zki(ji,jj)*( ztrb(ji+1,jj) - ztrb(ji,jj) ) zky(ji,jj) = zkj(ji,jj)*( ztrb(ji,jj+1) - ztrb(ji,jj) ) END DO END DO IF( cp_cfg == "orca" ) THEN SELECT CASE ( jp_cfg ) ! ! ======================= CASE ( 2 ) ! ORCA_R2 configuration ! ! ======================= ! Gibraltar enhancement of BBL zkx( mi0(139):mi1(140) , mj0(102):mj1(102) ) = 4.e0 * zkx( mi0(139):mi1(140) , mj0(102):mj1(102) ) zky( mi0(139):mi1(140) , mj0(102):mj1(102) ) = 4.e0 * zky( mi0(139):mi1(140) , mj0(102):mj1(102) ) ! Red Sea enhancement of BBL zkx( mi0(161):mi1(162) , mj0(88):mj1(88) ) = 10.e0 * zkx( mi0(161):mi1(162) , mj0(88):mj1(88) ) zky( mi0(161):mi1(162) , mj0(88):mj1(88) ) = 10.e0 * zky( mi0(161):mi1(162) , mj0(88):mj1(88) ) ! ! ======================= CASE ( 4 ) ! ORCA_R4 configuration ! ! ======================= ! Gibraltar enhancement of BBL zkx( mi0(70):mi1(71) , mj0(52):mj1(52) ) = 4.e0 * zkx( mi0(70):mi1(71) , mj0(52):mj1(52) ) zky( mi0(70):mi1(71) , mj0(52):mj1(52) ) = 4.e0 * zky( mi0(70):mi1(71) , mj0(52):mj1(52) ) END SELECT ENDIF ! ... second derivative (divergence) and add to the general tracer trend # if defined key_vectopt_loop && ! defined key_autotasking jj = 1 DO ji = jpi+2, jpij-jpi-1 ! vector opt. (forced unrolling) # else DO jj = 2, jpjm1 DO ji = 2, jpim1 # endif ik = mbkt(ji,jj) zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,ik) ) ztra = ( zkx(ji,jj) - zkx(ji-1,jj ) & & + zky(ji,jj) - zky(ji ,jj-1) ) * zbtr tra(ji,jj,ik,jn) = tra(ji,jj,ik,jn) + ztra #if ! defined key_vectopt_loop || defined key_autotasking END DO #endif END DO #endif ! 5. Along sigma advective trend ! ------------------------------- ! ... Second order centered tracer flux at u and v-points # if defined key_vectopt_loop && ! defined key_autotasking jj = 1 DO ji = 1, jpij-jpi ! vector opt. (forced unrolling) # else DO jj = 1, jpjm1 DO ji = 1, jpim1 # endif iku = mbku(ji,jj) ikv = mbkv(ji,jj) zfui = zalphax(ji,jj) *e2u(ji,jj) * fse3u(ji,jj,iku) * zunb(ji,jj) zfvj = zalphay(ji,jj) *e1v(ji,jj) * fse3v(ji,jj,ikv) * zvnb(ji,jj) ! upstream scheme zwx(ji,jj) = ( ( zfui + ABS( zfui ) ) * ztrb(ji ,jj ) & & +( zfui - ABS( zfui ) ) * ztrb(ji+1,jj ) ) * 0.5 zwy(ji,jj) = ( ( zfui + ABS( zfvj ) ) * ztrb(ji ,jj ) & & +( zfui - ABS( zfvj ) ) * ztrb(ji ,jj+1) ) * 0.5 #if ! defined key_vectopt_loop || defined key_autotasking END DO #endif END DO # if defined key_vectopt_loop && ! defined key_autotasking jj = 1 DO ji = jpi+2, jpij-jpi-1 ! vector opt. (forced unrolling) # else DO jj = 2, jpjm1 DO ji = 2, jpim1 # endif ik = mbkt(ji,jj) zbtr = 1. / ( e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,ik) ) ! horizontal advective trends ztra = - zbtr * ( zwx(ji,jj) - zwx(ji-1,jj ) & & + zwy(ji,jj) - zwy(ji ,jj-1) ) ! add it to the general tracer trends tra(ji,jj,ik,jn) = tra(ji,jj,ik,jn) + ztra #if ! defined key_vectopt_loop || defined key_autotasking END DO #endif END DO END DO IF(ln_ctl) THEN ! print mean trends (used for debugging) WRITE(charout, FMT="('bbl - adv')") CALL prt_ctl_trc_info(charout) CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm,clinfo2='trd') ENDIF ! 6. Vertical advection velocities ! -------------------------------- ! ... computes divergence perturbation (velocties to be removed from upper t boxes : DO jk= 1, jpkm1 DO jj=1, jpjm1 DO ji = 1, fs_jpim1 ! vertor opt. zwu(ji,jj) = -e2u(ji,jj) * u_trc_bbl(ji,jj,jk) zwv(ji,jj) = -e1v(ji,jj) * v_trc_bbl(ji,jj,jk) END DO END DO ! ... horizontal divergence DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. zbt = e1t(ji,jj) * e2t(ji,jj) zhdivn(ji,jj,jk) = ( zwu(ji,jj) - zwu(ji-1,jj ) & & + zwv(ji,jj) - zwv(ji ,jj-1) ) / zbt END DO END DO END DO ! ... horizontal bottom divergence # if defined key_vectopt_loop && ! defined key_autotasking jj = 1 DO ji = 1, jpij-jpi ! vector opt. (forced unrolling) # else DO jj = 1, jpjm1 DO ji = 1, jpim1 # endif iku = mbku(ji,jj) ikv = mbkv(ji,jj) zwu(ji,jj) = zalphax(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,iku) zwv(ji,jj) = zalphay(ji,jj) * e1v(ji,jj) * fse3v(ji,jj,ikv) #if ! defined key_vectopt_loop || defined key_autotasking END DO #endif END DO # if defined key_vectopt_loop && ! defined key_autotasking jj = 1 DO ji = jpi+2, jpij-jpi-1 ! vector opt. (forced unrolling) # else DO jj = 2, jpjm1 DO ji = 2, jpim1 # endif ik = mbkt(ji,jj) zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,ik) zhdivn(ji,jj,ik) = & & ( zwu(ji ,jj ) * ( zunb(ji ,jj ) - un(ji ,jj ,ik) *umask(ji ,jj ,1) ) & & - zwu(ji-1,jj ) * ( zunb(ji-1,jj ) - un(ji-1,jj ,ik) *umask(ji-1,jj ,1) ) & & + zwv(ji ,jj ) * ( zvnb(ji ,jj ) - vn(ji ,jj ,ik) *vmask(ji ,jj ,1) ) & & - zwv(ji ,jj-1) * ( zvnb(ji ,jj-1) - vn(ji ,jj-1,ik) *vmask(ji ,jj-1,1) ) & & ) / zbt # if ! defined key_vectopt_loop || defined key_autotasking END DO # endif END DO ! 7. compute additional vertical velocity to be used in t boxes ! ------------------------------------------------------------- ! ... Computation from the bottom ! Note that w_trc_bbl(:,:,jpk) has been set to 0 in trc_bbl_init DO jk = jpkm1, 1, -1 DO jj= 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. w_trc_bbl(ji,jj,jk) = w_trc_bbl(ji,jj,jk+1) - fse3t(ji,jj,jk)*zhdivn(ji,jj,jk) END DO END DO END DO ! Boundary condition on w_bbl (unchanged sign) CALL lbc_lnk( w_trc_bbl, 'W', 1. ) END SUBROUTINE trc_bbl_adv