Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

moldiff_red.F90 @ 3436

Last change on this file since 3436 was 3175, checked in by emillour, 11 months ago
Pluto PCM: Add the old Pluto LMDZ for reference (required prior step to making an LMDZ.PLUTO using the same framework as the other physics packages). TB+EM
File size: 39.8 KB

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

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

Download in other formats:

Original Format