Context Navigation

calltherm_interface.F90 @ 627

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

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

Download in other formats:

Original Format