Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

moldiff_red.F90 @ 3586

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

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

Download in other formats:

Original Format