Context Navigation

sugas_corrk.F90 @ 1980

Last change on this file since 1980 was 1521, checked in by emillour, 9 years ago

All GCMs: Updates to make planetary codes (+Earth) setups converge.

Made a "phy_common" directory in libf, to contain routines common (wrt structural nature of underlying code/grid) to all LMDZ-related physics packages.
moved all "mod_phys_*" and "mod_grid_phy_lmdz" files from dynlonlat_phylonlat to "phy_common"
moved "ioipsl_getincom_p.F90 from "misc" to "phy_common" and modified it to match Earth GCM version and renamed it ioipsl_getin_p_mod.F90
added an "abort_physics" (as in Earth GCM) in "phy_common"
added a "print_control_mod.F90 (as in Earth GCM) in phy_common
made similar changes in LMDZ.GENERIC and LMDZ.MARS

File size: 24.1 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	use comcstfi_mod, only: mugaz
29	use gases_h
30	use ioipsl_getin_p_mod, only: getin_p
31	use callkeys_mod, only: varactive,varfixed,graybody,callgasvis,&
32	continuum,H2Ocont_simple
33	implicit none
34
35	!==================================================================
36
37	logical file_ok
38
39	integer n, nt, np, nh, ng, nw, m, i
40	integer L_NGAUSScheck
41
42	character(len=200) :: file_id
43	character(len=500) :: file_path
44
45	! ALLOCATABLE ARRAYS -- AS 12/2011
46	REAL*8, DIMENSION(:,:,:,:,:), ALLOCATABLE,SAVE :: gasi8, gasv8 !read by master
47	character*20,allocatable,DIMENSION(:),SAVE :: gastype ! for check with gnom, read by master
48
49	real*8 x, xi(4), yi(4), ans, p
50	! For gray case (JL12)
51	real kappa_IR, kappa_VI, IR_VI_wnlimit
52	integer 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	!$OMP MASTER
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, adapted to uneven grids
221	do n=1,L_NPREF-1
222	do m=1,5
223	pfgasref((n-1)5+m) = pgasref(n)+(m-1)(pgasref(n+1) - pgasref(n))/5.
224	end do
225	end do
226	pfgasref(L_PINT) = pgasref(L_NPREF)
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	IF( .NOT. ALLOCATED( gasi8 ) ) ALLOCATE( gasi8(L_NTREF,L_NPREF,L_REFVAR,L_NSPECTI,L_NGAUSS) )
264	IF( .NOT. ALLOCATED( gasv8 ) ) ALLOCATE( gasv8(L_NTREF,L_NPREF,L_REFVAR,L_NSPECTV,L_NGAUSS) )
265	!$OMP END MASTER
266	!$OMP BARRIER
267
268	!-----------------------------------------------------------------------
269	! allocate the multidimensional arrays in radcommon_h
270	IF( .NOT. ALLOCATED( gasi ) ) ALLOCATE( gasi(L_NTREF,L_PINT,L_REFVAR,L_NSPECTI,L_NGAUSS) )
271	IF( .NOT. ALLOCATED( gasv ) ) ALLOCATE( gasv(L_NTREF,L_PINT,L_REFVAR,L_NSPECTV,L_NGAUSS) )
272
273	! display the values
274	print*,''
275	print*,'Correlated-k matrix size:'
276	print*,'[',L_NTREF,',',L_NPREF,',',L_REFVAR,',',L_NGAUSS,']'
277
278	!=======================================================================
279	! Get gaseous k-coefficients and interpolate onto finer pressure grid
280
281
282	! wavelength used to separate IR from VI in graybody. We will need that anyway
283	IR_VI_wnlimit=3000.
284	write(,)"graybody: Visible / Infrared separation set at",10000./IR_VI_wnlimit,"um"
285
286	nVI_limit=0
287	Do nw=1,L_NSPECTV
288	if ((WNOV(nw).gt.IR_VI_wnlimit).and.(L_NSPECTV.gt.1)) then
289	nVI_limit=nw-1
290	exit
291	endif
292	End do
293	nIR_limit=L_NSPECTI
294	Do nw=1,L_NSPECTI
295	if ((WNOI(nw).gt.IR_VI_wnlimit).and.(L_NSPECTI.gt.1)) then
296	nIR_limit=nw-1
297	exit
298	endif
299	End do
300
301	if (graybody) then
302	! constant absorption coefficient in visible
303	write(,)"graybody: constant absorption coefficient in visible:"
304	kappa_VI=-100000.
305	call getin_p("kappa_VI",kappa_VI)
306	write(,)" kappa_VI = ",kappa_VI
307	kappa_VI=kappa_VI1.e4 mugaz * 1.672621e-27 ! conversion from m^2/kg to cm^2/molecule
308
309	! constant absorption coefficient in IR
310	write(,)"graybody: constant absorption coefficient in InfraRed:"
311	kappa_IR=-100000.
312	call getin_p("kappa_IR",kappa_IR)
313	write(,)" kappa_IR = ",kappa_IR
314	kappa_IR=kappa_IR1.e4 mugaz * 1.672621e-27 ! conversion from m^2/kg to cm^2/molecule
315
316	write(,)"graybody: Visible / Infrared separation set at band: IR=",nIR_limit,", VI=",nVI_limit
317
318	Else
319	kappa_VI=1.e-30
320	kappa_IR=1.e-30
321	End if
322
323	!$OMP MASTER
324	! print*,corrkdir(1:4)
325	! VISIBLE
326	if (callgasvis) then
327	if ((corrkdir(1:4).eq.'null'))then !(TRIM(corrkdir).eq.'null_LowTeffStar')) then
328	gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTV,1:L_NGAUSS)=0.0
329	print*,'using no corrk data'
330	print*,'Visible corrk gaseous absorption is set to zero if graybody=F'
331	else
332	file_id='/corrk_data/'//trim(adjustl(banddir))//'/corrk_gcm_VI.dat'
333	file_path=TRIM(datadir)//TRIM(file_id)
334
335	! check that the file exists
336	inquire(FILE=file_path,EXIST=file_ok)
337	if(.not.file_ok) then
338	write(,)'The file ',TRIM(file_path)
339	write(,)'was not found by sugas_corrk.F90.'
340	write(,)'Are you sure you have absorption data for these bands?'
341	call abort
342	endif
343
344	open(111,file=TRIM(file_path),form='formatted')
345	read(111,*) gasv8
346	close(111)
347	end if
348
349	if(nVI_limit.eq.0) then
350	gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTV,1:L_NGAUSS)= &
351	gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTV,1:L_NGAUSS)+kappa_VI
352	else if (nVI_limit.eq.L_NSPECTV) then
353	gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTV,1:L_NGAUSS)= &
354	gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTV,1:L_NGAUSS)+kappa_IR
355	else
356	gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:nVI_limit,1:L_NGAUSS)= &
357	gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:nVI_limit,1:L_NGAUSS)+kappa_IR
358	gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,nVI_limit+1:L_NSPECTV,1:L_NGAUSS)= &
359	gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,nVI_limit+1:L_NSPECTV,1:L_NGAUSS)+kappa_VI
360	end if
361	else
362	print*,'Visible corrk gaseous absorption is set to zero.'
363	gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTV,1:L_NGAUSS)=0.0
364	endif
365	!$OMP END MASTER
366	!$OMP BARRIER
367
368	! INFRA-RED
369	if ((corrkdir(1:4).eq.'null'))then !.or.(TRIM(corrkdir).eq.'null_LowTeffStar')) then
370	print*,'Infrared corrk gaseous absorption is set to zero if graybody=F'
371	!$OMP MASTER
372	gasi8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTI,1:L_NGAUSS)=0.0
373	!$OMP END MASTER
374	!$OMP BARRIER
375	else
376	file_id='/corrk_data/'//trim(adjustl(banddir))//'/corrk_gcm_IR.dat'
377	file_path=TRIM(datadir)//TRIM(file_id)
378
379	! check that the file exists
380	inquire(FILE=file_path,EXIST=file_ok)
381	if(.not.file_ok) then
382	write(,)'The file ',TRIM(file_path)
383	write(,)'was not found by sugas_corrk.F90.'
384	write(,)'Are you sure you have absorption data for these bands?'
385	call abort
386	endif
387
388	!$OMP MASTER
389	open(111,file=TRIM(file_path),form='formatted')
390	read(111,*) gasi8
391	close(111)
392	!$OMP END MASTER
393	!$OMP BARRIER
394
395	! 'fzero' is a currently unused feature that allows optimisation
396	! of the radiative transfer by neglecting bands where absorption
397	! is close to zero. As it could be useful in the future, this
398	! section of the code has been kept commented and not erased.
399	! RW 7/3/12.
400
401	do nw=1,L_NSPECTI
402	fzeroi(nw) = 0.d0
403	! do nt=1,L_NTREF
404	! do np=1,L_NPREF
405	! do nh=1,L_REFVAR
406	! do ng = 1,L_NGAUSS
407	! if(gasi8(nt,np,nh,nw,ng).lt.1.0e-25)then
408	! fzeroi(nw)=fzeroi(nw)+1.d0
409	! endif
410	! end do
411	! end do
412	! end do
413	! end do
414	! fzeroi(nw)=fzeroi(nw)/dble(L_NTREFL_NPREFL_REFVAR*L_NGAUSS)
415	end do
416
417	do nw=1,L_NSPECTV
418	fzerov(nw) = 0.d0
419	! do nt=1,L_NTREF
420	! do np=1,L_NPREF
421	! do nh=1,L_REFVAR
422	! do ng = 1,L_NGAUSS
423	! if(gasv8(nt,np,nh,nw,ng).lt.1.0e-25)then
424	! fzerov(nw)=fzerov(nw)+1.d0
425	! endif
426	! end do
427	! end do
428	! end do
429	! end do
430	! fzerov(nw)=fzerov(nw)/dble(L_NTREFL_NPREFL_REFVAR*L_NGAUSS)
431	end do
432
433	endif
434
435	!$OMP MASTER
436	if(nIR_limit.eq.0) then
437	gasi8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTI,1:L_NGAUSS)= &
438	gasi8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTI,1:L_NGAUSS)+kappa_VI
439	else if (nIR_limit.eq.L_NSPECTI) then
440	gasi8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTI,1:L_NGAUSS)= &
441	gasi8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTI,1:L_NGAUSS)+kappa_IR
442	else
443	gasi8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:nIR_limit,1:L_NGAUSS)= &
444	gasi8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:nIR_limit,1:L_NGAUSS)+kappa_IR
445	gasi8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,nIR_limit+1:L_NSPECTI,1:L_NGAUSS)= &
446	gasi8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,nIR_limit+1:L_NSPECTI,1:L_NGAUSS)+kappa_VI
447	end if
448
449
450	! Take log10 of the values - this is what we will interpolate.
451	! Smallest value is 1.0E-200.
452
453	do nt=1,L_NTREF
454	do np=1,L_NPREF
455	do nh=1,L_REFVAR
456	do ng = 1,L_NGAUSS
457
458	do nw=1,L_NSPECTV
459	if(gasv8(nt,np,nh,nw,ng).gt.1.0d-200) then
460	gasv8(nt,np,nh,nw,ng) = log10(gasv8(nt,np,nh,nw,ng))
461	else
462	gasv8(nt,np,nh,nw,ng) = -200.0
463	end if
464	end do
465
466	do nw=1,L_NSPECTI
467	if(gasi8(nt,np,nh,nw,ng).gt.1.0d-200) then
468	gasi8(nt,np,nh,nw,ng) = log10(gasi8(nt,np,nh,nw,ng))
469	else
470	gasi8(nt,np,nh,nw,ng) = -200.0
471	end if
472	end do
473
474	end do
475	end do
476	end do
477	end do
478	!$OMP END MASTER
479	!$OMP BARRIER
480
481	! Interpolate the values: first the longwave
482
483	do nt=1,L_NTREF
484	do nh=1,L_REFVAR
485	do nw=1,L_NSPECTI
486	do ng=1,L_NGAUSS
487
488	! First, the initial interval
489
490	n = 1
491	do m=1,5
492	x = pfgasref(m)
493	xi(1) = pgasref(n)
494	xi(2) = pgasref(n+1)
495	xi(3) = pgasref(n+2)
496	xi(4) = pgasref(n+3)
497	yi(1) = gasi8(nt,n,nh,nw,ng)
498	yi(2) = gasi8(nt,n+1,nh,nw,ng)
499	yi(3) = gasi8(nt,n+2,nh,nw,ng)
500	yi(4) = gasi8(nt,n+3,nh,nw,ng)
501	call lagrange(x,xi,yi,ans)
502	gasi(nt,m,nh,nw,ng) = 10.0**ans
503	end do
504
505	do n=2,L_NPREF-2
506	do m=1,5
507	i = (n-1)*5+m
508	x = pfgasref(i)
509	xi(1) = pgasref(n-1)
510	xi(2) = pgasref(n)
511	xi(3) = pgasref(n+1)
512	xi(4) = pgasref(n+2)
513	yi(1) = gasi8(nt,n-1,nh,nw,ng)
514	yi(2) = gasi8(nt,n,nh,nw,ng)
515	yi(3) = gasi8(nt,n+1,nh,nw,ng)
516	yi(4) = gasi8(nt,n+2,nh,nw,ng)
517	call lagrange(x,xi,yi,ans)
518	gasi(nt,i,nh,nw,ng) = 10.0**ans
519	end do
520	end do
521
522	! Now, get the last interval
523
524	n = L_NPREF-1
525	do m=1,5
526	i = (n-1)*5+m
527	x = pfgasref(i)
528	xi(1) = pgasref(n-2)
529	xi(2) = pgasref(n-1)
530	xi(3) = pgasref(n)
531	xi(4) = pgasref(n+1)
532	yi(1) = gasi8(nt,n-2,nh,nw,ng)
533	yi(2) = gasi8(nt,n-1,nh,nw,ng)
534	yi(3) = gasi8(nt,n,nh,nw,ng)
535	yi(4) = gasi8(nt,n+1,nh,nw,ng)
536	call lagrange(x,xi,yi,ans)
537	gasi(nt,i,nh,nw,ng) = 10.0**ans
538	end do
539
540	! Fill the last pressure point
541
542	gasi(nt,L_PINT,nh,nw,ng) = &
543	10.0**gasi8(nt,L_NPREF,nh,nw,ng)
544
545	end do
546	end do
547	end do
548	end do
549
550	! Interpolate the values: now the shortwave
551
552	do nt=1,L_NTREF
553	do nh=1,L_REFVAR
554	do nw=1,L_NSPECTV
555	do ng=1,L_NGAUSS
556
557	! First, the initial interval
558
559	n = 1
560	do m=1,5
561	x = pfgasref(m)
562	xi(1) = pgasref(n)
563	xi(2) = pgasref(n+1)
564	xi(3) = pgasref(n+2)
565	xi(4) = pgasref(n+3)
566	yi(1) = gasv8(nt,n,nh,nw,ng)
567	yi(2) = gasv8(nt,n+1,nh,nw,ng)
568	yi(3) = gasv8(nt,n+2,nh,nw,ng)
569	yi(4) = gasv8(nt,n+3,nh,nw,ng)
570	call lagrange(x,xi,yi,ans)
571	gasv(nt,m,nh,nw,ng) = 10.0**ans
572	end do
573
574	do n=2,L_NPREF-2
575	do m=1,5
576	i = (n-1)*5+m
577	x = pfgasref(i)
578	xi(1) = pgasref(n-1)
579	xi(2) = pgasref(n)
580	xi(3) = pgasref(n+1)
581	xi(4) = pgasref(n+2)
582	yi(1) = gasv8(nt,n-1,nh,nw,ng)
583	yi(2) = gasv8(nt,n,nh,nw,ng)
584	yi(3) = gasv8(nt,n+1,nh,nw,ng)
585	yi(4) = gasv8(nt,n+2,nh,nw,ng)
586	call lagrange(x,xi,yi,ans)
587	gasv(nt,i,nh,nw,ng) = 10.0**ans
588	end do
589	end do
590
591	! Now, get the last interval
592
593	n = L_NPREF-1
594	do m=1,5
595	i = (n-1)*5+m
596	x = pfgasref(i)
597	xi(1) = pgasref(n-2)
598	xi(2) = pgasref(n-1)
599	xi(3) = pgasref(n)
600	xi(4) = pgasref(n+1)
601	yi(1) = gasv8(nt,n-2,nh,nw,ng)
602	yi(2) = gasv8(nt,n-1,nh,nw,ng)
603	yi(3) = gasv8(nt,n,nh,nw,ng)
604	yi(4) = gasv8(nt,n+1,nh,nw,ng)
605	call lagrange(x,xi,yi,ans)
606	gasv(nt,i,nh,nw,ng) = 10.0**ans
607	end do
608
609	! Fill the last pressure point
610
611	gasv(nt,L_PINT,nh,nw,ng) = &
612	10.0**gasv8(nt,L_NPREF,nh,nw,ng)
613
614	end do
615	end do
616	end do
617	end do
618
619
620	!=======================================================================
621	! Initialise the continuum absorption data
622	if(continuum)then
623	do igas=1,ngasmx
624
625	if (igas .eq. igas_N2) then
626
627	dummy = -9999
628	call interpolateN2N2(100.D+0,250.D+0,17500.D+0,testcont,.true.,dummy)
629
630	elseif (igas .eq. igas_H2) then
631
632	! first do self-induced absorption
633	dummy = -9999
634	call interpolateH2H2(500.D+0,250.D+0,17500.D+0,testcont,.true.,dummy)
635	! then cross-interactions with other gases
636	do jgas=1,ngasmx
637	if (jgas .eq. igas_N2) then
638	dummy = -9999
639	call interpolateN2H2(592.D+0,278.15D+0,200000.D+0,10000.D+0,testcont,.true.,dummy)
640	elseif (jgas .eq. igas_He) then
641	dummy = -9999
642	call interpolateH2He(500.D+0,250.D+0,200000.D+0,10000.D+0,testcont,.true.,dummy)
643	endif
644	enddo
645
646	elseif (igas .eq. igas_H2O) then
647
648	! H2O is special
649	if(H2Ocont_simple)then
650	call interpolateH2Ocont_PPC(990.D+0,296.D+0,683.2D+0*2,0.D+0,testcont,.true.)
651	else
652	dummy = -9999
653	call interpolateH2Ocont_CKD(990.D+0,296.D+0,683.2D+0*2,0.D+0,testcont,.true.,dummy)
654	endif
655
656	endif
657
658	enddo
659	endif
660
661	print*,'----------------------------------------------------'
662	print*,'And that`s all we have. It`s possible that other'
663	print*,'continuum absorption may be present, but if it is we'
664	print*,'don`t yet have data for it...'
665	print*,''
666
667	! Deallocate local arrays
668	!$OMP BARRIER
669	!$OMP MASTER
670	IF( ALLOCATED( gasi8 ) ) DEALLOCATE( gasi8 )
671	IF( ALLOCATED( gasv8 ) ) DEALLOCATE( gasv8 )
672	IF( ALLOCATED( pgasref ) ) DEALLOCATE( pgasref )
673	IF( ALLOCATED( gastype ) ) DEALLOCATE( gastype )
674	!$OMP END MASTER
675	!$OMP BARRIER
676
677	return
678	end subroutine sugas_corrk

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

Download in other formats:

Original Format