Context Navigation

optci.F90 @ 2552

Last change on this file since 2552 was 2543, checked in by aslmd, 4 years ago

Generic GCM:

Adding k-coefficients mixing on the fly
Working with MordernTrac?

JVO + YJ

File size: 15.0 KB

Rev	Line
[2032]	1	MODULE optci_mod
	2
	3	IMPLICIT NONE
	4
	5	CONTAINS
	6
[716]	7	subroutine optci(PLEV,TLEV,DTAUI,TAUCUMI, &
	8	QXIAER,QSIAER,GIAER,COSBI,WBARI,TAUAERO, &
	9	TMID,PMID,TAUGSURF,QVAR,MUVAR)
[135]	10
[2032]	11	use radinc_h, only: L_LEVELS, L_NLAYRAD, L_NSPECTI, L_NGAUSS, &
	12	L_NLEVRAD, L_REFVAR, naerkind
[2133]	13	use radcommon_h, only: gasi,tlimit,wrefVAR,Cmk,tgasref,pfgasref,wnoi,scalep,indi,glat_ig
[2032]	14	use gases_h, only: gfrac, ngasmx, igas_N2, igas_He, igas_H2O, igas_H2
[1384]	15	use comcstfi_mod, only: g, r, mugaz
[2520]	16	use callkeys_mod, only: kastprof,continuum,graybody
[2543]	17	use recombin_corrk_mod, only: corrk_recombin, gasi_recomb
[716]	18	implicit none
[135]	19
[716]	20	!==================================================================
	21	!
	22	! Purpose
	23	! -------
	24	! Calculates longwave optical constants at each level. For each
	25	! layer and spectral interval in the IR it calculates WBAR, DTAU
	26	! and COSBAR. For each level it calculates TAU.
	27	!
[2133]	28	! TAUCUMI(L,LW) is the cumulative optical depth at level L (or alternatively
[716]	29	! at the bottom of layer L), LW is the spectral wavelength interval.
	30	!
	31	! TLEV(L) - Temperature at the layer boundary (i.e., level)
	32	! PLEV(L) - Pressure at the layer boundary (i.e., level)
	33	!
	34	! Authors
	35	! -------
	36	! Adapted from the NASA Ames code by R. Wordsworth (2009)
	37	!
	38	!==================================================================
[135]	39
	40
[716]	41	real*8 DTAUI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
[1715]	42	real*8 DTAUKI(L_LEVELS,L_NSPECTI,L_NGAUSS)
[716]	43	real*8 TAUI(L_NLEVRAD,L_NSPECTI,L_NGAUSS)
	44	real*8 TAUCUMI(L_LEVELS,L_NSPECTI,L_NGAUSS)
	45	real*8 PLEV(L_LEVELS)
	46	real*8 TLEV(L_LEVELS)
	47	real*8 TMID(L_LEVELS), PMID(L_LEVELS)
	48	real*8 COSBI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
	49	real*8 WBARI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
[135]	50
[716]	51	! for aerosols
[1715]	52	real*8 QXIAER(L_LEVELS,L_NSPECTI,NAERKIND)
	53	real*8 QSIAER(L_LEVELS,L_NSPECTI,NAERKIND)
	54	real*8 GIAER(L_LEVELS,L_NSPECTI,NAERKIND)
	55	real*8 TAUAERO(L_LEVELS,NAERKIND)
	56	real*8 TAUAEROLK(L_LEVELS,L_NSPECTI,NAERKIND)
[716]	57	real*8 TAEROS(L_LEVELS,L_NSPECTI,NAERKIND)
[135]	58
[716]	59	integer L, NW, NG, K, LK, IAER
	60	integer MT(L_LEVELS), MP(L_LEVELS), NP(L_LEVELS)
	61	real*8 ANS, TAUGAS
	62	real*8 DPR(L_LEVELS), U(L_LEVELS)
	63	real*8 LCOEF(4), LKCOEF(L_LEVELS,4)
[135]	64
[716]	65	real*8 taugsurf(L_NSPECTI,L_NGAUSS-1)
[918]	66	real*8 DCONT,DAERO
[716]	67	double precision wn_cont, p_cont, p_air, T_cont, dtemp, dtempc
	68	double precision p_cross
[135]	69
[716]	70	! variable species mixing ratio variables
	71	real*8 QVAR(L_LEVELS), WRATIO(L_LEVELS), MUVAR(L_LEVELS)
	72	real*8 KCOEF(4)
	73	integer NVAR(L_LEVELS)
[1725]	74
	75	! temporary variables to reduce memory access time to gasi
	76	real*8 tmpk(2,2)
	77	real*8 tmpkvar(2,2,2)
[135]	78
[716]	79	! temporary variables for multiple aerosol calculation
[918]	80	real*8 atemp
	81	real*8 btemp(L_NLAYRAD,L_NSPECTI)
[135]	82
[716]	83	! variables for k in units m^-1
[873]	84	real*8 dz(L_LEVELS)
	85	!real*8 rho !! see test below
[135]	86
[716]	87	integer igas, jgas
[253]	88
[873]	89	integer interm
	90
[716]	91	!--- Kasting's CIA ----------------------------------------
	92	!real*8, parameter :: Ci(L_NSPECTI)=[ &
	93	! 3.8E-5, 1.2E-5, 2.8E-6, 7.6E-7, 4.5E-7, 2.3E-7, &
	94	! 5.4E-7, 1.6E-6, 0.0, &
	95	! 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
	96	! 0.0, 4.0E-7, 4.0E-6, 1.4E-5, &
	97	! 1.0E-5, 1.2E-6, 2.0E-7, 5.0E-8, 3.0E-8, 0.0 ]
	98	!real*8, parameter :: Ti(L_NSPECTI)=[ -2.2, -1.9, &
	99	! -1.7, -1.7, -1.7, -1.7, -1.7, -1.7, &
	100	! 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
	101	! -1.7,-1.7,-1.7,-1.7,-1.7,-1.7,-1.7, -1.7,0.0 ]
	102	!----------------------------------------------------------
[253]	103
[873]	104	!! AS: to save time in computing continuum (see bilinearbig)
	105	IF (.not.ALLOCATED(indi)) THEN
[878]	106	ALLOCATE(indi(L_NSPECTI,ngasmx,ngasmx))
[873]	107	indi = -9999 ! this initial value means "to be calculated"
	108	ENDIF
	109
[716]	110	!=======================================================================
	111	! Determine the total gas opacity throughout the column, for each
	112	! spectral interval, NW, and each Gauss point, NG.
[135]	113
[716]	114	taugsurf(:,:) = 0.0
	115	dpr(:) = 0.0
	116	lkcoef(:,:) = 0.0
[135]	117
[716]	118	do K=2,L_LEVELS
	119	DPR(k) = PLEV(K)-PLEV(K-1)
[135]	120
[716]	121	!--- Kasting's CIA ----------------------------------------
	122	!dz(k)=dpr(k)189.02TMID(K)/(0.03720*PMID(K))
	123	! this is CO2 path length (in cm) as written by Francois
	124	! delta_z = delta_p * R_specific * T / (g * P)
	125	! But Kasting states that W is in units of _atmosphere_ cm
	126	! So we do
	127	!dz(k)=dz(k)*(PMID(K)/1013.25)
	128	!dz(k)=dz(k)/100.0 ! in m for SI calc
	129	!----------------------------------------------------------
[135]	130
[716]	131	! if we have continuum opacities, we need dz
	132	if(kastprof)then
[961]	133	dz(k) = dpr(k)(1000.0d08.3145d0/muvar(k))TMID(K)/(gPMID(K))
[1016]	134	U(k) = CmkDPR(k)mugaz/muvar(k)
[716]	135	else
[1194]	136	dz(k) = dpr(k)RTMID(K)/(glat_igPMID(K))mugaz/muvar(k)
[1016]	137	U(k) = CmkDPR(k)mugaz/muvar(k) ! only Cmk line in optci.F
	138	!JL13 the mugaz/muvar factor takes into account water meanmolecular weight if water is present
[716]	139	endif
[135]	140
[716]	141	call tpindex(PMID(K),TMID(K),QVAR(K),pfgasref,tgasref,WREFVAR, &
	142	LCOEF,MT(K),MP(K),NVAR(K),WRATIO(K))
[135]	143
[716]	144	do LK=1,4
	145	LKCOEF(K,LK) = LCOEF(LK)
	146	end do
[918]	147	end do ! levels
[253]	148
[1715]	149	! Spectral dependance of aerosol absorption
[918]	150	do iaer=1,naerkind
[716]	151	DO NW=1,L_NSPECTI
[918]	152	do K=2,L_LEVELS
[716]	153	TAEROS(K,NW,IAER) = TAUAERO(K,IAER) * QXIAER(K,NW,IAER)
[918]	154	end do ! levels
[716]	155	END DO
[918]	156	end do
[135]	157
[918]	158	do NW=1,L_NSPECTI
[135]	159
[918]	160	do K=2,L_LEVELS
[1715]	161
	162	DAERO=SUM(TAEROS(K,NW,1:naerkind)) ! aerosol absorption
[873]	163
[1715]	164	DCONT = 0.0d0 ! continuum absorption
[135]	165
[873]	166	if(continuum.and.(.not.graybody))then
[716]	167	! include continua if necessary
	168	wn_cont = dble(wnoi(nw))
	169	T_cont = dble(TMID(k))
	170	do igas=1,ngasmx
[135]	171
[716]	172	if(gfrac(igas).eq.-1)then ! variable
	173	p_cont = dble(PMID(k)scalepQVAR(k)) ! qvar = mol/mol
	174	else
	175	p_cont = dble(PMID(k)scalepgfrac(igas)*(1.-QVAR(k)))
	176	endif
[253]	177
[961]	178	dtemp=0.0d0
[716]	179	if(igas.eq.igas_N2)then
[305]	180
[878]	181	interm = indi(nw,igas,igas)
	182	call interpolateN2N2(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
	183	indi(nw,igas,igas) = interm
[253]	184
[716]	185	elseif(igas.eq.igas_H2)then
[253]	186
[716]	187	! first do self-induced absorption
[878]	188	interm = indi(nw,igas,igas)
[873]	189	call interpolateH2H2(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
[878]	190	indi(nw,igas,igas) = interm
[253]	191
[716]	192	! then cross-interactions with other gases
	193	do jgas=1,ngasmx
	194	p_cross = dble(PMID(k)scalepgfrac(jgas)*(1.-QVAR(k)))
[961]	195	dtempc = 0.0d0
[716]	196	if(jgas.eq.igas_N2)then
[878]	197	interm = indi(nw,igas,jgas)
	198	call interpolateN2H2(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
	199	indi(nw,igas,jgas) = interm
[716]	200	elseif(jgas.eq.igas_He)then
[878]	201	interm = indi(nw,igas,jgas)
[873]	202	call interpolateH2He(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
[878]	203	indi(nw,igas,jgas) = interm
[716]	204	endif
	205	dtemp = dtemp + dtempc
	206	enddo
[135]	207
[2520]	208	elseif(igas.eq.igas_H2O.and.T_cont.gt.100.0)then
	209	! Compute self and foreign (with air) continuum of H2O
[716]	210	p_air = dble(PMID(k)*scalep) - p_cont ! note assumes background is air!
[2520]	211	interm = indi(nw,igas,igas)
	212	call interpolateH2O_self_foreign(wn_cont,T_cont,p_cont,p_air,dtemp,.false.,interm) ! MTCKD v3.3
	213	indi(nw,igas,igas) = interm
[135]	214
[716]	215	endif
[135]	216
[716]	217	DCONT = DCONT + dtemp
[135]	218
[716]	219	enddo
[135]	220
[716]	221	! Oobleck test
	222	!rho = PMID(k)scalep / (TMID(k)286.99)
	223	!if(WNOI(nw).gt.300.0 .and. WNOI(nw).lt.500.0)then
	224	! DCONT = rho * 0.125 * 4.6e-4
	225	!elseif(WNOI(nw).gt.500.0 .and. WNOI(nw).lt.700.0)then
	226	! DCONT = 1000dpr(k) 1.0 * 4.6e-4 / g
	227	! DCONT = rho * 1.0 * 4.6e-4
	228	!elseif(WNOI(nw).gt.700.0 .and. WNOI(nw).lt.900.0)then
	229	! DCONT = rho * 0.125 * 4.6e-4
	230	!endif
[135]	231
[716]	232	DCONT = DCONT*dz(k)
[135]	233
[716]	234	endif
[135]	235
[716]	236	do ng=1,L_NGAUSS-1
[135]	237
[716]	238	! Now compute TAUGAS
[135]	239
[716]	240	! Interpolate between water mixing ratios
	241	! WRATIO = 0.0 if the requested water amount is equal to, or outside the
	242	! the water data range
[253]	243
[716]	244	if(L_REFVAR.eq.1)then ! added by RW for special no variable case
[1725]	245
	246	! JVO 2017 : added tmpk because the repeated calls to gasi/v increased dramatically
	247	! the execution time of optci/v -> ~ factor 2 on the whole radiative
	248	! transfer on the tested simulations !
	249
[2543]	250	IF (corrk_recombin) THEN ! added by JVO
	251	tmpk = GASI_RECOMB(MT(K):MT(K)+1,MP(K):MP(K)+1,1,NW,NG) ! contains the mix of recombined species
	252	ELSE
	253	tmpk = GASI(MT(K):MT(K)+1,MP(K):MP(K)+1,1,NW,NG)
	254	ENDIF
[1725]	255
	256	KCOEF(1) = tmpk(1,1) ! KCOEF(1) = GASI(MT(K),MP(K),1,NW,NG)
	257	KCOEF(2) = tmpk(1,2) ! KCOEF(2) = GASI(MT(K),MP(K)+1,1,NW,NG)
	258	KCOEF(3) = tmpk(2,2) ! KCOEF(3) = GASI(MT(K)+1,MP(K)+1,1,NW,NG)
	259	KCOEF(4) = tmpk(2,1) ! KCOEF(4) = GASI(MT(K)+1,MP(K),1,NW,NG)
	260
[716]	261	else
[135]	262
[2543]	263	IF (corrk_recombin) THEN ! added by JVO
	264	tmpkvar = GASI_RECOMB(MT(K):MT(K)+1,MP(K):MP(K)+1,NVAR(K):NVAR(K)+1,NW,NG)
	265	ELSE
	266	tmpkvar = GASI(MT(K):MT(K)+1,MP(K):MP(K)+1,NVAR(K):NVAR(K)+1,NW,NG)
	267	ENDIF
[135]	268
[1725]	269	KCOEF(1) = tmpkvar(1,1,1) + WRATIO(K) * &
	270	( tmpkvar(1,1,2)-tmpkvar(1,1,1) )
[135]	271
[1725]	272	KCOEF(2) = tmpkvar(1,2,1) + WRATIO(K) * &
	273	( tmpkvar(1,2,2)-tmpkvar(1,2,1) )
[135]	274
[1725]	275	KCOEF(3) = tmpkvar(2,2,1) + WRATIO(K) * &
	276	( tmpkvar(2,2,2)-tmpkvar(2,2,1) )
	277
	278	KCOEF(4) = tmpkvar(2,1,1) + WRATIO(K) * &
	279	( tmpkvar(2,1,2)-tmpkvar(2,1,1) )
[873]	280
[716]	281	endif
[135]	282
[716]	283	! Interpolate the gaseous k-coefficients to the requested T,P values
[135]	284
[716]	285	ANS = LKCOEF(K,1)KCOEF(1) + LKCOEF(K,2)KCOEF(2) + &
	286	LKCOEF(K,3)KCOEF(3) + LKCOEF(K,4)KCOEF(4)
[135]	287
[716]	288	TAUGAS = U(k)*ANS
[135]	289
[716]	290	TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS + DCONT
[918]	291	DTAUKI(K,nw,ng) = TAUGAS &
	292	+ DCONT & ! For parameterized continuum absorption
	293	+ DAERO ! For aerosol absorption
[135]	294
[716]	295	end do
[135]	296
[716]	297	! Now fill in the "clear" part of the spectrum (NG = L_NGAUSS),
	298	! which holds continuum opacity only
[135]	299
[716]	300	NG = L_NGAUSS
[918]	301	DTAUKI(K,nw,ng) = 0.d0 &
	302	+ DCONT & ! For parameterized continuum absorption
	303	+ DAERO ! For aerosol absorption
[135]	304
[716]	305	end do
	306	end do
[135]	307
[716]	308	!=======================================================================
	309	! Now the full treatment for the layers, where besides the opacity
	310	! we need to calculate the scattering albedo and asymmetry factors
[135]	311
[873]	312	do iaer=1,naerkind
[918]	313	DO NW=1,L_NSPECTI
[1715]	314	DO K=2,L_LEVELS
	315	TAUAEROLK(K,NW,IAER) = TAUAERO(K,IAER)*QSIAER(K,NW,IAER) ! effect of scattering albedo
[716]	316	ENDDO
[918]	317	ENDDO
[873]	318	end do
[918]	319
	320	DO NW=1,L_NSPECTI
[1715]	321	DO L=1,L_NLAYRAD-1
[918]	322	K = 2*L+1
	323	btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind)) + SUM(TAUAEROLK(K+1,NW,1:naerkind))
	324	END DO ! L vertical loop
[1715]	325
	326	! Last level
	327	L = L_NLAYRAD
	328	K = 2*L+1
	329	btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind))
	330
[918]	331	END DO ! NW spectral loop
	332
[135]	333
[716]	334	DO NW=1,L_NSPECTI
	335	NG = L_NGAUSS
[1715]	336	DO L=1,L_NLAYRAD-1
[135]	337
[716]	338	K = 2*L+1
	339	DTAUI(L,nw,ng) = DTAUKI(K,NW,NG) + DTAUKI(K+1,NW,NG)! + 1.e-50
[135]	340
[716]	341	atemp = 0.
[961]	342	if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
[716]	343	do iaer=1,naerkind
	344	atemp = atemp + &
	345	GIAER(K,NW,IAER) * TAUAEROLK(K,NW,IAER) + &
	346	GIAER(K+1,NW,IAER) * TAUAEROLK(K+1,NW,IAER)
	347	end do
[918]	348	WBARI(L,nw,ng) = btemp(L,nw) / DTAUI(L,NW,NG)
[716]	349	else
	350	WBARI(L,nw,ng) = 0.0D0
[961]	351	DTAUI(L,NW,NG) = 1.0D-9
[716]	352	endif
[135]	353
[961]	354	if(btemp(L,nw) .GT. 0.0d0) then
[918]	355	cosbi(L,NW,NG) = atemp/btemp(L,nw)
[716]	356	else
	357	cosbi(L,NW,NG) = 0.0D0
	358	end if
[135]	359
[716]	360	END DO ! L vertical loop
[1715]	361
	362	! Last level
	363
	364	L = L_NLAYRAD
	365	K = 2*L+1
	366	DTAUI(L,nw,ng) = DTAUKI(K,NW,NG) ! + 1.e-50
[135]	367
[1715]	368	atemp = 0.
	369	if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
	370	do iaer=1,naerkind
	371	atemp = atemp + GIAER(K,NW,IAER) * TAUAEROLK(K,NW,IAER)
	372	end do
	373	WBARI(L,nw,ng) = btemp(L,nw) / DTAUI(L,NW,NG)
	374	else
	375	WBARI(L,nw,ng) = 0.0D0
	376	DTAUI(L,NW,NG) = 1.0D-9
	377	endif
	378
	379	if(btemp(L,nw) .GT. 0.0d0) then
	380	cosbi(L,NW,NG) = atemp/btemp(L,nw)
	381	else
	382	cosbi(L,NW,NG) = 0.0D0
	383	end if
	384
	385
[716]	386	! Now the other Gauss points, if needed.
[135]	387
[716]	388	DO NG=1,L_NGAUSS-1
	389	IF(TAUGSURF(NW,NG) .gt. TLIMIT) THEN
[135]	390
[1715]	391	DO L=1,L_NLAYRAD-1
[716]	392	K = 2*L+1
	393	DTAUI(L,nw,ng) = DTAUKI(K,NW,NG)+DTAUKI(K+1,NW,NG)! + 1.e-50
	394
[961]	395	if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
[716]	396
[918]	397	WBARI(L,nw,ng) = btemp(L,nw) / DTAUI(L,NW,NG)
[716]	398
	399	else
	400	WBARI(L,nw,ng) = 0.0D0
[961]	401	DTAUI(L,NW,NG) = 1.0D-9
[716]	402	endif
	403
	404	cosbi(L,NW,NG) = cosbi(L,NW,L_NGAUSS)
	405	END DO ! L vertical loop
[1715]	406
	407	! Last level
	408	L = L_NLAYRAD
	409	K = 2*L+1
	410	DTAUI(L,nw,ng) = DTAUKI(K,NW,NG)! + 1.e-50
	411
	412	if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
	413
	414	WBARI(L,nw,ng) = btemp(L,nw) / DTAUI(L,NW,NG)
	415
	416	else
	417	WBARI(L,nw,ng) = 0.0D0
	418	DTAUI(L,NW,NG) = 1.0D-9
	419	endif
	420
	421	cosbi(L,NW,NG) = cosbi(L,NW,L_NGAUSS)
	422
[716]	423	END IF
	424
	425	END DO ! NG Gauss loop
	426	END DO ! NW spectral loop
	427
	428	! Total extinction optical depths
	429
[918]	430	DO NG=1,L_NGAUSS ! full gauss loop
	431	DO NW=1,L_NSPECTI
[716]	432	TAUCUMI(1,NW,NG)=0.0D0
	433	DO K=2,L_LEVELS
	434	TAUCUMI(K,NW,NG)=TAUCUMI(K-1,NW,NG)+DTAUKI(K,NW,NG)
	435	END DO
	436	END DO ! end full gauss loop
	437	END DO
	438
	439	! be aware when comparing with textbook results
	440	! (e.g. Pierrehumbert p. 218) that
	441	! taucumi does not take the <cos theta>=0.5 factor into
	442	! account. It is the optical depth for a vertically
	443	! ascending ray with angle theta = 0.
	444
	445	!open(127,file='taucum.out')
	446	!do nw=1,L_NSPECTI
	447	! write(127,*) taucumi(L_LEVELS,nw,L_NGAUSS)
	448	!enddo
	449	!close(127)
[918]	450
	451	! print*,'WBARI'
	452	! print*,WBARI
	453	! print*,'DTAUI'
	454	! print*,DTAUI
	455	! call abort
[2131]	456
[716]	457	end subroutine optci
	458
[2032]	459	END MODULE optci_mod
[716]	460

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

Download in other formats:

Original Format