Context Navigation

lmdz_thermcell_flux2.f90 @ 5424

Last change on this file since 5424 was 5390, checked in by yann meurdesoif, 13 days ago

Remove UTF8 character that inihibit fortran parsing with GPU morphosis
Add missing END SUBROUTINE instead of simple END, that inhibit correct parsing with regulat expression parser (quick and dirty parsing)

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
[4590]	1	MODULE lmdz_thermcell_flux2
[1403]	2	!
	3	! $Id: lmdz_thermcell_flux2.f90 5390 2024-12-05 16:09:25Z fhourdin $
	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
[5390]	514	END SUBROUTINE thermcell_flux2
[4590]	515	END MODULE lmdz_thermcell_flux2

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

Download in other formats:

Original Format