Index: LMDZ6/trunk/libf/phylmd/atke_turbulence_ini_mod.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/atke_turbulence_ini_mod.F90 (revision 4543) +++ LMDZ6/trunk/libf/phylmd/atke_turbulence_ini_mod.F90 (revision 4545) @@ -5,6 +5,6 @@ save - integer :: iflag_atke - !$OMP THREADPRIVATE(iflag_atke) +integer :: iflag_atke, iflag_num_atke +!$OMP THREADPRIVATE(iflag_atke, iflag_num_atke) real :: kappa = 0.4 ! Von Karman constant !$OMP THREADPRIVATE(kappa) @@ -45,4 +45,8 @@ CALL getin_p('iflag_atke',iflag_atke) + ! 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 Index: LMDZ6/trunk/libf/phylmd/call_atke_mod.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/call_atke_mod.F90 (revision 4545) +++ LMDZ6/trunk/libf/phylmd/call_atke_mod.F90 (revision 4545) @@ -0,0 +1,160 @@ +module call_atke_mod + +USE atke_exchange_coeff_mod, 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,play,pinterf, & + tke,Km_out,Kh_out) + + + + +USE atke_turbulence_ini_mod, 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) :: 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 + + +REAL, DIMENSION(ngrid,nlay) :: wind_u_predict, wind_v_predict +REAL, DIMENSION(ngrid) :: wind1 +INTEGER i + + + +call atke_compute_km_kh(ngrid,nlay, & + wind_u,wind_v,temp,play,pinterf, & + tke,Km_out,Kh_out) + +if (iflag_num_atke .EQ. 1) then + + !! 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, & + wind_u_predict,wind_v_predict,temp,play,pinterf, & + 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 + +!K_big(:,1)=0. +!do k=1,nlay +! do i=1,ngrid +! print*, 'youhou', k, x_in(i,k), x_predict(i,k), K_big(i,k) +! end do +!enddo + + + +end subroutine atke_explicit_prediction + + + +end module call_atke_mod Index: LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90 (revision 4543) +++ LMDZ6/trunk/libf/phylmd/pbl_surface_mod.F90 (revision 4545) @@ -24,5 +24,5 @@ USE climb_wind_mod, ONLY : climb_wind_down, climb_wind_up USE coef_diff_turb_mod, ONLY : coef_diff_turb - USE atke_exchange_coeff_mod, ONLY : atke_compute_km_kh + USE call_atke_mod, ONLY : call_atke USE ioipsl_getin_p_mod, ONLY : getin_p USE cdrag_mod @@ -790,4 +790,5 @@ REAL, DIMENSION(klon) :: uzon_w, vmer_w, speed_w, zri1_w, pref_w !speed_w, zri1_w, pref_w, added by Fuxing WANG, 04/03/2015 REAL, DIMENSION(klon) :: zgeo1_w, tair1_w, qair1_w, tairsol_w + REAL, DIMENSION(klon) :: yus0, yvs0 !!! jyg le 25/03/2013 @@ -993,4 +994,5 @@ cdragh(:)=0. ; cdragm(:)=0. zu1(:)=0. ; zv1(:)=0. + yus0(:)=0. ; yvs0(:)=0. !albedo SB >>> alb_dir_m=0. ; alb_dif_m=0. ; alb3_lic(:)=0. @@ -1699,5 +1701,5 @@ IF (iflag_pbl>=50) THEN - CALL atke_compute_km_kh(knon,klev,yu,yv,yt, & + CALL call_atke(dtime,knon,klev,ycdragm, ycdragh,yus0,yvs0,yts,yu,yv,yt, & ypplay,ypaprs,ytke,ycoefm, ycoefh) @@ -1743,5 +1745,5 @@ IF (iflag_pbl>=50) THEN - CALL atke_compute_km_kh(knon,klev,yu_x,yv_x,yt_x, & + CALL call_atke(dtime,knon,klev,ycdragm_x,ycdragh_x,yus0,yvs0,yts_x,yu_x,yv_x,yt_x, & ypplay,ypaprs,ytke_x,ycoefm_x, ycoefh_x) @@ -1782,5 +1784,5 @@ IF (iflag_pbl>=50) THEN - CALL atke_compute_km_kh(knon,klev,yu_w,yv_w,yt_w, & + CALL call_atke(dtime,knon,klev,ycdragm_w,ycdragh_w,yus0,yvs0,yts_w,yu_w,yv_w,yt_w, & ypplay,ypaprs,ytke_w,ycoefm_w, ycoefh_w) Index: LMDZ6/trunk/libf/phylmdiso/call_atke_mod.F90 =================================================================== --- LMDZ6/trunk/libf/phylmdiso/call_atke_mod.F90 (revision 4545) +++ LMDZ6/trunk/libf/phylmdiso/call_atke_mod.F90 (revision 4545) @@ -0,0 +1,1 @@ +link ../phylmd/call_atke_mod.F90