Context Navigation

moldiff_MPF.F90 @ 3948

Last change on this file since 3948 was 3743, checked in by emillour, 6 months ago
Mars PCM: OpenMP bug fix: add some missing threadprivate statements in moldiff_MPF EM
File size: 42.5 KB

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

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

Download in other formats:

Original Format