Context Navigation

thermcell_main.F90 @ 950

Last change on this file since 950 was 940, checked in by Laurent Fairhead, 17 years ago
On remplace le fichier include dimphy.h par le module dimphy.F90i pour etre coherent avec le partout LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 23.2 KB

Rev	Line
[878]	1	!
	2	! $Header$
	3	!
	4	SUBROUTINE thermcell_main(ngrid,nlay,ptimestep &
	5	& ,pplay,pplev,pphi,debut &
	6	& ,pu,pv,pt,po &
	7	& ,pduadj,pdvadj,pdtadj,pdoadj &
	8	& ,fm0,entr0,zqla,lmax &
	9	& ,ratqscth,ratqsdiff,zqsatth &
[879]	10	& ,r_aspect,l_mix,w2di,tho &
[927]	11	& ,Ale_bl,Alp_bl,lalim_conv,wght_th &
	12	& ,zmax0, f0)
[878]	13
[940]	14	use dimphy
[878]	15	IMPLICIT NONE
	16
	17	!=======================================================================
	18	! Auteurs: Frederic Hourdin, Catherine Rio, Anne Mathieu
	19	! Version du 09.02.07
	20	! Calcul du transport vertical dans la couche limite en presence
	21	! de "thermiques" explicitement representes avec processus nuageux
	22	!
	23	! Réécriture à partir d'un listing papier à Habas, le 14/02/00
	24	!
	25	! le thermique est supposé homogène et dissipé par mélange avec
	26	! son environnement. la longueur l_mix contrôle l'efficacité du
	27	! mélange
	28	!
	29	! Le calcul du transport des différentes espèces se fait en prenant
	30	! en compte:
	31	! 1. un flux de masse montant
	32	! 2. un flux de masse descendant
	33	! 3. un entrainement
	34	! 4. un detrainement
	35	!
	36	!=======================================================================
	37
	38	!-----------------------------------------------------------------------
	39	! declarations:
	40	! -------------
	41
	42	#include "dimensions.h"
[940]	43	!#include "dimphy.h"
[878]	44	#include "YOMCST.h"
	45	#include "YOETHF.h"
	46	#include "FCTTRE.h"
[938]	47	#include "iniprint.h"
[878]	48
	49	! arguments:
	50	! ----------
	51
	52	INTEGER ngrid,nlay,w2di,tho
	53	real ptimestep,l_mix,r_aspect
	54	REAL pt(ngrid,nlay),pdtadj(ngrid,nlay)
	55	REAL pu(ngrid,nlay),pduadj(ngrid,nlay)
	56	REAL pv(ngrid,nlay),pdvadj(ngrid,nlay)
	57	REAL po(ngrid,nlay),pdoadj(ngrid,nlay)
	58	REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
	59	real pphi(ngrid,nlay)
	60
	61	! local:
	62	! ------
	63
[883]	64	integer,save :: igout=1
[938]	65	integer,save :: lunout1=6
[883]	66	integer,save :: lev_out=10
[878]	67
	68	INTEGER ig,k,l,ll
	69	real zsortie1d(klon)
	70	INTEGER lmax(klon),lmin(klon),lalim(klon)
	71	INTEGER lmix(klon)
	72	real linter(klon)
	73	real zmix(klon)
	74	real zmax(klon),zw2(klon,klev+1),ztva(klon,klev)
	75	real zmax_sec(klon)
	76	real w_est(klon,klev+1)
	77	!on garde le zmax du pas de temps precedent
	78	real zmax0(klon)
[927]	79	!FH/IM save zmax0
[878]	80
	81	real zlev(klon,klev+1),zlay(klon,klev)
	82	real deltaz(klon,klev)
	83	REAL zh(klon,klev),zdhadj(klon,klev)
	84	real zthl(klon,klev),zdthladj(klon,klev)
	85	REAL ztv(klon,klev)
	86	real zu(klon,klev),zv(klon,klev),zo(klon,klev)
	87	real zl(klon,klev)
	88	real zsortie(klon,klev)
	89	real zva(klon,klev)
	90	real zua(klon,klev)
	91	real zoa(klon,klev)
	92
	93	real zta(klon,klev)
	94	real zha(klon,klev)
	95	real fraca(klon,klev+1)
	96	real zf,zf2
	97	real thetath2(klon,klev),wth2(klon,klev),wth3(klon,klev)
	98	real q2(klon,klev)
[940]	99	! common/comtherm/thetath2,wth2
[878]	100
	101	real ratqscth(klon,klev)
	102	real var
	103	real vardiff
	104	real ratqsdiff(klon,klev)
	105	integer isplit,nsplit
	106	parameter (nsplit=10)
	107	data isplit/0/
	108	save isplit
	109
	110	logical sorties
	111	real rho(klon,klev),rhobarz(klon,klev+1),masse(klon,klev)
	112	real zpspsk(klon,klev)
	113
	114	real wmax(klon)
	115	real wmax_sec(klon)
	116	real fm0(klon,klev+1),entr0(klon,klev),detr(klon,klev)
	117	real detr0(klon,klev)
	118	real fm(klon,klev+1),entr(klon,klev)
	119
	120	real ztla(klon,klev),zqla(klon,klev),zqta(klon,klev)
	121	!niveau de condensation
[879]	122	integer nivcon(klon)
[878]	123	real zcon(klon)
	124	REAL CHI
	125	real zcon2(klon)
	126	real pcon(klon)
	127	real zqsat(klon,klev)
	128	real zqsatth(klon,klev)
	129
	130	real f_star(klon,klev+1),entr_star(klon,klev)
	131	real detr_star(klon,klev)
	132	real alim_star_tot(klon),alim_star2(klon)
	133	real alim_star(klon,klev)
	134	real f(klon), f0(klon)
[927]	135	!FH/IM save f0
[878]	136	real zlevinter(klon)
	137	logical debut
	138	real seuil
	139
	140	!
[879]	141	!nouvelles variables pour la convection
	142	real Ale_bl(klon)
	143	real Alp_bl(klon)
	144	real alp_int(klon)
	145	real ale_int(klon)
	146	integer n_int(klon)
	147	real fm_tot(klon)
	148	real wght_th(klon,klev)
	149	integer lalim_conv(klon)
[926]	150	!v1d logical therm
	151	!v1d save therm
[878]	152
	153	character*2 str2
	154	character*10 str10
	155
	156	EXTERNAL SCOPY
	157	!
	158
	159	!-----------------------------------------------------------------------
	160	! initialisation:
	161	! ---------------
	162	!
	163
	164	seuil=0.25
	165
[938]	166	if (prt_level.ge.1) print*,'thermcell_main V4'
[878]	167
	168	sorties=.true.
	169	IF(ngrid.NE.klon) THEN
	170	PRINT*
	171	PRINT*,'STOP dans convadj'
	172	PRINT*,'ngrid =',ngrid
	173	PRINT*,'klon =',klon
	174	ENDIF
	175	!
	176	!Initialisation
	177	!
	178	do ig=1,klon
[927]	179	!FH/IM 130308 if ((debut).or.((.not.debut).and.(f0(ig).lt.1.e-10))) then
	180	if ((.not.debut).and.(f0(ig).lt.1.e-10)) then
[878]	181	f0(ig)=1.e-5
	182	zmax0(ig)=40.
[926]	183	!v1d therm=.false.
[878]	184	endif
	185	enddo
	186
	187
	188	!-----------------------------------------------------------------------
	189	! Calcul de T,q,ql a partir de Tl et qT dans l environnement
	190	! --------------------------------------------------------------------
	191	!
	192	CALL thermcell_env(ngrid,nlay,po,pt,pu,pv,pplay, &
	193	& pplev,zo,zh,zl,ztv,zthl,zu,zv,zpspsk,zqsat,lev_out)
	194
[938]	195	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_env'
[878]	196
	197	!------------------------------------------------------------------------
	198	! --------------------
	199	!
	200	!
	201	! + + + + + + + + + + +
	202	!
	203	!
	204	! wa, fraca, wd, fracd -------------------- zlev(2), rhobarz
	205	! wh,wt,wo ...
	206	!
	207	! + + + + + + + + + + + zh,zu,zv,zo,rho
	208	!
	209	!
	210	! -------------------- zlev(1)
	211	! \\\\\\\\\\\\\\\\\\\\
	212	!
	213	!
	214
	215	!-----------------------------------------------------------------------
	216	! Calcul des altitudes des couches
	217	!-----------------------------------------------------------------------
	218
	219	do l=2,nlay
	220	zlev(:,l)=0.5*(pphi(:,l)+pphi(:,l-1))/RG
	221	enddo
	222	zlev(:,1)=0.
	223	zlev(:,nlay+1)=(2.*pphi(:,klev)-pphi(:,klev-1))/RG
	224	do l=1,nlay
	225	zlay(:,l)=pphi(:,l)/RG
	226	enddo
	227	!calcul de l epaisseur des couches
	228	do l=1,nlay
	229	deltaz(:,l)=zlev(:,l+1)-zlev(:,l)
	230	enddo
	231
	232	! print*,'2 OK convect8'
	233	!-----------------------------------------------------------------------
	234	! Calcul des densites
	235	!-----------------------------------------------------------------------
	236
	237	do l=1,nlay
	238	rho(:,l)=pplay(:,l)/(zpspsk(:,l)RDztv(:,l))
	239	enddo
	240
	241	do l=2,nlay
	242	rhobarz(:,l)=0.5*(rho(:,l)+rho(:,l-1))
	243	enddo
	244
	245	!calcul de la masse
	246	do l=1,nlay
	247	masse(:,l)=(pplev(:,l)-pplev(:,l+1))/RG
	248	enddo
	249
[938]	250	if (prt_level.ge.1) print*,'thermcell_main apres initialisation'
[878]	251
	252	!------------------------------------------------------------------
	253	!
	254	! /\|\
	255	! -------- \| F_k+1 -------
	256	! ----> D_k
	257	! /\|\ <---- E_k , A_k
	258	! -------- \| F_k ---------
	259	! ----> D_k-1
	260	! <---- E_k-1 , A_k-1
	261	!
	262	!
	263	!
	264	!
	265	!
	266	! ---------------------------
	267	!
	268	! ----- F_lmax+1=0 ---------- \
	269	! lmax (zmax) \|
	270	! --------------------------- \|
	271	! \|
	272	! --------------------------- \|
	273	! \|
	274	! --------------------------- \|
	275	! \|
	276	! --------------------------- \|
	277	! \|
	278	! --------------------------- \|
	279	! \| E
	280	! --------------------------- \| D
	281	! \|
	282	! --------------------------- \|
	283	! \|
	284	! --------------------------- \ \|
	285	! lalim \| \|
	286	! --------------------------- \| \|
	287	! \| \|
	288	! --------------------------- \| \|
	289	! \| A \|
	290	! --------------------------- \| \|
	291	! \| \|
	292	! --------------------------- \| \|
	293	! lmin (=1 pour le moment) \| \|
	294	! ----- F_lmin=0 ------------ / /
	295	!
	296	! ---------------------------
	297	! //////////////////////////
	298	!
	299	!
	300	!=============================================================================
	301	! Calculs initiaux ne faisant pas intervenir les changements de phase
	302	!=============================================================================
	303
	304	!------------------------------------------------------------------
	305	! 1. alim_star est le profil vertical de l'alimentation à la base du
	306	! panache thermique, calculé à partir de la flotabilité de l'air sec
	307	! 2. lmin et lalim sont les indices inferieurs et superieurs de alim_star
	308	!------------------------------------------------------------------
	309	!
	310	entr_star=0. ; detr_star=0. ; alim_star=0. ; alim_star_tot=0.
	311	CALL thermcell_init(ngrid,nlay,ztv,zlev, &
	312	& lalim,lmin,alim_star,alim_star_tot,lev_out)
	313
	314	call test_ltherm(ngrid,nlay,pplev,pplay,lmin,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_init lmin ')
	315	call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_init lalim ')
	316
	317
[938]	318	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_init'
	319	if (prt_level.ge.10) then
[878]	320	write(lunout,*) 'Dans thermcell_main 1'
	321	write(lunout,*) 'lmin ',lmin(igout)
	322	write(lunout,*) 'lalim ',lalim(igout)
	323	write(lunout,*) ' ig l alim_star thetav'
	324	write(lunout,'(i6,i4,2e15.5)') (igout,l,alim_star(igout,l) &
	325	& ,ztv(igout,l),l=1,lalim(igout)+4)
	326	endif
	327
[938]	328	!v1d do ig=1,klon
[926]	329	!v1d if (alim_star(ig,1).gt.1.e-10) then
	330	!v1d therm=.true.
	331	!v1d endif
[938]	332	!v1d enddo
[878]	333	!-----------------------------------------------------------------------------
	334	! 3. wmax_sec et zmax_sec sont les vitesses et altitudes maximum d'un
	335	! panache sec conservatif (e=d=0) alimente selon alim_star
	336	! Il s'agit d'un calcul de type CAPE
	337	! zmax_sec est utilisé pour déterminer la géométrie du thermique.
	338	!------------------------------------------------------------------------------
	339	!
	340	CALL thermcell_dry(ngrid,nlay,zlev,pphi,ztv,alim_star, &
	341	& lalim,lmin,zmax_sec,wmax_sec,lev_out)
	342
	343	call test_ltherm(ngrid,nlay,pplev,pplay,lmin,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry lmin ')
	344	call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry lalim ')
	345
[938]	346	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_dry'
	347	if (prt_level.ge.10) then
[878]	348	write(lunout,*) 'Dans thermcell_main 1b'
	349	write(lunout,*) 'lmin ',lmin(igout)
	350	write(lunout,*) 'lalim ',lalim(igout)
	351	write(lunout,*) ' ig l alim_star entr_star detr_star f_star '
	352	write(lunout,'(i6,i4,e15.5)') (igout,l,alim_star(igout,l) &
	353	& ,l=1,lalim(igout)+4)
	354	endif
	355
	356
	357
	358	!---------------------------------------------------------------------------------
	359	!calcul du melange et des variables dans le thermique
	360	!--------------------------------------------------------------------------------
	361	!
	362	CALL thermcell_plume(ngrid,nlay,ztv,zthl,po,zl,rhobarz, &
	363	& zlev,pplev,pphi,zpspsk,l_mix,r_aspect,alim_star, &
	364	& lalim,zmax_sec,f0,detr_star,entr_star,f_star,ztva, &
	365	& ztla,zqla,zqta,zha,zw2,zqsatth,lmix,linter,lev_out)
	366	call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lalim ')
	367	call test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lmix ')
	368
[938]	369	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_plume'
	370	if (prt_level.ge.10) then
[878]	371	write(lunout,*) 'Dans thermcell_main 2'
	372	write(lunout,*) 'lmin ',lmin(igout)
	373	write(lunout,*) 'lalim ',lalim(igout)
	374	write(lunout,*) ' ig l alim_star entr_star detr_star f_star '
	375	write(lunout,'(i6,i4,4e15.5)') (igout,l,alim_star(igout,l),entr_star(igout,l),detr_star(igout,l) &
	376	& ,f_star(igout,l+1),l=1,nint(linter(igout))+5)
	377	endif
	378
	379	!-------------------------------------------------------------------------------
	380	! Calcul des caracteristiques du thermique:zmax,zmix,wmax
	381	!-------------------------------------------------------------------------------
	382	!
	383	CALL thermcell_height(ngrid,nlay,lalim,lmin,linter,lmix,zw2, &
	384	& zlev,lmax,zmax,zmax0,zmix,wmax,lev_out)
	385
	386
	387	call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lalim ')
	388	call test_ltherm(ngrid,nlay,pplev,pplay,lmin ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmin ')
	389	call test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmix ')
	390	call test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmax ')
	391
[938]	392	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_height'
[878]	393
	394	!-------------------------------------------------------------------------------
	395	! Fermeture,determination de f
	396	!-------------------------------------------------------------------------------
	397
	398	CALL thermcell_closure(ngrid,nlay,r_aspect,ptimestep,rho, &
	399	& zlev,lalim,alim_star,zmax_sec,wmax_sec,zmax,wmax,f,f0,lev_out)
	400
[938]	401	if(prt_level.ge.1)print*,'thermcell_closure apres thermcell_closure'
[878]	402
	403	!-------------------------------------------------------------------------------
	404	!deduction des flux
	405	!-------------------------------------------------------------------------------
	406
	407	CALL thermcell_flux(ngrid,nlay,ptimestep,masse, &
	408	& lalim,lmax,alim_star, &
	409	& entr_star,detr_star,f,rhobarz,zlev,zw2,fm,entr, &
	410	& detr,zqla,zmax,lev_out,lunout,igout)
	411
[938]	412	if (prt_level.ge.1) print*,'thermcell_main apres thermcell_flux'
[878]	413	call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lalim ')
	414	call test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lmax ')
	415
	416	!c------------------------------------------------------------------
	417	! calcul du transport vertical
	418	!------------------------------------------------------------------
	419
	420	if (w2di.eq.1) then
	421	fm0=fm0+ptimestep*(fm-fm0)/float(tho)
	422	entr0=entr0+ptimestep*(entr-entr0)/float(tho)
	423	else
	424	fm0=fm
	425	entr0=entr
	426	detr0=detr
	427	endif
	428
	429	call thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,masse, &
	430	& zthl,zdthladj,zta,lev_out)
	431	call thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,masse, &
	432	& po,pdoadj,zoa,lev_out)
	433
[883]	434	!------------------------------------------------------------------
	435	! Calcul de la fraction de l'ascendance
	436	!------------------------------------------------------------------
	437	do ig=1,klon
	438	fraca(ig,1)=0.
	439	fraca(ig,nlay+1)=0.
	440	enddo
	441	do l=2,nlay
	442	do ig=1,klon
	443	if (zw2(ig,l).gt.1.e-10) then
	444	fraca(ig,l)=fm(ig,l)/(rhobarz(ig,l)*zw2(ig,l))
	445	else
	446	fraca(ig,l)=0.
	447	endif
	448	enddo
	449	enddo
	450
	451	!------------------------------------------------------------------
	452	! calcul du transport vertical du moment horizontal
	453	!------------------------------------------------------------------
[878]	454
[883]	455	if (1.eq.1) then
	456
	457
[878]	458	! Calcul du transport de V tenant compte d'echange par gradient
	459	! de pression horizontal avec l'environnement
	460
	461	call thermcell_dv2(ngrid,nlay,ptimestep,fm0,entr0,masse &
	462	& ,fraca,zmax &
[883]	463	& ,zu,zv,pduadj,pdvadj,zua,zva,igout,lev_out)
[878]	464	else
	465
	466	! calcul purement conservatif pour le transport de V
	467	call thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,masse &
	468	& ,zu,pduadj,zua,lev_out)
	469	call thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,masse &
	470	& ,zv,pdvadj,zva,lev_out)
	471	endif
	472
	473	! print*,'13 OK convect8'
	474	do l=1,nlay
	475	do ig=1,ngrid
	476	pdtadj(ig,l)=zdthladj(ig,l)*zpspsk(ig,l)
	477	enddo
	478	enddo
	479
[938]	480	! print*,'14 OK convect8'
[878]	481	!------------------------------------------------------------------
	482	! Calculs de diagnostiques pour les sorties
	483	!------------------------------------------------------------------
	484	!calcul de fraca pour les sorties
	485
	486	if (sorties) then
[938]	487	! print*,'14a OK convect8'
[878]	488	! calcul du niveau de condensation
	489	! initialisation
	490	do ig=1,ngrid
[879]	491	nivcon(ig)=0
[878]	492	zcon(ig)=0.
	493	enddo
	494	!nouveau calcul
	495	do ig=1,ngrid
	496	CHI=zh(ig,1)/(1669.0-122.0*zo(ig,1)/zqsat(ig,1)-zh(ig,1))
	497	pcon(ig)=pplay(ig,1)(zo(ig,1)/zqsat(ig,1))*CHI
	498	enddo
	499	do k=1,nlay
	500	do ig=1,ngrid
	501	if ((pcon(ig).le.pplay(ig,k)) &
	502	& .and.(pcon(ig).gt.pplay(ig,k+1))) then
	503	zcon2(ig)=zlay(ig,k)-(pcon(ig)-pplay(ig,k))/(RG*rho(ig,k))/100.
	504	endif
	505	enddo
	506	enddo
[938]	507	! print*,'14b OK convect8'
[878]	508	do k=nlay,1,-1
	509	do ig=1,ngrid
	510	if (zqla(ig,k).gt.1e-10) then
	511	nivcon(ig)=k
	512	zcon(ig)=zlev(ig,k)
	513	endif
	514	enddo
	515	enddo
[938]	516	! print*,'14c OK convect8'
[878]	517	!calcul des moments
	518	!initialisation
	519	do l=1,nlay
	520	do ig=1,ngrid
	521	q2(ig,l)=0.
	522	wth2(ig,l)=0.
	523	wth3(ig,l)=0.
	524	ratqscth(ig,l)=0.
	525	ratqsdiff(ig,l)=0.
	526	enddo
	527	enddo
[938]	528	! print*,'14d OK convect8'
[878]	529	do l=1,nlay
	530	do ig=1,ngrid
	531	zf=fraca(ig,l)
	532	zf2=zf/(1.-zf)
	533	thetath2(ig,l)=zf2(zha(ig,l)-zh(ig,l)/zpspsk(ig,l))*2
	534	wth2(ig,l)=zf2(zw2(ig,l))*2
	535	! print*,'wth2=',wth2(ig,l)
	536	wth3(ig,l)=zf2(1-2.fraca(ig,l))/(1-fraca(ig,l)) &
	537	& zw2(ig,l)zw2(ig,l)*zw2(ig,l)
	538	q2(ig,l)=zf2(zqta(ig,l)1000.-po(ig,l)1000.)*2
	539	!test: on calcul q2/po=ratqsc
	540	ratqscth(ig,l)=sqrt(max(q2(ig,l),1.e-6)/(po(ig,l)*1000.))
	541	enddo
	542	enddo
[879]	543	!calcul de ale_bl et alp_bl
	544	!pour le calcul d'une valeur intégrée entre la surface et lmax
	545	do ig=1,ngrid
	546	alp_int(ig)=0.
	547	ale_int(ig)=0.
	548	n_int(ig)=0
	549	enddo
	550	do ig=1,ngrid
	551	! do l=nivcon(ig),lmax(ig)
	552	do l=1,lmax(ig)
	553	alp_int(ig)=alp_int(ig)+0.5rhobarz(ig,l)wth3(ig,l)
	554	ale_int(ig)=ale_int(ig)+0.5zw2(ig,l)*2
	555	n_int(ig)=n_int(ig)+1
	556	enddo
	557	enddo
	558	! print*,'avant calcul ale et alp'
	559	!calcul de ALE et ALP pour la convection
	560	do ig=1,ngrid
	561	! Alp_bl(ig)=0.5rhobarz(ig,lmix_bis(ig))wth3(ig,lmix(ig))
	562	! Alp_bl(ig)=0.5rhobarz(ig,nivcon(ig))wth3(ig,nivcon(ig))
	563	! Alp_bl(ig)=0.5rhobarz(ig,lmix(ig))wth3(ig,lmix(ig))
	564	! & *0.1
	565	!valeur integree de alp_bl * 0.5:
	566	if (n_int(ig).gt.0) then
	567	Alp_bl(ig)=0.5*alp_int(ig)/n_int(ig)
	568	! if (Alp_bl(ig).lt.0.) then
	569	! Alp_bl(ig)=0.
	570	endif
	571	! endif
	572	! write(18,*),'rhobarz,wth3,Alp',rhobarz(ig,nivcon(ig)),
	573	! s wth3(ig,nivcon(ig)),Alp_bl(ig)
	574	! write(18,*),'ALP_BL',Alp_bl(ig),lmix(ig)
	575	! Ale_bl(ig)=0.5zw2(ig,lmix_bis(ig))*2
	576	! if (nivcon(ig).eq.1) then
	577	! Ale_bl(ig)=0.
	578	! else
	579	!valeur max de ale_bl:
	580	Ale_bl(ig)=0.5zw2(ig,lmix(ig))*2
	581	! & /2.
	582	! & *0.1
	583	! Ale_bl(ig)=0.5zw2(ig,lmix_bis(ig))*2
	584	! if (n_int(ig).gt.0) then
	585	! Ale_bl(ig)=ale_int(ig)/n_int(ig)
	586	! Ale_bl(ig)=4.
	587	! endif
	588	! endif
	589	! Ale_bl(ig)=0.5*wth2(ig,lmix_bis(ig))
	590	! Ale_bl(ig)=wth2(ig,nivcon(ig))
	591	! write(19,*),'wth2,ALE_BL',wth2(ig,nivcon(ig)),Ale_bl(ig)
	592	enddo
	593	!test:calcul de la ponderation des couches pour KE
	594	!initialisations
	595	! print*,'ponderation'
	596	do ig=1,ngrid
	597	fm_tot(ig)=0.
	598	enddo
	599	do ig=1,ngrid
	600	do k=1,klev
	601	wght_th(ig,k)=1.
	602	enddo
	603	enddo
	604	do ig=1,ngrid
	605	! lalim_conv(ig)=lmix_bis(ig)
	606	!la hauteur de la couche alim_conv = hauteur couche alim_therm
	607	lalim_conv(ig)=lalim(ig)
	608	! zentr(ig)=zlev(ig,lalim(ig))
	609	enddo
	610	do ig=1,ngrid
	611	do k=1,lalim_conv(ig)
	612	fm_tot(ig)=fm_tot(ig)+fm(ig,k)
	613	enddo
	614	enddo
	615	do ig=1,ngrid
	616	do k=1,lalim_conv(ig)
	617	if (fm_tot(ig).gt.1.e-10) then
	618	! wght_th(ig,k)=fm(ig,k)/fm_tot(ig)
	619	endif
	620	!on pondere chaque couche par a*
	621	if (alim_star(ig,k).gt.1.e-10) then
	622	wght_th(ig,k)=alim_star(ig,k)
	623	else
	624	wght_th(ig,k)=1.
	625	endif
	626	enddo
	627	enddo
	628	! print*,'apres wght_th'
	629	!test pour prolonger la convection
	630	do ig=1,ngrid
[926]	631	!v1d if ((alim_star(ig,1).lt.1.e-10).and.(therm)) then
	632	if ((alim_star(ig,1).lt.1.e-10)) then
[879]	633	lalim_conv(ig)=1
	634	wght_th(ig,1)=1.
	635	! print*,'lalim_conv ok',lalim_conv(ig),wght_th(ig,1)
	636	endif
	637	enddo
	638
[878]	639	!calcul du ratqscdiff
[938]	640	! print*,'14e OK convect8'
[878]	641	var=0.
	642	vardiff=0.
	643	ratqsdiff(:,:)=0.
	644	do ig=1,ngrid
	645	do l=1,lalim(ig)
	646	var=var+alim_star(ig,l)zqta(ig,l)1000.
	647	enddo
	648	enddo
[938]	649	! print*,'14f OK convect8'
[878]	650	do ig=1,ngrid
	651	do l=1,lalim(ig)
	652	zf=fraca(ig,l)
	653	zf2=zf/(1.-zf)
	654	vardiff=vardiff+alim_star(ig,l) &
	655	& (zqta(ig,l)1000.-var)**2
	656	! ratqsdiff=ratqsdiff+alim_star(ig,l)*
	657	! s (zqta(ig,l)1000.-po(ig,l)1000.)**2
	658	enddo
	659	enddo
[938]	660	! print*,'14g OK convect8'
[878]	661	do l=1,nlay
	662	do ig=1,ngrid
	663	ratqsdiff(ig,l)=sqrt(vardiff)/(po(ig,l)*1000.)
	664	! write(11,*)'ratqsdiff=',ratqsdiff(ig,l)
	665	enddo
	666	enddo
	667	!--------------------------------------------------------------------
	668	!
	669	!ecriture des fichiers sortie
	670	! print*,'15 OK convect8'
	671
	672	isplit=isplit+1
	673
	674
	675	#ifdef und
[938]	676	if (prt_level.ge.1) print*,'thermcell_main sorties 1D'
[878]	677	#include "thermcell_out1d.h"
	678	#endif
	679
	680
	681	! #define troisD
[938]	682	if (prt_level.ge.1) print*,'thermcell_main sorties 3D'
[878]	683	#ifdef troisD
	684	#include "thermcell_out3d.h"
	685	#endif
	686
	687	endif
	688
[938]	689	if (prt_level.ge.1) print*,'thermcell_main FIN OK'
[878]	690
	691	return
	692	end
	693
	694	!-----------------------------------------------------------------------------
	695
	696	subroutine test_ltherm(klon,klev,pplev,pplay,long,seuil,ztv,po,ztva,zqla,f_star,zw2,comment)
[938]	697	IMPLICIT NONE
	698	#include "iniprint.h"
[878]	699
[938]	700	integer i, k, klon,klev
[878]	701	real pplev(klon,klev+1),pplay(klon,klev)
	702	real ztv(klon,klev)
	703	real po(klon,klev)
	704	real ztva(klon,klev)
	705	real zqla(klon,klev)
	706	real f_star(klon,klev)
	707	real zw2(klon,klev)
	708	integer long(klon)
	709	real seuil
	710	character*21 comment
	711
[938]	712	if (prt_level.ge.1) THEN
	713	print*,'WARNING !!! TEST ',comment
	714	endif
[879]	715	return
	716
[878]	717	! test sur la hauteur des thermiques ...
	718	do i=1,klon
	719	if (pplay(i,long(i)).lt.seuil*pplev(i,1)) then
	720	print*,'WARNING ',comment,' au point ',i,' K= ',long(i)
	721	print,' K P(MB) THV(K) Qenv(g/kg)THVA QLA(g/kg) F W2'
	722	do k=1,klev
	723	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)
	724	enddo
	725	! stop
	726	endif
	727	enddo
	728
	729
	730	return
	731	end
	732

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

Download in other formats:

Original Format