GCC Code Coverage Report


Directory: ./
File: phys/tlift.f90
Date: 2022-01-11 19:19:34
Exec Total Coverage
Lines: 0 45 0.0%
Branches: 0 10 0.0%

Line Branch Exec Source
1
2 ! $Header$
3
4 SUBROUTINE tlift(p, t, rr, rs, gz, plcl, icb, nk, tvp, tpk, clw, nd, nl, &
5 dtvpdt1, dtvpdq1)
6 IMPLICIT NONE
7 ! Argument NK ajoute (jyg) = Niveau de depart de la
8 ! convection
9 INTEGER icb, nk, nd, nl
10 INTEGER,PARAMETER :: na=60
11 REAL gz(nd), tpk(nd), clw(nd), plcl
12 REAL t(nd), rr(nd), rs(nd), tvp(nd), p(nd)
13 REAL dtvpdt1(nd), dtvpdq1(nd) ! Derivatives of parcel virtual
14 ! temperature wrt T1 and Q1
15
16 REAL clw_new(na), qi(na)
17 REAL dtpdt1(na), dtpdq1(na) ! Derivatives of parcel temperature
18 ! wrt T1 and Q1
19 REAL gravity, cpd, cpv, cl, ci, cpvmcl, clmci, eps, alv0, alf0
20 REAL cpp, cpinv, ah0, alf, tg, s, ahg, tc, denom, alv, es, esi
21 REAL qsat_new, snew
22 INTEGER icbl, i, imin, j, icb1
23
24 LOGICAL ice_conv
25
26 ! *** ASSIGN VALUES OF THERMODYNAMIC CONSTANTS ***
27
28 ! sb CPD=1005.7
29 ! sb CPV=1870.0
30 ! sb CL=4190.0
31 ! sb CPVMCL=2320.0
32 ! sb RV=461.5
33 ! sb RD=287.04
34 ! sb EPS=RD/RV
35 ! sb ALV0=2.501E6
36 ! cccccccccccccccccccccc
37 ! constantes coherentes avec le modele du Centre Europeen
38 ! sb RD = 1000.0 * 1.380658E-23 * 6.0221367E+23 / 28.9644
39 ! sb RV = 1000.0 * 1.380658E-23 * 6.0221367E+23 / 18.0153
40 ! sb CPD = 3.5 * RD
41 ! sb CPV = 4.0 * RV
42 ! sb CL = 4218.0
43 ! sb CI=2090.0
44 ! sb CPVMCL=CL-CPV
45 ! sb CLMCI=CL-CI
46 ! sb EPS=RD/RV
47 ! sb ALV0=2.5008E+06
48 ! sb ALF0=3.34E+05
49
50 ! ccccccccccc
51 ! on utilise les constantes thermo du Centre Europeen: (SB)
52
53 include "YOMCST.h"
54 gravity = rg !sb: Pr que gravite ne devienne pas humidite!
55
56 cpd = rcpd
57 cpv = rcpv
58 cl = rcw
59 ci = rcs
60 cpvmcl = cl - cpv
61 clmci = cl - ci
62 eps = rd/rv
63 alv0 = rlvtt
64 alf0 = rlmlt ! (ALF0 = RLSTT-RLVTT)
65
66 ! ccccccccccccccccccccc
67
68 ! *** CALCULATE CERTAIN PARCEL QUANTITIES, INCLUDING STATIC ENERGY ***
69
70 icb1 = max(icb, 2)
71 icb1 = min(icb, nl)
72
73 ! jyg1
74 ! C CPP=CPD*(1.-RR(1))+RR(1)*CPV
75 cpp = cpd*(1.-rr(nk)) + rr(nk)*cpv
76 ! jyg2
77 cpinv = 1./cpp
78 ! jyg1
79 ! ICB may be below condensation level
80 ! CC DO 100 I=1,ICB1-1
81 ! CC TPK(I)=T(1)-GZ(I)*CPINV
82 ! CC TVP(I)=TPK(I)*(1.+RR(1)/EPS)
83 DO i = 1, icb1
84 clw(i) = 0.0
85 END DO
86
87 DO i = nk, icb1
88 tpk(i) = t(nk) - (gz(i)-gz(nk))*cpinv
89 ! jyg1
90 ! CC TVP(I)=TPK(I)*(1.+RR(NK)/EPS)
91 tvp(i) = tpk(i)*(1.+rr(nk)/eps-rr(nk))
92 ! jyg2
93 dtvpdt1(i) = 1. + rr(nk)/eps - rr(nk)
94 dtvpdq1(i) = tpk(i)*(1./eps-1.)
95
96 ! jyg2
97
98 END DO
99
100
101 ! *** FIND LIFTED PARCEL TEMPERATURE AND MIXING RATIO ***
102
103 ! jyg1
104 ! C AH0=(CPD*(1.-RR(1))+CL*RR(1))*T(1)
105 ! C $ +RR(1)*(ALV0-CPVMCL*(T(1)-273.15))
106 ah0 = (cpd*(1.-rr(nk))+cl*rr(nk))*t(nk) + rr(nk)*(alv0-cpvmcl*(t(nk)-273.15 &
107 )) + gz(nk)
108 ! jyg2
109
110 ! jyg1
111 imin = icb1
112 ! If ICB is below LCL, start loop at ICB+1
113 IF (plcl<p(icb1)) imin = min(imin+1, nl)
114
115 ! CC DO 300 I=ICB1,NL
116 DO i = imin, nl
117 ! jyg2
118 alv = alv0 - cpvmcl*(t(i)-273.15)
119 alf = alf0 + clmci*(t(i)-273.15)
120
121 rg = rs(i)
122 tg = t(i)
123 ! S=CPD+ALV*ALV*RG/(RV*T(I)*T(I))
124 ! jyg1
125 ! C S=CPD*(1.-RR(1))+CL*RR(1)+ALV*ALV*RG/(RV*T(I)*T(I))
126 s = cpd*(1.-rr(nk)) + cl*rr(nk) + alv*alv*rg/(rv*t(i)*t(i))
127 ! jyg2
128 s = 1./s
129
130 DO j = 1, 2
131 ! jyg1
132 ! C AHG=CPD*TG+(CL-CPD)*RR(1)*TG+ALV*RG+GZ(I)
133 ahg = cpd*tg + (cl-cpd)*rr(nk)*tg + alv*rg + gz(i)
134 ! jyg2
135 tg = tg + s*(ah0-ahg)
136 tc = tg - 273.15
137 denom = 243.5 + tc
138 denom = max(denom, 1.0)
139
140 ! FORMULE DE BOLTON POUR PSAT
141
142 es = 6.112*exp(17.67*tc/denom)
143 rg = eps*es/(p(i)-es*(1.-eps))
144
145
146 END DO
147
148 ! jyg1
149 ! C TPK(I)=(AH0-GZ(I)-ALV*RG)/(CPD+(CL-CPD)*RR(1))
150 tpk(i) = (ah0-gz(i)-alv*rg)/(cpd+(cl-cpd)*rr(nk))
151 ! jyg2
152 ! TPK(I)=(AH0-GZ(I)-ALV*RG-(CL-CPD)*T(I)*RR(1))/CPD
153
154 ! jyg1
155 ! C CLW(I)=RR(1)-RG
156 clw(i) = rr(nk) - rg
157 ! jyg2
158 clw(i) = max(0.0, clw(i))
159 ! jyg1
160 ! CC TVP(I)=TPK(I)*(1.+RG/EPS)
161 tvp(i) = tpk(i)*(1.+rg/eps-rr(nk))
162 ! jyg2
163
164 ! jyg1 Derivatives
165
166 dtpdt1(i) = cpd*s
167 dtpdq1(i) = alv*s
168
169 dtvpdt1(i) = dtpdt1(i)*(1.+rg/eps-rr(nk)+alv*rg/(rd*tpk(i)))
170 dtvpdq1(i) = dtpdq1(i)*(1.+rg/eps-rr(nk)+alv*rg/(rd*tpk(i))) - tpk(i)
171
172 ! jyg2
173
174 END DO
175
176 ice_conv = .FALSE.
177
178 IF (ice_conv) THEN
179
180 ! JAM
181 ! RAJOUT DE LA PROCEDURE ICEFRAC
182
183 ! sb CALL ICEFRAC(T,CLW,CLW_NEW,QI,ND,NL)
184
185 DO i = icb1, nl
186 IF (t(i)<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))
197 ! CC SNEW=
198 ! CPD*(1.-RR(1))+CL*RR(1)+ALV*ALV*QSAT_NEW/(RV*TPK(I)*TPK(I))
199 snew = cpd*(1.-rr(nk)) + cl*rr(nk) + alv*alv*qsat_new/(rv*tpk(i)* &
200 tpk(i))
201
202 snew = 1./snew
203 tpk(i) = tg + (alf*qi(i)+alv*rg*(1.-(esi/es)))*snew
204 ! @$$ PRINT*,'################################'
205 ! @$$ PRINT*,TPK(I)
206 ! @$$ PRINT*,(ALF*QI(I)+ALV*RG*(1.-(ESI/ES)))*SNEW
207 END DO
208 ! CC CLW(I)=RR(1)-QSAT_NEW
209 clw(i) = rr(nk) - qsat_new
210 clw(i) = max(0.0, clw(i))
211 ! jyg1
212 ! CC TVP(I)=TPK(I)*(1.+QSAT_NEW/EPS)
213 tvp(i) = tpk(i)*(1.+qsat_new/eps-rr(nk))
214 ! jyg2
215 ELSE
216 CONTINUE
217 END IF
218
219 END DO
220
221 END IF
222
223
224 ! *****************************************************
225 ! * BK : RAJOUT DE LA TEMPERATURE DES ASCENDANCES
226 ! * NON DILUES AU NIVEAU KLEV = ND
227 ! * POSONS LE ENVIRON EGAL A CELUI DE KLEV-1
228 ! *******************************************************
229
230 tpk(nl+1) = tpk(nl)
231
232 ! ******************************************************
233
234 rg = gravity ! RG redevient la gravite de YOMCST (sb)
235
236
237 RETURN
238 END SUBROUTINE tlift
239
240
241
242
243
244
245
246
247