Context Navigation

interpolateH2H2.F90 @ 3058

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

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

Download in other formats:

Original Format