Context Navigation

mztfsub_overlap.F @ 463

Last change on this file since 463 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: 44.5 KB

Rev	Line
[414]	1	c***********************************************************************
	2	c File with all subroutines required by mztf
	3	c
	4	c jan 98 malv basado en mztfsub_solar
	5	c jul 2011 malv+fgg adapted to LMD-MGCM
	6	c
	7	c contiene:
	8	c initial
	9	c intershape
	10	c interstrength
	11	c intz
	12	c rhist
	13	c interh
	14	c we
	15	c simrul
	16	c fi
	17	c f
	18	c findw
	19	c voigtf
	20	c***********************************************************************
	21
	22	c ****************************************************************
	23	subroutine initial
	24
	25	c ma & crs !evita troubles 16-july-96
	26	c ****************************************************************
	27
	28	implicit none
	29
	30	include 'nltedefs.h'
	31	include 'nlte_curtis.h'
	32
	33	c local variables
	34	integer i
	35
	36	c ***************
	37
	38	eqw = 0.0d00
	39	aa = 0.0d00
	40	bb = 0.0d00
	41	cc = 0.0d00
	42	dd = 0.0d00
	43
	44	do i=1,nbox
	45	ua(i) = 0.0d0
	46	ccbox(i) = 0.0d0
	47	ddbox(i) = 0.0d0
	48	end do
	49
	50	return
	51	end
	52
	53	c **********************************************************************
	54	subroutine intershape(alsx,alnx,adx,xtemp)
	55	c interpolates the line shape parameters at a temperature xtemp from
	56	c input histogram data.
	57	c **********************************************************************
	58
	59	implicit none
	60
	61	include 'nltedefs.h'
	62	include 'nlte_curtis.h'
	63
	64	c arguments
	65	real*8 alsx(nbox_max),alnx(nbox_max),adx(nbox_max),xtemp(nbox_max)
	66
	67	c local variables
	68	integer i, k
	69
	70	c ***********
	71
	72	! write (,) 'intershape xtemp =', xtemp
	73
	74	do 1, k=1,nbox
	75	if (xtemp(k).gt.tmax) then
	76	write (,) ' WARNING ! Tpath,tmax= ',xtemp(k),tmax
	77	xtemp(k) = tmax
	78	endif
	79	if (xtemp(k).lt.tmin) then
	80	write (,) ' WARNING ! Tpath,tmin= ',xtemp(k),tmin
	81	xtemp(k) = tmin
	82	endif
	83
	84	i = 1
	85	do while (i.le.mm)
	86	i = i + 1
	87
	88	if (abs(xtemp(k)-thist(i)) .lt. 1.0d-4) then !evita troubles
	89	alsx(k)=xls1(i,k) !16-july-1996
	90	alnx(k)=xln1(i,k)
	91	adx(k)=xld1(i,k)
	92	goto 1
	93	elseif ( thist(i) .le. xtemp(k) ) then
	94	alsx(k) = (( xls1(i,k)*(thist(i-1)-xtemp(k)) +
	95	@ xls1(i-1,k)*(xtemp(k)-thist(i)) )) / (thist(i-1)-thist(i))
	96	alnx(k) = (( xln1(i,k)*(thist(i-1)-xtemp(k)) +
	97	@ xln1(i-1,k)*(xtemp(k)-thist(i)) )) / (thist(i-1)-thist(i))
	98	adx(k) = (( xld1(i,k)*(thist(i-1)-xtemp(k)) +
	99	@ xld1(i-1,k)*(xtemp(k)-thist(i)) )) / (thist(i-1)-thist(i))
	100	goto 1
	101	end if
	102	end do
	103	write (,)
	104	@ ' error in xtemp(k). it should be between tmin and tmax'
	105	1 continue
	106
	107	return
	108	end
	109	c **********************************************************************
	110	subroutine interstrength (stx, ts, sx, xtemp)
	111	c interpolates the line strength at a temperature xtemp from
	112	c input histogram data.
	113	c **********************************************************************
	114
	115	implicit none
	116
	117	include 'nltedefs.h'
	118	include 'nlte_curtis.h'
	119
	120	c arguments
	121	real*8 stx ! output, total band strength
	122	real*8 ts ! input, temp for stx
	123	real*8 sx(nbox_max) ! output, strength for each box
	124	real*8 xtemp(nbox_max) ! input, temp for sx
	125
	126	c local variables
	127	integer i, k
	128
	129	c ***********
	130
	131	do 1, k=1,nbox
	132	! if(xtemp(k).lt.ts)then
	133	! write(,)'***********************'
	134	! write(,)'mztfsub_overlap/EEEEEEH!',xtemp(k),ts,k
	135	! write(,)'***********************'
	136	! endif
	137	if (xtemp(k).gt.tmax) xtemp(k) = tmax
	138	if (xtemp(k).lt.tmin) xtemp(k) = tmin
	139	i = 1
	140	do while (i.le.mm-1)
	141	i = i + 1
	142	! write(,)'mztfsub_overlap/136',i,xtemp(k),thist(i)
	143	if ( abs(xtemp(k)-thist(i)) .lt. 1.0d-4 ) then
	144	sx(k) = sk1(i,k)
	145	! write(,)'mztfsub_overlap/139',sx(k),k,i
	146	goto 1
	147	elseif ( thist(i) .le. xtemp(k) ) then
	148	sx(k) = ( sk1(i,k)(thist(i-1)-xtemp(k)) + sk1(i-1,k)
	149	@ (xtemp(k)-thist(i)) ) / (thist(i-1)-thist(i))
	150	! write(,)'mztfsub_overlap/144',sx(k),k,i
	151	goto 1
	152	end if
	153	end do
	154	write (,) ' error in xtemp(kr) =', xtemp(k),
	155	@ '. it should be between '
	156	write (,) ' tmin =',tmin, ' and tmax =',tmax
	157	stop
	158	1 continue
	159
	160	stx = 0.d0
	161	if (ts.gt.tmax) ts = dble( tmax )
	162	if (ts.lt.tmin) ts = dble( tmin )
	163	i = 1
	164	do while (i.le.mm-1)
	165	i = i + 1
	166	! write(,)'mztfsub_overlap/160',i,ts,thist(i)
	167	if ( abs(ts-thist(i)) .lt. 1.0d-4 ) then
	168	do k=1,nbox
	169	stx = stx + no(k) * sk1(i,k)
	170	! write(,)'mztfsub_overlap/164',stx
	171	end do
	172	return
	173	elseif ( thist(i) .le. ts ) then
	174	do k=1,nbox
	175	stx = stx + no(k) * (( sk1(i,k)*(thist(i-1)-ts) +
	176	@ sk1(i-1,k)*(ts-thist(i)) )) / (thist(i-1)-thist(i))
	177	! write(,)'mztfsub_overlap/171',stx
	178	end do
	179	! stop
	180	return
	181	end if
	182	end do
	183
	184	return
	185	end
	186	c **********************************************************************
	187	subroutine intz(h,aco2,ap,amr,at, con)
	188	c return interp. concentration, pressure,mixing ratio and temperature
	189	c for a input height h
	190	c **********************************************************************
	191
	192	implicit none
	193	include 'nltedefs.h'
	194	include 'nlte_atm.h'
	195	include 'nlte_curtis.h'
	196
	197	c arguments
	198	real h ! i
	199	real*8 con(nzy) ! i
	200	real*8 aco2, ap, at, amr ! o
	201
	202	c local variables
	203	integer k
	204
	205	c ************
	206
	207	if ( ( h.lt.zy(1) ).and.( h.le.-1.e-5 ) ) then
	208	write (,) ' zp= ',h,' zy(1)= ',zy(1)
	209	stop'from intz: error in interpolation, z < minimum height'
	210	elseif (h.gt.zy(nzy)) then
	211	write (,) ' zp= ',h,' zy(nzy)= ',zy(nzy)
	212	stop'from intz: error in interpolation, z > maximum height'
	213	end if
	214
	215	if (h.eq.zy(nzy)) then
	216	ap = dble( py(nzy) )
	217	aco2= con(nzy)
	218	at = dble( ty(nzy) )
	219	amr = dble( mr(nzy) )
	220	return
	221	end if
	222
	223	do k=1,nzy-1
	224	if( abs( h-zy(k) ).le.( 1.e-5 ) ) then
	225	ap = dble( py(k) )
	226	aco2= con(k)
	227	at = dble( ty(k) )
	228	amr = dble( mr(k) )
	229	return
	230	elseif(h.gt.zy(k).and.h.lt.zy(k+1))then
	231	ap = dble( exp( log(py(k)) + log(py(k+1)/py(k)) *
	232	@ (h-zy(k)) / (zy(k+1)-zy(k)) ) )
	233	aco2 = exp( log(con(k)) + log( con(k+1)/con(k) ) *
	234	@ (h-zy(k)) / (zy(k+1)-zy(k)) )
	235	at = dble( ty(k)+(ty(k+1)-ty(k))*(h-zy(k))/
	236	@ (zy(k+1)-zy(k)) )
	237	amr = dble( mr(k)+(mr(k+1)-mr(k))*(h-zy(k))/
	238	@ (zy(k+1)-zy(k)) )
	239	return
	240	end if
	241	end do
	242
	243	return
	244	end
	245
	246	c *******************************************************************
	247
	248	subroutine rhist (factor_comp)
	249
	250	c reads histogram data arrays created by ~/spectral/his.for
	251	c malv nov-98 add average distance between lines for overlapp
	252
	253	c *******************************************************************
	254
	255
	256	implicit none
	257
	258	include 'nltedefs.h'
	259	include 'nlte_curtis.h'
	260	include 'datafile.h'
	261
	262	c arguments
	263	real factor_comp
	264
	265	c local variables
	266	integer j, r
	267	real*8 sk1_aux, xls1_aux, xln1_aux, xld1_aux,weight,nu0
	268	character tonto*50
	269
	270	c formats
	271	! 100 format(80a1) ! Esto es si fuese byte dummy(80)
	272	100 format(a80) ! Esto es si fuese character dummy*80
	273	150 format(a50) ! Esto es si fuese character dummy*80
	274
	275	c ***************
	276
	277	open(unit=3,
	278	$ file=trim(datafile)//'/NLTEDAT/'
	279	$ //hisfile(1:len_trim(hisfile)),status='old')
	280	!read(3,100) dummy
	281	read(3,150) tonto
	282	read(3,*) weight
	283	read(3,*) mm
	284	read(3,*) nu0
	285	read(3,*) nbox
	286	!read(3,'(a)') dumm
	287	read(3,'(a)') tonto
	288
	289	if ( nbox .gt. nbox_max ) then
	290	write (,) ' nbox too large in input file ', hisfile
	291	stop ' Check maximum number nbox_max in mz1d.par '
	292	endif
	293	do 1 j=1,mm ! for each temperature
	294	read(3,*) thist(j)
	295	do r=1,nbox ! for each box
	296	read(3,*) no(r), sk1(j,r), xls1(j,r),xln1(j,r),xld1(j,r),
	297	@ dist(r)
	298	c xld1(j,r)=xld1(j,r)*0.83255 !0.83255=sqrt(log(2))
	299	enddo
	300	1 continue
	301	tmax=thist(1)
	302	tmin=thist(mm)
	303
	304	!close(unit=3,dispose='save')
	305	close(unit=3)
	306
	307
	308	do 2 j=1,mm
	309	do r=1,nbox
	310	sk1(j,r) = sk1(j,r) * factor_comp
	311	enddo
	312	2 continue
	313
	314
	315	return
	316	end
	317
	318	c ****************************************************************
	319	subroutine interh( sx, alsx, alnx, adx, xtemp, xtdop )
	320
	321	c interpolates a histogram at temperature xtemp from input histogr
	322
	323	c jan 98 malv version para mz1d basada en el inth de solar10:
	324	c mz5/curtis/mztfsub_solar.f
	325	c ****************************************************************
	326
	327	implicit none
	328
	329	include 'nltedefs.h'
	330	include 'nlte_curtis.h'
	331
	332	c arguments
	333	real*8 sx(nbox_max) ! o
	334	real*8 xtemp ! i
	335	real*8 alsx(nbox_max) ! o
	336	real*8 alnx(nbox_max) ! o
	337	real*8 adx(nbox_max) ! o
	338	real*8 xtdop ! i
	339
	340	c local variables
	341	integer i, j, k
	342
	343	c ************
	344
	345
	346	if (xtemp.gt.thist(1)) then
	347	! write (,) ' xtemp-path, thist(1)max: ',
	348	! @ xtemp,thist(1)
	349	xtemp = thist(1)
	350	elseif (xtemp.lt.thist(nhist)) then
	351	! write (,) ' xtemp-path, thist(nhist)min: ',
	352	! @ xtemp,thist(nhist)
	353	xtemp = thist(nhist)
	354	end if
	355
	356	i=0
	357	1 i=i+1
	358	if (abs(xtemp-thist(i)) .lt. 1.0d-4) goto 2
	359	if (thist(i).lt.xtemp) goto 2
	360	goto 1
	361	2 j=i
	362
	363	! write (,) 'InterH/ j, xtemp, th(1),th(nh)',
	364	! @ j, xtemp, thist(1),thist(nhist)
	365
	366	if (j.gt.1) then
	367	do k=1,nbox
	368	sx(k) = ( sk1(j,k) * (thist(j-1)-xtemp) + sk1(j-1,k) *
	369	@ (xtemp-thist(j)) ) / (thist(j-1)-thist(j))
	370	enddo
	371	elseif (j.eq.1) then
	372	do k=1,nbox
	373	sx(k) = sk1(1,k)
	374	enddo
	375	endif
	376
	377	if (xtdop.gt.thist(1)) then
	378	! write (,) ' xtdop-path, thist(1)max: ',
	379	! @ xtdop,thist(1)
	380	xtdop = thist(1)
	381	elseif (xtdop.lt.thist(nhist)) then
	382	! write (,) ' xtdop-path, thist(nhist)min: ',
	383	! @ xtdop,thist(nhist)
	384	xtdop = thist(nhist)
	385	end if
	386
	387	i=0
	388	4 i=i+1
	389	if (abs(xtdop-thist(i)) .lt. 1.0d-4) goto 5
	390	if (thist(i).lt.xtdop) goto 5
	391	goto 4
	392	5 j=i
	393
	394	! write (,) 'InterH/ j, xtdop',
	395	! @ j, xtdop
	396
	397	if (j.gt.1) then
	398	do k=1,nbox
	399	alsx(k) = ( xls1(j,k) * (thist(j-1)-xtdop) + xls1(j-1,k)*
	400	@ (xtdop-thist(j)) ) / (thist(j-1)-thist(j))
	401	alnx(k) = ( xln1(j,k) * (thist(j-1)-xtdop) + xln1(j-1,k)*
	402	@ (xtdop-thist(j)) ) / (thist(j-1)-thist(j))
	403	adx(k) = ( xld1(j,k) * (thist(j-1)-xtdop) + xld1(j-1,k)*
	404	@ (xtdop-thist(j)) ) / (thist(j-1)-thist(j))
	405	enddo
	406	elseif (j.eq.1) then
	407	do k=1,nbox
	408	alsx(k) = xls1(1,k)
	409	alnx(k) = xln1(1,k)
	410	adx(k) = xld1(j,k)
	411	enddo
	412	endif
	413
	414	c end
	415	return
	416	end
	417
	418	c **********************************************************************
	419	real*8 function we(ig,me,pe,pl, idummy, nt_local,p_local, Desp, wsL)
	420	c icls=5 -->para mztf
	421	c icls=1,2,3-->para fot, line shape (v=1,l=2,d=3) (only use if wr=2)
	422	c calculates an approximate equivalent width for an error estimate.
	423	c
	424	c ioverlap = 0 ....... no correction for overlaping
	425	c 1 ....... "lisat" first correction (see overlap_box.
	426	c 2 ....... " " " plus "supersaturation"
	427
	428	c idummy=0 do nothing
	429	c 1 write out some values for diagnostics
	430	c 2 correct the Strong Lorentz behaviour for SZA>90
	431	c 3 casos 1 & 2
	432
	433	c malv nov-98 add overlaping's corrections
	434	c **********************************************************************
	435
	436	implicit none
	437
	438	include 'nltedefs.h'
	439	include 'nlte_curtis.h'
	440	include 'tcr_15um.h'
	441
	442	c arguments
	443	integer ig ! ADDED FOR TRACEBACK
	444	real*8 me ! I. path's absorber amount
	445	real*8 pe ! I. path's presion total
	446	real*8 pl ! I. path's partial pressure of CO2
	447	real*8 nt_local ! I. needed for strong limit of Lorentz profil
	448	real*8 p_local ! I. " " "
	449	integer idummy ! I. indica varias opciones
	450	real*8 wsL ! O. need this for strong Lorentz correction
	451	real*8 Desp ! I. need this for strong Lorentz correction
	452
	453	c local variables
	454	integer i
	455	real*8 y,x,wl,wd
	456	real*8 cn(0:7),dn(0:7)
	457	real*8 pi, xx
	458	real*8 f_sat_box
	459	real*8 dv_sat_box, dv_corte_box
	460	real*8 area_core_box, area_wing_box
	461	real*8 wlgood , parentesis , xlor
	462	real*8 wsl_grad
	463
	464
	465	c data blocks
	466	data cn/9.99998291698d-1,-3.53508187098d-1,9.60267807976d-2,
	467	@ -2.04969011013d-2,3.43927368627d-3,-4.27593051557d-4,
	468	@ 3.42209457833d-5,-1.28380804108d-6/
	469	data dn/1.99999898289,5.774919878d-1,-5.05367549898d-1,
	470	@ 8.21896973657d-1,-2.5222672453,6.1007027481,
	471	@ -8.51001627836,4.6535116765/
	472
	473	c ***********
	474
	475	c equivalent width of atmospheric line.
	476
	477	pi = acos(-1.d0)
	478
	479	if ( idummy.gt.9 )
	480	@ write (,) ' S, m, alsa, pp =', ka(kr), me, alsa(kr), pl
	481
	482	y=ka(kr)*me
	483	! x=y/(2.0pi(alsa(kr)pl+alna(kr)(pe-pl)))
	484	x=y/(2.0d0pi alsa(kr)pl) !+alna(kr)(pe-pl)))
	485
	486	! Strong limit of Lorentz profile: WL = 2 SQRT( S * m * alsa*pl )
	487	! Para anular esto, comentar las siguientes 5 lineas
	488	! if ( x .gt. 1.e6 ) then
	489	! wl = 2.0sqrt( y alsa(kr)*pl )
	490	! else
	491	! wl=y/sqrt(1.0d0+pi*x/2.0d0)
	492	! endif
	493
	494	wl=y/sqrt(1.0d0+pi*x/2.0d0)
	495
	496	if (wl .le. 0.d0) then
	497	write(,)'mztfsub_overlap/496',ig,y,ka(kr),me,kr
	498	stop'WE/Lorentz EQW zero or negative!/498'!,ig
	499	endif
	500
	501	if ( idummy.gt.9 )
	502	@ write (,) ' y, x =', y, x
	503
	504	xlor = x
	505	if ( (idummy.eq.2 .or. idummy.eq.12) .and. xlor.gt.1e5 ) then
	506	! en caso que estemos en el regimen
	507	! Strong Lorentz y la presion local
	508	! vaya disminuyendo, corregimos la EQW
	509	! con un gradiente analitico (notebook)
	510	wsL = 2.0sqrt( y alsa(kr)*pl )
	511	wsl_grad = - 2.d0 * ka(kr)alsa(kr) nt_local*p_local / wsL
	512	wlgood = w_strongLor_prev(kr) + wsl_grad * Desp
	513	if (idummy.eq.12)
	514	@ write (,) ' W(wrong), W_SL, W_SL prev, W_SL corrected=',
	515	@ wl, wsL, w_strongLor_prev(kr), wlgood
	516	wl = wlgood
	517	endif
	518	! wsL = wl pero esto no lo hacemos todavia, porque necesitamos
	519	! el valor que ahora mismo tiene wsL para corregir la
	520	! expresion R&W below
	521
	522	! write (,) 'WE arguments me,pe,pl =', me,pe,pl
	523	! write (,) 'WE/ wl,ka(kr),alsa(kr) =',
	524	! @ wl, ka(kr),alsa(kr)
	525
	526
	527	!>>>>>>>
	528	500 format (a,i3,3(2x,1pe15.8))
	529	600 format (a,2(2x,1pe16.9))
	530	700 format (a,3(1x,1pe16.9))
	531	! if (kr.eq.8 .or. kr.eq.13) then
	532	! write (*,500) 'WE/kr,m,pt,pl=', kr, me, pe, pl
	533	! write (*,700) ' /aln,als,d_x=', alna(kr),alsa(kr),
	534	! @ 2.0pi( alsa(kr)pl + alna(kr)(pe-pl) )
	535	! write (,600) ' /alsap_CO2, alna*p_n2 :',
	536	! @ alsa(kr)pl, alna(kr)(pe-pl)
	537	! write (,600) ' ap, S =',
	538	! @ alsa(kr)pl + alna(kr)(pe-pl) , ka(kr)
	539	! write (,600) ' /Sm, x =', y, x
	540	! write (*,600) ' /aprox, WL =',
	541	! @ 2.sqrt( y( alsa(kr)pl+alna(kr)(pe-pl) ) ), WL
	542	! endif
	543	!
	544	! corrections to lorentz eqw due to overlaping and super-saturation
	545	!
	546
	547	i_supersat = 0
	548
	549	if ( icls.eq.5 .and. ioverlap.gt.0 ) then
	550	! for the moment, only consider overlaping for mztf.f, not fot.f
	551
	552	! definition of saturation in the lisat model
	553	!
	554	asat_box = 0.99d0
	555	f_sat_box = 2.d0 * x
	556	xx = f_sat_box / log( 1./(1-asat_box) )
	557	if ( xx .lt. 1.0d0 ) then
	558	dv_sat_box = 0.0d0
	559	asat_box = 1.0d0 - exp( - f_sat_box )
	560	else
	561	dv_sat_box = alsa(kr) * sqrt( xx - 1.0d0 )
	562	! approximation: only use of alsa in mars and venus
	563	endif
	564
	565	! area of saturated line
	566	!
	567	area_core_box = 2.d0 * dv_sat_box * asat_box
	568	area_wing_box = 0.5d0 * ( wl - area_core_box )
	569	dv_corte_box = dv_sat_box + 2.d0*area_wing_box/asat_box
	570
	571	! super-saturation or simple overlaping?
	572	!
	573	! i_supersat = 0
	574	xx = dv_sat_box - ( 0.5d0 * dist(kr) )
	575	if ( xx .ge. 0.0 ! definition of supersaturation
	576	@ .and. dv_sat_box .gt. 0.0 ! definition of saturation
	577	@ .and. (dist(kr).gt.0.0) ) ! box contains more than 1 line
	578	@ ! and not too far apart
	579	@ then
	580
	581	i_supersat = 1
	582
	583	else
	584	! no super-saturation, then use "lisat + first correction", i.e.,
	585	! correct for line products
	586	!
	587
	588	wl = wl
	589
	590	endif
	591
	592	end if ! end of overlaping loop
	593
	594	if (icls.eq.2) then
	595	we = wl
	596	return
	597	endif
	598
	599	cc doppler limit:
	600	if ( idummy.gt.9 )
	601	@ write (,) ' Sm, alf_dop =', y, alda(kr)sqrt(pi)
	602
	603	x = y / (alda(kr)*sqrt(pi))
	604	if ( x.lt.1.e-10 ) then ! to avoid underflow
	605	wd = y
	606	else
	607	wd=alda(kr)sqrt(4.0pix2(1.0+log(1.0+x))/(4.0+pix*2))
	608	endif
	609	if ( idummy.gt.9 )
	610	@ write (,) ' wd =', wd
	611
	612	cc doppler weak limit
	613	c wd = ka(kr) * me
	614
	615	cc good doppler
	616	if(icls.eq.5) then !para mztf
	617	!write (,) 'para mztf, icls=',icls
	618	if (x.lt.5.) then
	619	wd = 0.d0
	620	do i=0,7
	621	wd = wd + cn(i) * x**i
	622	end do
	623	wd = alda(kr) * x * sqrt(pi) * wd
	624	elseif (x.gt.5.) then
	625	wd = 0.d0
	626	do i=0,7
	627	wd = wd + dn(i) / (log(x))**i
	628	end do
	629	wd = alda(kr) * sqrt(log(x)) * wd
	630	else
	631	stop ' x should not be less than zero'
	632	end if
	633	end if
	634
	635
	636	if ( i_supersat .eq. 0 ) then
	637
	638	parentesis = wl2+wd2-(wdwl/y)*2
	639	! changed +(wdwl/y)*2 to -...14-3-84
	640
	641	if ( parentesis .lt. 0.0 ) then
	642	if ((idummy.eq.2 .or. idummy.eq.12) .and. xlor.gt.1e5) then
	643	parentesis = wl2+wd2-(wdwsL/y)*2
	644	! este cambio puede ser necesario cuando se hace
	645	! correccion Strong Lor, para evitar valores
	646	! negativos del parentesis en sqrt( )
	647	else
	648	stop ' WE/ Error en las EQW wl,wl,y '
	649	endif
	650	endif
	651
	652	we = sqrt( parentesis )
	653	! write (,) ' from we: xdop,alda,wd', sngl(x),alda(kr),sngl(wd)
	654	! write (,) ' from we: we', we
	655
	656	else
	657
	658	we = wl
	659	! if there is supersaturation we can ignore wd completely;
	660	! mztf.f will compute the eqw of the whole box afterwards
	661
	662	endif
	663
	664	if (icls.eq.3) we = wd
	665
	666	if ( idummy.gt.9 )
	667	@ write (,) ' wl,wd,w =', wl,wd,we
	668
	669	wsL = wl
	670
	671	return
	672	end
	673
	674
	675	c **********************************************************************
	676	real*8 function simrul(a,b,fsim,c,acc)
	677	c adaptively integrates fsim from a to b, within the criterion acc.
	678	c **********************************************************************
	679
	680	implicit none
	681
	682	real*8 res,a,b,g0,g1,g2,g3,g4,d,a0,a1,a2,h,x,acc,c,fsim
	683	real*8 s1(70),s2(70),s3(70)
	684	real*8 c1, c2
	685	integer*4 m,n,j
	686
	687	res=0.
	688	c=0.
	689	m=0
	690	n=0
	691	j=30
	692	g0=fsim(a)
	693	g2=fsim((a+b)/2.)
	694	g4=fsim(b)
	695	a0=(b-a)(g0+4.0g2+g4)/2.0
	696	1 d=2.0**n
	697	h=(b-a)/(4.0*d)
	698	x=a+(4.0m+1.0)h
	699	g1=fsim(x)
	700	g3=fsim(x+2.0*h)
	701	a1=h(g0+4.0g1+g2)
	702	a2=h(g2+4.0g3+g4)
	703	if ( abs(a1+a2-a0).gt.(acc/d)) goto 2
	704	res=res+(16.0*(a1+a2)-a0)/45.0
	705	m=m+1
	706	c=a+m*(b-a)/d
	707	6 if (m.eq.(2*(m/2))) goto 4
	708	if ((m.ne.1).or.(n.ne.0)) goto 5
	709	8 simrul=res
	710	return
	711	2 m=2*m
	712	n=n+1
	713	if (n.gt.j) goto 3
	714	a0=a1
	715	s1(n)=a2
	716	s2(n)=g3
	717	s3(n)=g4
	718	g4=g2
	719	g2=g1
	720	goto 1
	721	3 c1=c-(b-a)/d
	722	c2=c+(b-a)/d
	723	write(2,7) c1,c,c2,fsim(c1),fsim(c),fsim(c2)
	724	write(*,7) c1,c,c2,fsim(c1),fsim(c),fsim(c2)
	725	7 format(2x,17hsimrule fails at ,/,3e15.6,/,3e15.6)
	726	goto 8
	727	5 a0=s1(n)
	728	g0=g4
	729	g2=s2(n)
	730	g4=s3(n)
	731	goto 1
	732	4 m=m/2
	733	n=n-1
	734	goto 6
	735	end
	736
	737	c **********************************************************************
	738	subroutine findw(ig,iirw,idummy,c1,p1, Desp, wsL)
	739	c this routine sets up accuracy criteria and calls simrule between limit
	740	c that depend on the number of atmospheric and cell paths. it gives eqw.
	741
	742	c Add correction for EQW in Strong Lorentz regime and SZA>90
	743	c **********************************************************************
	744
	745	implicit none
	746	include 'nltedefs.h'
	747	include 'nlte_curtis.h'
	748
	749	c arguments
	750	integer ig ! ADDED FOR TRACEBACK
	751	integer iirw
	752	integer idummy ! I. indica varias opciones
	753	real*8 c1 ! I. needed for strong limit of Lorentz profil
	754	real*8 p1 ! I. " " "
	755	real*8 wsL ! O. need this for strong Lorentz correction
	756	real*8 Desp ! I. need this for strong Lorentz correction
	757
	758	c local variables
	759	real*8 ept,eps,xa
	760	real*8 acc, c
	761	real*8 we
	762	real*8 f, fi, simrul
	763
	764	external f,fi
	765
	766	c ******** ******* *******
	767
	768	if(icls.eq.5) then !para mztf
	769	! if(ig.eq.1682)write(,)'mztfsub_overlap/768',ua(kr),iirw
	770	if (iirw.eq.2) then !iirw=icf=2 ==> we use the w&r formula
	771	w = we(ig,ua(kr),pt,pp, idummy, c1,p1, Desp, wsL )
	772	return
	773	end if
	774	ept=we(ig,ua(kr),pt,pp, idummy,c1,p1, Desp, wsL)
	775	else !para fot
	776	if (iirw.eq.2) then ! icf=2 ==> we use the w&r formula
	777	w = we(ig,sl_ua,pt,pp, idummy,c1,p1, Desp, wsL)
	778	return
	779	end if
	780	ept=we(ig,sl_ua,pt,pp, idummy,c1,p1, Desp, wsL)
	781	end if
	782
	783	c the next block is a modification to avoid nul we.
	784	c this situation appears for weak lines and low path temperature, but
	785	c there is not any loss of accuracy. first july 1986
	786	if (ept.eq.0.) then ! for weak lines sometimes we=0
	787	ept=1.0e-18
	788	write (,) 'ept =',ept
	789	write (,) 'from we: we=0.0'
	790	return
	791	end if
	792
	793	acc = 4.d0
	794	acc = 10.d0**(-acc)
	795
	796	eps = acc * ept !accuracy 10-4 atmospheric eqw.
	797	xa=0.5*ept/f(0.d0) !width of doppler shifted atmospheric line.
	798	w=2.0*(simrul(0.0d0,xa,f,c,eps)+simrul(0.1d0,1.0/xa,fi,c,eps))
	799	!no shift.
	800
	801	return
	802	end
	803
	804
	805	c **********************************************************************
	806	double precision function fi(y)
	807	c returns the value of f(1/y)
	808	c **********************************************************************
	809
	810	implicit none
	811	real*8 f, y
	812
	813	fi=f(1.0/y)/y**2
	814	return
	815	end
	816
	817
	818	c **********************************************************************
	819	double precision function f(nu)
	820	c calculates 1-exp(-k(nu)u) for all series paths or combinations thereof
	821	c **********************************************************************
	822
	823	implicit none
	824	include 'nltedefs.h'
	825	include 'nlte_curtis.h'
	826
	827	double precision tra,xa,ya,za,yy,nu
	828	double precision voigtf
	829	tra=1.0d0
	830
	831	yy=1.0d0/alda(kr)
	832	xa=nu*yy
	833	ya= ( alsa(kr)pp + alna(kr)(pt-pp) ) * yy
	834	za=ka(kr)*yy
	835
	836	if(icls.eq.5) then !para mztf
	837	! write (,) 'icls=',icls
	838	tra=zaua(kr)voigtf(sngl(xa),sngl(ya))
	839	else
	840	tra=zasl_uavoigtf(sngl(xa),sngl(ya))
	841	end if
	842
	843	if (tra.gt.50.0) then
	844	tra=1.0 !2.0e-22 overflow cut-off.
	845	else if (tra.gt.1.0e-4) then
	846	tra=1.0-exp(-tra)
	847	end if
	848
	849	f=tra
	850	return
	851	end
	852
	853	c **********************************************************************
	854	double precision function voigtf(x1,y)
	855	c computes voigt function for any value of x1 and any +ve value of y.
	856	c where possible uses modified lorentz and modified doppler approximatio
	857	c otherwise uses a rearranged rybicki routine.
	858	c c(n) = exp(-(n/h)*2)/(pisqrt(pi)), with h = 2.5 .
	859	c accurate to better than 1 in 10000.
	860	c **********************************************************************
	861
	862	implicit none
	863
	864	real x1, y
	865	real x, xx, xxyy, xh,xhxh, yh,yhyh, f1,f2, p, q, xn,xnxn, voig
	866
	867	real*8 b,g0,g1,g2,g3,g4,d1,d2,d3,d4,c
	868	integer*4 n
	869
	870	dimension c(10)
	871	complex xp,xpp,z
	872
	873	data c(1)/0.15303405/
	874	data c(2)/0.94694928e-1/
	875	data c(3)/0.42549174e-1/
	876	data c(4)/0.13882935e-1/
	877	data c(5)/0.32892528e-2/
	878	data c(6)/0.56589906e-3/
	879	data c(7)/0.70697890e-4/
	880	data c(8)/0.64135678e-5/
	881	data c(9)/0.42249221e-6/
	882	data c(10)/0.20209868e-7/
	883
	884	x=abs(x1)
	885	if (x.gt.7.2) goto 1
	886	if ((y+x*0.3).gt.5.4) goto 1
	887	if (y.gt.0.01) goto 3
	888	if (x.lt.2.1) goto 2
	889	goto 3
	890	c here uses modified lorentz approx.
	891	1 xx=x*x
	892	xxyy=xx+y*y
	893	b=xx/xxyy
	894	voigtf=y(1.+(2.b-0.5+(0.75-(9.-12.b)b)/xxyy)/
	895	* xxyy)/(xxyy*3.141592654)
	896	return
	897	c here uses modified doppler approx.
	898	2 xx=x*x
	899	voigtf=0.56418958exp(-xx)(1.-y(1.-0.5y)(1.1289-xx(1.1623+
	900	* xx(0.080812+xx(0.13854-xx(0.033605-0.0073972xx))))))
	901	return
	902	c here uses a rearranged rybicki routine.
	903	3 xh=2.5*x
	904	xhxh=xh*xh
	905	yh=2.5*y
	906	yhyh=yh*yh
	907	f1=xhxh+yhyh
	908	f2=f1-0.5*yhyh
	909	if (y.lt.0.1) goto 20
	910	p=-y7.8539816 !7.8539816=2.5pi
	911	q=x*7.8539816
	912	xpp=cmplx(p,q)
	913	z=cexp(xpp)
	914	d1=xh*aimag(z)
	915	d2=-d1
	916	d3=yh*(1.-real(z))
	917	d4=-d3+2.*yh
	918	voig=0.17958712*(d1+d3)/f1
	919	goto 30
	920	20 p=x*7.8539816
	921	q=y*7.8539816
	922	xp=cmplx(p,q)
	923	z=ccos(xp)
	924	d1=xh*aimag(z)
	925	d2=-d1
	926	d3=yh*(1.-real(z))
	927	d4=-d3+2.*yh
	928	voig=0.56418958exp(yy-xx)cos(2.xy)+0.17958712*(d1+d3)/f1
	929	30 xn=0.
	930	do 55 n=1,10,2
	931	xn=xn+1.
	932	xnxn=xn*xn
	933	g1=xh-xn
	934	g2=g1*(xh+xn)
	935	g3=0.5g2g2
	936	voig=voig+c(n)(d2(g2+yhyh)+d4(f1+xnxn))/(g3+yhyh(f2+xnxn))
	937	xn=xn+1.
	938	xnxn=xn*xn
	939	g1=xh-xn
	940	g2=g1*(xh+xn)
	941	g3=0.5g2g2
	942	voig=voig+c(n+1)(d1(g2+yhyh)+d3*(f1+xnxn))/
	943	@ (g3+yhyh*(f2+xnxn))
	944	55 continue
	945	voigtf=voig
	946	return
	947	end

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

Download in other formats:

Original Format