Context Navigation

interpolateH2H2.F90 @ 2733

Last change on this file since 2733 was 2667, checked in by aslmd, 3 years ago
CIA : modified the loading of file to use high temperature dataset
Property svn:executable set to ``*
File size: 9.6 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	! print*,"We're in the Hot Jupiter case"
66	if (temp .gt. 3000.) then
67	if (strictboundcia) then
68	if (is_master) then
69	print*,'Your temperatures are too high for this H2-H2 CIA dataset (Hot Jupiter case). If you '
70	print*,'really want to run simulations with hydrogen at T > 400 K, contact'
71	print*,'Robin Wordsworth [rwordsworth@uchicago.edu].'
72	endif !is_master
73	stop
74	else
75	if (is_master) then
76	print*,'Your temperatures are too high for this H2-H2 CIA dataset (Hot Jupiter case)'
77	print*,'you have chosen strictboundcia = ', strictboundcia
78	print,'********************************************************'
79	print*,' we allow model to continue but with temp = 3000 '
80	print*,' ... for H2-H2 CIA dataset ... '
81	print*,' ... we assume we know what you are doing ... '
82	print,'********************************************************'
83	endif !is_master
84	temp = 3000.
85	endif !strictboundcia
86	endif !(temp .gt. 3000.)
87	else ! if not Hot Jupiter
88	if(temp.gt.400)then
89	if (strictboundcia) then
90	if (is_master) then
91	print*,'Your temperatures are too high for this H2-H2 CIA dataset. If you '
92	print*,'really want to run simulations with hydrogen at T > 400 K, contact'
93	print*,'Robin Wordsworth [rwordsworth@uchicago.edu].'
94	endif !is_master
95	stop
96	else
97	if (is_master) then
98	print*,'Your temperatures are too high for this H2-H2 CIA dataset'
99	print*,'you have chosen strictboundcia = ', strictboundcia
100	print,'********************************************************'
101	print*,' we allow model to continue but with temp = 400 '
102	print*,' ... for H2-H2 CIA dataset ... '
103	print*,' ... we assume we know what you are doing ... '
104	print,'********************************************************'
105	endif !is_master
106	temp = 400
107	endif !of strictboundcia
108	elseif(temp.lt.40)then
109	if (strictboundcia) then
110	if (is_master) then
111	print*,'Your temperatures are too low for this H2-H2 CIA dataset. If you '
112	print*,'really want to run simulations with hydrogen at T < 40 K, contact'
113	print*,'Robin Wordsworth [rwordsworth@uchicago.edu].'
114	endif ! is_master
115	stop
116	else
117	if (is_master) then
118	print*,'Your temperatures are too low for this H2-H2 CIA dataset'
119	print*,'you have chosen strictboundcia = ', strictboundcia
120	print,'********************************************************'
121	print*,' we allow model to continue but with temp = 40 '
122	print*,' ... for H2-H2 CIA dataset ... '
123	print*,' ... we assume we know what you are doing ... '
124	print,'********************************************************'
125	endif !is_master
126	temp = 40
127	endif !of strictboundcia
128	endif ! of (temp .gt. 400)
129	endif ! of ((H2orthopara_mixture .eq. "hot").and. (temp .gt. 3000.))
130	amagat = (273.15/temp)*(pres/101325.0)
131
132	if(firstcall)then ! called by sugas_corrk only
133	if (is_master) print*,'----------------------------------------------------'
134	if (is_master) print*,'Initialising H2-H2 continuum from HITRAN database...'
135
136	! 1.1 Open the ASCII files and set nS according to version
137	! Only two possible versions for now : 2011 or 2018 (sanity check in inifis_mod)
138	if (versH2H2cia.eq.2011) then
139	nS = 2428
140	if (H2orthopara_mixture.eq."normal") then
141	dt_file=TRIM(datadir)//'/continuum_data/H2-H2_norm_2011.cia'
142	else if (H2orthopara_mixture.eq."equilibrium") then
143	dt_file=TRIM(datadir)//'/continuum_data/H2-H2_eq_2011.cia'
144	else if (H2orthopara_mixture.eq."hot") then
145	dt_file=TRIM(datadir)//'/continuum_data/H2-H2_2011.cia'
146	ns = 9981
147	endif
148	else if (versH2H2cia.eq.2018) then
149	nS = 9600
150	if (H2orthopara_mixture.eq."normal") then
151	dt_file=TRIM(datadir)//'/continuum_data/H2-H2_norm_2018.cia'
152	else if (H2orthopara_mixture.eq."equilibrium") then
153	dt_file=TRIM(datadir)//'/continuum_data/H2-H2_eq_2018.cia'
154	endif
155	endif
156
157	if(.not.allocated(wn_arr)) allocate(wn_arr(nS))
158	if(.not.allocated(abs_arr)) allocate(abs_arr(nS,nT))
159
160	!$OMP MASTER
161	open(33,file=dt_file,form='formatted',status='old',iostat=ios)
162	if (ios.ne.0) then ! file not found
163	if (is_master) then
164	write(,) 'Error from interpolateH2H2'
165	write(,) 'Data file ',trim(dt_file),' not found.'
166	write(,) 'Check that your path to datagcm:',trim(datadir)
167	write(,) 'is correct. You can change it in callphys.def with:'
168	write(,) 'datadir = /absolute/path/to/datagcm'
169	write(,) 'Also check that the continuum data is there.'
170	endif
171	call abort
172	else
173
174	if(versH2H2cia.eq.2011) then
175	if (is_master) then
176	write(,) '... You are using H2-H2 CIA from 2011 but you should use more recent data available on HITRAN !'
177	write(,) '... (Especially if you are running a giant planet atmosphere)'
178	write(,) '... Just find out the H2-H2 CIA from 2018, put it in your datadir and have a look at interpolateH2H2.F90 ! .'
179	endif
180	endif
181
182	do iT=1,nT
183
184	! Only two possibles values for now : 2011 or 2018 (sanity check in inifis_mod)
185	if(versH2H2cia.eq.2011) then
186	read(33,fmat11) bleh,blah,blah,nres,Ttemp
187	else if (versH2H2cia.eq.2018) then
188	read(33,fmat18) bleh,bleh,bleh,blah,blah,nres,Ttemp
189	endif
190
191	if(nS.ne.nres)then
192	if (is_master) then
193	print*,'Resolution given in file: ',trim(dt_file)
194	print*,'is ',nres,', which does not match nS.'
195	print*,'Please adjust nS value in interpolateH2H2.F90'
196	endif
197	stop
198	endif
199	temp_arr(iT)=Ttemp
200
201	do k=1,nS
202	read(33,*) wn_arr(k),abs_arr(k,it)
203	end do
204
205	end do
206
207	endif
208	close(33)
209	!$OMP END MASTER
210	!$OMP BARRIER
211
212	if (is_master) then
213	print*,'interpolateH2H2: At wavenumber ',wn,' cm^-1'
214	print*,' temperature ',temp,' K'
215	print*,' pressure ',pres,' Pa'
216	endif
217	endif
218
219
220
221	call bilinearbig(nS,nT,wn_arr,temp_arr,abs_arr,wn,temp,abcoef,ind)
222
223	!print*,'the absorption is ',abcoef,' cm^5 molecule^-2'
224	!print,'or ',abcoeflosch**2,' cm^-1 amagat^-2'
225
226	abcoef=abcoeflosch2100.0amagat*2 ! convert to m^-1
227
228	!print*,'We have ',amagat,' amagats of H2'
229	!print*,'So the absorption is ',abcoef,' m^-1'
230
231	! unlike for Rayleigh scattering, we do not currently weight by the BB function
232	! however our bands are normally thin, so this is no big deal.
233
234	return
235	end subroutine interpolateH2H2

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

Download in other formats:

Original Format