Context Navigation

optci.F90 @ 3026

Last change on this file since 3026 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

File size: 14.0 KB

Line
1	subroutine optci(PQMO,NLAY,PLEV,TLEV,TMID,PMID, &
2	DTAUI,TAUCUMI,COSBI,WBARI,TAUGSURF,SEASHAZEFACT)
3
4	use radinc_h
5	use radcommon_h, only: gasi,gasi_recomb,tlimit,Cmk,gzlat_ig, &
6	tgasref,pfgasref,wnoi,scalep,indi
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,corrk_recombin, &
11	callclouds,callmufi,seashaze,uncoupl_optic_haze
12	use tracer_h, only : nmicro,nice
13
14	implicit none
15
16	!==================================================================
17	!
18	! Purpose
19	! -------
20	! Calculates longwave optical constants at each level. For each
21	! layer and spectral interval in the IR it calculates WBAR, DTAU
22	! and COSBAR. For each level it calculates TAU.
23	!
24	! TAUCUMI(L,LW) is the cumulative optical depth at level L (or alternatively
25	! at the bottom of layer L), LW is the spectral wavelength interval.
26	!
27	! TLEV(L) - Temperature at the layer boundary (i.e., level)
28	! PLEV(L) - Pressure at the layer boundary (i.e., level)
29	!
30	! Authors
31	! -------
32	! Adapted from the NASA Ames code by R. Wordsworth (2009)
33	! Clean and adaptation to Titan by J. Vatant d'Ollone (2016-17)
34	!
35	!==================================================================
36
37
38	!==========================================================
39	! Input/Output
40	!==========================================================
41	REAL*8, INTENT(IN) :: PQMO(nlay,nmicro) ! Tracers for microphysics optics (X/m2).
42	INTEGER, INTENT(IN) :: NLAY ! Number of pressure layers (for pqmo)
43	REAL*8, INTENT(IN) :: PLEV(L_LEVELS), TLEV(L_LEVELS)
44	REAL*8, INTENT(IN) :: TMID(L_LEVELS), PMID(L_LEVELS)
45	REAL*8, INTENT(IN) :: SEASHAZEFACT(L_LEVELS)
46
47	REAL*8, INTENT(OUT) :: DTAUI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
48	REAL*8, INTENT(OUT) :: TAUCUMI(L_LEVELS,L_NSPECTI,L_NGAUSS)
49	REAL*8, INTENT(OUT) :: COSBI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
50	REAL*8, INTENT(OUT) :: WBARI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
51	REAL*8, INTENT(OUT) :: TAUGSURF(L_NSPECTI,L_NGAUSS-1)
52	! ==========================================================
53
54	real*8 DTAUKI(L_LEVELS,L_NSPECTI,L_NGAUSS)
55
56	! Titan customisation
57	! J. Vatant d'Ollone (2016)
58	real*8 DHAZE_T(L_LEVELS,L_NSPECTI)
59	real*8 DHAZES_T(L_LEVELS,L_NSPECTI)
60	real*8 SSA_T(L_LEVELS,L_NSPECTI)
61	real*8 ASF_T(L_LEVELS,L_NSPECTI)
62	! ==========================
63
64	integer L, NW, NG, K, LK, IAER
65	integer MT(L_LEVELS), MP(L_LEVELS), NP(L_LEVELS)
66	real*8 ANS, TAUGAS
67	real*8 DPR(L_LEVELS), U(L_LEVELS)
68	real*8 LCOEF(4), LKCOEF(L_LEVELS,4)
69
70	real*8 DCONT
71	double precision wn_cont, p_cont, p_air, T_cont, dtemp, dtempc
72	double precision p_cross
73
74	real*8 KCOEF(4)
75
76	! temporary variable to reduce memory access time to gasi
77	real*8 tmpk(2,2)
78
79	! temporary variables for multiple aerosol calculation
80	real*8 atemp
81	real*8 btemp(L_NLAYRAD,L_NSPECTI)
82
83	! variables for k in units m^-1
84	real*8 dz(L_LEVELS)
85
86	integer igas, jgas, ilay
87
88	integer interm
89
90	! Variables for haze optics
91	character(len=200) file_path
92	logical file_ok
93	integer dumch
94	real*8 dumwvl
95
96	real*8 m3as,m3af
97	real*8 dtauaer_s,dtauaer_f
98	real*8,save :: rhoaer_s(L_NSPECTI),ssa_s(L_NSPECTI),asf_s(L_NSPECTI)
99	real*8,save :: rhoaer_f(L_NSPECTI),ssa_f(L_NSPECTI),asf_f(L_NSPECTI)
100	!$OMP THREADPRIVATE(rhoaer_s,rhoaer_f,ssa_s,ssa_f,asf_s,asf_f)
101
102	logical,save :: firstcall=.true.
103	!$OMP THREADPRIVATE(firstcall)
104
105
106	!! AS: to save time in computing continuum (see bilinearbig)
107	IF (.not.ALLOCATED(indi)) THEN
108	ALLOCATE(indi(L_NSPECTI,ngasmx,ngasmx))
109	indi = -9999 ! this initial value means "to be calculated"
110	ENDIF
111
112	! Some initialisation because there's a pb with disr_haze at the limits (nw=1)
113	! I should check this - For now we set vars to zero : better than nans - JVO 2017
114	DHAZE_T(:,:) = 0.0
115	SSA_T(:,:) = 0.0
116	ASF_T(:,:) = 0.0
117
118	! Load tabulated haze optical properties if needed.
119	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
120	IF (firstcall .AND. callmufi .AND. (.NOT. uncoupl_optic_haze)) THEN
121	OPEN(12,file=TRIM(haze_opt_file),form='formatted') ! The file has been inquired in physiq_mod firstcall
122	READ(12,*) ! dummy header
123	DO NW=1,L_NSPECTI
124	READ(12,*) dumch, dumwvl, rhoaer_f(nw), ssa_f(nw), asf_f(nw), rhoaer_s(nw), ssa_s(nw), asf_s(nw)
125	ENDDO
126	CLOSE(12)
127	ENDIF
128
129	!=======================================================================
130	! Determine the total gas opacity throughout the column, for each
131	! spectral interval, NW, and each Gauss point, NG.
132
133	taugsurf(:,:) = 0.0
134	dpr(:) = 0.0
135	lkcoef(:,:) = 0.0
136
137	do K=2,L_LEVELS
138
139	ilay = L_NLAYRAD+1 - k/2 ! int. arithmetic => gives the gcm layer index (reversed)
140
141	DPR(k) = PLEV(K)-PLEV(K-1)
142
143	! if we have continuum opacities, we need dz
144	dz(k) = dpr(k)RTMID(K)/(gzlat_ig(ilay)*PMID(K))
145	U(k) = Cmk(ilay)*DPR(k) ! only Cmk line in optci.F
146
147	call tpindex(PMID(K),TMID(K),pfgasref,tgasref,LCOEF,MT(K),MP(K))
148
149	do LK=1,4
150	LKCOEF(K,LK) = LCOEF(LK)
151	end do
152	end do ! levels
153
154	do NW=1,L_NSPECTI
155
156	do K=2,L_LEVELS
157
158	ilay = L_NLAYRAD+1 - k/2 ! int. arithmetic => gives the gcm layer index (reversed)
159
160	IF (callmufi .AND. (.NOT. uncoupl_optic_haze)) THEN
161	m3as = pqmo(ilay,2) / 2.0
162	m3af = pqmo(ilay,4) / 2.0
163
164	IF ( ilay .lt. 18 ) THEN
165	m3as = pqmo(18,2) / 2.0 *dz(k)/dz(18)
166	m3af = pqmo(18,4) / 2.0 *dz(k)/dz(18)
167	ENDIF
168
169	dtauaer_s = m3as*rhoaer_s(nw)
170	dtauaer_f = m3af*rhoaer_f(nw)
171	DHAZE_T(k,nw) = dtauaer_s + dtauaer_f
172
173	IF ( dtauaer_s + dtauaer_f .GT. 1.D-30 ) THEN
174	SSA_T(k,nw) = ( dtauaer_sssa_s(nw) + dtauaer_fssa_f(nw) ) / ( dtauaer_s+dtauaer_f )
175	ASF_T(k,nw) = ( dtauaer_sssa_s(nw)asf_s(nw) + dtauaer_fssa_f(nw)asf_f(nw) ) &
176	/ ( ssa_s(nw)dtauaer_s + ssa_f(nw)dtauaer_f )
177	ELSE
178	DHAZE_T(k,nw) = 0.D0
179	SSA_T(k,nw) = 1.0
180	ASF_T(k,nw) = 1.0
181	ENDIF
182
183	IF (callclouds.and.firstcall) &
184	WRITE(,) 'WARNING: In optci, optical properties &
185	&calculations are not implemented yet'
186	ELSE
187	! Call fixed vertical haze profile of extinction - same for all columns
188	call disr_haze(dz(k),plev(k),wnoi(nw),DHAZE_T(k,nw),SSA_T(k,nw),ASF_T(k,nw))
189	if (seashaze) DHAZE_T(k,nw) = DHAZE_T(k,nw)*seashazefact(k)
190	ENDIF
191
192	DCONT = 0.0d0 ! continuum absorption
193
194	if(continuum.and.(.not.graybody))then
195	! include continua if necessary
196	wn_cont = dble(wnoi(nw))
197	T_cont = dble(TMID(k))
198	do igas=1,ngasmx
199
200	p_cont = dble(PMID(k)scalepgfrac(igas,ilay))
201
202	dtemp=0.0d0
203	if(igas.eq.igas_N2)then
204
205	interm = indi(nw,igas,igas)
206	call interpolateN2N2(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
207	indi(nw,igas,igas) = interm
208
209	elseif(igas.eq.igas_H2)then
210
211	! first do self-induced absorption
212	interm = indi(nw,igas,igas)
213	call interpolateH2H2(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
214	indi(nw,igas,igas) = interm
215
216	! then cross-interactions with other gases
217	do jgas=1,ngasmx
218	p_cross = dble(PMID(k)scalepgfrac(jgas,ilay))
219	dtempc = 0.0d0
220	if(jgas.eq.igas_N2)then
221	interm = indi(nw,igas,jgas)
222	call interpolateN2H2(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
223	indi(nw,igas,jgas) = interm
224	endif
225	dtemp = dtemp + dtempc
226	enddo
227
228	elseif(igas.eq.igas_CH4)then
229
230	! first do self-induced absorption
231	interm = indi(nw,igas,igas)
232	call interpolateCH4CH4(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
233	indi(nw,igas,igas) = interm
234
235	! then cross-interactions with other gases
236	do jgas=1,ngasmx
237	p_cross = dble(PMID(k)scalepgfrac(jgas,ilay))
238	dtempc = 0.0d0
239	if(jgas.eq.igas_N2)then
240	interm = indi(nw,igas,jgas)
241	call interpolateN2CH4(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
242	indi(nw,igas,jgas) = interm
243	endif
244	dtemp = dtemp + dtempc
245	enddo
246
247	endif
248
249	DCONT = DCONT + dtemp
250
251	enddo
252
253	DCONT = DCONT*dz(k)
254
255	endif
256
257	do ng=1,L_NGAUSS-1
258
259	! Now compute TAUGAS
260
261	! JVO 2017 : added tmpk because the repeated calls to gasi/v increased dramatically
262	! the execution time of optci/v -> ~ factor 2 on the whole radiative
263	! transfer on the tested simulations !
264
265	if (corrk_recombin) then
266	tmpk = GASI_RECOMB(MT(K):MT(K)+1,MP(K):MP(K)+1,NW,NG)
267	else
268	tmpk = GASI(MT(K):MT(K)+1,MP(K):MP(K)+1,1,NW,NG)
269	endif
270
271	KCOEF(1) = tmpk(1,1) ! KCOEF(1) = GASI(MT(K),MP(K),1,NW,NG)
272	KCOEF(2) = tmpk(1,2) ! KCOEF(2) = GASI(MT(K),MP(K)+1,1,NW,NG)
273	KCOEF(3) = tmpk(2,2) ! KCOEF(3) = GASI(MT(K)+1,MP(K)+1,1,NW,NG)
274	KCOEF(4) = tmpk(2,1) ! KCOEF(4) = GASI(MT(K)+1,MP(K),1,NW,NG)
275
276
277	! Interpolate the gaseous k-coefficients to the requested T,P values
278
279	ANS = LKCOEF(K,1)KCOEF(1) + LKCOEF(K,2)KCOEF(2) + &
280	LKCOEF(K,3)KCOEF(3) + LKCOEF(K,4)KCOEF(4)
281
282	TAUGAS = U(k)*ANS
283
284	TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS + DCONT
285	DTAUKI(K,nw,ng) = TAUGAS &
286	+ DCONT & ! For parameterized continuum absorption
287	+ DHAZE_T(K,NW) ! For Titan haze
288
289	end do
290
291	! Now fill in the "clear" part of the spectrum (NG = L_NGAUSS),
292	! which holds continuum opacity only
293
294	NG = L_NGAUSS
295	DTAUKI(K,nw,ng) = 0.d0 &
296	+ DCONT & ! For parameterized continuum absorption
297	+ DHAZE_T(K,NW) ! For Titan Haze
298
299	end do
300	end do
301
302	!=======================================================================
303	! Now the full treatment for the layers, where besides the opacity
304	! we need to calculate the scattering albedo and asymmetry factors
305	! ======================================================================
306
307	! Haze scattering
308	DO NW=1,L_NSPECTI
309	DO K=2,L_LEVELS
310	DHAZES_T(K,NW) = DHAZE_T(K,NW) * SSA_T(K,NW)
311	ENDDO
312	ENDDO
313
314	DO NW=1,L_NSPECTI
315	DO L=1,L_NLAYRAD-1
316	K = 2*L+1
317	btemp(L,NW) = DHAZES_T(K,NW) + DHAZES_T(K+1,NW)
318	END DO ! L vertical loop
319
320	! Last level
321	L = L_NLAYRAD
322	K = 2*L+1
323	btemp(L,NW) = DHAZES_T(K,NW)
324
325	END DO ! NW spectral loop
326
327
328	DO NW=1,L_NSPECTI
329	NG = L_NGAUSS
330	DO L=1,L_NLAYRAD-1
331
332	K = 2*L+1
333	DTAUI(L,nw,ng) = DTAUKI(K,NW,NG) + DTAUKI(K+1,NW,NG)! + 1.e-50
334
335	atemp = 0.
336	if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
337	atemp = atemp + &
338	ASF_T(K,NW)*DHAZES_T(K,NW) + &
339	ASF_T(K+1,NW)*DHAZES_T(K+1,NW)
340
341	WBARI(L,nw,ng) = btemp(L,nw) / DTAUI(L,NW,NG)
342	else
343	WBARI(L,nw,ng) = 0.0D0
344	DTAUI(L,NW,NG) = 1.0D-9
345	endif
346
347	if(btemp(L,nw) .GT. 0.0d0) then
348	cosbi(L,NW,NG) = atemp/btemp(L,nw)
349	else
350	cosbi(L,NW,NG) = 0.0D0
351	end if
352
353	END DO ! L vertical loop
354
355	! Last level
356
357	L = L_NLAYRAD
358	K = 2*L+1
359	DTAUI(L,nw,ng) = DTAUKI(K,NW,NG) ! + 1.e-50
360
361	atemp = 0.
362	if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
363	atemp = atemp + ASF_T(K,NW)*DHAZES_T(K,NW)
364	WBARI(L,nw,ng) = btemp(L,nw) / DTAUI(L,NW,NG)
365	else
366	WBARI(L,nw,ng) = 0.0D0
367	DTAUI(L,NW,NG) = 1.0D-9
368	endif
369
370	if(btemp(L,nw) .GT. 0.0d0) then
371	cosbi(L,NW,NG) = atemp/btemp(L,nw)
372	else
373	cosbi(L,NW,NG) = 0.0D0
374	end if
375
376
377	! Now the other Gauss points, if needed.
378
379	DO NG=1,L_NGAUSS-1
380	IF(TAUGSURF(NW,NG) .gt. TLIMIT) THEN
381
382	DO L=1,L_NLAYRAD-1
383	K = 2*L+1
384	DTAUI(L,nw,ng) = DTAUKI(K,NW,NG)+DTAUKI(K+1,NW,NG)! + 1.e-50
385
386	if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
387
388	WBARI(L,nw,ng) = btemp(L,nw) / DTAUI(L,NW,NG)
389
390	else
391	WBARI(L,nw,ng) = 0.0D0
392	DTAUI(L,NW,NG) = 1.0D-9
393	endif
394
395	cosbi(L,NW,NG) = cosbi(L,NW,L_NGAUSS)
396	END DO ! L vertical loop
397
398	! Last level
399	L = L_NLAYRAD
400	K = 2*L+1
401	DTAUI(L,nw,ng) = DTAUKI(K,NW,NG)! + 1.e-50
402
403	if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
404
405	WBARI(L,nw,ng) = btemp(L,nw) / DTAUI(L,NW,NG)
406
407	else
408	WBARI(L,nw,ng) = 0.0D0
409	DTAUI(L,NW,NG) = 1.0D-9
410	endif
411
412	cosbi(L,NW,NG) = cosbi(L,NW,L_NGAUSS)
413
414	END IF
415
416	END DO ! NG Gauss loop
417	END DO ! NW spectral loop
418
419	! Total extinction optical depths
420
421	DO NG=1,L_NGAUSS ! full gauss loop
422	DO NW=1,L_NSPECTI
423	TAUCUMI(1,NW,NG)=0.0D0
424	DO K=2,L_LEVELS
425	TAUCUMI(K,NW,NG)=TAUCUMI(K-1,NW,NG)+DTAUKI(K,NW,NG)
426	END DO
427	END DO ! end full gauss loop
428	END DO
429
430	! be aware when comparing with textbook results
431	! (e.g. Pierrehumbert p. 218) that
432	! taucumi does not take the <cos theta>=0.5 factor into
433	! account. It is the optical depth for a vertically
434	! ascending ray with angle theta = 0.
435
436	if(firstcall) firstcall = .false.
437
438	return
439
440
441	end subroutine optci
442
443
444

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

Download in other formats:

Original Format