Context Navigation

sugas_corrk.F90 @ 3039

Last change on this file since 3039 was 2560, checked in by slebonnois, 4 years ago
SL: Implementation of SW computation based on generic model. Switch between this new SW module or old module that reads R. Haus tables implemented with a key (solarchoice)
File size: 18.2 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, only: L_NSPECTV, &
23	L_NGAUSS, L_NPREF, L_NTREF, L_REFVAR, L_PINT
24	use radcommon_h, only : pgasref,pfgasref,pgasmin,pgasmax
25	use radcommon_h, only : tgasref,tgasmin,tgasmax
26	use radcommon_h, only : gasv,FZEROV,gweight
27	use radcommon_h, only : WNOV
28	use datafile_mod, only: datadir,corrkdir,banddir
29	use gases_h, only: gnom, ngasmx, igas_H2O
30	implicit none
31	#include "YOMCST.h"
32	#include "clesphys.h"
33
34	!==================================================================
35
36	logical file_ok
37
38	integer n, nt, np, nh, ng, nw, m, i
39
40	character(len=200) :: file_id
41	character(len=500) :: file_path
42
43	! ALLOCATABLE ARRAYS -- AS 12/2011
44	REAL*8, DIMENSION(:,:,:,:,:), ALLOCATABLE,SAVE :: gasi8, gasv8 !read by master
45	character*20,allocatable,DIMENSION(:),SAVE :: gastype ! for check with gnom, read by master
46
47	real*8 x, xi(4), yi(4), ans, p
48	! For gray case (JL12)
49	real kappa_VI, kappa_IR, IR_VI_wnlimit
50	integer nVI_limit
51
52	integer ngas, igas, jgas
53
54	double precision testcont ! for continuum absorption initialisation
55
56	integer :: dummy
57	character(len=100) :: message
58	character(len=10),parameter :: subname="sugascorrk"
59
60	!=======================================================================
61	! Load variable species data, exit if we have wrong database
62	file_id='/corrk_data/' // TRIM(corrkdir) // '/Q.dat'
63	file_path=TRIM(datadir)//TRIM(file_id)
64
65	! check that the file exists
66	inquire(FILE=file_path,EXIST=file_ok)
67	if(.not.file_ok) then
68	write(,)'The file ',TRIM(file_path)
69	write(,)'was not found by sugas_corrk.F90, exiting.'
70	write(,)'Check that your path to datagcm:',trim(datadir)
71	write(,)' is correct. You can change it in callphys.def with:'
72	write(,)' datadir = /absolute/path/to/datagcm'
73	message='Also check that the corrkdir you chose in callphys.def exists.'
74	call abort_physic(subname,message,1)
75	endif
76
77	!$OMP MASTER
78	! check that database matches varactive toggle
79	open(111,file=TRIM(file_path),form='formatted')
80	read(111,*) ngas
81
82	if(ngas.ne.ngasmx)then
83	print*,'Number of gases in radiative transfer data (',ngas,') does not', &
84	'match that in gases_h (',ngasmx,')'
85	message='exiting.'
86	call abort_physic(subname,message,1)
87	endif
88
89	if(ngas.gt.5 .or. ngas.lt.1)then
90	print*,ngas,' species in database [', &
91	corrkdir(1:LEN_TRIM(corrkdir)),']'
92	message='radiative code cannot handle this.'
93	call abort_physic(subname,message,1)
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)then
112	print*,'The database [', &
113	corrkdir(1:LEN_TRIM(corrkdir)), &
114	'] has a variable species.'
115	message='Not possible with Venus...'
116	call abort_physic(subname,message,1)
117	endif
118
119	! Check that gastype and gnom match
120	do igas=1,ngas
121	print*,'Gas ',igas,' is ',trim(gnom(igas))
122	if (trim(gnom(igas)).ne.trim(gastype(igas))) then
123	print*,'Name of a gas in radiative transfer data (',trim(gastype(igas)),') does not ', &
124	'match that in gases.def (',trim(gnom(igas)),'). You should compare ', &
125	'gases.def with Q.dat in your radiative transfer directory.'
126	message='exiting.'
127	call abort_physic(subname,message,1)
128	endif
129	enddo
130	print*,'Confirmed gas match in radiative transfer and gases_h!'
131
132	! display the values
133	! print*,'Variable gas volume mixing ratios:'
134	! do n=1,L_REFVAR
135	!print*,n,'.',wrefvar(n),' kg/kg' ! pay attention!
136	! print*,n,'.',wrefvar(n),' mol/mol'
137	! end do
138	print*,''
139
140	!=======================================================================
141	! Set the weighting in g-space
142
143	file_id='/corrk_data/' // TRIM(corrkdir) // '/g.dat'
144	file_path=TRIM(datadir)//TRIM(file_id)
145
146	! check that the file exists
147	inquire(FILE=file_path,EXIST=file_ok)
148	if(.not.file_ok) then
149	write(,)'The file ',TRIM(file_path)
150	write(,)'was not found by sugas_corrk.F90, exiting.'
151	write(,)'Check that your path to datagcm:',trim(datadir)
152	write(,)' is correct. You can change it in callphys.def with:'
153	write(,)' datadir = /absolute/path/to/datagcm'
154	message='Also check that the corrkdir you chose in callphys.def exists.'
155	call abort_physic(subname,message,1)
156	endif
157
158	! check the array size is correct, load the coefficients
159	open(111,file=TRIM(file_path),form='formatted')
160	read(111,*) L_NGAUSS
161	IF( .NOT. ALLOCATED( gweight ) ) ALLOCATE( GWEIGHT(L_NGAUSS) )
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/' // TRIM(corrkdir) // '/p.dat'
180	file_path=TRIM(datadir)//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 ',TRIM(file_path)
186	write(,)'was not found by sugas_corrk.F90, exiting.'
187	write(,)'Check that your path to datagcm:',trim(datadir)
188	write(,)' is correct. You can change it in callphys.def with:'
189	write(,)' datadir = /absolute/path/to/datagcm'
190	message='Also check that the corrkdir you chose in callphys.def exists.'
191	call abort_physic(subname,message,1)
192	endif
193
194	! get array size, load the coefficients
195	open(111,file=TRIM(file_path),form='formatted')
196	read(111,*) L_NPREF
197	IF( .NOT. ALLOCATED( pgasref ) ) ALLOCATE( PGASREF(L_NPREF) )
198	read(111,*) pgasref
199	close(111)
200	L_PINT = (L_NPREF-1)*5+1
201	IF( .NOT. ALLOCATED( pfgasref ) ) ALLOCATE( PFGASREF(L_PINT) )
202
203	! display the values
204	print*,'Correlated-k pressure grid (mBar):'
205	do n=1,L_NPREF
206	print*,n,'. 1 x 10^',pgasref(n),' mBar'
207	end do
208	print*,''
209
210	! save the min / max matrix values
211	pgasmin = 10.0**pgasref(1)
212	pgasmax = 10.0**pgasref(L_NPREF)
213
214	! interpolate to finer grid, adapted to uneven grids
215	do n=1,L_NPREF-1
216	do m=1,5
217	pfgasref((n-1)5+m) = pgasref(n)+(m-1)(pgasref(n+1) - pgasref(n))/5.
218	end do
219	end do
220	pfgasref(L_PINT) = pgasref(L_NPREF)
221
222	!-----------------------------------------------------------------------
223	! temperature
224
225	file_id='/corrk_data/' // TRIM(corrkdir) // '/T.dat'
226	file_path=TRIM(datadir)//TRIM(file_id)
227
228	! check that the file exists
229	inquire(FILE=file_path,EXIST=file_ok)
230	if(.not.file_ok) then
231	write(,)'The file ',TRIM(file_path)
232	write(,)'was not found by sugas_corrk.F90, exiting.'
233	write(,)'Check that your path to datagcm:',trim(datadir)
234	write(,)' is correct. You can change it in callphys.def with:'
235	write(,)' datadir = /absolute/path/to/datagcm'
236	message='Also check that the corrkdir you chose in callphys.def exists.'
237	call abort_physic(subname,message,1)
238	endif
239
240	! get array size, load the coefficients
241	open(111,file=TRIM(file_path),form='formatted')
242	read(111,*) L_NTREF
243	IF( .NOT. ALLOCATED( tgasref ) ) ALLOCATE( TGASREF(L_NTREF) )
244	read(111,*) tgasref
245	close(111)
246
247	! display the values
248	print*,'Correlated-k temperature grid:'
249	do n=1,L_NTREF
250	print*,n,'.',tgasref(n),' K'
251	end do
252
253	! save the min / max matrix values
254	tgasmin = tgasref(1)
255	tgasmax = tgasref(L_NTREF)
256
257	IF( .NOT. ALLOCATED( gasv8 ) ) ALLOCATE( gasv8(L_NTREF,L_NPREF,L_REFVAR,L_NSPECTV,L_NGAUSS) )
258	!$OMP END MASTER
259	!$OMP BARRIER
260
261	!-----------------------------------------------------------------------
262	! allocate the multidimensional arrays in radcommon_h
263	IF( .NOT. ALLOCATED( gasv ) ) ALLOCATE( gasv(L_NTREF,L_PINT,L_REFVAR,L_NSPECTV,L_NGAUSS) )
264
265	! display the values
266	print*,''
267	print*,'Correlated-k matrix size:'
268	print*,'[',L_NTREF,',',L_NPREF,',',L_REFVAR,',',L_NGAUSS,']'
269
270	!=======================================================================
271	! Get gaseous k-coefficients and interpolate onto finer pressure grid
272
273
274	! wavelength used to separate IR from VI in graybody. We will need that anyway
275	IR_VI_wnlimit=3000.
276	! write(,)"graybody: Visible / Infrared separation set at",10000./IR_VI_wnlimit,"um"
277
278	nVI_limit=0
279	Do nw=1,L_NSPECTV
280	if ((WNOV(nw).gt.IR_VI_wnlimit).and.(L_NSPECTV.gt.1)) then
281	nVI_limit=nw-1
282	exit
283	endif
284	End do
285
286	kappa_VI=1.e-30
287	kappa_IR=1.e-30
288
289	!$OMP MASTER
290	! print*,corrkdir(1:4)
291	! VISIBLE
292	! if (callgasvis) then
293	if ((corrkdir(1:4).eq.'null'))then !(TRIM(corrkdir).eq.'null_LowTeffStar')) then
294	gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTV,1:L_NGAUSS)=0.0
295	print*,'using no corrk data'
296	print*,'Visible corrk gaseous absorption is set to zero if graybody=F'
297	else
298	file_id='/corrk_data/'//trim(adjustl(banddir))//'/corrk_gcm_VI.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	message='Are you sure you have absorption data for these bands?'
307	call abort_physic(subname,message,1)
308	endif
309
310	open(111,file=TRIM(file_path),form='formatted')
311	read(111,*) gasv8
312	close(111)
313	end if
314
315	if(nVI_limit.eq.0) then
316	gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTV,1:L_NGAUSS)= &
317	gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTV,1:L_NGAUSS)+kappa_VI
318	else if (nVI_limit.eq.L_NSPECTV) then
319	gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTV,1:L_NGAUSS)= &
320	gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTV,1:L_NGAUSS)+kappa_IR
321	else
322	gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:nVI_limit,1:L_NGAUSS)= &
323	gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:nVI_limit,1:L_NGAUSS)+kappa_IR
324	gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,nVI_limit+1:L_NSPECTV,1:L_NGAUSS)= &
325	gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,nVI_limit+1:L_NSPECTV,1:L_NGAUSS)+kappa_VI
326	end if
327	! else
328	! print*,'Visible corrk gaseous absorption is set to zero.'
329	! gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTV,1:L_NGAUSS)=0.0
330	! endif
331	!$OMP END MASTER
332	!$OMP BARRIER
333
334
335	! 'fzero' is a currently unused feature that allows optimisation
336	! of the radiative transfer by neglecting bands where absorption
337	! is close to zero. As it could be useful in the future, this
338	! section of the code has been kept commented and not erased.
339	! RW 7/3/12.
340
341	do nw=1,L_NSPECTV
342	fzerov(nw) = 0.d0
343	! do nt=1,L_NTREF
344	! do np=1,L_NPREF
345	! do nh=1,L_REFVAR
346	! do ng = 1,L_NGAUSS
347	! if(gasv8(nt,np,nh,nw,ng).lt.1.0e-25)then
348	! fzerov(nw)=fzerov(nw)+1.d0
349	! endif
350	! end do
351	! end do
352	! end do
353	! end do
354	! fzerov(nw)=fzerov(nw)/dble(L_NTREFL_NPREFL_REFVAR*L_NGAUSS)
355	end do
356
357	!$OMP MASTER
358
359	! Take log10 of the values - this is what we will interpolate.
360	! Smallest value is 1.0E-200.
361
362	do nt=1,L_NTREF
363	do np=1,L_NPREF
364	do nh=1,L_REFVAR
365	do ng = 1,L_NGAUSS
366
367	do nw=1,L_NSPECTV
368	if(gasv8(nt,np,nh,nw,ng).gt.1.0d-200) then
369	gasv8(nt,np,nh,nw,ng) = log10(gasv8(nt,np,nh,nw,ng))
370	else
371	gasv8(nt,np,nh,nw,ng) = -200.0
372	end if
373	end do
374
375	end do
376	end do
377	end do
378	end do
379	!$OMP END MASTER
380	!$OMP BARRIER
381
382	! Interpolate the values: now the shortwave
383
384	do nt=1,L_NTREF
385	do nh=1,L_REFVAR
386	do nw=1,L_NSPECTV
387	do ng=1,L_NGAUSS
388
389	! First, the initial interval
390
391	n = 1
392	do m=1,5
393	x = pfgasref(m)
394	xi(1) = pgasref(n)
395	xi(2) = pgasref(n+1)
396	xi(3) = pgasref(n+2)
397	xi(4) = pgasref(n+3)
398	yi(1) = gasv8(nt,n,nh,nw,ng)
399	yi(2) = gasv8(nt,n+1,nh,nw,ng)
400	yi(3) = gasv8(nt,n+2,nh,nw,ng)
401	yi(4) = gasv8(nt,n+3,nh,nw,ng)
402	call lagrange(x,xi,yi,ans)
403	gasv(nt,m,nh,nw,ng) = 10.0**ans
404	end do
405
406	do n=2,L_NPREF-2
407	do m=1,5
408	i = (n-1)*5+m
409	x = pfgasref(i)
410	xi(1) = pgasref(n-1)
411	xi(2) = pgasref(n)
412	xi(3) = pgasref(n+1)
413	xi(4) = pgasref(n+2)
414	yi(1) = gasv8(nt,n-1,nh,nw,ng)
415	yi(2) = gasv8(nt,n,nh,nw,ng)
416	yi(3) = gasv8(nt,n+1,nh,nw,ng)
417	yi(4) = gasv8(nt,n+2,nh,nw,ng)
418	call lagrange(x,xi,yi,ans)
419	gasv(nt,i,nh,nw,ng) = 10.0**ans
420	end do
421	end do
422
423	! Now, get the last interval
424
425	n = L_NPREF-1
426	do m=1,5
427	i = (n-1)*5+m
428	x = pfgasref(i)
429	xi(1) = pgasref(n-2)
430	xi(2) = pgasref(n-1)
431	xi(3) = pgasref(n)
432	xi(4) = pgasref(n+1)
433	yi(1) = gasv8(nt,n-2,nh,nw,ng)
434	yi(2) = gasv8(nt,n-1,nh,nw,ng)
435	yi(3) = gasv8(nt,n,nh,nw,ng)
436	yi(4) = gasv8(nt,n+1,nh,nw,ng)
437	call lagrange(x,xi,yi,ans)
438	gasv(nt,i,nh,nw,ng) = 10.0**ans
439	end do
440
441	! Fill the last pressure point
442
443	gasv(nt,L_PINT,nh,nw,ng) = &
444	10.0**gasv8(nt,L_NPREF,nh,nw,ng)
445
446	end do
447	end do
448	end do
449	end do
450
451
452	!=======================================================================
453	! Initialise the continuum absorption data
454	! if(continuum)then
455	do igas=1,ngasmx
456
457	! For Venus: only H2O, using CKD
458	if (igas .eq. igas_H2O) then
459
460	! H2O is special
461	! if(H2Ocont_simple)then
462	! call interpolateH2Ocont_PPC(990.D+0,296.D+0,683.2D+0*2,0.D+0,testcont,.true.)
463	! else
464	dummy = -9999
465	call interpolateH2Ocont_CKD(990.D+0,296.D+0,683.2D+0*2,0.D+0,testcont,.true.,dummy)
466	! endif
467
468	endif
469
470	enddo
471	! endif
472
473	! print*,'----------------------------------------------------'
474	! print*,'And that`s all we have. It`s possible that other'
475	! print*,'continuum absorption may be present, but if it is we'
476	! print*,'don`t yet have data for it...'
477	! print*,''
478
479	! Deallocate local arrays
480	!$OMP BARRIER
481	!$OMP MASTER
482	IF( ALLOCATED( gasv8 ) ) DEALLOCATE( gasv8 )
483	IF( ALLOCATED( pgasref ) ) DEALLOCATE( pgasref )
484	IF( ALLOCATED( gastype ) ) DEALLOCATE( gastype )
485	!$OMP END MASTER
486	!$OMP BARRIER
487
488	return
489	end subroutine sugas_corrk
490
491	!!===================================================================================
492	subroutine lagrange(x, xi, yi, ans)
493
494	! Lagrange interpolation - Polynomial interpolation at point x
495	! xi(1) <= x <= xi(4). Yi(n) is the functional value at XI(n).
496
497	implicit none
498
499	real*8 x, xi(4), yi(4), ans
500	real*8 fm1, fm2, fm3, fm4
501
502	!======================================================================!
503
504	fm1 = x - XI(1)
505	fm2 = x - XI(2)
506	fm3 = x - XI(3)
507	fm4 = x - XI(4)
508
509	! Get the answer at the requested X
510
511	ans = fm2fm3fm4*YI(1)/ &
512	((XI(1)-XI(2))(XI(1)-XI(3))(XI(1)-XI(4))) + &
513	fm1fm3fm4*YI(2)/ &
514	((XI(2)-XI(1))(XI(2)-XI(3))(XI(2)-XI(4))) + &
515	fm1fm2fm4*YI(3)/ &
516	((XI(3)-XI(1))(XI(3)-XI(2))(XI(3)-XI(4))) + &
517	fm1fm2fm3*YI(4)/ &
518	((XI(4)-XI(1))(XI(4)-XI(2))(XI(4)-XI(3)))
519
520	return
521	end

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

Download in other formats:

Original Format