Context Navigation

calltherm_interface.F90 @ 503

Last change on this file since 503 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

Line
1	!
2	! AC 2011-01-05
3	!
4	SUBROUTINE calltherm_interface (firstcall, &
5	& long,lati,zzlev,zzlay, &
6	& ptimestep,pu,pv,pt,pq,pdu,pdv,pdt,pdq,q2, &
7	& pplay,pplev,pphi,zpopsk, &
8	& pdu_th,pdv_th,pdt_th,pdq_th,lmax,zmaxth,pbl_dtke,hfmax,wstar)
9
10	USE ioipsl_getincom
11
12	implicit none
13	#include "callkeys.h"
14	#include "dimensions.h"
15	#include "dimphys.h"
16	#include "comcstfi.h"
17
18	!--------------------------------------------------------
19	! Input Variables
20	!--------------------------------------------------------
21
22	REAL, INTENT(IN) :: ptimestep
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)
29	LOGICAL, INTENT(IN) :: firstcall
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)
35
36	!--------------------------------------------------------
37	! Output Variables
38	!--------------------------------------------------------
39
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)
47	REAL, INTENT(OUT) :: wstar(ngridmx)
48
49	!--------------------------------------------------------
50	! Thermals local variables
51	!--------------------------------------------------------
52	REAL zu(ngridmx,nlayermx), zv(ngridmx,nlayermx)
53	REAL zt(ngridmx,nlayermx)
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)
65	REAL lmax_real(ngridmx)
66	REAL masse(ngridmx,nlayermx)
67	REAL zdz(ngridmx,nlayermx)
68	LOGICAL qtransport_thermals,dtke_thermals
69	INTEGER l,ig,iq,ii(1)
70	CHARACTER (LEN=20) :: modname
71
72	!--------------------------------------------------------
73	! Local variables for sub-timestep
74	!--------------------------------------------------------
75
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
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)
100	REAL pbl_teta(ngridmx),dteta(ngridmx,nlayermx)
101
102	!---------------------------------------------------------
103
104
105	! **********************************************************************
106	! Initialization
107	! **********************************************************************
108
109	lmax(:)=0.
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
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.
137
138
139	! **********************************************************************
140	! Preparing inputs for the thermals
141	! **********************************************************************
142
143	zu(:,:)=pu(:,:)+pdu(:,:)*ptimestep
144	zv(:,:)=pv(:,:)+pdv(:,:)*ptimestep
145	zt(:,:)=pt(:,:)+pdt(:,:)*ptimestep
146
147	pq_therm(:,:,:)=0.
148	qtransport_thermals=.true. !! default setting
149	!call getin("qtransport_thermals",qtransport_thermals)
150
151	if(qtransport_thermals) then
152	if(tracer) then
153	pq_therm(:,:,:)=pq(:,:,:)+pdq(:,:,:)*ptimestep
154	endif
155	endif
156
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.
217	endif
218
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)
228
229	fact=1./REAL(nsplit_thermals)
230
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
238
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
337	if(qtransport_thermals) then
338	if(tracer) then
339	pdq_th(:,:,:)=d_q_ajs(:,:,:)
340	endif
341	endif
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
351
352
353	! Temperature PDT is a TENDANCY
354	pdt_th(:,:)=d_t_ajs(:,:)/ptimestep
355
356
357	! **********************************************************************
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	! **********************************************************************
393	! Diagnostics
394	! **********************************************************************
395
396	if(outptherm) then
397	if (ngridmx .eq. 1) then
398	call WRITEDIAGFI(ngridmx,'entr_therm','entrainement thermique',&
399	& 'kg/m-2',1,entr_therm)
400	call WRITEDIAGFI(ngridmx,'detr_therm','detrainement thermique',&
401	& 'kg/m-2',1,detr_therm)
402	call WRITEDIAGFI(ngridmx,'fm_therm','flux masse thermique',&
403	& 'kg/m-2',1,fm_therm)
404	call WRITEDIAGFI(ngridmx,'zw2','vitesse verticale thermique',&
405	& 'm/s',1,zw2)
406	call WRITEDIAGFI(ngridmx,'heatFlux_up','heatFlux_updraft',&
407	& 'SI',1,heatFlux)
408	call WRITEDIAGFI(ngridmx,'heatFlux_down','heatFlux_downdraft',&
409	& 'SI',1,heatFlux_down)
410	call WRITEDIAGFI(ngridmx,'fraca','fraction coverage',&
411	& 'percent',1,fraca)
412	call WRITEDIAGFI(ngridmx,'buoyancyOut','buoyancyOut',&
413	& 'm.s-2',1,buoyancyOut)
414	call WRITEDIAGFI(ngridmx,'buoyancyEst','buoyancyEst',&
415	& 'm.s-2',1,buoyancyEst)
416	call WRITEDIAGFI(ngridmx,'d_t_th', &
417	& 'tendance temp TH','K',1,d_t_ajs)
418	call WRITEDIAGFI(ngridmx,'zmax', &
419	& 'pbl height','m',0,zmaxth)
420	else
421
422	call WRITEDIAGFI(ngridmx,'entr_therm','entrainement thermique',&
423	& 'kg/m-2',3,entr_therm)
424	call WRITEDIAGFI(ngridmx,'detr_therm','detrainement thermique',&
425	& 'kg/m-2',3,detr_therm)
426	call WRITEDIAGFI(ngridmx,'fm_therm','flux masse thermique',&
427	& 'kg/m-2',3,fm_therm)
428	call WRITEDIAGFI(ngridmx,'zw2','vitesse verticale thermique',&
429	& 'm/s',3,zw2)
430	call WRITEDIAGFI(ngridmx,'heatFlux','heatFlux',&
431	& 'SI',3,heatFlux)
432	call WRITEDIAGFI(ngridmx,'buoyancyOut','buoyancyOut',&
433	& 'SI',3,buoyancyOut)
434	call WRITEDIAGFI(ngridmx,'d_t_th', &
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