Context Navigation

callcorrk.F90 @ 3026

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

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

Download in other formats:

Original Format