Context Navigation

moldiff_red.F90 @ 1036

Last change on this file since 1036 was 1036, checked in by emillour, 11 years ago

Mars GCM: (a first step towards using parallel dynamics)

IMPORTANT CHANGE: Implemented dynamic tracers. It is no longer necessary to compile the model with the '-t #' option, number of tracers is simply read from tracer.def file (as before). Adapted makegcm_* scripts (and co.) accordingly. Technical aspects of the switch to dynamic tracers are:
- advtrac.h (in dyn3d) removed and replaced by module infotrac.F
- tracer.h (in phymars) removed and replaced by module tracer_mod.F90 (which contains nqmx, the number of tracers, etc. and can be used anywhere in the physics).
Included some side cleanups: removed unused files (in dyn3d) anldoppler2.F, anl_mcdstats.F and anl_stats-diag.F, and all the unecessary dimensions.* files in grid/dimension.
Checked that changes are clean and that GCM yields identical results (in debug mode) to previous svn version.

File size: 38.5 KB

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

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

Download in other formats:

Original Format