MODULE zdftke_crs !!====================================================================== !! *** MODULE zdftke *** !! Ocean physics: vertical mixing coefficient computed from the tke !! turbulent closure parameterization !!===================================================================== !! History : OPA ! 1991-03 (b. blanke) Original code !! 7.0 ! 1991-11 (G. Madec) bug fix !! 7.1 ! 1992-10 (G. Madec) new mixing length and eav !! 7.2 ! 1993-03 (M. Guyon) symetrical conditions !! 7.3 ! 1994-08 (G. Madec, M. Imbard) nn_pdl flag !! 7.5 ! 1996-01 (G. Madec) s-coordinates !! 8.0 ! 1997-07 (G. Madec) lbc !! 8.1 ! 1999-01 (E. Stretta) new option for the mixing length !! NEMO 1.0 ! 2002-06 (G. Madec) add tke_init routine !! - ! 2004-10 (C. Ethe ) 1D configuration !! 2.0 ! 2006-07 (S. Masson) distributed restart using iom !! 3.0 ! 2008-05 (C. Ethe, G.Madec) : update TKE physics: !! ! - tke penetration (wind steering) !! ! - suface condition for tke & mixing length !! ! - Langmuir cells !! - ! 2008-05 (J.-M. Molines, G. Madec) 2D form of avtb !! - ! 2008-06 (G. Madec) style + DOCTOR name for namelist parameters !! - ! 2008-12 (G. Reffray) stable discretization of the production term !! 3.2 ! 2009-06 (G. Madec, S. Masson) TKE restart compatible with key_cpl !! ! + cleaning of the parameters + bugs correction !! 3.3 ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase !!---------------------------------------------------------------------- #if defined key_zdftke || defined key_esopa !!---------------------------------------------------------------------- !! 'key_zdftke' TKE vertical physics !!---------------------------------------------------------------------- !! zdf_tke : update momentum and tracer Kz from a tke scheme !! tke_tke : tke time stepping: update tke at now time step (en) !! tke_avn : compute mixing length scale and deduce avm and avt !! zdf_tke_init : initialization, namelist read, and parameters control !! tke_rst : read/write tke restart in ocean restart file !!---------------------------------------------------------------------- USE crs USE zdf_oce , ONLY : avtb, avmb, avtb_2d USE zdftke USE crslbclnk !USE oce ! ocean: dynamics and active tracers variables USE phycst ! physical constants !USE dom_oce ! domain: ocean !USE domvvl ! domain: variable volume layer !USE sbc_oce ! surface boundary condition: ocean !USE zdf_oce ! vertical physics: ocean variables !USE zdfmxl ! vertical physics: mixed layer !USE lbclnk ! ocean lateral boundary conditions (or mpp link) !USE prtctl ! Print control !USE in_out_manager ! I/O manager !USE iom ! I/O manager library !USE lib_mpp ! MPP library USE wrk_nemo ! work arrays USE timing ! Timing !USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) IMPLICIT NONE PRIVATE PUBLIC tke_avn_crs PUBLIC tke_avn_ini_crs !LOGICAL , PUBLIC, PARAMETER :: lk_zdftke = .TRUE. !: TKE vertical mixing flag ! !!** Namelist namzdf_tke ** !LOGICAL :: ln_mxl0 ! mixing length scale surface value as function of wind stress or not !INTEGER :: nn_mxl ! type of mixing length (=0/1/2/3) !REAL(wp) :: rn_mxl0 ! surface min value of mixing length (kappa*z_o=0.4*0.1 m) [m] !INTEGER :: nn_pdl ! Prandtl number or not (ratio avt/avm) (=0/1) !REAL(wp) :: rn_ediff ! coefficient for avt: avt=rn_ediff*mxl*sqrt(e) !REAL(wp) :: rn_ediss ! coefficient of the Kolmogoroff dissipation !REAL(wp) :: rn_ebb ! coefficient of the surface input of tke !REAL(wp) :: rn_emin ! minimum value of tke [m2/s2] !REAL(wp) :: rn_emin0 ! surface minimum value of tke [m2/s2] !REAL(wp) :: rn_bshear ! background shear (>0) currently a numerical threshold (do not change it) !INTEGER :: nn_etau ! type of depth penetration of surface tke (=0/1/2/3) !INTEGER :: nn_htau ! type of tke profile of penetration (=0/1) !REAL(wp) :: rn_efr ! fraction of TKE surface value which penetrates in the ocean !LOGICAL :: ln_lc ! Langmuir cells (LC) as a source term of TKE or not !REAL(wp) :: rn_lc ! coef to compute vertical velocity of Langmuir cells !REAL(wp) :: ri_cri ! critic Richardson number (deduced from rn_ediff and rn_ediss values) !REAL(wp) :: rmxl_min ! minimum mixing length value (deduced from rn_ediff and rn_emin values) [m] !REAL(wp) :: rhftau_add = 1.e-3_wp ! add offset applied to HF part of taum (nn_etau=3) !REAL(wp) :: rhftau_scl = 1.0_wp ! scale factor applied to HF part of taum (nn_etau=3) !REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: en !: now turbulent kinetic energy [m2/s2] !REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: htau ! depth of tke penetration (nn_htau) !REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: dissl ! now mixing lenght of dissipation !REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: avt_k , avm_k ! not enhanced Kz !REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: avmu_k, avmv_k ! not enhanced Kz #if defined key_c1d ! !!** 1D cfg only ** ('key_c1d') !REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e_dis, e_mix !: dissipation and mixing turbulent lengh scales !REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e_pdl, e_ric !: prandl and local Richardson numbers #endif !! * Substitutions # include "domzgr_substitute.h90" # include "vectopt_loop_substitute.h90" !!---------------------------------------------------------------------- !! NEMO/OPA 4.0 , NEMO Consortium (2011) !! $Id: zdftke.F90 4990 2014-12-15 16:42:49Z timgraham $ !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) !!---------------------------------------------------------------------- CONTAINS ! INTEGER FUNCTION zdf_tke_alloc() ! !!---------------------------------------------------------------------- ! !! *** FUNCTION zdf_tke_alloc *** ! !!---------------------------------------------------------------------- ! ALLOCATE( & !#if defined key_c1d ! & e_dis(jpi,jpj,jpk) , e_mix(jpi,jpj,jpk) , & ! & e_pdl(jpi,jpj,jpk) , e_ric(jpi,jpj,jpk) , & !#endif ! & en (jpi,jpj,jpk) , htau (jpi,jpj) , dissl(jpi,jpj,jpk) , & ! & avt_k (jpi,jpj,jpk) , avm_k (jpi,jpj,jpk), & ! & avmu_k(jpi,jpj,jpk) , avmv_k(jpi,jpj,jpk), STAT= zdf_tke_alloc ) ! ! ! IF( lk_mpp ) CALL mpp_sum ( zdf_tke_alloc ) ! IF( zdf_tke_alloc /= 0 ) CALL ctl_warn('zdf_tke_alloc: failed to allocate arrays') ! ! ! END FUNCTION zdf_tke_alloc SUBROUTINE tke_avn_crs !!---------------------------------------------------------------------- !! *** ROUTINE tke_avn *** !! !! ** Purpose : Compute the vertical eddy viscosity and diffusivity !! !! ** Method : At this stage, en, the now TKE, is known (computed in !! the tke_tke routine). First, the now mixing lenth is !! computed from en and the strafification (N^2), then the mixings !! coefficients are computed. !! - Mixing length : a first evaluation of the mixing lengh !! scales is: !! mxl = sqrt(2*en) / N !! where N is the brunt-vaisala frequency, with a minimum value set !! to rmxl_min (rn_mxl0) in the interior (surface) ocean. !! The mixing and dissipative length scale are bound as follow : !! nn_mxl=0 : mxl bounded by the distance to surface and bottom. !! zmxld = zmxlm = mxl !! nn_mxl=1 : mxl bounded by the e3w and zmxld = zmxlm = mxl !! nn_mxl=2 : mxl bounded such that the vertical derivative of mxl is !! less than 1 (|d/dz(mxl)|<1) and zmxld = zmxlm = mxl !! nn_mxl=3 : mxl is bounded from the surface to the bottom usings !! |d/dz(xml)|<1 to obtain lup, and from the bottom to !! the surface to obtain ldown. the resulting length !! scales are: !! zmxld = sqrt( lup * ldown ) !! zmxlm = min ( lup , ldown ) !! - Vertical eddy viscosity and diffusivity: !! avm = max( avtb, rn_ediff * zmxlm * en^1/2 ) !! avt = max( avmb, pdlr * avm ) !! with pdlr=1 if nn_pdl=0, pdlr=1/pdl=F(Ri) otherwise. !! !! ** Action : - avt : now vertical eddy diffusivity (w-point) !! - avmu, avmv : now vertical eddy viscosity at uw- and vw-points !!---------------------------------------------------------------------- INTEGER :: ji, jj, jk ! dummy loop indices REAL(wp) :: zrn2, zraug, zcoef, zav ! local scalars REAL(wp) :: zdku, zdkv, zpdlr, zri, zsqen ! - - REAL(wp) :: zemlm, zemlp ! - - !REAL(wp), POINTER, DIMENSION(:,:,:) :: zmpdl, zmxlm, zmxld REAL(wp), POINTER, DIMENSION(:,:,:) :: zmxlm_crs,zmxld_crs REAL(wp), POINTER, DIMENSION(:,:,:) :: avm_crs, avmu_crs,avmv_crs !!-------------------------------------------------------------------- ! IF( nn_timing == 1 ) CALL timing_start('tke_avn') !CALL wrk_alloc( jpi,jpj,jpk, zmpdl, zmxlm, zmxld ) CALL wrk_alloc( jpi_crs,jpj_crs,jpk, zmxlm_crs, zmxld_crs, avm_crs, avmu_crs, avmv_crs ) ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< ! ! Mixing length ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< ! ! !* Buoyancy length scale: l=sqrt(2*e/n**2) ! IF( ln_mxl0 ) THEN ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rau0*g) DO jj = 2, jpj_crs-1 DO ji = 2, jpi_crs-1 !IF (mikt(ji,jj) .GT. 1) THEN ! zmxlm(ji,jj,mikt(ji,jj)) = rmxl_min !ELSE zraug = vkarmn * 2.e5_wp / ( rau0 * grav ) zmxlm_crs(ji,jj,1) = MAX( rn_mxl0, zraug * taum_crs(ji,jj) ) !END IF END DO END DO !ELSE ! DO jj = 2, jpjm1 ! DO ji = fs_2, fs_jpim1 ! surface set to the minimum value ! zmxlm(ji,jj,mikt(ji,jj)) = MAX( tmask(ji,jj,1) * rn_mxl0, rmxl_min) ! END DO ! END DO ENDIF zmxlm_crs(:,:,jpk) = rmxl_min ! last level set to the interior minium value ! DO jj = 2, jpj_crs-1 DO ji = 2, jpi_crs-1 ! vector opt. DO jk = 2, jpk-1 ! interior value : l=sqrt(2*e/n^2) zrn2 = MAX( rn2_crs(ji,jj,jk), rsmall ) zmxlm_crs(ji,jj,jk) = MAX( rmxl_min, SQRT( 2._wp * en_crs(ji,jj,jk) / zrn2 ) ) END DO zmxld_crs(ji,jj,1) = zmxlm_crs(ji,jj,1) ! surface set to the minimum value END DO END DO ! ! !* Physical limits for the mixing length ! zmxld_crs(:,:, 1 ) = zmxlm_crs(:,:,1) ! surface set to the zmxlm value zmxld_crs(:,:,jpk) = rmxl_min ! last level set to the minimum value ! SELECT CASE ( nn_mxl ) ! !CASE ( 0 ) ! bounded by the distance to surface and bottom ! DO jj = 2, jpjm1 ! DO ji = fs_2, fs_jpim1 ! vector opt. ! DO jk = mikt(ji,jj)+1, jpkm1 ! zemxl = MIN( fsdepw(ji,jj,jk) - fsdepw(ji,jj,mikt(ji,jj)), zmxlm(ji,jj,jk), & ! & fsdepw(ji,jj,mbkt(ji,jj)+1) - fsdepw(ji,jj,jk) ) ! zmxlm(ji,jj,jk) = zemxl ! zmxld(ji,jj,jk) = zemxl ! END DO ! END DO ! END DO ! ! !CASE ( 1 ) ! bounded by the vertical scale factor ! DO jj = 2, jpjm1 ! DO ji = fs_2, fs_jpim1 ! vector opt. ! DO jk = mikt(ji,jj)+1, jpkm1 ! zemxl = MIN( fse3w(ji,jj,jk), zmxlm(ji,jj,jk) ) ! zmxlm(ji,jj,jk) = zemxl ! zmxld(ji,jj,jk) = zemxl ! END DO ! END DO ! END DO ! ! !CASE ( 2 ) ! |dk[xml]| bounded by e3t : ! DO jj = 2, jpjm1 ! DO ji = fs_2, fs_jpim1 ! vector opt. ! DO jk = mikt(ji,jj)+1, jpkm1 ! from the surface to the bottom : ! zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk-1) + fse3t(ji,jj,jk-1), zmxlm(ji,jj,jk) ) ! END DO ! DO jk = jpkm1, mikt(ji,jj)+1, -1 ! from the bottom to the surface : ! zemxl = MIN( zmxlm(ji,jj,jk+1) + fse3t(ji,jj,jk+1), zmxlm(ji,jj,jk) ) ! zmxlm(ji,jj,jk) = zemxl ! zmxld(ji,jj,jk) = zemxl ! END DO ! END DO ! END DO ! ! CASE ( 3 ) ! lup and ldown, |dk[xml]| bounded by e3t : DO jj = 2, jpj_crs-1 DO ji = 2, jpi_crs-1 ! vector opt. DO jk = 2, jpkm1 ! from the surface to the bottom : lup zmxld_crs(ji,jj,jk) = MIN( zmxld_crs(ji,jj,jk-1) + e3t_crs(ji,jj,jk-1), zmxlm_crs(ji,jj,jk) ) END DO DO jk = jpkm1, 2 , -1 ! from the bottom to the surface : ldown zmxlm_crs(ji,jj,jk) = MIN( zmxlm_crs(ji,jj,jk+1) + e3t_crs(ji,jj,jk+1), zmxlm_crs(ji,jj,jk) ) END DO END DO END DO DO jk = 2, jpkm1 DO jj = 2, jpj_crs-1 DO ji = 2, jpi_crs-1 ! vector opt. zemlm = MIN ( zmxld_crs(ji,jj,jk), zmxlm_crs(ji,jj,jk) ) zemlp = SQRT( zmxld_crs(ji,jj,jk) * zmxlm_crs(ji,jj,jk) ) zmxlm_crs(ji,jj,jk) = zemlm zmxld_crs(ji,jj,jk) = zemlp END DO END DO END DO ! END SELECT ! ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< ! ! Vertical eddy viscosity and diffusivity (avmu, avmv, avt) ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< DO jk = 1, jpkm1 !* vertical eddy viscosity & diffivity at w-points DO jj = 2, jpj_crs-1 DO ji = 2, jpi_crs-1 ! vector opt. zsqen = SQRT( en_crs(ji,jj,jk) ) zav = rn_ediff * zmxlm_crs(ji,jj,jk) * zsqen avm_crs(ji,jj,jk) = MAX( zav, avmb(jk) ) * tmask_crs(ji,jj,jk) avt_crs(ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * tmask_crs(ji,jj,jk) !dissl(ji,jj,jk) = zsqen / zmxld(ji,jj,jk) END DO END DO END DO CALL crs_lbc_lnk( avm_crs, 'W', 1. ) ! Lateral boundary conditions (sign unchanged) ! DO jj = 2, jpj_crs-1 DO ji = 2, jpi_crs-1 ! vector opt. DO jk = 2, jpkm1 !* vertical eddy viscosity at u- and v-points avmu_crs(ji,jj,jk) = 0.5 * ( avm_crs(ji,jj,jk) + avm_crs(ji+1,jj ,jk) ) * umask_crs(ji,jj,jk) END DO DO jk = 2, jpkm1 avmv_crs(ji,jj,jk) = 0.5 * ( avm_crs(ji,jj,jk) + avm_crs(ji ,jj+1,jk) ) * vmask_crs(ji,jj,jk) END DO END DO END DO CALL crs_lbc_lnk( avmu_crs, 'U', 1. ) ; CALL crs_lbc_lnk( avmv_crs, 'V', 1. ) ! Lateral boundary conditions ! IF( nn_pdl == 1 ) THEN !* Prandtl number case: update avt DO jj = 2, jpj_crs-1 DO ji = 2, jpi_crs-1 ! vector opt. DO jk = 2, jpkm1 zcoef = avm_crs(ji,jj,jk) * 2._wp * e3w_crs(ji,jj,jk) * e3w_crs(ji,jj,jk) ! ! shear zdku = avmu_crs(ji-1,jj,jk) * ( un_crs(ji-1,jj,jk-1) - un_crs(ji-1,jj,jk) ) * ( ub_crs(ji-1,jj,jk-1) - ub_crs(ji-1,jj,jk) ) & & + avmu_crs(ji ,jj,jk) * ( un_crs(ji ,jj,jk-1) - un_crs(ji ,jj,jk) ) * ( ub_crs(ji ,jj,jk-1) - ub_crs(ji ,jj,jk) ) zdkv = avmv_crs(ji,jj-1,jk) * ( vn_crs(ji,jj-1,jk-1) - vn_crs(ji,jj-1,jk) ) * ( vb_crs(ji,jj-1,jk-1) - vb_crs(ji,jj-1,jk) ) & & + avmv_crs(ji,jj ,jk) * ( vn_crs(ji,jj ,jk-1) - vn_crs(ji,jj ,jk) ) * ( vb_crs(ji,jj ,jk-1) - vb_crs(ji,jj ,jk) ) ! ! local Richardson number zri = MAX( rb2_crs(ji,jj,jk), 0._wp ) * zcoef / (zdku + zdkv + rn_bshear ) zpdlr = MAX( 0.1_wp, 0.2 / MAX( 0.2 , zri ) ) avt_crs(ji,jj,jk) = MAX( zpdlr * avt_crs(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * tmask_crs(ji,jj,jk) END DO END DO END DO ENDIF CALL crs_lbc_lnk( avt_crs, 'W', 1. ) ! Lateral boundary conditions on avt (sign unchanged) !IF(ln_ctl) THEN ! CALL prt_ctl( tab3d_1=en , clinfo1=' tke - e: ', tab3d_2=avt, clinfo2=' t: ', ovlap=1, kdim=jpk) ! CALL prt_ctl( tab3d_1=avmu, clinfo1=' tke - u: ', mask1=umask, & ! & tab3d_2=avmv, clinfo2= ' v: ', mask2=vmask, ovlap=1, kdim=jpk ) !ENDIF ! CALL wrk_dealloc( jpi_crs,jpj_crs,jpk, zmxlm_crs, zmxld_crs, avm_crs, avmu_crs, avmv_crs ) ! IF( nn_timing == 1 ) CALL timing_stop('tke_avn') ! END SUBROUTINE tke_avn_crs SUBROUTINE tke_avn_ini_crs !!---------------------------------------------------------------------- !! !! !! !! !!---------------------------------------------------------------------- !IF( nn_avb == 0 ) THEN ! Define avmb, avtb from namelist parameter ! avtb(:) = rn_avt0 !ELSE ! Background profile of avt (fit a theoretical/observational profile (Krauss 1990) ! avmb(:) = rn_avm0 ! avtb(:) = rn_avt0 + ( 3.e-4_wp - 2._wp * rn_avt0 ) * 1.e-4_wp * gdepw_1d(:) ! m2/s ! IF(ln_sco .AND. lwp) CALL ctl_warn( 'avtb profile not valid in sco' ) !ENDIF avtb_2d_crs(:,:) = 1.e0 !IF( nn_havtb == 1 ) THEN ! decrease avtb in the equatorial band ! ! -15S -5S : linear decrease from avt0 to avt0/10. ! ! -5S +5N : cst value avt0/10. ! ! 5N 15N : linear increase from avt0/10, to avt0 ! WHERE(-15. <= gphit .AND. gphit < -5 ) avtb_2d = (1. - 0.09 * (gphit + 15.)) ! WHERE( -5. <= gphit .AND. gphit < 5 ) avtb_2d = 0.1 ! WHERE( 5. <= gphit .AND. gphit < 15 ) avtb_2d = (0.1 + 0.09 * (gphit - 5.)) !ENDIF END SUBROUTINE tke_avn_ini_crs #else !!---------------------------------------------------------------------- !! Dummy module : NO TKE scheme !!---------------------------------------------------------------------- LOGICAL, PUBLIC, PARAMETER :: lk_zdftke = .FALSE. !: TKE flag #endif !!====================================================================== END MODULE zdftke_crs