Context Navigation

callcorrk.F90 @ 1115

Last change on this file since 1115 was 1026, checked in by sglmd, 12 years ago
Added a flexible, 2-layer aerosol scenario (initially for Saturn, but can be simply tuned). Called aeroback2lay.
Property svn:executable set to ``*
File size: 31.6 KB

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

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

Download in other formats:

Original Format