Context Navigation

sugas_corrk.F90 @ 3477

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

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

Download in other formats:

Original Format