Context Navigation

optci.F90 @ 2870

Last change on this file since 2870 was 2861, checked in by mturbet, 3 years ago
add the CO2-CH4 CIA opacity from Turbet2020 in the model
File size: 16.5 KB

Line
1	MODULE optci_mod
2
3	IMPLICIT NONE
4
5	CONTAINS
6
7	subroutine optci(PLEV,TLEV,DTAUI,TAUCUMI, &
8	QXIAER,QSIAER,GIAER,COSBI,WBARI,TAUAERO, &
9	TMID,PMID,TAUGSURF,QVAR,MUVAR)
10
11	use radinc_h, only: L_LEVELS, L_NLAYRAD, L_NSPECTI, L_NGAUSS, &
12	L_NLEVRAD, L_REFVAR, naerkind
13	use radcommon_h, only: gasi,tlimit,wrefVAR,Cmk,tgasref,pfgasref,wnoi,scalep,indi,glat_ig
14	use gases_h, only: gfrac, ngasmx, igas_N2, igas_He, igas_H2O, igas_H2, igas_CH4
15	use comcstfi_mod, only: g, r, mugaz
16	use callkeys_mod, only: kastprof,continuum,graybody
17	use recombin_corrk_mod, only: corrk_recombin, gasi_recomb
18	use tpindex_mod, only: tpindex
19
20	implicit none
21
22	!==================================================================
23	!
24	! Purpose
25	! -------
26	! Calculates longwave optical constants at each level. For each
27	! layer and spectral interval in the IR it calculates WBAR, DTAU
28	! and COSBAR. For each level it calculates TAU.
29	!
30	! TAUCUMI(L,LW) is the cumulative optical depth at level L (or alternatively
31	! at the bottom of layer L), LW is the spectral wavelength interval.
32	!
33	! TLEV(L) - Temperature at the layer boundary (i.e., level)
34	! PLEV(L) - Pressure at the layer boundary (i.e., level)
35	!
36	! Authors
37	! -------
38	! Adapted from the NASA Ames code by R. Wordsworth (2009)
39	!
40	!==================================================================
41
42
43	real*8 DTAUI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
44	real*8 DTAUKI(L_LEVELS,L_NSPECTI,L_NGAUSS)
45	real*8 TAUI(L_NLEVRAD,L_NSPECTI,L_NGAUSS)
46	real*8 TAUCUMI(L_LEVELS,L_NSPECTI,L_NGAUSS)
47	real*8 PLEV(L_LEVELS)
48	real*8 TLEV(L_LEVELS)
49	real*8 TMID(L_LEVELS), PMID(L_LEVELS)
50	real*8 COSBI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
51	real*8 WBARI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
52
53	! for aerosols
54	real*8 QXIAER(L_LEVELS,L_NSPECTI,NAERKIND)
55	real*8 QSIAER(L_LEVELS,L_NSPECTI,NAERKIND)
56	real*8 GIAER(L_LEVELS,L_NSPECTI,NAERKIND)
57	real*8 TAUAERO(L_LEVELS,NAERKIND)
58	real*8 TAUAEROLK(L_LEVELS,L_NSPECTI,NAERKIND)
59	real*8 TAEROS(L_LEVELS,L_NSPECTI,NAERKIND)
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 taugsurf(L_NSPECTI,L_NGAUSS-1)
68	real*8 DCONT,DAERO
69	double precision wn_cont, p_cont, p_air, T_cont, dtemp, dtempc
70	double precision p_cross
71
72	! variable species mixing ratio variables
73	real*8 QVAR(L_LEVELS), WRATIO(L_LEVELS), MUVAR(L_LEVELS)
74	real*8 KCOEF(4)
75	integer NVAR(L_LEVELS)
76
77	! temporary variables to reduce memory access time to gasi
78	real*8 tmpk(2,2)
79	real*8 tmpkvar(2,2,2)
80
81	! temporary variables for multiple aerosol calculation
82	real*8 atemp
83	real*8 btemp(L_NLAYRAD,L_NSPECTI)
84
85	! variables for k in units m^-1
86	real*8 dz(L_LEVELS)
87	!real*8 rho !! see test below
88
89	integer igas, jgas
90
91	integer interm
92
93	!--- Kasting's CIA ----------------------------------------
94	!real*8, parameter :: Ci(L_NSPECTI)=[ &
95	! 3.8E-5, 1.2E-5, 2.8E-6, 7.6E-7, 4.5E-7, 2.3E-7, &
96	! 5.4E-7, 1.6E-6, 0.0, &
97	! 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
98	! 0.0, 4.0E-7, 4.0E-6, 1.4E-5, &
99	! 1.0E-5, 1.2E-6, 2.0E-7, 5.0E-8, 3.0E-8, 0.0 ]
100	!real*8, parameter :: Ti(L_NSPECTI)=[ -2.2, -1.9, &
101	! -1.7, -1.7, -1.7, -1.7, -1.7, -1.7, &
102	! 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
103	! -1.7,-1.7,-1.7,-1.7,-1.7,-1.7,-1.7, -1.7,0.0 ]
104	!----------------------------------------------------------
105
106	!! AS: to save time in computing continuum (see bilinearbig)
107	IF (.not.ALLOCATED(indi)) THEN
108	ALLOCATE(indi(L_NSPECTI,ngasmx,ngasmx))
109	indi = -9999 ! this initial value means "to be calculated"
110	ENDIF
111
112	!=======================================================================
113	! Determine the total gas opacity throughout the column, for each
114	! spectral interval, NW, and each Gauss point, NG.
115
116	taugsurf(:,:) = 0.0
117	dpr(:) = 0.0
118	lkcoef(:,:) = 0.0
119
120	do K=2,L_LEVELS
121	DPR(k) = PLEV(K)-PLEV(K-1)
122
123	!--- Kasting's CIA ----------------------------------------
124	!dz(k)=dpr(k)189.02TMID(K)/(0.03720*PMID(K))
125	! this is CO2 path length (in cm) as written by Francois
126	! delta_z = delta_p * R_specific * T / (g * P)
127	! But Kasting states that W is in units of _atmosphere_ cm
128	! So we do
129	!dz(k)=dz(k)*(PMID(K)/1013.25)
130	!dz(k)=dz(k)/100.0 ! in m for SI calc
131	!----------------------------------------------------------
132
133	! if we have continuum opacities, we need dz
134	if(kastprof)then
135	dz(k) = dpr(k)(1000.0d08.3145d0/muvar(k))TMID(K)/(gPMID(K))
136	U(k) = CmkDPR(k)mugaz/muvar(k)
137	else
138	dz(k) = dpr(k)RTMID(K)/(glat_igPMID(K))mugaz/muvar(k)
139	U(k) = CmkDPR(k)mugaz/muvar(k) ! only Cmk line in optci.F
140	!JL13 the mugaz/muvar factor takes into account water meanmolecular weight if water is present
141	endif
142
143	call tpindex(PMID(K),TMID(K),QVAR(K),pfgasref,tgasref,WREFVAR, &
144	LCOEF,MT(K),MP(K),NVAR(K),WRATIO(K))
145
146	do LK=1,4
147	LKCOEF(K,LK) = LCOEF(LK)
148	end do
149	end do ! levels
150
151	! Spectral dependance of aerosol absorption
152	do iaer=1,naerkind
153	DO NW=1,L_NSPECTI
154	do K=2,L_LEVELS
155	TAEROS(K,NW,IAER) = TAUAERO(K,IAER) * QXIAER(K,NW,IAER)
156	end do ! levels
157	END DO
158	end do
159
160	do NW=1,L_NSPECTI
161
162	do K=2,L_LEVELS
163
164	DAERO=SUM(TAEROS(K,NW,1:naerkind)) ! aerosol absorption
165
166	DCONT = 0.0d0 ! continuum absorption
167
168	if(continuum.and.(.not.graybody))then
169	! include continua if necessary
170	wn_cont = dble(wnoi(nw))
171	T_cont = dble(TMID(k))
172	do igas=1,ngasmx
173
174	if(gfrac(igas).eq.-1)then ! variable
175	p_cont = dble(PMID(k)scalepQVAR(k)) ! qvar = mol/mol
176	else
177	p_cont = dble(PMID(k)scalepgfrac(igas)*(1.-QVAR(k)))
178	endif
179
180	dtemp=0.0d0
181	if(igas.eq.igas_N2)then
182
183	interm = indi(nw,igas,igas)
184	call interpolateN2N2(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
185	indi(nw,igas,igas) = interm
186
187	elseif(igas.eq.igas_H2)then
188
189	! first do self-induced absorption
190	interm = indi(nw,igas,igas)
191	call interpolateH2H2(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
192	indi(nw,igas,igas) = interm
193
194	! then cross-interactions with other gases
195	do jgas=1,ngasmx
196	p_cross = dble(PMID(k)scalepgfrac(jgas)*(1.-QVAR(k)))
197	dtempc = 0.0d0
198	if(jgas.eq.igas_N2)then
199	interm = indi(nw,igas,jgas)
200	call interpolateN2H2(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
201	indi(nw,igas,jgas) = interm
202	elseif(jgas.eq.igas_CO2)then
203	interm = indi(nw,igas,jgas)
204	call interpolateCO2H2(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
205	indi(nw,igas,jgas) = interm
206	elseif(jgas.eq.igas_He)then
207	interm = indi(nw,igas,jgas)
208	call interpolateH2He(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
209	indi(nw,igas,jgas) = interm
210	endif
211	dtemp = dtemp + dtempc
212	enddo
213
214	elseif(igas.eq.igas_CH4)then
215
216	! first do self-induced absorption
217	interm = indi(nw,igas,igas)
218	call interpolateCH4CH4(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
219	indi(nw,igas,igas) = interm
220
221	! then cross-interactions with other gases
222	do jgas=1,ngasmx
223	p_cross = dble(PMID(k)scalepgfrac(jgas)*(1.-QVAR(k)))
224	dtempc = 0.0d0
225	if(jgas.eq.igas_H2)then
226	interm = indi(nw,igas,jgas)
227	call interpolateH2CH4(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
228	indi(nw,igas,jgas) = interm
229	elseif(jgas.eq.igas_CO2)then
230	interm = indi(nw,igas,jgas)
231	call interpolateCO2CH4(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
232	indi(nw,igas,jgas) = interm
233	elseif(jgas.eq.igas_He)then
234	interm = indi(nw,igas,jgas)
235	call interpolateHeCH4(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
236	indi(nw,igas,jgas) = interm
237	endif
238	dtemp = dtemp + dtempc
239	enddo
240
241	elseif(igas.eq.igas_H2O.and.T_cont.gt.100.0)then
242	! Compute self and foreign (with air) continuum of H2O
243	p_air = dble(PMID(k)*scalep) - p_cont ! note assumes background is air!
244	interm = indi(nw,igas,igas)
245	call interpolateH2O_self_foreign(wn_cont,T_cont,p_cont,p_air,dtemp,.false.,interm) ! MTCKD v3.3
246	indi(nw,igas,igas) = interm
247
248	endif
249
250	DCONT = DCONT + dtemp
251
252	enddo
253
254	! Oobleck test
255	!rho = PMID(k)scalep / (TMID(k)286.99)
256	!if(WNOI(nw).gt.300.0 .and. WNOI(nw).lt.500.0)then
257	! DCONT = rho * 0.125 * 4.6e-4
258	!elseif(WNOI(nw).gt.500.0 .and. WNOI(nw).lt.700.0)then
259	! DCONT = 1000dpr(k) 1.0 * 4.6e-4 / g
260	! DCONT = rho * 1.0 * 4.6e-4
261	!elseif(WNOI(nw).gt.700.0 .and. WNOI(nw).lt.900.0)then
262	! DCONT = rho * 0.125 * 4.6e-4
263	!endif
264
265	DCONT = DCONT*dz(k)
266
267	endif
268
269	do ng=1,L_NGAUSS-1
270
271	! Now compute TAUGAS
272
273	! Interpolate between water mixing ratios
274	! WRATIO = 0.0 if the requested water amount is equal to, or outside the
275	! the water data range
276
277	if(L_REFVAR.eq.1)then ! added by RW for special no variable case
278
279	! JVO 2017 : added tmpk because the repeated calls to gasi/v increased dramatically
280	! the execution time of optci/v -> ~ factor 2 on the whole radiative
281	! transfer on the tested simulations !
282
283	IF (corrk_recombin) THEN ! added by JVO
284	tmpk = GASI_RECOMB(MT(K):MT(K)+1,MP(K):MP(K)+1,1,NW,NG) ! contains the mix of recombined species
285	ELSE
286	tmpk = GASI(MT(K):MT(K)+1,MP(K):MP(K)+1,1,NW,NG)
287	ENDIF
288
289	KCOEF(1) = tmpk(1,1) ! KCOEF(1) = GASI(MT(K),MP(K),1,NW,NG)
290	KCOEF(2) = tmpk(1,2) ! KCOEF(2) = GASI(MT(K),MP(K)+1,1,NW,NG)
291	KCOEF(3) = tmpk(2,2) ! KCOEF(3) = GASI(MT(K)+1,MP(K)+1,1,NW,NG)
292	KCOEF(4) = tmpk(2,1) ! KCOEF(4) = GASI(MT(K)+1,MP(K),1,NW,NG)
293
294	else
295
296	IF (corrk_recombin) THEN ! added by JVO
297	tmpkvar = GASI_RECOMB(MT(K):MT(K)+1,MP(K):MP(K)+1,NVAR(K):NVAR(K)+1,NW,NG)
298	ELSE
299	tmpkvar = GASI(MT(K):MT(K)+1,MP(K):MP(K)+1,NVAR(K):NVAR(K)+1,NW,NG)
300	ENDIF
301
302	KCOEF(1) = tmpkvar(1,1,1) + WRATIO(K) * &
303	( tmpkvar(1,1,2)-tmpkvar(1,1,1) )
304
305	KCOEF(2) = tmpkvar(1,2,1) + WRATIO(K) * &
306	( tmpkvar(1,2,2)-tmpkvar(1,2,1) )
307
308	KCOEF(3) = tmpkvar(2,2,1) + WRATIO(K) * &
309	( tmpkvar(2,2,2)-tmpkvar(2,2,1) )
310
311	KCOEF(4) = tmpkvar(2,1,1) + WRATIO(K) * &
312	( tmpkvar(2,1,2)-tmpkvar(2,1,1) )
313
314	endif
315
316	! Interpolate the gaseous k-coefficients to the requested T,P values
317
318	ANS = LKCOEF(K,1)KCOEF(1) + LKCOEF(K,2)KCOEF(2) + &
319	LKCOEF(K,3)KCOEF(3) + LKCOEF(K,4)KCOEF(4)
320
321	TAUGAS = U(k)*ANS
322
323	TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS + DCONT
324	DTAUKI(K,nw,ng) = TAUGAS &
325	+ DCONT & ! For parameterized continuum absorption
326	+ DAERO ! For aerosol absorption
327
328	end do
329
330	! Now fill in the "clear" part of the spectrum (NG = L_NGAUSS),
331	! which holds continuum opacity only
332
333	NG = L_NGAUSS
334	DTAUKI(K,nw,ng) = 0.d0 &
335	+ DCONT & ! For parameterized continuum absorption
336	+ DAERO ! For aerosol absorption
337
338	end do
339	end do
340
341	!=======================================================================
342	! Now the full treatment for the layers, where besides the opacity
343	! we need to calculate the scattering albedo and asymmetry factors
344
345	do iaer=1,naerkind
346	DO NW=1,L_NSPECTI
347	DO K=2,L_LEVELS
348	TAUAEROLK(K,NW,IAER) = TAUAERO(K,IAER)*QSIAER(K,NW,IAER) ! effect of scattering albedo
349	ENDDO
350	ENDDO
351	end do
352
353	DO NW=1,L_NSPECTI
354	DO L=1,L_NLAYRAD-1
355	K = 2*L+1
356	btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind)) + SUM(TAUAEROLK(K+1,NW,1:naerkind))
357	END DO ! L vertical loop
358
359	! Last level
360	L = L_NLAYRAD
361	K = 2*L+1
362	btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind))
363
364	END DO ! NW spectral loop
365
366
367	DO NW=1,L_NSPECTI
368	NG = L_NGAUSS
369	DO L=1,L_NLAYRAD-1
370
371	K = 2*L+1
372	DTAUI(L,nw,ng) = DTAUKI(K,NW,NG) + DTAUKI(K+1,NW,NG)! + 1.e-50
373
374	atemp = 0.
375	if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
376	do iaer=1,naerkind
377	atemp = atemp + &
378	GIAER(K,NW,IAER) * TAUAEROLK(K,NW,IAER) + &
379	GIAER(K+1,NW,IAER) * TAUAEROLK(K+1,NW,IAER)
380	end do
381	WBARI(L,nw,ng) = btemp(L,nw) / DTAUI(L,NW,NG)
382	else
383	WBARI(L,nw,ng) = 0.0D0
384	DTAUI(L,NW,NG) = 1.0D-9
385	endif
386
387	if(btemp(L,nw) .GT. 0.0d0) then
388	cosbi(L,NW,NG) = atemp/btemp(L,nw)
389	else
390	cosbi(L,NW,NG) = 0.0D0
391	end if
392
393	END DO ! L vertical loop
394
395	! Last level
396
397	L = L_NLAYRAD
398	K = 2*L+1
399	DTAUI(L,nw,ng) = DTAUKI(K,NW,NG) ! + 1.e-50
400
401	atemp = 0.
402	if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
403	do iaer=1,naerkind
404	atemp = atemp + GIAER(K,NW,IAER) * TAUAEROLK(K,NW,IAER)
405	end do
406	WBARI(L,nw,ng) = btemp(L,nw) / DTAUI(L,NW,NG)
407	else
408	WBARI(L,nw,ng) = 0.0D0
409	DTAUI(L,NW,NG) = 1.0D-9
410	endif
411
412	if(btemp(L,nw) .GT. 0.0d0) then
413	cosbi(L,NW,NG) = atemp/btemp(L,nw)
414	else
415	cosbi(L,NW,NG) = 0.0D0
416	end if
417
418
419	! Now the other Gauss points, if needed.
420
421	DO NG=1,L_NGAUSS-1
422	IF(TAUGSURF(NW,NG) .gt. TLIMIT) THEN
423
424	DO L=1,L_NLAYRAD-1
425	K = 2*L+1
426	DTAUI(L,nw,ng) = DTAUKI(K,NW,NG)+DTAUKI(K+1,NW,NG)! + 1.e-50
427
428	if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
429
430	WBARI(L,nw,ng) = btemp(L,nw) / DTAUI(L,NW,NG)
431
432	else
433	WBARI(L,nw,ng) = 0.0D0
434	DTAUI(L,NW,NG) = 1.0D-9
435	endif
436
437	cosbi(L,NW,NG) = cosbi(L,NW,L_NGAUSS)
438	END DO ! L vertical loop
439
440	! Last level
441	L = L_NLAYRAD
442	K = 2*L+1
443	DTAUI(L,nw,ng) = DTAUKI(K,NW,NG)! + 1.e-50
444
445	if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
446
447	WBARI(L,nw,ng) = btemp(L,nw) / DTAUI(L,NW,NG)
448
449	else
450	WBARI(L,nw,ng) = 0.0D0
451	DTAUI(L,NW,NG) = 1.0D-9
452	endif
453
454	cosbi(L,NW,NG) = cosbi(L,NW,L_NGAUSS)
455
456	END IF
457
458	END DO ! NG Gauss loop
459	END DO ! NW spectral loop
460
461	! Total extinction optical depths
462
463	DO NG=1,L_NGAUSS ! full gauss loop
464	DO NW=1,L_NSPECTI
465	TAUCUMI(1,NW,NG)=0.0D0
466	DO K=2,L_LEVELS
467	TAUCUMI(K,NW,NG)=TAUCUMI(K-1,NW,NG)+DTAUKI(K,NW,NG)
468	END DO
469	END DO ! end full gauss loop
470	END DO
471
472	! be aware when comparing with textbook results
473	! (e.g. Pierrehumbert p. 218) that
474	! taucumi does not take the <cos theta>=0.5 factor into
475	! account. It is the optical depth for a vertically
476	! ascending ray with angle theta = 0.
477
478	!open(127,file='taucum.out')
479	!do nw=1,L_NSPECTI
480	! write(127,*) taucumi(L_LEVELS,nw,L_NGAUSS)
481	!enddo
482	!close(127)
483
484	! print*,'WBARI'
485	! print*,WBARI
486	! print*,'DTAUI'
487	! print*,DTAUI
488	! call abort
489
490	end subroutine optci
491
492	END MODULE optci_mod
493

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

Download in other formats:

Original Format