Context Navigation

nlte_aux.F @ 517

Last change on this file since 517 was 498, checked in by emillour, 14 years ago
Mars GCM: Cleanup of the NLTE routines which have been packed together to limit the number of files. Also enforced that file names are non-capitalized (needed by the create_make_gcm scripts to better evaluate dependencies when building the makefile). FGG+EM
File size: 78.5 KB

Rev	Line
[498]	1	c***********************************************************************
	2	c File with all subroutines required by mztf
	3	c Subroutines previously included in mztfsub_overlap.F
	4	c
	5	c jan 98 malv basado en mztfsub_solar
	6	c jul 2011 malv+fgg adapted to LMD-MGCM
	7	c
	8	c contiene:
	9	c initial
	10	c intershape
	11	c interstrength
	12	c intz
	13	c rhist
	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 'nlte_paramdef.h'
	31	include 'nlte_commons.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 'nlte_paramdef.h'
	62	include 'nlte_commons.h'
	63
	64	c arguments
	65	real*8 alsx(nbox_max),alnx(nbox_max),adx(nbox_max),
	66	& xtemp(nbox_max)
	67
	68	c local variables
	69	integer i, k
	70
	71	c ***********
	72
	73	! write (,) 'intershape xtemp =', xtemp
	74
	75	do 1, k=1,nbox
	76	if (xtemp(k).gt.tmax) then
	77	write (,) ' WARNING ! Tpath,tmax= ',xtemp(k),tmax
	78	xtemp(k) = tmax
	79	endif
	80	if (xtemp(k).lt.tmin) then
	81	write (,) ' WARNING ! Tpath,tmin= ',xtemp(k),tmin
	82	xtemp(k) = tmin
	83	endif
	84
	85	i = 1
	86	do while (i.le.mm)
	87	i = i + 1
	88
	89	if (abs(xtemp(k)-thist(i)) .lt. 1.0d-4) then !evita troubles
	90	alsx(k)=xls1(i,k) !16-july-1996
	91	alnx(k)=xln1(i,k)
	92	adx(k)=xld1(i,k)
	93	goto 1
	94	elseif ( thist(i) .le. xtemp(k) ) then
	95	alsx(k) = (( xls1(i,k)*(thist(i-1)-xtemp(k)) +
	96	@ xls1(i-1,k)*(xtemp(k)-thist(i)) )) /
	97	$ (thist(i-1)-thist(i))
	98	alnx(k) = (( xln1(i,k)*(thist(i-1)-xtemp(k)) +
	99	@ xln1(i-1,k)*(xtemp(k)-thist(i)) )) /
	100	$ (thist(i-1)-thist(i))
	101	adx(k) = (( xld1(i,k)*(thist(i-1)-xtemp(k)) +
	102	@ xld1(i-1,k)*(xtemp(k)-thist(i)) )) /
	103	$ (thist(i-1)-thist(i))
	104	goto 1
	105	end if
	106	end do
	107	write (,)
	108	@ ' error in xtemp(k). it should be between tmin and tmax'
	109	1 continue
	110
	111	return
	112	end
	113	c **********************************************************************
	114	subroutine interstrength (stx, ts, sx, xtemp)
	115	c interpolates the line strength at a temperature xtemp from
	116	c input histogram data.
	117	c **********************************************************************
	118
	119	implicit none
	120
	121	include 'nlte_paramdef.h'
	122	include 'nlte_commons.h'
	123
	124	c arguments
	125	real*8 stx ! output, total band strength
	126	real*8 ts ! input, temp for stx
	127	real*8 sx(nbox_max) ! output, strength for each box
	128	real*8 xtemp(nbox_max) ! input, temp for sx
	129
	130	c local variables
	131	integer i, k
	132
	133	c ***********
	134
	135	do 1, k=1,nbox
	136	! if(xtemp(k).lt.ts)then
	137	! write(,)'***********************'
	138	! write(,)'mztfsub_overlap/EEEEEEH!',xtemp(k),ts,k
	139	! write(,)'***********************'
	140	! endif
	141	if (xtemp(k).gt.tmax) xtemp(k) = tmax
	142	if (xtemp(k).lt.tmin) xtemp(k) = tmin
	143	i = 1
	144	do while (i.le.mm-1)
	145	i = i + 1
	146	! write(,)'mztfsub_overlap/136',i,xtemp(k),thist(i)
	147	if ( abs(xtemp(k)-thist(i)) .lt. 1.0d-4 ) then
	148	sx(k) = sk1(i,k)
	149	! write(,)'mztfsub_overlap/139',sx(k),k,i
	150	goto 1
	151	elseif ( thist(i) .le. xtemp(k) ) then
	152	sx(k) = ( sk1(i,k)(thist(i-1)-xtemp(k)) + sk1(i-1,k)
	153	@ (xtemp(k)-thist(i)) ) / (thist(i-1)-thist(i))
	154	! write(,)'mztfsub_overlap/144',sx(k),k,i
	155	goto 1
	156	end if
	157	end do
	158	write (,) ' error in xtemp(kr) =', xtemp(k),
	159	@ '. it should be between '
	160	write (,) ' tmin =',tmin, ' and tmax =',tmax
	161	stop
	162	1 continue
	163
	164	stx = 0.d0
	165	if (ts.gt.tmax) ts = dble( tmax )
	166	if (ts.lt.tmin) ts = dble( tmin )
	167	i = 1
	168	do while (i.le.mm-1)
	169	i = i + 1
	170	! write(,)'mztfsub_overlap/160',i,ts,thist(i)
	171	if ( abs(ts-thist(i)) .lt. 1.0d-4 ) then
	172	do k=1,nbox
	173	stx = stx + no(k) * sk1(i,k)
	174	! write(,)'mztfsub_overlap/164',stx
	175	end do
	176	return
	177	elseif ( thist(i) .le. ts ) then
	178	do k=1,nbox
	179	stx = stx + no(k) * (( sk1(i,k)*(thist(i-1)-ts) +
	180	@ sk1(i-1,k)*(ts-thist(i)) )) / (thist(i-1)-thist(i))
	181	! write(,)'mztfsub_overlap/171',stx
	182	end do
	183	! stop
	184	return
	185	end if
	186	end do
	187
	188	return
	189	end
	190
	191
	192	c **********************************************************************
	193	subroutine intz(h,aco2,ap,amr,at, con)
	194	c return interp. concentration, pressure,mixing ratio and temperature
	195	c for a input height h
	196	c **********************************************************************
	197
	198	implicit none
	199	include 'nlte_paramdef.h'
	200	include 'nlte_commons.h'
	201
	202	c arguments
	203	real h ! i
	204	real*8 con(nzy) ! i
	205	real*8 aco2, ap, at, amr ! o
	206
	207	c local variables
	208	integer k
	209
	210	c ************
	211
	212	if ( ( h.lt.zy(1) ).and.( h.le.-1.e-5 ) ) then
	213	write (,) ' zp= ',h,' zy(1)= ',zy(1)
	214	stop'from intz: error in interpolation, z < minimum height'
	215	elseif (h.gt.zy(nzy)) then
	216	write (,) ' zp= ',h,' zy(nzy)= ',zy(nzy)
	217	stop'from intz: error in interpolation, z > maximum height'
	218	end if
	219
	220	if (h.eq.zy(nzy)) then
	221	ap = dble( py(nzy) )
	222	aco2= con(nzy)
	223	at = dble( ty(nzy) )
	224	amr = dble( mr(nzy) )
	225	return
	226	end if
	227
	228	do k=1,nzy-1
	229	if( abs( h-zy(k) ).le.( 1.e-5 ) ) then
	230	ap = dble( py(k) )
	231	aco2= con(k)
	232	at = dble( ty(k) )
	233	amr = dble( mr(k) )
	234	return
	235	elseif(h.gt.zy(k).and.h.lt.zy(k+1))then
	236	ap = dble( exp( log(py(k)) + log(py(k+1)/py(k)) *
	237	@ (h-zy(k)) / (zy(k+1)-zy(k)) ) )
	238	aco2 = exp( log(con(k)) + log( con(k+1)/con(k) ) *
	239	@ (h-zy(k)) / (zy(k+1)-zy(k)) )
	240	at = dble( ty(k)+(ty(k+1)-ty(k))*(h-zy(k))/
	241	@ (zy(k+1)-zy(k)) )
	242	amr = dble( mr(k)+(mr(k+1)-mr(k))*(h-zy(k))/
	243	@ (zy(k+1)-zy(k)) )
	244	return
	245	end if
	246	end do
	247
	248	return
	249	end
	250
	251
	252
	253	c **********************************************************************
	254	real*8 function we(ig,me,pe,plaux,idummy,nt_local,p_local,
	255	$ Desp,wsL)
	256	c icls=5 -->para mztf
	257	c icls=1,2,3-->para fot, line shape (v=1,l=2,d=3) (only use if wr=2)
	258	c calculates an approximate equivalent width for an error estimate.
	259	c
	260	c ioverlap = 0 ....... no correction for overlaping
	261	c 1 ....... "lisat" first correction (see overlap_box.
	262	c 2 ....... " " " plus "supersaturation"
	263
	264	c idummy=0 do nothing
	265	c 1 write out some values for diagnostics
	266	c 2 correct the Strong Lorentz behaviour for SZA>90
	267	c 3 casos 1 & 2
	268
	269	c malv nov-98 add overlaping's corrections
	270	c **********************************************************************
	271
	272	implicit none
	273
	274	include 'nlte_paramdef.h'
	275	include 'nlte_commons.h'
	276
	277	c arguments
	278	integer ig ! ADDED FOR TRACEBACK
	279	real*8 me ! I. path's absorber amount
	280	real*8 pe ! I. path's presion total
	281	real*8 plaux ! I. path's partial pressure of CO2
	282	real*8 nt_local ! I. needed for strong limit of Lorentz profil
	283	real*8 p_local ! I. " " "
	284	integer idummy ! I. indica varias opciones
	285	real*8 wsL ! O. need this for strong Lorentz correction
	286	real*8 Desp ! I. need this for strong Lorentz correction
	287
	288	c local variables
	289	integer i
	290	real*8 y,x,wl,wd
	291	real*8 cn(0:7),dn(0:7)
	292	real*8 pi, xx
	293	real*8 f_sat_box
	294	real*8 dv_sat_box, dv_corte_box
	295	real*8 area_core_box, area_wing_box
	296	real*8 wlgood , parentesis , xlor
	297	real*8 wsl_grad
	298
	299
	300	c data blocks
	301	data cn/9.99998291698d-1,-3.53508187098d-1,9.60267807976d-2,
	302	@ -2.04969011013d-2,3.43927368627d-3,-4.27593051557d-4,
	303	@ 3.42209457833d-5,-1.28380804108d-6/
	304	data dn/1.99999898289,5.774919878d-1,-5.05367549898d-1,
	305	@ 8.21896973657d-1,-2.5222672453,6.1007027481,
	306	@ -8.51001627836,4.6535116765/
	307
	308	c ***********
	309
	310	c equivalent width of atmospheric line.
	311
	312	pi = acos(-1.d0)
	313
	314	if ( idummy.gt.9 )
	315	@ write (,) ' S, m, alsa, pp =', ka(kr), me, alsa(kr), plaux
	316
	317	y=ka(kr)*me
	318	! x=y/(2.0pi(alsa(kr)pl+alna(kr)(pe-pl)))
	319	x=y/(2.0d0pi alsa(kr)plaux) !+alna(kr)(pe-pl)))
	320
	321	! Strong limit of Lorentz profile: WL = 2 SQRT( S * m * alsa*pl )
	322	! Para anular esto, comentar las siguientes 5 lineas
	323	! if ( x .gt. 1.e6 ) then
	324	! wl = 2.0sqrt( y alsa(kr)*pl )
	325	! else
	326	! wl=y/sqrt(1.0d0+pi*x/2.0d0)
	327	! endif
	328
	329	wl=y/sqrt(1.0d0+pi*x/2.0d0)
	330
	331	if (wl .le. 0.d0) then
	332	write(,)'mztfsub_overlap/496',ig,y,ka(kr),me,kr
	333	stop'WE/Lorentz EQW zero or negative!/498' !,ig
	334	endif
	335
	336	if ( idummy.gt.9 )
	337	@ write (,) ' y, x =', y, x
	338
	339	xlor = x
	340	if ( (idummy.eq.2 .or. idummy.eq.12) .and. xlor.gt.1e5 ) then
	341	! en caso que estemos en el regimen
	342	! Strong Lorentz y la presion local
	343	! vaya disminuyendo, corregimos la EQW
	344	! con un gradiente analitico (notebook)
	345	wsL = 2.0sqrt( y alsa(kr)*plaux )
	346	wsl_grad = - 2.d0 * ka(kr)alsa(kr) nt_local*p_local / wsL
	347	wlgood = w_strongLor_prev(kr) + wsl_grad * Desp
	348	if (idummy.eq.12)
	349	@ write (,) ' W(wrong), W_SL, W_SL prev, W_SL corrected=',
	350	@ wl, wsL, w_strongLor_prev(kr), wlgood
	351	wl = wlgood
	352	endif
	353	! wsL = wl pero esto no lo hacemos todavia, porque necesitamos
	354	! el valor que ahora mismo tiene wsL para corregir la
	355	! expresion R&W below
	356
	357	! write (,) 'WE arguments me,pe,pl =', me,pe,pl
	358	! write (,) 'WE/ wl,ka(kr),alsa(kr) =',
	359	! @ wl, ka(kr),alsa(kr)
	360
	361
	362	!>>>>>>>
	363	500 format (a,i3,3(2x,1pe15.8))
	364	600 format (a,2(2x,1pe16.9))
	365	700 format (a,3(1x,1pe16.9))
	366	! if (kr.eq.8 .or. kr.eq.13) then
	367	! write (*,500) 'WE/kr,m,pt,pl=', kr, me, pe, pl
	368	! write (*,700) ' /aln,als,d_x=', alna(kr),alsa(kr),
	369	! @ 2.0pi( alsa(kr)pl + alna(kr)(pe-pl) )
	370	! write (,600) ' /alsap_CO2, alna*p_n2 :',
	371	! @ alsa(kr)pl, alna(kr)(pe-pl)
	372	! write (,600) ' ap, S =',
	373	! @ alsa(kr)pl + alna(kr)(pe-pl) , ka(kr)
	374	! write (,600) ' /Sm, x =', y, x
	375	! write (*,600) ' /aprox, WL =',
	376	! @ 2.sqrt( y( alsa(kr)pl+alna(kr)(pe-pl) ) ), WL
	377	! endif
	378	!
	379	! corrections to lorentz eqw due to overlaping and super-saturation
	380	!
	381
	382	i_supersat = 0
	383
	384	if ( icls.eq.5 .and. ioverlap.gt.0 ) then
	385	! for the moment, only consider overlaping for mztf.f, not fot.f
	386
	387	! definition of saturation in the lisat model
	388	!
	389	asat_box = 0.99d0
	390	f_sat_box = 2.d0 * x
	391	xx = f_sat_box / log( 1./(1-asat_box) )
	392	if ( xx .lt. 1.0d0 ) then
	393	dv_sat_box = 0.0d0
	394	asat_box = 1.0d0 - exp( - f_sat_box )
	395	else
	396	dv_sat_box = alsa(kr) * sqrt( xx - 1.0d0 )
	397	! approximation: only use of alsa in mars and venus
	398	endif
	399
	400	! area of saturated line
	401	!
	402	area_core_box = 2.d0 * dv_sat_box * asat_box
	403	area_wing_box = 0.5d0 * ( wl - area_core_box )
	404	dv_corte_box = dv_sat_box + 2.d0*area_wing_box/asat_box
	405
	406	! super-saturation or simple overlaping?
	407	!
	408	! i_supersat = 0
	409	xx = dv_sat_box - ( 0.5d0 * dist(kr) )
	410	if ( xx .ge. 0.0 ! definition of supersaturation
	411	@ .and. dv_sat_box .gt. 0.0 ! definition of saturation
	412	@ .and. (dist(kr).gt.0.0) ) ! box contains more than 1 line
	413	@ ! and not too far apart
	414	@ then
	415
	416	i_supersat = 1
	417
	418	else
	419	! no super-saturation, then use "lisat + first correction", i.e.,
	420	! correct for line products
	421	!
	422
	423	wl = wl
	424
	425	endif
	426
	427	end if ! end of overlaping loop
	428
	429	if (icls.eq.2) then
	430	we = wl
	431	return
	432	endif
	433
	434	cc doppler limit:
	435	if ( idummy.gt.9 )
	436	@ write (,) ' Sm, alf_dop =', y, alda(kr)sqrt(pi)
	437
	438	x = y / (alda(kr)*sqrt(pi))
	439	if ( x.lt.1.e-10 ) then ! to avoid underflow
	440	wd = y
	441	else
	442	wd=alda(kr)sqrt(4.0pix2(1.0+log(1.0+x))/(4.0+pix*2))
	443	endif
	444	if ( idummy.gt.9 )
	445	@ write (,) ' wd =', wd
	446
	447	cc doppler weak limit
	448	c wd = ka(kr) * me
	449
	450	cc good doppler
	451	if(icls.eq.5) then !para mztf
	452	!write (,) 'para mztf, icls=',icls
	453	if (x.lt.5.) then
	454	wd = 0.d0
	455	do i=0,7
	456	wd = wd + cn(i) * x**i
	457	end do
	458	wd = alda(kr) * x * sqrt(pi) * wd
	459	elseif (x.gt.5.) then
	460	wd = 0.d0
	461	do i=0,7
	462	wd = wd + dn(i) / (log(x))**i
	463	end do
	464	wd = alda(kr) * sqrt(log(x)) * wd
	465	else
	466	stop ' x should not be less than zero'
	467	end if
	468	end if
	469
	470
	471	if ( i_supersat .eq. 0 ) then
	472
	473	parentesis = wl2+wd2-(wdwl/y)*2
	474	! changed +(wdwl/y)*2 to -...14-3-84
	475
	476	if ( parentesis .lt. 0.0 ) then
	477	if ((idummy.eq.2 .or. idummy.eq.12) .and. xlor.gt.1e5) then
	478	parentesis = wl2+wd2-(wdwsL/y)*2
	479	! este cambio puede ser necesario cuando se hace
	480	! correccion Strong Lor, para evitar valores
	481	! negativos del parentesis en sqrt( )
	482	else
	483	stop ' WE/ Error en las EQW wl,wl,y '
	484	endif
	485	endif
	486
	487	we = sqrt( parentesis )
	488	! write (,) ' from we: xdop,alda,wd', sngl(x),alda(kr),sngl(wd)
	489	! write (,) ' from we: we', we
	490
	491	else
	492
	493	we = wl
	494	! if there is supersaturation we can ignore wd completely;
	495	! mztf.f will compute the eqw of the whole box afterwards
	496
	497	endif
	498
	499	if (icls.eq.3) we = wd
	500
	501	if ( idummy.gt.9 )
	502	@ write (,) ' wl,wd,w =', wl,wd,we
	503
	504	wsL = wl
	505
	506	return
	507	end
	508
	509
	510	c **********************************************************************
	511	real*8 function simrul(a,b,fsim,c,acc)
	512	c adaptively integrates fsim from a to b, within the criterion acc.
	513	c **********************************************************************
	514
	515	implicit none
	516
	517	real*8 res,a,b,g0,g1,g2,g3,g4,d,a0,a1,a2,h,x,acc,c,fsim
	518	real*8 s1(70),s2(70),s3(70)
	519	real*8 c1, c2
	520	integer*4 m,n,j
	521
	522	res=0.
	523	c=0.
	524	m=0
	525	n=0
	526	j=30
	527	g0=fsim(a)
	528	g2=fsim((a+b)/2.)
	529	g4=fsim(b)
	530	a0=(b-a)(g0+4.0g2+g4)/2.0
	531	1 d=2.0**n
	532	h=(b-a)/(4.0*d)
	533	x=a+(4.0m+1.0)h
	534	g1=fsim(x)
	535	g3=fsim(x+2.0*h)
	536	a1=h(g0+4.0g1+g2)
	537	a2=h(g2+4.0g3+g4)
	538	if ( abs(a1+a2-a0).gt.(acc/d)) goto 2
	539	res=res+(16.0*(a1+a2)-a0)/45.0
	540	m=m+1
	541	c=a+m*(b-a)/d
	542	6 if (m.eq.(2*(m/2))) goto 4
	543	if ((m.ne.1).or.(n.ne.0)) goto 5
	544	8 simrul=res
	545	return
	546	2 m=2*m
	547	n=n+1
	548	if (n.gt.j) goto 3
	549	a0=a1
	550	s1(n)=a2
	551	s2(n)=g3
	552	s3(n)=g4
	553	g4=g2
	554	g2=g1
	555	goto 1
	556	3 c1=c-(b-a)/d
	557	c2=c+(b-a)/d
	558	write(2,7) c1,c,c2,fsim(c1),fsim(c),fsim(c2)
	559	write(*,7) c1,c,c2,fsim(c1),fsim(c),fsim(c2)
	560	7 format(2x,17hsimrule fails at ,/,3e15.6,/,3e15.6)
	561	goto 8
	562	5 a0=s1(n)
	563	g0=g4
	564	g2=s2(n)
	565	g4=s3(n)
	566	goto 1
	567	4 m=m/2
	568	n=n-1
	569	goto 6
	570	end
	571
	572	c **********************************************************************
	573	subroutine findw(ig,iirw,idummy,c1,p1, Desp, wsL)
	574	c this routine sets up accuracy criteria and calls simrule between limit
	575	c that depend on the number of atmospheric and cell paths. it gives eqw.
	576
	577	c Add correction for EQW in Strong Lorentz regime and SZA>90
	578	c **********************************************************************
	579
	580	implicit none
	581	include 'nlte_paramdef.h'
	582	include 'nlte_commons.h'
	583
	584	c arguments
	585	integer ig ! ADDED FOR TRACEBACK
	586	integer iirw
	587	integer idummy ! I. indica varias opciones
	588	real*8 c1 ! I. needed for strong limit of Lorentz profil
	589	real*8 p1 ! I. " " "
	590	real*8 wsL ! O. need this for strong Lorentz correction
	591	real*8 Desp ! I. need this for strong Lorentz correction
	592
	593	c local variables
	594	real*8 ept,eps,xa
	595	real*8 acc, c
	596	real*8 we
	597	real*8 f, fi, simrul
	598
	599	external f,fi
	600
	601	c ******** ******* *******
	602
	603	if(icls.eq.5) then !para mztf
	604	! if(ig.eq.1682)write(,)'mztfsub_overlap/768',ua(kr),iirw
	605	if (iirw.eq.2) then !iirw=icf=2 ==> we use the w&r formula
	606	w = we(ig,ua(kr),pt,pp, idummy, c1,p1, Desp, wsL )
	607	return
	608	end if
	609	ept=we(ig,ua(kr),pt,pp, idummy,c1,p1, Desp, wsL)
	610	else !para fot
	611	if (iirw.eq.2) then ! icf=2 ==> we use the w&r formula
	612	w = we(ig,sl_ua,pt,pp, idummy,c1,p1, Desp, wsL)
	613	return
	614	end if
	615	ept=we(ig,sl_ua,pt,pp, idummy,c1,p1, Desp, wsL)
	616	end if
	617
	618	c the next block is a modification to avoid nul we.
	619	c this situation appears for weak lines and low path temperature, but
	620	c there is not any loss of accuracy. first july 1986
	621	if (ept.eq.0.) then ! for weak lines sometimes we=0
	622	ept=1.0e-18
	623	write (,) 'ept =',ept
	624	write (,) 'from we: we=0.0'
	625	return
	626	end if
	627
	628	acc = 4.d0
	629	acc = 10.d0**(-acc)
	630
	631	eps = acc * ept !accuracy 10-4 atmospheric eqw.
	632	xa=0.5*ept/f(0.d0) !width of doppler shifted atmospheric line.
	633	w=2.0*(simrul(0.0d0,xa,f,c,eps)+simrul(0.1d0,1.0/xa,fi,c,eps))
	634	!no shift.
	635
	636	return
	637	end
	638
	639
	640	c **********************************************************************
	641	double precision function fi(y)
	642	c returns the value of f(1/y)
	643	c **********************************************************************
	644
	645	implicit none
	646	real*8 f, y
	647
	648	fi=f(1.0/y)/y**2
	649	return
	650	end
	651
	652
	653	c **********************************************************************
	654	double precision function f(nuaux)
	655	c calculates 1-exp(-k(nu)u) for all series paths or combinations thereof
	656	c **********************************************************************
	657
	658	implicit none
	659	include 'nlte_paramdef.h'
	660	include 'nlte_commons.h'
	661
	662	double precision tra,xa,ya,za,yy,nuaux
	663	double precision voigtf
	664	tra=1.0d0
	665
	666	yy=1.0d0/alda(kr)
	667	xa=nuaux*yy
	668	ya= ( alsa(kr)pp + alna(kr)(pt-pp) ) * yy
	669	za=ka(kr)*yy
	670
	671	if(icls.eq.5) then !para mztf
	672	! write (,) 'icls=',icls
	673	tra=zaua(kr)voigtf(sngl(xa),sngl(ya))
	674	else
	675	tra=zasl_uavoigtf(sngl(xa),sngl(ya))
	676	end if
	677
	678	if (tra.gt.50.0) then
	679	tra=1.0 !2.0e-22 overflow cut-off.
	680	else if (tra.gt.1.0e-4) then
	681	tra=1.0-exp(-tra)
	682	end if
	683
	684	f=tra
	685	return
	686	end
	687
	688	c **********************************************************************
	689	double precision function voigtf(x1,y)
	690	c computes voigt function for any value of x1 and any +ve value of y.
	691	c where possible uses modified lorentz and modified doppler approximatio
	692	c otherwise uses a rearranged rybicki routine.
	693	c c(n) = exp(-(n/h)*2)/(pisqrt(pi)), with h = 2.5 .
	694	c accurate to better than 1 in 10000.
	695	c **********************************************************************
	696
	697	implicit none
	698
	699	real x1, y
	700	real x, xx, xxyy, xh,xhxh, yh,yhyh, f1,f2, p, q, xn,xnxn, voig
	701
	702	real*8 b,g0,g1,g2,g3,g4,d1,d2,d3,d4,c
	703	integer*4 n
	704
	705	dimension c(10)
	706	complex xp,xpp,z
	707
	708	data c(1)/0.15303405/
	709	data c(2)/0.94694928e-1/
	710	data c(3)/0.42549174e-1/
	711	data c(4)/0.13882935e-1/
	712	data c(5)/0.32892528e-2/
	713	data c(6)/0.56589906e-3/
	714	data c(7)/0.70697890e-4/
	715	data c(8)/0.64135678e-5/
	716	data c(9)/0.42249221e-6/
	717	data c(10)/0.20209868e-7/
	718
	719	x=abs(x1)
	720	if (x.gt.7.2) goto 1
	721	if ((y+x*0.3).gt.5.4) goto 1
	722	if (y.gt.0.01) goto 3
	723	if (x.lt.2.1) goto 2
	724	goto 3
	725	c here uses modified lorentz approx.
	726	1 xx=x*x
	727	xxyy=xx+y*y
	728	b=xx/xxyy
	729	voigtf=y(1.+(2.b-0.5+(0.75-(9.-12.b)b)/xxyy)/
	730	* xxyy)/(xxyy*3.141592654)
	731	return
	732	c here uses modified doppler approx.
	733	2 xx=x*x
	734	voigtf=0.56418958exp(-xx)(1.-y(1.-0.5y)(1.1289-xx(1.1623+
	735	* xx(0.080812+xx(0.13854-xx(0.033605-0.0073972xx))))))
	736	return
	737	c here uses a rearranged rybicki routine.
	738	3 xh=2.5*x
	739	xhxh=xh*xh
	740	yh=2.5*y
	741	yhyh=yh*yh
	742	f1=xhxh+yhyh
	743	f2=f1-0.5*yhyh
	744	if (y.lt.0.1) goto 20
	745	p=-y7.8539816 !7.8539816=2.5pi
	746	q=x*7.8539816
	747	xpp=cmplx(p,q)
	748	z=cexp(xpp)
	749	d1=xh*aimag(z)
	750	d2=-d1
	751	d3=yh*(1.-real(z))
	752	d4=-d3+2.*yh
	753	voig=0.17958712*(d1+d3)/f1
	754	goto 30
	755	20 p=x*7.8539816
	756	q=y*7.8539816
	757	xp=cmplx(p,q)
	758	z=ccos(xp)
	759	d1=xh*aimag(z)
	760	d2=-d1
	761	d3=yh*(1.-real(z))
	762	d4=-d3+2.*yh
	763	voig=0.56418958exp(yy-xx)cos(2.xy)+0.17958712*(d1+d3)/f1
	764	30 xn=0.
	765	do 55 n=1,10,2
	766	xn=xn+1.
	767	xnxn=xn*xn
	768	g1=xh-xn
	769	g2=g1*(xh+xn)
	770	g3=0.5g2g2
	771	voig=voig+c(n)(d2(g2+yhyh)+d4*(f1+xnxn))/
	772	& (g3+yhyh*(f2+xnxn))
	773	xn=xn+1.
	774	xnxn=xn*xn
	775	g1=xh-xn
	776	g2=g1*(xh+xn)
	777	g3=0.5g2g2
	778	voig=voig+c(n+1)(d1(g2+yhyh)+d3*(f1+xnxn))/
	779	@ (g3+yhyh*(f2+xnxn))
	780	55 continue
	781	voigtf=voig
	782	return
	783	end
	784
	785
	786
	787	c **********************************************************************
	788	c elimin_mz1d.F (includes smooth_cf)
	789	c ************************************************************************
	790	subroutine elimin_mz1d (c,vc, ilayer,nanaux,itblout, nwaux)
	791
	792	c Eliminate anomalous negative numbers in c(nl,nl) according to "nanaux":
	793
	794	c nanaux = 0 -> no eliminate
	795	c @ -> eliminate all numbers with absol.value<abs(max(c(n,r)))/300.
	796	c 2 -> eliminate all anomalous negative numbers in c(n,r).
	797	c 3 -> eliminate all anomalous negative numbers far from the main
	798	c diagonal.
	799	c 8 -> eliminate all non-zero numbers outside the main diagonal,
	800	c and the contibution of lower boundary.
	801	c 9 -> eliminate all non-zero numbers outside the main diagonal.
	802	c 4 -> hace un smoothing cuando la distancia de separacion entre
	803	c el valor maximo y el minimo de cf > 50 capas.
	804	c 5 -> elimina valores menores que 1.0d-19
	805	c 6 -> incluye los dos casos 4 y 5
	806	c 7 -> llama a lisa: smooth con width=nw & elimina mejorado
	807	c 78-> incluye los dos casos 7 y 8
	808	c 79-> incluye los dos casos 7 y 9
	809
	810	c itblout (itableout in calling program) is the option for writing
	811	c out or not the purged c(n,r) matrix:
	812	c itblout = 0 -> no write
	813	c = 1 -> write out in curtis***.out according to ilayer
	814
	815	c ilayer is the index for the layer selected to write out the matrix:
	816	c ilayer = 0 => matrix elements written out cover all the altitudes
	817	c with 5 layers steps
	818	c > 0 => " " " " are c(ilayer,*)
	819	c NOTA:
	820	c EXISTE LA POSIBILIDAD DE SACAR TODAS LAS CAPAS (TODA LA MATRIZ)
	821	c UTILIZANDO itableout=30 EN MZTUD
	822
	823	c jul 2011 malv+fgg adapted to LMD-MGCM
	824	c Sep-04 FGG+MALV Correct include and call parameters
	825	c cristina 25-sept-1996 y 27-ene-1997
	826	c JAN 98 MALV Version for mz1d
	827	c ************************************************************************
	828
	829	implicit none
	830
	831	include 'nlte_paramdef.h'
	832	include 'nlte_commons.h'
	833
	834	integer nanaux,j,i,itblout,kk,k,ir,in
	835	integer ilayer,jmin, jmax,np,nwaux,ntimes,ntimes2
	836	!* real*8 c(nl,nl), vc(nl), amax, cmax, cmin, cs(nl,nl), mini
	837	real*8 c(nl,nl), vc(nl), amax, cmax, cmin, mini
	838	real*8 aux(nl), auxs(nl)
	839	character layercode*3
	840
	841	ntimes=0
	842	ntimes2=0
	843	! type *,'from elimin_mz4: nan, nw',nan,nw
	844
	845	if (nanaux .eq. 0) goto 200
	846
	847	if(nanaux.eq.1)then
	848	do i=1,nl
	849	amax=1.0d-36
	850	do j=1,nl
	851	if(abs(c(i,j)).gt.amax)amax=abs(c(i,j))
	852	end do
	853	do j=1,nl
	854	if(abs(c(i,j)).lt.amax/300.0d0)c(i,j)=0.0d0
	855	end do
	856	enddo
	857	elseif(nanaux.eq.2)then
	858	do i=1,nl
	859	do j=1,nl
	860	if( ( j.le.(i-2) .or. j.gt.(i+2) ).and.
	861	@ ( c(i,j).lt.0.0d0 ) ) c(i,j)=0.0d0
	862	end do
	863	enddo
	864	elseif(nanaux.eq.3)then
	865	do i=1,nl
	866	do j=1,nl
	867	if (abs(i-j).ge.10) c(i,j)=0.0d0
	868	end do
	869	enddo
	870	elseif(nanaux.eq.8)then
	871	do i=1,nl
	872	do j=1,i-1
	873	c(i,j)=0.0d0
	874	enddo
	875	do j=i+1,nl
	876	c(i,j)=0.0d0
	877	enddo
	878	vc(i)= 0.d0
	879	enddo
	880	elseif(nanaux.eq.9)then
	881	do i=1,nl
	882	do j=1,i-1
	883	c(i,j)=0.0d0
	884	enddo
	885	do j=i+1,nl
	886	c(i,j)=0.0d0
	887	enddo
	888	enddo
	889	! elseif(nan.eq.7.or.nan.eq.78.or.nan.eq.79)then
	890	! call lisa(c, vc, nl, nw)
	891	end if
	892	if(nanaux.eq.78)then
	893	do i=1,nl
	894	do j=1,i-1
	895	c(i,j)=0.0d0
	896	enddo
	897	do j=i+1,nl
	898	c(i,j)=0.0d0
	899	enddo
	900	vc(i)= 0.d0
	901	enddo
	902	endif
	903	if(nanaux.eq.79)then
	904	do i=1,nl
	905	do j=1,i-1
	906	c(i,j)=0.0d0
	907	enddo
	908	do j=i+1,nl
	909	c(i,j)=0.0d0
	910	enddo
	911	enddo
	912	endif
	913
	914	if(nanaux.eq.5.or.nanaux.eq.6)then
	915	do i=1,nl
	916	mini = 1.0d-19
	917	do j=1,nl
	918	if(abs(c(i,j)).le.mini.and.c(i,j).ne.0.d0) then
	919	ntimes2=ntimes2+1
	920	end if
	921	if ( abs(c(i,j)).le.mini) c(i,j)=0.d0
	922	end do
	923	enddo
	924	end if
	925
	926	if(nanaux.eq.4.or.nanaux.eq.6)then
	927	do i=1,nl
	928	do j=1,nl
	929	aux(j)=c(i,j)
	930	auxs(j)=c(i,j)
	931	end do
	932	!call maxdp_2(aux,nl,cmax,jmax)
	933	cmax=maxval(aux)
	934	jmax=maxloc(aux,dim=1)
	935	if(abs(jmax-i).ge.50) then
	936	call smooth_cf(aux,auxs,i,nl,3)
	937	!!!call smooth_cf(aux,auxs,i,nl,5)
	938	ntimes=ntimes+1
	939	end if
	940	do j=1,nl
	941	c(i,j)=auxs(j)
	942	end do
	943	end do
	944	end if
	945
	946	! type *, 'elimin_mz4: c(n,r) procesed for elimination. '
	947	! type *, ' '
	948	! if(nan.eq.4.or.nan.eq.6) type *, ' call smoothing:',ntimes
	949	! if(nan.eq.5.or.nan.eq.6) type *, ' call elimina: ',ntimes2
	950	! if(nan.eq.7) type *, ' from elimin: lisa w=',nw
	951	! type *, ' '
	952
	953
	954	200 continue
	955
	956	c writting out of c(n,r) in ascii file
	957
	958	! if(itblout.eq.1) then
	959
	960	! if (ilayer.eq.0) then
	961
	962	! open (unit=2, status='new',
	963	! @ file=dircurtis//'curtis_gnu.out', recl=1024)
	964	! write(2,'(a)')
	965	! @ ' curtis matrix: table with 1.e+7 * acf(n,r) '
	966	! write(2,114) 'n,r', ( i, i=nl,1,-5 )
	967	! do in=nl,1,-5
	968	! write(2,*)
	969	! write(2,115) in, ( c(in,ir)*1.d7, ir=nl,1,-5 )
	970	! end do
	971	! close(2)
	972
	973
	974	! write (,) ' '
	975	! write (,) ' curtis.out has been created. '
	976	! write (,) ' '
	977
	978	! else
	979
	980	! write (layercode,132) ilayer
	981	! open (2, status='new',
	982	! @ file=dircurtis//'curtis'//layercode//'.out')
	983	! write(2,'(a)')
	984	! @ ' curtis matrix: table with 1.e+7 * acf(n,r) '
	985	! write(2,116) ' layer x c(',layercode,
	986	! @ ',x) c(x,', layercode,')'
	987	! do in=nl,1,-1
	988	! if (c(ilayer,ilayer).ne.0.d0) then
	989	! write(2,117) in, c(ilayer,in), c(in,ilayer),
	990	! @ c(ilayer,in)/c(ilayer,ilayer),
	991	! @ c(in,ilayer)/c(ilayer,ilayer)
	992	! else
	993	! write(2,118) in, c(ilayer,in), c(in,ilayer)
	994	! end if
	995	! end do
	996	! close(2)
	997	! write (,) ' '
	998	! write (,) dircurtis//'curtis'//layercode//'.out',
	999	! @ ' has been created.'
	1000	! write (,) ' '
	1001
	1002	! end if
	1003
	1004	! elseif(itblout.eq.0)then
	1005
	1006	! continue
	1007
	1008	! else
	1009
	1010	! write (,) ' error from elimin: ',
	1011	! @ ' itblout should be 1 or 0; itblout= ',itblout
	1012	! stop
	1013
	1014	! end if
	1015
	1016	return
	1017
	1018	112 format(10x,10(i3,9x))
	1019	113 format(1x,i3,2x,9(1pe9.2,2x))
	1020
	1021	114 format(1x,a3, 11(8x,i3))
	1022	115 format( 1x,i3, 2x, 11(1pe10.3))
	1023	116 format( 1x,a17,a2,a18,a2,a1 )
	1024	117 format( 3x,i3, 4(8x,1pe10.3) )
	1025	118 format( 3x,i3, 2(8x,1pe10.3) )
	1026	120 format( 1x,i3, 1x,i3, 2x, 11(1pe10.3))
	1027
	1028	132 format(i3)
	1029
	1030	! cambio: los formatos 114, 115 , 117 y 118
	1031	! cambio: al cambia nl de 51 a 140 hay que cambiar el formato i2-->i3
	1032	! y ahora en vez de 11 capas de 5 en 5, hay 28
	1033	!
	1034	end
	1035	c**************************************************************************
	1036	subroutine smooth_cf( c, cs, i, nl, w )
	1037	c hace un smoothing de c(i,*), de la contribucion de todas las capas
	1038	c menos de la capa en cuestion, la i.
	1039	c opcion w (width): el tamanho de la ventana del smoothing.
	1040	c output values: cs
	1041	c**************************************************************************
	1042
	1043	implicit none
	1044
	1045	integer j,np,i,nl,w
	1046	real*8 c(nl), cs(nl)
	1047
	1048	if(w.eq.0) then
	1049	do j=1,nl
	1050	cs(j)=c(j)
	1051	end do
	1052
	1053	elseif(w.eq.3) then
	1054
	1055	! write (,) 'smoothing w=3'
	1056	do j=1,i-4
	1057	if(j.eq.1) then
	1058	cs(j)=c(j)
	1059	else
	1060	cs(j)=1/3.d0*(c(j-1)+c(j)+c(j+1))
	1061	end if
	1062	end do
	1063	do j=i+4,nl-1
	1064	if(j.eq.nl) then
	1065	cs(j)=c(j)
	1066	else
	1067	cs(j)=1/3.d0*(c(j-1)+c(j)+c(j+1))
	1068	end if
	1069	end do
	1070	elseif(w.eq.5) then
	1071
	1072	! type *,'smoothing w=5'
	1073	do j=3,i-4
	1074	if(j.eq.1) then
	1075	cs(j)=c(j)
	1076	else
	1077	cs(j)=1/5.d0*(c(j-2)+c(j-1)+c(j)+c(j+1)+c(j+2))
	1078	end if
	1079	end do
	1080	do j=i+4,nl-2
	1081	if(j.eq.nl) then
	1082	cs(j)=c(j)
	1083	else
	1084	cs(j)=1/5.d0*(c(j-2)+c(j-1)+c(j)+c(j+1)+c(j+2))
	1085	end if
	1086	end do
	1087	end if
	1088	return
	1089	end
	1090
	1091
	1092
	1093	c*****************************************************************************
	1094	c suaviza
	1095	c*****************************************************************************
	1096	c
	1097	subroutine suaviza ( x, n, ismooth, y )
	1098	c
	1099	c x - input and return values
	1100	c y - auxiliary vector
	1101	c ismooth = 0 --> no smoothing is performed
	1102	c ismooth = 1 --> weak smoothing (5 points, centred weighted)
	1103	c ismooth = 2 --> normal smoothing (3 points, evenly weighted)
	1104	c ismooth = 3 --> strong smoothing (5 points, evenly weighted)
	1105
	1106
	1107	c malv august 1991
	1108	c*****************************************************************************
	1109
	1110	implicit none
	1111
	1112	integer n, imax, imin, i, ismooth
	1113	real*8 x(n), y(n)
	1114	c*****************************************************************************
	1115
	1116	imin=1
	1117	imax=n
	1118
	1119	if (ismooth.eq.0) then
	1120
	1121	return
	1122
	1123	elseif (ismooth.eq.1) then ! 5 points, with central weighting
	1124
	1125	do i=imin,imax
	1126	if(i.eq.imin)then
	1127	y(i)=x(imin)
	1128	elseif(i.eq.imax)then
	1129	y(i)=x(imax-1)+(x(imax-1)-x(imax-3))/2.d0
	1130	elseif(i.gt.(imin+1) .and. i.lt.(imax-1) )then
	1131	y(i) = ( x(i+2)/4.d0 + x(i+1)/2.d0 + 2.d0*x(i)/3.d0 +
	1132	& x(i-1)/2.d0 + x(i-2)/4.d0 )* 6.d0/13.d0
	1133	else
	1134	y(i)=(x(i+1)/2.d0+x(i)+x(i-1)/2.d0)/2.d0
	1135	end if
	1136	end do
	1137
	1138	elseif (ismooth.eq.2) then ! 3 points, evenly spaced
	1139
	1140	do i=imin,imax
	1141	if(i.eq.imin)then
	1142	y(i)=x(imin)
	1143	elseif(i.eq.imax)then
	1144	y(i)=x(imax-1)+(x(imax-1)-x(imax-3))/2.d0
	1145	else
	1146	y(i) = ( x(i+1)+x(i)+x(i-1) )/3.d0
	1147	end if
	1148	end do
	1149
	1150	elseif (ismooth.eq.3) then ! 5 points, evenly spaced
	1151
	1152	do i=imin,imax
	1153	if(i.eq.imin)then
	1154	y(i) = x(imin)
	1155	elseif(i.eq.(imin+1) .or. i.eq.(imax-1))then
	1156	y(i) = ( x(i+1)+x(i)+x(i-1) )/3.d0
	1157	elseif(i.eq.imax)then
	1158	y(i) = ( x(imax-1) + x(imax-1) + x(imax-2) ) / 3.d0
	1159	else
	1160	y(i) = ( x(i+2)+x(i+1)+x(i)+x(i-1)+x(i-2) )/5.d0
	1161	end if
	1162	end do
	1163
	1164	else
	1165
	1166	write (,) ' Error in suaviza.f Wrong ismooth value.'
	1167	stop
	1168
	1169	endif
	1170
	1171	c rehago el cambio, para devolver x(i)
	1172	do i=imin,imax
	1173	x(i)=y(i)
	1174	end do
	1175
	1176	return
	1177	end
	1178
	1179
	1180
	1181
	1182	c*****************************************************************************
	1183	c LUdec.F (includes lubksb_dp and ludcmp_dp subroutines)
	1184	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	1185	c
	1186	c Solution of linear equation without inverting matrix
	1187	c using LU decomposition:
	1188	c AA * xx = bb AA, bb: known
	1189	c xx: to be found
	1190	c AA and bb are not modified in this subroutine
	1191	c
	1192	c MALV , Sep 2007
	1193	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	1194
	1195	subroutine LUdec(xx,aa,bb,m,n)
	1196
	1197	implicit none
	1198
	1199	! Arguments
	1200	integer,intent(in) :: m, n
	1201	real*8,intent(in) :: aa(m,m), bb(m)
	1202	real*8,intent(out) :: xx(m)
	1203
	1204
	1205	! Local variables
	1206	real*8 a(n,n), b(n), x(n), d
	1207	integer i, j, indx(n)
	1208
	1209
	1210	! Subrutinas utilizadas
	1211	! ludcmp_dp, lubksb_dp
	1212
	1213	!!!!!!!!!!!!!!! Comienza el programa !!!!!!!!!!!!!!
	1214
	1215	do i=1,n
	1216	b(i) = bb(i+1)
	1217	do j=1,n
	1218	a(i,j) = aa(i+1,j+1)
	1219	enddo
	1220	enddo
	1221
	1222	! Descomposicion de auxm1
	1223	call ludcmp_dp ( a, n, n, indx, d)
	1224
	1225	! Sustituciones foward y backwards para hallar la solucion
	1226	do i=1,n
	1227	x(i) = b(i)
	1228	enddo
	1229	call lubksb_dp( a, n, n, indx, x )
	1230
	1231	do i=1,n
	1232	xx(i+1) = x(i)
	1233	enddo
	1234
	1235	return
	1236	end
	1237
	1238	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	1239
	1240	subroutine ludcmp_dp(a,n,np,indx,d)
	1241
	1242	c jul 2011 malv+fgg
	1243
	1244	implicit none
	1245
	1246	integer,intent(in) :: n, np
	1247	real*8,intent(inout) :: a(np,np)
	1248	real*8,intent(out) :: d
	1249	integer,intent(out) :: indx(n)
	1250
	1251	integer i, j, k, imax
	1252	real*8,parameter :: tiny=1.0d-20
	1253	real*8 vv(n), aamax, sum, dum
	1254
	1255
	1256	d=1.0d0
	1257	do 12 i=1,n
	1258	aamax=0.0d0
	1259	do 11 j=1,n
	1260	if (abs(a(i,j)).gt.aamax) aamax=abs(a(i,j))
	1261	11 continue
	1262	if (aamax.eq.0.0) then
	1263	write(,) 'ludcmp_dp: singular matrix!'
	1264	stop
	1265	endif
	1266	vv(i)=1.0d0/aamax
	1267	12 continue
	1268	do 19 j=1,n
	1269	if (j.gt.1) then
	1270	do 14 i=1,j-1
	1271	sum=a(i,j)
	1272	if (i.gt.1)then
	1273	do 13 k=1,i-1
	1274	sum=sum-a(i,k)*a(k,j)
	1275	13 continue
	1276	a(i,j)=sum
	1277	endif
	1278	14 continue
	1279	endif
	1280	aamax=0.0d0
	1281	do 16 i=j,n
	1282	sum=a(i,j)
	1283	if (j.gt.1)then
	1284	do 15 k=1,j-1
	1285	sum=sum-a(i,k)*a(k,j)
	1286	15 continue
	1287	a(i,j)=sum
	1288	endif
	1289	dum=vv(i)*abs(sum)
	1290	if (dum.ge.aamax) then
	1291	imax=i
	1292	aamax=dum
	1293	endif
	1294	16 continue
	1295	if (j.ne.imax)then
	1296	do 17 k=1,n
	1297	dum=a(imax,k)
	1298	a(imax,k)=a(j,k)
	1299	a(j,k)=dum
	1300	17 continue
	1301	d=-d
	1302	vv(imax)=vv(j)
	1303	endif
	1304	indx(j)=imax
	1305	if(j.ne.n)then
	1306	if(a(j,j).eq.0.0)a(j,j)=tiny
	1307	dum=1.0d0/a(j,j)
	1308	do 18 i=j+1,n
	1309	a(i,j)=a(i,j)*dum
	1310	18 continue
	1311	endif
	1312	19 continue
	1313	if(a(n,n).eq.0.0)a(n,n)=tiny
	1314	return
	1315	end
	1316
	1317	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	1318
	1319	subroutine lubksb_dp(a,n,np,indx,b)
	1320
	1321	c jul 2011 malv+fgg
	1322
	1323	implicit none
	1324
	1325	integer,intent(in) :: n,np
	1326	real*8,intent(in) :: a(np,np)
	1327	integer,intent(in) :: indx(n)
	1328	real*8,intent(out) :: b(n)
	1329
	1330	real*8 sum
	1331	integer ii, ll, i, j
	1332
	1333	ii=0
	1334	do 12 i=1,n
	1335	ll=indx(i)
	1336	sum=b(ll)
	1337	b(ll)=b(i)
	1338	if (ii.ne.0)then
	1339	do 11 j=ii,i-1
	1340	sum=sum-a(i,j)*b(j)
	1341	11 continue
	1342	else if (sum.ne.0.0) then
	1343	ii=i
	1344	endif
	1345	b(i)=sum
	1346	12 continue
	1347	do 14 i=n,1,-1
	1348	sum=b(i)
	1349	if(i.lt.n)then
	1350	do 13 j=i+1,n
	1351	sum=sum-a(i,j)*b(j)
	1352	13 continue
	1353	endif
	1354	b(i)=sum/a(i,i)
	1355	14 continue
	1356	return
	1357	end
	1358
	1359
	1360
	1361
	1362	c*****************************************************************************
	1363	c intersp
	1364	c ***********************************************************************
	1365	subroutine intersp(yy,zz,m,y,z,n,opt)
	1366	c interpolation soubroutine. input values: y(n) at z(n).
	1367	c output values: yy(m) at zz(m). options: 1 -> lineal; 2 -> logarithmic
	1368
	1369	c jul 2011 malv+fgg
	1370	c ***********************************************************************
	1371
	1372	implicit none
	1373
	1374	integer n,m,i,j,opt
	1375	real zz(m),yy(m),z(n),y(n)
	1376	real zmin,zzmin,zmax,zzmax
	1377
	1378	! write(,) ' interpolating'
	1379	! call minsp(z,n,zmin)
	1380	zmin=minval(z)
	1381	! call minsp(zz,m,zzmin)
	1382	zzmin=minval(zz)
	1383	! call maxsp(z,n,zmax)
	1384	zmax=maxval(z)
	1385	! call maxsp(zz,m,zzmax)
	1386	zzmax=maxval(zz)
	1387
	1388	if(zzmin.lt.zmin)then
	1389	write(,) 'from interp: new variable out of limits'
	1390	write(,) zzmin,'must be .ge. ',zmin
	1391	stop
	1392	! elseif(zzmax.gt.zmax)then
	1393	! write(,)'from interp: new variable out of limits'
	1394	! write(,)zzmax, 'must be .le. ',zmax
	1395	! stop
	1396	end if
	1397
	1398	do 1,i=1,m
	1399
	1400	do 2,j=1,n-1
	1401	if(zz(i).ge.z(j).and.zz(i).lt.z(j+1)) goto 3
	1402	2 continue
	1403	c in this case (zz(m).ge.z(n)) and j leaves the loop with j=n-1+1=n
	1404	if(opt.eq.1)then
	1405	yy(i)=y(n-1)+(y(n)-y(n-1))*(zz(i)-z(n-1))/(z(n)-z(n-1))
	1406	elseif(opt.eq.2)then
	1407	if(y(n).eq.0.0.or.y(n-1).eq.0.0)then
	1408	yy(i)=0.0
	1409	else
	1410	yy(i)=exp(log(y(n-1))+log(y(n)/y(n-1))*
	1411	@ (zz(i)-z(n-1))/(z(n)-z(n-1)))
	1412	end if
	1413	else
	1414	write(,)'from interp error: opt must be 1 or 2, opt= ',opt
	1415	end if
	1416	goto 1
	1417	3 continue
	1418	if(opt.eq.1)then
	1419	yy(i)=y(j)+(y(j+1)-y(j))*(zz(i)-z(j))/(z(j+1)-z(j))
	1420	elseif(opt.eq.2)then
	1421	if(y(j+1).eq.0.0.or.y(j).eq.0.0)then
	1422	yy(i)=0.0
	1423	else
	1424	yy(i)=exp(log(y(j))+log(y(j+1)/y(j))*
	1425	@ (zz(i)-z(j))/(z(j+1)-z(j)))
	1426	end if
	1427	else
	1428	write(,)'from interp error: opt must be 1 or 2, opt= ',opt
	1429	end if
	1430	1 continue
	1431
	1432	return
	1433	end
	1434
	1435
	1436
	1437	c*****************************************************************************
	1438	c interdp
	1439	c ***********************************************************************
	1440	subroutine interdp(yy,zz,m,y,z,n,opt)
	1441	c interpolation soubroutine. input values: y(n) at z(n).
	1442	c output values: yy(m) at zz(m). options: 1 -> lineal; 2 -> logarithmic
	1443	c jul 2011: malv+fgg Adapted to LMD-MGCM
	1444	c ***********************************************************************
	1445	implicit none
	1446	integer n,m,i,j,opt
	1447	real*8 zz(m),yy(m),z(n),y(n), zmin,zzmin,zmax,zzmax
	1448
	1449	! write (,) ' d interpolating '
	1450	! call mindp (z,n,zmin)
	1451	zmin=minval(z)
	1452	! call mindp (zz,m,zzmin)
	1453	zzmin=minval(zz)
	1454	! call maxdp (z,n,zmax)
	1455	zmax=maxval(z)
	1456	! call maxdp (zz,m,zzmax)
	1457	zzmax=maxval(zz)
	1458
	1459	if(zzmin.lt.zmin)then
	1460	write (,) 'from d interp: new variable out of limits'
	1461	write (,) zzmin,'must be .ge. ',zmin
	1462	stop
	1463	! elseif(zzmax.gt.zmax)then
	1464	! write (,) 'from interp: new variable out of limits'
	1465	! write (,) zzmax, 'must be .le. ',zmax
	1466	! stop
	1467	end if
	1468
	1469	do 1,i=1,m
	1470
	1471	do 2,j=1,n-1
	1472	if(zz(i).ge.z(j).and.zz(i).lt.z(j+1)) goto 3
	1473	2 continue
	1474	c in this case (zz(m).eq.z(n)) and j leaves the loop with j=n-1+1=n
	1475	if(opt.eq.1)then
	1476	yy(i)=y(n-1)+(y(n)-y(n-1))*(zz(i)-z(n-1))/(z(n)-z(n-1))
	1477	elseif(opt.eq.2)then
	1478	if(y(n).eq.0.0d0.or.y(n-1).eq.0.0d0)then
	1479	yy(i)=0.0d0
	1480	else
	1481	yy(i)=dexp(dlog(y(n-1))+dlog(y(n)/y(n-1))*
	1482	@ (zz(i)-z(n-1))/(z(n)-z(n-1)))
	1483	end if
	1484	else
	1485	write (,)
	1486	@ ' from d interp error: opt must be 1 or 2, opt= ',opt
	1487	end if
	1488	goto 1
	1489	3 continue
	1490	if(opt.eq.1)then
	1491	yy(i)=y(j)+(y(j+1)-y(j))*(zz(i)-z(j))/(z(j+1)-z(j))
	1492	! write (,) ' '
	1493	! write (,) ' z(j),z(j+1) =', z(j),z(j+1)
	1494	! write (,) ' t(j),t(j+1) =', y(j),y(j+1)
	1495	! write (,) ' zz, tt = ', zz(i), yy(i)
	1496	elseif(opt.eq.2)then
	1497	if(y(j+1).eq.0.0d0.or.y(j).eq.0.0d0)then
	1498	yy(i)=0.0d0
	1499	else
	1500	yy(i)=dexp(dlog(y(j))+dlog(y(j+1)/y(j))*
	1501	@ (zz(i)-z(j))/(z(j+1)-z(j)))
	1502	end if
	1503	else
	1504	write (,) ' from interp error: opt must be 1 or 2, opt= ',
	1505	@ opt
	1506	end if
	1507	1 continue
	1508	return
	1509	end
	1510
	1511
	1512	c*****************************************************************************
	1513	c interdp_limits.F
	1514	c ***********************************************************************
	1515
	1516	subroutine interdp_limits ( yy,zz,m, i1,i2, y,z,n, j1,j2, opt)
	1517
	1518	c Interpolation soubroutine.
	1519	c Returns values between indexes i1 & i2, donde 1 =< i1 =< i2 =< m
	1520	c Solo usan los indices de los inputs entre j1,j2, 1 =< j1 =< j2 =< n
	1521	c Input values: y(n) , z(n) (solo se usan los valores entre j1,j2)
	1522	c zz(m) (solo se necesita entre i1,i2)
	1523	c Output values: yy(m) (solo se calculan entre i1,i2)
	1524	c Options: opt=1 -> lineal ,, opt=2 -> logarithmic
	1525	c Difference with interdp:
	1526	c here interpolation proceeds between indexes i1,i2 only
	1527	c if i1=1 & i2=m, both subroutines are exactly the same
	1528	c thus previous calls to interdp or interdp2 could be easily replaced
	1529
	1530	c JAN 98 MALV Version for mz1d
	1531	c jul 2011 malv+fgg Adapted to LMD-MGCM
	1532	c ***********************************************************************
	1533
	1534	implicit none
	1535
	1536	! Arguments
	1537	integer n,m ! I. Dimensions
	1538	integer i1, i2, j1, j2, opt ! I
	1539	real*8 zz(m),yy(m) ! O
	1540	real*8 z(n),y(n) ! I
	1541
	1542	! Local variables
	1543	integer i,j
	1544	real*8 zmin,zzmin,zmax,zzmax
	1545
	1546	c *******************************
	1547
	1548	! type *, ' d interpolating '
	1549	! call mindp_limits (z,n,zmin, j1,j2)
	1550	zmin=minval(z(j1:j2))
	1551	! call mindp_limits (zz,m,zzmin, i1,i2)
	1552	zzmin=minval(zz(i1:i2))
	1553	! call maxdp_limits (z,n,zmax, j1,j2)
	1554	zmax=maxval(z(j1:j2))
	1555	! call maxdp_limits (zz,m,zzmax, i1,i2)
	1556	zzmax=maxval(zz(i1:i2))
	1557
	1558	if(zzmin.lt.zmin)then
	1559	write (,) 'from d interp: new variable out of limits'
	1560	write (,) zzmin,'must be .ge. ',zmin
	1561	stop
	1562	! elseif(zzmax.gt.zmax)then
	1563	! type *,'from interp: new variable out of limits'
	1564	! type *,zzmax, 'must be .le. ',zmax
	1565	! stop
	1566	end if
	1567
	1568	do 1,i=i1,i2
	1569
	1570	do 2,j=j1,j2-1
	1571	if(zz(i).ge.z(j).and.zz(i).lt.z(j+1)) goto 3
	1572	2 continue
	1573	c in this case (zz(i2).eq.z(j2)) and j leaves the loop with j=j2-1+1=j2
	1574	if(opt.eq.1)then
	1575	yy(i)=y(j2-1)+(y(j2)-y(j2-1))*
	1576	$ (zz(i)-z(j2-1))/(z(j2)-z(j2-1))
	1577	elseif(opt.eq.2)then
	1578	if(y(j2).eq.0.0d0.or.y(j2-1).eq.0.0d0)then
	1579	yy(i)=0.0d0
	1580	else
	1581	yy(i)=exp(log(y(j2-1))+log(y(j2)/y(j2-1))*
	1582	@ (zz(i)-z(j2-1))/(z(j2)-z(j2-1)))
	1583	end if
	1584	else
	1585	write (,) ' d interp : opt must be 1 or 2, opt= ',opt
	1586	end if
	1587	goto 1
	1588	3 continue
	1589	if(opt.eq.1)then
	1590	yy(i)=y(j)+(y(j+1)-y(j))*(zz(i)-z(j))/(z(j+1)-z(j))
	1591	! type *, ' '
	1592	! type *, ' z(j),z(j+1) =', z(j),z(j+1)
	1593	! type *, ' t(j),t(j+1) =', y(j),y(j+1)
	1594	! type *, ' zz, tt = ', zz(i), yy(i)
	1595	elseif(opt.eq.2)then
	1596	if(y(j+1).eq.0.0d0.or.y(j).eq.0.0d0)then
	1597	yy(i)=0.0d0
	1598	else
	1599	yy(i)=exp(log(y(j))+log(y(j+1)/y(j))*
	1600	@ (zz(i)-z(j))/(z(j+1)-z(j)))
	1601	end if
	1602	else
	1603	write (,) ' interp : opt must be 1 or 2, opt= ',opt
	1604	end if
	1605	1 continue
	1606	return
	1607	end
	1608
	1609
	1610
	1611
	1612	c*****************************************************************************
	1613	c Subroutines previously included in tcrco2_subrut.F
	1614	c***********************************************************************
	1615	c tcrco2_subrut.f
	1616	c
	1617	c jan 98 malv version for mz1d. copied from solar10/mz4sub.f
	1618	c jul 2011 malv+fgg adapted to LMD-MGCM
	1619	c***********************************************************************
	1620
	1621	************************************************************************
	1622
	1623	subroutine dinterconnection ( v, vt )
	1624
	1625	* input: vib. temp. from che*.for programs, vt(nl)
	1626	* output: test vibrational temp. for other che*.for, v(nl)
	1627	! iconex_smooth=1 ==> with smoothing
	1628	! iconex_smooth=0 ==> without smoothing
	1629	! iconex_tk=40 ==> with forced lte up to 40 km
	1630	! iconex_tk=20 ==> with forced lte up to 20 km
	1631	************************************************************************
	1632
	1633	implicit none
	1634	include 'nlte_paramdef.h'
	1635	include 'nlte_commons.h'
	1636
	1637	c argumentos
	1638	real*8 vt(nl), v(nl)
	1639
	1640	c local variables
	1641	integer i
	1642
	1643	c *************
	1644
	1645	do i=1,nl
	1646	v(i) = vt(i)
	1647	end do
	1648
	1649	! lo siguiente se utilizaba en solar10, pero es mejor introducirlo en
	1650	! la driver. por ahora no lo uso todavia.
	1651	! call fluctua(v,iconex_fluctua)
	1652	! call smooth_nl(v,iconex_smooth,nl)
	1653	! call forzar_tk(v,iconex_tk)
	1654
	1655	return
	1656	end
	1657
	1658	c***********************************************************************
	1659	function planckdp(tp,xnu)
	1660	c returns the black body function at wavenumber xnu and temperature t.
	1661	c***********************************************************************
	1662
	1663	implicit none
	1664
	1665	include 'nlte_paramdef.h'
	1666	include 'nlte_commons.h'
	1667
	1668	! common/datis/ pi, vlight, ee, hplanck, gamma, ab,
	1669	! @ n_avog, GG, R0, cte_sb, kboltzman, raddeg
	1670	! real*8 pi, vlight, ee, hplanck, gamma, ab,
	1671	! @ n_avog, GG, R0, cte_sb, kboltzman, raddeg
	1672
	1673	real*8 planckdp
	1674	real*8 xnu
	1675	real tp
	1676
	1677	planckdp = gammaxnu3.0 / exp( eexnu/dble(tp) )
	1678	!erg cm-2.sr-1/cm-1.
	1679
	1680	return
	1681	end
	1682
	1683	c ****************************************************************
	1684	function bandid (ib)
	1685	c returns the 2 character code of the band
	1686	c ****************************************************************
	1687	implicit none
	1688
	1689	integer ib
	1690	character*2 bandid
	1691
	1692	132 format(i2)
	1693	! encode (2,132,bandid) ib
	1694	write ( bandid, 132) ib
	1695
	1696	if ( ib .eq. 1 ) bandid = '01'
	1697	if ( ib .eq. 2 ) bandid = '02'
	1698	if ( ib .eq. 3 ) bandid = '03'
	1699	if ( ib .eq. 4 ) bandid = '04'
	1700	if ( ib .eq. 5 ) bandid = '05'
	1701	if ( ib .eq. 6 ) bandid = '06'
	1702	if ( ib .eq. 7 ) bandid = '07'
	1703	if ( ib .eq. 8 ) bandid = '08'
	1704	if ( ib .eq. 9 ) bandid = '09'
	1705	if ( ib .eq. 0 ) bandid = '00'
	1706
	1707	c end
	1708	return
	1709	end
	1710
	1711
	1712
	1713	c*****************************************************************************
	1714	c Subroutines previously included in mat_oper.F
	1715	c*****************************************************************************
	1716	c set of subroutines for the cz*.for programs:
	1717	! subroutine unit(a,n)
	1718	! subroutine diago(a,v,n) diagonal matrix with v
	1719	! subroutine invdiag(a,b,n) inverse of diagonal matrix
	1720	! subroutine sypvvv(a,b,c,d,n) suma y prod de 3 vectores, muy comun
	1721	! subroutine sypvmv(v,w,b,u,n) suma y prod de 3 vectores, muy comun
	1722	! subroutine mulmvv(w,b,u,v,n) prod matriz vector vector
	1723	! subroutine muymvv(w,b,u,v,n) prod matriz (inv.vector) vector
	1724	! subroutine samem (a,m,n)
	1725	! subroutine mulmv(a,b,c,n)
	1726	! subroutine mulmm(a,b,c,n)
	1727	! subroutine resmm(a,b,c,n)
	1728	! subroutine mulvv(a,b,c,n)
	1729	! subroutine sumvv(a,b,c,n)
	1730	! subroutine zerom(a,n)
	1731	! subroutine zero4m(a,b,c,d,n)
	1732	! subroutine zero3m(a,b,c,n)
	1733	! subroutine zero2m(a,b,n)
	1734	! subroutine zerov(a,n)
	1735	! subroutine zero4v(a,b,c,d,n)
	1736	! subroutine zero3v(a,b,c,n)
	1737	! subroutine zero2v(a,b,n)
	1738
	1739	!
	1740	!
	1741	! May-05 Sustituimos todos los zerojt de cristina por las subrutinas
	1742	! genericas zerov***
	1743	!
	1744	c ***********************************************************************
	1745	subroutine unit(a,n)
	1746	c store the unit value in the diagonal of a
	1747	c ***********************************************************************
	1748	real*8 a(n,n)
	1749	integer n,i,j,k
	1750	do 1,i=2,n-1
	1751	do 2,j=2,n-1
	1752	if(i.eq.j) then
	1753	a(i,j) = 1.d0
	1754	else
	1755	a(i,j)=0.0d0
	1756	end if
	1757	2 continue
	1758	1 continue
	1759	do k=1,n
	1760	a(n,k) = 0.0d0
	1761	a(1,k) = 0.0d0
	1762	a(k,1) = 0.0d0
	1763	a(k,n) = 0.0d0
	1764	end do
	1765	return
	1766	end
	1767
	1768	c ***********************************************************************
	1769	subroutine diago(a,v,n)
	1770	c store the vector v in the diagonal elements of the square matrix a
	1771	c ***********************************************************************
	1772	implicit none
	1773
	1774	integer n,i,j,k
	1775	real*8 a(n,n),v(n)
	1776
	1777	do 1,i=2,n-1
	1778	do 2,j=2,n-1
	1779	if(i.eq.j) then
	1780	a(i,j) = v(i)
	1781	else
	1782	a(i,j)=0.0d0
	1783	end if
	1784	2 continue
	1785	1 continue
	1786	do k=1,n
	1787	a(n,k) = 0.0d0
	1788	a(1,k) = 0.0d0
	1789	a(k,1) = 0.0d0
	1790	a(k,n) = 0.0d0
	1791	end do
	1792	return
	1793	end
	1794
	1795	c ***********************************************************************
	1796	subroutine samem (a,m,n)
	1797	c store the matrix m in the matrix a
	1798	c ***********************************************************************
	1799	real*8 a(n,n),m(n,n)
	1800	integer n,i,j,k
	1801	do 1,i=2,n-1
	1802	do 2,j=2,n-1
	1803	a(i,j) = m(i,j)
	1804	2 continue
	1805	1 continue
	1806	do k=1,n
	1807	a(n,k) = 0.0d0
	1808	a(1,k) = 0.0d0
	1809	a(k,1) = 0.0d0
	1810	a(k,n) = 0.0d0
	1811	end do
	1812	return
	1813	end
	1814	c ***********************************************************************
	1815	subroutine mulmv(a,b,c,n)
	1816	c do a(i)=b(i,j)*c(j). a, b, and c must be distint
	1817	c ***********************************************************************
	1818	implicit none
	1819
	1820	integer n,i,j
	1821	real*8 a(n),b(n,n),c(n),sum
	1822
	1823	do 1,i=2,n-1
	1824	sum=0.0d0
	1825	do 2,j=2,n-1
	1826	sum=sum+ (b(i,j)) * (c(j))
	1827	2 continue
	1828	a(i)=sum
	1829	1 continue
	1830	a(1) = 0.0d0
	1831	a(n) = 0.0d0
	1832	return
	1833	end
	1834
	1835	c ***********************************************************************
	1836	subroutine mulmm(a,b,c,n)
	1837	c ***********************************************************************
	1838	real*8 a(n,n), b(n,n), c(n,n)
	1839	integer n,i,j,k
	1840
	1841	! do i=2,n-1
	1842	! do j=2,n-1
	1843	! a(i,j)= 0.d00
	1844	! do k=2,n-1
	1845	! a(i,j) = a(i,j) + b(i,k) * c(k,j)
	1846	! end do
	1847	! end do
	1848	! end do
	1849	do j=2,n-1
	1850	do i=2,n-1
	1851	a(i,j)=0.d0
	1852	enddo
	1853	do k=2,n-1
	1854	do i=2,n-1
	1855	a(i,j)=a(i,j)+b(i,k)*c(k,j)
	1856	enddo
	1857	enddo
	1858	enddo
	1859	do k=1,n
	1860	a(n,k) = 0.0d0
	1861	a(1,k) = 0.0d0
	1862	a(k,1) = 0.0d0
	1863	a(k,n) = 0.0d0
	1864	end do
	1865
	1866	return
	1867	end
	1868
	1869	c ***********************************************************************
	1870	subroutine resmm(a,b,c,n)
	1871	c ***********************************************************************
	1872	real*8 a(n,n), b(n,n), c(n,n)
	1873	integer n,i,j,k
	1874
	1875	do i=2,n-1
	1876	do j=2,n-1
	1877	a(i,j)= b(i,j) - c(i,j)
	1878	end do
	1879	end do
	1880	do k=1,n
	1881	a(n,k) = 0.0d0
	1882	a(1,k) = 0.0d0
	1883	a(k,1) = 0.0d0
	1884	a(k,n) = 0.0d0
	1885	end do
	1886
	1887	return
	1888	end
	1889
	1890	c ***********************************************************************
	1891	subroutine sumvv(a,b,c,n)
	1892	c a(i)=b(i)+c(i)
	1893	c ***********************************************************************
	1894	implicit none
	1895
	1896	integer n,i
	1897	real*8 a(n),b(n),c(n)
	1898
	1899	do 1,i=2,n-1
	1900	a(i)= (b(i)) + (c(i))
	1901	1 continue
	1902	a(1) = 0.0d0
	1903	a(n) = 0.0d0
	1904	return
	1905	end
	1906
	1907	c ***********************************************************************
	1908	subroutine zerom(a,n)
	1909	c a(i,j)= 0.0
	1910	c ***********************************************************************
	1911
	1912	implicit none
	1913
	1914	integer n,i,j
	1915	real*8 a(n,n)
	1916
	1917	do 1,i=1,n
	1918	do 2,j=1,n
	1919	a(i,j) = 0.0d0
	1920	2 continue
	1921	1 continue
	1922	return
	1923	end
	1924
	1925	c ***********************************************************************
	1926	subroutine zero4m(a,b,c,d,n)
	1927	c a(i,j) = b(i,j) = c(i,j) = d(i,j) = 0.0
	1928	c ***********************************************************************
	1929	real*8 a(n,n), b(n,n), c(n,n), d(n,n)
	1930	integer n,i,j
	1931	do 1,i=1,n
	1932	do 2,j=1,n
	1933	a(i,j) = 0.0d0
	1934	b(i,j) = 0.0d0
	1935	c(i,j) = 0.0d0
	1936	d(i,j) = 0.0d0
	1937	2 continue
	1938	1 continue
	1939	return
	1940	end
	1941
	1942	c ***********************************************************************
	1943	subroutine zero3m(a,b,c,n)
	1944	c a(i,j) = b(i,j) = c(i,j) = 0.0
	1945	c **********************************************************************
	1946	real*8 a(n,n), b(n,n), c(n,n)
	1947	integer n,i,j
	1948	do 1,i=1,n
	1949	do 2,j=1,n
	1950	a(i,j) = 0.0d0
	1951	b(i,j) = 0.0d0
	1952	c(i,j) = 0.0d0
	1953	2 continue
	1954	1 continue
	1955	return
	1956	end
	1957
	1958	c ***********************************************************************
	1959	subroutine zero2m(a,b,n)
	1960	c a(i,j) = b(i,j) = 0.0
	1961	c ***********************************************************************
	1962	real*8 a(n,n), b(n,n)
	1963	integer n,i,j
	1964	do 1,i=1,n
	1965	do 2,j=1,n
	1966	a(i,j) = 0.0d0
	1967	b(i,j) = 0.0d0
	1968	2 continue
	1969	1 continue
	1970	return
	1971	end
	1972	c ***********************************************************************
	1973	subroutine zerov(a,n)
	1974	c a(i)= 0.0
	1975	c ***********************************************************************
	1976	real*8 a(n)
	1977	integer n,i
	1978	do 1,i=1,n
	1979	a(i) = 0.0d0
	1980	1 continue
	1981	return
	1982	end
	1983	c ***********************************************************************
	1984	subroutine zero4v(a,b,c,d,n)
	1985	c a(i) = b(i) = c(i) = d(i,j) = 0.0
	1986	c ***********************************************************************
	1987	real*8 a(n), b(n), c(n), d(n)
	1988	integer n,i
	1989	do 1,i=1,n
	1990	a(i) = 0.0d0
	1991	b(i) = 0.0d0
	1992	c(i) = 0.0d0
	1993	d(i) = 0.0d0
	1994	1 continue
	1995	return
	1996	end
	1997	c ***********************************************************************
	1998	subroutine zero3v(a,b,c,n)
	1999	c a(i) = b(i) = c(i) = 0.0
	2000	c ***********************************************************************
	2001	real*8 a(n), b(n), c(n)
	2002	integer n,i
	2003	do 1,i=1,n
	2004	a(i) = 0.0d0
	2005	b(i) = 0.0d0
	2006	c(i) = 0.0d0
	2007	1 continue
	2008	return
	2009	end
	2010	c ***********************************************************************
	2011	subroutine zero2v(a,b,n)
	2012	c a(i) = b(i) = 0.0
	2013	c ***********************************************************************
	2014	real*8 a(n), b(n)
	2015	integer n,i
	2016	do 1,i=1,n
	2017	a(i) = 0.0d0
	2018	b(i) = 0.0d0
	2019	1 continue
	2020	return
	2021	end
	2022
	2023

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

Download in other formats:

Original Format