Context Navigation

thermcell_flux2.F90 @ 4296

Last change on this file since 4296 was 4133, checked in by fhourdin, 3 years ago
Corrections thermiques pour replay
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.7 KB

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

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

Download in other formats:

Original Format