Context Navigation

calltherm_interface.F90 @ 577

Last change on this file since 577 was 561, checked in by acolaitis, 13 years ago

Thermals. Added a second instance of subtimestep variation depending on global timestep. Without active clouds without active dust, going from 48 to 96 timestep per days can be compensated in the thermals by cutting the subtimestep from 35 to 10. Resulting thermal structure is better, the model is also observed to be more stable. In this configuration, the cost of doubling the timestep is of 10% only.

File size: 16.5 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)
[561]	79	INTEGER isplit
	80	INTEGER,SAVE :: nsplit_thermals
	81	REAL, SAVE :: r_aspect_thermals
[342]	82	REAL fact
	83	REAL zfm_therm(ngridmx,nlayermx+1),zdt
	84	REAL zentr_therm(ngridmx,nlayermx),zdetr_therm(ngridmx,nlayermx)
	85	REAL zheatFlux(ngridmx,nlayermx)
	86	REAL zheatFlux_down(ngridmx,nlayermx)
	87	REAL zbuoyancyOut(ngridmx,nlayermx)
	88	REAL zbuoyancyEst(ngridmx,nlayermx)
	89	REAL zzw2(ngridmx,nlayermx+1)
	90	REAL zmax(ngridmx)
	91
	92	!--------------------------------------------------------
	93	! Diagnostics
	94	!--------------------------------------------------------
	95
[185]	96	REAL heatFlux(ngridmx,nlayermx)
	97	REAL heatFlux_down(ngridmx,nlayermx)
	98	REAL buoyancyOut(ngridmx,nlayermx)
	99	REAL buoyancyEst(ngridmx,nlayermx)
	100	REAL hfmax(ngridmx),wmax(ngridmx)
[499]	101	REAL pbl_teta(ngridmx),dteta(ngridmx,nlayermx)
[161]	102
[508]	103	!--------------------------------------------------------
	104	! Theta_m
	105	!--------------------------------------------------------
[342]	106
[508]	107	INTEGER ico2
	108	SAVE ico2
[342]	109
[161]	110	! **********************************************************************
[342]	111	! Initialization
[161]	112	! **********************************************************************
	113
[342]	114	lmax(:)=0.
[161]	115	pdu_th(:,:)=0.
	116	pdv_th(:,:)=0.
	117	pdt_th(:,:)=0.
	118	entr_therm(:,:)=0.
	119	detr_therm(:,:)=0.
	120	q2_therm(:,:)=0.
	121	dq2_therm(:,:)=0.
	122	ztla(:,:)=0.
	123	pbl_dtke(:,:)=0.
	124	fm_therm(:,:)=0.
	125	zw2(:,:)=0.
	126	fraca(:,:)=0.
[512]	127	zfraca(:,:)=0.
[161]	128	if (tracer) then
	129	pdq_th(:,:,:)=0.
	130	end if
[342]	131	d_t_ajs(:,:)=0.
	132	d_u_ajs(:,:)=0.
	133	d_v_ajs(:,:)=0.
	134	d_q_ajs(:,:,:)=0.
	135	heatFlux(:,:)=0.
	136	heatFlux_down(:,:)=0.
	137	buoyancyOut(:,:)=0.
	138	buoyancyEst(:,:)=0.
	139	zmaxth(:)=0.
	140	lmax_real(:)=0.
[161]	141
	142
[342]	143	! **********************************************************************
	144	! Preparing inputs for the thermals
	145	! **********************************************************************
[161]	146
[342]	147	zu(:,:)=pu(:,:)+pdu(:,:)*ptimestep
	148	zv(:,:)=pv(:,:)+pdv(:,:)*ptimestep
	149	zt(:,:)=pt(:,:)+pdt(:,:)*ptimestep
[161]	150
[342]	151	pq_therm(:,:,:)=0.
	152	qtransport_thermals=.true. !! default setting
	153	!call getin("qtransport_thermals",qtransport_thermals)
[161]	154
[342]	155	if(qtransport_thermals) then
	156	if(tracer) then
	157	pq_therm(:,:,:)=pq(:,:,:)+pdq(:,:,:)*ptimestep
	158	endif
	159	endif
[161]	160
[544]	161	dtke_thermals=.false. !! default setting
	162	call getin("dtke_thermals",dtke_thermals)
	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
[342]	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
[561]	194	! r_aspect_thermals ! Mainly control the shape of the temperature profile
	195	! in the surface layer. Decreasing it goes toward
	196	! a convective-adjustment like profile.
	197	! nsplit_thermals ! Sub-timestep for the thermals. Very dependant on the
	198	! chosen timestep for the radiative transfer.
	199	! It is recommended to run with 96 timestep per day and
	200	! iradia = 1., configuration in which thermals can run
	201	! very well with a sub-timestep of 10.
	202	IF (firstcall) THEN
[342]	203	#ifdef MESOSCALE
[561]	204	!! valid for timesteps < 200s
	205	nsplit_thermals=4
	206	r_aspect_thermals=0.7
[342]	207	#else
[561]	208	IF ((ptimestep .lt. 3699.*24./96.) .and. (iradia .eq. 1)) THEN
	209	nsplit_thermals=10
	210	ELSE
	211	nsplit_thermals=35
	212	ENDIF
	213	r_aspect_thermals=1.5
[342]	214	#endif
[561]	215	call getin("nsplit_thermals",nsplit_thermals)
	216	call getin("r_aspect_thermals",r_aspect_thermals)
	217	ENDIF
[342]	218
	219	! **********************************************************************
	220	! SUB-TIMESTEP LOOP
	221	! **********************************************************************
	222
	223	zdt=ptimestep/REAL(nsplit_thermals)
	224
	225	DO isplit=1,nsplit_thermals
	226
	227	! Initialization of intermediary variables
	228
	229	zfm_therm(:,:)=0.
	230	zentr_therm(:,:)=0.
	231	zdetr_therm(:,:)=0.
	232	zheatFlux(:,:)=0.
	233	zheatFlux_down(:,:)=0.
[508]	234	zbuoyancyOut(:,:)=0.
	235	zbuoyancyEst(:,:)=0.
[342]	236	zzw2(:,:)=0.
	237	zmax(:)=0.
	238	lmax(:)=0.
	239	d_t_the(:,:)=0.
	240	d_u_the(:,:)=0.
	241	d_v_the(:,:)=0.
	242	dq2_the(:,:)=0.
	243	if (nqmx .ne. 0) then
	244	d_q_the(:,:,:)=0.
[161]	245	endif
	246
[342]	247	CALL thermcell_main_mars(zdt &
	248	& ,pplay,pplev,pphi,zzlev,zzlay &
	249	& ,zu,zv,zt,pq_therm,q2_therm &
	250	& ,d_u_the,d_v_the,d_t_the,d_q_the,dq2_the &
	251	& ,zfm_therm,zentr_therm,zdetr_therm,lmax,zmax &
	252	& ,r_aspect_thermals &
	253	& ,zzw2,fraca,zpopsk &
	254	& ,ztla,zheatFlux,zheatFlux_down &
	255	& ,zbuoyancyOut,zbuoyancyEst)
[161]	256
[342]	257	fact=1./REAL(nsplit_thermals)
[161]	258
[342]	259	d_t_the(:,:)=d_t_the(:,:)ptimestepfact
	260	! d_u_the(:,:)=d_u_the(:,:)*fact
	261	! d_v_the(:,:)=d_v_the(:,:)*fact
	262	! dq2_the(:,:)=dq2_the(:,:)*fact
[508]	263	if (ico2 .ne. 0) then
	264	d_q_the(:,:,ico2)=d_q_the(:,:,ico2)*fact
	265	endif
[161]	266
[342]	267	zmaxth(:)=zmaxth(:)+zmax(:)*fact
	268	lmax_real(:)=lmax_real(:)+float(lmax(:))*fact
	269	fm_therm(:,:)=fm_therm(:,:) &
	270	& +zfm_therm(:,:)*fact
	271	entr_therm(:,:)=entr_therm(:,:) &
	272	& +zentr_therm(:,:)*fact
	273	detr_therm(:,:)=detr_therm(:,:) &
	274	& +zdetr_therm(:,:)*fact
[512]	275	zfraca(:,:)=zfraca(:,:) + fraca(:,:)*fact
[342]	276
	277	heatFlux(:,:)=heatFlux(:,:) &
	278	& +zheatFlux(:,:)*fact
	279	heatFlux_down(:,:)=heatFlux_down(:,:) &
	280	& +zheatFlux_down(:,:)*fact
[508]	281	buoyancyOut(:,:)=buoyancyOut(:,:) &
	282	& +zbuoyancyOut(:,:)*fact
	283	buoyancyEst(:,:)=buoyancyEst(:,:) &
	284	& +zbuoyancyEst(:,:)*fact
[512]	285
[342]	286
	287	zw2(:,:)=zw2(:,:) + zzw2(:,:)*fact
	288
	289	! accumulation de la tendance
	290
	291	d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_the(:,:)
	292	! d_u_ajs(:,:)=d_u_ajs(:,:)+d_u_the(:,:)
	293	! d_v_ajs(:,:)=d_v_ajs(:,:)+d_v_the(:,:)
[508]	294	if (ico2 .ne. 0) then
	295	d_q_ajs(:,:,ico2)=d_q_ajs(:,:,ico2)+d_q_the(:,:,ico2)
	296	endif
[342]	297	! dq2_therm(:,:)=dq2_therm(:,:)+dq2_the(:,:)
	298	! incrementation des variables meteo
	299
	300	zt(:,:) = zt(:,:) + d_t_the(:,:)
	301	! zu(:,:) = zu(:,:) + d_u_the(:,:)
	302	! zv(:,:) = zv(:,:) + d_v_the(:,:)
[508]	303	if (ico2 .ne. 0) then
	304	pq_therm(:,:,ico2) = &
	305	& pq_therm(:,:,ico2) + d_q_the(:,:,ico2)*ptimestep
	306	endif
[342]	307	! q2_therm(:,:) = q2_therm(:,:) + dq2_therm(:,:)
	308
	309
	310	ENDDO ! isplit
	311	!****************************************************************
	312
	313	! Now that we have computed total entrainment and detrainment, we can
	314	! advect u, v, and q in thermals. (theta already advected). We can do
	315	! that separatly because u,v,and q are not used in thermcell_main for
	316	! any thermals-related computation : they are purely passive.
	317
	318	! mass of cells
	319	do l=1,nlayermx
	320	masse(:,l)=(pplev(:,l)-pplev(:,l+1))/g
	321	enddo
	322
	323	! thickness of layers
	324	do l=1,nlayermx
	325	zdz(:,l)=zzlev(:,l+1)-zzlev(:,l)
	326	enddo
	327
	328	modname='momentum'
	329	call thermcell_dqup(ngridmx,nlayermx,ptimestep &
	330	& ,fm_therm,entr_therm,detr_therm, &
	331	& masse,zu,d_u_ajs,modname,zdz)
	332
	333	call thermcell_dqup(ngridmx,nlayermx,ptimestep &
	334	& ,fm_therm,entr_therm,detr_therm, &
	335	& masse,zv,d_v_ajs,modname,zdz)
	336
	337	if (nqmx .ne. 0.) then
	338	modname='tracer'
	339	DO iq=1,nqmx
[508]	340	if (iq .ne. ico2) then
[342]	341	call thermcell_dqup(ngridmx,nlayermx,ptimestep &
	342	& ,fm_therm,entr_therm,detr_therm, &
	343	& masse,pq_therm(:,:,iq),d_q_ajs(:,:,iq),modname,zdz)
[508]	344	endif
[342]	345	ENDDO
	346	endif
	347
[544]	348	if (dtke_thermals) then
	349	modname='tke'
	350	call thermcell_dqup(ngridmx,nlayermx,ptimestep &
	351	& ,fm_therm,entr_therm,detr_therm, &
	352	& masse,q2_therm,dq2_therm,modname,zdz)
	353	endif
	354
[342]	355	DO ig=1,ngridmx
	356	hfmax(ig)=MAXVAL(heatFlux(ig,:)+heatFlux_down(ig,:))
	357	wmax(ig)=MAXVAL(zw2(ig,:))
	358	ENDDO
	359
	360	lmax(:)=nint(lmax_real(:))
	361
	362	! **********************************************************************
	363	! **********************************************************************
	364	! **********************************************************************
	365	! CALLTHERM END
	366	! **********************************************************************
	367	! **********************************************************************
	368	! **********************************************************************
	369
	370
	371	! **********************************************************************
	372	! Preparing outputs
	373	! **********************************************************************
	374
	375	! Winds and tracers PDU, PDV, and PDQ are DERIVATIVES
	376
	377	pdu_th(:,:)=d_u_ajs(:,:)
	378	pdv_th(:,:)=d_v_ajs(:,:)
	379
[161]	380	if(qtransport_thermals) then
[342]	381	if(tracer) then
	382	pdq_th(:,:,:)=d_q_ajs(:,:,:)
	383	endif
[161]	384	endif
	385
[544]	386	IF(dtke_thermals) THEN
	387	DO l=2,nlayermx
	388	pbl_dtke(:,l)=0.5*(dq2_therm(:,l-1)+dq2_therm(:,l))
	389	ENDDO
	390
	391	pbl_dtke(:,1)=0.5*dq2_therm(:,1)
	392	pbl_dtke(:,nlayermx+1)=0.
	393	ENDIF
[161]	394
	395
[342]	396	! Temperature PDT is a TENDANCY
	397	pdt_th(:,:)=d_t_ajs(:,:)/ptimestep
	398
[499]	399
[342]	400	! **********************************************************************
[499]	401	! Compute the free convection velocity scale for vdifc
	402	! **********************************************************************
	403
	404
	405	! Potential temperature gradient
	406
	407	dteta(:,nlayermx)=0.
	408	DO l=1,nlayermx-1
	409	DO ig=1, ngridmx
	410	dteta(ig,l) = ((zt(ig,l+1)-zt(ig,l))/zpopsk(ig,l)) &
	411	& /(zzlay(ig,l+1)-zzlay(ig,l))
	412	ENDDO
	413	ENDDO
	414
	415	! Computation of the pbl mixed layer temperature
	416
	417	DO ig=1, ngridmx
	418	ii=MINLOC(abs(dteta(ig,1:lmax(ig))))
	419	pbl_teta(ig) = zt(ig,ii(1))/zpopsk(ig,ii(1))
	420	ENDDO
	421
	422	! We follow Spiga et. al 2010 (QJRMS)
	423	! ------------
	424
	425	DO ig=1, ngridmx
	426	IF (zmax(ig) .gt. 0.) THEN
	427	wstar(ig)=(gzmaxth(ig)hfmax(ig)/pbl_teta(ig))**(1./3.)
	428	ELSE
	429	wstar(ig)=0.
	430	ENDIF
	431	ENDDO
	432
	433
	434
	435	! **********************************************************************
[342]	436	! Diagnostics
	437	! **********************************************************************
[161]	438
	439	if(outptherm) then
[185]	440	if (ngridmx .eq. 1) then
	441	call WRITEDIAGFI(ngridmx,'entr_therm','entrainement thermique',&
[161]	442	& 'kg/m-2',1,entr_therm)
[185]	443	call WRITEDIAGFI(ngridmx,'detr_therm','detrainement thermique',&
[161]	444	& 'kg/m-2',1,detr_therm)
[185]	445	call WRITEDIAGFI(ngridmx,'fm_therm','flux masse thermique',&
[161]	446	& 'kg/m-2',1,fm_therm)
[185]	447	call WRITEDIAGFI(ngridmx,'zw2','vitesse verticale thermique',&
[161]	448	& 'm/s',1,zw2)
[185]	449	call WRITEDIAGFI(ngridmx,'heatFlux_up','heatFlux_updraft',&
[161]	450	& 'SI',1,heatFlux)
[185]	451	call WRITEDIAGFI(ngridmx,'heatFlux_down','heatFlux_downdraft',&
[161]	452	& 'SI',1,heatFlux_down)
[185]	453	call WRITEDIAGFI(ngridmx,'fraca','fraction coverage',&
[161]	454	& 'percent',1,fraca)
[185]	455	call WRITEDIAGFI(ngridmx,'buoyancyOut','buoyancyOut',&
[161]	456	& 'm.s-2',1,buoyancyOut)
[185]	457	call WRITEDIAGFI(ngridmx,'buoyancyEst','buoyancyEst',&
[161]	458	& 'm.s-2',1,buoyancyEst)
[185]	459	call WRITEDIAGFI(ngridmx,'d_t_th', &
[161]	460	& 'tendance temp TH','K',1,d_t_ajs)
[185]	461	call WRITEDIAGFI(ngridmx,'zmax', &
[342]	462	& 'pbl height','m',0,zmaxth)
[161]	463	else
	464
[185]	465	call WRITEDIAGFI(ngridmx,'entr_therm','entrainement thermique',&
[161]	466	& 'kg/m-2',3,entr_therm)
[185]	467	call WRITEDIAGFI(ngridmx,'detr_therm','detrainement thermique',&
[161]	468	& 'kg/m-2',3,detr_therm)
[185]	469	call WRITEDIAGFI(ngridmx,'fm_therm','flux masse thermique',&
[161]	470	& 'kg/m-2',3,fm_therm)
[185]	471	call WRITEDIAGFI(ngridmx,'zw2','vitesse verticale thermique',&
[161]	472	& 'm/s',3,zw2)
[185]	473	call WRITEDIAGFI(ngridmx,'heatFlux','heatFlux',&
[161]	474	& 'SI',3,heatFlux)
[185]	475	call WRITEDIAGFI(ngridmx,'buoyancyOut','buoyancyOut',&
[161]	476	& 'SI',3,buoyancyOut)
[185]	477	call WRITEDIAGFI(ngridmx,'d_t_th', &
[161]	478	& 'tendance temp TH','K',3,d_t_ajs)
	479
	480	endif
	481	endif
	482
	483	END

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

Download in other formats:

Original Format