ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SBDART version 1.4
c
c Created on Thu Jul 13 16:43:35 PDT 2000
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c sbdart
c
C============================================================================
c
C PURPOSE: COMPUTE PLANE-PARALLEL RADIATIVE TRANSFER IN THE EARTH'S ATMOSPHERE
C
C INPUT/OUTPUT VARIABLES described in rt.doc
C
c============================================================================
c
c
c The maximum number of radiation streams is set by the parameter NSTRMS
c while the maximum number of atmospheric layers is set by the parameter MXLYR
c To increase the number of streams or layers beyond the current limits, these
c parameters must be changed in each of the subroutines in which they appear
c (i.e., do a global replace, changing nstrms=40 to nstrms=??, or changing
c mxlyr=50 to mxlry=??). Remember to recompile after the changes are saved.
c
C SUBROUTINE CALLS
c
c absgas transfer gas absorption data
c absint compute path integrals for gaseous absorption
c atms load standard atmospheres (P,T,WH,WO)
c chkin check input parameters
c chkprn diagnostic print out
c cloudpar compute cloud rt properties (single scatter albedo...)
c depthscl compute radiative quantities for DISORT input
c disort solve multi-stream rt equation
c getmom compute phase functions
c kdistr compute k-distribution fits
c locate search on a 1-D table
c modatm modify atmosphere according to user inputs
c modmix modify constituent mix of atmosphere
c modrsc modify Rayleigh scattering strength for sensitivity studies
c openckf open diagnostic output files (when specified)
c rayleigh compute rayleigh scattering optical depth
c stdaer setup standard boundary layer aerosols
c saturate set relative humidity of cloud layer
c setfilt select sensor filter function
c stdout0 write output file header
c stdout1 write information at each frequency
c stdout2 write frequency integrated information
c suralb load surface albedo
c tauaero compute aerosol optical depth
c taugas compute absorption due to molecular species
c usraer setup user defined boundary layer aerosols.
c zensun compute solar zenith and azimuth given day and time
c zlayer find the computational layer with altitude z
c
c logical units files
c
c 11 INPUT
c
c 12 pmom direct access file
c (stays open throughout run)
c
c 13 usrcloud.dat, atms.dat, filter.dat, albedo.dat,
c aerosol.dat, aeroden.dat
c (each file is closed after read)
c
c 20-35 diagnostic output rtinfo.?? (?? = 00 --> 15)
c
C=======================================================================
program sbdart,57
parameter (lunc=20)
parameter (lun=11)
parameter (nstrms=40)
parameter (mxly=50)
parameter (maxulv=mxly+1)
parameter (naerz=5, ncldz=5)
parameter (re=6371.2)
parameter (nk=3)
real*8 zip
parameter (zip=-1.D0)
dimension z(mxly),p(mxly),t(mxly),wh(mxly),wo(mxly)
dimension taucld(mxly),gcld(mxly),wcld(mxly)
dimension dtaug(mxly),wreal(mxly)
dimension dtaur(mxly),taur(mxly)
dimension dtaua(mxly),waer(mxly),gaer(mxly),dtauab(mxly)
dimension dtcs(mxly),dtls(mxly)
dimension dtcv(mxly),dtlv(mxly)
dimension tauk(mxly,nk),dtauk(mxly,nk),wtk(mxly,nk),
& work(mxly,nk)
dimension umu(nstrms),utau(mxly)
parameter (nrstrm=nstrms)
dimension uzen(nrstrm),phi(nrstrm)
dimension sc(4)
parameter (mxq=63)
dimension uu(mxq,mxly)
parameter (nta=9)
parameter (naerw=47)
parameter (naerwm=naerw*nstrms)
dimension wlbaer(naerw),qbaer(naerw),wbaer(naerw),gbaer(naerw)
dimension pmaer(naerwm)
dimension airmass(mxly),gasabs(nta,mxly),trnsgas(nta,mxly)
dimension fdird(maxulv),flxdn(maxulv),flxup(maxulv),
& dfdt(maxulv),uavg(maxulv),uur(nrstrm,maxulv,nrstrm),
& u0u(nrstrm,maxulv),albmed(nstrms),trnmed(nstrms),hl(0:nstrms)
dimension zout(2)
dimension zcloud(ncldz),tcloud(ncldz),lcld(ncldz)
dimension zaer(naerz),taerst(naerz),jaer(naerz),laer(naerz)
dimension zbaer(mxly),dbaer(mxly)
dimension pmom(0:nstrms,mxly),temper(0:mxly)
parameter (npa=200) ! phase function arrays for high
dimension sphase(npa) ! angular resolution scattering
dimension sangle(npa) !
logical prnt(7),deltam,lamber,onlyfl,azmavg,usrtau,usrang,plank
real lwp(ncldz),nre(ncldz)
logical prnton
character*127 header
c=======================================================================
data pi/3.1415926536/,wvnmlo/0.0/,wvnmhi/0.0/,accur/0.0/
data nz, idatm,isat,nf,iaer,isalb,ickp,iout,imom
& / 33, 4, 0, 2, 0, 0, 0, 10, 0 /
data wlinf, wlsup, wlinc
& / .55, .55, 0./
data vis, rhaer, albcon, sza, csza, rhcld, amix
& / 23., zip, 0., 0., zip, zip, 0. /
data tbaer /zip/
data wlbaer /naerw*zip/
data qbaer /naerw*zip/
data wbaer /naerw*zip/
data gbaer /naerw*zip/
data dbaer /mxly*zip/
data zbaer /mxly*zip/
data abaer /zip/
data pmaer /naerwm*zip/
data sc/1.,0.,0.,0./
data iday,alat,alon,time /0,-64.767,-64.067,16./
data solfac /1.0/
data sclh2o, pbar, uw, uo3, o3trp, ztrp
& / zip, zip, zip, zip, zip, 0./
data xn2, xo2, xco2, xch4, xn2o
& / zip, zip, zip, zip, zip/
data xco, xno2, xso2, xnh3, xno, xhno3
& / zip, zip, zip, zip, zip, zip/
data xo4/1.0/
data nosct/0/
data xrsc/1.0/
data zpres/zip/
c
data zcloud, tcloud, nre, lcld
& /ncldz*0, ncldz*0., ncldz*8., ncldz*0/
data zaer, taerst, jaer, laer
& /naerz*0, naerz*0., naerz*0, naerz*0/
data prnt/ 7*.false./
data deltam, lamber, onlyfl, azmavg, usrang, usrtau
& / .true., .true., .true., .true., .false., .false./
data ntau,numu,nzen,nphi/4*0/
data umu/nstrms*0./
data uzen/nrstrm*zip/
data phi/nrstrm*zip/
data nothrm,ibcnd,phi0,fisot,temis/-1,0,0.,0.,0./
data btemp,ttemp /zip,zip/
data zout/0.,100./
data kdist/3/
data negprn/0/
data zgrid1,zgrid2,ngrid /0.,30.,mxly/
data krhclr /0/
data ipth /1/
data maxcmu/nstrms/
data maxumu,maxphi/2*nrstrm/
data nstr/4/
data lwp/ncldz*0./
data sphase/npa*zip/
data sangle/npa*zip/
c wavenumber resolution in CM-1
data wnres/20./
data ibeam/0/
c-----------------------------------------------------------------------
namelist /input/
& idatm, amix, isat, wlinf, wlsup, wlinc, sza, csza, solfac,
& nf, iday, time, alat, alon, zpres, pbar, sclh2o, uw, uo3,
& o3trp, ztrp, xrsc, xn2, xo2, xco2, xch4, xn2o, xco, xno2,
& xso2, xnh3, xno, xhno3, xo4, isalb, albcon, sc, zcloud,
& tcloud, lwp, nre, rhcld, krhclr, jaer, zaer, taerst, iaer,
& vis, rhaer, tbaer, wlbaer, qbaer, abaer, wbaer, gbaer, pmaer,
& zbaer, dbaer, nothrm, nosct, kdist, zgrid1, zgrid2, ngrid,
& ickp, zout, iout, deltam, lamber, ibcnd, phi0, prnt, ipth,
& fisot, temis, nstr, nzen, uzen, nphi, phi, imom, negprn,
& sphase, sangle, ttemp, btemp, ibeam
namelist /dinput/
& deltam,lamber,ibcnd,phi0,prnt,ipth,fisot,temis,nstr,
& nzen,uzen,nphi,phi,ttemp,btemp,imom
c-----------------------------------------------------------------------
c*****
iprint(kkk,nnn)=mod((kkk/2**nnn),2)
prnton(kkk,nnn)=iprint(kkk,nnn) .eq. 1
lunit(kkk,nnn)=iprint(kkk,nnn)*(lunc+nnn)
c*****
open(unit=11,file='INPUT',status='old',iostat=ios)
if(ios.eq.0) then
read(11,input,end=1)
read(11,dinput,end=2)
goto 2
1 continue
stop 'error: namelist block $INPUT not found'
2 continue
else
write(*,input)
stop
endif
zipp=zip
c assign user radiance angles
if(nphi.gt.1) then
if(numset
(zipp,phi,nrstrm).ne.2)
& write(*,'(a)') 'Error in MAIN -- ' //
& 'must specify exactly 2 values of phi when nphi is set'
if(nphi.gt.nrstrm)
& write(*,'(a)') 'Error in Main -- ' //
& 'specified nphi larger than nrstrm'
p1=min(phi(1),phi(2))
p2=max(phi(1),phi(2))
do i=1,nphi
phi(i)=p1+(i-1)*(p2-p1)/float(nphi-1)
enddo
else
nphi=numset(zipp,phi,nrstrm)
endif
if(nzen.ne.0) then
if(numset(zipp,uzen,nrstrm).ne.2)
& write(*,'(a)') 'Error in MAIN -- ' //
& 'must specify exactly 2 values of uzen when nzen is set'
if(nzen.gt.nrstrm)
& write(*,'(a)') 'Error in Main -- ' //
& 'specified nzen larger than nrstrm'
z1=min(uzen(1),uzen(2))
z2=max(uzen(1),uzen(2))
do i=1,nzen
uzen(i)=z1+(i-1)*(z2-z1)/float(nzen-1)
enddo
else
nzen=numset(zipp,uzen,nrstrm)
endif
c check input parameters
call chkin(isat,nf,wlinf,wlsup,solfac,isalb,albcon,
& sc,idatm,zpres,pbar,zcloud,tcloud,nre,jaer,zaer,
& taerst,iaer,zout,iout,prnt,lwp,nphi,nzen)
if(iout.le.0) stop
call openckf(ickp,lunc)
c solar spectrum, atmosphric profile, satellite sensor
c filtering, surface albedo and aerosol parameters
if(iday.ne.0) then
call zensun(abs(iday),time,alat,alon,sza,sazm,solfac)
if(iday.lt.0) then
write(*,'(i5,6f9.3)') abs(iday),time,alat,alon,sza,sazm,solfac
stop
endif
elseif(csza.ge.0.) then
sza=acos(csza)*180/pi
endif
if(abs(sza-90).lt..01) sza=95.
nosun=0
if(sza.ge.90.) nosun=1
call atms(idatm,amix,nz,z,p,t,wh,wo)
if(zgrid1.ne.0.) call zgrid(nz,z,p,t,wh,wo,zgrid1,zgrid2,ngrid)
if(zpres.gt.zip) then
call locate(z,nz,zpres,j)
fj=(zpres-z(j))/( pbar=p(j)*( endif
call modatm(nz,sclh2o,uw,uo3,o3trp,ztrp,pbar,z,p,wh,wo)
call modmix(xn2,xo2,xco2,xch4,xn2o,xco,xno2,xso2,xnh3,xno,xhno3)
c array dimensions
nstr1=nstr
mcldz=max(numset(0.,zcloud,ncldz),
& numset(0.,tcloud,ncldz),
& numset(0.,lwp,ncldz))
call zlayer(nz,z,mcldz,zcloud,lcld)
call zlayer(nz,z,numset(0.,taerst,naerz),zaer,laer)
lu1=lunit(ickp,1)
if(rhcld.ge.0) then
if(krhclr.eq.1) then
call satcloud(nz,ncldz,lcld,t,rhcld,wh,lu1)
else
call saturate(nz,ncldz,lcld,z,t,rhcld,wh,lu1)
endif
endif
call absint(uu,nz,z,p,t,wh,wo,lunit(ickp,0))
if(xrsc.ne.1.0) call modrsc(xrsc)
temper(0)=t(nz)
do 5 j=1,nz
temper(j)=t(nz+1-j)
5 continue
if(btemp.lt.0.) btemp=temper(nz)
if(ttemp.lt.0.) ttemp=temper(0)
c find nearest computational levels to given output altitudes
call nearest(z,nz,zout(1),nbot)
call nearest(z,nz,zout(2),ntop)
c write(*,'(4(10f10.3)/)') z
c write(*,*) 'ntop nbot ',ntop,nbot
ntop=nz-ntop+1
nbot=nz-nbot+2
c write(*,*) 'ntop nbot ',ntop,nbot
c set up for radiance calculation
if(nzen.ne.0) then
mpa=numset(zip,sangle,npa) ! number of angles at which aerosol
! phase function is specified
if(mpa.gt.0.and.tcloud(1).ne.0.) then
write(*,*) 'Error -- cloud optical depth cannot be ',
& ' treated in high precision radiance mode'
stop
endif
if(mpa.gt.0) then
if(numset(zip,sphase,npa).eq.0) then
g=gbaer(1)
c write(*,'(2a11)') 'sangle','sphase'
do i=1,mpa
am=cos(sangle(i)*pi/180)
sphase(i)=(1.-g**2)/(1.+g**2-2.*g*am)**1.5
c write(*,'(1p,2e11.3)') sangle(i),sphase(i)
enddo
endif
endif
usrang=.true.
onlyfl=.false.
numu=nzen
do 20 j=1,numu
umu(j)=min(1.,max(cos(uzen(numu+1-j)*pi/180.),-1.))
if(umu(j).eq.0.) then
if(j.eq.numu) then
umu(j)=-.0001
else
umu(j)=.0001
endif
endif
20 continue
if(nphi.ne.0) then
azmavg=.false.
else
nphi=1
phi(1)=0.
endif
endif
if(prnton(ickp,0))
& call chkprn(lunc,iday,time,alat,alon,sza,solfac,pbar,sclh2o,
& uw,uo3,o3trp,tcloud,zcloud,taerst,zaer,laer,
& zout,ntop,nbot,nz,z,p,t,wh,wo,uu)
call suralb(isalb,albcon,sc)
c load extinction, absorption, and asymmetry parameters for boundary
c layer aerosols, i.e., either rural, urban, oceanic or tropospheric aerosols
c and performs interpolation over relative humidity
if(rhaer.lt.0.) then
rh=relhum( else
rh=rhaer
endif
c set vertical profile
nzbaer=numset(zipp,zbaer,mxly)
ndbaer=numset(zipp,dbaer,mxly)
if(ndbaer.gt.0) then
if(nzbaer.eq.1) then
write(*,*) 'Error -- only one value of zbaer set'
stop
endif
if(nzbaer.ne.ndbaer.and.nzbaer.gt.1) then
write(*,*) 'Error -- number of elements must match:'
write(*,'(a,/,(1p10e11.3))') 'zbaer',(zbaer(i),i=1,nzbaer)
write(*,'(a,/,(1p10e11.3))') 'dbaer',(dbaer(i),i=1,ndbaer)
stop
endif
call uaerden(max(nzbaer,ndbaer),zbaer,dbaer)
endif
if(iaer.eq.5) then
call usraer(wlbaer,qbaer,wbaer,gbaer,pmaer,abaer,q55,imoma)
if(tbaer.eq.zip) tbaer=q55
else
call stdaer(iaer,rh)
endif
amu0=cos(sza*pi/180.0)
c load spectral response function
call setfilt(isat,wlinf,wlsup,wlinc,wl1,wl2,nwl)
call stdout0(iout,nwl,nz)
call zeroit(trnsgas,nta*mxly)
c beginning of wavelength looping
do 100 il=0,nwl-1
if(wlinc.gt.1) then
wl=(nwl-1)/(float(il)/wl2+float(nwl-1-il)/wl1)
dw1=wl-(nwl-1)/(float(il-1)/wl2+float(nwl-il)/wl1)
dw2=(nwl-1)/(float(il+1)/wl2+float(nwl-2-il)/wl1)-wl
elseif(wlinc.lt.0.) then
wr=wl2/wl1
wl=wl1*wr**(float(il)/(nwl-1))
dw1=wl-wl1*wr**(float(il-1)/(nwl-1))
dw2=wl1*wr**(float(il+1)/(nwl-1))-wl
else
wl=wl1+il*wlinc
dw1=wlinc
dw2=wlinc
endif
if(nothrm.lt.0) then
plank=wl.gt.2.
else
plank=nothrm.eq.0
endif
if(il.eq.0.and.il.ne.nwl-1) dw1=0.
if(il.eq.nwl-1.and.il.ne.0) dw2=0.
dw=.5*(dw1+dw2)
wvnm=(10000./wl)
dwvn=.0001*wvnm
wvnmlo=wvnm-dwvn
wvnmhi=wvnm+dwvn
c input flxin as solar energy (W/m2/um) or an arbitrary constant
ff=filter(wl)
if(nf.eq.0) then
if(amu0.ne.0.) flxin=pi/amu0
else
flxin=solfac*solirr(wl,nf)
endif
if(iout.eq.3) ff=ff*flxin
if(nosct.eq.1.or.nosun.eq.1) then
flxin=0.
amu0=1.
endif
rsfc=min(salbedo(wl),1.)
c calculate cloud optical depth and scattering parameters
call taucloud(nz,ncldz,wl,lcld,lwp,tcloud,
& nre,wcld,gcld,taucld,imom,pmom,lu1)
!write(*,'(a/1p(10e11.3))') 'wl',wl
!write(*,'(a/1p(10e11.3))') 'z',( !write(*,'(a/1p(10e11.3))') 'wcld',(wcld(jj),jj=1,nz)
!write(*,'(a/1p(10e11.3))') 'gcld',(gcld(jj),jj=1,nz)
!write(*,'(a/1p(10e11.3))') 'taucld',(taucld(jj),jj=1,nz)
!write(*,'(a,1p2e11.3)') 'wl,total ',wl,total(taucld,1,nz)
c boundary layer and stratospheric aerosols
call tauaero(wl,nz,z,vis,tbaer,naerz,taerst,laer,jaer,
& waer,gaer,dtaua,dtauab,imoma)
!write(*,'(a/1p(10e11.3))') 'wl,tbaer',wl,tbear
!write(*,'(a/1p(10e11.3))') 'z',( !write(*,'(a/1p(10e11.3))') 'dtaua',(dtaua(jj),jj=1,nz)
!write(*,'(a,1p2e11.3)') 'wl,tauaer ',wl,total(dtaua,1,nz)
c calculate optical thickness of gaseous absorption
c and rayleigh scattering (taur, dtaur)
call rayleigh(wl,uu,nz,taur,dtaur)
c TAUGAS returns -ln(molecular_transmission)
if(xo4.ne.1.0) call o4cfac(xo4)
call taugas(wl,uu,amu0,nz,z,dtcs,dtls,lunit(ickp,2))
call taugas(wl,uu,1.,nz,z,dtcv,dtlv,0)
call absgas(airmass,gasabs)
if(iout.eq.2) then
kdist=0
tauraytot=total(dtaur,1,nz)
tauaertot=total(dtaua,1,nz)
endif
c find number of passes, set up esft parameters for this wavelenth
if(kdist.eq.0.or.total(dtlv,1,nz).lt..01) then
kd=0
npass=1
else
kd=kdist
call kdistr(nz,mxly,mxq,uu,dtauk,tauk,wtk,work)
test=max(tauk(nz,1),tauk(nz,2),tauk(nz,3))
if(test.lt.0.01) then
npass=1
else
npass=nk
endif
c write(*,*) 'test = ',test,' npass = ',npass
endif
c loop through k-distribution terms, calculate total optical
c thickness (gases, rayleigh and clouds)
do 90 k=1,npass
c write(*,*) '0. maxulv=',maxulv
c skip call to DISORT if regular output not requested
if(.not. (iout.eq.2.or.iout.eq.3)) then
60 continue
call depthscl(kd,nz,k,npass,nstr,wl,amu0,z,dtls,
& dtcs,dtcv,dtlv,dtaur,dtaua,dtauab,taucld,dtauk,
& wtk,wcld,waer,gcld,gaer,dtaug,wreal,pmom,wt,
& imom,imoma,lunit(ickp,3))
if(npass.gt.1) ipth=0
c write(*,'(f10.3,3i4,1p7e12.4)')
c & wl,kd,k,npass,total(dtaur,1,nz),total(dtaua,1,nz),
c & total(taucld,1,nz),total(dtcv,1,nz),total(dtlv,1,nz),
c & total(dtaug,1,nz),dtaug(nz)
call disort(nz,dtaug,wreal,pmom,temper,wvnmlo,wvnmhi,
& usrtau,ntau,utau,nstr,usrang,numu,umu,nphi,phi,ibcnd,
& flxin,amu0,phi0,fisot,lamber,rsfc,hl,btemp,ttemp,
& temis,deltam,plank,onlyfl,accur,prnt,header,mxly,
& maxulv,maxumu,maxcmu,maxphi,fdird,flxdn,flxup,dfdt,
& uavg,uur,u0u,albmed,trnmed)
c add in atmospheric single scattering contribution when aerosol
c phase function (not moments) is specified.
c call zeroit(uur,nstrms*maxulv*nstrms)
if(ibeam.eq.2) call beamsct(nphi,phi,numu,umu,phi0,amu0,
& nz,dtaug,wreal,mpa,sangle,sphase,flxin,dtaur,iout,
& nbot,ntop,uur)
c stop
numu=nzen
c
c change nstr if a negative radiance is found within 5 layers
c of the output layer.
c
rmin=0.
if(iout.eq.20) then
call radrng(uur,nzen,nz+1,nphi,ntop,5,radmin,radmax)
elseif(iout.eq.21) then
call radrng(uur,nzen,nz+1,nphi,nbot,5,radmin,radmax)
elseif(iout.eq.22) then
call radrng(uur,nzen,nz+1,nphi,ntop,nz+1,radmin,radmax)
elseif(iout.eq.23) then
call radrng(uur,nzen,nz+1,nphi,ntop,5,radmin,radmax)
call radrng(uur,nzen,nz+1,nphi,nbot,5,radmn,radmx)
radmin=min(radmin,radmn)
radmax=max(radmax,radmx)
endif
if(negprn.eq.1.and.rmin.lt.0.) then
write(*,*) 'number of streams:',nstr
write(*,*) 'min,max radiance: ',radmin,radmax
do iz=1,nz+1,nz
write(*,'(x,a,i4)') 'iz = ',iz
do izen=1,nzen
write(*,'(f10.3,1p11e11.3)') uzen(izen),
& (uur(izen,iz,iphi),iphi=1,min(nphi,11))
enddo
enddo
endif
if(radmin.lt.-1.e-8*abs(radmax)) then
nstr=nstr-2
if(nstr.lt.max(4,nstr1-2)) nstr=min(nstrms,nstr1+2)
if(nstr.eq.nstr1) then
write(*,*)'%%% Could not repair negative radiance'
stop
endif
go to 60
else
nstr=nstr1
endif
c write(*,*) 'numu: ',numu
c write(*,*) 'umu: ',umu
c write(*,*) 'nphi: ',nphi
c write(*,*) 'phi: ',phi
c write(*,*) 'uur:'
c do ii=1,nzen
c write(*,'(1p10e11.3)') (uur(ii,nz+1,kk),kk=1,nphi)
c enddo
if(prnton(ickp,4)) then
write(lunc+4,'(a,i5,e11.3)') 'k, wt ',k,wt
write(lunc+4,'(3x,5a11)')
& 'dfdt','dflxup','fdird','flxnd','uavg'
write(lunc+4,'(i3,1p5e11.3)')
& (i,dfdt(i),flxup(i),fdird(i),flxdn(i),uavg(i),i=1,nz+1)
endif
endif
call stdout1(nz,z,ntop,nbot,iout,wl,dw,wt,
& fdird,flxdn,flxup,ff,nphi,nzen,sazm,phi,uzen,uur,
& gasabs,trnsgas,tauraytot,tauaertot,taucldtot,k,npass)
90 continue
100 continue
call stdout2(nz,iout,wl1,wl2,nphi,nzen,sazm,phi,uzen,z,
& airmass,trnsgas)
1000 format(/a/(1p10e11.3))
end