Context Navigation

optci.F90 @ 1716

Last change on this file since 1716 was 1715, checked in by jvatant, 9 years ago

Cleaning of the radiative code + coherence between optci and optcv :

+ Get rid of the actually dummy and confusing extra level L_LEVELS+1 for aerosols
+ Harmonisation between optci and optcv
+ Get rid of tauref (called nowhere) and csqtly of ini_radcommon subroutine

JVO

File size: 14.2 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,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,L_NSPECTI,NAERKIND)
45	real*8 QSIAER(L_LEVELS,L_NSPECTI,NAERKIND)
46	real*8 GIAER(L_LEVELS,L_NSPECTI,NAERKIND)
47	real*8 TAUAERO(L_LEVELS,NAERKIND)
48	real*8 TAUAEROLK(L_LEVELS,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	! Spectral dependance of aerosol absorption
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	DAERO=SUM(TAEROS(K,NW,1:naerkind)) ! aerosol absorption
151
152	DCONT = 0.0d0 ! continuum absorption
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	do ng=1,L_NGAUSS-1
229
230	! Now compute TAUGAS
231
232	! Interpolate between water mixing ratios
233	! WRATIO = 0.0 if the requested water amount is equal to, or outside the
234	! the water data range
235
236	if(L_REFVAR.eq.1)then ! added by RW for special no variable case
237	KCOEF(1) = GASI(MT(K),MP(K),1,NW,NG)
238	KCOEF(2) = GASI(MT(K),MP(K)+1,1,NW,NG)
239	KCOEF(3) = GASI(MT(K)+1,MP(K)+1,1,NW,NG)
240	KCOEF(4) = GASI(MT(K)+1,MP(K),1,NW,NG)
241	else
242
243	KCOEF(1) = GASI(MT(K),MP(K),NVAR(K),NW,NG) + WRATIO(K)* &
244	(GASI(MT(K),MP(K),NVAR(K)+1,NW,NG) - &
245	GASI(MT(K),MP(K),NVAR(K),NW,NG))
246
247	KCOEF(2) = GASI(MT(K),MP(K)+1,NVAR(K),NW,NG) + WRATIO(K)* &
248	(GASI(MT(K),MP(K)+1,NVAR(K)+1,NW,NG) - &
249	GASI(MT(K),MP(K)+1,NVAR(K),NW,NG))
250
251	KCOEF(3) = GASI(MT(K)+1,MP(K)+1,NVAR(K),NW,NG) + WRATIO(K)* &
252	(GASI(MT(K)+1,MP(K)+1,NVAR(K)+1,NW,NG) - &
253	GASI(MT(K)+1,MP(K)+1,NVAR(K),NW,NG))
254
255	KCOEF(4) = GASI(MT(K)+1,MP(K),NVAR(K),NW,NG) + WRATIO(K)* &
256	(GASI(MT(K)+1,MP(K),NVAR(K)+1,NW,NG) - &
257	GASI(MT(K)+1,MP(K),NVAR(K),NW,NG))
258
259	endif
260
261	! Interpolate the gaseous k-coefficients to the requested T,P values
262
263	ANS = LKCOEF(K,1)KCOEF(1) + LKCOEF(K,2)KCOEF(2) + &
264	LKCOEF(K,3)KCOEF(3) + LKCOEF(K,4)KCOEF(4)
265
266	TAUGAS = U(k)*ANS
267
268	TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS + DCONT
269	DTAUKI(K,nw,ng) = TAUGAS &
270	+ DCONT & ! For parameterized continuum absorption
271	+ DAERO ! For aerosol absorption
272
273	end do
274
275	! Now fill in the "clear" part of the spectrum (NG = L_NGAUSS),
276	! which holds continuum opacity only
277
278	NG = L_NGAUSS
279	DTAUKI(K,nw,ng) = 0.d0 &
280	+ DCONT & ! For parameterized continuum absorption
281	+ DAERO ! For aerosol absorption
282
283	end do
284	end do
285
286	!=======================================================================
287	! Now the full treatment for the layers, where besides the opacity
288	! we need to calculate the scattering albedo and asymmetry factors
289
290	do iaer=1,naerkind
291	DO NW=1,L_NSPECTI
292	DO K=2,L_LEVELS
293	TAUAEROLK(K,NW,IAER) = TAUAERO(K,IAER)*QSIAER(K,NW,IAER) ! effect of scattering albedo
294	ENDDO
295	ENDDO
296	end do
297
298	DO NW=1,L_NSPECTI
299	DO L=1,L_NLAYRAD-1
300	K = 2*L+1
301	btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind)) + SUM(TAUAEROLK(K+1,NW,1:naerkind))
302	END DO ! L vertical loop
303
304	! Last level
305	L = L_NLAYRAD
306	K = 2*L+1
307	btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind))
308
309	END DO ! NW spectral loop
310
311
312	DO NW=1,L_NSPECTI
313	NG = L_NGAUSS
314	DO L=1,L_NLAYRAD-1
315
316	K = 2*L+1
317	DTAUI(L,nw,ng) = DTAUKI(K,NW,NG) + DTAUKI(K+1,NW,NG)! + 1.e-50
318
319	atemp = 0.
320	if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
321	do iaer=1,naerkind
322	atemp = atemp + &
323	GIAER(K,NW,IAER) * TAUAEROLK(K,NW,IAER) + &
324	GIAER(K+1,NW,IAER) * TAUAEROLK(K+1,NW,IAER)
325	end do
326	WBARI(L,nw,ng) = btemp(L,nw) / DTAUI(L,NW,NG)
327	else
328	WBARI(L,nw,ng) = 0.0D0
329	DTAUI(L,NW,NG) = 1.0D-9
330	endif
331
332	if(btemp(L,nw) .GT. 0.0d0) then
333	cosbi(L,NW,NG) = atemp/btemp(L,nw)
334	else
335	cosbi(L,NW,NG) = 0.0D0
336	end if
337
338	END DO ! L vertical loop
339
340	! Last level
341
342	L = L_NLAYRAD
343	K = 2*L+1
344	DTAUI(L,nw,ng) = DTAUKI(K,NW,NG) ! + 1.e-50
345
346	atemp = 0.
347	if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
348	do iaer=1,naerkind
349	atemp = atemp + GIAER(K,NW,IAER) * TAUAEROLK(K,NW,IAER)
350	end do
351	WBARI(L,nw,ng) = btemp(L,nw) / DTAUI(L,NW,NG)
352	else
353	WBARI(L,nw,ng) = 0.0D0
354	DTAUI(L,NW,NG) = 1.0D-9
355	endif
356
357	if(btemp(L,nw) .GT. 0.0d0) then
358	cosbi(L,NW,NG) = atemp/btemp(L,nw)
359	else
360	cosbi(L,NW,NG) = 0.0D0
361	end if
362
363
364	! Now the other Gauss points, if needed.
365
366	DO NG=1,L_NGAUSS-1
367	IF(TAUGSURF(NW,NG) .gt. TLIMIT) THEN
368
369	DO L=1,L_NLAYRAD-1
370	K = 2*L+1
371	DTAUI(L,nw,ng) = DTAUKI(K,NW,NG)+DTAUKI(K+1,NW,NG)! + 1.e-50
372
373	if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
374
375	WBARI(L,nw,ng) = btemp(L,nw) / DTAUI(L,NW,NG)
376
377	else
378	WBARI(L,nw,ng) = 0.0D0
379	DTAUI(L,NW,NG) = 1.0D-9
380	endif
381
382	cosbi(L,NW,NG) = cosbi(L,NW,L_NGAUSS)
383	END DO ! L vertical loop
384
385	! Last level
386	L = L_NLAYRAD
387	K = 2*L+1
388	DTAUI(L,nw,ng) = DTAUKI(K,NW,NG)! + 1.e-50
389
390	if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
391
392	WBARI(L,nw,ng) = btemp(L,nw) / DTAUI(L,NW,NG)
393
394	else
395	WBARI(L,nw,ng) = 0.0D0
396	DTAUI(L,NW,NG) = 1.0D-9
397	endif
398
399	cosbi(L,NW,NG) = cosbi(L,NW,L_NGAUSS)
400
401	END IF
402
403	END DO ! NG Gauss loop
404	END DO ! NW spectral loop
405
406	! Total extinction optical depths
407
408	DO NG=1,L_NGAUSS ! full gauss loop
409	DO NW=1,L_NSPECTI
410	TAUCUMI(1,NW,NG)=0.0D0
411	DO K=2,L_LEVELS
412	TAUCUMI(K,NW,NG)=TAUCUMI(K-1,NW,NG)+DTAUKI(K,NW,NG)
413	END DO
414	END DO ! end full gauss loop
415	END DO
416
417	! be aware when comparing with textbook results
418	! (e.g. Pierrehumbert p. 218) that
419	! taucumi does not take the <cos theta>=0.5 factor into
420	! account. It is the optical depth for a vertically
421	! ascending ray with angle theta = 0.
422
423	!open(127,file='taucum.out')
424	!do nw=1,L_NSPECTI
425	! write(127,*) taucumi(L_LEVELS,nw,L_NGAUSS)
426	!enddo
427	!close(127)
428
429	! print*,'WBARI'
430	! print*,WBARI
431	! print*,'DTAUI'
432	! print*,DTAUI
433	! call abort
434
435
436	return
437
438
439	end subroutine optci
440
441
442

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

Download in other formats:

Original Format