Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

optci_pluto_mod.F90 @ 3504

Last change on this file since 3504 was 3329, checked in by afalco, 7 months ago
Pluto PCM: Include cooling, hazes in radiative module AF
File size: 10.1 KB

Line
1	MODULE optci_pluto_mod
2
3	IMPLICIT NONE
4
5	CONTAINS
6
7	subroutine optci_pluto(PLEV,TLEV,DTAUI,TAUCUMI, &
8	QXIAER,QSIAER,GIAER,COSBI,WBARI,TAUAERO, &
9	TMID,PMID,TAUGSURF,QVAR)
10
11	use radinc_h, only: L_LEVELS, L_NLAYRAD, L_NSPECTI, L_NGAUSS, &
12	L_NLEVRAD, L_REFVAR, naerkind
13	use radcommon_h, only: gasi,tlimit,wrefVAR,Cmk,tgasref,pfgasref,wnoi,scalep,indi,glat_ig
14	use gases_h, only: gfrac, ngasmx, igas_N2, igas_He, igas_H2O, igas_H2, &
15	igas_CH4, igas_N2
16	use comcstfi_mod, only: g, r, mugaz
17	use callkeys_mod, only: kastprof,continuum,graybody
18	use recombin_corrk_mod, only: corrk_recombin, gasi_recomb
19	use tpindex_mod, only: tpindex
20	implicit none
21
22	!==================================================================
23	!
24	! Purpose
25	! -------
26	! Calculates longwave optical constants at each level. For each
27	! layer and spectral interval in the IR it calculates WBAR, DTAU
28	! and COSBAR. For each level it calculates TAU.
29	!
30	! TAUI(L,LW) is the cumulative optical depth at level L (or alternatively
31	! at the bottom of layer L), LW is the spectral wavelength interval.
32	!
33	! TLEV(L) - Temperature at the layer boundary (i.e., level)
34	! PLEV(L) - Pressure at the layer boundary (i.e., level)
35	!
36	! Authors
37	! -------
38	! Adapted from the NASA Ames code by R. Wordsworth (2009)
39	!
40	!==================================================================
41
42
43
44	real*8 DTAUI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
45	real*8 DTAUKI(L_LEVELS+1,L_NSPECTI,L_NGAUSS)
46	real*8 TAUI(L_NLEVRAD,L_NSPECTI,L_NGAUSS)
47	real*8 TAUCUMI(L_LEVELS,L_NSPECTI,L_NGAUSS)
48	real*8 PLEV(L_LEVELS)
49	real*8 TLEV(L_LEVELS)
50	real*8 TMID(L_LEVELS), PMID(L_LEVELS)
51	real*8 COSBI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
52	real*8 WBARI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
53
54	! For aerosols
55	real*8 QXIAER(L_LEVELS+1,L_NSPECTI,NAERKIND)
56	real*8 QSIAER(L_LEVELS+1,L_NSPECTI,NAERKIND)
57	real*8 GIAER(L_LEVELS+1,L_NSPECTI,NAERKIND)
58	real*8 TAUAERO(L_LEVELS+1,NAERKIND)
59	real*8 TAUAEROLK(L_LEVELS+1,L_NSPECTI,NAERKIND)
60	real*8 TAEROS(L_LEVELS,L_NSPECTI,NAERKIND)
61
62	integer L, NW, NG, K, LK, IAER
63	integer MT(L_LEVELS), MP(L_LEVELS), NP(L_LEVELS)
64	real*8 ANS, TAUGAS
65	real*8 DPR(L_LEVELS), U(L_LEVELS)
66	real*8 LCOEF(4), LKCOEF(L_LEVELS,4)
67
68	real*8 taugsurf(L_NSPECTI,L_NGAUSS-1)
69	real*8 dco2
70
71	! mixing ratio variables
72	real*8 QVAR(L_LEVELS), WRATIO(L_LEVELS)
73	real*8 KCOEF(4)
74	integer NVAR(L_LEVELS)
75
76	! temporary variables for multiple aerosol calculation
77	real*8 atemp, btemp
78
79	! variables for k in units m^-1
80	real*8 rho, dz
81
82	!=======================================================================
83	! Determine the total gas opacity throughout the column, for each
84	! spectral interval, NW, and each Gauss point, NG.
85
86	! write(,)'L_LEVELS',L_LEVELS
87	! write(,)'L_NSPECTI',L_NSPECTI
88	DTAUI(:,:,:)=0.
89	DTAUKI(:,:,:)=0.
90
91	DO NG=1,L_NGAUSS-1
92	do NW=1,L_NSPECTI
93	TAUGSURF(NW,NG) = 0.0D0
94	end do
95	end do
96	do K=2,L_LEVELS
97	DPR(k) = PLEV(K)-PLEV(K-1)
98
99	! rho = PLEV(K)/(R*TMID(K))
100	rho = PMID(K)/(R*TMID(K))
101	dz = -DPR(k)/(g*rho)
102	!print*,'rho=',rho
103	!print*,'dz=',dz
104
105	U(k) = Cmk*DPR(k) ! only Cmk line in optci_pluto.F
106	! soon to be replaced by m^-1 !!!
107
108	call tpindex(PMID(K),TMID(K),QVAR(K),pfgasref,tgasref,WREFVAR, &
109	LCOEF,MT(K),MP(K),NVAR(K),WRATIO(K))
110
111	do LK=1,4
112	LKCOEF(K,LK) = LCOEF(LK)
113	end do
114
115	DO NW=1,L_NSPECTI
116	do iaer=1,naerkind
117	TAEROS(K,NW,IAER) = TAUAERO(K,IAER) * QXIAER(K,NW,IAER)
118	! write(22,*) 'TB17 Taero IR:',K,NW,IAER,TAEROS(K,NW,IAER)
119	end do
120	END DO
121	end do ! levels
122
123
124	do K=2,L_LEVELS
125	do nw=1,L_NSPECTI
126
127	DCO2 = 0.0 ! continuum absorption (no longer used)
128
129	do ng=1,L_NGAUSS-1
130
131	! Now compute TAUGAS
132	! Interpolate between mixing ratios
133	! WRATIO = 0.0 if the requested amount is equal to, or outside the
134	! the data range
135
136
137	if (L_REFVAR.eq.1)then ! added by RW for special no variable case
138	KCOEF(1) = GASI(MT(K),MP(K),1,NW,NG)
139	KCOEF(2) = GASI(MT(K),MP(K)+1,1,NW,NG)
140	KCOEF(3) = GASI(MT(K)+1,MP(K)+1,1,NW,NG)
141	KCOEF(4) = GASI(MT(K)+1,MP(K),1,NW,NG)
142	else
143
144	KCOEF(1) = GASI(MT(K),MP(K),NVAR(K),NW,NG) + WRATIO(K)* &
145	(GASI(MT(K),MP(K),NVAR(K)+1,NW,NG) - &
146	GASI(MT(K),MP(K),NVAR(K),NW,NG))
147
148	KCOEF(2) = GASI(MT(K),MP(K)+1,NVAR(K),NW,NG)+ WRATIO(K)* &
149	(GASI(MT(K),MP(K)+1,NVAR(K)+1,NW,NG) - &
150	GASI(MT(K),MP(K)+1,NVAR(K),NW,NG))
151
152	KCOEF(3)=GASI(MT(K)+1,MP(K)+1,NVAR(K),NW,NG)+WRATIO(K)* &
153	(GASI(MT(K)+1,MP(K)+1,NVAR(K)+1,NW,NG) - &
154	GASI(MT(K)+1,MP(K)+1,NVAR(K),NW,NG))
155
156	KCOEF(4) =GASI(MT(K)+1,MP(K),NVAR(K),NW,NG) + WRATIO(K)* &
157	(GASI(MT(K)+1,MP(K),NVAR(K)+1,NW,NG) - &
158	GASI(MT(K)+1,MP(K),NVAR(K),NW,NG))
159	endif
160
161	! Interpolate the gaseous k-coefficients to the requested T,P values
162
163	ANS = LKCOEF(K,1)KCOEF(1) + LKCOEF(K,2)KCOEF(2) + &
164	LKCOEF(K,3)KCOEF(3) + LKCOEF(K,4)KCOEF(4)
165
166	TAUGAS = U(k)*ANS
167
168	TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS
169
170	DTAUKI(K,nw,ng) = TAUGAS
171	do iaer=1,naerkind
172	DTAUKI(K,nw,ng) = DTAUKI(K,nw,ng) + TAEROS(K,NW,IAER) &
173	+ DCO2 ! For Kasting CIA
174	end do
175
176	end do
177
178	! Now fill in the "clear" part of the spectrum (NG = L_NGAUSS),
179	! which holds continuum opacity only
180
181	NG = L_NGAUSS
182	DTAUKI(K,nw,ng) = 0.0
183	do iaer=1,naerkind
184	DTAUKI(K,nw,ng) = DTAUKI(K,nw,ng) + TAEROS(K,NW,IAER) &
185	+ DCO2 ! For parameterized continuum absorption
186	end do ! a bug was found here!!
187
188	end do
189	end do
190
191
192	!=======================================================================
193	! Now the full treatment for the layers, where besides the opacity
194	! we need to calculate the scattering albedo and asymmetry factors
195
196	DO NW=1,L_NSPECTI
197	DO K=2,L_LEVELS
198	do iaer=1,naerkind
199	TAUAEROLK(K,NW,IAER) = TAUAERO(K,IAER)*QSIAER(K,NW,IAER)
200	end do
201	ENDDO
202	ENDDO
203
204	DO NW=1,L_NSPECTI
205	NG = L_NGAUSS
206	DO L=1,L_NLAYRAD-1
207
208	K = 2*L+1
209	DTAUI(L,nw,ng) = DTAUKI(K,NW,NG) + DTAUKI(K+1,NW,NG)! + 1.e-50
210
211	atemp = 0.
212	btemp = 0.
213	if(DTAUI(L,NW,NG) .GT. 1.0E-9) then
214	do iaer=1,naerkind
215	atemp = atemp + &
216	GIAER(K,NW,IAER) * TAUAEROLK(K,NW,IAER) + &
217	GIAER(K+1,NW,IAER) * TAUAEROLK(K+1,NW,IAER)
218	btemp = btemp + TAUAEROLK(K,NW,IAER) + TAUAEROLK(K+1,NW,IAER)
219	! * + 1.e-10
220	end do
221	WBARI(L,nw,ng) = btemp / DTAUI(L,NW,NG)
222	else
223	WBARI(L,nw,ng) = 0.0D0
224	DTAUI(L,NW,NG) = 1.0E-9
225	endif
226
227	if(btemp .GT. 0.0) then
228	cosbi(L,NW,NG) = atemp/btemp
229	else
230	cosbi(L,NW,NG) = 0.0D0
231	end if
232
233	END DO ! L vertical loop
234
235	! Last level
236
237	L = L_NLAYRAD
238	K = 2*L+1
239	DTAUI(L,nw,ng) = DTAUKI(K,NW,NG) ! + 1.e-50
240
241	atemp = 0.
242	if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
243	do iaer=1,naerkind
244	atemp = atemp + GIAER(K,NW,IAER) * TAUAEROLK(K,NW,IAER)
245	end do
246	WBARI(L,nw,ng) = btemp / DTAUI(L,NW,NG)
247	else
248	WBARI(L,nw,ng) = 0.0D0
249	DTAUI(L,NW,NG) = 1.0D-9
250	endif
251
252	if(btemp .GT. 0.0d0) then
253	cosbi(L,NW,NG) = atemp/btemp
254	else
255	cosbi(L,NW,NG) = 0.0D0
256	end if
257
258	! Now the other Gauss points, if needed.
259
260	DO NG=1,L_NGAUSS-1
261	IF(TAUGSURF(NW,NG) .gt. TLIMIT) THEN
262
263	DO L=1,L_NLAYRAD-1
264	K = 2*L+1
265	DTAUI(L,nw,ng) = DTAUKI(K,NW,NG)+DTAUKI(K+1,NW,NG)! + 1.e-50
266
267	btemp = 0.
268	if(DTAUI(L,NW,NG) .GT. 1.0E-9) then
269
270	do iaer=1,naerkind
271	btemp = btemp + TAUAEROLK(K,NW,IAER) + TAUAEROLK(K+1,NW,IAER)
272	end do
273	WBARI(L,nw,ng) = btemp / DTAUI(L,NW,NG)
274
275	else
276	WBARI(L,nw,ng) = 0.0D0
277	DTAUI(L,NW,NG) = 1.0E-9
278	endif
279
280	cosbi(L,NW,NG) = cosbi(L,NW,L_NGAUSS)
281	END DO ! L vertical loop
282
283	! Last level
284	L = L_NLAYRAD
285	K = 2*L+1
286	btemp = SUM(TAUAEROLK(K,NW,1:naerkind))
287	DTAUI(L,nw,ng) = DTAUKI(K,NW,NG)! + 1.e-50
288
289	if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
290
291	WBARI(L,nw,ng) = btemp / DTAUI(L,NW,NG)
292
293	else
294	WBARI(L,nw,ng) = 0.0D0
295	DTAUI(L,NW,NG) = 1.0D-9
296	endif
297
298	cosbi(L,NW,NG) = cosbi(L,NW,L_NGAUSS)
299
300	END IF
301
302	END DO ! NG Gauss loop
303	END DO ! NW spectral loop
304
305	! Total extinction optical depths
306
307	DO NW=1,L_NSPECTI
308	DO NG=1,L_NGAUSS ! full gauss loop
309	TAUI(1,NW,NG)=0.0D0
310	DO L=1,L_NLAYRAD
311	TAUI(L+1,NW,NG)=TAUI(L,NW,NG)+DTAUI(L,NW,NG)
312	END DO
313
314	TAUCUMI(1,NW,NG)=0.0D0
315	DO K=2,L_LEVELS
316	TAUCUMI(K,NW,NG)=TAUCUMI(K-1,NW,NG)+DTAUKI(K,NW,NG)
317	END DO
318	END DO ! end full gauss loop
319	END DO
320
321
322	end subroutine optci_pluto
323
324	END MODULE optci_pluto_mod
325

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

Download in other formats:

Original Format