Context Navigation

nlte_aux.F @ 498

Last change on this file since 498 was 498, checked in by emillour, 13 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

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 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