Context Navigation

optcv.F90 @ 1145

Last change on this file since 1145 was 1016, checked in by jleconte, 11 years ago

07/08/2013 == JL

Water cycle in double precision (largescale+moistadj)
Improved wate rayleigh.
First step for rayleigh with variable species. Now, just need to change optcv.
changed some interpolation indices in callcorrk to limit dependency of OLR on the number of layers

Property svn:executable set to *

File size: 11.7 KB

Line
1	SUBROUTINE OPTCV(DTAUV,TAUV,TAUCUMV,PLEV, &
2	QXVAER,QSVAER,GVAER,WBARV,COSBV, &
3	TAURAY,TAUAERO,TMID,PMID,TAUGSURF,QVAR,MUVAR)
4
5	use radinc_h
6	use radcommon_h, only: gasv, tlimit, wrefVAR, Cmk, tgasref, pfgasref,wnov,scalep,indv
7	use gases_h
8
9	implicit none
10
11	!==================================================================
12	!
13	! Purpose
14	! -------
15	! Calculates shortwave optical constants at each level.
16	!
17	! Authors
18	! -------
19	! Adapted from the NASA Ames code by R. Wordsworth (2009)
20	!
21	!==================================================================
22	!
23	! THIS SUBROUTINE SETS THE OPTICAL CONSTANTS IN THE VISUAL
24	! IT CALCUALTES FOR EACH LAYER, FOR EACH SPECRAL INTERVAL IN THE VISUAL
25	! LAYER: WBAR, DTAU, COSBAR
26	! LEVEL: TAU
27	!
28	! TAUV(L,NW,NG) is the cumulative optical depth at the top of radiation code
29	! layer L. NW is spectral wavelength interval, ng the Gauss point index.
30	!
31	! TLEV(L) - Temperature at the layer boundary
32	! PLEV(L) - Pressure at the layer boundary (i.e. level)
33	! GASV(NT,NPS,NW,NG) - Visible k-coefficients
34	!
35	!-------------------------------------------------------------------
36
37	#include "comcstfi.h"
38	#include "callkeys.h"
39
40
41	real*8 DTAUV(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
42	real*8 DTAUKV(L_LEVELS+1,L_NSPECTV,L_NGAUSS)
43	real*8 TAUV(L_NLEVRAD,L_NSPECTV,L_NGAUSS)
44	real*8 TAUCUMV(L_LEVELS,L_NSPECTV,L_NGAUSS)
45	real*8 PLEV(L_LEVELS)
46	real*8 TMID(L_LEVELS), PMID(L_LEVELS)
47	real*8 COSBV(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
48	real*8 WBARV(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
49
50	! for aerosols
51	real*8 QXVAER(L_LEVELS+1,L_NSPECTV,NAERKIND)
52	real*8 QSVAER(L_LEVELS+1,L_NSPECTV,NAERKIND)
53	real*8 GVAER(L_LEVELS+1,L_NSPECTV,NAERKIND)
54	real*8 TAUAERO(L_LEVELS+1,NAERKIND)
55	real*8 TAUAEROLK(L_LEVELS+1,L_NSPECTV,NAERKIND)
56	real*8 TAEROS(L_LEVELS,L_NSPECTV,NAERKIND)
57
58	integer L, NW, NG, K, LK, IAER
59	integer MT(L_LEVELS), MP(L_LEVELS), NP(L_LEVELS)
60	real*8 ANS, TAUGAS
61	real*8 TAURAY(L_NSPECTV)
62	real*8 TRAY(L_LEVELS,L_NSPECTV)
63	real*8 TRAYAER
64	real*8 DPR(L_LEVELS), U(L_LEVELS)
65	real*8 LCOEF(4), LKCOEF(L_LEVELS,4)
66
67	real*8 taugsurf(L_NSPECTV,L_NGAUSS-1)
68	real*8 DCONT,DAERO
69	double precision wn_cont, p_cont, p_air, T_cont, dtemp, dtempc
70	double precision p_cross
71
72	! variable species mixing ratio variables
73	real*8 QVAR(L_LEVELS), WRATIO(L_LEVELS), MUVAR(L_LEVELS)
74	real*8 KCOEF(4)
75	integer NVAR(L_LEVELS)
76
77	! temporary variables for multiple aerosol calculation
78	real*8 atemp(L_NLAYRAD,L_NSPECTV)
79	real*8 btemp(L_NLAYRAD,L_NSPECTV)
80	real*8 ctemp(L_NLAYRAD,L_NSPECTV)
81
82	! variables for k in units m^-1
83	real*8 dz(L_LEVELS)
84
85	integer igas, jgas
86
87	integer interm
88
89	!! AS: to save time in computing continuum (see bilinearbig)
90	IF (.not.ALLOCATED(indv)) THEN
91	ALLOCATE(indv(L_NSPECTV,ngasmx,ngasmx))
92	indv = -9999 ! this initial value means "to be calculated"
93	ENDIF
94
95	!=======================================================================
96	! Determine the total gas opacity throughout the column, for each
97	! spectral interval, NW, and each Gauss point, NG.
98	! Calculate the continuum opacities, i.e., those that do not depend on
99	! NG, the Gauss index.
100
101	taugsurf(:,:) = 0.0
102	dpr(:) = 0.0
103	lkcoef(:,:) = 0.0
104
105	do K=2,L_LEVELS
106	DPR(k) = PLEV(K)-PLEV(K-1)
107
108	! if we have continuum opacities, we need dz
109	if(kastprof)then
110	dz(k) = dpr(k)(1000.0d08.3145d0/muvar(k))TMID(K)/(gPMID(K))
111	U(k) = CmkDPR(k)mugaz/muvar(k)
112	else
113	dz(k) = dpr(k)RTMID(K)/(gPMID(K))mugaz/muvar(k)
114	U(k) = CmkDPR(k)mugaz/muvar(k) ! only Cmk line in optci.F
115	!JL13 the mugaz/muvar factor takes into account water meanmolecular weight if water is present
116	endif
117
118	call tpindex(PMID(K),TMID(K),QVAR(K),pfgasref,tgasref,WREFVAR, &
119	LCOEF,MT(K),MP(K),NVAR(K),WRATIO(K))
120
121	do LK=1,4
122	LKCOEF(K,LK) = LCOEF(LK)
123	end do
124	end do ! levels
125
126
127	do iaer=1,naerkind
128	do NW=1,L_NSPECTV
129	do K=2,L_LEVELS
130	TAEROS(K,NW,IAER) = TAUAERO(K,IAER) * QXVAER(K,NW,IAER)
131	end do ! levels
132	end do
133	end do
134	do NW=1,L_NSPECTV
135	do K=2,L_LEVELS
136	TRAY(K,NW) = TAURAY(NW) * DPR(K)
137	end do ! levels
138	end do
139
140	! we ignore K=1...
141	do K=2,L_LEVELS
142
143	do NW=1,L_NSPECTV
144
145	TRAYAER = TRAY(K,NW)
146	! TRAYAER is Tau RAYleigh scattering, plus AERosol opacity
147	do iaer=1,naerkind
148	TRAYAER = TRAYAER + TAEROS(K,NW,IAER)
149	end do
150
151	DCONT = 0.0 ! continuum absorption
152
153	if(continuum.and.(.not.graybody).and.callgasvis)then
154	! include continua if necessary
155	wn_cont = dble(wnov(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.0
166	if(igas.eq.igas_N2)then
167
168	interm = indv(nw,igas,igas)
169	! call interpolateN2N2(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
170	indv(nw,igas,igas) = interm
171	! only goes to 500 cm^-1, so unless we're around a cold brown dwarf, this is irrelevant in the visible
172
173	elseif(igas.eq.igas_H2)then
174
175	! first do self-induced absorption
176	interm = indv(nw,igas,igas)
177	call interpolateH2H2(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
178	indv(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.0
184	if(jgas.eq.igas_N2)then
185	interm = indv(nw,igas,jgas)
186	call interpolateN2H2(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
187	indv(nw,igas,jgas) = interm
188	! should be irrelevant in the visible
189	elseif(jgas.eq.igas_He)then
190	interm = indv(nw,igas,jgas)
191	call interpolateH2He(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
192	indv(nw,igas,jgas) = interm
193	endif
194	dtemp = dtemp + dtempc
195	enddo
196
197	elseif(igas.eq.igas_H2O.and.T_cont.gt.200.0)then
198
199	p_air = dble(PMID(k)*scalep) - p_cont ! note assumes background is air!
200	if(H2Ocont_simple)then
201	call interpolateH2Ocont_PPC(wn_cont,T_cont,p_cont,p_air,dtemp,.false.)
202	else
203	interm = indv(nw,igas,igas)
204	call interpolateH2Ocont_CKD(wn_cont,T_cont,p_cont,p_air,dtemp,.false.,interm)
205	indv(nw,igas,igas) = interm
206	endif
207
208	endif
209
210	DCONT = DCONT + dtemp
211
212	enddo
213
214	DCONT = DCONT*dz(k)
215
216	endif
217
218	do ng=1,L_NGAUSS-1
219
220	! Now compute TAUGAS
221
222	! Interpolate between water mixing ratios
223	! WRATIO = 0.0 if the requested water amount is equal to, or outside the
224	! the water data range
225
226	if(L_REFVAR.eq.1)then ! added by RW for special no variable case
227	KCOEF(1) = GASV(MT(K),MP(K),1,NW,NG)
228	KCOEF(2) = GASV(MT(K),MP(K)+1,1,NW,NG)
229	KCOEF(3) = GASV(MT(K)+1,MP(K)+1,1,NW,NG)
230	KCOEF(4) = GASV(MT(K)+1,MP(K),1,NW,NG)
231	else
232
233	KCOEF(1) = GASV(MT(K),MP(K),NVAR(K),NW,NG) + WRATIO(K)* &
234	(GASV(MT(K),MP(K),NVAR(K)+1,NW,NG) - &
235	GASV(MT(K),MP(K),NVAR(K),NW,NG))
236
237	KCOEF(2) = GASV(MT(K),MP(K)+1,NVAR(K),NW,NG) + WRATIO(K)* &
238	(GASV(MT(K),MP(K)+1,NVAR(K)+1,NW,NG) - &
239	GASV(MT(K),MP(K)+1,NVAR(K),NW,NG))
240
241	KCOEF(3) = GASV(MT(K)+1,MP(K)+1,NVAR(K),NW,NG) + WRATIO(K)*&
242	(GASV(MT(K)+1,MP(K)+1,NVAR(K)+1,NW,NG) - &
243	GASV(MT(K)+1,MP(K)+1,NVAR(K),NW,NG))
244
245	KCOEF(4) = GASV(MT(K)+1,MP(K),NVAR(K),NW,NG) + WRATIO(K)* &
246	(GASV(MT(K)+1,MP(K),NVAR(K)+1,NW,NG) - &
247	GASV(MT(K)+1,MP(K),NVAR(K),NW,NG))
248
249	endif
250
251	! Interpolate the gaseous k-coefficients to the requested T,P values
252
253	ANS = LKCOEF(K,1)KCOEF(1) + LKCOEF(K,2)KCOEF(2) + &
254	LKCOEF(K,3)KCOEF(3) + LKCOEF(K,4)KCOEF(4)
255
256	TAUGAS = U(k)*ANS
257
258	TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS + DCONT
259	DTAUKV(K,nw,ng) = TAUGAS &
260	+ TRAYAER & ! TRAYAER includes all scattering contributions
261	+ DCONT ! For parameterized continuum aborption
262
263	end do
264
265	! Now fill in the "clear" part of the spectrum (NG = L_NGAUSS),
266	! which holds continuum opacity only
267
268	NG = L_NGAUSS
269	DTAUKV(K,nw,ng) = TRAY(K,NW) + DCONT ! For parameterized continuum absorption
270
271	do iaer=1,naerkind
272	DTAUKV(K,nw,ng) = DTAUKV(K,nw,ng) + TAEROS(K,NW,IAER)
273	end do ! a bug was here!
274
275	end do
276	end do
277
278
279	!=======================================================================
280	! Now the full treatment for the layers, where besides the opacity
281	! we need to calculate the scattering albedo and asymmetry factors
282
283	do iaer=1,naerkind
284	DO NW=1,L_NSPECTV
285	DO K=2,L_LEVELS ! AS: shouldn't this be L_LEVELS+1 ? (see optci)
286	TAUAEROLK(K,NW,IAER) = TAUAERO(K,IAER) * QSVAER(K,NW,IAER)
287	ENDDO
288	ENDDO
289	end do
290
291	DO NW=1,L_NSPECTV
292	DO L=1,L_NLAYRAD-1
293	K = 2*L+1
294	atemp(L,NW) = SUM(GVAER(K,NW,1:naerkind) * TAUAEROLK(K,NW,1:naerkind))+SUM(GVAER(K+1,NW,1:naerkind) * TAUAEROLK(K+1,NW,1:naerkind))
295	btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind)) + SUM(TAUAEROLK(K+1,NW,1:naerkind))
296	ctemp(L,NW) = btemp(L,NW) + 0.9999*(TRAY(K,NW) + TRAY(K+1,NW))
297	btemp(L,NW) = btemp(L,NW) + TRAY(K,NW) + TRAY(K+1,NW)
298	COSBV(L,NW,1:L_NGAUSS) = atemp(L,NW)/btemp(L,NW)
299	END DO ! L vertical loop
300
301	!last level
302	L = L_NLAYRAD
303	K = 2*L+1
304	atemp(L,NW) = SUM(GVAER(K,NW,1:naerkind) * TAUAEROLK(K,NW,1:naerkind))
305	btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind))
306	ctemp(L,NW) = btemp(L,NW) + 0.9999*TRAY(K,NW)
307	btemp(L,NW) = btemp(L,NW) + TRAY(K,NW)
308	COSBV(L,NW,1:L_NGAUSS) = atemp(L,NW)/btemp(L,NW)
309
310
311	END DO ! NW spectral loop
312
313	DO NG=1,L_NGAUSS
314	DO NW=1,L_NSPECTV
315	DO L=1,L_NLAYRAD-1
316
317	K = 2*L+1
318	DTAUV(L,nw,ng) = DTAUKV(K,NW,NG) + DTAUKV(K+1,NW,NG)
319	WBARV(L,nw,ng) = ctemp(L,NW) / DTAUV(L,nw,ng)
320
321	END DO ! L vertical loop
322
323	! No vertical averaging on bottom layer
324
325	L = L_NLAYRAD
326	K = 2*L+1
327	DTAUV(L,nw,ng) = DTAUKV(K,NW,NG)
328
329	WBARV(L,NW,NG) = ctemp(L,NW) / DTAUV(L,NW,NG)
330	END DO ! NW spectral loop
331	END DO ! NG Gauss loop
332
333	! Total extinction optical depths
334
335	DO NG=1,L_NGAUSS ! full gauss loop
336	DO NW=1,L_NSPECTV
337	TAUV(1,NW,NG)=0.0D0
338	DO L=1,L_NLAYRAD
339	TAUV(L+1,NW,NG)=TAUV(L,NW,NG)+DTAUV(L,NW,NG)
340	END DO
341
342	TAUCUMV(1,NW,NG)=0.0D0
343	DO K=2,L_LEVELS
344	TAUCUMV(K,NW,NG)=TAUCUMV(K-1,NW,NG)+DTAUKV(K,NW,NG)
345	END DO
346	END DO
347	END DO ! end full gauss loop
348
349
350	return
351
352
353	end subroutine optcv

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

Download in other formats:

Original Format