Context Navigation

moldiff_red.F90 @ 1431

Last change on this file since 1431 was 1431, checked in by aslmd, 10 years ago
LMDZ.MARS -- commented redundant call to comgeomphy and corrected call to areas in moldiff_red
File size: 39.5 KB

Line
1	subroutine moldiff_red(ngrid,nlayer,nq,pplay,pplev,pt,pdt,pq,pdq,&
2	ptimestep,zzlay,pdteuv,pdtconduc,pdqdiff)
3
4	use tracer_mod, only: noms, mmol
5	use comgeomfi_h, only: area
6
7	implicit none
8
9	!#include "dimensions.h"
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	! July 2014 JYC ADD BALISTIC Transport coupling to compute wup for H and H2
20
21
22
23	!
24	! Input/Output
25	!
26	integer,intent(in) :: ngrid ! number of atmospheric columns
27	integer,intent(in) :: nlayer ! number of atmospheric layers
28	integer,intent(in) :: nq ! number of advected tracers
29	real ptimestep
30	real pplay(ngrid,nlayer)
31	real zzlay(ngrid,nlayer)
32	real pplev(ngrid,nlayer+1)
33	real pq(ngrid,nlayer,nq)
34	real pdq(ngrid,nlayer,nq)
35	real pt(ngrid,nlayer)
36	real pdt(ngrid,nlayer)
37	real pdteuv(ngrid,nlayer)
38	real pdtconduc(ngrid,nlayer)
39	real pdqdiff(ngrid,nlayer,nq)
40	!
41	! Local
42	!
43
44
45	! real hco2(ncompdiff),ho
46
47	integer,dimension(nq) :: indic_diff
48	integer ig,iq,nz,l,k,n,nn,p,ij0
49	integer istep,il,gcn,ntime,nlraf
50	real*8 masse
51	real*8 masseU,kBolt,RgazP,Rmars,Grav,Mmars
52	real*8 rho0,D0,T0,H0,time0,dZ,time,dZraf,tdiff,Zmin,Zmax
53	real*8 FacEsc,invsgmu,PhiEscH,PhiEscH2
54	real*8 hp(nlayer)
55	real*8 pp(nlayer)
56	real*8 pint(nlayer)
57	real*8 tt(nlayer),tnew(nlayer),tint(nlayer)
58	real*8 zz(nlayer)
59	real*8,dimension(:,:),allocatable,save :: qq,qnew,qint,FacMass
60	real*8,dimension(:,:),allocatable,save :: rhoK,rhokinit
61	real*8 rhoT(nlayer)
62	real*8 dmmeandz(nlayer)
63	real*8 massemoy(nlayer)
64	real*8,dimension(:),allocatable :: Praf,Traf,Rraf,Mraf,Nraf,Draf,Hraf,Wraf
65	real*8,dimension(:),allocatable :: Zraf,Tdiffraf
66	real*8,dimension(:),allocatable :: Prafold,Mrafold
67	real*8,dimension(:,:),allocatable :: Qraf,Rrafk,Nrafk
68	real*8,dimension(:,:),allocatable :: Rrafkold
69	real*8,dimension(:,:),allocatable :: Drafmol,Hrafmol,Wrafmol,Tdiffrafmol
70	real*8,dimension(:),allocatable :: Atri,Btri,Ctri,Dtri,Xtri,Tad,Dad,Zad,rhoad
71	real*8,dimension(:),allocatable :: alpha,beta,gama,delta,eps
72	real*8,dimension(:),allocatable,save :: wi,Wad,Uthermal,Lambdaexo,Hspecie
73	real*8,dimension(:),allocatable,save :: Mtot1,Mtot2,Mraf1,Mraf2
74	character(len=20),dimension(14) :: ListeDiff
75	real*8,parameter :: pi=3.141592653
76	real*8,parameter :: g=3.72d0
77	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
78	! tracer numbering in the molecular diffusion
79	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
80	! We need the index of escaping species
81
82	integer,save :: i_h2
83	integer,save :: i_h
84	! vertical index limit for the molecular diffusion
85	integer,save :: il0
86
87	! Tracer indexes in the GCM:
88	! integer,save :: g_co2=0
89	! integer,save :: g_co=0
90	! integer,save :: g_o=0
91	! integer,save :: g_o1d=0
92	! integer,save :: g_o2=0
93	! integer,save :: g_o3=0
94	! integer,save :: g_h=0
95	! integer,save :: g_h2=0
96	! integer,save :: g_oh=0
97	! integer,save :: g_ho2=0
98	! integer,save :: g_h2o2=0
99	! integer,save :: g_n2=0
100	! integer,save :: g_ar=0
101	! integer,save :: g_h2o=0
102	! integer,save :: g_n=0
103	! integer,save :: g_no=0
104	! integer,save :: g_no2=0
105	! integer,save :: g_n2d=0
106	! integer,save :: g_oplus=0
107	! integer,save :: g_hplus=0
108
109	integer,save :: ncompdiff
110	integer,dimension(:),allocatable,save :: gcmind ! array of GCM indexes
111	integer ierr,compteur
112
113	logical,save :: firstcall=.true.
114	! real abfac(ncompdiff)
115	real,dimension(:,:),allocatable,save :: dij
116	real,save :: step
117
118	! Initializations at first call
119	if (firstcall) then
120
121	! Liste des especes qui sont diffuses si elles sont presentes
122	! ListeDiff=['co2','o','n2','ar','co','h2','h','d2','d','hd','o2','oh','ho2','h2o_vap','h2o2','o1d','n','no','no2']
123	ListeDiff(1)='co2'
124	ListeDiff(2)='o'
125	ListeDiff(3)='n2'
126	ListeDiff(4)='ar'
127	ListeDiff(5)='co'
128	ListeDiff(6)='h2'
129	ListeDiff(7)='h'
130	ListeDiff(8)='d2'
131	ListeDiff(9)='hd'
132	ListeDiff(10)='d'
133	ListeDiff(11)='o2'
134	ListeDiff(12)='h2o_vap'
135	ListeDiff(13)='o3'
136	ListeDiff(14)='n'
137	i_h=1000
138	i_h2=1000
139	! On regarde quelles especes diffusables sont presentes
140
141	ncompdiff=0
142	indic_diff(1:nq)=0
143
144	do nn=1,nq
145	do n=1,14
146	if (ListeDiff(n) .eq. noms(nn)) then
147	indic_diff(nn)=1
148	! if (noms(nn) .eq. 'h') i_h=n
149	! if (noms(nn) .eq. 'h2') i_h2=n
150	endif
151
152	enddo
153	if (indic_diff(nn) .eq. 1) then
154	print*,'specie ', noms(nn), 'diffused in moldiff_red'
155	ncompdiff=ncompdiff+1
156	endif
157	enddo
158	print*,'number of diffused species:',ncompdiff
159	allocate(gcmind(ncompdiff))
160
161	! Store gcm indexes in gcmind
162	n=0
163	do nn=1,nq
164	if (indic_diff(nn) .eq. 1) then
165	n=n+1
166	gcmind(n)=nn
167	if (noms(nn) .eq. 'h') i_h=n
168	if (noms(nn) .eq. 'h2') i_h2=n
169	endif
170	enddo
171
172	! print*,'gcmind',gcmind,i_h,i_h2
173
174	! find vertical index above which diffusion is computed
175
176	do l=1,nlayer
177	if (pplay(1,l) .gt. Pdiff) then
178	il0=l
179	endif
180	enddo
181	il0=il0+1
182	print*,'vertical index for diffusion',il0,pplay(1,il0)
183
184	allocate(dij(ncompdiff,ncompdiff))
185
186	call moldiffcoeff_red(dij,indic_diff,gcmind,ncompdiff)
187	print*,'MOLDIFF EXO'
188
189	! allocatation des tableaux dependants du nombre d especes diffusees
190	allocate(qq(nlayer,ncompdiff))
191	allocate(qnew(nlayer,ncompdiff))
192	allocate(qint(nlayer,ncompdiff))
193	allocate(FacMass(nlayer,ncompdiff))
194	allocate(rhok(nlayer,ncompdiff))
195	allocate(rhokinit(nlayer,ncompdiff))
196
197	allocate(wi(ncompdiff))
198	allocate(wad(ncompdiff))
199	allocate(uthermal(ncompdiff))
200	allocate(lambdaexo(ncompdiff))
201	allocate(Hspecie(ncompdiff))
202	allocate(Mtot1(ncompdiff))
203	allocate(Mtot2(ncompdiff))
204	allocate(Mraf1(ncompdiff))
205	allocate(Mraf2(ncompdiff))
206
207	firstcall= .false.
208	step=1
209	endif ! of if (firstcall)
210
211	!
212	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
213
214	masseU=1.660538782d-27
215	kbolt=1.3806504d-23
216	RgazP=8.314472
217	! PI =3.141592653
218	! g=3.72d0
219	Rmars=3390000d0 ! Used to compute escape parameter no need to be more accurate
220	Grav=6.67d-11
221	Mmars=6.4d23
222	ij0=6000 ! For test
223	invsgmu=1d0/g/masseU
224
225	! Compute the wup(ig) for H and H2 using the balistic code from R Yelle
226
227	PhiEscH=0D0
228	PhiEscH2=0D0
229
230	do ig=1,ngrid
231	pp=dble(pplay(ig,:))
232
233	! Update the temperature
234
235	! CALL TMNEW(pt(ig,:),pdt(ig,:),pdtconduc(ig,:),pdteuv(ig,:) &
236	! & ,tt,ptimestep,nlayer,ig)
237	do l=1,nlayer
238	tt(l)=pt(ig,l)1D0+(pdt(ig,l)dble(ptimestep)+ &
239	pdtconduc(ig,l)*dble(ptimestep)+ &
240	pdteuv(ig,l)*dble(ptimestep))
241	! to cach Nans...
242	if (tt(l).ne.tt(l)) then
243	print*,'Err TMNEW',ig,l,tt(l),pt(ig,l), &
244	pdt(ig,l),pdtconduc(ig,l),pdteuv(ig,l),dble(ptimestep)
245	endif
246	enddo ! of do l=1,nlayer
247
248	! Update the mass mixing ratios modified by other processes
249
250	! CALL QMNEW(pq(ig,:,:),pdq(ig,:,:),qq,ptimestep,nlayer, &
251	! & ncompdiff,gcmind,ig)
252	do iq=1,ncompdiff
253	do l=1,nlayer
254	qq(l,iq)=pq(ig,l,gcmind(iq))*1D0+( &
255	pdq(ig,l,gcmind(iq))*dble(ptimestep))
256	qq(l,iq)=max(qq(l,iq),1d-30)
257	enddo ! of do l=1,nlayer
258	enddo ! of do iq=1,ncompdiff
259
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,mmol,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),pdq(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,mmol,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,mmol,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(mmol(gcmind(nn)))/masseU
378	! enddo
379	! Mraf(l)=1d0/sum(Qraf(l,:)/dble(mmol(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(mmol(gcmind(nn))))
407	Lambdaexo(nn)=masseUdble(mmol(gcmind(nn)))Grav*Mmars/ &
408	& (Rmars+Zraf(nlraf))/kbolt/Traf(nlraf)
409	wi(nn)=0D0
410	if (nn .eq. i_h .or. nn .eq. i_h2) then
411	wi(nn)=Uthermal(nn)/2d0/sqrt(PI)exp(-Lambdaexo(nn)) &
412	& (Lambdaexo(nn)+1d0)
413	endif
414	enddo
415
416	! print*,'wi',wi(i_h),wi(i_h2),Uthermal,Lambdaexo,mmol(gcmind(:))
417	! print*,'wi',wi
418	! stop
419
420	! Compute time step for diffusion
421
422	! Loop on species
423
424	T0=Traf(nlraf)
425	rho0=1d0
426
427	do nn=1,ncompdiff
428	masse=dble(mmol(gcmind(nn)))
429	! DIffusion coefficient
430	CALL DCOEFF(nn,dij,Praf,Traf,Nraf,Nrafk,Draf,nlraf,ncompdiff)
431	Drafmol(:,nn)=Draf(:)
432	! Scale height of the density of the specie
433	CALL HSCALEREAL(nn,Nrafk,Dzraf,Hraf,nlraf,ncompdiff)
434	Hrafmol(:,nn)=Hraf(:)
435	! Hspecie(nn)=kboltT0/masseinvsgmu
436	! Computation of the diffusion vertical velocity of the specie
437	CALL VELVERT(nn,Traf,Hraf,Draf,Dzraf,masse,Wraf,nlraf)
438	Wrafmol(:,nn)=Wraf(:)
439	! Computation of the diffusion time
440	CALL TIMEDIFF(nn,Hraf,Wraf,Tdiffraf,nlraf)
441	Tdiffrafmol(:,nn)=Tdiffraf
442	enddo
443	! We use a lower time step
444	Tdiff=minval(Tdiffrafmol)
445	Tdiff=minval(Tdiffrafmol(nlraf,:))*Mraf(nlraf)
446	! Some problems when H is dominant
447	! The time step is chosen function of atmospheric mass at higher level
448	! It is not very nice
449
450	! if (ig .eq. ij0) then
451	! print*,'test',ig,tdiff,tdiffmin,minloc(Tdiffrafmol),minloc(Tdiffrafmol(nlraf,:))
452	! endif
453	if (tdiff .lt. tdiffminMraf(nlraf)) tdiff=tdiffminMraf(nlraf)
454
455	! Number of time step
456	ntime=anint(dble(ptimestep)/tdiff)
457	! print*,'ptime',ig,ptimestep,tdiff,ntime,tdiffmin,Mraf(nlraf)
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(mmol(gcmind(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(mmol(gcmind(nn)))/masseU
572	enddo
573	Mraf(l)=1d0/sum(Qraf(l,:)/dble(mmol(gcmind(:))))
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,mmol,gcmind,nlraf,ncompdiff,nlayer,FacMass,ig)
595
596	CALL RHOTOT(pp,tt,massemoy,qnew,RHOT,RHOK,nlayer,ncompdiff)
597
598	! Update total escape flux of H and H2 (if q was really qnew, but not forget we will output
599	! the trend only at the end
600
601	PhiEscH=PhiEscH+wi(i_h)Nrafk(nlraf,i_h)area(ig) ! in s-1
602	PhiEscH2=PhiEscH2+wi(i_h2)Nrafk(nlraf,i_h2)area(ig) ! in s-1 (U in m/s, aire in m2, Nrafk in m-3)
603	! print*,'test',ig,wi(i_h),Nrafk(nlraf,i_h),wi(i_h2),Nrafk(nlraf,i_h2),area(ig),PhiEscH,PhiEscH2,i_h,i_h2
604	! stop
605
606
607	if (ig .eq. ij0) then
608	do l=il0,nlayer
609	write(*,'(i2,1x,19(e12.4,1x))') l,zz(l),tt(l),RHOK(l,1)/sum(RHOK(l,:)),RHOKINIT(l,1)/sum(RHOKINIT(l,:)),&
610	& RHOK(l,2)/sum(RHOK(l,:)),RHOKINIT(l,2)/sum(RHOKINIT(l,:)),&
611	& RHOK(l,6)/sum(RHOK(l,:)),RHOKINIT(l,6)/sum(RHOKINIT(l,:)),&
612	& RHOK(l,5)/sum(RHOK(l,:)),RHOKINIT(l,5)/sum(RHOKINIT(l,:)),&
613	& RHOK(l,7)/sum(RHOK(l,:)),RHOKINIT(l,7)/sum(RHOKINIT(l,:))
614	enddo
615	endif
616
617	! Compute total mass of each specie on the GCM vertical grid
618
619	Mtot2(1:ncompdiff)=0d0
620
621	do l=il0,nlayer
622	do nn=1,ncompdiff
623	Mtot2(nn)=Mtot2(nn)+1d0/gqnew(l,nn) &
624	& (dble(pplev(ig,l))-dble(pplev(ig,l+1)))
625	enddo
626	enddo
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	! print*,'Escape flux H, H2 (s-1)',PhiEscH,PhiEscH2
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 tracer_mod, only: nqmx
938	IMPLICIT NONE
939	!#include "dimensions.h"
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,nqmx) :: 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,mmol,qq,gc,nl,nq)
973	use tracer_mod, only: nqmx
974	IMPLICIT NONE
975	!#include "dimensions.h"
976
977	INTEGER :: nl,nq,l
978	INTEGER,dimension(nq) :: gc
979	REAL*8,DIMENSION(nl,nq) :: qq
980	REAL*8,DIMENSION(nl) :: massemoy
981	REAL,DIMENSION(nqmx) :: MMOL
982
983
984	do l=1,nl
985	massemoy(l)=1D0/sum(qq(l,:)/dble(mmol(gc(:))))
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	INTEGER :: nl,l,ig
1007	REAL*8,dimension(nl) :: P,T,M,Z,H
1008	REAL*8 :: z0
1009	REAL*8 :: kbolt,masseU,Konst,g,Hpm
1010	masseU=1.660538782d-27
1011	kbolt=1.3806504d-23
1012	Konst=kbolt/masseU
1013	g=3.72D0
1014
1015	z0=0d0
1016	Z(1)=z0
1017	H(1)=Konst*T(1)/M(1)/g
1018
1019	do l=2,nl
1020	H(l)=Konst*T(l)/M(l)/g
1021	Hpm=H(l-1)
1022	Z(l)=z(l-1)-Hpm*log(P(l)/P(l-1))
1023	if (Z(l) .ne. Z(l)) then
1024	print*,'pb',l,ig
1025	print*,'P',P
1026	print*,'T',T
1027	print*,'M',M
1028	print*,'Z',Z
1029	print*,'Hpm',Hpm
1030	endif
1031	enddo
1032
1033	END
1034
1035	SUBROUTINE RHOTOT(P,T,M,qq,rhoN,rhoK,nl,nq)
1036	IMPLICIT NONE
1037
1038	REAL*8 :: kbolt,masseU,Konst
1039	INTEGER :: nl,nq,l,iq
1040	REAL*8,DIMENSION(nl) :: P,T,M,RHON
1041	REAL*8,DIMENSION(nl,nq) :: RHOK,qq
1042	masseU=1.660538782d-27
1043	kbolt=1.3806504d-23
1044	Konst=Kbolt/masseU
1045
1046	do l=1,nl
1047	RHON(l)=P(l)*M(l)/T(l)/Konst
1048	do iq=1,nq
1049	RHOK(l,iq)=qq(l,iq)*RHON(l)
1050	enddo
1051	enddo
1052
1053	END
1054
1055	SUBROUTINE UPPER_RESOL(P,T,Z,M,R,Rk, &
1056	& qq,mmol,gc,Praf,Traf,Qraf,Mraf,Zraf, &
1057	& Nraf,Nrafk,Rraf,Rrafk,il,nl,nq,nlx,ig)
1058	use tracer_mod, only: nqmx
1059	IMPLICIT NONE
1060	!#include "dimensions.h"
1061
1062	INTEGER :: nl,nq,il,l,i,iq,nlx,iz,ig
1063	INTEGER :: gc(nq)
1064	INTEGER,DIMENSION(1) :: indz
1065	REAL*8, DIMENSION(nlx) :: P,T,Z,M,R
1066	REAL*8, DIMENSION(nlx,nq) :: qq,Rk
1067	REAL*8, DIMENSION(nl) :: Praf,Traf,Mraf,Zraf,Nraf,Rraf
1068	REAL*8 :: kbolt,masseU,Konst,g
1069	REAL*8, DIMENSION(nl,nq) :: Qraf,Rrafk,Nrafk
1070	REAL*8 :: facZ,dZ,H
1071	REAL,DIMENSION(nqmx) :: mmol
1072	masseU=1.660538782d-27
1073	kbolt=1.3806504d-23
1074	Konst=Kbolt/masseU
1075	g=3.72d0
1076
1077
1078	Zraf(1)=z(il)
1079	Praf(1)=P(il)
1080	Traf(1)=T(il)
1081	Nraf(1)=Praf(1)/kbolt/Traf(1)
1082	do iq=1,nq
1083	Qraf(1,iq)=qq(il,iq)
1084	enddo
1085	Mraf(1)=1d0/sum(Qraf(1,:)/dble(mmol(gc(:))))
1086	Rraf(1)=Mraf(1)masseUNraf(1)
1087	do iq=1,nq
1088	Rrafk(1,iq)=Rraf(1)*Qraf(1,iq)
1089	Nrafk(1,iq)=Rrafk(1,iq)/masseU/dble(mmol(gc(iq)))
1090	enddo
1091	Zraf(nl)=z(nlx)
1092
1093	do l=2,nl-1
1094	Zraf(l)=Zraf(1)+(Zraf(nl)-Zraf(1))/dble(nl-1)*dble((l-1))
1095	indz=maxloc(z,mask=z <= Zraf(l))
1096	iz=indz(1)
1097	if (iz .lt. 1 .or. iz .gt. nlx) then
1098	print*,'danger',iz,nl,Zraf(l),l,Zraf(1),Zraf(nl),z,P,T,ig
1099	stop
1100	endif
1101	dZ=Zraf(l)-Zraf(l-1)
1102	! dZ=Zraf(l)-z(iz)
1103	facz=(Zraf(l)-z(iz))/(z(iz+1)-z(iz))
1104	Traf(l)=T(iz)+(T(iz+1)-T(iz))*facz
1105	do iq=1,nq
1106	! Qraf(l,iq)=qq(iz,iq)+(qq(iz+1,iq)-qq(iz,iq))*facz
1107	Rrafk(l,iq)=Rk(iz,iq)+(Rk(iz+1,iq)-Rk(iz,iq))*facZ
1108	Rrafk(l,iq)=Rk(iz,iq)(Rk(iz+1,iq)/Rk(iz,iq))*facZ
1109	enddo
1110	! Mraf(l)=1D0/(sum(qraf(l,:)/dble(mmol(gc(:)))))
1111	Rraf(l)=sum(Rrafk(l,:))
1112	do iq=1,nq
1113	Qraf(l,iq)=Rrafk(l,iq)/Rraf(l)
1114	enddo
1115	Mraf(l)=1D0/(sum(qraf(l,:)/dble(mmol(gc(:)))))
1116	! H=Konst*Traf(l)/Mraf(l)/g
1117	! H=Konst*T(iz)/M(iz)/g
1118	! Praf(l)=P(iz)*exp(-dZ/H)
1119	! Praf(l)=Praf(l-1)*exp(-dZ/H)
1120	! print,'iz',l,iz,Praf(il-1)exp(-dZ/H),z(iz),z(iz+1),H
1121	Nraf(l)=Rraf(l)/Mraf(l)/masseU
1122	Praf(l)=Nraf(l)kboltTraf(l)
1123	! Rraf(l)=Nraf(l)Mraf(l)masseU
1124	do iq=1,nq
1125	! Rrafk(l,iq)=Rraf(l)*Qraf(l,iq)
1126	Nrafk(l,iq)=Rrafk(l,iq)/dble(mmol(gc(iq)))/masseU
1127	if (Nrafk(l,iq) .lt. 0. .or. &
1128	& Nrafk(l,iq) .ne. Nrafk(l,iq)) then
1129	print*,'pb interpolation',l,iq,Nrafk(l,iq),Rrafk(l,iq), &
1130	& Qraf(l,iq),Rk(iz,iq),Rk(iz+1,iq),facZ,Zraf(l),z(iz)
1131	stop
1132	endif
1133	enddo
1134	enddo
1135	Zraf(nl)=z(nlx)
1136	Traf(nl)=T(nlx)
1137	do iq=1,nq
1138	Rrafk(nl,iq)=Rk(nlx,iq)
1139	Qraf(nl,iq)=Rk(nlx,iq)/R(nlx)
1140	Nrafk(nl,iq)=Rk(nlx,iq)/dble(mmol(gc(iq)))/masseU
1141	enddo
1142	Nraf(nl)=sum(Nrafk(nl,:))
1143	Praf(nl)=Nraf(nl)kboltTraf(nl)
1144	Mraf(nl)=1D0/sum(Qraf(nl,:)/dble(mmol(gc(:))))
1145	END
1146
1147	SUBROUTINE CORRMASS(qq,qint,FacMass,nl,nq)
1148	IMPLICIT NONE
1149	INTEGER :: nl,nq,l,nn
1150	REAL*8,DIMENSION(nl,nq) :: qq,qint,FacMass
1151
1152	do nn=1,nq
1153	do l=1,nl
1154	FacMass(l,nn)=qq(l,nn)/qint(l,nn)
1155	enddo
1156	enddo
1157
1158	END
1159
1160
1161	SUBROUTINE DCOEFF(nn,Dij,P,T,N,Nk,D,nl,nq)
1162	IMPLICIT NONE
1163	REAL*8,DIMENSION(nl) :: N,T,D,P
1164	REAL*8,DIMENSION(nl,nq) :: Nk
1165	INTEGER :: nn,nl,nq,l,iq
1166	REAL,DIMENSION(nq,nq) :: Dij
1167	REAL*8 :: interm,P0,T0,ptfac,dfac
1168
1169	P0=1D5
1170	T0=273d0
1171
1172
1173	do l=1,nl
1174	ptfac=(P0/P(l))(T(l)/T0)*1.75d0
1175	D(l)=0d0
1176	interm=0d0
1177	do iq=1,nq
1178	if (iq .ne. nn) then
1179	dfac=dble(dij(nn,iq))*ptfac
1180	interm=interm+Nk(l,iq)/N(l)/dfac
1181	endif
1182	enddo
1183	D(l)=1d0/interm
1184	enddo
1185	END
1186
1187	SUBROUTINE HSCALEREAL(nn,Nk,Dz,H,nl,nq)
1188	IMPLICIT NONE
1189	INTEGER :: nn,nl,nq,l
1190	REAL*8,DIMENSION(nl) :: H
1191	REAL*8,DIMENSION(nl,nq) :: Nk
1192	REAL*8 :: Dz
1193
1194	H(1)=(-3D0Nk(1,nn)+4d0NK(2,nn)-Nk(3,nn))/(2D0*DZ)/ &
1195	& NK(1,nn)
1196
1197	H(1)=-1D0/H(1)
1198
1199	DO l=2,nl-1
1200	H(l)=(Nk(l+1,nn)-NK(l-1,nn))/(2D0*DZ)/NK(l,nn)
1201	H(l)=-1D0/H(l)
1202	ENDDO
1203
1204	H(nl)=(3D0Nk(nl,nn)-4D0Nk(nl-1,nn)+Nk(nl-2,nn))/(2D0*DZ)/ &
1205	& Nk(nl,nn)
1206	H(nl)=-1D0/H(nl)
1207
1208	! do l=1,nl
1209	! if (abs(H(l)) .lt. 100.) then
1210	! print*,'H',l,H(l),Nk(l,nn),nn
1211	! endif
1212	! enddo
1213
1214	END
1215
1216	SUBROUTINE VELVERT(nn,T,H,D,Dz,masse,W,nl)
1217	IMPLICIT NONE
1218	INTEGER :: l,nl,nn
1219	REAL*8,DIMENSION(nl) :: T,H,D,W,DT
1220	REAL*8 :: Dz,Hmol,masse
1221	REAL*8 :: kbolt,masseU,Konst,g
1222	masseU=1.660538782d-27
1223	kbolt=1.3806504d-23
1224	Konst=Kbolt/masseU
1225	g=3.72d0
1226
1227	DT(1)=1D0/T(1)(-3D0T(1)+4D0T(2)-T(3))/(2D0DZ)
1228	Hmol=Konst*T(1)/masse/g
1229	W(1)=-D(1)*(1D0/H(1)-1D0/Hmol-DT(1))
1230
1231	DO l=2,nl-1
1232	DT(l)=1D0/T(l)(T(l+1)-T(l-1))/(2D0DZ)
1233	Hmol=Konst*T(l)/masse/g
1234	W(l)=-D(l)*(1D0/H(l)-1D0/Hmol-DT(l))
1235	ENDDO
1236
1237	DT(nl)=1D0/T(nl)(3d0T(nl)-4D0T(nl-1)+T(nl-2))/(2D0DZ)
1238	Hmol=Konst*T(nl)/masse/g
1239	W(nl)=-D(nl)*(1D0/H(nl)-1D0/Hmol-DT(nl))
1240
1241	! do l=1,nl
1242	! print*,'W',W(l),D(l),H(l),DT(l)
1243	! enddo
1244
1245	END
1246
1247	SUBROUTINE TIMEDIFF(nn,H,W,TIME,nl)
1248	IMPLICIT NONE
1249	INTEGER :: nn,nl,l
1250	REAL*8,DIMENSION(nl) :: W,H,TIME
1251
1252	DO l=1,nl
1253	TIME(l)=abs(H(l)/W(l))
1254	if (TIME(l) .lt. 1.D-4) then
1255	! print*,'low tdiff',H(l),W(l),nn,l
1256	endif
1257	ENDDO
1258
1259	END
1260
1261
1262	SUBROUTINE DIFFPARAM(T,D,dz,alpha,beta,gama,delta,eps,nl)
1263	IMPLICIT NONE
1264	INTEGER :: nl,l
1265	REAL*8,DIMENSION(nl) :: T,D
1266	REAL*8 :: DZ,DZinv
1267	REAL*8,DIMENSION(nl) :: alpha,beta,gama,delta,eps
1268
1269	! Compute alpha,beta and delta
1270	! lower altitude values
1271	DZinv=1d0/(2D0*DZ)
1272
1273	beta(1)=1d0/T(1)
1274	alpha(1)=beta(1)(-3D0T(1)+4D0T(2)-T(3))Dzinv
1275	delta(1)=(-3D0D(1)+4D0D(2)-D(3))*Dzinv
1276
1277	beta(2)=1d0/T(2)
1278	alpha(2)=beta(2)(T(3)-T(1))Dzinv
1279	delta(2)=(D(3)-D(1))*Dzinv
1280
1281	! do l=2,nl-1
1282	! beta(l)=1D0/T(l)
1283	! alpha(l)=beta(l)(T(l+1)-T(l-1))Dzinv
1284	! delta(l)=(D(l+1)-D(l-1))*Dzinv
1285	! end do
1286
1287	do l=3,nl-1
1288	beta(l)=1D0/T(l)
1289	alpha(l)=beta(l)(T(l+1)-T(l-1))Dzinv
1290	delta(l)=(D(l+1)-D(l-1))*Dzinv
1291	gama(l-1)=(beta(l)-beta(l-2))*Dzinv
1292	eps(l-1)=(alpha(l)-alpha(l-2))*Dzinv
1293	enddo
1294
1295	! Upper altitude values
1296
1297	beta(nl)=1D0/T(nl)
1298	alpha(nl)=beta(nl)(3D0T(nl)-4D0T(nl-1)+T(nl-2))Dzinv
1299	delta(nl)=(3D0D(nl)-4D0D(nl-1)+D(nl-2))*Dzinv
1300
1301	! Compute the gama and eps coefficients
1302	! Lower altitude values
1303
1304	gama(1)=(-3D0beta(1)+4D0beta(2)-beta(3))*Dzinv
1305	eps(1)=(-3D0alpha(1)+4D0alpha(2)-alpha(3))*Dzinv
1306
1307	gama(nl-1)=(beta(nl)-beta(nl-2))*Dzinv
1308	eps(nl-1)=(alpha(nl)-alpha(nl-2))*Dzinv
1309
1310	! do l=2,nl-1
1311	! gama(l)=(beta(l+1)-beta(l-1))*Dzinv
1312	! eps(l)=(alpha(l+1)-alpha(l-1))*Dzinv
1313	! end do
1314
1315	gama(nl)=(3D0beta(nl)-4D0beta(nl-1)+beta(nl-2))*Dzinv
1316	eps(nl)=(3D0alpha(nl)-4D0alpha(nl-1)+alpha(nl-2))*Dzinv
1317
1318	! do l=1,nl
1319	! print*,'test diffparam',alpha(l),beta(l),delta(l),gama(l),eps(l)
1320	! enddo
1321	! stop
1322
1323	END
1324
1325
1326	SUBROUTINE MATCOEFF(alpha,beta,gama,delta,eps,Dad,rhoad, &
1327	& FacEsc,dz,dt,A,B,C,D,nl)
1328	IMPLICIT NONE
1329	INTEGER :: nl,l
1330	REAL*8, DIMENSION(nl) :: alpha,beta,gama,delta,eps,Dad,RHoad
1331	REAL*8 :: dz,dt,del1,del2,del3,FacEsc
1332	REAL*8, DIMENSION(nl) :: A,B,C,D
1333	del1=dt/2d0/dz
1334	del2=dt/dz/dz
1335	del3=dt
1336
1337	! lower boundary coefficients no change
1338	A(1)=0d0
1339	B(1)=1d0
1340	C(1)=0d0
1341	D(1)=rhoAd(1)
1342
1343	do l=2,nl-1
1344	A(l)=(delta(l)+(alpha(l)+beta(l))Dad(l))del1-Dad(l)*del2
1345	B(l)=-(delta(l)(alpha(l)+beta(l))+Dad(l)(gama(l)+eps(l)))*del3
1346	B(l)=B(l)+1d0+2d0Dad(l)del2
1347	C(l)=-(delta(l)+(alpha(l)+beta(l))Dad(l))del1-Dad(l)*del2
1348	D(l)=rhoAD(l)
1349	enddo
1350
1351
1352	! Upper boundary coefficients Diffusion profile
1353	C(nl)=0d0
1354	B(nl)=-1d0
1355	A(nl)=exp(-dZ(alpha(nl)+beta(nl)))FacEsc
1356	D(nl)=0D0
1357
1358
1359	END
1360
1361	SUBROUTINE Checkmass(X,Y,nl,nn)
1362	IMPLICIT NONE
1363
1364	INTEGER :: nl,nn
1365	REAL*8,DIMENSION(nl) :: X,Y
1366	REAL*8 Xtot,Ytot
1367
1368	Xtot=sum(X)
1369	Ytot=sum(Y)
1370
1371	if (abs((Xtot-Ytot)/Xtot) .gt. 1d-3) then
1372	print*,'no conservation for mass',Xtot,Ytot,nn
1373	endif
1374	END
1375
1376	SUBROUTINE Checkmass2(qold,qnew,P,il,nl,nn,nq)
1377	IMPLICIT NONE
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=3.72d0
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-2) THEN
1393	Print*,'We dont conserve mas',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 tracer_mod, only: nqmx
1401	IMPLICIT NONE
1402	!#include "dimensions.h"
1403	INTEGER :: nl,nq,nlx,il,nn,iP,ig,compteur
1404	INTEGER,DIMENSION(1) :: indP
1405	INTEGER,DIMENSION(nq) :: gc
1406	REAL*8,DIMENSION(nl) :: Z,P,T
1407	REAL*8,DIMENSION(nl,nq) :: Q,Nk,Rk
1408	REAL,DIMENSION(nqmx) :: M
1409	REAL*8,DIMENSION(nq) :: nNew
1410	REAL*8,DIMENSION(nlx) :: pp,tt,tnew
1411	REAL*8,DIMENSION(nlx) :: rhonew
1412	REAL*8,DIMENSION(nlx,nq) :: qq,qnew,rhoknew
1413	REAL*8 :: kbolt,masseU,Konst,g,Dz,facP,Hi
1414	REAL*8 :: Znew,Znew2,Pnew,Pnew2
1415	masseU=1.660538782d-27
1416	kbolt=1.3806504d-23
1417	Konst=Kbolt/masseU
1418	g=3.72d0
1419	Dz=Z(2)-Z(1)
1420	Znew=Z(nl)
1421	Znew2=Znew+dz
1422	! print*,'dz',Znew,Znew2,dz
1423	nNew(1:nq)=Nk(nl,1:nq)
1424	Pnew=P(nl)
1425
1426	do il=1,nlx
1427	! print*,'il',il,pp(il),P(1),P(nl)
1428	if (pp(il) .ge. P(1)) then
1429	qnew(il,:)=qq(il,:)
1430	tnew(il)=tt(il)
1431	endif
1432	if (pp(il) .lt. P(1)) then
1433	if (pp(il) .gt. P(nl)) then
1434	indP=maxloc(P,mask=P < pp(il))
1435	iP=indP(1)-1
1436	if (iP .lt. 1 .or. iP .gt. nl) then
1437	print*,'danger 2',iP,nl,pp(il)
1438	endif
1439	facP=(pp(il)-P(ip))/(P(ip+1)-P(ip))
1440	! print*,'P',P(ip),P(ip+1),facP,indP,iP
1441
1442	! do nn=1,nq
1443	! qnew(il,nn)=Q(iP,nn)+
1444	! & (Q(ip+1,nn)-Q(ip,nn))*facP
1445	! enddo
1446
1447	do nn=1,nq
1448	rhoknew(il,nn)=Rk(iP,nn)+ &
1449	& (Rk(ip+1,nn)-Rk(ip,nn))*facP
1450	enddo
1451	tnew(il)=T(iP)+(T(iP+1)-T(iP))*facP
1452	rhonew(il)=sum(rhoknew(il,:))
1453	do nn=1,nq
1454	qnew(il,nn)=rhoknew(il,nn)/rhonew(il)
1455	enddo
1456
1457	else ! pp < P(nl) need to extrapolate density of each specie
1458	Pnew2=Pnew
1459
1460	compteur=0
1461	do while (pnew2 .ge. pp(il))
1462	compteur=compteur+1
1463	do nn=1,nq
1464	Hi=Konst*T(nl)/dble(M(gc(nn)))/g
1465	Nnew(nn)=Nnew(nn)*exp(-dZ/Hi)
1466	enddo
1467	Pnew=Pnew2
1468	Pnew2=kboltT(nl)sum(Nnew(:))
1469	Znew=Znew2
1470	Znew2=Znew2+Dz
1471	if (compteur .ge. 100000) then
1472	print*,'error moldiff_red infinite loop'
1473	print*,ig,il,pp(il),tt(nl),pnew2,qnew(il,:),Znew2
1474	stop
1475	endif
1476	! print*,'test',Pnew2,Znew2,Nnew(nq),pp(il)
1477	enddo
1478
1479	facP=(pp(il)-Pnew)/(Pnew2-Pnew)
1480
1481	! do nn=1,nq
1482	! qnew(il,nn)=dble(M(gc(nn)))*Nnew(nn)
1483	! & /sum(dble(M(gc(:)))*Nnew(:))
1484	! enddo
1485
1486	do nn=1,nq
1487	rhoknew(il,nn)=dble(M(gc(nn)))*Nnew(nn)
1488	enddo
1489	rhonew(il)=sum(rhoknew(il,:))
1490	do nn=1,nq
1491	qnew(il,nn)=rhoknew(il,nn)/rhonew(il)
1492	enddo
1493	tnew(il)=T(nl)
1494	endif
1495	endif
1496	enddo
1497
1498	END
1499
1500	SUBROUTINE GCMGRID_P2(Z,P,Q,T,Nk,Rk,qq,qnew,tt,tnew &
1501	& ,pp,M,gc,nl,nq,nlx,facM,ig)
1502	use tracer_mod, only: nqmx
1503	IMPLICIT NONE
1504	!#include "dimensions.h"
1505	INTEGER :: nl,nq,nlx,il,nn,iP,ig,compteur
1506	INTEGER,DIMENSION(1) :: indP
1507	INTEGER,DIMENSION(nq) :: gc
1508	REAL*8,DIMENSION(nl) :: Z,P,T
1509	REAL*8,DIMENSION(nl,nq) :: Q,Nk,Rk
1510	REAL,DIMENSION(nqmx) :: M
1511	REAL*8,DIMENSION(nq) :: nNew
1512	REAL*8,DIMENSION(nlx) :: pp,rhonew,tt,tnew
1513	REAL*8,DIMENSION(nlx,nq) :: qq,qnew,facM,rhoknew
1514	REAL*8 :: kbolt,masseU,Konst,g,Dz,facP,Hi
1515	REAL*8 :: Znew,Znew2,Pnew,Pnew2
1516	masseU=1.660538782d-27
1517	kbolt=1.3806504d-23
1518	Konst=Kbolt/masseU
1519	g=3.72d0
1520	Dz=Z(2)-Z(1)
1521	Znew=Z(nl)
1522	Znew2=Znew+dz
1523	! print*,'dz',Znew,Znew2,dz
1524	nNew(1:nq)=Nk(nl,1:nq)
1525	Pnew=P(nl)
1526
1527	do il=1,nlx
1528	! print*,'il',il,pp(il),P(1),P(nl)
1529	if (pp(il) .ge. P(1)) then
1530	qnew(il,:)=qq(il,:)
1531	tnew(il)=tt(il)
1532	endif
1533	if (pp(il) .lt. P(1)) then
1534	if (pp(il) .gt. P(nl)) then
1535	indP=maxloc(P,mask=P < pp(il))
1536	iP=indP(1)-1
1537	if (iP .lt. 1 .or. iP .gt. nl) then
1538	print*,'danger 3',iP,nl,pp(il)
1539	endif
1540	facP=(pp(il)-P(ip))/(P(ip+1)-P(ip))
1541	! print*,'P',P(ip),P(ip+1),facP,indP,iP
1542
1543	! do nn=1,nq
1544	! qnew(il,nn)=Q(iP,nn)+
1545	! & (Q(ip+1,nn)-Q(ip,nn))*facP
1546	! enddo
1547
1548	do nn=1,nq
1549	rhoknew(il,nn)=(RK(iP,nn)+ &
1550	& (RK(iP+1,nn)-Rk(iP,nn))facP)facM(il,nn)
1551	enddo
1552	tnew(il)=T(iP)+(T(ip+1)-T(iP))*facP
1553	rhonew(il)=sum(rhoknew(il,:))
1554	do nn=1,nq
1555	qnew(il,nn)=rhoknew(il,nn)/rhonew(il)
1556	enddo
1557
1558	else ! pp < P(nl) need to extrapolate density of each specie
1559	Pnew2=Pnew
1560
1561	compteur=0
1562	do while (pnew2 .ge. pp(il))
1563	compteur=compteur+1
1564	do nn=1,nq
1565	Hi=Konst*T(nl)/dble(M(gc(nn)))/g
1566	Nnew(nn)=Nnew(nn)*exp(-dZ/Hi)
1567	enddo
1568	Pnew=Pnew2
1569	Pnew2=kboltT(nl)sum(Nnew(:))
1570	Znew=Znew2
1571	Znew2=Znew2+Dz
1572	if (compteur .ge. 100000) then
1573	print*,'pb moldiff_red infinite loop'
1574	print*,ig,nl,T(nl),pnew2,qnew(il,:),Znew2
1575	stop
1576	endif
1577
1578	! print*,'test',Pnew2,Znew2,Nnew(nq),pp(il)
1579	enddo
1580
1581	facP=(pp(il)-Pnew)/(Pnew2-Pnew)
1582
1583	! do nn=1,nq
1584	! qnew(il,nn)=dble(M(gc(nn)))*Nnew(nn)
1585	! & /sum(dble(M(gc(:)))*Nnew(:))
1586	! enddo
1587
1588	do nn=1,nq
1589	rhoknew(il,nn)=dble(M(gc(nn)))Nnew(nn)FacM(il,nn)
1590	enddo
1591	rhonew(il)=sum(rhoknew(il,:))
1592	do nn=1,nq
1593	qnew(il,nn)=rhoknew(il,nn)/rhonew(il)
1594	enddo
1595	tnew(il)=T(nl)
1596
1597	endif
1598	endif
1599	enddo
1600
1601	END
1602

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

Download in other formats:

Original Format