Context Navigation

moldiff_red.F90 @ 1358

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

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

Download in other formats:

Original Format