Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

optci.F90

Last change on this file was 3497, checked in by debatzbr, 2 weeks ago
Add AC6H6 condensation in the microphysics
File size: 22.2 KB

Line
1	subroutine optci(PQMO,NLAY,ZLEV,PLEV,TLEV,TMID,PMID, &
2	DTAUI,TAUCUMI,COSBI,WBARI,TAUGSURF,SEASHAZEFACT,&
3	DIAG_OPTH,DIAG_OPTT,CDCOLUMN)
4
5	use radinc_h
6	use radcommon_h, only: gasi,gasi_recomb,tlimit,Cmk,gzlat_ig, &
7	tgasref,pfgasref,wnoi,scalep,indi
8	use gases_h
9	use datafile_mod, only: haze_opt_file
10	use comcstfi_mod, only: pi,r
11	use callkeys_mod, only: continuum,graybody,corrk_recombin, &
12	callclouds,callmufi,seashaze,uncoupl_optic_haze,&
13	opt4clouds,FHIR,FCIR
14	use tracer_h, only: nmicro,nice,ices_indx
15
16	implicit none
17
18	!==================================================================
19	!
20	! Purpose
21	! -------
22	! Calculates longwave optical constants at each level. For each
23	! layer and spectral interval in the IR it calculates WBAR, DTAU
24	! and COSBAR. For each level it calculates TAU.
25	!
26	! TAUCUMI(L,LW) is the cumulative optical depth at level L (or alternatively
27	! at the bottom of layer L), LW is the spectral wavelength interval.
28	!
29	! TLEV(L) - Temperature at the layer boundary (i.e., level)
30	! PLEV(L) - Pressure at the layer boundary (i.e., level)
31	!
32	! Authors
33	! -------
34	! Adapted from the NASA Ames code by R. Wordsworth (2009)
35	!
36	! Modified
37	! --------
38	! J. Vatant d'Ollone (2016-17)
39	! --> Clean and adaptation to Titan
40	! B. de Batz de Trenquelléon (2022-2023)
41	! --> Clean and correction
42	! --> New optics added for Titan's clouds
43	!
44	!==================================================================
45
46
47	!==========================================================
48	! Input/Output
49	!==========================================================
50	REAL*8, INTENT(IN) :: PQMO(nlay,nmicro) ! Tracers for microphysics optics (X/m2).
51	INTEGER, INTENT(IN) :: NLAY ! Number of pressure layers (for pqmo)
52	REAL*8, INTENT(IN) :: ZLEV(NLAY+1)
53	REAL*8, INTENT(IN) :: PLEV(L_LEVELS), TLEV(L_LEVELS)
54	REAL*8, INTENT(IN) :: TMID(L_LEVELS), PMID(L_LEVELS)
55	REAL*8, INTENT(IN) :: SEASHAZEFACT(L_LEVELS)
56	INTEGER, INTENT(IN) :: CDCOLUMN
57
58	REAL*8, INTENT(OUT) :: DTAUI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
59	REAL*8, INTENT(OUT) :: TAUCUMI(L_LEVELS,L_NSPECTI,L_NGAUSS)
60	REAL*8, INTENT(OUT) :: COSBI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
61	REAL*8, INTENT(OUT) :: WBARI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
62	REAL*8, INTENT(OUT) :: TAUGSURF(L_NSPECTI,L_NGAUSS-1)
63	REAL*8, INTENT(OUT) :: DIAG_OPTH(L_LEVELS,L_NSPECTI,6)
64	REAL*8, INTENT(OUT) :: DIAG_OPTT(L_LEVELS,L_NSPECTI,6)
65	! ==========================================================
66
67	real*8 DTAUKI(L_LEVELS,L_NSPECTI,L_NGAUSS)
68
69	! Titan customisation
70	! J. Vatant d'Ollone (2016)
71	real*8 DHAZE_T(L_LEVELS,L_NSPECTI)
72	real*8 DHAZES_T(L_LEVELS,L_NSPECTI)
73	real*8 SSA_T(L_LEVELS,L_NSPECTI)
74	real*8 ASF_T(L_LEVELS,L_NSPECTI)
75	! ==========================
76
77	integer L, NW, NG, K, LK, IAER
78	integer MT(L_LEVELS), MP(L_LEVELS), NP(L_LEVELS)
79	real*8 ANS, TAUGAS
80	real*8 DPR(L_LEVELS), U(L_LEVELS)
81	real*8 LCOEF(4), LKCOEF(L_LEVELS,4)
82
83	real*8 DCONT
84	double precision wn_cont, p_cont, p_air, T_cont, dtemp, dtempc
85	double precision p_cross
86
87	real*8 KCOEF(4)
88
89	! temporary variable to reduce memory access time to gasi
90	real*8 tmpk(2,2)
91
92	! temporary variables for multiple aerosol calculation
93	real*8 atemp
94	real*8 btemp(L_NLAYRAD,L_NSPECTI)
95
96	! variables for k in units m^-1
97	real*8 dz(L_LEVELS)
98
99	integer igas, jgas, ilay
100
101	integer interm
102
103	! Variables for haze optics
104	character(len=200) file_path
105	logical file_ok
106	integer dumch
107	real*8 dumwvl
108	integer ilev_cutoff
109	real*8 corr_haze
110
111	! Variables for new optics
112	integer iq, iw, FTYPE, CTYPE
113	real*8 m0as,m0af,m0ccn,m3as,m3af,m3ccn,m3cld
114	real*8 dtauaer_s,dtauaer_f,dtau_ccn,dtau_cld
115	real*8,save :: rhoaer_s(L_NSPECTI),ssa_s(L_NSPECTI),asf_s(L_NSPECTI)
116	real*8,save :: rhoaer_f(L_NSPECTI),ssa_f(L_NSPECTI),asf_f(L_NSPECTI)
117	real*8,save :: ssa_ccn(L_NSPECTI),asf_ccn(L_NSPECTI)
118	real*8,save :: ssa_cld(L_NSPECTI),asf_cld(L_NSPECTI)
119	!$OMP THREADPRIVATE(rhoaer_s,rhoaer_f,ssa_s,ssa_f,ssa_cld,asf_s,asf_f,asf_cld)
120
121	logical,save :: firstcall=.true.
122	!$OMP THREADPRIVATE(firstcall)
123
124
125	!! AS: to save time in computing continuum (see bilinearbig)
126	IF (.not.ALLOCATED(indi)) THEN
127	ALLOCATE(indi(L_NSPECTI,ngasmx,ngasmx))
128	indi = -9999 ! this initial value means "to be calculated"
129	ENDIF
130
131	! Some initialisation because there's a pb with disr_haze at the limits (nw=1)
132	! I should check this - For now we set vars to zero : better than nans - JVO 2017
133	DHAZE_T(:,:) = 0.0
134	SSA_T(:,:) = 0.0
135	ASF_T(:,:) = 0.0
136
137	! Load tabulated haze optical properties if needed.
138	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
139	IF (firstcall .AND. callmufi .AND. (.NOT. uncoupl_optic_haze)) THEN
140	OPEN(12,file=TRIM(haze_opt_file),form='formatted') ! The file has been inquired in physiq_mod firstcall
141	READ(12,*) ! dummy header
142	DO NW=1,L_NSPECTI
143	READ(12,*) dumch, dumwvl, rhoaer_f(nw), ssa_f(nw), asf_f(nw), rhoaer_s(nw), ssa_s(nw), asf_s(nw)
144	ENDDO
145	CLOSE(12)
146	ENDIF
147
148	!=======================================================================
149	! Determine the total gas opacity throughout the column, for each
150	! spectral interval, NW, and each Gauss point, NG.
151
152	taugsurf(:,:) = 0.0
153	dpr(:) = 0.0
154	lkcoef(:,:) = 0.0
155
156	! Level of cutoff
157	!~~~~~~~~~~~~~~~~
158	ilev_cutoff = 26
159
160	do K=2,L_LEVELS
161	ilay = L_NLAYRAD+1 - k/2 ! int. arithmetic => gives the gcm layer index (reversed)
162	DPR(k) = PLEV(K)-PLEV(K-1)
163
164	! if we have continuum opacities, we need dz
165	dz(k) = dpr(k)RTMID(K)/(gzlat_ig(ilay)*PMID(K))
166	U(k) = Cmk(ilay)*DPR(k) ! only Cmk line in optci.F
167
168	call tpindex(PMID(K),TMID(K),pfgasref,tgasref,LCOEF,MT(K),MP(K))
169
170	do LK=1,4
171	LKCOEF(K,LK) = LCOEF(LK)
172	end do
173	end do ! L_LEVELS
174
175	do NW=1,L_NSPECTI
176	! We ignore K=1...
177	do K=2,L_LEVELS
178	! int. arithmetic => gives the gcm layer index (reversed)
179	ilay = L_NLAYRAD+1 - k/2
180
181	! Optics coupled with the microphysics :
182	IF (callmufi .AND. (.NOT. uncoupl_optic_haze)) THEN
183
184	!==========================================================================================
185	! Old optics (must have callclouds = .false.):
186	!==========================================================================================
187	IF (.NOT. opt4clouds) THEN
188	m3as = pqmo(ilay,2) / 2.0
189	m3af = pqmo(ilay,4) / 2.0
190	! Cut-off
191	IF (ilay .lt. ilev_cutoff) THEN
192	m3as = pqmo(ilev_cutoff,2) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(ilev_cutoff+1)-zlev(ilev_cutoff))
193	m3af = pqmo(ilev_cutoff,4) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(ilev_cutoff+1)-zlev(ilev_cutoff))
194	ENDIF
195
196	dtauaer_s = m3as*rhoaer_s(nw)
197	dtauaer_f = m3af*rhoaer_f(nw)
198
199	!==========================================================================================
200	! New optics :
201	!==========================================================================================
202	ELSE
203	iw = (L_NSPECTI + 1) - NW + L_NSPECTV ! Visible first and return
204	!-----------------------------
205	! HAZE (Spherical + Fractal) :
206	!-----------------------------
207	FTYPE = 1
208
209	! Spherical aerosols :
210	!---------------------
211	CTYPE = 5
212	m0as = pqmo(ilay,1) / 2.0
213	m3as = pqmo(ilay,2) / 2.0
214	! If not callclouds : must have a cut-off
215	IF (.NOT. callclouds) THEN
216	IF (ilay .lt. ilev_cutoff) THEN
217	m0as = pqmo(ilev_cutoff,1) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(ilev_cutoff+1)-zlev(ilev_cutoff))
218	m3as = pqmo(ilev_cutoff,2) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(ilev_cutoff+1)-zlev(ilev_cutoff))
219	ENDIF
220	ENDIF
221	call get_haze_and_cloud_opacity(FTYPE,CTYPE,m0as,m3as,iw,dtauaer_s,ssa_s(nw),asf_s(nw))
222
223	! Fractal aerosols :
224	!-------------------
225	CTYPE = FTYPE
226	m0af = pqmo(ilay,3) / 2.0
227	m3af = pqmo(ilay,4) / 2.0
228	! If not callclouds : must have a cut-off
229	IF (.NOT. callclouds) THEN
230	IF (ilay .lt. ilev_cutoff) THEN
231	m0af = pqmo(ilev_cutoff,3) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(ilev_cutoff+1)-zlev(ilev_cutoff))
232	m3af = pqmo(ilev_cutoff,4) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(ilev_cutoff+1)-zlev(ilev_cutoff))
233	ENDIF
234	ENDIF
235	call get_haze_and_cloud_opacity(FTYPE,CTYPE,m0af,m3af,iw,dtauaer_f,ssa_f(nw),asf_f(nw))
236	ENDIF
237
238	! Tuning of optical properties for haze :
239	!dtauaer_s = dtauaer_s * (FHIR * (1-ssa_s(nw)) + ssa_s(nw))
240	!ssa_s(nw) = ssa_s(nw) / (FHIR * (1-ssa_s(nw)) + ssa_s(nw))
241	dtauaer_f = dtauaer_f * (FHIR * (1-ssa_f(nw)) + ssa_f(nw))
242	ssa_f(nw) = ssa_f(nw) / (FHIR * (1-ssa_f(nw)) + ssa_f(nw))
243
244	! Total of Haze opacity (dtau), SSA (w) and ASF (COS) :
245	DHAZE_T(k,nw) = dtauaer_s + dtauaer_f
246	IF (dtauaer_s + dtauaer_f .GT. 1.D-30) THEN
247	SSA_T(k,nw) = ( dtauaer_sssa_s(nw) + dtauaer_fssa_f(nw) ) / ( dtauaer_s+dtauaer_f )
248	ASF_T(k,nw) = ( dtauaer_sssa_s(nw)asf_s(nw) + dtauaer_fssa_f(nw)asf_f(nw) ) &
249	/ ( ssa_s(nw)dtauaer_s + ssa_f(nw)dtauaer_f )
250	ELSE
251	DHAZE_T(k,nw) = 0.D0
252	SSA_T(k,nw) = 1.0
253	ASF_T(k,nw) = 1.0
254	ENDIF
255
256	! Opacity and albedo adjustement below cutoff :
257	IF (.NOT. callclouds) THEN
258	corr_haze=0.6-0.4TANH((PMID(K)100.-2500.)/250.)
259	IF (ilay .lt. ilev_cutoff) THEN
260	DHAZE_T(k,nw) = DHAZE_T(k,nw) * corr_haze
261	ENDIF
262	ENDIF
263
264	! Diagnostics for the haze :
265	DIAG_OPTH(k,nw,1) = DHAZE_T(k,nw) ! dtau
266	DIAG_OPTH(k,nw,2) = SSA_T(k,nw) ! wbar
267	DIAG_OPTH(k,nw,3) = ASF_T(k,nw) ! gbar
268
269	!---------------------
270	! CLOUDS (Spherical) :
271	!---------------------
272	IF (callclouds) THEN
273	CTYPE = 0
274	m0ccn = pqmo(ilay,5) / 2.0
275	m3ccn = pqmo(ilay,6) / 2.0
276	m3cld = pqmo(ilay,6) / 2.0
277
278	! Clear / Dark column method :
279	!-----------------------------
280
281	! CCN's SSA :
282	call get_haze_and_cloud_opacity(FTYPE,FTYPE,m0ccn,m3ccn,iw,dtau_ccn,ssa_ccn(nw),asf_ccn(nw))
283
284	! Clear column (CCN, C2H2, C2H6, HCN, AC6H6) :
285	IF (CDCOLUMN == 0) THEN
286	DO iq = 2, nice
287	m3cld = m3cld + (pqmo(ilay,ices_indx(iq)) / 2.0)
288	ENDDO
289	call get_haze_and_cloud_opacity(FTYPE,CTYPE,m0ccn,m3cld,iw,dtau_cld,ssa_cld(nw),asf_cld(nw))
290
291	! Dark column (CCN, CH4, C2H2, C2H6, HCN, AC6H6) :
292	ELSEIF (CDCOLUMN == 1) THEN
293	DO iq = 1, nice
294	m3cld = m3cld + (pqmo(ilay,ices_indx(iq)) / 2.0)
295	ENDDO
296	call get_haze_and_cloud_opacity(FTYPE,CTYPE,m0ccn,m3cld,iw,dtau_cld,ssa_cld(nw),asf_cld(nw))
297
298	ELSE
299	WRITE(,) 'WARNING OPTCI.F90 : CDCOLUMN MUST BE 0 OR 1'
300	WRITE(,) 'WE USE DARK COLUMN ...'
301	DO iq = 1, nice
302	m3cld = m3cld + (pqmo(ilay,ices_indx(iq)) / 2.0)
303	ENDDO
304	call get_haze_and_cloud_opacity(FTYPE,CTYPE,m0ccn,m3cld,iw,dtau_cld,ssa_cld(nw),asf_cld(nw))
305	ENDIF
306
307	! For small dropplets, opacity of nucleus dominates...
308	dtau_cld = (dtau_cldm3ccn + dtau_cldm3cld) / (m3ccn + m3cld)
309	ssa_cld(nw) = (ssa_ccn(nw)m3ccn + ssa_cld(nw)m3cld) / (m3ccn + m3cld)
310
311	! Tuning of optical properties for clouds :
312	dtau_cld = dtau_cld * (FCIR * (1-ssa_cld(nw)) + ssa_cld(nw))
313	ssa_cld(nw) = ssa_cld(nw) / (FCIR * (1-ssa_cld(nw)) + ssa_cld(nw))
314
315	! Total of Haze + Clouds opacity (dtau), SSA (w) and ASF (COS) :
316	DHAZE_T(k,nw) = dtauaer_s + dtauaer_f + dtau_cld
317	IF (DHAZE_T(k,nw) .GT. 1.D-30) THEN
318	SSA_T(k,nw) = ( dtauaer_sssa_s(nw) + dtauaer_fssa_f(nw) + dtau_cld*ssa_cld(nw) ) / ( dtauaer_s+dtauaer_f+dtau_cld )
319	ASF_T(k,nw) = ( dtauaer_sssa_s(nw)asf_s(nw) + dtauaer_fssa_f(nw)asf_f(nw) + dtau_cldssa_cld(nw)asf_cld(nw) ) &
320	/ ( ssa_s(nw)dtauaer_s + ssa_f(nw)dtauaer_f + ssa_cld(nw)*dtau_cld )
321	ELSE
322	DHAZE_T(k,nw) = 0.D0
323	SSA_T(k,nw) = 1.0
324	ASF_T(k,nw) = 1.0
325	ENDIF
326
327	! Diagnostics for clouds :
328	DIAG_OPTT(k,nw,1) = DHAZE_T(k,nw) ! dtau
329	DIAG_OPTT(k,nw,2) = SSA_T(k,nw) ! wbar
330	DIAG_OPTT(k,nw,3) = ASF_T(k,nw) ! gbar
331
332	ELSE
333	! Diagnostics for clouds :
334	DIAG_OPTT(k,nw,1) = 0.D0 ! dtau
335	DIAG_OPTT(k,nw,2) = 1.0 ! wbar
336	DIAG_OPTT(k,nw,3) = 1.0 ! gbar
337	ENDIF
338
339	! Optics and microphysics no coupled :
340	ELSE
341	! Call fixed vertical haze profile of extinction - same for all columns
342	call disr_haze(dz(k),plev(k),wnoi(nw),DHAZE_T(k,nw),SSA_T(k,nw),ASF_T(k,nw))
343	if (seashaze) DHAZE_T(k,nw) = DHAZE_T(k,nw)*seashazefact(k)
344	! Diagnostics for the haze :
345	DIAG_OPTH(k,nw,1) = DHAZE_T(k,nw) ! dtau
346	DIAG_OPTH(k,nw,2) = SSA_T(k,nw) ! wbar
347	DIAG_OPTH(k,nw,3) = ASF_T(k,nw) ! gbar
348	! Diagnostics for clouds :
349	DIAG_OPTT(k,nw,1) = 0.D0 ! dtau
350	DIAG_OPTT(k,nw,2) = 1.0 ! wbar
351	DIAG_OPTT(k,nw,3) = 1.0 ! gbar
352	ENDIF ! ENDIF callmufi
353
354	DCONT = 0.0d0 ! continuum absorption
355
356	if(continuum.and.(.not.graybody))then
357	! include continua if necessary
358	wn_cont = dble(wnoi(nw))
359	T_cont = dble(TMID(k))
360	do igas=1,ngasmx
361
362	p_cont = dble(PMID(k)scalepgfrac(igas,ilay))
363
364	dtemp=0.0d0
365	if(igas.eq.igas_N2)then
366
367	interm = indi(nw,igas,igas)
368	call interpolateN2N2(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
369	indi(nw,igas,igas) = interm
370
371	elseif(igas.eq.igas_H2)then
372
373	! first do self-induced absorption
374	interm = indi(nw,igas,igas)
375	call interpolateH2H2(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
376	indi(nw,igas,igas) = interm
377
378	! then cross-interactions with other gases
379	do jgas=1,ngasmx
380	p_cross = dble(PMID(k)scalepgfrac(jgas,ilay))
381	dtempc = 0.0d0
382	if(jgas.eq.igas_N2)then
383	interm = indi(nw,igas,jgas)
384	call interpolateN2H2(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
385	indi(nw,igas,jgas) = interm
386	endif
387	dtemp = dtemp + dtempc
388	enddo
389
390	elseif(igas.eq.igas_CH4)then
391
392	! first do self-induced absorption
393	interm = indi(nw,igas,igas)
394	call interpolateCH4CH4(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
395	indi(nw,igas,igas) = interm
396
397	! then cross-interactions with other gases
398	do jgas=1,ngasmx
399	p_cross = dble(PMID(k)scalepgfrac(jgas,ilay))
400	dtempc = 0.0d0
401	if(jgas.eq.igas_N2)then
402	interm = indi(nw,igas,jgas)
403	call interpolateN2CH4(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
404	indi(nw,igas,jgas) = interm
405	endif
406	dtemp = dtemp + dtempc
407	enddo
408
409	endif
410
411	DCONT = DCONT + dtemp
412
413	enddo
414
415	DCONT = DCONT*dz(k)
416
417	endif
418
419	do ng=1,L_NGAUSS-1
420
421	! Now compute TAUGAS
422
423	! JVO 2017 : added tmpk because the repeated calls to gasi/v increased dramatically
424	! the execution time of optci/v -> ~ factor 2 on the whole radiative
425	! transfer on the tested simulations !
426
427	if (corrk_recombin) then
428	tmpk = GASI_RECOMB(MT(K):MT(K)+1,MP(K):MP(K)+1,NW,NG)
429	else
430	tmpk = GASI(MT(K):MT(K)+1,MP(K):MP(K)+1,1,NW,NG)
431	endif
432
433	KCOEF(1) = tmpk(1,1) ! KCOEF(1) = GASI(MT(K),MP(K),1,NW,NG)
434	KCOEF(2) = tmpk(1,2) ! KCOEF(2) = GASI(MT(K),MP(K)+1,1,NW,NG)
435	KCOEF(3) = tmpk(2,2) ! KCOEF(3) = GASI(MT(K)+1,MP(K)+1,1,NW,NG)
436	KCOEF(4) = tmpk(2,1) ! KCOEF(4) = GASI(MT(K)+1,MP(K),1,NW,NG)
437
438
439	! Interpolate the gaseous k-coefficients to the requested T,P values
440
441	ANS = LKCOEF(K,1)KCOEF(1) + LKCOEF(K,2)KCOEF(2) + &
442	LKCOEF(K,3)KCOEF(3) + LKCOEF(K,4)KCOEF(4)
443
444	TAUGAS = U(k)*ANS
445
446	TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS + DCONT
447	DTAUKI(K,nw,ng) = TAUGAS &
448	+ DCONT & ! For parameterized continuum absorption
449	+ DHAZE_T(K,NW) ! For Titan haze
450
451	end do
452
453	! Now fill in the "clear" part of the spectrum (NG = L_NGAUSS),
454	! which holds continuum opacity only
455
456	NG = L_NGAUSS
457	DTAUKI(K,nw,ng) = 0.d0 &
458	+ DCONT & ! For parameterized continuum absorption
459	+ DHAZE_T(K,NW) ! For Titan Haze
460
461	DIAG_OPTH(K,nw,4) = 0.d0
462	DIAG_OPTH(K,nw,5) = TAUGAS
463	DIAG_OPTH(K,nw,6) = DCONT
464	DIAG_OPTT(K,nw,4) = 0.d0
465	DIAG_OPTT(K,nw,5) = TAUGAS
466	DIAG_OPTT(K,nw,6) = DCONT
467
468	end do ! k = L_LEVELS
469	end do ! nw = L_NSPECTI
470
471	!=======================================================================
472	! Now the full treatment for the layers, where besides the opacity
473	! we need to calculate the scattering albedo and asymmetry factors
474	! ======================================================================
475
476	! Haze scattering
477	DO NW=1,L_NSPECTI
478	DO K=2,L_LEVELS
479	DHAZES_T(K,NW) = DHAZE_T(K,NW) * SSA_T(K,NW)
480	ENDDO
481	ENDDO
482
483	DO NW=1,L_NSPECTI
484	DO L=1,L_NLAYRAD-1
485	K = 2*L+1
486	btemp(L,NW) = DHAZES_T(K,NW) + DHAZES_T(K+1,NW)
487	END DO ! L vertical loop
488
489	! Last level
490	L = L_NLAYRAD
491	K = 2*L+1
492	btemp(L,NW) = DHAZES_T(K,NW)
493
494	END DO ! NW spectral loop
495
496
497	DO NW=1,L_NSPECTI
498	NG = L_NGAUSS
499	DO L=1,L_NLAYRAD-1
500
501	K = 2*L+1
502	DTAUI(L,nw,ng) = DTAUKI(K,NW,NG) + DTAUKI(K+1,NW,NG)! + 1.e-50
503
504	atemp = 0.
505	if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
506	atemp = atemp + &
507	ASF_T(K,NW)*DHAZES_T(K,NW) + &
508	ASF_T(K+1,NW)*DHAZES_T(K+1,NW)
509
510	WBARI(L,nw,ng) = btemp(L,nw) / DTAUI(L,NW,NG)
511	else
512	WBARI(L,nw,ng) = 0.0D0
513	DTAUI(L,NW,NG) = 1.0D-9
514	endif
515
516	if(btemp(L,nw) .GT. 0.0d0) then
517	cosbi(L,NW,NG) = atemp/btemp(L,nw)
518	else
519	cosbi(L,NW,NG) = 0.0D0
520	end if
521
522	END DO ! L vertical loop
523
524	! Last level
525
526	L = L_NLAYRAD
527	K = 2*L+1
528	DTAUI(L,nw,ng) = DTAUKI(K,NW,NG) ! + 1.e-50
529
530	atemp = 0.
531	if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
532	atemp = atemp + ASF_T(K,NW)*DHAZES_T(K,NW)
533	WBARI(L,nw,ng) = btemp(L,nw) / DTAUI(L,NW,NG)
534	else
535	WBARI(L,nw,ng) = 0.0D0
536	DTAUI(L,NW,NG) = 1.0D-9
537	endif
538
539	if(btemp(L,nw) .GT. 0.0d0) then
540	cosbi(L,NW,NG) = atemp/btemp(L,nw)
541	else
542	cosbi(L,NW,NG) = 0.0D0
543	end if
544
545
546	! Now the other Gauss points, if needed.
547
548	DO NG=1,L_NGAUSS-1
549	IF(TAUGSURF(NW,NG) .gt. TLIMIT) THEN
550
551	DO L=1,L_NLAYRAD-1
552	K = 2*L+1
553	DTAUI(L,nw,ng) = DTAUKI(K,NW,NG)+DTAUKI(K+1,NW,NG)! + 1.e-50
554
555	if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
556
557	WBARI(L,nw,ng) = btemp(L,nw) / DTAUI(L,NW,NG)
558
559	else
560	WBARI(L,nw,ng) = 0.0D0
561	DTAUI(L,NW,NG) = 1.0D-9
562	endif
563
564	cosbi(L,NW,NG) = cosbi(L,NW,L_NGAUSS)
565	END DO ! L vertical loop
566
567	! Last level
568	L = L_NLAYRAD
569	K = 2*L+1
570	DTAUI(L,nw,ng) = DTAUKI(K,NW,NG)! + 1.e-50
571
572	if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
573
574	WBARI(L,nw,ng) = btemp(L,nw) / DTAUI(L,NW,NG)
575
576	else
577	WBARI(L,nw,ng) = 0.0D0
578	DTAUI(L,NW,NG) = 1.0D-9
579	endif
580
581	cosbi(L,NW,NG) = cosbi(L,NW,L_NGAUSS)
582
583	END IF
584
585	END DO ! NG Gauss loop
586	END DO ! NW spectral loop
587
588	! Total extinction optical depths
589	!DO NG=1,L_NGAUSS ! full gauss loop
590	! DO NW=1,L_NSPECTI
591	! TAUCUMI(1,NW,NG)=0.0D0
592	! DO K=2,L_LEVELS
593	! TAUCUMI(K,NW,NG)=TAUCUMI(K-1,NW,NG)+DTAUKI(K,NW,NG)
594	! END DO
595	! END DO ! end full gauss loop
596	!END DO
597
598	TAUCUMI(:,:,:) = DTAUKI(:,:,:)
599
600	! be aware when comparing with textbook results
601	! (e.g. Pierrehumbert p. 218) that
602	! taucumi does not take the <cos theta>=0.5 factor into
603	! account. It is the optical depth for a vertically
604	! ascending ray with angle theta = 0.
605
606	if(firstcall) firstcall = .false.
607
608	return
609
610
611	end subroutine optci
612
613
614

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

Download in other formats:

Original Format