Context Navigation

calltherm.F90 @ 1509

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

Modified closure for the new physics package, new values for the iflag_coupl parameter
FH

Fermeture modifiée pour la nouvelle physique, nouvelles valeurs définies pour
le paramètre iflag_coupl
FH

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 1496 2011-03-11 09:44:05Z idelkadi $
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,iflag_coupl &
220	& ,Ale,Alp,lalim_conv,wght_th &
221	& ,zmax0,f0,zw2,fraca,ztv,zpspsk &
222	& ,ztla,zthl)
223	if (prt_level.gt.10) write(lunout,*)'Apres thermcell_main OK'
224	else
225	abort_message = 'Cas des thermiques non prevu'
226	CALL abort_gcm (modname,abort_message,1)
227	endif
228
229	flag_bidouille_stratocu=iflag_thermals.eq.14.or.iflag_thermals.eq.16
230
231	fact(:)=0.
232	DO i=1,klon
233	logexpr1(i)=flag_bidouille_stratocu.or.weak_inversion(i).gt.0.5
234	IF(logexpr1(i)) fact(i)=1./REAL(nsplit_thermals)
235	ENDDO
236
237	DO k=1,klev
238	! transformation de la derivee en tendance
239	d_t_the(:,k)=d_t_the(:,k)dtimefact(:)
240	d_u_the(:,k)=d_u_the(:,k)dtimefact(:)
241	d_v_the(:,k)=d_v_the(:,k)dtimefact(:)
242	d_q_the(:,k)=d_q_the(:,k)dtimefact(:)
243	fm_therm(:,k)=fm_therm(:,k) &
244	& +zfm_therm(:,k)*fact(:)
245	entr_therm(:,k)=entr_therm(:,k) &
246	& +zentr_therm(:,k)*fact(:)
247	detr_therm(:,k)=detr_therm(:,k) &
248	& +zdetr_therm(:,k)*fact(:)
249	ENDDO
250	fm_therm(:,klev+1)=0.
251
252
253
254	! accumulation de la tendance
255	d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_the(:,:)
256	d_u_ajs(:,:)=d_u_ajs(:,:)+d_u_the(:,:)
257	d_v_ajs(:,:)=d_v_ajs(:,:)+d_v_the(:,:)
258	d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_the(:,:)
259
260	! incrementation des variables meteo
261	t_seri(:,:) = t_seri(:,:) + d_t_the(:,:)
262	u_seri(:,:) = u_seri(:,:) + d_u_the(:,:)
263	v_seri(:,:) = v_seri(:,:) + d_v_the(:,:)
264	qmemoire(:,:)=q_seri(:,:)
265	q_seri(:,:) = q_seri(:,:) + d_q_the(:,:)
266	if (prt_level.gt.10) write(lunout,*)'Apres apres thermcell_main OK'
267
268	DO i=1,klon
269	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)
270	fm_therm(i,klev+1)=0.
271	Ale_bl(i)=Ale_bl(i)+Ale(i)/REAL(nsplit_thermals)
272	! write(22,*)'ALE CALLTHERM',Ale_bl(i),Ale(i)
273	Alp_bl(i)=Alp_bl(i)+Alp(i)/REAL(nsplit_thermals)
274	! write(23,*)'ALP CALLTHERM',Alp_bl(i),Alp(i)
275	ENDDO
276
277	!IM 060508 marche pas comme cela !!! enddo ! isplit
278
279	! tests sur les valeurs negatives de l'eau
280	nbptspb=0
281	DO k = 1, klev
282	DO i = 1, klon
283	logexpr2(i,k)=.not.q_seri(i,k).ge.0.
284	if (logexpr2(i,k)) then
285	q_seri(i,k)=1.e-15
286	nbptspb=nbptspb+1
287	! if (prt_level.ge.10) then
288	! print*,'WARN eau<0 apres therm i=',i,' k=',k &
289	! & ,' dq,q',d_q_the(i,k),q_seri(i,k), &
290	! & 'fm=',zfm_therm(i,k),'entr=',entr_therm(i,k)
291	endif
292	ENDDO
293	ENDDO
294	IF(nbptspb.GT.0) print*,'Number of points with q_seri(i,k)<=0 ',nbptspb
295	! tests sur les valeurs de la temperature
296	nbptspb=0
297	DO k = 1, klev
298	DO i = 1, klon
299	logexpr2(i,k)=t_seri(i,k).lt.50..or.t_seri(i,k).gt.370.
300	if (logexpr2(i,k)) nbptspb=nbptspb+1
301	! if ((t_seri(i,k).lt.50.) .or. &
302	! & (t_seri(i,k).gt.370.)) then
303	! print*,'WARN temp apres therm i=',i,' k=',k &
304	! & ,' t_seri',t_seri(i,k)
305	! CALL abort
306	! endif
307	ENDDO
308	ENDDO
309	IF(nbptspb.GT.0) print*,'Number of points with q_seri(i,k)<=0 ',nbptspb
310	enddo ! isplit
311
312	!
313	!***************************************************************
314	! calcul du flux ascencant conservatif
315	! print*,'<<<<calcul flux ascendant conservatif'
316
317	fmc_therm=0.
318	do k=1,klev
319	do i=1,klon
320	if (entr_therm(i,k).gt.0.) then
321	fmc_therm(i,k+1)=fmc_therm(i,k)+entr_therm(i,k)
322	else
323	fmc_therm(i,k+1)=fmc_therm(i,k)
324	endif
325	detrc_therm(i,k)=(fmc_therm(i,k+1)-fm_therm(i,k+1)) &
326	& -(fmc_therm(i,k)-fm_therm(i,k))
327	enddo
328	enddo
329
330
331	!****************************************************************
332	! calcul de l'humidite dans l'ascendance
333	! print*,'<<<<calcul de lhumidite dans thermique'
334	!CR:on ne le calcule que pour le cas sec
335	if (iflag_thermals.le.11) then
336	do i=1,klon
337	zqasc(i,1)=q_seri(i,1)
338	do k=2,klev
339	if (fmc_therm(i,k+1).gt.1.e-6) then
340	zqasc(i,k)=(fmc_therm(i,k)*zqasc(i,k-1) &
341	& +entr_therm(i,k)*q_seri(i,k))/fmc_therm(i,k+1)
342	!CR:test on asseche le thermique
343	! zqasc(i,k)=zqasc(i,k)/2.
344	! else
345	! zqasc(i,k)=q_seri(i,k)
346	endif
347	enddo
348	enddo
349
350
351	! calcul de l'eau condensee dans l'ascendance
352	! print*,'<<<<calcul de leau condensee dans thermique'
353	do i=1,klon
354	do k=1,klev
355	clwcon0(i,k)=zqasc(i,k)-zqsat(i,k)
356	if (clwcon0(i,k).lt.0. .or. &
357	& (fm_therm(i,k+1)+detrc_therm(i,k)).lt.1.e-6) then
358	clwcon0(i,k)=0.
359	endif
360	enddo
361	enddo
362	else
363	do i=1,klon
364	do k=1,klev
365	clwcon0(i,k)=zqla(i,k)
366	if (clwcon0(i,k).lt.0. .or. &
367	& (fm_therm(i,k+1)+detrc_therm(i,k)).lt.1.e-6) then
368	clwcon0(i,k)=0.
369	endif
370	enddo
371	enddo
372	endif
373	!*******************************************************************
374
375
376	!jyg Protection contre les temperatures nulles
377	do i=1,klon
378	do k=1,klev
379	if (ztla(i,k) .lt. 1.e-10) fraca(i,k) =0.
380	enddo
381	enddo
382
383
384	return
385
386	end

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

Download in other formats:

Original Format