Context Navigation

nirco2abs.F @ 2266

Last change on this file since 2266 was 2240, checked in by cmathe, 5 years ago
GCM MARS: Initialization of mu0_int and ztim1 in nirco2abs.F to prevent crash during debug mode simulations.
File size: 8.6 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	! AS: OK firstcall absolute
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	mu0_int(1:ngrid) = 0.
186	ztim1 = 0.
187	do n=1,nstep
188	zday_int = (n-1)/float(nstep)
189	ztim2=COS(declin)COS(2.pi*(zday_int-.5))
190	ztim3=-COS(declin)SIN(2.pi*(zday_int-.5))
191	CALL solang(ngrid,sinlon,coslon,sinlat,coslat,
192	s ztim1,ztim2,ztim3,
193	s mu0_int,fract_int)
194	do ig=1,ngrid
195	zmu(ig)=sqrt(1224.mu0_int(ig)mu0_int(ig)+1.)/35.
196
197	if(nircorr.eq.1) then
198	do l=1,nlayer
199	pyy(l)=pplay(ig,l)
200	enddo
201	call interpnir(cor1,pyy,nlayer,corgcm,pres1d,npres)
202	call interpnir(oldoco2,pyy,nlayer,oco21d,pres1d,npres)
203	call interpnir(alfa2,pyy,nlayer,alfa,pres1d,npres)
204	endif
205
206	do l=1,nlayer
207	if(nircorr.eq.1) then
208	cor0=1./(1.+n_p0/pplay(ig,l))**n_b
209	oco2gcm=pq(ig,l,io)/pq(ig,l,ico2)
210	cociente1=oco2gcm/oldoco2(l)
211	merge=alog10(cociente1)alfa2(l)+alog10(cor0)
212	$ (1.-alfa2(l))
213	merge=10**merge
214	p2011=sqrt(merge)*cor0
215	else if (nircorr.eq.0) then
216	p2011=1.
217	cor1(l)=1.
218	endif
219
220	if(fract_int(ig).gt.0.) pdtnirco2(ig,l)=
221	& pdtnirco2(ig,l) + (1/float(nstep))*
222	& co2heat0sqrt((700.zmu(ig))/pplay(ig,l))
223	& /(1.+n_p0/pplay(ig,l))**n_b
224	! Corrections from tabulation
225	$ * cor1(l) * p2011
226	enddo
227	enddo
228	end do
229	END IF
230
231	return
232	end
233
234
235
236	subroutine interpnir(escout,p,nlayer,escin,pin,nl)
237	C
238	C subroutine to perform linear interpolation in pressure from 1D profile
239	C escin(nl) sampled on pressure grid pin(nl) to profile
240	C escout(nlayer) on pressure grid p(nlayer).
241	C
242	real escout(nlayer),p(nlayer)
243	real escin(nl),pin(nl),wm,wp
244	integer nl,nlayer,n1,n,nm,np
245	do n1=1,nlayer
246	if(p(n1) .gt. 1500. .or. p(n1) .lt. 1.0e-13) then
247	escout(n1) = 0.0
248	else
249	do n = 1,nl-1
250	if (p(n1).le.pin(n).and.p(n1).ge.pin(n+1)) then
251	nm=n
252	np=n+1
253	wm=abs(pin(np)-p(n1))/(pin(nm)-pin(np))
254	wp=1.0 - wm
255	endif
256	enddo
257	escout(n1) = escin(nm)wm + escin(np)wp
258	endif
259	enddo
260	return
261	end

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

Download in other formats:

Original Format