Context Navigation

moldiff_red.F90 @ 1357

Last change on this file since 1357 was 1357, checked in by emillour, 10 years ago

Mars GCM:

Update and bug correction on the escape fluxes of H and H2.

JYC

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

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

Download in other formats:

Original Format