Context Navigation

optci.F90 @ 3436

Last change on this file since 3436 was 3373, checked in by afalco, 6 months ago
Pluto.old PCM: Corrected some runtime issues with gfortran. AF
File size: 9.4 KB

Rev	Line
[3175]	1	subroutine optci(PLEV,TLEV,DTAUI,TAUCUMI, &
	2	QXIAER,QSIAER,GIAER,COSBI,WBARI,TAUAERO, &
	3	TMID,PMID,TAUGSURF,QVAR)
	4
	5
	6	use radinc_h
	7	use radcommon_h, only: gasi, tlimit, wrefVAR, Cmk,tgasref,pfgasref,wnoi
	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)
	60	real*8 dco2
	61
	62	! mixing ratio variables
	63	real*8 QVAR(L_LEVELS), WRATIO(L_LEVELS)
	64	real*8 KCOEF(4)
	65	integer NVAR(L_LEVELS)
	66
	67	! temporary variables for multiple aerosol calculation
[3373]	68	real*8 atemp
	69	real*8 btemp(L_NLAYRAD,L_NSPECTI)
[3175]	70
	71	! variables for k in units m^-1
	72	real*8 rho, dz
	73
	74	!=======================================================================
	75	! Determine the total gas opacity throughout the column, for each
	76	! spectral interval, NW, and each Gauss point, NG.
	77
	78	! write(,)'L_LEVELS',L_LEVELS
	79	! write(,)'L_NSPECTI',L_NSPECTI
	80	DTAUI(:,:,:)=0.
	81	DTAUKI(:,:,:)=0.
	82
	83	DO NG=1,L_NGAUSS-1
	84	do NW=1,L_NSPECTI
	85	TAUGSURF(NW,NG) = 0.0D0
	86	end do
	87	end do
	88	do K=2,L_LEVELS
	89	DPR(k) = PLEV(K)-PLEV(K-1)
	90
	91	! rho = PLEV(K)/(R*TMID(K))
	92	rho = PMID(K)/(R*TMID(K))
	93	dz = -DPR(k)/(g*rho)
	94	!print*,'rho=',rho
	95	!print*,'dz=',dz
	96
	97	U(k) = Cmk*DPR(k) ! only Cmk line in optci.F
	98	! soon to be replaced by m^-1 !!!
	99
	100	call tpindex(PMID(K),TMID(K),QVAR(K),pfgasref,tgasref,WREFVAR, &
	101	LCOEF,MT(K),MP(K),NVAR(K),WRATIO(K))
	102
	103	do LK=1,4
	104	LKCOEF(K,LK) = LCOEF(LK)
	105	end do
	106
	107	DO NW=1,L_NSPECTI
	108	do iaer=1,naerkind
	109	TAEROS(K,NW,IAER) = TAUAERO(K,IAER) * QXIAER(K,NW,IAER)
	110	! write(22,*) 'TB17 Taero IR:',K,NW,IAER,TAEROS(K,NW,IAER)
	111	end do
	112	END DO
	113	end do ! levels
	114
	115
	116	do K=2,L_LEVELS
	117	do nw=1,L_NSPECTI
	118
	119	DCO2 = 0.0 ! continuum absorption (no longer used)
	120
	121	do ng=1,L_NGAUSS-1
	122
	123	! Now compute TAUGAS
	124	! Interpolate between mixing ratios
	125	! WRATIO = 0.0 if the requested amount is equal to, or outside the
	126	! the data range
	127
	128
	129	if (L_REFVAR.eq.1)then ! added by RW for special no variable case
	130	KCOEF(1) = GASI(MT(K),MP(K),1,NW,NG)
	131	KCOEF(2) = GASI(MT(K),MP(K)+1,1,NW,NG)
	132	KCOEF(3) = GASI(MT(K)+1,MP(K)+1,1,NW,NG)
	133	KCOEF(4) = GASI(MT(K)+1,MP(K),1,NW,NG)
	134	else
	135
	136	KCOEF(1) = GASI(MT(K),MP(K),NVAR(K),NW,NG) + WRATIO(K)* &
	137	(GASI(MT(K),MP(K),NVAR(K)+1,NW,NG) - &
	138	GASI(MT(K),MP(K),NVAR(K),NW,NG))
	139
	140	KCOEF(2) = GASI(MT(K),MP(K)+1,NVAR(K),NW,NG)+ WRATIO(K)* &
	141	(GASI(MT(K),MP(K)+1,NVAR(K)+1,NW,NG) - &
	142	GASI(MT(K),MP(K)+1,NVAR(K),NW,NG))
	143
	144	KCOEF(3)=GASI(MT(K)+1,MP(K)+1,NVAR(K),NW,NG)+WRATIO(K)* &
	145	(GASI(MT(K)+1,MP(K)+1,NVAR(K)+1,NW,NG) - &
	146	GASI(MT(K)+1,MP(K)+1,NVAR(K),NW,NG))
	147
	148	KCOEF(4) =GASI(MT(K)+1,MP(K),NVAR(K),NW,NG) + WRATIO(K)* &
	149	(GASI(MT(K)+1,MP(K),NVAR(K)+1,NW,NG) - &
	150	GASI(MT(K)+1,MP(K),NVAR(K),NW,NG))
	151	endif
	152
	153	! Interpolate the gaseous k-coefficients to the requested T,P values
	154
	155	ANS = LKCOEF(K,1)KCOEF(1) + LKCOEF(K,2)KCOEF(2) + &
	156	LKCOEF(K,3)KCOEF(3) + LKCOEF(K,4)KCOEF(4)
	157
	158	TAUGAS = U(k)*ANS
	159
	160	TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS
	161
	162	DTAUKI(K,nw,ng) = TAUGAS
	163	do iaer=1,naerkind
	164	DTAUKI(K,nw,ng) = DTAUKI(K,nw,ng) + TAEROS(K,NW,IAER) &
	165	+ DCO2 ! For Kasting CIA
	166	end do
	167
	168	end do
	169
	170	! Now fill in the "clear" part of the spectrum (NG = L_NGAUSS),
	171	! which holds continuum opacity only
	172
	173	NG = L_NGAUSS
	174	DTAUKI(K,nw,ng) = 0.0
	175	do iaer=1,naerkind
	176	DTAUKI(K,nw,ng) = DTAUKI(K,nw,ng) + TAEROS(K,NW,IAER) &
	177	+ DCO2 ! For parameterized continuum absorption
	178	end do ! a bug was found here!!
	179
	180	end do
	181	end do
	182
	183
	184	!=======================================================================
	185	! Now the full treatment for the layers, where besides the opacity
	186	! we need to calculate the scattering albedo and asymmetry factors
	187
	188	DO NW=1,L_NSPECTI
[3373]	189	DO K=2,L_LEVELS
[3175]	190	do iaer=1,naerkind
	191	TAUAEROLK(K,NW,IAER) = TAUAERO(K,IAER)*QSIAER(K,NW,IAER)
	192	end do
	193	ENDDO
	194	ENDDO
	195
	196	DO NW=1,L_NSPECTI
	197	NG = L_NGAUSS
	198	DO L=1,L_NLAYRAD-1
	199
	200	K = 2*L+1
	201	DTAUI(L,nw,ng) = DTAUKI(K,NW,NG) + DTAUKI(K+1,NW,NG)! + 1.e-50
	202
	203	atemp = 0.
[3373]	204	btemp(L,NW) = 0.
	205	do iaer=1,naerkind
	206	atemp = atemp + &
	207	GIAER(K,NW,IAER) * TAUAEROLK(K,NW,IAER) + &
	208	GIAER(K+1,NW,IAER) * TAUAEROLK(K+1,NW,IAER)
	209	btemp(L,NW) = btemp(L,NW) + TAUAEROLK(K,NW,IAER) + TAUAEROLK(K+1,NW,IAER)
	210	! * + 1.e-10
	211	end do
	212
[3175]	213	if(DTAUI(L,NW,NG) .GT. 1.0E-9) then
[3373]	214	WBARI(L,nw,ng) = btemp(L,NW) / DTAUI(L,NW,NG)
[3175]	215	else
	216	WBARI(L,nw,ng) = 0.0D0
	217	DTAUI(L,NW,NG) = 1.0E-9
	218	endif
	219
[3373]	220	if(btemp(L,NW) .GT. 0.0) then
	221	cosbi(L,NW,NG) = atemp/btemp(L,NW)
[3175]	222	else
	223	cosbi(L,NW,NG) = 0.0D0
	224	end if
	225
	226	END DO ! L vertical loop
[3373]	227	! Last level
[3175]	228
[3373]	229	L = L_NLAYRAD
	230	K = 2*L+1
	231	DTAUI(L,nw,ng) = DTAUKI(K,NW,NG) ! + 1.e-50
	232	btemp(L,NW) = 0
	233	do iaer=1,naerkind
	234	btemp(L,NW) = btemp(L,NW) + TAUAEROLK(K,NW,IAER)
	235	enddo
	236
	237	atemp = 0.
	238	if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
	239	do iaer=1,naerkind
	240	atemp = atemp + GIAER(K,NW,IAER) * TAUAEROLK(K,NW,IAER)
	241	end do
	242	WBARI(L,nw,ng) = btemp(L,NW) / DTAUI(L,NW,NG)
	243	else
	244	WBARI(L,nw,ng) = 0.0D0
	245	DTAUI(L,NW,NG) = 1.0D-9
	246	endif
[3175]	247
[3373]	248	if(btemp(L,NW) .GT. 0.0d0) then
	249	cosbi(L,NW,NG) = atemp/btemp(L,NW)
	250	else
	251	cosbi(L,NW,NG) = 0.0D0
	252	end if
[3175]	253
[3373]	254	! Now the other Gauss points, if needed.
[3175]	255
	256	DO NG=1,L_NGAUSS-1
	257	IF(TAUGSURF(NW,NG) .gt. TLIMIT) THEN
	258
	259	DO L=1,L_NLAYRAD
	260	K = 2*L+1
	261	DTAUI(L,nw,ng) = DTAUKI(K,NW,NG)+DTAUKI(K+1,NW,NG)! + 1.e-50
	262
	263	if(DTAUI(L,NW,NG) .GT. 1.0E-9) then
[3373]	264	WBARI(L,nw,ng) = btemp(L,NW) / DTAUI(L,NW,NG)
[3175]	265	else
	266	WBARI(L,nw,ng) = 0.0D0
	267	DTAUI(L,NW,NG) = 1.0E-9
	268	endif
	269
	270	cosbi(L,NW,NG) = cosbi(L,NW,L_NGAUSS)
	271	END DO ! L vertical loop
	272	END IF
	273
	274	END DO ! NG Gauss loop
	275	END DO ! NW spectral loop
	276
	277	! Total extinction optical depths
	278
	279	DO NW=1,L_NSPECTI
	280	DO NG=1,L_NGAUSS ! full gauss loop
	281	TAUI(1,NW,NG)=0.0D0
	282	DO L=1,L_NLAYRAD
	283	TAUI(L+1,NW,NG)=TAUI(L,NW,NG)+DTAUI(L,NW,NG)
	284	END DO
	285
	286	TAUCUMI(1,NW,NG)=0.0D0
	287	DO K=2,L_LEVELS
	288	TAUCUMI(K,NW,NG)=TAUCUMI(K-1,NW,NG)+DTAUKI(K,NW,NG)
	289	END DO
	290	END DO ! end full gauss loop
	291	END DO
	292
	293
	294	return
	295	end subroutine optci
	296

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

Download in other formats:

Original Format