thermcell_qsat.F90

Go to the documentation of this file.

subroutine thermcell_qsat(klon,active,pplev,ztemp,zqta,zqsat)
implicit none

#include "YOMCST.h"
#include "YOETHF.h"
#include "FCTTRE.h"


!====================================================================
! DECLARATIONS
!====================================================================

! Arguments
INTEGER klon
REAL zpspsk(klon),pplev(klon)
REAL ztemp(klon),zqta(klon),zqsat(klon)
LOGICAL active(klon)

! Variables locales
INTEGER ig,iter
REAL tbef(klon),dt(klon)
REAL tdelta,qsatbef,zcor,qlbef,zdelta,zcvm5,dqsat,num,denom,dqsat_dt
logical zsat
REAL rlvcp
REAL, SAVE :: ddt0=.01
LOGICAL afaire(klon),tout_converge

!====================================================================
! INITIALISATIONS
!====================================================================

rlvcp = rlvtt/rcpd
tout_converge=.false.
afaire(:)=.false.
dt(:)=0.


!====================================================================
! Routine a vectoriser en copiant active dans converge et en mettant
! la boucle sur les iterations a l'exterieur est en mettant
! converge= false des que la convergence est atteinte.
!====================================================================

do ig=1,klon
   if (active(ig)) then
               tbef(ig)=ztemp(ig)
               zdelta=max(0.,sign(1.,rtt-tbef(ig)))
               qsatbef= r2es * foeew(tbef(ig),zdelta)/pplev(ig)
               qsatbef=min(0.5,qsatbef)
               zcor=1./(1.-retv*qsatbef)
               qsatbef=qsatbef*zcor
               qlbef=max(0.,zqta(ig)-qsatbef)
               dt(ig) = 0.5*rlvcp*qlbef
               zqsat(ig)=qsatbef
     endif
enddo

! Traitement du cas ou il y a condensation mais faible
! On ne condense pas mais on dit que le qsat est le qta
do ig=1,klon
   if (active(ig)) then
      if (0.<abs(dt(ig)).and.abs(dt(ig))<=ddt0) then
         zqsat(ig)=zqta(ig)
      endif
   endif
enddo

do iter=1,10
    afaire(:)=abs(dt(:)).gt.ddt0
    do ig=1,klon
               if (afaire(ig)) then
                 tbef(ig)=tbef(ig)+dt(ig)
                 zdelta=max(0.,sign(1.,rtt-tbef(ig)))
                 qsatbef= r2es * foeew(tbef(ig),zdelta)/pplev(ig)
                 qsatbef=min(0.5,qsatbef)
                 zcor=1./(1.-retv*qsatbef)
                 qsatbef=qsatbef*zcor
                 qlbef=zqta(ig)-qsatbef
                 zdelta=max(0.,sign(1.,rtt-tbef(ig)))
                 zcvm5=r5les*(1.-zdelta) + r5ies*zdelta
                 zcor=1./(1.-retv*qsatbef)
                 dqsat_dt=foede(tbef(ig),zdelta,zcvm5,qsatbef,zcor)
                 num=-tbef(ig)+ztemp(ig)+rlvcp*qlbef
                 denom=1.+rlvcp*dqsat_dt
                 zqsat(ig) = qsatbef
                 dt(ig)=num/denom
               endif
    enddo
enddo

return
end