source: LMDZ4/branches/LMDZ4-dev-20091210/libf/phylmd/tlift.F @ 3655

Last change on this file since 3655 was 776, checked in by Laurent Fairhead, 17 years ago

Suite du merge entre la version et la HEAD: quelques modifications
de Yann sur le

LF

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