Context Navigation

thermcell_main.F90 @ 5434

Last change on this file since 5434 was 1525, checked in by idelkadi, 14 years ago

Modifications concerning the cloud scheme:

In newmicro, it now possible to read a min and max effective radius of ice particles from physiq.def : rei_min and rei_max which were initially set to 3.5 and 61.29 microns in the code.

concerns : conf_phys.F90, nuage.h, newmicro.F

In physiq.F, in case of combination of iflag_cldcon>=5 (A. Jam cloud scheme) and iflag_coupl=5

concerns : physiq.F and thermcell_main.F90

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

File size: 27.3 KB

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

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

Download in other formats:

Original Format