Index: LMDZ6/trunk/libf/phylmd/atke_exchange_coeff_mod.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/atke_exchange_coeff_mod.F90 (revision 4630) +++ LMDZ6/trunk/libf/phylmd/atke_exchange_coeff_mod.F90 (revision 4631) @@ -5,5 +5,5 @@ contains -subroutine atke_compute_km_kh(ngrid,nlay, & +subroutine atke_compute_km_kh(ngrid,nlay,dtime, & wind_u,wind_v,temp,play,pinterf, & tke,Km_out,Kh_out) @@ -28,5 +28,5 @@ USE atke_turbulence_ini_mod, ONLY : iflag_atke, kappa, l0, ric, cinf, rpi, rcpd USE atke_turbulence_ini_mod, ONLY : cepsilon, pr_slope, pr_asym, pr_neut, rg, rd -USE atke_turbulence_ini_mod, ONLY : viscom, viscoh +USE atke_turbulence_ini_mod, ONLY : viscom, viscoh, clmix, iflag_atke_lmix, lmin implicit none @@ -37,6 +37,7 @@ INTEGER, INTENT(IN) :: ngrid ! number of horizontal index (flat grid) -INTEGER, INTENT(IN) :: nlay ! number of vertical index - +INTEGER, INTENT(IN) :: nlay ! number of vertical index + +REAL, INTENT(IN) :: dtime ! physics time step (s) REAL, DIMENSION(ngrid,nlay), INTENT(IN) :: wind_u ! zonal velocity (m/s) REAL, DIMENSION(ngrid,nlay), INTENT(IN) :: wind_v ! meridional velocity (m/s) @@ -60,8 +61,11 @@ REAL, DIMENSION(ngrid,nlay) :: dz_interf ! distance between two consecutive interfaces REAL, DIMENSION(ngrid,nlay) :: dz_lay ! distance between two layer middles (NB: first and last are half layers) +REAL, DIMENSION(ngrid,nlay+1) :: N2 ! square of Brunt Vaisala pulsation (at interface) +REAL, DIMENSION(ngrid,nlay+1) :: shear2 ! square of wind shear (at interface) REAL, DIMENSION(ngrid,nlay+1) :: Ri ! Richardson's number (at interface) REAL, DIMENSION(ngrid,nlay+1) :: Prandtl ! Turbulent Prandtl's number (at interface) REAL, DIMENSION(ngrid,nlay+1) :: Sm ! Stability function for momentum (at interface) REAL, DIMENSION(ngrid,nlay+1) :: Sh ! Stability function for heat (at interface) +LOGICAL, DIMENSION(ngrid,nlay+1) :: switch_num ! switch of numerical integration in very stable cases INTEGER :: igrid,ilay ! horizontal,vertical index (flat grid) @@ -69,5 +73,9 @@ REAL :: preff ! reference pressure for potential temperature calculations REAL :: thetam ! mean potential temperature at interface - +REAL :: delta ! discriminant of the second order polynomial +REAL :: qq ! tke=qq**2/2 +REAL :: shear ! wind shear +REAL :: lstrat ! mixing length depending on local stratification +REAL :: taustrat ! caracteristic timescale for turbulence in very stable conditions ! Initializations: @@ -109,15 +117,4 @@ -! Computing the mixing length: -! so far, we have neglected the effect of local stratification -!============================================================== - - -DO ilay=2,nlay+1 - DO igrid=1,ngrid - l_exchange(igrid,ilay) = kappa*l0*z_interf(igrid,ilay) / (kappa*z_interf(igrid,ilay) + l0) - ENDDO -ENDDO - ! Computes the gradient Richardson's number and stability functions: !=================================================================== @@ -136,7 +133,8 @@ dz_lay(igrid,ilay)=z_lay(igrid,ilay)-z_lay(igrid,ilay-1) thetam=0.5*(theta(igrid,ilay) + theta(igrid,ilay-1)) - Ri(igrid,ilay) = rg * (theta(igrid,ilay) - theta(igrid,ilay-1))/thetam / dz_lay(igrid,ilay) / & - MAX(((wind_u(igrid,ilay) - wind_u(igrid,ilay-1)) / dz_lay(igrid,ilay))**2 + & - ((wind_v(igrid,ilay) - wind_v(igrid,ilay-1)) / dz_lay(igrid,ilay))**2,1E-10) + N2(igrid,ilay) = rg * (theta(igrid,ilay) - theta(igrid,ilay-1))/thetam / dz_lay(igrid,ilay) + shear2(igrid,ilay)= (((wind_u(igrid,ilay) - wind_u(igrid,ilay-1)) / dz_lay(igrid,ilay))**2 + & + ((wind_v(igrid,ilay) - wind_v(igrid,ilay-1)) / dz_lay(igrid,ilay))**2 ) + Ri(igrid,ilay) = N2(igrid,ilay) / MAX(shear2(igrid,ilay),1E-10) IF (Ri(igrid,ilay) < 0.) THEN ! unstable cases @@ -157,22 +155,104 @@ ENDDO + +! Computing the mixing length: +!============================================================== + +switch_num(:,:)=.false. + +IF (iflag_atke_lmix .EQ. 1 ) THEN + + DO ilay=2,nlay + DO igrid=1,ngrid + l_exchange(igrid,ilay) = kappa*l0*z_interf(igrid,ilay) / (kappa*z_interf(igrid,ilay) + l0) + IF (N2(igrid,ilay) .GT. 0.) THEN + lstrat=clmix*sqrt(tke(igrid,ilay))/sqrt(N2(igrid,ilay)) + IF (lstrat .LT. l_exchange(igrid,ilay)) THEN + l_exchange(igrid,ilay)=max(lstrat,lmin) + switch_num(igrid,ilay)=.true. + ENDIF + ENDIF + ENDDO + ENDDO + +ELSE +! default: neglect effect of local stratification + + DO ilay=2,nlay+1 + DO igrid=1,ngrid + l_exchange(igrid,ilay) = kappa*l0*z_interf(igrid,ilay) / (kappa*z_interf(igrid,ilay) + l0) + ENDDO + + ENDDO +ENDIF + + ! Computing the TKE: !=================== IF (iflag_atke == 0) THEN -! stationary solution neglecting the vertical transport of TKE by turbulence + DO ilay=2,nlay + DO igrid=1,ngrid + tke(igrid,ilay) = cepsilon * l_exchange(igrid,ilay)**2 * Sm(igrid,ilay) * & + shear2(igrid,ilay) * (1. - Ri(igrid,ilay) / Prandtl(igrid,ilay)) + ENDDO + ENDDO + +ELSE IF (iflag_atke == 1) THEN + +! full implicit scheme resolved with a second order polynomial equation + DO ilay=2,nlay DO igrid=1,ngrid - tke(igrid,ilay) = cepsilon * l_exchange(igrid,ilay)**2 * Sm(igrid,ilay) * & - (((wind_u(igrid,ilay) - wind_u(igrid,ilay-1)) / dz_lay(igrid,ilay))**2 + & - ((wind_v(igrid,ilay) - wind_v(igrid,ilay-1)) / dz_lay(igrid,ilay))**2 ) * & - (1. - Ri(igrid,ilay) / Prandtl(igrid,ilay)) - ENDDO - ENDDO - -ELSE ! TO DO - + qq=sqrt(2.*tke(igrid,ilay)) + delta=1.+4.*dtime/cepsilon/l_exchange(igrid,ilay)/(2.**(3/2)) * & + (qq+dtime*l_exchange(igrid,ilay)/sqrt(2.)*Sm(igrid,ilay)*shear2(igrid,ilay) & + *(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay))) + qq=(-1. + sqrt(delta))/dtime*cepsilon*sqrt(2.)*l_exchange(igrid,ilay) + tke(igrid,ilay)=0.5*(qq**2) + ENDDO + ENDDO + +ELSE IF (iflag_atke == 2) THEN + +! full implicit scheme resolved with a second order polynomial equation +! + exact resolution for very stable cases + DO ilay=2,nlay + DO igrid=1,ngrid + qq=sqrt(2.*tke(igrid,ilay)) + IF (switch_num(igrid,ilay)) THEN + taustrat=clmix/2./sqrt(N2(igrid,ilay))*Sm(igrid,ilay)*shear2(igrid,ilay)*(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay)) & + - sqrt(N2(igrid,ilay))/2./cepsilon/clmix + taustrat=-1./taustrat + qq=qq*exp(-dtime/taustrat) + ELSE + delta=1.+4.*dtime/cepsilon/l_exchange(igrid,ilay)/(2.**(3/2)) * & + (qq+dtime*l_exchange(igrid,ilay)/sqrt(2.)*Sm(igrid,ilay) * shear2(igrid,ilay) & + *(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay))) + qq=(-1. + sqrt(delta))/dtime*cepsilon*sqrt(2.)*l_exchange(igrid,ilay) + ENDIF + tke(igrid,ilay)=0.5*(qq**2) + ENDDO + ENDDO + + + +ELSE IF (iflag_atke == 3) THEN + +! semi implicit scheme when l does not depend on tke +! positive-guaranteed if pr slope in stable condition >1 + DO ilay=2,nlay + DO igrid=1,ngrid + qq=sqrt(2.*tke(igrid,ilay)) + qq=(qq+dtime*l_exchange(igrid,ilay)/sqrt(2.)*Sm(igrid,ilay)*(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay))) & + /(1.+dtime*qq/cepsilon/l_exchange(igrid,ilay)/(2.**(3/2))) + tke(igrid,ilay)=0.5*(qq**2) + ENDDO + ENDDO + + +ELSE call abort_physic("atke_compute_km_kh", & - 'traitement non-stationnaire de la tke pas encore prevu', 1) + 'numerical treatment of TKE not possible yet', 1) END IF Index: LMDZ6/trunk/libf/phylmd/atke_turbulence_ini_mod.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/atke_turbulence_ini_mod.F90 (revision 4630) +++ LMDZ6/trunk/libf/phylmd/atke_turbulence_ini_mod.F90 (revision 4631) @@ -5,10 +5,10 @@ save -integer :: iflag_atke, iflag_num_atke -!$OMP THREADPRIVATE(iflag_atke, iflag_num_atke) +integer :: iflag_atke, iflag_num_atke, iflag_atke_lmix + !$OMP THREADPRIVATE(iflag_atke, iflag_num_atke, iflag_atke_lmix) real :: kappa = 0.4 ! Von Karman constant !$OMP THREADPRIVATE(kappa) - real :: l0, ric, ri0, cinf, cepsilon, pr_slope, pr_asym, pr_neut - !$OMP THREADPRIVATE(l0, ric, cinf, cepsilon, pr_slope, pr_asym, pr_neut) + real :: l0, ric, ri0, cinf, cepsilon, pr_slope, pr_asym, pr_neut, clmix + !$OMP THREADPRIVATE(l0, ric, cinf, cepsilon, pr_slope, pr_asym, pr_neut, clmix) integer :: lunout,prt_level !$OMP THREADPRIVATE(lunout,prt_level) @@ -19,4 +19,6 @@ !$OMP THREADPRIVATE(viscom,viscoh) + real :: lmin=0.001 ! minimum mixing length + !$OMP THREADPRIVATE(lmin) @@ -45,10 +47,23 @@ CALL getin_p('iflag_atke',iflag_atke) + ! flag that controls the calculation of mixing length in atke + iflag_atke_lmix=0 + CALL getin_p('iflag_atke_lmix',iflag_atke_lmix) + + if (iflag_atke .eq. 0 .and. iflag_atke_lmix>0) then + call abort_physic("atke_turbulence_ini", & + 'stationary scheme must use mixing length formulation not depending on tke', 1) + endif + + ! flag that controls the numerical treatment of diffusion coeffiient calculation iflag_num_atke=0 CALL getin_p('iflag_num_atke',iflag_num_atke) - ! asymptotic mixing length [m] - l0=150.0 + ! asymptotic mixing length in neutral conditions [m] + ! Sun et al 2011, JAMC + ! between 10 and 40 + + l0=15.0 CALL getin_p('atke_l0',l0) @@ -62,5 +77,5 @@ ! constant for tke dissipation calculation - cepsilon=16.6 ! default value as in yamada4 + cepsilon=16.6/2./sqrt(2.) ! default value as in yamada4 CALL getin_p('atke_cepsilon',cepsilon) @@ -68,4 +83,8 @@ pr_slope=5.0 ! default value from Zilitinkevich et al. 2005 CALL getin_p('atke_pr_slope',pr_slope) + if (pr_slope .le. 1) then + call abort_physic("atke_turbulence_ini", & + 'pr_slope has to be greater than 1 for consistency of the tke scheme', 1) + endif ! asymptotic turbulent prandt number value for Ri=-Inf @@ -77,4 +96,7 @@ CALL getin_p('atke_pr_neut',pr_neut) + ! coefficient for mixing length depending on local stratification + clmix=0.2 + CALL getin_p('atke_clmix',clmix) RETURN Index: LMDZ6/trunk/libf/phylmd/call_atke_mod.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/call_atke_mod.F90 (revision 4630) +++ LMDZ6/trunk/libf/phylmd/call_atke_mod.F90 (revision 4631) @@ -54,5 +54,5 @@ -call atke_compute_km_kh(ngrid,nlay, & +call atke_compute_km_kh(ngrid,nlay,dtime,& wind_u,wind_v,temp,play,pinterf, & tke,Km_out,Kh_out) @@ -70,5 +70,5 @@ - call atke_compute_km_kh(ngrid,nlay, & + call atke_compute_km_kh(ngrid,nlay,dtime,& wind_u_predict,wind_v_predict,temp,play,pinterf, & tke,Km_out,Kh_out)