Context Navigation

optcv.F90 @ 2142

Last change on this file since 2142 was 2133, checked in by jvatant, 6 years ago
Fic a bug in r2131 the writediagfi cannot be called in RT routines because of the iradia it goes crazy with outputs freq ... --JVO
Property svn:executable set to ``*
File size: 13.0 KB

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

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

Download in other formats:

Original Format