Context Navigation

thermcell.F @ 1160

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

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

Download in other formats:

Original Format