c2--567--1---------2---------3---------4---------5---------6---------712 c Program apriori c c -- Function: Given latitude and day of the year, this program c determines and outputs the 21-layer apriori ozone profile used c in SBUV Version 8 ozone retrieval package c PROGRAM CURRENT AS OF MAY 2004. c c -- Input & Output Variables c real lat !latitude (degree) (input) integer day !day of the year (input) real qu(13) !ozone amounts in Umkehr layers (12:0) (DU) ! (output) real q(81), q21(21), qtot !ozone amounts in fine layers (1:81), 21-layer ! scheme (1:21), and in total column, with ! layer 1 at the bottom of the atmosphere ! (DU) (output) real presly(22) ! Pressure at 21 layer boundaries (atmospheres) data presly / 1., 0.631, 0.398, 0.251, 0.158, 1 .1, 0.0631, 0.0398, 0.0251, 0.0158, 2 .01, 0.00631, 0.00398, 0.00251, 0.00158, 3 0.001,0.000631,0.000398,0.000251,0.000158, 4 0.0001,0./ c23456789112345678921234567893123456789412345678951234567896123456789712 c c c -- Local Variables c real*4 qu_ref(13,18,12) !climatological ozone amounts in 13 Umkehr layers ! arranged from layer 12 to layer 0, in 18 10-degree ! latitude bins from south to north, and ! 12 months, January to December (DU) (input) real cntr(18) data cntr /-85,-75,-65,-55,-45,-35,-25,-15,-5, 1 5,15,25,35,45,55,65,75,85/ real mon, fm, fl integer m1, m2, l1, l2, iu, i, i1 real ptop, slope, lam1, lam2 integer nu, n logical normal parameter (ptop = 1e-4, lam1 = 0.0, lam2 = 0.0, 1 nu = 13, n = 81, normal = .false.) c c Read Umkehr-layer ozone climatology data c open (unit=1, file='climat.txt', status='old', form='formatted') read (1, '(13f7.3)') qu_ref close (1) c c Input latitude and day c write (*, *) write (*, *) 'SBUV V8 Apriori' write (*, *) write (*, *) 1 'When prompted, input latitude (degree) and day of the year ', 2 '(in integer) and press enter: The output is latitude- and ', 3 'season-dependent 21-layer apriori ozone profile and total ', 4 'column ozone, in DU, with layer bottoms at 1.0, 0.631, 0.398, ', 5 '0.251, 0.158, 0.1, 0.0631, 0.0398, 0.0251, 0.0158, 0.01, ', 6 '0.00631, 0.00398, 0.00251, 0.00158, 0.001, 0.000631, ', 7 '0.000398, 0.000251, 0.000158, 0.0001 atmos, and top layer ', 8 'extending to infinity. ' write (*, *) write (*, *) 1 'For another latitude and day, repeat the process; to end ', 2 'the program, press cntl-c. ' write (*, *) 100 write (*,*) 'lat, day ?' read (*,*) lat, day if ((lat .lt. -90.0) .or. (lat .gt. 90.0)) then write (*, *) 'lat must lie between -90.0 and 90.0' go to 100 end if if ((day .lt. 1) .or. (day .gt. 366)) then write (*, *) 'day should be an integer between 1 and 366' go to 100 end if c c Perform bi-linear interpolation in latitude and month c l1 = int((lat + 85.0)/10.0) + 1 if (l1 .le. 0) l1 = 1 if (l1 .ge. 18) l1 = 17 l2 = l1 + 1 fl = (lat - cntr(l1))/(cntr(l2) - cntr(l1)) mon = (day + 15.25)/30.5 m1 = int(mon) m2 = m1 + 1 fm = mon - m1 if (m1 .eq. 0) m1 = 12 if (m2 .eq. 13) m2 = 1 do iu = 1,13 qu(iu)=(1.0 - fl)*(1.0 - fm)*qu_ref(iu,l1,m1) 1 + fl*(1.0 - fm)*qu_ref(iu,l2,m1) 2 + (1.0 - fl)*fm*qu_ref(iu,l1,m2) 3 + fl*fm*qu_ref(iu,l2,m2) end do c c Perform interpolation for 81-layer profile c call interpol_qo3(qu, ptop, q, slope, nu, n, lam1, lam2, normal) c c Convert 81-layer profile to 21-layer profile c qtot = 0.0 do i = 1,20 i1 = (i-1)*4 + 1 q21(i) = q(i1) + q(i1+1) + q(i1+2) + q(i1+3) qtot = qtot + q21(i) end do q21(21) = q(81) qtot = qtot + q21(21) c c Output 21-layer profile c write (*, *) 'profile (bottom-to-top):' write (*, '(2x, 10f7.3)') q21 write (*, *) 'column:' write (*, '(2x, f7.3)') qtot write (*, *) go to 100 c end c2--567--1---------2---------3---------4---------5---------6---------712 Subroutine interpol_qo3(q, p_top, q_fine, slope, n, n_fineplus, 1 lam1, lamn, normal) c c -- Interpol_qo3/Interpolation of layer ozone from Umkehr to fine c layers/N. Nath/June '01 c c -- Function: This routine interpolates layer ozone amount (q) c data from Umkehr layers to fine layers, as a preparation for c computing back-scattered solar radiation from the SBUV c experiment. The fine layers, like Umkehr layers, are of equal c delta[log(p)], where p is the pressure, and the output data are c in the same unit as Umkehr data (DU). The Umkehr data may be c given from bottom to top of the atmosphere (i.e., starting with c layer 0), or in the reverse order. The fine layer data are c from bottom to top. Spline interpolation of log(q) is used for c the lower layers, and linear interpolation is used for the top c layers starting at the middle of the next-to-top Umkehr layer. c Following the interpolation, the amount in the topmost layer, c which extends to infinity, is calculated based on the assumed c exponential drop of q (same as linear drop of log(q)) in top c fine layers. (Check if q's are positive before entering this c routine.) c c -- Input & Output Variables: c real q(n), p_top !layer ozone amount in Umkehr layers (DU), ! atmospheric pressure at top of the ! uppermost fine layer (atmos) (input) real q_fine(n_fineplus), slope !interpolated ozone amount in fine layers ! (DU) (includes amount in topmost layer, ! an extension of fine layers extending ! to infinity--computed at the end after ! interpolation), slope of log[q_fine] ! vs. z at the top layers (output) integer n, n_fineplus !number of Umkehr layers & fine layers ! (including the layer extending to ! infinity) (input) real lam1, lamn !end parameters for spline interpolation ! (input) logical normal!if true, Umkehr layer data are ordered ! from 0 upward; if false, the opposite ! (input) c c -- Local Variables: c integer nl, nl_fine, ns !maximum values of n and n_fine, & number ! of auxiliary q's parameter (nl = 13, nl_fine = 81, ns = 8) c real h !pressure scale height at z = 0 & standard ! temperature (cm) (i.e., p proportional ! to exp(-z/h)) parameter (h = 7.996e5) real q0(nl), z(nl), dz, z_fine(nl_fine), dz_fine !q save array, z array & dz for Umkehr ! layers, z array & dz for fine layers real logqs(nl), logq_fine(nl_fine) !log(scaled-q), log(q_fine) real r !ratio of successive q_fine's at the top real z_aux(nl*ns), logq_aux(nl*ns), q_aux(nl*ns), 1 dz_aux, qsum(nl), sum !quantities needed for auxiliary computns real devn2, devn, tol parameter (tol = 1e-4) integer n_fine, n_spline, i, is, n_aux, irpt c c Order data from bottom to top c if (.not. normal) then do i = 1,n q0(i) = q(i) end do do i = 1,n q(i) = q0(n+1-i) end do end if c c Prepare for interpolation c dz = h*log(2.0) !define dz & z array do i = 1,n z(i) = dz/2.0 + (i - 1)*dz end do c n_fine = n_fineplus - 1 !number of actual fine layers 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 = 0 !count how many z_fine's are below z(n-1), ! the extent of spline interpolation do i = 1,n_fine if (z_fine(i) .gt. z(n-1)) go to 100 n_spline = n_spline + 1 end do 100 continue c c Scale q to fine layers & convert to log c do i = 1,n-1 !note: layer n being excluded (see later) logqs(i) = log(q(i)*dz_fine/dz) end do c c o Iterate on the scaling so the layer ozone amount is c unchanged c dz_aux = dz/ns n_aux = (n - 1)*ns !note: layer n being excluded do i = 1,n_aux z_aux(i) = dz_aux/2.0 + (i - 1)*dz_aux end do do irpt = 1,8 !four iterations should suffice call spline(z, logqs, n-1, z_aux, logq_aux, n_aux, lam1, lamn) !first step in finding how the choice of ! qs affects q in Umkehr layers, as a ! result of spline interpolation do i = 1,n_aux q_aux(i) = exp(logq_aux(i)) end do do i = 1,n-1 sum = 0.0 do is = 1,ns sum = sum + q_aux((i-1)*ns + is) end do qsum(i) = sum*dz_aux/dz_fine !qsum(i) should be close to q(i) logqs(i) = logqs(i) + log(q(i)/qsum(i)) !correction so in the next step qsum(i) ! will be closer to q(i) end do devn2 = 0.0 do i = 1,n-1 devn = abs(1.0 - qsum(i)/q(i)) devn2 = devn2 + devn**2 end do devn = sqrt(devn2/(n-1)) if (devn .lt. tol) go to 200 end do 200 continue c c Perform spline interpolation for lower layers c call spline(z, logqs, n-1, z_fine, logq_fine, n_spline, lam1, 1 lamn) !note: q(n) is excluded, as that includes ! ozone up to pressure zero and is not ! representative of the layer itself c c Perform linear interpolation for top layers c slope = (logqs(n-1) - logqs(n-2))/(z(n-1) - z(n-2)) !note: q(n) is excluded do i = n_spline+1,n_fine logq_fine(i) = logqs(n-1) + slope*(z_fine(i) - z(n-1)) end do c c Convert log(q_fine) to q_fine c do i = 1,n_fine q_fine(i) = exp(logq_fine(i)) end do c c Compute q at the topmost layer c r = q_fine(n_fine)/q_fine(n_fine-1) q_fine(n_fineplus) = r*q_fine(n_fine)/(1.0 - r) !constant ratio r in successive layers, ! with (1 - r)**-1 coming from the ! series 1 + r + r**2 + ... c c Return original data c if (.not. normal) then do i = 1,n q(i) = q0(i) end do end if c return end c2--567--1---------2---------3---------4---------5---------6---------712 Subroutine spline(x0, y0, n0, x, y, n, lam1, lamn0) c c -- Spline/Spline interpolation/N. Nath/April '01. (Adapted from c a previous version, this routine has the following features: c (1) it incorporates different types of end conditions; (2) the c argument array {x0} may be an increasing or a decreasing c sequence; (3) {x} need not be ordered; (4) for x lying outside c {x0}, the cubic in the nearest interval is used; this feature c should be used with extreme caution.) c c -- Function: Performs a cubic spline interpolation for function c values (y) at an array of arguments (x) of dimension n, with c different end parameters lam1 & lamn0. x0 & y0 (of dimension c n0) are the input argument & associated function array. lam1 c [alternately lamn0] may take the following values: 0.0, -1.0, c -2.0, 1.0, which correspond respectively to y20(1) = (0.0, 0.5, c 1.0)*y20(2), & y10(1) = 0.0 [y20(n0) = (0.0, 0.5, 1.0)* c y20(n0-1), & y10(n0) = 0.0]; y20 and y10 are respectively the c second and first derivative at x0. c c -- Usage Notes: (1) Choice of lam1 and lamn0 depends on the c knowledge of the function. If nothing is known, use one of the c first two values. Note also that the second & third values c correspond respectively to a single cubic in the domain: c (2.0*x0(1) - x0(2), x0(2)), with y20 at the left boundary = c 0.0, & a quadratic in the domain (x0(1), x0(2)), and similarly c for the other end. (2) x0 must be monotonically increasing or c decreasing. (3) n0 must be > 2; if it exceeds 256, change the c value of n0max. c c -- Input & Output Variables: c real x0(n0), y0(n0), !base argument & function array (input), 1 x(n), y(n), !argument array at which function values ! are desired (input), & interpolated ! function array (output) 2 lam1, lamn0 !end parameters at x0(1) & x0(n0) (input) integer n0, n !dimensions of (x0, y0) & (x, y) (input) c c -- Local Variables: c integer n0max parameter (n0max = 256) c real y20(n0max) !second derivative array real h(n0max), a(n0max), b(n0max), d(n0max), 1 phi, lam, dy0, dy00, c, c1, f integer n01, l, u, m, i, j, k, i1 logical reverse c c -- Subprogram Called: order c c -- Ref.: Ahlberg et al, Theory Of Splines & Their Applications, c Academic Press, 1967; Hamming, Numerical Methods, Mcgraw-Hill, c 1973. c phi(lam) = lam*(lam + 1.0)*(lam + 2.0) c c Order x0 in an increasing sequence; change x accordingly c !note: this ordering, viz., x0 -> -x0 while y0 -> y0, c ! only changes the signs of the odd derivatives; thus c ! y20 remains unchanged c reverse = (x0(1) .gt. x0(n0)) if (reverse) call order(x0, n0, x, n) c c Determine y20 by solving tri-diagonal linear equations c c o Define h, a, & d (quantities independent of y0) c n01 = n0 - 1 do i = 1,n01 h(i) = x0(i+1) - x0(i) end do a(1) = -lam1/2.0 do i = 2,n01 i1 = i - 1 d(i) = h(i1)*a(i1) + 2.0*(h(i1) + h(i)) a(i) = -h(i)/d(i) end do d(n0) = h(n01)*(lamn0*a(n01) + 2.0) c c o Define b, & complete computations c dy0 = (y0(2) - y0(1))/h(1) b(1) = (phi(lam1)*dy0)/h(1)/2.0 do i = 2,n01 dy00 = dy0 dy0 = (y0(i+1) - y0(i))/h(i) b(i) = (6.0*(dy0 - dy00) - h(i-1)*b(i-1))/d(i) end do b(n0) = (-phi(lamn0)*dy0 - lamn0*h(n01)*b(n01))/d(n0) y20(n0) = b(n0) k = n0 do i = 2,n0 k = k - 1 y20(k) = a(k)*y20(k+1) + b(k) end do c c Interpolate using binary domain search c do j = 1,n l = 1 u = n0 m = (l + u)/2 do while (m .gt. l) if (x(j) .lt. x0(m)) then u = m else l = m end if m = (l + u)/2 end do c = (x0(l+1) - x(j))/h(l) c1 = 1.0 - c f = (h(l)**2)/6.0 y(j) = c*(y0(l) - f*y20(l)*(1.0 - c**2)) 1 + c1*(y0(l+1) - f*y20(l+1)*(1.0 - c1**2)) end do c c Reorder x0 and x c if (reverse) call order(x0, n0, x, n) c return end c --- Subroutine order(x0, n0, x, n) c c -- Order/Adjunct to Spline/N. Nath/April '01 c c -- Input & Output Variables: c real x0(n0), x(n) integer n0, n c c -- Local Variables: c integer i, j c c Change signs of x0 and x c do i = 1,n0 x0(i) = -x0(i) end do do j = 1,n x(j) = -x(j) end do c return end