!#define SW_DEBUG
module module_phot_tuv
use params_mod, only : dp, m2km, ppm2vmr, o2vmr, km2cm, m2s
IMPLICIT NONE
private
public :: tuv_driver, tuv_init, tuv_timestep_init
real, parameter :: conv1 = km2cm*.5
real, parameter :: conv2 = conv1*ppm2vmr
real, parameter :: bext340 = 5.E-6
real, parameter :: bexth2o = 5.E-6
integer :: nj
integer :: nlambda_start = 1
integer :: nlambda_af_start = 1
integer :: nlambda_af_end = 1
integer :: nconc, ntemp, nwave
integer :: n_temp_data, n_o3_data, n_air_dens_data
integer :: j_o2_ndx = -1
integer :: last, next
integer :: curjulday = 0
integer, allocatable :: rxn_ndx(:)
real :: dels
real :: esfact
logical :: has_exo_coldens
logical :: tuv_is_initialized = .false.
real, allocatable :: temp_tab(:)
real, allocatable :: conc_tab(:)
real, allocatable :: del_temp_tab(:)
real, allocatable :: del_conc_tab(:)
real, allocatable :: wl(:)
real, allocatable :: wc(:)
real, allocatable :: dw(:)
real, allocatable :: w_fac(:)
real, allocatable :: etfl(:)
real, allocatable :: albedo(:)
real, allocatable :: o2_xs(:)
real, allocatable :: so2_xs(:)
real, allocatable :: par_wght(:)
real, allocatable :: ery_wght(:)
real, allocatable :: z_temp_data(:), temp_data(:)
real, allocatable :: z_o3_data(:), o3_data(:)
real, allocatable :: z_air_dens_data(:), air_dens_data(:)
real, allocatable :: o3_xs_tab(:,:)
real, allocatable :: no2_xs_tab(:,:)
real, allocatable :: xsqy_tab(:,:,:,:)
character(len=32), allocatable :: tuv_jname(:)
logical :: has_so2, has_no2
logical :: is_full_tuv = .true.
! logical :: is_full_tuv = .false.
logical :: rxn_initialized
logical, allocatable :: xsqy_is_zdep(:)
type column_density
integer :: ncoldens_levs
integer :: ndays_of_year
real(8), pointer :: col_levs(:)
real(8), pointer :: day_of_year(:)
real(8), pointer :: o3_col_dens(:,:,:,:)
real(8), pointer :: o2_col_dens(:,:,:,:)
logical :: is_allocated
end type column_density
type(column_density), allocatable :: col_dens(:)
CONTAINS
subroutine tuv_driver( &
dm, curr_secs, ktau, config_flags, haveaer, &
gmt, julday, t_phy, moist, aerwrf, &
p8w, t8w, p_phy, chem, rho_phy, &
dz8w, xlat, xlong, z, z_at_w, gd_cloud, gd_cloud2, &
tauaer1,tauaer2,tauaer3,tauaer4, &
gaer1,gaer2,gaer3,gaer4, &
waer1,waer2,waer3,waer4, &
ph_macr,ph_o31d,ph_o33p,ph_no2,ph_no3o2,ph_no3o,ph_hno2, &
ph_hno3,ph_hno4,ph_h2o2,ph_ch2or,ph_ch2om,ph_ch3cho, &
ph_ch3coch3,ph_ch3coc2h5,ph_hcocho,ph_ch3cocho, &
ph_hcochest,ph_ch3o2h,ph_ch3coo2h,ph_ch3ono2,ph_hcochob, &
ph_n2o5,ph_o2,ph_n2o,ph_pooh,ph_pan,ph_mvk,ph_hyac, &
ph_glyald,ph_mek,ph_gly, &
pm2_5_dry,pm2_5_water,uvrad, &
dt_cld,af_dir,af_dn,af_up,par,erythema, &
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
USE module_data_radm2
USE tuv_subs, only : tuv_radfld, sundis, calc_zenith
USE module_params, only : kz
USE srb, only : sjo2
USE module_rxn, only : xsqy_table => xsqy_tab, the_subs, set_initialization
USE module_rxn, only : get_initialization
USE module_xsections, only : o3xs, no2xs_jpl06a
!---------------------------------------------------------------------
! ... dummy arguments
!---------------------------------------------------------------------
INTEGER, INTENT(IN ) :: dm,julday, &
ids,ide, jds,jde, kds,kde, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte
INTEGER, INTENT(IN ) :: ktau
REAL(KIND=8), INTENT(IN ) :: curr_secs
REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_moist ), &
INTENT(IN ) :: moist
REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
INTENT(INOUT ) :: &
pm2_5_dry,pm2_5_water
REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
INTENT(INOUT ) :: &
gd_cloud,gd_cloud2
REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
INTENT(IN ) :: tauaer1, tauaer2, tauaer3, tauaer4, &
waer1, waer2, waer3, waer4, &
gaer1, gaer2, gaer3, gaer4
REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
INTENT(INOUT ) :: &
ph_macr,ph_o31d,ph_o33p,ph_no2,ph_no3o2,ph_no3o,ph_hno2,&
ph_hno3,ph_hno4,ph_h2o2,ph_ch2or,ph_ch2om,ph_ch3cho, &
ph_ch3coch3,ph_ch3coc2h5,ph_hcocho,ph_ch3cocho, &
ph_hcochest,ph_ch3o2h,ph_ch3coo2h,ph_ch3ono2,ph_hcochob,&
ph_n2o5,ph_o2,ph_n2o,ph_pooh,ph_pan,ph_mvk,ph_hyac, &
ph_glyald,ph_mek,ph_gly
REAL, INTENT(INOUT) :: dt_cld(ims:ime,kms:kme,jms:jme), &
af_dir(ims:ime,kms:kme,jms:jme), &
af_dn(ims:ime,kms:kme,jms:jme), &
af_up(ims:ime,kms:kme,jms:jme), &
par(ims:ime,kms:kme,jms:jme), &
erythema(ims:ime,kms:kme,jms:jme)
REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ), &
INTENT(IN ) :: chem
REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
INTENT(IN ) :: &
z, &
t_phy, &
p_phy, &
dz8w, &
t8w,p8w,z_at_w , &
aerwrf , &
rho_phy
REAL, DIMENSION( ims:ime , jms:jme ) , &
INTENT(IN ) :: &
xlat, &
xlong
REAL, DIMENSION( ims:ime , jms:jme ) , &
INTENT(INOUT ) :: uvrad
REAL, INTENT(IN ) :: gmt
TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
LOGICAL, INTENT(IN) :: haveaer
!---------------------------------------------------------------------
! ... local variables
!---------------------------------------------------------------------
real, parameter :: boltz = 1.38065e-23 ! J/K/molecule
real, parameter :: xair_kg = 1.e3/28.97*6.023e23*1.e-6
real(dp), parameter :: Pa2hPa = .01_dp
integer :: astat, ierr
integer :: k, ki, kp1, i, j, n, wn
integer :: ktsm1, ktem1
integer :: nlev, nlyr
integer :: jndx
integer :: n_tuv_z, n_exo_z, last_n_tuv_z
integer :: min_ndx(2), max_ndx(2)
integer :: myproc
integer(kind=8) :: ixhour
integer :: minndx(2)
real :: xmin, gmtp, uvb_dd1, uvb_du1, uvb_dir1
real :: zenith, dobsi
real :: delz_top
real(kind=8) :: xtime, xhour
real :: max_tauaer
real :: Dobson(2)
real :: xsect(nwave)
real :: wrk(nwave), wrk_lam(nwave)
! real :: tlev(kts-1:kte)
! real :: tlyr(kts-1:kte)
! real :: o33(kts-1:kte)
real :: rhoa(kts-1:kte)
! real :: dens_air(kts-1:kte)
real :: aerext(kts-1:kte)
! real :: qll(kts-1:kte)
! real :: aircol(kts-1:kte)
! real :: o3col(kts-1:kte)
! real :: o2col(kts-1:kte)
! real :: so2col(kts-1:kte)
! real :: no2col(kts-1:kte)
real :: dtcld_col(kts:kte)
! real :: cldfrac(kts-1:kte)
real :: par_col(kts:kte)
! real :: tauaer300(kts-1:kte), tauaer400(kts-1:kte), &
! tauaer600(kts-1:kte), tauaer999(kts-1:kte)
! real :: waer300(kts-1:kte), waer400(kts-1:kte), &
! waer600(kts-1:kte), waer999(kts-1:kte)
! real :: gaer300(kts-1:kte), gaer400(kts-1:kte), &
! gaer600(kts-1:kte), gaer999(kts-1:kte)
real :: zen_angle(ims:ime,jms:jme)
real :: z_top(its:ite,jts:jte)
real :: z_exo(1), dens_exo(1), temp_exo(1)
real :: dummy(kz)
real(dp) :: p_top(its:ite,jts:jte)
real(dp) :: o2_exo_col(its:ite,jts:jte)
real(dp) :: o3_exo_col(its:ite,jts:jte)
real(dp) :: o3_exo_col_at_grnd(its:ite,jts:jte)
logical :: do_alloc
logical :: has_aer_ra_feedback
real, allocatable :: tuv_z(:)
real, allocatable :: dtuv_z(:)
real, allocatable :: cldfrac(:)
real, allocatable :: qll(:)
real, allocatable :: tlev(:)
real, allocatable :: tlyr(:)
real, allocatable :: dens_air(:)
real, allocatable :: o33(:)
real, allocatable :: aircol(:)
real, allocatable :: o3col(:)
real, allocatable :: o2col(:)
real, allocatable :: so2col(:)
real, allocatable :: no2col(:)
real, allocatable :: tauaer300(:), tauaer400(:), tauaer600(:), tauaer999(:)
real, allocatable :: waer300(:), waer400(:), waer600(:), waer999(:)
real, allocatable :: gaer300(:), gaer400(:), gaer600(:), gaer999(:)
real, allocatable :: dtcld(:,:), dtaer(:,:)
real, allocatable :: omcld(:,:), omaer(:,:)
real, allocatable :: gcld(:,:), gaer(:,:)
real, allocatable :: rad_fld(:,:)
real, allocatable :: e_fld(:,:)
real, allocatable :: tuv_prate(:,:)
real, allocatable :: xsqy(:,:)
real, allocatable :: srb_o2_xs(:,:)
real, allocatable :: o3_xs(:,:), o3_xs_tpose(:,:)
real, allocatable :: no2_xs(:,:), no2_xs_tpose(:,:)
real, allocatable :: rad_fld_tpose(:,:)
real, allocatable :: dir_fld(:,:), dwn_fld(:,:), up_fld(:,:)
real, allocatable :: e_dir(:,:), e_dn(:,:), e_up(:,:)
REAL :: zexo_grd(2*kte)
CHARACTER(len=256) :: msg
#ifdef SW_DEBUG
logical :: tuv_diags
#endif
call wrf_get_myproc( myproc )
has_aer_ra_feedback = config_flags%aer_ra_feedback == 1
xtime = curr_secs/60._8 ! min since simulation start
ixhour = int( gmt+.01,8 ) + int( xtime/60._8,8 )
xhour = real( ixhour,8 )
xmin = 60.*gmt + real( xtime-xhour*60._8,4 )
gmtp = real( mod(xhour,24._8),4 )
gmtp = gmtp + xmin/60.
call calc_zenith( xlat, -xlong, &
julday, real(gmtp,8), zen_angle, &
its, ite, jts, jte, &
ims, ime, jms, jme )
#ifdef SW_DEBUG
if( myproc == 6 ) then
write(*,'(''tuv_drvr: its,jts = '',2i5)') its,jts
minndx(:) = minloc( zen_angle(its:ite,jts:jte) )
write(*,'(''tuv_drvr: min zen_angle ndx = '',2i5)') minndx(:)
write(*,'(''tuv_drvr: min zen_angle = '',1p,g15.8)') minval( zen_angle(its:ite,jts:jte) )
minndx(1) = minndx(1) + its - 1
minndx(2) = minndx(2) + jts - 1
write(*,'(''tuv_drvr: minndx = '',2i5)') minndx(:)
write(*,'(''tuv_drvr: min zen_angle = '',1p,g15.8)') zen_angle(minndx(1),minndx(2))
endif
#endif
do j = jts,jte
where( zen_angle(its:ite,j) == 90. )
zen_angle(its:ite,j) = 89.9
endwhere
end do
ktsm1 = kts - 1 ; ktem1 = kte - 1
allocate( tuv_prate(kts:kte,nj),stat=ierr )
if( ierr /= 0 ) then
call wrf_error_fatal( 'tuv_driver: failed to allocate tuv_prate' )
endif
any_daylight: &
if( any( zen_angle(its:ite,jts:jte) < 90. ) ) then
if( config_flags%has_o3_exo_coldens ) then
! p_top(its:ite,jts:jte) = Pa2hPa * real( p_phy(its:ite,kte,jts:jte),dp )
z_exo(1) = 50.
call z_interp( z_air_dens_data, air_dens_data, n_air_dens_data, &
z_exo, dens_exo )
call z_interp( z_temp_data, temp_data, n_temp_data, z_exo, temp_exo )
p_top(its:ite,jts:jte) = temp_exo(1)*boltz*dens_exo(1)*1.e6*Pa2hPa
call tuv_timestep_init( dm, julday )
call p_interp( o2_exo_col, o3_exo_col, o3_exo_col_at_grnd, p_top, &
dm, its, ite, jts, jte )
endif
allocate( rad_fld(nwave,kts:kte),e_fld(kts:kte,nwave),stat=astat )
ierr = ierr + astat
allocate( dir_fld(kts:kte,nwave), dwn_fld(kts:kte,nwave), up_fld(kts:kte,nwave),stat=astat )
ierr = ierr + astat
allocate( e_dir(kts:kte,nwave), e_dn(kts:kte,nwave), e_up(kts:kte,nwave),stat=astat )
ierr = ierr + astat
if( .not. is_full_tuv ) then
if( any( .not. xsqy_is_zdep(:) ) ) then
allocate( rad_fld_tpose(kts:kte,nwave),stat=astat )
endif
elseif( any( xsqy_table(2:nj)%tpflag == 0 ) ) then
allocate( rad_fld_tpose(kts:kte,nwave),stat=astat )
endif
ierr = ierr + astat
allocate( srb_o2_xs(nwave,kts:kte) )
ierr = ierr + astat
allocate( xsqy(nwave,kts:kte) )
ierr = ierr + astat
if( ierr /= 0 ) then
call wrf_error_fatal( 'tuv_driver: failed to allocate rad_fld ... xsqy' )
endif
!-----------------------------------------------------------------------------
! set solar distance factor
!-----------------------------------------------------------------------------
if( curjulday /= julday ) then
curjulday = julday
esfact = sundis( julday )
endif
if( .not. config_flags%scale_o3_to_grnd_exo_coldens ) then
if( config_flags%scale_o3_to_du_at_grnd ) then
dobsi = max( 0.,config_flags%du_at_grnd )
else
dobsi = 0.
endif
endif
endif any_daylight
if( ierr /= 0 ) then
call wrf_error_fatal( 'tuv_driver: failed to allocate module variables' )
endif
z_top(:,:) = 0.
do j = jts,jte
do i = its,ite
if( zen_angle(i,j) < 90. ) then
wrk(1) = m2km*(z(i,kte,j) - z_at_w(i,kts,j))
z_top(i,j) = real( ceiling(wrk(1)) )
delz_top = z_top(i,j) - wrk(1)
if( delz_top < .3 ) then
z_top(i,j) = z_top(i,j) + 1.
endif
endif
end do
end do
if( is_full_tuv ) then
rxn_initialized = .not. get_initialization()
if( .not. rxn_initialized ) then
do n = 1,nj
jndx = rxn_ndx(n)
if( jndx /= -1 ) then
call the_subs(jndx)%xsqy_sub(nwave+1,wl,wc,kz,dummy,dummy,jndx)
endif
enddo
call set_initialization( status=.false. )
endif
endif
last_n_tuv_z = -1
lat_loop : &
do j = jts,jte
long_loop : &
do i = its,ite
#ifdef SW_DEBUG
tuv_diags = i == minndx(1) .and. j == minndx(2)
#endif
do n = 1,nj
tuv_prate(:,n) = 0.
end do
zenith = zen_angle(i,j)
!---------------------------------------------------------------------
! if night, skip radiative field calculation
!---------------------------------------------------------------------
has_daylight : &
if( zenith < 90. ) then
do k = kts,kte
rad_fld(:,k) = 0.
end do
do wn = 1,nwave
e_fld(:,wn) = 0.
end do
wrk(1) = m2km*(z(i,kte,j) - z_at_w(i,kts,j))
call setzgrid( wrk(1), n_exo_z, zexo_grd )
n_tuv_z = kte + n_exo_z
nlev = n_tuv_z - ktsm1 + 1
nlyr = nlev - 1
do_alloc = n_tuv_z /= last_n_tuv_z
last_n_tuv_z = n_tuv_z
!---------------------------------------------------------------------
! allocate column variables
!---------------------------------------------------------------------
if( do_alloc ) then
call tuv_deallocate
call tuv_allocate
endif
!---------------------------------------------------------------------
! column vertical grid including exo-model top
!---------------------------------------------------------------------
tuv_z(ktsm1) = 0.
tuv_z(kts:kte) = (z(i,kts:kte,j) - z_at_w(i,kts,j))*m2km
tuv_z(kte+1:kte+n_exo_z) = zexo_grd(1:n_exo_z)
! tuv_z(kte+1:kte+n_exo_z) = real( (/ (n,n=k,k+n_exo_z-1) /) )
! tuv_z(kte+n_exo_z+1:n_tuv_z) = real( (/ (25+5*n,n=1,5) /) )
!---------------------------------------------------------------------
! cloud fraction
!---------------------------------------------------------------------
cldfrac(:) = 0.
if( config_flags%cld_od_opt > 1 ) then
if( config_flags%pht_cldfrc_opt == 1 ) then
call cldfrac_binary( cldfrac(kts:kte), &
moist(i,kts:kte,j,p_qc), moist(i,kts:kte,j,p_qi), &
moist(i,kts:kte,j,p_qs), kts, kte )
elseif( config_flags%pht_cldfrc_opt == 2 ) then
if( config_flags%cu_diag == 1 ) then
call cldfrac_fractional( cldfrac(kts:kte), &
moist(i,kts:kte,j,p_qv), &
moist(i,kts:kte,j,p_qc) + gd_cloud(i,kts:kte,j), &
moist(i,kts:kte,j,p_qi) + gd_cloud2(i,kts:kte,j), &
moist(i,kts:kte,j,p_qs), &
p_phy(i,kts:kte,j), t_phy(i,kts:kte,j), kts, kte )
else
call cldfrac_fractional( cldfrac(kts:kte), &
moist(i,kts:kte,j,p_qv), moist(i,kts:kte,j,p_qc), &
moist(i,kts:kte,j,p_qi), moist(i,kts:kte,j,p_qs), &
p_phy(i,kts:kte,j), t_phy(i,kts:kte,j), kts, kte )
endif
endif
endif
!---------------------------------------------------------------------
! aerosols
!---------------------------------------------------------------------
tauaer300(:) = 0.
tauaer400(:) = 0.
tauaer600(:) = 0.
tauaer999(:) = 0.
waer300(:) = 0.
waer400(:) = 0.
waer600(:) = 0.
waer999(:) = 0.
gaer300(:) = 0.
gaer400(:) = 0.
gaer600(:) = 0.
gaer999(:) = 0.
if( has_aer_ra_feedback ) then
tauaer300(kts:kte) = tauaer1(i,kts:kte,j)
tauaer400(kts:kte) = tauaer2(i,kts:kte,j)
tauaer600(kts:kte) = tauaer3(i,kts:kte,j)
tauaer999(kts:kte) = tauaer4(i,kts:kte,j)
waer300(kts:kte) = waer1(i,kts:kte,j)
waer400(kts:kte) = waer2(i,kts:kte,j)
waer600(kts:kte) = waer3(i,kts:kte,j)
waer999(kts:kte) = waer4(i,kts:kte,j)
gaer300(kts:kte) = gaer1(i,kts:kte,j)
gaer400(kts:kte) = gaer2(i,kts:kte,j)
gaer600(kts:kte) = gaer3(i,kts:kte,j)
gaer999(kts:kte) = gaer4(i,kts:kte,j)
tauaer300(ktsm1) = tauaer300(kts)
tauaer400(ktsm1) = tauaer400(kts)
tauaer600(ktsm1) = tauaer600(kts)
tauaer999(ktsm1) = tauaer999(kts)
waer300(ktsm1) = waer300(kts)
waer400(ktsm1) = waer400(kts)
waer600(ktsm1) = waer600(kts)
waer999(ktsm1) = waer999(kts)
gaer300(ktsm1) = gaer300(kts)
gaer400(ktsm1) = gaer400(kts)
gaer600(ktsm1) = gaer600(kts)
gaer999(ktsm1) = gaer999(kts)
endif
tlev(ktsm1) = t8w(i,kts,j)
tlev(kts:kte) = t_phy(i,kts:kte,j)
call z_interp( z_temp_data, temp_data, n_temp_data, &
tuv_z(kte+1:n_tuv_z), tlev(kte+1:n_tuv_z) )
tlyr(ktsm1:n_tuv_z-1) = .5*(tlev(ktsm1:n_tuv_z-1) + tlev(kts:n_tuv_z))
rhoa(ktsm1) = p8w(i,kts,j)/(t8w(i,kts,j)*r_d)
rhoa(kts:kte) = rho_phy(i,kts:kte,j)
dens_air(ktsm1:kte) = xair_kg*rhoa(ktsm1:kte) ! air num.(molecules/cm^3)
call z_interp( z_air_dens_data, air_dens_data, n_air_dens_data, &
tuv_z(kte+1:n_tuv_z), dens_air(kte+1:n_tuv_z) )
o33(kts:kte) = ppm2vmr*chem(i,kts:kte,j,p_o3)*dens_air(kts:kte)
o33(ktsm1) = o33(kts)
call z_interp( z_o3_data, o3_data, n_o3_data, &
tuv_z(kte+1:n_tuv_z), o33(kte+1:n_tuv_z) )
qll(:) = 0.
qll(kts:kte) = moist(i,kts:kte,j,p_qc) + moist(i,kts:kte,j,p_qi)
if( config_flags%cu_diag == 1 ) then
qll(kts:kte) = qll(kts:kte) + gd_cloud(i,kts:kte,j) + gd_cloud2(i,kts:kte,j)
endif
qll(kts:kte) = 1.e3*rhoa(kts:kte)*qll(kts:kte)
where( qll(kts:kte) < 1.e-5 )
qll(kts:kte) = 0.
endwhere
if( haveaer .and. ktau > 1 )then
aerext(ktsm1) = aerext(kts)
else
aerext(ktsm1) = aerwrf(i,kts,j)
endif
if( .not. is_full_tuv ) then
call xs_int( o3_xs, tlyr, o3_xs_tab )
call xs_int( no2_xs, tlyr, no2_xs_tab )
else
call o3xs( nlyr,tlyr,nwave+1,wl,o3_xs_tpose )
call no2xs_jpl06a( nlyr,tlyr,nwave+1,wl,no2_xs_tpose )
o3_xs = transpose( o3_xs_tpose )
no2_xs = transpose( no2_xs_tpose )
endif
dtuv_z(ktsm1:n_tuv_z-1) = tuv_z(kts:n_tuv_z) - tuv_z(ktsm1:n_tuv_z-1)
aircol(ktsm1:n_tuv_z-1) = &
km2cm*dtuv_z(ktsm1:n_tuv_z-1) &
*real(sqrt(real(dens_air(kts:n_tuv_z),kind=8)*real(dens_air(ktsm1:n_tuv_z-1),kind=8)),kind=4)
! km2cm*dtuv_z(ktsm1:n_tuv_z-1)*(dens_air(kts:n_tuv_z) - dens_air(ktsm1:n_tuv_z-1)) &
! / log( dens_air(kts:n_tuv_z)/dens_air(ktsm1:n_tuv_z-1) )
o3col(ktsm1:n_tuv_z-1) = conv1*dtuv_z(ktsm1:n_tuv_z-1)*(o33(kts:n_tuv_z) + o33(ktsm1:n_tuv_z-1))
o2col(ktsm1:n_tuv_z-1) = o2vmr*aircol(ktsm1:n_tuv_z-1)
if( config_flags%has_o3_exo_coldens ) then
o3col(n_tuv_z-1) = o3col(n_tuv_z-1) + real( o3_exo_col(i,j),4 )
o2col(n_tuv_z-1) = o2col(n_tuv_z-1) + real( o2_exo_col(i,j),4 )
aircol(n_tuv_z-1) = aircol(n_tuv_z-1) + o2col(n_tuv_z-1)/o2vmr
endif
if( config_flags%scale_o3_to_grnd_exo_coldens ) then
dobsi = real( o3_exo_col_at_grnd(i,j),4 )/2.687e16
endif
so2col(:) = 0.
if( has_so2 ) then
so2col(ktsm1) = conv2*dtuv_z(ktsm1) &
*chem(i,kts,j,p_so2)*(dens_air(ktsm1) + dens_air(kts))
so2col(kts:ktem1) = conv2*dtuv_z(kts:ktem1) &
*(chem(i,kts:ktem1,j,p_so2)*dens_air(kts:ktem1) &
+ chem(i,kts+1:kte,j,p_so2)*dens_air(kts+1:kte))
endif
no2col(:) = 0.
if( has_no2 ) then
no2col(ktsm1) = conv2*dtuv_z(ktsm1) &
*chem(i,kts,j,p_no2)*(dens_air(ktsm1) + dens_air(kts))
no2col(kts:ktem1) = conv2*dtuv_z(kts:ktem1) &
*(chem(i,kts:ktem1,j,p_no2)*dens_air(kts:ktem1) &
+ chem(i,kts+1:kte,j,p_no2)*dens_air(kts+1:kte))
endif
call tuv_radfld( nlambda_start, config_flags%cld_od_opt, cldfrac, &
nlyr, nwave, zenith, tuv_z, albedo, &
aircol, o2col, o3col, so2col, no2col, &
tauaer300, tauaer400, tauaer600, tauaer999, &
waer300, waer400, waer600, waer999, &
gaer300, gaer400, gaer600, gaer999, &
dtaer, omaer, gaer, dtcld, omcld, gcld, &
has_aer_ra_feedback, &
qll, dobsi, o3_xs, no2_xs, o2_xs, &
so2_xs, wl(1), wc, tlev, srb_o2_xs, rad_fld, e_fld, &
e_dir, e_dn, e_up, &
dir_fld, dwn_fld, up_fld, dtcld_col )
do k = kts,kte
af_dir(i,k,j) = dot_product( dir_fld(k,nlambda_af_start:nlambda_af_end),dw(nlambda_af_start:nlambda_af_end) )
af_dn(i,k,j) = dot_product( dwn_fld(k,nlambda_af_start:nlambda_af_end),dw(nlambda_af_start:nlambda_af_end) )
af_up(i,k,j) = dot_product( up_fld(k,nlambda_af_start:nlambda_af_end),dw(nlambda_af_start:nlambda_af_end) )
end do
dt_cld(i,kts:kte,j) = dtcld_col(kts:kte)
wrk(nlambda_start:nwave) = dw(nlambda_start:nwave)*etfl(nlambda_start:nwave)
wrk_lam(nlambda_start:nwave) = wrk(nlambda_start:nwave)*par_wght(nlambda_start:nwave)
par(i,kts:kte,j) = matmul( e_fld(kts:kte,nlambda_start:nwave), wrk_lam(nlambda_start:nwave) )
wrk_lam(nlambda_start:nwave) = wrk(nlambda_start:nwave)*ery_wght(nlambda_start:nwave)
erythema(i,kts:kte,j) = matmul( e_fld(kts:kte,nlambda_start:nwave), wrk_lam(nlambda_start:nwave) )
if( .not. is_full_tuv ) then
if( any( .not. xsqy_is_zdep(:) ) ) then
rad_fld_tpose = transpose( rad_fld )
endif
elseif( any( xsqy_table(1:nj)%tpflag == 0 ) ) then
rad_fld_tpose = transpose( rad_fld )
endif
#ifdef SW_DEBUG
if( tuv_diags ) then
Dobson(1) = sum( o3col(ktsm1:ktem1) )/2.687e16
Dobson(2) = o3col(kte)/2.687e16
write(*,'(''tuv_drvr: o3col = '',1p,2g15.8)') Dobson(:)
open(unit=33,file='wrfchm_inp')
write(33,*) nlev
write(33,*) kte
write(33,*) tuv_z(ktsm1:n_tuv_z)
write(33,*) tlev(ktsm1:n_tuv_z)
write(33,*) o3col(ktsm1:n_tuv_z-1)
write(33,*) so2col(ktsm1:n_tuv_z-1)
write(33,*) no2col(ktsm1:n_tuv_z-1)
write(33,*) dens_air(ktsm1:n_tuv_z)
write(33,*) aircol(ktsm1:n_tuv_z-1)
do wn = 1,nwave
write(33,*) dtaer(ktsm1:n_tuv_z-1,wn)
end do
do wn = 1,nwave
write(33,*) omaer(ktsm1:n_tuv_z-1,wn)
end do
do wn = 1,nwave
write(33,*) gaer(ktsm1:n_tuv_z-1,wn)
end do
do wn = 1,nwave
write(33,*) dtcld(ktsm1:n_tuv_z-1,wn)
end do
do wn = 1,nwave
write(33,*) omcld(ktsm1:n_tuv_z-1,wn)
end do
do wn = 1,nwave
write(33,*) gcld(ktsm1:n_tuv_z-1,wn)
end do
write(33,*) zen_angle(i,j)
write(33,*) esfact
write(33,*) dobsi
close( 33 )
open(unit=33,file='WRF-TUV.flx.out' )
do n = nlambda_start,nwave
write(33,*) dir_fld(kts:kte,n)
end do
write(33,*) ' '
do n = nlambda_start,nwave
write(33,*) dwn_fld(kts:kte,n)
end do
write(33,*) ' '
do n = nlambda_start,nwave
write(33,*) up_fld(kts:kte,n)
end do
write(33,*) ' '
do n = nlambda_start,nwave
write(33,*) rad_fld_tpose(kts:kte,n)
end do
close( 33 )
endif
#endif
rate_loop: &
do n = 1,nj
!---------------------------------------------------------------------
! set cross-section x quantum yields
!---------------------------------------------------------------------
if( n /= j_o2_ndx ) then
if( .not. is_full_tuv ) then
if( xsqy_is_zdep(n) ) then
call xsqy_int( n, xsqy, tlev(kts:kte), dens_air(kts:kte) )
endif
else
jndx = rxn_ndx(n)
if( jndx /= -1 ) then
if( xsqy_table(jndx)%tpflag /= 0 ) then
call the_subs(jndx)%xsqy_sub(nwave+1,wl,wc,nlev,tlev,dens_air,jndx)
endif
endif
endif
elseif( .not. is_full_tuv ) then
call sjo2( kte, nwave, srb_o2_xs, xsqy )
endif
!---------------------------------------------------------------------
! compute tuv photorates
!---------------------------------------------------------------------
if( .not. is_full_tuv ) then
if( xsqy_is_zdep(n) ) then
do k = kts,kte
xsect(nlambda_start:nwave) = xsqy(nlambda_start:nwave,k)*w_fac(nlambda_start:nwave)*esfact
tuv_prate(k,n) = m2s*dot_product( rad_fld(nlambda_start:nwave,k),xsect(nlambda_start:nwave) )
end do
else
xsect(nlambda_start:nwave) = xsqy_tab(nlambda_start:nwave,1,1,n)*w_fac(nlambda_start:nwave)*esfact
tuv_prate(:,n) = m2s*matmul( rad_fld_tpose(:,nlambda_start:nwave),xsect(nlambda_start:nwave) )
endif
else
if( n /= j_o2_ndx ) then
if( xsqy_table(jndx)%tpflag > 0 ) then
do k = kts,kte
xsect(nlambda_start:nwave) = xsqy_table(jndx)%sq(k+1,nlambda_start:nwave)*w_fac(nlambda_start:nwave)*esfact
tuv_prate(k,n) = m2s*dot_product( rad_fld(nlambda_start:nwave,k),xsect(nlambda_start:nwave) )
end do
else
xsect(nlambda_start:nwave) = xsqy_table(jndx)%sq(nlambda_start:nwave,1)*w_fac(nlambda_start:nwave)*esfact
tuv_prate(:,n) = m2s*matmul( rad_fld_tpose(:,nlambda_start:nwave),xsect(nlambda_start:nwave) )
endif
else
do k = kts,kte
xsect(nlambda_start:nwave) = srb_o2_xs(nlambda_start:nwave,k)*w_fac(nlambda_start:nwave)*esfact
tuv_prate(k,n) = m2s*dot_product( rad_fld(nlambda_start:nwave,k),xsect(nlambda_start:nwave) )
end do
endif
endif
end do rate_loop
#ifdef SW_DEBUG
if( myproc == 6 ) then
if( tuv_diags ) then
open(unit=33,file='WRF-TUV.j.out' )
do n = 1,nj
select case( n )
case( 2:4,7:10,12:13 )
do k = kts,kte,5
write(33,'(1p,5g15.7)') tuv_prate(k:min(k+4,kte),n)/m2s
end do
if( n < 13 ) then
write(33,*) ' '
endif
end select
end do
close(33)
call wrf_abort
endif
endif
#endif
endif has_daylight
#if ( WRF_KPP == 1 )
#include "tuv2wrf_jvals.inc"
#endif
end do long_loop
end do lat_loop
#ifdef SW_DEBUG
if( any( zen_angle(its:ite,jts:jte) < 90. ) ) then
min_ndx(:) = minloc( z_top,mask=z_top>0. )
min_ndx(:) = (/ min_ndx(1) + its - 1,min_ndx(2) + jts - 1 /)
max_ndx(:) = maxloc( z_top,mask=z_top>0. )
max_ndx(:) = (/ max_ndx(1) + its - 1,max_ndx(2) + jts - 1 /)
write(msg,'(''tuv_driver: z_top min indices = '',2i6)') min_ndx(:)
call wrf_debug( 0,trim(msg) )
write(msg,'(''tuv_driver: z_top max indices = '',2i6)') max_ndx(:)
call wrf_debug( 0,trim(msg) )
write(msg,'(''tuv_driver: z_top min,max = '',1p2g15.7)') z_top(min_ndx(1),min_ndx(2)),z_top(max_ndx(1),max_ndx(2))
call wrf_debug( 0,trim(msg) )
endif
#endif
call tuv_deallocate
ierr = 0
if( allocated( rad_fld ) ) then
deallocate( rad_fld,stat=astat )
ierr = ierr + astat
deallocate( dir_fld, dwn_fld, up_fld,stat=astat )
ierr = ierr + astat
deallocate( e_dir, e_dn, e_up,stat=astat )
ierr = ierr + astat
endif
if( allocated( rad_fld_tpose ) ) then
deallocate( rad_fld_tpose,stat=astat )
ierr = ierr + astat
endif
if( allocated( e_fld ) ) then
deallocate( e_fld,stat=astat )
ierr = ierr + astat
endif
if( allocated( xsqy ) ) then
deallocate( xsqy,stat=astat )
ierr = ierr + astat
endif
if( allocated( srb_o2_xs ) ) then
deallocate( srb_o2_xs,stat=astat )
ierr = ierr + astat
endif
if( allocated( tuv_prate ) ) then
deallocate( tuv_prate,stat=astat )
ierr = ierr + astat
endif
if( ierr /= 0 ) then
call wrf_error_fatal( 'tuv_driver: failed to deallocate local variables' )
endif
CONTAINS
subroutine setzgrid( ztop, nexo, zexo_grd )
!---------------------------------------------------------------------
! set the vertical grid above model top up to 50 km for photorate
! calculation
!---------------------------------------------------------------------
integer, intent(out) :: nexo
real, intent(in) :: ztop
real, intent(out) :: zexo_grd(:)
real, parameter :: ztrop = 25.
real, parameter :: dztrop = 1.
real, parameter :: zlid = 50.
real, parameter :: dzlid = 5.
real :: z
real :: zBase
nexo = 0
if( ztop >= zlid ) then
return
elseif( ztop <= ztrop ) then
zBase = ztrop
else
zBase = dzlid*real( int(ztop/dzlid) )
endif
z = int( ztop )
do while( z < ztrop )
z = z + dztrop
nexo = nexo + 1
zexo_grd(nexo) = z
enddo
z = zBase
do while( z < zlid )
z = z + dzlid
nexo = nexo + 1
zexo_grd(nexo) = z
end do
end subroutine setzgrid
subroutine tuv_allocate
!---------------------------------------------------------------------
! allocate column variables
!---------------------------------------------------------------------
allocate( tuv_z(ktsm1:n_tuv_z),dtuv_z(ktsm1:n_tuv_z-1),stat=ierr )
allocate( cldfrac(ktsm1:n_tuv_z),stat=astat )
ierr = ierr + astat
allocate( qll(ktsm1:n_tuv_z),stat=astat )
ierr = ierr + astat
allocate( tlev(ktsm1:n_tuv_z),tlyr(ktsm1:n_tuv_z-1),stat=astat )
ierr = ierr + astat
allocate( dens_air(ktsm1:n_tuv_z),stat=astat )
ierr = ierr + astat
allocate( o33(ktsm1:n_tuv_z),stat=astat )
ierr = ierr + astat
allocate( o3_xs(nwave,ktsm1:n_tuv_z-1),stat=astat )
ierr = ierr + astat
if( is_full_tuv ) then
allocate( o3_xs_tpose(ktsm1:n_tuv_z-1,nwave),stat=astat )
ierr = ierr + astat
allocate( no2_xs_tpose(ktsm1:n_tuv_z-1,nwave),stat=astat )
ierr = ierr + astat
endif
allocate( no2_xs(nwave,ktsm1:n_tuv_z-1),stat=astat )
ierr = ierr + astat
allocate( aircol(ktsm1:n_tuv_z-1),stat=astat )
ierr = ierr + astat
allocate( o2col(ktsm1:n_tuv_z-1),stat=astat )
ierr = ierr + astat
allocate( o3col(ktsm1:n_tuv_z-1),stat=astat )
ierr = ierr + astat
allocate( no2col(ktsm1:n_tuv_z-1),stat=astat )
ierr = ierr + astat
allocate( so2col(ktsm1:n_tuv_z-1),stat=astat )
ierr = ierr + astat
allocate( tauaer300(ktsm1:n_tuv_z-1),tauaer400(ktsm1:n_tuv_z-1), &
tauaer600(ktsm1:n_tuv_z-1),tauaer999(ktsm1:n_tuv_z-1),stat=astat )
ierr = ierr + astat
allocate( waer300(ktsm1:n_tuv_z-1),waer400(ktsm1:n_tuv_z-1), &
waer600(ktsm1:n_tuv_z-1),waer999(ktsm1:n_tuv_z-1),stat=astat )
ierr = ierr + astat
allocate( gaer300(ktsm1:n_tuv_z-1),gaer400(ktsm1:n_tuv_z-1), &
gaer600(ktsm1:n_tuv_z-1),gaer999(ktsm1:n_tuv_z-1),stat=astat )
ierr = ierr + astat
allocate( dtaer(ktsm1:n_tuv_z-1,nwave),omaer(ktsm1:n_tuv_z-1,nwave), &
gaer(ktsm1:n_tuv_z-1,nwave),stat=astat )
ierr = ierr + astat
allocate( dtcld(ktsm1:n_tuv_z-1,nwave),omcld(ktsm1:n_tuv_z-1,nwave), &
gcld(ktsm1:n_tuv_z-1,nwave),stat=astat )
ierr = ierr + astat
if( ierr /= 0 ) then
call wrf_error_fatal( 'tuv_driver: failed to allocate column variables' )
endif
end subroutine tuv_allocate
subroutine tuv_deallocate
!---------------------------------------------------------------------
! deallocate column variables
!---------------------------------------------------------------------
integer :: astat, ierr
ierr = 0
if( allocated( tuv_z ) ) then
deallocate( tuv_z,stat=astat )
ierr = ierr + astat
endif
if( allocated( dtuv_z ) ) then
deallocate( dtuv_z,stat=astat )
ierr = ierr + astat
endif
if( allocated( cldfrac ) ) then
deallocate( cldfrac,stat=astat )
ierr = ierr + astat
endif
if( allocated( qll ) ) then
deallocate( qll,stat=astat )
ierr = ierr + astat
endif
if( allocated( tlev ) ) then
deallocate( tlev,stat=astat )
ierr = ierr + astat
endif
if( allocated( tlyr ) ) then
deallocate( tlyr,stat=astat )
ierr = ierr + astat
endif
if( allocated( dens_air ) ) then
deallocate( dens_air,stat=astat )
ierr = ierr + astat
endif
if( allocated( o33 ) ) then
deallocate( o33,stat=astat )
ierr = ierr + astat
endif
if( allocated( o3_xs ) ) then
deallocate( o3_xs,stat=astat )
ierr = ierr + astat
endif
if( allocated( o3_xs_tpose ) ) then
deallocate( o3_xs_tpose,stat=astat )
ierr = ierr + astat
endif
if( allocated( no2_xs ) ) then
deallocate( no2_xs,stat=astat )
ierr = ierr + astat
endif
if( allocated( no2_xs_tpose ) ) then
deallocate( no2_xs_tpose,stat=astat )
ierr = ierr + astat
endif
if( allocated( aircol ) ) then
deallocate( aircol,stat=astat )
ierr = ierr + astat
endif
if( allocated( o2col ) ) then
deallocate( o2col,stat=astat )
ierr = ierr + astat
endif
if( allocated( o3col ) ) then
deallocate( o3col,stat=astat )
ierr = ierr + astat
endif
if( allocated( no2col ) ) then
deallocate( no2col,stat=astat )
ierr = ierr + astat
endif
if( allocated( so2col ) ) then
deallocate( so2col,stat=astat )
ierr = ierr + astat
endif
if( allocated( tauaer300 ) ) then
deallocate( tauaer300, waer300, gaer300,stat=astat )
ierr = ierr + astat
endif
if( allocated( tauaer400 ) ) then
deallocate( tauaer400, waer400, gaer400,stat=astat )
ierr = ierr + astat
endif
if( allocated( tauaer600 ) ) then
deallocate( tauaer600, waer600, gaer600,stat=astat )
ierr = ierr + astat
endif
if( allocated( tauaer999 ) ) then
deallocate( tauaer999, waer999, gaer999,stat=astat )
ierr = ierr + astat
endif
if( allocated( dtcld ) ) then
deallocate( dtcld, omcld, gcld,stat=astat )
ierr = ierr + astat
endif
if( allocated( dtaer ) ) then
deallocate( dtaer, omaer, gaer,stat=astat )
ierr = ierr + astat
endif
if( ierr /= 0 ) then
call wrf_error_fatal( 'tuv_deallocate: failed to deallocate all variables' )
endif
end subroutine tuv_deallocate
end subroutine tuv_driver
subroutine tuv_init( &
domain, config_flags, z_at_w, aerwrf, g, &
af_lambda_start, af_lambda_end, start_lambda, &
ids,ide, jds,jde, kds,kde, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte )
USE module_configure
USE params_mod, only : lambda_cutoff
USE srb, only : init_srb
USE module_state_description, only : p_so2, p_no2, param_first_scalar
USE module_rxn, only : rxn_init, xsqy_table => xsqy_tab, npht_tab
USE module_rxn, only : get_initialization
! .. Scalar Arguments ..
INTEGER, INTENT (IN) :: domain
INTEGER, INTENT (IN) :: ide, ids, ime, ims, ite, its, jde, jds, &
jme, jms, jte, jts, kde, kds, kme, kms, kte, kts
REAL, INTENT (IN) :: g
REAL, INTENT (IN) :: af_lambda_start, af_lambda_end, start_lambda
! .. Array Arguments ..
REAL, INTENT (INOUT) :: aerwrf(ims:ime,kms:kme,jms:jme)
REAL, INTENT (IN) :: z_at_w(ims:ime,kms:kme,jms:jme)
TYPE(grid_config_rec_type), INTENT(IN) :: config_flags
! .. Local Scalars ..
INTEGER :: astat
INTEGER :: i, j, k, n, wn
! .. Local Arrays ..
CHARACTER(len=256) :: msg
#ifndef NETCDF
call wrf_error_fatal( 'tuv_init: requires netcdf' )
#endif
DO j = jts, min(jte,jde-1)
DO k = kts, min(kte,kde-1)
DO i = its, min(ite,ide-1)
aerwrf(i,k,j) = 0.
END DO
END DO
END DO
is_initialized: &
if( .not. tuv_is_initialized ) then
has_exo_coldens = config_flags%has_o3_exo_coldens .or. config_flags%scale_o3_to_grnd_exo_coldens
is_full_tuv = config_flags%is_full_tuv
#if ( WRF_KPP == 1 )
#include "tuvdef_jvals.inc"
#endif
call get_xsqy_tab
if( .not. is_full_tuv ) then
allocate( xsqy_is_zdep(nj),stat=astat )
if( astat /= 0 ) then
call wrf_error_fatal( 'tuv_init: failed to allocate xsqy_is_zdep' )
endif
xsqy_is_zdep(:) = .false.
if( j_o2_ndx > 0 ) then
xsqy_is_zdep(j_o2_ndx) = .true.
endif
do n = 1,nj
if( n /= j_o2_ndx ) then
t_loop : do j = 1,nconc
do i = 1,ntemp
if( any( xsqy_tab(:,i,j,n) /= xsqy_tab(:,1,1,n) ) ) then
xsqy_is_zdep(n) = .true.
exit t_loop
endif
end do
end do t_loop
endif
end do
endif
has_so2 = p_so2 >= param_first_scalar
has_no2 = p_no2 >= param_first_scalar
call init_srb
!---------------------------------------------------------------------
! ... locate starting wave bin
!---------------------------------------------------------------------
lambda_cutoff = start_lambda
do nlambda_start = 1,nwave
if( wc(nlambda_start) >= lambda_cutoff ) then
exit
endif
end do
if( nlambda_start > nwave ) then
write(msg,'(''tuv_init: '',1pg15.7,'' is not in photo wavelength interval ('',g15.7,'','',g15.7,'')'')') &
lambda_cutoff,wl(1),wl(nwave+1)
call wrf_error_fatal( trim(msg) )
endif
!---------------------------------------------------------------------
! ... locate af output wave bins
!---------------------------------------------------------------------
do nlambda_af_start = nlambda_start,nwave
if( wl(nlambda_af_start) <= af_lambda_start .and. &
wl(nlambda_af_start+1) >= af_lambda_start ) then
exit
endif
end do
do nlambda_af_end = nlambda_start,nwave
if( wl(nlambda_af_end) <= af_lambda_end .and. &
wl(nlambda_af_end+1) >= af_lambda_end ) then
exit
endif
end do
allocate( par_wght(nwave), ery_wght(nwave), stat=astat )
if( astat /= 0 ) then
call wrf_error_fatal( 'tuv_init: failed to allocate par_wght,ery_wght' )
endif
!---------------------------------------------------------------------
! ... set PAR weight
!---------------------------------------------------------------------
where (wc(:) > 400. .AND. wc(:) < 700.)
par_wght(:) = 8.36e-3 * wc(:)
elsewhere
par_wght(:) = 0.
end where
!---------------------------------------------------------------------
! ... set erythema weight
!---------------------------------------------------------------------
call fery( nwave, wc, ery_wght )
!---------------------------------------------------------------------
! ... set surface albedo
!---------------------------------------------------------------------
DO wn = 1,nwave
IF (wl(wn)<400.) then
albedo(wn) = 0.05
elseIF (wl(wn)>=400. .AND. wl(wn)<450.) then
albedo(wn) = 0.06
elseIF (wl(wn)>=450. .AND. wl(wn)<500.) then
albedo(wn) = 0.08
elseIF (wl(wn)>=500. .AND. wl(wn)<550.) then
albedo(wn) = 0.10
elseIF (wl(wn)>=550. .AND. wl(wn)<600.) then
albedo(wn) = 0.11
elseIF (wl(wn)>=600. .AND. wl(wn)<640.) then
albedo(wn) = 0.12
elseIF (wl(wn)>=640. .AND. wl(wn)<660.) then
albedo(wn) = 0.135
elseIF (wl(wn)>=660.) then
albedo(wn) = 0.15
endIF
END DO
!---------------------------------------------------------------------
! ... if full TUV then initialize
!---------------------------------------------------------------------
if( is_full_tuv ) then
call rxn_init( nwave+1,wl )
allocate( rxn_ndx(nj),stat=astat )
if( astat /= 0 ) then
call wrf_error_fatal( 'tuv_init: failed to allocate rxn_ndx' )
endif
rxn_ndx(1:nj) = -1
do j = 1,nj
if( j /= j_o2_ndx ) then
do n = 2,npht_tab
if( trim(xsqy_table(n)%wrf_label) == trim(tuv_jname(j)) ) then
rxn_ndx(j) = n
exit
endif
enddo
endif
enddo
rxn_initialized = .not. get_initialization()
endif
tuv_is_initialized = .true.
endif is_initialized
!---------------------------------------------------------------------
! ... get exo model column densities
!---------------------------------------------------------------------
if( has_exo_coldens ) then
call get_exo_coldens( domain, config_flags%exo_coldens_inname, &
ids,ide, jds,jde, kds,kde, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte )
endif
end subroutine tuv_init
subroutine get_xsqy_tab
!---------------------------------------------------------------------
! ... read in the cross section,quantum yield tables
!---------------------------------------------------------------------
use params_mod, only : hc
use srb, only : ila, isrb
use srb, only : nchebev_term, nchebev_wave
use srb, only : chebev_ac, chebev_bc
!---------------------------------------------------------------------
! ... local variables
!---------------------------------------------------------------------
integer :: astat, ierr
integer :: m
integer :: ncid, dimid, varid
character(len=132) :: filename
character(len=132) :: err_msg
character(len=64) :: varname
#ifdef NETCDF
include 'netcdf.inc'
!---------------------------------------------------------------------
! ... function declarations
!---------------------------------------------------------------------
LOGICAL, EXTERNAL :: wrf_dm_on_monitor
master_proc_a: &
if( wrf_dm_on_monitor() ) then
filename = 'wrf_tuv_xsqy.nc'
err_msg = 'get_xsqy_tab: failed to open file ' // trim(filename)
call handle_ncerr( nf_open( trim(filename), nf_noclobber, ncid ), trim(err_msg) )
!---------------------------------------------------------------------
! ... get dimensions
!---------------------------------------------------------------------
err_msg = 'get_xsqy_tab: failed to get nwave id'
call handle_ncerr( nf_inq_dimid( ncid, 'nwave', dimid ), trim(err_msg) )
err_msg = 'get_xsqy_tab: failed to get nwave'
call handle_ncerr( nf_inq_dimlen( ncid, dimid, nwave ), trim(err_msg) )
err_msg = 'get_xsqy_tab: failed to get ntemp id'
call handle_ncerr( nf_inq_dimid( ncid, 'ntemp', dimid ), trim(err_msg) )
err_msg = 'get_xsqy_tab: failed to get ntemp'
call handle_ncerr( nf_inq_dimlen( ncid, dimid, ntemp ), trim(err_msg) )
err_msg = 'get_xsqy_tab: failed to get nconc id'
call handle_ncerr( nf_inq_dimid( ncid, 'nconc', dimid ), trim(err_msg) )
err_msg = 'get_xsqy_tab: failed to get nconc'
call handle_ncerr( nf_inq_dimlen( ncid, dimid, nconc ), trim(err_msg) )
err_msg = 'get_xsqy_tab: failed to get nchebev_term id'
call handle_ncerr( nf_inq_dimid( ncid, 'nchebev_term', dimid ), trim(err_msg) )
err_msg = 'get_xsqy_tab: failed to get nchebev'
call handle_ncerr( nf_inq_dimlen( ncid, dimid, nchebev_term ), trim(err_msg) )
err_msg = 'get_xsqy_tab: failed to get nchebev_wave id'
call handle_ncerr( nf_inq_dimid( ncid, 'nchebev_wave', dimid ), trim(err_msg) )
err_msg = 'get_xsqy_tab: failed to get nchebev'
call handle_ncerr( nf_inq_dimlen( ncid, dimid, nchebev_wave ), trim(err_msg) )
err_msg = 'get_xsqy_tab: failed to get n_temp_data id'
call handle_ncerr( nf_inq_dimid( ncid, 'n_temp_data', dimid ), trim(err_msg) )
err_msg = 'get_xsqy_tab: failed to get n_temp_data'
call handle_ncerr( nf_inq_dimlen( ncid, dimid, n_temp_data ), trim(err_msg) )
err_msg = 'get_xsqy_tab: failed to get n_o3_data id'
call handle_ncerr( nf_inq_dimid( ncid, 'n_o3_data', dimid ), trim(err_msg) )
err_msg = 'get_xsqy_tab: failed to get n_o3_data'
call handle_ncerr( nf_inq_dimlen( ncid, dimid, n_o3_data ), trim(err_msg) )
err_msg = 'get_xsqy_tab: failed to get n_air_dens_data id'
call handle_ncerr( nf_inq_dimid( ncid, 'n_air_dens_data', dimid ), trim(err_msg) )
err_msg = 'get_xsqy_tab: failed to get n_air_dens_data'
call handle_ncerr( nf_inq_dimlen( ncid, dimid, n_air_dens_data ), trim(err_msg) )
endif master_proc_a
#ifdef DM_PARALLEL
!---------------------------------------------------------------------
! ... bcast the dimensions
!---------------------------------------------------------------------
CALL wrf_dm_bcast_bytes ( nwave, IWORDSIZE )
CALL wrf_dm_bcast_bytes ( ntemp, IWORDSIZE )
CALL wrf_dm_bcast_bytes ( nconc, IWORDSIZE )
CALL wrf_dm_bcast_bytes ( n_temp_data, IWORDSIZE )
CALL wrf_dm_bcast_bytes ( n_o3_data, IWORDSIZE )
CALL wrf_dm_bcast_bytes ( n_air_dens_data, IWORDSIZE )
CALL wrf_dm_bcast_bytes ( nchebev_term, IWORDSIZE )
CALL wrf_dm_bcast_bytes ( nchebev_wave, IWORDSIZE )
#endif
!---------------------------------------------------------------------
! ... allocate module arrays
!---------------------------------------------------------------------
ierr = 0
allocate( z_temp_data(n_temp_data), z_o3_data(n_o3_data), &
z_air_dens_data(n_air_dens_data),stat=astat )
ierr = astat + ierr
allocate( temp_data(n_temp_data), o3_data(n_o3_data), &
air_dens_data(n_air_dens_data),stat=astat )
ierr = astat + ierr
allocate( wl(nwave+1), wc(nwave), dw(nwave), w_fac(nwave), &
etfl(nwave), albedo(nwave), stat=astat )
ierr = astat + ierr
if( .not. is_full_tuv ) then
allocate( temp_tab(ntemp), conc_tab(nconc), stat=astat )
ierr = astat + ierr
allocate( del_temp_tab(ntemp-1), del_conc_tab(nconc-1), stat=astat )
ierr = astat + ierr
endif
allocate( chebev_ac(nchebev_term,nchebev_wave), stat=astat )
ierr = astat + ierr
allocate( chebev_bc(nchebev_term,nchebev_wave), stat=astat )
ierr = astat + ierr
allocate( o2_xs(nwave), so2_xs(nwave), stat=astat )
ierr = astat + ierr
allocate( o3_xs_tab(nwave,ntemp), no2_xs_tab(nwave,ntemp), stat=astat )
ierr = astat + ierr
if( .not. is_full_tuv ) then
allocate( xsqy_tab(nwave,ntemp,nconc,nj), stat=astat )
ierr = astat + ierr
endif
if( ierr /= 0 ) then
call wrf_error_fatal( 'tuv_init: failed to allocate z_temp_data ... xsqy_tab' )
endif
!---------------------------------------------------------------------
! ... read arrays
!---------------------------------------------------------------------
master_proc_b: &
if( wrf_dm_on_monitor() ) then
err_msg = 'get_xsqy_tab: failed to get z_temp_data variable id'
call handle_ncerr( nf_inq_varid( ncid, 'z_temp_data', varid ), trim(err_msg) )
err_msg = 'get_xsqy_tab: failed to read z_temp_data variable'
call handle_ncerr( nf_get_var_real( ncid, varid, z_temp_data ), trim(err_msg) )
err_msg = 'get_xsqy_tab: failed to get z_o3_data variable id'
call handle_ncerr( nf_inq_varid( ncid, 'z_o3_data', varid ), trim(err_msg) )
err_msg = 'get_xsqy_tab: failed to read z_o3_data variable'
call handle_ncerr( nf_get_var_real( ncid, varid, z_o3_data ), trim(err_msg) )
err_msg = 'get_xsqy_tab: failed to get z_air_dens_data variable id'
call handle_ncerr( nf_inq_varid( ncid, 'z_air_dens_data', varid ), trim(err_msg) )
err_msg = 'get_xsqy_tab: failed to read z_air_dens_data variable'
call handle_ncerr( nf_get_var_real( ncid, varid, z_air_dens_data ), trim(err_msg) )
err_msg = 'get_xsqy_tab: failed to get temp_data variable id'
call handle_ncerr( nf_inq_varid( ncid, 'temp_data', varid ), trim(err_msg) )
err_msg = 'get_xsqy_tab: failed to read temp_data variable'
call handle_ncerr( nf_get_var_real( ncid, varid, temp_data ), trim(err_msg) )
err_msg = 'get_xsqy_tab: failed to get o3_data variable id'
call handle_ncerr( nf_inq_varid( ncid, 'o3_data', varid ), trim(err_msg) )
err_msg = 'get_xsqy_tab: failed to read o3_data variable'
call handle_ncerr( nf_get_var_real( ncid, varid, o3_data ), trim(err_msg) )
err_msg = 'get_xsqy_tab: failed to get air_dens_data variable id'
call handle_ncerr( nf_inq_varid( ncid, 'air_dens_data', varid ), trim(err_msg) )
err_msg = 'get_xsqy_tab: failed to read air_dens_data variable'
call handle_ncerr( nf_get_var_real( ncid, varid, air_dens_data ), trim(err_msg) )
err_msg = 'get_xsqy_tab: failed to get wl variable id'
call handle_ncerr( nf_inq_varid( ncid, 'wl', varid ), trim(err_msg) )
err_msg = 'get_xsqy_tab: failed to read wl variable'
call handle_ncerr( nf_get_var_real( ncid, varid, wl ), trim(err_msg) )
err_msg = 'get_xsqy_tab: failed to get wc variable id'
call handle_ncerr( nf_inq_varid( ncid, 'wc', varid ), trim(err_msg) )
err_msg = 'get_xsqy_tab: failed to read wc variable'
call handle_ncerr( nf_get_var_real( ncid, varid, wc ), trim(err_msg) )
err_msg = 'get_xsqy_tab: failed to get etfl variable id'
call handle_ncerr( nf_inq_varid( ncid, 'etf', varid ), trim(err_msg) )
err_msg = 'get_xsqy_tab: failed to read etfl variable'
call handle_ncerr( nf_get_var_real( ncid, varid, etfl ), trim(err_msg) )
err_msg = 'get_xsqy_tab: failed to get chebev_ac variable id'
call handle_ncerr( nf_inq_varid( ncid, 'chebev_ac', varid ), trim(err_msg) )
err_msg = 'get_xsqy_tab: failed to read chebev_ac variable'
call handle_ncerr( nf_get_var_double( ncid, varid, chebev_ac ), trim(err_msg) )
err_msg = 'get_xsqy_tab: failed to get chebev_bc variable id'
call handle_ncerr( nf_inq_varid( ncid, 'chebev_bc', varid ), trim(err_msg) )
err_msg = 'get_xsqy_tab: failed to read chebev_bc variable'
call handle_ncerr( nf_get_var_double( ncid, varid, chebev_bc ), trim(err_msg) )
err_msg = 'get_xsqy_tab: failed to get ila variable id'
call handle_ncerr( nf_inq_varid( ncid, 'ila', varid ), trim(err_msg) )
err_msg = 'get_xsqy_tab: failed to read ila variable'
call handle_ncerr( nf_get_var_int( ncid, varid, ila ), trim(err_msg) )
err_msg = 'get_xsqy_tab: failed to get isrb variable id'
call handle_ncerr( nf_inq_varid( ncid, 'isrb', varid ), trim(err_msg) )
err_msg = 'get_xsqy_tab: failed to read isrb variable'
call handle_ncerr( nf_get_var_int( ncid, varid, isrb ), trim(err_msg) )
if( .not. is_full_tuv ) then
err_msg = 'get_xsqy_tab: failed to temp_tab variable id'
call handle_ncerr( nf_inq_varid( ncid, 'temps', varid ), trim(err_msg) )
err_msg = 'get_xsqy_tab: failed to read temp_tab variable'
call handle_ncerr( nf_get_var_real( ncid, varid, temp_tab ), trim(err_msg) )
err_msg = 'get_xsqy_tab: failed to conc_tab variable id'
call handle_ncerr( nf_inq_varid( ncid, 'concs', varid ), trim(err_msg) )
err_msg = 'get_xsqy_tab: failed to read conc_tab variable'
call handle_ncerr( nf_get_var_real( ncid, varid, conc_tab ), trim(err_msg) )
endif
err_msg = 'get_xsqy_tab: failed to o2_xs variable id'
call handle_ncerr( nf_inq_varid( ncid, 'o2_xs', varid ), trim(err_msg) )
err_msg = 'get_xsqy_tab: failed to read o2_xs variable'
call handle_ncerr( nf_get_var_real( ncid, varid, o2_xs ), trim(err_msg) )
err_msg = 'get_xsqy_tab: failed to so2_xs variable id'
call handle_ncerr( nf_inq_varid( ncid, 'so2_xs', varid ), trim(err_msg) )
err_msg = 'get_xsqy_tab: failed to read so2_xs variable'
call handle_ncerr( nf_get_var_real( ncid, varid, so2_xs ), trim(err_msg) )
err_msg = 'get_xsqy_tab: failed to o3_xs_tab variable id'
call handle_ncerr( nf_inq_varid( ncid, 'o3_xs', varid ), trim(err_msg) )
err_msg = 'get_xsqy_tab: failed to read o3_xs_tab variable'
call handle_ncerr( nf_get_var_real( ncid, varid, o3_xs_tab ), trim(err_msg) )
err_msg = 'get_xsqy_tab: failed to no2_xs_tab variable id'
call handle_ncerr( nf_inq_varid( ncid, 'no2_xs', varid ), trim(err_msg) )
err_msg = 'get_xsqy_tab: failed to read no2_xs_tab variable'
call handle_ncerr( nf_get_var_real( ncid, varid, no2_xs_tab ), trim(err_msg) )
if( .not. is_full_tuv ) then
do m = 1,nj
varname = trim(tuv_jname(m)) // '_xsqy'
err_msg = 'get_xsqy_tab: failed to ' // trim(varname) //' variable id'
call handle_ncerr( nf_inq_varid( ncid, trim(varname), varid ), trim(err_msg) )
err_msg = 'get_xsqy_tab: failed to read ' // trim(varname) // ' variable'
call handle_ncerr( nf_get_var_real( ncid, varid, xsqy_tab(:,:,:,m) ), trim(err_msg) )
end do
endif
endif master_proc_b
#ifdef DM_PARALLEL
!---------------------------------------------------------------------
! ... bcast the arrays
!---------------------------------------------------------------------
CALL wrf_dm_bcast_bytes ( ila, IWORDSIZE )
CALL wrf_dm_bcast_bytes ( isrb, IWORDSIZE )
CALL wrf_dm_bcast_bytes( wl, (nwave+1)*RWORDSIZE )
CALL wrf_dm_bcast_bytes( wc, nwave*RWORDSIZE )
CALL wrf_dm_bcast_bytes( etfl, nwave*RWORDSIZE )
if( .not. is_full_tuv ) then
CALL wrf_dm_bcast_bytes( temp_tab, ntemp*RWORDSIZE )
CALL wrf_dm_bcast_bytes( conc_tab, nconc*RWORDSIZE )
endif
CALL wrf_dm_bcast_bytes( o2_xs, nwave*RWORDSIZE )
CALL wrf_dm_bcast_bytes( so2_xs, nwave*RWORDSIZE )
CALL wrf_dm_bcast_bytes( z_temp_data, n_temp_data*RWORDSIZE )
CALL wrf_dm_bcast_bytes( z_o3_data, n_o3_data*RWORDSIZE )
CALL wrf_dm_bcast_bytes( z_air_dens_data, n_air_dens_data*RWORDSIZE )
CALL wrf_dm_bcast_bytes( temp_data, n_temp_data*RWORDSIZE )
CALL wrf_dm_bcast_bytes( o3_data, n_o3_data*RWORDSIZE )
CALL wrf_dm_bcast_bytes( air_dens_data, n_air_dens_data*RWORDSIZE )
#if RWORDSIZE == 4
CALL wrf_dm_bcast_bytes( chebev_ac, nchebev_term*nchebev_wave*2*RWORDSIZE )
CALL wrf_dm_bcast_bytes( chebev_bc, nchebev_term*nchebev_wave*2*RWORDSIZE )
#endif
#if RWORDSIZE == 8
CALL wrf_dm_bcast_bytes( chebev_ac, nchebev_term*nchebev_wave*RWORDSIZE )
CALL wrf_dm_bcast_bytes( chebev_bc, nchebev_term*nchebev_wave*RWORDSIZE )
#endif
CALL wrf_dm_bcast_bytes( o3_xs_tab, nwave*ntemp*RWORDSIZE )
CALL wrf_dm_bcast_bytes( no2_xs_tab, nwave*ntemp*RWORDSIZE )
if( .not. is_full_tuv ) then
CALL wrf_dm_bcast_bytes( xsqy_tab, nwave*ntemp*nconc*nj*RWORDSIZE )
endif
#endif
if( .not. is_full_tuv ) then
del_temp_tab(:ntemp-1) = 1./(temp_tab(2:ntemp) - temp_tab(1:ntemp-1))
del_conc_tab(:nconc-1) = 1./(conc_tab(2:nconc) - conc_tab(1:nconc-1))
endif
dw(:nwave) = wl(2:nwave+1) - wl(1:nwave)
w_fac(:nwave) = dw(:nwave)*etfl(:nwave)*1.e-13*wc(:nwave)/hc
if( wrf_dm_on_monitor() ) then
err_msg = 'get_xsqy_tab: failed to close file ' // trim(filename)
call handle_ncerr( nf_close( ncid ),trim(err_msg) )
endif
#endif
end subroutine get_xsqy_tab
subroutine get_exo_coldens( dm, exo_coldens_filename, &
ids,ide, jds,jde, kds,kde, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte )
!---------------------------------------------------------------------
! ... read in the exo column o2,o3 densities
!---------------------------------------------------------------------
!---------------------------------------------------------------------
! ... dummy arguments
!---------------------------------------------------------------------
integer, intent(in) :: dm
integer, intent(in) :: ide, ids, ime, ims, ite, its, jde, jds, &
jme, jms, jte, jts, kde, kds, kme, kms, kte, kts
character(len=*), intent(in) :: exo_coldens_filename
!---------------------------------------------------------------------
! ... local variables
!---------------------------------------------------------------------
INTEGER :: i, j, k
integer :: astat
integer :: ncid
integer :: dimid
integer :: varid
integer :: max_dom
integer :: cpos
integer :: iend, jend
integer :: lon_e, lat_e
integer :: ncoldens_levs
integer :: ndays_of_year
! real, allocatable :: coldens(:,:,:,:)
character(len=128) :: err_msg
character(len=64) :: filename
character(len=2) :: id_num
!---------------------------------------------------------------------
! ... function declarations
!---------------------------------------------------------------------
LOGICAL, EXTERNAL :: wrf_dm_on_monitor
#ifdef NETCDF
include 'netcdf.inc'
#define DM_BCAST_MACRO(A) CALL wrf_dm_bcast_bytes ( A , size ( A ) * DWORDSIZE )
!---------------------------------------------------------------------
! ... allocate column_density type
!---------------------------------------------------------------------
if( .not. allocated(col_dens) ) then
CALL nl_get_max_dom( 1,max_dom )
allocate( col_dens(max_dom),stat=astat )
if( astat /= 0 ) then
call wrf_error_fatal( 'get_exo_coldens: failed to allocate col_dens' )
endif
write(err_msg,'(''get_exo_coldens: intializing '',i2,'' domains'')') max_dom
call wrf_message( trim(err_msg) )
col_dens(:)%is_allocated = .false.
endif
!---------------------------------------------------------------------
! ... open column density netcdf file
!---------------------------------------------------------------------
col_dens_allocated : &
if( .not. col_dens(dm)%is_allocated ) then
! if( wrf_dm_on_monitor() ) then
cpos = index( exo_coldens_filename, '<domain>' )
if( cpos > 0 ) then
write(id_num,'(i2.2)') dm
filename = exo_coldens_filename(:cpos-1) // 'd' // id_num
else
filename = trim( exo_coldens_filename )
endif
err_msg = 'get_exo_coldens: intializing domain ' // id_num
call wrf_message( trim(err_msg) )
err_msg = 'get_exo_coldens: failed to open file ' // trim(filename)
call handle_ncerr( nf_open( trim(filename), nf_noclobber, ncid ), trim(err_msg) )
!---------------------------------------------------------------------
! ... get dimensions
!---------------------------------------------------------------------
err_msg = 'get_exo_coldens: failed to get col_dens levels id'
call handle_ncerr( nf_inq_dimid( ncid, 'coldens_levs', dimid ), trim(err_msg) )
err_msg = 'get_exo_coldens: failed to get col_dens levels'
call handle_ncerr( nf_inq_dimlen( ncid, dimid, ncoldens_levs ), trim(err_msg) )
err_msg = 'get_exo_coldens: failed to get number of days in year id'
call handle_ncerr( nf_inq_dimid( ncid, 'ndays_of_year', dimid ), trim(err_msg) )
err_msg = 'get_exo_coldens: failed to get number of days in year'
call handle_ncerr( nf_inq_dimlen( ncid, dimid, ndays_of_year ), trim(err_msg) )
err_msg = 'get_exo_coldens: failed to get west_east id'
call handle_ncerr( nf_inq_dimid( ncid, 'west_east', dimid ), trim(err_msg) )
err_msg = 'get_exo_coldens: failed to get west_east'
call handle_ncerr( nf_inq_dimlen( ncid, dimid, lon_e ), trim(err_msg) )
err_msg = 'get_exo_coldens: failed to get south_north id'
call handle_ncerr( nf_inq_dimid( ncid, 'south_north', dimid ), trim(err_msg) )
err_msg = 'get_exo_coldens: failed to get south_north'
call handle_ncerr( nf_inq_dimlen( ncid, dimid, lat_e ), trim(err_msg) )
! end IF
#ifdef DM_PARALLEL
!---------------------------------------------------------------------
! ... bcast the dimensions
!---------------------------------------------------------------------
! CALL wrf_dm_bcast_bytes ( ncoldens_levs , IWORDSIZE )
! CALL wrf_dm_bcast_bytes ( ndays_of_year , IWORDSIZE )
! CALL wrf_dm_bcast_bytes ( lon_e , IWORDSIZE )
! CALL wrf_dm_bcast_bytes ( lat_e , IWORDSIZE )
#endif
!---------------------------------------------------------------------
! ... allocate local arrays
!---------------------------------------------------------------------
iend = min( ite,ide-1 )
jend = min( jte,jde-1 )
! allocate( coldens(lon_e,lat_e,ncoldens_levs,ndays_of_year), stat=astat )
! if( astat /= 0 ) then
! call wrf_error_fatal( 'get_exo_coldens: failed to allocate coldens' )
! end if
!---------------------------------------------------------------------
! ... allocate column_density type component arrays
!---------------------------------------------------------------------
col_dens(dm)%ncoldens_levs = ncoldens_levs
col_dens(dm)%ndays_of_year = ndays_of_year
allocate( col_dens(dm)%col_levs(ncoldens_levs), &
col_dens(dm)%day_of_year(ndays_of_year), stat=astat )
if( astat /= 0 ) then
call wrf_error_fatal( 'get_exo_coldens: failed to allocate col_levs,day_of_year' )
end if
allocate( col_dens(dm)%o3_col_dens(its:iend,jts:jend,ncoldens_levs,ndays_of_year), &
col_dens(dm)%o2_col_dens(its:iend,jts:jend,ncoldens_levs,ndays_of_year), stat=astat )
if( astat /= 0 ) then
call wrf_error_fatal( 'get_exo_coldens: failed to allocate o3_col_dens,o2_col_dens' )
end if
col_dens(dm)%is_allocated = .true.
!---------------------------------------------------------------------
! ... read arrays
!---------------------------------------------------------------------
! IF ( wrf_dm_on_monitor() ) THEN
err_msg = 'get_exo_coldens: failed to get col_levs variable id'
call handle_ncerr( nf_inq_varid( ncid, 'coldens_levs', varid ), trim(err_msg) )
err_msg = 'get_exo_coldens: failed to read col_levs variable'
call handle_ncerr( nf_get_var_double( ncid, varid, col_dens(dm)%col_levs ), trim(err_msg) )
err_msg = 'get_exo_coldens: failed to get days_of_year variable id'
call handle_ncerr( nf_inq_varid( ncid, 'days_of_year', varid ), trim(err_msg) )
err_msg = 'get_exo_coldens: failed to read days_of_year variable'
call handle_ncerr( nf_get_var_double( ncid, varid, col_dens(dm)%day_of_year ), trim(err_msg) )
err_msg = 'get_exo_coldens: failed to get o3 col_dens variable id'
call handle_ncerr( nf_inq_varid( ncid, 'o3_column_density', varid ), trim(err_msg) )
err_msg = 'get_exo_coldens: failed to read o3 col_dens variable'
! call handle_ncerr( nf_get_var_real( ncid, varid, coldens ), trim(err_msg) )
call handle_ncerr( nf_get_vara_double( ncid, varid, (/its,jts,1,1/), &
(/iend-its+1,jend-jts+1,ncoldens_levs,ndays_of_year/), &
col_dens(dm)%o3_col_dens(its:iend,jts:jend,1:ncoldens_levs,1:ndays_of_year) ), trim(err_msg) )
! ENDIF
#ifdef DM_PARALLEL
! call wrf_debug( 0,'get_exo_coldens: bcast col_levs' )
! DM_BCAST_MACRO(col_dens(dm)%col_levs)
! call wrf_debug( 0,'get_exo_coldens: bcast day_of_year' )
! DM_BCAST_MACRO(col_dens(dm)%day_of_year)
! call wrf_debug( 0,'get_exo_coldens: bcast o3_col_dens' )
! CALL wrf_dm_bcast_bytes ( coldens, size(coldens)*RWORDSIZE )
#endif
! col_dens(dm)%o3_col_dens(its:iend,jts:jend,:ncoldens_levs,:ndays_of_year) = &
! coldens(its:iend,jts:jend,:ncoldens_levs,:ndays_of_year)
! IF ( wrf_dm_on_monitor() ) THEN
err_msg = 'get_exo_coldens: failed to get o2 col_dens variable id'
call handle_ncerr( nf_inq_varid( ncid, 'o2_column_density', varid ), trim(err_msg) )
err_msg = 'get_exo_coldens: failed to read o2 col_dens variable'
! call handle_ncerr( nf_get_var_real( ncid, varid, coldens ), trim(err_msg) )
call handle_ncerr( nf_get_vara_double( ncid, varid, (/its,jts,1,1/), &
(/iend-its+1,jend-jts+1,ncoldens_levs,ndays_of_year/), &
col_dens(dm)%o2_col_dens(its:iend,jts:jend,1:ncoldens_levs,1:ndays_of_year) ), trim(err_msg) )
!---------------------------------------------------------------------
! ... close column density netcdf file
!---------------------------------------------------------------------
err_msg = 'get_exo_coldens: failed to close file ' // trim(filename)
call handle_ncerr( nf_close( ncid ), trim(err_msg) )
! end if
#ifdef DM_PARALLEL
! call wrf_debug( 0,'get_exo_coldens: bcast o2_col_dens' )
! CALL wrf_dm_bcast_bytes ( coldens, size(coldens)*RWORDSIZE )
#endif
! col_dens(dm)%o2_col_dens(its:iend,jts:jend,:ncoldens_levs,:ndays_of_year) = &
! coldens(its:iend,jts:jend,:ncoldens_levs,:ndays_of_year)
! deallocate( coldens )
!---------------------------------------------------------------------
! ... some diagnostics
!---------------------------------------------------------------------
call wrf_debug( 100,' ' )
write(err_msg,'(''get_exo_coldens: coldens variables for domain '',i2)') dm
call wrf_debug( 100,trim(err_msg) )
call wrf_debug( 100,'get_exo_coldens: days_of_year' )
do k = 1,ndays_of_year,5
write(err_msg,'(1p,5g15.7)') col_dens(dm)%day_of_year(k:min(k+4,ndays_of_year))
call wrf_debug( 100,trim(err_msg) )
end do
call wrf_debug( 100,'get_exo_coldens: coldens levels' )
do k = 1,ncoldens_levs,5
write(err_msg,'(1p,5g15.7)') col_dens(dm)%col_levs(k:min(k+4,ncoldens_levs))
call wrf_debug( 100,trim(err_msg) )
end do
call wrf_debug( 100,' ' )
endif col_dens_allocated
#endif
end subroutine get_exo_coldens
subroutine tuv_timestep_init( dm, julday )
!-----------------------------------------------------------------------------
! ... setup the exo column time interpolation
!-----------------------------------------------------------------------------
implicit none
!-----------------------------------------------------------------------------
! ... dummy arguments
!-----------------------------------------------------------------------------
integer, intent(in) :: dm ! domain number
integer, intent(in) :: julday ! day of year at present time step
!-----------------------------------------------------------------------------
! ... local variables
!-----------------------------------------------------------------------------
integer :: m
real(dp) :: calday
calday = real( julday,kind=dp)
if( calday < col_dens(dm)%day_of_year(1) ) then
next = 1
last = 12
dels = (365._dp + calday - col_dens(dm)%day_of_year(12)) &
/ (365._dp + col_dens(dm)%day_of_year(1) - col_dens(dm)%day_of_year(12))
else if( calday >= col_dens(dm)%day_of_year(12) ) then
next = 1
last = 12
dels = (calday - col_dens(dm)%day_of_year(12)) &
/ (365. + col_dens(dm)%day_of_year(1) - col_dens(dm)%day_of_year(12))
else
do m = 11,1,-1
if( calday >= col_dens(dm)%day_of_year(m) ) then
exit
end if
end do
last = m
next = m + 1
dels = (calday - col_dens(dm)%day_of_year(m)) / (col_dens(dm)%day_of_year(m+1) - col_dens(dm)%day_of_year(m))
end if
end subroutine tuv_timestep_init
subroutine z_interp( z_tab, tab, ntab, z_out, out )
integer, intent(in) :: ntab
real, intent(in) :: z_tab(ntab)
real, intent(in) :: tab(ntab)
real, intent(in) :: z_out(:)
real, intent(out) :: out(:)
!---------------------------------------------------------------
! ... altitude interpolation
!---------------------------------------------------------------
integer :: k, kt, ktm1, n
real :: delz
n = size(out)
do k = 1,n
if( z_out(k) <= z_tab(1) ) then
out(k) = z_tab(1)
elseif( z_out(k) >= z_tab(ntab) ) then
out(k) = z_tab(ntab)
else
do kt = 1,ntab
if( z_tab(kt) >= z_out(k) ) then
ktm1 = kt - 1
delz = (z_out(k) - z_tab(ktm1))/(z_tab(kt) - z_tab(ktm1))
out(k) = tab(ktm1) + delz*(tab(kt) - tab(ktm1))
exit
endif
end do
endif
end do
end subroutine z_interp
subroutine p_interp( o2_exo_col, o3_exo_col, o3_exo_col_at_grnd, ptop, &
dm, its, ite, jts, jte )
!---------------------------------------------------------------
! ... pressure interpolation for exo col density
!---------------------------------------------------------------
implicit none
!---------------------------------------------------------------
! ... dummy arguments
!---------------------------------------------------------------
integer, intent(in) :: dm
integer, intent(in) :: its, ite
integer, intent(in) :: jts, jte
real(dp), intent(in) :: ptop(its:ite,jts:jte) ! pressure at photolysis top (hPa)
real(dp), intent(out) :: o2_exo_col(its:ite,jts:jte) ! exo model o2 column density (molecules/cm^2)
real(dp), intent(out) :: o3_exo_col(its:ite,jts:jte) ! exo model o3 column density (molecules/cm^2)
real(dp), intent(out) :: o3_exo_col_at_grnd(its:ite,jts:jte) ! exo model o3 column density at grnd (molecules/cm^2)
!---------------------------------------------------------------
! ... local variables
!---------------------------------------------------------------
integer :: i, j, k, ku, kl
integer :: Kgrnd
real(dp) :: pinterp
real(dp) :: delp
real(dp) :: tint_vals(2)
Kgrnd = col_dens(dm)%ncoldens_levs
lat_loop : &
do j = jts,jte
long_loop : &
do i = its,ite
pinterp = ptop(i,j)
if( pinterp < col_dens(dm)%col_levs(1) ) then
ku = 1
kl = 1
delp = 0._dp
else
do ku = 2,col_dens(dm)%ncoldens_levs
if( pinterp <= col_dens(dm)%col_levs(ku) ) then
kl = ku - 1
delp = log( pinterp/col_dens(dm)%col_levs(kl) )/log( col_dens(dm)%col_levs(ku)/col_dens(dm)%col_levs(kl) )
exit
end if
end do
end if
tint_vals(1) = col_dens(dm)%o2_col_dens(i,j,kl,last) &
+ delp * (col_dens(dm)%o2_col_dens(i,j,ku,last) &
- col_dens(dm)%o2_col_dens(i,j,kl,last))
tint_vals(2) = col_dens(dm)%o2_col_dens(i,j,kl,next) &
+ delp * (col_dens(dm)%o2_col_dens(i,j,ku,next) &
- col_dens(dm)%o2_col_dens(i,j,kl,next))
o2_exo_col(i,j) = tint_vals(1) + dels * (tint_vals(2) - tint_vals(1))
tint_vals(1) = col_dens(dm)%o3_col_dens(i,j,kl,last) &
+ delp * (col_dens(dm)%o3_col_dens(i,j,ku,last) &
- col_dens(dm)%o3_col_dens(i,j,kl,last))
tint_vals(2) = col_dens(dm)%o3_col_dens(i,j,kl,next) &
+ delp * (col_dens(dm)%o3_col_dens(i,j,ku,next) &
- col_dens(dm)%o3_col_dens(i,j,kl,next))
o3_exo_col(i,j) = tint_vals(1) + dels * (tint_vals(2) - tint_vals(1))
tint_vals(1) = col_dens(dm)%o3_col_dens(i,j,Kgrnd,last)
tint_vals(2) = col_dens(dm)%o3_col_dens(i,j,Kgrnd,next)
o3_exo_col_at_grnd(i,j) = tint_vals(1) + dels * (tint_vals(2) - tint_vals(1))
end do long_loop
end do lat_loop
end subroutine p_interp
subroutine xsqy_int( n, xsqy, tlev, dens_air )
!---------------------------------------------------------------------
! ... interpolate m,t tables for xs * qy
!---------------------------------------------------------------------
integer, intent(in) :: n
real, intent(in) :: tlev(:)
real, intent(in) :: dens_air(:)
real, intent(out) :: xsqy(:,:)
real, parameter :: m0 = 2.45e19
integer :: tndx, mndx, tndxp1, mndxp1
integer :: k, ku
real :: temp, dens
real :: w(4)
real :: del_t, del_d
ku = size( tlev )
do k = 1,ku
temp = tlev(k)
do tndx = 1,ntemp
if( temp_tab(tndx) > temp ) then
exit
endif
end do
tndx = max( min( tndx,ntemp ) - 1,1 )
tndxp1 = tndx + 1
del_t = max( 0.,min( 1.,(temp - temp_tab(tndx))*del_temp_tab(tndx) ) )
! dens = dens_air(k)
dens = dens_air(k)/m0
do mndx = 1,nconc
if( conc_tab(mndx) > dens ) then
exit
endif
end do
mndx = max( min( mndx,nconc ) - 1,1 )
mndxp1 = mndx + 1
del_d = max( 0.,min( 1.,(dens - conc_tab(mndx))*del_conc_tab(mndx) ) )
w(1) = (1. - del_t)*(1. - del_d)
w(2) = del_t*(1. - del_d)
w(3) = (1. - del_t)*del_d
w(4) = del_t*del_d
xsqy(1:nwave,k) = w(1)*xsqy_tab(1:nwave,tndx,mndx,n) &
+ w(2)*xsqy_tab(1:nwave,tndxp1,mndx,n) &
+ w(3)*xsqy_tab(1:nwave,tndx,mndxp1,n) &
+ w(4)*xsqy_tab(1:nwave,tndxp1,mndxp1,n)
end do
end subroutine xsqy_int
subroutine xs_int( xs, tlev, xs_tab )
!---------------------------------------------------------------------
! ... interpolate tables for xs
!---------------------------------------------------------------------
real, intent(in) :: tlev(:)
real, intent(in) :: xs_tab(:,:)
real, intent(out) :: xs(:,:)
integer :: tndx, tndxp1
integer :: k, ku
real :: temp
real :: w(2)
real :: del_t
ku = size( tlev )
do k = 1,ku
temp = tlev(k)
do tndx = 1,ntemp
if( temp_tab(tndx) > temp ) then
exit
endif
end do
tndx = max( min( tndx,ntemp ) - 1,1 )
tndxp1 = tndx + 1
del_t = max( 0.,min( 1.,(temp - temp_tab(tndx))*del_temp_tab(tndx) ) )
w(1) = (1. - del_t)
w(2) = del_t
xs(1:nwave,k) = w(1)*xs_tab(1:nwave,tndx) &
+ w(2)*xs_tab(1:nwave,tndxp1)
end do
end subroutine xs_int
FUNCTION chap(zeta)
! chapman function is used when the solar zenith angle exceeds
! 75 deg.
! interpolates between values given in, e.g., mccartney (1976).
! .. Scalar Arguments ..
REAL, intent(in) :: zeta
! .. Local Scalars ..
REAL :: rm
INTEGER :: i
LOGICAL :: fnd
! .. Local Arrays ..
REAL :: y(75:96)
! .. Function Return Value ..
REAL :: chap
! .. Data Statements ..
DATA (y(i),i=75,96) &
/3.800, 4.055, 4.348, 4.687, 5.083, 5.551, 6.113, &
6.799, 7.650, 8.732, 10.144, 12.051, 14.730, 18.686, 24.905, 35.466, &
55.211, 96.753, 197., 485., 1476., 9999./
fnd = .false.
DO i = 75, 96
rm = real(i)
IF (zeta < rm) then
chap = y(i) + (y(i+1) - y(i))*(zeta - (rm - 1.))
fnd = .true.
exit
ENDIF
END DO
IF( .not. fnd ) then
chap = y(96)
ENDIF
END FUNCTION chap
SUBROUTINE fery( nwave, w, wght )
!-----------------------------------------------------------------------------
!= PURPOSE:
!= Calculate the action spectrum value for erythema at a given wavelength
!= according to: McKinlay, A.F and B.L.Diffey, A reference action spectrum
!= for ultraviolet induced erythema in human skin, CIE Journal, vol 6,
!= pp 17-22, 1987.
!= Value at 300 nm = 0.6486
!-----------------------------------------------------------------------------
INTEGER, intent(in) :: nwave
REAL, intent(in) :: w(:)
REAL, intent(out) :: wght(:)
WHERE( w(1:nwave) < 298. )
wght(1:nwave) = 1.
ELSEWHERE( w(1:nwave) >= 298. .AND. w(1:nwave) < 328. )
wght(1:nwave) = 10.**(0.094*(298. - w(1:nwave)))
ELSEWHERE( w(1:nwave) >= 328. .AND. w(1:nwave) < 400. )
wght(1:nwave) = 10.**(0.015*(139. - w(1:nwave)))
ELSEWHERE( w(1:nwave) >= 400. )
wght(1:nwave) = 1.e-36
ENDWHERE
END SUBROUTINE fery
SUBROUTINE cldfrac_binary( CLDFRA,QC,QI, QS, kts, kte )
!---------------------------------------------------------------------
! !DESCRIPTION:
! Compute cloud fraction from input ice, snow, and cloud water fields
! if provided.
!
! Whether QI or QC is active or not is determined from the indices of
! the fields into the 4D scalar arrays in WRF. These indices are
! P_QI and P_QC, respectively, and they are passed in to the routine
! to enable testing to see if QI and QC represent active fields in
! the moisture 4D scalar array carried by WRF.
!
! If a field is active its index will have a value greater than or
! equal to PARAM_FIRST_SCALAR, which is also an input argument to
! this routine.
!EOP
!---------------------------------------------------------------------
IMPLICIT NONE
INTEGER, INTENT(IN) :: kts, kte
REAL, INTENT(OUT ) :: CLDFRA(kts:kte)
REAL, INTENT(IN) :: QI(kts:kte), &
QC(kts:kte), &
QS(kts:kte)
!----------------------------------------------------------------------
! local variables
!----------------------------------------------------------------------
REAL, parameter :: thresh = 1.e-9
INTEGER :: j
where( (qc(kts:kte) + qi(kts:kte) + qs(kts:kte)) > thresh )
cldfra(kts:kte) = 1.
elsewhere
cldfra(kts:kte) = 0.
endwhere
END SUBROUTINE cldfrac_binary
SUBROUTINE cldfrac_fractional( CLDFRA, QV, QC, QI, QS, &
p_phy, t_phy, &
kts, kte )
!----------------------------------------------------------------------
! dummy arguments
!----------------------------------------------------------------------
INTEGER, INTENT(IN) :: kts,kte
REAL, INTENT(OUT) :: CLDFRA(kts:kte)
REAL, INTENT(IN) :: QV(kts:kte), &
QI(kts:kte), &
QC(kts:kte), &
QS(kts:kte), &
t_phy(kts:kte), &
p_phy(kts:kte)
!----------------------------------------------------------------------
! local variables
!----------------------------------------------------------------------
REAL , PARAMETER :: ALPHA0 = 100.
REAL , PARAMETER :: GAMMA = 0.49
REAL , PARAMETER :: QCLDMIN = 1.E-12
REAL , PARAMETER :: PEXP = 0.25
REAL , PARAMETER :: RHGRID =1.0
REAL , PARAMETER :: SVP1 = 0.61078
REAL , PARAMETER :: SVP2 = 17.2693882
REAL , PARAMETER :: SVPI2 = 21.8745584
REAL , PARAMETER :: SVP3 = 35.86
REAL , PARAMETER :: SVPI3 = 7.66
REAL , PARAMETER :: SVPT0 = 273.15
REAL , PARAMETER :: r_d = 287.
REAL , PARAMETER :: r_v = 461.6
REAL , PARAMETER :: ep_2 = r_d/r_v
INTEGER :: i,j,k
INTEGER :: imax, jmax, kmax
REAL :: RHUM, tc, esw, esi, weight, qvsw, qvsi, qvs_weight
REAL :: QCLD, DENOM, ARG, SUBSAT, wrk
REAL :: relhum_max, wrk_max
! !DESCRIPTION:
!----------------------------------------------------------------------
! Compute cloud fraction from input ice and cloud water fields
! if provided.
!
! Whether QI or QC is active or not is determined from the indices of
! the fields into the 4D scalar arrays in WRF. These indices are
! P_QI and P_QC, respectively, and they are passed in to the routine
! to enable testing to see if QI and QC represent active fields in
! the moisture 4D scalar array carried by WRF.
!
! If a field is active its index will have a value greater than or
! equal to PARAM_FIRST_SCALAR, which is also an input argument to
! this routine.
!----------------------------------------------------------------------
!EOP
!-----------------------------------------------------------------------
! COMPUTE GRID-SCALE CLOUD COVER FOR RADIATION
! (modified by Ferrier, Feb '02)
!
! Cloud fraction parameterization follows Randall, 1994
! (see Hong et al., 1998)
!-----------------------------------------------------------------------
! Note: ep_2=287./461.6 Rd/Rv
! Note: R_D=287.
! Alternative calculation for critical RH for grid saturation
! RHGRID=0.90+.08*((100.-DX)/95.)**.5
! Calculate saturation mixing ratio weighted according to the fractions of
! water and ice.
! Following:
! Murray, F.W. 1966. ``On the computation of Saturation Vapor Pressure'' J. Appl. Meteor. 6 p.204
! es (in mb) = 6.1078 . exp[ a . (T-273.16)/ (T-b) ]
!
! over ice over water
! a = 21.8745584 17.2693882
! b = 7.66 35.86
!---------------------------------------------------------------------
relhum_max = -100.
wrk_max = -10000.
imax = 0; kmax = 0; jmax = 0
vert_loop: &
DO k = kts,kte
!---------------------------------------------------------------------
! Determine cloud fraction (modified from original algorithm)
!---------------------------------------------------------------------
QCLD = QI(k) + QC(k) + QS(k)
has_cloud : &
IF( QCLD >= QCLDMIN ) THEN
tc = t_phy(k) - SVPT0
esw = 1000.0 * SVP1 * EXP( SVP2 * tc / ( t_phy(k) - SVP3 ) )
esi = 1000.0 * SVP1 * EXP( SVPI2 * tc / ( t_phy(k) - SVPI3 ) )
QVSW = EP_2 * esw / ( p_phy(k) - esw )
QVSI = EP_2 * esi / ( p_phy(k) - esi )
weight = (QI(k) + QS(k)) / QCLD
QVS_WEIGHT = (1. - weight)*QVSW + weight*QVSI
RHUM = QV(k)/QVS_WEIGHT !--- Relative humidity
!---------------------------------------------------------------------
! Assume zero cloud fraction if there is no cloud mixing ratio
!---------------------------------------------------------------------
IF( RHUM >= RHGRID )THEN
!---------------------------------------------------------------------
! Assume cloud fraction of unity if near saturation and the cloud
! mixing ratio is at or above the minimum threshold
!---------------------------------------------------------------------
CLDFRA(k) = 1.
ELSE
!---------------------------------------------------------------------
! Adaptation of original algorithm (Randall, 1994; Zhao, 1995)
! modified based on assumed grid-scale saturation at RH=RHgrid.
!---------------------------------------------------------------------
SUBSAT = MAX( 1.E-10,RHGRID*QVS_WEIGHT - QV(k) )
DENOM = SUBSAT**GAMMA
ARG = MAX( -6.9,-ALPHA0*QCLD/DENOM ) ! <-- EXP(-6.9)=.001
!---------------------------------------------------------------------
! prevent negative values (new)
!---------------------------------------------------------------------
RHUM = MAX( 1.E-10, RHUM )
wrk = (RHUM/RHGRID)**PEXP*(1. - EXP( ARG ))
IF( wrk >= .01 ) then
CLDFRA(k) = wrk
ENDIF
ENDIF
ENDIF has_cloud
END DO vert_loop
END SUBROUTINE cldfrac_fractional
#ifdef NETCDF
subroutine handle_ncerr( ret, mes )
!---------------------------------------------------------------------
! ... netcdf error handling routine
!---------------------------------------------------------------------
implicit none
!---------------------------------------------------------------------
! ... dummy arguments
!---------------------------------------------------------------------
integer, intent(in) :: ret
character(len=*), intent(in) :: mes
include 'netcdf.inc'
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
#endif
end module module_phot_tuv