Context Navigation

moldiff_red.F90 @ 1198

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

Mars GCM:

IMPORTANT CHANGE: Removed all reference/use of ngridmx (dimphys.h) in routines (necessary prerequisite to using parallel dynamics); in most cases this just means adding 'ngrid' as routine argument, and making local saved variables allocatable (and allocated at first call). In the process, had to convert many *.h files to equivalent modules: yomaer.h => yomaer_h.F90 , surfdat.h => surfdat_h.F90 , comsaison.h => comsaison_h.F90 , yomlw.h => yomlw_h.F90 , comdiurn.h => comdiurn_h.F90 , dimradmars.h => dimradmars_mod.F90 , comgeomfi.h => comgeomfi_h.F90, comsoil.h => comsoil_h.F90 , slope.h => slope_mod.F90
Also updated EOF routines, everything is now in eofdump_mod.F90
Removed unused routine lectfux.F (in dyn3d)

File size: 38.6 KB

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

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

Download in other formats:

Original Format