Context Navigation

sugas_corrk.F90 @ 191

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

Line
1	subroutine sugas_corrk
2
3	!==================================================================
4	!
5	! Purpose
6	! -------
7	! Set up gaseous absorption parameters used by the radiation code.
8	! This subroutine is a replacement for the old 'setrad', which contained
9	! both absorption and scattering data.
10	!
11	! Authors
12	! -------
13	! Adapted and generalised from the NASA Ames code by Robin Wordsworth (2009)
14	!
15	! Summary
16	! -------
17	!
18	!==================================================================
19
20	use radinc_h
21	use radcommon_h, only : pgasref,pfgasref,pgasmin,pgasmax
22	use radcommon_h, only : tgasref,tgasmin,tgasmax
23	use radcommon_h, only : gasv,gasi,FZEROI,FZEROV,gweight
24	use radcommon_h, only : wrefvar,gastype
25
26	implicit none
27
28	#include "datafile.h"
29	#include "callkeys.h"
30
31	!==================================================================
32
33	logical file_ok
34
35	integer n, nt, np, nh, ng, nw, m, i
36	integer L_NGAUSScheck, L_NPREFcheck, L_NTREFcheck, L_REFVARcheck
37
38	character(len=100) :: file_id
39	character(len=100) :: file_path
40
41	real*8 gasi8(L_NTREF,L_NPREF,L_REFVAR,L_NSPECTI,L_NGAUSS)
42	real*8 gasv8(L_NTREF,L_NPREF,L_REFVAR,L_NSPECTV,L_NGAUSS)
43
44	real*8 x, xi(4), yi(4), ans, p
45
46	integer Nspecies
47
48
49	!=======================================================================
50	! Load variable species data, exit if we have wrong database
51	file_id='/corrk_data/' // corrkdir(1:LEN_TRIM(corrkdir)) // '/Q.dat'
52	file_path=datafile(1:LEN_TRIM(datafile))//file_id(1:LEN_TRIM(file_id))
53
54	! check that the file exists
55	inquire(FILE=file_path,EXIST=file_ok)
56	if(.not.file_ok) then
57	write(,)'The file ',file_path(1:LEN_TRIM(file_path))
58	write(,)'was not found by sugas_corrk.F90, exiting.'
59	call abort
60	endif
61
62	! check that database matches varactive toggle
63	open(111,file=file_path(1:LEN_TRIM(file_path)),form='formatted')
64	read(111,*) Nspecies
65
66	if(Nspecies.eq.1 .and. (varactive.or.varfixed))then
67	print*,'You have varactive/fixed=.true. but the database [', &
68	corrkdir(1:LEN_TRIM(corrkdir)), &
69	'] has no variable species, exiting.'
70	call abort
71	elseif(Nspecies.eq.2 .and. (.not.varactive) .and. (.not.varfixed))then
72	print*,'You have varactive and varfixed=.false. and the database [', &
73	corrkdir(1:LEN_TRIM(corrkdir)), &
74	'] has a variable species.'
75	call abort
76	elseif(Nspecies.gt.3 .or. Nspecies.lt.1)then
77	print*,Nspecies,' species in database [', &
78	corrkdir(1:LEN_TRIM(corrkdir)), &
79	'], radiative code cannot handle this.'
80	call abort
81	endif
82
83	do n=1,Nspecies
84	read(111,*) gastype(n)
85	print*,'Gas ',n,' is ',gastype(n)
86	enddo
87
88	! check the array size is correct, load the coefficients
89	open(111,file=file_path(1:LEN_TRIM(file_path)),form='formatted')
90	read(111,*) L_REFVARcheck
91	if(.not.(L_REFVARcheck.eq.L_REFVAR)) then
92	print*,'The size of your radiative transfer mixing ratio array does '
93	print*,'not match the value given in Q.dat, exiting.'
94	call abort
95	endif
96	read(111,*) wrefvar
97	close(111)
98
99	! display the values
100	print*,'Variable gas mixing ratios:'
101	do n=1,L_REFVAR
102	print*,n,'.',wrefvar(n),' kg/kg'
103	end do
104	print*,''
105
106
107	!=======================================================================
108	! Set the weighting in g-space
109
110	file_id='/corrk_data/' // corrkdir(1:LEN_TRIM(corrkdir)) // '/g.dat'
111	file_path=datafile(1:LEN_TRIM(datafile))//file_id(1:LEN_TRIM(file_id))
112
113	! check that the file exists
114	inquire(FILE=file_path,EXIST=file_ok)
115	if(.not.file_ok) then
116	write(,)'The file ',file_path(1:LEN_TRIM(file_path))
117	write(,)'was not found by sugas_corrk.F90, exiting.'
118	call abort
119	endif
120
121	! check the array size is correct, load the coefficients
122	open(111,file=file_path(1:LEN_TRIM(file_path)),form='formatted')
123	read(111,*) L_NGAUSScheck
124	if(.not.(L_NGAUSScheck.eq.L_NGAUSS)) then
125	print*,'The size of your radiative transfer g-space array does '
126	print*,'not match the value given in g.dat, exiting.'
127	call abort
128	endif
129	read(111,*) gweight
130	close(111)
131
132	! display the values
133	print*,'Correlated-k g-space grid:'
134	do n=1,L_NGAUSS
135	print*,n,'.',gweight(n)
136	end do
137	print*,''
138
139	!=======================================================================
140	! Set the reference pressure and temperature arrays. These are
141	! the pressures and temperatures at which we have k-coefficients.
142
143	!-----------------------------------------------------------------------
144	! pressure
145
146	file_id='/corrk_data/' // corrkdir(1:LEN_TRIM(corrkdir)) // '/p.dat'
147	file_path=datafile(1:LEN_TRIM(datafile))//file_id(1:LEN_TRIM(file_id))
148
149	! check that the file exists
150	inquire(FILE=file_path,EXIST=file_ok)
151	if(.not.file_ok) then
152	write(,)'The file ',file_path(1:LEN_TRIM(file_path))
153	write(,)'was not found by sugas_corrk.F90, exiting.'
154	call abort
155	endif
156
157	! check the array size is correct, load the coefficients
158	open(111,file=file_path(1:LEN_TRIM(file_path)),form='formatted')
159	read(111,*) L_NPREFcheck
160	if(.not.(L_NPREFcheck.eq.L_NPREF)) then
161	print*,'The size of your radiative transfer pressure array does '
162	print*,'not match the value given in p.dat, exiting.'
163	call abort
164	endif
165	read(111,*) pgasref
166	close(111)
167
168	! display the values
169	print*,'Correlated-k pressure grid (mBar):'
170	do n=1,L_NPREF
171	print*,n,'. 1 x 10^',pgasref(n),' mBar'
172	end do
173	print*,''
174
175	! save the min / max matrix values
176	pgasmin = 10.0**pgasref(1)
177	pgasmax = 10.0**pgasref(L_NPREF)
178
179	! interpolate to finer grid
180	do n=1,L_NPREF-1
181	do m=1,5
182	pfgasref((n-1)5+m) = pgasref(n)+(m-1)0.2
183	end do
184	end do
185	pfgasref(L_PINT) = pgasref(L_NPREF)
186	! Warning! this may need to be generalised if we want to use uneven grids!
187
188
189	!-----------------------------------------------------------------------
190	! temperature
191
192	file_id='/corrk_data/' // corrkdir(1:LEN_TRIM(corrkdir)) // '/T.dat'
193	file_path=datafile(1:LEN_TRIM(datafile))//file_id(1:LEN_TRIM(file_id))
194
195	! check that the file exists
196	inquire(FILE=file_path,EXIST=file_ok)
197	if(.not.file_ok) then
198	write(,)'The file ',file_path(1:LEN_TRIM(file_path))
199	write(,)'was not found by sugas_corrk.F90, exiting.'
200	call abort
201	endif
202
203	! check the array size is correct, load the coefficients
204	open(111,file=file_path(1:LEN_TRIM(file_path)),form='formatted')
205	read(111,*) L_NTREFcheck
206	if(.not.(L_NTREFcheck.eq.L_NTREF)) then
207	print*,'The size of your radiative transfer temperature array does '
208	print*,'not match the value given in T.dat, exiting.'
209	call abort
210	endif
211	read(111,*) tgasref
212	close(111)
213
214	! display the values
215	print*,'Correlated-k temperature grid:'
216	do n=1,L_NTREF
217	print*,n,'.',tgasref(n),' K'
218	end do
219
220	! save the min / max matrix values
221	tgasmin = tgasref(1)
222	tgasmax = tgasref(L_NTREF)
223
224	!=======================================================================
225	! Get gaseous k-coefficients and interpolate onto finer pressure grid
226
227	! VISIBLE
228	if (callgasvis) then
229	file_id='/corrk_data/'//trim(adjustl(banddir))//'/corrk_gcm_VI.dat'
230	file_path=datafile(1:LEN_TRIM(datafile))//file_id(1:LEN_TRIM(file_id))
231
232	! check that the file exists
233	inquire(FILE=file_path,EXIST=file_ok)
234	if(.not.file_ok) then
235	write(,)'The file ',file_path(1:LEN_TRIM(file_path))
236	write(,)'was not found by sugas_corrk.F90.'
237	write(,)'Are you sure you have absorption data for these bands?'
238	call abort
239	endif
240
241	open(111,file=file_path(1:LEN_TRIM(file_path)),form='formatted')
242	read(111,*) gasv8
243	close(111)
244
245	else
246	print*,'Visible gaseous absorption is set to zero.'
247	endif
248
249	! INFRA-RED
250	file_id='/corrk_data/'//trim(adjustl(banddir))//'/corrk_gcm_IR.dat'
251	file_path=datafile(1:LEN_TRIM(datafile))//file_id(1:LEN_TRIM(file_id))
252
253	! check that the file exists
254	inquire(FILE=file_path,EXIST=file_ok)
255	if(.not.file_ok) then
256	write(,)'The file ',file_path(1:LEN_TRIM(file_path))
257	write(,)'was not found by sugas_corrk.F90.'
258	write(,)'Are you sure you have absorption data for these bands?'
259	call abort
260	endif
261
262	open(111,file=file_path(1:LEN_TRIM(file_path)),form='formatted')
263	read(111,*) gasi8
264	close(111)
265
266	do nw=1,L_NSPECTI
267	fzeroi(nw) = 0.
268	end do
269	do nw=1,L_NSPECTV
270	fzerov(nw) = 0.
271	end do
272
273	! Take log10 of the values - this is what we will interpolate.
274	! Smallest value is 1.0E-200.
275
276	do nt=1,L_NTREF
277	do np=1,L_NPREF
278	do nh=1,L_REFVAR
279	do ng = 1,L_NGAUSS
280
281	do nw=1,L_NSPECTV
282	if(gasv8(nt,np,nh,nw,ng).gt.1.0d-200) then
283	gasv8(nt,np,nh,nw,ng) = &
284	log10(gasv8(nt,np,nh,nw,ng))
285	else
286	gasv8(nt,np,nh,nw,ng) = -200.0
287	end if
288	end do
289
290	do nw=1,L_NSPECTI
291	if(gasi8(nt,np,nh,nw,ng).gt.1.0d-200) then
292	gasi8(nt,np,nh,nw,ng) = &
293	log10(gasi8(nt,np,nh,nw,ng))
294	else
295	gasi8(nt,np,nh,nw,ng) = -200.0
296	end if
297	end do
298
299	end do
300	end do
301	end do
302	end do
303
304	! Interpolate the values: first the longwave
305
306	do nt=1,L_NTREF
307	do nh=1,L_REFVAR
308	do nw=1,L_NSPECTI
309	do ng=1,L_NGAUSS
310
311	! First, the initial interval
312
313	n = 1
314	do m=1,5
315	x = pfgasref(m)
316	xi(1) = pgasref(n)
317	xi(2) = pgasref(n+1)
318	xi(3) = pgasref(n+2)
319	xi(4) = pgasref(n+3)
320	yi(1) = gasi8(nt,n,nh,nw,ng)
321	yi(2) = gasi8(nt,n+1,nh,nw,ng)
322	yi(3) = gasi8(nt,n+2,nh,nw,ng)
323	yi(4) = gasi8(nt,n+3,nh,nw,ng)
324	call lagrange(x,xi,yi,ans)
325	gasi(nt,m,nh,nw,ng) = 10.0**ans
326	end do
327
328	do n=2,L_NPREF-2
329	do m=1,5
330	i = (n-1)*5+m
331	x = pfgasref(i)
332	xi(1) = pgasref(n-1)
333	xi(2) = pgasref(n)
334	xi(3) = pgasref(n+1)
335	xi(4) = pgasref(n+2)
336	yi(1) = gasi8(nt,n-1,nh,nw,ng)
337	yi(2) = gasi8(nt,n,nh,nw,ng)
338	yi(3) = gasi8(nt,n+1,nh,nw,ng)
339	yi(4) = gasi8(nt,n+2,nh,nw,ng)
340	call lagrange(x,xi,yi,ans)
341	gasi(nt,i,nh,nw,ng) = 10.0**ans
342	end do
343	end do
344
345	! Now, get the last interval
346
347	n = L_NPREF-1
348	do m=1,5
349	i = (n-1)*5+m
350	x = pfgasref(i)
351	xi(1) = pgasref(n-2)
352	xi(2) = pgasref(n-1)
353	xi(3) = pgasref(n)
354	xi(4) = pgasref(n+1)
355	yi(1) = gasi8(nt,n-2,nh,nw,ng)
356	yi(2) = gasi8(nt,n-1,nh,nw,ng)
357	yi(3) = gasi8(nt,n,nh,nw,ng)
358	yi(4) = gasi8(nt,n+1,nh,nw,ng)
359	call lagrange(x,xi,yi,ans)
360	gasi(nt,i,nh,nw,ng) = 10.0**ans
361	end do
362
363	! Fill the last pressure point
364
365	gasi(nt,L_PINT,nh,nw,ng) = &
366	10.0**gasi8(nt,L_NPREF,nh,nw,ng)
367
368	end do
369	end do
370	end do
371	end do
372
373	! Interpolate the values: now the shortwave
374
375	do nt=1,L_NTREF
376	do nh=1,L_REFVAR
377	do nw=1,L_NSPECTV
378	do ng=1,L_NGAUSS
379
380	! First, the initial interval
381
382	n = 1
383	do m=1,5
384	x = pfgasref(m)
385	xi(1) = pgasref(n)
386	xi(2) = pgasref(n+1)
387	xi(3) = pgasref(n+2)
388	xi(4) = pgasref(n+3)
389	yi(1) = gasv8(nt,n,nh,nw,ng)
390	yi(2) = gasv8(nt,n+1,nh,nw,ng)
391	yi(3) = gasv8(nt,n+2,nh,nw,ng)
392	yi(4) = gasv8(nt,n+3,nh,nw,ng)
393	call lagrange(x,xi,yi,ans)
394	gasv(nt,m,nh,nw,ng) = 10.0**ans
395	end do
396
397	do n=2,L_NPREF-2
398	do m=1,5
399	i = (n-1)*5+m
400	x = pfgasref(i)
401	xi(1) = pgasref(n-1)
402	xi(2) = pgasref(n)
403	xi(3) = pgasref(n+1)
404	xi(4) = pgasref(n+2)
405	yi(1) = gasv8(nt,n-1,nh,nw,ng)
406	yi(2) = gasv8(nt,n,nh,nw,ng)
407	yi(3) = gasv8(nt,n+1,nh,nw,ng)
408	yi(4) = gasv8(nt,n+2,nh,nw,ng)
409	call lagrange(x,xi,yi,ans)
410	gasv(nt,i,nh,nw,ng) = 10.0**ans
411	end do
412	end do
413
414	! Now, get the last interval
415
416	n = L_NPREF-1
417	do m=1,5
418	i = (n-1)*5+m
419	x = pfgasref(i)
420	xi(1) = pgasref(n-2)
421	xi(2) = pgasref(n-1)
422	xi(3) = pgasref(n)
423	xi(4) = pgasref(n+1)
424	yi(1) = gasv8(nt,n-2,nh,nw,ng)
425	yi(2) = gasv8(nt,n-1,nh,nw,ng)
426	yi(3) = gasv8(nt,n,nh,nw,ng)
427	yi(4) = gasv8(nt,n+1,nh,nw,ng)
428	call lagrange(x,xi,yi,ans)
429	gasv(nt,i,nh,nw,ng) = 10.0**ans
430	end do
431
432	! Fill the last pressure point
433
434	gasv(nt,L_PINT,nh,nw,ng) = &
435	10.0**gasv8(nt,L_NPREF,nh,nw,ng)
436
437	end do
438	end do
439	end do
440	end do
441
442	return
443	end subroutine sugas_corrk

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

Download in other formats:

Original Format