Context Navigation

optci.F90 @ 374

Last change on this file since 374 was 305, checked in by rwordsworth, 13 years ago
Several new files added as part of the climate evolution model (main program kcm.F90). Some general cleanup in physiq.F90 and callcorrk.F90. Bugs in dust radiative transfer and H2 Rayleigh scattering corrected.
File size: 11.4 KB

Rev	Line
[135]	1	subroutine optci(PLEV,TLEV,DTAUI,TAUCUMI, &
	2	QXIAER,QSIAER,GIAER,COSBI,WBARI,TAUAERO, &
[305]	3	TMID,PMID,TAUGSURF,QVAR,MUVAR)
[135]	4
	5	use radinc_h
[253]	6	use radcommon_h, only: gasi, tlimit, wrefVAR, Cmk,tgasref,pfgasref,wnoi,scalep
[135]	7	implicit none
	8
	9	!==================================================================
	10	!
	11	! Purpose
	12	! -------
	13	! Calculates longwave optical constants at each level. For each
	14	! layer and spectral interval in the IR it calculates WBAR, DTAU
	15	! and COSBAR. For each level it calculates TAU.
	16	!
	17	! TAUI(L,LW) is the cumulative optical depth at level L (or alternatively
	18	! at the bottom of layer L), LW is the spectral wavelength interval.
	19	!
	20	! TLEV(L) - Temperature at the layer boundary (i.e., level)
	21	! PLEV(L) - Pressure at the layer boundary (i.e., level)
	22	!
	23	! Authors
	24	! -------
	25	! Adapted from the NASA Ames code by R. Wordsworth (2009)
	26	!
	27	!==================================================================
	28
	29
	30	#include "comcstfi.h"
	31	#include "callkeys.h"
[253]	32	#include "gases.h"
[135]	33
	34
	35	real*8 DTAUI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
	36	real*8 DTAUKI(L_LEVELS+1,L_NSPECTI,L_NGAUSS)
	37	real*8 TAUI(L_NLEVRAD,L_NSPECTI,L_NGAUSS)
	38	real*8 TAUCUMI(L_LEVELS,L_NSPECTI,L_NGAUSS)
	39	real*8 PLEV(L_LEVELS)
	40	real*8 TLEV(L_LEVELS)
	41	real*8 TMID(L_LEVELS), PMID(L_LEVELS)
	42	real*8 COSBI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
	43	real*8 WBARI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
	44
	45	! For aerosols
	46	real*8 QXIAER(L_LEVELS+1,L_NSPECTI,NAERKIND)
	47	real*8 QSIAER(L_LEVELS+1,L_NSPECTI,NAERKIND)
	48	real*8 GIAER(L_LEVELS+1,L_NSPECTI,NAERKIND)
	49	real*8 TAUAERO(L_LEVELS+1,NAERKIND)
	50	real*8 TAUAEROLK(L_LEVELS+1,L_NSPECTI,NAERKIND)
	51	real*8 TAEROS(L_LEVELS,L_NSPECTI,NAERKIND)
	52
	53	integer L, NW, NG, K, LK, IAER
	54	integer MT(L_LEVELS), MP(L_LEVELS), NP(L_LEVELS)
	55	real*8 ANS, TAUGAS
	56	real*8 DPR(L_LEVELS), U(L_LEVELS)
	57	real*8 LCOEF(4), LKCOEF(L_LEVELS,4)
	58
	59	real*8 taugsurf(L_NSPECTI,L_NGAUSS-1)
[253]	60	real*8 DCONT
[305]	61	double precision wn_cont, p_cont, p_air, T_cont, dtemp
[135]	62
[305]	63	! Variable species mixing ratio variables
	64	real*8 QVAR(L_LEVELS), WRATIO(L_LEVELS), MUVAR(L_LEVELS)
[135]	65	real*8 KCOEF(4)
	66	integer NVAR(L_LEVELS)
	67
	68	! temporary variables for multiple aerosol calculation
	69	real*8 atemp, btemp
	70
	71	! variables for k in units m^-1
[253]	72	real*8 rho, dz(L_LEVELS)
[135]	73
[253]	74	integer igas
	75
	76	!--- Kasting's CIA ----------------------------------------
	77	!real*8, parameter :: Ci(L_NSPECTI)=[ &
	78	! 3.8E-5, 1.2E-5, 2.8E-6, 7.6E-7, 4.5E-7, 2.3E-7, &
	79	! 5.4E-7, 1.6E-6, 0.0, &
	80	! 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
	81	! 0.0, 4.0E-7, 4.0E-6, 1.4E-5, &
	82	! 1.0E-5, 1.2E-6, 2.0E-7, 5.0E-8, 3.0E-8, 0.0 ]
	83	!real*8, parameter :: Ti(L_NSPECTI)=[ -2.2, -1.9, &
	84	! -1.7, -1.7, -1.7, -1.7, -1.7, -1.7, &
	85	! 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
	86	! -1.7,-1.7,-1.7,-1.7,-1.7,-1.7,-1.7, -1.7,0.0 ]
	87	!----------------------------------------------------------
	88
[135]	89	!=======================================================================
	90	! Determine the total gas opacity throughout the column, for each
	91	! spectral interval, NW, and each Gauss point, NG.
	92
[253]	93	taugsurf(:,:) = 0.0
	94	dpr(:) = 0.0
	95	lkcoef(:,:) = 0.0
[135]	96
	97	do K=2,L_LEVELS
	98	DPR(k) = PLEV(K)-PLEV(K-1)
	99
[253]	100	!--- Kasting's CIA ----------------------------------------
	101	!dz(k)=dpr(k)189.02TMID(K)/(0.03720*PMID(K))
	102	! this is CO2 path length (in cm) as written by Francois
	103	! delta_z = delta_p * R_specific * T / (g * P)
	104	! But Kasting states that W is in units of _atmosphere_ cm
	105	! So we do
	106	!dz(k)=dz(k)*(PMID(K)/1013.25)
	107	!dz(k)=dz(k)/100.0 ! in m for SI calc
	108	!----------------------------------------------------------
[135]	109
[253]	110	! if we have continuum opacities, we need dz
[305]	111	if(kastprof)then
	112	dz(k) = dpr(k)(8314.5/muvar(k))TMID(K)/(g*PMID(K))
	113	U(k) = (Cmkmugaz/(muvar(k)))DPR(k)
	114	else
	115	dz(k) = dpr(k)RTMID(K)/(g*PMID(K))
	116	U(k) = Cmk*DPR(k) ! only Cmk line in optci.F
	117	endif
[135]	118
	119	call tpindex(PMID(K),TMID(K),QVAR(K),pfgasref,tgasref,WREFVAR, &
	120	LCOEF,MT(K),MP(K),NVAR(K),WRATIO(K))
[253]	121
[135]	122	do LK=1,4
	123	LKCOEF(K,LK) = LCOEF(LK)
	124	end do
	125
[253]	126
[135]	127	DO NW=1,L_NSPECTI
	128	do iaer=1,naerkind
[253]	129	TAEROS(K,NW,IAER) = TAUAERO(K,IAER) * QXIAER(K,NW,IAER)
[135]	130	end do
	131	END DO
	132	end do ! levels
	133
	134	do K=2,L_LEVELS
	135	do nw=1,L_NSPECTI
	136
[253]	137	DCONT = 0.0 ! continuum absorption
[135]	138
[305]	139	! include continua if necessary
	140	wn_cont = dble(wnoi(nw))
	141	T_cont = dble(TMID(k))
[253]	142	do igas=1,ngasmx
[135]	143
[305]	144	if(gfrac(igas).eq.-1)then ! variable
	145	p_cont = dble(PMID(k)scalepQVAR(k)) ! qvar = mol/mol
	146	else
	147	p_cont = dble(PMID(k)scalepgfrac(igas)*(1.-QVAR(k)))
	148	endif
[135]	149
[305]	150	dtemp=0.0
	151	if(gnom(igas).eq.'H2_')then
	152	call interpolateH2H2(wn_cont,T_cont,p_cont,dtemp,.false.)
	153	elseif(gnom(igas).eq.'H2O'.and.T_cont.gt.200.0)then
	154	p_air = dble(PMID(k)*scalep) - p_cont ! note assumes air!!
	155	call interpolateH2Ocont(wn_cont,T_cont,p_cont,p_air,dtemp,.false.)
[135]	156
[305]	157	endif
[253]	158
[305]	159	DCONT = DCONT + dtemp
[253]	160
	161	enddo
	162
[305]	163
	164	DCONT = DCONT*dz(k)
	165
[253]	166	!--- Kasting's CIA ----------------------------------------
	167	!DCO2 = dz(k)Ci(nw)(1.2859PMID(k)/1000.0)(TMID(k)/300.)**Ti(nw)
	168	!DCO2 = 130Ci(nw)(pmid(k)/1013.25)*2(tmid(k)/300.)*Ti(nw) dz(k)
	169	! these two have been verified to give the same results
	170	!----------------------------------------------------------
	171
	172
	173	do ng=1,L_NGAUSS-1
	174
[135]	175	! Now compute TAUGAS
	176
	177	! Interpolate between water mixing ratios
	178	! WRATIO = 0.0 if the requested water amount is equal to, or outside the
	179	! the water data range
	180
[253]	181	if(L_REFVAR.eq.1)then ! added by RW for special no variable case
[135]	182	KCOEF(1) = GASI(MT(K),MP(K),1,NW,NG)
	183	KCOEF(2) = GASI(MT(K),MP(K)+1,1,NW,NG)
	184	KCOEF(3) = GASI(MT(K)+1,MP(K)+1,1,NW,NG)
	185	KCOEF(4) = GASI(MT(K)+1,MP(K),1,NW,NG)
	186	else
	187
[253]	188	KCOEF(1) = GASI(MT(K),MP(K),NVAR(K),NW,NG) + WRATIO(K)* &
	189	(GASI(MT(K),MP(K),NVAR(K)+1,NW,NG) - &
	190	GASI(MT(K),MP(K),NVAR(K),NW,NG))
[135]	191
[253]	192	KCOEF(2) = GASI(MT(K),MP(K)+1,NVAR(K),NW,NG) + WRATIO(K)* &
	193	(GASI(MT(K),MP(K)+1,NVAR(K)+1,NW,NG) - &
	194	GASI(MT(K),MP(K)+1,NVAR(K),NW,NG))
[135]	195
[253]	196	KCOEF(3) = GASI(MT(K)+1,MP(K)+1,NVAR(K),NW,NG) + WRATIO(K)* &
	197	(GASI(MT(K)+1,MP(K)+1,NVAR(K)+1,NW,NG) - &
	198	GASI(MT(K)+1,MP(K)+1,NVAR(K),NW,NG))
[135]	199
[253]	200	KCOEF(4) = GASI(MT(K)+1,MP(K),NVAR(K),NW,NG) + WRATIO(K)* &
	201	(GASI(MT(K)+1,MP(K),NVAR(K)+1,NW,NG) - &
	202	GASI(MT(K)+1,MP(K),NVAR(K),NW,NG))
[135]	203	endif
	204
	205	! Interpolate the gaseous k-coefficients to the requested T,P values
	206
[253]	207	ANS = LKCOEF(K,1)KCOEF(1) + LKCOEF(K,2)KCOEF(2) + &
[135]	208	LKCOEF(K,3)KCOEF(3) + LKCOEF(K,4)KCOEF(4)
[253]	209
[135]	210	TAUGAS = U(k)*ANS
	211
[305]	212	if(graybody)then
	213	TAUGAS = 0.0
	214	DCONT = U(k)*3.3e-26
	215	endif
	216
[253]	217	TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS + DCONT
	218	!TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS
[135]	219
[253]	220	DTAUKI(K,nw,ng) = TAUGAS + DCONT ! For parameterized continuum absorption
	221
[135]	222	do iaer=1,naerkind
[253]	223	DTAUKI(K,nw,ng) = DTAUKI(K,nw,ng) + TAEROS(K,NW,IAER)
	224	end do ! a bug was here!
[135]	225
	226	end do
	227
	228	! Now fill in the "clear" part of the spectrum (NG = L_NGAUSS),
	229	! which holds continuum opacity only
	230
	231	NG = L_NGAUSS
[253]	232	DTAUKI(K,nw,ng) = 0.0 + DCONT ! For parameterized continuum absorption
	233
[135]	234	do iaer=1,naerkind
[253]	235	DTAUKI(K,nw,ng) = DTAUKI(K,nw,ng) + TAEROS(K,NW,IAER)
	236	end do ! a bug was here!
[135]	237
	238	end do
	239	end do
	240
	241
	242	!=======================================================================
	243	! Now the full treatment for the layers, where besides the opacity
	244	! we need to calculate the scattering albedo and asymmetry factors
	245
	246	DO NW=1,L_NSPECTI
	247	DO K=2,L_LEVELS+1
	248	do iaer=1,naerkind
	249	TAUAEROLK(K,NW,IAER) = TAUAERO(K,IAER)*QSIAER(K,NW,IAER)
	250	end do
	251	ENDDO
	252	ENDDO
	253
	254	DO NW=1,L_NSPECTI
	255	NG = L_NGAUSS
	256	DO L=1,L_NLAYRAD
	257
	258	K = 2*L+1
	259	DTAUI(L,nw,ng) = DTAUKI(K,NW,NG) + DTAUKI(K+1,NW,NG)! + 1.e-50
	260
	261	atemp = 0.
	262	btemp = 0.
	263	if(DTAUI(L,NW,NG) .GT. 1.0E-9) then
	264	do iaer=1,naerkind
	265	atemp = atemp + &
	266	GIAER(K,NW,IAER) * TAUAEROLK(K,NW,IAER) + &
	267	GIAER(K+1,NW,IAER) * TAUAEROLK(K+1,NW,IAER)
	268	btemp = btemp + TAUAEROLK(K,NW,IAER) + TAUAEROLK(K+1,NW,IAER)
	269	! * + 1.e-10
	270	end do
	271	WBARI(L,nw,ng) = btemp / DTAUI(L,NW,NG)
	272	else
	273	WBARI(L,nw,ng) = 0.0D0
	274	DTAUI(L,NW,NG) = 1.0E-9
	275	endif
	276
	277	if(btemp .GT. 0.0) then
	278	cosbi(L,NW,NG) = atemp/btemp
	279	else
	280	cosbi(L,NW,NG) = 0.0D0
	281	end if
	282
	283	END DO ! L vertical loop
	284
	285	! Now the other Gauss points, if needed.
	286
	287	DO NG=1,L_NGAUSS-1
	288	IF(TAUGSURF(NW,NG) .gt. TLIMIT) THEN
	289
	290	DO L=1,L_NLAYRAD
	291	K = 2*L+1
	292	DTAUI(L,nw,ng) = DTAUKI(K,NW,NG)+DTAUKI(K+1,NW,NG)! + 1.e-50
	293
	294	btemp = 0.
	295	if(DTAUI(L,NW,NG) .GT. 1.0E-9) then
	296
	297	do iaer=1,naerkind
	298	btemp = btemp + TAUAEROLK(K,NW,IAER) + TAUAEROLK(K+1,NW,IAER)
	299	end do
	300	WBARI(L,nw,ng) = btemp / DTAUI(L,NW,NG)
	301
	302	else
	303	WBARI(L,nw,ng) = 0.0D0
	304	DTAUI(L,NW,NG) = 1.0E-9
	305	endif
	306
	307	cosbi(L,NW,NG) = cosbi(L,NW,L_NGAUSS)
	308	END DO ! L vertical loop
	309	END IF
	310
	311	END DO ! NG Gauss loop
	312	END DO ! NW spectral loop
	313
	314	! Total extinction optical depths
	315
	316	DO NW=1,L_NSPECTI
	317	DO NG=1,L_NGAUSS ! full gauss loop
	318	TAUI(1,NW,NG)=0.0D0
	319	DO L=1,L_NLAYRAD
	320	TAUI(L+1,NW,NG)=TAUI(L,NW,NG)+DTAUI(L,NW,NG)
	321	END DO
	322
	323	TAUCUMI(1,NW,NG)=0.0D0
	324	DO K=2,L_LEVELS
	325	TAUCUMI(K,NW,NG)=TAUCUMI(K-1,NW,NG)+DTAUKI(K,NW,NG)
	326	END DO
	327	END DO ! end full gauss loop
	328	END DO
	329
	330	return
	331
	332
[305]	333	end subroutine optci
[135]	334
	335
	336

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

Download in other formats:

Original Format