Context Navigation

sugas_corrk.F90 @ 1150

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

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

Download in other formats:

Original Format