Context Navigation

moldiff_red.F90 @ 2325

Last change on this file since 2325 was 2157, checked in by emillour, 5 years ago

Mars GCM:

Use Wilke's formulation for molecular diffusion

JYC

File size: 39.3 KB

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

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

Download in other formats:

Original Format