Context Navigation

thermcell_dq.F90 @ 5225

Last change on this file since 5225 was 2311, checked in by Ehouarn Millour, 9 years ago

Further modifications to enforce physics/dynamics separation:

moved iniprint.h and misc_mod back to dyn3d_common, as these should only be used by dynamics.
created print_control_mod in the physics to store flags prt_level, lunout, debug to be local to physics (should be used rather than iniprint.h)
created abort_physic.F90 , which does the same job as abort_gcm() did, but should be used instead when in physics.
reactivated inifis (turned it into a module, inifis_mod.F90) to initialize physical constants and print_control_mod flags.

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.6 KB

Rev	Line
[1738]	1	subroutine thermcell_dq(ngrid,nlay,impl,ptimestep,fm,entr, &
[878]	2	& masse,q,dq,qa,lev_out)
[2311]	3	USE print_control_mod, ONLY: prt_level
[878]	4	implicit none
	5
	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
[1985]	42	if (impl<=-1) then
	43	call thermcell_dq_o(ngrid,nlay,impl,ptimestep,fm,entr, &
[1738]	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
[1943]	55	print,'entrdt>m,1',k,entr(ig,k)*ptimestep,masse(ig,k)
	56	abort_message = 'entr dt > m, 1st'
[2311]	57	CALL abort_physic (modname,abort_message,1)
[1738]	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
[1985]	126	! FH debut de modif : le calcul ci dessous modifiait numériquement
	127	! la concentration quand le flux de masse etait nul car on divisait
	128	! puis multipliait par masse/ptimestep.
	129	! q(:,k)=(masse(:,k)q(:,k)/ptimestep+fqa(:,k)-fqa(:,k+1)+fm(:,k+1)q(:,k+1)) &
	130	! & /(fm(:,k)+masse(:,k)/ptimestep)
	131	q(:,k)=(q(:,k)+ptimestep/masse(:,k)(fqa(:,k)-fqa(:,k+1)+fm(:,k+1)q(:,k+1))) &
	132	& /(1.+fm(:,k)*ptimestep/masse(:,k))
	133	! FH fin de modif.
[1738]	134	enddo
	135	endif
	136
	137	! Tendencies
	138	do k=1,nlay
	139	do ig=1,ngrid
	140	dq(ig,k)=(q(ig,k)-qold(ig,k))/ptimestep
	141	q(ig,k)=qold(ig,k)
	142	enddo
	143	enddo
	144
	145	return
	146	end
	147
	148	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	149	! Obsolete version kept for convergence with Cmip5 NPv3.1 simulations
	150	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	151
[1985]	152	subroutine thermcell_dq_o(ngrid,nlay,impl,ptimestep,fm,entr, &
[1738]	153	& masse,q,dq,qa,lev_out)
[2311]	154	USE print_control_mod, ONLY: prt_level
[1738]	155	implicit none
	156
	157	!=======================================================================
	158	!
	159	! Calcul du transport verticale dans la couche limite en presence
	160	! de "thermiques" explicitement representes
	161	! calcul du dq/dt une fois qu'on connait les ascendances
	162	!
	163	!=======================================================================
	164
[1985]	165	integer ngrid,nlay,impl
[878]	166
	167	real ptimestep
	168	real masse(ngrid,nlay),fm(ngrid,nlay+1)
	169	real entr(ngrid,nlay)
	170	real q(ngrid,nlay)
	171	real dq(ngrid,nlay)
	172	integer lev_out ! niveau pour les print
	173
	174	real qa(ngrid,nlay),detr(ngrid,nlay),wqd(ngrid,nlay+1)
	175
[972]	176	real zzm
	177
[878]	178	integer ig,k
[972]	179	real cfl
[878]	180
[972]	181	real qold(ngrid,nlay)
	182	real ztimestep
	183	integer niter,iter
[1403]	184	CHARACTER (LEN=20) :: modname='thermcell_dq'
	185	CHARACTER (LEN=80) :: abort_message
[878]	186
[972]	187
	188
	189	! Calcul du critere CFL pour l'advection dans la subsidence
[983]	190	cfl = 0.
[972]	191	do k=1,nlay
	192	do ig=1,ngrid
	193	zzm=masse(ig,k)/ptimestep
	194	cfl=max(cfl,fm(ig,k)/zzm)
	195	if (entr(ig,k).gt.zzm) then
[1943]	196	print,'entrdt>m,2',k,entr(ig,k)*ptimestep,masse(ig,k)
	197	abort_message = 'entr dt > m, 2nd'
[2311]	198	CALL abort_physic (modname,abort_message,1)
[972]	199	endif
	200	enddo
	201	enddo
	202
	203	!IM 090508 print*,'CFL CFL CFL CFL ',cfl
	204
	205	#undef CFL
	206	#ifdef CFL
	207	! On subdivise le calcul en niter pas de temps.
	208	niter=int(cfl)+1
	209	#else
	210	niter=1
	211	#endif
	212
	213	ztimestep=ptimestep/niter
	214	qold=q
	215
	216
	217	do iter=1,niter
[938]	218	if (prt_level.ge.1) print*,'Q2 THERMCEL_DQ 0'
[878]	219
[972]	220	! calcul du detrainement
[878]	221	do k=1,nlay
	222	do ig=1,ngrid
	223	detr(ig,k)=fm(ig,k)-fm(ig,k+1)+entr(ig,k)
	224	! print*,'Q2 DQ ',detr(ig,k),fm(ig,k),entr(ig,k)
	225	!test
	226	if (detr(ig,k).lt.0.) then
	227	entr(ig,k)=entr(ig,k)-detr(ig,k)
	228	detr(ig,k)=0.
	229	! print*,'detr2<0!!!','ig=',ig,'k=',k,'f=',fm(ig,k),
	230	! s 'f+1=',fm(ig,k+1),'e=',entr(ig,k),'d=',detr(ig,k)
	231	endif
	232	if (fm(ig,k+1).lt.0.) then
	233	! print*,'fm2<0!!!'
	234	endif
	235	if (entr(ig,k).lt.0.) then
	236	! print*,'entr2<0!!!'
	237	endif
	238	enddo
	239	enddo
	240
	241	! calcul de la valeur dans les ascendances
	242	do ig=1,ngrid
	243	qa(ig,1)=q(ig,1)
	244	enddo
	245
	246	do k=2,nlay
	247	do ig=1,ngrid
[972]	248	if ((fm(ig,k+1)+detr(ig,k))*ztimestep.gt. &
[878]	249	& 1.e-5*masse(ig,k)) then
	250	qa(ig,k)=(fm(ig,k)qa(ig,k-1)+entr(ig,k)q(ig,k)) &
	251	& /(fm(ig,k+1)+detr(ig,k))
	252	else
	253	qa(ig,k)=q(ig,k)
	254	endif
	255	if (qa(ig,k).lt.0.) then
	256	! print*,'qa<0!!!'
	257	endif
	258	if (q(ig,k).lt.0.) then
	259	! print*,'q<0!!!'
	260	endif
	261	enddo
	262	enddo
	263
[972]	264	! Calcul du flux subsident
	265
[878]	266	do k=2,nlay
	267	do ig=1,ngrid
[972]	268	#undef centre
	269	#ifdef centre
	270	wqd(ig,k)=fm(ig,k)0.5(q(ig,k-1)+q(ig,k))
	271	#else
	272
	273	#define plusqueun
	274	#ifdef plusqueun
	275	! Schema avec advection sur plus qu'une maille.
	276	zzm=masse(ig,k)/ztimestep
	277	if (fm(ig,k)>zzm) then
	278	wqd(ig,k)=zzmq(ig,k)+(fm(ig,k)-zzm)q(ig,k+1)
	279	else
	280	wqd(ig,k)=fm(ig,k)*q(ig,k)
	281	endif
	282	#else
[878]	283	wqd(ig,k)=fm(ig,k)*q(ig,k)
[972]	284	#endif
	285	#endif
	286
[878]	287	if (wqd(ig,k).lt.0.) then
	288	! print*,'wqd<0!!!'
	289	endif
	290	enddo
	291	enddo
	292	do ig=1,ngrid
	293	wqd(ig,1)=0.
	294	wqd(ig,nlay+1)=0.
	295	enddo
	296
[972]	297
	298	! Calcul des tendances
[878]	299	do k=1,nlay
	300	do ig=1,ngrid
[972]	301	q(ig,k)=q(ig,k)+(detr(ig,k)qa(ig,k)-entr(ig,k)q(ig,k) &
[878]	302	& -wqd(ig,k)+wqd(ig,k+1)) &
[972]	303	& *ztimestep/masse(ig,k)
[878]	304	! if (dq(ig,k).lt.0.) then
	305	! print*,'dq<0!!!'
	306	! endif
	307	enddo
	308	enddo
	309
[972]	310
	311	enddo
	312
	313
	314	! Calcul des tendances
	315	do k=1,nlay
	316	do ig=1,ngrid
	317	dq(ig,k)=(q(ig,k)-qold(ig,k))/ptimestep
	318	q(ig,k)=qold(ig,k)
	319	enddo
	320	enddo
	321
[878]	322	return
	323	end

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

Download in other formats:

Original Format