Context Navigation

nirco2abs.F @ 706

Last change on this file since 706 was 690, checked in by acolaitis, 13 years ago
Code re-organization in diverse parts of the GCM code. These are NOT cosmetic changes, but are needed for compilation of the Mesoscale model in NESTED configuration
File size: 8.5 KB

Line
1	SUBROUTINE nirco2abs(ngrid,nlayer,pplay,dist_sol,nq,pq,
2	$ mu0,fract,declin,pdtnirco2)
3
4	IMPLICIT NONE
5	c=======================================================================
6	c subject:
7	c --------
8	c Computing heating rate due to
9	c absorption by CO2 in the near-infrared
10	c This version includes NLTE effects
11	c
12	c (Scheme to be described in Forget et al., JGR, 2003)
13	c (old Scheme described in Forget et al., JGR, 1999)
14	c
15	c This version updated with a new functional fit,
16	c see NLTE correction-factor of Lopez-Valverde et al (1998)
17	c Stephen Lewis 2000
18	c
19	c jul 2011 malv+fgg New corrections for NLTE implemented
20	c 08/2002 : correction for bug when running with diurnal=F
21	c
22	c author: Frederic Hourdin 1996
23	c ------
24	c Francois Forget 1999
25	c
26	c input:
27	c -----
28	c ngrid number of gridpoint of horizontal grid
29	c nlayer Number of layer
30	c dist_sol sun-Mars distance (AU)
31	c mu0(ngridmx)
32	c fract(ngridmx) day fraction of the time interval
33	c declin latitude of subslar point
34	c
35	c output:
36	c -------
37	c
38	c pdtnirco2(ngrid,nlayer) Heating rate (K/s)
39	c
40	c
41	c=======================================================================
42	c
43	c 0. Declarations :
44	c ------------------
45	c
46	#include "dimensions.h"
47	#include "dimphys.h"
48	#include "comcstfi.h"
49	#include "callkeys.h"
50	#include "comdiurn.h"
51	#include "nirdata.h"
52	#include "tracer.h"
53
54	c-----------------------------------------------------------------------
55	c Input/Output
56	c ------------
57	integer,intent(in) :: ngrid ! number of (horizontal) grid points
58	integer,intent(in) :: nlayer ! number of atmospheric layers
59	real,intent(in) :: pplay(ngrid,nlayer) ! Pressure
60	real,intent(in) :: dist_sol ! Sun-Mars distance (in AU)
61	integer,intent(in) :: nq ! number of tracers
62	real,intent(in) :: pq(ngrid,nlayer,nq) ! tracers
63	real,intent(in) :: mu0(ngridmx) ! solar angle
64	real,intent(in) :: fract(ngridmx) ! day fraction of the time interval
65	real,intent(in) :: declin ! latitude of sub-solar point
66
67	real,intent(out) :: pdtnirco2(ngrid,nlayer) ! heating rate (K/s)
68	c
69	c Local variables :
70	c -----------------
71	INTEGER l,ig, n, nstep,i
72	REAL co2heat0, zmu(ngridmx)
73
74	c special diurnal=F
75	real mu0_int(ngridmx),fract_int(ngridmx),zday_int
76	real ztim1,ztim2,ztim3,step
77
78	c
79	c local saved variables
80	c ---------------------
81	logical,save :: firstcall=.true.
82	integer,save :: ico2=0 ! index of "co2" tracer
83	integer,save :: io=0 ! index of "o" tracer
84	c p0noonlte is a pressure below which non LTE effects are significant.
85	c REAL p0nonlte
86	c DATA p0nonlte/7.5e-3/
87	c SAVE p0nonlte
88
89	c parameters for CO2 heating fit
90	real n_a, n_p0, n_b
91	parameter (n_a=1.1956475)
92	parameter (n_b=1.9628251)
93	parameter (n_p0=0.0015888279)
94
95	c Variables added to implement NLTE correction factor (feb 2011)
96	real pyy(nlayer)
97	real cor1(nlayer),oldoco2(nlayer),alfa2(nlayer)
98	real p2011,cociente1,merge
99	real cor0,oco2gcm
100
101	c----------------------------------------------------------------------
102
103	c Initialisation
104	c --------------
105	if (firstcall) then
106	if (nircorr.eq.1) then
107	! we will need co2 and o tracers
108	ico2=igcm_co2
109	if (ico2==0) then
110	write(,) "nirco2abs error: I need a CO2 tracer"
111	write(,) " when running with nircorr==1"
112	stop
113	endif
114	io=igcm_o
115	if (io==0) then
116	write(,) "nirco2abs error: I need an O tracer"
117	write(,) " when running with nircorr==1"
118	stop
119	endif
120	endif
121	firstcall=.false.
122	endif
123
124
125	c co2heat is the heating by CO2 at 700Pa for a zero zenithal angle.
126	co2heat0=n_a(1.52/dist_sol)*2/daysec
127
128	c Simple calcul for a given sun incident angle (if diurnal=T)
129	c --------------------------------------------
130
131	IF (diurnal) THEN
132	do ig=1,ngrid
133	zmu(ig)=sqrt(1224.mu0(ig)mu0(ig)+1.)/35.
134
135	if(nircorr.eq.1) then
136	do l=1,nlayer
137	pyy(l)=pplay(ig,l)
138	enddo
139
140	call interpnir(cor1,pyy,nlayer,corgcm,pres1d,npres)
141	call interpnir(oldoco2,pyy,nlayer,oco21d,pres1d,npres)
142	call interpnir(alfa2,pyy,nlayer,alfa,pres1d,npres)
143	endif
144
145	do l=1,nlayer
146	! Calculations for the O/CO2 correction
147	if(nircorr.eq.1) then
148	cor0=1./(1.+n_p0/pplay(ig,l))**n_b
149	if(pq(ig,l,ico2).gt.1.e-6) then
150	oco2gcm=pq(ig,l,io)/pq(ig,l,ico2)
151	else
152	oco2gcm=1.e6
153	endif
154	cociente1=oco2gcm/oldoco2(l)
155	merge=alog10(cociente1)alfa2(l)+alog10(cor0)
156	$ (1.-alfa2(l))
157	merge=10**merge
158	p2011=sqrt(merge)*cor0
159	else if (nircorr.eq.0) then
160	p2011=1.
161	cor1(l)=1.
162	endif
163
164	if(fract(ig).gt.0.) pdtnirco2(ig,l)=
165	& co2heat0sqrt((700.zmu(ig))/pplay(ig,l))
166	& /(1.+n_p0/pplay(ig,l))**n_b
167	! Corrections from tabulation
168	$ * cor1(l) * p2011
169	c OLD SCHEME (forget et al. 1999)
170	c s co2heat0sqrt((700.zmu(ig))/pplay(ig,l))
171	c s / (1+p0nonlte/pplay(ig,l))
172	enddo
173	enddo
174
175
176	c Averaging over diurnal cycle (if diurnal=F)
177	c -------------------------------------------
178	c NIR CO2 abs is slightly non linear. To remove the diurnal
179	c cycle, it is better to average the heating rate over 1 day rather
180	c than using the mean mu0 computed by mucorr in physiq.F (FF, 1998)
181
182	ELSE ! if (.not.diurnal) then
183
184	nstep = 20 ! number of integration step /sol
185	do n=1,nstep
186	zday_int = (n-1)/float(nstep)
187	ztim2=COS(declin)COS(2.pi*(zday_int-.5))
188	ztim3=-COS(declin)SIN(2.pi*(zday_int-.5))
189	CALL solang(ngrid,sinlon,coslon,sinlat,coslat,
190	s ztim1,ztim2,ztim3,
191	s mu0_int,fract_int)
192	do ig=1,ngrid
193	zmu(ig)=sqrt(1224.mu0_int(ig)mu0_int(ig)+1.)/35.
194
195	if(nircorr.eq.1) then
196	do l=1,nlayer
197	pyy(l)=pplay(ig,l)
198	enddo
199	call interpnir(cor1,pyy,nlayer,corgcm,pres1d,npres)
200	call interpnir(oldoco2,pyy,nlayer,oco21d,pres1d,npres)
201	call interpnir(alfa2,pyy,nlayer,alfa,pres1d,npres)
202	endif
203
204	do l=1,nlayer
205	if(nircorr.eq.1) then
206	cor0=1./(1.+n_p0/pplay(ig,l))**n_b
207	oco2gcm=pq(ig,l,io)/pq(ig,l,ico2)
208	cociente1=oco2gcm/oldoco2(l)
209	merge=alog10(cociente1)alfa2(l)+alog10(cor0)
210	$ (1.-alfa2(l))
211	merge=10**merge
212	p2011=sqrt(merge)*cor0
213	else if (nircorr.eq.0) then
214	p2011=1.
215	cor1(l)=1.
216	endif
217
218	if(fract_int(ig).gt.0.) pdtnirco2(ig,l)=
219	& pdtnirco2(ig,l) + (1/float(nstep))*
220	& co2heat0sqrt((700.zmu(ig))/pplay(ig,l))
221	& /(1.+n_p0/pplay(ig,l))**n_b
222	! Corrections from tabulation
223	$ * cor1(l) * p2011
224	enddo
225	enddo
226	end do
227	END IF
228
229	return
230	end
231
232
233
234	subroutine interpnir(escout,p,nlayer,escin,pin,nl)
235	C
236	C subroutine to perform linear interpolation in pressure from 1D profile
237	C escin(nl) sampled on pressure grid pin(nl) to profile
238	C escout(nlayer) on pressure grid p(nlayer).
239	C
240	real escout(nlayer),p(nlayer)
241	real escin(nl),pin(nl),wm,wp
242	integer nl,nlayer,n1,n,nm,np
243	do n1=1,nlayer
244	if(p(n1) .gt. 1500. .or. p(n1) .lt. 1.0e-13) then
245	escout(n1) = 0.0
246	else
247	do n = 1,nl-1
248	if (p(n1).le.pin(n).and.p(n1).ge.pin(n+1)) then
249	nm=n
250	np=n+1
251	wm=abs(pin(np)-p(n1))/(pin(nm)-pin(np))
252	wp=1.0 - wm
253	endif
254	enddo
255	escout(n1) = escin(nm)wm + escin(np)wp
256	endif
257	enddo
258	return
259	end

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

Download in other formats:

Original Format