Context Navigation

moldiff_red.F90 @ 3461

Last change on this file since 3461 was 3221, checked in by flefevre, 9 months ago
Venus : a bit more cosmetics in molecular diffusion
File size: 39.8 KB

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

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

Download in other formats:

Original Format