! Description:
! Routines for setting upper boundary values for :
! o3, no, no2, hno3, ch4, co, n2o and n2o5
! From the default pressure level of 50 mb to the top level of the model,
! the mixing ratio was set to a value from the monthly climatology
! which varied only in latitude (from ub_vals_moz3.nc file).
! Between the tropopause level and the default pressure level of 50 mb,
! the mixing ratio was relaxed with a 10-day relaxation factor.
! The tropopause level was determined by the primary routine of
! the subroutine tropopause_twmo and the backup routine of
! the subroutine tropopause_climate (using data from clim_p_trop.nc file).
! Author:
! Adopted from module mo_fstrat of CAM-Chem
! and put into the current format by Jeff Lee, Feb. 2011
module module_upper_bc_driver
implicit none
private
public :: upper_bc_init
public :: upper_bc_driver
save
real, parameter :: mb2Pa = 100.
real, parameter :: vmr2ppmv = 1.e6
real, parameter :: tau_relax = 864000. ! 10 days relaxation
real :: fac_relax
real :: fixed_ubc_press
integer :: iend
integer :: jend
integer :: ub_month_n ! # of month in ubc file
integer :: ub_lev_n ! # of lev in ubc file
integer :: ub_lat_n ! # of lat in ubc file
integer :: ub_species_n ! # of species in ubc file
integer :: ub_nchar_n ! # of nchar in ubc file
integer :: max_dom ! wrf domain count
real , allocatable :: ub_p_m(:) ! ubc midpoint levels pressure(Pa)
real , allocatable :: ub_p_e(:) ! ubc edge levels pressure(Pa)
real , allocatable :: ub_vmr(:,:,:,:) ! values of vmr in ubc file
real , allocatable :: ub_lat(:) ! values of lat in ubc file
integer , allocatable :: ub_map(:) ! species indices for ubc species
integer :: nox_ndx = -1
integer :: ox_ndx = -1
!----------------------------------------------------
! private variables
!----------------------------------------------------
type chem_upper_bc
real, pointer :: vmr(:,:,:,:,:) ! values of vmr in ubc file
logical :: is_allocated
end type chem_upper_bc
type(chem_upper_bc), private, allocatable :: upper_bc(:)
contains
!==========================================================================
subroutine upper_bc_driver( id, dtstep, current_date_char, &
chem, p_phy, p8w, tropo_lev, &
ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
its, ite, jts, jte, kts, kte )
use module_configure
use module_state_description
use module_model_constants
!----------------------------------------------------
! input arguments
!----------------------------------------------------
integer, intent(in) :: id, &
ids,ide, jds,jde, kds,kde, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte
real, dimension(ims:ime, kms:kme, jms:jme ), &
intent(in) :: p_phy, & ! p at mid-level (Pa)
p8w ! p at interface (Pa)
integer, dimension(ims:ime, jms:jme ), &
intent(in) :: tropo_lev ! tropopause level
real, intent(in) :: dtstep
CHARACTER (LEN=256) , intent(in) :: current_date_char
real, dimension(ims:ime,kms:kme,jms:jme,num_chem ), &
intent(inout) :: chem
!----------------------------------------------------
! local variables
!----------------------------------------------------
integer :: next, last
integer :: kk, km
integer :: kmax
integer :: lev_relax
integer :: astat
integer :: i, j, k, m, n
integer :: ku(kts:kte)
integer :: kl(kts:kte)
integer :: del_p(kts:kte)
real :: del_time
real :: vmr_relax
real :: nox_ubc, xno, xno2, rno
real :: pint_vals(2)
real, allocatable :: table_ox(:)
real, allocatable :: p8w_temp(:)
real, allocatable :: chem_temp(:)
character(len=132) :: message
!-----------------------------------------------------------------------
! tile dimensions, needed for nest domains!!
!-------------------------------------------
iend = min(ite,ide-1)
jend = min(jte,jde-1)
!------------------------------------------------------------------------
!... setup for time interpolation
!------------------------------------------------------------------------
call get_time_interp_factors ( current_date_char, last, next, del_time )
!--------------------------------------------------------
! ... setup for ox conserving interp
!--------------------------------------------------------
if( ox_ndx > 0 ) then
allocate( table_ox(ub_lev_n-2),stat=astat )
if( astat /= 0 ) then
write(message,*) 'upper_bc_driver: table_ox allocation error = ',astat
CALL wrf_message( trim(message) )
call wrf_abort
end if
end if
j_tile_loop : &
do j = jts, jend
i_tile_loop : &
do i = its, iend
!--------------------------------------------------------
! setup for the pressure interpolation
!--------------------------------------------------------
lev_loop : &
do k = kte,kts,-1 ! from top to bottom
if( p_phy(i,k,j) <= ub_p_m(1) ) then
kl(k) = 1
ku(k) = 1
del_p(k) = 0.
else if( p_phy(i,k,j) >= ub_p_m(ub_lev_n) ) then
kl(k) = ub_lev_n
ku(k) = ub_lev_n
del_p(k) = 0.
else
do kk = 2,ub_lev_n ! from top to bottom of ub_p_m
if( p_phy(i,k,j) <= ub_p_m(kk) ) then
ku(k) = kk
kl(k) = kk - 1
del_p(k) = log( p_phy(i,k,j)/ub_p_m(kl(k)) ) &
/ log( ub_p_m(ku(k))/ub_p_m(kl(k)) )
exit
end if
end do
end if
end do lev_loop
!--------------------------------------------------------
! find max level less than fixed_ubc_press (default 50 mb)
! fixed UB vals from top of model to this level
!--------------------------------------------------------
if( p_phy(i,kte,j) > fixed_ubc_press ) then
kmax = kte + 1
else
do kmax = kte,kts+1,-1 ! from top to bottom
if( p_phy(i,kmax,j) > fixed_ubc_press ) then
exit
end if
end do
kmax = kmax + 1
endif
!--------------------------------------------------------
! above fixed_ubc_press (default is 50 mb),
! set the mixing ratios as fixed values
!--------------------------------------------------------
species_loop : &
do m = 1,ub_species_n
ub_overwrite : &
if( ub_map(m) > 0 .and. kmax <= kte ) then
!-----------------------------------------------------------
! o3
!-----------------------------------------------------------
if( m == ox_ndx ) then
table_ox(1:ub_lev_n-2) = & ! from top to bottom
upper_bc(id)%vmr(i,j,m,last,2:ub_lev_n-1) &
+ (upper_bc(id)%vmr(i,j,m,next,2:ub_lev_n-1) &
- upper_bc(id)%vmr(i,j,m,last,2:ub_lev_n-1))*del_time
km = kte - kmax + 1
allocate ( p8w_temp (km+1) )
allocate ( chem_temp(km) )
!-----------------------------------------------------------
! upper_bc_rebin requires "top-down" ordering
!-----------------------------------------------------------
p8w_temp (1:km+1) = p8w (i,kte:kmax-1:-1,j)
call upper_bc_rebin( ub_lev_n-2, km, ub_p_e, p8w_temp, &
table_ox, chem_temp )
chem(i,kmax:kte,j,p_o3) = chem_temp(km:1:-1)
deallocate( p8w_temp )
deallocate( chem_temp )
cycle species_loop
end if
!-----------------------------------------------------------
! others
!-----------------------------------------------------------
lev_loop_a : &
do k = kte,kmax,-1 ! from top to bottom
pint_vals(1) = upper_bc(id)%vmr(i,j,m,last,kl(k)) &
+ (upper_bc(id)%vmr(i,j,m,last,ku(k)) &
- upper_bc(id)%vmr(i,j,m,last,kl(k)))*del_p(k)
pint_vals(2) = upper_bc(id)%vmr(i,j,m,next,kl(k)) &
+ (upper_bc(id)%vmr(i,j,m,next,ku(k)) &
- upper_bc(id)%vmr(i,j,m,next,kl(k)))*del_p(k)
!-----------------------------------------------------------
! others: hno3, ch4, co, n2o, n2o5
!-----------------------------------------------------------
if( m /= nox_ndx ) then
chem(i,k,j,ub_map(m)) = pint_vals(1) &
+ (pint_vals(2) - pint_vals(1))*del_time
!-----------------------------------------------------------
! others: no and no2
!-----------------------------------------------------------
else
nox_ubc = pint_vals(1) + del_time * (pint_vals(2) - pint_vals(1))
if( p_no >= param_first_scalar ) then
xno = chem(i,k,j,p_no)
else
xno = 0.
end if
if( p_no2 >= param_first_scalar ) then
xno2 = chem(i,k,j,p_no2)
else
xno2 = 0.
end if
if( xno > 0. .or. xno2 > 0.) then
rno = xno / (xno + xno2)
end if
if( p_no >= param_first_scalar ) then
chem(i,k,j,p_no) = rno * nox_ubc
end if
if( p_no2 >= param_first_scalar ) then
chem(i,k,j,p_no2) = (1. - rno) * nox_ubc
end if
end if
end do lev_loop_a
end if ub_overwrite
end do species_loop
ub_relax : &
if( tropo_lev(i,j) > 0 .and. tropo_lev(i,j) < kmax ) then
!--------------------------------------------------------
!... relax lower stratosphere to extended ubc
! check to make sure ubc is not being imposed too low
! lev_relax = lowest model level (highest pressure)
! in which to relax to ubc
!--------------------------------------------------------
lev_relax = tropo_lev(i,j)
do while( p_phy(i,lev_relax,j) > ub_p_m(ub_lev_n) ) ! ub_p_m(ub_lev_n)= 440 mb
lev_relax = lev_relax + 1 ! move level upward
end do
if( lev_relax /= tropo_lev(i,j) ) then
write(message,*) 'upper_bc_driver:Warning,raised ubc: at j,i,= ',j,i
call wrf_message( trim(message) )
write(message,*) 'from ', tropo_lev(i,j),nint(p_phy(i,tropo_lev(i,j)-1,j)/mb2pa),' mb'
call wrf_message( trim(message) )
write(message,*) 'to ', lev_relax, nint(p_phy(i,lev_relax,j) /mb2pa),' mb'
call wrf_message( trim(message) )
endif
species_loop_a : &
do m = 1,ub_species_n
has_ubc : &
if( ub_map(m) > 0 ) then
lev_loop_b: do k = kmax-1,lev_relax,-1 ! from top to bottom
pint_vals(1) = upper_bc(id)%vmr(i,j,m,last,kl(k)) &
+ (upper_bc(id)%vmr(i,j,m,last,ku(k)) &
- upper_bc(id)%vmr(i,j,m,last,kl(k)))*del_p(k)
pint_vals(2) = upper_bc(id)%vmr(i,j,m,next,kl(k)) &
+ (upper_bc(id)%vmr(i,j,m,next,ku(k)) &
- upper_bc(id)%vmr(i,j,m,next,kl(k)))*del_p(k)
vmr_relax = pint_vals(1) &
+ (pint_vals(2) - pint_vals(1))*del_time
if( m /= nox_ndx ) then
chem(i,k,j,ub_map(m)) = chem(i,k,j,ub_map(m)) &
+ (vmr_relax - chem(i,k,j,ub_map(m))) * fac_relax
else
if( p_no >= param_first_scalar ) then
xno = chem(i,k,j,p_no)
else
xno = 0.
endif
if( p_no2 >= param_first_scalar ) then
xno2 = chem(i,k,j,p_no2)
else
xno2 = 0.
endif
if( xno > 0. .or. xno2 > 0. ) then
rno = xno / (xno + xno2)
endif
nox_ubc = xno + xno2
nox_ubc = nox_ubc + (vmr_relax - nox_ubc) * fac_relax
if( p_no >= param_first_scalar ) then
chem(i,k,j,p_no) = rno * nox_ubc
endif
if( p_no2 >= param_first_scalar ) then
chem(i,k,j,p_no2)= (1. - rno) * nox_ubc
endif
endif
end do lev_loop_b
endif has_ubc
end do species_loop_a
endif ub_relax
end do i_tile_loop
end do j_tile_loop
if( allocated( table_ox ) ) then
deallocate( table_ox )
endif
end subroutine upper_bc_driver
!===========================================================================
subroutine upper_bc_init(id, xlat, dt, config_flags, &
ids,ide, jds,jde, kds,kde, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte )
!---------------------------------------------------------------------
! ... new initialization routine for upper boundary condition
!---------------------------------------------------------------------
use module_domain ! needed for p_co, p_o3....
use module_state_description, only : param_first_scalar
use module_configure, only : grid_config_rec_type
implicit none
!---------------------------------------------------------------------
! ... dummy arguments
!---------------------------------------------------------------------
integer, intent(in) :: id, &
ids,ide, jds,jde, kds,kde, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte
real, intent(in) :: dt
real, intent(in) :: xlat(ims:ime, jms:jme) ! unit: degree
type(grid_config_rec_type), intent(in) :: config_flags
!---------------------------------------------------------------------
! ... local variables
!---------------------------------------------------------------------
integer :: astat
integer :: ncid
integer :: varid
integer :: dimid(4)
integer :: start(4)
integer :: count(4)
integer :: dimid_lev
integer :: dimid_lat
integer :: dimid_month
integer :: dimid_species
integer :: ndims
character(len=132) :: message
character(len=128) :: err_msg
character(len=64) :: filename
character(len=3) :: id_num
character(len=25), allocatable :: ub_specname(:)
character(len=80) :: attribute
integer :: lati
integer, allocatable :: lat_interp(:,:)
real, allocatable :: lat_del(:,:)
integer :: i, j, k, m
integer :: j1, j2
LOGICAL , EXTERNAL :: wrf_dm_on_monitor
!-------------------------------------------------------------------------
#ifdef NETCDF
include 'netcdf.inc'
#else
call wrf_message( 'upper_bc_init: requires netcdf' )
call wrf_abort
#endif
!-------------------------------------
!... note: grid%XLAT(ide,:) =0.
!... note: grid%XLONG(:,jde) =0.
!... => iend = min(ite,ide-1)
! jend = min(jte,jde-1)
!-------------------------------------
iend = min(ite,ide-1)
jend = min(jte,jde-1)
!---------------------------------------------------------------------
! ... allocate chem_upper_bc type
!---------------------------------------------------------------------
if( id == 1 .and. .not. allocated(upper_bc) ) then
CALL nl_get_max_dom( 1,max_dom )
allocate( upper_bc(max_dom),stat=astat )
if( astat /= 0 ) then
CALL wrf_message( 'upper_bc_init: failed to allocate upper_bc' )
CALL wrf_abort
end if
upper_bc(:)%is_allocated = .false.
endif
upper_bc_allocated : &
if( .not. upper_bc(id)%is_allocated ) then
!---------------------------------------------------------------------
!... allocate upper_bc type
!--------------------------------------------------------------------
is_d01 : &
IF( id == 1 ) then
master_proc : &
IF( wrf_dm_on_monitor() ) THEN
write(id_num,'(i3)') 100+id
write(message,*) 'upper_bc_init: intializing domain ' // id_num(2:3)
CALL wrf_message( trim(message) )
!---------------------------------------------------------------------
! ... open upper_bc netcdf file
!---------------------------------------------------------------------
filename = config_flags%fixed_ubc_inname
if( filename == ' ' ) then
call wrf_message( 'upper_bc_init: input filename not specified in namelist' )
call wrf_abort
endif
err_msg = 'upper_bc_init: failed to open file ' // trim(filename)
call handle_ncerr( nf_open( trim(filename), nf_noclobber, ncid ), trim(err_msg) )
write(message,*) 'upper_bc_init: id, open filename= ',id, filename
CALL wrf_message( trim(message) )
!---------------------------------------------------------------------
! ... get dimensions
!---------------------------------------------------------------------
err_msg = 'upper_bc_init: failed to get month id'
call handle_ncerr( nf_inq_dimid( ncid, 'month', dimid_month ), trim(err_msg) )
err_msg = 'upper_bc_init: failed to get month'
call handle_ncerr( nf_inq_dimlen( ncid, dimid_month, ub_month_n ), trim(err_msg) )
err_msg = 'upper_bc_init: failed to get lat id'
call handle_ncerr( nf_inq_dimid( ncid, 'lat', dimid_lat ), trim(err_msg) )
err_msg = 'upper_bc_init: failed to get lat'
call handle_ncerr( nf_inq_dimlen( ncid, dimid_lat, ub_lat_n ), trim(err_msg) )
err_msg = 'upper_bc_init: failed to get lev id'
call handle_ncerr( nf_inq_dimid( ncid, 'lev', dimid_lev ), trim(err_msg) )
err_msg = 'upper_bc_init: failed to get lev'
call handle_ncerr( nf_inq_dimlen( ncid, dimid_lev, ub_lev_n ), trim(err_msg) )
err_msg = 'upper_bc_init: failed to get species id'
call handle_ncerr( nf_inq_dimid( ncid, 'species', dimid_species ), trim(err_msg) )
err_msg = 'upper_bc_init: failed to get species'
call handle_ncerr( nf_inq_dimlen( ncid, dimid_species, ub_species_n ), trim(err_msg) )
err_msg = 'upper_bc_init: failed to get nchar id'
call handle_ncerr( nf_inq_dimid( ncid, 'nchar', dimid(1) ), trim(err_msg) )
err_msg = 'upper_bc_init: failed to get nchar'
call handle_ncerr( nf_inq_dimlen( ncid, dimid(1), ub_nchar_n ), trim(err_msg) )
END IF master_proc
!---------------------------------------------------------------------
! ... bcast the dimensions
!---------------------------------------------------------------------
#ifdef DM_PARALLEL
CALL wrf_dm_bcast_integer ( ub_month_n , 1 )
CALL wrf_dm_bcast_integer ( ub_lat_n , 1 )
CALL wrf_dm_bcast_integer ( ub_lev_n , 1 )
CALL wrf_dm_bcast_integer ( ub_species_n , 1 )
CALL wrf_dm_bcast_integer ( ub_nchar_n , 1 )
#endif
!---------------------------------------------------------------------
! ... allocate module arrays
!---------------------------------------------------------------------
allocate( ub_lat(ub_lat_n), stat=astat )
if( astat /= 0 ) then
call wrf_message( 'upper_bc_init: failed to allocate ub_lat' )
call wrf_abort
end if
allocate( ub_vmr(ub_lat_n,ub_species_n,ub_month_n,ub_lev_n), stat=astat )
if( astat /= 0 ) then
call wrf_message( 'upper_bc_init: failed to allocate ub_vmr' )
call wrf_abort
end if
allocate( ub_p_m(ub_lev_n), stat=astat )
if( astat /= 0 ) then
call wrf_message( 'upper_bc_init: failed to allocate ub_p_m' )
call wrf_abort
end if
allocate( ub_p_e(ub_lev_n-1), stat=astat )
if( astat /= 0 ) then
call wrf_message( 'upper_bc_init: failed to allocate ub_p_e' )
call wrf_abort
end if
allocate( ub_map(ub_species_n), stat=astat )
if( astat /= 0 ) then
call wrf_message( 'upper_bc_init: failed to allocate ub_map' )
call wrf_abort
end if
!---------------------------------------------------------------------
! ... read arrays
!---------------------------------------------------------------------
master_proc_a : &
IF ( wrf_dm_on_monitor() ) THEN
!lev
err_msg = 'upper_bc_init: failed to get lev variable id'
call handle_ncerr( nf_inq_varid( ncid, 'lev', varid ), trim(err_msg) )
err_msg = 'upper_bc_init: failed to read lev variable'
call handle_ncerr( nf_get_var_real( ncid, varid, ub_p_m ), trim(err_msg) )
!-------- check unit
attribute(:) = ' '
astat = nf_get_att_text( ncid, varid, 'units', attribute )
if( astat == nf_noerr )then
if( trim(attribute) == 'mb' .or. trim(attribute) == 'hpa' )then
write(message,*) 'upper_bc_init: units for lev = ',trim(attribute),'... converting to pa'
call wrf_message( trim(message) )
ub_p_m(:) = mb2Pa * ub_p_m(:)
else if( trim(attribute) /= 'pa' .and. trim(attribute) /= 'pa' )then
write(message,*) 'upper_bc_init: unknown units for lev, units=*',trim(attribute),'*'
call wrf_message( trim(message) )
call wrf_abort
end if
else
write(message,*) 'upper_bc_init: warning! units attribute for lev missing, assuming mb'
call wrf_message( trim(message) )
ub_p_m(:) = mb2Pa * ub_p_m(:)
end if
! ... set values of ub_p_e
ub_p_e(1:ub_lev_n-1) = .5*(ub_p_m(1:ub_lev_n-1) + ub_p_m(2:ub_lev_n))
!---------------------------
!lat
err_msg = 'upper_bc_init: failed to get lat variable id'
call handle_ncerr( nf_inq_varid( ncid, 'lat', varid ), trim(err_msg) )
err_msg = 'upper_bc_init: failed to read lat variable'
call handle_ncerr( nf_get_var_real( ncid, varid, ub_lat ), trim(err_msg) )
!---------------------------
!specname
allocate( ub_specname(ub_species_n), stat=astat )
if( astat /= 0 ) then
call wrf_message( 'upper_bc_init: failed to allocate ub_specname' )
call wrf_abort
end if
species_loop : &
do i = 1,ub_species_n
start(:2) = (/ 1, i /)
count(:2) = (/ ub_nchar_n, 1 /)
ub_specname(i)(:) = ' '
err_msg = 'upper_bc_init: failed to get specname variable id'
call handle_ncerr( nf_inq_varid( ncid, 'specname', varid ), trim(err_msg) )
err_msg = 'upper_bc_init: failed to read ub_specname variable'
call handle_ncerr( nf_get_vara_text( ncid, varid, start(:2), count(:2), ub_specname(i) ), trim(err_msg) )
ub_map(i) = 0
if( trim(ub_specname(i)) == 'HNO3' .and. p_hno3 >= param_first_scalar ) then
ub_map(i) = p_hno3
else if ( trim(ub_specname(i)) == 'CH4' .and. p_ch4 >= param_first_scalar ) then
ub_map(i) = p_ch4
else if ( trim(ub_specname(i)) == 'CO' .and. p_co >= param_first_scalar ) then
ub_map(i) = p_co
else if ( trim(ub_specname(i)) == 'N2O' .and. p_n2o >= param_first_scalar ) then
ub_map(i) = p_n2o
else if ( trim(ub_specname(i)) == 'N2O5' .and. p_n2o5 >= param_first_scalar ) then
ub_map(i) = p_n2o5
else if ( trim(ub_specname(i)) == 'OX' .and. p_o3 >= param_first_scalar ) then
ub_map(i) = p_o3
ox_ndx = i
else if ( trim(ub_specname(i)) == 'NOX' .and. &
(p_no >= param_first_scalar .or. p_no2 >= param_first_scalar) ) then
ub_map(i) = p_no
nox_ndx = i
endif
if( ub_map(i) == 0 ) then
write(message,*) 'upper_bc_init: ubc table species ',trim(ub_specname(i)), ' not used'
call wrf_message( trim(message) )
end if
end do species_loop
!--------------------------
!vmr
err_msg = 'upper_bc_init: failed to get vmr variable id'
call handle_ncerr( nf_inq_varid( ncid, 'vmr', varid ), trim(err_msg) )
! check dimensions
err_msg = 'upper_bc_init: failed to get ndims of vmr variable'
call handle_ncerr( nf_inq_varndims( ncid, varid, ndims ), trim(err_msg) )
if( ndims /= 4 ) then
write(message,*) 'upper_bc_init: error! variable vmr has ndims = ',ndims,', expecting 4'
call wrf_message( trim(message) )
call wrf_abort
end if
err_msg = 'upper_bc_init: failed to get dimid of vmr variable'
call handle_ncerr( nf_inq_vardimid( ncid, varid, dimid ), trim(err_msg) )
if( dimid(1) /= dimid_lat .or. dimid(2) /= dimid_species .or. &
dimid(3) /= dimid_month .or. dimid(4) /= dimid_lev )then
write(message,*) 'upper_bc_init: error! dimensions in wrong order for variable vmr,'// &
'expecting (lat,species,month,lev)'
call wrf_message( trim(message) )
call wrf_abort
end if
!------------------------------------------------------------------
! ... read in the ub mixing ratio values
!------------------------------------------------------------------
err_msg = 'upper_bc_init: failed to read vmr variable'
call handle_ncerr( nf_get_var_real( ncid, varid, ub_vmr ), trim(err_msg) )
!---------------------------------------------------------------------
!... close input netcdf file
!---------------------------------------------------------------------
err_msg = 'upper_bc_init: failed to close file ' // trim(filename)
call handle_ncerr( nf_close( ncid ), trim(err_msg) )
deallocate( ub_specname, stat=astat )
if( astat /= 0 ) then
write(message,*) ': failed to deallocate ub_specnmae; astat = ',astat
call wrf_message( trim(message) )
call wrf_abort
end if
ENDIF master_proc_a
!---------------------------------------------------------------------
! ... bcast the variables
!---------------------------------------------------------------------
#ifdef DM_PARALLEL
CALL wrf_dm_bcast_integer ( nox_ndx , 1 )
CALL wrf_dm_bcast_integer ( ox_ndx , 1 )
CALL wrf_dm_bcast_integer ( ub_map, size(ub_map) )
CALL wrf_dm_bcast_real ( ub_p_m, size(ub_p_m) )
CALL wrf_dm_bcast_real ( ub_p_e, size(ub_p_e) )
CALL wrf_dm_bcast_real ( ub_lat, size(ub_lat) )
CALL wrf_dm_bcast_real ( ub_vmr, size(ub_vmr) )
#endif
fixed_ubc_press = config_flags%fixed_ubc_press*mb2Pa
ENDIF is_d01
upper_bc(id)%is_allocated = .true.
write(message,*) 'upper_bc_init: ub_vmr(1,1,1,:)'
call wrf_message( trim(message) )
do k = 1,ub_lev_n,5
write(message,'(1p,5g15.7)') ub_vmr(1,1,1,k:min(k+4,ub_lev_n))
call wrf_message( trim(message) )
end do
!---------------------------------------------------------------------
! ... allocate domain specific arrays
!---------------------------------------------------------------------
allocate( lat_del(its:iend,jts:jend), stat=astat )
if( astat /= 0 ) then
call wrf_message( 'upper_bc_init: failed to allocate lat_del' )
call wrf_abort
end if
allocate( lat_interp(its:iend,jts:jend), stat=astat )
if( astat /= 0 ) then
call wrf_message( 'upper_bc_init: failed to allocate lat_interp' )
call wrf_abort
end if
!------------------------------------------------------------------
! ... regrid mixing ratio to model latitudes
!------------------------------------------------------------------
!-------------------------------------
!... get lat_del, lat_interp and lati
!-------------------------------------
call get_lat_interp_factors( id, xlat, jend, iend, &
ub_lat_n, ub_lat, &
lat_del, lat_interp, lati, &
ids,ide, jds,jde, kds,kde, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte )
!---------------------------------------------------------------------
! ... allocate upper_bc(id) component array
!---------------------------------------------------------------------
allocate( upper_bc(id)%vmr(its:iend,jts:jend,ub_species_n,ub_month_n,ub_lev_n), stat=astat )
if( astat /= 0 ) then
call wrf_message( 'upper_bc_init: failed to allocate upper_bc(id)%vrm' )
call wrf_abort
end if
!--------------------------------------------------------------------
!... regrid from ub_vmr to upper_bc(id)%vmr
!--------------------------------------------------------------------
do j = jts,jend
do i = its,iend
j1 = lat_interp(i,j)
j2 = j1 + lati
upper_bc(id)%vmr(i,j,:,:,:) = ub_vmr(j1,:,:,:) + &
lat_del(i,j) * (ub_vmr(j2,:,:,:) - ub_vmr(j1,:,:,:))
end do
end do
!----------------------------------------------------
! transform to ppm vmr
!----------------------------------------------------
upper_bc(id)%vmr(:,:,:,:,:) = upper_bc(id)%vmr(:,:,:,:,:)*vmr2ppmv
write(message,*) 'upper_bc_init: upper_bc(id)%vmr(its+5,jts+5,1,6,:)'
call wrf_message( trim(message) )
do k = 1,ub_lev_n,5
write(message,'(1p,5g15.7)') upper_bc(id)%vmr(its+5,jts+5,1,6,k:min(k+4,ub_lev_n))
call wrf_message( trim(message) )
end do
!--------------------------------------------------------
! ... set up the relaxation for lower stratosphere
!--------------------------------------------------------
! ... tau_relax = relaxation timescale (in sec)
! fac_relax = fractional relaxation towards ubc
! 1 => use ubc
! 0 => ignore ubc, use model concentrations
!--------------------------------------------------------
fac_relax = 1. - exp( -real(dt)/tau_relax )
!--------------------------------------------------------
if( id == max_dom ) then
deallocate( ub_vmr)
deallocate( ub_lat)
endif
deallocate( lat_del)
deallocate( lat_interp)
call wrf_message( ' ' )
write(message,*) 'upper_bc_init: DONE intialized domain ',id
call wrf_message( trim(message) )
call wrf_message( ' ' )
endif upper_bc_allocated
end subroutine upper_bc_init
!===========================================================================
subroutine handle_ncerr( ret, mes )
!---------------------------------------------------------------------
! ... netcdf error handling routine
!---------------------------------------------------------------------
implicit none
!---------------------------------------------------------------------
! ... dummy arguments
!---------------------------------------------------------------------
integer, intent(in) :: ret
character(len=*), intent(in) :: mes
#ifdef NETCDF
include 'netcdf.inc'
#endif
if( ret /= nf_noerr ) then
call wrf_message( trim(mes) )
call wrf_message( trim(nf_strerror(ret)) )
call wrf_abort
end if
end subroutine handle_ncerr
!===========================================================================
subroutine get_time_interp_factors ( current_date_char, last, next, del_time)
use module_date_time, only : get_julgmt
implicit none
!-----------------------------------------------------------------------------
! ... dummy arguments
!-----------------------------------------------------------------------------
CHARACTER (LEN=256),intent(in) :: current_date_char
integer, intent(out) :: next, last
real, intent(out) :: del_time
!-----------------------------------------------------------------------------
! ... local variables
!-----------------------------------------------------------------------------
integer, parameter :: day_of_year(12) = (/ 16, 45, 75, 105, 136, 166, 197, &
228, 258, 289, 319, 350 /)
INTEGER :: julyr , julday
REAL :: gmt
real :: calday
integer :: m
!------------------------------------------------------------------------------
call get_julgmt ( current_date_char, julyr, julday, gmt )
calday = real(julday) + gmt
if( calday < day_of_year(1) ) then
next = 1
last = 12
del_time = (365. + calday - day_of_year(12)) &
/ (365. + day_of_year(1) - day_of_year(12))
else if( calday >= day_of_year(12) ) then
next = 1
last = 12
del_time = (calday - day_of_year(12)) &
/ (365. + day_of_year(1) - day_of_year(12))
else
do m = 11,1,-1
if( calday >= day_of_year(m) ) then
exit
end if
end do
last = m
next = m + 1
del_time = (calday - day_of_year(m)) / (day_of_year(m+1) - day_of_year(m))
end if
del_time = max( min( 1.,del_time ),0. )
end subroutine get_time_interp_factors
!===========================================================================
subroutine upper_bc_rebin( nsrc, ntrg, src_x, trg_x, src, trg )
!---------------------------------------------------------------
! ... rebin src to trg
!---------------------------------------------------------------
implicit none
!---------------------------------------------------------------
! ... dummy arguments
!---------------------------------------------------------------
integer, intent(in) :: nsrc ! dimension source array
integer, intent(in) :: ntrg ! dimension target array
real, intent(in) :: src_x(nsrc+1) ! source coordinates
real, intent(in) :: trg_x(ntrg+1) ! target coordinates
real, intent(in) :: src(nsrc) ! source array
real, intent(out) :: trg(ntrg) ! target array
!---------------------------------------------------------------
! ... local variables
!---------------------------------------------------------------
integer :: i, l
integer :: si, si1
integer :: sil, siu
real :: y
real :: sl, su
real :: tl, tu
do i = 1,ntrg
tl = trg_x(i)
if( tl < src_x(nsrc+1) ) then
do sil = 1,nsrc+1
if( tl <= src_x(sil) ) then
exit
end if
end do
tu = trg_x(i+1)
do siu = 1,nsrc+1
if( tu <= src_x(siu) ) then
exit
end if
end do
y = 0.
sil = max( sil,2 )
siu = min( siu,nsrc+1 )
do si = sil,siu
si1 = si - 1
sl = max( tl,src_x(si1) )
su = min( tu,src_x(si) )
y = y + (su - sl)*src(si1)
end do
trg(i) = y/(trg_x(i+1) - trg_x(i))
else
trg(i) = 0.
end if
end do
end subroutine upper_bc_rebin
!===========================================================================
subroutine get_lat_interp_factors( id, xlat, jend, iend, &
from_lat_n, from_lat, &
lat_del, lat_interp, lati, &
ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
its, ite, jts, jte, kts, kte )
!---------------------------------------------------------------------
! ... Determine indicies and deltas for transform
! Note : it is assumed that the latitude and longitude
! arrays are monotonic
!--------------------------------------------------------------------
implicit none
!---------------------------------------------------------------------
! ... Dummy args
!---------------------------------------------------------------------
integer, intent(in) :: id, &
ids,ide, jds,jde, kds,kde, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte
real, intent(in) :: xlat(ims:ime, jms:jme) ! unit: degree
integer, intent(in) :: jend
integer, intent(in) :: iend
integer, intent(in) :: from_lat_n
real , intent(in) :: from_lat(from_lat_n)
integer, intent(out) :: lati
integer, intent(out) :: lat_interp(its:iend,jts:jend)
real, intent(out) :: lat_del(its:iend,jts:jend)
!---------------------------------------------------------------------
! ... Local variables
!---------------------------------------------------------------------
integer :: astat
real :: to_lat(jts:jend) ! values of wrf lat, 1d
integer :: to_lon_n
real :: countx
real :: target_lat
integer :: latl, latu
real :: max_lat, min_lat
logical :: from_lats_mono_pos
integer :: i, j
integer :: m, n
!---------------------------------------------------------------------
!... set "module" variable, converter(id)%lati
!---------------------------------------------------------------------
max_lat = from_lat(from_lat_n)
min_lat = from_lat(1)
if( from_lat(from_lat_n) >= from_lat(1)) then ! increasing lat
latl = 1
latu = from_lat_n
lati = 1
from_lats_mono_pos = .true.
else ! decreasing lat
latl = from_lat_n
latu = 1
lati = -1
from_lats_mono_pos = .false.
end if
!---------------------------------------------------------------------
! ... Set interpolation latitude indicies and del_lat
!---------------------------------------------------------------------
do j = jts,jend
do i = its,iend
target_lat = xlat(i,j)
if( target_lat <= min_lat ) then
lat_del(i,j) = 0.
lat_interp(i,j) = latl
else if( target_lat >= max_lat ) then
lat_del(i,j) = 1.
if( from_lats_mono_pos ) then
lat_interp(i,j) = latu - 1
else
lat_interp(i,j) = latu + 1
end if
else
do m = latl,latu,lati
if( target_lat < from_lat(m) ) then
n = m - lati
lat_interp(i,j) = min( from_lat_n,max( 1,n ) )
lat_del(i,j) = &
(target_lat - from_lat(n))/(from_lat(m) - from_lat(n))
exit
end if
end do
end if
end do
end do
end subroutine get_lat_interp_factors
end module module_upper_bc_driver