Context Navigation

thermcell_plume.F90 @ 3595

Last change on this file since 3595 was 3435, checked in by alesaux, 4 months ago
Generic PCM: Bug fix in the Thermal Plumes Model for generic tracers. Error in computation of potential temperature. ALS
File size: 14.3 KB

Rev	Line
[2060]	1	!
	2	!
	3	!
[2178]	4	SUBROUTINE thermcell_plume(ngrid,nlay,nq,ptimestep, &
	5	ztv,zhl,zqt,zql,zlev,pplev,zpopsk, &
[2127]	6	detr_star,entr_star,f_star, &
[2178]	7	ztva,zhla,zqta,zqla,zqsa, &
[2232]	8	zw2,lbot,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
[3342]	28	USE comcstfi_mod, ONLY: r, cpp, mugaz
	29	USE callkeys_mod, ONLY: water, generic_condensation
	30	USE generic_cloud_common_h, ONLY: Psat_generic, epsi_generic, RLVTT_generic
	31
[2060]	32	IMPLICIT NONE
	33
	34
[2127]	35	!===============================================================================
[2060]	36	! Declaration
[2127]	37	!===============================================================================
[2060]	38
[2113]	39	! Inputs:
	40	! -------
[2060]	41
[2178]	42	INTEGER, INTENT(in) :: ngrid
	43	INTEGER, INTENT(in) :: nlay
	44	INTEGER, INTENT(in) :: nq
[2060]	45
[2232]	46	INTEGER, INTENT(in) :: lbot(ngrid) ! First considered layer
	47
[2178]	48	REAL, INTENT(in) :: ptimestep
	49	REAL, INTENT(in) :: zlev(ngrid,nlay+1) ! Levels altitude
	50	REAL, INTENT(in) :: pplev(ngrid,nlay+1) ! Levels pressure
	51	REAL, INTENT(in) :: zpopsk(ngrid,nlay) ! Exner function
[2127]	52
[2178]	53	REAL, INTENT(in) :: ztv(ngrid,nlay) ! TRPV environment
	54	REAL, INTENT(in) :: zhl(ngrid,nlay) ! TP environment
	55	REAL, INTENT(in) :: zqt(ngrid,nlay) ! qt environment
	56	REAL, INTENT(in) :: zql(ngrid,nlay) ! ql environment
[2060]	57
[2113]	58	! Outputs:
	59	! --------
[2060]	60
[2178]	61	INTEGER, INTENT(out) :: lmin(ngrid) ! Plume bottom level (first unstable level)
[2060]	62
[2178]	63	REAL, INTENT(out) :: detr_star(ngrid,nlay) ! Normalized detrainment
	64	REAL, INTENT(out) :: entr_star(ngrid,nlay) ! Normalized entrainment
	65	REAL, INTENT(out) :: f_star(ngrid,nlay+1) ! Normalized mass flux
[2060]	66
[2178]	67	REAL, INTENT(out) :: ztva(ngrid,nlay) ! TRPV plume (after mixing)
	68	REAL, INTENT(out) :: zhla(ngrid,nlay) ! TP plume (after mixing)
	69	REAL, INTENT(out) :: zqla(ngrid,nlay) ! ql plume (after mixing)
	70	REAL, INTENT(out) :: zqta(ngrid,nlay) ! qt plume (after mixing)
	71	REAL, INTENT(out) :: zqsa(ngrid,nlay) ! qsat plume (after mixing)
[2232]	72	REAL, INTENT(out) :: zw2(ngrid,nlay+1) ! w plume (after mixing)
[2060]	73
[2113]	74	! Local:
	75	! ------
[2060]	76
	77	INTEGER ig, l, k
[2232]	78	INTEGER l_start
[2060]	79
[2178]	80	REAL ztva_est(ngrid,nlay) ! TRPV plume (before mixing)
	81	REAL zqla_est(ngrid,nlay) ! ql plume (before mixing)
	82	REAL zta_est(ngrid,nlay) ! TR plume (before mixing)
	83	REAL zqsa_est(ngrid) ! qsat plume (before mixing)
	84	REAL zw2_est(ngrid,nlay+1) ! w plume (before mixing)
[2060]	85
[2178]	86	REAL zta(ngrid,nlay) ! TR plume (after mixing)
[2060]	87
[2178]	88	REAL zbuoy(ngrid,nlay) ! Plume buoyancy
[2232]	89	REAL ztemp(ngrid) ! Temperature to compute saturation vapor pressure
[2178]	90	REAL zdz ! Layers heights
	91	REAL ztv2(ngrid,nlay) ! ztv + d_temp * Dirac(l=linf)
[2127]	92
[2178]	93	REAL zdw2 !
	94	REAL zw2fact !
	95	REAL zw2m ! Average vertical velocity between two successive levels
	96	REAL gamma ! Plume acceleration term (to compute vertical velocity)
	97	REAL test !
[2060]	98
[2178]	99	REAL psat ! Dummy argument for Psat_water()
[2071]	100
[2232]	101	LOGICAL active(ngrid) ! If the plume is active (speed and incoming mass flux > 0)
	102	LOGICAL activetmp(ngrid) ! If the plume is active (active=true and outgoing mass flux > 0)
[3342]	103
	104	REAL, SAVE :: metallicity ! metallicity of planet --- is not used here, but necessary to call function Psat_generic
	105	REAL :: RV_generic
	106	REAL :: RETV_comp, RLvCp_comp !values used for computation (depends if water or generic tracer)
	107
[2127]	108	!===============================================================================
[2060]	109	! Initialization
[2127]	110	!===============================================================================
[2060]	111
[2232]	112	ztva(:,:) = ztv(:,:) ! ztva is set to TPV environment
	113	zhla(:,:) = zhl(:,:) ! zhla is set to TP environment
	114	zqta(:,:) = zqt(:,:) ! zqta is set to qt environment
	115	zqla(:,:) = zql(:,:) ! zqla is set to ql environment
	116	zqsa(:,:) = 0.
	117	zw2(:,:) = 0.
[2060]	118
[2232]	119	ztva_est(:,:) = ztv(:,:) ! ztva_est is set to TPV environment
	120	zqla_est(:,:) = zql(:,:) ! zqla_est is set to ql environment
	121	zqsa_est(:) = 0.
	122	zw2_est(:,:) = 0.
[2060]	123
[2178]	124	zbuoy(:,:) = 0.
[2060]	125
[2178]	126	f_star(:,:) = 0.
	127	detr_star(:,:) = 0.
	128	entr_star(:,:) = 0.
[2060]	129
[2232]	130	lmin(:) = lbot(:)
[2060]	131
[2178]	132	ztv2(:,:) = ztv(:,:)
	133	ztv2(:,linf) = ztv(:,linf) + d_temp
[2060]	134
[2143]	135	active(:) = .false.
	136
[2232]	137	l_start = nlay
[3342]	138
	139	metallicity=0.0 ! default value --- is not used here but necessary to call function Psat_generic
	140
	141	! ALS24 for thermal plume model with generic tracer
	142	IF (water) THEN
	143	RETV_comp = RETV
	144	RLvCp_comp = RLvCp
	145	ELSEIF (generic_condensation .AND. .NOT. water ) THEN
	146	RV_generic = (8.3145111000.)/(epsi_genericmugaz)
[3435]	147	RETV_comp = RV_generic/r-1.
[3342]	148	RLvCp_comp = RLVTT_generic/cpp
	149	ENDIF
	150
[2127]	151	!===============================================================================
	152	! First layer computation
	153	!===============================================================================
[2060]	154
	155	DO ig=1,ngrid
[2232]	156	l = lbot(ig)
	157	l_start = MIN(l_start, lbot(ig)+1)
[2143]	158	DO WHILE (.not.active(ig).and.(pplev(ig,l+1) > pres_limit).and.(l < nlay))
	159	zbuoy(ig,l) = RG * (ztv2(ig,l) - ztv2(ig,l+1)) / ztv2(ig,l+1)
[2178]	160	IF (zbuoy(ig,l) > 0.) THEN
[2060]	161	lmin(ig) = l
[2178]	162	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	163	! AB: entrainement and mass flux initial values are set to 1. The physical value
	164	! will be computed thanks to the closure relation in thermcell_closure.
	165	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[2143]	166	entr_star(ig,l) = 1.
	167	f_star(ig,l+1) = 1.
[2178]	168	zdz = zlev(ig,l+1) - zlev(ig,l)
	169	zw2fact = 2. * fact_epsilon * zdz / (1. + betalpha)
	170	zdw2 = 2. * afact * zbuoy(ig,l) * zdz / (1. + betalpha)
	171	zw2_est(ig,l+1) = exp(-zw2fact) * zdw2
[2143]	172	zw2(ig,l+1) = zw2_est(ig,l+1)
[2145]	173	active(ig) = .true.
[2060]	174	ENDIF
	175	l = l + 1
	176	ENDDO
	177	ENDDO
	178
[2127]	179	!===============================================================================
	180	! Thermal plumes computations
	181	!===============================================================================
[2060]	182
[2232]	183	DO l=l_start,nlay-1
[2060]	184
[2127]	185	!-------------------------------------------------------------------------------
[2143]	186	! Is thermal plume (still) active ?
[2127]	187	!-------------------------------------------------------------------------------
	188
[2060]	189	DO ig=1,ngrid
[2232]	190	active(ig) = (zw2(ig,l) > 1.e-9).and.(f_star(ig,l) > 1.e-9)
[2060]	191	ENDDO
	192
[2143]	193	!-------------------------------------------------------------------------------
	194	! Latent heat release (before mixing)
	195	!-------------------------------------------------------------------------------
	196
[2127]	197	ztemp(:) = zpopsk(:,l) * zhla(:,l-1)
[2065]	198
[2060]	199	DO ig=1,ngrid
[2178]	200	IF (active(ig)) THEN
[3342]	201	IF (water) THEN
	202	CALL Psat_water(ztemp(ig), pplev(ig,l), psat, zqsa_est(ig))
	203	ELSEIF (generic_condensation .AND. .NOT. water) THEN
	204	CALL Psat_generic(ztemp(ig),pplev(ig,l),metallicity,psat,zqsa_est(ig))
	205	ENDIF
[2178]	206	ENDIF
[2060]	207	ENDDO
	208
[2127]	209	!-------------------------------------------------------------------------------
[2143]	210	! Vertical speed (before mixing)
[2127]	211	!-------------------------------------------------------------------------------
[2060]	212
	213	DO ig=1,ngrid
	214	IF (active(ig)) THEN
[2232]	215	zqla_est(ig,l) = MAX(0.,zqta(ig,l-1) - zqsa_est(ig)) ! zqla_est is set to ql plume
	216	zta_est(ig,l) = zhla(ig,l-1) * zpopsk(ig,l) & ! zta_est is set to TR plume
[3342]	217	& + RLvCp_comp * zqla_est(ig,l)
[2232]	218	ztva_est(ig,l) = zta_est(ig,l) / zpopsk(ig,l) & ! ztva_est is set to TRPV plume
[3342]	219	& * (1. + RETV_comp * (zqta(ig,l-1)-zqla_est(ig,l)) - zqla_est(ig,l))
[2060]	220
[2127]	221	zbuoy(ig,l) = RG * (ztva_est(ig,l) - ztv(ig,l)) / ztv(ig,l)
[2060]	222	zdz = zlev(ig,l+1) - zlev(ig,l)
	223
[2113]	224	zw2fact = 2. * fact_epsilon * zdz / (1. + betalpha)
[2178]	225	zdw2 = afact * zbuoy(ig,l) / fact_epsilon
[2232]	226	zw2_est(ig,l+1) = MAX(0., exp(-zw2fact) * (zw2_est(ig,l) - zdw2) + zdw2)
[2060]	227	ENDIF
	228	ENDDO
	229
[2127]	230	!-------------------------------------------------------------------------------
	231	! Mass flux, entrainment and detrainment
	232	!-------------------------------------------------------------------------------
[2060]	233
	234	DO ig=1,ngrid
	235	IF (active(ig)) THEN
	236
	237	zdz = zlev(ig,l+1) - zlev(ig,l)
[2232]	238	zw2m = (zw2_est(ig,l+1) + zw2(ig,l)) / 2.
[2127]	239	gamma = afact * zbuoy(ig,l) - fact_epsilon * zw2m
[2060]	240
[2232]	241	IF (zw2m > 0.) THEN
[2127]	242	test = gamma / zw2m - nu
[2060]	243	ELSE
[2232]	244	test = 0.
	245	print *, 'WARNING: vertical speed is negative while plume is active!'
[2127]	246	print *, 'ig,l', ig, l
[2232]	247	print *, 'zw2m', zw2m
[2060]	248	ENDIF
	249
[2143]	250	IF (test > 0.) THEN
	251	detr_star(ig,l) = zdz * f_star(ig,l) * nu
[2178]	252	entr_star(ig,l) = zdz * f_star(ig,l) * (betalpha * gamma / zw2m + nu) / (betalpha + 1)
[2060]	253	ELSE
[2178]	254	detr_star(ig,l) = zdz * f_star(ig,l) * ((betalpha + 1) * nu - betalpha * gamma / zw2m)
[2143]	255	entr_star(ig,l) = zdz * f_star(ig,l) * nu
[2060]	256	ENDIF
	257
[2127]	258	f_star(ig,l+1) = f_star(ig,l) + entr_star(ig,l) - detr_star(ig,l)
[2060]	259
	260	ENDIF
	261	ENDDO
	262
[2127]	263	!-------------------------------------------------------------------------------
	264	! Mixing between thermal plume and environment
	265	!-------------------------------------------------------------------------------
[2060]	266
[2232]	267	activetmp(:) = active(:).and.(f_star(:,l+1) > 1.e-9)
[2060]	268
	269	DO ig=1,ngrid
	270	IF (activetmp(ig)) THEN
[2127]	271	zhla(ig,l) = (f_star(ig,l) * zhla(ig,l-1) & ! zhla is set to TP in plume (mixed)
	272	& + entr_star(ig,l) * zhl(ig,l)) &
[2065]	273	& / (f_star(ig,l+1) + detr_star(ig,l))
[2470]	274	zqta(ig,l) = (f_star(ig,l) * zqta(ig,l-1) & ! zqta is set to qt in plume (mixed)
[2127]	275	& + entr_star(ig,l) * zqt(ig,l)) &
[2065]	276	& / (f_star(ig,l+1) + detr_star(ig,l))
[2060]	277	ENDIF
	278	ENDDO
	279
[2143]	280	!-------------------------------------------------------------------------------
	281	! Latent heat release (after mixing)
	282	!-------------------------------------------------------------------------------
	283
[2127]	284	ztemp(:) = zpopsk(:,l) * zhla(:,l)
[2060]	285
	286	DO ig=1,ngrid
	287	IF (activetmp(ig)) THEN
[3342]	288	IF (water) THEN
	289	CALL Psat_water(ztemp(ig), pplev(ig,l), psat, zqsa(ig,l))
	290	ELSEIF (generic_condensation .AND. .NOT. water ) THEN
	291	CALL Psat_generic(ztemp(ig),pplev(ig,l),metallicity,psat,zqsa(ig,l))
	292	ENDIF
[2060]	293	ENDIF
	294	ENDDO
	295
[2127]	296	!-------------------------------------------------------------------------------
[2143]	297	! Vertical speed (after mixing)
[2127]	298	!-------------------------------------------------------------------------------
[2060]	299
	300	DO ig=1,ngrid
	301	IF (activetmp(ig)) THEN
[2143]	302	zqla(ig,l) = MAX(0.,zqta(ig,l) - zqsa(ig,l)) ! zqla is set to ql plume (mixed)
	303	zta(ig,l) = zhla(ig,l) * zpopsk(ig,l) & ! ztva is set to TR plume (mixed)
[3342]	304	& + RLvCp_comp * zqla(ig,l)
[2232]	305	ztva(ig,l) = zta(ig,l) / zpopsk(ig,l) & ! ztva is set to TRPV plume (mixed)
[3342]	306	& * (1. + RETV_comp*(zqta(ig,l)-zqla(ig,l)) - zqla(ig,l))
[2127]	307
[2060]	308	zbuoy(ig,l) = RG * (ztva(ig,l) - ztv(ig,l)) / ztv(ig,l)
	309	zdz = zlev(ig,l+1) - zlev(ig,l)
	310
[2113]	311	zw2fact = 2. * fact_epsilon * zdz / (1. + betalpha)
[2178]	312	zdw2 = afact * zbuoy(ig,l) / fact_epsilon
[2232]	313	zw2(ig,l+1) = MAX(0., exp(-zw2fact) * (zw2(ig,l) - zdw2) + zdw2)
[2060]	314	ENDIF
	315	ENDDO
	316
	317	ENDDO
	318
	319
[2143]	320	RETURN
[2060]	321	END

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

Download in other formats:

Original Format