Context Navigation

sugas_corrk.F90 @ 304

Last change on this file since 304 was 304, checked in by rwordsworth, 13 years ago
Option added for FZERO but not implemented. May be useful for optimization in future.
File size: 17.3 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	#include "gases.h"
31
32	!==================================================================
33
34	logical file_ok
35
36	integer n, nt, np, nh, ng, nw, m, i
37	integer L_NGAUSScheck, L_NPREFcheck, L_NTREFcheck, L_REFVARcheck
38
39	character(len=100) :: file_id
40	character(len=100) :: file_path
41
42	real*8 gasi8(L_NTREF,L_NPREF,L_REFVAR,L_NSPECTI,L_NGAUSS)
43	real*8 gasv8(L_NTREF,L_NPREF,L_REFVAR,L_NSPECTV,L_NGAUSS)
44
45	real*8 x, xi(4), yi(4), ans, p
46
47	integer ngas, igas
48
49	double precision testcont ! for continuum absorption initialisation
50
51	!=======================================================================
52	! Load variable species data, exit if we have wrong database
53	file_id='/corrk_data/' // corrkdir(1:LEN_TRIM(corrkdir)) // '/Q.dat'
54	file_path=datafile(1:LEN_TRIM(datafile))//file_id(1:LEN_TRIM(file_id))
55
56
57	! check that the file exists
58	inquire(FILE=file_path,EXIST=file_ok)
59	if(.not.file_ok) then
60	write(,)'The file ',file_path(1:LEN_TRIM(file_path))
61	write(,)'was not found by sugas_corrk.F90, exiting.'
62	call abort
63	endif
64
65	! check that database matches varactive toggle
66	open(111,file=file_path(1:LEN_TRIM(file_path)),form='formatted')
67	read(111,*) ngas
68
69	if(ngas.ne.ngasmx)then
70	print*,'Number of gases in radiative transfer data (',ngas,') does not', &
71	'match that in gases.def (',ngasmx,'), exiting.'
72	call abort
73	endif
74
75	if(ngas.eq.1 .and. (varactive.or.varfixed))then
76	print*,'You have varactive/fixed=.true. but the database [', &
77	corrkdir(1:LEN_TRIM(corrkdir)), &
78	'] has no variable species, exiting.'
79	call abort
80	elseif(ngas.eq.2 .and. (.not.varactive) .and. (.not.varfixed))then
81	print*,'You have varactive and varfixed=.false. and the database [', &
82	corrkdir(1:LEN_TRIM(corrkdir)), &
83	'] has a variable species.'
84	call abort
85	elseif(ngas.gt.4 .or. ngas.lt.1)then
86	print*,ngas,' species in database [', &
87	corrkdir(1:LEN_TRIM(corrkdir)), &
88	'], radiative code cannot handle this.'
89	call abort
90	endif
91
92	if(ngas.gt.3)then
93	print*,'ngas>3, EXPERIMENTAL!'
94	endif
95
96
97	do n=1,ngas
98	read(111,*) gastype(n)
99	print*,'Gas ',n,' is ',gastype(n)
100	enddo
101
102	! check the array size is correct, load the coefficients
103	open(111,file=file_path(1:LEN_TRIM(file_path)),form='formatted')
104	read(111,*) L_REFVARcheck
105	if(.not.(L_REFVARcheck.eq.L_REFVAR)) then
106	print*, L_REFVARcheck
107	print*, L_REFVAR
108	print*,'The size of your radiative transfer mixing ratio array does '
109	print*,'not match the value given in Q.dat, exiting.'
110	call abort
111	endif
112	read(111,*) wrefvar
113	close(111)
114
115	! Check that gastype and gnom match
116	do n=1,ngas
117	print*,'Gas ',n,' is ',gnom(n)
118	if(gnom(n).ne.gastype(n))then
119	print*,'Name of a gas in radiative transfer data (',gastype(n),') does not', &
120	'match that in gases.def (',gnom(n),'), exiting.'
121	call abort
122	endif
123	enddo
124	print*,'Confirmed gas match in radiative transfer and gases.def!'
125
126	! display the values
127	print*,'Variable gas mixing ratios:'
128	do n=1,L_REFVAR
129	!print*,n,'.',wrefvar(n),' kg/kg' ! pay attention!
130	print*,n,'.',wrefvar(n),' mol/mol'
131	end do
132	print*,''
133
134	!=======================================================================
135	! Set the weighting in g-space
136
137	file_id='/corrk_data/' // corrkdir(1:LEN_TRIM(corrkdir)) // '/g.dat'
138	file_path=datafile(1:LEN_TRIM(datafile))//file_id(1:LEN_TRIM(file_id))
139
140	! check that the file exists
141	inquire(FILE=file_path,EXIST=file_ok)
142	if(.not.file_ok) then
143	write(,)'The file ',file_path(1:LEN_TRIM(file_path))
144	write(,)'was not found by sugas_corrk.F90, exiting.'
145	call abort
146	endif
147
148	! check the array size is correct, load the coefficients
149	open(111,file=file_path(1:LEN_TRIM(file_path)),form='formatted')
150	read(111,*) L_NGAUSScheck
151	if(.not.(L_NGAUSScheck.eq.L_NGAUSS)) then
152	print*,'The size of your radiative transfer g-space array does '
153	print*,'not match the value given in g.dat, exiting.'
154	call abort
155	endif
156	read(111,*) gweight
157	close(111)
158
159	! display the values
160	print*,'Correlated-k g-space grid:'
161	do n=1,L_NGAUSS
162	print*,n,'.',gweight(n)
163	end do
164	print*,''
165
166	!=======================================================================
167	! Set the reference pressure and temperature arrays. These are
168	! the pressures and temperatures at which we have k-coefficients.
169
170	!-----------------------------------------------------------------------
171	! pressure
172
173	file_id='/corrk_data/' // corrkdir(1:LEN_TRIM(corrkdir)) // '/p.dat'
174	file_path=datafile(1:LEN_TRIM(datafile))//file_id(1:LEN_TRIM(file_id))
175
176	! check that the file exists
177	inquire(FILE=file_path,EXIST=file_ok)
178	if(.not.file_ok) then
179	write(,)'The file ',file_path(1:LEN_TRIM(file_path))
180	write(,)'was not found by sugas_corrk.F90, exiting.'
181	call abort
182	endif
183
184	! check the array size is correct, load the coefficients
185	open(111,file=file_path(1:LEN_TRIM(file_path)),form='formatted')
186	read(111,*) L_NPREFcheck
187	if(.not.(L_NPREFcheck.eq.L_NPREF)) then
188	print*,'The size of your radiative transfer pressure array does '
189	print*,'not match the value given in p.dat, exiting.'
190	call abort
191	endif
192	read(111,*) pgasref
193	close(111)
194
195	! display the values
196	print*,'Correlated-k pressure grid (mBar):'
197	do n=1,L_NPREF
198	print*,n,'. 1 x 10^',pgasref(n),' mBar'
199	end do
200	print*,''
201
202	! save the min / max matrix values
203	pgasmin = 10.0**pgasref(1)
204	pgasmax = 10.0**pgasref(L_NPREF)
205
206	! interpolate to finer grid
207	do n=1,L_NPREF-1
208	do m=1,5
209	pfgasref((n-1)5+m) = pgasref(n)+(m-1)0.2
210	end do
211	end do
212	pfgasref(L_PINT) = pgasref(L_NPREF)
213	print*,'Warning, pfgasref needs generalisation to uneven grids!!'
214
215	!-----------------------------------------------------------------------
216	! temperature
217
218	file_id='/corrk_data/' // corrkdir(1:LEN_TRIM(corrkdir)) // '/T.dat'
219	file_path=datafile(1:LEN_TRIM(datafile))//file_id(1:LEN_TRIM(file_id))
220
221	! check that the file exists
222	inquire(FILE=file_path,EXIST=file_ok)
223	if(.not.file_ok) then
224	write(,)'The file ',file_path(1:LEN_TRIM(file_path))
225	write(,)'was not found by sugas_corrk.F90, exiting.'
226	call abort
227	endif
228
229	! check the array size is correct, load the coefficients
230	open(111,file=file_path(1:LEN_TRIM(file_path)),form='formatted')
231	read(111,*) L_NTREFcheck
232	if(.not.(L_NTREFcheck.eq.L_NTREF)) then
233	print*,'The size of your radiative transfer temperature array does '
234	print*,'not match the value given in T.dat, exiting.'
235	call abort
236	endif
237	read(111,*) tgasref
238	close(111)
239
240	! display the values
241	print*,'Correlated-k temperature grid:'
242	do n=1,L_NTREF
243	print*,n,'.',tgasref(n),' K'
244	end do
245
246	! save the min / max matrix values
247	tgasmin = tgasref(1)
248	tgasmax = tgasref(L_NTREF)
249
250	!=======================================================================
251	! Get gaseous k-coefficients and interpolate onto finer pressure grid
252
253	! VISIBLE
254	if (callgasvis.and..not.corrkdir(1:LEN_TRIM(corrkdir)).eq.'null') then
255	file_id='/corrk_data/'//trim(adjustl(banddir))//'/corrk_gcm_VI.dat'
256	file_path=datafile(1:LEN_TRIM(datafile))//file_id(1:LEN_TRIM(file_id))
257
258	! check that the file exists
259	inquire(FILE=file_path,EXIST=file_ok)
260	if(.not.file_ok) then
261	write(,)'The file ',file_path(1:LEN_TRIM(file_path))
262	write(,)'was not found by sugas_corrk.F90.'
263	write(,)'Are you sure you have absorption data for these bands?'
264	call abort
265	endif
266
267	open(111,file=file_path(1:LEN_TRIM(file_path)),form='formatted')
268	read(111,*) gasv8
269	close(111)
270
271	else
272	print*,'Visible corrk gaseous absorption is set to zero.'
273	gasv8(:,:,:,:,:)=0.0
274	endif
275
276	! INFRA-RED
277	if (.not.corrkdir(1:LEN_TRIM(corrkdir)).eq.'null') then
278	file_id='/corrk_data/'//trim(adjustl(banddir))//'/corrk_gcm_IR.dat'
279	file_path=datafile(1:LEN_TRIM(datafile))//file_id(1:LEN_TRIM(file_id))
280
281	! check that the file exists
282	inquire(FILE=file_path,EXIST=file_ok)
283	if(.not.file_ok) then
284	write(,)'The file ',file_path(1:LEN_TRIM(file_path))
285	write(,)'was not found by sugas_corrk.F90.'
286	write(,)'Are you sure you have absorption data for these bands?'
287	call abort
288	endif
289
290	open(111,file=file_path(1:LEN_TRIM(file_path)),form='formatted')
291	read(111,*) gasi8
292	close(111)
293
294	do nw=1,L_NSPECTI
295	fzeroi(nw) = 0.
296	! do nt=1,L_NTREF
297	! do np=1,L_NPREF
298	! do nh=1,L_REFVAR
299	! do ng = 1,L_NGAUSS
300	! if(gasi8(nt,np,nh,nw,ng).lt.1.0e-25)then
301	! fzeroi(nw)=fzeroi(nw)+1.
302	! endif
303	! end do
304	! end do
305	! end do
306	! end do
307	! fzeroi(nw)=fzeroi(nw)/dble(L_NTREFL_NPREFL_REFVAR*L_NGAUSS)
308	end do
309
310	do nw=1,L_NSPECTV
311	fzerov(nw) = 0.
312	! do nt=1,L_NTREF
313	! do np=1,L_NPREF
314	! do nh=1,L_REFVAR
315	! do ng = 1,L_NGAUSS
316	! if(gasv8(nt,np,nh,nw,ng).lt.1.0e-25)then
317	! fzerov(nw)=fzerov(nw)+1.
318	! endif
319	! end do
320	! end do
321	! end do
322	! end do
323	! fzerov(nw)=fzerov(nw)/dble(L_NTREFL_NPREFL_REFVAR*L_NGAUSS)
324	end do
325
326
327	do nw=1,L_NSPECTV
328	fzerov(nw) = 0.
329	end do
330
331	else
332	print*,'Infrared corrk gaseous absorption is set to zero.'
333	gasi8(:,:,:,:,:)=0.0
334	endif
335
336	! Take log10 of the values - this is what we will interpolate.
337	! Smallest value is 1.0E-200.
338
339	do nt=1,L_NTREF
340	do np=1,L_NPREF
341	do nh=1,L_REFVAR
342	do ng = 1,L_NGAUSS
343
344	do nw=1,L_NSPECTV
345	if(gasv8(nt,np,nh,nw,ng).gt.1.0d-200) then
346	gasv8(nt,np,nh,nw,ng) = &
347	log10(gasv8(nt,np,nh,nw,ng))
348	else
349	gasv8(nt,np,nh,nw,ng) = -200.0
350	end if
351	end do
352
353	do nw=1,L_NSPECTI
354	if(gasi8(nt,np,nh,nw,ng).gt.1.0d-200) then
355	gasi8(nt,np,nh,nw,ng) = &
356	log10(gasi8(nt,np,nh,nw,ng))
357	else
358	gasi8(nt,np,nh,nw,ng) = -200.0
359	end if
360	end do
361
362	end do
363	end do
364	end do
365	end do
366
367	! Interpolate the values: first the longwave
368
369	do nt=1,L_NTREF
370	do nh=1,L_REFVAR
371	do nw=1,L_NSPECTI
372	do ng=1,L_NGAUSS
373
374	! First, the initial interval
375
376	n = 1
377	do m=1,5
378	x = pfgasref(m)
379	xi(1) = pgasref(n)
380	xi(2) = pgasref(n+1)
381	xi(3) = pgasref(n+2)
382	xi(4) = pgasref(n+3)
383	yi(1) = gasi8(nt,n,nh,nw,ng)
384	yi(2) = gasi8(nt,n+1,nh,nw,ng)
385	yi(3) = gasi8(nt,n+2,nh,nw,ng)
386	yi(4) = gasi8(nt,n+3,nh,nw,ng)
387	call lagrange(x,xi,yi,ans)
388	gasi(nt,m,nh,nw,ng) = 10.0**ans
389	end do
390
391	do n=2,L_NPREF-2
392	do m=1,5
393	i = (n-1)*5+m
394	x = pfgasref(i)
395	xi(1) = pgasref(n-1)
396	xi(2) = pgasref(n)
397	xi(3) = pgasref(n+1)
398	xi(4) = pgasref(n+2)
399	yi(1) = gasi8(nt,n-1,nh,nw,ng)
400	yi(2) = gasi8(nt,n,nh,nw,ng)
401	yi(3) = gasi8(nt,n+1,nh,nw,ng)
402	yi(4) = gasi8(nt,n+2,nh,nw,ng)
403	call lagrange(x,xi,yi,ans)
404	gasi(nt,i,nh,nw,ng) = 10.0**ans
405	end do
406	end do
407
408	! Now, get the last interval
409
410	n = L_NPREF-1
411	do m=1,5
412	i = (n-1)*5+m
413	x = pfgasref(i)
414	xi(1) = pgasref(n-2)
415	xi(2) = pgasref(n-1)
416	xi(3) = pgasref(n)
417	xi(4) = pgasref(n+1)
418	yi(1) = gasi8(nt,n-2,nh,nw,ng)
419	yi(2) = gasi8(nt,n-1,nh,nw,ng)
420	yi(3) = gasi8(nt,n,nh,nw,ng)
421	yi(4) = gasi8(nt,n+1,nh,nw,ng)
422	call lagrange(x,xi,yi,ans)
423	gasi(nt,i,nh,nw,ng) = 10.0**ans
424	end do
425
426	! Fill the last pressure point
427
428	gasi(nt,L_PINT,nh,nw,ng) = &
429	10.0**gasi8(nt,L_NPREF,nh,nw,ng)
430
431	end do
432	end do
433	end do
434	end do
435
436	! Interpolate the values: now the shortwave
437
438	do nt=1,L_NTREF
439	do nh=1,L_REFVAR
440	do nw=1,L_NSPECTV
441	do ng=1,L_NGAUSS
442
443	! First, the initial interval
444
445	n = 1
446	do m=1,5
447	x = pfgasref(m)
448	xi(1) = pgasref(n)
449	xi(2) = pgasref(n+1)
450	xi(3) = pgasref(n+2)
451	xi(4) = pgasref(n+3)
452	yi(1) = gasv8(nt,n,nh,nw,ng)
453	yi(2) = gasv8(nt,n+1,nh,nw,ng)
454	yi(3) = gasv8(nt,n+2,nh,nw,ng)
455	yi(4) = gasv8(nt,n+3,nh,nw,ng)
456	call lagrange(x,xi,yi,ans)
457	gasv(nt,m,nh,nw,ng) = 10.0**ans
458	end do
459
460	do n=2,L_NPREF-2
461	do m=1,5
462	i = (n-1)*5+m
463	x = pfgasref(i)
464	xi(1) = pgasref(n-1)
465	xi(2) = pgasref(n)
466	xi(3) = pgasref(n+1)
467	xi(4) = pgasref(n+2)
468	yi(1) = gasv8(nt,n-1,nh,nw,ng)
469	yi(2) = gasv8(nt,n,nh,nw,ng)
470	yi(3) = gasv8(nt,n+1,nh,nw,ng)
471	yi(4) = gasv8(nt,n+2,nh,nw,ng)
472	call lagrange(x,xi,yi,ans)
473	gasv(nt,i,nh,nw,ng) = 10.0**ans
474	end do
475	end do
476
477	! Now, get the last interval
478
479	n = L_NPREF-1
480	do m=1,5
481	i = (n-1)*5+m
482	x = pfgasref(i)
483	xi(1) = pgasref(n-2)
484	xi(2) = pgasref(n-1)
485	xi(3) = pgasref(n)
486	xi(4) = pgasref(n+1)
487	yi(1) = gasv8(nt,n-2,nh,nw,ng)
488	yi(2) = gasv8(nt,n-1,nh,nw,ng)
489	yi(3) = gasv8(nt,n,nh,nw,ng)
490	yi(4) = gasv8(nt,n+1,nh,nw,ng)
491	call lagrange(x,xi,yi,ans)
492	gasv(nt,i,nh,nw,ng) = 10.0**ans
493	end do
494
495	! Fill the last pressure point
496
497	gasv(nt,L_PINT,nh,nw,ng) = &
498	10.0**gasv8(nt,L_NPREF,nh,nw,ng)
499
500	end do
501	end do
502	end do
503	end do
504
505
506
507	do igas=1,ngasmx
508	if(gnom(igas).eq.'H2_')then
509	call interpolateH2H2(500.D+0,250.D+0,17500.D+0,testcont,.true.)
510	elseif(gnom(igas).eq.'H2O')then
511	call interpolateH2Ocont(990.D+0,296.D+0,683.2D+0*2,0.D+0,testcont,.true.)
512	endif
513	enddo
514
515	return
516	end subroutine sugas_corrk

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

Download in other formats:

Original Format