Context Navigation

setspi.F90 @ 146

Last change on this file since 146 was 135, checked in by aslmd, 14 years ago
CHANGEMENT ARBORESCENCE ETAPE 2 -- NON COMPLET
File size: 6.5 KB

Line
1	subroutine setspi
2
3	!==================================================================
4	!
5	! Purpose
6	! -------
7	! Set up spectral intervals and Planck function in the longwave.
8	!
9	! Authors
10	! -------
11	! Adapted from setspi in the NASA Ames radiative code by
12	! Robin Wordsworth (2009).
13	!
14	! Called by
15	! ---------
16	! callcorrk.F
17	!
18	! Calls
19	! -----
20	! none
21	!
22	!==================================================================
23
24	use radinc_h, only: L_NSPECTI,corrkdir,banddir,NTstar,NTstop
25	use radcommon_h, only: BWNI,BLAMI,WNOI,DWNI,WAVEI,planckir,sigma
26
27	implicit none
28
29	#include "callkeys.h"
30	#include "comcstfi.h"
31	#include "datafile.h"
32
33	logical file_ok
34	integer nw, nt, m, mm, file_entries
35	real*8 a, b, ans, y, bpa, bma, T
36
37	character(len=30) :: temp1
38	character(len=100) :: file_id
39	character(len=100) :: file_path
40
41	! C1 and C2 values from Goody and Yung (2nd edition) MKS units
42	! These values lead to a "sigma" (sigma*T^4) of 5.67032E-8 W m^-2 K^-4
43
44	real*8 :: c1 = 3.741832D-16 ! W m^-2
45	real*8 :: c2 = 1.438786D-2 ! m K
46
47	real*8 :: lastband(2), plancksum
48
49	logical forceEC, planckcheck
50
51	real*8 :: x(12) = [ -0.981560634246719D0, -0.904117256370475D0, &
52	-0.769902674194305D0, -0.587317954286617D0, &
53	-0.367831498998180D0, -0.125233408511469D0, &
54	0.125233408511469D0, 0.367831498998180D0, &
55	0.587317954286617D0, 0.769902674194305D0, &
56	0.904117256370475D0, 0.981560634246719D0 ]
57
58	real*8 :: w(12) = [ 0.047175336386512D0, 0.106939325995318D0, &
59	0.160078328543346D0, 0.203167426723066D0, &
60	0.233492536538355D0, 0.249147045813403D0, &
61	0.249147045813403D0, 0.233492536538355D0, &
62	0.203167426723066D0, 0.160078328543346D0, &
63	0.106939325995318D0, 0.047175336386512D0 ]
64	mm=0
65
66	forceEC=.false.
67	planckcheck=.true.
68
69	!=======================================================================
70	! Set up spectral bands - wavenumber [cm^(-1)]. Go from smaller to
71	! larger wavenumbers.
72
73	write(temp1,'(i2.2)') L_NSPECTI
74	!file_id='/corrk_data/' // corrkdir(1:LEN_TRIM(corrkdir)) // '/narrowbands_IR.in'
75	file_id='/corrk_data/'//trim(adjustl(banddir))//'/narrowbands_IR.in'
76	file_path=datafile(1:LEN_TRIM(datafile))//file_id(1:LEN_TRIM(file_id))
77
78	! check that the file exists
79	inquire(FILE=file_path,EXIST=file_ok)
80	if(.not.file_ok) then
81	write(,)'The file ',file_path(1:LEN_TRIM(file_path))
82	write(,)'was not found by setspi.F90, exiting.'
83	call abort
84	endif
85
86	! check that the file contains the right number of bands
87	call system('wc -l '//file_path//' > bandlen.txt')
88	open(131,file='bandlen.txt', form='formatted')
89	read(131,*) file_entries
90	close(131)
91	call system('rm -f bandlen.txt')
92
93	if(file_entries-1.ne.L_NSPECTI) then
94	write(,) 'L_NSPECTI = ',L_NSPECTI, 'in the model, but there are '
95	write(,) file_entries-1,'entries in ', &
96	file_path(1:LEN_TRIM(file_path)),', exiting.'
97	call abort
98	endif
99
100	! load and display the data
101	open(111,file=file_path,form='formatted')
102	read(111,*)
103	do M=1,L_NSPECTI-1
104	read(111,*) BWNI(M)
105	end do
106	read(111,*) lastband
107	close(111)
108	BWNI(L_NSPECTI) =lastband(1)
109	BWNI(L_NSPECTI+1)=lastband(2)
110
111	print*,''
112	print*,'IR band limits:'
113	do M=1,L_NSPECTI+1
114	print*,m,'-->',BWNI(M),' cm^-1'
115	end do
116
117	! Set up mean wavenumbers and wavenumber deltas. Units of
118	! wavenumbers is cm^(-1); units of wavelengths is microns.
119
120	do M=1,L_NSPECTI
121	WNOI(M) = 0.5*(BWNI(M+1)+BWNI(M))
122	DWNI(M) = BWNI(M+1)-BWNI(M)
123	WAVEI(M) = 1.0E+4/WNOI(M)
124	BLAMI(M) = 0.01/BWNI(M)
125	end do
126	BLAMI(M) = 0.01/BWNI(M)
127	! note M=L_NSPECTI+1 after loop due to Fortran bizarreness
128
129	!=======================================================================
130	! For each IR wavelength interval, compute the integral of B(T), the
131	! Planck function, divided by the wavelength interval, in cm-1. The
132	! integration is in MKS units, the final answer is the same as the
133	! original planck.f; W m^-2 wavenumber^-1, where wavenumber is in CM^-1.
134
135	print*,''
136	print*,'Current Planck integration range:'
137	print*,'T = ',dble(NTstar)/1.0D+1, ' to ',dble(NTstop)/1.0D+1,' K.'
138
139	do NW=1,L_NSPECTI
140	a = 1.0D-2/BWNI(NW+1)
141	b = 1.0D-2/BWNI(NW)
142	bpa = (b+a)/2.0
143	bma = (b-a)/2.0
144	do nt=NTstar,NTstop
145	T = dble(NT)/1.0D+1
146	ans = 0.0D0
147
148	do mm=1,12
149	y = bma*x(mm)+bpa
150	ans = ans + w(mm)c1/(y5(exp(c2/(y*T))-1.0D0))
151	end do
152
153	planckir(NW,nt-NTstar+1) = ansbma/(PIDWNI(NW))
154	end do
155	end do
156
157	! force planck=sigmaepsT^4 for each temperature in array
158	if(forceEC)then
159	do nt=NTstar,NTstop
160	plancksum=0.0
161	T=dble(NT)/1.0D+1
162
163	do NW=1,L_NSPECTI
164	plancksum=plancksum+ &
165	planckir(NW,nt-NTstar+1)DWNI(NW)pi
166	end do
167
168	do NW=1,L_NSPECTI
169	planckir(NW,nt-NTstar+1)= &
170	planckir(NW,nt-NTstar+1)* &
171	sigma(dble(nt)/1.0D+1)*4/plancksum
172	end do
173	end do
174	endif
175
176	if(planckcheck)then
177	! check energy conservation at lower temperature boundary
178	plancksum=0.0
179	nt=NTstar
180	do NW=1,L_NSPECTI
181	plancksum=plancksum+planckir(NW,nt-NTstar+1)DWNI(NW)pi
182	end do
183	print*,'At lower limit:'
184	print,'in model sigT^4 = ',plancksum,' W m^-2'
185	print,'actual sigT^4 = ',sigma(dble(nt)/1.0D+1)*4,' W m^-2'
186
187	! check energy conservation at upper temperature boundary
188	plancksum=0.0
189	nt=NTstop
190	do NW=1,L_NSPECTI
191	plancksum=plancksum+planckir(NW,nt-NTstar+1)DWNI(NW)pi
192	end do
193	print*,'At upper limit:'
194	print,'in model sigT^4 = ',plancksum,' W m^-2'
195	print,'actual sigT^4 = ',sigma(dble(nt)/1.0D+1)*4,' W m^-2'
196	print*,''
197	endif
198
199	return
200	end subroutine setspi

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

Download in other formats:

Original Format