Context Navigation

sugas_corrk.F90 @ 298

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

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

Context Navigation

source: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/sugas_corrk.F90 @ 298

Download in other formats: