Context Navigation

optcv.F90 @ 486

Last change on this file since 486 was 471, checked in by aslmd, 13 years ago

LMDZ.GENERIC

13/12/2011 == AS

Same spirit as previous commit, but for ngasmx which is now read in gases.def -- before arrays w/ dim ngasmx are allocated dynamically
Allocation is done in su_gases.F90 which is called in inifis
Outside su_gases.F90, very few modifications to the code : the new module "gases_h.F90" simply replaces the old common "gases.h" !
Compiles fine. Tested with debugging options through pgdbg. Runs fine. Exact same results in Early Mars test case.

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))then
160	DCONT=0.0
161	endif
162
163
164	do NG=1,L_NGAUSS-1
165
166	!=======================================================================
167	! Now compute TAUGAS
168	! Interpolate between water mixing ratios
169	! WRATIO = 0.0 if the requested water amount is equal to, or outside the
170	! the water data range
171
172	if (L_REFVAR.eq.1)then ! added by RW for special no variable case
173	KCOEF(1) = GASV(MT(K),MP(K),1,NW,NG)
174	KCOEF(2) = GASV(MT(K),MP(K)+1,1,NW,NG)
175	KCOEF(3) = GASV(MT(K)+1,MP(K)+1,1,NW,NG)
176	KCOEF(4) = GASV(MT(K)+1,MP(K),1,NW,NG)
177	else
178	KCOEF(1) = GASV(MT(K),MP(K),NVAR(K),NW,NG) + WRATIO(K)* &
179	(GASV(MT(K),MP(K),NVAR(K)+1,NW,NG) - &
180	GASV(MT(K),MP(K),NVAR(K),NW,NG))
181
182	KCOEF(2) = GASV(MT(K),MP(K)+1,NVAR(K),NW,NG) + WRATIO(K)* &
183	(GASV(MT(K),MP(K)+1,NVAR(K)+1,NW,NG) - &
184	GASV(MT(K),MP(K)+1,NVAR(K),NW,NG))
185
186	KCOEF(3) = GASV(MT(K)+1,MP(K)+1,NVAR(K),NW,NG) + WRATIO(K)*&
187	(GASV(MT(K)+1,MP(K)+1,NVAR(K)+1,NW,NG) - &
188	GASV(MT(K)+1,MP(K)+1,NVAR(K),NW,NG))
189
190	KCOEF(4) = GASV(MT(K)+1,MP(K),NVAR(K),NW,NG) + WRATIO(K)* &
191	(GASV(MT(K)+1,MP(K),NVAR(K)+1,NW,NG) - &
192	GASV(MT(K)+1,MP(K),NVAR(K),NW,NG))
193	endif
194
195	! Interpolate the gaseous k-coefficients to the requested T,P values
196
197	ANS = LKCOEF(K,1)KCOEF(1) + LKCOEF(K,2)KCOEF(2) + &
198	LKCOEF(K,3)KCOEF(3) + LKCOEF(K,4)KCOEF(4)
199
200	TAUGAS = U(k)*ANS
201
202	if(graybody)then
203	TAUGAS = 0.0
204	DCONT = 0.0
205	endif
206
207
208
209	!TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS
210	TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS + DCONT
211	DTAUKV(K,nw,ng) = TAUGAS + TRAYAER & ! TRAYAER includes all scattering contributions
212	+ DCONT ! for continuum absorption
213
214	end do
215
216
217	! Now fill in the "clear" part of the spectrum (NG = L_NGAUSS),
218	! which holds continuum opacity only
219
220	NG = L_NGAUSS
221	DTAUKV(K,nw,ng) = TRAY(K,NW) + DCONT ! For parameterized continuum absorption
222	do iaer=1,naerkind
223	DTAUKV(K,nw,ng) = DTAUKV(K,nw,ng) + TAEROS(K,NW,IAER)
224	! & + DCONT ! For parameterized continuum absorption
225	end do ! a bug was here!
226
227	end do
228	end do
229
230
231	!=======================================================================
232	! Now the full treatment for the layers, where besides the opacity
233	! we need to calculate the scattering albedo and asymmetry factors
234
235	DO NW=1,L_NSPECTV
236	DO K=2,L_LEVELS
237	do iaer=1,naerkind
238	TAUAEROLK(K,NW,IAER) = TAUAERO(K,IAER) * QSVAER(K,NW,IAER)
239	end do
240	ENDDO
241	ENDDO
242
243
244	DO NW=1,L_NSPECTV
245	DO NG=1,L_NGAUSS
246	DO L=1,L_NLAYRAD-1
247	K = 2*L+1
248
249	DTAUV(L,nw,ng) = DTAUKV(K,NW,NG)+DTAUKV(K+1,NW,NG)
250
251	atemp=0.
252	btemp=TRAY(K,NW) + TRAY(K+1,NW)
253	ctemp=0.9999*(TRAY(K,NW) + TRAY(K+1,NW))
254	do iaer=1,naerkind
255	atemp = atemp + &
256	GVAER(K,NW,IAER) * TAUAEROLK(K,NW,IAER) + &
257	GVAER(K+1,NW,IAER) * TAUAEROLK(K+1,NW,IAER)
258	btemp = btemp + &
259	TAUAEROLK(K,NW,IAER) + TAUAEROLK(K+1,NW,IAER)
260	ctemp = ctemp + &
261	TAUAEROLK(K,NW,IAER) + TAUAEROLK(K+1,NW,IAER)
262	end do
263
264	COSBV(L,NW,NG) = atemp/btemp
265	WBARV(L,nw,ng) = ctemp/DTAUV(L,nw,ng)
266
267	END DO
268
269	! No vertical averaging on bottom layer
270
271	L = L_NLAYRAD
272	K = 2*L+1
273	DTAUV(L,nw,ng) = DTAUKV(K,NW,NG)
274
275	atemp=0.
276	btemp=TRAY(K,NW)
277	ctemp=0.9999*TRAY(K,NW)
278	do iaer=1,naerkind
279	atemp = atemp + GVAER(K,NW,IAER) * TAUAEROLK(K,NW,IAER)
280	btemp = btemp + TAUAEROLK(K,NW,IAER)
281	ctemp = ctemp + TAUAEROLK(K,NW,IAER)
282	end do
283	COSBV(L,NW,NG) = atemp/btemp
284	WBARV(L,nw,ng) = ctemp/DTAUV(L,nw,ng)
285
286	END DO ! NG gauss point loop
287	END DO ! NW spectral loop
288
289
290
291	! Total extinction optical depths
292
293	DO NW=1,L_NSPECTV
294	DO NG=1,L_NGAUSS ! full gauss loop
295	TAUV(1,NW,NG)=0.0D0
296	DO L=1,L_NLAYRAD
297	TAUV(L+1,NW,NG)=TAUV(L,NW,NG)+DTAUV(L,NW,NG)
298	END DO
299
300	TAUCUMV(1,NW,NG)=0.0D0
301	DO K=2,L_LEVELS
302	TAUCUMV(K,NW,NG)=TAUCUMV(K-1,NW,NG)+DTAUKV(K,NW,NG)
303	END DO
304	END DO ! end full gauss loop
305	END DO
306
307
308	RETURN
309	END SUBROUTINE OPTCV

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

Download in other formats:

Original Format