Context Navigation

moldiff_red.F90 @ 1543

Last change on this file since 1543 was 1530, checked in by emillour, 9 years ago

Venus and Titan GCMs:
Updates in the physics to keep up with updates in LMDZ5 (up to
LMDZ5 trunk, rev 2350) concerning dynamics/physics separation:

Adapted makelmdz and makelmdz_fcm script to stop if trying to compile 1d model or newstart or start2archive in parallel.
got rid of references to "dimensions.h" in physics. Within physics packages, use nbp_lon (=iim), nbp_lat (=jjmp1) and nbp_lev (=llm) from module mod_grid_phy_lmdz (in phy_common) instead. Only partially done for Titan, because of many hard-coded commons; a necessary first step will be to clean these up (using modules).

File size: 38.8 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
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 (ig .eq. ij0) then
452	! print*,'test',ig,Tdiff,tdiffmin,minloc(Tdiffrafmol),minloc(Tdiffrafmol(nlraf,:))
453	! endif
454	if (tdiff .lt. tdiffminMraf(nlraf)) tdiff=tdiffminMraf(nlraf)
455	ntime=anint(dble(ptimestep)/Tdiff)
456	! print*,'ptime',ig,ptimestep,Tdiff,ntime,tdiffmin,Mraf(nlraf)
457	! Adimensionned temperature
458
459	do l=1,nlraf
460	Tad(l)=Traf(l)/T0
461	enddo
462
463	do istep=1,ntime
464	do nn=1,ncompdiff
465
466	Draf(1:nlraf)=Drafmol(1:nlraf,nn)
467
468	! Parameters to adimension the problem
469
470	H0=kboltT0/dble(mol_tr(nn))invsgmu
471	D0=Draf(nlraf)
472	Time0=H0*H0/D0
473	Time=Tdiff/Time0
474
475	! Adimensions
476
477	do l=1,nlraf
478	Dad(l)=Draf(l)/D0
479	Zad(l)=Zraf(l)/H0
480	enddo
481	Wad(nn)=wi(nn)*Time0/H0
482	DZ=Zad(2)-Zad(1)
483	FacEsc=exp(-wad(nn)*DZ/Dad(nlraf-1))
484
485	do l=1,nlraf
486	RhoAd(l)=Rrafk(l,nn)/rho0
487	enddo
488
489	! Compute intermediary coefficients
490
491	CALL DIFFPARAM(Tad,Dad,DZ,alpha,beta,gama,delta,eps,nlraf)
492
493	! First way to include escape need to be validated
494	! Compute escape factor (exp(-ueff*dz/D(nz)))
495
496	! Compute matrix coefficients
497
498	CALL MATCOEFF(alpha,beta,gama,delta,eps,Dad,rhoAd,FacEsc, &
499	& dz,time,Atri,Btri,Ctri,Dtri,nlraf)
500
501	Xtri(:)=0D0
502
503	! SOLVE SYSTEM
504
505	CALL tridag(Atri,Btri,Ctri,Dtri,Xtri,nlraf)
506
507	! Xtri=rhoAd
508
509	! if (ig .eq. ij0 .and. (nn .eq. 1 .or. nn .eq. 5 .or. nn .eq. 6 .or. nn .eq. 16)) then
510	! do l=1,nlraf
511	! if (Xtri(l) .lt. 0.) then
512	! print,'l',l,rhoAd(l)rho0,Xtri(l)*rho0,nn,Tad(l),Zad(l),Dad(l)
513	! stop
514	! endif
515	! enddo
516	! endif
517
518	! Check mass conservation of speci
519
520	! CALL CheckMass(rhoAd,Xtri,nlraf,nn)
521
522	! Update mass density of the species
523
524	do l=1,nlraf
525	Rrafk(l,nn)=rho0*Xtri(l)
526	if (Rrafk(l,nn) .ne. Rrafk(l,nn) .or. &
527	& Rrafk(l,nn) .lt. 0 .and. nn .eq. 16) then
528
529	! Test if n(CO2) < 0 skip diffusion (should never happen)
530
531	print*,'pb moldiff',istep,ig,l,nn,Rrafk(l,nn),tdiff,&
532	& rho0*Rhoad(l),Zmin,Zmax,ntime
533	print*,'Atri',Atri
534	print*,'Btri',Btri
535	print*,'Ctri',Ctri
536	print*,'Dtri',Dtri
537	print*,'Xtri',Xtri
538	print*,'alpha',alpha
539	print*,'beta',beta
540	print*,'gama',gama
541	print*,'delta',delta
542	print*,'eps',eps
543	print*,'Dad',Dad
544	print*,'Temp',Traf
545	print*,'alt',Zraf
546	print*,'Mraf',Mraf
547	stop
548	! pdqdiff(1:ngrid,1:nlayer,1:nq)=0.
549	! return
550	! Rrafk(l,nn)=1D-30*Rraf(l)
551	Rrafk(l,nn)=rho0*Rhoad(l)
552
553	endif
554
555	enddo
556
557	enddo ! END Species Loop
558
559	! Update total mass density
560
561	do l=1,nlraf
562	Rraf(l)=sum(Rrafk(l,:))
563	enddo
564
565	! Compute new mass average at each altitude level and new pressure
566
567	do l=1,nlraf
568	do nn=1,ncompdiff
569	Qraf(l,nn)=Rrafk(l,nn)/Rraf(l)
570	Nrafk(l,nn)=Rrafk(l,nn)/dble(mol_tr(nn))/masseU
571	enddo
572	Mraf(l)=1d0/sum(Qraf(l,:)/dble(mol_tr(:)))
573	Nraf(l)=Rraf(l)/Mraf(l)/masseU
574	Praf(l)=Nraf(l)kboltTraf(l)
575	enddo
576
577	enddo ! END time Loop
578
579	! Compute the total mass of each species to check mass conservation
580
581	Mraf2(1:ncompdiff)=0d0
582	do nn=1,ncompdiff
583	do l=1,nlraf
584	Mraf2(nn)=Mraf2(nn)+Rrafk(l,nn)*Dzraf
585	enddo
586	enddo
587
588	! print*,'Mraf',Mraf2
589
590	! Reinterpolate values on the GCM pressure levels
591
592	CALL GCMGRID_P2(Zraf,Praf,Qraf,Traf,Nrafk,Rrafk,qq,qnew,tt,tnew,&
593	& pp,mol_tr,gcmind,nlraf,ncompdiff,nlayer,FacMass,ig)
594
595	CALL RHOTOT(pp,tt,massemoy,qnew,RHOT,RHOK,nlayer,ncompdiff)
596
597	if (ig .eq. ij0) then
598	do l=il0,nlayer
599	write(*,'(i2,1x,19(e12.4,1x))') l,zz(l),tt(l),RHOK(l,1)/sum(RHOK(l,:)),RHOKINIT(l,1)/sum(RHOKINIT(l,:)),&
600	& RHOK(l,2)/sum(RHOK(l,:)),RHOKINIT(l,2)/sum(RHOKINIT(l,:)),&
601	& RHOK(l,6)/sum(RHOK(l,:)),RHOKINIT(l,6)/sum(RHOKINIT(l,:)),&
602	& RHOK(l,5)/sum(RHOK(l,:)),RHOKINIT(l,5)/sum(RHOKINIT(l,:)),&
603	& RHOK(l,7)/sum(RHOK(l,:)),RHOKINIT(l,7)/sum(RHOKINIT(l,:))
604	enddo
605	endif
606
607	! Compute total mass of each specie on the GCM vertical grid
608
609	Mtot2(1:ncompdiff)=0d0
610
611	do l=il0,nlayer
612	do nn=1,ncompdiff
613	Mtot2(nn)=Mtot2(nn)+1d0/gqnew(l,nn) &
614	& (dble(pplev(ig,l))-dble(pplev(ig,l+1)))
615	enddo
616	enddo
617
618	! Check mass conservation of each specie on column
619
620	! do nn=1,ncompdiff
621	! CALL CheckMass2(qq,qnew,pplev(ig,:),il0,nlayer,nn,ncompdiff)
622	! enddo
623
624	! Compute the diffusion trends du to diffusion
625
626	do l=1,nlayer
627	do nn=1,ncompdiff
628	pdqdiff(ig,l,gcmind(nn))=(qnew(l,nn)-qq(l,nn))/ptimestep
629	enddo
630	enddo
631
632	! deallocation des tableaux
633
634	deallocate(Praf,Traf,Rraf,Mraf)
635	deallocate(Nraf,Draf,Hraf,Wraf)
636	deallocate(Zraf,Tdiffraf)
637	deallocate(Prafold,Mrafold)
638	deallocate(Qraf,Rrafk,Nrafk)
639	deallocate(Rrafkold)
640	deallocate(Drafmol,Hrafmol)
641	deallocate(Wrafmol,Tdiffrafmol)
642	deallocate(Atri,Btri,Ctri,Dtri,Xtri)
643	deallocate(Tad,Dad,Zad,rhoad)
644	deallocate(alpha,beta,gama,delta,eps)
645
646
647	enddo ! ig loop
648
649
650	return
651	end
652
653	! ********************************************************************
654	! ********************************************************************
655	! ********************************************************************
656
657	! JYC subtroutine solving MX = Y where M is defined as a block tridiagonal
658	! matrix (Thomas algorithm), tested on a example
659
660	subroutine tridagbloc(M,F,X,n1,n2)
661	parameter (nmax=100)
662	real8 M(n1n2,n1n2),F(n1n2),X(n1*n2)
663	real*8 A(n1,n1,n2),B(n1,n1,n2),C(n1,n1,n2),D(n1,n2)
664	real*8 at(n1,n1),bt(n1,n1),ct(n1,n1),dt(n1),gamt(n1,n1),y(n1,n1)
665	real*8 alf(n1,n1),gam(n1,n1,n2),alfinv(n1,n1)
666	real*8 uvec(n1,n2),uvect(n1),vvect(n1),xt(n1)
667	real*8 indx(n1)
668	real*8 rhu
669	integer n1,n2
670	integer i,p,q
671
672	X(:)=0.
673	! Define the bloc matrix A,B,C and the vector D
674	A(1:n1,1:n1,1)=M(1:n1,1:n1)
675	C(1:n1,1:n1,1)=M(1:n1,n1+1:2*n1)
676	D(1:n1,1)=F(1:n1)
677
678	do i=2,n2-1
679	A(1:n1,1:n1,i)=M((i-1)n1+1:in1,(i-1)n1+1:in1)
680	B(1:n1,1:n1,i)=M((i-1)n1+1:in1,(i-2)n1+1:(i-1)n1)
681	C(1:n1,1:n1,i)=M((i-1)n1+1:in1,in1+1:(i+1)n1)
682	D(1:n1,i)=F((i-1)n1+1:in1)
683	enddo
684	A(1:n1,1:n1,n2)=M((n2-1)n1+1:n2n1,(n2-1)n1+1:n2n1)
685	B(1:n1,1:n1,n2)=M((n2-1)n1+1:n2n1,(n2-2)n1+1:(n2-1)n1)
686	D(1:n1,n2)=F((n2-1)n1+1:n2n1)
687
688	! Initialization
689	y(:,:)=0.
690	do i=1,n1
691	y(i,i)=1.
692	enddo
693
694	at(:,:)=A(:,:,1)
695	ct(:,:)=C(:,:,1)
696	dt(:)=D(:,1)
697	call ludcmp(at,n1,n1,indx,rhu,ierr)
698	do p=1,n1
699	call lubksb(at,n1,n1,indx,y(1,p))
700	do q=1,n1
701	alfinv(q,p)=y(q,p)
702	enddo
703	enddo
704	gamt=matmul(alfinv,ct)
705	gam(:,:,1)=gamt(:,:)
706	uvect=matmul(alfinv,dt)
707	uvec(:,1)=uvect
708
709	do i=2,n2-1
710	y(:,:)=0.
711	do j=1,n1
712	y(j,j)=1.
713	enddo
714	bt(:,:)=B(:,:,i)
715	at(:,:)=A(:,:,i)-matmul(bt,gamt)
716	ct(:,:)=C(:,:,i)
717	dt(:)=D(:,i)
718	call ludcmp(at,n1,n1,indx,rhu,ierr)
719	do p=1,n1
720	call lubksb(at,n1,n1,indx,y(1,p))
721	do q=1,n1
722	alfinv(q,p)=y(q,p)
723	enddo
724	enddo
725	gamt=matmul(alfinv,ct)
726	gam(:,:,i)=gamt
727	vvect=dt-matmul(bt,uvect)
728	uvect=matmul(alfinv,vvect)
729	uvec(:,i)=uvect
730	enddo
731	bt=B(:,:,n2)
732	dt=D(:,n2)
733	at=A(:,:,n2)-matmul(bt,gamt)
734	vvect=dt-matmul(bt,uvect)
735	y(:,:)=0.
736	do j=1,n1
737	y(j,j)=1.
738	enddo
739	call ludcmp(at,n1,n1,indx,rhu,ierr)
740	do p=1,n1
741	call lubksb(at,n1,n1,indx,y(1,p))
742	do q=1,n1
743	alfinv(q,p)=y(q,p)
744	enddo
745	enddo
746	xt=matmul(alfinv,vvect)
747	X((n2-1)n1+1 :n1n2)=xt
748	do i=n2-1,1,-1
749	gamt=gam(:,:,i)
750	xt=X(in1+1:n1n2)
751	uvect=uvec(:,i)
752	vvect=matmul(gamt,xt)
753	X((i-1)n1+1:in1)=uvect-vvect
754	enddo
755
756	end
757
758	subroutine tridag(a,b,c,r,u,n)
759	! parameter (nmax=4000)
760	! dimension gam(nmax),a(n),b(n),c(n),r(n),u(n)
761	real*8 gam(n),a(n),b(n),c(n),r(n),u(n)
762	if(b(1).eq.0.)then
763	stop 'tridag: error: b(1)=0 !!! '
764	endif
765	bet=b(1)
766	u(1)=r(1)/bet
767	do 11 j=2,n
768	gam(j)=c(j-1)/bet
769	bet=b(j)-a(j)*gam(j)
770	if(bet.eq.0.) then
771	stop 'tridag: error: bet=0 !!! '
772	endif
773	u(j)=(r(j)-a(j)*u(j-1))/bet
774	11 continue
775	do 12 j=n-1,1,-1
776	u(j)=u(j)-gam(j+1)*u(j+1)
777	12 continue
778	return
779	end
780
781	! ********************************************************************
782	! ********************************************************************
783	! ********************************************************************
784
785	SUBROUTINE LUBKSB(A,N,NP,INDX,B)
786
787	implicit none
788
789	integer i,j,n,np,ii,ll
790	real*8 sum
791	real*8 a(np,np),indx(np),b(np)
792
793	! DIMENSION A(NP,NP),INDX(N),B(N)
794	II=0
795	DO 12 I=1,N
796	LL=INDX(I)
797	SUM=B(LL)
798	B(LL)=B(I)
799	IF (II.NE.0)THEN
800	DO 11 J=II,I-1
801	SUM=SUM-A(I,J)*B(J)
802	11 CONTINUE
803	ELSE IF (SUM.NE.0.) THEN
804	II=I
805	ENDIF
806	B(I)=SUM
807	12 CONTINUE
808	DO 14 I=N,1,-1
809	SUM=B(I)
810	IF(I.LT.N)THEN
811	DO 13 J=I+1,N
812	SUM=SUM-A(I,J)*B(J)
813	13 CONTINUE
814	ENDIF
815	B(I)=SUM/A(I,I)
816	14 CONTINUE
817	RETURN
818	END
819
820	! ********************************************************************
821	! ********************************************************************
822	! ********************************************************************
823
824	SUBROUTINE LUDCMP(A,N,NP,INDX,D,ierr)
825
826	implicit none
827
828	integer n,np,nmax,i,j,k,imax
829	real*8 d,tiny,aamax
830	real*8 a(np,np),indx(np)
831	integer ierr ! error =0 if OK, =1 if problem
832
833	PARAMETER (NMAX=100,TINY=1.0E-20)
834	! DIMENSION A(NP,NP),INDX(N),VV(NMAX)
835	real*8 sum,vv(nmax),dum
836
837	D=1.
838	DO 12 I=1,N
839	AAMAX=0.
840	DO 11 J=1,N
841	IF (ABS(A(I,J)).GT.AAMAX) AAMAX=ABS(A(I,J))
842	11 CONTINUE
843	IF (AAMAX.EQ.0.) then
844	write(,) 'In moldiff: Problem in LUDCMP with matrix A'
845	write(,) 'Singular matrix ?'
846	write(,) 'Matrix A = ', A
847	! stop
848	! TO DEBUG :
849	ierr =1
850	return
851	! stop
852	END IF
853
854	VV(I)=1./AAMAX
855	12 CONTINUE
856	DO 19 J=1,N
857	IF (J.GT.1) THEN
858	DO 14 I=1,J-1
859	SUM=A(I,J)
860	IF (I.GT.1)THEN
861	DO 13 K=1,I-1
862	SUM=SUM-A(I,K)*A(K,J)
863	13 CONTINUE
864	A(I,J)=SUM
865	ENDIF
866	14 CONTINUE
867	ENDIF
868	AAMAX=0.
869	DO 16 I=J,N
870	SUM=A(I,J)
871	IF (J.GT.1)THEN
872	DO 15 K=1,J-1
873	SUM=SUM-A(I,K)*A(K,J)
874	15 CONTINUE
875	A(I,J)=SUM
876	ENDIF
877	DUM=VV(I)*ABS(SUM)
878	IF (DUM.GE.AAMAX) THEN
879	IMAX=I
880	AAMAX=DUM
881	ENDIF
882	16 CONTINUE
883	IF (J.NE.IMAX)THEN
884	DO 17 K=1,N
885	DUM=A(IMAX,K)
886	A(IMAX,K)=A(J,K)
887	A(J,K)=DUM
888	17 CONTINUE
889	D=-D
890	VV(IMAX)=VV(J)
891	ENDIF
892	INDX(J)=IMAX
893	IF(J.NE.N)THEN
894	IF(A(J,J).EQ.0.)A(J,J)=TINY
895	DUM=1./A(J,J)
896	DO 18 I=J+1,N
897	A(I,J)=A(I,J)*DUM
898	18 CONTINUE
899	ENDIF
900	19 CONTINUE
901	IF(A(N,N).EQ.0.)A(N,N)=TINY
902	ierr =0
903	RETURN
904	END
905
906	SUBROUTINE TMNEW(T1,DT1,DT2,DT3,T2,dtime,nl,ig)
907	IMPLICIT NONE
908
909	INTEGER,INTENT(IN) :: nl,ig
910	REAL,INTENT(IN),DIMENSION(nl) :: T1,DT1,DT2,DT3
911	REAL*8,INTENT(OUT),DIMENSION(nl) :: T2
912	REAL,INTENT(IN) :: dtime
913	INTEGER :: l
914	DO l=1,nl
915	T2(l)=T1(l)1D0+(DT1(l)dble(dtime)+ &
916	& DT2(l)*dble(dtime)+ &
917	& DT3(l)dble(dtime))1D0
918	if (T2(l) .ne. T2(l)) then
919	print*,'Err TMNEW',ig,l,T2(l),T1(l),dT1(l),DT2(l), &
920	& DT3(l),dtime,dble(dtime)
921	endif
922
923	ENDDO
924	END
925
926	SUBROUTINE QMNEW(Q1,DQ,Q2,dtime,nl,nq,gc,ig)
927	use chemparam_mod
928	use infotrac
929	IMPLICIT NONE
930
931	INTEGER,INTENT(IN) :: nl,nq
932	INTEGER,INTENT(IN) :: ig
933	INTEGER,INTENT(IN),dimension(nq) :: gc
934	REAL,INTENT(IN),DIMENSION(nl,nqtot) :: Q1,DQ
935	REAL*8,INTENT(OUT),DIMENSION(nl,nq) :: Q2
936	REAL,INTENT(IN) :: dtime
937	INTEGER :: l,iq
938	DO l=1,nl
939	DO iq=1,nq
940	Q2(l,iq)=Q1(l,gc(iq))1D0+(DQ(l,gc(iq))dble(dtime))*1D0
941	Q2(l,iq)=max(Q2(l,iq),1d-30)
942	ENDDO
943	ENDDO
944	END
945
946	SUBROUTINE HSCALE(p,hp,nl)
947	IMPLICIT NONE
948
949	INTEGER :: nl
950	INTEGER :: l
951	REAL*8,dimension(nl) :: P
952	REAL*8,DIMENSION(nl) :: Hp
953
954	hp(1)=-log(P(2)/P(1))
955	hp(nl)=-log(P(nl)/P(nl-1))
956
957	DO l=2,nl-1
958	hp(l)=-log(P(l+1)/P(l-1))
959	ENDDO
960	END
961
962	SUBROUTINE MMOY(massemoy,mol_tr,qq,gc,nl,nq)
963	use chemparam_mod
964	use infotrac
965	IMPLICIT NONE
966
967	INTEGER :: nl,nq,l, nn
968	INTEGER,dimension(nq) :: gc
969	REAL*8,DIMENSION(nl,nq) :: qq
970	REAL*8,DIMENSION(nl) :: massemoy
971	REAL,DIMENSION(nq) :: mol_tr
972
973	do l=1,nl
974	massemoy(l)=1D0/sum(qq(l,:)/dble(mol_tr(:)))
975	! write(,),'L988 Masse molaire moy: ', massemoy(l)
976	enddo
977
978	END
979
980	SUBROUTINE DMMOY(M,H,DM,nl)
981	IMPLICIT NONE
982	INTEGER :: nl,l
983	REAL*8,DIMENSION(nl) :: M,H,DM
984
985	DM(1)=(-3D0M(1)+4D0M(2)-M(3))/2D0/H(1)
986	DM(nl)=(3D0M(nl)-4D0M(nl-1)+M(nl-2))/2D0/H(nl)
987
988	do l=2,nl-1
989	DM(l)=(M(l+1)-M(l-1))/H(l)
990	enddo
991
992	END
993
994	SUBROUTINE ZVERT(P,T,M,Z,nl,ig)
995	IMPLICIT NONE
996	INTEGER :: nl,l,ig
997	REAL*8,dimension(nl) :: P,T,M,Z,H
998	REAL*8 :: z0
999	REAL*8 :: kbolt,masseU,Konst,g,Hpm
1000	masseU=1.660538782d-27
1001	kbolt=1.3806504d-23
1002	Konst=kbolt/masseU
1003	g=3.72D0
1004
1005	z0=0d0
1006	Z(1)=z0
1007	H(1)=Konst*T(1)/M(1)/g
1008
1009	do l=2,nl
1010	H(l)=Konst*T(l)/M(l)/g
1011	Hpm=H(l-1)
1012	Z(l)=z(l-1)-Hpm*log(P(l)/P(l-1))
1013	if (Z(l) .ne. Z(l)) then
1014	print*,'pb',l,ig
1015	print*,'P',P
1016	print*,'T',T
1017	print*,'M',M
1018	print*,'Z',Z
1019	print*,'Hpm',Hpm
1020	endif
1021	enddo
1022
1023	END
1024
1025	SUBROUTINE RHOTOT(P,T,M,qq,rhoN,rhoK,nl,nq)
1026	IMPLICIT NONE
1027
1028	REAL*8 :: kbolt,masseU,Konst
1029	INTEGER :: nl,nq,l,iq
1030	REAL*8,DIMENSION(nl) :: P,T,M,RHON
1031	REAL*8,DIMENSION(nl,nq) :: RHOK,qq
1032	masseU=1.660538782d-27
1033	kbolt=1.3806504d-23
1034	Konst=Kbolt/masseU
1035
1036	do l=1,nl
1037	RHON(l)=P(l)*M(l)/T(l)/Konst
1038	do iq=1,nq
1039	RHOK(l,iq)=qq(l,iq)*RHON(l)
1040	enddo
1041	enddo
1042
1043	END
1044
1045	SUBROUTINE UPPER_RESOL(P,T,Z,M,R,Rk, &
1046	& qq,mol_tr,gc,Praf,Traf,Qraf,Mraf,Zraf, &
1047	& Nraf,Nrafk,Rraf,Rrafk,il,nl,nq,nlx,ig)
1048	use chemparam_mod
1049	use infotrac
1050	IMPLICIT NONE
1051
1052	INTEGER :: nl,nq,il,l,i,iq,nlx,iz,ig
1053	INTEGER :: gc(nq)
1054	INTEGER,DIMENSION(1) :: indz
1055	REAL*8, DIMENSION(nlx) :: P,T,Z,M,R
1056	REAL*8, DIMENSION(nlx,nq) :: qq,Rk
1057	REAL*8, DIMENSION(nl) :: Praf,Traf,Mraf,Zraf,Nraf,Rraf
1058	REAL*8 :: kbolt,masseU,Konst,g
1059	REAL*8, DIMENSION(nl,nq) :: Qraf,Rrafk,Nrafk
1060	REAL*8 :: facZ,dZ,H
1061	REAL,DIMENSION(nq) :: mol_tr
1062	masseU=1.660538782d-27
1063	kbolt=1.3806504d-23
1064	Konst=Kbolt/masseU
1065	g=3.72d0
1066
1067
1068	Zraf(1)=z(il)
1069	Praf(1)=P(il)
1070	Traf(1)=T(il)
1071	Nraf(1)=Praf(1)/kbolt/Traf(1)
1072	do iq=1,nq
1073	Qraf(1,iq)=qq(il,iq)
1074	enddo
1075	Mraf(1)=1d0/sum(Qraf(1,:)/dble(mol_tr(:)))
1076	Rraf(1)=Mraf(1)masseUNraf(1)
1077	do iq=1,nq
1078	Rrafk(1,iq)=Rraf(1)*Qraf(1,iq)
1079	Nrafk(1,iq)=Rrafk(1,iq)/masseU/dble(mol_tr(iq))
1080	enddo
1081	Zraf(nl)=z(nlx)
1082
1083	do l=2,nl-1
1084	Zraf(l)=Zraf(1)+(Zraf(nl)-Zraf(1))/dble(nl-1)*dble((l-1))
1085	indz=maxloc(z,mask=z <= Zraf(l))
1086	iz=indz(1)
1087	if (iz .lt. 1 .or. iz .gt. nlx) then
1088	print*,'danger',iz,nl,Zraf(l),l,Zraf(1),Zraf(nl),z,P,T,ig
1089	stop
1090	endif
1091	dZ=Zraf(l)-Zraf(l-1)
1092	! dZ=Zraf(l)-z(iz)
1093	facz=(Zraf(l)-z(iz))/(z(iz+1)-z(iz))
1094	Traf(l)=T(iz)+(T(iz+1)-T(iz))*facz
1095	do iq=1,nq
1096	! Qraf(l,iq)=qq(iz,iq)+(qq(iz+1,iq)-qq(iz,iq))*facz
1097	Rrafk(l,iq)=Rk(iz,iq)+(Rk(iz+1,iq)-Rk(iz,iq))*facZ
1098	Rrafk(l,iq)=Rk(iz,iq)(Rk(iz+1,iq)/Rk(iz,iq))*facZ
1099	enddo
1100	! Mraf(l)=1D0/(sum(qraf(l,:)/dble(mol_tr(:))))
1101	Rraf(l)=sum(Rrafk(l,:))
1102	do iq=1,nq
1103	Qraf(l,iq)=Rrafk(l,iq)/Rraf(l)
1104	enddo
1105	Mraf(l)=1D0/(sum(qraf(l,:)/dble(mol_tr(:))))
1106	! H=Konst*Traf(l)/Mraf(l)/g
1107	! H=Konst*T(iz)/M(iz)/g
1108	! Praf(l)=P(iz)*exp(-dZ/H)
1109	! Praf(l)=Praf(l-1)*exp(-dZ/H)
1110	! print,'iz',l,iz,Praf(il-1)exp(-dZ/H),z(iz),z(iz+1),H
1111	Nraf(l)=Rraf(l)/Mraf(l)/masseU
1112	Praf(l)=Nraf(l)kboltTraf(l)
1113	! Rraf(l)=Nraf(l)Mraf(l)masseU
1114	do iq=1,nq
1115	! Rrafk(l,iq)=Rraf(l)*Qraf(l,iq)
1116	Nrafk(l,iq)=Rrafk(l,iq)/dble(mol_tr(iq))/masseU
1117	if (Nrafk(l,iq) .lt. 0. .or. &
1118	& Nrafk(l,iq) .ne. Nrafk(l,iq)) then
1119	print*,'pb interpolation',l,iq,Nrafk(l,iq),Rrafk(l,iq), &
1120	& Qraf(l,iq),Rk(iz,iq),Rk(iz+1,iq),facZ,Zraf(l),z(iz)
1121	stop
1122	endif
1123	enddo
1124	enddo
1125	Zraf(nl)=z(nlx)
1126	Traf(nl)=T(nlx)
1127	do iq=1,nq
1128	Rrafk(nl,iq)=Rk(nlx,iq)
1129	Qraf(nl,iq)=Rk(nlx,iq)/R(nlx)
1130	Nrafk(nl,iq)=Rk(nlx,iq)/dble(mol_tr(iq))/masseU
1131	enddo
1132	Nraf(nl)=sum(Nrafk(nl,:))
1133	Praf(nl)=Nraf(nl)kboltTraf(nl)
1134	Mraf(nl)=1D0/sum(Qraf(nl,:)/dble(mol_tr(:)))
1135	END
1136
1137	SUBROUTINE CORRMASS(qq,qint,FacMass,nl,nq)
1138	IMPLICIT NONE
1139	INTEGER :: nl,nq,l,nn
1140	REAL*8,DIMENSION(nl,nq) :: qq,qint,FacMass
1141
1142	do nn=1,nq
1143	do l=1,nl
1144	FacMass(l,nn)=qq(l,nn)/qint(l,nn)
1145	enddo
1146	enddo
1147
1148	END
1149
1150
1151	SUBROUTINE DCOEFF(nn,Dij,P,T,N,Nk,D,nl,nq)
1152	IMPLICIT NONE
1153	REAL*8,DIMENSION(nl) :: N,T,D,P
1154	REAL*8,DIMENSION(nl,nq) :: Nk
1155	INTEGER :: nn,nl,nq,l,iq
1156	REAL,DIMENSION(nq,nq) :: Dij
1157	REAL*8 :: interm,P0,T0,ptfac,dfac
1158
1159	P0=1D5
1160	T0=273d0
1161
1162
1163	do l=1,nl
1164	ptfac=(P0/P(l))(T(l)/T0)*1.75d0
1165	D(l)=0d0
1166	interm=0d0
1167	do iq=1,nq
1168	if (iq .ne. nn) then
1169	dfac=dble(dij(nn,iq))*ptfac
1170	interm=interm+Nk(l,iq)/N(l)/dfac
1171	endif
1172	enddo
1173	D(l)=1d0/interm
1174	enddo
1175	END
1176
1177	SUBROUTINE HSCALEREAL(nn,Nk,Dz,H,nl,nq)
1178	IMPLICIT NONE
1179	INTEGER :: nn,nl,nq,l
1180	REAL*8,DIMENSION(nl) :: H
1181	REAL*8,DIMENSION(nl,nq) :: Nk
1182	REAL*8 :: Dz
1183
1184	H(1)=(-3D0Nk(1,nn)+4d0NK(2,nn)-Nk(3,nn))/(2D0*DZ)/ &
1185	& NK(1,nn)
1186
1187	H(1)=-1D0/H(1)
1188
1189	DO l=2,nl-1
1190	H(l)=(Nk(l+1,nn)-NK(l-1,nn))/(2D0*DZ)/NK(l,nn)
1191	H(l)=-1D0/H(l)
1192	ENDDO
1193
1194	H(nl)=(3D0Nk(nl,nn)-4D0Nk(nl-1,nn)+Nk(nl-2,nn))/(2D0*DZ)/ &
1195	& Nk(nl,nn)
1196	H(nl)=-1D0/H(nl)
1197
1198	! do l=1,nl
1199	! if (abs(H(l)) .lt. 100.) then
1200	! print*,'H',l,H(l),Nk(l,nn),nn
1201	! endif
1202	! enddo
1203
1204	END
1205
1206	SUBROUTINE VELVERT(nn,T,H,D,Dz,masse,W,nl)
1207	IMPLICIT NONE
1208	INTEGER :: l,nl,nn
1209	REAL*8,DIMENSION(nl) :: T,H,D,W,DT
1210	REAL*8 :: Dz,Hmol,masse
1211	REAL*8 :: kbolt,masseU,Konst,g
1212	masseU=1.660538782d-27
1213	kbolt=1.3806504d-23
1214	Konst=Kbolt/masseU
1215	g=3.72d0
1216
1217	DT(1)=1D0/T(1)(-3D0T(1)+4D0T(2)-T(3))/(2D0DZ)
1218	Hmol=Konst*T(1)/masse/g
1219	W(1)=-D(1)*(1D0/H(1)-1D0/Hmol-DT(1))
1220
1221	DO l=2,nl-1
1222	DT(l)=1D0/T(l)(T(l+1)-T(l-1))/(2D0DZ)
1223	Hmol=Konst*T(l)/masse/g
1224	W(l)=-D(l)*(1D0/H(l)-1D0/Hmol-DT(l))
1225	ENDDO
1226
1227	DT(nl)=1D0/T(nl)(3d0T(nl)-4D0T(nl-1)+T(nl-2))/(2D0DZ)
1228	Hmol=Konst*T(nl)/masse/g
1229	W(nl)=-D(nl)*(1D0/H(nl)-1D0/Hmol-DT(nl))
1230
1231	! do l=1,nl
1232	! print*,'W',W(l),D(l),H(l),DT(l)
1233	! enddo
1234
1235	END
1236
1237	SUBROUTINE TIMEDIFF(nn,H,W,TIME,nl)
1238	IMPLICIT NONE
1239	INTEGER :: nn,nl,l
1240	REAL*8,DIMENSION(nl) :: W,H,TIME
1241
1242	DO l=1,nl
1243	TIME(l)=abs(H(l)/W(l))
1244	if (TIME(l) .lt. 1.D-4) then
1245	! print*,'low tdiff',H(l),W(l),nn,l
1246	endif
1247	ENDDO
1248
1249	END
1250
1251
1252	SUBROUTINE DIFFPARAM(T,D,dz,alpha,beta,gama,delta,eps,nl)
1253	IMPLICIT NONE
1254	INTEGER :: nl,l
1255	REAL*8,DIMENSION(nl) :: T,D
1256	REAL*8 :: DZ,DZinv
1257	REAL*8,DIMENSION(nl) :: alpha,beta,gama,delta,eps
1258
1259	! Compute alpha,beta and delta
1260	! lower altitude values
1261	DZinv=1d0/(2D0*DZ)
1262
1263	beta(1)=1d0/T(1)
1264	alpha(1)=beta(1)(-3D0T(1)+4D0T(2)-T(3))Dzinv
1265	delta(1)=(-3D0D(1)+4D0D(2)-D(3))*Dzinv
1266
1267	beta(2)=1d0/T(2)
1268	alpha(2)=beta(2)(T(3)-T(1))Dzinv
1269	delta(2)=(D(3)-D(1))*Dzinv
1270
1271	! do l=2,nl-1
1272	! beta(l)=1D0/T(l)
1273	! alpha(l)=beta(l)(T(l+1)-T(l-1))Dzinv
1274	! delta(l)=(D(l+1)-D(l-1))*Dzinv
1275	! end do
1276
1277	do l=3,nl-1
1278	beta(l)=1D0/T(l)
1279	alpha(l)=beta(l)(T(l+1)-T(l-1))Dzinv
1280	delta(l)=(D(l+1)-D(l-1))*Dzinv
1281	gama(l-1)=(beta(l)-beta(l-2))*Dzinv
1282	eps(l-1)=(alpha(l)-alpha(l-2))*Dzinv
1283	enddo
1284
1285	! Upper altitude values
1286
1287	beta(nl)=1D0/T(nl)
1288	alpha(nl)=beta(nl)(3D0T(nl)-4D0T(nl-1)+T(nl-2))Dzinv
1289	delta(nl)=(3D0D(nl)-4D0D(nl-1)+D(nl-2))*Dzinv
1290
1291	! Compute the gama and eps coefficients
1292	! Lower altitude values
1293
1294	gama(1)=(-3D0beta(1)+4D0beta(2)-beta(3))*Dzinv
1295	eps(1)=(-3D0alpha(1)+4D0alpha(2)-alpha(3))*Dzinv
1296
1297	gama(nl-1)=(beta(nl)-beta(nl-2))*Dzinv
1298	eps(nl-1)=(alpha(nl)-alpha(nl-2))*Dzinv
1299
1300	! do l=2,nl-1
1301	! gama(l)=(beta(l+1)-beta(l-1))*Dzinv
1302	! eps(l)=(alpha(l+1)-alpha(l-1))*Dzinv
1303	! end do
1304
1305	gama(nl)=(3D0beta(nl)-4D0beta(nl-1)+beta(nl-2))*Dzinv
1306	eps(nl)=(3D0alpha(nl)-4D0alpha(nl-1)+alpha(nl-2))*Dzinv
1307
1308	! do l=1,nl
1309	! print*,'test diffparam',alpha(l),beta(l),delta(l),gama(l),eps(l)
1310	! enddo
1311	! stop
1312
1313	END
1314
1315
1316	SUBROUTINE MATCOEFF(alpha,beta,gama,delta,eps,Dad,rhoad, &
1317	& FacEsc,dz,dt,A,B,C,D,nl)
1318	IMPLICIT NONE
1319	INTEGER :: nl,l
1320	REAL*8, DIMENSION(nl) :: alpha,beta,gama,delta,eps,Dad,RHoad
1321	REAL*8 :: dz,dt,del1,del2,del3,FacEsc
1322	REAL*8, DIMENSION(nl) :: A,B,C,D
1323	del1=dt/2d0/dz
1324	del2=dt/dz/dz
1325	del3=dt
1326
1327	! lower boundary coefficients no change
1328	A(1)=0d0
1329	B(1)=1d0
1330	C(1)=0d0
1331	D(1)=rhoAd(1)
1332
1333	do l=2,nl-1
1334	A(l)=(delta(l)+(alpha(l)+beta(l))Dad(l))del1-Dad(l)*del2
1335	B(l)=-(delta(l)(alpha(l)+beta(l))+Dad(l)(gama(l)+eps(l)))*del3
1336	B(l)=B(l)+1d0+2d0Dad(l)del2
1337	C(l)=-(delta(l)+(alpha(l)+beta(l))Dad(l))del1-Dad(l)*del2
1338	D(l)=rhoAD(l)
1339	enddo
1340
1341
1342	! Upper boundary coefficients Diffusion profile
1343	C(nl)=0d0
1344	B(nl)=-1d0
1345	A(nl)=exp(-dZ(alpha(nl)+beta(nl)))FacEsc
1346	D(nl)=0D0
1347
1348
1349	END
1350
1351	SUBROUTINE Checkmass(X,Y,nl,nn)
1352	IMPLICIT NONE
1353
1354	INTEGER :: nl,nn
1355	REAL*8,DIMENSION(nl) :: X,Y
1356	REAL*8 Xtot,Ytot
1357
1358	Xtot=sum(X)
1359	Ytot=sum(Y)
1360
1361	if (abs((Xtot-Ytot)/Xtot) .gt. 1d-4) then
1362	print*,'no conservation for mass',Xtot,Ytot,nn
1363	endif
1364	END
1365
1366	SUBROUTINE Checkmass2(qold,qnew,P,il,nl,nn,nq)
1367	IMPLICIT NONE
1368	INTEGER :: nl,nn,l,nq,il
1369	REAL,DIMENSION(nl+1) :: P
1370	REAL*8,DIMENSION(nl,nq) :: qold,qnew
1371	REAL*8 :: DM,Mold,Mnew,g
1372	g=3.72d0
1373	DM=0d0
1374	Mold=0d0
1375	Mnew=0d0
1376	DO l=il,nl
1377	DM=DM+(qnew(l,nn)-qold(l,nn))*(dble(P(l))-dble(P(l+1)))/g
1378	Mold=Mold+qold(l,nn)*(dble(P(l))-dble(P(l+1)))/g
1379	Mnew=Mnew+qnew(l,nn)*(dble(P(l))-dble(P(l+1)))/g
1380	! print*,'l',l,qold(l,nn),qnew(l,nn),Mold,Mnew,DM,P(l),P(l+1)
1381	ENDDO
1382	IF (abs(DM/Mold) .gt. 1d-5) THEN
1383	Print*,'We dont conserve mas',nn,DM,Mold,Mnew,DM/Mold
1384	ENDIF
1385
1386	END
1387
1388	SUBROUTINE GCMGRID_P(Z,P,Q,T,Nk,Rk,qq,qnew,tt,tnew, &
1389	& pp,M,gc,nl,nq,nlx,ig)
1390	use chemparam_mod
1391	use infotrac
1392	IMPLICIT NONE
1393	INTEGER :: nl,nq,nlx,il,nn,iP,ig,compteur
1394	INTEGER,DIMENSION(1) :: indP
1395	INTEGER,DIMENSION(nq) :: gc
1396	REAL*8,DIMENSION(nl) :: Z,P,T
1397	REAL*8,DIMENSION(nl,nq) :: Q,Nk,Rk
1398	REAL,DIMENSION(nq) :: M
1399	REAL*8,DIMENSION(nq) :: nNew
1400	REAL*8,DIMENSION(nlx) :: pp,tt,tnew
1401	REAL*8,DIMENSION(nlx) :: rhonew
1402	REAL*8,DIMENSION(nlx,nq) :: qq,qnew,rhoknew
1403	REAL*8 :: kbolt,masseU,Konst,g,Dz,facP,Hi
1404	REAL*8 :: Znew,Znew2,Pnew,Pnew2
1405	masseU=1.660538782d-27
1406	kbolt=1.3806504d-23
1407	Konst=Kbolt/masseU
1408	g=3.72d0
1409	Dz=Z(2)-Z(1)
1410	Znew=Z(nl)
1411	Znew2=Znew+dz
1412	! print*,'dz',Znew,Znew2,dz
1413	nNew(1:nq)=Nk(nl,1:nq)
1414	Pnew=P(nl)
1415
1416	do il=1,nlx
1417	! print*,'il',il,pp(il),P(1),P(nl)
1418	if (pp(il) .ge. P(1)) then
1419	qnew(il,:)=qq(il,:)
1420	tnew(il)=tt(il)
1421	endif
1422	if (pp(il) .lt. P(1)) then
1423	if (pp(il) .gt. P(nl)) then
1424	indP=maxloc(P,mask=P < pp(il))
1425	iP=indP(1)-1
1426	if (iP .lt. 1 .or. iP .gt. nl) then
1427	print*,'danger 2',iP,nl,pp(il)
1428	endif
1429	facP=(pp(il)-P(ip))/(P(ip+1)-P(ip))
1430	! print*,'P',P(ip),P(ip+1),facP,indP,iP
1431
1432	! do nn=1,nq
1433	! qnew(il,nn)=Q(iP,nn)+
1434	! & (Q(ip+1,nn)-Q(ip,nn))*facP
1435	! enddo
1436
1437	do nn=1,nq
1438	rhoknew(il,nn)=Rk(iP,nn)+ &
1439	& (Rk(ip+1,nn)-Rk(ip,nn))*facP
1440	enddo
1441	tnew(il)=T(iP)+(T(iP+1)-T(iP))*facP
1442	rhonew(il)=sum(rhoknew(il,:))
1443	do nn=1,nq
1444	qnew(il,nn)=rhoknew(il,nn)/rhonew(il)
1445	enddo
1446
1447	else ! pp < P(nl) need to extrapolate density of each specie
1448	Pnew2=Pnew
1449
1450	compteur=0
1451	do while (pnew2 .ge. pp(il))
1452	compteur=compteur+1
1453	do nn=1,nq
1454	Hi=Konst*T(nl)/dble(M(nn))/g
1455	Nnew(nn)=Nnew(nn)*exp(-dZ/Hi)
1456	enddo
1457	Pnew=Pnew2
1458	Pnew2=kboltT(nl)sum(Nnew(:))
1459	Znew=Znew2
1460	Znew2=Znew2+Dz
1461	if (compteur .ge. 100000) then
1462	print*,'error moldiff_red infinite loop'
1463	print*,ig,il,pp(il),tt(nl),pnew2,qnew(il,:),Znew2
1464	stop
1465	endif
1466	! print*,'test',Pnew2,Znew2,Nnew(nq),pp(il)
1467	enddo
1468
1469	facP=(pp(il)-Pnew)/(Pnew2-Pnew)
1470
1471	! do nn=1,nq
1472	! qnew(il,nn)=dble(M(nn))*Nnew(nn)
1473	! & /sum(dble(M(:))*Nnew(:))
1474	! enddo
1475
1476	do nn=1,nq
1477	rhoknew(il,nn)=dble(M(nn))*Nnew(nn)
1478	enddo
1479	rhonew(il)=sum(rhoknew(il,:))
1480	do nn=1,nq
1481	qnew(il,nn)=rhoknew(il,nn)/rhonew(il)
1482	enddo
1483	tnew(il)=T(nl)
1484	endif
1485	endif
1486	enddo
1487
1488	END
1489
1490	SUBROUTINE GCMGRID_P2(Z,P,Q,T,Nk,Rk,qq,qnew,tt,tnew &
1491	& ,pp,M,gc,nl,nq,nlx,facM,ig)
1492	! use chemparam_mod
1493	! use infotrac
1494	IMPLICIT NONE
1495	INTEGER :: nl,nq,nlx,il,nn,iP,ig,compteur
1496	INTEGER,DIMENSION(1) :: indP
1497	INTEGER,DIMENSION(nq) :: gc
1498	REAL*8,DIMENSION(nl) :: Z,P,T
1499	REAL*8,DIMENSION(nl,nq) :: Q,Nk,Rk
1500	REAL,DIMENSION(nq) :: M
1501	REAL*8,DIMENSION(nq) :: nNew
1502	REAL*8,DIMENSION(nlx) :: pp,rhonew,tt,tnew
1503	REAL*8,DIMENSION(nlx,nq) :: qq,qnew,facM,rhoknew
1504	REAL*8 :: kbolt,masseU,Konst,g,Dz,facP,Hi
1505	REAL*8 :: Znew,Znew2,Pnew,Pnew2
1506	masseU=1.660538782d-27
1507	kbolt=1.3806504d-23
1508	Konst=Kbolt/masseU
1509	g=3.72d0
1510	Dz=Z(2)-Z(1)
1511	Znew=Z(nl)
1512	Znew2=Znew+dz
1513	! print*,'dz',Znew,Znew2,dz
1514	nNew(1:nq)=Nk(nl,1:nq)
1515	Pnew=P(nl)
1516
1517	do il=1,nlx
1518	! print*,'il',il,pp(il),P(1),P(nl)
1519	if (pp(il) .ge. P(1)) then
1520	qnew(il,:)=qq(il,:)
1521	tnew(il)=tt(il)
1522	endif
1523	if (pp(il) .lt. P(1)) then
1524	if (pp(il) .gt. P(nl)) then
1525	indP=maxloc(P,mask=P < pp(il))
1526	iP=indP(1)-1
1527	if (iP .lt. 1 .or. iP .gt. nl) then
1528	print*,'danger 3',iP,nl,pp(il)
1529	endif
1530	facP=(pp(il)-P(ip))/(P(ip+1)-P(ip))
1531	! print*,'P',P(ip),P(ip+1),facP,indP,iP
1532
1533	! do nn=1,nq
1534	! qnew(il,nn)=Q(iP,nn)+
1535	! & (Q(ip+1,nn)-Q(ip,nn))*facP
1536	! enddo
1537
1538	do nn=1,nq
1539	rhoknew(il,nn)=(RK(iP,nn)+ &
1540	& (RK(iP+1,nn)-Rk(iP,nn))facP)facM(il,nn)
1541	enddo
1542	tnew(il)=T(iP)+(T(ip+1)-T(iP))*facP
1543	rhonew(il)=sum(rhoknew(il,:))
1544	do nn=1,nq
1545	qnew(il,nn)=rhoknew(il,nn)/rhonew(il)
1546	enddo
1547
1548	else ! pp < P(nl) need to extrapolate density of each specie
1549	Pnew2=Pnew
1550
1551	compteur=0
1552	do while (pnew2 .ge. pp(il))
1553	compteur=compteur+1
1554	do nn=1,nq
1555	Hi=Konst*T(nl)/dble(M(nn))/g
1556	Nnew(nn)=Nnew(nn)*exp(-dZ/Hi)
1557	enddo
1558	Pnew=Pnew2
1559	Pnew2=kboltT(nl)sum(Nnew(:))
1560	Znew=Znew2
1561	Znew2=Znew2+Dz
1562	if (compteur .ge. 100000) then
1563	print*,'pb moldiff_red infinite loop'
1564	print*,ig,nl,T(nl),pnew2,qnew(il,:),Znew2
1565	stop
1566	endif
1567
1568	! print*,'test',Pnew2,Znew2,Nnew(nq),pp(il)
1569	enddo
1570
1571	facP=(pp(il)-Pnew)/(Pnew2-Pnew)
1572
1573	! do nn=1,nq
1574	! qnew(il,nn)=dble(M(nn))*Nnew(nn)
1575	! & /sum(dble(M(:))*Nnew(:))
1576	! enddo
1577
1578	do nn=1,nq
1579	rhoknew(il,nn)=dble(M(nn))Nnew(nn)FacM(il,nn)
1580	enddo
1581	rhonew(il)=sum(rhoknew(il,:))
1582	do nn=1,nq
1583	qnew(il,nn)=rhoknew(il,nn)/rhonew(il)
1584	enddo
1585	tnew(il)=T(nl)
1586
1587	endif
1588	endif
1589	enddo
1590
1591	! write(,), ' -- Sortie moldiff_red -- '
1592
1593	END
1594

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

Download in other formats:

Original Format