Context Navigation

moldiff_red.F90 @ 2808

Last change on this file since 2808 was 2611, checked in by romain.vande, 3 years ago
LMDZ_MARS RV : Open_MP; Reading files in parallel for moldiff parametrisation (callmoldiff = .true.)
File size: 40.6 KB

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

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

Download in other formats:

Original Format