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

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	#include "tracer.h"
18
19	!--------------------------------------------------------
20	! Input Variables
21	!--------------------------------------------------------
22
23	REAL, INTENT(IN) :: ptimestep
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)
30	LOGICAL, INTENT(IN) :: firstcall
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)
36
37	!--------------------------------------------------------
38	! Output Variables
39	!--------------------------------------------------------
40
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)
48	REAL, INTENT(OUT) :: wstar(ngridmx)
49
50	!--------------------------------------------------------
51	! Thermals local variables
52	!--------------------------------------------------------
53	REAL zu(ngridmx,nlayermx), zv(ngridmx,nlayermx)
54	REAL zt(ngridmx,nlayermx)
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)
61	REAL fraca(ngridmx,nlayermx+1),zfraca(ngridmx,nlayermx+1)
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)
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
80	INTEGER,SAVE :: nsplit_thermals
81	REAL, SAVE :: r_aspect_thermals
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
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)
101	REAL pbl_teta(ngridmx),dteta(ngridmx,nlayermx)
102
103	!--------------------------------------------------------
104	! Theta_m
105	!--------------------------------------------------------
106
107	INTEGER ico2
108	SAVE ico2
109
110	! **********************************************************************
111	! Initialization
112	! **********************************************************************
113
114	lmax(:)=0.
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.
127	zfraca(:,:)=0.
128	if (tracer) then
129	pdq_th(:,:,:)=0.
130	end if
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.
141
142
143	! **********************************************************************
144	! Preparing inputs for the thermals
145	! **********************************************************************
146
147	zu(:,:)=pu(:,:)+pdu(:,:)*ptimestep
148	zv(:,:)=pv(:,:)+pdv(:,:)*ptimestep
149	zt(:,:)=pt(:,:)+pdt(:,:)*ptimestep
150
151	pq_therm(:,:,:)=0.
152	qtransport_thermals=.true. !! default setting
153	!call getin("qtransport_thermals",qtransport_thermals)
154
155	if(qtransport_thermals) then
156	if(tracer) then
157	pq_therm(:,:,:)=pq(:,:,:)+pdq(:,:,:)*ptimestep
158	endif
159	endif
160
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
168
169	! **********************************************************************
170	! Polar night mixing : theta_m
171	! **********************************************************************
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
186	! **********************************************************************
187	! **********************************************************************
188	! **********************************************************************
189	! CALLTHERM
190	! **********************************************************************
191	! **********************************************************************
192	! **********************************************************************
193
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
203	#ifdef MESOSCALE
204	!! valid for timesteps < 200s
205	nsplit_thermals=4
206	r_aspect_thermals=0.7
207	#else
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
214	#endif
215	call getin("nsplit_thermals",nsplit_thermals)
216	call getin("r_aspect_thermals",r_aspect_thermals)
217	ENDIF
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.
234	zbuoyancyOut(:,:)=0.
235	zbuoyancyEst(:,:)=0.
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.
245	endif
246
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)
256
257	fact=1./REAL(nsplit_thermals)
258
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
263	if (ico2 .ne. 0) then
264	d_q_the(:,:,ico2)=d_q_the(:,:,ico2)*fact
265	endif
266
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
275	zfraca(:,:)=zfraca(:,:) + fraca(:,:)*fact
276
277	heatFlux(:,:)=heatFlux(:,:) &
278	& +zheatFlux(:,:)*fact
279	heatFlux_down(:,:)=heatFlux_down(:,:) &
280	& +zheatFlux_down(:,:)*fact
281	buoyancyOut(:,:)=buoyancyOut(:,:) &
282	& +zbuoyancyOut(:,:)*fact
283	buoyancyEst(:,:)=buoyancyEst(:,:) &
284	& +zbuoyancyEst(:,:)*fact
285
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(:,:)
294	if (ico2 .ne. 0) then
295	d_q_ajs(:,:,ico2)=d_q_ajs(:,:,ico2)+d_q_the(:,:,ico2)
296	endif
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(:,:)
303	if (ico2 .ne. 0) then
304	pq_therm(:,:,ico2) = &
305	& pq_therm(:,:,ico2) + d_q_the(:,:,ico2)*ptimestep
306	endif
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
340	if (iq .ne. ico2) then
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)
344	endif
345	ENDDO
346	endif
347
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
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
380	if(qtransport_thermals) then
381	if(tracer) then
382	pdq_th(:,:,:)=d_q_ajs(:,:,:)
383	endif
384	endif
385
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
394
395
396	! Temperature PDT is a TENDANCY
397	pdt_th(:,:)=d_t_ajs(:,:)/ptimestep
398
399
400	! **********************************************************************
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	! **********************************************************************
436	! Diagnostics
437	! **********************************************************************
438
439	if(outptherm) then
440	if (ngridmx .eq. 1) then
441	call WRITEDIAGFI(ngridmx,'entr_therm','entrainement thermique',&
442	& 'kg/m-2',1,entr_therm)
443	call WRITEDIAGFI(ngridmx,'detr_therm','detrainement thermique',&
444	& 'kg/m-2',1,detr_therm)
445	call WRITEDIAGFI(ngridmx,'fm_therm','flux masse thermique',&
446	& 'kg/m-2',1,fm_therm)
447	call WRITEDIAGFI(ngridmx,'zw2','vitesse verticale thermique',&
448	& 'm/s',1,zw2)
449	call WRITEDIAGFI(ngridmx,'heatFlux_up','heatFlux_updraft',&
450	& 'SI',1,heatFlux)
451	call WRITEDIAGFI(ngridmx,'heatFlux_down','heatFlux_downdraft',&
452	& 'SI',1,heatFlux_down)
453	call WRITEDIAGFI(ngridmx,'fraca','fraction coverage',&
454	& 'percent',1,fraca)
455	call WRITEDIAGFI(ngridmx,'buoyancyOut','buoyancyOut',&
456	& 'm.s-2',1,buoyancyOut)
457	call WRITEDIAGFI(ngridmx,'buoyancyEst','buoyancyEst',&
458	& 'm.s-2',1,buoyancyEst)
459	call WRITEDIAGFI(ngridmx,'d_t_th', &
460	& 'tendance temp TH','K',1,d_t_ajs)
461	call WRITEDIAGFI(ngridmx,'zmax', &
462	& 'pbl height','m',0,zmaxth)
463	else
464
465	call WRITEDIAGFI(ngridmx,'entr_therm','entrainement thermique',&
466	& 'kg/m-2',3,entr_therm)
467	call WRITEDIAGFI(ngridmx,'detr_therm','detrainement thermique',&
468	& 'kg/m-2',3,detr_therm)
469	call WRITEDIAGFI(ngridmx,'fm_therm','flux masse thermique',&
470	& 'kg/m-2',3,fm_therm)
471	call WRITEDIAGFI(ngridmx,'zw2','vitesse verticale thermique',&
472	& 'm/s',3,zw2)
473	call WRITEDIAGFI(ngridmx,'heatFlux','heatFlux',&
474	& 'SI',3,heatFlux)
475	call WRITEDIAGFI(ngridmx,'buoyancyOut','buoyancyOut',&
476	& 'SI',3,buoyancyOut)
477	call WRITEDIAGFI(ngridmx,'d_t_th', &
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