SUBROUTINE thermcell_closure(ngrid,nlay,r_aspect,ptimestep,rho, & & zlev,lalim,alim_star,zmax_sec,wmax_sec,zmax,wmax,f,f0,lev_out) !------------------------------------------------------------------------- !thermcell_closure: fermeture, determination de f !------------------------------------------------------------------------- IMPLICIT NONE 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 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) REAL f0(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.(1.eq.1)) then f(ig)=wmax_sec(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)) else f(ig)=wmax(ig)/zdenom f(ig)=f(ig)+(f0(ig)-f(ig))*exp((-ptimestep/ & & zmax(ig))*wmax(ig)) endif endif f0(ig)=f(ig) enddo if (lev_out.ge.1) print*,'apres fermeture' ! return end