MODULE module_ftuv_subs use module_wave_data, only : nw implicit none integer, parameter :: dp = selected_real_kind(14,300) !----------------------------------------------------------------------------- ! data for inter2, inter3 and interi !----------------------------------------------------------------------------- integer, private, save :: nintervals integer, private, allocatable, save :: xi(:), xcnt(:) real(dp), private, allocatable, save :: xfrac(:,:) real(dp), private, parameter :: km2cm = 1.e5_dp !----------------------------------------------------------------------------- ! data for lyman !----------------------------------------------------------------------------- integer, private, parameter :: nla = 2 integer, private, parameter :: ngast = 17 !----------------------------------------------------------------------------- ! data for schu_inti, schu !----------------------------------------------------------------------------- integer, private, parameter :: tdim = 501 !BSINGH(PNNL)-Lahey compiler doesn't allow use of 'real' in variable declaration real(dp), private, parameter :: tdim_real = 501.0 !real(dp), private, parameter :: t_del = 5._dp/real(tdim-1,kind=dp)!BSINGH(PNNL)-commented out !real(dp), private, parameter :: t_fac = real(tdim-1,kind=dp)/5._dp!BSINGH(PNNL)-commented out real(dp), private, parameter :: t_del = 5._dp/(tdim_real-1._dp)!BSINGH(PNNL)-redefined t_del and t_fac real(dp), private, parameter :: t_fac = (tdim_real-1._dp)/5._dp integer, private :: ii, jj real(dp), private, save :: d_table(0:tdim,ngast-1) real(dp), private, save :: x_table(0:tdim,ngast-1) real(dp), private, save :: o2_table(tdim) real(dp), private, dimension(12,ngast-1), save :: a_schu, b_schu !----------------------------------------------------------------------------- ! a_schu(16,12) coefficients for rj(m) (table 1 in kockarts 1994) ! b_schu(16,12) rj(o2)(table 2 in kockarts 1994) ! rjm attenuation coefficients rj(m) ! rjo2 rj(o2) !----------------------------------------------------------------------------- data ((a_schu(jj,ii),jj=1,12),ii=1,ngast-1) / & !a 57000-56500.5 cm-1 1.13402e-01,1.00088e-20,3.48747e-01,2.76282e-20,3.47322e-01,1.01267e-19, & 1.67351e-01,5.63588e-19,2.31433e-02,1.68267e-18,0.00000e+00,0.00000e+00, & !a 56500-56000.5 cm-1 2.55268e-03,1.64489e-21,1.85483e-01,2.03591e-21,2.60603e-01,4.62276e-21, & 2.50337e-01,1.45106e-20,1.92340e-01,7.57381e-20,1.06363e-01,7.89634e-19, & !a 56000-55500.5 cm-1 4.21594e-03,8.46639e-22,8.91886e-02,1.12935e-21,2.21334e-01,1.67868e-21, & 2.84446e-01,3.94782e-21,2.33442e-01,1.91554e-20,1.63433e-01,2.25346e-19, & !a 55500-55000.5 cm-1 3.93529e-03,6.79660e-22,4.46906e-02,9.00358e-22,1.33060e-01,1.55952e-21, & 3.25506e-01,3.43763e-21,2.79405e-01,1.62086e-20,2.10316e-01,1.53883e-19, & !a 55000-54500.5 cm-1 2.60939e-03,2.33791e-22,2.08101e-02,3.21734e-22,1.67186e-01,5.77191e-22, & 2.80694e-01,1.33362e-21,3.26867e-01,6.10533e-21,1.96539e-01,7.83142e-20, & !a 54500-54000.5 cm-1 9.33711e-03,1.32897e-22,3.63980e-02,1.78786e-22,1.46182e-01,3.38285e-22, & 3.81762e-01,8.93773e-22,2.58549e-01,4.28115e-21,1.64773e-01,4.67537e-20, & !a 54000-53500.5 cm-1 9.51799e-03,1.00252e-22,3.26320e-02,1.33766e-22,1.45962e-01,2.64831e-22, & 4.49823e-01,6.42879e-22,2.14207e-01,3.19594e-21,1.45616e-01,2.77182e-20, & !a 53500-53000.5 cm-1 7.87331e-03,3.38291e-23,6.91451e-02,4.77708e-23,1.29786e-01,8.30805e-23, & 3.05103e-01,2.36167e-22,3.35007e-01,8.59109e-22,1.49766e-01,9.63516e-21, & !a 53000-52500.5 cm-1 6.92175e-02,1.56323e-23,1.44403e-01,3.03795e-23,2.94489e-01,1.13219e-22, & 3.34773e-01,3.48121e-22,9.73632e-02,2.10693e-21,5.94308e-02,1.26195e-20, & !a 52500-52000.5 cm-1 1.47873e-01,8.62033e-24,3.15881e-01,3.51859e-23,4.08077e-01,1.90524e-22, & 8.08029e-02,9.93062e-22,3.90399e-02,6.38738e-21,8.13330e-03,9.93644e-22, & !a 52000-51500.5 cm-1 1.50269e-01,1.02621e-23,2.39823e-01,3.48120e-23,3.56408e-01,1.69494e-22, & 1.61277e-01,6.59294e-22,8.89713e-02,2.94571e-21,3.25063e-03,1.25548e-20, & !a 51500-51000.5 cm-1 2.55746e-01,8.49877e-24,2.94733e-01,2.06878e-23,2.86382e-01,9.30992e-23, & 1.21011e-01,3.66239e-22,4.21105e-02,1.75700e-21,0.00000e+00,0.00000e+00, & !a 51000-50500.5 cm-1 5.40111e-01,7.36085e-24,2.93263e-01,2.46742e-23,1.63417e-01,1.37832e-22, & 3.23781e-03,2.15052e-21,0.00000e+00,0.00000e+00,0.00000e+00,0.00000e+00, & !a 50500-50000.5 cm-1 8.18514e-01,7.17937e-24,1.82262e-01,4.17496e-23,0.00000e+00,0.00000e+00, & 0.00000e+00,0.00000e+00,0.00000e+00,0.00000e+00,0.00000e+00,0.00000e+00, & !a 50000-49500.5 cm-1 8.73680e-01,7.13444e-24,1.25583e-01,2.77819e-23,0.00000e+00,0.00000e+00, & 0.00000e+00,0.00000e+00,0.00000e+00,0.00000e+00,0.00000e+00,0.00000e+00, & !a 49500-49000.5 cm-1 3.32476e-04,7.00362e-24,9.89000e-01,6.99600e-24,0.00000e+00,0.00000e+00, & 0.00000e+00,0.00000e+00,0.00000e+00,0.00000e+00,0.00000e+00,0.00000e+00 / data ((b_schu(jj,ii),jj=1,12),ii=1,ngast-1) / & ! 57000-56500.5 cm-1 1.07382e-21,9.95029e-21,7.19430e-21,2.48960e-20,2.53735e-20,7.54467e-20, & 4.48987e-20,2.79981e-19,9.72535e-20,9.29745e-19,2.30892e-20,4.08009e-17, & ! 56500-56000.5 cm-1 3.16903e-22,1.98251e-21,5.87326e-22,3.44057e-21,2.53094e-21,8.81484e-21, & 8.82299e-21,4.17179e-20,2.64703e-20,2.43792e-19,8.73831e-20,1.46371e-18, & ! 56000-55500.5 cm-1 1.64421e-23,9.26011e-22,2.73137e-22,1.33640e-21,9.79188e-22,2.99706e-21, & 3.37768e-21,1.39438e-20,1.47898e-20,1.04322e-19,4.08014e-20,6.31023e-19, & ! 55500-55000.5 cm-1 8.68729e-24,7.31056e-22,8.78313e-23,1.07173e-21,8.28170e-22,2.54986e-21, & 2.57643e-21,9.42698e-21,9.92377e-21,5.21402e-20,3.34301e-20,2.91785e-19, & ! 55000-54500.5 cm-1 1.20679e-24,2.44092e-22,2.64326e-23,4.03998e-22,2.53514e-22,8.53166e-22, & 1.29834e-21,3.74482e-21,5.12103e-21,2.65798e-20,2.10948e-20,2.35315e-19, & ! 54500-54000.5 cm-1 2.79656e-24,1.40820e-22,3.60824e-23,2.69510e-22,4.02850e-22,8.83735e-22, & 1.77198e-21,6.60221e-21,9.60992e-21,8.13558e-20,4.95591e-21,1.22858e-17, & ! 54000-53500.5 cm-1 2.36959e-24,1.07535e-22,2.83333e-23,2.16789e-22,3.35242e-22,6.42753e-22, & 1.26395e-21,5.43183e-21,4.88083e-21,5.42670e-20,3.27481e-21,1.58264e-17, & ! 53500-53000.5 cm-1 8.65018e-25,3.70310e-23,1.04351e-23,6.43574e-23,1.17431e-22,2.70904e-22, & 4.88705e-22,1.65505e-21,2.19776e-21,2.71172e-20,2.65257e-21,2.13945e-17, & ! 53000-52500.5 cm-1 9.63263e-25,1.54249e-23,4.78065e-24,2.97642e-23,6.40637e-23,1.46464e-22, & 1.82634e-22,7.12786e-22,1.64805e-21,2.37376e-17,9.33059e-22,1.13741e-20, & ! 52500-52000.5 cm-1 1.08414e-24,8.37560e-24,9.15550e-24,2.99295e-23,9.38405e-23,1.95845e-22, & 2.84356e-22,3.39699e-21,1.94524e-22,2.72227e-19,1.18924e-21,3.20246e-17, & ! 52000-51500.5 cm-1 1.52817e-24,1.01885e-23,1.22946e-23,4.16517e-23,9.01287e-23,2.34869e-22, & 1.93510e-22,1.44956e-21,1.81051e-22,5.17773e-21,9.82059e-22,6.22768e-17, & ! 51500-51000.5 cm-1 2.12813e-24,8.48035e-24,5.23338e-24,1.93052e-23,1.99464e-23,7.48997e-23, & 4.96642e-22,6.15691e-17,4.47504e-23,2.76004e-22,8.26788e-23,1.65278e-21, & ! 51000-50500.5 cm-1 3.81336e-24,7.32307e-24,5.60549e-24,2.04651e-23,3.36883e-22,6.15708e-17, & 2.09877e-23,1.07474e-22,9.13562e-24,8.41252e-22,0.00000e+00,0.00000e+00, & ! 50500-50000.5 cm-1 5.75373e-24,7.15986e-24,5.90031e-24,3.05375e-23,2.97196e-22,8.92000e-17, & 8.55920e-24,1.66709e-17,0.00000e+00,0.00000e+00,0.00000e+00,0.00000e+00, & ! 50000-49500.5 cm-1 6.21281e-24,7.13108e-24,3.30780e-24,2.61196e-23,1.30783e-22,9.42550e-17, & 2.69241e-24,1.46500e-17,0.00000e+00,0.00000e+00,0.00000e+00,0.00000e+00, & ! 49500-49000.5 cm-1 6.81118e-24,6.98767e-24,7.55667e-25,2.75124e-23,1.94044e-22,1.45019e-16, & 1.92236e-24,3.73223e-17,0.00000e+00,0.00000e+00,0.00000e+00,0.00000e+00 / !----------------------------------------------------------------------------- ! data for set_o2_xsect !----------------------------------------------------------------------------- integer, public, parameter :: nwint = 105 real(dp), public, save :: wlint(nwint) real(dp), private, save :: xso2int(nwint) real(dp), private, save :: wlla(nla) real(dp), private, save :: wlgast(ngast) !----------------------------------------------------------------------------- ! o2 parameters !----------------------------------------------------------------------------- data (wlint(ii),ii=1,105)/ & 0.0000E+00, 0.1166E+03, 0.1167E+03, 0.1173E+03, 0.1179E+03, & 0.1187E+03, 0.1194E+03, 0.1202E+03, 0.1208E+03, 0.1214E+03, & 0.1219E+03, 0.1223E+03, 0.1231E+03, 0.1238E+03, 0.1246E+03, & 0.1254E+03, 0.1262E+03, 0.1270E+03, 0.1278E+03, 0.1286E+03, & 0.1294E+03, 0.1303E+03, 0.1311E+03, 0.1320E+03, 0.1329E+03, & 0.1338E+03, 0.1346E+03, 0.1356E+03, 0.1365E+03, 0.1374E+03, & 0.1384E+03, 0.1399E+03, 0.1418E+03, 0.1439E+03, 0.1459E+03, & 0.1481E+03, 0.1504E+03, 0.1526E+03, 0.1550E+03, 0.1574E+03, & 0.1600E+03, 0.1626E+03, 0.1653E+03, 0.1681E+03, 0.1709E+03, & 0.1731E+03, 0.1746E+03, 0.1754E+03, 0.1770E+03, 0.1786E+03, & 0.1802E+03, 0.1818E+03, 0.1835E+03, 0.1852E+03, 0.1869E+03, & 0.1887E+03, 0.1905E+03, 0.1923E+03, 0.1942E+03, 0.1961E+03, & 0.1980E+03, 0.2000E+03, 0.2020E+03, 0.2041E+03, 0.2050E+03, & 0.2060E+03, 0.2070E+03, 0.2080E+03, 0.2090E+03, 0.2100E+03, & 0.2110E+03, 0.2120E+03, 0.2130E+03, 0.2140E+03, 0.2150E+03, & 0.2160E+03, 0.2170E+03, 0.2180E+03, 0.2190E+03, 0.2200E+03, & 0.2210E+03, 0.2220E+03, 0.2230E+03, 0.2240E+03, 0.2250E+03, & 0.2260E+03, 0.2270E+03, 0.2280E+03, 0.2290E+03, 0.2300E+03, & 0.2310E+03, 0.2320E+03, 0.2330E+03, 0.2340E+03, 0.2350E+03, & 0.2360E+03, 0.2370E+03, 0.2380E+03, 0.2390E+03, 0.2400E+03, & 0.2410E+03, 0.2420E+03, 0.2430E+03, 0.2430E+03, 0.1000E+39 / !----------------------------------------------------------------------------- ! o2 parameters !----------------------------------------------------------------------------- data (xso2int(ii),ii=1,105)/ & 0.0000E+00, 0.1000E-19, 0.6350E-18, 0.7525E-18, 0.1425E-18, & 0.2025E-18, 0.2413E-17, 0.6400E-17, 0.5251E-17, 0.7329E-18, & 0.3445E-18, 0.3425E-18, 0.3925E-18, 0.8918E-17, 0.9197E-17, & 0.6625E-18, 0.2700E-18, 0.1575E-18, 0.3240E-18, 0.4990E-18, & 0.4875E-18, 0.5525E-18, 0.1068E-17, 0.1850E-17, 0.2275E-17, & 0.3425E-17, 0.5890E-17, 0.8365E-17, 0.1090E-16, 0.1275E-16, & 0.1340E-16, 0.1380E-16, 0.1440E-16, 0.1445E-16, 0.1350E-16, & 0.1220E-16, 0.1071E-16, 0.9075E-17, 0.7410E-17, 0.5775E-17, & 0.4210E-17, 0.2765E-17, 0.1655E-17, 0.9760E-18, 0.5900E-18, & 0.3660E-18, 0.2061E-18, 0.3889E-19, 0.4092E-20, 0.2273E-20, & 0.1256E-20, 0.6902E-21, 0.3750E-21, 0.2097E-21, 0.1226E-21, & 0.6508E-22, 0.4579E-22, 0.3220E-22, 0.2507E-22, 0.1941E-22, & 0.1609E-22, 0.1319E-22, 0.9907E-23, 0.7419E-23, 0.7240E-23, & 0.7090E-23, 0.6955E-23, 0.6770E-23, 0.6595E-23, 0.6375E-23, & 0.6145E-23, 0.5970E-23, 0.5805E-23, 0.5655E-23, 0.5470E-23, & 0.5240E-23, 0.5005E-23, 0.4760E-23, 0.4550E-23, 0.4360E-23, & 0.4175E-23, 0.3990E-23, 0.3780E-23, 0.3560E-23, 0.3330E-23, & 0.3095E-23, 0.2875E-23, 0.2875E-23, 0.2875E-23, 0.2700E-23, & 0.2530E-23, 0.2340E-23, 0.2175E-23, 0.2020E-23, 0.1860E-23, & 0.1705E-23, 0.1555E-23, 0.1410E-23, 0.1280E-23, 0.1160E-23, & 0.1055E-23, 0.9550E-24, 0.4500E-24, 0.0000E+00, 0.0000E+00 / !----------------------------------------------------------------------------- ! the lyman-alpha grid !----------------------------------------------------------------------------- data (wlla(ii),ii=1,2)/ 121.4_dp, 121.9_dp / !----------------------------------------------------------------------------- ! o2 parameters !----------------------------------------------------------------------------- data (wlgast(ii),ii=1,17)/ & 0.1754E+03_dp, 0.1770E+03_dp, 0.1786E+03_dp, 0.1802E+03_dp, 0.1818E+03_dp, & 0.1835E+03_dp, 0.1852E+03_dp, 0.1869E+03_dp, 0.1887E+03_dp, 0.1905E+03_dp, & 0.1923E+03_dp, 0.1942E+03_dp, 0.1961E+03_dp, 0.1980E+03_dp, 0.2000E+03_dp, & 0.2020E+03_dp, 0.2041E+03_dp / !----------------------------------------------------------------------------- ! setaer arrays !----------------------------------------------------------------------------- real(dp) :: op_cb1w(nw-1), om_cb1w(nw-1), g_cb1w(nw-1) real(dp) :: op_cb2w(nw-1), om_cb2w(nw-1), g_cb2w(nw-1) real(dp) :: op_oc1w(nw-1), g_oc1w(nw-1) real(dp) :: op_oc2w(nw-1), om_oc2w(nw-1), g_oc2w(nw-1) real(dp) :: op_antw(nw-1), g_antw(nw-1) real(dp) :: op_so4w(nw-1) real(dp) :: op_salw(nw-1) CONTAINS subroutine photoin( chem_opt, nz, zen, o3toms, esfact, & o3top, o2top, albedo, z, tlev, & tlay, airlev, rh, xlwc, o3, & acb1, acb2, aoc1, aoc2, aant, & aso4, asal, tauaer300, tauaer400, tauaer600, & tauaer999, waer300, waer400, waer600, & waer999, gaer300, gaer400, gaer600, & gaer999, aer_ra_feedback, prate, adjcoe, radfld ) use module_wave_data, only : nw, tuv_jmax, deltaw, sflx, & c20, c40, c60, c80, sq !----------------------------------------------- ! ... inter-face routine !----------------------------------------------- implicit none !---------------------------------------------------------- ! ... dummy arguments !---------------------------------------------------------- integer, intent(in) :: chem_opt integer, intent(in) :: nz real(dp), intent(in) :: zen real(dp), intent(in) :: o3toms real(dp), intent(in) :: esfact real(dp), intent(in) :: o3top real(dp), intent(in) :: o2top real(dp), intent(in) :: albedo(nw-1) real(dp), intent(in) :: z(nz) real(dp), intent(in) :: tlev(nz) real(dp), intent(in) :: tlay(nz-1) real(dp), intent(in) :: airlev(nz) real(dp), intent(in) :: rh(nz) real(dp), intent(in) :: xlwc(nz) real(dp), intent(in) :: o3(nz) real(dp), intent(in) :: acb1(nz) real(dp), intent(in) :: acb2(nz) real(dp), intent(in) :: aoc1(nz) real(dp), intent(in) :: aoc2(nz) real(dp), intent(in) :: aant(nz) real(dp), intent(in) :: aso4(nz) real(dp), intent(in) :: asal(nz) ! rajesh: declare aerosol optical properties real(dp), intent(in) :: tauaer300(nz-1) real(dp), intent(in) :: tauaer400(nz-1) real(dp), intent(in) :: tauaer600(nz-1) real(dp), intent(in) :: tauaer999(nz-1) real(dp), intent(in) :: waer300(nz-1) real(dp), intent(in) :: waer400(nz-1) real(dp), intent(in) :: waer600(nz-1) real(dp), intent(in) :: waer999(nz-1) real(dp), intent(in) :: gaer300(nz-1) real(dp), intent(in) :: gaer400(nz-1) real(dp), intent(in) :: gaer600(nz-1) real(dp), intent(in) :: gaer999(nz-1) INTEGER, INTENT(IN) :: aer_ra_feedback real(dp), intent(out) :: prate(nz,tuv_jmax) real(dp), intent(out) :: adjcoe(nz,tuv_jmax) real(dp), intent(out) :: radfld(nz,nw-1) !---------------------------------------------------------- ! ... local variables !---------------------------------------------------------- integer :: i, j, k, iw, n real(dp) :: factor, delzint character(len=8 ) :: cdate(4) character(len=10) :: ctime(4) !---------------------------------------------------------- ! ... altitude grid !---------------------------------------------------------- integer :: iz, izi real(dp) :: colinc(nz) real(dp) :: vcol(nz) real(dp) :: scol(nz) real(dp) :: to3(nz) !---------------------------------------------------------- ! ... fast-j adj !---------------------------------------------------------- real(dp), dimension(nz,tuv_jmax) :: adjcoe1, adjcoe2 !---------------------------------------------------------- ! ... solar zenith angle ! slant pathlengths in spherical geometry !---------------------------------------------------------- integer :: nid(0:nz-1) real(dp) :: dsdh(0:nz-1,nz-1) !---------------------------------------------------------- ! ... extra terrestrial solar flux and earth-sun distance ^-2 !---------------------------------------------------------- real(dp), dimension(nw-1) :: etf, delw, xsec !-------------------------------------------------------------- ! ... o2 absorption cross section !-------------------------------------------------------------- real(dp) :: xso2(nw-1,nz) !-------------------------------------------------------------- ! ... atmospheric optical parameters: !-------------------------------------------------------------- real(dp), dimension(nz-1,nw-1) :: dtrl, dto3, dto2 real(dp), dimension(nz-1,nw-1) :: dtcld, omcld, gcld real(dp), dimension(nz-1,nw-1) :: dtcb1, omcb1, gcb1 real(dp), dimension(nz-1,nw-1) :: dtcb2, omcb2, gcb2 real(dp), dimension(nz-1,nw-1) :: dtoc1, omoc1, goc1 real(dp), dimension(nz-1,nw-1) :: dtoc2, omoc2, goc2 real(dp), dimension(nz-1,nw-1) :: dtant, omant, gant real(dp), dimension(nz-1,nw-1) :: dtso4, omso4, gso4 real(dp), dimension(nz-1,nw-1) :: dtsal, omsal, gsal ! rajesh: declare optical property arrays for aerosols real(dp), dimension(nz-1,nw-1) :: dtaer, omaer, gaer !-------------------------------------------------------------- ! ... spectral irradiance and actinic flux (scalar irradiance): !-------------------------------------------------------------- real(dp), dimension(nz,nw-1) :: radxx !------------------------------------------------------------- ! ... j-values: !------------------------------------------------------------- integer :: m !------------------------------------------------------------- ! ... location and time !------------------------------------------------------------- integer :: iyear, imonth, iday real(dp) :: dtime, ut0 !------------------------------------------------------------- ! ... earth-sun distance effect !------------------------------------------------------------- etf(1:nw-1) = sflx(1:nw-1) * esfact !------------------------------------------------------- ! ... air profile and rayleigh optical depths (inter-face) !------------------------------------------------------- call setair( nz, z, airlev, dtrl, colinc, o2top ) !------------------------------------------------------------- ! ... ozone optical depths (must give temperature) (inter-face) !------------------------------------------------------------- call setozo( nz, z, tlay, dto3, to3, o3, airlev, o3top, o3toms ) !------------------------------------------------------------- ! ... cloud optical depths (must cloud water g/m3) (inter-face) !------------------------------------------------------------- call setcld( nz, z, xlwc, dtcld, omcld, gcld ) !------------------------------------------------------------- ! ... aerosol optical depths ... !------------------------------------------------------------- ! call setaer( nz, z, airlev, rh, acb1, & ! acb2, aoc1, aoc2, aant, aso4, & ! asal, & ! dtcb1, omcb1, gcb1, & ! dtcb2, omcb2, gcb2, & ! dtoc1, omoc1, goc1, & ! dtoc2, omoc2, goc2, & ! dtant, omant, gant, & ! dtso4, omso4, gso4, & ! dtsal, omsal, gsal ) !------------------------------------------------------------ ! rajesh: Conform Aerosol Optical Properties from 4 wavelengths to ! the entire spectra of FTUV !------------------------------------------------------------ ! Initialize aerosol optical properties to zero dtaer(:,:) = 0._dp omaer(:,:) = 0._dp gaer(:,:) = 0._dp if(aer_ra_feedback == 1) then call aer_wrf2ftuv(nz, z, tauaer300, tauaer400, & tauaer600, tauaer999, waer300, & waer400, waer600, waer999, & gaer300, gaer400, gaer600, & gaer999, dtaer, omaer, gaer ) endif !------------------------------------------------------------ ! ... photo-chemical and photo-biological weigting functions. ! for pchem, need to know temperature and pressure profiles. ! output: ! from pchem: sq(nw,nz,nj) - for each reaction jlabel(nj) !------------------------------------------------------------- call pchem( chem_opt, nz, tlev, airlev ) !------------------------------------------------------------- ! ... slant path lengths for spherical geometry !------------------------------------------------------------- call spheres( nz, z, zen, dsdh, nid ) call airmas( nz, z, zen, dsdh, nid, colinc, vcol, scol ) !--------------------------------------------------------------- ! ... modification of coefficent of j-vales function of to3 and zenith !--------------------------------------------------------------- if( zen <= 20.5_dp ) then call setz( nz, to3, tlev, c20, 20, adjcoe ) else if( zen > 20.5_dp .and. zen < 40.5_dp ) then call setz( nz, to3, tlev, c20, 20, adjcoe1 ) call setz( nz, to3, tlev, c40, 40, adjcoe2 ) factor = (zen - 20.5_dp) / 20._dp adjcoe(:nz,:tuv_jmax) = adjcoe1(:nz,:tuv_jmax) & + (adjcoe2(:nz,:tuv_jmax) - adjcoe1(:nz,:tuv_jmax)) * factor else if( zen >= 40.5_dp .and. zen < 60.5_dp ) then call setz( nz, to3, tlev, c40, 40, adjcoe1 ) call setz( nz, to3, tlev, c60, 60, adjcoe2 ) factor = (zen - 40.5_dp) / 20._dp adjcoe(:nz,:tuv_jmax) = adjcoe1(:nz,:tuv_jmax) & + (adjcoe2(:nz,:tuv_jmax) - adjcoe1(:nz,:tuv_jmax)) * factor else if( zen >= 60.5_dp .and. zen < 80._dp ) then call setz( nz, to3, tlev, c60, 60, adjcoe1 ) call setz( nz, to3, tlev, c80, 80, adjcoe2 ) factor = (zen - 60.5_dp) / 19.5_dp adjcoe(:nz,:tuv_jmax) = adjcoe1(:nz,:tuv_jmax) & + (adjcoe2(:nz,:tuv_jmax) - adjcoe1(:nz,:tuv_jmax)) * factor else if( zen >= 80. ) then call setz( nz, to3, tlev, c80, 80, adjcoe ) end if !------------------------------------------------------------------ ! ... effective o2 optical depth (sr bands, must know zenith angle!) ! assign o2 cross section to sq(*,*,1) !------------------------------------------------------------------ call set_o2_xsect( nz, z, colinc, vcol, scol, dto2, xso2 ) sq(1:nw-1,1:nz,1) = xso2(1:nw-1,1:nz) !----------------------------------------------- ! ... main wavelength loop !----------------------------------------------- delw(1:nw-1) = deltaw(1:nw-1) * etf(1:nw-1) !--------------------------------------------------- ! ... monochromatic radiative transfer: ! outputs are fdir, fdn, fup !--------------------------------------------------- call rtlink( nz, & nw, & zen, & albedo, & z, & nid, & dsdh, & dtrl, & ! opt depth of air (rayleigh) dto3, & ! opt depth of o3 dto2, & ! opt depth of o2 dtcld, omcld, gcld, & ! opt depth of cloud ! dtcb1, omcb1, gcb1, & ! dtcb2, omcb2, gcb2, & ! dtoc1, omoc1, goc1, & ! dtoc2, omoc2, goc2, & ! dtant, omant, gant, & ! dtso4, omso4, gso4, & ! dtsal, omsal, gsal, & dtaer, omaer, gaer, & ! pass aerosol optical properties radfld ) ! output of abs. scart flux. !---------------------------------------------------------- ! Interplation the top level !---------------------------------------------------------- radfld(:,:) = max( radfld(:,:),0._dp ) delzint = (z(nz-2) - z(nz-3))/(z(nz-1) - z(nz-2)) do iw = 1,nw-1 radfld(1,iw) = radfld(2,iw) + (radfld(2,iw)-radfld(3,iw))*delzint radfld(1,iw) = max(radfld(1,iw),radfld(2,iw)) end do !---------------------------------------------------------- ! ... j-val calculation ! spherical irradiance (actinic flux) ! as a function of altitude ! convert to quanta s-1 nm-1 cm-2 ! (1.e-4 * (wc*1e-9) / (hc = 6.62e-34 * 2.998e8)) !---------------------------------------------------------- do m = 1,tuv_jmax do iz = 1,nz izi = nz - iz + 1 xsec(:nw-1) = sq(:nw-1,iz,m) * delw(:nw-1) prate(iz,m) = dot_product( radfld(izi,:nw-1), xsec(:nw-1) ) end do prate(1:nz,m) = prate(1:nz,m) * adjcoe(1:nz,m) end do end subroutine photoin !---------------------------------------------------------------------------- !rajesh: subroutine to convert aerosol optical properties from 4 !wavelengths ! to entire spectra of FTUV !----------------------------------------------------------------------------- subroutine aer_wrf2ftuv( nzlev, z, tauaer300, tauaer400, & tauaer600, tauaer999, & waer300, waer400, waer600, waer999, & gaer300, gaer400, gaer600, gaer999, & dtaer, omaer, gaer ) use module_wave_data, only : nw, wc !---------------------------------------------------------------------- ! The routine is based on aerosol treatment in module_ra_rrtmg_sw.F ! INPUT: ! nzlev: number of specified altitude levels in the working grid ! z: specified altitude working grid ! Aerosol optical properties at 300, 400, 600 and 999 nm. ! tauaer300, tauaer400, tauaer600, tauaer999: Layer AODs ! waer300, waer400, waer600, waer999: Layer SSAs ! gaer300, gaer400, gaer600, gaer999: Layer asymmetry parameters ! OUTPUT: ! dtaer: Layer AOD at FTUV wavelengths ! omaer: Layer SSA at FTUV wavelengths ! gaer : Layer asymmetry parameters at FTUV wavelengths !------------------------------------------------------------------------ implicit none !----------------------------------------------------------------------------- ! ... Dummy arguments !----------------------------------------------------------------------------- integer, intent(in) :: nzlev real(dp), intent(in) :: z(nzlev) real(dp), intent(in) :: tauaer300(nzlev-1), tauaer400(nzlev-1), & tauaer600(nzlev-1), tauaer999(nzlev-1) real(dp), intent(in) :: waer300(nzlev-1), waer400(nzlev-1), & waer600(nzlev-1), waer999(nzlev-1) real(dp), intent(in) :: gaer300(nzlev-1), gaer400(nzlev-1), & gaer600(nzlev-1), gaer999(nzlev-1) ! Output arrays real(dp), intent(out) :: dtaer(nzlev-1,nw-1), & omaer(nzlev-1,nw-1), gaer(nzlev-1,nw-1) ! Local Variables integer :: k, wn, nloop real(dp) :: ang, slope real(dp), parameter :: thresh = 1.e-9_dp ang = 0._dp slope = 0._dp ! Start Calculation do wn = 1,nw-1 ! wavelength loop do k = 1,nzlev-1 ! level loop ! use angstrom exponent to calculate aerosol optical depth; wc is in nm. if(tauaer300(k) .gt. thresh .and. tauaer999(k) .gt. thresh) then ang = log(tauaer300(k)/tauaer999(k))/log(0.999_dp/0.3_dp) dtaer(k,wn) = tauaer400(k)*(0.4_dp/(wc(wn)*1.e-3_dp))**ang ! print *, ang, dtaer(k,wn), tauaer600(k)*(600./wc(wn))**ang ! ssa - use linear interpolation/extrapolation slope = (waer600(k)-waer400(k))/0.2_dp omaer(k,wn) = slope*((wc(wn)*1.e-3_dp)-0.6_dp)+waer600(k) if(omaer(k,wn) .lt. 0.4_dp) omaer(k,wn)=0.4_dp if(omaer(k,wn) .ge. 1.0_dp) omaer(k,wn)=1.0_dp ! asymmetry parameter - use linear interpolation/extrapolation slope = (gaer600(k)-gaer400(k))/0.2_dp gaer(k,wn) = slope*((wc(wn)*1.e-3_dp)-0.6_dp)+gaer600(k) if(gaer(k,wn) .lt. 0.5_dp) gaer(k,wn) = 0.5_dp if(gaer(k,wn) .ge. 1.0_dp) gaer(k,wn) = 1.0_dp endif enddo ! k enddo ! wn end subroutine aer_wrf2ftuv !-------------------------------------------------------------------- subroutine setaer( nzlev, z, airden, rh, acb1, & acb2, aoc1, aoc2, aant, aso4, & asal, & op_cb1, om_cb1, g_cb1, & op_cb2, om_cb2, g_cb2, & op_oc1, om_oc1, g_oc1, & op_oc2, om_oc2, g_oc2, & op_ant, om_ant, g_ant, & op_so4, om_so4, g_so4, & op_sal, om_sal, g_sal ) use module_wave_data, only : nw, wc !----------------------------------------------------------------------------- != PURPOSE: ! != Cal aer opt depth !----------------------------------------------------------------------------- ! != PARAMETERS: != NZLEV - INTEGER, number of specified altitude levels in the working (I) != grid != Z - real(dp), specified altitude working grid (km) (I) != axxx - aerosol mix ratio (I) ! != op - real(dp), optical depth (absorption) (O) != om - real(dp), single albedo != g - asysmetric factor (O) ! !----------------------------------------------------------------------------- ! ! VERTICAL DOMAIN is from bottom(1) to TOP (TOP=nzlev) ! CCM from top(1) to bottom(nzlev) !----------------------------------------------------------------------------- implicit none !----------------------------------------------------------------------------- ! ... Dummy arguments !----------------------------------------------------------------------------- integer, intent(in) :: nzlev real(dp), intent(in) :: z(nzlev) real(dp), intent(in) :: airden(nzlev) real(dp), intent(in) :: rh(nzlev) real(dp), intent(in) :: acb1(nzlev) real(dp), intent(in) :: acb2(nzlev) real(dp), intent(in) :: aoc1(nzlev) real(dp), intent(in) :: aoc2(nzlev) real(dp), intent(in) :: aant(nzlev) real(dp), intent(in) :: aso4(nzlev) real(dp), intent(in) :: asal(nzlev) real(dp), intent(out) :: op_cb1(nzlev-1,nw-1), om_cb1(nzlev-1,nw-1), g_cb1(nzlev-1,nw-1) real(dp), intent(out) :: op_cb2(nzlev-1,nw-1), om_cb2(nzlev-1,nw-1), g_cb2(nzlev-1,nw-1) real(dp), intent(out) :: op_oc1(nzlev-1,nw-1), om_oc1(nzlev-1,nw-1), g_oc1(nzlev-1,nw-1) real(dp), intent(out) :: op_oc2(nzlev-1,nw-1), om_oc2(nzlev-1,nw-1), g_oc2(nzlev-1,nw-1) real(dp), intent(out) :: op_ant(nzlev-1,nw-1), om_ant(nzlev-1,nw-1), g_ant(nzlev-1,nw-1) real(dp), intent(out) :: op_so4(nzlev-1,nw-1), om_so4(nzlev-1,nw-1), g_so4(nzlev-1,nw-1) real(dp), intent(out) :: op_sal(nzlev-1,nw-1), om_sal(nzlev-1,nw-1), g_sal(nzlev-1,nw-1) !----------------------------------------------------------------------------- ! ... local variables !----------------------------------------------------------------------------- integer :: k, wn, nz real(dp) :: wcen real(dp) :: dz(nzlev) real(dp) :: sig_oc2(nzlev-1), sig_cb2(nzlev-1) real(dp) :: opt_oc1(nzlev-1), opt_oc2(nzlev-1), opt_cb1(nzlev-1), opt_cb2(nzlev-1) real(dp) :: opt_ant(nzlev-1), opt_so4(nzlev-1), opt_sal(nzlev-1) real(dp) :: rw(nzlev-1), num(nzlev-1), wght(nzlev-1) !----------------------------------------------------------------------------- ! ... assign the constants !----------------------------------------------------------------------------- real(dp), parameter :: rm_cb1 = 0.035_dp real(dp), parameter :: rm_cb2 = 0.035_dp real(dp), parameter :: rm_oc1 = 0.10_dp real(dp), parameter :: rm_oc2 = 0.10_dp real(dp), parameter :: rm_ant = 0.15_dp real(dp), parameter :: rm_so4 = 0.24_dp real(dp), parameter :: rm_sal = 1.30_dp real(dp), parameter :: de_cb1 = 1.0_dp real(dp), parameter :: de_cb2 = 1.0_dp real(dp), parameter :: de_oc1 = 1.5_dp real(dp), parameter :: de_oc2 = 1.5_dp real(dp), parameter :: de_ant = 1.7_dp real(dp), parameter :: de_so4 = 1.7_dp real(dp), parameter :: de_sal = 2.2_dp real(dp), parameter :: ml_cb1 = 12._dp real(dp), parameter :: ml_cb2 = 12._dp real(dp), parameter :: ml_oc1 = 48._dp real(dp), parameter :: ml_oc2 = 48._dp real(dp), parameter :: ml_ant = 90._dp real(dp), parameter :: ml_so4 = 96._dp real(dp), parameter :: ml_sal = 22._dp real(dp), parameter :: qx_cb1 = 2.138_dp real(dp), parameter :: qx_cb2 = 2.138_dp real(dp), parameter :: qw_cb1 = 0.510_dp real(dp), parameter :: qw_cb2 = 0.510_dp real(dp), parameter :: sx_cb2 = 0.120_dp real(dp), parameter :: sw_cb2 = 1.000_dp real(dp), parameter :: qx_oc1 = 0.675_dp real(dp), parameter :: qx_oc2 = 0.675_dp real(dp), parameter :: qw_oc1 = 1.118_dp real(dp), parameter :: qw_oc2 = 1.118_dp real(dp), parameter :: sx_oc2 = 0.980_dp real(dp), parameter :: sw_oc2 = 1.000_dp real(dp), parameter :: qx_ant = 0.726_dp real(dp), parameter :: qw_ant = 1.770_dp real(dp), parameter :: qx_so4 = 2.830_dp real(dp), parameter :: qw_so4 = 3.605_dp real(dp), parameter :: qx_sal = 2.600_dp real(dp), parameter :: qw_sal = 2.575_dp real(dp), parameter :: pi = 3.1416_dp !----------------------------------------------------------------------------- ! ... hydroscropic growth !----------------------------------------------------------------------------- real(dp), parameter :: a1_cb2 = 1.000_dp real(dp), parameter :: a2_cb2 =-0.829_dp real(dp), parameter :: a3_cb2 = 1.350_dp real(dp), parameter :: a1_oc2 = 1.005_dp real(dp), parameter :: a2_oc2 =-0.0666_dp real(dp), parameter :: a3_oc2 = 0.7608_dp real(dp), parameter :: a1_ant = 1.000_dp real(dp), parameter :: a2_ant = 0.493_dp real(dp), parameter :: a3_ant = 0.527_dp real(dp), parameter :: a1_so4 = 1.000_dp real(dp), parameter :: a2_so4 = 0.493_dp real(dp), parameter :: a3_so4 = 0.528_dp real(dp), parameter :: a1_sal = 1.011_dp real(dp), parameter :: a2_sal = 0.457_dp real(dp), parameter :: a3_sal = 1.102_dp nz = nzlev - 1 !----------------------------------------------------------------------------- ! ... calculate vertical interveral dz !----------------------------------------------------------------------------- dz(1:nz) = abs( z(2:nz+1) - z(1:nz) ) * km2cm * pi * 1.e-8_dp !----------------------------------------------------------------------------- ! -------- CB1 ------------- !----------------------------------------------------------------------------- rw(:nz) = rm_cb1 call cnum( acb1, airden, de_cb1, ml_cb1, rm_cb1, & rw, num, wght, nzlev ) opt_cb1(:nz) = qx_cb1*wght(:nz) + (1._dp - wght(:nz))*qw_cb1 opt_cb1(:nz) = opt_cb1(:nz)*num(:nz)*rw(:nz)*rw(:nz)*dz(:nz) !----------------------------------------------------------------------------- ! -------- CB2 ------------- !----------------------------------------------------------------------------- rw(:nz) = rm_cb2*(a1_cb2 + rh(:nz)*(a2_cb2 + a3_cb2*rh(:nz))) call cnum( acb2, airden, de_cb2, ml_cb2, rm_cb2, & rw, num, wght, nzlev ) opt_cb2(:nz) = qx_cb2*wght(:nz) + (1._dp - wght(:nz))*qw_cb2 opt_cb2(:nz) = opt_cb2(:nz)*num(:nz)*rw(:nz)*rw(:nz)*dz(:nz) sig_cb2(:nz) = sx_cb2*wght(:nz) + (1._dp - wght(:nz))*sw_cb2 !----------------------------------------------------------------------------- ! -------- OC1 ------------- !----------------------------------------------------------------------------- rw(:nz) = rm_oc1 call cnum( aoc1, airden, de_oc1, ml_oc1, rm_oc1, & rw, num, wght, nzlev ) opt_oc1(:nz) = qx_oc1*wght(:nz) + (1._dp - wght(:nz))*qw_oc1 opt_oc1(:nz) = opt_oc1(:nz)*num(:nz)*rw(:nz)*rw(:nz)*dz(:nz) !----------------------------------------------------------------------------- ! -------- OC2 ------------- !----------------------------------------------------------------------------- rw(:nz) = rm_oc2*(a1_oc2 + rh(:nz)*(a2_oc2 + a3_oc2*rh(:nz))) call cnum( aoc2, airden, de_oc2, ml_oc2, rm_oc2, & rw, num, wght, nzlev ) opt_oc2(:nz) = qx_oc2*wght(:nz) + (1._dp - wght(:nz))*qw_oc2 opt_oc2(:nz) = opt_oc2(:nz)*num(:nz)*rw(:nz)*rw(:nz)*dz(:nz) sig_oc2(:nz) = sx_oc2*wght(:nz) + (1._dp - wght(:nz))*sw_oc2 !----------------------------------------------------------------------------- ! -------- ANT-NH4NO3 ------------- !----------------------------------------------------------------------------- rw(:nz) = rm_ant*(a1_ant + rh(:nz)*(a2_ant + a3_ant*rh(:nz))) call cnum( aant, airden, de_ant, ml_ant, rm_ant, & rw, num, wght, nzlev ) opt_ant(:nz) = qx_ant*wght(:nz) + (1._dp - wght(:nz))*qw_ant opt_ant(:nz) = opt_ant(:nz)*num(:nz)*rw(:nz)*rw(:nz)*dz(:nz) !----------------------------------------------------------------------------- ! -------- SO4 ------------- !----------------------------------------------------------------------------- rw(:nz) = rm_so4*(a1_so4 + rh(:nz)*(a2_so4 + a3_so4*rh(:nz))) call cnum( aso4, airden, de_so4, ml_so4, rm_so4, & rw, num, wght, nzlev ) opt_so4(:nz) = qx_so4*wght(:nz) + (1._dp - wght(:nz))*qw_so4 opt_so4(:nz) = opt_so4(:nz)*num(:nz)*rw(:nz)*rw(:nz)*dz(:nz) !----------------------------------------------------------------------------- ! -------- SALT ------------- !----------------------------------------------------------------------------- rw(:nz) = rm_sal*(a1_sal + rh(:nz)*(a2_sal + a3_sal*rh(:nz))) call cnum( asal, airden, de_sal, ml_sal, rm_sal, & rw, num, wght, nzlev ) opt_sal(:nz) = qx_sal*wght(:nz) + (1._dp - wght(:nz))*qw_sal opt_sal(:nz) = opt_sal(:nz)*num(:nz)*rw(:nz)*rw(:nz)*dz(:nz) WAVE_LENGTH_LOOP : & do wn = 1,nw-1 wcen = wc(wn) op_cb1(:nz,wn) = opt_cb1(:nz)*op_cb1w(wn) om_cb1(:nz,wn) = om_cb1w(wn) g_cb1(:nz,wn) = g_cb1w(wn) op_cb2(:nz,wn) = opt_cb2(:nz)*op_cb2w(wn) om_cb2(:nz,wn) = sig_cb2(:nz)*om_cb2w(wn) g_cb2(:nz,wn) = g_cb2w(wn) op_oc1(:nz,wn) = opt_oc1(:nz)*op_oc1w(wn) om_oc1(:nz,wn) = 0.98_dp g_oc1(:nz,wn) = g_oc1w(wn) op_oc2(:nz,wn) = opt_oc2(:nz)*op_oc2w(wn) om_oc2(:nz,wn) = sig_oc2(:nz) g_oc2(:nz,wn) = g_oc2w(wn) op_ant(:nz,wn) = opt_ant(:nz)*op_antw(wn) om_ant(:nz,wn) = 1.00_dp g_ant(:nz,wn) = g_antw(wn) op_so4(:nz,wn) = opt_so4(:nz)*op_so4w(wn) om_so4(:nz,wn) = 1.00_dp g_so4(:nz,wn) = 0.75_dp op_sal(:nz,wn) = opt_sal(:nz)*op_salw(wn) om_sal(:nz,wn) = 1.00_dp g_sal(:nz,wn) = 0.77_dp enddo WAVE_LENGTH_LOOP end subroutine setaer subroutine aer_init !----------------------------------------------------------------------------- ! ... initialize setaer !----------------------------------------------------------------------------- use module_wave_data, only : nw, wc implicit none !----------------------------------------------------------------------- ! Assign the wave-length dep values for aerososl !----------------------------------------------------------------------- real(dp), parameter :: d1_cb = 1.61892_dp real(dp), parameter :: d2_cb =-0.0006853_dp real(dp), parameter :: d3_cb =-1.07806e-6_dp real(dp), parameter :: o1_cb = 1.29185_dp real(dp), parameter :: o2_cb = 7.6863e-5_dp real(dp), parameter :: o3_cb =-1.28567e-6_dp real(dp), parameter :: g1_cb = 2.87326_dp real(dp), parameter :: g2_cb =-0.004482_dp real(dp), parameter :: g3_cb = 1.50361e-6_dp real(dp), parameter :: d1_oc = 10.5098_dp real(dp), parameter :: d2_oc =-0.0301455_dp real(dp), parameter :: d3_oc = 2.22837e-5_dp real(dp), parameter :: g1_oc = 1.67125_dp real(dp), parameter :: g2_oc =-0.0008341_dp real(dp), parameter :: g3_oc =-1.18801e-6_dp real(dp), parameter :: d1_ant = 9.44794_dp real(dp), parameter :: d2_ant =-0.0268029_dp real(dp), parameter :: d3_ant = 1.97106e-5_dp real(dp), parameter :: g1_ant = 1.34742_dp real(dp), parameter :: g2_ant =-1.85850e-5_dp real(dp), parameter :: g3_ant =-1.54297e-6_dp real(dp), parameter :: d1_so4 = 0.942351_dp real(dp), parameter :: d2_so4 = 0.00220168_dp real(dp), parameter :: d3_so4 =-3.9622e-6_dp real(dp), parameter :: d1_sal = 1.0_dp real(dp), parameter :: d2_sal = 0.0_dp real(dp), parameter :: d3_sal = 0.0_dp op_cb1w(:nw-1) = d1_cb + wc(:nw-1)*(d2_cb + d3_cb*wc(:nw-1)) om_cb1w(:nw-1) = .12_dp*(o1_cb + wc(:nw-1)*(o2_cb + o3_cb*wc(:nw-1))) g_cb1w(:nw-1) = .12_dp*(g1_cb + wc(:nw-1)*(g2_cb + g3_cb*wc(:nw-1))) op_cb2w(:nw-1) = op_cb1w(:nw-1) om_cb2w(:nw-1) = om_cb1w(:nw-1) g_cb2w(:nw-1) = g_cb1w(:nw-1) op_oc1w(:nw-1) = d1_oc + wc(:nw-1)*(d2_oc + d3_oc*wc(:nw-1)) g_oc1w(:nw-1) = 0.50_dp*(g1_oc + wc(:nw-1)*(g2_oc + g3_oc*wc(:nw-1))) op_oc2w(:nw-1) = op_oc1w(:nw-1) ! om_oc2w(:nw-1) = sig_oc2 g_oc2w(:nw-1) = g_oc1w(:nw-1) op_antw(:nw-1) = d1_ant + wc(:nw-1)*(d2_ant + d3_ant*wc(:nw-1)) g_antw(:nw-1) = 0.63_dp*(g1_ant + wc(:nw-1)*(g2_ant + g3_ant*wc(:nw-1))) op_so4w(:nw-1) = d1_so4 + wc(:nw-1)*(d2_so4 + d3_so4*wc(:nw-1)) op_salw(:nw-1) = d1_sal + wc(:nw-1)*(d2_sal + d3_sal*wc(:nw-1)) end subroutine aer_init subroutine cnum( mix, denx, cden, molw, rm, & rw, num, weight, nzlev ) !----------------------------------------------------------------------------- ! PURPOSE: ! ! Cal aer num and percent weigh !----------------------------------------------------------------------------- implicit none !----------------------------------------------------------------------------- ! ... Dummy arguments !----------------------------------------------------------------------------- integer, intent(in) :: nzlev real(dp), intent(in) :: molw, rm real(dp), intent(in) :: cden real(dp), intent(in) :: rw(nzlev-1) real(dp), intent(in) :: mix(nzlev-1), denx(nzlev-1) real(dp), intent(out) :: num(nzlev-1), weight(nzlev-1) !----------------------------------------------------------------------------- ! ... Local variables !----------------------------------------------------------------------------- REAL(dp), PARAMETER :: PI = 3.14159265358979324 REAL(dp), PARAMETER :: vol_fac = 4._dp * PI / 3._dp integer :: k real(dp) :: rm_cubic, factor real(dp) :: mas(nzlev-1), masw(nzlev-1) rm_cubic = rm**3 factor = 1._dp/(vol_fac * cden * 1.e-12_dp * rm_cubic) !----------------------------------------------------------------------------- ! Convert volumn mix ratio to mass !----------------------------------------------------------------------------- where( mix(:nzlev-1) >= 1.e-15_dp ) !=============================================================================== ! XUEXI ! UNLIKE MOZART mix = mix ratio ! mix = ug/m3 ! 1e-12 for ug/m3 --> gram/cm3 !=============================================================================== mas(:nzlev-1) = mix(:nzlev-1) * 1.e-12_dp ! 1e-12 for ug/m3 --> gram/cm3 num(:nzlev-1) = mas(:nzlev-1) * factor elsewhere num(:nzlev-1) = 0._dp weight(:nzlev-1) = 0._dp endwhere do k = 1,nzlev-1 if( abs(rw(k) - rm) >= 1.e-6_dp ) then masw(k) = vol_fac*(rw(k)**3 - rm_cubic)*1.e-12_dp * num(k) ! water mass gram/cm3 else masw(k) = 0._dp end if end do where( mix(:nzlev-1) >= 1.e-15_dp ) weight(:nzlev-1) = mas(:nzlev-1)/( mas(:nzlev-1) + masw(:nzlev-1) ) ! weight percentage (%*0.01) weight(:nzlev-1) = max( 0._dp, min( 1._dp,weight(:nzlev-1) ) ) endwhere end subroutine cnum subroutine setair( nz, z, airlev, dtrl, & cz, o2top ) use module_wave_data, only : nw, wc !----------------------------------------------------------------------------- ! purpose: ! set up an altitude profile of air molecules. subroutine includes a ! shape-conserving scaling method that allows scaling of the entire ! profile to a given sea-level pressure. !----------------------------------------------------------------------------- ! parameters: ! pmbnew - real(dp), sea-level pressure (mb) to which profile should be ! scaled. if pmbnew < 0, no scaling is done ! nz - integer, number of specified altitude levels in the working (i) ! grid ! z - real(dp), specified altitude working grid (km) (i) ! nw - integer, number of specified intervals + 1 in working (i) ! wavelength grid ! wl - real(dp), vector of lower limits of wavelength intervals in (i) ! working wavelength grid ! airlev - real(dp), air density (molec/cc) at each specified altitude (o) ! dtrl - real(dp), rayleigh optical depth at each specified altitude (o) ! and each specified wavelength ! cz - real(dp), number of air molecules per cm^2 at each specified (o) ! altitude layer !----------------------------------------------------------------------------- implicit none !----------------------------------------------------------------------------- ! ... dummy arguments !----------------------------------------------------------------------------- integer, intent(in) :: nz real(dp), intent(in) :: o2top real(dp), intent(in) :: z(nz) real(dp), intent(in) :: airlev(nz) !----------------------------------------------------------------------------- ! ... air density (molec cm-3) at each grid level ! rayleigh optical depths !----------------------------------------------------------------------------- real(dp), intent(out) :: dtrl(nz-1,nw-1) real(dp), intent(out) :: cz(nz) !----------------------------------------------------------------------------- ! ... local variables !----------------------------------------------------------------------------- integer :: i, ip1, iw real(dp) :: hscale real(dp) :: srayl real(dp) :: deltaz real(dp) :: wmicrn, xx !----------------------------------------------------- ! ... compute column increments (logarithmic integrals) !----------------------------------------------------- do i = 1,nz-1 ip1 = i + 1 deltaz = km2cm * (z(ip1) - z(i)) cz(i) = (airlev(ip1) - airlev(i)) / log(airlev(ip1)/airlev(i)) * deltaz end do !----------------------------------------------------- ! ... include exponential tail integral from infinity to 50 km, ! fold tail integral into top layer ! specify scale height near top of data. (scale height at 40km????) !----------------------------------------------------- ! hscale = 8.05e5 cz(nz) = o2top/.2095_dp !----------------------------------------------------- ! ... compute rayleigh cross sections and depths: !----------------------------------------------------- do iw = 1,nw-1 !----------------------------------------------------- ! ... rayleigh scattering cross section from wmo 1985 (originally from ! nicolet, m., on the molecular scattering in the terrestrial atmosphere: ! an empirical formula for its calculation in the homoshpere, planet. ! space sci., 32, 1467-1468, 1984. !----------------------------------------------------- wmicrn = 1.e-3_dp*wc(iw) if( wmicrn <= .55_dp ) then xx = 3.6772_dp + 0.389_dp*wmicrn + 0.09426_dp/wmicrn else xx = 4.04_dp end if srayl = 4.02e-28_dp/(wmicrn)**xx dtrl(1:nz-2,iw) = cz(1:nz-2)*srayl dtrl(nz-1,iw) = (cz(nz-1) + cz(nz))*srayl end do end subroutine setair subroutine zenith( lat, long, idate, ut, zen ) !----------------------------------------------------------------------------- ! purpose: ! calculate solar zenith angle for a given time and location. ! calculation is based on equations given in: paltridge and platt, radia- ! tive processes in meteorology and climatology, elsevier, pp. 62,63, 1976. ! fourier coefficients originally from: spencer, j.w., 1971, fourier ! series representation of the position of the sun, search, 2:172. ! note: this approximate program does not account fro changes from year ! to year. !----------------------------------------------------------------------------- ! parameters: ! lat - real(dp), latitude of location (degrees) (i) ! long - real(dp), longitude of location (degrees) (i) ! idate - integer, date in the form yymmdd (i) ! ut - real(dp), local time in decimal ut (e.g., 16.25 means 15 minutes (i) ! after 4 pm) ! zen - real(dp), solar zenith angle (degrees) (o) !----------------------------------------------------------------------------- implicit none !----------------------------------------------------------------------------- ! ... dummy arguments !----------------------------------------------------------------------------- integer, intent(in) :: idate real(dp), intent(in) :: lat,long real(dp), intent(in) :: ut real(dp), intent(out) :: zen !----------------------------------------------------------------------------- ! ... local variables !----------------------------------------------------------------------------- real(dp), parameter :: pi = 3.1415926 real(dp), parameter :: d2r = pi/180.0 real(dp), parameter :: r2d = 1.0/d2r integer :: i integer :: iiyear, imth, iday, ijd integer :: imn(12) = (/ 31,28,31,30,31,30,31,31,30,31,30,31 /) real(dp) :: lbut,lzut real(dp) :: rlt real(dp) :: d, tz, rdecl, eqr, eqh, zpt real(dp) :: csz, zr, caz, raz real(dp) :: sintz, costz, sin2tz, cos2tz, sin3tz, cos3tz !----------------------------------------------------------------------------- ! ... convert to radians !----------------------------------------------------------------------------- rlt = lat*d2r !----------------------------------------------------------------------------- ! ... parse date !----------------------------------------------------------------------------- iiyear = idate/10000 imth = (idate - iiyear*10000)/100 iday = idate - iiyear*10000 - imth*100 !----------------------------------------------------------------------------- ! ... identify and correct leap years !----------------------------------------------------------------------------- if( mod(iiyear,4) == 0 ) then imn(2) = 29 else imn(2) = 28 end if !----------------------------------------------------------------------------- ! ... compute current (julian) day of year ijd = 1 to 365 !----------------------------------------------------------------------------- ijd = 0 do i = 1,imth-1 ijd = ijd + imn(i) end do ijd = ijd + iday !----------------------------------------------------------------------------- ! ... calculate decimal julian day from start of year: !----------------------------------------------------------------------------- d = real(ijd-1,kind=dp) + ut/24._dp !----------------------------------------------------------------------------- ! ... equation 3.8 for "day-angle" !----------------------------------------------------------------------------- tz = 2._dp*pi*d/365._dp !----------------------------------------------------------------------------- ! ... calculate sine and cosine from addition theoremes for ! better performance; the computation of sin2tz, ! sin3tz, cos2tz and cos3tz is about 5-6 times faster ! than the evaluation of the intrinsic functions ! ! it is sin(x+y) = sin(x)*cos(y)+cos(x)*sin(y) ! and cos(x+y) = cos(x)*cos(y)-sin(x)*sin(y) ! ! sintz = sin(tz) costz = cos(tz) ! sin2tz = sin(2.*tz) cos2tz = sin(2.*tz) ! sin3tz = sin(3.*tz) cos3tz = cos(3.*tz) !----------------------------------------------------------------------------- sintz = sin( tz ) costz = cos( tz ) sin2tz = 2._dp*sintz*costz cos2tz = costz*costz - sintz*sintz sin3tz = sintz*cos2tz + costz*sin2tz cos3tz = costz*cos2tz - sintz*sin2tz !----------------------------------------------------------------------------- ! ... equation 3.7 for declination in radians !----------------------------------------------------------------------------- rdecl = 0.006918_dp - 0.399912_dp*costz + 0.070257_dp*sintz & - 0.006758_dp*cos2tz + 0.000907_dp*sin2tz & - 0.002697_dp*cos3tz + 0.001480_dp*sin3tz !----------------------------------------------------------------------------- ! ... equation 3.11 for equation of time in radians !----------------------------------------------------------------------------- eqr = 0.000075_dp + 0.001868_dp*costz - 0.032077_dp*sintz & - 0.014615_dp*cos2tz - 0.040849_dp*sin2tz !----------------------------------------------------------------------------- ! ... convert equation of time to hours !----------------------------------------------------------------------------- eqh = eqr*24._dp/(2._dp*pi) !----------------------------------------------------------------------------- ! ... calculate local hour angle (hours): !----------------------------------------------------------------------------- lbut = 12._dp - eqh - long*24._dp/360._dp !----------------------------------------------------------------------------- ! ... convert to angle from ut !----------------------------------------------------------------------------- lzut = 15._dp*(ut - lbut) zpt = lzut*d2r !----------------------------------------------------------------------------- ! ... equation 2.4 for cosine of zenith angle !----------------------------------------------------------------------------- csz = sin(rlt)*sin(rdecl) + cos(rlt)*cos(rdecl)*cos(zpt) zr = acos(csz) zen = zr*r2d end subroutine zenith subroutine calc_zenith( lat, long, ijd, gmt, azimuth, zenith ) !------------------------------------------------------------------- ! this subroutine calculates solar zenith and azimuth angles for a ! part time and location. must specify: ! input: ! lat - latitude in decimal degrees ! long - longitude in decimal degrees ! gmt - greenwich mean time - decimal military eg. ! 22.75 = 45 min after ten pm gmt ! output: ! zenith ! azimuth ! .. Scalar Arguments .. !------------------------------------------------------------------- integer, intent(in) :: ijd real(dp), intent(in) :: gmt, lat, long real(dp), intent(out) :: azimuth, zenith !------------------------------------------------------------------- ! .. Local variables !------------------------------------------------------------------- real(dp), parameter :: d2r = 3.1415926_dp/180.0_dp real(dp), parameter :: r2d = 1.0_dp/d2r integer :: jd real(dp) :: caz, csz, cw, d, decl, ec, epsi, eqt, eyt, feqt, feqt1, & feqt2, feqt3, feqt4, feqt5, feqt6, feqt7, lbgmt, lzgmt, ml, pepsi, & pi, ra, raz, rdecl, reqt, rlt, rml, rphi, rra, ssw, sw, tab, w, wr, & yt, zpt, zr !------------------------------------------------------------------- ! .. Intrinsic Functions .. !------------------------------------------------------------------- intrinsic acos, atan, cos, min, sin, tan !------------------------------------------------------------------- ! convert to radians !------------------------------------------------------------------- rlt = lat*d2r rphi = long*d2r jd = ijd d = jd + gmt/24.0_dp !------------------------------------------------------------------- ! calc geom mean longitude !------------------------------------------------------------------- ml = 279.2801988_dp + d*(.9856473354_dp + 2.267E-13_dp*d) rml = ml*d2r !------------------------------------------------------------------- ! calc equation of time in sec ! w = mean long of perigee ! e = eccentricity ! epsi = mean obliquity of ecliptic !------------------------------------------------------------------- w = 282.4932328_dp + d*(4.70684E-5_dp + 3.39E-13_dp*d) wr = w*d2r ec = 1.6720041E-2_dp - d*(1.1444E-9_dp + 9.4E-17_dp*d) epsi = 23.44266511_dp - d*(3.5626E-7_dp + 1.23E-15_dp*d) pepsi = epsi*d2r yt = (tan(pepsi/2.0_dp))**2 cw = cos(wr) sw = sin(wr) ssw = sin(2.0_dp*wr) eyt = 2._dp*ec*yt feqt1 = -sin(rml)*cw*(eyt + 2._dp*ec) feqt2 = cos(rml)*sw*(2._dp*ec - eyt) feqt3 = sin(2._dp*rml)*(yt - (5._dp*ec**2/4._dp)*(cw**2 - sw**2)) feqt4 = cos(2._dp*rml)*(5._dp*ec**2*ssw/4._dp) feqt5 = sin(3._dp*rml)*(eyt*cw) feqt6 = -cos(3._dp*rml)*(eyt*sw) feqt7 = -sin(4._dp*rml)*(.5_dp*yt**2) feqt = feqt1 + feqt2 + feqt3 + feqt4 + feqt5 + feqt6 + feqt7 eqt = feqt*13751.0_dp !------------------------------------------------------------------- ! convert eq of time from sec to deg !------------------------------------------------------------------- reqt = eqt/240._dp !------------------------------------------------------------------- ! calc right ascension in rads !------------------------------------------------------------------- ra = ml - reqt rra = ra*d2r !------------------------------------------------------------------- ! calc declination in rads, deg !------------------------------------------------------------------- tab = 0.43360_dp*sin(rra) rdecl = atan(tab) decl = rdecl/d2r !------------------------------------------------------------------- ! calc local hour angle !------------------------------------------------------------------- lbgmt = 12.0_dp - eqt/3600._dp + long*24._dp/360._dp lzgmt = 15.0_dp*(gmt-lbgmt) zpt = lzgmt*d2r csz = sin(rlt)*sin(rdecl) + cos(rlt)*cos(rdecl)*cos(zpt) if(csz.gt.1)print *,'calczen,csz ',csz csz = min( 1._dp,csz ) ! zr = acos(csz) ! zenith = zr/d2r zr = acos(csz) zenith = zr/d2r !------------------------------------------------------------------- ! calc local solar azimuth !------------------------------------------------------------------- caz = (sin(rdecl)-sin(rlt)*cos(zr))/(cos(rlt)*sin(zr)) if(caz < -0.999999_dp)then azimuth=180._dp elseif(caz > 0.999999_dp)then azimuth=0._dp else raz = acos(caz) azimuth = raz/d2r endif ! caz = min(1.,(sin(rdecl)-sin(rlt)*cos(zr))/(cos(rlt)*sin(zr))) ! if(caz.lt.-1)print *,calczen ,caz ! caz = max(-1.,caz) ! raz = acos(caz) ! azimuth = raz/d2r if( lzgmt > 0._dp ) azimuth = azimuth + (2._dp*(180._dp - azimuth)) end subroutine calc_zenith subroutine sundis( julday, esrm2 ) !----------------------------------------------------------------------------- ! purpose: ! calculate earth-sun distance variation for a given date. based on ! fourier coefficients originally from: spencer, j.w., 1971, fourier ! series representation of the position of the sun, search, 2:172 !----------------------------------------------------------------------------- ! parameters: ! idate - integer, specification of the date, from yymmdd (i) ! esrm2 - real(dp), variation of the earth-sun distance (o) ! esrm2 = (average e/s dist)^2 / (e/s dist on day idate)^2 !----------------------------------------------------------------------------- implicit none !----------------------------------------------------------------------------- ! ... dummy arguments !----------------------------------------------------------------------------- integer, intent(in) :: julday real(dp), intent(out) :: esrm2 !----------------------------------------------------------------------------- ! ... local variables !----------------------------------------------------------------------------- real(dp), parameter :: pi = 3.1415926_dp real(dp) :: dayn, thet0 real(dp) :: sinth, costh, sin2th, cos2th !----------------------------------------------------------------------------- ! ... parse date to find day number (julian day) !----------------------------------------------------------------------------- dayn = real(julday - 1,kind=dp) + .5_dp !----------------------------------------------------------------------------- ! ... define angular day number and compute esrm2: !----------------------------------------------------------------------------- thet0 = 2._dp*pi*dayn/365._dp !----------------------------------------------------------------------------- ! ... calculate sin(2*thet0), cos(2*thet0) from ! addition theorems for trig functions for better ! performance; the computation of sin2th, cos2th ! is about 5-6 times faster than the evaluation ! of the intrinsic functions sin and cos !----------------------------------------------------------------------------- sinth = sin( thet0 ) costh = cos( thet0 ) sin2th = 2._dp*sinth*costh cos2th = costh*costh - sinth*sinth esrm2 = 1.000110_dp + .034221_dp*costh + .001280_dp*sinth + .000719_dp*cos2th + .000077_dp*sin2th end subroutine sundis subroutine spheres( nz, z, zen, dsdh, nid ) !----------------------------------------------------------------------------- ! purpose: ! calculate slant path over vertical depth ds/dh in spherical geometry. ! calculation is based on: a.dahlback, and k.stamnes, a new spheric model ! for computing the radiation field available for photolysis and heating ! at twilight, planet.space sci., v39, n5, pp. 671-683, 1991 (appendix b) !----------------------------------------------------------------------------- ! parameters: ! nz - integer, number of specified altitude levels in the working (i) ! grid ! z - real(dp), specified altitude working grid (km) (i) ! zen - real(dp), solar zenith angle (degrees) (i) ! dsdh - real(dp), slant path of direct beam through each layer crossed (o) ! when travelling from the top of the atmosphere to layer i; ! dsdh(i,j), i = 0..nz-1, j = 1..nz-1 ! nid - integer, number of layers crossed by the direct beam when (o) ! travelling from the top of the atmosphere to layer i; ! nid(i), i = 0..nz-1 !----------------------------------------------------------------------------- implicit none !----------------------------------------------------------------------------- ! ... dummy arguments !----------------------------------------------------------------------------- integer, intent(in) :: nz integer, intent(out) :: nid(0:nz-1) real(dp), intent(in) :: zen real(dp), intent(in) :: z(nz) real(dp), intent(out) :: dsdh(0:nz-1,nz-1) !----------------------------------------------------------------------------- ! ... local variables !----------------------------------------------------------------------------- real(dp), parameter :: radius = 6.371e+3_dp ! radius earth (km) real(dp), parameter :: d2r = 3.1415926_dp/180.0_dp integer :: i, j, k integer :: id real(dp) :: re, ze(nz) real(dp) :: zd(0:nz-1) real(dp) :: zenrad, rpsinz, rj, rjp1, dsj, dhj, ga, gb, sm zenrad = zen*d2r !----------------------------------------------------------------------------- ! ... include the elevation above sea level to the radius of the earth: !----------------------------------------------------------------------------- re = radius + z(1) !----------------------------------------------------------------------------- ! ... correspondingly z changed to the elevation above earth surface: !----------------------------------------------------------------------------- ze(1:nz) = z(1:nz) - z(1) !----------------------------------------------------------------------------- ! ... inverse coordinate of z !----------------------------------------------------------------------------- zd(0) = ze(nz) do k = 1,nz-1 zd(k) = ze(nz - k) end do !----------------------------------------------------------------------------- ! ... initialize dsdh(i,j), nid(i) !----------------------------------------------------------------------------- dsdh(0:nz-1,1:nz-1) = 0._dp nid(0:nz-1) = 0 !----------------------------------------------------------------------------- ! ... calculate ds/dh of every layer !----------------------------------------------------------------------------- do i = 0,nz-1 rpsinz = (re + zd(i)) * sin(zenrad) if( zen > 90._dp .and. rpsinz < re ) then nid(i) = -1 else !----------------------------------------------------------------------------- ! ... find index of layer in which the screening height lies !----------------------------------------------------------------------------- id = i if( zen > 90._dp ) then do j = 1,nz-1 if( (rpsinz < ( zd(j-1) + re ) ) .and. (rpsinz >= ( zd(j) + re )) ) then id = j end if end do end if do j = 1,id if( j == id .and. id == i .and. zen > 90._dp ) then sm = -1._dp else sm = 1._dp end if rj = re + zd(j-1) rjp1 = re + zd(j) dhj = zd(j-1) - zd(j) ga = max( 0._dp,rj*rj - rpsinz*rpsinz ) gb = max( 0._dp,rjp1*rjp1 - rpsinz*rpsinz ) if( id > i .and. j == id ) then dsj = sqrt( ga ) else dsj = sqrt( ga ) - sm*sqrt( gb ) end if dsdh(i,j) = dsj / dhj end do nid(i) = id end if end do end subroutine spheres subroutine setz( nz, cz, tlev, c, ndx, adjcoe ) use module_wave_data, only : tuv_jmax !----------------------------------------------------------------------------- ! adjcoe - real(dp), coross section adjust coefficients !----------------------------------------------------------------------------- implicit none !----------------------------------------------------------------------------- ! ... dummy arguments !----------------------------------------------------------------------------- integer, intent(in) :: nz integer, intent(in) :: ndx real(dp), intent(in) :: cz(nz) real(dp), intent(in) :: tlev(nz) real(dp), intent(in) :: c(5,tuv_jmax) real(dp), intent(inout) :: adjcoe(nz,tuv_jmax) !----------------------------------------------------------------------------- ! ... local variables !----------------------------------------------------------------------------- integer :: k, m real(dp) :: tt, adjin real(dp) :: c0, c1, c2 real(dp) :: xz(nz) !----------------------------------------------------------------------------- ! 1 o2 + hv -> o + o ! 2 o3 -> o2 + o(1d) ! 3 o3 -> o2 + o(3p) ! 4 no2 -> no + o(3p) ! 5 no3 -> no + o2 ! 6 no3 -> no2 + o(3p) ! 7 n2o5 -> no3 + no + o(3p) ! 8 n2o5 -> no3 + no2 ! 9 n2o + hv -> n2 + o(1d) ! 10 ho2 + hv -> oh + o ! 11 h2o2 -> 2 oh ! 12 hno2 -> oh + no ! 13 hno3 -> oh + no2 ! 14 hno4 -> ho2 + no2 ! 15 ch2o -> h + hco ! 16 ch2o -> h2 + co ! 17 ch3cho -> ch3 + hco ! 18 ch3cho -> ch4 + co ! 19 ch3cho -> ch3co + h ! 20 c2h5cho -> c2h5 + hco ! 21 chocho -> products ! 22 ch3cocho -> products ! 23 ch3coch3 ! 24 ch3ooh -> ch3o + oh ! 25 ch3ono2 -> ch3o+no2 ! 26 pan + hv -> products ! 27 mvk + hv -> products ! 28 macr + hv -> products ! 29 hyac + hv -> products ! 30 glyald + hv -> products !----------------------------------------------------------------------------- xz(1:nz) = cz(1:nz)*1.e-18_dp do k = 1,tuv_jmax adjcoe(1:nz,k) = 1._dp end do tt = tlev(1)/281._dp do m = 1,tuv_jmax adjin = 1._dp if( m == 2 ) then !---------------------------------------------------------------------- ! ... temperature modification ! t0.9 (1.05) t0.95(1.025) t1.0(1.0) t1.15(1.02) t1.1(1.04) !---------------------------------------------------------------------- select case( ndx ) case( 20 ) c0 = 4.52372_dp ; c1 = -5.94317_dp ; c2 = 2.63156_dp case( 40 ) c0 = 4.99378_dp ; c1 = -7.92752_dp ; c2 = 3.94715_dp case( 60 ) c0 = .969867_dp ; c1 = -.841035_dp ; c2 = .878835_dp case( 80 ) c0 = 1.07801_dp ; c1 = -2.39580_dp ; c2 = 2.32632_dp end select adjin = c0 + tt*(c1 + c2*tt) else if( m == 11 ) then !---------------------------------------------------------------------- ! ... temperature modification ! t0.9 (1.05) t0.95(1.025) t1.0(1.0) t1.15(1.02) t1.1(1.04) !---------------------------------------------------------------------- select case( ndx ) case( 20 ) c0 = 2.43360_dp ; c1 = -3.61363_dp ; c2 = 2.19018_dp case( 40 ) c0 = 3.98265_dp ; c1 = -6.90516_dp ; c2 = 3.93602_dp case( 60 ) c0 = 3.49843_dp ; c1 = -5.98839_dp ; c2 = 3.50262_dp case( 80 ) c0 = 3.06312_dp ; c1 = -5.26281_dp ; c2 = 3.20980_dp end select adjin = c0 + tt*(c1 + c2*tt) end if call calcoe( nz, c(1,m), xz, tt, adjin, adjcoe(1,m) ) end do end subroutine setz subroutine setozo( nz, z, tlay, dto3, & to3, o3, airlev, o3top, o3toms ) use module_wave_data, only : nw, wl, xso3, s226, s263, s298 !----------------------------------------------------------------------------- ! purpose: ! set up an altitude profile of ozone, and corresponding absorption ! optical depths. subroutine includes a shape-conserving scaling method ! that allows scaling of the entire profile to a given overhead ozone ! column amount. !----------------------------------------------------------------------------- ! parameters: ! dobnew - real(dp), overhead ozone column amount (du) to which profile (i) ! should be scaled. if dobnew < 0, no scaling is done ! nz - integer, number of specified altitude levels in the working (i) ! grid ! z - real(dp), specified altitude working grid (km) ! nw - integer, number of specified intervals + 1 in working (i) ! wavelength grid ! wl - real(dp), vector of lower limits of wavelength intervals in (i) ! working wavelength grid ! xso3 - real(dp), molecular absoprtion cross section (cm^2) of o3 at (i) ! each specified wavelength (wmo value at 273) ! s226 - real(dp), molecular absoprtion cross section (cm^2) of o3 at (i) ! each specified wavelength (value from molina and molina at 226k) ! s263 - real(dp), molecular absoprtion cross section (cm^2) of o3 at (i) ! each specified wavelength (value from molina and molina at 263k) ! s298 - real(dp), molecular absoprtion cross section (cm^2) of o3 at (i) ! each specified wavelength (value from molina and molina at 298k) ! tlay - real(dp), temperature (k) at each specified altitude layer (i) ! dto3 - real(dp), optical depth due to ozone absorption at each (o) ! specified altitude at each specified wavelength ! to3 - real(dp), totol ozone column (o) !----------------------------------------------------------------------------- implicit none !----------------------------------------------------------------------------- ! ... dummy arguments !----------------------------------------------------------------------------- integer, intent(in) :: nz real(dp), intent(in) :: o3top real(dp), intent(in) :: o3toms real(dp), intent(in) :: z(nz) real(dp), intent(in) :: tlay(nz-1) real(dp), intent(in) :: airlev(nz) real(dp), intent(out) :: dto3(nz-1,nw-1) real(dp), intent(out) :: to3(nz) real(dp), intent(in) :: o3(nz) !----------------------------------------------------------------------------- ! ... local variables !----------------------------------------------------------------------------- integer :: k, iw, nd, ilat real(dp) :: so3 real(dp) :: scale real(dp) :: div1, div2 real(dp) :: o3den(nz) real(dp) :: cz(nz) !----------------------------------------------------------------------------- ! ... compute column increments !----------------------------------------------------------------------------- ! o3(1:nz) = o3(1:nz)*airlev(1:nz) ! cz(1:nz-1) = (o3(2:nz)) * 1.e5 * (z(2:nz) - z(1:nz-1)) o3den(1:nz) = o3(1:nz)*airlev(1:nz) cz(1:nz-1) = 0.5_dp*(o3den(2:nz) + o3den(1:nz-1))*km2cm*(z(2:nz) - z(1:nz-1)) to3(nz) = o3top do k = nz-1,1,-1 to3(k) = to3(k+1) + cz(k) end do ! Do not double count o3top cz(nz-1) = cz(nz-1) + o3top !----------------------------------------------------------------------------- ! ... scale o3 using toms data !----------------------------------------------------------------------------- if( o3toms > 0.0_dp ) then scale = o3toms/(to3(1)/2.687e16_dp) cz(1:nz) = cz(1:nz)*scale to3(1:nz) = to3(1:nz)*scale endif !----------------------------------------------------------------------------- ! ... calculate ozone optical depth for each layer, with temperature ! correction. output, dto3(kz,kw) !----------------------------------------------------------------------------- div1 = 1._dp / ( 263._dp - 226._dp) div2 = 1._dp / ( 298._dp - 263._dp) do iw = 1,nw-1 so3 = xso3(iw) do k = 1,nz-1 if( wl(iw) > 240.5_dp .and. wl(iw+1) < 350._dp ) then if( tlay(k) < 263._dp ) then so3 = s226(iw) + (s263(iw) - s226(iw)) * div1 * (tlay(k) - 226._dp) else so3 = s263(iw) + (s298(iw) - s263(iw)) * div2 * (tlay(k) - 263._dp) end if end if dto3(k,iw) = cz(k)*so3 end do end do end subroutine setozo subroutine set_o2_xsect( nz, z, cz, vcol, scol, dto2, xso2 ) use module_wave_data, only : nw, wl !----------------------------------------------------------------------------- ! purpose: ! compute equivalent optical depths for o2 absorption, parameterized in ! the sr bands and the lyman-alpha line. !----------------------------------------------------------------------------- ! parameters: ! nz - integer, number of specified altitude levels in the working (i) ! grid ! z - real(dp), specified altitude working grid (km) ! nw - integer, number of specified intervals + 1 in working ! wavelength grid ! wl - real(dp), vector of lower limits of wavelength intervals in ! working wavelength grid ! cz - real(dp), number of air molecules per cm^2 at each specified ! altitude layer ! zen - real(dp), solar zenith angle ! dto2 - real(dp), optical depth due to o2 absorption at each specified (o) ! vertical layer at each specified wavelength ! xso2 - real(dp), molecular absorption cross section in sr bands at (o) ! each specified altitude and wavelength. includes herzberg ! continuum. !----------------------------------------------------------------------------- ! edit history: ! 02/98 included lyman-alpha parameterization ! 03/97 fix dto2 problem at top level (nz) ! 02/97 changed offset for grid-end interpolation to relative number ! (x * (1 +- deltax)) ! 08/96 modified for early exit, no redundant read of data and smaller ! internal grid if possible; internal grid uses user grid points ! whenever possible ! 07/96 modified to work on internal grid and interpolate final values ! onto the user-defined grid !----------------------------------------------------------------------------- implicit none !----------------------------------------------------------------------------- ! ... dummy arguments !----------------------------------------------------------------------------- integer, intent(in) :: nz real(dp), intent(in) :: cz(nz) real(dp), intent(in) :: z(nz) real(dp), intent(in) :: vcol(nz) real(dp), intent(in) :: scol(nz) real(dp), intent(out) :: dto2(nz-1,nw-1) real(dp), intent(out) :: xso2(nw-1,nz) !----------------------------------------------------------------------------- ! ... local variables !----------------------------------------------------------------------------- integer :: iw, iz, igast real(dp) :: secchi(nz) !----------------------------------------------------------------------------- ! ... o2 optical depth and equivalent cross section on kockarts grid !----------------------------------------------------------------------------- real(dp) :: dto2k(nz-1,ngast-1) real(dp) :: xso2k(nz,ngast-1) !----------------------------------------------------------------------------- ! ... o2 optical depth and equivalent cross section in the lyman-alpha region !----------------------------------------------------------------------------- real(dp) :: dto2la(nz-1,nla-1) real(dp) :: xso2la(nz,nla-1) !----------------------------------------------------------------------------- ! ... temporary one-dimensional storage for optical depth and cross section values ! xxtmp - on internal grid ! xxuser - on user defined grid !----------------------------------------------------------------------------- real(dp), dimension(2*nwint) :: dttmp, xstmp ! CHANGED by Guohui real(dp), dimension(nwint) :: dtuser, xsuser ! CHANGED by Guohui real(dp) :: o2col(nz) real(dp) :: x, y real(dp) :: delo2 !----------------------------------------------------------------------------- ! ... check, whether user grid is in the o2 absorption band at all... ! if not, set cross section and optical depth values to zero and return !----------------------------------------------------------------------------- dto2(:nz-1,:nw-1) = 0._dp xso2(:nw-1,:nz) = 0._dp if( wl(1) > 243._dp ) then return end if !----------------------------------------------------------------------------- ! ... sec xhi or chapman calculation ! for zen > 95 degrees, use zen = 95. (this is only to compute effective o2 ! cross sections. still, better than setting dto2 = 0. as was done up to ! version 4.0) sm 1/2000 ! in future could replace with mu2(iz) (but mu2 is also wavelength-depenedent) ! or imporved chapman function !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- ! ... slant o2 column !----------------------------------------------------------------------------- o2col(1:nz) = 0.2095_dp * scol(1:nz) !----------------------------------------------------------------------------- ! ... effective secant of solar zenith angle. use 2.0 if no direct sun. ! for nz, use value at nz-1 !----------------------------------------------------------------------------- secchi(1:nz-1) = scol(1:nz-1)/vcol(1:nz-1) where( secchi(1:nz-1) == 0._dp ) secchi(1:nz-1) = 2._dp endwhere secchi(nz) = secchi(nz-1) !----------------------------------------------------------------------------- ! ... if necessary: ! kockarts parameterization of the sr bands, output values of o2 ! optical depth and o2 equivalent cross section are on his grid !----------------------------------------------------------------------------- if( wl(1) < wlgast(ngast) .and. wl(nw) > wlgast(1) ) then call schu( nz, o2col, secchi, dto2k, xso2k ) else dto2k(:,:) = 0._dp xso2k(:,:) = 0._dp end if !----------------------------------------------------------------------------- ! ... lyman-alpha parameterization, output values of o2 opticaldepth ! and o2 effective (equivalent) cross section !----------------------------------------------------------------------------- if( wl(1) <= wlla(nla) .and. wl(nw) >= wlla(1) ) then call lymana( nz, o2col, secchi, dto2la, xso2la ) else dto2la(:,:) = 0._dp xso2la(:,:) = 0._dp end if !----------------------------------------------------------------------------- ! ... loop through the altitude levels !----------------------------------------------------------------------------- do iz = 1,nz igast = 0 !----------------------------------------------------------------------------- ! ... loop through the internal wavelength grid !----------------------------------------------------------------------------- do iw = 1,nwint-1 !----------------------------------------------------------------------------- ! ... if outside kockarts grid and outside lyman-alpha, use the ! jpl/brasseur+solomon data, if inside ! kockarts grid, use the parameterized values from the call to schu, ! if inside lyman-alpha, use the paraemterized values from call to lymana !----------------------------------------------------------------------------- if( wlint(iw+1) <= wlgast(1) .or. wlint(iw) >= wlgast(ngast) ) then if( wlint(iw+1) <= wlla(1) .or. wlint(iw) >= wlla(nla) ) then xstmp(iw) = xso2int(iw) else xstmp(iw) = xso2la(iz,1) end if else igast = igast + 1 xstmp(iw) = xso2k(iz,igast) end if !----------------------------------------------------------------------------- ! ... compute the area in each bin (for correct interpolation purposes only!) !----------------------------------------------------------------------------- xstmp(iw) = xstmp(iw) * (wlint(iw+1) - wlint(iw)) end do !----------------------------------------------------------------------------- ! ... interpolate o2 cross section from the internal grid onto the user grid !----------------------------------------------------------------------------- call inter3( nw, wl, xsuser, nwint, wlint, xstmp ) xso2(1:nw-1,iz) = xsuser(1:nw-1)/(wl(2:nw) - wl(1:nw-1)) end do do iz = 1,nz-1 igast = 0 delo2 = .2095_dp * cz(iz) ! vertical o2 column !----------------------------------------------------------------------------- ! ... loop through the internal wavelength grid !----------------------------------------------------------------------------- do iw = 1,nwint-1 !----------------------------------------------------------------------------- ! ... if outside kockarts grid and outside lyman-alpha, use the ! jpl/brasseur+solomon data, if inside ! kockarts grid, use the parameterized values from the call to schu, ! if inside lyman-alpha, use the paraemterized values from call to lymana !----------------------------------------------------------------------------- if( wlint(iw+1) <= wlgast(1) .or. wlint(iw) >= wlgast(ngast) ) then if( wlint(iw+1) <= wlla(1) .or. wlint(iw) >= wlla(nla) ) then dttmp(iw) = xso2int(iw) * delo2 else dttmp(iw) = dto2la(iz,1) end if else igast = igast + 1 dttmp(iw) = dto2k(iz,igast) end if !----------------------------------------------------------------------------- ! ... compute the area in each bin (for correct interpolation purposes only!) !----------------------------------------------------------------------------- dttmp(iw) = dttmp(iw) * (wlint(iw+1) - wlint(iw)) end do !----------------------------------------------------------------------------- ! ... interpolate o2 optical depth from the internal grid onto the user grid !----------------------------------------------------------------------------- call inter3( nw, wl, dtuser, nwint, wlint, dttmp ) dto2(iz,1:nw-1) = dtuser(1:nw-1)/(wl(2:nw) - wl(1:nw-1)) end do end subroutine set_o2_xsect subroutine setcld( nz, z, xlwc, dtcld, omcld, gcld ) use module_wave_data, only : nw !----------------------------------------------------------------------------- != PURPOSE: != Set up cloud optical depth, single albedo and g !----------------------------------------------------------------------------- != PARAMETERS: != PARAMETERS: != NZ - INTEGER, number of specified altitude levels in the working (I) != grid != Z - real(dp), specified altitude working grid (km) (I) != XLWC Cloud water content g/M3 (I) != != dtcld - cloud optical depth != omcld - cloud droplet single albedo != gcld - g !----------------------------------------------------------------------------- ! ! VERTICAL DOMAIN is from bottom(1) to TOP (TOP=nz) ! CCM from top(1) to bottom(nz) !----------------------------------------------------------------------------- !**************************************************************************** implicit none !----------------------------------------------------------------------------- ! ... Dummy arguments !----------------------------------------------------------------------------- integer, intent(in) :: nz real(dp), intent(in) :: z(nz) real(dp), intent(in) :: xlwc(nz) real(dp), intent(out) :: dtcld(nz-1,nw-1) real(dp), intent(out) :: omcld(nz-1,nw-1) real(dp), intent(out) :: gcld(nz-1,nw-1) !----------------------------------------------------------------------------- ! ... Local variables !----------------------------------------------------------------------------- real(dp), parameter :: wden = 1.0_dp * 1.e6_dp ! g/m3 (1 m3 water = 1e6 g water) real(dp), parameter :: re = 10.0_dp * 1.e-6_dp ! assuming cloud drop radius = 10 um to M integer :: i ! vertical index for each blk real(dp) :: dz(nz-1) !----------------------------------------------------------------------------- ! ... calculate vertical interveral dz !----------------------------------------------------------------------------- do i=1,nz-1 dz(i) = (z(i+1)-z(i)) * 1.e3_dp ! dz in Meters end do !---------------------------------------------------- ! ....calculate cloud optical depth T ! following Liao et al. JGR, 104, 23697, 1999 !---------------------------------------------------- do i = 1,nz-1 dtcld(i,:) = 1.5_dp * xlwc(i)*dz(i)/ (wden * re) dtcld(i,:) = max( dtcld(i,:), 0._dp ) omcld(i,:) = .9999_dp gcld (i,:) = .85_dp end do end subroutine setcld subroutine schu_inti !----------------------------------------------------------------------------- ! ... initialize the tables !----------------------------------------------------------------------------- implicit none !----------------------------------------------------------------------------- ! ... local variables !----------------------------------------------------------------------------- integer :: i, iw, k, j1, jp1, j real(dp) :: col real(dp), dimension(6) :: a0, a1, b0, b1 do iw = 1,ngast-1 x_table(0,iw) = sum( a_schu(1:11:2,iw) ) d_table(0,iw) = sum( b_schu(1:11:2,iw) ) do k = 1,tdim col = 22._dp + t_del*real(k-1,kind=dp) o2_table(k) = col col = 10._dp**col a1(:) = a_schu(2:12:2,iw)*col b1(:) = b_schu(2:12:2,iw)*col where( a1(:) < 500._dp ) a0(:) = exp( -a1(:) ) elsewhere a0(:) = 0._dp endwhere where( b1(:) < 500._dp ) b0(:) = exp( -b1(:) ) elsewhere b0(:) = 0._dp endwhere x_table(k,iw) = dot_product( a_schu(1:11:2,iw),a0(:) ) d_table(k,iw) = dot_product( b_schu(1:11:2,iw),b0(:) ) end do end do end subroutine schu_inti subroutine schu( nz, o2col, secchi, dto2, xscho2 ) !----------------------------------------------------------------------------- ! purpose: ! calculate the equivalent absorption cross section of o2 in the sr bands. ! the algorithm is based on: g.kockarts, penetration of solar radiation ! in the schumann-runge bands of molecular oxygen: a robust approximation, ! annales geophysicae, v12, n12, pp. 1207ff, dec 1994. calculation is ! done on the wavelength grid used by kockarts (1994). final values do ! include effects from the herzberg continuum. !----------------------------------------------------------------------------- ! parameters: ! nz - integer, number of specified altitude levels in the working (i) ! grid ! o2col - real(dp), slant overhead o2 column (molec/cc) at each specified (i) ! altitude ! dto2 - real(dp), optical depth due to o2 absorption at each specified (o) ! vertical layer at each specified wavelength ! xscho2 - real(dp), molecular absorption cross section in sr bands at (o) ! each specified wavelength. includes herzberg continuum !----------------------------------------------------------------------------- implicit none !----------------------------------------------------------------------------- ! ... dummy arguments !----------------------------------------------------------------------------- integer, intent(in) :: nz real(dp), intent(inout) :: dto2(nz-1,ngast-1) real(dp), intent(inout) :: xscho2(nz,ngast-1) real(dp), intent(in) :: o2col(nz) real(dp), intent(in) :: secchi(nz) !----------------------------------------------------------------------------- ! ... local variables !----------------------------------------------------------------------------- integer :: i, iw, j, j1, jp1, k, ki, kp1 integer :: index(nz) integer :: minind(1), maxind(1) real(dp) :: a0, a1, b0, b1 real(dp), dimension(6) :: ac, bc real(dp), dimension(nz,6) :: aa, bb real(dp), dimension(nz) :: rjm, rjo2 real(dp), dimension(nz) :: rjmi, rjo2i, lo2col, dels !----------------------------------------------------------------------------- ! ... initialize r(m) !----------------------------------------------------------------------------- rjm(1:nz) = 0._dp rjo2(1:nz) = 0._dp !----------------------------------------------------------------------------- ! ... initialize the table interpolation !----------------------------------------------------------------------------- where( o2col(:) /= 0 ) lo2col(:) = log10( o2col(:) ) endwhere do ki = 1,nz if( o2col(ki) /= 0._dp ) then if( lo2col(ki) <= o2_table(1) ) then dels(ki) = 0._dp index(ki) = 1 else if( lo2col(ki) >= o2_table(tdim) ) then dels(ki) = 1._dp index(ki) = tdim-1 else do k = 2,tdim if( lo2col(ki) <= o2_table(k) ) then dels(ki) = t_fac*(lo2col(ki) - o2_table(k-1)) index(ki) = k-1 exit end if end do end if else index(ki) = 0 dels(ki) = 0._dp end if end do !----------------------------------------------------------------------------- ! ... calculate sum of exponentials (eqs 7 and 8 of kockarts 1994) !----------------------------------------------------------------------------- do iw = 1,ngast-1 do k = 1,nz ki = index(k) rjm(k) = x_table(ki,iw) + dels(k)*(x_table(ki+1,iw) - x_table(ki,iw)) rjo2(k) = d_table(ki,iw) + dels(k)*(d_table(ki+1,iw) - d_table(ki,iw)) end do do k = 1,nz-1 if( rjm(k) > 1.e-100_dp ) then kp1 = k + 1 if( rjm(kp1) > 0._dp ) then dto2(k,iw) = log( rjm(kp1) ) / secchi(kp1) - log( rjm(k) ) * secchi(k) else dto2(k,iw) = 1000._dp end if else dto2(k,iw) = 1000._dp end if end do do k = 1,nz if( rjm(k) > 1.e-100_dp ) then if( rjo2(k) > 1.e-100_dp ) then xscho2(k,iw) = rjo2(k)/rjm(k) else xscho2(k,iw) = 0._dp end if else xscho2(k,iw) = 0._dp end if end do end do end subroutine schu subroutine rtlink( nz, & nw, & zen, & albedo, & z, & nid, & dsdh, & dtrl, & dto3, & dto2, & dtcld, omcld, gcld, & ! dtcb1, omcb1, gcb1, & ! dtcb2, omcb2, gcb2, & ! dtoc1, omoc1, goc1, & ! dtoc2, omoc2, goc2, & ! dtant, omant, gant, & ! dtso4, omso4, gso4, & ! dtsal, omsal, gsal, & ! rajesh: pass aerosol optical properties dtaer, omaer, gaer, & radfld ) implicit none !----------------------------------------------------------------------- ! ... dummy arguments !----------------------------------------------------------------------- integer, intent(in) :: nz integer, intent(in) :: nw real(dp), intent(in) :: z(nz) real(dp), intent(in) :: zen real(dp), intent(in) :: albedo(nw-1) integer, intent(in) :: nid(0:nz-1) real(dp), intent(in) :: dsdh(0:nz-1,nz-1) real(dp), intent(in) :: dtrl(nz-1,nw-1) real(dp), intent(in) :: dto3(nz-1,nw-1) real(dp), intent(in) :: dto2(nz-1,nw-1) real(dp), intent(in) :: dtcld(nz-1,nw-1), omcld(nz-1,nw-1), gcld(nz-1,nw-1) ! real(dp), intent(in) :: dtcb1(nz-1,nw-1), omcb1(nz-1,nw-1), gcb1(nz-1,nw-1) ! real(dp), intent(in) :: dtcb2(nz-1,nw-1), omcb2(nz-1,nw-1), gcb2(nz-1,nw-1) ! real(dp), intent(in) :: dtoc1(nz-1,nw-1), omoc1(nz-1,nw-1), goc1(nz-1,nw-1) ! real(dp), intent(in) :: dtoc2(nz-1,nw-1), omoc2(nz-1,nw-1), goc2(nz-1,nw-1) ! real(dp), intent(in) :: dtant(nz-1,nw-1), omant(nz-1,nw-1), gant(nz-1,nw-1) ! real(dp), intent(in) :: dtso4(nz-1,nw-1), omso4(nz-1,nw-1), gso4(nz-1,nw-1) ! real(dp), intent(in) :: dtsal(nz-1,nw-1), omsal(nz-1,nw-1), gsal(nz-1,nw-1) ! rajesh: declare dust optical properties real(dp), intent(in) :: dtaer(nz-1,nw-1), omaer(nz-1,nw-1), gaer(nz-1,nw-1) real(dp), intent(out) :: radfld(nz,nw-1) !----------------------------------------------------------------------- ! ... local variables !----------------------------------------------------------------------- real(dp), parameter :: smallest = tiny( 1.0 ) integer :: i, ii, iw real(dp) :: dtsct, dtabs real(dp) :: dscld, dacld ! real(dp) :: dscb1, dacb1 ! real(dp) :: dscb2, dacb2 ! real(dp) :: dsoc1, daoc1 ! real(dp) :: dsoc2, daoc2 ! real(dp) :: dsant, daant ! real(dp) :: dsso4, daso4 ! real(dp) :: dssal, dasal ! rajesh: declare scattering and absorbing optical depths real(dp) :: dsaer, daaer real(dp), dimension(nz-1,nw-1) :: dt, om, g !----------------------------------------------------------------------- ! ... set any coefficients specific to rt scheme !----------------------------------------------------------------------- do iw = 1,nw-1 do i = 1,nz-1 dscld = dtcld(i,iw)*omcld(i,iw) dacld = dtcld(i,iw)*(1._dp - omcld(i,iw)) ! dscb1 = dtcb1(i,iw)*omcb1(i,iw) ! dacb1 = dtcb1(i,iw)*(1._dp - omcb1(i,iw)) ! dscb2 = dtcb2(i,iw)*omcb2(i,iw) ! dacb2 = dtcb2(i,iw)*(1._dp - omcb2(i,iw)) ! dsoc1 = dtoc1(i,iw)*omoc1(i,iw) ! daoc1 = dtoc1(i,iw)*(1._dp - omoc1(i,iw)) ! dsoc2 = dtoc2(i,iw)*omoc2(i,iw) ! daoc2 = dtoc2(i,iw)*(1._dp - omoc2(i,iw)) ! dsant = dtant(i,iw)*omant(i,iw) ! daant = dtant(i,iw)*(1._dp - omant(i,iw)) ! dsso4 = dtso4(i,iw)*omso4(i,iw) ! daso4 = dtso4(i,iw)*(1._dp - omso4(i,iw)) ! dssal = dtsal(i,iw)*omsal(i,iw) ! dasal = dtsal(i,iw)*(1._dp - omsal(i,iw)) ! rajesh: determine scattering and absorption optical depths for dust dsaer = dtaer(i,iw)*omaer(i,iw) daaer = dtaer(i,iw)*(1._dp - omaer(i,iw)) dtsct = dtrl(i,iw) + dscld + dsaer ! dscb1 + dscb2 + dsoc1 + dsoc2 + dsant + dsso4 + dssal dtabs = dto3(i,iw) + dto2(i,iw) + dacld + daaer ! dacb1 + dacb2 + daoc1 + daoc2 + daant + daso4 + dasal dtabs = max( dtabs,smallest ) dtsct = max( dtsct,smallest ) !----------------------------------------------------------------------- ! ... invert z-coordinate !----------------------------------------------------------------------- ii = nz - i dt(ii,iw) = dtsct + dtabs if( dtsct /= smallest ) then om(ii,iw) = dtsct/(dtsct + dtabs) g(ii,iw) = ( gcld(i,iw)*dscld + & ! gcb1(i,iw)*dscb1 + & ! gcb2(i,iw)*dscb2 + & ! goc1(i,iw)*dsoc1 + & ! goc2(i,iw)*dsoc2 + & ! gant(i,iw)*dsant + & ! gso4(i,iw)*dsso4 + & ! rajesh: add aerosol contribution gaer(i,iw)*dsaer ) / dtsct g(ii,iw) = max( smallest, g(ii,iw) ) g(ii,iw) = min( 1.d0, g(ii,iw) ) else om(ii,iw) = smallest g(ii,iw) = smallest end if ! print *, dt(ii,iw), om(ii,iw), g(ii,iw) end do end do !----------------------------------------------------------------------- ! ... call rt routine !----------------------------------------------------------------------- call ps2str( nz, nw, zen, albedo, dt, om, & g, dsdh, nid, radfld ) end subroutine rtlink subroutine tridec( syscnt, order, lower, main, upper ) !-------------------------------------------------------------------- ! dimension of lower(syscnt,order), main(syscnt,order), upper(syscnt,order) ! arguments ! ! latest revision june 1995 ! ! purpose trdi computes the solution of the tridiagonal ! linear system, ! b(1)*x(1)+c(1)*x(2) = y(1) ! a(i)*x(i-1)+b(i)*x(i)+c(i)*x(i+1) = y(i) ! i=2,3,...,n-1 ! a(n)*x(n-1)+b(n)*x(n) = y(n) ! ! usage call trdi (n, a, b, c, x ) ! ! arguments ! ! on input n ! the number of unknowns. this subroutine ! requires that n be greater than 2. ! ! a ! the subdiagonal of the matrix is stored in ! locations a(2) through a(n). ! ! b ! the diagonal of the matrix is stored in ! locations b(1) through b(n). ! ! c ! the super-diagonal of the matrix is stored in ! locations c(1) through c(n-1). ! ! x ! the right-hand side of the equations is ! stored in y(1) through y(n). ! ! on output x ! an array which contains the solution to the ! system of equations. ! ! special conditions this subroutine executes satisfactorily ! if the input matrix is diagonally dominant ! and non-singular. the diagonal elements ! are used to pivot, and no tests are made to ! determine singularity. if a singular, or ! nearly singular, matrix is used as input, ! a divide by zero or floating point overflow ! may result. ! ! precision compiler dependent ! ! language fortran ! ! history originally written by nancy werner at ncar ! in october, 1973. modified by s. walters at ! ncar in june 1989 to functionally replace ! the cal routine tridla. ! ! portability fortran 90 ! ! algorithm an lu-decomposition is obtained using the ! algorithm described in the reference below. ! ! to avoid unnecessary divisions, the alpha ! values used in the routine are the ! reciprocals of the alpha values described ! in the reference below. ! ! accuracy every component of the residual of the linear ! system (i.e. the difference between y and ! the matrix applied to x) should be less in ! magnitude than ten times the machine precision ! times the matrix order times the maximum ! absolute component of the solution vector ! times the largest absolute row sum of the ! input matrix. ! ! timing the timing is roughly proportional to the ! order n of the linear system. ! ! references analysis of numerical methods by ! e. isaacson and h. keller ! (john wiley and sons, 1966) pp. 55-58. !-------------------------------------------------------------------- implicit none !-------------------------------------------------------------------- ! ... dummy args !-------------------------------------------------------------------- integer, intent(in) :: syscnt, order real(dp), intent(in), dimension(syscnt,order) :: lower real(dp), intent(inout), dimension(syscnt,order) :: main, upper !-------------------------------------------------------------------- ! ... local variables !-------------------------------------------------------------------- integer :: i !---------------------------------------------------------------------- ! ... lu-decomposition !---------------------------------------------------------------------- main(:,1) = 1._dp / main(:,1) upper(:,1) = upper(:,1)*main(:,1) do i = 2,order-1 main(:,i) = 1._dp / (main(:,i) - lower(:,i)*upper(:,i-1)) upper(:,i) = upper(:,i)*main(:,i) end do end subroutine tridec subroutine trislv( syscnt, order, lower, main, upper, x ) implicit none !---------------------------------------------------------------------- ! ... dummy args !---------------------------------------------------------------------- integer, intent(in) :: syscnt, order real(dp), intent(in), dimension(syscnt,order) :: lower, & main, & upper real(dp), intent(inout), dimension(syscnt,order) :: x !---------------------------------------------------------------------- ! ... local variables !---------------------------------------------------------------------- integer :: i, im1, j, n, nm1 nm1 = order - 1 n = order !---------------------------------------------------------------------- ! ... solve the system !---------------------------------------------------------------------- x(:,1) = x(:,1)*main(:,1) do i = 2,nm1 x(:,i) = (x(:,i) - lower(:,i)*x(:,i-1))*main(:,i) end do x(:,n) = (x(:,n) - lower(:,n)*x(:,nm1))/(main(:,n) - lower(:,n)*upper(:,nm1)) do i = nm1,1,-1 x(:,i) = x(:,i) - upper(:,i)*x(:,i+1) end do end subroutine trislv subroutine ps2str( nz, nw, zen, rsfc, tauu, omu, & gu, dsdh, nid, radfld ) !----------------------------------------------------------------------------- ! purpose: ! solve two-stream equations for multiple layers. the subroutine is based ! on equations from: toon et al., j.geophys.res., v94 (d13), nov 20, 1989. ! it contains 9 two-stream methods to choose from. a pseudo-spherical ! correction has also been added. !----------------------------------------------------------------------------- ! parameters: ! nz - integer, number of specified altitude levels in the working (i) ! grid ! zen - real(dp), solar zenith angle (degrees) (i) ! rsfc - real(dp), surface albedo at current wavelength (i) ! tauu - real(dp), unscaled optical depth of each layer (i) ! omu - real(dp), unscaled single scattering albedo of each layer (i) ! gu - real(dp), unscaled asymmetry parameter of each layer (i) ! dsdh - real(dp), slant path of direct beam through each layer crossed (i) ! when travelling from the top of the atmosphere to layer i; ! dsdh(i,j), i = 0..nz-1, j = 1..nz-1 ! nid - integer, number of layers crossed by the direct beam when (i) ! travelling from the top of the atmosphere to layer i; ! nid(i), i = 0..nz-1 ! delta - logical, switch to use delta-scaling (i) ! .true. -> apply delta-scaling ! .false.-> do not apply delta-scaling ! fdr - real(dp), contribution of the direct component to the total (o) ! actinic flux at each altitude level ! fup - real(dp), contribution of the diffuse upwelling component to (o) ! the total actinic flux at each altitude level ! fdn - real(dp), contribution of the diffuse downwelling component to (o) ! the total actinic flux at each altitude level ! edr - real(dp), contribution of the direct component to the total (o) ! spectral irradiance at each altitude level ! eup - real(dp), contribution of the diffuse upwelling component to (o) ! the total spectral irradiance at each altitude level ! edn - real(dp), contribution of the diffuse downwelling component to (o) ! the total spectral irradiance at each altitude level !----------------------------------------------------------------------------- implicit none !----------------------------------------------------------------------------- ! ... dummy arguments !----------------------------------------------------------------------------- integer, intent(in) :: nz, nw integer, intent(in) :: nid(0:nz-1) real(dp), intent(in) :: zen real(dp), intent(in) :: rsfc(nw-1) real(dp), dimension(nz-1,nw-1), intent(in) :: tauu, omu, gu real(dp), intent(in) :: dsdh(0:nz-1,nz-1) real(dp), intent(out) :: radfld(nz,nw-1) !----------------------------------------------------------------------------- ! ... local variables !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- ! ... mu = cosine of solar zenith angle ! rsfc = surface albedo ! tauu = unscaled optical depth of each layer ! omu = unscaled single scattering albedo ! gu = unscaled asymmetry factor ! klev = max dimension of number of layers in atmosphere ! nlayer = number of layers in the atmosphere ! nlevel = nlayer + 1 = number of levels !----------------------------------------------------------------------------- real(dp), parameter :: smallest = tiny( 1.0_dp ) real(dp), parameter :: largest = huge( 1.0_dp ) real(dp), parameter :: d2r = 3.1415926_dp/180.0_dp integer :: mrows, nzm1, nzm2 real(dp), parameter :: eps = 1.e-3_dp real(dp), parameter :: pifs = 1._dp real(dp), parameter :: fdn0 = 0._dp integer :: row integer :: lev integer :: i, ip1, iw integer :: j, jl, ju real(dp) :: precis, wrk real(dp) :: tempg real(dp) :: mu, suma real(dp) :: g, om real(dp) :: gam1, gam2, gam3, gam4 real(dp), dimension(nz-1) :: f, gi, omi real(dp), dimension(0:nz) :: tauc, mu2 real(dp), dimension(nz-1) :: lam, taun, bgam real(dp), dimension(nz-1) :: cdn real(dp), dimension(0:nz,nw-1) :: tausla real(dp), dimension(nz-1,nw-1) :: cup, cuptn, cdntn real(dp), dimension(nz-1,nw-1) :: e1, e2, e3, e4 real(dp), dimension(2*(nz-1)) :: a, b, d, e real(dp), dimension(nw-1,2*(nz-1)) :: sub, main, super, y !----------------------------------------------------------------------------- ! ... for calculations of associated legendre polynomials for gama1,2,3,4 ! in delta-function, modified quadrature, hemispheric constant, ! hybrid modified eddington-delta function metods, p633,table1. ! w.e.meador and w.r.weaver, gas,1980,v37,p.630 ! w.j.wiscombe and g.w. grams, gas,1976,v33,p2440, ! uncomment the following two lines and the appropriate statements ! further down. !----------------------------------------------------------------------------- real(dp) :: expon, expon0, expon1, divisr, temp, up, dn real(dp) :: ssfc !----------------------------------------------------------------------------- ! ... initial conditions: pi*solar flux = 1; diffuse incidence = 0 !----------------------------------------------------------------------------- nzm1 = nz - 1 nzm2 = nz - 2 mrows = 2*nzm1 precis = epsilon( precis ) mu = cos( zen*d2r ) wave_loop : & do iw = 1,nw-1 !----------------------------------------------------------------------------- ! ... compute coefficients for each layer: ! gam1 - gam4 = 2-stream coefficients, different for different approximations ! expon0 = calculation of e when tau is zero ! expon1 = calculation of e when tau is taun ! cup and cdn = calculation when tau is zero ! cuptn and cdntn = calc. when tau is taun ! divisr = prevents division by zero !----------------------------------------------------------------------------- tauc(0:nz) = 0._dp tausla(0:nz,iw) = 0._dp mu2(0:nz) = sqrt( smallest ) !----------------------------------------------------------------------------- ! ... delta-scaling. have to be done for delta-eddington approximation, ! delta discrete ordinate, practical improved flux method, delta function, ! and hybrid modified eddington-delta function methods approximations !----------------------------------------------------------------------------- f(1:nzm1) = gu(:,iw)*gu(:,iw) gi(1:nzm1) = (gu(:,iw) - f(1:nzm1))/(1._dp - f(1:nzm1)) omi(1:nzm1) = (1._dp - f(1:nzm1))*omu(1:nzm1,iw)/(1._dp - omu(1:nzm1,iw)*f(1:nzm1)) taun(1:nzm1) = (1._dp - omu(1:nzm1,iw)*f(1:nzm1))*tauu(1:nzm1,iw) !----------------------------------------------------------------------------- ! ... calculate slant optical depth at the top of the atmosphere when zen>90. ! in this case, higher altitude of the top layer is recommended which can ! be easily changed in gridz.f. !----------------------------------------------------------------------------- if( zen > 90._dp ) then if( nid(0) < 0 ) then tausla(0,iw) = largest else ju = nid(0) tausla(0,iw) = 2._dp*dot_product( taun(1:ju),dsdh(0,1:ju) ) end if end if level_loop : & do i = 1,nzm1 g = gi(i) om = omi(i) tauc(i) = tauc(i-1) + taun(i) !----------------------------------------------------------------------------- ! ... stay away from 1 by precision. for g, also stay away from -1 !----------------------------------------------------------------------------- tempg = min( abs(g),1._dp - precis ) g = sign( tempg,g ) om = min( om,1._dp - precis ) !----------------------------------------------------------------------------- ! ... calculate slant optical depth !----------------------------------------------------------------------------- if( nid(i) < 0 ) then tausla(i,iw) = largest else ju = min( nid(i),i ) suma = dot_product( taun(1:ju),dsdh(i,1:ju) ) jl = min( nid(i),i ) + 1 tausla(i,iw) = suma + 2._dp*dot_product( taun(jl:nid(i)),dsdh(i,jl:nid(i)) ) if( abs(tausla(i,iw)-tausla(i-1,iw)) < 1.0e-30_dp ) then mu2(i) = sqrt( largest ) else mu2(i) = (tauc(i) - tauc(i-1))/(tausla(i,iw) - tausla(i-1,iw)) mu2(i) = sign( max( abs(mu2(i)),sqrt(smallest) ),mu2(i) ) end if end if !----------------------------------------------------------------------------- ! ... the following gamma equations are from pg 16,289, table 1 ! eddington approximation(joseph et al., 1976, jas, 33, 2452): !----------------------------------------------------------------------------- gam1 = .25_dp*(7._dp - om*(4._dp + 3._dp*g)) gam2 = -.25_dp*(1._dp - om*(4._dp - 3._dp*g)) gam3 = .25_dp*(2._dp - 3._dp*g*mu) gam4 = 1._dp - gam3 !----------------------------------------------------------------------------- ! ... lambda = pg 16,290 equation 21 ! big gamma = pg 16,290 equation 22 !----------------------------------------------------------------------------- lam(i) = sqrt( gam1*gam1 - gam2*gam2 ) bgam(i) = (gam1 - lam(i))/gam2 wrk = lam(i)*taun(i) if( wrk < 500._dp ) then expon = exp( -wrk ) else expon = 0._dp end if !----------------------------------------------------------------------------- ! ... e1 - e4 = pg 16,292 equation 44 !----------------------------------------------------------------------------- e1(i,iw) = 1._dp + bgam(i)*expon e2(i,iw) = 1._dp - bgam(i)*expon e3(i,iw) = bgam(i) + expon e4(i,iw) = bgam(i) - expon !----------------------------------------------------------------------------- ! ... the following sets up for the c equations 23, and 24 ! found on page 16,290 ! prevent division by zero (if lambda=1/mu, shift 1/mu^2 by eps = 1.e-3 ! which is approx equiv to shifting mu by 0.5*eps* (mu)**3 !----------------------------------------------------------------------------- if( tausla(i-1,iw) < 500._dp ) then expon0 = exp( -tausla(i-1,iw) ) else expon0 = 0._dp end if if( tausla(i,iw) < 500._dp ) then expon1 = exp( -tausla(i,iw) ) else expon1 = 0._dp end if divisr = lam(i)*lam(i) - 1._dp/(mu2(i)*mu2(i)) temp = max( eps,abs(divisr) ) divisr = 1._dp/sign( temp,divisr ) up = om*pifs*((gam1 - 1._dp/mu2(i))*gam3 + gam4*gam2)*divisr dn = om*pifs*((gam1 + 1._dp/mu2(i))*gam4 + gam2*gam3)*divisr !----------------------------------------------------------------------------- ! ... cup and cdn are when tau is equal to zero ! cuptn and cdntn are when tau is equal to taun !----------------------------------------------------------------------------- cup(i,iw) = up*expon0 cdn(i) = dn*expon0 cuptn(i,iw) = up*expon1 cdntn(i,iw) = dn*expon1 end do level_loop !----------------------------------------------------------------------------- ! ... set up matrix ! ssfc = pg 16,292 equation 37 where pi fs is one (unity). !----------------------------------------------------------------------------- if( tausla(nzm1,iw) < 500._dp ) then ssfc = rsfc(iw)*mu*exp( -tausla(nzm1,iw) )*pifs else ssfc = 0._dp end if !----------------------------------------------------------------------------- ! ... the following are from pg 16,292 equations 39 - 43. ! set up first row of matrix: !----------------------------------------------------------------------------- a(1) = 0._dp b(1) = e1(1,iw) d(1) = -e2(1,iw) e(1) = fdn0 - cdn(1) !----------------------------------------------------------------------------- ! ... set up odd rows 3 thru (mrows - 1): !----------------------------------------------------------------------------- a(3:mrows-1:2) = e2(1:nzm2,iw)*e3(1:nzm2,iw) - e4(1:nzm2,iw)*e1(1:nzm2,iw) b(3:mrows-1:2) = e1(1:nzm2,iw)*e1(2:nzm1,iw) - e3(1:nzm2,iw)*e3(2:nzm1,iw) d(3:mrows-1:2) = e3(1:nzm2,iw)*e4(2:nzm1,iw) - e1(1:nzm2,iw)*e2(2:nzm1,iw) e(3:mrows-1:2) = e3(1:nzm2,iw)*(cup(2:nzm1,iw) - cuptn(1:nzm2,iw)) + e1(1:nzm2,iw)*(cdntn(1:nzm2,iw) - cdn(2:nzm1)) !----------------------------------------------------------------------------- ! ... set up even rows 2 thru (mrows - 2): !----------------------------------------------------------------------------- a(2:mrows-2:2) = e2(2:nzm1,iw)*e1(1:nzm2,iw) - e3(1:nzm2,iw)*e4(2:nzm1,iw) b(2:mrows-2:2) = e2(1:nzm2,iw)*e2(2:nzm1,iw) - e4(1:nzm2,iw)*e4(2:nzm1,iw) d(2:mrows-2:2) = e1(2:nzm1,iw)*e4(2:nzm1,iw) - e2(2:nzm1,iw)*e3(2:nzm1,iw) e(2:mrows-2:2) = (cup(2:nzm1,iw) - cuptn(1:nzm2,iw))*e2(2:nzm1,iw) - (cdn(2:nzm1) - cdntn(1:nzm2,iw))*e4(2:nzm1,iw) !----------------------------------------------------------------------------- ! ... set up last row of matrix at mrows: !----------------------------------------------------------------------------- a(mrows) = e1(nzm1,iw) - rsfc(iw)*e3(nzm1,iw) b(mrows) = e2(nzm1,iw) - rsfc(iw)*e4(nzm1,iw) d(mrows) = 0._dp e(mrows) = ssfc - cuptn(nzm1,iw) + rsfc(iw)*cdntn(nzm1,iw) sub(iw,1:mrows) = a(1:mrows) main(iw,1:mrows) = b(1:mrows) super(iw,1:mrows) = d(1:mrows) y(iw,1:mrows) = e(1:mrows) end do wave_loop !----------------------------------------------------------------------------- ! ... solve the system !----------------------------------------------------------------------------- call tridec( nw-1, mrows, sub, main, super ) call trislv( nw-1, mrows, sub, main, super, y ) !----------------------------------------------------------------------------- ! ... unfold solution of matrix, compute output fluxes !----------------------------------------------------------------------------- do iw = 1,nw-1 !----------------------------------------------------------------------------- ! ... the following equations are from pg 16,291 equations 31 & 32 !----------------------------------------------------------------------------- e(:mrows) = y(iw,:mrows) if( tausla(0,iw) < 500._dp ) then radfld(1,iw) = 2._dp*(fdn0 + e(1)*e3(1,iw) - e(2)*e4(1,iw) + cup(1,iw)) + exp( -tausla(0,iw) ) else radfld(1,iw) = 2._dp*(fdn0 + e(1)*e3(1,iw) - e(2)*e4(1,iw) + cup(1,iw)) end if where( tausla(1:nzm1,iw) < 500._dp ) cdn(1:nzm1) = exp( -tausla(1:nzm1,iw) ) elsewhere cdn(1:nzm1) = 0._dp endwhere radfld(2:nz,iw) = 2._dp*(e(1:mrows-1:2)*(e3(1:nzm1,iw) + e1(1:nzm1,iw)) & + e(2:mrows:2)*(e4(1:nzm1,iw) + e2(1:nzm1,iw)) & + cdntn(1:nzm1,iw) + cuptn(1:nzm1,iw)) + cdn(1:nzm1) end do end subroutine ps2str subroutine pchem( chem_opt, nz, tlev, airlev ) use module_state_description, only : MOZART_KPP, MOZCART_KPP, & MOZART_MOSAIC_4BIN_KPP, MOZART_MOSAIC_4BIN_AQ_KPP use module_wave_data, only : r01, r04, r44, r08, r06, r10, & r11, r14, r15, r17, r18, & xs_mvk, xs_macr, xs_hyac, xs_glyald, nj !----------------------------------------------------------------------------- ! purpose: ! load various "weighting functions" (products of cross section and ! quantum yield at each altitude and for wavelength). the altitude ! dependence is necessary to ensure the consideration of pressure and ! temperature dependence of the cross sections or quantum yields. ! the actual reading, evaluation and interpolation is done is separate ! subroutines for ease of management and manipulation. please refer to ! the inline documentation of the specific subroutines for detail information. !----------------------------------------------------------------------------- ! parameters: ! nw - integer, number of specified intervals + 1 in working (i) ! wavelength grid ! wl - real(dp), vector of lower limits of wavelength intervals in (i) ! working wavelength grid ! nz - integer, number of altitude levels in working altitude grid (i) ! tlev - real(dp), temperature (k) at each specified altitude level (i) ! airlev - real(dp), air density (molec/cc) at each altitude level (i) ! j - integer, counter for number of weighting functions defined (io) ! sq - real(dp), cross section x quantum yield (cm^2) for each (o) ! photolysis reaction defined, at each defined wavelength and ! at each defined altitude level ! jlabel - character*40, string identifier for each photolysis reaction (o) ! defined !----------------------------------------------------------------------------- implicit none !----------------------------------------------------------------------------- ! ... dummy arguments !----------------------------------------------------------------------------- integer, intent(in) :: chem_opt integer, intent(in) :: nz real(dp), intent(in) :: tlev(nz) real(dp), intent(in) :: airlev(nz) !----------------------------------------------------------------------------- ! ... local variables !----------------------------------------------------------------------------- integer :: j character(len=40) :: jlabel(nj) !----------------------------------------------------------------------------- ! ... (1) o2 + hv -> o + o ! reserve first position. cross section parameterization in ! schumman-runge and lyman-alpha regions are zenith-angle dependent !----------------------------------------------------------------------------- j = 1 jlabel(j) = 'o2 + hv -> o + o ' !----------------------------------------------------------------------------- ! ... (2,3) o3 + hv -> (both channels) (tmp dep) !----------------------------------------------------------------------------- call r01( nz, tlev, airlev, j, jlabel ) !----------------------------------------------------------------------------- ! ... (4) no2 + hv -> no + o(3p) !----------------------------------------------------------------------------- j = j + 1 jlabel(j) = 'no2 + hv -> no + o(3p) ' !----------------------------------------------------------------------------- ! ... (5) no3 + hv -> no2 + o(3p) !----------------------------------------------------------------------------- j = j + 1 jlabel(j) = 'no3 + hv -> no2 + o(3p) ' !----------------------------------------------------------------------------- ! ... (6) no3 + hv -> no + o2 !----------------------------------------------------------------------------- j = j + 1 jlabel(j) = 'no3 + hv -> no + o2 ' !----------------------------------------------------------------------------- ! ... (7,8) n2o5 + hv -> (both channels) (tmp dep) !----------------------------------------------------------------------------- call r04( nz,tlev,airlev,j,jlabel ) !----------------------------------------------------------------------------- ! ... (9) n2o + hv -> n2 + o(1d) (tmp dep) !----------------------------------------------------------------------------- call r44( nz,tlev,airlev,j,jlabel ) !----------------------------------------------------------------------------- ! ... (10) ho2 + hv -> oh + o !----------------------------------------------------------------------------- j = j + 1 jlabel(j) = 'ho2 + hv -> oh + o ' !----------------------------------------------------------------------------- ! ... (11) h2o2 + hv -> 2 oh !----------------------------------------------------------------------------- call r08( nz,tlev,airlev,j,jlabel ) !----------------------------------------------------------------------------- ! ... (12) hno2 + hv -> oh + no !----------------------------------------------------------------------------- j = j + 1 jlabel(j) = 'hno2 + hv -> oh + no ' !----------------------------------------------------------------------------- ! ... (13) hno3 + hv -> oh + no2 !----------------------------------------------------------------------------- call r06( nz,tlev,airlev,j,jlabel ) !----------------------------------------------------------------------------- ! ... (14) hno4 + hv -> ho2 + no2 !----------------------------------------------------------------------------- j = j + 1 jlabel(j) = 'hno4 + hv -> ho2 + no2 ' !----------------------------------------------------------------------------- ! ... (15,16) ch2o + hv -> (both channels) !----------------------------------------------------------------------------- call r10( nz,tlev,airlev,j,jlabel ) !----------------------------------------------------------------------------- ! ... (17,18,19) ch3cho + hv -> (all three channels) !----------------------------------------------------------------------------- call r11( nz,tlev,airlev,j,jlabel ) !----------------------------------------------------------------------------- ! ... (20) c2h5cho + hv -> c2h5 + hco !----------------------------------------------------------------------------- j = j + 1 jlabel(j) = 'c2h5cho + hv -> c2h5 + hco ' !----------------------------------------------------------------------------- ! ... (21) chocho + hv -> products !----------------------------------------------------------------------------- j = j + 1 jlabel(j) = 'chocho + hv -> products ' !----------------------------------------------------------------------------- ! ... (22) ch3cocho + hv -> products !----------------------------------------------------------------------------- call r14( nz,tlev,airlev,j,jlabel ) !----------------------------------------------------------------------------- ! ... (23) ch3coch3 + hv -> products !----------------------------------------------------------------------------- call r15( nz,tlev,airlev,j,jlabel ) !----------------------------------------------------------------------------- ! ... (24) ch3ooh + hv -> ch3o + oh !----------------------------------------------------------------------------- j = j + 1 jlabel(j) = 'ch3ooh + hv -> ch3o + oh ' !----------------------------------------------------------------------------- ! ... (25) ch3ono2 + hv -> ch3o + no2 !----------------------------------------------------------------------------- call r17( nz,tlev,airlev,j,jlabel ) !----------------------------------------------------------------------------- ! ... (26) pan + hv -> products !----------------------------------------------------------------------------- call r18( nz,tlev,airlev,j,jlabel ) !----------------------------------------------------------------------------- ! ... (27) mvk + hv -> products ! (28) macr + hv -> products ! (29) hyac + hv -> products ! (30) glyald + hv -> products !----------------------------------------------------------------------------- is_mozart : & if( chem_opt == MOZART_KPP .or. chem_opt == MOZCART_KPP .or. & chem_opt == MOZART_MOSAIC_4BIN_KPP .or. & chem_opt == MOZART_MOSAIC_4BIN_AQ_KPP) then call xs_mvk( nz, tlev, airlev ) call xs_macr( nz, tlev, airlev ) call xs_hyac( nz, tlev, airlev ) call xs_glyald( nz, tlev, airlev ) else is_mozart !----------------------------------------------------------------------------- ! ... (27) cloo + hv -> products !----------------------------------------------------------------------------- j = j + 1 jlabel(j) = 'cloo + hv -> products ' !----------------------------------------------------------------------------- ! ... (28,29) clono2 + hv -> products !----------------------------------------------------------------------------- j = j + 1 jlabel(j) = 'clono2 + hv -> cl + no3' j = j + 1 jlabel(j) = 'clono2 + hv -> clo + no2' !----------------------------------------------------------------------------- ! ... (30) ch3cl + hv -> products !----------------------------------------------------------------------------- j = j + 1 jlabel(j) = 'ch3cl + hv -> products ' !----------------------------------------------------------------------------- ! ... (31) ccl2o + hv -> products !----------------------------------------------------------------------------- j = j + 1 jlabel(j) = 'ccl2o + hv -> products' !----------------------------------------------------------------------------- ! ... (32) ccl4 + hv -> products !----------------------------------------------------------------------------- j = j + 1 jlabel(j) = 'ccl4 + hv -> products ' !----------------------------------------------------------------------------- ! ... (33) cclfo + hv -> products !----------------------------------------------------------------------------- j = j + 1 jlabel(j) = 'cclfo + hv -> products ' !----------------------------------------------------------------------------- ! ... (34) ccf2o + hv -> products !----------------------------------------------------------------------------- j = j + 1 jlabel(j) = 'ccf2o + hv -> products ' !----------------------------------------------------------------------------- ! ... (35) cf2clcfcl2 (cfc-113) + hv -> products !----------------------------------------------------------------------------- j = j + 1 jlabel(j) = 'cf2clcfcl2 (cfc-113) + hv -> products' !----------------------------------------------------------------------------- ! ... (36) cf2clcf2cl (cfc-114) + hv -> products !----------------------------------------------------------------------------- j = j + 1 jlabel(j) = 'cf2clcf2cl (cfc-114) + hv -> products ' !----------------------------------------------------------------------------- ! ... (37) cf3cf2cl (cfc-115) + hv -> products !----------------------------------------------------------------------------- j = j + 1 jlabel(j) = 'cf3cf2cl (cfc-115) + hv -> products ' !----------------------------------------------------------------------------- ! ... (38) ccl3f (cfc-111) + hv -> products !----------------------------------------------------------------------------- j = j + 1 jlabel(j) = 'ccl3f (cfc-111) + hv -> products' !----------------------------------------------------------------------------- ! ... (39) ccl2f2 (cfc-112) + hv -> products !----------------------------------------------------------------------------- j = j + 1 jlabel(j) = 'ccl2f2 (cfc-112) + hv -> products' !----------------------------------------------------------------------------- ! (40) ch3ccl3 + hv -> products !----------------------------------------------------------------------------- j = j + 1 jlabel(j)='ch3ccl3 + hv -> products' !----------------------------------------------------------------------------- ! (41) cf3chcl2 (hcfc-123) + hv -> products !----------------------------------------------------------------------------- j = j + 1 jlabel(j)='cf3chcl2 (hcfc-123) + hv -> products' !----------------------------------------------------------------------------- ! (42) cf3chfcl (hcfc-124) + hv -> products !----------------------------------------------------------------------------- j = j + 1 jlabel(j)='cf3chfcl (hcfc-124) + hv -> products' !----------------------------------------------------------------------------- ! (43) ch3cfcl2 (hcfc-141b) + hv -> products !----------------------------------------------------------------------------- j = j + 1 jlabel(j)='ch3cfcl2 (hcfc-141b) + hv -> products' !----------------------------------------------------------------------------- ! (44) ch3cf2cl (hcfc-142b) + hv -> products !----------------------------------------------------------------------------- j = j + 1 jlabel(j)='ch3cf2cl (hcfc-142b) + hv -> products' !----------------------------------------------------------------------------- ! (45) cf3cf2chcl2 (hcfc-225ca) + hv -> products !----------------------------------------------------------------------------- j = j + 1 jlabel(j)='cf3cf2chcl2 (hcfc-225ca) + hv -> products' !----------------------------------------------------------------------------- ! (46) cf2clcf2chfcl (hcfc-225cb) + hv -> products !----------------------------------------------------------------------------- j = j + 1 jlabel(j)='cf2clcf2chfcl (hcfc-225cb) + hv -> products' !----------------------------------------------------------------------------- ! (47) chclf2 (hcfc-22) + hv -> products !----------------------------------------------------------------------------- j = j + 1 jlabel(j)='chclf2 (hcfc-22) + hv -> products' !----------------------------------------------------------------------------- ! (48) brono2 + hv -> br + o + no2 !----------------------------------------------------------------------------- j = j + 1 jlabel(j)='brono2 + hv -> br + o + no2' !----------------------------------------------------------------------------- ! (49) brono2 + hv -> bro + no2 !----------------------------------------------------------------------------- j = j + 1 jlabel(j)='brono2 + hv -> bro + no2' !----------------------------------------------------------------------------- ! (50) ch3br + hv -> products !----------------------------------------------------------------------------- j = j + 1 jlabel(j)=' ch3br + hv -> products' !----------------------------------------------------------------------------- ! (51) chbr3 + hv -> products !----------------------------------------------------------------------------- j = j + 1 jlabel(j)='chbr3 + hv -> products' !----------------------------------------------------------------------------- ! (52) cf3br (halon-1301) + hv -> products !----------------------------------------------------------------------------- j = j + 1 jlabel(j)='cf3br (halon-1301) + hv -> products' !----------------------------------------------------------------------------- ! (53) cf2brcf2br (halon-2402) + hv -> products !----------------------------------------------------------------------------- j = j + 1 jlabel(j)='cf2brcf2br (halon-2402) + hv -> products' !----------------------------------------------------------------------------- ! (54) cf2br2 (halon-1202) + hv -> products !----------------------------------------------------------------------------- j = j + 1 jlabel(j)='cf2br2 (halon-1202) + hv -> products' !----------------------------------------------------------------------------- ! (55) cf2brcl (halon-1211) + hv -> products !----------------------------------------------------------------------------- j = j + 1 jlabel(j)='cf2brcl (halon-1211) + hv -> products ' !----------------------------------------------------------------------------- ! (56) cl2 + hv -> cl + cl !----------------------------------------------------------------------------- j = j + 1 jlabel(j)='cl2 + hv -> cl + cl ' !----------------------------------------------------------------------------- ! (57) hocl + hv -> oh + cl !----------------------------------------------------------------------------- j = j + 1 jlabel(j)='hocl + hv -> oh + cl' !----------------------------------------------------------------------------- ! (58) fmcl + hv -> cl + co + ho2 !----------------------------------------------------------------------------- j = j + 1 jlabel(j)='fmcl + hv -> cl + co + ho2' end if is_mozart end subroutine pchem subroutine lymana( nz, o2col, secchi, dto2la, xso2la ) !----------------------------------------------------------------------------- ! purpose: ! calculate the effective absorption cross section of o2 in the lyman-alpha ! bands and an effective o2 optical depth at all altitudes. parameterized ! after: chabrillat, s., and g. kockarts, simple parameterization of the ! absorption of the solar lyman-alpha line, geophysical research letters, ! vol.24, no.21, pp 2659-2662, 1997. !----------------------------------------------------------------------------- ! parameters: ! nz - integer, number of specified altitude levels in the working (i) ! grid ! o2col - real(dp), slant overhead o2 column (molec/cc) at each specified (i) ! altitude ! dto2la - real(dp), optical depth due to o2 absorption at each specified (o) ! vertical layer ! xso2la - real(dp), molecular absorption cross section in la bands (o) !----------------------------------------------------------------------------- implicit none !----------------------------------------------------------------------------- ! ... dummy arguments !----------------------------------------------------------------------------- integer, intent(in) :: nz real(dp), intent(in) :: o2col(nz) real(dp), intent(in) :: secchi(nz) real(dp), intent(out) :: dto2la(nz-1,nla-1) real(dp), intent(out) :: xso2la(nz,nla-1) !----------------------------------------------------------------------------- ! ... local variables !----------------------------------------------------------------------------- integer :: i, k, kp1 real(dp), dimension(nz) :: rm, ro2 real(dp), save :: b(3), c(3), d(3), e(3) data b / 6.8431e-01_dp, 2.29841e-01_dp, 8.65412e-02_dp/, & c /8.22114e-21_dp, 1.77556e-20_dp, 8.22112e-21_dp/, & d / 6.0073e-21_dp, 4.28569e-21_dp, 1.28059e-20_dp/, & e /8.21666e-21_dp, 1.63296e-20_dp, 4.85121e-17_dp/ !----------------------------------------------------------------------------- ! ... calculate reduction factors at every altitude !----------------------------------------------------------------------------- rm(:) = 0._dp ro2(:) = 0._dp do k = 1,nz do i = 1,3 rm(k) = rm(k) + b(i) * exp( -c(i)*o2col(k) ) ro2(k) = ro2(k) + d(i) * exp( -e(i)*o2col(k) ) end do end do !----------------------------------------------------------------------------- ! ... calculate effective o2 optical depths and effective o2 cross sections !----------------------------------------------------------------------------- do k = 1,nz-1 if( rm(k) > 1.e-100_dp ) then kp1 = k + 1 if( rm(kp1) > 0._dp ) then dto2la(k,1) = log( rm(kp1) )/secchi(kp1) - log( rm(k) )/secchi(k) else dto2la(k,1) = 1000._dp end if else dto2la(k,1) = 1000._dp end if end do do k = 1,nz if( rm(k) > 1.e-100_dp ) then if( ro2(k) > 1.e-100_dp ) then xso2la(k,1) = ro2(k)/rm(k) else xso2la(k,1) = 0._dp end if else xso2la(k,1) = 0._dp end if end do end subroutine lymana subroutine inter2( ng, xg, yg, n, x, y, ierr ) !----------------------------------------------------------------------------- ! purpose: ! map input data given on single, discrete points onto a set of target ! bins. ! the original input data are given on single, discrete points of an ! arbitrary grid and are being linearly interpolated onto a specified set ! of target bins. in general, this is the case for most of the weighting ! functions (action spectra, molecular cross section, and quantum yield ! data), which have to be matched onto the specified wavelength intervals. ! the average value in each target bin is found by averaging the trapezoi- ! dal area underneath the input data curve (constructed by linearly connec- ! ting the discrete input values). ! some caution should be used near the endpoints of the grids. if the ! input data set does not span the range of the target grid, an error ! message is printed and the execution is stopped, as extrapolation of the ! data is not permitted. ! if the input data does not encompass the target grid, use addpnt to ! expand the input array. !----------------------------------------------------------------------------- ! parameters: ! ng - integer, number of bins + 1 in the target grid (i) ! xg - real(dp), target grid (e.g., wavelength grid); bin i is defined (i) ! as [xg(i),xg(i+1)] (i = 1..ng-1) ! yg - real(dp), y-data re-gridded onto xg, yg(i) specifies the value for (o) ! bin i (i = 1..ng-1) ! n - integer, number of points in input grid (i) ! x - real(dp), grid on which input data are defined (i) ! y - real(dp), input y-data (i) !----------------------------------------------------------------------------- implicit none !----------------------------------------------------------------------------- ! ... dummy arguments !----------------------------------------------------------------------------- integer, intent(in) :: ng, n integer, intent(out) :: ierr real(dp), intent(in) :: x(n) real(dp), intent(in) :: y(n) real(dp), intent(in) :: xg(ng) real(dp), intent(out) :: yg(ng) !----------------------------------------------------------------------------- ! ... local variables !----------------------------------------------------------------------------- integer :: ngintv integer :: i, k, jstart real(dp) :: area, xgl, xgu real(dp) :: darea, slope real(dp) :: a1, a2, b1, b2 ierr = 0 !----------------------------------------------------------------------------- ! ... test for correct ordering of data, by increasing value of x !----------------------------------------------------------------------------- do i = 2,n if( x(i) <= x(i-1) ) then ierr = 1 write(*,*) 'inter2: x coord not monotonically increasing' return end if end do do i = 2,ng if( xg(i) <= xg(i-1) ) then ierr = 2 write(*,*) 'inter2: xg coord not monotonically increasing' return end if end do !----------------------------------------------------------------------------- ! ... check for xg-values outside the x-range !----------------------------------------------------------------------------- if( x(1) > xg(1) .or. x(n) < xg(ng) ) then write(*,*) 'inter2: data does not span grid' write(*,*) ' use addpnt to expand data and re-run' ! call endrun end if !----------------------------------------------------------------------------- ! ... find the integral of each grid interval and use this to ! calculate the average y value for the interval ! xgl and xgu are the lower and upper limits of the grid interval !----------------------------------------------------------------------------- jstart = 1 ngintv = ng - 1 do i = 1,ngintv !----------------------------------------------------------------------------- ! ... initalize: !----------------------------------------------------------------------------- area = 0._dp xgl = xg(i) xgu = xg(i+1) !----------------------------------------------------------------------------- ! ... discard data before the first grid interval and after the ! last grid interval ! for internal grid intervals, start calculating area by interpolating ! between the last point which lies in the previous interval and the ! first point inside the current interval !----------------------------------------------------------------------------- k = jstart if( k <= n-1 ) then !----------------------------------------------------------------------------- ! ... if both points are before the first grid, go to the next point !----------------------------------------------------------------------------- do if( x(k+1) <= xgl ) then jstart = k - 1 k = k+1 if( k <= n-1 ) then cycle else exit end if else exit end if end do !----------------------------------------------------------------------------- ! ... if the last point is beyond the end of the grid, ! complete and go to the next grid !----------------------------------------------------------------------------- do if( k <= n-1 .and. x(k) < xgu ) then jstart = k-1 !----------------------------------------------------------------------------- ! ... compute x-coordinates of increment !----------------------------------------------------------------------------- a1 = max( x(k),xgl ) a2 = min( x(k+1),xgu ) !----------------------------------------------------------------------------- ! ... if points coincide, contribution is zero !----------------------------------------------------------------------------- if( x(k+1) == x(k) ) then darea = 0._dp else slope = (y(k+1) - y(k))/(x(k+1) - x(k)) b1 = y(k) + slope*(a1 - x(k)) b2 = y(k) + slope*(a2 - x(k)) darea = .5*(a2 - a1)*(b2 + b1) end if !----------------------------------------------------------------------------- ! ... find the area under the trapezoid from a1 to a2 !----------------------------------------------------------------------------- area = area + darea k = k+1 cycle else exit end if end do end if !----------------------------------------------------------------------------- ! ... calculate the average y after summing the areas in the interval !----------------------------------------------------------------------------- yg(i) = area/(xgu - xgl) end do end subroutine inter2 subroutine inter_inti( ng, xg, n, x ) !----------------------------------------------------------------------------- ! ... initialization !----------------------------------------------------------------------------- implicit none !----------------------------------------------------------------------------- ! ... dummy arguments !----------------------------------------------------------------------------- integer, intent(in) :: ng, n real(dp), intent(in) :: xg(ng) real(dp), intent(in) :: x(n) !----------------------------------------------------------------------------- ! ... local variables !----------------------------------------------------------------------------- integer :: i, ii, iil, astat integer :: ndim(1) allocate( xi(ng), xcnt(ng-1), stat=astat ) if( astat /= 0 ) then write(*,*) 'inter_inti: failed to allocate wrk arrays; error = ',astat stop else xi(:) = 0 xcnt(:) = 0 end if iil = 1 do i = 1,ng do ii = iil,n-1 if( xg(i) < x(ii) ) then xi(i) = ii - 1 iil = ii exit end if end do end do nintervals = count( xi(:) /= 0 ) if( nintervals == 0 ) then write(*,*) 'inter_inti: wavelength grids do not overlap' stop else nintervals = nintervals - 1 end if xcnt(1:nintervals) = xi(2:nintervals+1) - xi(1:nintervals) + 1 ndim(:) = maxval( xcnt(1:nintervals) ) allocate( xfrac(ndim(1),nintervals),stat=astat ) if( astat /= 0 ) then write(*,*) 'inter_inti: failed to allocate wrk array; error = ',astat stop else xfrac(:,:) = 1._dp end if do i = 1,nintervals iil = xi(i) xfrac(1,i) = (min( x(iil+1),xg(i+1) ) - xg(i))/(x(iil+1) - x(iil)) if( xcnt(i) > 1 ) then iil = xi(i) + xcnt(i) - 1 xfrac(xcnt(i),i) = (xg(i+1) - x(iil))/(x(iil+1) - x(iil)) end if end do ! write(*,*) 'inter_inti: diagnostics; ng,n,nintervals = ',ng,n,nintervals ! write(*,'('' xi'')') ! write(*,'(10i4)') xi(1:nintervals+1) ! write(*,'('' xcnt'')') ! write(*,'(10i4)') xcnt(1:nintervals) ! write(*,'('' xfrac'')') ! do i = 1,nintervals ! write(*,'(1p,5e21.13)') xfrac(1:xcnt(i),i) ! end do end subroutine inter_inti subroutine inter3( ng, xg, yg, n, x, y ) !----------------------------------------------------------------------------- ! purpose: ! map input data given on a set of bins onto a different set of target ! bins. ! the input data are given on a set of bins (representing the integral ! of the input quantity over the range of each bin) and are being matched ! onto another set of bins (target grid). a typical example would be an ! input data set spcifying the extra-terrestrial flux on wavelength inter- ! vals, that has to be matched onto the working wavelength grid. ! the resulting area in a given bin of the target grid is calculated by ! simply adding all fractional areas of the input data that cover that ! particular target bin. ! some caution should be used near the endpoints of the grids. if the ! input data do not span the full range of the target grid, the area in ! the "missing" bins will be assumed to be zero. if the input data extend ! beyond the upper limit of the target grid, the user has the option to ! integrate the "overhang" data and fold the remaining area back into the ! last target bin. using this option is recommended when re-gridding ! vertical profiles that directly affect the total optical depth of the ! model atmosphere. !----------------------------------------------------------------------------- ! parameters: ! ng - integer, number of bins + 1 in the target grid (i) ! xg - real(dp), target grid (e.g. working wavelength grid); bin i (i) ! is defined as [xg(i),xg(i+1)] (i = 1..ng-1) ! yg - real(dp), y-data re-gridded onto xg; yg(i) specifies the (o) ! y-value for bin i (i = 1..ng-1) ! n - integer, number of bins + 1 in the input grid (i) ! x - real(dp), input grid (e.g. data wavelength grid); bin i is (i) ! defined as [x(i),x(i+1)] (i = 1..n-1) ! y - real(dp), input y-data on grid x; y(i) specifies the (i) ! y-value for bin i (i = 1..n-1) !----------------------------------------------------------------------------- implicit none !----------------------------------------------------------------------------- ! ... dummy arguments !----------------------------------------------------------------------------- integer, intent(in) :: n, ng real(dp), intent(in) :: xg(ng) real(dp), intent(in) :: x(n) real(dp), intent(in) :: y(n) real(dp), intent(out) :: yg(ng) !----------------------------------------------------------------------------- ! ... local variables !----------------------------------------------------------------------------- integer :: i, ii, iil !----------------------------------------------------------------------------- ! ... do interpolation !----------------------------------------------------------------------------- yg(:) = 0. do i = 1,nintervals iil = xi(i) ii = xcnt(i) if( ii == 1 ) then yg(i) = xfrac(1,i)*y(iil) else yg(i) = dot_product( xfrac(1:ii,i),y(iil:iil+ii-1) ) end if end do end subroutine inter3 subroutine calcoe( nz, c, xz, tt, adjin, adjcoe ) !----------------------------------------------------------------------------- ! parameters: ! adjcoe - real(dp), coross section adjust coefficients (in and out) ! c(5,28)-polynomal coef ! tt -nomarlized temperature !-----------------------------------------------------------------------------* implicit none !----------------------------------------------------------------------------- ! ... dummy arguments !----------------------------------------------------------------------------- integer, intent(in) :: nz real(dp), intent(in) :: adjin real(dp), intent(in) :: tt real(dp), intent(in) :: c(5) real(dp), intent(in) :: xz(nz) real(dp), intent(inout) :: adjcoe(nz) !----------------------------------------------------------------------------- ! ... local variables !----------------------------------------------------------------------------- integer :: k real(dp) :: x do k = 1,nz x = xz(k) adjcoe(k) = adjin * (1._dp + .01_dp*(c(1) + x*(c(2) + x*(c(3) + x*(c(4) + x*c(5)))))) end do end subroutine calcoe subroutine airmas( nz, z, zen, dsdh, nid, cz, & vcol, scol ) !----------------------------------------------------------------------------- ! purpose: ! calculate vertical and slant air columns, in spherical geometry, as a ! function of altitude. !----------------------------------------------------------------------------- ! parameters: ! nz - integer, number of specified altitude levels in the working (i) ! grid ! z - real(dp), specified altitude working grid (km) (i) ! zen - real(dp), solar zenith angle (degrees) (i) ! dsdh - real(dp), slant path of direct beam through each layer crossed (o) ! when travelling from the top of the atmosphere to layer i; ! dsdh(i,j), i = 0..nz-1, j = 1..nz-1 ! nid - integer, number of layers crossed by the direct beam when (o) ! travelling from the top of the atmosphere to layer i; ! nid(i), i = 0..nz-1 ! vcol - real(dp), output, vertical air column, molec cm-2, above level iz ! scol - real(dp), output, slant air column in direction of sun, above iz ! also in molec cm-2 !----------------------------------------------------------------------------- implicit none !----------------------------------------------------------------------------- ! ... dummy arguments !----------------------------------------------------------------------------- integer, intent(in) :: nz integer, intent(in) :: nid(0:nz-1) real(dp), intent(in) :: z(nz) real(dp), intent(in) :: zen real(dp), intent(in) :: dsdh(0:nz-1,nz-1) real(dp), intent(in) :: cz(nz) real(dp), intent(out) :: vcol(nz) real(dp), intent(out) :: scol(nz) !----------------------------------------------------------------------------- ! ... parameters !----------------------------------------------------------------------------- real(dp), parameter :: largest = huge(1.0_dp) !----------------------------------------------------------------------------- ! ... local variables !----------------------------------------------------------------------------- integer :: id, j real(dp) :: sum, ssum, vsum, ratio !----------------------------------------------------------------------------- ! ... calculate vertical and slant column from each level: work downward !----------------------------------------------------------------------------- vsum = 0._dp ssum = 0._dp do id = 1,nz-1 vsum = vsum + cz(nz-id) vcol(nz-id) = vsum sum = 0. if( nid(id) < 0 ) then sum = largest else !----------------------------------------------------------------------------- ! ... single pass layers: !----------------------------------------------------------------------------- do j = 1,min( nid(id),id ) sum = sum + cz(nz-j)*dsdh(id,j) end do !----------------------------------------------------------------------------- ! ... double pass layers: !----------------------------------------------------------------------------- do j = min( nid(id),id )+1,nid(id) sum = sum + 2._dp*cz(nz-j)*dsdh(id,j) end do end if scol(nz - id) = sum end do !----------------------------------------------------------------------------- ! ... special section to set scol(nz) !----------------------------------------------------------------------------- if( scol(nz-2) /= 0._dp ) then ratio = scol(nz-1)/scol(nz-2) scol(nz) = ratio * scol(nz-1) else scol(nz) = 0._dp end if end subroutine airmas END MODULE module_ftuv_subs