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

Rev	Line
[185]	1	subroutine thermcell_dqupdown(ngrid,nlayer,ptimestep,fm0,entr0, &
[161]	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
[185]	15	#include "dimensions.h"
	16	#include "dimphys.h"
	17
[161]	18	! ============================ INPUTS ============================
	19
[185]	20	INTEGER, INTENT(IN) :: ngrid,nlayer
[161]	21	REAL, INTENT(IN) :: ptimestep
[185]	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)
[161]	28	CHARACTER (LEN=20), INTENT(IN) :: charvar
[185]	29	REAL, INTENT(IN) :: masse0(ngridmx,nlayermx)
	30	INTEGER, INTENT(IN) :: lmax(ngridmx)
[161]	31
	32	! ============================ OUTPUTS ===========================
	33
[185]	34	REAL, INTENT(OUT) :: dq_therm(ngridmx,nlayermx)
[161]	35
	36	! ============================ LOCAL =============================
	37
[185]	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)
[161]	45	INTEGER ig,k
[185]	46	LOGICAL active(ngridmx,nlayermx)
	47	INTEGER lmax_down(ngridmx),lmin_down(ngridmx)
[161]	48	INTEGER ncorec
	49
	50	! =========== Init ==============================================
	51
	52	entrd(:,:)=0.
	53	detrd(:,:)=0.
	54	qa(:,:)=q_therm(:,:)
	55	q(:,:)=q_therm(:,:)
[313]	56	! qd(:,:)=q_therm(:,:)
	57	! active(:,:)=.false.
[161]	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	!
[313]	63	! fmd(:,:)=-fm_down(:,:)
[161]	64	!
	65	!! ========== Entrainment, Detrainement and Mass =================
	66	!
	67	!! ========== DOWNDRAFT TRANSPORT DISABLED FOR NOW ===============
	68	!
[185]	69	! do ig=1,ngridmx
	70	! if (ztv(ig,nlayermx)-ztvd(ig,nlayermx) .gt. 0.5) then
[161]	71	! print*,"downdraft non nul derniere couche !!! (dqupdown)"
	72	! endif
[185]	73	! detrd(ig,nlayermx)=0.
	74	! entrd(ig,nlayermx)=0.
[161]	75	! enddo
	76	!
[185]	77	! do k=nlayermx-1,1,-1
	78	! do ig=1,ngridmx
[161]	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	!
[185]	98	! do k=1,nlayermx
	99	! do ig=1,ngridmx
[161]	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
[185]	106	! do k=nlayermx,1,-1
	107	! do ig=1,ngridmx
[161]	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	!
[185]	118	! do ig=1,ngridmx
[161]	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
[185]	142	! do k=nlayermx,2,-1
	143	! do ig=1,ngridmx
[161]	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	!
[185]	165	! do k=2,nlayermx
	166	! do ig=1,ngridmx
[161]	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	!
[185]	173	! do ig=1,ngridmx
[161]	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]	185	!! do ig=1,ngridmx
[161]	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	!!
[185]	190	!! do ig=1,ngridmx
[161]	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	!
[185]	204	do k=2,nlayermx
	205	do ig=1,ngridmx
[161]	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
[278]	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.
[161]	218	endif
	219	else
	220	qa(ig,k)=q(ig,k)
	221	endif
	222	enddo
	223	enddo
	224
	225	! ========== qd : q in downdraft =================================
	226	!
[185]	227	! do k=nlayermx-1,1,-1
	228	! do ig=1,ngridmx
[161]	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
[185]	249	do ig=1,ngridmx
[313]	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
[161]	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
[313]	263	! endif
[161]	264	enddo
	265
[185]	266	do k=2,nlayermx-1
	267	do ig=1, ngridmx
[161]	268
[313]	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)
[161]	276
	277
[313]	278	! else
[161]	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
[313]	284	! endif
[161]	285
	286	enddo
	287	enddo
	288
[185]	289	do ig=1, ngridmx
[161]	290
[313]	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)
[161]	298
[313]	299	! else
[185]	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)
[313]	303	! endif
[161]	304
	305	enddo
	306	return
	307	end

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

Download in other formats:

Original Format