Context Navigation

calltherm_interface.F90 @ 358

Last change on this file since 358 was 342, checked in by acolaitis, 14 years ago

M 341 libf/phymars/calltherm_interface.F90
D 341 libf/phymars/calltherm_mars.F90
^{---------------- Merged calltherm_mars with calltherm_interface for simplicity.}

Cleaned up variables and code

M 341 libf/phymars/thermcell_dqup.F90
^{---------------- dqup now output derivatives and not tendancies}

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

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

Download in other formats:

Original Format