Context Navigation

thermcell_dq.F90 @ 630

Last change on this file since 630 was 142, checked in by lfita, 10 years ago

Incorporing all the work done in LMDZ_WRFmeas about:

1.- Chekcing for NaNs?
2.- Instabilities in thermcell (negative SQRT & robustness IF)

File size: 8.2 KB

Rev	Line
[1]	1	subroutine thermcell_dq(ngrid,nlay,impl,ptimestep,fm,entr, &
	2	& masse,q,dq,qa,lev_out)
	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	!
	12	! Modif 2013/01/04 (FH hourdin@lmd.jussieu.fr)
	13	! Introduction of an implicit computation of vertical advection in
	14	! the environment of thermal plumes in thermcell_dq
	15	! impl = 0 : explicit, 1 : implicit, -1 : old version
	16	!
	17	!=======================================================================
	18
	19	integer ngrid,nlay,impl
	20
	21	real ptimestep
	22	real masse(ngrid,nlay),fm(ngrid,nlay+1)
	23	real entr(ngrid,nlay)
	24	real q(ngrid,nlay)
	25	real dq(ngrid,nlay)
	26	integer lev_out ! niveau pour les print
	27
	28	real qa(ngrid,nlay),detr(ngrid,nlay),wqd(ngrid,nlay+1)
	29
	30	real zzm
	31
	32	integer ig,k
	33	real cfl
	34
	35	real qold(ngrid,nlay),fqa(ngrid,nlay+1)
	36	integer niter,iter
	37	CHARACTER (LEN=20) :: modname='thermcell_dq'
	38	CHARACTER (LEN=80) :: abort_message
	39
	40	! Old explicite scheme
	41	if (impl==-1) then
	42	call thermcell_dq_o(ngrid,nlay,ptimestep,fm,entr, &
	43	& masse,q,dq,qa,lev_out)
	44	return
	45	endif
	46
	47	! Calcul du critere CFL pour l'advection dans la subsidence
	48	cfl = 0.
	49	do k=1,nlay
	50	do ig=1,ngrid
	51	zzm=masse(ig,k)/ptimestep
	52	cfl=max(cfl,fm(ig,k)/zzm)
	53	if (entr(ig,k).gt.zzm) then
	54	print,'entr dt > m ',entr(ig,k)ptimestep,masse(ig,k)
	55	abort_message = ''
	56	CALL abort_gcm (modname,abort_message,1)
	57	endif
	58	enddo
	59	enddo
	60
	61	qold=q
	62
	63
	64	if (prt_level.ge.1) print*,'Q2 THERMCEL_DQ 0'
	65
	66	! calcul du detrainement
	67	do k=1,nlay
	68	do ig=1,ngrid
	69	detr(ig,k)=fm(ig,k)-fm(ig,k+1)+entr(ig,k)
	70	! print*,'Q2 DQ ',detr(ig,k),fm(ig,k),entr(ig,k)
	71	!test
	72	if (detr(ig,k).lt.0.) then
	73	entr(ig,k)=entr(ig,k)-detr(ig,k)
	74	detr(ig,k)=0.
	75	! print*,'detr2<0!!!','ig=',ig,'k=',k,'f=',fm(ig,k),
	76	! s 'f+1=',fm(ig,k+1),'e=',entr(ig,k),'d=',detr(ig,k)
	77	endif
	78	if (fm(ig,k+1).lt.0.) then
	79	! print*,'fm2<0!!!'
	80	endif
	81	if (entr(ig,k).lt.0.) then
	82	! print*,'entr2<0!!!'
	83	endif
	84	enddo
	85	enddo
	86
	87	! Computation of tracer concentrations in the ascending plume
	88	do ig=1,ngrid
	89	qa(ig,1)=q(ig,1)
	90	enddo
	91
	92	do k=2,nlay
	93	do ig=1,ngrid
	94	if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt. &
	95	& 1.e-5*masse(ig,k)) then
	96	qa(ig,k)=(fm(ig,k)qa(ig,k-1)+entr(ig,k)q(ig,k)) &
	97	& /(fm(ig,k+1)+detr(ig,k))
	98	else
	99	qa(ig,k)=q(ig,k)
	100	endif
	101	if (qa(ig,k).lt.0.) then
	102	! print*,'qa<0!!!'
	103	endif
	104	if (q(ig,k).lt.0.) then
	105	! print*,'q<0!!!'
	106	endif
	107	enddo
	108	enddo
	109
	110	! Plume vertical flux
	111	do k=2,nlay-1
	112	fqa(:,k)=fm(:,k)*qa(:,k-1)
	113	enddo
	114	fqa(:,1)=0. ; fqa(:,nlay)=0.
	115
	116	! Trace species evolution
	117	if (impl==0) then
	118	do k=1,nlay-1
	119	q(:,k)=q(:,k)+(fqa(:,k)-fqa(:,k+1)-fm(:,k)q(:,k)+fm(:,k+1)q(:,k+1)) &
	120	& *ptimestep/masse(:,k)
	121	enddo
	122	else
	123	do k=nlay-1,1,-1
	124	q(:,k)=(masse(:,k)q(:,k)/ptimestep+fqa(:,k)-fqa(:,k+1)+fm(:,k+1)q(:,k+1)) &
	125	& /(fm(:,k)+masse(:,k)/ptimestep)
	126	enddo
	127	endif
	128
	129	! Tendencies
	130	do k=1,nlay
	131	do ig=1,ngrid
	132	dq(ig,k)=(q(ig,k)-qold(ig,k))/ptimestep
	133	q(ig,k)=qold(ig,k)
	134	enddo
	135	enddo
	136
	137	return
	138	end
	139
	140	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	141	! Obsolete version kept for convergence with Cmip5 NPv3.1 simulations
	142	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	143
	144	subroutine thermcell_dq_o(ngrid,nlay,ptimestep,fm,entr, &
	145	& masse,q,dq,qa,lev_out)
	146	implicit none
	147
	148	#include "iniprint.h"
	149	!=======================================================================
	150	!
	151	! Calcul du transport verticale dans la couche limite en presence
	152	! de "thermiques" explicitement representes
	153	! calcul du dq/dt une fois qu'on connait les ascendances
	154	!
	155	!=======================================================================
	156
	157	integer ngrid,nlay
	158
	159	real ptimestep
	160	real masse(ngrid,nlay),fm(ngrid,nlay+1)
	161	real entr(ngrid,nlay)
	162	real q(ngrid,nlay)
	163	real dq(ngrid,nlay)
	164	integer lev_out ! niveau pour les print
	165
	166	real qa(ngrid,nlay),detr(ngrid,nlay),wqd(ngrid,nlay+1)
	167
	168	real zzm
	169
	170	integer ig,k
	171	real cfl
	172
	173	real qold(ngrid,nlay)
	174	real ztimestep
	175	integer niter,iter
	176	CHARACTER (LEN=20) :: modname='thermcell_dq'
	177	CHARACTER (LEN=80) :: abort_message
	178
	179
	180
	181	! Calcul du critere CFL pour l'advection dans la subsidence
	182	cfl = 0.
	183	do k=1,nlay
	184	do ig=1,ngrid
	185	zzm=masse(ig,k)/ptimestep
	186	cfl=max(cfl,fm(ig,k)/zzm)
	187	if (entr(ig,k).gt.zzm) then
	188	print,'entr dt > m ',entr(ig,k)ptimestep,masse(ig,k)
	189	abort_message = ''
	190	CALL abort_gcm (modname,abort_message,1)
	191	endif
	192	enddo
	193	enddo
	194
	195	!IM 090508 print*,'CFL CFL CFL CFL ',cfl
	196
	197	#undef CFL
	198	#ifdef CFL
	199	! On subdivise le calcul en niter pas de temps.
	200	niter=int(cfl)+1
	201	#else
	202	niter=1
	203	#endif
	204
	205	ztimestep=ptimestep/niter
	206	qold=q
	207
	208
	209	do iter=1,niter
	210	if (prt_level.ge.1) print*,'Q2 THERMCEL_DQ 0'
	211
	212	! calcul du detrainement
	213	do k=1,nlay
	214	do ig=1,ngrid
	215	detr(ig,k)=fm(ig,k)-fm(ig,k+1)+entr(ig,k)
	216	! print*,'Q2 DQ ',detr(ig,k),fm(ig,k),entr(ig,k)
	217	!test
	218	if (detr(ig,k).lt.0.) then
	219	entr(ig,k)=entr(ig,k)-detr(ig,k)
	220	detr(ig,k)=0.
	221	! print*,'detr2<0!!!','ig=',ig,'k=',k,'f=',fm(ig,k),
	222	! s 'f+1=',fm(ig,k+1),'e=',entr(ig,k),'d=',detr(ig,k)
	223	endif
	224	if (fm(ig,k+1).lt.0.) then
	225	! print*,'fm2<0!!!'
	226	endif
	227	if (entr(ig,k).lt.0.) then
	228	! print*,'entr2<0!!!'
	229	endif
	230	enddo
	231	enddo
	232
	233	! calcul de la valeur dans les ascendances
	234	do ig=1,ngrid
	235	qa(ig,1)=q(ig,1)
	236	enddo
	237
	238	do k=2,nlay
	239	do ig=1,ngrid
	240	if ((fm(ig,k+1)+detr(ig,k))*ztimestep.gt. &
	241	& 1.e-5*masse(ig,k)) then
	242	qa(ig,k)=(fm(ig,k)qa(ig,k-1)+entr(ig,k)q(ig,k)) &
	243	& /(fm(ig,k+1)+detr(ig,k))
	244	else
	245	qa(ig,k)=q(ig,k)
	246	endif
	247	if (qa(ig,k).lt.0.) then
	248	! print*,'qa<0!!!'
	249	endif
	250	if (q(ig,k).lt.0.) then
	251	! print*,'q<0!!!'
	252	endif
	253	enddo
	254	enddo
	255
	256	! Calcul du flux subsident
	257
	258	do k=2,nlay
	259	do ig=1,ngrid
	260	#undef centre
	261	#ifdef centre
	262	wqd(ig,k)=fm(ig,k)0.5(q(ig,k-1)+q(ig,k))
	263	#else
	264
	265	#define plusqueun
	266	#ifdef plusqueun
	267	! Schema avec advection sur plus qu'une maille.
	268	zzm=masse(ig,k)/ztimestep
	269	if (fm(ig,k)>zzm) then
	270	wqd(ig,k)=zzmq(ig,k)+(fm(ig,k)-zzm)q(ig,k+1)
	271	else
	272	wqd(ig,k)=fm(ig,k)*q(ig,k)
	273	endif
	274	#else
	275	wqd(ig,k)=fm(ig,k)*q(ig,k)
	276	#endif
	277	#endif
	278
	279	if (wqd(ig,k).lt.0.) then
	280	! print*,'wqd<0!!!'
	281	endif
	282	enddo
	283	enddo
	284	do ig=1,ngrid
	285	wqd(ig,1)=0.
	286	wqd(ig,nlay+1)=0.
	287	enddo
	288
	289
	290	! Calcul des tendances
	291	do k=1,nlay
	292	do ig=1,ngrid
	293	q(ig,k)=q(ig,k)+(detr(ig,k)qa(ig,k)-entr(ig,k)q(ig,k) &
	294	& -wqd(ig,k)+wqd(ig,k+1)) &
	295	& *ztimestep/masse(ig,k)
	296	! if (dq(ig,k).lt.0.) then
	297	! print*,'dq<0!!!'
	298	! endif
	299	enddo
	300	enddo
	301
	302
	303	enddo
	304
	305
	306	! Calcul des tendances
	307	do k=1,nlay
	308	do ig=1,ngrid
	309	dq(ig,k)=(q(ig,k)-qold(ig,k))/ptimestep
	310	q(ig,k)=qold(ig,k)
	311	enddo
	312	enddo
	313
	314	return
	315	end

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

Download in other formats:

Original Format