module atke_exchange_coeff_mod implicit none contains subroutine atke_compute_km_kh(ngrid,nlay,dtime, & 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, clmix, iflag_atke_lmix, lmin implicit none ! Declarations: !============= INTEGER, INTENT(IN) :: ngrid ! number of horizontal index (flat grid) 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) 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) :: 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) 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 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: !================ 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 ! 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)) 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 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 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 IF (iflag_atke_lmix .EQ. 2 ) THEN ! add effect of wind shear on lstrat following grisogono and belusic 2008, qjrms, eq 2 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))/(2.*max(sqrt(shear2(igrid,ilay)),1E-10)*(1.+sqrt(Ri(igrid,ilay))/2.)) IF (lstrat .LT. l_exchange(igrid,ilay)) THEN l_exchange(igrid,ilay)=max(lstrat,lmin) ENDIF ENDIF ENDDO ENDDO ELSE ! default: neglect effect of local stratification and shear 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 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 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 (iflag_atke_lmix=1) 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", & 'numerical treatment of TKE not possible yet', 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