source: trunk/LMDZ.PLUTO/libf/phypluto/ave_stelspec.F90 @ 3504

Last change on this file since 3504 was 3504, checked in by afalco, 13 days ago

Pluto: print only on master process.
AF

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