Context Navigation

moldiff_red.F90 @ 3567

Last change on this file since 3567 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

Line
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
9	subroutine moldiff_red(ngrid,nlayer,nq,pplay,pplev,pt,pdt,pq,pdq,&
10	ptimestep,zzlay,pdteuv,pdtconduc,pdqdiff,&
11	PhiEscH,PhiEscH2,PhiEscD)
12
13	use tracer_mod, only: noms, mmol
14	use geometry_mod, only: cell_area
15	use planetwide_mod, only: planetwide_sumval
16	USE mod_phys_lmdz_para, only: is_master,bcast
17	use moldiffcoeff_red_mod, only: moldiffcoeff_red
18
19	implicit none
20
21	! July 2014 JYC ADD BALISTIC Transport coupling to compute wup for H and H2
22
23
24
25	!
26	! Input/Output
27	!
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
31	real ptimestep
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)
42	real*8 PhiEscH,PhiEscH2,PhiEscD
43	!
44	! Local
45	!
46
47
48	! real hco2(ncompdiff),ho
49
50	integer,dimension(nq) :: indic_diff
51	integer ig,iq,nz,l,k,n,nn,p,ij0
52	integer istep,il,gcn,ntime,nlraf
53	real*8 masse
54	real*8 masseU,kBolt,RgazP,Rmars,Grav,Mmars
55	real*8 rho0,D0,T0,H0,time0,dZ,time,dZraf,tdiff,Zmin,Zmax
56	real*8 FacEsc,invsgmu
57	real*8 PhiauxH(ngrid),PhiauxH2(ngrid),PhiauxD(ngrid)
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)
63	real*8,dimension(:,:),allocatable,save :: qq,qnew,qint,FacMass
64	real*8,dimension(:,:),allocatable,save :: rhoK,rhokinit
65	real*8 rhoT(nlayer)
66	real*8 dmmeandz(nlayer)
67	real*8 massemoy(nlayer)
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
76	real*8,dimension(:),allocatable,save :: wi,Wad,Uthermal,Lambdaexo,Hspecie
77	real*8,dimension(:),allocatable,save :: Mtot1,Mtot2,Mraf1,Mraf2
78	integer,parameter :: ListeDiffNb=16
79	character(len=20),dimension(ListeDiffNb) :: ListeDiff
80	real*8,parameter :: pi=3.141592653
81	real*8,parameter :: g=3.72d0
82	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
83	! tracer numbering in the molecular diffusion
84	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
85	! We need the index of escaping species
86
87	integer,save :: i_h2
88	integer,save :: i_h
89	integer,save :: i_d
90	! vertical index limit for the molecular diffusion
91	integer,save :: il0
92
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
95	! Tracer indexes in the GCM:
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
114	! integer,save :: g_oplus=0
115	! integer,save :: g_hplus=0
116
117	integer,save :: ncompdiff
118	integer,dimension(:),allocatable,save :: gcmind ! array of GCM indexes
119	integer ierr,compteur
120
121	logical,save :: firstcall=.true.
122	! real abfac(ncompdiff)
123	real,dimension(:,:),allocatable,save :: dij
124	real,save :: step
125
126	!$OMP THREADPRIVATE(ncompdiff,gcmind,firstcall,dij,step)
127
128	! Initializations at first call
129	if (firstcall) then
130
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'
147	ListeDiff(15)='he'
148	ListeDiff(16)='hdo_vap'
149	i_h=1000
150	i_h2=1000
151	i_d=1000
152	! On regarde quelles especes diffusables sont presentes
153
154	ncompdiff=0
155	indic_diff(1:nq)=0
156
157	do nn=1,nq
158	do n=1,ListeDiffNb
159	if (ListeDiff(n) .eq. noms(nn)) then
160	indic_diff(nn)=1
161	! if (noms(nn) .eq. 'h') i_h=n
162	! if (noms(nn) .eq. 'h2') i_h2=n
163	endif
164
165	enddo
166	if (indic_diff(nn) .eq. 1) then
167	print*,'specie ', noms(nn), 'diffused in moldiff_red'
168	ncompdiff=ncompdiff+1
169	endif
170	enddo
171	print*,'number of diffused species:',ncompdiff
172	allocate(gcmind(ncompdiff))
173
174	! Store gcm indexes in gcmind
175	n=0
176	do nn=1,nq
177	if (indic_diff(nn) .eq. 1) then
178	n=n+1
179	gcmind(n)=nn
180	if (noms(nn) .eq. 'h') i_h=n
181	if (noms(nn) .eq. 'h2') i_h2=n
182	if (noms(nn) .eq. 'd') i_d=n
183	endif
184	enddo
185
186	! print*,'gcmind',gcmind,i_h,i_h2
187
188	! find vertical index above which diffusion is computed
189	if(is_master) then
190	do l=1,nlayer
191	if (pplay(1,l) .gt. Pdiff) then
192	il0=l
193	endif
194	enddo
195	il0=il0+1
196	endif !(is_master)
197	CALL bcast(il0)
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
205	! allocatation des tableaux dependants du nombre d especes diffusees
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))
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
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
233	! PI =3.141592653
234	! g=3.72d0
235	Rmars=3390000d0 ! Used to compute escape parameter no need to be more accurate
236	Grav=6.67d-11
237	Mmars=6.4d23
238	ij0=6000 ! For test
239	invsgmu=1d0/g/masseU
240
241	! Compute the wup(ig) for H and H2 using the balistic code from R Yelle
242
243	PhiEscH=0D0
244	PhiEscH2=0D0
245	PhiEscD=0D0
246
247	do ig=1,ngrid
248	pp=dble(pplay(ig,:))
249
250	! Update the temperature
251
252	! CALL TMNEW(pt(ig,:),pdt(ig,:),pdtconduc(ig,:),pdteuv(ig,:) &
253	! & ,tt,ptimestep,nlayer,ig)
254	do l=1,nlayer
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
263	enddo ! of do l=1,nlayer
264
265	! Update the mass mixing ratios modified by other processes
266
267	! CALL QMNEW(pq(ig,:,:),pdq(ig,:,:),qq,ptimestep,nlayer, &
268	! & ncompdiff,gcmind,ig)
269	do iq=1,ncompdiff
270	do l=1,nlayer
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)
274	enddo ! of do l=1,nlayer
275	enddo ! of do iq=1,ncompdiff
276
277	! Compute the Pressure scale height
278
279	CALL HSCALE(pp,hp,nlayer)
280
281	! Compute the atmospheric mass (in Dalton)
282
283	CALL MMOY(massemoy,mmol,qq,gcmind,nlayer,ncompdiff)
284
285	! Compute the vertical gradient of atmospheric mass
286
287	CALL DMMOY(massemoy,hp,dmmeandz,nlayer)
288
289	! Compute the altitude of each layer
290
291	CALL ZVERT(pp,tt,massemoy,zz,nlayer,ig)
292
293	! Compute the total mass density (kg/m3)
294
295	CALL RHOTOT(pp,tt,massemoy,qq,RHOT,RHOK,nlayer,ncompdiff)
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
306	do l=il0,nlayer
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
313	Zmin=zz(il0)
314	Zmax=zz(nlayer)
315
316
317	! If Zmax > 4000 km there is a problem / stop
318
319	if (Zmax .gt. 4000000.) then
320	Print*,'Zmax too high',ig,zmax,zmin
321	do l=1,nlayer
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
329	! and fixed by the resolution (dzres) specified as module variable
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
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))
349
350	! before beginning, I use a better vertical resolution above il0,
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, &
356	& Nraf,Nrafk,Rraf,Rrafk,il0,nlraf,ncompdiff,nlayer,ig)
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 &
365	& ,pp,mmol,gcmind,nlraf,ncompdiff,nlayer,ig)
366
367	! We compute the mass correction factor of each specie at each pressure level
368
369	CALL CORRMASS(qq,qint,FacMass,nlayer,ncompdiff)
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
402	! do l=1,nlayer
403	! print*,'l',l,zz(l),pp(l),tt(l),sum(qq(l,:)),massemoy(l)
404	! enddo
405
406	! The diffusion is computed above il0 computed from Pdiff
407	! (defined above as module variable)
408	! No change below il0
409
410	do l=1,nlayer
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
428	if (nn .eq. i_h .or. nn .eq. i_h2 .or. nn .eq. i_d) then
429	wi(nn)=Uthermal(nn)/2d0/sqrt(PI)exp(-Lambdaexo(nn)) &
430	& (Lambdaexo(nn)+1d0)
431	endif
432	enddo
433
434	! print*,'wi',wi(i_h),wi(i_h2),wi(i_d),Uthermal,Lambdaexo,mmol(gcmind(:))
435	! print*,'wi',wi
436
437	! Compute time step for diffusion
438
439	! Loop on species
440
441	T0=Traf(nlraf)
442	rho0=1d0
443
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(:)
452	! Hspecie(nn)=kboltT0/masseinvsgmu
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
461	Tdiff=minval(Tdiffrafmol)
462	Tdiff=minval(Tdiffrafmol(nlraf,:))*Mraf(nlraf)
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
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)
471
472	! Number of time step
473	ntime=anint(dble(ptimestep)/tdiff)
474	! print*,'ptime',ig,ptimestep,tdiff,ntime,tdiffmin,Mraf(nlraf)
475	! Adimensionned temperature
476
477	do l=1,nlraf
478	Tad(l)=Traf(l)/T0
479	enddo
480
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
488	H0=kboltT0/dble(mmol(gcmind(nn)))invsgmu
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)
501	FacEsc=exp(-wad(nn)*DZ/Dad(nlraf-1))
502
503	do l=1,nlraf
504	RhoAd(l)=Rrafk(l,nn)/rho0
505	enddo
506
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
527	! if (ig .eq. ij0 .and. (nn .eq. 1 .or. nn .eq. 5 .or. nn .eq. 6 .or. nn .eq. 16)) then
528	! do l=1,nlraf
529	! if (Xtri(l) .lt. 0.) then
530	! print,'l',l,rhoAd(l)rho0,Xtri(l)*rho0,nn,Tad(l),Zad(l),Dad(l)
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
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
566	! pdqdiff(1:ngrid,1:nlayer,1:nq)=0.
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,&
611	& pp,mmol,gcmind,nlraf,ncompdiff,nlayer,FacMass,ig)
612
613	CALL RHOTOT(pp,tt,massemoy,qnew,RHOT,RHOK,nlayer,ncompdiff)
614
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
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)
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
625	! stop
626
627
628	if (ig .eq. ij0) then
629	do l=il0,nlayer
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,:)),&
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,:)),&
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
637
638	! Compute total mass of each specie on the GCM vertical grid
639
640	Mtot2(1:ncompdiff)=0d0
641
642	do l=il0,nlayer
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
652	! CALL CheckMass2(qq,qnew,pplev(ig,:),il0,nlayer,nn,ncompdiff)
653	! enddo
654
655	! Compute the diffusion trends du to diffusion
656
657	do l=1,nlayer
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
679	! print*,'Escape flux H, H2,D (s-1)',PhiEscH,PhiEscH2,PhiEscD
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)
683
684
685	end subroutine moldiff_red
686
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
790	end subroutine tridagbloc
791
792	subroutine tridag(a,b,c,r,u,n)
793	! parameter (nmax=4000)
794	! dimension gam(nmax),a(n),b(n),c(n),r(n),u(n)
795	real*8 gam(n),a(n),b(n),c(n),r(n),u(n)
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
813	end subroutine tridag
814
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
852	END SUBROUTINE LUBKSB
853
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
938	END SUBROUTINE LUDCMP
939
940	SUBROUTINE TMNEW(T1,DT1,DT2,DT3,T2,dtime,nl,ig)
941	IMPLICIT NONE
942
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
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
958	END SUBROUTINE TMNEW
959
960	SUBROUTINE QMNEW(Q1,DQ,Q2,dtime,nl,nq,gc,ig)
961	use tracer_mod, only: nqmx
962	IMPLICIT NONE
963
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
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
977	END SUBROUTINE QMNEW
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
993	END SUBROUTINE HSCALE
994
995	SUBROUTINE MMOY(massemoy,mmol,qq,gc,nl,nq)
996	use tracer_mod, only: nqmx
997	IMPLICIT NONE
998
999	INTEGER :: nl,nq,l
1000	INTEGER,dimension(nq) :: gc
1001	REAL*8,DIMENSION(nl,nq) :: qq
1002	REAL*8,DIMENSION(nl) :: massemoy
1003	REAL,DIMENSION(nqmx) :: MMOL
1004
1005
1006	do l=1,nl
1007	massemoy(l)=1D0/sum(qq(l,:)/dble(mmol(gc(:))))
1008	enddo
1009
1010	END SUBROUTINE MMOY
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
1024	END SUBROUTINE DMMOY
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
1055	END SUBROUTINE ZVERT
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
1075	END SUBROUTINE RHOTOT
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)
1080	use tracer_mod, only: nqmx
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
1092	REAL,DIMENSION(nqmx) :: mmol
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(:))))
1166	END SUBROUTINE UPPER_RESOL
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
1179	END SUBROUTINE CORRMASS
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
1204	!Temporary: eliminate modification to include Wilke's formulation
1205	!back to the old scheme to check effect
1206	!D(l)=1d0/interm
1207	D(l)=(1D0-Nk(l,nn)/N(l))/interm
1208	enddo
1209	END SUBROUTINE DCOEFF
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
1232	! do l=1,nl
1233	! if (abs(H(l)) .lt. 100.) then
1234	! print*,'H',l,H(l),Nk(l,nn),nn
1235	! endif
1236	! enddo
1237
1238	END SUBROUTINE HSCALEREAL
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
1269	END SUBROUTINE VELVERT
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
1283	END SUBROUTINE TIMEDIFF
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
1290	REAL*8 :: DZ,DZinv
1291	REAL*8,DIMENSION(nl) :: alpha,beta,gama,delta,eps
1292
1293	! Compute alpha,beta and delta
1294	! lower altitude values
1295	DZinv=1d0/(2D0*DZ)
1296
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
1300
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
1312	beta(l)=1D0/T(l)
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
1318
1319	! Upper altitude values
1320
1321	beta(nl)=1D0/T(nl)
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
1324
1325	! Compute the gama and eps coefficients
1326	! Lower altitude values
1327
1328	gama(1)=(-3D0beta(1)+4D0beta(2)-beta(3))*Dzinv
1329	eps(1)=(-3D0alpha(1)+4D0alpha(2)-alpha(3))*Dzinv
1330
1331	gama(nl-1)=(beta(nl)-beta(nl-2))*Dzinv
1332	eps(nl-1)=(alpha(nl)-alpha(nl-2))*Dzinv
1333
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
1338
1339	gama(nl)=(3D0beta(nl)-4D0beta(nl-1)+beta(nl-2))*Dzinv
1340	eps(nl)=(3D0alpha(nl)-4D0alpha(nl-1)+alpha(nl-2))*Dzinv
1341
1342	! do l=1,nl
1343	! print*,'test diffparam',alpha(l),beta(l),delta(l),gama(l),eps(l)
1344	! enddo
1345	! stop
1346
1347	END SUBROUTINE DIFFPARAM
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
1383	END SUBROUTINE MATCOEFF
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
1395	if (abs((Xtot-Ytot)/Xtot) .gt. 1d-3) then
1396	print*,'no conservation for mass',Xtot,Ytot,nn
1397	endif
1398	END SUBROUTINE Checkmass
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
1416	IF (abs(DM/Mold) .gt. 1d-2) THEN
1417	Print*,'We dont conserve mas',nn,DM,Mold,Mnew,DM/Mold
1418	ENDIF
1419
1420	END SUBROUTINE Checkmass2
1421
1422	SUBROUTINE GCMGRID_P(Z,P,Q,T,Nk,Rk,qq,qnew,tt,tnew, &
1423	& pp,M,gc,nl,nq,nlx,ig)
1424	use tracer_mod, only: nqmx
1425	IMPLICIT NONE
1426	INTEGER :: nl,nq,nlx,il,nn,iP,ig,compteur
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
1431	REAL,DIMENSION(nqmx) :: M
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
1483	compteur=0
1484	do while (pnew2 .ge. pp(il))
1485	compteur=compteur+1
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
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
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
1521	END SUBROUTINE GCMGRID_P
1522
1523	SUBROUTINE GCMGRID_P2(Z,P,Q,T,Nk,Rk,qq,qnew,tt,tnew &
1524	& ,pp,M,gc,nl,nq,nlx,facM,ig)
1525	use tracer_mod, only: nqmx
1526	IMPLICIT NONE
1527	INTEGER :: nl,nq,nlx,il,nn,iP,ig,compteur
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
1532	REAL,DIMENSION(nqmx) :: M
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
1583	compteur=0
1584	do while (pnew2 .ge. pp(il))
1585	compteur=compteur+1
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
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
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
1623	END SUBROUTINE GCMGRID_P2
1624
1625
1626	END MODULE moldiff_red_mod

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

Download in other formats:

Original Format