Context Navigation

optcv.F90 @ 2865

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

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

Download in other formats:

Original Format