Context Navigation

moldiff_red.F90 @ 1660

Last change on this file since 1660 was 1660, checked in by emillour, 8 years ago

Mars GCM:

Added possibility to run with an Helium "he" tracer (to be initialized at constant value of 3.6e-7 kg/kg_air, i.e. the 4ppm of Krasnopolsky 1996 EUVE satellite, using newstart).
corrected case for CH4 in aeronomars/photochemistry.F (was using undefined indexes when there is no CH4).
updated aki/cpi coefficients for Argon used to compute mean atmospheric Cp and thermal conductivity in aeronomars/concentrations.F

JYC+EM

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	enddo
1177	END
1178
1179	SUBROUTINE HSCALEREAL(nn,Nk,Dz,H,nl,nq)
1180	IMPLICIT NONE
1181	INTEGER :: nn,nl,nq,l
1182	REAL*8,DIMENSION(nl) :: H
1183	REAL*8,DIMENSION(nl,nq) :: Nk
1184	REAL*8 :: Dz
1185
1186	H(1)=(-3D0Nk(1,nn)+4d0NK(2,nn)-Nk(3,nn))/(2D0*DZ)/ &
1187	& NK(1,nn)
1188
1189	H(1)=-1D0/H(1)
1190
1191	DO l=2,nl-1
1192	H(l)=(Nk(l+1,nn)-NK(l-1,nn))/(2D0*DZ)/NK(l,nn)
1193	H(l)=-1D0/H(l)
1194	ENDDO
1195
1196	H(nl)=(3D0Nk(nl,nn)-4D0Nk(nl-1,nn)+Nk(nl-2,nn))/(2D0*DZ)/ &
1197	& Nk(nl,nn)
1198	H(nl)=-1D0/H(nl)
1199
1200	! do l=1,nl
1201	! if (abs(H(l)) .lt. 100.) then
1202	! print*,'H',l,H(l),Nk(l,nn),nn
1203	! endif
1204	! enddo
1205
1206	END
1207
1208	SUBROUTINE VELVERT(nn,T,H,D,Dz,masse,W,nl)
1209	IMPLICIT NONE
1210	INTEGER :: l,nl,nn
1211	REAL*8,DIMENSION(nl) :: T,H,D,W,DT
1212	REAL*8 :: Dz,Hmol,masse
1213	REAL*8 :: kbolt,masseU,Konst,g
1214	masseU=1.660538782d-27
1215	kbolt=1.3806504d-23
1216	Konst=Kbolt/masseU
1217	g=3.72d0
1218
1219	DT(1)=1D0/T(1)(-3D0T(1)+4D0T(2)-T(3))/(2D0DZ)
1220	Hmol=Konst*T(1)/masse/g
1221	W(1)=-D(1)*(1D0/H(1)-1D0/Hmol-DT(1))
1222
1223	DO l=2,nl-1
1224	DT(l)=1D0/T(l)(T(l+1)-T(l-1))/(2D0DZ)
1225	Hmol=Konst*T(l)/masse/g
1226	W(l)=-D(l)*(1D0/H(l)-1D0/Hmol-DT(l))
1227	ENDDO
1228
1229	DT(nl)=1D0/T(nl)(3d0T(nl)-4D0T(nl-1)+T(nl-2))/(2D0DZ)
1230	Hmol=Konst*T(nl)/masse/g
1231	W(nl)=-D(nl)*(1D0/H(nl)-1D0/Hmol-DT(nl))
1232
1233	! do l=1,nl
1234	! print*,'W',W(l),D(l),H(l),DT(l)
1235	! enddo
1236
1237	END
1238
1239	SUBROUTINE TIMEDIFF(nn,H,W,TIME,nl)
1240	IMPLICIT NONE
1241	INTEGER :: nn,nl,l
1242	REAL*8,DIMENSION(nl) :: W,H,TIME
1243
1244	DO l=1,nl
1245	TIME(l)=abs(H(l)/W(l))
1246	if (TIME(l) .lt. 1.D-4) then
1247	! print*,'low tdiff',H(l),W(l),nn,l
1248	endif
1249	ENDDO
1250
1251	END
1252
1253
1254	SUBROUTINE DIFFPARAM(T,D,dz,alpha,beta,gama,delta,eps,nl)
1255	IMPLICIT NONE
1256	INTEGER :: nl,l
1257	REAL*8,DIMENSION(nl) :: T,D
1258	REAL*8 :: DZ,DZinv
1259	REAL*8,DIMENSION(nl) :: alpha,beta,gama,delta,eps
1260
1261	! Compute alpha,beta and delta
1262	! lower altitude values
1263	DZinv=1d0/(2D0*DZ)
1264
1265	beta(1)=1d0/T(1)
1266	alpha(1)=beta(1)(-3D0T(1)+4D0T(2)-T(3))Dzinv
1267	delta(1)=(-3D0D(1)+4D0D(2)-D(3))*Dzinv
1268
1269	beta(2)=1d0/T(2)
1270	alpha(2)=beta(2)(T(3)-T(1))Dzinv
1271	delta(2)=(D(3)-D(1))*Dzinv
1272
1273	! do l=2,nl-1
1274	! beta(l)=1D0/T(l)
1275	! alpha(l)=beta(l)(T(l+1)-T(l-1))Dzinv
1276	! delta(l)=(D(l+1)-D(l-1))*Dzinv
1277	! end do
1278
1279	do l=3,nl-1
1280	beta(l)=1D0/T(l)
1281	alpha(l)=beta(l)(T(l+1)-T(l-1))Dzinv
1282	delta(l)=(D(l+1)-D(l-1))*Dzinv
1283	gama(l-1)=(beta(l)-beta(l-2))*Dzinv
1284	eps(l-1)=(alpha(l)-alpha(l-2))*Dzinv
1285	enddo
1286
1287	! Upper altitude values
1288
1289	beta(nl)=1D0/T(nl)
1290	alpha(nl)=beta(nl)(3D0T(nl)-4D0T(nl-1)+T(nl-2))Dzinv
1291	delta(nl)=(3D0D(nl)-4D0D(nl-1)+D(nl-2))*Dzinv
1292
1293	! Compute the gama and eps coefficients
1294	! Lower altitude values
1295
1296	gama(1)=(-3D0beta(1)+4D0beta(2)-beta(3))*Dzinv
1297	eps(1)=(-3D0alpha(1)+4D0alpha(2)-alpha(3))*Dzinv
1298
1299	gama(nl-1)=(beta(nl)-beta(nl-2))*Dzinv
1300	eps(nl-1)=(alpha(nl)-alpha(nl-2))*Dzinv
1301
1302	! do l=2,nl-1
1303	! gama(l)=(beta(l+1)-beta(l-1))*Dzinv
1304	! eps(l)=(alpha(l+1)-alpha(l-1))*Dzinv
1305	! end do
1306
1307	gama(nl)=(3D0beta(nl)-4D0beta(nl-1)+beta(nl-2))*Dzinv
1308	eps(nl)=(3D0alpha(nl)-4D0alpha(nl-1)+alpha(nl-2))*Dzinv
1309
1310	! do l=1,nl
1311	! print*,'test diffparam',alpha(l),beta(l),delta(l),gama(l),eps(l)
1312	! enddo
1313	! stop
1314
1315	END
1316
1317
1318	SUBROUTINE MATCOEFF(alpha,beta,gama,delta,eps,Dad,rhoad, &
1319	& FacEsc,dz,dt,A,B,C,D,nl)
1320	IMPLICIT NONE
1321	INTEGER :: nl,l
1322	REAL*8, DIMENSION(nl) :: alpha,beta,gama,delta,eps,Dad,RHoad
1323	REAL*8 :: dz,dt,del1,del2,del3,FacEsc
1324	REAL*8, DIMENSION(nl) :: A,B,C,D
1325	del1=dt/2d0/dz
1326	del2=dt/dz/dz
1327	del3=dt
1328
1329	! lower boundary coefficients no change
1330	A(1)=0d0
1331	B(1)=1d0
1332	C(1)=0d0
1333	D(1)=rhoAd(1)
1334
1335	do l=2,nl-1
1336	A(l)=(delta(l)+(alpha(l)+beta(l))Dad(l))del1-Dad(l)*del2
1337	B(l)=-(delta(l)(alpha(l)+beta(l))+Dad(l)(gama(l)+eps(l)))*del3
1338	B(l)=B(l)+1d0+2d0Dad(l)del2
1339	C(l)=-(delta(l)+(alpha(l)+beta(l))Dad(l))del1-Dad(l)*del2
1340	D(l)=rhoAD(l)
1341	enddo
1342
1343
1344	! Upper boundary coefficients Diffusion profile
1345	C(nl)=0d0
1346	B(nl)=-1d0
1347	A(nl)=exp(-dZ(alpha(nl)+beta(nl)))FacEsc
1348	D(nl)=0D0
1349
1350
1351	END
1352
1353	SUBROUTINE Checkmass(X,Y,nl,nn)
1354	IMPLICIT NONE
1355
1356	INTEGER :: nl,nn
1357	REAL*8,DIMENSION(nl) :: X,Y
1358	REAL*8 Xtot,Ytot
1359
1360	Xtot=sum(X)
1361	Ytot=sum(Y)
1362
1363	if (abs((Xtot-Ytot)/Xtot) .gt. 1d-3) then
1364	print*,'no conservation for mass',Xtot,Ytot,nn
1365	endif
1366	END
1367
1368	SUBROUTINE Checkmass2(qold,qnew,P,il,nl,nn,nq)
1369	IMPLICIT NONE
1370	INTEGER :: nl,nn,l,nq,il
1371	REAL,DIMENSION(nl+1) :: P
1372	REAL*8,DIMENSION(nl,nq) :: qold,qnew
1373	REAL*8 :: DM,Mold,Mnew,g
1374	g=3.72d0
1375	DM=0d0
1376	Mold=0d0
1377	Mnew=0d0
1378	DO l=il,nl
1379	DM=DM+(qnew(l,nn)-qold(l,nn))*(dble(P(l))-dble(P(l+1)))/g
1380	Mold=Mold+qold(l,nn)*(dble(P(l))-dble(P(l+1)))/g
1381	Mnew=Mnew+qnew(l,nn)*(dble(P(l))-dble(P(l+1)))/g
1382	! print*,'l',l,qold(l,nn),qnew(l,nn),Mold,Mnew,DM,P(l),P(l+1)
1383	ENDDO
1384	IF (abs(DM/Mold) .gt. 1d-2) THEN
1385	Print*,'We dont conserve mas',nn,DM,Mold,Mnew,DM/Mold
1386	ENDIF
1387
1388	END
1389
1390	SUBROUTINE GCMGRID_P(Z,P,Q,T,Nk,Rk,qq,qnew,tt,tnew, &
1391	& pp,M,gc,nl,nq,nlx,ig)
1392	use tracer_mod, only: nqmx
1393	IMPLICIT NONE
1394	INTEGER :: nl,nq,nlx,il,nn,iP,ig,compteur
1395	INTEGER,DIMENSION(1) :: indP
1396	INTEGER,DIMENSION(nq) :: gc
1397	REAL*8,DIMENSION(nl) :: Z,P,T
1398	REAL*8,DIMENSION(nl,nq) :: Q,Nk,Rk
1399	REAL,DIMENSION(nqmx) :: M
1400	REAL*8,DIMENSION(nq) :: nNew
1401	REAL*8,DIMENSION(nlx) :: pp,tt,tnew
1402	REAL*8,DIMENSION(nlx) :: rhonew
1403	REAL*8,DIMENSION(nlx,nq) :: qq,qnew,rhoknew
1404	REAL*8 :: kbolt,masseU,Konst,g,Dz,facP,Hi
1405	REAL*8 :: Znew,Znew2,Pnew,Pnew2
1406	masseU=1.660538782d-27
1407	kbolt=1.3806504d-23
1408	Konst=Kbolt/masseU
1409	g=3.72d0
1410	Dz=Z(2)-Z(1)
1411	Znew=Z(nl)
1412	Znew2=Znew+dz
1413	! print*,'dz',Znew,Znew2,dz
1414	nNew(1:nq)=Nk(nl,1:nq)
1415	Pnew=P(nl)
1416
1417	do il=1,nlx
1418	! print*,'il',il,pp(il),P(1),P(nl)
1419	if (pp(il) .ge. P(1)) then
1420	qnew(il,:)=qq(il,:)
1421	tnew(il)=tt(il)
1422	endif
1423	if (pp(il) .lt. P(1)) then
1424	if (pp(il) .gt. P(nl)) then
1425	indP=maxloc(P,mask=P < pp(il))
1426	iP=indP(1)-1
1427	if (iP .lt. 1 .or. iP .gt. nl) then
1428	print*,'danger 2',iP,nl,pp(il)
1429	endif
1430	facP=(pp(il)-P(ip))/(P(ip+1)-P(ip))
1431	! print*,'P',P(ip),P(ip+1),facP,indP,iP
1432
1433	! do nn=1,nq
1434	! qnew(il,nn)=Q(iP,nn)+
1435	! & (Q(ip+1,nn)-Q(ip,nn))*facP
1436	! enddo
1437
1438	do nn=1,nq
1439	rhoknew(il,nn)=Rk(iP,nn)+ &
1440	& (Rk(ip+1,nn)-Rk(ip,nn))*facP
1441	enddo
1442	tnew(il)=T(iP)+(T(iP+1)-T(iP))*facP
1443	rhonew(il)=sum(rhoknew(il,:))
1444	do nn=1,nq
1445	qnew(il,nn)=rhoknew(il,nn)/rhonew(il)
1446	enddo
1447
1448	else ! pp < P(nl) need to extrapolate density of each specie
1449	Pnew2=Pnew
1450
1451	compteur=0
1452	do while (pnew2 .ge. pp(il))
1453	compteur=compteur+1
1454	do nn=1,nq
1455	Hi=Konst*T(nl)/dble(M(gc(nn)))/g
1456	Nnew(nn)=Nnew(nn)*exp(-dZ/Hi)
1457	enddo
1458	Pnew=Pnew2
1459	Pnew2=kboltT(nl)sum(Nnew(:))
1460	Znew=Znew2
1461	Znew2=Znew2+Dz
1462	if (compteur .ge. 100000) then
1463	print*,'error moldiff_red infinite loop'
1464	print*,ig,il,pp(il),tt(nl),pnew2,qnew(il,:),Znew2
1465	stop
1466	endif
1467	! print*,'test',Pnew2,Znew2,Nnew(nq),pp(il)
1468	enddo
1469
1470	facP=(pp(il)-Pnew)/(Pnew2-Pnew)
1471
1472	! do nn=1,nq
1473	! qnew(il,nn)=dble(M(gc(nn)))*Nnew(nn)
1474	! & /sum(dble(M(gc(:)))*Nnew(:))
1475	! enddo
1476
1477	do nn=1,nq
1478	rhoknew(il,nn)=dble(M(gc(nn)))*Nnew(nn)
1479	enddo
1480	rhonew(il)=sum(rhoknew(il,:))
1481	do nn=1,nq
1482	qnew(il,nn)=rhoknew(il,nn)/rhonew(il)
1483	enddo
1484	tnew(il)=T(nl)
1485	endif
1486	endif
1487	enddo
1488
1489	END
1490
1491	SUBROUTINE GCMGRID_P2(Z,P,Q,T,Nk,Rk,qq,qnew,tt,tnew &
1492	& ,pp,M,gc,nl,nq,nlx,facM,ig)
1493	use tracer_mod, only: nqmx
1494	IMPLICIT NONE
1495	INTEGER :: nl,nq,nlx,il,nn,iP,ig,compteur
1496	INTEGER,DIMENSION(1) :: indP
1497	INTEGER,DIMENSION(nq) :: gc
1498	REAL*8,DIMENSION(nl) :: Z,P,T
1499	REAL*8,DIMENSION(nl,nq) :: Q,Nk,Rk
1500	REAL,DIMENSION(nqmx) :: M
1501	REAL*8,DIMENSION(nq) :: nNew
1502	REAL*8,DIMENSION(nlx) :: pp,rhonew,tt,tnew
1503	REAL*8,DIMENSION(nlx,nq) :: qq,qnew,facM,rhoknew
1504	REAL*8 :: kbolt,masseU,Konst,g,Dz,facP,Hi
1505	REAL*8 :: Znew,Znew2,Pnew,Pnew2
1506	masseU=1.660538782d-27
1507	kbolt=1.3806504d-23
1508	Konst=Kbolt/masseU
1509	g=3.72d0
1510	Dz=Z(2)-Z(1)
1511	Znew=Z(nl)
1512	Znew2=Znew+dz
1513	! print*,'dz',Znew,Znew2,dz
1514	nNew(1:nq)=Nk(nl,1:nq)
1515	Pnew=P(nl)
1516
1517	do il=1,nlx
1518	! print*,'il',il,pp(il),P(1),P(nl)
1519	if (pp(il) .ge. P(1)) then
1520	qnew(il,:)=qq(il,:)
1521	tnew(il)=tt(il)
1522	endif
1523	if (pp(il) .lt. P(1)) then
1524	if (pp(il) .gt. P(nl)) then
1525	indP=maxloc(P,mask=P < pp(il))
1526	iP=indP(1)-1
1527	if (iP .lt. 1 .or. iP .gt. nl) then
1528	print*,'danger 3',iP,nl,pp(il)
1529	endif
1530	facP=(pp(il)-P(ip))/(P(ip+1)-P(ip))
1531	! print*,'P',P(ip),P(ip+1),facP,indP,iP
1532
1533	! do nn=1,nq
1534	! qnew(il,nn)=Q(iP,nn)+
1535	! & (Q(ip+1,nn)-Q(ip,nn))*facP
1536	! enddo
1537
1538	do nn=1,nq
1539	rhoknew(il,nn)=(RK(iP,nn)+ &
1540	& (RK(iP+1,nn)-Rk(iP,nn))facP)facM(il,nn)
1541	enddo
1542	tnew(il)=T(iP)+(T(ip+1)-T(iP))*facP
1543	rhonew(il)=sum(rhoknew(il,:))
1544	do nn=1,nq
1545	qnew(il,nn)=rhoknew(il,nn)/rhonew(il)
1546	enddo
1547
1548	else ! pp < P(nl) need to extrapolate density of each specie
1549	Pnew2=Pnew
1550
1551	compteur=0
1552	do while (pnew2 .ge. pp(il))
1553	compteur=compteur+1
1554	do nn=1,nq
1555	Hi=Konst*T(nl)/dble(M(gc(nn)))/g
1556	Nnew(nn)=Nnew(nn)*exp(-dZ/Hi)
1557	enddo
1558	Pnew=Pnew2
1559	Pnew2=kboltT(nl)sum(Nnew(:))
1560	Znew=Znew2
1561	Znew2=Znew2+Dz
1562	if (compteur .ge. 100000) then
1563	print*,'pb moldiff_red infinite loop'
1564	print*,ig,nl,T(nl),pnew2,qnew(il,:),Znew2
1565	stop
1566	endif
1567
1568	! print*,'test',Pnew2,Znew2,Nnew(nq),pp(il)
1569	enddo
1570
1571	facP=(pp(il)-Pnew)/(Pnew2-Pnew)
1572
1573	! do nn=1,nq
1574	! qnew(il,nn)=dble(M(gc(nn)))*Nnew(nn)
1575	! & /sum(dble(M(gc(:)))*Nnew(:))
1576	! enddo
1577
1578	do nn=1,nq
1579	rhoknew(il,nn)=dble(M(gc(nn)))Nnew(nn)FacM(il,nn)
1580	enddo
1581	rhonew(il)=sum(rhoknew(il,:))
1582	do nn=1,nq
1583	qnew(il,nn)=rhoknew(il,nn)/rhonew(il)
1584	enddo
1585	tnew(il)=T(nl)
1586
1587	endif
1588	endif
1589	enddo
1590
1591	END
1592

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

Download in other formats:

Original Format