Line |
Branch |
Exec |
Source |
1 |
|
|
! |
2 |
|
|
! $Id: thermcell_plume.F90 3074 2017-11-15 13:31:44Z fhourdin $ |
3 |
|
|
! |
4 |
|
✗ |
SUBROUTINE thermcell_plume(itap,ngrid,klev,ptimestep,ztv,zthl,po,zl,rhobarz, & |
5 |
|
✗ |
& zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot, & |
6 |
|
✗ |
& lalim,f0,detr_star,entr_star,f_star,csc,ztva, & |
7 |
|
✗ |
& ztla,zqla,zqta,zha,zw2,w_est,ztva_est,zqsatth,lmix,lmix_bis,linter & |
8 |
|
|
& ,lev_out,lunout1,igout) |
9 |
|
|
! & ,lev_out,lunout1,igout,zbuoy,zbuoyjam) |
10 |
|
|
!-------------------------------------------------------------------------- |
11 |
|
|
! Auhtors : Catherine Rio, Frédéric Hourdin, Arnaud Jam |
12 |
|
|
! |
13 |
|
|
!thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance |
14 |
|
|
! This versions starts from a cleaning of thermcell_plume_6A (2019/01/20) |
15 |
|
|
! thermcell_plume_6A is activate for flag_thermas_ed < 10 |
16 |
|
|
! thermcell_plume_5B for flag_thermas_ed < 20 |
17 |
|
|
! thermcell_plume for flag_thermals_ed>= 20 |
18 |
|
|
! Various options are controled by the flag_thermals_ed parameter |
19 |
|
|
! = 20 : equivalent to thermcell_plume_6A with flag_thermals_ed=8 |
20 |
|
|
! = 21 : the Jam strato-cumulus modif is not activated in detrainment |
21 |
|
|
! = 29 : an other way to compute the modified buoyancy (to be tested) |
22 |
|
|
!-------------------------------------------------------------------------- |
23 |
|
|
USE IOIPSL, ONLY : getin |
24 |
|
|
USE ioipsl_getin_p_mod, ONLY : getin_p |
25 |
|
|
|
26 |
|
|
USE print_control_mod, ONLY: prt_level |
27 |
|
|
IMPLICIT NONE |
28 |
|
|
|
29 |
|
|
! |
30 |
|
|
! $Header$ |
31 |
|
|
! |
32 |
|
|
! ATTENTION!!!!: ce fichier include est compatible format fixe/format libre |
33 |
|
|
! veillez � n'utiliser que des ! pour les commentaires |
34 |
|
|
! et � bien positionner les & des lignes de continuation |
35 |
|
|
! (les placer en colonne 6 et en colonne 73) |
36 |
|
|
! |
37 |
|
|
! |
38 |
|
|
! A1.0 Fundamental constants |
39 |
|
|
REAL RPI,RCLUM,RHPLA,RKBOL,RNAVO |
40 |
|
|
! A1.1 Astronomical constants |
41 |
|
|
REAL RDAY,REA,REPSM,RSIYEA,RSIDAY,ROMEGA |
42 |
|
|
! A1.1.bis Constantes concernant l'orbite de la Terre: |
43 |
|
|
REAL R_ecc, R_peri, R_incl |
44 |
|
|
! A1.2 Geoide |
45 |
|
|
REAL RA,RG,R1SA |
46 |
|
|
! A1.3 Radiation |
47 |
|
|
! REAL RSIGMA,RI0 |
48 |
|
|
REAL RSIGMA |
49 |
|
|
! A1.4 Thermodynamic gas phase |
50 |
|
|
REAL RMO3,RMCO2,RMC,RMCH4,RMN2O,RMCFC11,RMCFC12 |
51 |
|
|
REAL R,RMD,RMV,RD,RV,RCPD,RCPV,RCVD,RCVV |
52 |
|
|
REAL RKAPPA,RETV, eps_w |
53 |
|
|
! A1.5,6 Thermodynamic liquid,solid phases |
54 |
|
|
REAL RCW,RCS |
55 |
|
|
! A1.7 Thermodynamic transition of phase |
56 |
|
|
REAL RLVTT,RLSTT,RLMLT,RTT,RATM |
57 |
|
|
! A1.8 Curve of saturation |
58 |
|
|
REAL RESTT,RALPW,RBETW,RGAMW,RALPS,RBETS,RGAMS |
59 |
|
|
REAL RALPD,RBETD,RGAMD |
60 |
|
|
! |
61 |
|
|
COMMON/YOMCST/RPI ,RCLUM ,RHPLA ,RKBOL ,RNAVO & |
62 |
|
|
& ,RDAY ,REA ,REPSM ,RSIYEA,RSIDAY,ROMEGA & |
63 |
|
|
& ,R_ecc, R_peri, R_incl & |
64 |
|
|
& ,RA ,RG ,R1SA & |
65 |
|
|
& ,RSIGMA & |
66 |
|
|
& ,R ,RMD ,RMV ,RD ,RV ,RCPD & |
67 |
|
|
& ,RMO3 ,RMCO2 ,RMC ,RMCH4 ,RMN2O ,RMCFC11 ,RMCFC12 & |
68 |
|
|
& ,RCPV ,RCVD ,RCVV ,RKAPPA,RETV, eps_w & |
69 |
|
|
& ,RCW ,RCS & |
70 |
|
|
& ,RLVTT ,RLSTT ,RLMLT ,RTT ,RATM & |
71 |
|
|
& ,RESTT ,RALPW ,RBETW ,RGAMW ,RALPS ,RBETS ,RGAMS & |
72 |
|
|
& ,RALPD ,RBETD ,RGAMD |
73 |
|
|
! ------------------------------------------------------------------ |
74 |
|
|
!$OMP THREADPRIVATE(/YOMCST/) |
75 |
|
|
! |
76 |
|
|
! $Id: YOETHF.h 2799 2017-02-24 18:50:33Z jyg $ |
77 |
|
|
! |
78 |
|
|
! ATTENTION!!!!: ce fichier include est compatible format fixe/format libre |
79 |
|
|
! veillez n'utiliser que des ! pour les commentaires |
80 |
|
|
! et bien positionner les & des lignes de continuation |
81 |
|
|
! (les placer en colonne 6 et en colonne 73) |
82 |
|
|
! |
83 |
|
|
!* COMMON *YOETHF* DERIVED CONSTANTS SPECIFIC TO ECMWF THERMODYNAMICS |
84 |
|
|
! |
85 |
|
|
! *R__ES* *CONSTANTS USED FOR COMPUTATION OF SATURATION |
86 |
|
|
! MIXING RATIO OVER LIQUID WATER(*R_LES*) OR |
87 |
|
|
! ICE(*R_IES*). |
88 |
|
|
! *RVTMP2* *RVTMP2=RCPV/RCPD-1. |
89 |
|
|
! *RHOH2O* *DENSITY OF LIQUID WATER. (RATM/100.) |
90 |
|
|
! |
91 |
|
|
REAL R2ES, R3LES, R3IES, R4LES, R4IES, R5LES, R5IES |
92 |
|
|
REAL RVTMP2, RHOH2O |
93 |
|
|
REAL R5ALVCP,R5ALSCP,RALVDCP,RALSDCP,RALFDCP,RTWAT,RTBER,RTBERCU |
94 |
|
|
REAL RTICE,RTICECU,RTWAT_RTICE_R,RTWAT_RTICECU_R,RKOOP1,RKOOP2 |
95 |
|
|
LOGICAL OK_BAD_ECMWF_THERMO ! If TRUE, then variables set by rrtm/suphec.F90 |
96 |
|
|
! If FALSE, then variables set by suphel.F90 |
97 |
|
|
COMMON /YOETHF/R2ES, R3LES, R3IES, R4LES, R4IES, R5LES, R5IES, & |
98 |
|
|
& RVTMP2, RHOH2O, & |
99 |
|
|
& R5ALVCP,R5ALSCP,RALVDCP,RALSDCP, & |
100 |
|
|
& RALFDCP,RTWAT,RTBER,RTBERCU, & |
101 |
|
|
& RTICE,RTICECU,RTWAT_RTICE_R,RTWAT_RTICECU_R,RKOOP1,& |
102 |
|
|
& RKOOP2, & |
103 |
|
|
& OK_BAD_ECMWF_THERMO |
104 |
|
|
|
105 |
|
|
!$OMP THREADPRIVATE(/YOETHF/) |
106 |
|
|
! |
107 |
|
|
! $Header$ |
108 |
|
|
! |
109 |
|
|
! |
110 |
|
|
! ATTENTION!!!!: ce fichier include est compatible format fixe/format libre |
111 |
|
|
! veillez n'utiliser que des ! pour les commentaires |
112 |
|
|
! et bien positionner les & des lignes de continuation |
113 |
|
|
! (les placer en colonne 6 et en colonne 73) |
114 |
|
|
! |
115 |
|
|
! ------------------------------------------------------------------ |
116 |
|
|
! This COMDECK includes the Thermodynamical functions for the cy39 |
117 |
|
|
! ECMWF Physics package. |
118 |
|
|
! Consistent with YOMCST Basic physics constants, assuming the |
119 |
|
|
! partial pressure of water vapour is given by a first order |
120 |
|
|
! Taylor expansion of Qs(T) w.r.t. to Temperature, using constants |
121 |
|
|
! in YOETHF |
122 |
|
|
! ------------------------------------------------------------------ |
123 |
|
|
REAL PTARG, PDELARG, P5ARG, PQSARG, PCOARG |
124 |
|
|
REAL FOEEW, FOEDE, qsats, qsatl, dqsats, dqsatl |
125 |
|
|
LOGICAL thermcep |
126 |
|
|
PARAMETER (thermcep=.TRUE.) |
127 |
|
|
! |
128 |
|
|
FOEEW ( PTARG,PDELARG ) = EXP ( & |
129 |
|
|
& (R3LES*(1.-PDELARG)+R3IES*PDELARG) * (PTARG-RTT) & |
130 |
|
|
& / (PTARG-(R4LES*(1.-PDELARG)+R4IES*PDELARG)) ) |
131 |
|
|
! |
132 |
|
|
FOEDE ( PTARG,PDELARG,P5ARG,PQSARG,PCOARG ) = PQSARG*PCOARG*P5ARG & |
133 |
|
|
& / (PTARG-(R4LES*(1.-PDELARG)+R4IES*PDELARG))**2 |
134 |
|
|
! |
135 |
|
|
qsats(ptarg) = 100.0 * 0.622 * 10.0 & |
136 |
|
|
& ** (2.07023 - 0.00320991 * ptarg & |
137 |
|
|
& - 2484.896 / ptarg + 3.56654 * LOG10(ptarg)) |
138 |
|
|
qsatl(ptarg) = 100.0 * 0.622 * 10.0 & |
139 |
|
|
& ** (23.8319 - 2948.964 / ptarg & |
140 |
|
|
& - 5.028 * LOG10(ptarg) & |
141 |
|
|
& - 29810.16 * EXP( - 0.0699382 * ptarg) & |
142 |
|
|
& + 25.21935 * EXP( - 2999.924 / ptarg)) |
143 |
|
|
! |
144 |
|
|
dqsats(ptarg,pqsarg) = RLVTT/RCPD*pqsarg * (3.56654/ptarg & |
145 |
|
|
& +2484.896*LOG(10.)/ptarg**2 & |
146 |
|
|
& -0.00320991*LOG(10.)) |
147 |
|
|
dqsatl(ptarg,pqsarg) = RLVTT/RCPD*pqsarg*LOG(10.)* & |
148 |
|
|
& (2948.964/ptarg**2-5.028/LOG(10.)/ptarg & |
149 |
|
|
& +25.21935*2999.924/ptarg**2*EXP(-2999.924/ptarg) & |
150 |
|
|
& +29810.16*0.0699382*EXP(-0.0699382*ptarg)) |
151 |
|
|
integer :: iflag_thermals,nsplit_thermals |
152 |
|
|
|
153 |
|
|
!!! nrlmd le 10/04/2012 |
154 |
|
|
integer :: iflag_trig_bl,iflag_clos_bl |
155 |
|
|
integer :: tau_trig_shallow,tau_trig_deep |
156 |
|
|
real :: s_trig |
157 |
|
|
!!! fin nrlmd le 10/04/2012 |
158 |
|
|
|
159 |
|
|
real,parameter :: r_aspect_thermals=2.,l_mix_thermals=30. |
160 |
|
|
real :: alp_bl_k |
161 |
|
|
real :: tau_thermals,fact_thermals_ed_dz |
162 |
|
|
integer,parameter :: w2di_thermals=0 |
163 |
|
|
integer :: isplit |
164 |
|
|
|
165 |
|
|
integer :: iflag_coupl,iflag_clos,iflag_wake |
166 |
|
|
integer :: iflag_thermals_ed,iflag_thermals_optflux,iflag_thermals_closure |
167 |
|
|
|
168 |
|
|
common/ctherm1/iflag_thermals,nsplit_thermals,iflag_thermals_closure |
169 |
|
|
common/ctherm2/tau_thermals,alp_bl_k,fact_thermals_ed_dz |
170 |
|
|
common/ctherm4/iflag_coupl,iflag_clos,iflag_wake |
171 |
|
|
common/ctherm5/iflag_thermals_ed,iflag_thermals_optflux |
172 |
|
|
|
173 |
|
|
!!! nrlmd le 10/04/2012 |
174 |
|
|
common/ctherm6/iflag_trig_bl,iflag_clos_bl |
175 |
|
|
common/ctherm7/tau_trig_shallow,tau_trig_deep |
176 |
|
|
common/ctherm8/s_trig |
177 |
|
|
!!! fin nrlmd le 10/04/2012 |
178 |
|
|
|
179 |
|
|
!$OMP THREADPRIVATE(/ctherm1/,/ctherm2/,/ctherm4/,/ctherm5/) |
180 |
|
|
!$OMP THREADPRIVATE(/ctherm6/,/ctherm7/,/ctherm8/) |
181 |
|
|
|
182 |
|
|
INTEGER itap |
183 |
|
|
INTEGER lunout1,igout |
184 |
|
|
INTEGER ngrid,klev |
185 |
|
|
REAL ptimestep |
186 |
|
|
REAL ztv(ngrid,klev) |
187 |
|
|
REAL zthl(ngrid,klev) |
188 |
|
|
REAL po(ngrid,klev) |
189 |
|
|
REAL zl(ngrid,klev) |
190 |
|
|
REAL rhobarz(ngrid,klev) |
191 |
|
|
REAL zlev(ngrid,klev+1) |
192 |
|
|
REAL pplev(ngrid,klev+1) |
193 |
|
|
REAL pphi(ngrid,klev) |
194 |
|
|
REAL zpspsk(ngrid,klev) |
195 |
|
|
REAL alim_star(ngrid,klev) |
196 |
|
|
REAL f0(ngrid) |
197 |
|
|
INTEGER lalim(ngrid) |
198 |
|
|
integer lev_out ! niveau pour les print |
199 |
|
|
integer nbpb |
200 |
|
|
|
201 |
|
|
real alim_star_tot(ngrid) |
202 |
|
|
|
203 |
|
|
REAL ztva(ngrid,klev) |
204 |
|
|
REAL ztla(ngrid,klev) |
205 |
|
|
REAL zqla(ngrid,klev) |
206 |
|
|
REAL zqta(ngrid,klev) |
207 |
|
|
REAL zha(ngrid,klev) |
208 |
|
|
|
209 |
|
|
REAL detr_star(ngrid,klev) |
210 |
|
|
REAL coefc |
211 |
|
|
REAL entr_star(ngrid,klev) |
212 |
|
✗ |
REAL detr(ngrid,klev) |
213 |
|
✗ |
REAL entr(ngrid,klev) |
214 |
|
|
|
215 |
|
|
REAL csc(ngrid,klev) |
216 |
|
|
|
217 |
|
|
REAL zw2(ngrid,klev+1) |
218 |
|
|
REAL w_est(ngrid,klev+1) |
219 |
|
|
REAL f_star(ngrid,klev+1) |
220 |
|
✗ |
REAL wa_moy(ngrid,klev+1) |
221 |
|
|
|
222 |
|
|
REAL ztva_est(ngrid,klev) |
223 |
|
✗ |
REAL ztv_est(ngrid,klev) |
224 |
|
✗ |
REAL zqla_est(ngrid,klev) |
225 |
|
|
REAL zqsatth(ngrid,klev) |
226 |
|
✗ |
REAL zta_est(ngrid,klev) |
227 |
|
✗ |
REAL ztemp(ngrid),zqsat(ngrid) |
228 |
|
|
REAL zdw2,zdw2bis |
229 |
|
|
REAL zw2modif |
230 |
|
|
REAL zw2fact,zw2factbis |
231 |
|
✗ |
REAL zeps(ngrid,klev) |
232 |
|
|
|
233 |
|
|
REAL linter(ngrid) |
234 |
|
|
INTEGER lmix(ngrid) |
235 |
|
|
INTEGER lmix_bis(ngrid) |
236 |
|
✗ |
REAL wmaxa(ngrid) |
237 |
|
|
|
238 |
|
|
INTEGER ig,l,k,lt,it,lm |
239 |
|
|
|
240 |
|
✗ |
real zdz,zbuoy(ngrid,klev),zalpha,gamma(ngrid,klev),zdqt(ngrid,klev),zw2m |
241 |
|
✗ |
real zbuoyjam(ngrid,klev),zdqtjam(ngrid,klev) |
242 |
|
|
real zdz2,zdz3,lmel,entrbis,zdzbis |
243 |
|
✗ |
real d_temp(ngrid) |
244 |
|
|
real ztv1,ztv2,factinv,zinv,zlmel |
245 |
|
|
real zlmelup,zlmeldwn,zlt,zltdwn,zltup |
246 |
|
|
real atv1,atv2,btv1,btv2 |
247 |
|
|
real ztv_est1,ztv_est2 |
248 |
|
|
real zcor,zdelta,zcvm5,qlbef |
249 |
|
|
real zbetalpha, coefzlmel |
250 |
|
|
real eps |
251 |
|
|
REAL REPS,RLvCp,DDT0 |
252 |
|
|
PARAMETER (DDT0=.01) |
253 |
|
|
logical Zsat |
254 |
|
✗ |
LOGICAL active(ngrid),activetmp(ngrid) |
255 |
|
|
REAL fact_gamma,fact_gamma2,fact_epsilon2 |
256 |
|
|
|
257 |
|
|
REAL, SAVE :: fact_epsilon=0.002 |
258 |
|
|
REAL, SAVE :: betalpha=0.9 |
259 |
|
|
REAL, SAVE :: afact=2./3. |
260 |
|
|
REAL, SAVE :: fact_shell=1. |
261 |
|
|
REAL,SAVE :: detr_min=1.e-5 |
262 |
|
|
REAL,SAVE :: entr_min=1.e-5 |
263 |
|
|
REAL,SAVE :: detr_q_coef=0.012 |
264 |
|
|
REAL,SAVE :: detr_q_power=0.5 |
265 |
|
|
REAL,SAVE :: mix0=0. |
266 |
|
|
INTEGER,SAVE :: thermals_flag_alim=0 |
267 |
|
|
|
268 |
|
|
!$OMP THREADPRIVATE(fact_epsilon, betalpha, afact, fact_shell) |
269 |
|
|
!$OMP THREADPRIVATE(detr_min, entr_min, detr_q_coef, detr_q_power) |
270 |
|
|
!$OMP THREADPRIVATE( mix0, thermals_flag_alim) |
271 |
|
|
|
272 |
|
|
LOGICAL, SAVE :: first=.true. |
273 |
|
|
!$OMP THREADPRIVATE(first) |
274 |
|
|
|
275 |
|
|
|
276 |
|
|
REAL c2(ngrid,klev) |
277 |
|
|
|
278 |
|
✗ |
if (ngrid==1) print*,'THERMCELL PLUME MODIFIE 2014/07/11' |
279 |
|
|
Zsat=.false. |
280 |
|
|
! Initialisation |
281 |
|
|
|
282 |
|
✗ |
RLvCp = RLVTT/RCPD |
283 |
|
✗ |
IF (first) THEN |
284 |
|
|
|
285 |
|
✗ |
CALL getin_p('thermals_fact_epsilon',fact_epsilon) |
286 |
|
✗ |
CALL getin_p('thermals_betalpha',betalpha) |
287 |
|
✗ |
CALL getin_p('thermals_afact',afact) |
288 |
|
✗ |
CALL getin_p('thermals_fact_shell',fact_shell) |
289 |
|
✗ |
CALL getin_p('thermals_detr_min',detr_min) |
290 |
|
✗ |
CALL getin_p('thermals_entr_min',entr_min) |
291 |
|
✗ |
CALL getin_p('thermals_detr_q_coef',detr_q_coef) |
292 |
|
✗ |
CALL getin_p('thermals_detr_q_power',detr_q_power) |
293 |
|
✗ |
CALL getin_p('thermals_mix0',mix0) |
294 |
|
✗ |
CALL getin_p('thermals_flag_alim',thermals_flag_alim) |
295 |
|
|
|
296 |
|
|
|
297 |
|
✗ |
first=.false. |
298 |
|
|
ENDIF |
299 |
|
|
|
300 |
|
✗ |
zbetalpha=betalpha/(1.+betalpha) |
301 |
|
|
|
302 |
|
|
|
303 |
|
|
! Initialisations des variables r?elles |
304 |
|
|
if (1==1) then |
305 |
|
✗ |
ztva(:,:)=ztv(:,:) |
306 |
|
✗ |
ztva_est(:,:)=ztva(:,:) |
307 |
|
✗ |
ztv_est(:,:)=ztv(:,:) |
308 |
|
✗ |
ztla(:,:)=zthl(:,:) |
309 |
|
✗ |
zqta(:,:)=po(:,:) |
310 |
|
✗ |
zqla(:,:)=0. |
311 |
|
✗ |
zha(:,:) = ztva(:,:) |
312 |
|
|
else |
313 |
|
|
ztva(:,:)=0. |
314 |
|
|
ztv_est(:,:)=0. |
315 |
|
|
ztva_est(:,:)=0. |
316 |
|
|
ztla(:,:)=0. |
317 |
|
|
zqta(:,:)=0. |
318 |
|
|
zha(:,:) =0. |
319 |
|
|
endif |
320 |
|
|
|
321 |
|
✗ |
zqla_est(:,:)=0. |
322 |
|
✗ |
zqsatth(:,:)=0. |
323 |
|
✗ |
zqla(:,:)=0. |
324 |
|
✗ |
detr_star(:,:)=0. |
325 |
|
✗ |
entr_star(:,:)=0. |
326 |
|
✗ |
alim_star(:,:)=0. |
327 |
|
✗ |
alim_star_tot(:)=0. |
328 |
|
✗ |
csc(:,:)=0. |
329 |
|
✗ |
detr(:,:)=0. |
330 |
|
✗ |
entr(:,:)=0. |
331 |
|
✗ |
zw2(:,:)=0. |
332 |
|
✗ |
zbuoy(:,:)=0. |
333 |
|
✗ |
zbuoyjam(:,:)=0. |
334 |
|
✗ |
gamma(:,:)=0. |
335 |
|
✗ |
zeps(:,:)=0. |
336 |
|
✗ |
w_est(:,:)=0. |
337 |
|
✗ |
f_star(:,:)=0. |
338 |
|
✗ |
wa_moy(:,:)=0. |
339 |
|
✗ |
linter(:)=1. |
340 |
|
|
! linter(:)=1. |
341 |
|
|
! Initialisation des variables entieres |
342 |
|
✗ |
lmix(:)=1 |
343 |
|
✗ |
lmix_bis(:)=2 |
344 |
|
✗ |
wmaxa(:)=0. |
345 |
|
|
|
346 |
|
|
|
347 |
|
|
!------------------------------------------------------------------------- |
348 |
|
|
! On ne considere comme actif que les colonnes dont les deux premieres |
349 |
|
|
! couches sont instables. |
350 |
|
|
!------------------------------------------------------------------------- |
351 |
|
|
|
352 |
|
✗ |
active(:)=ztv(:,1)>ztv(:,2) |
353 |
|
✗ |
d_temp(:)=0. ! Pour activer un contraste de temperature a la base |
354 |
|
|
! du panache |
355 |
|
|
! Cet appel pourrait être fait avant thermcell_plume dans thermcell_main |
356 |
|
✗ |
CALL thermcell_alim(thermals_flag_alim,ngrid,klev,ztv,d_temp,zlev,alim_star,lalim) |
357 |
|
|
|
358 |
|
|
!------------------------------------------------------------------------------ |
359 |
|
|
! Calcul dans la premiere couche |
360 |
|
|
! On decide dans cette version que le thermique n'est actif que si la premiere |
361 |
|
|
! couche est instable. |
362 |
|
|
! Pourrait etre change si on veut que le thermiques puisse se d??clencher |
363 |
|
|
! dans une couche l>1 |
364 |
|
|
!------------------------------------------------------------------------------ |
365 |
|
✗ |
do ig=1,ngrid |
366 |
|
|
! Le panache va prendre au debut les caracteristiques de l'air contenu |
367 |
|
|
! dans cette couche. |
368 |
|
✗ |
if (active(ig)) then |
369 |
|
✗ |
ztla(ig,1)=zthl(ig,1) |
370 |
|
✗ |
zqta(ig,1)=po(ig,1) |
371 |
|
✗ |
zqla(ig,1)=zl(ig,1) |
372 |
|
|
!cr: attention, prise en compte de f*(1)=1 |
373 |
|
✗ |
f_star(ig,2)=alim_star(ig,1) |
374 |
|
|
zw2(ig,2)=2.*RG*(ztv(ig,1)-ztv(ig,2))/ztv(ig,2) & |
375 |
|
|
& *(zlev(ig,2)-zlev(ig,1)) & |
376 |
|
✗ |
& *0.4*pphi(ig,1)/(pphi(ig,2)-pphi(ig,1)) |
377 |
|
✗ |
w_est(ig,2)=zw2(ig,2) |
378 |
|
|
endif |
379 |
|
|
enddo |
380 |
|
|
! |
381 |
|
|
|
382 |
|
|
!============================================================================== |
383 |
|
|
!boucle de calcul de la vitesse verticale dans le thermique |
384 |
|
|
!============================================================================== |
385 |
|
✗ |
do l=2,klev-1 |
386 |
|
|
!============================================================================== |
387 |
|
|
|
388 |
|
|
|
389 |
|
|
! On decide si le thermique est encore actif ou non |
390 |
|
|
! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test |
391 |
|
✗ |
do ig=1,ngrid |
392 |
|
|
active(ig)=active(ig) & |
393 |
|
|
& .and. zw2(ig,l)>1.e-10 & |
394 |
|
✗ |
& .and. f_star(ig,l)+alim_star(ig,l)>1.e-10 |
395 |
|
|
enddo |
396 |
|
|
|
397 |
|
|
|
398 |
|
|
|
399 |
|
|
!--------------------------------------------------------------------------- |
400 |
|
|
! calcul des proprietes thermodynamiques et de la vitesse de la couche l |
401 |
|
|
! sans tenir compte du detrainement et de l'entrainement dans cette |
402 |
|
|
! couche |
403 |
|
|
! C'est a dire qu'on suppose |
404 |
|
|
! ztla(l)=ztla(l-1) et zqta(l)=zqta(l-1) |
405 |
|
|
! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer |
406 |
|
|
! avant) a l'alimentation pour avoir un calcul plus propre |
407 |
|
|
!--------------------------------------------------------------------------- |
408 |
|
|
|
409 |
|
✗ |
ztemp(:)=zpspsk(:,l)*ztla(:,l-1) |
410 |
|
✗ |
call thermcell_qsat(ngrid,active,pplev(:,l),ztemp,zqta(:,l-1),zqsat(:)) |
411 |
|
✗ |
do ig=1,ngrid |
412 |
|
|
! print*,'active',active(ig),ig,l |
413 |
|
✗ |
if(active(ig)) then |
414 |
|
✗ |
zqla_est(ig,l)=max(0.,zqta(ig,l-1)-zqsat(ig)) |
415 |
|
✗ |
ztva_est(ig,l) = ztla(ig,l-1)*zpspsk(ig,l)+RLvCp*zqla_est(ig,l) |
416 |
|
✗ |
zta_est(ig,l)=ztva_est(ig,l) |
417 |
|
✗ |
ztva_est(ig,l) = ztva_est(ig,l)/zpspsk(ig,l) |
418 |
|
|
ztva_est(ig,l) = ztva_est(ig,l)*(1.+RETV*(zqta(ig,l-1) & |
419 |
|
✗ |
& -zqla_est(ig,l))-zqla_est(ig,l)) |
420 |
|
|
|
421 |
|
|
|
422 |
|
|
!Modif AJAM |
423 |
|
|
|
424 |
|
✗ |
zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l) |
425 |
|
✗ |
zdz=zlev(ig,l+1)-zlev(ig,l) |
426 |
|
✗ |
lmel=fact_thermals_ed_dz*zlev(ig,l) |
427 |
|
|
! lmel=0.09*zlev(ig,l) |
428 |
|
✗ |
zlmel=zlev(ig,l)+lmel |
429 |
|
✗ |
zlmelup=zlmel+(zdz/2) |
430 |
|
|
zlmeldwn=zlmel-(zdz/2) |
431 |
|
|
|
432 |
|
|
lt=l+1 |
433 |
|
|
zlt=zlev(ig,lt) |
434 |
|
|
zdz3=zlev(ig,lt+1)-zlt |
435 |
|
|
zltdwn=zlt-zdz3/2 |
436 |
|
|
zltup=zlt+zdz3/2 |
437 |
|
|
|
438 |
|
|
!========================================================================= |
439 |
|
|
! 3. Calcul de la flotabilite modifie par melange avec l'air au dessus |
440 |
|
|
!========================================================================= |
441 |
|
|
|
442 |
|
|
!-------------------------------------------------- |
443 |
|
|
lt=l+1 |
444 |
|
|
zlt=zlev(ig,lt) |
445 |
|
|
zdz2=zlev(ig,lt)-zlev(ig,l) |
446 |
|
|
|
447 |
|
✗ |
do while (lmel.gt.zdz2) |
448 |
|
✗ |
lt=lt+1 |
449 |
|
✗ |
zlt=zlev(ig,lt) |
450 |
|
✗ |
zdz2=zlt-zlev(ig,l) |
451 |
|
|
enddo |
452 |
|
✗ |
zdz3=zlev(ig,lt+1)-zlt |
453 |
|
✗ |
zltdwn=zlev(ig,lt)-zdz3/2 |
454 |
|
|
zlmelup=zlmel+(zdz/2) |
455 |
|
✗ |
coefzlmel=Min(1.,(zlmelup-zltdwn)/zdz) |
456 |
|
|
zbuoyjam(ig,l)=1.*RG*(coefzlmel*(ztva_est(ig,l)- & |
457 |
|
|
& ztv(ig,lt))/ztv(ig,lt)+(1.-coefzlmel)*(ztva_est(ig,l)- & |
458 |
|
✗ |
& ztv(ig,lt-1))/ztv(ig,lt-1))+0.*zbuoy(ig,l) |
459 |
|
|
|
460 |
|
|
!------------------------------------------------ |
461 |
|
|
!AJAM:nouveau calcul de w? |
462 |
|
|
!------------------------------------------------ |
463 |
|
|
zdz=zlev(ig,l+1)-zlev(ig,l) |
464 |
|
|
zdzbis=zlev(ig,l)-zlev(ig,l-1) |
465 |
|
|
zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l) |
466 |
|
✗ |
zw2fact=fact_epsilon*2.*zdz/(1.+betalpha) |
467 |
|
|
zw2factbis=fact_epsilon*2.*zdzbis/(1.+betalpha) |
468 |
|
✗ |
zdw2=afact*zbuoy(ig,l)/fact_epsilon |
469 |
|
|
zdw2bis=afact*zbuoy(ig,l-1)/fact_epsilon |
470 |
|
|
lm=Max(1,l-2) |
471 |
|
✗ |
w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2)+zdw2) |
472 |
|
|
endif |
473 |
|
|
enddo |
474 |
|
|
|
475 |
|
|
|
476 |
|
|
!------------------------------------------------- |
477 |
|
|
!calcul des taux d'entrainement et de detrainement |
478 |
|
|
!------------------------------------------------- |
479 |
|
|
|
480 |
|
✗ |
do ig=1,ngrid |
481 |
|
✗ |
if (active(ig)) then |
482 |
|
|
|
483 |
|
|
! zw2m=max(0.5*(w_est(ig,l)+w_est(ig,l+1)),0.1) |
484 |
|
✗ |
zw2m=w_est(ig,l+1) |
485 |
|
✗ |
zdz=zlev(ig,l+1)-zlev(ig,l) |
486 |
|
✗ |
zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l) |
487 |
|
✗ |
zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l+1))/rhobarz(ig,l) |
488 |
|
✗ |
zdqt(ig,l)=max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l) |
489 |
|
|
|
490 |
|
|
!========================================================================= |
491 |
|
|
! 4. Calcul de l'entrainement et du detrainement |
492 |
|
|
!========================================================================= |
493 |
|
|
|
494 |
|
|
detr_star(ig,l)=f_star(ig,l)*zdz & |
495 |
|
|
& *( mix0 * 0.1 / (zalpha+0.001) & |
496 |
|
|
& + MAX(detr_min, -afact*zbetalpha*zbuoyjam(ig,l)/zw2m & |
497 |
|
✗ |
& + detr_q_coef*(zdqt(ig,l)/zw2m)**detr_q_power)) |
498 |
|
|
|
499 |
|
✗ |
if ( iflag_thermals_ed == 20 ) then |
500 |
|
|
entr_star(ig,l)=f_star(ig,l)*zdz* ( & |
501 |
|
|
& mix0 * 0.1 / (zalpha+0.001) & |
502 |
|
|
& + zbetalpha*MAX(entr_min, & |
503 |
|
✗ |
& afact*zbuoyjam(ig,l)/zw2m - fact_epsilon)) |
504 |
|
|
else |
505 |
|
|
entr_star(ig,l)=f_star(ig,l)*zdz* ( & |
506 |
|
|
& mix0 * 0.1 / (zalpha+0.001) & |
507 |
|
|
& + zbetalpha*MAX(entr_min, & |
508 |
|
✗ |
& afact*zbuoy(ig,l)/zw2m - fact_epsilon)) |
509 |
|
|
endif |
510 |
|
|
|
511 |
|
|
! En dessous de lalim, on prend le max de alim_star et entr_star pour |
512 |
|
|
! alim_star et 0 sinon |
513 |
|
✗ |
if (l.lt.lalim(ig)) then |
514 |
|
✗ |
alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l)) |
515 |
|
✗ |
entr_star(ig,l)=0. |
516 |
|
|
endif |
517 |
|
|
f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l) & |
518 |
|
✗ |
& -detr_star(ig,l) |
519 |
|
|
|
520 |
|
|
endif |
521 |
|
|
enddo |
522 |
|
|
|
523 |
|
|
|
524 |
|
|
!============================================================================ |
525 |
|
|
! 5. calcul de la vitesse verticale en melangeant Tl et qt du thermique |
526 |
|
|
!=========================================================================== |
527 |
|
|
|
528 |
|
✗ |
activetmp(:)=active(:) .and. f_star(:,l+1)>1.e-10 |
529 |
|
✗ |
do ig=1,ngrid |
530 |
|
✗ |
if (activetmp(ig)) then |
531 |
|
|
Zsat=.false. |
532 |
|
|
ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+ & |
533 |
|
|
& (alim_star(ig,l)+entr_star(ig,l))*zthl(ig,l)) & |
534 |
|
✗ |
& /(f_star(ig,l+1)+detr_star(ig,l)) |
535 |
|
|
zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+ & |
536 |
|
|
& (alim_star(ig,l)+entr_star(ig,l))*po(ig,l)) & |
537 |
|
✗ |
& /(f_star(ig,l+1)+detr_star(ig,l)) |
538 |
|
|
|
539 |
|
|
endif |
540 |
|
|
enddo |
541 |
|
|
|
542 |
|
✗ |
ztemp(:)=zpspsk(:,l)*ztla(:,l) |
543 |
|
✗ |
call thermcell_qsat(ngrid,activetmp,pplev(:,l),ztemp,zqta(:,l),zqsatth(:,l)) |
544 |
|
✗ |
do ig=1,ngrid |
545 |
|
✗ |
if (activetmp(ig)) then |
546 |
|
|
! on ecrit de maniere conservative (sat ou non) |
547 |
|
|
! T = Tl +Lv/Cp ql |
548 |
|
✗ |
zqla(ig,l)=max(0.,zqta(ig,l)-zqsatth(ig,l)) |
549 |
|
✗ |
ztva(ig,l) = ztla(ig,l)*zpspsk(ig,l)+RLvCp*zqla(ig,l) |
550 |
|
✗ |
ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l) |
551 |
|
|
!on rajoute le calcul de zha pour diagnostiques (temp potentielle) |
552 |
|
✗ |
zha(ig,l) = ztva(ig,l) |
553 |
|
|
ztva(ig,l) = ztva(ig,l)*(1.+RETV*(zqta(ig,l) & |
554 |
|
✗ |
& -zqla(ig,l))-zqla(ig,l)) |
555 |
|
✗ |
zbuoy(ig,l)=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l) |
556 |
|
✗ |
zdz=zlev(ig,l+1)-zlev(ig,l) |
557 |
|
|
zdzbis=zlev(ig,l)-zlev(ig,l-1) |
558 |
|
✗ |
zeps(ig,l)=(entr_star(ig,l)+alim_star(ig,l))/(f_star(ig,l)*zdz) |
559 |
|
✗ |
zw2fact=fact_epsilon*2.*zdz/(1.+betalpha) |
560 |
|
|
zw2factbis=fact_epsilon*2.*zdzbis/(1.+betalpha) |
561 |
|
✗ |
zdw2= afact*zbuoy(ig,l)/(fact_epsilon) |
562 |
|
|
zdw2bis= afact*zbuoy(ig,l-1)/(fact_epsilon) |
563 |
|
✗ |
zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2)+zdw2) |
564 |
|
|
endif |
565 |
|
|
enddo |
566 |
|
|
|
567 |
|
✗ |
if (prt_level.ge.20) print*,'coucou calcul detr 460: ig, l',ig, l |
568 |
|
|
! |
569 |
|
|
!=========================================================================== |
570 |
|
|
! 6. initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max |
571 |
|
|
!=========================================================================== |
572 |
|
|
|
573 |
|
✗ |
nbpb=0 |
574 |
|
✗ |
do ig=1,ngrid |
575 |
|
✗ |
if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then |
576 |
|
|
! stop'On tombe sur le cas particulier de thermcell_dry' |
577 |
|
|
! print*,'On tombe sur le cas particulier de thermcell_plume' |
578 |
|
✗ |
nbpb=nbpb+1 |
579 |
|
✗ |
zw2(ig,l+1)=0. |
580 |
|
✗ |
linter(ig)=l+1 |
581 |
|
|
endif |
582 |
|
|
|
583 |
|
✗ |
if (zw2(ig,l+1).lt.0.) then |
584 |
|
|
linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l)) & |
585 |
|
✗ |
& -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l)) |
586 |
|
✗ |
zw2(ig,l+1)=0. |
587 |
|
|
!+CR:04/05/12:correction calcul linter pour calcul de zmax continu |
588 |
|
✗ |
elseif (f_star(ig,l+1).lt.0.) then |
589 |
|
|
linter(ig)=(l*(f_star(ig,l+1)-f_star(ig,l)) & |
590 |
|
✗ |
& -f_star(ig,l))/(f_star(ig,l+1)-f_star(ig,l)) |
591 |
|
✗ |
zw2(ig,l+1)=0. |
592 |
|
|
!fin CR:04/05/12 |
593 |
|
|
endif |
594 |
|
|
|
595 |
|
✗ |
wa_moy(ig,l+1)=sqrt(zw2(ig,l+1)) |
596 |
|
|
|
597 |
|
✗ |
if (wa_moy(ig,l+1).gt.wmaxa(ig)) then |
598 |
|
|
! lmix est le niveau de la couche ou w (wa_moy) est maximum |
599 |
|
|
!on rajoute le calcul de lmix_bis |
600 |
|
✗ |
if (zqla(ig,l).lt.1.e-10) then |
601 |
|
✗ |
lmix_bis(ig)=l+1 |
602 |
|
|
endif |
603 |
|
✗ |
lmix(ig)=l+1 |
604 |
|
✗ |
wmaxa(ig)=wa_moy(ig,l+1) |
605 |
|
|
endif |
606 |
|
|
enddo |
607 |
|
|
|
608 |
|
✗ |
if (nbpb>0) then |
609 |
|
✗ |
print*,'WARNING on tombe ',nbpb,' x sur un pb pour l=',l,' dans thermcell_plume' |
610 |
|
|
endif |
611 |
|
|
|
612 |
|
|
!========================================================================= |
613 |
|
|
! FIN DE LA BOUCLE VERTICALE |
614 |
|
|
enddo |
615 |
|
|
!========================================================================= |
616 |
|
|
|
617 |
|
|
!on recalcule alim_star_tot |
618 |
|
✗ |
do ig=1,ngrid |
619 |
|
✗ |
alim_star_tot(ig)=0. |
620 |
|
|
enddo |
621 |
|
✗ |
do ig=1,ngrid |
622 |
|
✗ |
do l=1,lalim(ig)-1 |
623 |
|
✗ |
alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l) |
624 |
|
|
enddo |
625 |
|
|
enddo |
626 |
|
|
|
627 |
|
|
|
628 |
|
✗ |
if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l |
629 |
|
|
|
630 |
|
|
|
631 |
|
|
|
632 |
|
✗ |
return |
633 |
|
|
end |
634 |
|
|
|