Context Navigation

calltherm.F90 @ 96

Last change on this file since 96 was 96, checked in by lfita, 10 years ago
Wrong check_var!!! on Ale_bl
File size: 18.2 KB

Line
1	!
2	! $Id: calltherm.F90 1795 2013-07-18 08:20:28Z emillour $
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,ztv,zpspsk,ztla,zthl &
11	!!! nrlmd le 10/04/2012
12	& ,pbl_tke,pctsrf,omega,airephy &
13	& ,zlcl_th,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &
14	& ,n2,s2,ale_bl_stat &
15	& ,therm_tke_max,env_tke_max &
16	& ,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &
17	& ,alp_bl_conv,alp_bl_stat &
18	!!! fin nrlmd le 10/04/2012
19	& ,zqla,ztva )
20
21	USE dimphy
22	USE indice_sol_mod
23
24	implicit none
25	#include "dimensions.h"
26	!#include "dimphy.h"
27	#include "thermcell.h"
28	#include "iniprint.h"
29
30
31	!IM 140508
32	INTEGER, SAVE :: itap
33	!$OMP THREADPRIVATE(itap)
34	REAL dtime
35	LOGICAL debut
36	LOGICAL logexpr0, logexpr2(klon,klev), logexpr1(klon)
37	REAL fact(klon)
38	INTEGER nbptspb
39
40	REAL u_seri(klon,klev),v_seri(klon,klev)
41	REAL t_seri(klon,klev),q_seri(klon,klev),qmemoire(klon,klev)
42	REAL weak_inversion(klon)
43	REAL paprs(klon,klev+1)
44	REAL pplay(klon,klev)
45	REAL pphi(klon,klev)
46	real zlev(klon,klev+1)
47	!test: on sort lentr et a* pour alimenter KE
48	REAL wght_th(klon,klev)
49	INTEGER lalim_conv(klon)
50	REAL zw2(klon,klev+1),fraca(klon,klev+1)
51
52	!FH Update Thermiques
53	REAL d_t_ajs(klon,klev), d_q_ajs(klon,klev)
54	REAL d_u_ajs(klon,klev),d_v_ajs(klon,klev)
55	real fm_therm(klon,klev+1)
56	real entr_therm(klon,klev),detr_therm(klon,klev)
57
58	!********************************************************
59	! declarations
60	LOGICAL flag_bidouille_stratocu
61	real fmc_therm(klon,klev+1),zqasc(klon,klev)
62	real zqla(klon,klev)
63	real zqta(klon,klev)
64	real ztv(klon,klev),ztva(klon,klev)
65	real zpspsk(klon,klev)
66	real ztla(klon,klev)
67	real zthl(klon,klev)
68	real wmax_sec(klon)
69	real zmax_sec(klon)
70	real f_sec(klon)
71	real detrc_therm(klon,klev)
72	! FH WARNING : il semble que ces save ne servent a rien
73	! save fmc_therm, detrc_therm
74	real clwcon0(klon,klev)
75	real zqsat(klon,klev)
76	real zw_sec(klon,klev+1)
77	integer lmix_sec(klon)
78	integer lmax(klon)
79	real ratqscth(klon,klev)
80	real ratqsdiff(klon,klev)
81	real zqsatth(klon,klev)
82	!nouvelles variables pour la convection
83	real Ale_bl(klon)
84	real Alp_bl(klon)
85	real Ale(klon)
86	real Alp(klon)
87	!RC
88	!on garde le zmax du pas de temps precedent
89	real zmax0(klon), f0(klon)
90
91	!!! nrlmd le 10/04/2012
92	real pbl_tke(klon,klev+1,nbsrf)
93	real pctsrf(klon,nbsrf)
94	real omega(klon,klev)
95	real airephy(klon)
96	real zlcl_th(klon),fraca0(klon),w0(klon),w_conv(klon)
97	real therm_tke_max0(klon),env_tke_max0(klon)
98	real n2(klon),s2(klon)
99	real ale_bl_stat(klon)
100	real therm_tke_max(klon,klev),env_tke_max(klon,klev)
101	real alp_bl_det(klon),alp_bl_fluct_m(klon),alp_bl_fluct_tke(klon),alp_bl_conv(klon),alp_bl_stat(klon)
102	!!! fin nrlmd le 10/04/2012
103
104	!********************************************************
105
106
107	! variables locales
108	REAL d_t_the(klon,klev), d_q_the(klon,klev)
109	REAL d_u_the(klon,klev),d_v_the(klon,klev)
110	!
111	real zfm_therm(klon,klev+1),zdt
112	real zentr_therm(klon,klev),zdetr_therm(klon,klev)
113	! FH A VERIFIER : SAVE INUTILES
114	! save zentr_therm,zfm_therm
115
116	character (len=20) :: modname='calltherm'
117	character (len=80) :: abort_message
118
119	integer i,k
120	logical, save :: first=.true.
121	!$OMP THREADPRIVATE(first)
122	!********************************************************
123
124	! Lluis
125	INTEGER :: llp
126	CHARACTER(LEN=50) :: lvarname, lfname
127	REAL :: largest
128
129	llp = 734
130	lfname = 'physiq'
131	largest = 10.e5
132	if (first) then
133	itap=0
134	first=.false.
135	endif
136
137	! L. Fita, LMD July 2014. Initializing tendencies and variables
138	d_t_the = 0.
139	d_q_the = 0.
140	d_u_the = 0.
141	d_v_the = 0.
142	Alp = 0.
143	Ale = 0.
144
145	PRINT *,' Lluis SUm(Alp): ', SUM(Alp),' SUM Ale: ',SUM(Ale)
146	IF (SUM(Alp) /= 0.) STOP
147	IF (SUM(Ale) /= 0.) STOP
148
149	! Incrementer le compteur de la physique
150	itap = itap + 1
151
152	! Modele du thermique
153	! ===================
154	! print*,'thermiques: WARNING on passe t au lieu de t_seri'
155
156	lfname='beginning calltherm'
157	lvarname = 't_seri'
158	CALL check_var3D(lfname, lvarname, t_seri, klon, klev, largest, .FALSE.)
159
160	! On prend comme valeur initiale des thermiques la valeur du pas
161	! de temps precedent
162	zfm_therm(:,:)=fm_therm(:,:)
163	zdetr_therm(:,:)=detr_therm(:,:)
164	zentr_therm(:,:)=entr_therm(:,:)
165
166	! On reinitialise les flux de masse a zero pour le cumul en
167	! cas de splitting
168	fm_therm(:,:)=0.
169	entr_therm(:,:)=0.
170	detr_therm(:,:)=0.
171
172	Ale_bl(:)=0.
173	Alp_bl(:)=0.
174	if (prt_level.ge.10) then
175	print*,'thermV4 nsplit: ',nsplit_thermals,' weak_inversion'
176	endif
177
178	! tests sur les valeurs negatives de l'eau
179	logexpr0=prt_level.ge.10
180	nbptspb=0
181	do k=1,klev
182	do i=1,klon
183	! Attention teste abderr 19-03-09
184	! logexpr2(i,k)=.not.q_seri(i,k).ge.0.
185	logexpr2(i,k)=.not.q_seri(i,k).ge.1.e-15
186	if (logexpr2(i,k)) then
187	q_seri(i,k)=1.e-15
188	nbptspb=nbptspb+1
189	endif
190	! if (logexpr0) &
191	! & print*,'WARN eau<0 avant therm i=',i,' k=',k &
192	! & ,' dq,q',d_q_the(i,k),q_seri(i,k)
193	enddo
194	enddo
195	if(nbptspb.GT.0) print*,'Number of points with q_seri(i,k)<=0 ',nbptspb
196	lfname='calltherm before thermcell_main'
197	lvarname = 't_seri'
198	CALL check_var3D(lfname, lvarname, t_seri, klon, klev, largest, .FALSE.)
199	lvarname = 'd_t_the'
200	CALL check_var3D(lfname, lvarname, d_t_the, klon, klev, largest, .FALSE.)
201	lvarname = 'Alp'
202	CALL check_var(lfname, lvarname, Alp, klon, largest*10.e5, .FALSE.)
203	lvarname = 'Ale'
204	CALL check_var(lfname, lvarname, Ale, klon, largest*10.e5, .FALSE.)
205
206	zdt=dtime/REAL(nsplit_thermals)
207	do isplit=1,nsplit_thermals
208
209	if (iflag_thermals.eq.1) then
210	CALL thermcell_2002(klon,klev,zdt &
211	& ,pplay,paprs,pphi &
212	& ,u_seri,v_seri,t_seri,q_seri &
213	& ,d_u_the,d_v_the,d_t_the,d_q_the &
214	& ,zfm_therm,zentr_therm &
215	& ,r_aspect_thermals,30.,w2di_thermals &
216	& ,tau_thermals)
217	else if (iflag_thermals.eq.2) then
218	CALL thermcell_sec(klon,klev,zdt &
219	& ,pplay,paprs,pphi,zlev &
220	& ,u_seri,v_seri,t_seri,q_seri &
221	& ,d_u_the,d_v_the,d_t_the,d_q_the &
222	& ,zfm_therm,zentr_therm &
223	& ,r_aspect_thermals,30.,w2di_thermals &
224	& ,tau_thermals)
225	else if (iflag_thermals.eq.3) then
226	CALL thermcell(klon,klev,zdt &
227	& ,pplay,paprs,pphi &
228	& ,u_seri,v_seri,t_seri,q_seri &
229	& ,d_u_the,d_v_the,d_t_the,d_q_the &
230	& ,zfm_therm,zentr_therm &
231	& ,r_aspect_thermals,l_mix_thermals,w2di_thermals &
232	& ,tau_thermals)
233	else if (iflag_thermals.eq.10) then
234	CALL thermcell_eau(klon,klev,zdt &
235	& ,pplay,paprs,pphi &
236	& ,u_seri,v_seri,t_seri,q_seri &
237	& ,d_u_the,d_v_the,d_t_the,d_q_the &
238	& ,zfm_therm,zentr_therm &
239	& ,r_aspect_thermals,l_mix_thermals,w2di_thermals &
240	& ,tau_thermals)
241	else if (iflag_thermals.eq.11) then
242	abort_message = 'cas non prevu dans calltherm'
243	CALL abort_gcm (modname,abort_message,1)
244
245	! CALL thermcell_pluie(klon,klev,zdt &
246	! & ,pplay,paprs,pphi,zlev &
247	! & ,u_seri,v_seri,t_seri,q_seri &
248	! & ,d_u_the,d_v_the,d_t_the,d_q_the &
249	! & ,zfm_therm,zentr_therm,zqla &
250	! & ,r_aspect_thermals,l_mix_thermals,w2di_thermals &
251	! & ,tau_thermals,3)
252	else if (iflag_thermals.eq.12) then
253	CALL calcul_sec(klon,klev,zdt &
254	& ,pplay,paprs,pphi,zlev &
255	& ,u_seri,v_seri,t_seri,q_seri &
256	& ,zmax_sec,wmax_sec,zw_sec,lmix_sec &
257	& ,r_aspect_thermals,l_mix_thermals,w2di_thermals &
258	& ,tau_thermals)
259	else if (iflag_thermals==13.or.iflag_thermals==14) then
260	CALL thermcellV0_main(itap,klon,klev,zdt &
261	& ,pplay,paprs,pphi,debut &
262	& ,u_seri,v_seri,t_seri,q_seri &
263	& ,d_u_the,d_v_the,d_t_the,d_q_the &
264	& ,zfm_therm,zentr_therm,zdetr_therm,zqasc,zqla,lmax &
265	& ,ratqscth,ratqsdiff,zqsatth &
266	& ,r_aspect_thermals,l_mix_thermals &
267	& ,tau_thermals,Ale,Alp,lalim_conv,wght_th &
268	& ,zmax0,f0,zw2,fraca)
269	else if (iflag_thermals>=15.and.iflag_thermals<=18) then
270
271	! print*,'THERM iflag_thermas_ed=',iflag_thermals_ed
272
273	CALL thermcell_main(itap,klon,klev,zdt &
274	& ,pplay,paprs,pphi,debut &
275	& ,u_seri,v_seri,t_seri,q_seri &
276	& ,d_u_the,d_v_the,d_t_the,d_q_the &
277	& ,zfm_therm,zentr_therm,zdetr_therm,zqasc,zqla,lmax &
278	& ,ratqscth,ratqsdiff,zqsatth &
279	! & ,r_aspect_thermals,l_mix_thermals &
280	! & ,tau_thermals,iflag_thermals_ed,iflag_coupl &
281	& ,Ale,Alp,lalim_conv,wght_th &
282	& ,zmax0,f0,zw2,fraca,ztv,zpspsk &
283	& ,ztla,zthl &
284	!!! nrlmd le 10/04/2012
285	& ,pbl_tke,pctsrf,omega,airephy &
286	& ,zlcl_th,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &
287	& ,n2,s2,ale_bl_stat &
288	& ,therm_tke_max,env_tke_max &
289	& ,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &
290	& ,alp_bl_conv,alp_bl_stat &
291	!!! fin nrlmd le 10/04/2012
292	& ,ztva )
293	if (prt_level.gt.10) write(lunout,*)'Apres thermcell_main OK'
294	else
295	abort_message = 'Cas des thermiques non prevu'
296	CALL abort_gcm (modname,abort_message,1)
297	endif
298
299	! Attention : les noms sont contre intuitif.
300	! flag_bidouille_stratocu est .true. si on ne fait pas de bidouille.
301	! Il aurait mieux valu avoir un nobidouille_stratocu
302	! Et pour simplifier :
303	! nobidouille_stratocu=.not.(iflag_thermals==13.or.iflag_thermals=15)
304	! Ce serait bien de changer, mai en prenant le temps de vÃ©rifier que ca
305	! fait bien ce qu'on croit.
306	lfname='calltherm after thermcell_main'
307	lvarname = 't_seri'
308	CALL check_var3D(lfname, lvarname, t_seri, klon, klev, largest, .FALSE.)
309	lvarname = 'd_t_the'
310	CALL check_var3D(lfname, lvarname, d_t_the, klon, klev, largest, .FALSE.)
311	lvarname = 'Alp'
312	CALL check_var(lfname, lvarname, Alp, klon, largest*10.e5, .FALSE.)
313	lvarname = 'Ale'
314	CALL check_var(lfname, lvarname, Ale, klon, largest*10.e5, .FALSE.)
315
316	flag_bidouille_stratocu=iflag_thermals<=12.or.iflag_thermals==14.or.iflag_thermals==16.or.iflag_thermals==18
317
318	if (iflag_thermals<=12) then
319	lmax=1
320	do k=1,klev-1
321	zdetr_therm(:,k)=zentr_therm(:,k)+zfm_therm(:,k)-zfm_therm(:,k+1)
322	enddo
323	endif
324
325	fact(:)=0.
326	DO i=1,klon
327	logexpr1(i)=flag_bidouille_stratocu.or.weak_inversion(i).gt.0.5
328	IF(logexpr1(i)) fact(i)=1./REAL(nsplit_thermals)
329	ENDDO
330
331	DO k=1,klev
332	! transformation de la derivee en tendance
333	d_t_the(:,k)=d_t_the(:,k)dtimefact(:)
334	d_u_the(:,k)=d_u_the(:,k)dtimefact(:)
335	d_v_the(:,k)=d_v_the(:,k)dtimefact(:)
336	d_q_the(:,k)=d_q_the(:,k)dtimefact(:)
337	fm_therm(:,k)=fm_therm(:,k) &
338	& +zfm_therm(:,k)*fact(:)
339	entr_therm(:,k)=entr_therm(:,k) &
340	& +zentr_therm(:,k)*fact(:)
341	detr_therm(:,k)=detr_therm(:,k) &
342	& +zdetr_therm(:,k)*fact(:)
343	ENDDO
344	fm_therm(:,klev+1)=0.
345
346	lfname='calltherm after transforming from derivate to tendency'
347	lvarname = 't_seri'
348	CALL check_var3D(lfname, lvarname, t_seri, klon, klev, largest, .FALSE.)
349	lvarname = 'd_t_the'
350	CALL check_var3D(lfname, lvarname, d_t_the, klon, klev, largest, .FALSE.)
351
352
353	! accumulation de la tendance
354	d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_the(:,:)
355	d_u_ajs(:,:)=d_u_ajs(:,:)+d_u_the(:,:)
356	d_v_ajs(:,:)=d_v_ajs(:,:)+d_v_the(:,:)
357	d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_the(:,:)
358
359	! incrementation des variables meteo
360	t_seri(:,:) = t_seri(:,:) + d_t_the(:,:)
361	u_seri(:,:) = u_seri(:,:) + d_u_the(:,:)
362	v_seri(:,:) = v_seri(:,:) + d_v_the(:,:)
363	qmemoire(:,:)=q_seri(:,:)
364	q_seri(:,:) = q_seri(:,:) + d_q_the(:,:)
365	if (prt_level.gt.10) write(lunout,*)'Apres apres thermcell_main OK'
366	lfname='calltherm after add tendencies'
367	lvarname = 't_seri'
368	CALL check_var3D(lfname, lvarname, t_seri, klon, klev, largest, .FALSE.)
369	lvarname = 'd_t_the'
370	CALL check_var3D(lfname, lvarname, d_t_the, klon, klev, largest, .FALSE.)
371	lvarname = 'Alp_bl'
372	CALL check_var(lfname, lvarname, Alp_bl, klon, largest, .FALSE.)
373	lvarname = 'Alp'
374	CALL check_var(lfname, lvarname, Alp, klon, largest*10.e5, .FALSE.)
375	lvarname = 'Ale_bl'
376	CALL check_var(lfname, lvarname, Ale_bl, klon, largest, .FALSE.)
377	lvarname = 'Ale'
378	CALL check_var(lfname, lvarname, Ale, klon, largest*10.e5, .FALSE.)
379
380	DO i=1,klon
381	fm_therm(i,klev+1)=0.
382	Ale_bl(i)=Ale_bl(i)+Ale(i)/REAL(nsplit_thermals)
383	! write(22,*)'ALE CALLTHERM',Ale_bl(i),Ale(i)
384	Alp_bl(i)=Alp_bl(i)+Alp(i)/REAL(nsplit_thermals)
385	! write(23,*)'ALP CALLTHERM',Alp_bl(i),Alp(i)
386	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)
387	ENDDO
388	lfname='calltherm after computing Al[p/e]_bl'
389	lvarname = 't_seri'
390	CALL check_var3D(lfname, lvarname, t_seri, klon, klev, largest, .FALSE.)
391	lvarname = 'd_t_the'
392	CALL check_var3D(lfname, lvarname, d_t_the, klon, klev, largest, .FALSE.)
393	lvarname = 'Alp_bl'
394	CALL check_var(lfname, lvarname, Alp_bl, klon, largest, .FALSE.)
395	lvarname = 'Alp'
396	CALL check_var(lfname, lvarname, Alp, klon, largest*10.e5, .FALSE.)
397	lvarname = 'Ale_bl'
398	CALL check_var(lfname, lvarname, Ale_bl, klon, largest, .FALSE.)
399	lvarname = 'Ale'
400	CALL check_var(lfname, lvarname, Ale, klon, largest*10.e5, .FALSE.)
401
402	!IM 060508 marche pas comme cela !!! enddo ! isplit
403
404	! tests sur les valeurs negatives de l'eau
405	nbptspb=0
406	DO k = 1, klev
407	DO i = 1, klon
408	logexpr2(i,k)=.not.q_seri(i,k).ge.0.
409	if (logexpr2(i,k)) then
410	q_seri(i,k)=1.e-15
411	nbptspb=nbptspb+1
412	! if (prt_level.ge.10) then
413	! print*,'WARN eau<0 apres therm i=',i,' k=',k &
414	! & ,' dq,q',d_q_the(i,k),q_seri(i,k), &
415	! & 'fm=',zfm_therm(i,k),'entr=',entr_therm(i,k)
416	endif
417	ENDDO
418	ENDDO
419	IF(nbptspb.GT.0) print*,'Number of points with q_seri(i,k)<=0 ',nbptspb
420	! tests sur les valeurs de la temperature
421	nbptspb=0
422	DO k = 1, klev
423	DO i = 1, klon
424	logexpr2(i,k)=t_seri(i,k).lt.50..or.t_seri(i,k).gt.370.
425	if (logexpr2(i,k)) nbptspb=nbptspb+1
426	! if ((t_seri(i,k).lt.50.) .or. &
427	! & (t_seri(i,k).gt.370.)) then
428	! print*,'WARN temp apres therm i=',i,' k=',k &
429	! & ,' t_seri',t_seri(i,k)
430	! CALL abort
431	! endif
432	ENDDO
433	ENDDO
434	IF(nbptspb.GT.0) print*,'Number of points with q_seri(i,k)<=0 ',nbptspb
435	enddo ! isplit
436
437	!
438	!***************************************************************
439	! calcul du flux ascencant conservatif
440	! print*,'<<<<calcul flux ascendant conservatif'
441
442	fmc_therm=0.
443	do k=1,klev
444	do i=1,klon
445	if (entr_therm(i,k).gt.0.) then
446	fmc_therm(i,k+1)=fmc_therm(i,k)+entr_therm(i,k)
447	else
448	fmc_therm(i,k+1)=fmc_therm(i,k)
449	endif
450	detrc_therm(i,k)=(fmc_therm(i,k+1)-fm_therm(i,k+1)) &
451	& -(fmc_therm(i,k)-fm_therm(i,k))
452	enddo
453	enddo
454
455
456	!****************************************************************
457	! calcul de l'humidite dans l'ascendance
458	! print*,'<<<<calcul de lhumidite dans thermique'
459	!CR:on ne le calcule que pour le cas sec
460	if (iflag_thermals.le.11) then
461	do i=1,klon
462	zqasc(i,1)=q_seri(i,1)
463	do k=2,klev
464	if (fmc_therm(i,k+1).gt.1.e-6) then
465	zqasc(i,k)=(fmc_therm(i,k)*zqasc(i,k-1) &
466	& +entr_therm(i,k)*q_seri(i,k))/fmc_therm(i,k+1)
467	!CR:test on asseche le thermique
468	! zqasc(i,k)=zqasc(i,k)/2.
469	! else
470	! zqasc(i,k)=q_seri(i,k)
471	endif
472	enddo
473	enddo
474
475
476	! calcul de l'eau condensee dans l'ascendance
477	! print*,'<<<<calcul de leau condensee dans thermique'
478	do i=1,klon
479	do k=1,klev
480	clwcon0(i,k)=zqasc(i,k)-zqsat(i,k)
481	if (clwcon0(i,k).lt.0. .or. &
482	& (fm_therm(i,k+1)+detrc_therm(i,k)).lt.1.e-6) then
483	clwcon0(i,k)=0.
484	endif
485	enddo
486	enddo
487	else
488	do i=1,klon
489	do k=1,klev
490	clwcon0(i,k)=zqla(i,k)
491	if (clwcon0(i,k).lt.0. .or. &
492	& (fm_therm(i,k+1)+detrc_therm(i,k)).lt.1.e-6) then
493	clwcon0(i,k)=0.
494	endif
495	enddo
496	enddo
497	endif
498	!*******************************************************************
499
500
501	!jyg Protection contre les temperatures nulles
502	do i=1,klon
503	do k=1,klev
504	if (ztla(i,k) .lt. 1.e-10) fraca(i,k) =0.
505	enddo
506	enddo
507
508
509	return
510
511	end

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

Download in other formats:

Original Format