trunk/phylmd/tlift.f90

SUBROUTINE TLIFT(P,T,RR,RS,GZ,PLCL,ICB,NK, &
     TVP,TPK,CLW,ND,NL, &
     DTVPDT1,DTVPDQ1)
  !
  ! From phylmd/tlift.F, v 1.1.1.1 2004/05/19 12:53:08

  !     Argument NK ajoute (jyg) = Niveau de depart de la
  !     convection
  !
  use YOMCST, only: rg, rcpd, RCPV, rcw, rcs, rv, rd, rlvtt, RLMLT

  implicit none

  integer, PARAMETER:: NA=60
  integer nd, icb, nk, nl
  real plcl
  REAL GZ(ND),TPK(ND),CLW(ND)
  REAL T(ND),RR(ND),RS(ND),TVP(ND),P(ND)
  REAL DTVPDT1(ND),DTVPDQ1(ND)   ! Derivatives of parcel virtual
  !                                   temperature wrt T1 and Q1
  !
  REAL QI(NA)
  REAL DTPDT1(NA),DTPDQ1(NA)      ! Derivatives of parcel temperature
  !                                   wrt T1 and Q1

  LOGICAL ICE_CONV
  real gravity, cpd, cpv, cl, ci, CPVMCL, CLMCI, EPS, alv0, alf0, CPP, cpinv
  real ah0, alv, alf, tg, s, ahg, tc, denom, es, esi, qsat_new, snew
  integer icb1, i, IMIN, j

  !--------------------------------------------------------------

  !   ***   ASSIGN VALUES OF THERMODYNAMIC CONSTANTS     ***
  ! on utilise les constantes thermo du Centre Europeen: (SB)
  !
  GRAVITY = RG !sb: Pr que gravite ne devienne pas humidite!
  !
  CPD = RCPD
  CPV = RCPV
  CL = RCW
  CI = RCS
  CPVMCL = CL-CPV
  CLMCI = CL-CI
  EPS = RD/RV
  ALV0 = RLVTT
  ALF0 = RLMLT ! (ALF0 = RLSTT-RLVTT)
  ! 
  !ccccccccccccccccccccc
  !
  !   ***  CALCULATE CERTAIN PARCEL QUANTITIES, INCLUDING STATIC ENERGY   ***
  !
  ICB1=MAX(ICB,2)
  ICB1=MIN(ICB,NL)
  !
  !jyg1
  !C      CPP=CPD*(1.-RR(1))+RR(1)*CPV
  CPP=CPD*(1.-RR(NK))+RR(NK)*CPV
  !jyg2
  CPINV=1./CPP
  !jyg1
  !         ICB may be below condensation level
  DO I=1,ICB1
     CLW(I)=0.0
  end DO
  !
  DO I=NK,ICB1
     TPK(I)=T(NK)-(GZ(I) - GZ(NK))*CPINV
     !jyg1
     !CC         TVP(I)=TPK(I)*(1.+RR(NK)/EPS)
     TVP(I)=TPK(I)*(1.+RR(NK)/EPS-RR(NK))
     !jyg2
     DTVPDT1(I) = 1.+RR(NK)/EPS-RR(NK)
     DTVPDQ1(I) = TPK(I)*(1./EPS-1.)
     !
     !jyg2

  end DO

  !
  !    ***  FIND LIFTED PARCEL TEMPERATURE AND MIXING RATIO    ***
  !
  !jyg1
  !C        AH0=(CPD*(1.-RR(1))+CL*RR(1))*T(1)
  !C     $     +RR(1)*(ALV0-CPVMCL*(T(1)-273.15))
  AH0=(CPD*(1.-RR(NK))+CL*RR(NK))*T(NK) &
       +RR(NK)*(ALV0-CPVMCL*(T(NK)-273.15)) + GZ(NK)
  !jyg2
  !
  !jyg1
  IMIN = ICB1
  !         If ICB is below LCL, start loop at ICB+1
  IF (PLCL .LT. P(ICB1)) IMIN = MIN(IMIN+1,NL)
  !
  DO  I=IMIN,NL
     !jyg2
     ALV=ALV0-CPVMCL*(T(I)-273.15)
     ALF=ALF0+CLMCI*(T(I)-273.15)

     RG=RS(I)
     TG=T(I)
     S=CPD*(1.-RR(NK))+CL*RR(NK)+ALV*ALV*RG/(RV*T(I)*T(I))
     !jyg2
     S=1./S

     DO  J=1,2
        !jyg1
        AHG=CPD*TG+(CL-CPD)*RR(NK)*TG+ALV*RG+GZ(I)
        !jyg2
        TG=TG+S*(AH0-AHG)
        TC=TG-273.15
        DENOM=243.5+TC
        DENOM=MAX(DENOM,1.0)
        !
        !       FORMULE DE BOLTON POUR PSAT
        !
        ES=6.112*EXP(17.67*TC/DENOM)
        RG=EPS*ES/(P(I)-ES*(1.-EPS))


     end DO

     !jyg1
     TPK(I)=(AH0-GZ(I)-ALV*RG)/(CPD+(CL-CPD)*RR(NK))
     !jyg2

     CLW(I)=RR(NK)-RG
     !jyg2
     CLW(I)=MAX(0.0,CLW(I))
     !jyg1
     TVP(I)=TPK(I)*(1.+RG/EPS-RR(NK))
     !jyg2
     !
     !jyg1       Derivatives
     !
     DTPDT1(I) = CPD*S
     DTPDQ1(I) = ALV*S
     !
     DTVPDT1(I) = DTPDT1(I)*(1. + RG/EPS - &
          RR(NK) + ALV*RG/(RD*TPK(I)) )
     DTVPDQ1(I) = DTPDQ1(I)*(1. + RG/EPS - &
          RR(NK) + ALV*RG/(RD*TPK(I)) ) - TPK(I)
     !
     !jyg2

  end DO
  !
  ICE_CONV = .FALSE.

  IF (ICE_CONV) THEN
     DO  I=ICB1,NL
        IF (T(I).LT.263.15) THEN
           TG=TPK(I)
           TC=TPK(I)-273.15
           DENOM=243.5+TC
           ES=6.112*EXP(17.67*TC/DENOM)
           ALV=ALV0-CPVMCL*(T(I)-273.15)
           ALF=ALF0+CLMCI*(T(I)-273.15)

           DO J=1,4
              ESI=EXP(23.33086-(6111.72784/TPK(I))+0.15215*LOG(TPK(I)))
              QSAT_NEW=EPS*ESI/(P(I)-ESI*(1.-EPS))

              SNEW= CPD*(1.-RR(NK))+CL*RR(NK) &
                   +ALV*ALV*QSAT_NEW/(RV*TPK(I)*TPK(I))
              !
              SNEW=1./SNEW
              TPK(I)=TG+(ALF*QI(I)+ALV*RG*(1.-(ESI/ES)))*SNEW
           ENDDO
           !CC        CLW(I)=RR(1)-QSAT_NEW
           CLW(I)=RR(NK)-QSAT_NEW
           CLW(I)=MAX(0.0,CLW(I))
           !jyg1
           !CC        TVP(I)=TPK(I)*(1.+QSAT_NEW/EPS)
           TVP(I)=TPK(I)*(1.+QSAT_NEW/EPS-RR(NK))
           !jyg2
        ELSE
           CONTINUE
        ENDIF

     end DO
     !
  ENDIF
  !

  !* BK :  RAJOUT DE LA TEMPERATURE DES ASCENDANCES
  !*   NON DILUES AU  NIVEAU KLEV = ND
  !*   POSONS LE ENVIRON EGAL A CELUI DE KLEV-1
  TPK(NL+1)=TPK(NL)

  RG = GRAVITY  ! RG redevient la gravite de YOMCST (sb)

END SUBROUTINE TLIFT
1	guez	3	SUBROUTINE TLIFT(P,T,RR,RS,GZ,PLCL,ICB,NK, &
2			TVP,TPK,CLW,ND,NL, &
3			DTVPDT1,DTVPDQ1)
4			!
5			! From phylmd/tlift.F, v 1.1.1.1 2004/05/19 12:53:08
6
7			! Argument NK ajoute (jyg) = Niveau de depart de la
8			! convection
9			!
10			use YOMCST, only: rg, rcpd, RCPV, rcw, rcs, rv, rd, rlvtt, RLMLT
11
12			implicit none
13
14			integer, PARAMETER:: NA=60
15			integer nd, icb, nk, nl
16			real plcl
17			REAL GZ(ND),TPK(ND),CLW(ND)
18			REAL T(ND),RR(ND),RS(ND),TVP(ND),P(ND)
19			REAL DTVPDT1(ND),DTVPDQ1(ND) ! Derivatives of parcel virtual
20			! temperature wrt T1 and Q1
21			!
22			REAL QI(NA)
23			REAL DTPDT1(NA),DTPDQ1(NA) ! Derivatives of parcel temperature
24			! wrt T1 and Q1
25
26			LOGICAL ICE_CONV
27			real gravity, cpd, cpv, cl, ci, CPVMCL, CLMCI, EPS, alv0, alf0, CPP, cpinv
28			real ah0, alv, alf, tg, s, ahg, tc, denom, es, esi, qsat_new, snew
29			integer icb1, i, IMIN, j
30
31			!--------------------------------------------------------------
32
33			! * ASSIGN VALUES OF THERMODYNAMIC CONSTANTS *
34			! on utilise les constantes thermo du Centre Europeen: (SB)
35			!
36			GRAVITY = RG !sb: Pr que gravite ne devienne pas humidite!
37			!
38			CPD = RCPD
39			CPV = RCPV
40			CL = RCW
41			CI = RCS
42			CPVMCL = CL-CPV
43			CLMCI = CL-CI
44			EPS = RD/RV
45			ALV0 = RLVTT
46			ALF0 = RLMLT ! (ALF0 = RLSTT-RLVTT)
47			!
48			!ccccccccccccccccccccc
49			!
50			! * CALCULATE CERTAIN PARCEL QUANTITIES, INCLUDING STATIC ENERGY *
51			!
52			ICB1=MAX(ICB,2)
53			ICB1=MIN(ICB,NL)
54			!
55			!jyg1
56			!C CPP=CPD(1.-RR(1))+RR(1)CPV
57			CPP=CPD(1.-RR(NK))+RR(NK)CPV
58			!jyg2
59			CPINV=1./CPP
60			!jyg1
61			! ICB may be below condensation level
62			DO I=1,ICB1
63			CLW(I)=0.0
64			end DO
65			!
66			DO I=NK,ICB1
67			TPK(I)=T(NK)-(GZ(I) - GZ(NK))*CPINV
68			!jyg1
69			!CC TVP(I)=TPK(I)*(1.+RR(NK)/EPS)
70			TVP(I)=TPK(I)*(1.+RR(NK)/EPS-RR(NK))
71			!jyg2
72			DTVPDT1(I) = 1.+RR(NK)/EPS-RR(NK)
73			DTVPDQ1(I) = TPK(I)*(1./EPS-1.)
74			!
75			!jyg2
76
77			end DO
78
79			!
80			! * FIND LIFTED PARCEL TEMPERATURE AND MIXING RATIO *
81			!
82			!jyg1
83			!C AH0=(CPD(1.-RR(1))+CLRR(1))*T(1)
84			!C $ +RR(1)(ALV0-CPVMCL(T(1)-273.15))
85			AH0=(CPD(1.-RR(NK))+CLRR(NK))*T(NK) &
86			+RR(NK)(ALV0-CPVMCL(T(NK)-273.15)) + GZ(NK)
87			!jyg2
88			!
89			!jyg1
90			IMIN = ICB1
91			! If ICB is below LCL, start loop at ICB+1
92			IF (PLCL .LT. P(ICB1)) IMIN = MIN(IMIN+1,NL)
93			!
94			DO I=IMIN,NL
95			!jyg2
96			ALV=ALV0-CPVMCL*(T(I)-273.15)
97			ALF=ALF0+CLMCI*(T(I)-273.15)
98
99			RG=RS(I)
100			TG=T(I)
101			S=CPD(1.-RR(NK))+CLRR(NK)+ALVALVRG/(RVT(I)T(I))
102			!jyg2
103			S=1./S
104
105			DO J=1,2
106			!jyg1
107			AHG=CPDTG+(CL-CPD)RR(NK)TG+ALVRG+GZ(I)
108			!jyg2
109			TG=TG+S*(AH0-AHG)
110			TC=TG-273.15
111			DENOM=243.5+TC
112			DENOM=MAX(DENOM,1.0)
113			!
114			! FORMULE DE BOLTON POUR PSAT
115			!
116			ES=6.112EXP(17.67TC/DENOM)
117			RG=EPSES/(P(I)-ES(1.-EPS))
118
119
120			end DO
121
122			!jyg1
123			TPK(I)=(AH0-GZ(I)-ALVRG)/(CPD+(CL-CPD)RR(NK))
124			!jyg2
125
126			CLW(I)=RR(NK)-RG
127			!jyg2
128			CLW(I)=MAX(0.0,CLW(I))
129			!jyg1
130			TVP(I)=TPK(I)*(1.+RG/EPS-RR(NK))
131			!jyg2
132			!
133			!jyg1 Derivatives
134			!
135			DTPDT1(I) = CPD*S
136			DTPDQ1(I) = ALV*S
137			!
138			DTVPDT1(I) = DTPDT1(I)*(1. + RG/EPS - &
139			RR(NK) + ALVRG/(RDTPK(I)) )
140			DTVPDQ1(I) = DTPDQ1(I)*(1. + RG/EPS - &
141			RR(NK) + ALVRG/(RDTPK(I)) ) - TPK(I)
142			!
143			!jyg2
144
145			end DO
146			!
147			ICE_CONV = .FALSE.
148
149			IF (ICE_CONV) THEN
150			DO I=ICB1,NL
151			IF (T(I).LT.263.15) THEN
152			TG=TPK(I)
153			TC=TPK(I)-273.15
154			DENOM=243.5+TC
155			ES=6.112EXP(17.67TC/DENOM)
156			ALV=ALV0-CPVMCL*(T(I)-273.15)
157			ALF=ALF0+CLMCI*(T(I)-273.15)
158
159			DO J=1,4
160			ESI=EXP(23.33086-(6111.72784/TPK(I))+0.15215*LOG(TPK(I)))
161			QSAT_NEW=EPSESI/(P(I)-ESI(1.-EPS))
162
163			SNEW= CPD(1.-RR(NK))+CLRR(NK) &
164			+ALVALVQSAT_NEW/(RVTPK(I)TPK(I))
165			!
166			SNEW=1./SNEW
167			TPK(I)=TG+(ALFQI(I)+ALVRG(1.-(ESI/ES)))SNEW
168			ENDDO
169			!CC CLW(I)=RR(1)-QSAT_NEW
170			CLW(I)=RR(NK)-QSAT_NEW
171			CLW(I)=MAX(0.0,CLW(I))
172			!jyg1
173			!CC TVP(I)=TPK(I)*(1.+QSAT_NEW/EPS)
174			TVP(I)=TPK(I)*(1.+QSAT_NEW/EPS-RR(NK))
175			!jyg2
176			ELSE
177			CONTINUE
178			ENDIF
179
180			end DO
181			!
182			ENDIF
183			!
184
185			!* BK : RAJOUT DE LA TEMPERATURE DES ASCENDANCES
186			!* NON DILUES AU NIVEAU KLEV = ND
187			!* POSONS LE ENVIRON EGAL A CELUI DE KLEV-1
188			TPK(NL+1)=TPK(NL)
189
190			RG = GRAVITY ! RG redevient la gravite de YOMCST (sb)
191
192			END SUBROUTINE TLIFT