Context Navigation

callcorrk.F90 @ 2099

Last change on this file since 2099 was 2051, checked in by jleconte, 7 years ago

12/12/2018 == JL

Correct a bug from commit 2032 in callcorrk: automatic arrays must be allocated

only if they haven't been before (callcorrk is called twice in physiq when cloud cover is not uniform).

Property svn:executable set to *

File size: 40.4 KB

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

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

Download in other formats:

Original Format