Context Navigation

optcv.F90 @ 2921

Last change on this file since 2921 was 2875, checked in by emillour, 2 years ago

Generic PCM:
Some minor fixes:

missing igas_CO2 in "use gases_h" in optci.F90, optcv.F90, sugas_corrk.F90
handle case when moist adjustement is not called in physiq: some of its outputs still need be set to zero as they are used later on.

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

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

Download in other formats:

Original Format