Context Navigation

calltherm.F90 @ 5394

Last change on this file since 5394 was 1311, checked in by Laurent Fairhead, 15 years ago

Modifications to thermals for TKE transport

Modifications aux thermiques pour le transport de la TKE

pbl_surface_mode.F90 : ok_flux_surf=.false. seulement pour klon>1
physiq.F : option iflag_pbl=10 pour transporter la TKE avec les thermiques.
calltherm.F90 : passage de iflag_thermals_ed en argument pour thermcell_plume
thermcell_main.F90 : Appel a plusieurs version de thermcell_plume en option
thermcell_plume.F90 : plusieurs versions dans le meme fichier (temporaire)
thermcell_height.F90 : verrue pour les cas ou les thermiques montent tout

en haut

yamada4 : inclusion de la diffusion verticale en option iflag_pbl=9

+ variables anciennement common, puis save/allocatable, remises en local

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 12.8 KB

Line
1	!
2	! $Id: calltherm.F90 1311 2010-02-18 13:14:02Z evignon $
3	!
4	subroutine calltherm(dtime &
5	& ,pplay,paprs,pphi,weak_inversion &
6	& ,u_seri,v_seri,t_seri,q_seri,zqsat,debut &
7	& ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs &
8	& ,fm_therm,entr_therm,detr_therm,zqasc,clwcon0,lmax,ratqscth, &
9	& ratqsdiff,zqsatth,Ale_bl,Alp_bl,lalim_conv,wght_th, &
10	& zmax0,f0,zw2,fraca)
11
12	USE dimphy
13	implicit none
14	#include "dimensions.h"
15	!#include "dimphy.h"
16	#include "thermcell.h"
17	#include "iniprint.h"
18
19	!IM 140508
20	INTEGER itap
21	REAL dtime
22	LOGICAL debut
23	LOGICAL logexpr0, logexpr2(klon,klev), logexpr1(klon)
24	REAL fact(klon)
25	INTEGER nbptspb
26
27	REAL u_seri(klon,klev),v_seri(klon,klev)
28	REAL t_seri(klon,klev),q_seri(klon,klev),qmemoire(klon,klev)
29	REAL weak_inversion(klon)
30	REAL paprs(klon,klev+1)
31	REAL pplay(klon,klev)
32	REAL pphi(klon,klev)
33	real zlev(klon,klev+1)
34	!test: on sort lentr et a* pour alimenter KE
35	REAL wght_th(klon,klev)
36	INTEGER lalim_conv(klon)
37	REAL zw2(klon,klev+1),fraca(klon,klev+1)
38
39	!FH Update Thermiques
40	REAL d_t_ajs(klon,klev), d_q_ajs(klon,klev)
41	REAL d_u_ajs(klon,klev),d_v_ajs(klon,klev)
42	real fm_therm(klon,klev+1)
43	real entr_therm(klon,klev),detr_therm(klon,klev)
44
45	!********************************************************
46	! declarations
47	LOGICAL flag_bidouille_stratocu
48	real fmc_therm(klon,klev+1),zqasc(klon,klev)
49	real zqla(klon,klev)
50	real zqta(klon,klev)
51	real wmax_sec(klon)
52	real zmax_sec(klon)
53	real f_sec(klon)
54	real detrc_therm(klon,klev)
55	! FH WARNING : il semble que ces save ne servent a rien
56	! save fmc_therm, detrc_therm
57	real clwcon0(klon,klev)
58	real zqsat(klon,klev)
59	real zw_sec(klon,klev+1)
60	integer lmix_sec(klon)
61	integer lmax(klon)
62	real ratqscth(klon,klev)
63	real ratqsdiff(klon,klev)
64	real zqsatth(klon,klev)
65	!nouvelles variables pour la convection
66	real Ale_bl(klon)
67	real Alp_bl(klon)
68	real Ale(klon)
69	real Alp(klon)
70	!RC
71	!on garde le zmax du pas de temps precedent
72	real zmax0(klon), f0(klon)
73	!********************************************************
74
75
76	! variables locales
77	REAL d_t_the(klon,klev), d_q_the(klon,klev)
78	REAL d_u_the(klon,klev),d_v_the(klon,klev)
79	!
80	real zfm_therm(klon,klev+1),zdt
81	real zentr_therm(klon,klev),zdetr_therm(klon,klev)
82	! FH A VERIFIER : SAVE INUTILES
83	! save zentr_therm,zfm_therm
84
85	character (len=20) :: modname='calltherm'
86	character (len=80) :: abort_message
87
88	integer i,k
89	logical, save :: first=.true.
90	!$OMP THREADPRIVATE(first)
91	!********************************************************
92	if (first) then
93	itap=0
94	first=.false.
95	endif
96
97	! Incrementer le compteur de la physique
98	itap = itap + 1
99
100	! Modele du thermique
101	! ===================
102	! print*,'thermiques: WARNING on passe t au lieu de t_seri'
103
104
105	! On prend comme valeur initiale des thermiques la valeur du pas
106	! de temps precedent
107	zfm_therm(:,:)=fm_therm(:,:)
108	zdetr_therm(:,:)=detr_therm(:,:)
109	zentr_therm(:,:)=entr_therm(:,:)
110
111	! On reinitialise les flux de masse a zero pour le cumul en
112	! cas de splitting
113	fm_therm(:,:)=0.
114	entr_therm(:,:)=0.
115	detr_therm(:,:)=0.
116
117	Ale_bl(:)=0.
118	Alp_bl(:)=0.
119	if (prt_level.ge.10) then
120	print*,'thermV4 nsplit: ',nsplit_thermals,' weak_inversion'
121	endif
122
123	! tests sur les valeurs negatives de l'eau
124	logexpr0=prt_level.ge.10
125	nbptspb=0
126	do k=1,klev
127	do i=1,klon
128	! Attention teste abderr 19-03-09
129	! logexpr2(i,k)=.not.q_seri(i,k).ge.0.
130	logexpr2(i,k)=.not.q_seri(i,k).ge.1.e-15
131	if (logexpr2(i,k)) then
132	q_seri(i,k)=1.e-15
133	nbptspb=nbptspb+1
134	endif
135	! if (logexpr0) &
136	! & print*,'WARN eau<0 avant therm i=',i,' k=',k &
137	! & ,' dq,q',d_q_the(i,k),q_seri(i,k)
138	enddo
139	enddo
140	if(nbptspb.GT.0) print*,'Number of points with q_seri(i,k)<=0 ',nbptspb
141
142	zdt=dtime/REAL(nsplit_thermals)
143	do isplit=1,nsplit_thermals
144
145	if (iflag_thermals.eq.1) then
146	CALL thermcell_2002(klon,klev,zdt &
147	& ,pplay,paprs,pphi &
148	& ,u_seri,v_seri,t_seri,q_seri &
149	& ,d_u_the,d_v_the,d_t_the,d_q_the &
150	& ,zfm_therm,zentr_therm &
151	& ,r_aspect_thermals,30.,w2di_thermals &
152	& ,tau_thermals,3)
153	else if (iflag_thermals.eq.2) then
154	CALL thermcell_sec(klon,klev,zdt &
155	& ,pplay,paprs,pphi,zlev &
156	& ,u_seri,v_seri,t_seri,q_seri &
157	& ,d_u_the,d_v_the,d_t_the,d_q_the &
158	& ,zfm_therm,zentr_therm &
159	& ,r_aspect_thermals,30.,w2di_thermals &
160	& ,tau_thermals,3)
161	else if (iflag_thermals.eq.3) then
162	CALL thermcell(klon,klev,zdt &
163	& ,pplay,paprs,pphi &
164	& ,u_seri,v_seri,t_seri,q_seri &
165	& ,d_u_the,d_v_the,d_t_the,d_q_the &
166	& ,zfm_therm,zentr_therm &
167	& ,r_aspect_thermals,l_mix_thermals,w2di_thermals &
168	& ,tau_thermals,3)
169	else if (iflag_thermals.eq.10) then
170	CALL thermcell_eau(klon,klev,zdt &
171	& ,pplay,paprs,pphi &
172	& ,u_seri,v_seri,t_seri,q_seri &
173	& ,d_u_the,d_v_the,d_t_the,d_q_the &
174	& ,zfm_therm,zentr_therm &
175	& ,r_aspect_thermals,l_mix_thermals,w2di_thermals &
176	& ,tau_thermals,3)
177	else if (iflag_thermals.eq.11) then
178	abort_message = 'cas non prevu dans calltherm'
179	CALL abort_gcm (modname,abort_message,1)
180
181	! CALL thermcell_pluie(klon,klev,zdt &
182	! & ,pplay,paprs,pphi,zlev &
183	! & ,u_seri,v_seri,t_seri,q_seri &
184	! & ,d_u_the,d_v_the,d_t_the,d_q_the &
185	! & ,zfm_therm,zentr_therm,zqla &
186	! & ,r_aspect_thermals,l_mix_thermals,w2di_thermals &
187	! & ,tau_thermals,3)
188	else if (iflag_thermals.eq.12) then
189	CALL calcul_sec(klon,klev,zdt &
190	& ,pplay,paprs,pphi,zlev &
191	& ,u_seri,v_seri,t_seri,q_seri &
192	& ,zmax_sec,wmax_sec,zw_sec,lmix_sec &
193	& ,r_aspect_thermals,l_mix_thermals,w2di_thermals &
194	& ,tau_thermals)
195	else if (iflag_thermals==13.or.iflag_thermals==14) then
196	CALL thermcellV0_main(itap,klon,klev,zdt &
197	& ,pplay,paprs,pphi,debut &
198	& ,u_seri,v_seri,t_seri,q_seri &
199	& ,d_u_the,d_v_the,d_t_the,d_q_the &
200	& ,zfm_therm,zentr_therm,zdetr_therm,zqasc,zqla,lmax &
201	& ,ratqscth,ratqsdiff,zqsatth &
202	& ,r_aspect_thermals,l_mix_thermals &
203	& ,tau_thermals,Ale,Alp,lalim_conv,wght_th &
204	& ,zmax0,f0,zw2,fraca)
205	else if (iflag_thermals==15.or.iflag_thermals==16) then
206
207	! print*,'THERM iflag_thermas_ed=',iflag_thermals_ed
208	CALL thermcell_main(itap,klon,klev,zdt &
209	& ,pplay,paprs,pphi,debut &
210	& ,u_seri,v_seri,t_seri,q_seri &
211	& ,d_u_the,d_v_the,d_t_the,d_q_the &
212	& ,zfm_therm,zentr_therm,zdetr_therm,zqasc,zqla,lmax &
213	& ,ratqscth,ratqsdiff,zqsatth &
214	& ,r_aspect_thermals,l_mix_thermals &
215	& ,tau_thermals,iflag_thermals_ed,Ale,Alp,lalim_conv,wght_th &
216	& ,zmax0,f0,zw2,fraca)
217	if (prt_level.gt.10) write(lunout,*)'Apres thermcell_main OK'
218	else
219	abort_message = 'Cas des thermiques non prevu'
220	CALL abort_gcm (modname,abort_message,1)
221	endif
222
223	flag_bidouille_stratocu=iflag_thermals.eq.14.or.iflag_thermals.eq.16
224
225	fact(:)=0.
226	DO i=1,klon
227	logexpr1(i)=flag_bidouille_stratocu.or.weak_inversion(i).gt.0.5
228	IF(logexpr1(i)) fact(i)=1./REAL(nsplit_thermals)
229	ENDDO
230
231	DO k=1,klev
232	! transformation de la derivee en tendance
233	d_t_the(:,k)=d_t_the(:,k)dtimefact(:)
234	d_u_the(:,k)=d_u_the(:,k)dtimefact(:)
235	d_v_the(:,k)=d_v_the(:,k)dtimefact(:)
236	d_q_the(:,k)=d_q_the(:,k)dtimefact(:)
237	fm_therm(:,k)=fm_therm(:,k) &
238	& +zfm_therm(:,k)*fact(:)
239	entr_therm(:,k)=entr_therm(:,k) &
240	& +zentr_therm(:,k)*fact(:)
241	detr_therm(:,k)=detr_therm(:,k) &
242	& +zdetr_therm(:,k)*fact(:)
243	ENDDO
244	fm_therm(:,klev+1)=0.
245
246
247
248	! accumulation de la tendance
249	d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_the(:,:)
250	d_u_ajs(:,:)=d_u_ajs(:,:)+d_u_the(:,:)
251	d_v_ajs(:,:)=d_v_ajs(:,:)+d_v_the(:,:)
252	d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_the(:,:)
253
254	! incrementation des variables meteo
255	t_seri(:,:) = t_seri(:,:) + d_t_the(:,:)
256	u_seri(:,:) = u_seri(:,:) + d_u_the(:,:)
257	v_seri(:,:) = v_seri(:,:) + d_v_the(:,:)
258	qmemoire(:,:)=q_seri(:,:)
259	q_seri(:,:) = q_seri(:,:) + d_q_the(:,:)
260	if (prt_level.gt.10) write(lunout,*)'Apres apres thermcell_main OK'
261
262	DO i=1,klon
263	if(prt_level.GE.10) print*,'calltherm i Alp_bl Alp Ale_bl Ale',i,Alp_bl(i),Alp(i),Ale_bl(i),Ale(i)
264	fm_therm(i,klev+1)=0.
265	Ale_bl(i)=Ale_bl(i)+Ale(i)/REAL(nsplit_thermals)
266	! write(22,*)'ALE CALLTHERM',Ale_bl(i),Ale(i)
267	Alp_bl(i)=Alp_bl(i)+Alp(i)/REAL(nsplit_thermals)
268	! write(23,*)'ALP CALLTHERM',Alp_bl(i),Alp(i)
269	ENDDO
270
271	!IM 060508 marche pas comme cela !!! enddo ! isplit
272
273	! tests sur les valeurs negatives de l'eau
274	nbptspb=0
275	DO k = 1, klev
276	DO i = 1, klon
277	logexpr2(i,k)=.not.q_seri(i,k).ge.0.
278	if (logexpr2(i,k)) then
279	q_seri(i,k)=1.e-15
280	nbptspb=nbptspb+1
281	! if (prt_level.ge.10) then
282	! print*,'WARN eau<0 apres therm i=',i,' k=',k &
283	! & ,' dq,q',d_q_the(i,k),q_seri(i,k), &
284	! & 'fm=',zfm_therm(i,k),'entr=',entr_therm(i,k)
285	endif
286	ENDDO
287	ENDDO
288	IF(nbptspb.GT.0) print*,'Number of points with q_seri(i,k)<=0 ',nbptspb
289	! tests sur les valeurs de la temperature
290	nbptspb=0
291	DO k = 1, klev
292	DO i = 1, klon
293	logexpr2(i,k)=t_seri(i,k).lt.50..or.t_seri(i,k).gt.370.
294	if (logexpr2(i,k)) nbptspb=nbptspb+1
295	! if ((t_seri(i,k).lt.50.) .or. &
296	! & (t_seri(i,k).gt.370.)) then
297	! print*,'WARN temp apres therm i=',i,' k=',k &
298	! & ,' t_seri',t_seri(i,k)
299	! CALL abort
300	! endif
301	ENDDO
302	ENDDO
303	IF(nbptspb.GT.0) print*,'Number of points with q_seri(i,k)<=0 ',nbptspb
304	enddo ! isplit
305
306	!
307	!***************************************************************
308	! calcul du flux ascencant conservatif
309	! print*,'<<<<calcul flux ascendant conservatif'
310
311	fmc_therm=0.
312	do k=1,klev
313	do i=1,klon
314	if (entr_therm(i,k).gt.0.) then
315	fmc_therm(i,k+1)=fmc_therm(i,k)+entr_therm(i,k)
316	else
317	fmc_therm(i,k+1)=fmc_therm(i,k)
318	endif
319	detrc_therm(i,k)=(fmc_therm(i,k+1)-fm_therm(i,k+1)) &
320	& -(fmc_therm(i,k)-fm_therm(i,k))
321	enddo
322	enddo
323
324
325	!****************************************************************
326	! calcul de l'humidite dans l'ascendance
327	! print*,'<<<<calcul de lhumidite dans thermique'
328	!CR:on ne le calcule que pour le cas sec
329	if (iflag_thermals.le.11) then
330	do i=1,klon
331	zqasc(i,1)=q_seri(i,1)
332	do k=2,klev
333	if (fmc_therm(i,k+1).gt.1.e-6) then
334	zqasc(i,k)=(fmc_therm(i,k)*zqasc(i,k-1) &
335	& +entr_therm(i,k)*q_seri(i,k))/fmc_therm(i,k+1)
336	!CR:test on asseche le thermique
337	! zqasc(i,k)=zqasc(i,k)/2.
338	! else
339	! zqasc(i,k)=q_seri(i,k)
340	endif
341	enddo
342	enddo
343
344
345	! calcul de l'eau condensee dans l'ascendance
346	! print*,'<<<<calcul de leau condensee dans thermique'
347	do i=1,klon
348	do k=1,klev
349	clwcon0(i,k)=zqasc(i,k)-zqsat(i,k)
350	if (clwcon0(i,k).lt.0. .or. &
351	& (fm_therm(i,k+1)+detrc_therm(i,k)).lt.1.e-6) then
352	clwcon0(i,k)=0.
353	endif
354	enddo
355	enddo
356	else
357	do i=1,klon
358	do k=1,klev
359	clwcon0(i,k)=zqla(i,k)
360	if (clwcon0(i,k).lt.0. .or. &
361	& (fm_therm(i,k+1)+detrc_therm(i,k)).lt.1.e-6) then
362	clwcon0(i,k)=0.
363	endif
364	enddo
365	enddo
366	endif
367	!*******************************************************************
368
369
370	return
371
372	end

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format