Context Navigation

calltherm.F90 @ 3627

Last change on this file since 3627 was 2346, checked in by Ehouarn Millour, 9 years ago

Physics/dynamics separation:

remove all references to dimensions.h from physics. nbp_lon (==iim) , nbp_lat (==jjm+1) and nbp_lev (==llm) from mod_grid_phy_lmdz should be used instead.
added module regular_lonlat_mod in phy_common to store information about the global (lon-lat) grid cell boundaries and centers.

Property copyright set to
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory
Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 15.3 KB

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

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

Download in other formats:

Original Format