Context Navigation

optci.F90 @ 3232

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

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

Download in other formats:

Original Format