Context Navigation

optcv.F90 @ 873

Last change on this file since 873 was 873, checked in by aslmd, 12 years ago

LMDZ.GENERIC (and LMDZ.UNIVERSAL)

Optimized calculations for continuum (done for H2 and He, to be done for others)
- new common bilinear interpolation routine (bilinearbig)
- optimization: only one calculation is actually needed

to find indexes of wavelength for bilinear interpolation
... because this will not change with level and integration step!

optimization: use while loop in bilinearbig
completely similar results obtained (test case for a gas giant, many simulated days)

NB: those changes really improve gcm speed (factor 2.2 for whole model!)

continuum was very expensive, now very cheap
--> e.g. 1 day, 25 dyn ts, 5 phys ts
--> before: 243 seconds (including 120 seconds for continuum bilinear interpolation)
--> after: 108 seconds

Corrected a bug: Continuum in inifis instead of continuum

... until now, most users (unbeknownst to them) were running with the continuum by default!

Cosmetic changes in optcv (mostly spaces and line breaks)

... so that comparisons with optci are easy e.g. through vimdiff

Property svn:executable set to *

File size: 10.9 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
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, btemp, ctemp
79
80	! variables for k in units m^-1
81	real*8 dz(L_LEVELS)
82
83	integer igas, jgas
84
85	integer interm
86
87	!! AS: to save time in computing continuum (see bilinearbig)
88	IF (.not.ALLOCATED(indv)) THEN
89	ALLOCATE(indv(L_NSPECTV,ngasmx))
90	indv = -9999 ! this initial value means "to be calculated"
91	ENDIF
92
93	!=======================================================================
94	! Determine the total gas opacity throughout the column, for each
95	! spectral interval, NW, and each Gauss point, NG.
96	! Calculate the continuum opacities, i.e., those that do not depend on
97	! NG, the Gauss index.
98
99	taugsurf(:,:) = 0.0
100	dpr(:) = 0.0
101	lkcoef(:,:) = 0.0
102
103	do K=2,L_LEVELS
104	DPR(k) = PLEV(K)-PLEV(K-1)
105
106	! if we have continuum opacities, we need dz
107	if(kastprof)then
108	dz(k) = dpr(k)(1000.08.3145/muvar(k))TMID(K)/(gPMID(K))
109	U(k) = (Cmkmugaz/(muvar(k)))DPR(k)
110	else
111	dz(k) = dpr(k)RTMID(K)/(g*PMID(K))
112	U(k) = Cmk*DPR(k) ! only Cmk line in optci.F
113	endif
114
115	call tpindex(PMID(K),TMID(K),QVAR(K),pfgasref,tgasref,WREFVAR, &
116	LCOEF,MT(K),MP(K),NVAR(K),WRATIO(K))
117
118	do LK=1,4
119	LKCOEF(K,LK) = LCOEF(LK)
120	end do
121
122
123	DO NW=1,L_NSPECTV
124	do iaer=1,naerkind
125	TAEROS(K,NW,IAER) = TAUAERO(K,IAER) * QXVAER(K,NW,IAER)
126	end do
127	TRAY(K,NW) = TAURAY(NW) * DPR(K)
128	END DO
129	end do ! levels
130
131	! we ignore K=1...
132	do K=2,L_LEVELS
133
134	do NW=1,L_NSPECTV
135
136	TRAYAER = TRAY(K,NW)
137	! TRAYAER is Tau RAYleigh scattering, plus AERosol opacity
138	do iaer=1,naerkind
139	TRAYAER = TRAYAER + TAEROS(K,NW,IAER)
140	end do
141
142	DCONT = 0.0 ! continuum absorption
143
144	if(continuum.and.(.not.graybody).and.callgasvis)then
145	! include continua if necessary
146	wn_cont = dble(wnov(nw))
147	T_cont = dble(TMID(k))
148	do igas=1,ngasmx
149
150	if(gfrac(igas).eq.-1)then ! variable
151	p_cont = dble(PMID(k)scalepQVAR(k)) ! qvar = mol/mol
152	else
153	p_cont = dble(PMID(k)scalepgfrac(igas)*(1.-QVAR(k)))
154	endif
155
156	dtemp=0.0
157	if(igas.eq.igas_N2)then
158
159	!call interpolateN2N2(wn_cont,T_cont,p_cont,dtemp,.false.)
160	! only goes to 500 cm^-1, so unless we're around a cold brown dwarf, this is irrelevant in the visible
161
162	elseif(igas.eq.igas_H2)then
163
164	! first do self-induced absorption
165	interm = indv(nw,igas)
166	call interpolateH2H2(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
167	indv(nw,igas) = interm
168
169	! then cross-interactions with other gases
170	do jgas=1,ngasmx
171	p_cross = dble(PMID(k)scalepgfrac(jgas)*(1.-QVAR(k)))
172	dtempc = 0.0
173	if(jgas.eq.igas_N2)then
174	call interpolateN2H2(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.)
175	! should be irrelevant in the visible
176	elseif(jgas.eq.igas_He)then
177	interm = indv(nw,jgas)
178	call interpolateH2He(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
179	indv(nw,jgas) = interm
180	endif
181	dtemp = dtemp + dtempc
182	enddo
183
184	elseif(igas.eq.igas_H2O.and.T_cont.gt.200.0)then
185
186	p_air = dble(PMID(k)*scalep) - p_cont ! note assumes background is air!
187	if(H2Ocont_simple)then
188	call interpolateH2Ocont_PPC(wn_cont,T_cont,p_cont,p_air,dtemp,.false.)
189	else
190	call interpolateH2Ocont_CKD(wn_cont,T_cont,p_cont,p_air,dtemp,.false.)
191	endif
192
193	endif
194
195	DCONT = DCONT + dtemp
196
197	enddo
198
199	DCONT = DCONT*dz(k)
200
201	endif
202
203	do ng=1,L_NGAUSS-1
204
205	! Now compute TAUGAS
206
207	! Interpolate between water mixing ratios
208	! WRATIO = 0.0 if the requested water amount is equal to, or outside the
209	! the water data range
210
211	if(L_REFVAR.eq.1)then ! added by RW for special no variable case
212	KCOEF(1) = GASV(MT(K),MP(K),1,NW,NG)
213	KCOEF(2) = GASV(MT(K),MP(K)+1,1,NW,NG)
214	KCOEF(3) = GASV(MT(K)+1,MP(K)+1,1,NW,NG)
215	KCOEF(4) = GASV(MT(K)+1,MP(K),1,NW,NG)
216	else
217
218	KCOEF(1) = GASV(MT(K),MP(K),NVAR(K),NW,NG) + WRATIO(K)* &
219	(GASV(MT(K),MP(K),NVAR(K)+1,NW,NG) - &
220	GASV(MT(K),MP(K),NVAR(K),NW,NG))
221
222	KCOEF(2) = GASV(MT(K),MP(K)+1,NVAR(K),NW,NG) + WRATIO(K)* &
223	(GASV(MT(K),MP(K)+1,NVAR(K)+1,NW,NG) - &
224	GASV(MT(K),MP(K)+1,NVAR(K),NW,NG))
225
226	KCOEF(3) = GASV(MT(K)+1,MP(K)+1,NVAR(K),NW,NG) + WRATIO(K)*&
227	(GASV(MT(K)+1,MP(K)+1,NVAR(K)+1,NW,NG) - &
228	GASV(MT(K)+1,MP(K)+1,NVAR(K),NW,NG))
229
230	KCOEF(4) = GASV(MT(K)+1,MP(K),NVAR(K),NW,NG) + WRATIO(K)* &
231	(GASV(MT(K)+1,MP(K),NVAR(K)+1,NW,NG) - &
232	GASV(MT(K)+1,MP(K),NVAR(K),NW,NG))
233
234	endif
235
236	! Interpolate the gaseous k-coefficients to the requested T,P values
237
238	ANS = LKCOEF(K,1)KCOEF(1) + LKCOEF(K,2)KCOEF(2) + &
239	LKCOEF(K,3)KCOEF(3) + LKCOEF(K,4)KCOEF(4)
240
241	TAUGAS = U(k)*ANS
242
243	TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS + DCONT
244	DTAUKV(K,nw,ng) = TAUGAS &
245	+ TRAYAER & ! TRAYAER includes all scattering contributions
246	+ DCONT ! For parameterized continuum aborption
247
248	end do
249
250	! Now fill in the "clear" part of the spectrum (NG = L_NGAUSS),
251	! which holds continuum opacity only
252
253	NG = L_NGAUSS
254	DTAUKV(K,nw,ng) = TRAY(K,NW) + DCONT ! For parameterized continuum absorption
255
256	do iaer=1,naerkind
257	DTAUKV(K,nw,ng) = DTAUKV(K,nw,ng) + TAEROS(K,NW,IAER)
258	end do ! a bug was here!
259
260	end do
261	end do
262
263
264	!=======================================================================
265	! Now the full treatment for the layers, where besides the opacity
266	! we need to calculate the scattering albedo and asymmetry factors
267
268	do iaer=1,naerkind
269	DO NW=1,L_NSPECTV
270	DO K=2,L_LEVELS ! AS: shouldn't this be L_LEVELS+1 ? (see optci)
271	TAUAEROLK(K,NW,IAER) = TAUAERO(K,IAER) * QSVAER(K,NW,IAER)
272	ENDDO
273	ENDDO
274	end do
275
276	DO NW=1,L_NSPECTV
277	DO NG=1,L_NGAUSS
278	DO L=1,L_NLAYRAD-1
279
280	K = 2*L+1
281	DTAUV(L,nw,ng) = DTAUKV(K,NW,NG) + DTAUKV(K+1,NW,NG)
282
283	atemp = 0.
284	btemp = TRAY(K,NW) + TRAY(K+1,NW)
285	ctemp=0.9999*(TRAY(K,NW) + TRAY(K+1,NW))
286	do iaer=1,naerkind
287	atemp = atemp + &
288	GVAER(K,NW,IAER) * TAUAEROLK(K,NW,IAER) + &
289	GVAER(K+1,NW,IAER) * TAUAEROLK(K+1,NW,IAER)
290	btemp = btemp + TAUAEROLK(K,NW,IAER) + TAUAEROLK(K+1,NW,IAER)
291	ctemp = ctemp + TAUAEROLK(K,NW,IAER) + TAUAEROLK(K+1,NW,IAER)
292	end do
293	WBARV(L,nw,ng) = ctemp / DTAUV(L,nw,ng)
294	COSBV(L,NW,NG) = atemp/btemp
295
296	END DO ! L vertical loop
297
298	! No vertical averaging on bottom layer
299
300	L = L_NLAYRAD
301	K = 2*L+1
302	DTAUV(L,nw,ng) = DTAUKV(K,NW,NG)
303
304	atemp=0.
305	btemp=TRAY(K,NW)
306	ctemp=0.9999*TRAY(K,NW)
307	do iaer=1,naerkind
308	atemp = atemp + GVAER(K,NW,IAER) * TAUAEROLK(K,NW,IAER)
309	btemp = btemp + TAUAEROLK(K,NW,IAER)
310	ctemp = ctemp + TAUAEROLK(K,NW,IAER)
311	end do
312	COSBV(L,NW,NG) = atemp/btemp
313	WBARV(L,nw,ng) = ctemp/DTAUV(L,nw,ng)
314
315	END DO ! NG Gauss loop
316	END DO ! NW spectral loop
317
318	! Total extinction optical depths
319
320	DO NW=1,L_NSPECTV
321	DO NG=1,L_NGAUSS ! full gauss loop
322	TAUV(1,NW,NG)=0.0D0
323	DO L=1,L_NLAYRAD
324	TAUV(L+1,NW,NG)=TAUV(L,NW,NG)+DTAUV(L,NW,NG)
325	END DO
326
327	TAUCUMV(1,NW,NG)=0.0D0
328	DO K=2,L_LEVELS
329	TAUCUMV(K,NW,NG)=TAUCUMV(K-1,NW,NG)+DTAUKV(K,NW,NG)
330	END DO
331	END DO ! end full gauss loop
332	END DO
333
334
335	return
336
337
338	end subroutine optcv

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

Download in other formats:

Original Format