Context Navigation

← Previous Revision
Next Revision →
Normal
Revision Log

optci.F90

Last change on this file was 3497, checked in by debatzbr, 2 weeks ago
Add AC6H6 condensation in the microphysics
File size: 22.2 KB

Rev	Line
[3083]	1	subroutine optci(PQMO,NLAY,ZLEV,PLEV,TLEV,TMID,PMID, &
	2	DTAUI,TAUCUMI,COSBI,WBARI,TAUGSURF,SEASHAZEFACT,&
[3318]	3	DIAG_OPTH,DIAG_OPTT,CDCOLUMN)
[135]	4
[716]	5	use radinc_h
[2050]	6	use radcommon_h, only: gasi,gasi_recomb,tlimit,Cmk,gzlat_ig, &
[2133]	7	tgasref,pfgasref,wnoi,scalep,indi
[716]	8	use gases_h
[2242]	9	use datafile_mod, only: haze_opt_file
[3083]	10	use comcstfi_mod, only: pi,r
	11	use callkeys_mod, only: continuum,graybody,corrk_recombin, &
	12	callclouds,callmufi,seashaze,uncoupl_optic_haze,&
[3318]	13	opt4clouds,FHIR,FCIR
[3083]	14	use tracer_h, only: nmicro,nice,ices_indx
[1897]	15
[716]	16	implicit none
[135]	17
[716]	18	!==================================================================
	19	!
	20	! Purpose
	21	! -------
	22	! Calculates longwave optical constants at each level. For each
	23	! layer and spectral interval in the IR it calculates WBAR, DTAU
	24	! and COSBAR. For each level it calculates TAU.
	25	!
[1822]	26	! TAUCUMI(L,LW) is the cumulative optical depth at level L (or alternatively
[716]	27	! at the bottom of layer L), LW is the spectral wavelength interval.
	28	!
	29	! TLEV(L) - Temperature at the layer boundary (i.e., level)
	30	! PLEV(L) - Pressure at the layer boundary (i.e., level)
	31	!
	32	! Authors
	33	! -------
	34	! Adapted from the NASA Ames code by R. Wordsworth (2009)
	35	!
[3090]	36	! Modified
	37	! --------
	38	! J. Vatant d'Ollone (2016-17)
	39	! --> Clean and adaptation to Titan
	40	! B. de Batz de Trenquelléon (2022-2023)
[3497]	41	! --> Clean and correction
[3090]	42	! --> New optics added for Titan's clouds
	43	!
[716]	44	!==================================================================
[135]	45
	46
[1822]	47	!==========================================================
	48	! Input/Output
	49	!==========================================================
[2050]	50	REAL*8, INTENT(IN) :: PQMO(nlay,nmicro) ! Tracers for microphysics optics (X/m2).
	51	INTEGER, INTENT(IN) :: NLAY ! Number of pressure layers (for pqmo)
[3083]	52	REAL*8, INTENT(IN) :: ZLEV(NLAY+1)
[1822]	53	REAL*8, INTENT(IN) :: PLEV(L_LEVELS), TLEV(L_LEVELS)
	54	REAL*8, INTENT(IN) :: TMID(L_LEVELS), PMID(L_LEVELS)
[2046]	55	REAL*8, INTENT(IN) :: SEASHAZEFACT(L_LEVELS)
[3083]	56	INTEGER, INTENT(IN) :: CDCOLUMN
[1822]	57
	58	REAL*8, INTENT(OUT) :: DTAUI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
	59	REAL*8, INTENT(OUT) :: TAUCUMI(L_LEVELS,L_NSPECTI,L_NGAUSS)
	60	REAL*8, INTENT(OUT) :: COSBI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
	61	REAL*8, INTENT(OUT) :: WBARI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
	62	REAL*8, INTENT(OUT) :: TAUGSURF(L_NSPECTI,L_NGAUSS-1)
[3318]	63	REAL*8, INTENT(OUT) :: DIAG_OPTH(L_LEVELS,L_NSPECTI,6)
	64	REAL*8, INTENT(OUT) :: DIAG_OPTT(L_LEVELS,L_NSPECTI,6)
[1822]	65	! ==========================================================
	66
[1722]	67	real*8 DTAUKI(L_LEVELS,L_NSPECTI,L_NGAUSS)
[135]	68
[1648]	69	! Titan customisation
	70	! J. Vatant d'Ollone (2016)
[1722]	71	real*8 DHAZE_T(L_LEVELS,L_NSPECTI)
	72	real*8 DHAZES_T(L_LEVELS,L_NSPECTI)
	73	real*8 SSA_T(L_LEVELS,L_NSPECTI)
	74	real*8 ASF_T(L_LEVELS,L_NSPECTI)
[1648]	75	! ==========================
	76
[716]	77	integer L, NW, NG, K, LK, IAER
	78	integer MT(L_LEVELS), MP(L_LEVELS), NP(L_LEVELS)
	79	real*8 ANS, TAUGAS
	80	real*8 DPR(L_LEVELS), U(L_LEVELS)
	81	real*8 LCOEF(4), LKCOEF(L_LEVELS,4)
[135]	82
[1788]	83	real*8 DCONT
[716]	84	double precision wn_cont, p_cont, p_air, T_cont, dtemp, dtempc
	85	double precision p_cross
[135]	86
[716]	87	real*8 KCOEF(4)
[1725]	88
	89	! temporary variable to reduce memory access time to gasi
	90	real*8 tmpk(2,2)
[1648]	91
[716]	92	! temporary variables for multiple aerosol calculation
[918]	93	real*8 atemp
	94	real*8 btemp(L_NLAYRAD,L_NSPECTI)
[135]	95
[716]	96	! variables for k in units m^-1
[873]	97	real*8 dz(L_LEVELS)
[135]	98
[1648]	99	integer igas, jgas, ilay
[253]	100
[873]	101	integer interm
	102
[2242]	103	! Variables for haze optics
	104	character(len=200) file_path
	105	logical file_ok
	106	integer dumch
	107	real*8 dumwvl
[3340]	108	integer ilev_cutoff
	109	real*8 corr_haze
[2242]	110
[3083]	111	! Variables for new optics
	112	integer iq, iw, FTYPE, CTYPE
[3090]	113	real*8 m0as,m0af,m0ccn,m3as,m3af,m3ccn,m3cld
[3083]	114	real*8 dtauaer_s,dtauaer_f,dtau_ccn,dtau_cld
[2242]	115	real*8,save :: rhoaer_s(L_NSPECTI),ssa_s(L_NSPECTI),asf_s(L_NSPECTI)
	116	real*8,save :: rhoaer_f(L_NSPECTI),ssa_f(L_NSPECTI),asf_f(L_NSPECTI)
[3083]	117	real*8,save :: ssa_ccn(L_NSPECTI),asf_ccn(L_NSPECTI)
	118	real*8,save :: ssa_cld(L_NSPECTI),asf_cld(L_NSPECTI)
	119	!$OMP THREADPRIVATE(rhoaer_s,rhoaer_f,ssa_s,ssa_f,ssa_cld,asf_s,asf_f,asf_cld)
[2242]	120
[1897]	121	logical,save :: firstcall=.true.
[2242]	122	!$OMP THREADPRIVATE(firstcall)
[1897]	123
[2242]	124
[873]	125	!! AS: to save time in computing continuum (see bilinearbig)
	126	IF (.not.ALLOCATED(indi)) THEN
[878]	127	ALLOCATE(indi(L_NSPECTI,ngasmx,ngasmx))
[873]	128	indi = -9999 ! this initial value means "to be calculated"
	129	ENDIF
	130
[2242]	131	! Some initialisation because there's a pb with disr_haze at the limits (nw=1)
[1792]	132	! I should check this - For now we set vars to zero : better than nans - JVO 2017
[2242]	133	DHAZE_T(:,:) = 0.0
	134	SSA_T(:,:) = 0.0
	135	ASF_T(:,:) = 0.0
[1792]	136
[2242]	137	! Load tabulated haze optical properties if needed.
	138	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	139	IF (firstcall .AND. callmufi .AND. (.NOT. uncoupl_optic_haze)) THEN
	140	OPEN(12,file=TRIM(haze_opt_file),form='formatted') ! The file has been inquired in physiq_mod firstcall
	141	READ(12,*) ! dummy header
	142	DO NW=1,L_NSPECTI
	143	READ(12,*) dumch, dumwvl, rhoaer_f(nw), ssa_f(nw), asf_f(nw), rhoaer_s(nw), ssa_s(nw), asf_s(nw)
	144	ENDDO
	145	CLOSE(12)
	146	ENDIF
	147
[3083]	148	!=======================================================================
	149	! Determine the total gas opacity throughout the column, for each
	150	! spectral interval, NW, and each Gauss point, NG.
[135]	151
[3083]	152	taugsurf(:,:) = 0.0
	153	dpr(:) = 0.0
	154	lkcoef(:,:) = 0.0
[135]	155
[3340]	156	! Level of cutoff
	157	!~~~~~~~~~~~~~~~~
	158	ilev_cutoff = 26
	159
[3083]	160	do K=2,L_LEVELS
	161	ilay = L_NLAYRAD+1 - k/2 ! int. arithmetic => gives the gcm layer index (reversed)
	162	DPR(k) = PLEV(K)-PLEV(K-1)
[135]	163
[3083]	164	! if we have continuum opacities, we need dz
[1947]	165	dz(k) = dpr(k)RTMID(K)/(gzlat_ig(ilay)*PMID(K))
	166	U(k) = Cmk(ilay)*DPR(k) ! only Cmk line in optci.F
[3083]	167
	168	call tpindex(PMID(K),TMID(K),pfgasref,tgasref,LCOEF,MT(K),MP(K))
[135]	169
[3083]	170	do LK=1,4
	171	LKCOEF(K,LK) = LCOEF(LK)
	172	end do
	173	end do ! L_LEVELS
[253]	174
[3083]	175	do NW=1,L_NSPECTI
	176	! We ignore K=1...
	177	do K=2,L_LEVELS
	178	! int. arithmetic => gives the gcm layer index (reversed)
	179	ilay = L_NLAYRAD+1 - k/2
[135]	180
[3083]	181	! Optics coupled with the microphysics :
	182	IF (callmufi .AND. (.NOT. uncoupl_optic_haze)) THEN
[1722]	183
[3083]	184	!==========================================================================================
	185	! Old optics (must have callclouds = .false.):
	186	!==========================================================================================
	187	IF (.NOT. opt4clouds) THEN
	188	m3as = pqmo(ilay,2) / 2.0
	189	m3af = pqmo(ilay,4) / 2.0
[3340]	190	! Cut-off
	191	IF (ilay .lt. ilev_cutoff) THEN
	192	m3as = pqmo(ilev_cutoff,2) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(ilev_cutoff+1)-zlev(ilev_cutoff))
	193	m3af = pqmo(ilev_cutoff,4) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(ilev_cutoff+1)-zlev(ilev_cutoff))
[3083]	194	ENDIF
[1897]	195
[3083]	196	dtauaer_s = m3as*rhoaer_s(nw)
	197	dtauaer_f = m3af*rhoaer_f(nw)
	198
	199	!==========================================================================================
	200	! New optics :
	201	!==========================================================================================
	202	ELSE
	203	iw = (L_NSPECTI + 1) - NW + L_NSPECTV ! Visible first and return
	204	!-----------------------------
	205	! HAZE (Spherical + Fractal) :
	206	!-----------------------------
	207	FTYPE = 1
[1722]	208
[3083]	209	! Spherical aerosols :
	210	!---------------------
	211	CTYPE = 5
	212	m0as = pqmo(ilay,1) / 2.0
	213	m3as = pqmo(ilay,2) / 2.0
[3340]	214	! If not callclouds : must have a cut-off
[3083]	215	IF (.NOT. callclouds) THEN
[3340]	216	IF (ilay .lt. ilev_cutoff) THEN
	217	m0as = pqmo(ilev_cutoff,1) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(ilev_cutoff+1)-zlev(ilev_cutoff))
	218	m3as = pqmo(ilev_cutoff,2) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(ilev_cutoff+1)-zlev(ilev_cutoff))
[3083]	219	ENDIF
	220	ENDIF
	221	call get_haze_and_cloud_opacity(FTYPE,CTYPE,m0as,m3as,iw,dtauaer_s,ssa_s(nw),asf_s(nw))
	222
	223	! Fractal aerosols :
	224	!-------------------
	225	CTYPE = FTYPE
	226	m0af = pqmo(ilay,3) / 2.0
	227	m3af = pqmo(ilay,4) / 2.0
[3340]	228	! If not callclouds : must have a cut-off
[3083]	229	IF (.NOT. callclouds) THEN
[3340]	230	IF (ilay .lt. ilev_cutoff) THEN
	231	m0af = pqmo(ilev_cutoff,3) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(ilev_cutoff+1)-zlev(ilev_cutoff))
	232	m3af = pqmo(ilev_cutoff,4) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(ilev_cutoff+1)-zlev(ilev_cutoff))
[3083]	233	ENDIF
	234	ENDIF
	235	call get_haze_and_cloud_opacity(FTYPE,CTYPE,m0af,m3af,iw,dtauaer_f,ssa_f(nw),asf_f(nw))
	236	ENDIF
[135]	237
[3318]	238	! Tuning of optical properties for haze :
	239	!dtauaer_s = dtauaer_s * (FHIR * (1-ssa_s(nw)) + ssa_s(nw))
	240	!ssa_s(nw) = ssa_s(nw) / (FHIR * (1-ssa_s(nw)) + ssa_s(nw))
	241	dtauaer_f = dtauaer_f * (FHIR * (1-ssa_f(nw)) + ssa_f(nw))
	242	ssa_f(nw) = ssa_f(nw) / (FHIR * (1-ssa_f(nw)) + ssa_f(nw))
	243
[3083]	244	! Total of Haze opacity (dtau), SSA (w) and ASF (COS) :
	245	DHAZE_T(k,nw) = dtauaer_s + dtauaer_f
	246	IF (dtauaer_s + dtauaer_f .GT. 1.D-30) THEN
	247	SSA_T(k,nw) = ( dtauaer_sssa_s(nw) + dtauaer_fssa_f(nw) ) / ( dtauaer_s+dtauaer_f )
	248	ASF_T(k,nw) = ( dtauaer_sssa_s(nw)asf_s(nw) + dtauaer_fssa_f(nw)asf_f(nw) ) &
	249	/ ( ssa_s(nw)dtauaer_s + ssa_f(nw)dtauaer_f )
	250	ELSE
	251	DHAZE_T(k,nw) = 0.D0
	252	SSA_T(k,nw) = 1.0
	253	ASF_T(k,nw) = 1.0
	254	ENDIF
[3318]	255
[3340]	256	! Opacity and albedo adjustement below cutoff :
	257	IF (.NOT. callclouds) THEN
	258	corr_haze=0.6-0.4TANH((PMID(K)100.-2500.)/250.)
	259	IF (ilay .lt. ilev_cutoff) THEN
	260	DHAZE_T(k,nw) = DHAZE_T(k,nw) * corr_haze
	261	ENDIF
	262	ENDIF
	263
[3083]	264	! Diagnostics for the haze :
[3318]	265	DIAG_OPTH(k,nw,1) = DHAZE_T(k,nw) ! dtau
	266	DIAG_OPTH(k,nw,2) = SSA_T(k,nw) ! wbar
	267	DIAG_OPTH(k,nw,3) = ASF_T(k,nw) ! gbar
[253]	268
[3083]	269	!---------------------
	270	! CLOUDS (Spherical) :
	271	!---------------------
	272	IF (callclouds) THEN
	273	CTYPE = 0
	274	m0ccn = pqmo(ilay,5) / 2.0
	275	m3ccn = pqmo(ilay,6) / 2.0
[3318]	276	m3cld = pqmo(ilay,6) / 2.0
[3083]	277
	278	! Clear / Dark column method :
	279	!-----------------------------
[305]	280
[3083]	281	! CCN's SSA :
	282	call get_haze_and_cloud_opacity(FTYPE,FTYPE,m0ccn,m3ccn,iw,dtau_ccn,ssa_ccn(nw),asf_ccn(nw))
	283
[3497]	284	! Clear column (CCN, C2H2, C2H6, HCN, AC6H6) :
[3083]	285	IF (CDCOLUMN == 0) THEN
	286	DO iq = 2, nice
	287	m3cld = m3cld + (pqmo(ilay,ices_indx(iq)) / 2.0)
	288	ENDDO
	289	call get_haze_and_cloud_opacity(FTYPE,CTYPE,m0ccn,m3cld,iw,dtau_cld,ssa_cld(nw),asf_cld(nw))
	290
[3497]	291	! Dark column (CCN, CH4, C2H2, C2H6, HCN, AC6H6) :
[3083]	292	ELSEIF (CDCOLUMN == 1) THEN
	293	DO iq = 1, nice
	294	m3cld = m3cld + (pqmo(ilay,ices_indx(iq)) / 2.0)
	295	ENDDO
	296	call get_haze_and_cloud_opacity(FTYPE,CTYPE,m0ccn,m3cld,iw,dtau_cld,ssa_cld(nw),asf_cld(nw))
	297
	298	ELSE
	299	WRITE(,) 'WARNING OPTCI.F90 : CDCOLUMN MUST BE 0 OR 1'
	300	WRITE(,) 'WE USE DARK COLUMN ...'
	301	DO iq = 1, nice
	302	m3cld = m3cld + (pqmo(ilay,ices_indx(iq)) / 2.0)
	303	ENDDO
	304	call get_haze_and_cloud_opacity(FTYPE,CTYPE,m0ccn,m3cld,iw,dtau_cld,ssa_cld(nw),asf_cld(nw))
	305	ENDIF
	306
	307	! For small dropplets, opacity of nucleus dominates...
[3318]	308	dtau_cld = (dtau_cldm3ccn + dtau_cldm3cld) / (m3ccn + m3cld)
[3083]	309	ssa_cld(nw) = (ssa_ccn(nw)m3ccn + ssa_cld(nw)m3cld) / (m3ccn + m3cld)
[3318]	310
	311	! Tuning of optical properties for clouds :
	312	dtau_cld = dtau_cld * (FCIR * (1-ssa_cld(nw)) + ssa_cld(nw))
	313	ssa_cld(nw) = ssa_cld(nw) / (FCIR * (1-ssa_cld(nw)) + ssa_cld(nw))
	314
[3083]	315	! Total of Haze + Clouds opacity (dtau), SSA (w) and ASF (COS) :
	316	DHAZE_T(k,nw) = dtauaer_s + dtauaer_f + dtau_cld
	317	IF (DHAZE_T(k,nw) .GT. 1.D-30) THEN
	318	SSA_T(k,nw) = ( dtauaer_sssa_s(nw) + dtauaer_fssa_f(nw) + dtau_cld*ssa_cld(nw) ) / ( dtauaer_s+dtauaer_f+dtau_cld )
	319	ASF_T(k,nw) = ( dtauaer_sssa_s(nw)asf_s(nw) + dtauaer_fssa_f(nw)asf_f(nw) + dtau_cldssa_cld(nw)asf_cld(nw) ) &
	320	/ ( ssa_s(nw)dtauaer_s + ssa_f(nw)dtauaer_f + ssa_cld(nw)*dtau_cld )
	321	ELSE
	322	DHAZE_T(k,nw) = 0.D0
	323	SSA_T(k,nw) = 1.0
	324	ASF_T(k,nw) = 1.0
	325	ENDIF
[253]	326
[3083]	327	! Diagnostics for clouds :
[3318]	328	DIAG_OPTT(k,nw,1) = DHAZE_T(k,nw) ! dtau
	329	DIAG_OPTT(k,nw,2) = SSA_T(k,nw) ! wbar
	330	DIAG_OPTT(k,nw,3) = ASF_T(k,nw) ! gbar
[135]	331
[3083]	332	ELSE
	333	! Diagnostics for clouds :
[3318]	334	DIAG_OPTT(k,nw,1) = 0.D0 ! dtau
	335	DIAG_OPTT(k,nw,2) = 1.0 ! wbar
	336	DIAG_OPTT(k,nw,3) = 1.0 ! gbar
[3083]	337	ENDIF
	338
	339	! Optics and microphysics no coupled :
	340	ELSE
	341	! Call fixed vertical haze profile of extinction - same for all columns
	342	call disr_haze(dz(k),plev(k),wnoi(nw),DHAZE_T(k,nw),SSA_T(k,nw),ASF_T(k,nw))
	343	if (seashaze) DHAZE_T(k,nw) = DHAZE_T(k,nw)*seashazefact(k)
	344	! Diagnostics for the haze :
[3318]	345	DIAG_OPTH(k,nw,1) = DHAZE_T(k,nw) ! dtau
	346	DIAG_OPTH(k,nw,2) = SSA_T(k,nw) ! wbar
	347	DIAG_OPTH(k,nw,3) = ASF_T(k,nw) ! gbar
[3083]	348	! Diagnostics for clouds :
[3318]	349	DIAG_OPTT(k,nw,1) = 0.D0 ! dtau
	350	DIAG_OPTT(k,nw,2) = 1.0 ! wbar
	351	DIAG_OPTT(k,nw,3) = 1.0 ! gbar
[3083]	352	ENDIF ! ENDIF callmufi
	353
	354	DCONT = 0.0d0 ! continuum absorption
[1648]	355
[3083]	356	if(continuum.and.(.not.graybody))then
	357	! include continua if necessary
	358	wn_cont = dble(wnoi(nw))
	359	T_cont = dble(TMID(k))
	360	do igas=1,ngasmx
[1648]	361
[3083]	362	p_cont = dble(PMID(k)scalepgfrac(igas,ilay))
[1648]	363
[3083]	364	dtemp=0.0d0
	365	if(igas.eq.igas_N2)then
[135]	366
[3083]	367	interm = indi(nw,igas,igas)
	368	call interpolateN2N2(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
	369	indi(nw,igas,igas) = interm
[135]	370
[3083]	371	elseif(igas.eq.igas_H2)then
[135]	372
[3083]	373	! first do self-induced absorption
	374	interm = indi(nw,igas,igas)
	375	call interpolateH2H2(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
	376	indi(nw,igas,igas) = interm
[135]	377
[3083]	378	! then cross-interactions with other gases
	379	do jgas=1,ngasmx
	380	p_cross = dble(PMID(k)scalepgfrac(jgas,ilay))
	381	dtempc = 0.0d0
	382	if(jgas.eq.igas_N2)then
	383	interm = indi(nw,igas,jgas)
	384	call interpolateN2H2(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
	385	indi(nw,igas,jgas) = interm
	386	endif
	387	dtemp = dtemp + dtempc
	388	enddo
[135]	389
[3083]	390	elseif(igas.eq.igas_CH4)then
[135]	391
[3083]	392	! first do self-induced absorption
	393	interm = indi(nw,igas,igas)
	394	call interpolateCH4CH4(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
	395	indi(nw,igas,igas) = interm
[135]	396
[3083]	397	! then cross-interactions with other gases
	398	do jgas=1,ngasmx
	399	p_cross = dble(PMID(k)scalepgfrac(jgas,ilay))
	400	dtempc = 0.0d0
	401	if(jgas.eq.igas_N2)then
	402	interm = indi(nw,igas,jgas)
	403	call interpolateN2CH4(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
	404	indi(nw,igas,jgas) = interm
	405	endif
	406	dtemp = dtemp + dtempc
	407	enddo
[253]	408
[3083]	409	endif
[1725]	410
[3083]	411	DCONT = DCONT + dtemp
[1725]	412
[3083]	413	enddo
[1725]	414
[3083]	415	DCONT = DCONT*dz(k)
[135]	416
[3083]	417	endif
[135]	418
[3083]	419	do ng=1,L_NGAUSS-1
[135]	420
[3083]	421	! Now compute TAUGAS
[135]	422
[3083]	423	! JVO 2017 : added tmpk because the repeated calls to gasi/v increased dramatically
	424	! the execution time of optci/v -> ~ factor 2 on the whole radiative
	425	! transfer on the tested simulations !
[135]	426
[3083]	427	if (corrk_recombin) then
	428	tmpk = GASI_RECOMB(MT(K):MT(K)+1,MP(K):MP(K)+1,NW,NG)
	429	else
	430	tmpk = GASI(MT(K):MT(K)+1,MP(K):MP(K)+1,1,NW,NG)
	431	endif
[135]	432
[3083]	433	KCOEF(1) = tmpk(1,1) ! KCOEF(1) = GASI(MT(K),MP(K),1,NW,NG)
	434	KCOEF(2) = tmpk(1,2) ! KCOEF(2) = GASI(MT(K),MP(K)+1,1,NW,NG)
	435	KCOEF(3) = tmpk(2,2) ! KCOEF(3) = GASI(MT(K)+1,MP(K)+1,1,NW,NG)
	436	KCOEF(4) = tmpk(2,1) ! KCOEF(4) = GASI(MT(K)+1,MP(K),1,NW,NG)
[135]	437
	438
[3083]	439	! Interpolate the gaseous k-coefficients to the requested T,P values
	440
	441	ANS = LKCOEF(K,1)KCOEF(1) + LKCOEF(K,2)KCOEF(2) + &
	442	LKCOEF(K,3)KCOEF(3) + LKCOEF(K,4)KCOEF(4)
	443
	444	TAUGAS = U(k)*ANS
	445
	446	TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS + DCONT
	447	DTAUKI(K,nw,ng) = TAUGAS &
	448	+ DCONT & ! For parameterized continuum absorption
	449	+ DHAZE_T(K,NW) ! For Titan haze
	450
	451	end do
	452
	453	! Now fill in the "clear" part of the spectrum (NG = L_NGAUSS),
	454	! which holds continuum opacity only
	455
	456	NG = L_NGAUSS
	457	DTAUKI(K,nw,ng) = 0.d0 &
	458	+ DCONT & ! For parameterized continuum absorption
	459	+ DHAZE_T(K,NW) ! For Titan Haze
	460
[3318]	461	DIAG_OPTH(K,nw,4) = 0.d0
	462	DIAG_OPTH(K,nw,5) = TAUGAS
	463	DIAG_OPTH(K,nw,6) = DCONT
	464	DIAG_OPTT(K,nw,4) = 0.d0
	465	DIAG_OPTT(K,nw,5) = TAUGAS
	466	DIAG_OPTT(K,nw,6) = DCONT
	467
[3083]	468	end do ! k = L_LEVELS
	469	end do ! nw = L_NSPECTI
	470
[716]	471	!=======================================================================
	472	! Now the full treatment for the layers, where besides the opacity
	473	! we need to calculate the scattering albedo and asymmetry factors
[1648]	474	! ======================================================================
[135]	475
[1648]	476	! Haze scattering
[918]	477	DO NW=1,L_NSPECTI
[1722]	478	DO K=2,L_LEVELS
[1648]	479	DHAZES_T(K,NW) = DHAZE_T(K,NW) * SSA_T(K,NW)
	480	ENDDO
	481	ENDDO
	482
	483	DO NW=1,L_NSPECTI
	484	DO L=1,L_NLAYRAD-1
[918]	485	K = 2*L+1
[1788]	486	btemp(L,NW) = DHAZES_T(K,NW) + DHAZES_T(K+1,NW)
[918]	487	END DO ! L vertical loop
[1648]	488
[1722]	489	! Last level
	490	L = L_NLAYRAD
	491	K = 2*L+1
[1788]	492	btemp(L,NW) = DHAZES_T(K,NW)
[1648]	493
[918]	494	END DO ! NW spectral loop
	495
[135]	496
[716]	497	DO NW=1,L_NSPECTI
	498	NG = L_NGAUSS
[1648]	499	DO L=1,L_NLAYRAD-1
[135]	500
[716]	501	K = 2*L+1
	502	DTAUI(L,nw,ng) = DTAUKI(K,NW,NG) + DTAUKI(K+1,NW,NG)! + 1.e-50
[135]	503
[716]	504	atemp = 0.
[961]	505	if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
[1722]	506	atemp = atemp + &
	507	ASF_T(K,NW)*DHAZES_T(K,NW) + &
	508	ASF_T(K+1,NW)*DHAZES_T(K+1,NW)
	509
[918]	510	WBARI(L,nw,ng) = btemp(L,nw) / DTAUI(L,NW,NG)
[716]	511	else
	512	WBARI(L,nw,ng) = 0.0D0
[961]	513	DTAUI(L,NW,NG) = 1.0D-9
[716]	514	endif
[135]	515
[961]	516	if(btemp(L,nw) .GT. 0.0d0) then
[918]	517	cosbi(L,NW,NG) = atemp/btemp(L,nw)
[716]	518	else
	519	cosbi(L,NW,NG) = 0.0D0
	520	end if
[135]	521
[716]	522	END DO ! L vertical loop
[1648]	523
[1722]	524	! Last level
[1648]	525
[1722]	526	L = L_NLAYRAD
	527	K = 2*L+1
	528	DTAUI(L,nw,ng) = DTAUKI(K,NW,NG) ! + 1.e-50
[1648]	529
[1722]	530	atemp = 0.
	531	if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
	532	atemp = atemp + ASF_T(K,NW)*DHAZES_T(K,NW)
	533	WBARI(L,nw,ng) = btemp(L,nw) / DTAUI(L,NW,NG)
	534	else
	535	WBARI(L,nw,ng) = 0.0D0
	536	DTAUI(L,NW,NG) = 1.0D-9
	537	endif
[1648]	538
[1722]	539	if(btemp(L,nw) .GT. 0.0d0) then
	540	cosbi(L,NW,NG) = atemp/btemp(L,nw)
	541	else
	542	cosbi(L,NW,NG) = 0.0D0
	543	end if
	544
	545
[716]	546	! Now the other Gauss points, if needed.
[135]	547
[716]	548	DO NG=1,L_NGAUSS-1
	549	IF(TAUGSURF(NW,NG) .gt. TLIMIT) THEN
[135]	550
[1648]	551	DO L=1,L_NLAYRAD-1
[716]	552	K = 2*L+1
	553	DTAUI(L,nw,ng) = DTAUKI(K,NW,NG)+DTAUKI(K+1,NW,NG)! + 1.e-50
	554
[961]	555	if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
[1722]	556
	557	WBARI(L,nw,ng) = btemp(L,nw) / DTAUI(L,NW,NG)
	558
[716]	559	else
	560	WBARI(L,nw,ng) = 0.0D0
[961]	561	DTAUI(L,NW,NG) = 1.0D-9
[716]	562	endif
	563
	564	cosbi(L,NW,NG) = cosbi(L,NW,L_NGAUSS)
	565	END DO ! L vertical loop
[1648]	566
[1722]	567	! Last level
	568	L = L_NLAYRAD
	569	K = 2*L+1
	570	DTAUI(L,nw,ng) = DTAUKI(K,NW,NG)! + 1.e-50
[1648]	571
[1722]	572	if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
	573
[1648]	574	WBARI(L,nw,ng) = btemp(L,nw) / DTAUI(L,NW,NG)
[1722]	575
	576	else
	577	WBARI(L,nw,ng) = 0.0D0
	578	DTAUI(L,NW,NG) = 1.0D-9
	579	endif
	580
	581	cosbi(L,NW,NG) = cosbi(L,NW,L_NGAUSS)
[1648]	582
[716]	583	END IF
	584
	585	END DO ! NG Gauss loop
	586	END DO ! NW spectral loop
	587
	588	! Total extinction optical depths
[3083]	589	!DO NG=1,L_NGAUSS ! full gauss loop
	590	! DO NW=1,L_NSPECTI
	591	! TAUCUMI(1,NW,NG)=0.0D0
	592	! DO K=2,L_LEVELS
	593	! TAUCUMI(K,NW,NG)=TAUCUMI(K-1,NW,NG)+DTAUKI(K,NW,NG)
	594	! END DO
	595	! END DO ! end full gauss loop
	596	!END DO
	597
	598	TAUCUMI(:,:,:) = DTAUKI(:,:,:)
[716]	599
	600	! be aware when comparing with textbook results
	601	! (e.g. Pierrehumbert p. 218) that
	602	! taucumi does not take the <cos theta>=0.5 factor into
	603	! account. It is the optical depth for a vertically
	604	! ascending ray with angle theta = 0.
	605
[1897]	606	if(firstcall) firstcall = .false.
	607
[716]	608	return
	609
	610
	611	end subroutine optci
	612
	613
	614

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

Download in other formats:

Original Format