Context Navigation

optci.F90 @ 3682

Last change on this file since 3682 was 3660, checked in by emoisan, 9 months ago
Titan physics: For debug: initialize optical variables for the first atmospheric layer. EMo
File size: 22.4 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	! initialize k=1 level (to avoid NaN values and error in debug mode)
176	DIAG_OPTT(1,:,1) = 0.0 ! dtau
177	DIAG_OPTT(1,:,2) = 1.0 ! wbar
178	DIAG_OPTT(1,:,3) = 1.0 ! gbar
179	DIAG_OPTT(1,:,4) = 0.d0
180	DIAG_OPTT(1,:,5) = 0.d0
181
182	do NW=1,L_NSPECTI
183	! We ignore K=1...
184	do K=2,L_LEVELS
185	! int. arithmetic => gives the gcm layer index (reversed)
186	ilay = L_NLAYRAD+1 - k/2
187
188	! Optics coupled with the microphysics :
189	IF (callmufi .AND. (.NOT. uncoupl_optic_haze)) THEN
190
191	!==========================================================================================
192	! Old optics (must have callclouds = .false.):
193	!==========================================================================================
194	IF (.NOT. opt4clouds) THEN
195	m3as = pqmo(ilay,2) / 2.0
196	m3af = pqmo(ilay,4) / 2.0
197	! Cut-off
198	IF (ilay .lt. ilev_cutoff) THEN
199	m3as = pqmo(ilev_cutoff,2) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(ilev_cutoff+1)-zlev(ilev_cutoff))
200	m3af = pqmo(ilev_cutoff,4) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(ilev_cutoff+1)-zlev(ilev_cutoff))
201	ENDIF
202
203	dtauaer_s = m3as*rhoaer_s(nw)
204	dtauaer_f = m3af*rhoaer_f(nw)
205
206	!==========================================================================================
207	! New optics :
208	!==========================================================================================
209	ELSE
210	iw = (L_NSPECTI + 1) - NW + L_NSPECTV ! Visible first and return
211	!-----------------------------
212	! HAZE (Spherical + Fractal) :
213	!-----------------------------
214	FTYPE = 1
215
216	! Spherical aerosols :
217	!---------------------
218	CTYPE = 5
219	m0as = pqmo(ilay,1) / 2.0
220	m3as = pqmo(ilay,2) / 2.0
221	! If not callclouds : must have a cut-off
222	IF (.NOT. callclouds) THEN
223	IF (ilay .lt. ilev_cutoff) THEN
224	m0as = pqmo(ilev_cutoff,1) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(ilev_cutoff+1)-zlev(ilev_cutoff))
225	m3as = pqmo(ilev_cutoff,2) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(ilev_cutoff+1)-zlev(ilev_cutoff))
226	ENDIF
227	ENDIF
228	call get_haze_and_cloud_opacity(FTYPE,CTYPE,m0as,m3as,iw,dtauaer_s,ssa_s(nw),asf_s(nw))
229
230	! Fractal aerosols :
231	!-------------------
232	CTYPE = FTYPE
233	m0af = pqmo(ilay,3) / 2.0
234	m3af = pqmo(ilay,4) / 2.0
235	! If not callclouds : must have a cut-off
236	IF (.NOT. callclouds) THEN
237	IF (ilay .lt. ilev_cutoff) THEN
238	m0af = pqmo(ilev_cutoff,3) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(ilev_cutoff+1)-zlev(ilev_cutoff))
239	m3af = pqmo(ilev_cutoff,4) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(ilev_cutoff+1)-zlev(ilev_cutoff))
240	ENDIF
241	ENDIF
242	call get_haze_and_cloud_opacity(FTYPE,CTYPE,m0af,m3af,iw,dtauaer_f,ssa_f(nw),asf_f(nw))
243	ENDIF
244
245	! Tuning of optical properties for haze :
246	!dtauaer_s = dtauaer_s * (FHIR * (1-ssa_s(nw)) + ssa_s(nw))
247	!ssa_s(nw) = ssa_s(nw) / (FHIR * (1-ssa_s(nw)) + ssa_s(nw))
248	dtauaer_f = dtauaer_f * (FHIR * (1-ssa_f(nw)) + ssa_f(nw))
249	ssa_f(nw) = ssa_f(nw) / (FHIR * (1-ssa_f(nw)) + ssa_f(nw))
250
251	! Total of Haze opacity (dtau), SSA (w) and ASF (COS) :
252	DHAZE_T(k,nw) = dtauaer_s + dtauaer_f
253	IF (dtauaer_s + dtauaer_f .GT. 1.D-30) THEN
254	SSA_T(k,nw) = ( dtauaer_sssa_s(nw) + dtauaer_fssa_f(nw) ) / ( dtauaer_s+dtauaer_f )
255	ASF_T(k,nw) = ( dtauaer_sssa_s(nw)asf_s(nw) + dtauaer_fssa_f(nw)asf_f(nw) ) &
256	/ ( ssa_s(nw)dtauaer_s + ssa_f(nw)dtauaer_f )
257	ELSE
258	DHAZE_T(k,nw) = 0.D0
259	SSA_T(k,nw) = 1.0
260	ASF_T(k,nw) = 1.0
261	ENDIF
262
263	! Opacity and albedo adjustement below cutoff :
264	IF (.NOT. callclouds) THEN
265	corr_haze=0.6-0.4TANH((PMID(K)100.-2500.)/250.)
266	IF (ilay .lt. ilev_cutoff) THEN
267	DHAZE_T(k,nw) = DHAZE_T(k,nw) * corr_haze
268	ENDIF
269	ENDIF
270
271	! Diagnostics for the haze :
272	DIAG_OPTH(k,nw,1) = DHAZE_T(k,nw) ! dtau
273	DIAG_OPTH(k,nw,2) = SSA_T(k,nw) ! wbar
274	DIAG_OPTH(k,nw,3) = ASF_T(k,nw) ! gbar
275
276	!---------------------
277	! CLOUDS (Spherical) :
278	!---------------------
279	IF (callclouds) THEN
280	CTYPE = 0
281	m0ccn = pqmo(ilay,5) / 2.0
282	m3ccn = pqmo(ilay,6) / 2.0
283	m3cld = pqmo(ilay,6) / 2.0
284
285	! Clear / Dark column method :
286	!-----------------------------
287
288	! CCN's SSA :
289	call get_haze_and_cloud_opacity(FTYPE,FTYPE,m0ccn,m3ccn,iw,dtau_ccn,ssa_ccn(nw),asf_ccn(nw))
290
291	! Clear column (CCN, C2H2, C2H6, HCN, AC6H6) :
292	IF (CDCOLUMN == 0) THEN
293	DO iq = 2, 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	! Dark column (CCN, CH4, C2H2, C2H6, HCN, AC6H6) :
299	ELSEIF (CDCOLUMN == 1) THEN
300	DO iq = 1, nice
301	m3cld = m3cld + (pqmo(ilay,ices_indx(iq)) / 2.0)
302	ENDDO
303	call get_haze_and_cloud_opacity(FTYPE,CTYPE,m0ccn,m3cld,iw,dtau_cld,ssa_cld(nw),asf_cld(nw))
304
305	ELSE
306	WRITE(,) 'WARNING OPTCI.F90 : CDCOLUMN MUST BE 0 OR 1'
307	WRITE(,) 'WE USE DARK COLUMN ...'
308	DO iq = 1, nice
309	m3cld = m3cld + (pqmo(ilay,ices_indx(iq)) / 2.0)
310	ENDDO
311	call get_haze_and_cloud_opacity(FTYPE,CTYPE,m0ccn,m3cld,iw,dtau_cld,ssa_cld(nw),asf_cld(nw))
312	ENDIF
313
314	! For small dropplets, opacity of nucleus dominates...
315	dtau_cld = (dtau_cldm3ccn + dtau_cldm3cld) / (m3ccn + m3cld)
316	ssa_cld(nw) = (ssa_ccn(nw)m3ccn + ssa_cld(nw)m3cld) / (m3ccn + m3cld)
317
318	! Tuning of optical properties for clouds :
319	dtau_cld = dtau_cld * (FCIR * (1-ssa_cld(nw)) + ssa_cld(nw))
320	ssa_cld(nw) = ssa_cld(nw) / (FCIR * (1-ssa_cld(nw)) + ssa_cld(nw))
321
322	! Total of Haze + Clouds opacity (dtau), SSA (w) and ASF (COS) :
323	DHAZE_T(k,nw) = dtauaer_s + dtauaer_f + dtau_cld
324	IF (DHAZE_T(k,nw) .GT. 1.D-30) THEN
325	SSA_T(k,nw) = ( dtauaer_sssa_s(nw) + dtauaer_fssa_f(nw) + dtau_cld*ssa_cld(nw) ) / ( dtauaer_s+dtauaer_f+dtau_cld )
326	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) ) &
327	/ ( ssa_s(nw)dtauaer_s + ssa_f(nw)dtauaer_f + ssa_cld(nw)*dtau_cld )
328	ELSE
329	DHAZE_T(k,nw) = 0.D0
330	SSA_T(k,nw) = 1.0
331	ASF_T(k,nw) = 1.0
332	ENDIF
333
334	! Diagnostics for clouds :
335	DIAG_OPTT(k,nw,1) = DHAZE_T(k,nw) ! dtau
336	DIAG_OPTT(k,nw,2) = SSA_T(k,nw) ! wbar
337	DIAG_OPTT(k,nw,3) = ASF_T(k,nw) ! gbar
338
339	ELSE
340	! Diagnostics for clouds :
341	DIAG_OPTT(k,nw,1) = 0.D0 ! dtau
342	DIAG_OPTT(k,nw,2) = 1.0 ! wbar
343	DIAG_OPTT(k,nw,3) = 1.0 ! gbar
344	ENDIF
345
346	! Optics and microphysics no coupled :
347	ELSE
348	! Call fixed vertical haze profile of extinction - same for all columns
349	call disr_haze(dz(k),plev(k),wnoi(nw),DHAZE_T(k,nw),SSA_T(k,nw),ASF_T(k,nw))
350	if (seashaze) DHAZE_T(k,nw) = DHAZE_T(k,nw)*seashazefact(k)
351	! Diagnostics for the haze :
352	DIAG_OPTH(k,nw,1) = DHAZE_T(k,nw) ! dtau
353	DIAG_OPTH(k,nw,2) = SSA_T(k,nw) ! wbar
354	DIAG_OPTH(k,nw,3) = ASF_T(k,nw) ! gbar
355	! Diagnostics for clouds :
356	DIAG_OPTT(k,nw,1) = 0.D0 ! dtau
357	DIAG_OPTT(k,nw,2) = 1.0 ! wbar
358	DIAG_OPTT(k,nw,3) = 1.0 ! gbar
359	ENDIF ! ENDIF callmufi
360
361	DCONT = 0.0d0 ! continuum absorption
362
363	if(continuum.and.(.not.graybody))then
364	! include continua if necessary
365	wn_cont = dble(wnoi(nw))
366	T_cont = dble(TMID(k))
367	do igas=1,ngasmx
368
369	p_cont = dble(PMID(k)scalepgfrac(igas,ilay))
370
371	dtemp=0.0d0
372	if(igas.eq.igas_N2)then
373
374	interm = indi(nw,igas,igas)
375	call interpolateN2N2(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
376	indi(nw,igas,igas) = interm
377
378	elseif(igas.eq.igas_H2)then
379
380	! first do self-induced absorption
381	interm = indi(nw,igas,igas)
382	call interpolateH2H2(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
383	indi(nw,igas,igas) = interm
384
385	! then cross-interactions with other gases
386	do jgas=1,ngasmx
387	p_cross = dble(PMID(k)scalepgfrac(jgas,ilay))
388	dtempc = 0.0d0
389	if(jgas.eq.igas_N2)then
390	interm = indi(nw,igas,jgas)
391	call interpolateN2H2(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
392	indi(nw,igas,jgas) = interm
393	endif
394	dtemp = dtemp + dtempc
395	enddo
396
397	elseif(igas.eq.igas_CH4)then
398
399	! first do self-induced absorption
400	interm = indi(nw,igas,igas)
401	call interpolateCH4CH4(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
402	indi(nw,igas,igas) = interm
403
404	! then cross-interactions with other gases
405	do jgas=1,ngasmx
406	p_cross = dble(PMID(k)scalepgfrac(jgas,ilay))
407	dtempc = 0.0d0
408	if(jgas.eq.igas_N2)then
409	interm = indi(nw,igas,jgas)
410	call interpolateN2CH4(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
411	indi(nw,igas,jgas) = interm
412	endif
413	dtemp = dtemp + dtempc
414	enddo
415
416	endif
417
418	DCONT = DCONT + dtemp
419
420	enddo
421
422	DCONT = DCONT*dz(k)
423
424	endif
425
426	do ng=1,L_NGAUSS-1
427
428	! Now compute TAUGAS
429
430	! JVO 2017 : added tmpk because the repeated calls to gasi/v increased dramatically
431	! the execution time of optci/v -> ~ factor 2 on the whole radiative
432	! transfer on the tested simulations !
433
434	if (corrk_recombin) then
435	tmpk = GASI_RECOMB(MT(K):MT(K)+1,MP(K):MP(K)+1,NW,NG)
436	else
437	tmpk = GASI(MT(K):MT(K)+1,MP(K):MP(K)+1,1,NW,NG)
438	endif
439
440	KCOEF(1) = tmpk(1,1) ! KCOEF(1) = GASI(MT(K),MP(K),1,NW,NG)
441	KCOEF(2) = tmpk(1,2) ! KCOEF(2) = GASI(MT(K),MP(K)+1,1,NW,NG)
442	KCOEF(3) = tmpk(2,2) ! KCOEF(3) = GASI(MT(K)+1,MP(K)+1,1,NW,NG)
443	KCOEF(4) = tmpk(2,1) ! KCOEF(4) = GASI(MT(K)+1,MP(K),1,NW,NG)
444
445
446	! Interpolate the gaseous k-coefficients to the requested T,P values
447
448	ANS = LKCOEF(K,1)KCOEF(1) + LKCOEF(K,2)KCOEF(2) + &
449	LKCOEF(K,3)KCOEF(3) + LKCOEF(K,4)KCOEF(4)
450
451	TAUGAS = U(k)*ANS
452
453	TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS + DCONT
454	DTAUKI(K,nw,ng) = TAUGAS &
455	+ DCONT & ! For parameterized continuum absorption
456	+ DHAZE_T(K,NW) ! For Titan haze
457
458	end do
459
460	! Now fill in the "clear" part of the spectrum (NG = L_NGAUSS),
461	! which holds continuum opacity only
462
463	NG = L_NGAUSS
464	DTAUKI(K,nw,ng) = 0.d0 &
465	+ DCONT & ! For parameterized continuum absorption
466	+ DHAZE_T(K,NW) ! For Titan Haze
467
468	DIAG_OPTH(K,nw,4) = 0.d0
469	DIAG_OPTH(K,nw,5) = TAUGAS
470	DIAG_OPTH(K,nw,6) = DCONT
471	DIAG_OPTT(K,nw,4) = 0.d0
472	DIAG_OPTT(K,nw,5) = TAUGAS
473	DIAG_OPTT(K,nw,6) = DCONT
474
475	end do ! k = L_LEVELS
476	end do ! nw = L_NSPECTI
477
478	!=======================================================================
479	! Now the full treatment for the layers, where besides the opacity
480	! we need to calculate the scattering albedo and asymmetry factors
481	! ======================================================================
482
483	! Haze scattering
484	DO NW=1,L_NSPECTI
485	DO K=2,L_LEVELS
486	DHAZES_T(K,NW) = DHAZE_T(K,NW) * SSA_T(K,NW)
487	ENDDO
488	ENDDO
489
490	DO NW=1,L_NSPECTI
491	DO L=1,L_NLAYRAD-1
492	K = 2*L+1
493	btemp(L,NW) = DHAZES_T(K,NW) + DHAZES_T(K+1,NW)
494	END DO ! L vertical loop
495
496	! Last level
497	L = L_NLAYRAD
498	K = 2*L+1
499	btemp(L,NW) = DHAZES_T(K,NW)
500
501	END DO ! NW spectral loop
502
503
504	DO NW=1,L_NSPECTI
505	NG = L_NGAUSS
506	DO L=1,L_NLAYRAD-1
507
508	K = 2*L+1
509	DTAUI(L,nw,ng) = DTAUKI(K,NW,NG) + DTAUKI(K+1,NW,NG)! + 1.e-50
510
511	atemp = 0.
512	if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
513	atemp = atemp + &
514	ASF_T(K,NW)*DHAZES_T(K,NW) + &
515	ASF_T(K+1,NW)*DHAZES_T(K+1,NW)
516
517	WBARI(L,nw,ng) = btemp(L,nw) / DTAUI(L,NW,NG)
518	else
519	WBARI(L,nw,ng) = 0.0D0
520	DTAUI(L,NW,NG) = 1.0D-9
521	endif
522
523	if(btemp(L,nw) .GT. 0.0d0) then
524	cosbi(L,NW,NG) = atemp/btemp(L,nw)
525	else
526	cosbi(L,NW,NG) = 0.0D0
527	end if
528
529	END DO ! L vertical loop
530
531	! Last level
532
533	L = L_NLAYRAD
534	K = 2*L+1
535	DTAUI(L,nw,ng) = DTAUKI(K,NW,NG) ! + 1.e-50
536
537	atemp = 0.
538	if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
539	atemp = atemp + ASF_T(K,NW)*DHAZES_T(K,NW)
540	WBARI(L,nw,ng) = btemp(L,nw) / DTAUI(L,NW,NG)
541	else
542	WBARI(L,nw,ng) = 0.0D0
543	DTAUI(L,NW,NG) = 1.0D-9
544	endif
545
546	if(btemp(L,nw) .GT. 0.0d0) then
547	cosbi(L,NW,NG) = atemp/btemp(L,nw)
548	else
549	cosbi(L,NW,NG) = 0.0D0
550	end if
551
552
553	! Now the other Gauss points, if needed.
554
555	DO NG=1,L_NGAUSS-1
556	IF(TAUGSURF(NW,NG) .gt. TLIMIT) THEN
557
558	DO L=1,L_NLAYRAD-1
559	K = 2*L+1
560	DTAUI(L,nw,ng) = DTAUKI(K,NW,NG)+DTAUKI(K+1,NW,NG)! + 1.e-50
561
562	if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
563
564	WBARI(L,nw,ng) = btemp(L,nw) / DTAUI(L,NW,NG)
565
566	else
567	WBARI(L,nw,ng) = 0.0D0
568	DTAUI(L,NW,NG) = 1.0D-9
569	endif
570
571	cosbi(L,NW,NG) = cosbi(L,NW,L_NGAUSS)
572	END DO ! L vertical loop
573
574	! Last level
575	L = L_NLAYRAD
576	K = 2*L+1
577	DTAUI(L,nw,ng) = DTAUKI(K,NW,NG)! + 1.e-50
578
579	if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
580
581	WBARI(L,nw,ng) = btemp(L,nw) / DTAUI(L,NW,NG)
582
583	else
584	WBARI(L,nw,ng) = 0.0D0
585	DTAUI(L,NW,NG) = 1.0D-9
586	endif
587
588	cosbi(L,NW,NG) = cosbi(L,NW,L_NGAUSS)
589
590	END IF
591
592	END DO ! NG Gauss loop
593	END DO ! NW spectral loop
594
595	! Total extinction optical depths
596	!DO NG=1,L_NGAUSS ! full gauss loop
597	! DO NW=1,L_NSPECTI
598	! TAUCUMI(1,NW,NG)=0.0D0
599	! DO K=2,L_LEVELS
600	! TAUCUMI(K,NW,NG)=TAUCUMI(K-1,NW,NG)+DTAUKI(K,NW,NG)
601	! END DO
602	! END DO ! end full gauss loop
603	!END DO
604
605	TAUCUMI(:,:,:) = DTAUKI(:,:,:)
606
607	! be aware when comparing with textbook results
608	! (e.g. Pierrehumbert p. 218) that
609	! taucumi does not take the <cos theta>=0.5 factor into
610	! account. It is the optical depth for a vertically
611	! ascending ray with angle theta = 0.
612
613	if(firstcall) firstcall = .false.
614
615	return
616
617
618	end subroutine optci
619
620
621

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

Download in other formats:

Original Format