Context Navigation

optci.F90 @ 1862

Last change on this file since 1862 was 1822, checked in by jvatant, 7 years ago
Preliminary modifs for the optical coupling of haze + Moved inits of setspi/v before init of mufi + Added access to tarcers in optci/v + Some coherence in call to directories JVO
File size: 12.7 KB

Line
1	subroutine optci(PQ,PLEV,TLEV,TMID,PMID, &
2	DTAUI,TAUCUMI,COSBI,WBARI,TAUGSURF)
3
4	use radinc_h
5	use radcommon_h, only: gasi,tlimit,Cmk,tgasref,pfgasref,wnoi,scalep,indi,glat_ig,gweight
6	use gases_h
7	use comcstfi_mod, only: g, r
8	use callkeys_mod, only: continuum,graybody
9	implicit none
10
11	!==================================================================
12	!
13	! Purpose
14	! -------
15	! Calculates longwave optical constants at each level. For each
16	! layer and spectral interval in the IR it calculates WBAR, DTAU
17	! and COSBAR. For each level it calculates TAU.
18	!
19	! TAUCUMI(L,LW) is the cumulative optical depth at level L (or alternatively
20	! at the bottom of layer L), LW is the spectral wavelength interval.
21	!
22	! TLEV(L) - Temperature at the layer boundary (i.e., level)
23	! PLEV(L) - Pressure at the layer boundary (i.e., level)
24	!
25	! Authors
26	! -------
27	! Adapted from the NASA Ames code by R. Wordsworth (2009)
28	! Clean and adaptation to Titan by J. Vatant d'Ollone (2016-17)
29	!
30	!==================================================================
31
32
33	!==========================================================
34	! Input/Output
35	!==========================================================
36	REAL*8, INTENT(IN) :: PQ ! Tracers (kg/kg_of_air).
37	REAL*8, INTENT(IN) :: PLEV(L_LEVELS), TLEV(L_LEVELS)
38	REAL*8, INTENT(IN) :: TMID(L_LEVELS), PMID(L_LEVELS)
39
40	REAL*8, INTENT(OUT) :: DTAUI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
41	REAL*8, INTENT(OUT) :: TAUCUMI(L_LEVELS,L_NSPECTI,L_NGAUSS)
42	REAL*8, INTENT(OUT) :: COSBI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
43	REAL*8, INTENT(OUT) :: WBARI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
44	REAL*8, INTENT(OUT) :: TAUGSURF(L_NSPECTI,L_NGAUSS-1)
45	! ==========================================================
46
47	real*8 DTAUKI(L_LEVELS,L_NSPECTI,L_NGAUSS)
48
49	! Titan customisation
50	! J. Vatant d'Ollone (2016)
51	real*8 DHAZE_T(L_LEVELS,L_NSPECTI)
52	real*8 DHAZES_T(L_LEVELS,L_NSPECTI)
53	real*8 SSA_T(L_LEVELS,L_NSPECTI)
54	real*8 ASF_T(L_LEVELS,L_NSPECTI)
55	real*8 INT_DTAU(L_NLAYRAD,L_NSPECTI)
56	real*8 K_HAZE(L_NLAYRAD,L_NSPECTI)
57
58	CHARACTER*2 str2
59	! ==========================
60
61	integer L, NW, NG, K, LK, IAER
62	integer MT(L_LEVELS), MP(L_LEVELS), NP(L_LEVELS)
63	real*8 ANS, TAUGAS
64	real*8 DPR(L_LEVELS), U(L_LEVELS)
65	real*8 LCOEF(4), LKCOEF(L_LEVELS,4)
66
67	real*8 DCONT
68	double precision wn_cont, p_cont, p_air, T_cont, dtemp, dtempc
69	double precision p_cross
70
71	real*8 KCOEF(4)
72
73	! temporary variable to reduce memory access time to gasi
74	real*8 tmpk(2,2)
75
76	! temporary variables for multiple aerosol calculation
77	real*8 atemp
78	real*8 btemp(L_NLAYRAD,L_NSPECTI)
79
80	! variables for k in units m^-1
81	real*8 dz(L_LEVELS)
82	!real*8 rho !! see test below
83
84	integer igas, jgas, ilay
85
86	integer interm
87
88	!! AS: to save time in computing continuum (see bilinearbig)
89	IF (.not.ALLOCATED(indi)) THEN
90	ALLOCATE(indi(L_NSPECTI,ngasmx,ngasmx))
91	indi = -9999 ! this initial value means "to be calculated"
92	ENDIF
93
94	! Some initialisation beacause there's a pb with disr_haze at the limits (nw=1)
95	! I should check this - For now we set vars to zero : better than nans - JVO 2017
96
97	dhaze_t(:,:) = 0.
98	ssa_t(:,:) = 0.
99	asf_t(:,:) = 0.
100
101	!=======================================================================
102	! Determine the total gas opacity throughout the column, for each
103	! spectral interval, NW, and each Gauss point, NG.
104
105	taugsurf(:,:) = 0.0
106	dpr(:) = 0.0
107	lkcoef(:,:) = 0.0
108
109	do K=2,L_LEVELS
110	DPR(k) = PLEV(K)-PLEV(K-1)
111
112	! if we have continuum opacities, we need dz
113	dz(k) = dpr(k)RTMID(K)/(glat_ig*PMID(K))
114	U(k) = Cmk*DPR(k) ! only Cmk line in optci.F
115
116	call tpindex(PMID(K),TMID(K),pfgasref,tgasref,LCOEF,MT(K),MP(K))
117
118	do LK=1,4
119	LKCOEF(K,LK) = LCOEF(LK)
120	end do
121	end do ! levels
122
123	do NW=1,L_NSPECTI
124
125	do K=2,L_LEVELS
126
127	ilay = k / 2 ! int. arithmetic => gives the gcm layer index
128
129	!================= Titan customisation ========================================
130	call disr_haze(dz(k),plev(k),wnoi(nw),dhaze_T(k,nw),SSA_T(k,nw),ASF_T(k,nw))
131	! =============================================================================
132
133	DCONT = 0.0d0 ! continuum absorption
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	do ng=1,L_NGAUSS-1
210
211	! Now compute TAUGAS
212
213	! JVO 2017 : added tmpk because the repeated calls to gasi/v increased dramatically
214	! the execution time of optci/v -> ~ factor 2 on the whole radiative
215	! transfer on the tested simulations !
216
217	tmpk = GASI(MT(K):MT(K)+1,MP(K):MP(K)+1,1,NW,NG)
218
219	KCOEF(1) = tmpk(1,1) ! KCOEF(1) = GASI(MT(K),MP(K),1,NW,NG)
220	KCOEF(2) = tmpk(1,2) ! KCOEF(2) = GASI(MT(K),MP(K)+1,1,NW,NG)
221	KCOEF(3) = tmpk(2,2) ! KCOEF(3) = GASI(MT(K)+1,MP(K)+1,1,NW,NG)
222	KCOEF(4) = tmpk(2,1) ! KCOEF(4) = GASI(MT(K)+1,MP(K),1,NW,NG)
223
224
225	! Interpolate the gaseous k-coefficients to the requested T,P values
226
227	ANS = LKCOEF(K,1)KCOEF(1) + LKCOEF(K,2)KCOEF(2) + &
228	LKCOEF(K,3)KCOEF(3) + LKCOEF(K,4)KCOEF(4)
229
230	TAUGAS = U(k)*ANS
231
232	TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS + DCONT
233	DTAUKI(K,nw,ng) = TAUGAS &
234	+ DCONT & ! For parameterized continuum absorption
235	+ DHAZE_T(K,NW) ! For Titan haze
236
237	end do
238
239	! Now fill in the "clear" part of the spectrum (NG = L_NGAUSS),
240	! which holds continuum opacity only
241
242	NG = L_NGAUSS
243	DTAUKI(K,nw,ng) = 0.d0 &
244	+ DCONT & ! For parameterized continuum absorption
245	+ DHAZE_T(K,NW) ! For Titan Haze
246
247	end do
248	end do
249
250	!=======================================================================
251	! Now the full treatment for the layers, where besides the opacity
252	! we need to calculate the scattering albedo and asymmetry factors
253	! ======================================================================
254
255	! Haze scattering
256	DO NW=1,L_NSPECTI
257	DO K=2,L_LEVELS
258	DHAZES_T(K,NW) = DHAZE_T(K,NW) * SSA_T(K,NW)
259	ENDDO
260	ENDDO
261
262	DO NW=1,L_NSPECTI
263	DO L=1,L_NLAYRAD-1
264	K = 2*L+1
265	btemp(L,NW) = DHAZES_T(K,NW) + DHAZES_T(K+1,NW)
266	END DO ! L vertical loop
267
268	! Last level
269	L = L_NLAYRAD
270	K = 2*L+1
271	btemp(L,NW) = DHAZES_T(K,NW)
272
273	END DO ! NW spectral loop
274
275
276	DO NW=1,L_NSPECTI
277	NG = L_NGAUSS
278	DO L=1,L_NLAYRAD-1
279
280	K = 2*L+1
281	DTAUI(L,nw,ng) = DTAUKI(K,NW,NG) + DTAUKI(K+1,NW,NG)! + 1.e-50
282
283	atemp = 0.
284	if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
285	atemp = atemp + &
286	ASF_T(K,NW)*DHAZES_T(K,NW) + &
287	ASF_T(K+1,NW)*DHAZES_T(K+1,NW)
288
289	WBARI(L,nw,ng) = btemp(L,nw) / DTAUI(L,NW,NG)
290	else
291	WBARI(L,nw,ng) = 0.0D0
292	DTAUI(L,NW,NG) = 1.0D-9
293	endif
294
295	if(btemp(L,nw) .GT. 0.0d0) then
296	cosbi(L,NW,NG) = atemp/btemp(L,nw)
297	else
298	cosbi(L,NW,NG) = 0.0D0
299	end if
300
301	END DO ! L vertical loop
302
303	! Last level
304
305	L = L_NLAYRAD
306	K = 2*L+1
307	DTAUI(L,nw,ng) = DTAUKI(K,NW,NG) ! + 1.e-50
308
309	atemp = 0.
310	if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
311	atemp = atemp + ASF_T(K,NW)*DHAZES_T(K,NW)
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
325	! Now the other Gauss points, if needed.
326
327	DO NG=1,L_NGAUSS-1
328	IF(TAUGSURF(NW,NG) .gt. TLIMIT) THEN
329
330	DO L=1,L_NLAYRAD-1
331	K = 2*L+1
332	DTAUI(L,nw,ng) = DTAUKI(K,NW,NG)+DTAUKI(K+1,NW,NG)! + 1.e-50
333
334	if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
335
336	WBARI(L,nw,ng) = btemp(L,nw) / DTAUI(L,NW,NG)
337
338	else
339	WBARI(L,nw,ng) = 0.0D0
340	DTAUI(L,NW,NG) = 1.0D-9
341	endif
342
343	cosbi(L,NW,NG) = cosbi(L,NW,L_NGAUSS)
344	END DO ! L vertical loop
345
346	! Last level
347	L = L_NLAYRAD
348	K = 2*L+1
349	DTAUI(L,nw,ng) = DTAUKI(K,NW,NG)! + 1.e-50
350
351	if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
352
353	WBARI(L,nw,ng) = btemp(L,nw) / DTAUI(L,NW,NG)
354
355	else
356	WBARI(L,nw,ng) = 0.0D0
357	DTAUI(L,NW,NG) = 1.0D-9
358	endif
359
360	cosbi(L,NW,NG) = cosbi(L,NW,L_NGAUSS)
361
362	END IF
363
364	END DO ! NG Gauss loop
365	END DO ! NW spectral loop
366
367	! Total extinction optical depths
368
369	DO NG=1,L_NGAUSS ! full gauss loop
370	DO NW=1,L_NSPECTI
371	TAUCUMI(1,NW,NG)=0.0D0
372	DO K=2,L_LEVELS
373	TAUCUMI(K,NW,NG)=TAUCUMI(K-1,NW,NG)+DTAUKI(K,NW,NG)
374	END DO
375	END DO ! end full gauss loop
376	END DO
377
378	! be aware when comparing with textbook results
379	! (e.g. Pierrehumbert p. 218) that
380	! taucumi does not take the <cos theta>=0.5 factor into
381	! account. It is the optical depth for a vertically
382	! ascending ray with angle theta = 0.
383
384
385	! Titan's outputs (J.V.O, 2016)===============================================
386	! do l=1,L_NLAYRAD
387	! do nw=1,L_NSPECTI
388	! INT_DTAU(L,NW) = 0.0d+0
389	! DO NG=1,L_NGAUSS
390	! INT_DTAU(L,NW)= INT_DTAU(L,NW) + dtaui(L,nw,ng)*gweight(NG)
391	! enddo
392	! enddo
393	! enddo
394
395	! do nw=1,L_NSPECTI
396	! write(str2,'(i2.2)') nw
397	! 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))
398	! 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))
399	! enddo
400
401	! ==============================================================================
402
403	return
404
405
406	end subroutine optci
407
408
409

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

Download in other formats:

Original Format