Context Navigation

sugas_corrk.F90 @ 728

Last change on this file since 728 was 716, checked in by rwordsworth, 13 years ago

Mainly updates to radiative transfer and gas management scheme.
Most CIA data now read from standard HITRAN datafiles. For the H2O
continuum, two options have been added: the standard CKD continuum,
and the empirical formula in PPC (Pierrehumbert 2010). Use the toggle
'H2Ocont_simple' in callphys.def to choose.

Note to Martians: I've changed the default values of 'sedimentation' and
'co2cond' in inifis to false. Both these are defined in the standard deftank
callphys.def file, so there shouldn't be any compatibility problems.

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

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

Download in other formats:

Original Format