Context Navigation

optcv.F90 @ 537

Last change on this file since 537 was 526, checked in by jleconte, 13 years ago

13/02/2012 == JL + AS

All outputs are now in netCDF format. Even in 1D (No more G1D)
Clean up of the call to callcorrk when CLFvarying=true
Corrects a bug in writediagspecIR/VI. Output are now in W/m2/cm-1 as a function of the wavenumber in cm-1
Enable writediagspecIR/V to work in the CLFvarying=true case (output now done in Physiq after writediagfi)
Add a simple treatment for the supersaturation of CO2 (see forget et al 2012)
corrects a small bug when no clouds are present in aeropacity

Property svn:executable set to *

File size: 10.3 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
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 "callkeys.h"
38	#include "comcstfi.h"
39
40	real*8 DTAUV(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
41	real*8 DTAUKV(L_LEVELS+1,L_NSPECTV,L_NGAUSS)
42	real*8 TAUV(L_NLEVRAD,L_NSPECTV,L_NGAUSS)
43	real*8 TAUCUMV(L_LEVELS,L_NSPECTV,L_NGAUSS)
44	real*8 PLEV(L_LEVELS)
45	real*8 TMID(L_LEVELS), PMID(L_LEVELS)
46	real*8 COSBV(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
47	real*8 WBARV(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
48	real*8 TAURAY(L_NSPECTV)
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, NG1(L_NSPECTV), LK, IAER
59	integer MT(L_LEVELS), MP(L_LEVELS), NP(L_LEVELS)
60	real*8 ANS, TAUGAS
61	real*8 TRAY(L_LEVELS,L_NSPECTV)
62	real*8 DPR(L_LEVELS), U(L_LEVELS)
63	real*8 LCOEF(4), LKCOEF(L_LEVELS,4)
64
65	real*8 taugsurf(L_NSPECTV,L_NGAUSS-1), TRAYAER
66
67	! Variable species mixing ratio variables
68	real*8 QVAR(L_LEVELS), WRATIO(L_LEVELS), MUVAR(L_LEVELS)
69	real*8 KCOEF(4)
70	integer NVAR(L_LEVELS)
71
72	! temporary variables for multiple aerosol calculation
73	real*8 atemp, btemp, ctemp
74
75	! variables for k in units m^-1
76	double precision wn_cont, p_cont, p_air, T_cont, dtemp
77	real*8 dz(L_LEVELS), DCONT
78
79	integer igas
80
81	!=======================================================================
82	! Determine the total gas opacity throughout the column, for each
83	! spectral interval, NW, and each Gauss point, NG.
84	! Calculate the continuum opacities, i.e., those that do not depend on
85	! NG, the Gauss index.
86
87	taugsurf(:,:) = 0.0
88	dpr(:) = 0.0
89	lkcoef(:,:) = 0.0
90
91	do K=2,L_LEVELS
92	DPR(k) = PLEV(K)-PLEV(K-1)
93
94	! if we have continuum opacities, we need dz
95	if(kastprof)then
96	dz(k) = dpr(k)(8314.5/muvar(k))TMID(K)/(g*PMID(K))
97	U(k) = (Cmkmugaz/(muvar(k)))DPR(k)
98	else
99	dz(k) = dpr(k)RTMID(K)/(g*PMID(K))
100	U(k) = Cmk*DPR(k) ! only Cmk line in optci.F
101	endif
102
103	call tpindex(PMID(K),TMID(K),QVAR(K),pfgasref,tgasref,WREFVAR, &
104	LCOEF,MT(K),MP(K),NVAR(K),WRATIO(K))
105
106	do LK=1,4
107	LKCOEF(K,LK) = LCOEF(LK)
108	end do
109
110	DO NW=1,L_NSPECTV
111	TRAY(K,NW) = TAURAY(NW) * DPR(K)
112
113	do iaer=1,naerkind
114	TAEROS(K,NW,IAER) = TAUAERO(K,IAER) * QXVAER(K,NW,IAER)
115	end do
116	!
117
118	END DO
119	end do
120
121	! TRAYAER is Tau RAYleigh scattering, plus AERosol opacity
122
123	! we ignore K=1...
124	do K=2,L_LEVELS
125	do NW=1,L_NSPECTV
126
127	TRAYAER = TRAY(K,NW)
128	do iaer=1,naerkind
129	TRAYAER = TRAYAER + TAEROS(K,NW,IAER)
130	end do
131
132	DCONT = 0.0 ! continuum absorption
133
134	! include continua if necessary
135	wn_cont = dble(wnov(nw))
136	T_cont = dble(TMID(k))
137	do igas=1,ngasmx
138
139	if(gfrac(igas).eq.-1)then ! variable
140	p_cont = dble(PMID(k)scalepQVAR(K)) ! qvar = mol/mol
141	else
142	p_cont = dble(PMID(k)scalepgfrac(igas)*(1.-QVAR(k)))
143	endif
144
145	dtemp=0.0
146	if(gnom(igas).eq.'H2_')then
147	call interpolateH2H2(wn_cont,T_cont,p_cont,dtemp,.false.)
148	elseif(gnom(igas).eq.'H2O'.and.T_cont.gt.200.0)then
149	p_air = dble(PMID(k)*scalep) - p_cont ! note assumes air!!
150	call interpolateH2Ocont(wn_cont,T_cont,p_cont,p_air,dtemp,.false.)
151	endif
152
153	DCONT = DCONT + dtemp
154
155	enddo
156
157	DCONT = DCONT*dz(k)
158
159	if(.not.(callgasvis.and.Continuum))then
160	DCONT=0.0
161	endif
162
163	do NG=1,L_NGAUSS-1
164
165	!=======================================================================
166	! Now compute TAUGAS
167	! Interpolate between water mixing ratios
168	! WRATIO = 0.0 if the requested water amount is equal to, or outside the
169	! the water data range
170
171	if (L_REFVAR.eq.1)then ! added by RW for special no variable case
172	KCOEF(1) = GASV(MT(K),MP(K),1,NW,NG)
173	KCOEF(2) = GASV(MT(K),MP(K)+1,1,NW,NG)
174	KCOEF(3) = GASV(MT(K)+1,MP(K)+1,1,NW,NG)
175	KCOEF(4) = GASV(MT(K)+1,MP(K),1,NW,NG)
176	else
177	KCOEF(1) = GASV(MT(K),MP(K),NVAR(K),NW,NG) + WRATIO(K)* &
178	(GASV(MT(K),MP(K),NVAR(K)+1,NW,NG) - &
179	GASV(MT(K),MP(K),NVAR(K),NW,NG))
180
181	KCOEF(2) = GASV(MT(K),MP(K)+1,NVAR(K),NW,NG) + WRATIO(K)* &
182	(GASV(MT(K),MP(K)+1,NVAR(K)+1,NW,NG) - &
183	GASV(MT(K),MP(K)+1,NVAR(K),NW,NG))
184
185	KCOEF(3) = GASV(MT(K)+1,MP(K)+1,NVAR(K),NW,NG) + WRATIO(K)*&
186	(GASV(MT(K)+1,MP(K)+1,NVAR(K)+1,NW,NG) - &
187	GASV(MT(K)+1,MP(K)+1,NVAR(K),NW,NG))
188
189	KCOEF(4) = GASV(MT(K)+1,MP(K),NVAR(K),NW,NG) + WRATIO(K)* &
190	(GASV(MT(K)+1,MP(K),NVAR(K)+1,NW,NG) - &
191	GASV(MT(K)+1,MP(K),NVAR(K),NW,NG))
192	endif
193
194	! Interpolate the gaseous k-coefficients to the requested T,P values
195
196	ANS = LKCOEF(K,1)KCOEF(1) + LKCOEF(K,2)KCOEF(2) + &
197	LKCOEF(K,3)KCOEF(3) + LKCOEF(K,4)KCOEF(4)
198
199	TAUGAS = U(k)*ANS
200
201	if(graybody)then
202	TAUGAS = 0.0
203	DCONT = 0.0
204	endif
205
206
207
208	!TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS
209	TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS + DCONT
210	DTAUKV(K,nw,ng) = TAUGAS + TRAYAER & ! TRAYAER includes all scattering contributions
211	+ DCONT ! for continuum absorption
212
213	end do
214
215
216	! Now fill in the "clear" part of the spectrum (NG = L_NGAUSS),
217	! which holds continuum opacity only
218
219	NG = L_NGAUSS
220	DTAUKV(K,nw,ng) = TRAY(K,NW) + DCONT ! For parameterized continuum absorption
221	do iaer=1,naerkind
222	DTAUKV(K,nw,ng) = DTAUKV(K,nw,ng) + TAEROS(K,NW,IAER)
223	! & + DCONT ! For parameterized continuum absorption
224	end do ! a bug was here!
225
226	end do
227	end do
228
229
230	!=======================================================================
231	! Now the full treatment for the layers, where besides the opacity
232	! we need to calculate the scattering albedo and asymmetry factors
233
234	DO NW=1,L_NSPECTV
235	DO K=2,L_LEVELS
236	do iaer=1,naerkind
237	TAUAEROLK(K,NW,IAER) = TAUAERO(K,IAER) * QSVAER(K,NW,IAER)
238	end do
239	ENDDO
240	ENDDO
241
242
243	DO NW=1,L_NSPECTV
244	DO NG=1,L_NGAUSS
245	DO L=1,L_NLAYRAD-1
246	K = 2*L+1
247
248	DTAUV(L,nw,ng) = DTAUKV(K,NW,NG)+DTAUKV(K+1,NW,NG)
249
250	atemp=0.
251	btemp=TRAY(K,NW) + TRAY(K+1,NW)
252	ctemp=0.9999*(TRAY(K,NW) + TRAY(K+1,NW))
253	do iaer=1,naerkind
254	atemp = atemp + &
255	GVAER(K,NW,IAER) * TAUAEROLK(K,NW,IAER) + &
256	GVAER(K+1,NW,IAER) * TAUAEROLK(K+1,NW,IAER)
257	btemp = btemp + &
258	TAUAEROLK(K,NW,IAER) + TAUAEROLK(K+1,NW,IAER)
259	ctemp = ctemp + &
260	TAUAEROLK(K,NW,IAER) + TAUAEROLK(K+1,NW,IAER)
261	end do
262
263	COSBV(L,NW,NG) = atemp/btemp
264	WBARV(L,nw,ng) = ctemp/DTAUV(L,nw,ng)
265
266	END DO
267
268	! No vertical averaging on bottom layer
269
270	L = L_NLAYRAD
271	K = 2*L+1
272	DTAUV(L,nw,ng) = DTAUKV(K,NW,NG)
273
274	atemp=0.
275	btemp=TRAY(K,NW)
276	ctemp=0.9999*TRAY(K,NW)
277	do iaer=1,naerkind
278	atemp = atemp + GVAER(K,NW,IAER) * TAUAEROLK(K,NW,IAER)
279	btemp = btemp + TAUAEROLK(K,NW,IAER)
280	ctemp = ctemp + TAUAEROLK(K,NW,IAER)
281	end do
282	COSBV(L,NW,NG) = atemp/btemp
283	WBARV(L,nw,ng) = ctemp/DTAUV(L,nw,ng)
284
285	END DO ! NG gauss point loop
286	END DO ! NW spectral loop
287
288
289
290	! Total extinction optical depths
291
292	DO NW=1,L_NSPECTV
293	DO NG=1,L_NGAUSS ! full gauss loop
294	TAUV(1,NW,NG)=0.0D0
295	DO L=1,L_NLAYRAD
296	TAUV(L+1,NW,NG)=TAUV(L,NW,NG)+DTAUV(L,NW,NG)
297	END DO
298
299	TAUCUMV(1,NW,NG)=0.0D0
300	DO K=2,L_LEVELS
301	TAUCUMV(K,NW,NG)=TAUCUMV(K-1,NW,NG)+DTAUKV(K,NW,NG)
302	END DO
303	END DO ! end full gauss loop
304	END DO
305
306
307	RETURN
308	END SUBROUTINE OPTCV

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

Download in other formats:

Original Format