Context Navigation

thermcell_main.F90 @ 1494

Last change on this file since 1494 was 1494, checked in by Laurent Fairhead, 13 years ago

Optimization of routines from the new physics

Optimisation de routines de la nouvelle physique

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

File size: 26.6 KB

Line
1	!
2	! $Id: thermcell_main.F90 1494 2011-03-09 10:05:02Z fairhead $
3	!
4	SUBROUTINE thermcell_main(itap,ngrid,nlay,ptimestep &
5	& ,pplay,pplev,pphi,debut &
6	& ,pu,pv,pt,po &
7	& ,pduadj,pdvadj,pdtadj,pdoadj &
8	& ,fm0,entr0,detr0,zqta,zqla,lmax &
9	& ,ratqscth,ratqsdiff,zqsatth &
10	& ,r_aspect,l_mix,tau_thermals,iflag_thermals_ed &
11	& ,Ale_bl,Alp_bl,lalim_conv,wght_th &
12	& ,zmax0, f0,zw2,fraca,ztv &
13	& ,zpspsk,ztla,zthl)
14
15	USE dimphy
16	USE comgeomphy , ONLY:rlond,rlatd
17	IMPLICIT NONE
18
19	!=======================================================================
20	! Auteurs: Frederic Hourdin, Catherine Rio, Anne Mathieu
21	! Version du 09.02.07
22	! Calcul du transport vertical dans la couche limite en presence
23	! de "thermiques" explicitement representes avec processus nuageux
24	!
25	! Reecriture a partir d'un listing papier a Habas, le 14/02/00
26	!
27	! le thermique est suppose homogene et dissipe par melange avec
28	! son environnement. la longueur l_mix controle l'efficacite du
29	! melange
30	!
31	! Le calcul du transport des differentes especes se fait en prenant
32	! en compte:
33	! 1. un flux de masse montant
34	! 2. un flux de masse descendant
35	! 3. un entrainement
36	! 4. un detrainement
37	!
38	!=======================================================================
39
40	!-----------------------------------------------------------------------
41	! declarations:
42	! -------------
43
44	#include "dimensions.h"
45	#include "YOMCST.h"
46	#include "YOETHF.h"
47	#include "FCTTRE.h"
48	#include "iniprint.h"
49
50	! arguments:
51	! ----------
52
53	!IM 140508
54	INTEGER itap
55
56	INTEGER ngrid,nlay,w2di
57	real tau_thermals
58	integer iflag_thermals_ed
59	real ptimestep,l_mix,r_aspect
60	REAL pt(ngrid,nlay),pdtadj(ngrid,nlay)
61	REAL pu(ngrid,nlay),pduadj(ngrid,nlay)
62	REAL pv(ngrid,nlay),pdvadj(ngrid,nlay)
63	REAL po(ngrid,nlay),pdoadj(ngrid,nlay)
64	REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
65	real pphi(ngrid,nlay)
66
67	! local:
68	! ------
69
70	integer icount
71	data icount/0/
72	save icount
73	!$OMP THREADPRIVATE(icount)
74
75	integer,save :: igout=1
76	!$OMP THREADPRIVATE(igout)
77	integer,save :: lunout1=6
78	!$OMP THREADPRIVATE(lunout1)
79	integer,save :: lev_out=10
80	!$OMP THREADPRIVATE(lev_out)
81
82	INTEGER ig,k,l,ll,ierr
83	real zsortie1d(klon)
84	INTEGER lmax(klon),lmin(klon),lalim(klon)
85	INTEGER lmix(klon)
86	INTEGER lmix_bis(klon)
87	real linter(klon)
88	real zmix(klon)
89	real zmax(klon),zw2(klon,klev+1),ztva(klon,klev),zw_est(klon,klev+1),ztva_est(klon,klev)
90	! real fraca(klon,klev)
91
92	real zmax_sec(klon)
93	!on garde le zmax du pas de temps precedent
94	real zmax0(klon)
95	!FH/IM save zmax0
96
97	real lambda
98
99	real zlev(klon,klev+1),zlay(klon,klev)
100	real deltaz(klon,klev)
101	REAL zh(klon,klev)
102	real zthl(klon,klev),zdthladj(klon,klev)
103	REAL ztv(klon,klev)
104	real zu(klon,klev),zv(klon,klev),zo(klon,klev)
105	real zl(klon,klev)
106	real zsortie(klon,klev)
107	real zva(klon,klev)
108	real zua(klon,klev)
109	real zoa(klon,klev)
110
111	real zta(klon,klev)
112	real zha(klon,klev)
113	real fraca(klon,klev+1)
114	real zf,zf2
115	real thetath2(klon,klev),wth2(klon,klev),wth3(klon,klev)
116	real q2(klon,klev)
117	! FH probleme de dimensionnement avec l'allocation dynamique
118	! common/comtherm/thetath2,wth2
119	real wq(klon,klev)
120	real wthl(klon,klev)
121	real wthv(klon,klev)
122
123	real ratqscth(klon,klev)
124	real var
125	real vardiff
126	real ratqsdiff(klon,klev)
127
128	logical sorties
129	real rho(klon,klev),rhobarz(klon,klev),masse(klon,klev)
130	real zpspsk(klon,klev)
131
132	real wmax(klon)
133	real wmax_tmp(klon)
134	real wmax_sec(klon)
135	real fm0(klon,klev+1),entr0(klon,klev),detr0(klon,klev)
136	real fm(klon,klev+1),entr(klon,klev),detr(klon,klev)
137
138	real ztla(klon,klev),zqla(klon,klev),zqta(klon,klev)
139	!niveau de condensation
140	integer nivcon(klon)
141	real zcon(klon)
142	REAL CHI
143	real zcon2(klon)
144	real pcon(klon)
145	real zqsat(klon,klev)
146	real zqsatth(klon,klev)
147
148	real f_star(klon,klev+1),entr_star(klon,klev)
149	real detr_star(klon,klev)
150	real alim_star_tot(klon)
151	real alim_star(klon,klev)
152	real alim_star_clos(klon,klev)
153	real f(klon), f0(klon)
154	!FH/IM save f0
155	real zlevinter(klon)
156	logical debut
157	real seuil
158	real csc(klon,klev)
159
160	!
161	!nouvelles variables pour la convection
162	real Ale_bl(klon)
163	real Alp_bl(klon)
164	real alp_int(klon)
165	real ale_int(klon)
166	integer n_int(klon)
167	real fm_tot(klon)
168	real wght_th(klon,klev)
169	integer lalim_conv(klon)
170	!v1d logical therm
171	!v1d save therm
172
173	character*2 str2
174	character*10 str10
175
176	character (len=20) :: modname='thermcell_main'
177	character (len=80) :: abort_message
178
179	EXTERNAL SCOPY
180	!
181
182	!-----------------------------------------------------------------------
183	! initialisation:
184	! ---------------
185	!
186
187	seuil=0.25
188
189	if (debut) then
190	fm0=0.
191	entr0=0.
192	detr0=0.
193
194
195	#undef wrgrads_thermcell
196	#ifdef wrgrads_thermcell
197	! Initialisation des sorties grads pour les thermiques.
198	! Pour l'instant en 1D sur le point igout.
199	! Utilise par thermcell_out3d.h
200	str10='therm'
201	call inigrads(1,1,rlond(igout),1.,-180.,180.,jjm, &
202	& rlatd(igout),-90.,90.,1.,llm,pplay(igout,:),1., &
203	& ptimestep,str10,'therm ')
204	#endif
205
206
207
208	endif
209
210	fm=0. ; entr=0. ; detr=0.
211
212
213	icount=icount+1
214
215	!IM 090508 beg
216	!print*,'====================================================================='
217	!print*,'====================================================================='
218	!print*,' PAS ',icount,' PAS ',icount,' PAS ',icount,' PAS ',icount
219	!print*,'====================================================================='
220	!print*,'====================================================================='
221	!IM 090508 end
222
223	if (prt_level.ge.1) print*,'thermcell_main V4'
224
225	sorties=.true.
226	IF(ngrid.NE.klon) THEN
227	PRINT*
228	PRINT*,'STOP dans convadj'
229	PRINT*,'ngrid =',ngrid
230	PRINT*,'klon =',klon
231	ENDIF
232	!
233	! write(lunout,*)'WARNING thermcell_main f0=max(f0,1.e-2)'
234	do ig=1,klon
235	f0(ig)=max(f0(ig),1.e-2)
236	zmax0(ig)=max(zmax0(ig),40.)
237	!IMmarche pas ?! if (f0(ig)<1.e-2) f0(ig)=1.e-2
238	enddo
239
240	if (prt_level.ge.20) then
241	do ig=1,ngrid
242	print*,'th_main ig f0',ig,f0(ig)
243	enddo
244	endif
245	!-----------------------------------------------------------------------
246	! Calcul de T,q,ql a partir de Tl et qT dans l environnement
247	! --------------------------------------------------------------------
248	!
249	CALL thermcell_env(ngrid,nlay,po,pt,pu,pv,pplay, &
250	& pplev,zo,zh,zl,ztv,zthl,zu,zv,zpspsk,zqsat,lev_out)
251
252	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_env'
253
254	!------------------------------------------------------------------------
255	! --------------------
256	!
257	!
258	! + + + + + + + + + + +
259	!
260	!
261	! wa, fraca, wd, fracd -------------------- zlev(2), rhobarz
262	! wh,wt,wo ...
263	!
264	! + + + + + + + + + + + zh,zu,zv,zo,rho
265	!
266	!
267	! -------------------- zlev(1)
268	! \\\\\\\\\\\\\\\\\\\\
269	!
270	!
271
272	!-----------------------------------------------------------------------
273	! Calcul des altitudes des couches
274	!-----------------------------------------------------------------------
275
276	do l=2,nlay
277	zlev(:,l)=0.5*(pphi(:,l)+pphi(:,l-1))/RG
278	enddo
279	zlev(:,1)=0.
280	zlev(:,nlay+1)=(2.*pphi(:,klev)-pphi(:,klev-1))/RG
281	do l=1,nlay
282	zlay(:,l)=pphi(:,l)/RG
283	enddo
284	!calcul de l epaisseur des couches
285	do l=1,nlay
286	deltaz(:,l)=zlev(:,l+1)-zlev(:,l)
287	enddo
288
289	! print*,'2 OK convect8'
290	!-----------------------------------------------------------------------
291	! Calcul des densites
292	!-----------------------------------------------------------------------
293
294	rho(:,:)=pplay(:,:)/(zpspsk(:,:)RDztv(:,:))
295
296	if (prt_level.ge.10)write(lunout,*) &
297	& 'WARNING thermcell_main rhobarz(:,1)=rho(:,1)'
298	rhobarz(:,1)=rho(:,1)
299
300	do l=2,nlay
301	rhobarz(:,l)=0.5*(rho(:,l)+rho(:,l-1))
302	enddo
303
304	!calcul de la masse
305	do l=1,nlay
306	masse(:,l)=(pplev(:,l)-pplev(:,l+1))/RG
307	enddo
308
309	if (prt_level.ge.1) print*,'thermcell_main apres initialisation'
310
311	!------------------------------------------------------------------
312	!
313	! /\|\
314	! -------- \| F_k+1 -------
315	! ----> D_k
316	! /\|\ <---- E_k , A_k
317	! -------- \| F_k ---------
318	! ----> D_k-1
319	! <---- E_k-1 , A_k-1
320	!
321	!
322	!
323	!
324	!
325	! ---------------------------
326	!
327	! ----- F_lmax+1=0 ---------- \
328	! lmax (zmax) \|
329	! --------------------------- \|
330	! \|
331	! --------------------------- \|
332	! \|
333	! --------------------------- \|
334	! \|
335	! --------------------------- \|
336	! \|
337	! --------------------------- \|
338	! \| E
339	! --------------------------- \| D
340	! \|
341	! --------------------------- \|
342	! \|
343	! --------------------------- \ \|
344	! lalim \| \|
345	! --------------------------- \| \|
346	! \| \|
347	! --------------------------- \| \|
348	! \| A \|
349	! --------------------------- \| \|
350	! \| \|
351	! --------------------------- \| \|
352	! lmin (=1 pour le moment) \| \|
353	! ----- F_lmin=0 ------------ / /
354	!
355	! ---------------------------
356	! //////////////////////////
357	!
358	!
359	!=============================================================================
360	! Calculs initiaux ne faisant pas intervenir les changements de phase
361	!=============================================================================
362
363	!------------------------------------------------------------------
364	! 1. alim_star est le profil vertical de l'alimentation a la base du
365	! panache thermique, calcule a partir de la flotabilite de l'air sec
366	! 2. lmin et lalim sont les indices inferieurs et superieurs de alim_star
367	!------------------------------------------------------------------
368	!
369	entr_star=0. ; detr_star=0. ; alim_star=0. ; alim_star_tot=0.
370	lmin=1
371
372	!-----------------------------------------------------------------------------
373	! 3. wmax_sec et zmax_sec sont les vitesses et altitudes maximum d'un
374	! panache sec conservatif (e=d=0) alimente selon alim_star
375	! Il s'agit d'un calcul de type CAPE
376	! zmax_sec est utilise pour determiner la geometrie du thermique.
377	!------------------------------------------------------------------------------
378	!---------------------------------------------------------------------------------
379	!calcul du melange et des variables dans le thermique
380	!--------------------------------------------------------------------------------
381	!
382	if (prt_level.ge.1) print*,'avant thermcell_plume ',lev_out
383	!IM 140508 CALL thermcell_plume(ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz, &
384
385	! Gestion temporaire de plusieurs appels à thermcell_plume au travers
386	! de la variable iflag_thermals
387
388	! print*,'THERM thermcell_main iflag_thermals_ed=',iflag_thermals_ed
389	if (iflag_thermals_ed<=9) then
390	! print*,'THERM NOUVELLE/NOUVELLE Arnaud'
391	CALL thermcell_plume(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,&
392	& zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot, &
393	& lalim,f0,detr_star,entr_star,f_star,csc,ztva, &
394	& ztla,zqla,zqta,zha,zw2,zw_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
395	& ,lev_out,lunout1,igout)
396
397	elseif (iflag_thermals_ed>9) then
398	! print*,'THERM RIO et al 2010, version d Arnaud'
399	CALL thermcellV1_plume(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,&
400	& zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot, &
401	& lalim,f0,detr_star,entr_star,f_star,csc,ztva, &
402	& ztla,zqla,zqta,zha,zw2,zw_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
403	& ,lev_out,lunout1,igout)
404
405	endif
406
407	if (prt_level.ge.1) print*,'apres thermcell_plume ',lev_out
408
409	call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lalim ')
410	call test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lmix ')
411
412	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_plume'
413	if (prt_level.ge.10) then
414	write(lunout1,*) 'Dans thermcell_main 2'
415	write(lunout1,*) 'lmin ',lmin(igout)
416	write(lunout1,*) 'lalim ',lalim(igout)
417	write(lunout1,*) ' ig l alim_star entr_star detr_star f_star '
418	write(lunout1,'(i6,i4,4e15.5)') (igout,l,alim_star(igout,l),entr_star(igout,l),detr_star(igout,l) &
419	& ,f_star(igout,l+1),l=1,nint(linter(igout))+5)
420	endif
421
422	!-------------------------------------------------------------------------------
423	! Calcul des caracteristiques du thermique:zmax,zmix,wmax
424	!-------------------------------------------------------------------------------
425	!
426	CALL thermcell_height(ngrid,nlay,lalim,lmin,linter,lmix,zw2, &
427	& zlev,lmax,zmax,zmax0,zmix,wmax,lev_out)
428	! Attention, w2 est transforme en sa racine carree dans cette routine
429	! Le probleme vient du fait que linter et lmix sont souvent égaux à 1.
430	wmax_tmp=0.
431	do l=1,nlay
432	wmax_tmp(:)=max(wmax_tmp(:),zw2(:,l))
433	enddo
434	! print*,"ZMAX ",lalim,lmin,linter,lmix,lmax,zmax,zmax0,zmix,wmax
435
436
437
438	call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lalim ')
439	call test_ltherm(ngrid,nlay,pplev,pplay,lmin ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmin ')
440	call test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmix ')
441	call test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmax ')
442
443	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_height'
444
445	!-------------------------------------------------------------------------------
446	! Fermeture,determination de f
447	!-------------------------------------------------------------------------------
448	!
449	!
450	!! write(lunout,*)'THERM NOUVEAU XXXXX'
451	CALL thermcell_dry(ngrid,nlay,zlev,pphi,ztv,alim_star, &
452	& lalim,lmin,zmax_sec,wmax_sec,lev_out)
453
454	call test_ltherm(ngrid,nlay,pplev,pplay,lmin,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry lmin ')
455	call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry lalim ')
456
457	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_dry'
458	if (prt_level.ge.10) then
459	write(lunout1,*) 'Dans thermcell_main 1b'
460	write(lunout1,*) 'lmin ',lmin(igout)
461	write(lunout1,*) 'lalim ',lalim(igout)
462	write(lunout1,*) ' ig l alim_star entr_star detr_star f_star '
463	write(lunout1,'(i6,i4,e15.5)') (igout,l,alim_star(igout,l) &
464	& ,l=1,lalim(igout)+4)
465	endif
466
467
468
469
470	! Choix de la fonction d'alimentation utilisee pour la fermeture.
471	! Apparemment sans importance
472	alim_star_clos(:,:)=alim_star(:,:)
473	alim_star_clos(:,:)=entr_star(:,:)+alim_star(:,:)
474
475	! Appel avec la version seche
476	CALL thermcell_closure(ngrid,nlay,r_aspect,ptimestep,rho, &
477	& zlev,lalim,alim_star_clos,f_star,zmax_sec,wmax_sec,f,lev_out)
478
479	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
480	! Appel avec les zmax et wmax tenant compte de la condensation
481	! Semble moins bien marcher
482	! CALL thermcell_closure(ngrid,nlay,r_aspect,ptimestep,rho, &
483	! & zlev,lalim,alim_star,f_star,zmax,wmax,f,lev_out)
484	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
485
486	if(prt_level.ge.1)print*,'thermcell_closure apres thermcell_closure'
487
488	if (tau_thermals>1.) then
489	lambda=exp(-ptimestep/tau_thermals)
490	f0=(1.-lambda)f+lambdaf0
491	else
492	f0=f
493	endif
494
495	! Test valable seulement en 1D mais pas genant
496	if (.not. (f0(1).ge.0.) ) then
497	abort_message = '.not. (f0(1).ge.0.)'
498	CALL abort_gcm (modname,abort_message,1)
499	endif
500
501	!-------------------------------------------------------------------------------
502	!deduction des flux
503	!-------------------------------------------------------------------------------
504
505	CALL thermcell_flux2(ngrid,nlay,ptimestep,masse, &
506	& lalim,lmax,alim_star, &
507	& entr_star,detr_star,f,rhobarz,zlev,zw2,fm,entr, &
508	& detr,zqla,lev_out,lunout1,igout)
509	!IM 060508 & detr,zqla,zmax,lev_out,lunout,igout)
510
511	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_flux'
512	call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lalim ')
513	call test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lmax ')
514
515	!------------------------------------------------------------------
516	! On ne prend pas directement les profils issus des calculs precedents
517	! mais on s'autorise genereusement une relaxation vers ceci avec
518	! une constante de temps tau_thermals (typiquement 1800s).
519	!------------------------------------------------------------------
520
521	if (tau_thermals>1.) then
522	lambda=exp(-ptimestep/tau_thermals)
523	fm0=(1.-lambda)fm+lambdafm0
524	entr0=(1.-lambda)entr+lambdaentr0
525	detr0=(1.-lambda)detr+lambdadetr0
526	else
527	fm0=fm
528	entr0=entr
529	detr0=detr
530	endif
531
532	!c------------------------------------------------------------------
533	! calcul du transport vertical
534	!------------------------------------------------------------------
535
536	call thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,masse, &
537	& zthl,zdthladj,zta,lev_out)
538	call thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,masse, &
539	& po,pdoadj,zoa,lev_out)
540
541	!------------------------------------------------------------------
542	! Calcul de la fraction de l'ascendance
543	!------------------------------------------------------------------
544	do ig=1,klon
545	fraca(ig,1)=0.
546	fraca(ig,nlay+1)=0.
547	enddo
548	do l=2,nlay
549	do ig=1,klon
550	if (zw2(ig,l).gt.1.e-10) then
551	fraca(ig,l)=fm(ig,l)/(rhobarz(ig,l)*zw2(ig,l))
552	else
553	fraca(ig,l)=0.
554	endif
555	enddo
556	enddo
557
558	!------------------------------------------------------------------
559	! calcul du transport vertical du moment horizontal
560	!------------------------------------------------------------------
561
562	!IM 090508
563	if (1.eq.1) then
564	!IM 070508 vers. _dq
565	! if (1.eq.0) then
566
567
568	! Calcul du transport de V tenant compte d'echange par gradient
569	! de pression horizontal avec l'environnement
570
571	call thermcell_dv2(ngrid,nlay,ptimestep,fm0,entr0,masse &
572	& ,fraca,zmax &
573	& ,zu,zv,pduadj,pdvadj,zua,zva,lev_out)
574
575	else
576
577	! calcul purement conservatif pour le transport de V
578	call thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,masse &
579	& ,zu,pduadj,zua,lev_out)
580	call thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,masse &
581	& ,zv,pdvadj,zva,lev_out)
582	endif
583
584	! print*,'13 OK convect8'
585	do l=1,nlay
586	do ig=1,ngrid
587	pdtadj(ig,l)=zdthladj(ig,l)*zpspsk(ig,l)
588	enddo
589	enddo
590
591	if (prt_level.ge.1) print*,'14 OK convect8'
592	!------------------------------------------------------------------
593	! Calculs de diagnostiques pour les sorties
594	!------------------------------------------------------------------
595	!calcul de fraca pour les sorties
596
597	if (sorties) then
598	if (prt_level.ge.1) print*,'14a OK convect8'
599	! calcul du niveau de condensation
600	! initialisation
601	do ig=1,ngrid
602	nivcon(ig)=0
603	zcon(ig)=0.
604	enddo
605	!nouveau calcul
606	do ig=1,ngrid
607	CHI=zh(ig,1)/(1669.0-122.0*zo(ig,1)/zqsat(ig,1)-zh(ig,1))
608	pcon(ig)=pplay(ig,1)(zo(ig,1)/zqsat(ig,1))*CHI
609	enddo
610	!IM do k=1,nlay
611	do k=1,nlay-1
612	do ig=1,ngrid
613	if ((pcon(ig).le.pplay(ig,k)) &
614	& .and.(pcon(ig).gt.pplay(ig,k+1))) then
615	zcon2(ig)=zlay(ig,k)-(pcon(ig)-pplay(ig,k))/(RG*rho(ig,k))/100.
616	endif
617	enddo
618	enddo
619	!IM
620	ierr=0
621	do ig=1,ngrid
622	if (pcon(ig).le.pplay(ig,nlay)) then
623	zcon2(ig)=zlay(ig,nlay)-(pcon(ig)-pplay(ig,nlay))/(RG*rho(ig,nlay))/100.
624	ierr=1
625	endif
626	enddo
627	if (ierr==1) then
628	abort_message = 'thermcellV0_main: les thermiques vont trop haut '
629	CALL abort_gcm (modname,abort_message,1)
630	endif
631
632	if (prt_level.ge.1) print*,'14b OK convect8'
633	do k=nlay,1,-1
634	do ig=1,ngrid
635	if (zqla(ig,k).gt.1e-10) then
636	nivcon(ig)=k
637	zcon(ig)=zlev(ig,k)
638	endif
639	enddo
640	enddo
641	if (prt_level.ge.1) print*,'14c OK convect8'
642	!calcul des moments
643	!initialisation
644	do l=1,nlay
645	do ig=1,ngrid
646	q2(ig,l)=0.
647	wth2(ig,l)=0.
648	wth3(ig,l)=0.
649	ratqscth(ig,l)=0.
650	ratqsdiff(ig,l)=0.
651	enddo
652	enddo
653	if (prt_level.ge.1) print*,'14d OK convect8'
654	if (prt_level.ge.10)write(lunout,*) &
655	& 'WARNING thermcell_main wth2=0. si zw2 > 1.e-10'
656	do l=1,nlay
657	do ig=1,ngrid
658	zf=fraca(ig,l)
659	zf2=zf/(1.-zf)
660	!
661	thetath2(ig,l)=zf2(ztla(ig,l)-zthl(ig,l))*2
662	if(zw2(ig,l).gt.1.e-10) then
663	wth2(ig,l)=zf2(zw2(ig,l))*2
664	else
665	wth2(ig,l)=0.
666	endif
667	wth3(ig,l)=zf2(1-2.fraca(ig,l))/(1-fraca(ig,l)) &
668	& zw2(ig,l)zw2(ig,l)*zw2(ig,l)
669	q2(ig,l)=zf2(zqta(ig,l)1000.-po(ig,l)1000.)*2
670	!test: on calcul q2/po=ratqsc
671	ratqscth(ig,l)=sqrt(max(q2(ig,l),1.e-6)/(po(ig,l)*1000.))
672	enddo
673	enddo
674	!calcul des flux: q, thetal et thetav
675	do l=1,nlay
676	do ig=1,ngrid
677	wq(ig,l)=fraca(ig,l)zw2(ig,l)(zqta(ig,l)1000.-po(ig,l)1000.)
678	wthl(ig,l)=fraca(ig,l)zw2(ig,l)(ztla(ig,l)-zthl(ig,l))
679	wthv(ig,l)=fraca(ig,l)zw2(ig,l)(ztva(ig,l)-ztv(ig,l))
680	enddo
681	enddo
682	!
683	if (prt_level.ge.10) then
684	ig=igout
685	do l=1,nlay
686	print*,'14f OK convect8 ig,l,zha zh zpspsk ',ig,l,zha(ig,l),zh(ig,l),zpspsk(ig,l)
687	print*,'14g OK convect8 ig,l,po',ig,l,po(ig,l)
688	enddo
689	endif
690
691	! print*,'avant calcul ale et alp'
692	!calcul de ALE et ALP pour la convection
693	Alp_bl(:)=0.
694	Ale_bl(:)=0.
695	! print*,'ALE,ALP ,l,zw2(ig,l),Ale_bl(ig),Alp_bl(ig)'
696	do l=1,nlay
697	do ig=1,ngrid
698	Alp_bl(ig)=max(Alp_bl(ig),0.5rhobarz(ig,l)wth3(ig,l) )
699	Ale_bl(ig)=max(Ale_bl(ig),0.5zw2(ig,l)*2)
700	! print*,'ALE,ALP',l,zw2(ig,l),Ale_bl(ig),Alp_bl(ig)
701	enddo
702	enddo
703
704	! print*,'AAAAAAA ',Alp_bl,Ale_bl,lmix
705
706
707	! TEST. IL FAUT REECRIRE LES ALE et ALP
708	! Ale_bl(:)=0.5wmax(:)wmax(:)
709	! Alp_bl(:)=0.1wmax(:)wmax(:)*wmax(:)
710
711	!test:calcul de la ponderation des couches pour KE
712	!initialisations
713	! print*,'ponderation'
714
715	fm_tot(:)=0.
716	wght_th(:,:)=1.
717	lalim_conv(:)=lalim(:)
718
719	do k=1,klev
720	do ig=1,ngrid
721	if (k<=lalim_conv(ig)) fm_tot(ig)=fm_tot(ig)+fm(ig,k)
722	enddo
723	enddo
724
725	! assez bizarre car, si on est dans la couche d'alim et que alim_star et
726	! plus petit que 1.e-10, on prend wght_th=1.
727	do k=1,klev
728	do ig=1,ngrid
729	if (k<=lalim_conv(ig).and.alim_star(ig,k)>1.e-10) then
730	wght_th(ig,k)=alim_star(ig,k)
731	endif
732	enddo
733	enddo
734
735	! print*,'apres wght_th'
736	!test pour prolonger la convection
737	do ig=1,ngrid
738	!v1d if ((alim_star(ig,1).lt.1.e-10).and.(therm)) then
739	if ((alim_star(ig,1).lt.1.e-10)) then
740	lalim_conv(ig)=1
741	wght_th(ig,1)=1.
742	! print*,'lalim_conv ok',lalim_conv(ig),wght_th(ig,1)
743	endif
744	enddo
745
746
747	!calcul du ratqscdiff
748	if (prt_level.ge.1) print*,'14e OK convect8'
749	var=0.
750	vardiff=0.
751	ratqsdiff(:,:)=0.
752
753	do l=1,klev
754	do ig=1,ngrid
755	if (l<=lalim(ig)) then
756	var=var+alim_star(ig,l)zqta(ig,l)1000.
757	endif
758	enddo
759	enddo
760
761	if (prt_level.ge.1) print*,'14f OK convect8'
762
763	do l=1,klev
764	do ig=1,ngrid
765	if (l<=lalim(ig)) then
766	zf=fraca(ig,l)
767	zf2=zf/(1.-zf)
768	vardiff=vardiff+alim_star(ig,l)(zqta(ig,l)1000.-var)**2
769	endif
770	enddo
771	enddo
772
773	if (prt_level.ge.1) print*,'14g OK convect8'
774	do l=1,nlay
775	do ig=1,ngrid
776	ratqsdiff(ig,l)=sqrt(vardiff)/(po(ig,l)*1000.)
777	! write(11,*)'ratqsdiff=',ratqsdiff(ig,l)
778	enddo
779	enddo
780	!--------------------------------------------------------------------
781	!
782	!ecriture des fichiers sortie
783	! print*,'15 OK convect8 CCCCCCCCCCCCCCCCCCc'
784
785	#ifdef wrgrads_thermcell
786	if (prt_level.ge.1) print*,'thermcell_main sorties 3D'
787	#include "thermcell_out3d.h"
788	#endif
789
790	endif
791
792	if (prt_level.ge.1) print*,'thermcell_main FIN OK'
793
794	return
795	end
796
797	!-----------------------------------------------------------------------------
798
799	subroutine test_ltherm(klon,klev,pplev,pplay,long,seuil,ztv,po,ztva,zqla,f_star,zw2,comment)
800	IMPLICIT NONE
801	#include "iniprint.h"
802
803	integer i, k, klon,klev
804	real pplev(klon,klev+1),pplay(klon,klev)
805	real ztv(klon,klev)
806	real po(klon,klev)
807	real ztva(klon,klev)
808	real zqla(klon,klev)
809	real f_star(klon,klev)
810	real zw2(klon,klev)
811	integer long(klon)
812	real seuil
813	character*21 comment
814
815	if (prt_level.ge.1) THEN
816	print*,'WARNING !!! TEST ',comment
817	endif
818	return
819
820	! test sur la hauteur des thermiques ...
821	do i=1,klon
822	!IMtemp if (pplay(i,long(i)).lt.seuil*pplev(i,1)) then
823	if (prt_level.ge.10) then
824	print*,'WARNING ',comment,' au point ',i,' K= ',long(i)
825	print,' K P(MB) THV(K) Qenv(g/kg)THVA QLA(g/kg) F W2'
826	do k=1,klev
827	write(6,'(i3,7f10.3)') k,pplay(i,k),ztv(i,k),1000po(i,k),ztva(i,k),1000zqla(i,k),f_star(i,k),zw2(i,k)
828	enddo
829	endif
830	enddo
831
832
833	return
834	end
835

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

Download in other formats:

Original Format