Context Navigation

thermcell.F @ 5456

Last change on this file since 5456 was 634, checked in by Laurent Fairhead, 20 years ago
Modifications faites à la physique pour la rendre parallele YM Une branche de travail LMDZ4_par_0 a été créée provisoirement afin de tester les modifs pleinement avant leurs inclusions dans le tronc principal LF
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 41.1 KB

Rev	Line
[542]	1	SUBROUTINE thermcell(ngrid,nlay,ptimestep
	2	s ,pplay,pplev,pphi
	3	s ,pu,pv,pt,po
	4	s ,pduadj,pdvadj,pdtadj,pdoadj
	5	s ,fm0,entr0
	6	c s ,pu_therm,pv_therm
	7	s ,r_aspect,l_mix,w2di,tho)
[634]	8	USE dimphy
[542]	9	IMPLICIT NONE
	10
	11	c=======================================================================
	12	c
	13	c Calcul du transport verticale dans la couche limite en presence
	14	c de "thermiques" explicitement representes
	15	c
	16	c Réécriture à partir d'un listing papier à Habas, le 14/02/00
	17	c
	18	c le thermique est supposé homogène et dissipé par mélange avec
	19	c son environnement. la longueur l_mix contrôle l'efficacité du
	20	c mélange
	21	c
	22	c Le calcul du transport des différentes espèces se fait en prenant
	23	c en compte:
	24	c 1. un flux de masse montant
	25	c 2. un flux de masse descendant
	26	c 3. un entrainement
	27	c 4. un detrainement
	28	c
	29	c=======================================================================
	30
	31	c-----------------------------------------------------------------------
	32	c declarations:
	33	c -------------
	34
[634]	35	cym#include "dimensions.h"
	36	cym#include "dimphy.h"
[542]	37	#include "YOMCST.h"
	38
	39	c arguments:
	40	c ----------
	41
	42	INTEGER ngrid,nlay,w2di,tho
	43	real ptimestep,l_mix,r_aspect
	44	REAL pt(ngrid,nlay),pdtadj(ngrid,nlay)
	45	REAL pu(ngrid,nlay),pduadj(ngrid,nlay)
	46	REAL pv(ngrid,nlay),pdvadj(ngrid,nlay)
	47	REAL po(ngrid,nlay),pdoadj(ngrid,nlay)
	48	REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
	49	real pphi(ngrid,nlay)
	50
	51	integer idetr
	52	save idetr
	53	data idetr/3/
	54
	55	c local:
	56	c ------
	57
	58	INTEGER ig,k,l,lmaxa(klon),lmix(klon)
	59	real zsortie1d(klon)
	60	c CR: on remplace lmax(klon,klev+1)
	61	INTEGER lmax(klon),lmin(klon),lentr(klon)
	62	real linter(klon)
	63	real zmix(klon), fracazmix(klon)
	64	c RC
	65	real zmax(klon),zw,zz,zw2(klon,klev+1),ztva(klon,klev),zzz
	66
	67	real zlev(klon,klev+1),zlay(klon,klev)
	68	REAL zh(klon,klev),zdhadj(klon,klev)
	69	REAL ztv(klon,klev)
	70	real zu(klon,klev),zv(klon,klev),zo(klon,klev)
	71	REAL wh(klon,klev+1)
	72	real wu(klon,klev+1),wv(klon,klev+1),wo(klon,klev+1)
	73	real zla(klon,klev+1)
	74	real zwa(klon,klev+1)
	75	real zld(klon,klev+1)
	76	real zwd(klon,klev+1)
	77	real zsortie(klon,klev)
	78	real zva(klon,klev)
	79	real zua(klon,klev)
	80	real zoa(klon,klev)
	81
	82	real zha(klon,klev)
	83	real wa_moy(klon,klev+1)
	84	real fraca(klon,klev+1)
	85	real fracc(klon,klev+1)
	86	real zf,zf2
[634]	87	real,allocatable,save :: thetath2(:,:),wth2(:,:)
	88	cym common/comtherm/thetath2,wth2
[542]	89
	90	real count_time
	91	integer isplit,nsplit,ialt
	92	parameter (nsplit=10)
	93	data isplit/0/
	94	save isplit
	95
	96	logical sorties
	97	real rho(klon,klev),rhobarz(klon,klev+1),masse(klon,klev)
	98	real zpspsk(klon,klev)
	99
	100	c real wmax(klon,klev),wmaxa(klon)
	101	real wmax(klon),wmaxa(klon)
	102	real wa(klon,klev,klev+1)
	103	real wd(klon,klev+1)
	104	real larg_part(klon,klev,klev+1)
	105	real fracd(klon,klev+1)
	106	real xxx(klon,klev+1)
	107	real larg_cons(klon,klev+1)
	108	real larg_detr(klon,klev+1)
	109	real fm0(klon,klev+1),entr0(klon,klev),detr(klon,klev)
	110	real pu_therm(klon,klev),pv_therm(klon,klev)
	111	real fm(klon,klev+1),entr(klon,klev)
	112	real fmc(klon,klev+1)
	113
	114	cCR:nouvelles variables
	115	real f_star(klon,klev+1),entr_star(klon,klev)
	116	real entr_star_tot(klon),entr_star2(klon)
	117	real f(klon), f0(klon)
	118	real zlevinter(klon)
	119	logical first
	120	data first /.false./
	121	save first
	122	cRC
	123
	124	character*2 str2
	125	character*10 str10
	126
	127	LOGICAL vtest(klon),down
	128
	129	EXTERNAL SCOPY
	130
	131	integer ncorrec,ll
	132	save ncorrec
	133	data ncorrec/0/
[634]	134	logical,save :: firstCall=.true.
[542]	135	c
	136	c-----------------------------------------------------------------------
	137	c initialisation:
	138	c ---------------
	139	c
[634]	140	if (firstcall) then
	141	allocate(thetath2(klon,klev),wth2(klon,klev))
	142	thetath2(:,:)=0.
	143	wth2(:,:)=0.
	144	firstcall=.false.
	145	endif
	146
[542]	147	sorties=.true.
	148	IF(ngrid.NE.klon) THEN
	149	PRINT*
	150	PRINT*,'STOP dans convadj'
	151	PRINT*,'ngrid =',ngrid
	152	PRINT*,'klon =',klon
	153	ENDIF
	154	c
	155	c-----------------------------------------------------------------------
	156	c incrementation eventuelle de tendances precedentes:
	157	c ---------------------------------------------------
	158
	159	print*,'0 OK convect8'
	160
	161	DO 1010 l=1,nlay
	162	DO 1015 ig=1,ngrid
	163	zpspsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**RKAPPA
	164	zh(ig,l)=pt(ig,l)/zpspsk(ig,l)
	165	zu(ig,l)=pu(ig,l)
	166	zv(ig,l)=pv(ig,l)
	167	zo(ig,l)=po(ig,l)
	168	ztv(ig,l)=zh(ig,l)(1.+0.61zo(ig,l))
	169	1015 CONTINUE
	170	1010 CONTINUE
	171
	172	print*,'1 OK convect8'
	173	c --------------------
	174	c
	175	c
	176	c + + + + + + + + + + +
	177	c
	178	c
	179	c wa, fraca, wd, fracd -------------------- zlev(2), rhobarz
	180	c wh,wt,wo ...
	181	c
	182	c + + + + + + + + + + + zh,zu,zv,zo,rho
	183	c
	184	c
	185	c -------------------- zlev(1)
	186	c \\\\\\\\\\\\\\\\\\\\
	187	c
	188	c
	189
	190	c-----------------------------------------------------------------------
	191	c Calcul des altitudes des couches
	192	c-----------------------------------------------------------------------
	193
	194	do l=2,nlay
	195	do ig=1,ngrid
	196	zlev(ig,l)=0.5*(pphi(ig,l)+pphi(ig,l-1))/RG
	197	enddo
	198	enddo
	199	do ig=1,ngrid
	200	zlev(ig,1)=0.
	201	zlev(ig,nlay+1)=(2.*pphi(ig,klev)-pphi(ig,klev-1))/RG
	202	enddo
	203	do l=1,nlay
	204	do ig=1,ngrid
	205	zlay(ig,l)=pphi(ig,l)/RG
	206	enddo
	207	enddo
	208
	209	c print*,'2 OK convect8'
	210	c-----------------------------------------------------------------------
	211	c Calcul des densites
	212	c-----------------------------------------------------------------------
	213
	214	do l=1,nlay
	215	do ig=1,ngrid
	216	rho(ig,l)=pplay(ig,l)/(zpspsk(ig,l)RDzh(ig,l))
	217	enddo
	218	enddo
	219
	220	do l=2,nlay
	221	do ig=1,ngrid
	222	rhobarz(ig,l)=0.5*(rho(ig,l)+rho(ig,l-1))
	223	enddo
	224	enddo
	225
	226	do k=1,nlay
	227	do l=1,nlay+1
	228	do ig=1,ngrid
	229	wa(ig,k,l)=0.
	230	enddo
	231	enddo
	232	enddo
	233
	234	c print*,'3 OK convect8'
	235	c------------------------------------------------------------------
	236	c Calcul de w2, quarre de w a partir de la cape
	237	c a partir de w2, on calcule wa, vitesse de l'ascendance
	238	c
	239	c ATTENTION: Dans cette version, pour cause d'economie de memoire,
	240	c w2 est stoke dans wa
	241	c
	242	c ATTENTION: dans convect8, on n'utilise le calcule des wa
	243	c independants par couches que pour calculer l'entrainement
	244	c a la base et la hauteur max de l'ascendance.
	245	c
	246	c Indicages:
	247	c l'ascendance provenant du niveau k traverse l'interface l avec
	248	c une vitesse wa(k,l).
	249	c
	250	c --------------------
	251	c
	252	c + + + + + + + + + +
	253	c
	254	c wa(k,l) ---- -------------------- l
	255	c /\
	256	c /\|\|\ + + + + + + + + + +
	257	c \|\|
	258	c \|\| --------------------
	259	c \|\|
	260	c \|\| + + + + + + + + + +
	261	c \|\|
	262	c \|\| --------------------
	263	c \|\|__
	264	c \|___ + + + + + + + + + + k
	265	c
	266	c --------------------
	267	c
	268	c
	269	c
	270	c------------------------------------------------------------------
	271
	272	cCR: ponderation entrainement des couches instables
	273	cdef des entr_star tels que entr=f*entr_star
	274	do l=1,klev
	275	do ig=1,ngrid
	276	entr_star(ig,l)=0.
	277	enddo
	278	enddo
	279	c determination de la longueur de la couche d entrainement
	280	do ig=1,ngrid
	281	lentr(ig)=1
	282	enddo
	283
	284	con ne considere que les premieres couches instables
	285	do k=nlay-2,1,-1
	286	do ig=1,ngrid
	287	if (ztv(ig,k).gt.ztv(ig,k+1).and.
	288	s ztv(ig,k+1).le.ztv(ig,k+2)) then
	289	lentr(ig)=k
	290	endif
	291	enddo
	292	enddo
	293
	294	c determination du lmin: couche d ou provient le thermique
	295	do ig=1,ngrid
	296	lmin(ig)=1
	297	enddo
	298	do ig=1,ngrid
	299	do l=nlay,2,-1
	300	if (ztv(ig,l-1).gt.ztv(ig,l)) then
	301	lmin(ig)=l-1
	302	endif
	303	enddo
	304	enddo
	305	c
	306	c definition de l'entrainement des couches
	307	do l=1,klev-1
	308	do ig=1,ngrid
	309	if (ztv(ig,l).gt.ztv(ig,l+1).and.
	310	s l.ge.lmin(ig).and.l.le.lentr(ig)) then
	311	entr_star(ig,l)=(ztv(ig,l)-ztv(ig,l+1))*
	312	s (zlev(ig,l+1)-zlev(ig,l))
	313	endif
	314	enddo
	315	enddo
	316	c pas de thermique si couches 1->5 stables
	317	do ig=1,ngrid
	318	if (lmin(ig).gt.5) then
	319	do l=1,klev
	320	entr_star(ig,l)=0.
	321	enddo
	322	endif
	323	enddo
	324	c calcul de l entrainement total
	325	do ig=1,ngrid
	326	entr_star_tot(ig)=0.
	327	enddo
	328	do ig=1,ngrid
	329	do k=1,klev
	330	entr_star_tot(ig)=entr_star_tot(ig)+entr_star(ig,k)
	331	enddo
	332	enddo
	333	c
	334	print*,'fin calcul entr_star'
	335	do k=1,klev
	336	do ig=1,ngrid
	337	ztva(ig,k)=ztv(ig,k)
	338	enddo
	339	enddo
	340	cRC
	341	c print*,'7 OK convect8'
	342	do k=1,klev+1
	343	do ig=1,ngrid
	344	zw2(ig,k)=0.
	345	fmc(ig,k)=0.
	346	cCR
	347	f_star(ig,k)=0.
	348	cRC
	349	larg_cons(ig,k)=0.
	350	larg_detr(ig,k)=0.
	351	wa_moy(ig,k)=0.
	352	enddo
	353	enddo
	354
	355	c print*,'8 OK convect8'
	356	do ig=1,ngrid
	357	linter(ig)=1.
	358	lmaxa(ig)=1
	359	lmix(ig)=1
	360	wmaxa(ig)=0.
	361	enddo
	362
	363	cCR:
	364	do l=1,nlay-2
	365	do ig=1,ngrid
	366	if (ztv(ig,l).gt.ztv(ig,l+1)
	367	s .and.entr_star(ig,l).gt.1.e-10
	368	s .and.zw2(ig,l).lt.1e-10) then
	369	f_star(ig,l+1)=entr_star(ig,l)
	370	ctest:calcul de dteta
	371	zw2(ig,l+1)=2.RG(ztv(ig,l)-ztv(ig,l+1))/ztv(ig,l+1)
	372	s *(zlev(ig,l+1)-zlev(ig,l))
	373	s 0.4pphi(ig,l)/(pphi(ig,l+1)-pphi(ig,l))
	374	larg_detr(ig,l)=0.
	375	else if ((zw2(ig,l).ge.1e-10).and.
	376	s (f_star(ig,l)+entr_star(ig,l).gt.1.e-10)) then
	377	f_star(ig,l+1)=f_star(ig,l)+entr_star(ig,l)
	378	ztva(ig,l)=(f_star(ig,l)*ztva(ig,l-1)+entr_star(ig,l)
	379	s *ztv(ig,l))/f_star(ig,l+1)
	380	zw2(ig,l+1)=zw2(ig,l)(f_star(ig,l)/f_star(ig,l+1))*2+
	381	s 2.RG(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
	382	s *(zlev(ig,l+1)-zlev(ig,l))
	383	endif
	384	c determination de zmax continu par interpolation lineaire
	385	if (zw2(ig,l+1).lt.0.) then
	386	ctest
	387	if (abs(zw2(ig,l+1)-zw2(ig,l)).lt.1e-10) then
	388	print*,'pb linter'
	389	endif
	390	linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l))
	391	s -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
	392	zw2(ig,l+1)=0.
	393	lmaxa(ig)=l
	394	else
	395	if (zw2(ig,l+1).lt.0.) then
	396	print*,'pb1 zw2<0'
	397	endif
	398	wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
	399	endif
	400	if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
	401	c lmix est le niveau de la couche ou w (wa_moy) est maximum
	402	lmix(ig)=l+1
	403	wmaxa(ig)=wa_moy(ig,l+1)
	404	endif
	405	enddo
	406	enddo
	407	print*,'fin calcul zw2'
	408	c
	409	c Calcul de la couche correspondant a la hauteur du thermique
	410	do ig=1,ngrid
	411	lmax(ig)=lentr(ig)
	412	enddo
	413	do ig=1,ngrid
	414	do l=nlay,lentr(ig)+1,-1
	415	if (zw2(ig,l).le.1.e-10) then
	416	lmax(ig)=l-1
	417	endif
	418	enddo
	419	enddo
	420	c pas de thermique si couches 1->5 stables
	421	do ig=1,ngrid
	422	if (lmin(ig).gt.5) then
	423	lmax(ig)=1
	424	lmin(ig)=1
	425	endif
	426	enddo
	427	c
	428	c Determination de zw2 max
	429	do ig=1,ngrid
	430	wmax(ig)=0.
	431	enddo
	432
	433	do l=1,nlay
	434	do ig=1,ngrid
	435	if (l.le.lmax(ig)) then
	436	if (zw2(ig,l).lt.0.)then
	437	print*,'pb2 zw2<0'
	438	endif
	439	zw2(ig,l)=sqrt(zw2(ig,l))
	440	wmax(ig)=max(wmax(ig),zw2(ig,l))
	441	else
	442	zw2(ig,l)=0.
	443	endif
	444	enddo
	445	enddo
	446
	447	c Longueur caracteristique correspondant a la hauteur des thermiques.
	448	do ig=1,ngrid
	449	zmax(ig)=0.
	450	zlevinter(ig)=zlev(ig,1)
	451	enddo
	452	do ig=1,ngrid
	453	c calcul de zlevinter
	454	zlevinter(ig)=(zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))*
	455	s linter(ig)+zlev(ig,lmax(ig))-lmax(ig)*(zlev(ig,lmax(ig)+1)
	456	s -zlev(ig,lmax(ig)))
	457	zmax(ig)=max(zmax(ig),zlevinter(ig)-zlev(ig,lmin(ig)))
	458	enddo
	459
	460	print*,'avant fermeture'
	461	c Fermeture,determination de f
	462	do ig=1,ngrid
	463	entr_star2(ig)=0.
	464	enddo
	465	do ig=1,ngrid
	466	if (entr_star_tot(ig).LT.1.e-10) then
	467	f(ig)=0.
	468	else
	469	do k=lmin(ig),lentr(ig)
	470	entr_star2(ig)=entr_star2(ig)+entr_star(ig,k)**2
	471	s /(rho(ig,k)*(zlev(ig,k+1)-zlev(ig,k)))
	472	enddo
	473	c Nouvelle fermeture
	474	f(ig)=wmax(ig)/(max(500.,zmax(ig))*r_aspect
	475	s entr_star2(ig))entr_star_tot(ig)
	476	ctest
	477	c if (first) then
	478	c f(ig)=f(ig)+(f0(ig)-f(ig))*exp(-ptimestep/zmax(ig)
	479	c s *wmax(ig))
	480	c endif
	481	endif
	482	c f0(ig)=f(ig)
	483	c first=.true.
	484	enddo
	485	print*,'apres fermeture'
	486
	487	c Calcul de l'entrainement
	488	do k=1,klev
	489	do ig=1,ngrid
	490	entr(ig,k)=f(ig)*entr_star(ig,k)
	491	enddo
	492	enddo
	493	c Calcul des flux
	494	do ig=1,ngrid
	495	do l=1,lmax(ig)-1
	496	fmc(ig,l+1)=fmc(ig,l)+entr(ig,l)
	497	enddo
	498	enddo
	499
	500	cRC
	501
	502
	503	c print*,'9 OK convect8'
	504	c print*,'WA1 ',wa_moy
	505
	506	c determination de l'indice du debut de la mixed layer ou w decroit
	507
	508	c calcul de la largeur de chaque ascendance dans le cas conservatif.
	509	c dans ce cas simple, on suppose que la largeur de l'ascendance provenant
	510	c d'une couche est égale à la hauteur de la couche alimentante.
	511	c La vitesse maximale dans l'ascendance est aussi prise comme estimation
	512	c de la vitesse d'entrainement horizontal dans la couche alimentante.
	513
	514	do l=2,nlay
	515	do ig=1,ngrid
	516	if (l.le.lmaxa(ig)) then
	517	zw=max(wa_moy(ig,l),1.e-10)
	518	larg_cons(ig,l)=zmax(ig)*r_aspect
	519	s fmc(ig,l)/(rhobarz(ig,l)zw)
	520	endif
	521	enddo
	522	enddo
	523
	524	do l=2,nlay
	525	do ig=1,ngrid
	526	if (l.le.lmaxa(ig)) then
	527	c if (idetr.eq.0) then
	528	c cette option est finalement en dur.
	529	if ((l_mix*zlev(ig,l)).lt.0.)then
	530	print,'pb l_mixzlev<0'
	531	endif
	532	larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
	533	c else if (idetr.eq.1) then
	534	c larg_detr(ig,l)=larg_cons(ig,l)
	535	c s sqrt(l_mixzlev(ig,l))/larg_cons(ig,lmix(ig))
	536	c else if (idetr.eq.2) then
	537	c larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
	538	c s *sqrt(wa_moy(ig,l))
	539	c else if (idetr.eq.4) then
	540	c larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
	541	c s *wa_moy(ig,l)
	542	c endif
	543	endif
	544	enddo
	545	enddo
	546
	547	c print*,'10 OK convect8'
	548	c print*,'WA2 ',wa_moy
	549	c calcul de la fraction de la maille concernée par l'ascendance en tenant
	550	c compte de l'epluchage du thermique.
	551	c
	552	cCR def de zmix continu (profil parabolique des vitesses)
	553	do ig=1,ngrid
	554	if (lmix(ig).gt.1.) then
	555	c test
	556	if (((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))
	557	s *((zlev(ig,lmix(ig)))-(zlev(ig,lmix(ig)+1)))
	558	s -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1))
	559	s *((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig))))).gt.1e-10)
	560	s then
	561	c
	562	zmix(ig)=((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))
	563	s ((zlev(ig,lmix(ig)))2-(zlev(ig,lmix(ig)+1))*2)
	564	s -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1))
	565	s ((zlev(ig,lmix(ig)-1))2-(zlev(ig,lmix(ig)))*2))
	566	s /(2.*((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))
	567	s *((zlev(ig,lmix(ig)))-(zlev(ig,lmix(ig)+1)))
	568	s -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1))
	569	s *((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig))))))
	570	else
	571	zmix(ig)=zlev(ig,lmix(ig))
	572	print*,'pb zmix'
	573	endif
	574	else
	575	zmix(ig)=0.
	576	endif
	577	ctest
	578	if ((zmax(ig)-zmix(ig)).lt.0.) then
	579	zmix(ig)=0.99*zmax(ig)
	580	c print*,'pb zmix>zmax'
	581	endif
	582	enddo
	583	c
	584	c calcul du nouveau lmix correspondant
	585	do ig=1,ngrid
	586	do l=1,klev
	587	if (zmix(ig).ge.zlev(ig,l).and.
	588	s zmix(ig).lt.zlev(ig,l+1)) then
	589	lmix(ig)=l
	590	endif
	591	enddo
	592	enddo
	593	c
	594	do l=2,nlay
	595	do ig=1,ngrid
	596	if(larg_cons(ig,l).gt.1.) then
	597	c print*,ig,l,lmix(ig),lmaxa(ig),larg_cons(ig,l),' KKK'
	598	fraca(ig,l)=(larg_cons(ig,l)-larg_detr(ig,l))
	599	s /(r_aspect*zmax(ig))
	600	c test
	601	fraca(ig,l)=max(fraca(ig,l),0.)
	602	fraca(ig,l)=min(fraca(ig,l),0.5)
	603	fracd(ig,l)=1.-fraca(ig,l)
	604	fracc(ig,l)=larg_cons(ig,l)/(r_aspect*zmax(ig))
	605	else
	606	c wa_moy(ig,l)=0.
	607	fraca(ig,l)=0.
	608	fracc(ig,l)=0.
	609	fracd(ig,l)=1.
	610	endif
	611	enddo
	612	enddo
	613	cCR: calcul de fracazmix
	614	do ig=1,ngrid
	615	fracazmix(ig)=(fraca(ig,lmix(ig)+1)-fraca(ig,lmix(ig)))/
	616	s (zlev(ig,lmix(ig)+1)-zlev(ig,lmix(ig)))*zmix(ig)
	617	s +fraca(ig,lmix(ig))-zlev(ig,lmix(ig))*(fraca(ig,lmix(ig)+1)
	618	s -fraca(ig,lmix(ig)))/(zlev(ig,lmix(ig)+1)-zlev(ig,lmix(ig)))
	619	enddo
	620	c
	621	do l=2,nlay
	622	do ig=1,ngrid
	623	if(larg_cons(ig,l).gt.1.) then
	624	if (l.gt.lmix(ig)) then
	625	ctest
	626	if (zmax(ig)-zmix(ig).lt.1.e-10) then
	627	c print*,'pb xxx'
	628	xxx(ig,l)=(lmaxa(ig)+1.-l)/(lmaxa(ig)+1.-lmix(ig))
	629	else
	630	xxx(ig,l)=(zmax(ig)-zlev(ig,l))/(zmax(ig)-zmix(ig))
	631	endif
	632	if (idetr.eq.0) then
	633	fraca(ig,l)=fracazmix(ig)
	634	else if (idetr.eq.1) then
	635	fraca(ig,l)=fracazmix(ig)*xxx(ig,l)
	636	else if (idetr.eq.2) then
	637	fraca(ig,l)=fracazmix(ig)(1.-(1.-xxx(ig,l))*2)
	638	else
	639	fraca(ig,l)=fracazmix(ig)xxx(ig,l)*2
	640	endif
	641	c print*,ig,l,lmix(ig),lmaxa(ig),xxx(ig,l),'LLLLLLL'
	642	fraca(ig,l)=max(fraca(ig,l),0.)
	643	fraca(ig,l)=min(fraca(ig,l),0.5)
	644	fracd(ig,l)=1.-fraca(ig,l)
	645	fracc(ig,l)=larg_cons(ig,l)/(r_aspect*zmax(ig))
	646	endif
	647	endif
	648	enddo
	649	enddo
	650
	651	print*,'fin calcul fraca'
	652	c print*,'11 OK convect8'
	653	c print*,'Ea3 ',wa_moy
	654	c------------------------------------------------------------------
	655	c Calcul de fracd, wd
	656	c somme wa - wd = 0
	657	c------------------------------------------------------------------
	658
	659
	660	do ig=1,ngrid
	661	fm(ig,1)=0.
	662	fm(ig,nlay+1)=0.
	663	enddo
	664
	665	do l=2,nlay
	666	do ig=1,ngrid
	667	fm(ig,l)=fraca(ig,l)wa_moy(ig,l)rhobarz(ig,l)
	668	cCR:test
	669	if (entr(ig,l-1).lt.1e-10.and.fm(ig,l).gt.fm(ig,l-1)
	670	s .and.l.gt.lmix(ig)) then
	671	fm(ig,l)=fm(ig,l-1)
	672	c write(1,*)'ajustement fm, l',l
	673	endif
	674	c write(1,*)'ig,l,fm(ig,l)',ig,l,fm(ig,l)
	675	cRC
	676	enddo
	677	do ig=1,ngrid
	678	if(fracd(ig,l).lt.0.1) then
	679	stop'fracd trop petit'
	680	else
	681	c vitesse descendante "diagnostique"
	682	wd(ig,l)=fm(ig,l)/(fracd(ig,l)*rhobarz(ig,l))
	683	endif
	684	enddo
	685	enddo
	686
	687	do l=1,nlay
	688	do ig=1,ngrid
	689	c masse(ig,l)=rho(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
	690	masse(ig,l)=(pplev(ig,l)-pplev(ig,l+1))/RG
	691	enddo
	692	enddo
	693
	694	print*,'12 OK convect8'
	695	c print*,'WA4 ',wa_moy
	696	cc------------------------------------------------------------------
	697	c calcul du transport vertical
	698	c------------------------------------------------------------------
	699
	700	go to 4444
	701	c print*,'XXXXXXXXXXXXXXX ptimestep= ',ptimestep
	702	do l=2,nlay-1
	703	do ig=1,ngrid
	704	if(fm(ig,l+1)*ptimestep.gt.masse(ig,l)
	705	s .and.fm(ig,l+1)*ptimestep.gt.masse(ig,l+1)) then
	706	c print*,'WARN!!! FM>M ig=',ig,' l=',l,' FM='
	707	c s ,fm(ig,l+1)*ptimestep
	708	c s ,' M=',masse(ig,l),masse(ig,l+1)
	709	endif
	710	enddo
	711	enddo
	712
	713	do l=1,nlay
	714	do ig=1,ngrid
	715	if(entr(ig,l)*ptimestep.gt.masse(ig,l)) then
	716	c print*,'WARN!!! E>M ig=',ig,' l=',l,' E=='
	717	c s ,entr(ig,l)*ptimestep
	718	c s ,' M=',masse(ig,l)
	719	endif
	720	enddo
	721	enddo
	722
	723	do l=1,nlay
	724	do ig=1,ngrid
	725	if(.not.fm(ig,l).ge.0..or..not.fm(ig,l).le.10.) then
	726	c print*,'WARN!!! fm exagere ig=',ig,' l=',l
	727	c s ,' FM=',fm(ig,l)
	728	endif
	729	if(.not.masse(ig,l).ge.1.e-10
	730	s .or..not.masse(ig,l).le.1.e4) then
	731	c print*,'WARN!!! masse exagere ig=',ig,' l=',l
	732	c s ,' M=',masse(ig,l)
	733	c print*,'rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)',
	734	c s rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)
	735	c print*,'zlev(ig,l+1),zlev(ig,l)'
	736	c s ,zlev(ig,l+1),zlev(ig,l)
	737	c print*,'pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)'
	738	c s ,pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)
	739	endif
	740	if(.not.entr(ig,l).ge.0..or..not.entr(ig,l).le.10.) then
	741	c print*,'WARN!!! entr exagere ig=',ig,' l=',l
	742	c s ,' E=',entr(ig,l)
	743	endif
	744	enddo
	745	enddo
	746
	747	4444 continue
	748
	749	cCR:redefinition du entr
	750	do l=1,nlay
	751	do ig=1,ngrid
	752	detr(ig,l)=fm(ig,l)+entr(ig,l)-fm(ig,l+1)
	753	if (detr(ig,l).lt.0.) then
	754	entr(ig,l)=entr(ig,l)-detr(ig,l)
	755	detr(ig,l)=0.
	756	c print*,'WARNING !!! detrainement negatif ',ig,l
	757	endif
	758	enddo
	759	enddo
	760	cRC
	761	if (w2di.eq.1) then
	762	fm0=fm0+ptimestep*(fm-fm0)/float(tho)
	763	entr0=entr0+ptimestep*(entr-entr0)/float(tho)
	764	else
	765	fm0=fm
	766	entr0=entr
	767	endif
	768
	769	if (1.eq.1) then
	770	call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
	771	. ,zh,zdhadj,zha)
	772	call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
	773	. ,zo,pdoadj,zoa)
	774	else
	775	call dqthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse,fraca
	776	. ,zh,zdhadj,zha)
	777	call dqthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse,fraca
	778	. ,zo,pdoadj,zoa)
	779	endif
	780
	781	if (1.eq.0) then
	782	call dvthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse
	783	. ,fraca,zmax
	784	. ,zu,zv,pduadj,pdvadj,zua,zva)
	785	else
	786	call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
	787	. ,zu,pduadj,zua)
	788	call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
	789	. ,zv,pdvadj,zva)
	790	endif
	791
	792	do l=1,nlay
	793	do ig=1,ngrid
	794	zf=0.5*(fracc(ig,l)+fracc(ig,l+1))
	795	zf2=zf/(1.-zf)
	796	thetath2(ig,l)=zf2(zha(ig,l)-zh(ig,l))*2
	797	wth2(ig,l)=zf2(0.5(wa_moy(ig,l)+wa_moy(ig,l+1)))**2
	798	enddo
	799	enddo
	800
	801
	802
	803	c print*,'13 OK convect8'
	804	c print*,'WA5 ',wa_moy
	805	do l=1,nlay
	806	do ig=1,ngrid
	807	pdtadj(ig,l)=zdhadj(ig,l)*zpspsk(ig,l)
	808	enddo
	809	enddo
	810
	811
	812	c do l=1,nlay
	813	c do ig=1,ngrid
	814	c if(abs(pdtadj(ig,l))*86400..gt.500.) then
	815	c print*,'WARN!!! ig=',ig,' l=',l
	816	c s ,' pdtadj=',pdtadj(ig,l)
	817	c endif
	818	c if(abs(pdoadj(ig,l))*86400..gt.1.) then
	819	c print*,'WARN!!! ig=',ig,' l=',l
	820	c s ,' pdoadj=',pdoadj(ig,l)
	821	c endif
	822	c enddo
	823	c enddo
	824
	825	print*,'14 OK convect8'
	826	c------------------------------------------------------------------
	827	c Calculs pour les sorties
	828	c------------------------------------------------------------------
	829
	830	if(sorties) then
	831	do l=1,nlay
	832	do ig=1,ngrid
	833	zla(ig,l)=(1.-fracd(ig,l))*zmax(ig)
	834	zld(ig,l)=fracd(ig,l)*zmax(ig)
	835	if(1.-fracd(ig,l).gt.1.e-10)
	836	s zwa(ig,l)=wd(ig,l)*fracd(ig,l)/(1.-fracd(ig,l))
	837	enddo
	838	enddo
	839
	840	cdeja fait
	841	c do l=1,nlay
	842	c do ig=1,ngrid
	843	c detr(ig,l)=fm(ig,l)+entr(ig,l)-fm(ig,l+1)
	844	c if (detr(ig,l).lt.0.) then
	845	c entr(ig,l)=entr(ig,l)-detr(ig,l)
	846	c detr(ig,l)=0.
	847	c print*,'WARNING !!! detrainement negatif ',ig,l
	848	c endif
	849	c enddo
	850	c enddo
	851
	852	c print*,'15 OK convect8'
	853
	854	isplit=isplit+1
	855
	856
	857	c #define und
	858	goto 123
	859	#ifdef und
	860	CALL writeg1d(1,nlay,wd,'wd ','wd ')
	861	CALL writeg1d(1,nlay,zwa,'wa ','wa ')
	862	CALL writeg1d(1,nlay,fracd,'fracd ','fracd ')
	863	CALL writeg1d(1,nlay,fraca,'fraca ','fraca ')
	864	CALL writeg1d(1,nlay,wa_moy,'wam ','wam ')
	865	CALL writeg1d(1,nlay,zla,'la ','la ')
	866	CALL writeg1d(1,nlay,zld,'ld ','ld ')
	867	CALL writeg1d(1,nlay,pt,'pt ','pt ')
	868	CALL writeg1d(1,nlay,zh,'zh ','zh ')
	869	CALL writeg1d(1,nlay,zha,'zha ','zha ')
	870	CALL writeg1d(1,nlay,zu,'zu ','zu ')
	871	CALL writeg1d(1,nlay,zv,'zv ','zv ')
	872	CALL writeg1d(1,nlay,zo,'zo ','zo ')
	873	CALL writeg1d(1,nlay,wh,'wh ','wh ')
	874	CALL writeg1d(1,nlay,wu,'wu ','wu ')
	875	CALL writeg1d(1,nlay,wv,'wv ','wv ')
	876	CALL writeg1d(1,nlay,wo,'w15uo ','wXo ')
	877	CALL writeg1d(1,nlay,zdhadj,'zdhadj ','zdhadj ')
	878	CALL writeg1d(1,nlay,pduadj,'pduadj ','pduadj ')
	879	CALL writeg1d(1,nlay,pdvadj,'pdvadj ','pdvadj ')
	880	CALL writeg1d(1,nlay,pdoadj,'pdoadj ','pdoadj ')
	881	CALL writeg1d(1,nlay,entr ,'entr ','entr ')
	882	CALL writeg1d(1,nlay,detr ,'detr ','detr ')
	883	CALL writeg1d(1,nlay,fm ,'fm ','fm ')
	884
	885	CALL writeg1d(1,nlay,pdtadj,'pdtadj ','pdtadj ')
	886	CALL writeg1d(1,nlay,pplay,'pplay ','pplay ')
	887	CALL writeg1d(1,nlay,pplev,'pplev ','pplev ')
	888
	889	c recalcul des flux en diagnostique...
	890	c print*,'PAS DE TEMPS ',ptimestep
	891	call dt2F(pplev,pplay,pt,pdtadj,wh)
	892	CALL writeg1d(1,nlay,wh,'wh2 ','wh2 ')
	893	#endif
	894	123 continue
	895	c#define troisD
	896	#ifdef troisD
	897	c if (sorties) then
	898	print*,'Debut des wrgradsfi'
	899
	900	c print*,'16 OK convect8'
	901	call wrgradsfi(1,nlay,wd,'wd ','wd ')
	902	call wrgradsfi(1,nlay,zwa,'wa ','wa ')
	903	call wrgradsfi(1,nlay,fracd,'fracd ','fracd ')
	904	call wrgradsfi(1,nlay,fraca,'fraca ','fraca ')
	905	call wrgradsfi(1,nlay,xxx,'xxx ','xxx ')
	906	call wrgradsfi(1,nlay,wa_moy,'wam ','wam ')
	907	c print*,'WA6 ',wa_moy
	908	call wrgradsfi(1,nlay,zla,'la ','la ')
	909	call wrgradsfi(1,nlay,zld,'ld ','ld ')
	910	call wrgradsfi(1,nlay,pt,'pt ','pt ')
	911	call wrgradsfi(1,nlay,zh,'zh ','zh ')
	912	call wrgradsfi(1,nlay,zha,'zha ','zha ')
	913	call wrgradsfi(1,nlay,zua,'zua ','zua ')
	914	call wrgradsfi(1,nlay,zva,'zva ','zva ')
	915	call wrgradsfi(1,nlay,zu,'zu ','zu ')
	916	call wrgradsfi(1,nlay,zv,'zv ','zv ')
	917	call wrgradsfi(1,nlay,zo,'zo ','zo ')
	918	call wrgradsfi(1,nlay,wh,'wh ','wh ')
	919	call wrgradsfi(1,nlay,wu,'wu ','wu ')
	920	call wrgradsfi(1,nlay,wv,'wv ','wv ')
	921	call wrgradsfi(1,nlay,wo,'wo ','wo ')
	922	call wrgradsfi(1,1,zmax,'zmax ','zmax ')
	923	call wrgradsfi(1,nlay,zdhadj,'zdhadj ','zdhadj ')
	924	call wrgradsfi(1,nlay,pduadj,'pduadj ','pduadj ')
	925	call wrgradsfi(1,nlay,pdvadj,'pdvadj ','pdvadj ')
	926	call wrgradsfi(1,nlay,pdoadj,'pdoadj ','pdoadj ')
	927	call wrgradsfi(1,nlay,entr,'entr ','entr ')
	928	call wrgradsfi(1,nlay,detr,'detr ','detr ')
	929	call wrgradsfi(1,nlay,fm,'fm ','fm ')
	930	call wrgradsfi(1,nlay,fmc,'fmc ','fmc ')
	931	call wrgradsfi(1,nlay,zw2,'zw2 ','zw2 ')
	932	call wrgradsfi(1,nlay,ztva,'ztva ','ztva ')
	933	call wrgradsfi(1,nlay,ztv,'ztv ','ztv ')
	934
	935	call wrgradsfi(1,nlay,zo,'zo ','zo ')
	936	call wrgradsfi(1,nlay,larg_cons,'Lc ','Lc ')
	937	call wrgradsfi(1,nlay,larg_detr,'Ldetr ','Ldetr ')
	938
	939	cCR:nouveaux diagnostiques
	940	call wrgradsfi(1,nlay,entr_star ,'entr_star ','entr_star ')
	941	call wrgradsfi(1,nlay,f_star ,'f_star ','f_star ')
	942	call wrgradsfi(1,1,zmax,'zmax ','zmax ')
	943	call wrgradsfi(1,1,zmix,'zmix ','zmix ')
	944	zsortie1d(:)=lmax(:)
	945	call wrgradsfi(1,1,zsortie1d,'lmax ','lmax ')
	946	call wrgradsfi(1,1,wmax,'wmax ','wmax ')
	947	zsortie1d(:)=lmix(:)
	948	call wrgradsfi(1,1,zsortie1d,'lmix ','lmix ')
	949	zsortie1d(:)=lentr(:)
	950	call wrgradsfi(1,1,zsortie1d,'lentr ','lentr ')
	951
	952	c print*,'17 OK convect8'
	953
	954	do k=1,klev/10
	955	write(str2,'(i2.2)') k
	956	str10='wa'//str2
	957	do l=1,nlay
	958	do ig=1,ngrid
	959	zsortie(ig,l)=wa(ig,k,l)
	960	enddo
	961	enddo
	962	CALL wrgradsfi(1,nlay,zsortie,str10,str10)
	963	do l=1,nlay
	964	do ig=1,ngrid
	965	zsortie(ig,l)=larg_part(ig,k,l)
	966	enddo
	967	enddo
	968	str10='la'//str2
	969	CALL wrgradsfi(1,nlay,zsortie,str10,str10)
	970	enddo
	971
	972
	973	c print*,'18 OK convect8'
	974	c endif
	975	print*,'Fin des wrgradsfi'
	976	#endif
	977
	978	endif
	979
	980	c if(wa_moy(1,4).gt.1.e-10) stop
	981
	982	print*,'19 OK convect8'
	983	return
	984	end
	985
	986	subroutine dqthermcell(ngrid,nlay,ptimestep,fm,entr,masse
	987	. ,q,dq,qa)
[634]	988	use dimphy
[542]	989	implicit none
	990
	991	c=======================================================================
	992	c
	993	c Calcul du transport verticale dans la couche limite en presence
	994	c de "thermiques" explicitement representes
	995	c calcul du dq/dt une fois qu'on connait les ascendances
	996	c
	997	c=======================================================================
	998
[634]	999	cym#include "dimensions.h"
	1000	cym#include "dimphy.h"
[542]	1001
	1002	integer ngrid,nlay
	1003
	1004	real ptimestep
	1005	real masse(ngrid,nlay),fm(ngrid,nlay+1)
	1006	real entr(ngrid,nlay)
	1007	real q(ngrid,nlay)
	1008	real dq(ngrid,nlay)
	1009
	1010	real qa(klon,klev),detr(klon,klev),wqd(klon,klev+1)
	1011
	1012	integer ig,k
	1013
	1014	c calcul du detrainement
	1015
	1016	do k=1,nlay
	1017	do ig=1,ngrid
	1018	detr(ig,k)=fm(ig,k)-fm(ig,k+1)+entr(ig,k)
	1019	enddo
	1020	enddo
	1021
	1022	c calcul de la valeur dans les ascendances
	1023	do ig=1,ngrid
	1024	qa(ig,1)=q(ig,1)
	1025	enddo
	1026
	1027	do k=2,nlay
	1028	do ig=1,ngrid
	1029	if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt.
	1030	s 1.e-5*masse(ig,k)) then
	1031	qa(ig,k)=(fm(ig,k)qa(ig,k-1)+entr(ig,k)q(ig,k))
	1032	s /(fm(ig,k+1)+detr(ig,k))
	1033	else
	1034	qa(ig,k)=q(ig,k)
	1035	endif
	1036	enddo
	1037	enddo
	1038
	1039	do k=2,nlay
	1040	do ig=1,ngrid
	1041	c wqd(ig,k)=fm(ig,k)0.5(q(ig,k-1)+q(ig,k))
	1042	wqd(ig,k)=fm(ig,k)*q(ig,k)
	1043	enddo
	1044	enddo
	1045	do ig=1,ngrid
	1046	wqd(ig,1)=0.
	1047	wqd(ig,nlay+1)=0.
	1048	enddo
	1049
	1050	do k=1,nlay
	1051	do ig=1,ngrid
	1052	dq(ig,k)=(detr(ig,k)qa(ig,k)-entr(ig,k)q(ig,k)
	1053	s -wqd(ig,k)+wqd(ig,k+1))
	1054	s /masse(ig,k)
	1055	enddo
	1056	enddo
	1057
	1058	return
	1059	end
	1060	subroutine dvthermcell(ngrid,nlay,ptimestep,fm,entr,masse
	1061	. ,fraca,larga
	1062	. ,u,v,du,dv,ua,va)
[634]	1063	use dimphy
[542]	1064	implicit none
	1065
	1066	c=======================================================================
	1067	c
	1068	c Calcul du transport verticale dans la couche limite en presence
	1069	c de "thermiques" explicitement representes
	1070	c calcul du dq/dt une fois qu'on connait les ascendances
	1071	c
	1072	c=======================================================================
	1073
[634]	1074	cym#include "dimensions.h"
	1075	cym#include "dimphy.h"
[542]	1076
	1077	integer ngrid,nlay
	1078
	1079	real ptimestep
	1080	real masse(ngrid,nlay),fm(ngrid,nlay+1)
	1081	real fraca(ngrid,nlay+1)
	1082	real larga(ngrid)
	1083	real entr(ngrid,nlay)
	1084	real u(ngrid,nlay)
	1085	real ua(ngrid,nlay)
	1086	real du(ngrid,nlay)
	1087	real v(ngrid,nlay)
	1088	real va(ngrid,nlay)
	1089	real dv(ngrid,nlay)
	1090
	1091	real qa(klon,klev),detr(klon,klev)
	1092	real wvd(klon,klev+1),wud(klon,klev+1)
	1093	real gamma0,gamma(klon,klev+1)
	1094	real dua,dva
	1095	integer iter
	1096
	1097	integer ig,k
	1098
	1099	c calcul du detrainement
	1100
	1101	do k=1,nlay
	1102	do ig=1,ngrid
	1103	detr(ig,k)=fm(ig,k)-fm(ig,k+1)+entr(ig,k)
	1104	enddo
	1105	enddo
	1106
	1107	c calcul de la valeur dans les ascendances
	1108	do ig=1,ngrid
	1109	ua(ig,1)=u(ig,1)
	1110	va(ig,1)=v(ig,1)
	1111	enddo
	1112
	1113	do k=2,nlay
	1114	do ig=1,ngrid
	1115	if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt.
	1116	s 1.e-5*masse(ig,k)) then
	1117	c On itère sur la valeur du coeff de freinage.
	1118	c gamma0=rho(ig,k)*(zlev(ig,k+1)-zlev(ig,k))
	1119	gamma0=masse(ig,k)
	1120	s sqrt( 0.5(fraca(ig,k+1)+fraca(ig,k)) )
	1121	s *0.5/larga(ig)
	1122	c gamma0=0.
	1123	c la première fois on multiplie le coefficient de freinage
	1124	c par le module du vent dans la couche en dessous.
	1125	dua=ua(ig,k-1)-u(ig,k-1)
	1126	dva=va(ig,k-1)-v(ig,k-1)
	1127	do iter=1,5
	1128	gamma(ig,k)=gamma0sqrt(dua2+dva*2)
	1129	ua(ig,k)=(fm(ig,k)*ua(ig,k-1)
	1130	s +(entr(ig,k)+gamma(ig,k))*u(ig,k))
	1131	s /(fm(ig,k+1)+detr(ig,k)+gamma(ig,k))
	1132	va(ig,k)=(fm(ig,k)*va(ig,k-1)
	1133	s +(entr(ig,k)+gamma(ig,k))*v(ig,k))
	1134	s /(fm(ig,k+1)+detr(ig,k)+gamma(ig,k))
	1135	c print*,k,ua(ig,k),va(ig,k),u(ig,k),v(ig,k),dua,dva
	1136	dua=ua(ig,k)-u(ig,k)
	1137	dva=va(ig,k)-v(ig,k)
	1138	enddo
	1139	else
	1140	ua(ig,k)=u(ig,k)
	1141	va(ig,k)=v(ig,k)
	1142	gamma(ig,k)=0.
	1143	endif
	1144	enddo
	1145	enddo
	1146
	1147	do k=2,nlay
	1148	do ig=1,ngrid
	1149	wud(ig,k)=fm(ig,k)*u(ig,k)
	1150	wvd(ig,k)=fm(ig,k)*v(ig,k)
	1151	enddo
	1152	enddo
	1153	do ig=1,ngrid
	1154	wud(ig,1)=0.
	1155	wud(ig,nlay+1)=0.
	1156	wvd(ig,1)=0.
	1157	wvd(ig,nlay+1)=0.
	1158	enddo
	1159
	1160	do k=1,nlay
	1161	do ig=1,ngrid
	1162	du(ig,k)=((detr(ig,k)+gamma(ig,k))*ua(ig,k)
	1163	s -(entr(ig,k)+gamma(ig,k))*u(ig,k)
	1164	s -wud(ig,k)+wud(ig,k+1))
	1165	s /masse(ig,k)
	1166	dv(ig,k)=((detr(ig,k)+gamma(ig,k))*va(ig,k)
	1167	s -(entr(ig,k)+gamma(ig,k))*v(ig,k)
	1168	s -wvd(ig,k)+wvd(ig,k+1))
	1169	s /masse(ig,k)
	1170	enddo
	1171	enddo
	1172
	1173	return
	1174	end
	1175	subroutine dqthermcell2(ngrid,nlay,ptimestep,fm,entr,masse,frac
	1176	. ,q,dq,qa)
[634]	1177	use dimphy
[542]	1178	implicit none
	1179
	1180	c=======================================================================
	1181	c
	1182	c Calcul du transport verticale dans la couche limite en presence
	1183	c de "thermiques" explicitement representes
	1184	c calcul du dq/dt une fois qu'on connait les ascendances
	1185	c
	1186	c=======================================================================
	1187
[634]	1188	cym#include "dimensions.h"
	1189	cym#include "dimphy.h"
[542]	1190
	1191	integer ngrid,nlay
	1192
	1193	real ptimestep
	1194	real masse(ngrid,nlay),fm(ngrid,nlay+1)
	1195	real entr(ngrid,nlay),frac(ngrid,nlay)
	1196	real q(ngrid,nlay)
	1197	real dq(ngrid,nlay)
	1198
	1199	real qa(klon,klev),detr(klon,klev),wqd(klon,klev+1)
	1200	real qe(klon,klev),zf,zf2
	1201
	1202	integer ig,k
	1203
	1204	c calcul du detrainement
	1205
	1206	do k=1,nlay
	1207	do ig=1,ngrid
	1208	detr(ig,k)=fm(ig,k)-fm(ig,k+1)+entr(ig,k)
	1209	enddo
	1210	enddo
	1211
	1212	c calcul de la valeur dans les ascendances
	1213	do ig=1,ngrid
	1214	qa(ig,1)=q(ig,1)
	1215	qe(ig,1)=q(ig,1)
	1216	enddo
	1217
	1218	do k=2,nlay
	1219	do ig=1,ngrid
	1220	if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt.
	1221	s 1.e-5*masse(ig,k)) then
	1222	zf=0.5*(frac(ig,k)+frac(ig,k+1))
	1223	zf2=1./(1.-zf)
	1224	qa(ig,k)=(fm(ig,k)qa(ig,k-1)+zf2entr(ig,k)*q(ig,k))
	1225	s /(fm(ig,k+1)+detr(ig,k)+entr(ig,k)zfzf2)
	1226	qe(ig,k)=(q(ig,k)-zfqa(ig,k))zf2
	1227	else
	1228	qa(ig,k)=q(ig,k)
	1229	qe(ig,k)=q(ig,k)
	1230	endif
	1231	enddo
	1232	enddo
	1233
	1234	do k=2,nlay
	1235	do ig=1,ngrid
	1236	c wqd(ig,k)=fm(ig,k)0.5(q(ig,k-1)+q(ig,k))
	1237	wqd(ig,k)=fm(ig,k)*qe(ig,k)
	1238	enddo
	1239	enddo
	1240	do ig=1,ngrid
	1241	wqd(ig,1)=0.
	1242	wqd(ig,nlay+1)=0.
	1243	enddo
	1244
	1245	do k=1,nlay
	1246	do ig=1,ngrid
	1247	dq(ig,k)=(detr(ig,k)qa(ig,k)-entr(ig,k)qe(ig,k)
	1248	s -wqd(ig,k)+wqd(ig,k+1))
	1249	s /masse(ig,k)
	1250	enddo
	1251	enddo
	1252
	1253	return
	1254	end
	1255	subroutine dvthermcell2(ngrid,nlay,ptimestep,fm,entr,masse
	1256	. ,fraca,larga
	1257	. ,u,v,du,dv,ua,va)
[634]	1258	use dimphy
[542]	1259	implicit none
	1260
	1261	c=======================================================================
	1262	c
	1263	c Calcul du transport verticale dans la couche limite en presence
	1264	c de "thermiques" explicitement representes
	1265	c calcul du dq/dt une fois qu'on connait les ascendances
	1266	c
	1267	c=======================================================================
	1268
[634]	1269	cym#include "dimensions.h"
	1270	cym#include "dimphy.h"
[542]	1271
	1272	integer ngrid,nlay
	1273
	1274	real ptimestep
	1275	real masse(ngrid,nlay),fm(ngrid,nlay+1)
	1276	real fraca(ngrid,nlay+1)
	1277	real larga(ngrid)
	1278	real entr(ngrid,nlay)
	1279	real u(ngrid,nlay)
	1280	real ua(ngrid,nlay)
	1281	real du(ngrid,nlay)
	1282	real v(ngrid,nlay)
	1283	real va(ngrid,nlay)
	1284	real dv(ngrid,nlay)
	1285
	1286	real qa(klon,klev),detr(klon,klev),zf,zf2
	1287	real wvd(klon,klev+1),wud(klon,klev+1)
	1288	real gamma0,gamma(klon,klev+1)
	1289	real ue(klon,klev),ve(klon,klev)
	1290	real dua,dva
	1291	integer iter
	1292
	1293	integer ig,k
	1294
	1295	c calcul du detrainement
	1296
	1297	do k=1,nlay
	1298	do ig=1,ngrid
	1299	detr(ig,k)=fm(ig,k)-fm(ig,k+1)+entr(ig,k)
	1300	enddo
	1301	enddo
	1302
	1303	c calcul de la valeur dans les ascendances
	1304	do ig=1,ngrid
	1305	ua(ig,1)=u(ig,1)
	1306	va(ig,1)=v(ig,1)
	1307	ue(ig,1)=u(ig,1)
	1308	ve(ig,1)=v(ig,1)
	1309	enddo
	1310
	1311	do k=2,nlay
	1312	do ig=1,ngrid
	1313	if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt.
	1314	s 1.e-5*masse(ig,k)) then
	1315	c On itère sur la valeur du coeff de freinage.
	1316	c gamma0=rho(ig,k)*(zlev(ig,k+1)-zlev(ig,k))
	1317	gamma0=masse(ig,k)
	1318	s sqrt( 0.5(fraca(ig,k+1)+fraca(ig,k)) )
	1319	s *0.5/larga(ig)
	1320	s *1.
	1321	c s *0.5
	1322	c gamma0=0.
	1323	zf=0.5*(fraca(ig,k)+fraca(ig,k+1))
	1324	zf=0.
	1325	zf2=1./(1.-zf)
	1326	c la première fois on multiplie le coefficient de freinage
	1327	c par le module du vent dans la couche en dessous.
	1328	dua=ua(ig,k-1)-u(ig,k-1)
	1329	dva=va(ig,k-1)-v(ig,k-1)
	1330	do iter=1,5
	1331	c On choisit une relaxation lineaire.
	1332	gamma(ig,k)=gamma0
	1333	c On choisit une relaxation quadratique.
	1334	gamma(ig,k)=gamma0sqrt(dua2+dva*2)
	1335	ua(ig,k)=(fm(ig,k)*ua(ig,k-1)
	1336	s +(zf2entr(ig,k)+gamma(ig,k))u(ig,k))
	1337	s /(fm(ig,k+1)+detr(ig,k)+entr(ig,k)zfzf2
	1338	s +gamma(ig,k))
	1339	va(ig,k)=(fm(ig,k)*va(ig,k-1)
	1340	s +(zf2entr(ig,k)+gamma(ig,k))v(ig,k))
	1341	s /(fm(ig,k+1)+detr(ig,k)+entr(ig,k)zfzf2
	1342	s +gamma(ig,k))
	1343	c print*,k,ua(ig,k),va(ig,k),u(ig,k),v(ig,k),dua,dva
	1344	dua=ua(ig,k)-u(ig,k)
	1345	dva=va(ig,k)-v(ig,k)
	1346	ue(ig,k)=(u(ig,k)-zfua(ig,k))zf2
	1347	ve(ig,k)=(v(ig,k)-zfva(ig,k))zf2
	1348	enddo
	1349	else
	1350	ua(ig,k)=u(ig,k)
	1351	va(ig,k)=v(ig,k)
	1352	ue(ig,k)=u(ig,k)
	1353	ve(ig,k)=v(ig,k)
	1354	gamma(ig,k)=0.
	1355	endif
	1356	enddo
	1357	enddo
	1358
	1359	do k=2,nlay
	1360	do ig=1,ngrid
	1361	wud(ig,k)=fm(ig,k)*ue(ig,k)
	1362	wvd(ig,k)=fm(ig,k)*ve(ig,k)
	1363	enddo
	1364	enddo
	1365	do ig=1,ngrid
	1366	wud(ig,1)=0.
	1367	wud(ig,nlay+1)=0.
	1368	wvd(ig,1)=0.
	1369	wvd(ig,nlay+1)=0.
	1370	enddo
	1371
	1372	do k=1,nlay
	1373	do ig=1,ngrid
	1374	du(ig,k)=((detr(ig,k)+gamma(ig,k))*ua(ig,k)
	1375	s -(entr(ig,k)+gamma(ig,k))*ue(ig,k)
	1376	s -wud(ig,k)+wud(ig,k+1))
	1377	s /masse(ig,k)
	1378	dv(ig,k)=((detr(ig,k)+gamma(ig,k))*va(ig,k)
	1379	s -(entr(ig,k)+gamma(ig,k))*ve(ig,k)
	1380	s -wvd(ig,k)+wvd(ig,k+1))
	1381	s /masse(ig,k)
	1382	enddo
	1383	enddo
	1384
	1385	return
	1386	end

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

Download in other formats:

Original Format