Context Navigation

mztvc.F @ 469

Last change on this file since 469 was 414, checked in by aslmd, 14 years ago

LMDZ.MARS : NEW NLTE MODEL FROM GRANADA AMIGOS

23/11/11 == FGG + MALV

New parameterization of the NLTE 15 micron cooling. The old parameterization is kept as an option, including or not variable atomic oxygen concentration. A new flag is introduced in callphys.def, nltemodel, to select which parameterization wants to be used (new one, old one with variable [O], or old one with fixed [O], see below). Includes many new subroutines and commons in phymars. Some existing routines are also modified:

-physiq.F. Call to the new subroutine NLTE_leedat in first call. Call to nltecool modified to allow for variable atomic oxygen. Depending on the value of nltemodel, the new subroutine NLTEdlvr09_TCOOL is called instead of nltecool.

-inifis.F. Reading of nltemodel is added.

-callkeys.h Declaration of nltemodel is added.

The following lines should be added to callphys.def (ideally after setting callnlte):

# NLTE 15um scheme to use.
# 0-> Old scheme, static oxygen
# 1-> Old scheme, dynamic oxygen
# 2-> New scheme
nltemodel = 2

A new directory, NLTEDAT, has to be included in datagcm.

Improvements into NLTE NIR heating parameterization to take into account variability with O/CO2 ratio and SZA. A new subroutine, NIR_leedat.F, and a new common, NIRdata.h, are included. A new flag, nircorr, is added in callphys.def, to include or not these corrections. The following files are modified:

-nirco2abs.F: nq and pq are added in the arguments. The corrections factors are interpolated to the GCM grid and included in the clculation. A new subroutine, interpnir, is added at the end of the file.

-physiq.F: Call to NIR_leedat added at first call. Modified call to nirco2abs

-inifis: Reading new flag nircorr.

-callkeys.h: Declaration of nircorr.

The following lines have to be added to callphys.def (ideally after callnirco2):

# NIR NLTE correction for variable SZA and O/CO2?
# matters only if callnirco2=T
# 0-> no correction
# 1-> include correction
nircorr=1

A new data file, NIRcorrection_feb2011.dat, has to be included in datagcm.

Small changes to the molecular diffusion scheme to fix the number of species considered, to avoid problems when compiling with more than 15 tracers (for example, when CH4 is included). Modified subroutines: aeronomars/moldiff.F and aeronomars/moldiffcoeff.F

File size: 14.4 KB

Line
1	c***********************************************************************
2
3	subroutine mztvc ( ig,vc, ib,isot,
4	@ iirw,iimu,itauout,icfout,itableout )
5
6	c jul 2011 malv+fgg
7	c***********************************************************************
8
9	implicit none
10
11	include 'comcstfi.h'
12	include 'nltedefs.h'
13	include 'nlte_atm.h'
14	include 'nlte_data.h'
15	include 'nlte_curtis.h'
16	include 'tcr_15um.h'
17	include 'nlte_results.h'
18
19	c arguments
20	integer ig ! ADDED FOR TRACEBACK
21	real*8 cf(nl,nl), cfup(nl,nl), cfdw(nl,nl) ! o
22	real*8 vc(nl), taugr(nl) ! o
23	integer ib ! i
24	integer isot ! i
25	integer iirw ! i
26	integer iimu ! i
27	integer itauout ! i
28	integer icfout ! i
29	integer itableout ! i
30
31	c local variables and constants
32	integer i, in, ir, im, k ,j
33	integer nmu
34	parameter (nmu = 8)
35	real*8 tau(nl,nl)
36	real*8 tauinf(nl)
37	real*8 con(nzy), coninf
38	real*8 c1, c2
39	real*8 t1, t2
40	real*8 p1, p2
41	real*8 mr1, mr2
42	real*8 st1, st2
43	real*8 c1box(70), c2box(70)
44	real*8 ff ! to avoid too small numbers
45	real*8 tvtbs(nzy)
46	real*8 st, beta, ts, eqwmu
47	real*8 mu(nmu), amu(nmu)
48	real*8 zld(nl), zyd(nzy)
49	real*8 correc
50	real deltanux ! width of vib-rot band (cm-1)
51	character isotcode*2
52	integer idummy
53	real*8 Desp,wsL
54
55	c formats
56	111 format(a1)
57	112 format(a2)
58	101 format(i1)
59	202 format(i2)
60	180 format(a80)
61	181 format(a80)
62	c***********************************************************************
63
64	c some needed values
65	! rl=sqrt(log(2.d0))
66	! pi2 = 3.14159265358989d0
67	beta = 1.8d0
68	! beta = 1.0d0
69	idummy = 0
70	Desp = 0.0d0
71	wsL = 0.0d0
72
73	!write (,) ' MZTUD/ iirw = ', iirw
74
75
76	c esto es para que las subroutines de mztfsub calculen we
77	c de la forma apropiada para mztf, no para fot
78	icls=icls_mztf
79
80	c codigos para filenames
81	! if (isot .eq. 1) isotcode = '26'
82	! if (isot .eq. 2) isotcode = '28'
83	! if (isot .eq. 3) isotcode = '36'
84	! if (isot .eq. 4) isotcode = '27'
85	! if (isot .eq. 5) isotcode = '62'
86	! if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5
87	! @ .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9) then
88	! write (ibcode1,101) ib
89	! else
90	! write (ibcode2,202) ib
91	! endif
92	! write (*,'( 30h calculating curtis matrix : ,2x,
93	! @ 8h band = ,i2,2x, 11h isotope = ,i2)') ib, isot
94
95	c integration in angle !!!!!!!!!!!!!!!!!!!!
96
97	c------- diffusivity approx.
98	if (iimu.eq.1) then
99	! write (,) ' diffusivity approx. beta = ',beta
100	mu(1) = 1.0d0
101	amu(1)= 1.0d0
102	c-------data for 8 points integration
103	elseif (iimu.eq.4) then
104	write (,)' 4 points for the gauss-legendre angle quadrature.'
105	mu(1)=(1.0d0+0.339981043584856)/2.0d0
106	mu(2)=(1.0d0-0.339981043584856)/2.0d0
107	mu(3)=(1.0d0+0.861136311594053)/2.0d0
108	mu(4)=(1.0d0-0.861136311594053)/2.0d0
109	amu(1)=0.652145154862546
110	amu(2)=amu(1)
111	amu(3)=0.347854845137454
112	amu(4)=amu(3)
113	beta=1.0d0
114	c-------data for 8 points integration
115	elseif(iimu.eq.8) then
116	write (,)' 8 points for the gauss-legendre angle quadrature.'
117	mu(1)=(1.0d0+0.183434642495650)/2.0d0
118	mu(2)=(1.0d0-0.183434642495650)/2.0d0
119	mu(3)=(1.0d0+0.525532409916329)/2.0d0
120	mu(4)=(1.0d0-0.525532409916329)/2.0d0
121	mu(5)=(1.0d0+0.796666477413627)/2.0d0
122	mu(6)=(1.0d0-0.796666477413627)/2.0d0
123	mu(7)=(1.0d0+0.960289856497536)/2.0d0
124	mu(8)=(1.0d0-0.960289856497536)/2.0d0
125	amu(1)=0.362683783378362
126	amu(2)=amu(1)
127	amu(3)=0.313706645877887
128	amu(4)=amu(3)
129	amu(5)=0.222381034453374
130	amu(6)=amu(5)
131	amu(7)=0.101228536290376
132	amu(8)=amu(7)
133	beta=1.0d0
134	end if
135	c!!!!!!!!!!!!!!!!!!!!!!!
136
137	ccc
138	ccc determine abundances included in the absorber amount
139	ccc
140
141	c first, set up the grid ready for interpolation.
142	do i=1,nzy
143	zyd(i) = dble(zy(i))
144	enddo
145	do i=1,nl
146	zld(i) = dble(zl(i))
147	enddo
148
149	c vibr. temp of the bending mode :
150	if (isot.eq.1) call interdp ( tvtbs,zyd,nzy, v626t1,zld,nl, 1 )
151	if (isot.eq.2) call interdp ( tvtbs,zyd,nzy, v628t1,zld,nl, 1 )
152	if (isot.eq.3) call interdp ( tvtbs,zyd,nzy, v636t1,zld,nl, 1 )
153	if (isot.eq.4) call interdp ( tvtbs,zyd,nzy, v627t1,zld,nl, 1 )
154	!if (isot.eq.5) call interdp ( tvtbs,zxd,nz, vcot1,zld,nl, 1 )
155
156	c 2nd: correccion a la n10(i) (cantidad de absorbente en el lower state)
157	c por similitud a la que se hace en cza.for ; esto solo se hace para CO2
158
159	!write (,) 'imr(isot) = ', isot, imr(isot)
160	do i=1,nzy
161	if (isot.eq.5) then
162	con(i) = dble( coy(i) * imrco )
163	else
164	con(i) = dble( co2y(i) * imr(isot) )
165	correc = 2.d0 * dexp( dble(-ee*elow(isot,2))/tvtbs(i) )
166	con(i) = con(i) * ( 1.d0 - correc )
167	! write (,) ' iz, correc, co2y(i), con(i) =',
168	! @ i,correc,co2y(i),con(i)
169	endif
170
171	!-----------------------------------------------------------------
172	! mlp & cristina. 17 july 1996 change the calculation of mr.
173	! it is used for calculating partial press
174	! alpha = alpha(self,co2)pp +alpha(n2)(pt-pp)
175	! for an isotope, if mr is obtained by
176	! co2*imr(iso)/nt
177	! we are considerin collisions with other co2 isotopes
178	! (including the major one, 626) as if they were with n2.
179	! assuming mr as co2/nt, we consider collisions
180	! of type 628-626 as of 626-626 instead of as 626-n2.
181	! mrx(i)=con(i)/ntx(i) ! old malv
182	! mrx(i)= dble(co2x(i)/ntx(i)) ! mlp & crs
183
184	! jan 98:
185	! esta modif de mlp implica anular el correc (deberia revisar esto)
186
187	mr(i) = dble(co2y(i)/nty(i)) ! malv, jan 98
188
189	!-----------------------------------------------------------------
190
191	end do
192
193	! como beta y 1.d5 son comunes a todas las weighted absorber amounts,
194	! los simplificamos:
195	! coninf = beta * 1.d5 * dble( con(n) / log( con(n-1) / con(n) ) )
196	!write (,) ' con(nz), con(nz-1) =', con(nz), con(nz-1)
197	coninf = dble( con(nzy) / log( con(nzy-1) / con(nzy) ) )
198	!write (,) ' coninf =', coninf
199
200	ccc
201	ccc temp dependence of the band strength and
202	ccc nlte correction factor for the absorber amount
203	ccc
204	call mztf_correccion ( coninf, con, ib, isot, itableout )
205
206	ccc
207	ccc reads histogrammed spectral data (strength for lte and vmr=1)
208	ccc
209	!hfile1 = dirspec//'hi'//dn !Ya no hacemos distincion d/n en esto
210	!! hfile1 = dirspec//'hid' !(see why in his.for)
211	! hfile1='hid'
212	!! if (ib.eq.13 .or. ib.eq.14 ) hfile1 = dirspec//'his'
213	! if (ib.eq.13 .or. ib.eq.14 ) hfile1 = 'his'
214
215	! if(ib.eq.1.or.ib.eq.2.or.ib.eq.3.or.ib.eq.4.or.ib.eq.5
216	! @ .or.ib.eq.6.or.ib.eq.7.or.ib.eq.8.or.ib.eq.9) then
217	! if (isot.eq.1) hisfile = hfile1//'26-'//ibcode1//'.dat'
218	! if (isot.eq.2) hisfile = hfile1//'28-'//ibcode1//'.dat'
219	! if (isot.eq.3) hisfile = hfile1//'36-'//ibcode1//'.dat'
220	! if (isot.eq.4) hisfile = hfile1//'27-'//ibcode1//'.dat'
221	! if (isot.eq.5) hisfile = hfile1//'62-'//ibcode1//'.dat'
222	! else
223	! if (isot.eq.1) hisfile = hfile1//'26-'//ibcode2//'.dat'
224	! if (isot.eq.2) hisfile = hfile1//'28-'//ibcode2//'.dat'
225	! if (isot.eq.3) hisfile = hfile1//'36-'//ibcode2//'.dat'
226	! if (isot.eq.4) hisfile = hfile1//'27-'//ibcode2//'.dat'
227	! if (isot.eq.5) hisfile = hfile1//'62-'//ibcode2//'.dat'
228	! endif
229	! write (,) 'hisfile: ', hisfile
230
231	! the argument to rhist is to make this compatible with mztf_comp.f,
232	! which is a useful modification of mztf.f (to change strengths of bands
233	! call rhist (1.0)
234	if(ib.eq.1) then
235	if(isot.eq.1) then !Case 1
236	mm=mm_c1
237	nbox=nbox_c1
238	tmin=tmin_c1
239	tmax=tmax_c1
240	do i=1,nbox_max
241	no(i)=no_c1(i)
242	dist(i)=dist_c1(i)
243	do j=1,nhist
244	sk1(j,i)=sk1_c1(j,i)
245	xls1(j,i)=xls1_c1(j,i)
246	xln1(j,i)=xln1_c1(j,i)
247	xld1(j,i)=xld1_c1(j,i)
248	enddo
249	enddo
250	do j=1,nhist
251	thist(j)=thist_c1(j)
252	enddo
253	else if(isot.eq.2) then !Case 2
254	mm=mm_c2
255	nbox=nbox_c2
256	tmin=tmin_c2
257	tmax=tmax_c2
258	do i=1,nbox_max
259	no(i)=no_c2(i)
260	dist(i)=dist_c2(i)
261	do j=1,nhist
262	sk1(j,i)=sk1_c2(j,i)
263	xls1(j,i)=xls1_c2(j,i)
264	xln1(j,i)=xln1_c2(j,i)
265	xld1(j,i)=xld1_c2(j,i)
266	enddo
267	enddo
268	do j=1,nhist
269	thist(j)=thist_c2(j)
270	enddo
271	else if(isot.eq.3) then !Case 3
272	mm=mm_c3
273	nbox=nbox_c3
274	tmin=tmin_c3
275	tmax=tmax_c3
276	do i=1,nbox_max
277	no(i)=no_c3(i)
278	dist(i)=dist_c3(i)
279	do j=1,nhist
280	sk1(j,i)=sk1_c3(j,i)
281	xls1(j,i)=xls1_c3(j,i)
282	xln1(j,i)=xln1_c3(j,i)
283	xld1(j,i)=xld1_c3(j,i)
284	enddo
285	enddo
286	do j=1,nhist
287	thist(j)=thist_c3(j)
288	enddo
289	else if(isot.eq.4) then !Case 4
290	mm=mm_c4
291	nbox=nbox_c4
292	tmin=tmin_c4
293	tmax=tmax_c4
294	do i=1,nbox_max
295	no(i)=no_c4(i)
296	dist(i)=dist_c4(i)
297	do j=1,nhist
298	sk1(j,i)=sk1_c4(j,i)
299	xls1(j,i)=xls1_c4(j,i)
300	xln1(j,i)=xln1_c4(j,i)
301	xld1(j,i)=xld1_c4(j,i)
302	enddo
303	enddo
304	do j=1,nhist
305	thist(j)=thist_c4(j)
306	enddo
307	else
308	write(,)'isot must be 2,3 or 4 for ib=1!!'
309	write(,)'stop at mztvc/310'
310	stop
311	endif
312	else if (ib.eq.2) then
313	if(isot.eq.1) then !Case 5
314	mm=mm_c5
315	nbox=nbox_c5
316	tmin=tmin_c5
317	tmax=tmax_c5
318	do i=1,nbox_max
319	no(i)=no_c5(i)
320	dist(i)=dist_c5(i)
321	do j=1,nhist
322	sk1(j,i)=sk1_c5(j,i)
323	xls1(j,i)=xls1_c5(j,i)
324	xln1(j,i)=xln1_c5(j,i)
325	xld1(j,i)=xld1_c5(j,i)
326	enddo
327	enddo
328	do j=1,nhist
329	thist(j)=thist_c5(j)
330	enddo
331	else
332	write(,)'isot must be 1 for ib=2!!'
333	write(,)'stop at mztvc/334'
334	stop
335	endif
336	else if (ib.eq.3) then
337	if(isot.eq.1) then !Case 6
338	mm=mm_c6
339	nbox=nbox_c6
340	tmin=tmin_c6
341	tmax=tmax_c6
342	do i=1,nbox_max
343	no(i)=no_c6(i)
344	dist(i)=dist_c6(i)
345	do j=1,nhist
346	sk1(j,i)=sk1_c6(j,i)
347	xls1(j,i)=xls1_c6(j,i)
348	xln1(j,i)=xln1_c6(j,i)
349	xld1(j,i)=xld1_c6(j,i)
350	enddo
351	enddo
352	do j=1,nhist
353	thist(j)=thist_c6(j)
354	enddo
355	else
356	write(,)'isot must be 1 for ib=3!!'
357	write(,)'stop at mztvc/358'
358	stop
359	endif
360	else if (ib.eq.4) then
361	if(isot.eq.1) then !Case 7
362	mm=mm_c7
363	nbox=nbox_c7
364	tmin=tmin_c7
365	tmax=tmax_c7
366	do i=1,nbox_max
367	no(i)=no_c7(i)
368	dist(i)=dist_c7(i)
369	do j=1,nhist
370	sk1(j,i)=sk1_c7(j,i)
371	xls1(j,i)=xls1_c7(j,i)
372	xln1(j,i)=xln1_c7(j,i)
373	xld1(j,i)=xld1_c7(j,i)
374	enddo
375	enddo
376	do j=1,nhist
377	thist(j)=thist_c7(j)
378	enddo
379	else
380	write(,)'isot must be 1 for ib=4!!'
381	write(,)'stop at mztvc/382'
382	stop
383	endif
384	else
385	write(,)'ib must be 1,2,3 or 4!!'
386	write(,)'stop at mztvc/387'
387	endif
388
389
390	c******
391	c****** calculation of tau(1,ir) for 1<=r
392	c******
393	call initial
394
395	ff=1.0e10
396
397	in=1
398
399	tau(in,1) = 1.d0
400
401	call initial
402	call intz (zl(in), c1,p1,mr1,t1, con)
403	do kr=1,nbox
404	ta(kr) = t1
405	end do
406	call interstrength (st1,t1,ka,ta)
407	do kr=1,nbox
408	c1box(kr) = c1 * ka(kr) * dble(deltaz)
409	end do
410	c1 = c1 * st1 * dble(deltaz)
411
412	do 2 ir=2,nl
413
414	call intz (zl(ir), c2,p2,mr2,t2, con)
415	do kr=1,nbox
416	ta(kr) = t2
417	end do
418	call interstrength (st2,t2,ka,ta)
419	do kr=1,nbox
420	c2box(kr) = c2 * ka(kr) * dble(deltaz)
421	end do
422	c2 = c2 * st2 * dble(deltaz)
423
424	aa = aa + ( p1mr1(c1ff) + p2mr2(c2ff)) / 2.d0
425	bb = bb + ( p1c1 + p2c2 ) / 2.d0
426	cc = cc + ( c1 + c2 ) / 2.d0
427	dd = dd + ( t1c1 + t2c2 ) / 2.d0
428	do kr=1,nbox
429	ccbox(kr) = ccbox(kr) + ( c1box(kr) + c2box(kr) ) /2.d0
430	ddbox(kr) = ddbox(kr) + ( t1c1box(kr) + t2c2box(kr) ) /2.d0
431	end do
432
433	mr1=mr2
434	t1=t2
435	c1=c2
436	p1=p2
437	do kr=1,nbox
438	c1box(kr) = c2box(kr)
439	end do
440
441	pt = bb / cc
442	pp = aa / (cc * ff)
443
444	ts = dd/cc
445	do kr=1,nbox
446	ta(kr) = ddbox(kr) / ccbox(kr)
447	end do
448	call interstrength(st,ts,ka,ta)
449	call intershape(alsa,alna,alda,ta)
450
451
452	eqwmu = 0.0d0
453	do im = 1,iimu
454	eqw=0.0d0
455	do kr=1,nbox
456	ua(kr) = ccbox(kr) / ka(kr) * beta * 1.0d5 / mu(im)
457	call findw(ig,iirw, idummy,c1,p1,Desp,wsL)
458	if ( i_supersat .eq. 0 ) then
459	eqw=eqw+no(kr)*w
460	else
461	eqw = w + (no(kr)-1) * ( asat_box*dist(kr) )
462	endif
463	end do
464	eqwmu = eqwmu + eqw * mu(im)*amu(im)
465	end do
466
467	tau(in,ir) = exp( - eqwmu / dble(deltanu(isot,ib)) )
468
469	2 continue
470
471
472
473	c
474	c due to the simmetry of the transmittances
475	c
476	do in=nl,2,-1
477	tau(in,1) = tau(1,in)
478	end do
479
480	vc(1) = 0.0d0
481	vc(nl) = 0.0d0
482	do in=2,nl-1 ! poner aqui nl-1 luego
483	vc(in) = pideltanu(isot,ib)/( 2.d0deltaz1.d5 )
484	@ ( tau(in-1,1) - tau(in+1,1) )
485	end do
486
487
488	c end
489	return
490	end

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

Download in other formats:

Original Format