Context Navigation

nlte_ch4.F90

Last change on this file was 3275, checked in by afalco, 8 months ago
Pluto PCM: Changed _vap to _gas; Included surfprop.F90; callcorrk includes methane AF
File size: 10.3 KB

Line
1	subroutine nlte_ch4(ngrid,nlayer,nq,pplay,pplev,pt,vmrch4,eps23,eps33,eps76)
2
3	use comgeomfi_h
4	use callkeys_mod, only: strobel
5	use tracer_h, only: igcm_n2,mmol
6	use comcstfi_mod, only: pi, g, cpp,r
7	implicit none
8
9	!==================================================================
10	! Purpose
11	! -------
12	! Calculation of epsilon - non LTE efficiency factor
13	! For CH4 heating and cooling rate calculation
14	! Cf Strobel et al 1996
15	!
16	! Inputs
17	! ------
18	! ngrid Number of vertical columns
19	! nlayer Number of layers
20	! pplay(ngrid,nlayer) Pressure layers
21	! pt
22	!
23	! Outputs
24	! -------
25	!
26	! eps_sw ! efficiency factor for heating rate
27	! eps_lw ! efficiency factor for cooling rate
28	!
29	! Authors
30	! -------
31	! Tanguy Bertrand (2015,2022)
32	!
33	!==================================================================
34
35	!-----------------------------------------------------------------------
36
37	! Arguments
38
39	INTEGER ngrid, nlayer, nq
40	REAL pplay(ngrid,nlayer)
41	REAL pplev(ngrid,nlayer)
42	REAL pt(ngrid,nlayer)
43	REAL vmrch4(ngrid,nlayer)
44
45	!-----------------------------------------------------------------------
46	! Local variables
47
48	INTEGER l,ig
49
50	REAL cbol ! Boltzman constant
51	REAL avog ! avogadro constant
52	REAL beta23
53	REAL beta33
54	REAL nd(ngrid,nlayer) ! density number cm-3
55	REAL na(ngrid,nlayer) ! column density number cm-2
56	REAL rho(ngrid,nlayer) ! density of atmosphere
57	REAL kn3,k3,k2n4,k3n4,k4
58	REAL An3,A3,A2n4,A3n4,A4
59	REAL dphin3,dphi3,dphi2n4,dphi3n4,dphi4
60	REAL phin3(ngrid,nlayer)
61	REAL phi3(ngrid,nlayer)
62	REAL phi2n4(ngrid,nlayer)
63	REAL phi3n4(ngrid,nlayer)
64	REAL phi4(ngrid,nlayer)
65	REAL phi4T(ngrid,nlayer)
66
67	! DATA from Strobel et al 1996
68	REAL, DIMENSION(22) :: gam2 ! band escape fonction 2.3
69	REAL, DIMENSION(22) :: Naval2 ! column density number 2.3
70	REAL, DIMENSION(22) :: gam3 ! band escape fonction 3.3
71	REAL, DIMENSION(22) :: Naval3 ! column density number 3.3
72	REAL gammafi23(nlayer),gammafi33(nlayer) ! Band espace fonction
73	REAL gammafi23mod(nlayer),gammafi33mod(nlayer) ! Band espace fonction
74	REAL vmr
75	REAL k4zT ! function to calculate k4 depending on T and vmr
76	REAL log_square ! log function
77	REAL log_affine ! log function
78	REAL ratio_affine ! linear approx function to calculate ratio for gamma calculation
79	REAL ratio_square !
80	REAL ratio_sq,ratio_af
81
82	! Output
83	REAL,intent(out) :: eps23(ngrid,nlayer)
84	REAL,intent(out) :: eps33(ngrid,nlayer)
85	REAL,intent(out) :: eps76(ngrid,nlayer)
86
87	! Option simplified coefficient
88	REAL klw, ksw, eps_sw1Pa
89
90	!-----------------------------------------------------------------------
91	eps23(:,:)=1.
92	eps33(:,:)=1.
93	eps76(:,:)=1.
94
95	IF (strobel) then
96
97	! Constantes
98	cbol=1.3806488e-23 !Boltzman constant
99	avog=6.022141e23
100
101	! band 3
102	k3=2.8e-11 !cm-3
103	A3=25.2
104	dphi3=k3/A3
105	! band n3
106	kn3=2.8e-11 !cm-3
107	An3=27.3
108	dphin3=kn3/An3
109	! band 3n4
110	k3n4=1.16e-11 !cm-3
111	A3n4=2.12
112	dphi3n4=k3n4/(3.*A3n4)
113	! band 4
114	!! Here we use function k4zT instead of constant k4
115	k4=4.5e-16 !cm-3
116	A4=2.12
117	dphi4=k4/A4
118	! band 2n4
119	k2n4=1.16e-11
120	A2n4=2.12
121	dphi2n4=k2n4/(2.*A2n4)
122
123	! 2.3 band :
124	beta23=0.09
125	! 3.3 band :
126	beta33=0.13
127
128	! Escape Function Gamma (Data Strobel et al) :
129	! 100 K 2.3
130	Naval2=(/1.00000000e+10,1.30000000e+15,2.92766055e+15,6.59322793e+15,1.48482564e+16,3.34389650e+16, &
131	7.53061067e+16,1.69592860e+17,3.81931020e+17,8.60126447e+17,1.93704482e+18,4.36231516e+18, &
132	9.82413694e+18,2.21244140e+19,4.98252108e+19,1.12208695e+20,2.52699209e+20,5.69090388e+20, &
133	1.28161806e+21,2.88626357e+21,6.50000000e+21,1.00000000e+24/)
134	gam2=(/1.00000000e+00,1.00000000e+00,1.00000000e+00,1.00000000e+00,1.00000000e+00,9.90000000e-01, &
135	9.80000000e-01,9.50000000e-0,9.00000000e-01,8.00000000e-01,7.00000000e-01,5.00000000e-01, &
136	3.50000000e-01,2.20000000e-01,1.20000000e-01,6.50000000e-02,3.20000000e-02,1.40000000e-02, &
137	6.20000000e-03,2.50000000e-03,1.05000000e-03,4.83227951e-06/)
138
139	! 100 K 3.3
140	Naval3=(/1.00000000e+10,1.30000000e+15,2.92766055e+15,6.59322793e+15,1.48482564e+16,3.34389650e+16, &
141	7.53061067e+16,1.69592860e+17,3.81931020e+17,8.60126447e+17,1.93704482e+18,4.36231516e+18, &
142	9.82413694e+18,2.21244140e+19,4.98252108e+19,1.12208695e+20,2.52699209e+20,5.69090388e+20, &
143	1.28161806e+21,2.88626357e+21,6.50000000e+21,1.00000000e+24/)
144	gam3=(/1.00000000e+00,9.80000000e-01,9.00000000e-01,8.50000000e-01,7.00000000e-01,5.00000000e-01, &
145	3.00000000e-01,1.50000000e-01,1.00000000e-01,6.00000000e-02,4.00000000e-02,2.30000000e-02, &
146	1.45000000e-02,9.00000000e-03,5.00000000e-03,3.00000000e-03,1.50000000e-03,8.00000000e-04, &
147	3.80000000e-04,1.80000000e-04,8.30000000e-05,6.81737544e-07/)
148
149	DO ig=1,ngrid
150
151	phi3n4(ig,:)=1.
152	phi2n4(ig,:)=1.
153	phin3(ig,:)=1.
154	phi4(ig,:)=1.
155	nd(ig,:)=1.
156	na(ig,:)=1.
157
158	!write(,) 'TB22 pplev=',pplev(1,:)
159	!write(,) 'TB22 tlev=',pt(1,:)
160
161	DO l=1,nlayer
162	!!calculation of density
163	rho(ig,l)=pplay(ig,l)/(r*pt(ig,l))
164
165	!!calculation of column number density for ch4 (cm-2)
166	na(ig,l)=vmrch4(ig,l)/100./(mmol(igcm_n2)1.e-3) &
167	pplev(ig,l)/gavog1.e-4
168
169	!!calculation of atmospheric number density (cm-3) for each layer
170	nd(ig,l)=rho(ig,l)/(mmol(igcm_n2)1.e-3)avog*1.e-6
171
172	!!calculation of the phi :
173	phin3(ig,l)=dphin3*nd(ig,l)
174	phi3(ig,l)=dphi3*nd(ig,l)
175	phi2n4(ig,l)=dphi2n4*nd(ig,l)
176	phi3n4(ig,l)=dphi3n4*nd(ig,l)
177	phi4(ig,l)=dphi4*nd(ig,l)
178	!! Getting phi4 depending on T and vmr
179	vmr=vmrch4(ig,l)/100.
180	phi4T(ig,l)=k4zT(pt(ig,l),vmr)/A4*nd(ig,l)
181	ENDDO
182	!write(,) 'TB22 nd=',nd(1,:)
183	!write(,) 'TB22 vmr=',vmrch4
184	!write(,) 'TB22 phi4T=',phi4T
185
186	!! interpolation of na values in cm-2 from Strobel values for
187	!the 100K model:
188	CALL interp_line(Naval2,gam2,22,na(ig,:),gammafi23,nlayer)
189	CALL interp_line(Naval3,gam3,22,na(ig,:),gammafi33,nlayer)
190
191	!write(,) 'TB22rad na = ',na(1,:)
192	!write(,) 'TB22rad gammafi = ',gammafi23
193
194	DO l=1,nlayer
195
196	!! Getting ratio depending on temperature and n (cm-2)
197	ratio_sq=ratio_square(pt(ig,l),na(ig,l))
198	ratio_af=ratio_affine(pt(ig,l),na(ig,l))
199	!! Modified gamma depending on T
200	gammafi23mod(l)=min(ratio_sq*gammafi23(l),1.)
201	gammafi33mod(l)=min(ratio_af*gammafi33(l),1.)
202
203	!! Heating : cf Strobel 1996 and Zalucha 2011
204	! eps 2.3
205	eps23(ig,l)=phin3(ig,l)/(phin3(ig,l)+gammafi23mod(l))* &
206	(beta23+phi3n4(ig,l)/(phi3n4(ig,l)+1)* &
207	phi4T(ig,l)/(phi4T(ig,l)+gammafi23mod(l))* &
208	(3.phi2n4(ig,l)+1.)/(1.+phi2n4(ig,l))(1-beta23)/3.)
209
210	! eps 3.3
211	eps33(ig,l)=phi3(ig,l)/(phi3(ig,l)+gammafi33mod(l))* &
212	(beta33+phi2n4(ig,l)/(phi2n4(ig,l)+1)* &
213	phi4T(ig,l)/(phi4T(ig,l)+gammafi33mod(l))* &
214	(1-beta33))
215
216	! Cooling from Zalucha et al 2013 eq10
217	eps76(ig,l)=1./(1.+gammafi23mod(l)/(2.*phi4T(ig,l)))
218
219	ENDDO
220	!write(,) 'TB22rad eps23 = ',eps23(1,:)
221	!write(,) 'TB22rad eps33 = ',eps33(1,:)
222	!write(,) 'TB22rad gamma33 = ',gammafi33mod
223	!write(,) 'TB22rad eps76 = ',eps76(1,:)
224
225	ENDDO ! ngrid
226
227	ELSE ! if (strobel)
228	!if strobel = false
229
230	klw = 3.6E-4 ! fit strobel
231	!eps_sw1Pa = 0.6 ! fit strobel
232	eps_sw1Pa = 0.90 ! to fit NH temperature en 1D [ch4]=0.5
233	!eps_sw1Pa = 0.98 ! to fit NH temperature en 3D = fit 1D [ch4]=0.25
234
235	ksw = 8.5E-5*(0.9/(eps_sw1Pa-0.1)-1.)
236	do l=1, nlayer
237	do ig=1, ngrid
238	rho(ig,l)=pplay(ig,l)/(r*pt(ig,l))
239	eps23(ig,l) = 0.1+ 0.9/(1.+ksw/rho(ig,l))
240	eps76(ig,l) = 1./(1.+klw/rho(ig,l))
241	end do
242	end do
243
244	END IF
245
246	end subroutine nlte_ch4
247
248	real function k4zT(temp,vmr)
249	implicit none
250	real, intent(in) :: temp,vmr
251	k4zT=1.9e-15exp((temp-240.)/70.)(1+10.*vmr)
252	return
253	end function k4zT
254
255	!! Log log affine and square function for ratio adjustment and extrapolation
256
257	real function log_square(x,y1,y2,x1,x2,coef)
258	implicit none
259	real, intent(in) :: x,y1,y2,x1,x2,coef
260	real aa,bb
261	aa=(log(y2)-log(y1))/((log(coefx2))4-(log(coefx1)**4))
262	bb=log(y2)-aa(log(coefx2))**4
263	log_square=exp(aa(log(coefx))**4+bb)
264	return
265	end function log_square
266
267	real function log_affine(x,y1,y2,x1,x2)
268	real, intent(in) :: x,y1,y2,x1,x2
269	real aa,bb
270	aa=(log(y2)-log(y1))/(log(x2)-log(x1))
271	bb=log(y2)-aa*log(x2)
272	log_affine=exp(aa*log(x)+bb)
273	return
274	end function log_affine
275
276	!! Ratio of gamma depending on T. Linear approximation
277	real function ratio_affine(temp,n)
278	implicit none
279	real, intent(in) :: temp,n
280	real aa,bb,ratio0
281	real log_affine
282	aa=(1.-0.2)/(100.-40.)
283	bb=1.-aa*100.
284	ratio0=aa*temp+bb
285	ratio_affine=log_affine(n,1.,ratio0,1e15,1e22)
286	return
287	end function ratio_affine
288
289	!! Ratio of gamma depending on T. Linear approximation
290	real function ratio_square(temp,n)
291	implicit none
292	real, intent(in) :: temp,n
293	real aa,bb,ratio0,coef
294	real log_square
295	coef=1.e-12
296	aa=(1.-0.2)/(100.-40.)
297	bb=1.-aa*100.
298	ratio0=aa*temp+bb
299	ratio_square=log_square(n,1.,ratio0,1e15,1e22,coef)
300	return
301	end function ratio_square
302
303
304

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

Download in other formats:

Original Format