c2--567--1---------2---------3---------4---------5---------6---------712 Subroutine interpol_kern(kern, p_top, p_bottom, kern_fine, nc, 1 nc1, n, n_fine, lam1, lamn) c c -- Interpol_kern/Interpolation of multiple scattering kernel from c Umkehr to fine layers/N. Nath/September '01 c c Modification: Incorporation of a lower boundary of the c atmosphere at p_bottom/N. Nath/October '01 c c -- Function: Given multiple scattering kernel, d(nval_msr)/dq, in c Umkehr layers (where nval_msr is proportional to log(rad_msr) c with total radiance rad = rad_ss + rad_msr, and q is the ozone c amount in those layers), this routine determines multiple c scattering kernel at equidistant, fine log(pressure) levels for c ozone profile retrieval calculations. Umkehr layers are to be c specified from bottom to top of the atmosphere. As written, c the routine is meant to be used for n = 11, but it can be used c as such for n = 13 (as in the routines 'interpol_qo3' & c 'interpol_t'). An extrapolation to a constant value is used at c the top, as the kernel is essentially constant there. c c -- Input/Output Note: The input kernel at the Umkehr layer c including p_bottom is a weighted average of a higher value for c the part of the layer above p_bottom and zero for the lower c part. Similarly, the output fine-layer kernel for the layer c including p_bottom is a weighted average of a higher value for c the top part and zero for the lower. Interpolation within this c Umkehr layer involves a transformation from the average kernel c to the local kernel for the upper part; a reverse trans- c formation yields the fine layer average kernel. c c -- Input & Output Variables: c real kern(nc,n), p_top !multiple scattering kernel in Umkehr ! layers (arbitrary unit), pressure at ! top of fine layers (atmos) (input) real p_bottom !pressure at lower boundary of atmosphere ! [the additional parameter] (atmos) ! (input) real kern_fine(nc,n_fine) !interpolated kernel in fine layers ! (output) integer nc, nc1, n, n_fine !total number of channels, starting ! channel number with shortest wavelength ! (nominal 5), number of Umkehr layers, ! number of fine layers excluding the top ! layer extending to infinity (input) real lam1, lamn !end parameters for spline interpolation ! (input) c c -- Local Variables: c integer nl, nl_fine, ns !maximum values of n & n_fine, & number of ! auxiliary kern's parameter (nl = 13, nl_fine = 80, ns = 8) c real h !pressure scale height at z = 0 & standard ! temperature (cm) parameter (h = 7.996e5) real z(nl), dz, z_fine(nl_fine), dz_fine, z_bottom !z array & dz for Umkehr layers, z array ! & dz for fine layers, z at p_bottom real kernmid(nl), kernx(nl_fine) !best guess of Umkehr mid-layer kernel ! (such that the spline average agrees with ! input kern), work variable real z_aux(nl*ns), kern_aux(nl*ns), dz_aux, f, 1 f_fine, kernav(nl), sum !quantities needed for auxiliary computns real devn2, devn, tol parameter (tol = 1e-5) c real maxdevn c integer imax integer ic, i, n_spline, nb, nbfine, is, n_aux, irpt c c Prepare for interpolation c dz = h*log(2.0) !define dz, z array, & z_bottom do i = 1,n z(i) = dz/2.0 + (i - 1)*dz end do z_bottom = -h*log(p_bottom) dz_fine = - h*log(p_top)/n_fine !define dz_fine & z_fine array do i = 1,n_fine z_fine(i) = dz_fine/2.0 + (i - 1)*dz_fine end do n_spline = int(z(n)/dz_fine) !number of z_fine's below z(n), the ! extent of spline interpolation nb = int(z_bottom/dz) + 1 !beginning Umkeyr layer where kern is ! non-zero do ic = nc1,nc do i = 1,nb-1 if (kern(ic,i) .ne. 0.0) then write (*, *) "Error: non-zero kern below lower boundary" write (*, *) " at channel ", ic, ", layer ", i end if end do end do nbfine = int(z_bottom/dz_fine) + 1 !beginning fine layer where kern is ! non-zero c c Set kern_fine = 0 for layers below lower boundary c do ic = nc1,nc do i = 1,nbfine-1 kern_fine(ic,i) = 0.0 end do end do c c Interpolate c dz_aux = dz/ns n_aux = (n-nb+1)*ns f = (z(nb) + dz/2.0 - z_bottom)/dz !fractional thickness of first Umkehr ! layer with non-zero kern if (f .eq. 0.0) then write (*, *) "Error w/ interpol: frac thickness equals zero" end if f_fine = (z_fine(nbfine) + dz_fine/2.0 - z_bottom)/dz_fine !fractional thickness of first fine layer ! with non-zero kern do i = 1,ns z_aux(i) = z_bottom + (i - 0.5)*f*dz_aux !z_aux starts at z_bottom, & first ns of ! them are in layer nb with reduced ! separation compared to the rest end do do i = ns+1,n_aux z_aux(i) = z(nb) + dz/2.0 + (i - ns - 0.5)*dz_aux !regular separation end do do ic = nc1,nc c c o Determine Umkehr mid-layer kernel such that the average c of the resulting spline agrees with kern c kernmid(nb) = kern(ic,nb)/f !average to local transformation ! (by definition f cannot be zero) do i = nb+1,n kernmid(i) = kern(ic,i) end do c c - Iterate until the desired accuracy is attained c c do i = nb,n c if (kernmid(i) .eq. 0.0) go to 100 c !skip iteration if any of kernmid is 0 c ! (this condition is symptomatic of layer c ! 5, for which kern is small, & therefore c ! refinement of kernmid is c ! inconsequential) c end do do irpt = 1,8 call spline(z(nb), kernmid(nb), n-nb+1, 1 z_aux, kern_aux, n_aux, lam1, lamn) do i = nb,n sum = 0.0 do is = 1,ns sum = sum + kern_aux((i-nb)*ns + is) end do kernav(i) = sum/ns if (i .eq. nb) kernav(i) = kernav(i)*f !kernav(i) should be close to kern(ic,i) c kernmid(i) = kernmid(i)*(kern(ic,i)/kernav(i)) if (i .eq. nb) then kernmid(i) = kernmid(i) + (kern(ic,i)-kernav(i))/f else kernmid(i) = kernmid(i) + (kern(ic,i)-kernav(i)) end if !correction so in the next step kernav(i) ! will be closer to kern(ic,i) !note: kern & kernav are layer average, ! whereas kernmid is local; division by f ! in modified form of kernmid correction ! is tied with that end do devn2 = 0.0 c imax = 0 c maxdevn = 0.0 do i = nb,n devn = abs(kernav(i) - kern(ic,i)) devn2 = devn2 + devn**2 c if (devn .gt. maxdevn) then c imax = i c maxdevn = devn c end if end do devn = sqrt(devn2/(n-nb+1)) c write (*, *) kernav c write (*, *) "devn, imax, maxdevn: ", devn, imax, maxdevn if (devn .lt. tol) go to 100 end do 100 continue c c o Perform spline interpolation up to middle of the top c Umkehr layer c call spline(z(nb), kernmid(nb), n-nb+1, 1 z_fine(nbfine), kernx(nbfine), n_spline-nbfine+1, 2 lam1, lamn) !extrapolation at lower end kern_fine(ic,nbfine) = kernx(nbfine)*f_fine !re-transformation from local to average do i = nbfine+1,n_spline kern_fine(ic,i) = kernx(i) end do c c o Set kernels for top fine layers equal to kern(ic,n) c do i = n_spline+1,n_fine kern_fine(ic,i) = kern(ic,n) end do end do c return end