Context Navigation

calltherm_interface.F90 @ 524

Last change on this file since 524 was 512, checked in by acolaitis, 13 years ago
Updated thermals version, that produces quite good overshoots in the inversion region. Series of zmax matches very well the LES and underline one of the weaknesses of convadj. This version will change, as this inversion is only produced in 222 levels. Some resolution tricks must be found for it to work at 32 levels.
File size: 15.9 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
[342]	160	! dtke_thermals=.false. !! default setting
	161	! !call getin("dtke_thermals",dtke_thermals)
	162	!
	163	! IF(dtke_thermals) THEN
	164	! DO l=1,nlayermx
	165	! q2_therm(:,l)=0.5*(q2(:,l)+q2(:,l+1))
	166	! ENDDO
	167	! ENDIF
	168
	169	! **********************************************************************
[508]	170	! Polar night mixing : theta_m
[342]	171	! **********************************************************************
[508]	172
	173	if(firstcall) then
	174	ico2=0
	175	if (tracer) then
	176	! Prepare Special treatment if one of the tracers is CO2 gas
	177	do iq=1,nqmx
	178	if (noms(iq).eq."co2") then
	179	ico2=iq
	180	end if
	181	enddo
	182	endif
	183	endif !of if firstcall
	184
	185
[342]	186	! **********************************************************************
[508]	187	! **********************************************************************
	188	! **********************************************************************
[342]	189	! CALLTHERM
	190	! **********************************************************************
	191	! **********************************************************************
	192	! **********************************************************************
	193
	194	! r_aspect_thermals ! ultimately conrols the amount of mass going through the thermals
	195	! decreasing it increases the thermals effect. Tests at gcm resolution
	196	! shows that too low values destabilize the model
	197	! when changing this value, one should check that the surface layer model
	198	! outputs the correct Cdu and Chu through changing the gustiness coefficient beta
	199
	200	#ifdef MESOSCALE
	201	!! valid for timesteps < 200s
	202	nsplit_thermals=2
	203	r_aspect_thermals=0.7
	204	#else
	205	nsplit_thermals=35
	206	r_aspect_thermals=1.5
	207	#endif
	208	call getin("nsplit_thermals",nsplit_thermals)
	209	call getin("r_aspect_thermals",r_aspect_thermals)
	210
	211	! **********************************************************************
	212	! SUB-TIMESTEP LOOP
	213	! **********************************************************************
	214
	215	zdt=ptimestep/REAL(nsplit_thermals)
	216
	217	DO isplit=1,nsplit_thermals
	218
	219	! Initialization of intermediary variables
	220
	221	zfm_therm(:,:)=0.
	222	zentr_therm(:,:)=0.
	223	zdetr_therm(:,:)=0.
	224	zheatFlux(:,:)=0.
	225	zheatFlux_down(:,:)=0.
[508]	226	zbuoyancyOut(:,:)=0.
	227	zbuoyancyEst(:,:)=0.
[342]	228	zzw2(:,:)=0.
	229	zmax(:)=0.
	230	lmax(:)=0.
	231	d_t_the(:,:)=0.
	232	d_u_the(:,:)=0.
	233	d_v_the(:,:)=0.
	234	dq2_the(:,:)=0.
	235	if (nqmx .ne. 0) then
	236	d_q_the(:,:,:)=0.
[161]	237	endif
	238
[342]	239	CALL thermcell_main_mars(zdt &
	240	& ,pplay,pplev,pphi,zzlev,zzlay &
	241	& ,zu,zv,zt,pq_therm,q2_therm &
	242	& ,d_u_the,d_v_the,d_t_the,d_q_the,dq2_the &
	243	& ,zfm_therm,zentr_therm,zdetr_therm,lmax,zmax &
	244	& ,r_aspect_thermals &
	245	& ,zzw2,fraca,zpopsk &
	246	& ,ztla,zheatFlux,zheatFlux_down &
	247	& ,zbuoyancyOut,zbuoyancyEst)
[161]	248
[342]	249	fact=1./REAL(nsplit_thermals)
[161]	250
[342]	251	d_t_the(:,:)=d_t_the(:,:)ptimestepfact
	252	! d_u_the(:,:)=d_u_the(:,:)*fact
	253	! d_v_the(:,:)=d_v_the(:,:)*fact
	254	! dq2_the(:,:)=dq2_the(:,:)*fact
[508]	255	if (ico2 .ne. 0) then
	256	d_q_the(:,:,ico2)=d_q_the(:,:,ico2)*fact
	257	endif
[161]	258
[342]	259	zmaxth(:)=zmaxth(:)+zmax(:)*fact
	260	lmax_real(:)=lmax_real(:)+float(lmax(:))*fact
	261	fm_therm(:,:)=fm_therm(:,:) &
	262	& +zfm_therm(:,:)*fact
	263	entr_therm(:,:)=entr_therm(:,:) &
	264	& +zentr_therm(:,:)*fact
	265	detr_therm(:,:)=detr_therm(:,:) &
	266	& +zdetr_therm(:,:)*fact
[512]	267	zfraca(:,:)=zfraca(:,:) + fraca(:,:)*fact
[342]	268
	269	heatFlux(:,:)=heatFlux(:,:) &
	270	& +zheatFlux(:,:)*fact
	271	heatFlux_down(:,:)=heatFlux_down(:,:) &
	272	& +zheatFlux_down(:,:)*fact
[508]	273	buoyancyOut(:,:)=buoyancyOut(:,:) &
	274	& +zbuoyancyOut(:,:)*fact
	275	buoyancyEst(:,:)=buoyancyEst(:,:) &
	276	& +zbuoyancyEst(:,:)*fact
[512]	277
[342]	278
	279	zw2(:,:)=zw2(:,:) + zzw2(:,:)*fact
	280
	281	! accumulation de la tendance
	282
	283	d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_the(:,:)
	284	! d_u_ajs(:,:)=d_u_ajs(:,:)+d_u_the(:,:)
	285	! d_v_ajs(:,:)=d_v_ajs(:,:)+d_v_the(:,:)
[508]	286	if (ico2 .ne. 0) then
	287	d_q_ajs(:,:,ico2)=d_q_ajs(:,:,ico2)+d_q_the(:,:,ico2)
	288	endif
[342]	289	! dq2_therm(:,:)=dq2_therm(:,:)+dq2_the(:,:)
	290	! incrementation des variables meteo
	291
	292	zt(:,:) = zt(:,:) + d_t_the(:,:)
	293	! zu(:,:) = zu(:,:) + d_u_the(:,:)
	294	! zv(:,:) = zv(:,:) + d_v_the(:,:)
[508]	295	if (ico2 .ne. 0) then
	296	pq_therm(:,:,ico2) = &
	297	& pq_therm(:,:,ico2) + d_q_the(:,:,ico2)*ptimestep
	298	endif
[342]	299	! q2_therm(:,:) = q2_therm(:,:) + dq2_therm(:,:)
	300
	301
	302	ENDDO ! isplit
	303	!****************************************************************
	304
	305	! Now that we have computed total entrainment and detrainment, we can
	306	! advect u, v, and q in thermals. (theta already advected). We can do
	307	! that separatly because u,v,and q are not used in thermcell_main for
	308	! any thermals-related computation : they are purely passive.
	309
	310	! mass of cells
	311	do l=1,nlayermx
	312	masse(:,l)=(pplev(:,l)-pplev(:,l+1))/g
	313	enddo
	314
	315	! thickness of layers
	316	do l=1,nlayermx
	317	zdz(:,l)=zzlev(:,l+1)-zzlev(:,l)
	318	enddo
	319
	320	modname='momentum'
	321	call thermcell_dqup(ngridmx,nlayermx,ptimestep &
	322	& ,fm_therm,entr_therm,detr_therm, &
	323	& masse,zu,d_u_ajs,modname,zdz)
	324
	325	call thermcell_dqup(ngridmx,nlayermx,ptimestep &
	326	& ,fm_therm,entr_therm,detr_therm, &
	327	& masse,zv,d_v_ajs,modname,zdz)
	328
	329	if (nqmx .ne. 0.) then
	330	modname='tracer'
	331	DO iq=1,nqmx
[508]	332	if (iq .ne. ico2) then
[342]	333	call thermcell_dqup(ngridmx,nlayermx,ptimestep &
	334	& ,fm_therm,entr_therm,detr_therm, &
	335	& masse,pq_therm(:,:,iq),d_q_ajs(:,:,iq),modname,zdz)
[508]	336	endif
[342]	337	ENDDO
	338	endif
	339
	340	DO ig=1,ngridmx
	341	hfmax(ig)=MAXVAL(heatFlux(ig,:)+heatFlux_down(ig,:))
	342	wmax(ig)=MAXVAL(zw2(ig,:))
	343	ENDDO
	344
	345	lmax(:)=nint(lmax_real(:))
	346
	347	! **********************************************************************
	348	! **********************************************************************
	349	! **********************************************************************
	350	! CALLTHERM END
	351	! **********************************************************************
	352	! **********************************************************************
	353	! **********************************************************************
	354
	355
	356	! **********************************************************************
	357	! Preparing outputs
	358	! **********************************************************************
	359
	360	! Winds and tracers PDU, PDV, and PDQ are DERIVATIVES
	361
	362	pdu_th(:,:)=d_u_ajs(:,:)
	363	pdv_th(:,:)=d_v_ajs(:,:)
	364
[161]	365	if(qtransport_thermals) then
[342]	366	if(tracer) then
	367	pdq_th(:,:,:)=d_q_ajs(:,:,:)
	368	endif
[161]	369	endif
	370
[342]	371	! IF(dtke_thermals) THEN
	372	! DO l=2,nlayermx
	373	! pbl_dtke(:,l)=0.5*(dq2_therm(:,l-1)+dq2_therm(:,l))
	374	! ENDDO
	375	!
	376	! pbl_dtke(:,1)=0.5*dq2_therm(:,1)
	377	! pbl_dtke(:,nlayermx+1)=0.
	378	! ENDIF
[161]	379
	380
[342]	381	! Temperature PDT is a TENDANCY
	382	pdt_th(:,:)=d_t_ajs(:,:)/ptimestep
	383
[499]	384
[342]	385	! **********************************************************************
[499]	386	! Compute the free convection velocity scale for vdifc
	387	! **********************************************************************
	388
	389
	390	! Potential temperature gradient
	391
	392	dteta(:,nlayermx)=0.
	393	DO l=1,nlayermx-1
	394	DO ig=1, ngridmx
	395	dteta(ig,l) = ((zt(ig,l+1)-zt(ig,l))/zpopsk(ig,l)) &
	396	& /(zzlay(ig,l+1)-zzlay(ig,l))
	397	ENDDO
	398	ENDDO
	399
	400	! Computation of the pbl mixed layer temperature
	401
	402	DO ig=1, ngridmx
	403	ii=MINLOC(abs(dteta(ig,1:lmax(ig))))
	404	pbl_teta(ig) = zt(ig,ii(1))/zpopsk(ig,ii(1))
	405	ENDDO
	406
	407	! We follow Spiga et. al 2010 (QJRMS)
	408	! ------------
	409
	410	DO ig=1, ngridmx
	411	IF (zmax(ig) .gt. 0.) THEN
	412	wstar(ig)=(gzmaxth(ig)hfmax(ig)/pbl_teta(ig))**(1./3.)
	413	ELSE
	414	wstar(ig)=0.
	415	ENDIF
	416	ENDDO
	417
	418
	419
	420	! **********************************************************************
[342]	421	! Diagnostics
	422	! **********************************************************************
[161]	423
	424	if(outptherm) then
[185]	425	if (ngridmx .eq. 1) then
	426	call WRITEDIAGFI(ngridmx,'entr_therm','entrainement thermique',&
[161]	427	& 'kg/m-2',1,entr_therm)
[185]	428	call WRITEDIAGFI(ngridmx,'detr_therm','detrainement thermique',&
[161]	429	& 'kg/m-2',1,detr_therm)
[185]	430	call WRITEDIAGFI(ngridmx,'fm_therm','flux masse thermique',&
[161]	431	& 'kg/m-2',1,fm_therm)
[185]	432	call WRITEDIAGFI(ngridmx,'zw2','vitesse verticale thermique',&
[161]	433	& 'm/s',1,zw2)
[185]	434	call WRITEDIAGFI(ngridmx,'heatFlux_up','heatFlux_updraft',&
[161]	435	& 'SI',1,heatFlux)
[185]	436	call WRITEDIAGFI(ngridmx,'heatFlux_down','heatFlux_downdraft',&
[161]	437	& 'SI',1,heatFlux_down)
[185]	438	call WRITEDIAGFI(ngridmx,'fraca','fraction coverage',&
[161]	439	& 'percent',1,fraca)
[185]	440	call WRITEDIAGFI(ngridmx,'buoyancyOut','buoyancyOut',&
[161]	441	& 'm.s-2',1,buoyancyOut)
[185]	442	call WRITEDIAGFI(ngridmx,'buoyancyEst','buoyancyEst',&
[161]	443	& 'm.s-2',1,buoyancyEst)
[185]	444	call WRITEDIAGFI(ngridmx,'d_t_th', &
[161]	445	& 'tendance temp TH','K',1,d_t_ajs)
[185]	446	call WRITEDIAGFI(ngridmx,'zmax', &
[342]	447	& 'pbl height','m',0,zmaxth)
[161]	448	else
	449
[185]	450	call WRITEDIAGFI(ngridmx,'entr_therm','entrainement thermique',&
[161]	451	& 'kg/m-2',3,entr_therm)
[185]	452	call WRITEDIAGFI(ngridmx,'detr_therm','detrainement thermique',&
[161]	453	& 'kg/m-2',3,detr_therm)
[185]	454	call WRITEDIAGFI(ngridmx,'fm_therm','flux masse thermique',&
[161]	455	& 'kg/m-2',3,fm_therm)
[185]	456	call WRITEDIAGFI(ngridmx,'zw2','vitesse verticale thermique',&
[161]	457	& 'm/s',3,zw2)
[185]	458	call WRITEDIAGFI(ngridmx,'heatFlux','heatFlux',&
[161]	459	& 'SI',3,heatFlux)
[185]	460	call WRITEDIAGFI(ngridmx,'buoyancyOut','buoyancyOut',&
[161]	461	& 'SI',3,buoyancyOut)
[185]	462	call WRITEDIAGFI(ngridmx,'d_t_th', &
[161]	463	& 'tendance temp TH','K',3,d_t_ajs)
	464
	465	endif
	466	endif
	467
	468	END

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

Download in other formats:

Original Format