Context Navigation

optcv.F90 @ 867

Last change on this file since 867 was 716, checked in by rwordsworth, 13 years ago

Mainly updates to radiative transfer and gas management scheme.
Most CIA data now read from standard HITRAN datafiles. For the H2O
continuum, two options have been added: the standard CKD continuum,
and the empirical formula in PPC (Pierrehumbert 2010). Use the toggle
'H2Ocont_simple' in callphys.def to choose.

Note to Martians: I've changed the default values of 'sedimentation' and
'co2cond' in inifis to false. Both these are defined in the standard deftank
callphys.def file, so there shouldn't be any compatibility problems.

Property svn:executable set to *

File size: 10.7 KB

Line
1	SUBROUTINE OPTCV(DTAUV,TAUV,TAUCUMV,PLEV, &
2	QXVAER,QSVAER,GVAER,WBARV,COSBV, &
3	TAURAY,TAUAERO,TMID,PMID,TAUGSURF,QVAR,MUVAR)
4
5	use radinc_h
6	use radcommon_h, only: gasv, tlimit, wrefVAR, Cmk, tgasref, pfgasref,wnov,scalep
7	use gases_h
8
9	implicit none
10
11	!==================================================================
12	!
13	! Purpose
14	! -------
15	! Calculates shortwave optical constants at each level.
16	!
17	! Authors
18	! -------
19	! Adapted from the NASA Ames code by R. Wordsworth (2009)
20	!
21	!==================================================================
22	!
23	! THIS SUBROUTINE SETS THE OPTICAL CONSTANTS IN THE VISUAL
24	! IT CALCUALTES FOR EACH LAYER, FOR EACH SPECRAL INTERVAL IN THE VISUAL
25	! LAYER: WBAR, DTAU, COSBAR
26	! LEVEL: TAU
27	!
28	! TAUV(L,NW,NG) is the cumulative optical depth at the top of radiation code
29	! layer L. NW is spectral wavelength interval, ng the Gauss point index.
30	!
31	! TLEV(L) - Temperature at the layer boundary
32	! PLEV(L) - Pressure at the layer boundary (i.e. level)
33	! GASV(NT,NPS,NW,NG) - Visible k-coefficients
34	!
35	!-------------------------------------------------------------------
36
37	#include "callkeys.h"
38	#include "comcstfi.h"
39
40	real*8 DTAUV(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
41	real*8 DTAUKV(L_LEVELS+1,L_NSPECTV,L_NGAUSS)
42	real*8 TAUV(L_NLEVRAD,L_NSPECTV,L_NGAUSS)
43	real*8 TAUCUMV(L_LEVELS,L_NSPECTV,L_NGAUSS)
44	real*8 PLEV(L_LEVELS)
45	real*8 TMID(L_LEVELS), PMID(L_LEVELS)
46	real*8 COSBV(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
47	real*8 WBARV(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
48	real*8 TAURAY(L_NSPECTV)
49
50	! for aerosols
51	real*8 QXVAER(L_LEVELS+1,L_NSPECTV,NAERKIND)
52	real*8 QSVAER(L_LEVELS+1,L_NSPECTV,NAERKIND)
53	real*8 GVAER(L_LEVELS+1,L_NSPECTV,NAERKIND)
54	real*8 TAUAERO(L_LEVELS+1,NAERKIND)
55	real*8 TAUAEROLK(L_LEVELS+1,L_NSPECTV,NAERKIND)
56	real*8 TAEROS(L_LEVELS,L_NSPECTV,NAERKIND)
57
58	integer L, NW, NG, K, NG1(L_NSPECTV), LK, IAER
59	integer MT(L_LEVELS), MP(L_LEVELS), NP(L_LEVELS)
60	real*8 ANS, TAUGAS
61	real*8 TRAY(L_LEVELS,L_NSPECTV)
62	real*8 DPR(L_LEVELS), U(L_LEVELS)
63	real*8 LCOEF(4), LKCOEF(L_LEVELS,4)
64
65	real*8 taugsurf(L_NSPECTV,L_NGAUSS-1), TRAYAER
66
67	! variable species mixing ratio variables
68	real*8 QVAR(L_LEVELS), WRATIO(L_LEVELS), MUVAR(L_LEVELS)
69	real*8 KCOEF(4)
70	integer NVAR(L_LEVELS)
71
72	! temporary variables for multiple aerosol calculation
73	real*8 atemp, btemp, ctemp
74
75	! variables for k in units m^-1
76	double precision wn_cont, p_cont, p_air, T_cont, dtemp
77	double precision p_cross
78	real*8 dz(L_LEVELS), DCONT
79
80	integer igas, jgas
81
82	!=======================================================================
83	! Determine the total gas opacity throughout the column, for each
84	! spectral interval, NW, and each Gauss point, NG.
85	! Calculate the continuum opacities, i.e., those that do not depend on
86	! NG, the Gauss index.
87
88	taugsurf(:,:) = 0.0
89	dpr(:) = 0.0
90	lkcoef(:,:) = 0.0
91
92	do K=2,L_LEVELS
93	DPR(k) = PLEV(K)-PLEV(K-1)
94
95	! if we have continuum opacities, we need dz
96	if(kastprof)then
97	dz(k) = dpr(k)(8314.5/muvar(k))TMID(K)/(g*PMID(K))
98	U(k) = (Cmkmugaz/(muvar(k)))DPR(k)
99	else
100	dz(k) = dpr(k)RTMID(K)/(g*PMID(K))
101	U(k) = Cmk*DPR(k) ! only Cmk line in optci.F
102	endif
103
104	call tpindex(PMID(K),TMID(K),QVAR(K),pfgasref,tgasref,WREFVAR, &
105	LCOEF,MT(K),MP(K),NVAR(K),WRATIO(K))
106
107	do LK=1,4
108	LKCOEF(K,LK) = LCOEF(LK)
109	end do
110
111	DO NW=1,L_NSPECTV
112	TRAY(K,NW) = TAURAY(NW) * DPR(K)
113
114	do iaer=1,naerkind
115	TAEROS(K,NW,IAER) = TAUAERO(K,IAER) * QXVAER(K,NW,IAER)
116	end do
117
118	END DO
119	end do
120
121	! TRAYAER is Tau RAYleigh scattering, plus AERosol opacity
122
123	! we ignore K=1...
124	do K=2,L_LEVELS
125	do NW=1,L_NSPECTV
126
127	TRAYAER = TRAY(K,NW)
128	do iaer=1,naerkind
129	TRAYAER = TRAYAER + TAEROS(K,NW,IAER)
130	end do
131
132	DCONT = 0.0 ! continuum absorption
133
134	if(callgasvis.and.continuum.and.(.not.graybody))then
135	! include continua if necessary
136	wn_cont = dble(wnov(nw))
137	T_cont = dble(TMID(k))
138	do igas=1,ngasmx
139
140	if(gfrac(igas).eq.-1)then ! variable
141	p_cont = dble(PMID(k)scalepQVAR(k)) ! qvar = mol/mol
142	else
143	p_cont = dble(PMID(k)scalepgfrac(igas)*(1.-QVAR(k)))
144	endif
145
146	dtemp=0.0
147	if(igas.eq.igas_N2)then
148
149	!call interpolateN2N2(wn_cont,T_cont,p_cont,dtemp,.false.)
150	! only goes to 500 cm^-1, so unless we're around a cold brown dwarf, this is irrelevant in the visible
151
152	elseif(igas.eq.igas_H2)then
153
154	! first do self-induced absorption
155	call interpolateH2H2(wn_cont,T_cont,p_cont,dtemp,.false.)
156
157	! then cross-interactions with other gases
158	do jgas=1,ngasmx
159	p_cross = dble(PMID(k)scalepgfrac(jgas)*(1.-QVAR(k)))
160	if(jgas.eq.igas_N2)then
161	call interpolateN2H2(wn_cont,T_cont,p_cross,p_cont,dtemp,.false.)
162	! should be irrelevant in the visible
163	elseif(jgas.eq.igas_He)then
164	call interpolateH2He(wn_cont,T_cont,p_cross,p_cont,dtemp,.false.)
165	endif
166	enddo
167
168	elseif(igas.eq.igas_H2O.and.T_cont.gt.200.0)then
169
170	p_air = dble(PMID(k)*scalep) - p_cont ! note assumes background is air!
171	if(H2Ocont_simple)then
172	call interpolateH2Ocont_PPC(wn_cont,T_cont,p_cont,p_air,dtemp,.false.)
173	else
174	call interpolateH2Ocont_CKD(wn_cont,T_cont,p_cont,p_air,dtemp,.false.)
175	endif
176
177	endif
178
179	DCONT = DCONT + dtemp
180
181	enddo
182
183	DCONT = DCONT*dz(k)
184	endif
185
186	do NG=1,L_NGAUSS-1
187
188	!=======================================================================
189	! Now compute TAUGAS
190	! Interpolate between water mixing ratios
191	! WRATIO = 0.0 if the requested water amount is equal to, or outside the
192	! the water data range
193
194	if (L_REFVAR.eq.1)then ! added by RW for special no variable case
195	KCOEF(1) = GASV(MT(K),MP(K),1,NW,NG)
196	KCOEF(2) = GASV(MT(K),MP(K)+1,1,NW,NG)
197	KCOEF(3) = GASV(MT(K)+1,MP(K)+1,1,NW,NG)
198	KCOEF(4) = GASV(MT(K)+1,MP(K),1,NW,NG)
199	else
200	KCOEF(1) = GASV(MT(K),MP(K),NVAR(K),NW,NG) + WRATIO(K)* &
201	(GASV(MT(K),MP(K),NVAR(K)+1,NW,NG) - &
202	GASV(MT(K),MP(K),NVAR(K),NW,NG))
203
204	KCOEF(2) = GASV(MT(K),MP(K)+1,NVAR(K),NW,NG) + WRATIO(K)* &
205	(GASV(MT(K),MP(K)+1,NVAR(K)+1,NW,NG) - &
206	GASV(MT(K),MP(K)+1,NVAR(K),NW,NG))
207
208	KCOEF(3) = GASV(MT(K)+1,MP(K)+1,NVAR(K),NW,NG) + WRATIO(K)*&
209	(GASV(MT(K)+1,MP(K)+1,NVAR(K)+1,NW,NG) - &
210	GASV(MT(K)+1,MP(K)+1,NVAR(K),NW,NG))
211
212	KCOEF(4) = GASV(MT(K)+1,MP(K),NVAR(K),NW,NG) + WRATIO(K)* &
213	(GASV(MT(K)+1,MP(K),NVAR(K)+1,NW,NG) - &
214	GASV(MT(K)+1,MP(K),NVAR(K),NW,NG))
215	endif
216
217	! Interpolate the gaseous k-coefficients to the requested T,P values
218
219	ANS = LKCOEF(K,1)KCOEF(1) + LKCOEF(K,2)KCOEF(2) + &
220	LKCOEF(K,3)KCOEF(3) + LKCOEF(K,4)KCOEF(4)
221
222	TAUGAS = U(k)*ANS
223
224	!TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS
225	TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS + DCONT
226	DTAUKV(K,nw,ng) = TAUGAS + TRAYAER & ! TRAYAER includes all scattering contributions
227	+ DCONT ! for continuum absorption
228
229	end do
230
231
232	! Now fill in the "clear" part of the spectrum (NG = L_NGAUSS),
233	! which holds continuum opacity only
234
235	NG = L_NGAUSS
236	DTAUKV(K,nw,ng) = TRAY(K,NW) + DCONT ! For parameterized continuum absorption
237	do iaer=1,naerkind
238	DTAUKV(K,nw,ng) = DTAUKV(K,nw,ng) + TAEROS(K,NW,IAER)
239	! & + DCONT ! For parameterized continuum absorption
240	end do ! a bug was here!
241
242	end do
243	end do
244
245
246	!=======================================================================
247	! Now the full treatment for the layers, where besides the opacity
248	! we need to calculate the scattering albedo and asymmetry factors
249
250	DO NW=1,L_NSPECTV
251	DO K=2,L_LEVELS
252	do iaer=1,naerkind
253	TAUAEROLK(K,NW,IAER) = TAUAERO(K,IAER) * QSVAER(K,NW,IAER)
254	end do
255	ENDDO
256	ENDDO
257
258
259	DO NW=1,L_NSPECTV
260	DO NG=1,L_NGAUSS
261	DO L=1,L_NLAYRAD-1
262	K = 2*L+1
263
264	DTAUV(L,nw,ng) = DTAUKV(K,NW,NG)+DTAUKV(K+1,NW,NG)
265
266	atemp=0.
267	btemp=TRAY(K,NW) + TRAY(K+1,NW)
268	ctemp=0.9999*(TRAY(K,NW) + TRAY(K+1,NW))
269	do iaer=1,naerkind
270	atemp = atemp + &
271	GVAER(K,NW,IAER) * TAUAEROLK(K,NW,IAER) + &
272	GVAER(K+1,NW,IAER) * TAUAEROLK(K+1,NW,IAER)
273	btemp = btemp + &
274	TAUAEROLK(K,NW,IAER) + TAUAEROLK(K+1,NW,IAER)
275	ctemp = ctemp + &
276	TAUAEROLK(K,NW,IAER) + TAUAEROLK(K+1,NW,IAER)
277	end do
278
279	COSBV(L,NW,NG) = atemp/btemp
280	WBARV(L,nw,ng) = ctemp/DTAUV(L,nw,ng)
281
282	END DO
283
284	! No vertical averaging on bottom layer
285
286	L = L_NLAYRAD
287	K = 2*L+1
288	DTAUV(L,nw,ng) = DTAUKV(K,NW,NG)
289
290	atemp=0.
291	btemp=TRAY(K,NW)
292	ctemp=0.9999*TRAY(K,NW)
293	do iaer=1,naerkind
294	atemp = atemp + GVAER(K,NW,IAER) * TAUAEROLK(K,NW,IAER)
295	btemp = btemp + TAUAEROLK(K,NW,IAER)
296	ctemp = ctemp + TAUAEROLK(K,NW,IAER)
297	end do
298	COSBV(L,NW,NG) = atemp/btemp
299	WBARV(L,nw,ng) = ctemp/DTAUV(L,nw,ng)
300
301	END DO ! NG gauss point loop
302	END DO ! NW spectral loop
303
304
305
306	! Total extinction optical depths
307
308	DO NW=1,L_NSPECTV
309	DO NG=1,L_NGAUSS ! full gauss loop
310	TAUV(1,NW,NG)=0.0D0
311	DO L=1,L_NLAYRAD
312	TAUV(L+1,NW,NG)=TAUV(L,NW,NG)+DTAUV(L,NW,NG)
313	END DO
314
315	TAUCUMV(1,NW,NG)=0.0D0
316	DO K=2,L_LEVELS
317	TAUCUMV(K,NW,NG)=TAUCUMV(K-1,NW,NG)+DTAUKV(K,NW,NG)
318	END DO
319	END DO ! end full gauss loop
320	END DO
321
322
323	RETURN
324	END SUBROUTINE OPTCV

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

Download in other formats:

Original Format