Context Navigation

← Previous Revision
Next Revision →
Normal
Revision Log

moldiff_red.F90

Last change on this file was 3015, checked in by emillour, 17 months ago
Mars PCM: Code cleanup in diffusion, turn variables from diffusion.h into module variables in moldiff_red.F90 (turn it into a module in the process). Also turn moldiffcoeff_red.F and thermosphere.F into modules. EM
File size: 41.3 KB

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

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

Download in other formats:

Original Format