Context Navigation

callcorrk.F90 @ 724

Last change on this file since 724 was 716, checked in by rwordsworth, 13 years ago

Mainly updates to radiative transfer and gas management scheme.
Most CIA data now read from standard HITRAN datafiles. For the H2O
continuum, two options have been added: the standard CKD continuum,
and the empirical formula in PPC (Pierrehumbert 2010). Use the toggle
'H2Ocont_simple' in callphys.def to choose.

Note to Martians: I've changed the default values of 'sedimentation' and
'co2cond' in inifis to false. Both these are defined in the standard deftank
callphys.def file, so there shouldn't be any compatibility problems.

Property svn:executable set to *

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

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

Download in other formats:

Original Format