Context Navigation

optci.F90 @ 258

Last change on this file since 258 was 253, checked in by emillour, 14 years ago

Generic GCM

Massive update to version 0.7

EM+RW

File size: 11.3 KB

Rev	Line
[135]	1	subroutine optci(PLEV,TLEV,DTAUI,TAUCUMI, &
	2	QXIAER,QSIAER,GIAER,COSBI,WBARI,TAUAERO, &
	3	TMID,PMID,TAUGSURF,QVAR)
	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
	61	double precision wn_cont, p_cont, T_cont
[135]	62
	63	! Water mixing ratio variables
	64	real*8 QVAR(L_LEVELS), WRATIO(L_LEVELS)
	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)
[253]	99	U(k) = Cmk*DPR(k) ! only Cmk line in optci.F
[135]	100
[253]	101	!--- Kasting's CIA ----------------------------------------
	102	!dz(k)=dpr(k)189.02TMID(K)/(0.03720*PMID(K))
	103	! this is CO2 path length (in cm) as written by Francois
	104	! delta_z = delta_p * R_specific * T / (g * P)
	105	! But Kasting states that W is in units of _atmosphere_ cm
	106	! So we do
	107	!dz(k)=dz(k)*(PMID(K)/1013.25)
	108	!dz(k)=dz(k)/100.0 ! in m for SI calc
	109	!----------------------------------------------------------
[135]	110
[253]	111	! if we have continuum opacities, we need dz
	112	dz(k)=dpr(k)RTMID(K)/(g*PMID(K))
[135]	113
	114	call tpindex(PMID(K),TMID(K),QVAR(K),pfgasref,tgasref,WREFVAR, &
	115	LCOEF,MT(K),MP(K),NVAR(K),WRATIO(K))
[253]	116
[135]	117	do LK=1,4
	118	LKCOEF(K,LK) = LCOEF(LK)
	119	end do
	120
[253]	121
[135]	122	DO NW=1,L_NSPECTI
	123	do iaer=1,naerkind
[253]	124	TAEROS(K,NW,IAER) = TAUAERO(K,IAER) * QXIAER(K,NW,IAER)
[135]	125	end do
	126	END DO
	127	end do ! levels
	128
	129	do K=2,L_LEVELS
	130	do nw=1,L_NSPECTI
	131
[253]	132	DCONT = 0.0 ! continuum absorption
[135]	133
[253]	134	! include H2 continuum
	135	do igas=1,ngasmx
	136	if(gnom(igas).eq.'H2_')then
[135]	137
[253]	138	wn_cont = dble(wnoi(nw))
	139	T_cont = dble(TMID(k))
	140	p_cont = dble(PMID(k)scalepgfrac(igas))
[135]	141
[253]	142	call interpolateH2H2(wn_cont,T_cont,p_cont,DCONT,.false.)
[135]	143
[253]	144	DCONT=DCONT*dz(k)
	145
	146
	147	endif
	148	enddo
	149
	150	!--- Kasting's CIA ----------------------------------------
	151	!DCO2 = dz(k)Ci(nw)(1.2859PMID(k)/1000.0)(TMID(k)/300.)**Ti(nw)
	152	!DCO2 = 130Ci(nw)(pmid(k)/1013.25)*2(tmid(k)/300.)*Ti(nw) dz(k)
	153	! these two have been verified to give the same results
	154	!----------------------------------------------------------
	155
	156	! Water continuum currently inactive!
	157	!if(varactive)then
	158	! call water_cont(PMID(K),TMID(K),WRATIO(K),WNOI(nw),DCO2)
	159	!endif
	160
	161	do ng=1,L_NGAUSS-1
	162
[135]	163	! Now compute TAUGAS
	164
	165	! Interpolate between water mixing ratios
	166	! WRATIO = 0.0 if the requested water amount is equal to, or outside the
	167	! the water data range
	168
[253]	169	if(L_REFVAR.eq.1)then ! added by RW for special no variable case
[135]	170	KCOEF(1) = GASI(MT(K),MP(K),1,NW,NG)
	171	KCOEF(2) = GASI(MT(K),MP(K)+1,1,NW,NG)
	172	KCOEF(3) = GASI(MT(K)+1,MP(K)+1,1,NW,NG)
	173	KCOEF(4) = GASI(MT(K)+1,MP(K),1,NW,NG)
	174	else
	175
[253]	176	KCOEF(1) = GASI(MT(K),MP(K),NVAR(K),NW,NG) + WRATIO(K)* &
	177	(GASI(MT(K),MP(K),NVAR(K)+1,NW,NG) - &
	178	GASI(MT(K),MP(K),NVAR(K),NW,NG))
[135]	179
[253]	180	KCOEF(2) = GASI(MT(K),MP(K)+1,NVAR(K),NW,NG) + WRATIO(K)* &
	181	(GASI(MT(K),MP(K)+1,NVAR(K)+1,NW,NG) - &
	182	GASI(MT(K),MP(K)+1,NVAR(K),NW,NG))
[135]	183
[253]	184	KCOEF(3) = GASI(MT(K)+1,MP(K)+1,NVAR(K),NW,NG) + WRATIO(K)* &
	185	(GASI(MT(K)+1,MP(K)+1,NVAR(K)+1,NW,NG) - &
	186	GASI(MT(K)+1,MP(K)+1,NVAR(K),NW,NG))
[135]	187
[253]	188	KCOEF(4) = GASI(MT(K)+1,MP(K),NVAR(K),NW,NG) + WRATIO(K)* &
	189	(GASI(MT(K)+1,MP(K),NVAR(K)+1,NW,NG) - &
	190	GASI(MT(K)+1,MP(K),NVAR(K),NW,NG))
[135]	191	endif
	192
	193	! Interpolate the gaseous k-coefficients to the requested T,P values
	194
[253]	195	ANS = LKCOEF(K,1)KCOEF(1) + LKCOEF(K,2)KCOEF(2) + &
[135]	196	LKCOEF(K,3)KCOEF(3) + LKCOEF(K,4)KCOEF(4)
[253]	197
[135]	198	TAUGAS = U(k)*ANS
	199
[253]	200	TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS + DCONT
	201	!TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS
[135]	202
[253]	203	DTAUKI(K,nw,ng) = TAUGAS + DCONT ! For parameterized continuum absorption
	204
[135]	205	do iaer=1,naerkind
[253]	206	DTAUKI(K,nw,ng) = DTAUKI(K,nw,ng) + TAEROS(K,NW,IAER)
	207	end do ! a bug was here!
[135]	208
	209	end do
	210
	211	! Now fill in the "clear" part of the spectrum (NG = L_NGAUSS),
	212	! which holds continuum opacity only
	213
	214	NG = L_NGAUSS
[253]	215	DTAUKI(K,nw,ng) = 0.0 + DCONT ! For parameterized continuum absorption
	216
[135]	217	do iaer=1,naerkind
[253]	218	DTAUKI(K,nw,ng) = DTAUKI(K,nw,ng) + TAEROS(K,NW,IAER)
	219	end do ! a bug was here!
[135]	220
	221	end do
	222	end do
	223
	224
	225	!=======================================================================
	226	! Now the full treatment for the layers, where besides the opacity
	227	! we need to calculate the scattering albedo and asymmetry factors
	228
	229	DO NW=1,L_NSPECTI
	230	DO K=2,L_LEVELS+1
	231	do iaer=1,naerkind
	232	TAUAEROLK(K,NW,IAER) = TAUAERO(K,IAER)*QSIAER(K,NW,IAER)
	233	end do
	234	ENDDO
	235	ENDDO
	236
	237	DO NW=1,L_NSPECTI
	238	NG = L_NGAUSS
	239	DO L=1,L_NLAYRAD
	240
	241	K = 2*L+1
	242	DTAUI(L,nw,ng) = DTAUKI(K,NW,NG) + DTAUKI(K+1,NW,NG)! + 1.e-50
	243
	244	atemp = 0.
	245	btemp = 0.
	246	if(DTAUI(L,NW,NG) .GT. 1.0E-9) then
	247	do iaer=1,naerkind
	248	atemp = atemp + &
	249	GIAER(K,NW,IAER) * TAUAEROLK(K,NW,IAER) + &
	250	GIAER(K+1,NW,IAER) * TAUAEROLK(K+1,NW,IAER)
	251	btemp = btemp + TAUAEROLK(K,NW,IAER) + TAUAEROLK(K+1,NW,IAER)
	252	! * + 1.e-10
	253	end do
	254	WBARI(L,nw,ng) = btemp / DTAUI(L,NW,NG)
	255	else
	256	WBARI(L,nw,ng) = 0.0D0
	257	DTAUI(L,NW,NG) = 1.0E-9
	258	endif
	259
	260	if(btemp .GT. 0.0) then
	261	cosbi(L,NW,NG) = atemp/btemp
	262	else
	263	cosbi(L,NW,NG) = 0.0D0
	264	end if
	265
	266	END DO ! L vertical loop
	267
	268	! Now the other Gauss points, if needed.
	269
	270	DO NG=1,L_NGAUSS-1
	271	IF(TAUGSURF(NW,NG) .gt. TLIMIT) THEN
	272
	273	DO L=1,L_NLAYRAD
	274	K = 2*L+1
	275	DTAUI(L,nw,ng) = DTAUKI(K,NW,NG)+DTAUKI(K+1,NW,NG)! + 1.e-50
	276
	277	btemp = 0.
	278	if(DTAUI(L,NW,NG) .GT. 1.0E-9) then
	279
	280	do iaer=1,naerkind
	281	btemp = btemp + TAUAEROLK(K,NW,IAER) + TAUAEROLK(K+1,NW,IAER)
	282	end do
	283	WBARI(L,nw,ng) = btemp / DTAUI(L,NW,NG)
	284
	285	else
	286	WBARI(L,nw,ng) = 0.0D0
	287	DTAUI(L,NW,NG) = 1.0E-9
	288	endif
	289
	290	cosbi(L,NW,NG) = cosbi(L,NW,L_NGAUSS)
	291	END DO ! L vertical loop
	292	END IF
	293
	294	END DO ! NG Gauss loop
	295	END DO ! NW spectral loop
	296
	297	! Total extinction optical depths
	298
	299	DO NW=1,L_NSPECTI
	300	DO NG=1,L_NGAUSS ! full gauss loop
	301	TAUI(1,NW,NG)=0.0D0
	302	DO L=1,L_NLAYRAD
	303	TAUI(L+1,NW,NG)=TAUI(L,NW,NG)+DTAUI(L,NW,NG)
	304	END DO
	305
	306	TAUCUMI(1,NW,NG)=0.0D0
	307	DO K=2,L_LEVELS
	308	TAUCUMI(K,NW,NG)=TAUCUMI(K-1,NW,NG)+DTAUKI(K,NW,NG)
	309	END DO
	310	END DO ! end full gauss loop
	311	END DO
	312
	313
	314	return
	315	end subroutine optci
	316
	317
	318	!-------------------------------------------------------------------------
[253]	319	subroutine water_cont(p,T,wratio,nu,DCONT)
[135]	320	! Calculates the continuum opacity for H2O
[253]	321	! NOT CURRENTLY USED
[135]	322
	323	implicit none
	324
[253]	325	real p, T, wratio, nu, DCONT
[135]	326	real h1, h2
	327
	328	h1 = exp(1800.*(1./T - 0.0034))
	329	h2 = 1.25e-22 + 1.67e-19exp(-2.62e-13nu)
	330
[253]	331	! DCONT = h1h2p(pwratio)*2/(RT)
	332	! DCONT=0.0 ! to be implemented...
[135]	333
	334	return
	335
	336	end subroutine water_cont
[253]	337

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

Download in other formats:

Original Format