Context Navigation

thermcell_plume.F90 @ 2176

Last change on this file since 2176 was 2145, checked in by aboissinot, 5 years ago
Some temporary outputs in thermcell_closure are removed and a bug in thermcell_plume is fixed.
File size: 14.5 KB

Rev	Line
[2060]	1	!
	2	!
	3	!
[2127]	4	SUBROUTINE thermcell_plume(ngrid,nlay,nq,ptimestep,ztv, &
	5	zhl,zqt,zql,rhobarz,zlev,pplev,pphi,zpopsk, &
	6	detr_star,entr_star,f_star, &
[2143]	7	ztva,zhla,zqla,zqta,zta,zqsa, &
	8	zw2,lmix,lmin)
[2060]	9
	10
[2127]	11	!===============================================================================
	12	! Purpose: calcule les valeurs de qt, thetal et w dans l ascendance
	13	!
	14	! Nota Bene
	15	! ql means "non-gaseous water mass mixing ratio" (liquid and solid)
	16	! qv means "vapor mass mixing ratio"
	17	! qt means "total water mass mixing ratio"
	18	! TP means "potential temperature"
	19	! TRPV means "virtual potential temperature with latent heat release"
	20	! TPV means "virtual potential temperature"
	21	! TR means "temperature with latent heat release"
	22	!===============================================================================
[2060]	23
	24	USE print_control_mod, ONLY: prt_level
[2071]	25	USE watercommon_h, ONLY: RLvCp, RETV, Psat_water
[2127]	26	USE tracer_h, ONLY: igcm_h2o_vap
[2060]	27	USE thermcell_mod
	28
	29	IMPLICIT NONE
	30
	31
[2127]	32	!===============================================================================
[2060]	33	! Declaration
[2127]	34	!===============================================================================
[2060]	35
[2113]	36	! Inputs:
	37	! -------
[2060]	38
[2127]	39	INTEGER ngrid, nlay, nq
[2060]	40
[2127]	41	REAL ptimestep
	42	REAL rhobarz(ngrid,nlay) ! Levels density
	43	REAL zlev(ngrid,nlay+1) ! Levels altitude
	44	REAL pplev(ngrid,nlay+1) ! Levels pressure
	45	REAL pphi(ngrid,nlay) ! Geopotential
	46	REAL zpopsk(ngrid,nlay) ! Exner function
	47
[2101]	48	REAL ztv(ngrid,nlay) ! TRPV environment
[2127]	49	REAL zhl(ngrid,nlay) ! TP environment
	50	REAL zqt(ngrid,nlay) ! qt environment
	51	REAL zql(ngrid,nlay) ! ql environment
[2060]	52
[2113]	53	! Outputs:
	54	! --------
[2060]	55
	56	INTEGER lmin(ngrid) ! plume base level (first unstable level)
	57	INTEGER lmix(ngrid) ! maximum vertical speed level
	58
[2101]	59	REAL detr_star(ngrid,nlay) ! normalized detrainment
	60	REAL entr_star(ngrid,nlay) ! normalized entrainment
	61	REAL f_star(ngrid,nlay+1) ! normalized mass flux
[2060]	62
[2101]	63	REAL ztva(ngrid,nlay) ! TRPV plume (after mixing)
[2143]	64	REAL zhla(ngrid,nlay) ! TP plume ?
[2101]	65	REAL zqla(ngrid,nlay) ! ql plume (after mixing)
[2143]	66	REAL zqta(ngrid,nlay) ! qt plume ?
[2127]	67	REAL zqsa(ngrid,nlay) ! qsat plume (after mixing)
[2143]	68	REAL zw2(ngrid,nlay+1) ! w plume (after mixing)
[2060]	69
[2113]	70	! Local:
	71	! ------
[2060]	72
	73	INTEGER ig, l, k
	74
[2127]	75	REAL ztva_est(ngrid,nlay) ! TRPV plume (before mixing)
[2101]	76	REAL zqla_est(ngrid,nlay) ! ql plume (before mixing)
	77	REAL zta_est(ngrid,nlay) ! TR plume (before mixing)
[2127]	78	REAL zqsa_est(ngrid) ! qsat plume (before mixing)
[2143]	79	REAL zw2_est(ngrid,nlay+1) ! w plume (before mixing)
[2060]	80
[2127]	81	REAL zta(ngrid,nlay) ! TR plume (after mixing)
[2060]	82
[2127]	83	REAL zbuoy(ngrid,nlay) ! Plume buoyancy
	84	REAL ztemp(ngrid) ! Temperature for saturation vapor pressure computation in plume
	85	REAL zdz ! Layers heights
	86	REAL ztv2(ngrid,nlay) ! ztv + d_temp * Dirac(l=linf)
	87
[2060]	88	REAL zbetalpha !
	89	REAL zdw2 !
	90	REAL zdw2bis !
	91	REAL zw2fact !
[2127]	92	REAL zw2m ! Average vertical velocity between two successive levels
	93	REAL gamma ! Plume acceleration term (to compute vertical velocity)
	94	REAL test ! Test to know how to compute entrainment and detrainment
	95
	96	REAL psat ! Dummy argument for Psat_water()
[2060]	97
[2127]	98	LOGICAL active(ngrid) ! If the plume is active at ig (speed and incoming mass flux > 0 or l=lmin)
	99	LOGICAL activetmp(ngrid) ! If the plume is active at ig (active=true and outgoing mass flux > 0)
[2071]	100
[2127]	101	!===============================================================================
[2060]	102	! Initialization
[2127]	103	!===============================================================================
[2060]	104
	105	zbetalpha = betalpha / (1. + betalpha)
	106
[2127]	107	ztva(:,:) = ztv(:,:) ! ztva is set to TPV environment
	108	zhla(:,:) = zhl(:,:) ! zhla is set to TP environment
	109	zqta(:,:) = zqt(:,:) ! zqta is set to qt environment
	110	zqla(:,:) = zql(:,:) ! zqla is set to ql environment
[2060]	111
[2127]	112	zqsa_est(:) = 0.
	113	zqsa(:,:) = 0.
[2060]	114
[2127]	115	zw2_est(:,:) = 0.
[2060]	116	zw2(:,:) = 0.
	117
	118	zbuoy(:,:) = 0.
	119
	120	f_star(:,:) = 0.
	121	detr_star(:,:) = 0.
	122	entr_star(:,:) = 0.
	123
	124	lmix(:) = 1
[2127]	125	lmin(:) = 1
[2060]	126
[2093]	127	ztv2(:,:) = ztv(:,:)
[2101]	128	ztv2(:,linf) = ztv(:,linf) + d_temp
[2093]	129
[2143]	130	active(:) = .false.
	131
[2127]	132	!===============================================================================
	133	! First layer computation
	134	!===============================================================================
[2060]	135
	136	DO ig=1,ngrid
	137	l = linf
[2143]	138	DO WHILE (.not.active(ig).and.(pplev(ig,l+1) > pres_limit).and.(l < nlay))
	139	zbuoy(ig,l) = RG * (ztv2(ig,l) - ztv2(ig,l+1)) / ztv2(ig,l+1)
	140	zdz = zlev(ig,l+1) - zlev(ig,l)
	141	zw2m = afact * zbuoy(ig,l) * zdz / (1. + betalpha)
	142	! gamma = afact * zbuoy(ig,l) - fact_epsilon * zw2m
	143	! test = gamma / zw2m - nu
	144	test = zbuoy(ig,l)
	145	IF (test > 0.) THEN
[2060]	146	lmin(ig) = l
[2143]	147	! entr_star(ig,l) = zdz * f_star(ig,l) * zbetalpha * gamma / zw2m - nu ! Problem because f*(ig,l) = 0
	148	! detr_star(ig,l) = f_star(ig,l) * nu ! Problem because f*(ig,l) = 0
	149	! f_star(ig,l+1) = entr_star(ig,l) - detr_star(ig,l)
	150	entr_star(ig,l) = 1.
	151	f_star(ig,l+1) = 1.
	152	zw2_est(ig,l+1) = zw2m * 2.
	153	zw2(ig,l+1) = zw2_est(ig,l+1)
[2145]	154	active(ig) = .true.
[2060]	155	ENDIF
	156	l = l + 1
	157	ENDDO
	158	ENDDO
	159
[2127]	160	!===============================================================================
	161	! Thermal plumes computations
	162	!===============================================================================
[2060]	163
[2101]	164	DO l=2,nlay-1
[2060]	165
[2127]	166	!-------------------------------------------------------------------------------
[2143]	167	! Is thermal plume (still) active ?
[2127]	168	!-------------------------------------------------------------------------------
	169
[2060]	170	DO ig=1,ngrid
[2143]	171	active(ig) = (active(ig).or.(l == lmin(ig)+1)) &
	172	& .and.(zw2(ig,l) > 1.e-10) &
	173	& .and.(f_star(ig,l) > 1.e-10)
[2060]	174	ENDDO
	175
[2143]	176	!-------------------------------------------------------------------------------
	177	! Latent heat release (before mixing)
	178	!-------------------------------------------------------------------------------
	179
[2127]	180	ztemp(:) = zpopsk(:,l) * zhla(:,l-1)
[2065]	181
[2060]	182	DO ig=1,ngrid
[2127]	183	CALL Psat_water(ztemp(ig), pplev(ig,l), psat, zqsa_est(ig))
[2060]	184	ENDDO
	185
[2127]	186	!-------------------------------------------------------------------------------
[2143]	187	! Vertical speed (before mixing)
[2127]	188	!-------------------------------------------------------------------------------
[2060]	189
	190	DO ig=1,ngrid
	191	IF (active(ig)) THEN
[2143]	192	zqla_est(ig,l) = MAX(0.,zqta(ig,l-1) - zqsa_est(ig)) ! zqla_est set to ql plume
[2127]	193	zta_est(ig,l) = zhla(ig,l-1) * zpopsk(ig,l) & ! zta_est set to TR plume
	194	& + RLvCp * zqla_est(ig,l)
	195	ztva_est(ig,l) = zta_est(ig,l) / zpopsk(ig,l) & ! ztva_est set to TRPV plume
	196	& * (1. + RETV * (zqta(ig,l-1)-zqla_est(ig,l)) - zqla_est(ig,l))
[2060]	197
[2127]	198	zbuoy(ig,l) = RG * (ztva_est(ig,l) - ztv(ig,l)) / ztv(ig,l)
[2060]	199	zdz = zlev(ig,l+1) - zlev(ig,l)
	200
[2127]	201	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	202	! AB: initial formulae
	203	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[2060]	204	! zw2fact = fact_epsilon * 2. * zdz / (1. + betalpha)
	205	! zdw2 = afact * zbuoy(ig,l) / fact_epsilon
	206	! zdw2bis = afact * zbuoy(ig,l-1) / fact_epsilon
[2127]	207	! zw2_est(ig,l+1) = Max(0.0001,exp(-zw2fact)*(zw2_est(ig,l)-zdw2)+zdw2)
	208	! zw2_est(ig,l+1) = Max(0.0001,exp(-zw2fact)*(zw2_est(ig,l)-zdw2bis)+zdw2)
	209	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	210	! AB: own derivation for zw2_est (Rio et al. 2010)
	211	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[2113]	212	! zw2fact = 2. * fact_epsilon * zdz
	213	! zdw2 = 2. * afact * zbuoy(ig,l) * zdz
	214	zw2fact = 2. * fact_epsilon * zdz / (1. + betalpha)
	215	zdw2 = 2. * afact * zbuoy(ig,l) * zdz / (1. + betalpha)
[2127]	216	zw2_est(ig,l+1) = Max(0., exp(-zw2fact) * zw2_est(ig,l) + zdw2)
[2060]	217	ENDIF
	218	ENDDO
	219
[2127]	220	!-------------------------------------------------------------------------------
	221	! Mass flux, entrainment and detrainment
	222	!-------------------------------------------------------------------------------
[2060]	223
	224	DO ig=1,ngrid
	225	IF (active(ig)) THEN
	226
	227	zdz = zlev(ig,l+1) - zlev(ig,l)
[2143]	228	zw2m = (zw2_est(ig,l+1) + zw2_est(ig,l)) / 2.
[2127]	229	gamma = afact * zbuoy(ig,l) - fact_epsilon * zw2m
[2060]	230
[2143]	231	IF (zw2_est(ig,l) > 0.) THEN
[2127]	232	test = gamma / zw2m - nu
[2060]	233	ELSE
[2127]	234	print *, 'ERROR: zw2_est is negative while plume is active!'
	235	print *, 'ig,l', ig, l
	236	print *, 'zw2_est', zw2_est(ig,l)
	237	call abort
[2060]	238	ENDIF
	239
[2143]	240	IF (test > 0.) THEN
	241	detr_star(ig,l) = zdz * f_star(ig,l) * nu
	242	entr_star(ig,l) = zdz * f_star(ig,l) * (zbetalpha * gamma / zw2m + nu)
[2060]	243	ELSE
[2143]	244	detr_star(ig,l) = zdz * f_star(ig,l) * (nu - betalpha * gamma / zw2m)
	245	entr_star(ig,l) = zdz * f_star(ig,l) * nu
[2060]	246	ENDIF
	247
[2127]	248	f_star(ig,l+1) = f_star(ig,l) + entr_star(ig,l) - detr_star(ig,l)
[2060]	249
	250	ENDIF
	251	ENDDO
	252
[2127]	253	!-------------------------------------------------------------------------------
	254	! Mixing between thermal plume and environment
	255	!-------------------------------------------------------------------------------
[2060]	256
[2143]	257	activetmp(:) = active(:).and.(f_star(:,l+1) > 1.e-10)
[2060]	258
	259	DO ig=1,ngrid
	260	IF (activetmp(ig)) THEN
[2127]	261	zhla(ig,l) = (f_star(ig,l) * zhla(ig,l-1) & ! zhla is set to TP in plume (mixed)
	262	& + entr_star(ig,l) * zhl(ig,l)) &
[2065]	263	& / (f_star(ig,l+1) + detr_star(ig,l))
	264	zqta(ig,l) = (f_star(ig,l) * zqta(ig,l-1) + & ! zqta is set to qt in plume (mixed)
[2127]	265	& + entr_star(ig,l) * zqt(ig,l)) &
[2065]	266	& / (f_star(ig,l+1) + detr_star(ig,l))
[2060]	267	ENDIF
	268	ENDDO
	269
[2143]	270	!-------------------------------------------------------------------------------
	271	! Latent heat release (after mixing)
	272	!-------------------------------------------------------------------------------
	273
[2127]	274	ztemp(:) = zpopsk(:,l) * zhla(:,l)
[2060]	275
	276	DO ig=1,ngrid
	277	IF (activetmp(ig)) THEN
[2127]	278	CALL Psat_water(ztemp(ig), pplev(ig,l), psat, zqsa(ig,l))
[2060]	279	ENDIF
	280	ENDDO
	281
[2127]	282	!-------------------------------------------------------------------------------
[2143]	283	! Vertical speed (after mixing)
[2127]	284	!-------------------------------------------------------------------------------
[2060]	285
	286	DO ig=1,ngrid
	287	IF (activetmp(ig)) THEN
[2143]	288	zqla(ig,l) = MAX(0.,zqta(ig,l) - zqsa(ig,l)) ! zqla is set to ql plume (mixed)
	289	zta(ig,l) = zhla(ig,l) * zpopsk(ig,l) & ! ztva is set to TR plume (mixed)
	290	& + RLvCp * zqla(ig,l)
[2127]	291	ztva(ig,l) = zta(ig,l) / zpopsk(ig,l) & ! ztva is set to TRPV plume (mixed)
	292	& * (1. + RETV*(zqta(ig,l)-zqla(ig,l)) - zqla(ig,l))
	293
[2060]	294	zbuoy(ig,l) = RG * (ztva(ig,l) - ztv(ig,l)) / ztv(ig,l)
	295	zdz = zlev(ig,l+1) - zlev(ig,l)
	296
[2127]	297	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	298	! AB: initial formula
	299	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[2060]	300	! zw2fact = fact_epsilon * 2. * zdz / (1. + betalpha)
	301	! zdw2 = afact * zbuoy(ig,l) / fact_epsilon
	302	! zw2(ig,l+1) = Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2)+zdw2)
[2127]	303	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	304	! AB: own derivation for zw2 (Rio et al. 2010)
	305	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[2113]	306	! zw2fact = 2. * (fact_epsilon * zdz + entr_star(ig,l) / f_star(ig,l))
	307	! zdw2 = 2. * afact * zbuoy(ig,l) * zdz
	308	zw2fact = 2. * fact_epsilon * zdz / (1. + betalpha)
	309	zdw2 = 2. * afact * zbuoy(ig,l) * zdz / (1. + betalpha)
[2060]	310	zw2(ig,l+1) = Max(0., exp(-zw2fact) * zw2(ig,l) + zdw2)
	311	ENDIF
	312	ENDDO
	313
	314	ENDDO
	315
	316
[2143]	317	RETURN
[2060]	318	END

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

Download in other formats:

Original Format