Context Navigation

optcv.F90 @ 3436

Last change on this file since 3436 was 3233, checked in by gmilcareck, 9 months ago
Add the possibility to use a fixed vertical molar fraction profile for the collision-induced absorption like the correlated-k.
Property svn:executable set to ``*
File size: 16.3 KB

Rev	Line
[2032]	1	MODULE optcv_mod
	2
	3	IMPLICIT NONE
	4
	5	CONTAINS
	6
[716]	7	SUBROUTINE OPTCV(DTAUV,TAUV,TAUCUMV,PLEV, &
	8	QXVAER,QSVAER,GVAER,WBARV,COSBV, &
[3233]	9	TAURAY,TAUAERO,TMID,PMID,TAUGSURF,QVAR,MUVAR,FRACVAR)
[253]	10
[2032]	11	use radinc_h, only: L_NLAYRAD, L_NLEVRAD, L_LEVELS, L_NSPECTV, L_NGAUSS, L_REFVAR, NAERKIND
[2133]	12	use radcommon_h, only: gasv, tlimit, wrefVAR, Cmk, tgasref, pfgasref,wnov,scalep,indv,glat_ig
[2875]	13	use gases_h, only: gfrac, ngasmx, igas_H2, igas_H2O, igas_He, igas_N2, &
	14	igas_CH4, igas_CO2
[1384]	15	use comcstfi_mod, only: g, r, mugaz
[3233]	16	use callkeys_mod, only: kastprof,continuum,graybody,callgasvis,varspec
[2543]	17	use recombin_corrk_mod, only: corrk_recombin, gasv_recomb
[2582]	18	use tpindex_mod, only: tpindex
[253]	19
[716]	20	implicit none
[253]	21
[716]	22	!==================================================================
	23	!
	24	! Purpose
	25	! -------
	26	! Calculates shortwave optical constants at each level.
	27	!
	28	! Authors
	29	! -------
	30	! Adapted from the NASA Ames code by R. Wordsworth (2009)
	31	!
	32	!==================================================================
	33	!
	34	! THIS SUBROUTINE SETS THE OPTICAL CONSTANTS IN THE VISUAL
[1715]	35	! IT CALCULATES FOR EACH LAYER, FOR EACH SPECTRAL INTERVAL IN THE VISUAL
[716]	36	! LAYER: WBAR, DTAU, COSBAR
	37	! LEVEL: TAU
	38	!
	39	! TAUV(L,NW,NG) is the cumulative optical depth at the top of radiation code
	40	! layer L. NW is spectral wavelength interval, ng the Gauss point index.
	41	!
	42	! TLEV(L) - Temperature at the layer boundary
	43	! PLEV(L) - Pressure at the layer boundary (i.e. level)
	44	! GASV(NT,NPS,NW,NG) - Visible k-coefficients
	45	!
	46	!-------------------------------------------------------------------
[253]	47
	48
[2972]	49	real*8,intent(out) :: DTAUV(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
[1715]	50	real*8 DTAUKV(L_LEVELS,L_NSPECTV,L_NGAUSS)
[2972]	51	real*8,intent(out) :: TAUV(L_NLEVRAD,L_NSPECTV,L_NGAUSS)
	52	real*8,intent(out) :: TAUCUMV(L_LEVELS,L_NSPECTV,L_NGAUSS)
	53	real*8,intent(in) :: PLEV(L_LEVELS)
	54	real*8,intent(in) :: TMID(L_LEVELS), PMID(L_LEVELS)
	55	real*8,intent(out) :: COSBV(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
	56	real*8,intent(out) :: WBARV(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
[253]	57
[716]	58	! for aerosols
[2972]	59	real*8,intent(in) :: QXVAER(L_LEVELS,L_NSPECTV,NAERKIND)
	60	real*8,intent(in) :: QSVAER(L_LEVELS,L_NSPECTV,NAERKIND)
	61	real*8,intent(in) :: GVAER(L_LEVELS,L_NSPECTV,NAERKIND)
	62	real*8,intent(in) :: TAUAERO(L_LEVELS,NAERKIND)
	63
	64	! local arrays (saved for convenience as need be allocated)
	65	real*8,save,allocatable :: TAUAEROLK(:,:,:)
	66	real*8,save,allocatable :: TAEROS(:,:,:)
	67	!$OMP THREADPRIVATE(TAUAEROLK,TAEROS)
[253]	68
[873]	69	integer L, NW, NG, K, LK, IAER
[716]	70	integer MT(L_LEVELS), MP(L_LEVELS), NP(L_LEVELS)
	71	real*8 ANS, TAUGAS
[2972]	72	real*8,intent(in) :: TAURAY(L_NSPECTV)
[716]	73	real*8 TRAY(L_LEVELS,L_NSPECTV)
	74	real*8 DPR(L_LEVELS), U(L_LEVELS)
	75	real*8 LCOEF(4), LKCOEF(L_LEVELS,4)
[253]	76
[2972]	77	real*8,intent(out) :: taugsurf(L_NSPECTV,L_NGAUSS-1)
[918]	78	real*8 DCONT,DAERO
[1715]	79	real*8 DRAYAER
[873]	80	double precision wn_cont, p_cont, p_air, T_cont, dtemp, dtempc
	81	double precision p_cross
[253]	82
[716]	83	! variable species mixing ratio variables
[2972]	84	real*8,intent(in) :: QVAR(L_LEVELS)
	85	real*8,intent(in) :: MUVAR(L_LEVELS)
[3233]	86	real*8,intent(in) :: FRACVAR(ngasmx,L_LEVELS)
[2972]	87	real*8 :: WRATIO(L_LEVELS)
[873]	88	real*8 KCOEF(4)
[716]	89	integer NVAR(L_LEVELS)
[1725]	90
	91	! temporary variables to reduce memory access time to gasv
	92	real*8 tmpk(2,2)
	93	real*8 tmpkvar(2,2,2)
[253]	94
[716]	95	! temporary variables for multiple aerosol calculation
[918]	96	real*8 atemp(L_NLAYRAD,L_NSPECTV)
	97	real*8 btemp(L_NLAYRAD,L_NSPECTV)
	98	real*8 ctemp(L_NLAYRAD,L_NSPECTV)
[253]	99
[716]	100	! variables for k in units m^-1
[873]	101	real*8 dz(L_LEVELS)
[253]	102
[2131]	103
[716]	104	integer igas, jgas
[253]	105
[873]	106	integer interm
	107
[2972]	108	logical :: firstcall=.true.
	109	!$OMP THREADPRIVATE(firstcall)
	110
	111	if (firstcall) then
	112	! allocate local arrays of size "naerkind" (which are also
	113	! "saved" so that this is done only once in for all even if
	114	! we don't need to store the value from a time step to the next)
	115	allocate(TAUAEROLK(L_LEVELS,L_NSPECTV,NAERKIND))
	116	allocate(TAEROS(L_LEVELS,L_NSPECTV,NAERKIND))
	117	firstcall=.false.
	118	endif ! of if (firstcall)
	119
[873]	120	!! AS: to save time in computing continuum (see bilinearbig)
	121	IF (.not.ALLOCATED(indv)) THEN
[878]	122	ALLOCATE(indv(L_NSPECTV,ngasmx,ngasmx))
[873]	123	indv = -9999 ! this initial value means "to be calculated"
	124	ENDIF
	125
[716]	126	!=======================================================================
	127	! Determine the total gas opacity throughout the column, for each
	128	! spectral interval, NW, and each Gauss point, NG.
	129	! Calculate the continuum opacities, i.e., those that do not depend on
	130	! NG, the Gauss index.
[253]	131
[716]	132	taugsurf(:,:) = 0.0
	133	dpr(:) = 0.0
	134	lkcoef(:,:) = 0.0
[253]	135
[716]	136	do K=2,L_LEVELS
	137	DPR(k) = PLEV(K)-PLEV(K-1)
[253]	138
[716]	139	! if we have continuum opacities, we need dz
	140	if(kastprof)then
[1016]	141	dz(k) = dpr(k)(1000.0d08.3145d0/muvar(k))TMID(K)/(gPMID(K))
	142	U(k) = CmkDPR(k)mugaz/muvar(k)
[716]	143	else
[1194]	144	dz(k) = dpr(k)RTMID(K)/(glat_igPMID(K))mugaz/muvar(k)
[1016]	145	U(k) = CmkDPR(k)mugaz/muvar(k) ! only Cmk line in optci.F
	146	!JL13 the mugaz/muvar factor takes into account water meanmolecular weight if water is present
[716]	147	endif
[253]	148
[716]	149	call tpindex(PMID(K),TMID(K),QVAR(K),pfgasref,tgasref,WREFVAR, &
	150	LCOEF,MT(K),MP(K),NVAR(K),WRATIO(K))
[253]	151
[716]	152	do LK=1,4
	153	LKCOEF(K,LK) = LCOEF(LK)
	154	end do
[918]	155	end do ! levels
[253]	156
[1715]	157	! Spectral dependance of aerosol absorption
[1987]	158	!JL18 It seems to be good to have aerosols in the first "radiative layer" of the gcm in the IR
	159	! but visible does not handle very well diffusion in first layer.
	160	! The tauaero and tauray are thus set to 0 (a small value for rayleigh because the code crashes otherwise)
	161	! in the 4 first semilayers in optcv, but not optci.
	162	! This solves random variations of the sw heating at the model top.
[918]	163	do iaer=1,naerkind
	164	do NW=1,L_NSPECTV
[1987]	165	TAEROS(1:4,NW,IAER)=0.d0
	166	do K=5,L_LEVELS
[873]	167	TAEROS(K,NW,IAER) = TAUAERO(K,IAER) * QXVAER(K,NW,IAER)
[918]	168	end do ! levels
	169	end do
	170	end do
[1715]	171
	172	! Rayleigh scattering
[918]	173	do NW=1,L_NSPECTV
[1987]	174	TRAY(1:4,NW) = 1d-30
	175	do K=5,L_LEVELS
[873]	176	TRAY(K,NW) = TAURAY(NW) * DPR(K)
[918]	177	end do ! levels
	178	end do
	179
[716]	180	! we ignore K=1...
	181	do K=2,L_LEVELS
[873]	182
[716]	183	do NW=1,L_NSPECTV
[253]	184
[1715]	185	DRAYAER = TRAY(K,NW)
	186	! DRAYAER is Tau RAYleigh scattering, plus AERosol opacity
[716]	187	do iaer=1,naerkind
[1715]	188	DRAYAER = DRAYAER + TAEROS(K,NW,IAER)
[716]	189	end do
[253]	190
[716]	191	DCONT = 0.0 ! continuum absorption
[253]	192
[873]	193	if(continuum.and.(.not.graybody).and.callgasvis)then
[716]	194	! include continua if necessary
	195	wn_cont = dble(wnov(nw))
	196	T_cont = dble(TMID(k))
	197	do igas=1,ngasmx
[305]	198
[716]	199	if(gfrac(igas).eq.-1)then ! variable
	200	p_cont = dble(PMID(k)scalepQVAR(k)) ! qvar = mol/mol
[3233]	201	elseif(varspec) then
	202	p_cont = dble(PMID(k)scalepFRACVAR(igas,k)*(1.-QVAR(k)))
[716]	203	else
	204	p_cont = dble(PMID(k)scalepgfrac(igas)*(1.-QVAR(k)))
	205	endif
[305]	206
[716]	207	dtemp=0.0
	208	if(igas.eq.igas_N2)then
[253]	209
[878]	210	interm = indv(nw,igas,igas)
	211	! call interpolateN2N2(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
	212	indv(nw,igas,igas) = interm
[716]	213	! only goes to 500 cm^-1, so unless we're around a cold brown dwarf, this is irrelevant in the visible
[253]	214
[716]	215	elseif(igas.eq.igas_H2)then
[253]	216
[716]	217	! first do self-induced absorption
[878]	218	interm = indv(nw,igas,igas)
[873]	219	call interpolateH2H2(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
[878]	220	indv(nw,igas,igas) = interm
[253]	221
[716]	222	! then cross-interactions with other gases
	223	do jgas=1,ngasmx
[3233]	224	if(varspec) then
	225	p_cross = dble(PMID(k)scalepFRACVAR(jgas,k)*(1.-QVAR(k)))
	226	else
	227	p_cross = dble(PMID(k)scalepgfrac(jgas)*(1.-QVAR(k)))
	228	endif
[873]	229	dtempc = 0.0
	230	if(jgas.eq.igas_N2)then
[878]	231	interm = indv(nw,igas,jgas)
	232	call interpolateN2H2(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
	233	indv(nw,igas,jgas) = interm
[716]	234	! should be irrelevant in the visible
[2860]	235	elseif(jgas.eq.igas_CO2)then
	236	interm = indv(nw,igas,jgas)
	237	call interpolateCO2H2(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
	238	indv(nw,igas,jgas) = interm
	239	! might not be relevant in the visible
[716]	240	elseif(jgas.eq.igas_He)then
[878]	241	interm = indv(nw,igas,jgas)
[873]	242	call interpolateH2He(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
[878]	243	indv(nw,igas,jgas) = interm
[716]	244	endif
[873]	245	dtemp = dtemp + dtempc
[716]	246	enddo
[2655]	247
	248	elseif(igas.eq.igas_CH4)then
[253]	249
[2655]	250	! first do self-induced absorption
	251	interm = indv(nw,igas,igas)
	252	call interpolateCH4CH4(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
	253	indv(nw,igas,igas) = interm
	254
	255	! then cross-interactions with other gases
	256	do jgas=1,ngasmx
[3233]	257	if(varspec) then
	258	p_cross = dble(PMID(k)scalepFRACVAR(jgas,k)*(1.-QVAR(k)))
	259	else
	260	p_cross = dble(PMID(k)scalepgfrac(jgas)*(1.-QVAR(k)))
	261	endif
[2655]	262	dtempc = 0.0d0
	263	if(jgas.eq.igas_H2)then
	264	interm = indv(nw,igas,jgas)
	265	call interpolateH2CH4(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
	266	indv(nw,igas,jgas) = interm
[2861]	267	elseif(jgas.eq.igas_CO2)then
	268	interm = indv(nw,igas,jgas)
	269	call interpolateCO2CH4(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
	270	indv(nw,igas,jgas) = interm
	271	! might not be relevant in the visible
[2655]	272	elseif(jgas.eq.igas_He)then
	273	interm = indv(nw,igas,jgas)
	274	call interpolateHeCH4(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
	275	indv(nw,igas,jgas) = interm
	276	endif
	277	dtemp = dtemp + dtempc
	278	enddo
	279
[2520]	280	elseif(igas.eq.igas_H2O.and.T_cont.gt.100.0)then
	281	! Compute self and foreign (with air) continuum of H2O
[716]	282	p_air = dble(PMID(k)*scalep) - p_cont ! note assumes background is air!
[2520]	283	interm = indv(nw,igas,igas)
	284	call interpolateH2O_self_foreign(wn_cont,T_cont,p_cont,p_air,dtemp,.false.,interm) ! MTCKD v3.3
	285	indv(nw,igas,igas) = interm
[253]	286
[716]	287	endif
[253]	288
[716]	289	DCONT = DCONT + dtemp
[253]	290
[716]	291	enddo
[253]	292
[873]	293	DCONT = DCONT*dz(k)
	294
[716]	295	endif
[253]	296
[873]	297	do ng=1,L_NGAUSS-1
[305]	298
[873]	299	! Now compute TAUGAS
[253]	300
[873]	301	! Interpolate between water mixing ratios
	302	! WRATIO = 0.0 if the requested water amount is equal to, or outside the
	303	! the water data range
	304
	305	if(L_REFVAR.eq.1)then ! added by RW for special no variable case
[1725]	306
	307	! JVO 2017 : added tmpk because the repeated calls to gasi/v increased dramatically
	308	! the execution time of optci/v -> ~ factor 2 on the whole radiative
	309	! transfer on the tested simulations !
	310
[2543]	311	IF (corrk_recombin) THEN ! Added by JVO
	312	tmpk = GASV_RECOMB(MT(K):MT(K)+1,MP(K):MP(K)+1,1,NW,NG) ! contains the mix of recombined species
	313	ELSE
	314	tmpk = GASV(MT(K):MT(K)+1,MP(K):MP(K)+1,1,NW,NG)
	315	ENDIF
[1725]	316
	317	KCOEF(1) = tmpk(1,1) ! KCOEF(1) = GASV(MT(K),MP(K),1,NW,NG)
	318	KCOEF(2) = tmpk(1,2) ! KCOEF(2) = GASV(MT(K),MP(K)+1,1,NW,NG)
	319	KCOEF(3) = tmpk(2,2) ! KCOEF(3) = GASV(MT(K)+1,MP(K)+1,1,NW,NG)
	320	KCOEF(4) = tmpk(2,1) ! KCOEF(4) = GASV(MT(K)+1,MP(K),1,NW,NG)
	321
[716]	322	else
[873]	323
[2543]	324	IF (corrk_recombin) THEN
	325	tmpkvar = GASV_RECOMB(MT(K):MT(K)+1,MP(K):MP(K)+1,NVAR(K):NVAR(K)+1,NW,NG)
	326	ELSE
	327	tmpkvar = GASV(MT(K):MT(K)+1,MP(K):MP(K)+1,NVAR(K):NVAR(K)+1,NW,NG)
	328	ENDIF
[253]	329
[1725]	330	KCOEF(1) = tmpkvar(1,1,1) + WRATIO(K) * &
	331	( tmpkvar(1,1,2)-tmpkvar(1,1,1) )
[253]	332
[1725]	333	KCOEF(2) = tmpkvar(1,2,1) + WRATIO(K) * &
	334	( tmpkvar(1,2,2)-tmpkvar(1,2,1) )
[253]	335
[1725]	336	KCOEF(3) = tmpkvar(2,2,1) + WRATIO(K) * &
	337	( tmpkvar(2,2,2)-tmpkvar(2,2,1) )
	338
	339	KCOEF(4) = tmpkvar(2,1,1) + WRATIO(K) * &
	340	( tmpkvar(2,1,2)-tmpkvar(2,1,1) )
[873]	341
[1725]	342
[716]	343	endif
[253]	344
[873]	345	! Interpolate the gaseous k-coefficients to the requested T,P values
[253]	346
[873]	347	ANS = LKCOEF(K,1)KCOEF(1) + LKCOEF(K,2)KCOEF(2) + &
[716]	348	LKCOEF(K,3)KCOEF(3) + LKCOEF(K,4)KCOEF(4)
[253]	349
[873]	350	TAUGAS = U(k)*ANS
[253]	351
[716]	352	TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS + DCONT
[873]	353	DTAUKV(K,nw,ng) = TAUGAS &
[1715]	354	+ DRAYAER & ! DRAYAER includes all scattering contributions
[873]	355	+ DCONT ! For parameterized continuum aborption
[253]	356
[716]	357	end do
[253]	358
[873]	359	! Now fill in the "clear" part of the spectrum (NG = L_NGAUSS),
	360	! which holds continuum opacity only
[253]	361
[873]	362	NG = L_NGAUSS
[1715]	363	DTAUKV(K,nw,ng) = DRAYAER + DCONT ! Scattering + parameterized continuum absorption
[253]	364
[716]	365	end do
	366	end do
[253]	367
	368
[716]	369	!=======================================================================
	370	! Now the full treatment for the layers, where besides the opacity
	371	! we need to calculate the scattering albedo and asymmetry factors
[253]	372
[1987]	373	!JL18 It seems to be good to have aerosols in the first "radiative layer" of the gcm in the IR
	374	! but not in the visible
	375	! The tauaero is thus set to 0 in the 4 first semilayers in optcv, but not optci.
	376	! This solves random variations of the sw heating at the model top.
[873]	377	do iaer=1,naerkind
[918]	378	DO NW=1,L_NSPECTV
[1987]	379	TAUAEROLK(1:4,NW,IAER)=0.d0
	380	DO K=5,L_LEVELS
[1715]	381	TAUAEROLK(K,NW,IAER) = TAUAERO(K,IAER) * QSVAER(K,NW,IAER) ! effect of scattering albedo
[918]	382	ENDDO
	383	ENDDO
[873]	384	end do
[253]	385
[716]	386	DO NW=1,L_NSPECTV
[919]	387	DO L=1,L_NLAYRAD-1
[918]	388	K = 2*L+1
	389	atemp(L,NW) = SUM(GVAER(K,NW,1:naerkind) * TAUAEROLK(K,NW,1:naerkind))+SUM(GVAER(K+1,NW,1:naerkind) * TAUAEROLK(K+1,NW,1:naerkind))
	390	btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind)) + SUM(TAUAEROLK(K+1,NW,1:naerkind))
[1715]	391	ctemp(L,NW) = btemp(L,NW) + 0.9999*(TRAY(K,NW) + TRAY(K+1,NW)) ! JVO 2017 : does this 0.999 is really meaningful ?
[918]	392	btemp(L,NW) = btemp(L,NW) + TRAY(K,NW) + TRAY(K+1,NW)
	393	COSBV(L,NW,1:L_NGAUSS) = atemp(L,NW)/btemp(L,NW)
	394	END DO ! L vertical loop
[919]	395
[1715]	396	! Last level
	397	L = L_NLAYRAD
	398	K = 2*L+1
	399	atemp(L,NW) = SUM(GVAER(K,NW,1:naerkind) * TAUAEROLK(K,NW,1:naerkind))
[919]	400	btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind))
[1715]	401	ctemp(L,NW) = btemp(L,NW) + 0.9999*TRAY(K,NW) ! JVO 2017 : does this 0.999 is really meaningful ?
[919]	402	btemp(L,NW) = btemp(L,NW) + TRAY(K,NW)
	403	COSBV(L,NW,1:L_NGAUSS) = atemp(L,NW)/btemp(L,NW)
	404
	405
[918]	406	END DO ! NW spectral loop
	407
	408	DO NG=1,L_NGAUSS
	409	DO NW=1,L_NSPECTV
[873]	410	DO L=1,L_NLAYRAD-1
[253]	411
[873]	412	K = 2*L+1
	413	DTAUV(L,nw,ng) = DTAUKV(K,NW,NG) + DTAUKV(K+1,NW,NG)
[918]	414	WBARV(L,nw,ng) = ctemp(L,NW) / DTAUV(L,nw,ng)
[253]	415
[873]	416	END DO ! L vertical loop
[253]	417
[1715]	418	! Last level
[253]	419
[716]	420	L = L_NLAYRAD
	421	K = 2*L+1
[919]	422	DTAUV(L,nw,ng) = DTAUKV(K,NW,NG)
	423
	424	WBARV(L,NW,NG) = ctemp(L,NW) / DTAUV(L,NW,NG)
[1722]	425
[918]	426	END DO ! NW spectral loop
	427	END DO ! NG Gauss loop
[716]	428
	429	! Total extinction optical depths
	430
[918]	431	DO NG=1,L_NGAUSS ! full gauss loop
	432	DO NW=1,L_NSPECTV
[716]	433	TAUCUMV(1,NW,NG)=0.0D0
	434	DO K=2,L_LEVELS
	435	TAUCUMV(K,NW,NG)=TAUCUMV(K-1,NW,NG)+DTAUKV(K,NW,NG)
	436	END DO
[1987]	437
[2004]	438	DO L=1,L_NLAYRAD
[1987]	439	TAUV(L,NW,NG)=TAUCUMV(2*L,NW,NG)
	440	END DO
[2004]	441	TAUV(L,NW,NG)=TAUCUMV(2*L_NLAYRAD+1,NW,NG)
[918]	442	END DO
	443	END DO ! end full gauss loop
[716]	444
	445
[2131]	446
	447
[2032]	448	end subroutine optcv
[873]	449
[2032]	450	END MODULE optcv_mod

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

Download in other formats:

Original Format