Context Navigation

thermcell_main.F90 @ 1299

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

Nettoyage general pour se rapprocher des normes et éviter des erreurs a la
compilation:

tous les FLOAT() sont remplacés par des REAL()
tous les STOP dans phylmd sont remplacés par des appels à abort_gcm
le common control défini dans le fichier control.h est remplacé par le module control_mod pour éviter des messages sur l'alignement des variables dans les déclarations
des $Header$ remplacés par des $Id$ pour svn

Quelques remplacements à faire ont pu m'échapper

General cleanup of the code to try and adhere to norms and to prevent some
compilation errors:

all FLOAT() instructions have been replaced by REAL() instructions
all STOP instructions in phylmd have been replaced by calls to abort_gcm
the common block control defined in the control.h file has been replaced by the control_mod to prevent compilation warnings on the alignement of declared variables
$Header$ replaced by $Id$ for svn

Some changes which should have been made might have escaped me

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 1299 2010-01-20 14:27:21Z 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 &
11	& ,Ale_bl,Alp_bl,lalim_conv,wght_th &
12	& ,zmax0, f0,zw2,fraca)
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
49	! arguments:
50	! ----------
51
52	!IM 140508
53	INTEGER itap
54
55	INTEGER ngrid,nlay,w2di
56	real tau_thermals
57	real ptimestep,l_mix,r_aspect
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
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_sec(klon)
132	real fm0(klon,klev+1),entr0(klon,klev),detr0(klon,klev)
133	real fm(klon,klev+1),entr(klon,klev),detr(klon,klev)
134
135	real ztla(klon,klev),zqla(klon,klev),zqta(klon,klev)
136	!niveau de condensation
137	integer nivcon(klon)
138	real zcon(klon)
139	REAL CHI
140	real zcon2(klon)
141	real pcon(klon)
142	real zqsat(klon,klev)
143	real zqsatth(klon,klev)
144
145	real f_star(klon,klev+1),entr_star(klon,klev)
146	real detr_star(klon,klev)
147	real alim_star_tot(klon)
148	real alim_star(klon,klev)
149	real alim_star_clos(klon,klev)
150	real f(klon), f0(klon)
151	!FH/IM save f0
152	real zlevinter(klon)
153	logical debut
154	real seuil
155	real csc(klon,klev)
156
157	!
158	!nouvelles variables pour la convection
159	real Ale_bl(klon)
160	real Alp_bl(klon)
161	real alp_int(klon)
162	real ale_int(klon)
163	integer n_int(klon)
164	real fm_tot(klon)
165	real wght_th(klon,klev)
166	integer lalim_conv(klon)
167	!v1d logical therm
168	!v1d save therm
169
170	character*2 str2
171	character*10 str10
172
173	character (len=20) :: modname='thermcell_main'
174	character (len=80) :: abort_message
175
176	EXTERNAL SCOPY
177	!
178
179	!-----------------------------------------------------------------------
180	! initialisation:
181	! ---------------
182	!
183
184	seuil=0.25
185
186	if (debut) then
187	fm0=0.
188	entr0=0.
189	detr0=0.
190
191
192	! #define wrgrads_thermcell
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	icount=icount+1
211
212	!IM 090508 beg
213	!print*,'====================================================================='
214	!print*,'====================================================================='
215	!print*,' PAS ',icount,' PAS ',icount,' PAS ',icount,' PAS ',icount
216	!print*,'====================================================================='
217	!print*,'====================================================================='
218	!IM 090508 end
219
220	if (prt_level.ge.1) print*,'thermcell_main V4'
221
222	sorties=.true.
223	IF(ngrid.NE.klon) THEN
224	PRINT*
225	PRINT*,'STOP dans convadj'
226	PRINT*,'ngrid =',ngrid
227	PRINT*,'klon =',klon
228	ENDIF
229	!
230	write(lunout,*)'WARNING thermcell_main f0=max(f0,1.e-2)'
231	do ig=1,klon
232	if (prt_level.ge.20) then
233	print*,'th_main ig f0',ig,f0(ig)
234	endif
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	!-----------------------------------------------------------------------
241	! Calcul de T,q,ql a partir de Tl et qT dans l environnement
242	! --------------------------------------------------------------------
243	!
244	CALL thermcell_env(ngrid,nlay,po,pt,pu,pv,pplay, &
245	& pplev,zo,zh,zl,ztv,zthl,zu,zv,zpspsk,zqsat,lev_out)
246
247	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_env'
248
249	!------------------------------------------------------------------------
250	! --------------------
251	!
252	!
253	! + + + + + + + + + + +
254	!
255	!
256	! wa, fraca, wd, fracd -------------------- zlev(2), rhobarz
257	! wh,wt,wo ...
258	!
259	! + + + + + + + + + + + zh,zu,zv,zo,rho
260	!
261	!
262	! -------------------- zlev(1)
263	! \\\\\\\\\\\\\\\\\\\\
264	!
265	!
266
267	!-----------------------------------------------------------------------
268	! Calcul des altitudes des couches
269	!-----------------------------------------------------------------------
270
271	do l=2,nlay
272	zlev(:,l)=0.5*(pphi(:,l)+pphi(:,l-1))/RG
273	enddo
274	zlev(:,1)=0.
275	zlev(:,nlay+1)=(2.*pphi(:,klev)-pphi(:,klev-1))/RG
276	do l=1,nlay
277	zlay(:,l)=pphi(:,l)/RG
278	enddo
279	!calcul de l epaisseur des couches
280	do l=1,nlay
281	deltaz(:,l)=zlev(:,l+1)-zlev(:,l)
282	enddo
283
284	! print*,'2 OK convect8'
285	!-----------------------------------------------------------------------
286	! Calcul des densites
287	!-----------------------------------------------------------------------
288
289	do l=1,nlay
290	rho(:,l)=pplay(:,l)/(zpspsk(:,l)RDztv(:,l))
291	enddo
292
293	!IM
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	CALL thermcell_plume(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz, &
383	& zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot, &
384	& lalim,f0,detr_star,entr_star,f_star,csc,ztva, &
385	& ztla,zqla,zqta,zha,zw2,zw_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
386	& ,lev_out,lunout1,igout)
387	if (prt_level.ge.1) print*,'apres thermcell_plume ',lev_out
388
389	call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lalim ')
390	call test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lmix ')
391
392	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_plume'
393	if (prt_level.ge.10) then
394	write(lunout1,*) 'Dans thermcell_main 2'
395	write(lunout1,*) 'lmin ',lmin(igout)
396	write(lunout1,*) 'lalim ',lalim(igout)
397	write(lunout1,*) ' ig l alim_star entr_star detr_star f_star '
398	write(lunout1,'(i6,i4,4e15.5)') (igout,l,alim_star(igout,l),entr_star(igout,l),detr_star(igout,l) &
399	& ,f_star(igout,l+1),l=1,nint(linter(igout))+5)
400	endif
401
402	!-------------------------------------------------------------------------------
403	! Calcul des caracteristiques du thermique:zmax,zmix,wmax
404	!-------------------------------------------------------------------------------
405	!
406	CALL thermcell_height(ngrid,nlay,lalim,lmin,linter,lmix,zw2, &
407	& zlev,lmax,zmax,zmax0,zmix,wmax,lev_out)
408
409
410
411	call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lalim ')
412	call test_ltherm(ngrid,nlay,pplev,pplay,lmin ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmin ')
413	call test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmix ')
414	call test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmax ')
415
416	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_height'
417
418	!-------------------------------------------------------------------------------
419	! Fermeture,determination de f
420	!-------------------------------------------------------------------------------
421	!
422	!
423	write(lunout,*)'THERM NOUVEAU RIO2009 31B'
424	CALL thermcell_dry(ngrid,nlay,zlev,pphi,ztv,alim_star, &
425	& lalim,lmin,zmax_sec,wmax_sec,lev_out)
426
427	call test_ltherm(ngrid,nlay,pplev,pplay,lmin,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry lmin ')
428	call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry lalim ')
429
430	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_dry'
431	if (prt_level.ge.10) then
432	write(lunout1,*) 'Dans thermcell_main 1b'
433	write(lunout1,*) 'lmin ',lmin(igout)
434	write(lunout1,*) 'lalim ',lalim(igout)
435	write(lunout1,*) ' ig l alim_star entr_star detr_star f_star '
436	write(lunout1,'(i6,i4,e15.5)') (igout,l,alim_star(igout,l) &
437	& ,l=1,lalim(igout)+4)
438	endif
439
440
441
442	print*,'THERM 26JJJ'
443
444	! Choix de la fonction d'alimentation utilisee pour la fermeture.
445	! Apparemment sans importance
446	alim_star_clos(:,:)=alim_star(:,:)
447	alim_star_clos(:,:)=entr_star(:,:)+alim_star(:,:)
448
449	! Appel avec la version seche
450	CALL thermcell_closure(ngrid,nlay,r_aspect,ptimestep,rho, &
451	& zlev,lalim,alim_star_clos,f_star,zmax_sec,wmax_sec,f,lev_out)
452
453	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
454	! Appel avec les zmax et wmax tenant compte de la condensation
455	! Semble moins bien marcher
456	! CALL thermcell_closure(ngrid,nlay,r_aspect,ptimestep,rho, &
457	! & zlev,lalim,alim_star,f_star,zmax,wmax,f,lev_out)
458	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
459
460	if(prt_level.ge.1)print*,'thermcell_closure apres thermcell_closure'
461
462	if (tau_thermals>1.) then
463	lambda=exp(-ptimestep/tau_thermals)
464	f0=(1.-lambda)f+lambdaf0
465	else
466	f0=f
467	endif
468
469	! Test valable seulement en 1D mais pas genant
470	if (.not. (f0(1).ge.0.) ) then
471	abort_message = '.not. (f0(1).ge.0.)'
472	CALL abort_gcm (modname,abort_message,1)
473	endif
474
475	!-------------------------------------------------------------------------------
476	!deduction des flux
477	!-------------------------------------------------------------------------------
478
479	CALL thermcell_flux2(ngrid,nlay,ptimestep,masse, &
480	& lalim,lmax,alim_star, &
481	& entr_star,detr_star,f,rhobarz,zlev,zw2,fm,entr, &
482	& detr,zqla,lev_out,lunout1,igout)
483	!IM 060508 & detr,zqla,zmax,lev_out,lunout,igout)
484
485	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_flux'
486	call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lalim ')
487	call test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lmax ')
488
489	!------------------------------------------------------------------
490	! On ne prend pas directement les profils issus des calculs precedents
491	! mais on s'autorise genereusement une relaxation vers ceci avec
492	! une constante de temps tau_thermals (typiquement 1800s).
493	!------------------------------------------------------------------
494
495	if (tau_thermals>1.) then
496	lambda=exp(-ptimestep/tau_thermals)
497	fm0=(1.-lambda)fm+lambdafm0
498	entr0=(1.-lambda)entr+lambdaentr0
499	detr0=(1.-lambda)detr+lambdadetr0
500	else
501	fm0=fm
502	entr0=entr
503	detr0=detr
504	endif
505
506	!c------------------------------------------------------------------
507	! calcul du transport vertical
508	!------------------------------------------------------------------
509
510	call thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,masse, &
511	& zthl,zdthladj,zta,lev_out)
512	call thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,masse, &
513	& po,pdoadj,zoa,lev_out)
514
515	!------------------------------------------------------------------
516	! Calcul de la fraction de l'ascendance
517	!------------------------------------------------------------------
518	do ig=1,klon
519	fraca(ig,1)=0.
520	fraca(ig,nlay+1)=0.
521	enddo
522	do l=2,nlay
523	do ig=1,klon
524	if (zw2(ig,l).gt.1.e-10) then
525	fraca(ig,l)=fm(ig,l)/(rhobarz(ig,l)*zw2(ig,l))
526	else
527	fraca(ig,l)=0.
528	endif
529	enddo
530	enddo
531
532	!------------------------------------------------------------------
533	! calcul du transport vertical du moment horizontal
534	!------------------------------------------------------------------
535
536	!IM 090508
537	if (1.eq.1) then
538	!IM 070508 vers. _dq
539	! if (1.eq.0) then
540
541
542	! Calcul du transport de V tenant compte d'echange par gradient
543	! de pression horizontal avec l'environnement
544
545	call thermcell_dv2(ngrid,nlay,ptimestep,fm0,entr0,masse &
546	& ,fraca,zmax &
547	& ,zu,zv,pduadj,pdvadj,zua,zva,lev_out)
548
549	else
550
551	! calcul purement conservatif pour le transport de V
552	call thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,masse &
553	& ,zu,pduadj,zua,lev_out)
554	call thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,masse &
555	& ,zv,pdvadj,zva,lev_out)
556	endif
557
558	! print*,'13 OK convect8'
559	do l=1,nlay
560	do ig=1,ngrid
561	pdtadj(ig,l)=zdthladj(ig,l)*zpspsk(ig,l)
562	enddo
563	enddo
564
565	if (prt_level.ge.1) print*,'14 OK convect8'
566	!------------------------------------------------------------------
567	! Calculs de diagnostiques pour les sorties
568	!------------------------------------------------------------------
569	!calcul de fraca pour les sorties
570
571	if (sorties) then
572	if (prt_level.ge.1) print*,'14a OK convect8'
573	! calcul du niveau de condensation
574	! initialisation
575	do ig=1,ngrid
576	nivcon(ig)=0
577	zcon(ig)=0.
578	enddo
579	!nouveau calcul
580	do ig=1,ngrid
581	CHI=zh(ig,1)/(1669.0-122.0*zo(ig,1)/zqsat(ig,1)-zh(ig,1))
582	pcon(ig)=pplay(ig,1)(zo(ig,1)/zqsat(ig,1))*CHI
583	enddo
584	do k=1,nlay
585	do ig=1,ngrid
586	if ((pcon(ig).le.pplay(ig,k)) &
587	& .and.(pcon(ig).gt.pplay(ig,k+1))) then
588	zcon2(ig)=zlay(ig,k)-(pcon(ig)-pplay(ig,k))/(RG*rho(ig,k))/100.
589	endif
590	enddo
591	enddo
592	if (prt_level.ge.1) print*,'14b OK convect8'
593	do k=nlay,1,-1
594	do ig=1,ngrid
595	if (zqla(ig,k).gt.1e-10) then
596	nivcon(ig)=k
597	zcon(ig)=zlev(ig,k)
598	endif
599	enddo
600	enddo
601	if (prt_level.ge.1) print*,'14c OK convect8'
602	!calcul des moments
603	!initialisation
604	do l=1,nlay
605	do ig=1,ngrid
606	q2(ig,l)=0.
607	wth2(ig,l)=0.
608	wth3(ig,l)=0.
609	ratqscth(ig,l)=0.
610	ratqsdiff(ig,l)=0.
611	enddo
612	enddo
613	if (prt_level.ge.1) print*,'14d OK convect8'
614	if (prt_level.ge.10)write(lunout,*) &
615	& 'WARNING thermcell_main wth2=0. si zw2 > 1.e-10'
616	do l=1,nlay
617	do ig=1,ngrid
618	zf=fraca(ig,l)
619	zf2=zf/(1.-zf)
620	!
621	if (prt_level.ge.10) print*,'14e OK convect8 ig,l,zf,zf2',ig,l,zf,zf2
622	!
623	if (prt_level.ge.10) print*,'14f OK convect8 ig,l,zha zh zpspsk ',ig,l,zha(ig,l),zh(ig,l),zpspsk(ig,l)
624	thetath2(ig,l)=zf2(ztla(ig,l)-zthl(ig,l))*2
625	if(zw2(ig,l).gt.1.e-10) then
626	wth2(ig,l)=zf2(zw2(ig,l))*2
627	else
628	wth2(ig,l)=0.
629	endif
630	! print*,'wth2=',wth2(ig,l)
631	wth3(ig,l)=zf2(1-2.fraca(ig,l))/(1-fraca(ig,l)) &
632	& zw2(ig,l)zw2(ig,l)*zw2(ig,l)
633	if (prt_level.ge.10) print*,'14g OK convect8 ig,l,po',ig,l,po(ig,l)
634	q2(ig,l)=zf2(zqta(ig,l)1000.-po(ig,l)1000.)*2
635	!test: on calcul q2/po=ratqsc
636	ratqscth(ig,l)=sqrt(max(q2(ig,l),1.e-6)/(po(ig,l)*1000.))
637	enddo
638	enddo
639	!calcul des flux: q, thetal et thetav
640	do l=1,nlay
641	do ig=1,ngrid
642	wq(ig,l)=fraca(ig,l)zw2(ig,l)(zqta(ig,l)1000.-po(ig,l)1000.)
643	wthl(ig,l)=fraca(ig,l)zw2(ig,l)(ztla(ig,l)-zthl(ig,l))
644	wthv(ig,l)=fraca(ig,l)zw2(ig,l)(ztva(ig,l)-ztv(ig,l))
645	enddo
646	enddo
647	!calcul de ale_bl et alp_bl
648	!pour le calcul d'une valeur integree entre la surface et lmax
649	do ig=1,ngrid
650	alp_int(ig)=0.
651	ale_int(ig)=0.
652	n_int(ig)=0
653	enddo
654	!
655	do l=1,nlay
656	do ig=1,ngrid
657	if(l.LE.lmax(ig)) THEN
658	alp_int(ig)=alp_int(ig)+0.5rhobarz(ig,l)wth3(ig,l)
659	ale_int(ig)=ale_int(ig)+0.5zw2(ig,l)*2
660	n_int(ig)=n_int(ig)+1
661	endif
662	enddo
663	enddo
664	! print*,'avant calcul ale et alp'
665	!calcul de ALE et ALP pour la convection
666	do ig=1,ngrid
667	! Alp_bl(ig)=0.5rhobarz(ig,lmix_bis(ig))wth3(ig,lmix(ig))
668	! Alp_bl(ig)=0.5rhobarz(ig,nivcon(ig))wth3(ig,nivcon(ig))
669	! Alp_bl(ig)=0.5rhobarz(ig,lmix(ig))wth3(ig,lmix(ig))
670	! & *0.1
671	!valeur integree de alp_bl * 0.5:
672	if (n_int(ig).gt.0) then
673	Alp_bl(ig)=0.5*alp_int(ig)/n_int(ig)
674	! if (Alp_bl(ig).lt.0.) then
675	! Alp_bl(ig)=0.
676	endif
677	! endif
678	! write(18,*),'rhobarz,wth3,Alp',rhobarz(ig,nivcon(ig)),
679	! s wth3(ig,nivcon(ig)),Alp_bl(ig)
680	! write(18,*),'ALP_BL',Alp_bl(ig),lmix(ig)
681	! Ale_bl(ig)=0.5zw2(ig,lmix_bis(ig))*2
682	! if (nivcon(ig).eq.1) then
683	! Ale_bl(ig)=0.
684	! else
685	!valeur max de ale_bl:
686	Ale_bl(ig)=0.5zw2(ig,lmix(ig))*2
687	! & /2.
688	! & *0.1
689	! Ale_bl(ig)=0.5zw2(ig,lmix_bis(ig))*2
690	! if (n_int(ig).gt.0) then
691	! Ale_bl(ig)=ale_int(ig)/n_int(ig)
692	! Ale_bl(ig)=4.
693	! endif
694	! endif
695	! Ale_bl(ig)=0.5*wth2(ig,lmix_bis(ig))
696	! Ale_bl(ig)=wth2(ig,nivcon(ig))
697	! write(19,*),'wth2,ALE_BL',wth2(ig,nivcon(ig)),Ale_bl(ig)
698	enddo
699	!test:calcul de la ponderation des couches pour KE
700	!initialisations
701	! print*,'ponderation'
702	do ig=1,ngrid
703	fm_tot(ig)=0.
704	enddo
705	do ig=1,ngrid
706	do k=1,klev
707	wght_th(ig,k)=1.
708	enddo
709	enddo
710	do ig=1,ngrid
711	! lalim_conv(ig)=lmix_bis(ig)
712	!la hauteur de la couche alim_conv = hauteur couche alim_therm
713	lalim_conv(ig)=lalim(ig)
714	! zentr(ig)=zlev(ig,lalim(ig))
715	enddo
716	do ig=1,ngrid
717	do k=1,lalim_conv(ig)
718	fm_tot(ig)=fm_tot(ig)+fm(ig,k)
719	enddo
720	enddo
721	do ig=1,ngrid
722	do k=1,lalim_conv(ig)
723	if (fm_tot(ig).gt.1.e-10) then
724	! wght_th(ig,k)=fm(ig,k)/fm_tot(ig)
725	endif
726	!on pondere chaque couche par a*
727	if (alim_star(ig,k).gt.1.e-10) then
728	wght_th(ig,k)=alim_star(ig,k)
729	else
730	wght_th(ig,k)=1.
731	endif
732	enddo
733	enddo
734	! print*,'apres wght_th'
735	!test pour prolonger la convection
736	do ig=1,ngrid
737	!v1d if ((alim_star(ig,1).lt.1.e-10).and.(therm)) then
738	if ((alim_star(ig,1).lt.1.e-10)) then
739	lalim_conv(ig)=1
740	wght_th(ig,1)=1.
741	! print*,'lalim_conv ok',lalim_conv(ig),wght_th(ig,1)
742	endif
743	enddo
744
745	!calcul du ratqscdiff
746	if (prt_level.ge.1) print*,'14e OK convect8'
747	var=0.
748	vardiff=0.
749	ratqsdiff(:,:)=0.
750	do ig=1,ngrid
751	do l=1,lalim(ig)
752	var=var+alim_star(ig,l)zqta(ig,l)1000.
753	enddo
754	enddo
755	if (prt_level.ge.1) print*,'14f OK convect8'
756	do ig=1,ngrid
757	do l=1,lalim(ig)
758	zf=fraca(ig,l)
759	zf2=zf/(1.-zf)
760	vardiff=vardiff+alim_star(ig,l) &
761	& (zqta(ig,l)1000.-var)**2
762	! ratqsdiff=ratqsdiff+alim_star(ig,l)*
763	! s (zqta(ig,l)1000.-po(ig,l)1000.)**2
764	enddo
765	enddo
766	if (prt_level.ge.1) print*,'14g OK convect8'
767	do l=1,nlay
768	do ig=1,ngrid
769	ratqsdiff(ig,l)=sqrt(vardiff)/(po(ig,l)*1000.)
770	! write(11,*)'ratqsdiff=',ratqsdiff(ig,l)
771	enddo
772	enddo
773	!--------------------------------------------------------------------
774	!
775	!ecriture des fichiers sortie
776	! print*,'15 OK convect8'
777
778	#ifdef wrgrads_thermcell
779	if (prt_level.ge.1) print*,'thermcell_main sorties 3D'
780	#include "thermcell_out3d.h"
781	#endif
782
783	endif
784
785	if (prt_level.ge.1) print*,'thermcell_main FIN OK'
786
787	return
788	end
789
790	!-----------------------------------------------------------------------------
791
792	subroutine test_ltherm(klon,klev,pplev,pplay,long,seuil,ztv,po,ztva,zqla,f_star,zw2,comment)
793	IMPLICIT NONE
794	#include "iniprint.h"
795
796	integer i, k, klon,klev
797	real pplev(klon,klev+1),pplay(klon,klev)
798	real ztv(klon,klev)
799	real po(klon,klev)
800	real ztva(klon,klev)
801	real zqla(klon,klev)
802	real f_star(klon,klev)
803	real zw2(klon,klev)
804	integer long(klon)
805	real seuil
806	character*21 comment
807
808	if (prt_level.ge.1) THEN
809	print*,'WARNING !!! TEST ',comment
810	endif
811	return
812
813	! test sur la hauteur des thermiques ...
814	do i=1,klon
815	!IMtemp if (pplay(i,long(i)).lt.seuil*pplev(i,1)) then
816	if (prt_level.ge.10) then
817	print*,'WARNING ',comment,' au point ',i,' K= ',long(i)
818	print,' K P(MB) THV(K) Qenv(g/kg)THVA QLA(g/kg) F W2'
819	do k=1,klev
820	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)
821	enddo
822	endif
823	enddo
824
825
826	return
827	end
828

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

Download in other formats:

Original Format