Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

nlte_ch4.F90_old_11_2022 @ 4103

Last change on this file since 4103 was 3175, checked in by emillour, 2 years ago
Pluto PCM: Add the old Pluto LMDZ for reference (required prior step to making an LMDZ.PLUTO using the same framework as the other physics packages). TB+EM
File size: 7.8 KB

Line
1	subroutine nlte_ch4(ngrid,nlayer,nq,pplay,pplev,pt,pq,eps23,eps76)
2
3	use comgeomfi_h
4	implicit none
5
6	!==================================================================
7	! Purpose
8	! -------
9	! Calculation of epsilon - non LTE efficiency factor
10	! For CH4 heating and cooling rate calculation
11	! Cf Strobel et al 1996
12	!
13	! Inputs
14	! ------
15	! ngrid Number of vertical columns
16	! nlayer Number of layers
17	! pplay(ngrid,nlayer) Pressure layers
18	! pt
19	!
20	! Outputs
21	! -------
22	!
23	! eps_sw ! efficiency factor for heating rate
24	! eps_lw ! efficiency factor for cooling rate
25	!
26	! Authors
27	! -------
28	! Tanguy Bertrand (2015)
29	!
30	!==================================================================
31
32	#include "dimensions.h"
33	#include "dimphys.h"
34	#include "comcstfi.h"
35	#include "surfdat.h"
36	#include "comvert.h"
37	#include "callkeys.h"
38	#include "tracer.h"
39
40	!-----------------------------------------------------------------------
41
42	! Arguments
43
44	INTEGER ngrid, nlayer, nq
45	REAL pplay(ngrid,nlayer)
46	REAL pplev(ngrid,nlayer)
47	REAL pt(ngrid,nlayer)
48	REAL pq(ngrid,nlayer,nq)
49
50	!-----------------------------------------------------------------------
51	! Local variables
52
53	INTEGER l,ig
54
55	REAL temper ! atmospheric temperature
56	REAL cbol ! Boltzman constant
57	REAL avog ! avogadro constant
58	REAL beta23
59	REAL beta33
60	REAL nd(ngrid,nlayer) ! density number cm-3
61	REAL na(ngrid,nlayer) ! column density number cm-2
62	REAL rho(ngrid,nlayer) ! density of atmosphere
63	REAL kn3,k3,k2n4,k3n4,k4
64	REAL An3,A3,A2n4,A3n4,A4
65	REAL dphin3,dphi3,dphi2n4,dphi3n4,dphi4
66	REAL phin3(ngrid,nlayer)
67	REAL phi3(ngrid,nlayer)
68	REAL phi2n4(ngrid,nlayer)
69	REAL phi3n4(ngrid,nlayer)
70	REAL phi4(ngrid,nlayer)
71
72	! DATA from Strobel et al 1996
73	REAL, DIMENSION(12) :: gam2 ! band escape fonction 2.3
74	REAL, DIMENSION(12) :: Naval2 ! column density number 2.3
75	REAL, DIMENSION(14) :: gam3 ! band escape fonction 3.3
76	REAL, DIMENSION(14) :: Naval3 ! column density number 3.3
77	REAL, DIMENSION(20) :: Navalcool
78	REAL, DIMENSION(20) :: gamcool
79	REAL gammafi23(nlayer),gammafi33(nlayer) ! Band espace fonction
80
81	! Output
82
83	REAL eps23(ngrid,nlayer)
84	REAL eps33(ngrid,nlayer)
85	REAL eps76(ngrid,nlayer)
86
87	! Option simplified coefficient
88	logical, save :: strobel=.true. !if false, coef are calculated from linear fit
89	REAL klw, ksw, eps_sw1Pa
90
91	!-----------------------------------------------------------------------
92
93	IF (strobel) then
94
95	! Constantes
96	cbol=1.3806488e-23 !Boltzman constant
97	avog = 6.022141e23
98	temper=100. ! cte Temperature atmosphere K
99
100	! Calculation phi
101	! A2=25.8 !Einstein coefficient sec-1
102	! P3=7.e-3 ! collisional probability
103
104	! band 3
105	k3=2.8e-11 !cm-3
106	A3=25.2
107	dphi3=k3/A3
108	! band n3
109	kn3=2.8e-11 !cm-3
110	An3=27.3
111	dphin3=kn3/An3
112	! band 3n4
113	k3n4=1.16e-11 !cm-3
114	A3n4=2.12
115	dphi3n4=k3n4/(3*A3n4)
116	! band 4
117	k4=4.5e-16 !cm-3
118	A4=2.12
119	dphi4=k4/A4
120	! band 2n4
121	k2n4=1.16e-11
122	A2n4=2.12
123	dphi2n4=k2n4/(2*A2n4)
124
125	! 2.3 band :
126	beta23=0.09
127	! 3.3 band :
128	beta23=0.13
129
130	! Escape Function Gamma (Data Strobel et al) :
131	! 100 K 2.3
132	Naval2=(/1.e10,1.e16,8.e16,4.e17,2.e18,1.e19,5.e19,2.4e20,1.3e21,6.e21,1.e23,1.e24/)
133	gam2=(/1.,1.,0.95,0.9,0.65,0.38,0.14,4.e-2,6.e-3,1.1e-3,6.e-5,5.e-6/)
134	! 40 K 2.3
135	! gam2=(/1.,1.,0.95,0.7,0.47,0.16,0.05,7.e-3,1.2e-3,2.e-4,1.e-5,7.e-7/)
136	! 100 K 3.3
137	Naval3=(/1e10,1e15,1.4e15,3e15,1.5e16,7e16,4e17,2e18,1e19,5e19,2.7e20,1.3e21,7e21,1e24/)
138	gam3=(/1.,1.,0.95,0.9,0.7,0.3,0.1,0.04,0.015,5.e-3,1.5e-3,4.e-4,8.e-5,1.e-6/)
139
140	! gam=(/1.,1.,1.,1.,0.95,0.9,0.8,0.7,0.6,0.5,0.2,0.12,0.07,0.03,0.014,0.006,0.0026,0.0007,1.e-4,1.e-5/) ! values of gamma
141	! Naval=(/1e10,1.2e15,1e16,2e16,2e17,4e17,9e17,2e18,4e18,1e19,2e19,5e19,1.3e20,2.4e20,6e20,1.4e20,3e21,1e22,1e23,1e24/) ! values of Na in cm-2
142	! gamcool=(/8.79e15,8.79e15,8.79e15,1.87e16,1.94e17,3.79e17,8.04e17,1.62e18,2.92e18,6.22e18,9.72e18,1.45e19,2.21e19, 2.76e19,3.55e19,4.01e19,5.24e19,6.39e19,9.99e19,1.49e20/)
143	! Navalcool=(/1e10,1.12e16,1.12e16, 3.12e16, 2.31e17, 6.31e17,1.53e18, 3.53e18, 7.53e18, 1.75e19, 3.75e19,8.75e19, 2.175e20, 4.575e20, 1.057e21, 1.197e21,4.19e21, 1.419e22, 1.14e23,1.114e24/)
144	! gamcool=gam
145
146	DO ig=1,ngrid
147
148	eps23(ig,:)=0.
149	eps33(ig,:)=0.
150	eps76(ig,:)=0.
151	phi3n4(ig,:)=0.
152	phi2n4(ig,:)=0.
153	phin3(ig,:)=0.
154	phi4(ig,:)=0.
155	nd(ig,:)=0.
156	na(ig,:)=0.
157
158	DO l=1,nlayer
159	! old
160	!! nd(ig,l)=pplay(ig,l)/(cbol*temper)
161	!! nd(ig,l)=pplay(ig,l)/(cbol*pt(ig,l)) ! number of mol density of the layer
162	!! na(ig,l)=na(ig,l)+nd(ig,l) ! density : sum of each layer
163
164	! calculation of density
165	! rho(ig,l)=pplay(ig,l)/(r*temper)
166	rho(ig,l)=pplay(ig,l)/(r*pt(ig,l))
167
168	! calculation of column number density for ch4 (cm-2)
169	na(ig,l)=pq(ig,l,igcm_ch4_gas)/(mmol(igcm_ch4_gas)1e-3) &
170	pplev(ig,l)/gavog1.e-4
171
172	! calculation of atmospheric number density (cm-3) for each layer
173	nd(ig,l)=rho(ig,l)/(mmol(igcm_n2)1.e-3)avog*1.e-6
174
175	! calculation of the phi :
176	phin3(ig,l)=dphin3*nd(ig,l)
177	phi3(ig,l)=dphi3*nd(ig,l)
178	phi2n4(ig,l)=dphi2n4*nd(ig,l)
179	phi3n4(ig,l)=dphi3n4*nd(ig,l)
180	phi4(ig,l)=dphi4*nd(ig,l)
181
182	! prendre en compte Temp
183	! diff entre les gammaf4 et f3 ?
184
185	ENDDO
186
187	! interpolation of na values in cm-2 from Strobel values :
188	CALL interp_line(Naval2,gam2,12,na(ig,:),gammafi23,nlayer)
189	CALL interp_line(Naval3,gam3,14,na(ig,:),gammafi33,nlayer)
190	! CALL interp_line(Naval,gam,20,na(ig,:),gammafi,nlayer)
191	! CALL interp_line(Naval,gamcool,20,na(ig,:),gammacool,nlayer)
192	! CALL interp_line(Navalcool,gam,20,na(ig,:),gammacool,nlayer)
193
194	! write(,) 'TBrad phi4 = ',phi4(1,:)
195	! write(,) 'TBrad gammafi = ',gammafi
196	DO l=1,nlayer
197	! Heating : cf Strobel 1996 and Zalucha 2011
198	! eps 2.3
199	eps23(ig,l)=phin3(ig,l)/(phin3(ig,l)+gammafi23(l))* &
200	(beta23+phi3n4(ig,l)/(phi3n4(ig,l)+1)* &
201	phi4(ig,l)/(phi4(ig,l)+gammafi23(l))* &
202	(3.phi2n4(ig,l))/(1.+phi2n4(ig,l))(1-beta23)/3.)
203
204	! eps 3.3
205	eps33(ig,l)=phi3(ig,l)/(phi3(ig,l)+gammafi33(l))* &
206	(beta33+phi2n4(ig,l)/(phi2n4(ig,l)+1)* &
207	phi4(ig,l)/(phi4(ig,l)+gammafi33(l))* &
208	(1-beta33))
209
210	! Cooling from Zalucha et al 2013 eq10
211	eps76(ig,l)=1./(1.+gammafi23(l)10./(2.phi4(ig,l)))
212
213	ENDDO
214	ENDDO ! ngrid
215	ELSE ! if (strobel)
216	! if strobel = false
217
218	klw = 3.6E-4 ! fit strobel
219	! eps_sw1Pa = 0.6 ! fit strobel
220	eps_sw1Pa = 0.90 ! to fit NH temperature en 1D [ch4]=0.5
221	! eps_sw1Pa = 0.98 ! to fit NH temperature en 3D = fit 1D [ch4]=0.25
222
223	ksw = 8.5E-5*(0.9/(eps_sw1Pa -0.1) -1.)
224	do l=1, nlayer
225	do ig=1, ngrid
226	rho(ig,l)=pplay(ig,l)/(r*pt(ig,l))
227	eps23(ig,l) = 0.1+ 0.9/(1+ ksw/rho(ig,l))
228	eps76(ig,l) = 1/(1+klw/rho(ig,l))
229	end do
230	end do
231	END IF
232	end subroutine nlte_ch4
233

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

Download in other formats:

Original Format