Context Navigation

callcorrk.F90 @ 2054

Last change on this file since 2054 was 2050, checked in by jvatant, 7 years ago

Major radiative transfer contribution : Add the 'corrk_recombin' option that allows to use
correlated-k for single species instead of pre-mix and enables more flexiblity for variable species.
-> Algorithm inspired from Lacis and Oinas 1991 and Amundsen et al 2016

+ Added 'recombin_corrk.F90' - Important improvements in 'sugas_corrk.90' and 'callcork.F90'

Must have the desired variable species as tracers -> TBD : Enable to force composition even if no tracers

To have decent CPU time recombining is not done on all gridpoints and wavelenghts but we calculate a gasi/v_recomb variable on the reference corrk-k T,P grid (only for T,P values who match the atmospheric conditions ) which is then processed as a standard pre-mix in optci/v routines, but updated every time tracers on the ref P grid have varied > 1%.

READ CAREFULY :

In case of 'corrk_recombin', the variable L_NREFVAR doesn't have the same meaning as before and doesn't stand for the different mixing ratios but the different species.
Input corr-k should be found in corrkdir within 'corrk_gcm_IR/VI_XXX.dat' and can contain a 'fixed' specie ( compulsory if you include self-broadening ) that MUST have been created with correct mixing ratios, or a variable specie for which mixing ratio MUST have been set to 1 ( no self-broadening then, assume it's a trace sepecie ) -> You can't neither have CIA of variable species included upstream in the corr-k

Property svn:executable set to *

File size: 25.4 KB

Rev	Line
[2046]	1	subroutine callcorrk(ngrid,nlayer,pq,nq,qsurf,zday, &
[1482]	2	albedo,albedo_equivalent,emis,mu0,pplev,pplay,pt, &
[1788]	3	tsurf,fract,dist_star, &
[253]	4	dtlw,dtsw,fluxsurf_lw, &
[1482]	5	fluxsurf_sw,fluxsurfabs_sw,fluxtop_lw, &
	6	fluxabs_sw,fluxtop_dn, &
[538]	7	OLR_nu,OSR_nu, &
[1822]	8	lastcall)
[253]	9
[1722]	10	use mod_phys_lmdz_para, only : is_master
[253]	11	use radinc_h
	12	use radcommon_h
[471]	13	use gases_h
[787]	14	USE tracer_h
[1822]	15	use callkeys_mod, only: global1d, szangle
[1384]	16	use comcstfi_mod, only: pi, mugaz, cpp
[2050]	17	use callkeys_mod, only: diurnal,tracer,seashaze,corrk_recombin, &
[1647]	18	strictboundcorrk,specOLR
[2046]	19	use geometry_mod, only: latitude
[787]	20
[253]	21	implicit none
	22
	23	!==================================================================
	24	!
	25	! Purpose
	26	! -------
	27	! Solve the radiative transfer using the correlated-k method for
	28	! the gaseous absorption and the Toon et al. (1989) method for
	29	! scatttering due to aerosols.
	30	!
	31	! Authors
	32	! -------
	33	! Emmanuel 01/2001, Forget 09/2001
	34	! Robin Wordsworth (2009)
[2050]	35	! Jan Vatant d'Ollone (2018) -> corrk recombining case
[253]	36	!
	37	!==================================================================
	38
	39	!-----------------------------------------------------------------------
	40	! Declaration of the arguments (INPUT - OUTPUT) on the LMD GCM grid
	41	! Layer #1 is the layer near the ground.
[1308]	42	! Layer #nlayer is the layer at the top.
[1483]	43	!-----------------------------------------------------------------------
[253]	44
[1483]	45
	46	! INPUT
	47	INTEGER,INTENT(IN) :: ngrid ! Number of atmospheric columns.
	48	INTEGER,INTENT(IN) :: nlayer ! Number of atmospheric layers.
[2050]	49	REAL,INTENT(IN) :: pq(ngrid,nlayer,nq) ! Tracers (X/kg).
[1483]	50	INTEGER,INTENT(IN) :: nq ! Number of tracers.
	51	REAL,INTENT(IN) :: qsurf(ngrid,nq) ! Tracers on surface (kg.m-2).
[2046]	52	REAL,INTENT(IN) :: zday ! Time elapsed since Ls=0 (sols).
[1483]	53	REAL,INTENT(IN) :: albedo(ngrid,L_NSPECTV) ! Spectral Short Wavelengths Albedo. By MT2015
	54	REAL,INTENT(IN) :: emis(ngrid) ! Long Wave emissivity.
	55	REAL,INTENT(IN) :: mu0(ngrid) ! Cosine of sun incident angle.
	56	REAL,INTENT(IN) :: pplev(ngrid,nlayer+1) ! Inter-layer pressure (Pa).
	57	REAL,INTENT(IN) :: pplay(ngrid,nlayer) ! Mid-layer pressure (Pa).
	58	REAL,INTENT(IN) :: pt(ngrid,nlayer) ! Air temperature (K).
	59	REAL,INTENT(IN) :: tsurf(ngrid) ! Surface temperature (K).
	60	REAL,INTENT(IN) :: fract(ngrid) ! Fraction of day.
	61	REAL,INTENT(IN) :: dist_star ! Distance star-planet (AU).
	62	logical,intent(in) :: lastcall ! Signals last call to physics.
	63
	64	! OUTPUT
	65	REAL,INTENT(OUT) :: dtlw(ngrid,nlayer) ! Heating rate (K/s) due to LW radiation.
	66	REAL,INTENT(OUT) :: dtsw(ngrid,nlayer) ! Heating rate (K/s) due to SW radiation.
	67	REAL,INTENT(OUT) :: fluxsurf_lw(ngrid) ! Incident LW flux to surf (W/m2).
	68	REAL,INTENT(OUT) :: fluxsurf_sw(ngrid) ! Incident SW flux to surf (W/m2)
	69	REAL,INTENT(OUT) :: fluxsurfabs_sw(ngrid) ! Absorbed SW flux by the surface (W/m2). By MT2015.
	70	REAL,INTENT(OUT) :: fluxtop_lw(ngrid) ! Outgoing LW flux to space (W/m2).
	71	REAL,INTENT(OUT) :: fluxabs_sw(ngrid) ! SW flux absorbed by the planet (W/m2).
	72	REAL,INTENT(OUT) :: fluxtop_dn(ngrid) ! Incident top of atmosphere SW flux (W/m2).
	73	REAL,INTENT(OUT) :: OLR_nu(ngrid,L_NSPECTI) ! Outgoing LW radition in each band (Normalized to the band width (W/m2/cm-1).
	74	REAL,INTENT(OUT) :: OSR_nu(ngrid,L_NSPECTV) ! Outgoing SW radition in each band (Normalized to the band width (W/m2/cm-1).
[1722]	75	REAL,INTENT(OUT) :: albedo_equivalent(ngrid) ! Spectrally Integrated Albedo. For Diagnostic. By MT2015
[1483]	76
	77
[253]	78	!-----------------------------------------------------------------------
	79	! Declaration of the variables required by correlated-k subroutines
[1483]	80	! Numbered from top to bottom (unlike in the GCM)
	81	!-----------------------------------------------------------------------
[253]	82
	83	REAL*8 tmid(L_LEVELS),pmid(L_LEVELS)
	84	REAL*8 tlevrad(L_LEVELS),plevrad(L_LEVELS)
	85
[1483]	86	! Optical values for the optci/cv subroutines
[253]	87	REAL*8 stel(L_NSPECTV),stel_fract(L_NSPECTV)
	88	REAL*8 dtaui(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
	89	REAL*8 dtauv(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
	90	REAL*8 cosbv(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
	91	REAL*8 cosbi(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
	92	REAL*8 wbari(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
	93	REAL*8 wbarv(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
	94	REAL*8 tauv(L_NLEVRAD,L_NSPECTV,L_NGAUSS)
	95	REAL*8 taucumv(L_LEVELS,L_NSPECTV,L_NGAUSS)
	96	REAL*8 taucumi(L_LEVELS,L_NSPECTI,L_NGAUSS)
	97
[961]	98	REAL*8 nfluxtopv,nfluxtopi,nfluxtop,fluxtopvdn
[1483]	99	REAL*8 nfluxoutv_nu(L_NSPECTV) ! Outgoing band-resolved VI flux at TOA (W/m2).
	100	REAL*8 nfluxtopi_nu(L_NSPECTI) ! Net band-resolved IR flux at TOA (W/m2).
	101	REAL*8 fluxupi_nu(L_NLAYRAD,L_NSPECTI) ! For 1D diagnostic.
[253]	102	REAL*8 fmneti(L_NLAYRAD),fmnetv(L_NLAYRAD)
	103	REAL*8 fluxupv(L_NLAYRAD),fluxupi(L_NLAYRAD)
	104	REAL*8 fluxdnv(L_NLAYRAD),fluxdni(L_NLAYRAD)
[1482]	105	REAL*8 albi,acosz
[1483]	106	REAL*8 albv(L_NSPECTV) ! Spectral Visible Albedo.
[253]	107
[2050]	108	INTEGER ig,l,k,nw,iq,ip,ilay,it
	109
	110	LOGICAL found
[253]	111
	112	real*8 taugsurf(L_NSPECTV,L_NGAUSS-1)
	113	real*8 taugsurfi(L_NSPECTI,L_NGAUSS-1)
	114
	115	logical OLRz
	116	real*8 NFLUXGNDV_nu(L_NSPECTV)
[1648]	117
[1483]	118	! Included by MT for albedo calculations.
[1482]	119	REAL*8 albedo_temp(L_NSPECTV) ! For equivalent albedo calculation.
[1526]	120	REAL*8 surface_stellar_flux ! Stellar flux reaching the surface. Useful for equivalent albedo calculation.
[2046]	121
	122	! For variable haze
	123	REAL*8 seashazefact(L_LEVELS)
[305]	124
[2050]	125	! For muphys optics
	126	REAL*8 pqmo(ngrid,nlayer,nmicro) ! Tracers for microphysics optics (X/m2).
	127	REAL*8 i2e(ngrid,nlayer) ! int 2 ext factor ( X.kg-1 -> X.m-2 for optics )
	128
	129	! For corr-k recombining
	130	REAL*8 pqr(ngrid,L_PINT,L_REFVAR) ! Tracers for corr-k recombining (mol/mol).
	131	REAL*8 fact, tmin, tmax
	132
	133	LOGICAL usept(L_PINT,L_NTREF) ! mask if pfref grid point will be used
	134	INTEGER inflay(L_PINT) ! nearest inferior GCM layer for pfgasref grid points
	135
[253]	136	!=======================================================================
[1822]	137	! I. Initialization on every call
[1483]	138	!=======================================================================
	139
[1529]	140
[1483]	141	! How much light do we get ?
[253]	142	do nw=1,L_NSPECTV
	143	stel(nw)=stellarf(nw)/(dist_star**2)
	144	end do
	145
[2050]	146	! Convert (microphysical) tracers for optics: X.kg-1 --> X.m-2
	147	! NOTE: it should be moved somewhere else: calmufi performs the same kind of
	148	! computations... waste of time...
	149	i2e(:,1:nlayer) = ( pplev(:,1:nlayer)-pplev(:,2:nlayer+1) ) / gzlat(:,1:nlayer)
	150	pqmo(:,:,:) = 0.0
	151	DO iq=1,nmicro
	152	pqmo(:,:,iq) = pq(:,:,iq)*i2e(:,:)
	153	ENDDO
[253]	154
[2050]	155	! Default value for fixed species for whom vmr has been taken
	156	! into account while computing high-resolution spectra
	157	if (corrk_recombin) pqr(:,:,:) = 1.0
	158
[1483]	159	!-----------------------------------------------------------------------
	160	do ig=1,ngrid ! Starting Big Loop over every GCM column
[253]	161	!-----------------------------------------------------------------------
[2050]	162
	163	! Recombine reference corr-k if needed
	164	if (corrk_recombin) then
	165
	166	! NB : To have decent CPU time recombining is not done on all gridpoints and wavelenghts but we
	167	! calculate a gasi/v_recomb variable on the reference corrk-k T,P grid (only for T,P values
	168	! who match the atmospheric conditions ) which is then processed as a standard pre-mix in
	169	! optci/v routines, but updated every time tracers on the ref P grid have varied > 1%.
[253]	170
[2050]	171	! Extract tracers for variable radiative species
	172	! Also find the nearest GCM layer under each ref pressure
	173	do ip=1,L_PINT
	174
	175	ilay=0
	176	found = .false.
	177	do l=1,nlayer
	178	if ( pplay(ig,l) .gt. 10.0**(pfgasref(ip)+2.0) ) then ! pfgasref=log(p[mbar])
	179	found=.true.
	180	ilay=l
	181	endif
	182	enddo
[1483]	183
[2050]	184	if (.not. found ) then ! set to min
	185	do iq=1,L_REFVAR
	186	if ( radvar_mask(iq) ) then
	187	pqr(ig,ip,iq) = pq(ig,1,radvar_indx(iq)) / rat_mmol(radvar_indx(iq)-nmicro) ! mol/mol
	188	endif
	189	enddo
	190	else
	191	if (ilay==nlayer) then ! set to max
	192	do iq=1,L_REFVAR
	193	if ( radvar_mask(iq) ) then
	194	pqr(ig,ip,iq) = pq(ig,nlayer,radvar_indx(iq)) / rat_mmol(radvar_indx(iq)-nmicro) ! mol/mol
	195	endif
	196	enddo
	197	else ! standard
	198	fact = ( 10.0**(pfgasref(ip)+2.0) - pplay(1,ilay+1) ) / ( pplay(1,ilay) - pplay(1,ilay+1) ) ! pfgasref=log(p[mbar])
	199	do iq=1,L_REFVAR
	200	if ( radvar_mask(iq) ) then
	201	pqr(ig,ip,iq) = pq(ig,ilay,radvar_indx(iq))*fact pq(ig,ilay+1,radvar_indx(iq))**(1.0-fact)
	202	pqr(ig,ip,iq) = pqr(ig,ip,iq) / rat_mmol(radvar_indx(iq)-nmicro) ! mol/mol
	203	endif
	204	enddo
	205	endif ! if ilay==nlayer
	206	endif ! if not found
	207
	208	inflay(ip) = ilay
	209
	210	enddo ! ip=1,L_PINT
	211
	212	! NB : The following usept is a trick to call recombine only for the reference T-P
	213	! grid points that are useful given the temperature range at this altitude
	214	! It saves a looot of time - JVO 18
	215	usept(:,:) = .true.
	216	do ip=1,L_PINT-1
	217	if ( inflay(ip+1)==nlayer ) then
	218	usept(ip,:) = .false.
	219	endif
	220	if ( inflay(ip) == 0 ) then
	221	usept(ip+1:,:) = .false.
	222	endif
	223	if ( usept(ip,1) ) then ! if not all false
	224	tmin = minval(pt(ig,min(inflay(ip+1)+1,nlayer):max(inflay(ip),1)))
	225	tmax = maxval(pt(ig,min(inflay(ip+1)+1,nlayer):max(inflay(ip),1)))
	226	do it=1,L_NTREF-1
	227	if ( tgasref(it+1) .lt. tmin ) then
	228	usept(ip,it) = .false.
	229	endif
	230	enddo
	231	do it=2,L_NTREF
	232	if ( tgasref(it-1) .gt. tmax ) then
	233	usept(ip,it) = .false.
	234	endif
	235	enddo
	236	! in case of out-of-bounds
	237	if ( tgasref(1) .lt. tmin ) usept(ip,1) = .true.
	238	if ( tgasref(L_NTREF) .gt. tmax ) usept(ip,L_NTREF) = .true.
	239	endif
	240	enddo ! ip=1,L_PINT-1
	241	! deal with last bound
	242	if ( inflay(L_PINT-1).ne.0 ) usept(L_PINT,:) = usept(L_PINT-1,:)
	243
	244
	245	do ip=1,L_PINT
	246
	247	! Recombine k at (useful only!) reference T-P values if tracers or T have enough varied
	248	do it=1,L_NTREF
	249
	250	if ( usept(ip,it) .eqv. .false. ) cycle
	251
	252	do l=1,L_REFVAR
	253	if ( abs( (pqr(ig,ip,l) - pqrold(ip,l)) / max(1.0e-30,pqrold(ip,l))) .GT. 0.01 & ! +- 1%
	254	.or. ( useptold(ip,it) .eqv. .false. ) ) then ! in case T change but not the tracers
	255	call recombin_corrk( pqr(ig,ip,:),ip,it )
	256	exit ! one is enough to trigger the update
	257	endif
	258	enddo
	259
	260	enddo
	261
	262	enddo ! ip=1,L_PINT
	263
	264	useptold(:,:)=usept(:,:)
	265
	266	endif ! if corrk_recombin
	267
[253]	268	!=======================================================================
[1483]	269	! II. Transformation of the GCM variables
	270	!=======================================================================
[253]	271
[1483]	272
	273	! Albedo and Emissivity.
	274	albi=1-emis(ig) ! Long Wave.
	275	DO nw=1,L_NSPECTV ! Short Wave loop.
[1482]	276	albv(nw)=albedo(ig,nw)
[1529]	277	ENDDO
[253]	278
[1483]	279	if ((ngrid.eq.1).and.(global1d)) then ! Fixed zenith angle 'szangle' in 1D simulations w/ globally-averaged sunlight.
[253]	280	acosz = cos(pi*szangle/180.0)
	281	print*,'acosz=',acosz,', szangle=',szangle
	282	else
[1483]	283	acosz=mu0(ig) ! Cosine of sun incident angle : 3D simulations or local 1D simulations using latitude.
[253]	284	endif
	285
	286	!-----------------------------------------------------------------------
	287	! Pressure and temperature
[1483]	288	!-----------------------------------------------------------------------
[253]	289
	290	DO l=1,nlayer
	291	plevrad(2*l) = pplay(ig,nlayer+1-l)/scalep
	292	plevrad(2*l+1) = pplev(ig,nlayer+1-l)/scalep
	293	tlevrad(2*l) = pt(ig,nlayer+1-l)
	294	tlevrad(2*l+1) = (pt(ig,nlayer+1-l)+pt(ig,max(nlayer-l,1)))/2
	295	END DO
	296
[600]	297	plevrad(1) = 0.
[1483]	298	plevrad(2) = 0. !! Trick to have correct calculations of fluxes in gflux(i/v).F, but the pmid levels are not impacted by this change.
[253]	299
	300	tlevrad(1) = tlevrad(2)
[1308]	301	tlevrad(2*nlayer+1)=tsurf(ig)
[253]	302
[1423]	303	pmid(1) = max(pgasmin,0.0001*plevrad(3))
	304	pmid(2) = pmid(1)
	305
[253]	306	tmid(1) = tlevrad(2)
[1423]	307	tmid(2) = tmid(1)
	308
	309	DO l=1,L_NLAYRAD-1
	310	tmid(2l+1) = tlevrad(2l+1)
	311	tmid(2l+2) = tlevrad(2l+1)
	312	pmid(2l+1) = plevrad(2l+1)
	313	pmid(2l+2) = plevrad(2l+1)
[253]	314	END DO
[1423]	315	pmid(L_LEVELS) = plevrad(L_LEVELS)
	316	tmid(L_LEVELS) = tlevrad(L_LEVELS)
[253]	317
[1423]	318	!!Alternative interpolation:
	319	! pmid(3) = pmid(1)
	320	! pmid(4) = pmid(1)
	321	! tmid(3) = tmid(1)
	322	! tmid(4) = tmid(1)
	323	! DO l=2,L_NLAYRAD-1
	324	! tmid(2l+1) = tlevrad(2l)
	325	! tmid(2l+2) = tlevrad(2l)
	326	! pmid(2l+1) = plevrad(2l)
	327	! pmid(2l+2) = plevrad(2l)
	328	! END DO
	329	! pmid(L_LEVELS) = plevrad(L_LEVELS-1)
	330	! tmid(L_LEVELS) = tlevrad(L_LEVELS-1)
	331
[1483]	332	! Test for out-of-bounds pressure.
[253]	333	if(plevrad(3).lt.pgasmin)then
	334	print*,'Minimum pressure is outside the radiative'
	335	print*,'transfer kmatrix bounds, exiting.'
	336	call abort
	337	elseif(plevrad(L_LEVELS).gt.pgasmax)then
	338	print*,'Maximum pressure is outside the radiative'
	339	print*,'transfer kmatrix bounds, exiting.'
	340	call abort
	341	endif
	342
[1483]	343	! Test for out-of-bounds temperature.
[253]	344	do k=1,L_LEVELS
	345	if(tlevrad(k).lt.tgasmin)then
	346	print*,'Minimum temperature is outside the radiative'
[1145]	347	print*,'transfer kmatrix bounds'
[858]	348	print*,"k=",k," tlevrad(k)=",tlevrad(k)
	349	print*,"tgasmin=",tgasmin
[1145]	350	if (strictboundcorrk) then
	351	call abort
	352	else
	353	print,'**********************************************'
	354	print*,'we allow model to continue with tlevrad=tgasmin'
	355	print*,' ... we assume we know what you are doing ... '
	356	print*,' ... but do not let this happen too often ... '
	357	print,'**********************************************'
[1147]	358	!tlevrad(k)=tgasmin
[1145]	359	endif
[253]	360	elseif(tlevrad(k).gt.tgasmax)then
[1648]	361	! print*,'Maximum temperature is outside the radiative'
	362	! print*,'transfer kmatrix bounds, exiting.'
	363	! print*,"k=",k," tlevrad(k)=",tlevrad(k)
	364	! print*,"tgasmax=",tgasmax
[1145]	365	if (strictboundcorrk) then
	366	call abort
	367	else
[1648]	368	! print,'**********************************************'
	369	! print*,'we allow model to continue with tlevrad=tgasmax'
	370	! print*,' ... we assume we know what you are doing ... '
	371	! print*,' ... but do not let this happen too often ... '
	372	! print,'**********************************************'
[1147]	373	!tlevrad(k)=tgasmax
[1145]	374	endif
[253]	375	endif
	376	enddo
[1016]	377	do k=1,L_NLAYRAD+1
	378	if(tmid(k).lt.tgasmin)then
	379	print*,'Minimum temperature is outside the radiative'
	380	print*,'transfer kmatrix bounds, exiting.'
[1145]	381	print*,"k=",k," tmid(k)=",tmid(k)
[1016]	382	print*,"tgasmin=",tgasmin
[1145]	383	if (strictboundcorrk) then
	384	call abort
	385	else
	386	print,'**********************************************'
	387	print*,'we allow model to continue with tmid=tgasmin'
	388	print*,' ... we assume we know what you are doing ... '
	389	print*,' ... but do not let this happen too often ... '
	390	print,'**********************************************'
	391	tmid(k)=tgasmin
	392	endif
[1016]	393	elseif(tmid(k).gt.tgasmax)then
[1648]	394	! print*,'Maximum temperature is outside the radiative'
	395	! print*,'transfer kmatrix bounds, exiting.'
	396	! print*,"k=",k," tmid(k)=",tmid(k)
	397	! print*,"tgasmax=",tgasmax
[1145]	398	if (strictboundcorrk) then
	399	call abort
	400	else
[1648]	401	! print,'**********************************************'
	402	! print*,'we allow model to continue with tmid=tgasmin'
	403	! print*,' ... we assume we know what you are doing ... '
	404	! print*,' ... but do not let this happen too often ... '
	405	! print,'**********************************************'
[1145]	406	tmid(k)=tgasmax
	407	endif
[1016]	408	endif
	409	enddo
[253]	410
	411	!=======================================================================
[1483]	412	! III. Calling the main radiative transfer subroutines
	413	!=======================================================================
[253]	414
[1947]	415	Cmk(:) = 0.01 * 1.0 / (gzlat(ig,:) * mugaz * 1.672621e-27) ! q_main=1.0 assumed.
	416	gzlat_ig(:) = gzlat(ig,:)
[2046]	417
	418	! Compute the haze seasonal modulation factor
	419	if (seashaze) call season_haze(zday,latitude(ig),plevrad,seashazefact)
[253]	420
	421	!-----------------------------------------------------------------------
[1483]	422	! Short Wave Part
	423	!-----------------------------------------------------------------------
[253]	424
[1483]	425	if(fract(ig) .ge. 1.0e-4) then ! Only during daylight.
[787]	426	if((ngrid.eq.1).and.(global1d))then
[253]	427	do nw=1,L_NSPECTV
[1483]	428	stel_fract(nw)= stel(nw)* 0.25 / acosz ! globally averaged = divide by 4, and we correct for solar zenith angle
[253]	429	end do
	430	else
	431	do nw=1,L_NSPECTV
[1161]	432	stel_fract(nw)= stel(nw) * fract(ig)
[253]	433	end do
[1483]	434	endif
	435
[2050]	436	call optcv(pqmo(ig,:,:),nlayer,plevrad,tmid,pmid, &
[2046]	437	dtauv,tauv,taucumv,wbarv,cosbv,tauray,taugsurf,seashazefact)
[253]	438
	439	call sfluxv(dtauv,tauv,taucumv,albv,dwnv,wbarv,cosbv, &
[1781]	440	acosz,stel_fract, &
	441	nfluxtopv,fluxtopvdn,nfluxoutv_nu,nfluxgndv_nu, &
[253]	442	fmnetv,fluxupv,fluxdnv,fzerov,taugsurf)
	443
[1483]	444	else ! During the night, fluxes = 0.
[962]	445	nfluxtopv = 0.0d0
[1529]	446	fluxtopvdn = 0.0d0
[962]	447	nfluxoutv_nu(:) = 0.0d0
	448	nfluxgndv_nu(:) = 0.0d0
[253]	449	do l=1,L_NLAYRAD
[962]	450	fmnetv(l)=0.0d0
	451	fluxupv(l)=0.0d0
	452	fluxdnv(l)=0.0d0
[253]	453	end do
	454	end if
	455
[1482]	456
[1526]	457	! Equivalent Albedo Calculation (for OUTPUT). MT2015
	458	if(fract(ig) .ge. 1.0e-4) then ! equivalent albedo makes sense only during daylight.
	459	surface_stellar_flux=sum(nfluxgndv_nu(1:L_NSPECTV))
	460	if(surface_stellar_flux .gt. 1.0e-3) then ! equivalent albedo makes sense only if the stellar flux received by the surface is positive.
[1529]	461	DO nw=1,L_NSPECTV
	462	albedo_temp(nw)=albedo(ig,nw)*nfluxgndv_nu(nw)
[1526]	463	ENDDO
[1529]	464	albedo_temp(1:L_NSPECTV)=albedo_temp(1:L_NSPECTV)/surface_stellar_flux
[1526]	465	albedo_equivalent(ig)=sum(albedo_temp(1:L_NSPECTV))
	466	else
	467	albedo_equivalent(ig)=0.0 ! Spectrally Integrated Albedo not defined for non-irradiated grid points. So we arbitrary set the equivalent albedo to 0.
	468	endif
[1529]	469	else
	470	albedo_equivalent(ig)=0.0 ! Spectrally Integrated Albedo not defined for non-irradiated grid points. So we arbitrary set the equivalent albedo to 0.
	471	endif
[1482]	472
	473
[253]	474	!-----------------------------------------------------------------------
[1483]	475	! Long Wave Part
	476	!-----------------------------------------------------------------------
[253]	477
[2050]	478	call optci(pqmo(ig,:,:),nlayer,plevrad,tlevrad,tmid,pmid, &
[2046]	479	dtaui,taucumi,cosbi,wbari,taugsurfi,seashazefact)
[538]	480
[253]	481	call sfluxi(plevrad,tlevrad,dtaui,taucumi,ubari,albi, &
[1781]	482	wnoi,dwni,cosbi,wbari,nfluxtopi,nfluxtopi_nu, &
[253]	483	fmneti,fluxupi,fluxdni,fluxupi_nu,fzeroi,taugsurfi)
	484
	485	!-----------------------------------------------------------------------
	486	! Transformation of the correlated-k code outputs
	487	! (into dtlw, dtsw, fluxsurf_lw, fluxsurf_sw, fluxtop_lw, fluxtop_sw)
	488
	489	! Flux incident at the top of the atmosphere
[961]	490	fluxtop_dn(ig)=fluxtopvdn
[253]	491
	492	fluxtop_lw(ig) = real(nfluxtopi)
	493	fluxabs_sw(ig) = real(-nfluxtopv)
	494	fluxsurf_lw(ig) = real(fluxdni(L_NLAYRAD))
	495	fluxsurf_sw(ig) = real(fluxdnv(L_NLAYRAD))
[1482]	496
	497	! Flux absorbed by the surface. By MT2015.
	498	fluxsurfabs_sw(ig) = fluxsurf_sw(ig)*(1.-albedo_equivalent(ig))
[253]	499
	500	if(fluxtop_dn(ig).lt.0.0)then
	501	print*,'Achtung! fluxtop_dn has lost the plot!'
	502	print*,'fluxtop_dn=',fluxtop_dn(ig)
	503	print*,'acosz=',acosz
	504	print*,'temp= ',pt(ig,:)
	505	print*,'pplay= ',pplay(ig,:)
	506	call abort
	507	endif
	508
	509	! Spectral output, for exoplanet observational comparison
	510	if(specOLR)then
	511	do nw=1,L_NSPECTI
[526]	512	OLR_nu(ig,nw)=nfluxtopi_nu(nw)/DWNI(nw) !JL Normalize to the bandwidth
[253]	513	end do
	514	do nw=1,L_NSPECTV
[366]	515	!GSR_nu(ig,nw)=nfluxgndv_nu(nw)
[526]	516	OSR_nu(ig,nw)=nfluxoutv_nu(nw)/DWNV(nw) !JL Normalize to the bandwidth
[253]	517	end do
	518	endif
	519
	520	! Finally, the heating rates
	521
[586]	522	DO l=2,L_NLAYRAD
	523	dtsw(ig,L_NLAYRAD+1-l)=(fmnetv(l)-fmnetv(l-1)) &
[1947]	524	gzlat(ig,L_NLAYRAD+1-l)/(cppscalep(plevrad(2l+1)-plevrad(2*l-1)))
[586]	525	dtlw(ig,L_NLAYRAD+1-l)=(fmneti(l)-fmneti(l-1)) &
[1947]	526	gzlat(ig,L_NLAYRAD+1-l)/(cppscalep(plevrad(2l+1)-plevrad(2*l-1)))
[586]	527	END DO
[253]	528
	529	! These are values at top of atmosphere
[586]	530	dtsw(ig,L_NLAYRAD)=(fmnetv(1)-nfluxtopv) &
[1947]	531	gzlat(ig,L_NLAYRAD)/(cppscalep*(plevrad(3)-plevrad(1)))
[586]	532	dtlw(ig,L_NLAYRAD)=(fmneti(1)-nfluxtopi) &
[1947]	533	gzlat(ig,L_NLAYRAD)/(cppscalep*(plevrad(3)-plevrad(1)))
[253]	534
	535
[1483]	536	!-----------------------------------------------------------------------
	537	end do ! End of big loop over every GCM column.
	538	!-----------------------------------------------------------------------
[253]	539
[1483]	540
	541
[253]	542	!-----------------------------------------------------------------------
	543	! Additional diagnostics
[1483]	544	!-----------------------------------------------------------------------
[253]	545
[1483]	546	! IR spectral output, for exoplanet observational comparison
	547	if(lastcall.and.(ngrid.eq.1))then ! could disable the 1D output, they are in the diagfi and diagspec... JL12
[253]	548
[1483]	549	print*,'Saving scalar quantities in surf_vals.out...'
	550	print*,'psurf = ', pplev(1,1),' Pa'
	551	open(116,file='surf_vals.out')
	552	write(116,*) tsurf(1),pplev(1,1),fluxtop_dn(1), &
	553	real(-nfluxtopv),real(nfluxtopi)
	554	close(116)
[253]	555
[526]	556
[1483]	557	! USEFUL COMMENT - Do Not Remove.
	558	!
[526]	559	! if(specOLR)then
	560	! open(117,file='OLRnu.out')
	561	! do nw=1,L_NSPECTI
	562	! write(117,*) OLR_nu(1,nw)
	563	! enddo
	564	! close(117)
	565	!
	566	! open(127,file='OSRnu.out')
	567	! do nw=1,L_NSPECTV
	568	! write(127,*) OSR_nu(1,nw)
	569	! enddo
	570	! close(127)
	571	! endif
[253]	572
[1483]	573	! OLR vs altitude: do it as a .txt file.
	574	OLRz=.false.
	575	if(OLRz)then
	576	print*,'saving IR vertical flux for OLRz...'
	577	open(118,file='OLRz_plevs.out')
	578	open(119,file='OLRz.out')
	579	do l=1,L_NLAYRAD
	580	write(118,) plevrad(2l)
	581	do nw=1,L_NSPECTI
	582	write(119,*) fluxupi_nu(l,nw)
	583	enddo
	584	enddo
	585	close(118)
	586	close(119)
	587	endif
[253]	588
[305]	589	endif
[253]	590
[1647]	591	if (lastcall) then
[470]	592	IF( ALLOCATED( gasi ) ) DEALLOCATE( gasi )
	593	IF( ALLOCATED( gasv ) ) DEALLOCATE( gasv )
[2050]	594	IF( ALLOCATED( gasi_recomb ) ) DEALLOCATE( gasi_recomb )
	595	IF( ALLOCATED( gasv_recomb ) ) DEALLOCATE( gasv_recomb )
	596	IF( ALLOCATED( pqrold ) ) DEALLOCATE( pqrold )
	597	IF( ALLOCATED( useptold ) ) DEALLOCATE( useptold )
[1315]	598	!$OMP BARRIER
	599	!$OMP MASTER
[470]	600	IF( ALLOCATED( pgasref ) ) DEALLOCATE( pgasref )
	601	IF( ALLOCATED( tgasref ) ) DEALLOCATE( tgasref )
	602	IF( ALLOCATED( pfgasref ) ) DEALLOCATE( pfgasref )
[2026]	603	IF( ALLOCATED( gweight ) ) DEALLOCATE( gweight )
[1315]	604	!$OMP END MASTER
[1529]	605	!$OMP BARRIER
[470]	606	endif
	607
[716]	608
[253]	609	end subroutine callcorrk

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

Download in other formats:

Original Format