chem/module_cam_mam_initmixrats.F

   1 ! module_cam_mam_initmixrats.F
   2 ! created by r.c.easter, june 2010
   3 !
   4 !--------------------------------------------------------------
   5 #include "MODAL_AERO_CPP_DEFINES.h"
   6
   7       module module_cam_mam_initmixrats
   8
   9       private
  10       public :: bdy_chem_value_cam_mam
  11       public :: cam_mam_initmixrats
  12
  13       contains
  14
  15
  16 !==============================================================
  17         subroutine cam_mam_initmixrats(         &
  18                 iflagaa, numgas, config_flags,  &
  19                 chem, convfac, alt, z_at_w, g,  &
  20                 ids,ide, jds,jde, kds,kde,      &
  21                 ims,ime, jms,jme, kms,kme,      &
  22                 its,ite, jts,jte, kts,kte       )
  23
  24 !
  25 !   initializes the species and number mixing ratios for each section
  26 !
  27 !   this top level routine simply calls other routines depending on value
  28 !       of config_flags%aer_ic_opt
  29 !
  30           !use module_state_description !Balwinder.Singh@gmail.com- Commented out to avoid multiple 'use' declarations of this module
  31         use module_configure, only:  grid_config_rec_type, num_chem, &
  32                 p_ncl_a1, p_ncl_a2, p_nh4_a1, p_nh4_a2, &
  33                 p_so4_a1, p_so4_a2, p_sulf
  34         use module_state_description, only:  num_chem, param_first_scalar,   &
  35                 aer_ic_default, aer_ic_pnnl
  36         use module_data_cam_mam_asect, only:  ai_phase, dens_aer, &
  37                 massptr_aer, ncomp_aer, nsize_aer, ntype_aer, numptr_aer, &
  38                 volumcen_sect
  39         use module_cam_mam_init, only:  pom_icbc_1p4_factor, so4_icbc_1p2_factor
  40
  41         use module_data_sorgam, only:  conmin
  42
  43         implicit none
  44
  45
  46 !   subr arguments
  47         type(grid_config_rec_type), intent(in) :: config_flags
  48
  49         integer, intent(in) ::   &
  50                 iflagaa, numgas,   &
  51                 ids, ide, jds, jde, kds, kde,   &
  52                 ims, ime, jms, jme, kms, kme,   &
  53                 its, ite, jts, jte, kts, kte
  54
  55         real, intent(inout),   &
  56                 dimension( ims:ime, kms:kme, jms:jme, 1:num_chem ) :: &
  57                 chem
  58
  59         real, intent(in),      &
  60                 dimension( ims:ime, kms:kme, jms:jme ) :: &
  61                 alt, convfac, z_at_w
  62 !       convfac = air molar density (mol/m3)
  63         real, intent(in) :: g
  64
  65 !   local variables
  66         integer :: i, ic, iphase, isize, itype
  67         integer :: j
  68         integer :: k
  69         integer :: l, ll
  70
  71         real, parameter :: splitfac = 0.98   ! fraction of initial so4 aerosol in accum mode
  72         real, parameter :: so4vaptoaer = 0.999   ! fraction of initial so4 that is aerosol
  73         real :: tmpa
  74         real :: zz
  75
  76
  77         do ic = numgas+1, num_chem
  78            chem(its:ite,kts:kte,jts:jte,ic) = conmin
  79         end do
  80
  81         do j = jts, jte
  82         do k = kts, kte-1
  83         do i = its, ite
  84
  85         if (config_flags%aer_ic_opt == aer_ic_default) then
  86 ! "default" ic values mimicing subr. aerosols_sorgam_init
  87            chem(i,k,j,p_so4_a1)=chem(i,k,j,p_sulf)*convfac(i,k,j)* &
  88               96.0*splitfac*so4vaptoaer * so4_icbc_1p2_factor
  89            chem(i,k,j,p_so4_a2)=chem(i,k,j,p_sulf)*convfac(i,k,j)* &
  90               96.0*(1.-splitfac)*so4vaptoaer * so4_icbc_1p2_factor
  91            chem(i,k,j,p_sulf)=chem(i,k,j,p_sulf)*(1.-so4vaptoaer)
  92 #if ( defined MODAL_AERO_7MODE )
  93            chem(i,k,j,p_nh4_a1) = 10.e-05
  94            chem(i,k,j,p_nh4_a2) = 10.e-05
  95 #endif
  96            chem(i,k,j,p_ncl_a1) = 20.e-05
  97            chem(i,k,j,p_ncl_a2) = 20.e-05
  98
  99         else if (config_flags%aer_ic_opt == aer_ic_pnnl) then
 100            zz = (z_at_w(i,k,j)+z_at_w(i,k+1,j))*0.5
 101            call cam_mam_init_aer_ic_pnnl(   &
 102                 chem, zz, i,k,j, ims,ime,jms,jme,kms,kme )
 103
 104         else
 105
 106            call wrf_error_fatal(   &
 107                 "cam_mam_initmixrats: bad value for aer_ic_opt" )
 108         end if
 109
 110 ! calculate aerosol number from aerosol volume and prescribed initial size
 111         iphase = ai_phase
 112         do itype = 1, ntype_aer
 113         do isize = 1, nsize_aer(itype)
 114            tmpa = 0.0
 115            do ll = 1, ncomp_aer(itype)
 116               l = massptr_aer(ll,isize,itype,iphase)
 117               if ((l >= param_first_scalar) .and. (l <= num_chem)) &
 118                  tmpa = tmpa + chem(i,k,j,l)/dens_aer(ll,itype)
 119            end do ! ll
 120            tmpa = tmpa*1.0e-6   ! chem is ug/kg but want (g/kg)/dens_aer = cm3/kg
 121            l = numptr_aer(isize,itype,iphase)
 122            chem(i,k,j,l) = tmpa/volumcen_sect(isize,itype)
 123         end do ! isize
 124         end do ! itype
 125
 126         end do ! i
 127         end do ! k
 128         end do ! j
 129
 130
 131 ! Aerosol species are returned from above in concentration units. Convert
 132 ! them to mixing ratio for use in advection, etc.
 133         do ic = numgas+1, num_chem
 134         do j = jts,jte
 135         do k = kts,kte-1
 136         do i = its,ite
 137             chem(i,k,j,ic) = chem(i,k,j,ic)*alt(i,k,j)
 138         end do ! i
 139         end do ! k
 140         end do ! j
 141         end do ! ic
 142
 143 ! Fill the top z-staggered location to prevent trouble during advection.
 144         do ic = numgas+1, num_chem
 145         do j = jts,jte
 146         do i = its,ite
 147             chem(i,kte,j,ic) = chem(i,kte-1,j,ic)
 148         end do ! i
 149         end do ! j
 150         end do ! ic
 151
 152
 153         return
 154         end subroutine cam_mam_initmixrats
 155
 156
 157 !==============================================================
 158 !
 159 !   subroutine to initialize aerosol values using the
 160 !   aer_ic_opt == aer_ic_pnnl option.
 161 !
 162 !   2010-11-17 rce - adapted from subr. sorgam_init_aer_ic_pnnl
 163 !
 164 !   called by: cam_mam_initmixrats
 165 !
 166       subroutine cam_mam_init_aer_ic_pnnl(                  &
 167           chem, z, i,k,j, ims,ime, jms,jme, kms,kme )
 168
 169       use module_configure, only:  grid_config_rec_type, num_chem, &
 170           p_bc_a1, p_bc_a3, p_dst_a1, p_dst_a3, p_dst_a5, p_dst_a7, &
 171           p_ncl_a1, p_ncl_a2, p_ncl_a3, p_ncl_a4, p_ncl_a6, &
 172           p_nh4_a1, p_nh4_a2, p_pom_a1, p_pom_a3, &
 173           p_soa_a1, p_soa_a2, p_so4_a1, p_so4_a2, p_sulf
 174       use module_data_sorgam, only:  conmin
 175       use module_cam_mam_init, only:  pom_icbc_1p4_factor, so4_icbc_1p2_factor
 176
 177       implicit none
 178
 179       integer,intent(in   ) :: i,k,j,                           &
 180                                ims,ime, jms,jme, kms,kme
 181       real,  dimension( ims:ime , kms:kme , jms:jme, num_chem ),&
 182            intent(inout   ) :: chem
 183
 184       real, intent(in     ) :: z
 185       real :: mult
 186
 187 !
 188 ! determine height multiplier...
 189 ! this should mimic the calculation in sorgam_set_aer_bc_pnnl,
 190 ! mosaic_init_wrf_mixrats_opt2, and bdy_chem_value_mosaic
 191 !
 192 !jdf comment these values and have another profile consistent with mosaic
 193         if( z <= 500. ) then
 194            mult = 1.0
 195         elseif( z > 500. &
 196              .and. z <= 1000. ) then
 197            mult = 1.0 - 0.001074*(z-500.)
 198         elseif( z > 1000. &
 199              .and. z <= 5000. ) then
 200            mult = 0.463 - 0.000111*(z-1000.)
 201         else
 202            mult = 0.019
 203         end if
 204
 205 ! these should match what is in sorgam_set_aer_bc_pnnl.
 206 ! values as of 2-dec-2004:
 207 !jdf comment these values and have another profile consistent with mosaic
 208       chem(i,k,j,p_sulf)     = mult*conmin
 209
 210 !     chem(i,k,j,p_so4aj)    = mult*0.300*0.97
 211 !     chem(i,k,j,p_so4ai)    = mult*0.300*0.03
 212       chem(i,k,j,p_so4_a1)   = mult*0.300*0.97 * so4_icbc_1p2_factor
 213       chem(i,k,j,p_so4_a2)   = mult*0.300*0.03 * so4_icbc_1p2_factor
 214
 215 !     chem(i,k,j,p_nh4aj)    = mult*0.094*0.97
 216 !     chem(i,k,j,p_nh4ai)    = mult*0.094*0.03
 217 #if ( defined MODAL_AERO_7MODE )
 218       chem(i,k,j,p_nh4_a1)   = mult*0.094*0.97
 219       chem(i,k,j,p_nh4_a2)   = mult*0.094*0.03
 220 #endif
 221
 222 !     chem(i,k,j,p_no3aj)    = mult*0.001*0.97
 223 !     chem(i,k,j,p_no3ai)    = mult*0.001*0.03
 224
 225 !     chem(i,k,j,p_naaj)     = 10.e-05
 226 !     chem(i,k,j,p_naai)     = 10.e-05
 227 !     chem(i,k,j,p_claj)     = 10.e-05
 228 !     chem(i,k,j,p_clai)     = 10.e-05
 229       chem(i,k,j,p_ncl_a2)   = 20.e-05
 230 #if ( defined MODAL_AERO_3MODE )
 231       chem(i,k,j,p_ncl_a1)   = 20.e-05
 232 #elif ( defined MODAL_AERO_7MODE )
 233       chem(i,k,j,p_ncl_a1)   = 20.e-05*0.2
 234       chem(i,k,j,p_ncl_a4)   = 20.e-05*0.8
 235 #endif
 236
 237 !     chem(i,k,j,p_ecj)      = mult*0.013*0.97
 238 !     chem(i,k,j,p_eci)      = mult*0.013*0.03
 239 #if ( defined MODAL_AERO_3MODE )
 240       chem(i,k,j,p_bc_a1)    = mult*0.013
 241 #elif ( defined MODAL_AERO_7MODE )
 242       chem(i,k,j,p_bc_a1)    = mult*0.013*0.99
 243       chem(i,k,j,p_bc_a3)    = mult*0.013*0.01
 244 #endif
 245
 246 !     chem(i,k,j,p_p25j)     = mult*4.500*0.97
 247 !     chem(i,k,j,p_p25i)     = mult*4.500*0.03
 248 !     chem(i,k,j,p_antha)    = mult*4.500/2.0
 249 #if ( defined MODAL_AERO_3MODE )
 250       chem(i,k,j,p_dst_a1)   = mult*4.500
 251       chem(i,k,j,p_dst_a3)   = mult*4.500/2.0
 252 #elif ( defined MODAL_AERO_7MODE )
 253       chem(i,k,j,p_dst_a5)   = mult*4.500
 254       chem(i,k,j,p_dst_a7)   = mult*4.500/2.0
 255 #endif
 256
 257 !     chem(i,k,j,p_orgpaj)   = mult*0.088*0.97
 258 !     chem(i,k,j,p_orgpai)   = mult*0.088*0.03
 259 #if ( defined MODAL_AERO_3MODE )
 260       chem(i,k,j,p_pom_a1)   = mult*0.088 * pom_icbc_1p4_factor
 261 #elif ( defined MODAL_AERO_7MODE )
 262       chem(i,k,j,p_pom_a1)   = mult*0.088*0.99 * pom_icbc_1p4_factor
 263       chem(i,k,j,p_pom_a3)   = mult*0.088*0.01 * pom_icbc_1p4_factor
 264 #endif
 265
 266 !     chem(i,k,j,p_orgaro1j) = conmin
 267 !     chem(i,k,j,p_orgaro1i) = conmin
 268 !     chem(i,k,j,p_orgaro2j) = conmin
 269 !     chem(i,k,j,p_orgaro2i) = conmin
 270 !     chem(i,k,j,p_orgalk1j) = conmin
 271 !     chem(i,k,j,p_orgalk1i) = conmin
 272 !     chem(i,k,j,p_orgole1j) = conmin
 273 !     chem(i,k,j,p_orgole1i) = conmin
 274 !     chem(i,k,j,p_orgba1j)  = conmin
 275 !     chem(i,k,j,p_orgba1i)  = conmin
 276 !     chem(i,k,j,p_orgba2j)  = conmin
 277 !     chem(i,k,j,p_orgba2i)  = conmin
 278 !     chem(i,k,j,p_orgba3j)  = conmin
 279 !     chem(i,k,j,p_orgba3i)  = conmin
 280 !     chem(i,k,j,p_orgba4j)  = conmin
 281 !     chem(i,k,j,p_orgba4i)  = conmin
 282       chem(i,k,j,p_soa_a1)   = conmin
 283       chem(i,k,j,p_soa_a2)   = conmin
 284
 285 !     chem(i,k,j,p_seas)     = mult*1.75
 286 #if ( defined MODAL_AERO_3MODE )
 287       chem(i,k,j,p_ncl_a3)   = mult*1.75
 288 #elif ( defined MODAL_AERO_7MODE )
 289       chem(i,k,j,p_ncl_a6)   = mult*1.75
 290 #endif
 291
 292
 293       end subroutine cam_mam_init_aer_ic_pnnl
 294
 295
 296 !==============================================================
 297 !
 298 !   subroutine to set aerosol inflow boundary values
 299 !
 300 !   2011-01-11 rce - adapted from subr. bdy_chem_value_sorgam
 301 !
 302 !   called by: flow_dep_bdy_chem
 303 !
 304       subroutine bdy_chem_value_cam_mam( chem, z, nch, config_flags, &
 305                                          alt, convfac, g )
 306
 307       use module_configure, only:  grid_config_rec_type
 308       use module_state_description, only:  aer_bc_default, aer_bc_pnnl
 309       use module_data_sorgam, only: conmin
 310
 311       implicit none
 312
 313       integer, intent(in)   :: nch        ! index number of chemical species
 314       real, intent(out) :: chem
 315       real, intent(in)  :: z          ! 3d height array
 316       real, intent(in)  :: alt, convfac, g
 317       type (grid_config_rec_type), intent(in) :: config_flags
 318
 319 !
 320 ! method for bc calculation is determined by aer_bc_opt
 321 !
 322       if (config_flags%aer_bc_opt == aer_bc_pnnl) then
 323          call cam_mam_set_aer_bc_pnnl( chem, z, nch, config_flags )
 324          return
 325       else if (config_flags%aer_bc_opt == aer_bc_default) then
 326          chem = conmin
 327          return
 328       else
 329          call wrf_error_fatal( &
 330             "bdy_chem_value_cam_mam -- unable to parse aer_bc_opt" )
 331       end if
 332
 333       end subroutine bdy_chem_value_cam_mam
 334
 335
 336 !==============================================================
 337 !
 338 !   subroutine to set aerosol inflow boundary values
 339 !      when aer_ic_opt == aer_ic_pnnl option.
 340 !
 341 !   2011-01-11 rce - adapted from subr. cam_mam_set_aer_bc_pnnl
 342 !
 343 !   called by bdy_chem_value_cam_mam
 344 !
 345       subroutine cam_mam_set_aer_bc_pnnl( chem, z, nch, config_flags )
 346
 347       use module_configure, only:  grid_config_rec_type, num_chem, &
 348           p_bc_a1, p_dst_a1, p_dst_a3, p_dst_a5, p_dst_a7, &
 349           p_ncl_a3, p_ncl_a6, p_nh4_a1, p_nh4_a2, &
 350           p_num_a3, p_num_a7, p_pom_a1, p_so4_a1, p_so4_a2
 351
 352       use module_state_description, only:  param_first_scalar
 353
 354       use module_data_sorgam, only:  conmin
 355
 356       use module_data_cam_mam_asect, only:  ai_phase, dens_aer, massptr_aer, &
 357           ncomp_aer, nsize_aer, ntype_aer, numptr_aer, volumcen_sect
 358
 359       use module_cam_mam_init, only:  pom_icbc_1p4_factor, so4_icbc_1p2_factor
 360
 361       implicit none
 362
 363       integer,intent(in) :: nch
 364       real,intent(in   ) :: z
 365       real,intent(inout) :: chem
 366       type (grid_config_rec_type), intent (in) :: config_flags
 367
 368 ! local variables
 369       integer :: inumber, isize, itype
 370       integer :: l, ll
 371       real :: bv_so4ai, bv_so4aj,         &
 372               bv_nh4ai, bv_nh4aj,         &
 373               bv_no3ai, bv_no3aj,         &
 374               bv_eci,   bv_ecj,           &
 375               bv_p25i,  bv_p25j,          &
 376               bv_orgpai,bv_orgpaj,        &
 377               bv_antha, bv_seas, bv_soila
 378       real :: mult
 379       real :: tmpchem(num_chem)
 380       real :: tmpvol
 381
 382       character(len=160) :: msg
 383
 384
 385       chem = conmin
 386 ! check that nch is an aerosol species
 387 #if ( defined MODAL_AERO_3MODE )
 388       if ((nch < p_so4_a1) .or. (nch > p_num_a3)) return
 389 #elif ( defined MODAL_AERO_7MODE )
 390       if ((nch < p_so4_a1) .or. (nch > p_num_a7)) return
 391 #else
 392       return
 393 #endif
 394
 395
 396 !
 397 ! determine height multiplier...
 398 ! this should mimic the calculation in sorgam_init_aer_ic_pnnl,
 399 ! mosaic_init_wrf_mixrats_opt2, and bdy_chem_value_mosaic
 400 !
 401 ! updated aerosol profile multiplier 1-apr-2005:
 402 !    height(m)     multiplier
 403 !    ---------     ----------
 404 !    <=2000        1.0
 405 !    2000<z<3000   linear transition zone to 0.25
 406 !    3000<z<5000   linear transision zone to 0.125
 407 !    >=5000        0.125
 408 !
 409 ! which translates to:
 410 !    2000<z<3000   mult = 1.00 + (z-2000.)*(0.25-1.0)/(3000.-2000.)
 411 !    3000<z<5000   mult = 0.25 + (z-3000.)*(0.125-0.25)/(5000.-3000.)
 412        if( z <= 500. ) then
 413           mult = 1.0
 414        elseif( z > 500. &
 415             .and. z <= 1000. ) then
 416           mult = 1.0 - 0.001074*(z-500.)
 417        elseif( z > 1000. &
 418             .and. z <= 5000. ) then
 419           mult = 0.463 - 0.000111*(z-1000.)
 420        else
 421           mult = 0.019
 422        end if
 423
 424 ! these should match what is in sorgam_init_aer_ic_pnnl.
 425       bv_so4aj = mult*0.300*0.97 * so4_icbc_1p2_factor
 426       bv_so4ai = mult*0.300*0.03 * so4_icbc_1p2_factor
 427       bv_nh4aj = mult*0.094*0.97
 428       bv_nh4ai = mult*0.094*0.03
 429       bv_no3aj = mult*0.001*0.97
 430       bv_no3ai = mult*0.001*0.03
 431       bv_ecj   = mult*0.013*0.97
 432       bv_eci   = mult*0.013*0.03
 433       bv_p25j  = mult*4.500*0.97
 434       bv_p25i  = mult*4.500*0.03
 435       bv_antha = mult*4.500/2.0
 436       bv_orgpaj = mult*0.088*0.97
 437       bv_orgpai = mult*0.088*0.03
 438       bv_seas   = mult*1.75
 439       bv_soila  = conmin
 440
 441
 442 ! set cam_mam mass mixing ratios
 443 ! (most of them are conmin ~= 0)
 444       tmpchem(:) = conmin
 445
 446       tmpchem(p_so4_a1) = bv_so4aj
 447       tmpchem(p_so4_a2) = bv_so4ai
 448
 449 #if ( defined MODAL_AERO_7MODE )
 450       tmpchem(p_nh4_a1) = bv_nh4aj
 451       tmpchem(p_nh4_a2) = bv_nh4ai
 452 #endif
 453
 454       tmpchem(p_pom_a1) = (bv_orgpai + bv_orgpaj) * pom_icbc_1p4_factor
 455
 456       tmpchem(p_bc_a1 ) = bv_eci + bv_ecj
 457
 458 #if ( defined MODAL_AERO_3MODE )
 459       tmpchem(p_dst_a1) = bv_p25i + bv_p25j
 460       tmpchem(p_dst_a3) = bv_antha + bv_soila
 461 #elif ( defined MODAL_AERO_7MODE )
 462       tmpchem(p_dst_a5) = bv_p25i + bv_p25j
 463       tmpchem(p_dst_a7) = bv_antha + bv_soila
 464 #endif
 465
 466 #if ( defined MODAL_AERO_3MODE )
 467       tmpchem(p_ncl_a3) = bv_seas
 468 #elif ( defined MODAL_AERO_7MODE )
 469       tmpchem(p_ncl_a6) = bv_seas
 470 #endif
 471
 472
 473 ! now set the chem value for species=nch
 474
 475 ! first check if species is an aerosol number by comparing nch to numptr_aer(...)
 476       inumber = 0
 477 itype_loop01: &
 478       do itype = 1, ntype_aer
 479       do isize = 1, nsize_aer(itype)
 480          if (nch == numptr_aer(isize,itype,ai_phase)) then
 481 ! calculate volume mixratio
 482             tmpvol = 0.0
 483             do ll = 1, ncomp_aer(itype)
 484                l = massptr_aer(ll,isize,itype,ai_phase)
 485                tmpvol = tmpchem(l)/dens_aer(ll,itype)
 486             end do
 487 ! calculate number mixratio from volume using default 1-particle volume
 488 ! number mixratio is particles/kg
 489 !    the 1.0e-6 factor is because the tmpchem values are ug/kg,
 490 !    dens_aer are g/cm3, and volumcen_sect are cm3/particle
 491             chem = tmpvol*1.0e-6/volumcen_sect(isize,itype)
 492             inumber = 1
 493             exit itype_loop01
 494          end if
 495       end do
 496       end do itype_loop01
 497
 498       if (inumber <= 0) then
 499 ! species must be an aerosol mass
 500          chem = tmpchem(nch)
 501       end if
 502
 503
 504       return
 505       end subroutine cam_mam_set_aer_bc_pnnl
 506
 507
 508
 509 !==============================================================
 510       end module module_cam_mam_initmixrats
 511