chem/module_phot_fastj.F

   1 !************************************************************************
   2 ! This computer software was prepared by Battelle Memorial Institute,
   3 ! hereinafter the Contractor, under Contract No. DE-AC05-76RL0 1830 with
   4 ! the Department of Energy (DOE). NEITHER THE GOVERNMENT NOR THE
   5 ! CONTRACTOR MAKES ANY WARRANTY, EXPRESS OR IMPLIED, OR ASSUMES ANY
   6 ! LIABILITY FOR THE USE OF THIS SOFTWARE.
   7 !
   8 ! Photolysis Option:  Fast-J
   9 ! * Primary investigators: Elaine G. Chapman and James C. Barnard
  10 ! * Co-investigators: Jerome D. Fast, William I. Gustafson Jr.
  11 ! Last update: February 2009
  12 !
  13 ! Contact:
  14 ! Jerome D. Fast, PhD
  15 ! Staff Scientist
  16 ! Pacific Northwest National Laboratory
  17 ! P.O. Box 999, MSIN K9-30
  18 ! Richland, WA, 99352
  19 ! Phone: (509) 372-6116
  20 ! Email: Jerome.Fast@pnl.gov
  21 !
  22 ! The original Fast-J code was provided by Oliver Wild (Univ. of Calif.
  23 ! Irvine); however, substantial modifications were necessary to make it
  24 ! compatible with WRF-Chem and to include the effect of prognostic
  25 ! aerosols on photolysis rates.
  26 !
  27 ! Please report any bugs or problems to Jerome Fast, the WRF-chem
  28 ! implmentation team leader for PNNL.
  29 !
  30 ! References:
  31 ! * Wild, O., X. Zhu, M.J. Prather, (2000), Accurate simulation of in-
  32 !   and below cloud photolysis in tropospheric chemical models, J. Atmos.
  33 !   Chem., 37, 245-282.
  34 ! * Barnard, J.C., E.G. Chapman, J.D. Fast, J.R. Schmelzer, J.R.
  35 !   Schulsser, and R.E. Shetter (2004), An evaluation of the FAST-J
  36 !   photolysis model for predicting nitrogen dioxide photolysis rates
  37 !   under clear and cloudy sky conditions, Atmos. Environ., 38,
  38 !   3393-3403.
  39 ! * Fast, J.D., W.I. Gustafson Jr., R.C. Easter, R.A. Zaveri, J.C.
  40 !   Barnard, E.G. Chapman, G.A. Grell, and S.E. Peckham (2005), Evolution
  41 !   of ozone, particulates, and aerosol direct radiative forcing in the
  42 !   vicinity of Houston using a fully-coupled meteorology-chemistry-
  43 !   aerosol model, J. Geophys. Res., 111, D21305,
  44 !   doi:10.1029/2005JD006721.
  45 ! * Barnard, J.C., J.D. Fast, G. Paredes-Miranda, W.P. Arnott, and
  46 !   A. Laskin (2010), Technical note: evaluation of the WRF-Chem
  47 !   "aerosol chemical to aerosol optical properties" module using data
  48 !   from the MILAGRO campaign, Atmos. Chem. Phys., 10, 7325-7340,
  49 !   doi:10.5194/acp-10-7325-2010.
  50 !
  51 ! Contact Jerome Fast for updates on the status of manuscripts under
  52 ! review.
  53 !
  54 ! Additional information:
  55 ! *  www.pnl.gov/atmospheric/research/wrf-chem
  56 !
  57 ! Support:
  58 ! Funding for adapting Fast-J was provided by the U.S. Department of
  59 ! Energy under the auspices of Atmospheric Science Program of the Office
  60 ! of Biological and Environmental Research the PNNL Laboratory Research
  61 ! and Directed Research and Development program.
  62 !************************************************************************
  63
  64 !WRF:MODEL_LAYER:CHEMISTRY
  65 !
  66         module module_phot_fastj
  67         integer, parameter :: lunerr = -1
  68
  69         contains
  70 !***********************************************************************
  71        subroutine fastj_driver(id,curr_secs,dtstep,config_flags,       &
  72                gmt,julday,t_phy,moist,p8w,p_phy,                       &
  73                chem,rho_phy,dz8w,xlat,xlong,z_at_w,                    &
  74                ph_o2,ph_o31d,ph_o33p,ph_no2,ph_no3o2,ph_no3o,ph_hno2,  &
  75                ph_hno3,ph_hno4,ph_h2o2,ph_ch2or,ph_ch2om,ph_ch3cho,    &
  76                ph_ch3coch3,ph_ch3coc2h5,ph_hcocho,ph_ch3cocho,         &
  77                ph_hcochest,ph_ch3o2h,ph_ch3coo2h,ph_ch3ono2,ph_hcochob,&
  78                ph_n2o5,                                                &
  79                tauaer1,tauaer2,tauaer3,tauaer4,                        &
  80                gaer1,gaer2,gaer3,gaer4,                                &
  81                waer1,waer2,waer3,waer4,                                &
  82                bscoef1,bscoef2,bscoef3,bscoef4,                        &
  83                l2aer,l3aer,l4aer,l5aer,l6aer,l7aer,                    &
  84                ids,ide, jds,jde, kds,kde,                              &
  85                ims,ime, jms,jme, kms,kme,                              &
  86                its,ite, jts,jte, kts,kte                               )
  87
  88
  89    USE module_configure
  90    USE module_state_description
  91    USE module_data_mosaic_therm, only: nbin_a, nbin_a_maxd
  92 !  USE module_data_mosaic_asect
  93    USE module_data_mosaic_other, only: kmaxd, nsubareas
  94    USE module_fastj_mie
  95 !  USE module_mosaic_therm, only: aerosol_optical_properties
  96 !  USE module_mosaic_driver, only: mapaer_tofrom_host
  97    USE module_fastj_data, only:  nb, nc
  98
  99    implicit none
 100 !jdf
 101 !  integer, parameter :: iprint = 0
 102    integer, parameter :: single = 4        !compiler dependent value real*4
 103    integer, parameter :: double = 8        !compiler dependent value real*8
 104    integer,parameter :: ipar_fastj=1,jpar=1
 105    integer,parameter :: jppj=14        !Number of photolytic reactions supplied
 106    logical,parameter :: ldeg45=.false. !Logical flag for degraded CTM resolution
 107    integer,save :: lpar           !Number of levels in CTM
 108    integer,save :: jpnl           !Number of levels requiring chemistry
 109    real(kind=double), dimension(ipar_fastj) :: xgrd !  Longitude (midpoint, radians)
 110    real(kind=double), dimension(jpar) :: ygrd           !  Latitude  (midpoint, radians)
 111    real(kind=double), dimension(jpar) :: ydgrd          !  Latitude  (midpoint, degrees)
 112    real(kind=double), dimension(kmaxd+1) :: etaa    !  Eta(a) value for level boundaries
 113    real(kind=double), dimension(kmaxd+1) :: etab    !  Eta(b) value for level boundaries
 114    real(kind=double) ::  tau_fastj          !  Time of Day (hours, GMT)
 115    integer month_fastj        !  Number of month (1-12)
 116    integer iday_fastj         !  Day of year
 117    integer nspint           ! Num of spectral intervals across solar
 118    parameter ( nspint = 4 ) ! spectrum for FAST-J
 119    real, dimension (nspint),save :: wavmid !cm
 120    real, dimension (nspint, kmaxd+1),save :: sizeaer,extaer,waer,gaer,tauaer,bscoef
 121    real, dimension (nspint, kmaxd+1),save :: l2,l3,l4,l5,l6,l7
 122    data wavmid     &
 123        / 0.30e-4, 0.40e-4, 0.60e-4 ,0.999e-04 /
 124    integer nphoto_fastj
 125    parameter (nphoto_fastj = 14)
 126    integer   &
 127    lfastj_no2,   lfastj_o3a,   lfastj_o3b,    lfastj_h2o2,   &
 128    lfastj_hchoa, lfastj_hchob, lfastj_ch3ooh, lfastj_no3x,   &
 129    lfastj_no3l,  lfastj_hono,  lfastj_n2o5,   lfastj_hno3,   &
 130    lfastj_hno4
 131    parameter( lfastj_no2   = 1 )
 132    parameter( lfastj_o3a   = 2 )
 133    parameter( lfastj_o3b   = 3 )
 134    parameter( lfastj_h2o2  = 4 )
 135    parameter( lfastj_hchoa = 5 )
 136    parameter( lfastj_hchob = 6 )
 137    parameter( lfastj_ch3ooh= 7 )
 138    parameter( lfastj_no3x  = 8 )
 139    parameter( lfastj_no3l  = 9 )
 140    parameter( lfastj_hono  = 10 )
 141    parameter( lfastj_n2o5  = 11 )
 142    parameter( lfastj_hno3  = 12 )
 143    parameter( lfastj_hno4  = 13 )
 144 !jdf
 145    INTEGER,      INTENT(IN   ) :: id,julday,                    &
 146                                   ids,ide, jds,jde, kds,kde,    &
 147                                   ims,ime, jms,jme, kms,kme,    &
 148                                   its,ite, jts,jte, kts,kte
 149    REAL(KIND=8), INTENT(IN   ) :: curr_secs
 150    REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_moist ),        &
 151          INTENT(IN ) ::                                   moist
 152    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ),                    &
 153          INTENT(INOUT ) ::                                          &
 154            ph_o2,ph_o31d,ph_o33p,ph_no2,ph_no3o2,ph_no3o,ph_hno2,   &
 155            ph_hno3,ph_hno4,ph_h2o2,ph_ch2or,ph_ch2om,ph_ch3cho,     &
 156            ph_ch3coch3,ph_ch3coc2h5,ph_hcocho,ph_ch3cocho,          &
 157            ph_hcochest,ph_ch3o2h,ph_ch3coo2h,ph_ch3ono2,ph_hcochob, &
 158            ph_n2o5
 159    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ),                    &
 160          INTENT(IN ) ::                                             &
 161            tauaer1,tauaer2,tauaer3,tauaer4,                         &
 162            gaer1,gaer2,gaer3,gaer4,                                 &
 163            waer1,waer2,waer3,waer4,                                 &
 164            bscoef1,bscoef2,bscoef3,bscoef4
 165    REAL, DIMENSION( ims:ime, kms:kme, jms:jme, 1:4 ),               &
 166          INTENT(IN ) ::                                             &
 167            l2aer,l3aer,l4aer,l5aer,l6aer,l7aer
 168
 169    REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ),         &
 170          INTENT(INOUT ) ::                                chem
 171    REAL,  DIMENSION( ims:ime , kms:kme , jms:jme )         ,    &
 172           INTENT(IN   ) ::                                      &
 173                                                       t_phy,    &
 174                                                       p_phy,    &
 175                                                        dz8w,    &
 176                                                         p8w,    &
 177                                                     rho_phy,    &
 178                                                       z_at_w
 179    REAL,  DIMENSION( ims:ime , jms:jme )                   ,    &
 180           INTENT(IN   ) ::                             xlat,    &
 181                                                       xlong
 182    REAL,      INTENT(IN   ) ::                                  &
 183                                                  dtstep,gmt
 184
 185    TYPE(grid_config_rec_type), INTENT(IN  ) :: config_flags
 186
 187
 188 !local variables
 189
 190       integer iclm, jclm
 191
 192       real number_bin(nbin_a_maxd,kmaxd) !These have to be the full kmaxd to
 193       real radius_wet(nbin_a_maxd,kmaxd) !match arrays in MOSAIC routines.
 194       complex refindx(nbin_a_maxd,kmaxd)
 195
 196       integer(kind=8) :: ixhour
 197       integer i,j, k, nsub
 198           real(kind=8) :: xtime, xhour
 199       real xmin, gmtp, tmidh
 200       real sla, slo
 201       real psfc
 202
 203       real cos_sza
 204       real, dimension(kts:kte) :: temp, ozone, dz
 205       real, dimension(0:kte) :: pbnd
 206       real, dimension(kts:kte) :: cloudmr, airdensity, relhum
 207       real, dimension(kts:kte+1) :: zatw
 208
 209       real valuej(kte,nphoto_fastj)
 210
 211       logical processingAerosols
 212
 213 !   set "pegasus" grid size variables
 214 !       itot = ite
 215 !       jtot = jte
 216     lpar = kte
 217     jpnl = kte
 218     nb   = lpar + 1 !for module module_fastj_data
 219     nc   = 2*nb     !ditto, and don't confuse this with nc in module_fastj_mie
 220         nsubareas = 1
 221         if ((kte > kmaxd) .or. (lpar <= 0)) then
 222             write( wrf_err_message, '(a,4i5)' )   &
 223                 '*** subr fastj_driver -- ' //   &
 224                 'lpar, kmaxd, kts, kte', lpar, kmaxd, kts, kte
 225         call wrf_message( trim(wrf_err_message) )
 226             wrf_err_message = '*** subr fastj_driver -- ' //   &
 227              'kte>kmaxd OR lpar<=0'
 228             call wrf_error_fatal( wrf_err_message )
 229         end if
 230
 231 ! Determine if aerosol data is provided in the chem array. Currently,
 232 ! only MOSAIC will work. The Mie routine does not know how to handle
 233 ! SORGAM aerosols.
 234     select case (config_flags%chem_opt)
 235     case ( CBMZ_MOSAIC_4BIN,    CBMZ_MOSAIC_8BIN, CBMZ_MOSAIC_KPP,  &
 236            CBMZ_MOSAIC_4BIN_AQ, CBMZ_MOSAIC_8BIN_AQ, &
 237            CBMZ_MOSAIC_DMS_4BIN, CBMZ_MOSAIC_DMS_8BIN, &
 238            CBMZ_MOSAIC_DMS_4BIN_AQ, CBMZ_MOSAIC_DMS_8BIN_AQ, &
 239            MOZART_MOSAIC_4BIN_KPP, MOZART_MOSAIC_4BIN_AQ_KPP, &
 240            CRI_MOSAIC_8BIN_AQ_KPP, CRI_MOSAIC_4BIN_AQ_KPP, SAPRC99_MOSAIC_8BIN_VBS2_AQ_KPP, &!BSINGH(12/05/2013): Added for SAPRC 8 bin vbs
 241            SAPRC99_MOSAIC_8BIN_VBS2_KPP )!BSINGH(04/07/2014): Added for SAPRC 8 bin vbs non-aq
 242        processingAerosols = .true.
 243     case default
 244        processingAerosols = .false.
 245     end select
 246
 247 ! Set nbin_a = ntot_amode and check nbin_a <= nbin_a_maxd.
 248 ! This duplicates the same assignment and check in module_mosaic_therm.F,
 249 ! but photolysis is called before aeorosols so this must be set too.
 250 !
 251 ! rce 2004-dec-07 -- nbin_a is initialized elsewhere
 252 !!$      nbin_a = ntot_amode
 253 !!$      if ((nbin_a .gt. nbin_a_maxd) .or. (nbin_a .le. 0)) then
 254 !!$          write( wrf_err_message, '(a,2(1x,i4))' )   &
 255 !!$              '*** subr fastj_driver -- nbin_a, nbin_a_maxd =',   &
 256 !!$              nbin_a, nbin_a_maxd
 257 !!$          call wrf_message( wrf_err_message )
 258 !!$          call wrf_error_fatal(   &
 259 !!$              '*** subr fastj_driver -- BAD VALUE for nbin_a' )
 260 !!$      end if
 261
 262 ! determine current time of day in Greenwich Mean Time at middle
 263 ! of current time step, tmidh.  do this by computing the number of minutes
 264 ! from beginning of run to middle of current time step
 265     xtime    = curr_secs/60._8 + real(dtstep,8)/120._8
 266     ixhour   = int(gmt + 0.01,8) + int(xtime/60._8,8)
 267     xhour    = real(ixhour,8)   !current hour
 268     xmin     = 60.*gmt + real(xtime-xhour*60_8,8)
 269     gmtp     = mod(xhour,24._8)
 270     tmidh    = gmtp + xmin/60.
 271
 272 ! execute for each i,j column and each nsub subarea
 273     do nsub = 1, nsubareas
 274         do jclm = jts, jte
 275         do iclm = its, ite
 276
 277        do k = kts, lpar
 278           dz(k) = dz8w(iclm, k, jclm)   ! cell depth (m)
 279        end do
 280
 281        if( processingAerosols ) then
 282          do k = kts, lpar
 283            l2(1,k)=l2aer(iclm,k,jclm,1)
 284            l2(2,k)=l2aer(iclm,k,jclm,2)
 285            l2(3,k)=l2aer(iclm,k,jclm,3)
 286            l2(4,k)=l2aer(iclm,k,jclm,4)
 287            l3(1,k)=l3aer(iclm,k,jclm,1)
 288            l3(2,k)=l3aer(iclm,k,jclm,2)
 289            l3(3,k)=l3aer(iclm,k,jclm,3)
 290            l3(4,k)=l3aer(iclm,k,jclm,4)
 291            l4(1,k)=l4aer(iclm,k,jclm,1)
 292            l4(2,k)=l4aer(iclm,k,jclm,2)
 293            l4(3,k)=l4aer(iclm,k,jclm,3)
 294            l4(4,k)=l4aer(iclm,k,jclm,4)
 295            l5(1,k)=l5aer(iclm,k,jclm,1)
 296            l5(2,k)=l5aer(iclm,k,jclm,2)
 297            l5(3,k)=l5aer(iclm,k,jclm,3)
 298            l5(4,k)=l5aer(iclm,k,jclm,4)
 299            l6(1,k)=l6aer(iclm,k,jclm,1)
 300            l6(2,k)=l6aer(iclm,k,jclm,2)
 301            l6(3,k)=l6aer(iclm,k,jclm,3)
 302            l6(4,k)=l6aer(iclm,k,jclm,4)
 303            l7(1,k)=l7aer(iclm,k,jclm,1)
 304            l7(2,k)=l7aer(iclm,k,jclm,2)
 305            l7(3,k)=l7aer(iclm,k,jclm,3)
 306            l7(4,k)=l7aer(iclm,k,jclm,4)
 307          enddo
 308 ! take chem data and extract 1 column to create rsub(l,k,m) array
 309 !         call mapaer_tofrom_host( 0,                     &
 310 !              ims,ime, jms,jme, kms,kme,                    &
 311 !              its,ite, jts,jte, kts,kte,                    &
 312 !              iclm, jclm, kts,lpar,                         &
 313 !              num_moist, num_chem, moist, chem,             &
 314 !              t_phy, p_phy, rho_phy                         )
 315 ! generate aerosol optical properties for cells in this column
 316 ! subroutine is located in file module_mosaic_therm
 317 !         call aerosol_optical_properties(iclm, jclm, lpar, refindx, &
 318 !              radius_wet, number_bin)
 319 ! execute mie code , located in file module_fastj_mie
 320 !         CALL wrf_debug(250,'fastj_driver: calling mieaer')
 321 !         call mieaer(                            &
 322 !              id, iclm, jclm, nbin_a,            &
 323 !              number_bin, radius_wet, refindx,   &
 324 !              dz, curr_secs, lpar,               &
 325 !              sizeaer,extaer,waer,gaer,tauaer,l2,l3,l4,l5,l6,l7,bscoef)
 326        end if
 327
 328 ! set lat, long
 329           sla = xlat(iclm,jclm)
 330           slo = xlong(iclm,jclm)
 331 ! set column pressures, temperature, and ozone
 332           psfc = p8w(iclm,1,jclm) * 10. ! convert pascals to dynes/cm2
 333           do k = kts, lpar
 334             pbnd(k) = p8w(iclm,k+1,jclm) *10.  ! convert pascals to dynes/cm2
 335             temp(k) = t_phy(iclm,k,jclm)
 336             ozone(k) = chem(iclm,k,jclm,p_o3) / 1.0e6   ! ppm->mol/mol air
 337         cloudmr(k) = moist(iclm,k,jclm,p_qc)/0.622
 338         airdensity(k) = rho_phy(iclm,k,jclm)/28.966e3
 339         relhum(k) = MIN( .95, moist(iclm,k,jclm,p_qv) / &
 340                        (3.80*exp(17.27*(t_phy(iclm,k,jclm)-273.)/ &
 341                        (t_phy(iclm,k,jclm)-36.))/(.01*p_phy(iclm,k,jclm))))
 342         relhum(k) = MAX(.001,relhum(k))
 343         zatw(k)=z_at_w(iclm,k,jclm)
 344           end do
 345       zatw(lpar+1)=z_at_w(iclm,lpar+1,jclm)
 346 ! call interface_fastj
 347           CALL wrf_debug(250,'fastj_driver: calling interface_fastj')
 348           call interface_fastj(tmidh,sla,slo,julday,           &
 349            pbnd, psfc, temp, ozone,                        &
 350            dz, cloudmr, airdensity, relhum, zatw,          &
 351            iclm, jclm, lpar, jpnl,                         &
 352            curr_secs, valuej, cos_sza, processingAerosols, &
 353            sizeaer,extaer,waer,gaer,tauaer,l2,l3,l4,l5,l6,l7)
 354 ! put column photolysis rates (valuej) into wrf photolysis (i,k,j) arrays
 355       CALL wrf_debug(250,'fastj_driver: calling mapJrates_tofrom_host')
 356       call mapJrates_tofrom_host( 0,                                &
 357            ims,ime, jms,jme, kms,kme,                                       &
 358            its,ite, jts,jte, kts,kte,                                       &
 359            iclm, jclm, kts,lpar,                                            &
 360            valuej,                                                                      &
 361            ph_o2,ph_o31d,ph_o33p,ph_no2,ph_no3o2,ph_no3o,ph_hno2,   &
 362            ph_hno3,ph_hno4,ph_h2o2,ph_ch2or,ph_ch2om,ph_ch3cho,         &
 363            ph_ch3coch3,ph_ch3coc2h5,ph_hcocho,ph_ch3cocho,              &
 364            ph_hcochest,ph_ch3o2h,ph_ch3coo2h,ph_ch3ono2,ph_hcochob, &
 365            ph_n2o5                                                  )
 366 ! put the aerosol optical properties into the wrf arrays (this is hard-
 367 ! coded to 4 spectral bins, nspint)
 368 !     do k=kts,kte
 369 !        tauaer1(iclm,k,jclm) = tauaer(1,k)
 370 !        tauaer2(iclm,k,jclm) = tauaer(2,k)
 371 !        tauaer3(iclm,k,jclm) = tauaer(3,k)
 372 !        tauaer4(iclm,k,jclm) = tauaer(4,k)
 373 !        gaer1(iclm,k,jclm)   = gaer(1,k)
 374 !        gaer2(iclm,k,jclm)   = gaer(2,k)
 375 !        gaer3(iclm,k,jclm)   = gaer(3,k)
 376 !        gaer4(iclm,k,jclm)   = gaer(4,k)
 377 !        waer1(iclm,k,jclm)   = waer(1,k)
 378 !        waer2(iclm,k,jclm)   = waer(2,k)
 379 !        waer3(iclm,k,jclm)   = waer(3,k)
 380 !        waer4(iclm,k,jclm)   = waer(4,k)
 381 !        bscoef1(iclm,k,jclm) = bscoef(1,k)
 382 !        bscoef2(iclm,k,jclm) = bscoef(2,k)
 383 !        bscoef3(iclm,k,jclm) = bscoef(3,k)
 384 !        bscoef4(iclm,k,jclm) = bscoef(4,k)
 385 !     end do
 386
 387     end do
 388     end do
 389     end do
 390
 391     return
 392     end subroutine  fastj_driver
 393
 394 !-----------------------------------------------------------------------
 395         subroutine mapJrates_tofrom_host( iflag,                        &
 396                 ims,ime, jms,jme, kms,kme,                                      &
 397                 its,ite, jts,jte, kts,kte,                                      &
 398         iclm, jclm, ktmaps,ktmape,                                          &
 399         valuej,                                                                             &
 400         ph_o2,ph_o31d,ph_o33p,ph_no2,ph_no3o2,ph_no3o,ph_hno2,      &
 401         ph_hno3,ph_hno4,ph_h2o2,ph_ch2or,ph_ch2om,ph_ch3cho,        &
 402         ph_ch3coch3,ph_ch3coc2h5,ph_hcocho,ph_ch3cocho,             &
 403         ph_hcochest,ph_ch3o2h,ph_ch3coo2h,ph_ch3ono2,ph_hcochob,    &
 404         ph_n2o5                                                     )
 405
 406         USE module_data_cbmz
 407
 408    implicit none
 409 !jdf
 410    integer nphoto_fastj
 411    parameter (nphoto_fastj = 14)
 412    integer   &
 413    lfastj_no2,   lfastj_o3a,   lfastj_o3b,    lfastj_h2o2,   &
 414    lfastj_hchoa, lfastj_hchob, lfastj_ch3ooh, lfastj_no3x,   &
 415    lfastj_no3l,  lfastj_hono,  lfastj_n2o5,   lfastj_hno3,   &
 416    lfastj_hno4
 417    parameter( lfastj_no2   = 1 )
 418    parameter( lfastj_o3a   = 2 )
 419    parameter( lfastj_o3b   = 3 )
 420    parameter( lfastj_h2o2  = 4 )
 421    parameter( lfastj_hchoa = 5 )
 422    parameter( lfastj_hchob = 6 )
 423    parameter( lfastj_ch3ooh= 7 )
 424    parameter( lfastj_no3x  = 8 )
 425    parameter( lfastj_no3l  = 9 )
 426    parameter( lfastj_hono  = 10 )
 427    parameter( lfastj_n2o5  = 11 )
 428    parameter( lfastj_hno3  = 12 )
 429    parameter( lfastj_hno4  = 13 )
 430 !jdf
 431    INTEGER,      INTENT(IN   ) :: iflag,                        &
 432                                   ims,ime, jms,jme, kms,kme,    &
 433                                   its,ite, jts,jte, kts,kte,    &
 434                                   iclm, jclm, ktmaps, ktmape
 435    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ),                   &
 436          INTENT(INOUT ) ::                                         &
 437            ph_o2,ph_o31d,ph_o33p,ph_no2,ph_no3o2,ph_no3o,ph_hno2,  &
 438            ph_hno3,ph_hno4,ph_h2o2,ph_ch2or,ph_ch2om,ph_ch3cho,    &
 439            ph_ch3coch3,ph_ch3coc2h5,ph_hcocho,ph_ch3cocho,         &
 440            ph_hcochest,ph_ch3o2h,ph_ch3coo2h,ph_ch3ono2,ph_hcochob,&
 441            ph_n2o5
 442
 443    REAL, DIMENSION( kte,nphoto_fastj ), INTENT(INOUT) :: valuej
 444
 445 ! local variables
 446         real ft
 447         integer kt
 448
 449         ft = 60.
 450
 451     if (iflag .gt. 0) go to 2000
 452 ! flag is <=0, put pegasus column J rates (in 1/sec) into WRF arrays (in 1/min)
 453         do kt = ktmaps, ktmape
 454           ph_no2(iclm,kt,jclm)       = valuej(kt,lfastj_no2) * ft
 455           ph_no3o(iclm,kt,jclm)      = valuej(kt,lfastj_no3x) * ft
 456           ph_no3o2(iclm,kt,jclm)     = valuej(kt,lfastj_no3l) * ft
 457           ph_o33p(iclm,kt,jclm)      = valuej(kt,lfastj_o3a) * ft
 458           ph_o31d(iclm,kt,jclm)      = valuej(kt,lfastj_o3b) * ft
 459           ph_hno2(iclm,kt,jclm)      = valuej(kt,lfastj_hono) * ft
 460           ph_hno3(iclm,kt,jclm)      = valuej(kt,lfastj_hno3) * ft
 461           ph_hno4(iclm,kt,jclm)      = valuej(kt,lfastj_hno4) * ft
 462           ph_h2o2(iclm,kt,jclm)      = valuej(kt,lfastj_h2o2) * ft
 463           ph_ch3o2h(iclm,kt,jclm)    = valuej(kt,lfastj_ch3ooh) * ft
 464           ph_ch2or(iclm,kt,jclm)     = valuej(kt,lfastj_hchoa) * ft
 465           ph_ch2om(iclm,kt,jclm)     = valuej(kt,lfastj_hchob) * ft
 466           ph_n2o5(iclm,kt,jclm)      = valuej(kt,lfastj_n2o5) * ft
 467
 468           ph_o2(iclm,kt,jclm)        = 0.0
 469           ph_ch3cho(iclm,kt,jclm)    = 0.0
 470           ph_ch3coch3(iclm,kt,jclm)  = 0.0
 471           ph_ch3coc2h5(iclm,kt,jclm) = 0.0
 472           ph_hcocho(iclm,kt,jclm)    = 0.0
 473           ph_ch3cocho(iclm,kt,jclm)  = 0.0
 474           ph_hcochest(iclm,kt,jclm)  = 0.0
 475           ph_ch3coo2h(iclm,kt,jclm)  = 0.0
 476           ph_ch3ono2(iclm,kt,jclm)   = 0.0
 477           ph_hcochob(iclm,kt,jclm)   = 0.0
 478
 479         end do
 480         return  !finished peg-> wrf mapping
 481
 482 2000    continue
 483 ! iflag > 0 ; put wrf ph_xxx Jrates (1/min) into pegasus column valuej (1/sec)
 484         do kt = ktmaps, ktmape
 485           valuej(kt,lfastj_no2)    = ph_no2(iclm,kt,jclm) /  ft
 486           valuej(kt,lfastj_no3x)   = ph_no3o(iclm,kt,jclm) / ft
 487           valuej(kt,lfastj_no3l)   = ph_no3o2(iclm,kt,jclm)/ ft
 488           valuej(kt,lfastj_o3a)    = ph_o33p(iclm,kt,jclm) / ft
 489           valuej(kt,lfastj_o3b)    = ph_o31d(iclm,kt,jclm) / ft
 490           valuej(kt,lfastj_hono)   = ph_hno2(iclm,kt,jclm) / ft
 491           valuej(kt,lfastj_hno3)   = ph_hno3(iclm,kt,jclm) / ft
 492           valuej(kt,lfastj_hno4)   = ph_hno4(iclm,kt,jclm) / ft
 493           valuej(kt,lfastj_h2o2)   = ph_h2o2(iclm,kt,jclm) / ft
 494           valuej(kt,lfastj_ch3ooh) = ph_ch3o2h(iclm,kt,jclm)/ft
 495           valuej(kt,lfastj_hchoa)  = ph_ch2or(iclm,kt,jclm)/ ft
 496           valuej(kt,lfastj_hchob)  = ph_ch2om(iclm,kt,jclm)/ ft
 497           valuej(kt,lfastj_n2o5)   = ph_n2o5(iclm,kt,jclm) / ft
 498         end do
 499
 500     return      !finished wrf->peg mapping
 501
 502     end subroutine mapJrates_tofrom_host
 503 !-----------------------------------------------------------------------
 504
 505
 506
 507 !***********************************************************************
 508         subroutine interface_fastj(tmidh,sla,slo,julian_day,   &
 509              pbnd, psfc, temp, ozone,                          &
 510              dz, cloudmr, airdensity, relhum, zatw,            &
 511              isvode, jsvode, lpar, jpnl,                       &
 512              curr_secs, valuej, cos_sza, processingAerosols,   &
 513              sizeaer,extaer,waer,gaer,tauaer,l2,l3,l4,l5,l6,l7)
 514 !-----------------------------------------------------------------
 515 ! sets parameters for fastj call.
 516 ! inputs
 517 !       tmidh -- GMT time in decimal hours at which to calculate
 518 !                photolysis rates
 519 !       sla -- latitude, decimal degrees in real*8
 520 !       slo -- negative of the longitude, decimal degrees in real*8
 521 !       julian_day -- day of the year in julian days
 522 !    pbnd(0:lpar) = pressure at top boundary of cell k (dynes/cm^2).
 523 !    psfc = surface pressure (dynes/cm^2).
 524 !    temp(lpar)= mid-cell temperature values (deg K)
 525 !    ozone(lpar) = mid-cell ozone mixing ratios
 526 !    surface_albedo -- broadband albedo (dimensionless)
 527 !    curr_secs -- elapsed time from start of simulation in seconds.
 528 !    isvode,jsvode  -- current column i,j.
 529 !
 530 !    lpar -- vertical extent of column (from module_fastj_cmnh)
 531 !
 532 ! outputs
 533 !       cos_sza -- cosine of solar zenith angle.
 534 !       valuej(lpar,nphoto_fastj-1) -- array of photolysis rates, s-1.
 535 !
 536 !
 537 ! local variables
 538 !    surface_pressure_mb -- surface pressure (mb). equal to col_press_mb(1).
 539 !    col_press_mb(lpar) -- for the column, grid cell boundary pressures
 540 !          (not at cell centers) up until the bottom pressure for the
 541 !          top cell (mb).
 542 !    col_temp_K(lpar+1) -- for the column, grid cell center temperature (deg K)
 543 !    col_ozone(lpar+1) -- for the column, grid cell center ozone mixing
 544 !          ratios (dimensionless)
 545 !    col_optical_depth(lpar+1) -- for the column, grid cell center cloud
 546 !          optical depths (dimensionless).SET TO ZERO IN THIS VERSION
 547 !    tauaer_550 -- aerosol optical thickness at 550 nm.
 548 !       note:  photolysis rates are calculated at centers of model layers
 549 !       the pressures are given at the boundaries defining
 550 !       the top and bottom of the layers
 551 !       so the number of pressure values is equal
 552 !       to the (number of layers) + 1 ; the last pressure is set = 0 in fastj code.
 553 !       pressures from the surface up to the bottom of the top (lpar<=kmaxd) cell
 554 !       ******** pressure 2
 555 !       layer 1 - temperature,optical depth, and O3 given here
 556 !       ******** pressure 1
 557 ! the optical depth is appropriate for the layer depth
 558 ! conversion factor:   1 dyne/cm2 = 0.001 mb
 559 !-----------------------------------------------------------------
 560         USE module_data_mosaic_other, only : kmaxd
 561         USE module_peg_util, only:  peg_message, peg_error_fatal
 562
 563         IMPLICIT NONE
 564
 565 !jdf
 566    integer, parameter :: iprint = 0
 567    integer, parameter :: single = 4        !compiler dependent value real*4
 568    integer, parameter :: double = 8        !compiler dependent value real*8
 569    integer,parameter :: ipar_fastj=1,jpar=1
 570    integer,parameter :: jppj=14        !Number of photolytic reactions supplied
 571    logical,parameter :: ldeg45=.false. !Logical flag for degraded CTM resolution
 572    integer lpar           !Number of levels in CTM
 573    integer jpnl           !Number of levels requiring chemistry
 574    real(kind=double), dimension(ipar_fastj) :: xgrd !  Longitude (midpoint, radians)
 575    real(kind=double), dimension(jpar) :: ygrd           !  Latitude  (midpoint, radians)
 576    real(kind=double), dimension(jpar) :: ydgrd          !  Latitude  (midpoint, degrees)
 577    real(kind=double), dimension(kmaxd+1) :: etaa    !  Eta(a) value for level boundaries
 578    real(kind=double), dimension(kmaxd+1) :: etab    !  Eta(b) value for level boundaries
 579    real(kind=double) ::  tau_fastj          !  Time of Day (hours, GMT)
 580    integer month_fastj        !  Number of month (1-12)
 581    integer iday_fastj         !  Day of year
 582    integer nphoto_fastj
 583    parameter (nphoto_fastj = 14)
 584    integer   &
 585    lfastj_no2,   lfastj_o3a,   lfastj_o3b,    lfastj_h2o2,   &
 586    lfastj_hchoa, lfastj_hchob, lfastj_ch3ooh, lfastj_no3x,   &
 587    lfastj_no3l,  lfastj_hono,  lfastj_n2o5,   lfastj_hno3,   &
 588    lfastj_hno4
 589    parameter( lfastj_no2   = 1 )
 590    parameter( lfastj_o3a   = 2 )
 591    parameter( lfastj_o3b   = 3 )
 592    parameter( lfastj_h2o2  = 4 )
 593    parameter( lfastj_hchoa = 5 )
 594    parameter( lfastj_hchob = 6 )
 595    parameter( lfastj_ch3ooh= 7 )
 596    parameter( lfastj_no3x  = 8 )
 597    parameter( lfastj_no3l  = 9 )
 598    parameter( lfastj_hono  = 10 )
 599    parameter( lfastj_n2o5  = 11 )
 600    parameter( lfastj_hno3  = 12 )
 601    parameter( lfastj_hno4  = 13 )
 602    integer nspint           ! Num of spectral intervals across solar
 603    parameter ( nspint = 4 ) ! spectrum for FAST-J
 604    real, dimension (nspint),save :: wavmid !cm
 605    real, dimension (nspint, kmaxd+1) :: sizeaer,extaer,waer,gaer,tauaer
 606    real, dimension (nspint, kmaxd+1) :: l2,l3,l4,l5,l6,l7
 607    data wavmid     &
 608        / 0.30e-4, 0.40e-4, 0.60e-4 ,0.999e-04 /
 609 !jdf
 610         real pbnd(0:lpar), psfc
 611         real temp(lpar), ozone(lpar), surface_albedo
 612         real dz(lpar), cloudmr(lpar), airdensity(lpar), relhum(lpar), zatw(lpar+1)
 613     real(kind=8) :: curr_secs
 614         integer isvode, jsvode
 615
 616         real cos_sza
 617     integer,parameter :: lunout=41
 618
 619         real valuej(lpar,nphoto_fastj)
 620
 621     real hl,rhl,factor1,part1,part2,cfrac,rhfrac
 622     real emziohl(lpar+1),clwp(lpar)
 623 !ec material to check output
 624         real valuej_no3rate(lpar)
 625
 626         real*8 lat,lon
 627     real*8 jvalue(lpar,nphoto_fastj)
 628     real sza
 629         real tau1
 630         real tmidh, sla, slo
 631
 632         integer julian_day,iozone1
 633     integer,parameter :: nfastj_rxns = 14
 634         integer k, l
 635
 636         real surface_pressure_mb, tauaer_550,   &
 637          col_press_mb,col_temp_K,col_ozone,col_optical_depth
 638         dimension col_press_mb(lpar+2),col_temp_K(lpar+1),   &
 639                 col_ozone(lpar+1),col_optical_depth(lpar+1)
 640     character*80 msg
 641
 642 ! define logical processingAerosols
 643 ! if processingAerosols = true, uses values calculated in subroutine
 644 ! mieaer for variables & arrays in common block mie.
 645 ! if processingAerosols = false, sets all variables & arrays in common
 646 ! block mie to zero.  (JCB-revised Fast-J requires common block mie info,
 647 ! regardless of whether aerosols are present or not.  Original Wild Fast-J
 648 ! did not use common block mie info.)
 649
 650       logical processingAerosols
 651
 652 ! set lat and longitude as real*8 for consistency with fastj code.
 653 ! variables lat and lon previously declared as reals
 654         lat  = sla
 655         lon  = slo
 656 !
 657 ! cloud optical depths currently treated by using fractional cloudiness
 658 ! based on relative humidity. cloudmr set up to use cloud liquid water
 659 ! but hooks into microphysics need to be tested - for now set cloudmr=0
 660 !
 661 ! parameters to calculate 'typical' liquid cloudwater path values for
 662 ! non convective clouds based on approximations in NCAR's CCM2
 663 ! 0.18 = reference liquid water concentration (gh2o/m3)
 664 ! hl = liquid water scale height (m)
 665 !
 666     hl=1080.+2000.0*cos(lat*0.017454329)
 667     rhl=1.0/hl
 668         do k =1, lpar+1
 669        emziohl(k)=exp(-zatw(k)*rhl)
 670     enddo
 671         do k =1, lpar
 672        clwp(k)=0.18*hl*(emziohl(k)-emziohl(k+1))
 673     enddo
 674 ! assume radius of cloud droplets is constant at 10 microns (0.001 cm) and
 675 ! that density of water is constant at 1 g2ho/cm3
 676 !       factor1=3./2./0.001/1.
 677     factor1=1500.0
 678         do k =1, lpar
 679            col_optical_depth(k) = 0.0
 680        cfrac=0.0
 681        cloudmr(k)=0.0
 682        if(cloudmr(k).gt.0.0) cfrac=1.0
 683 ! 18.0*airdensity converts mole h2o/mole air to g h2o/cm3, part1 is in g h2o/cm2
 684        part1=cloudmr(k)*cfrac*18.0*airdensity(k)*dz(k)*100.0
 685        if(relhum(k).lt.0.8) then
 686           rhfrac=0.0
 687        elseif(relhum(k).le.1.0.and.relhum(k).ge.0.8) then
 688 !            rhfrac=(relhum(k)-0.8)/(1.0-0.8)
 689           rhfrac=(relhum(k)-0.8)/0.2
 690        else
 691           rhfrac=1.0
 692        endif
 693        if(rhfrac.ge.0.01) then
 694 ! factor 1.0e4 converts clwp of g h2o/m2 to g h2o/cm2
 695           part2=rhfrac*clwp(k)/1.e4
 696        else
 697           part2=0.0
 698        endif
 699        if(cfrac.gt.0) part2=0.0
 700        col_optical_depth(k) = factor1*(part1+part2)
 701 !          col_optical_depth(k) = 0.0
 702 !       if(isvode.eq.33.and.jsvode.eq.34) &
 703 !          print*,'jdf opt',isvode,jsvode,k,col_optical_depth(k), &
 704 !          cfrac,rhfrac,relhum(k),clwp(k)
 705         end do
 706     col_optical_depth(lpar+1) = 0.0
 707         if (.not.processingAerosols) then
 708 ! set common block mie variables to 0 if
 709            call set_common_mie(lpar, &
 710            sizeaer,extaer,waer,gaer,tauaer,l2,l3,l4,l5,l6,l7)
 711         end if                           ! processingAerosols=false
 712
 713 ! set pressure, temperature, ozone of each cell in the column
 714 ! set iozone1 = lpar to allow replacement of climatological ozone with model
 715 ! predicted ozone to top of chemistry column; standard fastj climatological o3
 716 ! thereafter.
 717         surface_pressure_mb = psfc * 0.001
 718         tau1 = tmidh
 719         col_press_mb(1) = psfc * 0.001
 720         iozone1 = lpar
 721         do k =1, lpar
 722        col_press_mb(k+1) = pbnd(k) * 0.001
 723            col_temp_K(k) = temp(k)
 724            col_ozone(k) = ozone(k)
 725         end do
 726
 727         surface_albedo=0.055
 728
 729 ! set aerosol parameters needed by Fast-J
 730         if (processingAerosols) then
 731             tauaer_550 = 0.0    ! needed parameters already calculated by subroutine
 732                                 ! mieaer and passed into proper parts of fastj code
 733                                 ! via module_fastj_cmnmie
 734          else
 735             tauaer_550 = 0.05   ! no aerosols, assume typical constant aerosol optical thickness
 736          end if
 737
 738           CALL wrf_debug(250,'interface_fastj: calling fastj')
 739       call fastj(isvode,jsvode,lat,lon,surface_pressure_mb,surface_albedo,   &
 740            julian_day,  tau1,   &
 741           col_press_mb, col_temp_K, col_optical_depth, col_ozone,   &
 742           iozone1,tauaer_550,jvalue,sza,lpar,jpnl, &
 743           sizeaer,extaer,waer,gaer,tauaer,l2,l3,l4,l5,l6,l7)
 744
 745         !cos_sza = cosd(sza)
 746         cos_sza = cos(sza*3.141592653/180.)
 747
 748
 749 ! array jvalue (real*8) is returned from fastj.  array valuej(unspecified
 750 ! real, default of real*4) is sent on to
 751 ! other chemistry subroutines
 752            do k = 1, lpar
 753              valuej(k, lfastj_no2) = jvalue(k,lfastj_no2)
 754              valuej(k, lfastj_o3a) = jvalue(k,lfastj_o3a)
 755              valuej(k, lfastj_o3b) = jvalue(k,lfastj_o3b)
 756              valuej(k, lfastj_h2o2) = jvalue(k,lfastj_h2o2)
 757              valuej(k, lfastj_hchoa) = jvalue(k,lfastj_hchoa)
 758              valuej(k, lfastj_hchob) = jvalue(k,lfastj_hchob)
 759              valuej(k, lfastj_ch3ooh) = jvalue(k,lfastj_ch3ooh)
 760              valuej(k, lfastj_no3x) = jvalue(k,lfastj_no3x)
 761              valuej(k, lfastj_no3l) = jvalue(k,lfastj_no3l)
 762              valuej(k, lfastj_hono) = jvalue(k,lfastj_hono)
 763              valuej(k, lfastj_n2o5) = jvalue(k,lfastj_n2o5)
 764              valuej(k, lfastj_hno3) = jvalue(k,lfastj_hno3)
 765              valuej(k, lfastj_hno4) = jvalue(k,lfastj_hno4)
 766           end do
 767 ! diagnostic output and zeroed value if negative photolysis rates returned
 768           do k = 1, lpar
 769          valuej(k,nphoto_fastj)=0.0
 770          do l = 1, nphoto_fastj-1
 771             if (valuej(k,l) .lt. 0) then
 772                write( msg, '(a,f14.2,4i4,1x,e11.4)' )   &
 773                     'FASTJ negative Jrate ' //   &
 774                     'tsec i j k l J(k,l)', curr_secs,isvode,jsvode,k,l,valuej(k,l)
 775                call peg_message( lunerr, msg )
 776                valuej(k,l) = 0.0
 777 ! following code used if want run stopped with negative Jrate
 778 !                 msg = '*** subr interface_fastj -- ' //   &
 779 !                   'Negative J rate returned from Fast-J'
 780 !                 call peg_error_fatal( lunerr, msg )
 781             end if
 782          end do
 783       end do
 784 ! compute overall no3 photolysis rate
 785 ! wig: commented out since it is not used anywhere
 786 !       do k = 1, lpar
 787 !          valuej_no3rate(k) =   &
 788 !                valuej(k, lfastj_no3x) + valuej(k,lfastj_no3l)
 789 !       end do
 790
 791         end subroutine interface_fastj
 792
 793 !***********************************************************************
 794         subroutine set_common_mie(lpar, &
 795           sizeaer,extaer,waer,gaer,tauaer,l2,l3,l4,l5,l6,l7)
 796 ! for use when aerosols are not included in a model run.  sets variables
 797 ! in common block mie to zero, except for wavelengths.
 798 ! OUTPUT: in module_fastj_cmnmie
 799 !       wavmid ! fast-J wavelengths (cm)
 800 !       tauaer ! aerosol optical depth
 801 !       waer  ! aerosol single scattering albedo
 802 !       gaer  ! aerosol asymmetery factor
 803 !       extaer ! aerosol extinction
 804 !       l2,l3,l4,l5,l6,l7 ! Legendre coefficients, numbered 0,1,2,......
 805 !       sizeaer ! average wet radius
 806 ! INPUTS
 807 !       lpar = total number of vertical layers in chemistry model.  Passed
 808 !         via module_fastf_cmnh
 809 !       NB = total vertical layers + 1 considered by FastJ (lpar+1=kmaxd+1).
 810 !          passed via module_fast       j_data
 811 !------------------------------------------------------------------------
 812
 813         USE module_data_mosaic_other, only : kmaxd
 814
 815         IMPLICIT NONE
 816 !jdf
 817         integer lpar             ! Number of levels in CTM
 818         integer nspint           ! Num of spectral intervals across solar
 819         parameter ( nspint = 4 ) ! spectrum for FAST-J
 820         real, dimension (nspint),save :: wavmid !cm
 821         real, dimension (nspint, kmaxd+1) :: sizeaer,extaer,waer,gaer,tauaer
 822         real, dimension (nspint, kmaxd+1) :: l2,l3,l4,l5,l6,l7
 823         data wavmid     &
 824             / 0.30e-4, 0.40e-4, 0.60e-4 ,0.999e-04 /
 825 !jdf
 826
 827 ! LOCAL VARIABLES
 828         integer klevel   ! vertical level index
 829         integer ns       ! spectral loop index
 830
 831
 832 !  aerosol optical properties: set everything = 0 when no aerosol
 833         do 1000 ns=1,nspint
 834         do 1000 klevel = 1, lpar
 835           tauaer(ns,klevel)=0.
 836           waer(ns,klevel)=0.
 837           gaer(ns,klevel)=0.
 838           sizeaer(ns,klevel)=0.0
 839           extaer(ns,klevel)=0.0
 840           l2(ns,klevel)=0.0
 841           l3(ns,klevel)=0.0
 842           l4(ns,klevel)=0.0
 843           l5(ns,klevel)=0.0
 844           l6(ns,klevel)=0.0
 845           l7(ns,klevel)=0.0
 846 1000    continue
 847
 848         return
 849         end subroutine set_common_mie
 850 !***********************************************************************
 851     subroutine fastj(isvode,jsvode,lat,lon,surface_pressure,surface_albedo,   &
 852       julian_day,tau1, pressure, temperature, optical_depth, my_ozone1,   &
 853       iozone1,tauaer_550_1,jvalue,SZA_dum,lpar,jpnl, &
 854       sizeaer,extaer,waer,gaer,tauaer,l2,l3,l4,l5,l6,l7)
 855 ! input:
 856 !       lat = latitute; must be real*8
 857 !       lon = longitude; must be real*8
 858 !       surface_pressure (mb); real*4
 859 !       surface_albedo (broadband albedo); real*4
 860 !       julian_day; integer
 861 !       tau1 = time of calculation (GMT); real*4
 862 !       pressure (mb) = vector of pressure values, pressure(NB);
 863 !          real*4; NB is the number of model layers;
 864 !          pressure (NB+1) is defined as 0 mb in model
 865 !       temperature (degree K)= vector of temperature values, temperature(NB);
 866 !          real*4
 867 !       optical_depth (dimensionless) = vector of cloud optical depths,
 868 !          optical_depth(NB); real*4
 869 !       my_ozone1 (volume mixing ratio) = ozone at layer center
 870 !          ozone(iozone); real*4; if iozone1 <= NB-1, then climatology is
 871 !          used in the upper layers
 872 !       tauaer_550; real*4  aerosol optical thickness @ 550 nm
 873 ! input note: NB is the number of model layers -- photolysis rates are calculated
 874 !       at layer centers while pressures  are given at the boundaries defining
 875 !       the top and bottom of the layers.  The number of pressure values =
 876 !       (number of layers) + 1 ; see below
 877 !       ******** pressure 2
 878 !       layer 1 - optical depth, O3, and temperature given here
 879 !       ******** pressure 1
 880 !   temperature and o3 are defined at the layer center. optical depth is
 881 !   appropriate for the layer depth.
 882 ! output:
 883 !       jvalue = photolysis rates, an array of dimension jvalue(jpnl,jppj) where
 884 !         jpnl = # of models level at which photolysis rates are calculated
 885 !         note: level 1 = first level of model (adjacent to ground)
 886 !         jppj = # of chemical species for which photolysis rates are calculated;
 887 !         this is fixed and is not easy to change on the fly
 888 !         jpnl land jppl are defined in the common block "cmn_h.f"
 889 !       SZA_dum = solar zenith angle
 890 !-----------------------------------------------------------------------
 891
 892         USE module_data_mosaic_other, only : kmaxd
 893         USE module_fastj_data
 894
 895         IMPLICIT NONE
 896 !jdf
 897 ! Print Fast-J diagnostics if iprint /= 0
 898        integer, parameter :: iprint = 0
 899        integer, parameter :: single = 4        !compiler dependent value real*4
 900 !      integer, parameter :: double = 8        !compiler dependent value real*8
 901 ! following specific for fastj
 902 ! jppj will be gas phase mechanism dependent
 903        integer,parameter :: ipar_fastj=1,jpar=1
 904 !      integer,parameter :: jppj=14        !Number of photolytic reactions supplied
 905        logical,parameter :: ldeg45=.false. !Logical flag for degraded CTM resolution
 906 ! The vertical level variables are set in fastj_driver.
 907        integer lpar           !Number of levels in CTM
 908        integer jpnl           !Number of levels requiring chemistry
 909 ! following should be available from other wrf modules and passed into
 910 ! photodriver
 911        real(kind=double), dimension(ipar_fastj) :: xgrd !  Longitude (midpoint, radians)
 912        real(kind=double), dimension(jpar) :: ygrd           !  Latitude  (midpoint, radians)
 913        real(kind=double), dimension(jpar) :: ydgrd          !  Latitude  (midpoint, degrees)
 914        real(kind=double), dimension(kmaxd+1) :: etaa    !  Eta(a) value for level boundaries
 915        real(kind=double), dimension(kmaxd+1) :: etab    !  Eta(b) value for level boundaries
 916        real(kind=double) ::  tau_fastj          !  Time of Day (hours, GMT)
 917        integer month_fastj        !  Number of month (1-12)
 918        integer iday_fastj         !  Day of year
 919        real(kind=double), dimension(ipar_fastj,jpar) ::   P , SA
 920        real(kind=double), dimension(ipar_fastj,jpar,kmaxd+1) :: T, OD
 921        real(kind=double) my_ozone(kmaxd)       !real*8 version of ozone mixing ratios
 922        real tauaer_550
 923        integer iozone
 924        integer nslat               !  Latitude of current profile point
 925        integer nslon               !  Longitude of current profile point
 926        save :: nslat, nslon
 927        integer nspint           ! Num of spectral intervals across solar
 928        parameter ( nspint = 4 ) ! spectrum for FAST-J
 929        real, dimension (nspint),save :: wavmid !cm
 930        real, dimension (nspint, kmaxd+1) :: sizeaer,extaer,waer,gaer,tauaer
 931        real, dimension (nspint, kmaxd+1) :: l2,l3,l4,l5,l6,l7
 932        data wavmid     &
 933            / 0.30e-4, 0.40e-4, 0.60e-4 ,0.999e-04 /
 934 !jdf
 935
 936     integer julian_day
 937     real surface_pressure,surface_albedo,pressure(lpar+2),   &
 938          temperature(lpar+1)
 939         real optical_depth(lpar+1)
 940         real tau1
 941     real*8 pi_fastj,lat,lon,timej,jvalue(lpar,jppj)
 942     integer isvode,jsvode
 943
 944         integer iozone1,i
 945         real my_ozone1(lpar+1)
 946
 947         real tauaer_550_1
 948         real sza_dum
 949
 950         integer ientryno_fastj
 951         save ientryno_fastj
 952         data ientryno_fastj / 0 /
 953
 954
 955
 956 !  Just focus on one column
 957       nslat = 1
 958       nslon = 1
 959       pi_fastj=3.141592653589793D0
 960 !
 961
 962 ! JCB - note that pj(NB+1) = p and is defined such elsewhere
 963 ! wig: v3.0.0 release has do loop from 1:lpar+1 but T is never
 964 !      initialized at lpar+1 leading to model instabilities. T and OD
 965 !      are ultimately never used at lpar+1. So, I changed this loop to
 966 !      1:lpar and then just assign pj a value at lpar+1.
 967         do i=1,lpar
 968            pj(i)=pressure(i)
 969            T(nslon,nslat,i)=temperature(i)
 970            OD(nslon,nslat,i)=optical_depth(i)
 971         enddo
 972         pj(lpar+1) = pressure(i)
 973
 974 ! surface albedo
 975         SA(nslon,nslat)=surface_albedo
 976 !
 977         iozone=iozone1
 978         do i=1,iozone1
 979        my_ozone(i)=my_ozone1(i)
 980         enddo
 981 !
 982         tau_fastj=tau1 ! fix time
 983         iday_fastj=julian_day
 984 ! fix optical depth for situations where aerosols not considered
 985         tauaer_550=tauaer_550_1
 986 !
 987       month_fastj=int(dble(iday_fastj)*12.d0/365.d0)+1    !  Approximately
 988       xgrd(nslon)=lon*pi_fastj/180.d0
 989       ygrd(nslat)=lat*pi_fastj/180.d0
 990       ydgrd(nslat)=lat
 991
 992 !  Initial call to Fast-J to set things up--done only once
 993         if (ientryno_fastj .eq. 0) then
 994             call inphot2
 995             ientryno_fastj = 1
 996         end if
 997 !
 998 !  Now call fastj as appropriate
 999         timej=0.0 ! manually set offset to zero (JCB: 14 November 2001)
1000     call photoj(isvode,jsvode,jvalue,timej,nslat,nslon,iozone,tauaer_550, &
1001         my_ozone,p,t,od,sa,lpar,jpnl, &
1002         xgrd,ygrd,tau_fastj,month_fastj,iday_fastj,ydgrd, &
1003         sizeaer,extaer,waer,gaer,tauaer,l2,l3,l4,l5,l6,l7)
1004         sza_dum=SZA
1005
1006         return
1007       end subroutine fastj
1008 !**********************************************************************
1009 !---(pphot.f)-------generic CTM shell from UCIrvine (p-code 4.0, 7/99)
1010 !---------PPHOT calculates photolysis rates with the Fast-J scheme
1011 !---------subroutines:  inphot, photoj, Fast-J schemes...
1012 !-----------------------------------------------------------------------
1013 !
1014       subroutine inphot2
1015 !-----------------------------------------------------------------------
1016 !  Routine to initialise photolysis rate data, called directly from the
1017 !  cinit routine in ASAD. Currently use it to read the JPL spectral data
1018 !  and standard O3 and T profiles and to set the appropriate reaction index.
1019 !-----------------------------------------------------------------------
1020 !
1021 !     iph       Channel number for reading all data files
1022 !     rad       Radius of Earth (cm)
1023 !     zzht      Effective scale height above top of atmosphere (cm)
1024 !     dtaumax   Maximum opt.depth above which a new level should be inserted
1025 !     dtausub   No. of opt.depths at top of cloud requiring subdivision
1026 !     dsubdiv   Number of additional levels to add at top of cloud
1027 !     szamax    Solar zenith angle cut-off, above which to skip calculation
1028 !
1029 !-----------------------------------------------------------------------
1030 !
1031
1032         USE module_data_mosaic_other, only : kmaxd
1033         USE module_fastj_data
1034
1035         IMPLICIT NONE
1036 !jdf
1037 ! Print Fast-J diagnostics if iprint /= 0
1038         integer, parameter :: iprint = 0
1039         integer, parameter :: single = 4        !compiler dependent value real*4
1040 !       integer, parameter :: double = 8        !compiler dependent value real*8
1041        integer,parameter :: ipar_fastj=1,jpar=1
1042 !      integer,parameter :: jppj=14        !Number of photolytic reactions supplied
1043        logical,parameter :: ldeg45=.false. !Logical flag for degraded CTM resolution
1044        integer lpar           !Number of levels in CTM
1045        integer jpnl           !Number of levels requiring chemistry
1046        real(kind=double), dimension(ipar_fastj) :: xgrd !  Longitude (midpoint, radians)
1047        real(kind=double), dimension(jpar) :: ygrd           !  Latitude  (midpoint, radians)
1048        real(kind=double), dimension(jpar) :: ydgrd          !  Latitude  (midpoint, degrees)
1049        real(kind=double), dimension(kmaxd+1) :: etaa    !  Eta(a) value for level boundaries
1050        real(kind=double), dimension(kmaxd+1) :: etab    !  Eta(b) value for level boundaries
1051        real(kind=double) ::  tau_fastj          !  Time of Day (hours, GMT)
1052        integer month_fastj        !  Number of month (1-12)
1053        integer iday_fastj         !  Day of year
1054 !jdf
1055 !
1056 ! Set labels of photolysis rates required
1057 !ec032504      CALL RD_JS(iph,path_fastj_ratjd)
1058 !       call rd_js2
1059 !
1060 ! Read in JPL spectral data set
1061 !ec032504      CALL RD_TJPL(iph,path_fastj_jvspec)
1062         call rd_tjpl2
1063 !
1064 ! Read in T & O3 climatology
1065 !ec032504      CALL RD_PROF(iph,path_fastj_jvatms)
1066 !       call rd_prof2
1067 !
1068 ! Select Aerosol/Cloud types to be used
1069 !         call set_aer2
1070
1071       return
1072       end subroutine inphot2
1073 !*************************************************************************
1074
1075       subroutine photoj(isvode,jsvode,zpj,timej,nslat,nslon,iozone,tauaer_550_1, &
1076         my_ozone,p,t,od,sa,lpar,jpnl,xgrd,ygrd,tau_fastj,month_fastj,iday_fastj, &
1077         ydgrd,sizeaer,extaer,waer,gaer,tauaer,l2,l3,l4,l5,l6,l7)
1078 !-----------------------------------------------------------------------
1079 !----jv_trop.f:  new FAST J-Value code, troposphere only (mjprather 6/96)
1080 !----     uses special wavelength quadrature spectral data (jv_spec.dat)
1081 !---      that includes only 289 nm - 800 nm  (later a single 205 nm add-on)
1082 !---      uses special compact Mie code based on Feautrier/Auer/Prather vers.
1083 !-----------------------------------------------------------------------
1084 !
1085 !     zpj      External array providing J-values to main CTM code
1086 !     timej    Offset in hours from start of timestep to time J-values
1087 !              required for - take as half timestep for mid-step Js.
1088 !     solf     Solar distance factor, for scaling; normally given by:
1089 !                      1.0-(0.034*cos(real(iday_fastj-186)*2.0*pi_fastj/365.))
1090 !
1091 !----------basic common blocks:-----------------------------------------
1092
1093         USE module_data_mosaic_other, only : kmaxd
1094         USE module_fastj_data
1095
1096         IMPLICIT NONE
1097 !jdf
1098 ! Print Fast-J diagnostics if iprint /= 0
1099        integer, parameter :: iprint = 0
1100        integer, parameter :: single = 4        !compiler dependent value real*4
1101 !      integer, parameter :: double = 8        !compiler dependent value real*8
1102        integer,parameter :: ipar_fastj=1,jpar=1
1103 !      integer,parameter :: jppj=14        !Number of photolytic reactions supplied
1104        logical,parameter :: ldeg45=.false. !Logical flag for degraded CTM resolution
1105        integer lpar           !Number of levels in CTM
1106        integer jpnl           !Number of levels requiring chemistry
1107        real(kind=double), dimension(ipar_fastj) :: xgrd !  Longitude (midpoint, radians)
1108        real(kind=double), dimension(jpar) :: ygrd           !  Latitude  (midpoint, radians)
1109        real(kind=double), dimension(jpar) :: ydgrd          !  Latitude  (midpoint, degrees)
1110        real(kind=double), dimension(kmaxd+1) :: etaa    !  Eta(a) value for level boundaries
1111        real(kind=double), dimension(kmaxd+1) :: etab    !  Eta(b) value for level boundaries
1112        real(kind=double) ::  tau_fastj          !  Time of Day (hours, GMT)
1113        integer month_fastj        !  Number of month (1-12)
1114        integer iday_fastj         !  Day of year
1115        real(kind=double), dimension(ipar_fastj,jpar) ::   P , SA
1116        real(kind=double), dimension(ipar_fastj,jpar,kmaxd+1) :: T, OD
1117        real(kind=double) my_ozone(kmaxd)       !real*8 version of ozone mixing ratios
1118        real tauaer_550_1
1119        integer iozone
1120        integer nslat               !  Latitude of current profile point
1121        integer nslon               !  Longitude of current profile point
1122        integer nspint           ! Num of spectral intervals across solar
1123        parameter ( nspint = 4 ) ! spectrum for FAST-J
1124        real, dimension (nspint),save :: wavmid !cm
1125        real, dimension (nspint, kmaxd+1) :: sizeaer,extaer,waer,gaer,tauaer
1126        real, dimension (nspint, kmaxd+1) :: l2,l3,l4,l5,l6,l7
1127        data wavmid     &
1128            / 0.30e-4, 0.40e-4, 0.60e-4 ,0.999e-04 /
1129 !jdf
1130       real*8 zpj(lpar,jppj),timej,solf
1131         real*8 pi_fastj
1132
1133         integer i,j
1134     integer isvode,jsvode
1135
1136 !-----------------------------------------------------------------------
1137 !
1138       do i=1,jpnl
1139         do j=1,jppj
1140           zj(i,j)=0.D0
1141           zpj(i,j)=0.D0
1142         enddo
1143       enddo
1144 !
1145 !---Calculate new solar zenith angle
1146       CALL SOLAR2(timej,nslat,nslon, &
1147       xgrd,ygrd,tau_fastj,month_fastj,iday_fastj)
1148       if(SZA.gt.szamax) go to 10
1149 !
1150 !---Set up profiles on model levels
1151       CALL SET_PROF(isvode,jsvode,nslat,nslon,iozone,tauaer_550_1, &
1152         my_ozone,p,t,od,sa,lpar,jpnl,month_fastj,ydgrd)
1153 !
1154 !---Print out atmosphere
1155        if(iprint.ne.0)CALL PRTATM(3,nslat,nslon,tau_fastj,month_fastj,ydgrd)  ! code change jcb
1156 !       call prtatm(0)
1157 !
1158 !-----------------------------------------------------------------------
1159       CALL JVALUE(isvode,jsvode,lpar,jpnl, &
1160         ydgrd,sizeaer,extaer,waer,gaer,tauaer,l2,l3,l4,l5,l6,l7)
1161 !-----------------------------------------------------------------------
1162 !---Print solar flux terms
1163 !      WRITE(6,'(A16,I5,20I9)') '  wave (beg/end)',(i,i=1,jpnl)
1164 !      DO j=NW1,NW2
1165 !        WRITE(6,'(2F8.2,20F9.6)') WBIN(j),WBIN(j+1),
1166 !     $                            (FFF(j,i)/FL(j),i=1,jpnl)
1167 !      ENDDO
1168 !
1169 !---Include variation in distance to sun
1170         pi_fastj=3.1415926536d0
1171        solf=1.d0-(0.034d0*cos(dble(iday_fastj-186)*2.d0   &
1172              *pi_fastj/365.d0))
1173         if(iprint.ne.0)then
1174 !          write(6,'('' solf = '', f10.5)')solf
1175           write(*,'('' solf = '', f10.5)')solf
1176         endif
1177 !      solf=1.d0 ! code change jcb
1178 !-----------------------------------------------------------------------
1179       CALL JRATET(solf,nslat,nslon,p,t,od,sa,lpar,jpnl)
1180 !-----------------------------------------------------------------------
1181 !
1182
1183 !  "zj" updated in JRATET - pass this back to ASAD as "zpj"
1184       do i=1,jpnl
1185         do j=1,jppj
1186           zpj(i,j)= zj(i,j)
1187         enddo
1188       enddo
1189
1190 !
1191 !---Output selected values
1192   10  if((.not.ldeg45.and.nslon.eq.37.and.nslat.eq.36).or.   &
1193                (ldeg45.and.nslon.eq.19.and.nslat.eq.18)) then
1194         i=min(jppj,8)
1195 !        write(6,1000)iday_fastj,tau_fastj+timej,sza,jlabel(i),zpj(1,i)
1196       endif
1197 !
1198       return
1199 ! 1000 format(' Photolysis on day ',i4,' at ',f4.1,' hrs: SZA = ',f7.3,   &
1200 !                     ' J',a7,'= ',1pE10.3)
1201       end subroutine photoj
1202
1203 !*************************************************************************
1204       subroutine set_prof(isvode,jsvode,nslat,nslon,iozone,tauaer_550, &
1205         my_ozone,p,t,od,sa,lpar,jpnl,month_fastj,ydgrd)
1206 !-----------------------------------------------------------------------
1207 !  Routine to set up atmospheric profiles required by Fast-J using a
1208 !  doubled version of the level scheme used in the CTM. First pressure
1209 !  and z* altitude are defined, then O3 and T are taken from the supplied
1210 !  climatology and integrated to the CTM levels (may be overwritten with
1211 !  values directly from the CTM, if desired) and then black carbon and
1212 !  aerosol profiles are constructed.
1213 !                                                     Oliver (04/07/99)
1214 !-----------------------------------------------------------------------
1215 !
1216 !     pj       Pressure at boundaries of model levels (hPa)
1217 !     z        Altitude of boundaries of model levels (cm)
1218 !     odcol    Optical depth at each model level
1219 !     masfac   Conversion factor for pressure to column density
1220 !
1221 !     TJ       Temperature profile on model grid
1222 !     DM       Air column for each model level (molecules.cm-2)
1223 !     DO3      Ozone column for each model level (molecules.cm-2)
1224 !     DBC      Mass of Black Carbon at each model level (g.cm-3)  !  .....!
1225 !     PSTD     Approximate pressures of levels for supplied climatology
1226 !
1227 !-----------------------------------------------------------------------
1228
1229       USE module_data_mosaic_other, only : kmaxd
1230       USE module_fastj_data
1231
1232       IMPLICIT NONE
1233 !jdf
1234 ! Print Fast-J diagnostics if iprint /= 0
1235       integer, parameter :: iprint = 0
1236       integer, parameter :: single = 4        !compiler dependent value real*4
1237 !     integer, parameter :: double = 8        !compiler dependent value real*8
1238       integer,parameter :: ipar_fastj=1,jpar=1
1239 !     integer,parameter :: jppj=14        !Number of photolytic reactions supplied
1240       logical,parameter :: ldeg45=.false. !Logical flag for degraded CTM resolution
1241       integer lpar           !Number of levels in CTM
1242       integer jpnl           !Number of levels requiring chemistry
1243       real(kind=double), dimension(ipar_fastj) :: xgrd !  Longitude (midpoint, radians)
1244       real(kind=double), dimension(jpar) :: ygrd           !  Latitude  (midpoint, radians)
1245       real(kind=double), dimension(jpar) :: ydgrd          !  Latitude  (midpoint, degrees)
1246       real(kind=double), dimension(kmaxd+1) :: etaa    !  Eta(a) value for level boundaries
1247       real(kind=double), dimension(kmaxd+1) :: etab    !  Eta(b) value for level boundaries
1248       real(kind=double) ::  tau_fastj          !  Time of Day (hours, GMT)
1249       integer month_fastj        !  Number of month (1-12)
1250       integer iday_fastj         !  Day of year
1251       real(kind=double), dimension(ipar_fastj,jpar) ::   P , SA
1252       real(kind=double), dimension(ipar_fastj,jpar,kmaxd+1) :: T, OD
1253       real(kind=double) my_ozone(kmaxd)       !real*8 version of ozone mixing ratios
1254       real tauaer_550
1255       integer iozone
1256       integer nslat               !  Latitude of current profile point
1257       integer nslon               !  Longitude of current profile point
1258 !jdf
1259       real*8  pstd(52),oref2(51),tref2(51),bref2(51)
1260       real*8  odcol(lpar),dlogp,f0,t0,b0,pb,pc,xc,masfac,scaleh
1261       real vis, aerd1, aerd2
1262
1263       integer  i, k, l, m
1264       integer isvode,jsvode
1265
1266        pj(NB+1) = 0.d0 ! define top level
1267 !
1268 !  Set up cloud and surface properties
1269       call CLDSRF(isvode,jsvode,odcol,nslat,nslon,p,t,od,sa,lpar,jpnl)
1270
1271 !  Mass factor - delta-Pressure (mbars) to delta-Column (molecules.cm-2)
1272       masfac=100.d0*6.022d+23/(28.97d0*9.8d0*10.d0)
1273 !
1274 !  Set up pressure levels for O3/T climatology - assume that value
1275 !  given for each 2 km z* level applies from 1 km below to 1 km above,
1276 !  so select pressures at these boundaries. Surface level values at
1277 !  1000 mb are assumed to extend down to the actual P(nslon,nslat).
1278 !
1279       pstd(1) = max(pj(1),1000.d0)
1280       pstd(2) = 1000.d0*10.d0**(-1.d0/16.d0)
1281       dlogp = 10.d0**(-2.d0/16.d0)
1282       do i=3,51
1283         pstd(i) = pstd(i-1)*dlogp
1284       enddo
1285       pstd(52) = 0.d0
1286 !
1287 !  Select appropriate monthly and latitudinal profiles
1288       m = max(1,min(12,month_fastj))
1289       l = max(1,min(18,(int(ydgrd(nslat))+99)/10))
1290 !
1291 !  Temporary arrays for climatology data
1292       do i=1,51
1293         oref2(i)=oref(i,l,m)
1294         tref2(i)=tref(i,l,m)
1295         bref2(i)=bref(i)
1296       enddo
1297 !
1298 !  Apportion O3 and T on supplied climatology z* levels onto CTM levels
1299 !  with mass (pressure) weighting, assuming constant mixing ratio and
1300 !  temperature half a layer on either side of the point supplied.
1301 !
1302       do i = 1,NB
1303         F0 = 0.d0
1304         T0 = 0.d0
1305         B0 = 0.d0
1306         do k = 1,51
1307           PC = min(pj(i),pstd(k))
1308           PB = max(pj(i+1),pstd(k+1))
1309           if(PC.gt.PB) then
1310             XC = (PC-PB)/(pj(i)-pj(i+1))
1311             F0 = F0 + oref2(k)*XC
1312             T0 = T0 + tref2(k)*XC
1313             B0 = B0 + bref2(k)*XC
1314           endif
1315         enddo
1316         TJ(i) = T0
1317         DO3(i)= F0*1.d-6
1318         DBC(i)= B0
1319       enddo
1320 !
1321 !  Insert model values here to replace or supplement climatology.
1322 !  Note that CTM temperature is always used in x-section calculations
1323 !  (see JRATET); TJ is used in actinic flux calculation only.
1324 !
1325         do i=1,lpar ! JCB code change; just use climatlogy for upper levels
1326         if(i.le.iozone)DO3(i) = my_ozone(i) ! Volume Mixing Ratio
1327 !        TJ(i)  = T(nslon,nslat,I)   ! Kelvin
1328 ! JCB - overwrite climatology
1329 !       TJ(i)  = (T(nslon,nslat,i)+T(nslon,nslat,i+1))/2. ! JCB - take midpoint
1330 ! code change - now take temperature as appropriate for midpoint of layer
1331         TJ(i)=T(nslon,nslat,i)
1332         enddo
1333         if(lpar+1.le.iozone)then
1334         DO3(lpar+1) = my_ozone(lpar+1)  ! Above top of model (or use climatology)
1335         endif
1336 !       TJ(lpar+1)  = my_temp(lpar)   ! Above top of model (or use climatology)
1337 !wig 26-Aug-2000: Comment out following line so that climatology is used for
1338 !                 above the model top.
1339 !       TJ(lpar+1) = T(nslon,nslat,NB)  ! JCB - just use climatology or given temperature
1340 ! JCB read in O3
1341 !
1342 !
1343 !  Calculate effective altitudes using scale height at each level
1344       z(1) = 0.d0
1345       do i=1,lpar
1346         scaleh=1.3806d-19*masfac*TJ(i)
1347         z(i+1) = z(i)-(log(pj(i+1)/pj(i))*scaleh)
1348       enddo
1349 !
1350 !  Add Aerosol Column - include aerosol types here. Currently use soot
1351 !  water and ice; assume black carbon x-section of 10 m2/g, independent
1352 !  of wavelength; assume limiting temperature for ice of -40 deg C.
1353
1354       do i=1,lpar
1355 !        AER(1,i) = DBC(i)*10.d0*(z(i+1)-z(i)) ! DBC must be g/m^3
1356 ! calculate AER(1,i) according to aerosol density - use trap rule
1357         vis=23.0
1358         call aeroden(z(i)/100000.,vis,aerd1) ! convert cm to km
1359         call aeroden(z(i+1)/100000.,vis,aerd2)
1360 ! trap rule used here; convert cm to km; divide by 100000.
1361         AER(1,i)=(z(i+1)-z(i))/100000.*(aerd1+aerd2)/2./4287.55*tauaer_550
1362 !       write(6,*)i,z(i)/100000.,aerd1,aerd2,tauaer_550,AER(1,i)
1363 !
1364         if(T(nslon,nslat,I).gt.233.d0) then
1365           AER(2,i) = odcol(i)
1366           AER(3,i) = 0.d0
1367         else
1368           AER(2,i) = 0.d0
1369           AER(3,i) = odcol(i)
1370         endif
1371       enddo
1372       do k=1,MX
1373         AER(k,lpar+1) = 0.d0 ! just set equal to zero
1374       enddo
1375
1376         AER(1,lpar+1)=2.0*AER(1,lpar) ! kludge
1377 !
1378 !  Calculate column quantities for Fast-J
1379       do i=1,NB
1380         DM(i)  = (PJ(i)-PJ(i+1))*masfac
1381         DO3(i) = DO3(i)*DM(i)
1382       enddo
1383 !
1384       end subroutine set_prof
1385
1386 !******************************************************************
1387 !     SUBROUTINE CLDSRF(odcol)
1388       SUBROUTINE CLDSRF(isvode,jsvode,odcol,nslat,nslon,p,t,od,sa, &
1389          lpar,jpnl)
1390 !-----------------------------------------------------------------------
1391 !  Routine to set cloud and surface properties
1392 !-----------------------------------------------------------------------
1393 !     rflect   Surface albedo (Lambertian)
1394 !     odmax    Maximum allowed optical depth, above which they are scaled
1395 !     odcol    Optical depth at each model level
1396 !     odsum    Column optical depth
1397 !     nlbatm   Level of lower photolysis boundary - usually surface
1398 !-----------------------------------------------------------------------
1399
1400         USE module_data_mosaic_other, only : kmaxd
1401         USE module_fastj_data
1402
1403         IMPLICIT NONE
1404 !jdf
1405 ! Print Fast-J diagnostics if iprint /= 0
1406       integer, parameter :: iprint = 0
1407       integer, parameter :: single = 4        !compiler dependent value real*4
1408 !     integer, parameter :: double = 8        !compiler dependent value real*8
1409       integer,parameter :: ipar_fastj=1,jpar=1
1410 !     integer,parameter :: jppj=14        !Number of photolytic reactions supplied
1411       logical,parameter :: ldeg45=.false. !Logical flag for degraded CTM resolution
1412       integer lpar           !Number of levels in CTM
1413       integer jpnl           !Number of levels requiring chemistry
1414       real(kind=double), dimension(ipar_fastj) :: xgrd !  Longitude (midpoint, radians)
1415       real(kind=double), dimension(jpar) :: ygrd           !  Latitude  (midpoint, radians)
1416       real(kind=double), dimension(jpar) :: ydgrd          !  Latitude  (midpoint, degrees)
1417       real(kind=double), dimension(kmaxd+1) :: etaa    !  Eta(a) value for level boundaries
1418       real(kind=double), dimension(kmaxd+1) :: etab    !  Eta(b) value for level boundaries
1419       real(kind=double) ::  tau_fastj          !  Time of Day (hours, GMT)
1420       integer month_fastj        !  Number of month (1-12)
1421       integer iday_fastj         !  Day of year
1422       integer nslat               !  Latitude of current profile point
1423       integer nslon               !  Longitude of current profile point
1424       real(kind=double), dimension(ipar_fastj,jpar) ::   P , SA
1425       real(kind=double), dimension(ipar_fastj,jpar,kmaxd+1) :: T, OD
1426 !jdf
1427       integer i, j, k
1428       integer isvode, jsvode
1429       real*8  odcol(lpar), odsum, odmax, odtot
1430 !
1431 ! Default lower photolysis boundary as bottom of level 1
1432       nlbatm = 1
1433 !
1434 ! Set surface albedo
1435       RFLECT = dble(SA(nslon,nslat))
1436       RFLECT = max(0.d0,min(1.d0,RFLECT))
1437 !
1438 ! Zero aerosol column
1439       do k=1,MX
1440         do i=1,NB
1441           AER(k,i) = 0.d0
1442         enddo
1443       enddo
1444 !
1445 ! Scale optical depths as appropriate - limit column to 'odmax'
1446       odmax = 200.d0
1447       odsum =   0.d0
1448       do i=1,lpar
1449         odcol(i) = dble(OD(nslon,nslat,i))
1450         odsum    = odsum + odcol(i)
1451 !       if(isvode.eq.8.and.jsvode.eq.1) print*,'jdf odcol',odcol(i),odsum
1452       enddo
1453       if(odsum.gt.odmax) then
1454         odsum = odmax/odsum
1455         do i=1,lpar
1456           odcol(i) = odcol(i)*odsum
1457         enddo
1458         odsum = odmax
1459       endif
1460 ! Set sub-division switch if appropriate
1461       odtot=0.d0
1462       jadsub(nb)=0
1463       jadsub(nb-1)=0
1464       do i=nb-1,1,-1
1465         k=2*i
1466         jadsub(k)=0
1467         jadsub(k-1)=0
1468         odtot=odtot+odcol(i)
1469         if(odcol(i).gt.0.d0.and.dtausub.gt.0.d0) then
1470           if(odtot.le.dtausub) then
1471             jadsub(k)=1
1472             jadsub(k-1)=1
1473 !       if(isvode.eq.8.and.jsvode.eq.1) print*,'jdf in cldsf1',i,k,jadsub(k)
1474           elseif(odtot.gt.dtausub) then
1475             jadsub(k)=1
1476             jadsub(k-1)=0
1477             do j=1,2*(i-1)
1478               jadsub(j)=0
1479             enddo
1480 !       if(isvode.eq.8.and.jsvode.eq.1) print*,'jdf in cldsf2',i,k,jadsub(k)
1481             go to 20
1482           endif
1483         endif
1484       enddo
1485  20   continue
1486 !
1487       return
1488       end SUBROUTINE CLDSRF
1489
1490 !********************************************************************
1491       subroutine solar2(timej,nslat,nslon, &
1492       xgrd,ygrd,tau_fastj,month_fastj,iday_fastj)
1493 !-----------------------------------------------------------------------
1494 !  Routine to set up SZA for given lat, lon and time
1495 !-----------------------------------------------------------------------
1496 !     timej    Offset in hours from start of timestep to time J-values
1497 !              required for - take as half timestep for mid-step Js.
1498 !-----------------------------------------------------------------------
1499
1500         USE module_data_mosaic_other, only : kmaxd
1501         USE module_fastj_data
1502         IMPLICIT NONE
1503 !jdf
1504 ! Print Fast-J diagnostics if iprint /= 0
1505       integer, parameter :: iprint = 0
1506       integer, parameter :: single = 4        !compiler dependent value real*4
1507 !     integer, parameter :: double = 8        !compiler dependent value real*8
1508       integer,parameter :: ipar_fastj=1,jpar=1
1509 !     integer,parameter :: jppj=14        !Number of photolytic reactions supplied
1510       logical,parameter :: ldeg45=.false. !Logical flag for degraded CTM resolution
1511       integer lpar           !Number of levels in CTM
1512       integer jpnl           !Number of levels requiring chemistry
1513       real(kind=double), dimension(ipar_fastj) :: xgrd !  Longitude (midpoint, radians)
1514       real(kind=double), dimension(jpar) :: ygrd           !  Latitude  (midpoint, radians)
1515       real(kind=double), dimension(jpar) :: ydgrd          !  Latitude  (midpoint, degrees)
1516       real(kind=double), dimension(kmaxd+1) :: etaa    !  Eta(a) value for level boundaries
1517       real(kind=double), dimension(kmaxd+1) :: etab    !  Eta(b) value for level boundaries
1518       real(kind=double) ::  tau_fastj          !  Time of Day (hours, GMT)
1519       integer month_fastj        !  Number of month (1-12)
1520       integer iday_fastj         !  Day of year
1521       integer nslat               !  Latitude of current profile point
1522       integer nslon               !  Longitude of current profile point
1523 !jdf
1524       real*8 pi_fastj, pi180, loct, timej
1525       real*8 sindec, soldek, cosdec, sinlat, sollat, coslat, cosz
1526 !
1527       pi_fastj=3.141592653589793D0
1528       pi180=pi_fastj/180.d0
1529       sindec=0.3978d0*sin(0.9863d0*(dble(iday_fastj)-80.d0)*pi180)
1530       soldek=asin(sindec)
1531       cosdec=cos(soldek)
1532       sinlat=sin(ygrd(nslat))
1533       sollat=asin(sinlat)
1534       coslat=cos(sollat)
1535 !
1536       loct = (((tau_fastj+timej)*15.d0)-180.d0)*pi180 + xgrd(nslon)
1537       cosz = cosdec*coslat*cos(loct) + sindec*sinlat
1538       sza  = acos(cosz)/pi180
1539       U0 = cos(SZA*pi180)
1540 !
1541       return
1542       end subroutine solar2
1543
1544 !**********************************************************************
1545
1546
1547       SUBROUTINE JRATET(SOLF,nslat,nslon,p,t,od,sa,lpar,jpnl)
1548 !-----------------------------------------------------------------------
1549 !  Calculate and print J-values. Note that the loop in this routine
1550 !  only covers the jpnl levels actually needed by the CTM.
1551 !-----------------------------------------------------------------------
1552 !
1553 !     FFF    Actinic flux at each level for each wavelength bin
1554 !     QQQ    Cross sections for species (read in in RD_TJPL)
1555 !     SOLF   Solar distance factor, for scaling; normally given by:
1556 !                      1.0-(0.034*cos(real(iday_fastj-186)*2.0*pi_fastj/365.))
1557 !                      Assumes aphelion day 186, perihelion day 3.
1558 !     TQQ    Temperatures at which QQQ cross sections supplied
1559 !
1560 !-----------------------------------------------------------------------
1561
1562         USE module_data_mosaic_other, only : kmaxd
1563         USE module_fastj_data
1564
1565         IMPLICIT NONE
1566 !jdf
1567 ! Print Fast-J diagnostics if iprint /= 0
1568       integer, parameter :: iprint = 0
1569       integer, parameter :: single = 4        !compiler dependent value real*4
1570 !     integer, parameter :: double = 8        !compiler dependent value real*8
1571       integer,parameter :: ipar_fastj=1,jpar=1
1572 !     integer,parameter :: jppj=14        !Number of photolytic reactions supplied
1573       logical,parameter :: ldeg45=.false. !Logical flag for degraded CTM resolution
1574       integer lpar           !Number of levels in CTM
1575       integer jpnl           !Number of levels requiring chemistry
1576       real(kind=double), dimension(ipar_fastj) :: xgrd !  Longitude (midpoint, radians)
1577       real(kind=double), dimension(jpar) :: ygrd           !  Latitude  (midpoint, radians)
1578       real(kind=double), dimension(jpar) :: ydgrd          !  Latitude  (midpoint, degrees)
1579       real(kind=double), dimension(kmaxd+1) :: etaa    !  Eta(a) value for level boundaries
1580       real(kind=double), dimension(kmaxd+1) :: etab    !  Eta(b) value for level boundaries
1581       real(kind=double) ::  tau_fastj          !  Time of Day (hours, GMT)
1582       integer month_fastj        !  Number of month (1-12)
1583       integer iday_fastj         !  Day of year
1584       integer nslat               !  Latitude of current profile point
1585       integer nslon               !  Longitude of current profile point
1586       real(kind=double), dimension(ipar_fastj,jpar) ::   P , SA
1587       real(kind=double), dimension(ipar_fastj,jpar,kmaxd+1) :: T, OD
1588 !jdf
1589       integer i, j, k
1590       real*8 qo2tot, qo3tot, qo31d, qo33p, qqqt
1591 !      real*8 xseco2, xseco3, xsec1d, solf, tfact
1592       real*8  solf, tfact
1593 !
1594       do I=1,jpnl
1595        VALJ(1) = 0.d0
1596        VALJ(2) = 0.d0
1597        VALJ(3) = 0.d0
1598        do K=NW1,NW2                       ! Using model 'T's here
1599          QO2TOT= xseco2(K,dble(T(nslon,nslat,I)))
1600          VALJ(1) = VALJ(1) + QO2TOT*FFF(K,I)
1601          QO3TOT= xseco3(K,dble(T(nslon,nslat,I)))
1602          QO31D = xsec1d(K,dble(T(nslon,nslat,I)))*QO3TOT
1603          QO33P = QO3TOT - QO31D
1604          VALJ(2) = VALJ(2) + QO33P*FFF(K,I)
1605          VALJ(3) = VALJ(3) + QO31D*FFF(K,I)
1606        enddo
1607 !------Calculate remaining J-values with T-dep X-sections
1608        do J=4,NJVAL
1609          VALJ(J) = 0.d0
1610          TFACT = 0.d0
1611          if(TQQ(2,J).gt.TQQ(1,J)) TFACT = max(0.d0,min(1.d0,   &
1612               (T(nslon,nslat,I)-TQQ(1,J))/(TQQ(2,J)-TQQ(1,J)) ))
1613          do K=NW1,NW2
1614            QQQT = QQQ(K,1,J-3) + (QQQ(K,2,J-3) - QQQ(K,1,J-3))*TFACT
1615            VALJ(J) = VALJ(J) + QQQT*FFF(K,I)
1616 !------Additional code for pressure dependencies
1617 !           if(jpdep(J).ne.0) then
1618 !             VALJ(J) = VALJ(J) + QQQT*FFF(K,I)*
1619 !     $                   (zpdep(K,L)*(pj(i)+pj(i+1))*0.5d0)
1620 !           endif
1621          enddo
1622        enddo
1623        do j=1,jppj
1624          zj(i,j)=VALJ(jind(j))*jfacta(j)*SOLF
1625        enddo
1626 !  Herzberg bin
1627        do j=1,nhz
1628          zj(i,hzind(j))=hztoa(j)*fhz(i)*SOLF
1629        enddo
1630       enddo
1631       return
1632       end SUBROUTINE JRATET
1633
1634 !*********************************************************************
1635
1636
1637       SUBROUTINE PRTATM(N,nslat,nslon,tau_fastj,month_fastj,ydgrd)
1638 !-----------------------------------------------------------------------
1639 !  Print out the atmosphere and calculate appropriate columns
1640 !     N=1    Print out column totals only
1641 !     N=2    Print out full columns
1642 !     N=3    Print out full columns and climatology
1643 !-----------------------------------------------------------------------
1644
1645         USE module_data_mosaic_other, only : kmaxd
1646         USE module_fastj_data
1647
1648 !jdf
1649 ! Print Fast-J diagnostics if iprint /= 0
1650       integer, parameter :: iprint = 0
1651       integer, parameter :: single = 4        !compiler dependent value real*4
1652 !     integer, parameter :: double = 8        !compiler dependent value real*8
1653       integer,parameter :: ipar_fastj=1,jpar=1
1654 !     integer,parameter :: jppj=14        !Number of photolytic reactions supplied
1655       logical,parameter :: ldeg45=.false. !Logical flag for degraded CTM resolution
1656       integer lpar           !Number of levels in CTM
1657       integer jpnl           !Number of levels requiring chemistry
1658       real(kind=double), dimension(ipar_fastj) :: xgrd !  Longitude (midpoint, radians)
1659       real(kind=double), dimension(jpar) :: ygrd           !  Latitude  (midpoint, radians)
1660       real(kind=double), dimension(jpar) :: ydgrd          !  Latitude  (midpoint, degrees)
1661       real(kind=double), dimension(kmaxd+1) :: etaa    !  Eta(a) value for level boundaries
1662       real(kind=double), dimension(kmaxd+1) :: etab    !  Eta(b) value for level boundaries
1663       real(kind=double) ::  tau_fastj          !  Time of Day (hours, GMT)
1664       integer month_fastj        !  Number of month (1-12)
1665       integer iday_fastj         !  Day of year
1666       integer nslat               !  Latitude of current profile point
1667       integer nslon               !  Longitude of current profile point
1668 !jdf
1669       integer n, i, k, l, m
1670       real*8 COLO3(NB),COLO2(NB),COLAX(MX,NB),ZKM,ZSTAR,PJC
1671       real*8 climat(9),masfac,dlogp
1672       if(N.eq.0) return
1673 !---Calculate columns, for diagnostic output only:
1674       COLO3(NB) = DO3(NB)
1675       COLO2(NB) = DM(NB)*0.20948d0
1676       do K=1,MX
1677         COLAX(K,NB) = AER(K,NB)
1678       enddo
1679       do I=NB-1,1,-1
1680         COLO3(i) = COLO3(i+1)+DO3(i)
1681         COLO2(i) = COLO2(i+1)+DM(i)*0.20948d0
1682         do K=1,MX
1683           COLAX(k,i) = COLAX(k,i+1)+AER(k,i)
1684         enddo
1685       enddo
1686       write(*,1200) ' Tau=',tau_fastj,'  SZA=',sza
1687       write(*,1200) ' O3-column(DU)=',COLO3(1)/2.687d16,   &
1688                     '  column aerosol @1000nm=',(COLAX(K,1),K=1,MX)
1689 !---Print out atmosphere
1690       if(N.gt.1) then
1691         write(*,1000) (' AER-X ','col-AER',k=1,mx)
1692         do I=NB,1,-1
1693           PJC = PJ(I)
1694           ZKM =1.d-5*Z(I)
1695           ZSTAR = 16.d0*DLOG10(1013.d0/PJC)
1696           write(*,1100) I,ZKM,ZSTAR,DM(I),DO3(I),1.d6*DO3(I)/DM(I),   &
1697                TJ(I),PJC,COLO3(I),COLO2(I),(AER(K,I),COLAX(K,I),K=1,MX)
1698         enddo
1699       endif
1700 !
1701 !---Print out climatology
1702       if(N.gt.2) then
1703         do i=1,9
1704           climat(i)=0.d0
1705         enddo
1706         m = max(1,min(12,month_fastj))
1707         l = max(1,min(18,(int(ydgrd(nslat))+99)/10))
1708         masfac=100.d0*6.022d+23/(28.97d0*9.8d0*10.d0)
1709         write(*,*) 'Specified Climatology'
1710         write(*,1000)
1711         do i=51,1,-1
1712           dlogp=10.d0**(-1.d0/16.d0)
1713           PJC = 1000.d0*dlogp**(2*i-2)
1714           climat(1) = 16.d0*DLOG10(1000.D0/PJC)
1715           climat(2) = climat(1)
1716           climat(3) = PJC*(1.d0/dlogp-dlogp)*masfac
1717           if(i.eq.1) climat(3)=PJC*(1.d0-dlogp)*masfac
1718           climat(4)=climat(3)*oref(i,l,m)*1.d-6
1719           climat(5)=oref(i,l,m)
1720           climat(6)=tref(i,l,m)
1721           climat(7)=PJC
1722           climat(8)=climat(8)+climat(4)
1723           climat(9)=climat(9)+climat(3)*0.20948d0
1724           write(*,1100) I,(climat(k),k=1,9)
1725         enddo
1726         write(*,1200) ' O3-column(DU)=',climat(8)/2.687d16
1727       endif
1728       return
1729  1000 format(5X,'Zkm',3X,'Z*',8X,'M',8X,'O3',6X,'f-O3',5X,'T',7X,'P',6x,   &
1730           'col-O3',3X,'col-O2',2X,10(a7,2x))
1731  1100 format(1X,I2,0P,2F6.2,1P,2E10.3,0P,F7.3,F8.2,F10.4,1P,10E9.2)
1732  1200 format(A,F8.1,A,10(1pE10.3))
1733       end SUBROUTINE PRTATM
1734
1735       SUBROUTINE JVALUE(isvode,jsvode,lpar,jpnl, &
1736         ydgrd,sizeaer,extaer,waer,gaer,tauaer,l2,l3,l4,l5,l6,l7)
1737 !-----------------------------------------------------------------------
1738 !  Calculate the actinic flux at each level for the current SZA value.
1739 !        quit when SZA > 98.0 deg ==> tangent height = 63 km
1740 !             or         99.                           80 km
1741 !-----------------------------------------------------------------------
1742 !
1743 !     AVGF   Attenuation of beam at each level for each wavelength
1744 !     FFF    Actinic flux at each desired level
1745 !     FHZ    Actinic flux in Herzberg bin
1746 !     WAVE   Effective wavelength of each wavelength bin
1747 !     XQO2   Absorption cross-section of O2
1748 !     XQO3   Absorption cross-section of O3
1749 !
1750 !-----------------------------------------------------------------------
1751
1752         USE module_data_mosaic_other, only : kmaxd
1753         USE module_fastj_data
1754         USE module_peg_util, only:  peg_message
1755 !jdf
1756 ! Print Fast-J diagnostics if iprint /= 0
1757         integer, parameter :: iprint = 0
1758         integer, parameter :: single = 4        !compiler dependent value real*4
1759 !       integer, parameter :: double = 8        !compiler dependent value real*8
1760        integer,parameter :: ipar_fastj=1,jpar=1
1761 !      integer,parameter :: jppj=14        !Number of photolytic reactions supplied
1762        logical,parameter :: ldeg45=.false. !Logical flag for degraded CTM resolution
1763        integer lpar           !Number of levels in CTM
1764        integer jpnl           !Number of levels requiring chemistry
1765       real(kind=double), dimension(ipar_fastj) :: xgrd !  Longitude (midpoint, radians)
1766       real(kind=double), dimension(jpar) :: ygrd           !  Latitude  (midpoint, radians)
1767       real(kind=double), dimension(jpar) :: ydgrd          !  Latitude  (midpoint, degrees)
1768       real(kind=double), dimension(kmaxd+1) :: etaa    !  Eta(a) value for level boundaries
1769       real(kind=double), dimension(kmaxd+1) :: etab    !  Eta(b) value for level boundaries
1770       real(kind=double) ::  tau_fastj          !  Time of Day (hours, GMT)
1771       integer month_fastj        !  Number of month (1-12)
1772       integer iday_fastj         !  Day of year
1773       real(kind=double), dimension(ipar_fastj,jpar) ::   P , SA
1774       real(kind=double), dimension(ipar_fastj,jpar,kmaxd+1) :: T, OD
1775       integer nspint           ! Num of spectral intervals across solar
1776       parameter ( nspint = 4 ) ! spectrum for FAST-J
1777       real, dimension (nspint),save :: wavmid !cm
1778       real, dimension (nspint, kmaxd+1) :: sizeaer,extaer,waer,gaer,tauaer
1779       real, dimension (nspint, kmaxd+1) :: l2,l3,l4,l5,l6,l7
1780       data wavmid     &
1781         / 0.30e-4, 0.40e-4, 0.60e-4 ,0.999e-04 /
1782 !jdf
1783       integer j, k
1784 !      real*8  wave, xseco3, xseco2
1785       real*8  wave
1786       real*8  AVGF(lpar),XQO3(NB),XQO2(NB)
1787 ! diagnostics for error situations
1788 !      integer lunout
1789 !      parameter (lunout = 41)
1790     integer isvode,jsvode
1791         character*80 msg
1792 !
1793       do J=1,jpnl
1794         do K=NW1,NW2
1795           FFF(K,J) = 0.d0
1796         enddo
1797         FHZ(J) = 0.d0
1798       enddo
1799 !
1800 !---SZA check
1801       if(SZA.gt.szamax) GOTO 99
1802 !
1803 !---Calculate spherical weighting functions
1804       CALL SPHERE(lpar)
1805 !
1806 !---Loop over all wavelength bins
1807       do K=NW1,NW2
1808         WAVE = WL(K)
1809         do J=1,NB
1810           XQO3(J) = xseco3(K,dble(TJ(J)))
1811         enddo
1812         do J=1,NB
1813           XQO2(J) = xseco2(K,dble(TJ(J)))
1814         enddo
1815 !-----------------------------------------
1816         CALL OPMIE(isvode,jsvode,K,WAVE,XQO2,XQO3,AVGF,lpar,jpnl, &
1817           ydgrd,sizeaer,extaer,waer,gaer,tauaer,l2,l3,l4,l5,l6,l7)
1818 !-----------------------------------------
1819         do J=1,jpnl
1820           FFF(K,J) = FFF(K,J) + FL(K)*AVGF(J)
1821 ! diagnostic
1822           if ( FFF(K,J) .lt. 0) then
1823              write( msg, '(a,2i4,e14.6)' )   &
1824                   'FASTJ neg actinic flux ' //   &
1825                   'k j FFF(K,J) ', k, j, fff(k,j)
1826              call peg_message( lunerr, msg )
1827           end if
1828 ! end diagnostic
1829         enddo
1830       enddo
1831 !
1832 !---Herzberg continuum bin above 10 km, if required
1833       if(NHZ.gt.0) then
1834         K=NW2+1
1835         WAVE = 204.d0
1836         do J=1,NB
1837           XQO3(J) = HZO3
1838           XQO2(J) = HZO2
1839         enddo
1840         CALL OPMIE(isvode,jsvode,K,WAVE,XQO2,XQO3,AVGF,lpar,jpnl, &
1841           ydgrd,sizeaer,extaer,waer,gaer,tauaer,l2,l3,l4,l5,l6,l7)
1842         do J=1,jpnl
1843           if(z(j).gt.1.d6) FHZ(J)=AVGF(J)
1844         enddo
1845       endif
1846 !
1847    99 continue
1848  1000 format('  SZA=',f6.1,' Reflectvty=',f6.3,' OD=',10(1pe10.3))
1849
1850       return
1851       end SUBROUTINE JVALUE
1852
1853       FUNCTION xseco3(K,TTT)
1854 !-----------------------------------------------------------------------
1855 !  Cross-sections for O3 for all processes interpolated across 3 temps
1856 !-----------------------------------------------------------------------
1857
1858         USE module_fastj_data
1859
1860       integer k
1861 !      real*8 ttt, flint, xseco3
1862       real*8 ttt, xseco3
1863       xseco3  =   &
1864         flint(TTT,TQQ(1,2),TQQ(2,2),TQQ(3,2),QO3(K,1),QO3(K,2),QO3(K,3))
1865       return
1866       end FUNCTION xseco3
1867
1868       FUNCTION xsec1d(K,TTT)
1869 !-----------------------------------------------------------------------
1870 !  Quantum yields for O3 --> O2 + O(1D) interpolated across 3 temps
1871 !-----------------------------------------------------------------------
1872
1873         USE module_fastj_data
1874
1875       integer k
1876 !      real*8 ttt, flint, xsec1d
1877       real*8 ttt,  xsec1d
1878       xsec1d =   &
1879         flint(TTT,TQQ(1,3),TQQ(2,3),TQQ(3,3),Q1D(K,1),Q1D(K,2),Q1D(K,3))
1880       return
1881       end FUNCTION xsec1d
1882
1883       FUNCTION xseco2(K,TTT)
1884 !-----------------------------------------------------------------------
1885 !  Cross-sections for O2 interpolated across 3 temps; No S_R Bands yet!
1886 !-----------------------------------------------------------------------
1887
1888         USE module_fastj_data
1889
1890       integer k
1891 !      real*8 ttt, flint, xseco2
1892       real*8 ttt,  xseco2
1893       xseco2 =   &
1894         flint(TTT,TQQ(1,1),TQQ(2,1),TQQ(3,1),QO2(K,1),QO2(K,2),QO2(K,3))
1895       return
1896       end FUNCTION xseco2
1897
1898       REAL*8 FUNCTION flint (TINT,T1,T2,T3,F1,F2,F3)
1899 !-----------------------------------------------------------------------
1900 !  Three-point linear interpolation function
1901 !-----------------------------------------------------------------------
1902       real*8 TINT,T1,T2,T3,F1,F2,F3
1903       IF (TINT .LE. T2)  THEN
1904         IF (TINT .LE. T1)  THEN
1905           flint  = F1
1906         ELSE
1907           flint = F1 + (F2 - F1)*(TINT -T1)/(T2 -T1)
1908         ENDIF
1909       ELSE
1910         IF (TINT .GE. T3)  THEN
1911           flint  = F3
1912         ELSE
1913           flint = F2 + (F3 - F2)*(TINT -T2)/(T3 -T2)
1914         ENDIF
1915       ENDIF
1916       return
1917       end FUNCTION flint
1918
1919       SUBROUTINE SPHERE(lpar)
1920 !-----------------------------------------------------------------------
1921 !  Calculation of spherical geometry; derive tangent heights, slant path
1922 !  lengths and air mass factor for each layer. Not called when
1923 !  SZA > 98 degrees.  Beyond 90 degrees, include treatment of emergent
1924 !  beam (where tangent height is below altitude J-value desired at).
1925 !-----------------------------------------------------------------------
1926 !
1927 !     GMU     MU, cos(solar zenith angle)
1928 !     RZ      Distance from centre of Earth to each point (cm)
1929 !     RQ      Square of radius ratios
1930 !     TANHT   Tangent height for the current SZA
1931 !     XL      Slant path between points
1932 !     AMF     Air mass factor for slab between level and level above
1933 !
1934 !-----------------------------------------------------------------------
1935
1936         USE module_fastj_data
1937
1938 !jdf
1939       integer lpar
1940 !jdf
1941       integer i, j, k, ii
1942       real*8 airmas, gmu, xmu1, xmu2, xl, diff
1943       REAL*8 Ux,H,RZ(NB),RQ(NB),ZBYR
1944 !
1945 !  Inlined air mass factor function for top of atmosphere
1946       AIRMAS(Ux,H) = (1.0d0+H)/SQRT(Ux*Ux+2.0d0*H*(1.0d0-   &
1947                0.6817d0*EXP(-57.3d0*ABS(Ux)/SQRT(1.0d0+5500.d0*H))/   &
1948                                                    (1.0d0+0.625d0*H)))
1949 !
1950       GMU = U0
1951       RZ(1)=RAD+Z(1)
1952       ZBYR = ZZHT/RAD
1953       DO 2 II=2,NB
1954         RZ(II) = RAD + Z(II)
1955         RQ(II-1) = (RZ(II-1)/RZ(II))**2
1956     2 CONTINUE
1957       IF (GMU.LT.0.0D0) THEN
1958         TANHT = RZ(nlbatm)/DSQRT(1.0D0-GMU**2)
1959       ELSE
1960         TANHT = RZ(nlbatm)
1961       ENDIF
1962 !
1963 !  Go up from the surface calculating the slant paths between each level
1964 !  and the level above, and deriving the appropriate Air Mass Factor
1965       DO 16 J=1,NB
1966         DO K=1,NB
1967           AMF(K,J)=0.D0
1968         ENDDO
1969 !
1970 !  Air Mass Factors all zero if below the tangent height
1971         IF (RZ(J).LT.TANHT) GOTO 16
1972 !  Ascend from layer J calculating AMFs
1973         XMU1=ABS(GMU)
1974         DO 12 I=J,lpar
1975           XMU2=DSQRT(1.0D0-RQ(I)*(1.0D0-XMU1**2))
1976           XL=RZ(I+1)*XMU2-RZ(I)*XMU1
1977           AMF(I,J)=XL/(RZ(I+1)-RZ(I))
1978           XMU1=XMU2
1979    12   CONTINUE
1980 !  Use function and scale height to provide AMF above top of model
1981         AMF(NB,J)=AIRMAS(XMU1,ZBYR)
1982 !
1983 !  Twilight case - Emergent Beam
1984         IF (GMU.GE.0.0D0) GOTO 16
1985         XMU1=ABS(GMU)
1986 !  Descend from layer J
1987         DO 14 II=J-1,1,-1
1988           DIFF=RZ(II+1)*DSQRT(1.0D0-XMU1**2)-RZ(II)
1989           if(II.eq.1) DIFF=max(DIFF,0.d0)   ! filter
1990 !  Tangent height below current level - beam passes through twice
1991           IF (DIFF.LT.0.0D0) THEN
1992             XMU2=DSQRT(1.0D0-(1.0D0-XMU1**2)/RQ(II))
1993             XL=ABS(RZ(II+1)*XMU1-RZ(II)*XMU2)
1994             AMF(II,J)=2.d0*XL/(RZ(II+1)-RZ(II))
1995             XMU1=XMU2
1996 !  Lowest level intersected by emergent beam
1997           ELSE
1998             XL=RZ(II+1)*XMU1*2.0D0
1999 !            WTING=DIFF/(RZ(II+1)-RZ(II))
2000 !            AMF(II,J)=(1.0D0-WTING)*2.D0**XL/(RZ(II+1)-RZ(II))
2001             AMF(II,J)=XL/(RZ(II+1)-RZ(II))
2002             GOTO 16
2003           ENDIF
2004    14   CONTINUE
2005 !
2006    16 CONTINUE
2007       RETURN
2008       END SUBROUTINE SPHERE
2009
2010
2011       SUBROUTINE OPMIE(isvode,jsvode,KW,WAVEL,XQO2,XQO3,FMEAN,lpar,jpnl, &
2012         ydgrd,sizeaer,extaer,waer,gaer,tauaer,l2,l3,l4,l5,l6,l7)
2013 !-----------------------------------------------------------------------
2014 !  NEW Mie code for J's, only uses 8-term expansion, 4-Gauss pts
2015 !  Currently allow up to NP aerosol phase functions (at all altitudes) to
2016 !  be associated with optical depth AER(1:NC) = aerosol opt.depth @ 1000 nm
2017 !
2018 !  Pick Mie-wavelength with phase function and Qext:
2019 !
2020 !  01 RAYLE = Rayleigh phase
2021 !  02 ISOTR = isotropic
2022 !  03 ABSRB = fully absorbing 'soot', wavelength indep.
2023 !  04 S_Bkg = backgrnd stratospheric sulfate (n=1.46,log-norm:r=.09um/sigma=.6)
2024 !  05 S_Vol = volcanic stratospheric sulfate (n=1.46,log-norm:r=.08um/sigma=.8)
2025 !  06 W_H01 = water haze (H1/Deirm.) (n=1.335, gamma:  r-mode=0.1um /alpha=2)
2026 !  07 W_H04 = water haze (H1/Deirm.) (n=1.335, gamma:  r-mode=0.4um /alpha=2)
2027 !  08 W_C02 = water cloud (C1/Deirm.) (n=1.335, gamma:  r-mode=2.0um /alpha=6)
2028 !  09 W_C04 = water cloud (C1/Deirm.) (n=1.335, gamma:  r-mode=4.0um /alpha=6)
2029 !  10 W_C08 = water cloud (C1/Deirm.) (n=1.335, gamma:  r-mode=8.0um /alpha=6)
2030 !  11 W_C13 = water cloud (C1/Deirm.) (n=1.335, gamma:  r-mode=13.3um /alpha=6)
2031 !  12 W_L06 = water cloud (Lacis) (n=1.335, r-mode=5.5um / alpha=11/3)
2032 !  13 Ice-H = hexagonal ice cloud (Mishchenko)
2033 !  14 Ice-I = irregular ice cloud (Mishchenko)
2034 !
2035 !  Choice of aerosol index MIEDX is made in SET_AER; optical depths are
2036 !  apportioned to the AER array in SET_PROF
2037 !
2038 !-----------------------------------------------------------------------
2039 !  FUNCTION RAYLAY(WAVE)---RAYLEIGH CROSS-SECTION for wave > 170 nm
2040 !       WSQI = 1.E6/(WAVE*WAVE)
2041 !       REFRM1 = 1.0E-6*(64.328+29498.1/(146.-WSQI)+255.4/(41.-WSQI))
2042 !       RAYLAY = 5.40E-21*(REFRM1*WSQI)**2
2043 !-----------------------------------------------------------------------
2044 !
2045 !     DTAUX    Local optical depth of each CTM level
2046 !     PIRAY    Contribution of Rayleigh scattering to extinction
2047 !     PIAER    Contribution of Aerosol scattering to extinction
2048 !     TTAU     Optical depth of air vertically above each point (to top of atm)
2049 !     FTAU     Attenuation of solar beam
2050 !     POMEGA   Scattering phase function
2051 !     FMEAN    Mean actinic flux at desired levels
2052 !
2053 !-----------------------------------------------------------------------
2054
2055         USE module_data_mosaic_other, only : kmaxd
2056         USE module_fastj_data
2057         USE module_peg_util, only:  peg_message, peg_error_fatal
2058 !jdf
2059 ! Print Fast-J diagnostics if iprint /= 0
2060       integer, parameter :: iprint = 0
2061       integer, parameter :: single = 4        !compiler dependent value real*4
2062 !     integer, parameter :: double = 8        !compiler dependent value real*8
2063       integer,parameter :: ipar_fastj=1,jpar=1
2064 !     integer,parameter :: jppj=14        !Number of photolytic reactions supplied
2065       logical,parameter :: ldeg45=.false. !Logical flag for degraded CTM resolution
2066       integer lpar           !Number of levels in CTM
2067       integer jpnl           !Number of levels requiring chemistry
2068       real(kind=double), dimension(ipar_fastj) :: xgrd !  Longitude (midpoint, radians)
2069       real(kind=double), dimension(jpar) :: ygrd           !  Latitude  (midpoint, radians)
2070       real(kind=double), dimension(jpar) :: ydgrd          !  Latitude  (midpoint, degrees)
2071       real(kind=double), dimension(kmaxd+1) :: etaa    !  Eta(a) value for level boundaries
2072       real(kind=double), dimension(kmaxd+1) :: etab    !  Eta(b) value for level boundaries
2073       real(kind=double) ::  tau_fastj          !  Time of Day (hours, GMT)
2074       integer month_fastj        !  Number of month (1-12)
2075       integer iday_fastj         !  Day of year
2076       integer nspint           ! Num of spectral intervals across solar
2077       parameter ( nspint = 4 ) ! spectrum for FAST-J
2078       real, dimension (nspint),save :: wavmid !cm
2079       real, dimension (nspint, kmaxd+1) :: sizeaer,extaer,waer,gaer,tauaer
2080       real, dimension (nspint, kmaxd+1) :: l2,l3,l4,l5,l6,l7
2081       data wavmid     &
2082           / 0.30e-4, 0.40e-4, 0.60e-4 ,0.999e-04 /
2083       INTEGER    NL, N__, M__
2084       PARAMETER (NL=500, N__=2*NL, M__=4)  !wig increased nl from 350 to 500, 31-Oct-2005
2085       REAL(kind=double), dimension(M__) :: A,C1,H,V1,WT,EMU
2086       REAL(kind=double), dimension(M__,M__) :: B,AA,CC,S,W,U1
2087       REAL(kind=double), dimension(M__,2*M__) :: PM
2088       REAL(kind=double), dimension(2*M__) :: PM0
2089       REAL(kind=double), dimension(2*M__,N__) :: POMEGA
2090       REAL(kind=double), dimension(N__) :: ZTAU, FZ, FJ
2091       REAL(kind=double), dimension(M__,M__,N__) :: DD
2092       REAL(kind=double), dimension(M__,N__) :: RR
2093       REAL(kind=double) :: ZREFL,ZFLUX,RADIUS,ZU0
2094       INTEGER ND,N,M,MFIT
2095 !jdf
2096       integer jndlev(lpar),jaddlv(nc),jaddto(nc+1)
2097       integer KW,km,i,j,k,l,ix,j1
2098       integer isvode,jsvode
2099       real*8 QXMIE(MX),XLAER(MX),SSALB(MX)
2100       real*8 xlo2,xlo3,xlray,xltau,zk,taudn,tauup,zk2
2101       real*8 WAVEL,XQO2(NB),XQO3(NB),FMEAN(lpar),POMEGAJ(2*M__,NC+1)
2102       real*8 DTAUX(NB),PIRAY(NB),PIAER(MX,NB),TTAU(NC+1),FTAU(NC+1)
2103       real*8 ftaulog,dttau,dpomega(2*M__)
2104       real*8 ftaulog2,dttau2,dpomega2(2*M__)
2105 ! JJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJ
2106         real*8 PIAER_MX1(NB)
2107 ! BBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB
2108          character*80 msg
2109 !
2110 !---Pick nearest Mie wavelength, no interpolation--------------
2111                               KM=1
2112       if( WAVEL .gt. 355.d0 ) KM=2
2113       if( WAVEL .gt. 500.d0 ) KM=3
2114 !     if( WAVEL .gt. 800.d0 ) KM=4  !drop the 1000 nm wavelength
2115 !
2116 !---For Mie code scale extinction at 1000 nm to wavelength WAVEL (QXMIE)
2117 ! define angstrom exponent
2118 ! JJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJ
2119         ang=log(QAA(1,MIEDX(1))/QAA(4,MIEDX(1)))/log(300./999.)
2120        do I=1,MX
2121 ! QAA is extinction efficiency
2122          QXMIE(I) = QAA(KM,MIEDX(I))/QAA(4,MIEDX(I))
2123 ! scale to 550 nm using angstrom relationship
2124 ! note that this gives QXMIE at 550.0 nm = 1.0, aerosol optical thickness
2125 ! is defined at 550 nm
2126 ! convention -- I = 1 is aerosol, I > 1 are clouds
2127          if(I.eq.1) QXMIE(I) = (WAVEL/550.0)**ang
2128          SSALB(I) = SSA(KM,MIEDX(I)) ! single scattering albedo
2129       enddo
2130 ! BBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB
2131 !
2132 !---Reinitialize arrays
2133       do j=1,nc+1
2134         ttau(j)=0.d0
2135         ftau(j)=0.d0
2136       enddo
2137 !
2138 !---Set up total optical depth over each CTM level, DTAUX
2139       J1 = NLBATM
2140       do J=J1,NB
2141         XLO3=DO3(J)*XQO3(J)
2142         XLO2=DM(J)*XQO2(J)*0.20948d0
2143         XLRAY=DM(J)*QRAYL(KW)
2144 !  Zero absorption for testing purposes
2145 !        call NOABS(XLO3,XLO2,XLRAY,AER(1,j),RFLECT)
2146         do I=1,MX
2147 ! I is aerosol type, j is level, AER(I,J)*QXMIE(I) is the layer aerosol optical thickness
2148 ! at 1000 nm , AER(I,J), times extinction efficiency for the layer (normalized to be one at 1000 nm)
2149 ! therefore xlaer(i) is the layer optical depth at the wavelength index KM
2150           XLAER(I)=AER(I,J)*QXMIE(I)
2151         enddo
2152 !  Total optical depth from all elements
2153         DTAUX(J)=XLO3+XLO2+XLRAY
2154         do I=1,MX
2155           DTAUX(J)=DTAUX(J)+XLAER(I)
2156 !       if(isvode.eq.8.and.jsvode.eq.1) print*,'jdf xlaer',&
2157 !       i,j,km,dtaux(j),xlaer(i),aer(i,j),qxmie(i),xlo3,xlo2,xlray
2158         enddo
2159 ! JJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJ
2160 ! add in new aerosol information from Mie code
2161 ! layer aerosol optical thickness at wavelength index KM, layer j
2162 !       tauaer(km,j)=0.0
2163         dtaux(j)=dtaux(j)+tauaer(km,j)
2164 !       if(isvode.eq.8.and.jsvode.eq.1) print*,'jdf dtaux',&
2165 !       j,km,dtaux(j),tauaer(km,j)
2166 ! BBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB
2167 !  Fractional extinction for Rayleigh scattering and each aerosol type
2168         PIRAY(J)=XLRAY/DTAUX(J)
2169         do I=1,MX
2170           PIAER(I,J)=SSALB(I)*XLAER(I)/DTAUX(J)
2171         enddo
2172 ! JJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJ
2173 ! note the level is now important
2174         PIAER_MX1(J)=waer(km,j)*tauaer(km,j)/DTAUX(J)
2175 ! BBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB
2176       enddo ! of the level "j" loop
2177 !
2178 !---Define the scattering phase fn. with mix of Rayleigh(1) & Mie(MIEDX)
2179 !   No. of quadrature pts fixed at 4 (M__), expansion of phase fn @ 8
2180       N = M__
2181       MFIT = 2*M__
2182       do j=j1,NB  ! jcb: layer index
2183         do i=1,MFIT
2184           pomegaj(i,j) = PIRAY(J)*PAA(I,KM,1) ! jcb: paa are the expansion coefficients of the phase function
2185           do k=1,MX ! jcb: mx is # of aerosols
2186             pomegaj(i,j) = pomegaj(i,j) + PIAER(K,J)*PAA(I,KM,MIEDX(K))
2187           enddo
2188         enddo
2189 ! JJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJ
2190 ! i is the # of coefficients, KM is the wavelength index, j is the level
2191 ! note the level in now relevant because we allow aerosol properties to
2192 ! vary by level
2193            pomegaj(1,j) = pomegaj(1,j) + PIAER_MX1(J)*1.0 ! 1.0 is l0
2194            pomegaj(2,j) = pomegaj(2,j) + PIAER_MX1(J)*gaer(KM,j)*3.0 ! the three converts gear to l1
2195            pomegaj(3,j) = pomegaj(3,j) + PIAER_MX1(J)*l2(KM,j)
2196            pomegaj(4,j) = pomegaj(4,j) + PIAER_MX1(J)*l3(KM,j)
2197            pomegaj(5,j) = pomegaj(5,j) + PIAER_MX1(J)*l4(KM,j)
2198            pomegaj(6,j) = pomegaj(6,j) + PIAER_MX1(J)*l5(KM,j)
2199            pomegaj(7,j) = pomegaj(7,j) + PIAER_MX1(J)*l6(KM,j)
2200            pomegaj(8,j) = pomegaj(7,j) + PIAER_MX1(J)*l7(KM,j)
2201 ! BBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB
2202
2203       enddo
2204 !
2205 !---Calculate attenuated incident beam EXP(-TTAU/U0) and flux on surface
2206       do J=J1,NB
2207         if(AMF(J,J).gt.0.0D0) then
2208           XLTAU=0.0D0
2209           do I=1,NB
2210             XLTAU=XLTAU + DTAUX(I)*AMF(I,J)
2211           enddo
2212           if(XLTAU.gt.450.d0) then   ! for compilers with no underflow trapping
2213             FTAU(j)=0.d0
2214           else
2215             FTAU(J)=DEXP(-XLTAU)
2216           endif
2217         else
2218           FTAU(J)=0.0D0
2219         endif
2220       enddo
2221       if(U0.gt.0.D0) then
2222         ZFLUX = U0*FTAU(J1)*RFLECT/(1.d0+RFLECT)
2223       else
2224         ZFLUX = 0.d0
2225       endif
2226 !
2227 !------------------------------------------------------------------------
2228 !  Take optical properties on CTM layers and convert to a photolysis
2229 !  level grid corresponding to layer centres and boundaries. This is
2230 !  required so that J-values can be calculated for the centre of CTM
2231 !  layers; the index of these layers is kept in the jndlev array.
2232 !------------------------------------------------------------------------
2233 !
2234 !  Set lower boundary and levels to calculate J-values at
2235       J1=2*J1-1
2236       do j=1,lpar
2237         jndlev(j)=2*j
2238       enddo
2239 !
2240 !  Calculate column optical depths above each level, TTAU
2241       TTAU(NC+1)=0.0D0
2242       do J=NC,J1,-1
2243         I=(J+1)/2
2244         TTAU(J)=TTAU(J+1) + 0.5d0*DTAUX(I)
2245         jaddlv(j)=int(0.5d0*DTAUX(I)/dtaumax)
2246 !  Subdivide cloud-top levels if required
2247         if(jadsub(j).gt.0) then
2248           jadsub(j)=min(jaddlv(j)+1,nint(dtausub))*(nint(dsubdiv)-1)
2249           jaddlv(j)=jaddlv(j)+jadsub(j)
2250         endif
2251 !       if(isvode.eq.8.and.jsvode.eq.1) print*,'jdf in opmie',&
2252 !       j,jadsub(j),jaddlv(j),ttau(j),dtaux(i)
2253       enddo
2254 !
2255 !  Calculate attenuated beam, FTAU, level boundaries then level centres
2256       FTAU(NC+1)=1.0D0
2257       do J=NC-1,J1,-2
2258         I=(J+1)/2
2259         FTAU(J)=FTAU(I)
2260       enddo
2261       do J=NC,J1,-2
2262         FTAU(J)=sqrt(FTAU(J+1)*FTAU(J-1))
2263       enddo
2264 !
2265 !  Calculate scattering properties, level centres then level boundaries
2266 !  using an inverse interpolation to give correctly-weighted values
2267       do j=NC,J1,-2
2268         do i=1,MFIT
2269           pomegaj(i,j) = pomegaj(i,j/2)
2270         enddo
2271       enddo
2272       do j=J1+2,nc,2
2273         taudn = ttau(j-1)-ttau(j)
2274         tauup = ttau(j)-ttau(j+1)
2275         do i=1,MFIT
2276           pomegaj(i,j) = (pomegaj(i,j-1)*taudn +   &
2277                           pomegaj(i,j+1)*tauup) / (taudn+tauup)
2278         enddo
2279       enddo
2280 !  Define lower and upper boundaries
2281       do i=1,MFIT
2282         pomegaj(i,J1)   = pomegaj(i,J1+1)
2283         pomegaj(i,nc+1) = pomegaj(i,nc)
2284       enddo
2285 !
2286 !------------------------------------------------------------------------
2287 !  Calculate cumulative total and define levels we want J-values at.
2288 !  Sum upwards for levels, and then downwards for Mie code readjustments.
2289 !
2290 !     jaddlv(i)   Number of new levels to add between (i) and (i+1)
2291 !     jaddto(i)   Total number of new levels to add to and above level (i)
2292 !     jndlev(j)   Level needed for J-value for CTM layer (j)
2293 !
2294 !------------------------------------------------------------------------
2295 !
2296 !  Reinitialize level arrays
2297       do j=1,nc+1
2298         jaddto(j)=0
2299       enddo
2300 !
2301       jaddto(J1)=jaddlv(J1)
2302       do j=J1+1,nc
2303         jaddto(j)=jaddto(j-1)+jaddlv(j)
2304       enddo
2305       if((jaddto(nc)+nc).gt.nl) then
2306 !         print*,'jdf mie',isvode,jsvode,jaddto(nc),nc,nl
2307          write ( msg, '(a, 2i6)' )  &
2308            'FASTJ  Max NL exceeded '  //  &
2309            'jaddto(nc)+nc NL', jaddto(nc)+nc,NL
2310          call peg_message( lunerr, msg )
2311          msg = 'FASTJ subr OPMIE error. Max NL exceeded'
2312          call peg_error_fatal( lunerr, msg )
2313 !         write(6,1500)  jaddto(nc)+nc, 'NL',NL
2314 !         stop
2315       endif
2316       do i=1,lpar
2317         jndlev(i)=jndlev(i)+jaddto(jndlev(i)-1)
2318       enddo
2319       jaddto(nc)=jaddlv(nc)
2320       do j=nc-1,J1,-1
2321         jaddto(j)=jaddto(j+1)+jaddlv(j)
2322       enddo
2323 !
2324 !---------------------SET UP FOR MIE CODE-------------------------------
2325 !
2326 !  Transpose the ascending TTAU grid to a descending ZTAU grid.
2327 !  Double the resolution - TTAU points become the odd points on the
2328 !  ZTAU grid, even points needed for asymm phase fn soln, contain 'h'.
2329 !  Odd point added at top of grid for unattenuated beam   (Z='inf')
2330 !
2331 !        Surface:   TTAU(1)   now use ZTAU(2*NC+1)
2332 !        Top:       TTAU(NC)  now use ZTAU(3)
2333 !        Infinity:            now use ZTAU(1)
2334 !
2335 !  Mie scattering code only used from surface to level NC
2336 !------------------------------------------------------------------------
2337 !
2338 !  Initialise all Fast-J optical property arrays
2339       do k=1,N__
2340         do i=1,MFIT
2341           pomega(i,k) = 0.d0
2342         enddo
2343         ztau(k) = 0.d0
2344         fz(k)   = 0.d0
2345       enddo
2346 !
2347 !  Ascend through atmosphere transposing grid and adding extra points
2348       do j=J1,nc+1
2349         k = 2*(nc+1-j)+2*jaddto(j)+1
2350         ztau(k)= ttau(j)
2351         fz(k)  = ftau(j)
2352         do i=1,MFIT
2353           pomega(i,k) = pomegaj(i,j)
2354         enddo
2355       enddo
2356 !
2357 !  Check profiles if desired
2358 !      ND = 2*(NC+jaddto(J1)-J1)  + 3
2359 !      if(kw.eq.1) call CH_PROF
2360 !
2361 !------------------------------------------------------------------------
2362 !    Insert new levels, working downwards from the top of the atmosphere
2363 !  to the surface (down in 'j', up in 'k'). This allows ztau and pomega
2364 !  to be incremented linearly (in a +ve sense), and the flux fz to be
2365 !  attenuated top-down (avoiding problems where lower level fluxes are
2366 !  zero).
2367 !
2368 !    zk        fractional increment in level
2369 !    dttau     change in ttau per increment    (linear, positive)
2370 !    dpomega   change in pomega per increment  (linear)
2371 !    ftaulog   change in ftau per increment    (exponential, normally < 1)
2372 !
2373 !------------------------------------------------------------------------
2374 !
2375       do j=nc,J1,-1
2376           zk = 0.5d0/(1.d0+dble(jaddlv(j)-jadsub(j)))
2377           dttau = (ttau(j)-ttau(j+1))*zk
2378           do i=1,MFIT
2379             dpomega(i) = (pomegaj(i,j)-pomegaj(i,j+1))*zk
2380           enddo
2381 !  Filter attenuation factor - set minimum at 1.0d-05
2382           if(ftau(j+1).eq.0.d0) then
2383             ftaulog=0.d0
2384           else
2385             ftaulog = ftau(j)/ftau(j+1)
2386             if(ftaulog.lt.1.d-150) then
2387               ftaulog=1.0d-05
2388             else
2389               ftaulog=exp(log(ftaulog)*zk)
2390             endif
2391           endif
2392           k = 2*(nc-j+jaddto(j)-jaddlv(j))+1   !  k at level j+1
2393           l = 0
2394 !  Additional subdivision of first level if required
2395           if(jadsub(j).ne.0) then
2396             l=jadsub(j)/nint(dsubdiv-1)
2397             zk2=1.d0/dsubdiv
2398             dttau2=dttau*zk2
2399             ftaulog2=ftaulog**zk2
2400             do i=1,MFIT
2401               dpomega2(i)=dpomega(i)*zk2
2402             enddo
2403             do ix=1,2*(jadsub(j)+l)
2404               ztau(k+1) = ztau(k) + dttau2
2405               fz(k+1) = fz(k)*ftaulog2
2406               do i=1,MFIT
2407                 pomega(i,k+1) = pomega(i,k) + dpomega2(i)
2408               enddo
2409               k = k+1
2410             enddo
2411           endif
2412           l = 2*(jaddlv(j)-jadsub(j)-l)+1
2413 !
2414 !  Add values at all intermediate levels
2415           do ix=1,l
2416             ztau(k+1) = ztau(k) + dttau
2417             fz(k+1) = fz(k)*ftaulog
2418             do i=1,MFIT
2419               pomega(i,k+1) = pomega(i,k) + dpomega(i)
2420             enddo
2421             k = k+1
2422           enddo
2423 !
2424 !  Alternate method to attenuate fluxes, fz, using 2nd-order finite
2425 !  difference scheme - just need to comment in section below
2426 !          ix = 2*(jaddlv(j)-jadsub(j))+1
2427 !          if(l.le.0) then
2428 !            l=k-ix-1
2429 !          else
2430 !            l=k-ix
2431 !          endif
2432 !          call efold(ftau(j+1),ftau(j),ix+1,fz(l))
2433 !          if(jadsub(j).ne.0) then
2434 !            k = 2*(nc-j+jaddto(j)-jaddlv(j))+1 !  k at level j+1
2435 !            ix=2*(jadsub(j)+(jadsub(j)/nint(dsubdiv-1)))
2436 !            call efold(ftau(j+1),fz(k+ix),ix,fz(k))
2437 !          endif
2438 !
2439       enddo
2440 !
2441 !---Update total number of levels and check doesn't exceed N__
2442       ND = 2*(NC+jaddto(J1)-J1)  + 3
2443       if(nd.gt.N__) then
2444          write ( msg, '(a, 2i6)' )  &
2445            'FASTJ  Max N__ exceeded '  //  &
2446            'ND N__', ND, N__
2447          call peg_message( lunerr, msg )
2448          msg = 'FASTJ subr OPMIE error. Max N__ exceeded'
2449          call peg_error_fatal( lunerr, msg )
2450 !        write(6,1500) ND, 'N__',N__
2451 !        stop
2452       endif
2453 !
2454 !---Add boundary/ground layer to ensure no negative J's caused by
2455 !---too large a TTAU-step in the 2nd-order lower b.c.
2456       ZTAU(ND+1) = ZTAU(ND)*1.000005d0
2457       ZTAU(ND+2) = ZTAU(ND)*1.000010d0
2458       zk=max(abs(U0),0.01d0)
2459       zk=dexp(-ZTAU(ND)*5.d-6/zk)
2460       FZ(ND+1) = FZ(ND)*zk
2461       FZ(ND+2) = FZ(ND+1)*zk
2462       do I=1,MFIT
2463         POMEGA(I,ND+1)   = POMEGA(I,ND)
2464         POMEGA(I,ND+2)   = POMEGA(I,ND)
2465       enddo
2466       ND = ND+2
2467 !
2468       ZU0 = U0
2469       ZREFL = RFLECT
2470 !
2471 !-----------------------------------------
2472       CALL MIESCT(ND,N,M,MFIT,POMEGA,PM,PM0,FZ,WT,EMU,ZTAU,ZFLUX, &
2473         ZREFL,FJ,A,C1,H,V1,B,AA,CC,S,W,U1,DD,RR,RADIUS,ZU0)
2474 !-----------------------------------------
2475 !  Accumulate attenuation for selected levels
2476       l=2*(NC+jaddto(J1))+3
2477       do j=1,lpar
2478         k=l-(2*jndlev(j))
2479         if(k.gt.ND-2) then
2480           FMEAN(j) = 0.d0
2481         else
2482           FMEAN(j) = FJ(k)
2483         endif
2484       enddo
2485 !
2486       return
2487  1000 format(1x,i3,3(2x,1pe10.4),1x,i3)
2488  1300 format(1x,50(i3))
2489  1500 format(' Too many levels in photolysis code: need ',i3,' but ',a,   &
2490              ' dimensioned as ',i3)
2491       END SUBROUTINE OPMIE
2492
2493 !*********************************************************************
2494       subroutine EFOLD (F0, F1, N, F)
2495 !-----------------------------------------------------------------------
2496 !---  calculate the e-fold between two boundaries, given the value
2497 !---     at both boundaries F0(x=0) = top, F1(x=1) = bottom.
2498 !---  presume that F(x) proportional to exp[-A*x] for x=0 to x=1
2499 !---          d2F/dx2 = A*A*F  and thus expect F1 = F0 * exp[-A]
2500 !---           alternatively, could define A = ln[F0/F1]
2501 !---  let X = A*x, d2F/dX2 = F
2502 !---  assume equal spacing (not necessary, but makes this easier)
2503 !---      with N-1 intermediate points (and N layers of thickness dX = A/N)
2504 !---
2505 !---  2nd-order finite difference:  (F(i-1) - 2F(i) + F(i+1)) / dX*dX = F(i)
2506 !---      let D = 1 / dX*dX:
2507 !
2508 !  1  |   1        0        0        0        0        0   |    | F0 |
2509 !     |                                                    |    | 0  |
2510 !  2  |  -D      2D+1      -D        0        0        0   |    | 0  |
2511 !     |                                                    |    | 0  |
2512 !  3  |   0       -D      2D+1      -D        0        0   |    | 0  |
2513 !     |                                                    |    | 0  |
2514 !     |   0        0       -D      2D+1      -D        0   |    | 0  |
2515 !     |                                                    |    | 0  |
2516 !  N  |   0        0        0       -D      2D+1      -D   |    | 0  |
2517 !     |                                                    |    | 0  |
2518 ! N+1 |   0        0        0        0        0        1   |    | F1 |
2519 !
2520 !-----------------------------------------------------------------------
2521 !  Advantage of scheme over simple attenuation factor: conserves total
2522 !  number of photons - very useful when using scheme for heating rates.
2523 !  Disadvantage: although reproduces e-folds very well for small flux
2524 !  differences, starts to drift off when many orders of magnitude are
2525 !  involved.
2526 !-----------------------------------------------------------------------
2527       implicit none
2528       real*8 F0,F1,F(250)  !F(N+1)
2529       integer N
2530       integer I
2531       real*8 A,DX,D,DSQ,DDP1, B(101),R(101)
2532 !
2533       if(F0.eq.0.d0) then
2534         do I=1,N
2535           F(I)=0.d0
2536         enddo
2537         return
2538       elseif(F1.eq.0.d0) then
2539         A = DLOG(F0/1.d-250)
2540       else
2541         A = DLOG(F0/F1)
2542       endif
2543 !
2544       DX = float(N)/A
2545       D = DX*DX
2546       DSQ = D*D
2547       DDP1 = D+D+1.d0
2548 !
2549       B(2) = DDP1
2550       R(2) = +D*F0
2551       do I=3,N
2552         B(I) = DDP1 - DSQ/B(I-1)
2553         R(I) = +D*R(I-1)/B(I-1)
2554       enddo
2555       F(N+1) = F1
2556       do I=N,2,-1
2557         F(I) = (R(I) + D*F(I+1))/B(I)
2558       enddo
2559       F(1) = F0
2560       return
2561       end subroutine EFOLD
2562
2563
2564       subroutine CH_PROF(ND,N,M,MFIT,POMEGA,PM,PM0,FZ,WT,EMU,ZTAU, &
2565         ZFLUX,ZREFL,FJ,A,C1,H,V1,B,AA,CC,S,W,U1,DD,RR,RADIUS,ZU0)
2566 !-----------------------------------------------------------------------
2567 !  Check profiles to be passed to MIESCT
2568 !-----------------------------------------------------------------------
2569
2570       USE module_peg_util, only:  peg_message
2571
2572       implicit none
2573 !jdf
2574       integer, parameter :: single = 4      !compiler dependent value real*4
2575       integer, parameter :: double = 8      !compiler dependent value real*8
2576       INTEGER    NL, N__, M__
2577       PARAMETER (NL=350, N__=2*NL, M__=4)
2578       REAL(kind=double), dimension(M__) :: A,C1,H,V1,WT,EMU
2579       REAL(kind=double), dimension(M__,M__) :: B,AA,CC,S,W,U1
2580       REAL(kind=double), dimension(M__,2*M__) :: PM
2581       REAL(kind=double), dimension(2*M__) :: PM0
2582       REAL(kind=double), dimension(2*M__,N__) :: POMEGA
2583       REAL(kind=double), dimension(N__) :: ZTAU, FZ, FJ
2584       REAL(kind=double), dimension(M__,M__,N__) :: DD
2585       REAL(kind=double), dimension(M__,N__) :: RR
2586       REAL(kind=double) :: ZREFL,ZFLUX,RADIUS,ZU0
2587       INTEGER ND,N,M,MFIT
2588 !jdf
2589       integer i,j
2590       character*80 msg
2591 !      write(6,1100) 'lev','ztau','fz  ','pomega( )'
2592       do i=1,ND
2593         if(ztau(i).ne.0.d0) then
2594           write ( msg, '(a, i3, 2(1x,1pe9.3))' )        &
2595               'FASTJ subr CH_PROF ztau ne 0. check pomega. ' //  &
2596               'k ztau(k) fz(k) ', i,ztau(i),fz(i)
2597           call peg_message( lunerr, msg )
2598 !          write(6,1200) i,ztau(i),fz(i),(pomega(j,i),j=1,8)
2599         endif
2600       enddo
2601       return
2602  1100 format(1x,a3,4(a9,2x))
2603  1200 format(1x,i3,11(1x,1pe9.3))
2604       end subroutine CH_PROF
2605
2606
2607       SUBROUTINE MIESCT(ND,N,M,MFIT,POMEGA,PM,PM0,FZ,WT,EMU,ZTAU, &
2608         ZFLUX,ZREFL,FJ,A,C1,H,V1,B,AA,CC,S,W,U1,DD,RR,RADIUS,ZU0)
2609 !     SUBROUTINE MIESCT
2610 !-----------------------------------------------------------------------
2611 !   This is an adaption of the Prather radiative transfer code, (mjp, 10/95)
2612 !     Prather, 1974, Astrophys. J. 192, 787-792.
2613 !         Sol'n of inhomogeneous Rayleigh scattering atmosphere.
2614 !         (original Rayleigh w/ polarization)
2615 !     Cochran and Trafton, 1978, Ap.J., 219, 756-762.
2616 !         Raman scattering in the atmospheres of the major planets.
2617 !         (first use of anisotropic code)
2618 !     Jacob, Gottlieb and Prather, 1989, J.Geophys.Res., 94, 12975-13002.
2619 !         Chemistry of a polluted cloudy boundary layer,
2620 !         (documentation of extension to anisotropic scattering)
2621 !
2622 !    takes atmospheric structure and source terms from std J-code
2623 !    ALSO limited to 4 Gauss points, only calculates mean field!
2624 !
2625 !   mean rad. field ONLY (M=1)
2626 !   initialize variables FIXED/UNUSED in this special version:
2627 !   FTOP = 1.0 = astrophysical flux (unit of pi) at SZA, -ZU0, use for scaling
2628 !   FBOT = 0.0 = external isotropic flux on lower boundary
2629 !   SISOTP = 0.0 = Specific Intensity of isotropic radiation incident from top
2630 !
2631 !   SUBROUTINES:  MIESCT              needs 'jv_mie.cmn'
2632 !                 BLKSLV              needs 'jv_mie.cmn'
2633 !                 GEN (ID)            needs 'jv_mie.cmn'
2634 !                 LEGND0 (X,PL,N)
2635 !                 MATIN4 (A)
2636 !                 GAUSSP (N,XPT,XWT)
2637 !-----------------------------------------------------------------------
2638
2639 !      INCLUDE 'jv_mie.h'
2640
2641        implicit none
2642 !jdf
2643       integer, parameter :: single = 4      !compiler dependent value real*4
2644       integer, parameter :: double = 8      !compiler dependent value real*8
2645       INTEGER    NL, N__, M__
2646       PARAMETER (NL=350, N__=2*NL, M__=4)
2647       REAL(kind=double), dimension(M__) :: A,C1,H,V1,WT,EMU
2648       REAL(kind=double), dimension(M__,M__) :: B,AA,CC,S,W,U1
2649       REAL(kind=double), dimension(M__,2*M__) :: PM
2650       REAL(kind=double), dimension(2*M__) :: PM0
2651       REAL(kind=double), dimension(2*M__,N__) :: POMEGA
2652       REAL(kind=double), dimension(N__) :: ZTAU, FZ, FJ
2653       REAL(kind=double), dimension(M__,M__,N__) :: DD
2654       REAL(kind=double), dimension(M__,N__) :: RR
2655       REAL(kind=double) :: ZREFL,ZFLUX,RADIUS,ZU0
2656       INTEGER ND,N,M,MFIT
2657 !jdf
2658       integer i, id, im
2659       real*8  cmeq1
2660 !-----------------------------------------------------------------------
2661 !---fix scattering to 4 Gauss pts = 8-stream
2662       CALL GAUSSP (N,EMU,WT)
2663 !---solve eqn of R.T. only for first-order M=1
2664 !      ZFLUX = (ZU0*FZ(ND)*ZREFL+FBOT)/(1.0d0+ZREFL)
2665       ZFLUX = (ZU0*FZ(ND)*ZREFL)/(1.0d0+ZREFL)
2666       M=1
2667       DO I=1,N
2668         CALL LEGND0 (EMU(I),PM0,MFIT)
2669         DO IM=M,MFIT
2670           PM(I,IM) = PM0(IM)
2671         ENDDO
2672       ENDDO
2673 !
2674       CMEQ1 = 0.25D0
2675       CALL LEGND0 (-ZU0,PM0,MFIT)
2676       DO IM=M,MFIT
2677         PM0(IM) = CMEQ1*PM0(IM)
2678       ENDDO
2679 !
2680 !     CALL BLKSLV
2681       CALL BLKSLV(ND,N,M,MFIT,POMEGA,PM,PM0,FZ,WT,EMU,ZTAU,ZFLUX, &
2682         ZREFL,FJ,A,C1,H,V1,B,AA,CC,S,W,U1,DD,RR,RADIUS,ZU0)
2683 !
2684       DO ID=1,ND,2
2685         FJ(ID) = 4.0d0*FJ(ID) + FZ(ID)
2686       ENDDO
2687
2688       RETURN
2689       END SUBROUTINE MIESCT
2690
2691       SUBROUTINE BLKSLV(ND,N,M,MFIT,POMEGA,PM,PM0,FZ,WT,EMU,ZTAU, &
2692         ZFLUX,ZREFL,FJ,A,C1,H,V1,B,AA,CC,S,W,U1,DD,RR,RADIUS,ZU0)
2693 !-----------------------------------------------------------------------
2694 !  Solves the block tri-diagonal system:
2695 !               A(I)*X(I-1) + B(I)*X(I) + C(I)*X(I+1) = H(I)
2696 !-----------------------------------------------------------------------
2697 !      INCLUDE 'jv_mie.h'
2698
2699         implicit none
2700 !jdf
2701       integer, parameter :: single = 4      !compiler dependent value real*4
2702       integer, parameter :: double = 8      !compiler dependent value real*8
2703       INTEGER    NL, N__, M__
2704       PARAMETER (NL=350, N__=2*NL, M__=4)
2705       REAL(kind=double), dimension(M__) :: A,C1,H,V1,WT,EMU
2706       REAL(kind=double), dimension(M__,M__) :: B,AA,CC,S,W,U1
2707       REAL(kind=double), dimension(M__,2*M__) :: PM
2708       REAL(kind=double), dimension(2*M__) :: PM0
2709       REAL(kind=double), dimension(2*M__,N__) :: POMEGA
2710       REAL(kind=double), dimension(N__) :: ZTAU, FZ, FJ
2711       REAL(kind=double), dimension(M__,M__,N__) :: DD
2712       REAL(kind=double), dimension(M__,N__) :: RR
2713       REAL(kind=double) :: ZREFL,ZFLUX,RADIUS,ZU0
2714       INTEGER ND,N,M,MFIT
2715 !jdf
2716       integer i, j, k, id
2717       real*8  thesum
2718 !-----------UPPER BOUNDARY ID=1
2719       CALL GEN(1,ND,N,M,MFIT,POMEGA,PM,PM0,FZ,WT,EMU,ZTAU,ZFLUX, &
2720         ZREFL,FJ,A,C1,H,V1,B,AA,CC,S,W,U1,DD,RR,RADIUS,ZU0)
2721       CALL MATIN4 (B)
2722       DO I=1,N
2723          RR(I,1) = 0.0d0
2724         DO J=1,N
2725           THESUM = 0.0d0
2726          DO K=1,N
2727           THESUM = THESUM - B(I,K)*CC(K,J)
2728          ENDDO
2729          DD(I,J,1) = THESUM
2730          RR(I,1) = RR(I,1) + B(I,J)*H(J)
2731         ENDDO
2732       ENDDO
2733 !----------CONTINUE THROUGH ALL DEPTH POINTS ID=2 TO ID=ND-1
2734       DO ID=2,ND-1
2735         CALL GEN(ID,ND,N,M,MFIT,POMEGA,PM,PM0,FZ,WT,EMU,ZTAU,ZFLUX, &
2736           ZREFL,FJ,A,C1,H,V1,B,AA,CC,S,W,U1,DD,RR,RADIUS,ZU0)
2737         DO I=1,N
2738           DO J=1,N
2739           B(I,J) = B(I,J) + A(I)*DD(I,J,ID-1)
2740           ENDDO
2741           H(I) = H(I) - A(I)*RR(I,ID-1)
2742         ENDDO
2743         CALL MATIN4 (B)
2744         DO I=1,N
2745           RR(I,ID) = 0.0d0
2746           DO J=1,N
2747           RR(I,ID) = RR(I,ID) + B(I,J)*H(J)
2748           DD(I,J,ID) = - B(I,J)*C1(J)
2749           ENDDO
2750         ENDDO
2751       ENDDO
2752 !---------FINAL DEPTH POINT: ND
2753       CALL GEN(ND,ND,N,M,MFIT,POMEGA,PM,PM0,FZ,WT,EMU,ZTAU,ZFLUX, &
2754         ZREFL,FJ,A,C1,H,V1,B,AA,CC,S,W,U1,DD,RR,RADIUS,ZU0)
2755       DO I=1,N
2756         DO J=1,N
2757           THESUM = 0.0d0
2758           DO K=1,N
2759           THESUM = THESUM + AA(I,K)*DD(K,J,ND-1)
2760           ENDDO
2761         B(I,J) = B(I,J) + THESUM
2762         H(I) = H(I) - AA(I,J)*RR(J,ND-1)
2763         ENDDO
2764       ENDDO
2765       CALL MATIN4 (B)
2766       DO I=1,N
2767         RR(I,ND) = 0.0d0
2768         DO J=1,N
2769         RR(I,ND) = RR(I,ND) + B(I,J)*H(J)
2770         ENDDO
2771       ENDDO
2772 !-----------BACK SOLUTION
2773       DO ID=ND-1,1,-1
2774        DO I=1,N
2775         DO J=1,N
2776          RR(I,ID) = RR(I,ID) + DD(I,J,ID)*RR(J,ID+1)
2777         ENDDO
2778        ENDDO
2779       ENDDO
2780 !----------MEAN J & H
2781       DO ID=1,ND,2
2782         FJ(ID) = 0.0d0
2783        DO I=1,N
2784         FJ(ID) = FJ(ID) + RR(I,ID)*WT(I)
2785        ENDDO
2786       ENDDO
2787       DO ID=2,ND,2
2788         FJ(ID) = 0.0d0
2789        DO I=1,N
2790         FJ(ID) = FJ(ID) + RR(I,ID)*WT(I)*EMU(I)
2791        ENDDO
2792       ENDDO
2793 ! Output fluxes for testing purposes
2794 !      CALL CH_FLUX
2795 !      CALL CH_FLUX(ND,N,M,MFIT,POMEGA,PM,PM0,FZ,WT,EMU,ZTAU, &
2796 !       ZFLUX,ZREFL,FJ,A,C1,H,V1,B,AA,CC,S,W,U1,DD,RR,RADIUS,ZU0)
2797 !
2798       RETURN
2799       END SUBROUTINE BLKSLV
2800
2801
2802       SUBROUTINE CH_FLUX(ND,N,M,MFIT,POMEGA,PM,PM0,FZ,WT,EMU,ZTAU, &
2803         ZFLUX,ZREFL,FJ,A,C1,H,V1,B,AA,CC,S,W,U1,DD,RR,RADIUS,ZU0)
2804 !-----------------------------------------------------------------------
2805 !  Diagnostic routine to check fluxes at each level - makes most sense
2806 !  when running a conservative atmosphere (zero out absorption in
2807 !  OPMIE by calling the NOABS routine below)
2808 !-----------------------------------------------------------------------
2809
2810 !      INCLUDE 'jv_mie.h'
2811         IMPLICIT NONE
2812 !jdf
2813       integer, parameter :: single = 4      !compiler dependent value real*4
2814       integer, parameter :: double = 8      !compiler dependent value real*8
2815       INTEGER    NL, N__, M__
2816       PARAMETER (NL=350, N__=2*NL, M__=4)
2817       REAL(kind=double), dimension(M__) :: A,C1,H,V1,WT,EMU
2818       REAL(kind=double), dimension(M__,M__) :: B,AA,CC,S,W,U1
2819       REAL(kind=double), dimension(M__,2*M__) :: PM
2820       REAL(kind=double), dimension(2*M__) :: PM0
2821       REAL(kind=double), dimension(2*M__,N__) :: POMEGA
2822       REAL(kind=double), dimension(N__) :: ZTAU, FZ, FJ
2823       REAL(kind=double), dimension(M__,M__,N__) :: DD
2824       REAL(kind=double), dimension(M__,N__) :: RR
2825       REAL(kind=double) :: ZREFL,ZFLUX,RADIUS,ZU0
2826       INTEGER ND,N,M,MFIT
2827 !jdf
2828       integer I,ID
2829       real*8 FJCHEK(N__),FZMEAN
2830 !
2831 !  Odd (h) levels held as actinic flux, so recalculate irradiances
2832       DO ID=1,ND,2
2833         FJCHEK(ID) = 0.0d0
2834        DO I=1,N
2835         FJCHEK(ID) = FJCHEK(ID) + RR(I,ID)*WT(I)*EMU(i)
2836        ENDDO
2837       ENDDO
2838 !
2839 !  Even (j) levels are already held as irradiances
2840       DO ID=2,ND,2
2841        DO I=1,N
2842         FJCHEK(ID) = FJ(ID)
2843        ENDDO
2844       ENDDO
2845 !
2846 !  Output Downward and Upward fluxes down through atmosphere
2847 !     WRITE(6,1200)
2848       WRITE(34,1200)
2849       DO ID=2,ND,2
2850         FZMEAN=sqrt(FZ(ID)*FZ(ID-1))
2851 !       WRITE(6,1000) ID, ZU0*FZMEAN-2.0*(FJCHEK(id)-FJCHEK(id-1)),   &
2852         WRITE(34,1000) ID, ZU0*FZMEAN-2.0*(FJCHEK(id)-FJCHEK(id-1)),   &
2853                                      2.0*(FJCHEK(id)+FJCHEK(id-1)),   &
2854                                      2.0*(FJCHEK(id)+FJCHEK(id-1))/   &
2855                          (ZU0*FZMEAN-2.0*(FJCHEK(id)-FJCHEK(id-1)))
2856       ENDDO
2857       RETURN
2858  1000 FORMAT(1x,i3,1p,2E12.4,1x,0p,f9.4)
2859  1200 FORMAT(1x,'Lev',3x,'Downward',4x,'Upward',7x,'Ratio')
2860       END SUBROUTINE CH_FLUX
2861
2862       SUBROUTINE NOABS(XLO3,XLO2,XLRAY,BCAER,RFLECT)
2863 !-----------------------------------------------------------------------
2864 !  Zero out absorption terms to check scattering code. Leave a little
2865 !  Rayleigh to provide a minimal optical depth, and set surface albedo
2866 !  to unity.
2867 !-----------------------------------------------------------------------
2868       IMPLICIT NONE
2869       real*8 XLO3,XLO2,XLRAY,BCAER,RFLECT
2870       XLO3=0.d0
2871       XLO2=0.d0
2872       XLRAY=XLRAY*1.d-10
2873       BCAER=0.d0
2874       RFLECT=1.d0
2875       RETURN
2876       END SUBROUTINE NOABS
2877
2878       SUBROUTINE GEN(ID,ND,N,M,MFIT,POMEGA,PM,PM0,FZ,WT,EMU,ZTAU, &
2879         ZFLUX,ZREFL,FJ,A,C1,H,V1,B,AA,CC,S,W,U1,DD,RR,RADIUS,ZU0)
2880 !-----------------------------------------------------------------------
2881 !  Generates coefficient matrices for the block tri-diagonal system:
2882 !               A(I)*X(I-1) + B(I)*X(I) + C(I)*X(I+1) = H(I)
2883 !-----------------------------------------------------------------------
2884
2885 !      INCLUDE 'jv_mie.h'
2886         IMPLICIT NONE
2887 !jdf
2888       integer, parameter :: single = 4      !compiler dependent value real*4
2889       integer, parameter :: double = 8      !compiler dependent value real*8
2890       INTEGER    NL, N__, M__
2891       PARAMETER (NL=350, N__=2*NL, M__=4)
2892       REAL(kind=double), dimension(M__) :: A,C1,H,V1,WT,EMU
2893       REAL(kind=double), dimension(M__,M__) :: B,AA,CC,S,W,U1
2894       REAL(kind=double), dimension(M__,2*M__) :: PM
2895       REAL(kind=double), dimension(2*M__) :: PM0
2896       REAL(kind=double), dimension(2*M__,N__) :: POMEGA
2897       REAL(kind=double), dimension(N__) :: ZTAU, FZ, FJ
2898       REAL(kind=double), dimension(M__,M__,N__) :: DD
2899       REAL(kind=double), dimension(M__,N__) :: RR
2900       REAL(kind=double) :: ZREFL,ZFLUX,RADIUS,ZU0
2901       INTEGER ND,N,M,MFIT
2902 !jdf
2903       integer id, id0, id1, im, i, j, k, mstart
2904       real*8  sum0, sum1, sum2, sum3
2905       real*8  deltau, d1, d2, surfac
2906 !---------------------------------------------
2907       IF(ID.EQ.1 .OR. ID.EQ.ND) THEN
2908 !---------calculate generic 2nd-order terms for boundaries
2909        ID0 = ID
2910        ID1 = ID+1
2911        IF(ID.GE.ND) ID1 = ID-1
2912        DO 10 I=1,N
2913           SUM0 = 0.0d0
2914           SUM1 = 0.0d0
2915           SUM2 = 0.0d0
2916           SUM3 = 0.0d0
2917         DO IM=M,MFIT,2
2918           SUM0 = SUM0 + POMEGA(IM,ID0)*PM(I,IM)*PM0(IM)
2919           SUM2 = SUM2 + POMEGA(IM,ID1)*PM(I,IM)*PM0(IM)
2920         ENDDO
2921         DO IM=M+1,MFIT,2
2922           SUM1 = SUM1 + POMEGA(IM,ID0)*PM(I,IM)*PM0(IM)
2923           SUM3 = SUM3 + POMEGA(IM,ID1)*PM(I,IM)*PM0(IM)
2924         ENDDO
2925          H(I) = 0.5d0*(SUM0*FZ(ID0) + SUM2*FZ(ID1))
2926          A(I) = 0.5d0*(SUM1*FZ(ID0) + SUM3*FZ(ID1))
2927         DO J=1,I
2928           SUM0 = 0.0d0
2929           SUM1 = 0.0d0
2930           SUM2 = 0.0d0
2931           SUM3 = 0.0d0
2932          DO IM=M,MFIT,2
2933           SUM0 = SUM0 + POMEGA(IM,ID0)*PM(I,IM)*PM(J,IM)
2934           SUM2 = SUM2 + POMEGA(IM,ID1)*PM(I,IM)*PM(J,IM)
2935          ENDDO
2936          DO IM=M+1,MFIT,2
2937           SUM1 = SUM1 + POMEGA(IM,ID0)*PM(I,IM)*PM(J,IM)
2938           SUM3 = SUM3 + POMEGA(IM,ID1)*PM(I,IM)*PM(J,IM)
2939          ENDDO
2940          S(I,J) = - SUM2*WT(J)
2941          S(J,I) = - SUM2*WT(I)
2942          W(I,J) = - SUM1*WT(J)
2943          W(J,I) = - SUM1*WT(I)
2944          U1(I,J) = - SUM3*WT(J)
2945          U1(J,I) = - SUM3*WT(I)
2946           SUM0 = 0.5d0*(SUM0 + SUM2)
2947          B(I,J) = - SUM0*WT(J)
2948          B(J,I) = - SUM0*WT(I)
2949         ENDDO
2950          S(I,I) = S(I,I) + 1.0d0
2951          W(I,I) = W(I,I) + 1.0d0
2952          U1(I,I) = U1(I,I) + 1.0d0
2953          B(I,I) = B(I,I) + 1.0d0
2954    10  CONTINUE
2955        DO I=1,N
2956          SUM0 = 0.0d0
2957         DO J=1,N
2958          SUM0 = SUM0 + S(I,J)*A(J)/EMU(J)
2959         ENDDO
2960         C1(I) = SUM0
2961        ENDDO
2962        DO I=1,N
2963         DO J=1,N
2964           SUM0 = 0.0d0
2965           SUM2 = 0.0d0
2966          DO K=1,N
2967           SUM0 = SUM0 + S(J,K)*W(K,I)/EMU(K)
2968           SUM2 = SUM2 + S(J,K)*U1(K,I)/EMU(K)
2969          ENDDO
2970          A(J) = SUM0
2971          V1(J) = SUM2
2972         ENDDO
2973         DO J=1,N
2974          W(J,I) = A(J)
2975          U1(J,I) = V1(J)
2976         ENDDO
2977        ENDDO
2978        IF (ID.EQ.1) THEN
2979 !-------------upper boundary, 2nd-order, C-matrix is full (CC)
2980         DELTAU = ZTAU(2) - ZTAU(1)
2981         D2 = 0.25d0*DELTAU
2982         DO I=1,N
2983           D1 = EMU(I)/DELTAU
2984           DO J=1,N
2985            B(I,J) = B(I,J) + D2*W(I,J)
2986            CC(I,J) = D2*U1(I,J)
2987           ENDDO
2988           B(I,I) = B(I,I) + D1
2989           CC(I,I) = CC(I,I) - D1
2990 !         H(I) = H(I) + 2.0d0*D2*C1(I) + D1*SISOTP
2991           H(I) = H(I) + 2.0d0*D2*C1(I)
2992           A(I) = 0.0d0
2993         ENDDO
2994        ELSE
2995 !-------------lower boundary, 2nd-order, A-matrix is full (AA)
2996         DELTAU = ZTAU(ND) - ZTAU(ND-1)
2997         D2 = 0.25d0*DELTAU
2998         SURFAC = 4.0d0*ZREFL/(1.0d0 + ZREFL)
2999         DO I=1,N
3000           D1 = EMU(I)/DELTAU
3001           H(I) = H(I) - 2.0d0*D2*C1(I)
3002            SUM0 = 0.0d0
3003           DO J=1,N
3004            SUM0 = SUM0 + W(I,J)
3005           ENDDO
3006            SUM0 = D1 + D2*SUM0
3007            SUM1 = SURFAC*SUM0
3008           DO J=1,N
3009            B(I,J) = B(I,J) + D2*W(I,J) - SUM1*EMU(J)*WT(J)
3010           ENDDO
3011           B(I,I) = B(I,I) + D1
3012           H(I) = H(I) + SUM0*ZFLUX
3013           DO J=1,N
3014            AA(I,J) = - D2*U1(I,J)
3015           ENDDO
3016            AA(I,I) = AA(I,I) + D1
3017            C1(I) = 0.0d0
3018         ENDDO
3019        ENDIF
3020 !------------intermediate points:  can be even or odd, A & C diagonal
3021       ELSE
3022         DELTAU = ZTAU(ID+1) - ZTAU(ID-1)
3023         MSTART = M + MOD(ID+1,2)
3024         DO I=1,N
3025           A(I) = EMU(I)/DELTAU
3026           C1(I) = -A(I)
3027            SUM0 = 0.0d0
3028           DO IM=MSTART,MFIT,2
3029            SUM0 = SUM0 + POMEGA(IM,ID)*PM(I,IM)*PM0(IM)
3030           ENDDO
3031           H(I) = SUM0*FZ(ID)
3032           DO J=1,I
3033             SUM0 = 0.0d0
3034            DO IM=MSTART,MFIT,2
3035             SUM0 = SUM0 + POMEGA(IM,ID)*PM(I,IM)*PM(J,IM)
3036            ENDDO
3037             B(I,J) =  - SUM0*WT(J)
3038             B(J,I) =  - SUM0*WT(I)
3039           ENDDO
3040           B(I,I) = B(I,I) + 1.0d0
3041         ENDDO
3042       ENDIF
3043       RETURN
3044       END SUBROUTINE GEN
3045
3046       SUBROUTINE LEGND0 (X,PL,N)
3047 !---Calculates ORDINARY LEGENDRE fns of X (real)
3048 !---   from P[0] = PL(1) = 1,  P[1] = X, .... P[N-1] = PL(N)
3049       IMPLICIT NONE
3050       INTEGER N,I
3051       REAL*8 X,PL(N),DEN
3052 !---Always does PL(2) = P[1]
3053         PL(1) = 1.D0
3054         PL(2) = X
3055         DO I=3,N
3056          DEN = (I-1)
3057          PL(I) = PL(I-1)*X*(2.d0-1.D0/DEN) - PL(I-2)*(1.d0-1.D0/DEN)
3058         ENDDO
3059       RETURN
3060       END SUBROUTINE LEGND0
3061
3062       SUBROUTINE MATIN4 (A)
3063 !-----------------------------------------------------------------------
3064 !  invert 4x4 matrix A(4,4) in place with L-U decomposition (mjp, old...)
3065 !-----------------------------------------------------------------------
3066       IMPLICIT NONE
3067       REAL*8 A(4,4)
3068 !---SETUP L AND U
3069       A(2,1) = A(2,1)/A(1,1)
3070       A(2,2) = A(2,2)-A(2,1)*A(1,2)
3071       A(2,3) = A(2,3)-A(2,1)*A(1,3)
3072       A(2,4) = A(2,4)-A(2,1)*A(1,4)
3073       A(3,1) = A(3,1)/A(1,1)
3074       A(3,2) = (A(3,2)-A(3,1)*A(1,2))/A(2,2)
3075       A(3,3) = A(3,3)-A(3,1)*A(1,3)-A(3,2)*A(2,3)
3076       A(3,4) = A(3,4)-A(3,1)*A(1,4)-A(3,2)*A(2,4)
3077       A(4,1) = A(4,1)/A(1,1)
3078       A(4,2) = (A(4,2)-A(4,1)*A(1,2))/A(2,2)
3079       A(4,3) = (A(4,3)-A(4,1)*A(1,3)-A(4,2)*A(2,3))/A(3,3)
3080       A(4,4) = A(4,4)-A(4,1)*A(1,4)-A(4,2)*A(2,4)-A(4,3)*A(3,4)
3081 !---INVERT L
3082       A(4,3) = -A(4,3)
3083       A(4,2) = -A(4,2)-A(4,3)*A(3,2)
3084       A(4,1) = -A(4,1)-A(4,2)*A(2,1)-A(4,3)*A(3,1)
3085       A(3,2) = -A(3,2)
3086       A(3,1) = -A(3,1)-A(3,2)*A(2,1)
3087       A(2,1) = -A(2,1)
3088 !---INVERT U
3089       A(4,4) = 1.D0/A(4,4)
3090       A(3,4) = -A(3,4)*A(4,4)/A(3,3)
3091       A(3,3) = 1.D0/A(3,3)
3092       A(2,4) = -(A(2,3)*A(3,4)+A(2,4)*A(4,4))/A(2,2)
3093       A(2,3) = -A(2,3)*A(3,3)/A(2,2)
3094       A(2,2) = 1.D0/A(2,2)
3095       A(1,4) = -(A(1,2)*A(2,4)+A(1,3)*A(3,4)+A(1,4)*A(4,4))/A(1,1)
3096       A(1,3) = -(A(1,2)*A(2,3)+A(1,3)*A(3,3))/A(1,1)
3097       A(1,2) = -A(1,2)*A(2,2)/A(1,1)
3098       A(1,1) = 1.D0/A(1,1)
3099 !---MULTIPLY (U-INVERSE)*(L-INVERSE)
3100       A(1,1) = A(1,1)+A(1,2)*A(2,1)+A(1,3)*A(3,1)+A(1,4)*A(4,1)
3101       A(1,2) = A(1,2)+A(1,3)*A(3,2)+A(1,4)*A(4,2)
3102       A(1,3) = A(1,3)+A(1,4)*A(4,3)
3103       A(2,1) = A(2,2)*A(2,1)+A(2,3)*A(3,1)+A(2,4)*A(4,1)
3104       A(2,2) = A(2,2)+A(2,3)*A(3,2)+A(2,4)*A(4,2)
3105       A(2,3) = A(2,3)+A(2,4)*A(4,3)
3106       A(3,1) = A(3,3)*A(3,1)+A(3,4)*A(4,1)
3107       A(3,2) = A(3,3)*A(3,2)+A(3,4)*A(4,2)
3108       A(3,3) = A(3,3)+A(3,4)*A(4,3)
3109       A(4,1) = A(4,4)*A(4,1)
3110       A(4,2) = A(4,4)*A(4,2)
3111       A(4,3) = A(4,4)*A(4,3)
3112       RETURN
3113       END SUBROUTINE MATIN4
3114
3115       SUBROUTINE GAUSSP (N,XPT,XWT)
3116 !-----------------------------------------------------------------------
3117 !  Loads in pre-set Gauss points for 4 angles from 0 to +1 in cos(theta)=mu
3118 !-----------------------------------------------------------------------
3119       IMPLICIT NONE
3120       INTEGER N,I
3121       REAL*8  XPT(N),XWT(N)
3122       REAL*8 GPT4(4),GWT4(4)
3123       DATA GPT4/.06943184420297D0,.33000947820757D0,.66999052179243D0,   &
3124                 .93056815579703D0/
3125       DATA GWT4/.17392742256873D0,.32607257743127D0,.32607257743127D0,   &
3126                 .17392742256873D0/
3127       N = 4
3128       DO I=1,N
3129         XPT(I) = GPT4(I)
3130         XWT(I) = GWT4(I)
3131       ENDDO
3132       RETURN
3133       END SUBROUTINE GAUSSP
3134 !
3135       subroutine aeroden(zz,v,aerd)
3136
3137 ! purpose:   find number density of boundary layer aerosols, aerd,
3138 !            at a given altitude, zz, and for a specified visibility
3139 ! input:
3140 !   zz       altitude   (km)
3141 !   v        visibility for a horizontal surface path (km)
3142 ! output:
3143 !   aerd     aerosol density at altitude z
3144
3145 ! the vertical distribution of the boundary layer aerosol density is
3146 ! based on the 5s vertical profile models for 5 and 23 km visibility.
3147 ! above 5 km, the aden05 and aden23 models are the same
3148 ! below 5 km, the models differ as follows;
3149 ! aden05     0.99 km scale height (94% of extinction occurs below 5 km)
3150 ! aden23     1.45 km scale heigth (80% of extinction occurs below 5 km)
3151 !
3152
3153         implicit none
3154         integer mz, nz
3155         parameter (mz=33)
3156         real v,aerd
3157         real*8 zz ! compatability with fastj
3158         real alt, aden05, aden23, aer05,aer23
3159       dimension alt(mz),aden05(mz),aden23(mz)
3160 !jdf  dimension zbaer(*),dbaer(*)
3161
3162         real z, f, wth
3163         integer k, kp
3164       save alt,aden05,aden23,nz
3165
3166
3167       data nz/mz/
3168
3169       data alt/   &
3170               0.0,        1.0,        2.0,        3.0,        4.0,   &
3171               5.0,        6.0,        7.0,        8.0,        9.0,   &
3172              10.0,       11.0,       12.0,       13.0,       14.0,   &
3173              15.0,       16.0,       17.0,       18.0,       19.0,   &
3174              20.0,       21.0,       22.0,       23.0,       24.0,   &
3175              25.0,       30.0,       35.0,       40.0,       45.0,   &
3176              50.0,       70.0,      100.0/
3177       data aden05/   &
3178         1.378E+04,  5.030E+03,  1.844E+03,  6.731E+02,  2.453E+02,   &
3179         8.987E+01,  6.337E+01,  5.890E+01,  6.069E+01,  5.818E+01,   &
3180         5.675E+01,  5.317E+01,  5.585E+01,  5.156E+01,  5.048E+01,   &
3181         4.744E+01,  4.511E+01,  4.458E+01,  4.314E+01,  3.634E+01,   &
3182         2.667E+01,  1.933E+01,  1.455E+01,  1.113E+01,  8.826E+00,   &
3183         7.429E+00,  2.238E+00,  5.890E-01,  1.550E-01,  4.082E-02,   &
3184         1.078E-02,  5.550E-05,  1.969E-08/
3185       data aden23/   &
3186         2.828E+03,  1.244E+03,  5.371E+02,  2.256E+02,  1.192E+02,   &
3187         8.987E+01,  6.337E+01,  5.890E+01,  6.069E+01,  5.818E+01,   &
3188         5.675E+01,  5.317E+01,  5.585E+01,  5.156E+01,  5.048E+01,   &
3189         4.744E+01,  4.511E+01,  4.458E+01,  4.314E+01,  3.634E+01,   &
3190         2.667E+01,  1.933E+01,  1.455E+01,  1.113E+01,  8.826E+00,   &
3191         7.429E+00,  2.238E+00,  5.890E-01,  1.550E-01,  4.082E-02,   &
3192         1.078E-02,  5.550E-05,  1.969E-08/
3193 !
3194       z=max(0.,min(100.,real(zz)))
3195       aerd=0.
3196       if(z.gt.alt(nz)) return
3197
3198       call locate(alt,nz,z,k)
3199       kp=k+1
3200       f=(z-alt(k))/(alt(kp)-alt(k))
3201
3202       if(min(aden05(k),aden05(kp),aden23(k),aden23(kp)).le.0.) then
3203         aer05=aden05(k)*(1.-f)+aden05(kp)*f
3204         aer23=aden23(k)*(1.-f)+aden23(kp)*f
3205       else
3206         aer05=aden05(k)*(aden05(kp)/aden05(k))**f
3207         aer23=aden23(k)*(aden23(kp)/aden23(k))**f
3208       endif
3209
3210       wth=(1./v-1/5.)/(1./23.-1./5.)
3211       wth=max(0.,min(1.,wth))
3212
3213       aerd=(1.-wth)*aer05+wth*aer23
3214
3215 !      write(*,*) 'aeroden k,kp,z,aer05(k),aer05(kp),f,aerd'
3216 !      write(*,'(2i5,1p5e11.3)') k,kp,z,aden05(k),aden05(kp),f,aerd
3217
3218       return
3219         end subroutine aeroden
3220 !=======================================================================
3221       subroutine locate(xx,n,x,j)
3222 !
3223 ! purpose:  given an array xx of length n, and given a value X, returns
3224 !           a value J such that X is between xx(j) and xx(j+1). xx must
3225 !           be monotonic, either increasing of decreasing. this function
3226 !           returns j=1 or j=n-1 if x is out of range.
3227 !c
3228 ! input:
3229 !   xx      monitonic table
3230 !   n       size of xx
3231 !   x       single floating point value perhaps within the range of xx
3232 !
3233 ! output:
3234 !           function returns index value j, such that
3235 !
3236 !            for an increasing table
3237 !
3238 !                xx(j) .lt. x .le. xx(j+1),
3239 !                j=1    for x .lt. xx(1)
3240 !                j=n-1  for x .gt. xx(n)
3241 !
3242 !            for a decreasing table
3243 !                xx(j) .le. x .lt. xx(j+1)
3244 !                j=n-1  for x .lt. xx(n)
3245 !                j=1    for x .gt. xx(1)
3246 !
3247       implicit none
3248       integer j,n
3249       real x,xx(n)
3250       integer jl,jm,ju
3251
3252 ! . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3253
3254       if(x.eq.xx(1)) then
3255         j=1
3256         return
3257       endif
3258       if(x.eq.xx(n)) then
3259         j=n-1
3260         return
3261       endif
3262       jl=1
3263       ju=n
3264 10    if(ju-jl.gt.1) then
3265         jm=(ju+jl)/2
3266          if((xx(n).gt.xx(1)).eqv.(x.gt.xx(jm)))then
3267           jl=jm
3268         else
3269           ju=jm
3270         endif
3271       goto 10
3272       endif
3273       j=jl
3274       return
3275       end subroutine locate
3276 !************************************************************************
3277       subroutine rd_tjpl2
3278 !-----------------------------------------------------------------------
3279 !  set wavelength bins, solar fluxes, Rayleigh parameters, temperature-
3280 !  dependent cross sections and Rayleigh/aerosol scattering phase functions
3281 !  with temperature dependences. Current data originates from JPL 2000
3282 !-----------------------------------------------------------------------
3283 !
3284 !     NJVAL    Number of species to calculate J-values for
3285 !     NWWW     Number of wavelength bins, from NW1:NW2
3286 !     WBIN     Boundaries of wavelength bins
3287 !     WL       Centres of wavelength bins - 'effective wavelength'
3288 !     FL       Solar flux incident on top of atmosphere (cm-2.s-1)
3289 !     QRAYL    Rayleigh parameters (effective cross-section) (cm2)
3290 !     QBC      Black Carbon absorption extinct. (specific cross-sect.) (m2/g)
3291 !     QO2      O2 cross-sections
3292 !     QO3      O3 cross-sections
3293 !     Q1D      O3 => O(1D) quantum yield
3294 !     TQQ      Temperature for supplied cross sections
3295 !     QQQ      Supplied cross sections in each wavelength bin (cm2)
3296 !     NAA      Number of categories for scattering phase functions
3297 !     QAA      Aerosol scattering phase functions
3298 !     NK       Number of wavelengths at which functions supplied (set as 4)
3299 !     WAA      Wavelengths for the NK supplied phase functions
3300 !     PAA      Phase function: first 8 terms of expansion
3301 !     RAA      Effective radius associated with aerosol type
3302 !     SSA      Single scattering albedo
3303 !
3304 !     npdep    Number of pressure dependencies
3305 !     zpdep    Pressure dependencies by wavelength bin
3306 !     jpdep    Index of cross sections requiring pressure dependence
3307 !     lpdep    Label for pressure dependence
3308 !
3309 !-----------------------------------------------------------------------
3310
3311         USE module_data_mosaic_other, only : kmaxd
3312         USE module_fastj_data
3313         USE module_peg_util, only:  peg_message, peg_error_fatal
3314
3315         IMPLICIT NONE
3316 !jdf
3317 ! Print Fast-J diagnostics if iprint /= 0
3318         integer, parameter :: iprint = 0
3319         integer, parameter :: single = 4        !compiler dependent value real*4
3320 !       integer, parameter :: double = 8        !compiler dependent value real*8
3321         integer,parameter :: ipar_fastj=1,jpar=1
3322 !       integer,parameter :: jppj=14        !Number of photolytic reactions supplied
3323         logical,parameter :: ldeg45=.false. !Logical flag for degraded CTM resolution
3324         integer lpar           !Number of levels in CTM
3325         integer jpnl           !Number of levels requiring chemistry
3326         real(kind=double), dimension(ipar_fastj) :: xgrd !  Longitude (midpoint, radians)
3327         real(kind=double), dimension(jpar) :: ygrd           !  Latitude  (midpoint, radians)
3328         real(kind=double), dimension(jpar) :: ydgrd          !  Latitude  (midpoint, degrees)
3329         real(kind=double), dimension(kmaxd+1) :: etaa    !  Eta(a) value for level boundaries
3330         real(kind=double), dimension(kmaxd+1) :: etab    !  Eta(b) value for level boundaries
3331         real(kind=double) ::  tau_fastj          !  Time of Day (hours, GMT)
3332         integer month_fastj        !  Number of month (1-12)
3333         integer iday_fastj         !  Day of year
3334 !jdf
3335          integer i, j, k
3336          character*7  lpdep(3)
3337          character*80 msg
3338
3339          if(NJVAL.gt.NS) then
3340 ! fastj input files are not set up for current situation
3341             write ( msg, '(a, 2i6)' )  &
3342               'FASTJ  # xsect supplied > max allowed  '  //  &
3343               'NJVAL NS ', NJVAL, NS
3344             call peg_message( lunerr, msg )
3345             msg =               &
3346              'FASTJ Setup Error: # xsect supplied > max allowed. Increase NS'
3347             call peg_error_fatal( lunerr, msg )
3348 !           write(6,300) NJVAL,NS
3349 !           stop
3350          endif
3351
3352          if(NAA.gt.NP) then
3353             write ( msg, '(a, 2i6)' )  &
3354               'FASTJ  # aerosol/cloud types > NP  '  //  &
3355               'NAA NP ', NAA ,NP
3356             call peg_message( lunerr, msg )
3357             msg =               &
3358               'FASTJ Setup Error: Too many phase functions supplied. Increase NP'
3359             call peg_error_fatal( lunerr, msg )
3360 !            write(6,350) NAA
3361 !            stop
3362          endif
3363 !---Zero index arrays
3364       do j=1,jppj
3365         jind(j)=0
3366       enddo
3367       do j=1,NJVAL
3368         jpdep(j)=0
3369       enddo
3370       do j=1,nh
3371         hzind(j)=0
3372       enddo
3373 !
3374 !---Set mapping index
3375       do j=1,NJVAL
3376         do k=1,jppj
3377           if(jlabel(k).eq.titlej(1,j)) jind(k)=j
3378 !       write(6,*)k,jind(k)  ! jcb
3379 !       write(6,*)jlabel(k),titlej(1,j)  ! jcb
3380         enddo
3381         do k=1,npdep
3382           if(lpdep(k).eq.titlej(1,j)) jpdep(j)=k
3383         enddo
3384       enddo
3385       do k=1,jppj
3386         if(jfacta(k).eq.0.d0) then
3387 !            write(6,*) 'Not using photolysis reaction ',k
3388            write ( msg, '(a, i6)' )  &
3389              'FASTJ  Not using photolysis reaction ' , k
3390            call peg_message( lunerr, msg )
3391         end if
3392         if(jind(k).eq.0) then
3393           if(jfacta(k).eq.0.d0) then
3394             jind(k)=1
3395           else
3396               write ( msg, '(a, i6)' )  &
3397                'FASTJ  Which J-rate for photolysis reaction ' , k
3398               call peg_message( lunerr, msg )
3399 !              write(6,*) 'Which J-rate for photolysis reaction ',k,' ?'
3400 !              stop
3401               msg = 'FASTJ subr rd_tjpl2 Unknown Jrate. Incorrect FASTJ setup'
3402               call peg_error_fatal( lunerr, msg )
3403           endif
3404         endif
3405       enddo
3406 ! Herzberg index
3407       i=0
3408       do j=1,nhz
3409         do k=1,jppj
3410           if(jlabel(k).eq.hzlab(j)) then
3411             i=i+1
3412             hzind(i)=k
3413             hztoa(i)=hztmp(j)*jfacta(k)
3414           endif
3415         enddo
3416       enddo
3417       nhz=i
3418       if(nhz.eq.0) then
3419         if(iprint.ne.0) then
3420            write ( msg, '(a)' )  &
3421             'FASTJ  Not using Herzberg bin '
3422            call peg_message( lunerr, msg )
3423 !           write(6,400)
3424          end if
3425       else
3426         if(iprint.ne.0) then
3427            write ( msg, '(a)' )  &
3428               'FASTJ Using Herzberg bin for: '
3429            call peg_message( lunerr, msg )
3430            write( msg, '(a,10a7)' )     &
3431               'FASTJ ', (jlabel(hzind(i)),i=1,nhz)
3432 !          write(6,420) (jlabel(hzind(i)),i=1,nhz)
3433         end if
3434       endif
3435
3436 ! 300   format(' Number of x-sections supplied to Fast-J: ',i3,/,   &
3437 !             ' Maximum number allowed (NS) only set to: ',i3,   &
3438 !             ' - increase in jv_cmn.h',/,   &
3439 !             'RESULTS WILL BE IN ERROR' )
3440 ! 350   format(' Too many phase functions supplied; increase NP to ',i2,   &
3441 !             /,'RESULTS WILL BE IN ERROR'  )
3442 !  400 format(' Not using Herzberg bin')
3443 !  420 format(' Using Herzberg bin for: ',10a7)
3444
3445
3446          return
3447          end subroutine rd_tjpl2
3448 !********************************************************************
3449
3450
3451 end module module_phot_fastj