Context Navigation

callcorrk.F90 @ 537

Last change on this file since 537 was 526, checked in by jleconte, 13 years ago

13/02/2012 == JL + AS

All outputs are now in netCDF format. Even in 1D (No more G1D)
Clean up of the call to callcorrk when CLFvarying=true
Corrects a bug in writediagspecIR/VI. Output are now in W/m2/cm-1 as a function of the wavenumber in cm-1
Enable writediagspecIR/V to work in the CLFvarying=true case (output now done in Physiq after writediagfi)
Add a simple treatment for the supersaturation of CO2 (see forget et al 2012)
corrects a small bug when no clouds are present in aeropacity

Property svn:executable set to *

File size: 30.8 KB

Rev	Line
[526]	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, &
[526]	6	OLR_nu,OSR_nu, &
	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
	120	real*8 taugsurf(L_NSPECTV,L_NGAUSS-1)
	121	real*8 taugsurfi(L_NSPECTI,L_NGAUSS-1)
	122
[305]	123	real*8 qvar(L_LEVELS) ! mixing ratio of variable component (mol/mol)
[253]	124	REAL pq(ngridmx,nlayer,nq)
	125	REAL qsurf(ngridmx,nqmx) ! tracer on surface (e.g. kg.m-2)
	126	integer nq
	127
	128	! Local aerosol optical properties for each column on RADIATIVE grid
	129	real*8 QXVAER(L_LEVELS+1,L_NSPECTV,naerkind)
	130	real*8 QSVAER(L_LEVELS+1,L_NSPECTV,naerkind)
	131	real*8 GVAER(L_LEVELS+1,L_NSPECTV,naerkind)
	132	real*8 QXIAER(L_LEVELS+1,L_NSPECTI,naerkind)
	133	real*8 QSIAER(L_LEVELS+1,L_NSPECTI,naerkind)
	134	real*8 GIAER(L_LEVELS+1,L_NSPECTI,naerkind)
	135
	136	save qxvaer, qsvaer, gvaer
	137	save qxiaer, qsiaer, giaer
	138	save QREFvis3d, QREFir3d
	139
	140	REAL tau_col(ngrid) ! diagnostic from aeropacity
	141
	142	! Misc.
	143	logical firstcall, lastcall, nantest
	144	real*8 tempv(L_NSPECTV)
	145	real*8 tempi(L_NSPECTI)
	146	real*8 temp,temp1,temp2,pweight
	147	character(len=10) :: tmp1
	148	character(len=10) :: tmp2
	149
	150	! for fixed water vapour profiles
	151	integer i_var
	152	real RH
	153	real*8 pq_temp(nlayer)
	154	real ptemp, Ttemp, qsat
	155
	156	! real(KIND=r8) :: pq_temp(nlayer) ! better F90 way.. DOESNT PORT TO F77!!!
	157
[305]	158	!real ptime, pday
[253]	159	logical OLRz
	160	real*8 NFLUXGNDV_nu(L_NSPECTV)
	161
	162	! for H2O cloud fraction in aeropacity
	163	real cloudfrac(ngridmx,nlayermx)
	164	real totcloudfrac(ngridmx)
	165	logical clearsky
	166
	167	! Allow variations in cp with location
	168	real cpp3D(ngridmx,nlayermx) ! specific heat capacity at const. pressure
	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
	198	if(firstcall) then
	199
	200	call system('rm -f surf_vals_long.out')
	201
	202	!--------------------------------------------------
	203	! Effective radius and variance of the aerosols
	204
[486]	205	do iaer=1,naerkind
[253]	206	! these values will change once the microphysics gets to work
	207	! UNLESS tracer=.false., in which case we should be working with
	208	! a fixed aerosol layer, and be able to define reffrad in a
	209	! .def file. To be improved!
	210
[486]	211	if(iaer.eq.1)then ! CO2 ice
	212	do l=1,nlayer
	213	do ig=1,ngrid
	214	reffrad(ig,l,iaer) = 1.e-4
	215	nueffrad(ig,l,iaer) = 0.1
	216	enddo
	217	enddo
	218	endif
[253]	219
[486]	220	if(iaer.eq.2)then ! H2O ice
	221	do l=1,nlayer
	222	do ig=1,ngrid
	223	reffrad(ig,l,iaer) = 1.e-5
	224	nueffrad(ig,l,iaer) = 0.1
	225	enddo
	226	enddo
	227	endif
	228
	229	if(iaer.eq.3)then ! dust
	230	do l=1,nlayer
	231	do ig=1,ngrid
	232	reffrad(ig,l,iaer) = 1.e-5
	233	nueffrad(ig,l,iaer) = 0.1
	234	enddo
	235	enddo
	236	endif
	237
	238	if(iaer.gt.3)then
	239	print*,'Error in callcorrk, naerkind is too high.'
	240	print*,'The code still needs generalisation to arbitrary'
	241	print*,'aerosol kinds and number.'
	242	call abort
	243	endif
	244
	245	enddo
	246
[374]	247	print*, "callcorrk: Correlated-k data base folder:",trim(datadir)
[253]	248	call getin("corrkdir",corrkdir)
	249	print*, "corrkdir = ",corrkdir
	250	write( tmp1, '(i3)' ) L_NSPECTI
	251	write( tmp2, '(i3)' ) L_NSPECTV
	252	banddir=trim(adjustl(tmp1))//'x'//trim(adjustl(tmp2))
	253	banddir=trim(adjustl(corrkdir))//'/'//trim(adjustl(banddir))
	254
	255	call sugas_corrk ! set up gaseous absorption properties
	256	call setspi ! basic infrared properties
	257	call setspv ! basic visible properties
	258	call suaer_corrk ! set up aerosol optical properties
	259
	260	Cmk= 0.01 * 1.0 / (g * mugaz * 1.672621e-27) ! q_main=1.0 assumed
	261
	262	if((igcm_h2o_vap.eq.0) .and. varactive)then
	263	print*,'varactive in callcorrk but no h2o_vap tracer.'
	264	stop
	265	endif
[526]	266
	267	OLR_nu=0.
	268	OSR_nu=0.
[253]	269
	270	firstcall=.false.
	271
	272	end if
	273
	274	!=======================================================================
	275	! Initialization on every call
	276
	277	do l=1,nlayer
	278	do ig=1,ngrid
	279	do iaer=1,naerkind
	280	nueffrad(ig,l,iaer) = 0.1 ! stays at 0.1
	281	enddo
	282	enddo
	283	enddo
	284
[305]	285	if(kastprof)then
	286	radfixed=.true.
	287	endif
	288
[253]	289	if(radfixed)then
	290	do l=1,nlayer
	291	do ig=1,ngrid
	292	reffrad(ig,l,1) = 5.e-5 ! CO2 ice
	293	enddo
	294	enddo
	295	print*,'CO2 ice particle size = ',reffrad(1,1,1)/1.e-6,' um'
	296	if(naerkind.ge.2)then
	297	do l=1,nlayer
	298	do ig=1,ngrid
[305]	299	!reffrad(ig,l,2) = 2.e-5 ! H2O ice
	300	reffrad(ig,l,2) = 5.e-6 ! H2O ice
[253]	301	enddo
	302	enddo
	303	print*,'H2O ice particle size = ',reffrad(1,1,2)/1.e-6,' um'
	304	endif
	305	if(naerkind.eq.3)then
	306	do l=1,nlayer
	307	do ig=1,ngrid
	308	reffrad(ig,l,naerkind) = 2.e-6 ! dust
	309	enddo
	310	enddo
	311	print*,'Dust particle size = ',reffrad(1,1,naerkind)/1.e-6,' um'
	312	endif
	313	if(naerkind.gt.3)then
	314	print*,'Code not general enough to deal with naerkind > 3 yet.'
	315	call abort
	316	endif
	317
	318	else
	319
	320	maxrad=0.0
[305]	321	!averad=0.0
[253]	322	minrad=1.0
	323	do l=1,nlayer
[305]	324
	325	!masse = (pplev(ig,l) - pplev(ig,l+1))/g
	326
[253]	327	do ig=1,ngrid
[526]	328	if(tracer.and.igcm_co2_ice.gt.0)then
[253]	329	reffrad(ig,l,1) = CBRT( 3*pq(ig,l,igcm_co2_ice)/ &
	330	(4Nmix_co2pi*rho_co2) )
	331	endif
	332	reffrad(ig,l,1) = max(reffrad(ig,l,1),1.e-6)
	333	reffrad(ig,l,1) = min(reffrad(ig,l,1),500.e-6)
	334
[305]	335	!averad = averad + reffrad(ig,l,1)*area(ig)
[253]	336	maxrad = max(reffrad(ig,l,1),maxrad)
	337	minrad = min(reffrad(ig,l,1),minrad)
	338	enddo
	339	enddo
	340	if(igcm_co2_ice.gt.0)then
	341	print*,'Max. CO2 ice particle size = ',maxrad/1.e-6,' um'
	342	print*,'Min. CO2 ice particle size = ',minrad/1.e-6,' um'
	343	endif
	344
	345	if((naerkind.ge.2).and.water)then
	346	maxrad=0.0
	347	minrad=1.0
	348	do l=1,nlayer
	349	do ig=1,ngrid
	350	reffrad(ig,l,2) = CBRT( 3*pq(ig,l,igcm_h2o_ice)/ &
	351	(4Nmix_h2opi*rho_ice) )
	352	reffrad(ig,l,2) = max(reffrad(ig,l,2),1.e-6)
	353	reffrad(ig,l,2) = min(reffrad(ig,l,2),100.e-6)
	354
	355	maxrad = max(reffrad(ig,l,2),maxrad)
	356	minrad = min(reffrad(ig,l,2),minrad)
	357	enddo
	358	enddo
	359	print*,'Max. H2O ice particle size = ',maxrad/1.e-6,' um'
	360	print*,'Min. H2O ice particle size = ',minrad/1.e-6,' um'
	361	endif
	362
	363	if(naerkind.eq.3)then
	364	do l=1,nlayer
	365	do ig=1,ngrid
	366	reffrad(ig,l,naerkind) = 2.e-6 ! dust
	367	enddo
	368	enddo
	369	endif
	370
	371	endif
	372
[305]	373
[253]	374	! how much light we get
	375	fluxtoplanet=0
	376	do nw=1,L_NSPECTV
	377	stel(nw)=stellarf(nw)/(dist_star**2)
	378	fluxtoplanet=fluxtoplanet + stel(nw)
	379	end do
	380
	381	call aeroptproperties(ngrid,nlayer,reffrad,nueffrad, &
	382	QVISsQREF3d,omegaVIS3d,gVIS3d, &
	383	QIRsQREF3d,omegaIR3d,gIR3d, &
	384	QREFvis3d,QREFir3d) ! get 3D aerosol optical properties
	385
	386	call aeropacity(ngrid,nlayer,nq,pplay,pplev,pq,aerosol, &
	387	reffrad,QREFvis3d,QREFir3d, &
	388	tau_col,cloudfrac,totcloudfrac,clearsky) ! get aerosol optical depths
	389
[305]	390
[253]	391	!-----------------------------------------------------------------------
	392	! Starting Big Loop over every GCM column
	393	do ig=1,ngridmx
	394
	395	!=======================================================================
	396	! Transformation of the GCM variables
	397
	398	!-----------------------------------------------------------------------
	399	! Aerosol optical properties Qext, Qscat and g
	400	! The transformation in the vertical is the same as for temperature
	401
	402	! shortwave
	403	do iaer=1,naerkind
	404	DO nw=1,L_NSPECTV
	405	do l=1,nlayermx
	406
	407	temp1=QVISsQREF3d(ig,nlayermx+1-l,nw,iaer) &
	408	*QREFvis3d(ig,nlayermx+1-l,iaer)
	409
	410	temp2=QVISsQREF3d(ig,max(nlayermx-l,1),nw,iaer) &
	411	*QREFvis3d(ig,max(nlayermx-l,1),iaer)
	412
	413	qxvaer(2*l,nw,iaer) = temp1
	414	qxvaer(2*l+1,nw,iaer)=(temp1+temp2)/2
	415
	416	temp1=temp1*omegavis3d(ig,nlayermx+1-l,nw,iaer)
	417	temp2=temp2*omegavis3d(ig,max(nlayermx-l,1),nw,iaer)
	418
	419	qsvaer(2*l,nw,iaer) = temp1
	420	qsvaer(2*l+1,nw,iaer)=(temp1+temp2)/2
	421
	422	temp1=gvis3d(ig,nlayermx+1-l,nw,iaer)
	423	temp2=gvis3d(ig,max(nlayermx-l,1),nw,iaer)
	424
	425	gvaer(2*l,nw,iaer) = temp1
	426	gvaer(2*l+1,nw,iaer)=(temp1+temp2)/2
	427
	428	end do
	429
	430	qxvaer(1,nw,iaer)=qxvaer(2,nw,iaer)
	431	qxvaer(2*nlayermx+1,nw,iaer)=0.
	432
	433	qsvaer(1,nw,iaer)=qsvaer(2,nw,iaer)
	434	qsvaer(2*nlayermx+1,nw,iaer)=0.
	435
	436	gvaer(1,nw,iaer)=gvaer(2,nw,iaer)
	437	gvaer(2*nlayermx+1,nw,iaer)=0.
	438
	439	end do
	440
	441	! longwave
	442	DO nw=1,L_NSPECTI
	443	do l=1,nlayermx
	444
	445	temp1=QIRsQREF3d(ig,nlayermx+1-l,nw,iaer) &
	446	*QREFir3d(ig,nlayermx+1-l,iaer)
	447
	448	temp2=QIRsQREF3d(ig,max(nlayermx-l,1),nw,iaer) &
	449	*QREFir3d(ig,max(nlayermx-l,1),iaer)
	450
	451	qxiaer(2*l,nw,iaer) = temp1
	452	qxiaer(2*l+1,nw,iaer)=(temp1+temp2)/2
	453
	454	temp1=temp1*omegair3d(ig,nlayermx+1-l,nw,iaer)
	455	temp2=temp2*omegair3d(ig,max(nlayermx-l,1),nw,iaer)
	456
	457	qsiaer(2*l,nw,iaer) = temp1
	458	qsiaer(2*l+1,nw,iaer)=(temp1+temp2)/2
	459
	460	temp1=gir3d(ig,nlayermx+1-l,nw,iaer)
	461	temp2=gir3d(ig,max(nlayermx-l,1),nw,iaer)
	462
	463	giaer(2*l,nw,iaer) = temp1
	464	giaer(2*l+1,nw,iaer)=(temp1+temp2)/2
	465
	466	end do
	467
	468	qxiaer(1,nw,iaer)=qxiaer(2,nw,iaer)
	469	qxiaer(2*nlayermx+1,nw,iaer)=0.
	470
	471	qsiaer(1,nw,iaer)=qsiaer(2,nw,iaer)
	472	qsiaer(2*nlayermx+1,nw,iaer)=0.
	473
	474	giaer(1,nw,iaer)=giaer(2,nw,iaer)
	475	giaer(2*nlayermx+1,nw,iaer)=0.
	476
	477	end do
	478	end do
	479
	480	! test / correct for freaky s. s. albedo values
	481	do iaer=1,naerkind
	482	do k=1,L_LEVELS+1
	483
	484	do nw=1,L_NSPECTV
	485	if(qsvaer(k,nw,iaer).gt.1.05*qxvaer(k,nw,iaer))then
	486	print*,'Serious problems with qsvaer values in callcorrk'
	487	call abort
	488	endif
	489	if(qsvaer(k,nw,iaer).gt.qxvaer(k,nw,iaer))then
	490	qsvaer(k,nw,iaer)=qxvaer(k,nw,iaer)
	491	endif
	492	end do
	493
	494	do nw=1,L_NSPECTI
	495	if(qsiaer(k,nw,iaer).gt.1.05*qxiaer(k,nw,iaer))then
	496	print*,'Serious problems with qsiaer values in callcorrk'
	497	call abort
	498	endif
	499	if(qsiaer(k,nw,iaer).gt.qxiaer(k,nw,iaer))then
	500	qsiaer(k,nw,iaer)=qxiaer(k,nw,iaer)
	501	endif
	502	end do
	503
	504	end do
	505	end do
	506
	507	!-----------------------------------------------------------------------
	508	! Aerosol optical depths
	509
	510	do iaer=1,naerkind ! a bug was here
	511	do k=0,nlayer-1
	512
	513	pweight=(pplay(ig,L_NLAYRAD-k)-pplev(ig,L_NLAYRAD-k+1))/ &
	514	(pplev(ig,L_NLAYRAD-k)-pplev(ig,L_NLAYRAD-k+1))
	515
	516	temp=aerosol(ig,L_NLAYRAD-k,iaer)/QREFvis3d(ig,L_NLAYRAD-k,iaer)
	517
	518	tauaero(2k+2,iaer)=max(temppweight,0.0)
	519	tauaero(2k+3,iaer)=max(temp-tauaero(2k+2,iaer),0.0)
	520	!
	521	end do
	522	! boundary conditions
	523	tauaero(1,iaer) = tauaero(2,iaer)
	524	tauaero(L_LEVELS+1,iaer) = tauaero(L_LEVELS,iaer)
	525	!tauaero(1,iaer) = 0.
	526	!tauaero(L_LEVELS+1,iaer) = 0.
	527	end do
	528
	529	! Albedo and emissivity
	530	albi=1-emis(ig) ! longwave
	531	albv=albedo(ig) ! shortwave
	532
	533	if(noradsurf.and.(albv.gt.0.0))then
	534	print*,'For open lower boundary in callcorrk must'
	535	print*,'have surface albedo set to zero!'
	536	call abort
	537	endif
	538
[305]	539	if(ngridmx.eq.1) then ! fixed zenith angle 'szangle' in 1D
[253]	540	acosz = cos(pi*szangle/180.0)
	541	print*,'acosz=',acosz,', szangle=',szangle
	542	else
	543	acosz=mu0(ig) ! cosine of sun incident angle
	544	endif
	545
	546	!-----------------------------------------------------------------------
[305]	547	! Water vapour (to be generalised for other gases eventually)
[253]	548
[305]	549	if(varactive)then
[253]	550
	551	i_var=igcm_h2o_vap
	552	do l=1,nlayer
	553	qvar(2*l) = pq(ig,nlayer+1-l,i_var)
	554	qvar(2*l+1) = (pq(ig,nlayer+1-l,i_var)+pq(ig,max(nlayer-l,1),i_var))/2
	555	! Average approximation as for temperature...
	556	end do
	557	qvar(1)=qvar(2)
[526]	558	! qvar(2*nlayermx+1)=qsurf(ig,i_var) !JL12 not very good to compare kg/kg and kg/m2???
[253]	559
	560	elseif(varfixed)then
	561
	562	do l=1,nlayermx ! here we will assign fixed water vapour profiles globally
	563	RH = satval * ((pplay(ig,l)/pplev(ig,1) - 0.02) / 0.98)
	564	if(RH.lt.0.0) RH=0.0
	565
	566	ptemp=pplay(ig,l)
	567	Ttemp=pt(ig,l)
	568	call watersat(Ttemp,ptemp,qsat)
	569
	570	!pq_temp(l) = qsat ! fully saturated everywhere
	571	pq_temp(l) = RH * qsat ! ~realistic profile (e.g. 80% saturation at ground)
	572	end do
	573
	574	do l=1,nlayer
	575	qvar(2*l) = pq_temp(nlayer+1-l)
	576	qvar(2*l+1) = (pq_temp(nlayer+1-l)+pq_temp(max(nlayer-l,1)))/2
	577	end do
	578	qvar(1)=qvar(2)
	579
	580	! Lowest layer of atmosphere
	581	RH = satval * (1 - 0.02) / 0.98
	582	if(RH.lt.0.0) RH=0.0
	583
	584	ptemp = pplev(ig,1)
	585	Ttemp = tsurf(ig)
	586	call watersat(Ttemp,ptemp,qsat)
[526]	587
[253]	588
	589	!qvar(2*nlayermx+1)=qsat ! fully saturated everywhere
	590	qvar(2nlayermx+1)= RH qsat ! ~realistic profile (e.g. 80% saturation at ground)
	591
[526]	592
	593	!!!!!!!!!!!!!!!!!!!!!!!! JL: for completely constant water profile uncoment the following line
	594	! qvar=0.005
	595
[253]	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)
[526]	768
[253]	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
[526]	799	OLR_nu(ig,nw)=nfluxtopi_nu(nw)/DWNI(nw) !JL Normalize to the bandwidth
[253]	800	end do
	801	do nw=1,L_NSPECTV
[366]	802	!GSR_nu(ig,nw)=nfluxgndv_nu(nw)
[526]	803	OSR_nu(ig,nw)=nfluxoutv_nu(nw)/DWNV(nw) !JL Normalize to the bandwidth
[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
	851
[526]	852	if(lastcall.and.(ngrid.eq.1))then ! could disable the 1D output, they are in the diagfi and diagspec... JL12
	853
[305]	854	print*,'Saving scalar quantities in surf_vals.out...'
	855	print*,'psurf = ', pplev(1,1),' Pa'
[253]	856	open(116,file='surf_vals.out')
	857	write(116,*) tsurf(1),pplev(1,1),fluxtop_dn(1), &
	858	real(-nfluxtopv),real(nfluxtopi)
	859	close(116)
	860
[526]	861	! if(specOLR)then
	862	! open(117,file='OLRnu.out')
	863	! do nw=1,L_NSPECTI
	864	! write(117,*) OLR_nu(1,nw)
	865	! enddo
	866	! close(117)
	867	!
	868	! open(127,file='OSRnu.out')
	869	! do nw=1,L_NSPECTV
	870	! write(127,*) OSR_nu(1,nw)
	871	! enddo
	872	! close(127)
	873	! endif
[253]	874
	875	! OLR vs altitude: do it as a .txt file
	876	OLRz=.false.
	877	if(OLRz)then
	878	print*,'saving IR vertical flux for OLRz...'
	879	open(118,file='OLRz_plevs.out')
	880	open(119,file='OLRz.out')
	881	do l=1,L_NLAYRAD
	882	write(118,) plevrad(2l)
	883	do nw=1,L_NSPECTI
	884	write(119,*) fluxupi_nu(l,nw)
	885	enddo
	886	enddo
	887	close(118)
	888	close(119)
	889	endif
	890
[305]	891	endif
[253]	892
[486]	893	! see physiq.F for explanations about CLFvarying. This is temporary.
[470]	894	if (lastcall .and. .not.CLFvarying) then
	895	IF( ALLOCATED( gasi ) ) DEALLOCATE( gasi )
	896	IF( ALLOCATED( gasv ) ) DEALLOCATE( gasv )
	897	IF( ALLOCATED( pgasref ) ) DEALLOCATE( pgasref )
	898	IF( ALLOCATED( tgasref ) ) DEALLOCATE( tgasref )
	899	IF( ALLOCATED( wrefvar ) ) DEALLOCATE( wrefvar )
	900	IF( ALLOCATED( pfgasref ) ) DEALLOCATE( pfgasref )
	901	endif
	902
[253]	903	end subroutine callcorrk

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

Download in other formats:

Original Format