SUBROUTINE thermcell_env(ngrid,nlay,po,pt,pu,pv,pplay, & & pplev,zo,zh,zl,ztv,zthl,zu,zv,zpspsk,zqsat,lev_out) !-------------------------------------------------------------- !thermcell_env: calcule les caracteristiques de l environnement !necessaires au calcul des proprietes dans le thermique !-------------------------------------------------------------- IMPLICIT NONE #include "YOMCST.h" #include "YOETHF.h" #include "FCTTRE.h" #include "iniprint.h" INTEGER ngrid,nlay REAL po(ngrid,nlay) REAL pt(ngrid,nlay) REAL pu(ngrid,nlay) REAL pv(ngrid,nlay) REAL pplay(ngrid,nlay) REAL pplev(ngrid,nlay+1) integer lev_out ! niveau pour les print REAL zo(ngrid,nlay) REAL zl(ngrid,nlay) REAL zh(ngrid,nlay) REAL ztv(ngrid,nlay) REAL zthl(ngrid,nlay) REAL zpspsk(ngrid,nlay) REAL zu(ngrid,nlay) REAL zv(ngrid,nlay) REAL zqsat(ngrid,nlay) INTEGER ig,l,ll real zcor,zdelta,zcvm5,qlbef real Tbef,qsatbef real dqsat_dT,DT,num,denom REAL RLvCp,DDT0 PARAMETER (DDT0=.01) LOGICAL Zsat Zsat=.false. ! ! Pr Tprec=Tl calcul de qsat ! Si qsat>qT T=Tl, q=qT ! Sinon DDT=(-Tprec+Tl+RLVCP (qT-qsat(T')) / (1+RLVCP dqsat/dt) ! On cherche DDT < DDT0 ! ! calcul des caracteristiques de l environnement DO ll=1,nlay DO ig=1,ngrid zo(ig,ll)=po(ig,ll) zl(ig,ll)=0. zh(ig,ll)=pt(ig,ll) zqsat(ig,ll)=0. EndDO EndDO ! ! !recherche de saturation dans l environnement DO ll=1,nlay ! les points insatures sont definitifs DO ig=1,ngrid Tbef=pt(ig,ll) zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) qsatbef= R2ES * FOEEW(Tbef,zdelta)/pplev(ig,ll) qsatbef=MIN(0.5,qsatbef) zcor=1./(1.-retv*qsatbef) qsatbef=qsatbef*zcor Zsat = (max(0.,po(ig,ll)-qsatbef) .gt. 1.e-10) if (Zsat) then qlbef=max(0.,po(ig,ll)-qsatbef) ! si sature: ql est surestime, d'ou la sous-relax DT = 0.5*RLvCp*qlbef ! on pourra enchainer 2 ou 3 calculs sans Do while do while (abs(DT).gt.DDT0) ! il faut verifier si c,a conserve quand on repasse en insature ... Tbef=Tbef+DT zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) qsatbef= R2ES * FOEEW(Tbef,zdelta)/pplev(ig,ll) qsatbef=MIN(0.5,qsatbef) zcor=1./(1.-retv*qsatbef) qsatbef=qsatbef*zcor ! on veut le signe de qlbef qlbef=po(ig,ll)-qsatbef ! dqsat_dT zdelta=MAX(0.,SIGN(1.,RTT-Tbef)) zcvm5=R5LES*(1.-zdelta) + R5IES*zdelta zcor=1./(1.-retv*qsatbef) dqsat_dT=FOEDE(Tbef,zdelta,zcvm5,qsatbef,zcor) num=-Tbef+pt(ig,ll)+RLvCp*qlbef denom=1.+RLvCp*dqsat_dT if (denom.lt.1.e-10) then print*,'pb denom' endif DT=num/denom enddo ! on ecrit de maniere conservative (sat ou non) zl(ig,ll) = max(0.,qlbef) ! T = Tl +Lv/Cp ql zh(ig,ll) = pt(ig,ll)+RLvCp*zl(ig,ll) zo(ig,ll) = po(ig,ll)-zl(ig,ll) endif !on ecrit zqsat zqsat(ig,ll)=qsatbef EndDO EndDO ! ! !----------------------------------------------------------------------- ! incrementation eventuelle de tendances precedentes: ! --------------------------------------------------- if (prt_level.ge.1) print*,'0 OK convect8' DO 1010 l=1,nlay DO 1015 ig=1,ngrid zpspsk(ig,l)=(pplay(ig,l)/100000.)**RKAPPA zu(ig,l)=pu(ig,l) zv(ig,l)=pv(ig,l) !attention zh est maintenant le profil de T et plus le profil de theta ! ! ! T-> Theta ztv(ig,l)=zh(ig,l)/zpspsk(ig,l) !Theta_v ztv(ig,l)=ztv(ig,l)*(1.+RETV*(zo(ig,l)) & & -zl(ig,l)) !Thetal zthl(ig,l)=pt(ig,l)/zpspsk(ig,l) ! 1015 CONTINUE 1010 CONTINUE RETURN END