Context Navigation

optcv.F90 @ 2590

Last change on this file since 2590 was 2582, checked in by emillour, 4 years ago
Generic GCM: Fixed bug in tpindex (for low temperatures, between first and second reference temperatures, temperature was wrongly set to tref(1)). The input temperature was also allowed to be modified by the routine, which is probably not a good idea and no longer the case. Took this opportunity to turn tpindex into a module. GM+EM
Property svn:executable set to ``*
File size: 13.4 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,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_He)then
211	interm = indv(nw,igas,jgas)
212	call interpolateH2He(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
213	indv(nw,igas,jgas) = interm
214	endif
215	dtemp = dtemp + dtempc
216	enddo
217
218	elseif(igas.eq.igas_H2O.and.T_cont.gt.100.0)then
219	! Compute self and foreign (with air) continuum of H2O
220	p_air = dble(PMID(k)*scalep) - p_cont ! note assumes background is air!
221	interm = indv(nw,igas,igas)
222	call interpolateH2O_self_foreign(wn_cont,T_cont,p_cont,p_air,dtemp,.false.,interm) ! MTCKD v3.3
223	indv(nw,igas,igas) = interm
224
225	endif
226
227	DCONT = DCONT + dtemp
228
229	enddo
230
231	DCONT = DCONT*dz(k)
232
233	endif
234
235	do ng=1,L_NGAUSS-1
236
237	! Now compute TAUGAS
238
239	! Interpolate between water mixing ratios
240	! WRATIO = 0.0 if the requested water amount is equal to, or outside the
241	! the water data range
242
243	if(L_REFVAR.eq.1)then ! added by RW for special no variable case
244
245	! JVO 2017 : added tmpk because the repeated calls to gasi/v increased dramatically
246	! the execution time of optci/v -> ~ factor 2 on the whole radiative
247	! transfer on the tested simulations !
248
249	IF (corrk_recombin) THEN ! Added by JVO
250	tmpk = GASV_RECOMB(MT(K):MT(K)+1,MP(K):MP(K)+1,1,NW,NG) ! contains the mix of recombined species
251	ELSE
252	tmpk = GASV(MT(K):MT(K)+1,MP(K):MP(K)+1,1,NW,NG)
253	ENDIF
254
255	KCOEF(1) = tmpk(1,1) ! KCOEF(1) = GASV(MT(K),MP(K),1,NW,NG)
256	KCOEF(2) = tmpk(1,2) ! KCOEF(2) = GASV(MT(K),MP(K)+1,1,NW,NG)
257	KCOEF(3) = tmpk(2,2) ! KCOEF(3) = GASV(MT(K)+1,MP(K)+1,1,NW,NG)
258	KCOEF(4) = tmpk(2,1) ! KCOEF(4) = GASV(MT(K)+1,MP(K),1,NW,NG)
259
260	else
261
262	IF (corrk_recombin) THEN
263	tmpkvar = GASV_RECOMB(MT(K):MT(K)+1,MP(K):MP(K)+1,NVAR(K):NVAR(K)+1,NW,NG)
264	ELSE
265	tmpkvar = GASV(MT(K):MT(K)+1,MP(K):MP(K)+1,NVAR(K):NVAR(K)+1,NW,NG)
266	ENDIF
267
268	KCOEF(1) = tmpkvar(1,1,1) + WRATIO(K) * &
269	( tmpkvar(1,1,2)-tmpkvar(1,1,1) )
270
271	KCOEF(2) = tmpkvar(1,2,1) + WRATIO(K) * &
272	( tmpkvar(1,2,2)-tmpkvar(1,2,1) )
273
274	KCOEF(3) = tmpkvar(2,2,1) + WRATIO(K) * &
275	( tmpkvar(2,2,2)-tmpkvar(2,2,1) )
276
277	KCOEF(4) = tmpkvar(2,1,1) + WRATIO(K) * &
278	( tmpkvar(2,1,2)-tmpkvar(2,1,1) )
279
280
281	endif
282
283	! Interpolate the gaseous k-coefficients to the requested T,P values
284
285	ANS = LKCOEF(K,1)KCOEF(1) + LKCOEF(K,2)KCOEF(2) + &
286	LKCOEF(K,3)KCOEF(3) + LKCOEF(K,4)KCOEF(4)
287
288	TAUGAS = U(k)*ANS
289
290	TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS + DCONT
291	DTAUKV(K,nw,ng) = TAUGAS &
292	+ DRAYAER & ! DRAYAER includes all scattering contributions
293	+ DCONT ! For parameterized continuum aborption
294
295	end do
296
297	! Now fill in the "clear" part of the spectrum (NG = L_NGAUSS),
298	! which holds continuum opacity only
299
300	NG = L_NGAUSS
301	DTAUKV(K,nw,ng) = DRAYAER + DCONT ! Scattering + parameterized continuum absorption
302
303	end do
304	end do
305
306
307	!=======================================================================
308	! Now the full treatment for the layers, where besides the opacity
309	! we need to calculate the scattering albedo and asymmetry factors
310
311	!JL18 It seems to be good to have aerosols in the first "radiative layer" of the gcm in the IR
312	! but not in the visible
313	! The tauaero is thus set to 0 in the 4 first semilayers in optcv, but not optci.
314	! This solves random variations of the sw heating at the model top.
315	do iaer=1,naerkind
316	DO NW=1,L_NSPECTV
317	TAUAEROLK(1:4,NW,IAER)=0.d0
318	DO K=5,L_LEVELS
319	TAUAEROLK(K,NW,IAER) = TAUAERO(K,IAER) * QSVAER(K,NW,IAER) ! effect of scattering albedo
320	ENDDO
321	ENDDO
322	end do
323
324	DO NW=1,L_NSPECTV
325	DO L=1,L_NLAYRAD-1
326	K = 2*L+1
327	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))
328	btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind)) + SUM(TAUAEROLK(K+1,NW,1:naerkind))
329	ctemp(L,NW) = btemp(L,NW) + 0.9999*(TRAY(K,NW) + TRAY(K+1,NW)) ! JVO 2017 : does this 0.999 is really meaningful ?
330	btemp(L,NW) = btemp(L,NW) + TRAY(K,NW) + TRAY(K+1,NW)
331	COSBV(L,NW,1:L_NGAUSS) = atemp(L,NW)/btemp(L,NW)
332	END DO ! L vertical loop
333
334	! Last level
335	L = L_NLAYRAD
336	K = 2*L+1
337	atemp(L,NW) = SUM(GVAER(K,NW,1:naerkind) * TAUAEROLK(K,NW,1:naerkind))
338	btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind))
339	ctemp(L,NW) = btemp(L,NW) + 0.9999*TRAY(K,NW) ! JVO 2017 : does this 0.999 is really meaningful ?
340	btemp(L,NW) = btemp(L,NW) + TRAY(K,NW)
341	COSBV(L,NW,1:L_NGAUSS) = atemp(L,NW)/btemp(L,NW)
342
343
344	END DO ! NW spectral loop
345
346	DO NG=1,L_NGAUSS
347	DO NW=1,L_NSPECTV
348	DO L=1,L_NLAYRAD-1
349
350	K = 2*L+1
351	DTAUV(L,nw,ng) = DTAUKV(K,NW,NG) + DTAUKV(K+1,NW,NG)
352	WBARV(L,nw,ng) = ctemp(L,NW) / DTAUV(L,nw,ng)
353
354	END DO ! L vertical loop
355
356	! Last level
357
358	L = L_NLAYRAD
359	K = 2*L+1
360	DTAUV(L,nw,ng) = DTAUKV(K,NW,NG)
361
362	WBARV(L,NW,NG) = ctemp(L,NW) / DTAUV(L,NW,NG)
363
364	END DO ! NW spectral loop
365	END DO ! NG Gauss loop
366
367	! Total extinction optical depths
368
369	DO NG=1,L_NGAUSS ! full gauss loop
370	DO NW=1,L_NSPECTV
371	TAUCUMV(1,NW,NG)=0.0D0
372	DO K=2,L_LEVELS
373	TAUCUMV(K,NW,NG)=TAUCUMV(K-1,NW,NG)+DTAUKV(K,NW,NG)
374	END DO
375
376	DO L=1,L_NLAYRAD
377	TAUV(L,NW,NG)=TAUCUMV(2*L,NW,NG)
378	END DO
379	TAUV(L,NW,NG)=TAUCUMV(2*L_NLAYRAD+1,NW,NG)
380	END DO
381	END DO ! end full gauss loop
382
383
384
385
386	end subroutine optcv
387
388	END MODULE optcv_mod

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

Download in other formats:

Original Format