Context Navigation

callcorrk.F90 @ 726

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

17/07/2012 == JL for LK

Generalization of aerosol scheme:
- any number of aerosols can be used and id numbers are determined consistently by the code. Aerosol order not important anymore.
- addition of a module with the id numbers for aerosols (aerosol_mod.F90).
- initialization of aerosols id numbers in iniaerosol.F90
- compile with -s x where x *must* be equal to the number of aerosols turned on in callphys.def (either by a flag or by dusttau>0 for dust). => may have to erase object files when compiling with s option for the first time.
For no aerosols, run with aeroco2=.true. and aerofixco2=.true (the default distribution for fixed co2

aerosols is 1.e-9; can be changed in aeropacity).

If starting from an old start file, recreate start file with the q=0 option in newstart.e.
update callphys.def with aeroXXX and aerofixXXX options (only XXX=co2,h2o supported for

now). Dust is activated by setting dusttau>0. See the early mars case in deftank.

To add other aerosols, see Laura Kerber.

Property svn:executable set to *

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

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

Download in other formats:

Original Format