Context Navigation

nirco2abs.F @ 1692

Last change on this file since 1692 was 1545, checked in by emillour, 10 years ago

Venus and Titan GCMs:

Adaptation wrt previous changes for Titan and Venus where
longitude and latitude arrays (in phycommon/geometry_mod) were overwritten
with values from startphy.nc files, where values are given in degrees.
For the sake of homegeneity with other physics package, revert to "default"
behaviour: longitude/latitude are in radians and longitude_deg/latitude_deg
are in degrees.
Also added checking in phyetat0 that the longitude/latitude read in the
restartphy.nc files match the ones provided by the dynamics.

File size: 9.4 KB

Line
1	SUBROUTINE nirco2abs(nlon,nlev,nplay,dist_sol,nq,pq,
2	$ mu0,fract,pdtnirco2)
3
4	use dimphy
5	use geometry_mod, only: longitude_deg, latitude_deg
6	use chemparam_mod, only: i_co2, i_o
7	c use compo_hedin83_mod2
8
9
10	IMPLICIT NONE
11	c=======================================================================
12	c subject:
13	c --------
14	c Computing heating rate due to
15	c absorption by CO2 in the near-infrared
16	c This version includes NLTE effects
17	c
18	c (Scheme to be described in Forget et al., JGR, 2003)
19	c (old Scheme described in Forget et al., JGR, 1999)
20	c
21	c This version updated with a new functional fit,
22	c see NLTE correction-factor of Lopez-Valverde et al (1998)
23	c Stephen Lewis 2000
24	c
25	c oct 2014 g.gilli Coupling with photochemical model
26	C jan 2014 g.gilli Revision (following martian non-lte param)
27	C jun 2013 l.salmi First adaptation to Venus and NIR NLTE param
28	c jul 2011 malv+fgg New corrections for NLTE implemented
29	c 08/2002 : correction for bug when running with diurnal=F
30	c
31	c author: Frederic Hourdin 1996
32	c ------
33	c Francois Forget 1999
34	c
35	c input:
36	c -----
37	c nlon number of gridpoint of horizontal grid
38	c nlev Number of layer
39	c dist_sol sun-Venus distance (AU)
40	c mu0(nlon)
41	c fract(nlon) day fraction of the time interval
42	c declin latitude of subslar point
43	c
44	c output:
45	c -------
46	c
47	c pdtnirco2(nlon,nlev) Heating rate (K/sec)
48	c
49	c
50	c=======================================================================
51	c
52	c 0. Declarations :
53	c ------------------
54	c
55
56	#include "YOMCST.h"
57	#include "clesphys.h"
58	c#include "comdiurn.h"
59	#include "nirdata.h"
60	c#include "tracer.h"
61	#include "mmol.h"
62	c-----------------------------------------------------------------------
63	c Input/Output
64	c ------------
65	integer,intent(in) :: nlon ! number of (horizontal) grid points
66	integer,intent(in) :: nlev ! number of atmospheric layers
67
68	real,intent(in) :: nplay(nlon,nlev) ! Pressure
69	real,intent(in) :: dist_sol ! Sun-Venus distance (in AU)
70	integer,intent(in) :: nq ! number of tracers
71	real,intent(in) :: pq(nlon,nlev,nq) ! mass mixing ratio tracers
72	real,intent(in) :: mu0(nlon) ! solar angle
73	real,intent(in) :: fract(nlon) ! day fraction of the time interval
74	c real,intent(in) :: declin ! latitude of sub-solar point
75	real :: co2vmr_gcm(nlon,nlev), o3pvmr_gcm(nlon,nlev)
76
77	real,intent(out) :: pdtnirco2(nlon,nlev) ! heating rate (K/sec)
78
79	c
80	c Local variables :
81	c -----------------
82	INTEGER l,ig, n, nstep,i
83	REAL co2heat0, zmu(nlon)
84
85	c special diurnal=F
86	real mu0_int(nlon),fract_int(nlon),zday_int
87	real ztim1,ztim2,ztim3,step
88
89	c
90	c local saved variables
91	c ---------------------
92	logical,save :: firstcall=.true.
93	integer,save :: ico2=0 ! index of "co2" tracer
94	integer,save :: io=0 ! index of "o" tracer
95
96	cccc parameters for CO2 heating fit
97	c
98	c n_a = heating rate for Venusian day at p0, r0, mu =0 [K day-1]
99	c Here p0 = p_cloud top [Pa]
100	c n_p0 = is a pressure below which non LTE effects are significant [Pa]
101	c n_a Solar heating [K/Eday] at the cloud top, taken from Crisps table
102
103	real n_a, n_p0, n_b, p_ctop
104
105
106	cc "Nominal" values
107	parameter (n_a = 18.13/86400.0) !c K/Eday ---> K/sec
108	parameter (p_ctop=13.2e2)
109	c -- NLTE Param v3 --
110	parameter (n_p0=0.008)
111	parameter (n_b=1.362)
112
113	cc -- Varoxy5
114	C parameter (n_a = 20/86400.0)
115	C parameter (p_ctop=870) ! [Pa]
116	C parameter (n_b=1.98)
117	C parameter (n_p0=0.045)
118
119	c parameter (n_p0=0.1) !!!!cccc test varoxy5mod
120	c parameter (n_b=0.9)
121
122
123	c -- NLTE Param v2 --
124	C parameter (n_p0=0.01)
125	c parameter (n_b = 1.3)
126
127
128
129	c Variables added to implement NLTE correction factor (feb 2011)
130	real pyy(nlev)
131	real cor1(nlev),oldoco2(nlev),alfa2(nlev)
132	real p2011,cociente1,merge
133	real cor0,oco2gcm
134	!!!!
135	c real :: pic27(nlon,nlev), pic27b(nlon,nlev)
136	c real :: pic43(nlon,nlev), picnir(nlon,nlev)
137
138	c co2heat is the heating by CO2 at p_ctop=13.2e2 for a zero zenithal angle.
139
140	co2heat0=n_a(0.72/dist_sol)*2
141
142	CCCCCC TEST: reduce by X% nir Heating
143
144	c co2heat0 = co2heat0 * 0.8
145
146
147	c----------------------------------------------------------------------
148
149	c Initialisation
150	c --------------
151	if (firstcall) then
152	if (nircorr.eq.1) then
153	c ! we will need co2 and o tracers
154	ico2= i_co2
155	if (ico2==0) then
156	write(,) "nirco2abs error: I need a CO2 tracer"
157	write(,) " when running with nircorr==1"
158	stop
159	endif
160	io=i_o
161	if (io==0) then
162	write(,) "nirco2abs error: I need an O tracer"
163	write(,) " when running with nircorr==1"
164	stop
165	endif
166	endif
167	firstcall=.false.
168	endif
169
170
171	c
172	c Simple calcul for a given sun incident angle (if cycle_diurne=T)
173	c --------------------------------------------
174
175	IF (cycle_diurne) THEN
176
177	do ig=1,nlon
178	zmu(ig)=sqrt(1224.mu0(ig)mu0(ig)+1.)/35.
179
180
181	if(nircorr.eq.1) then
182	do l=1,nlev
183	pyy(l)=nplay(ig,l)
184	enddo
185
186	call interpnir(cor1,pyy,nlev,corgcm,pres1d,npres)
187	call interpnir(oldoco2,pyy,nlev,oco21d,pres1d,npres)
188	call interpnir(alfa2,pyy,nlev,alfa,pres1d,npres)
189
190	endif
191
192	do l=1,nlev
193
194	c Calculations for the O/CO2 correction
195	if(nircorr.eq.1) then
196	cor0=1./(1.+n_p0/nplay(ig,l))**n_b
197	if(pq(ig,l,ico2) .gt. 1.e-6) then
198	oco2gcm=pq(ig,l,io)/pq(ig,l,ico2)
199
200	else
201	oco2gcm=1.e6
202	endif
203	cociente1=oco2gcm/oldoco2(l)
204
205	c WRITE(,) "nirco2abs line 211", l, cociente1
206
207	merge=alog10(cociente1)alfa2(l)+alog10(cor0)
208	$ (1.-alfa2(l))
209	merge=10**merge
210	p2011=sqrt(merge)*cor0
211
212	else if (nircorr.eq.0) then
213	p2011=1.
214	cor1(l)=1.
215	endif
216
217	if(fract(ig).gt.0.) pdtnirco2(ig,l)=
218	& co2heat0sqrt((p_ctopzmu(ig))/nplay(ig,l))
219	& /(1.+n_p0/nplay(ig,l))**n_b
220	c Corrections from tabulation
221	$ * cor1(l) * p2011
222
223	enddo
224	enddo
225
226	c Averaging over diurnal cycle (if diurnal=F)
227	c -------------------------------------------
228	c NIR CO2 abs is slightly non linear. To remove the diurnal
229	c cycle, it is better to average the heating rate over 1 day rather
230	c than using the mean mu0 computed by mucorr in physiq.F (FF, 1998)
231
232	ELSE ! if (.not.diurnal) then
233	nstep = 20 ! number of integration step /sol
234	do n=1,nstep
235
236	zday_int = (n-1)/float(nstep)
237
238	CALL zenang(0.,zday_int,RDAY/nstep,
239	& latitude_deg,longitude_deg,
240	& mu0_int,fract_int)
241
242	do ig=1,nlon
243	zmu(ig)=sqrt(1224.mu0_int(ig)mu0_int(ig)+1.)/35.
244
245	if(nircorr.eq.1) then
246	do l=1,nlev
247	pyy(l)=nplay(ig,l)
248	enddo
249	call interpnir(cor1,pyy,nlev,corgcm,pres1d,npres)
250	call interpnir(oldoco2,pyy,nlev,oco21d,pres1d,npres)
251	call interpnir(alfa2,pyy,nlev,alfa,pres1d,npres)
252	endif
253	c
254
255	do l=1,nlev
256	c Calculations for the O/CO2 correction
257	if(nircorr.eq.1) then
258	cor0=1./(1.+n_p0/nplay(ig,l))**n_b
259	oco2gcm=pq(ig,l,io)/pq(ig,l,ico2)
260	cociente1=oco2gcm/oldoco2(l)
261	merge=alog10(cociente1)alfa2(l)+alog10(cor0)
262	$ (1.-alfa2(l))
263	merge=10**merge
264	p2011=sqrt(merge)*cor0
265
266	else if (nircorr.eq.0) then
267	p2011=1.
268	cor1(l)=1.
269	endif
270
271	if(fract_int(ig).gt.0.) pdtnirco2(ig,l)=
272	& pdtnirco2(ig,l) + (1/float(nstep))*
273	& co2heat0sqrt((p_ctopzmu(ig))/nplay(ig,l))
274	& /(1.+n_p0/nplay(ig,l))**n_b
275	! Corrections from tabulation
276	$ * cor1(l) * p2011
277
278	enddo
279	enddo
280	end do
281
282
283	END IF
284
285	return
286	end
287
288
289	subroutine interpnir(escout,p,nlev,escin,pin,nl)
290	C
291	C subroutine to perform linear interpolation in pressure from 1D profile
292	C escin(nl) sampled on pressure grid pin(nl) to profile
293	C escout(nlev) on pressure grid p(nlev).
294	C
295	real escout(nlev),p(nlev)
296	real escin(nl),pin(nl),wm,wp
297	integer nl,nlev,n1,n,nm,np
298	do n1=1,nlev
299	if(p(n1) .gt. 1500. .or. p(n1) .lt. 1.0e-13) then
300	c escout(n1) = 0.0
301	escout(n1) = 1.e-15
302	else
303	do n = 1,nl-1
304	if (p(n1).le.pin(n).and.p(n1).ge.pin(n+1)) then
305	nm=n
306	np=n+1
307	wm=abs(pin(np)-p(n1))/(pin(nm)-pin(np))
308	wp=1.0 - wm
309	endif
310	enddo
311	escout(n1) = escin(nm)wm + escin(np)wp
312	endif
313	enddo
314	return
315	end

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

Download in other formats:

Original Format