Context Navigation

moldiff_MPF.F90

Last change on this file was 3443, checked in by emillour, 3 months ago
Mars PCM: Add new molecular diffusion routine moldiff_MPF (modified pass flow scheme) But for intermetiate testing the legacy moldiff_red routine remains the default. The moldiff_scheme flag (1: legacy, 2: MPF) can be used to swich from one scheme to the other. JYC+EM
File size: 42.4 KB

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

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

Download in other formats:

Original Format