Context Navigation

optcv.F90 @ 3300

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

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

Download in other formats:

Original Format