Context Navigation

optci.F90 @ 486

Last change on this file since 486 was 471, checked in by aslmd, 13 years ago

LMDZ.GENERIC

13/12/2011 == AS

Same spirit as previous commit, but for ngasmx which is now read in gases.def -- before arrays w/ dim ngasmx are allocated dynamically
Allocation is done in su_gases.F90 which is called in inifis
Outside su_gases.F90, very few modifications to the code : the new module "gases_h.F90" simply replaces the old common "gases.h" !
Compiles fine. Tested with debugging options through pgdbg. Runs fine. Exact same results in Early Mars test case.

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
[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