Context Navigation

ave_stelspec.F90 @ 2613

Last change on this file since 2613 was 2422, checked in by yjaziri, 4 years ago
GENERIC GCM Add possibility to read stellar spectra from input file in "ave_stelspec.F90" with option "startype==11" in "callphys.def". Input file need to be in "datagcm/stellar_spectra/" and called in "callphys.def" by "stelspec_file=...". File format have one and only one header line. With this option you also need to specify in "callphys.def" "tstellar==..." YJ
File size: 8.4 KB

Line
1	subroutine ave_stelspec(STELLAR)
2
3	!==================================================================
4	!
5	! Purpose
6	! -------
7	! Average the chosen high resolution stellar spectrum over the
8	! visible bands in the model.
9	!
10	! Authors
11	! -------
12	! Robin Wordsworth (2010).
13	! Generalized to very late spectral types (and Brown dwarfs) Jeremy Leconte (2012)
14	!
15	! Called by
16	! ---------
17	! setspv.F
18	!
19	! Calls
20	! -----
21	! none
22	!
23	!==================================================================
24
25	use radinc_h, only: L_NSPECTV
26	use radcommon_h, only: BWNV, DWNV, tstellar
27	use datafile_mod, only: datadir
28	use callkeys_mod, only: stelbbody,stelTbb,startype
29	use ioipsl_getin_p_mod, only: getin_p
30
31	implicit none
32
33	real*8 STELLAR(L_NSPECTV)
34	! integer startype
35
36	integer Nfine
37	integer,parameter :: Nfineband=200
38	integer ifine,band
39
40	real,allocatable,save :: lam(:),stel_f(:) !read by master
41	real lamm,lamp
42	real dl
43
44	character(len=100) :: file_id,file_id_lam
45	character(len=200) :: file_path,file_path_lam
46	character(len=150) :: stelspec_file
47
48	real lam_temp
49	double precision stel_temp
50
51	integer :: ios ! file opening/reading status
52
53	STELLAR(:)=0.0
54
55	print*,'enter ave_stellspec'
56	if(stelbbody)then
57	tstellar=stelTbb
58	Nfine=L_NSPECTV*Nfineband
59	do band=1,L_NSPECTV
60	lamm=10000.0/BWNV(band+1)
61	lamp=10000.0/BWNV(band)
62	dl=(lamp-lamm)/(Nfineband)
63	do ifine=1,Nfineband
64	lam_temp=lamm+(lamp-lamm)*(ifine-1.)/(Nfineband)
65	call blackl(dble(lam_temp*1e-6),dble(tstellar),stel_temp)
66	STELLAR(band)=STELLAR(band)+stel_temp*dl
67	enddo
68	end do
69	STELLAR(1:L_NSPECTV)=STELLAR(1:L_NSPECTV)/sum(STELLAR(1:L_NSPECTV))
70	else
71	! load high resolution stellar data
72	Select Case(startype)
73	Case(1)
74	file_id='/stellar_spectra/sol.txt'
75	tstellar=5800.
76	file_id_lam='/stellar_spectra/lam.txt'
77	Nfine=5000
78	Case(2)
79	file_id='/stellar_spectra/gl581.txt'
80	tstellar=3200.
81	file_id_lam='/stellar_spectra/lam.txt'
82	Nfine=5000
83	Case(3)
84	file_id='/stellar_spectra/adleo.txt'
85	tstellar=3200.
86	file_id_lam='/stellar_spectra/lam.txt'
87	Nfine=5000
88	Case(4)
89	file_id='/stellar_spectra/gj644.txt'
90	print*,'Find out tstellar before using this star!'
91	call abort
92	file_id_lam='/stellar_spectra/lam.txt'
93	Nfine=5000
94	Case(5)
95	file_id='/stellar_spectra/hd128167.txt'
96	tstellar=6700. ! Segura et al. (2003)
97	file_id_lam='/stellar_spectra/lam.txt'
98	Nfine=5000
99	Case(6)
100	file_id='/stellar_spectra/BD_Teff-1600K.txt'
101	tstellar=1600.
102	file_id_lam='/stellar_spectra/lamBD.txt'
103	Nfine=5000
104	Case(7)
105	file_id='/stellar_spectra/BD_Teff-1000K.txt'
106	tstellar=1000.
107	file_id_lam='/stellar_spectra/lamBD.txt'
108	Nfine=5000
109	Case(8)
110	file_id='/stellar_spectra/Flux_K5_Teff4700_logg4.5_Met-0.5_BTsettle.dat'
111	tstellar=4700.
112	file_id_lam='/stellar_spectra/lambda_K5_Teff4700_logg4.5_Met-0.5_BTsettle.dat'
113	Nfine=3986
114	Case(9)
115	file_id='/stellar_spectra/Flux_TRAPPIST1.dat'
116	tstellar=2550.
117	file_id_lam='/stellar_spectra/lambda_TRAPPIST1.dat'
118	Nfine=5000
119	Case(10)
120	file_id='/stellar_spectra/Flux_Proxima.dat'
121	tstellar=3050.
122	file_id_lam='/stellar_spectra/lambda_Proxima.dat'
123	Nfine=5000
124	Case(11)
125	! look for a " stelspec_file= ..." option in def files
126	write(,) "Input stellar spectra files for radiative transfer is:"
127	stelspec_file = "sun.txt" ! default
128	call getin_p("stelspec_file",stelspec_file) ! default path
129	write(,) " stelspec_file = ",trim(stelspec_file)
130	write(,) 'Please use ',1,' and only ',1,' header line in ',trim(stelspec_file)
131	! look for a " tstellar= ..." option in def files
132	write(,) "Input stellar temperature is:"
133	tstellar = -1. ! default
134	call getin_p("tstellar",tstellar) ! default path
135	write(,) " tstellar = ",tstellar
136	if (tstellar.eq.-1.) then
137	write(,)'Error: with startype 11 tstellar need to be specified'
138	write(,)' Specified in callphys.def tstellar=...'
139	stop
140	end if
141
142	! Opening file
143	file_path = trim(datadir)//'/stellar_spectra/'//stelspec_file
144	print*, 'stellar flux : ', file_path
145	OPEN(UNIT=110,FILE=file_path,STATUS='old',iostat=ios)
146
147	if (ios /= 0) THEN
148	write(,)'Error: cannot open stelspec_file file ', trim(stelspec_file)
149	write(,)'It should be in :',trim(datadir),'/stellar_spectra/'
150	write(,)'1) You can change the data directory in callphys.def'
151	write(,)' with:'
152	write(,)' datadir=/path/to/the/directory'
153	write(,)'2) You can change the input stelspec_file file name in'
154	write(,)' callphys.def with:'
155	write(,)' stelspec_file=filename'
156	stop
157	end if
158
159	! Get number of line in the file
160	READ(110,*) ! skip header
161	Nfine = 0
162	do
163	read(110,*,iostat=ios)
164	if (ios<0) exit
165	Nfine = Nfine + 1
166	end do
167	rewind(110) ! Rewind file after counting lines
168	READ(110,*) ! skip header
169	Case Default
170	print*,'Error: unknown star type chosen'
171	call abort
172	End Select
173
174	!$OMP MASTER
175	allocate(lam(Nfine),stel_f(Nfine))
176
177	if (startype.eq.11) then
178	do ifine=1,Nfine
179	read(110,*) lam(ifine), stel_f(ifine) ! lam [um] stel_f [per unit of wavelength] (integrated and normalized by Fat1AU)
180	enddo
181	close(110)
182	else
183	file_path_lam=TRIM(datadir)//TRIM(file_id_lam)
184	open(110,file=file_path_lam,form='formatted',status='old',iostat=ios)
185	if (ios.ne.0) then ! file not found
186	write(,) 'Error from ave_stelspec'
187	write(,) 'Data file ',trim(file_id_lam),' not found.'
188	write(,)'Check that your path to datagcm:',trim(datadir)
189	write(,)' is correct. You can change it in callphys.def with:'
190	write(,)' datadir = /absolute/path/to/datagcm'
191	write(,)' Also check that there is a ',trim(file_id_lam),' there.'
192	call abort
193	else
194	do ifine=1,Nfine
195	read(110,*) lam(ifine)
196	enddo
197	close(110)
198	endif
199
200
201	! load high resolution wavenumber data
202	file_path=TRIM(datadir)//TRIM(file_id)
203	open(111,file=file_path,form='formatted',status='old',iostat=ios)
204	if (ios.ne.0) then ! file not found
205	write(,) 'Error from ave_stelspec'
206	write(,) 'Data file ',trim(file_id),' not found.'
207	write(,)'Check that your path to datagcm:',trim(datadir)
208	write(,)' is correct. You can change it in callphys.def with:'
209	write(,)' datadir = /absolute/path/to/datagcm'
210	write(,)' Also check that there is a ',trim(file_id),' there.'
211	call abort
212	else
213	do ifine=1,Nfine
214	read(111,*) stel_f(ifine)
215	enddo
216	close(111)
217	endif
218	end if
219	!$OMP END MASTER
220	!$OMP BARRIER
221
222	! sum data by band
223	band=1
224	Do while(lam(1).lt. real(10000.0/BWNV(band+1)))
225	if (band.gt.L_NSPECTV-1) exit
226	band=band+1
227	enddo
228	dl=lam(2)-lam(1)
229	STELLAR(band)=STELLAR(band)+stel_f(1)*dl
230	do ifine = 2,Nfine
231	if(lam(ifine) .gt. real(10000.0/BWNV(band)))then
232	band=band-1
233	endif
234	if(band .lt. 1) exit
235	dl=lam(ifine)-lam(ifine-1)
236	STELLAR(band)=STELLAR(band)+stel_f(ifine)*dl
237	end do
238
239
240	STELLAR(1:L_NSPECTV)=STELLAR(1:L_NSPECTV)/sum(STELLAR(1:L_NSPECTV))
241	!$OMP BARRIER
242	!$OMP MASTER
243	if (allocated(lam)) deallocate(lam)
244	if (allocated(stel_f)) deallocate(stel_f)
245	!$OMP END MASTER
246	!$OMP BARRIER
247	endif
248
249	end subroutine ave_stelspec

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

Download in other formats:

Original Format