Context Navigation

moldiff_red.F90 @ 552

Last change on this file since 552 was 552, checked in by emillour, 13 years ago

Mars GCM:

bug fix in "aeronomars/moldiff_red.F90" (wrong mmol array size in included

subroutine) + adapted extreme limit test for H to altitudes of ~4000km
(compatble with 50 layer model).

bug fix in "nirco2abs.F" index of "co2" and "o" tracers were hard coded as

1 and 3!

updates from FGG of euvheat.F, callkeys.h and inifis.F to have the

"euveff" parameter read from callphys.def

File size: 37.8 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 > 4000 km there is a problem / stop
317
318	if (Zmax .gt. 4000000.) 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	#include "dimensions.h"
921
922	INTEGER :: nl,nq
923	INTEGER :: l,iq,ig
924	INTEGER,dimension(nq) :: gc
925	REAL,DIMENSION(nl,nqmx) :: Q1,DQ
926	REAL*8,DIMENSION(nl,nq) :: Q2
927	REAL :: dtime
928	DO l=1,nl
929	DO iq=1,nq
930	Q2(l,iq)=Q1(l,gc(iq))1D0+(DQ(l,gc(iq))dble(dtime))*1D0
931	Q2(l,iq)=max(Q2(l,iq),1d-30)
932	ENDDO
933	ENDDO
934	END
935
936	SUBROUTINE HSCALE(p,hp,nl)
937	IMPLICIT NONE
938
939	INTEGER :: nl
940	INTEGER :: l
941	REAL*8,dimension(nl) :: P
942	REAL*8,DIMENSION(nl) :: Hp
943
944	hp(1)=-log(P(2)/P(1))
945	hp(nl)=-log(P(nl)/P(nl-1))
946
947	DO l=2,nl-1
948	hp(l)=-log(P(l+1)/P(l-1))
949	ENDDO
950	END
951
952	SUBROUTINE MMOY(massemoy,mmol,qq,gc,nl,nq)
953	IMPLICIT NONE
954	#include "dimensions.h"
955
956	INTEGER :: nl,nq,l
957	INTEGER,dimension(nq) :: gc
958	REAL*8,DIMENSION(nl,nq) :: qq
959	REAL*8,DIMENSION(nl) :: massemoy
960	REAL,DIMENSION(nqmx) :: MMOL
961
962
963	do l=1,nl
964	massemoy(l)=1D0/sum(qq(l,:)/dble(mmol(gc(:))))
965	enddo
966
967	END
968
969	SUBROUTINE DMMOY(M,H,DM,nl)
970	IMPLICIT NONE
971	INTEGER :: nl,l
972	REAL*8,DIMENSION(nl) :: M,H,DM
973
974	DM(1)=(-3D0M(1)+4D0M(2)-M(3))/2D0/H(1)
975	DM(nl)=(3D0M(nl)-4D0M(nl-1)+M(nl-2))/2D0/H(nl)
976
977	do l=2,nl-1
978	DM(l)=(M(l+1)-M(l-1))/H(l)
979	enddo
980
981	END
982
983	SUBROUTINE ZVERT(P,T,M,Z,nl,ig)
984	IMPLICIT NONE
985	INTEGER :: nl,l,ig
986	REAL*8,dimension(nl) :: P,T,M,Z,H
987	REAL*8 :: z0
988	REAL*8 :: kbolt,masseU,Konst,g,Hpm
989	masseU=1.660538782d-27
990	kbolt=1.3806504d-23
991	Konst=kbolt/masseU
992	g=3.72D0
993
994	z0=0d0
995	Z(1)=z0
996	H(1)=Konst*T(1)/M(1)/g
997
998	do l=2,nl
999	H(l)=Konst*T(l)/M(l)/g
1000	Hpm=H(l-1)
1001	Z(l)=z(l-1)-Hpm*log(P(l)/P(l-1))
1002	if (Z(l) .ne. Z(l)) then
1003	print*,'pb',l,ig
1004	print*,'P',P
1005	print*,'T',T
1006	print*,'M',M
1007	print*,'Z',Z
1008	print*,'Hpm',Hpm
1009	endif
1010	enddo
1011
1012	END
1013
1014	SUBROUTINE RHOTOT(P,T,M,qq,rhoN,rhoK,nl,nq)
1015	IMPLICIT NONE
1016
1017	REAL*8 :: kbolt,masseU,Konst
1018	INTEGER :: nl,nq,l,iq
1019	REAL*8,DIMENSION(nl) :: P,T,M,RHON
1020	REAL*8,DIMENSION(nl,nq) :: RHOK,qq
1021	masseU=1.660538782d-27
1022	kbolt=1.3806504d-23
1023	Konst=Kbolt/masseU
1024
1025	do l=1,nl
1026	RHON(l)=P(l)*M(l)/T(l)/Konst
1027	do iq=1,nq
1028	RHOK(l,iq)=qq(l,iq)*RHON(l)
1029	enddo
1030	enddo
1031
1032	END
1033
1034	SUBROUTINE UPPER_RESOL(P,T,Z,M,R,Rk, &
1035	& qq,mmol,gc,Praf,Traf,Qraf,Mraf,Zraf, &
1036	& Nraf,Nrafk,Rraf,Rrafk,il,nl,nq,nlx,ig)
1037	IMPLICIT NONE
1038	#include "dimensions.h"
1039
1040	INTEGER :: nl,nq,il,l,i,iq,nlx,iz,ig
1041	INTEGER :: gc(nq)
1042	INTEGER,DIMENSION(1) :: indz
1043	REAL*8, DIMENSION(nlx) :: P,T,Z,M,R
1044	REAL*8, DIMENSION(nlx,nq) :: qq,Rk
1045	REAL*8, DIMENSION(nl) :: Praf,Traf,Mraf,Zraf,Nraf,Rraf
1046	REAL*8 :: kbolt,masseU,Konst,g
1047	REAL*8, DIMENSION(nl,nq) :: Qraf,Rrafk,Nrafk
1048	REAL*8 :: facZ,dZ,H
1049	REAL,DIMENSION(nqmx) :: mmol
1050	masseU=1.660538782d-27
1051	kbolt=1.3806504d-23
1052	Konst=Kbolt/masseU
1053	g=3.72d0
1054
1055
1056	Zraf(1)=z(il)
1057	Praf(1)=P(il)
1058	Traf(1)=T(il)
1059	Nraf(1)=Praf(1)/kbolt/Traf(1)
1060	do iq=1,nq
1061	Qraf(1,iq)=qq(il,iq)
1062	enddo
1063	Mraf(1)=1d0/sum(Qraf(1,:)/dble(mmol(gc(:))))
1064	Rraf(1)=Mraf(1)masseUNraf(1)
1065	do iq=1,nq
1066	Rrafk(1,iq)=Rraf(1)*Qraf(1,iq)
1067	Nrafk(1,iq)=Rrafk(1,iq)/masseU/dble(mmol(gc(iq)))
1068	enddo
1069	Zraf(nl)=z(nlx)
1070
1071	do l=2,nl-1
1072	Zraf(l)=Zraf(1)+(Zraf(nl)-Zraf(1))/dble(nl-1)*dble((l-1))
1073	indz=maxloc(z,mask=z <= Zraf(l))
1074	iz=indz(1)
1075	if (iz .lt. 1 .or. iz .gt. nlx) then
1076	print*,'danger',iz,nl,Zraf(l),l,Zraf(1),Zraf(nl),z,P,T,ig
1077	stop
1078	endif
1079	dZ=Zraf(l)-Zraf(l-1)
1080	! dZ=Zraf(l)-z(iz)
1081	facz=(Zraf(l)-z(iz))/(z(iz+1)-z(iz))
1082	Traf(l)=T(iz)+(T(iz+1)-T(iz))*facz
1083	do iq=1,nq
1084	! Qraf(l,iq)=qq(iz,iq)+(qq(iz+1,iq)-qq(iz,iq))*facz
1085	Rrafk(l,iq)=Rk(iz,iq)+(Rk(iz+1,iq)-Rk(iz,iq))*facZ
1086	Rrafk(l,iq)=Rk(iz,iq)(Rk(iz+1,iq)/Rk(iz,iq))*facZ
1087	enddo
1088	! Mraf(l)=1D0/(sum(qraf(l,:)/dble(mmol(gc(:)))))
1089	Rraf(l)=sum(Rrafk(l,:))
1090	do iq=1,nq
1091	Qraf(l,iq)=Rrafk(l,iq)/Rraf(l)
1092	enddo
1093	Mraf(l)=1D0/(sum(qraf(l,:)/dble(mmol(gc(:)))))
1094	! H=Konst*Traf(l)/Mraf(l)/g
1095	! H=Konst*T(iz)/M(iz)/g
1096	! Praf(l)=P(iz)*exp(-dZ/H)
1097	! Praf(l)=Praf(l-1)*exp(-dZ/H)
1098	! print,'iz',l,iz,Praf(il-1)exp(-dZ/H),z(iz),z(iz+1),H
1099	Nraf(l)=Rraf(l)/Mraf(l)/masseU
1100	Praf(l)=Nraf(l)kboltTraf(l)
1101	! Rraf(l)=Nraf(l)Mraf(l)masseU
1102	do iq=1,nq
1103	! Rrafk(l,iq)=Rraf(l)*Qraf(l,iq)
1104	Nrafk(l,iq)=Rrafk(l,iq)/dble(mmol(gc(iq)))/masseU
1105	if (Nrafk(l,iq) .lt. 0. .or. &
1106	& Nrafk(l,iq) .ne. Nrafk(l,iq)) then
1107	print*,'pb interpolation',l,iq,Nrafk(l,iq),Rrafk(l,iq), &
1108	& Qraf(l,iq),Rk(iz,iq),Rk(iz+1,iq),facZ,Zraf(l),z(iz)
1109	stop
1110	endif
1111	enddo
1112	enddo
1113	Zraf(nl)=z(nlx)
1114	Traf(nl)=T(nlx)
1115	do iq=1,nq
1116	Rrafk(nl,iq)=Rk(nlx,iq)
1117	Qraf(nl,iq)=Rk(nlx,iq)/R(nlx)
1118	Nrafk(nl,iq)=Rk(nlx,iq)/dble(mmol(gc(iq)))/masseU
1119	enddo
1120	Nraf(nl)=sum(Nrafk(nl,:))
1121	Praf(nl)=Nraf(nl)kboltTraf(nl)
1122	Mraf(nl)=1D0/sum(Qraf(nl,:)/dble(mmol(gc(:))))
1123	END
1124
1125	SUBROUTINE CORRMASS(qq,qint,FacMass,nl,nq)
1126	IMPLICIT NONE
1127	INTEGER :: nl,nq,l,nn
1128	REAL*8,DIMENSION(nl,nq) :: qq,qint,FacMass
1129
1130	do nn=1,nq
1131	do l=1,nl
1132	FacMass(l,nn)=qq(l,nn)/qint(l,nn)
1133	enddo
1134	enddo
1135
1136	END
1137
1138
1139	SUBROUTINE DCOEFF(nn,Dij,P,T,N,Nk,D,nl,nq)
1140	IMPLICIT NONE
1141	REAL*8,DIMENSION(nl) :: N,T,D,P
1142	REAL*8,DIMENSION(nl,nq) :: Nk
1143	INTEGER :: nn,nl,nq,l,iq
1144	REAL,DIMENSION(nq,nq) :: Dij
1145	REAL*8 :: interm,P0,T0,ptfac,dfac
1146
1147	P0=1D5
1148	T0=273d0
1149
1150
1151	do l=1,nl
1152	ptfac=(P0/P(l))(T(l)/T0)*1.75d0
1153	D(l)=0d0
1154	interm=0d0
1155	do iq=1,nq
1156	if (iq .ne. nn) then
1157	dfac=dble(dij(nn,iq))*ptfac
1158	interm=interm+Nk(l,iq)/N(l)/dfac
1159	endif
1160	enddo
1161	D(l)=1d0/interm
1162	enddo
1163	END
1164
1165	SUBROUTINE HSCALEREAL(nn,Nk,Dz,H,nl,nq)
1166	IMPLICIT NONE
1167	INTEGER :: nn,nl,nq,l
1168	REAL*8,DIMENSION(nl) :: H
1169	REAL*8,DIMENSION(nl,nq) :: Nk
1170	REAL*8 :: Dz
1171
1172	H(1)=(-3D0Nk(1,nn)+4d0NK(2,nn)-Nk(3,nn))/(2D0*DZ)/ &
1173	& NK(1,nn)
1174
1175	H(1)=-1D0/H(1)
1176
1177	DO l=2,nl-1
1178	H(l)=(Nk(l+1,nn)-NK(l-1,nn))/(2D0*DZ)/NK(l,nn)
1179	H(l)=-1D0/H(l)
1180	ENDDO
1181
1182	H(nl)=(3D0Nk(nl,nn)-4D0Nk(nl-1,nn)+Nk(nl-2,nn))/(2D0*DZ)/ &
1183	& Nk(nl,nn)
1184	H(nl)=-1D0/H(nl)
1185
1186	do l=1,nl
1187	if (abs(H(l)) .lt. 100.) then
1188	! print*,'H',l,H(l),Nk(l,nn),nn
1189	endif
1190	enddo
1191
1192	END
1193
1194	SUBROUTINE VELVERT(nn,T,H,D,Dz,masse,W,nl)
1195	IMPLICIT NONE
1196	INTEGER :: l,nl,nn
1197	REAL*8,DIMENSION(nl) :: T,H,D,W,DT
1198	REAL*8 :: Dz,Hmol,masse
1199	REAL*8 :: kbolt,masseU,Konst,g
1200	masseU=1.660538782d-27
1201	kbolt=1.3806504d-23
1202	Konst=Kbolt/masseU
1203	g=3.72d0
1204
1205	DT(1)=1D0/T(1)(-3D0T(1)+4D0T(2)-T(3))/(2D0DZ)
1206	Hmol=Konst*T(1)/masse/g
1207	W(1)=-D(1)*(1D0/H(1)-1D0/Hmol-DT(1))
1208
1209	DO l=2,nl-1
1210	DT(l)=1D0/T(l)(T(l+1)-T(l-1))/(2D0DZ)
1211	Hmol=Konst*T(l)/masse/g
1212	W(l)=-D(l)*(1D0/H(l)-1D0/Hmol-DT(l))
1213	ENDDO
1214
1215	DT(nl)=1D0/T(nl)(3d0T(nl)-4D0T(nl-1)+T(nl-2))/(2D0DZ)
1216	Hmol=Konst*T(nl)/masse/g
1217	W(nl)=-D(nl)*(1D0/H(nl)-1D0/Hmol-DT(nl))
1218
1219	! do l=1,nl
1220	! print*,'W',W(l),D(l),H(l),DT(l)
1221	! enddo
1222
1223	END
1224
1225	SUBROUTINE TIMEDIFF(nn,H,W,TIME,nl)
1226	IMPLICIT NONE
1227	INTEGER :: nn,nl,l
1228	REAL*8,DIMENSION(nl) :: W,H,TIME
1229
1230	DO l=1,nl
1231	TIME(l)=abs(H(l)/W(l))
1232	if (TIME(l) .lt. 1.D-4) then
1233	! print*,'low tdiff',H(l),W(l),nn,l
1234	endif
1235	ENDDO
1236
1237	END
1238
1239
1240	SUBROUTINE DIFFPARAM(T,D,dz,alpha,beta,gama,delta,eps,nl)
1241	IMPLICIT NONE
1242	INTEGER :: nl,l
1243	REAL*8,DIMENSION(nl) :: T,D
1244	REAL*8 :: DZ
1245	REAL*8,DIMENSION(nl) :: alpha,beta,gama,delta,eps
1246
1247	! Compute alpha,beta and delta
1248	! lower altitude values
1249
1250	alpha(1)=1D0/T(1)(-3D0T(1)+4D0T(2)-T(3))/(2D0DZ)
1251	delta(1)=(-3D0D(1)+4D0D(2)-D(3))/(2D0*DZ)
1252	beta(1)=1D0/T(1)
1253
1254	do l=2,nl-1
1255	alpha(l)=1D0/T(l)(T(l+1)-T(l-1))/(2D0DZ)
1256	delta(l)=(D(l+1)-D(l-1))/(2D0*DZ)
1257	beta(l)=1D0/T(l)
1258	end do
1259
1260	! Upper altitude values
1261
1262	alpha(nl)=1D0/T(nl)(3D0T(nl)-4D0T(nl-1)+T(nl-2))/(2D0DZ)
1263	delta(nl)=(3D0D(nl)-4D0D(nl-1)+D(nl-2))/(2D0*DZ)
1264	beta(nl)=1D0/T(nl)
1265
1266	! Compute the gama and eps coefficients
1267	! Lower altitude values
1268
1269	gama(1)=(-3D0beta(1)+4D0beta(2)-beta(3))/(2d0*dz)
1270	eps(1)=(-3D0alpha(1)+4D0alpha(2)-alpha(3))/(2d0*dz)
1271
1272	do l=2,nl-1
1273	gama(l)=(beta(l+1)-beta(l-1))/(2d0*Dz)
1274	eps(l)=(alpha(l+1)-alpha(l-1))/(2d0*Dz)
1275	end do
1276
1277	gama(nl)=(3D0beta(nl)-4D0beta(nl-1)+beta(nl-2))/(2D0*DZ)
1278	eps(nl)=(3D0alpha(nl)-4D0alpha(nl-1)+alpha(nl-2))/(2D0*DZ)
1279
1280
1281	END
1282
1283
1284	SUBROUTINE MATCOEFF(alpha,beta,gama,delta,eps,Dad,rhoad, &
1285	& FacEsc,dz,dt,A,B,C,D,nl)
1286	IMPLICIT NONE
1287	INTEGER :: nl,l
1288	REAL*8, DIMENSION(nl) :: alpha,beta,gama,delta,eps,Dad,RHoad
1289	REAL*8 :: dz,dt,del1,del2,del3,FacEsc
1290	REAL*8, DIMENSION(nl) :: A,B,C,D
1291	del1=dt/2d0/dz
1292	del2=dt/dz/dz
1293	del3=dt
1294
1295	! lower boundary coefficients no change
1296	A(1)=0d0
1297	B(1)=1d0
1298	C(1)=0d0
1299	D(1)=rhoAd(1)
1300
1301	do l=2,nl-1
1302	A(l)=(delta(l)+(alpha(l)+beta(l))Dad(l))del1-Dad(l)*del2
1303	B(l)=-(delta(l)(alpha(l)+beta(l))+Dad(l)(gama(l)+eps(l)))*del3
1304	B(l)=B(l)+1d0+2d0Dad(l)del2
1305	C(l)=-(delta(l)+(alpha(l)+beta(l))Dad(l))del1-Dad(l)*del2
1306	D(l)=rhoAD(l)
1307	enddo
1308
1309
1310	! Upper boundary coefficients Diffusion profile
1311	C(nl)=0d0
1312	B(nl)=-1d0
1313	A(nl)=exp(-dZ(alpha(nl)+beta(nl)))FacEsc
1314	D(nl)=0D0
1315
1316
1317	END
1318
1319	SUBROUTINE Checkmass(X,Y,nl,nn)
1320	IMPLICIT NONE
1321
1322	INTEGER :: nl,nn
1323	REAL*8,DIMENSION(nl) :: X,Y
1324	REAL*8 Xtot,Ytot
1325
1326	Xtot=sum(X)
1327	Ytot=sum(Y)
1328
1329	if (abs((Xtot-Ytot)/Xtot) .gt. 1d-4) then
1330	print*,'no conservation for mass',Xtot,Ytot,nn
1331	endif
1332	END
1333
1334	SUBROUTINE Checkmass2(qold,qnew,P,il,nl,nn,nq)
1335	IMPLICIT NONE
1336	INTEGER :: nl,nn,l,nq,il
1337	REAL,DIMENSION(nl+1) :: P
1338	REAL*8,DIMENSION(nl,nq) :: qold,qnew
1339	REAL*8 :: DM,Mold,Mnew,g
1340	g=3.72d0
1341	DM=0d0
1342	Mold=0d0
1343	Mnew=0d0
1344	DO l=il,nl
1345	DM=DM+(qnew(l,nn)-qold(l,nn))*(dble(P(l))-dble(P(l+1)))/g
1346	Mold=Mold+qold(l,nn)*(dble(P(l))-dble(P(l+1)))/g
1347	Mnew=Mnew+qnew(l,nn)*(dble(P(l))-dble(P(l+1)))/g
1348	! print*,'l',l,qold(l,nn),qnew(l,nn),Mold,Mnew,DM,P(l),P(l+1)
1349	ENDDO
1350	IF (abs(DM/Mold) .gt. 1d-5) THEN
1351	Print*,'We dont conserve mas',nn,DM,Mold,Mnew,DM/Mold
1352	ENDIF
1353
1354	END
1355
1356	SUBROUTINE GCMGRID_P(Z,P,Q,T,Nk,Rk,qq,qnew,tt,tnew, &
1357	& pp,M,gc,nl,nq,nlx)
1358	IMPLICIT NONE
1359	#include "dimensions.h"
1360	INTEGER :: nl,nq,nlx,il,nn,iP
1361	INTEGER,DIMENSION(1) :: indP
1362	INTEGER,DIMENSION(nq) :: gc
1363	REAL*8,DIMENSION(nl) :: Z,P,T
1364	REAL*8,DIMENSION(nl,nq) :: Q,Nk,Rk
1365	REAL,DIMENSION(nqmx) :: M
1366	REAL*8,DIMENSION(nq) :: nNew
1367	REAL*8,DIMENSION(nlx) :: pp,tt,tnew
1368	REAL*8,DIMENSION(nlx) :: rhonew
1369	REAL*8,DIMENSION(nlx,nq) :: qq,qnew,rhoknew
1370	REAL*8 :: kbolt,masseU,Konst,g,Dz,facP,Hi
1371	REAL*8 :: Znew,Znew2,Pnew,Pnew2
1372	masseU=1.660538782d-27
1373	kbolt=1.3806504d-23
1374	Konst=Kbolt/masseU
1375	g=3.72d0
1376	Dz=Z(2)-Z(1)
1377	Znew=Z(nl)
1378	Znew2=Znew+dz
1379	! print*,'dz',Znew,Znew2,dz
1380	nNew(1:nq)=Nk(nl,1:nq)
1381	Pnew=P(nl)
1382
1383	do il=1,nlx
1384	! print*,'il',il,pp(il),P(1),P(nl)
1385	if (pp(il) .ge. P(1)) then
1386	qnew(il,:)=qq(il,:)
1387	tnew(il)=tt(il)
1388	endif
1389	if (pp(il) .lt. P(1)) then
1390	if (pp(il) .gt. P(nl)) then
1391	indP=maxloc(P,mask=P < pp(il))
1392	iP=indP(1)-1
1393	if (iP .lt. 1 .or. iP .gt. nl) then
1394	print*,'danger 2',iP,nl,pp(il)
1395	endif
1396	facP=(pp(il)-P(ip))/(P(ip+1)-P(ip))
1397	! print*,'P',P(ip),P(ip+1),facP,indP,iP
1398
1399	! do nn=1,nq
1400	! qnew(il,nn)=Q(iP,nn)+
1401	! & (Q(ip+1,nn)-Q(ip,nn))*facP
1402	! enddo
1403
1404	do nn=1,nq
1405	rhoknew(il,nn)=Rk(iP,nn)+ &
1406	& (Rk(ip+1,nn)-Rk(ip,nn))*facP
1407	enddo
1408	tnew(il)=T(iP)+(T(iP+1)-T(iP))*facP
1409	rhonew(il)=sum(rhoknew(il,:))
1410	do nn=1,nq
1411	qnew(il,nn)=rhoknew(il,nn)/rhonew(il)
1412	enddo
1413
1414	else ! pp < P(nl) need to extrapolate density of each specie
1415	Pnew2=Pnew
1416
1417	do while (pnew2 .ge. pp(il))
1418	do nn=1,nq
1419	Hi=Konst*T(nl)/dble(M(gc(nn)))/g
1420	Nnew(nn)=Nnew(nn)*exp(-dZ/Hi)
1421	enddo
1422	Pnew=Pnew2
1423	Pnew2=kboltT(nl)sum(Nnew(:))
1424	Znew=Znew2
1425	Znew2=Znew2+Dz
1426	! print*,'test',Pnew2,Znew2,Nnew(nq),pp(il)
1427	enddo
1428
1429	facP=(pp(il)-Pnew)/(Pnew2-Pnew)
1430
1431	! do nn=1,nq
1432	! qnew(il,nn)=dble(M(gc(nn)))*Nnew(nn)
1433	! & /sum(dble(M(gc(:)))*Nnew(:))
1434	! enddo
1435
1436	do nn=1,nq
1437	rhoknew(il,nn)=dble(M(gc(nn)))*Nnew(nn)
1438	enddo
1439	rhonew(il)=sum(rhoknew(il,:))
1440	do nn=1,nq
1441	qnew(il,nn)=rhoknew(il,nn)/rhonew(il)
1442	enddo
1443	tnew(il)=T(nl)
1444	endif
1445	endif
1446	enddo
1447
1448	END
1449
1450	SUBROUTINE GCMGRID_P2(Z,P,Q,T,Nk,Rk,qq,qnew,tt,tnew &
1451	& ,pp,M,gc,nl,nq,nlx,facM)
1452	IMPLICIT NONE
1453	#include "dimensions.h"
1454	INTEGER :: nl,nq,nlx,il,nn,iP
1455	INTEGER,DIMENSION(1) :: indP
1456	INTEGER,DIMENSION(nq) :: gc
1457	REAL*8,DIMENSION(nl) :: Z,P,T
1458	REAL*8,DIMENSION(nl,nq) :: Q,Nk,Rk
1459	REAL,DIMENSION(nqmx) :: M
1460	REAL*8,DIMENSION(nq) :: nNew
1461	REAL*8,DIMENSION(nlx) :: pp,rhonew,tt,tnew
1462	REAL*8,DIMENSION(nlx,nq) :: qq,qnew,facM,rhoknew
1463	REAL*8 :: kbolt,masseU,Konst,g,Dz,facP,Hi
1464	REAL*8 :: Znew,Znew2,Pnew,Pnew2
1465	masseU=1.660538782d-27
1466	kbolt=1.3806504d-23
1467	Konst=Kbolt/masseU
1468	g=3.72d0
1469	Dz=Z(2)-Z(1)
1470	Znew=Z(nl)
1471	Znew2=Znew+dz
1472	! print*,'dz',Znew,Znew2,dz
1473	nNew(1:nq)=Nk(nl,1:nq)
1474	Pnew=P(nl)
1475
1476	do il=1,nlx
1477	! print*,'il',il,pp(il),P(1),P(nl)
1478	if (pp(il) .ge. P(1)) then
1479	qnew(il,:)=qq(il,:)
1480	tnew(il)=tt(il)
1481	endif
1482	if (pp(il) .lt. P(1)) then
1483	if (pp(il) .gt. P(nl)) then
1484	indP=maxloc(P,mask=P < pp(il))
1485	iP=indP(1)-1
1486	if (iP .lt. 1 .or. iP .gt. nl) then
1487	print*,'danger 3',iP,nl,pp(il)
1488	endif
1489	facP=(pp(il)-P(ip))/(P(ip+1)-P(ip))
1490	! print*,'P',P(ip),P(ip+1),facP,indP,iP
1491
1492	! do nn=1,nq
1493	! qnew(il,nn)=Q(iP,nn)+
1494	! & (Q(ip+1,nn)-Q(ip,nn))*facP
1495	! enddo
1496
1497	do nn=1,nq
1498	rhoknew(il,nn)=(RK(iP,nn)+ &
1499	& (RK(iP+1,nn)-Rk(iP,nn))facP)facM(il,nn)
1500	enddo
1501	tnew(il)=T(iP)+(T(ip+1)-T(iP))*facP
1502	rhonew(il)=sum(rhoknew(il,:))
1503	do nn=1,nq
1504	qnew(il,nn)=rhoknew(il,nn)/rhonew(il)
1505	enddo
1506
1507	else ! pp < P(nl) need to extrapolate density of each specie
1508	Pnew2=Pnew
1509
1510	do while (pnew2 .ge. pp(il))
1511	do nn=1,nq
1512	Hi=Konst*T(nl)/dble(M(gc(nn)))/g
1513	Nnew(nn)=Nnew(nn)*exp(-dZ/Hi)
1514	enddo
1515	Pnew=Pnew2
1516	Pnew2=kboltT(nl)sum(Nnew(:))
1517	Znew=Znew2
1518	Znew2=Znew2+Dz
1519	! print*,'test',Pnew2,Znew2,Nnew(nq),pp(il)
1520	enddo
1521
1522	facP=(pp(il)-Pnew)/(Pnew2-Pnew)
1523
1524	! do nn=1,nq
1525	! qnew(il,nn)=dble(M(gc(nn)))*Nnew(nn)
1526	! & /sum(dble(M(gc(:)))*Nnew(:))
1527	! enddo
1528
1529	do nn=1,nq
1530	rhoknew(il,nn)=dble(M(gc(nn)))Nnew(nn)FacM(il,nn)
1531	enddo
1532	rhonew(il)=sum(rhoknew(il,:))
1533	do nn=1,nq
1534	qnew(il,nn)=rhoknew(il,nn)/rhonew(il)
1535	enddo
1536	tnew(il)=T(nl)
1537
1538	endif
1539	endif
1540	enddo
1541
1542	END

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

Download in other formats:

Original Format