Context Navigation

optci.F90 @ 3094

Last change on this file since 3094 was 3090, checked in by slebonnois, 15 months ago
BdeBatz? : Cleans microphysics and makes few corrections for physics
File size: 21.0 KB

Line
1	subroutine optci(PQMO,NLAY,ZLEV,PLEV,TLEV,TMID,PMID, &
2	DTAUI,TAUCUMI,COSBI,WBARI,TAUGSURF,SEASHAZEFACT,&
3	POPT_HAZE,POPT_CLOUDS,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
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 to Titan
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) :: POPT_HAZE(L_LEVELS,L_NSPECTI,3)
64	REAL*8, INTENT(OUT) :: POPT_CLOUDS(L_LEVELS,L_NSPECTI,3)
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
109	! Variables for new optics
110	integer iq, iw, FTYPE, CTYPE
111	real*8 m0as,m0af,m0ccn,m3as,m3af,m3ccn,m3cld
112	real*8 dtauaer_s,dtauaer_f,dtau_ccn,dtau_cld
113	real*8,save :: rhoaer_s(L_NSPECTI),ssa_s(L_NSPECTI),asf_s(L_NSPECTI)
114	real*8,save :: rhoaer_f(L_NSPECTI),ssa_f(L_NSPECTI),asf_f(L_NSPECTI)
115	real*8,save :: ssa_ccn(L_NSPECTI),asf_ccn(L_NSPECTI)
116	real*8,save :: ssa_cld(L_NSPECTI),asf_cld(L_NSPECTI)
117	!$OMP THREADPRIVATE(rhoaer_s,rhoaer_f,ssa_s,ssa_f,ssa_cld,asf_s,asf_f,asf_cld)
118
119	logical,save :: firstcall=.true.
120	!$OMP THREADPRIVATE(firstcall)
121
122
123	!! AS: to save time in computing continuum (see bilinearbig)
124	IF (.not.ALLOCATED(indi)) THEN
125	ALLOCATE(indi(L_NSPECTI,ngasmx,ngasmx))
126	indi = -9999 ! this initial value means "to be calculated"
127	ENDIF
128
129	! Some initialisation because there's a pb with disr_haze at the limits (nw=1)
130	! I should check this - For now we set vars to zero : better than nans - JVO 2017
131	DHAZE_T(:,:) = 0.0
132	SSA_T(:,:) = 0.0
133	ASF_T(:,:) = 0.0
134
135	! Load tabulated haze optical properties if needed.
136	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
137	IF (firstcall .AND. callmufi .AND. (.NOT. uncoupl_optic_haze)) THEN
138	OPEN(12,file=TRIM(haze_opt_file),form='formatted') ! The file has been inquired in physiq_mod firstcall
139	READ(12,*) ! dummy header
140	DO NW=1,L_NSPECTI
141	READ(12,*) dumch, dumwvl, rhoaer_f(nw), ssa_f(nw), asf_f(nw), rhoaer_s(nw), ssa_s(nw), asf_s(nw)
142	ENDDO
143	CLOSE(12)
144	ENDIF
145
146	!=======================================================================
147	! Determine the total gas opacity throughout the column, for each
148	! spectral interval, NW, and each Gauss point, NG.
149
150	taugsurf(:,:) = 0.0
151	dpr(:) = 0.0
152	lkcoef(:,:) = 0.0
153
154	do K=2,L_LEVELS
155	ilay = L_NLAYRAD+1 - k/2 ! int. arithmetic => gives the gcm layer index (reversed)
156	DPR(k) = PLEV(K)-PLEV(K-1)
157
158	! if we have continuum opacities, we need dz
159	dz(k) = dpr(k)RTMID(K)/(gzlat_ig(ilay)*PMID(K))
160	U(k) = Cmk(ilay)*DPR(k) ! only Cmk line in optci.F
161
162	call tpindex(PMID(K),TMID(K),pfgasref,tgasref,LCOEF,MT(K),MP(K))
163
164	do LK=1,4
165	LKCOEF(K,LK) = LCOEF(LK)
166	end do
167	end do ! L_LEVELS
168
169	do NW=1,L_NSPECTI
170	! We ignore K=1...
171	do K=2,L_LEVELS
172	! int. arithmetic => gives the gcm layer index (reversed)
173	ilay = L_NLAYRAD+1 - k/2
174
175	! Optics coupled with the microphysics :
176	IF (callmufi .AND. (.NOT. uncoupl_optic_haze)) THEN
177
178	!==========================================================================================
179	! Old optics (must have callclouds = .false.):
180	!==========================================================================================
181	IF (.NOT. opt4clouds) THEN
182	m3as = pqmo(ilay,2) / 2.0
183	m3af = pqmo(ilay,4) / 2.0
184	! Cut-off (here for p = 2.7e3Pa / alt = 70km)
185	IF (ilay .lt. 23) THEN
186	m3as = pqmo(23,2) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(24)-zlev(23))
187	m3af = pqmo(23,4) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(24)-zlev(23))
188	ENDIF
189
190	dtauaer_s = m3as*rhoaer_s(nw)
191	dtauaer_f = m3af*rhoaer_f(nw)
192
193	!==========================================================================================
194	! New optics :
195	!==========================================================================================
196	ELSE
197	iw = (L_NSPECTI + 1) - NW + L_NSPECTV ! Visible first and return
198	!-----------------------------
199	! HAZE (Spherical + Fractal) :
200	!-----------------------------
201	FTYPE = 1
202
203	! Spherical aerosols :
204	!---------------------
205	CTYPE = 5
206	m0as = pqmo(ilay,1) / 2.0
207	m3as = pqmo(ilay,2) / 2.0
208	! If not callclouds : must have a cut-off (here for p = 2.7e3Pa / alt = 70km)
209	IF (.NOT. callclouds) THEN
210	IF (ilay .lt. 23) THEN
211	m0as = pqmo(23,1) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(24)-zlev(23))
212	m3as = pqmo(23,2) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(24)-zlev(23))
213	ENDIF
214	ENDIF
215	call get_haze_and_cloud_opacity(FTYPE,CTYPE,m0as,m3as,iw,dtauaer_s,ssa_s(nw),asf_s(nw))
216
217	! Fractal aerosols :
218	!-------------------
219	CTYPE = FTYPE
220	m0af = pqmo(ilay,3) / 2.0
221	m3af = pqmo(ilay,4) / 2.0
222	! If not callclouds : must have a cut-off (here for p = 2.7e3Pa / alt = 70km)
223	IF (.NOT. callclouds) THEN
224	IF (ilay .lt. 23) THEN
225	m0af = pqmo(23,3) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(24)-zlev(23))
226	m3af = pqmo(23,4) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(24)-zlev(23))
227	ENDIF
228	ENDIF
229	call get_haze_and_cloud_opacity(FTYPE,CTYPE,m0af,m3af,iw,dtauaer_f,ssa_f(nw),asf_f(nw))
230	ENDIF
231
232	! Total of Haze opacity (dtau), SSA (w) and ASF (COS) :
233	DHAZE_T(k,nw) = dtauaer_s + dtauaer_f
234	IF (dtauaer_s + dtauaer_f .GT. 1.D-30) THEN
235	SSA_T(k,nw) = ( dtauaer_sssa_s(nw) + dtauaer_fssa_f(nw) ) / ( dtauaer_s+dtauaer_f )
236	ASF_T(k,nw) = ( dtauaer_sssa_s(nw)asf_s(nw) + dtauaer_fssa_f(nw)asf_f(nw) ) &
237	/ ( ssa_s(nw)dtauaer_s + ssa_f(nw)dtauaer_f )
238	ELSE
239	DHAZE_T(k,nw) = 0.D0
240	SSA_T(k,nw) = 1.0
241	ASF_T(k,nw) = 1.0
242	ENDIF
243	! Diagnostics for the haze :
244	POPT_HAZE(k,nw,1) = DHAZE_T(k,nw) ! dtau
245	POPT_HAZE(k,nw,2) = SSA_T(k,nw) ! wbar
246	POPT_HAZE(k,nw,3) = ASF_T(k,nw) ! gbar
247
248	!---------------------
249	! CLOUDS (Spherical) :
250	!---------------------
251	IF (callclouds) THEN
252	CTYPE = 0
253	m0ccn = pqmo(ilay,5) / 2.0
254	m3ccn = pqmo(ilay,6) / 2.0
255	m3cld = 0.0d0
256
257	! Clear / Dark column method :
258	!-----------------------------
259
260	! CCN's SSA :
261	call get_haze_and_cloud_opacity(FTYPE,FTYPE,m0ccn,m3ccn,iw,dtau_ccn,ssa_ccn(nw),asf_ccn(nw))
262
263	! Clear column (CCN, C2H2, C2H6, HCN) :
264	IF (CDCOLUMN == 0) THEN
265	DO iq = 2, nice
266	m3cld = m3cld + (pqmo(ilay,ices_indx(iq)) / 2.0)
267	ENDDO
268	call get_haze_and_cloud_opacity(FTYPE,CTYPE,m0ccn,m3cld,iw,dtau_cld,ssa_cld(nw),asf_cld(nw))
269
270	! Dark column (CCN, CH4, C2H2, C2H6, HCN) :
271	ELSEIF (CDCOLUMN == 1) THEN
272	DO iq = 1, nice
273	m3cld = m3cld + (pqmo(ilay,ices_indx(iq)) / 2.0)
274	ENDDO
275	call get_haze_and_cloud_opacity(FTYPE,CTYPE,m0ccn,m3cld,iw,dtau_cld,ssa_cld(nw),asf_cld(nw))
276
277	ELSE
278	WRITE(,) 'WARNING OPTCI.F90 : CDCOLUMN MUST BE 0 OR 1'
279	WRITE(,) 'WE USE DARK COLUMN ...'
280	DO iq = 1, nice
281	m3cld = m3cld + (pqmo(ilay,ices_indx(iq)) / 2.0)
282	ENDDO
283	call get_haze_and_cloud_opacity(FTYPE,CTYPE,m0ccn,m3cld,iw,dtau_cld,ssa_cld(nw),asf_cld(nw))
284	ENDIF
285
286	! For small dropplets, opacity of nucleus dominates...
287	ssa_cld(nw) = (ssa_ccn(nw)m3ccn + ssa_cld(nw)m3cld) / (m3ccn + m3cld)
288
289	! Total of Haze + Clouds opacity (dtau), SSA (w) and ASF (COS) :
290	DHAZE_T(k,nw) = dtauaer_s + dtauaer_f + dtau_cld
291	IF (DHAZE_T(k,nw) .GT. 1.D-30) THEN
292	SSA_T(k,nw) = ( dtauaer_sssa_s(nw) + dtauaer_fssa_f(nw) + dtau_cld*ssa_cld(nw) ) / ( dtauaer_s+dtauaer_f+dtau_cld )
293	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) ) &
294	/ ( ssa_s(nw)dtauaer_s + ssa_f(nw)dtauaer_f + ssa_cld(nw)*dtau_cld )
295	ELSE
296	DHAZE_T(k,nw) = 0.D0
297	SSA_T(k,nw) = 1.0
298	ASF_T(k,nw) = 1.0
299	ENDIF
300
301	! Tuning of optical properties (now useless...) :
302	DHAZE_T(k,nw) = DHAZE_T(k,nw) * (FHIR * (1-SSA_T(k,nw)) + SSA_T(k,nw))
303	SSA_T(k,nw) = SSA_T(k,nw) / (FHIR * (1-SSA_T(k,nw)) + SSA_T(k,nw))
304
305	! Diagnostics for clouds :
306	POPT_CLOUDS(k,nw,1) = DHAZE_T(k,nw) ! dtau
307	POPT_CLOUDS(k,nw,2) = SSA_T(k,nw) ! wbar
308	POPT_CLOUDS(k,nw,3) = ASF_T(k,nw) ! gbar
309
310	ELSE
311	! Diagnostics for clouds :
312	POPT_CLOUDS(k,nw,1) = 0.D0 ! dtau
313	POPT_CLOUDS(k,nw,2) = 1.0 ! wbar
314	POPT_CLOUDS(k,nw,3) = 1.0 ! gbar
315	ENDIF
316
317	! Optics and microphysics no coupled :
318	ELSE
319	! Call fixed vertical haze profile of extinction - same for all columns
320	call disr_haze(dz(k),plev(k),wnoi(nw),DHAZE_T(k,nw),SSA_T(k,nw),ASF_T(k,nw))
321	if (seashaze) DHAZE_T(k,nw) = DHAZE_T(k,nw)*seashazefact(k)
322	! Diagnostics for the haze :
323	POPT_HAZE(k,nw,1) = DHAZE_T(k,nw) ! dtau
324	POPT_HAZE(k,nw,2) = SSA_T(k,nw) ! wbar
325	POPT_HAZE(k,nw,3) = ASF_T(k,nw) ! gbar
326	! Diagnostics for clouds :
327	POPT_CLOUDS(k,nw,1) = 0.D0 ! dtau
328	POPT_CLOUDS(k,nw,2) = 1.0 ! wbar
329	POPT_CLOUDS(k,nw,3) = 1.0 ! gbar
330	ENDIF ! ENDIF callmufi
331
332	DCONT = 0.0d0 ! continuum absorption
333
334	if(continuum.and.(.not.graybody))then
335	! include continua if necessary
336	wn_cont = dble(wnoi(nw))
337	T_cont = dble(TMID(k))
338	do igas=1,ngasmx
339
340	p_cont = dble(PMID(k)scalepgfrac(igas,ilay))
341
342	dtemp=0.0d0
343	if(igas.eq.igas_N2)then
344
345	interm = indi(nw,igas,igas)
346	call interpolateN2N2(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
347	indi(nw,igas,igas) = interm
348
349	elseif(igas.eq.igas_H2)then
350
351	! first do self-induced absorption
352	interm = indi(nw,igas,igas)
353	call interpolateH2H2(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
354	indi(nw,igas,igas) = interm
355
356	! then cross-interactions with other gases
357	do jgas=1,ngasmx
358	p_cross = dble(PMID(k)scalepgfrac(jgas,ilay))
359	dtempc = 0.0d0
360	if(jgas.eq.igas_N2)then
361	interm = indi(nw,igas,jgas)
362	call interpolateN2H2(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
363	indi(nw,igas,jgas) = interm
364	endif
365	dtemp = dtemp + dtempc
366	enddo
367
368	elseif(igas.eq.igas_CH4)then
369
370	! first do self-induced absorption
371	interm = indi(nw,igas,igas)
372	call interpolateCH4CH4(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
373	indi(nw,igas,igas) = interm
374
375	! then cross-interactions with other gases
376	do jgas=1,ngasmx
377	p_cross = dble(PMID(k)scalepgfrac(jgas,ilay))
378	dtempc = 0.0d0
379	if(jgas.eq.igas_N2)then
380	interm = indi(nw,igas,jgas)
381	call interpolateN2CH4(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
382	indi(nw,igas,jgas) = interm
383	endif
384	dtemp = dtemp + dtempc
385	enddo
386
387	endif
388
389	DCONT = DCONT + dtemp
390
391	enddo
392
393	DCONT = DCONT*dz(k)
394
395	endif
396
397	do ng=1,L_NGAUSS-1
398
399	! Now compute TAUGAS
400
401	! JVO 2017 : added tmpk because the repeated calls to gasi/v increased dramatically
402	! the execution time of optci/v -> ~ factor 2 on the whole radiative
403	! transfer on the tested simulations !
404
405	if (corrk_recombin) then
406	tmpk = GASI_RECOMB(MT(K):MT(K)+1,MP(K):MP(K)+1,NW,NG)
407	else
408	tmpk = GASI(MT(K):MT(K)+1,MP(K):MP(K)+1,1,NW,NG)
409	endif
410
411	KCOEF(1) = tmpk(1,1) ! KCOEF(1) = GASI(MT(K),MP(K),1,NW,NG)
412	KCOEF(2) = tmpk(1,2) ! KCOEF(2) = GASI(MT(K),MP(K)+1,1,NW,NG)
413	KCOEF(3) = tmpk(2,2) ! KCOEF(3) = GASI(MT(K)+1,MP(K)+1,1,NW,NG)
414	KCOEF(4) = tmpk(2,1) ! KCOEF(4) = GASI(MT(K)+1,MP(K),1,NW,NG)
415
416
417	! Interpolate the gaseous k-coefficients to the requested T,P values
418
419	ANS = LKCOEF(K,1)KCOEF(1) + LKCOEF(K,2)KCOEF(2) + &
420	LKCOEF(K,3)KCOEF(3) + LKCOEF(K,4)KCOEF(4)
421
422	TAUGAS = U(k)*ANS
423
424	TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS + DCONT
425	DTAUKI(K,nw,ng) = TAUGAS &
426	+ DCONT & ! For parameterized continuum absorption
427	+ DHAZE_T(K,NW) ! For Titan haze
428
429	end do
430
431	! Now fill in the "clear" part of the spectrum (NG = L_NGAUSS),
432	! which holds continuum opacity only
433
434	NG = L_NGAUSS
435	DTAUKI(K,nw,ng) = 0.d0 &
436	+ DCONT & ! For parameterized continuum absorption
437	+ DHAZE_T(K,NW) ! For Titan Haze
438
439	end do ! k = L_LEVELS
440	end do ! nw = L_NSPECTI
441
442	!=======================================================================
443	! Now the full treatment for the layers, where besides the opacity
444	! we need to calculate the scattering albedo and asymmetry factors
445	! ======================================================================
446
447	! Haze scattering
448	DO NW=1,L_NSPECTI
449	DO K=2,L_LEVELS
450	DHAZES_T(K,NW) = DHAZE_T(K,NW) * SSA_T(K,NW)
451	ENDDO
452	ENDDO
453
454	DO NW=1,L_NSPECTI
455	DO L=1,L_NLAYRAD-1
456	K = 2*L+1
457	btemp(L,NW) = DHAZES_T(K,NW) + DHAZES_T(K+1,NW)
458	END DO ! L vertical loop
459
460	! Last level
461	L = L_NLAYRAD
462	K = 2*L+1
463	btemp(L,NW) = DHAZES_T(K,NW)
464
465	END DO ! NW spectral loop
466
467
468	DO NW=1,L_NSPECTI
469	NG = L_NGAUSS
470	DO L=1,L_NLAYRAD-1
471
472	K = 2*L+1
473	DTAUI(L,nw,ng) = DTAUKI(K,NW,NG) + DTAUKI(K+1,NW,NG)! + 1.e-50
474
475	atemp = 0.
476	if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
477	atemp = atemp + &
478	ASF_T(K,NW)*DHAZES_T(K,NW) + &
479	ASF_T(K+1,NW)*DHAZES_T(K+1,NW)
480
481	WBARI(L,nw,ng) = btemp(L,nw) / DTAUI(L,NW,NG)
482	else
483	WBARI(L,nw,ng) = 0.0D0
484	DTAUI(L,NW,NG) = 1.0D-9
485	endif
486
487	if(btemp(L,nw) .GT. 0.0d0) then
488	cosbi(L,NW,NG) = atemp/btemp(L,nw)
489	else
490	cosbi(L,NW,NG) = 0.0D0
491	end if
492
493	END DO ! L vertical loop
494
495	! Last level
496
497	L = L_NLAYRAD
498	K = 2*L+1
499	DTAUI(L,nw,ng) = DTAUKI(K,NW,NG) ! + 1.e-50
500
501	atemp = 0.
502	if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
503	atemp = atemp + ASF_T(K,NW)*DHAZES_T(K,NW)
504	WBARI(L,nw,ng) = btemp(L,nw) / DTAUI(L,NW,NG)
505	else
506	WBARI(L,nw,ng) = 0.0D0
507	DTAUI(L,NW,NG) = 1.0D-9
508	endif
509
510	if(btemp(L,nw) .GT. 0.0d0) then
511	cosbi(L,NW,NG) = atemp/btemp(L,nw)
512	else
513	cosbi(L,NW,NG) = 0.0D0
514	end if
515
516
517	! Now the other Gauss points, if needed.
518
519	DO NG=1,L_NGAUSS-1
520	IF(TAUGSURF(NW,NG) .gt. TLIMIT) THEN
521
522	DO L=1,L_NLAYRAD-1
523	K = 2*L+1
524	DTAUI(L,nw,ng) = DTAUKI(K,NW,NG)+DTAUKI(K+1,NW,NG)! + 1.e-50
525
526	if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
527
528	WBARI(L,nw,ng) = btemp(L,nw) / DTAUI(L,NW,NG)
529
530	else
531	WBARI(L,nw,ng) = 0.0D0
532	DTAUI(L,NW,NG) = 1.0D-9
533	endif
534
535	cosbi(L,NW,NG) = cosbi(L,NW,L_NGAUSS)
536	END DO ! L vertical loop
537
538	! Last level
539	L = L_NLAYRAD
540	K = 2*L+1
541	DTAUI(L,nw,ng) = DTAUKI(K,NW,NG)! + 1.e-50
542
543	if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
544
545	WBARI(L,nw,ng) = btemp(L,nw) / DTAUI(L,NW,NG)
546
547	else
548	WBARI(L,nw,ng) = 0.0D0
549	DTAUI(L,NW,NG) = 1.0D-9
550	endif
551
552	cosbi(L,NW,NG) = cosbi(L,NW,L_NGAUSS)
553
554	END IF
555
556	END DO ! NG Gauss loop
557	END DO ! NW spectral loop
558
559	! Total extinction optical depths
560	!DO NG=1,L_NGAUSS ! full gauss loop
561	! DO NW=1,L_NSPECTI
562	! TAUCUMI(1,NW,NG)=0.0D0
563	! DO K=2,L_LEVELS
564	! TAUCUMI(K,NW,NG)=TAUCUMI(K-1,NW,NG)+DTAUKI(K,NW,NG)
565	! END DO
566	! END DO ! end full gauss loop
567	!END DO
568
569	TAUCUMI(:,:,:) = DTAUKI(:,:,:)
570
571	! be aware when comparing with textbook results
572	! (e.g. Pierrehumbert p. 218) that
573	! taucumi does not take the <cos theta>=0.5 factor into
574	! account. It is the optical depth for a vertically
575	! ascending ray with angle theta = 0.
576
577	if(firstcall) firstcall = .false.
578
579	return
580
581
582	end subroutine optci
583
584
585

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

Download in other formats:

Original Format