Context Navigation

optcv.F90 @ 3641

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

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

Download in other formats:

Original Format