Context Navigation

callcorrk.F90 @ 486

Last change on this file since 486 was 486, checked in by rwordsworth, 13 years ago
Variable ice timestep added for ice evolution algorithm. Sedimentation improved to allow consistent dust (L. Kerber). Bug linked to allocatable matrices removed from rcm1d.F. Treatment of initial aerosol radii in callcorrk.F90 improved.
Property svn:executable set to ``*
File size: 30.6 KB

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

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

Download in other formats:

Original Format