Context Navigation

optci.F90 @ 1477

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

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

Download in other formats:

Original Format