chem/module_upper_bc_driver.F

   1 ! Description:
   2 !     Routines for setting upper boundary values for :
   3 !     o3, no, no2, hno3, ch4, co, n2o and n2o5
   4 !     From the default pressure level of 50 mb to the top level of the model,
   5 !     the mixing ratio was set to a value from the monthly climatology
   6 !     which varied only in latitude (from ub_vals_moz3.nc file).
   7 !     Between the tropopause level and the default pressure level of 50 mb,
   8 !     the mixing ratio was relaxed with a 10-day relaxation factor.
   9 !     The tropopause level was determined by the primary routine of
  10 !     the subroutine tropopause_twmo and the backup routine of
  11 !     the subroutine tropopause_climate (using data from clim_p_trop.nc file).
  12
  13 ! Author:
  14 !     Adopted from module mo_fstrat of CAM-Chem
  15 !     and put into the current format by Jeff Lee, Feb. 2011
  16
  17 module module_upper_bc_driver
  18
  19       implicit none
  20
  21       private
  22
  23       public  :: upper_bc_init
  24       public  :: upper_bc_driver
  25
  26       save
  27
  28       real,    parameter :: mb2Pa    = 100.
  29       real,    parameter :: vmr2ppmv = 1.e6
  30
  31       real,    parameter :: tau_relax = 864000.   ! 10 days relaxation
  32       real               :: fac_relax
  33       real               :: fixed_ubc_press
  34
  35       integer :: iend
  36       integer :: jend
  37
  38       integer :: ub_month_n      ! # of month   in ubc file
  39       integer :: ub_lev_n        ! # of lev     in ubc file
  40       integer :: ub_lat_n        ! # of lat     in ubc file
  41       integer :: ub_species_n    ! # of species in ubc file
  42       integer :: ub_nchar_n      ! # of nchar   in ubc file
  43       integer :: max_dom         ! wrf domain count
  44
  45       real    , allocatable :: ub_p_m(:)        ! ubc midpoint levels pressure(Pa)
  46       real    , allocatable :: ub_p_e(:)        ! ubc edge     levels pressure(Pa)
  47       real    , allocatable :: ub_vmr(:,:,:,:)  ! values of vmr in ubc file
  48       real    , allocatable :: ub_lat(:)        ! values of lat in ubc file
  49       integer , allocatable :: ub_map(:)        ! species indices for ubc species
  50
  51       integer :: nox_ndx   = -1
  52       integer :: ox_ndx    = -1
  53
  54 !----------------------------------------------------
  55 !     private variables
  56 !----------------------------------------------------
  57
  58       type chem_upper_bc
  59           real, pointer :: vmr(:,:,:,:,:)  ! values of vmr in ubc file
  60           logical       :: is_allocated
  61       end type chem_upper_bc
  62
  63       type(chem_upper_bc), private, allocatable :: upper_bc(:)
  64
  65       contains
  66
  67 !==========================================================================
  68       subroutine upper_bc_driver( id, dtstep, current_date_char,   &
  69                                   chem, p_phy, p8w, tropo_lev,     &
  70                                   ids, ide, jds, jde, kds, kde,    &
  71                                   ims, ime, jms, jme, kms, kme,    &
  72                                   its, ite, jts, jte, kts, kte     )
  73       use module_configure
  74       use module_state_description
  75       use module_model_constants
  76
  77 !----------------------------------------------------
  78 !     input arguments
  79 !----------------------------------------------------
  80       integer, intent(in) :: id,                           &
  81                              ids,ide, jds,jde, kds,kde,    &
  82                              ims,ime, jms,jme, kms,kme,    &
  83                              its,ite, jts,jte, kts,kte
  84
  85       real, dimension(ims:ime, kms:kme, jms:jme ),         &
  86                intent(in) :: p_phy,  & ! p at mid-level (Pa)
  87                              p8w       ! p at interface (Pa)
  88       integer, dimension(ims:ime, jms:jme ),               &
  89                intent(in) :: tropo_lev   ! tropopause level
  90
  91       real,    intent(in) :: dtstep
  92
  93       CHARACTER (LEN=256)       , intent(in) :: current_date_char
  94
  95       real, dimension(ims:ime,kms:kme,jms:jme,num_chem ),  &
  96                intent(inout) :: chem
  97
  98 !----------------------------------------------------
  99 !     local variables
 100 !----------------------------------------------------
 101       integer :: next, last
 102       integer :: kk, km
 103       integer :: kmax
 104       integer :: lev_relax
 105       integer :: astat
 106       integer :: i, j, k, m, n
 107       integer :: ku(kts:kte)
 108       integer :: kl(kts:kte)
 109       integer :: del_p(kts:kte)
 110
 111       real    :: del_time
 112       real    :: vmr_relax
 113       real    :: nox_ubc, xno, xno2, rno
 114       real    :: pint_vals(2)
 115       real, allocatable :: table_ox(:)
 116       real, allocatable :: p8w_temp(:)
 117       real, allocatable :: chem_temp(:)
 118       character(len=132) :: message
 119
 120 !-----------------------------------------------------------------------
 121 ! tile dimensions, needed for nest domains!!
 122 !-------------------------------------------
 123       iend = min(ite,ide-1)
 124       jend = min(jte,jde-1)
 125
 126 !------------------------------------------------------------------------
 127 !...  setup for time interpolation
 128 !------------------------------------------------------------------------
 129       call get_time_interp_factors ( current_date_char, last, next, del_time )
 130
 131 !--------------------------------------------------------
 132 ! ... setup for ox conserving interp
 133 !--------------------------------------------------------
 134       if( ox_ndx > 0 ) then
 135         allocate( table_ox(ub_lev_n-2),stat=astat )
 136         if( astat /= 0 ) then
 137            write(message,*) 'upper_bc_driver: table_ox allocation error = ',astat
 138            CALL wrf_message( trim(message) )
 139            call wrf_abort
 140         end if
 141       end if
 142
 143 j_tile_loop : &
 144    do j = jts, jend
 145 i_tile_loop : &
 146    do i = its, iend
 147
 148       !--------------------------------------------------------
 149       ! setup for the pressure interpolation
 150       !--------------------------------------------------------
 151 lev_loop : &
 152       do k = kte,kts,-1     ! from top to bottom
 153
 154          if( p_phy(i,k,j) <= ub_p_m(1) ) then
 155              kl(k) = 1
 156              ku(k) = 1
 157              del_p(k) = 0.
 158          else if(  p_phy(i,k,j) >= ub_p_m(ub_lev_n) ) then
 159              kl(k) = ub_lev_n
 160              ku(k) = ub_lev_n
 161              del_p(k) = 0.
 162          else
 163              do kk = 2,ub_lev_n    ! from top to bottom of ub_p_m
 164                 if(  p_phy(i,k,j) <= ub_p_m(kk) ) then
 165                    ku(k) = kk
 166                    kl(k) = kk - 1
 167                    del_p(k) = log( p_phy(i,k,j)/ub_p_m(kl(k)) ) &
 168                            / log( ub_p_m(ku(k))/ub_p_m(kl(k)) )
 169                    exit
 170                 end if
 171              end do
 172          end if
 173       end do lev_loop
 174
 175       !--------------------------------------------------------
 176       ! find max level less than fixed_ubc_press (default 50 mb)
 177       ! fixed UB vals from top of model to this level
 178       !--------------------------------------------------------
 179
 180       if( p_phy(i,kte,j) > fixed_ubc_press ) then
 181         kmax = kte + 1
 182       else
 183         do kmax = kte,kts+1,-1                        ! from top to bottom
 184           if( p_phy(i,kmax,j) > fixed_ubc_press ) then
 185             exit
 186           end if
 187         end do
 188         kmax = kmax + 1
 189       endif
 190
 191       !--------------------------------------------------------
 192       ! above fixed_ubc_press (default is 50 mb),
 193       ! set the mixing ratios as fixed values
 194       !--------------------------------------------------------
 195    species_loop : &
 196       do m = 1,ub_species_n
 197
 198    ub_overwrite : &
 199         if( ub_map(m) > 0 .and. kmax <= kte ) then
 200       !-----------------------------------------------------------
 201       ! o3
 202       !-----------------------------------------------------------
 203           if( m == ox_ndx ) then
 204             table_ox(1:ub_lev_n-2) =                             &  ! from top to bottom
 205               upper_bc(id)%vmr(i,j,m,last,2:ub_lev_n-1)          &
 206            + (upper_bc(id)%vmr(i,j,m,next,2:ub_lev_n-1)          &
 207            -  upper_bc(id)%vmr(i,j,m,last,2:ub_lev_n-1))*del_time
 208
 209             km = kte - kmax + 1
 210
 211             allocate ( p8w_temp (km+1) )
 212             allocate ( chem_temp(km) )
 213       !-----------------------------------------------------------
 214       ! upper_bc_rebin requires "top-down" ordering
 215       !-----------------------------------------------------------
 216             p8w_temp (1:km+1) = p8w (i,kte:kmax-1:-1,j)
 217
 218             call upper_bc_rebin( ub_lev_n-2, km, ub_p_e, p8w_temp, &
 219                                  table_ox, chem_temp               )
 220
 221             chem(i,kmax:kte,j,p_o3) = chem_temp(km:1:-1)
 222
 223             deallocate( p8w_temp )
 224             deallocate( chem_temp )
 225
 226             cycle species_loop
 227           end if
 228       !-----------------------------------------------------------
 229       ! others
 230       !-----------------------------------------------------------
 231 lev_loop_a : &
 232           do k = kte,kmax,-1    ! from top to bottom
 233
 234             pint_vals(1) = upper_bc(id)%vmr(i,j,m,last,kl(k))  &
 235                          + (upper_bc(id)%vmr(i,j,m,last,ku(k)) &
 236                          -  upper_bc(id)%vmr(i,j,m,last,kl(k)))*del_p(k)
 237             pint_vals(2) = upper_bc(id)%vmr(i,j,m,next,kl(k))  &
 238                          + (upper_bc(id)%vmr(i,j,m,next,ku(k)) &
 239                          -  upper_bc(id)%vmr(i,j,m,next,kl(k)))*del_p(k)
 240
 241       !-----------------------------------------------------------
 242       ! others: hno3, ch4, co, n2o, n2o5
 243       !-----------------------------------------------------------
 244             if( m /= nox_ndx ) then
 245               chem(i,k,j,ub_map(m)) = pint_vals(1)     &
 246                                     + (pint_vals(2) - pint_vals(1))*del_time
 247       !-----------------------------------------------------------
 248       ! others: no and no2
 249       !-----------------------------------------------------------
 250             else
 251               nox_ubc = pint_vals(1) + del_time * (pint_vals(2) - pint_vals(1))
 252               if( p_no >= param_first_scalar ) then
 253                 xno = chem(i,k,j,p_no)
 254               else
 255                 xno = 0.
 256               end if
 257
 258               if( p_no2 >= param_first_scalar ) then
 259                 xno2 = chem(i,k,j,p_no2)
 260               else
 261                 xno2 = 0.
 262               end if
 263
 264               if( xno > 0. .or. xno2 > 0.) then
 265                 rno = xno / (xno + xno2)
 266               end if
 267
 268               if( p_no >= param_first_scalar ) then
 269                 chem(i,k,j,p_no) = rno * nox_ubc
 270               end if
 271               if( p_no2 >= param_first_scalar ) then
 272                 chem(i,k,j,p_no2) = (1. - rno) * nox_ubc
 273               end if
 274             end if
 275           end do lev_loop_a
 276
 277         end if ub_overwrite
 278       end do species_loop
 279
 280    ub_relax : &
 281       if( tropo_lev(i,j) > 0 .and. tropo_lev(i,j) < kmax ) then
 282       !--------------------------------------------------------
 283       !... relax lower stratosphere to extended ubc
 284       !    check to make sure ubc is not being imposed too low
 285       !    lev_relax = lowest model level (highest pressure)
 286       !               in which to relax to ubc
 287       !--------------------------------------------------------
 288         lev_relax = tropo_lev(i,j)
 289
 290         do while( p_phy(i,lev_relax,j) > ub_p_m(ub_lev_n) ) ! ub_p_m(ub_lev_n)= 440 mb
 291           lev_relax = lev_relax + 1                         ! move level upward
 292         end do
 293
 294         if( lev_relax /= tropo_lev(i,j) ) then
 295           write(message,*) 'upper_bc_driver:Warning,raised ubc: at j,i,=  ',j,i
 296           call wrf_message( trim(message) )
 297           write(message,*) 'from ', tropo_lev(i,j),nint(p_phy(i,tropo_lev(i,j)-1,j)/mb2pa),' mb'
 298           call wrf_message( trim(message) )
 299           write(message,*) 'to   ', lev_relax,     nint(p_phy(i,lev_relax,j)     /mb2pa),' mb'
 300           call wrf_message( trim(message) )
 301         endif
 302
 303    species_loop_a : &
 304         do m = 1,ub_species_n
 305    has_ubc : &
 306           if( ub_map(m) > 0 ) then
 307 lev_loop_b: do k = kmax-1,lev_relax,-1      ! from top to bottom
 308
 309               pint_vals(1) = upper_bc(id)%vmr(i,j,m,last,kl(k))         &
 310                            + (upper_bc(id)%vmr(i,j,m,last,ku(k))        &
 311                            -  upper_bc(id)%vmr(i,j,m,last,kl(k)))*del_p(k)
 312               pint_vals(2) = upper_bc(id)%vmr(i,j,m,next,kl(k))         &
 313                            + (upper_bc(id)%vmr(i,j,m,next,ku(k))        &
 314                            -  upper_bc(id)%vmr(i,j,m,next,kl(k)))*del_p(k)
 315
 316               vmr_relax = pint_vals(1) &
 317                         + (pint_vals(2) - pint_vals(1))*del_time
 318
 319               if( m /= nox_ndx ) then
 320                 chem(i,k,j,ub_map(m)) = chem(i,k,j,ub_map(m)) &
 321                       + (vmr_relax - chem(i,k,j,ub_map(m))) * fac_relax
 322               else
 323                 if( p_no >= param_first_scalar ) then
 324                   xno  = chem(i,k,j,p_no)
 325                 else
 326                   xno  = 0.
 327                 endif
 328
 329                 if( p_no2 >= param_first_scalar ) then
 330                   xno2 = chem(i,k,j,p_no2)
 331                 else
 332                   xno2 = 0.
 333                 endif
 334
 335                 if( xno > 0. .or. xno2 > 0. ) then
 336                   rno  = xno / (xno + xno2)
 337                 endif
 338
 339                 nox_ubc = xno + xno2
 340                 nox_ubc = nox_ubc + (vmr_relax - nox_ubc) * fac_relax
 341
 342                 if( p_no >= param_first_scalar ) then
 343                   chem(i,k,j,p_no)  = rno * nox_ubc
 344                 endif
 345
 346                 if( p_no2 >= param_first_scalar ) then
 347                   chem(i,k,j,p_no2)= (1. - rno) * nox_ubc
 348                 endif
 349               endif
 350             end do lev_loop_b
 351           endif has_ubc
 352         end do species_loop_a
 353
 354       endif ub_relax
 355
 356      end do i_tile_loop
 357    end do j_tile_loop
 358
 359    if( allocated( table_ox ) ) then
 360      deallocate( table_ox )
 361    endif
 362
 363    end subroutine upper_bc_driver
 364
 365 !===========================================================================
 366
 367       subroutine upper_bc_init(id, xlat, dt, config_flags,    &
 368                                ids,ide, jds,jde, kds,kde,     &
 369                                ims,ime, jms,jme, kms,kme,     &
 370                                its,ite, jts,jte, kts,kte      )
 371 !---------------------------------------------------------------------
 372 !       ... new initialization routine for upper boundary condition
 373 !---------------------------------------------------------------------
 374
 375       use module_domain    ! needed for p_co, p_o3....
 376       use module_state_description, only : param_first_scalar
 377       use module_configure,         only : grid_config_rec_type
 378
 379       implicit none
 380
 381 !---------------------------------------------------------------------
 382 !       ... dummy arguments
 383 !---------------------------------------------------------------------
 384       integer, intent(in)  :: id,                          &
 385                               ids,ide, jds,jde, kds,kde,   &
 386                               ims,ime, jms,jme, kms,kme,   &
 387                               its,ite, jts,jte, kts,kte
 388       real, intent(in)     :: dt
 389       real,    intent(in)  :: xlat(ims:ime, jms:jme)       ! unit: degree
 390       type(grid_config_rec_type), intent(in) :: config_flags
 391
 392 !---------------------------------------------------------------------
 393 !       ... local variables
 394 !---------------------------------------------------------------------
 395
 396       integer :: astat
 397       integer :: ncid
 398       integer :: varid
 399       integer :: dimid(4)
 400       integer :: start(4)
 401       integer :: count(4)
 402       integer :: dimid_lev
 403       integer :: dimid_lat
 404       integer :: dimid_month
 405       integer :: dimid_species
 406       integer :: ndims
 407
 408       character(len=132) :: message
 409       character(len=128) :: err_msg
 410       character(len=64)  :: filename
 411       character(len=3)   :: id_num
 412       character(len=25), allocatable :: ub_specname(:)
 413
 414       character(len=80) :: attribute
 415
 416
 417       integer :: lati
 418       integer, allocatable :: lat_interp(:,:)
 419       real,    allocatable :: lat_del(:,:)
 420
 421       integer :: i, j, k, m
 422       integer :: j1, j2
 423
 424
 425       LOGICAL , EXTERNAL      :: wrf_dm_on_monitor
 426 !-------------------------------------------------------------------------
 427
 428 #ifdef NETCDF
 429 include 'netcdf.inc'
 430 #else
 431     call wrf_message( 'upper_bc_init: requires netcdf' )
 432     call wrf_abort
 433 #endif
 434
 435       !-------------------------------------
 436       !... note: grid%XLAT(ide,:)  =0.
 437       !... note: grid%XLONG(:,jde) =0.
 438       !...       => iend = min(ite,ide-1)
 439       !             jend = min(jte,jde-1)
 440       !-------------------------------------
 441
 442       iend = min(ite,ide-1)
 443       jend = min(jte,jde-1)
 444
 445 !---------------------------------------------------------------------
 446 ! ... allocate chem_upper_bc type
 447 !---------------------------------------------------------------------
 448   if( id == 1 .and. .not. allocated(upper_bc) ) then
 449     CALL nl_get_max_dom( 1,max_dom )
 450
 451     allocate( upper_bc(max_dom),stat=astat )
 452     if( astat /= 0 ) then
 453       CALL wrf_message( 'upper_bc_init: failed to allocate  upper_bc' )
 454       CALL wrf_abort
 455     end if
 456
 457     upper_bc(:)%is_allocated = .false.
 458   endif
 459
 460 upper_bc_allocated : &
 461   if( .not. upper_bc(id)%is_allocated ) then
 462 !---------------------------------------------------------------------
 463 !... allocate upper_bc type
 464 !--------------------------------------------------------------------
 465
 466 is_d01 : &
 467     IF( id == 1 ) then
 468 master_proc : &
 469       IF( wrf_dm_on_monitor() ) THEN
 470         write(id_num,'(i3)') 100+id
 471         write(message,*) 'upper_bc_init: intializing domain ' // id_num(2:3)
 472         CALL wrf_message( trim(message) )
 473
 474 !---------------------------------------------------------------------
 475 ! ... open upper_bc netcdf file
 476 !---------------------------------------------------------------------
 477         filename = config_flags%fixed_ubc_inname
 478         if( filename == ' ' ) then
 479           call wrf_message( 'upper_bc_init: input filename not specified in namelist' )
 480           call wrf_abort
 481         endif
 482
 483         err_msg = 'upper_bc_init: failed to open file ' // trim(filename)
 484         call handle_ncerr( nf_open( trim(filename), nf_noclobber, ncid ), trim(err_msg) )
 485         write(message,*) 'upper_bc_init: id, open filename= ',id, filename
 486         CALL wrf_message( trim(message) )
 487 !---------------------------------------------------------------------
 488 ! ... get dimensions
 489 !---------------------------------------------------------------------
 490        err_msg = 'upper_bc_init: failed to get month id'
 491        call handle_ncerr( nf_inq_dimid( ncid, 'month', dimid_month ), trim(err_msg) )
 492        err_msg = 'upper_bc_init: failed to get month'
 493        call handle_ncerr( nf_inq_dimlen( ncid, dimid_month, ub_month_n ), trim(err_msg) )
 494
 495       err_msg = 'upper_bc_init: failed to get lat id'
 496       call handle_ncerr( nf_inq_dimid( ncid, 'lat', dimid_lat ), trim(err_msg) )
 497       err_msg = 'upper_bc_init: failed to get lat'
 498       call handle_ncerr( nf_inq_dimlen( ncid, dimid_lat, ub_lat_n ), trim(err_msg) )
 499
 500       err_msg = 'upper_bc_init: failed to get lev id'
 501       call handle_ncerr( nf_inq_dimid( ncid, 'lev', dimid_lev ), trim(err_msg) )
 502       err_msg = 'upper_bc_init: failed to get lev'
 503       call handle_ncerr( nf_inq_dimlen( ncid, dimid_lev, ub_lev_n ), trim(err_msg) )
 504
 505       err_msg = 'upper_bc_init: failed to get species id'
 506       call handle_ncerr( nf_inq_dimid( ncid, 'species', dimid_species ), trim(err_msg) )
 507       err_msg = 'upper_bc_init: failed to get species'
 508       call handle_ncerr( nf_inq_dimlen( ncid, dimid_species, ub_species_n ), trim(err_msg) )
 509
 510       err_msg = 'upper_bc_init: failed to get nchar id'
 511       call handle_ncerr( nf_inq_dimid( ncid, 'nchar', dimid(1) ), trim(err_msg) )
 512       err_msg = 'upper_bc_init: failed to get nchar'
 513       call handle_ncerr( nf_inq_dimlen( ncid, dimid(1), ub_nchar_n ), trim(err_msg) )
 514
 515    END IF master_proc
 516 !---------------------------------------------------------------------
 517 ! ... bcast the dimensions
 518 !---------------------------------------------------------------------
 519 #ifdef DM_PARALLEL
 520       CALL wrf_dm_bcast_integer ( ub_month_n   , 1 )
 521       CALL wrf_dm_bcast_integer ( ub_lat_n     , 1 )
 522       CALL wrf_dm_bcast_integer ( ub_lev_n     , 1 )
 523       CALL wrf_dm_bcast_integer ( ub_species_n , 1 )
 524       CALL wrf_dm_bcast_integer ( ub_nchar_n   , 1 )
 525 #endif
 526 !---------------------------------------------------------------------
 527 ! ... allocate module arrays
 528 !---------------------------------------------------------------------
 529       allocate( ub_lat(ub_lat_n), stat=astat )
 530       if( astat /= 0 ) then
 531          call wrf_message( 'upper_bc_init: failed to allocate ub_lat' )
 532          call wrf_abort
 533       end if
 534
 535       allocate( ub_vmr(ub_lat_n,ub_species_n,ub_month_n,ub_lev_n), stat=astat )
 536       if( astat /= 0 ) then
 537          call wrf_message( 'upper_bc_init: failed to allocate ub_vmr' )
 538          call wrf_abort
 539       end if
 540
 541       allocate( ub_p_m(ub_lev_n), stat=astat )
 542       if( astat /= 0 ) then
 543          call wrf_message( 'upper_bc_init: failed to allocate ub_p_m' )
 544          call wrf_abort
 545       end if
 546
 547       allocate( ub_p_e(ub_lev_n-1), stat=astat )
 548       if( astat /= 0 ) then
 549          call wrf_message( 'upper_bc_init: failed to allocate ub_p_e' )
 550          call wrf_abort
 551       end if
 552
 553       allocate( ub_map(ub_species_n), stat=astat )
 554       if( astat /= 0 ) then
 555          call wrf_message( 'upper_bc_init: failed to allocate ub_map' )
 556          call wrf_abort
 557       end if
 558
 559 !---------------------------------------------------------------------
 560 ! ... read arrays
 561 !---------------------------------------------------------------------
 562 master_proc_a : &
 563    IF ( wrf_dm_on_monitor() ) THEN
 564
 565 !lev
 566       err_msg = 'upper_bc_init: failed to get lev variable id'
 567       call handle_ncerr( nf_inq_varid( ncid, 'lev', varid ), trim(err_msg) )
 568       err_msg = 'upper_bc_init: failed to read lev variable'
 569       call handle_ncerr( nf_get_var_real( ncid, varid, ub_p_m ), trim(err_msg) )
 570
 571       !-------- check unit
 572       attribute(:) = ' '
 573       astat = nf_get_att_text( ncid, varid, 'units', attribute )
 574       if( astat == nf_noerr )then
 575          if( trim(attribute) == 'mb' .or. trim(attribute) == 'hpa' )then
 576             write(message,*) 'upper_bc_init: units for lev = ',trim(attribute),'... converting to pa'
 577             call wrf_message( trim(message) )
 578             ub_p_m(:) = mb2Pa * ub_p_m(:)
 579          else if( trim(attribute) /= 'pa' .and. trim(attribute) /= 'pa' )then
 580             write(message,*) 'upper_bc_init: unknown units for lev, units=*',trim(attribute),'*'
 581             call wrf_message( trim(message) )
 582             call wrf_abort
 583          end if
 584       else
 585          write(message,*) 'upper_bc_init: warning! units attribute for lev missing, assuming mb'
 586          call wrf_message( trim(message) )
 587          ub_p_m(:) = mb2Pa * ub_p_m(:)
 588       end if
 589
 590 ! ... set values of ub_p_e
 591
 592       ub_p_e(1:ub_lev_n-1) = .5*(ub_p_m(1:ub_lev_n-1) + ub_p_m(2:ub_lev_n))
 593
 594 !---------------------------
 595 !lat
 596       err_msg = 'upper_bc_init: failed to get lat variable id'
 597       call handle_ncerr( nf_inq_varid( ncid, 'lat', varid ), trim(err_msg) )
 598       err_msg = 'upper_bc_init: failed to read lat variable'
 599       call handle_ncerr( nf_get_var_real( ncid, varid, ub_lat ), trim(err_msg) )
 600
 601 !---------------------------
 602 !specname
 603       allocate( ub_specname(ub_species_n), stat=astat )
 604       if( astat /= 0 ) then
 605          call wrf_message( 'upper_bc_init: failed to allocate ub_specname' )
 606          call wrf_abort
 607       end if
 608
 609 species_loop : &
 610       do i = 1,ub_species_n
 611         start(:2) = (/ 1, i /)
 612         count(:2) = (/ ub_nchar_n, 1 /)
 613
 614         ub_specname(i)(:) = ' '
 615
 616         err_msg = 'upper_bc_init: failed to get specname variable id'
 617         call handle_ncerr( nf_inq_varid( ncid, 'specname', varid ), trim(err_msg) )
 618         err_msg = 'upper_bc_init: failed to read ub_specname variable'
 619         call handle_ncerr( nf_get_vara_text( ncid, varid, start(:2), count(:2), ub_specname(i) ), trim(err_msg) )
 620
 621         ub_map(i) = 0
 622
 623         if( trim(ub_specname(i)) ==       'HNO3' .and. p_hno3 >= param_first_scalar ) then
 624           ub_map(i) = p_hno3
 625         else if ( trim(ub_specname(i)) == 'CH4'  .and. p_ch4  >= param_first_scalar ) then
 626           ub_map(i) = p_ch4
 627         else if ( trim(ub_specname(i)) == 'CO'   .and. p_co   >= param_first_scalar ) then
 628           ub_map(i) = p_co
 629         else if ( trim(ub_specname(i)) == 'N2O'  .and. p_n2o  >= param_first_scalar ) then
 630           ub_map(i) = p_n2o
 631         else if ( trim(ub_specname(i)) == 'N2O5' .and. p_n2o5 >= param_first_scalar ) then
 632           ub_map(i) = p_n2o5
 633         else if ( trim(ub_specname(i)) == 'OX'   .and. p_o3   >= param_first_scalar ) then
 634           ub_map(i) = p_o3
 635           ox_ndx    = i
 636         else if ( trim(ub_specname(i)) == 'NOX'  .and. &
 637                 (p_no >= param_first_scalar .or. p_no2 >= param_first_scalar) ) then
 638           ub_map(i) = p_no
 639           nox_ndx   = i
 640         endif
 641
 642         if( ub_map(i) == 0 ) then
 643          write(message,*) 'upper_bc_init: ubc table species ',trim(ub_specname(i)), ' not used'
 644          call wrf_message( trim(message) )
 645         end if
 646       end do species_loop
 647
 648 !--------------------------
 649 !vmr
 650       err_msg = 'upper_bc_init: failed to get vmr variable id'
 651       call handle_ncerr( nf_inq_varid( ncid, 'vmr', varid ), trim(err_msg) )
 652
 653       ! check dimensions
 654
 655       err_msg = 'upper_bc_init: failed to get ndims of vmr variable'
 656       call handle_ncerr( nf_inq_varndims( ncid, varid, ndims ), trim(err_msg) )
 657
 658       if( ndims /= 4 ) then
 659          write(message,*) 'upper_bc_init: error! variable vmr has ndims = ',ndims,', expecting 4'
 660          call wrf_message( trim(message) )
 661          call wrf_abort
 662       end if
 663
 664       err_msg = 'upper_bc_init: failed to get dimid of vmr variable'
 665       call handle_ncerr( nf_inq_vardimid( ncid, varid, dimid ), trim(err_msg) )
 666
 667       if( dimid(1) /= dimid_lat   .or. dimid(2) /= dimid_species .or. &
 668           dimid(3) /= dimid_month .or. dimid(4) /= dimid_lev )then
 669           write(message,*) 'upper_bc_init: error! dimensions in wrong order for variable vmr,'// &
 670                'expecting (lat,species,month,lev)'
 671           call wrf_message( trim(message) )
 672           call wrf_abort
 673       end if
 674
 675 !------------------------------------------------------------------
 676 ! ... read in the ub mixing ratio values
 677 !------------------------------------------------------------------
 678       err_msg = 'upper_bc_init: failed to read vmr variable'
 679       call handle_ncerr( nf_get_var_real( ncid, varid, ub_vmr ), trim(err_msg) )
 680
 681 !---------------------------------------------------------------------
 682 !... close input netcdf file
 683 !---------------------------------------------------------------------
 684       err_msg = 'upper_bc_init: failed to close file ' // trim(filename)
 685       call handle_ncerr( nf_close( ncid ), trim(err_msg) )
 686
 687      deallocate( ub_specname, stat=astat )
 688      if( astat /= 0 ) then
 689         write(message,*) ': failed to deallocate ub_specnmae; astat = ',astat
 690         call wrf_message( trim(message) )
 691         call wrf_abort
 692      end if
 693
 694    ENDIF master_proc_a
 695
 696 !---------------------------------------------------------------------
 697 ! ... bcast the variables
 698 !---------------------------------------------------------------------
 699 #ifdef DM_PARALLEL
 700       CALL wrf_dm_bcast_integer ( nox_ndx   , 1 )
 701       CALL wrf_dm_bcast_integer ( ox_ndx    , 1 )
 702
 703       CALL wrf_dm_bcast_integer ( ub_map, size(ub_map) )
 704       CALL wrf_dm_bcast_real    ( ub_p_m, size(ub_p_m) )
 705       CALL wrf_dm_bcast_real    ( ub_p_e, size(ub_p_e) )
 706       CALL wrf_dm_bcast_real    ( ub_lat, size(ub_lat) )
 707       CALL wrf_dm_bcast_real    ( ub_vmr, size(ub_vmr) )
 708 #endif
 709      fixed_ubc_press = config_flags%fixed_ubc_press*mb2Pa
 710    ENDIF is_d01
 711
 712    upper_bc(id)%is_allocated = .true.
 713
 714    write(message,*) 'upper_bc_init: ub_vmr(1,1,1,:)'
 715    call wrf_message( trim(message) )
 716    do k = 1,ub_lev_n,5
 717      write(message,'(1p,5g15.7)') ub_vmr(1,1,1,k:min(k+4,ub_lev_n))
 718      call wrf_message( trim(message) )
 719    end do
 720
 721 !---------------------------------------------------------------------
 722 ! ... allocate domain specific arrays
 723 !---------------------------------------------------------------------
 724    allocate( lat_del(its:iend,jts:jend), stat=astat )
 725    if( astat /= 0 ) then
 726      call wrf_message( 'upper_bc_init: failed to allocate lat_del' )
 727      call wrf_abort
 728    end if
 729
 730    allocate( lat_interp(its:iend,jts:jend), stat=astat )
 731    if( astat /= 0 ) then
 732      call wrf_message( 'upper_bc_init: failed to allocate lat_interp' )
 733      call wrf_abort
 734    end if
 735
 736 !------------------------------------------------------------------
 737 ! ... regrid  mixing ratio to model latitudes
 738 !------------------------------------------------------------------
 739 !-------------------------------------
 740 !... get lat_del, lat_interp and lati
 741 !-------------------------------------
 742    call get_lat_interp_factors( id, xlat, jend, iend,       &
 743                                 ub_lat_n, ub_lat,           &
 744                                 lat_del, lat_interp, lati,  &
 745                                 ids,ide, jds,jde, kds,kde,  &
 746                                 ims,ime, jms,jme, kms,kme,  &
 747                                 its,ite, jts,jte, kts,kte   )
 748
 749 !---------------------------------------------------------------------
 750 ! ... allocate upper_bc(id) component array
 751 !---------------------------------------------------------------------
 752       allocate( upper_bc(id)%vmr(its:iend,jts:jend,ub_species_n,ub_month_n,ub_lev_n), stat=astat )
 753       if( astat /= 0 ) then
 754          call wrf_message( 'upper_bc_init: failed to allocate upper_bc(id)%vrm' )
 755          call wrf_abort
 756       end if
 757
 758 !--------------------------------------------------------------------
 759 !... regrid from ub_vmr to upper_bc(id)%vmr
 760 !--------------------------------------------------------------------
 761    do j = jts,jend
 762      do i = its,iend
 763        j1 = lat_interp(i,j)
 764        j2 = j1 + lati
 765        upper_bc(id)%vmr(i,j,:,:,:) = ub_vmr(j1,:,:,:) + &
 766                         lat_del(i,j) * (ub_vmr(j2,:,:,:) - ub_vmr(j1,:,:,:))
 767      end do
 768    end do
 769
 770 !----------------------------------------------------
 771 !      transform to ppm vmr
 772 !----------------------------------------------------
 773     upper_bc(id)%vmr(:,:,:,:,:) = upper_bc(id)%vmr(:,:,:,:,:)*vmr2ppmv
 774
 775     write(message,*) 'upper_bc_init: upper_bc(id)%vmr(its+5,jts+5,1,6,:)'
 776     call wrf_message( trim(message) )
 777     do k = 1,ub_lev_n,5
 778       write(message,'(1p,5g15.7)') upper_bc(id)%vmr(its+5,jts+5,1,6,k:min(k+4,ub_lev_n))
 779       call wrf_message( trim(message) )
 780     end do
 781
 782 !--------------------------------------------------------
 783 ! ... set up the relaxation for lower stratosphere
 784 !--------------------------------------------------------
 785 ! ... tau_relax = relaxation timescale (in sec)
 786 !     fac_relax = fractional relaxation towards ubc
 787 !             1 => use ubc
 788 !             0 => ignore ubc, use model concentrations
 789 !--------------------------------------------------------
 790
 791     fac_relax = 1. - exp( -real(dt)/tau_relax )
 792
 793 !--------------------------------------------------------
 794
 795       if( id == max_dom ) then
 796         deallocate( ub_vmr)
 797         deallocate( ub_lat)
 798       endif
 799       deallocate( lat_del)
 800       deallocate( lat_interp)
 801
 802       call wrf_message( ' ' )
 803       write(message,*) 'upper_bc_init: DONE intialized domain ',id
 804       call wrf_message( trim(message) )
 805       call wrf_message( ' ' )
 806
 807   endif upper_bc_allocated
 808
 809   end subroutine upper_bc_init
 810
 811 !===========================================================================
 812
 813       subroutine handle_ncerr( ret, mes )
 814 !---------------------------------------------------------------------
 815 !       ... netcdf error handling routine
 816 !---------------------------------------------------------------------
 817
 818       implicit none
 819
 820 !---------------------------------------------------------------------
 821 !       ... dummy arguments
 822 !---------------------------------------------------------------------
 823       integer, intent(in) :: ret
 824       character(len=*), intent(in) :: mes
 825
 826 #ifdef NETCDF
 827 include 'netcdf.inc'
 828 #endif
 829
 830       if( ret /= nf_noerr ) then
 831          call wrf_message( trim(mes) )
 832          call wrf_message( trim(nf_strerror(ret)) )
 833          call wrf_abort
 834       end if
 835
 836       end subroutine handle_ncerr
 837
 838 !===========================================================================
 839
 840       subroutine get_time_interp_factors ( current_date_char, last, next, del_time)
 841
 842       use module_date_time, only : get_julgmt
 843
 844       implicit none
 845
 846 !-----------------------------------------------------------------------------
 847 ! ... dummy arguments
 848 !-----------------------------------------------------------------------------
 849       CHARACTER (LEN=256),intent(in) :: current_date_char
 850
 851       integer, intent(out) :: next, last
 852       real,    intent(out) :: del_time
 853 !-----------------------------------------------------------------------------
 854 ! ... local variables
 855 !-----------------------------------------------------------------------------
 856
 857       integer, parameter :: day_of_year(12) = (/ 16, 45, 75, 105, 136, 166, 197, &
 858                                                  228, 258, 289, 319, 350        /)
 859       INTEGER :: julyr , julday
 860       REAL    :: gmt
 861       real    :: calday
 862
 863       integer :: m
 864 !------------------------------------------------------------------------------
 865
 866       call get_julgmt ( current_date_char, julyr, julday, gmt )
 867
 868       calday = real(julday) + gmt
 869
 870       if( calday < day_of_year(1) ) then
 871          next = 1
 872          last = 12
 873          del_time = (365. + calday - day_of_year(12)) &
 874                 / (365. + day_of_year(1) - day_of_year(12))
 875       else if( calday >= day_of_year(12) ) then
 876          next = 1
 877          last = 12
 878          del_time = (calday - day_of_year(12)) &
 879                 / (365. + day_of_year(1) - day_of_year(12))
 880       else
 881          do m = 11,1,-1
 882             if( calday >= day_of_year(m) ) then
 883                exit
 884             end if
 885          end do
 886          last = m
 887          next = m + 1
 888          del_time = (calday - day_of_year(m)) / (day_of_year(m+1) - day_of_year(m))
 889       end if
 890
 891       del_time = max( min( 1.,del_time ),0. )
 892
 893       end subroutine get_time_interp_factors
 894
 895 !===========================================================================
 896
 897   subroutine upper_bc_rebin( nsrc, ntrg, src_x, trg_x, src, trg )
 898     !---------------------------------------------------------------
 899     !   ... rebin src to trg
 900     !---------------------------------------------------------------
 901
 902     implicit none
 903
 904     !---------------------------------------------------------------
 905     !   ... dummy arguments
 906     !---------------------------------------------------------------
 907     integer, intent(in)   :: nsrc                  ! dimension source array
 908     integer, intent(in)   :: ntrg                  ! dimension target array
 909     real, intent(in)      :: src_x(nsrc+1)         ! source coordinates
 910     real, intent(in)      :: trg_x(ntrg+1)         ! target coordinates
 911     real, intent(in)      :: src(nsrc)             ! source array
 912     real, intent(out)     :: trg(ntrg)             ! target array
 913
 914     !---------------------------------------------------------------
 915     !   ... local variables
 916     !---------------------------------------------------------------
 917     integer  :: i, l
 918     integer  :: si, si1
 919     integer  :: sil, siu
 920     real :: y
 921     real :: sl, su
 922     real :: tl, tu
 923
 924     do i = 1,ntrg
 925        tl = trg_x(i)
 926        if( tl < src_x(nsrc+1) ) then
 927           do sil = 1,nsrc+1
 928              if( tl <= src_x(sil) ) then
 929                 exit
 930              end if
 931           end do
 932           tu = trg_x(i+1)
 933           do siu = 1,nsrc+1
 934              if( tu <= src_x(siu) ) then
 935                 exit
 936              end if
 937           end do
 938           y   = 0.
 939           sil = max( sil,2 )
 940           siu = min( siu,nsrc+1 )
 941           do si = sil,siu
 942              si1 = si - 1
 943              sl  = max( tl,src_x(si1) )
 944              su  = min( tu,src_x(si) )
 945              y   = y + (su - sl)*src(si1)
 946           end do
 947           trg(i) = y/(trg_x(i+1) - trg_x(i))
 948        else
 949           trg(i) = 0.
 950        end if
 951     end do
 952
 953   end subroutine upper_bc_rebin
 954
 955 !===========================================================================
 956
 957   subroutine get_lat_interp_factors( id, xlat, jend, iend,          &
 958                                      from_lat_n, from_lat,          &
 959                                      lat_del, lat_interp, lati,     &
 960                                      ids, ide, jds, jde, kds, kde,  &
 961                                      ims, ime, jms, jme, kms, kme,  &
 962                                      its, ite, jts, jte, kts, kte   )
 963
 964 !---------------------------------------------------------------------
 965 !       ... Determine indicies and deltas for transform
 966 !           Note : it is assumed that the latitude and longitude
 967 !                  arrays are monotonic
 968 !--------------------------------------------------------------------
 969
 970       implicit none
 971
 972 !---------------------------------------------------------------------
 973 !       ... Dummy args
 974 !---------------------------------------------------------------------
 975       integer, intent(in) :: id,                           &
 976                              ids,ide, jds,jde, kds,kde,    &
 977                              ims,ime, jms,jme, kms,kme,    &
 978                              its,ite, jts,jte, kts,kte
 979
 980       real,    intent(in)  :: xlat(ims:ime, jms:jme)       ! unit: degree
 981       integer, intent(in)  :: jend
 982       integer, intent(in)  :: iend
 983       integer, intent(in)  :: from_lat_n
 984       real ,   intent(in)  :: from_lat(from_lat_n)
 985
 986       integer, intent(out) :: lati
 987       integer, intent(out) :: lat_interp(its:iend,jts:jend)
 988       real,    intent(out) :: lat_del(its:iend,jts:jend)
 989
 990 !---------------------------------------------------------------------
 991 !       ... Local variables
 992 !---------------------------------------------------------------------
 993       integer :: astat
 994
 995       real    :: to_lat(jts:jend)       ! values of wrf lat, 1d
 996
 997       integer :: to_lon_n
 998       real    :: countx
 999
1000       real    :: target_lat
1001       integer :: latl, latu
1002       real    :: max_lat, min_lat
1003       logical :: from_lats_mono_pos
1004
1005       integer :: i, j
1006       integer :: m, n
1007
1008 !---------------------------------------------------------------------
1009 !... set "module" variable, converter(id)%lati
1010 !---------------------------------------------------------------------
1011
1012      max_lat = from_lat(from_lat_n)
1013      min_lat = from_lat(1)
1014
1015      if( from_lat(from_lat_n) >=  from_lat(1)) then  ! increasing lat
1016         latl = 1
1017         latu = from_lat_n
1018         lati = 1
1019         from_lats_mono_pos = .true.
1020      else                                            ! decreasing lat
1021         latl = from_lat_n
1022         latu = 1
1023         lati = -1
1024         from_lats_mono_pos = .false.
1025      end if
1026
1027 !---------------------------------------------------------------------
1028 ! ... Set interpolation latitude indicies and del_lat
1029 !---------------------------------------------------------------------
1030       do j = jts,jend
1031       do i = its,iend
1032          target_lat = xlat(i,j)
1033
1034          if( target_lat <= min_lat ) then
1035             lat_del(i,j) = 0.
1036             lat_interp(i,j)  = latl
1037          else if( target_lat >= max_lat ) then
1038             lat_del(i,j) = 1.
1039             if( from_lats_mono_pos ) then
1040                lat_interp(i,j) = latu - 1
1041             else
1042                lat_interp(i,j) = latu + 1
1043             end if
1044          else
1045             do m = latl,latu,lati
1046                if( target_lat < from_lat(m) ) then
1047                   n = m - lati
1048                   lat_interp(i,j) = min( from_lat_n,max( 1,n ) )
1049                   lat_del(i,j)    = &
1050                       (target_lat - from_lat(n))/(from_lat(m) - from_lat(n))
1051                   exit
1052                end if
1053             end do
1054          end if
1055       end do
1056       end do
1057
1058       end subroutine get_lat_interp_factors
1059
1060 end module module_upper_bc_driver
1061