Context Navigation

optci.F90 @ 1715

Last change on this file since 1715 was 1648, checked in by jvatant, 8 years ago

Modifications to custom radiative transfer to Titan
+ Enables an altitude dependant gfrac for CIA computations

-> many radical changes in su_gases and co ..
-> read vertical CH4 profile with call_profilgases
-> Now you need a 'profile.def' that I will add in the deftank

+ Added interpolate CIA routines for CH4
+ Added temporary mean aerosol profile opacity routine (disr_haze)

File size: 13.2 KB

Line
1	subroutine optci(PLEV,TLEV,DTAUI,TAUCUMI, &
2	QXIAER,QSIAER,GIAER,COSBI,WBARI,TAUAERO, &
3	TMID,PMID,TAUGSURF,GWEIGHT)
4
5	use radinc_h
6	use radcommon_h, only: gasi,tlimit,Cmk,tgasref,pfgasref,wnoi,scalep,indi,glat_ig
7	use gases_h
8	use comcstfi_mod, only: g, r
9	use callkeys_mod, only: continuum,graybody
10	implicit none
11
12	!==================================================================
13	!
14	! Purpose
15	! -------
16	! Calculates longwave optical constants at each level. For each
17	! layer and spectral interval in the IR it calculates WBAR, DTAU
18	! and COSBAR. For each level it calculates TAU.
19	!
20	! TAUI(L,LW) is the cumulative optical depth at level L (or alternatively
21	! at the bottom of layer L), LW is the spectral wavelength interval.
22	!
23	! TLEV(L) - Temperature at the layer boundary (i.e., level)
24	! PLEV(L) - Pressure at the layer boundary (i.e., level)
25	!
26	! Authors
27	! -------
28	! Adapted from the NASA Ames code by R. Wordsworth (2009)
29	!
30	!==================================================================
31
32
33	real*8 DTAUI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
34	real*8 DTAUKI(L_LEVELS+1,L_NSPECTI,L_NGAUSS)
35	real*8 TAUI(L_NLEVRAD,L_NSPECTI,L_NGAUSS)
36	real*8 TAUCUMI(L_LEVELS,L_NSPECTI,L_NGAUSS)
37	real*8 PLEV(L_LEVELS)
38	real*8 TLEV(L_LEVELS)
39	real*8 TMID(L_LEVELS), PMID(L_LEVELS)
40	real*8 COSBI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
41	real*8 WBARI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
42
43	! for aerosols
44	real*8 QXIAER(L_LEVELS+1,L_NSPECTI,NAERKIND)
45	real*8 QSIAER(L_LEVELS+1,L_NSPECTI,NAERKIND)
46	real*8 GIAER(L_LEVELS+1,L_NSPECTI,NAERKIND)
47	real*8 TAUAERO(L_LEVELS+1,NAERKIND)
48	real*8 TAUAEROLK(L_LEVELS+1,L_NSPECTI,NAERKIND)
49	real*8 TAEROS(L_LEVELS,L_NSPECTI,NAERKIND)
50
51	! Titan customisation
52	! J. Vatant d'Ollone (2016)
53	real*8 GWEIGHT(L_NGAUSS)
54	real*8 DHAZE_T(L_LEVELS+1,L_NSPECTI)
55	real*8 DHAZES_T(L_LEVELS+1,L_NSPECTI)
56	real*8 SSA_T(L_LEVELS+1,L_NSPECTI)
57	real*8 ASF_T(L_LEVELS+1,L_NSPECTI)
58	real*8 INT_DTAU(L_NLAYRAD,L_NSPECTI)
59	real*8 K_HAZE(L_NLAYRAD,L_NSPECTI)
60
61	CHARACTER*2 str2
62	! ==========================
63
64	integer L, NW, NG, K, LK, IAER
65	integer MT(L_LEVELS), MP(L_LEVELS), NP(L_LEVELS)
66	real*8 ANS, TAUGAS
67	real*8 DPR(L_LEVELS), U(L_LEVELS)
68	real*8 LCOEF(4), LKCOEF(L_LEVELS,4)
69
70	real*8 taugsurf(L_NSPECTI,L_NGAUSS-1)
71	real*8 DCONT,DAERO
72	double precision wn_cont, p_cont, p_air, T_cont, dtemp, dtempc
73	double precision p_cross
74
75	real*8 KCOEF(4)
76
77	! temporary variables for multiple aerosol calculation
78	real*8 atemp
79	real*8 btemp(L_NLAYRAD,L_NSPECTI)
80
81	! variables for k in units m^-1
82	real*8 dz(L_LEVELS)
83	!real*8 rho !! see test below
84
85	integer igas, jgas, ilay
86
87	integer interm
88
89	!! AS: to save time in computing continuum (see bilinearbig)
90	IF (.not.ALLOCATED(indi)) THEN
91	ALLOCATE(indi(L_NSPECTI,ngasmx,ngasmx))
92	indi = -9999 ! this initial value means "to be calculated"
93	ENDIF
94
95	!=======================================================================
96	! Determine the total gas opacity throughout the column, for each
97	! spectral interval, NW, and each Gauss point, NG.
98
99	taugsurf(:,:) = 0.0
100	dpr(:) = 0.0
101	lkcoef(:,:) = 0.0
102
103	do K=2,L_LEVELS
104	DPR(k) = PLEV(K)-PLEV(K-1)
105
106	! if we have continuum opacities, we need dz
107	dz(k) = dpr(k)RTMID(K)/(glat_ig*PMID(K))
108	U(k) = Cmk*DPR(k) ! only Cmk line in optci.F
109
110	call tpindex(PMID(K),TMID(K),pfgasref,tgasref,LCOEF,MT(K),MP(K))
111
112	do LK=1,4
113	LKCOEF(K,LK) = LCOEF(LK)
114	end do
115	end do ! levels
116
117
118	do iaer=1,naerkind
119	DO NW=1,L_NSPECTI
120	do K=2,L_LEVELS
121	TAEROS(K,NW,IAER) = TAUAERO(K,IAER) * QXIAER(K,NW,IAER)
122	end do ! levels
123	END DO
124	end do
125
126	do NW=1,L_NSPECTI
127
128	do K=2,L_LEVELS
129
130	ilay = k / 2 ! int. arithmetic => gives the gcm layer index
131
132	! continuum absorption
133	DCONT = 0.0d0
134
135	if(continuum.and.(.not.graybody))then
136	! include continua if necessary
137	wn_cont = dble(wnoi(nw))
138	T_cont = dble(TMID(k))
139	do igas=1,ngasmx
140
141	p_cont = dble(PMID(k)scalepgfrac(igas,ilay))
142
143	dtemp=0.0d0
144	if(igas.eq.igas_N2)then
145
146	interm = indi(nw,igas,igas)
147	call interpolateN2N2(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
148	indi(nw,igas,igas) = interm
149
150	elseif(igas.eq.igas_H2)then
151
152	! first do self-induced absorption
153	interm = indi(nw,igas,igas)
154	call interpolateH2H2(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
155	indi(nw,igas,igas) = interm
156
157	! then cross-interactions with other gases
158	do jgas=1,ngasmx
159	p_cross = dble(PMID(k)scalepgfrac(jgas,ilay))
160	dtempc = 0.0d0
161	if(jgas.eq.igas_N2)then
162	interm = indi(nw,igas,jgas)
163	call interpolateN2H2(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
164	indi(nw,igas,jgas) = interm
165	endif
166	dtemp = dtemp + dtempc
167	enddo
168
169	elseif(igas.eq.igas_CH4)then
170
171	! first do self-induced absorption
172	interm = indi(nw,igas,igas)
173	call interpolateCH4CH4(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
174	indi(nw,igas,igas) = interm
175
176	! then cross-interactions with other gases
177	do jgas=1,ngasmx
178	p_cross = dble(PMID(k)scalepgfrac(jgas,ilay))
179	dtempc = 0.0d0
180	if(jgas.eq.igas_N2)then
181	interm = indi(nw,igas,jgas)
182	call interpolateN2CH4(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
183	indi(nw,igas,jgas) = interm
184	endif
185	dtemp = dtemp + dtempc
186	enddo
187
188	endif
189
190	DCONT = DCONT + dtemp
191
192	enddo
193
194	! Oobleck test
195	!rho = PMID(k)scalep / (TMID(k)286.99)
196	!if(WNOI(nw).gt.300.0 .and. WNOI(nw).lt.500.0)then
197	! DCONT = rho * 0.125 * 4.6e-4
198	!elseif(WNOI(nw).gt.500.0 .and. WNOI(nw).lt.700.0)then
199	! DCONT = 1000dpr(k) 1.0 * 4.6e-4 / g
200	! DCONT = rho * 1.0 * 4.6e-4
201	!elseif(WNOI(nw).gt.700.0 .and. WNOI(nw).lt.900.0)then
202	! DCONT = rho * 0.125 * 4.6e-4
203	!endif
204
205	DCONT = DCONT*dz(k)
206
207	endif
208
209	! aerosol absorption
210	DAERO=SUM(TAEROS(K,NW,1:naerkind))
211
212	!================= Titan customisation ========================================
213	call disr_haze(dz(k),plev(k),wnoi(nw),dhaze_T(k,nw),SSA_T(k,nw),ASF_T(k,nw))
214	! =============================================================================
215
216
217	do ng=1,L_NGAUSS-1
218
219	! Now compute TAUGAS
220
221	KCOEF(1) = GASI(MT(K),MP(K),1,NW,NG)
222	KCOEF(2) = GASI(MT(K),MP(K)+1,1,NW,NG)
223	KCOEF(3) = GASI(MT(K)+1,MP(K)+1,1,NW,NG)
224	KCOEF(4) = GASI(MT(K)+1,MP(K),1,NW,NG)
225
226	! Interpolate the gaseous k-coefficients to the requested T,P values
227
228	ANS = LKCOEF(K,1)KCOEF(1) + LKCOEF(K,2)KCOEF(2) + &
229	LKCOEF(K,3)KCOEF(3) + LKCOEF(K,4)KCOEF(4)
230
231	TAUGAS = U(k)*ANS
232
233	TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS + DCONT
234	DTAUKI(K,nw,ng) = TAUGAS &
235	+ DCONT & ! For parameterized continuum absorption
236	+ DAERO & ! For aerosol absorption
237	+ DHAZE_T(K,NW) ! For Titan haze
238
239	end do
240
241	! Now fill in the "clear" part of the spectrum (NG = L_NGAUSS),
242	! which holds continuum opacity only
243
244	NG = L_NGAUSS
245	DTAUKI(K,nw,ng) = 0.d0 &
246	+ DCONT & ! For parameterized continuum absorption
247	+ DAERO & ! For aerosol absorption
248	+ DHAZE_T(K,NW) ! For Titan Haze
249
250	end do
251	end do
252
253	DTAUKI(L_LEVELS+1,1:L_NSPECTI,1:L_NGAUSS)=0.d0
254
255
256	!=======================================================================
257	! Now the full treatment for the layers, where besides the opacity
258	! we need to calculate the scattering albedo and asymmetry factors
259	! ======================================================================
260
261	do iaer=1,naerkind
262	DO NW=1,L_NSPECTI
263	DO K=2,L_LEVELS+1
264	TAUAEROLK(K,NW,IAER) = TAUAERO(K,IAER)*QSIAER(K,NW,IAER)
265	ENDDO
266	ENDDO
267	end do
268
269	! Haze scattering
270	DO NW=1,L_NSPECTI
271	DO K=2,L_LEVELS+1
272	DHAZES_T(K,NW) = DHAZE_T(K,NW) * SSA_T(K,NW)
273	ENDDO
274	ENDDO
275
276	DO NW=1,L_NSPECTI
277	DO L=1,L_NLAYRAD-1
278	K = 2*L+1
279	btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind)) + SUM(TAUAEROLK(K+1,NW,1:naerkind)) &
280	+ DHAZES_T(K,NW) + DHAZES_T(K+1,NW)
281	END DO ! L vertical loop
282
283	!last level
284	L = L_NLAYRAD
285	K = 2*L+1
286
287	btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind)) + DHAZES_T(K,NW)
288
289	END DO ! NW spectral loop
290
291	! ======================================================================
292
293	DO NW=1,L_NSPECTI
294	NG = L_NGAUSS
295	DO L=1,L_NLAYRAD-1
296
297	K = 2*L+1
298
299	DTAUI(L,nw,ng) = DTAUKI(K,NW,NG) + DTAUKI(K+1,NW,NG)! + 1.e-50
300
301	atemp = 0.
302	if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
303	do iaer=1,naerkind
304	atemp = atemp + &
305	GIAER(K,NW,IAER) * TAUAEROLK(K,NW,IAER) + &
306	GIAER(K+1,NW,IAER) * TAUAEROLK(K+1,NW,IAER)
307	end do
308	atemp = atemp + ASF_T(K,NW)*DHAZES_T(K,NW)
309	WBARI(L,nw,ng) = btemp(L,nw) / DTAUI(L,NW,NG)
310	else
311	WBARI(L,nw,ng) = 0.0D0
312	DTAUI(L,NW,NG) = 1.0D-9
313	endif
314
315	if(btemp(L,nw) .GT. 0.0d0) then
316	cosbi(L,NW,NG) = atemp/btemp(L,nw)
317	else
318	cosbi(L,NW,NG) = 0.0D0
319	end if
320
321	END DO ! L vertical loop
322
323	! No vertical averaging on bottom layer
324
325	L = L_NLAYRAD
326	K = 2*L+1
327	DTAUI(L,nw,ng) = DTAUKI(K,NW,NG)
328
329	atemp = 0.
330	if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
331	do iaer=1,naerkind
332	atemp = atemp + &
333	GIAER(K,NW,IAER) * TAUAEROLK(K,NW,IAER)
334	end do
335	atemp = atemp + ASF_T(K,NW)*DHAZES_T(K,NW)
336	WBARI(L,nw,ng) = btemp(L,nw) / DTAUI(L,NW,NG)
337	else
338	WBARI(L,nw,ng) = 0.0D0
339	DTAUI(L,NW,NG) = 1.0D-9
340	endif
341
342	if(btemp(L,nw) .GT. 0.0d0) then
343	cosbi(L,NW,NG) = atemp/btemp(L,nw)
344	else
345	cosbi(L,NW,NG) = 0.0D0
346	end if
347
348	! Now the other Gauss points, if needed.
349
350	DO NG=1,L_NGAUSS-1
351
352	IF(TAUGSURF(NW,NG) .gt. TLIMIT) THEN
353
354	DO L=1,L_NLAYRAD-1
355	K = 2*L+1
356	DTAUI(L,nw,ng) = DTAUKI(K,NW,NG)+DTAUKI(K+1,NW,NG)! + 1.e-50
357
358	if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
359	WBARI(L,nw,ng) = btemp(L,nw) / DTAUI(L,NW,NG)
360	else
361	WBARI(L,nw,ng) = 0.0D0
362	DTAUI(L,NW,NG) = 1.0D-9
363	endif
364
365	cosbi(L,NW,NG) = cosbi(L,NW,L_NGAUSS)
366	END DO ! L vertical loop
367
368	! No vertical averaging on bottom layer
369
370	L = L_NLAYRAD
371	K = 2*L+1
372	DTAUI(L,nw,ng) = DTAUKI(K,NW,NG)
373	if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
374	WBARI(L,nw,ng) = btemp(L,nw) / DTAUI(L,NW,NG)
375	else
376	WBARI(L,nw,ng) = 0.0D0
377	DTAUI(L,NW,NG) = 1.0D-9
378	endif
379	cosbi(L,NW,NG) = cosbi(L,NW,L_NGAUSS)
380
381	END IF
382
383	END DO ! NG Gauss loop
384	END DO ! NW spectral loop
385
386	! Total extinction optical depths
387
388	DO NG=1,L_NGAUSS ! full gauss loop
389	DO NW=1,L_NSPECTI
390	TAUCUMI(1,NW,NG)=0.0D0
391	DO K=2,L_LEVELS
392	TAUCUMI(K,NW,NG)=TAUCUMI(K-1,NW,NG)+DTAUKI(K,NW,NG)
393	END DO
394	END DO ! end full gauss loop
395	END DO
396
397	! be aware when comparing with textbook results
398	! (e.g. Pierrehumbert p. 218) that
399	! taucumi does not take the <cos theta>=0.5 factor into
400	! account. It is the optical depth for a vertically
401	! ascending ray with angle theta = 0.
402
403
404	! Titan's outputs (J.V.O, 2016)===============================================
405	! do l=1,L_NLAYRAD
406	! do nw=1,L_NSPECTI
407	! INT_DTAU(L,NW) = 0.0d+0
408	! DO NG=1,L_NGAUSS
409	! INT_DTAU(L,NW)= INT_DTAU(L,NW) + dtaui(L,nw,ng)*gweight(NG)
410	! enddo
411	! enddo
412	! enddo
413
414	! do nw=1,L_NSPECTI
415	! write(str2,'(i2.2)') nw
416	! call writediagfi(1,'kgi'//str2,'Gaz extinction coefficient IR band '//str2,'m-1',1,int_dtau(L_NLAYRAD:1:-1,nw)/dz_lay(L_NLAYRAD:1:-1))
417	! call writediagfi(1,'khi'//str2,'Haze extinction coefficient IR band '//str2,'m-1',1,k_haze(L_NLAYRAD:1:-1,nw)/dz_lay(L_NLAYRAD:1:-1))
418	! enddo
419
420	! ==============================================================================
421
422	return
423
424
425	end subroutine optci
426
427
428

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

Download in other formats:

Original Format