Context Navigation

optcv.F90 @ 2001

Last change on this file since 2001 was 1987, checked in by jleconte, 7 years ago

28/08/2018 == JL

Start a series of commits to change the upper boundary conditions in the radiative transfer to solve some issues with the last two layers.
It seems to be good to have aerosols in the first "radiative layer" of the gcm in the IR but visible does not handle very well diffusion in first layer.
Tauaero and tauray are set to 0 (a small value for rayleigh because the code crashes otherwise) in the 4 first semilayers in optcv, but not optci.
This solves random variations of the sw heating at the model top.

Property svn:executable set to *

File size: 12.8 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,glat_ig
7	use gases_h
8	use comcstfi_mod, only: g, r, mugaz
9	use callkeys_mod, only: kastprof,continuum,graybody,H2Ocont_simple,callgasvis
10
11	implicit none
12
13	!==================================================================
14	!
15	! Purpose
16	! -------
17	! Calculates shortwave optical constants at each level.
18	!
19	! Authors
20	! -------
21	! Adapted from the NASA Ames code by R. Wordsworth (2009)
22	!
23	!==================================================================
24	!
25	! THIS SUBROUTINE SETS THE OPTICAL CONSTANTS IN THE VISUAL
26	! IT CALCULATES FOR EACH LAYER, FOR EACH SPECTRAL INTERVAL IN THE VISUAL
27	! LAYER: WBAR, DTAU, COSBAR
28	! LEVEL: TAU
29	!
30	! TAUV(L,NW,NG) is the cumulative optical depth at the top of radiation code
31	! layer L. NW is spectral wavelength interval, ng the Gauss point index.
32	!
33	! TLEV(L) - Temperature at the layer boundary
34	! PLEV(L) - Pressure at the layer boundary (i.e. level)
35	! GASV(NT,NPS,NW,NG) - Visible k-coefficients
36	!
37	!-------------------------------------------------------------------
38
39
40	real*8 DTAUV(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
41	real*8 DTAUKV(L_LEVELS,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
49	! for aerosols
50	real*8 QXVAER(L_LEVELS,L_NSPECTV,NAERKIND)
51	real*8 QSVAER(L_LEVELS,L_NSPECTV,NAERKIND)
52	real*8 GVAER(L_LEVELS,L_NSPECTV,NAERKIND)
53	real*8 TAUAERO(L_LEVELS,NAERKIND)
54	real*8 TAUAEROLK(L_LEVELS,L_NSPECTV,NAERKIND)
55	real*8 TAEROS(L_LEVELS,L_NSPECTV,NAERKIND)
56
57	integer L, NW, NG, K, LK, IAER
58	integer MT(L_LEVELS), MP(L_LEVELS), NP(L_LEVELS)
59	real*8 ANS, TAUGAS
60	real*8 TAURAY(L_NSPECTV)
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)
66	real*8 DCONT,DAERO
67	real*8 DRAYAER
68	double precision wn_cont, p_cont, p_air, T_cont, dtemp, dtempc
69	double precision p_cross
70
71	! variable species mixing ratio variables
72	real*8 QVAR(L_LEVELS), WRATIO(L_LEVELS), MUVAR(L_LEVELS)
73	real*8 KCOEF(4)
74	integer NVAR(L_LEVELS)
75
76	! temporary variables to reduce memory access time to gasv
77	real*8 tmpk(2,2)
78	real*8 tmpkvar(2,2,2)
79
80	! temporary variables for multiple aerosol calculation
81	real*8 atemp(L_NLAYRAD,L_NSPECTV)
82	real*8 btemp(L_NLAYRAD,L_NSPECTV)
83	real*8 ctemp(L_NLAYRAD,L_NSPECTV)
84
85	! variables for k in units m^-1
86	real*8 dz(L_LEVELS)
87
88	integer igas, jgas
89
90	integer interm
91
92	!! AS: to save time in computing continuum (see bilinearbig)
93	IF (.not.ALLOCATED(indv)) THEN
94	ALLOCATE(indv(L_NSPECTV,ngasmx,ngasmx))
95	indv = -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	! Calculate the continuum opacities, i.e., those that do not depend on
102	! NG, the Gauss index.
103
104	taugsurf(:,:) = 0.0
105	dpr(:) = 0.0
106	lkcoef(:,:) = 0.0
107
108	do K=2,L_LEVELS
109	DPR(k) = PLEV(K)-PLEV(K-1)
110
111	! if we have continuum opacities, we need dz
112	if(kastprof)then
113	dz(k) = dpr(k)(1000.0d08.3145d0/muvar(k))TMID(K)/(gPMID(K))
114	U(k) = CmkDPR(k)mugaz/muvar(k)
115	else
116	dz(k) = dpr(k)RTMID(K)/(glat_igPMID(K))mugaz/muvar(k)
117	U(k) = CmkDPR(k)mugaz/muvar(k) ! only Cmk line in optci.F
118	!JL13 the mugaz/muvar factor takes into account water meanmolecular weight if water is present
119	endif
120
121	call tpindex(PMID(K),TMID(K),QVAR(K),pfgasref,tgasref,WREFVAR, &
122	LCOEF,MT(K),MP(K),NVAR(K),WRATIO(K))
123
124	do LK=1,4
125	LKCOEF(K,LK) = LCOEF(LK)
126	end do
127	end do ! levels
128
129	! Spectral dependance of aerosol absorption
130	!JL18 It seems to be good to have aerosols in the first "radiative layer" of the gcm in the IR
131	! but visible does not handle very well diffusion in first layer.
132	! The tauaero and tauray are thus set to 0 (a small value for rayleigh because the code crashes otherwise)
133	! in the 4 first semilayers in optcv, but not optci.
134	! This solves random variations of the sw heating at the model top.
135	do iaer=1,naerkind
136	do NW=1,L_NSPECTV
137	TAEROS(1:4,NW,IAER)=0.d0
138	do K=5,L_LEVELS
139	TAEROS(K,NW,IAER) = TAUAERO(K,IAER) * QXVAER(K,NW,IAER)
140	end do ! levels
141	end do
142	end do
143
144	! Rayleigh scattering
145	do NW=1,L_NSPECTV
146	TRAY(1:4,NW) = 1d-30
147	do K=5,L_LEVELS
148	TRAY(K,NW) = TAURAY(NW) * DPR(K)
149	end do ! levels
150	end do
151
152	! we ignore K=1...
153	do K=2,L_LEVELS
154
155	do NW=1,L_NSPECTV
156
157	DRAYAER = TRAY(K,NW)
158	! DRAYAER is Tau RAYleigh scattering, plus AERosol opacity
159	do iaer=1,naerkind
160	DRAYAER = DRAYAER + TAEROS(K,NW,IAER)
161	end do
162
163	DCONT = 0.0 ! continuum absorption
164
165	if(continuum.and.(.not.graybody).and.callgasvis)then
166	! include continua if necessary
167	wn_cont = dble(wnov(nw))
168	T_cont = dble(TMID(k))
169	do igas=1,ngasmx
170
171	if(gfrac(igas).eq.-1)then ! variable
172	p_cont = dble(PMID(k)scalepQVAR(k)) ! qvar = mol/mol
173	else
174	p_cont = dble(PMID(k)scalepgfrac(igas)*(1.-QVAR(k)))
175	endif
176
177	dtemp=0.0
178	if(igas.eq.igas_N2)then
179
180	interm = indv(nw,igas,igas)
181	! call interpolateN2N2(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
182	indv(nw,igas,igas) = interm
183	! only goes to 500 cm^-1, so unless we're around a cold brown dwarf, this is irrelevant in the visible
184
185	elseif(igas.eq.igas_H2)then
186
187	! first do self-induced absorption
188	interm = indv(nw,igas,igas)
189	call interpolateH2H2(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
190	indv(nw,igas,igas) = interm
191
192	! then cross-interactions with other gases
193	do jgas=1,ngasmx
194	p_cross = dble(PMID(k)scalepgfrac(jgas)*(1.-QVAR(k)))
195	dtempc = 0.0
196	if(jgas.eq.igas_N2)then
197	interm = indv(nw,igas,jgas)
198	call interpolateN2H2(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
199	indv(nw,igas,jgas) = interm
200	! should be irrelevant in the visible
201	elseif(jgas.eq.igas_He)then
202	interm = indv(nw,igas,jgas)
203	call interpolateH2He(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
204	indv(nw,igas,jgas) = interm
205	endif
206	dtemp = dtemp + dtempc
207	enddo
208
209	elseif(igas.eq.igas_H2O.and.T_cont.gt.200.0)then
210
211	p_air = dble(PMID(k)*scalep) - p_cont ! note assumes background is air!
212	if(H2Ocont_simple)then
213	call interpolateH2Ocont_PPC(wn_cont,T_cont,p_cont,p_air,dtemp,.false.)
214	else
215	interm = indv(nw,igas,igas)
216	call interpolateH2Ocont_CKD(wn_cont,T_cont,p_cont,p_air,dtemp,.false.,interm)
217	indv(nw,igas,igas) = interm
218	endif
219
220	endif
221
222	DCONT = DCONT + dtemp
223
224	enddo
225
226	DCONT = DCONT*dz(k)
227
228	endif
229
230	do ng=1,L_NGAUSS-1
231
232	! Now compute TAUGAS
233
234	! Interpolate between water mixing ratios
235	! WRATIO = 0.0 if the requested water amount is equal to, or outside the
236	! the water data range
237
238	if(L_REFVAR.eq.1)then ! added by RW for special no variable case
239
240	! JVO 2017 : added tmpk because the repeated calls to gasi/v increased dramatically
241	! the execution time of optci/v -> ~ factor 2 on the whole radiative
242	! transfer on the tested simulations !
243
244	tmpk = GASV(MT(K):MT(K)+1,MP(K):MP(K)+1,1,NW,NG)
245
246	KCOEF(1) = tmpk(1,1) ! KCOEF(1) = GASV(MT(K),MP(K),1,NW,NG)
247	KCOEF(2) = tmpk(1,2) ! KCOEF(2) = GASV(MT(K),MP(K)+1,1,NW,NG)
248	KCOEF(3) = tmpk(2,2) ! KCOEF(3) = GASV(MT(K)+1,MP(K)+1,1,NW,NG)
249	KCOEF(4) = tmpk(2,1) ! KCOEF(4) = GASV(MT(K)+1,MP(K),1,NW,NG)
250
251	else
252
253	tmpkvar = GASV(MT(K):MT(K)+1,MP(K):MP(K)+1,NVAR(K):NVAR(K)+1,NW,NG)
254
255	KCOEF(1) = tmpkvar(1,1,1) + WRATIO(K) * &
256	( tmpkvar(1,1,2)-tmpkvar(1,1,1) )
257
258	KCOEF(2) = tmpkvar(1,2,1) + WRATIO(K) * &
259	( tmpkvar(1,2,2)-tmpkvar(1,2,1) )
260
261	KCOEF(3) = tmpkvar(2,2,1) + WRATIO(K) * &
262	( tmpkvar(2,2,2)-tmpkvar(2,2,1) )
263
264	KCOEF(4) = tmpkvar(2,1,1) + WRATIO(K) * &
265	( tmpkvar(2,1,2)-tmpkvar(2,1,1) )
266
267
268	endif
269
270	! Interpolate the gaseous k-coefficients to the requested T,P values
271
272	ANS = LKCOEF(K,1)KCOEF(1) + LKCOEF(K,2)KCOEF(2) + &
273	LKCOEF(K,3)KCOEF(3) + LKCOEF(K,4)KCOEF(4)
274
275	TAUGAS = U(k)*ANS
276
277	TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS + DCONT
278	DTAUKV(K,nw,ng) = TAUGAS &
279	+ DRAYAER & ! DRAYAER includes all scattering contributions
280	+ DCONT ! For parameterized continuum aborption
281
282	end do
283
284	! Now fill in the "clear" part of the spectrum (NG = L_NGAUSS),
285	! which holds continuum opacity only
286
287	NG = L_NGAUSS
288	DTAUKV(K,nw,ng) = DRAYAER + DCONT ! Scattering + parameterized continuum absorption
289
290	end do
291	end do
292
293
294	!=======================================================================
295	! Now the full treatment for the layers, where besides the opacity
296	! we need to calculate the scattering albedo and asymmetry factors
297
298	!JL18 It seems to be good to have aerosols in the first "radiative layer" of the gcm in the IR
299	! but not in the visible
300	! The tauaero is thus set to 0 in the 4 first semilayers in optcv, but not optci.
301	! This solves random variations of the sw heating at the model top.
302	do iaer=1,naerkind
303	DO NW=1,L_NSPECTV
304	TAUAEROLK(1:4,NW,IAER)=0.d0
305	DO K=5,L_LEVELS
306	TAUAEROLK(K,NW,IAER) = TAUAERO(K,IAER) * QSVAER(K,NW,IAER) ! effect of scattering albedo
307	ENDDO
308	ENDDO
309	end do
310
311	DO NW=1,L_NSPECTV
312	DO L=1,L_NLAYRAD-1
313	K = 2*L+1
314	atemp(L,NW) = SUM(GVAER(K,NW,1:naerkind) * TAUAEROLK(K,NW,1:naerkind))+SUM(GVAER(K+1,NW,1:naerkind) * TAUAEROLK(K+1,NW,1:naerkind))
315	btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind)) + SUM(TAUAEROLK(K+1,NW,1:naerkind))
316	ctemp(L,NW) = btemp(L,NW) + 0.9999*(TRAY(K,NW) + TRAY(K+1,NW)) ! JVO 2017 : does this 0.999 is really meaningful ?
317	btemp(L,NW) = btemp(L,NW) + TRAY(K,NW) + TRAY(K+1,NW)
318	COSBV(L,NW,1:L_NGAUSS) = atemp(L,NW)/btemp(L,NW)
319	END DO ! L vertical loop
320
321	! Last level
322	L = L_NLAYRAD
323	K = 2*L+1
324	atemp(L,NW) = SUM(GVAER(K,NW,1:naerkind) * TAUAEROLK(K,NW,1:naerkind))
325	btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind))
326	ctemp(L,NW) = btemp(L,NW) + 0.9999*TRAY(K,NW) ! JVO 2017 : does this 0.999 is really meaningful ?
327	btemp(L,NW) = btemp(L,NW) + TRAY(K,NW)
328	COSBV(L,NW,1:L_NGAUSS) = atemp(L,NW)/btemp(L,NW)
329
330
331	END DO ! NW spectral loop
332
333	DO NG=1,L_NGAUSS
334	DO NW=1,L_NSPECTV
335	DO L=1,L_NLAYRAD-1
336
337	K = 2*L+1
338	DTAUV(L,nw,ng) = DTAUKV(K,NW,NG) + DTAUKV(K+1,NW,NG)
339	WBARV(L,nw,ng) = ctemp(L,NW) / DTAUV(L,nw,ng)
340
341	END DO ! L vertical loop
342
343	! Last level
344
345	L = L_NLAYRAD
346	K = 2*L+1
347	DTAUV(L,nw,ng) = DTAUKV(K,NW,NG)
348
349	WBARV(L,NW,NG) = ctemp(L,NW) / DTAUV(L,NW,NG)
350
351	END DO ! NW spectral loop
352	END DO ! NG Gauss loop
353
354	! Total extinction optical depths
355
356	DO NG=1,L_NGAUSS ! full gauss loop
357	DO NW=1,L_NSPECTV
358	TAUCUMV(1,NW,NG)=0.0D0
359	DO K=2,L_LEVELS
360	TAUCUMV(K,NW,NG)=TAUCUMV(K-1,NW,NG)+DTAUKV(K,NW,NG)
361	END DO
362
363	DO L=1,L_NLAYRAD+1
364	TAUV(L,NW,NG)=TAUCUMV(2*L,NW,NG)
365	END DO
366	END DO
367	END DO ! end full gauss loop
368
369
370	return
371
372
373	end subroutine optcv

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

Download in other formats:

Original Format