source: LMDZ.3.3/branches/rel-LF/libf/phylmd/tlift.F @ 1546

Last change on this file since 1546 was 230, checked in by lmdzadmin, 23 years ago

Merge de la physique avec la branche principale
LF

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 5.8 KB
Line 
1c $Header$
2        SUBROUTINE TLIFT(P,T,RR,RS,GZ,PLCL,ICB,NK,
3     .                  TVP,TPK,CLW,ND,NL,
4     .                  DTVPDT1,DTVPDQ1)
5C
6C     Argument NK ajoute (jyg) = Niveau de depart de la
7C     convection
8C
9        PARAMETER (NA=60)
10        REAL GZ(ND),TPK(ND),CLW(ND)
11        REAL T(ND),RR(ND),RS(ND),TVP(ND),P(ND)
12        REAL DTVPDT1(ND),DTVPDQ1(ND)   ! Derivatives of parcel virtual
13C                                   temperature wrt T1 and Q1
14C
15        REAL CLW_NEW(NA),QI(NA)
16        REAL DTPDT1(NA),DTPDQ1(NA)      ! Derivatives of parcel temperature
17C                                   wrt T1 and Q1
18 
19C
20        LOGICAL ICE_CONV
21C
22C   ***   ASSIGN VALUES OF THERMODYNAMIC CONSTANTS     ***
23C
24c sb        CPD=1005.7
25c sb      CPV=1870.0
26c sb        CL=4190.0
27c sb        CPVMCL=2320.0
28c sb        RV=461.5
29c sb        RD=287.04
30c sb        EPS=RD/RV
31c sb        ALV0=2.501E6
32ccccccccccccccccccccccc
33c constantes coherentes avec le modele du Centre Europeen
34c sb      RD = 1000.0 * 1.380658E-23 * 6.0221367E+23 / 28.9644
35c sb      RV = 1000.0 * 1.380658E-23 * 6.0221367E+23 / 18.0153
36c sb      CPD = 3.5 * RD
37c sb      CPV = 4.0 * RV
38c sb      CL = 4218.0
39c sb      CI=2090.0
40c sb      CPVMCL=CL-CPV
41c sb      CLMCI=CL-CI
42c sb      EPS=RD/RV
43c sb      ALV0=2.5008E+06
44c sb      ALF0=3.34E+05
45 
46cccccccccccc
47c on utilise les constantes thermo du Centre Europeen: (SB)
48c
49#include "YOMCST.h"
50       GRAVITY = RG !sb: Pr que gravite ne devienne pas humidite!
51c
52       CPD = RCPD
53       CPV = RCPV
54       CL = RCW
55       CI = RCS
56       CPVMCL = CL-CPV
57       CLMCI = CL-CI
58       EPS = RD/RV
59       ALV0 = RLVTT
60       ALF0 = RLMLT ! (ALF0 = RLSTT-RLVTT)
61c
62cccccccccccccccccccccc
63C
64C   ***  CALCULATE CERTAIN PARCEL QUANTITIES, INCLUDING STATIC ENERGY   ***
65C
66        ICB1=MAX(ICB,2)
67        ICB1=MIN(ICB,NL)
68C
69Cjyg1
70CC      CPP=CPD*(1.-RR(1))+RR(1)*CPV
71      CPP=CPD*(1.-RR(NK))+RR(NK)*CPV
72Cjyg2
73      CPINV=1./CPP
74Cjyg1
75C         ICB may be below condensation level
76CCC        DO 100 I=1,ICB1-1
77CCC         TPK(I)=T(1)-GZ(I)*CPINV
78CCC         TVP(I)=TPK(I)*(1.+RR(1)/EPS)
79        DO 50 I=1,ICB1
80         CLW(I)=0.0
8150      CONTINUE
82C
83        DO 100 I=NK,ICB1
84         TPK(I)=T(NK)-(GZ(I) - GZ(NK))*CPINV
85Cjyg1
86CCC         TVP(I)=TPK(I)*(1.+RR(NK)/EPS)
87         TVP(I)=TPK(I)*(1.+RR(NK)/EPS-RR(NK))
88Cjyg2
89         DTVPDT1(I) = 1.+RR(NK)/EPS-RR(NK)
90         DTVPDQ1(I) = TPK(I)*(1./EPS-1.)
91C
92Cjyg2
93 
94  100   CONTINUE
95 
96C
97C    ***  FIND LIFTED PARCEL TEMPERATURE AND MIXING RATIO    ***
98C
99Cjyg1
100CC        AH0=(CPD*(1.-RR(1))+CL*RR(1))*T(1)
101CC     $     +RR(1)*(ALV0-CPVMCL*(T(1)-273.15))
102        AH0=(CPD*(1.-RR(NK))+CL*RR(NK))*T(NK)
103     $     +RR(NK)*(ALV0-CPVMCL*(T(NK)-273.15)) + GZ(NK)
104Cjyg2
105C
106Cjyg1
107        IMIN = ICB1
108C         If ICB is below LCL, start loop at ICB+1
109        IF (PLCL .LT. P(ICB1)) IMIN = MIN(IMIN+1,NL)
110C
111CCC        DO 300 I=ICB1,NL
112        DO 300 I=IMIN,NL
113Cjyg2
114         ALV=ALV0-CPVMCL*(T(I)-273.15)
115         ALF=ALF0+CLMCI*(T(I)-273.15)
116 
117        RG=RS(I)
118        TG=T(I)
119C       S=CPD+ALV*ALV*RG/(RV*T(I)*T(I))
120Cjyg1
121CC        S=CPD*(1.-RR(1))+CL*RR(1)+ALV*ALV*RG/(RV*T(I)*T(I))
122        S=CPD*(1.-RR(NK))+CL*RR(NK)+ALV*ALV*RG/(RV*T(I)*T(I))
123Cjyg2
124        S=1./S
125 
126        DO 200 J=1,2
127Cjyg1
128CC         AHG=CPD*TG+(CL-CPD)*RR(1)*TG+ALV*RG+GZ(I)
129         AHG=CPD*TG+(CL-CPD)*RR(NK)*TG+ALV*RG+GZ(I)
130Cjyg2
131        TG=TG+S*(AH0-AHG)
132        TC=TG-273.15
133        DENOM=243.5+TC
134        DENOM=MAX(DENOM,1.0)
135C
136C       FORMULE DE BOLTON POUR PSAT
137C
138        ES=6.112*EXP(17.67*TC/DENOM)
139        RG=EPS*ES/(P(I)-ES*(1.-EPS))
140 
141 
142  200   CONTINUE
143 
144Cjyg1
145CC        TPK(I)=(AH0-GZ(I)-ALV*RG)/(CPD+(CL-CPD)*RR(1))
146        TPK(I)=(AH0-GZ(I)-ALV*RG)/(CPD+(CL-CPD)*RR(NK))
147Cjyg2
148C       TPK(I)=(AH0-GZ(I)-ALV*RG-(CL-CPD)*T(I)*RR(1))/CPD
149 
150Cjyg1
151CC        CLW(I)=RR(1)-RG
152        CLW(I)=RR(NK)-RG
153Cjyg2
154        CLW(I)=MAX(0.0,CLW(I))
155Cjyg1
156CCC        TVP(I)=TPK(I)*(1.+RG/EPS)
157        TVP(I)=TPK(I)*(1.+RG/EPS-RR(NK))
158Cjyg2
159C
160Cjyg1       Derivatives
161C
162        DTPDT1(I) = CPD*S
163        DTPDQ1(I) = ALV*S
164C
165         DTVPDT1(I) = DTPDT1(I)*(1. + RG/EPS -
166     .           RR(NK) + ALV*RG/(RD*TPK(I)) )
167        DTVPDQ1(I) = DTPDQ1(I)*(1. + RG/EPS -
168     .           RR(NK) + ALV*RG/(RD*TPK(I)) ) - TPK(I)
169C
170Cjyg2
171 
172  300   CONTINUE
173C
174      ICE_CONV = .FALSE.
175
176      IF (ICE_CONV) THEN
177C
178CJAM
179C       RAJOUT DE LA PROCEDURE ICEFRAC
180C
181c sb        CALL ICEFRAC(T,CLW,CLW_NEW,QI,ND,NL)
182 
183        DO 400 I=ICB1,NL
184        IF (T(I).LT.263.15) THEN
185        TG=TPK(I)
186        TC=TPK(I)-273.15
187        DENOM=243.5+TC
188        ES=6.112*EXP(17.67*TC/DENOM)
189        ALV=ALV0-CPVMCL*(T(I)-273.15)
190        ALF=ALF0+CLMCI*(T(I)-273.15)
191 
192        DO J=1,4
193        ESI=EXP(23.33086-(6111.72784/TPK(I))+0.15215*LOG(TPK(I)))
194        QSAT_NEW=EPS*ESI/(P(I)-ESI*(1.-EPS))
195CCC        SNEW= CPD*(1.-RR(1))+CL*RR(1)+ALV*ALV*QSAT_NEW/(RV*TPK(I)*TPK(I))
196        SNEW= CPD*(1.-RR(NK))+CL*RR(NK)
197     .        +ALV*ALV*QSAT_NEW/(RV*TPK(I)*TPK(I))
198C
199        SNEW=1./SNEW
200        TPK(I)=TG+(ALF*QI(I)+ALV*RG*(1.-(ESI/ES)))*SNEW
201c$$$        PRINT*,'################################'
202c$$$        PRINT*,TPK(I)
203c$$$        PRINT*,(ALF*QI(I)+ALV*RG*(1.-(ESI/ES)))*SNEW
204        ENDDO
205CCC        CLW(I)=RR(1)-QSAT_NEW
206        CLW(I)=RR(NK)-QSAT_NEW
207        CLW(I)=MAX(0.0,CLW(I))
208Cjyg1
209CCC        TVP(I)=TPK(I)*(1.+QSAT_NEW/EPS)
210        TVP(I)=TPK(I)*(1.+QSAT_NEW/EPS-RR(NK))
211Cjyg2
212        ELSE
213        CONTINUE
214        ENDIF
215 
216  400   CONTINUE
217C
218      ENDIF
219C
220 
221******************************************************
222** BK :  RAJOUT DE LA TEMPERATURE DES ASCENDANCES
223**   NON DILUES AU  NIVEAU KLEV = ND
224**   POSONS LE ENVIRON EGAL A CELUI DE KLEV-1
225********************************************************
226 
227      TPK(NL+1)=TPK(NL)
228 
229*******************************************************
230
231      RG = GRAVITY  ! RG redevient la gravite de YOMCST (sb)
232 
233 
234        RETURN
235        END
236 
237 
238 
239 
240 
241 
242 
243 
Note: See TracBrowser for help on using the repository browser.