Context Navigation

optci.F90 @ 3700

Last change on this file since 3700 was 3373, checked in by afalco, 18 months ago
Pluto.old PCM: Corrected some runtime issues with gfortran. AF
File size: 9.4 KB

Line
1	subroutine optci(PLEV,TLEV,DTAUI,TAUCUMI, &
2	QXIAER,QSIAER,GIAER,COSBI,WBARI,TAUAERO, &
3	TMID,PMID,TAUGSURF,QVAR)
4
5
6	use radinc_h
7	use radcommon_h, only: gasi, tlimit, wrefVAR, Cmk,tgasref,pfgasref,wnoi
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
31	#include "comcstfi.h"
32	#include "callkeys.h"
33
34
35	real*8 DTAUI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
36	real*8 DTAUKI(L_LEVELS+1,L_NSPECTI,L_NGAUSS)
37	real*8 TAUI(L_NLEVRAD,L_NSPECTI,L_NGAUSS)
38	real*8 TAUCUMI(L_LEVELS,L_NSPECTI,L_NGAUSS)
39	real*8 PLEV(L_LEVELS)
40	real*8 TLEV(L_LEVELS)
41	real*8 TMID(L_LEVELS), PMID(L_LEVELS)
42	real*8 COSBI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
43	real*8 WBARI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
44
45	! For aerosols
46	real*8 QXIAER(L_LEVELS+1,L_NSPECTI,NAERKIND)
47	real*8 QSIAER(L_LEVELS+1,L_NSPECTI,NAERKIND)
48	real*8 GIAER(L_LEVELS+1,L_NSPECTI,NAERKIND)
49	real*8 TAUAERO(L_LEVELS+1,NAERKIND)
50	real*8 TAUAEROLK(L_LEVELS+1,L_NSPECTI,NAERKIND)
51	real*8 TAEROS(L_LEVELS,L_NSPECTI,NAERKIND)
52
53	integer L, NW, NG, K, LK, IAER
54	integer MT(L_LEVELS), MP(L_LEVELS), NP(L_LEVELS)
55	real*8 ANS, TAUGAS
56	real*8 DPR(L_LEVELS), U(L_LEVELS)
57	real*8 LCOEF(4), LKCOEF(L_LEVELS,4)
58
59	real*8 taugsurf(L_NSPECTI,L_NGAUSS-1)
60	real*8 dco2
61
62	! mixing ratio variables
63	real*8 QVAR(L_LEVELS), WRATIO(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 rho, dz
73
74	!=======================================================================
75	! Determine the total gas opacity throughout the column, for each
76	! spectral interval, NW, and each Gauss point, NG.
77
78	! write(,)'L_LEVELS',L_LEVELS
79	! write(,)'L_NSPECTI',L_NSPECTI
80	DTAUI(:,:,:)=0.
81	DTAUKI(:,:,:)=0.
82
83	DO NG=1,L_NGAUSS-1
84	do NW=1,L_NSPECTI
85	TAUGSURF(NW,NG) = 0.0D0
86	end do
87	end do
88	do K=2,L_LEVELS
89	DPR(k) = PLEV(K)-PLEV(K-1)
90
91	! rho = PLEV(K)/(R*TMID(K))
92	rho = PMID(K)/(R*TMID(K))
93	dz = -DPR(k)/(g*rho)
94	!print*,'rho=',rho
95	!print*,'dz=',dz
96
97	U(k) = Cmk*DPR(k) ! only Cmk line in optci.F
98	! soon to be replaced by m^-1 !!!
99
100	call tpindex(PMID(K),TMID(K),QVAR(K),pfgasref,tgasref,WREFVAR, &
101	LCOEF,MT(K),MP(K),NVAR(K),WRATIO(K))
102
103	do LK=1,4
104	LKCOEF(K,LK) = LCOEF(LK)
105	end do
106
107	DO NW=1,L_NSPECTI
108	do iaer=1,naerkind
109	TAEROS(K,NW,IAER) = TAUAERO(K,IAER) * QXIAER(K,NW,IAER)
110	! write(22,*) 'TB17 Taero IR:',K,NW,IAER,TAEROS(K,NW,IAER)
111	end do
112	END DO
113	end do ! levels
114
115
116	do K=2,L_LEVELS
117	do nw=1,L_NSPECTI
118
119	DCO2 = 0.0 ! continuum absorption (no longer used)
120
121	do ng=1,L_NGAUSS-1
122
123	! Now compute TAUGAS
124	! Interpolate between mixing ratios
125	! WRATIO = 0.0 if the requested amount is equal to, or outside the
126	! the data range
127
128
129	if (L_REFVAR.eq.1)then ! added by RW for special no variable case
130	KCOEF(1) = GASI(MT(K),MP(K),1,NW,NG)
131	KCOEF(2) = GASI(MT(K),MP(K)+1,1,NW,NG)
132	KCOEF(3) = GASI(MT(K)+1,MP(K)+1,1,NW,NG)
133	KCOEF(4) = GASI(MT(K)+1,MP(K),1,NW,NG)
134	else
135
136	KCOEF(1) = GASI(MT(K),MP(K),NVAR(K),NW,NG) + WRATIO(K)* &
137	(GASI(MT(K),MP(K),NVAR(K)+1,NW,NG) - &
138	GASI(MT(K),MP(K),NVAR(K),NW,NG))
139
140	KCOEF(2) = GASI(MT(K),MP(K)+1,NVAR(K),NW,NG)+ WRATIO(K)* &
141	(GASI(MT(K),MP(K)+1,NVAR(K)+1,NW,NG) - &
142	GASI(MT(K),MP(K)+1,NVAR(K),NW,NG))
143
144	KCOEF(3)=GASI(MT(K)+1,MP(K)+1,NVAR(K),NW,NG)+WRATIO(K)* &
145	(GASI(MT(K)+1,MP(K)+1,NVAR(K)+1,NW,NG) - &
146	GASI(MT(K)+1,MP(K)+1,NVAR(K),NW,NG))
147
148	KCOEF(4) =GASI(MT(K)+1,MP(K),NVAR(K),NW,NG) + WRATIO(K)* &
149	(GASI(MT(K)+1,MP(K),NVAR(K)+1,NW,NG) - &
150	GASI(MT(K)+1,MP(K),NVAR(K),NW,NG))
151	endif
152
153	! Interpolate the gaseous k-coefficients to the requested T,P values
154
155	ANS = LKCOEF(K,1)KCOEF(1) + LKCOEF(K,2)KCOEF(2) + &
156	LKCOEF(K,3)KCOEF(3) + LKCOEF(K,4)KCOEF(4)
157
158	TAUGAS = U(k)*ANS
159
160	TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS
161
162	DTAUKI(K,nw,ng) = TAUGAS
163	do iaer=1,naerkind
164	DTAUKI(K,nw,ng) = DTAUKI(K,nw,ng) + TAEROS(K,NW,IAER) &
165	+ DCO2 ! For Kasting CIA
166	end do
167
168	end do
169
170	! Now fill in the "clear" part of the spectrum (NG = L_NGAUSS),
171	! which holds continuum opacity only
172
173	NG = L_NGAUSS
174	DTAUKI(K,nw,ng) = 0.0
175	do iaer=1,naerkind
176	DTAUKI(K,nw,ng) = DTAUKI(K,nw,ng) + TAEROS(K,NW,IAER) &
177	+ DCO2 ! For parameterized continuum absorption
178	end do ! a bug was found here!!
179
180	end do
181	end do
182
183
184	!=======================================================================
185	! Now the full treatment for the layers, where besides the opacity
186	! we need to calculate the scattering albedo and asymmetry factors
187
188	DO NW=1,L_NSPECTI
189	DO K=2,L_LEVELS
190	do iaer=1,naerkind
191	TAUAEROLK(K,NW,IAER) = TAUAERO(K,IAER)*QSIAER(K,NW,IAER)
192	end do
193	ENDDO
194	ENDDO
195
196	DO NW=1,L_NSPECTI
197	NG = L_NGAUSS
198	DO L=1,L_NLAYRAD-1
199
200	K = 2*L+1
201	DTAUI(L,nw,ng) = DTAUKI(K,NW,NG) + DTAUKI(K+1,NW,NG)! + 1.e-50
202
203	atemp = 0.
204	btemp(L,NW) = 0.
205	do iaer=1,naerkind
206	atemp = atemp + &
207	GIAER(K,NW,IAER) * TAUAEROLK(K,NW,IAER) + &
208	GIAER(K+1,NW,IAER) * TAUAEROLK(K+1,NW,IAER)
209	btemp(L,NW) = btemp(L,NW) + TAUAEROLK(K,NW,IAER) + TAUAEROLK(K+1,NW,IAER)
210	! * + 1.e-10
211	end do
212
213	if(DTAUI(L,NW,NG) .GT. 1.0E-9) then
214	WBARI(L,nw,ng) = btemp(L,NW) / DTAUI(L,NW,NG)
215	else
216	WBARI(L,nw,ng) = 0.0D0
217	DTAUI(L,NW,NG) = 1.0E-9
218	endif
219
220	if(btemp(L,NW) .GT. 0.0) then
221	cosbi(L,NW,NG) = atemp/btemp(L,NW)
222	else
223	cosbi(L,NW,NG) = 0.0D0
224	end if
225
226	END DO ! L vertical loop
227	! Last level
228
229	L = L_NLAYRAD
230	K = 2*L+1
231	DTAUI(L,nw,ng) = DTAUKI(K,NW,NG) ! + 1.e-50
232	btemp(L,NW) = 0
233	do iaer=1,naerkind
234	btemp(L,NW) = btemp(L,NW) + TAUAEROLK(K,NW,IAER)
235	enddo
236
237	atemp = 0.
238	if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
239	do iaer=1,naerkind
240	atemp = atemp + GIAER(K,NW,IAER) * TAUAEROLK(K,NW,IAER)
241	end do
242	WBARI(L,nw,ng) = btemp(L,NW) / DTAUI(L,NW,NG)
243	else
244	WBARI(L,nw,ng) = 0.0D0
245	DTAUI(L,NW,NG) = 1.0D-9
246	endif
247
248	if(btemp(L,NW) .GT. 0.0d0) then
249	cosbi(L,NW,NG) = atemp/btemp(L,NW)
250	else
251	cosbi(L,NW,NG) = 0.0D0
252	end if
253
254	! Now the other Gauss points, if needed.
255
256	DO NG=1,L_NGAUSS-1
257	IF(TAUGSURF(NW,NG) .gt. TLIMIT) THEN
258
259	DO L=1,L_NLAYRAD
260	K = 2*L+1
261	DTAUI(L,nw,ng) = DTAUKI(K,NW,NG)+DTAUKI(K+1,NW,NG)! + 1.e-50
262
263	if(DTAUI(L,NW,NG) .GT. 1.0E-9) then
264	WBARI(L,nw,ng) = btemp(L,NW) / DTAUI(L,NW,NG)
265	else
266	WBARI(L,nw,ng) = 0.0D0
267	DTAUI(L,NW,NG) = 1.0E-9
268	endif
269
270	cosbi(L,NW,NG) = cosbi(L,NW,L_NGAUSS)
271	END DO ! L vertical loop
272	END IF
273
274	END DO ! NG Gauss loop
275	END DO ! NW spectral loop
276
277	! Total extinction optical depths
278
279	DO NW=1,L_NSPECTI
280	DO NG=1,L_NGAUSS ! full gauss loop
281	TAUI(1,NW,NG)=0.0D0
282	DO L=1,L_NLAYRAD
283	TAUI(L+1,NW,NG)=TAUI(L,NW,NG)+DTAUI(L,NW,NG)
284	END DO
285
286	TAUCUMI(1,NW,NG)=0.0D0
287	DO K=2,L_LEVELS
288	TAUCUMI(K,NW,NG)=TAUCUMI(K-1,NW,NG)+DTAUKI(K,NW,NG)
289	END DO
290	END DO ! end full gauss loop
291	END DO
292
293
294	return
295	end subroutine optci
296

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

Download in other formats:

Original Format