Context Navigation

optcv.F90 @ 601

Last change on this file since 601 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)

Property svn:executable set to *

File size: 10.2 KB

Rev	Line
[253]	1	SUBROUTINE OPTCV(DTAUV,TAUV,TAUCUMV,PLEV, &
	2	QXVAER,QSVAER,GVAER,WBARV,COSBV, &
[305]	3	TAURAY,TAUAERO,TMID,PMID,TAUGSURF,QVAR,MUVAR)
[253]	4
	5	use radinc_h
	6	use radcommon_h, only: gasv, tlimit, wrefVAR, Cmk, tgasref, pfgasref,wnov,scalep
[471]	7	use gases_h
[253]	8
	9	implicit none
	10
	11	!==================================================================
	12	!
	13	! Purpose
	14	! -------
	15	! Calculates shortwave optical constants at each level.
	16	!
	17	! Authors
	18	! -------
	19	! Adapted from the NASA Ames code by R. Wordsworth (2009)
	20	!
	21	!==================================================================
	22	!
	23	! THIS SUBROUTINE SETS THE OPTICAL CONSTANTS IN THE VISUAL
	24	! IT CALCUALTES FOR EACH LAYER, FOR EACH SPECRAL INTERVAL IN THE VISUAL
	25	! LAYER: WBAR, DTAU, COSBAR
	26	! LEVEL: TAU
	27	!
	28	! TAUV(L,NW,NG) is the cumulative optical depth at the top of radiation code
	29	! layer L. NW is spectral wavelength interval, ng the Gauss point index.
	30	!
	31	! TLEV(L) - Temperature at the layer boundary
	32	! PLEV(L) - Pressure at the layer boundary (i.e. level)
[305]	33	! GASV(NT,NPS,NW,NG) - Visible k-coefficients
[253]	34	!
	35	!-------------------------------------------------------------------
	36
	37	#include "callkeys.h"
	38	#include "comcstfi.h"
	39
	40	real*8 DTAUV(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
	41	real*8 DTAUKV(L_LEVELS+1,L_NSPECTV,L_NGAUSS)
	42	real*8 TAUV(L_NLEVRAD,L_NSPECTV,L_NGAUSS)
	43	real*8 TAUCUMV(L_LEVELS,L_NSPECTV,L_NGAUSS)
	44	real*8 PLEV(L_LEVELS)
	45	real*8 TMID(L_LEVELS), PMID(L_LEVELS)
	46	real*8 COSBV(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
	47	real*8 WBARV(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
	48	real*8 TAURAY(L_NSPECTV)
	49
	50	! For aerosols
	51	real*8 QXVAER(L_LEVELS+1,L_NSPECTV,NAERKIND)
	52	real*8 QSVAER(L_LEVELS+1,L_NSPECTV,NAERKIND)
	53	real*8 GVAER(L_LEVELS+1,L_NSPECTV,NAERKIND)
	54	real*8 TAUAERO(L_LEVELS+1,NAERKIND)
	55	real*8 TAUAEROLK(L_LEVELS+1,L_NSPECTV,NAERKIND)
	56	real*8 TAEROS(L_LEVELS,L_NSPECTV,NAERKIND)
	57
	58	integer L, NW, NG, K, NG1(L_NSPECTV), LK, IAER
	59	integer MT(L_LEVELS), MP(L_LEVELS), NP(L_LEVELS)
	60	real*8 ANS, TAUGAS
	61	real*8 TRAY(L_LEVELS,L_NSPECTV)
	62	real*8 DPR(L_LEVELS), U(L_LEVELS)
	63	real*8 LCOEF(4), LKCOEF(L_LEVELS,4)
	64
	65	real*8 taugsurf(L_NSPECTV,L_NGAUSS-1), TRAYAER
	66
[305]	67	! Variable species mixing ratio variables
	68	real*8 QVAR(L_LEVELS), WRATIO(L_LEVELS), MUVAR(L_LEVELS)
[253]	69	real*8 KCOEF(4)
	70	integer NVAR(L_LEVELS)
	71
	72	! temporary variables for multiple aerosol calculation
	73	real*8 atemp, btemp, ctemp
	74
	75	! variables for k in units m^-1
[305]	76	double precision wn_cont, p_cont, p_air, T_cont, dtemp
[253]	77	real*8 dz(L_LEVELS), DCONT
	78
	79	integer igas
	80
	81	!=======================================================================
	82	! Determine the total gas opacity throughout the column, for each
	83	! spectral interval, NW, and each Gauss point, NG.
	84	! Calculate the continuum opacities, i.e., those that do not depend on
	85	! NG, the Gauss index.
	86
	87	taugsurf(:,:) = 0.0
	88	dpr(:) = 0.0
	89	lkcoef(:,:) = 0.0
	90
	91	do K=2,L_LEVELS
	92	DPR(k) = PLEV(K)-PLEV(K-1)
	93
	94	! if we have continuum opacities, we need dz
[305]	95	if(kastprof)then
	96	dz(k) = dpr(k)(8314.5/muvar(k))TMID(K)/(g*PMID(K))
	97	U(k) = (Cmkmugaz/(muvar(k)))DPR(k)
	98	else
	99	dz(k) = dpr(k)RTMID(K)/(g*PMID(K))
	100	U(k) = Cmk*DPR(k) ! only Cmk line in optci.F
	101	endif
[253]	102
	103	call tpindex(PMID(K),TMID(K),QVAR(K),pfgasref,tgasref,WREFVAR, &
	104	LCOEF,MT(K),MP(K),NVAR(K),WRATIO(K))
	105
	106	do LK=1,4
	107	LKCOEF(K,LK) = LCOEF(LK)
	108	end do
	109
	110	DO NW=1,L_NSPECTV
	111	TRAY(K,NW) = TAURAY(NW) * DPR(K)
	112
	113	do iaer=1,naerkind
	114	TAEROS(K,NW,IAER) = TAUAERO(K,IAER) * QXVAER(K,NW,IAER)
	115	end do
	116	!
	117
	118	END DO
	119	end do
	120
	121	! TRAYAER is Tau RAYleigh scattering, plus AERosol opacity
	122
	123	! we ignore K=1...
	124	do K=2,L_LEVELS
	125	do NW=1,L_NSPECTV
	126
	127	TRAYAER = TRAY(K,NW)
	128	do iaer=1,naerkind
	129	TRAYAER = TRAYAER + TAEROS(K,NW,IAER)
	130	end do
	131
	132	DCONT = 0.0 ! continuum absorption
	133
[596]	134	if(callgasvis.and.Continuum.and.(.not.graybody))then
[305]	135	! include continua if necessary
[596]	136	wn_cont = dble(wnov(nw))
	137	T_cont = dble(TMID(k))
	138	do igas=1,ngasmx
	139
	140	if(gfrac(igas).eq.-1)then ! variable
	141	p_cont = dble(PMID(k)scalepQVAR(K)) ! qvar = mol/mol
	142	else
	143	p_cont = dble(PMID(k)scalepgfrac(igas)*(1.-QVAR(k)))
	144	endif
[305]	145
[596]	146	dtemp=0.0
	147	if(gnom(igas).eq.'H2_')then
	148	call interpolateH2H2(wn_cont,T_cont,p_cont,dtemp,.false.)
	149	elseif(gnom(igas).eq.'H2O'.and.T_cont.gt.200.0)then
	150	p_air = dble(PMID(k)*scalep) - p_cont ! note assumes air!!
	151	call interpolateH2Ocont(wn_cont,T_cont,p_cont,p_air,dtemp,.false.)
	152	endif
[305]	153
[596]	154	DCONT = DCONT + dtemp
[253]	155
[596]	156	enddo
[253]	157
[596]	158	DCONT = DCONT*dz(k)
[366]	159	endif
[253]	160
	161	do NG=1,L_NGAUSS-1
	162
	163	!=======================================================================
	164	! Now compute TAUGAS
	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
	169	if (L_REFVAR.eq.1)then ! added by RW for special no variable case
	170	KCOEF(1) = GASV(MT(K),MP(K),1,NW,NG)
	171	KCOEF(2) = GASV(MT(K),MP(K)+1,1,NW,NG)
	172	KCOEF(3) = GASV(MT(K)+1,MP(K)+1,1,NW,NG)
	173	KCOEF(4) = GASV(MT(K)+1,MP(K),1,NW,NG)
	174	else
	175	KCOEF(1) = GASV(MT(K),MP(K),NVAR(K),NW,NG) + WRATIO(K)* &
	176	(GASV(MT(K),MP(K),NVAR(K)+1,NW,NG) - &
	177	GASV(MT(K),MP(K),NVAR(K),NW,NG))
	178
	179	KCOEF(2) = GASV(MT(K),MP(K)+1,NVAR(K),NW,NG) + WRATIO(K)* &
	180	(GASV(MT(K),MP(K)+1,NVAR(K)+1,NW,NG) - &
	181	GASV(MT(K),MP(K)+1,NVAR(K),NW,NG))
	182
	183	KCOEF(3) = GASV(MT(K)+1,MP(K)+1,NVAR(K),NW,NG) + WRATIO(K)*&
	184	(GASV(MT(K)+1,MP(K)+1,NVAR(K)+1,NW,NG) - &
	185	GASV(MT(K)+1,MP(K)+1,NVAR(K),NW,NG))
	186
	187	KCOEF(4) = GASV(MT(K)+1,MP(K),NVAR(K),NW,NG) + WRATIO(K)* &
	188	(GASV(MT(K)+1,MP(K),NVAR(K)+1,NW,NG) - &
	189	GASV(MT(K)+1,MP(K),NVAR(K),NW,NG))
	190	endif
	191
	192	! Interpolate the gaseous k-coefficients to the requested T,P values
	193
	194	ANS = LKCOEF(K,1)KCOEF(1) + LKCOEF(K,2)KCOEF(2) + &
	195	LKCOEF(K,3)KCOEF(3) + LKCOEF(K,4)KCOEF(4)
	196
	197	TAUGAS = U(k)*ANS
	198
[305]	199
[253]	200	!TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS
	201	TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS + DCONT
	202	DTAUKV(K,nw,ng) = TAUGAS + TRAYAER & ! TRAYAER includes all scattering contributions
	203	+ DCONT ! for continuum absorption
	204
	205	end do
	206
	207
	208	! Now fill in the "clear" part of the spectrum (NG = L_NGAUSS),
	209	! which holds continuum opacity only
	210
	211	NG = L_NGAUSS
	212	DTAUKV(K,nw,ng) = TRAY(K,NW) + DCONT ! For parameterized continuum absorption
	213	do iaer=1,naerkind
	214	DTAUKV(K,nw,ng) = DTAUKV(K,nw,ng) + TAEROS(K,NW,IAER)
	215	! & + DCONT ! For parameterized continuum absorption
	216	end do ! a bug was here!
	217
	218	end do
	219	end do
	220
	221
	222	!=======================================================================
	223	! Now the full treatment for the layers, where besides the opacity
	224	! we need to calculate the scattering albedo and asymmetry factors
	225
	226	DO NW=1,L_NSPECTV
	227	DO K=2,L_LEVELS
	228	do iaer=1,naerkind
	229	TAUAEROLK(K,NW,IAER) = TAUAERO(K,IAER) * QSVAER(K,NW,IAER)
	230	end do
	231	ENDDO
	232	ENDDO
	233
	234
	235	DO NW=1,L_NSPECTV
	236	DO NG=1,L_NGAUSS
	237	DO L=1,L_NLAYRAD-1
	238	K = 2*L+1
	239
	240	DTAUV(L,nw,ng) = DTAUKV(K,NW,NG)+DTAUKV(K+1,NW,NG)
	241
	242	atemp=0.
	243	btemp=TRAY(K,NW) + TRAY(K+1,NW)
	244	ctemp=0.9999*(TRAY(K,NW) + TRAY(K+1,NW))
	245	do iaer=1,naerkind
	246	atemp = atemp + &
	247	GVAER(K,NW,IAER) * TAUAEROLK(K,NW,IAER) + &
	248	GVAER(K+1,NW,IAER) * TAUAEROLK(K+1,NW,IAER)
	249	btemp = btemp + &
	250	TAUAEROLK(K,NW,IAER) + TAUAEROLK(K+1,NW,IAER)
	251	ctemp = ctemp + &
	252	TAUAEROLK(K,NW,IAER) + TAUAEROLK(K+1,NW,IAER)
	253	end do
	254
	255	COSBV(L,NW,NG) = atemp/btemp
	256	WBARV(L,nw,ng) = ctemp/DTAUV(L,nw,ng)
	257
	258	END DO
	259
	260	! No vertical averaging on bottom layer
	261
	262	L = L_NLAYRAD
	263	K = 2*L+1
	264	DTAUV(L,nw,ng) = DTAUKV(K,NW,NG)
	265
	266	atemp=0.
	267	btemp=TRAY(K,NW)
	268	ctemp=0.9999*TRAY(K,NW)
	269	do iaer=1,naerkind
	270	atemp = atemp + GVAER(K,NW,IAER) * TAUAEROLK(K,NW,IAER)
	271	btemp = btemp + TAUAEROLK(K,NW,IAER)
	272	ctemp = ctemp + TAUAEROLK(K,NW,IAER)
	273	end do
	274	COSBV(L,NW,NG) = atemp/btemp
	275	WBARV(L,nw,ng) = ctemp/DTAUV(L,nw,ng)
	276
	277	END DO ! NG gauss point loop
	278	END DO ! NW spectral loop
	279
	280
	281
	282	! Total extinction optical depths
	283
	284	DO NW=1,L_NSPECTV
	285	DO NG=1,L_NGAUSS ! full gauss loop
	286	TAUV(1,NW,NG)=0.0D0
	287	DO L=1,L_NLAYRAD
	288	TAUV(L+1,NW,NG)=TAUV(L,NW,NG)+DTAUV(L,NW,NG)
	289	END DO
	290
	291	TAUCUMV(1,NW,NG)=0.0D0
	292	DO K=2,L_LEVELS
	293	TAUCUMV(K,NW,NG)=TAUCUMV(K-1,NW,NG)+DTAUKV(K,NW,NG)
	294	END DO
	295	END DO ! end full gauss loop
	296	END DO
	297
	298
	299	RETURN
	300	END SUBROUTINE OPTCV

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

Download in other formats:

Original Format