Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

moldiff_MPF.F90

Last change on this file was 3956, checked in by emillour, 3 weeks ago
Venus PCM: Follow-up of r3954: comsctfi.h and YOCST.h no longer exist. Use modules comcstfi_mod and YOMCST_mod instead. EM
File size: 41.7 KB

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

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

Download in other formats:

Original Format