Context Navigation

moldiff_red.F90 @ 2203

Last change on this file since 2203 was 1675, checked in by slebonnois, 8 years ago
SL: update of xml file for xios output + a few modifications for run with uppr atmosphere
File size: 39.1 KB

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

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

Download in other formats:

Original Format