subroutine drybicki(tmin,dtau,u,s,ds,ntau,xu,wu,nu,opp,dopp,wpp & ,npp,b1,oppsum,f) implicit double precision (o-z,a-h) include "parameters.dat" c --------------------------------------------------------------- c Input arrays and variables c --------------------------------------------------------------- dimension dtau(nd1),u(nd1,nd2),s(nd1),ds(nd1,nd2),xu(nd6),wu(nd6) & , opp(nd1,nd2), dopp(nd1,nd2), wpp(nd2), b1(nd2) ,oppsum(nd1) c --------------------------------------------------------------- c Output arrays c --------------------------------------------------------------- dimension f(nd1,nd1) c --------------------------------------------------------------- c Internal arrays and variables c --------------------------------------------------------------- dimension gf(nd1), dgf(nd1), dt(nd1), w(nd1), ta(nd1) & , tb(nd1), tc(nd1), tu(nd1,nd1), e(nd1,nd1), fj(nd1) & , va(nd1), vb(nd1), vc(nd1), v(nd1,nd1), indx(nd1) c --------------------------------------------------------------- c Initialize arrays c --------------------------------------------------------------- do j1 = 1, ntau do j2 = 1, ntau e(j2,j1) = zero f(j2,j1) = zero enddo e(j1,j1) = one end do c --------------------------------------------------------------- c Loop through frequencies c --------------------------------------------------------------- do jf = 1, npp bf1 = b1(jf) wgf = wpp(jf) do jt = 1, ntau gf(jt) = opp(jt,jf) dgf(jt) = dopp(jt,jf) end do c ---------------------------------------------- c Calculate monochromatic optical depth c ---------------------------------------------- do jt = 1, ntau - 1 dt(jt) = half*(gf(jt)+gf(jt+1))*dtau(jt) c dt(jt) = dmax1(half*(gf(jt)+gf(jt+1))*dtau(jt),tmin) end do c ---------------------------------------------- c Find optically thin level c ---------------------------------------------- tmax = zero do k = 1, ntau -1 tmax = tmax + dt(k) end do k = 1 tot = dt(k) do while((tot .lt. tmin).and.(k.lt.ntau)) k = k + 1 tot = tot + dt(k) end do ntmn = max(k,1) if ((tmax .lt. tmin).or.(ntmn .eq. ntau)) then fac = half*wgf*bf1 do jt = 1, ntau f(jt,jt)=f(jt,jt)+fac*dgf(jt) end do else c ------------------------------------------------ c Loop through zenith angles c ------------------------------------------------ do ju = 1, nu c ------------------------------------- c Calculate W matrix (Diagonal) c ------------------------------------- wght=wu(ju)*wgf do jt = 1, ntau w(jt)=wght*gf(jt) end do c ------------------------------------- c Calculate T matrix (Tridiagonal) c ------------------------------------- do jt = 1, ntmn-1 ta(jt) = zero c tb(jt) = one c tc(jt) = -one tb(jt) = one/tmin/tmin tc(jt) = -one/tmin/tmin end do call tmatrix(xu(ju), dt(ntmn), ntau-ntmn+1 & , ta(ntmn), tb(ntmn), tc(ntmn)) c ------------------------------------- c Calculate LU decomposition of T c ------------------------------------- c call ludcmpt(ta,tb,tc,ntau) c ------------------------------------- c *** Calculate E Matrix *** c ------------------------------------- c Calculate U matrix (diagonal) c ------------------------------------- do j1 = 1, ntau do j2 = 1, ntau tu(j2,j1)=zero end do end do do jt = ntmn, ntau tu(jt,jt) = u(jt,jf)/oppsum(jt) end do c ------------------------------------- c Calculate T(-1)*U store in TU c ------------------------------------- c do jt = 1, ntau c call lubksbt(ta,tb,tc,tu(1,jt),ntau) c enddo do jt = 1, ntau call tridag(ta,tb,tc,tu(1,jt),ntau,ierr) if(ierr .eq. 1) then write(unit=*,fmt=*) ' Error in drybicki' do k = 1, ntau write(unit=*,fmt=*) ta(k),tb(k),tc(k) end do stop end if end do c ------------------------------------- c Calculate W*T(-1)*U add into E c ------------------------------------- do j1 = 1, ntau do j2 = 1, ntau e(j2,j1)=e(j2,j1)-w(j2)*tu(j2,j1) end do end do C ------------------------------------- C *** Calculate F Matrix *** C ------------------------------------- C Calculate Feautrier Intensities C ------------------------------------- do jt = 1, ntmn-1 fj(jt) = zero end do do jt = ntmn, ntau-1 fj(jt) = s(jt) end do fj(ntau) = s(ntau) + (two*xu(ju)/dt(ntau-1))*bf1 c call lubksbt(ta,tb,tc,fj,ntau) call tridag(ta,tb,tc,fj,ntau,ierr) if(ierr .eq. 1) then write(unit=*,fmt=*) ' Error in drybicki' do k = 1, ntau write(unit=*,fmt=*) ta(k),tb(k),tc(k) end do stop end if c ------------------------------------- c Add derivative of W terms into F c ------------------------------------- do jt = 1, ntau f(jt,jt)=f(jt,jt)+wght*fj(jt)*dgf(jt) end do c ------------------------------------- c Calculate V matrix (Tridiagonal) c ------------------------------------- fup = bf1 fdn = zero call vmatrix(xu(ju),dtau(ntmn),gf(ntmn),dgf(ntmn) & ,fj(ntmn),ntau-ntmn+1,fup,fdn,va(ntmn),vb(ntmn) & ,vc(ntmn)) do j1 = 1, ntau do j2 = 1, ntau v(j2,j1) = zero end do end do v(ntmn,ntmn) = ds(ntmn,jf)-vb(ntmn) v(ntmn,ntmn+1) = -vc(ntmn) do j1 = ntmn+1, ntau-1 v(j1,j1-1) = -va(j1) v(j1,j1) = ds(j1,jf)-vb(j1) v(j1,j1+1) = -vc(j1) end do v(ntau,ntau-1) = -va(ntau) v(ntau,ntau) = ds(ntau,jf)-vb(ntau) c ------------------------------------- c Calculate T(-1)*V store in V c ------------------------------------- c do jt = 1, ntau c call lubksbt(ta,tb,tc,v(1,jt),ntau) c end do do jt = 1, ntau call tridag(ta,tb,tc,v(1,jt),ntau,ierr) if(ierr .eq. 1) then write(unit=*,fmt=*) ' Error in drybicki' do k = 1, ntau write(unit=*,fmt=*) ta(k),tb(k),tc(k) end do stop end if end do c ------------------------------------- c Calculate W*T(-1)*V add into F c ------------------------------------- do j1 = 1, ntau do j2 = 1, ntau f(j2,j1)=f(j2,j1)+w(j2)*v(j2,j1) end do end do enddo end if enddo c ---------------------------------------- c Get LU decompostion of E c ---------------------------------------- call ludcmp(e,ntau,nd1,indx,dinx,ierr) if(ierr .eq. 1) then write(unit=*,fmt=*) ' Error in drybicki ' stop 200 end if c ---------------------------------------- c Calculate E(-1)*F store in F c ---------------------------------------- do jt = 1, ntau call lubksb(e,ntau,nd1,indx,f(1,jt)) end do return end