module atke_exchange_coeff_mod implicit none contains subroutine atke_compute_km_kh(ngrid,nlay, & wind_u,wind_v,temp,play,pinterf, & tke,Km_out,Kh_out) !======================================================================== ! Routine that computes turbulent Km / Kh coefficients with a ! 1.5 order closure scheme (TKE) with or without stationarity assumption ! ! This parameterization has been constructed in the framework of a ! collective and collaborative workshop, ! the so-called 'Atelier TKE (ATKE)' with ! K. Arjdal, L. Raillard, C. Dehondt, P. Tiengou, A. Spiga, F. Cheruy, T Dubos, ! M. Coulon-Decorzens, S. Fromang, G. Riviere, A. Sima, F. Hourdin, E. Vignon ! ! Main assumptions of the model : ! (1) dry atmosphere ! (2) horizontal homogeneity (Dx=Dy=0.) !======================================================================= 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 implicit none ! Declarations: !============= INTEGER, INTENT(IN) :: ngrid ! number of horizontal index (flat grid) INTEGER, INTENT(IN) :: nlay ! number of vertical index REAL, DIMENSION(ngrid,nlay), INTENT(IN) :: wind_u ! zonal velocity (m/s) REAL, DIMENSION(ngrid,nlay), INTENT(IN) :: wind_v ! meridional velocity (m/s) REAL, DIMENSION(ngrid,nlay), INTENT(IN) :: temp ! temperature (K) REAL, DIMENSION(ngrid,nlay), INTENT(IN) :: play ! pressure (Pa) REAL, DIMENSION(ngrid,nlay+1), INTENT(IN) :: pinterf ! pressure at interfaces(Pa) REAL, DIMENSION(ngrid,nlay+1), INTENT(INOUT) :: tke ! turbulent kinetic energy at interface between layers REAL, DIMENSION(ngrid,nlay), INTENT(OUT) :: Km_out ! output: Exchange coefficient for momentum at interface between layers REAL, DIMENSION(ngrid,nlay), INTENT(OUT) :: Kh_out ! output: Exchange coefficient for heat flux at interface between layers ! Local variables REAL, DIMENSION(ngrid,nlay+1) :: Km ! Exchange coefficient for momentum at interface between layers REAL, DIMENSION(ngrid,nlay+1) :: Kh ! Exchange coefficient for heat flux at interface between layers REAL, DIMENSION(ngrid,nlay) :: theta ! Potential temperature REAL, DIMENSION(ngrid,nlay+1) :: l_exchange ! Length of exchange (at interface) REAL, DIMENSION(ngrid,nlay+1) :: z_interf ! Altitude at the interface REAL, DIMENSION(ngrid,nlay) :: z_lay ! Altitude of layers 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) :: 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) INTEGER :: igrid,ilay ! horizontal,vertical index (flat grid) REAL :: cn,Ri0,Ri1 ! parameter for Sm stability function and Prandlt REAL :: preff ! reference pressure for potential temperature calculations REAL :: thetam ! mean potential temperature at interface ! Initializations: !================ DO igrid=1,ngrid dz_interf(igrid,1) = 0.0 z_interf(igrid,1) = 0.0 END DO ! Calculation of potential temperature: (if vapor -> todo virtual potential temperature) !===================================== preff=100000. ! The result should not depend on the choice of preff DO ilay=1,nlay DO igrid = 1, ngrid theta(igrid,ilay)=temp(igrid,ilay)*(preff/play(igrid,ilay))**(rd/rcpd) END DO END DO ! Calculation of altitude of layers' middle and bottom interfaces: !================================================================= DO ilay=2,nlay+1 DO igrid=1,ngrid dz_interf(igrid,ilay-1) = rd*temp(igrid,ilay-1)/rg/play(igrid,ilay-1)*(pinterf(igrid,ilay-1)-pinterf(igrid,ilay)) z_interf(igrid,ilay) = z_interf(igrid,ilay-1) + dz_interf(igrid,ilay-1) ENDDO ENDDO DO ilay=1,nlay DO igrid=1,ngrid z_lay(igrid,ilay)=0.5*(z_interf(igrid, ilay+1) + z_interf(igrid, ilay)) ENDDO ENDDO ! 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: !=================================================================== ! calculation of cn = Sm value at Ri=0 ! direct dependance on cepsilon to guarantee Fm=1 (first-order like stability function) at Ri=0 cn=(1./sqrt(cepsilon))**(2/3) ! calculation of Ri0 such that continuity in slope of Sm at Ri=0 Ri0=2./rpi*(cinf - cn)*ric/cn ! calculation of Ri1 to guarantee continuity in slope of Prandlt number at Ri=0 Ri1 = -2./rpi * (pr_asym - pr_neut) / pr_slope DO ilay=2,nlay DO igrid=1,ngrid 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) IF (Ri(igrid,ilay) < 0.) THEN ! unstable cases Sm(igrid,ilay) = 2./rpi * (cinf-cn) * atan(-Ri(igrid,ilay)/Ri0) + cn Prandtl(igrid,ilay) = -2./rpi * (pr_asym - pr_neut) * atan(Ri(igrid,ilay)/Ri1) + pr_neut ELSE ! stable cases Sm(igrid,ilay) = max(0.,cn*(1.-Ri(igrid,ilay)/Ric)) Prandtl(igrid,ilay) = pr_neut + Ri(igrid,ilay) * pr_slope IF (Ri(igrid,ilay) .GE. Prandtl(igrid,ilay)) THEN call abort_physic("atke_compute_km_kh", & 'Ri>=Pr in stable conditions -> violates energy conservation principles, change pr_neut or slope', 1) ENDIF END IF Sh(igrid,ilay) = Sm(igrid,ilay) / Prandtl(igrid,ilay) ENDDO ENDDO ! 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) * & (((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 call abort_physic("atke_compute_km_kh", & 'traitement non-stationnaire de la tke pas encore prevu', 1) END IF ! Computing eddy diffusivity coefficients: !======================================== DO ilay=2,nlay ! TODO: also calculate for nlay+1 ? DO igrid=1,ngrid ! we add the molecular viscosity to Km,h Km(igrid,ilay) = viscom + l_exchange(igrid,ilay) * Sm(igrid,ilay) * tke(igrid,ilay)**0.5 Kh(igrid,ilay) = viscoh + l_exchange(igrid,ilay) * Sh(igrid,ilay) * tke(igrid,ilay)**0.5 END DO END DO ! for output: !=========== Km_out(1:ngrid,2:nlay)=Km(1:ngrid,2:nlay) Kh_out(1:ngrid,2:nlay)=Kh(1:ngrid,2:nlay) end subroutine atke_compute_km_kh end module atke_exchange_coeff_mod