Context Navigation

calltherm_interface.F90 @ 506

Last change on this file since 506 was 499, checked in by acolaitis, 14 years ago

M 498 libf/phymars/calltherm_interface.F90
^{------------------> Added computation of wstar. Wmax is no longer an output.}

M 498 libf/phymars/physiq.F
^{------------------> Minor changes.}

M 498 libf/phymars/vdif_cd.F
^{------------------> Changed computation of bulk Richardson number to allow more unstable configurations.}

M 498 libf/phymars/vdifc.F
^{------------------> Changed computation of surface conductance coefficients following a new approach. Comparisons with LES are much better.}

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

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

Download in other formats:

Original Format