Context Navigation

nlte_aux.F @ 694

Last change on this file since 694 was 694, checked in by aslmd, 12 years ago
LMDZ.MARS : forgot something in the previous commit.
File size: 78.6 KB

Line
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 iaa_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	iaa_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	iaa_we = sqrt( parentesis )
488	! write (,) ' from iaa_we: xdop,alda,wd', sngl(x),alda(kr),sngl(wd)
489	! write (,) ' from iaa_we: we', iaa_we
490
491	else
492
493	iaa_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) iaa_we = wd
500
501	if ( idummy.gt.9 )
502	@ write (,) ' wl,wd,w =', wl,wd,iaa_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 iaa_we
597	real*8 iaa_f, iaa_fi, simrul
598
599	external iaa_f,iaa_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 = iaa_we(ig,ua(kr),pt,pp, idummy, c1,p1, Desp, wsL )
607	return
608	end if
609	ept=iaa_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 = iaa_we(ig,sl_ua,pt,pp, idummy,c1,p1, Desp, wsL)
613	return
614	end if
615	ept=iaa_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/iaa_f(0.d0) !width of doppler shifted atmospheric line.
633	w = 2.0*( simrul(0.0d0,xa,iaa_f,c,eps)
634	. + simrul(0.1d0,1.0/xa,iaa_fi,c,eps) )
635	!no shift.
636
637	return
638	end
639
640
641	c **********************************************************************
642	double precision function iaa_fi(y)
643	c returns the value of f(1/y)
644	c **********************************************************************
645
646	implicit none
647	real*8 iaa_f, y
648
649	iaa_fi=iaa_f(1.0/y)/y**2
650	return
651	end
652
653
654	c **********************************************************************
655	double precision function iaa_f(nuaux)
656	c calculates 1-exp(-k(nu)u) for all series paths or combinations thereof
657	c **********************************************************************
658
659	implicit none
660	include 'nlte_paramdef.h'
661	include 'nlte_commons.h'
662
663	double precision tra,xa,ya,za,yy,nuaux
664	double precision voigtf
665	tra=1.0d0
666
667	yy=1.0d0/alda(kr)
668	xa=nuaux*yy
669	ya= ( alsa(kr)pp + alna(kr)(pt-pp) ) * yy
670	za=ka(kr)*yy
671
672	if(icls.eq.5) then !para mztf
673	! write (,) 'icls=',icls
674	tra=zaua(kr)voigtf(sngl(xa),sngl(ya))
675	else
676	tra=zasl_uavoigtf(sngl(xa),sngl(ya))
677	end if
678
679	if (tra.gt.50.0) then
680	tra=1.0 !2.0e-22 overflow cut-off.
681	else if (tra.gt.1.0e-4) then
682	tra=1.0-exp(-tra)
683	end if
684
685	iaa_f=tra
686	return
687	end
688
689	c **********************************************************************
690	double precision function voigtf(x1,y)
691	c computes voigt function for any value of x1 and any +ve value of y.
692	c where possible uses modified lorentz and modified doppler approximatio
693	c otherwise uses a rearranged rybicki routine.
694	c c(n) = exp(-(n/h)*2)/(pisqrt(pi)), with h = 2.5 .
695	c accurate to better than 1 in 10000.
696	c **********************************************************************
697
698	implicit none
699
700	real x1, y
701	real x, xx, xxyy, xh,xhxh, yh,yhyh, f1,f2, p, q, xn,xnxn, voig
702
703	real*8 b,g0,g1,g2,g3,g4,d1,d2,d3,d4,c
704	integer*4 n
705
706	dimension c(10)
707	complex xp,xpp,z
708
709	data c(1)/0.15303405/
710	data c(2)/0.94694928e-1/
711	data c(3)/0.42549174e-1/
712	data c(4)/0.13882935e-1/
713	data c(5)/0.32892528e-2/
714	data c(6)/0.56589906e-3/
715	data c(7)/0.70697890e-4/
716	data c(8)/0.64135678e-5/
717	data c(9)/0.42249221e-6/
718	data c(10)/0.20209868e-7/
719
720	x=abs(x1)
721	if (x.gt.7.2) goto 1
722	if ((y+x*0.3).gt.5.4) goto 1
723	if (y.gt.0.01) goto 3
724	if (x.lt.2.1) goto 2
725	goto 3
726	c here uses modified lorentz approx.
727	1 xx=x*x
728	xxyy=xx+y*y
729	b=xx/xxyy
730	voigtf=y(1.+(2.b-0.5+(0.75-(9.-12.b)b)/xxyy)/
731	* xxyy)/(xxyy*3.141592654)
732	return
733	c here uses modified doppler approx.
734	2 xx=x*x
735	voigtf=0.56418958exp(-xx)(1.-y(1.-0.5y)(1.1289-xx(1.1623+
736	* xx(0.080812+xx(0.13854-xx(0.033605-0.0073972xx))))))
737	return
738	c here uses a rearranged rybicki routine.
739	3 xh=2.5*x
740	xhxh=xh*xh
741	yh=2.5*y
742	yhyh=yh*yh
743	f1=xhxh+yhyh
744	f2=f1-0.5*yhyh
745	if (y.lt.0.1) goto 20
746	p=-y7.8539816 !7.8539816=2.5pi
747	q=x*7.8539816
748	xpp=cmplx(p,q)
749	z=cexp(xpp)
750	d1=xh*aimag(z)
751	d2=-d1
752	d3=yh*(1.-real(z))
753	d4=-d3+2.*yh
754	voig=0.17958712*(d1+d3)/f1
755	goto 30
756	20 p=x*7.8539816
757	q=y*7.8539816
758	xp=cmplx(p,q)
759	z=ccos(xp)
760	d1=xh*aimag(z)
761	d2=-d1
762	d3=yh*(1.-real(z))
763	d4=-d3+2.*yh
764	voig=0.56418958exp(yy-xx)cos(2.xy)+0.17958712*(d1+d3)/f1
765	30 xn=0.
766	do 55 n=1,10,2
767	xn=xn+1.
768	xnxn=xn*xn
769	g1=xh-xn
770	g2=g1*(xh+xn)
771	g3=0.5g2g2
772	voig=voig+c(n)(d2(g2+yhyh)+d4*(f1+xnxn))/
773	& (g3+yhyh*(f2+xnxn))
774	xn=xn+1.
775	xnxn=xn*xn
776	g1=xh-xn
777	g2=g1*(xh+xn)
778	g3=0.5g2g2
779	voig=voig+c(n+1)(d1(g2+yhyh)+d3*(f1+xnxn))/
780	@ (g3+yhyh*(f2+xnxn))
781	55 continue
782	voigtf=voig
783	return
784	end
785
786
787
788	c **********************************************************************
789	c elimin_mz1d.F (includes smooth_cf)
790	c ************************************************************************
791	subroutine elimin_mz1d (c,vc, ilayer,nanaux,itblout, nwaux)
792
793	c Eliminate anomalous negative numbers in c(nl,nl) according to "nanaux":
794
795	c nanaux = 0 -> no eliminate
796	c @ -> eliminate all numbers with absol.value<abs(max(c(n,r)))/300.
797	c 2 -> eliminate all anomalous negative numbers in c(n,r).
798	c 3 -> eliminate all anomalous negative numbers far from the main
799	c diagonal.
800	c 8 -> eliminate all non-zero numbers outside the main diagonal,
801	c and the contibution of lower boundary.
802	c 9 -> eliminate all non-zero numbers outside the main diagonal.
803	c 4 -> hace un smoothing cuando la distancia de separacion entre
804	c el valor maximo y el minimo de cf > 50 capas.
805	c 5 -> elimina valores menores que 1.0d-19
806	c 6 -> incluye los dos casos 4 y 5
807	c 7 -> llama a lisa: smooth con width=nw & elimina mejorado
808	c 78-> incluye los dos casos 7 y 8
809	c 79-> incluye los dos casos 7 y 9
810
811	c itblout (itableout in calling program) is the option for writing
812	c out or not the purged c(n,r) matrix:
813	c itblout = 0 -> no write
814	c = 1 -> write out in curtis***.out according to ilayer
815
816	c ilayer is the index for the layer selected to write out the matrix:
817	c ilayer = 0 => matrix elements written out cover all the altitudes
818	c with 5 layers steps
819	c > 0 => " " " " are c(ilayer,*)
820	c NOTA:
821	c EXISTE LA POSIBILIDAD DE SACAR TODAS LAS CAPAS (TODA LA MATRIZ)
822	c UTILIZANDO itableout=30 EN MZTUD
823
824	c jul 2011 malv+fgg adapted to LMD-MGCM
825	c Sep-04 FGG+MALV Correct include and call parameters
826	c cristina 25-sept-1996 y 27-ene-1997
827	c JAN 98 MALV Version for mz1d
828	c ************************************************************************
829
830	implicit none
831
832	include 'nlte_paramdef.h'
833	include 'nlte_commons.h'
834
835	integer nanaux,j,i,itblout,kk,k,ir,in
836	integer ilayer,jmin, jmax,np,nwaux,ntimes,ntimes2
837	!* real*8 c(nl,nl), vc(nl), amax, cmax, cmin, cs(nl,nl), mini
838	real*8 c(nl,nl), vc(nl), amax, cmax, cmin, mini
839	real*8 aux(nl), auxs(nl)
840	character layercode*3
841
842	ntimes=0
843	ntimes2=0
844	! type *,'from elimin_mz4: nan, nw',nan,nw
845
846	if (nanaux .eq. 0) goto 200
847
848	if(nanaux.eq.1)then
849	do i=1,nl
850	amax=1.0d-36
851	do j=1,nl
852	if(abs(c(i,j)).gt.amax)amax=abs(c(i,j))
853	end do
854	do j=1,nl
855	if(abs(c(i,j)).lt.amax/300.0d0)c(i,j)=0.0d0
856	end do
857	enddo
858	elseif(nanaux.eq.2)then
859	do i=1,nl
860	do j=1,nl
861	if( ( j.le.(i-2) .or. j.gt.(i+2) ).and.
862	@ ( c(i,j).lt.0.0d0 ) ) c(i,j)=0.0d0
863	end do
864	enddo
865	elseif(nanaux.eq.3)then
866	do i=1,nl
867	do j=1,nl
868	if (abs(i-j).ge.10) c(i,j)=0.0d0
869	end do
870	enddo
871	elseif(nanaux.eq.8)then
872	do i=1,nl
873	do j=1,i-1
874	c(i,j)=0.0d0
875	enddo
876	do j=i+1,nl
877	c(i,j)=0.0d0
878	enddo
879	vc(i)= 0.d0
880	enddo
881	elseif(nanaux.eq.9)then
882	do i=1,nl
883	do j=1,i-1
884	c(i,j)=0.0d0
885	enddo
886	do j=i+1,nl
887	c(i,j)=0.0d0
888	enddo
889	enddo
890	! elseif(nan.eq.7.or.nan.eq.78.or.nan.eq.79)then
891	! call lisa(c, vc, nl, nw)
892	end if
893	if(nanaux.eq.78)then
894	do i=1,nl
895	do j=1,i-1
896	c(i,j)=0.0d0
897	enddo
898	do j=i+1,nl
899	c(i,j)=0.0d0
900	enddo
901	vc(i)= 0.d0
902	enddo
903	endif
904	if(nanaux.eq.79)then
905	do i=1,nl
906	do j=1,i-1
907	c(i,j)=0.0d0
908	enddo
909	do j=i+1,nl
910	c(i,j)=0.0d0
911	enddo
912	enddo
913	endif
914
915	if(nanaux.eq.5.or.nanaux.eq.6)then
916	do i=1,nl
917	mini = 1.0d-19
918	do j=1,nl
919	if(abs(c(i,j)).le.mini.and.c(i,j).ne.0.d0) then
920	ntimes2=ntimes2+1
921	end if
922	if ( abs(c(i,j)).le.mini) c(i,j)=0.d0
923	end do
924	enddo
925	end if
926
927	if(nanaux.eq.4.or.nanaux.eq.6)then
928	do i=1,nl
929	do j=1,nl
930	aux(j)=c(i,j)
931	auxs(j)=c(i,j)
932	end do
933	!call maxdp_2(aux,nl,cmax,jmax)
934	cmax=maxval(aux)
935	jmax=maxloc(aux,dim=1)
936	if(abs(jmax-i).ge.50) then
937	call smooth_cf(aux,auxs,i,nl,3)
938	!!!call smooth_cf(aux,auxs,i,nl,5)
939	ntimes=ntimes+1
940	end if
941	do j=1,nl
942	c(i,j)=auxs(j)
943	end do
944	end do
945	end if
946
947	! type *, 'elimin_mz4: c(n,r) procesed for elimination. '
948	! type *, ' '
949	! if(nan.eq.4.or.nan.eq.6) type *, ' call smoothing:',ntimes
950	! if(nan.eq.5.or.nan.eq.6) type *, ' call elimina: ',ntimes2
951	! if(nan.eq.7) type *, ' from elimin: lisa w=',nw
952	! type *, ' '
953
954
955	200 continue
956
957	c writting out of c(n,r) in ascii file
958
959	! if(itblout.eq.1) then
960
961	! if (ilayer.eq.0) then
962
963	! open (unit=2, status='new',
964	! @ file=dircurtis//'curtis_gnu.out', recl=1024)
965	! write(2,'(a)')
966	! @ ' curtis matrix: table with 1.e+7 * acf(n,r) '
967	! write(2,114) 'n,r', ( i, i=nl,1,-5 )
968	! do in=nl,1,-5
969	! write(2,*)
970	! write(2,115) in, ( c(in,ir)*1.d7, ir=nl,1,-5 )
971	! end do
972	! close(2)
973
974
975	! write (,) ' '
976	! write (,) ' curtis.out has been created. '
977	! write (,) ' '
978
979	! else
980
981	! write (layercode,132) ilayer
982	! open (2, status='new',
983	! @ file=dircurtis//'curtis'//layercode//'.out')
984	! write(2,'(a)')
985	! @ ' curtis matrix: table with 1.e+7 * acf(n,r) '
986	! write(2,116) ' layer x c(',layercode,
987	! @ ',x) c(x,', layercode,')'
988	! do in=nl,1,-1
989	! if (c(ilayer,ilayer).ne.0.d0) then
990	! write(2,117) in, c(ilayer,in), c(in,ilayer),
991	! @ c(ilayer,in)/c(ilayer,ilayer),
992	! @ c(in,ilayer)/c(ilayer,ilayer)
993	! else
994	! write(2,118) in, c(ilayer,in), c(in,ilayer)
995	! end if
996	! end do
997	! close(2)
998	! write (,) ' '
999	! write (,) dircurtis//'curtis'//layercode//'.out',
1000	! @ ' has been created.'
1001	! write (,) ' '
1002
1003	! end if
1004
1005	! elseif(itblout.eq.0)then
1006
1007	! continue
1008
1009	! else
1010
1011	! write (,) ' error from elimin: ',
1012	! @ ' itblout should be 1 or 0; itblout= ',itblout
1013	! stop
1014
1015	! end if
1016
1017	return
1018
1019	112 format(10x,10(i3,9x))
1020	113 format(1x,i3,2x,9(1pe9.2,2x))
1021
1022	114 format(1x,a3, 11(8x,i3))
1023	115 format( 1x,i3, 2x, 11(1pe10.3))
1024	116 format( 1x,a17,a2,a18,a2,a1 )
1025	117 format( 3x,i3, 4(8x,1pe10.3) )
1026	118 format( 3x,i3, 2(8x,1pe10.3) )
1027	120 format( 1x,i3, 1x,i3, 2x, 11(1pe10.3))
1028
1029	132 format(i3)
1030
1031	! cambio: los formatos 114, 115 , 117 y 118
1032	! cambio: al cambia nl de 51 a 140 hay que cambiar el formato i2-->i3
1033	! y ahora en vez de 11 capas de 5 en 5, hay 28
1034	!
1035	end
1036	c**************************************************************************
1037	subroutine smooth_cf( c, cs, i, nl, w )
1038	c hace un smoothing de c(i,*), de la contribucion de todas las capas
1039	c menos de la capa en cuestion, la i.
1040	c opcion w (width): el tamanho de la ventana del smoothing.
1041	c output values: cs
1042	c**************************************************************************
1043
1044	implicit none
1045
1046	integer j,np,i,nl,w
1047	real*8 c(nl), cs(nl)
1048
1049	if(w.eq.0) then
1050	do j=1,nl
1051	cs(j)=c(j)
1052	end do
1053
1054	elseif(w.eq.3) then
1055
1056	! write (,) 'smoothing w=3'
1057	do j=1,i-4
1058	if(j.eq.1) then
1059	cs(j)=c(j)
1060	else
1061	cs(j)=1/3.d0*(c(j-1)+c(j)+c(j+1))
1062	end if
1063	end do
1064	do j=i+4,nl-1
1065	if(j.eq.nl) then
1066	cs(j)=c(j)
1067	else
1068	cs(j)=1/3.d0*(c(j-1)+c(j)+c(j+1))
1069	end if
1070	end do
1071	elseif(w.eq.5) then
1072
1073	! type *,'smoothing w=5'
1074	do j=3,i-4
1075	if(j.eq.1) then
1076	cs(j)=c(j)
1077	else
1078	cs(j)=1/5.d0*(c(j-2)+c(j-1)+c(j)+c(j+1)+c(j+2))
1079	end if
1080	end do
1081	do j=i+4,nl-2
1082	if(j.eq.nl) then
1083	cs(j)=c(j)
1084	else
1085	cs(j)=1/5.d0*(c(j-2)+c(j-1)+c(j)+c(j+1)+c(j+2))
1086	end if
1087	end do
1088	end if
1089	return
1090	end
1091
1092
1093
1094	c*****************************************************************************
1095	c suaviza
1096	c*****************************************************************************
1097	c
1098	subroutine suaviza ( x, n, ismooth, y )
1099	c
1100	c x - input and return values
1101	c y - auxiliary vector
1102	c ismooth = 0 --> no smoothing is performed
1103	c ismooth = 1 --> weak smoothing (5 points, centred weighted)
1104	c ismooth = 2 --> normal smoothing (3 points, evenly weighted)
1105	c ismooth = 3 --> strong smoothing (5 points, evenly weighted)
1106
1107
1108	c malv august 1991
1109	c*****************************************************************************
1110
1111	implicit none
1112
1113	integer n, imax, imin, i, ismooth
1114	real*8 x(n), y(n)
1115	c*****************************************************************************
1116
1117	imin=1
1118	imax=n
1119
1120	if (ismooth.eq.0) then
1121
1122	return
1123
1124	elseif (ismooth.eq.1) then ! 5 points, with central weighting
1125
1126	do i=imin,imax
1127	if(i.eq.imin)then
1128	y(i)=x(imin)
1129	elseif(i.eq.imax)then
1130	y(i)=x(imax-1)+(x(imax-1)-x(imax-3))/2.d0
1131	elseif(i.gt.(imin+1) .and. i.lt.(imax-1) )then
1132	y(i) = ( x(i+2)/4.d0 + x(i+1)/2.d0 + 2.d0*x(i)/3.d0 +
1133	& x(i-1)/2.d0 + x(i-2)/4.d0 )* 6.d0/13.d0
1134	else
1135	y(i)=(x(i+1)/2.d0+x(i)+x(i-1)/2.d0)/2.d0
1136	end if
1137	end do
1138
1139	elseif (ismooth.eq.2) then ! 3 points, evenly spaced
1140
1141	do i=imin,imax
1142	if(i.eq.imin)then
1143	y(i)=x(imin)
1144	elseif(i.eq.imax)then
1145	y(i)=x(imax-1)+(x(imax-1)-x(imax-3))/2.d0
1146	else
1147	y(i) = ( x(i+1)+x(i)+x(i-1) )/3.d0
1148	end if
1149	end do
1150
1151	elseif (ismooth.eq.3) then ! 5 points, evenly spaced
1152
1153	do i=imin,imax
1154	if(i.eq.imin)then
1155	y(i) = x(imin)
1156	elseif(i.eq.(imin+1) .or. i.eq.(imax-1))then
1157	y(i) = ( x(i+1)+x(i)+x(i-1) )/3.d0
1158	elseif(i.eq.imax)then
1159	y(i) = ( x(imax-1) + x(imax-1) + x(imax-2) ) / 3.d0
1160	else
1161	y(i) = ( x(i+2)+x(i+1)+x(i)+x(i-1)+x(i-2) )/5.d0
1162	end if
1163	end do
1164
1165	else
1166
1167	write (,) ' Error in suaviza.f Wrong ismooth value.'
1168	stop
1169
1170	endif
1171
1172	c rehago el cambio, para devolver x(i)
1173	do i=imin,imax
1174	x(i)=y(i)
1175	end do
1176
1177	return
1178	end
1179
1180
1181
1182
1183	c*****************************************************************************
1184	c LUdec.F (includes lubksb_dp and ludcmp_dp subroutines)
1185	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1186	c
1187	c Solution of linear equation without inverting matrix
1188	c using LU decomposition:
1189	c AA * xx = bb AA, bb: known
1190	c xx: to be found
1191	c AA and bb are not modified in this subroutine
1192	c
1193	c MALV , Sep 2007
1194	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1195
1196	subroutine LUdec(xx,aa,bb,m,n)
1197
1198	implicit none
1199
1200	! Arguments
1201	integer,intent(in) :: m, n
1202	real*8,intent(in) :: aa(m,m), bb(m)
1203	real*8,intent(out) :: xx(m)
1204
1205
1206	! Local variables
1207	real*8 a(n,n), b(n), x(n), d
1208	integer i, j, indx(n)
1209
1210
1211	! Subrutinas utilizadas
1212	! ludcmp_dp, lubksb_dp
1213
1214	!!!!!!!!!!!!!!! Comienza el programa !!!!!!!!!!!!!!
1215
1216	do i=1,n
1217	b(i) = bb(i+1)
1218	do j=1,n
1219	a(i,j) = aa(i+1,j+1)
1220	enddo
1221	enddo
1222
1223	! Descomposicion de auxm1
1224	call ludcmp_dp ( a, n, n, indx, d)
1225
1226	! Sustituciones foward y backwards para hallar la solucion
1227	do i=1,n
1228	x(i) = b(i)
1229	enddo
1230	call lubksb_dp( a, n, n, indx, x )
1231
1232	do i=1,n
1233	xx(i+1) = x(i)
1234	enddo
1235
1236	return
1237	end
1238
1239	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1240
1241	subroutine ludcmp_dp(a,n,np,indx,d)
1242
1243	c jul 2011 malv+fgg
1244
1245	implicit none
1246
1247	integer,intent(in) :: n, np
1248	real*8,intent(inout) :: a(np,np)
1249	real*8,intent(out) :: d
1250	integer,intent(out) :: indx(n)
1251
1252	integer i, j, k, imax
1253	real*8,parameter :: tiny=1.0d-20
1254	real*8 vv(n), aamax, sum, dum
1255
1256
1257	d=1.0d0
1258	do 12 i=1,n
1259	aamax=0.0d0
1260	do 11 j=1,n
1261	if (abs(a(i,j)).gt.aamax) aamax=abs(a(i,j))
1262	11 continue
1263	if (aamax.eq.0.0) then
1264	write(,) 'ludcmp_dp: singular matrix!'
1265	stop
1266	endif
1267	vv(i)=1.0d0/aamax
1268	12 continue
1269	do 19 j=1,n
1270	if (j.gt.1) then
1271	do 14 i=1,j-1
1272	sum=a(i,j)
1273	if (i.gt.1)then
1274	do 13 k=1,i-1
1275	sum=sum-a(i,k)*a(k,j)
1276	13 continue
1277	a(i,j)=sum
1278	endif
1279	14 continue
1280	endif
1281	aamax=0.0d0
1282	do 16 i=j,n
1283	sum=a(i,j)
1284	if (j.gt.1)then
1285	do 15 k=1,j-1
1286	sum=sum-a(i,k)*a(k,j)
1287	15 continue
1288	a(i,j)=sum
1289	endif
1290	dum=vv(i)*abs(sum)
1291	if (dum.ge.aamax) then
1292	imax=i
1293	aamax=dum
1294	endif
1295	16 continue
1296	if (j.ne.imax)then
1297	do 17 k=1,n
1298	dum=a(imax,k)
1299	a(imax,k)=a(j,k)
1300	a(j,k)=dum
1301	17 continue
1302	d=-d
1303	vv(imax)=vv(j)
1304	endif
1305	indx(j)=imax
1306	if(j.ne.n)then
1307	if(a(j,j).eq.0.0)a(j,j)=tiny
1308	dum=1.0d0/a(j,j)
1309	do 18 i=j+1,n
1310	a(i,j)=a(i,j)*dum
1311	18 continue
1312	endif
1313	19 continue
1314	if(a(n,n).eq.0.0)a(n,n)=tiny
1315	return
1316	end
1317
1318	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1319
1320	subroutine lubksb_dp(a,n,np,indx,b)
1321
1322	c jul 2011 malv+fgg
1323
1324	implicit none
1325
1326	integer,intent(in) :: n,np
1327	real*8,intent(in) :: a(np,np)
1328	integer,intent(in) :: indx(n)
1329	real*8,intent(out) :: b(n)
1330
1331	real*8 sum
1332	integer ii, ll, i, j
1333
1334	ii=0
1335	do 12 i=1,n
1336	ll=indx(i)
1337	sum=b(ll)
1338	b(ll)=b(i)
1339	if (ii.ne.0)then
1340	do 11 j=ii,i-1
1341	sum=sum-a(i,j)*b(j)
1342	11 continue
1343	else if (sum.ne.0.0) then
1344	ii=i
1345	endif
1346	b(i)=sum
1347	12 continue
1348	do 14 i=n,1,-1
1349	sum=b(i)
1350	if(i.lt.n)then
1351	do 13 j=i+1,n
1352	sum=sum-a(i,j)*b(j)
1353	13 continue
1354	endif
1355	b(i)=sum/a(i,i)
1356	14 continue
1357	return
1358	end
1359
1360
1361
1362
1363	c*****************************************************************************
1364	c intersp
1365	c ***********************************************************************
1366	subroutine intersp(yy,zz,m,y,z,n,opt)
1367	c interpolation soubroutine. input values: y(n) at z(n).
1368	c output values: yy(m) at zz(m). options: 1 -> lineal; 2 -> logarithmic
1369
1370	c jul 2011 malv+fgg
1371	c ***********************************************************************
1372
1373	implicit none
1374
1375	integer n,m,i,j,opt
1376	real zz(m),yy(m),z(n),y(n)
1377	real zmin,zzmin,zmax,zzmax
1378
1379	! write(,) ' interpolating'
1380	! call minsp(z,n,zmin)
1381	zmin=minval(z)
1382	! call minsp(zz,m,zzmin)
1383	zzmin=minval(zz)
1384	! call maxsp(z,n,zmax)
1385	zmax=maxval(z)
1386	! call maxsp(zz,m,zzmax)
1387	zzmax=maxval(zz)
1388
1389	if(zzmin.lt.zmin)then
1390	write(,) 'from interp: new variable out of limits'
1391	write(,) zzmin,'must be .ge. ',zmin
1392	stop
1393	! elseif(zzmax.gt.zmax)then
1394	! write(,)'from interp: new variable out of limits'
1395	! write(,)zzmax, 'must be .le. ',zmax
1396	! stop
1397	end if
1398
1399	do 1,i=1,m
1400
1401	do 2,j=1,n-1
1402	if(zz(i).ge.z(j).and.zz(i).lt.z(j+1)) goto 3
1403	2 continue
1404	c in this case (zz(m).ge.z(n)) and j leaves the loop with j=n-1+1=n
1405	if(opt.eq.1)then
1406	yy(i)=y(n-1)+(y(n)-y(n-1))*(zz(i)-z(n-1))/(z(n)-z(n-1))
1407	elseif(opt.eq.2)then
1408	if(y(n).eq.0.0.or.y(n-1).eq.0.0)then
1409	yy(i)=0.0
1410	else
1411	yy(i)=exp(log(y(n-1))+log(y(n)/y(n-1))*
1412	@ (zz(i)-z(n-1))/(z(n)-z(n-1)))
1413	end if
1414	else
1415	write(,)'from interp error: opt must be 1 or 2, opt= ',opt
1416	end if
1417	goto 1
1418	3 continue
1419	if(opt.eq.1)then
1420	yy(i)=y(j)+(y(j+1)-y(j))*(zz(i)-z(j))/(z(j+1)-z(j))
1421	elseif(opt.eq.2)then
1422	if(y(j+1).eq.0.0.or.y(j).eq.0.0)then
1423	yy(i)=0.0
1424	else
1425	yy(i)=exp(log(y(j))+log(y(j+1)/y(j))*
1426	@ (zz(i)-z(j))/(z(j+1)-z(j)))
1427	end if
1428	else
1429	write(,)'from interp error: opt must be 1 or 2, opt= ',opt
1430	end if
1431	1 continue
1432
1433	return
1434	end
1435
1436
1437
1438	c*****************************************************************************
1439	c interdp
1440	c ***********************************************************************
1441	subroutine interdp(yy,zz,m,y,z,n,opt)
1442	c interpolation soubroutine. input values: y(n) at z(n).
1443	c output values: yy(m) at zz(m). options: 1 -> lineal; 2 -> logarithmic
1444	c jul 2011: malv+fgg Adapted to LMD-MGCM
1445	c ***********************************************************************
1446	implicit none
1447	integer n,m,i,j,opt
1448	real*8 zz(m),yy(m),z(n),y(n), zmin,zzmin,zmax,zzmax
1449
1450	! write (,) ' d interpolating '
1451	! call mindp (z,n,zmin)
1452	zmin=minval(z)
1453	! call mindp (zz,m,zzmin)
1454	zzmin=minval(zz)
1455	! call maxdp (z,n,zmax)
1456	zmax=maxval(z)
1457	! call maxdp (zz,m,zzmax)
1458	zzmax=maxval(zz)
1459
1460	if(zzmin.lt.zmin)then
1461	write (,) 'from d interp: new variable out of limits'
1462	write (,) zzmin,'must be .ge. ',zmin
1463	stop
1464	! elseif(zzmax.gt.zmax)then
1465	! write (,) 'from interp: new variable out of limits'
1466	! write (,) zzmax, 'must be .le. ',zmax
1467	! stop
1468	end if
1469
1470	do 1,i=1,m
1471
1472	do 2,j=1,n-1
1473	if(zz(i).ge.z(j).and.zz(i).lt.z(j+1)) goto 3
1474	2 continue
1475	c in this case (zz(m).eq.z(n)) and j leaves the loop with j=n-1+1=n
1476	if(opt.eq.1)then
1477	yy(i)=y(n-1)+(y(n)-y(n-1))*(zz(i)-z(n-1))/(z(n)-z(n-1))
1478	elseif(opt.eq.2)then
1479	if(y(n).eq.0.0d0.or.y(n-1).eq.0.0d0)then
1480	yy(i)=0.0d0
1481	else
1482	yy(i)=dexp(dlog(y(n-1))+dlog(y(n)/y(n-1))*
1483	@ (zz(i)-z(n-1))/(z(n)-z(n-1)))
1484	end if
1485	else
1486	write (,)
1487	@ ' from d interp error: opt must be 1 or 2, opt= ',opt
1488	end if
1489	goto 1
1490	3 continue
1491	if(opt.eq.1)then
1492	yy(i)=y(j)+(y(j+1)-y(j))*(zz(i)-z(j))/(z(j+1)-z(j))
1493	! write (,) ' '
1494	! write (,) ' z(j),z(j+1) =', z(j),z(j+1)
1495	! write (,) ' t(j),t(j+1) =', y(j),y(j+1)
1496	! write (,) ' zz, tt = ', zz(i), yy(i)
1497	elseif(opt.eq.2)then
1498	if(y(j+1).eq.0.0d0.or.y(j).eq.0.0d0)then
1499	yy(i)=0.0d0
1500	else
1501	yy(i)=dexp(dlog(y(j))+dlog(y(j+1)/y(j))*
1502	@ (zz(i)-z(j))/(z(j+1)-z(j)))
1503	end if
1504	else
1505	write (,) ' from interp error: opt must be 1 or 2, opt= ',
1506	@ opt
1507	end if
1508	1 continue
1509	return
1510	end
1511
1512
1513	c*****************************************************************************
1514	c interdp_limits.F
1515	c ***********************************************************************
1516
1517	subroutine interdp_limits ( yy,zz,m, i1,i2, y,z,n, j1,j2, opt)
1518
1519	c Interpolation soubroutine.
1520	c Returns values between indexes i1 & i2, donde 1 =< i1 =< i2 =< m
1521	c Solo usan los indices de los inputs entre j1,j2, 1 =< j1 =< j2 =< n
1522	c Input values: y(n) , z(n) (solo se usan los valores entre j1,j2)
1523	c zz(m) (solo se necesita entre i1,i2)
1524	c Output values: yy(m) (solo se calculan entre i1,i2)
1525	c Options: opt=1 -> lineal ,, opt=2 -> logarithmic
1526	c Difference with interdp:
1527	c here interpolation proceeds between indexes i1,i2 only
1528	c if i1=1 & i2=m, both subroutines are exactly the same
1529	c thus previous calls to interdp or interdp2 could be easily replaced
1530
1531	c JAN 98 MALV Version for mz1d
1532	c jul 2011 malv+fgg Adapted to LMD-MGCM
1533	c ***********************************************************************
1534
1535	implicit none
1536
1537	! Arguments
1538	integer n,m ! I. Dimensions
1539	integer i1, i2, j1, j2, opt ! I
1540	real*8 zz(m),yy(m) ! O
1541	real*8 z(n),y(n) ! I
1542
1543	! Local variables
1544	integer i,j
1545	real*8 zmin,zzmin,zmax,zzmax
1546
1547	c *******************************
1548
1549	! type *, ' d interpolating '
1550	! call mindp_limits (z,n,zmin, j1,j2)
1551	zmin=minval(z(j1:j2))
1552	! call mindp_limits (zz,m,zzmin, i1,i2)
1553	zzmin=minval(zz(i1:i2))
1554	! call maxdp_limits (z,n,zmax, j1,j2)
1555	zmax=maxval(z(j1:j2))
1556	! call maxdp_limits (zz,m,zzmax, i1,i2)
1557	zzmax=maxval(zz(i1:i2))
1558
1559	if(zzmin.lt.zmin)then
1560	write (,) 'from d interp: new variable out of limits'
1561	write (,) zzmin,'must be .ge. ',zmin
1562	stop
1563	! elseif(zzmax.gt.zmax)then
1564	! type *,'from interp: new variable out of limits'
1565	! type *,zzmax, 'must be .le. ',zmax
1566	! stop
1567	end if
1568
1569	do 1,i=i1,i2
1570
1571	do 2,j=j1,j2-1
1572	if(zz(i).ge.z(j).and.zz(i).lt.z(j+1)) goto 3
1573	2 continue
1574	c in this case (zz(i2).eq.z(j2)) and j leaves the loop with j=j2-1+1=j2
1575	if(opt.eq.1)then
1576	yy(i)=y(j2-1)+(y(j2)-y(j2-1))*
1577	$ (zz(i)-z(j2-1))/(z(j2)-z(j2-1))
1578	elseif(opt.eq.2)then
1579	if(y(j2).eq.0.0d0.or.y(j2-1).eq.0.0d0)then
1580	yy(i)=0.0d0
1581	else
1582	yy(i)=exp(log(y(j2-1))+log(y(j2)/y(j2-1))*
1583	@ (zz(i)-z(j2-1))/(z(j2)-z(j2-1)))
1584	end if
1585	else
1586	write (,) ' d interp : opt must be 1 or 2, opt= ',opt
1587	end if
1588	goto 1
1589	3 continue
1590	if(opt.eq.1)then
1591	yy(i)=y(j)+(y(j+1)-y(j))*(zz(i)-z(j))/(z(j+1)-z(j))
1592	! type *, ' '
1593	! type *, ' z(j),z(j+1) =', z(j),z(j+1)
1594	! type *, ' t(j),t(j+1) =', y(j),y(j+1)
1595	! type *, ' zz, tt = ', zz(i), yy(i)
1596	elseif(opt.eq.2)then
1597	if(y(j+1).eq.0.0d0.or.y(j).eq.0.0d0)then
1598	yy(i)=0.0d0
1599	else
1600	yy(i)=exp(log(y(j))+log(y(j+1)/y(j))*
1601	@ (zz(i)-z(j))/(z(j+1)-z(j)))
1602	end if
1603	else
1604	write (,) ' interp : opt must be 1 or 2, opt= ',opt
1605	end if
1606	1 continue
1607	return
1608	end
1609
1610
1611
1612
1613	c*****************************************************************************
1614	c Subroutines previously included in tcrco2_subrut.F
1615	c***********************************************************************
1616	c tcrco2_subrut.f
1617	c
1618	c jan 98 malv version for mz1d. copied from solar10/mz4sub.f
1619	c jul 2011 malv+fgg adapted to LMD-MGCM
1620	c***********************************************************************
1621
1622	************************************************************************
1623
1624	subroutine dinterconnection ( v, vt )
1625
1626	* input: vib. temp. from che*.for programs, vt(nl)
1627	* output: test vibrational temp. for other che*.for, v(nl)
1628	! iconex_smooth=1 ==> with smoothing
1629	! iconex_smooth=0 ==> without smoothing
1630	! iconex_tk=40 ==> with forced lte up to 40 km
1631	! iconex_tk=20 ==> with forced lte up to 20 km
1632	************************************************************************
1633
1634	implicit none
1635	include 'nlte_paramdef.h'
1636	include 'nlte_commons.h'
1637
1638	c argumentos
1639	real*8 vt(nl), v(nl)
1640
1641	c local variables
1642	integer i
1643
1644	c *************
1645
1646	do i=1,nl
1647	v(i) = vt(i)
1648	end do
1649
1650	! lo siguiente se utilizaba en solar10, pero es mejor introducirlo en
1651	! la driver. por ahora no lo uso todavia.
1652	! call fluctua(v,iconex_fluctua)
1653	! call smooth_nl(v,iconex_smooth,nl)
1654	! call forzar_tk(v,iconex_tk)
1655
1656	return
1657	end
1658
1659	c***********************************************************************
1660	function planckdp(tp,xnu)
1661	c returns the black body function at wavenumber xnu and temperature t.
1662	c***********************************************************************
1663
1664	implicit none
1665
1666	include 'nlte_paramdef.h'
1667	include 'nlte_commons.h'
1668
1669	! common/datis/ pi, vlight, ee, hplanck, gamma, ab,
1670	! @ n_avog, GG, R0, cte_sb, kboltzman, raddeg
1671	! real*8 pi, vlight, ee, hplanck, gamma, ab,
1672	! @ n_avog, GG, R0, cte_sb, kboltzman, raddeg
1673
1674	real*8 planckdp
1675	real*8 xnu
1676	real tp
1677
1678	planckdp = gammaxnu3.0 / exp( eexnu/dble(tp) )
1679	!erg cm-2.sr-1/cm-1.
1680
1681	return
1682	end
1683
1684	c ****************************************************************
1685	function bandid (ib)
1686	c returns the 2 character code of the band
1687	c ****************************************************************
1688	implicit none
1689
1690	integer ib
1691	character*2 bandid
1692
1693	132 format(i2)
1694	! encode (2,132,bandid) ib
1695	write ( bandid, 132) ib
1696
1697	if ( ib .eq. 1 ) bandid = '01'
1698	if ( ib .eq. 2 ) bandid = '02'
1699	if ( ib .eq. 3 ) bandid = '03'
1700	if ( ib .eq. 4 ) bandid = '04'
1701	if ( ib .eq. 5 ) bandid = '05'
1702	if ( ib .eq. 6 ) bandid = '06'
1703	if ( ib .eq. 7 ) bandid = '07'
1704	if ( ib .eq. 8 ) bandid = '08'
1705	if ( ib .eq. 9 ) bandid = '09'
1706	if ( ib .eq. 0 ) bandid = '00'
1707
1708	c end
1709	return
1710	end
1711
1712
1713
1714	c*****************************************************************************
1715	c Subroutines previously included in mat_oper.F
1716	c*****************************************************************************
1717	c set of subroutines for the cz*.for programs:
1718	! subroutine unit(a,n)
1719	! subroutine diago(a,v,n) diagonal matrix with v
1720	! subroutine invdiag(a,b,n) inverse of diagonal matrix
1721	! subroutine sypvvv(a,b,c,d,n) suma y prod de 3 vectores, muy comun
1722	! subroutine sypvmv(v,w,b,u,n) suma y prod de 3 vectores, muy comun
1723	! subroutine mulmvv(w,b,u,v,n) prod matriz vector vector
1724	! subroutine muymvv(w,b,u,v,n) prod matriz (inv.vector) vector
1725	! subroutine samem (a,m,n)
1726	! subroutine mulmv(a,b,c,n)
1727	! subroutine mulmm(a,b,c,n)
1728	! subroutine resmm(a,b,c,n)
1729	! subroutine mulvv(a,b,c,n)
1730	! subroutine sumvv(a,b,c,n)
1731	! subroutine zerom(a,n)
1732	! subroutine zero4m(a,b,c,d,n)
1733	! subroutine zero3m(a,b,c,n)
1734	! subroutine zero2m(a,b,n)
1735	! subroutine zerov(a,n)
1736	! subroutine zero4v(a,b,c,d,n)
1737	! subroutine zero3v(a,b,c,n)
1738	! subroutine zero2v(a,b,n)
1739
1740	!
1741	!
1742	! May-05 Sustituimos todos los zerojt de cristina por las subrutinas
1743	! genericas zerov***
1744	!
1745	c ***********************************************************************
1746	subroutine unit(a,n)
1747	c store the unit value in the diagonal of a
1748	c ***********************************************************************
1749	real*8 a(n,n)
1750	integer n,i,j,k
1751	do 1,i=2,n-1
1752	do 2,j=2,n-1
1753	if(i.eq.j) then
1754	a(i,j) = 1.d0
1755	else
1756	a(i,j)=0.0d0
1757	end if
1758	2 continue
1759	1 continue
1760	do k=1,n
1761	a(n,k) = 0.0d0
1762	a(1,k) = 0.0d0
1763	a(k,1) = 0.0d0
1764	a(k,n) = 0.0d0
1765	end do
1766	return
1767	end
1768
1769	c ***********************************************************************
1770	subroutine diago(a,v,n)
1771	c store the vector v in the diagonal elements of the square matrix a
1772	c ***********************************************************************
1773	implicit none
1774
1775	integer n,i,j,k
1776	real*8 a(n,n),v(n)
1777
1778	do 1,i=2,n-1
1779	do 2,j=2,n-1
1780	if(i.eq.j) then
1781	a(i,j) = v(i)
1782	else
1783	a(i,j)=0.0d0
1784	end if
1785	2 continue
1786	1 continue
1787	do k=1,n
1788	a(n,k) = 0.0d0
1789	a(1,k) = 0.0d0
1790	a(k,1) = 0.0d0
1791	a(k,n) = 0.0d0
1792	end do
1793	return
1794	end
1795
1796	c ***********************************************************************
1797	subroutine samem (a,m,n)
1798	c store the matrix m in the matrix a
1799	c ***********************************************************************
1800	real*8 a(n,n),m(n,n)
1801	integer n,i,j,k
1802	do 1,i=2,n-1
1803	do 2,j=2,n-1
1804	a(i,j) = m(i,j)
1805	2 continue
1806	1 continue
1807	do k=1,n
1808	a(n,k) = 0.0d0
1809	a(1,k) = 0.0d0
1810	a(k,1) = 0.0d0
1811	a(k,n) = 0.0d0
1812	end do
1813	return
1814	end
1815	c ***********************************************************************
1816	subroutine mulmv(a,b,c,n)
1817	c do a(i)=b(i,j)*c(j). a, b, and c must be distint
1818	c ***********************************************************************
1819	implicit none
1820
1821	integer n,i,j
1822	real*8 a(n),b(n,n),c(n),sum
1823
1824	do 1,i=2,n-1
1825	sum=0.0d0
1826	do 2,j=2,n-1
1827	sum=sum+ (b(i,j)) * (c(j))
1828	2 continue
1829	a(i)=sum
1830	1 continue
1831	a(1) = 0.0d0
1832	a(n) = 0.0d0
1833	return
1834	end
1835
1836	c ***********************************************************************
1837	subroutine mulmm(a,b,c,n)
1838	c ***********************************************************************
1839	real*8 a(n,n), b(n,n), c(n,n)
1840	integer n,i,j,k
1841
1842	! do i=2,n-1
1843	! do j=2,n-1
1844	! a(i,j)= 0.d00
1845	! do k=2,n-1
1846	! a(i,j) = a(i,j) + b(i,k) * c(k,j)
1847	! end do
1848	! end do
1849	! end do
1850	do j=2,n-1
1851	do i=2,n-1
1852	a(i,j)=0.d0
1853	enddo
1854	do k=2,n-1
1855	do i=2,n-1
1856	a(i,j)=a(i,j)+b(i,k)*c(k,j)
1857	enddo
1858	enddo
1859	enddo
1860	do k=1,n
1861	a(n,k) = 0.0d0
1862	a(1,k) = 0.0d0
1863	a(k,1) = 0.0d0
1864	a(k,n) = 0.0d0
1865	end do
1866
1867	return
1868	end
1869
1870	c ***********************************************************************
1871	subroutine resmm(a,b,c,n)
1872	c ***********************************************************************
1873	real*8 a(n,n), b(n,n), c(n,n)
1874	integer n,i,j,k
1875
1876	do i=2,n-1
1877	do j=2,n-1
1878	a(i,j)= b(i,j) - c(i,j)
1879	end do
1880	end do
1881	do k=1,n
1882	a(n,k) = 0.0d0
1883	a(1,k) = 0.0d0
1884	a(k,1) = 0.0d0
1885	a(k,n) = 0.0d0
1886	end do
1887
1888	return
1889	end
1890
1891	c ***********************************************************************
1892	subroutine sumvv(a,b,c,n)
1893	c a(i)=b(i)+c(i)
1894	c ***********************************************************************
1895	implicit none
1896
1897	integer n,i
1898	real*8 a(n),b(n),c(n)
1899
1900	do 1,i=2,n-1
1901	a(i)= (b(i)) + (c(i))
1902	1 continue
1903	a(1) = 0.0d0
1904	a(n) = 0.0d0
1905	return
1906	end
1907
1908	c ***********************************************************************
1909	subroutine zerom(a,n)
1910	c a(i,j)= 0.0
1911	c ***********************************************************************
1912
1913	implicit none
1914
1915	integer n,i,j
1916	real*8 a(n,n)
1917
1918	do 1,i=1,n
1919	do 2,j=1,n
1920	a(i,j) = 0.0d0
1921	2 continue
1922	1 continue
1923	return
1924	end
1925
1926	c ***********************************************************************
1927	subroutine zero4m(a,b,c,d,n)
1928	c a(i,j) = b(i,j) = c(i,j) = d(i,j) = 0.0
1929	c ***********************************************************************
1930	real*8 a(n,n), b(n,n), c(n,n), d(n,n)
1931	integer n,i,j
1932	do 1,i=1,n
1933	do 2,j=1,n
1934	a(i,j) = 0.0d0
1935	b(i,j) = 0.0d0
1936	c(i,j) = 0.0d0
1937	d(i,j) = 0.0d0
1938	2 continue
1939	1 continue
1940	return
1941	end
1942
1943	c ***********************************************************************
1944	subroutine zero3m(a,b,c,n)
1945	c a(i,j) = b(i,j) = c(i,j) = 0.0
1946	c **********************************************************************
1947	real*8 a(n,n), b(n,n), c(n,n)
1948	integer n,i,j
1949	do 1,i=1,n
1950	do 2,j=1,n
1951	a(i,j) = 0.0d0
1952	b(i,j) = 0.0d0
1953	c(i,j) = 0.0d0
1954	2 continue
1955	1 continue
1956	return
1957	end
1958
1959	c ***********************************************************************
1960	subroutine zero2m(a,b,n)
1961	c a(i,j) = b(i,j) = 0.0
1962	c ***********************************************************************
1963	real*8 a(n,n), b(n,n)
1964	integer n,i,j
1965	do 1,i=1,n
1966	do 2,j=1,n
1967	a(i,j) = 0.0d0
1968	b(i,j) = 0.0d0
1969	2 continue
1970	1 continue
1971	return
1972	end
1973	c ***********************************************************************
1974	subroutine zerov(a,n)
1975	c a(i)= 0.0
1976	c ***********************************************************************
1977	real*8 a(n)
1978	integer n,i
1979	do 1,i=1,n
1980	a(i) = 0.0d0
1981	1 continue
1982	return
1983	end
1984	c ***********************************************************************
1985	subroutine zero4v(a,b,c,d,n)
1986	c a(i) = b(i) = c(i) = d(i,j) = 0.0
1987	c ***********************************************************************
1988	real*8 a(n), b(n), c(n), d(n)
1989	integer n,i
1990	do 1,i=1,n
1991	a(i) = 0.0d0
1992	b(i) = 0.0d0
1993	c(i) = 0.0d0
1994	d(i) = 0.0d0
1995	1 continue
1996	return
1997	end
1998	c ***********************************************************************
1999	subroutine zero3v(a,b,c,n)
2000	c a(i) = b(i) = c(i) = 0.0
2001	c ***********************************************************************
2002	real*8 a(n), b(n), c(n)
2003	integer n,i
2004	do 1,i=1,n
2005	a(i) = 0.0d0
2006	b(i) = 0.0d0
2007	c(i) = 0.0d0
2008	1 continue
2009	return
2010	end
2011	c ***********************************************************************
2012	subroutine zero2v(a,b,n)
2013	c a(i) = b(i) = 0.0
2014	c ***********************************************************************
2015	real*8 a(n), b(n)
2016	integer n,i
2017	do 1,i=1,n
2018	a(i) = 0.0d0
2019	b(i) = 0.0d0
2020	1 continue
2021	return
2022	end
2023
2024

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

Download in other formats:

Original Format