Context Navigation

optcv.F90 @ 3749

Last change on this file since 3749 was 3695, checked in by debatzbr, 9 months ago
Calculate the visible opacity everywhere BBT
Property svn:executable set to ``*
File size: 13.2 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_N2
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) ! Rayleigh scattering
73	real*8 TRAY(L_LEVELS,L_NSPECTV) ! Rayleigh scattering
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	dtauv(:,:,:) = 0.0
132	tauv(:,:,:) = 0.0
133	taucumv(:,:,:) = 0.0
134
135	taugsurf(:,:) = 0.0
136	dpr(:) = 0.0 ! pressure difference between levels
137	lkcoef(:,:) = 0.0
138	DTAUKV(:,:,:) = 0.0
139
140	do K=2,L_LEVELS
141	DPR(k) = PLEV(K)-PLEV(K-1)
142
143	! if we have continuum opacities, we need dz
144	if(kastprof)then
145	dz(k) = dpr(k)(1000.0d08.3145d0/muvar(k))TMID(K)/(gPMID(K))
146	U(k) = CmkDPR(k)mugaz/muvar(k)
147	else
148	dz(k) = dpr(k)RTMID(K)/(glat_igPMID(K))mugaz/muvar(k)
149	U(k) = CmkDPR(k)mugaz/muvar(k) ! only Cmk line in optci.F
150	!JL13 the mugaz/muvar factor takes into account water meanmolecular weight if water is present
151	endif
152
153	call tpindex(PMID(K),TMID(K),QVAR(K),pfgasref,tgasref,WREFVAR, &
154	LCOEF,MT(K),MP(K),NVAR(K),WRATIO(K))
155
156	do LK=1,4
157	LKCOEF(K,LK) = LCOEF(LK)
158	end do
159	end do ! levels
160
161	! Spectral dependance of aerosol absorption
162	!JL18 It seems to be good to have aerosols in the first "radiative layer" of the gcm in the IR
163	! but visible does not handle very well diffusion in first layer.
164	! The tauaero and tauray are thus set to 0 (a small value for rayleigh because the code crashes otherwise)
165	! in the 4 first semilayers in optcv, but not optci.
166	! This solves random variations of the sw heating at the model top.
167	do iaer=1,naerkind
168	do NW=1,L_NSPECTV
169	TAEROS(1:4,NW,IAER)=0.d0
170	do K=5,L_LEVELS
171	TAEROS(K,NW,IAER) = TAUAERO(K,IAER) * QXVAER(K,NW,IAER)
172	end do ! levels
173	end do
174	end do
175
176	! Rayleigh scattering
177	do NW=1,L_NSPECTV
178	TRAY(1:4,NW) = 1d-30
179	do K=5,L_LEVELS
180	TRAY(K,NW) = TAURAY(NW) * DPR(K)
181	end do ! levels
182	end do
183
184	! we ignore K=1...
185	do K=2,L_LEVELS
186
187	do NW=1,L_NSPECTV
188
189	DRAYAER = TRAY(K,NW)
190	! DRAYAER is Tau RAYleigh scattering, plus AERosol opacity
191	do iaer=1,naerkind
192	DRAYAER = DRAYAER + TAEROS(K,NW,IAER)
193	end do
194
195	DCONT = 0.0 ! continuum absorption
196
197	if(continuum.and.(.not.graybody).and.callgasvis)then
198	! include continua if necessary
199	wn_cont = dble(wnov(nw))
200	T_cont = dble(TMID(k))
201	do igas=1,ngasmx
202
203	if(gfrac(igas).eq.-1)then ! variable
204	p_cont = dble(PMID(k)scalepQVAR(k)) ! qvar = mol/mol
205	else
206	p_cont = dble(PMID(k)scalepgfrac(igas)*(1.-QVAR(k)))
207	endif
208
209	dtemp=0.0
210	if(igas.eq.igas_N2)then
211
212	interm = indv(nw,igas,igas)
213	! call interpolateN2N2(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
214	indv(nw,igas,igas) = interm
215	! only goes to 500 cm^-1, so unless we're around a cold brown dwarf, this is irrelevant in the visible
216
217	! elseif(igas.eq.igas_H2)then !AF24: removed
218
219	elseif(igas.eq.igas_CH4)then
220
221	! first do self-induced absorption
222	interm = indv(nw,igas,igas)
223	call interpolateCH4CH4(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
224	indv(nw,igas,igas) = interm
225
226	! then cross-interactions with other gases !AF24: removed
227
228	! elseif(igas.eq.igas_H2O.and.T_cont.gt.100.0)then !AF24: removed
229
230	endif
231	DCONT = DCONT + dtemp
232
233	enddo
234
235	DCONT = DCONT*dz(k)
236
237	endif
238
239	do ng=1,L_NGAUSS-1
240
241	! Now compute TAUGAS
242
243	! Interpolate between water mixing ratios
244	! WRATIO = 0.0 if the requested water amount is equal to, or outside the
245	! the water data range
246
247	if(L_REFVAR.eq.1)then ! added by RW for special no variable case
248
249	! JVO 2017 : added tmpk because the repeated calls to gasi/v increased dramatically
250	! the execution time of optci/v -> ~ factor 2 on the whole radiative
251	! transfer on the tested simulations !
252
253	IF (corrk_recombin) THEN ! Added by JVO
254	tmpk = GASV_RECOMB(MT(K):MT(K)+1,MP(K):MP(K)+1,1,NW,NG) ! contains the mix of recombined species
255	ELSE
256	tmpk = GASV(MT(K):MT(K)+1,MP(K):MP(K)+1,1,NW,NG)
257	ENDIF
258
259	KCOEF(1) = tmpk(1,1) ! KCOEF(1) = GASV(MT(K),MP(K),1,NW,NG)
260	KCOEF(2) = tmpk(1,2) ! KCOEF(2) = GASV(MT(K),MP(K)+1,1,NW,NG)
261	KCOEF(3) = tmpk(2,2) ! KCOEF(3) = GASV(MT(K)+1,MP(K)+1,1,NW,NG)
262	KCOEF(4) = tmpk(2,1) ! KCOEF(4) = GASV(MT(K)+1,MP(K),1,NW,NG)
263
264	else
265
266	IF (corrk_recombin) THEN
267	tmpkvar = GASV_RECOMB(MT(K):MT(K)+1,MP(K):MP(K)+1,NVAR(K):NVAR(K)+1,NW,NG)
268	ELSE
269	tmpkvar = GASV(MT(K):MT(K)+1,MP(K):MP(K)+1,NVAR(K):NVAR(K)+1,NW,NG)
270	ENDIF
271
272	KCOEF(1) = tmpkvar(1,1,1) + WRATIO(K) * &
273	( tmpkvar(1,1,2)-tmpkvar(1,1,1) )
274
275	KCOEF(2) = tmpkvar(1,2,1) + WRATIO(K) * &
276	( tmpkvar(1,2,2)-tmpkvar(1,2,1) )
277
278	KCOEF(3) = tmpkvar(2,2,1) + WRATIO(K) * &
279	( tmpkvar(2,2,2)-tmpkvar(2,2,1) )
280
281	KCOEF(4) = tmpkvar(2,1,1) + WRATIO(K) * &
282	( tmpkvar(2,1,2)-tmpkvar(2,1,1) )
283
284
285	endif
286
287	! Interpolate the gaseous k-coefficients to the requested T,P values
288
289	ANS = LKCOEF(K,1)KCOEF(1) + LKCOEF(K,2)KCOEF(2) + &
290	LKCOEF(K,3)KCOEF(3) + LKCOEF(K,4)KCOEF(4)
291
292	TAUGAS = U(k)*ANS
293
294	TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS + DCONT
295	DTAUKV(K,nw,ng) = TAUGAS &
296	+ DRAYAER & ! DRAYAER includes all scattering contributions
297	+ DCONT ! For parameterized continuum aborption
298
299	end do
300
301	! Now fill in the "clear" part of the spectrum (NG = L_NGAUSS),
302	! which holds continuum opacity only
303
304	NG = L_NGAUSS
305	DTAUKV(K,nw,ng) = DRAYAER + DCONT ! Scattering + parameterized continuum absorption
306
307	end do
308	end do
309
310
311	!=======================================================================
312	! Now the full treatment for the layers, where besides the opacity
313	! we need to calculate the scattering albedo and asymmetry factors
314
315	!JL18 It seems to be good to have aerosols in the first "radiative layer" of the gcm in the IR
316	! but not in the visible
317	! The tauaero is thus set to 0 in the 4 first semilayers in optcv, but not optci.
318	! This solves random variations of the sw heating at the model top.
319	do iaer=1,naerkind
320	DO NW=1,L_NSPECTV
321	TAUAEROLK(1:4,NW,IAER)=0.d0
322	DO K=5,L_LEVELS
323	TAUAEROLK(K,NW,IAER) = TAUAERO(K,IAER) * QSVAER(K,NW,IAER) ! effect of scattering albedo
324	ENDDO
325	ENDDO
326	end do
327
328	DO NW=1,L_NSPECTV
329	DO L=1,L_NLAYRAD-1
330	K = 2*L+1
331	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))
332	btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind)) + SUM(TAUAEROLK(K+1,NW,1:naerkind))
333	ctemp(L,NW) = btemp(L,NW) + 0.9999*(TRAY(K,NW) + TRAY(K+1,NW)) ! JVO 2017 : does this 0.999 is really meaningful ?
334	btemp(L,NW) = btemp(L,NW) + TRAY(K,NW) + TRAY(K+1,NW)
335	COSBV(L,NW,1:L_NGAUSS) = atemp(L,NW)/btemp(L,NW)
336	END DO ! L vertical loop
337
338	! Last level
339	L = L_NLAYRAD
340	K = 2*L+1
341	atemp(L,NW) = SUM(GVAER(K,NW,1:naerkind) * TAUAEROLK(K,NW,1:naerkind))
342	btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind))
343	ctemp(L,NW) = btemp(L,NW) + 0.9999*TRAY(K,NW) ! JVO 2017 : does this 0.999 is really meaningful ?
344	btemp(L,NW) = btemp(L,NW) + TRAY(K,NW)
345	COSBV(L,NW,1:L_NGAUSS) = atemp(L,NW)/btemp(L,NW)
346
347
348	END DO ! NW spectral loop
349
350	DO NG=1,L_NGAUSS
351	DO NW=1,L_NSPECTV
352	DO L=1,L_NLAYRAD-1
353
354	K = 2*L+1
355	DTAUV(L,nw,ng) = DTAUKV(K,NW,NG) + DTAUKV(K+1,NW,NG)
356	WBARV(L,nw,ng) = ctemp(L,NW) / DTAUV(L,nw,ng)
357
358	END DO ! L vertical loop
359
360	! Last level
361
362	L = L_NLAYRAD
363	K = 2*L+1
364	DTAUV(L,nw,ng) = DTAUKV(K,NW,NG)
365
366	WBARV(L,NW,NG) = ctemp(L,NW) / DTAUV(L,NW,NG)
367
368	END DO ! NW spectral loop
369	END DO ! NG Gauss loop
370
371	! Total extinction optical depths
372	DO NG=1,L_NGAUSS ! full gauss loop
373	DO NW=1,L_NSPECTV
374	TAUV(1,NW,NG)=0.0D0
375	TAUCUMV(1,NW,NG)=0.0D0
376
377	DO K=2,L_LEVELS
378	TAUCUMV(K,NW,NG)=TAUCUMV(K-1,NW,NG)+DTAUKV(K,NW,NG)
379	END DO
380
381	DO L=1,L_NLAYRAD
382	TAUV(L,NW,NG)=TAUCUMV(2*L,NW,NG)
383	END DO
384	TAUV(L,NW,NG)=TAUCUMV(2*L_NLAYRAD+1,NW,NG)
385	END DO
386	END DO ! end full gauss loop
387
388	end subroutine optcv
389
390	END MODULE optcv_mod

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

Download in other formats:

Original Format