MODULE trabbl !!============================================================================== !! *** MODULE trabbl *** !! Ocean physics : advective and/or diffusive bottom boundary layer scheme !!============================================================================== #if defined key_trabbl_dif || defined key_trabbl_adv || defined key_esopa !!---------------------------------------------------------------------- !! 'key_trabbl_dif' or diffusive bottom boundary layer !! 'key_trabbl_adv' advective bottom boundary layer !!---------------------------------------------------------------------- !! tra_bbl_dif : update the active tracer trends due to the bottom !! boundary layer (diffusive only) !! tra_bbl_adv : update the active tracer trends due to the bottom !! boundary layer (advective and/or diffusive) !! tra_bbl_init : initialization, namlist read, parameters control !!---------------------------------------------------------------------- !! * Modules used USE oce ! ocean dynamics and active tracers USE dom_oce ! ocean space and time domain USE trdtra_oce ! ocean active tracer trends USE in_out_manager ! I/O manager IMPLICIT NONE PRIVATE !! * Routine accessibility PUBLIC tra_bbl_dif ! routine called by step.F90 PUBLIC tra_bbl_adv ! routine called by step.F90 !! * Shared module variables LOGICAL, PUBLIC, PARAMETER :: & lk_trabbl_dif = .TRUE. ! diffusive bottom boundary layer flag # if defined key_trabbl_adv LOGICAL, PUBLIC, PARAMETER :: & lk_trabbl_adv = .TRUE. ! bottom boundary layer flag REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: & w_bbl, & ! vertical increment of velocity due to advective BBL ! ! only affect tracer vertical advection u_bbl, v_bbl ! velocity involved in exhanges in the advective BBL # else LOGICAL, PUBLIC, PARAMETER :: & lk_trabbl_adv = .FALSE. ! advective bottom boundary layer flag # endif !! * Module variables INTEGER, DIMENSION(jpi,jpj) :: & mbkt, mbku, mbkv ! ??? REAL(wp) :: & !!! * bbl namelist * atrbbl = 1.e+3 ! lateral coeff. for bottom boundary layer scheme (m2/s) !! * Substitutions # include "domzgr_substitute.h90" # include "vectopt_loop_substitute.h90" !!---------------------------------------------------------------------- !! OPA 9.0 , LODYC-IPSL (2003) !!---------------------------------------------------------------------- CONTAINS SUBROUTINE tra_bbl_dif( kt ) !!---------------------------------------------------------------------- !! *** ROUTINE tra_bbl_dif *** !! !! ** Purpose : Compute the before tracer (t & s) 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 !! a purely diffusive bottom boundary layer. !! !! ** Method : When the product grad( rho) * grad(h) < 0 (where grad !! is an along bottom slope gradient) an additional lateral diffu- !! sive trend along the bottom slope is added to the general tracer !! trend, otherwise nothing is done. !! 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 temperature 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 (ta,sa) of the !! botton ocean tracer point: !! ta = ta + difft !! !! ** Action : - update (ta,sa) at the bottom level with the bottom !! boundary layer trend !! - save the trends in bbltrd ('key_diatrends') !! !! References : !! Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591. !! !! History : !! 8.0 ! 96-06 (L. Mortier) Original code !! 8.0 ! 97-11 (G. Madec) Optimization !! 8.5 ! 02-08 (G. Madec) free form + modules !!---------------------------------------------------------------------- !! * Arguments INTEGER, INTENT( in ) :: kt ! ocean time-step !! * Local declarations INTEGER :: ji, jj ! dummy loop indices INTEGER :: ik # if defined key_partial_steps INTEGER :: iku1, iku2, ikv1,ikv2 ! temporary intergers REAL(wp) :: ze3u, ze3v ! temporary scalars # else INTEGER :: iku, ikv # endif REAL(wp) :: & zsign, zt, zs, zh, zalbet, & ! temporary scalars zgdrho, zbtr, zta, zsa REAL(wp), DIMENSION(jpi,jpj) :: & zki, zkj, zkw, zkx, zky, zkz, & ! temporary workspace arrays ztnb, zsnb, zdep, & ztbb, zsbb, zahu, zahv REAL(wp) :: & fsalbt, pft, pfs, pfh ! statement function !!---------------------------------------------------------------------- ! 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 == nit000 ) CALL tra_bbl_init ! 0. 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 and S at ocean bottom zsnb(ji,jj) = sn(ji,jj,ik) * tmask(ji,jj,1) ztbb(ji,jj) = tb(ji,jj,ik) * tmask(ji,jj,1) ! masked before T and S at ocean bottom zsbb(ji,jj) = sb(ji,jj,ik) * tmask(ji,jj,1) 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_partial_steps ! partial steps correction # 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 iku1 = MAX( mbathy(ji+1,jj )-1, 1 ) iku2 = MAX( mbathy(ji ,jj )-1, 1 ) ikv1 = MAX( mbathy(ji ,jj+1)-1, 1 ) ikv2 = MAX( mbathy(ji ,jj )-1, 1 ) ze3u = MIN( fse3u(ji,jj,iku1), fse3u(ji,jj,iku2) ) ze3v = MIN( fse3v(ji,jj,ikv1), fse3v(ji,jj,ikv2) ) zahu(ji,jj) = atrbbl * e2u(ji,jj) * ze3u / e1u(ji,jj) * umask(ji,jj,1) zahv(ji,jj) = atrbbl * e1v(ji,jj) * ze3v / e2v(ji,jj) * vmask(ji,jj,1) # if ! defined key_vectopt_loop || defined key_autotasking END DO # endif END DO # else # 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 ! 1. 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 # 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 ! 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) ) ! 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) # 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) # else DO jj = 1, jpjm1 DO ji = 1, jpim1 # endif ! 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) ) ! 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) # if ! defined key_vectopt_loop || defined key_autotasking END DO # endif END DO ! 2. Additional second order diffusive trends ! ------------------------------------------- ! first derivative (gradient) # 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 zkx(ji,jj) = zki(ji,jj) * ( ztbb(ji+1,jj) - ztbb(ji,jj) ) zkz(ji,jj) = zki(ji,jj) * ( zsbb(ji+1,jj) - zsbb(ji,jj) ) zky(ji,jj) = zkj(ji,jj) * ( ztbb(ji,jj+1) - ztbb(ji,jj) ) zkw(ji,jj) = zkj(ji,jj) * ( zsbb(ji,jj+1) - zsbb(ji,jj) ) # if ! defined key_vectopt_loop || defined key_autotasking END DO # endif 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 = max( mbathy(ji,jj)-1, 1 ) zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,ik) ) zta = ( zkx(ji,jj) - zkx(ji-1,jj ) & + zky(ji,jj) - zky(ji ,jj-1) ) * zbtr zsa = ( zkz(ji,jj) - zkz(ji-1,jj ) & + zkw(ji,jj) - zkw(ji ,jj-1) ) * zbtr ta(ji,jj,ik) = ta(ji,jj,ik) + zta sa(ji,jj,ik) = sa(ji,jj,ik) + zsa # if defined key_trdtra || defined key_trdmld bbltrd(ji,jj,1) = + zta bbltrd(ji,jj,2) = + zsa # endif # if ! defined key_vectopt_loop || defined key_autotasking END DO # endif END DO IF( l_ctl .AND. lwp ) THEN zta = SUM( ta(2:jpim1,2:jpjm1,1:jpkm1) * tmask(2:jpim1,2:jpjm1,1:jpkm1) ) zsa = SUM( sa(2:jpim1,2:jpjm1,1:jpkm1) * tmask(2:jpim1,2:jpjm1,1:jpkm1) ) WRITE(numout,*) ' bbl - Ta: ', zta-t_ctl, ' Sa: ', zsa-s_ctl t_ctl = zta ; s_ctl = zsa ENDIF END SUBROUTINE tra_bbl_dif # if defined key_trabbl_adv !!---------------------------------------------------------------------- !! 'key_trabbl_adv' advective bottom boundary layer !!---------------------------------------------------------------------- # include "trabbl_adv.h90" # else !!---------------------------------------------------------------------- !! Default option : NO advective bottom boundary layer !!---------------------------------------------------------------------- SUBROUTINE tra_bbl_adv (kt ) ! Empty routine INTEGER, INTENT(in) :: kt WRITE(*,*) kt END SUBROUTINE tra_bbl_adv # endif SUBROUTINE tra_bbl_init !!---------------------------------------------------------------------- !! *** ROUTINE tra_bbl_init *** !! !! ** Purpose : Initialization for the bottom boundary layer scheme. !! !! ** Method : Read the nambbl namelist and check the parameters !! called by tra_bbl at the first timestep (nit000) !! !! History : !! 8.5 ! 02-08 (G. Madec) Original code !!---------------------------------------------------------------------- !! * Local declarations INTEGER :: ji, jj ! dummy loop indices NAMELIST/nambbl/ atrbbl !!---------------------------------------------------------------------- ! Read Namelist nambbl : bottom boundary layer scheme ! -------------------- REWIND ( numnam ) READ ( numnam, nambbl ) ! Parameter control and print ! --------------------------- IF(lwp) THEN WRITE(numout,*) WRITE(numout,*) 'tra_bbl_init : * Diffusive Bottom Boundary Layer' WRITE(numout,*) '~~~~~~~~~~~~' IF( lk_trabbl_adv ) THEN WRITE(numout,*) ' * Advective Bottom Boundary Layer' ENDIF WRITE(numout,*) ' Namelist nambbl : set bbl parameters' WRITE(numout,*) WRITE(numout,*) ' bottom boundary layer coef. atrbbl = ', atrbbl WRITE(numout,*) ENDIF DO jj = 1, jpj DO ji = 1, jpi mbkt(ji,jj) = MAX( mbathy(ji,jj) - 1, 1 ) ! vertical index of the bottom ocean T-level END DO END DO DO jj = 1, jpjm1 DO ji = 1, jpim1 mbku(ji,jj) = MAX( MIN( mbathy(ji+1,jj ), mbathy(ji,jj) ) - 1, 1 ) mbkv(ji,jj) = MAX( MIN( mbathy(ji ,jj+1), mbathy(ji,jj) ) - 1, 1 ) END DO END DO !!bug ??? !!bug Caution : define the vakue of mbku & mbkv everywhere!!! but lbc mpp lnk : pb when closed (0) # if defined key_trabbl_adv w_bbl(:,:,:) = 0.e0 ! initialisation of w_bbl to zero # endif END SUBROUTINE tra_bbl_init #else !!---------------------------------------------------------------------- !! Dummy module : No bottom boundary layer scheme !!---------------------------------------------------------------------- LOGICAL, PUBLIC, PARAMETER :: lk_trabbl_dif = .FALSE. LOGICAL, PUBLIC, PARAMETER :: lk_trabbl_adv = .FALSE. CONTAINS SUBROUTINE tra_bbl_dif (kt ) ! Empty routine INTEGER, INTENT(in) :: kt WRITE(*,*) kt END SUBROUTINE tra_bbl_dif SUBROUTINE tra_bbl_adv (kt ) ! Empty routine INTEGER, INTENT(in) :: kt WRITE(*,*) kt END SUBROUTINE tra_bbl_adv #endif !!====================================================================== END MODULE trabbl