!-------------------------------------------------------------- ! ! MODULE OSS_IR_MODULE_SCAT: contains items needed for the ! OSS Forward Model (IR). AER Inc. 2004. ! !-------------------------------------------------------------- MODULE oss_ir_module_scat USE CloudModule USE asolar_source_function USE oss_addbl IMPLICIT NONE PRIVATE !------------------------------------------------------------------------ ! Public items made available to calling program(s) !------------------------------------------------------------------------ PUBLIC :: oss_init_ir,ossdrv_ir PUBLIC :: oss_destroy_ir !------------------------------------------------------------------------ ! Parameters/Dimensions !------------------------------------------------------------------------ INTEGER, PARAMETER :: mxlev=101,mxlay=mxlev-1,mxchan=2500,mxfsmp=9000 INTEGER, PARAMETER :: mxchsel=100,mxSf=50,mxHmol=10,MxmolS=7 INTEGER, PARAMETER :: mxtmp=10,mxsize1=mxlay*mxtmp,mxsize2=mxsize1*mxmols INTEGER, PARAMETER :: mxParG=(mxHmol+1)*mxlev+10 REAL, PARAMETER :: deg2rad = 0.017453293,pi=3.14159265 REAL, PARAMETER :: rair=1.020408163!= 10./g ->converts p(mb) in g/cm**2 !------------------------------------------------------------------------ ! Pressure Grid / Molecules / Geophysical Pointers !------------------------------------------------------------------------ INTEGER :: Nlev INTEGER*2 :: NmolFix,Nmol INTEGER*2 :: MolID(MxHmol),MolIDFix(MxHmol) INTEGER :: NparG,ITempG,ITskinG,IPsfcG,ImolInd(MxHmol),ICldIrG INTEGER :: IcldLiqG,IcldIceG REAL,DIMENSION(:),POINTER :: spectPtsGrd INTEGER :: nSpcPts COMMON/ossPgrd/NLev COMMON/MolCntrl/MolID,MolIDFix,NmolFix,Nmol COMMON/XPntrsG/NparG,ITempG,ITskinG,IPsfcG,ImolInd,ICldIrG !------------------------------------------------------------------------ ! Instrument parameters !------------------------------------------------------------------------ INTEGER :: nchan,nsf,nfsmp_ir REAL, dimension(:) , allocatable :: SfGrd INTEGER*2, dimension(:), allocatable :: nch_ir REAL*8, dimension(:) , allocatable :: vFreq REAL, dimension(:,:) , allocatable :: coef_ir REAL, dimension(:) , allocatable :: cFreq,pref,sunrad,nsmpPerChan INTEGER, dimension(:,:), allocatable :: ichmap_ir,Nsmp2NchanMap COMMON/Instr/nchan,nfsmp_ir COMMON/EmRfGrd/nsf ! Spectral grid on which surface properties ! are specified !------------------------------------------------------------------------ ! OD tables !------------------------------------------------------------------------ INTEGER :: NLayOD,NTmpOD,NWvpOD INTEGER*2, dimension(:) , allocatable :: NmolS INTEGER*2, dimension(:,:) , allocatable :: ImolS real, dimension (:,:) , allocatable :: Tmptab real, dimension (:,:) , allocatable :: Wvptab real, dimension (:,:) , allocatable :: kfix_ir,kh2o_ir,dkh2o_ir,kvar_ir !------------------------------------------------------------------------ ! Variables saved for interpolation of cloud properties !------------------------------------------------------------------------ INTEGER :: ip2!=1 REAL :: v2o REAL, DIMENSION(mxLay) :: atabs,btabs,atscat,btscat,aasym,basym REAL, DIMENSION(:,:), ALLOCATABLE :: tabs2d,tscat2d,asym2d CONTAINS !---------------------------------------------------------------------------- ! PURPOSE: Performs the 1st interface between the module and the ! calling program. ! It reads in the unit numbers, the files names, the geophysical vector ! pointers and returns back information that is useful to the user, ! pertaining to the frequency, the pressure grid, etc. ! !---------------------------------------------------------------------------- SUBROUTINE oss_init_ir(ius,iuo,sel_file,od_file,iuSol,sol_file,ctab_file, & nmol_in,MolID_in, & NparG_in,ITempG_in,ITskinG_in,IPsfcG_in,Imolind_in,IcldLiq_in,IcldIce_in, & nsf_in,SfGrd_in,dataType_cld,cldTyp,nHyd,nchan_out,cFreq_out, & nlev_out,Pref_out) !---Input variables CHARACTER(len=*) :: sel_file,od_file,sol_file,ctab_file CHARACTER(LEN=*),DIMENSION(:),INTENT(IN) :: dataType_cld CHARACTER(LEN=*),INTENT(IN) :: cldTyp INTEGER,INTENT(IN) :: nHyd INTEGER :: ius,iuo,NparG_in,iuSol INTEGER :: ITempG_in,ITskinG_in,IPsfcG_in,nsf_in INTEGER :: IcldLiq_in,IcldIce_in INTEGER :: nmol_in INTEGER, DIMENSION(:) :: molid_in INTEGER, DIMENSION(:) :: Imolind_in REAL, DIMENSION(:) :: SfGrd_in !---Output variables INTEGER :: nchan_out,nlev_out REAL, DIMENSION(:) :: cfreq_out,pref_out !---read solar file call getSolar(iuSol,sol_file) !---read OD and SEL data CALL GetOD(ius,iuo,sel_file,od_file,nmol_in,molid_in) !---Size conformity checks (safer interface with calling program/subroutine) IF (SIZE(cfreq_out) < nchan) STOP 'Insufficient Size for CFREQ_OUT' IF (SIZE(Pref_out) < nlev) STOP 'Insufficient Size for PREF_OUT' IF (SIZE(SfGrd_in) > MxSf) STOP 'SfGrd Vector too large' IF (SIZE(MolID_in) > mxHmol)STOP 'MolID Vector too large' IF (SIZE(Imolind_in)> mxHmol)STOP 'Imolind Vector too large' ITempG = ITempG_in ITskinG = ITskinG_in IPsfcG = IPsfcG_in Imolind(1:nmol) = Imolind_in(1:nmol) IcldLiqG = IcldLiq_in IcldIceG = IcldIce_in NparG = NparG_in Nsf = Nsf_in allocate (sunrad(nfsmp_ir),sfgrd(nsf)) SfGrd(1:nsf) = SfGrd_in(1:nsf) !---interpolate solar call InterpSolar(NFSmp_ir,vfreq,Sunrad) nchan_out = nchan cFreq_out(1:nchan) = cFreq(1:nchan) nlev_out = nlev Pref_out(1:nlev) = Pref(1:nlev) CALL initCloud(ctab_file,dataType_cld,spectPtsGrd,IcldLiqG,IcldIceG, & cldTyp,nHyd) nSpcPts=size(spectPtsGrd) allocate(tabs2d(mxLay,nSpcPts), tscat2d(mxLay,nSpcPts),& asym2d(mxLay,nSpcPts)) RETURN END SUBROUTINE oss_init_ir !---------------------------------------------------------------------------- ! PURPOSE: Reads the pre-computed Look Up Tables needed later in the ! OSS radiative transfer model. !---------------------------------------------------------------------------- SUBROUTINE GetOD(ius,iuo,selfile,odfile,nmol_in,molid_in) !---Input variables CHARACTER(len=*) :: selfile,odfile INTEGER :: ius,iuo,nmol_in INTEGER, DIMENSION(:) :: molid_in !---Local variables CHARACTER(len=100) :: filename,instr_info INTEGER*2 :: NmolS_tmp,ImolS_tmp(MxmolS),imol INTEGER :: imols_indx(MxmolS) INTEGER :: uid_sel,uid_od,k,nf_sel,nchmax,kk,ismp INTEGER :: isels,nsize1,nsize2,ks REAL :: odfac = 0. REAL , dimension (:) , allocatable :: kvar_tmp INTEGER, DIMENSION(MxHmol):: map,mapS,iflag = 1 !---open files filename=TRIM(selfile) OPEN(ius,file=filename,form='unformatted',status='old') filename=TRIM(odfile) OPEN(iuo,file=filename,form='unformatted',status='old') !---Check input molecular selection DO k=2,nmol_in IF(molid_in(k).LE.molid_in(k-1))STOP 'GetOD: HITRAN IDs must be in ascending order' END DO IF(molid_in(1).NE.1)STOP 'GetOD: Water vapor must be specified in molecular selection' !---Read sel_file header READ(ius)uid_sel READ(ius) READ(ius)instr_info READ(ius)nchan,nf_sel,nchmax nfsmp_ir=nf_sel !should read nfsmp_ir from OD file and compare allocate(cFreq(nchan)) READ(ius)cFreq(1:nchan) !---Sanity Checks (of Sel File) IF (nchan.GT.mxchan) STOP 'GetOD: Nchan > Mxchan ' IF (nf_sel.GT.mxfsmp) STOP 'GetOD: Nfsmp_ir > MxFSmp ' IF (nchmax.GT.mxchsel) STOP 'GetOD: NChmax > Mxchsel ' !---Get OD file header READ(iuo)uid_od READ(iuo) READ(iuo)nmolfix,nmol READ(iuo)molidfix(1:nmolfix),molid(1:nmol) READ(iuo)nlayod,ntmpod nsize1=nlayod*ntmpod !---Sanity Checks (of OD LUT) allocate (kfix_ir(nsize1,nfsmp_ir),dkh2o_ir(nsize1,nfsmp_ir)) allocate (kvar_ir(mxmols*nsize1,nfsmp_ir),kh2o_ir(nsize1,nfsmp_ir)) allocate (wvptab(nmol+2,nlayod)) allocate (tmptab(ntmpod,nlayod),pref(nlayod+1)) allocate (vfreq(nfsmp_ir)) allocate (coef_ir(nchmax,nfsmp_ir),ichmap_ir(nchmax,nfsmp_ir)) allocate (nch_ir(nfsmp_ir)) allocate (nmols(nfsmp_ir),imols(mxmols,nfsmp_ir)) IF (nlayod.GT.mxlay) STOP 'GetOD: NLay > Mxlay ' IF (ntmpod.GT.mxtmp) STOP 'GetOD: Ntmp > Mxtmp ' !IF (uid_sel.NE.uid_od) STOP 'Inconsistent IDs in Sel & OD LUTs' nlev=nlayod+1 READ(iuo)pref(1:nlev),tmptab(1:ntmpod,1:nlayod),wvptab(1:nmol+2,1:nlayod) !---Build mapping vector for user-selected variable molecules kk=0 map(1:nmol)=0 DO k=1,nmol !---map(k) contains new relative index for selected molecule and is zero for others IF(ANY (molid_in(1:nmol_in) == molid(k))) THEN kk=kk+1 map(k)=kk END IF END DO IF(kk.NE.nmol_in)THEN PRINT *,'This database contains only the molecules with' PRINT *,'the following HITRAN ID as variable molecules: ' PRINT *,molid(1:nmol) STOP 'GetOD: ' END IF !---Read oss parameters DO ismp=1,nfsmp_ir READ(ius)iselS,nch_ir(ismp) ! iselS not needed in present application READ(ius)coef_ir(1:nch_ir(ismp),ismp),ichmap_ir(1:nch_ir(ismp),ismp) END DO !---Read absorption coefficient tables PRINT *, 'GetOD: Loading absorption tables for '//TRIM(instr_info) DO ISmp=1,Nfsmp_ir READ(iuo)vfreq(ismp),nmols_tmp READ(iuo)Imols_tmp(1:nmols_tmp) nsize2=nsize1*(nmols_tmp-1) allocate (kvar_tmp(max(1,nsize2))) READ(iuo)kfix_ir(1:nsize1,ismp),dkh2o_ir(1:nsize1,ismp),& kh2o_ir(1:nsize1,ismp),kvar_tmp(1:nsize2) !---code uses kfix_ir*rfix instead of kfix_ir CALL init_kfix(kfix_ir(1,ismp),wvptab) ! Keep only user-selected molecules as variables molecules - ! merge remaining molecules with fixed gases ! First apply OD threshold to identify weakly absorbing constituents iflag(2:nmol) = 0 CALL odthresh(kfix_ir(1,ismp),kvar_tmp,imols_tmp,nmols_tmp,wvptab,odfac,iflag) maps(1:nmol) = map(1:nmol) * iflag(1:nmol) kk = 0 DO ks = 1,nmols_tmp imol = imols_tmp(ks) IF(mapS(imol) > 0)THEN kk = kk+1 imols_indx(kk) = ks ELSE CALL cum_fix(kfix_ir(1,ismp),kvar_tmp,imol,nmols_tmp,ks,wvptab) END IF END DO nmols(ismp) = kk CALL shrink_var(kvar_tmp,imols_tmp,nmols_tmp,mapS,imols_indx,& kvar_ir(1,ISmp),imols(1,ISmp),nmols(ISmp)) deallocate(kvar_tmp) END DO nmol = nmol_in molid(1:nmol) = molid_in(1:nmol) CLOSE(ius) CLOSE(iuo) RETURN END SUBROUTINE GetOD SUBROUTINE odthresh(kfix0,kvar0,imols,nmols,w,odfac,iflag) !---Input variables INTEGER*2 :: NmolS,ImolS(NmolS) REAL, DIMENSION(NlayOD,NtmpOD) :: kfix0 REAL, DIMENSION(NmolS-1,NlayOD,NtmpOD) :: kvar0 REAL :: odfac REAL, DIMENSION(:,:) :: w !---Output variables INTEGER, DIMENSION(:) :: iflag !---Local variables REAL :: kbuf(50),ktot INTEGER :: nt,l,ks ! Compare optical depth of individual dry constituents to ! total optical depth and flag weakly absorbing molecules DO nt = 1,NtmpOD DO l = 1,NlayOD ktot = kfix0(l,nt) DO ks = 2,nmols kbuf(ks-1) = kvar0(ks-1,l,nt) * w(imols(ks)+2,l) ktot = ktot + kbuf(ks-1) END DO DO ks=2,nmols ! could exit at first occurence of iflag = 1 IF(kbuf(ks-1) > odfac * ktot) iflag(imols(ks)) = 1 END DO END DO END DO RETURN END SUBROUTINE odthresh SUBROUTINE cum_fix(kfix0,kvar0,imol,nmols,ks,w) !---In/Out variables INTEGER :: ks INTEGER*2 :: NmolS,Imol REAL :: kfix0(NlayOD,NtmpOD),kvar0(NmolS-1,NlayOD,NtmpOD) REAL, DIMENSION(:,:) :: w !---Local variables INTEGER :: nt DO nt=1,ntmpOD kfix0(1:NLayOD,nt) = kfix0(1:NLayOD,nt) +& kvar0(ks-1,1:NLayOD,nt) * w(imol+2,1:nlayod) END DO RETURN ENTRY init_kfix(kfix0,w) DO nt=1,ntmpOD kfix0(1:NLayOD,nt) = kfix0(1:NLayOD,nt) * w(1,1:NlayOD) END DO RETURN END SUBROUTINE cum_fix SUBROUTINE shrink_var(kvar0,imols0,nmols0,maps,imols_indx,kvar1,imols1,nmols1) !---In/Out variables INTEGER*2 :: nmols0,imols0(nmols0) INTEGER*2 :: nmols1,imols1(nmols1) INTEGER, DIMENSION(:) :: maps,imols_indx REAL :: kvar0(nmols0-1,NlayOD,NtmpOD) REAL :: kvar1(nmols1-1,NlayOD,NtmpOD) kvar1(1:nmols1-1,1:NLayOD,1:NtmpOD) = & kvar0(imols_indx(2:nmols1)-1,1:NLayOD,1:NtmpOD) imols1(1:nmols1) = maps(imols0(imols_indx(1:nmols1))) RETURN END SUBROUTINE shrink_var !---------------------------------------------------------------------------- ! PURPOSE: Driver for the IR radiative transfer model. ! Computes both radiances and their jacobians wrt geophysical ! parameters. ! xG: Profile vector of geophysical parmaters (containing ! temperature and constiutuents profiles, surface ! pressure and skin temperature, and cloud parameters ! EmRf: vector of MW surface emissivities. ! Path: Array of quantities characterizing the atmospehric ! path in plane-parallel atmosphere. ! Y: Vector of radiances. ! xkt: Array of derivatives of the radiances wrt geophysical parameters. ! xkEmRf: Array of radiance derivatives wrt surface emissivities. !---------------------------------------------------------------------------- SUBROUTINE ossdrv_ir(xG,EmRf,obsang,sunang,y,xkt,xkEmRf,pland) !---Input variables REAL, DIMENSION(:) :: xG REAL, DIMENSION(2,Nsf) :: EmRf REAL :: pland,obsang,sunang !INTEGER, DIMENSION(:) :: kchan !---Output variables REAL, DIMENSION(:) :: y REAL, DIMENSION(NparG,nchan) :: xkt REAL, DIMENSION(2,Nsf,nchan) :: XkEmRf !---Local variables LOGICAL :: lookup=.false.,sun REAL, DIMENSION(MxLay) :: tavl,wfix,dtu,dtl REAL, DIMENSION(MxHmol,MxLay) :: q,w REAL, DIMENSION(MxLay,MxHmol) :: dwqu,dwql INTEGER, DIMENSION(MxLay) :: indx REAL, DIMENSION(MxLay) :: a1,a2,ad1,ad2,tautot,dtaudtmp,amt REAL, DIMENSION(MxmolS,MxLay) :: abso REAL, DIMENSION(mxLay) :: tabs,tscat,asym REAL, DIMENSION(mxLay) :: flxnetMono REAL :: vEmRf(2),xkt_tmp(MxParG),xkEmRf_tmp(2) REAL :: umu,umu0,delphi,vn,a,fbeam,rad,xx,xkEm INTEGER :: nsurf,nobs,n1,n2,ip0,nn,ich,i,ich0 INTEGER,PARAMETER :: iflux=0,nStr=4 !====================================================================== ! Initialize radiance vector and k-matrix !====================================================================== y(1:Nchan) = 0. xkt(1:NparG,1:Nchan) = 0. xkEmRf(1:2,1:Nsf,1:Nchan) = 0. !====================================================================== ! Compute path variables !====================================================================== ! Compute path geometry CALL setpath_ir(xG,obsang,sunang,nsurf,nobs,lookup,n1,n2,sun,umu,umu0,delphi) !---Compute average temperature and molecular amounts for the layers CALL fpath_ir(xG,N2,tavl,wfix,q,w,dtu,dtl,dwqu,dwql) !---Compute coefficients for temperature interpolation of ODs CALL settabindx(tavl,N2,indx,a1,a2,ad1,ad2) !====================================================================== ! Loop over spectral points !====================================================================== ip0 = 2 vn = real(vFreq(1))-0.01 CALL planck_set(vn,tavl,xG(ITskinG),n2) !---Get profiles of optical properties for input to scattering module: CALL setCloudOptLayer(xG,tavl,nSurf,pref,xG(iPsfcG),iCldLiqG,iCldIceG, & tabs2d,tscat2d,asym2d) ip2=2 v2o=0. NfsmpLoop: DO nn=1,Nfsmp_ir !if(ANY(kchan(ichMap(1:nch(nn),nn)) == 1)) then !---Determine optical properties for nnth spectral point: CALL optInt(vFreq(nn),nSurf,tabs,tscat,asym) !---Compute molecular optical depth for all atmospheric layers CALL OSStran(kfix_ir(1,nn),kh2o_ir(1,nn),dkh2o_ir(1,nn),kvar_ir(1,nn),indx,a1,a2,& ad1,ad2,N2,NmolS(nn),ImolS(1,nn),wfix,q,w,tautot,abso,dtaudtmp) !---Interpolate input surface emissivity to node wavenumber CALL vinterp(EmRf,2,nn,vEmrf,ip0,a) vn=vFreq(nn) fbeam=sunrad(nn) IF (.NOT.sun)fbeam=0. !---Perform RT calculations !clear sky model !CALL OSSrad(tautot,abso,dtaudtmp,tavl,dtu,dtl,dwqu,dwql,nmols(nn),& ! imols(1,nn),xG,vEmRf,fbeam,vn,N1,N2,sun,umu,umu0,lookup,rad, & ! xkt_tmp,xkEmRf_tmp) CALL scattRT(tautot,xG,xG(iTskinG),tabs,tscat,asym,vEmRf(1),fBeam,vn,nSurf-1,& nobs,sun,nStr,umu,umu0,lookup,delphi,iflux,flxNetMono, & rad,xkt_tmp,xkEm) nchLoop: DO ich=1,nch_ir(nn) ich0 = ichMap_ir(ich,nn) y(ich0) = y(ich0)+rad*coef_ir(ich,nn) xkt(1:NParG,ich0) = xkt(1:NParG,ich0)+xkt_tmp(1:NParG)*coef_ir(ich,nn) DO i=1,2 xx = xkEmRf_tmp(i)*coef_ir(ich,nn) xkEmRf(i,ip0-1,ich0) = xkEmRf(i,ip0-1,ich0)+ xx*(1-a) xkEmRf(i,ip0,ich0) = xkEmRf(i,ip0,ich0) + xx*a END DO END DO nchLoop !EINDIF END DO NfsmpLoop RETURN END SUBROUTINE ossdrv_ir !---------------------------------------------------------------------------- ! PURPOSE: Interpolation in wavenumber !---------------------------------------------------------------------------- SUBROUTINE vinterp(datin,nx,nn,datout,ip0,a) !---Input variables INTEGER :: nx,nn,ip0 REAL :: datin(nx,nfsmp_ir) !---Output variables REAL :: a,datout(nx) !---Local variables REAL :: vn INTEGER :: ip,i vn = vFreq(nn) IpLoop: DO ip=ip0,nsf-1 IF(SfGrd(ip).GT.vn) EXIT IpLoop END DO IpLoop ip0=ip a=(vn-SfGrd(ip0-1))/(SfGrd(ip0)-SfGrd(ip0-1)) DO i=1,nx datout(i)=a*(datin(i,ip0)-datin(i,ip0-1))+datin(i,ip0-1) END DO RETURN END SUBROUTINE vinterp !---------------------------------------------------------------------------- ! PURPOSE: Computes items related to the viewing geometry. !---------------------------------------------------------------------------- SUBROUTINE setpath_ir(xG,obsang,sunang,nsurf,nobs,lookup,n1,n2,sun,umu,umu0,delphi) !---Input variables REAL, DIMENSION(:) :: xG REAL :: obsang,sunang LOGICAL :: lookup !---Output variables INTEGER :: Nobs,Nsurf,N1,N2 LOGICAL :: sun REAL :: umu,umu0,delphi !---Local variables REAL :: Psurf,viewang INTEGER :: i !---Compute the surface level nobs=1 Psurf=xG(IPsfcG) NlevLoop:DO i=2,nlev-1 IF(pref(i).GE.psurf) EXIT NlevLoop ENDDO NlevLoop nsurf=i !---Set viewing geometry parameters delphi = 0. viewang = obsang IF(lookup)THEN N1=NSurf N2=Nobs-1 ELSE IF(viewang > 90.)STOP 'EIA must be less than 90' N1=Nobs N2=NSurf-1 END IF umu =COS(viewang*deg2rad) IF(sunang.LT.80.)THEN sun=.TRUE. umu0=COS(sunang*deg2rad) ELSE sun=.FALSE. umu0=1. END IF RETURN END SUBROUTINE setpath_ir !---------------------------------------------------------------------------- ! PURPOSE: This subroutine calculates the average temperature and molecular ! amounts for all layers for given profiles of temperature and ! molecular concentrations. It also calculates the derivatives of ! tavl with respect to a change in the lower and upper boundary ! temperatures and the derivatives of the molecular amounts with ! respect to a change in the mixing ratios at the layer ! boundaries. Molecular amounts are in molec./cm**2. ! Integration assumes that T is linear in z (LnT linear in LnP) ! and LnQ linear in LnP. !---------------------------------------------------------------------------- SUBROUTINE fpath_ir(xG,Nlast,tavl,wfix,q,w,dtu,dtl,dwqu,dwql) !---Input variables INTEGER :: nlast REAL, DIMENSION(:) :: xG !---Output variables REAL, DIMENSION(MxLay) :: tavl,wfix,dtu,dtl REAL, DIMENSION(MxHmol,MxLay) :: q,w REAL, DIMENSION(MxLay,MxHmol) :: dwqu,dwql !---Local variables INTEGER :: nsurf,l,IH2o,k,IXoff,n REAL :: psfc,tsfc,scal,scal2,qsfc,wtot,wdry REAL, DIMENSION(mxlev) :: dp,qcor,qr NSurf=Nlast+1 Psfc=xG(IPsfcG) !------------------------------------------------------------------------ ! Compute Average Temperature for the layers !------------------------------------------------------------------------ DO l=1,Nlast-1 dp(l)=(pref(l+1)-pref(l)) scal=1./dp(l) CALL lpsum_log(pref(l),pref(l+1),xG(ITempG+l-1),xG(ITempG+l),& scal,tavl(l),dtu(l),dtl(l)) END DO !---Surface layer dp(Nlast)=(PSfc-pref(NSurf-1)) scal=1./dp(Nlast) scal2=alog(PSfc/pref(NSurf))/alog(pref(NSurf-1)/pref(NSurf)) CALL lint_log(xG,pref,Nsurf,Psfc,Tsfc) CALL lpsum_log(pref(NSurf-1),PSfc,xG(NSurf-1),TSfc,& scal,tavl(Nlast),dtu(Nlast),dtl(Nlast)) dtu(Nlast)=dtu(Nlast)+dtl(Nlast)*scal2*tsfc/xG(NSurf-1) dtl(Nlast)=dtl(Nlast)*(1-scal2)*tsfc/xG(NSurf) !------------------------------------------------------------------------ ! Calculate amounts for individual species and derivatives wrt mixing ! ratios for retrieved constituents. !------------------------------------------------------------------------ scal=rair IH2o=ImolInd(1)-1 qcor(1:Nsurf)=1/(1.+xG(IH2o+1:IH2o+Nsurf)) DO k=1,nmol IXoff=ImolInd(k)-1 !---Transform mix. ratios into mass fractions (relative to total air mass) qr(1:Nsurf)=xG(IXoff+1:IXoff+NSurf)*qcor(1:Nsurf) DO l=1,Nlast-1 CALL lpsum_log(pref(l),pref(l+1),qr(l),qr(l+1),& scal,w(k,l),dwqu(l,k),dwql(l,k)) END DO !---Surface Layer CALL lint_log(qr,pref,nsurf,psfc,qsfc) CALL lpsum_log(pref(NSurf-1),psfc,qr(NSurf-1),qsfc,& scal,w(k,Nlast),dwqu(Nlast,k),dwql(Nlast,k)) dwqu(Nlast,k)=dwqu(Nlast,k) +dwql(Nlast,k)*scal2*qsfc/qr(NSurf-1) dwql(Nlast,k)=dwql(Nlast,k)*(1-scal2)*qsfc/qr(NSurf) END DO !---Derivatives of amount wrt dry mixing ratios DO n=1,Nlast dwqu(n,1)=dwqu(n,1)*qcor(n)**2 dwql(n,1)=dwql(n,1)*qcor(n+1)**2 dwqu(n,2:nmol)=dwqu(n,2:nmol)*qcor(n) dwql(n,2:nmol)=dwql(n,2:nmol)*qcor(n+1) END DO !---Compute fix gas amount and average layer mixing ratios DO l=1,Nlast wtot=dp(l)*rair wdry=wtot-w(1,l) wfix(l)=wdry q(1,l)=w(1,l)/wtot q(2:nmol,l)=w(2:nmol,l)/wdry END DO RETURN END SUBROUTINE fpath_ir !---------------------------------------------------------------------------- ! PURPOSE: Computes average layer quantities (or integrated amount) ! using a linear dependence on P. !---------------------------------------------------------------------------- SUBROUTINE lpsum(pu,pl,xu,xl,scal,xint,dxu,dxl) !---Input variables REAL, INTENT(IN) :: pu,pl,xu,xl,scal !---Output variables REAL, INTENT(OUT) :: xint,dxu,dxl xint = 0.5*(pl-pu)*scal dxu = xint dxl = xint xint = xint*(xu+xl) RETURN END SUBROUTINE lpsum !---------------------------------------------------------------------------- ! PURPOSE: Computes average layer quantities (or integrated amount) ! using a log-x dependence on log-p. !---------------------------------------------------------------------------- SUBROUTINE lpsum_log(pu,pl,xu,xl,scal,xint,dxu,dxl) REAL, PARAMETER :: epsiln=1.e-05 !---Input variables REAL, INTENT(IN) :: pu,pl,xu,xl,scal !---Output variables REAL, INTENT(OUT) :: xint,dxu,dxl !---Local variables REAL :: hp,x0,zeta,alza,alpha hp = alog(pl/pu) x0 = pl*xl*hp zeta = pu*xu/(pl*xl) IF(ABS(zeta-1.).GT.epsiln)THEN alza = alog(zeta) xint = x0*(zeta-1.)/alza alpha = zeta/(zeta-1.)-1./alza ELSE xint = x0*2./(3.-zeta) alpha = zeta/(3.-zeta) END IF xint = xint*scal dxu = xint*alpha/xu dxl = xint*(1.-alpha)/xl RETURN END SUBROUTINE lpsum_log !---------------------------------------------------------------------------- ! PURPOSE: Interpolation in linear pressure !---------------------------------------------------------------------------- SUBROUTINE lint(x,p,N0,p0,x0) INTEGER :: n0 REAL :: p0 REAL, DIMENSION(:), INTENT(IN) :: x,p !---Output variable REAL :: x0 !---Local variable REAL :: xx xx = (x(N0)-x(N0-1))/(p(N0)-p(N0-1)) x0 = x(N0-1)+(p0-p(N0-1))*xx RETURN END SUBROUTINE lint !---------------------------------------------------------------------------- ! PURPOSE: Interpolation in Log pressure !---------------------------------------------------------------------------- SUBROUTINE lint_log(x,p,N0,p0,x0) INTEGER :: n0 REAL :: p0 REAL, DIMENSION(:), INTENT(IN) :: x,p !---Output variable REAL :: x0 !---Local variable REAL :: xx xx = alog(x(N0)/x(N0-1))/alog(p(N0)/p(N0-1)) x0 = x(N0-1)*(p0/p(N0-1))**xx RETURN END SUBROUTINE lint_log !---------------------------------------------------------------------------- ! PURPOSE: Given the LUTs, the function computes the layer ! optical depths. Refer to ossrad_mw for the transmittance and ! brightness temperature calculation !---------------------------------------------------------------------------- SUBROUTINE OSStran(kfix,kh2o,dkh2o,kvar,indx,a1,a2,ad1,ad2,N2,& NmolS,ImolS,wfix,q,w,tautot,abso,dtaudtmp) !---Input variables INTEGER*2 :: NmolS,ImolS(NmolS) INTEGER :: n2 INTEGER, DIMENSION(:) :: indx REAL, DIMENSION(:) :: a1,a2,ad1,ad2,wfix REAL, DIMENSION(:,:) :: q,w REAL, DIMENSION(NlayOD,NtmpOD) :: kfix,kh2o,dkh2o REAL :: kvar(NmolS-1,NlayOD,NtmpOD) !---Output variables REAL, DIMENSION(:) :: tautot,dtaudtmp REAL, DIMENSION(:,:) :: abso !---Local variables INTEGER :: l,indx1,indx2,indx3,ks,imol REAL :: al,bl,adl,bdl,df1,df2,abs0,dabs0 REAL :: absh2o1,absh2o2,absh2o3,absh2o REAL :: dabsh2odq,dabsh2o DO l=1,N2 indx1 = indx(l)-1 indx2 = indx1+1 indx3 = indx1+2 al=a1(l) bl=a2(l) adl=ad1(l) bdl=ad2(l) !---Fixed gases df1 = kfix(l,indx1) - kfix(l,indx3) df2 = kfix(l,indx2) - kfix(l,indx3) abs0 = al*df1 + bl*df2 + kfix(l,indx3) tautot(l) = abs0 * wfix(l) dabs0 = adl*df1 + bdl*df2 dtaudtmp(l) = dabs0 * wfix(l) abso(1,l) = -abs0 !---Water vapor absh2o1 = q(1,l) * dkh2o(l,indx1) + kh2o(l,indx1) absh2o2 = q(1,l) * dkh2o(l,indx2) + kh2o(l,indx2) absh2o3 = q(1,l) * dkh2o(l,indx3) + kh2o(l,indx3) df1 = absh2o1 - absh2o3 df2 = absh2o2 - absh2o3 absh2o = al*df1 + bl*df2 + absh2o3 tautot(l) = tautot(l) + absh2o * w(1,l) dabsh2odq = al * (dkh2o(l,indx1)-dkh2o(l,indx3)) & + bl * (dkh2o(l,indx2)-dkh2o(l,indx3)) + dkh2o(l,indx3) abso(1,l) = dabsh2odq * q(1,l) + absh2o + abso(1,l) dabsh2o = adl*df1 + bdl*df2 dtaudtmp(l) = dtaudtmp(l) + dabsh2o * w(1,l) !---Variable gases DO kS=2,nmolS Imol=ImolS(kS) df1 = kvar(kS-1,l,indx1) - kvar(kS-1,l,indx3) df2 = kvar(kS-1,l,indx2) - kvar(kS-1,l,indx3) abs0 = al*df1 + bl*df2 + kvar(kS-1,l,indx3) dabs0 = adl*df1 + bdl*df2 abso(kS,l) = abs0 tautot(l) = tautot(l) + abs0 * w(Imol,l) dtaudtmp(l) = dtaudtmp(l) + dabs0 * w(Imol,l) abso(1,l) = abso(1,l) - abs0 * q(Imol,l) END DO END DO RETURN END SUBROUTINE OSStran !---------------------------------------------------------------------------- ! PURPOSE: Computes temperature/Water vapor indexes and interpolation ! coefficients. !---------------------------------------------------------------------------- SUBROUTINE settabindx(tavl,N2,mindx,a1,a2,ad1,ad2) !---Input variables INTEGER :: n2 REAL, DIMENSION(:) :: tavl !---Output variables INTEGER, DIMENSION(:) :: mindx REAL, DIMENSION(:) :: a1,a2,ad1,ad2 !---Local variables INTEGER :: l,i,j REAL :: den1,den2 !---Compute coefficients for temperature interpolation of ODs DO l=1,N2 TempLoop: DO i=3,ntmpod-2 IF(tmptab(i,l).GE.tavl(l)) exit TempLoop END DO TempLoop IF (ABS(tavl(l)-tmptab(i-1,l)).LE.ABS(tavl(l)-tmptab(i,l))) THEN j=i-1 ELSE j=i END IF mindx(l)=j den1=(tmptab(j-1,l)-tmptab(j,l))*(tmptab(j-1,l)-tmptab(j+1,l)) den2=(tmptab(j,l)-tmptab(j-1,l))*(tmptab(j,l)-tmptab(j+1,l)) a1(l)=(tavl(l)-tmptab(j,l))*(tavl(l)-tmptab(j+1,l))/den1 a2(l)=(tavl(l)-tmptab(j-1,l))*(tavl(l)-tmptab(j+1,l))/den2 ad1(l)=(2*tavl(l)-tmptab(j,l)-tmptab(j+1,l))/den1 ad2(l)=(2*tavl(l)-tmptab(j-1,l)-tmptab(j+1,l))/den2 END DO RETURN END SUBROUTINE settabindx !---------------------------------------------------------------------------- ! PURPOSE: Compute radiances (in mw/m2/str/cm-1) and derivatives of radiances ! with respect to atmospheric and surface parameters !---------------------------------------------------------------------------- SUBROUTINE ossrad(tautot,abso,dtaudtmp,tavl,dtu,dtl,dwqu,dwql,nmols,& imols,xG,EmRf,bs_sun,vn,N1,N2,sun,umu,umu0,lookup,rad,xkt,xkEmRf) !---Input variables LOGICAL :: sun,lookup INTEGER :: N1,N2 REAL :: umu,umu0,bs_sun,vn INTEGER*2 :: NmolS,ImolS(NmolS) REAL, DIMENSION(:) :: xG,EmRf,tautot,dtaudtmp,tavl,dtu,dtl REAL, DIMENSION(:,:) :: abso,dwqu,dwql !---Output variables: REAL :: rad REAL, DIMENSION(:) :: xkt,xkEmRf !---Local variables: INTEGER :: l,ks,k,ixoff REAL :: sec,secdif,sec0,em,sun_rsfc,bs,dbs,sumtau_dwn REAL :: tausfc,radsun,draddrsfc,draddtau_sun,draddemis REAL :: draddtskn,sumtau0_up,dtran,txsun,rsfc REAL, DIMENSION(MxLev) :: txdn,txup,bbar,dbbar,draddtmp,draddtau,drdw sec = 1/umu sec0 = 1/umu0 secdif = sec em = EmRf(1) sun_rsfc = EmRf(2) !----------------------------------------------------------------------- ! Compute Planck function and its derivative wrt temperature !----------------------------------------------------------------------- CALL planck_int(vn,tavl,xG(ITskinG),n2,bbar,dbbar,bs,dbs) !----------------------------------------------------------------------- ! Compute transmittance profile along viewing path down to surface !----------------------------------------------------------------------- txdn(N1) = 1. sumtau_dwn = 0. DO l=N1,N2 sumtau_dwn = sumtau_dwn+tautot(l)*sec txdn(l+1) = EXP(-sumtau_dwn) END DO tausfc = txdn(N2+1) !---Initialize radiance and derivative arrays rad = 0. radsun = 0. draddrsfc = 0. draddtau_sun = 0. draddtmp(1:N2) = 0. draddtau(1:N2) = 0. draddemis = 0. draddtskn = 0. !----------------------------------------------------------------------- ! 1- Downwelling thermal radiance calculation: !----------------------------------------------------------------------- IF(tausfc.GT.1.e-06)THEN txup(N2+1) = tausfc sumtau0_up = 0. DO l=N2,1,-1 sumtau0_up = sumtau0_up+tautot(l) txup(l) = EXP(-(sumtau_dwn+sumtau0_up*secdif)) END DO DO l=1,N2 dtran = txup(l+1)-txup(l) !---Derivative of downwelling emission wrt to temperature and constituents draddtmp(l) = dtran*dbbar(l) draddtau(l) = (txup(l)*bbar(l)-rad)*secdif rad = rad+dtran*bbar(l) ENDDO END IF !----------------------------------------------------------------------- ! 2- Add solar component !----------------------------------------------------------------------- IF(sun)THEN txsun = EXP(-(sumtau_dwn+sumtau0_up*sec0)) !---Derivatives wrt sfc solar reflectance draddrsfc = txsun*bs_sun*umu0/pi radsun = sun_rsfc*draddrsfc draddtau_sun = -radsun*(sec+sec0) END IF !----------------------------------------------------------------------- ! 3- Add Surface terms !----------------------------------------------------------------------- rsfc = (1.-em) !---Derivatives wrt emissivity and sfc skin temperature: draddemis = tausfc*bs-rad draddtskn = em*tausfc*dbs rad = rad*rsfc+em*tausfc*bs+radsun !----------------------------------------------------------------------- ! 4- Upwelling thermal radiance calculation !----------------------------------------------------------------------- DO l=N2,1,-1 dtran = txdn(l)-txdn(l+1) draddtau(l) = draddtau(l)*rsfc+(txdn(l+1)*bbar(l)-rad)*sec+draddtau_sun draddtmp(l) = draddtmp(l)*rsfc+dtran*dbbar(l)+draddtau(l)*dtaudtmp(l) rad = rad+dtran*bbar(l) ENDDO !----------------------------------------------------------------------- ! Compute level derivatives and and map to array XKT: !----------------------------------------------------------------------- xkt=0. ! Air Temperature xkt(ITempG+1:ITempG+N2-1)=draddtmp(2:N2)*dtu(2:N2)+draddtmp(1:N2-1)*dtl(1:N2-1) !---level 1 xkt(ITempG)=draddtmp(1)*dtu(1) !---bottom level (N2+1) xkt(ITempG+N2)=draddtmp(N2)*dtl(N2) !---Molecular concentrations DO kS=1,nmolS k=ImolS(kS) IXoff=ImolInd(k)-1 drdw(1:N2)=draddtau(1:N2)*abso(k,1:N2) xkt(IXoff+2:IXoff+N2)=drdw(2:N2)*dwqu(2:N2,k)+drdw(1:N2-1)*dwql(1:N2-1,k) !---level 1 xkt(IXoff+1)=drdw(1)*dwqu(1,k) !---bottom level (N2+1) xkt(IXoff+N2+1)=drdw(N2)*dwql(N2,k) END DO !---Surface terms xkt(ITskinG)=draddtskn !--Tskin xkEmRf(1)=draddemis xkEmRf(2)=draddrsfc RETURN END SUBROUTINE ossrad SUBROUTINE planck_int(vn,t,ts,nlay,b,db,bs,dbs) !---Input variables INTEGER :: nlay REAL :: vn,ts REAL, DIMENSION(:) :: t !---Output variables REAL :: bs,dbs REAL, DIMENSION(:) :: b,db !---Local variables REAL, PARAMETER :: dvint=10. REAL :: bs1,bs2,dbs1,dbs2,ars,brs,ads,bds REAL :: v1,v2,v2n REAL, DIMENSION(mxlay) :: b1,b2,db1,db2,ar,br,ad,bd COMMON//v2,b2,db2,bs2,dbs2 SAVE ar,br,ad,bd,ars,brs,ads,bds IF(vn.GE.v2)THEN v2n = v2+dvint IF(vn.GE.v2n)THEN CALL planck_set(vn,t,ts,nlay) v2 = vn v2n = v2+dvint END IF v1 = v2 b1(1:nlay) = b2(1:nlay) db1(1:nlay) = db2(1:nlay) bs1 = bs2 dbs1 = dbs2 v2 = v2n CALL planck_set(v2,t,ts,nlay) ar(1:nlay) = (b2(1:nlay)-b1(1:nlay))/dvint br(1:nlay) = (b1(1:nlay)*v2-b2(1:nlay)*v1)/dvint ad(1:nlay) = (db2(1:nlay)-db1(1:nlay))/dvint bd(1:nlay) = (db1(1:nlay)*v2-db2(1:nlay)*v1)/dvint ars = (bs2-bs1)/dvint brs = (bs1*v2-bs2*v1)/dvint ads = (dbs2-dbs1)/dvint bds = (dbs1*v2-dbs2*v1)/dvint END IF b(1:nlay) = vn*ar(1:nlay)+br(1:nlay) db(1:nlay) = vn*ad(1:nlay)+bd(1:nlay) bs = vn*ars+brs dbs = vn*ads+bds RETURN END SUBROUTINE planck_int SUBROUTINE planck_set(vn,t,ts,nlay) !---Input variables INTEGER :: nlay REAL :: vn,ts REAL, DIMENSION(:) :: t !---Local variables INTEGER :: l REAL, DIMENSION(mxlay) :: b2,db2 REAL :: v2,bs2,dbs2 COMMON//v2,b2,db2,bs2,dbs2 v2=vn DO l=1,nlay CALL planck(vn,t(l),b2(l),db2(l)) END DO CALL planck(vn,ts,bs2,dbs2) RETURN END SUBROUTINE planck_set !---------------------------------------------------------------------------- ! PURPOSE: Calculates planck function and its derivative wrt temperature !---------------------------------------------------------------------------- SUBROUTINE planck(vn,t,rad,draddt) REAL, PARAMETER :: h=6.626176e-27,c=2.997925e+10,b=1.380662e-16 REAL, PARAMETER :: c1 = 2.d0*h*c*c,c2 = h*c/b !---Input variables REAL :: t,vn !---Output variables REAL :: rad,draddt !---Local variables REAL :: f1,f2,f3 f1=c1*vn**3 f2=c2*vn/t f3=EXP(f2) rad=f1/(f3-1.) draddt=(rad*rad)*f2*f3/(f1*t) RETURN END SUBROUTINE planck !--------------------------------------------------------------------------------- ! PURPOSE: Given a node frequency vn, perform spectral interpolation to determine ! cloud optical properties at vn based on stored values at bounding ! hingepoints v1 and v2. !--------------------------------------------------------------------------------- SUBROUTINE optInt(vn,nSurf,tabs,tscat,asym) !---Input variables REAL*8, INTENT(IN) :: vn INTEGER, INTENT(IN) :: nSurf !---Output variables REAL, DIMENSION(:), INTENT(INOUT) :: tabs,tscat,asym !---Local variables INTEGER :: nlay,ip,ip1 REAL :: v1,dvint nlay = nSurf-1 if(vn >= v2o)then DO ip=ip2,nSpcPts-1 IF(spectPtsGrd(ip).GT.vn) EXIT END DO ip1=ip-1 ip2=ip v1=spectPtsGrd(ip1) v2o=spectPtsGrd(ip2) dvint=v2o-v1 atabs(1:nlay) = (tabs2d(1:nlay,ip2)-tabs2d(1:nlay,ip1)) / dvint btabs(1:nlay) = ((v2o*tabs2d(1:nlay,ip1))-(v1*tabs2d(1:nlay,ip2))) / dvint atscat(1:nlay) = (tscat2d(1:nlay,ip2)-tscat2d(1:nlay,ip1)) / dvint btscat(1:nlay) = ((v2o*tscat2d(1:nlay,ip1))-(v1*tscat2d(1:nlay,ip2))) / dvint aasym(1:nlay) = (asym2d(1:nlay,ip2)-asym2d(1:nlay,ip1)) / dvint basym(1:nlay) = ((v2o*asym2d(1:nlay,ip1))-(v1*asym2d(1:nlay,ip2))) / dvint endif tabs(1:nlay) = vn * atabs(1:nlay) + btabs(1:nlay) tscat(1:nlay) = vn * atscat(1:nlay) + btscat(1:nlay) asym(1:nlay) = vn * aasym(1:nlay) + basym(1:nlay) RETURN END SUBROUTINE optInt !-- SUBROUTINE OSS_DESTROY_IR() if (allocated(SfGrd)) deallocate(SfGrd) if (allocated(nch_ir)) deallocate(nch_ir) if (allocated(vFreq)) deallocate(vFreq) if (allocated(coef_ir)) deallocate(coef_ir) if (allocated(cFreq)) deallocate(cFreq) if (allocated(pref)) deallocate(pref) if (allocated(sunrad)) deallocate(sunrad) if (allocated(nsmpPerChan)) deallocate(nsmpPerChan) if (allocated(ichmap_ir)) deallocate(ichmap_ir) if (allocated(Nsmp2NchanMap)) deallocate(Nsmp2NchanMap) if (allocated(NmolS)) deallocate(NmolS) if (allocated(ImolS)) deallocate(ImolS) if (allocated(Tmptab)) deallocate(Tmptab) if (allocated(Wvptab)) deallocate(Wvptab) if (allocated(kfix_ir)) deallocate(kfix_ir) if (allocated(kh2o_ir)) deallocate(kh2o_ir) if (allocated(dkh2o_ir)) deallocate(dkh2o_ir) if (allocated(kvar_ir)) deallocate(kvar_ir) END SUBROUTINE OSS_DESTROY_IR END MODULE oss_ir_module_scat