Context Navigation

callcorrk.F90 @ 2613

Last change on this file since 2613 was 2543, checked in by aslmd, 3 years ago

Generic GCM:

Adding k-coefficients mixing on the fly
Working with MordernTrac?

JVO + YJ

Property svn:executable set to *

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

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

Download in other formats:

Original Format