afrt.cpp

Go to the documentation of this file.

/*
 * afrt.cpp
 *
 *  Created on: Mar 22, 2019
 *  Author: Samuel Anderson
 */
 
#include <math.h>
#include <afrt.h>
#include <array>
 
/**************************************************************************
 * NAME: Afrt()
 *
 * DESCRIPTION: Constructor
 *
 *************************************************************************/
 
Afrt::Afrt()
{
}
 
/**************************************************************************
 * NAME: ~Afrt()
 *
 * DESCRIPTION: Destructor
 *
 *************************************************************************/
 
Afrt::~Afrt()
{
}
 
/**************************************************************************
 * NAME: Initialize()
 *
 * DESCRIPTION: Initialize object class should it be required
 *
 *************************************************************************/
 
int Afrt::initialize()
{
    return RT_SUCCESS;
}
 
/**************************************************************************
 * NAME: ocn()
 *
 * DESCRIPTION: Compute ocean surface reflectance
 *
 *************************************************************************/
 
int Afrt::ocn()
{
    return RT_SUCCESS;
}
 
/**************************************************************************
 * NAME: rt1()
 *
 * DESCRIPTION: First part of the radiative transfer program. Compute
 * optical thickness profiles.
 *
 *************************************************************************/
 
int Afrt::rt1()
{
    return RT_SUCCESS;
}
 
/**************************************************************************
 * NAME: rt1()
 *
 * DESCRIPTION: First part of the radiative transfer program. Compute
 * optical thickness profiles.
 *
 *************************************************************************/
 
int Afrt::rt2(int ilm, int isd, int itau, int iwnd,
                rt2_in* rt2in, rt2_out* rt2out,
                rt1_in* rt1in, rt1_out* rt1out,
                ocn_in* ocnin, ocn_out* ocnout,
                phs_in* phsin, phs_out* phsout)
{
    int status = 0;
 
    ifc = rt1out->ifc[0][0][0];
    iref = rt2in->iref;
    iglint = rt2in->iglint;
    itrans = rt2in->itrans;
    iwatr = rt2in->iwatr;
    ifoam = rt2in->ifoam;
    ipol = rt2in->ipol;
    nx = rt2in->inx;
    dphi = rt2in->dphi;
    alw = rt2in->albwat[ilm];
    nmu = 2*nx-2;
    nophi = 360.0/dphi;
    jpart = nophi/2+1;
    ddphi = dphi*D2R;
    conr = 3.0/(8.0*M_PI);
// rt1 inputs
    rho = rt1out->rho[0][0][0];
    nolyr = rt1out->nolyr[0][0][0];
    wvlth = rt1out->wvlth[0][0][0];
    tautot = rt1out->tautot[0][0][0];
    for (int i=0; i<nolyr; i++) {
        dtmm[i] = rt1out->dtmm[0][0][0][i];
        dtrr[i] = rt1out->dtrr[0][0][0][i];
        dtaa[i] = rt1out->dtaa[0][0][0][i];
        dtot[i] = rt1out->dtot[0][0][0][i];
    }
// phs inputs
    int olm = ilm-rt2in->ilm1;
    int osd = isd-rt2in->isd1;
    cnst = pow(wvlth,2.0)/(4.0*pow(M_PI,2.0)*phsout->qscat[olm][osd]);
    qsqt = phsout->qscat[olm][osd]/phsout->qext[olm][osd];
    for (int i=0; i<ntf; i++) {
        thd[i] = phsout->thd[olm][osd][i];
        for (int k=0; k<nstk; k++) {
            t[k][i] = phsout->t[olm][osd][k][i];
        }
    }
// ocn inputs
    xrw = ocnin->xr;
    xiw = ocnin->xi;
    vz = ocnin->v;
    for (int i=0; i<nsza; i++) {
        for (int j=0; j<nthe; j++) {
            for (int l=0; l<nphi; l++) {
                for (int k=0; k<16; k++) {
                    txx[i][j][l][k] = ocnout->txx[0][0][i][j][l][k];
                }
            }
        }
    }
 
    memset(&fnull, 0, sizeof(fdat));
    unit53.assign(nolyr+1, fnull);
    unit54.assign(nolyr+1, fnull);
    unit55.assign(nolyr+1, fnull);
    unit64.assign(nolyr+1, fnull);
    unit65.assign(nolyr+1, fnull);
    unit71.assign(nolyr+1, fnull);
    unit72.assign(nolyr+1, fnull);
 
    insz = 0;
// loop over solar zenith angle
    for (ksza = rt2in->ithe01; ksza <= rt2in->ithe02; ksza++) {
        int isza = ksza-rt2in->ithe01;
        msza[isza] = ksza;
        insz += 1;
        calb = 0.0;
        jpass = 0;
        for (int i=0; i<nmu; i++) {
            c[i] = 1.0;
        }
        for (int i=0; i < 2*nmu; i++) {
            for (int k=1; k<nstk; k++ ) {
                pp[k][i] = 0.0;
                qq[k][i] = 0.0;
            }
        }
        for (int i=0; i<2; i++) {
            ei[i] = 0.50;
            ei[i+2] = 0.0;
        }
        for (int i = 0; i < nx; i++) {
            rmu[i] = rt2in->rx[i];
            rmu[nmu-i] = 180.0 - rt2in->rx[i];
        }
        for (int i = 0; i < nmu+1; i++) {
            cmu[i] = cos(rmu[i]*D2R);
        }
        for (int i = 0; i < nmu; i++) {
            double xmu = (rmu[i] + rmu[i+1])/2.0;
            the[i] = xmu;
            if (i > nx-2) {
                the[i] = the[nmu-1-i];
            }
            if (abs(xmu - rt2in->the0in[ksza-1]) <= 1.0e-3) {
                amuo = cos(xmu*D2R);
                kkx = i;
            }
            cosmu[i] = cos(xmu*D2R);
            sinmu[i] = sin(xmu*D2R);
        }
        for (int l=0; l<nophi+1; l++) {
            phi[l] = l*dphi;
        }
        for (int i = 0; i < jpart; i++) {
            costh[i] = cos(((double)i)*dphi*D2R);
            sinth[i] = sin(((double)i)*dphi*D2R);
        }
        for (int i = 0; i < nmu; i++) {
            dcmu2[i] = (pow(cmu[i], 2.0) - pow(cmu[i+1], 2.0))/2.0;
            dcmu[i] = cmu[i] - cmu[i+1];
        }
// determine an average value of phase matrix in the forward/
// backward direction
        if (ifc == 1) {
            hump(cnst, t, pp, nmu, dphi+0.001, rmu, thd, 0);
            hump(cnst, t, qq, nmu, dphi+0.001, rmu, thd, 1);
        }
        if (iref == 2) {
            anglw();
        }
// compute the attenuation parameters
// determine the total optical thickness at every level in the atmosphere
        taupl[0] = 0.0;
        for (int i=0; i<nolyr; i++) {
            taupl[i+1] = taupl[i] + dtot[i];
        }
// compute the solar attenuation up to the middle as well as
// bottom of the layer (plane parallel case)
        efactb[0] = 1.0;
        for (int i = 0; i < nolyr; i++) {
            efact[i] = exp(-(0.50*dtot[i]+taupl[i])/amuo);
            efactb[i+1] = exp(-taupl[i+1]/amuo);
        }
        eo[0] = 0.5*exp(-taupl[nolyr]/amuo);
        eo[1] = 0.5*exp(-taupl[nolyr]/amuo);
        totl[0] = 1.0e-10;
        for (int i = 0; i < nolyr; i++) {
            totl[i+1] = totl[i] + dtot[i];
        }
        for (int i = 0; i < nolyr; i++) {
            for (int j=0; j<nmu; j++) {
                emdtm[j][i] = exp(-dtot[i]/abs(cosmu[j]));
                emtm[j][i] = exp(-(totl[i+1]-0.5*dtot[i])/abs(cosmu[j]));
            }
        }
// compute the scattering phase matrices (rayleigh/aerosol)
// compute phase matrices for first scatter
        int kk = kkx;
        for (int ii = 0; ii < nmu; ii++) {
            matrx(ii, kk);
            for (int i = 0; i < jpart; i++) {
                for (int k = 0; k < 32; k++) {
                    pfs[ii][i][k] = p[i][k];
                }
            }
        }
// compute phase matrices for second scatter
        for (int i = 0; i < nx-1; i++) {
            qsp[i] = 0.0;
        }
        for (int ii = 0; ii < nmu; ii++) {
            for (int kk = 0; kk < nmu; kk++) {
                matrx(ii, kk);
                for (int i = 0; i < jpart; i++) {
                    for (int k = 0; k < 32; k++) {
                        ppin[kk][ii][i][k] = p[i][k];
                    }
                }
            }
            if (ifc != 0) {
                ee[0] = 0.5;
                ee[1] = 0.5;
                for (int k = 0; k < nx-1; k++) {
                    double qspp = 0.0;
                    for (int m = 0; m < jpart; m++) {
                        double xfot = cnst*(ee[0]*(ppin[k][ii][m][0] +
                                        ppin[k][ii][m][3]) + ee[1]*
                                       (ppin[k][ii][m][1] + ppin[k][ii][m][5]));
                        if (m==0 || m==jpart-1) {
                            qspp = qspp + xfot;
                        } else {
                            qspp = qspp + 2.0*xfot;
                        }
                    }
                    qsp[k] = qsp[k] + qspp*dcmu[ii]*ddphi;
                }
            }
        }
        if (ifc != 0) {
            for (int i = 0; i < nx-1; i++) {
                c[i] = 1.0/qsp[i];
                c[nmu-i] = c[i];
            }
        }
// begin radiative transfer calculations
        int ibgn = 0;
        int iend = 1;
        if (insz == 1 && (iref == 0 || itrans == 1)) {
            iend = 2;
        }
// illumination from the top of the atmosphere(kzz = 0)
// illumination from the bottom of the atmosphere(kzz = 1)
        for (kzz = ibgn; kzz < iend; kzz++) {
            jpass = 0;
            int nump = rt2in->inpass1;
            double d1px[nphi][nsza][2];
            double d2px[nphi][nsza][2];
            for (int k = 0; k < 2; k++) {
                for (int i = 0; i < nx-1; i++) {
                    for (int j = 0; j < jpart; j++) {
                        d1px[j][i][k] = 0.0;
                        d2px[j][i][k] = 0.0;
                    }
                }
            }
            if (insz == 1 && (iref == 0 || itrans == 1)) {
                nump = rt2in->inpass2;
            }
            if (kzz == 0) {
                single_down();
            }
            if (kzz == 1 && iref == 0) {
                // initialize unit 55
                unit55.assign(nolyr+1, fnull);
                if(iref==0) {
                    for (int i=nx-1; i<nmu; i++) {
                        for (int j=0; j<jpart; j++) {
                            for (int k=0; k<2; k++) {
                                unit55[nolyr].f[j][i][k] = 0.5/M_PI;
                                unit55[nolyr].f[j][i][k+2] = 0.0;
                            }
                        }
                    }
                }
                if(iref==1 || iref==2 || iref==3) {
                    fio = unit53[nolyr];
                }
                single_up();
            }
            if (kzz == 1 && itrans == 1) {
//              focn_below();
//              snglup();
            }
            getrad();
            for (int ka = 0; ka < nump; ka++) {
// initialize the disk units
                kdx = ka - (ka/2)*2;
                if (kdx == 0) {
                    irad = &unit64;
                    iwrt = &unit65;
                } else {
                    irad = &unit65;
                    iwrt = &unit64;
                }
                if (kzz == 0) {
                    multp1d();
                } else if (kzz == 1 && iref == 0) {
                    multp2d();
                } else if (kzz == 1 && itrans == 1) {
//                  multp2d_trans();
                }
                jpass++;
                getrad();
                double d3 = 1.0;
                if (jpass > 2) {
                    for (int i=0; i<nx-1; i++ ) {
                        for (int j=0; j<jpart; j++) {
                            d1px[j][i][0]=d1px[j][i][1];
                            d2px[j][i][0]=d2px[j][i][1];
                        }
                    }
                    for (int i=0; i<nx-1; i++) {
                        for (int j=0; j<jpart; j++) {
                            double r1=geoiup[j][i][jpass-3];
                            double r2=geoiup[j][i][jpass-2];
                            double r3=geoiup[j][i][jpass-1];
                            double rat=(r3-r2)/(r2-r1);
                            double r4=r1+(r2-r1)/(1.0-rat);
                            d1px[j][i][1]=r4;
                        }
                    }
                    for (int i=0; i<nx-1; i++) {
                        for (int j=0; j<jpart; j++) {
                            double u1=geoidn[j][i][jpass-3];
                            double u2=geoidn[j][i][jpass-2];
                            double u3=geoidn[j][i][jpass-1];
                            double uat=(u3-u2)/(u2-u1);
                            double u4=u1+(u2-u1)/(1.0-uat);
                            d2px[j][i][1]=u4;
                        }
                    }
                    double df1[nphi][nsza];
                    double df2[nphi][nsza];
                    // compute the difference
                    for (int i=0; i<nx-1; i++ ) {
                        for (int j=0; j<jpart; j++) {
                            df1[j][i]=100.0*(d1px[j][i][1]-d1px[j][i][0])/d1px[j][i][1];
                            df2[j][i]=100.0*(d2px[j][i][1]-d2px[j][i][0])/d2px[j][i][1];
                        }
                    }
                    // compute the max value
                    double xmax1 = 0.0;
                    double xmax2 = 0.0;
                    for (int i=0; i<nx-1; i++ ) {
                        for (int j=0; j<jpart; j++) {
                            if(df1[j][i]>xmax1) xmax1=df1[j][i];
                            if(df2[j][i]>xmax2) xmax2=df2[j][i];
                        }
                    }
                    d3 = max(xmax1, xmax2);
                }
                unit71[jpass] = iwrt->at(0);
                unit72[jpass] = iwrt->at(nolyr);
                if (jpass > 2) {
                    if (d3 <= 0.1 || jpass >= 20) break;
                }
            }
            geocor();
            //      illumination from bottom is done only once
            if (kzz == 1) break;
        }
    }
    for (int i=0; i<(rt2in->ithe02-rt2in->ithe01+1); i++) {
        int m=msza[i];
        rt2out->tma[0][0][0][0][i]=rt2in->the0in[m];
        rt2out->tmb[0][0][0][0][i]=fdirc[m];
        rt2out->tmc[0][0][0][0][i]=sbarz[m];
        rt2out->tmfd[0][0][0][0][i]=fdown[m];
        rt2out->tmfu[0][0][0][0][i]=fup[m];
        rt2out->tms[0][0][0][0][i]=oalb[m];
        rt2out->tmg[0][0][0][0][i]=albtdr[m];
        rt2out->tmh[0][0][0][0][i]=albtdf[m];
        rt2out->tmp[0][0][0][0][i]=albtrf[m];
        rt2out->tmq[0][0][0][0][i]=albwl[m];
        rt2out->tmpp[0][0][0][0][i]=albwl[m]/(fdirc[m]+fdown[m]);
        rt2out->tmqq[0][0][0][0][i]=oalb[m]-rt2out->tmpp[0][0][0][0][i];
        rt2out->tmrr[0][0][0][0][i]=albrdr[m];
        rt2out->tmss[0][0][0][0][i]=albrdf[m];
        for (int is=0; is<jpart; is++) {
            for (int ir=0; ir<nx-1; ir++) {
                for (int s=0; s<nstk; s++) {
                    rt2out->xzerod[0][0][0][0][i][ir][is][s] = xzerod[is][ir][m][s];
                    rt2out->xzeroz[0][0][0][0][i][ir][is][s] = xzeroz[is][ir][m][s];
                }
            }
        }
    }
 
    return status;
}
 
/**************************************************************************
 * NAME: hump()
 *
 * DESCRIPTION: program to compute pbar for each mup when mu=mup and delfi=0
 *
 *************************************************************************/
 
int Afrt::hump(double cnst, double t[][ntf], double tp[][nlyr],
        int nm, double lp, double rmu[], double thd[], int ixy)
{
    double p[nstk][nstk];
    memset(p, 0, nstk*nstk*sizeof(double));
    double dp=0.10;
    double dt=0.10;
    double fi1=lp/2;
    int ni1=fi1/dp+1;
    dp=fi1/((double) ni1);
    double costh[nlyr];
    double sinth[nlyr];
    for (int k=0; k<ni1; k++) {
        double angle;
        if(ixy!=0) {
            angle=(180.0-lp/2.0+((double)k)*dp+dp/2.0)*D2R;
        } else {
            angle=(((double)k)*dp+dp/2.0)*D2R;
        }
        sinth[k] = sin(angle);
        costh[k] = cos(angle);
    }
    double deldeg = dt;
    double dcmu[nlyr];
    for (int kl=0; kl<nm/2; kl++) {
        double deg = rmu[kl+1] - rmu[kl];
        double nodeg = deg/deldeg;
        deldeg = deg / nodeg;
        double thetd = (rmu[kl]+rmu[kl+1])/2.0;
        double amup = cos(thetd*D2R);
        for (int i=0; i<nodeg+1; i++) {
            dcmu[i] = cos((thetd-deg/2.0+((double)i)*deldeg)*D2R)-
                       cos((thetd-deg/2.0+((double)(i+1))*deldeg)*D2R);
        }
        double ddmu = cos(rmu[kl]*D2R)-cos(rmu[kl+1]*D2R);
        double amup2=pow(amup,2.0);
        for (int i=0; i<nstk; i++) {
            int k=kl*nstk+i;
            for (int j=0; j<nstk; j++) {
                tp[j][k]=0.0;
            }
        }
        for (int i=0; i<nodeg; i++) {
            double amu=cos((thetd-(deg-deldeg)/2.0+((double)i)*deldeg)*D2R);
            if(ixy!=0) {
                amu=-amu;
            }
            memset(p, 0, nstk*nstk*sizeof(double));
            double amusq=pow(amu,2.0);
            double amumu=amu*amup;
            double anunu = sqrt((1.0-amup2)*(1.0-amusq));
            for (int j=0; j<ni1; j++) {
                double copsi=anunu+amumu*costh[j];
                double sifi=copsi*costh[j];
                double sfisq=pow(sinth[j],2.0);
                double cfisq=pow(costh[j],2.0);
                double sisq=pow(copsi,2.0);
                //double a=amup*costh[j];
                //double b=amup*copsi;
                //double c=amu*costh[j];
                //double d=amu*copsi;
                double e=sfisq*(amusq+amup2);
                double g=amumu*sfisq;
                double h=sfisq*(amusq-amup2);
                //double x=amumu+anunu*costh[j];
                double yz= acos(amumu+anunu*costh[j])/D2R;
                int m=10.0*yz;
                if(m==0) m=0;
                double t1, t2, t3, t4;
                xntpln(yz,thd[m],thd[m+1],t[0][m],t[0][m+1],t1);
                xntpln(yz,thd[m],thd[m+1],t[1][m],t[1][m+1],t2);
                xntpln(yz,thd[m],thd[m+1],t[2][m],t[2][m+1],t3);
                xntpln(yz,thd[m],thd[m+1],t[3][m],t[3][m+1],t4);
                p[2][2]=p[2][2]+(sifi-g)*(t1+t2)+(cfisq+sisq-e)*t3;
                p[1][1]=p[1][1]+sisq*t1+cfisq*t2+2.0 *sifi * t3;
                p[1][0]=p[1][0]+sfisq*(amup2*t1+amusq*t2+2.0 * amumu *t3 );
                p[0][1]=p[0][1]+sfisq*(amup2*t2+amusq*t1+2.0 * amumu *t3);
                p[0][0]=p[0][0]+cfisq*t1+sisq*t2+2.0 * sifi *t3;
                p[2][3]=p[2][3]+(sisq-cfisq-h)*t4;
                p[3][2]=p[3][2]+(cfisq-sisq-h)*t4;
                p[3][3]=p[3][3]+(sifi+g)*(t1+t2)+(cfisq+sisq+e)*t3;
            }
            for (int k=0; k<nstk; k++) {
                int kk = kl*4+k;
                for (int l=0; l<nstk; l++) {
                    //double ttpp1=cnst*p[0][k]*dcmu[i]*dp/fi1;
                    //double ttpp2=cnst*p[1][k]*dcmu[i]*dp/fi1;
                    //double ttpp=0.5*(ttpp1+ttpp2);
                    tp[l][kk]=p[l][k]*dcmu[i]*dp/fi1+tp[l][kk];
                }
            }
        }
        for (int k=0; k<nstk; k++) {
            int kk = kl*4+k;
            for (int l=0; l<nstk; l++) {
                tp[l][kk]=tp[l][kk]/ddmu;
            }
        }
        //int id = kl*4;
        //double ppbar=cnst*(tp[0][id]+tp[1][id]+tp[0][id+1]+tp[1][id+1])*0.5;
    }
 
    return RT_SUCCESS;
}
 
/**************************************************************************
 * NAME: anglw()
 *
 * DESCRIPTION: compute angles of refraction and sin/cos assuming flat ocean
 *
 *************************************************************************/
 
int Afrt::anglw()
{
    double rmuw[2*nthe];
    for (int i=0; i<nx; i++) {
        double cnuw = sqrt(1.0-pow(cmu[i],2.0))/xrw;
        rmuw[i]=asin(cnuw)/D2R;
    }
// compute the center of the bin and its sin/cos
    for (int i=0; i<nx-1; i++) {
        double thew = 0.5*(rmuw[i]+rmuw[i+1]);
        double dnxy = pow(xrw,2.0)+pow(xiw,2.0);
        double yrw=xrw/dnxy;
        double yiw=xiw/dnxy;
        frs(the[i],xrw,xiw,rfair[i]);
        frs(thew,yrw,yiw,rfwat[i]);
    }
    return RT_SUCCESS;
}
 
/**************************************************************************
 * NAME: frs()
 *
 * DESCRIPTION: compute the reflectivity of a flat ocean
 *
 *************************************************************************/
 
int Afrt::frs(double xx, double xr, double xi, double& rfsea)
{
    double cost = cos(xx*D2R);
    // calculate the elements of fresnel matrix
    double sinsq = (1.0-pow(cost,2.0));
    double sinsq2 = pow(sinsq,2.0);
    double a = (pow(xr,2.0)-pow(xi,2.0)-1.0 + pow(cost,2.0));
    double b = (-2.0 * xr * xi);
    double r = sqrt(pow(a,2.0) + pow(b,2.0));
    double tmr = pow(cost,2.0)*r;
    double xp = sqrt(2.0*r + 2.0*a);
    double xmm = (2.0*r - 2.0*a);
    if(xmm < 0.0 && xmm > -9.0) xmm = -xmm;
    double xm=sqrt(xmm);
    double dnmr1 = sinsq2 + tmr + cost*sinsq*xp;
    double qrmu = (sinsq2 - tmr)/dnmr1;
    double qimu = (cost*sinsq*xm)/dnmr1;
    double dnmr2 = pow(cost,2.0) + r + cost*xp;
    double rrr = (pow(cost,2.0) - r)/dnmr2;
    double rri = cost*xm/dnmr2;
    double rer = qrmu*rrr - qimu*rri;
    double rei = qimu*rrr + qrmu*rri;
    double r11 = pow(rer,2.0) + pow(rei,2.0);
    double r22 = pow(rrr,2.0) + pow(rri,2.0);
    //double r33 = rer*rrr + rei*rri;
    //double r34 = rri*rer - rei*rrr;
    rfsea = 0.5*(r11+r22);
 
    return RT_SUCCESS;
}
 
/**************************************************************************
 * NAME: matrx()
 *
 * DESCRIPTION: Compute scattering matrices
 *
 *************************************************************************/
 
int Afrt::matrx(int ii, int kk)
{
    amusq=pow(cosmu[ii],2.0);
    amumu=cosmu[ii]*cosmu[kk];
    amups=pow(cosmu[kk],2.0);
    for (int l=0; l<jpart; l++) {
        copsi=sinmu[ii]*sinmu[kk]+cosmu[ii]*cosmu[kk]*costh[l];
        copsq=pow(copsi,2.0);
        cfisq=pow(costh[l],2.0);
        copcs=copsi*costh[l];
        sfisq=pow(sinth[l],2.0);
        int iq, ip, ir, jp, mn;
        if((cosmu[kk]-cosmu[ii])==0) {
            if(l==0) {
                if((ii-nx)<0) {
                    iq=4*ii;
                } else {
                    iq=4*(nmu-ii-1);
                }
                mn=0;
                for (int ik=0; ik<16; ik++) {
                    jp=4*(ik/4);
                    ip=iq+jp/4;
                    ir=ik-jp;
                    p[l][ik]=pp[ir][ip];
                }
                depol(ii, kk, l);
                if((ii-nx)<0) {
                    if(l==0) {
                        mats(ii, kk, jpart);
                    }
                }
            } else {
                    mats(ii, kk, l);
            }
        } else if((cosmu[ii]+cosmu[kk])==0) {
            if((l-jpart)==0) {
                if((ii-nx)<0) {
                    iq=4*ii;
                } else {
                    iq=4*(nmu-ii-1);
                }
                mn=jpart*32-31;
                for (int ik=mn; ik<=mn+15; ik++) {
                    jp=4*((ik-mn)/4);
                    ip=iq+jp/4;
                    ir=ik-jp-mn;
                    p[l][ik-mn]=qq[ir][ip];
                }
                mn=mn-1;
                depol(ii, kk, l);
                if((ii-nx)<0) {
                    if(l==0) {
                        mats(ii, kk, jpart);
                    }
                }
            } else {
                mats(ii, kk, l);
            }
        } else {
            mats(ii, kk, l);
        }
    }
 
    return RT_SUCCESS;
}
/**************************************************************************
 * NAME: mats()
 *
 * DESCRIPTION: Compute scattering matrices subroutine
 *
 *************************************************************************/
 
int Afrt::mats(int ii, int kk, int l)
{
    double x = cosmu[ii] * cosmu[kk] + sinmu[ii] * sinmu[kk] * costh[l];
    if(x>1.0) x = 1.0;
    if(x<-1.0) x = -1.0;
    if(abs(x) < 1.0e-6) x = 0.0;
    double tf = 10.0*(acos(x)/D2R);
    double tfm1 = floor(tf + 0.01);
    double tfp1 = tfm1 + 1.0;
    int mt = 0;
    if ((l - jpart) <= 0) {
        mt = tfm1;
    }
    double tmt1, tmt2, tmt3, tmt4;
    xntpln(tf, tfm1, tfp1, t[0][mt], t[0][mt+1], tmt1);
    xntpln(tf, tfm1, tfp1, t[1][mt], t[1][mt+1], tmt2);
    xntpln(tf, tfm1, tfp1, t[2][mt], t[2][mt+1], tmt3);
    xntpln(tf, tfm1, tfp1, t[3][mt], t[3][mt+1], tmt4);
    p[l][0] = cfisq*tmt1 + copsq*tmt2 + 2.0*copcs*tmt3;
    p[l][1] = sfisq*(amups*tmt1 + amusq*tmt2 + 2.0*amumu*tmt3);
    p[l][2] = -sinth[l]*(cosmu[kk]*costh[l]*tmt1 + cosmu[ii]*copsi*tmt2 +
              (cosmu[kk]*copsi + cosmu[ii]*costh[l])*tmt3);
    p[l][3] = sinth[l]*(cosmu[kk]*copsi - cosmu[ii]*costh[l])*tmt4;
    p[l][4] = sfisq*(amups*tmt2 + amusq*tmt1 + 2.0*amumu*tmt3);
    p[l][5] = copsq * tmt1 + cfisq * tmt2 + 2.0 * copcs * tmt3;
    //double pmats = cnst*0.5*(p[l][0]+p[l][1]+p[l][4]+p[l][5])*4.0*M_PI;
    p[l][6]=sinth[l]*(cosmu[kk]*costh[l]*tmt2+cosmu[ii]*copsi*tmt1+
            (cosmu[kk]*copsi+cosmu[ii]*costh[l])*tmt3);
    p[l][7]=-p[l][3];
    p[l][8]=2.0*sinth[l]*(cosmu[ii]*costh[l]*tmt1+cosmu[kk]*copsi*tmt2+
            (cosmu[kk]*costh[l]+cosmu[ii]*copsi)*tmt3);
    p[l][9]=-2.0*sinth[l]*(cosmu[ii]*costh[l]*tmt2+cosmu[kk]*copsi*tmt1+
            (cosmu[kk]*costh[l]+cosmu[ii]*copsi)*tmt3);
    p[l][10]=(copcs-amumu*sfisq)*(tmt1+tmt2)+(cfisq+copsq-sfisq*(amusq+amups))*tmt3;
    p[l][11]=(cfisq-copsq-sfisq*(amusq-amups))*tmt4;
    p[l][12]=2.0*sinth[l]*(cosmu[ii]*copsi-cosmu[kk]*costh[l])*tmt4;
    p[l][13]=-p[l][12];
    p[l][14]=(copsq-cfisq-sfisq*(amusq-amups))*tmt4;
    p[l][15]=(copcs+amumu*sfisq)*(tmt1+tmt2)+(cfisq+copsq+sfisq*(amusq+amups))*tmt3;
    depol(ii, kk, l);
 
    return RT_SUCCESS;
}
 
/**************************************************************************
 * NAME: depol()
 *
 * DESCRIPTION: Compute scattering matrices subroutine
 *
 *************************************************************************/
 
int Afrt::depol(int ii, int kk, int l)
{
    p[l][16] = copsq;
    p[l][17] = amusq *sfisq;
    p[l][18] = -cosmu[ii]*copsi*sinth[l];
    p[l][19] = 0.0;
    p[l][20] = amups*sfisq;
    p[l][21] = cfisq;
    p[l][22] = cosmu[kk]*costh[l]*sinth[l];
    p[l][23] = 0.0;
    p[l][24] = 2.0*cosmu[kk]*copsi*sinth[l];
    p[l][25] = -2.0*cosmu[ii]*costh[l]*sinth[l];
    p[l][26] = -amumu*sfisq + copcs;
    p[l][27] = 0.0;
    p[l][28] = 0.0;
    p[l][29] = 0.0;
    p[l][30] = 0.0;
    p[l][31] = amumu*sfisq + copcs;
    // apply molecular depolarization correction
    if (ipol == 1) {
        double gam = rho/(2.0-rho);
        double agm = (1.0-gam)/(1+2.0*gam);
        double bgm = gam/(1.0+2.0*gam);
        double cgm = (1.0-3.0*gam)/(1.0+2.0*gam);
        for (int i=0; i<3; i++) {
            for (int j=0; j<3; j++) {
                int k = 16+i*4+j;
                double dgm = bgm;
                if(i==2 || j==2) dgm = 0.0;
                p[l][k] = agm * p[l][k] + dgm;
            }
        }
        p[l][31] = p[l][31] + cgm;
    }
    return RT_SUCCESS;
}
 
/**************************************************************************
 * NAME: single_down()
 *
 * DESCRIPTION: Compute scattering matrices subroutine
 *
 *************************************************************************/
 
int Afrt::single_down()
{
// initialize the buffers
    ftmp = fnull;
// initialize all levels
    unit53.assign(nolyr+1, fnull);
    unit64.assign(nolyr+1, fnull);
// compute the downward diffuse radiation at each level
    double qfi[nstk];
    double gfot[nstk];
    for (int il=0; il<nolyr; il++) {
        fiib = unit53[il];
        fio = unit53[il+1];
        tmsl=dtmm[il]*qsqt*cnst;
        trsl=dtrr[il]*conr;
        for (int is=0; is<nstk; is++) {
            qfi[is]=ei[is]*efact[il]/dtot[il];
        }
        for (int it=0; it<(nx-1); it++) {
            for (int ip=0; ip<jpart; ip++) {
                for (int is=0; is<nstk; is++) {
                    gfot[is]=0.0;
                }
                for (int i=0; i<nstk; i++) {
                    for (int j=0; j<nstk; j++) {
                        int ij=i*4+j;
                        gfot[i]=gfot[i]+(pfs[it][ip][ij]*tmsl +
                                pfs[it][ip][ij+16]*trsl)*qfi[j];
                    }
                    ftmp.f[ip][it][i]=gfot[i]*(1.0-emdtm[it][il]);
                    ftmp.f[ip][it][2]= (i>=2) ? -ftmp.f[ip][it][2] : ftmp.f[ip][it][2];
                    ftmp.f[ip][it][3]= (i>=2) ? -ftmp.f[ip][it][3] : ftmp.f[ip][it][3];
                    fio.f[ip][it][i] = fiib.f[ip][it][i]*emdtm[it][il] +
                            ftmp.f[ip][it][i];
                }
            }
        }
        unit53[il+1] = fio;
        unit54[il] = ftmp;
    }
    if(iref==1) {
//        fltocn_new();
        unit53[nolyr] = fio;
    }
// save the reflected radiation
    ftmpb = fio;
// compute the upward diffuse radiation at each level
    for (int il=0; il<nolyr; il++) {
        int im=nolyr-il-1;
        fiib = unit53[im+1];
        fio = unit53[im];
        ftmp = unit54[im];
        tmsl=dtmm[im]*qsqt*cnst;
        trsl=dtrr[im]*conr;
        for (int is=0; is<nstk; is++) {
            qfi[is]=ei[is]*efact[im]/dtot[im];
        }
        for (int it=nx-1; it<nmu; it++) {
            for (int ip=0; ip<jpart; ip++) {
                for (int is=0; is<nstk; is++) {
                    gfot[is]=0.0;
                }
                for (int i=0; i<nstk; i++) {
                    for (int j=0; j<nstk; j++) {
                        int ij=i*4+j;
                        gfot[i]=gfot[i]+(pfs[it][ip][ij]*tmsl +
                                pfs[it][ip][ij+16]*trsl)*qfi[j];
                    }
                    ftmp.f[ip][it][i] = gfot[i]*(1.0-emdtm[it][im]);
                    ftmp.f[ip][it][2]= (i>=2) ? -ftmp.f[ip][it][2] : ftmp.f[ip][it][2];
                    ftmp.f[ip][it][3]= (i>=2) ? -ftmp.f[ip][it][3] : ftmp.f[ip][it][3];
                    fio.f[ip][it][i]=fiib.f[ip][it][i]*emdtm[it][im] +
                            ftmp.f[ip][it][i];
                }
            }
        }
// add the contribution from the specular reflection
/*
        for (int is=0; is<nstk; is++) {
            qfi[is]=eox[is]*efact[il]/dtot[il]
        }
 
        for (int it=nx-1; it<nmu; it++) {
            int itp=nmu-it;
            for (int ip=0; ip<jpart; ip++) {
                for (int is=0; is<nstk; is++) {
                    gfot[is]=0.0
                }
                for (int i=0; i<nstk; i++) {
                    for (int j=0; j<nstk; j++) {
                        int ij=i*4+j;
                        gfot[i]=gfot[i]+(pfs[itp][ip][ij]*tmsl +
                                pfs[itp][ip][ij+16]*trsl)*qfi[j];
                    }
                    ftmp.f[ip][it][i] = gfot[i]*(1.0-emdtm[itp][im]);
                    ftmp.f[ip][it][2]= (i>=2) ? -ftmp.f[ip][it][2]: ftmp.f[ip][it][2];
                    ftmp.f[ip][it][3]= (i>=2) ? -ftmp.f[ip][it][3]: ftmp.f[ip][it][3];
                    fio.f[ip][it][i]=fio.f[ip][it][i]+ftmp.f[ip][it][i]+
                            ftmpb.f[ip][it][i]*atnflx[it][il+1];
                }
            }
        }
*/
        unit53[im] = fio;
        unit54[im] = ftmp;
    }
// copy all records to unit 64 for multiple scattering calculations
    for (int il=0; il<=nolyr; il++) {
        unit64[il] = unit53[il];
    }
 
    return RT_SUCCESS;
}
 
/**************************************************************************
 * NAME: single_up()
 *
 * DESCRIPTION: Compute scattering matrices subroutine
 *
 *************************************************************************/
 
int Afrt::single_up()
{
    memset(&fio, 0, sizeof(fdat));
    memset(&ftmp, 0, sizeof(fdat));
    memset(&ftmpa, 0, sizeof(fdat));
// compute the upwelling single scattered diffused radiation
// save the upwelling diffuse radiation for level nolyr+1 in the buffer ftmpa
    fiib = unit55[nolyr];
    ftmpa = unit55[nolyr];
    for (int il=0; il<nolyr; il++) {
        int im=nolyr-il-1;
        tmsl=dtmm[im]*qsqt*cnst;
        trsl=dtrr[im]*conr;
        dlyr=dtot[im];
        for (int i=0; i<nmu; i++) {
            for (int j=0; j<jpart; j++) {
                for (int k=0; k<nstk; k++) {
                    ftmp.f[j][i][k]=0.0;
                    fiic.f[j][i][k]=ftmpa.f[j][i][k]*emtm[i][im];
                }
            }
        }
        mdiffn(nx-1, nmu-1, im, dlyr);
        for (int i=0; i<nmu; i++) {
            for (int j=0; j<jpart; j++) {
                for (int k=0; k<nstk; k++) {
                    fio.f[j][i][k] = fiib.f[j][i][k]*emdtm[i][im] +
                                    ftmp.f[j][i][k]*(1.0-emdtm[i][im]);
                }
            }
        }
        fio = unit55[im];
        fiib = fio;
    }
// compute the downward radiation
    fiib = unit55[0];
    for (int il=0; il<nolyr; il++) {
        fio = unit55[il+1];
        tmsl=dtmm[il]*qsqt*cnst;
        trsl=dtrr[il]*conr;
        dlyr=dtot[il];
        for (int i=0; i<(nx-1); i++) {
            for (int j=0; j<jpart; j++) {
                for (int k=0; k<nstk; k++) {
                    fiic.f[j][i][k]=ftmpa.f[j][i][k]*emtm[i][il];
                }
            }
        }
        mdiffn(0, nx-2, il, dlyr);
        for (int i=0; i<(nx-1); i++) {
            for (int j=0; j<jpart; j++) {
                for (int k=0; k<nstk; k++) {
                    fio.f[j][i][k]=fiib.f[j][i][k]*emdtm[i][il] +
                                    ftmp.f[j][i][k]*(1.0-emdtm[i][il]);
                }
            }
        }
        fiib = fio;
        unit55[il+1] = fio;
    }
// if iref=0 then copy all records to unit 64 for multiple scattering
// calculations
    if(iref == 0 || itrans == 1) {
        for (int il=0; il<nolyr; il++) {
            unit64[il] = unit55[il];
        }
        eo[0]=0.0;
        eo[1]=0.0;
    }
    return RT_SUCCESS;
}
 
/**************************************************************************
 * NAME: getrad()
 *
 * DESCRIPTION: Compute scattering matrices subroutine
 *
 *************************************************************************/
 
int Afrt::getrad()
{
// get the radiances at the top and bottom of the atmosphere
    if (jpass == 0) {
        fio = unit53[nolyr];
    } else {
        fio = iwrt->at(nolyr);
    }
    double suma = 0.0;
    for (int i=0; i<nx-1; i++ ) {
        double psum = 0.0;
        for (int k=0; k<2; k++ ) {
            psum = fio.f[0][i][k] + fio.f[jpart-1][i][k] + psum;
            for (int j=1; j<nophi/2; j++ ) {
                psum = psum + 2.0*fio.f[j][i][k];
            }
        }
        suma = suma + psum*dcmu2[i]*6.2831854/double(nophi);
    }
    double sumdwn = suma;
    factr = eo[0] + eo[1];
    fluxd[nolyr][jpass] = sumdwn;
    for (int i=0; i<nx-1; i++ ) {
        for (int j=0; j<jpart; j++) {
            geoidn[j][i][jpass] = fio.f[j][i][0]+fio.f[j][i][1];
        }
    }
    if (jpass == 0) {
        fio = unit53[0];
    } else {
        fio = iwrt->at(0);
    }
    for (int i=nx-1; i<nmu; i++) {
        int ij = i-(nx-1);
        for (int j=0; j<jpart; j++) {
            geoiup[j][ij][jpass] = fio.f[j][i][0]+fio.f[j][i][1];
        }
    }
    suma = 0.0;
    for (int i=nx-1; i<nmu; i++) {
        double psum = 0.0;
        for (int k=0; k<2; k++) {
            psum = fio.f[0][i][k]+fio.f[jpart-1][i][k] + psum;
        }
        for (int k=0; k<2; k++) {
            for (int j=1; j<nophi/2; j++) {
                psum = psum + 2.*fio.f[j][i][k];
            }
        }
        suma=suma+psum*dcmu2[i]*6.2831854/nophi;
    }
    double sumup = abs(suma);
    fluxu[0][jpass]=sumup;
 
    return RT_SUCCESS;
}
 
/**************************************************************************
 * NAME: multp1d()
 *
 * DESCRIPTION: Compute scattering matrices subroutine
 *
 *************************************************************************/
 
int Afrt::multp1d()
{
// compute the downward diffuse radiation at each level
    sngla = fnull;
    snglb = fnull;
    iwrt->at(0) = irad->at(0);
    for (int il=0; il<nolyr; il++) {
        if(il==0) {
            fiib = iwrt->at(il);
            fio = irad->at(il+1);
            sngla = unit53[il];
            snglb = unit53[il+1];
// determine average intensity at the center of the layer il
            for (int i=0; i<nmu; i++) {
                for (int j=0; j<jpart; j++) {
                    for (int k=0; k<nstk; k++) {
                        fiic.f[j][i][k] = 0.50 * (fio.f[j][i][k] + fiib.f[j][i][k]);
                        ftmp.f[j][i][k] = 0.0;
                    }
                }
            }
            tmsl=dtmm[il]*qsqt*cnst;
            trsl=dtrr[il]*conr;
            dlyr=dtot[il];
            if(ifc==0) tmsl=0.0;
            mdiffn(0,nx-2,il,dlyr);
            for (int i=0; i<nx-1; i++) {
                for (int j=0; j<jpart; j++) {
                    for (int k=0; k<nstk; k++) {
                        fio.f[j][i][k]=fiib.f[j][i][k]*emdtm[i][il] +
                                ftmp.f[j][i][k]*(1.0-emdtm[i][il]) +
                                snglb.f[j][i][k]-sngla.f[j][i][k]*emdtm[i][il];
                    }
                }
            }
            iwrt->at(il+1) = fio;
        } else {
            fiib = iwrt->at(il-1);
            sngla = unit53[il-1];
            snglb = unit53[il+1];
            fiic = fio;
            ftmp = fnull;
            fio = irad->at(il+1);
            tmsl=(dtmm[il-1]+dtmm[il])*qsqt*cnst;
            trsl=(dtrr[il-1]+dtrr[il])*conr;
            dlyr=(dtot[il-1]+dtot[il]);
            if(ifc==0) tmsl=0.0;
            mdiffn(0,nx-2,il,dlyr);
            for (int i=0; i<nx-1; i++) {
                for (int j=0; j<jpart; j++) {
                    for (int k=0; k<nstk; k++) {
                        fio.f[j][i][k] = fiib.f[j][i][k]*emdtm[i][il-1]*emdtm[i][il]+
                           ftmp.f[j][i][k]*(1.0-emdtm[i][il-1]*emdtm[i][il])+
                           snglb.f[j][i][k]-sngla.f[j][i][k]*emdtm[i][il-1]*emdtm[i][il];
                    }
                }
            }
            iwrt->at(il+1) = fio;
        }
    }
    if(iref==1) {
//        fltocn_new();
        iwrt->at(nolyr) = fio;
    }
 
    if(iref==2) {
        rsea_new();
        iwrt->at(nolyr) = fio;
    }
 
    if(iref==3) {
//        brdfg(fio,brdfx,cosmu,dcmu,dcmu2,eo,ddphi,amuo,M_PI,sumc,
//                sumcpi,sumdwn,calb,kkx,nx-1,nmu,jpart);
        iwrt->at(nolyr) = fio;
    }
// compute the upward diffuse radiation at each level
    for (int il=0; il<nolyr; il++) {
        int im=nolyr-il-1;
        if(il==0) {
            fiib = iwrt->at(im+1);
            fio = iwrt->at(im);
            sngla = unit53[im+1];
            snglb = unit53[im];
            for (int i=0; i<nmu; i++) {
                for (int j=0; j<jpart; j++) {
                    for (int k=0; k<nstk; k++) {
                        fiic.f[j][i][k] = 0.50*(fio.f[j][i][k]+fiib.f[j][i][k]);
                        ftmp.f[j][i][k] = 0.0;
                    }
                }
            }
            tmsl=dtmm[im]*qsqt*cnst;
            trsl=dtrr[im]*conr;
            dlyr=dtot[im];
            if(ifc==0) tmsl=0.0;
            mdiffn(nx-1,nmu-1,im,dlyr);
            for (int i=nx-1; i<nmu; i++) {
                for (int j=0; j<jpart; j++) {
                    for (int k=0; k<nstk; k++) {
                        fio.f[j][i][k] = fiib.f[j][i][k]*emdtm[i][im] +
                                ftmp.f[j][i][k]*(1.0-emdtm[i][im]) +
                                snglb.f[j][i][k]-sngla.f[j][i][k]*emdtm[i][im];
                    }
                }
            }
            iwrt->at(im) = fio;
        } else {
            fiib = iwrt->at(im+2);
            sngla = unit53[im+2];
            snglb = unit53[im];
            for (int i=0; i<nmu; i++) {
                for (int j=0; j<jpart; j++) {
                    for (int k=0; k<nstk; k++) {
                        fiic.f[j][i][k]=fio.f[j][i][k];
                        ftmp.f[j][i][k]=0.0;
                    }
                }
            }
            fio = iwrt->at(im);
            tmsl=(dtmm[im+1]+dtmm[im])*qsqt*cnst;
            trsl=(dtrr[im+1]+dtrr[im])*conr;
            dlyr=(dtot[im+1]+dtot[im]);
            if(ifc==0) tmsl=0.0;
            mdiffn(nx-1,nmu-1,im,dlyr);
            for (int i=nx-1; i<nmu; i++) {
                for (int j=0; j<jpart; j++) {
                    for (int k=0; k<nstk; k++) {
                        fio.f[j][i][k] = fiib.f[j][i][k]*emdtm[i][im+1]*emdtm[i][im] +
                                ftmp.f[j][i][k]*(1.0-emdtm[i][im+1]*emdtm[i][im])+
                                snglb.f[j][i][k]-sngla.f[j][i][k]*emdtm[i][im+1]*emdtm[i][im];
                    }
                }
            }
            iwrt->at(im) = fio;
        }
    }
    return RT_SUCCESS;
}
 
/**************************************************************************
 * NAME: multp2d()
 *
 * DESCRIPTION: Compute scattering matrices subroutine
 *
 *************************************************************************/
 
int Afrt::multp2d()
{
    // compute the downward diffuse radiation at each level c
    iwrt->at(0) = irad->at(0);
    for (int il=0; il<nolyr; il++) {
        if (il == 1) {
            fiib = iwrt->at(il);
            fio = irad->at(il+1);
// determine average intensity at the center of the layer il
            for (int i=0; i<nmu; i++) {
                for (int j=0; j<jpart; j++) {
                    for (int k=0; k<nstk; k++) {
                        fiic.f[j][i][k] = 0.50 * (fio.f[j][i][k] + fiib.f[j][i][k]);
                        ftmp.f[j][i][k] = 0.0;
                    }
                }
            }
            tmsl=dtmm[il]*qsqt*cnst;
            trsl=dtrr[il]*conr;
            dlyr=dtot[il];
            if(ifc==0) tmsl=0.0;
            mdiffn(0,nx-2,il,dlyr);
            for (int i=0; i<nx-1; i++) {
                for (int j=0; j<jpart; j++) {
                    for (int k=0; k<nstk; k++) {
                        fio.f[j][i][k] = fiib.f[j][i][k]*emdtm[i][il] +
                                ftmp.f[j][i][k]*(1.0-emdtm[i][il]);
                    }
                }
            }
            iwrt->at(il+1) = fio;
        } else {
            fiib = iwrt->at(il-1);
            fiic = fio;
            ftmp = fnull;
            fio = irad->at(il+1);
            tmsl=(dtmm[il-1]+dtmm[il])*qsqt*cnst;
            trsl=(dtrr[il-1]+dtrr[il])*conr;
            dlyr=(dtot[il-1]+dtot[il]);
            if(ifc==0) tmsl=0.0;
            mdiffn(0,nx-2,il,dlyr);
            for (int i=0; i<nx-1; i++) {
                for (int j=0; j<jpart; j++) {
                    for (int k=0; k<nstk; k++) {
                        fio.f[j][i][k] = fiib.f[j][i][k]*emdtm[i][il-1]*emdtm[i][il] +
                                ftmp.f[j][i][k]*(1.0-emdtm[i][il-1]*emdtm[i][il]);
                    }
                }
            }
            iwrt->at(il+1) = fio;
        }
    }
// compute the upward diffuse radiation at each level
    for (int il=0; il<nolyr; il++) {
        int im=nolyr-il-1;
        if(il==1) {
            fiib = iwrt->at(im+1);
            fio = iwrt->at(im);
            for (int i=0; i<nmu; i++) {
                for (int j=0; j<jpart; j++) {
                    for (int k=0; k<nstk; k++) {
                        fiic.f[j][i][k]=0.50*(fio.f[j][i][k]+fiib.f[j][i][k]);
                        ftmp.f[j][i][k]=0.0;
                    }
                }
            }
            tmsl=dtmm[im]*qsqt*cnst;
            trsl=dtrr[im]*conr;
            dlyr=dtot[im];
            if(ifc==0) tmsl=0.0;
            mdiffn(nx-1,nmu-1,im,dlyr);
            for (int i=nx-1; i<nmu; i++) {
                for (int j=0; j<jpart; j++) {
                    for (int k=0; k<nstk; k++) {
                        fio.f[j][i][k] = fiib.f[j][i][k]*emdtm[i][im] +
                                ftmp.f[j][i][k]*(1.0-emdtm[i][im]);
                    }
                }
            }
            iwrt->at(im) = fio;
        } else {
            fiib = iwrt->at(im+2);
            for (int i=0; i<nmu; i++) {
                for (int j=0; j<jpart; j++) {
                    for (int k=0; k<nstk; k++) {
                        fiic.f[j][i][k]=fio.f[j][i][k];
                        ftmp.f[j][i][k]=0.0;
                    }
                }
            }
            fio = iwrt->at(im);
            tmsl=(dtmm[im+1]+dtmm[im])*qsqt*cnst;
            trsl=(dtrr[im+1]+dtrr[im])*conr;
            dlyr=(dtot[im+1]+dtot[im]);
            if(ifc==0) tmsl=0.0;
            mdiffn(nx-1,nmu-1,im,dlyr);
            for (int i=nx-1; i<nmu; i++) {
                for (int j=0; j<jpart; j++) {
                    for (int k=0; k<nstk; k++) {
                        fio.f[j][i][k] = fiib.f[j][i][k]*emdtm[i][im+1]*emdtm[i][im] +
                                ftmp.f[j][i][k]*(1.0-emdtm[i][im+1]*emdtm[i][im]);
                    }
                }
            }
            iwrt->at(im) = fio;
        }
    }
    return RT_SUCCESS;
}
 
/**************************************************************************
 * NAME: geocor()
 *
 * DESCRIPTION: Compute scattering matrices subroutine
 *
 *************************************************************************/
int Afrt::geocor()
{
 
    int m1 = jpass-2;
    int m2 = jpass-1;
    int m3 = jpass;
    double tup[nphi][nthe][nsza];
    double tdn[nphi][nthe][nsza];
// apply geo series correction to the fluxes leaving the bottom of the atmosphere
    double temp1[nlyr];
    double temp2[nlyr];
    double rr1, rr2;
    geom(fluxd[nolyr][m1],fluxd[nolyr][m2],fluxd[nolyr][m3],rr1,temp1[nolyr]);
    geom(fluxu[0][m1],fluxu[0][m2],fluxu[0][m3],rr2,temp2[0]);
    //
    // if(kzz==0) {
    //     double totflxgs=temp1[nolyr]+temp2[0]+amuo*factr;
    //     double ctest=totflxgs/amuo;
    // }
    double gz, sb;
    if(kzz==0) gz = temp1[nolyr];
    if(kzz==1 && itrans==0) sb = temp1[nolyr];
// apply geo series correction to the radiances leaving the top &
// bottom of the atmosphere
    ftmp = unit71[m3];
    ftmpa = unit71[m2];
    ftmpb = unit71[m1];
// if iglint=1, remove the direct component
    if(kzz==0 && iglint==1) {
        for (int i=nx-1; i<nmu; i++) {
            for (int j=0; j<jpart; j++) {
                for (int k=0; k<nstk; k++) {
                    double abcz=exp(-tautot/abs(cosmu[i]));
                    ftmp.f[j][i][k]=ftmp.f[j][i][k]-fglint.f[j][i][k]*abcz;
                    ftmpa.f[j][i][k]=ftmpa.f[j][i][k]-fglint.f[j][i][k]*abcz;
                    ftmpb.f[j][i][k]=ftmpb.f[j][i][k]-fglint.f[j][i][k]*abcz;
                }
            }
        }
    }
    //
    double ratiog;
    for (int i=nx-1; i<nmu; i++) {
        int m = nmu-i-1;
        for (int j=0; j<jpart; j++) {
            for (int k=0; k<nstk; k++) {
                if(ftmp.f[j][i][k] <= 1.0e-15) {
                    fioup.f[j][m][k]=ftmp.f[j][i][k];
                } else {
                    geom(ftmpb.f[j][i][k],ftmpa.f[j][i][k],ftmp.f[j][i][k],
                            ratiog,fioup.f[j][m][k]);
                }
            }
        }
    }
    ftmp = unit72[m3];
    ftmpa = unit72[m2];
    ftmpb = unit72[m1];
    for (int i=nx-1; i<nmu; i++) {
        int m = nmu-i-1;
        for (int j=0; j<jpart; j++) {
            for (int k=0; k<nstk; k++) {
                if(ftmp.f[j][i][k]<=1.0e-15) {
                    fioup_btm.f[j][m][k]=ftmp.f[j][i][k];
                } else {
                    geom(ftmpb.f[j][i][k],ftmpa.f[j][i][k],ftmp.f[j][i][k],
                            ratiog,fioup_btm.f[j][m][k]);
                }
            }
        }
    }
    for (int i=0; i<nx-1; i++) {
        for (int j=0; j<jpart; j++) {
            for (int k=0; k<nstk; k++) {
                if(ftmp.f[j][i][k]<=1.0e-15) {
                    fiodn.f[j][i][k]=ftmp.f[j][i][k];
                } else {
                    geom(ftmpb.f[j][i][k],ftmpa.f[j][i][k],ftmp.f[j][i][k],
                            ratiog,fiodn.f[j][i][k]);
                }
            }
        }
    }
    if(kzz==0) {
        fdirc[ksza-1]=amuo*factr*M_PI;
        fdown[ksza-1]=temp1[nolyr]*M_PI;
        fup[ksza-1]=temp2[0]*M_PI;
        for (int i=0; i<nx-1; i++) {
            for (int j=0; j<jpart; j++) {
                xzeroz[j][i][ksza-1][0]=fioup.f[j][i][0]+fioup.f[j][i][1];
                xzeroz[j][i][ksza-1][1]=fioup.f[j][i][0]-fioup.f[j][i][1];
                xzeroz[j][i][ksza-1][2]=fioup.f[j][i][2];
                xzeroz[j][i][ksza-1][3]=fioup.f[j][i][3];
                xzerod[j][i][ksza-1][0]=fiodn.f[j][i][0]+fiodn.f[j][i][1];
                xzerod[j][i][ksza-1][1]=fiodn.f[j][i][0]-fiodn.f[j][i][1];
                xzerod[j][i][ksza-1][2]=fiodn.f[j][i][2];
                xzerod[j][i][ksza-1][3]=fiodn.f[j][i][3];
            }
        }
        if(iref==1 || iref==2 || iref==3) {
            oalb[ksza-1]=calb;
        }
    }
    for (int i=0; i<nx-1; i++) {
        for (int j=0; j<jpart; j++) {
            xzero_up[j][i][ksza-1]=fioup.f[j][i][0]+fioup.f[j][i][2];
            xzero_btm[j][i][ksza-1]=fioup_btm.f[j][i][0]+fioup_btm.f[j][i][2];
        }
    }
    if(kzz==1 && itrans==1) {
        for (int i=0; i<nx-1; i++) {
            for (int j=0; j<jpart; j++) {
                xzero_up[j][i][ksza-1]=fioup.f[j][i][0]+fioup.f[j][i][2];
                xzero_btm[j][i][ksza-1]=fioup_btm.f[j][i][0]+fioup_btm.f[j][i][2];
            }
        }
    } else if (kzz==1 && itrans==0) {
        double ef=amuo*efactb[nolyr];
        double ftot=(gz+ef);
        for (int i=0; i<nx-1; i++) {
            for (int j=0; j<jpart; j++) {
                tdn[j][i][0]=(fiodn.f[j][i][0]+fiodn.f[j][i][1]);
                tdn[j][i][1]=(fiodn.f[j][i][0]-fiodn.f[j][i][1]);
                tdn[j][i][2]=(fiodn.f[j][i][2]);
                tdn[j][i][3]=(fiodn.f[j][i][3]);
                tup[j][i][0]=(fioup.f[j][i][0]+fioup.f[j][i][1]);
                tup[j][i][1]=(fioup.f[j][i][0]-fioup.f[j][i][1]);
                tup[j][i][2]=(fioup.f[j][i][2]);
                tup[j][i][3]=(fioup.f[j][i][3]);
            }
        }
        sbarz[ksza-1]=sb;
        for (int i=0; i<nx-1; i++) {
            for (int j=0; j<jpart; j++) {
                for (int s=0; s<nstk; s++) {
                    tupz[j][i][ksza-1][s]=tup[j][i][s]*ftot;
                    tdwnz[j][i][ksza-1][s]=tdn[j][i][s]*ftot;
                }
            }
        }
    }
    if(insz>1) {
        double ef=amuo*efactb[nolyr];
        double ftot=(gz+ef);
        sbarz[ksza-1]=sb;
        for (int i=0; i<nx-1; i++) {
            for (int j=0; j<jpart; j++) {
                for (int s=0; s<nstk; s++) {
                    tupz[j][i][ksza-1][s]=tup[j][i][s]*ftot;
                    tdwnz[j][i][ksza-1][s]=tdn[j][i][s]*ftot;
                }
            }
        }
    }
    return RT_SUCCESS;
}
 
/**************************************************************************
 * NAME: rsea_new()
 *
 * DESCRIPTION: Treat the lower boundary of the atmosphere as non-lambertian
 *   surface (rough ocean surface) and reflect all incident light
 *   (diffuse and direct light) using the reflection matrix of the surface
 *
 *************************************************************************/
int Afrt::rsea_new()
{
    double tempa[nphi][2*nsza][nstk];
    double tempc[nphi][2*nsza][nstk];
    double tss[nphi][2*nsza][nstk];
    for (int i=0; i<nmu; i++) {
        for (int j=0; j<jpart; j++) {
            for (int k=0; k<nstk; k++) {
                tempa[j][i][k]=0.0;
                tempc[j][i][k]=0.0;
                tss[j][i][k]=0.0;
            }
        }
    }
// total downward flux
    double sumdnz;
    fluxlvl(fio.f,sumdnz,0);
    double pflx0 = (eo[0]+eo[1])*cosmu[kkx];
    double pflx = pflx0+sumdnz;
// contribution from the direct beam
    for (int it=0; it<(nx-1); it++) {
        int itt=nmu-it-1;
        for (int ip=0; ip<jpart; ip++) {
            for (int is=0; is<nstk; is++) {
                for (int k=0; k<nstk; k++) {
                    int ms=is*4+k;
                    tempa[ip][itt][is]=tempa[ip][itt][is] +
                            txx[it][kkx][ip][ms]*eo[k];
                }
            }
        }
    }
    double flxtempa;
    fluxlvl(tempa,flxtempa,1);
// contribution from diffuse beam (integrate over thetaprime & phiprime)
    for (int it=0; it<(nx-1); it++) {
        int itt=nmu-it-1;
        for (int ip=0; ip<jpart; ip++) {
            for (int is=0; is<nstk; is++) {
                for (int itp=0; itp<(nx-1); itp++) {
                    double sangi = abs(dcmu[itp])*ddphi;
                    for (int ipp=0; ipp<nophi; ipp++) {
                        int jpp=ipp;
                        int jpx=abs(ipp-ip);
                        double fiit[nphi][2*nsza][nstk];
                        if(ipp>jpart-1)jpp=nophi-jpp;
                        if(jpx>jpart-1)jpx=nophi-jpx;
                        for (int i=0; i<2; i++) {
                            fiit[jpp][itp][i]=fio.f[jpp][itp][i];
                            fiit[jpp][itp][i+2]=fio.f[jpp][itp][i+2];
                            if(ipp>jpart-1) {
                                fiit[jpp][itp][i+2]= -fiit[jpp][itp][i+2];
                            }
                        }
                        double prodt3a=0.0;
                        for (int k=0; k<nstk; k++) {
                            int m = is*4+k;
                            prodt3a=prodt3a+txx[it][itp][jpx][m]*fiit[jpp][itp][k];
                        }
                        tempc[ip][itt][is]=tempc[ip][itt][is]+prodt3a*sangi;
                    }
                }
            }
        }
    }
// save the glint contribution from the direct beam
// into fglint buffer
    for (int i=0; i<nmu; i++) {
        for (int j=0; j<jpart; j++) {
            for (int k=0; k<nstk; k++) {
                fglint.f[j][i][k]=tempa[j][i][k];
            }
        }
    }
    double flxtempc;
    fluxlvl(tempc,flxtempc,1);
    for (int i=nx-1; i<nmu; i++) {
        int ii=i-(nx-1);
        //double sang = abs(dcmu2[i])*ddphi;
        for (int j=0; j<jpart; j++) {
            for (int k=0; k<nstk; k++) {
                tss[j][i][k]=tempa[j][i][k]+tempc[j][i][k];
                tss[j][ii][k]=0.0;
            }
        }
    }
 
// compute the upwelling fluxes due to each component
    double flxza, flxzc;
    fluxlvl(tempa,flxza,1);
    fluxlvl(tempc,flxzc,1);
    //double albdtot=(flxza+flxzc)/pflx;
    albrdr[ksza-1]=flxza/pflx;
    albrdf[ksza-1]=flxzc/pflx;
 
// correct the upwelling radiance for foam and underwater radiances
// foam correction effective reflectance of the whitecaps (Koepke, 1984)
    double ref[39] = {
            0.220,0.220,0.220,0.220,0.220,0.220,0.215,0.210,0.200,0.190,
            0.175,0.155,0.130,0.080,0.100,0.105,0.100,0.080,0.045,0.055,
            0.065,0.060,0.055,0.040,0.000,0.000,0.000,0.000,0.000,0.000,
            0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000};
    double radxi[nphi][2*nsza];
    memset(radxi, 0.0, nphi*2*nsza*sizeof(double));
    double wl=wvlth*1.0e4;
// compute whitecaps reflectance 'afoam' for wavelength 'wl'
    int iwl=(int)((wl-0.2)/0.1);
    double wlp=0.5+(iwl-1)*0.1;
    double afoam=ref[iwl+1]+(wl-wlp)/0.1*(ref[iwl]-ref[iwl+1]);
    double fracfm, ffoam, xifm;
    if(ifoam==1) {
        fracfm=2.95e-6*pow(vz,3.52);
        ffoam=afoam*fracfm;
        xifm=(ffoam*pflx)/M_PI;
    } else {
        ffoam=0.0;
        fracfm=0.0;
        xifm=0.0;
    }
// underwater contribution compute the upwelling radiance from the
// albedo of water just below the ocean surface
    if(iwatr==1) {
// compute the downward diffused flux
        double fzz;
        fluxlvl(fio.f,fzz,0);
        double suma = 0.0;
        double sumtt=0.0;
        double sumd=0.0;
        double sumu=0.0;
        for (int i=0; i<(nx-1); i++) {
            double psum = 0.0;
            for (int k=0; k<2; k++) {
                psum=(fio.f[0][i][k]+fio.f[jpart][i][k])+psum;
                for (int j=1; j<nophi/2; j++) {
                    psum=psum+2.*fio.f[j][i][k];
                }
            }
            suma=suma+psum*dcmu2[i]*(1.0-rfair[i])*ddphi;
        }
        sumtt=(eo[0]+eo[1])*cosmu[kkx]*(1.0-rfair[kkx]);
        sumd=suma+sumtt;
// compute the upwelling flux and radiance just below the ocean surface
        sumu=alw*sumd;
        double radu=sumu/M_PI;
        albtdr[ksza-1]=sumtt;
        albtdf[ksza-1]=suma;
        albtdt[ksza-1]=sumd;
        albtrf[ksza-1]=sumu;
// compute the radiance just above the ocean surface
        for (int i=0; i<(nx-1); i++) {
            int ii=nmu-i-1;
            for (int j=0; j<jpart; j++) {
                radxi[j][ii]=((1.0-rfwat[i])*radu)/(pow(xrw,2.0)+xiw*2);
                radxi[j][i]=0.0;
            }
        }
    }
// compute the flux of the water leaving radiance
    double qsumzz=0.0;
    for (int i=nx-1; i<nmu; i++) {
        double qsumz=0.0;
        for (int j=0; j<jpart; j++) {
            if(j==0 || j==jpart-1) {
                qsumz=qsumz+radxi[j][i];
            } else {
                qsumz=qsumz+2.0*radxi[j][i];
            }
        }
        qsumzz=qsumzz+qsumz*abs(dcmu2[i])*ddphi;
    }
    albwl[ksza-1]=qsumzz;
    for (int it=nx-1; it<nmu; it++) {
        for (int ip=0; ip<jpart; ip++) {
            for (int ks=0; ks<nstk; ks++) {
                if(ks<=2) {
                    fio.f[ip][it][ks]=(1.0-fracfm)*tss[ip][it][ks] +
                                       0.5*xifm+0.5*radxi[ip][it];
                } else {
                    fio.f[ip][it][ks]=(1.0-fracfm)*tss[ip][it][ks];
                }
            }
        }
    }
    //
    for (int it=nx-1; it<nmu; it++) {
        int itp=nmu-it-1;
        for (int ip=0; ip<jpart; ip++) {
            raddir[ip][itp][ksza-1]=tempa[ip][it][0]+tempa[ip][it][0];
            radocn[ip][itp][ksza-1]=radxi[ip][it];
            radsky[ip][itp][ksza-1]=tempc[ip][it][0]+tempc[ip][it][0];
        }
    }
    double sumc;
    fluxlvl(fio.f,sumc,1);
    //double sumcpi=sumc*M_PI;
    calb=sumc/pflx;
 
    return RT_SUCCESS;
}
 
/**************************************************************************
 * NAME: fluxlvl()
 *
 * DESCRIPTION: Compute fluxes
 *
 *************************************************************************/
int Afrt::fluxlvl(double buft[][2*nsza][nstk], double& sumg, int iflag)
{
    int m = 0;
    int n = nx-1;
    if(iflag==1) {
        m=nx-1;
        n=nmu;
    }
    double anophi = (double) nophi;
    sumg = 0.0;
    for (int i=m; i<n; i++) {
        double psumg = 0.0;
        for (int k=0; k<2; k++) {
            psumg=buft[0][i][k]+buft[jpart-1][i][k]+psumg;
            for (int j=1; j<nophi/2; j++) {
                psumg=psumg+2.0*buft[j][i][k];
            }
        }
        sumg=sumg+psumg*dcmu2[i]*6.2831854/anophi;
    }
    sumg=abs(sumg);
 
    return RT_SUCCESS;
}
 
/**************************************************************************
 * NAME: geom()
 *
 * DESCRIPTION: Compute scattering matrices subroutine
 *
 *************************************************************************/
int Afrt::geom(double a, double b, double c, double& r, double& g)
{
    g = c;
    r = 0.0;
    if (g > 1.0e-15) {
        double dif1 = b - a;
        double dif2 = c - b;
        double dif2p = abs(100.0 * dif2 / c);
        if (dif2p >= 1.0e-2) {
            r = dif2 / dif1;
            g = a + dif1 / (1.0 - r);
        }
    }
    return RT_SUCCESS;
}
 
/**************************************************************************
 * NAME: xntpln()
 *
 * DESCRIPTION: Compute scattering matrices subroutine
 *
 *************************************************************************/
int Afrt::xntpln(double x,double x1,double x2,double y1,double y2,double &y)
{
    double slope=(y1-y2)/(x1-x2);
    y = y1 + slope*(x-x1);
 
    return RT_SUCCESS;
}
 
/**************************************************************************
 * NAME: mdiffn()
 *
 * DESCRIPTION: Compute scattering matrices subroutine
 *
 *************************************************************************/
int Afrt::mdiffn(int ib, int ie, int il, double dlyr)
{
    double fiit[nphi][2*nsza][nstk];
    double fiiz[nphi][nthe][2*nsza][nstk];
 
    for (int kk=0; kk<nmu; kk++) {
        for (int ll=0; ll<jpart; ll++) {
            for (int is=0; is<nstk; is++) {
                fiit[ll][kk][is]=fiic.f[ll][kk][is];
            }
        }
        for (int ll=jpart; ll<=nophi; ll++) {
            fiit[ll][kk][0]=fiic.f[nophi-ll][kk][0];
            fiit[ll][kk][1]=fiic.f[nophi-ll][kk][1];
            fiit[ll][kk][2]=-fiic.f[nophi-ll][kk][2];
            fiit[ll][kk][3]=-fiic.f[nophi-ll][kk][3];
        }
    }
    for (int kk=0; kk<nmu; kk++) {
        for (int ip=0; ip<jpart; ip++) {
            for (int ll=0; ll<=nophi; ll++) {
                int mmp=ll+ip;
                if(ll>jpart-1) {
                    mmp=nophi*(1+mmp/nophi)-mmp;
                }
                for (int is=0; is<nstk; is++) {
                    fiiz[ll][ip][kk][is]=fiit[mmp][kk][is];
                }
            }
        }
    }
    for (int ii=0; ii<nmu; ii++) {
        for (int kk=0; kk<nmu; kk++) {
            for (int ll=jpart; ll<=nophi; ll++) {
                for (int is=0; is<nstk; is++) {
                    for (int j=0; j<nstk; j++) {
                      int ij=is*4+j;
                      if (is<2) {
                          ppin[kk][ii][ll][ij]=ppin[kk][ii][nophi-ll][ij];
                          ppin[kk][ii][ll][ij+16]=ppin[kk][ii][nophi-ll][ij+16];
                      } else {
                          ppin[kk][ii][ll][ij]=-ppin[kk][ii][nophi-ll][ij];
                          ppin[kk][ii][ll][ij+16]=-ppin[kk][ii][nophi-ll][ij+16];
                      }
                   }
               }
            }
        }
    }
 
    for (int it=ib; it<=ie; it++) {
        for (int ip=0; ip<jpart; ip++) {
            for (int is=0; is<nstk; is++) {
                double sumta=0.0;
                for (int kk=0; kk<nmu; kk++) {
                    double sumtb=0.0;
                    for (int ll=0; ll<nophi; ll++) {
                        double prod1=0.0;
                        for (int j=0; j<nstk; j++) {
                            int ij=is*4+j;
                            prod1=prod1+(c[kk]*tmsl*ppin[kk][it][ll][ij]+
                              trsl*ppin[kk][it][ll][ij+16])*fiiz[ll][ip][kk][j]/dlyr;
                        }
                        sumtb=sumtb+prod1;
                     }
                     sumta=sumta+sumtb*dcmu[kk];
                 }
                 ftmp.f[ip][it][is]=sumta*ddphi;
            }
        }
    }
    return RT_SUCCESS;
}