Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

sugas_corrk.F90 @ 3232

Last change on this file since 3232 was 3175, checked in by emillour, 2 years ago
Pluto PCM: Add the old Pluto LMDZ for reference (required prior step to making an LMDZ.PLUTO using the same framework as the other physics packages). TB+EM
File size: 14.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	!
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	use datafile_mod, only: datadir
26
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, L_NPREFcheck, L_NTREFcheck, L_REFVARcheck
37
38	character(len=100) :: file_id
39	character(len=100) :: file_path
40
41	real*8 gasi8(L_NTREF,L_NPREF,L_REFVAR,L_NSPECTI,L_NGAUSS)
42	real*8 gasv8(L_NTREF,L_NPREF,L_REFVAR,L_NSPECTV,L_NGAUSS)
43
44	real*8 x, xi(4), yi(4), ans, p
45
46	integer Nspecies
47
48	real*8 Xtemp
49
50	!=======================================================================
51	! Load variable species data, exit if we have wrong database
52	file_id='/corrk_data/' // corrkdir(1:LEN_TRIM(corrkdir)) // '/Q.dat'
53	file_path=TRIM(datadir)//TRIM(file_id)
54	! check that the file exists
55	inquire(FILE=file_path,EXIST=file_ok)
56	if(.not.file_ok) then
57	write(,)'The file ',file_path(1:LEN_TRIM(file_path))
58	write(,)'was not found by sugas_corrk.F90, exiting.'
59	call abort
60	endif
61
62	! check that database matches varactive toggle
63	open(111,file=file_path(1:LEN_TRIM(file_path)),form='formatted')
64	read(111,*) Nspecies
65
66	! if(Nspecies.eq.1 .and. (varactive.or.varfixed))then
67	! print*,'You have varactive/fixed=.true. but the database [', &
68	! corrkdir(1:LEN_TRIM(corrkdir)), &
69	! '] has no variable species, exiting.'
70	! call abort
71	! elseif(Nspecies.eq.2 .and. (.not.varactive) .and. (.not.varfixed))then
72	! print*,'You have varactive and varfixed=.false. and the database [', &
73	! corrkdir(1:LEN_TRIM(corrkdir)), &
74	! '] has a variable species.'
75	! call abort
76	! elseif(Nspecies.gt.3 .or. Nspecies.lt.1)then
77	! print*,Nspecies,' species in database [', &
78	! corrkdir(1:LEN_TRIM(corrkdir)), &
79	! '], radiative code cannot handle this.'
80	! call abort
81	! endif
82
83	do n=1,Nspecies
84	read(111,*) gastype(n)
85	print*,'Gas ',n,' is ',gastype(n)
86	enddo
87
88	! check the array size is correct, load the coefficients
89	open(111,file=file_path(1:LEN_TRIM(file_path)),form='formatted')
90	read(111,*) L_REFVARcheck
91	if(.not.(L_REFVARcheck.eq.L_REFVAR)) then
92	print*,'The size of your radiative transfer mixing ratio array does '
93	print*,'not match the value given in Q.dat, exiting.'
94	call abort
95	endif
96	read(111,*) wrefvar
97	close(111)
98
99	! display the values
100	print*,'Variable gas mixing ratios:'
101	do n=1,L_REFVAR
102	print*,n,'.',wrefvar(n),' kg/kg'
103	end do
104	print*,''
105
106
107	!=======================================================================
108	! Set the weighting in g-space
109
110	file_id='/corrk_data/' // corrkdir(1:LEN_TRIM(corrkdir)) // '/g.dat'
111	file_path=TRIM(datadir)//TRIM(file_id)
112
113	! check that the file exists
114	inquire(FILE=file_path,EXIST=file_ok)
115	if(.not.file_ok) then
116	write(,)'The file ',file_path(1:LEN_TRIM(file_path))
117	write(,)'was not found by sugas_corrk.F90, exiting.'
118	call abort
119	endif
120
121	! check the array size is correct, load the coefficients
122	open(111,file=file_path(1:LEN_TRIM(file_path)),form='formatted')
123	read(111,*) L_NGAUSScheck
124	if(.not.(L_NGAUSScheck.eq.L_NGAUSS)) then
125	print*,'The size of your radiative transfer g-space array does '
126	print*,'not match the value given in g.dat, exiting.'
127	call abort
128	endif
129	read(111,*) gweight
130	close(111)
131
132	! display the values
133	print*,'Correlated-k g-space grid:'
134	do n=1,L_NGAUSS
135	print*,n,'.',gweight(n)
136	end do
137	print*,''
138
139	!=======================================================================
140	! Set the reference pressure and temperature arrays. These are
141	! the pressures and temperatures at which we have k-coefficients.
142
143	!-----------------------------------------------------------------------
144	! pressure
145
146	file_id='/corrk_data/' // corrkdir(1:LEN_TRIM(corrkdir)) // '/p.dat'
147	file_path=TRIM(datadir)//TRIM(file_id)
148
149	! check that the file exists
150	inquire(FILE=file_path,EXIST=file_ok)
151	if(.not.file_ok) then
152	write(,)'The file ',file_path(1:LEN_TRIM(file_path))
153	write(,)'was not found by sugas_corrk.F90, exiting.'
154	call abort
155	endif
156
157	! check the array size is correct, load the coefficients
158	open(111,file=file_path(1:LEN_TRIM(file_path)),form='formatted')
159	read(111,*) L_NPREFcheck
160	if(.not.(L_NPREFcheck.eq.L_NPREF)) then
161	print*,'file_path = ',file_path, L_NPREFcheck,L_NPREF
162	print*,'check in radinc_h.F90'
163	print*,'The size of your radiative transfer pressure array does '
164	print*,'not match the value given in p.dat, exiting.'
165	call abort
166	endif
167	read(111,*) pgasref
168	close(111)
169
170	! display the values
171	print*,'Correlated-k pressure grid (mBar):'
172	do n=1,L_NPREF
173	print*,n,'. 1 x 10^',pgasref(n),' mBar'
174	end do
175	print*,''
176
177	! save the min / max matrix values
178	pgasmin = 10.0**pgasref(1)
179	pgasmax = 10.0**pgasref(L_NPREF)
180
181	! interpolate to finer grid
182	do n=1,L_NPREF-1
183	Xtemp=pgasref(n+1)-pgasref(n)
184	do m=1,5
185	! pfgasref((n-1)5+m) = pgasref(n)+(m-1)0.2
186
187	pfgasref((n-1)5+m) = pgasref(n)+(m-1)(Xtemp/5)
188
189	end do
190	end do
191	pfgasref(L_PINT) = pgasref(L_NPREF)
192	! Warning! this may need to be generalised if we want to use uneven grids!
193	! print*,'pgasref',pgasref
194	! print*,'pfgasref',pfgasref
195
196
197	!-----------------------------------------------------------------------
198	! temperature
199
200	file_id='/corrk_data/' // corrkdir(1:LEN_TRIM(corrkdir)) // '/T.dat'
201	file_path=TRIM(datadir)//TRIM(file_id)
202
203	! check that the file exists
204	inquire(FILE=file_path,EXIST=file_ok)
205	if(.not.file_ok) then
206	write(,)'The file ',file_path(1:LEN_TRIM(file_path))
207	write(,)'was not found by sugas_corrk.F90, exiting.'
208	call abort
209	endif
210
211	! check the array size is correct, load the coefficients
212	open(111,file=file_path(1:LEN_TRIM(file_path)),form='formatted')
213	read(111,*) L_NTREFcheck
214	if(.not.(L_NTREFcheck.eq.L_NTREF)) then
215	print*,'The size of your radiative transfer temperature array does '
216	print*,'not match the value given in T.dat, exiting.'
217	call abort
218	endif
219	read(111,*) tgasref
220	close(111)
221
222	! display the values
223	print*,'Correlated-k temperature grid:'
224	do n=1,L_NTREF
225	print*,n,'.',tgasref(n),' K'
226	end do
227
228	! save the min / max matrix values
229	tgasmin = tgasref(1)
230	tgasmax = tgasref(L_NTREF)
231
232	!=======================================================================
233	! Get gaseous k-coefficients and interpolate onto finer pressure grid
234
235	! VISIBLE
236	if (callgasvis) then
237	file_id='/corrk_data/'//trim(adjustl(banddir))//'/corrk_gcm_VI.dat'
238	file_path=TRIM(datadir)//TRIM(file_id)
239
240	! check that the file exists
241	inquire(FILE=file_path,EXIST=file_ok)
242	if(.not.file_ok) then
243	write(,)'The file ',file_path(1:LEN_TRIM(file_path))
244	write(,)'was not found by sugas_corrk.F90.'
245	write(,)'Are you sure you have absorption data for these bands?'
246	call abort
247	endif
248
249	open(111,file=file_path(1:LEN_TRIM(file_path)),form='formatted')
250	read(111,*) gasv8
251	close(111)
252
253	else
254	print*,'Visible gaseous absorption is set to zero.'
255	endif
256
257	! INFRA-RED
258	file_id='/corrk_data/'//trim(adjustl(banddir))//'/corrk_gcm_IR.dat'
259	file_path=TRIM(datadir)//TRIM(file_id)
260	! check that the file exists
261	inquire(FILE=file_path,EXIST=file_ok)
262	if(.not.file_ok) then
263	write(,)'The file ',file_path(1:LEN_TRIM(file_path))
264	write(,)'was not found by sugas_corrk.F90.'
265	write(,)'Are you sure you have absorption data for these bands?'
266	call abort
267	endif
268
269	open(111,file=file_path(1:LEN_TRIM(file_path)),form='formatted')
270	read(111,*) gasi8
271	close(111)
272
273	do nw=1,L_NSPECTI
274	fzeroi(nw) = 0.
275	end do
276	do nw=1,L_NSPECTV
277	fzerov(nw) = 0.
278	end do
279
280	! Take log10 of the values - this is what we will interpolate.
281	! Smallest value is 1.0E-200.
282
283	do nt=1,L_NTREF
284	do np=1,L_NPREF
285	do nh=1,L_REFVAR
286	do ng = 1,L_NGAUSS
287
288	do nw=1,L_NSPECTV
289	if(gasv8(nt,np,nh,nw,ng).gt.1.0d-200) then
290	gasv8(nt,np,nh,nw,ng) = &
291	log10(gasv8(nt,np,nh,nw,ng))
292	else
293	gasv8(nt,np,nh,nw,ng) = -200.0
294	end if
295	end do
296
297	do nw=1,L_NSPECTI
298	if(gasi8(nt,np,nh,nw,ng).gt.1.0d-200) then
299	gasi8(nt,np,nh,nw,ng) = &
300	log10(gasi8(nt,np,nh,nw,ng))
301	else
302	gasi8(nt,np,nh,nw,ng) = -200.0
303	end if
304	end do
305
306	end do
307	end do
308	end do
309	end do
310
311	! Interpolate the values: first the longwave
312
313	do nt=1,L_NTREF
314	do nh=1,L_REFVAR
315	do nw=1,L_NSPECTI
316	do ng=1,L_NGAUSS
317
318	! First, the initial interval
319
320	n = 1
321	do m=1,5
322	x = pfgasref(m)
323	xi(1) = pgasref(n)
324	xi(2) = pgasref(n+1)
325	xi(3) = pgasref(n+2)
326	xi(4) = pgasref(n+3)
327	yi(1) = gasi8(nt,n,nh,nw,ng)
328	yi(2) = gasi8(nt,n+1,nh,nw,ng)
329	yi(3) = gasi8(nt,n+2,nh,nw,ng)
330	yi(4) = gasi8(nt,n+3,nh,nw,ng)
331	call lagrange(x,xi,yi,ans)
332	gasi(nt,m,nh,nw,ng) = 10.0**ans
333	end do
334
335	do n=2,L_NPREF-2
336	do m=1,5
337	i = (n-1)*5+m
338	x = pfgasref(i)
339	xi(1) = pgasref(n-1)
340	xi(2) = pgasref(n)
341	xi(3) = pgasref(n+1)
342	xi(4) = pgasref(n+2)
343	yi(1) = gasi8(nt,n-1,nh,nw,ng)
344	yi(2) = gasi8(nt,n,nh,nw,ng)
345	yi(3) = gasi8(nt,n+1,nh,nw,ng)
346	yi(4) = gasi8(nt,n+2,nh,nw,ng)
347	call lagrange(x,xi,yi,ans)
348	gasi(nt,i,nh,nw,ng) = 10.0**ans
349	end do
350	end do
351
352	! Now, get the last interval
353
354	n = L_NPREF-1
355	do m=1,5
356	i = (n-1)*5+m
357	x = pfgasref(i)
358	xi(1) = pgasref(n-2)
359	xi(2) = pgasref(n-1)
360	xi(3) = pgasref(n)
361	xi(4) = pgasref(n+1)
362	yi(1) = gasi8(nt,n-2,nh,nw,ng)
363	yi(2) = gasi8(nt,n-1,nh,nw,ng)
364	yi(3) = gasi8(nt,n,nh,nw,ng)
365	yi(4) = gasi8(nt,n+1,nh,nw,ng)
366	call lagrange(x,xi,yi,ans)
367	gasi(nt,i,nh,nw,ng) = 10.0**ans
368	end do
369
370	! Fill the last pressure point
371
372	gasi(nt,L_PINT,nh,nw,ng) = &
373	10.0**gasi8(nt,L_NPREF,nh,nw,ng)
374
375	end do
376	end do
377	end do
378	end do
379
380	! Interpolate the values: now the shortwave
381
382	do nt=1,L_NTREF
383	do nh=1,L_REFVAR
384	do nw=1,L_NSPECTV
385	do ng=1,L_NGAUSS
386
387	! First, the initial interval
388
389	n = 1
390	do m=1,5
391	x = pfgasref(m)
392	xi(1) = pgasref(n)
393	xi(2) = pgasref(n+1)
394	xi(3) = pgasref(n+2)
395	xi(4) = pgasref(n+3)
396	yi(1) = gasv8(nt,n,nh,nw,ng)
397	yi(2) = gasv8(nt,n+1,nh,nw,ng)
398	yi(3) = gasv8(nt,n+2,nh,nw,ng)
399	yi(4) = gasv8(nt,n+3,nh,nw,ng)
400	call lagrange(x,xi,yi,ans)
401	gasv(nt,m,nh,nw,ng) = 10.0**ans
402	end do
403
404	do n=2,L_NPREF-2
405	do m=1,5
406	i = (n-1)*5+m
407	x = pfgasref(i)
408	xi(1) = pgasref(n-1)
409	xi(2) = pgasref(n)
410	xi(3) = pgasref(n+1)
411	xi(4) = pgasref(n+2)
412	yi(1) = gasv8(nt,n-1,nh,nw,ng)
413	yi(2) = gasv8(nt,n,nh,nw,ng)
414	yi(3) = gasv8(nt,n+1,nh,nw,ng)
415	yi(4) = gasv8(nt,n+2,nh,nw,ng)
416	call lagrange(x,xi,yi,ans)
417	gasv(nt,i,nh,nw,ng) = 10.0**ans
418	end do
419	end do
420
421	! Now, get the last interval
422
423	n = L_NPREF-1
424	do m=1,5
425	i = (n-1)*5+m
426	x = pfgasref(i)
427	xi(1) = pgasref(n-2)
428	xi(2) = pgasref(n-1)
429	xi(3) = pgasref(n)
430	xi(4) = pgasref(n+1)
431	yi(1) = gasv8(nt,n-2,nh,nw,ng)
432	yi(2) = gasv8(nt,n-1,nh,nw,ng)
433	yi(3) = gasv8(nt,n,nh,nw,ng)
434	yi(4) = gasv8(nt,n+1,nh,nw,ng)
435	call lagrange(x,xi,yi,ans)
436	gasv(nt,i,nh,nw,ng) = 10.0**ans
437	end do
438
439	! Fill the last pressure point
440
441	gasv(nt,L_PINT,nh,nw,ng) = &
442	10.0**gasv8(nt,L_NPREF,nh,nw,ng)
443
444	end do
445	end do
446	end do
447	end do
448
449	return
450	end subroutine sugas_corrk

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

Download in other formats:

Original Format