Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

interpolate_continuum.F90 @ 3794

Last change on this file since 3794 was 3794, checked in by slebonnois, 3 days ago
SL: VenusPCM - update of CIA in solar module (adapted from Generic PCM) + potential output of band-resolved net flux
File size: 36.4 KB

Line
1	module interpolate_continuum_mod
2
3	implicit none
4
5	contains
6
7	subroutine interpolate_continuum(filename,igas_X,igas_Y,c_WN,ind_WN,temp,pres_X,pres_Y,abs_coef,firstcall)
8
9	!==================================================================
10	!
11	! Purpose
12	! -------
13	! Generic routine to calculate continuum opacities, using lookup tables provided here: https://web.lmd.jussieu.fr/~lmdz/planets/generic/datagcm/continuum_data/
14	! More information on the data here: https://lmdz-forge.lmd.jussieu.fr/mediawiki/Planets/index.php/Continuum_Database
15	!
16	! Author
17	! -------
18	! M. Turbet (2025)
19	!
20	!==================================================================
21
22	use datafile_mod, only: datadir
23	use mod_phys_lmdz_para, only : is_master
24
25	use gases_h, only: ngasmx, gnom, &
26	igas_H2, igas_H2O, igas_He, igas_N2, &
27	igas_CH4, igas_CO2, igas_O2
28
29	use radinc_h, only: L_NSPECTI, L_NSPECTV
30
31	use radcommon_h, only : BWNV,BWNI,WNOI,WNOV
32
33
34	implicit none
35
36	! input
37	integer,intent(in) :: ind_WN ! wavenumber index
38	integer,intent(in) :: igas_X ! index of molecule X
39	integer,intent(in) :: igas_Y ! index of molecule Y
40	double precision,intent(in) :: temp ! temperature (Kelvin)
41	double precision,intent(in) :: pres_X ! partial pressure of molecule X (Pascals)
42	double precision,intent(in) :: pres_Y ! partial pressure of molecule Y (Pascals)
43	character(len=*),intent(in) :: filename ! name of the lookup table
44	character(len=2),intent(in) :: c_WN ! wavelength chanel: infrared (IR) or visible (VI)
45	logical,intent(in) :: firstcall
46
47	! output
48	double precision,intent(out) :: abs_coef ! absorption coefficient (m^-1)
49
50	! intermediate variables
51	double precision amagat_X ! density of molecule X (in amagat units)
52	double precision amagat_Y ! density of molecule Y (in amagat units)
53
54	character(len=512) :: line
55
56	integer i, pos, iT, iW, iB, count_norm, igas
57
58	double precision temp_value, temp_abs, temp_wn
59
60	double precision z_temp
61
62	integer num_wn, num_T
63
64	double precision, dimension(:), allocatable :: temp_arr
65	double precision, dimension(:), allocatable :: wn_arr
66	double precision, dimension(:,:), allocatable :: abs_arr
67
68	integer ios
69
70	! Temperature array, continuum absorption grid for the pair N2-N2
71	integer,save :: num_T_N2N2
72	double precision,save,dimension(:),allocatable :: temp_arr_N2N2
73	double precision,save,dimension(:,:),allocatable :: abs_arr_N2N2_IR
74	double precision,save,dimension(:,:),allocatable :: abs_arr_N2N2_VI
75	! None of these saved variables are THREADPRIVATE because read by master
76	! and then only accessed but never modified and thus can be shared
77
78	! Temperature array, continuum absorption grid for the pair O2-O2
79	integer,save :: num_T_O2O2
80	double precision,save,dimension(:),allocatable :: temp_arr_O2O2
81	double precision,save,dimension(:,:),allocatable :: abs_arr_O2O2_IR
82	double precision,save,dimension(:,:),allocatable :: abs_arr_O2O2_VI
83	! None of these saved variables are THREADPRIVATE because read by master
84	! and then only accessed but never modified and thus can be shared
85
86	! Temperature array, continuum absorption grid for the pair H2-H2
87	integer,save :: num_T_H2H2
88	double precision,save,dimension(:),allocatable :: temp_arr_H2H2
89	double precision,save,dimension(:,:),allocatable :: abs_arr_H2H2_IR
90	double precision,save,dimension(:,:),allocatable :: abs_arr_H2H2_VI
91	! None of these saved variables are THREADPRIVATE because read by master
92	! and then only accessed but never modified and thus can be shared
93
94	! Temperature array, continuum absorption grid for the pair CO2-CO2
95	integer,save :: num_T_CO2CO2
96	double precision,save,dimension(:),allocatable :: temp_arr_CO2CO2
97	double precision,save,dimension(:,:),allocatable :: abs_arr_CO2CO2_IR
98	double precision,save,dimension(:,:),allocatable :: abs_arr_CO2CO2_VI
99	! None of these saved variables are THREADPRIVATE because read by master
100	! and then only accessed but never modified and thus can be shared
101
102	! Temperature array, continuum absorption grid for the pair CH4-CH4
103	integer,save :: num_T_CH4CH4
104	double precision,save,dimension(:),allocatable :: temp_arr_CH4CH4
105	double precision,save,dimension(:,:),allocatable :: abs_arr_CH4CH4_IR
106	double precision,save,dimension(:,:),allocatable :: abs_arr_CH4CH4_VI
107	! None of these saved variables are THREADPRIVATE because read by master
108	! and then only accessed but never modified and thus can be shared
109
110	! Temperature array, continuum absorption grid for the pair H2O-H2O
111	integer,save :: num_T_H2OH2O
112	double precision,save,dimension(:),allocatable :: temp_arr_H2OH2O
113	double precision,save,dimension(:,:),allocatable :: abs_arr_H2OH2O_IR
114	double precision,save,dimension(:,:),allocatable :: abs_arr_H2OH2O_VI
115	! None of these saved variables are THREADPRIVATE because read by master
116	! and then only accessed but never modified and thus can be shared
117
118	! Temperature array, continuum absorption grid for the pair H2-He
119	integer,save :: num_T_H2He
120	double precision,save,dimension(:),allocatable :: temp_arr_H2He
121	double precision,save,dimension(:,:),allocatable :: abs_arr_H2He_IR
122	double precision,save,dimension(:,:),allocatable :: abs_arr_H2He_VI
123	! None of these saved variables are THREADPRIVATE because read by master
124	! and then only accessed but never modified and thus can be shared
125
126	! Temperature array, continuum absorption grid for the pair H2-CH4
127	integer,save :: num_T_H2CH4
128	double precision,save,dimension(:),allocatable :: temp_arr_H2CH4
129	double precision,save,dimension(:,:),allocatable :: abs_arr_H2CH4_IR
130	double precision,save,dimension(:,:),allocatable :: abs_arr_H2CH4_VI
131	! None of these saved variables are THREADPRIVATE because read by master
132	! and then only accessed but never modified and thus can be shared
133
134	! Temperature array, continuum absorption grid for the pair CO2-H2
135	integer,save :: num_T_CO2H2
136	double precision,save,dimension(:),allocatable :: temp_arr_CO2H2
137	double precision,save,dimension(:,:),allocatable :: abs_arr_CO2H2_IR
138	double precision,save,dimension(:,:),allocatable :: abs_arr_CO2H2_VI
139	! None of these saved variables are THREADPRIVATE because read by master
140	! and then only accessed but never modified and thus can be shared
141
142	! Temperature array, continuum absorption grid for the pair CO2-CH4
143	integer,save :: num_T_CO2CH4
144	double precision,save,dimension(:),allocatable :: temp_arr_CO2CH4
145	double precision,save,dimension(:,:),allocatable :: abs_arr_CO2CH4_IR
146	double precision,save,dimension(:,:),allocatable :: abs_arr_CO2CH4_VI
147	! None of these saved variables are THREADPRIVATE because read by master
148	! and then only accessed but never modified and thus can be shared
149
150	! Temperature array, continuum absorption grid for the pair N2-H2
151	integer,save :: num_T_N2H2
152	double precision,save,dimension(:),allocatable :: temp_arr_N2H2
153	double precision,save,dimension(:,:),allocatable :: abs_arr_N2H2_IR
154	double precision,save,dimension(:,:),allocatable :: abs_arr_N2H2_VI
155	! None of these saved variables are THREADPRIVATE because read by master
156	! and then only accessed but never modified and thus can be shared
157
158	! Temperature array, continuum absorption grid for the pair N2-CH4
159	integer,save :: num_T_N2CH4
160	double precision,save,dimension(:),allocatable :: temp_arr_N2CH4
161	double precision,save,dimension(:,:),allocatable :: abs_arr_N2CH4_IR
162	double precision,save,dimension(:,:),allocatable :: abs_arr_N2CH4_VI
163	! None of these saved variables are THREADPRIVATE because read by master
164	! and then only accessed but never modified and thus can be shared
165
166	! Temperature array, continuum absorption grid for the pair CO2-O2
167	integer,save :: num_T_CO2O2
168	double precision,save,dimension(:),allocatable :: temp_arr_CO2O2
169	double precision,save,dimension(:,:),allocatable :: abs_arr_CO2O2_IR
170	double precision,save,dimension(:,:),allocatable :: abs_arr_CO2O2_VI
171	! None of these saved variables are THREADPRIVATE because read by master
172	! and then only accessed but never modified and thus can be shared
173
174	! Temperature array, continuum absorption grid for the pair N2-O2
175	integer,save :: num_T_N2O2
176	double precision,save,dimension(:), allocatable :: temp_arr_N2O2
177	double precision,save,dimension(:,:), allocatable :: abs_arr_N2O2_IR
178	double precision,save,dimension(:,:), allocatable :: abs_arr_N2O2_VI
179	! None of these saved variables are THREADPRIVATE because read by master
180	! and then only accessed but never modified and thus can be shared
181
182	! Temperature array, continuum absorption grid for the pair H2O-N2
183	integer,save :: num_T_H2ON2
184	double precision,save,dimension(:),allocatable :: temp_arr_H2ON2
185	double precision,save,dimension(:,:),allocatable :: abs_arr_H2ON2_IR
186	double precision,save,dimension(:,:),allocatable :: abs_arr_H2ON2_VI
187	! None of these saved variables are THREADPRIVATE because read by master
188	! and then only accessed but never modified and thus can be shared
189
190	! Temperature array, continuum absorption grid for the pair H2O-O2
191	integer,save :: num_T_H2OO2
192	double precision,save,dimension(:),allocatable :: temp_arr_H2OO2
193	double precision,save,dimension(:,:),allocatable :: abs_arr_H2OO2_IR
194	double precision,save,dimension(:,:),allocatable :: abs_arr_H2OO2_VI
195	! None of these saved variables are THREADPRIVATE because read by master
196	! and then only accessed but never modified and thus can be shared
197
198	! Temperature array, continuum absorption grid for the pair H2O-CO2
199	integer,save :: num_T_H2OCO2
200	double precision,save,dimension(:),allocatable :: temp_arr_H2OCO2
201	double precision,save,dimension(:,:),allocatable :: abs_arr_H2OCO2_IR
202	double precision,save,dimension(:,:),allocatable :: abs_arr_H2OCO2_VI
203	! None of these saved variables are THREADPRIVATE because read by master
204	! and then only accessed but never modified and thus can be shared
205
206
207	if(firstcall)then ! called by sugas_corrk only
208	if (is_master) print*,'----------------------------------------------------'
209	if (is_master) print*,'Initialising continuum (interpolate_continuum routine) from ', trim(filename)
210
211	!$OMP MASTER
212
213	open(unit=33, file=trim(filename), status="old", action="read",iostat=ios)
214
215	if (ios.ne.0) then ! file not found
216	if (is_master) then
217	write(,) 'Error from interpolate_continuum routine'
218	write(,) 'Data file ',trim(filename),' not found.'
219	write(,) 'Check that your path to datagcm:',trim(datadir)
220	write(,) 'is correct. You can change it in callphys.def with:'
221	write(,) 'datadir = /absolute/path/to/datagcm'
222	write(,) 'Also check that the continuum data is there.'
223	write(,) 'Latest continuum data can be downloaded here:'
224	write(,) 'https://web.lmd.jussieu.fr/~lmdz/planets/generic/datagcm/continuum_data/'
225	endif
226	call abort_physic("interpolate_continuum","missing input file",1)
227	endif
228
229	! We read the first line of the file to get the number of temperatures provided in the data file
230	read(33, '(A)') line
231
232	i = 1
233	iT = 0
234
235	do while (i .lt. len_trim(line))
236	pos = index(line(i:), 'T=')
237	if (pos == 0) exit
238	i = i + pos
239	iT = iT + 1
240	read(line(i+2:i+10), '(E9.2)') temp_value
241	end do
242
243	num_T=iT ! num_T is the number of temperatures provided in the data file
244
245	! We read all the remaining lines of the file to get the number of wavenumbers provided in the data file
246	iW = 0
247	do
248	read(33,*, end=501) line
249	iW = iW + 1
250	end do
251
252	501 continue
253
254	num_wn=iW ! num_wn is the number of wavenumbers provided in the data file
255
256	close(33)
257
258	allocate(temp_arr(num_T))
259	allocate(wn_arr(num_wn))
260	allocate(abs_arr(num_wn,num_T))
261
262	! We now open and read the file a second time to extract the temperature array, wavenumber array and continuum absorption data
263
264	open(unit=33, file=trim(filename), status="old", action="read")
265
266	! We extract the temperature array (temp_arr)
267
268	read(33, '(A)') line
269
270	i = 1
271	iT = 0
272
273	do while (i .lt. len_trim(line))
274	pos = index(line(i:), 'T=')
275	if (pos == 0) exit
276	i = i + pos
277	iT = iT + 1
278	read(line(i+2:i+10), '(E9.2)') temp_arr(iT)
279	end do
280
281	! We extract the wavenumber array (wn_arr) and continuum absorption (abs_arr)
282
283	do iW=1,num_wn
284	read(33,*) wn_arr(iW), (abs_arr(iW, iT), iT=1,num_T)
285	end do
286
287	close(33)
288
289	print*,'We read continuum absorption data for the pair ', trim(gnom(igas_X)),'-',trim(gnom(igas_Y))
290	print*,'Temperature grid of the dataset: ', temp_arr(:)
291
292	! We loop on all molecular pairs with available continuum data and fill the corresponding array
293
294	if ((igas_X .eq. igas_CO2) .and. (igas_Y .eq. igas_CO2)) then
295	num_T_CO2CO2=num_T
296	allocate(temp_arr_CO2CO2(num_T_CO2CO2))
297	allocate(abs_arr_CO2CO2_VI(L_NSPECTV,num_T_CO2CO2))
298	allocate(abs_arr_CO2CO2_IR(L_NSPECTI,num_T_CO2CO2))
299	temp_arr_CO2CO2(:)=temp_arr(:)
300	abs_arr_CO2CO2_VI(:,:)=0.
301	abs_arr_CO2CO2_IR(:,:)=0.
302	call interpolate_wn_abs_coeff(wn_arr,num_wn,abs_arr,abs_arr_CO2CO2_VI,abs_arr_CO2CO2_IR,num_T_CO2CO2)
303	elseif ((igas_X .eq. igas_N2) .and. (igas_Y .eq. igas_N2)) then
304	num_T_N2N2=num_T
305	allocate(temp_arr_N2N2(num_T_N2N2))
306	allocate(abs_arr_N2N2_VI(L_NSPECTV,num_T_N2N2))
307	allocate(abs_arr_N2N2_IR(L_NSPECTI,num_T_N2N2))
308	temp_arr_N2N2(:)=temp_arr(:)
309	abs_arr_N2N2_VI(:,:)=0.
310	abs_arr_N2N2_IR(:,:)=0.
311	call interpolate_wn_abs_coeff(wn_arr,num_wn,abs_arr,abs_arr_N2N2_VI,abs_arr_N2N2_IR,num_T_N2N2)
312	elseif ((igas_X .eq. igas_O2) .and. (igas_Y .eq. igas_O2)) then
313	num_T_O2O2=num_T
314	allocate(temp_arr_O2O2(num_T_O2O2))
315	allocate(abs_arr_O2O2_VI(L_NSPECTV,num_T_O2O2))
316	allocate(abs_arr_O2O2_IR(L_NSPECTI,num_T_O2O2))
317	temp_arr_O2O2(:)=temp_arr(:)
318	abs_arr_O2O2_VI(:,:)=0.
319	abs_arr_O2O2_IR(:,:)=0.
320	call interpolate_wn_abs_coeff(wn_arr,num_wn,abs_arr,abs_arr_O2O2_VI,abs_arr_O2O2_IR,num_T_O2O2)
321	elseif ((igas_X .eq. igas_CH4) .and. (igas_Y .eq. igas_CH4)) then
322	num_T_CH4CH4=num_T
323	allocate(temp_arr_CH4CH4(num_T_CH4CH4))
324	allocate(abs_arr_CH4CH4_VI(L_NSPECTV,num_T_CH4CH4))
325	allocate(abs_arr_CH4CH4_IR(L_NSPECTI,num_T_CH4CH4))
326	temp_arr_CH4CH4(:)=temp_arr(:)
327	abs_arr_CH4CH4_VI(:,:)=0.
328	abs_arr_CH4CH4_IR(:,:)=0.
329	call interpolate_wn_abs_coeff(wn_arr,num_wn,abs_arr,abs_arr_CH4CH4_VI,abs_arr_CH4CH4_IR,num_T_CH4CH4)
330	elseif ((igas_X .eq. igas_H2) .and. (igas_Y .eq. igas_H2)) then
331	num_T_H2H2=num_T
332	allocate(temp_arr_H2H2(num_T_H2H2))
333	allocate(abs_arr_H2H2_VI(L_NSPECTV,num_T_H2H2))
334	allocate(abs_arr_H2H2_IR(L_NSPECTI,num_T_H2H2))
335	temp_arr_H2H2(:)=temp_arr(:)
336	abs_arr_H2H2_VI(:,:)=0.
337	abs_arr_H2H2_IR(:,:)=0.
338	call interpolate_wn_abs_coeff(wn_arr,num_wn,abs_arr,abs_arr_H2H2_VI,abs_arr_H2H2_IR,num_T_H2H2)
339	elseif ((igas_X .eq. igas_H2O) .and. (igas_Y .eq. igas_H2O)) then
340	num_T_H2OH2O=num_T
341	allocate(temp_arr_H2OH2O(num_T_H2OH2O))
342	allocate(abs_arr_H2OH2O_VI(L_NSPECTV,num_T_H2OH2O))
343	allocate(abs_arr_H2OH2O_IR(L_NSPECTI,num_T_H2OH2O))
344	temp_arr_H2OH2O(:)=temp_arr(:)
345	abs_arr_H2OH2O_VI(:,:)=0.
346	abs_arr_H2OH2O_IR(:,:)=0.
347	call interpolate_wn_abs_coeff(wn_arr,num_wn,abs_arr,abs_arr_H2OH2O_VI,abs_arr_H2OH2O_IR,num_T_H2OH2O)
348	elseif ((igas_X .eq. igas_N2) .and. (igas_Y .eq. igas_H2)) then
349	num_T_N2H2=num_T
350	allocate(temp_arr_N2H2(num_T_N2H2))
351	allocate(abs_arr_N2H2_VI(L_NSPECTV,num_T_N2H2))
352	allocate(abs_arr_N2H2_IR(L_NSPECTI,num_T_N2H2))
353	temp_arr_N2H2(:)=temp_arr(:)
354	abs_arr_N2H2_VI(:,:)=0.
355	abs_arr_N2H2_IR(:,:)=0.
356	call interpolate_wn_abs_coeff(wn_arr,num_wn,abs_arr,abs_arr_N2H2_VI,abs_arr_N2H2_IR,num_T_N2H2)
357	elseif ((igas_X .eq. igas_N2) .and. (igas_Y .eq. igas_O2)) then
358	num_T_N2O2=num_T
359	allocate(temp_arr_N2O2(num_T_N2O2))
360	allocate(abs_arr_N2O2_VI(L_NSPECTV,num_T_N2O2))
361	allocate(abs_arr_N2O2_IR(L_NSPECTI,num_T_N2O2))
362	temp_arr_N2O2(:)=temp_arr(:)
363	abs_arr_N2O2_VI(:,:)=0.
364	abs_arr_N2O2_IR(:,:)=0.
365	call interpolate_wn_abs_coeff(wn_arr,num_wn,abs_arr,abs_arr_N2O2_VI,abs_arr_N2O2_IR,num_T_N2O2)
366	elseif ((igas_X .eq. igas_N2) .and. (igas_Y .eq. igas_CH4)) then
367	num_T_N2CH4=num_T
368	allocate(temp_arr_N2CH4(num_T_N2CH4))
369	allocate(abs_arr_N2CH4_VI(L_NSPECTV,num_T_N2CH4))
370	allocate(abs_arr_N2CH4_IR(L_NSPECTI,num_T_N2CH4))
371	temp_arr_N2CH4(:)=temp_arr(:)
372	abs_arr_N2CH4_VI(:,:)=0.
373	abs_arr_N2CH4_IR(:,:)=0.
374	call interpolate_wn_abs_coeff(wn_arr,num_wn,abs_arr,abs_arr_N2CH4_VI,abs_arr_N2CH4_IR,num_T_N2CH4)
375	elseif ((igas_X .eq. igas_CO2) .and. (igas_Y .eq. igas_O2)) then
376	num_T_CO2O2=num_T
377	allocate(temp_arr_CO2O2(num_T_CO2O2))
378	allocate(abs_arr_CO2O2_VI(L_NSPECTV,num_T_CO2O2))
379	allocate(abs_arr_CO2O2_IR(L_NSPECTI,num_T_CO2O2))
380	temp_arr_CO2O2(:)=temp_arr(:)
381	abs_arr_CO2O2_VI(:,:)=0.
382	abs_arr_CO2O2_IR(:,:)=0.
383	call interpolate_wn_abs_coeff(wn_arr,num_wn,abs_arr,abs_arr_CO2O2_VI,abs_arr_CO2O2_IR,num_T_CO2O2)
384	elseif ((igas_X .eq. igas_H2) .and. (igas_Y .eq. igas_CH4)) then
385	num_T_H2CH4=num_T
386	allocate(temp_arr_H2CH4(num_T_H2CH4))
387	allocate(abs_arr_H2CH4_VI(L_NSPECTV,num_T_H2CH4))
388	allocate(abs_arr_H2CH4_IR(L_NSPECTI,num_T_H2CH4))
389	temp_arr_H2CH4(:)=temp_arr(:)
390	abs_arr_H2CH4_VI(:,:)=0.
391	abs_arr_H2CH4_IR(:,:)=0.
392	call interpolate_wn_abs_coeff(wn_arr,num_wn,abs_arr,abs_arr_H2CH4_VI,abs_arr_H2CH4_IR,num_T_H2CH4)
393	elseif ((igas_X .eq. igas_H2) .and. (igas_Y .eq. igas_He)) then
394	num_T_H2He=num_T
395	allocate(temp_arr_H2He(num_T_H2He))
396	allocate(abs_arr_H2He_VI(L_NSPECTV,num_T_H2He))
397	allocate(abs_arr_H2He_IR(L_NSPECTI,num_T_H2He))
398	temp_arr_H2He(:)=temp_arr(:)
399	abs_arr_H2He_VI(:,:)=0.
400	abs_arr_H2He_IR(:,:)=0.
401	call interpolate_wn_abs_coeff(wn_arr,num_wn,abs_arr,abs_arr_H2He_VI,abs_arr_H2He_IR,num_T_H2He)
402	elseif ((igas_X .eq. igas_H2O) .and. (igas_Y .eq. igas_N2)) then
403	num_T_H2ON2=num_T
404	allocate(temp_arr_H2ON2(num_T_H2ON2))
405	allocate(abs_arr_H2ON2_VI(L_NSPECTV,num_T_H2ON2))
406	allocate(abs_arr_H2ON2_IR(L_NSPECTI,num_T_H2ON2))
407	temp_arr_H2ON2(:)=temp_arr(:)
408	abs_arr_H2ON2_VI(:,:)=0.
409	abs_arr_H2ON2_IR(:,:)=0.
410	call interpolate_wn_abs_coeff(wn_arr,num_wn,abs_arr,abs_arr_H2ON2_VI,abs_arr_H2ON2_IR,num_T_H2ON2)
411	elseif ((igas_X .eq. igas_H2O) .and. (igas_Y .eq. igas_O2)) then
412	num_T_H2OO2=num_T
413	allocate(temp_arr_H2OO2(num_T_H2OO2))
414	allocate(abs_arr_H2OO2_VI(L_NSPECTV,num_T_H2OO2))
415	allocate(abs_arr_H2OO2_IR(L_NSPECTI,num_T_H2OO2))
416	temp_arr_H2OO2(:)=temp_arr(:)
417	abs_arr_H2OO2_VI(:,:)=0.
418	abs_arr_H2OO2_IR(:,:)=0.
419	call interpolate_wn_abs_coeff(wn_arr,num_wn,abs_arr,abs_arr_H2OO2_VI,abs_arr_H2OO2_IR,num_T_H2OO2)
420	elseif ((igas_X .eq. igas_H2O) .and. (igas_Y .eq. igas_CO2)) then
421	num_T_H2OCO2=num_T
422	allocate(temp_arr_H2OCO2(num_T_H2OCO2))
423	allocate(abs_arr_H2OCO2_VI(L_NSPECTV,num_T_H2OCO2))
424	allocate(abs_arr_H2OCO2_IR(L_NSPECTI,num_T_H2OCO2))
425	temp_arr_H2OCO2(:)=temp_arr(:)
426	abs_arr_H2OCO2_VI(:,:)=0.
427	abs_arr_H2OCO2_IR(:,:)=0.
428	call interpolate_wn_abs_coeff(wn_arr,num_wn,abs_arr,abs_arr_H2OCO2_VI,abs_arr_H2OCO2_IR,num_T_H2OCO2)
429	elseif ((igas_X .eq. igas_CO2) .and. (igas_Y .eq. igas_CO2)) then
430	num_T_CO2CO2=num_T
431	allocate(temp_arr_CO2CO2(num_T_CO2CO2))
432	allocate(abs_arr_CO2CO2_VI(L_NSPECTV,num_T_CO2CO2))
433	allocate(abs_arr_CO2CO2_IR(L_NSPECTI,num_T_CO2CO2))
434	temp_arr_CO2CO2(:)=temp_arr(:)
435	abs_arr_CO2CO2_VI(:,:)=0.
436	abs_arr_CO2CO2_IR(:,:)=0.
437	call interpolate_wn_abs_coeff(wn_arr,num_wn,abs_arr,abs_arr_CO2CO2_VI,abs_arr_CO2CO2_IR,num_T_CO2CO2)
438	elseif ((igas_X .eq. igas_CO2) .and. (igas_Y .eq. igas_H2)) then
439	num_T_CO2H2=num_T
440	allocate(temp_arr_CO2H2(num_T_CO2H2))
441	allocate(abs_arr_CO2H2_VI(L_NSPECTV,num_T_CO2H2))
442	allocate(abs_arr_CO2H2_IR(L_NSPECTI,num_T_CO2H2))
443	temp_arr_CO2H2(:)=temp_arr(:)
444	abs_arr_CO2H2_VI(:,:)=0.
445	abs_arr_CO2H2_IR(:,:)=0.
446	call interpolate_wn_abs_coeff(wn_arr,num_wn,abs_arr,abs_arr_CO2H2_VI,abs_arr_CO2H2_IR,num_T_CO2H2)
447	elseif ((igas_X .eq. igas_CO2) .and. (igas_Y .eq. igas_CH4)) then
448	num_T_CO2CH4=num_T
449	allocate(temp_arr_CO2CH4(num_T_CO2CH4))
450	allocate(abs_arr_CO2CH4_VI(L_NSPECTV,num_T_CO2CH4))
451	allocate(abs_arr_CO2CH4_IR(L_NSPECTI,num_T_CO2CH4))
452	temp_arr_CO2CH4(:)=temp_arr(:)
453	abs_arr_CO2CH4_VI(:,:)=0.
454	abs_arr_CO2CH4_IR(:,:)=0.
455	call interpolate_wn_abs_coeff(wn_arr,num_wn,abs_arr,abs_arr_CO2CH4_VI,abs_arr_CO2CH4_IR,num_T_CO2CH4)
456	endif ! igas_X / igas_Y condition
457
458
459	!$OMP END MASTER
460	!$OMP BARRIER
461
462
463	endif ! firstcall
464
465	! We loop on all molecular pairs with available continuum data and interpolate in the temperature field
466	! Two options: we call visible (VI) or infrared (IR) tables, depending on the value of c_WN
467
468	if ((igas_X .eq. igas_CO2) .and. (igas_Y .eq. igas_CO2)) then
469	call T_boundaries_continuum(z_temp,temp,temp_arr_CO2CO2,num_T_CO2CO2)
470	if(c_WN .eq. 'IR') then
471	call interpolate_T_abs_coeff(z_temp,temp_arr_CO2CO2,num_T_CO2CO2,abs_coef,abs_arr_CO2CO2_IR(ind_WN,:))
472	elseif(c_WN .eq. 'VI') then
473	call interpolate_T_abs_coeff(z_temp,temp_arr_CO2CO2,num_T_CO2CO2,abs_coef,abs_arr_CO2CO2_VI(ind_WN,:))
474	else
475	print*,'You must select visible (VI) or infrared (IR) canal.'
476	stop
477	endif
478	elseif ((igas_X .eq. igas_N2) .and. (igas_Y .eq. igas_N2)) then
479	call T_boundaries_continuum(z_temp,temp,temp_arr_N2N2,num_T_N2N2)
480	if(c_WN .eq. 'IR') then
481	call interpolate_T_abs_coeff(z_temp,temp_arr_N2N2,num_T_N2N2,abs_coef,abs_arr_N2N2_IR(ind_WN,:))
482	elseif(c_WN .eq. 'VI') then
483	call interpolate_T_abs_coeff(z_temp,temp_arr_N2N2,num_T_N2N2,abs_coef,abs_arr_N2N2_VI(ind_WN,:))
484	else
485	print*,'You must select visible (VI) or infrared (IR) canal.'
486	stop
487	endif
488	elseif ((igas_X .eq. igas_O2) .and. (igas_Y .eq. igas_O2)) then
489	call T_boundaries_continuum(z_temp,temp,temp_arr_O2O2,num_T_O2O2)
490	if(c_WN .eq. 'IR') then
491	call interpolate_T_abs_coeff(z_temp,temp_arr_O2O2,num_T_O2O2,abs_coef,abs_arr_O2O2_IR(ind_WN,:))
492	elseif(c_WN .eq. 'VI') then
493	call interpolate_T_abs_coeff(z_temp,temp_arr_O2O2,num_T_O2O2,abs_coef,abs_arr_O2O2_VI(ind_WN,:))
494	else
495	print*,'You must select visible (VI) or infrared (IR) canal.'
496	stop
497	endif
498	elseif ((igas_X .eq. igas_CH4) .and. (igas_Y .eq. igas_CH4)) then
499	call T_boundaries_continuum(z_temp,temp,temp_arr_CH4CH4,num_T_CH4CH4)
500	if(c_WN .eq. 'IR') then
501	call interpolate_T_abs_coeff(z_temp,temp_arr_CH4CH4,num_T_CH4CH4,abs_coef,abs_arr_CH4CH4_IR(ind_WN,:))
502	elseif(c_WN .eq. 'VI') then
503	call interpolate_T_abs_coeff(z_temp,temp_arr_CH4CH4,num_T_CH4CH4,abs_coef,abs_arr_CH4CH4_VI(ind_WN,:))
504	else
505	print*,'You must select visible (VI) or infrared (IR) canal.'
506	stop
507	endif
508	elseif ((igas_X .eq. igas_H2) .and. (igas_Y .eq. igas_H2)) then
509	call T_boundaries_continuum(z_temp,temp,temp_arr_H2H2,num_T_H2H2)
510	if(c_WN .eq. 'IR') then
511	call interpolate_T_abs_coeff(z_temp,temp_arr_H2H2,num_T_H2H2,abs_coef,abs_arr_H2H2_IR(ind_WN,:))
512	elseif(c_WN .eq. 'VI') then
513	call interpolate_T_abs_coeff(z_temp,temp_arr_H2H2,num_T_H2H2,abs_coef,abs_arr_H2H2_VI(ind_WN,:))
514	else
515	print*,'You must select visible (VI) or infrared (IR) canal.'
516	stop
517	endif
518	elseif ((igas_X .eq. igas_H2O) .and. (igas_Y .eq. igas_H2O)) then
519	call T_boundaries_continuum(z_temp,temp,temp_arr_H2OH2O,num_T_H2OH2O)
520	if(c_WN .eq. 'IR') then
521	call interpolate_T_abs_coeff(z_temp,temp_arr_H2OH2O,num_T_H2OH2O,abs_coef,abs_arr_H2OH2O_IR(ind_WN,:))
522	elseif(c_WN .eq. 'VI') then
523	call interpolate_T_abs_coeff(z_temp,temp_arr_H2OH2O,num_T_H2OH2O,abs_coef,abs_arr_H2OH2O_VI(ind_WN,:))
524	else
525	print*,'You must select visible (VI) or infrared (IR) canal.'
526	stop
527	endif
528	elseif ((igas_X .eq. igas_N2) .and. (igas_Y .eq. igas_H2)) then
529	call T_boundaries_continuum(z_temp,temp,temp_arr_N2H2,num_T_N2H2)
530	if(c_WN .eq. 'IR') then
531	call interpolate_T_abs_coeff(z_temp,temp_arr_N2H2,num_T_N2H2,abs_coef,abs_arr_N2H2_IR(ind_WN,:))
532	elseif(c_WN .eq. 'VI') then
533	call interpolate_T_abs_coeff(z_temp,temp_arr_N2H2,num_T_N2H2,abs_coef,abs_arr_N2H2_VI(ind_WN,:))
534	else
535	print*,'You must select visible (VI) or infrared (IR) canal.'
536	stop
537	endif
538	elseif ((igas_X .eq. igas_N2) .and. (igas_Y .eq. igas_O2)) then
539	call T_boundaries_continuum(z_temp,temp,temp_arr_N2O2,num_T_N2O2)
540	if(c_WN .eq. 'IR') then
541	call interpolate_T_abs_coeff(z_temp,temp_arr_N2O2,num_T_N2O2,abs_coef,abs_arr_N2O2_IR(ind_WN,:))
542	elseif(c_WN .eq. 'VI') then
543	call interpolate_T_abs_coeff(z_temp,temp_arr_N2O2,num_T_N2O2,abs_coef,abs_arr_N2O2_VI(ind_WN,:))
544	else
545	print*,'You must select visible (VI) or infrared (IR) canal.'
546	stop
547	endif
548	elseif ((igas_X .eq. igas_N2) .and. (igas_Y .eq. igas_CH4)) then
549	call T_boundaries_continuum(z_temp,temp,temp_arr_N2CH4,num_T_N2CH4)
550	if(c_WN .eq. 'IR') then
551	call interpolate_T_abs_coeff(z_temp,temp_arr_N2CH4,num_T_N2CH4,abs_coef,abs_arr_N2CH4_IR(ind_WN,:))
552	elseif(c_WN .eq. 'VI') then
553	call interpolate_T_abs_coeff(z_temp,temp_arr_N2CH4,num_T_N2CH4,abs_coef,abs_arr_N2CH4_VI(ind_WN,:))
554	else
555	print*,'You must select visible (VI) or infrared (IR) canal.'
556	stop
557	endif
558	elseif ((igas_X .eq. igas_CO2) .and. (igas_Y .eq. igas_O2)) then
559	call T_boundaries_continuum(z_temp,temp,temp_arr_CO2O2,num_T_CO2O2)
560	if(c_WN .eq. 'IR') then
561	call interpolate_T_abs_coeff(z_temp,temp_arr_CO2O2,num_T_CO2O2,abs_coef,abs_arr_CO2O2_IR(ind_WN,:))
562	elseif(c_WN .eq. 'VI') then
563	call interpolate_T_abs_coeff(z_temp,temp_arr_CO2O2,num_T_CO2O2,abs_coef,abs_arr_CO2O2_VI(ind_WN,:))
564	else
565	print*,'You must select visible (VI) or infrared (IR) canal.'
566	stop
567	endif
568	elseif ((igas_X .eq. igas_H2) .and. (igas_Y .eq. igas_CH4)) then
569	call T_boundaries_continuum(z_temp,temp,temp_arr_H2CH4,num_T_H2CH4)
570	if(c_WN .eq. 'IR') then
571	call interpolate_T_abs_coeff(z_temp,temp_arr_H2CH4,num_T_H2CH4,abs_coef,abs_arr_H2CH4_IR(ind_WN,:))
572	elseif(c_WN .eq. 'VI') then
573	call interpolate_T_abs_coeff(z_temp,temp_arr_H2CH4,num_T_H2CH4,abs_coef,abs_arr_H2CH4_VI(ind_WN,:))
574	else
575	print*,'You must select visible (VI) or infrared (IR) canal.'
576	stop
577	endif
578	elseif ((igas_X .eq. igas_H2) .and. (igas_Y .eq. igas_He)) then
579	call T_boundaries_continuum(z_temp,temp,temp_arr_H2He,num_T_H2He)
580	if(c_WN .eq. 'IR') then
581	call interpolate_T_abs_coeff(z_temp,temp_arr_H2He,num_T_H2He,abs_coef,abs_arr_H2He_IR(ind_WN,:))
582	elseif(c_WN .eq. 'VI') then
583	call interpolate_T_abs_coeff(z_temp,temp_arr_H2He,num_T_H2He,abs_coef,abs_arr_H2He_VI(ind_WN,:))
584	else
585	print*,'You must select visible (VI) or infrared (IR) canal.'
586	stop
587	endif
588	elseif ((igas_X .eq. igas_H2O) .and. (igas_Y .eq. igas_N2)) then
589	call T_boundaries_continuum(z_temp,temp,temp_arr_H2ON2,num_T_H2ON2)
590	if(c_WN .eq. 'IR') then
591	call interpolate_T_abs_coeff(z_temp,temp_arr_H2ON2,num_T_H2ON2,abs_coef,abs_arr_H2ON2_IR(ind_WN,:))
592	elseif(c_WN .eq. 'VI') then
593	call interpolate_T_abs_coeff(z_temp,temp_arr_H2ON2,num_T_H2ON2,abs_coef,abs_arr_H2ON2_VI(ind_WN,:))
594	else
595	print*,'You must select visible (VI) or infrared (IR) canal.'
596	stop
597	endif
598	elseif ((igas_X .eq. igas_H2O) .and. (igas_Y .eq. igas_O2)) then
599	call T_boundaries_continuum(z_temp,temp,temp_arr_H2OO2,num_T_H2OO2)
600	if(c_WN .eq. 'IR') then
601	call interpolate_T_abs_coeff(z_temp,temp_arr_H2OO2,num_T_H2OO2,abs_coef,abs_arr_H2OO2_IR(ind_WN,:))
602	elseif(c_WN .eq. 'VI') then
603	call interpolate_T_abs_coeff(z_temp,temp_arr_H2OO2,num_T_H2OO2,abs_coef,abs_arr_H2OO2_VI(ind_WN,:))
604	else
605	print*,'You must select visible (VI) or infrared (IR) canal.'
606	stop
607	endif
608	elseif ((igas_X .eq. igas_H2O) .and. (igas_Y .eq. igas_CO2)) then
609	call T_boundaries_continuum(z_temp,temp,temp_arr_H2OCO2,num_T_H2OCO2)
610	if(c_WN .eq. 'IR') then
611	call interpolate_T_abs_coeff(z_temp,temp_arr_H2OCO2,num_T_H2OCO2,abs_coef,abs_arr_H2OCO2_IR(ind_WN,:))
612	elseif(c_WN .eq. 'VI') then
613	call interpolate_T_abs_coeff(z_temp,temp_arr_H2OCO2,num_T_H2OCO2,abs_coef,abs_arr_H2OCO2_VI(ind_WN,:))
614	else
615	print*,'You must select visible (VI) or infrared (IR) canal.'
616	stop
617	endif
618	elseif ((igas_X .eq. igas_CO2) .and. (igas_Y .eq. igas_H2)) then
619	call T_boundaries_continuum(z_temp,temp,temp_arr_CO2H2,num_T_CO2H2)
620	if(c_WN .eq. 'IR') then
621	call interpolate_T_abs_coeff(z_temp,temp_arr_CO2H2,num_T_CO2H2,abs_coef,abs_arr_CO2H2_IR(ind_WN,:))
622	elseif(c_WN .eq. 'VI') then
623	call interpolate_T_abs_coeff(z_temp,temp_arr_CO2H2,num_T_CO2H2,abs_coef,abs_arr_CO2H2_VI(ind_WN,:))
624	else
625	print*,'You must select visible (VI) or infrared (IR) canal.'
626	stop
627	endif
628	elseif ((igas_X .eq. igas_CO2) .and. (igas_Y .eq. igas_CH4)) then
629	call T_boundaries_continuum(z_temp,temp,temp_arr_CO2CH4,num_T_CO2CH4)
630	if(c_WN .eq. 'IR') then
631	call interpolate_T_abs_coeff(z_temp,temp_arr_CO2CH4,num_T_CO2CH4,abs_coef,abs_arr_CO2CH4_IR(ind_WN,:))
632	elseif(c_WN .eq. 'VI') then
633	call interpolate_T_abs_coeff(z_temp,temp_arr_CO2CH4,num_T_CO2CH4,abs_coef,abs_arr_CO2CH4_VI(ind_WN,:))
634	else
635	print*,'You must select visible (VI) or infrared (IR) canal.'
636	stop
637	endif
638	endif ! igas_X / igas_Y condition
639
640	! We compute the values of amagat for molecules X and Y
641	amagat_X = (273.15/temp)*(pres_X/101325.0)
642	amagat_Y = (273.15/temp)*(pres_Y/101325.0)
643
644	! We convert the absorption coefficient from cm^-1 amagat^-2 into m^-1
645	abs_coef=abs_coef100.0amagat_X*amagat_Y
646
647	!!!! TEST !!!
648	!abs_coef=abs_coef*30
649
650	!print*,'We have ',amagat_X,' amagats of molecule ', trim(gnom(igas_X))
651	!print*,'We have ',amagat_X,' amagats of molecule ', trim(gnom(igas_Y))
652	!print*,'So the absorption is ',abs_coef,' m^-1'
653
654	end subroutine interpolate_continuum
655
656
657	subroutine interpolate_wn_abs_coeff(wn_arr,num_wn,abs_arr_in,abs_arr_out_VI,abs_arr_out_IR,num_T)
658
659	!==================================================================
660	!
661	! Purpose
662	! -------
663	! Interpolate the continuum data into the visible (VI) and infrared (IR) spectral chanels.
664	!
665	! Author
666	! -------
667	! M. Turbet (2025)
668	!
669	!==================================================================
670
671	use radcommon_h, only : BWNV,BWNI,WNOI,WNOV
672	use radinc_h, only: L_NSPECTI, L_NSPECTV
673	use mod_phys_lmdz_para, only : is_master
674
675	implicit none
676
677	integer iW, iB, count_norm
678	integer,intent(in) :: num_T
679	integer,intent(in) :: num_wn
680	double precision,intent(in) :: wn_arr(num_wn)
681	double precision,intent(in) :: abs_arr_in(num_wn,num_T)
682	double precision,intent(out) :: abs_arr_out_IR(L_NSPECTI,num_T)
683	double precision,intent(out) :: abs_arr_out_VI(L_NSPECTV,num_T)
684
685	! First visible (VI) chanel
686
687	! We get read of all the wavenumbers lower than the minimum wavenumber in the visible wavenumber grid
688	iW=1
689	do while((wn_arr(iW) .lt. BWNV(1)) .and. (iW .lt. num_wn))
690	iW=iW+1
691	enddo
692
693	! We compute the mean of the continuum absorption inside each wavenumber visible (VI) chanel
694	do iB = 1, L_NSPECTV
695	count_norm=0
696	do while((wn_arr(iW) .lt. BWNV(iB+1)) .and. (iW .lt. num_wn))
697	abs_arr_out_VI(iB,:)=abs_arr_out_VI(iB,:)+abs_arr_in(iW,:)
698	count_norm=count_norm+1
699	iW=iW+1
700	enddo
701	if(count_norm .ge. 1) abs_arr_out_VI(iB,:)=abs_arr_out_VI(iB,:)/count_norm
702	end do
703
704	! Then infrared (IR) chanel
705
706	! We get read of all the wavenumbers lower than the minimum wavenumber in the infrared wavenumber grid
707	iW=1
708	do while((wn_arr(iW) .lt. BWNI(1)) .and. (iW .lt. num_wn))
709	iW=iW+1
710	enddo
711
712	! We compute the mean of the continuum absorption inside each wavenumber infrared (IR) chanel
713	do iB = 1, L_NSPECTI
714	count_norm=0
715	do while((wn_arr(iW) .lt. BWNI(iB+1)) .and. (iW .lt. num_wn))
716	abs_arr_out_IR(iB,:)=abs_arr_out_IR(iB,:)+abs_arr_in(iW,:)
717	count_norm=count_norm+1
718	iW=iW+1
719	enddo
720	if(count_norm .ge. 1) abs_arr_out_IR(iB,:)=abs_arr_out_IR(iB,:)/count_norm
721	end do
722
723	if (is_master) then
724	print*, 'Continuum absorption, first temperature, visible (VI):'
725	do iB = 1, L_NSPECTV
726	print*,WNOV(iB),' cm-1',abs_arr_out_VI(iB,1), ' cm-1 amagat-2'
727	end do
728
729	print*, 'Continuum absorption, first temperature, infrared (IR):'
730	do iB = 1, L_NSPECTI
731	print*,WNOI(iB),' cm-1',abs_arr_out_IR(iB,1), ' cm-1 amagat-2'
732	end do
733	endif
734
735	end subroutine interpolate_wn_abs_coeff
736
737
738	subroutine T_boundaries_continuum(z_temp,temp,temp_arr,num_T)
739
740	!==================================================================
741	!
742	! Purpose
743	! -------
744	! Check if the temperature is outside the boundaries of the continuum data temperatures.
745	!
746	! Author
747	! -------
748	! M. Turbet (2025)
749	!
750	!==================================================================
751
752	! use callkeys_mod, only: strictboundcia
753	use mod_phys_lmdz_para, only : is_master
754
755	implicit none
756
757	double precision,intent(out) :: z_temp
758	double precision,intent(in) :: temp
759	integer,intent(in) :: num_T
760	double precision,intent(in) :: temp_arr(num_T)
761
762	character(len=22) :: rname = "T_boundaries_continuum"
763
764	z_temp=temp
765
766	if(z_temp .lt. minval(temp_arr)) then
767	! if (strictboundcia) then
768	if (is_master) then
769	print*,'Your temperatures are too low for this continuum dataset'
770	print*, 'Minimum temperature is ', minval(temp_arr), ' K'
771	endif
772	call abort_physic(rname,"temperature too low",1)
773	! else
774	! z_temp=minval(temp_arr)
775	! endif
776	elseif(z_temp .gt. maxval(temp_arr)) then
777	! if (strictboundcia) then
778	if (is_master) then
779	print*,'Your temperatures are too high for this continuum dataset'
780	print*, 'Maximum temperature is ', maxval(temp_arr), ' K'
781	endif
782	call abort_physic(rname,"temperature too low",1)
783	! else
784	! z_temp=maxval(temp_arr)
785	! endif
786	endif
787
788	end subroutine T_boundaries_continuum
789
790
791	subroutine interpolate_T_abs_coeff(z_temp,temp_arr,num_T,abs_coef,abs_arr)
792
793	!==================================================================
794	!
795	! Purpose
796	! -------
797	! Interpolate in the continuum data using the temperature field
798	!
799	! Author
800	! -------
801	! M. Turbet (2025)
802	!
803	!==================================================================
804
805	implicit none
806
807	integer iT
808	double precision z_temp
809	integer num_T
810	double precision temp_arr(num_T)
811
812	double precision abs_coef
813	double precision abs_arr(num_T)
814
815	! Check where to interpolate
816	iT=1
817	do while ( z_temp .gt. temp_arr(iT) )
818	iT=iT+1
819	end do
820
821	! We proceed to a simple linear interpolation using the two most nearby temperatures
822	if(iT .lt. num_T) then
823	abs_coef=abs_arr(iT-1)+(abs_arr(iT)-abs_arr(iT-1))*(z_temp-temp_arr(iT-1))/(temp_arr(iT)-temp_arr(iT-1))
824	else
825	abs_coef=abs_arr(iT)
826	endif
827
828	!print*,'the absorption is ',abs_coef,' cm^-1 amagat^-2'
829
830
831	end subroutine interpolate_T_abs_coeff
832
833	end module interpolate_continuum_mod

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

Download in other formats:

Original Format