Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

lmdz_thermcell_dq.F90 @ 5065

Last change on this file since 5065 was 4590, checked in by fhourdin, 17 months ago
Passage des thermiques a la nouvelle norme.
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.8 KB

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

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

Download in other formats:

Original Format