Context Navigation

optci.F90 @ 2239

Last change on this file since 2239 was 2133, checked in by jvatant, 6 years ago
Fic a bug in r2131 the writediagfi cannot be called in RT routines because of the iradia it goes crazy with outputs freq ... --JVO
File size: 13.2 KB

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

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

Download in other formats:

Original Format