phys/module_ra_cam_support.F

   1 MODULE module_ra_cam_support
   2   use module_cam_support, only: endrun
   3   implicit none
   4
   5 #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
   6   include 'mpif.h'
   7 #else
   8   integer, parameter :: MPI_UNDEFINED = -1
   9 #endif
  10
  11       integer, parameter :: r8 = 8
  12       real(r8), parameter:: inf = 1.e20 ! CAM sets this differently in infnan.F90
  13       integer, parameter:: bigint = O'17777777777'           ! largest possible 32-bit integer
  14
  15       integer :: ixcldliq
  16       integer :: ixcldice
  17 !     integer :: levsiz    ! size of level dimension on dataset
  18       integer, parameter :: nbands = 2          ! Number of spectral bands
  19       integer, parameter :: naer_all = 12 + 1
  20       integer, parameter :: naer = 10 + 1
  21       integer, parameter :: bnd_nbr_LW=7
  22       integer, parameter :: ndstsz = 4    ! number of dust size bins
  23       integer :: idxSUL
  24       integer :: idxSSLT
  25       integer :: idxDUSTfirst
  26       integer :: idxCARBONfirst
  27       integer :: idxOCPHO
  28       integer :: idxBCPHO
  29       integer :: idxOCPHI
  30       integer :: idxBCPHI
  31       integer :: idxBG
  32       integer :: idxVOLC
  33       real, pointer :: ozmixin_save(:,:,:,:), lat_ozone_save(:), plev_ozone_save(:)
  34       integer :: levsiz_ozone_save=-1
  35       logical :: have_ozone=.false.
  36   integer :: mxaerl                            ! Maximum level of background aerosol
  37
  38 ! indices to sections of array that represent
  39 ! groups of aerosols
  40
  41   integer, parameter :: &
  42       numDUST         = 4, &
  43       numCARBON      = 4
  44
  45 ! portion of each species group to use in computation
  46 ! of relative radiative forcing.
  47
  48   real(r8) :: sulscl_rf  = 0._r8 !
  49   real(r8) :: carscl_rf  = 0._r8
  50   real(r8) :: ssltscl_rf = 0._r8
  51   real(r8) :: dustscl_rf = 0._r8
  52   real(r8) :: bgscl_rf   = 0._r8
  53   real(r8) :: volcscl_rf = 0._r8
  54
  55 ! "background" aerosol species mmr.
  56   real(r8) :: tauback = 0._r8
  57
  58 ! portion of each species group to use in computation
  59 ! of aerosol forcing in driving the climate
  60   real(r8) :: sulscl  = 1._r8
  61   real(r8) :: carscl  = 1._r8
  62   real(r8) :: ssltscl = 1._r8
  63   real(r8) :: dustscl = 1._r8
  64   real(r8) :: volcscl = 1._r8
  65
  66 !From volcrad.F90 module
  67      integer, parameter :: idx_LW_0500_0650=3
  68      integer, parameter :: idx_LW_0650_0800=4
  69      integer, parameter :: idx_LW_0800_1000=5
  70      integer, parameter :: idx_LW_1000_1200=6
  71      integer, parameter :: idx_LW_1200_2000=7
  72
  73 ! First two values represent the overlap of volcanics with the non-window
  74 ! (0-800, 1200-2200 cm^-1) and window (800-1200 cm^-1) regions.|  Coefficients
  75 ! were derived using crm_volc_minimize.pro with spectral flux optimization
  76 ! on first iteration, total heating rate on subsequent iterations (2-9).
  77 ! Five profiles for HLS, HLW, MLS, MLW, and TRO conditions were given equal
  78 ! weight.  RMS heating rate errors for a visible stratospheric optical
  79 ! depth of 1.0 are 0.02948 K/day.
  80 !
  81       real(r8) :: abs_cff_mss_aer(bnd_nbr_LW) = &
  82          (/ 70.257384, 285.282943, &
  83          1.0273851e+02, 6.3073303e+01, 1.2039569e+02, &
  84          3.6343643e+02, 2.7138528e+02 /)
  85
  86 !From radae.F90 module
  87       real(r8), parameter:: min_tp_h2o = 160.0        ! min T_p for pre-calculated abs/emis
  88       real(r8), parameter:: max_tp_h2o = 349.999999   ! max T_p for pre-calculated abs/emis
  89       real(r8), parameter:: dtp_h2o = 21.111111111111 ! difference in adjacent elements of tp_h2o
  90       real(r8), parameter:: min_te_h2o = -120.0       ! min T_e-T_p for pre-calculated abs/emis
  91       real(r8), parameter:: max_te_h2o = 79.999999    ! max T_e-T_p for pre-calculated abs/emis
  92       real(r8), parameter:: dte_h2o  = 10.0           ! difference in adjacent elements of te_h2o
  93       real(r8), parameter:: min_rh_h2o = 0.0          ! min RH for pre-calculated abs/emis
  94       real(r8), parameter:: max_rh_h2o = 1.19999999   ! max RH for pre-calculated abs/emis
  95       real(r8), parameter:: drh_h2o = 0.2             ! difference in adjacent elements of RH
  96       real(r8), parameter:: min_lu_h2o = -8.0         ! min log_10(U) for pre-calculated abs/emis
  97       real(r8), parameter:: min_u_h2o  = 1.0e-8       ! min pressure-weighted path-length
  98       real(r8), parameter:: max_lu_h2o =  3.9999999   ! max log_10(U) for pre-calculated abs/emis
  99       real(r8), parameter:: dlu_h2o  = 0.5            ! difference in adjacent elements of lu_h2o
 100       real(r8), parameter:: min_lp_h2o = -3.0         ! min log_10(P) for pre-calculated abs/emis
 101       real(r8), parameter:: min_p_h2o = 1.0e-3        ! min log_10(P) for pre-calculated abs/emis
 102       real(r8), parameter:: max_lp_h2o = -0.0000001   ! max log_10(P) for pre-calculated abs/emis
 103       real(r8), parameter:: dlp_h2o = 0.3333333333333 ! difference in adjacent elements of lp_h2o
 104       integer, parameter :: n_u = 25   ! Number of U in abs/emis tables
 105       integer, parameter :: n_p = 10   ! Number of P in abs/emis tables
 106       integer, parameter :: n_tp = 10  ! Number of T_p in abs/emis tables
 107       integer, parameter :: n_te = 21  ! Number of T_e in abs/emis tables
 108       integer, parameter :: n_rh = 7   ! Number of RH in abs/emis tables
 109       real(r8):: c16,c17,c26,c27,c28,c29,c30,c31
 110       real(r8):: fwcoef      ! Farwing correction constant
 111       real(r8):: fwc1,fwc2   ! Farwing correction constants
 112       real(r8):: fc1         ! Farwing correction constant
 113       real(r8):: amco2 ! Molecular weight of co2   (g/mol)
 114       real(r8):: amd   ! Molecular weight of dry air (g/mol)
 115       real(r8):: p0    ! Standard pressure (dynes/cm**2)
 116
 117 ! These are now allocatable. JM 20090612
 118   real(r8), allocatable, dimension(:,:,:,:,:)  :: ah2onw   ! (n_p, n_tp, n_u, n_te, n_rh)   ! absorptivity (non-window)
 119   real(r8), allocatable, dimension(:,:,:,:,:)  :: eh2onw   ! (n_p, n_tp, n_u, n_te, n_rh)   ! emissivity   (non-window)
 120   real(r8), allocatable, dimension(:,:,:,:,:)  :: ah2ow    ! (n_p, n_tp, n_u, n_te, n_rh)    ! absorptivity (window, for adjacent layers)
 121   real(r8), allocatable, dimension(:,:,:,:,:)  :: cn_ah2ow ! (n_p, n_tp, n_u, n_te, n_rh)    ! continuum transmission for absorptivity (window)
 122   real(r8), allocatable, dimension(:,:,:,:,:)  :: cn_eh2ow ! (n_p, n_tp, n_u, n_te, n_rh)    ! continuum transmission for emissivity   (window)
 123   real(r8), allocatable, dimension(:,:,:,:,:)  :: ln_ah2ow ! (n_p, n_tp, n_u, n_te, n_rh)    ! line-only transmission for absorptivity (window)
 124   real(r8), allocatable, dimension(:,:,:,:,:)  :: ln_eh2ow ! (n_p, n_tp, n_u, n_te, n_rh)    ! line-only transmission for emissivity   (window)
 125
 126 !
 127 ! Constant coefficients for water vapor overlap with trace gases.
 128 ! Reference: Ramanathan, V. and  P.Downey, 1986: A Nonisothermal
 129 !            Emissivity and Absorptivity Formulation for Water Vapor
 130 !            Journal of Geophysical Research, vol. 91., D8, pp 8649-8666
 131 !
 132   real(r8):: coefh(2,4) = reshape(  &
 133          (/ (/5.46557e+01,-7.30387e-02/), &
 134             (/1.09311e+02,-1.46077e-01/), &
 135             (/5.11479e+01,-6.82615e-02/), &
 136             (/1.02296e+02,-1.36523e-01/) /), (/2,4/) )
 137 !
 138   real(r8):: coefj(3,2) = reshape( &
 139             (/ (/2.82096e-02,2.47836e-04,1.16904e-06/), &
 140                (/9.27379e-02,8.04454e-04,6.88844e-06/) /), (/3,2/) )
 141 !
 142   real(r8):: coefk(3,2) = reshape( &
 143             (/ (/2.48852e-01,2.09667e-03,2.60377e-06/) , &
 144                (/1.03594e+00,6.58620e-03,4.04456e-06/) /), (/3,2/) )
 145
 146   integer, parameter :: ntemp = 192 ! Number of temperatures in H2O sat. table for Tp
 147   real(r8) :: estblh2o(0:ntemp)       ! saturation vapor pressure for H2O for Tp rang
 148   integer, parameter :: o_fa = 6   ! Degree+1 of poly of T_e for absorptivity as U->inf.
 149   integer, parameter :: o_fe = 6   ! Degree+1 of poly of T_e for emissivity as U->inf.
 150
 151 !-----------------------------------------------------------------------------
 152 ! Data for f in C/H/E fit -- value of A and E as U->infinity
 153 ! New C/LT/E fit (Hitran 2K, CKD 2.4) -- no change
 154 !     These values are determined by integrals of Planck functions or
 155 !     derivatives of Planck functions only.
 156 !-----------------------------------------------------------------------------
 157 !
 158 ! fa/fe coefficients for 2 bands (0-800 & 1200-2200, 800-1200 cm^-1)
 159 !
 160 ! Coefficients of polynomial for f_a in T_e
 161 !
 162   real(r8), parameter:: fat(o_fa,nbands) = reshape( (/ &
 163        (/-1.06665373E-01,  2.90617375E-02, -2.70642049E-04,   &   ! 0-800&1200-2200 cm^-1
 164           1.07595511E-06, -1.97419681E-09,  1.37763374E-12/), &   !   0-800&1200-2200 cm^-1
 165        (/ 1.10666537E+00, -2.90617375E-02,  2.70642049E-04,   &   ! 800-1200 cm^-1
 166          -1.07595511E-06,  1.97419681E-09, -1.37763374E-12/) /) & !   800-1200 cm^-1
 167        , (/o_fa,nbands/) )
 168 !
 169 ! Coefficients of polynomial for f_e in T_e
 170 !
 171   real(r8), parameter:: fet(o_fe,nbands) = reshape( (/ &
 172       (/3.46148163E-01,  1.51240299E-02, -1.21846479E-04,   &   ! 0-800&1200-2200 cm^-1
 173         4.04970123E-07, -6.15368936E-10,  3.52415071E-13/), &   !   0-800&1200-2200 cm^-1
 174       (/6.53851837E-01, -1.51240299E-02,  1.21846479E-04,   &   ! 800-1200 cm^-1
 175        -4.04970123E-07,  6.15368936E-10, -3.52415071E-13/) /) & !   800-1200 cm^-1
 176       , (/o_fa,nbands/) )
 177
 178
 179       real(r8) ::  gravit     ! Acceleration of gravity (cgs)
 180       real(r8) ::  rga        ! 1./gravit
 181       real(r8) ::  gravmks    ! Acceleration of gravity (mks)
 182       real(r8) ::  cpair      ! Specific heat of dry air
 183       real(r8) ::  epsilo     ! Ratio of mol. wght of H2O to dry air
 184       real(r8) ::  epsqs      ! Ratio of mol. wght of H2O to dry air
 185       real(r8) ::  sslp       ! Standard sea-level pressure
 186       real(r8) ::  stebol     ! Stefan-Boltzmann's constant
 187       real(r8) ::  rgsslp     ! 0.5/(gravit*sslp)
 188       real(r8) ::  dpfo3      ! Voigt correction factor for O3
 189       real(r8) ::  dpfco2     ! Voigt correction factor for CO2
 190       real(r8) ::  dayspy     ! Number of days per 1 year
 191       real(r8) ::  pie        ! 3.14.....
 192       real(r8) ::  mwdry      ! molecular weight dry air ~ kg/kmole (shr_const_mwdair)
 193       real(r8) ::  scon       ! solar constant (not used in WRF)
 194       real(r8) ::  co2mmr
 195 real(r8) ::   mwco2              ! molecular weight of carbon dioxide
 196 real(r8) ::   mwh2o              ! molecular weight water vapor (shr_const_mwwv)
 197 real(r8) ::   mwch4              ! molecular weight ch4
 198 real(r8) ::   mwn2o              ! molecular weight n2o
 199 real(r8) ::   mwf11              ! molecular weight cfc11
 200 real(r8) ::   mwf12              ! molecular weight cfc12
 201 real(r8) ::   cappa              ! R/Cp
 202 real(r8) ::   rair               ! Gas constant for dry air (J/K/kg)
 203 real(r8) ::   tmelt              ! freezing T of fresh water ~ K
 204 real(r8) ::   r_universal        ! Universal gas constant ~ J/K/kmole
 205 real(r8) ::   latvap             ! latent heat of evaporation ~ J/kg
 206 real(r8) ::   latice             ! latent heat of fusion ~ J/kg
 207 real(r8) ::   zvir               ! R_V/R_D - 1.
 208   integer plenest  ! length of saturation vapor pressure table
 209   parameter (plenest=250)
 210 !
 211 ! Table of saturation vapor pressure values es from tmin degrees
 212 ! to tmax+1 degrees k in one degree increments.  ttrice defines the
 213 ! transition region where es is a combination of ice & water values
 214 !
 215 real(r8) estbl(plenest)      ! table values of saturation vapor pressure
 216 real(r8) tmin       ! min temperature (K) for table
 217 real(r8) tmax       ! max temperature (K) for table
 218 real(r8) pcf(6)     ! polynomial coeffs -> es transition water to ice
 219 !real(r8), allocatable :: pin(:)           ! ozone pressure level (levsiz)
 220 !real(r8), allocatable :: ozmix(:,:,:)     ! mixing ratio
 221 !real(r8), allocatable, target :: abstot_3d(:,:,:,:) ! Non-adjacent layer absorptivites
 222 !real(r8), allocatable, target :: absnxt_3d(:,:,:,:) ! Nearest layer absorptivities
 223 !real(r8), allocatable, target :: emstot_3d(:,:,:)   ! Total emissivity
 224
 225 !From aer_optics.F90 module
 226 integer, parameter :: idxVIS = 8     ! index to visible band
 227 integer, parameter :: nrh = 1000   ! number of relative humidity values for look-up-table
 228 integer, parameter :: nspint = 19   ! number of spectral intervals
 229
 230 ! These are now allocatable,  JM 20090612
 231 real(r8), allocatable, dimension(:,:) :: ksul    ! (nrh, nspint)    ! sulfate specific extinction  ( m^2 g-1 )
 232 real(r8), allocatable, dimension(:,:) :: wsul    ! (nrh, nspint)    ! sulfate single scattering albedo
 233 real(r8), allocatable, dimension(:,:) :: gsul    ! (nrh, nspint)    ! sulfate asymmetry parameter
 234 real(r8), allocatable, dimension(:,:) :: ksslt   ! (nrh, nspint)   ! sea-salt specific extinction  ( m^2 g-1 )
 235 real(r8), allocatable, dimension(:,:) :: wsslt   ! (nrh, nspint)   ! sea-salt single scattering albedo
 236 real(r8), allocatable, dimension(:,:) :: gsslt   ! (nrh, nspint)   ! sea-salt asymmetry parameter
 237 real(r8), allocatable, dimension(:,:) :: kcphil  ! (nrh, nspint)  ! hydrophilic carbon specific extinction  ( m^2 g-1 )
 238 real(r8), allocatable, dimension(:,:) :: wcphil  ! (nrh, nspint)  ! hydrophilic carbon single scattering albedo
 239 real(r8), allocatable, dimension(:,:) :: gcphil  ! (nrh, nspint)  ! hydrophilic carbon asymmetry parameter
 240
 241 real(r8) :: kbg(nspint)          ! background specific extinction  ( m^2 g-1 )
 242 real(r8) :: wbg(nspint)          ! background single scattering albedo
 243 real(r8) :: gbg(nspint)          ! background asymmetry parameter
 244 real(r8) :: kcphob(nspint)       ! hydrophobic carbon specific extinction  ( m^2 g-1 )
 245 real(r8) :: wcphob(nspint)       ! hydrophobic carbon single scattering albedo
 246 real(r8) :: gcphob(nspint)       ! hydrophobic carbon asymmetry parameter
 247 real(r8) :: kcb(nspint)          ! black carbon specific extinction  ( m^2 g-1 )
 248 real(r8) :: wcb(nspint)          ! black carbon single scattering albedo
 249 real(r8) :: gcb(nspint)          ! black carbon asymmetry parameter
 250 real(r8) :: kvolc(nspint)        ! volcanic specific extinction  ( m^2 g-1)
 251 real(r8) :: wvolc(nspint)        ! volcanic single scattering albedo
 252 real(r8) :: gvolc(nspint)        ! volcanic asymmetry parameter
 253
 254 real(r8) :: kdst(ndstsz, nspint) ! dust specific extinction  ( m^2 g-1 )
 255 real(r8) :: wdst(ndstsz, nspint) ! dust single scattering albedo
 256 real(r8) :: gdst(ndstsz, nspint) ! dust asymmetry parameter
 257 !
 258 !From comozp.F90 module
 259       real(r8) cplos    ! constant for ozone path length integral
 260       real(r8) cplol    ! constant for ozone path length integral
 261
 262    real(r8) :: co2vmr = 3.550e-4         ! co2   volume mixing ratio
 263    real(r8) :: n2ovmr = 0.311e-6         ! n2o   volume mixing ratio
 264    real(r8) :: ch4vmr = 1.714e-6         ! ch4   volume mixing ratio
 265    real(r8) :: f11vmr = 0.280e-9         ! cfc11 volume mixing ratio
 266    real(r8) :: f12vmr = 0.503e-9         ! cfc12 volume mixing ratio
 267
 268    integer, parameter  :: cyr = 233  ! number of years of co2 data
 269
 270    integer  :: yrdata(cyr) = &
 271  (/ 1869, 1870, 1871, 1872, 1873, 1874, 1875, &
 272     1876, 1877, 1878, 1879, 1880, 1881, 1882, &
 273     1883, 1884, 1885, 1886, 1887, 1888, 1889, &
 274     1890, 1891, 1892, 1893, 1894, 1895, 1896, &
 275     1897, 1898, 1899, 1900, 1901, 1902, 1903, &
 276     1904, 1905, 1906, 1907, 1908, 1909, 1910, &
 277     1911, 1912, 1913, 1914, 1915, 1916, 1917, &
 278     1918, 1919, 1920, 1921, 1922, 1923, 1924, &
 279     1925, 1926, 1927, 1928, 1929, 1930, 1931, &
 280     1932, 1933, 1934, 1935, 1936, 1937, 1938, &
 281     1939, 1940, 1941, 1942, 1943, 1944, 1945, &
 282     1946, 1947, 1948, 1949, 1950, 1951, 1952, &
 283     1953, 1954, 1955, 1956, 1957, 1958, 1959, &
 284     1960, 1961, 1962, 1963, 1964, 1965, 1966, &
 285     1967, 1968, 1969, 1970, 1971, 1972, 1973, &
 286     1974, 1975, 1976, 1977, 1978, 1979, 1980, &
 287     1981, 1982, 1983, 1984, 1985, 1986, 1987, &
 288     1988, 1989, 1990, 1991, 1992, 1993, 1994, &
 289     1995, 1996, 1997, 1998, 1999, 2000, 2001, &
 290     2002, 2003, 2004, 2005, 2006, 2007, 2008, &
 291     2009, 2010, 2011, 2012, 2013, 2014, 2015, &
 292     2016, 2017, 2018, 2019, 2020, 2021, 2022, &
 293     2023, 2024, 2025, 2026, 2027, 2028, 2029, &
 294     2030, 2031, 2032, 2033, 2034, 2035, 2036, &
 295     2037, 2038, 2039, 2040, 2041, 2042, 2043, &
 296     2044, 2045, 2046, 2047, 2048, 2049, 2050, &
 297     2051, 2052, 2053, 2054, 2055, 2056, 2057, &
 298     2058, 2059, 2060, 2061, 2062, 2063, 2064, &
 299     2065, 2066, 2067, 2068, 2069, 2070, 2071, &
 300     2072, 2073, 2074, 2075, 2076, 2077, 2078, &
 301     2079, 2080, 2081, 2082, 2083, 2084, 2085, &
 302     2086, 2087, 2088, 2089, 2090, 2091, 2092, &
 303     2093, 2094, 2095, 2096, 2097, 2098, 2099, &
 304     2100, 2101                               /)
 305
 306 ! A2 future scenario
 307     real(r8)  :: co2(cyr) = &
 308  (/ 289.263, 289.263, 289.416, 289.577, 289.745, 289.919, 290.102, &
 309     290.293, 290.491, 290.696, 290.909, 291.129, 291.355, 291.587, 291.824, &
 310     292.066, 292.313, 292.563, 292.815, 293.071, 293.328, 293.586, 293.843, &
 311     294.098, 294.35, 294.598, 294.842, 295.082, 295.32, 295.558, 295.797,   &
 312     296.038, 296.284, 296.535, 296.794, 297.062, 297.338, 297.62, 297.91,   &
 313     298.204, 298.504, 298.806, 299.111, 299.419, 299.729, 300.04, 300.352,  &
 314     300.666, 300.98, 301.294, 301.608, 301.923, 302.237, 302.551, 302.863,  &
 315     303.172, 303.478, 303.779, 304.075, 304.366, 304.651, 304.93, 305.206,  &
 316     305.478, 305.746, 306.013, 306.28, 306.546, 306.815, 307.087, 307.365,  &
 317     307.65, 307.943, 308.246, 308.56, 308.887, 309.228, 309.584, 309.956,   &
 318     310.344, 310.749, 311.172, 311.614, 312.077, 312.561, 313.068, 313.599, &
 319     314.154, 314.737, 315.347, 315.984, 316.646, 317.328, 318.026, 318.742, &
 320     319.489, 320.282, 321.133, 322.045, 323.021, 324.06, 325.155, 326.299,  &
 321     327.484, 328.698, 329.933, 331.194, 332.499, 333.854, 335.254, 336.69,  &
 322     338.15, 339.628, 341.125, 342.65, 344.206, 345.797, 347.397, 348.98,    &
 323     350.551, 352.1, 354.3637, 355.7772, 357.1601, 358.5306, 359.9046,       &
 324     361.4157, 363.0445, 364.7761, 366.6064, 368.5322, 370.534, 372.5798,    &
 325     374.6564, 376.7656, 378.9087, 381.0864, 383.2994, 385.548, 387.8326,    &
 326     390.1536, 392.523, 394.9625, 397.4806, 400.075, 402.7444, 405.4875,     &
 327     408.3035, 411.1918, 414.1518, 417.1831, 420.2806, 423.4355, 426.6442,   &
 328     429.9076, 433.2261, 436.6002, 440.0303, 443.5168, 447.06, 450.6603,     &
 329     454.3059, 457.9756, 461.6612, 465.3649, 469.0886, 472.8335, 476.6008,   &
 330     480.3916, 484.2069, 488.0473, 491.9184, 495.8295, 499.7849, 503.7843,   &
 331     507.8278, 511.9155, 516.0476, 520.2243, 524.4459, 528.7127, 533.0213,   &
 332     537.3655, 541.7429, 546.1544, 550.6005, 555.0819, 559.5991, 564.1525,   &
 333     568.7429, 573.3701, 578.0399, 582.7611, 587.5379, 592.3701, 597.2572,   &
 334     602.1997, 607.1975, 612.2507, 617.3596, 622.524, 627.7528, 633.0616,    &
 335     638.457, 643.9384, 649.505, 655.1568, 660.8936, 666.7153, 672.6219,     &
 336     678.6133, 684.6945, 690.8745, 697.1569, 703.5416, 710.0284, 716.6172,   &
 337     723.308, 730.1008, 736.9958, 743.993, 751.0975, 758.3183, 765.6594,     &
 338     773.1207, 780.702, 788.4033, 796.2249, 804.1667, 812.2289, 820.4118,    &
 339     828.6444, 828.6444 /)
 340
 341       integer  :: ntoplw      ! top level to solve for longwave cooling (WRF sets this to 1 for model top below 10 mb)
 342
 343       logical :: masterproc = .true.
 344       logical :: ozncyc            ! true => cycle ozone dataset
 345 !     logical :: dosw              ! True => shortwave calculation this timestep
 346 !     logical :: dolw              ! True => longwave calculation this timestep
 347       logical :: indirect          ! True => include indirect radiative effects of sulfate aerosols
 348 !     logical :: doabsems          ! True => abs/emiss calculation this timestep
 349       logical :: radforce   = .false.          ! True => calculate aerosol shortwave forcing
 350       logical :: trace_gas=.false.             ! set true for chemistry
 351       logical :: strat_volcanic   = .false.    ! True => volcanic aerosol mass available
 352
 353     real(r8) retab(95)
 354     !
 355     !       Tabulated values of re(T) in the temperature interval
 356     !       180 K -- 274 K; hexagonal columns assumed:
 357     !
 358     data retab /                                                &
 359          5.92779, 6.26422, 6.61973, 6.99539, 7.39234,   &
 360          7.81177, 8.25496, 8.72323, 9.21800, 9.74075, 10.2930,  &
 361          10.8765, 11.4929, 12.1440, 12.8317, 13.5581, 14.2319,  &
 362          15.0351, 15.8799, 16.7674, 17.6986, 18.6744, 19.6955,  &
 363          20.7623, 21.8757, 23.0364, 24.2452, 25.5034, 26.8125,  &
 364          27.7895, 28.6450, 29.4167, 30.1088, 30.7306, 31.2943,  &
 365          31.8151, 32.3077, 32.7870, 33.2657, 33.7540, 34.2601,  &
 366          34.7892, 35.3442, 35.9255, 36.5316, 37.1602, 37.8078,  &
 367          38.4720, 39.1508, 39.8442, 40.5552, 41.2912, 42.0635,  &
 368          42.8876, 43.7863, 44.7853, 45.9170, 47.2165, 48.7221,  &
 369          50.4710, 52.4980, 54.8315, 57.4898, 60.4785, 63.7898,  &
 370          65.5604, 71.2885, 75.4113, 79.7368, 84.2351, 88.8833,  &
 371          93.6658, 98.5739, 103.603, 108.752, 114.025, 119.424,  &
 372          124.954, 130.630, 136.457, 142.446, 148.608, 154.956,  &
 373          161.503, 168.262, 175.248, 182.473, 189.952, 197.699,  &
 374          205.728, 214.055, 222.694, 231.661, 240.971, 250.639/
 375     !
 376     save retab
 377 contains
 378
 379
 380
 381 subroutine sortarray(n, ain, indxa)
 382 !-----------------------------------------------
 383 !
 384 ! Purpose:
 385 !       Sort an array
 386 ! Alogrithm:
 387 !       Based on Shell's sorting method.
 388 !
 389 ! Author: T. Craig
 390 !-----------------------------------------------
 391 !  use shr_kind_mod, only: r8 => shr_kind_r8
 392    implicit none
 393 !
 394 !  Arguments
 395 !
 396    integer , intent(in) :: n             ! total number of elements
 397    integer , intent(inout) :: indxa(n)   ! array of integers
 398    real(r8), intent(inout) :: ain(n)     ! array to sort
 399 !
 400 !  local variables
 401 !
 402    integer :: i, j                ! Loop indices
 403    integer :: ni                  ! Starting increment
 404    integer :: itmp                ! Temporary index
 405    real(r8):: atmp                ! Temporary value to swap
 406
 407    ni = 1
 408    do while(.TRUE.)
 409       ni = 3*ni + 1
 410       if (ni <= n) cycle
 411       exit
 412    end do
 413
 414    do while(.TRUE.)
 415       ni = ni/3
 416       do i = ni + 1, n
 417          atmp = ain(i)
 418          itmp = indxa(i)
 419          j = i
 420          do while(.TRUE.)
 421             if (ain(j-ni) <= atmp) exit
 422             ain(j) = ain(j-ni)
 423             indxa(j) = indxa(j-ni)
 424             j = j - ni
 425             if (j > ni) cycle
 426             exit
 427          end do
 428          ain(j) = atmp
 429          indxa(j) = itmp
 430       end do
 431       if (ni > 1) cycle
 432       exit
 433    end do
 434    return
 435
 436 end subroutine sortarray
 437 subroutine trcab(lchnk   ,ncol    ,pcols, pverp,               &
 438                  k1      ,k2      ,ucfc11  ,ucfc12  ,un2o0   , &
 439                  un2o1   ,uch4    ,uco211  ,uco212  ,uco213  , &
 440                  uco221  ,uco222  ,uco223  ,bn2o0   ,bn2o1   , &
 441                  bch4    ,to3co2  ,pnm     ,dw      ,pnew    , &
 442                  s2c     ,uptype  ,dplh2o  ,abplnk1 ,tco2    , &
 443                  th2o    ,to3     ,abstrc  , &
 444                  aer_trn_ttl)
 445 !-----------------------------------------------------------------------
 446 !
 447 ! Purpose:
 448 ! Calculate absorptivity for non nearest layers for CH4, N2O, CFC11 and
 449 ! CFC12.
 450 !
 451 ! Method:
 452 ! See CCM3 description for equations.
 453 !
 454 ! Author: J. Kiehl
 455 !
 456 !-----------------------------------------------------------------------
 457 !  use shr_kind_mod, only: r8 => shr_kind_r8
 458 !  use ppgrid
 459 !  use volcrad
 460
 461    implicit none
 462
 463 !------------------------------Arguments--------------------------------
 464 !
 465 ! Input arguments
 466 !
 467    integer, intent(in) :: lchnk                    ! chunk identifier
 468    integer, intent(in) :: ncol                     ! number of atmospheric columns
 469    integer, intent(in) :: pcols, pverp
 470    integer, intent(in) :: k1,k2                    ! level indices
 471 !
 472    real(r8), intent(in) :: to3co2(pcols)           ! pressure weighted temperature
 473    real(r8), intent(in) :: pnm(pcols,pverp)        ! interface pressures
 474    real(r8), intent(in) :: ucfc11(pcols,pverp)     ! CFC11 path length
 475    real(r8), intent(in) :: ucfc12(pcols,pverp)     ! CFC12 path length
 476    real(r8), intent(in) :: un2o0(pcols,pverp)      ! N2O path length
 477 !
 478    real(r8), intent(in) :: un2o1(pcols,pverp)      ! N2O path length (hot band)
 479    real(r8), intent(in) :: uch4(pcols,pverp)       ! CH4 path length
 480    real(r8), intent(in) :: uco211(pcols,pverp)     ! CO2 9.4 micron band path length
 481    real(r8), intent(in) :: uco212(pcols,pverp)     ! CO2 9.4 micron band path length
 482    real(r8), intent(in) :: uco213(pcols,pverp)     ! CO2 9.4 micron band path length
 483 !
 484    real(r8), intent(in) :: uco221(pcols,pverp)     ! CO2 10.4 micron band path length
 485    real(r8), intent(in) :: uco222(pcols,pverp)     ! CO2 10.4 micron band path length
 486    real(r8), intent(in) :: uco223(pcols,pverp)     ! CO2 10.4 micron band path length
 487    real(r8), intent(in) :: bn2o0(pcols,pverp)      ! pressure factor for n2o
 488    real(r8), intent(in) :: bn2o1(pcols,pverp)      ! pressure factor for n2o
 489 !
 490    real(r8), intent(in) :: bch4(pcols,pverp)       ! pressure factor for ch4
 491    real(r8), intent(in) :: dw(pcols)               ! h2o path length
 492    real(r8), intent(in) :: pnew(pcols)             ! pressure
 493    real(r8), intent(in) :: s2c(pcols,pverp)        ! continuum path length
 494    real(r8), intent(in) :: uptype(pcols,pverp)     ! p-type h2o path length
 495 !
 496    real(r8), intent(in) :: dplh2o(pcols)           ! p squared h2o path length
 497    real(r8), intent(in) :: abplnk1(14,pcols,pverp) ! Planck factor
 498    real(r8), intent(in) :: tco2(pcols)             ! co2 transmission factor
 499    real(r8), intent(in) :: th2o(pcols)             ! h2o transmission factor
 500    real(r8), intent(in) :: to3(pcols)              ! o3 transmission factor
 501
 502    real(r8), intent(in) :: aer_trn_ttl(pcols,pverp,pverp,bnd_nbr_LW) ! aer trn.
 503
 504 !
 505 !  Output Arguments
 506 !
 507    real(r8), intent(out) :: abstrc(pcols)           ! total trace gas absorptivity
 508 !
 509 !--------------------------Local Variables------------------------------
 510 !
 511    integer  i,l                     ! loop counters
 512
 513    real(r8) sqti(pcols)             ! square root of mean temp
 514    real(r8) du1                     ! cfc11 path length
 515    real(r8) du2                     ! cfc12 path length
 516    real(r8) acfc1                   ! cfc11 absorptivity 798 cm-1
 517    real(r8) acfc2                   ! cfc11 absorptivity 846 cm-1
 518 !
 519    real(r8) acfc3                   ! cfc11 absorptivity 933 cm-1
 520    real(r8) acfc4                   ! cfc11 absorptivity 1085 cm-1
 521    real(r8) acfc5                   ! cfc12 absorptivity 889 cm-1
 522    real(r8) acfc6                   ! cfc12 absorptivity 923 cm-1
 523    real(r8) acfc7                   ! cfc12 absorptivity 1102 cm-1
 524 !
 525    real(r8) acfc8                   ! cfc12 absorptivity 1161 cm-1
 526    real(r8) du01                    ! n2o path length
 527    real(r8) dbeta01                 ! n2o pressure factor
 528    real(r8) dbeta11                 !         "
 529    real(r8) an2o1                   ! absorptivity of 1285 cm-1 n2o band
 530 !
 531    real(r8) du02                    ! n2o path length
 532    real(r8) dbeta02                 ! n2o pressure factor
 533    real(r8) an2o2                   ! absorptivity of 589 cm-1 n2o band
 534    real(r8) du03                    ! n2o path length
 535    real(r8) dbeta03                 ! n2o pressure factor
 536 !
 537    real(r8) an2o3                   ! absorptivity of 1168 cm-1 n2o band
 538    real(r8) duch4                   ! ch4 path length
 539    real(r8) dbetac                  ! ch4 pressure factor
 540    real(r8) ach4                    ! absorptivity of 1306 cm-1 ch4 band
 541    real(r8) du11                    ! co2 path length
 542 !
 543    real(r8) du12                    !       "
 544    real(r8) du13                    !       "
 545    real(r8) dbetc1                  ! co2 pressure factor
 546    real(r8) dbetc2                  ! co2 pressure factor
 547    real(r8) aco21                   ! absorptivity of 1064 cm-1 band
 548 !
 549    real(r8) du21                    ! co2 path length
 550    real(r8) du22                    !       "
 551    real(r8) du23                    !       "
 552    real(r8) aco22                   ! absorptivity of 961 cm-1 band
 553    real(r8) tt(pcols)               ! temp. factor for h2o overlap factor
 554 !
 555    real(r8) psi1                    !                 "
 556    real(r8) phi1                    !                 "
 557    real(r8) p1                      ! h2o overlap factor
 558    real(r8) w1                      !        "
 559    real(r8) ds2c(pcols)             ! continuum path length
 560 !
 561    real(r8) duptyp(pcols)           ! p-type path length
 562    real(r8) tw(pcols,6)             ! h2o transmission factor
 563    real(r8) g1(6)                   !         "
 564    real(r8) g2(6)                   !         "
 565    real(r8) g3(6)                   !         "
 566 !
 567    real(r8) g4(6)                   !         "
 568    real(r8) ab(6)                   ! h2o temp. factor
 569    real(r8) bb(6)                   !         "
 570    real(r8) abp(6)                  !         "
 571    real(r8) bbp(6)                  !         "
 572 !
 573    real(r8) tcfc3                   ! transmission for cfc11 band
 574    real(r8) tcfc4                   ! transmission for cfc11 band
 575    real(r8) tcfc6                   ! transmission for cfc12 band
 576    real(r8) tcfc7                   ! transmission for cfc12 band
 577    real(r8) tcfc8                   ! transmission for cfc12 band
 578 !
 579    real(r8) tlw                     ! h2o transmission
 580    real(r8) tch4                    ! ch4 transmission
 581 !
 582 !--------------------------Data Statements------------------------------
 583 !
 584    data g1 /0.0468556,0.0397454,0.0407664,0.0304380,0.0540398,0.0321962/
 585    data g2 /14.4832,4.30242,5.23523,3.25342,0.698935,16.5599/
 586    data g3 /26.1898,18.4476,15.3633,12.1927,9.14992,8.07092/
 587    data g4 /0.0261782,0.0369516,0.0307266,0.0243854,0.0182932,0.0161418/
 588    data ab /3.0857e-2,2.3524e-2,1.7310e-2,2.6661e-2,2.8074e-2,2.2915e-2/
 589    data bb /-1.3512e-4,-6.8320e-5,-3.2609e-5,-1.0228e-5,-9.5743e-5,-1.0304e-4/
 590    data abp/2.9129e-2,2.4101e-2,1.9821e-2,2.6904e-2,2.9458e-2,1.9892e-2/
 591    data bbp/-1.3139e-4,-5.5688e-5,-4.6380e-5,-8.0362e-5,-1.0115e-4,-8.8061e-5/
 592 !
 593 !--------------------------Statement Functions--------------------------
 594 !
 595    real(r8) func, u, b
 596    func(u,b) = u/sqrt(4.0 + u*(1.0 + 1.0 / b))
 597 !
 598 !------------------------------------------------------------------------
 599 !
 600    do i = 1,ncol
 601       sqti(i) = sqrt(to3co2(i))
 602 !
 603 ! h2o transmission
 604 !
 605       tt(i) = abs(to3co2(i) - 250.0)
 606       ds2c(i) = abs(s2c(i,k1) - s2c(i,k2))
 607       duptyp(i) = abs(uptype(i,k1) - uptype(i,k2))
 608    end do
 609 !
 610    do l = 1,6
 611       do i = 1,ncol
 612          psi1 = exp(abp(l)*tt(i) + bbp(l)*tt(i)*tt(i))
 613          phi1 = exp(ab(l)*tt(i) + bb(l)*tt(i)*tt(i))
 614          p1 = pnew(i)*(psi1/phi1)/sslp
 615          w1 = dw(i)*phi1
 616          tw(i,l) = exp(-g1(l)*p1*(sqrt(1.0 + g2(l)*(w1/p1)) - 1.0) - &
 617                    g3(l)*ds2c(i)-g4(l)*duptyp(i))
 618       end do
 619    end do
 620 !
 621    do i=1,ncol
 622       tw(i,1)=tw(i,1)*(0.7*aer_trn_ttl(i,k1,k2,idx_LW_0650_0800)+&! l=1: 0750--0820 cm-1
 623                        0.3*aer_trn_ttl(i,k1,k2,idx_LW_0800_1000))
 624       tw(i,2)=tw(i,2)*aer_trn_ttl(i,k1,k2,idx_LW_0800_1000) ! l=2: 0820--0880 cm-1
 625       tw(i,3)=tw(i,3)*aer_trn_ttl(i,k1,k2,idx_LW_0800_1000) ! l=3: 0880--0900 cm-1
 626       tw(i,4)=tw(i,4)*aer_trn_ttl(i,k1,k2,idx_LW_0800_1000) ! l=4: 0900--1000 cm-1
 627       tw(i,5)=tw(i,5)*aer_trn_ttl(i,k1,k2,idx_LW_1000_1200) ! l=5: 1000--1120 cm-1
 628       tw(i,6)=tw(i,6)*aer_trn_ttl(i,k1,k2,idx_LW_1000_1200) ! l=6: 1120--1170 cm-1
 629    end do                    ! end loop over lon
 630    do i = 1,ncol
 631       du1 = abs(ucfc11(i,k1) - ucfc11(i,k2))
 632       du2 = abs(ucfc12(i,k1) - ucfc12(i,k2))
 633 !
 634 ! cfc transmissions
 635 !
 636       tcfc3 = exp(-175.005*du1)
 637       tcfc4 = exp(-1202.18*du1)
 638       tcfc6 = exp(-5786.73*du2)
 639       tcfc7 = exp(-2873.51*du2)
 640       tcfc8 = exp(-2085.59*du2)
 641 !
 642 ! Absorptivity for CFC11 bands
 643 !
 644       acfc1 =  50.0*(1.0 - exp(-54.09*du1))*tw(i,1)*abplnk1(7,i,k2)
 645       acfc2 =  60.0*(1.0 - exp(-5130.03*du1))*tw(i,2)*abplnk1(8,i,k2)
 646       acfc3 =  60.0*(1.0 - tcfc3)*tw(i,4)*tcfc6*abplnk1(9,i,k2)
 647       acfc4 = 100.0*(1.0 - tcfc4)*tw(i,5)*abplnk1(10,i,k2)
 648 !
 649 ! Absorptivity for CFC12 bands
 650 !
 651       acfc5 = 45.0*(1.0 - exp(-1272.35*du2))*tw(i,3)*abplnk1(11,i,k2)
 652       acfc6 = 50.0*(1.0 - tcfc6)* tw(i,4) * abplnk1(12,i,k2)
 653       acfc7 = 80.0*(1.0 - tcfc7)* tw(i,5) * tcfc4*abplnk1(13,i,k2)
 654       acfc8 = 70.0*(1.0 - tcfc8)* tw(i,6) * abplnk1(14,i,k2)
 655 !
 656 ! Emissivity for CH4 band 1306 cm-1
 657 !
 658       tlw = exp(-1.0*sqrt(dplh2o(i)))
 659       tlw=tlw*aer_trn_ttl(i,k1,k2,idx_LW_1200_2000)
 660       duch4 = abs(uch4(i,k1) - uch4(i,k2))
 661       dbetac = abs(bch4(i,k1) - bch4(i,k2))/duch4
 662       ach4 = 6.00444*sqti(i)*log(1.0 + func(duch4,dbetac))*tlw*abplnk1(3,i,k2)
 663       tch4 = 1.0/(1.0 + 0.02*func(duch4,dbetac))
 664 !
 665 ! Absorptivity for N2O bands
 666 !
 667       du01 = abs(un2o0(i,k1) - un2o0(i,k2))
 668       du11 = abs(un2o1(i,k1) - un2o1(i,k2))
 669       dbeta01 = abs(bn2o0(i,k1) - bn2o0(i,k2))/du01
 670       dbeta11 = abs(bn2o1(i,k1) - bn2o1(i,k2))/du11
 671 !
 672 ! 1285 cm-1 band
 673 !
 674       an2o1 = 2.35558*sqti(i)*log(1.0 + func(du01,dbeta01) &
 675               + func(du11,dbeta11))*tlw*tch4*abplnk1(4,i,k2)
 676       du02 = 0.100090*du01
 677       du12 = 0.0992746*du11
 678       dbeta02 = 0.964282*dbeta01
 679 !
 680 ! 589 cm-1 band
 681 !
 682       an2o2 = 2.65581*sqti(i)*log(1.0 + func(du02,dbeta02) + &
 683               func(du12,dbeta02))*th2o(i)*tco2(i)*abplnk1(5,i,k2)
 684       du03 = 0.0333767*du01
 685       dbeta03 = 0.982143*dbeta01
 686 !
 687 ! 1168 cm-1 band
 688 !
 689       an2o3 = 2.54034*sqti(i)*log(1.0 + func(du03,dbeta03))* &
 690               tw(i,6)*tcfc8*abplnk1(6,i,k2)
 691 !
 692 ! Emissivity for 1064 cm-1 band of CO2
 693 !
 694       du11 = abs(uco211(i,k1) - uco211(i,k2))
 695       du12 = abs(uco212(i,k1) - uco212(i,k2))
 696       du13 = abs(uco213(i,k1) - uco213(i,k2))
 697       dbetc1 = 2.97558*abs(pnm(i,k1) + pnm(i,k2))/(2.0*sslp*sqti(i))
 698       dbetc2 = 2.0*dbetc1
 699       aco21 = 3.7571*sqti(i)*log(1.0 + func(du11,dbetc1) &
 700               + func(du12,dbetc2) + func(du13,dbetc2)) &
 701               *to3(i)*tw(i,5)*tcfc4*tcfc7*abplnk1(2,i,k2)
 702 !
 703 ! Emissivity for 961 cm-1 band
 704 !
 705       du21 = abs(uco221(i,k1) - uco221(i,k2))
 706       du22 = abs(uco222(i,k1) - uco222(i,k2))
 707       du23 = abs(uco223(i,k1) - uco223(i,k2))
 708       aco22 = 3.8443*sqti(i)*log(1.0 + func(du21,dbetc1) &
 709               + func(du22,dbetc1) + func(du23,dbetc2)) &
 710               *tw(i,4)*tcfc3*tcfc6*abplnk1(1,i,k2)
 711 !
 712 ! total trace gas absorptivity
 713 !
 714       abstrc(i) = acfc1 + acfc2 + acfc3 + acfc4 + acfc5 + acfc6 + &
 715                   acfc7 + acfc8 + an2o1 + an2o2 + an2o3 + ach4 + &
 716                   aco21 + aco22
 717    end do
 718 !
 719    return
 720 !
 721 end subroutine trcab
 722
 723
 724
 725 subroutine trcabn(lchnk   ,ncol    ,pcols, pverp,               &
 726                   k2      ,kn      ,ucfc11  ,ucfc12  ,un2o0   , &
 727                   un2o1   ,uch4    ,uco211  ,uco212  ,uco213  , &
 728                   uco221  ,uco222  ,uco223  ,tbar    ,bplnk   , &
 729                   winpl   ,pinpl   ,tco2    ,th2o    ,to3     , &
 730                   uptype  ,dw      ,s2c     ,up2     ,pnew    , &
 731                   abstrc  ,uinpl   , &
 732                   aer_trn_ngh)
 733 !-----------------------------------------------------------------------
 734 !
 735 ! Purpose:
 736 ! Calculate nearest layer absorptivity due to CH4, N2O, CFC11 and CFC12
 737 !
 738 ! Method:
 739 ! Equations in CCM3 description
 740 !
 741 ! Author: J. Kiehl
 742 !
 743 !-----------------------------------------------------------------------
 744 !
 745 !  use shr_kind_mod, only: r8 => shr_kind_r8
 746 !  use ppgrid
 747 !  use volcrad
 748
 749    implicit none
 750
 751 !------------------------------Arguments--------------------------------
 752 !
 753 ! Input arguments
 754 !
 755    integer, intent(in) :: lchnk                 ! chunk identifier
 756    integer, intent(in) :: ncol                  ! number of atmospheric columns
 757    integer, intent(in) :: pcols, pverp
 758    integer, intent(in) :: k2                    ! level index
 759    integer, intent(in) :: kn                    ! level index
 760 !
 761    real(r8), intent(in) :: tbar(pcols,4)        ! pressure weighted temperature
 762    real(r8), intent(in) :: ucfc11(pcols,pverp)  ! CFC11 path length
 763    real(r8), intent(in) :: ucfc12(pcols,pverp)  ! CFC12 path length
 764    real(r8), intent(in) :: un2o0(pcols,pverp)   ! N2O path length
 765    real(r8), intent(in) :: un2o1(pcols,pverp)   ! N2O path length (hot band)
 766 !
 767    real(r8), intent(in) :: uch4(pcols,pverp)    ! CH4 path length
 768    real(r8), intent(in) :: uco211(pcols,pverp)  ! CO2 9.4 micron band path length
 769    real(r8), intent(in) :: uco212(pcols,pverp)  ! CO2 9.4 micron band path length
 770    real(r8), intent(in) :: uco213(pcols,pverp)  ! CO2 9.4 micron band path length
 771    real(r8), intent(in) :: uco221(pcols,pverp)  ! CO2 10.4 micron band path length
 772 !
 773    real(r8), intent(in) :: uco222(pcols,pverp)  ! CO2 10.4 micron band path length
 774    real(r8), intent(in) :: uco223(pcols,pverp)  ! CO2 10.4 micron band path length
 775    real(r8), intent(in) :: bplnk(14,pcols,4)    ! weighted Planck fnc. for absorptivity
 776    real(r8), intent(in) :: winpl(pcols,4)       ! fractional path length
 777    real(r8), intent(in) :: pinpl(pcols,4)       ! pressure factor for subdivided layer
 778 !
 779    real(r8), intent(in) :: tco2(pcols)          ! co2 transmission
 780    real(r8), intent(in) :: th2o(pcols)          ! h2o transmission
 781    real(r8), intent(in) :: to3(pcols)           ! o3 transmission
 782    real(r8), intent(in) :: dw(pcols)            ! h2o path length
 783    real(r8), intent(in) :: pnew(pcols)          ! pressure factor
 784 !
 785    real(r8), intent(in) :: s2c(pcols,pverp)     ! h2o continuum factor
 786    real(r8), intent(in) :: uptype(pcols,pverp)  ! p-type path length
 787    real(r8), intent(in) :: up2(pcols)           ! p squared path length
 788    real(r8), intent(in) :: uinpl(pcols,4)       ! Nearest layer subdivision factor
 789    real(r8), intent(in) :: aer_trn_ngh(pcols,bnd_nbr_LW)
 790                              ! [fraction] Total transmission between
 791                              !            nearest neighbor sub-levels
 792 !
 793 !  Output Arguments
 794 !
 795    real(r8), intent(out) :: abstrc(pcols)        ! total trace gas absorptivity
 796
 797 !
 798 !--------------------------Local Variables------------------------------
 799 !
 800    integer i,l                   ! loop counters
 801 !
 802    real(r8) sqti(pcols)          ! square root of mean temp
 803    real(r8) rsqti(pcols)         ! reciprocal of sqti
 804    real(r8) du1                  ! cfc11 path length
 805    real(r8) du2                  ! cfc12 path length
 806    real(r8) acfc1                ! absorptivity of cfc11 798 cm-1 band
 807 !
 808    real(r8) acfc2                ! absorptivity of cfc11 846 cm-1 band
 809    real(r8) acfc3                ! absorptivity of cfc11 933 cm-1 band
 810    real(r8) acfc4                ! absorptivity of cfc11 1085 cm-1 band
 811    real(r8) acfc5                ! absorptivity of cfc11 889 cm-1 band
 812    real(r8) acfc6                ! absorptivity of cfc11 923 cm-1 band
 813 !
 814    real(r8) acfc7                ! absorptivity of cfc11 1102 cm-1 band
 815    real(r8) acfc8                ! absorptivity of cfc11 1161 cm-1 band
 816    real(r8) du01                 ! n2o path length
 817    real(r8) dbeta01              ! n2o pressure factors
 818    real(r8) dbeta11              !        "
 819 !
 820    real(r8)  an2o1               ! absorptivity of the 1285 cm-1 n2o band
 821    real(r8) du02                 ! n2o path length
 822    real(r8) dbeta02              ! n2o pressure factor
 823    real(r8) an2o2                ! absorptivity of the 589 cm-1 n2o band
 824    real(r8) du03                 ! n2o path length
 825 !
 826    real(r8) dbeta03              ! n2o pressure factor
 827    real(r8) an2o3                ! absorptivity of the 1168 cm-1 n2o band
 828    real(r8) duch4                ! ch4 path length
 829    real(r8) dbetac               ! ch4 pressure factor
 830    real(r8) ach4                 ! absorptivity of the 1306 cm-1 ch4 band
 831 !
 832    real(r8) du11                 ! co2 path length
 833    real(r8) du12                 !       "
 834    real(r8) du13                 !       "
 835    real(r8) dbetc1               ! co2 pressure factor
 836    real(r8) dbetc2               ! co2 pressure factor
 837 !
 838    real(r8) aco21                ! absorptivity of the 1064 cm-1 co2 band
 839    real(r8) du21                 ! co2 path length
 840    real(r8) du22                 !       "
 841    real(r8) du23                 !       "
 842    real(r8) aco22                ! absorptivity of the 961 cm-1 co2 band
 843 !
 844    real(r8) tt(pcols)            ! temp. factor for h2o overlap
 845    real(r8) psi1                 !          "
 846    real(r8) phi1                 !          "
 847    real(r8) p1                   ! factor for h2o overlap
 848    real(r8) w1                   !          "
 849 !
 850    real(r8) ds2c(pcols)          ! continuum path length
 851    real(r8) duptyp(pcols)        ! p-type path length
 852    real(r8) tw(pcols,6)          ! h2o transmission overlap
 853    real(r8) g1(6)                ! h2o overlap factor
 854    real(r8) g2(6)                !         "
 855 !
 856    real(r8) g3(6)                !         "
 857    real(r8) g4(6)                !         "
 858    real(r8) ab(6)                ! h2o temp. factor
 859    real(r8) bb(6)                !         "
 860    real(r8) abp(6)               !         "
 861 !
 862    real(r8) bbp(6)               !         "
 863    real(r8) tcfc3                ! transmission of cfc11 band
 864    real(r8) tcfc4                ! transmission of cfc11 band
 865    real(r8) tcfc6                ! transmission of cfc12 band
 866    real(r8) tcfc7                !         "
 867 !
 868    real(r8) tcfc8                !         "
 869    real(r8) tlw                  ! h2o transmission
 870    real(r8) tch4                 ! ch4 transmission
 871 !
 872 !--------------------------Data Statements------------------------------
 873 !
 874    data g1 /0.0468556,0.0397454,0.0407664,0.0304380,0.0540398,0.0321962/
 875    data g2 /14.4832,4.30242,5.23523,3.25342,0.698935,16.5599/
 876    data g3 /26.1898,18.4476,15.3633,12.1927,9.14992,8.07092/
 877    data g4 /0.0261782,0.0369516,0.0307266,0.0243854,0.0182932,0.0161418/
 878    data ab /3.0857e-2,2.3524e-2,1.7310e-2,2.6661e-2,2.8074e-2,2.2915e-2/
 879    data bb /-1.3512e-4,-6.8320e-5,-3.2609e-5,-1.0228e-5,-9.5743e-5,-1.0304e-4/
 880    data abp/2.9129e-2,2.4101e-2,1.9821e-2,2.6904e-2,2.9458e-2,1.9892e-2/
 881    data bbp/-1.3139e-4,-5.5688e-5,-4.6380e-5,-8.0362e-5,-1.0115e-4,-8.8061e-5/
 882 !
 883 !--------------------------Statement Functions--------------------------
 884 !
 885    real(r8) func, u, b
 886    func(u,b) = u/sqrt(4.0 + u*(1.0 + 1.0 / b))
 887 !
 888 !------------------------------------------------------------------
 889 !
 890    do i = 1,ncol
 891       sqti(i) = sqrt(tbar(i,kn))
 892       rsqti(i) = 1. / sqti(i)
 893 !
 894 ! h2o transmission
 895 !
 896       tt(i) = abs(tbar(i,kn) - 250.0)
 897       ds2c(i) = abs(s2c(i,k2+1) - s2c(i,k2))*uinpl(i,kn)
 898       duptyp(i) = abs(uptype(i,k2+1) - uptype(i,k2))*uinpl(i,kn)
 899    end do
 900 !
 901    do l = 1,6
 902       do i = 1,ncol
 903          psi1 = exp(abp(l)*tt(i)+bbp(l)*tt(i)*tt(i))
 904          phi1 = exp(ab(l)*tt(i)+bb(l)*tt(i)*tt(i))
 905          p1 = pnew(i) * (psi1/phi1) / sslp
 906          w1 = dw(i) * winpl(i,kn) * phi1
 907          tw(i,l) = exp(- g1(l)*p1*(sqrt(1.0+g2(l)*(w1/p1))-1.0) &
 908                    - g3(l)*ds2c(i)-g4(l)*duptyp(i))
 909       end do
 910    end do
 911 !
 912    do i=1,ncol
 913       tw(i,1)=tw(i,1)*(0.7*aer_trn_ngh(i,idx_LW_0650_0800)+&! l=1: 0750--0820 cm-1
 914                        0.3*aer_trn_ngh(i,idx_LW_0800_1000))
 915       tw(i,2)=tw(i,2)*aer_trn_ngh(i,idx_LW_0800_1000) ! l=2: 0820--0880 cm-1
 916       tw(i,3)=tw(i,3)*aer_trn_ngh(i,idx_LW_0800_1000) ! l=3: 0880--0900 cm-1
 917       tw(i,4)=tw(i,4)*aer_trn_ngh(i,idx_LW_0800_1000) ! l=4: 0900--1000 cm-1
 918       tw(i,5)=tw(i,5)*aer_trn_ngh(i,idx_LW_1000_1200) ! l=5: 1000--1120 cm-1
 919       tw(i,6)=tw(i,6)*aer_trn_ngh(i,idx_LW_1000_1200) ! l=6: 1120--1170 cm-1
 920    end do                    ! end loop over lon
 921
 922    do i = 1,ncol
 923 !
 924       du1 = abs(ucfc11(i,k2+1) - ucfc11(i,k2)) * winpl(i,kn)
 925       du2 = abs(ucfc12(i,k2+1) - ucfc12(i,k2)) * winpl(i,kn)
 926 !
 927 ! cfc transmissions
 928 !
 929       tcfc3 = exp(-175.005*du1)
 930       tcfc4 = exp(-1202.18*du1)
 931       tcfc6 = exp(-5786.73*du2)
 932       tcfc7 = exp(-2873.51*du2)
 933       tcfc8 = exp(-2085.59*du2)
 934 !
 935 ! Absorptivity for CFC11 bands
 936 !
 937       acfc1 = 50.0*(1.0 - exp(-54.09*du1)) * tw(i,1)*bplnk(7,i,kn)
 938       acfc2 = 60.0*(1.0 - exp(-5130.03*du1))*tw(i,2)*bplnk(8,i,kn)
 939       acfc3 = 60.0*(1.0 - tcfc3)*tw(i,4)*tcfc6 * bplnk(9,i,kn)
 940       acfc4 = 100.0*(1.0 - tcfc4)* tw(i,5) * bplnk(10,i,kn)
 941 !
 942 ! Absorptivity for CFC12 bands
 943 !
 944       acfc5 = 45.0*(1.0 - exp(-1272.35*du2))*tw(i,3)*bplnk(11,i,kn)
 945       acfc6 = 50.0*(1.0 - tcfc6)*tw(i,4)*bplnk(12,i,kn)
 946       acfc7 = 80.0*(1.0 - tcfc7)* tw(i,5)*tcfc4 *bplnk(13,i,kn)
 947       acfc8 = 70.0*(1.0 - tcfc8)*tw(i,6)*bplnk(14,i,kn)
 948 !
 949 ! Absorptivity for CH4 band 1306 cm-1
 950 !
 951       tlw = exp(-1.0*sqrt(up2(i)))
 952       tlw=tlw*aer_trn_ngh(i,idx_LW_1200_2000)
 953       duch4 = abs(uch4(i,k2+1) - uch4(i,k2)) * winpl(i,kn)
 954       dbetac = 2.94449 * pinpl(i,kn) * rsqti(i) / sslp
 955       ach4 = 6.00444*sqti(i)*log(1.0 + func(duch4,dbetac)) * tlw * bplnk(3,i,kn)
 956       tch4 = 1.0/(1.0 + 0.02*func(duch4,dbetac))
 957 !
 958 ! Absorptivity for N2O bands
 959 !
 960       du01 = abs(un2o0(i,k2+1) - un2o0(i,k2)) * winpl(i,kn)
 961       du11 = abs(un2o1(i,k2+1) - un2o1(i,k2)) * winpl(i,kn)
 962       dbeta01 = 19.399 *  pinpl(i,kn) * rsqti(i) / sslp
 963       dbeta11 = dbeta01
 964 !
 965 ! 1285 cm-1 band
 966 !
 967       an2o1 = 2.35558*sqti(i)*log(1.0 + func(du01,dbeta01) &
 968               + func(du11,dbeta11)) * tlw * tch4 * bplnk(4,i,kn)
 969       du02 = 0.100090*du01
 970       du12 = 0.0992746*du11
 971       dbeta02 = 0.964282*dbeta01
 972 !
 973 ! 589 cm-1 band
 974 !
 975       an2o2 = 2.65581*sqti(i)*log(1.0 + func(du02,dbeta02) &
 976               +  func(du12,dbeta02)) * tco2(i) * th2o(i) * bplnk(5,i,kn)
 977       du03 = 0.0333767*du01
 978       dbeta03 = 0.982143*dbeta01
 979 !
 980 ! 1168 cm-1 band
 981 !
 982       an2o3 = 2.54034*sqti(i)*log(1.0 + func(du03,dbeta03)) * &
 983               tw(i,6) * tcfc8 * bplnk(6,i,kn)
 984 !
 985 ! Absorptivity for 1064 cm-1 band of CO2
 986 !
 987       du11 = abs(uco211(i,k2+1) - uco211(i,k2)) * winpl(i,kn)
 988       du12 = abs(uco212(i,k2+1) - uco212(i,k2)) * winpl(i,kn)
 989       du13 = abs(uco213(i,k2+1) - uco213(i,k2)) * winpl(i,kn)
 990       dbetc1 = 2.97558 * pinpl(i,kn) * rsqti(i) / sslp
 991       dbetc2 = 2.0 * dbetc1
 992       aco21 = 3.7571*sqti(i)*log(1.0 + func(du11,dbetc1) &
 993               + func(du12,dbetc2) + func(du13,dbetc2)) &
 994               * to3(i) * tw(i,5) * tcfc4 * tcfc7 * bplnk(2,i,kn)
 995 !
 996 ! Absorptivity for 961 cm-1 band of co2
 997 !
 998       du21 = abs(uco221(i,k2+1) - uco221(i,k2)) * winpl(i,kn)
 999       du22 = abs(uco222(i,k2+1) - uco222(i,k2)) * winpl(i,kn)
1000       du23 = abs(uco223(i,k2+1) - uco223(i,k2)) * winpl(i,kn)
1001       aco22 = 3.8443*sqti(i)*log(1.0 + func(du21,dbetc1) &
1002               + func(du22,dbetc1) + func(du23,dbetc2)) &
1003               * tw(i,4) * tcfc3 * tcfc6 * bplnk(1,i,kn)
1004 !
1005 ! total trace gas absorptivity
1006 !
1007       abstrc(i) = acfc1 + acfc2 + acfc3 + acfc4 + acfc5 + acfc6 + &
1008                   acfc7 + acfc8 + an2o1 + an2o2 + an2o3 + ach4 + &
1009                   aco21 + aco22
1010    end do
1011 !
1012    return
1013 !
1014 end subroutine trcabn
1015
1016
1017
1018 subroutine trcems(lchnk   ,ncol    ,pcols, pverp,               &
1019                   k       ,co2t    ,pnm     ,ucfc11  ,ucfc12  , &
1020                   un2o0   ,un2o1   ,bn2o0   ,bn2o1   ,uch4    , &
1021                   bch4    ,uco211  ,uco212  ,uco213  ,uco221  , &
1022                   uco222  ,uco223  ,uptype  ,w       ,s2c     , &
1023                   up2     ,emplnk  ,th2o    ,tco2    ,to3     , &
1024                   emstrc  , &
1025                  aer_trn_ttl)
1026 !-----------------------------------------------------------------------
1027 !
1028 ! Purpose:
1029 !  Calculate emissivity for CH4, N2O, CFC11 and CFC12 bands.
1030 !
1031 ! Method:
1032 !  See CCM3 Description for equations.
1033 !
1034 ! Author: J. Kiehl
1035 !
1036 !-----------------------------------------------------------------------
1037 !  use shr_kind_mod, only: r8 => shr_kind_r8
1038 !  use ppgrid
1039 !  use volcrad
1040
1041    implicit none
1042
1043 !
1044 !------------------------------Arguments--------------------------------
1045 !
1046 ! Input arguments
1047 !
1048    integer, intent(in) :: lchnk                 ! chunk identifier
1049    integer, intent(in) :: ncol                  ! number of atmospheric columns
1050    integer, intent(in) :: pcols, pverp
1051
1052    real(r8), intent(in) :: co2t(pcols,pverp)    ! pressure weighted temperature
1053    real(r8), intent(in) :: pnm(pcols,pverp)     ! interface pressure
1054    real(r8), intent(in) :: ucfc11(pcols,pverp)  ! CFC11 path length
1055    real(r8), intent(in) :: ucfc12(pcols,pverp)  ! CFC12 path length
1056    real(r8), intent(in) :: un2o0(pcols,pverp)   ! N2O path length
1057 !
1058    real(r8), intent(in) :: un2o1(pcols,pverp)   ! N2O path length (hot band)
1059    real(r8), intent(in) :: uch4(pcols,pverp)    ! CH4 path length
1060    real(r8), intent(in) :: uco211(pcols,pverp)  ! CO2 9.4 micron band path length
1061    real(r8), intent(in) :: uco212(pcols,pverp)  ! CO2 9.4 micron band path length
1062    real(r8), intent(in) :: uco213(pcols,pverp)  ! CO2 9.4 micron band path length
1063 !
1064    real(r8), intent(in) :: uco221(pcols,pverp)  ! CO2 10.4 micron band path length
1065    real(r8), intent(in) :: uco222(pcols,pverp)  ! CO2 10.4 micron band path length
1066    real(r8), intent(in) :: uco223(pcols,pverp)  ! CO2 10.4 micron band path length
1067    real(r8), intent(in) :: uptype(pcols,pverp)  ! continuum path length
1068    real(r8), intent(in) :: bn2o0(pcols,pverp)   ! pressure factor for n2o
1069 !
1070    real(r8), intent(in) :: bn2o1(pcols,pverp)   ! pressure factor for n2o
1071    real(r8), intent(in) :: bch4(pcols,pverp)    ! pressure factor for ch4
1072    real(r8), intent(in) :: emplnk(14,pcols)     ! emissivity Planck factor
1073    real(r8), intent(in) :: th2o(pcols)          ! water vapor overlap factor
1074    real(r8), intent(in) :: tco2(pcols)          ! co2 overlap factor
1075 !
1076    real(r8), intent(in) :: to3(pcols)           ! o3 overlap factor
1077    real(r8), intent(in) :: s2c(pcols,pverp)     ! h2o continuum path length
1078    real(r8), intent(in) :: w(pcols,pverp)       ! h2o path length
1079    real(r8), intent(in) :: up2(pcols)           ! pressure squared h2o path length
1080 !
1081    integer, intent(in) :: k                 ! level index
1082
1083    real(r8), intent(in) :: aer_trn_ttl(pcols,pverp,pverp,bnd_nbr_LW) ! aer trn.
1084
1085 !
1086 !  Output Arguments
1087 !
1088    real(r8), intent(out) :: emstrc(pcols,pverp)  ! total trace gas emissivity
1089
1090 !
1091 !--------------------------Local Variables------------------------------
1092 !
1093    integer i,l               ! loop counters
1094 !
1095    real(r8) sqti(pcols)          ! square root of mean temp
1096    real(r8) ecfc1                ! emissivity of cfc11 798 cm-1 band
1097    real(r8) ecfc2                !     "      "    "   846 cm-1 band
1098    real(r8) ecfc3                !     "      "    "   933 cm-1 band
1099    real(r8) ecfc4                !     "      "    "   1085 cm-1 band
1100 !
1101    real(r8) ecfc5                !     "      "  cfc12 889 cm-1 band
1102    real(r8) ecfc6                !     "      "    "   923 cm-1 band
1103    real(r8) ecfc7                !     "      "    "   1102 cm-1 band
1104    real(r8) ecfc8                !     "      "    "   1161 cm-1 band
1105    real(r8) u01                  ! n2o path length
1106 !
1107    real(r8) u11                  ! n2o path length
1108    real(r8) beta01               ! n2o pressure factor
1109    real(r8) beta11               ! n2o pressure factor
1110    real(r8) en2o1                ! emissivity of the 1285 cm-1 N2O band
1111    real(r8) u02                  ! n2o path length
1112 !
1113    real(r8) u12                  ! n2o path length
1114    real(r8) beta02               ! n2o pressure factor
1115    real(r8) en2o2                ! emissivity of the 589 cm-1 N2O band
1116    real(r8) u03                  ! n2o path length
1117    real(r8) beta03               ! n2o pressure factor
1118 !
1119    real(r8) en2o3                ! emissivity of the 1168 cm-1 N2O band
1120    real(r8) betac                ! ch4 pressure factor
1121    real(r8) ech4                 ! emissivity of 1306 cm-1 CH4 band
1122    real(r8) betac1               ! co2 pressure factor
1123    real(r8) betac2               ! co2 pressure factor
1124 !
1125    real(r8) eco21                ! emissivity of 1064 cm-1 CO2 band
1126    real(r8) eco22                ! emissivity of 961 cm-1 CO2 band
1127    real(r8) tt(pcols)            ! temp. factor for h2o overlap factor
1128    real(r8) psi1                 ! narrow band h2o temp. factor
1129    real(r8) phi1                 !             "
1130 !
1131    real(r8) p1                   ! h2o line overlap factor
1132    real(r8) w1                   !          "
1133    real(r8) tw(pcols,6)          ! h2o transmission overlap
1134    real(r8) g1(6)                ! h2o overlap factor
1135    real(r8) g2(6)                !          "
1136 !
1137    real(r8) g3(6)                !          "
1138    real(r8) g4(6)                !          "
1139    real(r8) ab(6)                !          "
1140    real(r8) bb(6)                !          "
1141    real(r8) abp(6)               !          "
1142 !
1143    real(r8) bbp(6)               !          "
1144    real(r8) tcfc3                ! transmission for cfc11 band
1145    real(r8) tcfc4                !          "
1146    real(r8) tcfc6                ! transmission for cfc12 band
1147    real(r8) tcfc7                !          "
1148 !
1149    real(r8) tcfc8                !          "
1150    real(r8) tlw                  ! h2o overlap factor
1151    real(r8) tch4                 ! ch4 overlap factor
1152 !
1153 !--------------------------Data Statements------------------------------
1154 !
1155    data g1 /0.0468556,0.0397454,0.0407664,0.0304380,0.0540398,0.0321962/
1156    data g2 /14.4832,4.30242,5.23523,3.25342,0.698935,16.5599/
1157    data g3 /26.1898,18.4476,15.3633,12.1927,9.14992,8.07092/
1158    data g4 /0.0261782,0.0369516,0.0307266,0.0243854,0.0182932,0.0161418/
1159    data ab /3.0857e-2,2.3524e-2,1.7310e-2,2.6661e-2,2.8074e-2,2.2915e-2/
1160    data bb /-1.3512e-4,-6.8320e-5,-3.2609e-5,-1.0228e-5,-9.5743e-5,-1.0304e-4/
1161    data abp/2.9129e-2,2.4101e-2,1.9821e-2,2.6904e-2,2.9458e-2,1.9892e-2/
1162    data bbp/-1.3139e-4,-5.5688e-5,-4.6380e-5,-8.0362e-5,-1.0115e-4,-8.8061e-5/
1163 !
1164 !--------------------------Statement Functions--------------------------
1165 !
1166    real(r8) func, u, b
1167    func(u,b) = u/sqrt(4.0 + u*(1.0 + 1.0 / b))
1168 !
1169 !-----------------------------------------------------------------------
1170 !
1171    do i = 1,ncol
1172       sqti(i) = sqrt(co2t(i,k))
1173 !
1174 ! Transmission for h2o
1175 !
1176       tt(i) = abs(co2t(i,k) - 250.0)
1177    end do
1178 !
1179    do l = 1,6
1180       do i = 1,ncol
1181          psi1 = exp(abp(l)*tt(i)+bbp(l)*tt(i)*tt(i))
1182          phi1 = exp(ab(l)*tt(i)+bb(l)*tt(i)*tt(i))
1183          p1 = pnm(i,k) * (psi1/phi1) / sslp
1184          w1 = w(i,k) * phi1
1185          tw(i,l) = exp(- g1(l)*p1*(sqrt(1.0+g2(l)*(w1/p1))-1.0) &
1186                    - g3(l)*s2c(i,k)-g4(l)*uptype(i,k))
1187       end do
1188    end do
1189
1190 !     Overlap H2O tranmission with STRAER continuum in 6 trace gas
1191 !                 subbands
1192
1193       do i=1,ncol
1194          tw(i,1)=tw(i,1)*(0.7*aer_trn_ttl(i,k,1,idx_LW_0650_0800)+&! l=1: 0750--0820 cm-1
1195                           0.3*aer_trn_ttl(i,k,1,idx_LW_0800_1000))
1196          tw(i,2)=tw(i,2)*aer_trn_ttl(i,k,1,idx_LW_0800_1000) ! l=2: 0820--0880 cm-1
1197          tw(i,3)=tw(i,3)*aer_trn_ttl(i,k,1,idx_LW_0800_1000) ! l=3: 0880--0900 cm-1
1198          tw(i,4)=tw(i,4)*aer_trn_ttl(i,k,1,idx_LW_0800_1000) ! l=4: 0900--1000 cm-1
1199          tw(i,5)=tw(i,5)*aer_trn_ttl(i,k,1,idx_LW_1000_1200) ! l=5: 1000--1120 cm-1
1200          tw(i,6)=tw(i,6)*aer_trn_ttl(i,k,1,idx_LW_1000_1200) ! l=6: 1120--1170 cm-1
1201       end do                    ! end loop over lon
1202 !
1203    do i = 1,ncol
1204 !
1205 ! transmission due to cfc bands
1206 !
1207       tcfc3 = exp(-175.005*ucfc11(i,k))
1208       tcfc4 = exp(-1202.18*ucfc11(i,k))
1209       tcfc6 = exp(-5786.73*ucfc12(i,k))
1210       tcfc7 = exp(-2873.51*ucfc12(i,k))
1211       tcfc8 = exp(-2085.59*ucfc12(i,k))
1212 !
1213 ! Emissivity for CFC11 bands
1214 !
1215       ecfc1 = 50.0*(1.0 - exp(-54.09*ucfc11(i,k))) * tw(i,1) * emplnk(7,i)
1216       ecfc2 = 60.0*(1.0 - exp(-5130.03*ucfc11(i,k)))* tw(i,2) * emplnk(8,i)
1217       ecfc3 = 60.0*(1.0 - tcfc3)*tw(i,4)*tcfc6*emplnk(9,i)
1218       ecfc4 = 100.0*(1.0 - tcfc4)*tw(i,5)*emplnk(10,i)
1219 !
1220 ! Emissivity for CFC12 bands
1221 !
1222       ecfc5 = 45.0*(1.0 - exp(-1272.35*ucfc12(i,k)))*tw(i,3)*emplnk(11,i)
1223       ecfc6 = 50.0*(1.0 - tcfc6)*tw(i,4)*emplnk(12,i)
1224       ecfc7 = 80.0*(1.0 - tcfc7)*tw(i,5)* tcfc4 * emplnk(13,i)
1225       ecfc8 = 70.0*(1.0 - tcfc8)*tw(i,6) * emplnk(14,i)
1226 !
1227 ! Emissivity for CH4 band 1306 cm-1
1228 !
1229       tlw = exp(-1.0*sqrt(up2(i)))
1230
1231 !     Overlap H2O vibration rotation band with STRAER continuum
1232 !             for CH4 1306 cm-1 and N2O 1285 cm-1 bands
1233
1234             tlw=tlw*aer_trn_ttl(i,k,1,idx_LW_1200_2000)
1235       betac = bch4(i,k)/uch4(i,k)
1236       ech4 = 6.00444*sqti(i)*log(1.0 + func(uch4(i,k),betac)) *tlw * emplnk(3,i)
1237       tch4 = 1.0/(1.0 + 0.02*func(uch4(i,k),betac))
1238 !
1239 ! Emissivity for N2O bands
1240 !
1241       u01 = un2o0(i,k)
1242       u11 = un2o1(i,k)
1243       beta01 = bn2o0(i,k)/un2o0(i,k)
1244       beta11 = bn2o1(i,k)/un2o1(i,k)
1245 !
1246 ! 1285 cm-1 band
1247 !
1248       en2o1 = 2.35558*sqti(i)*log(1.0 + func(u01,beta01) + &
1249               func(u11,beta11))*tlw*tch4*emplnk(4,i)
1250       u02 = 0.100090*u01
1251       u12 = 0.0992746*u11
1252       beta02 = 0.964282*beta01
1253 !
1254 ! 589 cm-1 band
1255 !
1256       en2o2 = 2.65581*sqti(i)*log(1.0 + func(u02,beta02) + &
1257               func(u12,beta02)) * tco2(i) * th2o(i) * emplnk(5,i)
1258       u03 = 0.0333767*u01
1259       beta03 = 0.982143*beta01
1260 !
1261 ! 1168 cm-1 band
1262 !
1263       en2o3 = 2.54034*sqti(i)*log(1.0 + func(u03,beta03)) * &
1264               tw(i,6) * tcfc8 * emplnk(6,i)
1265 !
1266 ! Emissivity for 1064 cm-1 band of CO2
1267 !
1268       betac1 = 2.97558*pnm(i,k) / (sslp*sqti(i))
1269       betac2 = 2.0 * betac1
1270       eco21 = 3.7571*sqti(i)*log(1.0 + func(uco211(i,k),betac1) &
1271               + func(uco212(i,k),betac2) + func(uco213(i,k),betac2)) &
1272               * to3(i) * tw(i,5) * tcfc4 * tcfc7 * emplnk(2,i)
1273 !
1274 ! Emissivity for 961 cm-1 band
1275 !
1276       eco22 = 3.8443*sqti(i)*log(1.0 + func(uco221(i,k),betac1) &
1277               + func(uco222(i,k),betac1) + func(uco223(i,k),betac2)) &
1278               * tw(i,4) * tcfc3 * tcfc6 * emplnk(1,i)
1279 !
1280 ! total trace gas emissivity
1281 !
1282       emstrc(i,k) = ecfc1 + ecfc2 + ecfc3 + ecfc4 + ecfc5 +ecfc6 + &
1283                     ecfc7 + ecfc8 + en2o1 + en2o2 + en2o3 + ech4 + &
1284                     eco21 + eco22
1285    end do
1286 !
1287    return
1288 !
1289 end subroutine trcems
1290
1291 subroutine trcmix(lchnk   ,ncol     ,pcols, pver, &
1292                   pmid    ,clat, n2o      ,ch4     ,          &
1293                   cfc11   , cfc12   )
1294 !-----------------------------------------------------------------------
1295 !
1296 ! Purpose:
1297 ! Specify zonal mean mass mixing ratios of CH4, N2O, CFC11 and
1298 ! CFC12
1299 !
1300 ! Method:
1301 ! Distributions assume constant mixing ratio in the troposphere
1302 ! and a decrease of mixing ratio in the stratosphere. Tropopause
1303 ! defined by ptrop. The scale height of the particular trace gas
1304 ! depends on latitude. This assumption produces a more realistic
1305 ! stratospheric distribution of the various trace gases.
1306 !
1307 ! Author: J. Kiehl
1308 !
1309 !-----------------------------------------------------------------------
1310 !  use shr_kind_mod, only: r8 => shr_kind_r8
1311 !  use ppgrid
1312 !  use phys_grid,    only: get_rlat_all_p
1313 !  use physconst,    only: mwdry, mwch4, mwn2o, mwf11, mwf12
1314 !  use ghg_surfvals, only: ch4vmr, n2ovmr, f11vmr, f12vmr
1315
1316    implicit none
1317
1318 !-----------------------------Arguments---------------------------------
1319 !
1320 ! Input
1321 !
1322    integer, intent(in) :: lchnk                    ! chunk identifier
1323    integer, intent(in) :: ncol                     ! number of atmospheric columns
1324    integer, intent(in) :: pcols, pver
1325
1326    real(r8), intent(in) :: pmid(pcols,pver)        ! model pressures
1327    real(r8), intent(in) :: clat(pcols)             ! latitude in radians for columns
1328 !
1329 ! Output
1330 !
1331    real(r8), intent(out) :: n2o(pcols,pver)         ! nitrous oxide mass mixing ratio
1332    real(r8), intent(out) :: ch4(pcols,pver)         ! methane mass mixing ratio
1333    real(r8), intent(out) :: cfc11(pcols,pver)       ! cfc11 mass mixing ratio
1334    real(r8), intent(out) :: cfc12(pcols,pver)       ! cfc12 mass mixing ratio
1335
1336 !
1337 !--------------------------Local Variables------------------------------
1338
1339    real(r8) :: rmwn2o       ! ratio of molecular weight n2o   to dry air
1340    real(r8) :: rmwch4       ! ratio of molecular weight ch4   to dry air
1341    real(r8) :: rmwf11       ! ratio of molecular weight cfc11 to dry air
1342    real(r8) :: rmwf12       ! ratio of molecular weight cfc12 to dry air
1343 !
1344    integer i                ! longitude loop index
1345    integer k                ! level index
1346 !
1347 !  real(r8) clat(pcols)         ! latitude in radians for columns
1348    real(r8) coslat(pcols)       ! cosine of latitude
1349    real(r8) dlat                ! latitude in degrees
1350    real(r8) ptrop               ! pressure level of tropopause
1351    real(r8) pratio              ! pressure divided by ptrop
1352 !
1353    real(r8) xn2o                ! pressure scale height for n2o
1354    real(r8) xch4                ! pressure scale height for ch4
1355    real(r8) xcfc11              ! pressure scale height for cfc11
1356    real(r8) xcfc12              ! pressure scale height for cfc12
1357 !
1358    real(r8) ch40                ! tropospheric mass mixing ratio for ch4
1359    real(r8) n2o0                ! tropospheric mass mixing ratio for n2o
1360    real(r8) cfc110              ! tropospheric mass mixing ratio for cfc11
1361    real(r8) cfc120              ! tropospheric mass mixing ratio for cfc12
1362 !
1363 !-----------------------------------------------------------------------
1364    rmwn2o = mwn2o/mwdry      ! ratio of molecular weight n2o   to dry air
1365    rmwch4 = mwch4/mwdry      ! ratio of molecular weight ch4   to dry air
1366    rmwf11 = mwf11/mwdry      ! ratio of molecular weight cfc11 to dry air
1367    rmwf12 = mwf12/mwdry      ! ratio of molecular weight cfc12 to dry air
1368 !
1369 ! get latitudes
1370 !
1371 !  call get_rlat_all_p(lchnk, ncol, clat)
1372    do i = 1, ncol
1373       coslat(i) = cos(clat(i))
1374    end do
1375 !
1376 ! set tropospheric mass mixing ratios
1377 !
1378    ch40   = rmwch4 * ch4vmr
1379    n2o0   = rmwn2o * n2ovmr
1380    cfc110 = rmwf11 * f11vmr
1381    cfc120 = rmwf12 * f12vmr
1382
1383    do i = 1, ncol
1384       coslat(i) = cos(clat(i))
1385    end do
1386 !
1387    do k = 1,pver
1388       do i = 1,ncol
1389 !
1390 !        set stratospheric scale height factor for gases
1391          dlat = abs(57.2958 * clat(i))
1392          if(dlat.le.45.0) then
1393             xn2o = 0.3478 + 0.00116 * dlat
1394             xch4 = 0.2353
1395             xcfc11 = 0.7273 + 0.00606 * dlat
1396             xcfc12 = 0.4000 + 0.00222 * dlat
1397          else
1398             xn2o = 0.4000 + 0.013333 * (dlat - 45)
1399             xch4 = 0.2353 + 0.0225489 * (dlat - 45)
1400             xcfc11 = 1.00 + 0.013333 * (dlat - 45)
1401             xcfc12 = 0.50 + 0.024444 * (dlat - 45)
1402          end if
1403 !
1404 !        pressure of tropopause
1405          ptrop = 250.0e2 - 150.0e2*coslat(i)**2.0
1406 !
1407 !        determine output mass mixing ratios
1408          if (pmid(i,k) >= ptrop) then
1409             ch4(i,k) = ch40
1410             n2o(i,k) = n2o0
1411             cfc11(i,k) = cfc110
1412             cfc12(i,k) = cfc120
1413          else
1414             pratio = pmid(i,k)/ptrop
1415             ch4(i,k) = ch40 * (pratio)**xch4
1416             n2o(i,k) = n2o0 * (pratio)**xn2o
1417             cfc11(i,k) = cfc110 * (pratio)**xcfc11
1418             cfc12(i,k) = cfc120 * (pratio)**xcfc12
1419          end if
1420       end do
1421    end do
1422 !
1423    return
1424 !
1425 end subroutine trcmix
1426
1427 subroutine trcmix_clwrf(lchnk   ,ncol     ,pcols, pver,       &
1428                   pmid    ,clat, n2ovmr, ch4vmr, f11vmr,      &
1429                   f12vmr, n2o      ,ch4     ,                 &
1430                   cfc11   , cfc12   )
1431 !-----------------------------------------------------------------------
1432 !
1433 ! Purpose:
1434 ! Specify zonal mean mass mixing ratios of CH4, N2O, CFC11 and
1435 ! CFC12
1436 !
1437 ! Method:
1438 ! Distributions assume constant mixing ratio in the troposphere
1439 ! and a decrease of mixing ratio in the stratosphere. Tropopause
1440 ! defined by ptrop. The scale height of the particular trace gas
1441 ! depends on latitude. This assumption produces a more realistic
1442 ! stratospheric distribution of the various trace gases.
1443 !
1444 ! Author: J. Kiehl
1445 !
1446 !-----------------------------------------------------------------------
1447 !  use shr_kind_mod, only: r8 => shr_kind_r8
1448 !  use ppgrid
1449 !  use phys_grid,    only: get_rlat_all_p
1450 !  use physconst,    only: mwdry, mwch4, mwn2o, mwf11, mwf12
1451
1452    implicit none
1453
1454 !-----------------------------Arguments---------------------------------
1455 !
1456 ! Input
1457 !
1458    integer, intent(in) :: lchnk                    ! chunk identifier
1459    integer, intent(in) :: ncol                     ! number of atmospheric columns
1460    integer, intent(in) :: pcols, pver
1461
1462    real(r8), intent(in) :: pmid(pcols,pver)        ! model pressures
1463    real(r8), intent(in) :: clat(pcols)             ! latitude in radians for columns
1464    real(r8), intent(in) :: n2ovmr, ch4vmr, f11vmr, f12vmr
1465 !
1466 ! Output
1467 !
1468    real(r8), intent(out) :: n2o(pcols,pver)         ! nitrous oxide mass mixing ratio
1469    real(r8), intent(out) :: ch4(pcols,pver)         ! methane mass mixing ratio
1470    real(r8), intent(out) :: cfc11(pcols,pver)       ! cfc11 mass mixing ratio
1471    real(r8), intent(out) :: cfc12(pcols,pver)       ! cfc12 mass mixing ratio
1472
1473 !
1474 !--------------------------Local Variables------------------------------
1475
1476    real(r8) :: rmwn2o       ! ratio of molecular weight n2o   to dry air
1477    real(r8) :: rmwch4       ! ratio of molecular weight ch4   to dry air
1478    real(r8) :: rmwf11       ! ratio of molecular weight cfc11 to dry air
1479    real(r8) :: rmwf12       ! ratio of molecular weight cfc12 to dry air
1480 !
1481    integer i                ! longitude loop index
1482    integer k                ! level index
1483 !
1484 !  real(r8) clat(pcols)         ! latitude in radians for columns
1485    real(r8) coslat(pcols)       ! cosine of latitude
1486    real(r8) dlat                ! latitude in degrees
1487    real(r8) ptrop               ! pressure level of tropopause
1488    real(r8) pratio              ! pressure divided by ptrop
1489 !
1490    real(r8) xn2o                ! pressure scale height for n2o
1491    real(r8) xch4                ! pressure scale height for ch4
1492    real(r8) xcfc11              ! pressure scale height for cfc11
1493    real(r8) xcfc12              ! pressure scale height for cfc12
1494 !
1495    real(r8) ch40                ! tropospheric mass mixing ratio for ch4
1496    real(r8) n2o0                ! tropospheric mass mixing ratio for n2o
1497    real(r8) cfc110              ! tropospheric mass mixing ratio for cfc11
1498    real(r8) cfc120              ! tropospheric mass mixing ratio for cfc12
1499   CHARACTER (LEN=256)           :: message
1500 !
1501 !-----------------------------------------------------------------------
1502    rmwn2o = mwn2o/mwdry      ! ratio of molecular weight n2o   to dry air
1503    rmwch4 = mwch4/mwdry      ! ratio of molecular weight ch4   to dry air
1504    rmwf11 = mwf11/mwdry      ! ratio of molecular weight cfc11 to dry air
1505    rmwf12 = mwf12/mwdry      ! ratio of molecular weight cfc12 to dry air
1506 !
1507 ! get latitudes
1508 !
1509 !  call get_rlat_all_p(lchnk, ncol, clat)
1510    do i = 1, ncol
1511       coslat(i) = cos(clat(i))
1512    end do
1513 !
1514 ! set tropospheric mass mixing ratios
1515 !
1516    ch40   = rmwch4 * ch4vmr
1517    n2o0   = rmwn2o * n2ovmr
1518    cfc110 = rmwf11 * f11vmr
1519    cfc120 = rmwf12 * f12vmr
1520
1521    do i = 1, ncol
1522       coslat(i) = cos(clat(i))
1523    end do
1524 !
1525    do k = 1,pver
1526       do i = 1,ncol
1527 !
1528 !        set stratospheric scale height factor for gases
1529          dlat = abs(57.2958 * clat(i))
1530          if(dlat.le.45.0) then
1531             xn2o = 0.3478 + 0.00116 * dlat
1532             xch4 = 0.2353
1533             xcfc11 = 0.7273 + 0.00606 * dlat
1534             xcfc12 = 0.4000 + 0.00222 * dlat
1535          else
1536             xn2o = 0.4000 + 0.013333 * (dlat - 45)
1537             xch4 = 0.2353 + 0.0225489 * (dlat - 45)
1538             xcfc11 = 1.00 + 0.013333 * (dlat - 45)
1539             xcfc12 = 0.50 + 0.024444 * (dlat - 45)
1540          end if
1541 !
1542 !        pressure of tropopause
1543          ptrop = 250.0e2 - 150.0e2*coslat(i)**2.0
1544 !
1545 !        determine output mass mixing ratios
1546          if (pmid(i,k) >= ptrop) then
1547             ch4(i,k) = ch40
1548             n2o(i,k) = n2o0
1549             cfc11(i,k) = cfc110
1550             cfc12(i,k) = cfc120
1551          else
1552             pratio = pmid(i,k)/ptrop
1553             ch4(i,k) = ch40 * (pratio)**xch4
1554             n2o(i,k) = n2o0 * (pratio)**xn2o
1555             cfc11(i,k) = cfc110 * (pratio)**xcfc11
1556             cfc12(i,k) = cfc120 * (pratio)**xcfc12
1557          end if
1558       end do
1559    end do
1560 !
1561    WRITE(message,*)'CLWRF-trcmix. ch4vmr:',ch4vmr,' n2ovmr:',n2ovmr,' cfc11vmr:',f11vmr,' cfc12vmr:',f12vmr
1562    CALL wrf_debug(0, message)
1563    WRITE(message,*)'CLWRF-trcmix. for i= ncol/2', ncol/2,' k=pver/2',pver/2,'__________'
1564    CALL wrf_debug(0, message)
1565
1566    ptrop = 250.0e2 - 150.0e2*coslat(ncol/2)**2.0
1567    if (any(pmid >= ptrop)) then
1568       WRITE(message,*)'CLWRF-trcmix. pmid: ',pmid(ncol/2,pver/2),' ch4:',ch4(ncol/2,pver/2),' n2o:',n2o(ncol/2,pver/2),                 &
1569            ' cfc11:',cfc11(ncol/2,pver/2),' cfc12:',cfc12(ncol/2,pver/2)
1570    else
1571       WRITE(message,*)'CLWRF-trcmix. pratio: ',pmid(ncol/2,pver/2)/ptrop,' ch4:',ch4(ncol/2,pver/2),' n2o:',n2o(ncol/2,pver/2),                 &
1572            ' cfc11:',cfc11(ncol/2,pver/2),' cfc12:',cfc12(ncol/2,pver/2)
1573    endif
1574    CALL wrf_debug(0, message)
1575
1576    return
1577 !
1578 end subroutine trcmix_clwrf
1579
1580 subroutine trcplk(lchnk   ,ncol    ,pcols, pver, pverp,         &
1581                   tint    ,tlayr   ,tplnke  ,emplnk  ,abplnk1 , &
1582                   abplnk2 )
1583 !-----------------------------------------------------------------------
1584 !
1585 ! Purpose:
1586 !   Calculate Planck factors for absorptivity and emissivity of
1587 !   CH4, N2O, CFC11 and CFC12
1588 !
1589 ! Method:
1590 !   Planck function and derivative evaluated at the band center.
1591 !
1592 ! Author: J. Kiehl
1593 !
1594 !-----------------------------------------------------------------------
1595 !  use shr_kind_mod, only: r8 => shr_kind_r8
1596 !  use ppgrid
1597
1598    implicit none
1599 !------------------------------Arguments--------------------------------
1600 !
1601 ! Input arguments
1602 !
1603    integer, intent(in) :: lchnk                ! chunk identifier
1604    integer, intent(in) :: ncol                 ! number of atmospheric columns
1605    integer, intent(in) :: pcols, pver, pverp
1606
1607    real(r8), intent(in) :: tint(pcols,pverp)   ! interface temperatures
1608    real(r8), intent(in) :: tlayr(pcols,pverp)  ! k-1 level temperatures
1609    real(r8), intent(in) :: tplnke(pcols)       ! Top Layer temperature
1610 !
1611 ! output arguments
1612 !
1613    real(r8), intent(out) :: emplnk(14,pcols)         ! emissivity Planck factor
1614    real(r8), intent(out) :: abplnk1(14,pcols,pverp)  ! non-nearest layer Plack factor
1615    real(r8), intent(out) :: abplnk2(14,pcols,pverp)  ! nearest layer factor
1616
1617 !
1618 !--------------------------Local Variables------------------------------
1619 !
1620    integer wvl                   ! wavelength index
1621    integer i,k                   ! loop counters
1622 !
1623    real(r8) f1(14)                   ! Planck function factor
1624    real(r8) f2(14)                   !        "
1625    real(r8) f3(14)                   !        "
1626 !
1627 !--------------------------Data Statements------------------------------
1628 !
1629    data f1 /5.85713e8,7.94950e8,1.47009e9,1.40031e9,1.34853e8, &
1630             1.05158e9,3.35370e8,3.99601e8,5.35994e8,8.42955e8, &
1631             4.63682e8,5.18944e8,8.83202e8,1.03279e9/
1632    data f2 /2.02493e11,3.04286e11,6.90698e11,6.47333e11, &
1633             2.85744e10,4.41862e11,9.62780e10,1.21618e11, &
1634             1.79905e11,3.29029e11,1.48294e11,1.72315e11, &
1635             3.50140e11,4.31364e11/
1636    data f3 /1383.0,1531.0,1879.0,1849.0,848.0,1681.0, &
1637             1148.0,1217.0,1343.0,1561.0,1279.0,1328.0, &
1638             1586.0,1671.0/
1639 !
1640 !-----------------------------------------------------------------------
1641 !
1642 ! Calculate emissivity Planck factor
1643 !
1644    do wvl = 1,14
1645       do i = 1,ncol
1646          emplnk(wvl,i) = f1(wvl)/(tplnke(i)**4.0*(exp(f3(wvl)/tplnke(i))-1.0))
1647       end do
1648    end do
1649 !
1650 ! Calculate absorptivity Planck factor for tint and tlayr temperatures
1651 !
1652    do wvl = 1,14
1653       do k = ntoplw, pverp
1654          do i = 1, ncol
1655 !
1656 ! non-nearlest layer function
1657 !
1658             abplnk1(wvl,i,k) = (f2(wvl)*exp(f3(wvl)/tint(i,k)))  &
1659                                /(tint(i,k)**5.0*(exp(f3(wvl)/tint(i,k))-1.0)**2.0)
1660 !
1661 ! nearest layer function
1662 !
1663             abplnk2(wvl,i,k) = (f2(wvl)*exp(f3(wvl)/tlayr(i,k))) &
1664                                /(tlayr(i,k)**5.0*(exp(f3(wvl)/tlayr(i,k))-1.0)**2.0)
1665          end do
1666       end do
1667    end do
1668 !
1669    return
1670 end subroutine trcplk
1671
1672 subroutine trcpth(lchnk   ,ncol    ,pcols, pver, pverp,         &
1673                   tnm     ,pnm     ,cfc11   ,cfc12   ,n2o     , &
1674                   ch4     ,qnm     ,ucfc11  ,ucfc12  ,un2o0   , &
1675                   un2o1   ,uch4    ,uco211  ,uco212  ,uco213  , &
1676                   uco221  ,uco222  ,uco223  ,bn2o0   ,bn2o1   , &
1677                   bch4    ,uptype  )
1678 !-----------------------------------------------------------------------
1679 !
1680 ! Purpose:
1681 ! Calculate path lengths and pressure factors for CH4, N2O, CFC11
1682 ! and CFC12.
1683 !
1684 ! Method:
1685 ! See CCM3 description for details
1686 !
1687 ! Author: J. Kiehl
1688 !
1689 !-----------------------------------------------------------------------
1690 !  use shr_kind_mod, only: r8 => shr_kind_r8
1691 !  use ppgrid
1692 !  use ghg_surfvals, only: co2mmr
1693
1694    implicit none
1695
1696 !------------------------------Arguments--------------------------------
1697 !
1698 ! Input arguments
1699 !
1700    integer, intent(in) :: lchnk                 ! chunk identifier
1701    integer, intent(in) :: ncol                  ! number of atmospheric columns
1702    integer, intent(in) :: pcols, pver, pverp
1703
1704    real(r8), intent(in) :: tnm(pcols,pver)      ! Model level temperatures
1705    real(r8), intent(in) :: pnm(pcols,pverp)     ! Pres. at model interfaces (dynes/cm2)
1706    real(r8), intent(in) :: qnm(pcols,pver)      ! h2o specific humidity
1707    real(r8), intent(in) :: cfc11(pcols,pver)    ! CFC11 mass mixing ratio
1708 !
1709    real(r8), intent(in) :: cfc12(pcols,pver)    ! CFC12 mass mixing ratio
1710    real(r8), intent(in) :: n2o(pcols,pver)      ! N2O mass mixing ratio
1711    real(r8), intent(in) :: ch4(pcols,pver)      ! CH4 mass mixing ratio
1712
1713 !
1714 ! Output arguments
1715 !
1716    real(r8), intent(out) :: ucfc11(pcols,pverp)  ! CFC11 path length
1717    real(r8), intent(out) :: ucfc12(pcols,pverp)  ! CFC12 path length
1718    real(r8), intent(out) :: un2o0(pcols,pverp)   ! N2O path length
1719    real(r8), intent(out) :: un2o1(pcols,pverp)   ! N2O path length (hot band)
1720    real(r8), intent(out) :: uch4(pcols,pverp)    ! CH4 path length
1721 !
1722    real(r8), intent(out) :: uco211(pcols,pverp)  ! CO2 9.4 micron band path length
1723    real(r8), intent(out) :: uco212(pcols,pverp)  ! CO2 9.4 micron band path length
1724    real(r8), intent(out) :: uco213(pcols,pverp)  ! CO2 9.4 micron band path length
1725    real(r8), intent(out) :: uco221(pcols,pverp)  ! CO2 10.4 micron band path length
1726    real(r8), intent(out) :: uco222(pcols,pverp)  ! CO2 10.4 micron band path length
1727 !
1728    real(r8), intent(out) :: uco223(pcols,pverp)  ! CO2 10.4 micron band path length
1729    real(r8), intent(out) :: bn2o0(pcols,pverp)   ! pressure factor for n2o
1730    real(r8), intent(out) :: bn2o1(pcols,pverp)   ! pressure factor for n2o
1731    real(r8), intent(out) :: bch4(pcols,pverp)    ! pressure factor for ch4
1732    real(r8), intent(out) :: uptype(pcols,pverp)  ! p-type continuum path length
1733
1734 !
1735 !---------------------------Local variables-----------------------------
1736 !
1737    integer   i               ! Longitude index
1738    integer   k               ! Level index
1739 !
1740    real(r8) co2fac(pcols,1)      ! co2 factor
1741    real(r8) alpha1(pcols)        ! stimulated emission term
1742    real(r8) alpha2(pcols)        ! stimulated emission term
1743    real(r8) rt(pcols)            ! reciprocal of local temperature
1744    real(r8) rsqrt(pcols)         ! reciprocal of sqrt of temp
1745 !
1746    real(r8) pbar(pcols)          ! mean pressure
1747    real(r8) dpnm(pcols)          ! difference in pressure
1748    real(r8) diff                 ! diffusivity factor
1749 !
1750 !--------------------------Data Statements------------------------------
1751 !
1752    data diff /1.66/
1753 !
1754 !-----------------------------------------------------------------------
1755 !
1756 !  Calculate path lengths for the trace gases at model top
1757 !
1758    do i = 1,ncol
1759       ucfc11(i,ntoplw) = 1.8 * cfc11(i,ntoplw) * pnm(i,ntoplw) * rga
1760       ucfc12(i,ntoplw) = 1.8 * cfc12(i,ntoplw) * pnm(i,ntoplw) * rga
1761       un2o0(i,ntoplw) = diff * 1.02346e5 * n2o(i,ntoplw) * pnm(i,ntoplw) * rga / sqrt(tnm(i,ntoplw))
1762       un2o1(i,ntoplw) = diff * 2.01909 * un2o0(i,ntoplw) * exp(-847.36/tnm(i,ntoplw))
1763       uch4(i,ntoplw)  = diff * 8.60957e4 * ch4(i,ntoplw) * pnm(i,ntoplw) * rga / sqrt(tnm(i,ntoplw))
1764       co2fac(i,1)     = diff * co2mmr * pnm(i,ntoplw) * rga
1765       alpha1(i) = (1.0 - exp(-1540.0/tnm(i,ntoplw)))**3.0/sqrt(tnm(i,ntoplw))
1766       alpha2(i) = (1.0 - exp(-1360.0/tnm(i,ntoplw)))**3.0/sqrt(tnm(i,ntoplw))
1767       uco211(i,ntoplw) = 3.42217e3 * co2fac(i,1) * alpha1(i) * exp(-1849.7/tnm(i,ntoplw))
1768       uco212(i,ntoplw) = 6.02454e3 * co2fac(i,1) * alpha1(i) * exp(-2782.1/tnm(i,ntoplw))
1769       uco213(i,ntoplw) = 5.53143e3 * co2fac(i,1) * alpha1(i) * exp(-3723.2/tnm(i,ntoplw))
1770       uco221(i,ntoplw) = 3.88984e3 * co2fac(i,1) * alpha2(i) * exp(-1997.6/tnm(i,ntoplw))
1771       uco222(i,ntoplw) = 3.67108e3 * co2fac(i,1) * alpha2(i) * exp(-3843.8/tnm(i,ntoplw))
1772       uco223(i,ntoplw) = 6.50642e3 * co2fac(i,1) * alpha2(i) * exp(-2989.7/tnm(i,ntoplw))
1773       bn2o0(i,ntoplw) = diff * 19.399 * pnm(i,ntoplw)**2.0 * n2o(i,ntoplw) * &
1774                    1.02346e5 * rga / (sslp*tnm(i,ntoplw))
1775       bn2o1(i,ntoplw) = bn2o0(i,ntoplw) * exp(-847.36/tnm(i,ntoplw)) * 2.06646e5
1776       bch4(i,ntoplw) = diff * 2.94449 * ch4(i,ntoplw) * pnm(i,ntoplw)**2.0 * rga * &
1777                   8.60957e4 / (sslp*tnm(i,ntoplw))
1778       uptype(i,ntoplw) = diff * qnm(i,ntoplw) * pnm(i,ntoplw)**2.0 *  &
1779                     exp(1800.0*(1.0/tnm(i,ntoplw) - 1.0/296.0)) * rga / sslp
1780    end do
1781 !
1782 ! Calculate trace gas path lengths through model atmosphere
1783 !
1784    do k = ntoplw,pver
1785       do i = 1,ncol
1786          rt(i) = 1./tnm(i,k)
1787          rsqrt(i) = sqrt(rt(i))
1788          pbar(i) = 0.5 * (pnm(i,k+1) + pnm(i,k)) / sslp
1789          dpnm(i) = (pnm(i,k+1) - pnm(i,k)) * rga
1790          alpha1(i) = diff * rsqrt(i) * (1.0 - exp(-1540.0/tnm(i,k)))**3.0
1791          alpha2(i) = diff * rsqrt(i) * (1.0 - exp(-1360.0/tnm(i,k)))**3.0
1792          ucfc11(i,k+1) = ucfc11(i,k) +  1.8 * cfc11(i,k) * dpnm(i)
1793          ucfc12(i,k+1) = ucfc12(i,k) +  1.8 * cfc12(i,k) * dpnm(i)
1794          un2o0(i,k+1) = un2o0(i,k) + diff * 1.02346e5 * n2o(i,k) * rsqrt(i) * dpnm(i)
1795          un2o1(i,k+1) = un2o1(i,k) + diff * 2.06646e5 * n2o(i,k) * &
1796                         rsqrt(i) * exp(-847.36/tnm(i,k)) * dpnm(i)
1797          uch4(i,k+1) = uch4(i,k) + diff * 8.60957e4 * ch4(i,k) * rsqrt(i) * dpnm(i)
1798          uco211(i,k+1) = uco211(i,k) + 1.15*3.42217e3 * alpha1(i) * &
1799                          co2mmr * exp(-1849.7/tnm(i,k)) * dpnm(i)
1800          uco212(i,k+1) = uco212(i,k) + 1.15*6.02454e3 * alpha1(i) * &
1801                          co2mmr * exp(-2782.1/tnm(i,k)) * dpnm(i)
1802          uco213(i,k+1) = uco213(i,k) + 1.15*5.53143e3 * alpha1(i) * &
1803                          co2mmr * exp(-3723.2/tnm(i,k)) * dpnm(i)
1804          uco221(i,k+1) = uco221(i,k) + 1.15*3.88984e3 * alpha2(i) * &
1805                          co2mmr * exp(-1997.6/tnm(i,k)) * dpnm(i)
1806          uco222(i,k+1) = uco222(i,k) + 1.15*3.67108e3 * alpha2(i) * &
1807                          co2mmr * exp(-3843.8/tnm(i,k)) * dpnm(i)
1808          uco223(i,k+1) = uco223(i,k) + 1.15*6.50642e3 * alpha2(i) * &
1809                          co2mmr * exp(-2989.7/tnm(i,k)) * dpnm(i)
1810          bn2o0(i,k+1) = bn2o0(i,k) + diff * 19.399 * pbar(i) * rt(i) &
1811                         * 1.02346e5 * n2o(i,k) * dpnm(i)
1812          bn2o1(i,k+1) = bn2o1(i,k) + diff * 19.399 * pbar(i) * rt(i) &
1813                         * 2.06646e5 * exp(-847.36/tnm(i,k)) * n2o(i,k)*dpnm(i)
1814          bch4(i,k+1) = bch4(i,k) + diff * 2.94449 * rt(i) * pbar(i) &
1815                        * 8.60957e4 * ch4(i,k) * dpnm(i)
1816          uptype(i,k+1) = uptype(i,k) + diff *qnm(i,k) * &
1817                          exp(1800.0*(1.0/tnm(i,k) - 1.0/296.0)) * pbar(i) * dpnm(i)
1818       end do
1819    end do
1820 !
1821    return
1822 end subroutine trcpth
1823
1824
1825
1826 subroutine aqsat(t       ,p       ,es      ,qs        ,ii      , &
1827                  ilen    ,kk      ,kstart  ,kend      )
1828 !-----------------------------------------------------------------------
1829 !
1830 ! Purpose:
1831 ! Utility procedure to look up and return saturation vapor pressure from
1832 ! precomputed table, calculate and return saturation specific humidity
1833 ! (g/g),for input arrays of temperature and pressure (dimensioned ii,kk)
1834 ! This routine is useful for evaluating only a selected region in the
1835 ! vertical.
1836 !
1837 ! Method:
1838 ! <Describe the algorithm(s) used in the routine.>
1839 ! <Also include any applicable external references.>
1840 !
1841 ! Author: J. Hack
1842 !
1843 !------------------------------Arguments--------------------------------
1844 !
1845 ! Input arguments
1846 !
1847    integer, intent(in) :: ii             ! I dimension of arrays t, p, es, qs
1848    integer, intent(in) :: kk             ! K dimension of arrays t, p, es, qs
1849    integer, intent(in) :: ilen           ! Length of vectors in I direction which
1850    integer, intent(in) :: kstart         ! Starting location in K direction
1851    integer, intent(in) :: kend           ! Ending location in K direction
1852    real(r8), intent(in) :: t(ii,kk)          ! Temperature
1853    real(r8), intent(in) :: p(ii,kk)          ! Pressure
1854 !
1855 ! Output arguments
1856 !
1857    real(r8), intent(out) :: es(ii,kk)         ! Saturation vapor pressure
1858    real(r8), intent(out) :: qs(ii,kk)         ! Saturation specific humidity
1859 !
1860 !---------------------------Local workspace-----------------------------
1861 !
1862    real(r8) omeps             ! 1 - 0.622
1863    integer i, k           ! Indices
1864 !
1865 !-----------------------------------------------------------------------
1866 !
1867    omeps = 1.0 - epsqs
1868    do k=kstart,kend
1869       do i=1,ilen
1870          es(i,k) = estblf(t(i,k))
1871 !
1872 ! Saturation specific humidity
1873 !
1874          qs(i,k) = epsqs*es(i,k)/(p(i,k) - omeps*es(i,k))
1875 !
1876 ! The following check is to avoid the generation of negative values
1877 ! that can occur in the upper stratosphere and mesosphere
1878 !
1879          qs(i,k) = min(1.0_r8,qs(i,k))
1880 !
1881          if (qs(i,k) < 0.0) then
1882             qs(i,k) = 1.0
1883             es(i,k) = p(i,k)
1884          end if
1885       end do
1886    end do
1887 !
1888    return
1889 end subroutine aqsat
1890 !===============================================================================
1891   subroutine cldefr(lchnk   ,ncol    ,pcols, pver, pverp, &
1892        landfrac,t       ,rel     ,rei     ,ps      ,pmid    , landm, icefrac, snowh)
1893 !-----------------------------------------------------------------------
1894 !
1895 ! Purpose:
1896 ! Compute cloud water and ice particle size
1897 !
1898 ! Method:
1899 ! use empirical formulas to construct effective radii
1900 !
1901 ! Author: J.T. Kiehl, B. A. Boville, P. Rasch
1902 !
1903 !-----------------------------------------------------------------------
1904
1905     implicit none
1906 !------------------------------Arguments--------------------------------
1907 !
1908 ! Input arguments
1909 !
1910     integer, intent(in) :: lchnk                 ! chunk identifier
1911     integer, intent(in) :: ncol                  ! number of atmospheric columns
1912     integer, intent(in) :: pcols, pver, pverp
1913
1914     real(r8), intent(in) :: landfrac(pcols)      ! Land fraction
1915     real(r8), intent(in) :: icefrac(pcols)       ! Ice fraction
1916     real(r8), intent(in) :: t(pcols,pver)        ! Temperature
1917     real(r8), intent(in) :: ps(pcols)            ! Surface pressure
1918     real(r8), intent(in) :: pmid(pcols,pver)     ! Midpoint pressures
1919     real(r8), intent(in) :: landm(pcols)
1920     real(r8), intent(in) :: snowh(pcols)         ! Snow depth over land, water equivalent (m)
1921 !
1922 ! Output arguments
1923 !
1924     real(r8), intent(out) :: rel(pcols,pver)      ! Liquid effective drop size (microns)
1925     real(r8), intent(out) :: rei(pcols,pver)      ! Ice effective drop size (microns)
1926 !
1927
1928 !++pjr
1929 ! following Kiehl
1930          call reltab(ncol, pcols, pver, t, landfrac, landm, icefrac, rel, snowh)
1931
1932 ! following Kristjansson and Mitchell
1933          call reitab(ncol, pcols, pver, t, rei)
1934 !--pjr
1935 !
1936 !
1937     return
1938   end subroutine cldefr
1939
1940
1941 subroutine background(lchnk, ncol, pint, pcols, pverr, pverrp, mmr)
1942 !-----------------------------------------------------------------------
1943 !
1944 ! Purpose:
1945 ! Set global mean tropospheric aerosol background (or tuning) field
1946 !
1947 ! Method:
1948 ! Specify aerosol mixing ratio.
1949 ! Aerosol mass mixing ratio
1950 ! is specified so that the column visible aerosol optical depth is a
1951 ! specified global number (tauback). This means that the actual mixing
1952 ! ratio depends on pressure thickness of the lowest three atmospheric
1953 ! layers near the surface.
1954 !
1955 !-----------------------------------------------------------------------
1956 !  use shr_kind_mod, only: r8 => shr_kind_r8
1957 !  use aer_optics, only: kbg,idxVIS
1958 !  use physconst, only: gravit
1959 !-----------------------------------------------------------------------
1960    implicit none
1961 !-----------------------------------------------------------------------
1962 !#include "ptrrgrid.h"
1963 !------------------------------Arguments--------------------------------
1964 !
1965 ! Input arguments
1966 !
1967    integer, intent(in) :: lchnk                 ! chunk identifier
1968    integer, intent(in) :: ncol                  ! number of atmospheric columns
1969    integer, intent(in) :: pcols,pverr,pverrp
1970
1971    real(r8), intent(in) :: pint(pcols,pverrp)   ! Interface pressure (mks)
1972 !
1973 ! Output arguments
1974 !
1975    real(r8), intent(out) :: mmr(pcols,pverr)    ! "background" aerosol mass mixing ratio
1976 !
1977 !---------------------------Local variables-----------------------------
1978 !
1979    integer i          ! Longitude index
1980    integer k          ! Level index
1981 !
1982    real(r8) mass2mmr  ! Factor to convert mass to mass mixing ratio
1983    real(r8) mass      ! Mass of "background" aerosol as specified by tauback
1984 !
1985 !-----------------------------------------------------------------------
1986 !
1987    do i=1,ncol
1988       mass2mmr =  gravmks / (pint(i,pverrp)-pint(i,pverrp-mxaerl))
1989       do k=1,pverr
1990 !
1991 ! Compute aerosol mass mixing ratio for specified levels (1.e3 factor is
1992 ! for units conversion of the extinction coefficiant from m2/g to m2/kg)
1993 !
1994         if ( k >= pverrp-mxaerl ) then
1995 ! kaervs is not consistent with the values in aer_optics
1996 ! this ?should? be changed.
1997 ! rhfac is also implemented differently
1998             mass = tauback / (1.e3 * kbg(idxVIS))
1999             mmr(i,k) = mass2mmr*mass
2000          else
2001             mmr(i,k) = 0._r8
2002          endif
2003 !
2004       enddo
2005    enddo
2006 !
2007    return
2008 end subroutine background
2009
2010 subroutine scale_aerosols(AEROSOLt, pcols, pver, ncol, lchnk, scale)
2011 !-----------------------------------------------------------------
2012 ! scale each species as determined by scale factors
2013 !-----------------------------------------------------------------
2014   integer, intent(in) :: ncol, lchnk ! number of columns and chunk index
2015   integer, intent(in) :: pcols, pver
2016   real(r8), intent(in) :: scale(naer_all) ! scale each aerosol by this amount
2017   real(r8), intent(inout) :: AEROSOLt(pcols, pver, naer_all) ! aerosols
2018   integer m
2019
2020   do m = 1, naer_all
2021      AEROSOLt(:ncol, :, m) = scale(m)*AEROSOLt(:ncol, :, m)
2022   end do
2023
2024   return
2025 end subroutine scale_aerosols
2026
2027 subroutine get_int_scales(scales)
2028   real(r8), intent(out)::scales(naer_all)  ! scale each aerosol by this amount
2029   integer i                                  ! index through species
2030
2031 !initialize
2032   scales = 1.
2033
2034   scales(idxBG) = 1._r8
2035   scales(idxSUL) = sulscl
2036   scales(idxSSLT) = ssltscl
2037
2038   do i = idxCARBONfirst, idxCARBONfirst+numCARBON-1
2039     scales(i) = carscl
2040   enddo
2041
2042   do i = idxDUSTfirst, idxDUSTfirst+numDUST-1
2043     scales(i) = dustscl
2044   enddo
2045
2046   scales(idxVOLC) = volcscl
2047
2048   return
2049 end subroutine get_int_scales
2050
2051 subroutine vert_interpolate (Match_ps, aerosolc, m_hybi, paerlev, naer_c, pint, n, AEROSOL_mmr, pcols, pver, pverp, ncol, c)
2052 !--------------------------------------------------------------------
2053 ! Input: match surface pressure, cam interface pressure,
2054 !        month index, number of columns, chunk index
2055 !
2056 ! Output: Aerosol mass mixing ratio (AEROSOL_mmr)
2057 !
2058 ! Method:
2059 !         interpolate column mass (cumulative) from match onto
2060 !           cam's vertical grid (pressure coordinate)
2061 !         convert back to mass mixing ratio
2062 !
2063 !--------------------------------------------------------------------
2064
2065 !  use physconst,     only: gravit
2066
2067    integer, intent(in)  :: paerlev,naer_c,pcols,pver,pverp
2068    real(r8), intent(out) :: AEROSOL_mmr(pcols,pver,naer)  ! aerosol mmr from MATCH
2069    real(r8), intent(in) :: Match_ps(pcols)                ! surface pressure at a particular month
2070    real(r8), intent(in) :: pint(pcols,pverp)              ! interface pressure from CAM
2071    real(r8), intent(in) :: aerosolc(pcols,paerlev,naer_c)
2072    real(r8), intent(in) :: m_hybi(paerlev)
2073
2074    integer, intent(in) :: ncol,c                          ! chunk index and number of columns
2075    integer, intent(in) :: n                               ! prv or nxt month index
2076 !
2077 ! Local workspace
2078 !
2079    integer m                           ! index to aerosol species
2080    integer kupper(pcols)               ! last upper bound for interpolation
2081    integer i, k, kk, kkstart, kount    ! loop vars for interpolation
2082    integer isv, ksv, msv               ! loop indices to save
2083
2084    logical bad                         ! indicates a bad point found
2085    logical lev_interp_comp             ! interpolation completed for a level
2086
2087    real(r8) AEROSOL(pcols,pverp,naer)  ! cumulative mass of aerosol in column beneath upper
2088                                        ! interface of level in column at particular month
2089    real(r8) dpl, dpu                   ! lower and upper intepolation factors
2090    real(r8) v_coord                    ! vertical coordinate
2091    real(r8) m_to_mmr                   ! mass to mass mixing ratio conversion factor
2092    real(r8) AER_diff                   ! temp var for difference between aerosol masses
2093 !
2094 ! Initialize index array
2095 !
2096    do i=1,ncol
2097       kupper(i) = 1
2098    end do
2099 !
2100 ! assign total mass to topmost level
2101 !
2102
2103    do i=1,ncol
2104    do m=1,naer
2105    AEROSOL(i,1,m) = AEROSOLc(i,1,m)
2106    enddo
2107    enddo
2108 !
2109 ! At every pressure level, interpolate onto that pressure level
2110 !
2111    do k=2,pver
2112 !
2113 ! Top level we need to start looking is the top level for the previous k
2114 ! for all longitude points
2115 !
2116       kkstart = paerlev
2117       do i=1,ncol
2118          kkstart = min0(kkstart,kupper(i))
2119       end do
2120       kount = 0
2121 !
2122 ! Store level indices for interpolation
2123 !
2124 ! for the pressure interpolation should be comparing
2125 ! pint(column,lev) with M_hybi(lev)*M_ps_cam_col(month,column,chunk)
2126 !
2127       lev_interp_comp = .false.
2128       do kk=kkstart,paerlev-1
2129          if(.not.lev_interp_comp) then
2130          do i=1,ncol
2131             v_coord = pint(i,k)
2132             if (M_hybi(kk)*Match_ps(i) .lt. v_coord .and. v_coord .le. M_hybi(kk+1)*Match_ps(i)) then
2133                kupper(i) = kk
2134                kount = kount + 1
2135             end if
2136          end do
2137 !
2138 ! If all indices for this level have been found, do the interpolation and
2139 ! go to the next level
2140 !
2141 ! Interpolate in pressure.
2142 !
2143          if (kount.eq.ncol) then
2144             do i=1,ncol
2145              do m=1,naer
2146                dpu = pint(i,k) - M_hybi(kupper(i))*Match_ps(i)
2147                dpl = M_hybi(kupper(i)+1)*Match_ps(i) - pint(i,k)
2148                AEROSOL(i,k,m) = &
2149                     (AEROSOLc(i,kupper(i)  ,m)*dpl + &
2150                      AEROSOLc(i,kupper(i)+1,m)*dpu)/(dpl + dpu)
2151              enddo
2152             enddo !i
2153             lev_interp_comp = .true.
2154          end if
2155          end if
2156       end do
2157 !
2158 ! If we've fallen through the kk=1,levsiz-1 loop, we cannot interpolate and
2159
2160 ! must extrapolate from the bottom or top pressure level for at least some
2161 ! of the longitude points.
2162 !
2163
2164       if(.not.lev_interp_comp) then
2165          do i=1,ncol
2166           do m=1,naer
2167             if (pint(i,k) .lt. M_hybi(1)*Match_ps(i)) then
2168                AEROSOL(i,k,m) =  AEROSOLc(i,1,m)
2169             else if (pint(i,k) .gt. M_hybi(paerlev)*Match_ps(i)) then
2170                AEROSOL(i,k,m) = 0.0
2171             else
2172                dpu = pint(i,k) - M_hybi(kupper(i))*Match_ps(i)
2173                dpl = M_hybi(kupper(i)+1)*Match_ps(i) - pint(i,k)
2174                AEROSOL(i,k,m) = &
2175                     (AEROSOLc(i,kupper(i)  ,m)*dpl + &
2176                      AEROSOLc(i,kupper(i)+1,m)*dpu)/(dpl + dpu)
2177             end if
2178           enddo
2179          end do
2180
2181          if (kount.gt.ncol) then
2182             call endrun ('VERT_INTERPOLATE: Bad data: non-monotonicity suspected in dependent variable')
2183          end if
2184       end if
2185    end do
2186
2187 !  call t_startf ('vi_checks')
2188 !
2189 ! aerosol mass beneath lowest interface (pverp) must be 0
2190 !
2191    AEROSOL(1:ncol,pverp,:) = 0.
2192 !
2193 ! Set mass in layer to zero whenever it is less than
2194 !   1.e-40 kg/m^2 in the layer
2195 !
2196    do m = 1, naer
2197       do k = 1, pver
2198          do i = 1, ncol
2199             if (AEROSOL(i,k,m) < 1.e-40_r8) AEROSOL(i,k,m) = 0.
2200          end do
2201       end do
2202    end do
2203 !
2204 ! Set mass in layer to zero whenever it is less than
2205 !   10^-15 relative to column total mass
2206 ! convert back to mass mixing ratios.
2207 ! exit if mmr is negative
2208 !
2209    do m = 1, naer
2210       do k = 1, pver
2211          do i = 1, ncol
2212             AER_diff = AEROSOL(i,k,m) - AEROSOL(i,k+1,m)
2213             if( abs(AER_diff) < 1e-15*AEROSOL(i,1,m)) then
2214                AER_diff = 0.
2215             end if
2216             m_to_mmr = gravmks / (pint(i,k+1)-pint(i,k))
2217             AEROSOL_mmr(i,k,m)= AER_diff * m_to_mmr
2218             if (AEROSOL_mmr(i,k,m) < 0) then
2219                write(6,*)'vert_interpolate: mmr < 0, m, col, lev, mmr',m, i, k, AEROSOL_mmr(i,k,m)
2220                write(6,*)'vert_interpolate: aerosol(k),(k+1)',AEROSOL(i,k,m),AEROSOL(i,k+1,m)
2221                write(6,*)'vert_interpolate: pint(k+1),(k)',pint(i,k+1),pint(i,k)
2222                write(6,*)'n,c',n,c
2223                CALL wrf_error_fatal('CLWRF-module_ra_cam_support. vert_interpolate: ERROR -- error -- Error of computation pint=NaN')
2224                call endrun()
2225             end if
2226          end do
2227       end do
2228    end do
2229
2230 !  call t_stopf ('vi_checks')
2231 !  call t_stopf ('vert_interpolate')
2232
2233    return
2234 end subroutine vert_interpolate
2235
2236
2237 !===============================================================================
2238   subroutine cldems(lchnk   ,ncol    ,pcols, pver, pverp, clwp    ,fice    ,rei     ,emis    )
2239 !-----------------------------------------------------------------------
2240 !
2241 ! Purpose:
2242 ! Compute cloud emissivity using cloud liquid water path (g/m**2)
2243 !
2244 ! Method:
2245 ! <Describe the algorithm(s) used in the routine.>
2246 ! <Also include any applicable external references.>
2247 !
2248 ! Author: J.T. Kiehl
2249 !
2250 !-----------------------------------------------------------------------
2251
2252     implicit none
2253 !------------------------------Parameters-------------------------------
2254 !
2255     real(r8) kabsl                  ! longwave liquid absorption coeff (m**2/g)
2256     parameter (kabsl = 0.090361)
2257 !
2258 !------------------------------Arguments--------------------------------
2259 !
2260 ! Input arguments
2261 !
2262     integer, intent(in) :: lchnk                   ! chunk identifier
2263     integer, intent(in) :: ncol                    ! number of atmospheric columns
2264     integer, intent(in) :: pcols, pver, pverp
2265
2266     real(r8), intent(in) :: clwp(pcols,pver)       ! cloud liquid water path (g/m**2)
2267     real(r8), intent(in) :: rei(pcols,pver)        ! ice effective drop size (microns)
2268     real(r8), intent(in) :: fice(pcols,pver)       ! fractional ice content within cloud
2269 !
2270 ! Output arguments
2271 !
2272     real(r8), intent(out) :: emis(pcols,pver)       ! cloud emissivity (fraction)
2273 !
2274 !---------------------------Local workspace-----------------------------
2275 !
2276     integer i,k                 ! longitude, level indices
2277     real(r8) kabs                   ! longwave absorption coeff (m**2/g)
2278     real(r8) kabsi                  ! ice absorption coefficient
2279 !
2280 !-----------------------------------------------------------------------
2281 !
2282     do k=1,pver
2283        do i=1,ncol
2284           kabsi = 0.005 + 1./rei(i,k)
2285           kabs = kabsl*(1.-fice(i,k)) + kabsi*fice(i,k)
2286           emis(i,k) = 1. - exp(-1.66*kabs*clwp(i,k))
2287        end do
2288     end do
2289 !
2290     return
2291   end subroutine cldems
2292
2293 !===============================================================================
2294   subroutine cldovrlap(lchnk   ,ncol    ,pcols, pver, pverp, pint    ,cld     ,nmxrgn  ,pmxrgn  )
2295 !-----------------------------------------------------------------------
2296 !
2297 ! Purpose:
2298 ! Partitions each column into regions with clouds in neighboring layers.
2299 ! This information is used to implement maximum overlap in these regions
2300 ! with random overlap between them.
2301 ! On output,
2302 !    nmxrgn contains the number of regions in each column
2303 !    pmxrgn contains the interface pressures for the lower boundaries of
2304 !           each region!
2305 ! Method:
2306
2307 !
2308 ! Author: W. Collins
2309 !
2310 !-----------------------------------------------------------------------
2311
2312     implicit none
2313 !
2314 ! Input arguments
2315 !
2316     integer, intent(in) :: lchnk                ! chunk identifier
2317     integer, intent(in) :: ncol                 ! number of atmospheric columns
2318     integer, intent(in) :: pcols, pver, pverp
2319
2320     real(r8), intent(in) :: pint(pcols,pverp)   ! Interface pressure
2321     real(r8), intent(in) :: cld(pcols,pver)     ! Fractional cloud cover
2322 !
2323 ! Output arguments
2324 !
2325     real(r8), intent(out) :: pmxrgn(pcols,pverp)! Maximum values of pressure for each
2326 !    maximally overlapped region.
2327 !    0->pmxrgn(i,1) is range of pressure for
2328 !    1st region,pmxrgn(i,1)->pmxrgn(i,2) for
2329 !    2nd region, etc
2330     integer nmxrgn(pcols)                    ! Number of maximally overlapped regions
2331 !
2332 !---------------------------Local variables-----------------------------
2333 !
2334     integer i                    ! Longitude index
2335     integer k                    ! Level index
2336     integer n                    ! Max-overlap region counter
2337
2338     real(r8) pnm(pcols,pverp)    ! Interface pressure
2339
2340     logical cld_found            ! Flag for detection of cloud
2341     logical cld_layer(pver)      ! Flag for cloud in layer
2342 !
2343 !------------------------------------------------------------------------
2344 !
2345
2346     do i = 1, ncol
2347        cld_found = .false.
2348        cld_layer(:) = cld(i,:) > 0.0_r8
2349        pmxrgn(i,:) = 0.0
2350        pnm(i,:)=pint(i,:)*10.
2351        n = 1
2352        do k = 1, pver
2353           if (cld_layer(k) .and.  .not. cld_found) then
2354              cld_found = .true.
2355           else if ( .not. cld_layer(k) .and. cld_found) then
2356              cld_found = .false.
2357              if (count(cld_layer(k:pver)) == 0) then
2358                 exit
2359              endif
2360              pmxrgn(i,n) = pnm(i,k)
2361              n = n + 1
2362           endif
2363        end do
2364        pmxrgn(i,n) = pnm(i,pverp)
2365        nmxrgn(i) = n
2366     end do
2367
2368     return
2369   end subroutine cldovrlap
2370
2371 !===============================================================================
2372   subroutine cldclw(lchnk   ,ncol    ,pcols, pver, pverp, zi      ,clwp    ,tpw     ,hl      )
2373 !-----------------------------------------------------------------------
2374 !
2375 ! Purpose:
2376 ! Evaluate cloud liquid water path clwp (g/m**2)
2377 !
2378 ! Method:
2379 ! <Describe the algorithm(s) used in the routine.>
2380 ! <Also include any applicable external references.>
2381 !
2382 ! Author: J.T. Kiehl
2383 !
2384 !-----------------------------------------------------------------------
2385
2386     implicit none
2387
2388 !
2389 ! Input arguments
2390 !
2391     integer, intent(in) :: lchnk                 ! chunk identifier
2392     integer, intent(in) :: ncol                  ! number of atmospheric columns
2393     integer, intent(in) :: pcols, pver, pverp
2394
2395     real(r8), intent(in) :: zi(pcols,pverp)      ! height at layer interfaces(m)
2396     real(r8), intent(in) :: tpw(pcols)           ! total precipitable water (mm)
2397 !
2398 ! Output arguments
2399 !
2400     real(r8) clwp(pcols,pver)     ! cloud liquid water path (g/m**2)
2401     real(r8) hl(pcols)            ! liquid water scale height
2402     real(r8) rhl(pcols)           ! 1/hl
2403
2404 !
2405 !---------------------------Local workspace-----------------------------
2406 !
2407     integer i,k               ! longitude, level indices
2408     real(r8) clwc0                ! reference liquid water concentration (g/m**3)
2409     real(r8) emziohl(pcols,pverp) ! exp(-zi/hl)
2410 !
2411 !-----------------------------------------------------------------------
2412 !
2413 ! Set reference liquid water concentration
2414 !
2415     clwc0 = 0.21
2416 !
2417 ! Diagnose liquid water scale height from precipitable water
2418 !
2419     do i=1,ncol
2420        hl(i)  = 700.0*log(max(tpw(i)+1.0_r8,1.0_r8))
2421        rhl(i) = 1.0/hl(i)
2422     end do
2423 !
2424 ! Evaluate cloud liquid water path (vertical integral of exponential fn)
2425 !
2426     do k=1,pverp
2427        do i=1,ncol
2428           emziohl(i,k) = exp(-zi(i,k)*rhl(i))
2429        end do
2430     end do
2431     do k=1,pver
2432        do i=1,ncol
2433           clwp(i,k) = clwc0*hl(i)*(emziohl(i,k+1) - emziohl(i,k))
2434        end do
2435     end do
2436 !
2437     return
2438   end subroutine cldclw
2439
2440
2441 !===============================================================================
2442   subroutine reltab(ncol, pcols, pver, t, landfrac, landm, icefrac, rel, snowh)
2443 !-----------------------------------------------------------------------
2444 !
2445 ! Purpose:
2446 ! Compute cloud water size
2447 !
2448 ! Method:
2449 ! analytic formula following the formulation originally developed by J. T. Kiehl
2450 !
2451 ! Author: Phil Rasch
2452 !
2453 !-----------------------------------------------------------------------
2454 !   use physconst,          only: tmelt
2455     implicit none
2456 !------------------------------Arguments--------------------------------
2457 !
2458 ! Input arguments
2459 !
2460     integer, intent(in) :: ncol
2461     integer, intent(in) :: pcols, pver
2462     real(r8), intent(in) :: landfrac(pcols)      ! Land fraction
2463     real(r8), intent(in) :: icefrac(pcols)       ! Ice fraction
2464     real(r8), intent(in) :: snowh(pcols)         ! Snow depth over land, water equivalent (m)
2465     real(r8), intent(in) :: landm(pcols)         ! Land fraction ramping to zero over ocean
2466     real(r8), intent(in) :: t(pcols,pver)        ! Temperature
2467
2468 !
2469 ! Output arguments
2470 !
2471     real(r8), intent(out) :: rel(pcols,pver)      ! Liquid effective drop size (microns)
2472 !
2473 !---------------------------Local workspace-----------------------------
2474 !
2475     integer i,k               ! Lon, lev indices
2476     real(r8) rliqland         ! liquid drop size if over land
2477     real(r8) rliqocean        ! liquid drop size if over ocean
2478     real(r8) rliqice          ! liquid drop size if over sea ice
2479 !
2480 !-----------------------------------------------------------------------
2481 !
2482     rliqocean = 14.0_r8
2483     rliqice   = 14.0_r8
2484     rliqland  = 8.0_r8
2485     do k=1,pver
2486        do i=1,ncol
2487 ! jrm Reworked effective radius algorithm
2488           ! Start with temperature-dependent value appropriate for continental air
2489           ! Note: findmcnew has a pressure dependence here
2490           rel(i,k) = rliqland + (rliqocean-rliqland) * min(1.0_r8,max(0.0_r8,(tmelt-t(i,k))*0.05))
2491           ! Modify for snow depth over land
2492           rel(i,k) = rel(i,k) + (rliqocean-rel(i,k)) * min(1.0_r8,max(0.0_r8,snowh(i)*10.))
2493           ! Ramp between polluted value over land to clean value over ocean.
2494           rel(i,k) = rel(i,k) + (rliqocean-rel(i,k)) * min(1.0_r8,max(0.0_r8,1.0-landm(i)))
2495           ! Ramp between the resultant value and a sea ice value in the presence of ice.
2496           rel(i,k) = rel(i,k) + (rliqice-rel(i,k)) * min(1.0_r8,max(0.0_r8,icefrac(i)))
2497 ! end jrm
2498        end do
2499     end do
2500   end subroutine reltab
2501 !===============================================================================
2502   subroutine reitab(ncol, pcols, pver, t, re)
2503     !
2504
2505     integer, intent(in) :: ncol, pcols, pver
2506     real(r8), intent(out) :: re(pcols,pver)
2507     real(r8), intent(in) :: t(pcols,pver)
2508     real(r8) corr
2509     integer i
2510     integer k
2511     integer index
2512     !
2513     do k=1,pver
2514        do i=1,ncol
2515           index = int(t(i,k)-179.)
2516           index = min(max(index,1),94)
2517           corr = t(i,k) - int(t(i,k))
2518           re(i,k) = retab(index)*(1.-corr)              &
2519                +retab(index+1)*corr
2520           !           re(i,k) = amax1(amin1(re(i,k),30.),10.)
2521        end do
2522     end do
2523     !
2524     return
2525   end subroutine reitab
2526
2527   function exp_interpol(x, f, y) result(g)
2528
2529     ! Purpose:
2530     !   interpolates f(x) to point y
2531     !   assuming f(x) = f(x0) exp a(x - x0)
2532     !   where a = ( ln f(x1) - ln f(x0) ) / (x1 - x0)
2533     !   x0 <= x <= x1
2534     !   assumes x is monotonically increasing
2535
2536     ! Author: D. Fillmore
2537
2538 !   use shr_kind_mod, only: r8 => shr_kind_r8
2539
2540     implicit none
2541
2542     real(r8), intent(in), dimension(:) :: x  ! grid points
2543     real(r8), intent(in), dimension(:) :: f  ! grid function values
2544     real(r8), intent(in) :: y                ! interpolation point
2545     real(r8) :: g                            ! interpolated function value
2546
2547     integer :: k  ! interpolation point index
2548     integer :: n  ! length of x
2549     real(r8) :: a
2550
2551     n = size(x)
2552
2553     ! find k such that x(k) < y =< x(k+1)
2554     ! set k = 1 if y <= x(1)  and  k = n-1 if y > x(n)
2555
2556     if (y <= x(1)) then
2557       k = 1
2558     else if (y >= x(n)) then
2559       k = n - 1
2560     else
2561       k = 1
2562       do while (y > x(k+1) .and. k < n)
2563         k = k + 1
2564       end do
2565     end if
2566
2567     ! interpolate
2568     a = (  log( f(k+1) / f(k) )  ) / ( x(k+1) - x(k) )
2569     g = f(k) * exp( a * (y - x(k)) )
2570
2571   end function exp_interpol
2572
2573   function lin_interpol(x, f, y) result(g)
2574
2575     ! Purpose:
2576     !   interpolates f(x) to point y
2577     !   assuming f(x) = f(x0) + a * (x - x0)
2578     !   where a = ( f(x1) - f(x0) ) / (x1 - x0)
2579     !   x0 <= x <= x1
2580     !   assumes x is monotonically increasing
2581
2582     ! Author: D. Fillmore
2583
2584 !   use shr_kind_mod, only: r8 => shr_kind_r8
2585
2586     implicit none
2587
2588     real(r8), intent(in), dimension(:) :: x  ! grid points
2589     real(r8), intent(in), dimension(:) :: f  ! grid function values
2590     real(r8), intent(in) :: y                ! interpolation point
2591     real(r8) :: g                            ! interpolated function value
2592
2593     integer :: k  ! interpolation point index
2594     integer :: n  ! length of x
2595     real(r8) :: a
2596
2597     n = size(x)
2598
2599     ! find k such that x(k) < y =< x(k+1)
2600     ! set k = 1 if y <= x(1)  and  k = n-1 if y > x(n)
2601
2602     if (y <= x(1)) then
2603       k = 1
2604     else if (y >= x(n)) then
2605       k = n - 1
2606     else
2607       k = 1
2608       do while (y > x(k+1) .and. k < n)
2609         k = k + 1
2610       end do
2611     end if
2612
2613     ! interpolate
2614     a = (  f(k+1) - f(k) ) / ( x(k+1) - x(k) )
2615     g = f(k) + a * (y - x(k))
2616
2617   end function lin_interpol
2618
2619   function lin_interpol2(x, f, y) result(g)
2620
2621     ! Purpose:
2622     !   interpolates f(x) to point y
2623     !   assuming f(x) = f(x0) + a * (x - x0)
2624     !   where a = ( f(x1) - f(x0) ) / (x1 - x0)
2625     !   x0 <= x <= x1
2626     !   assumes x is monotonically increasing
2627
2628     ! Author: D. Fillmore ::  J. Done changed from r8 to r4
2629
2630     implicit none
2631
2632     real, intent(in), dimension(:) :: x  ! grid points
2633     real, intent(in), dimension(:) :: f  ! grid function values
2634     real, intent(in) :: y                ! interpolation point
2635     real :: g                            ! interpolated function value
2636
2637     integer :: k  ! interpolation point index
2638     integer :: n  ! length of x
2639     real    :: a
2640
2641     n = size(x)
2642
2643     ! find k such that x(k) < y =< x(k+1)
2644     ! set k = 1 if y <= x(1)  and  k = n-1 if y > x(n)
2645
2646     if (y <= x(1)) then
2647       k = 1
2648     else if (y >= x(n)) then
2649       k = n - 1
2650     else
2651       k = 1
2652       do while (y > x(k+1) .and. k < n)
2653         k = k + 1
2654       end do
2655     end if
2656
2657     ! interpolate
2658     a = (  f(k+1) - f(k) ) / ( x(k+1) - x(k) )
2659     g = f(k) + a * (y - x(k))
2660
2661   end function lin_interpol2
2662
2663
2664 subroutine getfactors (cycflag, np1, cdayminus, cdayplus, cday, &
2665                        fact1, fact2)
2666 !---------------------------------------------------------------------------
2667 !
2668 ! Purpose: Determine time interpolation factors (normally for a boundary dataset)
2669 !          for linear interpolation.
2670 !
2671 ! Method:  Assume 365 days per year.  Output variable fact1 will be the weight to
2672 !          apply to data at calendar time "cdayminus", and fact2 the weight to apply
2673 !          to data at time "cdayplus".  Combining these values will produce a result
2674 !          valid at time "cday".  Output arguments fact1 and fact2 will be between
2675 !          0 and 1, and fact1 + fact2 = 1 to roundoff.
2676 !
2677 ! Author:  Jim Rosinski
2678 !
2679 !---------------------------------------------------------------------------
2680    implicit none
2681 !
2682 ! Arguments
2683 !
2684    logical, intent(in) :: cycflag             ! flag indicates whether dataset is being cycled yearly
2685
2686    integer, intent(in) :: np1                 ! index points to forward time slice matching cdayplus
2687
2688    real(r8), intent(in) :: cdayminus          ! calendar day of rearward time slice
2689    real(r8), intent(in) :: cdayplus           ! calendar day of forward time slice
2690    real(r8), intent(in) :: cday               ! calenar day to be interpolated to
2691    real(r8), intent(out) :: fact1             ! time interpolation factor to apply to rearward time slice
2692    real(r8), intent(out) :: fact2             ! time interpolation factor to apply to forward time slice
2693
2694 !  character(len=*), intent(in) :: str        ! string to be added to print in case of error (normally the callers name)
2695 !
2696 ! Local workspace
2697 !
2698    real(r8) :: deltat                         ! time difference (days) between cdayminus and cdayplus
2699    real(r8), parameter :: daysperyear = 365.  ! number of days in a year
2700 !
2701 ! Initial sanity checks
2702 !
2703 !  if (np1 == 1 .and. .not. cycflag) then
2704 !     call endrun ('GETFACTORS:'//str//' cycflag false and forward month index = Jan. not allowed')
2705 !  end if
2706
2707 !  if (np1 < 1) then
2708 !     call endrun ('GETFACTORS:'//str//' input arg np1 must be > 0')
2709 !  end if
2710
2711    if (cycflag) then
2712       if ((cday < 1.) .or. (cday > (daysperyear+1.))) then
2713          write(6,*) 'GETFACTORS:', ' bad cday=',cday
2714          call endrun ()
2715       end if
2716    else
2717       if (cday < 1.) then
2718          write(6,*) 'GETFACTORS:',  ' bad cday=',cday
2719          call endrun ()
2720       end if
2721    end if
2722 !
2723 ! Determine time interpolation factors.  Account for December-January
2724 ! interpolation if dataset is being cycled yearly.
2725 !
2726    if (cycflag .and. np1 == 1) then                     ! Dec-Jan interpolation
2727       deltat = cdayplus + daysperyear - cdayminus
2728       if (cday > cdayplus) then                         ! We are in December
2729          fact1 = (cdayplus + daysperyear - cday)/deltat
2730          fact2 = (cday - cdayminus)/deltat
2731       else                                              ! We are in January
2732          fact1 = (cdayplus - cday)/deltat
2733          fact2 = (cday + daysperyear - cdayminus)/deltat
2734       end if
2735    else
2736       deltat = cdayplus - cdayminus
2737       fact1 = (cdayplus - cday)/deltat
2738       fact2 = (cday - cdayminus)/deltat
2739    end if
2740
2741    if (.not. validfactors (fact1, fact2)) then
2742       write(6,*) 'GETFACTORS: ', ' bad fact1 and/or fact2=', fact1, fact2
2743       call endrun ()
2744    end if
2745
2746    return
2747 end subroutine getfactors
2748
2749 logical function validfactors (fact1, fact2)
2750 !---------------------------------------------------------------------------
2751 !
2752 ! Purpose: check sanity of time interpolation factors to within 32-bit roundoff
2753 !
2754 !---------------------------------------------------------------------------
2755    implicit none
2756
2757    real(r8), intent(in) :: fact1, fact2           ! time interpolation factors
2758
2759    validfactors = .true.
2760    if (abs(fact1+fact2-1.) > 1.e-6 .or. &
2761        fact1 > 1.000001 .or. fact1 < -1.e-6 .or. &
2762        fact2 > 1.000001 .or. fact2 < -1.e-6) then
2763
2764       validfactors = .false.
2765    end if
2766
2767    return
2768 end function validfactors
2769
2770 subroutine get_rf_scales(scales)
2771
2772   real(r8), intent(out)::scales(naer_all)  ! scale aerosols by this amount
2773
2774   integer i                                  ! loop index
2775
2776   scales(idxBG) = bgscl_rf
2777   scales(idxSUL) = sulscl_rf
2778   scales(idxSSLT) = ssltscl_rf
2779
2780   do i = idxCARBONfirst, idxCARBONfirst+numCARBON-1
2781     scales(i) = carscl_rf
2782   enddo
2783
2784   do i = idxDUSTfirst, idxDUSTfirst+numDUST-1
2785     scales(i) = dustscl_rf
2786   enddo
2787
2788   scales(idxVOLC) = volcscl_rf
2789
2790 end subroutine get_rf_scales
2791
2792 function psi(tpx,iband)
2793 !
2794 ! History: First version for Hitran 1996 (C/H/E)
2795 !          Current version for Hitran 2000 (C/LT/E)
2796 ! Short function for Hulst-Curtis-Godson temperature factors for
2797 !   computing effective H2O path
2798 ! Line data for H2O: Hitran 2000, plus H2O patches v11.0 for 1341 missing
2799 !                    lines between 500 and 2820 cm^-1.
2800 !                    See cfa-www.harvard.edu/HITRAN
2801 ! Isotopes of H2O: all
2802 ! Line widths: air-broadened only (self set to 0)
2803 ! Code for line strengths and widths: GENLN3
2804 ! Reference: Edwards, D.P., 1992: GENLN2, A General Line-by-Line Atmospheric
2805 !                     Transmittance and Radiance Model, Version 3.0 Description
2806 !                     and Users Guide, NCAR/TN-367+STR, 147 pp.
2807 !
2808 ! Note: functions have been normalized by dividing by their values at
2809 !       a path temperature of 160K
2810 !
2811 ! spectral intervals:
2812 !   1 = 0-800 cm^-1 and 1200-2200 cm^-1
2813 !   2 = 800-1200 cm^-1
2814 !
2815 ! Formulae: Goody and Yung, Atmospheric Radiation: Theoretical Basis,
2816 !           2nd edition, Oxford University Press, 1989.
2817 ! Psi: function for pressure along path
2818 !      eq. 6.30, p. 228
2819 !
2820    real(r8),intent(in):: tpx      ! path temperature
2821    integer, intent(in):: iband    ! band to process
2822    real(r8) psi                   ! psi for given band
2823    real(r8),parameter ::  psi_r0(nbands) = (/ 5.65308452E-01, -7.30087891E+01/)
2824    real(r8),parameter ::  psi_r1(nbands) = (/ 4.07519005E-03,  1.22199547E+00/)
2825    real(r8),parameter ::  psi_r2(nbands) = (/-1.04347237E-05, -7.12256227E-03/)
2826    real(r8),parameter ::  psi_r3(nbands) = (/ 1.23765354E-08,  1.47852825E-05/)
2827
2828    psi = (((psi_r3(iband) * tpx) + psi_r2(iband)) * tpx + psi_r1(iband)) * tpx + psi_r0(iband)
2829 end function psi
2830
2831 function phi(tpx,iband)
2832 !
2833 ! History: First version for Hitran 1996 (C/H/E)
2834 !          Current version for Hitran 2000 (C/LT/E)
2835 ! Short function for Hulst-Curtis-Godson temperature factors for
2836 !   computing effective H2O path
2837 ! Line data for H2O: Hitran 2000, plus H2O patches v11.0 for 1341 missing
2838 !                    lines between 500 and 2820 cm^-1.
2839 !                    See cfa-www.harvard.edu/HITRAN
2840 ! Isotopes of H2O: all
2841 ! Line widths: air-broadened only (self set to 0)
2842 ! Code for line strengths and widths: GENLN3
2843 ! Reference: Edwards, D.P., 1992: GENLN2, A General Line-by-Line Atmospheric
2844 !                     Transmittance and Radiance Model, Version 3.0 Description
2845 !                     and Users Guide, NCAR/TN-367+STR, 147 pp.
2846 !
2847 ! Note: functions have been normalized by dividing by their values at
2848 !       a path temperature of 160K
2849 !
2850 ! spectral intervals:
2851 !   1 = 0-800 cm^-1 and 1200-2200 cm^-1
2852 !   2 = 800-1200 cm^-1
2853 !
2854 ! Formulae: Goody and Yung, Atmospheric Radiation: Theoretical Basis,
2855 !           2nd edition, Oxford University Press, 1989.
2856 ! Phi: function for H2O path
2857 !      eq. 6.25, p. 228
2858 !
2859    real(r8),intent(in):: tpx      ! path temperature
2860    integer, intent(in):: iband    ! band to process
2861    real(r8) phi                   ! phi for given band
2862    real(r8),parameter ::  phi_r0(nbands) = (/ 9.60917711E-01, -2.21031342E+01/)
2863    real(r8),parameter ::  phi_r1(nbands) = (/ 4.86076751E-04,  4.24062610E-01/)
2864    real(r8),parameter ::  phi_r2(nbands) = (/-1.84806265E-06, -2.95543415E-03/)
2865    real(r8),parameter ::  phi_r3(nbands) = (/ 2.11239959E-09,  7.52470896E-06/)
2866
2867    phi = (((phi_r3(iband) * tpx) + phi_r2(iband)) * tpx + phi_r1(iband)) &
2868           * tpx + phi_r0(iband)
2869 end function phi
2870
2871 function fh2oself( temp )
2872 !
2873 ! Short function for H2O self-continuum temperature factor in
2874 !   calculation of effective H2O self-continuum path length
2875 !
2876 ! H2O Continuum: CKD 2.4
2877 ! Code for continuum: GENLN3
2878 ! Reference: Edwards, D.P., 1992: GENLN2, A General Line-by-Line Atmospheric
2879 !                     Transmittance and Radiance Model, Version 3.0 Description
2880 !                     and Users Guide, NCAR/TN-367+STR, 147 pp.
2881 !
2882 ! In GENLN, the temperature scaling of the self-continuum is handled
2883 !    by exponential interpolation/extrapolation from observations at
2884 !    260K and 296K by:
2885 !
2886 !         TFAC =  (T(IPATH) - 296.0)/(260.0 - 296.0)
2887 !         CSFFT = CSFF296*(CSFF260/CSFF296)**TFAC
2888 !
2889 ! For 800-1200 cm^-1, (CSFF260/CSFF296) ranges from ~2.1 to ~1.9
2890 !     with increasing wavenumber.  The ratio <CSFF260>/<CSFF296>,
2891 !     where <> indicates average over wavenumber, is ~2.07
2892 !
2893 ! fh2oself is (<CSFF260>/<CSFF296>)**TFAC
2894 !
2895    real(r8),intent(in) :: temp     ! path temperature
2896    real(r8) fh2oself               ! mean ratio of self-continuum at temp and 296K
2897
2898    fh2oself = 2.0727484**((296.0 - temp) / 36.0)
2899 end function fh2oself
2900
2901 ! from wv_saturation.F90
2902
2903 subroutine esinti(epslon  ,latvap  ,latice  ,rh2o    ,cpair   ,tmelt   )
2904 !-----------------------------------------------------------------------
2905 !
2906 ! Purpose:
2907 ! Initialize es lookup tables
2908 !
2909 ! Method:
2910 ! <Describe the algorithm(s) used in the routine.>
2911 ! <Also include any applicable external references.>
2912 !
2913 ! Author: J. Hack
2914 !
2915 !-----------------------------------------------------------------------
2916 !  use shr_kind_mod, only: r8 => shr_kind_r8
2917 !  use wv_saturation, only: gestbl
2918    implicit none
2919 !------------------------------Arguments--------------------------------
2920 !
2921 ! Input arguments
2922 !
2923    real(r8), intent(in) :: epslon          ! Ratio of h2o to dry air molecular weights
2924    real(r8), intent(in) :: latvap          ! Latent heat of vaporization
2925    real(r8), intent(in) :: latice          ! Latent heat of fusion
2926    real(r8), intent(in) :: rh2o            ! Gas constant for water vapor
2927    real(r8), intent(in) :: cpair           ! Specific heat of dry air
2928    real(r8), intent(in) :: tmelt           ! Melting point of water (K)
2929 !
2930 !---------------------------Local workspace-----------------------------
2931 !
2932    real(r8) tmn             ! Minimum temperature entry in table
2933    real(r8) tmx             ! Maximum temperature entry in table
2934    real(r8) trice           ! Trans range from es over h2o to es over ice
2935    logical ip           ! Ice phase (true or false)
2936 !
2937 !-----------------------------------------------------------------------
2938 !
2939 ! Specify control parameters first
2940 !
2941    tmn   = 173.16
2942    tmx   = 375.16
2943    trice =  20.00
2944    ip    = .true.
2945 !
2946 ! Call gestbl to build saturation vapor pressure table.
2947 !
2948    call gestbl(tmn     ,tmx     ,trice   ,ip      ,epslon  , &
2949                latvap  ,latice  ,rh2o    ,cpair   ,tmelt )
2950 !
2951    return
2952 end subroutine esinti
2953
2954 subroutine gestbl(tmn     ,tmx     ,trice   ,ip      ,epsil   , &
2955                   latvap  ,latice  ,rh2o    ,cpair   ,tmeltx   )
2956 !-----------------------------------------------------------------------
2957 !
2958 ! Purpose:
2959 ! Builds saturation vapor pressure table for later lookup procedure.
2960 !
2961 ! Method:
2962 ! Uses Goff & Gratch (1946) relationships to generate the table
2963 ! according to a set of free parameters defined below.  Auxiliary
2964 ! routines are also included for making rapid estimates (well with 1%)
2965 ! of both es and d(es)/dt for the particular table configuration.
2966 !
2967 ! Author: J. Hack
2968 !
2969 !-----------------------------------------------------------------------
2970 !  use pmgrid, only: masterproc
2971    implicit none
2972 !------------------------------Arguments--------------------------------
2973 !
2974 ! Input arguments
2975 !
2976    real(r8), intent(in) :: tmn           ! Minimum temperature entry in es lookup table
2977    real(r8), intent(in) :: tmx           ! Maximum temperature entry in es lookup table
2978    real(r8), intent(in) :: epsil         ! Ratio of h2o to dry air molecular weights
2979    real(r8), intent(in) :: trice         ! Transition range from es over range to es over ice
2980    real(r8), intent(in) :: latvap        ! Latent heat of vaporization
2981    real(r8), intent(in) :: latice        ! Latent heat of fusion
2982    real(r8), intent(in) :: rh2o          ! Gas constant for water vapor
2983    real(r8), intent(in) :: cpair         ! Specific heat of dry air
2984    real(r8), intent(in) :: tmeltx        ! Melting point of water (K)
2985 !
2986 !---------------------------Local variables-----------------------------
2987 !
2988    real(r8) t             ! Temperature
2989    real(r8) rgasv
2990    real(r8) cp
2991    real(r8) hlatf
2992    real(r8) ttrice
2993    real(r8) hlatv
2994    integer n          ! Increment counter
2995    integer lentbl     ! Calculated length of lookup table
2996    integer itype      ! Ice phase: 0 -> no ice phase
2997 !            1 -> ice phase, no transition
2998 !           -x -> ice phase, x degree transition
2999    logical ip         ! Ice phase logical flag
3000    logical icephs
3001 !
3002 !-----------------------------------------------------------------------
3003 !
3004 ! Set es table parameters
3005 !
3006    tmin   = tmn       ! Minimum temperature entry in table
3007    tmax   = tmx       ! Maximum temperature entry in table
3008    ttrice = trice     ! Trans. range from es over h2o to es over ice
3009    icephs = ip        ! Ice phase (true or false)
3010 !
3011 ! Set physical constants required for es calculation
3012 !
3013    epsqs  = epsil
3014    hlatv  = latvap
3015    hlatf  = latice
3016    rgasv  = rh2o
3017    cp     = cpair
3018    tmelt  = tmeltx
3019 !
3020    lentbl = INT(tmax-tmin+2.000001)
3021    if (lentbl .gt. plenest) then
3022       write(6,9000) tmax, tmin, plenest
3023       call endrun ('GESTBL')    ! Abnormal termination
3024    end if
3025 !
3026 ! Begin building es table.
3027 ! Check whether ice phase requested.
3028 ! If so, set appropriate transition range for temperature
3029 !
3030    if (icephs) then
3031       if (ttrice /= 0.0) then
3032          itype = -ttrice
3033       else
3034          itype = 1
3035       end if
3036    else
3037       itype = 0
3038    end if
3039 !
3040    t = tmin - 1.0
3041    do n=1,lentbl
3042       t = t + 1.0
3043       call gffgch(t,estbl(n),itype)
3044    end do
3045 !
3046    do n=lentbl+1,plenest
3047       estbl(n) = -99999.0
3048    end do
3049 !
3050 ! Table complete -- Set coefficients for polynomial approximation of
3051 ! difference between saturation vapor press over water and saturation
3052 ! pressure over ice for -ttrice < t < 0 (degrees C). NOTE: polynomial
3053 ! is valid in the range -40 < t < 0 (degrees C).
3054 !
3055 !                  --- Degree 5 approximation ---
3056 !
3057    pcf(1) =  5.04469588506e-01
3058    pcf(2) = -5.47288442819e+00
3059    pcf(3) = -3.67471858735e-01
3060    pcf(4) = -8.95963532403e-03
3061    pcf(5) = -7.78053686625e-05
3062 !
3063 !                  --- Degree 6 approximation ---
3064 !
3065 !-----pcf(1) =  7.63285250063e-02
3066 !-----pcf(2) = -5.86048427932e+00
3067 !-----pcf(3) = -4.38660831780e-01
3068 !-----pcf(4) = -1.37898276415e-02
3069 !-----pcf(5) = -2.14444472424e-04
3070 !-----pcf(6) = -1.36639103771e-06
3071 !
3072    if (masterproc) then
3073       write(6,*)' *** SATURATION VAPOR PRESSURE TABLE COMPLETED ***'
3074    end if
3075
3076    return
3077 !
3078 9000 format('GESTBL: FATAL ERROR *********************************',/, &
3079             ' TMAX AND TMIN REQUIRE A LARGER DIMENSION ON THE LENGTH', &
3080             ' OF THE SATURATION VAPOR PRESSURE TABLE ESTBL(PLENEST)',/, &
3081             ' TMAX, TMIN, AND PLENEST => ', 2f7.2, i3)
3082 !
3083 end subroutine gestbl
3084
3085 subroutine gffgch(t       ,es      ,itype   )
3086 !-----------------------------------------------------------------------
3087 !
3088 ! Purpose:
3089 ! Computes saturation vapor pressure over water and/or over ice using
3090 ! Goff & Gratch (1946) relationships.
3091 ! <Say what the routine does>
3092 !
3093 ! Method:
3094 ! T (temperature), and itype are input parameters, while es (saturation
3095 ! vapor pressure) is an output parameter.  The input parameter itype
3096 ! serves two purposes: a value of zero indicates that saturation vapor
3097 ! pressures over water are to be returned (regardless of temperature),
3098 ! while a value of one indicates that saturation vapor pressures over
3099 ! ice should be returned when t is less than freezing degrees.  If itype
3100 ! is negative, its absolute value is interpreted to define a temperature
3101 ! transition region below freezing in which the returned
3102 ! saturation vapor pressure is a weighted average of the respective ice
3103 ! and water value.  That is, in the temperature range 0 => -itype
3104 ! degrees c, the saturation vapor pressures are assumed to be a weighted
3105 ! average of the vapor pressure over supercooled water and ice (all
3106 ! water at 0 c; all ice at -itype c).  Maximum transition range => 40 c
3107 !
3108 ! Author: J. Hack
3109 !
3110 !-----------------------------------------------------------------------
3111 !  use shr_kind_mod, only: r8 => shr_kind_r8
3112 !  use physconst, only: tmelt
3113 !  use abortutils, only: endrun
3114
3115    implicit none
3116 !------------------------------Arguments--------------------------------
3117 !
3118 ! Input arguments
3119 !
3120    real(r8), intent(in) :: t          ! Temperature
3121 !
3122 ! Output arguments
3123 !
3124    integer, intent(inout) :: itype   ! Flag for ice phase and associated transition
3125
3126    real(r8), intent(out) :: es         ! Saturation vapor pressure
3127 !
3128 !---------------------------Local variables-----------------------------
3129 !
3130    real(r8) e1         ! Intermediate scratch variable for es over water
3131    real(r8) e2         ! Intermediate scratch variable for es over water
3132    real(r8) eswtr      ! Saturation vapor pressure over water
3133    real(r8) f          ! Intermediate scratch variable for es over water
3134    real(r8) f1         ! Intermediate scratch variable for es over water
3135    real(r8) f2         ! Intermediate scratch variable for es over water
3136    real(r8) f3         ! Intermediate scratch variable for es over water
3137    real(r8) f4         ! Intermediate scratch variable for es over water
3138    real(r8) f5         ! Intermediate scratch variable for es over water
3139    real(r8) ps         ! Reference pressure (mb)
3140    real(r8) t0         ! Reference temperature (freezing point of water)
3141    real(r8) term1      ! Intermediate scratch variable for es over ice
3142    real(r8) term2      ! Intermediate scratch variable for es over ice
3143    real(r8) term3      ! Intermediate scratch variable for es over ice
3144    real(r8) tr         ! Transition range for es over water to es over ice
3145    real(r8) ts         ! Reference temperature (boiling point of water)
3146    real(r8) weight     ! Intermediate scratch variable for es transition
3147    integer itypo   ! Intermediate scratch variable for holding itype
3148 !
3149 !-----------------------------------------------------------------------
3150 !
3151 ! Check on whether there is to be a transition region for es
3152 !
3153    if (itype < 0) then
3154       tr    = abs(float(itype))
3155       itypo = itype
3156       itype = 1
3157    else
3158       tr    = 0.0
3159       itypo = itype
3160    end if
3161    if (tr > 40.0) then
3162       write(6,900) tr
3163       call endrun ('GFFGCH')                ! Abnormal termination
3164    end if
3165 !
3166    if(t < (tmelt - tr) .and. itype == 1) go to 10
3167 !
3168 ! Water
3169 !
3170    ps = 1013.246
3171    ts = 373.16
3172    e1 = 11.344*(1.0 - t/ts)
3173    e2 = -3.49149*(ts/t - 1.0)
3174    f1 = -7.90298*(ts/t - 1.0)
3175    f2 = 5.02808*log10(ts/t)
3176    f3 = -1.3816*(10.0**e1 - 1.0)/10000000.0
3177    f4 = 8.1328*(10.0**e2 - 1.0)/1000.0
3178    f5 = log10(ps)
3179    f  = f1 + f2 + f3 + f4 + f5
3180    es = (10.0**f)*100.0
3181    eswtr = es
3182 !
3183    if(t >= tmelt .or. itype == 0) go to 20
3184 !
3185 ! Ice
3186 !
3187 10 continue
3188    t0    = tmelt
3189    term1 = 2.01889049/(t0/t)
3190    term2 = 3.56654*log(t0/t)
3191    term3 = 20.947031*(t0/t)
3192    es    = 575.185606e10*exp(-(term1 + term2 + term3))
3193 !
3194    if (t < (tmelt - tr)) go to 20
3195 !
3196 ! Weighted transition between water and ice
3197 !
3198    weight = min((tmelt - t)/tr,1.0_r8)
3199    es = weight*es + (1.0 - weight)*eswtr
3200 !
3201 20 continue
3202    itype = itypo
3203    return
3204 !
3205 900 format('GFFGCH: FATAL ERROR ******************************',/, &
3206            'TRANSITION RANGE FOR WATER TO ICE SATURATION VAPOR', &
3207            ' PRESSURE, TR, EXCEEDS MAXIMUM ALLOWABLE VALUE OF', &
3208            ' 40.0 DEGREES C',/, ' TR = ',f7.2)
3209 !
3210 end subroutine gffgch
3211
3212    real(r8) function estblf( td )
3213 !
3214 ! Saturation vapor pressure table lookup
3215 !
3216    real(r8), intent(in) :: td         ! Temperature for saturation lookup
3217 !
3218    real(r8) :: e       ! intermediate variable for es look-up
3219    real(r8) :: ai
3220    integer  :: i
3221 !
3222    e = max(min(td,tmax),tmin)   ! partial pressure
3223    i = int(e-tmin)+1
3224    ai = aint(e-tmin)
3225    estblf = (tmin+ai-e+1.)* &
3226             estbl(i)-(tmin+ai-e)* &
3227             estbl(i+1)
3228    end function estblf
3229
3230
3231 function findvalue(ix,n,ain,indxa)
3232 !-----------------------------------------------------------------------
3233 !
3234 ! Purpose:
3235 ! Subroutine for finding ix-th smallest value in the array
3236 ! The elements are rearranged so that the ix-th smallest
3237 ! element is in the ix place and all smaller elements are
3238 ! moved to the elements up to ix (with random order).
3239 !
3240 ! Algorithm: Based on the quicksort algorithm.
3241 !
3242 ! Author:       T. Craig
3243 !
3244 !-----------------------------------------------------------------------
3245 !  use shr_kind_mod, only: r8 => shr_kind_r8
3246    implicit none
3247 !
3248 ! arguments
3249 !
3250    integer, intent(in) :: ix                ! element to search for
3251    integer, intent(in) :: n                 ! total number of elements
3252    integer, intent(inout):: indxa(n)        ! array of integers
3253    real(r8), intent(in) :: ain(n)           ! array to search
3254 !
3255    integer findvalue                        ! return value
3256 !
3257 ! local variables
3258 !
3259    integer i,j
3260    integer il,im,ir
3261
3262    integer ia
3263    integer itmp
3264 !
3265 !---------------------------Routine-----------------------------
3266 !
3267    il=1
3268    ir=n
3269    do
3270       if (ir-il <= 1) then
3271          if (ir-il == 1) then
3272             if (ain(indxa(ir)) < ain(indxa(il))) then
3273                itmp=indxa(il)
3274                indxa(il)=indxa(ir)
3275                indxa(ir)=itmp
3276             endif
3277          endif
3278          findvalue=indxa(ix)
3279          return
3280       else
3281          im=(il+ir)/2
3282          itmp=indxa(im)
3283          indxa(im)=indxa(il+1)
3284          indxa(il+1)=itmp
3285          if (ain(indxa(il+1)) > ain(indxa(ir))) then
3286             itmp=indxa(il+1)
3287             indxa(il+1)=indxa(ir)
3288             indxa(ir)=itmp
3289          endif
3290          if (ain(indxa(il)) > ain(indxa(ir))) then
3291             itmp=indxa(il)
3292             indxa(il)=indxa(ir)
3293             indxa(ir)=itmp
3294          endif
3295          if (ain(indxa(il+1)) > ain(indxa(il))) then
3296             itmp=indxa(il+1)
3297             indxa(il+1)=indxa(il)
3298             indxa(il)=itmp
3299          endif
3300          i=il+1
3301          j=ir
3302          ia=indxa(il)
3303          do
3304             do
3305                i=i+1
3306                if (ain(indxa(i)) >= ain(ia)) exit
3307             end do
3308             do
3309                j=j-1
3310                if (ain(indxa(j)) <= ain(ia)) exit
3311             end do
3312             if (j < i) exit
3313             itmp=indxa(i)
3314             indxa(i)=indxa(j)
3315             indxa(j)=itmp
3316          end do
3317          indxa(il)=indxa(j)
3318          indxa(j)=ia
3319          if (j >= ix)ir=j-1
3320          if (j <= ix)il=i
3321       endif
3322    end do
3323 end function findvalue
3324
3325
3326 subroutine radini(gravx   ,cpairx  ,epsilox ,stebolx, pstdx )
3327 !-----------------------------------------------------------------------
3328 !
3329 ! Purpose:
3330 ! Initialize various constants for radiation scheme; note that
3331 ! the radiation scheme uses cgs units.
3332 !
3333 ! Method:
3334 ! <Describe the algorithm(s) used in the routine.>
3335 ! <Also include any applicable external references.>
3336 !
3337 ! Author: W. Collins (H2O parameterization) and J. Kiehl
3338 !
3339 !-----------------------------------------------------------------------
3340 !  use shr_kind_mod, only: r8 => shr_kind_r8
3341 !  use ppgrid,       only: pver, pverp
3342 !  use comozp,       only: cplos, cplol
3343 !  use pmgrid,       only: masterproc, plev, plevp
3344 !  use radae,        only: radaeini
3345 !  use physconst,    only: mwdry, mwco2
3346 #if ( defined SPMD )
3347 !   use mpishorthand
3348 #endif
3349    implicit none
3350
3351 !------------------------------Arguments--------------------------------
3352 !
3353 ! Input arguments
3354 !
3355    real, intent(in) :: gravx      ! Acceleration of gravity (MKS)
3356    real, intent(in) :: cpairx     ! Specific heat of dry air (MKS)
3357    real, intent(in) :: epsilox    ! Ratio of mol. wght of H2O to dry air
3358    real, intent(in) :: stebolx    ! Stefan-Boltzmann's constant (MKS)
3359    real(r8), intent(in) :: pstdx      ! Standard pressure (Pascals)
3360 !
3361 !---------------------------Local variables-----------------------------
3362 !
3363    integer k       ! Loop variable
3364
3365    real(r8) v0         ! Volume of a gas at stp (m**3/kmol)
3366    real(r8) p0         ! Standard pressure (pascals)
3367    real(r8) amd        ! Effective molecular weight of dry air (kg/kmol)
3368    real(r8) goz        ! Acceleration of gravity (m/s**2)
3369 !
3370 !-----------------------------------------------------------------------
3371 !
3372 ! Set general radiation consts; convert to cgs units where appropriate:
3373 !
3374    gravit  =  100.*gravx
3375    rga     =  1./gravit
3376    gravmks =  gravx
3377    cpair   =  1.e4*cpairx
3378    epsilo  =  epsilox
3379    sslp    =  1.013250e6
3380    stebol  =  1.e3*stebolx
3381    rgsslp  =  0.5/(gravit*sslp)
3382    dpfo3   =  2.5e-3
3383    dpfco2  =  5.0e-3
3384    dayspy  =  365.
3385    pie     =  4.*atan(1.)
3386 !
3387 ! Initialize ozone data.
3388 !
3389    v0  = 22.4136         ! Volume of a gas at stp (m**3/kmol)
3390    p0  = 0.1*sslp        ! Standard pressure (pascals)
3391    amd = 28.9644         ! Molecular weight of dry air (kg/kmol)
3392    goz = gravx           ! Acceleration of gravity (m/s**2)
3393 !
3394 ! Constants for ozone path integrals (multiplication by 100 for unit
3395 ! conversion to cgs from mks):
3396 !
3397    cplos = v0/(amd*goz)       *100.0
3398    cplol = v0/(amd*goz*p0)*0.5*100.0
3399 !
3400 ! Derived constants
3401 ! If the top model level is above ~90 km (0.1 Pa), set the top level to compute
3402 ! longwave cooling to about 80 km (1 Pa)
3403 ! WRF: assume top level > 0.1 mb
3404 !  if (hypm(1) .lt. 0.1) then
3405 !     do k = 1, pver
3406 !        if (hypm(k) .lt. 1.) ntoplw  = k
3407 !     end do
3408 !  else
3409       ntoplw = 1
3410 !  end if
3411 !   if (masterproc) then
3412 !     write (6,*) 'RADINI: ntoplw =',ntoplw, ' pressure:',hypm(ntoplw)
3413 !   endif
3414
3415    call radaeini( pstdx, mwdry, mwco2 )
3416    return
3417 end subroutine radini
3418
3419 subroutine oznini(ozmixm,pin,levsiz,num_months,XLAT,                &
3420                      ids, ide, jds, jde, kds, kde,                  &
3421                      ims, ime, jms, jme, kms, kme,                  &
3422                      its, ite, jts, jte, kts, kte)
3423 !
3424 ! This subroutine assumes uniform distribution of ozone concentration.
3425 ! It should be replaced by monthly climatology that varies latitudinally and vertically
3426 !
3427
3428 #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
3429   use module_dm, only: local_communicator
3430 #endif
3431       IMPLICIT NONE
3432
3433    INTEGER,      INTENT(IN   )    ::   ids,ide, jds,jde, kds,kde, &
3434                                        ims,ime, jms,jme, kms,kme, &
3435                                        its,ite, jts,jte, kts,kte
3436
3437    INTEGER,      INTENT(IN   )    ::   levsiz, num_months
3438
3439    REAL,  DIMENSION( ims:ime, jms:jme ), INTENT(IN   )  ::     XLAT
3440
3441    REAL,  DIMENSION( ims:ime, levsiz, jms:jme, num_months ),      &
3442           INTENT(OUT   ) ::                                  OZMIXM
3443
3444    REAL,  DIMENSION(levsiz), INTENT(OUT )  ::                   PIN
3445
3446 ! Local
3447    INTEGER, PARAMETER :: latsiz = 64
3448    INTEGER, PARAMETER :: lonsiz = 1
3449    INTEGER :: i, j, k, itf, jtf, ktf, m, pin_unit, lat_unit, oz_unit, ierr
3450    REAL    :: interp_pt
3451    CHARACTER*255 :: message
3452    real, pointer :: ozmixin(:,:,:,:), lat_ozone(:), plev(:)
3453 !   REAL, DIMENSION( lonsiz, levsiz, latsiz, num_months )    ::   &
3454 !                                                            OZMIXIN
3455
3456    logical, external :: wrf_dm_on_monitor
3457
3458    jtf=min0(jte,jde-1)
3459    ktf=min0(kte,kde-1)
3460    itf=min0(ite,ide-1)
3461
3462    if_have_ozone: if(.not.have_ozone) then
3463     call wrf_debug(1,'Do not have ozone.  Must read it in.')
3464     ! Allocate and set local aliases:
3465     levsiz_ozone_save=levsiz
3466     allocate(plev_ozone_save(levsiz),lat_ozone_save(latsiz))
3467     allocate(ozmixin_save(lonsiz, levsiz, latsiz, num_months))
3468     plev=>plev_ozone_save
3469     lat_ozone=>lat_ozone_save
3470     ozmixin=>ozmixin_save
3471 #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
3472     if_master: if(wrf_dm_on_monitor()) then
3473     call wrf_debug(1,'Master rank reads ozone.')
3474 #endif
3475 !-- read in ozone pressure data
3476
3477      WRITE(message,*)'num_months = ',num_months
3478      CALL wrf_debug(50,message)
3479
3480       pin_unit = 27
3481         OPEN(pin_unit, FILE='ozone_plev.formatted',FORM='FORMATTED',STATUS='OLD')
3482         do k = 1,levsiz
3483         READ (pin_unit,*)plev(k)
3484         end do
3485       close(27)
3486
3487       do k=1,levsiz
3488          plev(k) = plev(k)*100.
3489       end do
3490       pin=plev ! copy to grid array
3491
3492 !-- read in ozone lat data
3493
3494       lat_unit = 28
3495         OPEN(lat_unit, FILE='ozone_lat.formatted',FORM='FORMATTED',STATUS='OLD')
3496         do j = 1,latsiz
3497         READ (lat_unit,*)lat_ozone(j)
3498         end do
3499       close(28)
3500
3501
3502 !-- read in ozone data
3503
3504       oz_unit = 29
3505       OPEN(oz_unit, FILE='ozone.formatted',FORM='FORMATTED',STATUS='OLD')
3506
3507       do m=2,num_months
3508       do j=1,latsiz ! latsiz=64
3509       do k=1,levsiz ! levsiz=59
3510       do i=1,lonsiz ! lonsiz=1
3511         READ (oz_unit,*)ozmixin(i,k,j,m)
3512       enddo
3513       enddo
3514       enddo
3515       enddo
3516       close(29)
3517 #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
3518      endif if_master
3519      call wrf_debug(1,"Broadcast ozone to other ranks.")
3520 # if ( RWORDSIZE == DWORDSIZE )
3521      call MPI_Bcast(ozmixin,size(ozmixin),MPI_DOUBLE_PRECISION,0,local_communicator,ierr)
3522      call MPI_Bcast(pin,size(pin),MPI_DOUBLE_PRECISION,0,local_communicator,ierr)
3523      plev=pin
3524      call MPI_Bcast(lat_ozone,size(lat_ozone),MPI_DOUBLE_PRECISION,0,local_communicator,ierr)
3525 # else
3526      call MPI_Bcast(ozmixin,size(ozmixin),MPI_REAL,0,local_communicator,ierr)
3527      call MPI_Bcast(pin,size(pin),MPI_REAL,0,local_communicator,ierr)
3528      plev=pin
3529      call MPI_Bcast(lat_ozone,size(lat_ozone),MPI_REAL,0,local_communicator,ierr)
3530 # endif
3531 #endif
3532    else ! already read in ozone data
3533     ! Make sure, first:
3534     if(levsiz/=levsiz_ozone_save) then
3535 3081   format('Logic error in caller: levsiz=',I0,' but prior call used ',I0,'.')
3536        write(message,3081) levsiz,levsiz_ozone_save
3537        call wrf_error_fatal(message)
3538     endif
3539     if(.not.(associated(plev_ozone_save) .and. &
3540              associated(lat_ozone_save) .and. &
3541              associated(ozmixin_save))) then
3542        call wrf_error_fatal('Ozone save arrays are not allocated.')
3543     endif
3544     ! Recover the pointers to allocated data:
3545     plev=>plev_ozone_save
3546     lat_ozone=>lat_ozone_save
3547     ozmixin=>ozmixin_save
3548    endif if_have_ozone
3549
3550 !-- latitudinally interpolate ozone data (and extend longitudinally)
3551 !-- using function lin_interpol2(x, f, y) result(g)
3552 ! Purpose:
3553 !   interpolates f(x) to point y
3554 !   assuming f(x) = f(x0) + a * (x - x0)
3555 !   where a = ( f(x1) - f(x0) ) / (x1 - x0)
3556 !   x0 <= x <= x1
3557 !   assumes x is monotonically increasing
3558 !    real, intent(in), dimension(:) :: x  ! grid points
3559 !    real, intent(in), dimension(:) :: f  ! grid function values
3560 !    real, intent(in) :: y                ! interpolation point
3561 !    real :: g                            ! interpolated function value
3562 !---------------------------------------------------------------------------
3563
3564       do m=2,num_months
3565       do j=jts,jtf
3566       do k=1,levsiz
3567       do i=its,itf
3568          interp_pt=XLAT(i,j)
3569          ozmixm(i,k,j,m)=lin_interpol2(lat_ozone(:),ozmixin(1,k,:,m),interp_pt)
3570       enddo
3571       enddo
3572       enddo
3573       enddo
3574
3575 ! Old code for fixed ozone
3576
3577 !     pin(1)=70.
3578 !     DO k=2,levsiz
3579 !     pin(k)=pin(k-1)+16.
3580 !     ENDDO
3581
3582 !     DO k=1,levsiz
3583 !         pin(k) = pin(k)*100.
3584 !     end do
3585
3586 !     DO m=1,num_months
3587 !     DO j=jts,jtf
3588 !     DO i=its,itf
3589 !     DO k=1,2
3590 !      ozmixm(i,k,j,m)=1.e-6
3591 !     ENDDO
3592 !     DO k=3,levsiz
3593 !      ozmixm(i,k,j,m)=1.e-7
3594 !     ENDDO
3595 !     ENDDO
3596 !     ENDDO
3597 !     ENDDO
3598
3599 END SUBROUTINE oznini
3600
3601
3602 subroutine aerosol_init(m_psp,m_psn,m_hybi,aerosolcp,aerosolcn,paerlev,naer_c,shalf,pptop,    &
3603                      ids, ide, jds, jde, kds, kde,                  &
3604                      ims, ime, jms, jme, kms, kme,                  &
3605                      its, ite, jts, jte, kts, kte)
3606 !
3607 !  This subroutine assumes a uniform aerosol distribution in both time and space.
3608 !  It should be modified if aerosol data are available from WRF-CHEM or other sources
3609 !
3610       IMPLICIT NONE
3611
3612    INTEGER,      INTENT(IN   )    ::   ids,ide, jds,jde, kds,kde, &
3613                                        ims,ime, jms,jme, kms,kme, &
3614                                        its,ite, jts,jte, kts,kte
3615
3616    INTEGER,      INTENT(IN   )    ::   paerlev,naer_c
3617
3618    REAL,     intent(in)                        :: pptop
3619    REAL,     DIMENSION( kms:kme ), intent(in)  :: shalf
3620
3621    REAL,  DIMENSION( ims:ime, paerlev, jms:jme, naer_c ),      &
3622           INTENT(INOUT   ) ::                                  aerosolcn , aerosolcp
3623
3624    REAL,  DIMENSION(paerlev), INTENT(OUT )  ::                m_hybi
3625    REAL,  DIMENSION( ims:ime, jms:jme),  INTENT(OUT )  ::       m_psp,m_psn
3626
3627    REAL ::                                                      psurf
3628    real, dimension(29) :: hybi
3629    integer k ! index through vertical levels
3630
3631    INTEGER :: i, j, itf, jtf, ktf,m
3632
3633    data hybi/0, 0.0065700002014637, 0.0138600002974272, 0.023089999333024, &
3634     0.0346900001168251, 0.0491999983787537, 0.0672300010919571,      &
3635      0.0894500017166138, 0.116539999842644, 0.149159997701645,       &
3636     0.187830001115799, 0.232859998941422, 0.284209996461868,         &
3637     0.341369986534119, 0.403340011835098, 0.468600004911423,         &
3638     0.535290002822876, 0.601350009441376, 0.66482001543045,          &
3639     0.724009990692139, 0.777729988098145, 0.825269997119904,         &
3640     0.866419970989227, 0.901350021362305, 0.930540025234222,         &
3641     0.954590022563934, 0.974179983139038, 0.990000009536743, 1/
3642
3643    jtf=min0(jte,jde-1)
3644    ktf=min0(kte,kde-1)
3645    itf=min0(ite,ide-1)
3646
3647     do k=1,paerlev
3648       m_hybi(k)=hybi(k)
3649     enddo
3650
3651 !
3652 ! mxaerl = max number of levels (from bottom) for background aerosol
3653 ! Limit background aerosol height to regions below 900 mb
3654 !
3655
3656    psurf = 1.e05
3657    mxaerl = 0
3658 !  do k=pver,1,-1
3659    do k=kms,kme-1
3660 !     if (hypm(k) >= 9.e4) mxaerl = mxaerl + 1
3661       if (shalf(k)*psurf+pptop  >= 9.e4) mxaerl = mxaerl + 1
3662    end do
3663    mxaerl = max(mxaerl,1)
3664 !  if (masterproc) then
3665       write(6,*)'AEROSOLS:  Background aerosol will be limited to ', &
3666                 'bottom ',mxaerl,' model interfaces.'
3667 !               'bottom ',mxaerl,' model interfaces. Top interface is ', &
3668 !               hypi(pverp-mxaerl),' pascals'
3669 !  end if
3670
3671      DO j=jts,jtf
3672      DO i=its,itf
3673       m_psp(i,j)=psurf
3674       m_psn(i,j)=psurf
3675      ENDDO
3676      ENDDO
3677
3678      DO j=jts,jtf
3679      DO i=its,itf
3680      DO k=1,paerlev
3681 ! aerosolc arrays are upward cumulative (kg/m2) at each level
3682 ! Here we assume uniform vertical distribution (aerosolc linear with hybi)
3683       aerosolcp(i,k,j,idxSUL)=1.e-7*(1.-hybi(k))
3684       aerosolcn(i,k,j,idxSUL)=1.e-7*(1.-hybi(k))
3685       aerosolcp(i,k,j,idxSSLT)=1.e-22*(1.-hybi(k))
3686       aerosolcn(i,k,j,idxSSLT)=1.e-22*(1.-hybi(k))
3687       aerosolcp(i,k,j,idxDUSTfirst)=1.e-7*(1.-hybi(k))
3688       aerosolcn(i,k,j,idxDUSTfirst)=1.e-7*(1.-hybi(k))
3689       aerosolcp(i,k,j,idxDUSTfirst+1)=1.e-7*(1.-hybi(k))
3690       aerosolcn(i,k,j,idxDUSTfirst+1)=1.e-7*(1.-hybi(k))
3691       aerosolcp(i,k,j,idxDUSTfirst+2)=1.e-7*(1.-hybi(k))
3692       aerosolcn(i,k,j,idxDUSTfirst+2)=1.e-7*(1.-hybi(k))
3693       aerosolcp(i,k,j,idxDUSTfirst+3)=1.e-7*(1.-hybi(k))
3694       aerosolcn(i,k,j,idxDUSTfirst+3)=1.e-7*(1.-hybi(k))
3695       aerosolcp(i,k,j,idxOCPHO)=1.e-7*(1.-hybi(k))
3696       aerosolcn(i,k,j,idxOCPHO)=1.e-7*(1.-hybi(k))
3697       aerosolcp(i,k,j,idxBCPHO)=1.e-9*(1.-hybi(k))
3698       aerosolcn(i,k,j,idxBCPHO)=1.e-9*(1.-hybi(k))
3699       aerosolcp(i,k,j,idxOCPHI)=1.e-7*(1.-hybi(k))
3700       aerosolcn(i,k,j,idxOCPHI)=1.e-7*(1.-hybi(k))
3701       aerosolcp(i,k,j,idxBCPHI)=1.e-8*(1.-hybi(k))
3702       aerosolcn(i,k,j,idxBCPHI)=1.e-8*(1.-hybi(k))
3703      ENDDO
3704      ENDDO
3705      ENDDO
3706
3707      call aer_optics_initialize
3708
3709
3710 END subroutine aerosol_init
3711
3712   subroutine aer_optics_initialize
3713
3714 USE module_wrf_error
3715
3716 !   use shr_kind_mod, only: r8 => shr_kind_r8
3717 !   use pmgrid  ! masterproc is here
3718 !   use ioFileMod, only: getfil
3719
3720 !#if ( defined SPMD )
3721 !    use mpishorthand
3722 !#endif
3723     implicit none
3724
3725 !   include 'netcdf.inc'
3726
3727
3728     integer :: nrh_opac  ! number of relative humidity values for OPAC data
3729     integer :: nbnd      ! number of spectral bands, should be identical to nspint
3730     real(r8), parameter :: wgt_sscm = 6.0 / 7.0
3731     integer :: krh_opac  ! rh index for OPAC rh grid
3732     integer :: krh       ! another rh index
3733     integer :: ksz       ! dust size bin index
3734     integer :: kbnd      ! band index
3735
3736     real(r8) :: rh   ! local relative humidity variable
3737
3738     integer, parameter :: irh=8
3739     real(r8) :: rh_opac(irh)        ! OPAC relative humidity grid
3740     real(r8) :: ksul_opac(irh,nspint)    ! sulfate  extinction
3741     real(r8) :: wsul_opac(irh,nspint)    !          single scattering albedo
3742     real(r8) :: gsul_opac(irh,nspint)    !          asymmetry parameter
3743     real(r8) :: ksslt_opac(irh,nspint)   ! sea-salt
3744     real(r8) :: wsslt_opac(irh,nspint)
3745     real(r8) :: gsslt_opac(irh,nspint)
3746     real(r8) :: kssam_opac(irh,nspint)   ! sea-salt accumulation mode
3747     real(r8) :: wssam_opac(irh,nspint)
3748     real(r8) :: gssam_opac(irh,nspint)
3749     real(r8) :: ksscm_opac(irh,nspint)   ! sea-salt coarse mode
3750     real(r8) :: wsscm_opac(irh,nspint)
3751     real(r8) :: gsscm_opac(irh,nspint)
3752     real(r8) :: kcphil_opac(irh,nspint)  ! hydrophilic organic carbon
3753     real(r8) :: wcphil_opac(irh,nspint)
3754     real(r8) :: gcphil_opac(irh,nspint)
3755     real(r8) :: dummy(nspint)
3756
3757       LOGICAL                 :: opened
3758       LOGICAL , EXTERNAL      :: wrf_dm_on_monitor
3759
3760       CHARACTER*80 errmess
3761       INTEGER cam_aer_unit
3762       integer :: i
3763
3764 !   read aerosol optics data
3765
3766       IF ( wrf_dm_on_monitor() ) THEN
3767         DO i = 10,99
3768           INQUIRE ( i , OPENED = opened )
3769           IF ( .NOT. opened ) THEN
3770             cam_aer_unit = i
3771             GOTO 2010
3772           ENDIF
3773         ENDDO
3774         cam_aer_unit = -1
3775  2010   CONTINUE
3776       ENDIF
3777       CALL wrf_dm_bcast_bytes ( cam_aer_unit , IWORDSIZE )
3778       IF ( cam_aer_unit < 0 ) THEN
3779         CALL wrf_error_fatal ( 'module_ra_cam: aer_optics_initialize: Can not find unused fortran unit to read in lookup table.' )
3780       ENDIF
3781
3782         IF ( wrf_dm_on_monitor() ) THEN
3783           OPEN(cam_aer_unit,FILE='CAM_AEROPT_DATA',                  &
3784                FORM='UNFORMATTED',STATUS='OLD',ERR=9010)
3785           call wrf_debug(50,'reading CAM_AEROPT_DATA')
3786         ENDIF
3787
3788 #define DM_BCAST_MACRO(A) CALL wrf_dm_bcast_bytes ( A , size ( A ) * r8 )
3789
3790          IF ( wrf_dm_on_monitor() ) then
3791          READ (cam_aer_unit,ERR=9010) dummy
3792          READ (cam_aer_unit,ERR=9010) rh_opac
3793          READ (cam_aer_unit,ERR=9010) ksul_opac
3794          READ (cam_aer_unit,ERR=9010) wsul_opac
3795          READ (cam_aer_unit,ERR=9010) gsul_opac
3796          READ (cam_aer_unit,ERR=9010) kssam_opac
3797          READ (cam_aer_unit,ERR=9010) wssam_opac
3798          READ (cam_aer_unit,ERR=9010) gssam_opac
3799          READ (cam_aer_unit,ERR=9010) ksscm_opac
3800          READ (cam_aer_unit,ERR=9010) wsscm_opac
3801          READ (cam_aer_unit,ERR=9010) gsscm_opac
3802          READ (cam_aer_unit,ERR=9010) kcphil_opac
3803          READ (cam_aer_unit,ERR=9010) wcphil_opac
3804          READ (cam_aer_unit,ERR=9010) gcphil_opac
3805          READ (cam_aer_unit,ERR=9010) kcb
3806          READ (cam_aer_unit,ERR=9010) wcb
3807          READ (cam_aer_unit,ERR=9010) gcb
3808          READ (cam_aer_unit,ERR=9010) kdst
3809          READ (cam_aer_unit,ERR=9010) wdst
3810          READ (cam_aer_unit,ERR=9010) gdst
3811          READ (cam_aer_unit,ERR=9010) kbg
3812          READ (cam_aer_unit,ERR=9010) wbg
3813          READ (cam_aer_unit,ERR=9010) gbg
3814          READ (cam_aer_unit,ERR=9010) kvolc
3815          READ (cam_aer_unit,ERR=9010) wvolc
3816          READ (cam_aer_unit,ERR=9010) gvolc
3817          endif
3818
3819          DM_BCAST_MACRO(rh_opac)
3820          DM_BCAST_MACRO(ksul_opac)
3821          DM_BCAST_MACRO(wsul_opac)
3822          DM_BCAST_MACRO(gsul_opac)
3823          DM_BCAST_MACRO(kssam_opac)
3824          DM_BCAST_MACRO(wssam_opac)
3825          DM_BCAST_MACRO(gssam_opac)
3826          DM_BCAST_MACRO(ksscm_opac)
3827          DM_BCAST_MACRO(wsscm_opac)
3828          DM_BCAST_MACRO(gsscm_opac)
3829          DM_BCAST_MACRO(kcphil_opac)
3830          DM_BCAST_MACRO(wcphil_opac)
3831          DM_BCAST_MACRO(gcphil_opac)
3832          DM_BCAST_MACRO(kcb)
3833          DM_BCAST_MACRO(wcb)
3834          DM_BCAST_MACRO(gcb)
3835          DM_BCAST_MACRO(kvolc)
3836          DM_BCAST_MACRO(wvolc)
3837          DM_BCAST_MACRO(gvolc)
3838          DM_BCAST_MACRO(kdst)
3839          DM_BCAST_MACRO(wdst)
3840          DM_BCAST_MACRO(gdst)
3841          DM_BCAST_MACRO(kbg)
3842          DM_BCAST_MACRO(wbg)
3843          DM_BCAST_MACRO(gbg)
3844
3845          IF ( wrf_dm_on_monitor() ) CLOSE (cam_aer_unit)
3846
3847     ! map OPAC aerosol species onto CAM aerosol species
3848     ! CAM name             OPAC name
3849     ! sul   or SO4         = suso                  sulfate soluble
3850     ! sslt  or SSLT        = 1/7 ssam + 6/7 sscm   sea-salt accumulation/coagulation mode
3851     ! cphil or CPHI        = waso                  water soluble (carbon)
3852     ! cphob or CPHO        = waso @ rh = 0
3853     ! cb    or BCPHI/BCPHO = soot
3854
3855     ksslt_opac(:,:) = (1.0 - wgt_sscm) * kssam_opac(:,:) + wgt_sscm * ksscm_opac(:,:)
3856
3857     wsslt_opac(:,:) = ( (1.0 - wgt_sscm) * kssam_opac(:,:) * wssam_opac(:,:) &
3858                   + wgt_sscm * ksscm_opac(:,:) * wsscm_opac(:,:) ) &
3859                   / ksslt_opac(:,:)
3860
3861     gsslt_opac(:,:) = ( (1.0 - wgt_sscm) * kssam_opac(:,:) * wssam_opac(:,:) * gssam_opac(:,:) &
3862                   + wgt_sscm * ksscm_opac(:,:) * wsscm_opac(:,:) * gsscm_opac(:,:) ) &
3863                    / ( ksslt_opac(:,:) * wsslt_opac(:,:) )
3864
3865     do i=1,nspint
3866     kcphob(i) = kcphil_opac(1,i)
3867     wcphob(i) = wcphil_opac(1,i)
3868     gcphob(i) = gcphil_opac(1,i)
3869     end do
3870
3871     ! interpolate optical properties of hygrospopic aerosol species
3872     !   onto a uniform relative humidity grid
3873
3874     nbnd = nspint
3875
3876     do krh = 1, nrh
3877       rh = 1.0_r8 / nrh * (krh - 1)
3878       do kbnd = 1, nbnd
3879         ksul(krh, kbnd) = exp_interpol( rh_opac, &
3880           ksul_opac(:, kbnd) / ksul_opac(1, kbnd), rh ) * ksul_opac(1, kbnd)
3881         wsul(krh, kbnd) = lin_interpol( rh_opac, &
3882           wsul_opac(:, kbnd) / wsul_opac(1, kbnd), rh ) * wsul_opac(1, kbnd)
3883         gsul(krh, kbnd) = lin_interpol( rh_opac, &
3884           gsul_opac(:, kbnd) / gsul_opac(1, kbnd), rh ) * gsul_opac(1, kbnd)
3885         ksslt(krh, kbnd) = exp_interpol( rh_opac, &
3886           ksslt_opac(:, kbnd) / ksslt_opac(1, kbnd), rh ) * ksslt_opac(1, kbnd)
3887         wsslt(krh, kbnd) = lin_interpol( rh_opac, &
3888           wsslt_opac(:, kbnd) / wsslt_opac(1, kbnd), rh ) * wsslt_opac(1, kbnd)
3889         gsslt(krh, kbnd) = lin_interpol( rh_opac, &
3890           gsslt_opac(:, kbnd) / gsslt_opac(1, kbnd), rh ) * gsslt_opac(1, kbnd)
3891         kcphil(krh, kbnd) = exp_interpol( rh_opac, &
3892           kcphil_opac(:, kbnd) / kcphil_opac(1, kbnd), rh ) * kcphil_opac(1, kbnd)
3893         wcphil(krh, kbnd) = lin_interpol( rh_opac, &
3894           wcphil_opac(:, kbnd) / wcphil_opac(1, kbnd), rh ) * wcphil_opac(1, kbnd)
3895         gcphil(krh, kbnd) = lin_interpol( rh_opac, &
3896           gcphil_opac(:, kbnd) / gcphil_opac(1, kbnd), rh )  * gcphil_opac(1, kbnd)
3897       end do
3898     end do
3899
3900      RETURN
3901 9010 CONTINUE
3902      WRITE( errmess , '(A35,I4)' ) 'module_ra_cam: error reading unit ',cam_aer_unit
3903      CALL wrf_error_fatal(errmess)
3904
3905 END subroutine aer_optics_initialize
3906
3907
3908 subroutine radaeini( pstdx, mwdryx, mwco2x )
3909
3910 USE module_wrf_error
3911
3912 !
3913 ! Initialize radae module data
3914 !
3915 !
3916 ! Input variables
3917 !
3918    real(r8), intent(in) :: pstdx   ! Standard pressure (dynes/cm^2)
3919    real(r8), intent(in) :: mwdryx  ! Molecular weight of dry air
3920    real(r8), intent(in) :: mwco2x  ! Molecular weight of carbon dioxide
3921 !
3922 !      Variables for loading absorptivity/emissivity
3923 !
3924    integer ncid_ae                ! NetCDF file id for abs/ems file
3925
3926    integer pdimid                 ! pressure dimension id
3927    integer psize                  ! pressure dimension size
3928
3929    integer tpdimid                ! path temperature dimension id
3930    integer tpsize                 ! path temperature size
3931
3932    integer tedimid                ! emission temperature dimension id
3933    integer tesize                 ! emission temperature size
3934
3935    integer udimid                 ! u (H2O path) dimension id
3936    integer usize                  ! u (H2O path) dimension size
3937
3938    integer rhdimid                ! relative humidity dimension id
3939    integer rhsize                 ! relative humidity dimension size
3940
3941    integer    ah2onwid            ! var. id for non-wndw abs.
3942    integer    eh2onwid            ! var. id for non-wndw ems.
3943    integer    ah2owid             ! var. id for wndw abs. (adjacent layers)
3944    integer cn_ah2owid             ! var. id for continuum trans. for wndw abs.
3945    integer cn_eh2owid             ! var. id for continuum trans. for wndw ems.
3946    integer ln_ah2owid             ! var. id for line trans. for wndw abs.
3947    integer ln_eh2owid             ! var. id for line trans. for wndw ems.
3948
3949 !  character*(NF_MAX_NAME) tmpname! dummy variable for var/dim names
3950    character(len=256) locfn       ! local filename
3951    integer tmptype                ! dummy variable for variable type
3952    integer ndims                  ! number of dimensions
3953 !  integer dims(NF_MAX_VAR_DIMS)  ! vector of dimension ids
3954    integer natt                   ! number of attributes
3955 !
3956 ! Variables for setting up H2O table
3957 !
3958    integer t                     ! path temperature
3959    integer tmin                  ! mininum path temperature
3960    integer tmax                  ! maximum path temperature
3961    integer itype                 ! type of sat. pressure (=0 -> H2O only)
3962    integer i
3963    real(r8) tdbl
3964
3965       LOGICAL                 :: opened
3966       LOGICAL , EXTERNAL      :: wrf_dm_on_monitor
3967
3968       CHARACTER*80 errmess
3969       INTEGER cam_abs_unit
3970
3971 !
3972 ! Constants to set
3973 !
3974    p0     = pstdx
3975    amd    = mwdryx
3976    amco2  = mwco2x
3977 !
3978 ! Coefficients for h2o emissivity and absorptivity for overlap of H2O
3979 !    and trace gases.
3980 !
3981    c16  = coefj(3,1)/coefj(2,1)
3982    c17  = coefk(3,1)/coefk(2,1)
3983    c26  = coefj(3,2)/coefj(2,2)
3984    c27  = coefk(3,2)/coefk(2,2)
3985    c28  = .5
3986    c29  = .002053
3987    c30  = .1
3988    c31  = 3.0e-5
3989 !
3990 ! Initialize further longwave constants referring to far wing
3991 ! correction for overlap of H2O and trace gases; R&D refers to:
3992 !
3993 !            Ramanathan, V. and  P.Downey, 1986: A Nonisothermal
3994 !            Emissivity and Absorptivity Formulation for Water Vapor
3995 !            Journal of Geophysical Research, vol. 91., D8, pp 8649-8666
3996 !
3997    fwcoef = .1           ! See eq(33) R&D
3998    fwc1   = .30          ! See eq(33) R&D
3999    fwc2   = 4.5          ! See eq(33) and eq(34) in R&D
4000    fc1    = 2.6          ! See eq(34) R&D
4001
4002       IF ( wrf_dm_on_monitor() ) THEN
4003         DO i = 10,99
4004           INQUIRE ( i , OPENED = opened )
4005           IF ( .NOT. opened ) THEN
4006             cam_abs_unit = i
4007             GOTO 2010
4008           ENDIF
4009         ENDDO
4010         cam_abs_unit = -1
4011  2010   CONTINUE
4012       ENDIF
4013       CALL wrf_dm_bcast_bytes ( cam_abs_unit , IWORDSIZE )
4014       IF ( cam_abs_unit < 0 ) THEN
4015         CALL wrf_error_fatal ( 'module_ra_cam: radaeinit: Can not find unused fortran unit to read in lookup table.' )
4016       ENDIF
4017
4018         IF ( wrf_dm_on_monitor() ) THEN
4019           OPEN(cam_abs_unit,FILE='CAM_ABS_DATA',                  &
4020                FORM='UNFORMATTED',STATUS='OLD',ERR=9010)
4021           call wrf_debug(50,'reading CAM_ABS_DATA')
4022         ENDIF
4023
4024 #define DM_BCAST_MACRO(A) CALL wrf_dm_bcast_bytes ( A , size ( A ) * r8 )
4025
4026          IF ( wrf_dm_on_monitor() ) then
4027          READ (cam_abs_unit,ERR=9010) ah2onw
4028          READ (cam_abs_unit,ERR=9010) eh2onw
4029          READ (cam_abs_unit,ERR=9010) ah2ow
4030          READ (cam_abs_unit,ERR=9010) cn_ah2ow
4031          READ (cam_abs_unit,ERR=9010) cn_eh2ow
4032          READ (cam_abs_unit,ERR=9010) ln_ah2ow
4033          READ (cam_abs_unit,ERR=9010) ln_eh2ow
4034
4035          endif
4036
4037          DM_BCAST_MACRO(ah2onw)
4038          DM_BCAST_MACRO(eh2onw)
4039          DM_BCAST_MACRO(ah2ow)
4040          DM_BCAST_MACRO(cn_ah2ow)
4041          DM_BCAST_MACRO(cn_eh2ow)
4042          DM_BCAST_MACRO(ln_ah2ow)
4043          DM_BCAST_MACRO(ln_eh2ow)
4044
4045          IF ( wrf_dm_on_monitor() ) CLOSE (cam_abs_unit)
4046
4047 ! Set up table of H2O saturation vapor pressures for use in calculation
4048 !     effective path RH.  Need separate table from table in wv_saturation
4049 !     because:
4050 !     (1. Path temperatures can fall below minimum of that table; and
4051 !     (2. Abs/Emissivity tables are derived with RH for water only.
4052 !
4053       tmin = nint(min_tp_h2o)
4054       tmax = nint(max_tp_h2o)+1
4055       itype = 0
4056       do t = tmin, tmax
4057 !        call gffgch(dble(t),estblh2o(t-tmin),itype)
4058          tdbl = t
4059          call gffgch(tdbl,estblh2o(t-tmin),itype)
4060       end do
4061
4062      RETURN
4063 9010 CONTINUE
4064      WRITE( errmess , '(A35,I4)' ) 'module_ra_cam: error reading unit ',cam_abs_unit
4065      CALL wrf_error_fatal(errmess)
4066 end subroutine radaeini
4067
4068
4069 end MODULE module_ra_cam_support