Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

sugas_corrk.F90

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

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

Download in other formats:

Original Format