Context Navigation

optci.F90 @ 3663

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

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

Download in other formats:

Original Format