module lmdz_atke_exchange_coeff implicit none contains subroutine atke_compute_km_kh(ngrid,nlay,dtime, & wind_u,wind_v,temp,qvap,play,pinterf,cdrag_uv, & tke,eps,tke_shear,tke_buoy,tke_trans,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) horizontal homogeneity (Dx=Dy=0.) !======================================================================= USE lmdz_atke_turbulence_ini, ONLY : iflag_atke, kappa, l0, ric, cinf, rpi, rcpd, atke_ok_virtual, ri0, ri1 USE lmdz_atke_turbulence_ini, ONLY : cepsilon, pr_slope, pr_asym, pr_neut, ctkes, rg, rd, rv, atke_ok_vdiff USE lmdz_atke_turbulence_ini, ONLY : viscom, viscoh, clmix, clmixshear, iflag_atke_lmix, lmin, smmin, cn !!------------------------------------------------------------------------------------------------------------- ! integer :: iflag_atke ! flag that controls options in atke_compute_km_kh ! integer :: iflag_atke_lmix ! flag that controls the calculation of mixing length in atke ! integer :: iflag_num_atke ! flag that controls the numerical treatment of diffusion coeffiient calculation ! logical :: atke_ok_vdiff ! activate vertical diffusion of TKE or not ! logical :: atke_ok_virtual ! account for vapor for flottability ! real :: kappa = 0.4 ! Von Karman constant ! real :: cn ! Sm value at Ri=0 ! real :: cinf ! Sm value at Ri=-Inf ! real :: ri0 ! Richardson number near zero to guarantee continuity in slope of Sm (stability function) at Ri=0 ! real :: ri1 ! Richardson number near zero to guarantee continuity in slope of Pr (Prandlt's number) at Ri=0 ! real :: lmin ! minimum mixing length corresponding to the Kolmogov dissipation scale in planetary atmospheres (Chen et al 2016, JGR Atmos) ! real :: ctkes ! coefficient for surface TKE ! real :: clmixshear ! coefficient for mixing length depending on local wind shear ! Tunable parameters for the ATKE scheme and their range of values !!------------------------------------------------------------------------------------------------------------- ! real :: cepsilon ! controls the value of the dissipation length scale, range [1.2 - 10] ! real :: cke ! controls the value of the diffusion coefficient of TKE, range [1 - 5] ! real :: l0 ! asymptotic mixing length far from the ground [m] (Sun et al 2011, JAMC), range [15 - 75] ! real :: clmix ! controls the value of the mixing length in stratified conditions, range [0.1 - 2] ! real :: ric ! critical Richardson number controlling the slope of Sm in stable conditions, range [0.19 - 0.25] ! real :: smmin ! minimum value of Sm in very stable conditions (defined here as minsqrt(Ez/Ek)) at large Ri, range [0.025 - 0.1] ! real :: pr_neut ! neutral value of the Prandtl number (Ri=0), range [0.7 - 1] ! real :: pr_slope ! linear slope of Pr with Ri in the very stable regime, range [3 - 5] ! real :: cinffac ! ratio between cinf and cn controlling the convective limit of Sm, range [1.2 - 5.0] ! real :: pr_asym ! value of Prandlt in the convective limit(Ri=-Inf), range [0.3 - 0.5] !!------------------------------------------------------------------------------------------------------------- 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) :: qvap ! specific humidity (kg/kg) REAL, DIMENSION(ngrid,nlay), INTENT(IN) :: play ! pressure (Pa) REAL, DIMENSION(ngrid,nlay+1), INTENT(IN) :: pinterf ! pressure at interfaces(Pa) REAL, DIMENSION(ngrid), INTENT(IN) :: cdrag_uv ! surface drag coefficient for momentum REAL, DIMENSION(ngrid,nlay+1), INTENT(INOUT) :: tke ! turbulent kinetic energy at interface between layers (m2/s2) REAL, DIMENSION(ngrid,nlay+1), INTENT(OUT) :: eps ! output: TKE dissipation rate at interface between layers (m2/s3) REAL, DIMENSION(ngrid,nlay+1), INTENT(OUT) :: tke_shear! output: TKE shear production rate (m2/s3) REAL, DIMENSION(ngrid,nlay+1), INTENT(OUT) :: tke_buoy ! output: TKE buoyancy production rate (m2/s3) REAL, DIMENSION(ngrid,nlay+1), INTENT(OUT) :: tke_trans! output: TKE transport (diffusion) term (m2/s3) REAL, DIMENSION(ngrid,nlay), INTENT(OUT) :: Km_out ! output: Exchange coefficient for momentum at interface between layers (m2/s) REAL, DIMENSION(ngrid,nlay), INTENT(OUT) :: Kh_out ! output: Exchange coefficient for heat flux at interface between layers (m2/s) ! 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) INTEGER :: igrid,ilay ! horizontal,vertical index (flat grid) 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 REAL :: netloss ! net loss term of tke REAL :: netsource ! net source term of tke REAL :: ustar ! friction velocity estimation REAL :: invtau REAL :: rvap ! 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 -> virtual potential temperature) !===================================== preff=100000. ! results 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 ! account for water vapor mass for buoyancy calculation IF (atke_ok_virtual) THEN DO ilay=1,nlay DO igrid = 1, ngrid rvap=max(0.,qvap(igrid,ilay)/(1.-qvap(igrid,ilay))) theta(igrid,ilay)=theta(igrid,ilay)*(1.+rvap/(RD/RV))/(1.+rvap) END DO END DO ENDIF ! 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: !=================================================================== 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(smmin,cn*(1.-Ri(igrid,ilay)/Ric)) ! prandlt expression from venayagamoorthy and stretch 2010, Li et al 2019 Prandtl(igrid,ilay) = pr_neut*exp(-pr_slope/pr_neut*Ri(igrid,ilay)+Ri(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: !============================================================== IF (iflag_atke_lmix .EQ. 1 ) THEN ! Blackadar formulation (~kappa l) + buoyancy length scale (Deardoff 1980) for very stable conditions 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)) lstrat=max(lstrat,lmin) !Inverse interpolation, Van de Wiel et al. 2010 l_exchange(igrid,ilay)=(1./(l_exchange(igrid,ilay))+1./(lstrat))**(-1.0) ENDIF ENDDO ENDDO ELSE IF (iflag_atke_lmix .EQ. 2 ) THEN ! add effect of wind shear on lstrat following grisogono and belusic 2008, qjrms ! implies 2 tuning coefficients clmix and clmixshear 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. .AND. shear2(igrid,ilay) .GT. 0.) THEN lstrat=min(clmix*sqrt(tke(igrid,ilay))/sqrt(N2(igrid,ilay)), & clmixshear*sqrt(tke(igrid,ilay))/sqrt(shear2(igrid,ilay))) lstrat=max(lstrat,lmin) !Inverse interpolation, Van de Wiel et al. 2010 l_exchange(igrid,ilay)=(1./(l_exchange(igrid,ilay))+1./(lstrat))**(-1.0) ENDIF ENDDO ENDDO ELSE IF (iflag_atke_lmix .EQ. 3 ) THEN ! add effect of wind shear on lstrat following grisogono 2010, qjrms ! keeping a single tuning coefficient clmix 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. .AND. shear2(igrid,ilay) .GT. 0.) THEN lstrat=clmix*sqrt(tke(igrid,ilay))/sqrt(shear2(igrid,ilay))*(1.0+Ri(igrid,ilay)/(2.*Prandtl(igrid,ilay))) lstrat=max(lstrat,lmin) !Inverse interpolation, Van de Wiel et al. 2010 l_exchange(igrid,ilay)=(1./(l_exchange(igrid,ilay))+1./(lstrat))**(-1.0) ENDIF ENDDO ENDDO ELSE ! default Blackadar formulation: 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 k>=2: !======================== IF (iflag_atke == 0) THEN ! stationary solution (dtke/dt=0) 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)) eps(igrid,ilay) = (tke(igrid,ilay)**(3./2))/(cepsilon*l_exchange(igrid,ilay)) tke_shear(igrid,ilay)=l_exchange(igrid,ilay)*Sm(igrid,ilay)*sqrt(tke(igrid,ilay))*shear2(igrid,ilay) tke_buoy(igrid,ilay)=-l_exchange(igrid,ilay)*Sm(igrid,ilay)*sqrt(tke(igrid,ilay))*shear2(igrid,ilay) & *(Ri(igrid,ilay) / Prandtl(igrid,ilay)) ENDDO ENDDO ELSE IF (iflag_atke == 1) THEN ! full implicit scheme resolved with a second order polynomial equation ! default solution which shows fair convergence properties DO ilay=2,nlay DO igrid=1,ngrid qq=max(sqrt(2.*tke(igrid,ilay)),1.e-10) delta=(2.*sqrt(2.)*cepsilon*l_exchange(igrid,ilay)/dtime)**2. & +4.*(2.*sqrt(2.)*cepsilon*l_exchange(igrid,ilay)/dtime*qq + & 2.*l_exchange(igrid,ilay)*l_exchange(igrid,ilay)*cepsilon*Sm(igrid,ilay) & *shear2(igrid,ilay) * (1. - Ri(igrid,ilay) / Prandtl(igrid,ilay))) qq=(-2.*sqrt(2.)*cepsilon*l_exchange(igrid,ilay)/dtime + sqrt(delta))/2. qq=max(0.,qq) tke(igrid,ilay)=0.5*(qq**2) eps(igrid,ilay) = (tke(igrid,ilay)**(3./2))/(cepsilon*l_exchange(igrid,ilay)) tke_shear(igrid,ilay)=l_exchange(igrid,ilay)*Sm(igrid,ilay)*sqrt(tke(igrid,ilay))*shear2(igrid,ilay) tke_buoy(igrid,ilay)=-l_exchange(igrid,ilay)*Sm(igrid,ilay)*sqrt(tke(igrid,ilay))*shear2(igrid,ilay) & *(Ri(igrid,ilay) / Prandtl(igrid,ilay)) ENDDO ENDDO ELSE IF (iflag_atke == 2) 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=max(sqrt(2.*tke(igrid,ilay)),1.e-10) qq=(qq+l_exchange(igrid,ilay)*Sm(igrid,ilay)*dtime/sqrt(2.) & *shear2(igrid,ilay)*(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay))) & /(1.+qq*dtime/(cepsilon*l_exchange(igrid,ilay)*2.*sqrt(2.))) tke_shear(igrid,ilay)=l_exchange(igrid,ilay)*Sm(igrid,ilay)*qq/sqrt(2.)*shear2(igrid,ilay) tke_buoy(igrid,ilay)=-l_exchange(igrid,ilay)*Sm(igrid,ilay)*qq/sqrt(2.)*shear2(igrid,ilay) & *(Ri(igrid,ilay) / Prandtl(igrid,ilay)) tke(igrid,ilay)=0.5*(qq**2) eps(igrid,ilay) = (tke(igrid,ilay)**(3./2))/(cepsilon*l_exchange(igrid,ilay)) ENDDO ENDDO ELSE IF (iflag_atke == 3) THEN ! numerical resolution adapted from that in MAR (Deleersnijder 1992) ! positively defined by construction DO ilay=2,nlay DO igrid=1,ngrid eps(igrid,ilay) = (tke(igrid,ilay)**(3./2))/(cepsilon*l_exchange(igrid,ilay)) qq=max(sqrt(2.*tke(igrid,ilay)),1.e-10) tke_shear(igrid,ilay)=l_exchange(igrid,ilay)*Sm(igrid,ilay)*qq/sqrt(2.)*shear2(igrid,ilay) tke_buoy(igrid,ilay)=-l_exchange(igrid,ilay)*Sm(igrid,ilay)*qq/sqrt(2.)*shear2(igrid,ilay) & *(Ri(igrid,ilay) / Prandtl(igrid,ilay)) IF (Ri(igrid,ilay) .LT. 0.) THEN netloss=qq/(2.*sqrt(2.)*cepsilon*l_exchange(igrid,ilay)) netsource=l_exchange(igrid,ilay)*Sm(igrid,ilay)/sqrt(2.)*shear2(igrid,ilay)*(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay)) ELSE netloss=qq/(2.*sqrt(2.)*cepsilon*l_exchange(igrid,ilay))+ & l_exchange(igrid,ilay)*Sm(igrid,ilay)/sqrt(2.)*N2(igrid,ilay)/Prandtl(igrid,ilay) netsource=l_exchange(igrid,ilay)*Sm(igrid,ilay)/sqrt(2.)*shear2(igrid,ilay) ENDIF qq=((qq**2)/dtime+qq*netsource)/(qq/dtime+netloss) tke(igrid,ilay)=0.5*(qq**2) ENDDO ENDDO ELSE IF (iflag_atke == 4) THEN ! semi implicit scheme from Arpege (V. Masson methodology with ! Taylor expansion of the dissipation term) DO ilay=2,nlay DO igrid=1,ngrid qq=max(sqrt(2.*tke(igrid,ilay)),1.e-10) tke_shear(igrid,ilay)=l_exchange(igrid,ilay)*Sm(igrid,ilay)*qq/sqrt(2.)*shear2(igrid,ilay) tke_buoy(igrid,ilay)=-l_exchange(igrid,ilay)*Sm(igrid,ilay)*qq/sqrt(2.)*shear2(igrid,ilay) & *(Ri(igrid,ilay) / Prandtl(igrid,ilay)) qq=(l_exchange(igrid,ilay)*Sm(igrid,ilay)/sqrt(2.)*shear2(igrid,ilay)*(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay)) & +qq*(1.+dtime*qq/(cepsilon*l_exchange(igrid,ilay)*2.*sqrt(2.)))) & /(1.+2.*qq*dtime/(cepsilon*l_exchange(igrid,ilay)*2.*sqrt(2.))) qq=max(0.,qq) tke(igrid,ilay)=0.5*(qq**2) eps(igrid,ilay) = (tke(igrid,ilay)**(3./2))/(cepsilon*l_exchange(igrid,ilay)) ENDDO ENDDO ELSE call abort_physic("atke_compute_km_kh", & 'numerical treatment of TKE not possible yet', 1) END IF ! We impose a 0 tke at nlay+1 !============================== DO igrid=1,ngrid tke(igrid,nlay+1)=0. eps(igrid,nlay+1)=0. tke_shear(igrid,nlay+1)=0. tke_buoy(igrid,nlay+1)=0. END DO ! Calculation of surface TKE (k=1) !================================= ! surface TKE calculation inspired from what is done in Arpege (see E. Bazile note) DO igrid=1,ngrid ustar=sqrt(cdrag_uv(igrid)*(wind_u(igrid,1)**2+wind_v(igrid,1)**2)) tke(igrid,1)=ctkes*(ustar**2) eps(igrid,1)=0. ! arbitrary as TKE is not properly defined at the surface tke_shear(igrid,1)=0. tke_buoy(igrid,1)=0. END DO ! vertical diffusion of TKE !========================== tke_trans(:,:)=0. IF (atke_ok_vdiff) THEN CALL atke_vdiff_tke(ngrid,nlay,dtime,z_lay,z_interf,temp,play,l_exchange,Sm,tke,tke_trans) ENDIF ! 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 !=============================================================================================== subroutine atke_vdiff_tke(ngrid,nlay,dtime,z_lay,z_interf,temp,play,l_exchange,Sm,tke,tke_trans) ! routine that computes the vertical diffusion of TKE by the turbulence ! using an implicit resolution (See note by Dufresne and Ghattas (2009)) ! E Vignon, July 2023 USE lmdz_atke_turbulence_ini, ONLY : rd, cke, viscom 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) :: z_lay ! altitude of mid-layers (m) REAL, DIMENSION(ngrid,nlay+1), INTENT(IN) :: z_interf ! altitude of bottom interfaces (m) 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) :: l_exchange ! mixing length at interfaces between layers REAL, DIMENSION(ngrid,nlay+1), INTENT(IN) :: Sm ! stability function for eddy diffusivity for momentum at interface between layers REAL, DIMENSION(ngrid,nlay+1), INTENT(INOUT) :: tke ! turbulent kinetic energy at interface between layers REAL, DIMENSION(ngrid,nlay+1), INTENT(INOUT) :: tke_trans ! turbulent kinetic energy transport term (m2/s3) INTEGER :: igrid,ilay REAL, DIMENSION(ngrid,nlay+1) :: Ke ! eddy diffusivity for TKE REAL, DIMENSION(ngrid,nlay+1) :: dtke REAL, DIMENSION(ngrid,nlay+1) :: ak, bk, ck, CCK, DDK REAL :: gammak,Kem,KKb,KKt ! Few initialisations CCK(:,:)=0. DDK(:,:)=0. dtke(:,:)=0. ! Eddy diffusivity for TKE DO ilay=2,nlay DO igrid=1,ngrid Ke(igrid,ilay)=(viscom+l_exchange(igrid,ilay)*Sm(igrid,ilay)*sqrt(tke(igrid,ilay)))*cke ENDDO ENDDO ! at the top of the atmosphere set to 0 Ke(:,nlay+1)=0. ! at the surface, set it equal to that at the first model level Ke(:,1)=Ke(:,2) ! calculate intermediary variables DO ilay=2,nlay DO igrid=1,ngrid Kem=0.5*(Ke(igrid,ilay+1)+Ke(igrid,ilay)) KKt=Kem*play(igrid,ilay)/rd/temp(igrid,ilay)/(z_interf(igrid,ilay+1)-z_interf(igrid,ilay)) Kem=0.5*(Ke(igrid,ilay)+Ke(igrid,ilay-1)) KKb=Kem*play(igrid,ilay-1)/rd/temp(igrid,ilay-1)/(z_interf(igrid,ilay)-z_interf(igrid,ilay-1)) gammak=1./(z_lay(igrid,ilay)-z_lay(igrid,ilay-1)) ak(igrid,ilay)=-gammak*dtime*KKb ck(igrid,ilay)=-gammak*dtime*KKt bk(igrid,ilay)=1.+gammak*dtime*(KKt+KKb) ENDDO ENDDO ! calculate CCK and DDK coefficients ! downhill phase DO igrid=1,ngrid CCK(igrid,nlay)=tke(igrid,nlay)/bk(igrid,nlay) DDK(igrid,nlay)=-ak(igrid,nlay)/bk(igrid,nlay) ENDDO DO ilay=nlay-1,2,-1 DO igrid=1,ngrid CCK(igrid,ilay)=(tke(igrid,ilay)/bk(igrid,ilay)-ck(igrid,ilay)/bk(igrid,ilay)*CCK(igrid,ilay+1)) & / (1.+ck(igrid,ilay)/bk(igrid,ilay)*DDK(igrid,ilay+1)) DDK(igrid,ilay)=-ak(igrid,ilay)/bk(igrid,ilay)/(1+ck(igrid,ilay)/bk(igrid,ilay)*DDK(igrid,ilay+1)) ENDDO ENDDO ! calculate TKE ! uphill phase DO ilay=2,nlay+1 DO igrid=1,ngrid dtke(igrid,ilay)=CCK(igrid,ilay)+DDK(igrid,ilay)*tke(igrid,ilay-1)-tke(igrid,ilay) ENDDO ENDDO ! update TKE tke(:,:)=tke(:,:)+dtke(:,:) tke_trans(:,:)=dtke(:,:)/dtime end subroutine atke_vdiff_tke end module lmdz_atke_exchange_coeff