Context Navigation

callcorrk.F90 @ 3318

Last change on this file since 3318 was 3318, checked in by slebonnois, 7 months ago
Titan PCM update : optics + microphysics
Property svn:executable set to ``*
File size: 40.0 KB

Rev	Line
[2046]	1	subroutine callcorrk(ngrid,nlayer,pq,nq,qsurf,zday, &
[3083]	2	albedo,albedo_equivalent,emis,mu0,pplev,pplay,zzlev, &
	3	pt,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, &
[3318]	8	int_dtaui,int_dtauv,popthi,popthv,poptti,popttv, &
[1822]	9	lastcall)
[253]	10
[1722]	11	use mod_phys_lmdz_para, only : is_master
[253]	12	use radinc_h
	13	use radcommon_h
[471]	14	use gases_h
[787]	15	USE tracer_h
[1822]	16	use callkeys_mod, only: global1d, szangle
[1384]	17	use comcstfi_mod, only: pi, mugaz, cpp
[2050]	18	use callkeys_mod, only: diurnal,tracer,seashaze,corrk_recombin, &
[2366]	19	strictboundcorrk,specOLR,diagdtau, &
[3083]	20	tplanckmin,tplanckmax,callclouds,Fcloudy
[2046]	21	use geometry_mod, only: latitude
[787]	22
[253]	23	implicit none
	24
	25	!==================================================================
	26	!
	27	! Purpose
	28	! -------
	29	! Solve the radiative transfer using the correlated-k method for
	30	! the gaseous absorption and the Toon et al. (1989) method for
	31	! scatttering due to aerosols.
	32	!
	33	! Authors
	34	! -------
	35	! Emmanuel 01/2001, Forget 09/2001
	36	! Robin Wordsworth (2009)
[3090]	37	!
	38	! Modified
	39	! --------
	40	! Jan Vatant d'Ollone (2018)
	41	! --> corrk recombining case
	42	! B. de Batz de Trenquelléon (2023)
	43	! --> Titan's haze and could optics
	44	! --> clear/dark columns method
	45	! --> optical diagnostics
[253]	46	!
	47	!==================================================================
	48
	49	!-----------------------------------------------------------------------
	50	! Declaration of the arguments (INPUT - OUTPUT) on the LMD GCM grid
	51	! Layer #1 is the layer near the ground.
[1308]	52	! Layer #nlayer is the layer at the top.
[1483]	53	!-----------------------------------------------------------------------
[253]	54
[1483]	55
	56	! INPUT
	57	INTEGER,INTENT(IN) :: ngrid ! Number of atmospheric columns.
	58	INTEGER,INTENT(IN) :: nlayer ! Number of atmospheric layers.
[2050]	59	REAL,INTENT(IN) :: pq(ngrid,nlayer,nq) ! Tracers (X/kg).
[1483]	60	INTEGER,INTENT(IN) :: nq ! Number of tracers.
	61	REAL,INTENT(IN) :: qsurf(ngrid,nq) ! Tracers on surface (kg.m-2).
[2046]	62	REAL,INTENT(IN) :: zday ! Time elapsed since Ls=0 (sols).
[1483]	63	REAL,INTENT(IN) :: albedo(ngrid,L_NSPECTV) ! Spectral Short Wavelengths Albedo. By MT2015
	64	REAL,INTENT(IN) :: emis(ngrid) ! Long Wave emissivity.
	65	REAL,INTENT(IN) :: mu0(ngrid) ! Cosine of sun incident angle.
	66	REAL,INTENT(IN) :: pplev(ngrid,nlayer+1) ! Inter-layer pressure (Pa).
	67	REAL,INTENT(IN) :: pplay(ngrid,nlayer) ! Mid-layer pressure (Pa).
[3083]	68	REAL,INTENT(IN) :: zzlev(ngrid,nlayer+1) ! Altitude at the atmospheric layer boundaries (ref : local surf).
[1483]	69	REAL,INTENT(IN) :: pt(ngrid,nlayer) ! Air temperature (K).
	70	REAL,INTENT(IN) :: tsurf(ngrid) ! Surface temperature (K).
	71	REAL,INTENT(IN) :: fract(ngrid) ! Fraction of day.
	72	REAL,INTENT(IN) :: dist_star ! Distance star-planet (AU).
	73	logical,intent(in) :: lastcall ! Signals last call to physics.
	74
	75	! OUTPUT
	76	REAL,INTENT(OUT) :: dtlw(ngrid,nlayer) ! Heating rate (K/s) due to LW radiation.
	77	REAL,INTENT(OUT) :: dtsw(ngrid,nlayer) ! Heating rate (K/s) due to SW radiation.
	78	REAL,INTENT(OUT) :: fluxsurf_lw(ngrid) ! Incident LW flux to surf (W/m2).
	79	REAL,INTENT(OUT) :: fluxsurf_sw(ngrid) ! Incident SW flux to surf (W/m2)
	80	REAL,INTENT(OUT) :: fluxsurfabs_sw(ngrid) ! Absorbed SW flux by the surface (W/m2). By MT2015.
	81	REAL,INTENT(OUT) :: fluxtop_lw(ngrid) ! Outgoing LW flux to space (W/m2).
	82	REAL,INTENT(OUT) :: fluxabs_sw(ngrid) ! SW flux absorbed by the planet (W/m2).
	83	REAL,INTENT(OUT) :: fluxtop_dn(ngrid) ! Incident top of atmosphere SW flux (W/m2).
	84	REAL,INTENT(OUT) :: OLR_nu(ngrid,L_NSPECTI) ! Outgoing LW radition in each band (Normalized to the band width (W/m2/cm-1).
	85	REAL,INTENT(OUT) :: OSR_nu(ngrid,L_NSPECTV) ! Outgoing SW radition in each band (Normalized to the band width (W/m2/cm-1).
[1722]	86	REAL,INTENT(OUT) :: albedo_equivalent(ngrid) ! Spectrally Integrated Albedo. For Diagnostic. By MT2015
[3083]	87	REAL,INTENT(OUT) :: int_dtaui(ngrid,nlayer,L_NSPECTI) ! IR optical thickness of layers within narrowbands for diags ().
	88	REAL,INTENT(OUT) :: int_dtauv(ngrid,nlayer,L_NSPECTV) ! VI optical thickness of layers within narrowbands for diags ().
[3090]	89	! Diagnostics
[3318]	90	REAL,INTENT(OUT) :: popthi(ngrid,nlayer,L_NSPECTI,8) ! IR optical properties of haze within narrowbands (dtau,tau,k,w,g,drayaer,taugaz,dcont).
	91	REAL,INTENT(OUT) :: popthv(ngrid,nlayer,L_NSPECTV,8) ! VI optical properties of haze within narrowbands (dtau,tau,k,w,g,drayaer,taugaz,dcont).
	92	REAL,INTENT(OUT) :: poptti(ngrid,nlayer,L_NSPECTI,8) ! IR optical properties of haze and clouds within narrowbands (dtau,tau,k,w,g,drayaer,taugaz,dcont).
	93	REAL,INTENT(OUT) :: popttv(ngrid,nlayer,L_NSPECTV,8) ! VI optical properties of haze and clouds within narrowbands (dtau,tau,k,w,g,drayaer,taugaz,dcont).
[1483]	94
	95
[3083]	96
[253]	97	!-----------------------------------------------------------------------
	98	! Declaration of the variables required by correlated-k subroutines
[1483]	99	! Numbered from top to bottom (unlike in the GCM)
	100	!-----------------------------------------------------------------------
[253]	101
	102	REAL*8 tmid(L_LEVELS),pmid(L_LEVELS)
	103	REAL*8 tlevrad(L_LEVELS),plevrad(L_LEVELS)
	104
[1483]	105	! Optical values for the optci/cv subroutines
[253]	106	REAL*8 stel(L_NSPECTV),stel_fract(L_NSPECTV)
[3090]	107	! IR
[253]	108	REAL*8 dtaui(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
[3090]	109	REAL*8 taucumi(L_LEVELS,L_NSPECTI,L_NGAUSS)
	110	REAL*8 wbari(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
	111	REAL*8 cosbi(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
	112	! VI
[253]	113	REAL*8 dtauv(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
[3090]	114	REAL*8 tauv(L_NLEVRAD,L_NSPECTV,L_NGAUSS)
	115	REAL*8 taucumv(L_LEVELS,L_NSPECTV,L_NGAUSS)
	116	REAL*8 wbarv(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
[253]	117	REAL*8 cosbv(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
[3083]	118
[3090]	119	! Temporary optical values for the optci/cv subroutines (clear column)
	120	REAL*8 dtaui_cc(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
	121	REAL*8 dtauv_cc(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
	122	REAL*8 tauv_cc(L_NLEVRAD,L_NSPECTV,L_NGAUSS)
	123	REAL*8 taucumi_cc(L_LEVELS,L_NSPECTI,L_NGAUSS)
	124	REAL*8 taucumv_cc(L_LEVELS,L_NSPECTV,L_NGAUSS)
[3083]	125	REAL*8 wbari_cc(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
	126	REAL*8 wbarv_cc(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
[3090]	127	REAL*8 cosbi_cc(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
	128	REAL*8 cosbv_cc(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
	129	! Temporary optical values for the optci/cv subroutines (dark column)
	130	REAL*8 dtaui_dc(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
	131	REAL*8 dtauv_dc(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
[3083]	132	REAL*8 tauv_dc(L_NLEVRAD,L_NSPECTV,L_NGAUSS)
[3090]	133	REAL*8 taucumi_dc(L_LEVELS,L_NSPECTI,L_NGAUSS)
[3083]	134	REAL*8 taucumv_dc(L_LEVELS,L_NSPECTV,L_NGAUSS)
[3090]	135	REAL*8 wbari_dc(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
	136	REAL*8 wbarv_dc(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
	137	REAL*8 cosbi_dc(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
	138	REAL*8 cosbv_dc(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
	139
	140	! Optical diagnostics
[3318]	141	! Haze
	142	REAL*8 diag_opthi(L_LEVELS,L_NSPECTI,6)
	143	REAL*8 diag_opthv(L_LEVELS,L_NSPECTV,6)
	144	! Clouds
	145	REAL*8 diag_optti(L_LEVELS,L_NSPECTI,6)
	146	REAL*8 diag_opttv(L_LEVELS,L_NSPECTV,6)
[3090]	147	! Temporary optical diagnostics (clear column)
[3318]	148	REAL*8 diag_optti_cc(L_LEVELS,L_NSPECTI,6)
	149	REAL*8 diag_opttv_cc(L_LEVELS,L_NSPECTV,6)
[3090]	150	! Temporary optical diagnostics (dark column)
[3318]	151	REAL*8 diag_optti_dc(L_LEVELS,L_NSPECTI,6)
	152	REAL*8 diag_opttv_dc(L_LEVELS,L_NSPECTV,6)
[3083]	153
[3090]	154
[961]	155	REAL*8 nfluxtopv,nfluxtopi,nfluxtop,fluxtopvdn
[1483]	156	REAL*8 nfluxoutv_nu(L_NSPECTV) ! Outgoing band-resolved VI flux at TOA (W/m2).
	157	REAL*8 nfluxtopi_nu(L_NSPECTI) ! Net band-resolved IR flux at TOA (W/m2).
	158	REAL*8 fluxupi_nu(L_NLAYRAD,L_NSPECTI) ! For 1D diagnostic.
[253]	159	REAL*8 fmneti(L_NLAYRAD),fmnetv(L_NLAYRAD)
	160	REAL*8 fluxupv(L_NLAYRAD),fluxupi(L_NLAYRAD)
	161	REAL*8 fluxdnv(L_NLAYRAD),fluxdni(L_NLAYRAD)
[1482]	162	REAL*8 albi,acosz
[1483]	163	REAL*8 albv(L_NSPECTV) ! Spectral Visible Albedo.
[253]	164
[3083]	165	INTEGER ig,l,k,nw,iq,ip,ilay,it,lev2lay,cdcolumn,ng
[2050]	166
	167	LOGICAL found
[253]	168
	169	real*8 taugsurf(L_NSPECTV,L_NGAUSS-1)
	170	real*8 taugsurfi(L_NSPECTI,L_NGAUSS-1)
	171
[2366]	172	! Miscellaneous :
	173	character(len=100) :: message
	174	character(len=10),parameter :: subname="callcorrk"
	175
[253]	176	logical OLRz
	177	real*8 NFLUXGNDV_nu(L_NSPECTV)
[1648]	178
[1483]	179	! Included by MT for albedo calculations.
[1482]	180	REAL*8 albedo_temp(L_NSPECTV) ! For equivalent albedo calculation.
[1526]	181	REAL*8 surface_stellar_flux ! Stellar flux reaching the surface. Useful for equivalent albedo calculation.
[2046]	182
	183	! For variable haze
	184	REAL*8 seashazefact(L_LEVELS)
[305]	185
[2050]	186	! For muphys optics
	187	REAL*8 pqmo(ngrid,nlayer,nmicro) ! Tracers for microphysics optics (X/m2).
	188	REAL*8 i2e(ngrid,nlayer) ! int 2 ext factor ( X.kg-1 -> X.m-2 for optics )
	189
	190	! For corr-k recombining
	191	REAL*8 pqr(ngrid,L_PINT,L_REFVAR) ! Tracers for corr-k recombining (mol/mol).
	192	REAL*8 fact, tmin, tmax
	193
	194	LOGICAL usept(L_PINT,L_NTREF) ! mask if pfref grid point will be used
	195	INTEGER inflay(L_PINT) ! nearest inferior GCM layer for pfgasref grid points
	196
[253]	197	!=======================================================================
[1822]	198	! I. Initialization on every call
[1483]	199	!=======================================================================
	200
[1529]	201
[1483]	202	! How much light do we get ?
[253]	203	do nw=1,L_NSPECTV
	204	stel(nw)=stellarf(nw)/(dist_star**2)
	205	end do
	206
[2050]	207	! Convert (microphysical) tracers for optics: X.kg-1 --> X.m-2
	208	! NOTE: it should be moved somewhere else: calmufi performs the same kind of
	209	! computations... waste of time...
	210	i2e(:,1:nlayer) = ( pplev(:,1:nlayer)-pplev(:,2:nlayer+1) ) / gzlat(:,1:nlayer)
	211	pqmo(:,:,:) = 0.0
	212	DO iq=1,nmicro
	213	pqmo(:,:,iq) = pq(:,:,iq)*i2e(:,:)
	214	ENDDO
[253]	215
[2050]	216	! Default value for fixed species for whom vmr has been taken
	217	! into account while computing high-resolution spectra
	218	if (corrk_recombin) pqr(:,:,:) = 1.0
	219
[1483]	220	!-----------------------------------------------------------------------
	221	do ig=1,ngrid ! Starting Big Loop over every GCM column
[253]	222	!-----------------------------------------------------------------------
[2050]	223
	224	! Recombine reference corr-k if needed
	225	if (corrk_recombin) then
	226
	227	! NB : To have decent CPU time recombining is not done on all gridpoints and wavelenghts but we
	228	! calculate a gasi/v_recomb variable on the reference corrk-k T,P grid (only for T,P values
	229	! who match the atmospheric conditions ) which is then processed as a standard pre-mix in
	230	! optci/v routines, but updated every time tracers on the ref P grid have varied > 1%.
[253]	231
[2050]	232	! Extract tracers for variable radiative species
	233	! Also find the nearest GCM layer under each ref pressure
	234	do ip=1,L_PINT
	235
	236	ilay=0
	237	found = .false.
	238	do l=1,nlayer
	239	if ( pplay(ig,l) .gt. 10.0**(pfgasref(ip)+2.0) ) then ! pfgasref=log(p[mbar])
	240	found=.true.
	241	ilay=l
	242	endif
	243	enddo
[1483]	244
[2050]	245	if (.not. found ) then ! set to min
	246	do iq=1,L_REFVAR
	247	if ( radvar_mask(iq) ) then
	248	pqr(ig,ip,iq) = pq(ig,1,radvar_indx(iq)) / rat_mmol(radvar_indx(iq)-nmicro) ! mol/mol
	249	endif
	250	enddo
	251	else
	252	if (ilay==nlayer) then ! set to max
	253	do iq=1,L_REFVAR
	254	if ( radvar_mask(iq) ) then
	255	pqr(ig,ip,iq) = pq(ig,nlayer,radvar_indx(iq)) / rat_mmol(radvar_indx(iq)-nmicro) ! mol/mol
	256	endif
	257	enddo
	258	else ! standard
[2379]	259	fact = ( 10.0**(pfgasref(ip)+2.0) - pplay(ig,ilay+1) ) / ( pplay(ig,ilay) - pplay(ig,ilay+1) ) ! pfgasref=log(p[mbar])
[2050]	260	do iq=1,L_REFVAR
	261	if ( radvar_mask(iq) ) then
	262	pqr(ig,ip,iq) = pq(ig,ilay,radvar_indx(iq))*fact pq(ig,ilay+1,radvar_indx(iq))**(1.0-fact)
	263	pqr(ig,ip,iq) = pqr(ig,ip,iq) / rat_mmol(radvar_indx(iq)-nmicro) ! mol/mol
	264	endif
	265	enddo
	266	endif ! if ilay==nlayer
	267	endif ! if not found
	268
	269	inflay(ip) = ilay
	270
	271	enddo ! ip=1,L_PINT
	272
	273	! NB : The following usept is a trick to call recombine only for the reference T-P
	274	! grid points that are useful given the temperature range at this altitude
	275	! It saves a looot of time - JVO 18
	276	usept(:,:) = .true.
	277	do ip=1,L_PINT-1
	278	if ( inflay(ip+1)==nlayer ) then
	279	usept(ip,:) = .false.
	280	endif
	281	if ( inflay(ip) == 0 ) then
	282	usept(ip+1:,:) = .false.
	283	endif
	284	if ( usept(ip,1) ) then ! if not all false
[2408]	285	tmin = minval(pt(ig,max(inflay(ip+1)+1,1):min(inflay(ip)+1,nlayer)))
	286	tmax = maxval(pt(ig,max(inflay(ip+1)+1,1):min(inflay(ip)+1,nlayer)))
[2050]	287	do it=1,L_NTREF-1
	288	if ( tgasref(it+1) .lt. tmin ) then
	289	usept(ip,it) = .false.
	290	endif
	291	enddo
	292	do it=2,L_NTREF
	293	if ( tgasref(it-1) .gt. tmax ) then
	294	usept(ip,it) = .false.
	295	endif
	296	enddo
	297	! in case of out-of-bounds
	298	if ( tgasref(1) .lt. tmin ) usept(ip,1) = .true.
	299	if ( tgasref(L_NTREF) .gt. tmax ) usept(ip,L_NTREF) = .true.
	300	endif
	301	enddo ! ip=1,L_PINT-1
	302	! deal with last bound
	303	if ( inflay(L_PINT-1).ne.0 ) usept(L_PINT,:) = usept(L_PINT-1,:)
	304
	305
	306	do ip=1,L_PINT
	307
	308	! Recombine k at (useful only!) reference T-P values if tracers or T have enough varied
	309	do it=1,L_NTREF
	310
	311	if ( usept(ip,it) .eqv. .false. ) cycle
	312
	313	do l=1,L_REFVAR
	314	if ( abs( (pqr(ig,ip,l) - pqrold(ip,l)) / max(1.0e-30,pqrold(ip,l))) .GT. 0.01 & ! +- 1%
	315	.or. ( useptold(ip,it) .eqv. .false. ) ) then ! in case T change but not the tracers
	316	call recombin_corrk( pqr(ig,ip,:),ip,it )
	317	exit ! one is enough to trigger the update
	318	endif
	319	enddo
	320
	321	enddo
	322
	323	enddo ! ip=1,L_PINT
	324
	325	useptold(:,:)=usept(:,:)
	326
	327	endif ! if corrk_recombin
	328
[253]	329	!=======================================================================
[1483]	330	! II. Transformation of the GCM variables
	331	!=======================================================================
[253]	332
[1483]	333
	334	! Albedo and Emissivity.
	335	albi=1-emis(ig) ! Long Wave.
	336	DO nw=1,L_NSPECTV ! Short Wave loop.
[1482]	337	albv(nw)=albedo(ig,nw)
[1529]	338	ENDDO
[253]	339
[1483]	340	if ((ngrid.eq.1).and.(global1d)) then ! Fixed zenith angle 'szangle' in 1D simulations w/ globally-averaged sunlight.
[253]	341	acosz = cos(pi*szangle/180.0)
	342	print*,'acosz=',acosz,', szangle=',szangle
	343	else
[1483]	344	acosz=mu0(ig) ! Cosine of sun incident angle : 3D simulations or local 1D simulations using latitude.
[253]	345	endif
	346
	347	!-----------------------------------------------------------------------
	348	! Pressure and temperature
[1483]	349	!-----------------------------------------------------------------------
[253]	350
	351	DO l=1,nlayer
	352	plevrad(2*l) = pplay(ig,nlayer+1-l)/scalep
	353	plevrad(2*l+1) = pplev(ig,nlayer+1-l)/scalep
	354	tlevrad(2*l) = pt(ig,nlayer+1-l)
	355	tlevrad(2*l+1) = (pt(ig,nlayer+1-l)+pt(ig,max(nlayer-l,1)))/2
	356	END DO
	357
[600]	358	plevrad(1) = 0.
[1483]	359	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]	360
	361	tlevrad(1) = tlevrad(2)
[1308]	362	tlevrad(2*nlayer+1)=tsurf(ig)
[253]	363
[1423]	364	pmid(1) = max(pgasmin,0.0001*plevrad(3))
	365	pmid(2) = pmid(1)
	366
[253]	367	tmid(1) = tlevrad(2)
[1423]	368	tmid(2) = tmid(1)
	369
	370	DO l=1,L_NLAYRAD-1
	371	tmid(2l+1) = tlevrad(2l+1)
	372	tmid(2l+2) = tlevrad(2l+1)
	373	pmid(2l+1) = plevrad(2l+1)
	374	pmid(2l+2) = plevrad(2l+1)
[253]	375	END DO
[1423]	376	pmid(L_LEVELS) = plevrad(L_LEVELS)
	377	tmid(L_LEVELS) = tlevrad(L_LEVELS)
[253]	378
[1423]	379	!!Alternative interpolation:
	380	! pmid(3) = pmid(1)
	381	! pmid(4) = pmid(1)
	382	! tmid(3) = tmid(1)
	383	! tmid(4) = tmid(1)
	384	! DO l=2,L_NLAYRAD-1
	385	! tmid(2l+1) = tlevrad(2l)
	386	! tmid(2l+2) = tlevrad(2l)
	387	! pmid(2l+1) = plevrad(2l)
	388	! pmid(2l+2) = plevrad(2l)
	389	! END DO
	390	! pmid(L_LEVELS) = plevrad(L_LEVELS-1)
	391	! tmid(L_LEVELS) = tlevrad(L_LEVELS-1)
	392
[1483]	393	! Test for out-of-bounds pressure.
[253]	394	if(plevrad(3).lt.pgasmin)then
	395	print*,'Minimum pressure is outside the radiative'
	396	print*,'transfer kmatrix bounds, exiting.'
	397	call abort
	398	elseif(plevrad(L_LEVELS).gt.pgasmax)then
	399	print*,'Maximum pressure is outside the radiative'
	400	print*,'transfer kmatrix bounds, exiting.'
	401	call abort
	402	endif
	403
[1483]	404	! Test for out-of-bounds temperature.
[253]	405	do k=1,L_LEVELS
	406	if(tlevrad(k).lt.tgasmin)then
	407	print*,'Minimum temperature is outside the radiative'
[1145]	408	print*,'transfer kmatrix bounds'
[858]	409	print*,"k=",k," tlevrad(k)=",tlevrad(k)
	410	print*,"tgasmin=",tgasmin
[1145]	411	if (strictboundcorrk) then
	412	call abort
	413	else
	414	print,'**********************************************'
	415	print*,'we allow model to continue with tlevrad=tgasmin'
	416	print*,' ... we assume we know what you are doing ... '
	417	print*,' ... but do not let this happen too often ... '
	418	print,'**********************************************'
[1147]	419	!tlevrad(k)=tgasmin
[1145]	420	endif
[253]	421	elseif(tlevrad(k).gt.tgasmax)then
[1648]	422	! print*,'Maximum temperature is outside the radiative'
	423	! print*,'transfer kmatrix bounds, exiting.'
	424	! print*,"k=",k," tlevrad(k)=",tlevrad(k)
	425	! print*,"tgasmax=",tgasmax
[1145]	426	if (strictboundcorrk) then
	427	call abort
	428	else
[1648]	429	! print,'**********************************************'
	430	! print*,'we allow model to continue with tlevrad=tgasmax'
	431	! print*,' ... we assume we know what you are doing ... '
	432	! print*,' ... but do not let this happen too often ... '
	433	! print,'**********************************************'
[1147]	434	!tlevrad(k)=tgasmax
[1145]	435	endif
[253]	436	endif
[2366]	437
	438	if (tlevrad(k).lt.tplanckmin) then
	439	print*,'Minimum temperature is outside the boundaries for'
	440	print*,'Planck function integration set in callphys.def, aborting.'
	441	print*,"k=",k," tlevrad(k)=",tlevrad(k)
	442	print*,"tplanckmin=",tplanckmin
	443	message="Minimum temperature outside Planck function bounds - Change tplanckmin in callphys.def"
	444	call abort_physic(subname,message,1)
	445	else if (tlevrad(k).gt.tplanckmax) then
	446	print*,'Maximum temperature is outside the boundaries for'
	447	print*,'Planck function integration set in callphys.def, aborting.'
	448	print*,"k=",k," tlevrad(k)=",tlevrad(k)
	449	print*,"tplanckmax=",tplanckmax
	450	message="Maximum temperature outside Planck function bounds - Change tplanckmax in callphys.def"
	451	call abort_physic(subname,message,1)
	452	endif
	453
[253]	454	enddo
[2366]	455
[1016]	456	do k=1,L_NLAYRAD+1
	457	if(tmid(k).lt.tgasmin)then
	458	print*,'Minimum temperature is outside the radiative'
	459	print*,'transfer kmatrix bounds, exiting.'
[1145]	460	print*,"k=",k," tmid(k)=",tmid(k)
[1016]	461	print*,"tgasmin=",tgasmin
[1145]	462	if (strictboundcorrk) then
	463	call abort
	464	else
	465	print,'**********************************************'
	466	print*,'we allow model to continue with tmid=tgasmin'
	467	print*,' ... we assume we know what you are doing ... '
	468	print*,' ... but do not let this happen too often ... '
	469	print,'**********************************************'
	470	tmid(k)=tgasmin
	471	endif
[1016]	472	elseif(tmid(k).gt.tgasmax)then
[1648]	473	! print*,'Maximum temperature is outside the radiative'
	474	! print*,'transfer kmatrix bounds, exiting.'
	475	! print*,"k=",k," tmid(k)=",tmid(k)
	476	! print*,"tgasmax=",tgasmax
[1145]	477	if (strictboundcorrk) then
	478	call abort
	479	else
[1648]	480	! print,'**********************************************'
	481	! print*,'we allow model to continue with tmid=tgasmin'
	482	! print*,' ... we assume we know what you are doing ... '
	483	! print*,' ... but do not let this happen too often ... '
	484	! print,'**********************************************'
[1145]	485	tmid(k)=tgasmax
	486	endif
[1016]	487	endif
	488	enddo
[253]	489
	490	!=======================================================================
[1483]	491	! III. Calling the main radiative transfer subroutines
	492	!=======================================================================
[253]	493
[1947]	494	Cmk(:) = 0.01 * 1.0 / (gzlat(ig,:) * mugaz * 1.672621e-27) ! q_main=1.0 assumed.
	495	gzlat_ig(:) = gzlat(ig,:)
[2046]	496
	497	! Compute the haze seasonal modulation factor
	498	if (seashaze) call season_haze(zday,latitude(ig),plevrad,seashazefact)
[253]	499
	500	!-----------------------------------------------------------------------
[1483]	501	! Short Wave Part
	502	!-----------------------------------------------------------------------
[253]	503
[3083]	504	! Clear column :
	505	cdcolumn = 0
	506	call optcv(pqmo(ig,:,:),nlayer,zzlev(ig,:),plevrad,tmid,pmid, &
	507	dtauv_cc,tauv_cc,taucumv_cc,wbarv_cc,cosbv_cc,tauray,taugsurf,seashazefact,&
[3318]	508	diag_opthv,diag_opttv_cc,cdcolumn)
[3083]	509	! Dark column :
	510	cdcolumn = 1
	511	call optcv(pqmo(ig,:,:),nlayer,zzlev(ig,:),plevrad,tmid,pmid, &
	512	dtauv_dc,tauv_dc,taucumv_dc,wbarv_dc,cosbv_dc,tauray,taugsurf,seashazefact,&
[3318]	513	diag_opthv,diag_opttv_dc,cdcolumn)
[3083]	514
	515	! Mean opacity, ssa and asf :
	516	where ((exp(-dtauv_cc(:,:,:)).ge.1.d-40) .and. (exp(-dtauv_dc(:,:,:)).ge.1.d-40))
	517	dtauv(:,:,:) = -log((1-Fcloudy)exp(-dtauv_cc(:,:,:)) + Fcloudyexp(-dtauv_dc(:,:,:)))
	518	elsewhere
	519	dtauv(:,:,:) = dtauv_dc(:,:,:) ! No need to average...
	520	endwhere
	521	do ng = 1, L_NGAUSS
	522	do nw = 1, L_NSPECTV
	523	taucumv(1,nw,ng) = 0.0d0
	524	do k = 2, L_LEVELS
	525	if ((exp(-taucumv_cc(k,nw,ng)).ge.1.d-40) .and. (exp(-taucumv_dc(k,nw,ng)).ge.1.d-40)) then
	526	taucumv(k,nw,ng) = taucumv(k-1,nw,ng) - log((1-Fcloudy)exp(-taucumv_cc(k,nw,ng)) + Fcloudyexp(-taucumv_dc(k,nw,ng)))
	527	else
	528	taucumv(k,nw,ng) = taucumv(k-1,nw,ng) + taucumv_dc(k,nw,ng) ! No need to average...
	529	end if
	530	end do
	531	do l = 1, L_NLAYRAD
	532	tauv(l,nw,ng) = taucumv(2*l,nw,ng)
	533	end do
	534	tauv(l,nw,ng) = taucumv(2*L_NLAYRAD+1,nw,ng)
	535	end do
	536	end do
[3318]	537	wbarv = (1-Fcloudy) * wbarv_cc + Fcloudy * wbarv_dc
	538	cosbv = (1-Fcloudy) * cosbv_cc + Fcloudy * cosbv_dc
[3083]	539
	540	! Diagnostics for clouds :
[3090]	541	if (callclouds) then
[3318]	542	where (diag_opttv_cc(:,:,1) .lt. 1.d-30)
	543	diag_opttv_cc(:,:,1) = 1.d-30
[3090]	544	endwhere
[3318]	545	where (diag_opttv_dc(:,:,1) .lt. 1.d-30)
	546	diag_opttv_dc(:,:,1) = 1.d-30
[3090]	547	endwhere
[3318]	548	diag_opttv(:,:,1) = -log( (1-Fcloudy)exp(-diag_opttv_cc(:,:,1)) + Fcloudyexp(-diag_opttv_dc(:,:,1)) )
	549	diag_opttv(:,:,2) = (1-Fcloudy) * diag_opttv_cc(:,:,2) + Fcloudy * diag_opttv_dc(:,:,2)
	550	diag_opttv(:,:,3) = (1-Fcloudy) * diag_opttv_cc(:,:,3) + Fcloudy * diag_opttv_dc(:,:,3)
	551	diag_opttv(:,:,4) = -log( (1-Fcloudy)exp(-diag_opttv_cc(:,:,4)) + Fcloudyexp(-diag_opttv_dc(:,:,4)) )
	552	diag_opttv(:,:,5) = -log( (1-Fcloudy)exp(-diag_opttv_cc(:,:,5)) + Fcloudyexp(-diag_opttv_dc(:,:,5)) )
	553	diag_opttv(:,:,6) = -log( (1-Fcloudy)exp(-diag_opttv_cc(:,:,6)) + Fcloudyexp(-diag_opttv_dc(:,:,6)) )
[3090]	554	endif
[3083]	555
[1483]	556	if(fract(ig) .ge. 1.0e-4) then ! Only during daylight.
[787]	557	if((ngrid.eq.1).and.(global1d))then
[253]	558	do nw=1,L_NSPECTV
[1483]	559	stel_fract(nw)= stel(nw)* 0.25 / acosz ! globally averaged = divide by 4, and we correct for solar zenith angle
[253]	560	end do
	561	else
	562	do nw=1,L_NSPECTV
[1161]	563	stel_fract(nw)= stel(nw) * fract(ig)
[253]	564	end do
[1483]	565	endif
[253]	566	call sfluxv(dtauv,tauv,taucumv,albv,dwnv,wbarv,cosbv, &
[1781]	567	acosz,stel_fract, &
	568	nfluxtopv,fluxtopvdn,nfluxoutv_nu,nfluxgndv_nu, &
[253]	569	fmnetv,fluxupv,fluxdnv,fzerov,taugsurf)
	570
[1483]	571	else ! During the night, fluxes = 0.
[3083]	572	nfluxtopv = 0.0d0
	573	fluxtopvdn = 0.0d0
	574	nfluxoutv_nu(:) = 0.0d0
	575	nfluxgndv_nu(:) = 0.0d0
[253]	576	do l=1,L_NLAYRAD
[962]	577	fmnetv(l)=0.0d0
	578	fluxupv(l)=0.0d0
	579	fluxdnv(l)=0.0d0
[253]	580	end do
	581	end if
	582
[1482]	583
[1526]	584	! Equivalent Albedo Calculation (for OUTPUT). MT2015
	585	if(fract(ig) .ge. 1.0e-4) then ! equivalent albedo makes sense only during daylight.
	586	surface_stellar_flux=sum(nfluxgndv_nu(1:L_NSPECTV))
	587	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]	588	DO nw=1,L_NSPECTV
	589	albedo_temp(nw)=albedo(ig,nw)*nfluxgndv_nu(nw)
[1526]	590	ENDDO
[1529]	591	albedo_temp(1:L_NSPECTV)=albedo_temp(1:L_NSPECTV)/surface_stellar_flux
[1526]	592	albedo_equivalent(ig)=sum(albedo_temp(1:L_NSPECTV))
	593	else
	594	albedo_equivalent(ig)=0.0 ! Spectrally Integrated Albedo not defined for non-irradiated grid points. So we arbitrary set the equivalent albedo to 0.
	595	endif
[1529]	596	else
	597	albedo_equivalent(ig)=0.0 ! Spectrally Integrated Albedo not defined for non-irradiated grid points. So we arbitrary set the equivalent albedo to 0.
	598	endif
[1482]	599
	600
[253]	601	!-----------------------------------------------------------------------
[1483]	602	! Long Wave Part
	603	!-----------------------------------------------------------------------
[253]	604
[538]	605
[3083]	606	! Clear column :
	607	cdcolumn = 0
	608	call optci(pqmo(ig,:,:),nlayer,zzlev(ig,:),plevrad,tlevrad,tmid,pmid,&
	609	dtaui_cc,taucumi_cc,cosbi_cc,wbari_cc,taugsurfi,seashazefact, &
[3318]	610	diag_opthi,diag_optti_cc,cdcolumn)
[3083]	611	! Dark column :
	612	cdcolumn = 1
	613	call optci(pqmo(ig,:,:),nlayer,zzlev(ig,:),plevrad,tlevrad,tmid,pmid,&
	614	dtaui_dc,taucumi_dc,cosbi_dc,wbari_dc,taugsurfi,seashazefact, &
[3318]	615	diag_opthi,diag_optti_dc,cdcolumn)
[3083]	616
	617	! Mean opacity, ssa and asf :
	618	where ((exp(-dtaui_cc(:,:,:)).ge.1.d-40) .and. (exp(-dtaui_dc(:,:,:)).ge.1.d-40))
	619	dtaui(:,:,:) = -log((1-Fcloudy)exp(-dtaui_cc(:,:,:)) + Fcloudyexp(-dtaui_dc(:,:,:)))
	620	elsewhere
	621	dtaui(:,:,:) = dtaui_dc(:,:,:) ! No need to average...
	622	endwhere
	623	do ng = 1, L_NGAUSS
	624	do nw = 1, L_NSPECTI
	625	taucumi(1,nw,ng) = 0.0d0
	626	do k = 2, L_LEVELS
	627	if ((exp(-taucumi_cc(k,nw,ng)).ge.1.d-40) .and. (exp(-taucumi_dc(k,nw,ng)).ge.1.d-40)) then
	628	taucumi(k,nw,ng) = taucumi(k-1,nw,ng) - log((1-Fcloudy)exp(-taucumi_cc(k,nw,ng)) + Fcloudyexp(-taucumi_dc(k,nw,ng)))
	629	else
	630	taucumi(k,nw,ng) = taucumi(k-1,nw,ng) + taucumi_dc(k,nw,ng) ! No need to average...
	631	end if
	632	end do
	633	end do
	634	end do
[3318]	635	wbari = (1-Fcloudy) * wbari_cc + Fcloudy * wbari_dc
	636	cosbi = (1-Fcloudy) * cosbi_cc + Fcloudy * cosbi_dc
[3083]	637
	638	! Diagnostics for clouds :
[3090]	639	if (callclouds) then
[3318]	640	where (diag_optti_cc(:,:,1) .lt. 1.d-30)
	641	diag_optti_cc(:,:,1) = 1.d-30
[3090]	642	endwhere
[3318]	643	where (diag_optti_dc(:,:,1) .lt. 1.d-30)
	644	diag_optti_dc(:,:,1) = 1.d-30
[3090]	645	endwhere
[3318]	646	diag_optti(:,:,1) = -log( (1-Fcloudy)exp(-diag_optti_cc(:,:,1)) + Fcloudyexp(-diag_optti_dc(:,:,1)) )
	647	diag_optti(:,:,2) = (1-Fcloudy) * diag_optti_cc(:,:,2) + Fcloudy * diag_optti_dc(:,:,2)
	648	diag_optti(:,:,3) = (1-Fcloudy) * diag_optti_cc(:,:,3) + Fcloudy * diag_optti_dc(:,:,3)
	649	diag_optti(:,:,4) = -log( (1-Fcloudy)exp(-diag_optti_cc(:,:,4)) + Fcloudyexp(-diag_optti_dc(:,:,4)) )
	650	diag_optti(:,:,5) = -log( (1-Fcloudy)exp(-diag_optti_cc(:,:,5)) + Fcloudyexp(-diag_optti_dc(:,:,5)) )
	651	diag_optti(:,:,6) = -log( (1-Fcloudy)exp(-diag_optti_cc(:,:,6)) + Fcloudyexp(-diag_optti_dc(:,:,6)) )
[3090]	652	endif
[3083]	653
[253]	654	call sfluxi(plevrad,tlevrad,dtaui,taucumi,ubari,albi, &
[1781]	655	wnoi,dwni,cosbi,wbari,nfluxtopi,nfluxtopi_nu, &
[253]	656	fmneti,fluxupi,fluxdni,fluxupi_nu,fzeroi,taugsurfi)
	657
	658	!-----------------------------------------------------------------------
	659	! Transformation of the correlated-k code outputs
	660	! (into dtlw, dtsw, fluxsurf_lw, fluxsurf_sw, fluxtop_lw, fluxtop_sw)
	661
	662	! Flux incident at the top of the atmosphere
[961]	663	fluxtop_dn(ig)=fluxtopvdn
[253]	664
	665	fluxtop_lw(ig) = real(nfluxtopi)
	666	fluxabs_sw(ig) = real(-nfluxtopv)
	667	fluxsurf_lw(ig) = real(fluxdni(L_NLAYRAD))
	668	fluxsurf_sw(ig) = real(fluxdnv(L_NLAYRAD))
[1482]	669
	670	! Flux absorbed by the surface. By MT2015.
	671	fluxsurfabs_sw(ig) = fluxsurf_sw(ig)*(1.-albedo_equivalent(ig))
[253]	672
	673	if(fluxtop_dn(ig).lt.0.0)then
	674	print*,'Achtung! fluxtop_dn has lost the plot!'
	675	print*,'fluxtop_dn=',fluxtop_dn(ig)
	676	print*,'acosz=',acosz
	677	print*,'temp= ',pt(ig,:)
	678	print*,'pplay= ',pplay(ig,:)
	679	call abort
	680	endif
	681
	682	! Spectral output, for exoplanet observational comparison
	683	if(specOLR)then
	684	do nw=1,L_NSPECTI
[526]	685	OLR_nu(ig,nw)=nfluxtopi_nu(nw)/DWNI(nw) !JL Normalize to the bandwidth
[253]	686	end do
	687	do nw=1,L_NSPECTV
[366]	688	!GSR_nu(ig,nw)=nfluxgndv_nu(nw)
[526]	689	OSR_nu(ig,nw)=nfluxoutv_nu(nw)/DWNV(nw) !JL Normalize to the bandwidth
[253]	690	end do
	691	endif
	692
	693	! Finally, the heating rates
	694
[586]	695	DO l=2,L_NLAYRAD
	696	dtsw(ig,L_NLAYRAD+1-l)=(fmnetv(l)-fmnetv(l-1)) &
[1947]	697	gzlat(ig,L_NLAYRAD+1-l)/(cppscalep(plevrad(2l+1)-plevrad(2*l-1)))
[586]	698	dtlw(ig,L_NLAYRAD+1-l)=(fmneti(l)-fmneti(l-1)) &
[1947]	699	gzlat(ig,L_NLAYRAD+1-l)/(cppscalep(plevrad(2l+1)-plevrad(2*l-1)))
[586]	700	END DO
[253]	701
	702	! These are values at top of atmosphere
[586]	703	dtsw(ig,L_NLAYRAD)=(fmnetv(1)-nfluxtopv) &
[1947]	704	gzlat(ig,L_NLAYRAD)/(cppscalep*(plevrad(3)-plevrad(1)))
[586]	705	dtlw(ig,L_NLAYRAD)=(fmneti(1)-nfluxtopi) &
[1947]	706	gzlat(ig,L_NLAYRAD)/(cppscalep*(plevrad(3)-plevrad(1)))
[253]	707
	708
[2133]	709	! Optical thickness diagnostics (added by JVO)
	710	if (diagdtau) then
[3083]	711	do l=1,L_NLAYRAD
	712	do nw=1,L_NSPECTV
	713	int_dtauv(ig,l,nw) = 0.0d0
	714	DO k=1,L_NGAUSS
	715	! Output exp(-tau) because gweight ponderates exp and not tau itself
	716	int_dtauv(ig,l,nw)= int_dtauv(ig,l,nw) + exp(-dtauv(l,nw,k))*gweight(k)
	717	ENDDO
	718	enddo
	719	do nw=1,L_NSPECTI
	720	int_dtaui(ig,l,nw) = 0.0d0
	721	DO k=1,L_NGAUSS
	722	! Output exp(-tau) because gweight ponderates exp and not tau itself
	723	int_dtaui(ig,l,nw)= int_dtaui(ig,l,nw) + exp(-dtaui(l,nw,k))*gweight(k)
	724	ENDDO
	725	enddo
	726	enddo
	727	endif
	728
	729
[3090]	730	! Diagnostics : optical properties of haze and clouds (dtau, tau, k, w, g) :
[3083]	731	do l = 2, L_LEVELS, 2
	732	lev2lay = L_NLAYRAD+1 - l/2
	733	! Visible :
	734	do nw = 1, L_NSPECTV
	735	popthv(ig,lev2lay,nw,:) = 0.0d0
[3318]	736	popttv(ig,lev2lay,nw,:) = 0.0d0
[3083]	737	! Optical thickness (dtau) :
[3318]	738	popthv(ig,lev2lay,nw,1) = (diag_opthv(l,nw,1) + diag_opthv(l+1,nw,1)) / 2.0
[3090]	739	if (callclouds) then
[3318]	740	popttv(ig,lev2lay,nw,1) = (diag_opttv(l,nw,1) + diag_opttv(l+1,nw,1)) / 2.0
[3090]	741	endif
[3083]	742	! Opacity (tau) :
	743	do k = L_NLAYRAD, lev2lay, -1
	744	popthv(ig,lev2lay,nw,2) = popthv(ig,lev2lay,nw,2) + popthv(ig,k,nw,1)
[3090]	745	if (callclouds) then
[3318]	746	popttv(ig,lev2lay,nw,2) = popttv(ig,lev2lay,nw,2) + popttv(ig,k,nw,1)
[3090]	747	endif
[3083]	748	enddo
	749	! Extinction (k) :
	750	popthv(ig,lev2lay,nw,3) = popthv(ig,lev2lay,nw,1) / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay))
[3090]	751	if (callclouds) then
[3318]	752	popttv(ig,lev2lay,nw,3) = popttv(ig,lev2lay,nw,1) / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay))
[3090]	753	endif
[3083]	754	! Simple Scattering Albedo (w) and Asymmetry Parameter (g) :
[3318]	755	popthv(ig,lev2lay,nw,4) = (diag_opthv(l,nw,2) + diag_opthv(l+1,nw,2)) / 2.0
	756	popthv(ig,lev2lay,nw,5) = (diag_opthv(l,nw,3) + diag_opthv(l+1,nw,3)) / 2.0
[3090]	757	if (callclouds) then
[3318]	758	popttv(ig,lev2lay,nw,4) = (diag_opttv(l,nw,2) + diag_opttv(l+1,nw,2)) / 2.0
	759	popttv(ig,lev2lay,nw,5) = (diag_opttv(l,nw,3) + diag_opttv(l+1,nw,3)) / 2.0
[3090]	760	endif
[3318]	761	! DRAYAER, TAUGAS, DCONT :
	762	popthv(ig,lev2lay,nw,6) = (diag_opthv(l,nw,4) + diag_opthv(l+1,nw,4)) / 2.0 / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay))
	763	popthv(ig,lev2lay,nw,7) = (diag_opthv(l,nw,5) + diag_opthv(l+1,nw,5)) / 2.0 / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay))
	764	popthv(ig,lev2lay,nw,8) = (diag_opthv(l,nw,6) + diag_opthv(l+1,nw,6)) / 2.0 / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay))
	765	if (callclouds) then
	766	popttv(ig,lev2lay,nw,6) = (diag_opttv(l,nw,4) + diag_opttv(l+1,nw,4)) / 2.0 / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay))
	767	popttv(ig,lev2lay,nw,7) = (diag_opttv(l,nw,5) + diag_opttv(l+1,nw,5)) / 2.0 / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay))
	768	popttv(ig,lev2lay,nw,8) = (diag_opttv(l,nw,6) + diag_opttv(l+1,nw,6)) / 2.0 / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay))
	769	endif
[3083]	770	enddo
	771	! Infra-Red
	772	do nw=1,L_NSPECTI
	773	popthi(ig,lev2lay,nw,:) = 0.0d0
[3318]	774	poptti(ig,lev2lay,nw,:) = 0.0d0
[3083]	775	! Optical thickness (dtau) :
[3318]	776	popthi(ig,lev2lay,nw,1) = (diag_opthi(l,nw,1) + diag_opthi(l+1,nw,1)) / 2.0
[3090]	777	if (callclouds) then
[3318]	778	poptti(ig,lev2lay,nw,1) = (diag_optti(l,nw,1) + diag_optti(l+1,nw,1)) / 2.0
[3090]	779	endif
[3083]	780	! Opacity (tau) :
	781	do k = L_NLAYRAD, lev2lay, -1
	782	popthi(ig,lev2lay,nw,2) = popthi(ig,lev2lay,nw,2) + popthi(ig,k,nw,1)
[3090]	783	if (callclouds) then
[3318]	784	poptti(ig,lev2lay,nw,2) = poptti(ig,lev2lay,nw,2) + poptti(ig,k,nw,1)
[3090]	785	endif
[3083]	786	enddo
	787	! Extinction (k) :
	788	popthi(ig,lev2lay,nw,3) = popthi(ig,lev2lay,nw,1) / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay))
[3090]	789	if (callclouds) then
[3318]	790	poptti(ig,lev2lay,nw,3) = poptti(ig,lev2lay,nw,1) / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay))
[3090]	791	endif
[3083]	792	! Simple Scattering Albedo (w) and Asymmetry Parameter (g) :
[3318]	793	popthi(ig,lev2lay,nw,4) = (diag_opthi(l,nw,2) + diag_opthi(l+1,nw,2)) / 2.0
	794	popthi(ig,lev2lay,nw,5) = (diag_opthi(l,nw,3) + diag_opthi(l+1,nw,3)) / 2.0
[3090]	795	if (callclouds) then
[3318]	796	poptti(ig,lev2lay,nw,4) = (diag_optti(l,nw,2) + diag_optti(l+1,nw,2)) / 2.0
	797	poptti(ig,lev2lay,nw,5) = (diag_optti(l,nw,3) + diag_optti(l+1,nw,3)) / 2.0
[3090]	798	endif
[3318]	799	! DRAYAER, TAUGAS, DCONT :
	800	popthi(ig,lev2lay,nw,6) = (diag_opthi(l,nw,4) + diag_opthi(l+1,nw,4)) / 2.0 / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay))
	801	popthi(ig,lev2lay,nw,7) = (diag_opthi(l,nw,5) + diag_opthi(l+1,nw,5)) / 2.0 / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay))
	802	popthi(ig,lev2lay,nw,8) = (diag_opthi(l,nw,6) + diag_opthi(l+1,nw,6)) / 2.0 / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay))
	803	if (callclouds) then
	804	poptti(ig,lev2lay,nw,6) = (diag_optti(l,nw,4) + diag_optti(l+1,nw,4)) / 2.0 / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay))
	805	poptti(ig,lev2lay,nw,7) = (diag_optti(l,nw,5) + diag_optti(l+1,nw,5)) / 2.0 / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay))
	806	poptti(ig,lev2lay,nw,8) = (diag_optti(l,nw,6) + diag_optti(l+1,nw,6)) / 2.0 / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay))
	807	endif
[3083]	808	enddo
	809	enddo
[2133]	810
[3083]	811
[1483]	812	!-----------------------------------------------------------------------
	813	end do ! End of big loop over every GCM column.
	814	!-----------------------------------------------------------------------
[253]	815
[1483]	816
[253]	817	!-----------------------------------------------------------------------
	818	! Additional diagnostics
[1483]	819	!-----------------------------------------------------------------------
[253]	820
[1483]	821	! IR spectral output, for exoplanet observational comparison
	822	if(lastcall.and.(ngrid.eq.1))then ! could disable the 1D output, they are in the diagfi and diagspec... JL12
[253]	823
[1483]	824	print*,'Saving scalar quantities in surf_vals.out...'
	825	print*,'psurf = ', pplev(1,1),' Pa'
	826	open(116,file='surf_vals.out')
	827	write(116,*) tsurf(1),pplev(1,1),fluxtop_dn(1), &
	828	real(-nfluxtopv),real(nfluxtopi)
	829	close(116)
[253]	830
[1483]	831	! USEFUL COMMENT - Do Not Remove.
	832	!
[526]	833	! if(specOLR)then
	834	! open(117,file='OLRnu.out')
	835	! do nw=1,L_NSPECTI
	836	! write(117,*) OLR_nu(1,nw)
	837	! enddo
	838	! close(117)
	839	!
	840	! open(127,file='OSRnu.out')
	841	! do nw=1,L_NSPECTV
	842	! write(127,*) OSR_nu(1,nw)
	843	! enddo
	844	! close(127)
	845	! endif
[253]	846
[1483]	847	! OLR vs altitude: do it as a .txt file.
	848	OLRz=.false.
	849	if(OLRz)then
	850	print*,'saving IR vertical flux for OLRz...'
	851	open(118,file='OLRz_plevs.out')
	852	open(119,file='OLRz.out')
	853	do l=1,L_NLAYRAD
	854	write(118,) plevrad(2l)
	855	do nw=1,L_NSPECTI
	856	write(119,*) fluxupi_nu(l,nw)
	857	enddo
	858	enddo
	859	close(118)
	860	close(119)
	861	endif
[253]	862
[305]	863	endif
[253]	864
[1647]	865	if (lastcall) then
[470]	866	IF( ALLOCATED( gasi ) ) DEALLOCATE( gasi )
	867	IF( ALLOCATED( gasv ) ) DEALLOCATE( gasv )
[2050]	868	IF( ALLOCATED( gasi_recomb ) ) DEALLOCATE( gasi_recomb )
	869	IF( ALLOCATED( gasv_recomb ) ) DEALLOCATE( gasv_recomb )
	870	IF( ALLOCATED( pqrold ) ) DEALLOCATE( pqrold )
	871	IF( ALLOCATED( useptold ) ) DEALLOCATE( useptold )
[1315]	872	!$OMP BARRIER
	873	!$OMP MASTER
[470]	874	IF( ALLOCATED( pgasref ) ) DEALLOCATE( pgasref )
	875	IF( ALLOCATED( tgasref ) ) DEALLOCATE( tgasref )
	876	IF( ALLOCATED( pfgasref ) ) DEALLOCATE( pfgasref )
[2026]	877	IF( ALLOCATED( gweight ) ) DEALLOCATE( gweight )
[1315]	878	!$OMP END MASTER
[1529]	879	!$OMP BARRIER
[470]	880	endif
	881
[716]	882
[253]	883	end subroutine callcorrk

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

Download in other formats:

Original Format