Context Navigation

← Previous Revision
Next Revision →
Normal
Revision Log

optci.F90

Last change on this file was 3700, checked in by emoisan, 3 months ago

Titan:

Add safeguards not to enter microphisics when there are no clouds
Removed a comment

EMo

File size: 22.6 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
[3660]	175	! initialize k=1 level (to avoid NaN values and error in debug mode)
	176	DIAG_OPTT(1,:,1) = 0.0 ! dtau
	177	DIAG_OPTT(1,:,2) = 1.0 ! wbar
	178	DIAG_OPTT(1,:,3) = 1.0 ! gbar
	179	DIAG_OPTT(1,:,4) = 0.d0
	180	DIAG_OPTT(1,:,5) = 0.d0
	181
[3083]	182	do NW=1,L_NSPECTI
	183	! We ignore K=1...
	184	do K=2,L_LEVELS
	185	! int. arithmetic => gives the gcm layer index (reversed)
	186	ilay = L_NLAYRAD+1 - k/2
[135]	187
[3083]	188	! Optics coupled with the microphysics :
	189	IF (callmufi .AND. (.NOT. uncoupl_optic_haze)) THEN
[1722]	190
[3083]	191	!==========================================================================================
	192	! Old optics (must have callclouds = .false.):
	193	!==========================================================================================
	194	IF (.NOT. opt4clouds) THEN
	195	m3as = pqmo(ilay,2) / 2.0
	196	m3af = pqmo(ilay,4) / 2.0
[3340]	197	! Cut-off
	198	IF (ilay .lt. ilev_cutoff) THEN
	199	m3as = pqmo(ilev_cutoff,2) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(ilev_cutoff+1)-zlev(ilev_cutoff))
	200	m3af = pqmo(ilev_cutoff,4) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(ilev_cutoff+1)-zlev(ilev_cutoff))
[3083]	201	ENDIF
[1897]	202
[3083]	203	dtauaer_s = m3as*rhoaer_s(nw)
	204	dtauaer_f = m3af*rhoaer_f(nw)
	205
	206	!==========================================================================================
	207	! New optics :
	208	!==========================================================================================
	209	ELSE
	210	iw = (L_NSPECTI + 1) - NW + L_NSPECTV ! Visible first and return
	211	!-----------------------------
	212	! HAZE (Spherical + Fractal) :
	213	!-----------------------------
	214	FTYPE = 1
[1722]	215
[3083]	216	! Spherical aerosols :
	217	!---------------------
	218	CTYPE = 5
	219	m0as = pqmo(ilay,1) / 2.0
	220	m3as = pqmo(ilay,2) / 2.0
[3340]	221	! If not callclouds : must have a cut-off
[3083]	222	IF (.NOT. callclouds) THEN
[3340]	223	IF (ilay .lt. ilev_cutoff) THEN
	224	m0as = pqmo(ilev_cutoff,1) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(ilev_cutoff+1)-zlev(ilev_cutoff))
	225	m3as = pqmo(ilev_cutoff,2) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(ilev_cutoff+1)-zlev(ilev_cutoff))
[3083]	226	ENDIF
	227	ENDIF
	228	call get_haze_and_cloud_opacity(FTYPE,CTYPE,m0as,m3as,iw,dtauaer_s,ssa_s(nw),asf_s(nw))
	229
	230	! Fractal aerosols :
	231	!-------------------
	232	CTYPE = FTYPE
	233	m0af = pqmo(ilay,3) / 2.0
	234	m3af = pqmo(ilay,4) / 2.0
[3340]	235	! If not callclouds : must have a cut-off
[3083]	236	IF (.NOT. callclouds) THEN
[3340]	237	IF (ilay .lt. ilev_cutoff) THEN
	238	m0af = pqmo(ilev_cutoff,3) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(ilev_cutoff+1)-zlev(ilev_cutoff))
	239	m3af = pqmo(ilev_cutoff,4) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(ilev_cutoff+1)-zlev(ilev_cutoff))
[3083]	240	ENDIF
	241	ENDIF
	242	call get_haze_and_cloud_opacity(FTYPE,CTYPE,m0af,m3af,iw,dtauaer_f,ssa_f(nw),asf_f(nw))
	243	ENDIF
[135]	244
[3318]	245	! Tuning of optical properties for haze :
	246	!dtauaer_s = dtauaer_s * (FHIR * (1-ssa_s(nw)) + ssa_s(nw))
	247	!ssa_s(nw) = ssa_s(nw) / (FHIR * (1-ssa_s(nw)) + ssa_s(nw))
	248	dtauaer_f = dtauaer_f * (FHIR * (1-ssa_f(nw)) + ssa_f(nw))
	249	ssa_f(nw) = ssa_f(nw) / (FHIR * (1-ssa_f(nw)) + ssa_f(nw))
	250
[3083]	251	! Total of Haze opacity (dtau), SSA (w) and ASF (COS) :
	252	DHAZE_T(k,nw) = dtauaer_s + dtauaer_f
	253	IF (dtauaer_s + dtauaer_f .GT. 1.D-30) THEN
	254	SSA_T(k,nw) = ( dtauaer_sssa_s(nw) + dtauaer_fssa_f(nw) ) / ( dtauaer_s+dtauaer_f )
	255	ASF_T(k,nw) = ( dtauaer_sssa_s(nw)asf_s(nw) + dtauaer_fssa_f(nw)asf_f(nw) ) &
	256	/ ( ssa_s(nw)dtauaer_s + ssa_f(nw)dtauaer_f )
	257	ELSE
	258	DHAZE_T(k,nw) = 0.D0
	259	SSA_T(k,nw) = 1.0
	260	ASF_T(k,nw) = 1.0
	261	ENDIF
[3318]	262
[3340]	263	! Opacity and albedo adjustement below cutoff :
	264	IF (.NOT. callclouds) THEN
	265	corr_haze=0.6-0.4TANH((PMID(K)100.-2500.)/250.)
	266	IF (ilay .lt. ilev_cutoff) THEN
	267	DHAZE_T(k,nw) = DHAZE_T(k,nw) * corr_haze
	268	ENDIF
	269	ENDIF
	270
[3083]	271	! Diagnostics for the haze :
[3318]	272	DIAG_OPTH(k,nw,1) = DHAZE_T(k,nw) ! dtau
	273	DIAG_OPTH(k,nw,2) = SSA_T(k,nw) ! wbar
	274	DIAG_OPTH(k,nw,3) = ASF_T(k,nw) ! gbar
[253]	275
[3083]	276	!---------------------
	277	! CLOUDS (Spherical) :
	278	!---------------------
	279	IF (callclouds) THEN
	280	CTYPE = 0
	281	m0ccn = pqmo(ilay,5) / 2.0
	282	m3ccn = pqmo(ilay,6) / 2.0
[3318]	283	m3cld = pqmo(ilay,6) / 2.0
[3083]	284
	285	! Clear / Dark column method :
	286	!-----------------------------
[305]	287
[3083]	288	! CCN's SSA :
	289	call get_haze_and_cloud_opacity(FTYPE,FTYPE,m0ccn,m3ccn,iw,dtau_ccn,ssa_ccn(nw),asf_ccn(nw))
	290
[3497]	291	! Clear column (CCN, C2H2, C2H6, HCN, AC6H6) :
[3083]	292	IF (CDCOLUMN == 0) THEN
	293	DO iq = 2, 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
[3497]	298	! Dark column (CCN, CH4, C2H2, C2H6, HCN, AC6H6) :
[3083]	299	ELSEIF (CDCOLUMN == 1) THEN
	300	DO iq = 1, nice
	301	m3cld = m3cld + (pqmo(ilay,ices_indx(iq)) / 2.0)
	302	ENDDO
	303	call get_haze_and_cloud_opacity(FTYPE,CTYPE,m0ccn,m3cld,iw,dtau_cld,ssa_cld(nw),asf_cld(nw))
	304
	305	ELSE
	306	WRITE(,) 'WARNING OPTCI.F90 : CDCOLUMN MUST BE 0 OR 1'
	307	WRITE(,) 'WE USE DARK COLUMN ...'
	308	DO iq = 1, nice
	309	m3cld = m3cld + (pqmo(ilay,ices_indx(iq)) / 2.0)
	310	ENDDO
	311	call get_haze_and_cloud_opacity(FTYPE,CTYPE,m0ccn,m3cld,iw,dtau_cld,ssa_cld(nw),asf_cld(nw))
	312	ENDIF
	313
	314	! For small dropplets, opacity of nucleus dominates...
[3700]	315	IF ((m3ccn + m3cld) .le. tiny(m3ccn)) THEN !no cloud !emoisan tests
	316	dtau_cld = 0.
	317	ssa_cld(nw) = 0.
	318	ELSE
	319	dtau_cld = (dtau_cldm3ccn + dtau_cldm3cld) / (m3ccn + m3cld)
	320	ssa_cld(nw) = (ssa_ccn(nw)m3ccn + ssa_cld(nw)m3cld) / (m3ccn + m3cld)
	321	ENDIF
[3318]	322
	323	! Tuning of optical properties for clouds :
	324	dtau_cld = dtau_cld * (FCIR * (1-ssa_cld(nw)) + ssa_cld(nw))
	325	ssa_cld(nw) = ssa_cld(nw) / (FCIR * (1-ssa_cld(nw)) + ssa_cld(nw))
	326
[3083]	327	! Total of Haze + Clouds opacity (dtau), SSA (w) and ASF (COS) :
	328	DHAZE_T(k,nw) = dtauaer_s + dtauaer_f + dtau_cld
	329	IF (DHAZE_T(k,nw) .GT. 1.D-30) THEN
	330	SSA_T(k,nw) = ( dtauaer_sssa_s(nw) + dtauaer_fssa_f(nw) + dtau_cld*ssa_cld(nw) ) / ( dtauaer_s+dtauaer_f+dtau_cld )
	331	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) ) &
	332	/ ( ssa_s(nw)dtauaer_s + ssa_f(nw)dtauaer_f + ssa_cld(nw)*dtau_cld )
	333	ELSE
	334	DHAZE_T(k,nw) = 0.D0
	335	SSA_T(k,nw) = 1.0
	336	ASF_T(k,nw) = 1.0
	337	ENDIF
[253]	338
[3083]	339	! Diagnostics for clouds :
[3318]	340	DIAG_OPTT(k,nw,1) = DHAZE_T(k,nw) ! dtau
	341	DIAG_OPTT(k,nw,2) = SSA_T(k,nw) ! wbar
	342	DIAG_OPTT(k,nw,3) = ASF_T(k,nw) ! gbar
[135]	343
[3083]	344	ELSE
	345	! Diagnostics for clouds :
[3318]	346	DIAG_OPTT(k,nw,1) = 0.D0 ! dtau
	347	DIAG_OPTT(k,nw,2) = 1.0 ! wbar
	348	DIAG_OPTT(k,nw,3) = 1.0 ! gbar
[3083]	349	ENDIF
[3660]	350
[3083]	351	! Optics and microphysics no coupled :
	352	ELSE
	353	! Call fixed vertical haze profile of extinction - same for all columns
	354	call disr_haze(dz(k),plev(k),wnoi(nw),DHAZE_T(k,nw),SSA_T(k,nw),ASF_T(k,nw))
	355	if (seashaze) DHAZE_T(k,nw) = DHAZE_T(k,nw)*seashazefact(k)
	356	! Diagnostics for the haze :
[3318]	357	DIAG_OPTH(k,nw,1) = DHAZE_T(k,nw) ! dtau
	358	DIAG_OPTH(k,nw,2) = SSA_T(k,nw) ! wbar
	359	DIAG_OPTH(k,nw,3) = ASF_T(k,nw) ! gbar
[3083]	360	! Diagnostics for clouds :
[3318]	361	DIAG_OPTT(k,nw,1) = 0.D0 ! dtau
	362	DIAG_OPTT(k,nw,2) = 1.0 ! wbar
	363	DIAG_OPTT(k,nw,3) = 1.0 ! gbar
[3083]	364	ENDIF ! ENDIF callmufi
	365
	366	DCONT = 0.0d0 ! continuum absorption
[1648]	367
[3083]	368	if(continuum.and.(.not.graybody))then
	369	! include continua if necessary
	370	wn_cont = dble(wnoi(nw))
	371	T_cont = dble(TMID(k))
	372	do igas=1,ngasmx
[1648]	373
[3083]	374	p_cont = dble(PMID(k)scalepgfrac(igas,ilay))
[1648]	375
[3083]	376	dtemp=0.0d0
	377	if(igas.eq.igas_N2)then
[135]	378
[3083]	379	interm = indi(nw,igas,igas)
	380	call interpolateN2N2(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
	381	indi(nw,igas,igas) = interm
[135]	382
[3083]	383	elseif(igas.eq.igas_H2)then
[135]	384
[3083]	385	! first do self-induced absorption
	386	interm = indi(nw,igas,igas)
	387	call interpolateH2H2(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
	388	indi(nw,igas,igas) = interm
[135]	389
[3083]	390	! then cross-interactions with other gases
	391	do jgas=1,ngasmx
	392	p_cross = dble(PMID(k)scalepgfrac(jgas,ilay))
	393	dtempc = 0.0d0
	394	if(jgas.eq.igas_N2)then
	395	interm = indi(nw,igas,jgas)
	396	call interpolateN2H2(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
	397	indi(nw,igas,jgas) = interm
	398	endif
	399	dtemp = dtemp + dtempc
	400	enddo
[135]	401
[3083]	402	elseif(igas.eq.igas_CH4)then
[135]	403
[3083]	404	! first do self-induced absorption
	405	interm = indi(nw,igas,igas)
	406	call interpolateCH4CH4(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
	407	indi(nw,igas,igas) = interm
[135]	408
[3083]	409	! then cross-interactions with other gases
	410	do jgas=1,ngasmx
	411	p_cross = dble(PMID(k)scalepgfrac(jgas,ilay))
	412	dtempc = 0.0d0
	413	if(jgas.eq.igas_N2)then
	414	interm = indi(nw,igas,jgas)
	415	call interpolateN2CH4(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
	416	indi(nw,igas,jgas) = interm
	417	endif
	418	dtemp = dtemp + dtempc
	419	enddo
[253]	420
[3083]	421	endif
[1725]	422
[3083]	423	DCONT = DCONT + dtemp
[1725]	424
[3083]	425	enddo
[1725]	426
[3083]	427	DCONT = DCONT*dz(k)
[135]	428
[3083]	429	endif
[135]	430
[3083]	431	do ng=1,L_NGAUSS-1
[135]	432
[3083]	433	! Now compute TAUGAS
[135]	434
[3083]	435	! JVO 2017 : added tmpk because the repeated calls to gasi/v increased dramatically
	436	! the execution time of optci/v -> ~ factor 2 on the whole radiative
	437	! transfer on the tested simulations !
[135]	438
[3083]	439	if (corrk_recombin) then
	440	tmpk = GASI_RECOMB(MT(K):MT(K)+1,MP(K):MP(K)+1,NW,NG)
	441	else
	442	tmpk = GASI(MT(K):MT(K)+1,MP(K):MP(K)+1,1,NW,NG)
	443	endif
[135]	444
[3083]	445	KCOEF(1) = tmpk(1,1) ! KCOEF(1) = GASI(MT(K),MP(K),1,NW,NG)
	446	KCOEF(2) = tmpk(1,2) ! KCOEF(2) = GASI(MT(K),MP(K)+1,1,NW,NG)
	447	KCOEF(3) = tmpk(2,2) ! KCOEF(3) = GASI(MT(K)+1,MP(K)+1,1,NW,NG)
	448	KCOEF(4) = tmpk(2,1) ! KCOEF(4) = GASI(MT(K)+1,MP(K),1,NW,NG)
[135]	449
	450
[3083]	451	! Interpolate the gaseous k-coefficients to the requested T,P values
	452
	453	ANS = LKCOEF(K,1)KCOEF(1) + LKCOEF(K,2)KCOEF(2) + &
	454	LKCOEF(K,3)KCOEF(3) + LKCOEF(K,4)KCOEF(4)
	455
	456	TAUGAS = U(k)*ANS
	457
	458	TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS + DCONT
	459	DTAUKI(K,nw,ng) = TAUGAS &
	460	+ DCONT & ! For parameterized continuum absorption
	461	+ DHAZE_T(K,NW) ! For Titan haze
	462
	463	end do
	464
	465	! Now fill in the "clear" part of the spectrum (NG = L_NGAUSS),
	466	! which holds continuum opacity only
	467
	468	NG = L_NGAUSS
	469	DTAUKI(K,nw,ng) = 0.d0 &
	470	+ DCONT & ! For parameterized continuum absorption
	471	+ DHAZE_T(K,NW) ! For Titan Haze
	472
[3318]	473	DIAG_OPTH(K,nw,4) = 0.d0
	474	DIAG_OPTH(K,nw,5) = TAUGAS
	475	DIAG_OPTH(K,nw,6) = DCONT
	476	DIAG_OPTT(K,nw,4) = 0.d0
	477	DIAG_OPTT(K,nw,5) = TAUGAS
	478	DIAG_OPTT(K,nw,6) = DCONT
	479
[3083]	480	end do ! k = L_LEVELS
	481	end do ! nw = L_NSPECTI
	482
[716]	483	!=======================================================================
	484	! Now the full treatment for the layers, where besides the opacity
	485	! we need to calculate the scattering albedo and asymmetry factors
[1648]	486	! ======================================================================
[135]	487
[1648]	488	! Haze scattering
[918]	489	DO NW=1,L_NSPECTI
[1722]	490	DO K=2,L_LEVELS
[1648]	491	DHAZES_T(K,NW) = DHAZE_T(K,NW) * SSA_T(K,NW)
	492	ENDDO
	493	ENDDO
	494
	495	DO NW=1,L_NSPECTI
	496	DO L=1,L_NLAYRAD-1
[918]	497	K = 2*L+1
[1788]	498	btemp(L,NW) = DHAZES_T(K,NW) + DHAZES_T(K+1,NW)
[918]	499	END DO ! L vertical loop
[1648]	500
[1722]	501	! Last level
	502	L = L_NLAYRAD
	503	K = 2*L+1
[1788]	504	btemp(L,NW) = DHAZES_T(K,NW)
[1648]	505
[918]	506	END DO ! NW spectral loop
	507
[135]	508
[716]	509	DO NW=1,L_NSPECTI
	510	NG = L_NGAUSS
[1648]	511	DO L=1,L_NLAYRAD-1
[135]	512
[716]	513	K = 2*L+1
	514	DTAUI(L,nw,ng) = DTAUKI(K,NW,NG) + DTAUKI(K+1,NW,NG)! + 1.e-50
[135]	515
[716]	516	atemp = 0.
[961]	517	if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
[1722]	518	atemp = atemp + &
	519	ASF_T(K,NW)*DHAZES_T(K,NW) + &
	520	ASF_T(K+1,NW)*DHAZES_T(K+1,NW)
	521
[918]	522	WBARI(L,nw,ng) = btemp(L,nw) / DTAUI(L,NW,NG)
[716]	523	else
	524	WBARI(L,nw,ng) = 0.0D0
[961]	525	DTAUI(L,NW,NG) = 1.0D-9
[716]	526	endif
[135]	527
[961]	528	if(btemp(L,nw) .GT. 0.0d0) then
[918]	529	cosbi(L,NW,NG) = atemp/btemp(L,nw)
[716]	530	else
	531	cosbi(L,NW,NG) = 0.0D0
	532	end if
[135]	533
[716]	534	END DO ! L vertical loop
[1648]	535
[1722]	536	! Last level
[1648]	537
[1722]	538	L = L_NLAYRAD
	539	K = 2*L+1
	540	DTAUI(L,nw,ng) = DTAUKI(K,NW,NG) ! + 1.e-50
[1648]	541
[1722]	542	atemp = 0.
	543	if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
	544	atemp = atemp + ASF_T(K,NW)*DHAZES_T(K,NW)
	545	WBARI(L,nw,ng) = btemp(L,nw) / DTAUI(L,NW,NG)
	546	else
	547	WBARI(L,nw,ng) = 0.0D0
	548	DTAUI(L,NW,NG) = 1.0D-9
	549	endif
[1648]	550
[1722]	551	if(btemp(L,nw) .GT. 0.0d0) then
	552	cosbi(L,NW,NG) = atemp/btemp(L,nw)
	553	else
	554	cosbi(L,NW,NG) = 0.0D0
	555	end if
	556
	557
[716]	558	! Now the other Gauss points, if needed.
[135]	559
[716]	560	DO NG=1,L_NGAUSS-1
	561	IF(TAUGSURF(NW,NG) .gt. TLIMIT) THEN
[135]	562
[1648]	563	DO L=1,L_NLAYRAD-1
[716]	564	K = 2*L+1
	565	DTAUI(L,nw,ng) = DTAUKI(K,NW,NG)+DTAUKI(K+1,NW,NG)! + 1.e-50
	566
[961]	567	if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
[1722]	568
	569	WBARI(L,nw,ng) = btemp(L,nw) / DTAUI(L,NW,NG)
	570
[716]	571	else
	572	WBARI(L,nw,ng) = 0.0D0
[961]	573	DTAUI(L,NW,NG) = 1.0D-9
[716]	574	endif
	575
	576	cosbi(L,NW,NG) = cosbi(L,NW,L_NGAUSS)
	577	END DO ! L vertical loop
[1648]	578
[1722]	579	! Last level
	580	L = L_NLAYRAD
	581	K = 2*L+1
	582	DTAUI(L,nw,ng) = DTAUKI(K,NW,NG)! + 1.e-50
[1648]	583
[1722]	584	if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
	585
[1648]	586	WBARI(L,nw,ng) = btemp(L,nw) / DTAUI(L,NW,NG)
[1722]	587
	588	else
	589	WBARI(L,nw,ng) = 0.0D0
	590	DTAUI(L,NW,NG) = 1.0D-9
	591	endif
	592
	593	cosbi(L,NW,NG) = cosbi(L,NW,L_NGAUSS)
[1648]	594
[716]	595	END IF
	596
	597	END DO ! NG Gauss loop
	598	END DO ! NW spectral loop
	599
	600	! Total extinction optical depths
[3083]	601	!DO NG=1,L_NGAUSS ! full gauss loop
	602	! DO NW=1,L_NSPECTI
	603	! TAUCUMI(1,NW,NG)=0.0D0
	604	! DO K=2,L_LEVELS
	605	! TAUCUMI(K,NW,NG)=TAUCUMI(K-1,NW,NG)+DTAUKI(K,NW,NG)
	606	! END DO
	607	! END DO ! end full gauss loop
	608	!END DO
	609
	610	TAUCUMI(:,:,:) = DTAUKI(:,:,:)
[716]	611
	612	! be aware when comparing with textbook results
	613	! (e.g. Pierrehumbert p. 218) that
	614	! taucumi does not take the <cos theta>=0.5 factor into
	615	! account. It is the optical depth for a vertically
	616	! ascending ray with angle theta = 0.
	617
[1897]	618	if(firstcall) firstcall = .false.
	619
[716]	620	return
	621
	622
	623	end subroutine optci
	624
	625
	626

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

Download in other formats:

Original Format