Context Navigation

thermcell_flux2.F90 @ 5220

Last change on this file since 5220 was 2408, checked in by Laurent Fairhead, 9 years ago
Merged trunk changes r2298:2396 into testing branch
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 16.4 KB

Rev	Line
[1403]	1	!
	2	! $Id: thermcell_flux2.F90 2408 2015-12-14 10:43:09Z abarral $
	3	!
[972]	4	SUBROUTINE thermcell_flux2(ngrid,klev,ptimestep,masse, &
	5	& lalim,lmax,alim_star, &
	6	& entr_star,detr_star,f,rhobarz,zlev,zw2,fm,entr, &
	7	& detr,zqla,lev_out,lunout1,igout)
	8	!IM 060508 & detr,zqla,zmax,lev_out,lunout,igout)
	9
	10
	11	!---------------------------------------------------------------------------
	12	!thermcell_flux: deduction des flux
	13	!---------------------------------------------------------------------------
	14
[2408]	15	USE print_control_mod, ONLY: prt_level
[972]	16	IMPLICIT NONE
[1026]	17	#include "thermcell.h"
[972]	18
	19	INTEGER ig,l
	20	INTEGER ngrid,klev
	21
	22	REAL alim_star(ngrid,klev),entr_star(ngrid,klev)
	23	REAL detr_star(ngrid,klev)
	24	REAL zw2(ngrid,klev+1)
	25	REAL zlev(ngrid,klev+1)
	26	REAL masse(ngrid,klev)
	27	REAL ptimestep
	28	REAL rhobarz(ngrid,klev)
	29	REAL f(ngrid)
	30	INTEGER lmax(ngrid)
	31	INTEGER lalim(ngrid)
	32	REAL zqla(ngrid,klev)
	33	REAL zmax(ngrid)
	34
	35	integer ncorecfm1,ncorecfm2,ncorecfm3,ncorecalpha
	36	integer ncorecfm4,ncorecfm5,ncorecfm6,ncorecfm7,ncorecfm8
	37
	38
	39	REAL entr(ngrid,klev),detr(ngrid,klev)
	40	REAL fm(ngrid,klev+1)
	41	REAL zfm
	42
[1403]	43	integer igout,lout
[972]	44	integer lev_out
	45	integer lunout1
	46
	47	REAL f_old,ddd0,eee0,ddd,eee,zzz
	48
	49	REAL fomass_max,alphamax
	50	save fomass_max,alphamax
	51
[2408]	52	logical check_debug,labort_physic
[1403]	53
	54	character (len=20) :: modname='thermcell_flux2'
	55	character (len=80) :: abort_message
	56
[972]	57	fomass_max=0.5
	58	alphamax=0.7
	59
	60	ncorecfm1=0
	61	ncorecfm2=0
	62	ncorecfm3=0
	63	ncorecfm4=0
	64	ncorecfm5=0
	65	ncorecfm6=0
	66	ncorecfm7=0
	67	ncorecfm8=0
	68	ncorecalpha=0
	69
	70	!initialisation
	71	fm(:,:)=0.
	72
	73	if (prt_level.ge.10) then
	74	write(lunout1,*) 'Dans thermcell_flux 0'
	75	write(lunout1,*) 'flux base ',f(igout)
	76	write(lunout1,*) 'lmax ',lmax(igout)
	77	write(lunout1,*) 'lalim ',lalim(igout)
	78	write(lunout1,*) 'ig= ',igout
	79	write(lunout1,) ' l E A* D* '
	80	write(lunout1,'(i4,3e15.5)') (l,entr_star(igout,l),alim_star(igout,l),detr_star(igout,l) &
	81	& ,l=1,lmax(igout))
	82	endif
	83
	84
	85	!-------------------------------------------------------------------------
	86	! Verification de la nullite des entrainement et detrainement au dessus
	87	! de lmax(ig)
[1403]	88	! Active uniquement si check_debug=.true. ou prt_level>=10
[972]	89	!-------------------------------------------------------------------------
	90
[1403]	91	check_debug=.false..or.prt_level>=10
	92
	93	if (check_debug) then
[972]	94	do l=1,klev
	95	do ig=1,ngrid
	96	if (l.le.lmax(ig)) then
	97	if (entr_star(ig,l).gt.1.) then
	98	print*,'WARNING thermcell_flux 1 ig,l,lmax(ig)',ig,l,lmax(ig)
	99	print*,'entr_star(ig,l)',entr_star(ig,l)
	100	print*,'alim_star(ig,l)',alim_star(ig,l)
	101	print*,'detr_star(ig,l)',detr_star(ig,l)
	102	endif
	103	else
	104	if (abs(entr_star(ig,l))+abs(alim_star(ig,l))+abs(detr_star(ig,l)).gt.0.) then
	105	print*,'cas 1 : ig,l,lmax(ig)',ig,l,lmax(ig)
	106	print*,'entr_star(ig,l)',entr_star(ig,l)
	107	print*,'alim_star(ig,l)',alim_star(ig,l)
	108	print*,'detr_star(ig,l)',detr_star(ig,l)
[1403]	109	abort_message = ''
[2408]	110	labort_physic=.true.
	111	CALL abort_physic (modname,abort_message,1)
[972]	112	endif
	113	endif
	114	enddo
	115	enddo
[1403]	116	endif
[972]	117
	118	!-------------------------------------------------------------------------
	119	! Multiplication par le flux de masse issu de la femreture
	120	!-------------------------------------------------------------------------
	121
	122	do l=1,klev
	123	entr(:,l)=f(:)*(entr_star(:,l)+alim_star(:,l))
	124	detr(:,l)=f(:)*detr_star(:,l)
	125	enddo
	126
	127	if (prt_level.ge.10) then
	128	write(lunout1,*) 'Dans thermcell_flux 1'
	129	write(lunout1,*) 'flux base ',f(igout)
	130	write(lunout1,*) 'lmax ',lmax(igout)
	131	write(lunout1,*) 'lalim ',lalim(igout)
	132	write(lunout1,*) 'ig= ',igout
	133	write(lunout1,*) ' l E D W2'
	134	write(lunout1,'(i4,3e15.5)') (l,entr(igout,l),detr(igout,l) &
	135	& ,zw2(igout,l+1),l=1,lmax(igout))
	136	endif
	137
	138	fm(:,1)=0.
	139	do l=1,klev
	140	do ig=1,ngrid
	141	if (l.lt.lmax(ig)) then
	142	fm(ig,l+1)=fm(ig,l)+entr(ig,l)-detr(ig,l)
	143	elseif(l.eq.lmax(ig)) then
	144	fm(ig,l+1)=0.
	145	detr(ig,l)=fm(ig,l)+entr(ig,l)
	146	else
	147	fm(ig,l+1)=0.
	148	endif
	149	enddo
	150	enddo
	151
	152
	153
	154	! Test provisoire : pour comprendre pourquoi on corrige plein de fois
	155	! le cas fm6, on commence par regarder une premiere fois avant les
	156	! autres corrections.
	157
	158	do l=1,klev
	159	do ig=1,ngrid
	160	if (detr(ig,l).gt.fm(ig,l)) then
	161	ncorecfm8=ncorecfm8+1
	162	! igout=ig
	163	endif
	164	enddo
	165	enddo
	166
	167	! if (prt_level.ge.10) &
	168	! & call printflux(ngrid,klev,lunout1,igout,f,lmax,lalim, &
	169	! & ptimestep,masse,entr,detr,fm,'2 ')
	170
	171
	172
	173	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	174	! FH Version en cours de test;
	175	! par rapport a thermcell_flux, on fait une grande boucle sur "l"
[1146]	176	! et on modifie le flux avec tous les contr�les appliques d'affilee
[972]	177	! pour la meme couche
	178	! Momentanement, on duplique le calcule du flux pour pouvoir comparer
	179	! les flux avant et apres modif
	180	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	181
	182	do l=1,klev
	183
	184	do ig=1,ngrid
	185	if (l.lt.lmax(ig)) then
	186	fm(ig,l+1)=fm(ig,l)+entr(ig,l)-detr(ig,l)
	187	elseif(l.eq.lmax(ig)) then
	188	fm(ig,l+1)=0.
	189	detr(ig,l)=fm(ig,l)+entr(ig,l)
	190	else
	191	fm(ig,l+1)=0.
	192	endif
	193	enddo
	194
	195
	196	!-------------------------------------------------------------------------
	197	! Verification de la positivite des flux de masse
	198	!-------------------------------------------------------------------------
	199
	200	! do l=1,klev
	201	do ig=1,ngrid
	202	if (fm(ig,l+1).lt.0.) then
	203	! print*,'fm1<0',l+1,lmax(ig),fm(ig,l+1)
	204	ncorecfm1=ncorecfm1+1
	205	fm(ig,l+1)=fm(ig,l)
	206	detr(ig,l)=entr(ig,l)
	207	endif
	208	enddo
	209	! enddo
	210
	211	if (prt_level.ge.10) &
	212	& write(lunout1,'(i4,4e14.4)') l,masse(igout,l)/ptimestep, &
	213	& entr(igout,l),detr(igout,l),fm(igout,l+1)
	214
	215	!-------------------------------------------------------------------------
	216	!Test sur fraca croissant
	217	!-------------------------------------------------------------------------
[1026]	218	if (iflag_thermals_optflux==0) then
[972]	219	! do l=1,klev
	220	do ig=1,ngrid
	221	if (l.ge.lalim(ig).and.l.le.lmax(ig) &
	222	& .and.(zw2(ig,l+1).gt.1.e-10).and.(zw2(ig,l).gt.1.e-10) ) then
	223	! zzz est le flux en l+1 a frac constant
	224	zzz=fm(ig,l)rhobarz(ig,l+1)zw2(ig,l+1) &
	225	& /(rhobarz(ig,l)*zw2(ig,l))
	226	if (fm(ig,l+1).gt.zzz) then
	227	detr(ig,l)=detr(ig,l)+fm(ig,l+1)-zzz
	228	fm(ig,l+1)=zzz
	229	ncorecfm4=ncorecfm4+1
	230	endif
	231	endif
	232	enddo
	233	! enddo
	234	endif
	235
	236	if (prt_level.ge.10) &
	237	& write(lunout1,'(i4,4e14.4)') l,masse(igout,l)/ptimestep, &
	238	& entr(igout,l),detr(igout,l),fm(igout,l+1)
	239
	240
	241	!-------------------------------------------------------------------------
	242	!test sur flux de masse croissant
	243	!-------------------------------------------------------------------------
[1026]	244	if (iflag_thermals_optflux==0) then
[972]	245	! do l=1,klev
	246	do ig=1,ngrid
	247	if ((fm(ig,l+1).gt.fm(ig,l)).and.(l.gt.lalim(ig))) then
	248	f_old=fm(ig,l+1)
	249	fm(ig,l+1)=fm(ig,l)
	250	detr(ig,l)=detr(ig,l)+f_old-fm(ig,l+1)
	251	ncorecfm5=ncorecfm5+1
	252	endif
	253	enddo
	254	! enddo
	255	endif
	256
	257	if (prt_level.ge.10) &
	258	& write(lunout1,'(i4,4e14.4)') l,masse(igout,l)/ptimestep, &
	259	& entr(igout,l),detr(igout,l),fm(igout,l+1)
	260
	261	!fin 1.eq.0
	262	!-------------------------------------------------------------------------
	263	!detr ne peut pas etre superieur a fm
	264	!-------------------------------------------------------------------------
	265
	266	if(1.eq.1) then
	267
	268	! do l=1,klev
[1403]	269
	270
	271
[2408]	272	labort_physic=.false.
[972]	273	do ig=1,ngrid
	274	if (entr(ig,l)<0.) then
[2408]	275	labort_physic=.true.
[1403]	276	igout=ig
	277	lout=l
[972]	278	endif
[1403]	279	enddo
	280
[2408]	281	if (labort_physic) then
[1403]	282	print*,'N1 ig,l,entr',igout,lout,entr(igout,lout)
	283	abort_message = 'entr negatif'
[2408]	284	CALL abort_physic (modname,abort_message,1)
[1403]	285	endif
	286
	287	do ig=1,ngrid
[972]	288	if (detr(ig,l).gt.fm(ig,l)) then
	289	ncorecfm6=ncorecfm6+1
	290	detr(ig,l)=fm(ig,l)
	291	entr(ig,l)=fm(ig,l+1)
	292
	293	! Dans le cas ou on est au dessus de la couche d'alimentation et que le
	294	! detrainement est plus fort que le flux de masse, on stope le thermique.
	295	!test:on commente
	296	! if (l.gt.lalim(ig)) then
	297	! lmax(ig)=l
	298	! fm(ig,l+1)=0.
	299	! entr(ig,l)=0.
	300	! else
	301	! ncorecfm7=ncorecfm7+1
	302	! endif
	303	endif
	304
	305	if(l.gt.lmax(ig)) then
	306	detr(ig,l)=0.
	307	fm(ig,l+1)=0.
	308	entr(ig,l)=0.
	309	endif
[1403]	310	enddo
[972]	311
[2408]	312	labort_physic=.false.
[1403]	313	do ig=1,ngrid
[972]	314	if (entr(ig,l).lt.0.) then
[2408]	315	labort_physic=.true.
[1403]	316	igout=ig
[972]	317	endif
	318	enddo
[2408]	319	if (labort_physic) then
[1403]	320	ig=igout
	321	print*,'ig,l,lmax(ig)',ig,l,lmax(ig)
	322	print*,'entr(ig,l)',entr(ig,l)
	323	print*,'fm(ig,l)',fm(ig,l)
	324	abort_message = 'probleme dans thermcell flux'
[2408]	325	CALL abort_physic (modname,abort_message,1)
[1403]	326	endif
	327
	328
[972]	329	! enddo
	330	endif
	331
	332
	333	if (prt_level.ge.10) &
	334	& write(lunout1,'(i4,4e14.4)') l,masse(igout,l)/ptimestep, &
	335	& entr(igout,l),detr(igout,l),fm(igout,l+1)
	336
	337	!-------------------------------------------------------------------------
	338	!fm ne peut pas etre negatif
	339	!-------------------------------------------------------------------------
	340
	341	! do l=1,klev
	342	do ig=1,ngrid
	343	if (fm(ig,l+1).lt.0.) then
	344	detr(ig,l)=detr(ig,l)+fm(ig,l+1)
	345	fm(ig,l+1)=0.
	346	ncorecfm2=ncorecfm2+1
	347	endif
[1403]	348	enddo
	349
[2408]	350	labort_physic=.false.
[1403]	351	do ig=1,ngrid
[972]	352	if (detr(ig,l).lt.0.) then
[2408]	353	labort_physic=.true.
[1403]	354	igout=ig
	355	endif
	356	enddo
[2408]	357	if (labort_physic) then
[1403]	358	ig=igout
[972]	359	print*,'cas 2 : ig,l,lmax(ig)',ig,l,lmax(ig)
	360	print*,'detr(ig,l)',detr(ig,l)
	361	print*,'fm(ig,l)',fm(ig,l)
[1403]	362	abort_message = 'probleme dans thermcell flux'
[2408]	363	CALL abort_physic (modname,abort_message,1)
[1403]	364	endif
[972]	365	! enddo
	366
	367	if (prt_level.ge.10) &
	368	& write(lunout1,'(i4,4e14.4)') l,masse(igout,l)/ptimestep, &
	369	& entr(igout,l),detr(igout,l),fm(igout,l+1)
	370
	371	!-----------------------------------------------------------------------
	372	!la fraction couverte ne peut pas etre superieure a 1
	373	!-----------------------------------------------------------------------
	374
	375
	376	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	377	! FH Partie a revisiter.
	378	! Il semble qu'etaient codees ici deux optiques dans le cas
	379	! F/ (rho *w) > 1
	380	! soit limiter la hauteur du thermique en considerant que c'est
	381	! la derniere chouche, soit limiter F a rho w.
	382	! Dans le second cas, il faut en fait limiter a un peu moins
	383	! que ca parce qu'on a des 1 / ( 1 -alpha) un peu plus loin
	384	! dans thermcell_main et qu'il semble de toutes facons deraisonable
	385	! d'avoir des fractions de 1..
	386	! Ci dessous, et dans l'etat actuel, le premier des deux if est
	387	! sans doute inutile.
	388	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	389
	390	! do l=1,klev
	391	do ig=1,ngrid
	392	if (zw2(ig,l+1).gt.1.e-10) then
	393	zfm=rhobarz(ig,l+1)zw2(ig,l+1)alphamax
	394	if ( fm(ig,l+1) .gt. zfm) then
	395	f_old=fm(ig,l+1)
	396	fm(ig,l+1)=zfm
	397	! zw2(ig,l+1)=0.
	398	! zqla(ig,l+1)=0.
	399	detr(ig,l)=detr(ig,l)+f_old-fm(ig,l+1)
	400	! lmax(ig)=l+1
	401	! zmax(ig)=zlev(ig,lmax(ig))
	402	! print*,'alpha>1',l+1,lmax(ig)
	403	ncorecalpha=ncorecalpha+1
	404	endif
	405	endif
	406	enddo
	407	! enddo
	408	!
	409
	410
	411	if (prt_level.ge.10) &
	412	& write(lunout1,'(i4,4e14.4)') l,masse(igout,l)/ptimestep, &
	413	& entr(igout,l),detr(igout,l),fm(igout,l+1)
	414
	415	! Fin de la grande boucle sur les niveaux verticaux
	416	enddo
	417
	418	! if (prt_level.ge.10) &
	419	! & call printflux(ngrid,klev,lunout1,igout,f,lmax,lalim, &
	420	! & ptimestep,masse,entr,detr,fm,'8 ')
	421
	422
	423	!-----------------------------------------------------------------------
	424	! On fait en sorte que la quantite totale d'air entraine dans le
	425	! panache ne soit pas trop grande comparee a la masse de la maille
	426	!-----------------------------------------------------------------------
	427
	428	if (1.eq.1) then
[2408]	429	labort_physic=.false.
[972]	430	do l=1,klev-1
	431	do ig=1,ngrid
	432	eee0=entr(ig,l)
	433	ddd0=detr(ig,l)
	434	eee=entr(ig,l)-masse(ig,l)*fomass_max/ptimestep
	435	ddd=detr(ig,l)-eee
	436	if (eee.gt.0.) then
	437	ncorecfm3=ncorecfm3+1
	438	entr(ig,l)=entr(ig,l)-eee
	439	if ( ddd.gt.0.) then
	440	! l'entrainement est trop fort mais l'exces peut etre compense par une
	441	! diminution du detrainement)
	442	detr(ig,l)=ddd
	443	else
	444	! l'entrainement est trop fort mais l'exces doit etre compense en partie
	445	! par un entrainement plus fort dans la couche superieure
	446	if(l.eq.lmax(ig)) then
	447	detr(ig,l)=fm(ig,l)+entr(ig,l)
	448	else
	449	if(l.ge.lmax(ig).and.0.eq.1) then
[1403]	450	igout=ig
	451	lout=l
[2408]	452	labort_physic=.true.
[1403]	453	endif
	454	entr(ig,l+1)=entr(ig,l+1)-ddd
	455	detr(ig,l)=0.
	456	fm(ig,l+1)=fm(ig,l)+entr(ig,l)
	457	detr(ig,l)=0.
	458	endif
	459	endif
	460	endif
	461	enddo
	462	enddo
[2408]	463	if (labort_physic) then
[1403]	464	ig=igout
	465	l=lout
[972]	466	print*,'ig,l',ig,l
	467	print*,'eee0',eee0
	468	print*,'ddd0',ddd0
	469	print*,'eee',eee
	470	print*,'ddd',ddd
	471	print*,'entr',entr(ig,l)
	472	print*,'detr',detr(ig,l)
	473	print*,'masse',masse(ig,l)
	474	print*,'fomass_max',fomass_max
	475	print,'masse(ig,l)fomass_max/ptimestep',masse(ig,l)*fomass_max/ptimestep
	476	print*,'ptimestep',ptimestep
	477	print*,'lmax(ig)',lmax(ig)
	478	print*,'fm(ig,l+1)',fm(ig,l+1)
	479	print*,'fm(ig,l)',fm(ig,l)
[1403]	480	abort_message = 'probleme dans thermcell_flux'
[2408]	481	CALL abort_physic (modname,abort_message,1)
[972]	482	endif
[1403]	483	endif
[972]	484	!
	485	! ddd=detr(ig)-entre
	486	!on s assure que tout s annule bien en zmax
	487	do ig=1,ngrid
	488	fm(ig,lmax(ig)+1)=0.
	489	entr(ig,lmax(ig))=0.
	490	detr(ig,lmax(ig))=fm(ig,lmax(ig))+entr(ig,lmax(ig))
	491	enddo
	492
	493	!-----------------------------------------------------------------------
	494	! Impression du nombre de bidouilles qui ont ete necessaires
	495	!-----------------------------------------------------------------------
	496
	497	!IM 090508 beg
	498	! if (ncorecfm1+ncorecfm2+ncorecfm3+ncorecfm4+ncorecfm5+ncorecalpha > 0 ) then
	499	!
	500	! print*,'PB thermcell : on a du coriger ',ncorecfm1,'x fm1',&
	501	! & ncorecfm2,'x fm2',ncorecfm3,'x fm3 et', &
	502	! & ncorecfm4,'x fm4',ncorecfm5,'x fm5 et', &
	503	! & ncorecfm6,'x fm6', &
	504	! & ncorecfm7,'x fm7', &
	505	! & ncorecfm8,'x fm8', &
	506	! & ncorecalpha,'x alpha'
	507	! endif
	508	!IM 090508 end
	509
	510	! if (prt_level.ge.10) &
	511	! & call printflux(ngrid,klev,lunout1,igout,f,lmax,lalim, &
	512	! & ptimestep,masse,entr,detr,fm,'fin')
	513
	514
	515	return
	516	end

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

Download in other formats:

Original Format