Context Navigation

optci.F90 @ 1351

Last change on this file since 1351 was 1194, checked in by sglmd, 11 years ago
Latitude-dependent gravity field added. Option oblate = .true. in callphys.def, and three additional variables needed: J2, Rmean and MassPlanet?.
File size: 13.0 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,indi,glat_ig
7	use gases_h
8	implicit none
9
10	!==================================================================
11	!
12	! Purpose
13	! -------
14	! Calculates longwave optical constants at each level. For each
15	! layer and spectral interval in the IR it calculates WBAR, DTAU
16	! and COSBAR. For each level it calculates TAU.
17	!
18	! TAUI(L,LW) is the cumulative optical depth at level L (or alternatively
19	! at the bottom of layer L), LW is the spectral wavelength interval.
20	!
21	! TLEV(L) - Temperature at the layer boundary (i.e., level)
22	! PLEV(L) - Pressure at the layer boundary (i.e., level)
23	!
24	! Authors
25	! -------
26	! Adapted from the NASA Ames code by R. Wordsworth (2009)
27	!
28	!==================================================================
29
30	#include "comcstfi.h"
31	#include "callkeys.h"
32
33
34	real*8 DTAUI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
35	real*8 DTAUKI(L_LEVELS+1,L_NSPECTI,L_NGAUSS)
36	real*8 TAUI(L_NLEVRAD,L_NSPECTI,L_NGAUSS)
37	real*8 TAUCUMI(L_LEVELS,L_NSPECTI,L_NGAUSS)
38	real*8 PLEV(L_LEVELS)
39	real*8 TLEV(L_LEVELS)
40	real*8 TMID(L_LEVELS), PMID(L_LEVELS)
41	real*8 COSBI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
42	real*8 WBARI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
43
44	! for aerosols
45	real*8 QXIAER(L_LEVELS+1,L_NSPECTI,NAERKIND)
46	real*8 QSIAER(L_LEVELS+1,L_NSPECTI,NAERKIND)
47	real*8 GIAER(L_LEVELS+1,L_NSPECTI,NAERKIND)
48	real*8 TAUAERO(L_LEVELS+1,NAERKIND)
49	real*8 TAUAEROLK(L_LEVELS+1,L_NSPECTI,NAERKIND)
50	real*8 TAEROS(L_LEVELS,L_NSPECTI,NAERKIND)
51
52	integer L, NW, NG, K, LK, IAER
53	integer MT(L_LEVELS), MP(L_LEVELS), NP(L_LEVELS)
54	real*8 ANS, TAUGAS
55	real*8 DPR(L_LEVELS), U(L_LEVELS)
56	real*8 LCOEF(4), LKCOEF(L_LEVELS,4)
57
58	real*8 taugsurf(L_NSPECTI,L_NGAUSS-1)
59	real*8 DCONT,DAERO
60	double precision wn_cont, p_cont, p_air, T_cont, dtemp, dtempc
61	double precision p_cross
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
70	real*8 btemp(L_NLAYRAD,L_NSPECTI)
71
72	! variables for k in units m^-1
73	real*8 dz(L_LEVELS)
74	!real*8 rho !! see test below
75
76	integer igas, jgas
77
78	integer interm
79
80	!--- Kasting's CIA ----------------------------------------
81	!real*8, parameter :: Ci(L_NSPECTI)=[ &
82	! 3.8E-5, 1.2E-5, 2.8E-6, 7.6E-7, 4.5E-7, 2.3E-7, &
83	! 5.4E-7, 1.6E-6, 0.0, &
84	! 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
85	! 0.0, 4.0E-7, 4.0E-6, 1.4E-5, &
86	! 1.0E-5, 1.2E-6, 2.0E-7, 5.0E-8, 3.0E-8, 0.0 ]
87	!real*8, parameter :: Ti(L_NSPECTI)=[ -2.2, -1.9, &
88	! -1.7, -1.7, -1.7, -1.7, -1.7, -1.7, &
89	! 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
90	! -1.7,-1.7,-1.7,-1.7,-1.7,-1.7,-1.7, -1.7,0.0 ]
91	!----------------------------------------------------------
92
93	!! AS: to save time in computing continuum (see bilinearbig)
94	IF (.not.ALLOCATED(indi)) THEN
95	ALLOCATE(indi(L_NSPECTI,ngasmx,ngasmx))
96	indi = -9999 ! this initial value means "to be calculated"
97	ENDIF
98
99	!=======================================================================
100	! Determine the total gas opacity throughout the column, for each
101	! spectral interval, NW, and each Gauss point, NG.
102
103	taugsurf(:,:) = 0.0
104	dpr(:) = 0.0
105	lkcoef(:,:) = 0.0
106
107	do K=2,L_LEVELS
108	DPR(k) = PLEV(K)-PLEV(K-1)
109
110	!--- Kasting's CIA ----------------------------------------
111	!dz(k)=dpr(k)189.02TMID(K)/(0.03720*PMID(K))
112	! this is CO2 path length (in cm) as written by Francois
113	! delta_z = delta_p * R_specific * T / (g * P)
114	! But Kasting states that W is in units of _atmosphere_ cm
115	! So we do
116	!dz(k)=dz(k)*(PMID(K)/1013.25)
117	!dz(k)=dz(k)/100.0 ! in m for SI calc
118	!----------------------------------------------------------
119
120	! if we have continuum opacities, we need dz
121	if(kastprof)then
122	dz(k) = dpr(k)(1000.0d08.3145d0/muvar(k))TMID(K)/(gPMID(K))
123	U(k) = CmkDPR(k)mugaz/muvar(k)
124	else
125	dz(k) = dpr(k)RTMID(K)/(glat_igPMID(K))mugaz/muvar(k)
126	U(k) = CmkDPR(k)mugaz/muvar(k) ! only Cmk line in optci.F
127	!JL13 the mugaz/muvar factor takes into account water meanmolecular weight if water is present
128	endif
129
130	call tpindex(PMID(K),TMID(K),QVAR(K),pfgasref,tgasref,WREFVAR, &
131	LCOEF,MT(K),MP(K),NVAR(K),WRATIO(K))
132
133	do LK=1,4
134	LKCOEF(K,LK) = LCOEF(LK)
135	end do
136	end do ! levels
137
138
139	do iaer=1,naerkind
140	DO NW=1,L_NSPECTI
141	do K=2,L_LEVELS
142	TAEROS(K,NW,IAER) = TAUAERO(K,IAER) * QXIAER(K,NW,IAER)
143	end do ! levels
144	END DO
145	end do
146
147	do NW=1,L_NSPECTI
148
149	do K=2,L_LEVELS
150
151	! continuum absorption
152	DCONT = 0.0d0
153
154	if(continuum.and.(.not.graybody))then
155	! include continua if necessary
156	wn_cont = dble(wnoi(nw))
157	T_cont = dble(TMID(k))
158	do igas=1,ngasmx
159
160	if(gfrac(igas).eq.-1)then ! variable
161	p_cont = dble(PMID(k)scalepQVAR(k)) ! qvar = mol/mol
162	else
163	p_cont = dble(PMID(k)scalepgfrac(igas)*(1.-QVAR(k)))
164	endif
165
166	dtemp=0.0d0
167	if(igas.eq.igas_N2)then
168
169	interm = indi(nw,igas,igas)
170	call interpolateN2N2(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
171	indi(nw,igas,igas) = interm
172
173	elseif(igas.eq.igas_H2)then
174
175	! first do self-induced absorption
176	interm = indi(nw,igas,igas)
177	call interpolateH2H2(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
178	indi(nw,igas,igas) = interm
179
180	! then cross-interactions with other gases
181	do jgas=1,ngasmx
182	p_cross = dble(PMID(k)scalepgfrac(jgas)*(1.-QVAR(k)))
183	dtempc = 0.0d0
184	if(jgas.eq.igas_N2)then
185	interm = indi(nw,igas,jgas)
186	call interpolateN2H2(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
187	indi(nw,igas,jgas) = interm
188	elseif(jgas.eq.igas_He)then
189	interm = indi(nw,igas,jgas)
190	call interpolateH2He(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
191	indi(nw,igas,jgas) = interm
192	endif
193	dtemp = dtemp + dtempc
194	enddo
195
196	elseif(igas.eq.igas_H2O.and.T_cont.gt.200.0)then
197
198	p_air = dble(PMID(k)*scalep) - p_cont ! note assumes background is air!
199	if(H2Ocont_simple)then
200	call interpolateH2Ocont_PPC(wn_cont,T_cont,p_cont,p_air,dtemp,.false.)
201	else
202	interm = indi(nw,igas,igas)
203	call interpolateH2Ocont_CKD(wn_cont,T_cont,p_cont,p_air,dtemp,.false.,interm)
204	indi(nw,igas,igas) = interm
205	endif
206
207	endif
208
209	DCONT = DCONT + dtemp
210
211	enddo
212
213	! Oobleck test
214	!rho = PMID(k)scalep / (TMID(k)286.99)
215	!if(WNOI(nw).gt.300.0 .and. WNOI(nw).lt.500.0)then
216	! DCONT = rho * 0.125 * 4.6e-4
217	!elseif(WNOI(nw).gt.500.0 .and. WNOI(nw).lt.700.0)then
218	! DCONT = 1000dpr(k) 1.0 * 4.6e-4 / g
219	! DCONT = rho * 1.0 * 4.6e-4
220	!elseif(WNOI(nw).gt.700.0 .and. WNOI(nw).lt.900.0)then
221	! DCONT = rho * 0.125 * 4.6e-4
222	!endif
223
224	DCONT = DCONT*dz(k)
225
226	endif
227
228	! aerosol absorption
229	DAERO=SUM(TAEROS(K,NW,1:naerkind))
230
231	do ng=1,L_NGAUSS-1
232
233	! Now compute TAUGAS
234
235	! Interpolate between water mixing ratios
236	! WRATIO = 0.0 if the requested water amount is equal to, or outside the
237	! the water data range
238
239	if(L_REFVAR.eq.1)then ! added by RW for special no variable case
240	KCOEF(1) = GASI(MT(K),MP(K),1,NW,NG)
241	KCOEF(2) = GASI(MT(K),MP(K)+1,1,NW,NG)
242	KCOEF(3) = GASI(MT(K)+1,MP(K)+1,1,NW,NG)
243	KCOEF(4) = GASI(MT(K)+1,MP(K),1,NW,NG)
244	else
245
246	KCOEF(1) = GASI(MT(K),MP(K),NVAR(K),NW,NG) + WRATIO(K)* &
247	(GASI(MT(K),MP(K),NVAR(K)+1,NW,NG) - &
248	GASI(MT(K),MP(K),NVAR(K),NW,NG))
249
250	KCOEF(2) = GASI(MT(K),MP(K)+1,NVAR(K),NW,NG) + WRATIO(K)* &
251	(GASI(MT(K),MP(K)+1,NVAR(K)+1,NW,NG) - &
252	GASI(MT(K),MP(K)+1,NVAR(K),NW,NG))
253
254	KCOEF(3) = GASI(MT(K)+1,MP(K)+1,NVAR(K),NW,NG) + WRATIO(K)* &
255	(GASI(MT(K)+1,MP(K)+1,NVAR(K)+1,NW,NG) - &
256	GASI(MT(K)+1,MP(K)+1,NVAR(K),NW,NG))
257
258	KCOEF(4) = GASI(MT(K)+1,MP(K),NVAR(K),NW,NG) + WRATIO(K)* &
259	(GASI(MT(K)+1,MP(K),NVAR(K)+1,NW,NG) - &
260	GASI(MT(K)+1,MP(K),NVAR(K),NW,NG))
261
262	endif
263
264	! Interpolate the gaseous k-coefficients to the requested T,P values
265
266	ANS = LKCOEF(K,1)KCOEF(1) + LKCOEF(K,2)KCOEF(2) + &
267	LKCOEF(K,3)KCOEF(3) + LKCOEF(K,4)KCOEF(4)
268
269	TAUGAS = U(k)*ANS
270
271	TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS + DCONT
272	DTAUKI(K,nw,ng) = TAUGAS &
273	+ DCONT & ! For parameterized continuum absorption
274	+ DAERO ! For aerosol absorption
275
276	end do
277
278	! Now fill in the "clear" part of the spectrum (NG = L_NGAUSS),
279	! which holds continuum opacity only
280
281	NG = L_NGAUSS
282	DTAUKI(K,nw,ng) = 0.d0 &
283	+ DCONT & ! For parameterized continuum absorption
284	+ DAERO ! For aerosol absorption
285
286	end do
287	end do
288
289	DTAUKI(L_LEVELS+1,1:L_NSPECTI,1:L_NGAUSS)=0.d0
290	!=======================================================================
291	! Now the full treatment for the layers, where besides the opacity
292	! we need to calculate the scattering albedo and asymmetry factors
293
294	do iaer=1,naerkind
295	DO NW=1,L_NSPECTI
296	DO K=2,L_LEVELS+1
297	TAUAEROLK(K,NW,IAER) = TAUAERO(K,IAER)*QSIAER(K,NW,IAER)
298	ENDDO
299	ENDDO
300	end do
301
302	DO NW=1,L_NSPECTI
303	DO L=1,L_NLAYRAD
304	K = 2*L+1
305	btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind)) + SUM(TAUAEROLK(K+1,NW,1:naerkind))
306	END DO ! L vertical loop
307	END DO ! NW spectral loop
308
309
310	DO NW=1,L_NSPECTI
311	NG = L_NGAUSS
312	DO L=1,L_NLAYRAD
313
314	K = 2*L+1
315	DTAUI(L,nw,ng) = DTAUKI(K,NW,NG) + DTAUKI(K+1,NW,NG)! + 1.e-50
316
317	atemp = 0.
318	if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
319	do iaer=1,naerkind
320	atemp = atemp + &
321	GIAER(K,NW,IAER) * TAUAEROLK(K,NW,IAER) + &
322	GIAER(K+1,NW,IAER) * TAUAEROLK(K+1,NW,IAER)
323	end do
324	WBARI(L,nw,ng) = btemp(L,nw) / DTAUI(L,NW,NG)
325	else
326	WBARI(L,nw,ng) = 0.0D0
327	DTAUI(L,NW,NG) = 1.0D-9
328	endif
329
330	if(btemp(L,nw) .GT. 0.0d0) then
331	cosbi(L,NW,NG) = atemp/btemp(L,nw)
332	else
333	cosbi(L,NW,NG) = 0.0D0
334	end if
335
336	END DO ! L vertical loop
337
338	! Now the other Gauss points, if needed.
339
340	DO NG=1,L_NGAUSS-1
341	IF(TAUGSURF(NW,NG) .gt. TLIMIT) THEN
342
343	DO L=1,L_NLAYRAD
344	K = 2*L+1
345	DTAUI(L,nw,ng) = DTAUKI(K,NW,NG)+DTAUKI(K+1,NW,NG)! + 1.e-50
346
347	if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
348
349	WBARI(L,nw,ng) = btemp(L,nw) / DTAUI(L,NW,NG)
350
351	else
352	WBARI(L,nw,ng) = 0.0D0
353	DTAUI(L,NW,NG) = 1.0D-9
354	endif
355
356	cosbi(L,NW,NG) = cosbi(L,NW,L_NGAUSS)
357	END DO ! L vertical loop
358	END IF
359
360	END DO ! NG Gauss loop
361	END DO ! NW spectral loop
362
363	! Total extinction optical depths
364
365	DO NG=1,L_NGAUSS ! full gauss loop
366	DO NW=1,L_NSPECTI
367	TAUCUMI(1,NW,NG)=0.0D0
368	DO K=2,L_LEVELS
369	TAUCUMI(K,NW,NG)=TAUCUMI(K-1,NW,NG)+DTAUKI(K,NW,NG)
370	END DO
371	END DO ! end full gauss loop
372	END DO
373
374	! be aware when comparing with textbook results
375	! (e.g. Pierrehumbert p. 218) that
376	! taucumi does not take the <cos theta>=0.5 factor into
377	! account. It is the optical depth for a vertically
378	! ascending ray with angle theta = 0.
379
380	!open(127,file='taucum.out')
381	!do nw=1,L_NSPECTI
382	! write(127,*) taucumi(L_LEVELS,nw,L_NGAUSS)
383	!enddo
384	!close(127)
385
386	! print*,'WBARI'
387	! print*,WBARI
388	! print*,'DTAUI'
389	! print*,DTAUI
390	! call abort
391
392
393	return
394
395
396	end subroutine optci
397
398
399

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

Download in other formats:

Original Format