Context Navigation

interpolateH2H2.F90 @ 2860

Last change on this file since 2860 was 2837, checked in by aslmd, 3 years ago

Changed the H2H2 CIA file in the Hot Jupiter case

Modified the default H2H2 CIA file. The file used now
(H2-H2_2011_extended.cia) is a combinaison of multiple files existing
for hot temperature. It is very similar to the H2H2 CIA used in Baudino
et al, 2015, but has been extended to 7000 K instead of 5000 K.

Property svn:executable set to *

File size: 9.4 KB

Line
1
2	subroutine interpolateH2H2(wn,temp,pres,abcoef,firstcall,ind)
3
4	!==================================================================
5	!
6	! Purpose
7	! -------
8	! Calculates the H2-H2 CIA opacity, using a lookup table from
9	! HITRAN (2011 or later)
10	!
11	! Authors
12	! -------
13	! R. Wordsworth (2011)
14	!
15	! + J.Vatant d'Ollone (2019)
16	! - Enable updated HITRAN file (Karman2019,Fletcher2018)
17	! (2018 one should be default for giant planets)
18	!==================================================================
19
20	use callkeys_mod, only: versH2H2cia, strictboundcia, H2orthopara_mixture
21	use datafile_mod, only: datadir
22	use mod_phys_lmdz_para, only : is_master
23
24	implicit none
25
26	! input
27	double precision wn ! wavenumber (cm^-1)
28	double precision temp ! temperature (Kelvin)
29	double precision pres ! pressure (Pascals)
30
31	! output
32	double precision abcoef ! absorption coefficient (m^-1)
33
34	integer nS,nT
35	parameter(nT=10)
36
37	double precision, parameter :: losch = 2.6867774e19
38	! Loschmit's number (molecule cm^-3 at STP)
39	! converts cm^5 molecule^-2 --> cm^-1 amagat^-2
40	! see Richard et al. 2011, JQSRT for details
41
42	double precision amagat
43	double precision temp_arr(nT)
44
45	double precision, dimension(:), allocatable :: wn_arr
46	double precision, dimension(:,:), allocatable :: abs_arr
47
48	integer k,iT
49	logical firstcall
50
51	save nS, wn_arr, temp_arr, abs_arr !read by master
52
53	character*100 dt_file
54	integer ios
55
56	character(LEN=*), parameter :: fmat11 = "(A20,F10.3,F10.3,I7,F7.1,E10.3,F5.3)"
57	character(LEN=*), parameter :: fmat18 = "(A12,A3,A5,F10.6,F10.4,I7,F7.3,E10.3,F5.3)"
58
59	character*20 bleh
60	double precision blah, Ttemp
61	integer nres
62
63	integer ind
64	if ((H2orthopara_mixture .eq. "hot")) then
65	if (temp .gt. 7000.) then
66	if (strictboundcia) then
67	if (is_master) then
68	print*,'Your temperatures are too high for this H2-H2 CIA dataset (>7000K, Hot Jupiter case)'
69	endif !is_master
70	stop
71	else
72	if (is_master) then
73	print*,'Your temperatures are too high for this H2-H2 CIA dataset (Hot Jupiter case)'
74	print*,'you have chosen strictboundcia = ', strictboundcia
75	print,'********************************************************'
76	print*,' we allow model to continue but with temp = 7000 '
77	print*,' ... for H2-H2 CIA dataset ... '
78	print*,' ... we assume we know what you are doing ... '
79	print,'********************************************************'
80	endif !is_master
81	temp = 7000.
82	endif !strictboundcia
83	endif !(temp .gt. 7000.)
84	else ! if not "hot"
85	if(temp.gt.400)then
86	if (strictboundcia) then
87	if (is_master) then
88	print*,'Your temperatures are too high for this H2-H2 CIA dataset. If you '
89	print*,'really want to run simulations with hydrogen at T > 400 K, contact'
90	print*,'Robin Wordsworth [rwordsworth@uchicago.edu].'
91	endif !is_master
92	stop
93	else
94	if (is_master) then
95	print*,'Your temperatures are too high for this H2-H2 CIA dataset'
96	print*,'you have chosen strictboundcia = ', strictboundcia
97	print,'********************************************************'
98	print*,' we allow model to continue but with temp = 400 '
99	print*,' ... for H2-H2 CIA dataset ... '
100	print*,' ... we assume we know what you are doing ... '
101	print,'********************************************************'
102	endif !is_master
103	temp = 400
104	endif !of strictboundcia
105	elseif(temp.lt.40)then
106	if (strictboundcia) then
107	if (is_master) then
108	print*,'Your temperatures are too low for this H2-H2 CIA dataset. If you '
109	print*,'really want to run simulations with hydrogen at T < 40 K, contact'
110	print*,'Robin Wordsworth [rwordsworth@uchicago.edu].'
111	endif ! is_master
112	stop
113	else
114	if (is_master) then
115	print*,'Your temperatures are too low for this H2-H2 CIA dataset'
116	print*,'you have chosen strictboundcia = ', strictboundcia
117	print,'********************************************************'
118	print*,' we allow model to continue but with temp = 40 '
119	print*,' ... for H2-H2 CIA dataset ... '
120	print*,' ... we assume we know what you are doing ... '
121	print,'********************************************************'
122	endif !is_master
123	temp = 40
124	endif !of strictboundcia
125	endif ! of (temp .gt. 400)
126	endif ! of ((H2orthopara_mixture .eq. "hot").and. (temp .gt. 3000.))
127	amagat = (273.15/temp)*(pres/101325.0)
128
129	if(firstcall)then ! called by sugas_corrk only
130	if (is_master) print*,'----------------------------------------------------'
131	if (is_master) print*,'Initialising H2-H2 continuum from HITRAN database...'
132
133	! 1.1 Open the ASCII files and set nS according to version
134	! Only two possible versions for now : 2011 or 2018 (sanity check in inifis_mod)
135	if (versH2H2cia.eq.2011) then
136	nS = 2428
137	if (H2orthopara_mixture.eq."normal") then
138	dt_file=TRIM(datadir)//'/continuum_data/H2-H2_norm_2011.cia'
139	else if (H2orthopara_mixture.eq."equilibrium") then
140	dt_file=TRIM(datadir)//'/continuum_data/H2-H2_eq_2011.cia'
141	else if (H2orthopara_mixture.eq."hot") then
142	dt_file=TRIM(datadir)//'/continuum_data/H2-H2_2011_extended.cia'
143	ns = 800
144	endif
145	else if (versH2H2cia.eq.2018) then
146	nS = 9600
147	if (H2orthopara_mixture.eq."normal") then
148	dt_file=TRIM(datadir)//'/continuum_data/H2-H2_norm_2018.cia'
149	else if (H2orthopara_mixture.eq."equilibrium") then
150	dt_file=TRIM(datadir)//'/continuum_data/H2-H2_eq_2018.cia'
151	endif
152	endif
153
154	if(.not.allocated(wn_arr)) allocate(wn_arr(nS))
155	if(.not.allocated(abs_arr)) allocate(abs_arr(nS,nT))
156
157	!$OMP MASTER
158	open(33,file=dt_file,form='formatted',status='old',iostat=ios)
159	if (ios.ne.0) then ! file not found
160	if (is_master) then
161	write(,) 'Error from interpolateH2H2'
162	write(,) 'Data file ',trim(dt_file),' not found.'
163	write(,) 'Check that your path to datagcm:',trim(datadir)
164	write(,) 'is correct. You can change it in callphys.def with:'
165	write(,) 'datadir = /absolute/path/to/datagcm'
166	write(,) 'Also check that the continuum data is there.'
167	endif
168	call abort
169	else
170
171	if(versH2H2cia.eq.2011) then
172	if (is_master) then
173	write(,) '... You are using H2-H2 CIA from 2011 but you should use more recent data available on HITRAN !'
174	write(,) '... (Especially if you are running a giant planet atmosphere)'
175	write(,) '... Just find out the H2-H2 CIA from 2018, put it in your datadir and have a look at interpolateH2H2.F90 ! .'
176	endif
177	endif
178
179	do iT=1,nT
180
181	! Only two possibles values for now : 2011 or 2018 (sanity check in inifis_mod)
182	if(versH2H2cia.eq.2011) then
183	read(33,fmat11) bleh,blah,blah,nres,Ttemp
184	else if (versH2H2cia.eq.2018) then
185	read(33,fmat18) bleh,bleh,bleh,blah,blah,nres,Ttemp
186	endif
187
188	if(nS.ne.nres)then
189	if (is_master) then
190	print*,'Resolution given in file: ',trim(dt_file)
191	print*,'is ',nres,', which does not match nS.'
192	print*,'Please adjust nS value in interpolateH2H2.F90'
193	endif
194	stop
195	endif
196	temp_arr(iT)=Ttemp
197
198	do k=1,nS
199	read(33,*) wn_arr(k),abs_arr(k,it)
200	end do
201
202	end do
203
204	endif
205	close(33)
206	!$OMP END MASTER
207	!$OMP BARRIER
208
209	if (is_master) then
210	print*,'interpolateH2H2: At wavenumber ',wn,' cm^-1'
211	print*,' temperature ',temp,' K'
212	print*,' pressure ',pres,' Pa'
213	endif
214	endif
215
216
217
218	call bilinearbig(nS,nT,wn_arr,temp_arr,abs_arr,wn,temp,abcoef,ind)
219
220	!print*,'the absorption is ',abcoef,' cm^5 molecule^-2'
221	!print,'or ',abcoeflosch**2,' cm^-1 amagat^-2'
222
223	abcoef=abcoeflosch2100.0amagat*2 ! convert to m^-1
224
225	!print*,'We have ',amagat,' amagats of H2'
226	!print*,'So the absorption is ',abcoef,' m^-1'
227
228	! unlike for Rayleigh scattering, we do not currently weight by the BB function
229	! however our bands are normally thin, so this is no big deal.
230
231	return
232	end subroutine interpolateH2H2

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

Download in other formats:

Original Format