Context Navigation

sugas_corrk.F90 @ 3529

Last change on this file since 3529 was 2054, checked in by jvatant, 6 years ago
Minor correct a lacking condition --JVO
File size: 27.8 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
23	use radcommon_h, only : pgasref,pfgasref,pgasmin,pgasmax
24	use radcommon_h, only : tgasref,tgasmin,tgasmax
25	use radcommon_h, only : gasv,gasi,FZEROI,FZEROV,gweight, w_cum
26	use radcommon_h, only : WNOI,WNOV
27	use radcommon_h, only : gasi_recomb, gasv_recomb, pqrold, useptold
28	use radcommon_h, only : radvar_mask, radvar_indx, permut_idx
29	use datafile_mod, only: datadir, corrkdir, banddir
30	use comcstfi_mod, only: mugaz
31	use gases_h
32	use ioipsl_getin_p_mod, only: getin_p
33	use callkeys_mod, only: graybody,callgasvis, continuum, tracer, corrk_recombin
34	use tracer_h, only: noms, nqtot_p
35	implicit none
36
37	!==================================================================
38
39	logical file_ok
40
41	integer n, nt, np, nh, ng, nw, m, i
42
43	character(len=200) :: file_id
44	character(len=500) :: file_path
45
46	! ALLOCATABLE ARRAYS -- AS 12/2011
47	REAL*8, DIMENSION(:,:,:,:,:), ALLOCATABLE,SAVE :: gasi8, gasv8 !read by master
48	character*20,allocatable,DIMENSION(:),SAVE :: gastype ! for name check, read by master
49
50	real*8 x, xi(4), yi(4), ans, p
51	! For gray case (JL12)
52	real kappa_IR, kappa_VI, IR_VI_wnlimit
53	integer nVI_limit,nIR_limit
54
55	integer ngas, igas, jgas
56
57	double precision testcont ! for continuum absorption initialisation
58
59	integer :: dummy
60
61	integer ngasvar
62	logical found
63
64	!=======================================================================
65	! Load variable species data, exit if we have wrong database
66	file_id='/corrk_data/' // TRIM(corrkdir) // '/Q.dat'
67	file_path=TRIM(datadir)//TRIM(file_id)
68
69	! check that the file exists
70	inquire(FILE=file_path,EXIST=file_ok)
71	if(.not.file_ok) then
72	write(,)'The file ',TRIM(file_path)
73	write(,)'was not found by sugas_corrk.F90, exiting.'
74	write(,)'Check that your path to datagcm:',trim(datadir)
75	write(,)' is correct. You can change it in callphys.def with:'
76	write(,)' datadir = /absolute/path/to/datagcm'
77	write(,)'Also check that the corrkdir you chose in callphys.def exists.'
78	call abort
79	endif
80
81	!$OMP MASTER
82	! check that database matches varactive toggle
83	open(111,file=TRIM(file_path),form='formatted')
84	read(111,*) ngas
85
86	if(ngas.lt.1)then
87	print*,ngas,' species in database [', &
88	corrkdir(1:LEN_TRIM(corrkdir)), &
89	'], radiative code cannot handle this.'
90	call abort
91	endif
92
93	if ( corrk_recombin ) then
94	L_REFVAR = ngas
95
96	! dynamical allocations and read from Q.dat
97	IF ( .NOT. ALLOCATED( gastype ) ) ALLOCATE( gastype( L_REFVAR ) )
98	IF ( .NOT. ALLOCATED( radvar_mask ) ) ALLOCATE( radvar_mask( L_REFVAR ) )
99	IF ( .NOT. ALLOCATED( radvar_indx ) ) ALLOCATE( radvar_indx( L_REFVAR ) )
100
101	ngasvar = 0
102
103	write(,) ' '
104	do igas=1,L_REFVAR
105	read(111,*) gastype(igas)
106	read(111,*) radvar_mask(igas) ! boolean
107	print *, 'Found radiatively active specie : ',trim(gastype(igas)), ' : Variable =', radvar_mask(igas)
108	if ( radvar_mask(igas) ) ngasvar = ngasvar+1
109	enddo
110
111	print *, ' '
112	print *, 'TOTAL number of species found for radiative transfer = ', L_REFVAR
113	print *, ' including : ', ngasvar, 'variable species.'
114	print *, ' '
115
116	if (.not. tracer .and. ngasvar .gt. 0) then
117	write(,) 'sugas_corrk:error: You try to run variable species for corrk recombining'
118	write(,) 'in the radiative transfer but without tracers : I will stop here ! '
119	call abort
120	endif
121
122	! NB (JVO 18) We should enable here to force composition even if no associated tracers
123
124	! Check that we have matching tracers for the variable radiative active species
125	if ( ngasvar .GT. 0 ) then
126
127	do igas=1,L_REFVAR
128	if ( radvar_mask(igas) ) then
129
130	found = .false.
131	do jgas=1,nqtot_p
132	if ( trim(noms(jgas)) == trim(gastype(igas)) .OR. &
133	trim(noms(jgas)) == trim(gastype(igas))//'_') then
134	radvar_indx(igas) = jgas
135	found=.true.
136	exit
137	endif
138	enddo
139
140	if (.not. found) then
141	write(,) "sugas_corrk:error: "//trim(gastype(igas))//" is missing from tracers."
142	call abort
143	endif
144
145	endif
146	enddo ! igas=1,L_REFVAR
147
148	endif ! if ngasvar .gt. 0
149
150	else
151	L_REFVAR = 1 ! JVO 2017 : those will not be used, just set to 1 to keep the code running
152	endif ! if corrk_recombin
153
154	!=======================================================================
155	! Set the weighting in g-space
156
157	file_id='/corrk_data/' // TRIM(corrkdir) // '/g.dat'
158	file_path=TRIM(datadir)//TRIM(file_id)
159
160	! check that the file exists
161	inquire(FILE=file_path,EXIST=file_ok)
162	if(.not.file_ok) then
163	write(,)'The file ',TRIM(file_path)
164	write(,)'was not found by sugas_corrk.F90, exiting.'
165	write(,)'Check that your path to datagcm:',trim(datadir)
166	write(,)' is correct. You can change it in callphys.def with:'
167	write(,)' datadir = /absolute/path/to/datagcm'
168	write(,)'Also check that the corrkdir you chose in callphys.def exists.'
169	call abort
170	endif
171
172	! check the array size is correct, load the coefficients
173	open(111,file=TRIM(file_path),form='formatted')
174	read(111,*) L_NGAUSS
175	IF( .NOT. ALLOCATED( gweight ) ) ALLOCATE( GWEIGHT(L_NGAUSS) )
176	read(111,*) gweight
177	close(111)
178
179	! display the values
180	print*,'Correlated-k g-space grid:'
181	do n=1,L_NGAUSS
182	print*,n,'.',gweight(n)
183	end do
184	print*,''
185
186	IF (.NOT. ALLOCATED(permut_idx)) ALLOCATE(permut_idx(L_NGAUSS*L_NGAUSS)) ! for corr-k recombining
187	permut_idx = (/(i, i=1,L_NGAUSS*L_NGAUSS)/) ! for the recombin_corrk firstcall
188
189	IF (.NOT. ALLOCATED(w_cum)) ALLOCATE(w_cum(L_NGAUSS)) ! for corr-k recombining
190	w_cum(1)= gweight(1)
191	DO n=2,L_NGAUSS
192	w_cum(n) = w_cum(n-1)+gweight(n)
193	ENDDO
194
195	!=======================================================================
196	! Set the reference pressure and temperature arrays. These are
197	! the pressures and temperatures at which we have k-coefficients.
198
199	!-----------------------------------------------------------------------
200	! pressure
201
202	file_id='/corrk_data/' // TRIM(corrkdir) // '/p.dat'
203	file_path=TRIM(datadir)//TRIM(file_id)
204
205	! check that the file exists
206	inquire(FILE=file_path,EXIST=file_ok)
207	if(.not.file_ok) then
208	write(,)'The file ',TRIM(file_path)
209	write(,)'was not found by sugas_corrk.F90, exiting.'
210	write(,)'Check that your path to datagcm:',trim(datadir)
211	write(,)' is correct. You can change it in callphys.def with:'
212	write(,)' datadir = /absolute/path/to/datagcm'
213	write(,)'Also check that the corrkdir you chose in callphys.def exists.'
214	call abort
215	endif
216
217	! get array size, load the coefficients
218	open(111,file=TRIM(file_path),form='formatted')
219	read(111,*) L_NPREF
220	IF( .NOT. ALLOCATED( pgasref ) ) ALLOCATE( PGASREF(L_NPREF) )
221	read(111,*) pgasref
222	close(111)
223	L_PINT = (L_NPREF-1)*5+1
224	IF( .NOT. ALLOCATED( pfgasref ) ) ALLOCATE( PFGASREF(L_PINT) )
225
226	! display the values
227	print*,'Correlated-k pressure grid (mBar):'
228	do n=1,L_NPREF
229	print*,n,'. 1 x 10^',pgasref(n),' mBar'
230	end do
231	print*,''
232
233	! save the min / max matrix values
234	pgasmin = 10.0**pgasref(1)
235	pgasmax = 10.0**pgasref(L_NPREF)
236
237	! interpolate to finer grid, adapted to uneven grids
238	do n=1,L_NPREF-1
239	do m=1,5
240	pfgasref((n-1)5+m) = pgasref(n)+(m-1)(pgasref(n+1) - pgasref(n))/5.
241	end do
242	end do
243	pfgasref(L_PINT) = pgasref(L_NPREF)
244
245	!-----------------------------------------------------------------------
246	! temperature
247
248	file_id='/corrk_data/' // TRIM(corrkdir) // '/T.dat'
249	file_path=TRIM(datadir)//TRIM(file_id)
250
251	! check that the file exists
252	inquire(FILE=file_path,EXIST=file_ok)
253	if(.not.file_ok) then
254	write(,)'The file ',TRIM(file_path)
255	write(,)'was not found by sugas_corrk.F90, exiting.'
256	write(,)'Check that your path to datagcm:',trim(datadir)
257	write(,)' is correct. You can change it in callphys.def with:'
258	write(,)' datadir = /absolute/path/to/datagcm'
259	write(,)'Also check that the corrkdir you chose in callphys.def exists.'
260	call abort
261	endif
262
263	! get array size, load the coefficients
264	open(111,file=TRIM(file_path),form='formatted')
265	read(111,*) L_NTREF
266	IF( .NOT. ALLOCATED( tgasref ) ) ALLOCATE( TGASREF(L_NTREF) )
267	read(111,*) tgasref
268	close(111)
269
270	! display the values
271	print*,'Correlated-k temperature grid:'
272	do n=1,L_NTREF
273	print*,n,'.',tgasref(n),' K'
274	end do
275
276	! save the min / max matrix values
277	tgasmin = tgasref(1)
278	tgasmax = tgasref(L_NTREF)
279
280	IF( .NOT. ALLOCATED( gasi8 ) ) ALLOCATE( gasi8(L_NTREF,L_NPREF,L_REFVAR,L_NSPECTI,L_NGAUSS) )
281	IF( .NOT. ALLOCATED( gasv8 ) ) ALLOCATE( gasv8(L_NTREF,L_NPREF,L_REFVAR,L_NSPECTV,L_NGAUSS) )
282	!$OMP END MASTER
283	!$OMP BARRIER
284
285	!-----------------------------------------------------------------------
286	! allocate the multidimensional arrays in radcommon_h
287	IF( .NOT. ALLOCATED( gasi ) ) ALLOCATE( gasi(L_NTREF,L_PINT,L_REFVAR,L_NSPECTI,L_NGAUSS) )
288	IF( .NOT. ALLOCATED( gasv ) ) ALLOCATE( gasv(L_NTREF,L_PINT,L_REFVAR,L_NSPECTV,L_NGAUSS) )
289
290	if (corrk_recombin) then
291	IF( .NOT. ALLOCATED( gasi_recomb ) ) ALLOCATE( gasi_recomb(L_NTREF,L_PINT,L_NSPECTI,L_NGAUSS) )
292	IF( .NOT. ALLOCATED( gasv_recomb ) ) ALLOCATE( gasv_recomb(L_NTREF,L_PINT,L_NSPECTV,L_NGAUSS) )
293	IF( .NOT. ALLOCATED( pqrold ) ) ALLOCATE( pqrold(L_PINT,L_REFVAR) )
294	IF( .NOT. ALLOCATED( useptold ) ) ALLOCATE( useptold(L_PINT,L_NTREF) )
295	pqrold(:,:) = 0.0
296	useptold(:,:) = .false.
297	endif
298
299	! display the values
300	print*,''
301	print*,'Correlated-k matrix size:'
302	print*,'[',L_NTREF,',',L_NPREF,',',L_REFVAR,',',L_NGAUSS,']'
303
304	!=======================================================================
305	! Get gaseous k-coefficients and interpolate onto finer pressure grid
306
307
308	! wavelength used to separate IR from VI in graybody. We will need that anyway
309	IR_VI_wnlimit=3000.
310	write(,)"graybody: Visible / Infrared separation set at",10000./IR_VI_wnlimit,"um"
311
312	nVI_limit=0
313	Do nw=1,L_NSPECTV
314	if ((WNOV(nw).gt.IR_VI_wnlimit).and.(L_NSPECTV.gt.1)) then
315	nVI_limit=nw-1
316	exit
317	endif
318	End do
319	nIR_limit=L_NSPECTI
320	Do nw=1,L_NSPECTI
321	if ((WNOI(nw).gt.IR_VI_wnlimit).and.(L_NSPECTI.gt.1)) then
322	nIR_limit=nw-1
323	exit
324	endif
325	End do
326
327	if (graybody) then
328	! constant absorption coefficient in visible
329	write(,)"graybody: constant absorption coefficient in visible:"
330	kappa_VI=-100000.
331	call getin_p("kappa_VI",kappa_VI)
332	write(,)" kappa_VI = ",kappa_VI
333	kappa_VI=kappa_VI1.e4 mugaz * 1.672621e-27 ! conversion from m^2/kg to cm^2/molecule
334
335	! constant absorption coefficient in IR
336	write(,)"graybody: constant absorption coefficient in InfraRed:"
337	kappa_IR=-100000.
338	call getin_p("kappa_IR",kappa_IR)
339	write(,)" kappa_IR = ",kappa_IR
340	kappa_IR=kappa_IR1.e4 mugaz * 1.672621e-27 ! conversion from m^2/kg to cm^2/molecule
341
342	write(,)"graybody: Visible / Infrared separation set at band: IR=",nIR_limit,", VI=",nVI_limit
343
344	Else
345	kappa_VI=1.e-30
346	kappa_IR=1.e-30
347	End if
348
349	!$OMP MASTER
350	! print*,corrkdir(1:4)
351	! VISIBLE
352	if (callgasvis) then
353	if ((corrkdir(1:4).eq.'null'))then !(TRIM(corrkdir).eq.'null_LowTeffStar')) then
354	gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTV,1:L_NGAUSS)=0.0
355	print*,'using no corrk data'
356	print*,'Visible corrk gaseous absorption is set to zero if graybody=F'
357	else
358
359	if (corrk_recombin .eqv. .false. ) then
360	file_id='/corrk_data/'//trim(adjustl(banddir))//'/corrk_gcm_VI.dat'
361	file_path=TRIM(datadir)//TRIM(file_id)
362
363	! check that the file exists
364	inquire(FILE=file_path,EXIST=file_ok)
365	if(.not.file_ok) then
366	write(,)'The file ',TRIM(file_path)
367	write(,)'was not found by sugas_corrk.F90.'
368	write(,)'Are you sure you have absorption data for these bands?'
369	call abort
370	endif
371
372	open(111,file=TRIM(file_path),form='formatted')
373	read(111,*) gasv8
374	close(111)
375
376	else
377
378	do igas=1,L_REFVAR
379	! JVO 18 : To comply with Q.dat watch out for 2-letters gas trailing underscore :
380	! e.g. if ever you have N2_ (symmetric so shouldn't happen) file should be ..._N2_.dat
381	file_id='/corrk_data/'//trim(adjustl(banddir))//'/corrk_gcm_VI_'//trim(adjustl(gastype(igas)))//'.dat'
382	file_path=TRIM(datadir)//TRIM(file_id)
383
384	! check that the file exists
385	inquire(FILE=file_path,EXIST=file_ok)
386	if(.not.file_ok) then
387	write(,)'The file ',TRIM(file_path)
388	write(,)'was not found by sugas_corrk.F90.'
389	write(,)'Are you sure you have absorption data for these bands and this specie?'
390	call abort
391	endif
392
393	open(111,file=TRIM(file_path),form='formatted')
394	read(111,*) gasv8(:,:,igas,:,:)
395	close(111)
396	enddo
397
398	endif ! if corrk_recombin
399
400	end if
401
402	if(nVI_limit.eq.0) then
403	gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTV,1:L_NGAUSS)= &
404	gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTV,1:L_NGAUSS)+kappa_VI
405	else if (nVI_limit.eq.L_NSPECTV) then
406	gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTV,1:L_NGAUSS)= &
407	gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTV,1:L_NGAUSS)+kappa_IR
408	else
409	gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:nVI_limit,1:L_NGAUSS)= &
410	gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:nVI_limit,1:L_NGAUSS)+kappa_IR
411	gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,nVI_limit+1:L_NSPECTV,1:L_NGAUSS)= &
412	gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,nVI_limit+1:L_NSPECTV,1:L_NGAUSS)+kappa_VI
413	end if
414	else
415	print*,'Visible corrk gaseous absorption is set to zero.'
416	gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTV,1:L_NGAUSS)=0.0
417	endif
418	!$OMP END MASTER
419	!$OMP BARRIER
420
421	! INFRA-RED
422	if ((corrkdir(1:4).eq.'null'))then !.or.(TRIM(corrkdir).eq.'null_LowTeffStar')) then
423	print*,'Infrared corrk gaseous absorption is set to zero if graybody=F'
424	!$OMP MASTER
425	gasi8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTI,1:L_NGAUSS)=0.0
426	!$OMP END MASTER
427	!$OMP BARRIER
428	else
429
430	if ( corrk_recombin .eqv. .false. ) then
431	file_id='/corrk_data/'//trim(adjustl(banddir))//'/corrk_gcm_IR.dat'
432	file_path=TRIM(datadir)//TRIM(file_id)
433
434	! check that the file exists
435	inquire(FILE=file_path,EXIST=file_ok)
436	if(.not.file_ok) then
437	write(,)'The file ',TRIM(file_path)
438	write(,)'was not found by sugas_corrk.F90.'
439	write(,)'Are you sure you have absorption data for these bands?'
440	call abort
441	endif
442
443	!$OMP MASTER
444	open(111,file=TRIM(file_path),form='formatted')
445	read(111,*) gasi8
446	close(111)
447	!$OMP END MASTER
448	!$OMP BARRIER
449
450	else
451
452	do igas=1,L_REFVAR
453	! JVO 18 : To comply with Q.dat watch out for 2-letters gas trailing underscore :
454	! e.g. if ever tou have N2_ (symmetric so shouldn't happen) file should be ..._N2_.dat
455	file_id='/corrk_data/'//trim(adjustl(banddir))//'/corrk_gcm_IR_'//trim(adjustl(gastype(igas)))//'.dat'
456	file_path=TRIM(datadir)//TRIM(file_id)
457
458	! check that the file exists
459	inquire(FILE=file_path,EXIST=file_ok)
460	if(.not.file_ok) then
461	write(,)'The file ',TRIM(file_path)
462	write(,)'was not found by sugas_corrk.F90.'
463	write(,)'Are you sure you have absorption data for these bands and this specie?'
464	call abort
465	endif
466
467	!$OMP MASTER
468	open(111,file=TRIM(file_path),form='formatted')
469	read(111,*) gasi8(:,:,igas,:,:)
470	close(111)
471	!$OMP END MASTER
472	!$OMP BARRIER
473	enddo ! do igas=1,L_REFVAR
474
475	endif! if corrk_recombin
476
477
478	! 'fzero' is a currently unused feature that allows optimisation
479	! of the radiative transfer by neglecting bands where absorption
480	! is close to zero. As it could be useful in the future, this
481	! section of the code has been kept commented and not erased.
482	! RW 7/3/12.
483
484	do nw=1,L_NSPECTI
485	fzeroi(nw) = 0.d0
486	! do nt=1,L_NTREF
487	! do np=1,L_NPREF
488	! do nh=1,L_REFVAR
489	! do ng = 1,L_NGAUSS
490	! if(gasi8(nt,np,nh,nw,ng).lt.1.0e-25)then
491	! fzeroi(nw)=fzeroi(nw)+1.d0
492	! endif
493	! end do
494	! end do
495	! end do
496	! end do
497	! fzeroi(nw)=fzeroi(nw)/dble(L_NTREFL_NPREFL_REFVAR*L_NGAUSS)
498	end do
499
500	do nw=1,L_NSPECTV
501	fzerov(nw) = 0.d0
502	! do nt=1,L_NTREF
503	! do np=1,L_NPREF
504	! do nh=1,L_REFVAR
505	! do ng = 1,L_NGAUSS
506	! if(gasv8(nt,np,nh,nw,ng).lt.1.0e-25)then
507	! fzerov(nw)=fzerov(nw)+1.d0
508	! endif
509	! end do
510	! end do
511	! end do
512	! end do
513	! fzerov(nw)=fzerov(nw)/dble(L_NTREFL_NPREFL_REFVAR*L_NGAUSS)
514	end do
515
516	endif
517
518	!$OMP MASTER
519	if(nIR_limit.eq.0) then
520	gasi8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTI,1:L_NGAUSS)= &
521	gasi8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTI,1:L_NGAUSS)+kappa_VI
522	else if (nIR_limit.eq.L_NSPECTI) then
523	gasi8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTI,1:L_NGAUSS)= &
524	gasi8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTI,1:L_NGAUSS)+kappa_IR
525	else
526	gasi8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:nIR_limit,1:L_NGAUSS)= &
527	gasi8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:nIR_limit,1:L_NGAUSS)+kappa_IR
528	gasi8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,nIR_limit+1:L_NSPECTI,1:L_NGAUSS)= &
529	gasi8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,nIR_limit+1:L_NSPECTI,1:L_NGAUSS)+kappa_VI
530	end if
531
532
533	! Take log10 of the values - this is what we will interpolate.
534	! Smallest value is 1.0E-200.
535
536	do nt=1,L_NTREF
537	do np=1,L_NPREF
538	do nh=1,L_REFVAR
539	do ng = 1,L_NGAUSS
540
541	do nw=1,L_NSPECTV
542	if(gasv8(nt,np,nh,nw,ng).gt.1.0d-200) then
543	gasv8(nt,np,nh,nw,ng) = log10(gasv8(nt,np,nh,nw,ng))
544	else
545	gasv8(nt,np,nh,nw,ng) = -200.0
546	end if
547	end do
548
549	do nw=1,L_NSPECTI
550	if(gasi8(nt,np,nh,nw,ng).gt.1.0d-200) then
551	gasi8(nt,np,nh,nw,ng) = log10(gasi8(nt,np,nh,nw,ng))
552	else
553	gasi8(nt,np,nh,nw,ng) = -200.0
554	end if
555	end do
556
557	end do
558	end do
559	end do
560	end do
561	!$OMP END MASTER
562	!$OMP BARRIER
563
564	! Interpolate the values: first the longwave
565
566	do nt=1,L_NTREF
567	do nh=1,L_REFVAR
568	do nw=1,L_NSPECTI
569	do ng=1,L_NGAUSS
570
571	! First, the initial interval
572
573	n = 1
574	do m=1,5
575	x = pfgasref(m)
576	xi(1) = pgasref(n)
577	xi(2) = pgasref(n+1)
578	xi(3) = pgasref(n+2)
579	xi(4) = pgasref(n+3)
580	yi(1) = gasi8(nt,n,nh,nw,ng)
581	yi(2) = gasi8(nt,n+1,nh,nw,ng)
582	yi(3) = gasi8(nt,n+2,nh,nw,ng)
583	yi(4) = gasi8(nt,n+3,nh,nw,ng)
584	call lagrange(x,xi,yi,ans)
585	gasi(nt,m,nh,nw,ng) = 10.0**ans
586	end do
587
588	do n=2,L_NPREF-2
589	do m=1,5
590	i = (n-1)*5+m
591	x = pfgasref(i)
592	xi(1) = pgasref(n-1)
593	xi(2) = pgasref(n)
594	xi(3) = pgasref(n+1)
595	xi(4) = pgasref(n+2)
596	yi(1) = gasi8(nt,n-1,nh,nw,ng)
597	yi(2) = gasi8(nt,n,nh,nw,ng)
598	yi(3) = gasi8(nt,n+1,nh,nw,ng)
599	yi(4) = gasi8(nt,n+2,nh,nw,ng)
600	call lagrange(x,xi,yi,ans)
601	gasi(nt,i,nh,nw,ng) = 10.0**ans
602	end do
603	end do
604
605	! Now, get the last interval
606
607	n = L_NPREF-1
608	do m=1,5
609	i = (n-1)*5+m
610	x = pfgasref(i)
611	xi(1) = pgasref(n-2)
612	xi(2) = pgasref(n-1)
613	xi(3) = pgasref(n)
614	xi(4) = pgasref(n+1)
615	yi(1) = gasi8(nt,n-2,nh,nw,ng)
616	yi(2) = gasi8(nt,n-1,nh,nw,ng)
617	yi(3) = gasi8(nt,n,nh,nw,ng)
618	yi(4) = gasi8(nt,n+1,nh,nw,ng)
619	call lagrange(x,xi,yi,ans)
620	gasi(nt,i,nh,nw,ng) = 10.0**ans
621	end do
622
623	! Fill the last pressure point
624
625	gasi(nt,L_PINT,nh,nw,ng) = &
626	10.0**gasi8(nt,L_NPREF,nh,nw,ng)
627
628	end do
629	end do
630	end do
631	end do
632
633	! Interpolate the values: now the shortwave
634
635	do nt=1,L_NTREF
636	do nh=1,L_REFVAR
637	do nw=1,L_NSPECTV
638	do ng=1,L_NGAUSS
639
640	! First, the initial interval
641
642	n = 1
643	do m=1,5
644	x = pfgasref(m)
645	xi(1) = pgasref(n)
646	xi(2) = pgasref(n+1)
647	xi(3) = pgasref(n+2)
648	xi(4) = pgasref(n+3)
649	yi(1) = gasv8(nt,n,nh,nw,ng)
650	yi(2) = gasv8(nt,n+1,nh,nw,ng)
651	yi(3) = gasv8(nt,n+2,nh,nw,ng)
652	yi(4) = gasv8(nt,n+3,nh,nw,ng)
653	call lagrange(x,xi,yi,ans)
654	gasv(nt,m,nh,nw,ng) = 10.0**ans
655	end do
656
657	do n=2,L_NPREF-2
658	do m=1,5
659	i = (n-1)*5+m
660	x = pfgasref(i)
661	xi(1) = pgasref(n-1)
662	xi(2) = pgasref(n)
663	xi(3) = pgasref(n+1)
664	xi(4) = pgasref(n+2)
665	yi(1) = gasv8(nt,n-1,nh,nw,ng)
666	yi(2) = gasv8(nt,n,nh,nw,ng)
667	yi(3) = gasv8(nt,n+1,nh,nw,ng)
668	yi(4) = gasv8(nt,n+2,nh,nw,ng)
669	call lagrange(x,xi,yi,ans)
670	gasv(nt,i,nh,nw,ng) = 10.0**ans
671	end do
672	end do
673
674	! Now, get the last interval
675
676	n = L_NPREF-1
677	do m=1,5
678	i = (n-1)*5+m
679	x = pfgasref(i)
680	xi(1) = pgasref(n-2)
681	xi(2) = pgasref(n-1)
682	xi(3) = pgasref(n)
683	xi(4) = pgasref(n+1)
684	yi(1) = gasv8(nt,n-2,nh,nw,ng)
685	yi(2) = gasv8(nt,n-1,nh,nw,ng)
686	yi(3) = gasv8(nt,n,nh,nw,ng)
687	yi(4) = gasv8(nt,n+1,nh,nw,ng)
688	call lagrange(x,xi,yi,ans)
689	gasv(nt,i,nh,nw,ng) = 10.0**ans
690	end do
691
692	! Fill the last pressure point
693
694	gasv(nt,L_PINT,nh,nw,ng) = &
695	10.0**gasv8(nt,L_NPREF,nh,nw,ng)
696
697	end do
698	end do
699	end do
700	end do
701
702	if (corrk_recombin) then
703	gasi_recomb(:,:,:,:) = kappa_IR ! non-zero init
704	gasv_recomb(:,:,:,:) = kappa_VI ! non-zero init
705	endif
706
707	!=======================================================================
708	! Initialise the continuum absorption data
709	if(continuum)then
710	do igas=1,ngasmx
711
712	if (igas .eq. igas_N2) then
713
714	dummy = -9999
715	call interpolateN2N2(100.D+0,250.D+0,17500.D+0,testcont,.true.,dummy)
716
717	elseif (igas .eq. igas_H2) then
718
719	! first do self-induced absorption
720	dummy = -9999
721	call interpolateH2H2(500.D+0,250.D+0,17500.D+0,testcont,.true.,dummy)
722	! then cross-interactions with other gases
723	do jgas=1,ngasmx
724	if (jgas .eq. igas_N2) then
725	dummy = -9999
726	call interpolateN2H2(592.D+0,278.15D+0,200000.D+0,10000.D+0,testcont,.true.,dummy)
727	endif
728	enddo
729
730	elseif (igas .eq. igas_CH4) then
731
732	! first do self-induced absorption
733	dummy = -9999
734	call interpolateCH4CH4(200.D+0,200.D+0,7500.D+0,testcont,.true.,dummy)
735	! then cross-interactions with other gases
736	do jgas=1,ngasmx
737	if (jgas .eq. igas_N2) then
738	dummy = -9999
739	call interpolateN2CH4(200.D+0,250.0D+0,100000.D+0,5000.D+0,testcont,.true.,dummy)
740	endif
741	enddo
742
743	endif
744
745	enddo
746	endif
747
748	print*,'----------------------------------------------------'
749	print*,'And that`s all we have. It`s possible that other'
750	print*,'continuum absorption may be present, but if it is we'
751	print*,'don`t yet have data for it...'
752	print*,''
753
754	! Deallocate local arrays
755	!$OMP BARRIER
756	!$OMP MASTER
757	IF( ALLOCATED( gasi8 ) ) DEALLOCATE( gasi8 )
758	IF( ALLOCATED( gasv8 ) ) DEALLOCATE( gasv8 )
759	IF( ALLOCATED( pgasref ) ) DEALLOCATE( pgasref )
760	IF( ALLOCATED( gastype ) ) DEALLOCATE( gastype )
761	!$OMP END MASTER
762	!$OMP BARRIER
763
764	return
765	end subroutine sugas_corrk

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

Download in other formats:

Original Format