Context Navigation

optci.F90 @ 1145

Last change on this file since 1145 was 1016, checked in by jleconte, 11 years ago

07/08/2013 == JL

Water cycle in double precision (largescale+moistadj)
Improved wate rayleigh.
First step for rayleigh with variable species. Now, just need to change optcv.
changed some interpolation indices in callcorrk to limit dependency of OLR on the number of layers

File size: 12.9 KB

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

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

Download in other formats:

Original Format