Context Navigation

optci.F90 @ 1781

Last change on this file since 1781 was 1781, checked in by jvatant, 7 years ago
Useless argument Gweight in rad. tr. routines ( present in radcommon.h ) -JVO
File size: 13.5 KB

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

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

Download in other formats:

Original Format