Context Navigation

moldiff_red.F90 @ 3461

Last change on this file since 3461 was 3221, checked in by flefevre, 9 months ago
Venus : a bit more cosmetics in molecular diffusion
File size: 39.8 KB

Rev	Line
[1442]	1	subroutine moldiff_red(ngrid,nlayer,nq,pplay,pplev,pt,pq,&
[2464]	2	ptimestep,pdteuv,pdtconduc,pdqdiff)
[1442]	3
	4	USE chemparam_mod
[1621]	5	USE infotrac_phy
[1442]	6	USE dimphy
	7
	8	implicit none
	9
	10	#include "comcstfi.h"
	11	#include "diffusion.h"
	12
	13
	14	!
	15	! Input/Output
	16	!
	17	integer,intent(in) :: ngrid ! number of atmospheric columns
	18	integer,intent(in) :: nlayer ! number of atmospheric layers
	19	integer,intent(in) :: nq ! number of advected tracers
	20	real ptimestep
	21	real pplay(ngrid,nlayer)
	22	real pplev(ngrid,nlayer+1)
	23	real pq(ngrid,nlayer,nq)
	24	! real pdq(ngrid,nlayer,nq)
	25	real pt(ngrid,nlayer)
	26	! real pdt(ngrid,nlayer)
	27	real pdteuv(ngrid,nlayer)
	28	real pdtconduc(ngrid,nlayer)
	29	real pdqdiff(ngrid,nlayer,nq)
	30	!
	31	! Local
	32	!
	33
	34	! real hco2(ncompdiff),ho
	35	integer,dimension(nq) :: indic_diff
	36	integer ig,iq,nz,l,k,n,nn,p,ij0
	37	integer istep,il,gcn,ntime,nlraf
	38	real*8 masse
	39	real*8 masseU,kBolt,RgazP,Rvenus,Grav,Mvenus,PI
	40	real*8 rho0,D0,T0,H0,time0,dZ,time,dZraf,tdiff,Zmin,Zmax
	41	real*8 FacEsc,invsgmu
	42	real*8 hp(nlayer)
	43	real*8 pp(nlayer)
	44	real*8 pint(nlayer)
	45	real*8 tt(nlayer),tnew(nlayer),tint(nlayer)
	46	real*8 zz(nlayer)
	47	real*8,dimension(:,:),allocatable,save :: qq,qnew,qint,FacMass
	48	real*8,dimension(:,:),allocatable,save :: rhoK,rhokinit
	49	real*8 rhoT(nlayer)
	50	real*8 dmmeandz(nlayer)
	51	real*8 massemoy(nlayer)
	52	real*8,dimension(:),allocatable :: Praf,Traf,Rraf,Mraf,Nraf,Draf,Hraf,Wraf
	53	real*8,dimension(:),allocatable :: Zraf,Tdiffraf
	54	real*8,dimension(:),allocatable :: Prafold,Mrafold
	55	real*8,dimension(:,:),allocatable :: Qraf,Rrafk,Nrafk
	56	real*8,dimension(:,:),allocatable :: Rrafkold
	57	real*8,dimension(:,:),allocatable :: Drafmol,Hrafmol,Wrafmol,Tdiffrafmol
	58	real*8,dimension(:),allocatable :: Atri,Btri,Ctri,Dtri,Xtri,Tad,Dad,Zad,rhoad
	59	real*8,dimension(:),allocatable :: alpha,beta,gama,delta,eps
	60	real*8,dimension(:),allocatable,save :: wi,Wad,Uthermal,Lambdaexo,Hspecie
	61	real*8,dimension(:),allocatable,save :: Mtot1,Mtot2,Mraf1,Mraf2
[2795]	62	integer,parameter :: nb_esp_diff=16
[2464]	63	character(len=20),dimension(nb_esp_diff) :: ListeDiff
[1442]	64	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	65	! tracer numbering in the molecular diffusion
	66	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	67	! We need the index of escaping species
	68
	69	integer,save :: i_esc_h2
	70	integer,save :: i_esc_h
	71	! vertical index limit for the molecular diffusion
	72	integer,save :: il0
	73
	74	! Tracer indexes in the GCM:
	75	! integer,save :: g_co2=0
	76	! integer,save :: g_co=0
	77	! integer,save :: g_o=0
	78	! integer,save :: g_o1d=0
	79	! integer,save :: g_o2=0
	80	! integer,save :: g_o3=0
	81	! integer,save :: g_h=0
	82	! integer,save :: g_h2=0
	83	! integer,save :: g_oh=0
	84	! integer,save :: g_ho2=0
	85	! integer,save :: g_h2o2=0
	86	! integer,save :: g_n2=0
	87	! integer,save :: g_ar=0
	88	! integer,save :: g_h2o=0
	89	! integer,save :: g_n=0
	90	! integer,save :: g_no=0
	91	! integer,save :: g_no2=0
	92	! integer,save :: g_n2d=0
	93	! integer,save :: g_oplus=0
	94	! integer,save :: g_hplus=0
	95
	96	integer,save :: ncompdiff
	97	integer,dimension(:),allocatable,save :: gcmind ! array of GCM indexes
	98	! Gab
	99	real,dimension(:),allocatable,save :: mol_tr ! array of mol mass traceurs
	100	integer ierr,compteur
	101
	102	logical,save :: firstcall=.true.
	103	! real abfac(ncompdiff)
	104	real,dimension(:,:),allocatable,save :: dij
	105	real,save :: step
	106
	107	! Initializations at first call
[3218]	108
[1442]	109	if (firstcall) then
	110
[3218]	111	! list of diffused species (if present)
[1442]	112
[3218]	113	ListeDiff(1) = 'co2'
	114	ListeDiff(2) = 'o'
	115	ListeDiff(3) = 'n2'
	116	ListeDiff(4) = 'ar'
	117	ListeDiff(5) = 'co'
	118	ListeDiff(6) = 'h2'
	119	ListeDiff(7) = 'h'
	120	ListeDiff(8) = 'd2'
	121	ListeDiff(9) = 'hd'
	122	ListeDiff(10) = 'd'
	123	ListeDiff(11) = 'o2'
	124	ListeDiff(12) = 'h2o'
	125	ListeDiff(13) = 'o3'
	126	ListeDiff(14) = 'n'
	127	ListeDiff(15) = 'he'
	128	ListeDiff(16) = 'n2d'
[1442]	129
[3218]	130	i_esc_h = 1000
	131	i_esc_h2 = 1000
[1442]	132
[3218]	133	! search for species that can be diffused
[1442]	134
[3218]	135	ncompdiff = 0
	136	indic_diff(1:nq) = 0
[1442]	137
[3218]	138	do nn = 1,nq
	139	do n = 1,nb_esp_diff
	140	if (ListeDiff(n) == tname(nn)) then
	141	indic_diff(nn) = 1
	142	end if
	143	end do
[1442]	144
[3218]	145	if (indic_diff(nn) == 1) then
	146	print*,'species ', tname(nn), 'diffused in moldiff_red'
	147	ncompdiff = ncompdiff + 1
	148	end if
	149	end do
[1442]	150
[3218]	151	print*,'number of diffused species:',ncompdiff
	152	allocate(gcmind(ncompdiff))
	153	allocate(mol_tr(ncompdiff))
	154
	155	! store gcm indexes in gcmind and molar masses in mol_tr
	156
	157	n = 0
	158	do nn = 1,nq
	159	if (indic_diff(nn) == 1) then
	160	n = n + 1
	161	gcmind(n) = nn
	162	mol_tr(n) = M_tr(nn)
	163	if (tname(nn) == 'h') i_esc_h = n
	164	if (tname(nn) == 'h2') i_esc_h2 = n
	165	end if
	166	end do
	167
[3221]	168	! find vertical index above which diffusion is computed
[1442]	169
[3221]	170	do l = 1,nlayer
	171	if (pplay(1,l) > Pdiff) then
	172	il0 = l
	173	end if
	174	end do
[1442]	175
[3221]	176	il0 = il0 + 1
[1442]	177
[3221]	178	print*,'vertical index for diffusion',il0,pplay(1,il0)
[1442]	179
[3221]	180	allocate(dij(ncompdiff,ncompdiff))
[1442]	181
[3221]	182	call moldiffcoeff_red(dij,indic_diff,gcmind,ncompdiff)
[1442]	183
[3221]	184	print*,'MOLDIFF EXO'
[1442]	185
[3221]	186	! allocation of arrays depending on the number of diffused species
	187
	188	allocate(qq(nlayer,ncompdiff))
	189	allocate(qnew(nlayer,ncompdiff))
	190	allocate(qint(nlayer,ncompdiff))
	191	allocate(FacMass(nlayer,ncompdiff))
	192	allocate(rhok(nlayer,ncompdiff))
	193	allocate(rhokinit(nlayer,ncompdiff))
	194
	195	allocate(wi(ncompdiff))
	196	allocate(wad(ncompdiff))
	197	allocate(uthermal(ncompdiff))
	198	allocate(lambdaexo(ncompdiff))
	199	allocate(Hspecie(ncompdiff))
	200	allocate(Mtot1(ncompdiff))
	201	allocate(Mtot2(ncompdiff))
	202	allocate(Mraf1(ncompdiff))
	203	allocate(Mraf2(ncompdiff))
	204
	205	firstcall = .false.
	206	step = 1
	207
	208	end if ! of if (firstcall)
	209
[1442]	210	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	211
	212	masseU=1.660538782d-27
	213	kbolt=1.3806504d-23
	214	RgazP=8.314472
	215	PI =3.141592653
	216	g=8.87d0
	217	Rvenus=6051800d0 ! Used to compute escape parameter no need to be more accurate
	218	Grav=6.67d-11
	219	Mvenus=4.86d24
	220	ij0=6000 ! For test
	221	invsgmu=1d0/g/masseU
	222
	223	! print*,'moldiff',i_esc_h2,i_esc_h,ncompdiff
	224	do ig=1,ngrid
	225	pp=dble(pplay(ig,:))
	226
	227	!!!
	228	!! =======> Gab 29 oct 2014
[1591]	229	!!$$$ THIS UPDATE IS ALREADY DONE in physique (t_seri & tr_seri) $$$$$$
[1442]	230
	231	! Update the temperature modified by other processes
	232
	233	do l=1,nlayer
	234	tt(l)=pt(ig,l)
	235	enddo ! of do l=1,nlayer
	236
	237
	238	! CALL TMNEW(pt(ig,:),pdt(ig,:),pdtconduc(ig,:),pdteuv(ig,:) &
	239	! & ,tt,ptimestep,nlayer,ig)
	240	! do l=1,nlayer
	241	! tt(l)=pt(ig,l)1D0+(pdt(ig,l)dble(ptimestep)+ &
	242	! pdtconduc(ig,l)*dble(ptimestep)+ &
	243	! pdteuv(ig,l)*dble(ptimestep))
	244	! ! to cach Nans...
	245	! if (tt(l).ne.tt(l)) then
	246	! print*,'Err TMNEW',ig,l,tt(l),pt(ig,l), &
	247	! pdt(ig,l),pdtconduc(ig,l),pdteuv(ig,l),dble(ptimestep)
	248	! endif
	249	! enddo ! of do l=1,nlayer
	250
	251	! Update the mass mixing ratios modified by other processes
	252
	253	do iq=1,ncompdiff
	254	do l=1,nlayer
	255	qq(l,iq)=pq(ig,l,gcmind(iq))
	256	qq(l,iq)=max(qq(l,iq),1d-30)
	257	enddo ! of do l=1,nlayer
	258	enddo ! of do iq=1,ncompdiff
	259	! STOP
	260
	261	! Compute the Pressure scale height
	262
	263	CALL HSCALE(pp,hp,nlayer)
	264
	265	! Compute the atmospheric mass (in Dalton)
	266
	267	CALL MMOY(massemoy,mol_tr,qq,gcmind,nlayer,ncompdiff)
	268
	269	! Compute the vertical gradient of atmospheric mass
	270
	271	CALL DMMOY(massemoy,hp,dmmeandz,nlayer)
	272
	273	! Compute the altitude of each layer
	274
	275	CALL ZVERT(pp,tt,massemoy,zz,nlayer,ig)
	276
	277	! Compute the total mass density (kg/m3)
	278
	279	CALL RHOTOT(pp,tt,massemoy,qq,RHOT,RHOK,nlayer,ncompdiff)
	280	RHOKINIT=RHOK
	281
	282	! Compute total mass of each specie before new grid
	283	! For conservation tests
	284	! The conservation is not always fulfilled especially
	285	! for species very far from diffusion equilibrium "photochemical species"
	286	! e.g. OH, O(1D)
	287
	288	Mtot1(1:ncompdiff)=0d0
	289
	290	do l=il0,nlayer
	291	do nn=1,ncompdiff
	292	Mtot1(nn)=Mtot1(nn)+1d0/gqq(l,nn) &
	293	& (dble(pplev(ig,l))-dble(pplev(ig,l+1)))
	294	enddo
	295	enddo
	296
	297	Zmin=zz(il0)
	298	Zmax=zz(nlayer)
	299
	300
	301	! If Zmax > 4000 km there is a problem / stop
	302
	303	if (Zmax .gt. 4000000.) then
	304	Print*,'Zmax too high',ig,zmax,zmin
	305	do l=1,nlayer
	306	print*,'old',zz(l),pt(ig,l),pdteuv(ig,l)
	307	print*,'l',l,rhot(l),tt(l),pp(l),massemoy(l),qq(l,:)
	308	enddo
	309	stop
	310	endif
	311
	312	! The number of diffusion layers is variable
	313	! and fixed by the resolution (dzres) specified in diffusion.h
	314	! I fix a minimum number of layers 40
	315
	316	nlraf=int((Zmax-Zmin)/1000./dzres)+1
	317	if (nlraf .le. 40) nlraf=40
	318
	319	! if (nlraf .ge. 200) print*,ig,nlraf,Zmin,Zmax
	320
	321	! allocate arrays:
	322	allocate(Praf(nlraf),Traf(nlraf),Rraf(nlraf),Mraf(nlraf))
	323	allocate(Nraf(nlraf),Draf(nlraf),Hraf(nlraf),Wraf(nlraf))
	324	allocate(Zraf(nlraf),Tdiffraf(nlraf))
	325	allocate(Prafold(nlraf),Mrafold(nlraf))
	326	allocate(Qraf(nlraf,ncompdiff),Rrafk(nlraf,ncompdiff),Nrafk(nlraf,ncompdiff))
	327	allocate(Rrafkold(nlraf,ncompdiff))
	328	allocate(Drafmol(nlraf,ncompdiff),Hrafmol(nlraf,ncompdiff))
	329	allocate(Wrafmol(nlraf,ncompdiff),Tdiffrafmol(nlraf,ncompdiff))
	330	allocate(Atri(nlraf),Btri(nlraf),Ctri(nlraf),Dtri(nlraf),Xtri(nlraf))
	331	allocate(Tad(nlraf),Dad(nlraf),Zad(nlraf),rhoad(nlraf))
	332	allocate(alpha(nlraf),beta(nlraf),gama(nlraf),delta(nlraf),eps(nlraf))
	333
	334	! before beginning, I use a better vertical resolution above il0,
	335	! altitude grid reinterpolation
	336	! The diffusion is solved on an altitude grid with constant step dz
	337
	338	CALL UPPER_RESOL(pp,tt,zz,massemoy,RHOT,RHOK, &
	339	& qq,mol_tr,gcmind,Praf,Traf,Qraf,Mraf,Zraf, &
	340	& Nraf,Nrafk,Rraf,Rrafk,il0,nlraf,ncompdiff,nlayer,ig)
	341
	342	Prafold=Praf
	343	Rrafkold=Rrafk
	344	Mrafold=Mraf
	345
	346	! We reddo interpolation of the gcm grid to estimate mass loss due to interpolation processes.
	347
	348	CALL GCMGRID_P(Zraf,Praf,Qraf,Traf,Nrafk,Rrafk,qq,qint,tt,tint &
	349	& ,pp,mol_tr,gcmind,nlraf,ncompdiff,nlayer,ig)
	350
	351	! We compute the mass correction factor of each specie at each pressure level
	352
	353	CALL CORRMASS(qq,qint,FacMass,nlayer,ncompdiff)
	354
	355	! Altitude step
	356
	357	Dzraf=Zraf(2)-Zraf(1)
	358
	359	! Total mass computed on the altitude grid
	360
	361	Mraf1(1:ncompdiff)=0d0
	362	do nn=1,ncompdiff
	363	do l=1,nlraf
	364	Mraf1(nn)=Mraf1(nn)+Rrafk(l,nn)*Dzraf
	365	enddo
	366	enddo
	367
	368	! Reupdate values for mass conservation
	369
	370	! do l=1,nlraf
	371	! print*,'test',l,Nraf(l),Praf(l)
	372	! do nn=1,ncompdiff
	373	! Rrafk(l,nn)=RrafK(l,nn)*Mtot1(nn)/Mraf1(nn)
	374	! enddo
	375	! Rraf(l)=sum(Rrafk(l,:))
	376	! do nn=1,ncompdiff
	377	! Qraf(l,nn)=Rrafk(l,nn)/Rraf(l)
	378	! Nrafk(l,nn)=Rrafk(l,nn)/dble(mol_tr(gcmind(nn)))/masseU
	379	! enddo
	380	! Mraf(l)=1d0/sum(Qraf(l,:)/dble(mol_tr(gcmind(:))))
	381	! Nraf(l)=Rraf(l)/Mraf(l)/masseU
	382	! Praf(l)=kboltTraf(l)Nraf(l)
	383	! print*,'test',l,Nraf(l),Praf(l)
	384	! enddo
	385
	386	! do l=1,nlayer
	387	! print*,'l',l,zz(l),pp(l),tt(l),sum(qq(l,:)),massemoy(l)
	388	! enddo
	389
	390	! The diffusion is computed above il0 computed from Pdiff in diffusion.h
	391	! No change below il0
	392
	393	do l=1,nlayer
	394	qnew(l,:)=qq(l,:) ! No effet below il0
	395	enddo
	396
	397	! all species treated independently
	398
	399	! Upper boundary velocity
	400	! Jeans escape for H and H2
	401	! Comparison with and without escape still to be done
	402	! No escape for other species
	403
	404
	405	do nn=1,ncompdiff
	406	Uthermal(nn)=sqrt(2d0kboltTraf(nlraf)/masseU/ &
	407	& dble(mol_tr(nn)))
	408	Lambdaexo(nn)=masseUdble(mol_tr(nn))Grav*Mvenus/ &
	409	& (Rvenus+Zraf(nlraf))/kbolt/Traf(nlraf)
	410	wi(nn)=0D0
	411	if (nn .eq. i_esc_h .or. nn .eq. i_esc_h2) then
	412	wi(nn)=Uthermal(nn)/2d0/sqrt(PI)exp(-Lambdaexo(nn)) &
	413	& (Lambdaexo(nn)+1d0)
	414	endif
	415	enddo
	416
	417	! Compute time step for diffusion
	418
	419	! Loop on species
	420
	421	T0=Traf(nlraf)
	422	rho0=1d0
	423
	424	do nn=1,ncompdiff
	425	masse=dble(mol_tr(nn))
	426	! DIffusion coefficient
	427	CALL DCOEFF(nn,dij,Praf,Traf,Nraf,Nrafk,Draf,nlraf,ncompdiff)
[2836]	428	!Draf(:) = max(100.,Draf(:))
	429	!Draf(:) = 0.01*Draf(:)
[1442]	430	Drafmol(:,nn)=Draf(:)
	431	! Scale height of the density of the specie
	432	CALL HSCALEREAL(nn,Nrafk,Dzraf,Hraf,nlraf,ncompdiff)
	433	Hrafmol(:,nn)=Hraf(:)
	434	! Hspecie(nn)=kboltT0/masseinvsgmu
	435	! Computation of the diffusion vertical velocity of the specie
	436	CALL VELVERT(nn,Traf,Hraf,Draf,Dzraf,masse,Wraf,nlraf)
	437	Wrafmol(:,nn)=Wraf(:)
	438	! Computation of the diffusion time
	439	CALL TIMEDIFF(nn,Hraf,Wraf,Tdiffraf,nlraf)
	440	Tdiffrafmol(:,nn)=Tdiffraf
	441	enddo
	442	! We use a lower time step
	443	Tdiff=minval(Tdiffrafmol)
	444	Tdiff=minval(Tdiffrafmol(nlraf,:))*Mraf(nlraf)
	445	! Number of time step
	446
	447	! Some problems when H is dominant
	448	! The time step is chosen function of atmospheric mass at higher level
	449	! It is not very nice
	450
	451	!!!! TEST GAB 14jan15 : increase/decrease number of time steps for diffusion (reduce/augmente tdiff)
	452	! Tdiff=Tdiff/2.
	453
[1591]	454	if (tdiff .lt. tdiffminMraf(nlraf)) tdiff=tdiffminMraf(nlraf)
	455
	456	! if (tdiff .lt. tdiffmin*Mraf(nlraf)) then
	457	! print, 'Moldiff L454', tdiff, ptimestep, tdiffminMraf(nlraf)
	458	! endif
	459
	460	! JYC + GG aout 2015 add this condition below
	461	if (tdiff .ge. ptimestep) tdiff=ptimestep
	462	! Number of time step
	463	ntime=anint(dble(ptimestep)/tdiff)
	464	!print*,'ptime',ig,tdiff,ntime,Mraf(nlraf)
	465
[1442]	466	! Adimensionned temperature
	467
	468	do l=1,nlraf
	469	Tad(l)=Traf(l)/T0
	470	enddo
	471
	472	do istep=1,ntime
	473	do nn=1,ncompdiff
	474
	475	Draf(1:nlraf)=Drafmol(1:nlraf,nn)
	476
	477	! Parameters to adimension the problem
	478
	479	H0=kboltT0/dble(mol_tr(nn))invsgmu
	480	D0=Draf(nlraf)
	481	Time0=H0*H0/D0
	482	Time=Tdiff/Time0
	483
	484	! Adimensions
	485
	486	do l=1,nlraf
	487	Dad(l)=Draf(l)/D0
	488	Zad(l)=Zraf(l)/H0
	489	enddo
	490	Wad(nn)=wi(nn)*Time0/H0
	491	DZ=Zad(2)-Zad(1)
	492	FacEsc=exp(-wad(nn)*DZ/Dad(nlraf-1))
	493
	494	do l=1,nlraf
	495	RhoAd(l)=Rrafk(l,nn)/rho0
	496	enddo
	497
	498	! Compute intermediary coefficients
	499
	500	CALL DIFFPARAM(Tad,Dad,DZ,alpha,beta,gama,delta,eps,nlraf)
	501
	502	! First way to include escape need to be validated
	503	! Compute escape factor (exp(-ueff*dz/D(nz)))
	504
	505	! Compute matrix coefficients
	506
	507	CALL MATCOEFF(alpha,beta,gama,delta,eps,Dad,rhoAd,FacEsc, &
	508	& dz,time,Atri,Btri,Ctri,Dtri,nlraf)
	509
	510	Xtri(:)=0D0
	511
	512	! SOLVE SYSTEM
	513
	514	CALL tridag(Atri,Btri,Ctri,Dtri,Xtri,nlraf)
	515
	516	! Xtri=rhoAd
	517
	518	! if (ig .eq. ij0 .and. (nn .eq. 1 .or. nn .eq. 5 .or. nn .eq. 6 .or. nn .eq. 16)) then
	519	! do l=1,nlraf
	520	! if (Xtri(l) .lt. 0.) then
	521	! print,'l',l,rhoAd(l)rho0,Xtri(l)*rho0,nn,Tad(l),Zad(l),Dad(l)
	522	! stop
	523	! endif
	524	! enddo
	525	! endif
	526
	527	! Check mass conservation of speci
	528
	529	! CALL CheckMass(rhoAd,Xtri,nlraf,nn)
	530
	531	! Update mass density of the species
	532
	533	do l=1,nlraf
	534	Rrafk(l,nn)=rho0*Xtri(l)
	535	if (Rrafk(l,nn) .ne. Rrafk(l,nn) .or. &
	536	& Rrafk(l,nn) .lt. 0 .and. nn .eq. 16) then
	537
	538	! Test if n(CO2) < 0 skip diffusion (should never happen)
	539
	540	print*,'pb moldiff',istep,ig,l,nn,Rrafk(l,nn),tdiff,&
	541	& rho0*Rhoad(l),Zmin,Zmax,ntime
	542	print*,'Atri',Atri
	543	print*,'Btri',Btri
	544	print*,'Ctri',Ctri
	545	print*,'Dtri',Dtri
	546	print*,'Xtri',Xtri
	547	print*,'alpha',alpha
	548	print*,'beta',beta
	549	print*,'gama',gama
	550	print*,'delta',delta
	551	print*,'eps',eps
	552	print*,'Dad',Dad
	553	print*,'Temp',Traf
	554	print*,'alt',Zraf
	555	print*,'Mraf',Mraf
	556	stop
	557	! pdqdiff(1:ngrid,1:nlayer,1:nq)=0.
	558	! return
	559	! Rrafk(l,nn)=1D-30*Rraf(l)
	560	Rrafk(l,nn)=rho0*Rhoad(l)
	561
	562	endif
	563
	564	enddo
	565
	566	enddo ! END Species Loop
	567
	568	! Update total mass density
	569
	570	do l=1,nlraf
	571	Rraf(l)=sum(Rrafk(l,:))
	572	enddo
	573
	574	! Compute new mass average at each altitude level and new pressure
	575
	576	do l=1,nlraf
	577	do nn=1,ncompdiff
	578	Qraf(l,nn)=Rrafk(l,nn)/Rraf(l)
	579	Nrafk(l,nn)=Rrafk(l,nn)/dble(mol_tr(nn))/masseU
	580	enddo
	581	Mraf(l)=1d0/sum(Qraf(l,:)/dble(mol_tr(:)))
	582	Nraf(l)=Rraf(l)/Mraf(l)/masseU
	583	Praf(l)=Nraf(l)kboltTraf(l)
	584	enddo
	585
	586	enddo ! END time Loop
	587
	588	! Compute the total mass of each species to check mass conservation
	589
	590	Mraf2(1:ncompdiff)=0d0
	591	do nn=1,ncompdiff
	592	do l=1,nlraf
	593	Mraf2(nn)=Mraf2(nn)+Rrafk(l,nn)*Dzraf
	594	enddo
	595	enddo
	596
	597	! print*,'Mraf',Mraf2
	598
	599	! Reinterpolate values on the GCM pressure levels
	600
	601	CALL GCMGRID_P2(Zraf,Praf,Qraf,Traf,Nrafk,Rrafk,qq,qnew,tt,tnew,&
	602	& pp,mol_tr,gcmind,nlraf,ncompdiff,nlayer,FacMass,ig)
	603
[1591]	604	!!! Call added by JY+GG Aout 2014
	605	CALL MMOY(massemoy,mol_tr,qnew,gcmind,nlayer,ncompdiff)
	606
[1442]	607	CALL RHOTOT(pp,tt,massemoy,qnew,RHOT,RHOK,nlayer,ncompdiff)
	608
	609	if (ig .eq. ij0) then
	610	do l=il0,nlayer
	611	write(*,'(i2,1x,19(e12.4,1x))') l,zz(l),tt(l),RHOK(l,1)/sum(RHOK(l,:)),RHOKINIT(l,1)/sum(RHOKINIT(l,:)),&
	612	& RHOK(l,2)/sum(RHOK(l,:)),RHOKINIT(l,2)/sum(RHOKINIT(l,:)),&
	613	& RHOK(l,6)/sum(RHOK(l,:)),RHOKINIT(l,6)/sum(RHOKINIT(l,:)),&
	614	& RHOK(l,5)/sum(RHOK(l,:)),RHOKINIT(l,5)/sum(RHOKINIT(l,:)),&
	615	& RHOK(l,7)/sum(RHOK(l,:)),RHOKINIT(l,7)/sum(RHOKINIT(l,:))
	616	enddo
	617	endif
	618
	619	! Compute total mass of each specie on the GCM vertical grid
	620
	621	Mtot2(1:ncompdiff)=0d0
	622
	623	do l=il0,nlayer
	624	do nn=1,ncompdiff
	625	Mtot2(nn)=Mtot2(nn)+1d0/gqnew(l,nn) &
	626	& (dble(pplev(ig,l))-dble(pplev(ig,l+1)))
	627	enddo
	628	enddo
	629
[2464]	630	! DEBUG
	631	! print*,"rapport MASSE: ",Mtot2(:)/Mtot1(:)
	632
[1442]	633	! Check mass conservation of each specie on column
	634
	635	! do nn=1,ncompdiff
	636	! CALL CheckMass2(qq,qnew,pplev(ig,:),il0,nlayer,nn,ncompdiff)
	637	! enddo
	638
	639	! Compute the diffusion trends du to diffusion
	640
	641	do l=1,nlayer
	642	do nn=1,ncompdiff
	643	pdqdiff(ig,l,gcmind(nn))=(qnew(l,nn)-qq(l,nn))/ptimestep
	644	enddo
	645	enddo
	646
	647	! deallocation des tableaux
	648
	649	deallocate(Praf,Traf,Rraf,Mraf)
	650	deallocate(Nraf,Draf,Hraf,Wraf)
	651	deallocate(Zraf,Tdiffraf)
	652	deallocate(Prafold,Mrafold)
	653	deallocate(Qraf,Rrafk,Nrafk)
	654	deallocate(Rrafkold)
	655	deallocate(Drafmol,Hrafmol)
	656	deallocate(Wrafmol,Tdiffrafmol)
	657	deallocate(Atri,Btri,Ctri,Dtri,Xtri)
	658	deallocate(Tad,Dad,Zad,rhoad)
	659	deallocate(alpha,beta,gama,delta,eps)
	660
	661
	662	enddo ! ig loop
	663
	664
	665	return
	666	end
	667
	668	! ********************************************************************
	669	! ********************************************************************
	670	! ********************************************************************
	671
	672	! JYC subtroutine solving MX = Y where M is defined as a block tridiagonal
	673	! matrix (Thomas algorithm), tested on a example
	674
	675	subroutine tridagbloc(M,F,X,n1,n2)
	676	parameter (nmax=100)
	677	real8 M(n1n2,n1n2),F(n1n2),X(n1*n2)
	678	real*8 A(n1,n1,n2),B(n1,n1,n2),C(n1,n1,n2),D(n1,n2)
	679	real*8 at(n1,n1),bt(n1,n1),ct(n1,n1),dt(n1),gamt(n1,n1),y(n1,n1)
	680	real*8 alf(n1,n1),gam(n1,n1,n2),alfinv(n1,n1)
	681	real*8 uvec(n1,n2),uvect(n1),vvect(n1),xt(n1)
	682	real*8 indx(n1)
	683	real*8 rhu
	684	integer n1,n2
	685	integer i,p,q
	686
	687	X(:)=0.
	688	! Define the bloc matrix A,B,C and the vector D
	689	A(1:n1,1:n1,1)=M(1:n1,1:n1)
	690	C(1:n1,1:n1,1)=M(1:n1,n1+1:2*n1)
	691	D(1:n1,1)=F(1:n1)
	692
	693	do i=2,n2-1
	694	A(1:n1,1:n1,i)=M((i-1)n1+1:in1,(i-1)n1+1:in1)
	695	B(1:n1,1:n1,i)=M((i-1)n1+1:in1,(i-2)n1+1:(i-1)n1)
	696	C(1:n1,1:n1,i)=M((i-1)n1+1:in1,in1+1:(i+1)n1)
	697	D(1:n1,i)=F((i-1)n1+1:in1)
	698	enddo
	699	A(1:n1,1:n1,n2)=M((n2-1)n1+1:n2n1,(n2-1)n1+1:n2n1)
	700	B(1:n1,1:n1,n2)=M((n2-1)n1+1:n2n1,(n2-2)n1+1:(n2-1)n1)
	701	D(1:n1,n2)=F((n2-1)n1+1:n2n1)
	702
	703	! Initialization
	704	y(:,:)=0.
	705	do i=1,n1
	706	y(i,i)=1.
	707	enddo
	708
	709	at(:,:)=A(:,:,1)
	710	ct(:,:)=C(:,:,1)
	711	dt(:)=D(:,1)
	712	call ludcmp(at,n1,n1,indx,rhu,ierr)
	713	do p=1,n1
	714	call lubksb(at,n1,n1,indx,y(1,p))
	715	do q=1,n1
	716	alfinv(q,p)=y(q,p)
	717	enddo
	718	enddo
	719	gamt=matmul(alfinv,ct)
	720	gam(:,:,1)=gamt(:,:)
	721	uvect=matmul(alfinv,dt)
	722	uvec(:,1)=uvect
	723
	724	do i=2,n2-1
	725	y(:,:)=0.
	726	do j=1,n1
	727	y(j,j)=1.
	728	enddo
	729	bt(:,:)=B(:,:,i)
	730	at(:,:)=A(:,:,i)-matmul(bt,gamt)
	731	ct(:,:)=C(:,:,i)
	732	dt(:)=D(:,i)
	733	call ludcmp(at,n1,n1,indx,rhu,ierr)
	734	do p=1,n1
	735	call lubksb(at,n1,n1,indx,y(1,p))
	736	do q=1,n1
	737	alfinv(q,p)=y(q,p)
	738	enddo
	739	enddo
	740	gamt=matmul(alfinv,ct)
	741	gam(:,:,i)=gamt
	742	vvect=dt-matmul(bt,uvect)
	743	uvect=matmul(alfinv,vvect)
	744	uvec(:,i)=uvect
	745	enddo
	746	bt=B(:,:,n2)
	747	dt=D(:,n2)
	748	at=A(:,:,n2)-matmul(bt,gamt)
	749	vvect=dt-matmul(bt,uvect)
	750	y(:,:)=0.
	751	do j=1,n1
	752	y(j,j)=1.
	753	enddo
	754	call ludcmp(at,n1,n1,indx,rhu,ierr)
	755	do p=1,n1
	756	call lubksb(at,n1,n1,indx,y(1,p))
	757	do q=1,n1
	758	alfinv(q,p)=y(q,p)
	759	enddo
	760	enddo
	761	xt=matmul(alfinv,vvect)
	762	X((n2-1)n1+1 :n1n2)=xt
	763	do i=n2-1,1,-1
	764	gamt=gam(:,:,i)
	765	xt=X(in1+1:n1n2)
	766	uvect=uvec(:,i)
	767	vvect=matmul(gamt,xt)
	768	X((i-1)n1+1:in1)=uvect-vvect
	769	enddo
	770
	771	end
	772
	773	subroutine tridag(a,b,c,r,u,n)
	774	! parameter (nmax=4000)
	775	! dimension gam(nmax),a(n),b(n),c(n),r(n),u(n)
	776	real*8 gam(n),a(n),b(n),c(n),r(n),u(n)
	777	if(b(1).eq.0.)then
	778	stop 'tridag: error: b(1)=0 !!! '
	779	endif
	780	bet=b(1)
	781	u(1)=r(1)/bet
	782	do 11 j=2,n
	783	gam(j)=c(j-1)/bet
	784	bet=b(j)-a(j)*gam(j)
	785	if(bet.eq.0.) then
	786	stop 'tridag: error: bet=0 !!! '
	787	endif
	788	u(j)=(r(j)-a(j)*u(j-1))/bet
	789	11 continue
	790	do 12 j=n-1,1,-1
	791	u(j)=u(j)-gam(j+1)*u(j+1)
	792	12 continue
	793	return
	794	end
	795
	796	! ********************************************************************
	797	! ********************************************************************
	798	! ********************************************************************
	799
	800	SUBROUTINE LUBKSB(A,N,NP,INDX,B)
	801
	802	implicit none
	803
	804	integer i,j,n,np,ii,ll
	805	real*8 sum
	806	real*8 a(np,np),indx(np),b(np)
	807
	808	! DIMENSION A(NP,NP),INDX(N),B(N)
	809	II=0
	810	DO 12 I=1,N
	811	LL=INDX(I)
	812	SUM=B(LL)
	813	B(LL)=B(I)
	814	IF (II.NE.0)THEN
	815	DO 11 J=II,I-1
	816	SUM=SUM-A(I,J)*B(J)
	817	11 CONTINUE
	818	ELSE IF (SUM.NE.0.) THEN
	819	II=I
	820	ENDIF
	821	B(I)=SUM
	822	12 CONTINUE
	823	DO 14 I=N,1,-1
	824	SUM=B(I)
	825	IF(I.LT.N)THEN
	826	DO 13 J=I+1,N
	827	SUM=SUM-A(I,J)*B(J)
	828	13 CONTINUE
	829	ENDIF
	830	B(I)=SUM/A(I,I)
	831	14 CONTINUE
	832	RETURN
	833	END
	834
	835	! ********************************************************************
	836	! ********************************************************************
	837	! ********************************************************************
	838
	839	SUBROUTINE LUDCMP(A,N,NP,INDX,D,ierr)
	840
	841	implicit none
	842
	843	integer n,np,nmax,i,j,k,imax
	844	real*8 d,tiny,aamax
	845	real*8 a(np,np),indx(np)
	846	integer ierr ! error =0 if OK, =1 if problem
	847
	848	PARAMETER (NMAX=100,TINY=1.0E-20)
	849	! DIMENSION A(NP,NP),INDX(N),VV(NMAX)
	850	real*8 sum,vv(nmax),dum
	851
	852	D=1.
	853	DO 12 I=1,N
	854	AAMAX=0.
	855	DO 11 J=1,N
	856	IF (ABS(A(I,J)).GT.AAMAX) AAMAX=ABS(A(I,J))
	857	11 CONTINUE
	858	IF (AAMAX.EQ.0.) then
	859	write(,) 'In moldiff: Problem in LUDCMP with matrix A'
	860	write(,) 'Singular matrix ?'
	861	write(,) 'Matrix A = ', A
	862	! stop
	863	! TO DEBUG :
	864	ierr =1
	865	return
	866	! stop
	867	END IF
	868
	869	VV(I)=1./AAMAX
	870	12 CONTINUE
	871	DO 19 J=1,N
	872	IF (J.GT.1) THEN
	873	DO 14 I=1,J-1
	874	SUM=A(I,J)
	875	IF (I.GT.1)THEN
	876	DO 13 K=1,I-1
	877	SUM=SUM-A(I,K)*A(K,J)
	878	13 CONTINUE
	879	A(I,J)=SUM
	880	ENDIF
	881	14 CONTINUE
	882	ENDIF
	883	AAMAX=0.
	884	DO 16 I=J,N
	885	SUM=A(I,J)
	886	IF (J.GT.1)THEN
	887	DO 15 K=1,J-1
	888	SUM=SUM-A(I,K)*A(K,J)
	889	15 CONTINUE
	890	A(I,J)=SUM
	891	ENDIF
	892	DUM=VV(I)*ABS(SUM)
	893	IF (DUM.GE.AAMAX) THEN
	894	IMAX=I
	895	AAMAX=DUM
	896	ENDIF
	897	16 CONTINUE
	898	IF (J.NE.IMAX)THEN
	899	DO 17 K=1,N
	900	DUM=A(IMAX,K)
	901	A(IMAX,K)=A(J,K)
	902	A(J,K)=DUM
	903	17 CONTINUE
	904	D=-D
	905	VV(IMAX)=VV(J)
	906	ENDIF
	907	INDX(J)=IMAX
	908	IF(J.NE.N)THEN
	909	IF(A(J,J).EQ.0.)A(J,J)=TINY
	910	DUM=1./A(J,J)
	911	DO 18 I=J+1,N
	912	A(I,J)=A(I,J)*DUM
	913	18 CONTINUE
	914	ENDIF
	915	19 CONTINUE
	916	IF(A(N,N).EQ.0.)A(N,N)=TINY
	917	ierr =0
	918	RETURN
	919	END
	920
	921	SUBROUTINE TMNEW(T1,DT1,DT2,DT3,T2,dtime,nl,ig)
	922	IMPLICIT NONE
	923
	924	INTEGER,INTENT(IN) :: nl,ig
	925	REAL,INTENT(IN),DIMENSION(nl) :: T1,DT1,DT2,DT3
	926	REAL*8,INTENT(OUT),DIMENSION(nl) :: T2
	927	REAL,INTENT(IN) :: dtime
	928	INTEGER :: l
	929	DO l=1,nl
	930	T2(l)=T1(l)1D0+(DT1(l)dble(dtime)+ &
	931	& DT2(l)*dble(dtime)+ &
	932	& DT3(l)dble(dtime))1D0
	933	if (T2(l) .ne. T2(l)) then
	934	print*,'Err TMNEW',ig,l,T2(l),T1(l),dT1(l),DT2(l), &
	935	& DT3(l),dtime,dble(dtime)
	936	endif
	937
	938	ENDDO
	939	END
	940
	941	SUBROUTINE QMNEW(Q1,DQ,Q2,dtime,nl,nq,gc,ig)
	942	use chemparam_mod
[1621]	943	use infotrac_phy
[1442]	944	IMPLICIT NONE
	945
	946	INTEGER,INTENT(IN) :: nl,nq
	947	INTEGER,INTENT(IN) :: ig
	948	INTEGER,INTENT(IN),dimension(nq) :: gc
	949	REAL,INTENT(IN),DIMENSION(nl,nqtot) :: Q1,DQ
	950	REAL*8,INTENT(OUT),DIMENSION(nl,nq) :: Q2
	951	REAL,INTENT(IN) :: dtime
	952	INTEGER :: l,iq
	953	DO l=1,nl
	954	DO iq=1,nq
	955	Q2(l,iq)=Q1(l,gc(iq))1D0+(DQ(l,gc(iq))dble(dtime))*1D0
	956	Q2(l,iq)=max(Q2(l,iq),1d-30)
	957	ENDDO
	958	ENDDO
	959	END
	960
	961	SUBROUTINE HSCALE(p,hp,nl)
	962	IMPLICIT NONE
	963
	964	INTEGER :: nl
	965	INTEGER :: l
	966	REAL*8,dimension(nl) :: P
	967	REAL*8,DIMENSION(nl) :: Hp
	968
	969	hp(1)=-log(P(2)/P(1))
	970	hp(nl)=-log(P(nl)/P(nl-1))
	971
	972	DO l=2,nl-1
	973	hp(l)=-log(P(l+1)/P(l-1))
	974	ENDDO
	975	END
	976
	977	SUBROUTINE MMOY(massemoy,mol_tr,qq,gc,nl,nq)
	978	use chemparam_mod
[1621]	979	use infotrac_phy
[1442]	980	IMPLICIT NONE
	981
	982	INTEGER :: nl,nq,l, nn
	983	INTEGER,dimension(nq) :: gc
	984	REAL*8,DIMENSION(nl,nq) :: qq
	985	REAL*8,DIMENSION(nl) :: massemoy
	986	REAL,DIMENSION(nq) :: mol_tr
	987
	988	do l=1,nl
	989	massemoy(l)=1D0/sum(qq(l,:)/dble(mol_tr(:)))
	990	! write(,),'L988 Masse molaire moy: ', massemoy(l)
	991	enddo
	992
	993	END
	994
	995	SUBROUTINE DMMOY(M,H,DM,nl)
	996	IMPLICIT NONE
	997	INTEGER :: nl,l
	998	REAL*8,DIMENSION(nl) :: M,H,DM
	999
	1000	DM(1)=(-3D0M(1)+4D0M(2)-M(3))/2D0/H(1)
	1001	DM(nl)=(3D0M(nl)-4D0M(nl-1)+M(nl-2))/2D0/H(nl)
	1002
	1003	do l=2,nl-1
	1004	DM(l)=(M(l+1)-M(l-1))/H(l)
	1005	enddo
	1006
	1007	END
	1008
	1009	SUBROUTINE ZVERT(P,T,M,Z,nl,ig)
	1010	IMPLICIT NONE
[1591]	1011	#include "YOMCST.h"
[1442]	1012	INTEGER :: nl,l,ig
	1013	REAL*8,dimension(nl) :: P,T,M,Z,H
	1014	REAL*8 :: z0
	1015	REAL*8 :: kbolt,masseU,Konst,g,Hpm
[1591]	1016	masseU=1.e-3/RNAVO
	1017	kbolt=RKBOL
[1442]	1018	Konst=kbolt/masseU
[1591]	1019	g=RG
[1442]	1020
	1021	z0=0d0
	1022	Z(1)=z0
	1023	H(1)=Konst*T(1)/M(1)/g
	1024
	1025	do l=2,nl
	1026	H(l)=Konst*T(l)/M(l)/g
	1027	Hpm=H(l-1)
	1028	Z(l)=z(l-1)-Hpm*log(P(l)/P(l-1))
	1029	if (Z(l) .ne. Z(l)) then
	1030	print*,'pb',l,ig
	1031	print*,'P',P
	1032	print*,'T',T
	1033	print*,'M',M
	1034	print*,'Z',Z
	1035	print*,'Hpm',Hpm
	1036	endif
	1037	enddo
	1038
	1039	END
	1040
	1041	SUBROUTINE RHOTOT(P,T,M,qq,rhoN,rhoK,nl,nq)
	1042	IMPLICIT NONE
[1591]	1043	#include "YOMCST.h"
[1442]	1044	REAL*8 :: kbolt,masseU,Konst
	1045	INTEGER :: nl,nq,l,iq
	1046	REAL*8,DIMENSION(nl) :: P,T,M,RHON
	1047	REAL*8,DIMENSION(nl,nq) :: RHOK,qq
[1591]	1048	masseU=1.e-3/RNAVO
	1049	kbolt=RKBOL
[1442]	1050	Konst=Kbolt/masseU
	1051
	1052	do l=1,nl
	1053	RHON(l)=P(l)*M(l)/T(l)/Konst
	1054	do iq=1,nq
	1055	RHOK(l,iq)=qq(l,iq)*RHON(l)
	1056	enddo
	1057	enddo
	1058
	1059	END
	1060
[1609]	1061	SUBROUTINE UPPER_RESOL(P,T,Z,M,RR,Rk, &
[1442]	1062	& qq,mol_tr,gc,Praf,Traf,Qraf,Mraf,Zraf, &
	1063	& Nraf,Nrafk,Rraf,Rrafk,il,nl,nq,nlx,ig)
	1064	use chemparam_mod
[1621]	1065	use infotrac_phy
[1442]	1066	IMPLICIT NONE
[1591]	1067	#include "YOMCST.h"
[1442]	1068	INTEGER :: nl,nq,il,l,i,iq,nlx,iz,ig
	1069	INTEGER :: gc(nq)
	1070	INTEGER,DIMENSION(1) :: indz
[1609]	1071	REAL*8, DIMENSION(nlx) :: P,T,Z,M,RR
[1442]	1072	REAL*8, DIMENSION(nlx,nq) :: qq,Rk
	1073	REAL*8, DIMENSION(nl) :: Praf,Traf,Mraf,Zraf,Nraf,Rraf
	1074	REAL*8 :: kbolt,masseU,Konst,g
	1075	REAL*8, DIMENSION(nl,nq) :: Qraf,Rrafk,Nrafk
	1076	REAL*8 :: facZ,dZ,H
	1077	REAL,DIMENSION(nq) :: mol_tr
[1591]	1078	masseU=1.e-3/RNAVO
	1079	kbolt=RKBOL
[1442]	1080	Konst=Kbolt/masseU
[1591]	1081	g=RG
[1442]	1082
	1083
	1084	Zraf(1)=z(il)
	1085	Praf(1)=P(il)
	1086	Traf(1)=T(il)
	1087	Nraf(1)=Praf(1)/kbolt/Traf(1)
	1088	do iq=1,nq
	1089	Qraf(1,iq)=qq(il,iq)
	1090	enddo
	1091	Mraf(1)=1d0/sum(Qraf(1,:)/dble(mol_tr(:)))
	1092	Rraf(1)=Mraf(1)masseUNraf(1)
	1093	do iq=1,nq
	1094	Rrafk(1,iq)=Rraf(1)*Qraf(1,iq)
	1095	Nrafk(1,iq)=Rrafk(1,iq)/masseU/dble(mol_tr(iq))
	1096	enddo
	1097	Zraf(nl)=z(nlx)
	1098
	1099	do l=2,nl-1
	1100	Zraf(l)=Zraf(1)+(Zraf(nl)-Zraf(1))/dble(nl-1)*dble((l-1))
	1101	indz=maxloc(z,mask=z <= Zraf(l))
	1102	iz=indz(1)
	1103	if (iz .lt. 1 .or. iz .gt. nlx) then
	1104	print*,'danger',iz,nl,Zraf(l),l,Zraf(1),Zraf(nl),z,P,T,ig
	1105	stop
	1106	endif
	1107	dZ=Zraf(l)-Zraf(l-1)
	1108	! dZ=Zraf(l)-z(iz)
	1109	facz=(Zraf(l)-z(iz))/(z(iz+1)-z(iz))
	1110	Traf(l)=T(iz)+(T(iz+1)-T(iz))*facz
	1111	do iq=1,nq
	1112	! Qraf(l,iq)=qq(iz,iq)+(qq(iz+1,iq)-qq(iz,iq))*facz
	1113	Rrafk(l,iq)=Rk(iz,iq)+(Rk(iz+1,iq)-Rk(iz,iq))*facZ
	1114	Rrafk(l,iq)=Rk(iz,iq)(Rk(iz+1,iq)/Rk(iz,iq))*facZ
	1115	enddo
	1116	! Mraf(l)=1D0/(sum(qraf(l,:)/dble(mol_tr(:))))
	1117	Rraf(l)=sum(Rrafk(l,:))
	1118	do iq=1,nq
	1119	Qraf(l,iq)=Rrafk(l,iq)/Rraf(l)
	1120	enddo
	1121	Mraf(l)=1D0/(sum(qraf(l,:)/dble(mol_tr(:))))
	1122	! H=Konst*Traf(l)/Mraf(l)/g
	1123	! H=Konst*T(iz)/M(iz)/g
	1124	! Praf(l)=P(iz)*exp(-dZ/H)
	1125	! Praf(l)=Praf(l-1)*exp(-dZ/H)
	1126	! print,'iz',l,iz,Praf(il-1)exp(-dZ/H),z(iz),z(iz+1),H
	1127	Nraf(l)=Rraf(l)/Mraf(l)/masseU
	1128	Praf(l)=Nraf(l)kboltTraf(l)
	1129	! Rraf(l)=Nraf(l)Mraf(l)masseU
	1130	do iq=1,nq
	1131	! Rrafk(l,iq)=Rraf(l)*Qraf(l,iq)
	1132	Nrafk(l,iq)=Rrafk(l,iq)/dble(mol_tr(iq))/masseU
	1133	if (Nrafk(l,iq) .lt. 0. .or. &
	1134	& Nrafk(l,iq) .ne. Nrafk(l,iq)) then
	1135	print*,'pb interpolation',l,iq,Nrafk(l,iq),Rrafk(l,iq), &
	1136	& Qraf(l,iq),Rk(iz,iq),Rk(iz+1,iq),facZ,Zraf(l),z(iz)
	1137	stop
	1138	endif
	1139	enddo
	1140	enddo
	1141	Zraf(nl)=z(nlx)
	1142	Traf(nl)=T(nlx)
	1143	do iq=1,nq
	1144	Rrafk(nl,iq)=Rk(nlx,iq)
[1609]	1145	Qraf(nl,iq)=Rk(nlx,iq)/RR(nlx)
[1442]	1146	Nrafk(nl,iq)=Rk(nlx,iq)/dble(mol_tr(iq))/masseU
	1147	enddo
	1148	Nraf(nl)=sum(Nrafk(nl,:))
	1149	Praf(nl)=Nraf(nl)kboltTraf(nl)
	1150	Mraf(nl)=1D0/sum(Qraf(nl,:)/dble(mol_tr(:)))
	1151	END
	1152
	1153	SUBROUTINE CORRMASS(qq,qint,FacMass,nl,nq)
	1154	IMPLICIT NONE
	1155	INTEGER :: nl,nq,l,nn
	1156	REAL*8,DIMENSION(nl,nq) :: qq,qint,FacMass
	1157
	1158	do nn=1,nq
	1159	do l=1,nl
	1160	FacMass(l,nn)=qq(l,nn)/qint(l,nn)
	1161	enddo
	1162	enddo
	1163
	1164	END
	1165
	1166
	1167	SUBROUTINE DCOEFF(nn,Dij,P,T,N,Nk,D,nl,nq)
	1168	IMPLICIT NONE
	1169	REAL*8,DIMENSION(nl) :: N,T,D,P
	1170	REAL*8,DIMENSION(nl,nq) :: Nk
	1171	INTEGER :: nn,nl,nq,l,iq
	1172	REAL,DIMENSION(nq,nq) :: Dij
	1173	REAL*8 :: interm,P0,T0,ptfac,dfac
	1174
	1175	P0=1D5
	1176	T0=273d0
	1177
	1178
	1179	do l=1,nl
	1180	ptfac=(P0/P(l))(T(l)/T0)*1.75d0
	1181	D(l)=0d0
	1182	interm=0d0
	1183	do iq=1,nq
	1184	if (iq .ne. nn) then
	1185	dfac=dble(dij(nn,iq))*ptfac
	1186	interm=interm+Nk(l,iq)/N(l)/dfac
	1187	endif
	1188	enddo
[2791]	1189	D(l)=(1d0-Nk(l,nn)/N(l))/interm
[1442]	1190	enddo
	1191	END
	1192
	1193	SUBROUTINE HSCALEREAL(nn,Nk,Dz,H,nl,nq)
	1194	IMPLICIT NONE
	1195	INTEGER :: nn,nl,nq,l
	1196	REAL*8,DIMENSION(nl) :: H
	1197	REAL*8,DIMENSION(nl,nq) :: Nk
	1198	REAL*8 :: Dz
	1199
	1200	H(1)=(-3D0Nk(1,nn)+4d0NK(2,nn)-Nk(3,nn))/(2D0*DZ)/ &
	1201	& NK(1,nn)
	1202
	1203	H(1)=-1D0/H(1)
	1204
	1205	DO l=2,nl-1
	1206	H(l)=(Nk(l+1,nn)-NK(l-1,nn))/(2D0*DZ)/NK(l,nn)
	1207	H(l)=-1D0/H(l)
	1208	ENDDO
	1209
	1210	H(nl)=(3D0Nk(nl,nn)-4D0Nk(nl-1,nn)+Nk(nl-2,nn))/(2D0*DZ)/ &
	1211	& Nk(nl,nn)
	1212	H(nl)=-1D0/H(nl)
	1213
	1214	! do l=1,nl
	1215	! if (abs(H(l)) .lt. 100.) then
	1216	! print*,'H',l,H(l),Nk(l,nn),nn
	1217	! endif
	1218	! enddo
	1219
	1220	END
	1221
	1222	SUBROUTINE VELVERT(nn,T,H,D,Dz,masse,W,nl)
	1223	IMPLICIT NONE
[1591]	1224	#include "YOMCST.h"
[1442]	1225	INTEGER :: l,nl,nn
	1226	REAL*8,DIMENSION(nl) :: T,H,D,W,DT
	1227	REAL*8 :: Dz,Hmol,masse
	1228	REAL*8 :: kbolt,masseU,Konst,g
[1591]	1229	masseU=1.e-3/RNAVO
	1230	kbolt=RKBOL
[1442]	1231	Konst=Kbolt/masseU
[1591]	1232	g=RG
[1442]	1233
	1234	DT(1)=1D0/T(1)(-3D0T(1)+4D0T(2)-T(3))/(2D0DZ)
	1235	Hmol=Konst*T(1)/masse/g
	1236	W(1)=-D(1)*(1D0/H(1)-1D0/Hmol-DT(1))
	1237
	1238	DO l=2,nl-1
	1239	DT(l)=1D0/T(l)(T(l+1)-T(l-1))/(2D0DZ)
	1240	Hmol=Konst*T(l)/masse/g
	1241	W(l)=-D(l)*(1D0/H(l)-1D0/Hmol-DT(l))
	1242	ENDDO
	1243
	1244	DT(nl)=1D0/T(nl)(3d0T(nl)-4D0T(nl-1)+T(nl-2))/(2D0DZ)
	1245	Hmol=Konst*T(nl)/masse/g
	1246	W(nl)=-D(nl)*(1D0/H(nl)-1D0/Hmol-DT(nl))
	1247
	1248	! do l=1,nl
	1249	! print*,'W',W(l),D(l),H(l),DT(l)
	1250	! enddo
	1251
	1252	END
	1253
	1254	SUBROUTINE TIMEDIFF(nn,H,W,TIME,nl)
	1255	IMPLICIT NONE
	1256	INTEGER :: nn,nl,l
	1257	REAL*8,DIMENSION(nl) :: W,H,TIME
	1258
	1259	DO l=1,nl
	1260	TIME(l)=abs(H(l)/W(l))
	1261	if (TIME(l) .lt. 1.D-4) then
	1262	! print*,'low tdiff',H(l),W(l),nn,l
	1263	endif
	1264	ENDDO
	1265
	1266	END
	1267
	1268
	1269	SUBROUTINE DIFFPARAM(T,D,dz,alpha,beta,gama,delta,eps,nl)
	1270	IMPLICIT NONE
	1271	INTEGER :: nl,l
	1272	REAL*8,DIMENSION(nl) :: T,D
	1273	REAL*8 :: DZ,DZinv
	1274	REAL*8,DIMENSION(nl) :: alpha,beta,gama,delta,eps
	1275
	1276	! Compute alpha,beta and delta
	1277	! lower altitude values
	1278	DZinv=1d0/(2D0*DZ)
	1279
	1280	beta(1)=1d0/T(1)
	1281	alpha(1)=beta(1)(-3D0T(1)+4D0T(2)-T(3))Dzinv
	1282	delta(1)=(-3D0D(1)+4D0D(2)-D(3))*Dzinv
	1283
	1284	beta(2)=1d0/T(2)
	1285	alpha(2)=beta(2)(T(3)-T(1))Dzinv
	1286	delta(2)=(D(3)-D(1))*Dzinv
	1287
	1288	! do l=2,nl-1
	1289	! beta(l)=1D0/T(l)
	1290	! alpha(l)=beta(l)(T(l+1)-T(l-1))Dzinv
	1291	! delta(l)=(D(l+1)-D(l-1))*Dzinv
	1292	! end do
	1293
	1294	do l=3,nl-1
	1295	beta(l)=1D0/T(l)
	1296	alpha(l)=beta(l)(T(l+1)-T(l-1))Dzinv
	1297	delta(l)=(D(l+1)-D(l-1))*Dzinv
	1298	gama(l-1)=(beta(l)-beta(l-2))*Dzinv
	1299	eps(l-1)=(alpha(l)-alpha(l-2))*Dzinv
	1300	enddo
	1301
	1302	! Upper altitude values
	1303
	1304	beta(nl)=1D0/T(nl)
	1305	alpha(nl)=beta(nl)(3D0T(nl)-4D0T(nl-1)+T(nl-2))Dzinv
	1306	delta(nl)=(3D0D(nl)-4D0D(nl-1)+D(nl-2))*Dzinv
	1307
	1308	! Compute the gama and eps coefficients
	1309	! Lower altitude values
	1310
	1311	gama(1)=(-3D0beta(1)+4D0beta(2)-beta(3))*Dzinv
	1312	eps(1)=(-3D0alpha(1)+4D0alpha(2)-alpha(3))*Dzinv
	1313
	1314	gama(nl-1)=(beta(nl)-beta(nl-2))*Dzinv
	1315	eps(nl-1)=(alpha(nl)-alpha(nl-2))*Dzinv
	1316
	1317	! do l=2,nl-1
	1318	! gama(l)=(beta(l+1)-beta(l-1))*Dzinv
	1319	! eps(l)=(alpha(l+1)-alpha(l-1))*Dzinv
	1320	! end do
	1321
	1322	gama(nl)=(3D0beta(nl)-4D0beta(nl-1)+beta(nl-2))*Dzinv
	1323	eps(nl)=(3D0alpha(nl)-4D0alpha(nl-1)+alpha(nl-2))*Dzinv
	1324
	1325	! do l=1,nl
	1326	! print*,'test diffparam',alpha(l),beta(l),delta(l),gama(l),eps(l)
	1327	! enddo
	1328	! stop
	1329
	1330	END
	1331
	1332
	1333	SUBROUTINE MATCOEFF(alpha,beta,gama,delta,eps,Dad,rhoad, &
	1334	& FacEsc,dz,dt,A,B,C,D,nl)
	1335	IMPLICIT NONE
	1336	INTEGER :: nl,l
	1337	REAL*8, DIMENSION(nl) :: alpha,beta,gama,delta,eps,Dad,RHoad
	1338	REAL*8 :: dz,dt,del1,del2,del3,FacEsc
	1339	REAL*8, DIMENSION(nl) :: A,B,C,D
	1340	del1=dt/2d0/dz
	1341	del2=dt/dz/dz
	1342	del3=dt
	1343
	1344	! lower boundary coefficients no change
	1345	A(1)=0d0
	1346	B(1)=1d0
	1347	C(1)=0d0
	1348	D(1)=rhoAd(1)
	1349
	1350	do l=2,nl-1
	1351	A(l)=(delta(l)+(alpha(l)+beta(l))Dad(l))del1-Dad(l)*del2
	1352	B(l)=-(delta(l)(alpha(l)+beta(l))+Dad(l)(gama(l)+eps(l)))*del3
	1353	B(l)=B(l)+1d0+2d0Dad(l)del2
	1354	C(l)=-(delta(l)+(alpha(l)+beta(l))Dad(l))del1-Dad(l)*del2
	1355	D(l)=rhoAD(l)
	1356	enddo
	1357
	1358
	1359	! Upper boundary coefficients Diffusion profile
	1360	C(nl)=0d0
	1361	B(nl)=-1d0
	1362	A(nl)=exp(-dZ(alpha(nl)+beta(nl)))FacEsc
	1363	D(nl)=0D0
	1364
	1365
	1366	END
	1367
	1368	SUBROUTINE Checkmass(X,Y,nl,nn)
	1369	IMPLICIT NONE
	1370
	1371	INTEGER :: nl,nn
	1372	REAL*8,DIMENSION(nl) :: X,Y
	1373	REAL*8 Xtot,Ytot
	1374
	1375	Xtot=sum(X)
	1376	Ytot=sum(Y)
	1377
	1378	if (abs((Xtot-Ytot)/Xtot) .gt. 1d-4) then
	1379	print*,'no conservation for mass',Xtot,Ytot,nn
	1380	endif
	1381	END
	1382
	1383	SUBROUTINE Checkmass2(qold,qnew,P,il,nl,nn,nq)
	1384	IMPLICIT NONE
[1591]	1385	#include "YOMCST.h"
[1442]	1386	INTEGER :: nl,nn,l,nq,il
	1387	REAL,DIMENSION(nl+1) :: P
	1388	REAL*8,DIMENSION(nl,nq) :: qold,qnew
	1389	REAL*8 :: DM,Mold,Mnew,g
[1591]	1390	g=RG
[1442]	1391	DM=0d0
	1392	Mold=0d0
	1393	Mnew=0d0
	1394	DO l=il,nl
	1395	DM=DM+(qnew(l,nn)-qold(l,nn))*(dble(P(l))-dble(P(l+1)))/g
	1396	Mold=Mold+qold(l,nn)*(dble(P(l))-dble(P(l+1)))/g
	1397	Mnew=Mnew+qnew(l,nn)*(dble(P(l))-dble(P(l+1)))/g
	1398	! print*,'l',l,qold(l,nn),qnew(l,nn),Mold,Mnew,DM,P(l),P(l+1)
	1399	ENDDO
	1400	IF (abs(DM/Mold) .gt. 1d-5) THEN
[1591]	1401	Print*,'We dont conserve mass',nn,DM,Mold,Mnew,DM/Mold
[1442]	1402	ENDIF
	1403
	1404	END
	1405
	1406	SUBROUTINE GCMGRID_P(Z,P,Q,T,Nk,Rk,qq,qnew,tt,tnew, &
	1407	& pp,M,gc,nl,nq,nlx,ig)
	1408	use chemparam_mod
[1621]	1409	use infotrac_phy
[1442]	1410	IMPLICIT NONE
[1591]	1411	#include "YOMCST.h"
[1442]	1412	INTEGER :: nl,nq,nlx,il,nn,iP,ig,compteur
	1413	INTEGER,DIMENSION(1) :: indP
	1414	INTEGER,DIMENSION(nq) :: gc
	1415	REAL*8,DIMENSION(nl) :: Z,P,T
	1416	REAL*8,DIMENSION(nl,nq) :: Q,Nk,Rk
	1417	REAL,DIMENSION(nq) :: M
	1418	REAL*8,DIMENSION(nq) :: nNew
	1419	REAL*8,DIMENSION(nlx) :: pp,tt,tnew
	1420	REAL*8,DIMENSION(nlx) :: rhonew
	1421	REAL*8,DIMENSION(nlx,nq) :: qq,qnew,rhoknew
	1422	REAL*8 :: kbolt,masseU,Konst,g,Dz,facP,Hi
	1423	REAL*8 :: Znew,Znew2,Pnew,Pnew2
[1591]	1424	masseU=1.e-3/RNAVO
	1425	kbolt=RKBOL
[1442]	1426	Konst=Kbolt/masseU
[1591]	1427	g=RG
[1442]	1428	Dz=Z(2)-Z(1)
	1429	Znew=Z(nl)
	1430	Znew2=Znew+dz
	1431	! print*,'dz',Znew,Znew2,dz
	1432	nNew(1:nq)=Nk(nl,1:nq)
	1433	Pnew=P(nl)
	1434
	1435	do il=1,nlx
	1436	! print*,'il',il,pp(il),P(1),P(nl)
	1437	if (pp(il) .ge. P(1)) then
	1438	qnew(il,:)=qq(il,:)
	1439	tnew(il)=tt(il)
	1440	endif
	1441	if (pp(il) .lt. P(1)) then
	1442	if (pp(il) .gt. P(nl)) then
	1443	indP=maxloc(P,mask=P < pp(il))
	1444	iP=indP(1)-1
	1445	if (iP .lt. 1 .or. iP .gt. nl) then
	1446	print*,'danger 2',iP,nl,pp(il)
	1447	endif
	1448	facP=(pp(il)-P(ip))/(P(ip+1)-P(ip))
	1449	! print*,'P',P(ip),P(ip+1),facP,indP,iP
	1450
	1451	! do nn=1,nq
	1452	! qnew(il,nn)=Q(iP,nn)+
	1453	! & (Q(ip+1,nn)-Q(ip,nn))*facP
	1454	! enddo
	1455
	1456	do nn=1,nq
	1457	rhoknew(il,nn)=Rk(iP,nn)+ &
	1458	& (Rk(ip+1,nn)-Rk(ip,nn))*facP
	1459	enddo
	1460	tnew(il)=T(iP)+(T(iP+1)-T(iP))*facP
	1461	rhonew(il)=sum(rhoknew(il,:))
	1462	do nn=1,nq
	1463	qnew(il,nn)=rhoknew(il,nn)/rhonew(il)
	1464	enddo
	1465
	1466	else ! pp < P(nl) need to extrapolate density of each specie
	1467	Pnew2=Pnew
	1468
	1469	compteur=0
[1675]	1470	do while (Pnew2 .ge. pp(il))
[1442]	1471	compteur=compteur+1
	1472	do nn=1,nq
	1473	Hi=Konst*T(nl)/dble(M(nn))/g
	1474	Nnew(nn)=Nnew(nn)*exp(-dZ/Hi)
	1475	enddo
	1476	Pnew=Pnew2
	1477	Pnew2=kboltT(nl)sum(Nnew(:))
	1478	Znew=Znew2
	1479	Znew2=Znew2+Dz
	1480	if (compteur .ge. 100000) then
	1481	print*,'error moldiff_red infinite loop'
[1675]	1482	print*,ig,il,pp(il),tt(nl),Pnew2,qnew(il,:),Znew2
[1442]	1483	stop
	1484	endif
	1485	! print*,'test',Pnew2,Znew2,Nnew(nq),pp(il)
	1486	enddo
	1487
	1488	facP=(pp(il)-Pnew)/(Pnew2-Pnew)
	1489
	1490	! do nn=1,nq
	1491	! qnew(il,nn)=dble(M(nn))*Nnew(nn)
	1492	! & /sum(dble(M(:))*Nnew(:))
	1493	! enddo
	1494
	1495	do nn=1,nq
	1496	rhoknew(il,nn)=dble(M(nn))*Nnew(nn)
	1497	enddo
	1498	rhonew(il)=sum(rhoknew(il,:))
	1499	do nn=1,nq
	1500	qnew(il,nn)=rhoknew(il,nn)/rhonew(il)
	1501	enddo
	1502	tnew(il)=T(nl)
	1503	endif
	1504	endif
	1505	enddo
	1506
	1507	END
	1508
	1509	SUBROUTINE GCMGRID_P2(Z,P,Q,T,Nk,Rk,qq,qnew,tt,tnew &
	1510	& ,pp,M,gc,nl,nq,nlx,facM,ig)
	1511	! use chemparam_mod
	1512	! use infotrac
	1513	IMPLICIT NONE
[1591]	1514	#include "YOMCST.h"
[1442]	1515	INTEGER :: nl,nq,nlx,il,nn,iP,ig,compteur
	1516	INTEGER,DIMENSION(1) :: indP
	1517	INTEGER,DIMENSION(nq) :: gc
	1518	REAL*8,DIMENSION(nl) :: Z,P,T
	1519	REAL*8,DIMENSION(nl,nq) :: Q,Nk,Rk
	1520	REAL,DIMENSION(nq) :: M
	1521	REAL*8,DIMENSION(nq) :: nNew
	1522	REAL*8,DIMENSION(nlx) :: pp,rhonew,tt,tnew
	1523	REAL*8,DIMENSION(nlx,nq) :: qq,qnew,facM,rhoknew
	1524	REAL*8 :: kbolt,masseU,Konst,g,Dz,facP,Hi
	1525	REAL*8 :: Znew,Znew2,Pnew,Pnew2
[1591]	1526	masseU=1.e-3/RNAVO
	1527	kbolt=RKBOL
[1442]	1528	Konst=Kbolt/masseU
[1591]	1529	g=RG
[1442]	1530	Dz=Z(2)-Z(1)
	1531	Znew=Z(nl)
	1532	Znew2=Znew+dz
	1533	! print*,'dz',Znew,Znew2,dz
	1534	nNew(1:nq)=Nk(nl,1:nq)
	1535	Pnew=P(nl)
	1536
	1537	do il=1,nlx
	1538	! print*,'il',il,pp(il),P(1),P(nl)
	1539	if (pp(il) .ge. P(1)) then
	1540	qnew(il,:)=qq(il,:)
	1541	tnew(il)=tt(il)
	1542	endif
	1543	if (pp(il) .lt. P(1)) then
	1544	if (pp(il) .gt. P(nl)) then
	1545	indP=maxloc(P,mask=P < pp(il))
	1546	iP=indP(1)-1
	1547	if (iP .lt. 1 .or. iP .gt. nl) then
	1548	print*,'danger 3',iP,nl,pp(il)
	1549	endif
	1550	facP=(pp(il)-P(ip))/(P(ip+1)-P(ip))
	1551	! print*,'P',P(ip),P(ip+1),facP,indP,iP
	1552
	1553	! do nn=1,nq
	1554	! qnew(il,nn)=Q(iP,nn)+
	1555	! & (Q(ip+1,nn)-Q(ip,nn))*facP
	1556	! enddo
	1557
	1558	do nn=1,nq
	1559	rhoknew(il,nn)=(RK(iP,nn)+ &
	1560	& (RK(iP+1,nn)-Rk(iP,nn))facP)facM(il,nn)
	1561	enddo
	1562	tnew(il)=T(iP)+(T(ip+1)-T(iP))*facP
	1563	rhonew(il)=sum(rhoknew(il,:))
	1564	do nn=1,nq
	1565	qnew(il,nn)=rhoknew(il,nn)/rhonew(il)
	1566	enddo
	1567
	1568	else ! pp < P(nl) need to extrapolate density of each specie
	1569	Pnew2=Pnew
	1570
	1571	compteur=0
[1675]	1572	do while (Pnew2 .ge. pp(il))
[1442]	1573	compteur=compteur+1
	1574	do nn=1,nq
	1575	Hi=Konst*T(nl)/dble(M(nn))/g
	1576	Nnew(nn)=Nnew(nn)*exp(-dZ/Hi)
	1577	enddo
	1578	Pnew=Pnew2
	1579	Pnew2=kboltT(nl)sum(Nnew(:))
	1580	Znew=Znew2
	1581	Znew2=Znew2+Dz
	1582	if (compteur .ge. 100000) then
	1583	print*,'pb moldiff_red infinite loop'
[1675]	1584	print*,ig,nl,T(nl),Pnew2,qnew(il,:),Znew2
[1442]	1585	stop
	1586	endif
	1587
	1588	! print*,'test',Pnew2,Znew2,Nnew(nq),pp(il)
	1589	enddo
	1590
	1591	facP=(pp(il)-Pnew)/(Pnew2-Pnew)
	1592
	1593	! do nn=1,nq
	1594	! qnew(il,nn)=dble(M(nn))*Nnew(nn)
	1595	! & /sum(dble(M(:))*Nnew(:))
	1596	! enddo
	1597
	1598	do nn=1,nq
	1599	rhoknew(il,nn)=dble(M(nn))Nnew(nn)FacM(il,nn)
	1600	enddo
	1601	rhonew(il)=sum(rhoknew(il,:))
	1602	do nn=1,nq
	1603	qnew(il,nn)=rhoknew(il,nn)/rhonew(il)
	1604	enddo
	1605	tnew(il)=T(nl)
	1606
	1607	endif
	1608	endif
	1609	enddo
	1610
	1611	! write(,), ' -- Sortie moldiff_red -- '
	1612
	1613	END
	1614

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

Download in other formats:

Original Format