Context Navigation

moldiff_red.F90 @ 2732

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

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

Download in other formats:

Original Format