Context Navigation

ave_stelspec.F90 @ 3922

Last change on this file since 3922 was 3893, checked in by gmilcareck, 6 months ago
Remove all "call abort" and "stop" statement in LMDZ.GENERIC and replacing them by call abort_physic().
File size: 8.6 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	call abort_physic("GJ644", "Tstellar unknown", 1)
91	file_id_lam='/stellar_spectra/lam.txt'
92	Nfine=5000
93	Case(5)
94	file_id='/stellar_spectra/Wasp43_flux.txt'
95	tstellar=4520.
96	file_id_lam='/stellar_spectra/Wasp43_lambda.txt'
97	Nfine=395562
98	Case(6)
99	file_id='/stellar_spectra/BD_Teff-1600K.txt'
100	tstellar=1600.
101	file_id_lam='/stellar_spectra/lamBD.txt'
102	Nfine=5000
103	Case(7)
104	file_id='/stellar_spectra/BD_Teff-1000K.txt'
105	tstellar=1000.
106	file_id_lam='/stellar_spectra/lamBD.txt'
107	Nfine=5000
108	Case(8)
109	file_id='/stellar_spectra/Flux_K5_Teff4700_logg4.5_Met-0.5_BTsettle.dat'
110	tstellar=4700.
111	file_id_lam='/stellar_spectra/lambda_K5_Teff4700_logg4.5_Met-0.5_BTsettle.dat'
112	Nfine=3986
113	Case(9)
114	file_id='/stellar_spectra/Flux_TRAPPIST1.dat'
115	tstellar=2550.
116	file_id_lam='/stellar_spectra/lambda_TRAPPIST1.dat'
117	Nfine=5000
118	Case(10)
119	file_id='/stellar_spectra/Flux_Proxima.dat'
120	tstellar=3050.
121	file_id_lam='/stellar_spectra/lambda_Proxima.dat'
122	Nfine=5000
123	Case(11)
124	! look for a " stelspec_file= ..." option in def files
125	write(,) "Input stellar spectra files for radiative transfer is:"
126	stelspec_file = "sun.txt" ! default
127	call getin_p("stelspec_file",stelspec_file) ! default path
128	write(,) " stelspec_file = ",trim(stelspec_file)
129	write(,) 'Please use ',1,' and only ',1,' header line in ',trim(stelspec_file)
130	! look for a " tstellar= ..." option in def files
131	write(,) "Input stellar temperature is:"
132	tstellar = -1. ! default
133	call getin_p("tstellar",tstellar) ! default path
134	write(,) " tstellar = ",tstellar
135	if (tstellar.eq.-1.) then
136	write(,)'Error: with startype 11 tstellar need to be specified'
137	write(,)' Specified in callphys.def tstellar=...'
138	call abort_physic("ave_stelspec", "startype 11 tstellar needs to be specified",1)
139	end if
140
141	! Opening file
142	file_path = trim(datadir)//'/stellar_spectra/'//stelspec_file
143	print*, 'stellar flux : ', file_path
144	OPEN(UNIT=110,FILE=file_path,STATUS='old',iostat=ios)
145
146	if (ios /= 0) THEN
147	write(,)'Error: cannot open stelspec_file file ', trim(stelspec_file)
148	write(,)'It should be in :',trim(datadir),'/stellar_spectra/'
149	write(,)'1) You can change the data directory in callphys.def'
150	write(,)' with:'
151	write(,)' datadir=/path/to/the/directory'
152	write(,)'2) You can change the input stelspec_file file name in'
153	write(,)' callphys.def with:'
154	write(,)' stelspec_file=filename'
155	call abort_physic("ave_stelspec", "Unable to read stellar flux file", 1)
156	end if
157
158	! Get number of line in the file
159	READ(110,*) ! skip header
160	Nfine = 0
161	do
162	read(110,*,iostat=ios)
163	if (ios<0) exit
164	Nfine = Nfine + 1
165	end do
166	rewind(110) ! Rewind file after counting lines
167	READ(110,*) ! skip header
168	Case Default
169	call abort_physic("ave_stelspec", "Unknown star type chosen", 1)
170	End Select
171
172	!$OMP MASTER
173	allocate(lam(Nfine),stel_f(Nfine))
174
175	if (startype.eq.11) then
176	do ifine=1,Nfine
177	read(110,*) lam(ifine), stel_f(ifine) ! lam [um] stel_f [per unit of wavelength] (integrated and normalized by Fat1AU)
178	enddo
179	close(110)
180	else
181	file_path_lam=TRIM(datadir)//TRIM(file_id_lam)
182	open(110,file=file_path_lam,form='formatted',status='old',iostat=ios)
183	if (ios.ne.0) then ! file not found
184	write(,) 'Error from ave_stelspec'
185	write(,) 'Data file ',trim(file_id_lam),' not found.'
186	write(,)'Check that your path to datagcm:',trim(datadir)
187	write(,)' is correct. You can change it in callphys.def with:'
188	write(,)' datadir = /absolute/path/to/datagcm'
189	write(,)' Also check that there is a ',trim(file_id_lam),' there.'
190	call abort_physic("ave_stelspec", "Unable to find file", 1)
191	else
192	do ifine=1,Nfine
193	read(110,*) lam(ifine)
194	enddo
195	close(110)
196	endif
197
198
199	! load high resolution wavenumber data
200	file_path=TRIM(datadir)//TRIM(file_id)
201	open(111,file=file_path,form='formatted',status='old',iostat=ios)
202	if (ios.ne.0) then ! file not found
203	write(,) 'Error from ave_stelspec'
204	write(,) 'Data file ',trim(file_id),' not found.'
205	write(,)'Check that your path to datagcm:',trim(datadir)
206	write(,)' is correct. You can change it in callphys.def with:'
207	write(,)' datadir = /absolute/path/to/datagcm'
208	write(,)' Also check that there is a ',trim(file_id),' there.'
209	call abort_physic("ave_stelspec", "Unable to find file", 1)
210	else
211	do ifine=1,Nfine
212	read(111,*) stel_f(ifine)
213	enddo
214	close(111)
215	endif
216	end if
217	!$OMP END MASTER
218	!$OMP BARRIER
219
220	! sum data by band
221	band=1
222	Do while(lam(1).lt. real(10000.0/BWNV(band+1)))
223	if (band.gt.L_NSPECTV-1) exit
224	band=band+1
225	enddo
226	dl=lam(2)-lam(1)
227	STELLAR(band)=STELLAR(band)+stel_f(1)*dl
228	do ifine = 2,Nfine
229	if(lam(ifine) .gt. real(10000.0/BWNV(band)))then
230	band=band-1
231	endif
232	if(band .lt. 1) exit
233	dl=lam(ifine)-lam(ifine-1)
234	STELLAR(band)=STELLAR(band)+stel_f(ifine)*dl
235	end do
236
237
238	STELLAR(1:L_NSPECTV)=STELLAR(1:L_NSPECTV)/sum(STELLAR(1:L_NSPECTV))
239	!$OMP BARRIER
240	!$OMP MASTER
241	if (allocated(lam)) deallocate(lam)
242	if (allocated(stel_f)) deallocate(stel_f)
243	!$OMP END MASTER
244	!$OMP BARRIER
245	endif
246
247	end subroutine ave_stelspec

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

Download in other formats:

Original Format