Context Navigation

thermcell_plume.F90 @ 2613

Last change on this file since 2613 was 2470, checked in by yjaziri, 4 years ago
Generic GCM: global1d and szangle for 1D simulation moved from callcorrk to callkeys to defined a consistent 1D sza in physiq_mod used also in chemistry + typo YJ
File size: 13.0 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
	28
	29	IMPLICIT NONE
	30
	31
[2127]	32	!===============================================================================
[2060]	33	! Declaration
[2127]	34	!===============================================================================
[2060]	35
[2113]	36	! Inputs:
	37	! -------
[2060]	38
[2178]	39	INTEGER, INTENT(in) :: ngrid
	40	INTEGER, INTENT(in) :: nlay
	41	INTEGER, INTENT(in) :: nq
[2060]	42
[2232]	43	INTEGER, INTENT(in) :: lbot(ngrid) ! First considered layer
	44
[2178]	45	REAL, INTENT(in) :: ptimestep
	46	REAL, INTENT(in) :: zlev(ngrid,nlay+1) ! Levels altitude
	47	REAL, INTENT(in) :: pplev(ngrid,nlay+1) ! Levels pressure
	48	REAL, INTENT(in) :: zpopsk(ngrid,nlay) ! Exner function
[2127]	49
[2178]	50	REAL, INTENT(in) :: ztv(ngrid,nlay) ! TRPV environment
	51	REAL, INTENT(in) :: zhl(ngrid,nlay) ! TP environment
	52	REAL, INTENT(in) :: zqt(ngrid,nlay) ! qt environment
	53	REAL, INTENT(in) :: zql(ngrid,nlay) ! ql environment
[2060]	54
[2113]	55	! Outputs:
	56	! --------
[2060]	57
[2178]	58	INTEGER, INTENT(out) :: lmin(ngrid) ! Plume bottom level (first unstable level)
[2060]	59
[2178]	60	REAL, INTENT(out) :: detr_star(ngrid,nlay) ! Normalized detrainment
	61	REAL, INTENT(out) :: entr_star(ngrid,nlay) ! Normalized entrainment
	62	REAL, INTENT(out) :: f_star(ngrid,nlay+1) ! Normalized mass flux
[2060]	63
[2178]	64	REAL, INTENT(out) :: ztva(ngrid,nlay) ! TRPV plume (after mixing)
	65	REAL, INTENT(out) :: zhla(ngrid,nlay) ! TP plume (after mixing)
	66	REAL, INTENT(out) :: zqla(ngrid,nlay) ! ql plume (after mixing)
	67	REAL, INTENT(out) :: zqta(ngrid,nlay) ! qt plume (after mixing)
	68	REAL, INTENT(out) :: zqsa(ngrid,nlay) ! qsat plume (after mixing)
[2232]	69	REAL, INTENT(out) :: zw2(ngrid,nlay+1) ! w plume (after mixing)
[2060]	70
[2113]	71	! Local:
	72	! ------
[2060]	73
	74	INTEGER ig, l, k
[2232]	75	INTEGER l_start
[2060]	76
[2178]	77	REAL ztva_est(ngrid,nlay) ! TRPV plume (before mixing)
	78	REAL zqla_est(ngrid,nlay) ! ql plume (before mixing)
	79	REAL zta_est(ngrid,nlay) ! TR plume (before mixing)
	80	REAL zqsa_est(ngrid) ! qsat plume (before mixing)
	81	REAL zw2_est(ngrid,nlay+1) ! w plume (before mixing)
[2060]	82
[2178]	83	REAL zta(ngrid,nlay) ! TR plume (after mixing)
[2060]	84
[2178]	85	REAL zbuoy(ngrid,nlay) ! Plume buoyancy
[2232]	86	REAL ztemp(ngrid) ! Temperature to compute saturation vapor pressure
[2178]	87	REAL zdz ! Layers heights
	88	REAL ztv2(ngrid,nlay) ! ztv + d_temp * Dirac(l=linf)
[2127]	89
[2178]	90	REAL zdw2 !
	91	REAL zw2fact !
	92	REAL zw2m ! Average vertical velocity between two successive levels
	93	REAL gamma ! Plume acceleration term (to compute vertical velocity)
	94	REAL test !
[2060]	95
[2178]	96	REAL psat ! Dummy argument for Psat_water()
[2071]	97
[2232]	98	LOGICAL active(ngrid) ! If the plume is active (speed and incoming mass flux > 0)
	99	LOGICAL activetmp(ngrid) ! If the plume is active (active=true and outgoing mass flux > 0)
[2178]	100
[2127]	101	!===============================================================================
[2060]	102	! Initialization
[2127]	103	!===============================================================================
[2060]	104
[2232]	105	ztva(:,:) = ztv(:,:) ! ztva is set to TPV environment
	106	zhla(:,:) = zhl(:,:) ! zhla is set to TP environment
	107	zqta(:,:) = zqt(:,:) ! zqta is set to qt environment
	108	zqla(:,:) = zql(:,:) ! zqla is set to ql environment
	109	zqsa(:,:) = 0.
	110	zw2(:,:) = 0.
[2060]	111
[2232]	112	ztva_est(:,:) = ztv(:,:) ! ztva_est is set to TPV environment
	113	zqla_est(:,:) = zql(:,:) ! zqla_est is set to ql environment
	114	zqsa_est(:) = 0.
	115	zw2_est(:,:) = 0.
[2060]	116
[2178]	117	zbuoy(:,:) = 0.
[2060]	118
[2178]	119	f_star(:,:) = 0.
	120	detr_star(:,:) = 0.
	121	entr_star(:,:) = 0.
[2060]	122
[2232]	123	lmin(:) = lbot(:)
[2060]	124
[2178]	125	ztv2(:,:) = ztv(:,:)
	126	ztv2(:,linf) = ztv(:,linf) + d_temp
[2060]	127
[2143]	128	active(:) = .false.
	129
[2232]	130	l_start = nlay
	131
[2127]	132	!===============================================================================
	133	! First layer computation
	134	!===============================================================================
[2060]	135
	136	DO ig=1,ngrid
[2232]	137	l = lbot(ig)
	138	l_start = MIN(l_start, lbot(ig)+1)
[2143]	139	DO WHILE (.not.active(ig).and.(pplev(ig,l+1) > pres_limit).and.(l < nlay))
	140	zbuoy(ig,l) = RG * (ztv2(ig,l) - ztv2(ig,l+1)) / ztv2(ig,l+1)
[2178]	141	IF (zbuoy(ig,l) > 0.) THEN
[2060]	142	lmin(ig) = l
[2178]	143	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	144	! AB: entrainement and mass flux initial values are set to 1. The physical value
	145	! will be computed thanks to the closure relation in thermcell_closure.
	146	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[2143]	147	entr_star(ig,l) = 1.
	148	f_star(ig,l+1) = 1.
[2178]	149	zdz = zlev(ig,l+1) - zlev(ig,l)
	150	zw2fact = 2. * fact_epsilon * zdz / (1. + betalpha)
	151	zdw2 = 2. * afact * zbuoy(ig,l) * zdz / (1. + betalpha)
	152	zw2_est(ig,l+1) = exp(-zw2fact) * zdw2
[2143]	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
[2232]	164	DO l=l_start,nlay-1
[2060]	165
[2127]	166	!-------------------------------------------------------------------------------
[2143]	167	! Is thermal plume (still) active ?
[2127]	168	!-------------------------------------------------------------------------------
	169
[2060]	170	DO ig=1,ngrid
[2232]	171	active(ig) = (zw2(ig,l) > 1.e-9).and.(f_star(ig,l) > 1.e-9)
[2060]	172	ENDDO
	173
[2143]	174	!-------------------------------------------------------------------------------
	175	! Latent heat release (before mixing)
	176	!-------------------------------------------------------------------------------
	177
[2127]	178	ztemp(:) = zpopsk(:,l) * zhla(:,l-1)
[2065]	179
[2060]	180	DO ig=1,ngrid
[2178]	181	IF (active(ig)) THEN
	182	CALL Psat_water(ztemp(ig), pplev(ig,l), psat, zqsa_est(ig))
	183	ENDIF
[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
[2232]	192	zqla_est(ig,l) = MAX(0.,zqta(ig,l-1) - zqsa_est(ig)) ! zqla_est is set to ql plume
	193	zta_est(ig,l) = zhla(ig,l-1) * zpopsk(ig,l) & ! zta_est is set to TR plume
[2127]	194	& + RLvCp * zqla_est(ig,l)
[2232]	195	ztva_est(ig,l) = zta_est(ig,l) / zpopsk(ig,l) & ! ztva_est is set to TRPV plume
[2127]	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
[2113]	201	zw2fact = 2. * fact_epsilon * zdz / (1. + betalpha)
[2178]	202	zdw2 = afact * zbuoy(ig,l) / fact_epsilon
[2232]	203	zw2_est(ig,l+1) = MAX(0., exp(-zw2fact) * (zw2_est(ig,l) - zdw2) + zdw2)
[2060]	204	ENDIF
	205	ENDDO
	206
[2127]	207	!-------------------------------------------------------------------------------
	208	! Mass flux, entrainment and detrainment
	209	!-------------------------------------------------------------------------------
[2060]	210
	211	DO ig=1,ngrid
	212	IF (active(ig)) THEN
	213
	214	zdz = zlev(ig,l+1) - zlev(ig,l)
[2232]	215	zw2m = (zw2_est(ig,l+1) + zw2(ig,l)) / 2.
[2127]	216	gamma = afact * zbuoy(ig,l) - fact_epsilon * zw2m
[2060]	217
[2232]	218	IF (zw2m > 0.) THEN
[2127]	219	test = gamma / zw2m - nu
[2060]	220	ELSE
[2232]	221	test = 0.
	222	print *, 'WARNING: vertical speed is negative while plume is active!'
[2127]	223	print *, 'ig,l', ig, l
[2232]	224	print *, 'zw2m', zw2m
[2060]	225	ENDIF
	226
[2143]	227	IF (test > 0.) THEN
	228	detr_star(ig,l) = zdz * f_star(ig,l) * nu
[2178]	229	entr_star(ig,l) = zdz * f_star(ig,l) * (betalpha * gamma / zw2m + nu) / (betalpha + 1)
[2060]	230	ELSE
[2178]	231	detr_star(ig,l) = zdz * f_star(ig,l) * ((betalpha + 1) * nu - betalpha * gamma / zw2m)
[2143]	232	entr_star(ig,l) = zdz * f_star(ig,l) * nu
[2060]	233	ENDIF
	234
[2127]	235	f_star(ig,l+1) = f_star(ig,l) + entr_star(ig,l) - detr_star(ig,l)
[2060]	236
	237	ENDIF
	238	ENDDO
	239
[2127]	240	!-------------------------------------------------------------------------------
	241	! Mixing between thermal plume and environment
	242	!-------------------------------------------------------------------------------
[2060]	243
[2232]	244	activetmp(:) = active(:).and.(f_star(:,l+1) > 1.e-9)
[2060]	245
	246	DO ig=1,ngrid
	247	IF (activetmp(ig)) THEN
[2127]	248	zhla(ig,l) = (f_star(ig,l) * zhla(ig,l-1) & ! zhla is set to TP in plume (mixed)
	249	& + entr_star(ig,l) * zhl(ig,l)) &
[2065]	250	& / (f_star(ig,l+1) + detr_star(ig,l))
[2470]	251	zqta(ig,l) = (f_star(ig,l) * zqta(ig,l-1) & ! zqta is set to qt in plume (mixed)
[2127]	252	& + entr_star(ig,l) * zqt(ig,l)) &
[2065]	253	& / (f_star(ig,l+1) + detr_star(ig,l))
[2060]	254	ENDIF
	255	ENDDO
	256
[2143]	257	!-------------------------------------------------------------------------------
	258	! Latent heat release (after mixing)
	259	!-------------------------------------------------------------------------------
	260
[2127]	261	ztemp(:) = zpopsk(:,l) * zhla(:,l)
[2060]	262
	263	DO ig=1,ngrid
	264	IF (activetmp(ig)) THEN
[2127]	265	CALL Psat_water(ztemp(ig), pplev(ig,l), psat, zqsa(ig,l))
[2060]	266	ENDIF
	267	ENDDO
	268
[2127]	269	!-------------------------------------------------------------------------------
[2143]	270	! Vertical speed (after mixing)
[2127]	271	!-------------------------------------------------------------------------------
[2060]	272
	273	DO ig=1,ngrid
	274	IF (activetmp(ig)) THEN
[2143]	275	zqla(ig,l) = MAX(0.,zqta(ig,l) - zqsa(ig,l)) ! zqla is set to ql plume (mixed)
	276	zta(ig,l) = zhla(ig,l) * zpopsk(ig,l) & ! ztva is set to TR plume (mixed)
	277	& + RLvCp * zqla(ig,l)
[2232]	278	ztva(ig,l) = zta(ig,l) / zpopsk(ig,l) & ! ztva is set to TRPV plume (mixed)
[2127]	279	& * (1. + RETV*(zqta(ig,l)-zqla(ig,l)) - zqla(ig,l))
	280
[2060]	281	zbuoy(ig,l) = RG * (ztva(ig,l) - ztv(ig,l)) / ztv(ig,l)
	282	zdz = zlev(ig,l+1) - zlev(ig,l)
	283
[2113]	284	zw2fact = 2. * fact_epsilon * zdz / (1. + betalpha)
[2178]	285	zdw2 = afact * zbuoy(ig,l) / fact_epsilon
[2232]	286	zw2(ig,l+1) = MAX(0., exp(-zw2fact) * (zw2(ig,l) - zdw2) + zdw2)
[2060]	287	ENDIF
	288	ENDDO
	289
	290	ENDDO
	291
	292
[2143]	293	RETURN
[2060]	294	END

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

Download in other formats:

Original Format