libf/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: 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 with regard to T1 and Q1

  REAL QI(NA)

  REAL DTPDT1(NA),DTPDQ1(NA)
  ! Derivatives of parcel temperature with regard to 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)

  CPD = RCPD
  CPV = RCPV
  CL = RCW
  CI = RCS
  CPVMCL = CL-CPV
  CLMCI = CL-CI
  EPS = RD/RV
  ALV0 = RLVTT
  ALF0 = RLMLT ! (ALF0 = RLSTT-RLVTT)

  ! ***  CALCULATE CERTAIN PARCEL QUANTITIES, INCLUDING STATIC ENERGY   ***

  ICB1=MIN(ICB,NL)

  CPP=CPD*(1.-RR(NK))+RR(NK)*CPV
  CPINV=1./CPP

  !         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
     TVP(I)=TPK(I)*(1.+RR(NK)/EPS-RR(NK))
     DTVPDT1(I) = 1.+RR(NK)/EPS-RR(NK)
     DTVPDQ1(I) = TPK(I)*(1./EPS-1.)
  end DO

  !    ***  FIND LIFTED PARCEL TEMPERATURE AND MIXING RATIO    ***

  AH0=(CPD*(1.-RR(NK))+CL*RR(NK))*T(NK) &
       +RR(NK)*(ALV0-CPVMCL*(T(NK)-273.15)) + GZ(NK)
  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
     ALV=ALV0-CPVMCL*(T(I)-273.15)
     ALF=ALF0+CLMCI*(T(I)-273.15)

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

     S=1./S

     DO  J=1,2
        AHG=CPD*TG+(CL-CPD)*RR(NK)*TG+ALV*GRAVITY+GZ(I)

        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)
        GRAVITY=EPS*ES/(P(I)-ES*(1.-EPS))
     end DO

     TPK(I)=(AH0-GZ(I)-ALV*GRAVITY)/(CPD+(CL-CPD)*RR(NK))
     CLW(I)=RR(NK)-GRAVITY
     CLW(I)=MAX(0.0,CLW(I))
     TVP(I)=TPK(I)*(1.+GRAVITY/EPS-RR(NK))

     !jyg1       Derivatives
     DTPDT1(I) = CPD*S
     DTPDQ1(I) = ALV*S

     DTVPDT1(I) = DTPDT1(I)*(1. + GRAVITY/EPS - &
          RR(NK) + ALV*GRAVITY/(RD*TPK(I)) )
     DTVPDQ1(I) = DTPDQ1(I)*(1. + GRAVITY/EPS - &
          RR(NK) + ALV*GRAVITY/(RD*TPK(I)) ) - TPK(I)
  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*GRAVITY*(1.-(ESI/ES)))*SNEW
           ENDDO
           CLW(I)=RR(NK)-QSAT_NEW
           CLW(I)=MAX(0.0,CLW(I))
           TVP(I)=TPK(I)*(1.+QSAT_NEW/EPS-RR(NK))
        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)

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