Context Navigation

sugas_corrk.F90 @ 937

Last change on this file since 937 was 879, checked in by jleconte, 12 years ago
forgot two files in previous commit
File size: 22.6 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	! Added double gray case by Jeremy Leconte (2012)
15	! New HITRAN continuum data section by RW (2012)
16	!
17	! Summary
18	! -------
19	!
20	!==================================================================
21
22	use radinc_h
23	use radcommon_h, only : pgasref,pfgasref,pgasmin,pgasmax
24	use radcommon_h, only : tgasref,tgasmin,tgasmax
25	use radcommon_h, only : gasv,gasi,FZEROI,FZEROV,gweight
26	use radcommon_h, only : wrefvar,WNOI,WNOV
27	use datafile_mod, only: datadir
28
29	use gases_h
30	use ioipsl_getincom
31	implicit none
32
33	#include "callkeys.h"
34	#include "comcstfi.h"
35
36	!==================================================================
37
38	logical file_ok
39
40	integer n, nt, np, nh, ng, nw, m, i
41	integer L_NGAUSScheck
42
43	character(len=200) :: file_id
44	character(len=500) :: file_path
45
46	! ALLOCATABLE ARRAYS -- AS 12/2011
47	REAL*8, DIMENSION(:,:,:,:,:), ALLOCATABLE :: gasi8, gasv8
48	character*20,allocatable,DIMENSION(:) :: gastype ! for check with gnom
49
50	real*8 x, xi(4), yi(4), ans, p
51	! For gray case (JL12)
52	real kappa_IR, kappa_VI, IR_VI_wnlimit, nVI_limit,nIR_limit
53
54	integer ngas, igas, jgas
55
56	double precision testcont ! for continuum absorption initialisation
57
58	integer :: dummy
59
60	!=======================================================================
61	! Load variable species data, exit if we have wrong database
62	file_id='/corrk_data/' // TRIM(corrkdir) // '/Q.dat'
63	file_path=TRIM(datadir)//TRIM(file_id)
64
65	! check that the file exists
66	inquire(FILE=file_path,EXIST=file_ok)
67	if(.not.file_ok) then
68	write(,)'The file ',TRIM(file_path)
69	write(,)'was not found by sugas_corrk.F90, exiting.'
70	write(,)'Check that your path to datagcm:',trim(datadir)
71	write(,)' is correct. You can change it in callphys.def with:'
72	write(,)' datadir = /absolute/path/to/datagcm'
73	write(,)'Also check that the corrkdir you chose in callphys.def exists.'
74	call abort
75	endif
76
77	! check that database matches varactive toggle
78	open(111,file=TRIM(file_path),form='formatted')
79	read(111,*) ngas
80
81	if(ngas.ne.ngasmx)then
82	print*,'Number of gases in radiative transfer data (',ngas,') does not', &
83	'match that in gases.def (',ngasmx,'), exiting.'
84	call abort
85	endif
86
87	if(ngas.eq.1 .and. (varactive.or.varfixed))then
88	print*,'You have varactive/fixed=.true. but the database [', &
89	corrkdir(1:LEN_TRIM(corrkdir)), &
90	'] has no variable species, exiting.'
91	call abort
92	elseif(ngas.gt.5 .or. ngas.lt.1)then
93	print*,ngas,' species in database [', &
94	corrkdir(1:LEN_TRIM(corrkdir)), &
95	'], radiative code cannot handle this.'
96	call abort
97	endif
98
99	! dynamically allocate gastype and read from Q.dat
100	IF ( .NOT. ALLOCATED( gastype ) ) ALLOCATE( gastype( ngas ) )
101
102	do igas=1,ngas
103	read(111,*) gastype(igas)
104	print*,'Gas ',igas,' is ',gastype(igas)
105	enddo
106
107	! get array size, load the coefficients
108	open(111,file=TRIM(file_path),form='formatted')
109	read(111,*) L_REFVAR
110	IF( .NOT. ALLOCATED( wrefvar ) ) ALLOCATE( WREFVAR(L_REFVAR) )
111	read(111,*) wrefvar
112	close(111)
113
114	if(L_REFVAR.gt.1 .and. (.not.varactive) .and. (.not.varfixed))then
115	print*,'You have varactive and varfixed=.false. and the database [', &
116	corrkdir(1:LEN_TRIM(corrkdir)), &
117	'] has a variable species.'
118	call abort
119	endif
120
121	! Check that gastype and gnom match
122	do igas=1,ngas
123	print*,'Gas ',igas,' is ',trim(gnom(igas))
124	if (trim(gnom(igas)).ne.trim(gastype(igas))) then
125	print*,'Name of a gas in radiative transfer data (',trim(gastype(igas)),') does not ', &
126	'match that in gases.def (',trim(gnom(igas)),'), exiting. You should compare ', &
127	'gases.def with Q.dat in your radiative transfer directory.'
128	call abort
129	endif
130	enddo
131	print*,'Confirmed gas match in radiative transfer and gases.def!'
132
133	! display the values
134	print*,'Variable gas volume mixing ratios:'
135	do n=1,L_REFVAR
136	!print*,n,'.',wrefvar(n),' kg/kg' ! pay attention!
137	print*,n,'.',wrefvar(n),' mol/mol'
138	end do
139	print*,''
140
141	!=======================================================================
142	! Set the weighting in g-space
143
144	file_id='/corrk_data/' // TRIM(corrkdir) // '/g.dat'
145	file_path=TRIM(datadir)//TRIM(file_id)
146
147	! check that the file exists
148	inquire(FILE=file_path,EXIST=file_ok)
149	if(.not.file_ok) then
150	write(,)'The file ',TRIM(file_path)
151	write(,)'was not found by sugas_corrk.F90, exiting.'
152	write(,)'Check that your path to datagcm:',trim(datadir)
153	write(,)' is correct. You can change it in callphys.def with:'
154	write(,)' datadir = /absolute/path/to/datagcm'
155	write(,)'Also check that the corrkdir you chose in callphys.def exists.'
156	call abort
157	endif
158
159	! check the array size is correct, load the coefficients
160	open(111,file=TRIM(file_path),form='formatted')
161	read(111,*) L_NGAUSScheck
162	if(.not.(L_NGAUSScheck.eq.L_NGAUSS)) then
163	print*,'The size of your radiative transfer g-space array does '
164	print*,'not match the value given in g.dat, exiting.'
165	call abort
166	endif
167	read(111,*) gweight
168	close(111)
169
170	! display the values
171	print*,'Correlated-k g-space grid:'
172	do n=1,L_NGAUSS
173	print*,n,'.',gweight(n)
174	end do
175	print*,''
176
177	!=======================================================================
178	! Set the reference pressure and temperature arrays. These are
179	! the pressures and temperatures at which we have k-coefficients.
180
181	!-----------------------------------------------------------------------
182	! pressure
183
184	file_id='/corrk_data/' // TRIM(corrkdir) // '/p.dat'
185	file_path=TRIM(datadir)//TRIM(file_id)
186
187	! check that the file exists
188	inquire(FILE=file_path,EXIST=file_ok)
189	if(.not.file_ok) then
190	write(,)'The file ',TRIM(file_path)
191	write(,)'was not found by sugas_corrk.F90, exiting.'
192	write(,)'Check that your path to datagcm:',trim(datadir)
193	write(,)' is correct. You can change it in callphys.def with:'
194	write(,)' datadir = /absolute/path/to/datagcm'
195	write(,)'Also check that the corrkdir you chose in callphys.def exists.'
196	call abort
197	endif
198
199	! get array size, load the coefficients
200	open(111,file=TRIM(file_path),form='formatted')
201	read(111,*) L_NPREF
202	IF( .NOT. ALLOCATED( pgasref ) ) ALLOCATE( PGASREF(L_NPREF) )
203	read(111,*) pgasref
204	close(111)
205	L_PINT = (L_NPREF-1)*5+1
206	IF( .NOT. ALLOCATED( pfgasref ) ) ALLOCATE( PFGASREF(L_PINT) )
207
208	! display the values
209	print*,'Correlated-k pressure grid (mBar):'
210	do n=1,L_NPREF
211	print*,n,'. 1 x 10^',pgasref(n),' mBar'
212	end do
213	print*,''
214
215	! save the min / max matrix values
216	pgasmin = 10.0**pgasref(1)
217	pgasmax = 10.0**pgasref(L_NPREF)
218
219	! interpolate to finer grid
220	do n=1,L_NPREF-1
221	do m=1,5
222	pfgasref((n-1)5+m) = pgasref(n)+(m-1)0.2
223	end do
224	end do
225	pfgasref(L_PINT) = pgasref(L_NPREF)
226	! Warning, pfgasref needs generalisation to uneven grids!
227
228	!-----------------------------------------------------------------------
229	! temperature
230
231	file_id='/corrk_data/' // TRIM(corrkdir) // '/T.dat'
232	file_path=TRIM(datadir)//TRIM(file_id)
233
234	! check that the file exists
235	inquire(FILE=file_path,EXIST=file_ok)
236	if(.not.file_ok) then
237	write(,)'The file ',TRIM(file_path)
238	write(,)'was not found by sugas_corrk.F90, exiting.'
239	write(,)'Check that your path to datagcm:',trim(datadir)
240	write(,)' is correct. You can change it in callphys.def with:'
241	write(,)' datadir = /absolute/path/to/datagcm'
242	write(,)'Also check that the corrkdir you chose in callphys.def exists.'
243	call abort
244	endif
245
246	! get array size, load the coefficients
247	open(111,file=TRIM(file_path),form='formatted')
248	read(111,*) L_NTREF
249	IF( .NOT. ALLOCATED( tgasref ) ) ALLOCATE( TGASREF(L_NTREF) )
250	read(111,*) tgasref
251	close(111)
252
253	! display the values
254	print*,'Correlated-k temperature grid:'
255	do n=1,L_NTREF
256	print*,n,'.',tgasref(n),' K'
257	end do
258
259	! save the min / max matrix values
260	tgasmin = tgasref(1)
261	tgasmax = tgasref(L_NTREF)
262
263	!-----------------------------------------------------------------------
264	! allocate the multidimensional arrays in radcommon_h
265
266	IF( .NOT. ALLOCATED( gasi8 ) ) ALLOCATE( gasi8(L_NTREF,L_NPREF,L_REFVAR,L_NSPECTI,L_NGAUSS) )
267	IF( .NOT. ALLOCATED( gasv8 ) ) ALLOCATE( gasv8(L_NTREF,L_NPREF,L_REFVAR,L_NSPECTV,L_NGAUSS) )
268	IF( .NOT. ALLOCATED( gasi ) ) ALLOCATE( gasi(L_NTREF,L_PINT,L_REFVAR,L_NSPECTI,L_NGAUSS) )
269	IF( .NOT. ALLOCATED( gasv ) ) ALLOCATE( gasv(L_NTREF,L_PINT,L_REFVAR,L_NSPECTV,L_NGAUSS) )
270
271	! display the values
272	print*,''
273	print*,'Correlated-k matrix size:'
274	print*,'[',L_NTREF,',',L_NPREF,',',L_REFVAR,',',L_NGAUSS,']'
275
276	!=======================================================================
277	! Get gaseous k-coefficients and interpolate onto finer pressure grid
278
279	if (graybody) then
280	! constant absorption coefficient in visible
281	write(,)"graybody: constant absorption coefficient in visible:"
282	kappa_VI=-100000.
283	call getin("kappa_VI",kappa_VI)
284	write(,)" kappa_VI = ",kappa_VI
285	kappa_VI=kappa_VI1.e4 mugaz * 1.672621e-27 ! conversion from m^2/kg to cm^2/molecule
286
287	! constant absorption coefficient in IR
288	write(,)"graybody: constant absorption coefficient in InfraRed:"
289	kappa_IR=-100000.
290	call getin("kappa_IR",kappa_IR)
291	write(,)" kappa_IR = ",kappa_IR
292	kappa_IR=kappa_IR1.e4 mugaz * 1.672621e-27 ! conversion from m^2/kg to cm^2/molecule
293
294	! wavelength used to separate IR from VI
295	IR_VI_wnlimit=3000.
296	write(,)"graybody: Visible / Infrared separation set at",10000./IR_VI_wnlimit,"um"
297
298	nVI_limit=L_NSPECTV
299	Do nw=1,L_NSPECTV
300	if (WNOV(nw).gt.IR_VI_wnlimit) then
301	nVI_limit=nw-1
302	exit
303	endif
304	End do
305	nIR_limit=L_NSPECTI
306	Do nw=1,L_NSPECTI
307	if (WNOI(nw).gt.IR_VI_wnlimit) then
308	nIR_limit=nw-1
309	exit
310	endif
311	End do
312	write(,)"graybody: Visible / Infrared separation set at band ",nIR_limit,nVI_limit
313
314	End if
315
316
317	! VISIBLE
318	if (callgasvis.and.(.not.TRIM(corrkdir).eq.'null').and.(.not.TRIM(corrkdir).eq.'null_LowTeffStar').and.(.not.graybody)) then
319	file_id='/corrk_data/'//trim(adjustl(banddir))//'/corrk_gcm_VI.dat'
320	file_path=TRIM(datadir)//TRIM(file_id)
321
322	! check that the file exists
323	inquire(FILE=file_path,EXIST=file_ok)
324	if(.not.file_ok) then
325	write(,)'The file ',TRIM(file_path)
326	write(,)'was not found by sugas_corrk.F90.'
327	write(,)'Are you sure you have absorption data for these bands?'
328	call abort
329	endif
330
331	open(111,file=TRIM(file_path),form='formatted')
332	read(111,*) gasv8
333	close(111)
334
335	else if (callgasvis.and.graybody) then
336	if(nVI_limit.eq.0) then
337	gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTV,1:L_NGAUSS)=kappa_VI
338	else if (nVI_limit.eq.L_NSPECTV) then
339	gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTV,1:L_NGAUSS)=kappa_IR
340	else
341	gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:nVI_limit,1:L_NGAUSS)=kappa_IR
342	gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,nVI_limit+1:L_NSPECTV,1:L_NGAUSS)=kappa_VI
343	end if
344	else
345	print*,'Visible corrk gaseous absorption is set to zero.'
346	gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTV,1:L_NGAUSS)=0.0
347	endif
348
349	! INFRA-RED
350	if ((.not.TRIM(corrkdir).eq.'null').and.(.not.TRIM(corrkdir).eq.'null_LowTeffStar').and.(.not.graybody)) then
351	file_id='/corrk_data/'//trim(adjustl(banddir))//'/corrk_gcm_IR.dat'
352	file_path=TRIM(datadir)//TRIM(file_id)
353
354	! check that the file exists
355	inquire(FILE=file_path,EXIST=file_ok)
356	if(.not.file_ok) then
357	write(,)'The file ',TRIM(file_path)
358	write(,)'was not found by sugas_corrk.F90.'
359	write(,)'Are you sure you have absorption data for these bands?'
360	call abort
361	endif
362
363	open(111,file=TRIM(file_path),form='formatted')
364	read(111,*) gasi8
365	close(111)
366
367	! 'fzero' is a currently unused feature that allows optimisation
368	! of the radiative transfer by neglecting bands where absorption
369	! is close to zero. As it could be useful in the future, this
370	! section of the code has been kept commented and not erased.
371	! RW 7/3/12.
372
373	do nw=1,L_NSPECTI
374	fzeroi(nw) = 0.
375	! do nt=1,L_NTREF
376	! do np=1,L_NPREF
377	! do nh=1,L_REFVAR
378	! do ng = 1,L_NGAUSS
379	! if(gasi8(nt,np,nh,nw,ng).lt.1.0e-25)then
380	! fzeroi(nw)=fzeroi(nw)+1.
381	! endif
382	! end do
383	! end do
384	! end do
385	! end do
386	! fzeroi(nw)=fzeroi(nw)/dble(L_NTREFL_NPREFL_REFVAR*L_NGAUSS)
387	end do
388
389	do nw=1,L_NSPECTV
390	fzerov(nw) = 0.
391	! do nt=1,L_NTREF
392	! do np=1,L_NPREF
393	! do nh=1,L_REFVAR
394	! do ng = 1,L_NGAUSS
395	! if(gasv8(nt,np,nh,nw,ng).lt.1.0e-25)then
396	! fzerov(nw)=fzerov(nw)+1.
397	! endif
398	! end do
399	! end do
400	! end do
401	! end do
402	! fzerov(nw)=fzerov(nw)/dble(L_NTREFL_NPREFL_REFVAR*L_NGAUSS)
403	end do
404
405	else if (graybody) then
406	if(nIR_limit.eq.0) then
407	gasi8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTI,1:L_NGAUSS)=kappa_VI
408	else if (nIR_limit.eq.L_NSPECTI) then
409	gasi8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTI,1:L_NGAUSS)=kappa_IR
410	else
411	gasi8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:nIR_limit,1:L_NGAUSS)=kappa_IR
412	gasi8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,nIR_limit+1:L_NSPECTI,1:L_NGAUSS)=kappa_VI
413	end if
414	else
415	print*,'Infrared corrk gaseous absorption is set to zero.'
416	gasi8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTI,1:L_NGAUSS)=0.0
417	endif
418
419
420	! Take log10 of the values - this is what we will interpolate.
421	! Smallest value is 1.0E-200.
422
423	do nt=1,L_NTREF
424	do np=1,L_NPREF
425	do nh=1,L_REFVAR
426	do ng = 1,L_NGAUSS
427
428	do nw=1,L_NSPECTV
429	if(gasv8(nt,np,nh,nw,ng).gt.1.0d-200) then
430	gasv8(nt,np,nh,nw,ng) = log10(gasv8(nt,np,nh,nw,ng))
431	else
432	gasv8(nt,np,nh,nw,ng) = -200.0
433	end if
434	end do
435
436	do nw=1,L_NSPECTI
437	if(gasi8(nt,np,nh,nw,ng).gt.1.0d-200) then
438	gasi8(nt,np,nh,nw,ng) = log10(gasi8(nt,np,nh,nw,ng))
439	else
440	gasi8(nt,np,nh,nw,ng) = -200.0
441	end if
442	end do
443
444	end do
445	end do
446	end do
447	end do
448
449	! Interpolate the values: first the longwave
450
451	do nt=1,L_NTREF
452	do nh=1,L_REFVAR
453	do nw=1,L_NSPECTI
454	do ng=1,L_NGAUSS
455
456	! First, the initial interval
457
458	n = 1
459	do m=1,5
460	x = pfgasref(m)
461	xi(1) = pgasref(n)
462	xi(2) = pgasref(n+1)
463	xi(3) = pgasref(n+2)
464	xi(4) = pgasref(n+3)
465	yi(1) = gasi8(nt,n,nh,nw,ng)
466	yi(2) = gasi8(nt,n+1,nh,nw,ng)
467	yi(3) = gasi8(nt,n+2,nh,nw,ng)
468	yi(4) = gasi8(nt,n+3,nh,nw,ng)
469	call lagrange(x,xi,yi,ans)
470	gasi(nt,m,nh,nw,ng) = 10.0**ans
471	end do
472
473	do n=2,L_NPREF-2
474	do m=1,5
475	i = (n-1)*5+m
476	x = pfgasref(i)
477	xi(1) = pgasref(n-1)
478	xi(2) = pgasref(n)
479	xi(3) = pgasref(n+1)
480	xi(4) = pgasref(n+2)
481	yi(1) = gasi8(nt,n-1,nh,nw,ng)
482	yi(2) = gasi8(nt,n,nh,nw,ng)
483	yi(3) = gasi8(nt,n+1,nh,nw,ng)
484	yi(4) = gasi8(nt,n+2,nh,nw,ng)
485	call lagrange(x,xi,yi,ans)
486	gasi(nt,i,nh,nw,ng) = 10.0**ans
487	end do
488	end do
489
490	! Now, get the last interval
491
492	n = L_NPREF-1
493	do m=1,5
494	i = (n-1)*5+m
495	x = pfgasref(i)
496	xi(1) = pgasref(n-2)
497	xi(2) = pgasref(n-1)
498	xi(3) = pgasref(n)
499	xi(4) = pgasref(n+1)
500	yi(1) = gasi8(nt,n-2,nh,nw,ng)
501	yi(2) = gasi8(nt,n-1,nh,nw,ng)
502	yi(3) = gasi8(nt,n,nh,nw,ng)
503	yi(4) = gasi8(nt,n+1,nh,nw,ng)
504	call lagrange(x,xi,yi,ans)
505	gasi(nt,i,nh,nw,ng) = 10.0**ans
506	end do
507
508	! Fill the last pressure point
509
510	gasi(nt,L_PINT,nh,nw,ng) = &
511	10.0**gasi8(nt,L_NPREF,nh,nw,ng)
512
513	end do
514	end do
515	end do
516	end do
517
518	! Interpolate the values: now the shortwave
519
520	do nt=1,L_NTREF
521	do nh=1,L_REFVAR
522	do nw=1,L_NSPECTV
523	do ng=1,L_NGAUSS
524
525	! First, the initial interval
526
527	n = 1
528	do m=1,5
529	x = pfgasref(m)
530	xi(1) = pgasref(n)
531	xi(2) = pgasref(n+1)
532	xi(3) = pgasref(n+2)
533	xi(4) = pgasref(n+3)
534	yi(1) = gasv8(nt,n,nh,nw,ng)
535	yi(2) = gasv8(nt,n+1,nh,nw,ng)
536	yi(3) = gasv8(nt,n+2,nh,nw,ng)
537	yi(4) = gasv8(nt,n+3,nh,nw,ng)
538	call lagrange(x,xi,yi,ans)
539	gasv(nt,m,nh,nw,ng) = 10.0**ans
540	end do
541
542	do n=2,L_NPREF-2
543	do m=1,5
544	i = (n-1)*5+m
545	x = pfgasref(i)
546	xi(1) = pgasref(n-1)
547	xi(2) = pgasref(n)
548	xi(3) = pgasref(n+1)
549	xi(4) = pgasref(n+2)
550	yi(1) = gasv8(nt,n-1,nh,nw,ng)
551	yi(2) = gasv8(nt,n,nh,nw,ng)
552	yi(3) = gasv8(nt,n+1,nh,nw,ng)
553	yi(4) = gasv8(nt,n+2,nh,nw,ng)
554	call lagrange(x,xi,yi,ans)
555	gasv(nt,i,nh,nw,ng) = 10.0**ans
556	end do
557	end do
558
559	! Now, get the last interval
560
561	n = L_NPREF-1
562	do m=1,5
563	i = (n-1)*5+m
564	x = pfgasref(i)
565	xi(1) = pgasref(n-2)
566	xi(2) = pgasref(n-1)
567	xi(3) = pgasref(n)
568	xi(4) = pgasref(n+1)
569	yi(1) = gasv8(nt,n-2,nh,nw,ng)
570	yi(2) = gasv8(nt,n-1,nh,nw,ng)
571	yi(3) = gasv8(nt,n,nh,nw,ng)
572	yi(4) = gasv8(nt,n+1,nh,nw,ng)
573	call lagrange(x,xi,yi,ans)
574	gasv(nt,i,nh,nw,ng) = 10.0**ans
575	end do
576
577	! Fill the last pressure point
578
579	gasv(nt,L_PINT,nh,nw,ng) = &
580	10.0**gasv8(nt,L_NPREF,nh,nw,ng)
581
582	end do
583	end do
584	end do
585	end do
586
587
588	!=======================================================================
589	! Initialise the continuum absorption data
590
591	do igas=1,ngasmx
592
593	if (igas .eq. igas_N2) then
594
595	dummy = -9999
596	call interpolateN2N2(100.D+0,250.D+0,17500.D+0,testcont,.true.,dummy)
597
598	elseif (igas .eq. igas_H2) then
599
600	! first do self-induced absorption
601	dummy = -9999
602	call interpolateH2H2(500.D+0,250.D+0,17500.D+0,testcont,.true.,dummy)
603	! then cross-interactions with other gases
604	do jgas=1,ngasmx
605	if (jgas .eq. igas_N2) then
606	dummy = -9999
607	call interpolateN2H2(592.D+0,278.15D+0,200000.D+0,10000.D+0,testcont,.true.,dummy)
608	elseif (jgas .eq. igas_He) then
609	dummy = -9999
610	call interpolateH2He(500.D+0,250.D+0,200000.D+0,10000.D+0,testcont,.true.,dummy)
611	endif
612	enddo
613
614	elseif (igas .eq. igas_H2O) then
615
616	! H2O is special
617	if(H2Ocont_simple)then
618	call interpolateH2Ocont_PPC(990.D+0,296.D+0,683.2D+0*2,0.D+0,testcont,.true.)
619	else
620	dummy = -9999
621	call interpolateH2Ocont_CKD(990.D+0,296.D+0,683.2D+0*2,0.D+0,testcont,.true.,dummy)
622	endif
623
624	endif
625
626	enddo
627
628	print*,'----------------------------------------------------'
629	print*,'And that`s all we have. It`s possible that other'
630	print*,'continuum absorption may be present, but if it is we'
631	print*,'don`t yet have data for it...'
632	print*,''
633
634	! Deallocate local arrays
635	IF( ALLOCATED( gasi8 ) ) DEALLOCATE( gasi8 )
636	IF( ALLOCATED( gasv8 ) ) DEALLOCATE( gasv8 )
637	IF( ALLOCATED( pgasref ) ) DEALLOCATE( pgasref )
638	IF( ALLOCATED( gastype ) ) DEALLOCATE( gastype )
639
640	return
641	end subroutine sugas_corrk

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

Download in other formats:

Original Format