Context Navigation

moldiff_red.F90 @ 2203

Last change on this file since 2203 was 1675, checked in by slebonnois, 8 years ago
SL: update of xml file for xios output + a few modifications for run with uppr atmosphere
File size: 39.1 KB

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

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

Download in other formats:

Original Format