Context Navigation

optcv.F90 @ 3696

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

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

Download in other formats:

Original Format