Context Navigation

calltherm_interface.F90 @ 545

Last change on this file since 545 was 544, checked in by acolaitis, 13 years ago

27/02/12 == AC

Continuation of thermals setting, comparisons with mesoscale results on Case C
Added possibility to call gcm (or 1d) with constant prescribed sensible heat flux, in the spirit of direct comparisons with LES

... This is directly comparable to the variable tke_heat_flux in namelist.input
... Requires the use of yamada4.F from terrestrial GCM (mixes more, seem more numerically stable)
... Usually requires high timesteps (>1000) to avoids crashes. Best approach is to compare
height of first model level z1 and teta_1 between LES and 1D, and increase the timesteps until results
between the two models are comparable (might require a slitghly different tke_heat_flux between the two models
due to difference in vertical diffusion schemes and subgrid effects)
... Syntax for use is to add "tke_heat_flux = ..." in callphys.def

Corrected some stuff with tke transport in thermals (which is anyway deactivated but one day maybe...)

File size: 16.1 KB

Rev	Line
[161]	1	!
	2	! AC 2011-01-05
	3	!
[185]	4	SUBROUTINE calltherm_interface (firstcall, &
[161]	5	& long,lati,zzlev,zzlay, &
	6	& ptimestep,pu,pv,pt,pq,pdu,pdv,pdt,pdq,q2, &
[185]	7	& pplay,pplev,pphi,zpopsk, &
[499]	8	& pdu_th,pdv_th,pdt_th,pdq_th,lmax,zmaxth,pbl_dtke,hfmax,wstar)
[161]	9
[342]	10	USE ioipsl_getincom
[161]	11
	12	implicit none
	13	#include "callkeys.h"
[185]	14	#include "dimensions.h"
	15	#include "dimphys.h"
[342]	16	#include "comcstfi.h"
[508]	17	#include "tracer.h"
[185]	18
[161]	19	!--------------------------------------------------------
[342]	20	! Input Variables
[161]	21	!--------------------------------------------------------
	22
	23	REAL, INTENT(IN) :: ptimestep
[185]	24	REAL, INTENT(IN) :: pplev(ngridmx,nlayermx+1),pplay(ngridmx,nlayermx)
	25	REAL, INTENT(IN) :: pphi(ngridmx,nlayermx)
	26	REAL, INTENT(IN) :: pu(ngridmx,nlayermx),pv(ngridmx,nlayermx)
	27	REAL, INTENT(IN) :: pt(ngridmx,nlayermx),pq(ngridmx,nlayermx,nqmx)
	28	REAL, INTENT(IN) :: zzlay(ngridmx,nlayermx)
	29	REAL, INTENT(IN) :: zzlev(ngridmx,nlayermx+1)
[161]	30	LOGICAL, INTENT(IN) :: firstcall
[185]	31	REAL, INTENT(IN) :: pdu(ngridmx,nlayermx),pdv(ngridmx,nlayermx)
	32	REAL, INTENT(IN) :: pdq(ngridmx,nlayermx,nqmx),pdt(ngridmx,nlayermx)
	33	REAL, INTENT(IN) :: q2(ngridmx,nlayermx+1)
	34	REAL, INTENT(IN) :: long(ngridmx),lati(ngridmx)
	35	REAL, INTENT(IN) :: zpopsk(ngridmx,nlayermx)
[161]	36
	37	!--------------------------------------------------------
[342]	38	! Output Variables
[161]	39	!--------------------------------------------------------
	40
[342]	41	REAL, INTENT(OUT) :: pdu_th(ngridmx,nlayermx)
	42	REAL, INTENT(OUT) :: pdv_th(ngridmx,nlayermx)
	43	REAL, INTENT(OUT) :: pdt_th(ngridmx,nlayermx)
	44	REAL, INTENT(OUT) :: pdq_th(ngridmx,nlayermx,nqmx)
	45	INTEGER, INTENT(OUT) :: lmax(ngridmx)
	46	REAL, INTENT(OUT) :: zmaxth(ngridmx)
	47	REAL, INTENT(OUT) :: pbl_dtke(ngridmx,nlayermx+1)
[499]	48	REAL, INTENT(OUT) :: wstar(ngridmx)
[161]	49
	50	!--------------------------------------------------------
[342]	51	! Thermals local variables
[161]	52	!--------------------------------------------------------
[342]	53	REAL zu(ngridmx,nlayermx), zv(ngridmx,nlayermx)
	54	REAL zt(ngridmx,nlayermx)
[185]	55	REAL d_t_ajs(ngridmx,nlayermx)
	56	REAL d_u_ajs(ngridmx,nlayermx), d_q_ajs(ngridmx,nlayermx,nqmx)
	57	REAL d_v_ajs(ngridmx,nlayermx)
	58	REAL fm_therm(ngridmx,nlayermx+1), entr_therm(ngridmx,nlayermx)
	59	REAL detr_therm(ngridmx,nlayermx)
	60	REAL zw2(ngridmx,nlayermx+1)
[512]	61	REAL fraca(ngridmx,nlayermx+1),zfraca(ngridmx,nlayermx+1)
[185]	62	REAL ztla(ngridmx,nlayermx)
	63	REAL q_therm(ngridmx,nlayermx), pq_therm(ngridmx,nlayermx,nqmx)
	64	REAL q2_therm(ngridmx,nlayermx), dq2_therm(ngridmx,nlayermx)
[342]	65	REAL lmax_real(ngridmx)
	66	REAL masse(ngridmx,nlayermx)
	67	REAL zdz(ngridmx,nlayermx)
[161]	68	LOGICAL qtransport_thermals,dtke_thermals
[499]	69	INTEGER l,ig,iq,ii(1)
[342]	70	CHARACTER (LEN=20) :: modname
[161]	71
[342]	72	!--------------------------------------------------------
	73	! Local variables for sub-timestep
	74	!--------------------------------------------------------
[161]	75
[342]	76	REAL d_t_the(ngridmx,nlayermx), d_q_the(ngridmx,nlayermx,nqmx)
	77	REAL d_u_the(ngridmx,nlayermx),d_v_the(ngridmx,nlayermx)
	78	REAL dq2_the(ngridmx,nlayermx)
	79	INTEGER isplit,nsplit_thermals
	80	REAL r_aspect_thermals
	81	REAL fact
	82	REAL zfm_therm(ngridmx,nlayermx+1),zdt
	83	REAL zentr_therm(ngridmx,nlayermx),zdetr_therm(ngridmx,nlayermx)
	84	REAL zheatFlux(ngridmx,nlayermx)
	85	REAL zheatFlux_down(ngridmx,nlayermx)
	86	REAL zbuoyancyOut(ngridmx,nlayermx)
	87	REAL zbuoyancyEst(ngridmx,nlayermx)
	88	REAL zzw2(ngridmx,nlayermx+1)
	89	REAL zmax(ngridmx)
	90
	91	!--------------------------------------------------------
	92	! Diagnostics
	93	!--------------------------------------------------------
	94
[185]	95	REAL heatFlux(ngridmx,nlayermx)
	96	REAL heatFlux_down(ngridmx,nlayermx)
	97	REAL buoyancyOut(ngridmx,nlayermx)
	98	REAL buoyancyEst(ngridmx,nlayermx)
	99	REAL hfmax(ngridmx),wmax(ngridmx)
[499]	100	REAL pbl_teta(ngridmx),dteta(ngridmx,nlayermx)
[161]	101
[508]	102	!--------------------------------------------------------
	103	! Theta_m
	104	!--------------------------------------------------------
[342]	105
[508]	106	INTEGER ico2
	107	SAVE ico2
[342]	108
[161]	109	! **********************************************************************
[342]	110	! Initialization
[161]	111	! **********************************************************************
	112
[342]	113	lmax(:)=0.
[161]	114	pdu_th(:,:)=0.
	115	pdv_th(:,:)=0.
	116	pdt_th(:,:)=0.
	117	entr_therm(:,:)=0.
	118	detr_therm(:,:)=0.
	119	q2_therm(:,:)=0.
	120	dq2_therm(:,:)=0.
	121	ztla(:,:)=0.
	122	pbl_dtke(:,:)=0.
	123	fm_therm(:,:)=0.
	124	zw2(:,:)=0.
	125	fraca(:,:)=0.
[512]	126	zfraca(:,:)=0.
[161]	127	if (tracer) then
	128	pdq_th(:,:,:)=0.
	129	end if
[342]	130	d_t_ajs(:,:)=0.
	131	d_u_ajs(:,:)=0.
	132	d_v_ajs(:,:)=0.
	133	d_q_ajs(:,:,:)=0.
	134	heatFlux(:,:)=0.
	135	heatFlux_down(:,:)=0.
	136	buoyancyOut(:,:)=0.
	137	buoyancyEst(:,:)=0.
	138	zmaxth(:)=0.
	139	lmax_real(:)=0.
[161]	140
	141
[342]	142	! **********************************************************************
	143	! Preparing inputs for the thermals
	144	! **********************************************************************
[161]	145
[342]	146	zu(:,:)=pu(:,:)+pdu(:,:)*ptimestep
	147	zv(:,:)=pv(:,:)+pdv(:,:)*ptimestep
	148	zt(:,:)=pt(:,:)+pdt(:,:)*ptimestep
[161]	149
[342]	150	pq_therm(:,:,:)=0.
	151	qtransport_thermals=.true. !! default setting
	152	!call getin("qtransport_thermals",qtransport_thermals)
[161]	153
[342]	154	if(qtransport_thermals) then
	155	if(tracer) then
	156	pq_therm(:,:,:)=pq(:,:,:)+pdq(:,:,:)*ptimestep
	157	endif
	158	endif
[161]	159
[544]	160	dtke_thermals=.false. !! default setting
	161	call getin("dtke_thermals",dtke_thermals)
	162	IF(dtke_thermals) THEN
	163	DO l=1,nlayermx
	164	q2_therm(:,l)=0.5*(q2(:,l)+q2(:,l+1))
	165	ENDDO
	166	ENDIF
[342]	167
	168	! **********************************************************************
[508]	169	! Polar night mixing : theta_m
[342]	170	! **********************************************************************
[508]	171
	172	if(firstcall) then
	173	ico2=0
	174	if (tracer) then
	175	! Prepare Special treatment if one of the tracers is CO2 gas
	176	do iq=1,nqmx
	177	if (noms(iq).eq."co2") then
	178	ico2=iq
	179	end if
	180	enddo
	181	endif
	182	endif !of if firstcall
	183
	184
[342]	185	! **********************************************************************
[508]	186	! **********************************************************************
	187	! **********************************************************************
[342]	188	! CALLTHERM
	189	! **********************************************************************
	190	! **********************************************************************
	191	! **********************************************************************
	192
	193	! r_aspect_thermals ! ultimately conrols the amount of mass going through the thermals
	194	! decreasing it increases the thermals effect. Tests at gcm resolution
	195	! shows that too low values destabilize the model
	196	! when changing this value, one should check that the surface layer model
	197	! outputs the correct Cdu and Chu through changing the gustiness coefficient beta
	198
	199	#ifdef MESOSCALE
	200	!! valid for timesteps < 200s
	201	nsplit_thermals=2
	202	r_aspect_thermals=0.7
	203	#else
	204	nsplit_thermals=35
	205	r_aspect_thermals=1.5
	206	#endif
	207	call getin("nsplit_thermals",nsplit_thermals)
	208	call getin("r_aspect_thermals",r_aspect_thermals)
	209
	210	! **********************************************************************
	211	! SUB-TIMESTEP LOOP
	212	! **********************************************************************
	213
	214	zdt=ptimestep/REAL(nsplit_thermals)
	215
	216	DO isplit=1,nsplit_thermals
	217
	218	! Initialization of intermediary variables
	219
	220	zfm_therm(:,:)=0.
	221	zentr_therm(:,:)=0.
	222	zdetr_therm(:,:)=0.
	223	zheatFlux(:,:)=0.
	224	zheatFlux_down(:,:)=0.
[508]	225	zbuoyancyOut(:,:)=0.
	226	zbuoyancyEst(:,:)=0.
[342]	227	zzw2(:,:)=0.
	228	zmax(:)=0.
	229	lmax(:)=0.
	230	d_t_the(:,:)=0.
	231	d_u_the(:,:)=0.
	232	d_v_the(:,:)=0.
	233	dq2_the(:,:)=0.
	234	if (nqmx .ne. 0) then
	235	d_q_the(:,:,:)=0.
[161]	236	endif
	237
[342]	238	CALL thermcell_main_mars(zdt &
	239	& ,pplay,pplev,pphi,zzlev,zzlay &
	240	& ,zu,zv,zt,pq_therm,q2_therm &
	241	& ,d_u_the,d_v_the,d_t_the,d_q_the,dq2_the &
	242	& ,zfm_therm,zentr_therm,zdetr_therm,lmax,zmax &
	243	& ,r_aspect_thermals &
	244	& ,zzw2,fraca,zpopsk &
	245	& ,ztla,zheatFlux,zheatFlux_down &
	246	& ,zbuoyancyOut,zbuoyancyEst)
[161]	247
[342]	248	fact=1./REAL(nsplit_thermals)
[161]	249
[342]	250	d_t_the(:,:)=d_t_the(:,:)ptimestepfact
	251	! d_u_the(:,:)=d_u_the(:,:)*fact
	252	! d_v_the(:,:)=d_v_the(:,:)*fact
	253	! dq2_the(:,:)=dq2_the(:,:)*fact
[508]	254	if (ico2 .ne. 0) then
	255	d_q_the(:,:,ico2)=d_q_the(:,:,ico2)*fact
	256	endif
[161]	257
[342]	258	zmaxth(:)=zmaxth(:)+zmax(:)*fact
	259	lmax_real(:)=lmax_real(:)+float(lmax(:))*fact
	260	fm_therm(:,:)=fm_therm(:,:) &
	261	& +zfm_therm(:,:)*fact
	262	entr_therm(:,:)=entr_therm(:,:) &
	263	& +zentr_therm(:,:)*fact
	264	detr_therm(:,:)=detr_therm(:,:) &
	265	& +zdetr_therm(:,:)*fact
[512]	266	zfraca(:,:)=zfraca(:,:) + fraca(:,:)*fact
[342]	267
	268	heatFlux(:,:)=heatFlux(:,:) &
	269	& +zheatFlux(:,:)*fact
	270	heatFlux_down(:,:)=heatFlux_down(:,:) &
	271	& +zheatFlux_down(:,:)*fact
[508]	272	buoyancyOut(:,:)=buoyancyOut(:,:) &
	273	& +zbuoyancyOut(:,:)*fact
	274	buoyancyEst(:,:)=buoyancyEst(:,:) &
	275	& +zbuoyancyEst(:,:)*fact
[512]	276
[342]	277
	278	zw2(:,:)=zw2(:,:) + zzw2(:,:)*fact
	279
	280	! accumulation de la tendance
	281
	282	d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_the(:,:)
	283	! d_u_ajs(:,:)=d_u_ajs(:,:)+d_u_the(:,:)
	284	! d_v_ajs(:,:)=d_v_ajs(:,:)+d_v_the(:,:)
[508]	285	if (ico2 .ne. 0) then
	286	d_q_ajs(:,:,ico2)=d_q_ajs(:,:,ico2)+d_q_the(:,:,ico2)
	287	endif
[342]	288	! dq2_therm(:,:)=dq2_therm(:,:)+dq2_the(:,:)
	289	! incrementation des variables meteo
	290
	291	zt(:,:) = zt(:,:) + d_t_the(:,:)
	292	! zu(:,:) = zu(:,:) + d_u_the(:,:)
	293	! zv(:,:) = zv(:,:) + d_v_the(:,:)
[508]	294	if (ico2 .ne. 0) then
	295	pq_therm(:,:,ico2) = &
	296	& pq_therm(:,:,ico2) + d_q_the(:,:,ico2)*ptimestep
	297	endif
[342]	298	! q2_therm(:,:) = q2_therm(:,:) + dq2_therm(:,:)
	299
	300
	301	ENDDO ! isplit
	302	!****************************************************************
	303
	304	! Now that we have computed total entrainment and detrainment, we can
	305	! advect u, v, and q in thermals. (theta already advected). We can do
	306	! that separatly because u,v,and q are not used in thermcell_main for
	307	! any thermals-related computation : they are purely passive.
	308
	309	! mass of cells
	310	do l=1,nlayermx
	311	masse(:,l)=(pplev(:,l)-pplev(:,l+1))/g
	312	enddo
	313
	314	! thickness of layers
	315	do l=1,nlayermx
	316	zdz(:,l)=zzlev(:,l+1)-zzlev(:,l)
	317	enddo
	318
	319	modname='momentum'
	320	call thermcell_dqup(ngridmx,nlayermx,ptimestep &
	321	& ,fm_therm,entr_therm,detr_therm, &
	322	& masse,zu,d_u_ajs,modname,zdz)
	323
	324	call thermcell_dqup(ngridmx,nlayermx,ptimestep &
	325	& ,fm_therm,entr_therm,detr_therm, &
	326	& masse,zv,d_v_ajs,modname,zdz)
	327
	328	if (nqmx .ne. 0.) then
	329	modname='tracer'
	330	DO iq=1,nqmx
[508]	331	if (iq .ne. ico2) then
[342]	332	call thermcell_dqup(ngridmx,nlayermx,ptimestep &
	333	& ,fm_therm,entr_therm,detr_therm, &
	334	& masse,pq_therm(:,:,iq),d_q_ajs(:,:,iq),modname,zdz)
[508]	335	endif
[342]	336	ENDDO
	337	endif
	338
[544]	339	if (dtke_thermals) then
	340	modname='tke'
	341	call thermcell_dqup(ngridmx,nlayermx,ptimestep &
	342	& ,fm_therm,entr_therm,detr_therm, &
	343	& masse,q2_therm,dq2_therm,modname,zdz)
	344	endif
	345
[342]	346	DO ig=1,ngridmx
	347	hfmax(ig)=MAXVAL(heatFlux(ig,:)+heatFlux_down(ig,:))
	348	wmax(ig)=MAXVAL(zw2(ig,:))
	349	ENDDO
	350
	351	lmax(:)=nint(lmax_real(:))
	352
	353	! **********************************************************************
	354	! **********************************************************************
	355	! **********************************************************************
	356	! CALLTHERM END
	357	! **********************************************************************
	358	! **********************************************************************
	359	! **********************************************************************
	360
	361
	362	! **********************************************************************
	363	! Preparing outputs
	364	! **********************************************************************
	365
	366	! Winds and tracers PDU, PDV, and PDQ are DERIVATIVES
	367
	368	pdu_th(:,:)=d_u_ajs(:,:)
	369	pdv_th(:,:)=d_v_ajs(:,:)
	370
[161]	371	if(qtransport_thermals) then
[342]	372	if(tracer) then
	373	pdq_th(:,:,:)=d_q_ajs(:,:,:)
	374	endif
[161]	375	endif
	376
[544]	377	IF(dtke_thermals) THEN
	378	DO l=2,nlayermx
	379	pbl_dtke(:,l)=0.5*(dq2_therm(:,l-1)+dq2_therm(:,l))
	380	ENDDO
	381
	382	pbl_dtke(:,1)=0.5*dq2_therm(:,1)
	383	pbl_dtke(:,nlayermx+1)=0.
	384	ENDIF
[161]	385
	386
[342]	387	! Temperature PDT is a TENDANCY
	388	pdt_th(:,:)=d_t_ajs(:,:)/ptimestep
	389
[499]	390
[342]	391	! **********************************************************************
[499]	392	! Compute the free convection velocity scale for vdifc
	393	! **********************************************************************
	394
	395
	396	! Potential temperature gradient
	397
	398	dteta(:,nlayermx)=0.
	399	DO l=1,nlayermx-1
	400	DO ig=1, ngridmx
	401	dteta(ig,l) = ((zt(ig,l+1)-zt(ig,l))/zpopsk(ig,l)) &
	402	& /(zzlay(ig,l+1)-zzlay(ig,l))
	403	ENDDO
	404	ENDDO
	405
	406	! Computation of the pbl mixed layer temperature
	407
	408	DO ig=1, ngridmx
	409	ii=MINLOC(abs(dteta(ig,1:lmax(ig))))
	410	pbl_teta(ig) = zt(ig,ii(1))/zpopsk(ig,ii(1))
	411	ENDDO
	412
	413	! We follow Spiga et. al 2010 (QJRMS)
	414	! ------------
	415
	416	DO ig=1, ngridmx
	417	IF (zmax(ig) .gt. 0.) THEN
	418	wstar(ig)=(gzmaxth(ig)hfmax(ig)/pbl_teta(ig))**(1./3.)
	419	ELSE
	420	wstar(ig)=0.
	421	ENDIF
	422	ENDDO
	423
	424
	425
	426	! **********************************************************************
[342]	427	! Diagnostics
	428	! **********************************************************************
[161]	429
	430	if(outptherm) then
[185]	431	if (ngridmx .eq. 1) then
	432	call WRITEDIAGFI(ngridmx,'entr_therm','entrainement thermique',&
[161]	433	& 'kg/m-2',1,entr_therm)
[185]	434	call WRITEDIAGFI(ngridmx,'detr_therm','detrainement thermique',&
[161]	435	& 'kg/m-2',1,detr_therm)
[185]	436	call WRITEDIAGFI(ngridmx,'fm_therm','flux masse thermique',&
[161]	437	& 'kg/m-2',1,fm_therm)
[185]	438	call WRITEDIAGFI(ngridmx,'zw2','vitesse verticale thermique',&
[161]	439	& 'm/s',1,zw2)
[185]	440	call WRITEDIAGFI(ngridmx,'heatFlux_up','heatFlux_updraft',&
[161]	441	& 'SI',1,heatFlux)
[185]	442	call WRITEDIAGFI(ngridmx,'heatFlux_down','heatFlux_downdraft',&
[161]	443	& 'SI',1,heatFlux_down)
[185]	444	call WRITEDIAGFI(ngridmx,'fraca','fraction coverage',&
[161]	445	& 'percent',1,fraca)
[185]	446	call WRITEDIAGFI(ngridmx,'buoyancyOut','buoyancyOut',&
[161]	447	& 'm.s-2',1,buoyancyOut)
[185]	448	call WRITEDIAGFI(ngridmx,'buoyancyEst','buoyancyEst',&
[161]	449	& 'm.s-2',1,buoyancyEst)
[185]	450	call WRITEDIAGFI(ngridmx,'d_t_th', &
[161]	451	& 'tendance temp TH','K',1,d_t_ajs)
[185]	452	call WRITEDIAGFI(ngridmx,'zmax', &
[342]	453	& 'pbl height','m',0,zmaxth)
[161]	454	else
	455
[185]	456	call WRITEDIAGFI(ngridmx,'entr_therm','entrainement thermique',&
[161]	457	& 'kg/m-2',3,entr_therm)
[185]	458	call WRITEDIAGFI(ngridmx,'detr_therm','detrainement thermique',&
[161]	459	& 'kg/m-2',3,detr_therm)
[185]	460	call WRITEDIAGFI(ngridmx,'fm_therm','flux masse thermique',&
[161]	461	& 'kg/m-2',3,fm_therm)
[185]	462	call WRITEDIAGFI(ngridmx,'zw2','vitesse verticale thermique',&
[161]	463	& 'm/s',3,zw2)
[185]	464	call WRITEDIAGFI(ngridmx,'heatFlux','heatFlux',&
[161]	465	& 'SI',3,heatFlux)
[185]	466	call WRITEDIAGFI(ngridmx,'buoyancyOut','buoyancyOut',&
[161]	467	& 'SI',3,buoyancyOut)
[185]	468	call WRITEDIAGFI(ngridmx,'d_t_th', &
[161]	469	& 'tendance temp TH','K',3,d_t_ajs)
	470
	471	endif
	472	endif
	473
	474	END

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

Download in other formats:

Original Format