Context Navigation

sugas_corrk.F90 @ 477

Last change on this file since 477 was 471, checked in by aslmd, 14 years ago

LMDZ.GENERIC

13/12/2011 == AS

Same spirit as previous commit, but for ngasmx which is now read in gases.def -- before arrays w/ dim ngasmx are allocated dynamically
Allocation is done in su_gases.F90 which is called in inifis
Outside su_gases.F90, very few modifications to the code : the new module "gases_h.F90" simply replaces the old common "gases.h" !
Compiles fine. Tested with debugging options through pgdbg. Runs fine. Exact same results in Early Mars test case.

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

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

Download in other formats:

Original Format