Context Navigation

mztvc.F @ 455

Last change on this file since 455 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

Rev	Line
[414]	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