Context Navigation

moldiff_red.F90 @ 535

Last change on this file since 535 was 517, checked in by emillour, 14 years ago
Mars GCM: Adding the new molecular diffusion routines (old "moldiff" and "moldiffcoeff" are kept for now and if further tests are needed) by JYC: moldiffcoeff_red.F, moldiff_red.F90, diffusion.h Also updated "phymars/comdiurn.h", "aeronomars/conc.h" and "aeronomars/chimiedata.h" to be fixed/free form compatible.
File size: 37.7 KB

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

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

Download in other formats:

Original Format