Context Navigation

optci.F90 @ 2006

Last change on this file since 2006 was 1947, checked in by jvatant, 8 years ago
Enables altitude-depending gravity field g(z) (glat->gzlat) in physics + Can be dangerous ( disagreement with dyn) but important (compulsive !) to have correct altitudes in the chemistry + Can be activated with eff_gz flag in callphys.def -- JVO
File size: 13.6 KB

Line
1	subroutine optci(PQO,NLAY,PLEV,TLEV,TMID,PMID, &
2	DTAUI,TAUCUMI,COSBI,WBARI,TAUGSURF)
3
4	use radinc_h
5	use radcommon_h, only: gasi,tlimit,Cmk,gzlat_ig,tgasref,pfgasref,wnoi,scalep,indi,gweight
6	use gases_h
7	use comcstfi_mod, only: r
8	use callkeys_mod, only: continuum,graybody,callclouds,callmufi, uncoupl_optic_haze
9	use tracer_h, only : nmicro,nice
10	use MMP_OPTICS
11
12	implicit none
13
14	!==================================================================
15	!
16	! Purpose
17	! -------
18	! Calculates longwave optical constants at each level. For each
19	! layer and spectral interval in the IR it calculates WBAR, DTAU
20	! and COSBAR. For each level it calculates TAU.
21	!
22	! TAUCUMI(L,LW) is the cumulative optical depth at level L (or alternatively
23	! at the bottom of layer L), LW is the spectral wavelength interval.
24	!
25	! TLEV(L) - Temperature at the layer boundary (i.e., level)
26	! PLEV(L) - Pressure at the layer boundary (i.e., level)
27	!
28	! Authors
29	! -------
30	! Adapted from the NASA Ames code by R. Wordsworth (2009)
31	! Clean and adaptation to Titan by J. Vatant d'Ollone (2016-17)
32	!
33	!==================================================================
34
35
36	!==========================================================
37	! Input/Output
38	!==========================================================
39	REAL*8, INTENT(IN) :: PQO(nlay,nmicro) ! Tracers (X/m2).
40	INTEGER, INTENT(IN) :: NLAY ! Number of pressure layers (for pqo)
41	REAL*8, INTENT(IN) :: PLEV(L_LEVELS), TLEV(L_LEVELS)
42	REAL*8, INTENT(IN) :: TMID(L_LEVELS), PMID(L_LEVELS)
43
44	REAL*8, INTENT(OUT) :: DTAUI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
45	REAL*8, INTENT(OUT) :: TAUCUMI(L_LEVELS,L_NSPECTI,L_NGAUSS)
46	REAL*8, INTENT(OUT) :: COSBI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
47	REAL*8, INTENT(OUT) :: WBARI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
48	REAL*8, INTENT(OUT) :: TAUGSURF(L_NSPECTI,L_NGAUSS-1)
49	! ==========================================================
50
51	real*8 DTAUKI(L_LEVELS,L_NSPECTI,L_NGAUSS)
52
53	! Titan customisation
54	! J. Vatant d'Ollone (2016)
55	real*8 DHAZE_T(L_LEVELS,L_NSPECTI)
56	real*8 DHAZES_T(L_LEVELS,L_NSPECTI)
57	real*8 SSA_T(L_LEVELS,L_NSPECTI)
58	real*8 ASF_T(L_LEVELS,L_NSPECTI)
59	real*8 INT_DTAU(L_NLAYRAD,L_NSPECTI)
60	real*8 K_HAZE(L_NLAYRAD,L_NSPECTI)
61
62	CHARACTER*2 str2
63	! ==========================
64
65	integer L, NW, NG, K, LK, IAER
66	integer MT(L_LEVELS), MP(L_LEVELS), NP(L_LEVELS)
67	real*8 ANS, TAUGAS
68	real*8 DPR(L_LEVELS), U(L_LEVELS)
69	real*8 LCOEF(4), LKCOEF(L_LEVELS,4)
70
71	real*8 DCONT
72	double precision wn_cont, p_cont, p_air, T_cont, dtemp, dtempc
73	double precision p_cross
74
75	real*8 KCOEF(4)
76
77	! temporary variable to reduce memory access time to gasi
78	real*8 tmpk(2,2)
79
80	! temporary variables for multiple aerosol calculation
81	real*8 atemp
82	real*8 btemp(L_NLAYRAD,L_NSPECTI)
83
84	! variables for k in units m^-1
85	real*8 dz(L_LEVELS)
86	!real*8 rho !! see test below
87
88	integer igas, jgas, ilay
89
90	integer interm
91
92	real*8 m0as,m3as,m0af,m3af
93	real*8 ext_s,sca_s,ssa_s,asf_s
94	real*8 ext_f,sca_f,ssa_f,asf_f
95	logical,save :: firstcall=.true.
96	!$OMP THREADPRIVATE(firstcall)
97
98	!! AS: to save time in computing continuum (see bilinearbig)
99	IF (.not.ALLOCATED(indi)) THEN
100	ALLOCATE(indi(L_NSPECTI,ngasmx,ngasmx))
101	indi = -9999 ! this initial value means "to be calculated"
102	ENDIF
103
104	! Some initialisation beacause there's a pb with disr_haze at the limits (nw=1)
105	! I should check this - For now we set vars to zero : better than nans - JVO 2017
106
107	dhaze_t(:,:) = 0.
108	ssa_t(:,:) = 0.
109	asf_t(:,:) = 0.
110
111	!=======================================================================
112	! Determine the total gas opacity throughout the column, for each
113	! spectral interval, NW, and each Gauss point, NG.
114
115	taugsurf(:,:) = 0.0
116	dpr(:) = 0.0
117	lkcoef(:,:) = 0.0
118
119	do K=2,L_LEVELS
120
121	ilay = k / 2 ! int. arithmetic => gives the gcm layer index
122
123	DPR(k) = PLEV(K)-PLEV(K-1)
124
125	! if we have continuum opacities, we need dz
126	dz(k) = dpr(k)RTMID(K)/(gzlat_ig(ilay)*PMID(K))
127	U(k) = Cmk(ilay)*DPR(k) ! only Cmk line in optci.F
128
129	call tpindex(PMID(K),TMID(K),pfgasref,tgasref,LCOEF,MT(K),MP(K))
130
131	do LK=1,4
132	LKCOEF(K,LK) = LCOEF(LK)
133	end do
134	end do ! levels
135
136	do NW=1,L_NSPECTI
137
138	do K=2,L_LEVELS
139
140	ilay = k / 2 ! int. arithmetic => gives the gcm layer index
141
142	! Optical coupling of YAMMS is plugged but inactivated for now
143	! as long as the microphysics only isn't fully debugged -- JVO 01/18
144	IF (callmufi .AND. (.NOT. uncoupl_optic_haze)) THEN
145	m0as = pqo(ilay,1)
146	m3as = pqo(ilay,2)
147	m0af = pqo(ilay,3)
148	m3af = pqo(ilay,4)
149
150	IF (.NOT.mmp_sph_optics_ir(m0as,m3as,nw,ext_s,sca_s,ssa_s,asf_s)) &
151	CALL abort_gcm("optcv", "Fatal error in mmp_sph_optics_ir", 12)
152	IF (.NOT.mmp_fra_optics_ir(m0af,m3af,nw,ext_f,sca_f,ssa_f,asf_f)) &
153	CALL abort_gcm("optcv", "Fatal error in mmp_fra_optics_ir", 12)
154	dhaze_T(k,nw) = ext_s+ext_f
155	SSA_T(k,nw) = (sca_s+sca_f)/dhaze_T(k,nw)
156	ASF_T(k,nw) = (asf_ssca_s + asf_fsca_f) /(sca_s+sca_f)
157	IF (callclouds.and.firstcall) &
158	WRITE(,) 'WARNING: In optci, optical properties &
159	&calculations are not implemented yet'
160	ELSE
161	! Call fixed vertical haze profile of extinction - same for all columns
162	call disr_haze(dz(k),plev(k),wnoi(nw),dhaze_T(k,nw),SSA_T(k,nw),ASF_T(k,nw))
163	ENDIF
164
165	DCONT = 0.0d0 ! continuum absorption
166
167	if(continuum.and.(.not.graybody))then
168	! include continua if necessary
169	wn_cont = dble(wnoi(nw))
170	T_cont = dble(TMID(k))
171	do igas=1,ngasmx
172
173	p_cont = dble(PMID(k)scalepgfrac(igas,ilay))
174
175	dtemp=0.0d0
176	if(igas.eq.igas_N2)then
177
178	interm = indi(nw,igas,igas)
179	call interpolateN2N2(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
180	indi(nw,igas,igas) = interm
181
182	elseif(igas.eq.igas_H2)then
183
184	! first do self-induced absorption
185	interm = indi(nw,igas,igas)
186	call interpolateH2H2(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
187	indi(nw,igas,igas) = interm
188
189	! then cross-interactions with other gases
190	do jgas=1,ngasmx
191	p_cross = dble(PMID(k)scalepgfrac(jgas,ilay))
192	dtempc = 0.0d0
193	if(jgas.eq.igas_N2)then
194	interm = indi(nw,igas,jgas)
195	call interpolateN2H2(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
196	indi(nw,igas,jgas) = interm
197	endif
198	dtemp = dtemp + dtempc
199	enddo
200
201	elseif(igas.eq.igas_CH4)then
202
203	! first do self-induced absorption
204	interm = indi(nw,igas,igas)
205	call interpolateCH4CH4(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
206	indi(nw,igas,igas) = interm
207
208	! then cross-interactions with other gases
209	do jgas=1,ngasmx
210	p_cross = dble(PMID(k)scalepgfrac(jgas,ilay))
211	dtempc = 0.0d0
212	if(jgas.eq.igas_N2)then
213	interm = indi(nw,igas,jgas)
214	call interpolateN2CH4(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
215	indi(nw,igas,jgas) = interm
216	endif
217	dtemp = dtemp + dtempc
218	enddo
219
220	endif
221
222	DCONT = DCONT + dtemp
223
224	enddo
225
226	DCONT = DCONT*dz(k)
227
228	endif
229
230	do ng=1,L_NGAUSS-1
231
232	! Now compute TAUGAS
233
234	! JVO 2017 : added tmpk because the repeated calls to gasi/v increased dramatically
235	! the execution time of optci/v -> ~ factor 2 on the whole radiative
236	! transfer on the tested simulations !
237
238	tmpk = GASI(MT(K):MT(K)+1,MP(K):MP(K)+1,1,NW,NG)
239
240	KCOEF(1) = tmpk(1,1) ! KCOEF(1) = GASI(MT(K),MP(K),1,NW,NG)
241	KCOEF(2) = tmpk(1,2) ! KCOEF(2) = GASI(MT(K),MP(K)+1,1,NW,NG)
242	KCOEF(3) = tmpk(2,2) ! KCOEF(3) = GASI(MT(K)+1,MP(K)+1,1,NW,NG)
243	KCOEF(4) = tmpk(2,1) ! KCOEF(4) = GASI(MT(K)+1,MP(K),1,NW,NG)
244
245
246	! Interpolate the gaseous k-coefficients to the requested T,P values
247
248	ANS = LKCOEF(K,1)KCOEF(1) + LKCOEF(K,2)KCOEF(2) + &
249	LKCOEF(K,3)KCOEF(3) + LKCOEF(K,4)KCOEF(4)
250
251	TAUGAS = U(k)*ANS
252
253	TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS + DCONT
254	DTAUKI(K,nw,ng) = TAUGAS &
255	+ DCONT & ! For parameterized continuum absorption
256	+ DHAZE_T(K,NW) ! For Titan haze
257
258	end do
259
260	! Now fill in the "clear" part of the spectrum (NG = L_NGAUSS),
261	! which holds continuum opacity only
262
263	NG = L_NGAUSS
264	DTAUKI(K,nw,ng) = 0.d0 &
265	+ DCONT & ! For parameterized continuum absorption
266	+ DHAZE_T(K,NW) ! For Titan Haze
267
268	end do
269	end do
270
271	!=======================================================================
272	! Now the full treatment for the layers, where besides the opacity
273	! we need to calculate the scattering albedo and asymmetry factors
274	! ======================================================================
275
276	! Haze scattering
277	DO NW=1,L_NSPECTI
278	DO K=2,L_LEVELS
279	DHAZES_T(K,NW) = DHAZE_T(K,NW) * SSA_T(K,NW)
280	ENDDO
281	ENDDO
282
283	DO NW=1,L_NSPECTI
284	DO L=1,L_NLAYRAD-1
285	K = 2*L+1
286	btemp(L,NW) = DHAZES_T(K,NW) + DHAZES_T(K+1,NW)
287	END DO ! L vertical loop
288
289	! Last level
290	L = L_NLAYRAD
291	K = 2*L+1
292	btemp(L,NW) = DHAZES_T(K,NW)
293
294	END DO ! NW spectral loop
295
296
297	DO NW=1,L_NSPECTI
298	NG = L_NGAUSS
299	DO L=1,L_NLAYRAD-1
300
301	K = 2*L+1
302	DTAUI(L,nw,ng) = DTAUKI(K,NW,NG) + DTAUKI(K+1,NW,NG)! + 1.e-50
303
304	atemp = 0.
305	if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
306	atemp = atemp + &
307	ASF_T(K,NW)*DHAZES_T(K,NW) + &
308	ASF_T(K+1,NW)*DHAZES_T(K+1,NW)
309
310	WBARI(L,nw,ng) = btemp(L,nw) / DTAUI(L,NW,NG)
311	else
312	WBARI(L,nw,ng) = 0.0D0
313	DTAUI(L,NW,NG) = 1.0D-9
314	endif
315
316	if(btemp(L,nw) .GT. 0.0d0) then
317	cosbi(L,NW,NG) = atemp/btemp(L,nw)
318	else
319	cosbi(L,NW,NG) = 0.0D0
320	end if
321
322	END DO ! L vertical loop
323
324	! Last level
325
326	L = L_NLAYRAD
327	K = 2*L+1
328	DTAUI(L,nw,ng) = DTAUKI(K,NW,NG) ! + 1.e-50
329
330	atemp = 0.
331	if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
332	atemp = atemp + ASF_T(K,NW)*DHAZES_T(K,NW)
333	WBARI(L,nw,ng) = btemp(L,nw) / DTAUI(L,NW,NG)
334	else
335	WBARI(L,nw,ng) = 0.0D0
336	DTAUI(L,NW,NG) = 1.0D-9
337	endif
338
339	if(btemp(L,nw) .GT. 0.0d0) then
340	cosbi(L,NW,NG) = atemp/btemp(L,nw)
341	else
342	cosbi(L,NW,NG) = 0.0D0
343	end if
344
345
346	! Now the other Gauss points, if needed.
347
348	DO NG=1,L_NGAUSS-1
349	IF(TAUGSURF(NW,NG) .gt. TLIMIT) THEN
350
351	DO L=1,L_NLAYRAD-1
352	K = 2*L+1
353	DTAUI(L,nw,ng) = DTAUKI(K,NW,NG)+DTAUKI(K+1,NW,NG)! + 1.e-50
354
355	if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
356
357	WBARI(L,nw,ng) = btemp(L,nw) / DTAUI(L,NW,NG)
358
359	else
360	WBARI(L,nw,ng) = 0.0D0
361	DTAUI(L,NW,NG) = 1.0D-9
362	endif
363
364	cosbi(L,NW,NG) = cosbi(L,NW,L_NGAUSS)
365	END DO ! L vertical loop
366
367	! Last level
368	L = L_NLAYRAD
369	K = 2*L+1
370	DTAUI(L,nw,ng) = DTAUKI(K,NW,NG)! + 1.e-50
371
372	if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
373
374	WBARI(L,nw,ng) = btemp(L,nw) / DTAUI(L,NW,NG)
375
376	else
377	WBARI(L,nw,ng) = 0.0D0
378	DTAUI(L,NW,NG) = 1.0D-9
379	endif
380
381	cosbi(L,NW,NG) = cosbi(L,NW,L_NGAUSS)
382
383	END IF
384
385	END DO ! NG Gauss loop
386	END DO ! NW spectral loop
387
388	! Total extinction optical depths
389
390	DO NG=1,L_NGAUSS ! full gauss loop
391	DO NW=1,L_NSPECTI
392	TAUCUMI(1,NW,NG)=0.0D0
393	DO K=2,L_LEVELS
394	TAUCUMI(K,NW,NG)=TAUCUMI(K-1,NW,NG)+DTAUKI(K,NW,NG)
395	END DO
396	END DO ! end full gauss loop
397	END DO
398
399	! be aware when comparing with textbook results
400	! (e.g. Pierrehumbert p. 218) that
401	! taucumi does not take the <cos theta>=0.5 factor into
402	! account. It is the optical depth for a vertically
403	! ascending ray with angle theta = 0.
404
405
406	! Titan's outputs (J.V.O, 2016)===============================================
407	! do l=1,L_NLAYRAD
408	! do nw=1,L_NSPECTI
409	! INT_DTAU(L,NW) = 0.0d+0
410	! DO NG=1,L_NGAUSS
411	! INT_DTAU(L,NW)= INT_DTAU(L,NW) + dtaui(L,nw,ng)*gweight(NG)
412	! enddo
413	! enddo
414	! enddo
415
416	! do nw=1,L_NSPECTI
417	! write(str2,'(i2.2)') nw
418	! 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))
419	! 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))
420	! enddo
421
422	! ==============================================================================
423
424	if(firstcall) firstcall = .false.
425
426	return
427
428
429	end subroutine optci
430
431
432

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

Download in other formats:

Original Format