1 |
|
|
MODULE lmdz_thermcell_plume_6A |
2 |
|
|
! |
3 |
|
|
! $Id: lmdz_thermcell_plume_6A.F90 4590 2023-06-29 01:03:15Z fhourdin $ |
4 |
|
|
! |
5 |
|
|
CONTAINS |
6 |
|
|
|
7 |
|
288 |
SUBROUTINE thermcell_plume_6A(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz, & |
8 |
|
288 |
& zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot, & |
9 |
|
288 |
& lalim,f0,detr_star,entr_star,f_star,csc,ztva, & |
10 |
|
288 |
& 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 |
|
|
!thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance |
15 |
|
|
!-------------------------------------------------------------------------- |
16 |
|
|
|
17 |
|
|
USE lmdz_thermcell_ini, ONLY: prt_level,fact_thermals_ed_dz,iflag_thermals_ed,RLvCP,RETV,RG |
18 |
|
|
USE lmdz_thermcell_ini, ONLY: fact_epsilon, betalpha, afact, fact_shell |
19 |
|
|
USE lmdz_thermcell_ini, ONLY: detr_min, entr_min, detr_q_coef, detr_q_power |
20 |
|
|
USE lmdz_thermcell_ini, ONLY: mix0, thermals_flag_alim |
21 |
|
|
USE lmdz_thermcell_alim, ONLY : thermcell_alim |
22 |
|
|
USE lmdz_thermcell_qsat, ONLY : thermcell_qsat |
23 |
|
|
|
24 |
|
|
|
25 |
|
|
IMPLICIT NONE |
26 |
|
|
|
27 |
|
|
integer,intent(in) :: itap,lev_out,lunout1,igout,ngrid,nlay |
28 |
|
|
real,intent(in) :: ptimestep |
29 |
|
|
real,intent(in),dimension(ngrid,nlay) :: ztv |
30 |
|
|
real,intent(in),dimension(ngrid,nlay) :: zthl |
31 |
|
|
real,intent(in),dimension(ngrid,nlay) :: po |
32 |
|
|
real,intent(in),dimension(ngrid,nlay) :: zl |
33 |
|
|
real,intent(in),dimension(ngrid,nlay) :: rhobarz |
34 |
|
|
real,intent(in),dimension(ngrid,nlay+1) :: zlev |
35 |
|
|
real,intent(in),dimension(ngrid,nlay+1) :: pplev |
36 |
|
|
real,intent(in),dimension(ngrid,nlay) :: pphi |
37 |
|
|
real,intent(in),dimension(ngrid,nlay) :: zpspsk |
38 |
|
|
real,intent(in),dimension(ngrid) :: f0 |
39 |
|
|
|
40 |
|
|
integer,intent(out) :: lalim(ngrid) |
41 |
|
|
real,intent(out),dimension(ngrid,nlay) :: alim_star |
42 |
|
|
real,intent(out),dimension(ngrid) :: alim_star_tot |
43 |
|
|
real,intent(out),dimension(ngrid,nlay) :: detr_star |
44 |
|
|
real,intent(out),dimension(ngrid,nlay) :: entr_star |
45 |
|
|
real,intent(out),dimension(ngrid,nlay+1) :: f_star |
46 |
|
|
real,intent(out),dimension(ngrid,nlay) :: csc |
47 |
|
|
real,intent(out),dimension(ngrid,nlay) :: ztva |
48 |
|
|
real,intent(out),dimension(ngrid,nlay) :: ztla |
49 |
|
|
real,intent(out),dimension(ngrid,nlay) :: zqla |
50 |
|
|
real,intent(out),dimension(ngrid,nlay) :: zqta |
51 |
|
|
real,intent(out),dimension(ngrid,nlay) :: zha |
52 |
|
|
real,intent(out),dimension(ngrid,nlay+1) :: zw2 |
53 |
|
|
real,intent(out),dimension(ngrid,nlay+1) :: w_est |
54 |
|
|
real,intent(out),dimension(ngrid,nlay) :: ztva_est |
55 |
|
|
real,intent(out),dimension(ngrid,nlay) :: zqsatth |
56 |
|
|
integer,intent(out),dimension(ngrid) :: lmix |
57 |
|
|
integer,intent(out),dimension(ngrid) :: lmix_bis |
58 |
|
|
real,intent(out),dimension(ngrid) :: linter |
59 |
|
|
|
60 |
|
|
REAL zdw2,zdw2bis |
61 |
|
|
REAL zw2modif |
62 |
|
|
REAL zw2fact,zw2factbis |
63 |
|
576 |
REAL,dimension(ngrid,nlay) :: zeps |
64 |
|
|
|
65 |
|
576 |
REAL, dimension(ngrid) :: wmaxa(ngrid) |
66 |
|
|
|
67 |
|
|
INTEGER ig,l,k,lt,it,lm |
68 |
|
|
integer nbpb |
69 |
|
|
|
70 |
|
576 |
real,dimension(ngrid,nlay) :: detr |
71 |
|
576 |
real,dimension(ngrid,nlay) :: entr |
72 |
|
576 |
real,dimension(ngrid,nlay+1) :: wa_moy |
73 |
|
576 |
real,dimension(ngrid,nlay) :: ztv_est |
74 |
|
576 |
real,dimension(ngrid) :: ztemp,zqsat |
75 |
|
576 |
real,dimension(ngrid,nlay) :: zqla_est |
76 |
|
576 |
real,dimension(ngrid,nlay) :: zta_est |
77 |
|
|
|
78 |
|
576 |
real,dimension(ngrid,nlay) :: zbuoy,gamma,zdqt |
79 |
|
|
real zdz,zalpha,zw2m |
80 |
|
576 |
real,dimension(ngrid,nlay) :: zbuoyjam,zdqtjam |
81 |
|
|
real zbuoybis,zdz2,zdz3,lmel,entrbis,zdzbis |
82 |
|
576 |
real, dimension(ngrid) :: d_temp |
83 |
|
|
real ztv1,ztv2,factinv,zinv,zlmel |
84 |
|
|
real zlmelup,zlmeldwn,zlt,zltdwn,zltup |
85 |
|
|
real atv1,atv2,btv1,btv2 |
86 |
|
|
real ztv_est1,ztv_est2 |
87 |
|
|
real zcor,zdelta,zcvm5,qlbef |
88 |
|
|
real zbetalpha, coefzlmel |
89 |
|
|
real eps |
90 |
|
|
logical Zsat |
91 |
|
576 |
LOGICAL,dimension(ngrid) :: active,activetmp |
92 |
|
|
REAL fact_gamma,fact_gamma2,fact_epsilon2 |
93 |
|
|
REAL coefc |
94 |
|
576 |
REAL,dimension(ngrid,nlay) :: c2 |
95 |
|
|
|
96 |
✗✓ |
288 |
if (ngrid==1) print*,'THERMCELL PLUME MODIFIE 2014/07/11' |
97 |
|
|
Zsat=.false. |
98 |
|
|
! Initialisation |
99 |
|
|
|
100 |
|
|
|
101 |
|
288 |
zbetalpha=betalpha/(1.+betalpha) |
102 |
|
|
|
103 |
|
|
|
104 |
|
|
! Initialisations des variables r?elles |
105 |
|
|
if (1==1) then |
106 |
✓✓✓✓
|
11176128 |
ztva(:,:)=ztv(:,:) |
107 |
✓✓✓✓
|
11176128 |
ztva_est(:,:)=ztva(:,:) |
108 |
✓✓✓✓
|
11176128 |
ztv_est(:,:)=ztv(:,:) |
109 |
✓✓✓✓
|
11176128 |
ztla(:,:)=zthl(:,:) |
110 |
✓✓✓✓
|
11176128 |
zqta(:,:)=po(:,:) |
111 |
✓✓✓✓
|
11176128 |
zqla(:,:)=0. |
112 |
✓✓✓✓
|
11176128 |
zha(:,:) = ztva(:,:) |
113 |
|
|
else |
114 |
|
|
ztva(:,:)=0. |
115 |
|
|
ztv_est(:,:)=0. |
116 |
|
|
ztva_est(:,:)=0. |
117 |
|
|
ztla(:,:)=0. |
118 |
|
|
zqta(:,:)=0. |
119 |
|
|
zha(:,:) =0. |
120 |
|
|
endif |
121 |
|
|
|
122 |
✓✓✓✓
|
11176128 |
zqla_est(:,:)=0. |
123 |
✓✓✓✓
|
11176128 |
zqsatth(:,:)=0. |
124 |
✓✓✓✓
|
11176128 |
zqla(:,:)=0. |
125 |
✓✓✓✓
|
11176128 |
detr_star(:,:)=0. |
126 |
✓✓✓✓
|
11176128 |
entr_star(:,:)=0. |
127 |
✓✓✓✓
|
11176128 |
alim_star(:,:)=0. |
128 |
✓✓ |
286560 |
alim_star_tot(:)=0. |
129 |
✓✓✓✓
|
11176128 |
csc(:,:)=0. |
130 |
✓✓✓✓
|
11176128 |
detr(:,:)=0. |
131 |
✓✓✓✓
|
11176128 |
entr(:,:)=0. |
132 |
✓✓✓✓
|
11462688 |
zw2(:,:)=0. |
133 |
✓✓✓✓
|
11176128 |
zbuoy(:,:)=0. |
134 |
✓✓✓✓
|
11176128 |
zbuoyjam(:,:)=0. |
135 |
✓✓✓✓
|
11176128 |
gamma(:,:)=0. |
136 |
✓✓✓✓
|
11176128 |
zeps(:,:)=0. |
137 |
✓✓✓✓
|
11462688 |
w_est(:,:)=0. |
138 |
✓✓✓✓
|
11462688 |
f_star(:,:)=0. |
139 |
✓✓✓✓
|
11462688 |
wa_moy(:,:)=0. |
140 |
✓✓ |
286560 |
linter(:)=1. |
141 |
|
|
! linter(:)=1. |
142 |
|
|
! Initialisation des variables entieres |
143 |
✓✓ |
286560 |
lmix(:)=1 |
144 |
✓✓ |
286560 |
lmix_bis(:)=2 |
145 |
✓✓ |
286560 |
wmaxa(:)=0. |
146 |
|
|
|
147 |
|
|
! Initialisation a 0 en cas de sortie dans replay |
148 |
✓✓ |
286560 |
zqsat(:)=0. |
149 |
✓✓✓✓
|
11176128 |
zta_est(:,:)=0. |
150 |
✓✓✓✓
|
11176128 |
zdqt(:,:)=0. |
151 |
✓✓✓✓
|
11176128 |
zdqtjam(:,:)=0. |
152 |
✓✓✓✓
|
11176128 |
c2(:,:)=0. |
153 |
|
|
|
154 |
|
|
|
155 |
|
|
!------------------------------------------------------------------------- |
156 |
|
|
! On ne considere comme actif que les colonnes dont les deux premieres |
157 |
|
|
! couches sont instables. |
158 |
|
|
!------------------------------------------------------------------------- |
159 |
|
|
|
160 |
✓✓ |
286560 |
active(:)=ztv(:,1)>ztv(:,2) |
161 |
✓✓ |
286560 |
d_temp(:)=0. ! Pour activer un contraste de temperature a la base |
162 |
|
|
! du panache |
163 |
|
|
! Cet appel pourrait être fait avant thermcell_plume dans thermcell_main |
164 |
|
288 |
CALL thermcell_alim(thermals_flag_alim,ngrid,nlay,ztv,d_temp,zlev,alim_star,lalim) |
165 |
|
|
|
166 |
|
|
!------------------------------------------------------------------------------ |
167 |
|
|
! Calcul dans la premiere couche |
168 |
|
|
! On decide dans cette version que le thermique n'est actif que si la premiere |
169 |
|
|
! couche est instable. |
170 |
|
|
! Pourrait etre change si on veut que le thermiques puisse se d??clencher |
171 |
|
|
! dans une couche l>1 |
172 |
|
|
!------------------------------------------------------------------------------ |
173 |
✓✓ |
286560 |
do ig=1,ngrid |
174 |
|
|
! Le panache va prendre au debut les caracteristiques de l'air contenu |
175 |
|
|
! dans cette couche. |
176 |
✓✓ |
286560 |
if (active(ig)) then |
177 |
|
137337 |
ztla(ig,1)=zthl(ig,1) |
178 |
|
137337 |
zqta(ig,1)=po(ig,1) |
179 |
|
137337 |
zqla(ig,1)=zl(ig,1) |
180 |
|
|
!cr: attention, prise en compte de f*(1)=1 |
181 |
|
137337 |
f_star(ig,2)=alim_star(ig,1) |
182 |
|
|
zw2(ig,2)=2.*RG*(ztv(ig,1)-ztv(ig,2))/ztv(ig,2) & |
183 |
|
|
& *(zlev(ig,2)-zlev(ig,1)) & |
184 |
|
137337 |
& *0.4*pphi(ig,1)/(pphi(ig,2)-pphi(ig,1)) |
185 |
|
137337 |
w_est(ig,2)=zw2(ig,2) |
186 |
|
|
endif |
187 |
|
|
enddo |
188 |
|
|
! |
189 |
|
|
|
190 |
|
|
!============================================================================== |
191 |
|
|
!boucle de calcul de la vitesse verticale dans le thermique |
192 |
|
|
!============================================================================== |
193 |
✓✓ |
10944 |
do l=2,nlay-1 |
194 |
|
|
!============================================================================== |
195 |
|
|
|
196 |
|
|
|
197 |
|
|
! On decide si le thermique est encore actif ou non |
198 |
|
|
! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test |
199 |
✓✓ |
10602720 |
do ig=1,ngrid |
200 |
|
|
active(ig)=active(ig) & |
201 |
|
|
& .and. zw2(ig,l)>1.e-10 & |
202 |
✓✓✓✓
✗✓ |
20453429 |
& .and. f_star(ig,l)+alim_star(ig,l)>1.e-10 |
203 |
|
|
enddo |
204 |
|
|
|
205 |
|
|
|
206 |
|
|
|
207 |
|
|
!--------------------------------------------------------------------------- |
208 |
|
|
! calcul des proprietes thermodynamiques et de la vitesse de la couche l |
209 |
|
|
! sans tenir compte du detrainement et de l'entrainement dans cette |
210 |
|
|
! couche |
211 |
|
|
! C'est a dire qu'on suppose |
212 |
|
|
! ztla(l)=ztla(l-1) et zqta(l)=zqta(l-1) |
213 |
|
|
! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer |
214 |
|
|
! avant) a l'alimentation pour avoir un calcul plus propre |
215 |
|
|
!--------------------------------------------------------------------------- |
216 |
|
|
|
217 |
✓✓ |
10602720 |
ztemp(:)=zpspsk(:,l)*ztla(:,l-1) |
218 |
|
10656 |
call thermcell_qsat(ngrid,active,pplev(:,l),ztemp,zqta(:,l-1),zqsat(:)) |
219 |
✓✓ |
10602720 |
do ig=1,ngrid |
220 |
|
|
! print*,'active',active(ig),ig,l |
221 |
✓✓ |
10602720 |
if(active(ig)) then |
222 |
|
741355 |
zqla_est(ig,l)=max(0.,zqta(ig,l-1)-zqsat(ig)) |
223 |
|
741355 |
ztva_est(ig,l) = ztla(ig,l-1)*zpspsk(ig,l)+RLvCp*zqla_est(ig,l) |
224 |
|
741355 |
zta_est(ig,l)=ztva_est(ig,l) |
225 |
|
741355 |
ztva_est(ig,l) = ztva_est(ig,l)/zpspsk(ig,l) |
226 |
|
|
ztva_est(ig,l) = ztva_est(ig,l)*(1.+RETV*(zqta(ig,l-1) & |
227 |
|
741355 |
& -zqla_est(ig,l))-zqla_est(ig,l)) |
228 |
|
|
|
229 |
|
|
|
230 |
|
|
!Modif AJAM |
231 |
|
|
|
232 |
|
741355 |
zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l) |
233 |
|
741355 |
zdz=zlev(ig,l+1)-zlev(ig,l) |
234 |
|
741355 |
lmel=fact_thermals_ed_dz*zlev(ig,l) |
235 |
|
|
! lmel=0.09*zlev(ig,l) |
236 |
|
741355 |
zlmel=zlev(ig,l)+lmel |
237 |
|
741355 |
zlmelup=zlmel+(zdz/2) |
238 |
|
741355 |
zlmeldwn=zlmel-(zdz/2) |
239 |
|
|
|
240 |
|
|
lt=l+1 |
241 |
|
|
zlt=zlev(ig,lt) |
242 |
|
741355 |
zdz3=zlev(ig,lt+1)-zlt |
243 |
|
741355 |
zltdwn=zlt-zdz3/2 |
244 |
|
741355 |
zltup=zlt+zdz3/2 |
245 |
|
|
|
246 |
|
|
!========================================================================= |
247 |
|
|
! 3. Calcul de la flotabilite modifie par melange avec l'air au dessus |
248 |
|
|
!========================================================================= |
249 |
|
|
|
250 |
|
|
!-------------------------------------------------- |
251 |
✗✓ |
741355 |
if (iflag_thermals_ed.lt.8) then |
252 |
|
|
!-------------------------------------------------- |
253 |
|
|
!AJ052014: J'ai remplac?? la boucle do par un do while |
254 |
|
|
! afin de faire moins de calcul dans la boucle |
255 |
|
|
!-------------------------------------------------- |
256 |
|
|
do while (zlmelup.gt.zltup) |
257 |
|
|
lt=lt+1 |
258 |
|
|
zlt=zlev(ig,lt) |
259 |
|
|
zdz3=zlev(ig,lt+1)-zlt |
260 |
|
|
zltdwn=zlt-zdz3/2 |
261 |
|
|
zltup=zlt+zdz3/2 |
262 |
|
|
enddo |
263 |
|
|
!-------------------------------------------------- |
264 |
|
|
!AJ052014: Si iflag_thermals_ed<8 (par ex 6), alors |
265 |
|
|
! on cherche o?? se trouve l'altitude d'inversion |
266 |
|
|
! en calculant ztv1 (interpolation de la valeur de |
267 |
|
|
! theta au niveau lt en utilisant les niveaux lt-1 et |
268 |
|
|
! lt-2) et ztv2 (interpolation avec les niveaux lt+1 |
269 |
|
|
! et lt+2). Si theta r??ellement calcul??e au niveau lt |
270 |
|
|
! comprise entre ztv1 et ztv2, alors il y a inversion |
271 |
|
|
! et on calcule son altitude zinv en supposant que ztv(lt) |
272 |
|
|
! est une combinaison lineaire de ztv1 et ztv2. |
273 |
|
|
! Ensuite, on calcule la flottabilite en comparant |
274 |
|
|
! la temperature de la couche l a celle de l'air situe |
275 |
|
|
! l+lmel plus haut, ce qui necessite de savoir quel fraction |
276 |
|
|
! de cet air est au-dessus ou en-dessous de l'inversion |
277 |
|
|
!-------------------------------------------------- |
278 |
|
|
atv1=(ztv(ig,lt-1)-ztv(ig,lt-2))/(zlev(ig,lt-1)-zlev(ig,lt-2)) |
279 |
|
|
btv1=(ztv(ig,lt-2)*zlev(ig,lt-1)-ztv(ig,lt-1)*zlev(ig,lt-2)) & |
280 |
|
|
& /(zlev(ig,lt-1)-zlev(ig,lt-2)) |
281 |
|
|
atv2=(ztv(ig,lt+2)-ztv(ig,lt+1))/(zlev(ig,lt+2)-zlev(ig,lt+1)) |
282 |
|
|
btv2=(ztv(ig,lt+1)*zlev(ig,lt+2)-ztv(ig,lt+2)*zlev(ig,lt+1)) & |
283 |
|
|
& /(zlev(ig,lt+2)-zlev(ig,lt+1)) |
284 |
|
|
|
285 |
|
|
ztv1=atv1*zlt+btv1 |
286 |
|
|
ztv2=atv2*zlt+btv2 |
287 |
|
|
|
288 |
|
|
if (ztv(ig,lt).gt.ztv1.and.ztv(ig,lt).lt.ztv2) then |
289 |
|
|
|
290 |
|
|
!-------------------------------------------------- |
291 |
|
|
!AJ052014: D??calage de zinv qui est entre le haut |
292 |
|
|
! et le bas de la couche lt |
293 |
|
|
!-------------------------------------------------- |
294 |
|
|
factinv=(ztv2-ztv(ig,lt))/(ztv2-ztv1) |
295 |
|
|
zinv=zltdwn+zdz3*factinv |
296 |
|
|
|
297 |
|
|
|
298 |
|
|
if (zlmeldwn.ge.zinv) then |
299 |
|
|
ztv_est(ig,l)=atv2*zlmel+btv2 |
300 |
|
|
zbuoyjam(ig,l)=fact_shell*RG*(ztva_est(ig,l)-ztv_est(ig,l))/ztv_est(ig,l) & |
301 |
|
|
& +(1.-fact_shell)*zbuoy(ig,l) |
302 |
|
|
elseif (zlmelup.ge.zinv) then |
303 |
|
|
ztv_est2=atv2*0.5*(zlmelup+zinv)+btv2 |
304 |
|
|
ztv_est1=atv1*0.5*(zinv+zlmeldwn)+btv1 |
305 |
|
|
ztv_est(ig,l)=((zlmelup-zinv)/zdz)*ztv_est2+((zinv-zlmeldwn)/zdz)*ztv_est1 |
306 |
|
|
|
307 |
|
|
zbuoyjam(ig,l)=fact_shell*RG*(((zlmelup-zinv)/zdz)*(ztva_est(ig,l)- & |
308 |
|
|
& ztv_est2)/ztv_est2+((zinv-zlmeldwn)/zdz)*(ztva_est(ig,l)- & |
309 |
|
|
& ztv_est1)/ztv_est1)+(1.-fact_shell)*zbuoy(ig,l) |
310 |
|
|
|
311 |
|
|
else |
312 |
|
|
ztv_est(ig,l)=atv1*zlmel+btv1 |
313 |
|
|
zbuoyjam(ig,l)=fact_shell*RG*(ztva_est(ig,l)-ztv_est(ig,l))/ztv_est(ig,l) & |
314 |
|
|
& +(1.-fact_shell)*zbuoy(ig,l) |
315 |
|
|
endif |
316 |
|
|
|
317 |
|
|
else ! if (ztv(ig,lt).gt.ztv1.and.ztv(ig,lt).lt.ztv2) then |
318 |
|
|
|
319 |
|
|
if (zlmeldwn.gt.zltdwn) then |
320 |
|
|
zbuoyjam(ig,l)=fact_shell*RG*((ztva_est(ig,l)- & |
321 |
|
|
& ztv(ig,lt))/ztv(ig,lt))+(1.-fact_shell)*zbuoy(ig,l) |
322 |
|
|
else |
323 |
|
|
zbuoyjam(ig,l)=fact_shell*RG*(((zlmelup-zltdwn)/zdz)*(ztva_est(ig,l)- & |
324 |
|
|
& ztv(ig,lt))/ztv(ig,lt)+((zltdwn-zlmeldwn)/zdz)*(ztva_est(ig,l)- & |
325 |
|
|
& ztv(ig,lt-1))/ztv(ig,lt-1))+(1.-fact_shell)*zbuoy(ig,l) |
326 |
|
|
|
327 |
|
|
endif |
328 |
|
|
|
329 |
|
|
! zbuoyjam(ig,l)=fact_shell*RG*(((zlmelup-zltdwn)/zdz)*(ztva_est(ig,l)- & |
330 |
|
|
! & ztv1)/ztv1+((zltdwn-zlmeldwn)/zdz)*(ztva_est(ig,l)- & |
331 |
|
|
! & ztv(ig,lt-1))/ztv(ig,lt-1))+(1.-fact_shell)*zbuoy(ig,l) |
332 |
|
|
! zdqt(ig,l)=Max(0.,((lmel+zdz3-zdz2)/zdz3)*(zqta(ig,l-1)- & |
333 |
|
|
! & po(ig,lt))/po(ig,lt)+((zdz2-lmel)/zdz3)*(zqta(ig,l-1)- & |
334 |
|
|
! & po(ig,lt-1))/po(ig,lt-1)) |
335 |
|
|
endif ! if (ztv(ig,lt).gt.ztv1.and.ztv(ig,lt).lt.ztv2) then |
336 |
|
|
|
337 |
|
|
else ! if (iflag_thermals_ed.lt.8) then |
338 |
|
|
lt=l+1 |
339 |
|
|
zlt=zlev(ig,lt) |
340 |
|
|
zdz2=zlev(ig,lt)-zlev(ig,l) |
341 |
|
|
|
342 |
✗✓ |
741355 |
do while (lmel.gt.zdz2) |
343 |
|
|
lt=lt+1 |
344 |
|
|
zlt=zlev(ig,lt) |
345 |
|
|
zdz2=zlt-zlev(ig,l) |
346 |
|
|
enddo |
347 |
|
741355 |
zdz3=zlev(ig,lt+1)-zlt |
348 |
|
741355 |
zltdwn=zlev(ig,lt)-zdz3/2 |
349 |
|
|
zlmelup=zlmel+(zdz/2) |
350 |
|
741355 |
coefzlmel=Min(1.,(zlmelup-zltdwn)/zdz) |
351 |
|
|
zbuoyjam(ig,l)=1.*RG*(coefzlmel*(ztva_est(ig,l)- & |
352 |
|
|
& ztv(ig,lt))/ztv(ig,lt)+(1.-coefzlmel)*(ztva_est(ig,l)- & |
353 |
|
741355 |
& ztv(ig,lt-1))/ztv(ig,lt-1))+0.*zbuoy(ig,l) |
354 |
|
|
endif ! if (iflag_thermals_ed.lt.8) then |
355 |
|
|
|
356 |
|
|
!------------------------------------------------ |
357 |
|
|
!AJAM:nouveau calcul de w? |
358 |
|
|
!------------------------------------------------ |
359 |
|
|
zdz=zlev(ig,l+1)-zlev(ig,l) |
360 |
|
|
zdzbis=zlev(ig,l)-zlev(ig,l-1) |
361 |
|
|
zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l) |
362 |
|
|
|
363 |
|
741355 |
zw2fact=fact_epsilon*2.*zdz/(1.+betalpha) |
364 |
|
|
zw2factbis=fact_epsilon*2.*zdzbis/(1.+betalpha) |
365 |
|
741355 |
zdw2=afact*zbuoy(ig,l)/fact_epsilon |
366 |
|
741355 |
zdw2bis=afact*zbuoy(ig,l-1)/fact_epsilon |
367 |
|
|
! zdw2bis=0.5*(zdw2+zdw2bis) |
368 |
|
|
lm=Max(1,l-2) |
369 |
|
|
! zdw2=(afact/fact_epsilon)*((zdz/zdzbis)*zbuoy(ig,l) & |
370 |
|
|
! & +((zdzbis-zdz)/zdzbis)*zbuoy(ig,l-1)) |
371 |
|
|
! zdw2bis=(afact/fact_epsilon)*((zdz/zdzbis)*zbuoy(ig,l-1) & |
372 |
|
|
! & +((zdzbis-zdz)/zdzbis)*zbuoy(ig,l-1)) |
373 |
|
|
! w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2)+zdw2) |
374 |
|
|
! w_est(ig,l+1)=(zdz/zdzbis)*Max(0.0001,exp(-zw2fact)* & |
375 |
|
|
! & (w_est(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* & |
376 |
|
|
! & Max(0.0001,exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2) |
377 |
|
|
! w_est(ig,l+1)=Max(0.0001,(1-exp(-zw2fact))*zdw2+w_est(ig,l)*exp(-zw2fact)) |
378 |
|
|
|
379 |
|
|
!-------------------------------------------------- |
380 |
|
|
!AJ052014: J'ai remplac? w_est(ig,l) par zw2(ig,l) |
381 |
|
|
!-------------------------------------------------- |
382 |
✓✗ |
741355 |
if (iflag_thermals_ed==8) then |
383 |
|
|
! Ancienne version |
384 |
|
|
! w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* & |
385 |
|
|
! & (w_est(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* & |
386 |
|
|
! & (exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2)) |
387 |
|
|
|
388 |
|
741355 |
w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2)+zdw2) |
389 |
|
|
|
390 |
|
|
! Nouvelle version Arnaud |
391 |
|
|
else |
392 |
|
|
! w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* & |
393 |
|
|
! & (w_est(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* & |
394 |
|
|
! & (exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2)) |
395 |
|
|
|
396 |
|
|
w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2bis)+zdw2) |
397 |
|
|
|
398 |
|
|
! w_est(ig,l+1)=Max(0.0001,(zdz/(zdzbis+zdz))*(exp(-zw2fact)* & |
399 |
|
|
! & (w_est(ig,l)-zdw2bis)+zdw2)+(zdzbis/(zdzbis+zdz))* & |
400 |
|
|
! & (exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2bis)) |
401 |
|
|
|
402 |
|
|
|
403 |
|
|
|
404 |
|
|
! w_est(ig,l+1)=Max(0.0001,(w_est(ig,l)+zdw2bis*zw2fact)*exp(-zw2fact)) |
405 |
|
|
|
406 |
|
|
! w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact)+ & |
407 |
|
|
! & (zdzbis-zdz)/zdzbis*(zw2(ig,l-1)+zdw2bis*zw2factbis)*exp(-zw2factbis)) |
408 |
|
|
|
409 |
|
|
! w_est(ig,l+1)=Max(0.0001,exp(-zw2factbis)*(w_est(ig,l-1)-zdw2bis)+zdw2) |
410 |
|
|
|
411 |
|
|
endif |
412 |
|
|
|
413 |
|
|
|
414 |
✗✓ |
741355 |
if (iflag_thermals_ed<6) then |
415 |
|
|
zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l+1))/rhobarz(ig,l) |
416 |
|
|
! fact_epsilon=0.0005/(zalpha+0.025)**0.5 |
417 |
|
|
! fact_epsilon=Min(0.003,0.0004/(zalpha)**0.5) |
418 |
|
|
fact_epsilon=0.0002/(zalpha+0.1) |
419 |
|
|
zw2fact=fact_epsilon*2.*zdz/(1.+betalpha) |
420 |
|
|
zw2factbis=fact_epsilon*2.*zdzbis/(1.+betalpha) |
421 |
|
|
zdw2=afact*zbuoy(ig,l)/fact_epsilon |
422 |
|
|
zdw2bis=afact*zbuoy(ig,l-1)/fact_epsilon |
423 |
|
|
! w_est(ig,l+1)=Max(0.0001,(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact)) |
424 |
|
|
|
425 |
|
|
! w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* & |
426 |
|
|
! & (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* & |
427 |
|
|
! & (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2)) |
428 |
|
|
|
429 |
|
|
w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2bis)+zdw2) |
430 |
|
|
|
431 |
|
|
|
432 |
|
|
endif |
433 |
|
|
!-------------------------------------------------- |
434 |
|
|
!AJ052014: J'ai comment? ce if plus n?cessaire puisqu' |
435 |
|
|
!on fait max(0.0001,.....) |
436 |
|
|
!-------------------------------------------------- |
437 |
|
|
|
438 |
|
|
! if (w_est(ig,l+1).lt.0.) then |
439 |
|
|
! w_est(ig,l+1)=zw2(ig,l) |
440 |
|
|
! w_est(ig,l+1)=0.0001 |
441 |
|
|
! endif |
442 |
|
|
|
443 |
|
|
endif |
444 |
|
|
enddo |
445 |
|
|
|
446 |
|
|
|
447 |
|
|
!------------------------------------------------- |
448 |
|
|
!calcul des taux d'entrainement et de detrainement |
449 |
|
|
!------------------------------------------------- |
450 |
|
|
|
451 |
✓✓ |
10602720 |
do ig=1,ngrid |
452 |
✓✓ |
10602720 |
if (active(ig)) then |
453 |
|
|
|
454 |
|
|
! zw2m=max(0.5*(w_est(ig,l)+w_est(ig,l+1)),0.1) |
455 |
|
741355 |
zw2m=w_est(ig,l+1) |
456 |
|
|
! zw2m=zw2(ig,l) |
457 |
|
741355 |
zdz=zlev(ig,l+1)-zlev(ig,l) |
458 |
|
741355 |
zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l) |
459 |
|
|
! zbuoybis=zbuoy(ig,l)+RG*0.1/300. |
460 |
|
|
zbuoybis=zbuoy(ig,l) |
461 |
|
741355 |
zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l+1))/rhobarz(ig,l) |
462 |
|
741355 |
zdqt(ig,l)=max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l) |
463 |
|
|
|
464 |
|
|
|
465 |
|
|
! entr_star(ig,l)=f_star(ig,l)*zdz*zbetalpha*MAX(0., & |
466 |
|
|
! & afact*zbuoybis/zw2m - fact_epsilon ) |
467 |
|
|
|
468 |
|
|
! entr_star(ig,l)=MAX(0.,f_star(ig,l)*zdz*zbetalpha* & |
469 |
|
|
! & afact*zbuoybis/zw2m - fact_epsilon ) |
470 |
|
|
|
471 |
|
|
|
472 |
|
|
|
473 |
|
|
! zbuoyjam(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l) |
474 |
|
|
|
475 |
|
|
!========================================================================= |
476 |
|
|
! 4. Calcul de l'entrainement et du detrainement |
477 |
|
|
!========================================================================= |
478 |
|
|
|
479 |
|
|
! entr_star(ig,l)=f_star(ig,l)*zdz*zbetalpha*MAX(0., & |
480 |
|
|
! & afact*zbuoyjam(ig,l)/zw2m - fact_epsilon ) |
481 |
|
|
! entrbis=entr_star(ig,l) |
482 |
|
|
|
483 |
✗✓ |
741355 |
if (iflag_thermals_ed.lt.6) then |
484 |
|
|
fact_epsilon=0.0002/(zalpha+0.1) |
485 |
|
|
endif |
486 |
|
|
|
487 |
|
|
|
488 |
|
|
|
489 |
|
|
detr_star(ig,l)=f_star(ig,l)*zdz & |
490 |
|
|
& *( mix0 * 0.1 / (zalpha+0.001) & |
491 |
|
|
& + MAX(detr_min, -afact*zbetalpha*zbuoyjam(ig,l)/zw2m & |
492 |
|
741355 |
& + detr_q_coef*(zdqt(ig,l)/zw2m)**detr_q_power)) |
493 |
|
|
|
494 |
|
|
! detr_star(ig,l)=(zdz/zdzbis)*detr_star(ig,l)+ & |
495 |
|
|
! & ((zdzbis-zdz)/zdzbis)*detr_star(ig,l-1) |
496 |
|
|
|
497 |
|
|
zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l) |
498 |
|
|
|
499 |
|
|
entr_star(ig,l)=f_star(ig,l)*zdz* ( & |
500 |
|
|
& mix0 * 0.1 / (zalpha+0.001) & |
501 |
|
|
& + zbetalpha*MAX(entr_min, & |
502 |
|
741355 |
& afact*zbuoyjam(ig,l)/zw2m - fact_epsilon)) |
503 |
|
|
|
504 |
|
|
|
505 |
|
|
! entr_star(ig,l)=f_star(ig,l)*zdz* ( & |
506 |
|
|
! & mix0 * 0.1 / (zalpha+0.001) & |
507 |
|
|
! & + MAX(entr_min, & |
508 |
|
|
! & zbetalpha*afact*zbuoyjam(ig,l)/zw2m - fact_epsilon + & |
509 |
|
|
! & detr_q_coef*(zdqt(ig,l)/zw2m)**detr_q_power)) |
510 |
|
|
|
511 |
|
|
|
512 |
|
|
! entr_star(ig,l)=(zdz/zdzbis)*entr_star(ig,l)+ & |
513 |
|
|
! & ((zdzbis-zdz)/zdzbis)*entr_star(ig,l-1) |
514 |
|
|
|
515 |
|
|
! entr_star(ig,l)=Max(0.,f_star(ig,l)*zdz*zbetalpha* & |
516 |
|
|
! & afact*zbuoy(ig,l)/zw2m & |
517 |
|
|
! & - 1.*fact_epsilon) |
518 |
|
|
|
519 |
|
|
|
520 |
|
|
! En dessous de lalim, on prend le max de alim_star et entr_star pour |
521 |
|
|
! alim_star et 0 sinon |
522 |
✓✓ |
741355 |
if (l.lt.lalim(ig)) then |
523 |
|
203435 |
alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l)) |
524 |
|
203435 |
entr_star(ig,l)=0. |
525 |
|
|
endif |
526 |
|
|
! if (l.lt.lalim(ig).and.alim_star(ig,l)>alim_star(ig,l-1)) then |
527 |
|
|
! alim_star(ig,l)=entrbis |
528 |
|
|
! endif |
529 |
|
|
|
530 |
|
|
! print*,'alim0',zlev(ig,l),entr_star(ig,l),detr_star(ig,l),zw2m,zbuoy(ig,l),f_star(ig,l) |
531 |
|
|
! Calcul du flux montant normalise |
532 |
|
|
f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l) & |
533 |
|
741355 |
& -detr_star(ig,l) |
534 |
|
|
|
535 |
|
|
endif |
536 |
|
|
enddo |
537 |
|
|
|
538 |
|
|
|
539 |
|
|
!============================================================================ |
540 |
|
|
! 5. calcul de la vitesse verticale en melangeant Tl et qt du thermique |
541 |
|
|
!=========================================================================== |
542 |
|
|
|
543 |
✓✓✓✓
✓✓ |
10602720 |
activetmp(:)=active(:) .and. f_star(:,l+1)>1.e-10 |
544 |
✓✓ |
10602720 |
do ig=1,ngrid |
545 |
✓✓ |
10602720 |
if (activetmp(ig)) then |
546 |
|
|
Zsat=.false. |
547 |
|
|
ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+ & |
548 |
|
|
& (alim_star(ig,l)+entr_star(ig,l))*zthl(ig,l)) & |
549 |
|
604018 |
& /(f_star(ig,l+1)+detr_star(ig,l)) |
550 |
|
|
zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+ & |
551 |
|
|
& (alim_star(ig,l)+entr_star(ig,l))*po(ig,l)) & |
552 |
|
604018 |
& /(f_star(ig,l+1)+detr_star(ig,l)) |
553 |
|
|
|
554 |
|
|
endif |
555 |
|
|
enddo |
556 |
|
|
|
557 |
✓✓ |
10602720 |
ztemp(:)=zpspsk(:,l)*ztla(:,l) |
558 |
|
10656 |
call thermcell_qsat(ngrid,activetmp,pplev(:,l),ztemp,zqta(:,l),zqsatth(:,l)) |
559 |
✓✓ |
10602720 |
do ig=1,ngrid |
560 |
✓✓ |
10602720 |
if (activetmp(ig)) then |
561 |
|
|
! on ecrit de maniere conservative (sat ou non) |
562 |
|
|
! T = Tl +Lv/Cp ql |
563 |
|
604018 |
zqla(ig,l)=max(0.,zqta(ig,l)-zqsatth(ig,l)) |
564 |
|
604018 |
ztva(ig,l) = ztla(ig,l)*zpspsk(ig,l)+RLvCp*zqla(ig,l) |
565 |
|
604018 |
ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l) |
566 |
|
|
!on rajoute le calcul de zha pour diagnostiques (temp potentielle) |
567 |
|
604018 |
zha(ig,l) = ztva(ig,l) |
568 |
|
|
ztva(ig,l) = ztva(ig,l)*(1.+RETV*(zqta(ig,l) & |
569 |
|
604018 |
& -zqla(ig,l))-zqla(ig,l)) |
570 |
|
604018 |
zbuoy(ig,l)=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l) |
571 |
|
604018 |
zdz=zlev(ig,l+1)-zlev(ig,l) |
572 |
|
604018 |
zdzbis=zlev(ig,l)-zlev(ig,l-1) |
573 |
|
604018 |
zeps(ig,l)=(entr_star(ig,l)+alim_star(ig,l))/(f_star(ig,l)*zdz) |
574 |
|
|
!!!!!!! fact_epsilon=0.002 |
575 |
|
604018 |
zw2fact=fact_epsilon*2.*zdz/(1.+betalpha) |
576 |
|
|
zw2factbis=fact_epsilon*2.*zdzbis/(1.+betalpha) |
577 |
|
604018 |
zdw2= afact*zbuoy(ig,l)/(fact_epsilon) |
578 |
|
604018 |
zdw2bis= afact*zbuoy(ig,l-1)/(fact_epsilon) |
579 |
|
|
! zdw2=(afact/fact_epsilon)*((zdz/zdzbis)*zbuoy(ig,l) & |
580 |
|
|
! & +((zdzbis-zdz)/zdzbis)*zbuoy(ig,l-1)) |
581 |
|
|
! lm=Max(1,l-2) |
582 |
|
|
! zdw2bis=(afact/fact_epsilon)*((zdz/zdzbis)*zbuoy(ig,l-1) & |
583 |
|
|
! & +((zdzbis-zdz)/zdzbis)*zbuoy(ig,l-1)) |
584 |
|
|
! zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2bis)+zdw2) |
585 |
|
|
! zw2(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact)+ & |
586 |
|
|
! & (zdzbis-zdz)/zdzbis*(zw2(ig,l-1)+zdw2bis*zw2factbis)*exp(-zw2factbis)) |
587 |
|
|
! zw2(ig,l+1)=Max(0.0001,(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact)) |
588 |
|
|
! zw2(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* & |
589 |
|
|
! & (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* & |
590 |
|
|
! & (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2)) |
591 |
✓✗ |
604018 |
if (iflag_thermals_ed==8) then |
592 |
|
604018 |
zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2)+zdw2) |
593 |
|
|
else |
594 |
|
|
zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2bis)+zdw2) |
595 |
|
|
endif |
596 |
|
|
! zw2(ig,l+1)=Max(0.0001,(zdz/(zdz+zdzbis))*(exp(-zw2fact)* & |
597 |
|
|
! & (zw2(ig,l)-zdw2)+zdw2bis)+(zdzbis/(zdz+zdzbis))* & |
598 |
|
|
! & (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2bis)) |
599 |
|
|
|
600 |
|
|
|
601 |
✗✓ |
604018 |
if (iflag_thermals_ed.lt.6) then |
602 |
|
|
zalpha=f0(ig)*f_star(ig,l)/sqrt(zw2(ig,l+1))/rhobarz(ig,l) |
603 |
|
|
! fact_epsilon=0.0005/(zalpha+0.025)**0.5 |
604 |
|
|
! fact_epsilon=Min(0.003,0.0004/(zalpha)**0.5) |
605 |
|
|
fact_epsilon=0.0002/(zalpha+0.1)**1 |
606 |
|
|
zw2fact=fact_epsilon*2.*zdz/(1.+betalpha) |
607 |
|
|
zw2factbis=fact_epsilon*2.*zdzbis/(1.+betalpha) |
608 |
|
|
zdw2= afact*zbuoy(ig,l)/(fact_epsilon) |
609 |
|
|
zdw2bis= afact*zbuoy(ig,l-1)/(fact_epsilon) |
610 |
|
|
|
611 |
|
|
! zw2(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* & |
612 |
|
|
! & (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* & |
613 |
|
|
! & (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2)) |
614 |
|
|
! zw2(ig,l+1)=Max(0.0001,(zw2(ig,l)+zdw2*zw2fact)*exp(-zw2fact)) |
615 |
|
|
zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2bis)+zdw2) |
616 |
|
|
|
617 |
|
|
endif |
618 |
|
|
|
619 |
|
|
|
620 |
|
|
endif |
621 |
|
|
enddo |
622 |
|
|
|
623 |
✗✓ |
10656 |
if (prt_level.ge.20) print*,'coucou calcul detr 460: ig, l',ig, l |
624 |
|
|
! |
625 |
|
|
!=========================================================================== |
626 |
|
|
! 6. initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max |
627 |
|
|
!=========================================================================== |
628 |
|
|
|
629 |
|
10656 |
nbpb=0 |
630 |
✓✓ |
10602720 |
do ig=1,ngrid |
631 |
✓✓✗✓
|
10592064 |
if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then |
632 |
|
|
! stop'On tombe sur le cas particulier de thermcell_dry' |
633 |
|
|
! print*,'On tombe sur le cas particulier de thermcell_plume' |
634 |
|
|
nbpb=nbpb+1 |
635 |
|
|
zw2(ig,l+1)=0. |
636 |
|
|
linter(ig)=l+1 |
637 |
|
|
endif |
638 |
|
|
|
639 |
✗✓ |
10592064 |
if (zw2(ig,l+1).lt.0.) then |
640 |
|
|
linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l)) & |
641 |
|
|
& -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l)) |
642 |
|
|
zw2(ig,l+1)=0. |
643 |
|
|
!+CR:04/05/12:correction calcul linter pour calcul de zmax continu |
644 |
✓✓ |
10592064 |
elseif (f_star(ig,l+1).lt.0.) then |
645 |
|
|
linter(ig)=(l*(f_star(ig,l+1)-f_star(ig,l)) & |
646 |
|
137337 |
& -f_star(ig,l))/(f_star(ig,l+1)-f_star(ig,l)) |
647 |
|
137337 |
zw2(ig,l+1)=0. |
648 |
|
|
!fin CR:04/05/12 |
649 |
|
|
endif |
650 |
|
|
|
651 |
|
10592064 |
wa_moy(ig,l+1)=sqrt(zw2(ig,l+1)) |
652 |
|
|
|
653 |
✓✓ |
10602720 |
if (wa_moy(ig,l+1).gt.wmaxa(ig)) then |
654 |
|
|
! lmix est le niveau de la couche ou w (wa_moy) est maximum |
655 |
|
|
!on rajoute le calcul de lmix_bis |
656 |
✓✓ |
568355 |
if (zqla(ig,l).lt.1.e-10) then |
657 |
|
370846 |
lmix_bis(ig)=l+1 |
658 |
|
|
endif |
659 |
|
568355 |
lmix(ig)=l+1 |
660 |
|
568355 |
wmaxa(ig)=wa_moy(ig,l+1) |
661 |
|
|
endif |
662 |
|
|
enddo |
663 |
|
|
|
664 |
✗✓ |
10944 |
if (nbpb>0) then |
665 |
|
|
print*,'WARNING on tombe ',nbpb,' x sur un pb pour l=',l,' dans thermcell_plume' |
666 |
|
|
endif |
667 |
|
|
|
668 |
|
|
!========================================================================= |
669 |
|
|
! FIN DE LA BOUCLE VERTICALE |
670 |
|
|
enddo |
671 |
|
|
!========================================================================= |
672 |
|
|
|
673 |
|
|
!on recalcule alim_star_tot |
674 |
✓✓ |
286560 |
do ig=1,ngrid |
675 |
|
286560 |
alim_star_tot(ig)=0. |
676 |
|
|
enddo |
677 |
✓✓ |
286560 |
do ig=1,ngrid |
678 |
✓✓ |
628030 |
do l=1,lalim(ig)-1 |
679 |
|
627742 |
alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l) |
680 |
|
|
enddo |
681 |
|
|
enddo |
682 |
|
|
|
683 |
|
|
|
684 |
✗✓ |
288 |
if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l |
685 |
|
|
|
686 |
|
|
#undef wrgrads_thermcell |
687 |
|
|
#ifdef wrgrads_thermcell |
688 |
|
|
call wrgradsfi(1,nlay,entr_star(igout,1:nlay),'esta ','esta ') |
689 |
|
|
call wrgradsfi(1,nlay,detr_star(igout,1:nlay),'dsta ','dsta ') |
690 |
|
|
call wrgradsfi(1,nlay,zbuoy(igout,1:nlay),'buoy ','buoy ') |
691 |
|
|
call wrgradsfi(1,nlay,zdqt(igout,1:nlay),'dqt ','dqt ') |
692 |
|
|
call wrgradsfi(1,nlay,w_est(igout,1:nlay),'w_est ','w_est ') |
693 |
|
|
call wrgradsfi(1,nlay,w_est(igout,2:nlay+1),'w_es2 ','w_es2 ') |
694 |
|
|
call wrgradsfi(1,nlay,zw2(igout,1:nlay),'zw2A ','zw2A ') |
695 |
|
|
#endif |
696 |
|
|
|
697 |
|
|
|
698 |
|
288 |
RETURN |
699 |
|
|
end |
700 |
|
|
|
701 |
|
|
|
702 |
|
|
|
703 |
|
|
|
704 |
|
|
|
705 |
|
|
|
706 |
|
|
|
707 |
|
|
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
708 |
|
|
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
709 |
|
|
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
710 |
|
|
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
711 |
|
|
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
712 |
|
|
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
713 |
|
|
SUBROUTINE thermcell_plume_5B(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz, & |
714 |
|
|
& zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot, & |
715 |
|
|
& lalim,f0,detr_star,entr_star,f_star,csc,ztva, & |
716 |
|
|
& ztla,zqla,zqta,zha,zw2,w_est,ztva_est,zqsatth,lmix,lmix_bis,linter & |
717 |
|
|
& ,lev_out,lunout1,igout) |
718 |
|
|
!& ,lev_out,lunout1,igout,zbuoy,zbuoyjam) |
719 |
|
|
|
720 |
|
|
!-------------------------------------------------------------------------- |
721 |
|
|
!thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance |
722 |
|
|
! Version conforme a l'article de Rio et al. 2010. |
723 |
|
|
! Code ecrit par Catherine Rio, Arnaud Jam et Frederic Hourdin |
724 |
|
|
!-------------------------------------------------------------------------- |
725 |
|
|
|
726 |
|
|
USE lmdz_thermcell_ini, ONLY: prt_level,fact_thermals_ed_dz,iflag_thermals_ed,RLvCP,RETV,RG |
727 |
|
|
USE lmdz_thermcell_qsat, ONLY : thermcell_qsat |
728 |
|
|
IMPLICIT NONE |
729 |
|
|
|
730 |
|
|
INTEGER itap |
731 |
|
|
INTEGER lunout1,igout |
732 |
|
|
INTEGER ngrid,nlay |
733 |
|
|
REAL ptimestep |
734 |
|
|
REAL ztv(ngrid,nlay) |
735 |
|
|
REAL zthl(ngrid,nlay) |
736 |
|
|
REAL po(ngrid,nlay) |
737 |
|
|
REAL zl(ngrid,nlay) |
738 |
|
|
REAL rhobarz(ngrid,nlay) |
739 |
|
|
REAL zlev(ngrid,nlay+1) |
740 |
|
|
REAL pplev(ngrid,nlay+1) |
741 |
|
|
REAL pphi(ngrid,nlay) |
742 |
|
|
REAL zpspsk(ngrid,nlay) |
743 |
|
|
REAL alim_star(ngrid,nlay) |
744 |
|
|
REAL f0(ngrid) |
745 |
|
|
INTEGER lalim(ngrid) |
746 |
|
|
integer lev_out ! niveau pour les print |
747 |
|
|
integer nbpb |
748 |
|
|
|
749 |
|
|
real alim_star_tot(ngrid) |
750 |
|
|
|
751 |
|
|
REAL ztva(ngrid,nlay) |
752 |
|
|
REAL ztla(ngrid,nlay) |
753 |
|
|
REAL zqla(ngrid,nlay) |
754 |
|
|
REAL zqta(ngrid,nlay) |
755 |
|
|
REAL zha(ngrid,nlay) |
756 |
|
|
|
757 |
|
|
REAL detr_star(ngrid,nlay) |
758 |
|
|
REAL coefc |
759 |
|
|
REAL entr_star(ngrid,nlay) |
760 |
|
|
REAL detr(ngrid,nlay) |
761 |
|
|
REAL entr(ngrid,nlay) |
762 |
|
|
|
763 |
|
|
REAL csc(ngrid,nlay) |
764 |
|
|
|
765 |
|
|
REAL zw2(ngrid,nlay+1) |
766 |
|
|
REAL w_est(ngrid,nlay+1) |
767 |
|
|
REAL f_star(ngrid,nlay+1) |
768 |
|
|
REAL wa_moy(ngrid,nlay+1) |
769 |
|
|
|
770 |
|
|
REAL ztva_est(ngrid,nlay) |
771 |
|
|
REAL zqla_est(ngrid,nlay) |
772 |
|
|
REAL zqsatth(ngrid,nlay) |
773 |
|
|
REAL zta_est(ngrid,nlay) |
774 |
|
|
REAL zbuoyjam(ngrid,nlay) |
775 |
|
|
REAL ztemp(ngrid),zqsat(ngrid) |
776 |
|
|
REAL zdw2 |
777 |
|
|
REAL zw2modif |
778 |
|
|
REAL zw2fact |
779 |
|
|
REAL zeps(ngrid,nlay) |
780 |
|
|
|
781 |
|
|
REAL linter(ngrid) |
782 |
|
|
INTEGER lmix(ngrid) |
783 |
|
|
INTEGER lmix_bis(ngrid) |
784 |
|
|
REAL wmaxa(ngrid) |
785 |
|
|
|
786 |
|
|
INTEGER ig,l,k |
787 |
|
|
|
788 |
|
|
real zdz,zbuoy(ngrid,nlay),zalpha,gamma(ngrid,nlay),zdqt(ngrid,nlay),zw2m |
789 |
|
|
real zbuoybis |
790 |
|
|
real zcor,zdelta,zcvm5,qlbef,zdz2 |
791 |
|
|
real betalpha,zbetalpha |
792 |
|
|
real eps, afact |
793 |
|
|
logical Zsat |
794 |
|
|
LOGICAL active(ngrid),activetmp(ngrid) |
795 |
|
|
REAL fact_gamma,fact_epsilon,fact_gamma2,fact_epsilon2 |
796 |
|
|
REAL c2(ngrid,nlay) |
797 |
|
|
Zsat=.false. |
798 |
|
|
! Initialisation |
799 |
|
|
|
800 |
|
|
fact_epsilon=0.002 |
801 |
|
|
betalpha=0.9 |
802 |
|
|
afact=2./3. |
803 |
|
|
|
804 |
|
|
zbetalpha=betalpha/(1.+betalpha) |
805 |
|
|
|
806 |
|
|
|
807 |
|
|
! Initialisations des variables reeles |
808 |
|
|
if (1==1) then |
809 |
|
|
ztva(:,:)=ztv(:,:) |
810 |
|
|
ztva_est(:,:)=ztva(:,:) |
811 |
|
|
ztla(:,:)=zthl(:,:) |
812 |
|
|
zqta(:,:)=po(:,:) |
813 |
|
|
zha(:,:) = ztva(:,:) |
814 |
|
|
else |
815 |
|
|
ztva(:,:)=0. |
816 |
|
|
ztva_est(:,:)=0. |
817 |
|
|
ztla(:,:)=0. |
818 |
|
|
zqta(:,:)=0. |
819 |
|
|
zha(:,:) =0. |
820 |
|
|
endif |
821 |
|
|
|
822 |
|
|
zqla_est(:,:)=0. |
823 |
|
|
zqsatth(:,:)=0. |
824 |
|
|
zqla(:,:)=0. |
825 |
|
|
detr_star(:,:)=0. |
826 |
|
|
entr_star(:,:)=0. |
827 |
|
|
alim_star(:,:)=0. |
828 |
|
|
alim_star_tot(:)=0. |
829 |
|
|
csc(:,:)=0. |
830 |
|
|
detr(:,:)=0. |
831 |
|
|
entr(:,:)=0. |
832 |
|
|
zw2(:,:)=0. |
833 |
|
|
zbuoy(:,:)=0. |
834 |
|
|
zbuoyjam(:,:)=0. |
835 |
|
|
gamma(:,:)=0. |
836 |
|
|
zeps(:,:)=0. |
837 |
|
|
w_est(:,:)=0. |
838 |
|
|
f_star(:,:)=0. |
839 |
|
|
wa_moy(:,:)=0. |
840 |
|
|
linter(:)=1. |
841 |
|
|
! linter(:)=1. |
842 |
|
|
! Initialisation des variables entieres |
843 |
|
|
lmix(:)=1 |
844 |
|
|
lmix_bis(:)=2 |
845 |
|
|
wmaxa(:)=0. |
846 |
|
|
lalim(:)=1 |
847 |
|
|
|
848 |
|
|
|
849 |
|
|
!------------------------------------------------------------------------- |
850 |
|
|
! On ne considere comme actif que les colonnes dont les deux premieres |
851 |
|
|
! couches sont instables. |
852 |
|
|
!------------------------------------------------------------------------- |
853 |
|
|
active(:)=ztv(:,1)>ztv(:,2) |
854 |
|
|
|
855 |
|
|
!------------------------------------------------------------------------- |
856 |
|
|
! Definition de l'alimentation |
857 |
|
|
!------------------------------------------------------------------------- |
858 |
|
|
do l=1,nlay-1 |
859 |
|
|
do ig=1,ngrid |
860 |
|
|
if (ztv(ig,l)> ztv(ig,l+1) .and. ztv(ig,1)>=ztv(ig,l) ) then |
861 |
|
|
alim_star(ig,l)=MAX((ztv(ig,l)-ztv(ig,l+1)),0.) & |
862 |
|
|
& *sqrt(zlev(ig,l+1)) |
863 |
|
|
lalim(ig)=l+1 |
864 |
|
|
alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l) |
865 |
|
|
endif |
866 |
|
|
enddo |
867 |
|
|
enddo |
868 |
|
|
do l=1,nlay |
869 |
|
|
do ig=1,ngrid |
870 |
|
|
if (alim_star_tot(ig) > 1.e-10 ) then |
871 |
|
|
alim_star(ig,l)=alim_star(ig,l)/alim_star_tot(ig) |
872 |
|
|
endif |
873 |
|
|
enddo |
874 |
|
|
enddo |
875 |
|
|
alim_star_tot(:)=1. |
876 |
|
|
|
877 |
|
|
|
878 |
|
|
|
879 |
|
|
!------------------------------------------------------------------------------ |
880 |
|
|
! Calcul dans la premiere couche |
881 |
|
|
! On decide dans cette version que le thermique n'est actif que si la premiere |
882 |
|
|
! couche est instable. |
883 |
|
|
! Pourrait etre change si on veut que le thermiques puisse se d??clencher |
884 |
|
|
! dans une couche l>1 |
885 |
|
|
!------------------------------------------------------------------------------ |
886 |
|
|
do ig=1,ngrid |
887 |
|
|
! Le panache va prendre au debut les caracteristiques de l'air contenu |
888 |
|
|
! dans cette couche. |
889 |
|
|
if (active(ig)) then |
890 |
|
|
ztla(ig,1)=zthl(ig,1) |
891 |
|
|
zqta(ig,1)=po(ig,1) |
892 |
|
|
zqla(ig,1)=zl(ig,1) |
893 |
|
|
!cr: attention, prise en compte de f*(1)=1 |
894 |
|
|
f_star(ig,2)=alim_star(ig,1) |
895 |
|
|
zw2(ig,2)=2.*RG*(ztv(ig,1)-ztv(ig,2))/ztv(ig,2) & |
896 |
|
|
& *(zlev(ig,2)-zlev(ig,1)) & |
897 |
|
|
& *0.4*pphi(ig,1)/(pphi(ig,2)-pphi(ig,1)) |
898 |
|
|
w_est(ig,2)=zw2(ig,2) |
899 |
|
|
endif |
900 |
|
|
enddo |
901 |
|
|
! |
902 |
|
|
|
903 |
|
|
!============================================================================== |
904 |
|
|
!boucle de calcul de la vitesse verticale dans le thermique |
905 |
|
|
!============================================================================== |
906 |
|
|
do l=2,nlay-1 |
907 |
|
|
!============================================================================== |
908 |
|
|
|
909 |
|
|
|
910 |
|
|
! On decide si le thermique est encore actif ou non |
911 |
|
|
! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test |
912 |
|
|
do ig=1,ngrid |
913 |
|
|
active(ig)=active(ig) & |
914 |
|
|
& .and. zw2(ig,l)>1.e-10 & |
915 |
|
|
& .and. f_star(ig,l)+alim_star(ig,l)>1.e-10 |
916 |
|
|
enddo |
917 |
|
|
|
918 |
|
|
|
919 |
|
|
|
920 |
|
|
!--------------------------------------------------------------------------- |
921 |
|
|
! calcul des proprietes thermodynamiques et de la vitesse de la couche l |
922 |
|
|
! sans tenir compte du detrainement et de l'entrainement dans cette |
923 |
|
|
! couche |
924 |
|
|
! C'est a dire qu'on suppose |
925 |
|
|
! ztla(l)=ztla(l-1) et zqta(l)=zqta(l-1) |
926 |
|
|
! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer |
927 |
|
|
! avant) a l'alimentation pour avoir un calcul plus propre |
928 |
|
|
!--------------------------------------------------------------------------- |
929 |
|
|
|
930 |
|
|
ztemp(:)=zpspsk(:,l)*ztla(:,l-1) |
931 |
|
|
call thermcell_qsat(ngrid,active,pplev(:,l),ztemp,zqta(:,l-1),zqsat(:)) |
932 |
|
|
|
933 |
|
|
do ig=1,ngrid |
934 |
|
|
! print*,'active',active(ig),ig,l |
935 |
|
|
if(active(ig)) then |
936 |
|
|
zqla_est(ig,l)=max(0.,zqta(ig,l-1)-zqsat(ig)) |
937 |
|
|
ztva_est(ig,l) = ztla(ig,l-1)*zpspsk(ig,l)+RLvCp*zqla_est(ig,l) |
938 |
|
|
zta_est(ig,l)=ztva_est(ig,l) |
939 |
|
|
ztva_est(ig,l) = ztva_est(ig,l)/zpspsk(ig,l) |
940 |
|
|
ztva_est(ig,l) = ztva_est(ig,l)*(1.+RETV*(zqta(ig,l-1) & |
941 |
|
|
& -zqla_est(ig,l))-zqla_est(ig,l)) |
942 |
|
|
|
943 |
|
|
!------------------------------------------------ |
944 |
|
|
!AJAM:nouveau calcul de w? |
945 |
|
|
!------------------------------------------------ |
946 |
|
|
zdz=zlev(ig,l+1)-zlev(ig,l) |
947 |
|
|
zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l) |
948 |
|
|
|
949 |
|
|
zw2fact=fact_epsilon*2.*zdz/(1.+betalpha) |
950 |
|
|
zdw2=(afact)*zbuoy(ig,l)/(fact_epsilon) |
951 |
|
|
w_est(ig,l+1)=Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2)+zdw2) |
952 |
|
|
|
953 |
|
|
|
954 |
|
|
if (w_est(ig,l+1).lt.0.) then |
955 |
|
|
w_est(ig,l+1)=zw2(ig,l) |
956 |
|
|
endif |
957 |
|
|
endif |
958 |
|
|
enddo |
959 |
|
|
|
960 |
|
|
|
961 |
|
|
!------------------------------------------------- |
962 |
|
|
!calcul des taux d'entrainement et de detrainement |
963 |
|
|
!------------------------------------------------- |
964 |
|
|
|
965 |
|
|
do ig=1,ngrid |
966 |
|
|
if (active(ig)) then |
967 |
|
|
|
968 |
|
|
zw2m=max(0.5*(w_est(ig,l)+w_est(ig,l+1)),0.1) |
969 |
|
|
zw2m=w_est(ig,l+1) |
970 |
|
|
zdz=zlev(ig,l+1)-zlev(ig,l) |
971 |
|
|
zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l) |
972 |
|
|
! zbuoybis=zbuoy(ig,l)+RG*0.1/300. |
973 |
|
|
zbuoybis=zbuoy(ig,l) |
974 |
|
|
zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l+1))/rhobarz(ig,l) |
975 |
|
|
zdqt(ig,l)=max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l) |
976 |
|
|
|
977 |
|
|
|
978 |
|
|
entr_star(ig,l)=f_star(ig,l)*zdz* zbetalpha*MAX(0., & |
979 |
|
|
& afact*zbuoybis/zw2m - fact_epsilon ) |
980 |
|
|
|
981 |
|
|
|
982 |
|
|
detr_star(ig,l)=f_star(ig,l)*zdz & |
983 |
|
|
& *MAX(1.e-3, -afact*zbetalpha*zbuoy(ig,l)/zw2m & |
984 |
|
|
& + 0.012*(zdqt(ig,l)/zw2m)**0.5 ) |
985 |
|
|
|
986 |
|
|
! En dessous de lalim, on prend le max de alim_star et entr_star pour |
987 |
|
|
! alim_star et 0 sinon |
988 |
|
|
if (l.lt.lalim(ig)) then |
989 |
|
|
alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l)) |
990 |
|
|
entr_star(ig,l)=0. |
991 |
|
|
endif |
992 |
|
|
|
993 |
|
|
! Calcul du flux montant normalise |
994 |
|
|
f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l) & |
995 |
|
|
& -detr_star(ig,l) |
996 |
|
|
|
997 |
|
|
endif |
998 |
|
|
enddo |
999 |
|
|
|
1000 |
|
|
|
1001 |
|
|
!---------------------------------------------------------------------------- |
1002 |
|
|
!calcul de la vitesse verticale en melangeant Tl et qt du thermique |
1003 |
|
|
!--------------------------------------------------------------------------- |
1004 |
|
|
activetmp(:)=active(:) .and. f_star(:,l+1)>1.e-10 |
1005 |
|
|
do ig=1,ngrid |
1006 |
|
|
if (activetmp(ig)) then |
1007 |
|
|
Zsat=.false. |
1008 |
|
|
ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+ & |
1009 |
|
|
& (alim_star(ig,l)+entr_star(ig,l))*zthl(ig,l)) & |
1010 |
|
|
& /(f_star(ig,l+1)+detr_star(ig,l)) |
1011 |
|
|
zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+ & |
1012 |
|
|
& (alim_star(ig,l)+entr_star(ig,l))*po(ig,l)) & |
1013 |
|
|
& /(f_star(ig,l+1)+detr_star(ig,l)) |
1014 |
|
|
|
1015 |
|
|
endif |
1016 |
|
|
enddo |
1017 |
|
|
|
1018 |
|
|
ztemp(:)=zpspsk(:,l)*ztla(:,l) |
1019 |
|
|
call thermcell_qsat(ngrid,activetmp,pplev(:,l),ztemp,zqta(:,l),zqsatth(:,l)) |
1020 |
|
|
|
1021 |
|
|
do ig=1,ngrid |
1022 |
|
|
if (activetmp(ig)) then |
1023 |
|
|
! on ecrit de maniere conservative (sat ou non) |
1024 |
|
|
! T = Tl +Lv/Cp ql |
1025 |
|
|
zqla(ig,l)=max(0.,zqta(ig,l)-zqsatth(ig,l)) |
1026 |
|
|
ztva(ig,l) = ztla(ig,l)*zpspsk(ig,l)+RLvCp*zqla(ig,l) |
1027 |
|
|
ztva(ig,l) = ztva(ig,l)/zpspsk(ig,l) |
1028 |
|
|
!on rajoute le calcul de zha pour diagnostiques (temp potentielle) |
1029 |
|
|
zha(ig,l) = ztva(ig,l) |
1030 |
|
|
ztva(ig,l) = ztva(ig,l)*(1.+RETV*(zqta(ig,l) & |
1031 |
|
|
& -zqla(ig,l))-zqla(ig,l)) |
1032 |
|
|
zbuoy(ig,l)=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l) |
1033 |
|
|
zdz=zlev(ig,l+1)-zlev(ig,l) |
1034 |
|
|
zeps(ig,l)=(entr_star(ig,l)+alim_star(ig,l))/(f_star(ig,l)*zdz) |
1035 |
|
|
|
1036 |
|
|
zw2fact=fact_epsilon*2.*zdz/(1.+betalpha) |
1037 |
|
|
zdw2=afact*zbuoy(ig,l)/(fact_epsilon) |
1038 |
|
|
zw2(ig,l+1)=Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2)+zdw2) |
1039 |
|
|
endif |
1040 |
|
|
enddo |
1041 |
|
|
|
1042 |
|
|
if (prt_level.ge.20) print*,'coucou calcul detr 460: ig, l',ig, l |
1043 |
|
|
! |
1044 |
|
|
!--------------------------------------------------------------------------- |
1045 |
|
|
!initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max |
1046 |
|
|
!--------------------------------------------------------------------------- |
1047 |
|
|
|
1048 |
|
|
nbpb=0 |
1049 |
|
|
do ig=1,ngrid |
1050 |
|
|
if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then |
1051 |
|
|
! stop'On tombe sur le cas particulier de thermcell_dry' |
1052 |
|
|
! print*,'On tombe sur le cas particulier de thermcell_plume' |
1053 |
|
|
nbpb=nbpb+1 |
1054 |
|
|
zw2(ig,l+1)=0. |
1055 |
|
|
linter(ig)=l+1 |
1056 |
|
|
endif |
1057 |
|
|
|
1058 |
|
|
if (zw2(ig,l+1).lt.0.) then |
1059 |
|
|
linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l)) & |
1060 |
|
|
& -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l)) |
1061 |
|
|
zw2(ig,l+1)=0. |
1062 |
|
|
elseif (f_star(ig,l+1).lt.0.) then |
1063 |
|
|
linter(ig)=(l*(f_star(ig,l+1)-f_star(ig,l)) & |
1064 |
|
|
& -f_star(ig,l))/(f_star(ig,l+1)-f_star(ig,l)) |
1065 |
|
|
! print*,"linter plume", linter(ig) |
1066 |
|
|
zw2(ig,l+1)=0. |
1067 |
|
|
endif |
1068 |
|
|
|
1069 |
|
|
wa_moy(ig,l+1)=sqrt(zw2(ig,l+1)) |
1070 |
|
|
|
1071 |
|
|
if (wa_moy(ig,l+1).gt.wmaxa(ig)) then |
1072 |
|
|
! lmix est le niveau de la couche ou w (wa_moy) est maximum |
1073 |
|
|
!on rajoute le calcul de lmix_bis |
1074 |
|
|
if (zqla(ig,l).lt.1.e-10) then |
1075 |
|
|
lmix_bis(ig)=l+1 |
1076 |
|
|
endif |
1077 |
|
|
lmix(ig)=l+1 |
1078 |
|
|
wmaxa(ig)=wa_moy(ig,l+1) |
1079 |
|
|
endif |
1080 |
|
|
enddo |
1081 |
|
|
|
1082 |
|
|
if (nbpb>0) then |
1083 |
|
|
print*,'WARNING on tombe ',nbpb,' x sur un pb pour l=',l,' dans thermcell_plume' |
1084 |
|
|
endif |
1085 |
|
|
|
1086 |
|
|
!========================================================================= |
1087 |
|
|
! FIN DE LA BOUCLE VERTICALE |
1088 |
|
|
enddo |
1089 |
|
|
!========================================================================= |
1090 |
|
|
|
1091 |
|
|
!on recalcule alim_star_tot |
1092 |
|
|
do ig=1,ngrid |
1093 |
|
|
alim_star_tot(ig)=0. |
1094 |
|
|
enddo |
1095 |
|
|
do ig=1,ngrid |
1096 |
|
|
do l=1,lalim(ig)-1 |
1097 |
|
|
alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l) |
1098 |
|
|
enddo |
1099 |
|
|
enddo |
1100 |
|
|
|
1101 |
|
|
|
1102 |
|
|
if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l |
1103 |
|
|
|
1104 |
|
|
#undef wrgrads_thermcell |
1105 |
|
|
#ifdef wrgrads_thermcell |
1106 |
|
|
call wrgradsfi(1,nlay,entr_star(igout,1:nlay),'esta ','esta ') |
1107 |
|
|
call wrgradsfi(1,nlay,detr_star(igout,1:nlay),'dsta ','dsta ') |
1108 |
|
|
call wrgradsfi(1,nlay,zbuoy(igout,1:nlay),'buoy ','buoy ') |
1109 |
|
|
call wrgradsfi(1,nlay,zdqt(igout,1:nlay),'dqt ','dqt ') |
1110 |
|
|
call wrgradsfi(1,nlay,w_est(igout,1:nlay),'w_est ','w_est ') |
1111 |
|
|
call wrgradsfi(1,nlay,w_est(igout,2:nlay+1),'w_es2 ','w_es2 ') |
1112 |
|
|
call wrgradsfi(1,nlay,zw2(igout,1:nlay),'zw2A ','zw2A ') |
1113 |
|
|
#endif |
1114 |
|
|
|
1115 |
|
|
|
1116 |
|
|
return |
1117 |
|
|
end |
1118 |
|
|
END MODULE lmdz_thermcell_plume_6A |