Context Navigation

thermcell_flux.F90 @ 3562

Last change on this file since 3562 was 2647, checked in by dbardet, 3 years ago
Thermal plume model: Fix an index error for calculating the maximal authorized interlayer mass flux of the plume.
File size: 15.8 KB

Rev	Line
[2060]	1	!
	2	!
	3	!
[2101]	4	SUBROUTINE thermcell_flux(ngrid,nlay,ptimestep,masse, &
[2132]	5	lmin,lmax,entr_star,detr_star, &
[2143]	6	f,rhobarz,zlev,zw2,fm,entr,detr)
[2060]	7
	8
[2127]	9	!===============================================================================
[2143]	10	! Purpose: deduction des flux
[2127]	11	!
	12	! Modif 2019/04 (AB alexandre.boissinot@lmd.jussieu.fr)
	13	!
	14	!===============================================================================
[2060]	15
	16	USE print_control_mod, ONLY: prt_level
[2132]	17	USE thermcell_mod, ONLY: fomass_max, alpha_max
[2060]	18
	19	IMPLICIT NONE
	20
	21
[2127]	22	!===============================================================================
[2060]	23	! Declaration
[2127]	24	!===============================================================================
[2060]	25
[2127]	26	! Inputs:
	27	! -------
[2060]	28
[2177]	29	INTEGER, INTENT(in) :: ngrid
	30	INTEGER, INTENT(in) :: nlay
	31	INTEGER, INTENT(in) :: lmin(ngrid)
[2060]	32
[2177]	33	REAL, INTENT(in) :: entr_star(ngrid,nlay)
	34	REAL, INTENT(in) :: detr_star(ngrid,nlay)
	35	REAL, INTENT(in) :: zw2(ngrid,nlay+1)
	36	REAL, INTENT(in) :: zlev(ngrid,nlay+1)
	37	REAL, INTENT(in) :: masse(ngrid,nlay)
	38	REAL, INTENT(in) :: ptimestep
	39	REAL, INTENT(in) :: rhobarz(ngrid,nlay)
	40	REAL, INTENT(in) :: f(ngrid)
[2060]	41
[2127]	42	! Outputs:
[2177]	43	! --------
[2060]	44
[2177]	45	INTEGER, INTENT(inout) :: lmax(ngrid)
[2060]	46
[2177]	47	REAL, INTENT(out) :: entr(ngrid,nlay)
	48	REAL, INTENT(out) :: detr(ngrid,nlay)
	49	REAL, INTENT(out) :: fm(ngrid,nlay+1)
	50
[2127]	51	! Local:
	52	! ------
[2060]	53
[2127]	54	INTEGER ig, l, k
[2143]	55	INTEGER igout, lout ! Error grid point and level
[2060]	56
[2177]	57	REAL fmax ! Maximal authorized mass flux (alpha < alpha_max)
	58	REAL fff0 ! Save initial value of mass flux
	59	REAL emax ! Maximal authorized entrainment (entr*dt < mass_max)
[2127]	60	REAL eee0 ! Save initial value of entrainment
	61	REAL ddd0 ! Save initial value of detrainment
[2060]	62	REAL eee ! eee0 - layer mass * maximal authorized mass fraction / dt
	63	REAL ddd ! ddd0 - eee
[2143]	64	REAL fact
	65	REAL test
[2060]	66
[2143]	67	INTEGER ncorecentr
	68	INTEGER ncorecdetr
	69	INTEGER nerrorequa
	70	INTEGER ncorecfact
	71	INTEGER ncorecalpha
[2060]	72
	73	LOGICAL labort_physic
	74
[2127]	75	!===============================================================================
[2060]	76	! Initialization
[2127]	77	!===============================================================================
[2060]	78
[2143]	79	nerrorequa = 0
	80	ncorecentr = 0
	81	ncorecdetr = 0
	82	ncorecfact = 0
[2060]	83	ncorecalpha = 0
	84
	85	entr(:,:) = 0.
	86	detr(:,:) = 0.
	87	fm(:,:) = 0.
	88
[2143]	89	labort_physic = .false.
	90
	91	fact = 0.
	92
[2127]	93	!===============================================================================
[2143]	94	! Calcul de l'entrainement, du detrainement et du flux de masse
[2127]	95	!===============================================================================
[2060]	96
[2127]	97	!-------------------------------------------------------------------------------
[2060]	98	! Multiplication par la norme issue de la relation de fermeture
[2127]	99	!-------------------------------------------------------------------------------
[2060]	100
[2101]	101	DO l=1,nlay
[2127]	102	entr(:,l) = f(:) * entr_star(:,l)
[2060]	103	detr(:,l) = f(:) * detr_star(:,l)
	104	ENDDO
	105
[2127]	106	!-------------------------------------------------------------------------------
[2177]	107	! Mass flux
[2127]	108	!-------------------------------------------------------------------------------
[2060]	109
[2101]	110	DO l=1,nlay
[2060]	111	DO ig=1,ngrid
[2143]	112	IF (l < lmax(ig) .and. l >= lmin(ig)) THEN
[2060]	113	fm(ig,l+1) = fm(ig,l) + entr(ig,l) - detr(ig,l)
[2143]	114	ELSEIF (l == lmax(ig)) THEN
[2060]	115	fm(ig,l+1) = 0.
[2177]	116	entr(ig,l) = 0.
[2060]	117	detr(ig,l) = fm(ig,l) + entr(ig,l)
	118	ELSE
	119	fm(ig,l+1) = 0.
	120	entr(ig,l) = 0.
	121	detr(ig,l) = 0.
	122	ENDIF
	123	ENDDO
	124	ENDDO
	125
[2127]	126	!===============================================================================
[2143]	127	! Checking
[2127]	128	!===============================================================================
[2060]	129
[2101]	130	DO l=1,nlay
[2060]	131
[2127]	132	!-------------------------------------------------------------------------------
[2132]	133	! Is incoming mass flux positive ?
[2127]	134	!-------------------------------------------------------------------------------
[2060]	135
	136	DO ig=1,ngrid
[2143]	137	IF (fm(ig,l) < 0.) THEN
[2060]	138	labort_physic = .true.
	139	igout = ig
	140	lout = l
	141	ENDIF
	142	ENDDO
	143
[2127]	144	!-------------------------------------------------------------------------------
[2132]	145	! Is entrainment positive ?
[2127]	146	!-------------------------------------------------------------------------------
[2101]	147
[2060]	148	DO ig=1,ngrid
[2143]	149	IF (entr(ig,l) < 0.) THEN
[2060]	150	labort_physic = .true.
	151	igout = ig
	152	lout = l
	153	ENDIF
	154	ENDDO
	155
[2127]	156	!-------------------------------------------------------------------------------
[2132]	157	! Is detrainment positive ?
[2127]	158	!-------------------------------------------------------------------------------
[2103]	159
[2060]	160	DO ig=1,ngrid
[2143]	161	IF (detr(ig,l) < 0.) THEN
[2060]	162	labort_physic = .true.
	163	igout = ig
	164	lout = l
	165	ENDIF
	166	ENDDO
	167
[2143]	168	!-------------------------------------------------------------------------------
	169	! Abort
	170	!-------------------------------------------------------------------------------
	171
[2060]	172	IF (labort_physic) THEN
[2143]	173	print *, '---------------------------------------------------------'
	174	print *, 'ERROR: mass flux has negative value(s)!'
[2177]	175	print *, 'ig,l,norm', igout, lout, f(igout)
[2103]	176	print *, 'lmin,lmax', lmin(igout), lmax(igout)
[2143]	177	print *, '- - - - - - - - - - - - - - - - - - - - - - - - - - - - -'
[2127]	178	DO k=nlay,1,-1
[2480]	179	print *, 'fm,w ', fm(igout,k+1), zw2(igout,k+1)
[2127]	180	print *, 'entr,detr', entr(igout,k), detr(igout,k)
	181	ENDDO
[2480]	182	print *, 'fm,w ', fm(igout,1), zw2(igout,1)
[2143]	183	print *, '---------------------------------------------------------'
	184	CALL abort
[2060]	185	ENDIF
	186
[2127]	187	!-------------------------------------------------------------------------------
[2143]	188	! Is entrained mass lesser than fomass_max ?
[2127]	189	!-------------------------------------------------------------------------------
[2060]	190
[2127]	191	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[2060]	192	! AB : Entrainment is bigger than the maximal authorized value.
	193	! If we consider that the excess entrainement is in fact plume air which
	194	! is not detrained then we compensate it by decreasing detr.
	195	! If it's not enough, we can increase entr in the layer above and decrease
	196	! the outgoing mass flux in the current layer.
[2177]	197	! If it's still insufficient, code will abort (now commented).
[2103]	198	! Else we reset entr to its intial value and we renormalize entrainment,
	199	! detrainment and mass flux profiles such as we do not exceed the maximal
	200	! authorized entrained mass.
[2127]	201	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[2060]	202
	203	DO ig=1,ngrid
	204	eee0 = entr(ig,l)
	205	ddd0 = detr(ig,l)
[2177]	206	emax = masse(ig,l) * fomass_max / ptimestep
	207	IF (emax < 0.) THEN
	208	print *, 'ERROR: layer mass is negative!'
	209	print *, 'mass,emax', masse(ig,l), emax
	210	print *, 'ig,l', ig, l
	211	ENDIF
	212	IF (eee0 > emax) THEN
	213	entr(ig,l) = emax
	214	ddd = ddd0 + emax - eee0
[2132]	215	ncorecentr = ncorecentr + 1
[2143]	216	IF (ddd > 0.) THEN
[2060]	217	detr(ig,l) = ddd
[2143]	218	ELSEIF (l == lmax(ig)) THEN
[2103]	219	detr(ig,l) = fm(ig,l) + entr(ig,l)
[2143]	220	ELSEIF (entr(ig,l+1) > ABS(ddd)) THEN
[2103]	221	detr(ig,l) = 0.
	222	fm(ig,l+1) = fm(ig,l) + entr(ig,l)
	223	entr(ig,l+1) = entr(ig,l+1) + ddd
[2480]	224	ELSE
[2177]	225	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	226	! AB: Simulation abort (try to reduce the physical time step).
	227	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[2480]	228	! entr(ig,l) = entr(ig,l) + eee
	229	! igout = ig
	230	! lout = l
	231	! labort_physic = .true.
[2177]	232	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	233	! AB: We can renormalize the plume mass fluxes. I think it does not work.
	234	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	235	! fact = max(fact, eee0 / emax)
[2480]	236	fact = eee0 / emax
	237	entr(ig,l) = eee0
	238	ncorecfact = ncorecfact + 1
[2177]	239	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	240	! AB: The renormalization can be just applied in the local plume.
	241	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[2480]	242	DO k=lmin(ig),lmax(ig)
	243	entr(ig,k) = entr(ig,k) * emax / eee0
	244	detr(ig,k) = detr(ig,k) * emax / eee0
	245	fm(ig,k) = fm(ig,k) * emax / eee0
	246	ENDDO
[2060]	247	ENDIF
	248	ENDIF
	249	ENDDO
	250
[2143]	251	IF (labort_physic) THEN
	252	print *, '---------------------------------------------------------'
	253	print *, 'ERROR: Entrainment is greater than maximal authorized value!'
	254	print *, ' Nor detrainment neither entrainment can compensate it!'
	255	print *, 'ig,l,entr', igout, lout, entr(igout,lout)
	256	print *, 'lmin,lmax', lmin(igout), lmax(igout)
	257	print *, '- - - - - - - - - - - - - - - - - - - - - - - - - - - - -'
	258	print , 'e_max :', masse(igout,lout)fomass_max/ptimestep
	259	print *, ' fomass_max :', fomass_max
	260	print *, ' masse :', masse(igout,lout)
	261	print *, ' ptimestep :', ptimestep
	262	print *, 'norm :', f(igout)
	263	print , 'entr :', entr_star(igout,lout)
	264	print *, '- - - - - - - - - - - - - - - - - - - - - - - - - - - - -'
	265	DO k=nlay,1,-1
[2480]	266	print *, 'fm,w ', fm(igout,k+1), zw2(igout,k+1)
[2143]	267	print *, 'entr,detr', entr(igout,k), detr(igout,k)
	268	ENDDO
[2480]	269	print *, 'fm,w ', fm(igout,1), zw2(igout,1)
[2143]	270	print *, '---------------------------------------------------------'
	271	CALL abort
	272	ENDIF
[2060]	273
[2127]	274	!-------------------------------------------------------------------------------
[2177]	275	! Is updraft fraction lesser than alpha_max ?
[2127]	276	!-------------------------------------------------------------------------------
[2060]	277
	278	DO ig=1,ngrid
[2177]	279	fff0 = fm(ig,l+1)
[2647]	280	fmax = rhobarz(ig,l) * zw2(ig,l+1) * alpha_max
	281	! rhobarz is an arithmetic mean of the mid-layer rho, and thus evaluated
	282	! at the interlayer
	283	! so rhobarz(ig,l) can be taken to evaluate the inter-layer maximal
	284	! authorized mass flux fmax (with the plume vertical velocity zw2).
[2177]	285	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	286	! AB: The plume mass flux can be reduced.
	287	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	288	! IF (fff0 > fmax) THEN
	289	! fm(ig,l+1) = fmax
	290	! detr(ig,l) = detr(ig,l) + fff0 - fmax
	291	! ncorecalpha = ncorecalpha + 1
	292	! entr(ig,l+1) = entr(ig,l+1) + fff0 - fmax
	293	! ENDIF
	294	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	295	! AB: The plume can be stopped here.
	296	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	297	IF (fff0 > fmax) THEN
[2143]	298	ncorecalpha = ncorecalpha + 1
[2177]	299	DO k=l+1,lmax(ig)
	300	entr(ig,k) = 0.
	301	detr(ig,k) = 0.
	302	fm(ig,k) = 0.
	303	ENDDO
	304	lmax(ig) = l
	305	entr(ig,l) = 0.
	306	detr(ig,l) = fm(ig,l)
	307	ENDIF
	308	ENDDO
	309
	310	!-------------------------------------------------------------------------------
	311	! Is detrainment lesser than incoming mass flux ?
	312	!-------------------------------------------------------------------------------
	313
[2127]	314	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[2177]	315	! AB : Even if fm has no negative value, it can be lesser than detr.
	316	! That is not suitable because when we will mix the plume with the
	317	! environment, it will detrain more mass than it is physically able to do.
	318	! When it occures, that imply that entr + fm is greater than detr,
	319	! otherwise we get a negative outgoing mass flux (cf. above).
	320	! That is why we decrease entrainment and detrainment as follows.
[2127]	321	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[2177]	322
	323	DO ig=1,ngrid
	324	IF (detr(ig,l) > fm(ig,l)) THEN
	325	detr(ig,l) = fm(ig,l)
	326	entr(ig,l) = fm(ig,l+1)
	327	ncorecdetr = ncorecdetr + 1
[2060]	328	ENDIF
	329	ENDDO
	330
	331	ENDDO
	332
[2177]	333	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	334	! AB: The renormalization can be applied everywhere.
	335	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	336	! IF (fact > 0.) THEN
	337	! entr(:,:) = entr(:,:) / fact
	338	! detr(:,:) = detr(:,:) / fact
	339	! fm(:,:) = fm(:,:) / fact
	340	! ENDIF
[2143]	341
[2127]	342	!-------------------------------------------------------------------------------
[2143]	343	! Is equation df/dz = e - d still verified ?
[2127]	344	!-------------------------------------------------------------------------------
[2060]	345
[2143]	346	! DO l=1,nlay
	347	! DO ig=1,ngrid
	348	! test = abs(fm(ig,l) + entr(ig,l) - detr(ig,l) - fm(ig,l+1))
	349	! IF (test > 1.e-10) THEN
	350	! nerrorequa = nerrorequa + 1
	351	! ENDIF
	352	! ENDDO
	353	! ENDDO
	354
	355	!-------------------------------------------------------------------------------
	356	! Reset top boundary conditions
	357	!-------------------------------------------------------------------------------
	358
[2103]	359	DO ig=1,ngrid
[2177]	360	IF (lmax(ig) > 0) THEN
	361	detr(ig,lmax(ig)) = fm(ig,lmax(ig))
[2060]	362	fm(ig,lmax(ig)+1) = 0.
	363	entr(ig,lmax(ig)) = 0.
[2103]	364	ENDIF
	365	ENDDO
[2060]	366
[2143]	367	!===============================================================================
	368	! Outputs
	369	!===============================================================================
[2060]	370
[2143]	371	IF (prt_level > 0) THEN
	372
[2177]	373	IF (ncorecdetr > 0) THEN
[2060]	374	print *, 'WARNING: Detrainment is greater than mass flux!'
[2143]	375	print *, 'In', ncorecdetr, 'grid point(s) over', nlay, 'x', ngrid
[2060]	376	ENDIF
	377
[2177]	378	IF (ncorecentr > 0) THEN
[2103]	379	print *, 'WARNING: Entrained mass is greater than maximal authorized value!'
[2177]	380	print *, 'In', ncorecentr, 'grid point(s) over', nlay, 'x', ngrid
[2060]	381	ENDIF
	382
[2177]	383	IF (ncorecfact > 0) THEN
[2143]	384	print *, 'WARNING: Entrained mass needs renormalization to be ok!'
[2177]	385	print *, 'In', ncorecfact, 'grid point(s) over', nlay, 'x', ngrid
	386	! print *, 'WARNING: Entr fact:', fact
[2143]	387	ENDIF
	388
[2177]	389	! IF (nerrorequa > 0) THEN
[2143]	390	! print *, 'WARNING: !'
	391	! print *, 'in', nerrorequa, 'grid point(s) over', nlay, 'x', ngrid
	392	! ENDIF
	393
[2177]	394	IF (ncorecalpha > 0) THEN
[2103]	395	print *, 'WARNING: Updraft fraction is greater than maximal authorized value!'
[2177]	396	print *, 'In', ncorecalpha, 'grid point(s) over', nlay, 'x', ngrid
[2060]	397	ENDIF
	398
	399	ENDIF
	400
	401
	402	RETURN
	403	END

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

Download in other formats:

Original Format