SUBROUTINE thermcell_closure(ngrid,nlay,r_aspect,ptimestep,rho, & & zlev,lalim,alim_star,alim_star_tot,zmax_sec,wmax_sec,zmax,wmax,f,lev_out) !------------------------------------------------------------------------- !thermcell_closure: fermeture, determination de f !------------------------------------------------------------------------- IMPLICIT NONE #include "iniprint.h" #include "thermcell.h" INTEGER ngrid,nlay INTEGER ig,k REAL r_aspect,ptimestep integer lev_out ! niveau pour les print INTEGER lalim(ngrid) REAL alim_star(ngrid,nlay) REAL alim_star_tot(ngrid) REAL rho(ngrid,nlay) REAL zlev(ngrid,nlay) REAL zmax(ngrid),zmax_sec(ngrid) REAL wmax(ngrid),wmax_sec(ngrid) real zdenom REAL alim_star2(ngrid) REAL f(ngrid) do ig=1,ngrid alim_star2(ig)=0. enddo do ig=1,ngrid if (alim_star(ig,1).LT.1.e-10) then f(ig)=0. else do k=1,lalim(ig) alim_star2(ig)=alim_star2(ig)+alim_star(ig,k)**2 & & /(rho(ig,k)*(zlev(ig,k+1)-zlev(ig,k))) enddo zdenom=max(500.,zmax(ig))*r_aspect*alim_star2(ig) if (zdenom<1.e-14) then print*,'ig=',ig print*,'alim_star2',alim_star2(ig) print*,'zmax',zmax(ig) print*,'r_aspect',r_aspect print*,'zdenom',zdenom print*,'alim_star',alim_star(ig,:) print*,'zmax_sec',zmax_sec(ig) print*,'wmax_sec',wmax_sec(ig) stop endif if ((zmax_sec(ig).gt.1.e-10).and.(iflag_thermals_ed.eq.0)) then f(ig)=wmax_sec(ig)*alim_star_tot(ig)/(max(500.,zmax_sec(ig))*r_aspect & & *alim_star2(ig)) ! f(ig)=f(ig)+(f0(ig)-f(ig))*exp((-ptimestep/ & ! & zmax_sec(ig))*wmax_sec(ig)) print*,'closure dry',f(ig),wmax_sec(ig),alim_star_tot(ig),zmax_sec(ig) else f(ig)=wmax(ig)*alim_star_tot(ig)/zdenom ! f(ig)=f(ig)+(f0(ig)-f(ig))*exp((-ptimestep/ & ! & zmax(ig))*wmax(ig)) if(prt_level.GE.10) print*,'closure moist',f(ig),wmax(ig),alim_star_tot(ig),zmax(ig) endif endif ! f0(ig)=f(ig) enddo if (prt_level.ge.1) print*,'apres fermeture' ! return end