Context Navigation

optcv.F90 @ 2987

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

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

Download in other formats:

Original Format