Context Navigation

optci.F90 @ 3580

Last change on this file since 3580 was 3233, checked in by gmilcareck, 11 months ago
Add the possibility to use a fixed vertical molar fraction profile for the collision-induced absorption like the correlated-k.
File size: 17.8 KB

Rev	Line
[2032]	1	MODULE optci_mod
	2
	3	IMPLICIT NONE
	4
	5	CONTAINS
	6
[716]	7	subroutine optci(PLEV,TLEV,DTAUI,TAUCUMI, &
	8	QXIAER,QSIAER,GIAER,COSBI,WBARI,TAUAERO, &
[3233]	9	TMID,PMID,TAUGSURF,QVAR,MUVAR,FRACVAR)
[135]	10
[2032]	11	use radinc_h, only: L_LEVELS, L_NLAYRAD, L_NSPECTI, L_NGAUSS, &
	12	L_NLEVRAD, L_REFVAR, naerkind
[2133]	13	use radcommon_h, only: gasi,tlimit,wrefVAR,Cmk,tgasref,pfgasref,wnoi,scalep,indi,glat_ig
[2875]	14	use gases_h, only: gfrac, ngasmx, igas_N2, igas_He, igas_H2O, igas_H2, &
	15	igas_CH4, igas_CO2
[1384]	16	use comcstfi_mod, only: g, r, mugaz
[3233]	17	use callkeys_mod, only: kastprof,continuum,graybody,varspec
[2543]	18	use recombin_corrk_mod, only: corrk_recombin, gasi_recomb
[2582]	19	use tpindex_mod, only: tpindex
	20
[716]	21	implicit none
[135]	22
[716]	23	!==================================================================
	24	!
	25	! Purpose
	26	! -------
	27	! Calculates longwave optical constants at each level. For each
	28	! layer and spectral interval in the IR it calculates WBAR, DTAU
	29	! and COSBAR. For each level it calculates TAU.
	30	!
[2133]	31	! TAUCUMI(L,LW) is the cumulative optical depth at level L (or alternatively
[716]	32	! at the bottom of layer L), LW is the spectral wavelength interval.
	33	!
	34	! TLEV(L) - Temperature at the layer boundary (i.e., level)
	35	! PLEV(L) - Pressure at the layer boundary (i.e., level)
	36	!
	37	! Authors
	38	! -------
	39	! Adapted from the NASA Ames code by R. Wordsworth (2009)
	40	!
	41	!==================================================================
[135]	42
	43
[2957]	44	real*8,intent(out) :: DTAUI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
[1715]	45	real*8 DTAUKI(L_LEVELS,L_NSPECTI,L_NGAUSS)
[716]	46	real*8 TAUI(L_NLEVRAD,L_NSPECTI,L_NGAUSS)
[2957]	47	real*8,intent(out) :: TAUCUMI(L_LEVELS,L_NSPECTI,L_NGAUSS)
	48	real*8,intent(in) :: PLEV(L_LEVELS)
	49	real*8,intent(in) :: TLEV(L_LEVELS) ! not used
	50	real*8,intent(in) :: TMID(L_LEVELS)
	51	real*8,intent(in) :: PMID(L_LEVELS)
	52	real*8,intent(out) :: COSBI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
	53	real*8,intent(out) :: WBARI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
[135]	54
[716]	55	! for aerosols
[2957]	56	real*8,intent(in) :: QXIAER(L_LEVELS,L_NSPECTI,NAERKIND)
	57	real*8,intent(in) :: QSIAER(L_LEVELS,L_NSPECTI,NAERKIND)
	58	real*8,intent(in) :: GIAER(L_LEVELS,L_NSPECTI,NAERKIND)
	59	real*8,intent(in) :: TAUAERO(L_LEVELS,NAERKIND)
[135]	60
[2972]	61	! local variables (saved for convenience as need be allocated)
	62	real*8,save,allocatable :: TAUAEROLK(:,:,:)
	63	real*8,save,allocatable :: TAEROS(:,:,:)
	64	!$OMP THREADPRIVATE(TAUAEROLK,TAEROS)
	65
[716]	66	integer L, NW, NG, K, LK, IAER
	67	integer MT(L_LEVELS), MP(L_LEVELS), NP(L_LEVELS)
	68	real*8 ANS, TAUGAS
	69	real*8 DPR(L_LEVELS), U(L_LEVELS)
	70	real*8 LCOEF(4), LKCOEF(L_LEVELS,4)
[135]	71
[2957]	72	real*8,intent(out) :: taugsurf(L_NSPECTI,L_NGAUSS-1)
[918]	73	real*8 DCONT,DAERO
[716]	74	double precision wn_cont, p_cont, p_air, T_cont, dtemp, dtempc
	75	double precision p_cross
[135]	76
[716]	77	! variable species mixing ratio variables
[2957]	78	real*8,intent(in) :: QVAR(L_LEVELS)
	79	real*8,intent(in) :: MUVAR(L_LEVELS)
[3233]	80	real*8,intent(in) :: FRACVAR(ngasmx,L_LEVELS)
[2957]	81	real*8 WRATIO(L_LEVELS)
[716]	82	real*8 KCOEF(4)
	83	integer NVAR(L_LEVELS)
[1725]	84
	85	! temporary variables to reduce memory access time to gasi
	86	real*8 tmpk(2,2)
	87	real*8 tmpkvar(2,2,2)
[135]	88
[716]	89	! temporary variables for multiple aerosol calculation
[918]	90	real*8 atemp
	91	real*8 btemp(L_NLAYRAD,L_NSPECTI)
[135]	92
[716]	93	! variables for k in units m^-1
[873]	94	real*8 dz(L_LEVELS)
	95	!real*8 rho !! see test below
[135]	96
[716]	97	integer igas, jgas
[253]	98
[873]	99	integer interm
[2972]	100
	101	logical :: firstcall=.true.
	102	!$OMP THREADPRIVATE(firstcall)
[873]	103
[716]	104	!--- Kasting's CIA ----------------------------------------
	105	!real*8, parameter :: Ci(L_NSPECTI)=[ &
	106	! 3.8E-5, 1.2E-5, 2.8E-6, 7.6E-7, 4.5E-7, 2.3E-7, &
	107	! 5.4E-7, 1.6E-6, 0.0, &
	108	! 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
	109	! 0.0, 4.0E-7, 4.0E-6, 1.4E-5, &
	110	! 1.0E-5, 1.2E-6, 2.0E-7, 5.0E-8, 3.0E-8, 0.0 ]
	111	!real*8, parameter :: Ti(L_NSPECTI)=[ -2.2, -1.9, &
	112	! -1.7, -1.7, -1.7, -1.7, -1.7, -1.7, &
	113	! 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
	114	! -1.7,-1.7,-1.7,-1.7,-1.7,-1.7,-1.7, -1.7,0.0 ]
	115	!----------------------------------------------------------
[253]	116
[2972]	117	if (firstcall) then
	118	! allocate local arrays of size "naerkind" (which are also
	119	! "saved" so that this is done only once in for all even if
	120	! we don't need to store the value from a time step to the next)
	121	allocate(TAUAEROLK(L_LEVELS,L_NSPECTI,NAERKIND))
	122	allocate(TAEROS(L_LEVELS,L_NSPECTI,NAERKIND))
	123	firstcall=.false.
	124	endif ! of if (firstcall)
	125
[873]	126	!! AS: to save time in computing continuum (see bilinearbig)
	127	IF (.not.ALLOCATED(indi)) THEN
[878]	128	ALLOCATE(indi(L_NSPECTI,ngasmx,ngasmx))
[873]	129	indi = -9999 ! this initial value means "to be calculated"
	130	ENDIF
	131
[716]	132	!=======================================================================
	133	! Determine the total gas opacity throughout the column, for each
	134	! spectral interval, NW, and each Gauss point, NG.
[135]	135
[716]	136	taugsurf(:,:) = 0.0
	137	dpr(:) = 0.0
	138	lkcoef(:,:) = 0.0
[135]	139
[716]	140	do K=2,L_LEVELS
	141	DPR(k) = PLEV(K)-PLEV(K-1)
[135]	142
[716]	143	!--- Kasting's CIA ----------------------------------------
	144	!dz(k)=dpr(k)189.02TMID(K)/(0.03720*PMID(K))
	145	! this is CO2 path length (in cm) as written by Francois
	146	! delta_z = delta_p * R_specific * T / (g * P)
	147	! But Kasting states that W is in units of _atmosphere_ cm
	148	! So we do
	149	!dz(k)=dz(k)*(PMID(K)/1013.25)
	150	!dz(k)=dz(k)/100.0 ! in m for SI calc
	151	!----------------------------------------------------------
[135]	152
[716]	153	! if we have continuum opacities, we need dz
	154	if(kastprof)then
[961]	155	dz(k) = dpr(k)(1000.0d08.3145d0/muvar(k))TMID(K)/(gPMID(K))
[1016]	156	U(k) = CmkDPR(k)mugaz/muvar(k)
[716]	157	else
[1194]	158	dz(k) = dpr(k)RTMID(K)/(glat_igPMID(K))mugaz/muvar(k)
[1016]	159	U(k) = CmkDPR(k)mugaz/muvar(k) ! only Cmk line in optci.F
	160	!JL13 the mugaz/muvar factor takes into account water meanmolecular weight if water is present
[716]	161	endif
[135]	162
[716]	163	call tpindex(PMID(K),TMID(K),QVAR(K),pfgasref,tgasref,WREFVAR, &
	164	LCOEF,MT(K),MP(K),NVAR(K),WRATIO(K))
[135]	165
[716]	166	do LK=1,4
	167	LKCOEF(K,LK) = LCOEF(LK)
	168	end do
[918]	169	end do ! levels
[253]	170
[1715]	171	! Spectral dependance of aerosol absorption
[918]	172	do iaer=1,naerkind
[716]	173	DO NW=1,L_NSPECTI
[918]	174	do K=2,L_LEVELS
[716]	175	TAEROS(K,NW,IAER) = TAUAERO(K,IAER) * QXIAER(K,NW,IAER)
[918]	176	end do ! levels
[716]	177	END DO
[918]	178	end do
[135]	179
[918]	180	do NW=1,L_NSPECTI
[135]	181
[918]	182	do K=2,L_LEVELS
[1715]	183
	184	DAERO=SUM(TAEROS(K,NW,1:naerkind)) ! aerosol absorption
[873]	185
[1715]	186	DCONT = 0.0d0 ! continuum absorption
[135]	187
[873]	188	if(continuum.and.(.not.graybody))then
[716]	189	! include continua if necessary
	190	wn_cont = dble(wnoi(nw))
	191	T_cont = dble(TMID(k))
	192	do igas=1,ngasmx
[135]	193
[716]	194	if(gfrac(igas).eq.-1)then ! variable
	195	p_cont = dble(PMID(k)scalepQVAR(k)) ! qvar = mol/mol
[3233]	196	elseif(varspec) then
	197	p_cont = dble(PMID(k)scalepFRACVAR(igas,k)*(1.-QVAR(k)))
[716]	198	else
	199	p_cont = dble(PMID(k)scalepgfrac(igas)*(1.-QVAR(k)))
	200	endif
[253]	201
[961]	202	dtemp=0.0d0
[716]	203	if(igas.eq.igas_N2)then
[305]	204
[878]	205	interm = indi(nw,igas,igas)
	206	call interpolateN2N2(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
	207	indi(nw,igas,igas) = interm
[253]	208
[716]	209	elseif(igas.eq.igas_H2)then
[253]	210
[716]	211	! first do self-induced absorption
[878]	212	interm = indi(nw,igas,igas)
[873]	213	call interpolateH2H2(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
[878]	214	indi(nw,igas,igas) = interm
[253]	215
[716]	216	! then cross-interactions with other gases
	217	do jgas=1,ngasmx
[3233]	218	if(varspec) then
	219	p_cross = dble(PMID(k)scalepFRACVAR(jgas,k)*(1.-QVAR(k)))
	220	else
	221	p_cross = dble(PMID(k)scalepgfrac(jgas)*(1.-QVAR(k)))
	222	endif
[961]	223	dtempc = 0.0d0
[716]	224	if(jgas.eq.igas_N2)then
[878]	225	interm = indi(nw,igas,jgas)
	226	call interpolateN2H2(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
	227	indi(nw,igas,jgas) = interm
[2860]	228	elseif(jgas.eq.igas_CO2)then
	229	interm = indi(nw,igas,jgas)
	230	call interpolateCO2H2(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
	231	indi(nw,igas,jgas) = interm
[716]	232	elseif(jgas.eq.igas_He)then
[878]	233	interm = indi(nw,igas,jgas)
[873]	234	call interpolateH2He(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
[878]	235	indi(nw,igas,jgas) = interm
[716]	236	endif
	237	dtemp = dtemp + dtempc
	238	enddo
[135]	239
[2655]	240	elseif(igas.eq.igas_CH4)then
	241
	242	! first do self-induced absorption
	243	interm = indi(nw,igas,igas)
	244	call interpolateCH4CH4(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
	245	indi(nw,igas,igas) = interm
	246
	247	! then cross-interactions with other gases
	248	do jgas=1,ngasmx
[3233]	249	if(varspec) then
	250	p_cross = dble(PMID(k)scalepFRACVAR(jgas,k)*(1.-QVAR(k)))
	251	else
	252	p_cross = dble(PMID(k)scalepgfrac(jgas)*(1.-QVAR(k)))
	253	endif
[2655]	254	dtempc = 0.0d0
	255	if(jgas.eq.igas_H2)then
	256	interm = indi(nw,igas,jgas)
	257	call interpolateH2CH4(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
	258	indi(nw,igas,jgas) = interm
[2861]	259	elseif(jgas.eq.igas_CO2)then
	260	interm = indi(nw,igas,jgas)
	261	call interpolateCO2CH4(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
	262	indi(nw,igas,jgas) = interm
[2655]	263	elseif(jgas.eq.igas_He)then
	264	interm = indi(nw,igas,jgas)
	265	call interpolateHeCH4(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
	266	indi(nw,igas,jgas) = interm
	267	endif
	268	dtemp = dtemp + dtempc
	269	enddo
	270
[2520]	271	elseif(igas.eq.igas_H2O.and.T_cont.gt.100.0)then
	272	! Compute self and foreign (with air) continuum of H2O
[716]	273	p_air = dble(PMID(k)*scalep) - p_cont ! note assumes background is air!
[2520]	274	interm = indi(nw,igas,igas)
	275	call interpolateH2O_self_foreign(wn_cont,T_cont,p_cont,p_air,dtemp,.false.,interm) ! MTCKD v3.3
	276	indi(nw,igas,igas) = interm
[135]	277
[716]	278	endif
[135]	279
[716]	280	DCONT = DCONT + dtemp
[135]	281
[716]	282	enddo
[135]	283
[716]	284	! Oobleck test
	285	!rho = PMID(k)scalep / (TMID(k)286.99)
	286	!if(WNOI(nw).gt.300.0 .and. WNOI(nw).lt.500.0)then
	287	! DCONT = rho * 0.125 * 4.6e-4
	288	!elseif(WNOI(nw).gt.500.0 .and. WNOI(nw).lt.700.0)then
	289	! DCONT = 1000dpr(k) 1.0 * 4.6e-4 / g
	290	! DCONT = rho * 1.0 * 4.6e-4
	291	!elseif(WNOI(nw).gt.700.0 .and. WNOI(nw).lt.900.0)then
	292	! DCONT = rho * 0.125 * 4.6e-4
	293	!endif
[135]	294
[716]	295	DCONT = DCONT*dz(k)
[135]	296
[716]	297	endif
[135]	298
[716]	299	do ng=1,L_NGAUSS-1
[135]	300
[716]	301	! Now compute TAUGAS
[135]	302
[716]	303	! Interpolate between water mixing ratios
	304	! WRATIO = 0.0 if the requested water amount is equal to, or outside the
	305	! the water data range
[253]	306
[716]	307	if(L_REFVAR.eq.1)then ! added by RW for special no variable case
[1725]	308
	309	! JVO 2017 : added tmpk because the repeated calls to gasi/v increased dramatically
	310	! the execution time of optci/v -> ~ factor 2 on the whole radiative
	311	! transfer on the tested simulations !
	312
[2543]	313	IF (corrk_recombin) THEN ! added by JVO
	314	tmpk = GASI_RECOMB(MT(K):MT(K)+1,MP(K):MP(K)+1,1,NW,NG) ! contains the mix of recombined species
	315	ELSE
	316	tmpk = GASI(MT(K):MT(K)+1,MP(K):MP(K)+1,1,NW,NG)
	317	ENDIF
[1725]	318
	319	KCOEF(1) = tmpk(1,1) ! KCOEF(1) = GASI(MT(K),MP(K),1,NW,NG)
	320	KCOEF(2) = tmpk(1,2) ! KCOEF(2) = GASI(MT(K),MP(K)+1,1,NW,NG)
	321	KCOEF(3) = tmpk(2,2) ! KCOEF(3) = GASI(MT(K)+1,MP(K)+1,1,NW,NG)
	322	KCOEF(4) = tmpk(2,1) ! KCOEF(4) = GASI(MT(K)+1,MP(K),1,NW,NG)
	323
[716]	324	else
[135]	325
[2543]	326	IF (corrk_recombin) THEN ! added by JVO
	327	tmpkvar = GASI_RECOMB(MT(K):MT(K)+1,MP(K):MP(K)+1,NVAR(K):NVAR(K)+1,NW,NG)
	328	ELSE
	329	tmpkvar = GASI(MT(K):MT(K)+1,MP(K):MP(K)+1,NVAR(K):NVAR(K)+1,NW,NG)
	330	ENDIF
[135]	331
[1725]	332	KCOEF(1) = tmpkvar(1,1,1) + WRATIO(K) * &
	333	( tmpkvar(1,1,2)-tmpkvar(1,1,1) )
[135]	334
[1725]	335	KCOEF(2) = tmpkvar(1,2,1) + WRATIO(K) * &
	336	( tmpkvar(1,2,2)-tmpkvar(1,2,1) )
[135]	337
[1725]	338	KCOEF(3) = tmpkvar(2,2,1) + WRATIO(K) * &
	339	( tmpkvar(2,2,2)-tmpkvar(2,2,1) )
	340
	341	KCOEF(4) = tmpkvar(2,1,1) + WRATIO(K) * &
	342	( tmpkvar(2,1,2)-tmpkvar(2,1,1) )
[873]	343
[716]	344	endif
[135]	345
[716]	346	! Interpolate the gaseous k-coefficients to the requested T,P values
[135]	347
[716]	348	ANS = LKCOEF(K,1)KCOEF(1) + LKCOEF(K,2)KCOEF(2) + &
	349	LKCOEF(K,3)KCOEF(3) + LKCOEF(K,4)KCOEF(4)
[135]	350
[716]	351	TAUGAS = U(k)*ANS
[135]	352
[716]	353	TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS + DCONT
[918]	354	DTAUKI(K,nw,ng) = TAUGAS &
	355	+ DCONT & ! For parameterized continuum absorption
	356	+ DAERO ! For aerosol absorption
[135]	357
[716]	358	end do
[135]	359
[716]	360	! Now fill in the "clear" part of the spectrum (NG = L_NGAUSS),
	361	! which holds continuum opacity only
[135]	362
[716]	363	NG = L_NGAUSS
[918]	364	DTAUKI(K,nw,ng) = 0.d0 &
	365	+ DCONT & ! For parameterized continuum absorption
	366	+ DAERO ! For aerosol absorption
[135]	367
[716]	368	end do
	369	end do
[135]	370
[716]	371	!=======================================================================
	372	! Now the full treatment for the layers, where besides the opacity
	373	! we need to calculate the scattering albedo and asymmetry factors
[135]	374
[873]	375	do iaer=1,naerkind
[918]	376	DO NW=1,L_NSPECTI
[1715]	377	DO K=2,L_LEVELS
	378	TAUAEROLK(K,NW,IAER) = TAUAERO(K,IAER)*QSIAER(K,NW,IAER) ! effect of scattering albedo
[716]	379	ENDDO
[918]	380	ENDDO
[873]	381	end do
[918]	382
	383	DO NW=1,L_NSPECTI
[1715]	384	DO L=1,L_NLAYRAD-1
[918]	385	K = 2*L+1
	386	btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind)) + SUM(TAUAEROLK(K+1,NW,1:naerkind))
	387	END DO ! L vertical loop
[1715]	388
	389	! Last level
	390	L = L_NLAYRAD
	391	K = 2*L+1
	392	btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind))
	393
[918]	394	END DO ! NW spectral loop
	395
[135]	396
[716]	397	DO NW=1,L_NSPECTI
	398	NG = L_NGAUSS
[1715]	399	DO L=1,L_NLAYRAD-1
[135]	400
[716]	401	K = 2*L+1
	402	DTAUI(L,nw,ng) = DTAUKI(K,NW,NG) + DTAUKI(K+1,NW,NG)! + 1.e-50
[135]	403
[716]	404	atemp = 0.
[961]	405	if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
[716]	406	do iaer=1,naerkind
	407	atemp = atemp + &
	408	GIAER(K,NW,IAER) * TAUAEROLK(K,NW,IAER) + &
	409	GIAER(K+1,NW,IAER) * TAUAEROLK(K+1,NW,IAER)
	410	end do
[918]	411	WBARI(L,nw,ng) = btemp(L,nw) / DTAUI(L,NW,NG)
[716]	412	else
	413	WBARI(L,nw,ng) = 0.0D0
[961]	414	DTAUI(L,NW,NG) = 1.0D-9
[716]	415	endif
[135]	416
[961]	417	if(btemp(L,nw) .GT. 0.0d0) then
[918]	418	cosbi(L,NW,NG) = atemp/btemp(L,nw)
[716]	419	else
	420	cosbi(L,NW,NG) = 0.0D0
	421	end if
[135]	422
[716]	423	END DO ! L vertical loop
[1715]	424
	425	! Last level
	426
	427	L = L_NLAYRAD
	428	K = 2*L+1
	429	DTAUI(L,nw,ng) = DTAUKI(K,NW,NG) ! + 1.e-50
[135]	430
[1715]	431	atemp = 0.
	432	if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
	433	do iaer=1,naerkind
	434	atemp = atemp + GIAER(K,NW,IAER) * TAUAEROLK(K,NW,IAER)
	435	end do
	436	WBARI(L,nw,ng) = btemp(L,nw) / DTAUI(L,NW,NG)
	437	else
	438	WBARI(L,nw,ng) = 0.0D0
	439	DTAUI(L,NW,NG) = 1.0D-9
	440	endif
	441
	442	if(btemp(L,nw) .GT. 0.0d0) then
	443	cosbi(L,NW,NG) = atemp/btemp(L,nw)
	444	else
	445	cosbi(L,NW,NG) = 0.0D0
	446	end if
	447
	448
[716]	449	! Now the other Gauss points, if needed.
[135]	450
[716]	451	DO NG=1,L_NGAUSS-1
	452	IF(TAUGSURF(NW,NG) .gt. TLIMIT) THEN
[135]	453
[1715]	454	DO L=1,L_NLAYRAD-1
[716]	455	K = 2*L+1
	456	DTAUI(L,nw,ng) = DTAUKI(K,NW,NG)+DTAUKI(K+1,NW,NG)! + 1.e-50
	457
[961]	458	if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
[716]	459
[918]	460	WBARI(L,nw,ng) = btemp(L,nw) / DTAUI(L,NW,NG)
[716]	461
	462	else
	463	WBARI(L,nw,ng) = 0.0D0
[961]	464	DTAUI(L,NW,NG) = 1.0D-9
[716]	465	endif
	466
	467	cosbi(L,NW,NG) = cosbi(L,NW,L_NGAUSS)
	468	END DO ! L vertical loop
[1715]	469
	470	! Last level
	471	L = L_NLAYRAD
	472	K = 2*L+1
	473	DTAUI(L,nw,ng) = DTAUKI(K,NW,NG)! + 1.e-50
	474
	475	if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
	476
	477	WBARI(L,nw,ng) = btemp(L,nw) / DTAUI(L,NW,NG)
	478
	479	else
	480	WBARI(L,nw,ng) = 0.0D0
	481	DTAUI(L,NW,NG) = 1.0D-9
	482	endif
	483
	484	cosbi(L,NW,NG) = cosbi(L,NW,L_NGAUSS)
	485
[716]	486	END IF
	487
	488	END DO ! NG Gauss loop
	489	END DO ! NW spectral loop
	490
	491	! Total extinction optical depths
	492
[918]	493	DO NG=1,L_NGAUSS ! full gauss loop
	494	DO NW=1,L_NSPECTI
[716]	495	TAUCUMI(1,NW,NG)=0.0D0
	496	DO K=2,L_LEVELS
	497	TAUCUMI(K,NW,NG)=TAUCUMI(K-1,NW,NG)+DTAUKI(K,NW,NG)
	498	END DO
	499	END DO ! end full gauss loop
	500	END DO
	501
	502	! be aware when comparing with textbook results
	503	! (e.g. Pierrehumbert p. 218) that
	504	! taucumi does not take the <cos theta>=0.5 factor into
	505	! account. It is the optical depth for a vertically
	506	! ascending ray with angle theta = 0.
	507
	508	!open(127,file='taucum.out')
	509	!do nw=1,L_NSPECTI
	510	! write(127,*) taucumi(L_LEVELS,nw,L_NGAUSS)
	511	!enddo
	512	!close(127)
[918]	513
	514	! print*,'WBARI'
	515	! print*,WBARI
	516	! print*,'DTAUI'
	517	! print*,DTAUI
	518	! call abort
[2131]	519
[716]	520	end subroutine optci
	521
[2032]	522	END MODULE optci_mod
[716]	523

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

Download in other formats:

Original Format