tilft43.F

Go to the documentation of this file.

!
! $Header$
!
        SUBROUTINE tlift43(P,T,Q,QS,GZ,ICB,NK,TVP,TPK,CLW,ND,NL,KK)
        REAL gz(nd),tpk(nd),clw(nd),p(nd)
        REAL t(nd),q(nd),qs(nd),tvp(nd),lv0
c
c   ***   assign values of thermodynamic constants     ***
c
c -- sb:
c!      CPD=1005.7
c!      CPV=1870.0
c!      CL=4190.0
c!      RV=461.5
c!      RD=287.04
c!      LV0=2.501E6
c!      G=9.8
c!      ROWL=1000.0
c ajouts:
#include "YOMCST.h"
        cpd = rcpd
        cpv = rcpv
        cl = rcw
        lv0 = rlvtt
        g = rg
        rowl= ratm/100.
        gravity = rg !sb: Pr que gravite ne devienne pas humidite!
c sb --
c
        cpvmcl=cl-cpv
        eps=rd/rv
        epsi=1./eps
c
c   ***  calculate certain parcel quantities, including static energy   ***
c
        ah0=(cpd*(1.-q(nk))+cl*q(nk))*t(nk)+q(nk)*(lv0-cpvmcl*(
     1   t(nk)-273.15))+gz(nk)
        cpp=cpd*(1.-q(nk))+q(nk)*cpv
        cpinv=1./cpp
c
        IF(kk.EQ.1)THEN
c
c   ***   calculate lifted parcel quantities below cloud base   ***
c
        DO 50 i=1,icb-1
         clw(i)=0.0
   50   CONTINUE
        DO 100 i=nk,icb-1
         tpk(i)=t(nk)-(gz(i)-gz(nk))*cpinv
         tvp(i)=tpk(i)*(1.+q(nk)*epsi)
  100   CONTINUE
        END IF
c
c    ***  find lifted parcel quantities above cloud base    ***
c
        nst=icb
        nsb=icb
        IF(kk.EQ.2)THEN  
         nst=nl
         nsb=icb+1
        END IF
        DO 300 i=nsb,nst
         tg=t(i)
         qg=qs(i)
         alv=lv0-cpvmcl*(t(i)-273.15)
         DO 200 j=1,2
          s=cpd+alv*alv*qg/(rv*t(i)*t(i))
          s=1./s
          ahg=cpd*tg+(cl-cpd)*q(nk)*t(i)+alv*qg+gz(i)
          tg=tg+s*(ah0-ahg)
          tg=max(tg,35.0)
          tc=tg-273.15
          denom=243.5+tc
          IF(tc.GE.0.0)THEN  
           es=6.112*exp(17.67*tc/denom)
          ELSE  
           es=exp(23.33086-6111.72784/tg+0.15215*log(tg))
          END IF  
          qg=eps*es/(p(i)-es*(1.-eps))
  200    CONTINUE
         alv=lv0-cpvmcl*(t(i)-273.15)
         tpk(i)=(ah0-(cl-cpd)*q(nk)*t(i)-gz(i)-alv*qg)/cpd
         clw(i)=q(nk)-qg
         clw(i)=max(0.0,clw(i))
         rg=qg/(1.-q(nk))
         tvp(i)=tpk(i)*(1.+rg*epsi)
  300   CONTINUE

c -- sb:
        rg = gravity  ! RG redevient la gravite de YOMCST (sb)
c sb --

        RETURN
        END