Context Navigation

moldiff_red.F90 @ 3094

Last change on this file since 3094 was 2836, checked in by abierjon, 2 years ago

VENUS GCM:

INCLUDING THE IONOSPHERE CODE IN VENUS GCM

ATTENTION: INCREASED MEMORY DEMAND

NEEDS AT LEAST 135 GB of allocated memory

==============================================================
===> LIST OF APPENDED FILES AND INTERNAL ADDITIONS <==========
==============================================================

NEW MODEL SPECIES

deftank/traceur-chemistry-IONOSPHERE.def

Neutrals: 4 Species: N, NO, NO2, N(2D)

Ions: 15 species:

CO2+, CO+, O+, O2+, H+,

N2+, H2O+, OH+, C+, HCO+,

H3O+, HCO2+, N+, NO+, elec

FROM 36 to 55 chemical species

NEW KEYWORD OF PHYSIQ.DEF

deftank/physiq-96x96x90-chemistry-IONOSPHERE.def

ok_ionchem: keyword supposed to activate ion chemistry.

(be careful that n, no, no2, n2d and the ion species are in the deftank/traceur-chemistry-IONOSPHERE.def)

ok_jonline: keyword supposed to activate the online photochemistry

(be careful that n, no, no2 and n2d are in the traceur-chemistry-IONOSPHERE.def)

==============================================================
===> LIST OF MODIFIED PROGRAMS <=========================
==============================================================

nonoro_gwd_ran_mod.F90

Change EPFLUXMAX value from 5.E-3 to 1.E-3

photochemistry_venus.F90 (krates; photolysis_online; indices;)

Import of keywords: ok_ionchem, tuneupperatm & ok_jonline
addition of ion species in order
Forcing electroneutrality
Update of the reactions a001 and a002 with taking into account the other species
Change of formula for a002 with the formulation of Baulch et al., 1976 (confirmed by Smith and Robertson, 2008)

photolysis_mod.F90

Modification of the rdsolarflux subroutine to include interpolation with Atlas1 and Atlas3 in connection with E10.7

concentration2.F90

Added chemical species n2d, no, no2 and n in the conductivity calculation.

The ions have been excluded because their sum is 10^{5 times less dense than the neutrals and

their thermal conductivity is unknown}

iono_h.F90

Addition of the phdisrate routine
replace the electronic temperature of Mars by that of

origin = 1: Theis et al. 1980 (Venus) with bilinear interpolation altitude/cos(SZA)
origin = 2: Theis et al. 1984 (Venus) with the formula of the electronic temperature at high solar activity

addition of an ion temperature model based on VIRA

cleshphys.h

added ok_ionchem & ok_jonline in COMMON/clesphys_l/

conf_phys.f90

add ok_ionchem & ok_jonline as parameters to read from physiq.def file set to .false. by default

chemparam_mod.F90

Add species in M_tr and corresponding i_X. Set all i_X to zero before reading traceur-chemistry-IONOSPHERE.def
Added Type_tr table to differentiate species: 1 == neutral, 2 == ION, 3 == liquid, 10 == others

euvheat.F90; hrtherm.F; jthermcalc_e107.F; param_read_e107.F

Normalization with Mars

A.M

File size: 39.1 KB

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

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

Download in other formats:

Original Format