Context Navigation

optci.F90 @ 671

Last change on this file since 671 was 596, checked in by jleconte, 13 years ago

Added double gray case (if graybody=true in callphys.def):
- opacities are set to a constant value in sugas_corrk.
- the values are kappa_IR m^{2/kg in the infrared (to be read in callphys.def)}

kappa_VI m^{2/kg in the visible (to be read in callphys.def)}

Cleaned continuum part in optc*
Added .def files for a typical 1d earth case in deftank (dry case for the moment)

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
[471]	7	use gases_h
[135]	8	implicit none
	9
	10	!==================================================================
	11	!
	12	! Purpose
	13	! -------
	14	! Calculates longwave optical constants at each level. For each
	15	! layer and spectral interval in the IR it calculates WBAR, DTAU
	16	! and COSBAR. For each level it calculates TAU.
	17	!
	18	! TAUI(L,LW) is the cumulative optical depth at level L (or alternatively
	19	! at the bottom of layer L), LW is the spectral wavelength interval.
	20	!
	21	! TLEV(L) - Temperature at the layer boundary (i.e., level)
	22	! PLEV(L) - Pressure at the layer boundary (i.e., level)
	23	!
	24	! Authors
	25	! -------
	26	! Adapted from the NASA Ames code by R. Wordsworth (2009)
	27	!
	28	!==================================================================
	29
	30
	31	#include "comcstfi.h"
	32	#include "callkeys.h"
	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
[596]	139	if(Continuum.and.(.not.graybody))then
[305]	140	! include continua if necessary
[596]	141	wn_cont = dble(wnoi(nw))
	142	T_cont = dble(TMID(k))
	143	do igas=1,ngasmx
[135]	144
[596]	145	if(gfrac(igas).eq.-1)then ! variable
	146	p_cont = dble(PMID(k)scalepQVAR(k)) ! qvar = mol/mol
	147	else
	148	p_cont = dble(PMID(k)scalepgfrac(igas)*(1.-QVAR(k)))
	149	endif
[135]	150
[596]	151	dtemp=0.0
	152	if(gnom(igas).eq.'H2_')then
	153	call interpolateH2H2(wn_cont,T_cont,p_cont,dtemp,.false.)
	154	elseif(gnom(igas).eq.'H2O'.and.T_cont.gt.200.0)then
	155	p_air = dble(PMID(k)*scalep) - p_cont ! note assumes air!!
	156	call interpolateH2Ocont(wn_cont,T_cont,p_cont,p_air,dtemp,.false.)
	157	endif
[135]	158
[596]	159	DCONT = DCONT + dtemp
	160	enddo
[253]	161
[596]	162	DCONT = DCONT*dz(k)
[526]	163	endif
[305]	164
[253]	165	!--- Kasting's CIA ----------------------------------------
	166	!DCO2 = dz(k)Ci(nw)(1.2859PMID(k)/1000.0)(TMID(k)/300.)**Ti(nw)
	167	!DCO2 = 130Ci(nw)(pmid(k)/1013.25)*2(tmid(k)/300.)*Ti(nw) dz(k)
	168	! these two have been verified to give the same results
	169	!----------------------------------------------------------
	170
	171
	172	do ng=1,L_NGAUSS-1
	173
[135]	174	! Now compute TAUGAS
	175
	176	! Interpolate between water mixing ratios
	177	! WRATIO = 0.0 if the requested water amount is equal to, or outside the
	178	! the water data range
	179
[253]	180	if(L_REFVAR.eq.1)then ! added by RW for special no variable case
[135]	181	KCOEF(1) = GASI(MT(K),MP(K),1,NW,NG)
	182	KCOEF(2) = GASI(MT(K),MP(K)+1,1,NW,NG)
	183	KCOEF(3) = GASI(MT(K)+1,MP(K)+1,1,NW,NG)
	184	KCOEF(4) = GASI(MT(K)+1,MP(K),1,NW,NG)
	185	else
	186
[253]	187	KCOEF(1) = GASI(MT(K),MP(K),NVAR(K),NW,NG) + WRATIO(K)* &
	188	(GASI(MT(K),MP(K),NVAR(K)+1,NW,NG) - &
	189	GASI(MT(K),MP(K),NVAR(K),NW,NG))
[135]	190
[253]	191	KCOEF(2) = GASI(MT(K),MP(K)+1,NVAR(K),NW,NG) + WRATIO(K)* &
	192	(GASI(MT(K),MP(K)+1,NVAR(K)+1,NW,NG) - &
	193	GASI(MT(K),MP(K)+1,NVAR(K),NW,NG))
[135]	194
[253]	195	KCOEF(3) = GASI(MT(K)+1,MP(K)+1,NVAR(K),NW,NG) + WRATIO(K)* &
	196	(GASI(MT(K)+1,MP(K)+1,NVAR(K)+1,NW,NG) - &
	197	GASI(MT(K)+1,MP(K)+1,NVAR(K),NW,NG))
[135]	198
[253]	199	KCOEF(4) = GASI(MT(K)+1,MP(K),NVAR(K),NW,NG) + WRATIO(K)* &
	200	(GASI(MT(K)+1,MP(K),NVAR(K)+1,NW,NG) - &
	201	GASI(MT(K)+1,MP(K),NVAR(K),NW,NG))
[135]	202	endif
	203
	204	! Interpolate the gaseous k-coefficients to the requested T,P values
	205
[253]	206	ANS = LKCOEF(K,1)KCOEF(1) + LKCOEF(K,2)KCOEF(2) + &
[135]	207	LKCOEF(K,3)KCOEF(3) + LKCOEF(K,4)KCOEF(4)
[253]	208
[135]	209	TAUGAS = U(k)*ANS
	210
[253]	211	TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS + DCONT
	212	!TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS
	213	DTAUKI(K,nw,ng) = TAUGAS + DCONT ! For parameterized continuum absorption
	214
[135]	215	do iaer=1,naerkind
[253]	216	DTAUKI(K,nw,ng) = DTAUKI(K,nw,ng) + TAEROS(K,NW,IAER)
	217	end do ! a bug was here!
[135]	218
	219	end do
	220
	221	! Now fill in the "clear" part of the spectrum (NG = L_NGAUSS),
	222	! which holds continuum opacity only
	223
	224	NG = L_NGAUSS
[253]	225	DTAUKI(K,nw,ng) = 0.0 + DCONT ! For parameterized continuum absorption
	226
[135]	227	do iaer=1,naerkind
[253]	228	DTAUKI(K,nw,ng) = DTAUKI(K,nw,ng) + TAEROS(K,NW,IAER)
	229	end do ! a bug was here!
[135]	230
	231	end do
	232	end do
	233
	234
	235	!=======================================================================
	236	! Now the full treatment for the layers, where besides the opacity
	237	! we need to calculate the scattering albedo and asymmetry factors
	238
	239	DO NW=1,L_NSPECTI
	240	DO K=2,L_LEVELS+1
	241	do iaer=1,naerkind
	242	TAUAEROLK(K,NW,IAER) = TAUAERO(K,IAER)*QSIAER(K,NW,IAER)
	243	end do
	244	ENDDO
	245	ENDDO
	246
	247	DO NW=1,L_NSPECTI
	248	NG = L_NGAUSS
	249	DO L=1,L_NLAYRAD
	250
	251	K = 2*L+1
	252	DTAUI(L,nw,ng) = DTAUKI(K,NW,NG) + DTAUKI(K+1,NW,NG)! + 1.e-50
	253
	254	atemp = 0.
	255	btemp = 0.
	256	if(DTAUI(L,NW,NG) .GT. 1.0E-9) then
	257	do iaer=1,naerkind
	258	atemp = atemp + &
	259	GIAER(K,NW,IAER) * TAUAEROLK(K,NW,IAER) + &
	260	GIAER(K+1,NW,IAER) * TAUAEROLK(K+1,NW,IAER)
	261	btemp = btemp + TAUAEROLK(K,NW,IAER) + TAUAEROLK(K+1,NW,IAER)
	262	! * + 1.e-10
	263	end do
	264	WBARI(L,nw,ng) = btemp / DTAUI(L,NW,NG)
	265	else
	266	WBARI(L,nw,ng) = 0.0D0
	267	DTAUI(L,NW,NG) = 1.0E-9
	268	endif
	269
	270	if(btemp .GT. 0.0) then
	271	cosbi(L,NW,NG) = atemp/btemp
	272	else
	273	cosbi(L,NW,NG) = 0.0D0
	274	end if
	275
	276	END DO ! L vertical loop
	277
	278	! Now the other Gauss points, if needed.
	279
	280	DO NG=1,L_NGAUSS-1
	281	IF(TAUGSURF(NW,NG) .gt. TLIMIT) THEN
	282
	283	DO L=1,L_NLAYRAD
	284	K = 2*L+1
	285	DTAUI(L,nw,ng) = DTAUKI(K,NW,NG)+DTAUKI(K+1,NW,NG)! + 1.e-50
	286
	287	btemp = 0.
	288	if(DTAUI(L,NW,NG) .GT. 1.0E-9) then
	289
	290	do iaer=1,naerkind
	291	btemp = btemp + TAUAEROLK(K,NW,IAER) + TAUAEROLK(K+1,NW,IAER)
	292	end do
	293	WBARI(L,nw,ng) = btemp / DTAUI(L,NW,NG)
	294
	295	else
	296	WBARI(L,nw,ng) = 0.0D0
	297	DTAUI(L,NW,NG) = 1.0E-9
	298	endif
	299
	300	cosbi(L,NW,NG) = cosbi(L,NW,L_NGAUSS)
	301	END DO ! L vertical loop
	302	END IF
	303
	304	END DO ! NG Gauss loop
	305	END DO ! NW spectral loop
	306
	307	! Total extinction optical depths
	308
	309	DO NW=1,L_NSPECTI
	310	DO NG=1,L_NGAUSS ! full gauss loop
	311	TAUI(1,NW,NG)=0.0D0
	312	DO L=1,L_NLAYRAD
	313	TAUI(L+1,NW,NG)=TAUI(L,NW,NG)+DTAUI(L,NW,NG)
	314	END DO
	315
	316	TAUCUMI(1,NW,NG)=0.0D0
	317	DO K=2,L_LEVELS
	318	TAUCUMI(K,NW,NG)=TAUCUMI(K-1,NW,NG)+DTAUKI(K,NW,NG)
	319	END DO
	320	END DO ! end full gauss loop
	321	END DO
	322
	323	return
	324
	325
[305]	326	end subroutine optci
[135]	327
	328
	329

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

Download in other formats:

Original Format