Context Navigation

optci.F90 @ 3073

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

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

Download in other formats:

Original Format