Context Navigation

thermcell_dq.F90 @ 1910

Last change on this file since 1910 was 1907, checked in by lguez, 11 years ago

Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:

$ svn propget copyright libf/phylmd/physiq.F90
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory

Also added the files defining the CeCILL version 2 license, in French
and English, at the top of the LMDZ tree.

Property copyright set to
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory
Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 8.2 KB

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

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

Download in other formats:

Original Format