Context Navigation

moldiff_red.F90 @ 815

Last change on this file since 815 was 710, checked in by emillour, 13 years ago

Mars GCM:

Decreased time step for molecular diffusion (diffusion.h) for better stability
Added test in moldiff_red.F90 to enforce that the system can't be stuck in an infinite loop.
The creation and output of coefficients in file 'coeffs.dat' in moldiff_red.F is now controled by an internal flag (by default set to false).

JYC + EM

File size: 37.5 KB

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

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

Download in other formats:

Original Format