Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

optcv.F_old @ 3700

Last change on this file since 3700 was 3175, checked in by emillour, 23 months ago
Pluto PCM: Add the old Pluto LMDZ for reference (required prior step to making an LMDZ.PLUTO using the same framework as the other physics packages). TB+EM
File size: 8.6 KB

Line
1	SUBROUTINE OPTCV(DTAUV,TAUV,TAUCUMV,PLEV,
2	& QXVAER,QSVAER,GVAER,WBARV,COSBV,
3	& TAURAY,TAUAERO,TMID,PMID,TAUGSURF,QVAR)
4
5	use radinc_h
6	use radcommon_h, only: gasv, tlimit, wrefVAR, Cmk, tgasref,
7	& pfgasref
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
24	C
25	C THIS SUBROUTINE SETS THE OPTICAL CONSTANTS IN THE VISUAL
26	C IT CALCUALTES FOR EACH LAYER, FOR EACH SPECRAL INTERVAL IN THE VISUAL
27	C LAYER: WBAR, DTAU, COSBAR
28	C LEVEL: TAU
29	C
30	C TAUV(L,NW,NG) is the cumulative optical depth at the top of radiation code
31	C layer L. NW is spectral wavelength interval, ng the Gauss point index.
32	C
33	C TLEV(L) - Temperature at the layer boundary
34	C PLEV(L) - Pressure at the layer boundary (i.e. level)
35	C GASV(NT,NPS,NW,NG) - Visual CO2 k-coefficients
36	C
37	C----------------------------------------------------------------------C
38
39	#include "callkeys.h"
40	#include "comcstfi.h"
41
42	real*8 DTAUV(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
43	real*8 DTAUKV(L_LEVELS+1,L_NSPECTV,L_NGAUSS)
44	real*8 TAUV(L_NLEVRAD,L_NSPECTV,L_NGAUSS)
45	real*8 TAUCUMV(L_LEVELS,L_NSPECTV,L_NGAUSS)
46	real*8 PLEV(L_LEVELS)
47	real*8 TMID(L_LEVELS), PMID(L_LEVELS)
48	real*8 COSBV(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
49	real*8 WBARV(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
50	real*8 TAURAY(L_NSPECTV)
51
52	! For aerosols
53	real*8 QXVAER(L_LEVELS+1,L_NSPECTV,NAERKIND)
54	real*8 QSVAER(L_LEVELS+1,L_NSPECTV,NAERKIND)
55	real*8 GVAER(L_LEVELS+1,L_NSPECTV,NAERKIND)
56	real*8 TAUAERO(L_LEVELS+1,NAERKIND)
57	real*8 TAUAEROLK(L_LEVELS+1,L_NSPECTV,NAERKIND)
58	real*8 TAEROS(L_LEVELS,L_NSPECTV,NAERKIND)
59
60	integer L, NW, NG, K, NG1(L_NSPECTV), LK, IAER
61	integer MT(L_LEVELS), MP(L_LEVELS), NP(L_LEVELS)
62	real*8 ANS, TAUGAS
63	real*8 TRAY(L_LEVELS,L_NSPECTV)
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), TRAYAER
68
69	! mixing ratio variables
70	real*8 QVAR(L_LEVELS), WRATIO(L_LEVELS)
71	real*8 KCOEF(4)
72	integer NVAR(L_LEVELS)
73
74	! temporary variables for multiple aerosol calculation
75	real*8 atemp, btemp, ctemp
76
77	! variables for k in units m^-1
78	real*8 rho, dz
79
80	!=======================================================================
81	! Determine the total gas opacity throughout the column, for each
82	! spectral interval, NW, and each Gauss point, NG.
83	! Calculate the continuum opacities, i.e., those that do not depend on
84	! NG, the Gauss index.
85
86
87	do K=2,L_LEVELS
88	DPR(K) = PLEV(K)-PLEV(K-1)
89
90
91	! rho = PLEV(K)/(R*TMID(K))
92	rho = PMID(K)/(R*TMID(K))
93	dz = -DPR(k)/(g*rho)
94
95	U(k) = Cmk*DPR(k)
96
97	call tpindex(PMID(K),TMID(K),QVAR(K),pfgasref,tgasref,WREFVAR,
98	* LCOEF,MT(K),MP(K),NVAR(K),WRATIO(K))
99
100	do LK=1,4
101	LKCOEF(K,LK) = LCOEF(LK)
102	end do
103
104	DO NW=1,L_NSPECTV
105	TRAY(K,NW) = TAURAY(NW) * DPR(K)
106
107	do iaer=1,naerkind
108	TAEROS(K,NW,IAER) = TAUAERO(K,IAER) * QXVAER(K,NW,IAER)
109	! write(20,*) 'TB17 taero',K,NW,IAER,TAEROS(K,NW,IAER)
110	end do
111	!
112
113	END DO
114	end do
115
116	! TRAYAER is Tau RAYleigh scattering, plus AERosol opacity
117
118
119	! we ignore K=1... hope this is ok...
120	do K=2,L_LEVELS
121	do NW=1,L_NSPECTV
122
123	TRAYAER = TRAY(K,NW)
124	do iaer=1,naerkind
125	TRAYAER = TRAYAER + TAEROS(K,NW,IAER)
126	end do
127
128	do NG=1,L_NGAUSS-1
129
130	!=======================================================================
131	! Now compute TAUGAS
132	! Interpolate between mixing ratios
133	! WRATIO = 0.0 if the requested amount is equal to, or outside the
134	! the data range
135
136	if (L_REFVAR.eq.1)then ! added by RW for special no variable case
137	KCOEF(1) = GASV(MT(K),MP(K),1,NW,NG)
138	KCOEF(2) = GASV(MT(K),MP(K)+1,1,NW,NG)
139	KCOEF(3) = GASV(MT(K)+1,MP(K)+1,1,NW,NG)
140	KCOEF(4) = GASV(MT(K)+1,MP(K),1,NW,NG)
141	else
142	KCOEF(1) = GASV(MT(K),MP(K),NVAR(K),NW,NG) + WRATIO(K)*
143	* (GASV(MT(K),MP(K),NVAR(K)+1,NW,NG) -
144	* GASV(MT(K),MP(K),NVAR(K),NW,NG))
145
146	KCOEF(2) = GASV(MT(K),MP(K)+1,NVAR(K),NW,NG) + WRATIO(K)*
147	* (GASV(MT(K),MP(K)+1,NVAR(K)+1,NW,NG) -
148	* GASV(MT(K),MP(K)+1,NVAR(K),NW,NG))
149
150	KCOEF(3) = GASV(MT(K)+1,MP(K)+1,NVAR(K),NW,NG) + WRATIO(K)*
151	* (GASV(MT(K)+1,MP(K)+1,NVAR(K)+1,NW,NG) -
152	* GASV(MT(K)+1,MP(K)+1,NVAR(K),NW,NG))
153
154	KCOEF(4) = GASV(MT(K)+1,MP(K),NVAR(K),NW,NG) + WRATIO(K)*
155	* (GASV(MT(K)+1,MP(K),NVAR(K)+1,NW,NG) -
156	* GASV(MT(K)+1,MP(K),NVAR(K),NW,NG))
157	endif
158
159	! Interpolate the gaseous k-coefficients to the requested T,P values
160
161	ANS = LKCOEF(K,1)KCOEF(1) + LKCOEF(K,2)KCOEF(2) +
162	* LKCOEF(K,3)KCOEF(3) + LKCOEF(K,4)KCOEF(4)
163
164	TAUGAS = U(k)*ANS
165	TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS
166	DTAUKV(K,nw,ng) = TAUGAS + TRAYAER
167	! write(21,*) 'TB17 taugas',K,NW,ng,TAUGAS
168	end do
169
170
171	! Now fill in the "clear" part of the spectrum (NG = L_NGAUSS),
172	! which holds continuum opacity only
173
174	NG = L_NGAUSS
175	DTAUKV(K,nw,ng) = TRAY(K,NW)
176	do iaer=1,naerkind
177	DTAUKV(K,nw,ng) = DTAUKV(K,nw,ng) + TAEROS(K,NW,IAER)
178	end do
179
180	end do
181	end do
182
183
184	!=======================================================================
185	! Now the full treatment for the layers, where besides the opacity
186	! we need to calculate the scattering albedo and asymmetry factors
187	TAUAEROLK(:,:,:) = 1.e-12 ! TB17
188
189	DO NW=1,L_NSPECTV
190	DO K=2,L_LEVELS
191	do iaer=1,naerkind
192	TAUAEROLK(K,NW,IAER) = TAUAERO(K,IAER) * QSVAER(K,NW,IAER)
193	TAUAEROLK(K,NW,IAER) = max(TAUAEROLK(K,NW,IAER),1.e-12) ! TB17
194	end do
195	ENDDO
196	ENDDO
197	print*, 'TBbug TAUAEROLK =', TAUAEROLK
198
199
200	DO NW=1,L_NSPECTV
201	DO NG=1,L_NGAUSS
202	DO L=1,L_NLAYRAD-1
203	K = 2*L+1
204
205	DTAUV(L,nw,ng) = DTAUKV(K,NW,NG)+DTAUKV(K+1,NW,NG)
206
207	atemp=0.
208	btemp=TRAY(K,NW) + TRAY(K+1,NW)
209	ctemp=0.9999*(TRAY(K,NW) + TRAY(K+1,NW))
210	do iaer=1,naerkind
211	atemp = atemp +
212	* GVAER(K,NW,IAER) * TAUAEROLK(K,NW,IAER) +
213	* GVAER(K+1,NW,IAER) * TAUAEROLK(K+1,NW,IAER)
214	btemp = btemp +
215	* TAUAEROLK(K,NW,IAER) + TAUAEROLK(K+1,NW,IAER)
216	ctemp = ctemp +
217	* TAUAEROLK(K,NW,IAER) + TAUAEROLK(K+1,NW,IAER)
218	end do
219
220	COSBV(L,NW,NG) = atemp/btemp
221	WBARV(L,nw,ng) = ctemp/DTAUV(L,nw,ng)
222
223	END DO
224
225	! No vertical averaging on bottom layer
226
227	L = L_NLAYRAD
228	K = 2*L+1
229	DTAUV(L,nw,ng) = DTAUKV(K,NW,NG)
230
231	atemp=0.
232	btemp=TRAY(K,NW)
233	ctemp=0.9999*TRAY(K,NW)
234	do iaer=1,naerkind
235	atemp = atemp + GVAER(K,NW,IAER) * TAUAEROLK(K,NW,IAER)
236	btemp = btemp + TAUAEROLK(K,NW,IAER)
237	ctemp = ctemp + TAUAEROLK(K,NW,IAER)
238	end do
239	COSBV(L,NW,NG) = atemp/btemp
240	WBARV(L,nw,ng) = ctemp/DTAUV(L,nw,ng)
241	!WBARV(L,nw,ng) = WBARV(L-1,nw,ng) !TB17 !! bug
242
243	END DO ! NG gauss point loop
244	END DO ! NW spectral loop
245	print*, 'TBbug WBARV =', WBARV
246
247
248
249	! Total extinction optical depths
250
251	DO NW=1,L_NSPECTV
252	DO NG=1,L_NGAUSS ! full gauss loop
253	TAUV(1,NW,NG)=0.0D0
254	DO L=1,L_NLAYRAD
255	TAUV(L+1,NW,NG)=TAUV(L,NW,NG)+DTAUV(L,NW,NG)
256	END DO
257
258	TAUCUMV(1,NW,NG)=0.0D0
259	DO K=2,L_LEVELS
260	TAUCUMV(K,NW,NG)=TAUCUMV(K-1,NW,NG)+DTAUKV(K,NW,NG)
261	END DO
262	END DO ! end full gauss loop
263	END DO
264
265
266	RETURN
267	END

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

Download in other formats:

Original Format