Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

lmdz_thermcell_flux2.f90 @ 5322

Last change on this file since 5322 was 5268, checked in by abarral, 3 weeks ago
.f90 <-> .F90 depending on cpp key use
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
[4590]	1	MODULE lmdz_thermcell_flux2
[1403]	2	!
	3	! $Id: lmdz_thermcell_flux2.f90 5268 2024-10-23 17:02:39Z abarral $
	4	!
[4590]	5	CONTAINS
	6
[4133]	7	SUBROUTINE thermcell_flux2(ngrid,nlay,ptimestep,masse, &
[972]	8	& lalim,lmax,alim_star, &
	9	& entr_star,detr_star,f,rhobarz,zlev,zw2,fm,entr, &
	10	& detr,zqla,lev_out,lunout1,igout)
	11	!IM 060508 & detr,zqla,zmax,lev_out,lunout,igout)
	12
	13
	14	!---------------------------------------------------------------------------
	15	!thermcell_flux: deduction des flux
	16	!---------------------------------------------------------------------------
	17
[4590]	18	USE lmdz_thermcell_ini, ONLY : prt_level,iflag_thermals_optflux
[972]	19	IMPLICIT NONE
	20
[4133]	21	! arguments
	22	INTEGER, intent(in) :: ngrid,nlay
	23	REAL, intent(in) :: ptimestep
	24	REAL, intent(in), dimension(ngrid,nlay) :: masse
	25	INTEGER, intent(in), dimension(ngrid) :: lalim,lmax
	26	REAL, intent(in), dimension(ngrid,nlay) :: alim_star,entr_star,detr_star
	27	REAL, intent(in), dimension(ngrid) :: f
	28	REAL, intent(in), dimension(ngrid,nlay) :: rhobarz
	29	REAL, intent(in), dimension(ngrid,nlay+1) :: zw2,zlev
	30	! FH : laisser ca le temps de verifier qu'on a bien fait de commenter les
	31	! lignes faisant apparaitre zqla, zmax ...
	32	! REAL, intent(in), dimension(ngrid) :: zmax(ngrid)
	33	! enlever aussi zqla
	34	REAL, intent(in), dimension(ngrid,nlay) :: zqla ! not used
	35	integer, intent(in) :: lev_out, lunout1
	36
	37	REAL,intent(out), dimension(ngrid,nlay) :: entr,detr
	38	REAL,intent(out), dimension(ngrid,nlay+1) :: fm
	39
	40	! local
[972]	41	INTEGER ig,l
[4133]	42	integer igout,lout
	43	REAL zfm
[972]	44	integer ncorecfm1,ncorecfm2,ncorecfm3,ncorecalpha
	45	integer ncorecfm4,ncorecfm5,ncorecfm6,ncorecfm7,ncorecfm8
	46
	47
	48	REAL f_old,ddd0,eee0,ddd,eee,zzz
	49
[3035]	50	REAL,SAVE :: fomass_max=0.5
	51	REAL,SAVE :: alphamax=0.7
	52	!$OMP THREADPRIVATE(fomass_max,alphamax)
[972]	53
[2311]	54	logical check_debug,labort_physic
[1403]	55
	56	character (len=20) :: modname='thermcell_flux2'
	57	character (len=80) :: abort_message
	58
[972]	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
[4133]	94	do l=1,nlay
[972]	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 = ''
[2311]	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
[4133]	122	do l=1,nlay
[972]	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.
[4133]	139	do l=1,nlay
[972]	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
[4133]	158	do l=1,nlay
[972]	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) &
[4133]	168	! & call printflux(ngrid,nlay,lunout1,igout,f,lmax,lalim, &
[972]	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
[4133]	182	do l=1,nlay
[972]	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
[4133]	200	! do l=1,nlay
[972]	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
[4133]	219	! do l=1,nlay
[972]	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
[4133]	245	! do l=1,nlay
[972]	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
[4133]	268	! do l=1,nlay
[1403]	269
	270
	271
[2311]	272	labort_physic=.false.
[972]	273	do ig=1,ngrid
	274	if (entr(ig,l)<0.) then
[2311]	275	labort_physic=.true.
[1403]	276	igout=ig
	277	lout=l
[972]	278	endif
[1403]	279	enddo
	280
[2311]	281	if (labort_physic) then
[1403]	282	print*,'N1 ig,l,entr',igout,lout,entr(igout,lout)
	283	abort_message = 'entr negatif'
[2311]	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
[2311]	312	labort_physic=.false.
[1403]	313	do ig=1,ngrid
[972]	314	if (entr(ig,l).lt.0.) then
[2311]	315	labort_physic=.true.
[1403]	316	igout=ig
[972]	317	endif
	318	enddo
[2311]	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'
[2311]	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
[4133]	341	! do l=1,nlay
[972]	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
[2311]	350	labort_physic=.false.
[1403]	351	do ig=1,ngrid
[972]	352	if (detr(ig,l).lt.0.) then
[2311]	353	labort_physic=.true.
[1403]	354	igout=ig
	355	endif
	356	enddo
[2311]	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'
[2311]	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
[4133]	390	! do l=1,nlay
[972]	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	detr(ig,l)=detr(ig,l)+f_old-fm(ig,l+1)
	398	! lmax(ig)=l+1
	399	! zmax(ig)=zlev(ig,lmax(ig))
	400	! print*,'alpha>1',l+1,lmax(ig)
	401	ncorecalpha=ncorecalpha+1
	402	endif
	403	endif
	404	enddo
	405	! enddo
	406	!
	407
	408
	409	if (prt_level.ge.10) &
	410	& write(lunout1,'(i4,4e14.4)') l,masse(igout,l)/ptimestep, &
	411	& entr(igout,l),detr(igout,l),fm(igout,l+1)
	412
	413	! Fin de la grande boucle sur les niveaux verticaux
	414	enddo
	415
	416	! if (prt_level.ge.10) &
[4133]	417	! & call printflux(ngrid,nlay,lunout1,igout,f,lmax,lalim, &
[972]	418	! & ptimestep,masse,entr,detr,fm,'8 ')
	419
	420
	421	!-----------------------------------------------------------------------
	422	! On fait en sorte que la quantite totale d'air entraine dans le
	423	! panache ne soit pas trop grande comparee a la masse de la maille
	424	!-----------------------------------------------------------------------
	425
	426	if (1.eq.1) then
[2311]	427	labort_physic=.false.
[4133]	428	do l=1,nlay-1
[972]	429	do ig=1,ngrid
	430	eee0=entr(ig,l)
	431	ddd0=detr(ig,l)
	432	eee=entr(ig,l)-masse(ig,l)*fomass_max/ptimestep
	433	ddd=detr(ig,l)-eee
	434	if (eee.gt.0.) then
	435	ncorecfm3=ncorecfm3+1
	436	entr(ig,l)=entr(ig,l)-eee
	437	if ( ddd.gt.0.) then
	438	! l'entrainement est trop fort mais l'exces peut etre compense par une
	439	! diminution du detrainement)
	440	detr(ig,l)=ddd
	441	else
	442	! l'entrainement est trop fort mais l'exces doit etre compense en partie
	443	! par un entrainement plus fort dans la couche superieure
	444	if(l.eq.lmax(ig)) then
	445	detr(ig,l)=fm(ig,l)+entr(ig,l)
	446	else
	447	if(l.ge.lmax(ig).and.0.eq.1) then
[1403]	448	igout=ig
	449	lout=l
[2311]	450	labort_physic=.true.
[1403]	451	endif
	452	entr(ig,l+1)=entr(ig,l+1)-ddd
	453	detr(ig,l)=0.
	454	fm(ig,l+1)=fm(ig,l)+entr(ig,l)
	455	detr(ig,l)=0.
	456	endif
	457	endif
	458	endif
	459	enddo
	460	enddo
[2311]	461	if (labort_physic) then
[1403]	462	ig=igout
	463	l=lout
[972]	464	print*,'ig,l',ig,l
	465	print*,'eee0',eee0
	466	print*,'ddd0',ddd0
	467	print*,'eee',eee
	468	print*,'ddd',ddd
	469	print*,'entr',entr(ig,l)
	470	print*,'detr',detr(ig,l)
	471	print*,'masse',masse(ig,l)
	472	print*,'fomass_max',fomass_max
	473	print,'masse(ig,l)fomass_max/ptimestep',masse(ig,l)*fomass_max/ptimestep
	474	print*,'ptimestep',ptimestep
	475	print*,'lmax(ig)',lmax(ig)
	476	print*,'fm(ig,l+1)',fm(ig,l+1)
	477	print*,'fm(ig,l)',fm(ig,l)
[1403]	478	abort_message = 'probleme dans thermcell_flux'
[2311]	479	CALL abort_physic (modname,abort_message,1)
[972]	480	endif
[1403]	481	endif
[972]	482	!
	483	! ddd=detr(ig)-entre
	484	!on s assure que tout s annule bien en zmax
	485	do ig=1,ngrid
	486	fm(ig,lmax(ig)+1)=0.
	487	entr(ig,lmax(ig))=0.
	488	detr(ig,lmax(ig))=fm(ig,lmax(ig))+entr(ig,lmax(ig))
	489	enddo
	490
	491	!-----------------------------------------------------------------------
	492	! Impression du nombre de bidouilles qui ont ete necessaires
	493	!-----------------------------------------------------------------------
	494
	495	!IM 090508 beg
	496	! if (ncorecfm1+ncorecfm2+ncorecfm3+ncorecfm4+ncorecfm5+ncorecalpha > 0 ) then
	497	!
	498	! print*,'PB thermcell : on a du coriger ',ncorecfm1,'x fm1',&
	499	! & ncorecfm2,'x fm2',ncorecfm3,'x fm3 et', &
	500	! & ncorecfm4,'x fm4',ncorecfm5,'x fm5 et', &
	501	! & ncorecfm6,'x fm6', &
	502	! & ncorecfm7,'x fm7', &
	503	! & ncorecfm8,'x fm8', &
	504	! & ncorecalpha,'x alpha'
	505	! endif
	506	!IM 090508 end
	507
	508	! if (prt_level.ge.10) &
[4133]	509	! & call printflux(ngrid,nlay,lunout1,igout,f,lmax,lalim, &
[972]	510	! & ptimestep,masse,entr,detr,fm,'fin')
	511
	512
[4133]	513	RETURN
[972]	514	end
[4590]	515	END MODULE lmdz_thermcell_flux2

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

Download in other formats:

Original Format