Context Navigation

optci.F90 @ 374

Last change on this file since 374 was 305, checked in by rwordsworth, 13 years ago
Several new files added as part of the climate evolution model (main program kcm.F90). Some general cleanup in physiq.F90 and callcorrk.F90. Bugs in dust radiative transfer and H2 Rayleigh scattering corrected.
File size: 11.4 KB

Line
1	subroutine optci(PLEV,TLEV,DTAUI,TAUCUMI, &
2	QXIAER,QSIAER,GIAER,COSBI,WBARI,TAUAERO, &
3	TMID,PMID,TAUGSURF,QVAR,MUVAR)
4
5	use radinc_h
6	use radcommon_h, only: gasi, tlimit, wrefVAR, Cmk,tgasref,pfgasref,wnoi,scalep
7	implicit none
8
9	!==================================================================
10	!
11	! Purpose
12	! -------
13	! Calculates longwave optical constants at each level. For each
14	! layer and spectral interval in the IR it calculates WBAR, DTAU
15	! and COSBAR. For each level it calculates TAU.
16	!
17	! TAUI(L,LW) is the cumulative optical depth at level L (or alternatively
18	! at the bottom of layer L), LW is the spectral wavelength interval.
19	!
20	! TLEV(L) - Temperature at the layer boundary (i.e., level)
21	! PLEV(L) - Pressure at the layer boundary (i.e., level)
22	!
23	! Authors
24	! -------
25	! Adapted from the NASA Ames code by R. Wordsworth (2009)
26	!
27	!==================================================================
28
29
30	#include "comcstfi.h"
31	#include "callkeys.h"
32	#include "gases.h"
33
34
35	real*8 DTAUI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
36	real*8 DTAUKI(L_LEVELS+1,L_NSPECTI,L_NGAUSS)
37	real*8 TAUI(L_NLEVRAD,L_NSPECTI,L_NGAUSS)
38	real*8 TAUCUMI(L_LEVELS,L_NSPECTI,L_NGAUSS)
39	real*8 PLEV(L_LEVELS)
40	real*8 TLEV(L_LEVELS)
41	real*8 TMID(L_LEVELS), PMID(L_LEVELS)
42	real*8 COSBI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
43	real*8 WBARI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
44
45	! For aerosols
46	real*8 QXIAER(L_LEVELS+1,L_NSPECTI,NAERKIND)
47	real*8 QSIAER(L_LEVELS+1,L_NSPECTI,NAERKIND)
48	real*8 GIAER(L_LEVELS+1,L_NSPECTI,NAERKIND)
49	real*8 TAUAERO(L_LEVELS+1,NAERKIND)
50	real*8 TAUAEROLK(L_LEVELS+1,L_NSPECTI,NAERKIND)
51	real*8 TAEROS(L_LEVELS,L_NSPECTI,NAERKIND)
52
53	integer L, NW, NG, K, LK, IAER
54	integer MT(L_LEVELS), MP(L_LEVELS), NP(L_LEVELS)
55	real*8 ANS, TAUGAS
56	real*8 DPR(L_LEVELS), U(L_LEVELS)
57	real*8 LCOEF(4), LKCOEF(L_LEVELS,4)
58
59	real*8 taugsurf(L_NSPECTI,L_NGAUSS-1)
60	real*8 DCONT
61	double precision wn_cont, p_cont, p_air, T_cont, dtemp
62
63	! Variable species mixing ratio variables
64	real*8 QVAR(L_LEVELS), WRATIO(L_LEVELS), MUVAR(L_LEVELS)
65	real*8 KCOEF(4)
66	integer NVAR(L_LEVELS)
67
68	! temporary variables for multiple aerosol calculation
69	real*8 atemp, btemp
70
71	! variables for k in units m^-1
72	real*8 rho, dz(L_LEVELS)
73
74	integer igas
75
76	!--- Kasting's CIA ----------------------------------------
77	!real*8, parameter :: Ci(L_NSPECTI)=[ &
78	! 3.8E-5, 1.2E-5, 2.8E-6, 7.6E-7, 4.5E-7, 2.3E-7, &
79	! 5.4E-7, 1.6E-6, 0.0, &
80	! 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
81	! 0.0, 4.0E-7, 4.0E-6, 1.4E-5, &
82	! 1.0E-5, 1.2E-6, 2.0E-7, 5.0E-8, 3.0E-8, 0.0 ]
83	!real*8, parameter :: Ti(L_NSPECTI)=[ -2.2, -1.9, &
84	! -1.7, -1.7, -1.7, -1.7, -1.7, -1.7, &
85	! 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
86	! -1.7,-1.7,-1.7,-1.7,-1.7,-1.7,-1.7, -1.7,0.0 ]
87	!----------------------------------------------------------
88
89	!=======================================================================
90	! Determine the total gas opacity throughout the column, for each
91	! spectral interval, NW, and each Gauss point, NG.
92
93	taugsurf(:,:) = 0.0
94	dpr(:) = 0.0
95	lkcoef(:,:) = 0.0
96
97	do K=2,L_LEVELS
98	DPR(k) = PLEV(K)-PLEV(K-1)
99
100	!--- Kasting's CIA ----------------------------------------
101	!dz(k)=dpr(k)189.02TMID(K)/(0.03720*PMID(K))
102	! this is CO2 path length (in cm) as written by Francois
103	! delta_z = delta_p * R_specific * T / (g * P)
104	! But Kasting states that W is in units of _atmosphere_ cm
105	! So we do
106	!dz(k)=dz(k)*(PMID(K)/1013.25)
107	!dz(k)=dz(k)/100.0 ! in m for SI calc
108	!----------------------------------------------------------
109
110	! if we have continuum opacities, we need dz
111	if(kastprof)then
112	dz(k) = dpr(k)(8314.5/muvar(k))TMID(K)/(g*PMID(K))
113	U(k) = (Cmkmugaz/(muvar(k)))DPR(k)
114	else
115	dz(k) = dpr(k)RTMID(K)/(g*PMID(K))
116	U(k) = Cmk*DPR(k) ! only Cmk line in optci.F
117	endif
118
119	call tpindex(PMID(K),TMID(K),QVAR(K),pfgasref,tgasref,WREFVAR, &
120	LCOEF,MT(K),MP(K),NVAR(K),WRATIO(K))
121
122	do LK=1,4
123	LKCOEF(K,LK) = LCOEF(LK)
124	end do
125
126
127	DO NW=1,L_NSPECTI
128	do iaer=1,naerkind
129	TAEROS(K,NW,IAER) = TAUAERO(K,IAER) * QXIAER(K,NW,IAER)
130	end do
131	END DO
132	end do ! levels
133
134	do K=2,L_LEVELS
135	do nw=1,L_NSPECTI
136
137	DCONT = 0.0 ! continuum absorption
138
139	! include continua if necessary
140	wn_cont = dble(wnoi(nw))
141	T_cont = dble(TMID(k))
142	do igas=1,ngasmx
143
144	if(gfrac(igas).eq.-1)then ! variable
145	p_cont = dble(PMID(k)scalepQVAR(k)) ! qvar = mol/mol
146	else
147	p_cont = dble(PMID(k)scalepgfrac(igas)*(1.-QVAR(k)))
148	endif
149
150	dtemp=0.0
151	if(gnom(igas).eq.'H2_')then
152	call interpolateH2H2(wn_cont,T_cont,p_cont,dtemp,.false.)
153	elseif(gnom(igas).eq.'H2O'.and.T_cont.gt.200.0)then
154	p_air = dble(PMID(k)*scalep) - p_cont ! note assumes air!!
155	call interpolateH2Ocont(wn_cont,T_cont,p_cont,p_air,dtemp,.false.)
156
157	endif
158
159	DCONT = DCONT + dtemp
160
161	enddo
162
163
164	DCONT = DCONT*dz(k)
165
166	!--- Kasting's CIA ----------------------------------------
167	!DCO2 = dz(k)Ci(nw)(1.2859PMID(k)/1000.0)(TMID(k)/300.)**Ti(nw)
168	!DCO2 = 130Ci(nw)(pmid(k)/1013.25)*2(tmid(k)/300.)*Ti(nw) dz(k)
169	! these two have been verified to give the same results
170	!----------------------------------------------------------
171
172
173	do ng=1,L_NGAUSS-1
174
175	! Now compute TAUGAS
176
177	! Interpolate between water mixing ratios
178	! WRATIO = 0.0 if the requested water amount is equal to, or outside the
179	! the water data range
180
181	if(L_REFVAR.eq.1)then ! added by RW for special no variable case
182	KCOEF(1) = GASI(MT(K),MP(K),1,NW,NG)
183	KCOEF(2) = GASI(MT(K),MP(K)+1,1,NW,NG)
184	KCOEF(3) = GASI(MT(K)+1,MP(K)+1,1,NW,NG)
185	KCOEF(4) = GASI(MT(K)+1,MP(K),1,NW,NG)
186	else
187
188	KCOEF(1) = GASI(MT(K),MP(K),NVAR(K),NW,NG) + WRATIO(K)* &
189	(GASI(MT(K),MP(K),NVAR(K)+1,NW,NG) - &
190	GASI(MT(K),MP(K),NVAR(K),NW,NG))
191
192	KCOEF(2) = GASI(MT(K),MP(K)+1,NVAR(K),NW,NG) + WRATIO(K)* &
193	(GASI(MT(K),MP(K)+1,NVAR(K)+1,NW,NG) - &
194	GASI(MT(K),MP(K)+1,NVAR(K),NW,NG))
195
196	KCOEF(3) = GASI(MT(K)+1,MP(K)+1,NVAR(K),NW,NG) + WRATIO(K)* &
197	(GASI(MT(K)+1,MP(K)+1,NVAR(K)+1,NW,NG) - &
198	GASI(MT(K)+1,MP(K)+1,NVAR(K),NW,NG))
199
200	KCOEF(4) = GASI(MT(K)+1,MP(K),NVAR(K),NW,NG) + WRATIO(K)* &
201	(GASI(MT(K)+1,MP(K),NVAR(K)+1,NW,NG) - &
202	GASI(MT(K)+1,MP(K),NVAR(K),NW,NG))
203	endif
204
205	! Interpolate the gaseous k-coefficients to the requested T,P values
206
207	ANS = LKCOEF(K,1)KCOEF(1) + LKCOEF(K,2)KCOEF(2) + &
208	LKCOEF(K,3)KCOEF(3) + LKCOEF(K,4)KCOEF(4)
209
210	TAUGAS = U(k)*ANS
211
212	if(graybody)then
213	TAUGAS = 0.0
214	DCONT = U(k)*3.3e-26
215	endif
216
217	TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS + DCONT
218	!TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS
219
220	DTAUKI(K,nw,ng) = TAUGAS + DCONT ! For parameterized continuum absorption
221
222	do iaer=1,naerkind
223	DTAUKI(K,nw,ng) = DTAUKI(K,nw,ng) + TAEROS(K,NW,IAER)
224	end do ! a bug was here!
225
226	end do
227
228	! Now fill in the "clear" part of the spectrum (NG = L_NGAUSS),
229	! which holds continuum opacity only
230
231	NG = L_NGAUSS
232	DTAUKI(K,nw,ng) = 0.0 + DCONT ! For parameterized continuum absorption
233
234	do iaer=1,naerkind
235	DTAUKI(K,nw,ng) = DTAUKI(K,nw,ng) + TAEROS(K,NW,IAER)
236	end do ! a bug was here!
237
238	end do
239	end do
240
241
242	!=======================================================================
243	! Now the full treatment for the layers, where besides the opacity
244	! we need to calculate the scattering albedo and asymmetry factors
245
246	DO NW=1,L_NSPECTI
247	DO K=2,L_LEVELS+1
248	do iaer=1,naerkind
249	TAUAEROLK(K,NW,IAER) = TAUAERO(K,IAER)*QSIAER(K,NW,IAER)
250	end do
251	ENDDO
252	ENDDO
253
254	DO NW=1,L_NSPECTI
255	NG = L_NGAUSS
256	DO L=1,L_NLAYRAD
257
258	K = 2*L+1
259	DTAUI(L,nw,ng) = DTAUKI(K,NW,NG) + DTAUKI(K+1,NW,NG)! + 1.e-50
260
261	atemp = 0.
262	btemp = 0.
263	if(DTAUI(L,NW,NG) .GT. 1.0E-9) then
264	do iaer=1,naerkind
265	atemp = atemp + &
266	GIAER(K,NW,IAER) * TAUAEROLK(K,NW,IAER) + &
267	GIAER(K+1,NW,IAER) * TAUAEROLK(K+1,NW,IAER)
268	btemp = btemp + TAUAEROLK(K,NW,IAER) + TAUAEROLK(K+1,NW,IAER)
269	! * + 1.e-10
270	end do
271	WBARI(L,nw,ng) = btemp / DTAUI(L,NW,NG)
272	else
273	WBARI(L,nw,ng) = 0.0D0
274	DTAUI(L,NW,NG) = 1.0E-9
275	endif
276
277	if(btemp .GT. 0.0) then
278	cosbi(L,NW,NG) = atemp/btemp
279	else
280	cosbi(L,NW,NG) = 0.0D0
281	end if
282
283	END DO ! L vertical loop
284
285	! Now the other Gauss points, if needed.
286
287	DO NG=1,L_NGAUSS-1
288	IF(TAUGSURF(NW,NG) .gt. TLIMIT) THEN
289
290	DO L=1,L_NLAYRAD
291	K = 2*L+1
292	DTAUI(L,nw,ng) = DTAUKI(K,NW,NG)+DTAUKI(K+1,NW,NG)! + 1.e-50
293
294	btemp = 0.
295	if(DTAUI(L,NW,NG) .GT. 1.0E-9) then
296
297	do iaer=1,naerkind
298	btemp = btemp + TAUAEROLK(K,NW,IAER) + TAUAEROLK(K+1,NW,IAER)
299	end do
300	WBARI(L,nw,ng) = btemp / DTAUI(L,NW,NG)
301
302	else
303	WBARI(L,nw,ng) = 0.0D0
304	DTAUI(L,NW,NG) = 1.0E-9
305	endif
306
307	cosbi(L,NW,NG) = cosbi(L,NW,L_NGAUSS)
308	END DO ! L vertical loop
309	END IF
310
311	END DO ! NG Gauss loop
312	END DO ! NW spectral loop
313
314	! Total extinction optical depths
315
316	DO NW=1,L_NSPECTI
317	DO NG=1,L_NGAUSS ! full gauss loop
318	TAUI(1,NW,NG)=0.0D0
319	DO L=1,L_NLAYRAD
320	TAUI(L+1,NW,NG)=TAUI(L,NW,NG)+DTAUI(L,NW,NG)
321	END DO
322
323	TAUCUMI(1,NW,NG)=0.0D0
324	DO K=2,L_LEVELS
325	TAUCUMI(K,NW,NG)=TAUCUMI(K-1,NW,NG)+DTAUKI(K,NW,NG)
326	END DO
327	END DO ! end full gauss loop
328	END DO
329
330	return
331
332
333	end subroutine optci
334
335
336

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

Download in other formats:

Original Format