Context Navigation

sugas_corrk.F90 @ 2551

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

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

Download in other formats:

Original Format