Context Navigation

optci.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

File size: 12.9 KB

Line
1	subroutine optci(PLEV,TLEV,DTAUI,TAUCUMI, &
2	QXIAER,QSIAER,GIAER,COSBI,WBARI,TAUAERO, &
3	TMID,PMID,TAUGSURF,QVAR,MUVAR)
4
5	use radinc_h
6	use radcommon_h, only: gasi,tlimit,wrefVAR,Cmk,tgasref,pfgasref,wnoi,scalep,indi
7	use gases_h
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	#include "comcstfi.h"
31	#include "callkeys.h"
32
33
34	real*8 DTAUI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
35	real*8 DTAUKI(L_LEVELS+1,L_NSPECTI,L_NGAUSS)
36	real*8 TAUI(L_NLEVRAD,L_NSPECTI,L_NGAUSS)
37	real*8 TAUCUMI(L_LEVELS,L_NSPECTI,L_NGAUSS)
38	real*8 PLEV(L_LEVELS)
39	real*8 TLEV(L_LEVELS)
40	real*8 TMID(L_LEVELS), PMID(L_LEVELS)
41	real*8 COSBI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
42	real*8 WBARI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
43
44	! for aerosols
45	real*8 QXIAER(L_LEVELS+1,L_NSPECTI,NAERKIND)
46	real*8 QSIAER(L_LEVELS+1,L_NSPECTI,NAERKIND)
47	real*8 GIAER(L_LEVELS+1,L_NSPECTI,NAERKIND)
48	real*8 TAUAERO(L_LEVELS+1,NAERKIND)
49	real*8 TAUAEROLK(L_LEVELS+1,L_NSPECTI,NAERKIND)
50	real*8 TAEROS(L_LEVELS,L_NSPECTI,NAERKIND)
51
52	integer L, NW, NG, K, LK, IAER
53	integer MT(L_LEVELS), MP(L_LEVELS), NP(L_LEVELS)
54	real*8 ANS, TAUGAS
55	real*8 DPR(L_LEVELS), U(L_LEVELS)
56	real*8 LCOEF(4), LKCOEF(L_LEVELS,4)
57
58	real*8 taugsurf(L_NSPECTI,L_NGAUSS-1)
59	real*8 DCONT
60	double precision wn_cont, p_cont, p_air, T_cont, dtemp, dtempc
61	double precision p_cross
62
63	! variable species mixing ratio variables
64	real*8 QVAR(L_LEVELS), WRATIO(L_LEVELS), MUVAR(L_LEVELS)
65	real*8 KCOEF(4)
66	integer NVAR(L_LEVELS)
67
68	! temporary variables for multiple aerosol calculation
69	real*8 atemp, btemp
70
71	! variables for k in units m^-1
72	real*8 dz(L_LEVELS)
73	!real*8 rho !! see test below
74
75	integer igas, jgas
76
77	integer interm
78
79	!--- Kasting's CIA ----------------------------------------
80	!real*8, parameter :: Ci(L_NSPECTI)=[ &
81	! 3.8E-5, 1.2E-5, 2.8E-6, 7.6E-7, 4.5E-7, 2.3E-7, &
82	! 5.4E-7, 1.6E-6, 0.0, &
83	! 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
84	! 0.0, 4.0E-7, 4.0E-6, 1.4E-5, &
85	! 1.0E-5, 1.2E-6, 2.0E-7, 5.0E-8, 3.0E-8, 0.0 ]
86	!real*8, parameter :: Ti(L_NSPECTI)=[ -2.2, -1.9, &
87	! -1.7, -1.7, -1.7, -1.7, -1.7, -1.7, &
88	! 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
89	! -1.7,-1.7,-1.7,-1.7,-1.7,-1.7,-1.7, -1.7,0.0 ]
90	!----------------------------------------------------------
91
92	!! AS: to save time in computing continuum (see bilinearbig)
93	IF (.not.ALLOCATED(indi)) THEN
94	ALLOCATE(indi(L_NSPECTI,ngasmx))
95	indi = -9999 ! this initial value means "to be calculated"
96	ENDIF
97
98	!=======================================================================
99	! Determine the total gas opacity throughout the column, for each
100	! spectral interval, NW, and each Gauss point, NG.
101
102	taugsurf(:,:) = 0.0
103	dpr(:) = 0.0
104	lkcoef(:,:) = 0.0
105
106	do K=2,L_LEVELS
107	DPR(k) = PLEV(K)-PLEV(K-1)
108
109	!--- Kasting's CIA ----------------------------------------
110	!dz(k)=dpr(k)189.02TMID(K)/(0.03720*PMID(K))
111	! this is CO2 path length (in cm) as written by Francois
112	! delta_z = delta_p * R_specific * T / (g * P)
113	! But Kasting states that W is in units of _atmosphere_ cm
114	! So we do
115	!dz(k)=dz(k)*(PMID(K)/1013.25)
116	!dz(k)=dz(k)/100.0 ! in m for SI calc
117	!----------------------------------------------------------
118
119	! if we have continuum opacities, we need dz
120	if(kastprof)then
121	dz(k) = dpr(k)(1000.08.3145/muvar(k))TMID(K)/(gPMID(K))
122	U(k) = (Cmkmugaz/(muvar(k)))DPR(k)
123	else
124	dz(k) = dpr(k)RTMID(K)/(g*PMID(K))
125	U(k) = Cmk*DPR(k) ! only Cmk line in optci.F
126	endif
127
128	call tpindex(PMID(K),TMID(K),QVAR(K),pfgasref,tgasref,WREFVAR, &
129	LCOEF,MT(K),MP(K),NVAR(K),WRATIO(K))
130
131	do LK=1,4
132	LKCOEF(K,LK) = LCOEF(LK)
133	end do
134
135
136	DO NW=1,L_NSPECTI
137	do iaer=1,naerkind
138	TAEROS(K,NW,IAER) = TAUAERO(K,IAER) * QXIAER(K,NW,IAER)
139	end do
140	END DO
141	end do ! levels
142
143	do K=2,L_LEVELS
144
145	do NW=1,L_NSPECTI
146
147	DCONT = 0.0 ! continuum absorption
148
149	if(continuum.and.(.not.graybody))then
150	! include continua if necessary
151	wn_cont = dble(wnoi(nw))
152	T_cont = dble(TMID(k))
153	do igas=1,ngasmx
154
155	if(gfrac(igas).eq.-1)then ! variable
156	p_cont = dble(PMID(k)scalepQVAR(k)) ! qvar = mol/mol
157	else
158	p_cont = dble(PMID(k)scalepgfrac(igas)*(1.-QVAR(k)))
159	endif
160
161	dtemp=0.0
162	if(igas.eq.igas_N2)then
163
164	call interpolateN2N2(wn_cont,T_cont,p_cont,dtemp,.false.)
165
166	elseif(igas.eq.igas_H2)then
167
168	! first do self-induced absorption
169	interm = indi(nw,igas)
170	call interpolateH2H2(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
171	indi(nw,igas) = interm
172
173	! then cross-interactions with other gases
174	do jgas=1,ngasmx
175	p_cross = dble(PMID(k)scalepgfrac(jgas)*(1.-QVAR(k)))
176	dtempc = 0.0
177	if(jgas.eq.igas_N2)then
178	call interpolateN2H2(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.)
179	elseif(jgas.eq.igas_He)then
180	interm = indi(nw,jgas)
181	call interpolateH2He(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
182	indi(nw,jgas) = interm
183	endif
184	dtemp = dtemp + dtempc
185	enddo
186
187	elseif(igas.eq.igas_H2O.and.T_cont.gt.200.0)then
188
189	p_air = dble(PMID(k)*scalep) - p_cont ! note assumes background is air!
190	if(H2Ocont_simple)then
191	call interpolateH2Ocont_PPC(wn_cont,T_cont,p_cont,p_air,dtemp,.false.)
192	else
193	call interpolateH2Ocont_CKD(wn_cont,T_cont,p_cont,p_air,dtemp,.false.)
194	endif
195
196	endif
197
198	DCONT = DCONT + dtemp
199
200	enddo
201
202	! Oobleck test
203	!rho = PMID(k)scalep / (TMID(k)286.99)
204	!if(WNOI(nw).gt.300.0 .and. WNOI(nw).lt.500.0)then
205	! DCONT = rho * 0.125 * 4.6e-4
206	!elseif(WNOI(nw).gt.500.0 .and. WNOI(nw).lt.700.0)then
207	! DCONT = 1000dpr(k) 1.0 * 4.6e-4 / g
208	! DCONT = rho * 1.0 * 4.6e-4
209	!elseif(WNOI(nw).gt.700.0 .and. WNOI(nw).lt.900.0)then
210	! DCONT = rho * 0.125 * 4.6e-4
211	!endif
212
213	DCONT = DCONT*dz(k)
214
215	endif
216
217	! RW 7/3/12: already done above
218	!if(.not.Continuum)then
219	! DCONT=0.0
220	!endif
221
222	!--- Kasting's CIA ----------------------------------------
223	!DCO2 = dz(k)Ci(nw)(1.2859PMID(k)/1000.0)(TMID(k)/300.)**Ti(nw)
224	!DCO2 = 130Ci(nw)(pmid(k)/1013.25)*2(tmid(k)/300.)*Ti(nw) dz(k)
225	! these two have been verified to give the same results
226	!----------------------------------------------------------
227
228	do ng=1,L_NGAUSS-1
229
230	! Now compute TAUGAS
231
232	! Interpolate between water mixing ratios
233	! WRATIO = 0.0 if the requested water amount is equal to, or outside the
234	! the water data range
235
236	if(L_REFVAR.eq.1)then ! added by RW for special no variable case
237	KCOEF(1) = GASI(MT(K),MP(K),1,NW,NG)
238	KCOEF(2) = GASI(MT(K),MP(K)+1,1,NW,NG)
239	KCOEF(3) = GASI(MT(K)+1,MP(K)+1,1,NW,NG)
240	KCOEF(4) = GASI(MT(K)+1,MP(K),1,NW,NG)
241	else
242
243	KCOEF(1) = GASI(MT(K),MP(K),NVAR(K),NW,NG) + WRATIO(K)* &
244	(GASI(MT(K),MP(K),NVAR(K)+1,NW,NG) - &
245	GASI(MT(K),MP(K),NVAR(K),NW,NG))
246
247	KCOEF(2) = GASI(MT(K),MP(K)+1,NVAR(K),NW,NG) + WRATIO(K)* &
248	(GASI(MT(K),MP(K)+1,NVAR(K)+1,NW,NG) - &
249	GASI(MT(K),MP(K)+1,NVAR(K),NW,NG))
250
251	KCOEF(3) = GASI(MT(K)+1,MP(K)+1,NVAR(K),NW,NG) + WRATIO(K)* &
252	(GASI(MT(K)+1,MP(K)+1,NVAR(K)+1,NW,NG) - &
253	GASI(MT(K)+1,MP(K)+1,NVAR(K),NW,NG))
254
255	KCOEF(4) = GASI(MT(K)+1,MP(K),NVAR(K),NW,NG) + WRATIO(K)* &
256	(GASI(MT(K)+1,MP(K),NVAR(K)+1,NW,NG) - &
257	GASI(MT(K)+1,MP(K),NVAR(K),NW,NG))
258
259	endif
260
261	! Interpolate the gaseous k-coefficients to the requested T,P values
262
263	ANS = LKCOEF(K,1)KCOEF(1) + LKCOEF(K,2)KCOEF(2) + &
264	LKCOEF(K,3)KCOEF(3) + LKCOEF(K,4)KCOEF(4)
265
266	TAUGAS = U(k)*ANS
267
268	TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS + DCONT
269	DTAUKI(K,nw,ng) = TAUGAS &
270	+ DCONT ! For parameterized continuum absorption
271
272	do iaer=1,naerkind
273	DTAUKI(K,nw,ng) = DTAUKI(K,nw,ng) + TAEROS(K,NW,IAER)
274	end do ! a bug was here!
275
276	end do
277
278	! Now fill in the "clear" part of the spectrum (NG = L_NGAUSS),
279	! which holds continuum opacity only
280
281	NG = L_NGAUSS
282	DTAUKI(K,nw,ng) = 0.0 + DCONT ! For parameterized continuum absorption
283
284	do iaer=1,naerkind
285	DTAUKI(K,nw,ng) = DTAUKI(K,nw,ng) + TAEROS(K,NW,IAER)
286	end do ! a bug was here!
287
288	end do
289	end do
290
291
292	!=======================================================================
293	! Now the full treatment for the layers, where besides the opacity
294	! we need to calculate the scattering albedo and asymmetry factors
295
296	do iaer=1,naerkind
297	DO NW=1,L_NSPECTI
298	DO K=2,L_LEVELS+1
299	TAUAEROLK(K,NW,IAER) = TAUAERO(K,IAER)*QSIAER(K,NW,IAER)
300	ENDDO
301	ENDDO
302	end do
303
304	DO NW=1,L_NSPECTI
305	NG = L_NGAUSS
306	DO L=1,L_NLAYRAD
307
308	K = 2*L+1
309	DTAUI(L,nw,ng) = DTAUKI(K,NW,NG) + DTAUKI(K+1,NW,NG)! + 1.e-50
310
311	atemp = 0.
312	btemp = 0.
313	if(DTAUI(L,NW,NG) .GT. 1.0E-9) then
314	do iaer=1,naerkind
315	atemp = atemp + &
316	GIAER(K,NW,IAER) * TAUAEROLK(K,NW,IAER) + &
317	GIAER(K+1,NW,IAER) * TAUAEROLK(K+1,NW,IAER)
318	btemp = btemp + TAUAEROLK(K,NW,IAER) + TAUAEROLK(K+1,NW,IAER)
319	! * + 1.e-10
320	end do
321	WBARI(L,nw,ng) = btemp / DTAUI(L,NW,NG)
322	else
323	WBARI(L,nw,ng) = 0.0D0
324	DTAUI(L,NW,NG) = 1.0E-9
325	endif
326
327	if(btemp .GT. 0.0) then
328	cosbi(L,NW,NG) = atemp/btemp
329	else
330	cosbi(L,NW,NG) = 0.0D0
331	end if
332
333	END DO ! L vertical loop
334
335	! Now the other Gauss points, if needed.
336
337	DO NG=1,L_NGAUSS-1
338	IF(TAUGSURF(NW,NG) .gt. TLIMIT) THEN
339
340	DO L=1,L_NLAYRAD
341	K = 2*L+1
342	DTAUI(L,nw,ng) = DTAUKI(K,NW,NG)+DTAUKI(K+1,NW,NG)! + 1.e-50
343
344	btemp = 0.
345	if(DTAUI(L,NW,NG) .GT. 1.0E-9) then
346
347	do iaer=1,naerkind
348	btemp = btemp + TAUAEROLK(K,NW,IAER) + TAUAEROLK(K+1,NW,IAER)
349	end do
350	WBARI(L,nw,ng) = btemp / DTAUI(L,NW,NG)
351
352	else
353	WBARI(L,nw,ng) = 0.0D0
354	DTAUI(L,NW,NG) = 1.0E-9
355	endif
356
357	cosbi(L,NW,NG) = cosbi(L,NW,L_NGAUSS)
358	END DO ! L vertical loop
359	END IF
360
361	END DO ! NG Gauss loop
362	END DO ! NW spectral loop
363
364	! Total extinction optical depths
365
366	DO NW=1,L_NSPECTI
367	DO NG=1,L_NGAUSS ! full gauss loop
368	TAUI(1,NW,NG)=0.0D0
369	DO L=1,L_NLAYRAD
370	TAUI(L+1,NW,NG)=TAUI(L,NW,NG)+DTAUI(L,NW,NG)
371	END DO
372
373	TAUCUMI(1,NW,NG)=0.0D0
374	DO K=2,L_LEVELS
375	TAUCUMI(K,NW,NG)=TAUCUMI(K-1,NW,NG)+DTAUKI(K,NW,NG)
376	END DO
377	END DO ! end full gauss loop
378	END DO
379
380	! be aware when comparing with textbook results
381	! (e.g. Pierrehumbert p. 218) that
382	! taucumi does not take the <cos theta>=0.5 factor into
383	! account. It is the optical depth for a vertically
384	! ascending ray with angle theta = 0.
385
386	!open(127,file='taucum.out')
387	!do nw=1,L_NSPECTI
388	! write(127,*) taucumi(L_LEVELS,nw,L_NGAUSS)
389	!enddo
390	!close(127)
391
392	return
393
394
395	end subroutine optci
396
397
398

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

Download in other formats:

Original Format