Context Navigation

sugas_corrk.F90 @ 537

Last change on this file since 537 was 484, checked in by aslmd, 13 years ago
GENERIC: Allocatable gastype in sugas_corrk instead of hardcoded (it was a problem for more than 4 gases) JL + AS.
File size: 18.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	!
15	! Summary
16	! -------
17	!
18	!==================================================================
19
20	use radinc_h
21	use radcommon_h, only : pgasref,pfgasref,pgasmin,pgasmax
22	use radcommon_h, only : tgasref,tgasmin,tgasmax
23	use radcommon_h, only : gasv,gasi,FZEROI,FZEROV,gweight
24	use radcommon_h, only : wrefvar
25	use datafile_mod, only: datadir
26	use gases_h
27	implicit none
28
29	#include "callkeys.h"
30
31	!==================================================================
32
33	logical file_ok
34
35	integer n, nt, np, nh, ng, nw, m, i
36	integer L_NGAUSScheck
37
38	character(len=100) :: file_id
39	character(len=300) :: file_path
40
41	!! ALLOCATABLE ARRAYS -- AS 12/2011
42	REAL*8, DIMENSION(:,:,:,:,:), ALLOCATABLE :: gasi8, gasv8
43	character*3,allocatable,DIMENSION(:) :: gastype ! for check with gnom
44
45	real*8 x, xi(4), yi(4), ans, p
46
47	integer ngas, igas
48
49	double precision testcont ! for continuum absorption initialisation
50
51	!=======================================================================
52	! Load variable species data, exit if we have wrong database
53	file_id='/corrk_data/' // TRIM(corrkdir) // '/Q.dat'
54	file_path=TRIM(datadir)//TRIM(file_id)
55
56
57	! check that the file exists
58	inquire(FILE=file_path,EXIST=file_ok)
59	if(.not.file_ok) then
60	write(,)'The file ',TRIM(file_path)
61	write(,)'was not found by sugas_corrk.F90, exiting.'
62	write(,)'Check that your path to datagcm:',trim(datadir)
63	write(,)' is correct. You can change it in callphys.def with:'
64	write(,)' datadir = /absolute/path/to/datagcm'
65	write(,)'Also check that the corrkdir you chose in callphys.def exists.'
66	call abort
67	endif
68
69	! check that database matches varactive toggle
70	open(111,file=TRIM(file_path),form='formatted')
71	read(111,*) ngas
72
73	if(ngas.ne.ngasmx)then
74	print*,'Number of gases in radiative transfer data (',ngas,') does not', &
75	'match that in gases.def (',ngasmx,'), exiting.'
76	call abort
77	endif
78
79	if(ngas.eq.1 .and. (varactive.or.varfixed))then
80	print*,'You have varactive/fixed=.true. but the database [', &
81	corrkdir(1:LEN_TRIM(corrkdir)), &
82	'] has no variable species, exiting.'
83	call abort
84	elseif(ngas.eq.2 .and. (.not.varactive) .and. (.not.varfixed))then
85	print*,'You have varactive and varfixed=.false. and the database [', &
86	corrkdir(1:LEN_TRIM(corrkdir)), &
87	'] has a variable species.'
88	call abort
89	elseif(ngas.gt.4 .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	if(ngas.gt.3)then
97	print*,'ngas>3, EXPERIMENTAL!'
98	endif
99
100	IF ( .NOT. ALLOCATED( gastype ) ) ALLOCATE( gastype( ngas ) )
101	do n=1,ngas
102	read(111,*) gastype(n)
103	print*,'Gas ',n,' is ',gastype(n)
104	enddo
105
106	! GET array size, load the coefficients
107	open(111,file=TRIM(file_path),form='formatted')
108	read(111,*) L_REFVAR
109	IF( .NOT. ALLOCATED( wrefvar ) ) ALLOCATE( WREFVAR(L_REFVAR) )
110	read(111,*) wrefvar
111	close(111)
112
113	! Check that gastype and gnom match
114	do n=1,ngas
115	print*,'Gas ',n,' is ',gnom(n)
116	if(gnom(n).ne.gastype(n))then
117	print*,'Name of a gas in radiative transfer data (',gastype(n),') does not ', &
118	'match that in gases.def (',gnom(n),'), exiting. You should compare ', &
119	'gases.def with Q.dat in your radiative transfer directory.'
120	call abort
121	endif
122	enddo
123	print*,'Confirmed gas match in radiative transfer and gases.def!'
124
125	! display the values
126	print*,'Variable gas mixing ratios:'
127	do n=1,L_REFVAR
128	!print*,n,'.',wrefvar(n),' kg/kg' ! pay attention!
129	print*,n,'.',wrefvar(n),' mol/mol'
130	end do
131	print*,''
132
133	!=======================================================================
134	! Set the weighting in g-space
135
136	file_id='/corrk_data/' // TRIM(corrkdir) // '/g.dat'
137	file_path=TRIM(datadir)//TRIM(file_id)
138
139	! check that the file exists
140	inquire(FILE=file_path,EXIST=file_ok)
141	if(.not.file_ok) then
142	write(,)'The file ',TRIM(file_path)
143	write(,)'was not found by sugas_corrk.F90, exiting.'
144	write(,)'Check that your path to datagcm:',trim(datadir)
145	write(,)' is correct. You can change it in callphys.def with:'
146	write(,)' datadir = /absolute/path/to/datagcm'
147	write(,)'Also check that the corrkdir you chose in callphys.def exists.'
148	call abort
149	endif
150
151	! check the array size is correct, load the coefficients
152	open(111,file=TRIM(file_path),form='formatted')
153	read(111,*) L_NGAUSScheck
154	if(.not.(L_NGAUSScheck.eq.L_NGAUSS)) then
155	print*,'The size of your radiative transfer g-space array does '
156	print*,'not match the value given in g.dat, exiting.'
157	call abort
158	endif
159	read(111,*) gweight
160	close(111)
161
162	! display the values
163	print*,'Correlated-k g-space grid:'
164	do n=1,L_NGAUSS
165	print*,n,'.',gweight(n)
166	end do
167	print*,''
168
169	!=======================================================================
170	! Set the reference pressure and temperature arrays. These are
171	! the pressures and temperatures at which we have k-coefficients.
172
173	!-----------------------------------------------------------------------
174	! pressure
175
176	file_id='/corrk_data/' // TRIM(corrkdir) // '/p.dat'
177	file_path=TRIM(datadir)//TRIM(file_id)
178
179	! check that the file exists
180	inquire(FILE=file_path,EXIST=file_ok)
181	if(.not.file_ok) then
182	write(,)'The file ',TRIM(file_path)
183	write(,)'was not found by sugas_corrk.F90, exiting.'
184	write(,)'Check that your path to datagcm:',trim(datadir)
185	write(,)' is correct. You can change it in callphys.def with:'
186	write(,)' datadir = /absolute/path/to/datagcm'
187	write(,)'Also check that the corrkdir you chose in callphys.def exists.'
188	call abort
189	endif
190
191	! GET array size, load the coefficients
192	open(111,file=TRIM(file_path),form='formatted')
193	read(111,*) L_NPREF
194	IF( .NOT. ALLOCATED( pgasref ) ) ALLOCATE( PGASREF(L_NPREF) )
195	read(111,*) pgasref
196	close(111)
197	L_PINT = (L_NPREF-1)*5+1
198	IF( .NOT. ALLOCATED( pfgasref ) ) ALLOCATE( PFGASREF(L_PINT) )
199
200
201	! display the values
202	print*,'Correlated-k pressure grid (mBar):'
203	do n=1,L_NPREF
204	print*,n,'. 1 x 10^',pgasref(n),' mBar'
205	end do
206	print*,''
207
208	! save the min / max matrix values
209	pgasmin = 10.0**pgasref(1)
210	pgasmax = 10.0**pgasref(L_NPREF)
211
212	! interpolate to finer grid
213	do n=1,L_NPREF-1
214	do m=1,5
215	pfgasref((n-1)5+m) = pgasref(n)+(m-1)0.2
216	end do
217	end do
218	pfgasref(L_PINT) = pgasref(L_NPREF)
219	print*,'Warning, pfgasref needs generalisation to uneven grids!!'
220
221	!-----------------------------------------------------------------------
222	! temperature
223
224	file_id='/corrk_data/' // TRIM(corrkdir) // '/T.dat'
225	file_path=TRIM(datadir)//TRIM(file_id)
226
227	! check that the file exists
228	inquire(FILE=file_path,EXIST=file_ok)
229	if(.not.file_ok) then
230	write(,)'The file ',TRIM(file_path)
231	write(,)'was not found by sugas_corrk.F90, exiting.'
232	write(,)'Check that your path to datagcm:',trim(datadir)
233	write(,)' is correct. You can change it in callphys.def with:'
234	write(,)' datadir = /absolute/path/to/datagcm'
235	write(,)'Also check that the corrkdir you chose in callphys.def exists.'
236	call abort
237	endif
238
239	! GET array size, load the coefficients
240	open(111,file=TRIM(file_path),form='formatted')
241	read(111,*) L_NTREF
242	IF( .NOT. ALLOCATED( tgasref ) ) ALLOCATE( TGASREF(L_NTREF) )
243	read(111,*) tgasref
244	close(111)
245
246	! display the values
247	print*,'Correlated-k temperature grid:'
248	do n=1,L_NTREF
249	print*,n,'.',tgasref(n),' K'
250	end do
251
252	! save the min / max matrix values
253	tgasmin = tgasref(1)
254	tgasmax = tgasref(L_NTREF)
255
256
257	!-----------------------------------------------------------------------
258	!-----------------------------------------------------------------------
259	! ALLOCATE THE MULTIDIM ARRAYS IN radcommon_h
260	PRINT *, L_NTREF,L_NPREF,L_REFVAR
261	IF( .NOT. ALLOCATED( gasi8 ) ) ALLOCATE( gasi8(L_NTREF,L_NPREF,L_REFVAR,L_NSPECTI,L_NGAUSS) )
262	IF( .NOT. ALLOCATED( gasv8 ) ) ALLOCATE( gasv8(L_NTREF,L_NPREF,L_REFVAR,L_NSPECTV,L_NGAUSS) )
263	IF( .NOT. ALLOCATED( gasi ) ) ALLOCATE( gasi(L_NTREF,L_PINT,L_REFVAR,L_NSPECTI,L_NGAUSS) )
264	IF( .NOT. ALLOCATED( gasv ) ) ALLOCATE( gasv(L_NTREF,L_PINT,L_REFVAR,L_NSPECTV,L_NGAUSS) )
265	!-----------------------------------------------------------------------
266	!-----------------------------------------------------------------------
267
268
269
270	!=======================================================================
271	! Get gaseous k-coefficients and interpolate onto finer pressure grid
272
273	! VISIBLE
274	if (callgasvis.and..not.TRIM(corrkdir).eq.'null') then
275	file_id='/corrk_data/'//trim(adjustl(banddir))//'/corrk_gcm_VI.dat'
276	file_path=TRIM(datadir)//TRIM(file_id)
277
278	! check that the file exists
279	inquire(FILE=file_path,EXIST=file_ok)
280	if(.not.file_ok) then
281	write(,)'The file ',TRIM(file_path)
282	write(,)'was not found by sugas_corrk.F90.'
283	write(,)'Are you sure you have absorption data for these bands?'
284	call abort
285	endif
286
287	open(111,file=TRIM(file_path),form='formatted')
288	read(111,*) gasv8
289	close(111)
290
291	else
292	print*,'Visible corrk gaseous absorption is set to zero.'
293	gasv8(:,:,:,:,:)=0.0
294	endif
295
296	! INFRA-RED
297	if (.not.TRIM(corrkdir).eq.'null') then
298	file_id='/corrk_data/'//trim(adjustl(banddir))//'/corrk_gcm_IR.dat'
299	file_path=TRIM(datadir)//TRIM(file_id)
300
301	! check that the file exists
302	inquire(FILE=file_path,EXIST=file_ok)
303	if(.not.file_ok) then
304	write(,)'The file ',TRIM(file_path)
305	write(,)'was not found by sugas_corrk.F90.'
306	write(,)'Are you sure you have absorption data for these bands?'
307	call abort
308	endif
309
310	open(111,file=TRIM(file_path),form='formatted')
311	read(111,*) gasi8
312	close(111)
313
314	do nw=1,L_NSPECTI
315	fzeroi(nw) = 0.
316	! do nt=1,L_NTREF
317	! do np=1,L_NPREF
318	! do nh=1,L_REFVAR
319	! do ng = 1,L_NGAUSS
320	! if(gasi8(nt,np,nh,nw,ng).lt.1.0e-25)then
321	! fzeroi(nw)=fzeroi(nw)+1.
322	! endif
323	! end do
324	! end do
325	! end do
326	! end do
327	! fzeroi(nw)=fzeroi(nw)/dble(L_NTREFL_NPREFL_REFVAR*L_NGAUSS)
328	end do
329
330	do nw=1,L_NSPECTV
331	fzerov(nw) = 0.
332	! do nt=1,L_NTREF
333	! do np=1,L_NPREF
334	! do nh=1,L_REFVAR
335	! do ng = 1,L_NGAUSS
336	! if(gasv8(nt,np,nh,nw,ng).lt.1.0e-25)then
337	! fzerov(nw)=fzerov(nw)+1.
338	! endif
339	! end do
340	! end do
341	! end do
342	! end do
343	! fzerov(nw)=fzerov(nw)/dble(L_NTREFL_NPREFL_REFVAR*L_NGAUSS)
344	end do
345
346
347	do nw=1,L_NSPECTV
348	fzerov(nw) = 0.
349	end do
350
351	else
352	print*,'Infrared corrk gaseous absorption is set to zero.'
353	gasi8(:,:,:,:,:)=0.0
354	endif
355
356	! Take log10 of the values - this is what we will interpolate.
357	! Smallest value is 1.0E-200.
358
359	do nt=1,L_NTREF
360	do np=1,L_NPREF
361	do nh=1,L_REFVAR
362	do ng = 1,L_NGAUSS
363
364	do nw=1,L_NSPECTV
365	if(gasv8(nt,np,nh,nw,ng).gt.1.0d-200) then
366	gasv8(nt,np,nh,nw,ng) = &
367	log10(gasv8(nt,np,nh,nw,ng))
368	else
369	gasv8(nt,np,nh,nw,ng) = -200.0
370	end if
371	end do
372
373	do nw=1,L_NSPECTI
374	if(gasi8(nt,np,nh,nw,ng).gt.1.0d-200) then
375	gasi8(nt,np,nh,nw,ng) = &
376	log10(gasi8(nt,np,nh,nw,ng))
377	else
378	gasi8(nt,np,nh,nw,ng) = -200.0
379	end if
380	end do
381
382	end do
383	end do
384	end do
385	end do
386
387	! Interpolate the values: first the longwave
388
389	do nt=1,L_NTREF
390	do nh=1,L_REFVAR
391	do nw=1,L_NSPECTI
392	do ng=1,L_NGAUSS
393
394	! First, the initial interval
395
396	n = 1
397	do m=1,5
398	x = pfgasref(m)
399	xi(1) = pgasref(n)
400	xi(2) = pgasref(n+1)
401	xi(3) = pgasref(n+2)
402	xi(4) = pgasref(n+3)
403	yi(1) = gasi8(nt,n,nh,nw,ng)
404	yi(2) = gasi8(nt,n+1,nh,nw,ng)
405	yi(3) = gasi8(nt,n+2,nh,nw,ng)
406	yi(4) = gasi8(nt,n+3,nh,nw,ng)
407	call lagrange(x,xi,yi,ans)
408	gasi(nt,m,nh,nw,ng) = 10.0**ans
409	end do
410
411	do n=2,L_NPREF-2
412	do m=1,5
413	i = (n-1)*5+m
414	x = pfgasref(i)
415	xi(1) = pgasref(n-1)
416	xi(2) = pgasref(n)
417	xi(3) = pgasref(n+1)
418	xi(4) = pgasref(n+2)
419	yi(1) = gasi8(nt,n-1,nh,nw,ng)
420	yi(2) = gasi8(nt,n,nh,nw,ng)
421	yi(3) = gasi8(nt,n+1,nh,nw,ng)
422	yi(4) = gasi8(nt,n+2,nh,nw,ng)
423	call lagrange(x,xi,yi,ans)
424	gasi(nt,i,nh,nw,ng) = 10.0**ans
425	end do
426	end do
427
428	! Now, get the last interval
429
430	n = L_NPREF-1
431	do m=1,5
432	i = (n-1)*5+m
433	x = pfgasref(i)
434	xi(1) = pgasref(n-2)
435	xi(2) = pgasref(n-1)
436	xi(3) = pgasref(n)
437	xi(4) = pgasref(n+1)
438	yi(1) = gasi8(nt,n-2,nh,nw,ng)
439	yi(2) = gasi8(nt,n-1,nh,nw,ng)
440	yi(3) = gasi8(nt,n,nh,nw,ng)
441	yi(4) = gasi8(nt,n+1,nh,nw,ng)
442	call lagrange(x,xi,yi,ans)
443	gasi(nt,i,nh,nw,ng) = 10.0**ans
444	end do
445
446	! Fill the last pressure point
447
448	gasi(nt,L_PINT,nh,nw,ng) = &
449	10.0**gasi8(nt,L_NPREF,nh,nw,ng)
450
451	end do
452	end do
453	end do
454	end do
455
456	! Interpolate the values: now the shortwave
457
458	do nt=1,L_NTREF
459	do nh=1,L_REFVAR
460	do nw=1,L_NSPECTV
461	do ng=1,L_NGAUSS
462
463	! First, the initial interval
464
465	n = 1
466	do m=1,5
467	x = pfgasref(m)
468	xi(1) = pgasref(n)
469	xi(2) = pgasref(n+1)
470	xi(3) = pgasref(n+2)
471	xi(4) = pgasref(n+3)
472	yi(1) = gasv8(nt,n,nh,nw,ng)
473	yi(2) = gasv8(nt,n+1,nh,nw,ng)
474	yi(3) = gasv8(nt,n+2,nh,nw,ng)
475	yi(4) = gasv8(nt,n+3,nh,nw,ng)
476	call lagrange(x,xi,yi,ans)
477	gasv(nt,m,nh,nw,ng) = 10.0**ans
478	end do
479
480	do n=2,L_NPREF-2
481	do m=1,5
482	i = (n-1)*5+m
483	x = pfgasref(i)
484	xi(1) = pgasref(n-1)
485	xi(2) = pgasref(n)
486	xi(3) = pgasref(n+1)
487	xi(4) = pgasref(n+2)
488	yi(1) = gasv8(nt,n-1,nh,nw,ng)
489	yi(2) = gasv8(nt,n,nh,nw,ng)
490	yi(3) = gasv8(nt,n+1,nh,nw,ng)
491	yi(4) = gasv8(nt,n+2,nh,nw,ng)
492	call lagrange(x,xi,yi,ans)
493	gasv(nt,i,nh,nw,ng) = 10.0**ans
494	end do
495	end do
496
497	! Now, get the last interval
498
499	n = L_NPREF-1
500	do m=1,5
501	i = (n-1)*5+m
502	x = pfgasref(i)
503	xi(1) = pgasref(n-2)
504	xi(2) = pgasref(n-1)
505	xi(3) = pgasref(n)
506	xi(4) = pgasref(n+1)
507	yi(1) = gasv8(nt,n-2,nh,nw,ng)
508	yi(2) = gasv8(nt,n-1,nh,nw,ng)
509	yi(3) = gasv8(nt,n,nh,nw,ng)
510	yi(4) = gasv8(nt,n+1,nh,nw,ng)
511	call lagrange(x,xi,yi,ans)
512	gasv(nt,i,nh,nw,ng) = 10.0**ans
513	end do
514
515	! Fill the last pressure point
516
517	gasv(nt,L_PINT,nh,nw,ng) = &
518	10.0**gasv8(nt,L_NPREF,nh,nw,ng)
519
520	end do
521	end do
522	end do
523	end do
524
525
526
527	do igas=1,ngasmx
528	if(gnom(igas).eq.'H2_')then
529	call interpolateH2H2(500.D+0,250.D+0,17500.D+0,testcont,.true.)
530	elseif(gnom(igas).eq.'H2O')then
531	call interpolateH2Ocont(990.D+0,296.D+0,683.2D+0*2,0.D+0,testcont,.true.)
532	endif
533	enddo
534
535	!!! DEALLOCATE LOCAL ARRAYS
536	IF( ALLOCATED( gasi8 ) ) DEALLOCATE( gasi8 )
537	IF( ALLOCATED( gasv8 ) ) DEALLOCATE( gasv8 )
538	IF( ALLOCATED( pgasref ) ) DEALLOCATE( pgasref )
539	IF( ALLOCATED( gastype ) ) DEALLOCATE( gastype )
540
541	return
542	end subroutine sugas_corrk

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

Download in other formats:

Original Format