Context Navigation

calltherm.F90 @ 5327

Last change on this file since 5327 was 1428, checked in by Laurent Fairhead, 14 years ago

Problems with some initialisations in the new physics package. These are
temporary fixes. Further investigations needed. Trac ticket issued. JYG

Problèmes d'initialisations dans la nouvelle physique. Ces corrections
sont temporaires, il faut trouver la vraie raison des plantages. Un ticket
trac est créé. JYG

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

File size: 13.2 KB

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

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

Download in other formats:

Original Format