1 |
|
|
MODULE lmdz_thermcell_plume |
2 |
|
|
! |
3 |
|
|
! $Id: thermcell_plume.F90 3074 2017-11-15 13:31:44Z fhourdin $ |
4 |
|
|
! |
5 |
|
|
CONTAINS |
6 |
|
|
|
7 |
|
|
SUBROUTINE thermcell_plume(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz, & |
8 |
|
|
& zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot, & |
9 |
|
|
& lalim,f0,detr_star,entr_star,f_star,csc,ztva, & |
10 |
|
|
& ztla,zqla,zqta,zha,zw2,w_est,ztva_est,zqsatth,lmix,lmix_bis,linter & |
11 |
|
|
& ,lev_out,lunout1,igout) |
12 |
|
|
! & ,lev_out,lunout1,igout,zbuoy,zbuoyjam) |
13 |
|
|
!-------------------------------------------------------------------------- |
14 |
|
|
! Auhtors : Catherine Rio, Frédéric Hourdin, Arnaud Jam |
15 |
|
|
! |
16 |
|
|
!thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance |
17 |
|
|
! This versions starts from a cleaning of thermcell_plume_6A (2019/01/20) |
18 |
|
|
! thermcell_plume_6A is activate for flag_thermas_ed < 10 |
19 |
|
|
! thermcell_plume_5B for flag_thermas_ed < 20 |
20 |
|
|
! thermcell_plume for flag_thermals_ed>= 20 |
21 |
|
|
! Various options are controled by the flag_thermals_ed parameter |
22 |
|
|
! = 20 : equivalent to thermcell_plume_6A with flag_thermals_ed=8 |
23 |
|
|
! = 21 : the Jam strato-cumulus modif is not activated in detrainment |
24 |
|
|
! = 29 : an other way to compute the modified buoyancy (to be tested) |
25 |
|
|
!-------------------------------------------------------------------------- |
26 |
|
|
|
27 |
|
|
USE lmdz_thermcell_ini, ONLY: prt_level,fact_thermals_ed_dz,iflag_thermals_ed,RLvCP,RETV,RG |
28 |
|
|
USE lmdz_thermcell_ini, ONLY: fact_epsilon, betalpha, afact, fact_shell |
29 |
|
|
USE lmdz_thermcell_ini, ONLY: detr_min, entr_min, detr_q_coef, detr_q_power |
30 |
|
|
USE lmdz_thermcell_ini, ONLY: mix0, thermals_flag_alim |
31 |
|
|
USE lmdz_thermcell_alim, ONLY : thermcell_alim |
32 |
|
|
USE lmdz_thermcell_qsat, ONLY : thermcell_qsat |
33 |
|
|
|
34 |
|
|
|
35 |
|
|
IMPLICIT NONE |
36 |
|
|
|
37 |
|
|
integer,intent(in) :: itap,lev_out,lunout1,igout,ngrid,nlay |
38 |
|
|
real,intent(in) :: ptimestep |
39 |
|
|
real,intent(in),dimension(ngrid,nlay) :: ztv |
40 |
|
|
real,intent(in),dimension(ngrid,nlay) :: zthl |
41 |
|
|
real,intent(in),dimension(ngrid,nlay) :: po |
42 |
|
|
real,intent(in),dimension(ngrid,nlay) :: zl |
43 |
|
|
real,intent(in),dimension(ngrid,nlay) :: rhobarz |
44 |
|
|
real,intent(in),dimension(ngrid,nlay+1) :: zlev |
45 |
|
|
real,intent(in),dimension(ngrid,nlay+1) :: pplev |
46 |
|
|
real,intent(in),dimension(ngrid,nlay) :: pphi |
47 |
|
|
real,intent(in),dimension(ngrid,nlay) :: zpspsk |
48 |
|
|
real,intent(in),dimension(ngrid) :: f0 |
49 |
|
|
|
50 |
|
|
integer,intent(out) :: lalim(ngrid) |
51 |
|
|
real,intent(out),dimension(ngrid,nlay) :: alim_star |
52 |
|
|
real,intent(out),dimension(ngrid) :: alim_star_tot |
53 |
|
|
real,intent(out),dimension(ngrid,nlay) :: detr_star |
54 |
|
|
real,intent(out),dimension(ngrid,nlay) :: entr_star |
55 |
|
|
real,intent(out),dimension(ngrid,nlay+1) :: f_star |
56 |
|
|
real,intent(out),dimension(ngrid,nlay) :: csc |
57 |
|
|
real,intent(out),dimension(ngrid,nlay) :: ztva |
58 |
|
|
real,intent(out),dimension(ngrid,nlay) :: ztla |
59 |
|
|
real,intent(out),dimension(ngrid,nlay) :: zqla |
60 |
|
|
real,intent(out),dimension(ngrid,nlay) :: zqta |
61 |
|
|
real,intent(out),dimension(ngrid,nlay) :: zha |
62 |
|
|
real,intent(out),dimension(ngrid,nlay+1) :: zw2 |
63 |
|
|
real,intent(out),dimension(ngrid,nlay+1) :: w_est |
64 |
|
|
real,intent(out),dimension(ngrid,nlay) :: ztva_est |
65 |
|
|
real,intent(out),dimension(ngrid,nlay) :: zqsatth |
66 |
|
|
integer,intent(out),dimension(ngrid) :: lmix(ngrid) |
67 |
|
|
integer,intent(out),dimension(ngrid) :: lmix_bis(ngrid) |
68 |
|
|
real,intent(out),dimension(ngrid) :: linter(ngrid) |
69 |
|
|
|
70 |
|
|
|
71 |
|
|
REAL,dimension(ngrid,nlay+1) :: wa_moy |
72 |
|
|
REAL,dimension(ngrid,nlay) :: entr,detr |
73 |
|
|
REAL,dimension(ngrid,nlay) :: ztv_est |
74 |
|
|
REAL,dimension(ngrid,nlay) :: zqla_est |
75 |
|
|
REAL,dimension(ngrid,nlay) :: zta_est |
76 |
|
|
REAL,dimension(ngrid) :: ztemp,zqsat |
77 |
|
|
REAL zdw2,zdw2bis |
78 |
|
|
REAL zw2modif |
79 |
|
|
REAL zw2fact,zw2factbis |
80 |
|
|
REAL,dimension(ngrid,nlay) :: zeps |
81 |
|
|
|
82 |
|
|
REAL,dimension(ngrid) :: wmaxa |
83 |
|
|
|
84 |
|
|
INTEGER ig,l,k,lt,it,lm,nbpb |
85 |
|
|
|
86 |
|
|
real,dimension(ngrid,nlay) :: zbuoy,gamma,zdqt |
87 |
|
|
real zdz,zalpha,zw2m |
88 |
|
|
real,dimension(ngrid,nlay) :: zbuoyjam,zdqtjam |
89 |
|
|
real zdz2,zdz3,lmel,entrbis,zdzbis |
90 |
|
|
real,dimension(ngrid) :: d_temp |
91 |
|
|
real ztv1,ztv2,factinv,zinv,zlmel |
92 |
|
|
real zlmelup,zlmeldwn,zlt,zltdwn,zltup |
93 |
|
|
real atv1,atv2,btv1,btv2 |
94 |
|
|
real ztv_est1,ztv_est2 |
95 |
|
|
real zcor,zdelta,zcvm5,qlbef |
96 |
|
|
real zbetalpha, coefzlmel |
97 |
|
|
real eps |
98 |
|
|
logical Zsat |
99 |
|
|
LOGICAL,dimension(ngrid) :: active,activetmp |
100 |
|
|
REAL fact_gamma,fact_gamma2,fact_epsilon2 |
101 |
|
|
|
102 |
|
|
|
103 |
|
|
REAL,dimension(ngrid,nlay) :: c2 |
104 |
|
|
|
105 |
|
|
if (ngrid==1) print*,'THERMCELL PLUME MODIFIE 2014/07/11' |
106 |
|
|
Zsat=.false. |
107 |
|
|
! Initialisation |
108 |
|
|
|
109 |
|
|
|
110 |
|
|
zbetalpha=betalpha/(1.+betalpha) |
111 |
|
|
|
112 |
|
|
|
113 |
|
|
! Initialisations des variables r?elles |
114 |
|
|
if (1==1) then |
115 |
|
|
ztva(:,:)=ztv(:,:) |
116 |
|
|
ztva_est(:,:)=ztva(:,:) |
117 |
|
|
ztv_est(:,:)=ztv(:,:) |
118 |
|
|
ztla(:,:)=zthl(:,:) |
119 |
|
|
zqta(:,:)=po(:,:) |
120 |
|
|
zqla(:,:)=0. |
121 |
|
|
zha(:,:) = ztva(:,:) |
122 |
|
|
else |
123 |
|
|
ztva(:,:)=0. |
124 |
|
|
ztv_est(:,:)=0. |
125 |
|
|
ztva_est(:,:)=0. |
126 |
|
|
ztla(:,:)=0. |
127 |
|
|
zqta(:,:)=0. |
128 |
|
|
zha(:,:) =0. |
129 |
|
|
endif |
130 |
|
|
|
131 |
|
|
zqla_est(:,:)=0. |
132 |
|
|
zqsatth(:,:)=0. |
133 |
|
|
zqla(:,:)=0. |
134 |
|
|
detr_star(:,:)=0. |
135 |
|
|
entr_star(:,:)=0. |
136 |
|
|
alim_star(:,:)=0. |
137 |
|
|
alim_star_tot(:)=0. |
138 |
|
|
csc(:,:)=0. |
139 |
|
|
detr(:,:)=0. |
140 |
|
|
entr(:,:)=0. |
141 |
|
|
zw2(:,:)=0. |
142 |
|
|
zbuoy(:,:)=0. |
143 |
|
|
zbuoyjam(:,:)=0. |
144 |
|
|
gamma(:,:)=0. |
145 |
|
|
zeps(:,:)=0. |
146 |
|
|
w_est(:,:)=0. |
147 |
|
|
f_star(:,:)=0. |
148 |
|
|
wa_moy(:,:)=0. |
149 |
|
|
linter(:)=1. |
150 |
|
|
! linter(:)=1. |
151 |
|
|
! Initialisation des variables entieres |
152 |
|
|
lmix(:)=1 |
153 |
|
|
lmix_bis(:)=2 |
154 |
|
|
wmaxa(:)=0. |
155 |
|
|
|
156 |
|
|
|
157 |
|
|
!------------------------------------------------------------------------- |
158 |
|
|
! On ne considere comme actif que les colonnes dont les deux premieres |
159 |
|
|
! couches sont instables. |
160 |
|
|
!------------------------------------------------------------------------- |
161 |
|
|
|
162 |
|
|
active(:)=ztv(:,1)>ztv(:,2) |
163 |
|
|
d_temp(:)=0. ! Pour activer un contraste de temperature a la base |
164 |
|
|
! du panache |
165 |
|
|
! Cet appel pourrait être fait avant thermcell_plume dans thermcell_main |
166 |
|
|
CALL thermcell_alim(thermals_flag_alim,ngrid,nlay,ztv,d_temp,zlev,alim_star,lalim) |
167 |
|
|
|
168 |
|
|
!------------------------------------------------------------------------------ |
169 |
|
|
! Calcul dans la premiere couche |
170 |
|
|
! On decide dans cette version que le thermique n'est actif que si la premiere |
171 |
|
|
! couche est instable. |
172 |
|
|
! Pourrait etre change si on veut que le thermiques puisse se d??clencher |
173 |
|
|
! dans une couche l>1 |
174 |
|
|
!------------------------------------------------------------------------------ |
175 |
|
|
do ig=1,ngrid |
176 |
|
|
! Le panache va prendre au debut les caracteristiques de l'air contenu |
177 |
|
|
! dans cette couche. |
178 |
|
|
if (active(ig)) then |
179 |
|
|
ztla(ig,1)=zthl(ig,1) |
180 |
|
|
zqta(ig,1)=po(ig,1) |
181 |
|
|
zqla(ig,1)=zl(ig,1) |
182 |
|
|
!cr: attention, prise en compte de f*(1)=1 |
183 |
|
|
f_star(ig,2)=alim_star(ig,1) |
184 |
|
|
zw2(ig,2)=2.*RG*(ztv(ig,1)-ztv(ig,2))/ztv(ig,2) & |
185 |
|
|
& *(zlev(ig,2)-zlev(ig,1)) & |
186 |
|
|
& *0.4*pphi(ig,1)/(pphi(ig,2)-pphi(ig,1)) |
187 |
|
|
w_est(ig,2)=zw2(ig,2) |
188 |
|
|
endif |
189 |
|
|
enddo |
190 |
|
|
! |
191 |
|
|
|
192 |
|
|
!============================================================================== |
193 |
|
|
!boucle de calcul de la vitesse verticale dans le thermique |
194 |
|
|
!============================================================================== |
195 |
|
|
do l=2,nlay-1 |
196 |
|
|
!============================================================================== |
197 |
|
|
|
198 |
|
|
|
199 |
|
|
! On decide si le thermique est encore actif ou non |
200 |
|
|
! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test |
201 |
|
|
do ig=1,ngrid |
202 |
|
|
active(ig)=active(ig) & |
203 |
|
|
& .and. zw2(ig,l)>1.e-10 & |
204 |
|
|
& .and. f_star(ig,l)+alim_star(ig,l)>1.e-10 |
205 |
|
|
enddo |
206 |
|
|
|
207 |
|
|
|
208 |
|
|
|
209 |
|
|
!--------------------------------------------------------------------------- |
210 |
|
|
! calcul des proprietes thermodynamiques et de la vitesse de la couche l |
211 |
|
|
! sans tenir compte du detrainement et de l'entrainement dans cette |
212 |
|
|
! couche |
213 |
|
|
! C'est a dire qu'on suppose |
214 |
|
|
! ztla(l)=ztla(l-1) et zqta(l)=zqta(l-1) |
215 |
|
|
! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer |
216 |
|
|
! avant) a l'alimentation pour avoir un calcul plus propre |
217 |
|
|
!--------------------------------------------------------------------------- |
218 |
|
|
|
219 |
|
|
ztemp(:)=zpspsk(:,l)*ztla(:,l-1) |
220 |
|
|
call thermcell_qsat(ngrid,active,pplev(:,l),ztemp,zqta(:,l-1),zqsat(:)) |
221 |
|
|
do ig=1,ngrid |
222 |
|
|
! print*,'active',active(ig),ig,l |
223 |
|
|
if(active(ig)) then |
224 |
|
|
zqla_est(ig,l)=max(0.,zqta(ig,l-1)-zqsat(ig)) |
225 |
|
|
ztva_est(ig,l) = ztla(ig,l-1)*zpspsk(ig,l)+RLvCp*zqla_est(ig,l) |
226 |
|
|
zta_est(ig,l)=ztva_est(ig,l) |
227 |
|
|
ztva_est(ig,l) = ztva_est(ig,l)/zpspsk(ig,l) |
228 |
|
|
ztva_est(ig,l) = ztva_est(ig,l)*(1.+RETV*(zqta(ig,l-1) & |
229 |
|
|
& -zqla_est(ig,l))-zqla_est(ig,l)) |
230 |
|
|
|
231 |
|
|
|
232 |
|
|
!Modif AJAM |
233 |
|
|
|
234 |
|
|
zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l) |
235 |
|
|
zdz=zlev(ig,l+1)-zlev(ig,l) |
236 |
|
|
lmel=fact_thermals_ed_dz*zlev(ig,l) |
237 |
|
|
! lmel=0.09*zlev(ig,l) |
238 |
|
|
zlmel=zlev(ig,l)+lmel |
239 |
|
|
zlmelup=zlmel+(zdz/2) |
240 |
|
|
zlmeldwn=zlmel-(zdz/2) |
241 |
|
|
|
242 |
|
|
lt=l+1 |
243 |
|
|
zlt=zlev(ig,lt) |
244 |
|
|
zdz3=zlev(ig,lt+1)-zlt |
245 |
|
|
zltdwn=zlt-zdz3/2 |
246 |
|
|
zltup=zlt+zdz3/2 |
247 |
|
|
|
248 |
|
|
!========================================================================= |
249 |
|
|
! 3. Calcul de la flotabilite modifie par melange avec l'air au dessus |
250 |
|
|
!========================================================================= |
251 |
|
|
|
252 |
|
|
!-------------------------------------------------- |
253 |
|
|
lt=l+1 |
254 |
|
|
zlt=zlev(ig,lt) |
255 |
|
|
zdz2=zlev(ig,lt)-zlev(ig,l) |
256 |
|
|
|
257 |
|
|
do while (lmel.gt.zdz2) |
258 |
|
|
lt=lt+1 |
259 |
|
|
zlt=zlev(ig,lt) |
260 |
|
|
zdz2=zlt-zlev(ig,l) |
261 |
|
|
enddo |
262 |
|
|
zdz3=zlev(ig,lt+1)-zlt |
263 |
|
|
zltdwn=zlev(ig,lt)-zdz3/2 |
264 |
|
|
zlmelup=zlmel+(zdz/2) |
265 |
|
|
coefzlmel=Min(1.,(zlmelup-zltdwn)/zdz) |
266 |
|
|
zbuoyjam(ig,l)=1.*RG*(coefzlmel*(ztva_est(ig,l)- & |
267 |
|
|
& ztv(ig,lt))/ztv(ig,lt)+(1.-coefzlmel)*(ztva_est(ig,l)- & |
268 |
|
|
& ztv(ig,lt-1))/ztv(ig,lt-1))+0.*zbuoy(ig,l) |
269 |
|
|
|
270 |
|
|
!------------------------------------------------ |
271 |
|
|
!AJAM:nouveau calcul de w? |
272 |
|
|
!------------------------------------------------ |
273 |
|
|
zdz=zlev(ig,l+1)-zlev(ig,l) |
274 |
|
|
zdzbis=zlev(ig,l)-zlev(ig,l-1) |
275 |
|
|
zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l) |
276 |
|
|
zw2fact=fact_epsilon*2.*zdz/(1.+betalpha) |
277 |
|
|
zw2factbis=fact_epsilon*2.*zdzbis/(1.+betalpha) |
278 |
|
|
zdw2=afact*zbuoy(ig,l)/fact_epsilon |
279 |
|
|
zdw2bis=afact*zbuoy(ig,l-1)/fact_epsilon |
280 |
|
|
lm=Max(1,l-2) |
281 |
|
|
w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2)+zdw2) |
282 |
|
|
endif |
283 |
|
|
enddo |
284 |
|
|
|
285 |
|
|
|
286 |
|
|
!------------------------------------------------- |
287 |
|
|
!calcul des taux d'entrainement et de detrainement |
288 |
|
|
!------------------------------------------------- |
289 |
|
|
|
290 |
|
|
do ig=1,ngrid |
291 |
|
|
if (active(ig)) then |
292 |
|
|
|
293 |
|
|
! zw2m=max(0.5*(w_est(ig,l)+w_est(ig,l+1)),0.1) |
294 |
|
|
zw2m=w_est(ig,l+1) |
295 |
|
|
zdz=zlev(ig,l+1)-zlev(ig,l) |
296 |
|
|
zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l) |
297 |
|
|
zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l+1))/rhobarz(ig,l) |
298 |
|
|
zdqt(ig,l)=max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l) |
299 |
|
|
|
300 |
|
|
!========================================================================= |
301 |
|
|
! 4. Calcul de l'entrainement et du detrainement |
302 |
|
|
!========================================================================= |
303 |
|
|
|
304 |
|
|
detr_star(ig,l)=f_star(ig,l)*zdz & |
305 |
|
|
& *( mix0 * 0.1 / (zalpha+0.001) & |
306 |
|
|
& + MAX(detr_min, -afact*zbetalpha*zbuoyjam(ig,l)/zw2m & |
307 |
|
|
& + detr_q_coef*(zdqt(ig,l)/zw2m)**detr_q_power)) |
308 |
|
|
|
309 |
|
|
if ( iflag_thermals_ed == 20 ) then |
310 |
|
|
entr_star(ig,l)=f_star(ig,l)*zdz* ( & |
311 |
|
|
& mix0 * 0.1 / (zalpha+0.001) & |
312 |
|
|
& + zbetalpha*MAX(entr_min, & |
313 |
|
|
& afact*zbuoyjam(ig,l)/zw2m - fact_epsilon)) |
314 |
|
|
else |
315 |
|
|
entr_star(ig,l)=f_star(ig,l)*zdz* ( & |
316 |
|
|
& mix0 * 0.1 / (zalpha+0.001) & |
317 |
|
|
& + zbetalpha*MAX(entr_min, & |
318 |
|
|
& afact*zbuoy(ig,l)/zw2m - fact_epsilon)) |
319 |
|
|
endif |
320 |
|
|
|
321 |
|
|
! En dessous de lalim, on prend le max de alim_star et entr_star pour |
322 |
|
|
! alim_star et 0 sinon |
323 |
|
|
if (l.lt.lalim(ig)) then |
324 |
|
|
alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l)) |
325 |
|
|
entr_star(ig,l)=0. |
326 |
|
|
endif |
327 |
|
|
f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l) & |
328 |
|
|
& -detr_star(ig,l) |
329 |
|
|
|
330 |
|
|
endif |
331 |
|
|
enddo |
332 |
|
|
|
333 |
|
|
|
334 |
|
|
!============================================================================ |
335 |
|
|
! 5. calcul de la vitesse verticale en melangeant Tl et qt du thermique |
336 |
|
|
!=========================================================================== |
337 |
|
|
|
338 |
|
|
activetmp(:)=active(:) .and. f_star(:,l+1)>1.e-10 |
339 |
|
|
do ig=1,ngrid |
340 |
|
|
if (activetmp(ig)) then |
341 |
|
|
Zsat=.false. |
342 |
|
|
ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+ & |
343 |
|
|
& (alim_star(ig,l)+entr_star(ig,l))*zthl(ig,l)) & |
344 |
|
|
& /(f_star(ig,l+1)+detr_star(ig,l)) |
345 |
|
|
zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+ & |
346 |
|
|
& (alim_star(ig,l)+entr_star(ig,l))*po(ig,l)) & |
347 |
|
|
& /(f_star(ig,l+1)+detr_star(ig,l)) |
348 |
|
|
|
349 |
|
|
endif |
350 |
|
|
enddo |
351 |
|
|
|
352 |
|
|
ztemp(:)=zpspsk(:,l)*ztla(:,l) |
353 |
|
|
call thermcell_qsat(ngrid,activetmp,pplev(:,l),ztemp,zqta(:,l),zqsatth(:,l)) |
354 |
|
|
do ig=1,ngrid |
355 |
|
|
if (activetmp(ig)) then |
356 |
|
|
! on ecrit de maniere conservative (sat ou non) |
357 |
|
|
! T = Tl +Lv/Cp ql |
358 |
|
|
zqla(ig,l)=max(0.,zqta(ig,l)-zqsatth(ig,l)) |
359 |
|
|
ztva(ig,l) = ztla(ig,l)*zpspsk(ig,l)+RLvCp*zqla(ig,l) |
360 |
|
|
ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l) |
361 |
|
|
!on rajoute le calcul de zha pour diagnostiques (temp potentielle) |
362 |
|
|
zha(ig,l) = ztva(ig,l) |
363 |
|
|
ztva(ig,l) = ztva(ig,l)*(1.+RETV*(zqta(ig,l) & |
364 |
|
|
& -zqla(ig,l))-zqla(ig,l)) |
365 |
|
|
zbuoy(ig,l)=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l) |
366 |
|
|
zdz=zlev(ig,l+1)-zlev(ig,l) |
367 |
|
|
zdzbis=zlev(ig,l)-zlev(ig,l-1) |
368 |
|
|
zeps(ig,l)=(entr_star(ig,l)+alim_star(ig,l))/(f_star(ig,l)*zdz) |
369 |
|
|
zw2fact=fact_epsilon*2.*zdz/(1.+betalpha) |
370 |
|
|
zw2factbis=fact_epsilon*2.*zdzbis/(1.+betalpha) |
371 |
|
|
zdw2= afact*zbuoy(ig,l)/(fact_epsilon) |
372 |
|
|
zdw2bis= afact*zbuoy(ig,l-1)/(fact_epsilon) |
373 |
|
|
zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2)+zdw2) |
374 |
|
|
endif |
375 |
|
|
enddo |
376 |
|
|
|
377 |
|
|
if (prt_level.ge.20) print*,'coucou calcul detr 460: ig, l',ig, l |
378 |
|
|
! |
379 |
|
|
!=========================================================================== |
380 |
|
|
! 6. initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max |
381 |
|
|
!=========================================================================== |
382 |
|
|
|
383 |
|
|
nbpb=0 |
384 |
|
|
do ig=1,ngrid |
385 |
|
|
if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then |
386 |
|
|
! stop'On tombe sur le cas particulier de thermcell_dry' |
387 |
|
|
! print*,'On tombe sur le cas particulier de thermcell_plume' |
388 |
|
|
nbpb=nbpb+1 |
389 |
|
|
zw2(ig,l+1)=0. |
390 |
|
|
linter(ig)=l+1 |
391 |
|
|
endif |
392 |
|
|
|
393 |
|
|
if (zw2(ig,l+1).lt.0.) then |
394 |
|
|
linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l)) & |
395 |
|
|
& -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l)) |
396 |
|
|
zw2(ig,l+1)=0. |
397 |
|
|
!+CR:04/05/12:correction calcul linter pour calcul de zmax continu |
398 |
|
|
elseif (f_star(ig,l+1).lt.0.) then |
399 |
|
|
linter(ig)=(l*(f_star(ig,l+1)-f_star(ig,l)) & |
400 |
|
|
& -f_star(ig,l))/(f_star(ig,l+1)-f_star(ig,l)) |
401 |
|
|
zw2(ig,l+1)=0. |
402 |
|
|
!fin CR:04/05/12 |
403 |
|
|
endif |
404 |
|
|
|
405 |
|
|
wa_moy(ig,l+1)=sqrt(zw2(ig,l+1)) |
406 |
|
|
|
407 |
|
|
if (wa_moy(ig,l+1).gt.wmaxa(ig)) then |
408 |
|
|
! lmix est le niveau de la couche ou w (wa_moy) est maximum |
409 |
|
|
!on rajoute le calcul de lmix_bis |
410 |
|
|
if (zqla(ig,l).lt.1.e-10) then |
411 |
|
|
lmix_bis(ig)=l+1 |
412 |
|
|
endif |
413 |
|
|
lmix(ig)=l+1 |
414 |
|
|
wmaxa(ig)=wa_moy(ig,l+1) |
415 |
|
|
endif |
416 |
|
|
enddo |
417 |
|
|
|
418 |
|
|
if (nbpb>0) then |
419 |
|
|
print*,'WARNING on tombe ',nbpb,' x sur un pb pour l=',l,' dans thermcell_plume' |
420 |
|
|
endif |
421 |
|
|
|
422 |
|
|
!========================================================================= |
423 |
|
|
! FIN DE LA BOUCLE VERTICALE |
424 |
|
|
enddo |
425 |
|
|
!========================================================================= |
426 |
|
|
|
427 |
|
|
!on recalcule alim_star_tot |
428 |
|
|
do ig=1,ngrid |
429 |
|
|
alim_star_tot(ig)=0. |
430 |
|
|
enddo |
431 |
|
|
do ig=1,ngrid |
432 |
|
|
do l=1,lalim(ig)-1 |
433 |
|
|
alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l) |
434 |
|
|
enddo |
435 |
|
|
enddo |
436 |
|
|
|
437 |
|
|
|
438 |
|
|
if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l |
439 |
|
|
|
440 |
|
|
#undef wrgrads_thermcell |
441 |
|
|
#ifdef wrgrads_thermcell |
442 |
|
|
call wrgradsfi(1,nlay,entr_star(igout,1:nlay),'esta ','esta ') |
443 |
|
|
call wrgradsfi(1,nlay,detr_star(igout,1:nlay),'dsta ','dsta ') |
444 |
|
|
call wrgradsfi(1,nlay,zbuoy(igout,1:nlay),'buoy ','buoy ') |
445 |
|
|
call wrgradsfi(1,nlay,zdqt(igout,1:nlay),'dqt ','dqt ') |
446 |
|
|
call wrgradsfi(1,nlay,w_est(igout,1:nlay),'w_est ','w_est ') |
447 |
|
|
call wrgradsfi(1,nlay,w_est(igout,2:nlay+1),'w_es2 ','w_es2 ') |
448 |
|
|
call wrgradsfi(1,nlay,zw2(igout,1:nlay),'zw2A ','zw2A ') |
449 |
|
|
#endif |
450 |
|
|
|
451 |
|
|
|
452 |
|
|
RETURN |
453 |
|
|
end |
454 |
|
|
END MODULE lmdz_thermcell_plume |