Context Navigation

callcorrk.F90 @ 606

Last change on this file since 606 was 600, checked in by jleconte, 13 years ago

Corrects the computation of planck function at the surface in sfluxi

so that its integral is equal to sigma Tsurf^4.

This ensure that no flux is lost due to:

-truncation of the planck function at high/low wavenumber
-numerical error during first spectral computation of the planck function
-discrepancy between Tsurf and NTS/NTfac in sfluxi

OLR now equal to LW net heating/cooling at equilibrium!
As much as possible, only the value of the stephan boltzmann constant defined in racommon_h (and the

corresponding variable, sigma) should be used. Now done in physics, vdifc and turbdiff.

Property svn:executable set to *

File size: 31.2 KB

Rev	Line
[526]	1	subroutine callcorrk(ngrid,nlayer,pq,nq,qsurf, &
[253]	2	albedo,emis,mu0,pplev,pplay,pt, &
[586]	3	tsurf,fract,dist_star,aerosol,muvar, &
[253]	4	dtlw,dtsw,fluxsurf_lw, &
	5	fluxsurf_sw,fluxtop_lw,fluxabs_sw,fluxtop_dn, &
[538]	6	OLR_nu,OSR_nu, &
[526]	7	reffrad,tau_col,cloudfrac,totcloudfrac, &
[253]	8	clearsky,firstcall,lastcall)
	9
[305]	10
[253]	11	use radinc_h
	12	use radcommon_h
	13	use watercommon_h
[374]	14	use datafile_mod, only: datadir
[253]	15	use ioipsl_getincom
[471]	16	use gases_h
[253]	17
	18	implicit none
	19
	20	!==================================================================
	21	!
	22	! Purpose
	23	! -------
	24	! Solve the radiative transfer using the correlated-k method for
	25	! the gaseous absorption and the Toon et al. (1989) method for
	26	! scatttering due to aerosols.
	27	!
	28	! Authors
	29	! -------
	30	! Emmanuel 01/2001, Forget 09/2001
	31	! Robin Wordsworth (2009)
	32	!
	33	!==================================================================
	34
	35	#include "dimphys.h"
	36	#include "comcstfi.h"
	37	#include "callkeys.h"
	38	#include "tracer.h"
	39
	40	!-----------------------------------------------------------------------
	41	! Declaration of the arguments (INPUT - OUTPUT) on the LMD GCM grid
	42	! Layer #1 is the layer near the ground.
	43	! Layer #nlayermx is the layer at the top.
	44
	45	! INPUT
	46	INTEGER icount
	47	INTEGER ngrid,nlayer
	48	REAL aerosol(ngrid,nlayermx,naerkind) ! aerosol tau (kg/kg)
	49	REAL albedo(ngrid) ! SW albedo
	50	REAL emis(ngrid) ! LW emissivity
	51	REAL pplay(ngrid,nlayermx) ! pres. level in GCM mid of layer
	52	REAL pplev(ngrid,nlayermx+1) ! pres. level at GCM layer boundaries
	53
	54	REAL pt(ngrid,nlayermx) ! air temperature (K)
	55	REAL tsurf(ngrid) ! surface temperature (K)
	56	REAL dist_star,mu0(ngrid) ! distance star-planet (AU)
	57	REAL fract(ngrid) ! fraction of day
	58
	59	! Globally varying aerosol optical properties on GCM grid
	60	! Not needed everywhere so not in radcommon_h
	61	REAL :: QVISsQREF3d(ngridmx,nlayermx,L_NSPECTV,naerkind)
	62	REAL :: omegaVIS3d(ngridmx,nlayermx,L_NSPECTV,naerkind)
	63	REAL :: gVIS3d(ngridmx,nlayermx,L_NSPECTV,naerkind)
	64
	65	REAL :: QIRsQREF3d(ngridmx,nlayermx,L_NSPECTI,naerkind)
	66	REAL :: omegaIR3d(ngridmx,nlayermx,L_NSPECTI,naerkind)
	67	REAL :: gIR3d(ngridmx,nlayermx,L_NSPECTI,naerkind)
	68
	69	REAL :: QREFvis3d(ngridmx,nlayermx,naerkind)
	70	REAL :: QREFir3d(ngridmx,nlayermx,naerkind)
	71
	72	! REAL :: omegaREFvis3d(ngridmx,nlayermx,naerkind)
	73	! REAL :: omegaREFir3d(ngridmx,nlayermx,naerkind) ! not sure of the point of these...
	74
	75	REAL reffrad(ngrid,nlayer,naerkind)
	76	REAL nueffrad(ngrid,nlayer,naerkind)
	77
	78	! OUTPUT
	79	REAL dtsw(ngridmx,nlayermx) ! heating rate (K/s) due to SW
	80	REAL dtlw(ngridmx,nlayermx) ! heating rate (K/s) due to LW
	81	REAL fluxsurf_lw(ngridmx) ! incident LW flux to surf (W/m2)
	82	REAL fluxtop_lw(ngridmx) ! outgoing LW flux to space (W/m2)
	83	REAL fluxsurf_sw(ngridmx) ! incident SW flux to surf (W/m2)
	84	REAL fluxabs_sw(ngridmx) ! SW flux absorbed by planet (W/m2)
	85	REAL fluxtop_dn(ngridmx) ! incident top of atmosphere SW flux (W/m2)
[526]	86	REAL OLR_nu(ngrid,L_NSPECTI)! Outgoing LW radition in each band (Normalized to the band width (W/m2/cm-1)
	87	REAL OSR_nu(ngrid,L_NSPECTV)! Outgoing SW radition in each band (Normalized to the band width (W/m2/cm-1)
[253]	88	!-----------------------------------------------------------------------
	89	! Declaration of the variables required by correlated-k subroutines
	90	! Numbered from top to bottom unlike in the GCM!
	91
	92	REAL*8 tmid(L_LEVELS),pmid(L_LEVELS)
	93	REAL*8 tlevrad(L_LEVELS),plevrad(L_LEVELS)
	94
	95	! Optical values for the optci/cv subroutines
	96	REAL*8 stel(L_NSPECTV),stel_fract(L_NSPECTV)
	97	REAL*8 dtaui(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
	98	REAL*8 dtauv(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
	99	REAL*8 cosbv(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
	100	REAL*8 cosbi(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
	101	REAL*8 wbari(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
	102	REAL*8 wbarv(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
	103	REAL*8 tauv(L_NLEVRAD,L_NSPECTV,L_NGAUSS)
	104	REAL*8 taucumv(L_LEVELS,L_NSPECTV,L_NGAUSS)
	105	REAL*8 taucumi(L_LEVELS,L_NSPECTI,L_NGAUSS)
	106
	107	REAL*8 tauaero(L_LEVELS+1,naerkind)
	108	REAL*8 nfluxtopv,nfluxtopi,nfluxtop
[366]	109	real*8 nfluxoutv_nu(L_NSPECTV) ! outgoing band-resolved VI flux at TOA (W/m2)
	110	real*8 nfluxtopi_nu(L_NSPECTI) ! net band-resolved IR flux at TOA (W/m2)
[253]	111	real*8 fluxupi_nu(L_NLAYRAD,L_NSPECTI) ! for 1D diagnostic
	112	REAL*8 fmneti(L_NLAYRAD),fmnetv(L_NLAYRAD)
	113	REAL*8 fluxupv(L_NLAYRAD),fluxupi(L_NLAYRAD)
	114	REAL*8 fluxdnv(L_NLAYRAD),fluxdni(L_NLAYRAD)
	115	REAL*8 albi,albv,acosz
	116
	117	INTEGER ig,l,k,nw,iaer,irad
	118
	119	real fluxtoplanet
[590]	120	real szangle
	121	logical global1d
	122	save szangle,global1d
[253]	123	real*8 taugsurf(L_NSPECTV,L_NGAUSS-1)
	124	real*8 taugsurfi(L_NSPECTI,L_NGAUSS-1)
	125
[305]	126	real*8 qvar(L_LEVELS) ! mixing ratio of variable component (mol/mol)
[253]	127	REAL pq(ngridmx,nlayer,nq)
	128	REAL qsurf(ngridmx,nqmx) ! tracer on surface (e.g. kg.m-2)
	129	integer nq
	130
	131	! Local aerosol optical properties for each column on RADIATIVE grid
	132	real*8 QXVAER(L_LEVELS+1,L_NSPECTV,naerkind)
	133	real*8 QSVAER(L_LEVELS+1,L_NSPECTV,naerkind)
	134	real*8 GVAER(L_LEVELS+1,L_NSPECTV,naerkind)
	135	real*8 QXIAER(L_LEVELS+1,L_NSPECTI,naerkind)
	136	real*8 QSIAER(L_LEVELS+1,L_NSPECTI,naerkind)
	137	real*8 GIAER(L_LEVELS+1,L_NSPECTI,naerkind)
	138
	139	save qxvaer, qsvaer, gvaer
	140	save qxiaer, qsiaer, giaer
	141	save QREFvis3d, QREFir3d
	142
	143	REAL tau_col(ngrid) ! diagnostic from aeropacity
	144
	145	! Misc.
	146	logical firstcall, lastcall, nantest
	147	real*8 tempv(L_NSPECTV)
	148	real*8 tempi(L_NSPECTI)
	149	real*8 temp,temp1,temp2,pweight
	150	character(len=10) :: tmp1
	151	character(len=10) :: tmp2
	152
	153	! for fixed water vapour profiles
	154	integer i_var
	155	real RH
	156	real*8 pq_temp(nlayer)
	157	real ptemp, Ttemp, qsat
	158
	159	! real(KIND=r8) :: pq_temp(nlayer) ! better F90 way.. DOESNT PORT TO F77!!!
	160
[305]	161	!real ptime, pday
[253]	162	logical OLRz
	163	real*8 NFLUXGNDV_nu(L_NSPECTV)
	164
	165	! for H2O cloud fraction in aeropacity
	166	real cloudfrac(ngridmx,nlayermx)
	167	real totcloudfrac(ngridmx)
	168	logical clearsky
	169
	170	! for weird cloud test
	171	real pqtest(ngridmx,nlayer,nq)
	172
	173	! are particle radii fixed?
	174	logical radfixed
	175	real maxrad, minrad
	176
	177	real CBRT
	178	external CBRT
	179
[366]	180	! included by RW for runaway greenhouse 1D study
[305]	181	real muvar(ngridmx,nlayermx+1)
	182	real vtmp(nlayermx)
	183	REAL*8 muvarrad(L_LEVELS)
	184
[253]	185	radfixed=.false.
	186
	187	!=======================================================================
	188	! Initialization on first call
	189
	190	qxvaer(:,:,:)=0.0
	191	qsvaer(:,:,:)=0.0
	192	gvaer(:,:,:) =0.0
	193
	194	qxiaer(:,:,:)=0.0
	195	qsiaer(:,:,:)=0.0
	196	giaer(:,:,:) =0.0
	197
[538]	198
[253]	199	if(firstcall) then
	200
	201	call system('rm -f surf_vals_long.out')
	202
	203	!--------------------------------------------------
	204	! Effective radius and variance of the aerosols
	205
[486]	206	do iaer=1,naerkind
[253]	207	! these values will change once the microphysics gets to work
	208	! UNLESS tracer=.false., in which case we should be working with
	209	! a fixed aerosol layer, and be able to define reffrad in a
	210	! .def file. To be improved!
	211
[486]	212	if(iaer.eq.1)then ! CO2 ice
	213	do l=1,nlayer
	214	do ig=1,ngrid
	215	reffrad(ig,l,iaer) = 1.e-4
	216	nueffrad(ig,l,iaer) = 0.1
	217	enddo
	218	enddo
	219	endif
[253]	220
[486]	221	if(iaer.eq.2)then ! H2O ice
	222	do l=1,nlayer
	223	do ig=1,ngrid
	224	reffrad(ig,l,iaer) = 1.e-5
	225	nueffrad(ig,l,iaer) = 0.1
	226	enddo
	227	enddo
	228	endif
	229
	230	if(iaer.eq.3)then ! dust
	231	do l=1,nlayer
	232	do ig=1,ngrid
	233	reffrad(ig,l,iaer) = 1.e-5
	234	nueffrad(ig,l,iaer) = 0.1
	235	enddo
	236	enddo
	237	endif
	238
	239	if(iaer.gt.3)then
	240	print*,'Error in callcorrk, naerkind is too high.'
	241	print*,'The code still needs generalisation to arbitrary'
	242	print*,'aerosol kinds and number.'
	243	call abort
	244	endif
	245
	246	enddo
	247
[374]	248	print*, "callcorrk: Correlated-k data base folder:",trim(datadir)
[253]	249	call getin("corrkdir",corrkdir)
	250	print*, "corrkdir = ",corrkdir
	251	write( tmp1, '(i3)' ) L_NSPECTI
	252	write( tmp2, '(i3)' ) L_NSPECTV
	253	banddir=trim(adjustl(tmp1))//'x'//trim(adjustl(tmp2))
	254	banddir=trim(adjustl(corrkdir))//'/'//trim(adjustl(banddir))
	255
	256	call sugas_corrk ! set up gaseous absorption properties
	257	call setspi ! basic infrared properties
	258	call setspv ! basic visible properties
	259	call suaer_corrk ! set up aerosol optical properties
	260
	261	Cmk= 0.01 * 1.0 / (g * mugaz * 1.672621e-27) ! q_main=1.0 assumed
	262
	263	if((igcm_h2o_vap.eq.0) .and. varactive)then
	264	print*,'varactive in callcorrk but no h2o_vap tracer.'
	265	stop
	266	endif
	267
[538]	268	OLR_nu=0.
	269	OSR_nu=0.
	270
[590]	271	if (ngridmx.eq.1) then
	272	PRINT*, 'Simulate global averaged conditions ?'
	273	global1d = .false. ! default value
	274	call getin("global1d",global1d)
	275	write(,) "global1d = ",global1d
	276	if (global1d) then
	277	PRINT *,'Solar Zenith angle (deg.) ?'
	278	PRINT *,'(assumed for averaged solar flux S/4)'
	279	szangle=60.0 ! default value
	280	call getin("szangle",szangle)
	281	write(,) "szangle = ",szangle
	282	endif
	283	endif
	284
[253]	285	firstcall=.false.
	286
	287	end if
	288
	289	!=======================================================================
	290	! Initialization on every call
	291
	292	do l=1,nlayer
	293	do ig=1,ngrid
	294	do iaer=1,naerkind
	295	nueffrad(ig,l,iaer) = 0.1 ! stays at 0.1
	296	enddo
	297	enddo
	298	enddo
	299
[305]	300	if(kastprof)then
	301	radfixed=.true.
	302	endif
	303
[253]	304	if(radfixed)then
	305	do l=1,nlayer
	306	do ig=1,ngrid
	307	reffrad(ig,l,1) = 5.e-5 ! CO2 ice
	308	enddo
	309	enddo
	310	print*,'CO2 ice particle size = ',reffrad(1,1,1)/1.e-6,' um'
	311	if(naerkind.ge.2)then
	312	do l=1,nlayer
	313	do ig=1,ngrid
[305]	314	!reffrad(ig,l,2) = 2.e-5 ! H2O ice
	315	reffrad(ig,l,2) = 5.e-6 ! H2O ice
[253]	316	enddo
	317	enddo
	318	print*,'H2O ice particle size = ',reffrad(1,1,2)/1.e-6,' um'
	319	endif
	320	if(naerkind.eq.3)then
	321	do l=1,nlayer
	322	do ig=1,ngrid
	323	reffrad(ig,l,naerkind) = 2.e-6 ! dust
	324	enddo
	325	enddo
	326	print*,'Dust particle size = ',reffrad(1,1,naerkind)/1.e-6,' um'
	327	endif
	328	if(naerkind.gt.3)then
	329	print*,'Code not general enough to deal with naerkind > 3 yet.'
	330	call abort
	331	endif
	332
	333	else
	334
	335	maxrad=0.0
[305]	336	!averad=0.0
[253]	337	minrad=1.0
	338	do l=1,nlayer
[305]	339
	340	!masse = (pplev(ig,l) - pplev(ig,l+1))/g
	341
[253]	342	do ig=1,ngrid
[596]	343	!if(tracer)then
	344	if(tracer.and.igcm_co2_ice.gt.0)then
[538]	345
	346	if(igcm_co2_ice.lt.1)then
	347	print*,'Tracers but no CO2 ice still seems to be a problem...'
	348	print*,'Aborting in callcorrk.'
	349	stop
	350	endif
	351
[253]	352	reffrad(ig,l,1) = CBRT( 3*pq(ig,l,igcm_co2_ice)/ &
	353	(4Nmix_co2pi*rho_co2) )
	354	endif
	355	reffrad(ig,l,1) = max(reffrad(ig,l,1),1.e-6)
	356	reffrad(ig,l,1) = min(reffrad(ig,l,1),500.e-6)
	357
[305]	358	!averad = averad + reffrad(ig,l,1)*area(ig)
[253]	359	maxrad = max(reffrad(ig,l,1),maxrad)
	360	minrad = min(reffrad(ig,l,1),minrad)
	361	enddo
	362	enddo
	363	if(igcm_co2_ice.gt.0)then
	364	print*,'Max. CO2 ice particle size = ',maxrad/1.e-6,' um'
	365	print*,'Min. CO2 ice particle size = ',minrad/1.e-6,' um'
	366	endif
	367
	368	if((naerkind.ge.2).and.water)then
	369	maxrad=0.0
	370	minrad=1.0
	371	do l=1,nlayer
	372	do ig=1,ngrid
	373	reffrad(ig,l,2) = CBRT( 3*pq(ig,l,igcm_h2o_ice)/ &
	374	(4Nmix_h2opi*rho_ice) )
	375	reffrad(ig,l,2) = max(reffrad(ig,l,2),1.e-6)
	376	reffrad(ig,l,2) = min(reffrad(ig,l,2),100.e-6)
	377
	378	maxrad = max(reffrad(ig,l,2),maxrad)
	379	minrad = min(reffrad(ig,l,2),minrad)
	380	enddo
	381	enddo
	382	print*,'Max. H2O ice particle size = ',maxrad/1.e-6,' um'
	383	print*,'Min. H2O ice particle size = ',minrad/1.e-6,' um'
	384	endif
	385
	386	if(naerkind.eq.3)then
	387	do l=1,nlayer
	388	do ig=1,ngrid
	389	reffrad(ig,l,naerkind) = 2.e-6 ! dust
	390	enddo
	391	enddo
	392	endif
	393
	394	endif
	395
[305]	396
[253]	397	! how much light we get
	398	fluxtoplanet=0
	399	do nw=1,L_NSPECTV
	400	stel(nw)=stellarf(nw)/(dist_star**2)
	401	fluxtoplanet=fluxtoplanet + stel(nw)
	402	end do
	403
	404	call aeroptproperties(ngrid,nlayer,reffrad,nueffrad, &
	405	QVISsQREF3d,omegaVIS3d,gVIS3d, &
	406	QIRsQREF3d,omegaIR3d,gIR3d, &
	407	QREFvis3d,QREFir3d) ! get 3D aerosol optical properties
	408
	409	call aeropacity(ngrid,nlayer,nq,pplay,pplev,pq,aerosol, &
	410	reffrad,QREFvis3d,QREFir3d, &
	411	tau_col,cloudfrac,totcloudfrac,clearsky) ! get aerosol optical depths
	412
[305]	413
[253]	414	!-----------------------------------------------------------------------
	415	! Starting Big Loop over every GCM column
	416	do ig=1,ngridmx
	417
	418	!=======================================================================
	419	! Transformation of the GCM variables
	420
	421	!-----------------------------------------------------------------------
	422	! Aerosol optical properties Qext, Qscat and g
	423	! The transformation in the vertical is the same as for temperature
	424
	425	! shortwave
	426	do iaer=1,naerkind
	427	DO nw=1,L_NSPECTV
	428	do l=1,nlayermx
	429
	430	temp1=QVISsQREF3d(ig,nlayermx+1-l,nw,iaer) &
	431	*QREFvis3d(ig,nlayermx+1-l,iaer)
	432
	433	temp2=QVISsQREF3d(ig,max(nlayermx-l,1),nw,iaer) &
	434	*QREFvis3d(ig,max(nlayermx-l,1),iaer)
	435
	436	qxvaer(2*l,nw,iaer) = temp1
	437	qxvaer(2*l+1,nw,iaer)=(temp1+temp2)/2
	438
	439	temp1=temp1*omegavis3d(ig,nlayermx+1-l,nw,iaer)
	440	temp2=temp2*omegavis3d(ig,max(nlayermx-l,1),nw,iaer)
	441
	442	qsvaer(2*l,nw,iaer) = temp1
	443	qsvaer(2*l+1,nw,iaer)=(temp1+temp2)/2
	444
	445	temp1=gvis3d(ig,nlayermx+1-l,nw,iaer)
	446	temp2=gvis3d(ig,max(nlayermx-l,1),nw,iaer)
	447
	448	gvaer(2*l,nw,iaer) = temp1
	449	gvaer(2*l+1,nw,iaer)=(temp1+temp2)/2
	450
	451	end do
	452
	453	qxvaer(1,nw,iaer)=qxvaer(2,nw,iaer)
	454	qxvaer(2*nlayermx+1,nw,iaer)=0.
	455
	456	qsvaer(1,nw,iaer)=qsvaer(2,nw,iaer)
	457	qsvaer(2*nlayermx+1,nw,iaer)=0.
	458
	459	gvaer(1,nw,iaer)=gvaer(2,nw,iaer)
	460	gvaer(2*nlayermx+1,nw,iaer)=0.
	461
	462	end do
	463
	464	! longwave
	465	DO nw=1,L_NSPECTI
	466	do l=1,nlayermx
	467
	468	temp1=QIRsQREF3d(ig,nlayermx+1-l,nw,iaer) &
	469	*QREFir3d(ig,nlayermx+1-l,iaer)
	470
	471	temp2=QIRsQREF3d(ig,max(nlayermx-l,1),nw,iaer) &
	472	*QREFir3d(ig,max(nlayermx-l,1),iaer)
	473
	474	qxiaer(2*l,nw,iaer) = temp1
	475	qxiaer(2*l+1,nw,iaer)=(temp1+temp2)/2
	476
	477	temp1=temp1*omegair3d(ig,nlayermx+1-l,nw,iaer)
	478	temp2=temp2*omegair3d(ig,max(nlayermx-l,1),nw,iaer)
	479
	480	qsiaer(2*l,nw,iaer) = temp1
	481	qsiaer(2*l+1,nw,iaer)=(temp1+temp2)/2
	482
	483	temp1=gir3d(ig,nlayermx+1-l,nw,iaer)
	484	temp2=gir3d(ig,max(nlayermx-l,1),nw,iaer)
	485
	486	giaer(2*l,nw,iaer) = temp1
	487	giaer(2*l+1,nw,iaer)=(temp1+temp2)/2
	488
	489	end do
	490
	491	qxiaer(1,nw,iaer)=qxiaer(2,nw,iaer)
	492	qxiaer(2*nlayermx+1,nw,iaer)=0.
	493
	494	qsiaer(1,nw,iaer)=qsiaer(2,nw,iaer)
	495	qsiaer(2*nlayermx+1,nw,iaer)=0.
	496
	497	giaer(1,nw,iaer)=giaer(2,nw,iaer)
	498	giaer(2*nlayermx+1,nw,iaer)=0.
	499
	500	end do
	501	end do
	502
	503	! test / correct for freaky s. s. albedo values
	504	do iaer=1,naerkind
	505	do k=1,L_LEVELS+1
	506
	507	do nw=1,L_NSPECTV
	508	if(qsvaer(k,nw,iaer).gt.1.05*qxvaer(k,nw,iaer))then
	509	print*,'Serious problems with qsvaer values in callcorrk'
	510	call abort
	511	endif
	512	if(qsvaer(k,nw,iaer).gt.qxvaer(k,nw,iaer))then
	513	qsvaer(k,nw,iaer)=qxvaer(k,nw,iaer)
	514	endif
	515	end do
	516
	517	do nw=1,L_NSPECTI
	518	if(qsiaer(k,nw,iaer).gt.1.05*qxiaer(k,nw,iaer))then
	519	print*,'Serious problems with qsiaer values in callcorrk'
	520	call abort
	521	endif
	522	if(qsiaer(k,nw,iaer).gt.qxiaer(k,nw,iaer))then
	523	qsiaer(k,nw,iaer)=qxiaer(k,nw,iaer)
	524	endif
	525	end do
	526
	527	end do
	528	end do
	529
	530	!-----------------------------------------------------------------------
	531	! Aerosol optical depths
	532
	533	do iaer=1,naerkind ! a bug was here
	534	do k=0,nlayer-1
	535
	536	pweight=(pplay(ig,L_NLAYRAD-k)-pplev(ig,L_NLAYRAD-k+1))/ &
	537	(pplev(ig,L_NLAYRAD-k)-pplev(ig,L_NLAYRAD-k+1))
	538
	539	temp=aerosol(ig,L_NLAYRAD-k,iaer)/QREFvis3d(ig,L_NLAYRAD-k,iaer)
	540
[588]	541	tauaero(2k+2,iaer)=max(temppweight,0.d0)
	542	tauaero(2k+3,iaer)=max(temp-tauaero(2k+2,iaer),0.d0)
[253]	543	!
	544	end do
	545	! boundary conditions
	546	tauaero(1,iaer) = tauaero(2,iaer)
	547	tauaero(L_LEVELS+1,iaer) = tauaero(L_LEVELS,iaer)
	548	!tauaero(1,iaer) = 0.
	549	!tauaero(L_LEVELS+1,iaer) = 0.
	550	end do
	551
	552	! Albedo and emissivity
	553	albi=1-emis(ig) ! longwave
	554	albv=albedo(ig) ! shortwave
	555
	556	if(noradsurf.and.(albv.gt.0.0))then
	557	print*,'For open lower boundary in callcorrk must'
	558	print*,'have surface albedo set to zero!'
	559	call abort
	560	endif
	561
[590]	562	if ((ngridmx.eq.1).and.(global1d)) then ! fixed zenith angle 'szangle' in 1D simulations w/ globally-averaged sunlight
[253]	563	acosz = cos(pi*szangle/180.0)
	564	print*,'acosz=',acosz,', szangle=',szangle
	565	else
[590]	566	acosz=mu0(ig) ! cosine of sun incident angle : 3D simulations or local 1D simulations using latitude
[253]	567	endif
	568
	569	!-----------------------------------------------------------------------
[305]	570	! Water vapour (to be generalised for other gases eventually)
[253]	571
[305]	572	if(varactive)then
[253]	573
	574	i_var=igcm_h2o_vap
	575	do l=1,nlayer
	576	qvar(2*l) = pq(ig,nlayer+1-l,i_var)
	577	qvar(2*l+1) = (pq(ig,nlayer+1-l,i_var)+pq(ig,max(nlayer-l,1),i_var))/2
	578	! Average approximation as for temperature...
	579	end do
	580	qvar(1)=qvar(2)
	581
	582	elseif(varfixed)then
	583
	584	do l=1,nlayermx ! here we will assign fixed water vapour profiles globally
	585	RH = satval * ((pplay(ig,l)/pplev(ig,1) - 0.02) / 0.98)
	586	if(RH.lt.0.0) RH=0.0
	587
	588	ptemp=pplay(ig,l)
	589	Ttemp=pt(ig,l)
	590	call watersat(Ttemp,ptemp,qsat)
	591
	592	!pq_temp(l) = qsat ! fully saturated everywhere
	593	pq_temp(l) = RH * qsat ! ~realistic profile (e.g. 80% saturation at ground)
	594	end do
	595
	596	do l=1,nlayer
	597	qvar(2*l) = pq_temp(nlayer+1-l)
	598	qvar(2*l+1) = (pq_temp(nlayer+1-l)+pq_temp(max(nlayer-l,1)))/2
	599	end do
	600	qvar(1)=qvar(2)
	601
	602	! Lowest layer of atmosphere
	603	RH = satval * (1 - 0.02) / 0.98
	604	if(RH.lt.0.0) RH=0.0
	605
	606	ptemp = pplev(ig,1)
	607	Ttemp = tsurf(ig)
	608	call watersat(Ttemp,ptemp,qsat)
	609
	610	!qvar(2*nlayermx+1)=qsat ! fully saturated everywhere
	611	qvar(2nlayermx+1)= RH qsat ! ~realistic profile (e.g. 80% saturation at ground)
[538]	612	!qvar=0.005 ! completely constant profile (JL)
[253]	613
	614	else
	615	do k=1,L_LEVELS
	616	qvar(k) = 1.0D-7
	617	end do
	618	end if
	619
[538]	620	if(.not.kastprof)then
[305]	621	! IMPORTANT: Now convert from kg/kg to mol/mol
[253]	622	do k=1,L_LEVELS
	623	qvar(k) = qvar(k)/epsi
	624	end do
[538]	625	end if
[253]	626
[366]	627	!-----------------------------------------------------------------------
	628	! kcm mode only
[305]	629	if(kastprof)then
	630
[486]	631	! initial values equivalent to mugaz
[305]	632	DO l=1,nlayer
[366]	633	muvarrad(2*l) = mugaz
	634	muvarrad(2*l+1) = mugaz
	635	END DO
	636
[486]	637	!do k=1,L_LEVELS
	638	! qvar(k) = 0.0
	639	!end do
	640	!print*,'ASSUMING qH2O=0 EVERYWHERE IN CALLCORRK!'
[366]	641	endif
	642
	643
	644	if(kastprof.and.(ngasmx.gt.1))then
	645
	646	DO l=1,nlayer
[305]	647	muvarrad(2*l) = muvar(ig,nlayer+2-l)
	648	muvarrad(2*l+1) = (muvar(ig,nlayer+2-l) + &
	649	muvar(ig,max(nlayer+1-l,1)))/2
	650	END DO
	651
	652	muvarrad(1) = muvarrad(2)
	653	muvarrad(2*nlayermx+1)=muvar(ig,1)
	654
	655	print*,'Recalculating qvar with VARIABLE epsi for kastprof'
[538]	656	print*,'Assumes that the variable gas is H2O!!!'
	657	print*,'Assumes that there is only one tracer'
	658	!i_var=igcm_h2o_vap
	659	i_var=1
	660	if(nqmx.gt.1)then
	661	print*,'Need 1 tracer only to run kcm1d.e'
	662	stop
	663	endif
[305]	664	do l=1,nlayer
	665	vtmp(l)=pq(ig,l,i_var)*muvar(ig,l+1)/mH2O
	666	end do
	667
	668	do l=1,nlayer
	669	qvar(2*l) = vtmp(nlayer+1-l)
	670	qvar(2*l+1) = ( vtmp(nlayer+1-l) + vtmp(max(nlayer-l,1)) )/2
	671	end do
	672	qvar(1)=qvar(2)
	673
[538]	674	print*,'Warning: reducing qvar in callcorrk.'
	675	print*,'Temperature profile no longer consistent ', &
	676	'with saturated H2O.'
	677	do k=1,L_LEVELS
	678	qvar(k) = qvar(k)*satval
	679	end do
	680
[305]	681	endif
	682
[253]	683	! Keep values inside limits for which we have radiative transfer coefficients
	684	if(L_REFVAR.gt.1)then ! there was a bug here!
	685	do k=1,L_LEVELS
	686	if(qvar(k).lt.wrefvar(1))then
	687	qvar(k)=wrefvar(1)+1.0e-8
	688	elseif(qvar(k).gt.wrefvar(L_REFVAR))then
	689	qvar(k)=wrefvar(L_REFVAR)-1.0e-8
	690	endif
	691	end do
	692	endif
	693
	694	!-----------------------------------------------------------------------
	695	! Pressure and temperature
	696
	697	DO l=1,nlayer
	698	plevrad(2*l) = pplay(ig,nlayer+1-l)/scalep
	699	plevrad(2*l+1) = pplev(ig,nlayer+1-l)/scalep
	700	tlevrad(2*l) = pt(ig,nlayer+1-l)
	701	tlevrad(2*l+1) = (pt(ig,nlayer+1-l)+pt(ig,max(nlayer-l,1)))/2
	702	END DO
	703
[600]	704	plevrad(1) = 0.
[253]	705	plevrad(2) = max(pgasmin,0.0001*plevrad(3))
	706
	707	tlevrad(1) = tlevrad(2)
	708	tlevrad(2*nlayermx+1)=tsurf(ig)
	709
	710	tmid(1) = tlevrad(2)
	711	tmid(2) = tlevrad(2)
	712	pmid(1) = plevrad(2)
	713	pmid(2) = plevrad(2)
	714
	715	DO l=1,L_NLAYRAD-1
	716	tmid(2l+1) = tlevrad(2l+1)
	717	tmid(2l+2) = tlevrad(2l+1)
	718	pmid(2l+1) = plevrad(2l+1)
	719	pmid(2l+2) = plevrad(2l+1)
	720	END DO
	721	pmid(L_LEVELS) = plevrad(L_LEVELS)
	722	tmid(L_LEVELS) = tlevrad(L_LEVELS)
	723
	724	! test for out-of-bounds pressure
	725	if(plevrad(3).lt.pgasmin)then
	726	print*,'Minimum pressure is outside the radiative'
	727	print*,'transfer kmatrix bounds, exiting.'
	728	call abort
	729	elseif(plevrad(L_LEVELS).gt.pgasmax)then
	730	print*,'Maximum pressure is outside the radiative'
	731	print*,'transfer kmatrix bounds, exiting.'
	732	call abort
	733	endif
	734
	735	! test for out-of-bounds temperature
	736	do k=1,L_LEVELS
	737	if(tlevrad(k).lt.tgasmin)then
	738	print*,'Minimum temperature is outside the radiative'
	739	print*,'transfer kmatrix bounds, exiting.'
[486]	740	!print*,'WARNING, OVERRIDING FOR TEST'
	741	call abort
[253]	742	elseif(tlevrad(k).gt.tgasmax)then
	743	print*,'Maximum temperature is outside the radiative'
	744	print*,'transfer kmatrix bounds, exiting.'
[486]	745	!print*,'WARNING, OVERRIDING FOR TEST'
	746	call abort
[253]	747	endif
	748	enddo
	749
	750	!=======================================================================
	751	! Calling the main radiative transfer subroutines
	752
	753
	754	!-----------------------------------------------------------------------
	755	! Shortwave
	756
	757	if(fract(ig) .ge. 1.0e-4) then ! only during daylight!
	758
	759	fluxtoplanet=0.
	760
[590]	761	if((ngridmx.eq.1).and.(global1d))then
[253]	762	do nw=1,L_NSPECTV
	763	stel_fract(nw)= stel(nw) * 0.25 / acosz
	764	fluxtoplanet=fluxtoplanet + stel_fract(nw)
	765	! globally averaged = divide by 4
	766	! but we correct for solar zenith angle
	767	end do
	768	else
	769	do nw=1,L_NSPECTV
	770	stel_fract(nw)= stel(nw) * fract(ig)
	771	fluxtoplanet=fluxtoplanet + stel_fract(nw)
	772	end do
	773	endif
	774
	775	call optcv(dtauv,tauv,taucumv,plevrad, &
	776	qxvaer,qsvaer,gvaer,wbarv,cosbv,tauray,tauaero, &
[305]	777	tmid,pmid,taugsurf,qvar,muvarrad)
[253]	778
	779	call sfluxv(dtauv,tauv,taucumv,albv,dwnv,wbarv,cosbv, &
[366]	780	acosz,stel_fract,gweight, &
	781	nfluxtopv,nfluxoutv_nu,nfluxgndv_nu, &
[253]	782	fmnetv,fluxupv,fluxdnv,fzerov,taugsurf)
	783
	784	else ! during the night, fluxes = 0
[366]	785	nfluxtopv = 0.0
	786	nfluxoutv_nu(:) = 0.0
	787	nfluxgndv_nu(:) = 0.0
[253]	788	do l=1,L_NLAYRAD
	789	fmnetv(l)=0.0
	790	fluxupv(l)=0.0
	791	fluxdnv(l)=0.0
	792	end do
	793	end if
	794
	795	!-----------------------------------------------------------------------
	796	! Longwave
	797
	798	call optci(plevrad,tlevrad,dtaui,taucumi, &
	799	qxiaer,qsiaer,giaer,cosbi,wbari,tauaero,tmid,pmid, &
[305]	800	taugsurfi,qvar,muvarrad)
[538]	801
[253]	802	call sfluxi(plevrad,tlevrad,dtaui,taucumi,ubari,albi, &
	803	wnoi,dwni,cosbi,wbari,gweight,nfluxtopi,nfluxtopi_nu, &
	804	fmneti,fluxupi,fluxdni,fluxupi_nu,fzeroi,taugsurfi)
	805
	806	!-----------------------------------------------------------------------
	807	! Transformation of the correlated-k code outputs
	808	! (into dtlw, dtsw, fluxsurf_lw, fluxsurf_sw, fluxtop_lw, fluxtop_sw)
	809
	810	! Flux incident at the top of the atmosphere
	811	fluxtop_dn(ig)=fluxdnv(1)
	812
	813	fluxtop_lw(ig) = real(nfluxtopi)
	814	fluxabs_sw(ig) = real(-nfluxtopv)
	815	fluxsurf_lw(ig) = real(fluxdni(L_NLAYRAD))
	816	fluxsurf_sw(ig) = real(fluxdnv(L_NLAYRAD))
	817
	818	if(fluxtop_dn(ig).lt.0.0)then
	819	print*,'Achtung! fluxtop_dn has lost the plot!'
	820	print*,'fluxtop_dn=',fluxtop_dn(ig)
	821	print*,'acosz=',acosz
	822	print*,'aerosol=',aerosol(ig,:,:)
	823	print*,'temp= ',pt(ig,:)
	824	print*,'pplay= ',pplay(ig,:)
	825	call abort
	826	endif
	827
	828	! Spectral output, for exoplanet observational comparison
	829	if(specOLR)then
	830	do nw=1,L_NSPECTI
[526]	831	OLR_nu(ig,nw)=nfluxtopi_nu(nw)/DWNI(nw) !JL Normalize to the bandwidth
[253]	832	end do
	833	do nw=1,L_NSPECTV
[366]	834	!GSR_nu(ig,nw)=nfluxgndv_nu(nw)
[526]	835	OSR_nu(ig,nw)=nfluxoutv_nu(nw)/DWNV(nw) !JL Normalize to the bandwidth
[253]	836	end do
	837	endif
	838
	839	! Finally, the heating rates
	840
[586]	841	DO l=2,L_NLAYRAD
	842	dtsw(ig,L_NLAYRAD+1-l)=(fmnetv(l)-fmnetv(l-1)) &
	843	g/(cppscalep(plevrad(2l+1)-plevrad(2*l-1)))
	844	dtlw(ig,L_NLAYRAD+1-l)=(fmneti(l)-fmneti(l-1)) &
	845	g/(cppscalep(plevrad(2l+1)-plevrad(2*l-1)))
	846	END DO
[253]	847
	848	! These are values at top of atmosphere
[586]	849	dtsw(ig,L_NLAYRAD)=(fmnetv(1)-nfluxtopv) &
	850	g/(cppscalep*(plevrad(3)-plevrad(1)))
	851	dtlw(ig,L_NLAYRAD)=(fmneti(1)-nfluxtopi) &
	852	g/(cppscalep*(plevrad(3)-plevrad(1)))
[253]	853
	854	! ---------------------------------------------------------------
	855	end do ! end of big loop over every GCM column (ig = 1:ngrid)
	856
	857
	858	!-----------------------------------------------------------------------
	859	! Additional diagnostics
	860
	861	! IR spectral output, for exoplanet observational comparison
	862
	863
[526]	864	if(lastcall.and.(ngrid.eq.1))then ! could disable the 1D output, they are in the diagfi and diagspec... JL12
	865
[305]	866	print*,'Saving scalar quantities in surf_vals.out...'
	867	print*,'psurf = ', pplev(1,1),' Pa'
[253]	868	open(116,file='surf_vals.out')
	869	write(116,*) tsurf(1),pplev(1,1),fluxtop_dn(1), &
	870	real(-nfluxtopv),real(nfluxtopi)
	871	close(116)
	872
[538]	873	! I am useful, please don`t remove me!
[526]	874	! if(specOLR)then
	875	! open(117,file='OLRnu.out')
	876	! do nw=1,L_NSPECTI
	877	! write(117,*) OLR_nu(1,nw)
	878	! enddo
	879	! close(117)
	880	!
	881	! open(127,file='OSRnu.out')
	882	! do nw=1,L_NSPECTV
	883	! write(127,*) OSR_nu(1,nw)
	884	! enddo
	885	! close(127)
	886	! endif
[253]	887
	888	! OLR vs altitude: do it as a .txt file
	889	OLRz=.false.
	890	if(OLRz)then
	891	print*,'saving IR vertical flux for OLRz...'
	892	open(118,file='OLRz_plevs.out')
	893	open(119,file='OLRz.out')
	894	do l=1,L_NLAYRAD
	895	write(118,) plevrad(2l)
	896	do nw=1,L_NSPECTI
	897	write(119,*) fluxupi_nu(l,nw)
	898	enddo
	899	enddo
	900	close(118)
	901	close(119)
	902	endif
	903
[305]	904	endif
[253]	905
[486]	906	! see physiq.F for explanations about CLFvarying. This is temporary.
[470]	907	if (lastcall .and. .not.CLFvarying) then
	908	IF( ALLOCATED( gasi ) ) DEALLOCATE( gasi )
	909	IF( ALLOCATED( gasv ) ) DEALLOCATE( gasv )
	910	IF( ALLOCATED( pgasref ) ) DEALLOCATE( pgasref )
	911	IF( ALLOCATED( tgasref ) ) DEALLOCATE( tgasref )
	912	IF( ALLOCATED( wrefvar ) ) DEALLOCATE( wrefvar )
	913	IF( ALLOCATED( pfgasref ) ) DEALLOCATE( pfgasref )
	914	endif
	915
[253]	916	end subroutine callcorrk

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

Download in other formats:

Original Format