Context Navigation

callcorrk.F90 @ 1680

Last change on this file since 1680 was 1648, checked in by jvatant, 8 years ago

Modifications to custom radiative transfer to Titan
+ Enables an altitude dependant gfrac for CIA computations

-> many radical changes in su_gases and co ..
-> read vertical CH4 profile with call_profilgases
-> Now you need a 'profile.def' that I will add in the deftank

+ Added interpolate CIA routines for CH4
+ Added temporary mean aerosol profile opacity routine (disr_haze)

Property svn:executable set to *

File size: 31.0 KB

Rev	Line
[526]	1	subroutine callcorrk(ngrid,nlayer,pq,nq,qsurf, &
[1482]	2	albedo,albedo_equivalent,emis,mu0,pplev,pplay,pt, &
[1648]	3	tsurf,fract,dist_star,aerosol, &
[253]	4	dtlw,dtsw,fluxsurf_lw, &
[1482]	5	fluxsurf_sw,fluxsurfabs_sw,fluxtop_lw, &
	6	fluxabs_sw,fluxtop_dn, &
[538]	7	OLR_nu,OSR_nu, &
[1647]	8	tau_col,firstcall,lastcall)
[253]	9
	10	use radinc_h
	11	use radcommon_h
[374]	12	use datafile_mod, only: datadir
[1521]	13	use ioipsl_getin_p_mod, only: getin_p
[471]	14	use gases_h
[1647]	15	use radii_mod, only : su_aer_radii,back2lay_reffrad
	16	use aerosol_mod, only : iaero_back2lay
[787]	17	USE tracer_h
[1384]	18	use comcstfi_mod, only: pi, mugaz, cpp
[1648]	19	use callkeys_mod, only: diurnal,tracer,nosurf, &
[1647]	20	strictboundcorrk,specOLR
[787]	21
[253]	22	implicit none
	23
	24	!==================================================================
	25	!
	26	! Purpose
	27	! -------
	28	! Solve the radiative transfer using the correlated-k method for
	29	! the gaseous absorption and the Toon et al. (1989) method for
	30	! scatttering due to aerosols.
	31	!
	32	! Authors
	33	! -------
	34	! Emmanuel 01/2001, Forget 09/2001
	35	! Robin Wordsworth (2009)
	36	!
	37	!==================================================================
	38
	39	!-----------------------------------------------------------------------
	40	! Declaration of the arguments (INPUT - OUTPUT) on the LMD GCM grid
	41	! Layer #1 is the layer near the ground.
[1308]	42	! Layer #nlayer is the layer at the top.
[1483]	43	!-----------------------------------------------------------------------
[253]	44
[1483]	45
	46	! INPUT
	47	INTEGER,INTENT(IN) :: ngrid ! Number of atmospheric columns.
	48	INTEGER,INTENT(IN) :: nlayer ! Number of atmospheric layers.
	49	REAL,INTENT(IN) :: pq(ngrid,nlayer,nq) ! Tracers (kg/kg_of_air).
	50	INTEGER,INTENT(IN) :: nq ! Number of tracers.
	51	REAL,INTENT(IN) :: qsurf(ngrid,nq) ! Tracers on surface (kg.m-2).
	52	REAL,INTENT(IN) :: albedo(ngrid,L_NSPECTV) ! Spectral Short Wavelengths Albedo. By MT2015
	53	REAL,INTENT(IN) :: emis(ngrid) ! Long Wave emissivity.
	54	REAL,INTENT(IN) :: mu0(ngrid) ! Cosine of sun incident angle.
	55	REAL,INTENT(IN) :: pplev(ngrid,nlayer+1) ! Inter-layer pressure (Pa).
	56	REAL,INTENT(IN) :: pplay(ngrid,nlayer) ! Mid-layer pressure (Pa).
	57	REAL,INTENT(IN) :: pt(ngrid,nlayer) ! Air temperature (K).
	58	REAL,INTENT(IN) :: tsurf(ngrid) ! Surface temperature (K).
	59	REAL,INTENT(IN) :: fract(ngrid) ! Fraction of day.
	60	REAL,INTENT(IN) :: dist_star ! Distance star-planet (AU).
	61	logical,intent(in) :: firstcall ! Signals first call to physics.
	62	logical,intent(in) :: lastcall ! Signals last call to physics.
	63
	64	! OUTPUT
	65	REAL,INTENT(OUT) :: aerosol(ngrid,nlayer,naerkind) ! Aerosol tau (kg/kg).
	66	REAL,INTENT(OUT) :: dtlw(ngrid,nlayer) ! Heating rate (K/s) due to LW radiation.
	67	REAL,INTENT(OUT) :: dtsw(ngrid,nlayer) ! Heating rate (K/s) due to SW radiation.
	68	REAL,INTENT(OUT) :: fluxsurf_lw(ngrid) ! Incident LW flux to surf (W/m2).
	69	REAL,INTENT(OUT) :: fluxsurf_sw(ngrid) ! Incident SW flux to surf (W/m2)
	70	REAL,INTENT(OUT) :: fluxsurfabs_sw(ngrid) ! Absorbed SW flux by the surface (W/m2). By MT2015.
	71	REAL,INTENT(OUT) :: fluxtop_lw(ngrid) ! Outgoing LW flux to space (W/m2).
	72	REAL,INTENT(OUT) :: fluxabs_sw(ngrid) ! SW flux absorbed by the planet (W/m2).
	73	REAL,INTENT(OUT) :: fluxtop_dn(ngrid) ! Incident top of atmosphere SW flux (W/m2).
	74	REAL,INTENT(OUT) :: OLR_nu(ngrid,L_NSPECTI) ! Outgoing LW radition in each band (Normalized to the band width (W/m2/cm-1).
	75	REAL,INTENT(OUT) :: OSR_nu(ngrid,L_NSPECTV) ! Outgoing SW radition in each band (Normalized to the band width (W/m2/cm-1).
	76	REAL,INTENT(OUT) :: tau_col(ngrid) ! Diagnostic from aeropacity.
[1647]	77	REAL,INTENT(OUT) :: albedo_equivalent(ngrid) ! Spectrally Integrated Albedo. For Diagnostic. By MT2015
[1483]	78
	79
	80	! Globally varying aerosol optical properties on GCM grid ; not needed everywhere so not in radcommon_h.
[1308]	81	REAL :: QVISsQREF3d(ngrid,nlayer,L_NSPECTV,naerkind)
	82	REAL :: omegaVIS3d(ngrid,nlayer,L_NSPECTV,naerkind)
	83	REAL :: gVIS3d(ngrid,nlayer,L_NSPECTV,naerkind)
	84	REAL :: QIRsQREF3d(ngrid,nlayer,L_NSPECTI,naerkind)
	85	REAL :: omegaIR3d(ngrid,nlayer,L_NSPECTI,naerkind)
	86	REAL :: gIR3d(ngrid,nlayer,L_NSPECTI,naerkind)
[253]	87
[1308]	88	! REAL :: omegaREFvis3d(ngrid,nlayer,naerkind)
	89	! REAL :: omegaREFir3d(ngrid,nlayer,naerkind) ! not sure of the point of these...
[253]	90
[1483]	91	REAL,ALLOCATABLE,SAVE :: reffrad(:,:,:) ! aerosol effective radius (m)
[858]	92	REAL,ALLOCATABLE,SAVE :: nueffrad(:,:,:) ! aerosol effective variance
[1315]	93	!$OMP THREADPRIVATE(reffrad,nueffrad)
[253]	94
	95	!-----------------------------------------------------------------------
	96	! Declaration of the variables required by correlated-k subroutines
[1483]	97	! Numbered from top to bottom (unlike in the GCM)
	98	!-----------------------------------------------------------------------
[253]	99
	100	REAL*8 tmid(L_LEVELS),pmid(L_LEVELS)
	101	REAL*8 tlevrad(L_LEVELS),plevrad(L_LEVELS)
	102
[1483]	103	! Optical values for the optci/cv subroutines
[253]	104	REAL*8 stel(L_NSPECTV),stel_fract(L_NSPECTV)
	105	REAL*8 dtaui(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
	106	REAL*8 dtauv(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
	107	REAL*8 cosbv(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
	108	REAL*8 cosbi(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
	109	REAL*8 wbari(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
	110	REAL*8 wbarv(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
	111	REAL*8 tauv(L_NLEVRAD,L_NSPECTV,L_NGAUSS)
	112	REAL*8 taucumv(L_LEVELS,L_NSPECTV,L_NGAUSS)
	113	REAL*8 taucumi(L_LEVELS,L_NSPECTI,L_NGAUSS)
	114
	115	REAL*8 tauaero(L_LEVELS+1,naerkind)
[961]	116	REAL*8 nfluxtopv,nfluxtopi,nfluxtop,fluxtopvdn
[1483]	117	REAL*8 nfluxoutv_nu(L_NSPECTV) ! Outgoing band-resolved VI flux at TOA (W/m2).
	118	REAL*8 nfluxtopi_nu(L_NSPECTI) ! Net band-resolved IR flux at TOA (W/m2).
	119	REAL*8 fluxupi_nu(L_NLAYRAD,L_NSPECTI) ! For 1D diagnostic.
[253]	120	REAL*8 fmneti(L_NLAYRAD),fmnetv(L_NLAYRAD)
	121	REAL*8 fluxupv(L_NLAYRAD),fluxupi(L_NLAYRAD)
	122	REAL*8 fluxdnv(L_NLAYRAD),fluxdni(L_NLAYRAD)
[1482]	123	REAL*8 albi,acosz
[1483]	124	REAL*8 albv(L_NSPECTV) ! Spectral Visible Albedo.
[253]	125
[1483]	126	INTEGER ig,l,k,nw,iaer
[253]	127
[590]	128	real szangle
	129	logical global1d
	130	save szangle,global1d
[1315]	131	!$OMP THREADPRIVATE(szangle,global1d)
[253]	132	real*8 taugsurf(L_NSPECTV,L_NGAUSS-1)
	133	real*8 taugsurfi(L_NSPECTI,L_NGAUSS-1)
	134
[1483]	135	! Local aerosol optical properties for each column on RADIATIVE grid.
[1529]	136	real*8,save,allocatable :: QXVAER(:,:,:)
	137	real*8,save,allocatable :: QSVAER(:,:,:)
	138	real*8,save,allocatable :: GVAER(:,:,:)
	139	real*8,save,allocatable :: QXIAER(:,:,:)
	140	real*8,save,allocatable :: QSIAER(:,:,:)
	141	real*8,save,allocatable :: GIAER(:,:,:)
[253]	142
[787]	143	real, dimension(:,:,:), save, allocatable :: QREFvis3d
	144	real, dimension(:,:,:), save, allocatable :: QREFir3d
[1315]	145	!$OMP THREADPRIVATE(QXVAER,QSVAER,GVAER,QXIAER,QSIAER,GIAER,QREFvis3d,QREFir3d)
[787]	146
[253]	147
[1483]	148	! Miscellaneous :
[253]	149	real*8 temp,temp1,temp2,pweight
	150	character(len=10) :: tmp1
	151	character(len=10) :: tmp2
	152
	153	logical OLRz
	154	real*8 NFLUXGNDV_nu(L_NSPECTV)
[1648]	155
[1483]	156	! Included by MT for albedo calculations.
[1482]	157	REAL*8 albedo_temp(L_NSPECTV) ! For equivalent albedo calculation.
[1526]	158	REAL*8 surface_stellar_flux ! Stellar flux reaching the surface. Useful for equivalent albedo calculation.
[305]	159
[1483]	160
[726]	161	!===============================================================
[1483]	162	! I.a Initialization on first call
	163	!===============================================================
[253]	164
[1483]	165
[1529]	166	if(firstcall) then
[253]	167
[1529]	168	! test on allocated necessary because of CLFvarying (two calls to callcorrk in physiq)
	169	if(.not.allocated(QXVAER)) allocate(QXVAER(L_LEVELS+1,L_NSPECTV,naerkind))
	170	if(.not.allocated(QSVAER)) allocate(QSVAER(L_LEVELS+1,L_NSPECTV,naerkind))
	171	if(.not.allocated(GVAER)) allocate(GVAER(L_LEVELS+1,L_NSPECTV,naerkind))
	172	if(.not.allocated(QXIAER)) allocate(QXIAER(L_LEVELS+1,L_NSPECTI,naerkind))
	173	if(.not.allocated(QSIAER)) allocate(QSIAER(L_LEVELS+1,L_NSPECTI,naerkind))
	174	if(.not.allocated(GIAER)) allocate(GIAER(L_LEVELS+1,L_NSPECTI,naerkind))
[253]	175
[1483]	176	!!! ALLOCATED instances are necessary because of CLFvarying (strategy to call callcorrk twice in physiq...)
[1308]	177	IF(.not.ALLOCATED(QREFvis3d)) ALLOCATE(QREFvis3d(ngrid,nlayer,naerkind))
	178	IF(.not.ALLOCATED(QREFir3d)) ALLOCATE(QREFir3d(ngrid,nlayer,naerkind))
[861]	179	! Effective radius and variance of the aerosols
	180	IF(.not.ALLOCATED(reffrad)) allocate(reffrad(ngrid,nlayer,naerkind))
	181	IF(.not.ALLOCATED(nueffrad)) allocate(nueffrad(ngrid,nlayer,naerkind))
[787]	182
[253]	183	call system('rm -f surf_vals_long.out')
	184
[728]	185	if(naerkind.gt.4)then
	186	print*,'Code not general enough to deal with naerkind > 4 yet.'
	187	call abort
	188	endif
[1308]	189	call su_aer_radii(ngrid,nlayer,reffrad,nueffrad)
[1483]	190
[1529]	191
[728]	192	!--------------------------------------------------
[1483]	193	! Set up correlated k
	194	!--------------------------------------------------
	195
	196
[374]	197	print*, "callcorrk: Correlated-k data base folder:",trim(datadir)
[1315]	198	call getin_p("corrkdir",corrkdir)
[253]	199	print*, "corrkdir = ",corrkdir
	200	write( tmp1, '(i3)' ) L_NSPECTI
	201	write( tmp2, '(i3)' ) L_NSPECTV
	202	banddir=trim(adjustl(tmp1))//'x'//trim(adjustl(tmp2))
	203	banddir=trim(adjustl(corrkdir))//'/'//trim(adjustl(banddir))
	204
[1483]	205	call setspi ! Basic infrared properties.
	206	call setspv ! Basic visible properties.
	207	call sugas_corrk ! Set up gaseous absorption properties.
	208	call suaer_corrk ! Set up aerosol optical properties.
[1498]	209
[253]	210
[716]	211	OLR_nu(:,:) = 0.
	212	OSR_nu(:,:) = 0.
[538]	213
[787]	214	if (ngrid.eq.1) then
[1483]	215	PRINT*, 'Simulate global averaged conditions ?'
	216	global1d = .false. ! default value
	217	call getin_p("global1d",global1d)
	218	write(,) "global1d = ",global1d
	219
	220	! Test of incompatibility : if global1d is true, there should not be any diurnal cycle.
	221	if (global1d.and.diurnal) then
	222	print*,'if global1d is true, diurnal must be set to false'
	223	stop
	224	endif
[622]	225
[1483]	226	if (global1d) then
	227	PRINT *,'Solar Zenith angle (deg.) ?'
	228	PRINT *,'(assumed for averaged solar flux S/4)'
	229	szangle=60.0 ! default value
	230	call getin_p("szangle",szangle)
	231	write(,) "szangle = ",szangle
	232	endif
[590]	233	endif
	234
[858]	235	end if ! of if (firstcall)
[253]	236
	237	!=======================================================================
[1483]	238	! I.b Initialization on every call
	239	!=======================================================================
	240
[1529]	241	qxvaer(:,:,:)=0.0
	242	qsvaer(:,:,:)=0.0
	243	gvaer(:,:,:) =0.0
	244
	245	qxiaer(:,:,:)=0.0
	246	qsiaer(:,:,:)=0.0
	247	giaer(:,:,:) =0.0
	248
[728]	249	!--------------------------------------------------
	250	! Effective radius and variance of the aerosols
[1483]	251	!--------------------------------------------------
	252
[726]	253	do iaer=1,naerkind
[1647]	254
[1026]	255	if(iaer.eq.iaero_back2lay)then
[1529]	256	call back2lay_reffrad(ngrid,reffrad(1,1,iaero_back2lay),nlayer,pplev)
[1026]	257	endif
[1483]	258
	259	end do !iaer=1,naerkind.
[253]	260
	261
[1483]	262	! How much light do we get ?
[253]	263	do nw=1,L_NSPECTV
	264	stel(nw)=stellarf(nw)/(dist_star**2)
	265	end do
	266
[1483]	267	! Get 3D aerosol optical properties.
[253]	268	call aeroptproperties(ngrid,nlayer,reffrad,nueffrad, &
	269	QVISsQREF3d,omegaVIS3d,gVIS3d, &
	270	QIRsQREF3d,omegaIR3d,gIR3d, &
[1483]	271	QREFvis3d,QREFir3d)
[253]	272
[1483]	273	! Get aerosol optical depths.
[253]	274	call aeropacity(ngrid,nlayer,nq,pplay,pplev,pq,aerosol, &
	275	reffrad,QREFvis3d,QREFir3d, &
[1647]	276	tau_col)
[1529]	277
[1483]	278
	279
	280	!-----------------------------------------------------------------------
	281	do ig=1,ngrid ! Starting Big Loop over every GCM column
[253]	282	!-----------------------------------------------------------------------
	283
[1483]	284
[253]	285	!=======================================================================
[1483]	286	! II. Transformation of the GCM variables
	287	!=======================================================================
[253]	288
[1483]	289
[253]	290	!-----------------------------------------------------------------------
[1483]	291	! Aerosol optical properties Qext, Qscat and g.
	292	! The transformation in the vertical is the same as for temperature.
	293	!-----------------------------------------------------------------------
[253]	294
[1483]	295
[253]	296	do iaer=1,naerkind
[1483]	297	! Shortwave.
	298	do nw=1,L_NSPECTV
	299
[1308]	300	do l=1,nlayer
[253]	301
[1308]	302	temp1=QVISsQREF3d(ig,nlayer+1-l,nw,iaer) &
	303	*QREFvis3d(ig,nlayer+1-l,iaer)
[253]	304
[1308]	305	temp2=QVISsQREF3d(ig,max(nlayer-l,1),nw,iaer) &
	306	*QREFvis3d(ig,max(nlayer-l,1),iaer)
[253]	307
	308	qxvaer(2*l,nw,iaer) = temp1
	309	qxvaer(2*l+1,nw,iaer)=(temp1+temp2)/2
	310
[1308]	311	temp1=temp1*omegavis3d(ig,nlayer+1-l,nw,iaer)
	312	temp2=temp2*omegavis3d(ig,max(nlayer-l,1),nw,iaer)
[253]	313
	314	qsvaer(2*l,nw,iaer) = temp1
	315	qsvaer(2*l+1,nw,iaer)=(temp1+temp2)/2
	316
[1308]	317	temp1=gvis3d(ig,nlayer+1-l,nw,iaer)
	318	temp2=gvis3d(ig,max(nlayer-l,1),nw,iaer)
[253]	319
	320	gvaer(2*l,nw,iaer) = temp1
	321	gvaer(2*l+1,nw,iaer)=(temp1+temp2)/2
	322
[1483]	323	end do ! nlayer
[253]	324
	325	qxvaer(1,nw,iaer)=qxvaer(2,nw,iaer)
[1308]	326	qxvaer(2*nlayer+1,nw,iaer)=0.
[253]	327
	328	qsvaer(1,nw,iaer)=qsvaer(2,nw,iaer)
[1308]	329	qsvaer(2*nlayer+1,nw,iaer)=0.
[253]	330
	331	gvaer(1,nw,iaer)=gvaer(2,nw,iaer)
[1308]	332	gvaer(2*nlayer+1,nw,iaer)=0.
[253]	333
[1483]	334	end do ! L_NSPECTV
	335
	336	do nw=1,L_NSPECTI
	337	! Longwave
[1308]	338	do l=1,nlayer
[253]	339
[1308]	340	temp1=QIRsQREF3d(ig,nlayer+1-l,nw,iaer) &
	341	*QREFir3d(ig,nlayer+1-l,iaer)
[253]	342
[1308]	343	temp2=QIRsQREF3d(ig,max(nlayer-l,1),nw,iaer) &
	344	*QREFir3d(ig,max(nlayer-l,1),iaer)
[253]	345
	346	qxiaer(2*l,nw,iaer) = temp1
	347	qxiaer(2*l+1,nw,iaer)=(temp1+temp2)/2
	348
[1308]	349	temp1=temp1*omegair3d(ig,nlayer+1-l,nw,iaer)
	350	temp2=temp2*omegair3d(ig,max(nlayer-l,1),nw,iaer)
[253]	351
	352	qsiaer(2*l,nw,iaer) = temp1
	353	qsiaer(2*l+1,nw,iaer)=(temp1+temp2)/2
	354
[1308]	355	temp1=gir3d(ig,nlayer+1-l,nw,iaer)
	356	temp2=gir3d(ig,max(nlayer-l,1),nw,iaer)
[253]	357
	358	giaer(2*l,nw,iaer) = temp1
	359	giaer(2*l+1,nw,iaer)=(temp1+temp2)/2
	360
[1483]	361	end do ! nlayer
[253]	362
	363	qxiaer(1,nw,iaer)=qxiaer(2,nw,iaer)
[1308]	364	qxiaer(2*nlayer+1,nw,iaer)=0.
[253]	365
	366	qsiaer(1,nw,iaer)=qsiaer(2,nw,iaer)
[1308]	367	qsiaer(2*nlayer+1,nw,iaer)=0.
[253]	368
	369	giaer(1,nw,iaer)=giaer(2,nw,iaer)
[1308]	370	giaer(2*nlayer+1,nw,iaer)=0.
[253]	371
[1483]	372	end do ! L_NSPECTI
	373
	374	end do ! naerkind
[253]	375
[1483]	376	! Test / Correct for freaky s. s. albedo values.
[253]	377	do iaer=1,naerkind
	378	do k=1,L_LEVELS+1
	379
	380	do nw=1,L_NSPECTV
	381	if(qsvaer(k,nw,iaer).gt.1.05*qxvaer(k,nw,iaer))then
[726]	382	print*,'Serious problems with qsvaer values'
	383	print*,'in callcorrk'
[253]	384	call abort
	385	endif
	386	if(qsvaer(k,nw,iaer).gt.qxvaer(k,nw,iaer))then
	387	qsvaer(k,nw,iaer)=qxvaer(k,nw,iaer)
	388	endif
	389	end do
	390
	391	do nw=1,L_NSPECTI
	392	if(qsiaer(k,nw,iaer).gt.1.05*qxiaer(k,nw,iaer))then
[726]	393	print*,'Serious problems with qsiaer values'
	394	print*,'in callcorrk'
[253]	395	call abort
	396	endif
	397	if(qsiaer(k,nw,iaer).gt.qxiaer(k,nw,iaer))then
	398	qsiaer(k,nw,iaer)=qxiaer(k,nw,iaer)
	399	endif
	400	end do
	401
[1483]	402	end do ! L_LEVELS
	403	end do ! naerkind
[253]	404
	405	!-----------------------------------------------------------------------
	406	! Aerosol optical depths
[1483]	407	!-----------------------------------------------------------------------
[253]	408
	409	do iaer=1,naerkind ! a bug was here
	410	do k=0,nlayer-1
	411
	412	pweight=(pplay(ig,L_NLAYRAD-k)-pplev(ig,L_NLAYRAD-k+1))/ &
[1483]	413	(pplev(ig,L_NLAYRAD-k)-pplev(ig,L_NLAYRAD-k+1))
[253]	414	temp=aerosol(ig,L_NLAYRAD-k,iaer)/QREFvis3d(ig,L_NLAYRAD-k,iaer)
[588]	415	tauaero(2k+2,iaer)=max(temppweight,0.d0)
	416	tauaero(2k+3,iaer)=max(temp-tauaero(2k+2,iaer),0.d0)
[1483]	417
[253]	418	end do
	419	! boundary conditions
	420	tauaero(1,iaer) = tauaero(2,iaer)
	421	tauaero(L_LEVELS+1,iaer) = tauaero(L_LEVELS,iaer)
	422	!tauaero(1,iaer) = 0.
	423	!tauaero(L_LEVELS+1,iaer) = 0.
[1483]	424
	425	end do ! naerkind
[253]	426
[1483]	427	! Albedo and Emissivity.
	428	albi=1-emis(ig) ! Long Wave.
	429	DO nw=1,L_NSPECTV ! Short Wave loop.
[1482]	430	albv(nw)=albedo(ig,nw)
[1529]	431	ENDDO
[253]	432
[1483]	433	if (nosurf) then ! Case with no surface.
[1482]	434	DO nw=1,L_NSPECTV
	435	if(albv(nw).gt.0.0) then
[1529]	436	print*,'For open lower boundary in callcorrk must'
[1482]	437	print*,'have spectral surface band albedos all set to zero!'
	438	call abort
[1529]	439	endif
	440	ENDDO
[1482]	441	endif
[253]	442
[1483]	443	if ((ngrid.eq.1).and.(global1d)) then ! Fixed zenith angle 'szangle' in 1D simulations w/ globally-averaged sunlight.
[253]	444	acosz = cos(pi*szangle/180.0)
	445	print*,'acosz=',acosz,', szangle=',szangle
	446	else
[1483]	447	acosz=mu0(ig) ! Cosine of sun incident angle : 3D simulations or local 1D simulations using latitude.
[253]	448	endif
	449
	450	!-----------------------------------------------------------------------
	451	! Pressure and temperature
[1483]	452	!-----------------------------------------------------------------------
[253]	453
	454	DO l=1,nlayer
	455	plevrad(2*l) = pplay(ig,nlayer+1-l)/scalep
	456	plevrad(2*l+1) = pplev(ig,nlayer+1-l)/scalep
	457	tlevrad(2*l) = pt(ig,nlayer+1-l)
	458	tlevrad(2*l+1) = (pt(ig,nlayer+1-l)+pt(ig,max(nlayer-l,1)))/2
	459	END DO
	460
[600]	461	plevrad(1) = 0.
[1483]	462	plevrad(2) = 0. !! Trick to have correct calculations of fluxes in gflux(i/v).F, but the pmid levels are not impacted by this change.
[253]	463
	464	tlevrad(1) = tlevrad(2)
[1308]	465	tlevrad(2*nlayer+1)=tsurf(ig)
[253]	466
[1423]	467	pmid(1) = max(pgasmin,0.0001*plevrad(3))
	468	pmid(2) = pmid(1)
	469
[253]	470	tmid(1) = tlevrad(2)
[1423]	471	tmid(2) = tmid(1)
	472
	473	DO l=1,L_NLAYRAD-1
	474	tmid(2l+1) = tlevrad(2l+1)
	475	tmid(2l+2) = tlevrad(2l+1)
	476	pmid(2l+1) = plevrad(2l+1)
	477	pmid(2l+2) = plevrad(2l+1)
[253]	478	END DO
[1423]	479	pmid(L_LEVELS) = plevrad(L_LEVELS)
	480	tmid(L_LEVELS) = tlevrad(L_LEVELS)
[253]	481
[1423]	482	!!Alternative interpolation:
	483	! pmid(3) = pmid(1)
	484	! pmid(4) = pmid(1)
	485	! tmid(3) = tmid(1)
	486	! tmid(4) = tmid(1)
	487	! DO l=2,L_NLAYRAD-1
	488	! tmid(2l+1) = tlevrad(2l)
	489	! tmid(2l+2) = tlevrad(2l)
	490	! pmid(2l+1) = plevrad(2l)
	491	! pmid(2l+2) = plevrad(2l)
	492	! END DO
	493	! pmid(L_LEVELS) = plevrad(L_LEVELS-1)
	494	! tmid(L_LEVELS) = tlevrad(L_LEVELS-1)
	495
[1483]	496	! Test for out-of-bounds pressure.
[253]	497	if(plevrad(3).lt.pgasmin)then
	498	print*,'Minimum pressure is outside the radiative'
	499	print*,'transfer kmatrix bounds, exiting.'
	500	call abort
	501	elseif(plevrad(L_LEVELS).gt.pgasmax)then
	502	print*,'Maximum pressure is outside the radiative'
	503	print*,'transfer kmatrix bounds, exiting.'
	504	call abort
	505	endif
	506
[1483]	507	! Test for out-of-bounds temperature.
[253]	508	do k=1,L_LEVELS
	509	if(tlevrad(k).lt.tgasmin)then
	510	print*,'Minimum temperature is outside the radiative'
[1145]	511	print*,'transfer kmatrix bounds'
[858]	512	print*,"k=",k," tlevrad(k)=",tlevrad(k)
	513	print*,"tgasmin=",tgasmin
[1145]	514	if (strictboundcorrk) then
	515	call abort
	516	else
	517	print,'**********************************************'
	518	print*,'we allow model to continue with tlevrad=tgasmin'
	519	print*,' ... we assume we know what you are doing ... '
	520	print*,' ... but do not let this happen too often ... '
	521	print,'**********************************************'
[1147]	522	!tlevrad(k)=tgasmin
[1145]	523	endif
[253]	524	elseif(tlevrad(k).gt.tgasmax)then
[1648]	525	! print*,'Maximum temperature is outside the radiative'
	526	! print*,'transfer kmatrix bounds, exiting.'
	527	! print*,"k=",k," tlevrad(k)=",tlevrad(k)
	528	! print*,"tgasmax=",tgasmax
[1145]	529	if (strictboundcorrk) then
	530	call abort
	531	else
[1648]	532	! print,'**********************************************'
	533	! print*,'we allow model to continue with tlevrad=tgasmax'
	534	! print*,' ... we assume we know what you are doing ... '
	535	! print*,' ... but do not let this happen too often ... '
	536	! print,'**********************************************'
[1147]	537	!tlevrad(k)=tgasmax
[1145]	538	endif
[253]	539	endif
	540	enddo
[1016]	541	do k=1,L_NLAYRAD+1
	542	if(tmid(k).lt.tgasmin)then
	543	print*,'Minimum temperature is outside the radiative'
	544	print*,'transfer kmatrix bounds, exiting.'
[1145]	545	print*,"k=",k," tmid(k)=",tmid(k)
[1016]	546	print*,"tgasmin=",tgasmin
[1145]	547	if (strictboundcorrk) then
	548	call abort
	549	else
	550	print,'**********************************************'
	551	print*,'we allow model to continue with tmid=tgasmin'
	552	print*,' ... we assume we know what you are doing ... '
	553	print*,' ... but do not let this happen too often ... '
	554	print,'**********************************************'
	555	tmid(k)=tgasmin
	556	endif
[1016]	557	elseif(tmid(k).gt.tgasmax)then
[1648]	558	! print*,'Maximum temperature is outside the radiative'
	559	! print*,'transfer kmatrix bounds, exiting.'
	560	! print*,"k=",k," tmid(k)=",tmid(k)
	561	! print*,"tgasmax=",tgasmax
[1145]	562	if (strictboundcorrk) then
	563	call abort
	564	else
[1648]	565	! print,'**********************************************'
	566	! print*,'we allow model to continue with tmid=tgasmin'
	567	! print*,' ... we assume we know what you are doing ... '
	568	! print*,' ... but do not let this happen too often ... '
	569	! print,'**********************************************'
[1145]	570	tmid(k)=tgasmax
	571	endif
[1016]	572	endif
	573	enddo
[253]	574
	575	!=======================================================================
[1483]	576	! III. Calling the main radiative transfer subroutines
	577	!=======================================================================
[253]	578
	579
[1483]	580	Cmk= 0.01 * 1.0 / (glat(ig) * mugaz * 1.672621e-27) ! q_main=1.0 assumed.
[1529]	581	glat_ig=glat(ig)
[1194]	582
[253]	583	!-----------------------------------------------------------------------
[1483]	584	! Short Wave Part
	585	!-----------------------------------------------------------------------
[253]	586
[1483]	587	if(fract(ig) .ge. 1.0e-4) then ! Only during daylight.
[787]	588	if((ngrid.eq.1).and.(global1d))then
[253]	589	do nw=1,L_NSPECTV
[1483]	590	stel_fract(nw)= stel(nw)* 0.25 / acosz ! globally averaged = divide by 4, and we correct for solar zenith angle
[253]	591	end do
	592	else
	593	do nw=1,L_NSPECTV
[1161]	594	stel_fract(nw)= stel(nw) * fract(ig)
[253]	595	end do
[1483]	596	endif
	597
[253]	598	call optcv(dtauv,tauv,taucumv,plevrad, &
	599	qxvaer,qsvaer,gvaer,wbarv,cosbv,tauray,tauaero, &
[1648]	600	tmid,pmid,taugsurf,gweight)
[253]	601
	602	call sfluxv(dtauv,tauv,taucumv,albv,dwnv,wbarv,cosbv, &
[366]	603	acosz,stel_fract,gweight, &
[961]	604	nfluxtopv,fluxtopvdn,nfluxoutv_nu,nfluxgndv_nu, &
[253]	605	fmnetv,fluxupv,fluxdnv,fzerov,taugsurf)
	606
[1483]	607	else ! During the night, fluxes = 0.
[962]	608	nfluxtopv = 0.0d0
[1529]	609	fluxtopvdn = 0.0d0
[962]	610	nfluxoutv_nu(:) = 0.0d0
	611	nfluxgndv_nu(:) = 0.0d0
[253]	612	do l=1,L_NLAYRAD
[962]	613	fmnetv(l)=0.0d0
	614	fluxupv(l)=0.0d0
	615	fluxdnv(l)=0.0d0
[253]	616	end do
	617	end if
	618
[1482]	619
[1526]	620	! Equivalent Albedo Calculation (for OUTPUT). MT2015
	621	if(fract(ig) .ge. 1.0e-4) then ! equivalent albedo makes sense only during daylight.
	622	surface_stellar_flux=sum(nfluxgndv_nu(1:L_NSPECTV))
	623	if(surface_stellar_flux .gt. 1.0e-3) then ! equivalent albedo makes sense only if the stellar flux received by the surface is positive.
[1529]	624	DO nw=1,L_NSPECTV
	625	albedo_temp(nw)=albedo(ig,nw)*nfluxgndv_nu(nw)
[1526]	626	ENDDO
[1529]	627	albedo_temp(1:L_NSPECTV)=albedo_temp(1:L_NSPECTV)/surface_stellar_flux
[1526]	628	albedo_equivalent(ig)=sum(albedo_temp(1:L_NSPECTV))
	629	else
	630	albedo_equivalent(ig)=0.0 ! Spectrally Integrated Albedo not defined for non-irradiated grid points. So we arbitrary set the equivalent albedo to 0.
	631	endif
[1529]	632	else
	633	albedo_equivalent(ig)=0.0 ! Spectrally Integrated Albedo not defined for non-irradiated grid points. So we arbitrary set the equivalent albedo to 0.
	634	endif
[1482]	635
	636
[253]	637	!-----------------------------------------------------------------------
[1483]	638	! Long Wave Part
	639	!-----------------------------------------------------------------------
[253]	640
	641	call optci(plevrad,tlevrad,dtaui,taucumi, &
	642	qxiaer,qsiaer,giaer,cosbi,wbari,tauaero,tmid,pmid, &
[1648]	643	taugsurfi,gweight)
[538]	644
[253]	645	call sfluxi(plevrad,tlevrad,dtaui,taucumi,ubari,albi, &
	646	wnoi,dwni,cosbi,wbari,gweight,nfluxtopi,nfluxtopi_nu, &
	647	fmneti,fluxupi,fluxdni,fluxupi_nu,fzeroi,taugsurfi)
	648
	649	!-----------------------------------------------------------------------
	650	! Transformation of the correlated-k code outputs
	651	! (into dtlw, dtsw, fluxsurf_lw, fluxsurf_sw, fluxtop_lw, fluxtop_sw)
	652
	653	! Flux incident at the top of the atmosphere
[961]	654	fluxtop_dn(ig)=fluxtopvdn
[253]	655
	656	fluxtop_lw(ig) = real(nfluxtopi)
	657	fluxabs_sw(ig) = real(-nfluxtopv)
	658	fluxsurf_lw(ig) = real(fluxdni(L_NLAYRAD))
	659	fluxsurf_sw(ig) = real(fluxdnv(L_NLAYRAD))
[1482]	660
	661	! Flux absorbed by the surface. By MT2015.
	662	fluxsurfabs_sw(ig) = fluxsurf_sw(ig)*(1.-albedo_equivalent(ig))
[253]	663
	664	if(fluxtop_dn(ig).lt.0.0)then
	665	print*,'Achtung! fluxtop_dn has lost the plot!'
	666	print*,'fluxtop_dn=',fluxtop_dn(ig)
	667	print*,'acosz=',acosz
	668	print*,'aerosol=',aerosol(ig,:,:)
	669	print*,'temp= ',pt(ig,:)
	670	print*,'pplay= ',pplay(ig,:)
	671	call abort
	672	endif
	673
	674	! Spectral output, for exoplanet observational comparison
	675	if(specOLR)then
	676	do nw=1,L_NSPECTI
[526]	677	OLR_nu(ig,nw)=nfluxtopi_nu(nw)/DWNI(nw) !JL Normalize to the bandwidth
[253]	678	end do
	679	do nw=1,L_NSPECTV
[366]	680	!GSR_nu(ig,nw)=nfluxgndv_nu(nw)
[526]	681	OSR_nu(ig,nw)=nfluxoutv_nu(nw)/DWNV(nw) !JL Normalize to the bandwidth
[253]	682	end do
	683	endif
	684
	685	! Finally, the heating rates
	686
[586]	687	DO l=2,L_NLAYRAD
	688	dtsw(ig,L_NLAYRAD+1-l)=(fmnetv(l)-fmnetv(l-1)) &
[1194]	689	glat(ig)/(cppscalep(plevrad(2l+1)-plevrad(2*l-1)))
[586]	690	dtlw(ig,L_NLAYRAD+1-l)=(fmneti(l)-fmneti(l-1)) &
[1194]	691	glat(ig)/(cppscalep(plevrad(2l+1)-plevrad(2*l-1)))
[586]	692	END DO
[253]	693
	694	! These are values at top of atmosphere
[586]	695	dtsw(ig,L_NLAYRAD)=(fmnetv(1)-nfluxtopv) &
[1194]	696	glat(ig)/(cppscalep*(plevrad(3)-plevrad(1)))
[586]	697	dtlw(ig,L_NLAYRAD)=(fmneti(1)-nfluxtopi) &
[1194]	698	glat(ig)/(cppscalep*(plevrad(3)-plevrad(1)))
[253]	699
	700
[1483]	701	!-----------------------------------------------------------------------
	702	end do ! End of big loop over every GCM column.
	703	!-----------------------------------------------------------------------
[253]	704
[1483]	705
	706
[253]	707	!-----------------------------------------------------------------------
	708	! Additional diagnostics
[1483]	709	!-----------------------------------------------------------------------
[253]	710
[1483]	711	! IR spectral output, for exoplanet observational comparison
	712	if(lastcall.and.(ngrid.eq.1))then ! could disable the 1D output, they are in the diagfi and diagspec... JL12
[253]	713
[1483]	714	print*,'Saving scalar quantities in surf_vals.out...'
	715	print*,'psurf = ', pplev(1,1),' Pa'
	716	open(116,file='surf_vals.out')
	717	write(116,*) tsurf(1),pplev(1,1),fluxtop_dn(1), &
	718	real(-nfluxtopv),real(nfluxtopi)
	719	close(116)
[253]	720
[526]	721
[1483]	722	! USEFUL COMMENT - Do Not Remove.
	723	!
[526]	724	! if(specOLR)then
	725	! open(117,file='OLRnu.out')
	726	! do nw=1,L_NSPECTI
	727	! write(117,*) OLR_nu(1,nw)
	728	! enddo
	729	! close(117)
	730	!
	731	! open(127,file='OSRnu.out')
	732	! do nw=1,L_NSPECTV
	733	! write(127,*) OSR_nu(1,nw)
	734	! enddo
	735	! close(127)
	736	! endif
[253]	737
[1483]	738	! OLR vs altitude: do it as a .txt file.
	739	OLRz=.false.
	740	if(OLRz)then
	741	print*,'saving IR vertical flux for OLRz...'
	742	open(118,file='OLRz_plevs.out')
	743	open(119,file='OLRz.out')
	744	do l=1,L_NLAYRAD
	745	write(118,) plevrad(2l)
	746	do nw=1,L_NSPECTI
	747	write(119,*) fluxupi_nu(l,nw)
	748	enddo
	749	enddo
	750	close(118)
	751	close(119)
	752	endif
[253]	753
[305]	754	endif
[253]	755
[1483]	756	! See physiq.F for explanations about CLFvarying. This is temporary.
[1647]	757	if (lastcall) then
[470]	758	IF( ALLOCATED( gasi ) ) DEALLOCATE( gasi )
	759	IF( ALLOCATED( gasv ) ) DEALLOCATE( gasv )
[1315]	760	!$OMP BARRIER
	761	!$OMP MASTER
[470]	762	IF( ALLOCATED( pgasref ) ) DEALLOCATE( pgasref )
	763	IF( ALLOCATED( tgasref ) ) DEALLOCATE( tgasref )
	764	IF( ALLOCATED( pfgasref ) ) DEALLOCATE( pfgasref )
[1315]	765	!$OMP END MASTER
[1529]	766	!$OMP BARRIER
[861]	767	IF ( ALLOCATED(reffrad)) DEALLOCATE(reffrad)
	768	IF ( ALLOCATED(nueffrad)) DEALLOCATE(nueffrad)
[470]	769	endif
	770
[716]	771
[253]	772	end subroutine callcorrk

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

Download in other formats:

Original Format