Context Navigation

lmdz_thermcell_flux2.f90 @ 5522

Last change on this file since 5522 was 5512, checked in by yann meurdesoif, 8 days ago

Implement GPU automatic port for :

Thermics
acama_gwd_rando
flott_gwd_rando

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

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

Download in other formats:

Original Format