Context Navigation

optcv.F90 @ 2245

Last change on this file since 2245 was 2242, checked in by jvatant, 5 years ago

+ Update the interface with YAMMS so it now correctly handles the small values of the moments,
requiring dynamics to have a threshold quite low (set to 1D-200) and a sanity check in calmufi
corresponding to this value. Thus it removes 'most' of the unphysical radius obtained in
YAMMS. There are still some, but at least there is no more problem of model stability for the moments
and the code can run.
Still, take care the day you want to calculate opacities from the radii and not the moments.
Although, note that there are some negative values, in the output files, but theses negatives are

harmless, as they are only present in output files, dynamics reseting to epsilon after.

+ Given theses pbs with the radii, also update the optics so that it computes the opacities in
a simpler way, directly for M3, through look-up tables, M3 being a good proxy.
--JVO

Property svn:executable set to *

File size: 14.7 KB

Line
1	SUBROUTINE OPTCV(PQMO,NLAY,PLEV,TMID,PMID, &
2	DTAUV,TAUV,TAUCUMV,WBARV,COSBV,TAURAY,TAUGSURF,SEASHAZEFACT)
3
4	use radinc_h
5	use radcommon_h, only: gasv,gasv_recomb,tlimit,Cmk,gzlat_ig, &
6	tgasref,pfgasref,wnov,scalep,indv
7	use gases_h
8	use datafile_mod, only: haze_opt_file
9	use comcstfi_mod, only: r
10	use callkeys_mod, only: continuum,graybody,callgasvis,corrk_recombin, &
11	callclouds,callmufi,seashaze,uncoupl_optic_haze
12	use tracer_h, only: nmicro,nice
13	use MMP_OPTICS
14
15	implicit none
16
17	!==================================================================
18	!
19	! Purpose
20	! -------
21	! Calculates shortwave optical constants at each level.
22	!
23	! Authors
24	! -------
25	! Adapted from the NASA Ames code by R. Wordsworth (2009)
26	! Clean and adaptation to Titan by J. Vatant d'Ollone (2016-17)
27	!
28	!==================================================================
29	!
30	! THIS SUBROUTINE SETS THE OPTICAL CONSTANTS IN THE VISUAL
31	! IT CALCULATES FOR EACH LAYER, FOR EACH SPECTRAL INTERVAL IN THE VISUAL
32	! LAYER: WBAR, DTAU, COSBAR
33	! LEVEL: TAU
34	!
35	! TAUV(L,NW,NG) is the cumulative optical depth at the top of radiation code
36	! layer L. NW is spectral wavelength interval, ng the Gauss point index.
37	!
38	! TLEV(L) - Temperature at the layer boundary
39	! PLEV(L) - Pressure at the layer boundary (i.e. level)
40	! GASV(NT,NPS,NW,NG) - Visible k-coefficients
41	!
42	!-------------------------------------------------------------------
43
44
45	!==========================================================
46	! Input/Output
47	!==========================================================
48	REAL*8, INTENT(IN) :: PQMO(nlay,nmicro) ! Tracers for microphysics optics (X/m2).
49	INTEGER, INTENT(IN) :: NLAY ! Number of pressure layers (for pqmo)
50	REAL*8, INTENT(IN) :: PLEV(L_LEVELS)
51	REAL*8, INTENT(IN) :: TMID(L_LEVELS), PMID(L_LEVELS)
52	REAL*8, INTENT(IN) :: TAURAY(L_NSPECTV)
53	REAL*8, INTENT(IN) :: SEASHAZEFACT(L_LEVELS)
54
55	REAL*8, INTENT(OUT) :: DTAUV(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
56	REAL*8, INTENT(OUT) :: TAUV(L_NLEVRAD,L_NSPECTV,L_NGAUSS)
57	REAL*8, INTENT(OUT) :: TAUCUMV(L_LEVELS,L_NSPECTV,L_NGAUSS)
58	REAL*8, INTENT(OUT) :: COSBV(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
59	REAL*8, INTENT(OUT) :: WBARV(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
60	REAL*8, INTENT(OUT) :: TAUGSURF(L_NSPECTV,L_NGAUSS-1)
61	! ==========================================================
62
63	real*8 DTAUKV(L_LEVELS,L_NSPECTV,L_NGAUSS)
64
65	! Titan customisation
66	! J. Vatant d'Ollone (2016)
67	real*8 DHAZE_T(L_LEVELS,L_NSPECTI)
68	real*8 DHAZES_T(L_LEVELS,L_NSPECTI)
69	real*8 SSA_T(L_LEVELS,L_NSPECTI)
70	real*8 ASF_T(L_LEVELS,L_NSPECTI)
71	! ==========================
72
73	integer L, NW, NG, K, LK, IAER
74	integer MT(L_LEVELS), MP(L_LEVELS), NP(L_LEVELS)
75	real*8 ANS, TAUGAS
76	real*8 TRAY(L_LEVELS,L_NSPECTV)
77	real*8 DPR(L_LEVELS), U(L_LEVELS)
78	real*8 LCOEF(4), LKCOEF(L_LEVELS,4)
79
80	real*8 DCONT
81	real*8 DRAYAER
82	double precision wn_cont, p_cont, p_air, T_cont, dtemp, dtempc
83	double precision p_cross
84
85	real*8 KCOEF(4)
86
87	! temporary variable to reduce memory access time to gasv
88	real*8 tmpk(2,2)
89
90	! temporary variables for multiple aerosol calculation
91	real*8 atemp(L_NLAYRAD,L_NSPECTV)
92	real*8 btemp(L_NLAYRAD,L_NSPECTV)
93	real*8 ctemp(L_NLAYRAD,L_NSPECTV)
94
95	! variables for k in units m^-1
96	real*8 dz(L_LEVELS)
97
98	integer igas, jgas, ilay
99
100	integer interm
101
102	! Variables for haze optics
103	character(len=200) file_path
104	logical file_ok
105	integer dumch
106	real*8 dumwvl
107
108	real*8 m3as,m3af
109	real*8 dtauaer_s,dtauaer_f
110	real*8,save :: rhoaer_s(L_NSPECTV),ssa_s(L_NSPECTV),asf_s(L_NSPECTV)
111	real*8,save :: rhoaer_f(L_NSPECTV),ssa_f(L_NSPECTV),asf_f(L_NSPECTV)
112	!$OMP THREADPRIVATE(rhoaer_s,rhoaer_f,ssa_s,ssa_f,asf_s,asf_f)
113
114	logical,save :: firstcall=.true.
115	!$OMP THREADPRIVATE(firstcall)
116
117
118	!! AS: to save time in computing continuum (see bilinearbig)
119	IF (.not.ALLOCATED(indv)) THEN
120	ALLOCATE(indv(L_NSPECTV,ngasmx,ngasmx))
121	indv = -9999 ! this initial value means "to be calculated"
122	ENDIF
123
124	! Some initialisation beacause there's a pb with disr_haze at the limits (nw=1)
125	! I should check this - For now we set vars to zero : better than nans - JVO 2017
126	DHAZE_T(:,:) = 0.0
127	SSA_T(:,:) = 0.0
128	ASF_T(:,:) = 0.0
129
130	! Load tabulated haze optical properties if needed.
131	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
132	IF (firstcall .AND. callmufi .AND. (.NOT. uncoupl_optic_haze)) THEN
133	OPEN(12,file=TRIM(haze_opt_file),form='formatted') ! The file has been inquired in physiq_mod firstcall
134	READ(12,*) ! dummy header
135	DO NW=1,L_NSPECTI
136	READ(12,*) ! there's IR 1st
137	ENDDO
138	DO NW=1,L_NSPECTV
139	READ(12,*) dumch, dumwvl, rhoaer_f(nw), ssa_f(nw), asf_f(nw), rhoaer_s(nw), ssa_s(nw), asf_s(nw)
140	ENDDO
141	CLOSE(12)
142	ENDIF
143
144	!=======================================================================
145	! Determine the total gas opacity throughout the column, for each
146	! spectral interval, NW, and each Gauss point, NG.
147	! Calculate the continuum opacities, i.e., those that do not depend on
148	! NG, the Gauss index.
149
150	taugsurf(:,:) = 0.0
151	dpr(:) = 0.0
152	lkcoef(:,:) = 0.0
153
154	do K=2,L_LEVELS
155
156	ilay = L_NLAYRAD+1 - k/2 ! int. arithmetic => gives the gcm layer index (reversed)
157
158	DPR(k) = PLEV(K)-PLEV(K-1)
159
160	! if we have continuum opacities, we need dz
161
162	dz(k) = dpr(k)RTMID(K)/(gzlat_ig(ilay)*PMID(K))
163	U(k) = Cmk(ilay)*DPR(k) ! only Cmk line in optcv.F
164
165	call tpindex(PMID(K),TMID(K),pfgasref,tgasref,LCOEF,MT(K),MP(K))
166
167	do LK=1,4
168	LKCOEF(K,LK) = LCOEF(LK)
169	end do
170	end do ! levels
171
172	! Rayleigh scattering
173	do NW=1,L_NSPECTV
174	TRAY(1:4,NW) = 1.d-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	ilay = L_NLAYRAD+1 - k/2 ! int. arithmetic => gives the gcm layer index (reversed)
184
185	do NW=1,L_NSPECTV
186
187	IF (callmufi .AND. (.NOT. uncoupl_optic_haze)) THEN
188	m3as = pqmo(ilay,2) / 2.0
189	m3af = pqmo(ilay,4) / 2.0
190
191	IF ( ilay .lt. 18 ) THEN
192	m3as = pqmo(18,2) / 2.0 * dz(k) / dz(18)
193	m3af = pqmo(18,4) / 2.0 * dz(k) / dz(18)
194	ENDIF
195
196	dtauaer_s = m3as*rhoaer_s(nw)
197	dtauaer_f = m3af*rhoaer_f(nw)
198	DHAZE_T(k,nw) = dtauaer_s + dtauaer_f
199
200	IF ( dtauaer_s + dtauaer_f .GT. 1.D-30 ) THEN
201	SSA_T(k,nw) = ( dtauaer_sssa_s(nw) + dtauaer_fssa_f(nw) ) / ( dtauaer_s+dtauaer_f )
202	ASF_T(k,nw) = ( dtauaer_sssa_s(nw)asf_s(nw) + dtauaer_fssa_f(nw)asf_f(nw) ) &
203	/ ( ssa_s(nw)dtauaer_s + ssa_f(nw)dtauaer_f )
204	ELSE
205	DHAZE_T(k,nw) = 0.D0
206	SSA_T(k,nw) = 1.0
207	ASF_T(k,nw) = 1.0
208	ENDIF
209
210	IF (callclouds.and.firstcall) &
211	WRITE(,) 'WARNING: In optcv, optical properties &
212	&calculations are not implemented yet'
213	ELSE
214	! Call fixed vertical haze profile of extinction - same for all columns
215	call disr_haze(dz(k),plev(k),wnov(nw),DHAZE_T(k,nw),SSA_T(k,nw),ASF_T(k,nw))
216	if (seashaze) DHAZE_T(k,nw) = DHAZE_T(k,nw)*seashazefact(k)
217	ENDIF
218
219	!JL18 It seems to be good to have aerosols in the first "radiative layer" of the gcm in the IR
220	! but visible does not handle very well diffusion in first layer.
221	! The tauaero and tauray are thus set to 0 (a small value for rayleigh because the code crashes otherwise)
222	! in the 4 first semilayers in optcv, but not optci.
223	! This solves random variations of the sw heating at the model top.
224	if (k<5) DHAZE_T(K,:) = 0.0
225
226	DRAYAER = TRAY(K,NW)
227	! DRAYAER is Tau RAYleigh scattering, plus AERosol opacity
228	DRAYAER = DRAYAER + DHAZE_T(K,NW) ! Titan's aerosol
229
230	DCONT = 0.0 ! continuum absorption
231
232	if(continuum.and.(.not.graybody).and.callgasvis)then
233	! include continua if necessary
234	wn_cont = dble(wnov(nw))
235	T_cont = dble(TMID(k))
236	do igas=1,ngasmx
237
238	p_cont = dble(PMID(k)scalepgfrac(igas,ilay))
239
240	dtemp=0.0
241	if(igas.eq.igas_N2)then
242
243	interm = indv(nw,igas,igas)
244	! call interpolateN2N2(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
245	indv(nw,igas,igas) = interm
246	! only goes to 500 cm^-1, so unless we're around a cold brown dwarf, this is irrelevant in the visible
247
248	elseif(igas.eq.igas_H2)then
249
250	! first do self-induced absorption
251	interm = indv(nw,igas,igas)
252	call interpolateH2H2(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	p_cross = dble(PMID(k)scalepgfrac(jgas,ilay))
258	dtempc = 0.0
259	if(jgas.eq.igas_N2)then
260	interm = indv(nw,igas,jgas)
261	call interpolateN2H2(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
262	indv(nw,igas,jgas) = interm
263	! should be irrelevant in the visible
264	endif
265	dtemp = dtemp + dtempc
266	enddo
267
268	elseif(igas.eq.igas_CH4)then
269
270	! first do self-induced absorption
271	interm = indv(nw,igas,igas)
272	call interpolateCH4CH4(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
273	indv(nw,igas,igas) = interm
274
275	! then cross-interactions with other gases
276	do jgas=1,ngasmx
277	p_cross = dble(PMID(k)scalepgfrac(jgas,ilay))
278	dtempc = 0.0
279	if(jgas.eq.igas_N2)then
280	interm = indv(nw,igas,jgas)
281	call interpolateN2CH4(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
282	indv(nw,igas,jgas) = interm
283	endif
284	dtemp = dtemp + dtempc
285	enddo
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	! JVO 2017 : added tmpk because the repeated calls to gasi/v increased dramatically
302	! the execution time of optci/v -> ~ factor 2 on the whole radiative
303	! transfer on the tested simulations !
304
305	if (corrk_recombin) then
306	tmpk = GASV_RECOMB(MT(K):MT(K)+1,MP(K):MP(K)+1,NW,NG)
307	else
308	tmpk = GASV(MT(K):MT(K)+1,MP(K):MP(K)+1,1,NW,NG)
309	endif
310
311	KCOEF(1) = tmpk(1,1) ! KCOEF(1) = GASV(MT(K),MP(K),1,NW,NG)
312	KCOEF(2) = tmpk(1,2) ! KCOEF(2) = GASV(MT(K),MP(K)+1,1,NW,NG)
313	KCOEF(3) = tmpk(2,2) ! KCOEF(3) = GASV(MT(K)+1,MP(K)+1,1,NW,NG)
314	KCOEF(4) = tmpk(2,1) ! KCOEF(4) = GASV(MT(K)+1,MP(K),1,NW,NG)
315
316	! Interpolate the gaseous k-coefficients to the requested T,P values
317
318	ANS = LKCOEF(K,1)KCOEF(1) + LKCOEF(K,2)KCOEF(2) + &
319	LKCOEF(K,3)KCOEF(3) + LKCOEF(K,4)KCOEF(4)
320
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, including Titan's haze
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
346	! Haze scattering
347	!JL18 It seems to be good to have aerosols in the first "radiative layer" of the gcm in the IR
348	! but not in the visible
349	! The dhaze_s is thus set to 0 in the 4 first semilayers in optcv, but not optci.
350	! This solves random variations of the sw heating at the model top.
351	DO NW=1,L_NSPECTV
352	DHAZES_T(1:4,NW) = 0.d0
353	DO K=5,L_LEVELS
354	DHAZES_T(K,NW) = DHAZE_T(K,NW) * SSA_T(K,NW) ! effect of scattering albedo on haze
355	ENDDO
356	ENDDO
357
358
359	DO NW=1,L_NSPECTV
360	DO L=1,L_NLAYRAD-1
361	K = 2*L+1
362	atemp(L,NW) = ASF_T(K,NW)DHAZES_T(K,NW) + ASF_T(K+1,NW)DHAZES_T(K+1,NW)
363	btemp(L,NW) = DHAZES_T(K,NW) + DHAZES_T(K+1,NW)
364	ctemp(L,NW) = btemp(L,NW) + 0.9999*(TRAY(K,NW) + TRAY(K+1,NW)) ! JVO 2017 : does this 0.999 is really meaningful ?
365	btemp(L,NW) = btemp(L,NW) + TRAY(K,NW) + TRAY(K+1,NW)
366	COSBV(L,NW,1:L_NGAUSS) = atemp(L,NW)/btemp(L,NW)
367	END DO ! L vertical loop
368
369	! Last level
370	L = L_NLAYRAD
371	K = 2*L+1
372	atemp(L,NW) = ASF_T(K,NW)*DHAZES_T(K,NW)
373	btemp(L,NW) = DHAZES_T(K,NW)
374	ctemp(L,NW) = btemp(L,NW) + 0.9999*TRAY(K,NW) ! JVO 2017 : does this 0.999 is really meaningful ?
375	btemp(L,NW) = btemp(L,NW) + TRAY(K,NW)
376	COSBV(L,NW,1:L_NGAUSS) = atemp(L,NW)/btemp(L,NW)
377
378
379	END DO ! NW spectral loop
380
381	DO NG=1,L_NGAUSS
382	DO NW=1,L_NSPECTV
383	DO L=1,L_NLAYRAD-1
384
385	K = 2*L+1
386	DTAUV(L,nw,ng) = DTAUKV(K,NW,NG) + DTAUKV(K+1,NW,NG)
387	WBARV(L,nw,ng) = ctemp(L,NW) / DTAUV(L,nw,ng)
388
389	END DO ! L vertical loop
390
391	! Last level
392
393	L = L_NLAYRAD
394	K = 2*L+1
395	DTAUV(L,nw,ng) = DTAUKV(K,NW,NG)
396
397	WBARV(L,NW,NG) = ctemp(L,NW) / DTAUV(L,NW,NG)
398
399	END DO ! NW spectral loop
400	END DO ! NG Gauss loop
401
402	! Total extinction optical depths
403
404	DO NG=1,L_NGAUSS ! full gauss loop
405	DO NW=1,L_NSPECTV
406	TAUCUMV(1,NW,NG)=0.0D0
407	DO K=2,L_LEVELS
408	TAUCUMV(K,NW,NG)=TAUCUMV(K-1,NW,NG)+DTAUKV(K,NW,NG)
409	END DO
410
411	DO L=1,L_NLAYRAD
412	TAUV(L,NW,NG)=TAUCUMV(2*L,NW,NG)
413	END DO
414	TAUV(L,NW,NG)=TAUCUMV(2*L_NLAYRAD+1,NW,NG)
415	END DO
416	END DO ! end full gauss loop
417
418	if(firstcall) firstcall = .false.
419
420	return
421
422
423	end subroutine optcv

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

Download in other formats:

Original Format