Context Navigation

sugas_corrk.F90 @ 1371

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

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

Download in other formats:

Original Format