Context Navigation

thermcell_flux2.F90 @ 3626

Last change on this file since 3626 was 3102, checked in by oboucher, 7 years ago
Removing x permission from these files
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:keywords set to `Author Date Id Revision`
File size: 16.4 KB

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

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

Download in other formats:

Original Format