wrfv2_fire/chem/module_input_chem_data.F

   1 !dis
   2 !dis    Open Source License/Disclaimer, Forecast Systems Laboratory
   3 !dis    NOAA/OAR/FSL, 325 Broadway Boulder, CO 80305
   4 !dis
   5 !dis    This software is distributed under the Open Source Definition,
   6 !dis    which may be found at http://www.opensource.org/osd.html.
   7 !dis
   8 !dis    In particular, redistribution and use in source and binary forms,
   9 !dis    with or without modification, are permitted provided that the
  10 !dis    following conditions are met:
  11 !dis
  12 !dis    - Redistributions of source code must retain this notice, this
  13 !dis    list of conditions and the following disclaimer.
  14 !dis
  15 !dis    - Redistributions in binary form must provide access to this
  16 !dis    notice, this list of conditions and the following disclaimer, and
  17 !dis    the underlying source code.
  18 !dis
  19 !dis    - All modifications to this software must be clearly documented,
  20 !dis    and are solely the responsibility of the agent making the
  21 !dis    modifications.
  22 !dis
  23 !dis    - If significant modifications or enhancements are made to this
  24 !dis    software, the FSL Software Policy Manager
  25 !dis    (softwaremgr@fsl.noaa.gov) should be notified.
  26 !dis
  27 !dis    THIS SOFTWARE AND ITS DOCUMENTATION ARE IN THE PUBLIC DOMAIN
  28 !dis    AND ARE FURNISHED "AS IS."  THE AUTHORS, THE UNITED STATES
  29 !dis    GOVERNMENT, ITS INSTRUMENTALITIES, OFFICERS, EMPLOYEES, AND
  30 !dis    AGENTS MAKE NO WARRANTY, EXPRESS OR IMPLIED, AS TO THE USEFULNESS
  31 !dis    OF THE SOFTWARE AND DOCUMENTATION FOR ANY PURPOSE.  THEY ASSUME
  32 !dis    NO RESPONSIBILITY (1) FOR THE USE OF THE SOFTWARE AND
  33 !dis    DOCUMENTATION; OR (2) TO PROVIDE TECHNICAL SUPPORT TO USERS.
  34 !dis
  35 !dis
  36
  37 !WRF:PACKAGE:IO
  38
  39
  40 !CPP directives to control ic/bc conditions...
  41 !(The directive in module_mosaic_initmixrats also needs to be set.)
  42 !  CASENAME = 0   Uses Texas August 2004 case values and profiles
  43 !             1   Uses same concentrations as TX, but uses different
  44 !                 profiles depending on the species. (NEAQS2004 case)
  45 #define CASENAME 0
  46
  47
  48 MODULE module_input_chem_data
  49
  50    USE module_io_domain
  51    USE module_domain
  52    USE module_data_sorgam, ONLY : conmin, rgasuniv, epsilc, grav
  53    USE module_get_file_names, ONLY : eligible_file_name, number_of_eligible_files, unix_ls
  54
  55    IMPLICIT NONE
  56
  57 !  REAL :: grav = 9.8104
  58    REAL, PARAMETER :: mwso4 = 96.0576
  59
  60 ! Variables for adaptive time steps...
  61 #if ( EM_CORE == 1 )
  62       TYPE(WRFU_Time), DIMENSION(max_domains) :: last_chem_time
  63 #endif
  64 #if ( NMM_CORE == 1)
  65       TYPE(WRFU_Time), DIMENSION(1) :: last_chem_time
  66 #endif
  67
  68 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  69 ! Initial atmospheric chemistry profile data
  70     INTEGER :: k_loop          ! Used for loop index
  71     INTEGER :: lo              ! number of chemicals in inital profile
  72     INTEGER :: logg            ! number of final chemical species (nch-1)
  73     INTEGER :: kx              ! number of vertical levels in temp profile
  74     INTEGER :: kxm1
  75
  76     PARAMETER( kx=16, kxm1=kx-1, logg=100, lo=34)
  77
  78     INTEGER, DIMENSION(logg)                     :: iref
  79
  80     real, allocatable :: ch4_lbc(:,:)
  81     real, allocatable :: n2o_lbc(:,:)
  82     real, allocatable :: h2_lbc(:,:)
  83     REAL, DIMENSION(logg)                        :: fracref
  84     REAL, DIMENSION(kx)                          :: dens
  85     REAL, DIMENSION(kx+1)                        :: zfa
  86     REAL, DIMENSION(kx+1)                        :: zfa_bdy
  87     REAL, DIMENSION(lo  ,kx)                     :: xl
  88     REAL                                         :: so4vaptoaer
  89     DATA so4vaptoaer/.999/
  90
  91     CHARACTER (LEN=4), DIMENSION(logg) :: ggnam
  92
  93 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  94 !
  95 !  The idealized profile is converted from the NALROM chemistry model.
  96 !
  97 ! The variable "iref" is the reference index
  98 !    "fracref" is the reference fraction correpsonding to iref
  99 !
 100 !  For example: WRF-Chem species number 1 for RADM is SO2.  iref(1)
 101 !  is 12, and XL(12,K) is the profile for SO2.
 102 !
 103 !  Note: NALROM calculates lumped OX (XL(1) =O3+NO2+HNO3+...) and a
 104 !  lumped NOX (XL(2)=NO+NO2+NO3+2N2O5+HO2NO2+HONO).  But XL(33) is
 105 !  strictly O3, and XL(34)=NOx=(NO+NO2 only).
 106 !
 107 !  Short-lived species are initialized to steady-state equilibrium -
 108 !  since they are short-lived.  The short-lived species within a lumped category
 109 !  (Ox , NOx, or NO3+N2O5 in our case) would be renormalized to the lumped class
 110 !  after the steady-state equilibrium concentrations are determined.
 111 !
 112 !  The following is a list of long-lived species
 113 !
 114 !      NAMEL( 1)='OX        '
 115 !      NAMEL( 2)='NOX       '
 116 !      NAMEL( 3)='HNO3      '
 117 !      NAMEL( 4)='H2O2      '
 118 !      NAMEL( 5)='CH3OOH    '
 119 !      NAMEL( 6)='CO        '
 120 !      NAMEL( 7)='ISOPRENE  '
 121 !      NAMEL( 8)='CH2O      '
 122 !      NAMEL( 9)='CH3CHO    '
 123 !      NAMEL(10)='PAN       '
 124 !      NAMEL(11)='OTHER ALKA'
 125 !      NAMEL(12)='SO2       '
 126 !      NAMEL(13)='BUTANE    '
 127 !      NAMEL(14)='ETHENE    '
 128 !      NAMEL(15)='PROPENE+   '
 129 !      NAMEL(16)='PPN       '
 130 !      NAMEL(17)='MEK       '
 131 !      NAMEL(18)='RCHO      '
 132 !      NAMEL(19)='SO4=      '
 133 !      NAMEL(20)='MVK       '
 134 !      NAMEL(21)='MACR      '
 135 !      NAMEL(22)='HAC       '
 136 !      NAMEL(23)='MGLY      '
 137 !      NAMEL(24)='HPAN      '
 138 !      NAMEL(25)='MPAN      '
 139 !      NAMEL(26)='PROPANE   '
 140 !      NAMEL(27)='ACETYLENE '
 141 !      NAMEL(28)='OH        '
 142 !      NAMEL(29)='HO2       '
 143 !      NAMEL(30)='NO3+N2O5  '
 144 !      NAMEL(31)='HO2NO2    '
 145 !      NAMEL(32)='SUM RO2   '
 146 !      NAMEL(33)='OZONE     '
 147 !      NAMEL(34)='NOX       '
 148 !
 149 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1
 150
 151      DATA dens/ 2.738E+18, 5.220E+18, 7.427E+18, 9.202E+18, &
 152                 1.109E+19, 1.313E+19, 1.525E+19, 1.736E+19, &
 153                 1.926E+19, 2.074E+19, 2.188E+19, 2.279E+19, &
 154                 2.342E+19, 2.384E+19, 2.414E+19, 2.434E+19  /
 155
 156 !     Profile heights in meters
 157 !     DATA ZFA/    0.,   85.,  212.,  385.,  603.,  960., 1430., 2010., &
 158 !               2850., 4010., 5340., 6900., 8510.,10200.,12100.,16000., &
 159 !              21000./
 160 #if ( ! EM_CORE == 0 )
 161       DATA ZFA_BDY/    0.,   85.,  212.,  385.,  603.,  960., 1430., 2010., &
 162                 2850., 4010., 5340., 6900., 8510.,10200.,12100.,16000., &
 163                21000./
 164
 165 !     Profile heights in meters
 166       DATA ZFA/    0.,   85.,  212.,  385.,  603.,  960., 1430., 2010., &
 167                 2850., 4010., 5340., 6900., 8510.,10200.,12100.,16000., &
 168                21000./
 169 #endif
 170 #if ( ! NMM_CORE == 0 )
 171
 172       DATA ZFA_BDY/    0.,   85.,  212.,  385.,  603.,  960., 1430., 2010., &
 173                 2850., 4010., 5340., 6900., 8510.,10200.,12100.,16000., &
 174                21000./
 175
 176 !     Profile pressure in hpa
 177       DATA ZFA/ 100000.,  98500.,  98000.,  96000.,  94000.,  90000., 85000., 75000., &
 178                  71000.,  65000.,  52000.,  48000.,  45000.,  30000., 25000., 20000., &
 179                   5000./
 180 #endif
 181
 182 !wig: To match the xl profile to the correct species, match WRF's p_<species>
 183 !     flag with iref(p_<species>-1) to get the value of the first index in xl,
 184 !     e.g. p_o3=6, iref(6-1)=1, so xl(1,:) is the ozone profile.
 185 !     See gasprofile_init_pnnl for an explination of what height
 186 !     each index represents.
 187       DATA (xl(1,k_loop),k_loop=1,kx) &
 188       / 1.68E-07, 1.68E-07, 5.79E-08, 5.24E-08, 5.26E-08, &
 189        5.16E-08, 4.83E-08, 4.50E-08, 4.16E-08, 3.80E-08, 3.56E-08, &
 190        3.35E-08, 3.15E-08, 3.08E-08, 3.06E-08, 3.00E-08/
 191
 192       DATA (xl(2,k_loop),k_loop=1,kx) &
 193       / 4.06E-10, 4.06E-10, 2.16E-10, 1.37E-10, 9.47E-11, &
 194        6.95E-11, 5.31E-11, 4.19E-11, 3.46E-11, 3.01E-11, 2.71E-11, &
 195        2.50E-11, 2.35E-11, 2.26E-11, 2.20E-11, 2.16E-11/
 196
 197       DATA (xl(3,k_loop),k_loop=1,kx) &
 198       / 9.84E-10, 9.84E-10, 5.66E-10, 4.24E-10, 3.26E-10, &
 199        2.06E-10, 1.12E-10, 7.33E-11, 7.03E-11, 7.52E-11, 7.96E-11, &
 200        7.56E-11, 7.27E-11, 7.07E-11, 7.00E-11, 7.00E-11/
 201
 202       DATA (xl(4,k_loop),k_loop=1,kx) &
 203       / 8.15E-10, 8.15E-10, 8.15E-10, 8.15E-10, 8.15E-10, &
 204        8.65E-10, 1.07E-09, 1.35E-09, 1.47E-09, 1.47E-09, 1.47E-09, &
 205        1.47E-09, 1.45E-09, 1.43E-09, 1.40E-09, 1.38E-09/
 206
 207       DATA (xl(5,k_loop),k_loop=1,kx) &
 208       / 4.16E-10, 4.16E-10, 4.16E-10, 4.16E-10, 4.16E-10, &
 209        4.46E-10, 5.57E-10, 1.11E-09, 1.63E-09, 1.63E-09, 1.63E-09, &
 210        1.63E-09, 1.61E-09, 1.59E-09, 1.57E-09, 1.54E-09/
 211
 212 !  CO is 70 ppbv at top, 80 throughout troposphere
 213       DATA (xl(6,k_loop),k_loop=1,kx)  / 7.00E-08, kxm1*8.00E-08/
 214
 215       DATA (xl(7,k_loop),k_loop=1,kx) &
 216       / 8.33E-29, 8.33E-29, 8.33E-29, 8.33E-29, 8.33E-29, &
 217        1.33E-28, 3.54E-28, 1.85E-28, 1.29E-29, 1.03E-30, 1.72E-31, &
 218        7.56E-32, 1.22E-31, 2.14E-31, 2.76E-31, 2.88E-31/
 219
 220       DATA (xl(8,k_loop),k_loop=1,kx) &
 221       / 9.17E-11, 9.17E-11, 9.17E-11, 9.17E-11, 9.17E-11, &
 222        1.03E-10, 1.55E-10, 2.68E-10, 4.47E-10, 4.59E-10, 4.72E-10, &
 223        4.91E-10, 5.05E-10, 5.13E-10, 5.14E-10, 5.11E-10/
 224       DATA (xl(9,k_loop),k_loop=1,kx) &
 225       / 7.10E-12, 7.10E-12, 7.10E-12, 7.10E-12, 7.10E-12, &
 226        7.36E-12, 1.02E-11, 2.03E-11, 2.98E-11, 3.01E-11, 3.05E-11, &
 227        3.08E-11, 3.08E-11, 3.06E-11, 3.03E-11, 2.99E-11/
 228       DATA (xl(10,k_loop),k_loop=1,kx) &
 229       / 4.00E-11, 4.00E-11, 4.00E-11, 3.27E-11, 2.51E-11, &
 230        2.61E-11, 2.20E-11, 1.69E-11, 1.60E-11, 1.47E-11, 1.37E-11, &
 231        1.30E-11, 1.24E-11, 1.20E-11, 1.18E-11, 1.17E-11/
 232       DATA (xl(11,k_loop),k_loop=1,kx) &
 233       / 1.15E-16, 1.15E-16, 2.46E-15, 2.30E-14, 1.38E-13, &
 234        6.25E-13, 2.31E-12, 7.32E-12, 1.87E-11, 3.68E-11, 6.10E-11, &
 235        9.05E-11, 1.22E-10, 1.50E-10, 1.70E-10, 1.85E-10/
 236       DATA (xl(12,k_loop),k_loop=1,kx) &
 237       / 1.00E-10, 1.00E-10, 1.00E-10, 1.00E-10, 1.00E-10, &
 238        1.00E-10, 1.00E-10, 1.00E-10, 1.00E-10, 1.00E-10, 1.00E-10, &
 239        1.00E-10, 1.00E-10, 1.00E-10, 1.00E-10, 1.00E-10/
 240       DATA (xl(13,k_loop),k_loop=1,kx) &
 241       / 1.26E-11, 1.26E-11, 2.02E-11, 2.50E-11, 3.02E-11, &
 242        4.28E-11, 6.62E-11, 1.08E-10, 1.54E-10, 2.15E-10, 2.67E-10, &
 243        3.24E-10, 3.67E-10, 3.97E-10, 4.16E-10, 4.31E-10/
 244       DATA (xl(14,k_loop),k_loop=1,kx) &
 245       / 1.15E-16, 1.15E-16, 2.46E-15, 2.30E-14, 1.38E-13, &
 246        6.25E-13, 2.31E-12, 7.32E-12, 1.87E-11, 3.68E-11, 6.10E-11, &
 247        9.05E-11, 1.22E-10, 1.50E-10, 1.70E-10, 1.85E-10/
 248       DATA (xl(15,k_loop),k_loop=1,kx) &
 249       / 1.00E-20, 1.00E-20, 6.18E-20, 4.18E-18, 1.23E-16, &
 250        2.13E-15, 2.50E-14, 2.21E-13, 1.30E-12, 4.66E-12, 1.21E-11, &
 251        2.54E-11, 4.47E-11, 6.63E-11, 8.37E-11, 9.76E-11/
 252       DATA (xl(16,k_loop),k_loop=1,kx) &
 253       / 1.23E-11, 1.23E-11, 1.23E-11, 1.23E-11, 1.23E-11, &
 254        1.20E-11, 9.43E-12, 3.97E-12, 1.19E-12, 1.11E-12, 9.93E-13, &
 255        8.66E-13, 7.78E-13, 7.26E-13, 7.04E-13, 6.88E-13/
 256       DATA (xl(17,k_loop),k_loop=1,kx) &
 257       / 1.43E-12, 1.43E-12, 1.43E-12, 1.43E-12, 1.43E-12, &
 258        1.50E-12, 2.64E-12, 8.90E-12, 1.29E-11, 1.30E-11, 1.32E-11, &
 259        1.32E-11, 1.31E-11, 1.30E-11, 1.29E-11, 1.27E-11/
 260       DATA (xl(18,k_loop),k_loop=1,kx) &
 261        / 3.61E-13, 3.61E-13, 3.61E-13, 3.61E-13, 3.61E-13, &
 262        3.58E-13, 5.22E-13, 1.75E-12, 2.59E-12, 2.62E-12, 2.64E-12, &
 263        2.66E-12, 2.65E-12, 2.62E-12, 2.60E-12, 2.57E-12/
 264       DATA (xl(19,k_loop),k_loop=1,kx) &
 265        / 5.00E-11, 5.00E-11, 5.00E-11, 5.00E-11, 5.00E-11, &
 266        5.00E-11, 5.00E-11, 5.00E-11, 5.00E-11, 5.00E-11, 5.00E-11, &
 267        5.00E-11, 5.00E-11, 5.00E-11, 5.00E-11, 5.00E-11/
 268
 269       DATA (xl(20,k_loop),k_loop=1,kx)/kx*1.E-20/
 270       DATA (xl(21,k_loop),k_loop=1,kx)/kx*1.E-20/
 271       DATA (xl(22,k_loop),k_loop=1,kx)/kx*1.E-20/
 272       DATA (xl(23,k_loop),k_loop=1,kx)/kx*1.E-20/
 273       DATA (xl(24,k_loop),k_loop=1,kx)/kx*1.E-20/
 274       DATA (xl(25,k_loop),k_loop=1,kx)/kx*1.E-20/
 275
 276 ! Propane - Gregory PEM-West A 25 ppt median marine boundary layer
 277       DATA (xl(26,k_loop),k_loop=1,kx) &
 278       /5.00E-13, 1.24E-12, 2.21E-12, 3.27E-12, 4.71E-12, &
 279        6.64E-12, 9.06E-12, 1.19E-11, 1.47E-11, 1.72E-11, &
 280        1.93E-11, 2.11E-11, 2.24E-11, 2.34E-11, 2.42E-11, 2.48E-11/
 281 ! Acetylene - Gregory PEM-West A 53 ppt median marine boundary layer
 282       DATA (xl(27,k_loop),k_loop=1,kx) &
 283       /1.00E-12, 2.48E-12, 4.42E-12, 6.53E-12, 9.42E-12, &
 284        1.33E-11, 1.81E-11, 2.37E-11, 2.95E-11, 3.44E-11, &
 285        3.85E-11, 4.22E-11, 4.49E-11, 4.69E-11, 4.84E-11, 4.95E-11/
 286 ! OH
 287       DATA (xl(28,k_loop),k_loop=1,kx) &
 288        / 9.80E+06, 9.80E+06, 4.89E+06, 2.42E+06, 1.37E+06, &
 289        9.18E+05, 7.29E+05, 6.26E+05, 5.01E+05, 4.33E+05, 4.05E+05, &
 290        3.27E+05, 2.54E+05, 2.03E+05, 1.74E+05, 1.52E+05/
 291 ! HO2
 292       DATA (xl(29,k_loop),k_loop=1,kx) &
 293        / 5.74E+07, 5.74E+07, 7.42E+07, 8.38E+07, 8.87E+07, &
 294        9.76E+07, 1.15E+08, 1.34E+08, 1.46E+08, 1.44E+08, 1.40E+08, &
 295        1.36E+08, 1.31E+08, 1.28E+08, 1.26E+08, 1.26E+08/
 296 ! NO3+N2O5
 297       DATA (xl(30,k_loop),k_loop=1,kx) &
 298        / 5.52E+05, 5.52E+05, 3.04E+05, 2.68E+05, 2.32E+05, &
 299        1.66E+05, 1.57E+05, 1.72E+05, 1.98E+05, 2.22E+05, 2.43E+05, &
 300        2.75E+05, 3.00E+05, 3.18E+05, 3.32E+05, 3.39E+05/
 301 ! HO2NO2
 302       DATA (xl(31,k_loop),k_loop=1,kx) &
 303        / 7.25E+07, 7.25E+07, 6.36E+07, 5.55E+07, 4.94E+07, &
 304        3.66E+07, 2.01E+07, 9.57E+06, 4.75E+06, 2.37E+06, 1.62E+06, &
 305        9.86E+05, 7.05E+05, 5.63E+05, 4.86E+05, 4.41E+05/
 306 ! Sum of RO2 &
 307       DATA (xl(32,k_loop),k_loop=1,kx) &
 308        / 9.14E+06, 9.14E+06, 1.46E+07, 2.14E+07, 2.76E+07, &
 309        3.62E+07, 5.47E+07, 1.19E+08, 2.05E+08, 2.25E+08, 2.39E+08, &
 310        2.58E+08, 2.82E+08, 2.99E+08, 3.08E+08, 3.15E+08/
 311 ! O3 <--This is not the O3 used for RADM2 or CBMZ (wig)
 312       DATA (xl(33,k_loop),k_loop=1,kx) &
 313        / 8.36E+11, 8.36E+11, 4.26E+11, 4.96E+11, 6.05E+11, &
 314        6.93E+11, 7.40E+11, 7.74E+11, 7.82E+11, 7.75E+11, 7.69E+11, &
 315        7.59E+11, 7.54E+11, 7.50E+11, 7.47E+11, 7.45E+11/
 316 ! NOx (NO+NO2)
 317       DATA (xl(34,k_loop),k_loop=1,kx) &
 318        / 1.94E+09, 1.94E+09, 1.53E+09, 1.24E+09, 1.04E+09, &
 319        8.96E+08, 7.94E+08, 7.11E+08, 6.44E+08, 6.00E+08, 5.70E+08, &
 320        5.49E+08, 5.35E+08, 5.28E+08, 5.24E+08, 5.23E+08/
 321
 322 CONTAINS
 323
 324 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 325 ! Sets up the cross reference mapping indices and fractional
 326 ! apportionment of the default species profiles for use with
 327 ! ICs and BCs.
 328 ! William.Gustafson@pnl.gov; 2-May-2007
 329 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 330 SUBROUTINE setup_gasprofile_maps(chem_opt, numgas)
 331   integer, intent(in) :: chem_opt, numgas
 332
 333   select case(chem_opt)
 334   case (RADM2, RADM2_KPP, RADM2SORG, RADM2SORG_AQ, RADM2SORG_KPP, &
 335         RACM_KPP, RACMPM_KPP, RACM_MIM_KPP, RACMSORG_AQ, RACMSORG_KPP, &
 336         GOCARTRACM_KPP, GOCARTRADM2,GOCARTRADM2_KPP,CHEM_TRACER, CHEM_TRACE2)
 337      call setup_gasprofile_map_radm_racm
 338
 339   case (CBMZ, CBMZ_BB, CBMZ_BB_KPP, CBMZ_MOSAIC_4BIN, CBMZ_MOSAIC_8BIN, &
 340         CBMZ_MOSAIC_4BIN_AQ, CBMZ_MOSAIC_8BIN_AQ)
 341      call setup_gasprofile_map_cbmz(numgas)
 342
 343   case (CBM4_KPP)
 344      call setup_gasprofile_map_cbm4(numgas)
 345
 346   case (GOCART_SIMPLE)
 347      call wrf_debug("setup_profile_maps: nothing done for gocart simple")
 348
 349   case (MOZART_KPP)
 350      call wrf_debug("setup_profile_maps: nothing done for mozart_kpp")
 351
 352   case (MOZCART_KPP)
 353      call wrf_debug("setup_profile_maps: nothing done for mozcart_kpp")
 354
 355   case default
 356      call wrf_error_fatal("setup_profile_maps: could not decipher chem_opt value")
 357
 358   end select
 359
 360 END SUBROUTINE setup_gasprofile_maps
 361
 362 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 363 ! Maps for indices and fractional apportionment of the default
 364 ! species profiles to the RADM2 and RACM species for ICs and BCs.
 365 ! William.Gustafson@pnl.gov; 2-May-2007
 366 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 367 SUBROUTINE setup_gasprofile_map_radm_racm
 368
 369   iref(:)    = 7 !default value
 370   iref(1:41) = (/12,19,2,2,1,3,4,9,8,5,5,32,6,6,6,30,30,10,26,13,11,6,6, &
 371                  14,15,15,23,23,32,16,23,31,17,23,23,23,23,23,7,28,29/)
 372
 373   fracref(:)    = 1. !default value
 374   fracref(1:41) = (/1.,1.,.75,.25,1.,1.,1.,1.,1.,1., &
 375                     .5,.5,6.25E-4,7.5E-4,6.25E-5,.1, &
 376                     .9,1.,1.,1.,1.,8.E-3,1.,1.,1.,.5,&
 377                     1.,1.,.5,1.,1.,1.,1.,1.,1.,1.,1.,&
 378                     1.,1.,1.,1./)
 379
 380   ggnam(:) = 'JUNK' !default value
 381   ggnam(1:41) = (/ 'SO2 ','SULF','NO2 ','NO  ','O3  ','HNO3',    &
 382                    'H2O2','ALD ','HCHO','OP1 ','OP2 ','PAA ',    &
 383                    'ORA1','ORA2','NH3 ','N2O5','NO3 ','PAN ',    &
 384                    'HC3 ','HC5 ','HC8 ','ETH ','CO  ','OL2 ',    &
 385                    'OLT ','OLI ','TOL ','XYL ','ACO3','TPAN',    &
 386                    'HONO','HNO4','KET ','GLY ','MGLY','DCB ',    &
 387                    'ONIT','CSL ','ISO ','HO  ','HO2 '           /)
 388
 389 END SUBROUTINE setup_gasprofile_map_radm_racm
 390
 391 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 392 ! Maps for indices and fractional apportionment of the default
 393 ! species profiles to the CBMZ species for ICs and BCs. The profiles
 394 ! for HC3, HC5, and HC8 used for RADM are tacked onto the end of the
 395 ! CBMZ list in order to construct the PAR species.
 396 ! William.Gustafson@pnl.gov; 2-May-2007
 397 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 398 SUBROUTINE setup_gasprofile_map_cbmz(numgas)
 399   integer, intent(in) :: numgas
 400   integer, parameter :: listlast = 33
 401   iref(:)    = 7 !default value
 402   iref(1:listlast) = (/12,19, 2, 2, 1, 3, &
 403                         4, 9, 8, 5, 5, 6, &
 404                         6, 6,30,30,10, 6, &
 405                         6,14,15,15,23,23, &
 406                        23,31,17,23,23,23, &
 407                         7,28,29          /)
 408
 409   fracref(:)    = 1. !default value
 410   fracref(1:listlast) = (/1.,1.,.75,.25,1.,1., &
 411                           1.,1.,1.,1.,.5,6.25E-4, &
 412                           7.5E-4,6.25E-5,.1,.9,1.,8.E-3, &
 413                           1.,1.,1.,.5,1.,1.,   &
 414                           1.,1.,1.,1.,1.,1.,   &
 415                           1.,1.,1.            /)
 416
 417   ggnam(:) = 'JUNK' !default value - not really all junk, the
 418                     !last species all just goto the default
 419   ggnam(1:listlast) = (/ 'SO2 ','SULF','NO2 ','NO  ','O3  ','HNO3', &
 420                          'H2O2','ALD ','HCHO','OP1 ','OP2 ','ORA1', &
 421                          'ORA2','NH3 ','N2O5','NO3 ','PAN ','ETH ', &
 422                          'CO  ','OL2 ','OLT ','OLI ','TOL ','XYL ', &
 423                          'HONO','HNO4','KET ','MGLY','ONIT','CSL ', &
 424                          'ISO ','HO  ','HO2 '                      /)
 425 !
 426 ! After the CBMZ species, add the RADM HC3, HC5, and HC8 to be used
 427 ! for constructing PAR..
 428 !
 429   if( numgas < listlast ) &
 430        call wrf_error_fatal("numgas < listlast in setup_gasprofile_map_cbmz")
 431   iref(numgas+1:numgas+3)    = (/   26,     13,    11/)
 432   fracref(numgas+1:numgas+3) = (/   1.,     1.,    1./)
 433   ggnam(numgas+1:numgas+3)   = (/'HC3 ','HC5 ','HC8 '/)
 434
 435 END SUBROUTINE setup_gasprofile_map_cbmz
 436
 437 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 438
 439 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 440 ! Maps for indices and fractional apportionment of the default
 441 ! species profiles to the CBM4 species for ICs and BCs. The profiles
 442 ! for HC3, HC5, and HC8 used for RADM are tacked onto the end of the
 443 ! CBM4 list in order to construct the PAR species.
 444 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 445 SUBROUTINE setup_gasprofile_map_cbm4(numgas)
 446   integer, intent(in) :: numgas
 447   integer, parameter :: listlast = 24
 448   iref(:)    = 7 !default value
 449   iref(1:listlast) = (/12,19, 2, 2, 1, 3, &
 450                         4, 9, 8,  &
 451                         6,30,30,10, 6, &
 452                         6,15,23,23, &
 453                         23,23,23, &
 454                         7,28,29          /)
 455
 456   fracref(:)    = 1. !default value
 457   fracref(1:listlast) = (/1.,1.,.75,.25,1.,1., &
 458                           1.,1.,1., &
 459                           6.25E-5,.1,.9,1.,1., &
 460                           1.,7.,1.,1.,   &
 461                           1.,1.,1.,   &
 462                           1.,1.,1.            /)
 463
 464   ggnam(:) = 'JUNK' !default value - not really all junk, the
 465                     !last species all just goto the default
 466   ggnam(1:listlast) = (/ 'SO2 ','SULF','NO2 ','NO  ','O3  ','HNO3', &
 467                          'H2O2','ALD2','HCHO', &
 468                          'NH3 ','N2O5','NO3 ','PAN ','ETH ', &
 469                          'CO  ','OLE ','TOL ','XYL ', &
 470                          'HONO','ONIT','CRES', &
 471                          'ISO ','HO  ','HO2 '                      /)
 472 !
 473 ! After the CBM4 species, add the RADM HC3, HC5, and HC8 to be used
 474 ! for constructing PAR..
 475 !
 476   if( numgas < listlast ) &
 477        call wrf_error_fatal("numgas < listlast in setup_gasprofile_map_cbm4")
 478   iref(numgas+1:numgas+7)    = (/   26,     13,    11,  17,   15,     15,   6/)
 479   fracref(numgas+1:numgas+7) = (/   1.,     1.,    1.,   1.,   .5,   1.,   7.5E-4/)
 480   ggnam(numgas+1:numgas+7)   = (/'HC3 ','HC5 ','HC8 ','KET ','OLI ','OLT ','ORA2'/)
 481
 482 END SUBROUTINE setup_gasprofile_map_cbm4
 483
 484 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 485
 486 SUBROUTINE input_ext_chem_file (si_grid)
 487 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 488
 489    IMPLICIT NONE
 490
 491    TYPE(domain)           ::  si_grid
 492
 493    INTEGER :: i , j , k, l, &
 494               ids, ide, jds, jde, kds, kde,    &
 495               ims, ime, jms, jme, kms, kme,    &
 496               ips, ipe, jps, jpe, kps, kpe
 497    INTEGER :: si_jday
 498    INTEGER :: dat_jday,dat_year,dat_month,dat_day,dat_hour,dat_min,dat_sec
 499    INTEGER :: time_loop_max , time_loop
 500    INTEGER, DIMENSION(2) :: num_steps
 501    INTEGER :: fid, ierr, rc
 502    INTEGER :: status_next_var
 503    INTEGER :: debug_level
 504    INTEGER :: si_year,si_month,si_day,si_hour,si_min,si_sec
 505    INTEGER :: total_time_sec , file_counter
 506 #if ( ! NMM_CORE == 0 )
 507    REAL, ALLOCATABLE, DIMENSION(:,:,:) :: pint
 508    REAL, ALLOCATABLE, DIMENSION(:,:) :: pdsl
 509 #endif
 510
 511    LOGICAL :: input_from_file , need_new_file
 512
 513    REAL, ALLOCATABLE, DIMENSION(:,:,:) :: si_zsigf, si_zsig
 514    REAL, ALLOCATABLE, DIMENSION(:,:,:) :: ch_zsigf, ch_zsig, ozone
 515    REAL    :: ext_time, dat_time
 516    REAL    :: wgt0,wgt_time1,wgt_time2
 517
 518    CHARACTER (LEN=80) :: inpname, message
 519    CHARACTER (LEN=19) :: date_string
 520    CHARACTER (LEN=19) :: extract_date, use_date, current_date_char
 521    CHARACTER*80                           :: timestr
 522
 523    TYPE (grid_config_rec_type)              :: config_flags
 524    TYPE(domain) , POINTER ::  null_domain, chem_grid, chgrid
 525    TYPE(domain) , POINTER ::  chem_grid2, chgrid2
 526 !   TYPE(ESMF_Time)                        :: CurrTime
 527
 528    !  Interface block for routine that passes pointers and needs to know that they
 529    !  are receiving pointers.
 530
 531
 532      CALL model_to_grid_config_rec ( si_grid%id , model_config_rec , config_flags )
 533       !  After current_date has been set, fill in the julgmt stuff.
 534
 535       CALL geth_julgmt ( config_flags%julyr , config_flags%julday , config_flags%gmt )
 536
 537    WRITE ( extract_date , FMT = '(I4.4,"-",I2.2,"-",I2.2,"_",I2.2,":",I2.2,":",I2.2)' ) &
 538            model_config_rec%start_year  (si_grid%id) , &
 539            model_config_rec%start_month (si_grid%id) , &
 540            model_config_rec%start_day   (si_grid%id) , &
 541            model_config_rec%start_hour  (si_grid%id) , &
 542            model_config_rec%start_minute(si_grid%id) , &
 543            model_config_rec%start_second(si_grid%id)
 544
 545    write(message,'(A,A)') 'Subroutine input_chem: finding data at date/time: ',extract_date
 546    CALL  wrf_message ( TRIM(message) )
 547
 548
 549    !  And here is an instance of using the information in the NAMELIST.
 550
 551    CALL nl_get_debug_level ( 1,debug_level )
 552    CALL set_wrf_debug_level ( debug_level )
 553
 554
 555    !  Allocated and configure the mother domain.  Since we are in the nesting down
 556    !  mode, we know a) we got a nest, and b) we only got 1 nest.
 557
 558    NULLIFY( null_domain )
 559
 560    CALL       wrf_debug ( 100 , 'wrfchem:input_chem: calling alloc_and_configure_domain 1' )
 561    ! Note that this is *not* the intended use of alloc_and_configure_domain()
 562    ! It does not seem to hurt anything, yet...
 563
 564 !   if( si_grid%id .EQ. 1) then
 565    CALL alloc_and_configure_domain ( domain_id  = 1           , &
 566                                      grid       = chem_grid   , &
 567                                      parent     = null_domain , &
 568                                      kid        = -1            )
 569
 570 !    else
 571 !     CALL alloc_and_configure_domain ( domain_id  = 2 ,                  &
 572 !                                       grid       = chem_grid ,        &
 573 !                                       parent     = parent_grid ,        &
 574 !                                       kid        = 1                   )
 575 !    endif
 576
 577
 578
 579    CALL       wrf_debug ( 100 , 'wrfchem:input_chem: set pointer for domain 1' )
 580    chgrid => chem_grid
 581
 582    !  Get a list of available file names to try.  This fills up the eligible_file_name
 583    !  array with number_of_eligible_files entries.  This routine issues a nonstandard
 584    !  call (system).
 585
 586    file_counter = 1
 587    need_new_file = .FALSE.
 588
 589    CALL unix_ls ( 'wrf_chem_input' , chem_grid%id )
 590    write(message,'(A,A)') 'number of eligible files ', number_of_eligible_files
 591    CALL wrf_message( TRIM(message) )
 592
 593 !   !  Open the input data (chemin_d01_000000) for reading.
 594 !   CALL wrf_debug          ( 100 , 'subroutine input_chem: calling open_r_dataset for wrfout' )
 595 !   CALL construct_filename ( inpname , 'chemin' , chgrid%id, 2 , 0 , 6 )
 596
 597    CALL construct_filename2a (inpname , chgrid%input_chem_inname, chgrid%id , 2, extract_date)
 598    write(message,'(A,A)') 'Subroutine input_chem: Opening data file ',TRIM(inpname)
 599    CALL wrf_message( TRIM(message) )
 600
 601    CALL open_r_dataset ( fid, TRIM(inpname) , chgrid, config_flags, "DATASET=INPUT", ierr )
 602
 603    IF ( ierr .NE. 0 ) THEN
 604       WRITE( wrf_err_message , FMT='(A,A,A,I8)' ) 'subroutine chemin: error opening ',TRIM(inpname),' for reading ierr=',ierr
 605       CALL WRF_ERROR_FATAL ( wrf_err_message )
 606    ENDIF
 607
 608    !  How many data time levels in the input file?
 609
 610    num_steps = -1
 611    time_loop_max = 0
 612    CALL wrf_debug    ( 100, 'subroutine input_chem: find time_loop_max' )
 613
 614      !  Which times are in this file, and more importantly, are any of them the
 615      !  ones that we want?  We need to loop over times in each files, loop
 616      !  over files.
 617
 618    get_the_right_time : DO
 619 !     CALL ext_ncd_get_next_time ( fid, date_string, status_next_var )
 620       CALL wrf_get_next_time ( fid, date_string, status_next_var )
 621
 622       write(message,'(19A)') 'Subroutine input_chem: file date/time = ',date_string
 623       CALL  wrf_message ( TRIM(message) )
 624
 625       IF ( status_next_var == 0 ) THEN
 626          CALL wrf_debug          ( 100 , 'input_ext_chem_file: calling close_dataset  for ' // TRIM(eligible_file_name(file_counter)) )
 627          CALL close_dataset      ( fid , config_flags , "DATASET=INPUT" )
 628          time_loop_max = time_loop_max + 1
 629             IF ( time_loop_max .GT. number_of_eligible_files ) THEN
 630                WRITE( wrf_err_message , FMT='(A,A,A,I8)' ) 'program input_chem_data: opening too many files'
 631                CALL WRF_ERROR_FATAL ( wrf_err_message )
 632             END IF
 633
 634             IF ( time_loop_max .EQ. number_of_eligible_files ) THEN
 635               num_steps(1) = time_loop_max
 636               num_steps(2) = time_loop_max+1
 637               use_date = date_string
 638               wgt_time1 = dat_time
 639
 640               EXIT get_the_right_time
 641             END IF
 642             CYCLE get_the_right_time
 643       ELSE
 644 !        ELSE IF ( TRIM(date_string) .LT. TRIM(current_date_char) ) THEN
 645 !           CYCLE get_the_right_time
 646 !        ELSE IF ( TRIM(date_string) .EQ. TRIM(current_date_char) ) THEN
 647 !           EXIT get_the_right_time
 648 !        ELSE IF ( TRIM(date_string) .GT. TRIM(current_date_char) ) THEN
 649 !           WRITE( wrf_err_message , FMT='(A,A,A,A,A)' ) 'Found ',TRIM(date_string),' before I found ',TRIM(current_date_char),' .'
 650 !           CALL WRF_ERROR_FATAL ( wrf_err_message )
 651 !        END IF
 652
 653 !    For now, the input date and time MUST match
 654 !
 655 !    Put the time check here and set num_steps
 656
 657         num_steps(1) = time_loop_max
 658         num_steps(2) = time_loop_max+1
 659         use_date = date_string
 660         wgt_time1 = dat_time
 661
 662         EXIT get_the_right_time
 663
 664       ENDIF
 665       if( num_steps(2) == time_loop_max ) then
 666         wgt_time2 = dat_time
 667       endif
 668    END DO get_the_right_time
 669
 670    num_steps(2) = MIN(num_steps(2),time_loop_max)
 671
 672 !  wgt0 = (ext_time  - wgt_time1) / (wgt_time2 - wgt_time1)
 673    wgt0 = 0.
 674
 675 ! Make sure the right date and time for the chemin data has been found
 676
 677    write(message,'(A,A20,A,I9)') 'Subroutine input_chem: use_date, num_steps(1) = ',use_date,num_steps(1)
 678    if( num_steps(1) > 0 ) then
 679       write(message,'(A,A)') 'Subroutine input_chem: extracting data at date/time: ',use_date
 680       CALL  wrf_message ( TRIM( message ) )
 681    else
 682       WRITE( wrf_err_message, FMT='(A)' ) 'subroutine input_chem: error finding chemin date/time #2'
 683       CALL WRF_ERROR_FATAL ( wrf_err_message )
 684    endif
 685
 686    !  There has to be a more elegant way to get to the beginning of the file, but this will do.
 687
 688    CALL close_dataset      ( fid , config_flags , "DATASET=INPUT" )
 689
 690    CALL construct_filename2a (inpname , chgrid%input_chem_inname, chgrid%id , 2, extract_date)
 691    write(message,'(A,A)') 'Subroutine input_chem: Opening data file ',TRIM(inpname)
 692    CALL wrf_message( TRIM(message) )
 693
 694    CALL open_r_dataset     ( fid, TRIM(inpname) , chgrid , config_flags , "DATASET=INPUT", ierr )
 695    IF ( ierr .NE. 0 ) THEN
 696       WRITE( wrf_err_message , FMT='(A,A,A,I8)' ) 'subroutine chemin: error opening ',TRIM(inpname),' for reading ierr=',ierr
 697       CALL WRF_ERROR_FATAL ( wrf_err_message )
 698    ENDIF
 699
 700    !  We know how many time periods to process (right now - all of them), we have the input data
 701    !  (re-)opened, so we begin.
 702
 703    big_time_loop_thingy : DO time_loop = 1 , time_loop_max
 704
 705         CALL wrf_debug          ( 100 , 'input_chem: calling input_history' )
 706         CALL input_history      ( fid , chgrid , config_flags, ierr )
 707         CALL wrf_debug          ( 100 , 'input_chem: back from input_history' )
 708
 709         IF( time_loop .EQ. num_steps(1) ) THEN
 710
 711           ! Get grid dimensions
 712           CALL get_ijk_from_grid (  si_grid ,                        &
 713                                     ids, ide, jds, jde, kds, kde,    &
 714                                     ims, ime, jms, jme, kms, kme,    &
 715                                     ips, ipe, jps, jpe, kps, kpe    )
 716
 717           ! Get scalar grid point heights
 718           ALLOCATE( si_zsigf(ims:ime,kms:kme,jms:jme) )
 719           ALLOCATE( ch_zsigf(ims:ime,kms:kme,jms:jme) )
 720           ALLOCATE(  si_zsig(ims:ime,kms:kme,jms:jme) )
 721           ALLOCATE(  ch_zsig(ims:ime,kms:kme,jms:jme) )
 722 #if ( EM_CORE == 1 )
 723           si_zsigf = (si_grid%ph_1 + si_grid%phb) / grav
 724           ch_zsigf = ( chgrid%ph_1 +  chgrid%phb) / grav
 725 #endif
 726
 727 #if ( NMM_CORE == 1 )
 728
 729    ! Get scalar grid point heights
 730    ALLOCATE(  pint(ips:ipe,kps:kpe,jps:jpe) )
 731    ALLOCATE(  pdsl(ips:ipe,jps:jpe) )
 732
 733        IF(chgrid%sigma.EQ. 1)THEN
 734          do j=jps,jpe
 735          do i=ips,ipe
 736            pdsl(i,j)=chgrid%pd(i,j)
 737          ENDDO
 738          ENDDO
 739        ELSE
 740          do j=jps,jpe
 741          do i=ips,ipe
 742            pdsl(i,j)=chgrid%res(i,j)*chgrid%pd(i,j)
 743          enddo
 744          enddO
 745        ENDIF
 746 !!
 747 !!!!!    SHOULD PINT,Z,W BE REDEFINED IF RESTRT?
 748 !
 749       do j=jps,jpe
 750         do k=kps,kpe
 751         do i=ips,ipe
 752           pint(i,k,j)=si_grid%eta1(k)*chgrid%pdtop+si_grid%eta2(k)*pdsl(i,j)+chgrid%pt
 753           ch_zsigf(i,k,j)=pint(i,k,j)
 754         ENDDO
 755         ENDDO
 756       ENDDO
 757
 758           CALL wrf_debug (0, message)
 759
 760        IF(si_grid%sigma.EQ. 1)THEN
 761          do j=jps,jpe
 762          do i=ips,ipe
 763            pdsl(i,j)=si_grid%pd(i,j)
 764          ENDDO
 765          ENDDO
 766        ELSE
 767          do j=jps,jpe
 768          do i=ips,ipe
 769            pdsl(i,j)=si_grid%res(i,j)*si_grid%pd(i,j)
 770          enddo
 771          enddO
 772        ENDIF
 773 !         write(message,'(1e15.6,i6)') pdsl(ips,jps), si_grid%sigma
 774 !!
 775 !!!!!    SHOULD PINT,Z,W BE REDEFINED IF RESTRT?
 776 !
 777       do j=jps,jpe
 778         do k=kps,kpe
 779         do i=ips,ipe
 780           pint(i,k,j)=si_grid%eta1(k)*si_grid%pdtop+si_grid%eta2(k)*pdsl(i,j)+si_grid%pt
 781 !         if (i.EQ. ips .and. j .EQ. jps) then
 782 !         print *,k,pint(i,k,j),si_grid%eta1(k),si_grid%pdtop,si_grid%eta2(k),pdsl(i,j),si_grid%pt
 783 !         endif
 784           si_zsigf(i,k,j)=pint(i,k,j)
 785         ENDDO
 786         ENDDO
 787       ENDDO
 788
 789 !         write(message,'(4e15.6)') si_zsigf(1,1:4,1)
 790 !         CALL wrf_error_fatal (message)
 791
 792           DEALLOCATE( pint); DEALLOCATE( pdsl )
 793 #endif
 794
 795
 796           do k=1,kde-1
 797             si_zsig(:,k,:) = 0.5 * ( si_zsigf(:,k,:) + si_zsigf(:,k+1,:) )
 798             ch_zsig(:,k,:) = 0.5 * ( ch_zsigf(:,k,:) + ch_zsigf(:,k+1,:) )
 799           enddo
 800           si_zsig(:,kde,:) = 0.5 * ( 3. * si_zsigf(:,kde,:) - si_zsigf(:,kde-1,:) )
 801           ch_zsig(:,kde,:) = 0.5 * ( 3. * ch_zsigf(:,kde,:) - ch_zsigf(:,kde-1,:) )
 802
 803       ! check size of si_grid vs chgrid
 804
 805         IF( size(si_grid%chem,1) .NE. ime-ims+1 .OR.  &
 806             size(si_grid%chem,2) .NE. kme-kms+1 .OR.  &
 807             size(si_grid%chem,3) .NE. jme-jms+1 .OR.  &
 808             size(si_grid%chem,4) .NE. num_chem  ) then
 809
 810           CALL wrf_debug (100, ' SI grid dimensions ' )
 811           write(message,'(4i8.8)') size(si_grid%chem,1),size(si_grid%chem,2), &
 812                                    size(si_grid%chem,3),size(si_grid%chem,4)
 813           CALL wrf_debug (100, message)
 814           CALL wrf_debug (100, ' Input data dimensions ' )
 815           write(message,'(4i8.8)') ime-ims+1,kme-kms+1,jme-jms+1,num_chem
 816           CALL wrf_debug (100, message)
 817           write(wrf_err_message,'(A)') 'ERROR IN MODULE_INPUT_CHEM: bad dimensions in input data '
 818           CALL WRF_ERROR_FATAL ( wrf_err_message )
 819         ENDIF
 820
 821         ! Fill top level to prevent spurious interpolation results (no extrapolation)
 822         chgrid%chem(:,kde,:,:) = chgrid%chem(:,kde-1,:,:)
 823
 824       ! Interpolate the chemistry data to the SI grid (holding place for time interpolation)
 825
 826         call vinterp_chem(ims, ime, jms, jme, kms, kme, kme, num_chem, ch_zsig, si_zsig, &
 827                           chgrid%chem, si_grid%chem, .false.)
 828
 829           if(wgt0 == 0.) EXIT big_time_loop_thingy
 830         ENDIF
 831
 832          IF( time_loop .EQ. num_steps(2) ) THEN
 833
 834 !        ! input chemistry sigma levels
 835 !         ch_zsigf = ( chgrid%ph_1 +  chgrid%phb) / grav
 836 !         do k=1,kde-1
 837 !           ch_zsig(:,k,:) = 0.5 * ( ch_zsigf(:,k,:) + ch_zsigf(:,k+1,:) )
 838 !         enddo
 839 !         ch_zsig(:,kde,:) = 0.5 * ( 3. * ch_zsigf(:,kde,:) - ch_zsigf(:,kde-1,:) )
 840
 841 !       ! Fill top level to prevent spurious interpolation results (no extrapolation)
 842 !         chgrid%chem(:,kde,:,:) = chgrid%chem(:,kde-1,:,:)
 843
 844 !     ! Interpolate the chemistry data to the temp chgrid2 structure
 845
 846 !         call vinterp_chem(ims, ime, jms, jme, kms, kme, kme, num_chem, ch_zsig, si_zsig, &
 847 !                           chgrid%chem, chgrid2%chem, .false.)
 848
 849 !     ! use linear interpolation in time to get new chem arrays
 850 !         si_grid%chem = (1. - wgt0) * si_grid%chem + &
 851 !                                wgt0  * chgrid2%chem
 852
 853           DEALLOCATE( si_zsigf); DEALLOCATE( si_zsig )
 854           DEALLOCATE( ch_zsigf); DEALLOCATE( ch_zsig )
 855
 856           EXIT big_time_loop_thingy
 857         ENDIF
 858    END DO big_time_loop_thingy
 859
 860 !  Check for errors in chemin data set
 861
 862    do l=2,num_chem
 863    do j=jds,jde
 864    do k=kds,kde
 865    do i=ids,ide
 866      si_grid%chem(i,k,j,l) = MAX(si_grid%chem(i,k,j,l),epsilc)
 867    enddo
 868    enddo
 869    enddo
 870    enddo
 871
 872
 873 !  Close chemin data set
 874    CALL close_dataset ( fid , config_flags , "DATASET=INPUT" )
 875
 876 ! free memory
 877    CALL domain_destroy( chem_grid )
 878
 879    CALL wrf_debug ( 100,' input_ext_chem_data: exit subroutine ')
 880
 881    RETURN
 882
 883   END SUBROUTINE input_ext_chem_file
 884 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 885   SUBROUTINE vinterp_chem(nx1, nx2, ny1, ny2, nz1, nz_in, nz_out, nch, z_in, z_out, &
 886                  data_in, data_out, extrapolate)
 887
 888     ! Interpolates columns of chemistry data from one set of height surfaces to
 889     ! another
 890
 891     INTEGER, INTENT(IN)                :: nx1, nx2
 892     INTEGER, INTENT(IN)                :: ny1, ny2
 893     INTEGER, INTENT(IN)                :: nz1
 894     INTEGER, INTENT(IN)                :: nz_in
 895     INTEGER, INTENT(IN)                :: nz_out
 896     INTEGER, INTENT(IN)                :: nch
 897     REAL, INTENT(IN)                   :: z_in (nx1:nx2,nz1:nz_in ,ny1:ny2)
 898     REAL, INTENT(IN)                   :: z_out(nx1:nx2,nz1:nz_out,ny1:ny2)
 899     REAL, INTENT(IN)                   :: data_in (nx1:nx2,nz1:nz_in ,ny1:ny2,nch)
 900     REAL, INTENT(OUT)                  :: data_out(nx1:nx2,nz1:nz_out,ny1:ny2,nch)
 901     LOGICAL, INTENT(IN)                :: extrapolate
 902
 903     INTEGER                            :: i,j,l
 904     INTEGER                            :: k,kk
 905     REAL                               :: desired_z
 906     REAL                               :: dvaldz
 907     REAL                               :: wgt0
 908
 909 !   Loop over the number of chemical species
 910     chem_loop: DO l = 2, nch
 911
 912       data_out(:,:,:,l) = -99999.9
 913
 914       DO j = ny1, ny2
 915         DO i = nx1, nx2
 916
 917           output_loop: DO k = nz1, nz_out
 918 #if ( EM_CORE == 1 )
 919
 920             desired_z = z_out(i,k,j)
 921             IF (desired_z .LT. z_in(i,1,j)) THEN
 922
 923               IF ((desired_z - z_in(i,1,j)).LT. 0.0001) THEN
 924                  data_out(i,k,j,l) = data_in(i,1,j,l)
 925               ELSE
 926                 IF (extrapolate) THEN
 927                   ! Extrapolate downward because desired height level is below
 928                   ! the lowest level in our input data.  Extrapolate using simple
 929                   ! 1st derivative of value with respect to height for the bottom 2
 930                   ! input layers.
 931
 932                   ! Add a check to make sure we are not using the gradient of
 933                   ! a very thin layer
 934
 935                   IF ( (z_in(i,1,j) - z_in(i,2,j)) .GT. 0.001) THEN
 936                     dvaldz = (data_in(i,1,j,l) - data_in(i,2,j,l)) / &
 937                               (z_in(i,1,j)  - z_in(i,2,j) )
 938                   ELSE
 939                     dvaldz = (data_in(i,1,j,l) - data_in(i,3,j,l)) / &
 940                               (z_in(i,1,j)  - z_in(i,3,j) )
 941                   ENDIF
 942                   data_out(i,k,j,l) = MAX( data_in(i,1,j,l) + &
 943                                 dvaldz * (desired_z-z_in(i,1,j)), 0.)
 944                 ELSE
 945                   data_out(i,k,j,l) = data_in(i,1,j,l)
 946                 ENDIF
 947               ENDIF
 948             ELSE IF (desired_z .GT. z_in(i,nz_in,j)) THEN
 949               IF ( (z_in(i,nz_in,j) - desired_z) .LT. 0.0001) THEN
 950                  data_out(i,k,j,l) = data_in(i,nz_in,j,l)
 951               ELSE
 952                 IF (extrapolate) THEN
 953                   ! Extrapolate upward
 954                   IF ( (z_in(i,nz_in-1,j)-z_in(i,nz_in,j)) .GT. 0.0005) THEN
 955                     dvaldz = (data_in(i,nz_in,j,l) - data_in(i,nz_in-1,j,l)) / &
 956                                (z_in(i,nz_in,j)  - z_in(i,nz_in-1,j))
 957                   ELSE
 958                     dvaldz = (data_in(i,nz_in,j,l) - data_in(i,nz_in-2,j,l)) / &
 959                                (z_in(i,nz_in,j)  - z_in(i,nz_in-2,j))
 960                   ENDIF
 961                   data_out(i,k,j,l) =  MAX( data_in(i,nz_in,j,l) + &
 962                            dvaldz * (desired_z-z_in(i,nz_in,j)), 0.)
 963                 ELSE
 964                   data_out(i,k,j,l) = data_in(i,nz_in,j,l)
 965                 ENDIF
 966               ENDIF
 967             ELSE
 968               ! We can trap between two levels and linearly interpolate
 969
 970               input_loop:  DO kk = 1, nz_in-1
 971                 IF (desired_z .EQ. z_in(i,kk,j) )THEN
 972                   data_out(i,k,j,l) = data_in(i,kk,j,l)
 973                   EXIT input_loop
 974                 ELSE IF (desired_z .EQ. z_in(i,kk+1,j) )THEN
 975                   data_out(i,k,j,l) = data_in(i,kk+1,j,l)
 976                   EXIT input_loop
 977                 ELSE IF ( (desired_z .GT. z_in(i,kk,j)) .AND. &
 978                           (desired_z .LT. z_in(i,kk+1,j)) ) THEN
 979                   wgt0 = (desired_z - z_in(i,kk+1,j)) / &
 980                          (z_in(i,kk,j)-z_in(i,kk+1,j))
 981                   data_out(i,k,j,l) = MAX( wgt0*data_in(i,kk,j,l) + &
 982                                     (1.-wgt0)*data_in(i,kk+1,j,l), 0.)
 983                   EXIT input_loop
 984                 ENDIF
 985               ENDDO input_loop
 986
 987             ENDIF
 988 #endif
 989 #if ( NMM_CORE == 1 )
 990
 991             desired_z = z_out(i,k,j)
 992             IF (desired_z .GT. z_in(i,1,j)) THEN
 993
 994               IF ((desired_z - z_in(i,1,j)).GT. 0.0001) THEN
 995                  data_out(i,k,j,l) = data_in(i,1,j,l)
 996               ELSE
 997                 IF (extrapolate) THEN
 998                   ! Extrapolate upward because desired pressure level is above
 999                   ! the highest level in our input data.  Extrapolate using simple
1000                   ! 1st derivative of value with respect to height for the bottom 2
1001                   ! input layers.
1002
1003                   ! Add a check to make sure we are not using the gradient of
1004                   ! a very thin layer
1005
1006                   IF ( (z_in(i,1,j) - z_in(i,2,j)) .LT. 0.001) THEN
1007                     dvaldz = (data_in(i,2,j,l) - data_in(i,1,j,l)) / &
1008                               (z_in(i,2,j)  - z_in(i,1,j) )
1009                   ELSE
1010                     dvaldz = (data_in(i,3,j,l) - data_in(i,1,j,l)) / &
1011                               (z_in(i,3,j)  - z_in(i,1,j) )
1012                   ENDIF
1013                   data_out(i,k,j,l) = MAX( data_in(i,1,j,l) + &
1014                                  dvaldz * (desired_z-z_in(i,1,j)), 0.)
1015                 ELSE
1016                   data_out(i,k,j,l) = data_in(i,1,j,l)
1017                 ENDIF
1018               ENDIF
1019             ELSE IF (desired_z .LT. z_in(i,nz_in,j)) THEN
1020               IF ( (z_in(i,nz_in,j) - desired_z) .LT. 0.0001) THEN
1021                  data_out(i,k,j,l) = data_in(i,nz_in,j,l)
1022               ELSE
1023                 IF (extrapolate) THEN
1024                   ! Extrapolate upward
1025                   IF ( (z_in(i,nz_in-1,j)-z_in(i,nz_in,j)) .LT. 0.0005) THEN
1026                     dvaldz = (data_in(i,nz_in,j,l) - data_in(i,nz_in-1,j,l)) / &
1027                                (z_in(i,nz_in,j)  - z_in(i,nz_in-1,j))
1028                   ELSE
1029                     dvaldz = (data_in(i,nz_in,j,l) - data_in(i,nz_in-2,j,l)) / &
1030                                (z_in(i,nz_in,j)  - z_in(i,nz_in-2,j))
1031                   ENDIF
1032                   data_out(i,k,j,l) =  MAX( data_in(i,nz_in,j,l) + &
1033                            dvaldz * (z_in(i,nz_in,j) - desired_z), 0.)
1034                 ELSE
1035                   data_out(i,k,j,l) = data_in(i,nz_in,j,l)
1036                 ENDIF
1037               ENDIF
1038             ELSE
1039               ! We can trap between two levels and linearly interpolate
1040
1041               input_loop:  DO kk = 1, nz_in-1
1042                 IF (desired_z .EQ. z_in(i,kk,j) )THEN
1043                   data_out(i,k,j,l) = data_in(i,kk,j,l)
1044                   EXIT input_loop
1045                 ELSE IF ( (desired_z .LT. z_in(i,kk,j)) .AND. &
1046                           (desired_z .GT. z_in(i,kk+1,j)) ) THEN
1047                   wgt0 = (desired_z - z_in(i,kk+1,j)) / &
1048                          (z_in(i,kk,j)-z_in(i,kk+1,j))
1049                   data_out(i,k,j,l) = MAX( wgt0*data_in(i,kk,j,l) + &
1050                                     (1.-wgt0)*data_in(i,kk+1,j,l), 0.)
1051                   EXIT input_loop
1052                 ENDIF
1053               ENDDO input_loop
1054
1055             ENDIF
1056 #endif
1057           ENDDO output_loop
1058         ENDDO
1059       ENDDO
1060     ENDDO chem_loop
1061
1062     RETURN
1063   END SUBROUTINE vinterp_chem
1064 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1065 SUBROUTINE input_chem_profile (si_grid)
1066
1067    IMPLICIT NONE
1068
1069    TYPE(domain)           ::  si_grid
1070
1071    INTEGER :: i , j , k, &
1072               ids, ide, jds, jde, kds, kde,    &
1073               ims, ime, jms, jme, kms, kme,    &
1074               ips, ipe, jps, jpe, kps, kpe
1075    INTEGER :: fid, ierr, numgas
1076    INTEGER :: debug_level
1077
1078    REAL, ALLOCATABLE, DIMENSION(:,:,:) :: si_zsigf, si_zsig
1079
1080 #if ( ! NMM_CORE == 0 )
1081    REAL, ALLOCATABLE, DIMENSION(:,:,:) :: pint
1082    REAL, ALLOCATABLE, DIMENSION(:,:) :: pdsl
1083 #endif
1084
1085    CHARACTER (LEN=80) :: inpname, message
1086
1087    write(message,'(A)') 'Subroutine input_chem_profile: '
1088    CALL  wrf_message ( TRIM(message) )
1089
1090    !  And here is an instance of using the information in the NAMELIST.
1091
1092    CALL nl_get_debug_level ( 1,debug_level )
1093    CALL set_wrf_debug_level ( debug_level )
1094
1095    ! Get grid dimensions
1096    CALL get_ijk_from_grid (  si_grid ,                        &
1097                              ids, ide, jds, jde, kds, kde,    &
1098                              ims, ime, jms, jme, kms, kme,    &
1099                              ips, ipe, jps, jpe, kps, kpe    )
1100
1101    ! Get scalar grid point heights
1102    ALLOCATE( si_zsigf(ims:ime,kms:kme,jms:jme) )
1103    ALLOCATE(  si_zsig(ims:ime,kms:kme,jms:jme) )
1104
1105 #if ( ! EM_CORE == 0 )
1106      write(message,'(A)') 'WRF_EM_CORE  '
1107      si_zsigf = (si_grid%ph_1 + si_grid%phb) / grav
1108 #endif
1109 #if ( ! NMM_CORE == 0 )
1110    ! Get scalar grid point heights
1111    ALLOCATE(  pint(ims:ime,kms:kme,jms:jme) )
1112    ALLOCATE(  pdsl(ims:ime,jms:jme) )
1113
1114      write(message,'(A)') 'WRF_NMM_CORE  '
1115      CALL  wrf_message ( message )
1116
1117        IF(si_grid%sigma.EQ. 1)THEN
1118          do j=jps,jpe
1119          do i=ips,ipe
1120            pdsl(i,j)=si_grid%pd(i,j)
1121          ENDDO
1122          ENDDO
1123        ELSE
1124          do j=jps,jpe
1125          do i=ips,ipe
1126            pdsl(i,j)=si_grid%res(i,j)*si_grid%pd(i,j)
1127          enddo
1128          enddO
1129        ENDIF
1130 !!
1131 !!***
1132 !!
1133 !!
1134 !!!!!    SHOULD PINT,Z,W BE REDEFINED IF RESTRT?
1135 !
1136 !     print *,' ips=',ips,' ipe=',ipe
1137 !     print *,' jps=',jps,' jpe=',jpe
1138 !     print *,' kps=',kps,' kpe=',kpe
1139 !     print *,' sigma=',si_grid%sigma
1140 !     print *,' pdtop=',si_grid%pdtop,' pt=',si_grid%pt
1141
1142       do j=jps,jpe
1143         do k=kps,kpe
1144         do i=ips,ipe
1145           pint(i,k,j)=si_grid%eta1(k)*si_grid%pdtop+si_grid%eta2(k)*pdsl(i,j)+si_grid%pt
1146           si_zsigf(i,k,j)=pint(i,k,j)
1147         ENDDO
1148         ENDDO
1149       ENDDO
1150 !       do k=kps,kpe
1151 !          print *,k,pint(1,k,1),si_grid%eta1(k),si_grid%pdtop,si_grid%eta2(k),pdsl(1,1),si_grid%pt
1152 !       enddo
1153 !
1154 !  si_zsigf = si_grid%z
1155 #endif
1156
1157 !  si_zsigf = (si_grid%ph_1 + si_grid%phb) / grav
1158
1159    do k=1,kde-1
1160      si_zsig(:,k,:) = 0.5 * ( si_zsigf(:,k,:) + si_zsigf(:,k+1,:) )
1161    enddo
1162    si_zsig(:,kde,:) = 0.5 * ( 3. * si_zsigf(:,kde,:) - si_zsigf(:,kde-1,:) )
1163
1164    ! Determine the index of the last gas species
1165    numgas = get_last_gas(si_grid%chem_opt)
1166
1167    ! Setup the cross reference mappings between the default profiles and
1168    ! the gas mechanism species (wig, 2-May-2007)
1169    call setup_gasprofile_maps(si_grid%chem_opt, numgas)
1170
1171    ! An alternative ozone profile option for initialization
1172    if( si_grid%gas_ic_opt == GAS_IC_PNNL ) &
1173         call gasprofile_init_pnnl( si_grid%chem_opt )
1174
1175    ! Interpolate the chemistry data to the SI grid. These values should typically
1176    ! be set to match the values in bdy_chem_value_tracer so that the boundaries
1177    ! and interior match each other.
1178    IF ( si_grid%chem_opt == CHEM_TRACER ) THEN
1179       si_grid%chem(ims:ime,kms:kme,jms:jme,1:numgas) = 0.0001
1180 !     si_grid%chem(ims:ime,kms:kme,jms:jme,p_so2) = 0.0001
1181       si_grid%chem(ims:ime,kms:kme,jms:jme,p_co ) = 0.08
1182    ELSE IF ( si_grid%chem_opt == CHEM_TRACE2 ) THEN
1183       si_grid%chem(ims:ime,kms:kme,jms:jme,p_TRACER_1 ) = 0.08
1184    ELSE IF ( si_grid%chem_opt == GOCART_SIMPLE ) THEN
1185       si_grid%chem(ims:ime,kms:kme,jms:jme,1:num_chem) = 1.e-12
1186       si_grid%chem(ims:ime,kms:kme,jms:jme,p_so2) = 1.e-6
1187       si_grid%chem(ims:ime,kms:kme,jms:jme,p_sulf) = 3.e-6
1188       si_grid%chem(ims:ime,kms:kme,jms:jme,p_dms) = 1.e-6
1189       si_grid%chem(ims:ime,kms:kme,jms:jme,p_msa) = 1.e-6
1190       si_grid%chem(ims:ime,kms:kme,jms:jme,p_bc1 ) = 1.e-2
1191       si_grid%chem(ims:ime,kms:kme,jms:jme,p_bc2 ) = 1.e-2
1192       si_grid%chem(ims:ime,kms:kme,jms:jme,p_oc1 ) = 1.e-2
1193       si_grid%chem(ims:ime,kms:kme,jms:jme,p_oc2 ) = 1.e-2
1194       si_grid%chem(ims:ime,kms:kme,jms:jme,p_p25 ) = 1.
1195    ELSE
1196       CALL make_chem_profile (ims, ime, jms, jme, kms, kme, num_chem, numgas, &
1197                               si_grid%chem_opt, si_zsig, si_grid%chem)
1198    END IF
1199
1200    CALL wrf_debug       ( 100,' input_chem_profile: exit subroutine ')
1201
1202    DEALLOCATE( si_zsigf ); DEALLOCATE( si_zsig )
1203 #if ( ! NMM_CORE == 0 )
1204    DEALLOCATE( pdsl ); DEALLOCATE( pint )
1205 #endif
1206    RETURN
1207
1208   END SUBROUTINE input_chem_profile
1209 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1210   SUBROUTINE make_chem_profile ( nx1, nx2, ny1, ny2, nz1, nz2, nch, numgas, &
1211                                  chem_opt, zgrid, chem )
1212     IMPLICIT NONE
1213
1214     INTEGER, INTENT(IN) :: nx1, ny1, nz1
1215     INTEGER, INTENT(IN) :: nx2, ny2, nz2
1216     INTEGER, INTENT(IN) :: nch, numgas, chem_opt
1217     !REAL, INTENT(IN), DIMENSION(nx1:nx2,nz1:nz2,ny1:ny2) :: zgrid
1218     REAL, DIMENSION(nx1:nx2,nz1:nz2,ny1:ny2) :: zgrid
1219
1220     CHARACTER (LEN=80) :: message
1221     INTEGER :: i, j, k, l, is
1222
1223     REAL, DIMENSION(nx1:nx2,nz1:kx ,ny1:ny2,lo+1):: chprof
1224     REAL, DIMENSION(nx1:nx2,nz1:kx ,ny1:ny2)     :: zprof
1225
1226     REAL, DIMENSION(nx1:nx2,nz1:nz2,ny1:ny2,nch) :: chem
1227     REAL, DIMENSION(nx1:nx2,nz1:nz2,ny1:ny2,lo ) :: stor
1228
1229     REAL :: hc358 !wig, 2-May-2007
1230     REAL :: olit
1231
1232 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1233 !
1234 ! Check the number of species
1235
1236      if( nch .NE. num_chem) then
1237        message = ' Input_chem_profile: wrong number of chemical species'
1238 !       CALL WRF_ERROR_FATAL ( message )
1239      endif
1240
1241       ! Vertically flip the chemistry data as it is given top down and
1242       ! heights are bottom up. Fill temp 3D chemical and profile array,
1243       ! keep chem slot 1 open as vinterp_chem assumes there is no data.
1244       DO j=ny1,ny2
1245       DO k=  1,kx
1246       DO i=nx1,nx2
1247          chprof(i,k,j,2:lo+1) = xl(1:lo,kx-k+1)
1248          zprof(i,k,j) = 0.5*(zfa(k)+zfa(k+1))
1249       ENDDO
1250       ENDDO
1251       ENDDO
1252 !
1253 ! return xl to previous value for next time...
1254 !   34 chemicals (lo), 16 vertical levels (kx)
1255 !     DO i=lo-6,lo
1256 !       xl(i,1:kx)=xl(i,1:kx)*dens(1:kx)
1257 !     ENDDO
1258
1259 ! Change number concentrations to mixing ratios for short-lived NALROM species
1260       do k=1,kx
1261          chprof(:,k,:,lo-5:lo+1) = chprof(:,k,:,lo-5:lo+1)/dens(k)
1262       end do
1263
1264       ! Interpolate temp 3D chemical and profile array to WRF grid
1265       call vinterp_chem(nx1, nx2, ny1, ny2, nz1, kx, nz2, lo, zprof, zgrid, &
1266                           chprof, chem, .false.)
1267
1268       ! place interpolated data into temp storage array
1269       stor(nx1:nx2,nz1:nz2,ny1:ny2,1:lo) = chem(nx1:nx2,nz1:nz2,ny1:ny2,2:lo+1)
1270
1271       ! Here is where the chemistry profile is constructed
1272       !chem(:,:,:,1) = stor(:,:,:,1) * 0.
1273       chem(nx1:nx2,nz1:nz2,ny1:ny2,1) = -999.
1274
1275 !      DO l=2,nch
1276       DO l=2,numgas
1277         is=iref(l-1)
1278         DO j=ny1,ny2
1279         DO k=nz1,nz2
1280         DO i=nx1,nx2
1281           chem(i,k,j,l)=fracref(l-1)*stor(i,k,j,is)*1.E6
1282         ENDDO
1283         ENDDO
1284         ENDDO
1285       ENDDO
1286 !
1287 ! For CBMZ, we need to construct PAR based on a combination of other
1288 ! species. This cannot be done with the looping construct above so
1289 ! we have to treat it separately. (wig, 2-May-2007)
1290 !
1291       SELECT CASE(chem_opt)
1292       CASE (CBMZ,CBMZ_BB,CBMZ_BB_KPP,CBMZ_MOSAIC_4BIN,CBMZ_MOSAIC_8BIN, &
1293             CBMZ_MOSAIC_4BIN_AQ,CBMZ_MOSAIC_8BIN_AQ)
1294          do j = ny1,ny2
1295          do k = nz1,nz2
1296          do i = nx1,nx2
1297             !Construct the sum of the profiles for hc3, hc5, & hc8
1298             hc358 = ( 2.9*fracref(numgas+1)*stor(i,k,j,iref(numgas+1)) &
1299                      +4.8*fracref(numgas+2)*stor(i,k,j,iref(numgas+2)) &
1300                      +7.9*fracref(numgas+3)*stor(i,k,j,iref(numgas+3)) &
1301                     )*1.E6
1302             chem(i,k,j,p_par) =                                    &
1303                  0.4*chem(i,k,j,p_ald) + hc358                     &
1304                  + 0.9*chem(i,k,j,p_ket) + 2.8*chem(i,k,j,p_oli)   &
1305                  + 1.8*chem(i,k,j,p_olt) + 1.0*chem(i,k,j,p_ora2)
1306          end do
1307          end do
1308          end do
1309
1310
1311
1312       CASE (CBM4_KPP)
1313          do j = ny1,ny2
1314          do k = nz1,nz2
1315          do i = nx1,nx2
1316             !Construct the sum of the profiles for hc3, hc5, & hc8
1317             hc358 = ( 2.9*fracref(numgas+1)*stor(i,k,j,iref(numgas+1)) &
1318                      +4.8*fracref(numgas+2)*stor(i,k,j,iref(numgas+2)) &
1319                      +7.9*fracref(numgas+3)*stor(i,k,j,iref(numgas+3)) &
1320                     )*1.E6
1321             olit = ( 0.9*fracref(numgas+4)*stor(i,k,j,iref(numgas+4)) &
1322                      +2.8*fracref(numgas+5)*stor(i,k,j,iref(numgas+5)) &
1323                      +1.8*fracref(numgas+6)*stor(i,k,j,iref(numgas+6)) &
1324                      +1.0*fracref(numgas+7)*stor(i,k,j,iref(numgas+7)) &
1325                     )*1.E6
1326             chem(i,k,j,p_par) =  0.4*chem(i,k,j,p_ald2) + hc358  + olit
1327          end do
1328          end do
1329          end do
1330       END SELECT
1331
1332       RETURN
1333   END SUBROUTINE make_chem_profile
1334 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1335 !
1336 ! this is a kludge routine as of now....
1337 !
1338   SUBROUTINE bdy_chem_value_sorgam (chem, z, nch, config_flags, &
1339                                       alt,convfac,g)
1340   USE module_data_sorgam
1341
1342     IMPLICIT NONE
1343
1344     REAL,    intent(OUT)  :: chem
1345     REAL,    intent(IN)   :: z          ! 3D height array
1346     INTEGER, intent(IN)   :: nch        ! index number of chemical species
1347     REAL,  INTENT(IN   ) ::   alt, convfac
1348     real, INTENT (IN) :: g
1349     TYPE (grid_config_rec_type), intent(in) :: config_flags
1350
1351     INTEGER :: i, k, l
1352     REAL, DIMENSION(lo+1,1:kx):: cprof  ! chemical profile, diff. index order
1353
1354     REAL, DIMENSION(1:kx):: zprof
1355     REAL, DIMENSION(lo ) :: stor
1356     REAL                 :: wgt0
1357
1358     real :: chemsulf_radm,chem_so4aj,chem_so4ai
1359      real tempfac
1360       REAL :: splitfac
1361                         !between gas and aerosol phase
1362
1363 !factor for splitting initial conc. of SO4
1364 !3rd moment i-mode [3rd moment/m^3]
1365       REAL :: m3nuc
1366 !3rd MOMENT j-mode [3rd moment/m^3]
1367       REAL :: m3acc
1368 !       REAL ESN36
1369       REAL :: m3cor
1370       DATA splitfac/.98/
1371
1372 !
1373 ! method for bc calculation is determined by aer_bc_opt
1374 !
1375        if (config_flags%aer_bc_opt == AER_BC_PNNL) then
1376            call sorgam_set_aer_bc_pnnl( chem, z, nch, config_flags )
1377            return
1378        else if (config_flags%aer_bc_opt == AER_BC_DEFAULT) then
1379            continue
1380        else
1381            call wrf_error_fatal(   &
1382                "bdy_chem_value_sorgam -- unable to parse aer_bc_opt" )
1383        end if
1384
1385 ! do default calculation of sorgam aerosol bc values
1386        chem=conmin
1387 !      tempfac=(t+t0)*((p+pb)/p1000mb)**rcp
1388 !      convfac=(p+pb)/rgasuniv/tempfac
1389 !
1390 !--- units for advection....
1391 !
1392        if(nch.eq.p_nu0)chem=1.e8*alt
1393        if(nch.eq.p_ac0)chem=1.e8*alt
1394        if(nch.eq.p_nh4aj)chem=10.e-5*alt
1395        if(nch.eq.p_nh4ai)chem=10.e-5*alt
1396        if(nch.eq.p_no3aj)chem=10.e-5*alt
1397        if(nch.eq.p_no3ai)chem=10.e-5*alt
1398 !
1399 ! recalculate sulf profile for aerosols
1400 !
1401      if   ( nch .eq. p_so4aj.or.nch.eq.p_so4ai                        &
1402         .or.nch .eq. p_nu0  .or.nch.eq.p_ac0                          &
1403         .or.nch .eq. p_corn                    ) then
1404
1405       ! Vertically flip the chemistry data as it is given top down
1406       !     and heights in zfa are bottom up
1407       ! Fill chemical profile array cprof
1408       ! Keep chem slot 1 open as vinterp_chem assumes there is no data
1409       !     (this isn't really needed in this subr)
1410       ! Convert species 28-34 (lo-6:lo) from (molecules/cm3) to (mol/mol)
1411       DO k = 1,kx
1412         zprof(k) = 0.5*(zfa_bdy(k)+zfa_bdy(k+1))
1413         DO l = 1,lo-7
1414            cprof(l+1,k) = xl(l,kx+1-k)
1415         END DO
1416 ! Fix number concentrations to mixing ratios for short-lived NALROM species
1417         DO l = lo-6,lo
1418             cprof(l+1,k) = xl(l,kx+1-k)/dens(kx+1-k)
1419         ENDDO
1420       ENDDO
1421
1422       ! Interpolate temp 1D chemical profile array to WRF field
1423       IF (z .LT. zprof(1)) THEN
1424         stor(1:lo) = cprof(2:lo+1,1)
1425       ELSE IF (z .GT. zprof(kx)) THEN
1426         stor(1:lo) = cprof(2:lo+1,kx)
1427       ELSE
1428         ! We can trap between two levels and linearly interpolate
1429         input_loop:  DO k = 1, kx-1
1430           IF (z .EQ. zprof(k) )THEN
1431             stor(1:lo) = cprof(2:lo+1,k)
1432             EXIT input_loop
1433           ELSE IF ( (z .GT. zprof(k)) .AND. &
1434                     (z .LT. zprof(k+1)) ) THEN
1435             wgt0 = (z   - zprof(k+1)) / &
1436                    (zprof(k) - zprof(k+1))
1437             stor(1:lo) = MAX( wgt0 *cprof(2:lo+1,k  ) + &
1438                           (1.-wgt0)*cprof(2:lo+1,k+1), 0.)
1439             EXIT input_loop
1440           ENDIF
1441         ENDDO input_loop
1442       ENDIF
1443
1444       ! Here is where the chemistry value is constructed
1445       chemsulf_radm = fracref(p_sulf-1)*stor( iref(p_sulf-1) )*1.E6
1446 !
1447 ! now have sulf
1448 !
1449        chem_so4aj=chemsulf_radm*CONVFAC*MWSO4*splitfac*so4vaptoaer
1450        chem_so4ai=chemsulf_radm*CONVFAC*MWSO4*(1.-splitfac)*so4vaptoaer
1451        if(nch.eq.p_so4aj)chem=chem_so4aj*alt
1452        if(nch.eq.p_so4ai)chem=chem_so4ai*alt
1453        m3nuc=so4fac*chem_so4ai+conmin*(nh4fac+no3fac+orgfac*9+2*anthfac)
1454        m3acc=so4fac*chem_so4aj+conmin*(nh4fac+no3fac+orgfac*9+2*anthfac)
1455        m3cor=conmin*(soilfac+seasfac+anthfac)
1456 !
1457 ! compute values for aerosol input data
1458 !
1459        if(nch.eq.p_nu0.or.nch.eq.p_ac0.or.nch.eq.p_corn)then
1460          xxlsgn = log(sginin)
1461         xxlsga = log(sginia)
1462         xxlsgc = log(sginic)
1463
1464         l2sginin = xxlsgn**2
1465         l2sginia = xxlsga**2
1466         l2sginic = xxlsgc**2
1467
1468         en1 = exp(0.125*l2sginin)
1469         ea1 = exp(0.125*l2sginia)
1470         ec1 = exp(0.125*l2sginic)
1471
1472
1473         esn04 = en1**4
1474         esa04 = ea1**4
1475         esc04 = ec1**4
1476
1477         esn05 = esn04*en1
1478         esa05 = esa04*ea1
1479
1480         esn08 = esn04*esn04
1481         esa08 = esa04*esa04
1482         esc08 = esc04*esc04
1483
1484         esn09 = esn04*esn05
1485         esa09 = esa04*esa05
1486
1487         esn12 = esn04*esn04*esn04
1488         esa12 = esa04*esa04*esa04
1489         esc12 = esc04*esc04*esc04
1490
1491         esn16 = esn08*esn08
1492         esa16 = esa08*esa08
1493         esc16 = esc08*esc08
1494
1495         esn20 = esn16*esn04
1496         esa20 = esa16*esa04
1497         esc20 = esc16*esc04
1498
1499         esn24 = esn12*esn12
1500         esa24 = esa12*esa12
1501         esc24 = esc12*esc12
1502
1503         esn25 = esn16*esn09
1504         esa25 = esa16*esa09
1505
1506         esn28 = esn20*esn08
1507         esa28 = esa20*esa08
1508         esc28 = esc20*esc08
1509
1510
1511         esn32 = esn16*esn16
1512         esa32 = esa16*esa16
1513         esc32 = esc16*esc16
1514
1515         esn36 = esn16*esn20
1516         esa36 = esa16*esa20
1517         esc36 = esc16*esc20
1518        endif
1519 !
1520 ! Units are something like number concentration
1521 !
1522        if(nch.eq.p_nu0)chem=m3nuc/((dginin**3)*esn36)*alt
1523        if(nch.eq.p_ac0)chem=m3acc/((dginia**3)*esa36)*alt
1524        if(nch.eq.p_corn)chem=m3cor/((dginic**3)*esc36)*alt
1525      endif
1526
1527
1528   END SUBROUTINE bdy_chem_value_sorgam
1529
1530   SUBROUTINE bdy_chem_value_gocart ( chem, nch )
1531
1532 ! This subroutine is called to set the boundary values of chemistry
1533
1534     IMPLICIT NONE
1535
1536     REAL,    intent(OUT)  :: chem
1537     INTEGER, intent(IN)   :: nch        ! index number of chemical species
1538 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1539     if( nch == p_so2  ) then
1540        chem = 5.e-6
1541     else if( nch == p_sulf ) then
1542        chem = 3.e-6
1543     else if( nch == p_dms ) then
1544        chem = 1.e-6
1545     else if( nch == p_msa ) then
1546        chem = 1.e-6
1547     else if( nch == p_bc1 ) then
1548        chem = 1.e-2
1549     else if( nch == p_bc2 ) then
1550        chem = 1.e-2
1551     else if( nch == p_oc1 ) then
1552        chem = 1.e-2
1553     else if( nch == p_oc2 ) then
1554        chem = 1.e-2
1555     else if( nch == p_p25 ) then
1556        chem = 1.
1557     else
1558        chem = 1.e-12
1559     end if
1560
1561   END SUBROUTINE bdy_chem_value_GOCART
1562   SUBROUTINE bdy_chem_value_tracer ( chem, nch )
1563
1564 ! This subroutine is called to set the boundary values of chemistry
1565 ! species when chem_opt==CHEM_TRACER. Typically, the boundary values
1566 ! here should be set to match those in input_chem_profile so that the
1567 ! interior and boundary values are the same.
1568 ! William.Gustafson@pnl.gov; 16-Jun-2005
1569
1570     IMPLICIT NONE
1571
1572     REAL,    intent(OUT)  :: chem
1573     INTEGER, intent(IN)   :: nch        ! index number of chemical species
1574 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1575
1576     if( nch .ne. p_co  ) then
1577        chem = 0.0001
1578     else if( nch == p_co ) then
1579        chem = 0.08
1580     else
1581        chem = conmin
1582     end if
1583     if( nch .eq. p_tracer_1  ) then
1584        chem = 0.08
1585     endif
1586
1587   END SUBROUTINE bdy_chem_value_tracer
1588 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1589
1590 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1591   SUBROUTINE bdy_chem_value_racm ( chem, z, nch, numgas,p_co2 )
1592
1593     IMPLICIT NONE
1594
1595     REAL,    intent(OUT)  :: chem
1596     REAL,    intent(IN)   :: z          ! 3D height array
1597     INTEGER, intent(IN)   :: nch,p_co2  ! index number of chemical species
1598     INTEGER, intent(IN)   :: numgas     ! index number of last gas species
1599
1600     INTEGER :: i, k, irefcur
1601
1602     REAL, DIMENSION(kx):: cprof         ! chemical profile, diff. index order
1603
1604     REAL, DIMENSION(1:kx):: zprof
1605     REAL                 :: stor
1606     REAL                 :: wgt0
1607
1608     CHARACTER (LEN=80) :: message
1609 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1610
1611 ! Check the number of species
1612 !     if((nch-1).gt.logg)return
1613      if (nch.eq.p_co2)then
1614        chem=370.
1615        return
1616      endif
1617      if (nch.eq.p_co2+1)then
1618        chem=1.7
1619        return
1620      endif
1621      if (nch.ge.p_co2+2)return
1622
1623
1624 !    if( nch .GT. logg+1) then
1625      if( nch .GT. numgas) then
1626        message = ' Input_chem_profile: wrong number of chemical species'
1627        return
1628 !       CALL WRF_ERROR_FATAL ( message )
1629      endif
1630
1631       ! Vertically flip the chemistry data as it is given top down
1632       !     and heights in zfa are bottom up
1633       ! Fill 1D chemical profile array cprof
1634       ! Convert species 28-34 (lo-6:lo) from (molecules/cm3) to (mol/mol)
1635       irefcur = iref(nch-1)
1636       DO k = 1,kx
1637         zprof(k) = 0.5*(zfa_bdy(k)+zfa_bdy(k+1))
1638         if (irefcur .lt. lo-6) then
1639           cprof(k) = xl(irefcur,kx+1-k)
1640         else
1641           cprof(k) = xl(irefcur,kx+1-k)/dens(kx+1-k)
1642         end if
1643       ENDDO
1644
1645       ! Interpolate temp 3D chemical profile array to WRF field
1646       IF (z .LT. zprof(1)) THEN
1647         stor = cprof(1)
1648       ELSE IF (z .GT. zprof(kx)) THEN
1649         stor = cprof(kx)
1650       ELSE
1651         ! We can trap between two levels and linearly interpolate
1652         input_loop:  DO k = 1, kx-1
1653           IF (z .EQ. zprof(k) )THEN
1654             stor = cprof(k)
1655             EXIT input_loop
1656           ELSE IF ( (z .GT. zprof(k)) .AND. &
1657                     (z .LT. zprof(k+1)) ) THEN
1658             wgt0 = (z   - zprof(k+1)) / &
1659                    (zprof(k) - zprof(k+1))
1660             stor = MAX( wgt0 *cprof(k  ) + &
1661                      (1.-wgt0)*cprof(k+1), 0.)
1662             EXIT input_loop
1663           ENDIF
1664         ENDDO input_loop
1665       ENDIF
1666
1667       ! Here is where the chemistry value is constructed
1668       chem = fracref(nch-1)*stor*1.E6
1669
1670       ! special code for sulfate/h2so4
1671       if(nch.eq.p_sulf.and.p_nu0.gt.1)then
1672         chem=chem*(1.-so4vaptoaer)
1673       endif
1674
1675       RETURN
1676   END SUBROUTINE bdy_chem_value_racm
1677
1678 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1679
1680   SUBROUTINE bdy_chem_value ( chem, z, nch, numgas )
1681
1682     IMPLICIT NONE
1683
1684     REAL,    intent(OUT)  :: chem
1685     REAL,    intent(IN)   :: z          ! 3D height array
1686     INTEGER, intent(IN)   :: nch        ! index number of chemical species
1687     INTEGER, intent(IN)   :: numgas     ! index number of last gas species
1688
1689     INTEGER :: i, k, irefcur
1690
1691     REAL, DIMENSION(kx):: cprof         ! chemical profile, diff. index order
1692
1693     REAL, DIMENSION(1:kx):: zprof
1694     REAL                 :: stor
1695     REAL                 :: wgt0
1696
1697     CHARACTER (LEN=80) :: message
1698 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1699
1700 ! Check the number of species
1701 !     if((nch-1).gt.logg)return
1702 ! for radmkpp there is co2 and ch4 in the variable list
1703 !
1704      if(p_co2.gt.1)then
1705      if (nch.eq.p_co2)then
1706        chem=370.
1707        return
1708      endif
1709      if (nch.eq.p_ch4)then
1710        chem=1.7
1711        return
1712      endif
1713      endif
1714
1715 !     if( nch .GT. numgas) then
1716 !       message = ' Input_chem_profile: wrong number of chemical species'
1717 !       return
1718 !       CALL WRF_ERROR_FATAL ( message )
1719 !     endif
1720
1721       ! Vertically flip the chemistry data as it is given top down
1722       !     and heights in zfa are bottom up
1723       ! Fill 1D chemical profile array cprof
1724       ! Convert species 28-34 (lo-6:lo) from (molecules/cm3) to (mol/mol)
1725       irefcur = iref(nch-1)
1726       DO k = 1,kx
1727         zprof(k) = 0.5*(zfa_bdy(k)+zfa_bdy(k+1))
1728         if (irefcur .lt. lo-6) then
1729           cprof(k) = xl(irefcur,kx+1-k)
1730         else
1731           cprof(k) = xl(irefcur,kx+1-k)/dens(kx+1-k)
1732         end if
1733       ENDDO
1734
1735       ! Interpolate temp 3D chemical profile array to WRF field
1736       IF (z .LT. zprof(1)) THEN
1737         stor = cprof(1)
1738       ELSE IF (z .GT. zprof(kx)) THEN
1739         stor = cprof(kx)
1740       ELSE
1741         ! We can trap between two levels and linearly interpolate
1742         input_loop:  DO k = 1, kx-1
1743           IF (z .EQ. zprof(k) )THEN
1744             stor = cprof(k)
1745             EXIT input_loop
1746           ELSE IF (z .EQ. zprof(k+1) )THEN
1747             stor = cprof(k+1)
1748             EXIT input_loop
1749           ELSE IF ( (z .GT. zprof(k)) .AND. &
1750                     (z .LT. zprof(k+1)) ) THEN
1751             wgt0 = (z   - zprof(k+1)) / &
1752                    (zprof(k) - zprof(k+1))
1753             stor = MAX( wgt0 *cprof(k  ) + &
1754                      (1.-wgt0)*cprof(k+1), 0.)
1755             EXIT input_loop
1756           ENDIF
1757         ENDDO input_loop
1758       ENDIF
1759
1760       ! Here is where the chemistry value is constructed
1761       chem = fracref(nch-1)*stor*1.E6
1762
1763       ! special code for sulfate/h2so4
1764       if(nch.eq.p_sulf.and.p_nu0.gt.1)then
1765         chem=chem*(1.-so4vaptoaer)
1766       endif
1767
1768       RETURN
1769   END SUBROUTINE bdy_chem_value
1770 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1771
1772 #if (EM_CORE == 1 )
1773    SUBROUTINE flow_dep_bdy_chem  (  chem,                                       &
1774                                chem_bxs,chem_btxs,                                  &
1775                                chem_bxe,chem_btxe,                                  &
1776                                chem_bys,chem_btys,                                  &
1777                                chem_bye,chem_btye,                                  &
1778                                dt,                                              &
1779                                spec_bdy_width,z,                                &
1780                                have_bcs_chem,                        &
1781                                u, v, config_flags, alt, &
1782                                t,pb,p,t0,p1000mb,rcp,ph,phb,g, &
1783                                spec_zone, ic,           &
1784                                ids,ide, jds,jde, kds,kde,  & ! domain dims
1785                                ims,ime, jms,jme, kms,kme,  & ! memory dims
1786                                ips,ipe, jps,jpe, kps,kpe,  & ! patch  dims
1787                                its,ite, jts,jte, kts,kte )
1788
1789 !  This subroutine sets zero gradient conditions for outflow and a set profile value
1790 !  for inflow in the boundary specified region. Note that field must be unstaggered.
1791 !  The velocities, u and v, will only be used to check their sign (coupled vels OK)
1792 !  spec_zone is the width of the outer specified b.c.s that are set here.
1793 !  (JD August 2000)
1794
1795       IMPLICIT NONE
1796
1797       INTEGER,      INTENT(IN   )    :: ids,ide, jds,jde, kds,kde
1798       INTEGER,      INTENT(IN   )    :: ims,ime, jms,jme, kms,kme
1799       INTEGER,      INTENT(IN   )    :: ips,ipe, jps,jpe, kps,kpe
1800       INTEGER,      INTENT(IN   )    :: its,ite, jts,jte, kts,kte
1801       INTEGER,      INTENT(IN   )    :: spec_zone,spec_bdy_width,ic
1802       REAL,         INTENT(IN   )    :: dt
1803
1804
1805       REAL,  DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(INOUT) :: chem
1806       REAL,  DIMENSION( jms:jme , kds:kde , spec_bdy_width), INTENT(IN   ) :: chem_bxs, chem_bxe, chem_btxs, chem_btxe
1807       REAL,  DIMENSION( ims:ime , kds:kde , spec_bdy_width), INTENT(IN   ) :: chem_bys, chem_bye, chem_btys, chem_btye
1808       REAL,  DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(IN   ) :: z
1809       REAL,  DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(IN   ) :: alt
1810       REAL,  DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(IN   ) :: u
1811       REAL,  DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(IN   ) :: v
1812    REAL,  DIMENSION( ims:ime , kms:kme , jms:jme )         ,         &
1813           INTENT(IN   ) ::                                           &
1814                                ph,phb,t,pb,p
1815    real, INTENT (IN) :: g,rcp,t0,p1000mb
1816       TYPE( grid_config_rec_type ) config_flags
1817
1818       INTEGER    :: i, j, k, numgas
1819       INTEGER    :: ibs, ibe, jbs, jbe, itf, jtf, ktf
1820       INTEGER    :: i_inner, j_inner
1821       INTEGER    :: b_dist
1822       integer    :: itestbc, i_bdy_method
1823       real tempfac,convfac
1824       real       :: chem_bv_def
1825       logical    :: have_bcs_chem
1826
1827       chem_bv_def = conmin
1828       numgas = get_last_gas(config_flags%chem_opt)
1829       itestbc=0
1830       if(p_nu0.gt.1)itestbc=1
1831       ibs = ids
1832       ibe = ide-1
1833       itf = min(ite,ide-1)
1834       jbs = jds
1835       jbe = jde-1
1836       jtf = min(jte,jde-1)
1837       ktf = kde-1
1838
1839 ! i_bdy_method determines which "bdy_chem_value" routine to use
1840 !   1=radm2 or racm gas for  p_so2     <= ic <= p_ho2
1841 !   2=sorgam aerosol    for  p_so4aj   <= ic <= p_corn
1842 !   3=cbmz gas          for  p_hcl     <= ic <= p_isopo2
1843 !                        OR  p_dms     <= ic <= p_mtf
1844 !   4=mosaic aerosol    for  p_so4_a01 <= ic <= p_num_a01
1845 !                        OR  p_so4_a02 <= ic <= p_num_a02
1846 !                        OR  ...
1847 !   5=tracer mode
1848 !   0=none              for all other ic values
1849 ! (note:  some cbmz packages use dms,...,mtf while others do not)
1850 ! (note:  different mosaic packages use different number of sections)
1851       i_bdy_method = 0
1852       if ((ic .ge. p_so2) .and. (ic .le. p_ho2)) then
1853           i_bdy_method = 1
1854
1855         if (config_flags%chem_opt == RACM_KPP .or.          &
1856             config_flags%chem_opt == GOCARTRACM_KPP .or.      &
1857             config_flags%chem_opt == RACMSORG_KPP .or.      &
1858             config_flags%chem_opt == RACM_MIM_KPP .or.      &
1859             config_flags%chem_opt == RACMSORG_KPP ) then
1860           i_bdy_method = 9
1861         end if
1862         if (config_flags%chem_opt == RACMPM_KPP ) then
1863           i_bdy_method = 9
1864         end if
1865
1866
1867       else if ((ic .ge. p_so4aj) .and. (ic .le. p_corn)) then
1868           i_bdy_method = 2
1869       else if ((ic .ge. p_hcl) .and. (ic .le. p_isopo2)) then
1870           i_bdy_method = 3
1871       else if ((ic .ge. p_dms) .and. (ic .le. p_mtf)) then
1872           i_bdy_method = 3
1873       else if ((ic .ge. p_so4_a01) .and. (ic .le. p_num_a01)) then
1874           i_bdy_method = 4
1875       else if ((ic .ge. p_so4_a02) .and. (ic .le. p_num_a02)) then
1876           i_bdy_method = 4
1877       else if ((ic .ge. p_so4_a03) .and. (ic .le. p_num_a03)) then
1878           i_bdy_method = 4
1879       else if ((ic .ge. p_so4_a04) .and. (ic .le. p_num_a04)) then
1880           i_bdy_method = 4
1881       else if ((ic .ge. p_so4_a05) .and. (ic .le. p_num_a05)) then
1882           i_bdy_method = 4
1883       else if ((ic .ge. p_so4_a06) .and. (ic .le. p_num_a06)) then
1884           i_bdy_method = 4
1885       else if ((ic .ge. p_so4_a07) .and. (ic .le. p_num_a07)) then
1886           i_bdy_method = 4
1887       else if ((ic .ge. p_so4_a08) .and. (ic .le. p_num_a08)) then
1888           i_bdy_method = 4
1889       else if (config_flags%chem_opt == CHEM_TRACER) then
1890           i_bdy_method = 5
1891       else if (config_flags%chem_opt == CHEM_TRACE2) then
1892           i_bdy_method = 5
1893       else if (config_flags%chem_opt == GOCART_SIMPLE) then
1894           i_bdy_method = 7
1895       end if
1896       if (have_bcs_chem) i_bdy_method =6
1897       if (ic .lt. param_first_scalar) i_bdy_method = 0
1898
1899 !----------------------------------------------------------------------
1900 !      if (i_bdy_method .eq. 1) then
1901 !          print 90010, '_bdy_radm2  for ic=', ic, i_bdy_method
1902 !      else if (i_bdy_method .eq. 2) then
1903 !          print 90010, '_bdy_sorgam for ic=', ic, i_bdy_method
1904 !      else if (i_bdy_method .eq. 3) then
1905 !          print 90010, '_bdy_cbmz   for ic=', ic, i_bdy_method
1906 !      else if (i_bdy_method .eq. 4) then
1907 !          print 90010, '_bdy_mosaic for ic=', ic, i_bdy_method
1908 !      else if (i_bdy_method .eq. 5) then
1909 !          print 90010, '_bdy_tracer for ic=', ic, i_bdy_method
1910 !          print *,numgas,num_chem
1911 !      else
1912 !          print 90010, '_bdy_NONE** for ic=', ic, i_bdy_method
1913 !      end if
1914 !90010 format( a, 2(1x,i5) )
1915 !90020 format( a, 1p, 2e12.2 )
1916 !----------------------------------------------------------------------
1917
1918       IF (jts - jbs .lt. spec_zone) THEN
1919 ! Y-start boundary
1920         DO j = jts, min(jtf,jbs+spec_zone-1)
1921           b_dist = j - jbs
1922           DO k = kts, ktf
1923             DO i = max(its,b_dist+ibs), min(itf,ibe-b_dist)
1924               i_inner = max(i,ibs+spec_zone)
1925               i_inner = min(i_inner,ibe-spec_zone)
1926               IF(v(i,k,j) .lt. 0.)THEN
1927                 chem(i,k,j) = chem(i_inner,k,jbs+spec_zone)
1928               ELSE
1929                 if (i_bdy_method .eq. 1) then
1930                    CALL bdy_chem_value (   &
1931                         chem(i,k,j), z(i,k,j), ic, numgas )
1932                 else if (i_bdy_method .eq. 9) then
1933                    CALL bdy_chem_value_racm(   &
1934                         chem(i,k,j), z(i,k,j), ic, numgas,p_co2 )
1935                 else if (i_bdy_method .eq. 2) then
1936                    tempfac=(t(i,k,j)+t0)*((p(i,k,j) + pb(i,k,j))/p1000mb)**rcp
1937                    convfac=(p(i,k,j)+pb(i,k,j))/rgasuniv/tempfac
1938                    CALL bdy_chem_value_sorgam (   &
1939                         chem(i,k,j), z(i,k,j), ic, config_flags,   &
1940                         alt(i,k,j),convfac,g)
1941                 else if (i_bdy_method .eq. 3) then
1942                    CALL bdy_chem_value_cbmz (   &
1943                         chem(i,k,j), z(i,k,j), ic, numgas )
1944                 else if (i_bdy_method .eq. 4) then
1945                    CALL bdy_chem_value_mosaic (   &
1946                         chem(i,k,j), alt(i,k,j), z(i,k,j), ic, config_flags )
1947                 else if (i_bdy_method .eq. 5) then
1948                    CALL bdy_chem_value_tracer ( chem(i,k,j), ic )
1949                 else if (i_bdy_method .eq. 7) then
1950                    CALL bdy_chem_value_gocart ( chem(i,k,j), ic )
1951                 else if (i_bdy_method .eq. 6) then
1952                    CALL bdy_chem_value_gcm ( chem(i,k,j),chem_bys(i,k,1),chem_btys(i,k,1),dt,ic)
1953                 else
1954                    chem(i,k,j) = chem_bv_def
1955                 endif
1956               ENDIF
1957             ENDDO
1958           ENDDO
1959         ENDDO
1960       ENDIF
1961       IF (jbe - jtf .lt. spec_zone) THEN
1962 ! Y-end boundary
1963         DO j = max(jts,jbe-spec_zone+1), jtf
1964           b_dist = jbe - j
1965           DO k = kts, ktf
1966             DO i = max(its,b_dist+ibs), min(itf,ibe-b_dist)
1967               i_inner = max(i,ibs+spec_zone)
1968               i_inner = min(i_inner,ibe-spec_zone)
1969               IF(v(i,k,j+1) .gt. 0.)THEN
1970                 chem(i,k,j) = chem(i_inner,k,jbe-spec_zone)
1971               ELSE
1972                 if (i_bdy_method .eq. 1) then
1973                    CALL bdy_chem_value (   &
1974                         chem(i,k,j), z(i,k,j), ic, numgas )
1975                 else if (i_bdy_method .eq. 9) then
1976                    CALL bdy_chem_value_racm (   &
1977                         chem(i,k,j), z(i,k,j), ic, numgas,p_co2 )
1978                 else if (i_bdy_method .eq. 2) then
1979                    tempfac=(t(i,k,j)+t0)*((p(i,k,j) + pb(i,k,j))/p1000mb)**rcp
1980                    convfac=(p(i,k,j)+pb(i,k,j))/rgasuniv/tempfac
1981                    CALL bdy_chem_value_sorgam (   &
1982                         chem(i,k,j), z(i,k,j), ic, config_flags,   &
1983                         alt(i,k,j),convfac,g)
1984                 else if (i_bdy_method .eq. 3) then
1985                    CALL bdy_chem_value_cbmz (   &
1986                         chem(i,k,j), z(i,k,j), ic, numgas )
1987                 else if (i_bdy_method .eq. 4) then
1988                    CALL bdy_chem_value_mosaic (   &
1989                         chem(i,k,j), alt(i,k,j), z(i,k,j), ic, config_flags )
1990                 else if (i_bdy_method .eq. 5) then
1991                    CALL bdy_chem_value_tracer ( chem(i,k,j), ic )
1992                 else if (i_bdy_method .eq. 6) then
1993                    CALL bdy_chem_value_gcm ( chem(i,k,j),chem_bye(i,k,1),chem_btye(i,k,1),dt,ic)
1994                 else if (i_bdy_method .eq. 7) then
1995                    CALL bdy_chem_value_gocart ( chem(i,k,j), ic )
1996                 else
1997                    chem(i,k,j) = chem_bv_def
1998                 endif
1999               ENDIF
2000             ENDDO
2001           ENDDO
2002         ENDDO
2003       ENDIF
2004
2005       IF (its - ibs .lt. spec_zone) THEN
2006 ! X-start boundary
2007         DO i = its, min(itf,ibs+spec_zone-1)
2008           b_dist = i - ibs
2009           DO k = kts, ktf
2010             DO j = max(jts,b_dist+jbs+1), min(jtf,jbe-b_dist-1)
2011               j_inner = max(j,jbs+spec_zone)
2012               j_inner = min(j_inner,jbe-spec_zone)
2013               IF(u(i,k,j) .lt. 0.)THEN
2014                 chem(i,k,j) = chem(ibs+spec_zone,k,j_inner)
2015               ELSE
2016                 if (i_bdy_method .eq. 1) then
2017                    CALL bdy_chem_value (   &
2018                         chem(i,k,j), z(i,k,j), ic, numgas )
2019                 else if (i_bdy_method .eq. 9) then
2020                    CALL bdy_chem_value_racm (   &
2021                         chem(i,k,j), z(i,k,j), ic, numgas,p_co2 )
2022                 else if (i_bdy_method .eq. 2) then
2023                    tempfac=(t(i,k,j)+t0)*((p(i,k,j) + pb(i,k,j))/p1000mb)**rcp
2024                    convfac=(p(i,k,j)+pb(i,k,j))/rgasuniv/tempfac
2025                    CALL bdy_chem_value_sorgam (   &
2026                         chem(i,k,j), z(i,k,j), ic, config_flags,   &
2027                         alt(i,k,j),convfac,g)
2028                 else if (i_bdy_method .eq. 3) then
2029                    CALL bdy_chem_value_cbmz (   &
2030                         chem(i,k,j), z(i,k,j), ic, numgas )
2031                 else if (i_bdy_method .eq. 4) then
2032                    CALL bdy_chem_value_mosaic (   &
2033                         chem(i,k,j), alt(i,k,j), z(i,k,j), ic, config_flags )
2034                 else if (i_bdy_method .eq. 5) then
2035                    CALL bdy_chem_value_tracer ( chem(i,k,j), ic )
2036                 else if (i_bdy_method .eq. 6) then
2037                    CALL bdy_chem_value_gcm ( chem(i,k,j),chem_bxs(j,k,1),chem_btxs(j,k,1),dt,ic)
2038                 else if (i_bdy_method .eq. 7) then
2039                    CALL bdy_chem_value_gocart ( chem(i,k,j), ic )
2040                 else
2041                    chem(i,k,j) = chem_bv_def
2042                 endif
2043               ENDIF
2044             ENDDO
2045           ENDDO
2046         ENDDO
2047       ENDIF
2048
2049       IF (ibe - itf .lt. spec_zone) THEN
2050 ! X-end boundary
2051         DO i = max(its,ibe-spec_zone+1), itf
2052           b_dist = ibe - i
2053           DO k = kts, ktf
2054             DO j = max(jts,b_dist+jbs+1), min(jtf,jbe-b_dist-1)
2055               j_inner = max(j,jbs+spec_zone)
2056               j_inner = min(j_inner,jbe-spec_zone)
2057               IF(u(i+1,k,j) .gt. 0.)THEN
2058                 chem(i,k,j) = chem(ibe-spec_zone,k,j_inner)
2059               ELSE
2060                 if (i_bdy_method .eq. 1) then
2061                    CALL bdy_chem_value (   &
2062                         chem(i,k,j), z(i,k,j), ic, numgas )
2063                 else if (i_bdy_method .eq. 9) then
2064                    CALL bdy_chem_value_racm (   &
2065                         chem(i,k,j), z(i,k,j), ic, numgas,p_co2 )
2066                 else if (i_bdy_method .eq. 2) then
2067                    tempfac=(t(i,k,j)+t0)*((p(i,k,j) + pb(i,k,j))/p1000mb)**rcp
2068                    convfac=(p(i,k,j)+pb(i,k,j))/rgasuniv/tempfac
2069                    CALL bdy_chem_value_sorgam (   &
2070                         chem(i,k,j), z(i,k,j), ic, config_flags,   &
2071                         alt(i,k,j),convfac,g)
2072                 else if (i_bdy_method .eq. 3) then
2073                    CALL bdy_chem_value_cbmz (   &
2074                         chem(i,k,j), z(i,k,j), ic, numgas )
2075                 else if (i_bdy_method .eq. 4) then
2076                    CALL bdy_chem_value_mosaic (   &
2077                         chem(i,k,j), alt(i,k,j), z(i,k,j), ic, config_flags )
2078                 else if (i_bdy_method .eq. 5) then
2079                    CALL bdy_chem_value_tracer ( chem(i,k,j), ic )
2080                 else if (i_bdy_method .eq. 6) then
2081                    CALL bdy_chem_value_gcm ( chem(i,k,j),chem_bxe(j,k,1),chem_btxe(j,k,1),dt,ic)
2082                 else if (i_bdy_method .eq. 7) then
2083                    CALL bdy_chem_value_gocart ( chem(i,k,j), ic )
2084                 else
2085                    chem(i,k,j) = chem_bv_def
2086                 endif
2087               ENDIF
2088             ENDDO
2089           ENDDO
2090         ENDDO
2091       ENDIF
2092
2093    END SUBROUTINE flow_dep_bdy_chem
2094 #else
2095    SUBROUTINE flow_dep_bdy_chem  (  chem, chem_b,chem_bt,dt,                    &
2096                                spec_bdy_width,z,                                &
2097                                ijds, ijde,have_bcs_chem,                        &
2098                                u, v, config_flags, alt, &
2099                                t,pb,p,t0,p1000mb,rcp,ph,phb,g, &
2100                                spec_zone, ic,           &
2101                                ids,ide, jds,jde, kds,kde,  & ! domain dims
2102                                ims,ime, jms,jme, kms,kme,  & ! memory dims
2103                                ips,ipe, jps,jpe, kps,kpe,  & ! patch  dims
2104                                its,ite, jts,jte, kts,kte )
2105
2106 !  This subroutine sets zero gradient conditions for outflow and a set profile value
2107 !  for inflow in the boundary specified region. Note that field must be unstaggered.
2108 !  The velocities, u and v, will only be used to check their sign (coupled vels OK)
2109 !  spec_zone is the width of the outer specified b.c.s that are set here.
2110 !  (JD August 2000)
2111
2112       IMPLICIT NONE
2113
2114       INTEGER,      INTENT(IN   )    :: ids,ide, jds,jde, kds,kde
2115       INTEGER,      INTENT(IN   )    :: ims,ime, jms,jme, kms,kme
2116       INTEGER,      INTENT(IN   )    :: ips,ipe, jps,jpe, kps,kpe
2117       INTEGER,      INTENT(IN   )    :: its,ite, jts,jte, kts,kte
2118       INTEGER,      INTENT(IN   )    :: ijds,ijde
2119       INTEGER,      INTENT(IN   )    :: spec_zone,spec_bdy_width,ic
2120       REAL,         INTENT(IN   )    :: dt
2121
2122
2123       REAL,  DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(INOUT) :: chem
2124       REAL,  DIMENSION( ijds:ijde , kds:kde , spec_bdy_width, 4 ), INTENT(IN   ) :: chem_b
2125       REAL,  DIMENSION( ijds:ijde , kds:kde , spec_bdy_width, 4 ), INTENT(IN   ) :: chem_bt
2126       REAL,  DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(IN   ) :: z
2127       REAL,  DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(IN   ) :: alt
2128       REAL,  DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(IN   ) :: u
2129       REAL,  DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(IN   ) :: v
2130    REAL,  DIMENSION( ims:ime , kms:kme , jms:jme )         ,         &
2131           INTENT(IN   ) ::                                           &
2132                                ph,phb,t,pb,p
2133    real, INTENT (IN) :: g,rcp,t0,p1000mb
2134       TYPE( grid_config_rec_type ) config_flags
2135
2136       INTEGER    :: i, j, k, numgas
2137       INTEGER    :: ibs, ibe, jbs, jbe, itf, jtf, ktf
2138       INTEGER    :: i_inner, j_inner
2139       INTEGER    :: b_dist
2140       integer    :: itestbc, i_bdy_method
2141       real tempfac,convfac
2142       real       :: chem_bv_def
2143       logical    :: have_bcs_chem
2144
2145       chem_bv_def = conmin
2146       numgas = get_last_gas(config_flags%chem_opt)
2147       itestbc=0
2148       if(p_nu0.gt.1)itestbc=1
2149       ibs = ids
2150       ibe = ide-1
2151       itf = min(ite,ide-1)
2152       jbs = jds
2153       jbe = jde-1
2154       jtf = min(jte,jde-1)
2155       ktf = kde-1
2156
2157 ! i_bdy_method determines which "bdy_chem_value" routine to use
2158 !   1=radm2 or racm gas for  p_so2     <= ic <= p_ho2
2159 !   2=sorgam aerosol    for  p_so4aj   <= ic <= p_corn
2160 !   3=cbmz gas          for  p_hcl     <= ic <= p_isopo2
2161 !                        OR  p_dms     <= ic <= p_mtf
2162 !   4=mosaic aerosol    for  p_so4_a01 <= ic <= p_num_a01
2163 !                        OR  p_so4_a02 <= ic <= p_num_a02
2164 !                        OR  ...
2165 !   5=tracer mode
2166 !   0=none              for all other ic values
2167 ! (note:  some cbmz packages use dms,...,mtf while others do not)
2168 ! (note:  different mosaic packages use different number of sections)
2169       i_bdy_method = 0
2170       if ((ic .ge. p_so2) .and. (ic .le. p_ho2)) then
2171           i_bdy_method = 1
2172
2173         if (config_flags%chem_opt == RACM_KPP .or.          &
2174             config_flags%chem_opt == GOCARTRACM_KPP .or.      &
2175             config_flags%chem_opt == RACMSORG_KPP .or.      &
2176             config_flags%chem_opt == RACM_MIM_KPP .or.      &
2177             config_flags%chem_opt == RACMSORG_KPP ) then
2178           i_bdy_method = 9
2179         end if
2180         if (config_flags%chem_opt == RACMPM_KPP ) then
2181           i_bdy_method = 9
2182         end if
2183
2184
2185       else if ((ic .ge. p_so4aj) .and. (ic .le. p_corn)) then
2186           i_bdy_method = 2
2187       else if ((ic .ge. p_hcl) .and. (ic .le. p_isopo2)) then
2188           i_bdy_method = 3
2189       else if ((ic .ge. p_dms) .and. (ic .le. p_mtf)) then
2190           i_bdy_method = 3
2191       else if ((ic .ge. p_so4_a01) .and. (ic .le. p_num_a01)) then
2192           i_bdy_method = 4
2193       else if ((ic .ge. p_so4_a02) .and. (ic .le. p_num_a02)) then
2194           i_bdy_method = 4
2195       else if ((ic .ge. p_so4_a03) .and. (ic .le. p_num_a03)) then
2196           i_bdy_method = 4
2197       else if ((ic .ge. p_so4_a04) .and. (ic .le. p_num_a04)) then
2198           i_bdy_method = 4
2199       else if ((ic .ge. p_so4_a05) .and. (ic .le. p_num_a05)) then
2200           i_bdy_method = 4
2201       else if ((ic .ge. p_so4_a06) .and. (ic .le. p_num_a06)) then
2202           i_bdy_method = 4
2203       else if ((ic .ge. p_so4_a07) .and. (ic .le. p_num_a07)) then
2204           i_bdy_method = 4
2205       else if ((ic .ge. p_so4_a08) .and. (ic .le. p_num_a08)) then
2206           i_bdy_method = 4
2207       else if (config_flags%chem_opt == CHEM_TRACER) then
2208           i_bdy_method = 5
2209       else if (config_flags%chem_opt == CHEM_TRACE2) then
2210           i_bdy_method = 5
2211       end if
2212       if (have_bcs_chem) i_bdy_method =6
2213       if (ic .lt. param_first_scalar) i_bdy_method = 0
2214
2215 !----------------------------------------------------------------------
2216 !      if (i_bdy_method .eq. 1) then
2217 !          print 90010, '_bdy_radm2  for ic=', ic, i_bdy_method
2218 !      else if (i_bdy_method .eq. 2) then
2219 !          print 90010, '_bdy_sorgam for ic=', ic, i_bdy_method
2220 !      else if (i_bdy_method .eq. 3) then
2221 !          print 90010, '_bdy_cbmz   for ic=', ic, i_bdy_method
2222 !      else if (i_bdy_method .eq. 4) then
2223 !          print 90010, '_bdy_mosaic for ic=', ic, i_bdy_method
2224 !      else if (i_bdy_method .eq. 5) then
2225 !          print 90010, '_bdy_tracer for ic=', ic, i_bdy_method
2226 !      else
2227 !          print 90010, '_bdy_NONE** for ic=', ic, i_bdy_method
2228 !      end if
2229 !90010 format( a, 2(1x,i5) )
2230 !90020 format( a, 1p, 2e12.2 )
2231 !----------------------------------------------------------------------
2232
2233 !     if(ic.eq.p_O3)THEN
2234 !     write(0,*)'in flow_chem ',jts,jbs,spec_zone
2235 !     write(0,*)'in flow_chem ',its,ibs,b_dist,i_bdy_method,ic
2236 !     endif
2237       IF (jts - jbs .lt. spec_zone) THEN
2238 ! Y-start boundary
2239         DO j = jts, min(jtf,jbs+spec_zone-1)
2240           b_dist = j - jbs
2241           DO k = kts, ktf
2242             DO i = max(its,b_dist+ibs), min(itf,ibe-b_dist)
2243               i_inner = max(i,ibs+spec_zone)
2244               i_inner = min(i_inner,ibe-spec_zone)
2245               IF(v(i,k,j) .lt. 0.)THEN
2246                 chem(i,k,j) = chem(i_inner,k,jbs+spec_zone)
2247 !               if(j.eq.jts+1.and.k.eq.kts.and.ic.eq.p_o3)then
2248 !                  write(0,*)'Yflow',i,j,k,i_bdy_method
2249 !                  write(0,*)chem(i_inner,k,jbs+spec_zone),v(i,k,j)
2250 !               endif
2251               ELSE
2252                 if (i_bdy_method .eq. 1) then
2253                    CALL bdy_chem_value (   &
2254                         chem(i,k,j), z(i,k,j), ic, numgas )
2255                 else if (i_bdy_method .eq. 9) then
2256                    CALL bdy_chem_value_racm(   &
2257                         chem(i,k,j), z(i,k,j), ic, numgas,p_co2 )
2258                 else if (i_bdy_method .eq. 2) then
2259                    tempfac=(t(i,k,j)+t0)*((p(i,k,j) + pb(i,k,j))/p1000mb)**rcp
2260                    convfac=(p(i,k,j)+pb(i,k,j))/rgasuniv/tempfac
2261                    CALL bdy_chem_value_sorgam (   &
2262                         chem(i,k,j), z(i,k,j), ic, config_flags,   &
2263                         alt(i,k,j),convfac,g)
2264                 else if (i_bdy_method .eq. 3) then
2265                    CALL bdy_chem_value_cbmz (   &
2266                         chem(i,k,j), z(i,k,j), ic, numgas )
2267                 else if (i_bdy_method .eq. 4) then
2268                    CALL bdy_chem_value_mosaic (   &
2269                         chem(i,k,j), alt(i,k,j), z(i,k,j), ic, config_flags )
2270                 else if (i_bdy_method .eq. 5) then
2271                    CALL bdy_chem_value_tracer ( chem(i,k,j), ic )
2272                 else if (i_bdy_method .eq. 6) then
2273                    CALL bdy_chem_value_gcm ( chem(i,k,j),chem_b(i,k,1,P_YSB),chem_bt(i,k,1,P_YSB),dt,ic)
2274 !               if(k.eq.kts.and.ic.eq.p_o3)then
2275 !                  write(0,*)'Ygcm',i,j,k,i_bdy_method
2276 !                  write(0,*)chem(i,k,j),chem_b(i,k,1,P_YSB),chem_bt(i,k,1,P_YSB),dt
2277 !               endif
2278                 else
2279                    chem(i,k,j) = chem_bv_def
2280                 endif
2281               ENDIF
2282             ENDDO
2283           ENDDO
2284         ENDDO
2285       ENDIF
2286       IF (jbe - jtf .lt. spec_zone) THEN
2287 ! Y-end boundary
2288         DO j = max(jts,jbe-spec_zone+1), jtf
2289           b_dist = jbe - j
2290           DO k = kts, ktf
2291             DO i = max(its,b_dist+ibs), min(itf,ibe-b_dist)
2292               i_inner = max(i,ibs+spec_zone)
2293               i_inner = min(i_inner,ibe-spec_zone)
2294               IF(v(i,k,j+1) .gt. 0.)THEN
2295                 chem(i,k,j) = chem(i_inner,k,jbe-spec_zone)
2296               ELSE
2297                 if (i_bdy_method .eq. 1) then
2298                    CALL bdy_chem_value (   &
2299                         chem(i,k,j), z(i,k,j), ic, numgas )
2300                 else if (i_bdy_method .eq. 9) then
2301                    CALL bdy_chem_value_racm (   &
2302                         chem(i,k,j), z(i,k,j), ic, numgas,p_co2 )
2303                 else if (i_bdy_method .eq. 2) then
2304                    tempfac=(t(i,k,j)+t0)*((p(i,k,j) + pb(i,k,j))/p1000mb)**rcp
2305                    convfac=(p(i,k,j)+pb(i,k,j))/rgasuniv/tempfac
2306                    CALL bdy_chem_value_sorgam (   &
2307                         chem(i,k,j), z(i,k,j), ic, config_flags,   &
2308                         alt(i,k,j),convfac,g)
2309                 else if (i_bdy_method .eq. 3) then
2310                    CALL bdy_chem_value_cbmz (   &
2311                         chem(i,k,j), z(i,k,j), ic, numgas )
2312                 else if (i_bdy_method .eq. 4) then
2313                    CALL bdy_chem_value_mosaic (   &
2314                         chem(i,k,j), alt(i,k,j), z(i,k,j), ic, config_flags )
2315                 else if (i_bdy_method .eq. 5) then
2316                    CALL bdy_chem_value_tracer ( chem(i,k,j), ic )
2317                 else if (i_bdy_method .eq. 6) then
2318                    CALL bdy_chem_value_gcm ( chem(i,k,j),chem_b(i,k,1,P_YEB),chem_bt(i,k,1,P_YEB),dt,ic)
2319                 else
2320                    chem(i,k,j) = chem_bv_def
2321                 endif
2322               ENDIF
2323             ENDDO
2324           ENDDO
2325         ENDDO
2326       ENDIF
2327
2328       IF (its - ibs .lt. spec_zone) THEN
2329 ! X-start boundary
2330         DO i = its, min(itf,ibs+spec_zone-1)
2331           b_dist = i - ibs
2332           DO k = kts, ktf
2333             DO j = max(jts,b_dist+jbs+1), min(jtf,jbe-b_dist-1)
2334               j_inner = max(j,jbs+spec_zone)
2335               j_inner = min(j_inner,jbe-spec_zone)
2336               IF(u(i,k,j) .lt. 0.)THEN
2337                 chem(i,k,j) = chem(ibs+spec_zone,k,j_inner)
2338               ELSE
2339                 if (i_bdy_method .eq. 1) then
2340                    CALL bdy_chem_value (   &
2341                         chem(i,k,j), z(i,k,j), ic, numgas )
2342                 else if (i_bdy_method .eq. 9) then
2343                    CALL bdy_chem_value_racm (   &
2344                         chem(i,k,j), z(i,k,j), ic, numgas,p_co2 )
2345                 else if (i_bdy_method .eq. 2) then
2346                    tempfac=(t(i,k,j)+t0)*((p(i,k,j) + pb(i,k,j))/p1000mb)**rcp
2347                    convfac=(p(i,k,j)+pb(i,k,j))/rgasuniv/tempfac
2348                    CALL bdy_chem_value_sorgam (   &
2349                         chem(i,k,j), z(i,k,j), ic, config_flags,   &
2350                         alt(i,k,j),convfac,g)
2351                 else if (i_bdy_method .eq. 3) then
2352                    CALL bdy_chem_value_cbmz (   &
2353                         chem(i,k,j), z(i,k,j), ic, numgas )
2354                 else if (i_bdy_method .eq. 4) then
2355                    CALL bdy_chem_value_mosaic (   &
2356                         chem(i,k,j), alt(i,k,j), z(i,k,j), ic, config_flags )
2357                 else if (i_bdy_method .eq. 5) then
2358                    CALL bdy_chem_value_tracer ( chem(i,k,j), ic )
2359                 else if (i_bdy_method .eq. 6) then
2360                    CALL bdy_chem_value_gcm ( chem(i,k,j),chem_b(j,k,1,P_XSB),chem_bt(j,k,1,P_XSB),dt,ic)
2361                 else
2362                    chem(i,k,j) = chem_bv_def
2363                 endif
2364               ENDIF
2365             ENDDO
2366           ENDDO
2367         ENDDO
2368       ENDIF
2369
2370       IF (ibe - itf .lt. spec_zone) THEN
2371 ! X-end boundary
2372         DO i = max(its,ibe-spec_zone+1), itf
2373           b_dist = ibe - i
2374           DO k = kts, ktf
2375             DO j = max(jts,b_dist+jbs+1), min(jtf,jbe-b_dist-1)
2376               j_inner = max(j,jbs+spec_zone)
2377               j_inner = min(j_inner,jbe-spec_zone)
2378               IF(u(i+1,k,j) .gt. 0.)THEN
2379                 chem(i,k,j) = chem(ibe-spec_zone,k,j_inner)
2380               ELSE
2381                 if (i_bdy_method .eq. 1) then
2382                    CALL bdy_chem_value (   &
2383                         chem(i,k,j), z(i,k,j), ic, numgas )
2384                 else if (i_bdy_method .eq. 9) then
2385                    CALL bdy_chem_value_racm (   &
2386                         chem(i,k,j), z(i,k,j), ic, numgas,p_co2 )
2387                 else if (i_bdy_method .eq. 2) then
2388                    tempfac=(t(i,k,j)+t0)*((p(i,k,j) + pb(i,k,j))/p1000mb)**rcp
2389                    convfac=(p(i,k,j)+pb(i,k,j))/rgasuniv/tempfac
2390                    CALL bdy_chem_value_sorgam (   &
2391                         chem(i,k,j), z(i,k,j), ic, config_flags,   &
2392                         alt(i,k,j),convfac,g)
2393                 else if (i_bdy_method .eq. 3) then
2394                    CALL bdy_chem_value_cbmz (   &
2395                         chem(i,k,j), z(i,k,j), ic, numgas )
2396                 else if (i_bdy_method .eq. 4) then
2397                    CALL bdy_chem_value_mosaic (   &
2398                         chem(i,k,j), alt(i,k,j), z(i,k,j), ic, config_flags )
2399                 else if (i_bdy_method .eq. 5) then
2400                    CALL bdy_chem_value_tracer ( chem(i,k,j), ic )
2401                 else if (i_bdy_method .eq. 6) then
2402                    CALL bdy_chem_value_gcm ( chem(i,k,j),chem_b(j,k,1,P_XEB),chem_bt(j,k,1,P_XEB),dt,ic)
2403                 else
2404                    chem(i,k,j) = chem_bv_def
2405                 endif
2406               ENDIF
2407             ENDDO
2408           ENDDO
2409         ENDDO
2410       ENDIF
2411
2412    END SUBROUTINE flow_dep_bdy_chem
2413 #endif
2414 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2415 !
2416 ! this is a kludge routine as of now....
2417 !
2418
2419
2420    SUBROUTINE flow_dep_bdy_s1 (  field,                     &
2421                                u, v, config_flags, &
2422                                spec_zone,                  &
2423                                ids,ide, jds,jde, kds,kde,  & ! domain dims
2424                                ims,ime, jms,jme, kms,kme,  & ! memory dims
2425                                ips,ipe, jps,jpe, kps,kpe,  & ! patch  dims
2426                                its,ite, jts,jte, kts,kte)
2427
2428 !  This subroutine sets zero gradient conditions for outflow and zero value
2429 !  for inflow in the boundary specified region. Note that field must be unstaggered.
2430 !  The velocities, u and v, will only be used to check their sign (coupled vels OK)
2431 !  spec_zone is the width of the outer specified b.c.s that are set here.
2432 !  (JD August 2000)
2433
2434       IMPLICIT NONE
2435
2436       INTEGER,      INTENT(IN   )    :: ids,ide, jds,jde, kds,kde
2437       INTEGER,      INTENT(IN   )    :: ims,ime, jms,jme, kms,kme
2438       INTEGER,      INTENT(IN   )    :: ips,ipe, jps,jpe, kps,kpe
2439       INTEGER,      INTENT(IN   )    :: its,ite, jts,jte, kts,kte
2440       INTEGER,      INTENT(IN   )    :: spec_zone
2441
2442
2443       REAL,  DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(INOUT) :: field
2444       REAL,  DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(IN   ) :: u
2445       REAL,  DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(IN   ) :: v
2446       TYPE( grid_config_rec_type ) config_flags
2447
2448       INTEGER    :: i, j, k, ibs, ibe, jbs, jbe, itf, jtf, ktf, i_inner, j_inner
2449       INTEGER    :: b_dist, b_limit
2450       LOGICAL    :: periodic_x
2451 !----------------------------------------------
2452 !jeff
2453 !--   hardwire outmost lateral boundary value at constant value
2454
2455       real value_bc
2456
2457       value_bc = 1.0
2458
2459 !-----------------------------------------------
2460
2461       periodic_x = config_flags%periodic_x
2462
2463       ibs = ids
2464       ibe = ide-1
2465       itf = min(ite,ide-1)
2466       jbs = jds
2467       jbe = jde-1
2468       jtf = min(jte,jde-1)
2469       ktf = kde-1
2470
2471       IF (jts - jbs .lt. spec_zone) THEN
2472 ! Y-start boundary
2473         DO j = jts, min(jtf,jbs+spec_zone-1)
2474           b_dist = j - jbs
2475           b_limit = b_dist
2476           IF(periodic_x)b_limit = 0
2477           DO k = kts, ktf
2478             DO i = max(its,b_limit+ibs), min(itf,ibe-b_limit)
2479               i_inner = max(i,ibs+spec_zone)
2480               i_inner = min(i_inner,ibe-spec_zone)
2481               IF(periodic_x)i_inner = i
2482               IF(v(i,k,j) .lt. 0.)THEN
2483                 field(i,k,j) = field(i_inner,k,jbs+spec_zone)
2484               ELSE
2485                 field(i,k,j) = value_bc
2486               ENDIF
2487             ENDDO
2488           ENDDO
2489         ENDDO
2490       ENDIF
2491       IF (jbe - jtf .lt. spec_zone) THEN
2492 ! Y-end boundary
2493         DO j = max(jts,jbe-spec_zone+1), jtf
2494           b_dist = jbe - j
2495           b_limit = b_dist
2496           IF(periodic_x)b_limit = 0
2497           DO k = kts, ktf
2498             DO i = max(its,b_limit+ibs), min(itf,ibe-b_limit)
2499               i_inner = max(i,ibs+spec_zone)
2500               i_inner = min(i_inner,ibe-spec_zone)
2501               IF(periodic_x)i_inner = i
2502               IF(v(i,k,j+1) .gt. 0.)THEN
2503                 field(i,k,j) = field(i_inner,k,jbe-spec_zone)
2504               ELSE
2505                 field(i,k,j) = value_bc
2506               ENDIF
2507             ENDDO
2508           ENDDO
2509         ENDDO
2510       ENDIF
2511
2512     IF(.NOT.periodic_x)THEN
2513       IF (its - ibs .lt. spec_zone) THEN
2514 ! X-start boundary
2515         DO i = its, min(itf,ibs+spec_zone-1)
2516           b_dist = i - ibs
2517           DO k = kts, ktf
2518             DO j = max(jts,b_dist+jbs+1), min(jtf,jbe-b_dist-1)
2519               j_inner = max(j,jbs+spec_zone)
2520               j_inner = min(j_inner,jbe-spec_zone)
2521               IF(u(i,k,j) .lt. 0.)THEN
2522                 field(i,k,j) = field(ibs+spec_zone,k,j_inner)
2523               ELSE
2524                 field(i,k,j) = value_bc
2525               ENDIF
2526             ENDDO
2527           ENDDO
2528         ENDDO
2529       ENDIF
2530
2531       IF (ibe - itf .lt. spec_zone) THEN
2532 ! X-end boundary
2533         DO i = max(its,ibe-spec_zone+1), itf
2534           b_dist = ibe - i
2535           DO k = kts, ktf
2536             DO j = max(jts,b_dist+jbs+1), min(jtf,jbe-b_dist-1)
2537               j_inner = max(j,jbs+spec_zone)
2538               j_inner = min(j_inner,jbe-spec_zone)
2539               IF(u(i+1,k,j) .gt. 0.)THEN
2540                 field(i,k,j) = field(ibe-spec_zone,k,j_inner)
2541               ELSE
2542                 field(i,k,j) = value_bc
2543               ENDIF
2544             ENDDO
2545           ENDDO
2546         ENDDO
2547       ENDIF
2548     ENDIF
2549
2550    END SUBROUTINE flow_dep_bdy_s1
2551
2552 !------------------------------------------------------------------------------
2553
2554    SUBROUTINE flow_dep_bdy_s2 (  field,                     &
2555                                u, v, config_flags, &
2556                                spec_zone,                  &
2557                                ids,ide, jds,jde, kds,kde,  & ! domain dims
2558                                ims,ime, jms,jme, kms,kme,  & ! memory dims
2559                                ips,ipe, jps,jpe, kps,kpe,  & ! patch  dims
2560                                its,ite, jts,jte, kts,kte, dtstep, ktau)
2561
2562 !  This subroutine sets zero gradient conditions for outflow and zero value
2563 !  for inflow in the boundary specified region. Note that field must be unstaggered.
2564 !  The velocities, u and v, will only be used to check their sign (coupled vels OK)
2565 !  spec_zone is the width of the outer specified b.c.s that are set here.
2566 !  (JD August 2000)
2567
2568       IMPLICIT NONE
2569
2570       INTEGER,      INTENT(IN   )    :: ids,ide, jds,jde, kds,kde
2571       INTEGER,      INTENT(IN   )    :: ims,ime, jms,jme, kms,kme
2572       INTEGER,      INTENT(IN   )    :: ips,ipe, jps,jpe, kps,kpe
2573       INTEGER,      INTENT(IN   )    :: its,ite, jts,jte, kts,kte
2574       INTEGER,      INTENT(IN   )    :: spec_zone
2575
2576
2577       REAL,  DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(INOUT) :: field
2578       REAL,  DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(IN   ) :: u
2579       REAL,  DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(IN   ) :: v
2580       TYPE( grid_config_rec_type ) config_flags
2581
2582       INTEGER    :: i, j, k, ibs, ibe, jbs, jbe, itf, jtf, ktf, i_inner, j_inner
2583       INTEGER    :: b_dist, b_limit
2584       LOGICAL    :: periodic_x
2585 !----------------------------------------------
2586 !jeff
2587 !--   hardwire outmost lateral boundary value with 1-day efolding decay
2588
2589       real,    INTENT(IN   ) :: dtstep
2590       integer, INTENT(IN   ) :: ktau
2591
2592       real value_bc
2593       real factor_decay
2594
2595 !--  initial value
2596
2597       value_bc = 1.0
2598
2599 !-- decay factor, with efolding time of one day
2600
2601       factor_decay = 1./(86400./dtstep)
2602
2603 !-----------------------------------------------
2604
2605       periodic_x = config_flags%periodic_x
2606
2607       ibs = ids
2608       ibe = ide-1
2609       itf = min(ite,ide-1)
2610       jbs = jds
2611       jbe = jde-1
2612       jtf = min(jte,jde-1)
2613       ktf = kde-1
2614
2615       IF (jts - jbs .lt. spec_zone) THEN
2616 ! Y-start boundary
2617         DO j = jts, min(jtf,jbs+spec_zone-1)
2618           b_dist = j - jbs
2619           b_limit = b_dist
2620           IF(periodic_x)b_limit = 0
2621           DO k = kts, ktf
2622             DO i = max(its,b_limit+ibs), min(itf,ibe-b_limit)
2623               i_inner = max(i,ibs+spec_zone)
2624               i_inner = min(i_inner,ibe-spec_zone)
2625               IF(periodic_x)i_inner = i
2626               IF(v(i,k,j) .lt. 0.)THEN
2627                 field(i,k,j) = field(i_inner,k,jbs+spec_zone)
2628               ELSE
2629                 if (ktau .eq. 1) then
2630                    field(i,k,j) = value_bc
2631                 else
2632                    field(i,k,j) = field(i,k,j) * (1. - factor_decay)
2633                 endif
2634               ENDIF
2635             ENDDO
2636           ENDDO
2637         ENDDO
2638       ENDIF
2639       IF (jbe - jtf .lt. spec_zone) THEN
2640 ! Y-end boundary
2641         DO j = max(jts,jbe-spec_zone+1), jtf
2642           b_dist = jbe - j
2643           b_limit = b_dist
2644           IF(periodic_x)b_limit = 0
2645           DO k = kts, ktf
2646             DO i = max(its,b_limit+ibs), min(itf,ibe-b_limit)
2647               i_inner = max(i,ibs+spec_zone)
2648               i_inner = min(i_inner,ibe-spec_zone)
2649               IF(periodic_x)i_inner = i
2650               IF(v(i,k,j+1) .gt. 0.)THEN
2651                 field(i,k,j) = field(i_inner,k,jbe-spec_zone)
2652               ELSE
2653                 if (ktau .eq. 1) then
2654                    field(i,k,j) = value_bc
2655                 else
2656                    field(i,k,j) = field(i,k,j) * (1. - factor_decay)
2657                 endif
2658               ENDIF
2659             ENDDO
2660           ENDDO
2661         ENDDO
2662       ENDIF
2663
2664     IF(.NOT.periodic_x)THEN
2665       IF (its - ibs .lt. spec_zone) THEN
2666 ! X-start boundary
2667         DO i = its, min(itf,ibs+spec_zone-1)
2668           b_dist = i - ibs
2669           DO k = kts, ktf
2670             DO j = max(jts,b_dist+jbs+1), min(jtf,jbe-b_dist-1)
2671               j_inner = max(j,jbs+spec_zone)
2672               j_inner = min(j_inner,jbe-spec_zone)
2673               IF(u(i,k,j) .lt. 0.)THEN
2674                 field(i,k,j) = field(ibs+spec_zone,k,j_inner)
2675               ELSE
2676                 if (ktau .eq. 1) then
2677                    field(i,k,j) = value_bc
2678                 else
2679                    field(i,k,j) = field(i,k,j) * (1. - factor_decay)
2680                 endif
2681               ENDIF
2682             ENDDO
2683           ENDDO
2684         ENDDO
2685       ENDIF
2686
2687       IF (ibe - itf .lt. spec_zone) THEN
2688 ! X-end boundary
2689         DO i = max(its,ibe-spec_zone+1), itf
2690           b_dist = ibe - i
2691           DO k = kts, ktf
2692             DO j = max(jts,b_dist+jbs+1), min(jtf,jbe-b_dist-1)
2693               j_inner = max(j,jbs+spec_zone)
2694               j_inner = min(j_inner,jbe-spec_zone)
2695               IF(u(i+1,k,j) .gt. 0.)THEN
2696                 field(i,k,j) = field(ibe-spec_zone,k,j_inner)
2697               ELSE
2698                 if (ktau .eq. 1) then
2699                    field(i,k,j) = value_bc
2700                 else
2701                    field(i,k,j) = field(i,k,j) * (1. - factor_decay)
2702                 endif
2703               ENDIF
2704             ENDDO
2705           ENDDO
2706         ENDDO
2707       ENDIF
2708     ENDIF
2709
2710    END SUBROUTINE flow_dep_bdy_s2
2711
2712 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2713   SUBROUTINE bdy_chem_value_gcm ( chem, chem_b, chem_bt, dt,ic)
2714
2715     IMPLICIT NONE
2716
2717     REAL,    intent(OUT)  :: chem
2718     REAL,    intent(IN)   :: chem_b
2719     REAL,    intent(IN)   :: chem_bt
2720     REAL,    intent(IN)   :: dt
2721     INTEGER, intent(IN)   :: ic
2722
2723
2724     CHARACTER (LEN=80) :: message
2725 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2726
2727 !     if( nch .GT. numgas) then
2728 !       message = ' Input_chem_profile: wrong number of chemical species'
2729 !       return
2730 !       CALL WRF_ERROR_FATAL ( message )
2731 !     endif
2732
2733
2734       !print*,'before',chem,chem_bt ,dt, ic
2735
2736       chem=max(epsilc,chem_b + chem_bt * dt)
2737       !print*,'after',chem
2738       RETURN
2739   END SUBROUTINE bdy_chem_value_gcm
2740 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2741
2742 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2743    SUBROUTINE cv_mmdd_jday ( YY, MM, DD, JDAY)
2744 !
2745 !     Subroutine to compute the julian day given the month and day
2746 !
2747 !
2748       INTEGER,      INTENT(IN )    :: YY, MM, DD
2749       INTEGER,      INTENT(OUT)    :: JDAY
2750
2751       INTEGER, DIMENSION(12) :: imon, imon_a
2752       INTEGER                :: i
2753
2754       DATA imon_a /0,31,59,90,120,151,181,212,243,273,304,334/
2755 !
2756 !..... Check for leap year.
2757 !
2758       do i=1,12
2759          imon(i) = imon_a(i)
2760       enddo
2761       if(YY .eq. (YY/4)*4) then
2762          do i=3,12
2763             imon(i) = imon(i) + 1
2764          enddo
2765       endif
2766 !
2767 !..... Convert month, day to julian day.
2768 !
2769       jday = imon(mm) + dd
2770
2771
2772    END SUBROUTINE cv_mmdd_jday
2773
2774 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2775
2776    integer FUNCTION get_last_gas(chem_opt)
2777      implicit none
2778      integer, intent(in) :: chem_opt
2779
2780      ! Determine the index of the last gas species, which depends
2781      ! upon the gas mechanism.
2782
2783      select case (chem_opt)
2784      case (0)
2785         get_last_gas = 0
2786
2787      case (RADM2, RADM2_KPP, RADM2SORG, RADM2SORG_AQ, RADM2SORG_KPP, &
2788            RACM_KPP, RACMPM_KPP, RACM_MIM_KPP, RACMSORG_AQ, RACMSORG_KPP, &
2789                GOCARTRACM_KPP,GOCARTRADM2,GOCARTRADM2_KPP)
2790         get_last_gas = p_ho2
2791
2792      case (CBMZ)
2793         get_last_gas = p_mtf
2794
2795      case (CBMZ_BB,CBMZ_BB_KPP,CBMZ_MOSAIC_4BIN,CBMZ_MOSAIC_8BIN,CBMZ_MOSAIC_4BIN_AQ,CBMZ_MOSAIC_8BIN_AQ)
2796         get_last_gas = p_isopo2
2797
2798      case (CHEM_TRACER)
2799         get_last_gas = p_co
2800
2801      case (CHEM_TRACE2)
2802         get_last_gas = p_tracer_1
2803
2804      case (GOCART_SIMPLE)
2805         get_last_gas = p_msa
2806
2807      case (CBM4_KPP)
2808         get_last_gas = p_ho2
2809
2810      case (MOZART_KPP)
2811         get_last_gas = p_meko2
2812
2813      case (MOZCART_KPP)
2814         get_last_gas = p_meko2
2815
2816      case default
2817         call wrf_error_fatal("get_last_gas: could not decipher chem_opt value")
2818
2819      end select
2820
2821    END FUNCTION get_last_gas
2822 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2823
2824 !****************************************************************
2825 !                                                               *
2826 !   SUBROUTINE TO SET AEROSOL BC VALUES USING THE               *
2827 !   aer_bc_opt == aer_bc_pnnl OPTION.                           *
2828 !                                                               *
2829 !   wig 22-Apr-2004, original routine                           *
2830 !       rce 25-apr-2004 - changed name to                       *
2831 !                         "sorgam_set_aer_bc_pnnl"              *
2832 !       wig  7-May-2004, added height dependance                *
2833 !                                                               *
2834 !   CALLS THE FOLLOWING SUBROUTINES:  NONE                      *
2835 !                                                               *
2836 !   CALLED BY:                        bdy_chem_value_sorgam     *
2837 !                                                               *
2838 !****************************************************************
2839     SUBROUTINE sorgam_set_aer_bc_pnnl( chem, z, nch, config_flags )
2840       USE module_data_sorgam, ONLY : dginia, dginin, dginic, esn36, esc36, esa36, seasfac, no3fac, nh4fac, so4fac, soilfac, anthfac, orgfac
2841
2842       implicit none
2843
2844       INTEGER,INTENT(IN   ) :: nch
2845       real,intent(in      ) :: z
2846       REAL,INTENT(INOUT   ) :: chem
2847       TYPE (grid_config_rec_type) , INTENT (in) :: config_flags
2848
2849       REAL :: mult,                       &
2850               m3acc, m3cor, m3nuc,        &
2851               bv_so4ai, bv_so4aj,         &
2852               bv_nh4ai, bv_nh4aj,         &
2853               bv_no3ai, bv_no3aj,         &
2854               bv_eci,   bv_ecj,           &
2855               bv_p25i,  bv_p25j,          &
2856               bv_orgpai,bv_orgpaj,        &
2857               bv_antha, bv_seas, bv_soila
2858
2859 !
2860 ! Determine height multiplier...
2861 ! This should mimic the calculation in sorgam_init_aer_ic_pnnl,
2862 ! mosaic_init_wrf_mixrats_opt2, and bdy_chem_value_mosaic
2863 !!$!    Height(m)     Multiplier
2864 !!$!    ---------     ----------
2865 !!$!    <=2000        1.0
2866 !!$!    2000<z<3000   linear transition zone to 0.5
2867 !!$!    3000<z<5000   linear transision zone to 0.25
2868 !!$!    >=5000        0.25
2869 !!$!
2870 !!$! which translates to:
2871 !!$!    2000<z<3000   mult = 1.0 + (z-2000.)*(0.5-1.0)/(3000.-2000.)
2872 !!$!    3000<z<5000   mult = 0.5 + (z-3000.)*(0.25-0.5)/(5000.-3000.)
2873 !!$! or in reduced form:
2874 !!$      if( z <= 2000. ) then
2875 !!$         mult = 1.0
2876 !!$      elseif( z > 2000. &
2877 !!$           .and. z <= 3000. ) then
2878 !!$         mult = 1.0 - 0.0005*(z-2000.)
2879 !!$      elseif( z > 3000. &
2880 !!$           .and. z <= 5000. ) then
2881 !!$         mult = 0.5 - 1.25e-4*(z-3000.)
2882 !!$      else
2883 !!$         mult = 0.25
2884 !!$      end if
2885 ! Updated aerosol profile multiplier 1-Apr-2005:
2886 !    Height(m)     Multiplier
2887 !    ---------     ----------
2888 !    <=2000        1.0
2889 !    2000<z<3000   linear transition zone to 0.25
2890 !    3000<z<5000   linear transision zone to 0.125
2891 !    >=5000        0.125
2892 !
2893 ! which translates to:
2894 !    2000<z<3000   mult = 1.00 + (z-2000.)*(0.25-1.0)/(3000.-2000.)
2895 !    3000<z<5000   mult = 0.25 + (z-3000.)*(0.125-0.25)/(5000.-3000.)
2896 ! or in reduced form:
2897 !       if( z <= 2000. ) then
2898 !          mult = 1.0
2899 !       elseif( z > 2000. &
2900 !            .and. z <= 3000. ) then
2901 !          mult = 1.0 - 0.00075*(z-2000.)
2902 !       elseif( z > 3000. &
2903 !            .and. z <= 5000. ) then
2904 !          mult = 0.25 - 4.166666667e-5*(z-3000.)
2905 !       else
2906 !          mult = 0.125
2907 !       end if
2908         if( z <= 500. ) then
2909            mult = 1.0
2910         elseif( z > 500. &
2911              .and. z <= 1000. ) then
2912            mult = 1.0 - 0.001074*(z-500.)
2913         elseif( z > 1000. &
2914              .and. z <= 5000. ) then
2915            mult = 0.463 - 0.000111*(z-1000.)
2916         else
2917            mult = 0.019
2918         end if
2919
2920 ! These should match what is in sorgam_init_aer_ic_pnnl.
2921 ! Boundary values as of 2-Dec-2004:
2922 !     bv_so4aj  = mult*2.375
2923 !     bv_so4ai  = mult*0.179
2924 !     bv_nh4aj  = mult*0.9604
2925 !     bv_nh4ai  = mult*0.0196
2926 !     bv_no3aj  = mult*0.0650
2927 !     bv_no3ai  = mult*0.0050
2928 !     bv_ecj    = mult*0.1630
2929 !     bv_eci    = mult*0.0120
2930 !     bv_p25j   = mult*0.6350
2931 !     bv_p25i   = mult*0.0490
2932 !     bv_orgpaj = mult*0.9300
2933 !     bv_orgpai = mult*0.0700
2934 !     bv_antha  = mult*2.2970
2935 !     bv_seas   = mult*0.2290
2936 !     bv_soila  = conmin
2937       bv_so4aj = mult*0.300*0.97
2938       bv_so4ai = mult*0.300*0.03
2939       bv_nh4aj = mult*0.094*0.97
2940       bv_nh4ai = mult*0.094*0.03
2941       bv_no3aj = mult*0.001*0.97
2942       bv_no3ai = mult*0.001*0.03
2943       bv_ecj   = mult*0.013*0.97
2944       bv_eci   = mult*0.013*0.03
2945       bv_p25j  = mult*4.500*0.97
2946       bv_p25i  = mult*4.500*0.03
2947       bv_antha = mult*4.500/2.0
2948       bv_orgpaj = mult*0.088*0.97
2949       bv_orgpai = mult*0.088*0.03
2950       bv_seas   = mult*1.75
2951       bv_soila  = conmin
2952
2953 ! m3... calculations should match the very end of module_aerosols_sorgam.F
2954 !... i-mode (note that the 8 SOA species have bv=conmin)
2955       m3nuc = so4fac*bv_so4ai + nh4fac*bv_nh4ai + &
2956         no3fac*bv_no3ai + &
2957         orgfac*8.0*conmin + orgfac*bv_orgpai + &
2958         anthfac*bv_p25i + anthfac*bv_eci
2959
2960 !... j-mode (note that the 8 SOA species have bv=conmin)
2961       m3acc = so4fac*bv_so4aj + nh4fac*bv_nh4aj + &
2962         no3fac*bv_no3aj + &
2963         orgfac*8.0*conmin + orgfac*bv_orgpaj + &
2964         anthfac*bv_p25j + anthfac*bv_ecj
2965
2966 !...c-mode
2967       m3cor = soilfac*bv_soila + seasfac*bv_seas + &
2968         anthfac*bv_antha
2969
2970 ! Cannot set_sulf here because it is a "radm2" species whose bc value
2971 ! is set via bdy_chem_value. Instead, xl(iref(p_sulf-1),:) is set to
2972 ! the value conmin in subroutine gasprofile_init_pnnl
2973 !      if( nch == p_sulf    ) chem = conmin !as per rce's 0 recommendation
2974
2975       if( nch == p_so4aj   ) chem = bv_so4aj
2976       if( nch == p_so4ai   ) chem = bv_so4ai
2977       if( nch == p_nh4aj   ) chem = bv_nh4aj
2978       if( nch == p_nh4ai   ) chem = bv_nh4ai
2979       if( nch == p_no3aj   ) chem = bv_no3aj
2980       if( nch == p_no3ai   ) chem = bv_no3ai
2981       if( nch == p_ecj     ) chem = bv_ecj
2982       if( nch == p_eci     ) chem = bv_eci
2983       if( nch == p_p25j    ) chem = bv_p25j
2984       if( nch == p_p25i    ) chem = bv_p25i
2985       if( nch == p_orgpaj  ) chem = bv_orgpaj
2986       if( nch == p_orgpai  ) chem = bv_orgpai
2987
2988       if( nch == p_orgaro1j) chem = conmin
2989       if( nch == p_orgaro1i) chem = conmin
2990       if( nch == p_orgaro2j) chem = conmin
2991       if( nch == p_orgaro2i) chem = conmin
2992       if( nch == p_orgalk1j) chem = conmin
2993       if( nch == p_orgalk1i) chem = conmin
2994       if( nch == p_orgole1j) chem = conmin
2995       if( nch == p_orgole1i) chem = conmin
2996       if( nch == p_orgba1j ) chem = conmin
2997       if( nch == p_orgba1i ) chem = conmin
2998       if( nch == p_orgba2j ) chem = conmin
2999       if( nch == p_orgba2i ) chem = conmin
3000       if( nch == p_orgba3j ) chem = conmin
3001       if( nch == p_orgba3i ) chem = conmin
3002       if( nch == p_orgba4j ) chem = conmin
3003       if( nch == p_orgba4i ) chem = conmin
3004
3005       if( nch == p_antha   ) chem = bv_antha
3006       if( nch == p_soila   ) chem = bv_soila
3007       if( nch == p_seas    ) chem = bv_seas
3008
3009       if( nch == p_nu0     ) chem = m3nuc/((dginin**3)*esn36)
3010       if( nch == p_ac0     ) chem = m3acc/((dginia**3)*esa36)
3011       if( nch == p_corn    ) chem = m3cor/((dginic**3)*esc36)
3012
3013     END SUBROUTINE sorgam_set_aer_bc_pnnl
3014 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3015
3016 !****************************************************************
3017 !                                                               *
3018 !   SUBROUTINE TO OVERWRITE THE PREDEFINED OZONE PROFILE        *
3019 !   WHEN gas_ic_opt == gas_ic_pnnl                              *
3020 !         OR gas_bc_opt == gas_bc_pnnl                              *
3021 !                                                               *
3022 !   wig, 21-Apr-2004                                            *
3023 !       rce 25-apr-2004 - changed name to                       *
3024 !                         "gasprofile_init_pnnl"                *
3025 !                                                               *
3026 !   CALLS THE FOLLOWING SUBROUTINES:  NONE                      *
3027 !                                                               *
3028 !   CALLED BY:                        chem_init                 *
3029 !                                     input_chem_profile        *
3030 !                                                               *
3031 !****************************************************************
3032     SUBROUTINE gasprofile_init_pnnl( chem_opt )
3033       use module_data_sorgam,only:  conmin
3034       implicit none
3035       INTEGER, INTENT (in) :: chem_opt
3036
3037       integer :: k
3038
3039       call wrf_debug ( 500 , 'wrfchem:gasprofile_init_pnnl' )
3040 !     print*,'gasprofile_init_pnnl redefining o3 and sulf profiles.'
3041
3042 ! Original O3 profile values:
3043 !      / 1.68E-07, 1.68E-07, 5.79E-08, 5.24E-08, 5.26E-08, &
3044 !       5.16E-08, 4.83E-08, 4.50E-08, 4.16E-08, 3.80E-08, 3.56E-08, &
3045 !       3.35E-08, 3.15E-08, 3.08E-08, 3.06E-08, 3.00E-08/
3046
3047 ! Note that heights associated with 2nd index of xl correspond to upside-down
3048 ! zfa values that have been "de-staggered".
3049 ! Height = 0.5*(zfa(1:kx) + zfa(2:kx+1)) and then flipped:
3050 !      / 18500., 14050., 11150., 9355., 7705., 6120., 4675., 3430.,
3051 !         2430.,  1720.,  1195., 781.5,  494., 298.5, 148.5,  42.5
3052
3053 #if (CASENAME != 4)
3054
3055 #if (CASENAME == 0)
3056       if( p_o3 > 1 ) then
3057                                           !Rounded to closest level:
3058          xl(iref(p_o3-1),11:16) = 4.00e-8 !40 ppbv below 1 km
3059          xl(iref(p_o3-1),3:10)  = 6.50e-8 !65 ppbv > 2 km and < stratosphere
3060                                           !  Changed from 70 ppbv 1-Apr-2005
3061       end if
3062 #endif
3063
3064 #if (CASENAME == 1)
3065       if( p_o3 > 1 ) then
3066          xl(iref(p_o3-1),11:16) = 3.50e-8 !35 ppbv below 1 km
3067          xl(iref(p_o3-1),3:10)  = 6.00e-8 !60 ppbv > 2 km and < stratosphere
3068       end if
3069 #endif
3070 #if (CASENAME == 2)
3071       if( p_o3 > 1 ) then
3072          xl(iref(p_o3-1),12:16) = 2.50e-8
3073          xl(iref(p_o3-1),11:11) = 3.10e-8
3074          xl(iref(p_o3-1),10:10) = 3.20e-8
3075          xl(iref(p_o3-1),9:9)   = 3.60e-8
3076          xl(iref(p_o3-1),8:8)   = 4.20e-8
3077          xl(iref(p_o3-1),7:7)   = 4.20e-8
3078          xl(iref(p_o3-1),6:6)   = 4.20e-8
3079          xl(iref(p_o3-1),5:5)   = 3.80e-8
3080          xl(iref(p_o3-1),4:4)   = 3.80e-8
3081          xl(iref(p_o3-1),3:3)   = 5.80e-8
3082          xl(iref(p_o3-1),2:2)   = 1.04e-7
3083          xl(iref(p_o3-1),1:1)   = 2.00e-7
3084       end if
3085       if( p_co > 1 ) then
3086          xl(iref(p_co-1),12:16) = 1.00e-7
3087          xl(iref(p_co-1),11:11) = 1.00e-7
3088          xl(iref(p_co-1),10:10) = 9.80e-8
3089          xl(iref(p_co-1),9:9)   = 8.80e-8
3090          xl(iref(p_co-1),8:8)   = 7.80e-8
3091          xl(iref(p_co-1),7:7)   = 7.70e-8
3092          xl(iref(p_co-1),6:6)   = 8.00e-8
3093          xl(iref(p_co-1),5:5)   = 8.20e-8
3094          xl(iref(p_co-1),4:4)   = 8.20e-8
3095          xl(iref(p_co-1),3:3)   = 8.20e-8
3096          xl(iref(p_co-1),2:2)   = 6.50e-8
3097          xl(iref(p_co-1),1:1)   = 6.10e-8
3098       end if
3099       if( p_so2 > 1 ) then
3100          xl(iref(p_so2-1),12:16) = 6.00e-11
3101          xl(iref(p_so2-1),11:11) = 1.00e-11
3102          xl(iref(p_so2-1),10:10) = 4.00e-12
3103          xl(iref(p_so2-1),9:9)   = 5.00e-12
3104          xl(iref(p_so2-1),8:8)   = 2.00e-12
3105          xl(iref(p_so2-1),7:7)   = 2.00e-12
3106          xl(iref(p_so2-1),6:6)   = 5.00e-12
3107          xl(iref(p_so2-1),5:5)   = 4.00e-12
3108          xl(iref(p_so2-1),4:4)   = 5.00e-12
3109          xl(iref(p_so2-1),3:3)   = 1.00e-11
3110          xl(iref(p_so2-1),2:2)   = 2.00e-12
3111          xl(iref(p_so2-1),1:1)   = 5.00e-13
3112       end if
3113       if( p_no2 > 1 ) then
3114          xl(iref(p_no2-1),12:16) = 1.00e-11
3115          xl(iref(p_no2-1),11:11) = 1.50e-11
3116          xl(iref(p_no2-1),10:10) = 1.80e-11
3117          xl(iref(p_no2-1),9:9)   = 1.90e-11
3118          xl(iref(p_no2-1),8:8)   = 2.80e-11
3119          xl(iref(p_no2-1),7:7)   = 3.20e-11
3120          xl(iref(p_no2-1),6:6)   = 2.60e-11
3121          xl(iref(p_no2-1),5:5)   = 1.80e-11
3122          xl(iref(p_no2-1),4:4)   = 1.70e-11
3123          xl(iref(p_no2-1),3:3)   = 4.30e-11
3124          xl(iref(p_no2-1),2:2)   = 2.50e-10
3125          xl(iref(p_no2-1),1:1)   = 2.80e-10
3126       end if
3127       if( p_h2o2 > 1 ) then
3128          xl(iref(p_h2o2-1),12:16) = 0.18e-8
3129          xl(iref(p_h2o2-1),11:11) = 0.13e-8
3130          xl(iref(p_h2o2-1),10:10) = 0.30e-9
3131          xl(iref(p_h2o2-1),9:9)   = 0.77e-9
3132          xl(iref(p_h2o2-1),8:8)   = 0.30e-9
3133          xl(iref(p_h2o2-1),7:7)   = 0.25e-9
3134          xl(iref(p_h2o2-1),6:6)   = 0.21e-9
3135          xl(iref(p_h2o2-1),5:5)   = 0.23e-9
3136          xl(iref(p_h2o2-1),4:4)   = 0.24e-9
3137          xl(iref(p_h2o2-1),3:3)   = 0.36e-9
3138          xl(iref(p_h2o2-1),2:2)   = 0.15e-10
3139          xl(iref(p_h2o2-1),1:1)   = 0.17e-10
3140       end if
3141       if( p_hno3 > 1 ) then
3142          xl(iref(p_hno3-1),12:16) = 0.98e-11
3143          xl(iref(p_hno3-1),11:11) = 0.58e-11
3144          xl(iref(p_hno3-1),10:10) = 0.53e-11
3145          xl(iref(p_hno3-1),9:9)   = 0.53e-11
3146          xl(iref(p_hno3-1),8:8)   = 0.67e-11
3147          xl(iref(p_hno3-1),7:7)   = 0.82e-11
3148          xl(iref(p_hno3-1),6:6)   = 0.82e-11
3149          xl(iref(p_hno3-1),5:5)   = 0.21e-11
3150          xl(iref(p_hno3-1),4:4)   = 0.39e-11
3151          xl(iref(p_hno3-1),3:3)   = 0.57e-10
3152          xl(iref(p_hno3-1),2:2)   = 0.17e-9
3153          xl(iref(p_hno3-1),1:1)   = 0.12e-9
3154       end if
3155       if( p_no3 > 1 ) then
3156          xl(iref(p_no3-1),12:16) = 0.98e-11
3157          xl(iref(p_no3-1),11:11) = 0.58e-11
3158          xl(iref(p_no3-1),10:10) = 0.53e-11
3159          xl(iref(p_no3-1),9:9)   = 0.53e-11
3160          xl(iref(p_no3-1),8:8)   = 0.68e-11
3161          xl(iref(p_no3-1),7:7)   = 0.25e-11
3162          xl(iref(p_no3-1),6:6)   = 0.82e-11
3163          xl(iref(p_no3-1),5:5)   = 0.21e-11
3164          xl(iref(p_no3-1),4:4)   = 0.39e-11
3165          xl(iref(p_no3-1),3:3)   = 0.57e-10
3166          xl(iref(p_no3-1),2:2)   = 0.17e-9
3167          xl(iref(p_no3-1),1:1)   = 0.12e-9
3168       end if
3169       if( p_n2o5 > 1 ) then
3170          xl(iref(p_n2o5-1),12:16) = 0.12e-14
3171          xl(iref(p_n2o5-1),11:11) = 0.11e-12
3172          xl(iref(p_n2o5-1),10:10) = 0.20e-12
3173          xl(iref(p_n2o5-1),9:9)   = 0.25e-12
3174          xl(iref(p_n2o5-1),8:8)   = 0.98e-12
3175          xl(iref(p_n2o5-1),7:7)   = 0.12e-11
3176          xl(iref(p_n2o5-1),6:6)   = 0.13e-11
3177          xl(iref(p_n2o5-1),5:5)   = 0.45e-12
3178          xl(iref(p_n2o5-1),4:4)   = 0.38e-12
3179          xl(iref(p_n2o5-1),3:3)   = 0.94e-12
3180          xl(iref(p_n2o5-1),2:2)   = 0.16e-11
3181          xl(iref(p_n2o5-1),1:1)   = 0.11e-11
3182       end if
3183       if( p_nh3 > 1 ) then
3184          xl(iref(p_nh3-1),12:16) = 0.14e-9
3185          xl(iref(p_nh3-1),11:11) = 0.60e-10
3186          xl(iref(p_nh3-1),10:10) = 0.30e-10
3187          xl(iref(p_nh3-1),9:9)   = 0.24e-10
3188          xl(iref(p_nh3-1),8:8)   = 0.49e-11
3189          xl(iref(p_nh3-1),7:7)   = 0.19e-11
3190          xl(iref(p_nh3-1),6:6)   = 0.41e-11
3191          xl(iref(p_nh3-1),5:5)   = 0.11e-11
3192          xl(iref(p_nh3-1),4:4)   = 0.52e-11
3193          xl(iref(p_nh3-1),3:3)   = 0.35e-10
3194          xl(iref(p_nh3-1),2:2)   = 0.86e-11
3195          xl(iref(p_nh3-1),1:1)   = 0.84e-11
3196       end if
3197       if( p_hcho > 1 ) then
3198          xl(iref(p_hcho-1),12:16) = 4.00e-10
3199          xl(iref(p_hcho-1),11:11) = 4.20e-10
3200          xl(iref(p_hcho-1),10:10) = 3.60e-10
3201          xl(iref(p_hcho-1),9:9)   = 3.00e-10
3202          xl(iref(p_hcho-1),8:8)   = 1.80e-10
3203          xl(iref(p_hcho-1),7:7)   = 1.40e-10
3204          xl(iref(p_hcho-1),6:6)   = 7.60e-11
3205          xl(iref(p_hcho-1),5:5)   = 6.20e-11
3206          xl(iref(p_hcho-1),4:4)   = 5.20e-11
3207          xl(iref(p_hcho-1),3:3)   = 8.80e-11
3208          xl(iref(p_hcho-1),2:2)   = 1.30e-11
3209          xl(iref(p_hcho-1),1:1)   = 1.40e-11
3210       end if
3211       if( p_ket > 1 ) then
3212          xl(iref(p_ket-1),12:16) = 0.19e-9
3213          xl(iref(p_ket-1),11:11) = 0.15e-9
3214          xl(iref(p_ket-1),10:10) = 0.13e-9
3215          xl(iref(p_ket-1),9:9)   = 0.98e-10
3216          xl(iref(p_ket-1),8:8)   = 0.54e-10
3217          xl(iref(p_ket-1),7:7)   = 0.46e-10
3218          xl(iref(p_ket-1),6:6)   = 0.66e-10
3219          xl(iref(p_ket-1),5:5)   = 0.64e-10
3220          xl(iref(p_ket-1),4:4)   = 0.68e-10
3221          xl(iref(p_ket-1),3:3)   = 0.63e-10
3222          xl(iref(p_ket-1),2:2)   = 0.91e-11
3223          xl(iref(p_ket-1),1:1)   = 0.55e-11
3224       end if
3225       if( p_par > 1 ) then
3226          xl(iref(p_par-1),12:16) = 0.75e-10+0.68e-10
3227          xl(iref(p_par-1),11:11) = 0.50e-10+0.21e-10
3228          xl(iref(p_par-1),10:10) = 0.41e-10+0.11e-10
3229          xl(iref(p_par-1),9:9)   = 0.37e-10+0.17e-10
3230          xl(iref(p_par-1),8:8)   = 0.24e-10+0.53e-11
3231          xl(iref(p_par-1),7:7)   = 0.22e-10+0.45e-11
3232          xl(iref(p_par-1),6:6)   = 0.31e-10+0.79e-11
3233          xl(iref(p_par-1),5:5)   = 0.28e-10+0.22e-11
3234          xl(iref(p_par-1),4:4)   = 0.29e-10+0.26e-11
3235          xl(iref(p_par-1),3:3)   = 0.29e-10+0.43e-11
3236          xl(iref(p_par-1),2:2)   = 0.16e-10+0.18e-12
3237          xl(iref(p_par-1),1:1)   = 0.17e-10+0.26e-12
3238       end if
3239       if( p_olt > 1 ) then
3240          xl(iref(p_olt-1),12:16) = 0.54e-13
3241          xl(iref(p_olt-1),11:11) = 0.52e-13
3242          xl(iref(p_olt-1),10:10) = 0.33e-13
3243          xl(iref(p_olt-1),9:9)   = 0.16e-13
3244          xl(iref(p_olt-1),8:8)   = 0.20e-14
3245          xl(iref(p_olt-1),7:7)   = 0.14e-14
3246          xl(iref(p_olt-1),6:6)   = 0.37e-14
3247          xl(iref(p_olt-1),5:5)   = 0.26e-15
3248          xl(iref(p_olt-1),4:4)   = 0.40e-15
3249          xl(iref(p_olt-1),3:3)   = 0.15e-14
3250          xl(iref(p_olt-1),2:2)   = 0.33e-17
3251          xl(iref(p_olt-1),1:1)   = 0.27e-15
3252       end if
3253       if( p_ol2 > 1 ) then
3254          xl(iref(p_ol2-1),12:16) = 0.87e-12
3255          xl(iref(p_ol2-1),11:11) = 0.10e-11
3256          xl(iref(p_ol2-1),10:10) = 0.68e-12
3257          xl(iref(p_ol2-1),9:9)   = 0.38e-12
3258          xl(iref(p_ol2-1),8:8)   = 0.38e-12
3259          xl(iref(p_ol2-1),7:7)   = 0.45e-13
3260          xl(iref(p_ol2-1),6:6)   = 0.58e-13
3261          xl(iref(p_ol2-1),5:5)   = 0.52e-14
3262          xl(iref(p_ol2-1),4:4)   = 0.12e-13
3263          xl(iref(p_ol2-1),3:3)   = 0.66e-13
3264          xl(iref(p_ol2-1),2:2)   = 0.14e-14
3265          xl(iref(p_ol2-1),1:1)   = 0.51e-13
3266       end if
3267       if( p_ald > 1 ) then
3268          xl(iref(p_ald-1),12:16) = 0.30e-10
3269          xl(iref(p_ald-1),11:11) = 0.28e-10
3270          xl(iref(p_ald-1),10:10) = 0.22e-10
3271          xl(iref(p_ald-1),9:9)   = 0.16e-10
3272          xl(iref(p_ald-1),8:8)   = 0.10e-10
3273          xl(iref(p_ald-1),7:7)   = 0.83e-11
3274          xl(iref(p_ald-1),6:6)   = 0.89e-11
3275          xl(iref(p_ald-1),5:5)   = 0.93e-11
3276          xl(iref(p_ald-1),4:4)   = 0.98e-11
3277          xl(iref(p_ald-1),3:3)   = 0.10e-10
3278          xl(iref(p_ald-1),2:2)   = 0.93e-12
3279          xl(iref(p_ald-1),1:1)   = 0.64e-12
3280       end if
3281       if( p_tol > 1 ) then
3282          xl(iref(p_tol-1),12:16) = 0.29e-11/2.0
3283          xl(iref(p_tol-1),11:11) = 0.16e-11/2.0
3284          xl(iref(p_tol-1),10:10) = 0.12e-11/2.0
3285          xl(iref(p_tol-1),9:9)   = 0.97e-12/2.0
3286          xl(iref(p_tol-1),8:8)   = 0.16e-12/2.0
3287          xl(iref(p_tol-1),7:7)   = 0.13e-12/2.0
3288          xl(iref(p_tol-1),6:6)   = 0.13e-12/2.0
3289          xl(iref(p_tol-1),5:5)   = 0.28e-13/2.0
3290          xl(iref(p_tol-1),4:4)   = 0.42e-13/2.0
3291          xl(iref(p_tol-1),3:3)   = 0.12e-12/2.0
3292          xl(iref(p_tol-1),2:2)   = 0.39e-15/2.0
3293          xl(iref(p_tol-1),1:1)   = 0.11e-14/2.0
3294       end if
3295       if( p_xyl > 1 ) then
3296          xl(iref(p_xyl-1),12:16) = 0.29e-11/2.0
3297          xl(iref(p_xyl-1),11:11) = 0.16e-11/2.0
3298          xl(iref(p_xyl-1),10:10) = 0.12e-11/2.0
3299          xl(iref(p_xyl-1),9:9)   = 0.97e-12/2.0
3300          xl(iref(p_xyl-1),8:8)   = 0.16e-12/2.0
3301          xl(iref(p_xyl-1),7:7)   = 0.13e-12/2.0
3302          xl(iref(p_xyl-1),6:6)   = 0.13e-12/2.0
3303          xl(iref(p_xyl-1),5:5)   = 0.28e-13/2.0
3304          xl(iref(p_xyl-1),4:4)   = 0.42e-13/2.0
3305          xl(iref(p_xyl-1),3:3)   = 0.12e-12/2.0
3306          xl(iref(p_xyl-1),2:2)   = 0.39e-15/2.0
3307          xl(iref(p_xyl-1),1:1)   = 0.11e-14/2.0
3308       end if
3309       if( p_eth > 1 ) then
3310          xl(iref(p_eth-1),12:16) = 0.60e-9
3311          xl(iref(p_eth-1),11:11) = 0.63e-9
3312          xl(iref(p_eth-1),10:10) = 0.49e-9
3313          xl(iref(p_eth-1),9:9)   = 0.51e-9
3314          xl(iref(p_eth-1),8:8)   = 0.44e-9
3315          xl(iref(p_eth-1),7:7)   = 0.43e-9
3316          xl(iref(p_eth-1),6:6)   = 0.45e-9
3317          xl(iref(p_eth-1),5:5)   = 0.47e-9
3318          xl(iref(p_eth-1),4:4)   = 0.48e-9
3319          xl(iref(p_eth-1),3:3)   = 0.47e-9
3320          xl(iref(p_eth-1),2:2)   = 0.38e-9
3321          xl(iref(p_eth-1),1:1)   = 0.36e-9
3322       end if
3323 #endif
3324
3325 #endif
3326
3327 #if (CASENAME == 1)
3328 ! so2 profile based on mirage 2 output, used for neaqs case, 7-20-05 egc
3329 ! decreased by one magnitude, 27-oct-2005 wig
3330       if( p_so2 > 1 ) then
3331          xl(iref(p_so2-1), 1:2) = 0.035e-10
3332          xl(iref(p_so2-1),   3) = 0.081e-10
3333          xl(iref(p_so2-1), 4:8) = 0.10e-10
3334          xl(iref(p_so2-1),   9) = 0.60e-10
3335          xl(iref(p_so2-1), 10) = 1.1e-10
3336          xl(iref(p_so2-1), 11) = 1.46e-10
3337          xl(iref(p_so2-1), 12) = 1.74e-10
3338          xl(iref(p_so2-1), 13) = 1.94e-10
3339          xl(iref(p_so2-1), 14) = 2.80e-10
3340          xl(iref(p_so2-1), 15:16) = 3.0e-10
3341       end if
3342 #endif
3343
3344       if( p_sulf > 1 ) then
3345          xl(iref(p_sulf-1),:)   = conmin
3346       end if
3347
3348     end SUBROUTINE gasprofile_init_pnnl
3349
3350 #ifdef CHEM_DBG_I
3351 !-----------------------------------------------------------------------
3352 subroutine chem_dbg(i,j,k,dtstep,itimestep,                           &
3353      dz8w,t_phy,p_phy,rho_phy,chem,                                   &
3354      emis_ant,                                                        &
3355      ids,ide, jds,jde, kds,kde,                                       &
3356      ims,ime, jms,jme, kms,kme,                                       &
3357      its,ite, jts,jte, kts,kte,                                       &
3358      kemit,                                                           &
3359      ph_macr,ph_o31d,ph_o33p,ph_no2,ph_no3o2,ph_no3o,ph_hno2,         &
3360      ph_hno3,ph_hno4,ph_h2o2,ph_ch2or,ph_ch2om,ph_ch3cho,             &
3361      ph_ch3coch3,ph_ch3coc2h5,ph_hcocho,ph_ch3cocho,                  &
3362      ph_hcochest,ph_ch3o2h,ph_ch3coo2h,ph_ch3ono2,ph_hcochob,ph_n2o5, &
3363      ph_o2                                                            )
3364
3365   IMPLICIT NONE
3366   INTEGER,      INTENT(IN   ) :: i,j,k,                        &
3367                                  ids,ide, jds,jde, kds,kde,    &
3368                                  ims,ime, jms,jme, kms,kme,    &
3369                                  its,ite, jts,jte, kts,kte,    &
3370                                  kemit
3371   real,         intent(in   ) :: dtstep
3372   integer,      intent(in   ) :: itimestep
3373   REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ),         &
3374        INTENT(INOUT ) :: chem
3375   REAL,  DIMENSION( ims:ime , kms:kme , jms:jme ),    &
3376        INTENT(IN   ) ::  dz8w,t_phy,p_phy,rho_phy
3377   REAL, DIMENSION( ims:ime, kms:kemit, jms:jme,num_emis_ant ),          &
3378        INTENT(IN ) ::                                                   &
3379                 emis_ant
3380   REAL,  DIMENSION( ims:ime , kms:kme , jms:jme ),                      &
3381        INTENT(IN   ), OPTIONAL ::                                       &
3382        ph_macr,ph_o31d,ph_o33p,ph_no2,ph_no3o2,ph_no3o,ph_hno2,         &
3383        ph_hno3,ph_hno4,ph_h2o2,ph_ch2or,ph_ch2om,ph_ch3cho,             &
3384        ph_ch3coch3,ph_ch3coc2h5,ph_hcocho,ph_ch3cocho,                  &
3385        ph_hcochest,ph_ch3o2h,ph_ch3coo2h,ph_ch3ono2,ph_hcochob,ph_n2o5, &
3386        ph_o2
3387
3388   integer :: n
3389   real :: conva,convg
3390
3391   print*,"itimestep =",itimestep
3392
3393   print*,"MET DATA AT (i,k,j):",i,k,j
3394   print*,"t_phy,p_phy,rho_phy=",t_phy(i,k,j),p_phy(i,k,j),rho_phy(i,k,j)
3395
3396   if(dz8w(i,k,j) /= 0.) then
3397      conva = dtstep/(dz8w(i,k,j)*60.)
3398      convg = 4.828e-4/rho_phy(i,k,j)*dtstep/(dz8w(i,k,j)*60.)
3399      print*,"ADJUSTED EMISSIONS (PPM) AT (i,k,j):",i,k,j
3400      print*,"dtstep,dz8w(i,k,j):",dtstep,dz8w(i,k,j)
3401      print*,"e_pm25 i,j:",emis_ant(i,k,j,p_e_pm25i)*conva, &
3402           emis_ant(i,k,j,p_e_pm25j)*conva
3403      print*,"e_ec i,j:",emis_ant(i,k,j,p_e_eci)*conva, &
3404           emis_ant(i,k,j,p_e_ecj)*conva
3405      print*,"e_org i,j:",emis_ant(i,k,j,p_e_orgi)*conva, &
3406           emis_ant(i,k,j,p_e_orgj)*conva
3407      print*,"e_so2:",emis_ant(i,k,j,p_e_so2)*convg
3408      print*,"e_no:",emis_ant(i,k,j,p_e_no)*convg
3409      print*,"e_co:",emis_ant(i,k,j,p_e_co)*convg
3410      print*,"e_eth:",emis_ant(i,k,j,p_e_eth)*convg
3411      print*,"e_hc3:",emis_ant(i,k,j,p_e_hc3)*convg
3412      print*,"e_hc5:",emis_ant(i,k,j,p_e_hc5)*convg
3413      print*,"e_hc8:",emis_ant(i,k,j,p_e_hc8)*convg
3414      print*,"e_xyl:",emis_ant(i,k,j,p_e_xyl)*convg
3415      print*,"e_ol2:",emis_ant(i,k,j,p_e_ol2)*convg
3416      print*,"e_olt:",emis_ant(i,k,j,p_e_olt)*convg
3417      print*,"e_oli:",emis_ant(i,k,j,p_e_oli)*convg
3418      print*,"e_tol:",emis_ant(i,k,j,p_e_tol)*convg
3419      print*,"e_csl:",emis_ant(i,k,j,p_e_csl)*convg
3420      print*,"e_hcho:",emis_ant(i,k,j,p_e_hcho)*convg
3421      print*,"e_ald:",emis_ant(i,k,j,p_e_ald)*convg
3422      print*,"e_ket:",emis_ant(i,k,j,p_e_ket)*convg
3423      print*,"e_ora2:",emis_ant(i,k,j,p_e_ora2)*convg
3424      print*,"e_pm25:",emis_ant(i,k,j,p_e_pm_25)*conva
3425      print*,"e_pm10:",emis_ant(i,k,j,p_e_pm_10)*conva
3426      print*,"e_nh3:",emis_ant(i,k,j,p_e_nh3)*convg
3427      print*,"e_no2:",emis_ant(i,k,j,p_e_no2)*convg
3428      print*,"e_ch3oh:",emis_ant(i,k,j,p_e_ch3oh)*convg
3429      print*,"e_c2h5oh:",emis_ant(i,k,j,p_e_c2h5oh)*convg
3430      print*,"e_iso:",emis_ant(i,k,j,p_e_iso)*convg
3431      print*,"e_so4 f,c:",emis_ant(i,k,j,p_e_so4j)*conva, &
3432           emis_ant(i,k,j,p_e_so4c(i,k,j)*conva
3433      print*,"e_no3 f,c:",emis_ant(i,k,j,p_e_no3j)*conva, &
3434           emis_ant(i,k,j,p_e_no3c(i,k,j)*conva
3435      print*,"e_orgc:",emis_ant(i,k,j,p_e_orgc)*conva
3436      print*,"e_ecc:",emis_ant(i,k,j,p_e_ecc)*conva
3437      print*
3438   else
3439      print*,"dz8w=0 so cannot show adjusted emissions"
3440   end if
3441   print*,"CHEM_DBG PRINT (PPM or ug/m^3) AT (i,k,j):",i,k,j
3442   do n=1,num_chem
3443      print*,n,chem(i,k,j,n)
3444   end do
3445   if( present(ph_macr) ) then
3446      print*,"PHOTOLYSIS DATA:"
3447      print*,"ph_macr:",ph_macr(i,:,j)
3448      print*,"ph_o31d:",ph_o31d(i,:,j)
3449      print*,"ph_o33p:",ph_o33p(i,:,j)
3450      print*,"ph_no2:",ph_no2(i,:,j)
3451      print*,"ph_no3o2:",ph_no3o2(i,:,j)
3452      print*,"ph_no3o:",ph_no3o(i,:,j)
3453      print*,"ph_hno2:",ph_hno2(i,:,j)
3454      print*,"ph_hno3:",ph_hno3(i,:,j)
3455      print*,"ph_hno4:",ph_hno4(i,:,j)
3456      print*,"ph_h2o2:",ph_h2o2(i,:,j)
3457      print*,"ph_ch2or:",ph_ch2or(i,:,j)
3458      print*,"ph_ch2om:",ph_ch2om(i,:,j)
3459      print*,"ph_ch3cho:",ph_ch3cho(i,:,j)
3460      print*,"ph_ch3coch3:",ph_ch3coch3(i,:,j)
3461      print*,"ph_ch3coc2h5:",ph_ch3coc2h5(i,:,j)
3462      print*,"ph_hcocho:",ph_hcocho(i,:,j)
3463      print*,"ph_ch3cocho:",ph_ch3cocho(i,:,j)
3464      print*,"ph_hcochest:",ph_hcochest(i,:,j)
3465      print*,"ph_ch3o2h:",ph_ch3o2h(i,:,j)
3466      print*,"ph_ch3coo2h:",ph_ch3coo2h(i,:,j)
3467      print*,"ph_ch3ono2:",ph_ch3ono2(i,:,j)
3468      print*,"ph_hcochob:",ph_hcochob(i,:,j)
3469      print*,"ph_n2o5:",ph_n2o5(i,:,j)
3470      print*,"ph_o2:",ph_o2(i,:,j)
3471   end if
3472   print*
3473 end subroutine chem_dbg
3474 #endif
3475
3476 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3477 SUBROUTINE med_read_bin_chem_emiss ( grid , config_flags ,intime, itime_max)
3478   ! Driver layer
3479 !  USE module_domain
3480 !  USE module_io_domain
3481 !  USE module_timing
3482   ! Model layer
3483 !  USE module_configure
3484 !  USE module_bc_time_utilities
3485 #ifdef DM_PARALLEL
3486 !  USE module_dm
3487 #endif
3488 !  USE module_date_time
3489 !  USE module_utility
3490
3491
3492    IMPLICIT NONE
3493
3494   ! Arguments
3495    TYPE(domain)                               :: grid
3496    TYPE (grid_config_rec_type) , INTENT(IN)   :: config_flags
3497   !INTEGER , INTENT(IN)                       :: start_step , step , end_step
3498 !  Type (ESMF_Time )                          :: start_time, stop_time, CurrTime
3499 !   TYPE(WRFU_TimeInterval) :: time_interval
3500
3501
3502   ! Local data
3503    LOGICAL, EXTERNAL                      :: wrf_dm_on_monitor
3504    LOGICAL                                :: emiss_opened
3505    INTEGER                                :: intime, itime,itime_max, ierr , open_status , fid
3506    REAL                                   :: time, btime, bfrq
3507    REAL, ALLOCATABLE :: dumc0(:,:,:)
3508    CHARACTER (LEN=256)                    :: message
3509    CHARACTER (LEN=80)                     :: bdyname
3510
3511    CHARACTER (LEN=9 ),DIMENSION(30)       :: ename
3512    INTEGER :: nv,i , j , k, &
3513               ids, ide, jds, jde, kds, kde,    &
3514               ims, ime, jms, jme, kms, kme,    &
3515               ips, ipe, jps, jpe, kps, kpe
3516
3517
3518 #include <wrf_io_flags.h>
3519
3520    write(message, '(A,I9)') 'call read emissions', intime
3521    call wrf_message( TRIM( message ) )
3522
3523    IF(intime == 0 ) THEN
3524      CALL construct_filename1 ( bdyname , '../../run/wrfem12k_00to12z' , grid%id , 2 )
3525
3526      IF (wrf_dm_on_monitor()) THEN
3527         open (91,file=bdyname,form='unformatted')
3528      ENDIf
3529    write(message, '(A,A)') ' OPENED FILE: ',bdyname
3530    call wrf_message( TRIM( message ) )
3531    ENDIF
3532    IF(intime == 12 ) THEN
3533      CALL construct_filename1 ( bdyname , '../../run/wrfem12k_12to24z' , grid%id , 2 )
3534
3535      IF (wrf_dm_on_monitor()) THEN
3536         open (91,file=bdyname,form='unformatted')
3537      ENDIf
3538      write(message, '(A,A)') ' OPENED FILE: ',bdyname
3539      call wrf_message( TRIM( message ) )
3540    ENDIF
3541    CALL wrf_debug( 100 , 'med_read_bin_chem_emiss: calling emissions' )
3542
3543    CALL get_ijk_from_grid (  grid ,                        &
3544                              ids, ide, jds, jde, kds, kde,    &
3545                              ims, ime, jms, jme, kms, kme,    &
3546                              ips, ipe, jps, jpe, kps, kpe    )
3547
3548    ALLOCATE (dumc0(ids:ide-1,kds:grid%kemit,jds:jde-1))
3549
3550      write(message, '(A,6I6)') ' I am reading emissions, dims: =',ids, ide-1, jds, jde-1, kds, grid%kemit
3551      call wrf_message( TRIM( message ) )
3552
3553    IF(intime == 0 .or. intime == 12) then
3554      read(91)nv
3555      read(91)ename
3556      write(message, '(A,I10)') ' Number of emissions: ',nv
3557      call wrf_message( TRIM( message ) )
3558
3559 !    write(message, '(A,30A10)') ' Array names : ',ename
3560 !    call wrf_message( TRIM( message ) )
3561    ENDIF
3562        read(91)itime
3563      write(message, '(A,I8,A,I8)') ' EMISSIONS INPUT FILE TIME PERIOD (GMT): ',itime-1,' TO ',itime
3564      call wrf_message( TRIM( message ) )
3565
3566          read(91)dumc0
3567          grid%emis_ant(ids:ide-1,kds:grid%kemit,jds:jde-1,p_e_so2)=dumc0
3568          read(91)dumc0
3569          grid%emis_ant(ids:ide-1,kds:grid%kemit,jds:jde-1,p_e_no)=dumc0
3570          read(91)dumc0
3571          grid%emis_ant(ids:ide-1,kds:grid%kemit,jds:jde-1,p_e_ald)=dumc0
3572          read(91)dumc0
3573          grid%emis_ant(ids:ide-1,kds:grid%kemit,jds:jde-1,p_e_hcho)=dumc0
3574          read(91)dumc0
3575          grid%emis_ant(ids:ide-1,kds:grid%kemit,jds:jde-1,p_e_ora2)=dumc0
3576          read(91)dumc0
3577          grid%emis_ant(ids:ide-1,kds:grid%kemit,jds:jde-1,p_e_nh3)=dumc0
3578          read(91)dumc0
3579          grid%emis_ant(ids:ide-1,kds:grid%kemit,jds:jde-1,p_e_hc3)=dumc0
3580          read(91)dumc0
3581          grid%emis_ant(ids:ide-1,kds:grid%kemit,jds:jde-1,p_e_hc5)=dumc0
3582          read(91)dumc0
3583          grid%emis_ant(ids:ide-1,kds:grid%kemit,jds:jde-1,p_e_hc8)=dumc0
3584          read(91)dumc0
3585          grid%emis_ant(ids:ide-1,kds:grid%kemit,jds:jde-1,p_e_eth)=dumc0
3586          read(91)dumc0
3587          grid%emis_ant(ids:ide-1,kds:grid%kemit,jds:jde-1,p_e_co)=dumc0
3588          read(91)dumc0
3589          grid%emis_ant(ids:ide-1,kds:grid%kemit,jds:jde-1,p_e_ol2)=dumc0
3590          read(91)dumc0
3591          grid%emis_ant(ids:ide-1,kds:grid%kemit,jds:jde-1,p_e_olt)=dumc0
3592          read(91)dumc0
3593          grid%emis_ant(ids:ide-1,kds:grid%kemit,jds:jde-1,p_e_oli)=dumc0
3594          read(91)dumc0
3595          grid%emis_ant(ids:ide-1,kds:grid%kemit,jds:jde-1,p_e_tol)=dumc0
3596          read(91)dumc0
3597          grid%emis_ant(ids:ide-1,kds:grid%kemit,jds:jde-1,p_e_xyl)=dumc0
3598          read(91)dumc0
3599          grid%emis_ant(ids:ide-1,kds:grid%kemit,jds:jde-1,p_e_ket)=dumc0
3600          read(91)dumc0
3601          grid%emis_ant(ids:ide-1,kds:grid%kemit,jds:jde-1,p_e_csl)=dumc0
3602          read(91)dumc0
3603          grid%emis_ant(ids:ide-1,kds:grid%kemit,jds:jde-1,p_e_iso)=dumc0
3604          read(91)dumc0
3605          grid%emis_ant(ids:ide-1,kds:grid%kemit,jds:jde-1,p_e_pm25i)=dumc0
3606          read(91)dumc0
3607          grid%emis_ant(ids:ide-1,kds:grid%kemit,jds:jde-1,p_e_pm25j)=dumc0
3608          read(91)dumc0
3609          grid%emis_ant(ids:ide-1,kds:grid%kemit,jds:jde-1,p_e_so4i)=dumc0
3610          read(91)dumc0
3611          grid%emis_ant(ids:ide-1,kds:grid%kemit,jds:jde-1,p_e_so4j)=dumc0
3612          read(91)dumc0
3613          grid%emis_ant(ids:ide-1,kds:grid%kemit,jds:jde-1,p_e_no3i)=dumc0
3614          read(91)dumc0
3615          grid%emis_ant(ids:ide-1,kds:grid%kemit,jds:jde-1,p_e_no3j)=dumc0
3616          read(91)dumc0
3617          grid%emis_ant(ids:ide-1,kds:grid%kemit,jds:jde-1,p_e_orgi)=dumc0
3618          read(91)dumc0
3619          grid%emis_ant(ids:ide-1,kds:grid%kemit,jds:jde-1,p_e_orgj)=dumc0
3620          read(91)dumc0
3621          grid%emis_ant(ids:ide-1,kds:grid%kemit,jds:jde-1,p_e_eci)=dumc0
3622          read(91)dumc0
3623          grid%emis_ant(ids:ide-1,kds:grid%kemit,jds:jde-1,p_e_ecj)=dumc0
3624          read(91)dumc0
3625          grid%emis_ant(ids:ide-1,kds:grid%kemit,jds:jde-1,p_e_pm_10)=dumc0
3626
3627     DEALLOCATE ( dumc0 )
3628    RETURN
3629 END SUBROUTINE med_read_bin_chem_emiss
3630
3631 SUBROUTINE med_read_bin_chem_fireemiss ( grid , config_flags ,intime, itime_max)
3632 !   USE module_bc_time_utilities
3633 #ifdef DM_PARALLEL
3634 !   USE module_dm
3635 #endif
3636
3637
3638    IMPLICIT NONE
3639
3640   ! Arguments
3641    TYPE(domain)                               :: grid
3642    TYPE (grid_config_rec_type) , INTENT(IN)   :: config_flags
3643   !INTEGER , INTENT(IN)                       :: start_step , step , end_step
3644 !  Type (ESMF_Time )                          :: start_time, stop_time, CurrTime
3645 !   TYPE(WRFU_TimeInterval) :: time_interval
3646
3647
3648   ! Local data
3649    LOGICAL, EXTERNAL                      :: wrf_dm_on_monitor
3650    LOGICAL                                :: emiss_opened
3651    INTEGER                                :: intime, itime,itime_max, ierr , open_status , fid
3652    REAL                                   :: time, btime, bfrq
3653    REAL, ALLOCATABLE :: dumc0(:,:)
3654    CHARACTER (LEN=256)                    :: message
3655    CHARACTER (LEN=80)                     :: bdyname
3656
3657    CHARACTER (LEN=9 ),DIMENSION(30)       :: ename
3658    INTEGER :: nv,i , j , k, &
3659               ids, ide, jds, jde, kds, kde,    &
3660               ims, ime, jms, jme, kms, kme,    &
3661               ips, ipe, jps, jpe, kps, kpe
3662
3663
3664 #include <wrf_io_flags.h>
3665
3666    write(message, '(A,I9)') 'call read fire emissions', intime
3667    call wrf_message( TRIM( message ) )
3668
3669    IF(intime == 0 ) THEN
3670      CALL construct_filename1 ( bdyname , '../../run/fireemiss' , grid%id , 2 )
3671
3672      IF (wrf_dm_on_monitor()) THEN
3673         open (91,file=bdyname,form='unformatted')
3674      ENDIf
3675    write(message, '(A,A)') ' OPENED FILE: ',bdyname
3676    call wrf_message( TRIM( message ) )
3677    ENDIF
3678    CALL wrf_debug( 100 , 'med_read_bin_chem_fireemiss: calling emissions' )
3679
3680    CALL get_ijk_from_grid (  grid ,                        &
3681                              ids, ide, jds, jde, kds, kde,    &
3682                              ims, ime, jms, jme, kms, kme,    &
3683                              ips, ipe, jps, jpe, kps, kpe    )
3684
3685    ALLOCATE (dumc0(ids:ide-1,jds:jde-1))
3686
3687      write(message, '(A,4I6)') ' I am reading fire emissions, dims: =',ids, ide-1, jds, jde-1
3688      call wrf_message( TRIM( message ) )
3689
3690    IF(intime == 0 .or. intime == 12) then
3691      read(92)nv
3692      read(92)ename
3693      write(message, '(A,I10)') ' Number of emissions: ',nv
3694      call wrf_message( TRIM( message ) )
3695
3696    ENDIF
3697        read(92)itime
3698      write(message, '(A,I8,A,I8)') ' EMISSIONS INPUT FILE TIME PERIOD (GMT): ',itime-1,' TO ',itime
3699      call wrf_message( TRIM( message ) )
3700
3701          read(92)dumc0
3702 !        grid%e_so2(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
3703          read(92)dumc0
3704 !        grid%e_no(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
3705          read(92)dumc0
3706 !        grid%e_ald(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
3707          read(92)dumc0
3708 !        grid%e_hcho(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
3709          read(92)dumc0
3710 !        grid%e_ora2(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
3711          read(92)dumc0
3712 !        grid%e_nh3(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
3713          read(92)dumc0
3714 !        grid%e_hc3(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
3715          read(92)dumc0
3716 !        grid%e_hc5(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
3717          read(92)dumc0
3718 !        grid%e_hc8(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
3719          read(92)dumc0
3720 !        grid%e_eth(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
3721          read(92)dumc0
3722 !        grid%e_co(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
3723          read(92)dumc0
3724 !        grid%e_ol2(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
3725          read(92)dumc0
3726 !        grid%e_olt(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
3727          read(92)dumc0
3728 !        grid%e_oli(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
3729          read(92)dumc0
3730 !        grid%e_tol(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
3731          read(92)dumc0
3732 !        grid%e_xyl(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
3733          read(92)dumc0
3734 !        grid%e_ket(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
3735          read(92)dumc0
3736 !        grid%e_csl(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
3737          read(92)dumc0
3738 !        grid%e_iso(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
3739          read(92)dumc0
3740 !        grid%e_pm25i(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
3741          read(92)dumc0
3742 !        grid%e_pm25j(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
3743          read(92)dumc0
3744 !        grid%e_so4i(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
3745          read(92)dumc0
3746 !        grid%e_so4j(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
3747          read(92)dumc0
3748 !        grid%e_no3i(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
3749          read(92)dumc0
3750 !        grid%e_no3j(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
3751          read(92)dumc0
3752 !        grid%e_orgi(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
3753          read(92)dumc0
3754 !        grid%e_orgj(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
3755          read(92)dumc0
3756 !        grid%e_eci(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
3757          read(92)dumc0
3758 !        grid%e_ecj(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
3759          read(92)dumc0
3760 !        grid%e_pm10(ids:ide-1,kds:grid%kemit,jds:jde-1)=dumc0
3761
3762     DEALLOCATE ( dumc0 )
3763    RETURN
3764 END SUBROUTINE med_read_bin_chem_fireemiss
3765
3766 SUBROUTINE mozcart_lbc_init( chem, num_chem,  &
3767                              ims, ime, jms, jme, kms, kme,    &
3768                              its, ite, jts, jte, kts )
3769
3770     USE module_state_description, only : p_ch4, p_n2o, p_h2
3771
3772 integer, intent(in) :: num_chem
3773 integer, intent(in) :: ims, ime, jms, jme, kms, kme,    &
3774                        its, ite, jts, jte, kts
3775 real, intent(in)    :: chem(ims:ime,kms:kme,jms:jme,num_chem)
3776
3777    integer           :: astat
3778
3779    if( p_ch4 > 1 ) then
3780       allocate( ch4_lbc(its:ite,jts:jte),stat=astat )
3781       if( astat /= 0 ) then
3782          call wrf_error_fatal("mozcart_lbc_init: failed to allocate ch4 lbc")
3783       end if
3784       ch4_lbc(its:ite,jts:jte) = chem(its:ite,kts,jts:jte,p_ch4)
3785    end if
3786    if( p_n2o > 1 ) then
3787       allocate( n2o_lbc(its:ite,jts:jte),stat=astat )
3788       if( astat /= 0 ) then
3789          call wrf_error_fatal("mozcart_lbc_init: failed to allocate n2o lbc")
3790       end if
3791       n2o_lbc(its:ite,jts:jte) = chem(its:ite,kts,jts:jte,p_n2o)
3792    end if
3793    if( p_h2 > 1 ) then
3794       allocate( h2_lbc(its:ite,jts:jte),stat=astat )
3795       if( astat /= 0 ) then
3796          call wrf_error_fatal("mozcart_lbc_init: failed to allocate h2 lbc")
3797       end if
3798    end if
3799    h2_lbc(its:ite,jts:jte) = chem(its:ite,kts,jts:jte,p_h2)
3800
3801 END SUBROUTINE mozcart_lbc_init
3802
3803 SUBROUTINE mozcart_lbc_set( chem, num_chem,  &
3804                             ims, ime, jms, jme, kms, kme,    &
3805                             its, ite, jts, jte, kts )
3806
3807     USE module_state_description, only : p_ch4, p_n2o, p_h2
3808
3809 integer, intent(in) :: num_chem
3810 integer, intent(in) :: ims, ime, jms, jme, kms, kme,    &
3811                        its, ite, jts, jte, kts
3812 real, intent(inout) :: chem(ims:ime,kms:kme,jms:jme,num_chem)
3813
3814    if( p_ch4 > 1 ) then
3815       chem(its:ite,kts,jts:jte,p_ch4) = ch4_lbc(its:ite,jts:jte)
3816    end if
3817    if( p_n2o > 1 ) then
3818       chem(its:ite,kts,jts:jte,p_n2o) = n2o_lbc(its:ite,jts:jte)
3819    end if
3820    if( p_h2 > 1 ) then
3821       chem(its:ite,kts,jts:jte,p_h2) = h2_lbc(its:ite,jts:jte)
3822    end if
3823
3824 END SUBROUTINE mozcart_lbc_set
3825
3826 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3827 END MODULE module_input_chem_data