Context Navigation

thermcell_main.F90 @ 2176

Last change on this file since 2176 was 2144, checked in by aboissinot, 5 years ago
Thermal plume model parameters are now set in callphys.def and (re)add flag to choose between thermcell_dv2 or thermcell_dq to mix horizontal momentum.
File size: 18.4 KB

Rev	Line
[2060]	1	!
	2	!
	3	!
[2127]	4	SUBROUTINE thermcell_main(ngrid,nlay,nq,ptimestep,firstcall, &
	5	pplay,pplev,pphi,zpopsk, &
	6	pu,pv,pt,pq, &
	7	pduadj,pdvadj,pdtadj,pdqadj, &
	8	f0,fm0,entr0,detr0,zw2,fraca, &
	9	zqta,zqla,ztv,ztva,zhla,zhl,zqsa, &
[2132]	10	lmin,lmix,lmax)
[2060]	11
[2101]	12
[2127]	13	!===============================================================================
[2060]	14	! Auteurs: Frederic Hourdin, Catherine Rio, Anne Mathieu
	15	! Version du 09.02.07
	16	! Calcul du transport vertical dans la couche limite en presence
	17	! de "thermiques" explicitement representes avec processus nuageux
	18	!
	19	! Reecriture a partir d'un listing papier a Habas, le 14/02/00
	20	!
	21	! le thermique est suppose homogene et dissipe par melange avec
	22	! son environnement. la longueur l_mix controle l'efficacite du
	23	! melange
	24	!
	25	! Le calcul du transport des differentes especes se fait en prenant
	26	! en compte:
	27	! 1. un flux de masse montant
	28	! 2. un flux de masse descendant
	29	! 3. un entrainement
	30	! 4. un detrainement
	31	!
	32	! Modif 2013/01/04 (FH hourdin@lmd.jussieu.fr)
	33	! Introduction of an implicit computation of vertical advection in
	34	! the environment of thermal plumes in thermcell_dq
[2127]	35	! impl = 0 : explicit ; impl = 1 : implicit ; impl =-1 : old version
[2060]	36	! controled by iflag_thermals =
	37	! 15, 16 run with impl=-1 : numerical convergence with NPv3
	38	! 17, 18 run with impl=1 : more stable
	39	! 15 and 17 correspond to the activation of the stratocumulus "bidouille"
[2127]	40	!
	41	! Major changes 2018-19 (AB alexandre.boissinot@lmd.jussieu.fr)
	42	! New detr and entre formulae (no longer alimentation)
	43	! lmin can be greater than 1
	44	! Mix every tracer (EN COURS)
	45	! Old version of thermcell_dq is removed
	46	! Alternative version thermcell_dv2 is removed
	47	!
	48	!===============================================================================
[2060]	49
	50	USE thermcell_mod
[2143]	51	USE print_control_mod, ONLY: prt_level
[2060]	52
	53	IMPLICIT NONE
	54
	55
[2127]	56	!===============================================================================
[2060]	57	! Declaration
[2127]	58	!===============================================================================
[2060]	59
[2127]	60	! Inputs:
	61	! -------
[2060]	62
[2127]	63	INTEGER ngrid, nlay, nq
[2060]	64
	65	REAL ptimestep
[2127]	66	REAL pplay(ngrid,nlay) ! Layer pressure
	67	REAL pplev(ngrid,nlay+1) ! Level pressure
	68	REAL pphi(ngrid,nlay) ! Geopotential
[2060]	69
[2127]	70	REAL pu(ngrid,nlay) ! Zonal wind
	71	REAL pv(ngrid,nlay) ! Meridional wind
	72	REAL pt(ngrid,nlay) ! Temperature
	73	REAL pq(ngrid,nlay,nq) ! Tracers mass mixing ratio
[2060]	74
[2101]	75	LOGICAL firstcall
[2060]	76
[2127]	77	! Outputs:
	78	! --------
[2060]	79
	80	REAL pduadj(ngrid,nlay) ! u convective variations
	81	REAL pdvadj(ngrid,nlay) ! v convective variations
[2127]	82	REAL pdtadj(ngrid,nlay) ! t convective variations
	83	REAL pdqadj(ngrid,nlay,nq) ! q convective variations
[2060]	84
[2127]	85	REAL f0(ngrid) ! mass flux norm (after possible time relaxation)
[2101]	86	REAL fm0(ngrid,nlay+1) ! mass flux (after possible time relaxation)
	87	REAL entr0(ngrid,nlay) ! entrainment (after possible time relaxation)
	88	REAL detr0(ngrid,nlay) ! detrainment (after possible time relaxation)
[2060]	89
[2127]	90	! Local:
	91	! ------
[2060]	92
[2127]	93	INTEGER ig, k, l, iq
	94	INTEGER lmax(ngrid) ! Highest layer reached by the plume
	95	INTEGER lmix(ngrid) ! Layer in which plume vertical speed is maximal
	96	INTEGER lmin(ngrid) ! First unstable layer
[2060]	97
[2143]	98	REAL zmix(ngrid) ! Altitude of maximal vertical speed
	99	REAL zmax(ngrid) ! Maximal altitudes where plumes are active
	100	REAL zmin(ngrid) ! Minimal altitudes where plumes are active
	101
[2127]	102	REAL zlay(ngrid,nlay) ! Layers altitudes
	103	REAL zlev(ngrid,nlay+1) ! Levels altitudes
	104	REAL rho(ngrid,nlay) ! Layers densities
	105	REAL rhobarz(ngrid,nlay) ! Levels densities
	106	REAL masse(ngrid,nlay) ! Layers masses
	107	REAL zpopsk(ngrid,nlay) ! Exner function
[2060]	108
[2127]	109	REAL zu(ngrid,nlay) ! u environment
	110	REAL zv(ngrid,nlay) ! v environment
	111	REAL zt(ngrid,nlay) ! TR environment
	112	REAL zqt(ngrid,nlay) ! qt environment
	113	REAL zql(ngrid,nlay) ! ql environment
	114	REAL zhl(ngrid,nlay) ! TP environment
	115	REAL ztv(ngrid,nlay) ! TRPV environment
	116	REAL zqs(ngrid,nlay) ! qsat environment
[2060]	117
[2127]	118	REAL zua(ngrid,nlay) ! u plume
	119	REAL zva(ngrid,nlay) ! v plume
	120	REAL zta(ngrid,nlay) ! TR plume
[2101]	121	REAL zqla(ngrid,nlay) ! qv plume
	122	REAL zqta(ngrid,nlay) ! qt plume
[2127]	123	REAL zhla(ngrid,nlay) ! TP plume
	124	REAL ztva(ngrid,nlay) ! TRPV plume
	125	REAL zqsa(ngrid,nlay) ! qsat plume
[2060]	126
[2127]	127	REAL zqa(ngrid,nlay,nq) ! q plume (ql=0, qv=qt)
[2060]	128
[2127]	129	REAL f_star(ngrid,nlay+1) ! Normalized mass flux
[2143]	130	REAL entr_star(ngrid,nlay) ! Normalized entrainment (E* = e* dz)
	131	REAL detr_star(ngrid,nlay) ! Normalized detrainment (D* = d* dz)
[2060]	132
[2127]	133	REAL fm(ngrid,nlay+1) ! Mass flux
[2143]	134	REAL entr(ngrid,nlay) ! Entrainment (E = e dz)
	135	REAL detr(ngrid,nlay) ! Detrainment (D = d dz)
[2060]	136
[2143]	137	REAL f(ngrid) ! Mass flux norm
[2127]	138	REAL lambda ! Time relaxation coefficent
[2143]	139	REAL fraca(ngrid,nlay+1) ! Updraft fraction
	140	REAL linter(ngrid) ! Level (continuous) of maximal vertical speed
	141	REAL wmax(ngrid) ! Maximal vertical speed
	142	REAL zw2(ngrid,nlay+1) ! Plume vertical speed
[2127]	143	REAL zdthladj(ngrid,nlay) ! Potential temperature variations
	144	REAL dummy(ngrid,nlay) ! Dummy argument for thermcell_dq()
[2060]	145
[2127]	146	!===============================================================================
[2060]	147	! Initialization
[2127]	148	!===============================================================================
[2060]	149
[2101]	150	IF (firstcall) THEN
[2060]	151	fm0(:,:) = 0.
	152	entr0(:,:) = 0.
	153	detr0(:,:) = 0.
	154	ENDIF
	155
[2127]	156	f_star(:,:) = 0.
	157	entr_star(:,:) = 0.
	158	detr_star(:,:) = 0.
	159
	160	f(:) = 0.
	161
[2060]	162	fm(:,:) = 0.
	163	entr(:,:) = 0.
	164	detr(:,:) = 0.
	165
[2127]	166	lmax(:) = 1
	167	lmix(:) = 1
	168	lmin(:) = 1
	169
	170	pduadj(:,:) = 0.0
	171	pdvadj(:,:) = 0.0
	172	pdtadj(:,:) = 0.0
	173	pdqadj(:,:,:) = 0.0
	174
[2101]	175	DO ig=1,ngrid
[2143]	176	! AB: Careful: Hard-coded value from Earth version!
	177	! f0(ig) = max(f0(ig), 1.e-2)
	178	! AB: No pescribed minimal value for f0
	179	f0(ig) = max(f0(ig), 0.)
[2060]	180	ENDDO
	181
[2127]	182	!===============================================================================
	183	! Environment settings
	184	!===============================================================================
[2060]	185
[2127]	186	!-------------------------------------------------------------------------------
[2060]	187	! Calcul de T,q,ql a partir de Tl et qt dans l environnement
[2127]	188	!-------------------------------------------------------------------------------
[2060]	189
[2127]	190	CALL thermcell_env(ngrid,nlay,nq,pq,pt,pu,pv,pplay,pplev, &
	191	& zqt,zql,zt,ztv,zhl,zu,zv,zpopsk,zqs)
[2060]	192
[2127]	193	!-------------------------------------------------------------------------------
	194	! Levels and layers altitudes
	195	!-------------------------------------------------------------------------------
[2060]	196
	197	DO l=2,nlay
	198	zlev(:,l) = 0.5 * (pphi(:,l) + pphi(:,l-1)) / RG
	199	ENDDO
	200
	201	zlev(:,1) = 0.
[2101]	202	zlev(:,nlay+1) = (2. * pphi(:,nlay) - pphi(:,nlay-1)) / RG
[2060]	203
	204	DO l=1,nlay
[2143]	205	zlay(:,l) = pphi(:,l) / RG
[2060]	206	ENDDO
	207
[2127]	208	!-------------------------------------------------------------------------------
	209	! Levels and layers densities
	210	!-------------------------------------------------------------------------------
[2060]	211
[2127]	212	rho(:,:) = pplay(:,:) / (zpopsk(:,:) * RD * ztv(:,:))
[2060]	213
	214	IF (prt_level.ge.10) THEN
[2143]	215	print *, 'WARNING: density in the first layer is equal to density at the first level!'
	216	print *, 'rhobarz(:,1)', rhobarz(:,1)
	217	print *, 'rho(:,1) ', rho(:,1)
[2060]	218	ENDIF
	219
	220	rhobarz(:,1) = rho(:,1)
	221
	222	DO l=2,nlay
	223	rhobarz(:,l) = 0.5 * (rho(:,l) + rho(:,l-1))
	224	ENDDO
	225
[2127]	226	!-------------------------------------------------------------------------------
	227	! Layers masses
	228	!-------------------------------------------------------------------------------
[2060]	229
	230	DO l=1,nlay
	231	masse(:,l) = (pplev(:,l) - pplev(:,l+1)) / RG
	232	ENDDO
	233
[2127]	234	!===============================================================================
	235	! Explicative schemes
	236	!===============================================================================
	237
	238	!-------------------------------------------------------------------------------
	239	! Thermal plume variables
	240	!-------------------------------------------------------------------------------
[2060]	241
[2127]	242	! top of the model
	243	! ===========================
	244	!
	245	! ---------------------------
	246	! _
	247	! ----- F_lmax+1=0 ------zmax \
	248	! lmax \|
	249	! ------F_lmax>0------------- \|
	250	! \|
	251	! --------------------------- \|
	252	! \|
	253	! --------------------------- \|
	254	! \|
	255	! ------------------wmax,zmix \|
	256	! lmix \|
	257	! --------------------------- \|
	258	! \|
	259	! --------------------------- \|
	260	! \| E, D
	261	! --------------------------- \|
	262	! \|
	263	! --------------------------- rhobarz, f_star, fm, fm0, zw2, fraca
	264	! zt, zu, zv, zo, rho \|
	265	! --------------------------- \|
	266	! \|
	267	! --------------------------- \|
	268	! \|
	269	! --------------------------- \|
	270	! \|
	271	! ------F_lmin+1>0----------- \|
	272	! lmin \|
	273	! ----- F_lmin=0 ------------ _/
	274	!
	275	! ---------------------------
	276	!
	277	! ===========================
	278	! bottom of the model
[2060]	279
[2127]	280	!-------------------------------------------------------------------------------
	281	! Zoom on layers k and k-1
	282	!-------------------------------------------------------------------------------
[2060]	283
[2127]	284	! \| /\|\ \| \|
	285	! \|---- \| F_k+1 -----------\|--------------------------\| level k+1
	286	! \| \| w_k+1 \| \|
	287	! \| --\|--> D_k \|
	288	! \| \| \| layer k
	289	! \| <--\|-- E_k \|
	290	! \| /\|\ \| \|
	291	! \|---- \| F_k ----------\|-----------------------------\| level k
	292	! \| \| w_k \| \|
	293	! \| --\|--> D_k-1 \|
	294	! \| \| \| layer k-1
	295	! \| <--\|-- E_k-1 \|
	296	! \| /\|\ \| \|
	297	! \|---- \| F_k-1 -----\|--------------------------------\| level k-1
	298	! \| w_k-1
	299	! 0 fraca 1
	300	! \__________________/ \______________________________/
	301	! plume (fraca) environment (1-fraca)
[2060]	302
[2127]	303	!===============================================================================
	304	! Thermal plumes computation
	305	!===============================================================================
[2060]	306
[2127]	307	!-------------------------------------------------------------------------------
	308	! Thermal plumes speeds, fluxes, tracers and temperatures
	309	!-------------------------------------------------------------------------------
[2060]	310
[2127]	311	CALL thermcell_plume(ngrid,nlay,nq,ptimestep,ztv, &
	312	& zhl,zqt,zql,rhobarz,zlev,pplev,pphi,zpopsk, &
	313	& detr_star,entr_star,f_star, &
[2143]	314	& ztva,zhla,zqla,zqta,zta,zqsa, &
	315	& zw2,lmix,lmin)
	316
[2127]	317	!-------------------------------------------------------------------------------
[2101]	318	! Thermal plumes characteristics: zmax, zmix, wmax
[2127]	319	!-------------------------------------------------------------------------------
[2060]	320
[2127]	321	! AB: Careful, zw2 became its square root in thermcell_height!
[2143]	322	CALL thermcell_height(ngrid,nlay,lmin,linter,lmix,lmax,zw2, &
	323	& zlev,zmin,zmix,zmax,wmax,f_star)
[2060]	324
[2127]	325	!===============================================================================
[2101]	326	! Closure and mass fluxes computation
[2127]	327	!===============================================================================
[2060]	328
[2127]	329	!-------------------------------------------------------------------------------
	330	! Closure
	331	!-------------------------------------------------------------------------------
[2060]	332
[2127]	333	CALL thermcell_closure(ngrid,nlay,ptimestep,rho,zlev, &
[2143]	334	& lmax,entr_star,zmin,zmax,wmax,f)
[2060]	335
	336	IF (tau_thermals>1.) THEN
	337	lambda = exp(-ptimestep/tau_thermals)
[2102]	338	f0(:) = (1.-lambda) * f(:) + lambda * f0(:)
[2060]	339	ELSE
[2101]	340	f0(:) = f(:)
[2060]	341	ENDIF
	342
[2127]	343	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[2060]	344	! Test valable seulement en 1D mais pas genant
[2127]	345	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[2060]	346	IF (.not. (f0(1).ge.0.) ) THEN
[2143]	347	print *, 'ERROR: mass flux norm is not positive!'
[2101]	348	print *, 'f0 =', f0(1)
[2143]	349	CALL abort
[2060]	350	ENDIF
	351
[2127]	352	!-------------------------------------------------------------------------------
[2101]	353	! Mass fluxes
[2127]	354	!-------------------------------------------------------------------------------
[2060]	355
	356	CALL thermcell_flux(ngrid,nlay,ptimestep,masse, &
[2132]	357	& lmin,lmax,entr_star,detr_star, &
[2127]	358	& f,rhobarz,zlev,zw2,fm,entr,detr,zqla)
[2060]	359
[2127]	360	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[2060]	361	! On ne prend pas directement les profils issus des calculs precedents mais on
	362	! s'autorise genereusement une relaxation vers ceci avec une constante de temps
[2127]	363	! tau_thermals (typiquement 1800s sur Terre).
	364	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[2060]	365
	366	IF (tau_thermals>1.) THEN
	367	lambda = exp(-ptimestep/tau_thermals)
[2064]	368	fm0 = (1.-lambda) * fm + lambda * fm0
	369	entr0 = (1.-lambda) * entr + lambda * entr0
	370	detr0 = (1.-lambda) * detr + lambda * detr0
[2060]	371	ELSE
[2101]	372	fm0(:,:) = fm(:,:)
	373	entr0(:,:) = entr(:,:)
	374	detr0(:,:) = detr(:,:)
[2060]	375	ENDIF
	376
[2127]	377	!-------------------------------------------------------------------------------
[2101]	378	! Updraft fraction
[2127]	379	!-------------------------------------------------------------------------------
[2101]	380
	381	DO ig=1,ngrid
	382	fraca(ig,1) = 0.
	383	fraca(ig,nlay+1) = 0.
	384	ENDDO
	385
	386	DO l=2,nlay
	387	DO ig=1,ngrid
[2143]	388	IF (zw2(ig,l) > 0.) THEN
[2101]	389	fraca(ig,l) = fm(ig,l) / (rhobarz(ig,l) * zw2(ig,l))
	390	ELSE
	391	fraca(ig,l) = 0.
	392	ENDIF
	393	ENDDO
	394	ENDDO
	395
[2127]	396	!===============================================================================
[2101]	397	! Transport vertical
[2127]	398	!===============================================================================
[2060]	399
[2127]	400	!-------------------------------------------------------------------------------
	401	! Calcul du transport vertical de la temperature potentielle
	402	!-------------------------------------------------------------------------------
[2101]	403
[2127]	404	CALL thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,detr0,masse, &
	405	& zhl,zdthladj,dummy,lmin)
[2060]	406
	407	DO l=1,nlay
	408	DO ig=1,ngrid
[2127]	409	pdtadj(ig,l) = zdthladj(ig,l) * zpopsk(ig,l)
[2060]	410	ENDDO
	411	ENDDO
	412
[2127]	413	!-------------------------------------------------------------------------------
	414	! Calcul du transport vertical des traceurs
	415	!-------------------------------------------------------------------------------
[2060]	416
[2127]	417	DO iq=1,nq
	418	CALL thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,detr0,masse, &
	419	& pq(:,:,iq),pdqadj(:,:,iq),zqa(:,:,iq),lmin)
[2102]	420	ENDDO
[2060]	421
[2127]	422	!-------------------------------------------------------------------------------
	423	! Calcul du transport vertical du moment horizontal
	424	!-------------------------------------------------------------------------------
[2060]	425
[2144]	426	IF (dvimpl) THEN
	427	CALL thermcell_dv2(ngrid,nlay,ptimestep,fm0,entr0,detr0,masse,fraca, &
	428	& zmax,zmin,pu,pv,pduadj,pdvadj,zua,zva)
	429	ELSE
	430	CALL thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,detr0,masse, &
	431	& zu,pduadj,zua,lmin)
	432	CALL thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,detr0,masse, &
	433	& zv,pdvadj,zva,lmin)
	434	ENDIF
[2060]	435
	436
	437	RETURN
	438	END

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

Download in other formats:

Original Format