wrfv2_fire/chem/module_ftuv_driver.F

   1 ! FTUV has been updated to conform with the v3 vertical indexing
   2 ! method with kte=kde-1. "Saves" were also added to the module level
   3 ! variables that are needed between calls since they could in theory be
   4 ! lost between calls to FTUV without them. In reality, I do not think
   5 ! this was happening due to the "use" statements, but better safe then
   6 ! sorry. Finally, a critical bug was fixed where the aerosol arrays were
   7 ! not set to 0 with non-SORGAM aerosol settings leading to junk J rates.
   8 !
   9 !!!! NOTE TO FTUV DEVELOPERS:
  10 !         The initialization routine needs to be changed to work using
  11 !         variable terrain height since the model top is not constant
  12 !         in time. It also is not constant in space, which is a bigger
  13 !         issue for most cases. At the very least, one needs to consider
  14 !         if the model top has been set at 50 vs. 100 mb, or something
  15 !         else entirely.
  16 !
  17 ! William.Gustafson@pnl.gov; 18-Sep-2008
  18 !
  19       module module_ftuv_driver
  20
  21       use module_wave_data, only : nw
  22
  23       implicit none
  24
  25 !     private variables
  26       integer, private, parameter :: nref = 51
  27       integer, private, save :: nz, kcon
  28       integer, private, save :: curjulday = 0
  29       real(8), private, save :: esfact = 1.0
  30
  31       real(8), private, save :: o3top = 2.97450e+16
  32       real(8), private, save :: o2top = 3.60906e+21
  33
  34       real(8), private, save :: zref(nref), tref(nref), airref(nref), o3ref(nref)
  35       real(8), private, save :: albedo(nw-1,2)
  36
  37       real(8), private, allocatable, save :: z_ph(:), tlev_ph(:), tlay_ph(:), airlev_ph(:)
  38       real(8), private, allocatable, save :: rh_ph(:), xlwc_ph(:)
  39       real(8), private, allocatable, save :: o3_ph(:)
  40       real(8), private, allocatable, save :: acb1_ph(:), acb2_ph(:)
  41       real(8), private, allocatable, save :: aoc1_ph(:), aoc2_ph(:)
  42       real(8), private, allocatable, save :: aant_ph(:), aso4_ph(:), asal_ph(:)
  43
  44       real(8), private, allocatable, save :: ftuv(:,:)
  45
  46
  47       data zref/ 00.0, 01.0, 02.0, 03.0, 04.0, 05.0, 06.0, 07.0, 08.0, 09.0,   &
  48                  10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 17.0, 18.0, 19.0,   &
  49                  20.0, 21.0, 22.0, 23.0, 24.0, 25.0, 26.0, 27.0, 28.0, 29.0,   &
  50                  30.0, 31.0, 32.0, 33.0, 34.0, 35.0, 36.0, 37.0, 38.0, 39.0,   &
  51                  40.0, 41.0, 42.0, 43.0, 44.0, 45.0, 46.0, 47.0, 48.0, 49.0,   &
  52                  50.0 /
  53
  54      data tref/ 288.150, 281.651, 275.154, 268.659, 262.166, 255.676, 249.187, &
  55                 242.700, 236.215, 229.733, 223.252, 216.774, 216.650, 216.650, &
  56                 216.650, 216.650, 216.650, 216.650, 216.650, 216.650, 216.650, &
  57                 217.581, 218.574, 219.567, 220.560, 221.552, 222.544, 223.536, &
  58                 224.527, 225.518, 226.509, 227.500, 228.490, 230.973, 233.743, &
  59                 236.513, 239.282, 242.050, 244.818, 247.584, 250.350, 253.114, &
  60                 255.878, 258.641, 261.403, 264.164, 266.925, 269.684, 270.650, &
  61                 270.650, 270.650 /
  62
  63      data airref/ 2.55E+19, 2.31E+19, 2.09E+19, 1.89E+19, 1.70E+19, 1.53E+19,  &
  64                   1.37E+19, 1.23E+19, 1.09E+19, 9.71E+18, 8.60E+18, 7.59E+18,  &
  65                   6.49E+18, 5.54E+18, 4.74E+18, 4.05E+18, 3.46E+18, 2.96E+18,  &
  66                   2.53E+18, 2.16E+18, 1.85E+18, 1.57E+18, 1.34E+18, 1.14E+18,  &
  67                   9.76E+17, 8.33E+17, 7.12E+17, 6.09E+17, 5.21E+17, 4.47E+17,  &
  68                   3.83E+17, 3.28E+17, 2.82E+17, 2.41E+17, 2.06E+17, 1.76E+17,  &
  69                   1.51E+17, 1.30E+17, 1.12E+17, 9.62E+16, 8.31E+16, 7.19E+16,  &
  70                   6.23E+16, 5.40E+16, 4.70E+16, 4.09E+16, 3.56E+16, 3.11E+16,  &
  71                   2.74E+16, 2.42E+16, 2.14E+16 /
  72
  73       data o3ref/ 8.00E+11, 7.39E+11, 6.90E+11, 6.22E+11, 5.80E+11, 5.67E+11,  &
  74                   5.70E+11, 5.86E+11, 6.50E+11, 8.23E+11, 1.13E+12, 1.57E+12,  &
  75                   2.02E+12, 2.24E+12, 2.35E+12, 2.57E+12, 2.95E+12, 3.47E+12,  &
  76                   4.04E+12, 4.49E+12, 4.77E+12, 4.88E+12, 4.86E+12, 4.73E+12,  &
  77                   4.54E+12, 4.32E+12, 4.03E+12, 3.65E+12, 3.24E+12, 2.85E+12,  &
  78                   2.52E+12, 2.26E+12, 2.03E+12, 1.80E+12, 1.58E+12, 1.40E+12,  &
  79                   1.22E+12, 1.04E+12, 8.73E+11, 7.31E+11, 6.07E+11, 4.95E+11,  &
  80                   3.98E+11, 3.18E+11, 2.54E+11, 2.03E+11, 1.62E+11, 1.30E+11,  &
  81                   1.04E+11, 8.27E+10, 6.61E+10/
  82
  83
  84 !     private parameters
  85       integer, private, parameter :: pid_no2   = 4
  86       integer, private, parameter :: pid_o31d  = 2
  87       integer, private, parameter :: pid_o33p  = 3
  88       integer, private, parameter :: pid_hono  = 12
  89       integer, private, parameter :: pid_hno3  = 13
  90       integer, private, parameter :: pid_hno4  = 14
  91       integer, private, parameter :: pid_no3o2 = 5
  92       integer, private, parameter :: pid_no3o  = 6
  93       integer, private, parameter :: pid_h2o2  = 11
  94       integer, private, parameter :: pid_ch2om = 16
  95       integer, private, parameter :: pid_ch2or = 15
  96       integer, private, parameter :: pid_ald   = 17
  97       integer, private, parameter :: pid_op1   = 24
  98       integer, private, parameter :: pid_op2   = 1    ! unknown
  99       integer, private, parameter :: pid_paa   = 16
 100       integer, private, parameter :: pid_ket   = 23
 101       integer, private, parameter :: pid_glya  = 21
 102       integer, private, parameter :: pid_glyb  = 21
 103       integer, private, parameter :: pid_mgly  = 22
 104       integer, private, parameter :: pid_dcb   = 1     ! unknown
 105       integer, private, parameter :: pid_onit  = 25
 106
 107
 108       contains
 109       subroutine ftuv_radm2_driver(                                    &
 110                               id, curr_secs, dtstep, config_flags,     &
 111                               gmt, julday,                             &
 112                               p_phy, t_phy, rho_phy, p8w, t8w,         &
 113                               xlat, xlong,    z_at_w,                  &
 114                               moist, chem,gd_cloud,gd_cloud2,          &
 115                               ph_no2,ph_o31d,ph_o33p,ph_hono,          &
 116                               ph_hno3,ph_hno4,ph_no3o2,ph_no3o,        &
 117                               ph_h2o2,ph_ch2om,ph_ch2or,ph_ald,        &
 118                               ph_op1,ph_op2,ph_paa,ph_ket,ph_glya,     &
 119                               ph_glyb,ph_mgly,ph_dcb,ph_onit,          &
 120                               ph_macr,ph_ch3coc2h5,                    &
 121                               ids,ide, jds,jde, kds,kde,               &
 122                               ims,ime, jms,jme, kms,kme,               &
 123                               its,ite, jts,jte, kts,kte                )
 124
 125       use module_configure
 126       use module_state_description
 127       use module_model_constants
 128
 129       use module_wave_data, only : tuv_jmax
 130
 131       implicit none
 132
 133 !     input arguments
 134       integer, intent(in   ) :: id, julday,                            &
 135                                 ids,ide, jds,jde, kds,kde,             &
 136                                 ims,ime, jms,jme, kms,kme,             &
 137                                 its,ite, jts,jte, kts,kte
 138
 139       real(kind=8), intent(in) :: curr_secs
 140
 141       real,    intent(in   ) :: dtstep, gmt
 142
 143       type(grid_config_rec_type),  intent(in   ) :: config_flags
 144
 145       real, dimension( ims:ime, kms:kme, jms:jme ),                    &
 146             intent(in   ) :: p_phy, t_phy, rho_phy, p8w, t8w, z_at_w
 147 !
 148 ! these arrays are for cloudwater/ice coming from grell convection, optional
 149       real, dimension( ims:ime, kms:kme, jms:jme ),                    &
 150             intent(in   ) :: gd_cloud,gd_cloud2
 151
 152       real, dimension( ims:ime, jms:jme ),                             &
 153             intent(in   ) :: xlat, xlong
 154
 155       real, dimension( ims:ime, kms:kme, jms:jme, num_moist ),         &
 156             intent(in   ) :: moist
 157
 158       real, dimension( ims:ime, kms:kme, jms:jme, num_chem ),          &
 159             intent(in   ) :: chem
 160
 161 !     output arguments
 162       real, dimension( ims:ime, kms:kme, jms:jme ),                    &
 163             intent(out  ) :: ph_no2, ph_o31d, ph_o33p, ph_hono,        &
 164                              ph_hno3, ph_hno4, ph_no3o2, ph_no3o,      &
 165                              ph_h2o2, ph_ch2om, ph_ch2or, ph_ald,      &
 166                              ph_op1, ph_op2, ph_paa, ph_ket, ph_glya,  &
 167                              ph_glyb, ph_mgly, ph_dcb,                 &
 168                              ph_onit, ph_macr, ph_ch3coc2h5
 169
 170 !     local parameters
 171       real(8), parameter :: xair_kg = 1.0e3/28.97*6.023e23*1.0e-6
 172
 173 !     local scalars
 174       integer :: i, j, k
 175       integer :: isorg
 176       real(8)    :: xtime, xhour, xmin, gmtp
 177       real(8)    :: dobson
 178       integer :: lu
 179
 180       real(8)    :: wrk, qs, es
 181
 182       logical, save :: first = .true.
 183
 184 !     local arrays
 185       real(8) :: temp(kts:kte), o3(kts:kte), air(kts:kte)
 186       real(8) :: zh(kts:kte), rh(kts:kte), xlwc(kts:kte)
 187
 188       real(8) :: acb1(kts:kte), acb2(kts:kte), aoc1(kts:kte), aoc2(kts:kte)
 189       real(8) :: aant(kts:kte), aso4(kts:kte), asal(kts:kte)
 190
 191       real(8) :: prate(kts:kte,tuv_jmax)
 192
 193 !     if the procedure is called first time, the initial subroutine is called
 194       if( first ) then
 195         call photo_inti( kte, 20.d0 )
 196         first = .false.
 197       endif
 198
 199 !     calculate the GMT
 200       isorg=0
 201       aer_select: SELECT CASE(config_flags%chem_opt)
 202          CASE (RADM2SORG,RADM2SORG_KPP,RACMSORG_KPP,RACMSORG_AQ)
 203              isorg=1
 204              CALL wrf_debug(15,'SORGAM aerosols initialization ')
 205          CASE DEFAULT
 206              CALL wrf_debug(15,'no aerosols initialization yet')
 207    END SELECT aer_select
 208
 209       xtime = curr_secs/60._8
 210       xhour = real( int(gmt+0.01,8) + int(xtime/60._8,8),8 )
 211       xmin  = 60.0*gmt + ( xtime - xhour*60._8 )
 212       gmtp  = mod( xhour, 24.0d0 )
 213       gmtp  = gmtp + xmin/60.0
 214
 215 !     set total ozone and land use type
 216       dobson = 265.0; lu = 1
 217
 218       xlwc(:) = 0.  !initialize in case no cloud water is present
 219
 220 !     get photolysis rates
 221       J_TILE_LOOP : do j = jts, jte
 222       I_TILE_LOOP : do i = its, ite
 223
 224         do k = kts, kte-1
 225
 226           temp(k+1) = t_phy(i,k,j)                            ! temperature(K)
 227           air(k+1)  = xair_kg*rho_phy(i,k,j)                  ! air num.(molecules/cm^3)
 228           o3(k+1)   = chem(i,k,j,p_o3)*1.0e-6                 ! ozone(VMR)
 229
 230           wrk       = t_phy(i,k,j) - 273.0
 231           es        = 6.11*10.0**(7.63*wrk/(241.9 + wrk))     ! Magnus EQ
 232           qs        = 0.622*es/(p_phy(i,k,j)*0.01)            ! sat mass mix (H2O)
 233           wrk       = moist(i,k,j,p_qv)/qs
 234           rh(k+1)   = max( 0.0d0, min( 1.0d0, wrk) )              ! relative huminity
 235
 236           if(p_qc.gt.1)xlwc(k+1) = 1.0e3*moist(i,k,j,p_qc)*rho_phy(i,k,j)  ! cloud water(g/m^3)
 237           if( xlwc(k+1) < 1.0e-6 ) xlwc(k+1) = 0.0
 238
 239           zh(k+1) = .5*(z_at_w(i,k+1,j)+z_at_w(i,k,j))*0.001 - z_at_w(i,kts,j)*0.001      ! height(km)
 240 !         zh(k+1) = z(i,k,j)*0.001 - z_at_w(i,1,j)*0.001      ! height(km)
 241
 242 !===============================================================================
 243 ! XUEXI
 244 !      ADD the aerosol effect on J
 245 !      chem is in the unit of ug/m3
 246 !      cb1 = p_ecj * 1.0
 247 !      cb2 = p_ecj * 0.0
 248 !      oc1 = p_orgpaj      ( Primary organic conc. Aitken mode)
 249 !          + p_p25aj       ( Unknown PM25 conc. Acc. mode )
 250 !      ant = p_no3aj       ( Nitrate conc. Acc. mode )
 251 !          + p_nh4aj       ( Ammonium conc. Acc. mode )
 252 !      so4 = p_so4aj       ( Sulfate conc. Acc. mode )
 253 !      sa1 = p_seas        ( Coarse soil-derived aero sols )
 254 !
 255 !=======================================================
 256           if(isorg == 1 ) then
 257
 258           acb1(k+1) = chem(i,k,j,p_ecj)
 259           acb2(k+1) = 0.0
 260 !         aoc1(k+1) = chem(i,k,j,p_orgpaj) + chem(i,k,j,p_orgbaj) + chem(i,k,j,p_orgaj)
 261           aoc1(k+1) = chem(i,k,j,p_orgpaj)
 262           aoc2(k+1) = 0.0
 263           aant(k+1) = chem(i,k,j,p_no3aj) + chem(i,k,j,p_nh4aj)
 264           aso4(k+1) = chem(i,k,j,p_so4aj)
 265           asal(k+1) = chem(i,k,j,p_seas)
 266
 267           else !Cannot do aerosol effect for non-SORGAM aerosol modules
 268                !This is a temporary hack; William.Gustafson@pnl.gov, 17-Sep-2008
 269
 270           acb1(k+1) = 0.0
 271           acb2(k+1) = 0.0
 272           aoc1(k+1) = 0.0
 273           aoc2(k+1) = 0.0
 274           aant(k+1) = 0.0
 275           aso4(k+1) = 0.0
 276           asal(k+1) = 0.0
 277
 278           endif
 279
 280 !==================================================================================
 281
 282         enddo
 283
 284         temp(1) = t8w(i,kts,j)
 285         air(1)  = xair_kg*( p8w(i,kts,j)/t8w(i,kts,j)/r_d )
 286         o3(1)   = o3(2)
 287         rh(1)   = rh(2)
 288         xlwc(1) = 0.
 289         zh(1)   = 0.
 290         acb1(1) = acb1(2)
 291         acb2(1) = acb2(2)
 292         aoc1(1) = aoc1(2)
 293         aoc2(1) = aoc2(2)
 294         aant(1) = aant(2)
 295         aso4(1) = aso4(2)
 296         asal(1) = asal(2)
 297
 298 !     smooth air density ...
 299         do k = kts+1, kte-1
 300           air(k) = .5*air(k) + .25*( air(k-1) + air(k+1) )
 301         enddo
 302
 303         call photo( kte,                  &
 304                     tuv_jmax,             &
 305                     julday,               &
 306                     gmtp,                 &
 307                     dble(xlat(i,j)),      &
 308                    -dble(xlong(i,j)),     &
 309                     dobson,               &
 310                     lu,                   &
 311                     zh,                   &
 312                     temp,                 &
 313                     air,                  &
 314                     rh,                   &
 315                     xlwc,                 &
 316                     o3,                   &
 317                     acb1,                 &
 318                     acb2,                 &
 319                     aoc1,                 &
 320                     aoc2,                 &
 321                     aant,                 &
 322                     aso4,                 &
 323                     asal,                 &
 324                     prate )
 325
 326         do k = kts, kte-1
 327           ph_no2(i,k,j)   = prate(k,pid_no2)*60.
 328           ph_o31d(i,k,j)  = prate(k,pid_o31d)*60.
 329           ph_o33p(i,k,j)  = prate(k,pid_o33p)*60.
 330           ph_hono(i,k,j)  = prate(k,pid_hono)*60.
 331           ph_hno3(i,k,j)  = prate(k,pid_hno3)*60.
 332           ph_hno4(i,k,j)  = prate(k,pid_hno4)*60.
 333           ph_no3o2(i,k,j) = prate(k,pid_no3o2)*60.
 334           ph_no3o(i,k,j)  = prate(k,pid_no3o)*60.
 335           ph_h2o2(i,k,j)  = prate(k,pid_h2o2)*60.
 336           ph_ch2om(i,k,j) = prate(k,pid_ch2om)*60.
 337           ph_ch2or(i,k,j) = prate(k,pid_ch2or)*60.
 338           ph_ald(i,k,j)   = prate(k,pid_ald)*60.
 339           ph_op1(i,k,j)   = prate(k,pid_op1)*60.
 340           ph_op2(i,k,j)   = prate(k,pid_op2)*60.
 341           ph_paa(i,k,j)   = prate(k,pid_paa)*60.
 342           ph_ket(i,k,j)   = prate(k,pid_ket)*60.
 343           ph_glya(i,k,j)  = prate(k,pid_glya)*60.
 344           ph_glyb(i,k,j)  = prate(k,pid_glyb)*60.
 345           ph_mgly(i,k,j)  = prate(k,pid_mgly)*60.
 346           ph_dcb(i,k,j)   = prate(k,pid_dcb)*60.
 347           ph_onit(i,k,j)  = prate(k,pid_onit)*60.
 348           ph_macr(i,k,j)  = prate(k,28) *60.  !
 349           ph_ch3coc2h5(i,k,j)=prate(k,23)*60. !
 350         enddo
 351
 352       enddo I_TILE_LOOP
 353       enddo J_TILE_LOOP
 354
 355       end subroutine ftuv_radm2_driver
 356
 357       subroutine photo_inti( nlev, ztopin )
 358
 359       use module_wave_data, only : nw, wl, tuv_jmax, wave_data_inti
 360       use module_ftuv_subs, only : nwint, wlint, schu_inti, inter_inti
 361
 362       implicit none
 363
 364 !     input arguments
 365       integer, intent(in) :: nlev
 366       real(8), intent(in) :: ztopin
 367
 368 !     local arguments
 369       integer :: k, kdis, nabv, iw
 370       real(8) :: ztop
 371
 372 !     change unit from m to km
 373       ztop = ztopin
 374
 375 !     first, set the top levels height, temperature and o3 profile
 376       do k = 1, nref
 377         if( ztop < zref(k) ) exit
 378       enddo
 379       if( k==1 .or. k>nref ) then
 380         print*, 'Cannot process these two conditions ...', ztop
 381         stop 'in photo_inti'
 382       endif
 383
 384       kcon = k + 1
 385       kdis = nref - kcon
 386       nabv = kdis/2 + 1
 387       if( mod(kdis,2) /= 0 ) then
 388         nabv = nabv + 1
 389       endif
 390       nz = nlev + nabv
 391
 392 !     allocate memory space
 393       allocate( z_ph(nz), tlev_ph(nz), tlay_ph(nz-1), airlev_ph(nz) )
 394       allocate( rh_ph(nz), xlwc_ph(nz) )
 395       allocate( o3_ph(nz) )
 396       allocate( acb1_ph(nz), acb2_ph(nz), aoc1_ph(nz), aoc2_ph(nz) )
 397       allocate( aant_ph(nz), aso4_ph(nz), asal_ph(nz) )
 398       allocate( ftuv(nz,tuv_jmax) )
 399
 400       z_ph(nlev+1:nz-1) = zref(kcon:nref:2)
 401       tlev_ph(nlev+1:nz-1) = tref(kcon:nref:2)
 402       airlev_ph(nlev+1:nz-1) = airref(kcon:nref:2)
 403       o3_ph(nlev+1:nz-1) = o3ref(kcon:nref:2)/airref(kcon:nref:2)  ! Changed to vmr
 404
 405       z_ph(nz) = zref(nref)
 406       tlev_ph(nz) = tref(nref)
 407       airlev_ph(nz) = airref(nref)
 408       o3_ph(nz) = o3ref(nref)/airref(nref)
 409
 410       tlay_ph(nlev+1:) = 0.5*( tlev_ph(nlev+1:nz-1) + tlev_ph(nlev+2:) )
 411
 412       rh_ph(:)   = 0.0
 413       xlwc_ph(:) = 0.0
 414       acb1_ph(:) = 0.0
 415       acb2_ph(:) = 0.0
 416       aoc1_ph(:) = 0.0
 417       aoc2_ph(:) = 0.0
 418       aant_ph(:) = 0.0
 419       aso4_ph(:) = 0.0
 420       asal_ph(:) = 0.0
 421
 422 !     second, set albedo for land
 423       do iw = 1, nw-1
 424         if( wl(iw)<400.0 ) albedo(iw,1) = 0.05
 425         if( (wl(iw)>=400.0 ) .and. (wl(iw)<450.0) ) albedo(iw,1) = 0.06
 426         if( (wl(iw)>=450.0 ) .and. (wl(iw)<500.0) ) albedo(iw,1) = 0.08
 427         if( (wl(iw)>=500.0 ) .and. (wl(iw)<550.0) ) albedo(iw,1) = 0.10
 428         if( (wl(iw)>=550.0 ) .and. (wl(iw)<600.0) ) albedo(iw,1) = 0.11
 429         if( (wl(iw)>=600.0 ) .and. (wl(iw)<640.0) ) albedo(iw,1) = 0.12
 430         if( (wl(iw)>=640.0 ) .and. (wl(iw)<660.0) ) albedo(iw,1) = 0.135
 431         if( wl(iw)>=660.0 ) albedo(iw,1) = 0.15
 432       end do
 433
 434 !     set albedo for sea water
 435       albedo(1:17,2) = (/ 0.8228, 0.8140, 0.8000, 0.7820, 0.7600,   &
 436                           0.7340, 0.7040, 0.6700, 0.6340, 0.4782,   &
 437                           0.3492, 0.3000, 0.2700, 0.2400, 0.1700,   &
 438                           0.0800, 0.0600 /)
 439
 440 !     third, intialize other modules
 441       call schu_inti
 442       call wave_data_inti( nz )
 443       call inter_inti( nw, wl, nwint, wlint )
 444
 445       print*, 'finish initialization ...'
 446
 447       end subroutine photo_inti
 448
 449       subroutine photo( nlev, njout, julday, gmtp, alat, along, o3toms, lu,     &
 450                         zin, tlevin, airlevin, rhin, xlwcin, o3in,              &
 451                         acb1in, acb2in, aoc1in, aoc2in,                         &
 452                         aantin, aso4in, asalin,                                 &
 453                         prate )
 454
 455       use module_ftuv_subs, only : sundis, calc_zenith, photoin
 456
 457       implicit none
 458
 459 !     input arguments
 460       integer, intent(in) :: nlev, njout
 461       integer, intent(in) :: julday
 462       real(8), intent(in) :: gmtp
 463       real(8), intent(in) :: alat, along
 464       real(8), intent(in) :: o3toms
 465       integer, intent(in) :: lu
 466
 467       real(8), intent(in) :: zin(nlev)
 468       real(8), intent(in) :: tlevin(nlev)
 469       real(8), intent(in) :: airlevin(nlev)
 470       real(8), intent(in) :: rhin(nlev)
 471       real(8), intent(in) :: xlwcin(nlev)
 472       real(8), intent(in) :: o3in(nlev)
 473       real(8), intent(in) :: acb1in(nlev)
 474       real(8), intent(in) :: acb2in(nlev)
 475       real(8), intent(in) :: aoc1in(nlev)
 476       real(8), intent(in) :: aoc2in(nlev)
 477       real(8), intent(in) :: aantin(nlev)
 478       real(8), intent(in) :: aso4in(nlev)
 479       real(8), intent(in) :: asalin(nlev)
 480
 481 !     output arguments
 482       real(8), intent(out) :: prate(nlev,njout)
 483
 484 !     local arugments
 485       integer :: k, n
 486       real(8) :: azim, zen
 487
 488 !-----------------------------------------------------------------------------
 489 ! Now xx-ftuv only considers the following 28 photolysis reactions ...
 490 !-----------------------------------------------------------------------------
 491 ! 1 o2 + hv -> o + o
 492 ! 2 o3 -> o2 + o(1d)
 493 ! 3 o3 -> o2 + o(3p)
 494 ! 4 no2 -> no + o(3p)
 495 ! 5 no3 -> no + o2
 496 ! 6 no3 -> no2 + o(3p)
 497 ! 7 n2o5 -> no3 + no + o(3p)
 498 ! 8 n2o5 -> no3 + no2
 499 ! 9 n2o + hv -> n2 + o(1d)
 500 ! 10 ho2 + hv -> oh + o
 501 ! 11 h2o2 -> 2 oh
 502 ! 12 hno2 -> oh + no
 503 ! 13 hno3 -> oh + no2
 504 ! 14 hno4 -> ho2 + no2
 505 ! 15 ch2o -> h + hco
 506 ! 16 ch2o -> h2 + co
 507 ! 17 ch3cho -> ch3 + hco
 508 ! 18 ch3cho -> ch4 + co
 509 ! 19 ch3cho -> ch3co + h
 510 ! 20 c2h5cho -> c2h5 + hco
 511 ! 21 chocho -> products
 512 ! 22 ch3cocho -> products
 513 ! 23 ch3coch3
 514 ! 24 ch3ooh -> ch3o + oh
 515 ! 25 ch3ono2 -> ch3o+no2
 516 ! 26 pan + hv -> products
 517 ! 27 mvk + hv -> products
 518 ! 28 macr + hv -> products
 519 !-----------------------------------------------------------------------------
 520       if( curjulday /= julday ) then
 521         curjulday = julday
 522         call sundis( curjulday, esfact )
 523       endif
 524
 525 !     solar zenith angle calculation:
 526       call calc_zenith( alat, along, julday, gmtp, azim, zen )
 527       if( zen == 90.0 ) then
 528         zen = 89.0
 529       elseif( zen >= 95.0 ) then
 530         prate(:,:) = 0.0
 531         return
 532       endif
 533
 534 !     assgin vertical profiles
 535       z_ph(1:nlev) = zin(1:nlev)
 536       tlev_ph(1:nlev) = tlevin(1:nlev)
 537       airlev_ph(1:nlev) = airlevin(1:nlev)
 538       rh_ph(1:nlev) = rhin(1:nlev)
 539       xlwc_ph(1:nlev) = xlwcin(1:nlev)
 540       o3_ph(1:nlev) = o3in(1:nlev)
 541       acb1_ph(1:nlev) = acb1in(1:nlev)
 542       acb2_ph(1:nlev) = acb2in(1:nlev)
 543       aoc1_ph(1:nlev) = aoc1in(1:nlev)
 544       aoc2_ph(1:nlev) = aoc2in(1:nlev)
 545       aant_ph(1:nlev) = aantin(1:nlev)
 546       aso4_ph(1:nlev) = aso4in(1:nlev)
 547       asal_ph(1:nlev) = asalin(1:nlev)
 548
 549       tlay_ph(1:nlev-1) = 0.5*( tlev_ph(1:nlev-1) + tlev_ph(2:nlev) )
 550       tlay_ph(nlev) = 0.5*( tlev_ph(nlev) + tref(kcon) )
 551
 552       call photoin( nz,             &
 553                     zen,            &
 554                     o3toms,         &
 555                     esfact,         &
 556                     o3top,          &
 557                     o2top,          &
 558                     albedo(:,lu),   &
 559                     z_ph,           &
 560                     tlev_ph,        &
 561                     tlay_ph,        &
 562                     airlev_ph,      &
 563                     rh_ph,          &
 564                     xlwc_ph,        &
 565                     o3_ph,          &
 566                     acb1_ph,        &
 567                     acb2_ph,        &
 568                     aoc1_ph,        &
 569                     aoc2_ph,        &
 570                     aant_ph,        &
 571                     aso4_ph,        &
 572                     asal_ph,        &
 573                     ftuv )
 574
 575       do k = 1, nlev-1
 576         prate(k,1:njout) = ftuv(nz-k,1:njout)
 577       enddo
 578
 579       end subroutine photo
 580
 581       end module module_ftuv_driver