Context Navigation

optci.F90 @ 1899

Last change on this file since 1899 was 1897, checked in by jvatant, 7 years ago

Making Titan's hazy again - part II
+ Major updates of J.Burgalat YAMMS library and optical coupling, including :
++ Added the routines for haze optics inside YAMMS
++ Calling rad. transf. with interactive haze is plugged
in but should stay unactive as long as the microphysics is
in test phase : cf "uncoupl_optic_haze" flag : true for now !
++ Also some sanity checks for negative tendencies and
some others upkeep of YAMMS model
+ Also added a temporary CPP key USE_QTEST in physiq_mod
that enables to have microphysical tendencies separated
from dynamics for debugging and test phases
-- JVO and JB

File size: 14.0 KB

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

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

Download in other formats:

Original Format