Context Navigation

callcorrk.F90 @ 1990

Last change on this file since 1990 was 1988, checked in by jleconte, 7 years ago

28/08/2018 == JL

We now shift the radiative model top from p=0 to the middle of the last physical layer. This is done by changing pmid and plevrad in callcorrk and some corrections need to be done in gfluxv.
This seems to get rid of the aratic temperature behavior in the last two layers of the model (especially on the night side on synchronous planets).
Additional speedup corrections have been made in gfluxi that change nothing to the result.
Finally, if aerosols are present in the last layer we must account for them. Provides better upper boundary condition in the IR. They must however be put to zero in the sw (see optcv and changes in last commit.)
This has been done for water ice in aeropacity, but same correction should probably be done for other aerosol types.

Property svn:executable set to *

File size: 38.8 KB

Rev	Line
[526]	1	subroutine callcorrk(ngrid,nlayer,pq,nq,qsurf, &
[1482]	2	albedo,albedo_equivalent,emis,mu0,pplev,pplay,pt, &
[858]	3	tsurf,fract,dist_star,aerosol,muvar, &
[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, &
[858]	8	tau_col,cloudfrac,totcloudfrac, &
[253]	9	clearsky,firstcall,lastcall)
	10
[1699]	11	use mod_phys_lmdz_para, only : is_master
[253]	12	use radinc_h
	13	use radcommon_h
	14	use watercommon_h
[374]	15	use datafile_mod, only: datadir
[1521]	16	use ioipsl_getin_p_mod, only: getin_p
[471]	17	use gases_h
[1026]	18	use radii_mod, only : su_aer_radii,co2_reffrad,h2o_reffrad,dust_reffrad,h2so4_reffrad,back2lay_reffrad
[1677]	19	use aerosol_mod, only : iaero_co2,iaero_h2o,iaero_dust,iaero_h2so4, iaero_back2lay, iaero_nh3, iaero_aurora
[787]	20	USE tracer_h
[1384]	21	use comcstfi_mod, only: pi, mugaz, cpp
[1709]	22	use callkeys_mod, only: varactive,diurnal,tracer,water,varfixed,satval, &
[1529]	23	kastprof,strictboundcorrk,specOLR,CLFvarying
[787]	24
[253]	25	implicit none
	26
	27	!==================================================================
	28	!
	29	! Purpose
	30	! -------
	31	! Solve the radiative transfer using the correlated-k method for
	32	! the gaseous absorption and the Toon et al. (1989) method for
	33	! scatttering due to aerosols.
	34	!
	35	! Authors
	36	! -------
	37	! Emmanuel 01/2001, Forget 09/2001
	38	! Robin Wordsworth (2009)
	39	!
	40	!==================================================================
	41
	42	!-----------------------------------------------------------------------
	43	! Declaration of the arguments (INPUT - OUTPUT) on the LMD GCM grid
	44	! Layer #1 is the layer near the ground.
[1308]	45	! Layer #nlayer is the layer at the top.
[1483]	46	!-----------------------------------------------------------------------
[253]	47
[1483]	48
	49	! INPUT
	50	INTEGER,INTENT(IN) :: ngrid ! Number of atmospheric columns.
	51	INTEGER,INTENT(IN) :: nlayer ! Number of atmospheric layers.
	52	REAL,INTENT(IN) :: pq(ngrid,nlayer,nq) ! Tracers (kg/kg_of_air).
	53	INTEGER,INTENT(IN) :: nq ! Number of tracers.
	54	REAL,INTENT(IN) :: qsurf(ngrid,nq) ! Tracers on surface (kg.m-2).
	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	REAL,INTENT(IN) :: muvar(ngrid,nlayer+1)
	65	REAL,INTENT(IN) :: cloudfrac(ngrid,nlayer) ! Fraction of clouds (%).
[858]	66	logical,intent(in) :: clearsky
[1483]	67	logical,intent(in) :: firstcall ! Signals first call to physics.
	68	logical,intent(in) :: lastcall ! Signals last call to physics.
	69
	70	! OUTPUT
	71	REAL,INTENT(OUT) :: aerosol(ngrid,nlayer,naerkind) ! Aerosol tau (kg/kg).
	72	REAL,INTENT(OUT) :: dtlw(ngrid,nlayer) ! Heating rate (K/s) due to LW radiation.
	73	REAL,INTENT(OUT) :: dtsw(ngrid,nlayer) ! Heating rate (K/s) due to SW radiation.
	74	REAL,INTENT(OUT) :: fluxsurf_lw(ngrid) ! Incident LW flux to surf (W/m2).
	75	REAL,INTENT(OUT) :: fluxsurf_sw(ngrid) ! Incident SW flux to surf (W/m2)
	76	REAL,INTENT(OUT) :: fluxsurfabs_sw(ngrid) ! Absorbed SW flux by the surface (W/m2). By MT2015.
	77	REAL,INTENT(OUT) :: fluxtop_lw(ngrid) ! Outgoing LW flux to space (W/m2).
	78	REAL,INTENT(OUT) :: fluxabs_sw(ngrid) ! SW flux absorbed by the planet (W/m2).
	79	REAL,INTENT(OUT) :: fluxtop_dn(ngrid) ! Incident top of atmosphere SW flux (W/m2).
	80	REAL,INTENT(OUT) :: OLR_nu(ngrid,L_NSPECTI) ! Outgoing LW radition in each band (Normalized to the band width (W/m2/cm-1).
	81	REAL,INTENT(OUT) :: OSR_nu(ngrid,L_NSPECTV) ! Outgoing SW radition in each band (Normalized to the band width (W/m2/cm-1).
	82	REAL,INTENT(OUT) :: tau_col(ngrid) ! Diagnostic from aeropacity.
	83	REAL,INTENT(OUT) :: albedo_equivalent(ngrid) ! Spectrally Integrated Albedo. For Diagnostic. By MT2015
	84	REAL,INTENT(OUT) :: totcloudfrac(ngrid) ! Column Fraction of clouds (%).
	85
	86
	87
	88
[253]	89
[1483]	90	! Globally varying aerosol optical properties on GCM grid ; not needed everywhere so not in radcommon_h.
[1308]	91	REAL :: QVISsQREF3d(ngrid,nlayer,L_NSPECTV,naerkind)
	92	REAL :: omegaVIS3d(ngrid,nlayer,L_NSPECTV,naerkind)
	93	REAL :: gVIS3d(ngrid,nlayer,L_NSPECTV,naerkind)
	94	REAL :: QIRsQREF3d(ngrid,nlayer,L_NSPECTI,naerkind)
	95	REAL :: omegaIR3d(ngrid,nlayer,L_NSPECTI,naerkind)
	96	REAL :: gIR3d(ngrid,nlayer,L_NSPECTI,naerkind)
[253]	97
[1308]	98	! REAL :: omegaREFvis3d(ngrid,nlayer,naerkind)
	99	! REAL :: omegaREFir3d(ngrid,nlayer,naerkind) ! not sure of the point of these...
[253]	100
[1483]	101	REAL,ALLOCATABLE,SAVE :: reffrad(:,:,:) ! aerosol effective radius (m)
[858]	102	REAL,ALLOCATABLE,SAVE :: nueffrad(:,:,:) ! aerosol effective variance
[1315]	103	!$OMP THREADPRIVATE(reffrad,nueffrad)
[253]	104
	105	!-----------------------------------------------------------------------
	106	! Declaration of the variables required by correlated-k subroutines
[1483]	107	! Numbered from top to bottom (unlike in the GCM)
	108	!-----------------------------------------------------------------------
[253]	109
	110	REAL*8 tmid(L_LEVELS),pmid(L_LEVELS)
	111	REAL*8 tlevrad(L_LEVELS),plevrad(L_LEVELS)
	112
[1483]	113	! Optical values for the optci/cv subroutines
[253]	114	REAL*8 stel(L_NSPECTV),stel_fract(L_NSPECTV)
	115	REAL*8 dtaui(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
	116	REAL*8 dtauv(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
	117	REAL*8 cosbv(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
	118	REAL*8 cosbi(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
	119	REAL*8 wbari(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
	120	REAL*8 wbarv(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
	121	REAL*8 tauv(L_NLEVRAD,L_NSPECTV,L_NGAUSS)
	122	REAL*8 taucumv(L_LEVELS,L_NSPECTV,L_NGAUSS)
	123	REAL*8 taucumi(L_LEVELS,L_NSPECTI,L_NGAUSS)
	124
[1715]	125	REAL*8 tauaero(L_LEVELS,naerkind)
[961]	126	REAL*8 nfluxtopv,nfluxtopi,nfluxtop,fluxtopvdn
[1483]	127	REAL*8 nfluxoutv_nu(L_NSPECTV) ! Outgoing band-resolved VI flux at TOA (W/m2).
	128	REAL*8 nfluxtopi_nu(L_NSPECTI) ! Net band-resolved IR flux at TOA (W/m2).
	129	REAL*8 fluxupi_nu(L_NLAYRAD,L_NSPECTI) ! For 1D diagnostic.
[253]	130	REAL*8 fmneti(L_NLAYRAD),fmnetv(L_NLAYRAD)
	131	REAL*8 fluxupv(L_NLAYRAD),fluxupi(L_NLAYRAD)
	132	REAL*8 fluxdnv(L_NLAYRAD),fluxdni(L_NLAYRAD)
[1482]	133	REAL*8 albi,acosz
[1483]	134	REAL*8 albv(L_NSPECTV) ! Spectral Visible Albedo.
[253]	135
[1483]	136	INTEGER ig,l,k,nw,iaer
[253]	137
[590]	138	real szangle
	139	logical global1d
	140	save szangle,global1d
[1315]	141	!$OMP THREADPRIVATE(szangle,global1d)
[253]	142	real*8 taugsurf(L_NSPECTV,L_NGAUSS-1)
	143	real*8 taugsurfi(L_NSPECTI,L_NGAUSS-1)
[1483]	144	real*8 qvar(L_LEVELS) ! Mixing ratio of variable component (mol/mol).
[253]	145
[1483]	146	! Local aerosol optical properties for each column on RADIATIVE grid.
[1529]	147	real*8,save,allocatable :: QXVAER(:,:,:)
	148	real*8,save,allocatable :: QSVAER(:,:,:)
	149	real*8,save,allocatable :: GVAER(:,:,:)
	150	real*8,save,allocatable :: QXIAER(:,:,:)
	151	real*8,save,allocatable :: QSIAER(:,:,:)
	152	real*8,save,allocatable :: GIAER(:,:,:)
[253]	153
[787]	154	real, dimension(:,:,:), save, allocatable :: QREFvis3d
	155	real, dimension(:,:,:), save, allocatable :: QREFir3d
[1315]	156	!$OMP THREADPRIVATE(QXVAER,QSVAER,GVAER,QXIAER,QSIAER,GIAER,QREFvis3d,QREFir3d)
[787]	157
[253]	158
[1483]	159	! Miscellaneous :
[253]	160	real*8 temp,temp1,temp2,pweight
	161	character(len=10) :: tmp1
	162	character(len=10) :: tmp2
	163
[1483]	164	! For fixed water vapour profiles.
[253]	165	integer i_var
	166	real RH
	167	real*8 pq_temp(nlayer)
[1483]	168	! real(KIND=r8) :: pq_temp(nlayer) ! better F90 way.. DOESNT PORT TO F77!!!
[253]	169	real ptemp, Ttemp, qsat
	170
	171	logical OLRz
	172	real*8 NFLUXGNDV_nu(L_NSPECTV)
	173
[1483]	174	! Included by RW for runaway greenhouse 1D study.
[1308]	175	real vtmp(nlayer)
[305]	176	REAL*8 muvarrad(L_LEVELS)
[1482]	177
[1483]	178	! Included by MT for albedo calculations.
[1482]	179	REAL*8 albedo_temp(L_NSPECTV) ! For equivalent albedo calculation.
[1526]	180	REAL*8 surface_stellar_flux ! Stellar flux reaching the surface. Useful for equivalent albedo calculation.
[305]	181
[1483]	182
[726]	183	!===============================================================
[1483]	184	! I.a Initialization on first call
	185	!===============================================================
[253]	186
[1483]	187
[1529]	188	if(firstcall) then
[253]	189
[1529]	190	! test on allocated necessary because of CLFvarying (two calls to callcorrk in physiq)
[1715]	191	if(.not.allocated(QXVAER)) allocate(QXVAER(L_LEVELS,L_NSPECTV,naerkind))
	192	if(.not.allocated(QSVAER)) allocate(QSVAER(L_LEVELS,L_NSPECTV,naerkind))
	193	if(.not.allocated(GVAER)) allocate(GVAER(L_LEVELS,L_NSPECTV,naerkind))
	194	if(.not.allocated(QXIAER)) allocate(QXIAER(L_LEVELS,L_NSPECTI,naerkind))
	195	if(.not.allocated(QSIAER)) allocate(QSIAER(L_LEVELS,L_NSPECTI,naerkind))
	196	if(.not.allocated(GIAER)) allocate(GIAER(L_LEVELS,L_NSPECTI,naerkind))
[253]	197
[1483]	198	!!! ALLOCATED instances are necessary because of CLFvarying (strategy to call callcorrk twice in physiq...)
[1308]	199	IF(.not.ALLOCATED(QREFvis3d)) ALLOCATE(QREFvis3d(ngrid,nlayer,naerkind))
	200	IF(.not.ALLOCATED(QREFir3d)) ALLOCATE(QREFir3d(ngrid,nlayer,naerkind))
[861]	201	! Effective radius and variance of the aerosols
	202	IF(.not.ALLOCATED(reffrad)) allocate(reffrad(ngrid,nlayer,naerkind))
	203	IF(.not.ALLOCATED(nueffrad)) allocate(nueffrad(ngrid,nlayer,naerkind))
[787]	204
[1829]	205	#ifndef MESOSCALE
[253]	206	call system('rm -f surf_vals_long.out')
[1829]	207	#endif
[253]	208
[728]	209	if(naerkind.gt.4)then
	210	print*,'Code not general enough to deal with naerkind > 4 yet.'
	211	call abort
	212	endif
[1308]	213	call su_aer_radii(ngrid,nlayer,reffrad,nueffrad)
[1483]	214
[1529]	215
[728]	216	!--------------------------------------------------
[1483]	217	! Set up correlated k
	218	!--------------------------------------------------
	219
	220
[374]	221	print*, "callcorrk: Correlated-k data base folder:",trim(datadir)
[1315]	222	call getin_p("corrkdir",corrkdir)
[253]	223	print*, "corrkdir = ",corrkdir
	224	write( tmp1, '(i3)' ) L_NSPECTI
	225	write( tmp2, '(i3)' ) L_NSPECTV
	226	banddir=trim(adjustl(tmp1))//'x'//trim(adjustl(tmp2))
	227	banddir=trim(adjustl(corrkdir))//'/'//trim(adjustl(banddir))
	228
[1483]	229	call setspi ! Basic infrared properties.
	230	call setspv ! Basic visible properties.
	231	call sugas_corrk ! Set up gaseous absorption properties.
	232	call suaer_corrk ! Set up aerosol optical properties.
[1498]	233
[253]	234
	235	if((igcm_h2o_vap.eq.0) .and. varactive)then
	236	print*,'varactive in callcorrk but no h2o_vap tracer.'
	237	stop
	238	endif
	239
[716]	240	OLR_nu(:,:) = 0.
	241	OSR_nu(:,:) = 0.
[538]	242
[787]	243	if (ngrid.eq.1) then
[1483]	244	PRINT*, 'Simulate global averaged conditions ?'
	245	global1d = .false. ! default value
	246	call getin_p("global1d",global1d)
	247	write(,) "global1d = ",global1d
	248
	249	! Test of incompatibility : if global1d is true, there should not be any diurnal cycle.
	250	if (global1d.and.diurnal) then
	251	print*,'if global1d is true, diurnal must be set to false'
	252	stop
	253	endif
[622]	254
[1483]	255	if (global1d) then
	256	PRINT *,'Solar Zenith angle (deg.) ?'
	257	PRINT *,'(assumed for averaged solar flux S/4)'
	258	szangle=60.0 ! default value
	259	call getin_p("szangle",szangle)
	260	write(,) "szangle = ",szangle
	261	endif
[590]	262	endif
	263
[858]	264	end if ! of if (firstcall)
[253]	265
	266	!=======================================================================
[1483]	267	! I.b Initialization on every call
	268	!=======================================================================
	269
[1529]	270	qxvaer(:,:,:)=0.0
	271	qsvaer(:,:,:)=0.0
	272	gvaer(:,:,:) =0.0
	273
	274	qxiaer(:,:,:)=0.0
	275	qsiaer(:,:,:)=0.0
	276	giaer(:,:,:) =0.0
	277
[728]	278	!--------------------------------------------------
	279	! Effective radius and variance of the aerosols
[1483]	280	!--------------------------------------------------
	281
[726]	282	do iaer=1,naerkind
[650]	283
[1483]	284	if ((iaer.eq.iaero_co2).and.tracer.and.(igcm_co2_ice.gt.0)) then ! Treat condensed co2 particles.
[1529]	285	call co2_reffrad(ngrid,nlayer,nq,pq,reffrad(1,1,iaero_co2))
[1699]	286	if (is_master) then
	287	print*,'Max. CO2 ice particle size = ',maxval(reffrad(1:ngrid,1:nlayer,iaer))/1.e-6,' um'
	288	print*,'Min. CO2 ice particle size = ',minval(reffrad(1:ngrid,1:nlayer,iaer))/1.e-6,' um'
	289	end if
	290	end if
[1483]	291
	292	if ((iaer.eq.iaero_h2o).and.water) then ! Treat condensed water particles. To be generalized for other aerosols ...
[1529]	293	call h2o_reffrad(ngrid,nlayer,pq(1,1,igcm_h2o_ice),pt, &
[858]	294	reffrad(1,1,iaero_h2o),nueffrad(1,1,iaero_h2o))
[1699]	295	if (is_master) then
	296	print*,'Max. H2O cloud particle size = ',maxval(reffrad(1:ngrid,1:nlayer,iaer))/1.e-6,' um'
	297	print*,'Min. H2O cloud particle size = ',minval(reffrad(1:ngrid,1:nlayer,iaer))/1.e-6,' um'
	298	end if
[253]	299	endif
[1483]	300
[726]	301	if(iaer.eq.iaero_dust)then
[1529]	302	call dust_reffrad(ngrid,nlayer,reffrad(1,1,iaero_dust))
[1699]	303	if (is_master) then
	304	print*,'Dust particle size = ',reffrad(1,1,iaer)/1.e-6,' um'
	305	end if
[253]	306	endif
[1483]	307
[726]	308	if(iaer.eq.iaero_h2so4)then
[1529]	309	call h2so4_reffrad(ngrid,nlayer,reffrad(1,1,iaero_h2so4))
[1699]	310	if (is_master) then
	311	print*,'H2SO4 particle size =',reffrad(1,1,iaer)/1.e-6,' um'
	312	end if
[253]	313	endif
[1483]	314
[1026]	315	if(iaer.eq.iaero_back2lay)then
[1529]	316	call back2lay_reffrad(ngrid,reffrad(1,1,iaero_back2lay),nlayer,pplev)
[1026]	317	endif
[1677]	318	! if(iaer.eq.iaero_nh3)then
	319	! call nh3_reffrad(ngrid,nlayer,reffrad(1,1,iaero_nh3))
	320	! endif
	321	! if(iaer.eq.iaero_aurora)then
	322	! call aurora_reffrad(ngrid,nlayer,reffrad(1,1,iaero_aurora))
	323	! endif
	324
[1483]	325	end do !iaer=1,naerkind.
[253]	326
[1715]	327
[1483]	328	! How much light do we get ?
[253]	329	do nw=1,L_NSPECTV
	330	stel(nw)=stellarf(nw)/(dist_star**2)
	331	end do
	332
[1483]	333	! Get 3D aerosol optical properties.
[253]	334	call aeroptproperties(ngrid,nlayer,reffrad,nueffrad, &
	335	QVISsQREF3d,omegaVIS3d,gVIS3d, &
	336	QIRsQREF3d,omegaIR3d,gIR3d, &
[1483]	337	QREFvis3d,QREFir3d)
[253]	338
[1483]	339	! Get aerosol optical depths.
[253]	340	call aeropacity(ngrid,nlayer,nq,pplay,pplev,pq,aerosol, &
	341	reffrad,QREFvis3d,QREFir3d, &
[1483]	342	tau_col,cloudfrac,totcloudfrac,clearsky)
[1529]	343
[1483]	344
	345
	346	!-----------------------------------------------------------------------
	347	do ig=1,ngrid ! Starting Big Loop over every GCM column
[253]	348	!-----------------------------------------------------------------------
	349
[1483]	350
[253]	351	!=======================================================================
[1483]	352	! II. Transformation of the GCM variables
	353	!=======================================================================
[253]	354
[1483]	355
[253]	356	!-----------------------------------------------------------------------
[1483]	357	! Aerosol optical properties Qext, Qscat and g.
	358	! The transformation in the vertical is the same as for temperature.
	359	!-----------------------------------------------------------------------
[253]	360
[1483]	361
[253]	362	do iaer=1,naerkind
[1483]	363	! Shortwave.
	364	do nw=1,L_NSPECTV
	365
[1308]	366	do l=1,nlayer
[253]	367
[1308]	368	temp1=QVISsQREF3d(ig,nlayer+1-l,nw,iaer) &
	369	*QREFvis3d(ig,nlayer+1-l,iaer)
[253]	370
[1308]	371	temp2=QVISsQREF3d(ig,max(nlayer-l,1),nw,iaer) &
	372	*QREFvis3d(ig,max(nlayer-l,1),iaer)
[253]	373
	374	qxvaer(2*l,nw,iaer) = temp1
	375	qxvaer(2*l+1,nw,iaer)=(temp1+temp2)/2
	376
[1308]	377	temp1=temp1*omegavis3d(ig,nlayer+1-l,nw,iaer)
	378	temp2=temp2*omegavis3d(ig,max(nlayer-l,1),nw,iaer)
[253]	379
	380	qsvaer(2*l,nw,iaer) = temp1
	381	qsvaer(2*l+1,nw,iaer)=(temp1+temp2)/2
	382
[1308]	383	temp1=gvis3d(ig,nlayer+1-l,nw,iaer)
	384	temp2=gvis3d(ig,max(nlayer-l,1),nw,iaer)
[253]	385
	386	gvaer(2*l,nw,iaer) = temp1
	387	gvaer(2*l+1,nw,iaer)=(temp1+temp2)/2
	388
[1483]	389	end do ! nlayer
[253]	390
	391	qxvaer(1,nw,iaer)=qxvaer(2,nw,iaer)
[1308]	392	qxvaer(2*nlayer+1,nw,iaer)=0.
[253]	393
	394	qsvaer(1,nw,iaer)=qsvaer(2,nw,iaer)
[1308]	395	qsvaer(2*nlayer+1,nw,iaer)=0.
[253]	396
	397	gvaer(1,nw,iaer)=gvaer(2,nw,iaer)
[1308]	398	gvaer(2*nlayer+1,nw,iaer)=0.
[253]	399
[1483]	400	end do ! L_NSPECTV
	401
	402	do nw=1,L_NSPECTI
	403	! Longwave
[1308]	404	do l=1,nlayer
[253]	405
[1308]	406	temp1=QIRsQREF3d(ig,nlayer+1-l,nw,iaer) &
	407	*QREFir3d(ig,nlayer+1-l,iaer)
[253]	408
[1308]	409	temp2=QIRsQREF3d(ig,max(nlayer-l,1),nw,iaer) &
	410	*QREFir3d(ig,max(nlayer-l,1),iaer)
[253]	411
	412	qxiaer(2*l,nw,iaer) = temp1
	413	qxiaer(2*l+1,nw,iaer)=(temp1+temp2)/2
	414
[1308]	415	temp1=temp1*omegair3d(ig,nlayer+1-l,nw,iaer)
	416	temp2=temp2*omegair3d(ig,max(nlayer-l,1),nw,iaer)
[253]	417
	418	qsiaer(2*l,nw,iaer) = temp1
	419	qsiaer(2*l+1,nw,iaer)=(temp1+temp2)/2
	420
[1308]	421	temp1=gir3d(ig,nlayer+1-l,nw,iaer)
	422	temp2=gir3d(ig,max(nlayer-l,1),nw,iaer)
[253]	423
	424	giaer(2*l,nw,iaer) = temp1
	425	giaer(2*l+1,nw,iaer)=(temp1+temp2)/2
	426
[1483]	427	end do ! nlayer
[253]	428
	429	qxiaer(1,nw,iaer)=qxiaer(2,nw,iaer)
[1308]	430	qxiaer(2*nlayer+1,nw,iaer)=0.
[253]	431
	432	qsiaer(1,nw,iaer)=qsiaer(2,nw,iaer)
[1308]	433	qsiaer(2*nlayer+1,nw,iaer)=0.
[253]	434
	435	giaer(1,nw,iaer)=giaer(2,nw,iaer)
[1308]	436	giaer(2*nlayer+1,nw,iaer)=0.
[253]	437
[1483]	438	end do ! L_NSPECTI
	439
	440	end do ! naerkind
[253]	441
[1483]	442	! Test / Correct for freaky s. s. albedo values.
[253]	443	do iaer=1,naerkind
[1715]	444	do k=1,L_LEVELS
[253]	445
	446	do nw=1,L_NSPECTV
	447	if(qsvaer(k,nw,iaer).gt.1.05*qxvaer(k,nw,iaer))then
[726]	448	print*,'Serious problems with qsvaer values'
	449	print*,'in callcorrk'
[253]	450	call abort
	451	endif
	452	if(qsvaer(k,nw,iaer).gt.qxvaer(k,nw,iaer))then
	453	qsvaer(k,nw,iaer)=qxvaer(k,nw,iaer)
	454	endif
	455	end do
	456
	457	do nw=1,L_NSPECTI
	458	if(qsiaer(k,nw,iaer).gt.1.05*qxiaer(k,nw,iaer))then
[726]	459	print*,'Serious problems with qsiaer values'
	460	print*,'in callcorrk'
[253]	461	call abort
	462	endif
	463	if(qsiaer(k,nw,iaer).gt.qxiaer(k,nw,iaer))then
	464	qsiaer(k,nw,iaer)=qxiaer(k,nw,iaer)
	465	endif
	466	end do
	467
[1483]	468	end do ! L_LEVELS
	469	end do ! naerkind
[253]	470
	471	!-----------------------------------------------------------------------
	472	! Aerosol optical depths
[1483]	473	!-----------------------------------------------------------------------
[253]	474
	475	do iaer=1,naerkind ! a bug was here
	476	do k=0,nlayer-1
	477
	478	pweight=(pplay(ig,L_NLAYRAD-k)-pplev(ig,L_NLAYRAD-k+1))/ &
[1483]	479	(pplev(ig,L_NLAYRAD-k)-pplev(ig,L_NLAYRAD-k+1))
[253]	480	temp=aerosol(ig,L_NLAYRAD-k,iaer)/QREFvis3d(ig,L_NLAYRAD-k,iaer)
[588]	481	tauaero(2k+2,iaer)=max(temppweight,0.d0)
	482	tauaero(2k+3,iaer)=max(temp-tauaero(2k+2,iaer),0.d0)
[1483]	483
[253]	484	end do
	485	! boundary conditions
	486	tauaero(1,iaer) = tauaero(2,iaer)
	487	!tauaero(1,iaer) = 0.
[1988]	488	!JL18 at time of testing, the two above conditions gave the same results bit for bit.
	489
[1483]	490	end do ! naerkind
[253]	491
[1483]	492	! Albedo and Emissivity.
	493	albi=1-emis(ig) ! Long Wave.
	494	DO nw=1,L_NSPECTV ! Short Wave loop.
[1482]	495	albv(nw)=albedo(ig,nw)
[1529]	496	ENDDO
[253]	497
[1483]	498	if ((ngrid.eq.1).and.(global1d)) then ! Fixed zenith angle 'szangle' in 1D simulations w/ globally-averaged sunlight.
[253]	499	acosz = cos(pi*szangle/180.0)
	500	print*,'acosz=',acosz,', szangle=',szangle
	501	else
[1483]	502	acosz=mu0(ig) ! Cosine of sun incident angle : 3D simulations or local 1D simulations using latitude.
[253]	503	endif
	504
[1483]	505	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	506	!!! Note by JL13 : In the following, some indices were changed in the interpolations,
	507	!!! so that the model results are less dependent on the number of layers !
	508	!!!
	509	!!! --- The older versions are commented with the comment !JL13index ---
	510	!!!
	511	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[1016]	512
	513
[253]	514	!-----------------------------------------------------------------------
[1483]	515	! Water vapour (to be generalised for other gases eventually ...)
	516	!-----------------------------------------------------------------------
[253]	517
[305]	518	if(varactive)then
[253]	519
	520	i_var=igcm_h2o_vap
	521	do l=1,nlayer
	522	qvar(2*l) = pq(ig,nlayer+1-l,i_var)
[1016]	523	qvar(2*l+1) = pq(ig,nlayer+1-l,i_var)
	524	!JL13index qvar(2*l+1) = (pq(ig,nlayer+1-l,i_var)+pq(ig,max(nlayer-l,1),i_var))/2
	525	!JL13index ! Average approximation as for temperature...
[253]	526	end do
	527	qvar(1)=qvar(2)
	528
	529	elseif(varfixed)then
	530
[1483]	531	do l=1,nlayer ! Here we will assign fixed water vapour profiles globally.
[253]	532	RH = satval * ((pplay(ig,l)/pplev(ig,1) - 0.02) / 0.98)
	533	if(RH.lt.0.0) RH=0.0
	534
	535	ptemp=pplay(ig,l)
	536	Ttemp=pt(ig,l)
	537	call watersat(Ttemp,ptemp,qsat)
	538
	539	!pq_temp(l) = qsat ! fully saturated everywhere
	540	pq_temp(l) = RH * qsat ! ~realistic profile (e.g. 80% saturation at ground)
	541	end do
	542
	543	do l=1,nlayer
	544	qvar(2*l) = pq_temp(nlayer+1-l)
	545	qvar(2*l+1) = (pq_temp(nlayer+1-l)+pq_temp(max(nlayer-l,1)))/2
	546	end do
[1483]	547
[253]	548	qvar(1)=qvar(2)
	549
	550	! Lowest layer of atmosphere
	551	RH = satval * (1 - 0.02) / 0.98
	552	if(RH.lt.0.0) RH=0.0
	553
[1016]	554	! ptemp = pplev(ig,1)
	555	! Ttemp = tsurf(ig)
	556	! call watersat(Ttemp,ptemp,qsat)
[253]	557
[1308]	558	qvar(2nlayer+1)= RH qsat ! ~realistic profile (e.g. 80% saturation at ground)
[1016]	559
[253]	560	else
	561	do k=1,L_LEVELS
	562	qvar(k) = 1.0D-7
	563	end do
[1483]	564	end if ! varactive/varfixed
[253]	565
[538]	566	if(.not.kastprof)then
[1483]	567	! IMPORTANT: Now convert from kg/kg to mol/mol.
[728]	568	do k=1,L_LEVELS
	569	qvar(k) = qvar(k)/(epsi+qvar(k)*(1.-epsi))
	570	end do
[538]	571	end if
[253]	572
[366]	573	!-----------------------------------------------------------------------
[1483]	574	! kcm mode only !
	575	!-----------------------------------------------------------------------
	576
[305]	577	if(kastprof)then
[1716]	578
	579	if(.not.global1d)then ! garde-fou/safeguard added by MT (to be removed in the future)
	580	write(,) 'You have to fix mu0, '
	581	write(,) 'the cosinus of the solar angle'
	582	stop
	583	endif
	584
[1483]	585	! Initial values equivalent to mugaz.
[305]	586	DO l=1,nlayer
[366]	587	muvarrad(2*l) = mugaz
	588	muvarrad(2*l+1) = mugaz
	589	END DO
	590
[1016]	591	if(ngasmx.gt.1)then
[366]	592
[1016]	593	DO l=1,nlayer
[1483]	594	muvarrad(2*l) = muvar(ig,nlayer+2-l)
[1016]	595	muvarrad(2*l+1) = (muvar(ig,nlayer+2-l) + &
[1483]	596	muvar(ig,max(nlayer+1-l,1)))/2
[1016]	597	END DO
	598
	599	muvarrad(1) = muvarrad(2)
[1483]	600	muvarrad(2*nlayer+1) = muvar(ig,1)
[366]	601
[1016]	602	print*,'Recalculating qvar with VARIABLE epsi for kastprof'
	603	print*,'Assumes that the variable gas is H2O!!!'
	604	print*,'Assumes that there is only one tracer'
[1483]	605
[1016]	606	!i_var=igcm_h2o_vap
	607	i_var=1
[1483]	608
[1016]	609	if(nq.gt.1)then
	610	print*,'Need 1 tracer only to run kcm1d.e'
	611	stop
	612	endif
[1483]	613
[1016]	614	do l=1,nlayer
	615	vtmp(l)=pq(ig,l,i_var)/(epsi+pq(ig,l,i_var)*(1.-epsi))
	616	!vtmp(l)=pq(ig,l,i_var)*muvar(ig,l+1)/mH2O !JL to be changed
	617	end do
[366]	618
[1016]	619	do l=1,nlayer
	620	qvar(2*l) = vtmp(nlayer+1-l)
	621	qvar(2*l+1) = vtmp(nlayer+1-l)
	622	! qvar(2*l+1) = ( vtmp(nlayer+1-l) + vtmp(max(nlayer-l,1)) )/2
	623	end do
	624	qvar(1)=qvar(2)
	625
	626	print*,'Warning: reducing qvar in callcorrk.'
	627	print*,'Temperature profile no longer consistent ', &
[1483]	628	'with saturated H2O. qsat=',satval
	629
[1016]	630	do k=1,L_LEVELS
	631	qvar(k) = qvar(k)*satval
	632	end do
	633
	634	endif
	635	else ! if kastprof
[366]	636	DO l=1,nlayer
[305]	637	muvarrad(2*l) = muvar(ig,nlayer+2-l)
[1016]	638	muvarrad(2*l+1) = (muvar(ig,nlayer+2-l)+muvar(ig,max(nlayer+1-l,1)))/2
[305]	639	END DO
	640
	641	muvarrad(1) = muvarrad(2)
[1308]	642	muvarrad(2*nlayer+1)=muvar(ig,1)
[1483]	643	endif ! if kastprof
[1016]	644
[1483]	645	! Keep values inside limits for which we have radiative transfer coefficients !!!
	646	if(L_REFVAR.gt.1)then ! (there was a bug here)
[253]	647	do k=1,L_LEVELS
	648	if(qvar(k).lt.wrefvar(1))then
	649	qvar(k)=wrefvar(1)+1.0e-8
	650	elseif(qvar(k).gt.wrefvar(L_REFVAR))then
	651	qvar(k)=wrefvar(L_REFVAR)-1.0e-8
	652	endif
	653	end do
	654	endif
	655
	656	!-----------------------------------------------------------------------
	657	! Pressure and temperature
[1483]	658	!-----------------------------------------------------------------------
[253]	659
	660	DO l=1,nlayer
	661	plevrad(2*l) = pplay(ig,nlayer+1-l)/scalep
	662	plevrad(2*l+1) = pplev(ig,nlayer+1-l)/scalep
	663	tlevrad(2*l) = pt(ig,nlayer+1-l)
	664	tlevrad(2*l+1) = (pt(ig,nlayer+1-l)+pt(ig,max(nlayer-l,1)))/2
	665	END DO
	666
[600]	667	plevrad(1) = 0.
[1988]	668	! plevrad(2) = 0. !! JL18 enabling this line puts the radiative top at p=0 which was the idea before, but does not seem to perform best after all.
[253]	669
	670	tlevrad(1) = tlevrad(2)
[1308]	671	tlevrad(2*nlayer+1)=tsurf(ig)
[253]	672
[1988]	673	pmid(1) = pplay(ig,nlayer)/scalep
[1423]	674	pmid(2) = pmid(1)
	675
[253]	676	tmid(1) = tlevrad(2)
[1423]	677	tmid(2) = tmid(1)
	678
	679	DO l=1,L_NLAYRAD-1
	680	tmid(2l+1) = tlevrad(2l+1)
	681	tmid(2l+2) = tlevrad(2l+1)
	682	pmid(2l+1) = plevrad(2l+1)
	683	pmid(2l+2) = plevrad(2l+1)
[253]	684	END DO
[1423]	685	pmid(L_LEVELS) = plevrad(L_LEVELS)
	686	tmid(L_LEVELS) = tlevrad(L_LEVELS)
[253]	687
[1423]	688	!!Alternative interpolation:
	689	! pmid(3) = pmid(1)
	690	! pmid(4) = pmid(1)
	691	! tmid(3) = tmid(1)
	692	! tmid(4) = tmid(1)
	693	! DO l=2,L_NLAYRAD-1
	694	! tmid(2l+1) = tlevrad(2l)
	695	! tmid(2l+2) = tlevrad(2l)
	696	! pmid(2l+1) = plevrad(2l)
	697	! pmid(2l+2) = plevrad(2l)
	698	! END DO
	699	! pmid(L_LEVELS) = plevrad(L_LEVELS-1)
	700	! tmid(L_LEVELS) = tlevrad(L_LEVELS-1)
	701
[1483]	702	! Test for out-of-bounds pressure.
[253]	703	if(plevrad(3).lt.pgasmin)then
	704	print*,'Minimum pressure is outside the radiative'
	705	print*,'transfer kmatrix bounds, exiting.'
	706	call abort
	707	elseif(plevrad(L_LEVELS).gt.pgasmax)then
	708	print*,'Maximum pressure is outside the radiative'
	709	print*,'transfer kmatrix bounds, exiting.'
	710	call abort
	711	endif
	712
[1483]	713	! Test for out-of-bounds temperature.
[253]	714	do k=1,L_LEVELS
	715	if(tlevrad(k).lt.tgasmin)then
	716	print*,'Minimum temperature is outside the radiative'
[1145]	717	print*,'transfer kmatrix bounds'
[858]	718	print*,"k=",k," tlevrad(k)=",tlevrad(k)
	719	print*,"tgasmin=",tgasmin
[1145]	720	if (strictboundcorrk) then
	721	call abort
	722	else
	723	print,'**********************************************'
[1940]	724	print*,'we allow model to continue with tlevrad<tgasmin'
[1145]	725	print*,' ... we assume we know what you are doing ... '
	726	print*,' ... but do not let this happen too often ... '
	727	print,'**********************************************'
[1940]	728	!tlevrad(k)=tgasmin ! Used in the source function !
[1145]	729	endif
[253]	730	elseif(tlevrad(k).gt.tgasmax)then
	731	print*,'Maximum temperature is outside the radiative'
	732	print*,'transfer kmatrix bounds, exiting.'
[1145]	733	print*,"k=",k," tlevrad(k)=",tlevrad(k)
	734	print*,"tgasmax=",tgasmax
	735	if (strictboundcorrk) then
	736	call abort
	737	else
	738	print,'**********************************************'
[1940]	739	print*,'we allow model to continue with tlevrad<tgasmax'
[1145]	740	print*,' ... we assume we know what you are doing ... '
	741	print*,' ... but do not let this happen too often ... '
	742	print,'**********************************************'
[1940]	743	!tlevrad(k)=tgasmax ! Used in the source function !
[1145]	744	endif
[253]	745	endif
	746	enddo
[1016]	747	do k=1,L_NLAYRAD+1
	748	if(tmid(k).lt.tgasmin)then
	749	print*,'Minimum temperature is outside the radiative'
	750	print*,'transfer kmatrix bounds, exiting.'
[1145]	751	print*,"k=",k," tmid(k)=",tmid(k)
[1016]	752	print*,"tgasmin=",tgasmin
[1145]	753	if (strictboundcorrk) then
	754	call abort
	755	else
	756	print,'**********************************************'
[1940]	757	print*,'we allow model to continue but with tmid=tgasmin'
[1145]	758	print*,' ... we assume we know what you are doing ... '
	759	print*,' ... but do not let this happen too often ... '
	760	print,'**********************************************'
	761	tmid(k)=tgasmin
	762	endif
[1016]	763	elseif(tmid(k).gt.tgasmax)then
	764	print*,'Maximum temperature is outside the radiative'
	765	print*,'transfer kmatrix bounds, exiting.'
[1145]	766	print*,"k=",k," tmid(k)=",tmid(k)
	767	print*,"tgasmax=",tgasmax
	768	if (strictboundcorrk) then
	769	call abort
	770	else
	771	print,'**********************************************'
[1940]	772	print*,'we allow model to continue but with tmid=tgasmin'
[1145]	773	print*,' ... we assume we know what you are doing ... '
	774	print*,' ... but do not let this happen too often ... '
	775	print,'**********************************************'
	776	tmid(k)=tgasmax
	777	endif
[1016]	778	endif
	779	enddo
[253]	780
	781	!=======================================================================
[1483]	782	! III. Calling the main radiative transfer subroutines
	783	!=======================================================================
[253]	784
	785
[1483]	786	Cmk= 0.01 * 1.0 / (glat(ig) * mugaz * 1.672621e-27) ! q_main=1.0 assumed.
[1529]	787	glat_ig=glat(ig)
[1194]	788
[253]	789	!-----------------------------------------------------------------------
[1483]	790	! Short Wave Part
	791	!-----------------------------------------------------------------------
[253]	792
[1483]	793	if(fract(ig) .ge. 1.0e-4) then ! Only during daylight.
[787]	794	if((ngrid.eq.1).and.(global1d))then
[253]	795	do nw=1,L_NSPECTV
[1483]	796	stel_fract(nw)= stel(nw)* 0.25 / acosz ! globally averaged = divide by 4, and we correct for solar zenith angle
[253]	797	end do
	798	else
	799	do nw=1,L_NSPECTV
[1161]	800	stel_fract(nw)= stel(nw) * fract(ig)
[253]	801	end do
[1483]	802	endif
	803
[253]	804	call optcv(dtauv,tauv,taucumv,plevrad, &
	805	qxvaer,qsvaer,gvaer,wbarv,cosbv,tauray,tauaero, &
[305]	806	tmid,pmid,taugsurf,qvar,muvarrad)
[253]	807
	808	call sfluxv(dtauv,tauv,taucumv,albv,dwnv,wbarv,cosbv, &
[1781]	809	acosz,stel_fract, &
	810	nfluxtopv,fluxtopvdn,nfluxoutv_nu,nfluxgndv_nu, &
[253]	811	fmnetv,fluxupv,fluxdnv,fzerov,taugsurf)
	812
[1483]	813	else ! During the night, fluxes = 0.
[962]	814	nfluxtopv = 0.0d0
[1529]	815	fluxtopvdn = 0.0d0
[962]	816	nfluxoutv_nu(:) = 0.0d0
	817	nfluxgndv_nu(:) = 0.0d0
[253]	818	do l=1,L_NLAYRAD
[962]	819	fmnetv(l)=0.0d0
	820	fluxupv(l)=0.0d0
	821	fluxdnv(l)=0.0d0
[253]	822	end do
	823	end if
	824
[1482]	825
[1526]	826	! Equivalent Albedo Calculation (for OUTPUT). MT2015
	827	if(fract(ig) .ge. 1.0e-4) then ! equivalent albedo makes sense only during daylight.
	828	surface_stellar_flux=sum(nfluxgndv_nu(1:L_NSPECTV))
	829	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]	830	DO nw=1,L_NSPECTV
	831	albedo_temp(nw)=albedo(ig,nw)*nfluxgndv_nu(nw)
[1526]	832	ENDDO
[1529]	833	albedo_temp(1:L_NSPECTV)=albedo_temp(1:L_NSPECTV)/surface_stellar_flux
[1526]	834	albedo_equivalent(ig)=sum(albedo_temp(1:L_NSPECTV))
	835	else
	836	albedo_equivalent(ig)=0.0 ! Spectrally Integrated Albedo not defined for non-irradiated grid points. So we arbitrary set the equivalent albedo to 0.
	837	endif
[1529]	838	else
	839	albedo_equivalent(ig)=0.0 ! Spectrally Integrated Albedo not defined for non-irradiated grid points. So we arbitrary set the equivalent albedo to 0.
	840	endif
[1482]	841
	842
[253]	843	!-----------------------------------------------------------------------
[1483]	844	! Long Wave Part
	845	!-----------------------------------------------------------------------
[253]	846
	847	call optci(plevrad,tlevrad,dtaui,taucumi, &
	848	qxiaer,qsiaer,giaer,cosbi,wbari,tauaero,tmid,pmid, &
[305]	849	taugsurfi,qvar,muvarrad)
[538]	850
[253]	851	call sfluxi(plevrad,tlevrad,dtaui,taucumi,ubari,albi, &
[1781]	852	wnoi,dwni,cosbi,wbari,nfluxtopi,nfluxtopi_nu, &
[253]	853	fmneti,fluxupi,fluxdni,fluxupi_nu,fzeroi,taugsurfi)
	854
	855	!-----------------------------------------------------------------------
	856	! Transformation of the correlated-k code outputs
	857	! (into dtlw, dtsw, fluxsurf_lw, fluxsurf_sw, fluxtop_lw, fluxtop_sw)
	858
	859	! Flux incident at the top of the atmosphere
[961]	860	fluxtop_dn(ig)=fluxtopvdn
[253]	861
	862	fluxtop_lw(ig) = real(nfluxtopi)
	863	fluxabs_sw(ig) = real(-nfluxtopv)
	864	fluxsurf_lw(ig) = real(fluxdni(L_NLAYRAD))
	865	fluxsurf_sw(ig) = real(fluxdnv(L_NLAYRAD))
[1482]	866
	867	! Flux absorbed by the surface. By MT2015.
	868	fluxsurfabs_sw(ig) = fluxsurf_sw(ig)*(1.-albedo_equivalent(ig))
[253]	869
	870	if(fluxtop_dn(ig).lt.0.0)then
	871	print*,'Achtung! fluxtop_dn has lost the plot!'
	872	print*,'fluxtop_dn=',fluxtop_dn(ig)
	873	print*,'acosz=',acosz
	874	print*,'aerosol=',aerosol(ig,:,:)
	875	print*,'temp= ',pt(ig,:)
	876	print*,'pplay= ',pplay(ig,:)
	877	call abort
	878	endif
	879
	880	! Spectral output, for exoplanet observational comparison
	881	if(specOLR)then
	882	do nw=1,L_NSPECTI
[526]	883	OLR_nu(ig,nw)=nfluxtopi_nu(nw)/DWNI(nw) !JL Normalize to the bandwidth
[253]	884	end do
	885	do nw=1,L_NSPECTV
[366]	886	!GSR_nu(ig,nw)=nfluxgndv_nu(nw)
[526]	887	OSR_nu(ig,nw)=nfluxoutv_nu(nw)/DWNV(nw) !JL Normalize to the bandwidth
[253]	888	end do
	889	endif
	890
	891	! Finally, the heating rates
	892
[586]	893	DO l=2,L_NLAYRAD
	894	dtsw(ig,L_NLAYRAD+1-l)=(fmnetv(l)-fmnetv(l-1)) &
[1194]	895	glat(ig)/(cppscalep(plevrad(2l+1)-plevrad(2*l-1)))
[586]	896	dtlw(ig,L_NLAYRAD+1-l)=(fmneti(l)-fmneti(l-1)) &
[1194]	897	glat(ig)/(cppscalep(plevrad(2l+1)-plevrad(2*l-1)))
[586]	898	END DO
[253]	899
	900	! These are values at top of atmosphere
[586]	901	dtsw(ig,L_NLAYRAD)=(fmnetv(1)-nfluxtopv) &
[1988]	902	glat(ig)/(cppscalep*(plevrad(3)-plevrad(2)))
[586]	903	dtlw(ig,L_NLAYRAD)=(fmneti(1)-nfluxtopi) &
[1988]	904	glat(ig)/(cppscalep*(plevrad(3)-plevrad(2)))
[253]	905
	906
[1483]	907	!-----------------------------------------------------------------------
	908	end do ! End of big loop over every GCM column.
	909	!-----------------------------------------------------------------------
[253]	910
[1483]	911
	912
[253]	913	!-----------------------------------------------------------------------
	914	! Additional diagnostics
[1483]	915	!-----------------------------------------------------------------------
[253]	916
[1483]	917	! IR spectral output, for exoplanet observational comparison
	918	if(lastcall.and.(ngrid.eq.1))then ! could disable the 1D output, they are in the diagfi and diagspec... JL12
[253]	919
[1483]	920	print*,'Saving scalar quantities in surf_vals.out...'
	921	print*,'psurf = ', pplev(1,1),' Pa'
	922	open(116,file='surf_vals.out')
	923	write(116,*) tsurf(1),pplev(1,1),fluxtop_dn(1), &
	924	real(-nfluxtopv),real(nfluxtopi)
	925	close(116)
[253]	926
[526]	927
[1483]	928	! USEFUL COMMENT - Do Not Remove.
	929	!
[526]	930	! if(specOLR)then
	931	! open(117,file='OLRnu.out')
	932	! do nw=1,L_NSPECTI
	933	! write(117,*) OLR_nu(1,nw)
	934	! enddo
	935	! close(117)
	936	!
	937	! open(127,file='OSRnu.out')
	938	! do nw=1,L_NSPECTV
	939	! write(127,*) OSR_nu(1,nw)
	940	! enddo
	941	! close(127)
	942	! endif
[253]	943
[1483]	944	! OLR vs altitude: do it as a .txt file.
	945	OLRz=.false.
	946	if(OLRz)then
	947	print*,'saving IR vertical flux for OLRz...'
	948	open(118,file='OLRz_plevs.out')
	949	open(119,file='OLRz.out')
	950	do l=1,L_NLAYRAD
	951	write(118,) plevrad(2l)
	952	do nw=1,L_NSPECTI
	953	write(119,*) fluxupi_nu(l,nw)
	954	enddo
	955	enddo
	956	close(118)
	957	close(119)
	958	endif
[253]	959
[305]	960	endif
[253]	961
[1483]	962	! See physiq.F for explanations about CLFvarying. This is temporary.
[470]	963	if (lastcall .and. .not.CLFvarying) then
	964	IF( ALLOCATED( gasi ) ) DEALLOCATE( gasi )
	965	IF( ALLOCATED( gasv ) ) DEALLOCATE( gasv )
[1315]	966	!$OMP BARRIER
	967	!$OMP MASTER
[470]	968	IF( ALLOCATED( pgasref ) ) DEALLOCATE( pgasref )
	969	IF( ALLOCATED( tgasref ) ) DEALLOCATE( tgasref )
	970	IF( ALLOCATED( wrefvar ) ) DEALLOCATE( wrefvar )
	971	IF( ALLOCATED( pfgasref ) ) DEALLOCATE( pfgasref )
[1315]	972	!$OMP END MASTER
[1529]	973	!$OMP BARRIER
[861]	974	IF ( ALLOCATED(reffrad)) DEALLOCATE(reffrad)
	975	IF ( ALLOCATED(nueffrad)) DEALLOCATE(nueffrad)
[470]	976	endif
	977
[716]	978
[253]	979	end subroutine callcorrk

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

Download in other formats:

Original Format