Context Navigation

calltherm_interface.F90 @ 623

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

Thermals: since the correction of the advection scheme, it seems like momentum is too efficiently transported in the surface layer, leading to an instability that will crash the model in 3D. I am very sorry that I did not see that before ! So, I deactived the momentum transport in thermals temporarily so that I can find a fix. It seems like this propblem is not present on Earth, as the solution they adopt for that problem is not working for Mars.

File size: 16.3 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	LOGICAL qtransport_thermals,dtke_thermals
68	INTEGER l,ig,iq,ii(1)
69
70	!--------------------------------------------------------
71	! Local variables for sub-timestep
72	!--------------------------------------------------------
73
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)
77	INTEGER isplit
78	INTEGER,SAVE :: nsplit_thermals
79	REAL, SAVE :: 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
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)
99	REAL pbl_teta(ngridmx),dteta(ngridmx,nlayermx)
100
101	!--------------------------------------------------------
102	! Theta_m
103	!--------------------------------------------------------
104
105	INTEGER ico2
106	SAVE ico2
107
108	! **********************************************************************
109	! Initialization
110	! **********************************************************************
111
112	lmax(:)=0
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.
125	zfraca(:,:)=0.
126	if (tracer) then
127	pdq_th(:,:,:)=0.
128	end if
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.
139
140
141	! **********************************************************************
142	! Preparing inputs for the thermals
143	! **********************************************************************
144
145	zu(:,:)=pu(:,:)+pdu(:,:)*ptimestep
146	zv(:,:)=pv(:,:)+pdv(:,:)*ptimestep
147	zt(:,:)=pt(:,:)+pdt(:,:)*ptimestep
148
149	pq_therm(:,:,:)=0.
150	qtransport_thermals=.true. !! default setting
151	!call getin("qtransport_thermals",qtransport_thermals)
152
153	if(qtransport_thermals) then
154	if(tracer) then
155	pq_therm(:,:,:)=pq(:,:,:)+pdq(:,:,:)*ptimestep
156	endif
157	endif
158
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
166
167	! **********************************************************************
168	! Polar night mixing : theta_m
169	! **********************************************************************
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
184	! **********************************************************************
185	! **********************************************************************
186	! **********************************************************************
187	! CALLTHERM
188	! **********************************************************************
189	! **********************************************************************
190	! **********************************************************************
191
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
201	#ifdef MESOSCALE
202	!! valid for timesteps < 200s
203	nsplit_thermals=4
204	r_aspect_thermals=0.7
205	#else
206	IF ((ptimestep .le. 3699.*24./96.) .and. (iradia .eq. 1)) THEN
207	nsplit_thermals=10
208	ELSE
209	nsplit_thermals=35
210	ENDIF
211	r_aspect_thermals=1.
212	#endif
213	call getin("nsplit_thermals",nsplit_thermals)
214	call getin("r_aspect_thermals",r_aspect_thermals)
215	ENDIF
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.
232	zbuoyancyOut(:,:)=0.
233	zbuoyancyEst(:,:)=0.
234	zzw2(:,:)=0.
235	zmax(:)=0.
236	lmax(:)=0
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.
243	endif
244
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)
254
255	fact=1./REAL(nsplit_thermals)
256
257	d_t_the(:,:)=d_t_the(:,:)ptimestepfact
258	! d_u_the(:,:)=d_u_the(:,:)*fact
259	! d_v_the(:,:)=d_v_the(:,:)*fact
260	! dq2_the(:,:)=dq2_the(:,:)*fact
261	if (ico2 .ne. 0) then
262	d_q_the(:,:,ico2)=d_q_the(:,:,ico2)*fact
263	endif
264
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
273	zfraca(:,:)=zfraca(:,:) + fraca(:,:)*fact
274
275	heatFlux(:,:)=heatFlux(:,:) &
276	& +zheatFlux(:,:)*fact
277	heatFlux_down(:,:)=heatFlux_down(:,:) &
278	& +zheatFlux_down(:,:)*fact
279	buoyancyOut(:,:)=buoyancyOut(:,:) &
280	& +zbuoyancyOut(:,:)*fact
281	buoyancyEst(:,:)=buoyancyEst(:,:) &
282	& +zbuoyancyEst(:,:)*fact
283
284
285	zw2(:,:)=zw2(:,:) + zzw2(:,:)*fact
286
287	! accumulation de la tendance
288
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(:,:)
292	if (ico2 .ne. 0) then
293	d_q_ajs(:,:,ico2)=d_q_ajs(:,:,ico2)+d_q_the(:,:,ico2)
294	endif
295	! dq2_therm(:,:)=dq2_therm(:,:)+dq2_the(:,:)
296	! incrementation des variables meteo
297
298	zt(:,:) = zt(:,:) + d_t_the(:,:)
299	! zu(:,:) = zu(:,:) + d_u_the(:,:)
300	! zv(:,:) = zv(:,:) + d_v_the(:,:)
301	if (ico2 .ne. 0) then
302	pq_therm(:,:,ico2) = &
303	& pq_therm(:,:,ico2) + d_q_the(:,:,ico2)*ptimestep
304	endif
305	! q2_therm(:,:) = q2_therm(:,:) + dq2_therm(:,:)
306
307
308	ENDDO ! isplit
309	!****************************************************************
310
311	lmax(:)=nint(lmax_real(:))
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	! 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)
330
331	if (nqmx .ne. 0.) then
332	DO iq=1,nqmx
333	if (iq .ne. ico2) then
334	call thermcell_dqup(ngridmx,nlayermx,ptimestep &
335	& ,fm_therm,entr_therm, &
336	& masse,pq_therm(:,:,iq),d_q_ajs(:,:,iq))
337	endif
338	ENDDO
339	endif
340
341	if (dtke_thermals) then
342	call thermcell_dqup(ngridmx,nlayermx,ptimestep &
343	& ,fm_therm,entr_therm, &
344	& masse,q2_therm,dq2_therm)
345	endif
346
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
365	! Winds and tracers PDU, PDV, and PDQ are DERIVATIVES
366
367	pdu_th(:,:)=d_u_ajs(:,:)
368	pdv_th(:,:)=d_v_ajs(:,:)
369
370	if(qtransport_thermals) then
371	if(tracer) then
372	pdq_th(:,:,:)=d_q_ajs(:,:,:)
373	endif
374	endif
375
376	IF(dtke_thermals) THEN
377	DO l=2,nlayermx
378	pbl_dtke(:,l)=0.5*(dq2_therm(:,l-1)+dq2_therm(:,l))
379	ENDDO
380
381	pbl_dtke(:,1)=0.5*dq2_therm(:,1)
382	pbl_dtke(:,nlayermx+1)=0.
383	ENDIF
384
385
386	! Temperature PDT is a TENDANCY
387	pdt_th(:,:)=d_t_ajs(:,:)/ptimestep
388
389
390	! **********************************************************************
391	! Compute the free convection velocity scale for vdifc
392	! **********************************************************************
393
394
395	! Potential temperature gradient
396
397	dteta(:,nlayermx)=0.
398	DO l=1,nlayermx-1
399	DO ig=1, ngridmx
400	dteta(ig,l) = ((zt(ig,l+1)-zt(ig,l))/zpopsk(ig,l)) &
401	& /(zzlay(ig,l+1)-zzlay(ig,l))
402	ENDDO
403	ENDDO
404
405	! Computation of the pbl mixed layer temperature
406
407	DO ig=1, ngridmx
408	ii=MINLOC(abs(dteta(ig,1:lmax(ig))))
409	pbl_teta(ig) = zt(ig,ii(1))/zpopsk(ig,ii(1))
410	ENDDO
411
412	! We follow Spiga et. al 2010 (QJRMS)
413	! ------------
414
415	DO ig=1, ngridmx
416	IF (zmax(ig) .gt. 0.) THEN
417	wstar(ig)=(gzmaxth(ig)hfmax(ig)/pbl_teta(ig))**(1./3.)
418	ELSE
419	wstar(ig)=0.
420	ENDIF
421	ENDDO
422
423
424
425	! **********************************************************************
426	! Diagnostics
427	! **********************************************************************
428
429	if(outptherm) then
430	if (ngridmx .eq. 1) then
431	call WRITEDIAGFI(ngridmx,'entr_therm','entrainement thermique',&
432	& 'kg/m-2',1,entr_therm)
433	call WRITEDIAGFI(ngridmx,'detr_therm','detrainement thermique',&
434	& 'kg/m-2',1,detr_therm)
435	call WRITEDIAGFI(ngridmx,'fm_therm','flux masse thermique',&
436	& 'kg/m-2',1,fm_therm)
437	call WRITEDIAGFI(ngridmx,'zw2','vitesse verticale thermique',&
438	& 'm/s',1,zw2)
439	call WRITEDIAGFI(ngridmx,'heatFlux_up','heatFlux_updraft',&
440	& 'SI',1,heatFlux)
441	call WRITEDIAGFI(ngridmx,'heatFlux_down','heatFlux_downdraft',&
442	& 'SI',1,heatFlux_down)
443	call WRITEDIAGFI(ngridmx,'fraca','fraction coverage',&
444	& 'percent',1,fraca)
445	call WRITEDIAGFI(ngridmx,'buoyancyOut','buoyancyOut',&
446	& 'm.s-2',1,buoyancyOut)
447	call WRITEDIAGFI(ngridmx,'buoyancyEst','buoyancyEst',&
448	& 'm.s-2',1,buoyancyEst)
449	call WRITEDIAGFI(ngridmx,'d_t_th', &
450	& 'tendance temp TH','K',1,d_t_ajs)
451	call WRITEDIAGFI(ngridmx,'d_q_th', &
452	& 'tendance traceur TH','kg/kg',1,d_q_ajs)
453	call WRITEDIAGFI(ngridmx,'zmax', &
454	& 'pbl height','m',0,zmaxth)
455	else
456
457	call WRITEDIAGFI(ngridmx,'entr_therm','entrainement thermique',&
458	& 'kg/m-2',3,entr_therm)
459	call WRITEDIAGFI(ngridmx,'detr_therm','detrainement thermique',&
460	& 'kg/m-2',3,detr_therm)
461	call WRITEDIAGFI(ngridmx,'fm_therm','flux masse thermique',&
462	& 'kg/m-2',3,fm_therm)
463	call WRITEDIAGFI(ngridmx,'zw2','vitesse verticale thermique',&
464	& 'm/s',3,zw2)
465	call WRITEDIAGFI(ngridmx,'heatFlux','heatFlux',&
466	& 'SI',3,heatFlux)
467	call WRITEDIAGFI(ngridmx,'buoyancyOut','buoyancyOut',&
468	& 'SI',3,buoyancyOut)
469	call WRITEDIAGFI(ngridmx,'d_t_th', &
470	& 'tendance temp TH','K',3,d_t_ajs)
471
472	endif
473	endif
474
475	END

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

Download in other formats:

Original Format