Context Navigation

calltherm_mars.F90 @ 338

Last change on this file since 338 was 337, checked in by acolaitis, 14 years ago

10/10/2011 == AC

*
This commit aims at increasing the thermals speed, especially for large tracer number configurations. The idea behind this commit is to advect non-active conserved variables outside of the sub-timestep of the thermals. Because these variables are not used in thermals computation, we can decouple them:

momentum: can be decoupled because we assume a constant ratio between horizontal velocity in alimentation layer and maximum vertical velocity in the thermals.

tracer: can be decoupled because we do not take condensation of any tracer into account and hence do not liberate latent heat nor form clouds in the thermals.

temperature: cannot be decoupled (of course)
*

D 336 libf/phymars/thermcell_dqupdown.F90
^{---------------- Deleted and replaced by a simpler version. Notes about downdraft advection are still available from revision 336 of SVN in thermcell_dqupdown.}

A 0 libf/phymars/thermcell_dqup.F90
^{---------------- New upward advection for tracers and momentum in thermals. Several changes are done compared to the old approach:}

Updraft quantities are not longer computed by making hypothesis on the amount of advected air.
- In general, the formalism for updraft computation is much simpler and clearer.
Tracer tendancies are no longer computed using the conservation equation. Instead, we use the divergence of an approximated turbulent flux of the concerned quantity, where downdraft are also neglected.

M 336 libf/phymars/thermcell_main_mars.F90
^{---------------- The Main does not call anymore thermcell_dqupdown, which it was doing 2+tracer_number times per subtimestep (140 times per physical step for a 2 tracer config)

M 336 libf/phymars/calltherm_mars.F90}---------------- Entrainment, detrainment and mass-fluxes are recomputed in the sub-timestep loop. Their final value after iterations is used by the new advection routine to compute tracer and momentum fluxes.

* Results

Conservation of tracers has been assessed over 1 yr in 1D and found to be comparable to that obtained with the simple convective adjustment. (it actually seems to be better by a factor of 10%!)
GCM speed-up is of about 20% compared to the old thermal configuration, for a 2 tracer case.
Advection of sharp tracer profiles has been successfully observed, similar to the old method.

File size: 8.2 KB

Line
1	!
2	! $Id: calltherm.F90 1428 2010-09-13 08:43:37Z fairhead $
3	!
4	subroutine calltherm_mars(ptimestep,zzlev,zzlay &
5	& ,pplay,pplev,pphi &
6	& ,u_seri,v_seri,t_seri,pq_therm,q2_therm &
7	& ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs,dq2_therm &
8	& ,fm_therm,entr_therm,detr_therm,lmax,zmaxth,&
9	& zw2,fraca,zpopsk,ztla,heatFlux,heatFlux_down,&
10	& buoyancyOut,buoyancyEst,hfmax,wmax)
11
12	USE ioipsl_getincom
13	implicit none
14
15	#include "dimensions.h"
16	#include "dimphys.h"
17	#include "comcstfi.h"
18
19	REAL ptimestep
20	LOGICAL logexpr0, logexpr2(ngridmx,nlayermx), logexpr1(ngridmx)
21	REAL fact
22	INTEGER nbptspb,iq,l
23
24	REAL, INTENT(IN) :: zzlay(ngridmx,nlayermx)
25	REAL, INTENT(IN) :: zzlev(ngridmx,nlayermx+1)
26
27	REAL u_seri(ngridmx,nlayermx),v_seri(ngridmx,nlayermx)
28	REAL t_seri(ngridmx,nlayermx),pq_therm(ngridmx,nlayermx,nqmx)
29	REAL q2_therm(ngridmx,nlayermx)
30	REAL pplev(ngridmx,nlayermx+1)
31	REAL pplay(ngridmx,nlayermx)
32	REAL pphi(ngridmx,nlayermx)
33	real zlev(ngridmx,nlayermx+1)
34	!test: on sort lentr et a* pour alimenter KE
35	REAL zw2(ngridmx,nlayermx+1),fraca(ngridmx,nlayermx+1)
36	REAL zzw2(ngridmx,nlayermx+1)
37
38	!FH Update Thermiques
39	REAL d_t_ajs(ngridmx,nlayermx), d_q_ajs(ngridmx,nlayermx,nqmx)
40	REAL d_u_ajs(ngridmx,nlayermx),d_v_ajs(ngridmx,nlayermx)
41	REAL dq2_therm(ngridmx,nlayermx), dq2_the(ngridmx,nlayermx)
42	real fm_therm(ngridmx,nlayermx+1)
43	real entr_therm(ngridmx,nlayermx),detr_therm(ngridmx,nlayermx)
44	REAL masse(ngridmx,nlayermx)
45
46	!********************************************************
47	! declarations
48	real zpopsk(ngridmx,nlayermx)
49	real ztla(ngridmx,nlayermx)
50	real wmax(ngridmx)
51	real hfmax(ngridmx)
52	integer lmax(ngridmx)
53	real lmax_real(ngridmx)
54	real zmax(ngridmx),zmaxth(ngridmx)
55	REAL zdz(ngridmx,nlayermx)
56
57
58	!nouvelles variables pour la convection
59	!RC
60	!on garde le zmax du pas de temps precedent
61	!********************************************************
62
63
64	! variables locales
65	REAL d_t_the(ngridmx,nlayermx), d_q_the(ngridmx,nlayermx,nqmx)
66	REAL d_u_the(ngridmx,nlayermx),d_v_the(ngridmx,nlayermx)
67	!
68	integer isplit,nsplit_thermals
69	real r_aspect_thermals
70
71	real zfm_therm(ngridmx,nlayermx+1),zdt
72	real zentr_therm(ngridmx,nlayermx),zdetr_therm(ngridmx,nlayermx)
73	real heatFlux(ngridmx,nlayermx)
74	real heatFlux_down(ngridmx,nlayermx)
75	real buoyancyOut(ngridmx,nlayermx)
76	real buoyancyEst(ngridmx,nlayermx)
77	real zheatFlux(ngridmx,nlayermx)
78	real zheatFlux_down(ngridmx,nlayermx)
79	real zbuoyancyOut(ngridmx,nlayermx)
80	real zbuoyancyEst(ngridmx,nlayermx)
81
82	character (len=20) :: modname
83	character (len=80) :: abort_message
84
85	integer i,k
86	logical, save :: first=.true.
87
88	REAL tstart,tstop
89
90
91	! Modele du thermique
92	! ===================
93
94	! r_aspect_thermals ! ultimately conrols the amount of mass going through the thermals
95	! decreasing it increases the thermals effect. Tests at gcm resolution
96	! shows that too low values destabilize the model
97	! when changing this value, one should check that the surface layer model
98	! outputs the correct Cdu and Chu through changing the gustiness coefficient beta
99
100
101	#ifdef MESOSCALE
102	!! valid for timesteps < 200s
103	nsplit_thermals=2
104	r_aspect_thermals=0.7
105	#else
106	nsplit_thermals=35
107	r_aspect_thermals=1.5
108	#endif
109
110	call getin("nsplit_thermals",nsplit_thermals)
111	call getin("r_aspect_thermals",r_aspect_thermals)
112
113	fm_therm(:,:)=0.
114	detr_therm(:,:)=0.
115	entr_therm(:,:)=0.
116
117	heatFlux(:,:)=0.
118	heatFlux_down(:,:)=0.
119	! buoyancyOut(:,:)=0.
120	! buoyancyEst(:,:)=0.
121
122	zw2(:,:)=0.
123	zmaxth(:)=0.
124	lmax_real(:)=0.
125
126	zdt=ptimestep/REAL(nsplit_thermals)
127
128	do isplit=1,nsplit_thermals
129
130	! call cpu_time(tstart)
131
132
133	! On reinitialise les flux de masse a zero pour le cumul en
134	! cas de splitting
135
136	zfm_therm(:,:)=0.
137	zentr_therm(:,:)=0.
138	zdetr_therm(:,:)=0.
139	!
140	zheatFlux(:,:)=0.
141	zheatFlux_down(:,:)=0.
142	! zbuoyancyOut(:,:)=0.
143	! zbuoyancyEst(:,:)=0.
144
145	zzw2(:,:)=0.
146	zmax(:)=0.
147	lmax(:)=0.
148
149	d_t_the(:,:)=0.
150	d_u_the(:,:)=0.
151	d_v_the(:,:)=0.
152	! dq2_the(:,:)=0.
153	if (nqmx .ne. 0) then
154	d_q_the(:,:,:)=0.
155	endif
156
157	CALL thermcell_main_mars(zdt &
158	! CALL thermcell_main_mars_coupled_v2(zdt &
159	& ,pplay,pplev,pphi,zzlev,zzlay &
160	& ,u_seri,v_seri,t_seri,pq_therm,q2_therm &
161	& ,d_u_the,d_v_the,d_t_the,d_q_the,dq2_the &
162	& ,zfm_therm,zentr_therm,zdetr_therm,lmax,zmax &
163	& ,r_aspect_thermals &
164	& ,zzw2,fraca,zpopsk &
165	& ,ztla,zheatFlux,zheatFlux_down &
166	& ,zbuoyancyOut,zbuoyancyEst)
167
168	fact=1./REAL(nsplit_thermals)
169	! transformation de la derivee en tendance
170
171	d_t_the(:,:)=d_t_the(:,:)ptimestepfact
172	! d_u_the(:,:)=d_u_the(:,:)*fact
173	! d_v_the(:,:)=d_v_the(:,:)*fact
174	! dq2_the(:,:)=dq2_the(:,:)*fact
175
176	! if (nqmx .ne. 0) then
177	! d_q_the(:,:,:)=d_q_the(:,:,:)*fact
178	! endif
179
180	zmaxth(:)=zmaxth(:)+zmax(:)*fact
181	lmax_real(:)=lmax_real(:)+float(lmax(:))*fact
182	fm_therm(:,:)=fm_therm(:,:) &
183	& +zfm_therm(:,:)*fact
184	entr_therm(:,:)=entr_therm(:,:) &
185	& +zentr_therm(:,:)*fact
186	detr_therm(:,:)=detr_therm(:,:) &
187	& +zdetr_therm(:,:)*fact
188
189	heatFlux(:,:)=heatFlux(:,:) &
190	& +zheatFlux(:,:)*fact
191	heatFlux_down(:,:)=heatFlux_down(:,:) &
192	& +zheatFlux_down(:,:)*fact
193	! buoyancyOut(:,:)=buoyancyOut(:,:) &
194	! & +zbuoyancyOut(:,:)*fact
195	! buoyancyEst(:,:)=buoyancyEst(:,:) &
196	! & +zbuoyancyEst(:,:)*fact
197
198	zw2(:,:)=zw2(:,:) + zzw2(:,:)*fact
199
200	! accumulation de la tendance
201
202	d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_the(:,:)
203	! d_u_ajs(:,:)=d_u_ajs(:,:)+d_u_the(:,:)
204	! d_v_ajs(:,:)=d_v_ajs(:,:)+d_v_the(:,:)
205	! d_q_ajs(:,:,:)=d_q_ajs(:,:,:)+d_q_the(:,:,:)
206	! dq2_therm(:,:)=dq2_therm(:,:)+dq2_the(:,:)
207	! incrementation des variables meteo
208
209	t_seri(:,:) = t_seri(:,:) + d_t_the(:,:)
210	! u_seri(:,:) = u_seri(:,:) + d_u_the(:,:)
211	! v_seri(:,:) = v_seri(:,:) + d_v_the(:,:)
212	! pq_therm(:,:,:) = pq_therm(:,:,:) + d_q_the(:,:,:)
213	! q2_therm(:,:) = q2_therm(:,:) + dq2_therm(:,:)
214
215
216	! call cpu_time(tstop)
217	! print*,'elapsed time in thermals : ',tstop-tstart
218
219	enddo ! isplit
220
221
222	!****************************************************************
223
224	! Now that we have computed total entrainment and detrainment, we can
225	! advect u, v, and q in thermals. (theta already advected). We can do
226	! that separatly because u,v,and q are not used in thermcell_main for
227	! any thermals-related computation : they are purely passive.
228
229	!calcul de la masse
230	do l=1,nlayermx
231	masse(:,l)=(pplev(:,l)-pplev(:,l+1))/g
232	enddo
233
234	!calcul de l'epaisseur des couches
235	do l=1,nlayermx
236	zdz(:,l)=zzlev(:,l+1)-zzlev(:,l)
237	enddo
238
239
240	modname='momentum'
241	call thermcell_dqup(ngridmx,nlayermx,ptimestep &
242	& ,fm_therm,entr_therm,detr_therm, &
243	& masse,u_seri,d_u_ajs,modname,zdz)
244
245	call thermcell_dqup(ngridmx,nlayermx,ptimestep &
246	& ,fm_therm,entr_therm,detr_therm, &
247	& masse,v_seri,d_v_ajs,modname,zdz)
248
249	if (nqmx .ne. 0.) then
250	modname='tracer'
251	DO iq=1,nqmx
252	call thermcell_dqup(ngridmx,nlayermx,ptimestep &
253	& ,fm_therm,entr_therm,detr_therm, &
254	& masse,pq_therm(:,:,iq),d_q_ajs(:,:,iq),modname,zdz)
255
256	ENDDO
257	endif
258
259	DO i=1,ngridmx
260	hfmax(i)=MAXVAL(heatFlux(i,:)+heatFlux_down(i,:))
261	wmax(i)=MAXVAL(zw2(i,:))
262	ENDDO
263
264	lmax(:)=nint(lmax_real(:))
265
266	return
267
268	end

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

Download in other formats:

Original Format