Context Navigation

nirco2abs.F @ 1684

Last change on this file since 1684 was 1524, checked in by emillour, 9 years ago

All GCMS:
More updates to enforce dynamics/physics separation:

get rid of references to "temps_mod" from physics packages;
make a "time_phylmdz_mod.F90" module to store that
information and fill it via "iniphysiq".

File size: 8.5 KB

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

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

Download in other formats:

Original Format