Context Navigation

moldiff_red.F90 @ 1528

Last change on this file since 1528 was 1528, checked in by emillour, 9 years ago

Mars GCM:

Got rid of references to "dimensions.h" from physics packages: use nbp_lon (=iim), nbp_lat (==jjp1) and nbp_lev from module mod_grid_phy_lmdz (in phy_common) instead.
Added "ioipsl_getin_p_mod.F90" (getin_p routine) in phy_common to correctly read in parameters from *.def files in a parallel environment.

File size: 39.2 KB

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

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

Download in other formats:

Original Format