Index: trunk/LMDZ.MARS/libf/phymars/lmdz_atke_exchange_coeff.F90 =================================================================== --- trunk/LMDZ.MARS/libf/phymars/lmdz_atke_exchange_coeff.F90 (revision 3168) +++ trunk/LMDZ.MARS/libf/phymars/lmdz_atke_exchange_coeff.F90 (revision 3168) @@ -0,0 +1,450 @@ +module lmdz_atke_exchange_coeff + +implicit none + +contains + +subroutine atke_compute_km_kh(ngrid,nlay,dtime, & + wind_u,wind_v,temp,qvap,play,pinterf,z_lay,z_interf,cdrag_uv, & + 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 +! +! Transposed in the Mars PCM by Lucas Lange +! +! 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 + +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,nlay), INTENT(IN) :: z_lay ! altitude at the middle of the vertical layers (m) +REAL, DIMENSION(ngrid,nlay+1), INTENT(IN) :: z_interf ! altitude at interfaces (m) +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 + +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) :: 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 +END DO + +! Calculation of potential temperature: (if vapor -> virtual potential temperature) +!===================================== + +preff=610. +! 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 the layer's distance at the middle and bottom interfaces: +!================================================================= + +dz_interf(:,1) = 0. ! first layer +dz_lay(:,1) = 0. +DO ilay = 2, nlay + DO igrid=1,ngrid + dz_lay(igrid,ilay)=z_lay(igrid,ilay)-z_lay(igrid,ilay-1) + dz_interf(igrid,ilay)=z_interf(igrid,ilay)-z_interf(igrid,ilay-1) + ENDDO +END DO + + +! 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)) + 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) + 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(igrid,ilay)=0.5*(qq**2) + 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 + qq=max(sqrt(2.*tke(igrid,ilay)),1.e-10) + 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) + 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) + 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. +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) +END DO + + +! vertical diffusion of TKE +!========================== +IF (atke_ok_vdiff) THEN + CALL atke_vdiff_tke(ngrid,nlay,dtime,z_lay,z_interf,temp,play,l_exchange,Sm,tke) +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) + +! 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 + + + +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(:,:) + + +end subroutine atke_vdiff_tke + + + +end module lmdz_atke_exchange_coeff Index: trunk/LMDZ.MARS/libf/phymars/lmdz_atke_turbulence_ini.F90 =================================================================== --- trunk/LMDZ.MARS/libf/phymars/lmdz_atke_turbulence_ini.F90 (revision 3168) +++ trunk/LMDZ.MARS/libf/phymars/lmdz_atke_turbulence_ini.F90 (revision 3168) @@ -0,0 +1,177 @@ +MODULE lmdz_atke_turbulence_ini + +implicit none + +! declaration of constants and parameters +save + +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 :: cinffac + !$OMP THREADPRIVATE(cinffac) + real :: l0,ric,ri0,ri1,cinf,cn,cepsilon,pr_slope,pr_asym,pr_neut,clmix,clmixshear,smmin,ctkes,cke + !$OMP THREADPRIVATE(l0,ri0,ri1,ric,cinf,cn,cepsilon,pr_slope,pr_asym,pr_neut,clmix,clmixshear,smmin,ctkes,cke) + integer :: lunout,prt_level + !$OMP THREADPRIVATE(lunout,prt_level) + real :: rg, rd, rpi, rcpd, rv + !$OMP THREADPRIVATE(rg, rd, rpi, rcpd, rv) + real :: viscom, viscoh + !$OMP THREADPRIVATE(viscom,viscoh) + real :: lmin=0.01 ! minimum mixing length corresponding to the Kolmogov dissipation scale + ! in planetary atmospheres (Chen et al 2016, JGR Atmos) + !$OMP THREADPRIVATE(lmin) + logical :: atke_ok_vdiff, atke_ok_virtual + !$OMP THREADPRIVATE(atke_ok_vdiff,atke_ok_virtual) + +CONTAINS + +SUBROUTINE atke_ini(rg_in, rd_in, rpi_in, rcpd_in, rv_in, viscom_in, viscoh_in) + + USE ioipsl_getin_p_mod, ONLY : getin_p + + real, intent(in) :: rg_in, rd_in, rpi_in, rcpd_in, rv_in, viscom_in, viscoh_in + + + ! input arguments (universal constants for planet) + !------------------------------------------------- + + ! gravity acceleration + rg=rg_in + ! dry gas constant (R/M, R=perfect gas constant and M is the molar mass of the fluid) + rd=rd_in + ! Pi number + rpi=rpi_in + ! cp per unit mass of the gas + rcpd=rcpd_in + ! water vapor constant (for simulations in Earth's atmosphere) + rv=rv_in + ! kinematic molecular viscosity for momentum + viscom=viscom_in + ! kinematic molecular viscosity for heat + viscoh=viscoh_in + + + ! Read flag values in .def files + !------------------------------- + + + ! flag that controls options in atke_compute_km_kh + iflag_atke=0 + 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 + + ! activate vertical diffusion of TKE or not + atke_ok_vdiff=.false. + CALL getin_p('atke_ok_vdiff',atke_ok_vdiff) + + + ! account for vapor for flottability + atke_ok_virtual=.false. + CALL getin_p('atke_ok_virtual',atke_ok_virtual) + + + ! 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 in neutral conditions [m] + ! Sun et al 2011, JAMC + ! between 10 and 40 + l0=15.0 + CALL getin_p('atke_l0',l0) + + ! critical Richardson number + ric=0.25 + CALL getin_p('atke_ric',ric) + + + ! constant for tke dissipation calculation + cepsilon=5.87 ! default value as in yamada4 + CALL getin_p('atke_cepsilon',cepsilon) + + + ! 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) + + ! asymptotic value of Sm for Ri=-Inf + cinffac=2.0 + CALL getin_p('atke_cinffac',cinffac) + cinf=cinffac*cn + if (cinf .le. cn) then + call abort_physic("atke_turbulence_ini", & + 'cinf cannot be lower than cn', 1) + endif + + + ! coefficient for surface TKE + ! following Lenderink & Holtslag 2004, ctkes=(cepsilon)**(2/3) + ! (provided by limit condition in neutral conditions) + ctkes=(cepsilon)**(2./3.) + + ! slope of Pr=f(Ri) for stable conditions + 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 + + ! value of turbulent prandtl number in neutral conditions (Ri=0) + pr_neut=0.8 + CALL getin_p('atke_pr_neut',pr_neut) + + + ! asymptotic turbulent prandt number value for Ri=-Inf + pr_asym=0.4 + CALL getin_p('atke_pr_asym',pr_asym) + if (pr_asym .ge. pr_neut) then + call abort_physic("atke_turbulence_ini", & + 'pr_asym must be be lower than pr_neut', 1) + endif + + + + ! coefficient for mixing length depending on local stratification + clmix=0.5 + CALL getin_p('atke_clmix',clmix) + + ! coefficient for mixing length depending on local wind shear + clmixshear=0.5 + CALL getin_p('atke_clmixshear',clmixshear) + + + ! minimum anisotropy coefficient (defined here as minsqrt(Ez/Ek)) at large Ri. + ! From Zilitinkevich et al. 2013, it equals sqrt(0.03)~0.17 + smmin=0.17 + CALL getin_p('atke_smmin',smmin) + + + ! ratio between the eddy diffusivity coeff for tke wrt that for momentum + ! default value from Lenderink et al. 2004 + cke=2. + CALL getin_p('atke_cke',cke) + + + ! 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) + + + RETURN + +END SUBROUTINE atke_ini + +END MODULE lmdz_atke_turbulence_ini Index: trunk/LMDZ.MARS/libf/phymars/lmdz_call_atke.F90 =================================================================== --- trunk/LMDZ.MARS/libf/phymars/lmdz_call_atke.F90 (revision 3168) +++ trunk/LMDZ.MARS/libf/phymars/lmdz_call_atke.F90 (revision 3168) @@ -0,0 +1,162 @@ +module lmdz_call_atke + +USE lmdz_atke_exchange_coeff, ONLY : atke_compute_km_kh + +implicit none + + +contains + +subroutine call_atke(dtime,ngrid,nlay,cdrag_uv,cdrag_t,u_surf,v_surf,temp_surf, & + wind_u,wind_v,temp,qvap,play,pinterf,zlay,zinterf, & + tke,Km_out,Kh_out) + + + + +USE lmdz_atke_turbulence_ini, ONLY : iflag_num_atke, rg, rd + +implicit none + + +! Declarations: +!============= + +REAL, INTENT(IN) :: dtime ! timestep (s) +INTEGER, INTENT(IN) :: ngrid ! number of horizontal index (flat grid) +INTEGER, INTENT(IN) :: nlay ! number of vertical index + + +REAL, DIMENSION(ngrid), INTENT(IN) :: cdrag_uv ! drag coefficient for wind +REAL, DIMENSION(ngrid), INTENT(IN) :: cdrag_t ! drag coefficient for temperature + +REAL, DIMENSION(ngrid), INTENT(IN) :: u_surf ! x wind velocity at the surface +REAL, DIMENSION(ngrid), INTENT(IN) :: v_surf ! y wind velocity at the surface +REAL, DIMENSION(ngrid), INTENT(IN) :: temp_surf ! surface temperature + +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,nlay), INTENT(IN) :: zlay ! altitude at the middle of the vertical layers (m) +REAL, DIMENSION(ngrid,nlay+1), INTENT(IN) :: zinterf ! altitude at interfaces (m) + + +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 + + +REAL, DIMENSION(ngrid,nlay) :: wind_u_predict, wind_v_predict +REAL, DIMENSION(ngrid) :: wind1 +INTEGER i + + +call atke_compute_km_kh(ngrid,nlay,dtime,& + wind_u,wind_v,temp,qvap,play,pinterf,zlay,zinterf,cdrag_uv,& + tke,Km_out,Kh_out) + + + +if (iflag_num_atke .EQ. 1) then + !! In this case, we make an explicit prediction of the wind shear to calculate the tke in a + !! forward backward way + !! pay attention that the treatment of the TKE + !! has to be adapted when solving the TKE with a prognostic equation + + do i=1,ngrid + wind1(i)=sqrt(wind_u(i,1)**2+wind_v(i,1)**2) + enddo + call atke_explicit_prediction(ngrid,nlay,rg,rd,dtime,pinterf,play,temp,wind1,wind_u,Km_out,u_surf,cdrag_uv,wind_u_predict) + call atke_explicit_prediction(ngrid,nlay,rg,rd,dtime,pinterf,play,temp,wind1,wind_v,Km_out,v_surf,cdrag_uv,wind_v_predict) + + + call atke_compute_km_kh(ngrid,nlay,dtime,& + wind_u_predict,wind_v_predict,temp,qvap,play,pinterf,zlay,zinterf,cdrag_uv, & + tke,Km_out,Kh_out) + +end if + + + + + + +end subroutine call_atke + +!---------------------------------------------------------------------------------------- + +subroutine atke_explicit_prediction(ngrid,nlay,rg,rd,dtime,pinterf,play,temp,wind1,x_in,K_in,x_surf,cdrag,x_predict) + +INTEGER, INTENT(IN) :: ngrid ! number of horizontal index (flat grid) +INTEGER, INTENT(IN) :: nlay ! number of vertical index +REAL, INTENT(IN) :: rg,rd,dtime ! gravity, R dry air and timestep +REAL, DIMENSION(ngrid,nlay), INTENT(IN) :: play ! pressure middle of layers (Pa) +REAL, DIMENSION(ngrid,nlay), INTENT(IN) :: temp ! temperature (K) +REAL, DIMENSION(ngrid), INTENT(IN) :: wind1 ! wind speed first level (m/s) +REAL, DIMENSION(ngrid,nlay+1), INTENT(IN) :: pinterf ! pressure at interfaces(Pa) +REAL, DIMENSION(ngrid,nlay), INTENT(IN) :: x_in ! variable at the beginning of timestep +REAL, DIMENSION(ngrid,nlay+1), INTENT(IN) :: K_in ! eddy diffusivity coef at the beginning of time step +REAL, DIMENSION(ngrid), INTENT(IN) :: x_surf ! surface variable +REAL, DIMENSION(ngrid), INTENT(IN) :: cdrag ! drag coefficient +REAL, DIMENSION(ngrid,nlay), INTENT(OUT) :: x_predict ! variable at the end of time step after explicit prediction + + +integer i,k +real ml,F1,rho +real, dimension(ngrid) :: play1,temp1 +real, dimension(ngrid,nlay+1) :: K_big + +! computation of K_big + +play1(:)=play(:,1) +temp1(:)=temp(:,1) + +! "big K" calculation +do k=2,nlay-1 + do i=1,ngrid + rho=pinterf(i,k)/rd/(0.5*(temp(i,k-1)+temp(i,k))) + K_big(i,k)=rg*K_in(i,k)/(play(i,k)-play(i,k+1))*(rho**2) + enddo +enddo +! speficic treatment for k=nlay +do i=1,ngrid + rho=pinterf(i,nlay)/rd/temp(i,nlay) + K_big(i,nlay)=rg*K_in(i,nlay)/(2*(play(i,nlay)-pinterf(i,nlay+1)))*(rho**2) +enddo + + + +! x_predict calculation for 2<=k<=nlay-1 +do k=2,nlay-1 + do i=1,ngrid + ml=(pinterf(i,k)-pinterf(i,k+1))/rg + x_predict(i,k)=x_in(i,k)-dtime/ml*(-K_big(i,k+1)*x_in(i,k+1) & + + (K_big(i,k)+K_big(i,k+1))*x_in(i,k) & + - K_big(i,k)*x_in(i,k-1)) + enddo +enddo + +! Specific treatment for k=1 +do i=1,ngrid + ml=(pinterf(i,1)-pinterf(i,2))/rg + F1=-play1(i)/rd/temp1(i)*wind1(i)*cdrag(i)*(x_in(i,1)-x_surf(i)) ! attention convention sens du flux + x_predict(i,1)=x_in(i,1)-dtime/ml*(-K_big(i,2)*(x_in(i,2) - x_in(i,1))-F1) +enddo + +! Specific treatment for k=nlay +! flux at the top of the atmosphere=0 +do i=1,ngrid + ml=0.5*(pinterf(i,nlay)-pinterf(i,nlay+1))/rg + x_predict(i,nlay)=x_in(i,nlay)+dtime/ml*(K_big(i,nlay)*(x_in(i,nlay)-x_in(i,nlay-1))) +enddo + + +end subroutine atke_explicit_prediction + + + +end module lmdz_call_atke