Context Navigation

sugas_corrk.F90 @ 370

Last change on this file since 370 was 368, checked in by rwordsworth, 14 years ago
.nc saves in physiq.F90 adapted to new diagfi.def format. More info added to gases.def warnings.
File size: 17.7 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,gastype
25
26	implicit none
27
28	#include "datafile.h"
29	#include "callkeys.h"
30	#include "gases.h"
31
32	!==================================================================
33
34	logical file_ok
35
36	integer n, nt, np, nh, ng, nw, m, i
37	integer L_NGAUSScheck, L_NPREFcheck, L_NTREFcheck, L_REFVARcheck
38
39	character(len=100) :: file_id
40	character(len=100) :: file_path
41
42	real*8 gasi8(L_NTREF,L_NPREF,L_REFVAR,L_NSPECTI,L_NGAUSS)
43	real*8 gasv8(L_NTREF,L_NPREF,L_REFVAR,L_NSPECTV,L_NGAUSS)
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/' // corrkdir(1:LEN_TRIM(corrkdir)) // '/Q.dat'
54	file_path=datafile(1:LEN_TRIM(datafile))//file_id(1:LEN_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 ',file_path(1:LEN_TRIM(file_path))
61	write(,)'was not found by sugas_corrk.F90, exiting.'
62	write(,)'Check that your path to datagcm in datafile.h is correct and that'
63	write(,)'the corrkdir you chose in callphys.def exists.'
64	call abort
65	endif
66
67	! check that database matches varactive toggle
68	open(111,file=file_path(1:LEN_TRIM(file_path)),form='formatted')
69	read(111,*) ngas
70
71	if(ngas.ne.ngasmx)then
72	print*,'Number of gases in radiative transfer data (',ngas,') does not', &
73	'match that in gases.def (',ngasmx,'), exiting.'
74	call abort
75	endif
76
77	if(ngas.eq.1 .and. (varactive.or.varfixed))then
78	print*,'You have varactive/fixed=.true. but the database [', &
79	corrkdir(1:LEN_TRIM(corrkdir)), &
80	'] has no variable species, exiting.'
81	call abort
82	elseif(ngas.eq.2 .and. (.not.varactive) .and. (.not.varfixed))then
83	print*,'You have varactive and varfixed=.false. and the database [', &
84	corrkdir(1:LEN_TRIM(corrkdir)), &
85	'] has a variable species.'
86	call abort
87	elseif(ngas.gt.4 .or. ngas.lt.1)then
88	print*,ngas,' species in database [', &
89	corrkdir(1:LEN_TRIM(corrkdir)), &
90	'], radiative code cannot handle this.'
91	call abort
92	endif
93
94	if(ngas.gt.3)then
95	print*,'ngas>3, EXPERIMENTAL!'
96	endif
97
98
99	do n=1,ngas
100	read(111,*) gastype(n)
101	print*,'Gas ',n,' is ',gastype(n)
102	enddo
103
104	! check the array size is correct, load the coefficients
105	open(111,file=file_path(1:LEN_TRIM(file_path)),form='formatted')
106	read(111,*) L_REFVARcheck
107	if(.not.(L_REFVARcheck.eq.L_REFVAR)) then
108	print*, L_REFVARcheck
109	print*, L_REFVAR
110	print*,'The size of your radiative transfer mixing ratio array does '
111	print*,'not match the value given in Q.dat, exiting.'
112	print*,'Check the value of L_NREFVAR in radinc_h.F90.'
113	call abort
114	endif
115	read(111,*) wrefvar
116	close(111)
117
118	! Check that gastype and gnom match
119	do n=1,ngas
120	print*,'Gas ',n,' is ',gnom(n)
121	if(gnom(n).ne.gastype(n))then
122	print*,'Name of a gas in radiative transfer data (',gastype(n),') does not ', &
123	'match that in gases.def (',gnom(n),'), 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 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/' // corrkdir(1:LEN_TRIM(corrkdir)) // '/g.dat'
142	file_path=datafile(1:LEN_TRIM(datafile))//file_id(1:LEN_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 ',file_path(1:LEN_TRIM(file_path))
148	write(,)'was not found by sugas_corrk.F90, exiting.'
149	write(,)'Check that your path to datagcm in datafile.h is correct and that'
150	write(,)'the corrkdir you chose in callphys.def exists.'
151	call abort
152	endif
153
154	! check the array size is correct, load the coefficients
155	open(111,file=file_path(1:LEN_TRIM(file_path)),form='formatted')
156	read(111,*) L_NGAUSScheck
157	if(.not.(L_NGAUSScheck.eq.L_NGAUSS)) then
158	print*,'The size of your radiative transfer g-space array does '
159	print*,'not match the value given in g.dat, exiting.'
160	call abort
161	endif
162	read(111,*) gweight
163	close(111)
164
165	! display the values
166	print*,'Correlated-k g-space grid:'
167	do n=1,L_NGAUSS
168	print*,n,'.',gweight(n)
169	end do
170	print*,''
171
172	!=======================================================================
173	! Set the reference pressure and temperature arrays. These are
174	! the pressures and temperatures at which we have k-coefficients.
175
176	!-----------------------------------------------------------------------
177	! pressure
178
179	file_id='/corrk_data/' // corrkdir(1:LEN_TRIM(corrkdir)) // '/p.dat'
180	file_path=datafile(1:LEN_TRIM(datafile))//file_id(1:LEN_TRIM(file_id))
181
182	! check that the file exists
183	inquire(FILE=file_path,EXIST=file_ok)
184	if(.not.file_ok) then
185	write(,)'The file ',file_path(1:LEN_TRIM(file_path))
186	write(,)'was not found by sugas_corrk.F90, exiting.'
187	call abort
188	endif
189
190	! check the array size is correct, load the coefficients
191	open(111,file=file_path(1:LEN_TRIM(file_path)),form='formatted')
192	read(111,*) L_NPREFcheck
193	if(.not.(L_NPREFcheck.eq.L_NPREF)) then
194	print*,'The size of your radiative transfer pressure array does '
195	print*,'not match the value given in p.dat, exiting.'
196	call abort
197	endif
198	read(111,*) pgasref
199	close(111)
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/' // corrkdir(1:LEN_TRIM(corrkdir)) // '/T.dat'
225	file_path=datafile(1:LEN_TRIM(datafile))//file_id(1:LEN_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 ',file_path(1:LEN_TRIM(file_path))
231	write(,)'was not found by sugas_corrk.F90, exiting.'
232	call abort
233	endif
234
235	! check the array size is correct, load the coefficients
236	open(111,file=file_path(1:LEN_TRIM(file_path)),form='formatted')
237	read(111,*) L_NTREFcheck
238	if(.not.(L_NTREFcheck.eq.L_NTREF)) then
239	print*,'The size of your radiative transfer temperature array does '
240	print*,'not match the value given in T.dat, exiting.'
241	call abort
242	endif
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	! Get gaseous k-coefficients and interpolate onto finer pressure grid
258
259	! VISIBLE
260	if (callgasvis.and..not.corrkdir(1:LEN_TRIM(corrkdir)).eq.'null') then
261	file_id='/corrk_data/'//trim(adjustl(banddir))//'/corrk_gcm_VI.dat'
262	file_path=datafile(1:LEN_TRIM(datafile))//file_id(1:LEN_TRIM(file_id))
263
264	! check that the file exists
265	inquire(FILE=file_path,EXIST=file_ok)
266	if(.not.file_ok) then
267	write(,)'The file ',file_path(1:LEN_TRIM(file_path))
268	write(,)'was not found by sugas_corrk.F90.'
269	write(,)'Are you sure you have absorption data for these bands?'
270	call abort
271	endif
272
273	open(111,file=file_path(1:LEN_TRIM(file_path)),form='formatted')
274	read(111,*) gasv8
275	close(111)
276
277	else
278	print*,'Visible corrk gaseous absorption is set to zero.'
279	gasv8(:,:,:,:,:)=0.0
280	endif
281
282	! INFRA-RED
283	if (.not.corrkdir(1:LEN_TRIM(corrkdir)).eq.'null') then
284	file_id='/corrk_data/'//trim(adjustl(banddir))//'/corrk_gcm_IR.dat'
285	file_path=datafile(1:LEN_TRIM(datafile))//file_id(1:LEN_TRIM(file_id))
286
287	! check that the file exists
288	inquire(FILE=file_path,EXIST=file_ok)
289	if(.not.file_ok) then
290	write(,)'The file ',file_path(1:LEN_TRIM(file_path))
291	write(,)'was not found by sugas_corrk.F90.'
292	write(,)'Are you sure you have absorption data for these bands?'
293	call abort
294	endif
295
296	open(111,file=file_path(1:LEN_TRIM(file_path)),form='formatted')
297	read(111,*) gasi8
298	close(111)
299
300	do nw=1,L_NSPECTI
301	fzeroi(nw) = 0.
302	! do nt=1,L_NTREF
303	! do np=1,L_NPREF
304	! do nh=1,L_REFVAR
305	! do ng = 1,L_NGAUSS
306	! if(gasi8(nt,np,nh,nw,ng).lt.1.0e-25)then
307	! fzeroi(nw)=fzeroi(nw)+1.
308	! endif
309	! end do
310	! end do
311	! end do
312	! end do
313	! fzeroi(nw)=fzeroi(nw)/dble(L_NTREFL_NPREFL_REFVAR*L_NGAUSS)
314	end do
315
316	do nw=1,L_NSPECTV
317	fzerov(nw) = 0.
318	! do nt=1,L_NTREF
319	! do np=1,L_NPREF
320	! do nh=1,L_REFVAR
321	! do ng = 1,L_NGAUSS
322	! if(gasv8(nt,np,nh,nw,ng).lt.1.0e-25)then
323	! fzerov(nw)=fzerov(nw)+1.
324	! endif
325	! end do
326	! end do
327	! end do
328	! end do
329	! fzerov(nw)=fzerov(nw)/dble(L_NTREFL_NPREFL_REFVAR*L_NGAUSS)
330	end do
331
332
333	do nw=1,L_NSPECTV
334	fzerov(nw) = 0.
335	end do
336
337	else
338	print*,'Infrared corrk gaseous absorption is set to zero.'
339	gasi8(:,:,:,:,:)=0.0
340	endif
341
342	! Take log10 of the values - this is what we will interpolate.
343	! Smallest value is 1.0E-200.
344
345	do nt=1,L_NTREF
346	do np=1,L_NPREF
347	do nh=1,L_REFVAR
348	do ng = 1,L_NGAUSS
349
350	do nw=1,L_NSPECTV
351	if(gasv8(nt,np,nh,nw,ng).gt.1.0d-200) then
352	gasv8(nt,np,nh,nw,ng) = &
353	log10(gasv8(nt,np,nh,nw,ng))
354	else
355	gasv8(nt,np,nh,nw,ng) = -200.0
356	end if
357	end do
358
359	do nw=1,L_NSPECTI
360	if(gasi8(nt,np,nh,nw,ng).gt.1.0d-200) then
361	gasi8(nt,np,nh,nw,ng) = &
362	log10(gasi8(nt,np,nh,nw,ng))
363	else
364	gasi8(nt,np,nh,nw,ng) = -200.0
365	end if
366	end do
367
368	end do
369	end do
370	end do
371	end do
372
373	! Interpolate the values: first the longwave
374
375	do nt=1,L_NTREF
376	do nh=1,L_REFVAR
377	do nw=1,L_NSPECTI
378	do ng=1,L_NGAUSS
379
380	! First, the initial interval
381
382	n = 1
383	do m=1,5
384	x = pfgasref(m)
385	xi(1) = pgasref(n)
386	xi(2) = pgasref(n+1)
387	xi(3) = pgasref(n+2)
388	xi(4) = pgasref(n+3)
389	yi(1) = gasi8(nt,n,nh,nw,ng)
390	yi(2) = gasi8(nt,n+1,nh,nw,ng)
391	yi(3) = gasi8(nt,n+2,nh,nw,ng)
392	yi(4) = gasi8(nt,n+3,nh,nw,ng)
393	call lagrange(x,xi,yi,ans)
394	gasi(nt,m,nh,nw,ng) = 10.0**ans
395	end do
396
397	do n=2,L_NPREF-2
398	do m=1,5
399	i = (n-1)*5+m
400	x = pfgasref(i)
401	xi(1) = pgasref(n-1)
402	xi(2) = pgasref(n)
403	xi(3) = pgasref(n+1)
404	xi(4) = pgasref(n+2)
405	yi(1) = gasi8(nt,n-1,nh,nw,ng)
406	yi(2) = gasi8(nt,n,nh,nw,ng)
407	yi(3) = gasi8(nt,n+1,nh,nw,ng)
408	yi(4) = gasi8(nt,n+2,nh,nw,ng)
409	call lagrange(x,xi,yi,ans)
410	gasi(nt,i,nh,nw,ng) = 10.0**ans
411	end do
412	end do
413
414	! Now, get the last interval
415
416	n = L_NPREF-1
417	do m=1,5
418	i = (n-1)*5+m
419	x = pfgasref(i)
420	xi(1) = pgasref(n-2)
421	xi(2) = pgasref(n-1)
422	xi(3) = pgasref(n)
423	xi(4) = pgasref(n+1)
424	yi(1) = gasi8(nt,n-2,nh,nw,ng)
425	yi(2) = gasi8(nt,n-1,nh,nw,ng)
426	yi(3) = gasi8(nt,n,nh,nw,ng)
427	yi(4) = gasi8(nt,n+1,nh,nw,ng)
428	call lagrange(x,xi,yi,ans)
429	gasi(nt,i,nh,nw,ng) = 10.0**ans
430	end do
431
432	! Fill the last pressure point
433
434	gasi(nt,L_PINT,nh,nw,ng) = &
435	10.0**gasi8(nt,L_NPREF,nh,nw,ng)
436
437	end do
438	end do
439	end do
440	end do
441
442	! Interpolate the values: now the shortwave
443
444	do nt=1,L_NTREF
445	do nh=1,L_REFVAR
446	do nw=1,L_NSPECTV
447	do ng=1,L_NGAUSS
448
449	! First, the initial interval
450
451	n = 1
452	do m=1,5
453	x = pfgasref(m)
454	xi(1) = pgasref(n)
455	xi(2) = pgasref(n+1)
456	xi(3) = pgasref(n+2)
457	xi(4) = pgasref(n+3)
458	yi(1) = gasv8(nt,n,nh,nw,ng)
459	yi(2) = gasv8(nt,n+1,nh,nw,ng)
460	yi(3) = gasv8(nt,n+2,nh,nw,ng)
461	yi(4) = gasv8(nt,n+3,nh,nw,ng)
462	call lagrange(x,xi,yi,ans)
463	gasv(nt,m,nh,nw,ng) = 10.0**ans
464	end do
465
466	do n=2,L_NPREF-2
467	do m=1,5
468	i = (n-1)*5+m
469	x = pfgasref(i)
470	xi(1) = pgasref(n-1)
471	xi(2) = pgasref(n)
472	xi(3) = pgasref(n+1)
473	xi(4) = pgasref(n+2)
474	yi(1) = gasv8(nt,n-1,nh,nw,ng)
475	yi(2) = gasv8(nt,n,nh,nw,ng)
476	yi(3) = gasv8(nt,n+1,nh,nw,ng)
477	yi(4) = gasv8(nt,n+2,nh,nw,ng)
478	call lagrange(x,xi,yi,ans)
479	gasv(nt,i,nh,nw,ng) = 10.0**ans
480	end do
481	end do
482
483	! Now, get the last interval
484
485	n = L_NPREF-1
486	do m=1,5
487	i = (n-1)*5+m
488	x = pfgasref(i)
489	xi(1) = pgasref(n-2)
490	xi(2) = pgasref(n-1)
491	xi(3) = pgasref(n)
492	xi(4) = pgasref(n+1)
493	yi(1) = gasv8(nt,n-2,nh,nw,ng)
494	yi(2) = gasv8(nt,n-1,nh,nw,ng)
495	yi(3) = gasv8(nt,n,nh,nw,ng)
496	yi(4) = gasv8(nt,n+1,nh,nw,ng)
497	call lagrange(x,xi,yi,ans)
498	gasv(nt,i,nh,nw,ng) = 10.0**ans
499	end do
500
501	! Fill the last pressure point
502
503	gasv(nt,L_PINT,nh,nw,ng) = &
504	10.0**gasv8(nt,L_NPREF,nh,nw,ng)
505
506	end do
507	end do
508	end do
509	end do
510
511
512
513	do igas=1,ngasmx
514	if(gnom(igas).eq.'H2_')then
515	call interpolateH2H2(500.D+0,250.D+0,17500.D+0,testcont,.true.)
516	elseif(gnom(igas).eq.'H2O')then
517	call interpolateH2Ocont(990.D+0,296.D+0,683.2D+0*2,0.D+0,testcont,.true.)
518	endif
519	enddo
520
521	return
522	end subroutine sugas_corrk

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

Download in other formats:

Original Format