Context Navigation

thermcell_dqupdown.F90 @ 328

Last change on this file since 328 was 313, checked in by acolaitis, 14 years ago

* AC 10/10/2011 *

*
This commit aims at increasing the thermals speed. Using these corrections, gcm performances in 64x48x32 using 1 tracer goes from 27.9% elapsed time in thermals to 18.76%.

*
Additional work needs to be done in tracer advection to gain speed in high tracer number configuration. (tracer advection (but not momentum nor temperature) could be decoupled from sub-timestep, as they do not act on the thermals scheme (water vapor is neglected as we use theta and not theta_v, and radiative effect of dust is not computed in the thermals.))

*
=> TOP 5 of routine contributions to gcm runtime :

Each sample counts as 0.01 seconds.

% cumulative self self total

time seconds seconds calls s/call s/call name
18.76 6.33 6.33 960 0.01 0.01 thermcell_main_mars_
17.19 12.13 5.80 svml_powf4.A
13.72 16.76 4.63 10369 0.00 0.00 filtreg_

3.94 18.09 1.33 intel_new_memset
3.73 19.35 1.26 2880 0.00 0.00 thermcell_dqupdown_

note: thermcell_main_mars_ does call quite a lot power computations (see svml_powf4.A), but this number will not increase with tracer numbers.

*
=> LOG:

M 312 libf/phymars/thermcell_main_mars.F90
^{------------------- removed (commented) computations on buoyancy which is purely diagnostic}

tuned internal convergence loop and added convergence criterion

M 312 libf/phymars/thermcell_dqupdown.F90
^{------------------- removed (commented) downdraft-related if-loops (as we do not advect tracers and momentum in downdrafts for now)}

M 312 libf/phymars/calltherm_mars.F90
^{------------------- removed (commented) diagnostic-related computations}

changed default thermals spliting and aspect ratio
corrected a bug where maximum height was not correctly computed and could result in convective adjustment used in place of thermals
when using certains sets of nsplit and r_aspect (was not happening with the baseline version, so that this correction is transparent to
users)

Property svn:executable set to *

File size: 10.0 KB

Line
1	subroutine thermcell_dqupdown(ngrid,nlayer,ptimestep,fm0,entr0, &
2	& detr0,masse0,q_therm,dq_therm,ztvd,fm_down,ztv,charvar,lmax)
3	implicit none
4
5	! #include "iniprint.h"
6	!=======================================================================
7	!
8	! Calcul du transport verticale dans la couche limite en presence
9	! de "thermiques" explicitement representes
10	! calcul du dq/dt une fois qu'on connait les ascendances
11	! Version modifiee pour prendre les downdrafts a la place de la
12	! subsidence compensatoire
13	!=======================================================================
14
15	#include "dimensions.h"
16	#include "dimphys.h"
17
18	! ============================ INPUTS ============================
19
20	INTEGER, INTENT(IN) :: ngrid,nlayer
21	REAL, INTENT(IN) :: ptimestep
22	REAL, INTENT(IN) :: fm0(ngridmx,nlayermx+1)
23	REAL, INTENT(IN) :: entr0(ngridmx,nlayermx),detr0(ngridmx,nlayermx)
24	REAL, INTENT(IN) :: q_therm(ngridmx,nlayermx)
25	REAL, INTENT(IN) :: fm_down(ngridmx,nlayermx+1)
26	REAL, INTENT(IN) :: ztvd(ngridmx,nlayermx)
27	REAL, INTENT(IN) :: ztv(ngridmx,nlayermx)
28	CHARACTER (LEN=20), INTENT(IN) :: charvar
29	REAL, INTENT(IN) :: masse0(ngridmx,nlayermx)
30	INTEGER, INTENT(IN) :: lmax(ngridmx)
31
32	! ============================ OUTPUTS ===========================
33
34	REAL, INTENT(OUT) :: dq_therm(ngridmx,nlayermx)
35
36	! ============================ LOCAL =============================
37
38	! REAL detr0(ngridmx,nlayermx)
39	REAL detrd(ngridmx,nlayermx)
40	REAL entrd(ngridmx,nlayermx)
41	REAL fmd(ngridmx,nlayermx+1)
42	REAL q(ngridmx,nlayermx)
43	REAL qa(ngridmx,nlayermx)
44	REAL qd(ngridmx,nlayermx)
45	INTEGER ig,k
46	LOGICAL active(ngridmx,nlayermx)
47	INTEGER lmax_down(ngridmx),lmin_down(ngridmx)
48	INTEGER ncorec
49
50	! =========== Init ==============================================
51
52	entrd(:,:)=0.
53	detrd(:,:)=0.
54	qa(:,:)=q_therm(:,:)
55	q(:,:)=q_therm(:,:)
56	! qd(:,:)=q_therm(:,:)
57	! active(:,:)=.false.
58
59	! previous calculation of zdthl_down uses the divergence of fmd
60	! so it can be negative without problem. Here we include the sign
61	! of fmd in the equations, so it has to be positive
62	!
63	! fmd(:,:)=-fm_down(:,:)
64	!
65	!! ========== Entrainment, Detrainement and Mass =================
66	!
67	!! ========== DOWNDRAFT TRANSPORT DISABLED FOR NOW ===============
68	!
69	! do ig=1,ngridmx
70	! if (ztv(ig,nlayermx)-ztvd(ig,nlayermx) .gt. 0.5) then
71	! print*,"downdraft non nul derniere couche !!! (dqupdown)"
72	! endif
73	! detrd(ig,nlayermx)=0.
74	! entrd(ig,nlayermx)=0.
75	! enddo
76	!
77	! do k=nlayermx-1,1,-1
78	! do ig=1,ngridmx
79	!
80	! if (ztv(ig,k)-ztvd(ig,k) .gt. 0.0001) then
81	! detrd(ig,k)=MAX(0.,(fmd(ig,k+1)*(ztv(ig,k)-ztvd(ig,k+1))) &
82	! & /(ztv(ig,k)-ztvd(ig,k)) - fmd(ig,k))
83	! entrd(ig,k)=MAX(0.,(fmd(ig,k+1)*(ztvd(ig,k)-ztvd(ig,k+1))) &
84	! & /(ztv(ig,k)-ztvd(ig,k)))
85	!
86	! endif
87	! enddo
88	! enddo
89	!
90	!! ======= We have computed entrainment and detrainment from a prescribed
91	!! mass flux and potential temp profile. Due to the way downdraft are parametrized,
92	!! this can yield negative entr and detr. We force it to be positive, but in
93	!! order to conserve tracers, we need to recompute an adequate mass flux
94	!! and modify interface rates, to preserve consistency.
95	! lmax_down(:)=1
96	! lmin_down(:)=1
97	!
98	! do k=1,nlayermx
99	! do ig=1,ngridmx
100	! if ((entrd(ig,k).gt.0.) .or. (detrd(ig,k).gt.0.)) then
101	!! if (entrd(ig,k).gt.detrd(ig,k)) then
102	! lmax_down(ig)=min(k,lmax(ig))
103	! endif
104	! enddo
105	! enddo
106	! do k=nlayermx,1,-1
107	! do ig=1,ngridmx
108	! if ((entrd(ig,k).gt.0.) .or. (detrd(ig,k).gt.0.)) then
109	!! if (detrd(ig,k).gt.entrd(ig,k)) then
110	! lmin_down(ig)=k
111	! endif
112	! enddo
113	! enddo
114	!
115	!
116	! fmd(:,:)=0.
117	!
118	! do ig=1,ngridmx
119	! if ((lmax_down(ig).gt.1) .and. ((lmax_down(ig)-lmin_down(ig)).gt.1)) then
120	!! fmd(ig,lmax_down(ig))=0.
121	!! entrd(ig,lmax_down(ig))=detrd(ig,lmax_down(ig))
122	!! detrd(ig,lmax_down(ig))=0.
123	!! print*,lmin_down(ig),lmax_down(ig),lmax(ig)
124	!
125	! fmd(ig,lmax_down(ig)+1)=0.
126	!
127	! do k=lmax_down(ig),lmin_down(ig)+1,-1
128	! fmd(ig,k)=fmd(ig,k+1)+entrd(ig,k)-detrd(ig,k)
129	! enddo
130	!
131	! fmd(ig,lmin_down(ig))=0.
132	! detrd(ig,lmin_down(ig))=fmd(ig,lmin_down(ig)+1)+entrd(ig,lmin_down(ig))
133	!
134	! else
135	! entrd(ig,:)=0.
136	! detrd(ig,:)=0.
137	! active(ig,:)=.false.
138	! endif
139	!
140	! enddo
141	! ncorec=0
142	! do k=nlayermx,2,-1
143	! do ig=1,ngridmx
144	! if (fmd(ig,k).lt.0.) then
145	!! detrd(ig,k)=max(0.,detrd(ig,k)+fmd(ig,k-1))
146	!! fmd(ig,k-1)=0.
147	!! entrd(ig,k-1)=0.
148	!! detrd(ig,k-1)=0.
149	!! lmin_down(ig)=k-1
150	! fmd(ig,k)=fmd(ig,k+1)
151	! detrd(ig,k)=entrd(ig,k)
152	! ncorec=ncorec+1
153	!! fmd(ig,k)=0.
154	!! detrd(ig,k)=entrd(ig,k)+fmd(ig,k+1)
155	!
156	! endif
157	! enddo
158	! enddo
159	!
160	! if (ncorec .ne. 0) then
161	! print*, 'corrections for negative downward mass flux :',ncorec
162	! endif
163	! print*, lmin_down(:),lmax_down(:)
164	!
165	! do k=2,nlayermx
166	! do ig=1,ngridmx
167	! active(ig,k)=(k.ge.lmin_down(ig)).and.(k.le.lmax_down(ig)) &
168	! & .and.(((fmd(ig,k)+detrd(ig,k))ptimestep).gt.1.e-6masse0(ig,k))
169	! enddo
170	! enddo
171	!
172	!
173	! do ig=1,ngridmx
174	! do k=lmin_down(ig),lmax_down(ig)
175	! if(.not.active(ig,k)) then
176	! active(ig,:)=.false.
177	! endif
178	! enddo
179	! enddo
180	!
181	! if(charvar .eq. 'tke') then
182	! active(:,:)=.false.
183	! endif
184	!
185	!! do ig=1,ngridmx
186	!! active(ig,lmax_down(ig))=(((fmd(ig,lmax_down(ig))+detrd(ig,lmax_down(ig)))* &
187	!! & ptimestep).gt.1.e-8*masse0(ig,lmax_down(ig)))
188	!! enddo
189	!!
190	!! do ig=1,ngridmx
191	!! if (lmax_down(ig).gt.1) then
192	!! do k=lmax_down(ig)-1,lmin_down(ig),-1
193	!! active(ig,k)=(((fmd(ig,k)+detrd(ig,k))* &
194	!! & ptimestep).gt.1.e-8*masse0(ig,k)) &
195	!! & .and. active(ig,k+1)
196	!! enddo
197	!! else
198	!! active(ig,:)=.false.
199	!! endif
200	!! enddo
201	!!
202	!! ========== qa : q in updraft ==================================
203	!
204	do k=2,nlayermx
205	do ig=1,ngridmx
206	if ((fm0(ig,k+1)+detr0(ig,k))*ptimestep.gt. &
207	& 1.e-5*masse0(ig,k)) then
208	qa(ig,k)=(fm0(ig,k)qa(ig,k-1)+entr0(ig,k)q(ig,k)) &
209	& /(fm0(ig,k+1)+detr0(ig,k))
210
211	if ((qa(ig,k).lt.0.) .and. (charvar .ne. 'momentum')) then
212	if ((q(ig,k) .gt. 0.) .and. ( q(ig,k) .gt. 1.e-15 )) then
213	!! only print if it is the thermal scheme which makes qa<0
214	print*,'qa<0 created by thermals!!!',charvar,ig,k,fm0(ig,k),qa(ig,k-1),entr0(ig,k),q(ig,k),fm0(ig,k+1),detr0(ig,k)
215	print*,'---------> Cancelling qa'
216	endif
217	qa(ig,k)=q(ig,k) !! no action of thermals in this case.
218	endif
219	else
220	qa(ig,k)=q(ig,k)
221	endif
222	enddo
223	enddo
224
225	! ========== qd : q in downdraft =================================
226	!
227	! do k=nlayermx-1,1,-1
228	! do ig=1,ngridmx
229	! if (active(ig,k)) then
230	! qd(ig,k)=(fmd(ig,k+1)qd(ig,k+1)+entrd(ig,k)q(ig,k)) &
231	! & /(fmd(ig,k)+detrd(ig,k))
232	!
233	! if ((qd(ig,k).lt.0.) .and. (charvar .ne. 'momentum')) then
234	! print*,'qd<0!!!',charvar,ig,k,fmd(ig,k),qd(ig,k),entrd(ig,k),q(ig,k),fmd(ig,k+1),detrd(ig,k),lmin_down(ig),lmax_down(ig)
235	! print*, '---------> cancelling qd, no downdraft for this gridpoint'
236	! qd(ig,k)=q(ig,k)
237	! active(ig,:)=.false.
238	! endif
239	! else
240	! qd(ig,k)=q(ig,k)
241	! endif
242	!! print*,'active,k,entr,detr,q,qd (down) :',active(ig,k),k,entrd(ig,k),detrd(ig,k),q(ig,k),qd(ig,k)
243	! enddo
244	! enddo
245	!
246	!
247	! ====== dq ======================================================
248
249	do ig=1,ngridmx
250	! if(active(ig,1)) then
251	!
252	! dq_therm(ig,1)=(detr0(ig,1)qa(ig,1)+detrd(ig,1)qd(ig,1) &
253	! & +fm0(ig,2)*q(ig,2) &
254	! & -entr0(ig,1)q(ig,1)-entrd(ig,1)q(ig,1) &
255	! & -fmd(ig,2)*q(ig,1)) &
256	! & *ptimestep/masse0(ig,1)
257	!
258	! else
259	dq_therm(ig,1)=(detr0(ig,1)qa(ig,1)+fm0(ig,2)q(ig,2) &
260	& -entr0(ig,1)*q(ig,1)) &
261	& *ptimestep/masse0(ig,1)
262
263	! endif
264	enddo
265
266	do k=2,nlayermx-1
267	do ig=1, ngridmx
268
269	! if(active(ig,k)) then
270	!
271	! dq_therm(ig,k)=(detr0(ig,k)qa(ig,k)+detrd(ig,k)qd(ig,k) &
272	! & +fm0(ig,k+1)q(ig,k+1)+fmd(ig,k)q(ig,k-1) &
273	! & -entr0(ig,k)q(ig,k)-entrd(ig,k)q(ig,k) &
274	! & -fm0(ig,k)q(ig,k)-fmd(ig,k+1)q(ig,k)) &
275	! & *ptimestep/masse0(ig,k)
276
277
278	! else
279	dq_therm(ig,k)=(detr0(ig,k)qa(ig,k)+fm0(ig,k+1)q(ig,k+1) &
280	& -entr0(ig,k)q(ig,k)-fm0(ig,k)q(ig,k)) &
281	& *ptimestep/masse0(ig,k)
282
283
284	! endif
285
286	enddo
287	enddo
288
289	do ig=1, ngridmx
290
291	! if(active(ig,nlayermx)) then
292	!
293	! dq_therm(ig,nlayermx)=(detr0(ig,nlayermx)qa(ig,nlayermx)+detrd(ig,nlayermx)qd(ig,nlayermx) &
294	! & +fmd(ig,nlayermx)*q(ig,nlayermx-1) &
295	! & -entr0(ig,nlayermx)q(ig,nlayermx)-entrd(ig,nlayermx)q(ig,nlayermx) &
296	! & -fm0(ig,nlayermx)*q(ig,nlayermx)) &
297	! & *ptimestep/masse0(ig,nlayermx)
298
299	! else
300	dq_therm(ig,nlayermx)=(detr0(ig,nlayermx)*qa(ig,nlayermx) &
301	& -entr0(ig,nlayermx)q(ig,nlayermx)-fm0(ig,nlayermx)q(ig,nlayermx)) &
302	& *ptimestep/masse0(ig,nlayermx)
303	! endif
304
305	enddo
306	return
307	end

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

Download in other formats:

Original Format