Context Navigation

thermcell_flux.F90 @ 4799

Last change on this file since 4799 was 2311, checked in by Ehouarn Millour, 10 years ago

Further modifications to enforce physics/dynamics separation:

moved iniprint.h and misc_mod back to dyn3d_common, as these should only be used by dynamics.
created print_control_mod in the physics to store flags prt_level, lunout, debug to be local to physics (should be used rather than iniprint.h)
created abort_physic.F90 , which does the same job as abort_gcm() did, but should be used instead when in physics.
reactivated inifis (turned it into a module, inifis_mod.F90) to initialize physical constants and print_control_mod flags.

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.8 KB

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

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

Download in other formats:

Original Format