Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

thermcell_dq.F90 @ 5507

Last change on this file since 5507 was 3331, checked in by acozic, 7 years ago
Add modification for isotopes
Property svn:executable set to ``*
File size: 22.9 KB

Rev	Line
[3331]	1	subroutine thermcell_dq(ngrid,nlay,impl,ptimestep,fm,entr, &
	2	& masse,q,dq,qa,lev_out)
	3	USE print_control_mod, ONLY: prt_level
	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	!
	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
	41	! Old explicite scheme
	42	if (impl<=-1) then
	43	call thermcell_dq_o(ngrid,nlay,impl,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,'entrdt>m,1',k,entr(ig,k)*ptimestep,masse(ig,k)
	56	abort_message = 'entr dt > m, 1st'
	57	CALL abort_physic (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	! 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.
	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	#ifdef ISO
	149	subroutine thermcell_dq_iso(ngrid,nlay,impl,ptimestep,fm,entr, &
	150	& masse,q,dq,qa,lev_out,xt,dxt,xta)
	151	USE infotrac_phy, ONLY : ntraciso
	152	#ifdef ISOVERIF
	153	USE isotopes_mod, ONLY: iso_eau,iso_HDO
	154	USE isotopes_verif_mod, ONLY: iso_verif_egalite, &
	155	iso_verif_aberrant_enc_vect2D,iso_verif_egalite_vect2D
	156	#endif
	157	implicit none
	158
	159
	160	#include "iniprint.h"
	161	!=======================================================================
	162	!
	163	! Calcul du transport verticale dans la couche limite en presence
	164	! de "thermiques" explicitement representes
	165	! calcul du dq/dt une fois qu'on connait les ascendances
	166	!
	167	! Modif 2013/01/04 (FH hourdin@lmd.jussieu.fr)
	168	! Introduction of an implicit computation of vertical advection in
	169	! the environment of thermal plumes in thermcell_dq
	170	! impl = 0 : explicit, 1 : implicit, -1 : old version
	171	!
	172	!=======================================================================
	173
	174	integer ngrid,nlay,impl
	175
	176	real ptimestep
	177	real masse(ngrid,nlay),fm(ngrid,nlay+1)
	178	real entr(ngrid,nlay)
	179	real q(ngrid,nlay)
	180	real dq(ngrid,nlay)
	181	integer lev_out ! niveau pour les print
	182
	183	real qa(ngrid,nlay),detr(ngrid,nlay),wqd(ngrid,nlay+1)
	184
	185	real zzm
	186
	187	integer ig,k
	188	real cfl
	189
	190	real qold(ngrid,nlay),fqa(ngrid,nlay+1)
	191	integer niter,iter
	192	CHARACTER (LEN=20) :: modname='thermcell_dq'
	193	CHARACTER (LEN=80) :: abort_message
	194
	195	! ifdef ISO
	196	real xt(ntraciso,ngrid,nlay)
	197	real dxt(ntraciso,ngrid,nlay)
	198	real xta(ntraciso,ngrid,nlay)
	199	real q_avant_entr(ngrid,nlay)
	200	integer ixt
	201	real wxtd(ntraciso,ngrid,nlay+1)
	202	real xtold(ntraciso,ngrid,nlay)
	203	real fxta(ntraciso,ngrid,nlay+1)
	204	! endif
	205
	206	! Old explicite scheme
	207	if (impl==-1) then
	208	call thermcell_dq_o_iso(ngrid,nlay,ptimestep,fm,entr, &
	209	& masse,q,dq,qa,lev_out,xt,dxt,xta)
	210	return
	211	endif
	212
	213	! Calcul du critere CFL pour l'advection dans la subsidence
	214	cfl = 0.
	215	do k=1,nlay
	216	do ig=1,ngrid
	217	zzm=masse(ig,k)/ptimestep
	218	cfl=max(cfl,fm(ig,k)/zzm)
	219	if (entr(ig,k).gt.zzm) then
	220	print,'entr dt > m ',entr(ig,k)ptimestep,masse(ig,k)
	221	abort_message = ''
	222	CALL abort_gcm (modname,abort_message,1)
	223	endif
	224	enddo
	225	enddo
	226
	227	qold=q
	228	!#ifdef ISO
	229	xtold=xt
	230	!#endif
	231
	232	if (prt_level.ge.1) print*,'Q2 THERMCEL_DQ 0'
	233
	234	! calcul du detrainement
	235	do k=1,nlay
	236	do ig=1,ngrid
	237	detr(ig,k)=fm(ig,k)-fm(ig,k+1)+entr(ig,k)
	238	! print*,'Q2 DQ ',detr(ig,k),fm(ig,k),entr(ig,k)
	239	!test
	240	if (detr(ig,k).lt.0.) then
	241	entr(ig,k)=entr(ig,k)-detr(ig,k)
	242	detr(ig,k)=0.
	243	! print*,'detr2<0!!!','ig=',ig,'k=',k,'f=',fm(ig,k),
	244	! s 'f+1=',fm(ig,k+1),'e=',entr(ig,k),'d=',detr(ig,k)
	245	endif
	246	if (fm(ig,k+1).lt.0.) then
	247	! print*,'fm2<0!!!'
	248	endif
	249	if (entr(ig,k).lt.0.) then
	250	! print*,'entr2<0!!!'
	251	endif
	252	enddo
	253	enddo
	254	! bla90
	255
	256	! Computation of tracer concentrations in the ascending plume
	257	do ig=1,ngrid
	258	qa(ig,1)=q(ig,1)
	259	enddo
	260	!#ifdef ISO
	261	do ig=1,ngrid
	262	do ixt=1,ntraciso
	263	xta(ixt,ig,1)=xt(ixt,ig,1)
	264	enddo
	265	enddo
	266	!#endif
	267
	268	do k=2,nlay
	269	do ig=1,ngrid
	270	if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt. &
	271	& 1.e-5*masse(ig,k)) then
	272	qa(ig,k)=(fm(ig,k)qa(ig,k-1)+entr(ig,k)q(ig,k)) &
	273	& /(fm(ig,k+1)+detr(ig,k))
	274	!#ifdef ISO
	275	do ixt=1,ntraciso
	276	xta(ixt,ig,k)=(fm(ig,k)*xta(ixt,ig,k-1) &
	277	& +entr(ig,k)*xt(ixt,ig,k)) &
	278	& /(fm(ig,k+1)+detr(ig,k))
	279	enddo
	280	!#endif
	281	else
	282	qa(ig,k)=q(ig,k)
	283	!#ifdef ISO
	284	do ixt=1,ntraciso
	285	xta(ixt,ig,k)=xt(ixt,ig,k)
	286	enddo
	287	!#endif
	288	endif
	289	if (qa(ig,k).lt.0.) then
	290	! print*,'qa<0!!!'
	291	endif
	292	if (q(ig,k).lt.0.) then
	293	! print*,'q<0!!!'
	294	endif
	295	enddo
	296	enddo
	297
	298	!#ifdef ISO
	299	#ifdef ISOVERIF
	300	if (iso_HDO.gt.0) then
	301	call iso_verif_aberrant_enc_vect2D( &
	302	& xta,qa, &
	303	& 'thermcell_dq_iso 320, qa',ntraciso,ngrid,nlay)
	304	endif
	305	#endif
	306	!#endif
	307	!
	308	! Plume vertical flux
	309	do k=2,nlay-1
	310	fqa(:,k)=fm(:,k)*qa(:,k-1)
	311	enddo
	312	fqa(:,1)=0. ; fqa(:,nlay)=0.
	313	!#ifdef ISO
	314	do k=2,nlay-1
	315	do ixt=1,ntraciso
	316	fxta(ixt,:,k)=fm(:,k)*xta(ixt,:,k-1)
	317	enddo
	318	enddo
	319	fxta(ixt,:,1)=0. ; fxta(ixt,:,nlay)=0.
	320	!#endif
	321
	322
	323	! Trace species evolution
	324	if (impl==0) then
	325	do k=1,nlay-1
	326	q(:,k)=q(:,k)+(fqa(:,k)-fqa(:,k+1)-fm(:,k)q(:,k)+fm(:,k+1)q(:,k+1)) &
	327	& *ptimestep/masse(:,k)
	328	!#ifdef ISO
	329	do ixt=1,ntraciso
	330	xt(ixt,:,k)=xt(ixt,:,k)+(fxta(ixt,:,k)-fxta(ixt,:,k+1) &
	331	& -fm(:,k)xt(ixt,:,k)+fm(:,k+1)xt(ixt,:,k+1)) &
	332	& *ptimestep/masse(:,k)
	333	enddo ! do ixt=1,ntraciso
	334	!#endif
	335	enddo !do k=1,nlay-1
	336	else !if (impl==0) then
	337	do k=nlay-1,1,-1
	338	q(:,k)=(masse(:,k)q(:,k)/ptimestep+fqa(:,k)-fqa(:,k+1)+fm(:,k+1)q(:,k+1)) &
	339	& /(fm(:,k)+masse(:,k)/ptimestep)
	340	!#ifdef ISO
	341	do ixt=1,ntraciso
	342	xt(ixt,:,k)=xt(ixt,:,k)+(fxta(ixt,:,k)-fxta(ixt,:,k+1) &
	343	& -fm(:,k)xt(ixt,:,k)+fm(:,k+1)xt(ixt,:,k+1)) &
	344	& *ptimestep/masse(:,k)
	345	enddo ! do ixt=1,ntraciso
	346	!#endif
	347	enddo !do k=nlay-1,1,-1
	348	endif !if (impl==0) then
	349
	350	! Tendencies
	351	do k=1,nlay
	352	do ig=1,ngrid
	353	dq(ig,k)=(q(ig,k)-qold(ig,k))/ptimestep
	354	q(ig,k)=qold(ig,k)
	355	!#ifdef ISO
	356	do ixt=1,ntraciso
	357	dxt(ixt,ig,k)=(xt(ixt,ig,k)-xtold(ixt,ig,k))/ptimestep
	358	xt(ixt,ig,k)=xtold(ixt,ig,k)
	359	enddo ! do ixt=1,ntraciso
	360	!#endif
	361	enddo
	362	enddo
	363
	364	#ifdef ISOVERIF
	365	if (iso_HDO.gt.0) then
	366	call iso_verif_aberrant_enc_vect2D( &
	367	& xt,q, &
	368	& 'thermcell_dq_iso 219, q',ntraciso,ngrid,nlay)
	369	endif
	370	if (iso_eau.gt.0) then
	371	call iso_verif_egalite_vect2D(xt,q, &
	372	& 'thermcell_dq_iso 223, q',ntraciso,ngrid,nlay)
	373	call iso_verif_egalite_vect2D(dxt,dq, &
	374	& 'thermcell_dq_iso 224, dq',ntraciso,ngrid,nlay)
	375	endif
	376	#endif
	377
	378	return
	379	end
	380
	381	#endif
	382	! ifdef iso
	383	!
	384	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	385	! Obsolete version kept for convergence with Cmip5 NPv3.1 simulations
	386	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	387
	388	subroutine thermcell_dq_o(ngrid,nlay,impl,ptimestep,fm,entr, &
	389	& masse,q,dq,qa,lev_out)
	390	USE print_control_mod, ONLY: prt_level
	391	implicit none
	392
	393	!=======================================================================
	394	!
	395	! Calcul du transport verticale dans la couche limite en presence
	396	! de "thermiques" explicitement representes
	397	! calcul du dq/dt une fois qu'on connait les ascendances
	398	!
	399	!=======================================================================
	400
	401	integer ngrid,nlay,impl
	402
	403	real ptimestep
	404	real masse(ngrid,nlay),fm(ngrid,nlay+1)
	405	real entr(ngrid,nlay)
	406	real q(ngrid,nlay)
	407	real dq(ngrid,nlay)
	408	integer lev_out ! niveau pour les print
	409
	410	real qa(ngrid,nlay),detr(ngrid,nlay),wqd(ngrid,nlay+1)
	411
	412	real zzm
	413
	414	integer ig,k
	415	real cfl
	416
	417	real qold(ngrid,nlay)
	418	real ztimestep
	419	integer niter,iter
	420	CHARACTER (LEN=20) :: modname='thermcell_dq'
	421	CHARACTER (LEN=80) :: abort_message
	422
	423
	424
	425	! Calcul du critere CFL pour l'advection dans la subsidence
	426	cfl = 0.
	427	do k=1,nlay
	428	do ig=1,ngrid
	429	zzm=masse(ig,k)/ptimestep
	430	cfl=max(cfl,fm(ig,k)/zzm)
	431	if (entr(ig,k).gt.zzm) then
	432	print,'entrdt>m,2',k,entr(ig,k)*ptimestep,masse(ig,k)
	433	abort_message = 'entr dt > m, 2nd'
	434	CALL abort_physic (modname,abort_message,1)
	435	endif
	436	enddo
	437	enddo
	438
	439	!IM 090508 print*,'CFL CFL CFL CFL ',cfl
	440
	441	#undef CFL
	442	#ifdef CFL
	443	! On subdivise le calcul en niter pas de temps.
	444	niter=int(cfl)+1
	445	#else
	446	niter=1
	447	#endif
	448
	449	ztimestep=ptimestep/niter
	450	qold=q
	451
	452
	453	do iter=1,niter
	454	if (prt_level.ge.1) print*,'Q2 THERMCEL_DQ 0'
	455
	456	! calcul du detrainement
	457	do k=1,nlay
	458	do ig=1,ngrid
	459	detr(ig,k)=fm(ig,k)-fm(ig,k+1)+entr(ig,k)
	460	! print*,'Q2 DQ ',detr(ig,k),fm(ig,k),entr(ig,k)
	461	!test
	462	if (detr(ig,k).lt.0.) then
	463	entr(ig,k)=entr(ig,k)-detr(ig,k)
	464	detr(ig,k)=0.
	465	! print*,'detr2<0!!!','ig=',ig,'k=',k,'f=',fm(ig,k),
	466	! s 'f+1=',fm(ig,k+1),'e=',entr(ig,k),'d=',detr(ig,k)
	467	endif
	468	if (fm(ig,k+1).lt.0.) then
	469	! print*,'fm2<0!!!'
	470	endif
	471	if (entr(ig,k).lt.0.) then
	472	! print*,'entr2<0!!!'
	473	endif
	474	enddo
	475	enddo
	476
	477	! calcul de la valeur dans les ascendances
	478	do ig=1,ngrid
	479	qa(ig,1)=q(ig,1)
	480	enddo
	481
	482	do k=2,nlay
	483	do ig=1,ngrid
	484	if ((fm(ig,k+1)+detr(ig,k))*ztimestep.gt. &
	485	& 1.e-5*masse(ig,k)) then
	486	qa(ig,k)=(fm(ig,k)qa(ig,k-1)+entr(ig,k)q(ig,k)) &
	487	& /(fm(ig,k+1)+detr(ig,k))
	488	else
	489	qa(ig,k)=q(ig,k)
	490	endif
	491	if (qa(ig,k).lt.0.) then
	492	! print*,'qa<0!!!'
	493	endif
	494	if (q(ig,k).lt.0.) then
	495	! print*,'q<0!!!'
	496	endif
	497	enddo
	498	enddo
	499
	500	! Calcul du flux subsident
	501
	502	do k=2,nlay
	503	do ig=1,ngrid
	504	#undef centre
	505	#ifdef centre
	506	wqd(ig,k)=fm(ig,k)0.5(q(ig,k-1)+q(ig,k))
	507	#else
	508
	509	#define plusqueun
	510	#ifdef plusqueun
	511	! Schema avec advection sur plus qu'une maille.
	512	zzm=masse(ig,k)/ztimestep
	513	if (fm(ig,k)>zzm) then
	514	wqd(ig,k)=zzmq(ig,k)+(fm(ig,k)-zzm)q(ig,k+1)
	515	else
	516	wqd(ig,k)=fm(ig,k)*q(ig,k)
	517	endif
	518	#else
	519	wqd(ig,k)=fm(ig,k)*q(ig,k)
	520	#endif
	521	#endif
	522
	523	if (wqd(ig,k).lt.0.) then
	524	! print*,'wqd<0!!!'
	525	endif
	526	enddo
	527	enddo
	528	do ig=1,ngrid
	529	wqd(ig,1)=0.
	530	wqd(ig,nlay+1)=0.
	531	enddo
	532
	533
	534	! Calcul des tendances
	535	do k=1,nlay
	536	do ig=1,ngrid
	537	q(ig,k)=q(ig,k)+(detr(ig,k)qa(ig,k)-entr(ig,k)q(ig,k) &
	538	& -wqd(ig,k)+wqd(ig,k+1)) &
	539	& *ztimestep/masse(ig,k)
	540	! if (dq(ig,k).lt.0.) then
	541	! print*,'dq<0!!!'
	542	! endif
	543	enddo
	544	enddo
	545
	546
	547	enddo
	548
	549
	550	! Calcul des tendances
	551	do k=1,nlay
	552	do ig=1,ngrid
	553	dq(ig,k)=(q(ig,k)-qold(ig,k))/ptimestep
	554	q(ig,k)=qold(ig,k)
	555	enddo
	556	enddo
	557
	558	return
	559	end
	560
	561
	562	#ifdef ISO
	563	subroutine thermcell_dq_o_iso(ngrid,nlay,ptimestep,fm,entr, &
	564	& masse,q,dq,qa,lev_out,xt,dxt,xta)
	565	#ifdef ISO
	566	use infotrac_phy, ONLY: ntraciso
	567	#ifdef ISOVERIF
	568	USE isotopes_mod, ONLY: iso_eau,iso_HDO
	569	USE isotopes_verif_mod, ONLY: iso_verif_egalite, &
	570	iso_verif_aberrant_enc_vect2D,iso_verif_egalite_vect2D
	571	#endif
	572	#endif
	573	implicit none
	574
	575	#include "iniprint.h"
	576	!=======================================================================
	577	!
	578	! Calcul du transport verticale dans la couche limite en presence
	579	! de "thermiques" explicitement representes
	580	! calcul du dq/dt une fois qu'on connait les ascendances
	581	!
	582	!=======================================================================
	583
	584	integer ngrid,nlay
	585
	586	real ptimestep
	587	real masse(ngrid,nlay),fm(ngrid,nlay+1)
	588	real entr(ngrid,nlay)
	589	real q(ngrid,nlay)
	590	real dq(ngrid,nlay)
	591	integer lev_out ! niveau pour les print
	592
	593	real qa(ngrid,nlay),detr(ngrid,nlay),wqd(ngrid,nlay+1)
	594
	595	real zzm
	596
	597	integer ig,k
	598	real cfl
	599
	600	real qold(ngrid,nlay)
	601	real ztimestep
	602	integer niter,iter
	603	CHARACTER (LEN=20) :: modname='thermcell_dq'
	604	CHARACTER (LEN=80) :: abort_message
	605
	606	! ifdef ISO
	607	real xt(ntraciso,ngrid,nlay)
	608	real dxt(ntraciso,ngrid,nlay)
	609	real xta(ntraciso,ngrid,nlay)
	610	real xtold(ntraciso,ngrid,nlay)
	611	real q_avant_entr(ngrid,nlay)
	612	integer ixt
	613	real wxtd(ntraciso,ngrid,nlay+1)
	614	! endif
	615
	616	! Calcul du critere CFL pour l'advection dans la subsidence
	617	cfl = 0.
	618	do k=1,nlay
	619	do ig=1,ngrid
	620	zzm=masse(ig,k)/ptimestep
	621	cfl=max(cfl,fm(ig,k)/zzm)
	622	if (entr(ig,k).gt.zzm) then
	623	print,'entr dt > m ',entr(ig,k)ptimestep,masse(ig,k)
	624	abort_message = ''
	625	CALL abort_gcm (modname,abort_message,1)
	626	endif
	627	enddo
	628	enddo
	629
	630	!IM 090508 print*,'CFL CFL CFL CFL ',cfl
	631
	632	#undef CFL
	633	#ifdef CFL
	634	! On subdivise le calcul en niter pas de temps.
	635	niter=int(cfl)+1
	636	#else
	637	niter=1
	638	#endif
	639
	640	ztimestep=ptimestep/niter
	641	qold=q
	642	!#ifdef ISO
	643	xtold=xt
	644	!#endif
	645
	646	do iter=1,niter
	647	if (prt_level.ge.1) print*,'Q2 THERMCEL_DQ 0'
	648
	649	! calcul du detrainement
	650	do k=1,nlay
	651	do ig=1,ngrid
	652	detr(ig,k)=fm(ig,k)-fm(ig,k+1)+entr(ig,k)
	653	! print*,'Q2 DQ ',detr(ig,k),fm(ig,k),entr(ig,k)
	654	!test
	655	if (detr(ig,k).lt.0.) then
	656	entr(ig,k)=entr(ig,k)-detr(ig,k)
	657	detr(ig,k)=0.
	658	! print*,'detr2<0!!!','ig=',ig,'k=',k,'f=',fm(ig,k),
	659	! s 'f+1=',fm(ig,k+1),'e=',entr(ig,k),'d=',detr(ig,k)
	660	endif
	661	if (fm(ig,k+1).lt.0.) then
	662	! print*,'fm2<0!!!'
	663	endif
	664	if (entr(ig,k).lt.0.) then
	665	! print*,'entr2<0!!!'
	666	endif
	667	enddo
	668	enddo
	669
	670	! calcul de la valeur dans les ascendances
	671	do ig=1,ngrid
	672	qa(ig,1)=q(ig,1)
	673	enddo
	674	!#ifdef ISO
	675	do ig=1,ngrid
	676	do ixt=1,ntraciso
	677	xta(ixt,ig,1)=xt(ixt,ig,1)
	678	enddo
	679	enddo
	680	!#endif
	681
	682	do k=2,nlay
	683	do ig=1,ngrid
	684	if ((fm(ig,k+1)+detr(ig,k))*ztimestep.gt. &
	685	& 1.e-5*masse(ig,k)) then
	686	qa(ig,k)=(fm(ig,k)qa(ig,k-1)+entr(ig,k)q(ig,k)) &
	687	& /(fm(ig,k+1)+detr(ig,k))
	688	!#ifdef ISO
	689	do ixt=1,ntraciso
	690	xta(ixt,ig,k)=(fm(ig,k)*xta(ixt,ig,k-1) &
	691	& +entr(ig,k)*xt(ixt,ig,k)) &
	692	& /(fm(ig,k+1)+detr(ig,k))
	693	enddo
	694	!#endif
	695	else
	696	qa(ig,k)=q(ig,k)
	697	!#ifdef ISO
	698	do ixt=1,ntraciso
	699	xta(ixt,ig,k)=xt(ixt,ig,k)
	700	enddo
	701	!#endif
	702	endif
	703	if (qa(ig,k).lt.0.) then
	704	! print*,'qa<0!!!'
	705	endif
	706	if (q(ig,k).lt.0.) then
	707	! print*,'q<0!!!'
	708	endif
	709	enddo
	710	enddo
	711
	712	!#ifdef ISO
	713	#ifdef ISOVERIF
	714	if (iso_HDO.gt.0) then
	715	call iso_verif_aberrant_enc_vect2D( &
	716	& xta,qa, &
	717	& 'thermcell_dq_iso 320, qa',ntraciso,ngrid,nlay)
	718	endif
	719	#endif
	720	!#endif
	721
	722	! Calcul du flux subsident
	723
	724	do k=2,nlay
	725	do ig=1,ngrid
	726	#undef centre
	727	#ifdef centre
	728	wqd(ig,k)=fm(ig,k)0.5(q(ig,k-1)+q(ig,k))
	729	!#ifdef ISO
	730	! attention à l'advection: on suppose que R varie linéairement entre les
	731	! couches, comme q.
	732	do ixt=1,ntraciso
	733	wxtd(ixt,ig,k)=wqd(ig,k)0.5( &
	734	& xt(ixt,ig,k-1)/q(ig,k-1)+xt(ixt,ig,k)/q(ig,k))
	735	enddo !do ixt=1,ntraciso
	736	!#endif
	737	#else
	738
	739	#define plusqueun
	740	#ifdef plusqueun
	741	! Schema avec advection sur plus qu'une maille.
	742	zzm=masse(ig,k)/ztimestep
	743	if (fm(ig,k)>zzm) then
	744	wqd(ig,k)=zzmq(ig,k)+(fm(ig,k)-zzm)q(ig,k+1)
	745	!#ifdef ISO
	746	do ixt=1,ntraciso
	747	wxtd(ixt,ig,k)=zzmxt(ixt,ig,k)+(fm(ig,k)-zzm)xt(ixt,ig,k+1)
	748	enddo !do ixt=1,ntraciso
	749	!#endif
	750	else
	751	wqd(ig,k)=fm(ig,k)*q(ig,k)
	752	!#ifdef ISO
	753	do ixt=1,ntraciso
	754	wxtd(ixt,ig,k)=fm(ig,k)*xt(ixt,ig,k)
	755	enddo !do ixt=1,ntraciso
	756	!#endif
	757	endif
	758	#else
	759	wqd(ig,k)=fm(ig,k)*q(ig,k)
	760	!#ifdef ISO
	761	do ixt=1,ntraciso
	762	wxtd(ixt,ig,k)=fm(ig,k)*xt(ixt,ig,k)
	763	enddo !do ixt=1,ntraciso
	764	!#endif
	765	#endif
	766	#endif
	767
	768	if (wqd(ig,k).lt.0.) then
	769	! print*,'wqd<0!!!'
	770	endif
	771	enddo
	772	enddo
	773	do ig=1,ngrid
	774	wqd(ig,1)=0.
	775	wqd(ig,nlay+1)=0.
	776	!#ifdef ISO
	777	do ixt=1,ntraciso
	778	wxtd(ixt,ig,1)=0.0
	779	wxtd(ixt,ig,nlay+1)=0.
	780	enddo !do ixt=1,ntraciso
	781	!#endif
	782	enddo
	783
	784
	785	! Calcul des tendances
	786	do k=1,nlay
	787	do ig=1,ngrid
	788	!#ifdef ISO
	789	q_avant_entr(ig,k)=q(ig,k)+(detr(ig,k)*qa(ig,k) &
	790	& -wqd(ig,k)+wqd(ig,k+1)) &
	791	& *ztimestep/masse(ig,k)
	792	!#endif
	793	q(ig,k)=q(ig,k)+(detr(ig,k)qa(ig,k)-entr(ig,k)q(ig,k) &
	794	& -wqd(ig,k)+wqd(ig,k+1)) &
	795	& *ztimestep/masse(ig,k)
	796	! if (dq(ig,k).lt.0.) then
	797	! print*,'dq<0!!!'
	798	! endif
	799	enddo
	800	enddo
	801
	802	!#ifdef ISO
	803	! Hypothèse: on ajoute le détrainement, et on fait la subsidence. Ensuite, on
	804	! retranche la perte vers l'entrainement en supposant que c'est un puit non
	805	! fractionnant
	806	do k=1,nlay
	807	do ig=1,ngrid
	808	do ixt=1,ntraciso
	809	xt(ixt,ig,k)=xt(ixt,ig,k)+(detr(ig,k)*xta(ixt,ig,k) &
	810	& -wxtd(ixt,ig,k)+wxtd(ixt,ig,k+1)) &
	811	& *ztimestep/masse(ig,k)
	812	xt(ixt,ig,k)=xt(ixt,ig,k)/q_avant_entr(ig,k)*q(ig,k)
	813	enddo !do ixt=1,ntraciso
	814	enddo !do ig=1,ngrid
	815	enddo !do k=1,nlay
	816	#ifdef ISOVERIF
	817	if (iso_eau.gt.0) then
	818	call iso_verif_egalite_vect2D( &
	819	& xt,q, &
	820	& 'thermcell_dq_iso 421',ntraciso,ngrid,nlay)
	821	endif
	822	if (iso_HDO.gt.0) then
	823	call iso_verif_aberrant_enc_vect2D( &
	824	& xt,q, &
	825	& 'thermcell_dq_iso 426',ntraciso,ngrid,nlay)
	826	endif
	827	#endif
	828	!#endif
	829
	830	enddo
	831
	832
	833	! Calcul des tendances
	834	do k=1,nlay
	835	do ig=1,ngrid
	836	dq(ig,k)=(q(ig,k)-qold(ig,k))/ptimestep
	837	q(ig,k)=qold(ig,k)
	838	!#ifdef ISO
	839	do ixt=1,ntraciso
	840	dxt(ixt,ig,k)=(xt(ixt,ig,k)-xtold(ixt,ig,k))/ptimestep
	841	xt(ixt,ig,k)=xtold(ixt,ig,k)
	842	enddo ! do ixt=1,ntraciso
	843	!#endif
	844	enddo
	845	enddo
	846
	847	return
	848	end
	849
	850	#endif
	851	! ifdef iso

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

Download in other formats:

Original Format