Context Navigation

callcorrk.F90 @ 485

Last change on this file since 485 was 471, checked in by aslmd, 14 years ago

LMDZ.GENERIC

13/12/2011 == AS

Same spirit as previous commit, but for ngasmx which is now read in gases.def -- before arrays w/ dim ngasmx are allocated dynamically
Allocation is done in su_gases.F90 which is called in inifis
Outside su_gases.F90, very few modifications to the code : the new module "gases_h.F90" simply replaces the old common "gases.h" !
Compiles fine. Tested with debugging options through pgdbg. Runs fine. Exact same results in Early Mars test case.

Property svn:executable set to *

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

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

Download in other formats:

Original Format