Context Navigation

calltherm.F90 @ 3577

Last change on this file since 3577 was 1790, checked in by Ehouarn Millour, 11 years ago

Ré-implémentation des commissions écrasées par la rev 1785. (Notamment lmdz1d.F)
UG
...................................................

Reimplementation of the commitions undonned by rev 1785. (Especially lmdy1d.F)
UG

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

File size: 15.0 KB

Line
1	!
2	! $Id: calltherm.F90 1790 2013-07-17 09:24:04Z fhourdin $
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	if (first) then
124	itap=0
125	first=.false.
126	endif
127
128	! Incrementer le compteur de la physique
129	itap = itap + 1
130
131	! Modele du thermique
132	! ===================
133	! print*,'thermiques: WARNING on passe t au lieu de t_seri'
134
135
136	! On prend comme valeur initiale des thermiques la valeur du pas
137	! de temps precedent
138	zfm_therm(:,:)=fm_therm(:,:)
139	zdetr_therm(:,:)=detr_therm(:,:)
140	zentr_therm(:,:)=entr_therm(:,:)
141
142	! On reinitialise les flux de masse a zero pour le cumul en
143	! cas de splitting
144	fm_therm(:,:)=0.
145	entr_therm(:,:)=0.
146	detr_therm(:,:)=0.
147
148	Ale_bl(:)=0.
149	Alp_bl(:)=0.
150	if (prt_level.ge.10) then
151	print*,'thermV4 nsplit: ',nsplit_thermals,' weak_inversion'
152	endif
153
154	! tests sur les valeurs negatives de l'eau
155	logexpr0=prt_level.ge.10
156	nbptspb=0
157	do k=1,klev
158	do i=1,klon
159	! Attention teste abderr 19-03-09
160	! logexpr2(i,k)=.not.q_seri(i,k).ge.0.
161	logexpr2(i,k)=.not.q_seri(i,k).ge.1.e-15
162	if (logexpr2(i,k)) then
163	q_seri(i,k)=1.e-15
164	nbptspb=nbptspb+1
165	endif
166	! if (logexpr0) &
167	! & print*,'WARN eau<0 avant therm i=',i,' k=',k &
168	! & ,' dq,q',d_q_the(i,k),q_seri(i,k)
169	enddo
170	enddo
171	if(nbptspb.GT.0) print*,'Number of points with q_seri(i,k)<=0 ',nbptspb
172
173	zdt=dtime/REAL(nsplit_thermals)
174	do isplit=1,nsplit_thermals
175
176	if (iflag_thermals.eq.1) then
177	CALL thermcell_2002(klon,klev,zdt &
178	& ,pplay,paprs,pphi &
179	& ,u_seri,v_seri,t_seri,q_seri &
180	& ,d_u_the,d_v_the,d_t_the,d_q_the &
181	& ,zfm_therm,zentr_therm &
182	& ,r_aspect_thermals,30.,w2di_thermals &
183	& ,tau_thermals)
184	else if (iflag_thermals.eq.2) then
185	CALL thermcell_sec(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 &
190	& ,r_aspect_thermals,30.,w2di_thermals &
191	& ,tau_thermals)
192	else if (iflag_thermals.eq.3) then
193	CALL thermcell(klon,klev,zdt &
194	& ,pplay,paprs,pphi &
195	& ,u_seri,v_seri,t_seri,q_seri &
196	& ,d_u_the,d_v_the,d_t_the,d_q_the &
197	& ,zfm_therm,zentr_therm &
198	& ,r_aspect_thermals,l_mix_thermals,w2di_thermals &
199	& ,tau_thermals)
200	else if (iflag_thermals.eq.10) then
201	CALL thermcell_eau(klon,klev,zdt &
202	& ,pplay,paprs,pphi &
203	& ,u_seri,v_seri,t_seri,q_seri &
204	& ,d_u_the,d_v_the,d_t_the,d_q_the &
205	& ,zfm_therm,zentr_therm &
206	& ,r_aspect_thermals,l_mix_thermals,w2di_thermals &
207	& ,tau_thermals)
208	else if (iflag_thermals.eq.11) then
209	abort_message = 'cas non prevu dans calltherm'
210	CALL abort_gcm (modname,abort_message,1)
211
212	! CALL thermcell_pluie(klon,klev,zdt &
213	! & ,pplay,paprs,pphi,zlev &
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,zqla &
217	! & ,r_aspect_thermals,l_mix_thermals,w2di_thermals &
218	! & ,tau_thermals,3)
219	else if (iflag_thermals.eq.12) then
220	CALL calcul_sec(klon,klev,zdt &
221	& ,pplay,paprs,pphi,zlev &
222	& ,u_seri,v_seri,t_seri,q_seri &
223	& ,zmax_sec,wmax_sec,zw_sec,lmix_sec &
224	& ,r_aspect_thermals,l_mix_thermals,w2di_thermals &
225	& ,tau_thermals)
226	else if (iflag_thermals==13.or.iflag_thermals==14) then
227	CALL thermcellV0_main(itap,klon,klev,zdt &
228	& ,pplay,paprs,pphi,debut &
229	& ,u_seri,v_seri,t_seri,q_seri &
230	& ,d_u_the,d_v_the,d_t_the,d_q_the &
231	& ,zfm_therm,zentr_therm,zdetr_therm,zqasc,zqla,lmax &
232	& ,ratqscth,ratqsdiff,zqsatth &
233	& ,r_aspect_thermals,l_mix_thermals &
234	& ,tau_thermals,Ale,Alp,lalim_conv,wght_th &
235	& ,zmax0,f0,zw2,fraca)
236	else if (iflag_thermals>=15.and.iflag_thermals<=18) then
237
238	! print*,'THERM iflag_thermas_ed=',iflag_thermals_ed
239	CALL thermcell_main(itap,klon,klev,zdt &
240	& ,pplay,paprs,pphi,debut &
241	& ,u_seri,v_seri,t_seri,q_seri &
242	& ,d_u_the,d_v_the,d_t_the,d_q_the &
243	& ,zfm_therm,zentr_therm,zdetr_therm,zqasc,zqla,lmax &
244	& ,ratqscth,ratqsdiff,zqsatth &
245	! & ,r_aspect_thermals,l_mix_thermals &
246	! & ,tau_thermals,iflag_thermals_ed,iflag_coupl &
247	& ,Ale,Alp,lalim_conv,wght_th &
248	& ,zmax0,f0,zw2,fraca,ztv,zpspsk &
249	& ,ztla,zthl &
250	!!! nrlmd le 10/04/2012
251	& ,pbl_tke,pctsrf,omega,airephy &
252	& ,zlcl_th,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &
253	& ,n2,s2,ale_bl_stat &
254	& ,therm_tke_max,env_tke_max &
255	& ,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &
256	& ,alp_bl_conv,alp_bl_stat &
257	!!! fin nrlmd le 10/04/2012
258	& ,ztva )
259	if (prt_level.gt.10) write(lunout,*)'Apres thermcell_main OK'
260	else
261	abort_message = 'Cas des thermiques non prevu'
262	CALL abort_gcm (modname,abort_message,1)
263	endif
264
265	! Attention : les noms sont contre intuitif.
266	! flag_bidouille_stratocu est .true. si on ne fait pas de bidouille.
267	! Il aurait mieux valu avoir un nobidouille_stratocu
268	! Et pour simplifier :
269	! nobidouille_stratocu=.not.(iflag_thermals==13.or.iflag_thermals=15)
270	! Ce serait bien de changer, mai en prenant le temps de vérifier que ca
271	! fait bien ce qu'on croit.
272
273	flag_bidouille_stratocu=iflag_thermals<=12.or.iflag_thermals==14.or.iflag_thermals==16.or.iflag_thermals==18
274
275	if (iflag_thermals<=12) then
276	lmax=1
277	do k=1,klev-1
278	zdetr_therm(:,k)=zentr_therm(:,k)+zfm_therm(:,k)-zfm_therm(:,k+1)
279	enddo
280	endif
281
282	fact(:)=0.
283	DO i=1,klon
284	logexpr1(i)=flag_bidouille_stratocu.or.weak_inversion(i).gt.0.5
285	IF(logexpr1(i)) fact(i)=1./REAL(nsplit_thermals)
286	ENDDO
287
288	DO k=1,klev
289	! transformation de la derivee en tendance
290	d_t_the(:,k)=d_t_the(:,k)dtimefact(:)
291	d_u_the(:,k)=d_u_the(:,k)dtimefact(:)
292	d_v_the(:,k)=d_v_the(:,k)dtimefact(:)
293	d_q_the(:,k)=d_q_the(:,k)dtimefact(:)
294	fm_therm(:,k)=fm_therm(:,k) &
295	& +zfm_therm(:,k)*fact(:)
296	entr_therm(:,k)=entr_therm(:,k) &
297	& +zentr_therm(:,k)*fact(:)
298	detr_therm(:,k)=detr_therm(:,k) &
299	& +zdetr_therm(:,k)*fact(:)
300	ENDDO
301	fm_therm(:,klev+1)=0.
302
303
304
305	! accumulation de la tendance
306	d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_the(:,:)
307	d_u_ajs(:,:)=d_u_ajs(:,:)+d_u_the(:,:)
308	d_v_ajs(:,:)=d_v_ajs(:,:)+d_v_the(:,:)
309	d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_the(:,:)
310
311	! incrementation des variables meteo
312	t_seri(:,:) = t_seri(:,:) + d_t_the(:,:)
313	u_seri(:,:) = u_seri(:,:) + d_u_the(:,:)
314	v_seri(:,:) = v_seri(:,:) + d_v_the(:,:)
315	qmemoire(:,:)=q_seri(:,:)
316	q_seri(:,:) = q_seri(:,:) + d_q_the(:,:)
317	if (prt_level.gt.10) write(lunout,*)'Apres apres thermcell_main OK'
318
319	DO i=1,klon
320	fm_therm(i,klev+1)=0.
321	Ale_bl(i)=Ale_bl(i)+Ale(i)/REAL(nsplit_thermals)
322	! write(22,*)'ALE CALLTHERM',Ale_bl(i),Ale(i)
323	Alp_bl(i)=Alp_bl(i)+Alp(i)/REAL(nsplit_thermals)
324	! write(23,*)'ALP CALLTHERM',Alp_bl(i),Alp(i)
325	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)
326	ENDDO
327
328	!IM 060508 marche pas comme cela !!! enddo ! isplit
329
330	! tests sur les valeurs negatives de l'eau
331	nbptspb=0
332	DO k = 1, klev
333	DO i = 1, klon
334	logexpr2(i,k)=.not.q_seri(i,k).ge.0.
335	if (logexpr2(i,k)) then
336	q_seri(i,k)=1.e-15
337	nbptspb=nbptspb+1
338	! if (prt_level.ge.10) then
339	! print*,'WARN eau<0 apres therm i=',i,' k=',k &
340	! & ,' dq,q',d_q_the(i,k),q_seri(i,k), &
341	! & 'fm=',zfm_therm(i,k),'entr=',entr_therm(i,k)
342	endif
343	ENDDO
344	ENDDO
345	IF(nbptspb.GT.0) print*,'Number of points with q_seri(i,k)<=0 ',nbptspb
346	! tests sur les valeurs de la temperature
347	nbptspb=0
348	DO k = 1, klev
349	DO i = 1, klon
350	logexpr2(i,k)=t_seri(i,k).lt.50..or.t_seri(i,k).gt.370.
351	if (logexpr2(i,k)) nbptspb=nbptspb+1
352	! if ((t_seri(i,k).lt.50.) .or. &
353	! & (t_seri(i,k).gt.370.)) then
354	! print*,'WARN temp apres therm i=',i,' k=',k &
355	! & ,' t_seri',t_seri(i,k)
356	! CALL abort
357	! endif
358	ENDDO
359	ENDDO
360	IF(nbptspb.GT.0) print*,'Number of points with q_seri(i,k)<=0 ',nbptspb
361	enddo ! isplit
362
363	!
364	!***************************************************************
365	! calcul du flux ascencant conservatif
366	! print*,'<<<<calcul flux ascendant conservatif'
367
368	fmc_therm=0.
369	do k=1,klev
370	do i=1,klon
371	if (entr_therm(i,k).gt.0.) then
372	fmc_therm(i,k+1)=fmc_therm(i,k)+entr_therm(i,k)
373	else
374	fmc_therm(i,k+1)=fmc_therm(i,k)
375	endif
376	detrc_therm(i,k)=(fmc_therm(i,k+1)-fm_therm(i,k+1)) &
377	& -(fmc_therm(i,k)-fm_therm(i,k))
378	enddo
379	enddo
380
381
382	!****************************************************************
383	! calcul de l'humidite dans l'ascendance
384	! print*,'<<<<calcul de lhumidite dans thermique'
385	!CR:on ne le calcule que pour le cas sec
386	if (iflag_thermals.le.11) then
387	do i=1,klon
388	zqasc(i,1)=q_seri(i,1)
389	do k=2,klev
390	if (fmc_therm(i,k+1).gt.1.e-6) then
391	zqasc(i,k)=(fmc_therm(i,k)*zqasc(i,k-1) &
392	& +entr_therm(i,k)*q_seri(i,k))/fmc_therm(i,k+1)
393	!CR:test on asseche le thermique
394	! zqasc(i,k)=zqasc(i,k)/2.
395	! else
396	! zqasc(i,k)=q_seri(i,k)
397	endif
398	enddo
399	enddo
400
401
402	! calcul de l'eau condensee dans l'ascendance
403	! print*,'<<<<calcul de leau condensee dans thermique'
404	do i=1,klon
405	do k=1,klev
406	clwcon0(i,k)=zqasc(i,k)-zqsat(i,k)
407	if (clwcon0(i,k).lt.0. .or. &
408	& (fm_therm(i,k+1)+detrc_therm(i,k)).lt.1.e-6) then
409	clwcon0(i,k)=0.
410	endif
411	enddo
412	enddo
413	else
414	do i=1,klon
415	do k=1,klev
416	clwcon0(i,k)=zqla(i,k)
417	if (clwcon0(i,k).lt.0. .or. &
418	& (fm_therm(i,k+1)+detrc_therm(i,k)).lt.1.e-6) then
419	clwcon0(i,k)=0.
420	endif
421	enddo
422	enddo
423	endif
424	!*******************************************************************
425
426
427	!jyg Protection contre les temperatures nulles
428	do i=1,klon
429	do k=1,klev
430	if (ztla(i,k) .lt. 1.e-10) fraca(i,k) =0.
431	enddo
432	enddo
433
434
435	return
436
437	end

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

Download in other formats:

Original Format