module esfsw_driver_mod
!
!
! fil
!
!
! smf
!
!
!
!
! Code that initializes and calculates shortwave radiative quantities
! such as flux and heating rate.
!
!
! This code initializes the necessary shortwave radiative parameters
! in the initialization subroutine. It then uses delta-eddington approximation
! and doubling and adding technique to calculate solar flux and
! heating rate.
!
!
! shared modules:
use mpp_mod, only: input_nml_file
use fms_mod, only: open_namelist_file, fms_init, &
mpp_pe, mpp_root_pe, stdlog, &
file_exist, write_version_number, &
check_nml_error, error_mesg, &
FATAL, close_file
use constants_mod, only: PI, GRAV, radcon_mks, o2mixrat, &
rhoair, pstd_mks, WTMAIR, &
constants_init
! shared radiation package modules:
use esfsw_parameters_mod, only: Solar_spect, esfsw_parameters_init, &
esfsw_parameters_end
use rad_utilities_mod, only: Rad_control, rad_utilities_init, &
cldrad_properties_type, &
cld_specification_type, &
astronomy_type, &
aerosol_diagnostics_type, &
radiative_gases_type, &
aerosol_type, aerosol_properties_type,&
Cldrad_control, &
atmos_input_type, surface_type, &
sw_output_type, Sw_control
use tracer_manager_mod, only : get_tracer_index, &
NO_TRACER
use field_manager_mod, only : MODEL_ATMOS
!---------------------------------------------------------------------
implicit none
private
!--------------------------------------------------------------------
! esfsw_driver_mod is the internal driver for the esf shortwave
! package.
!------------------------------------------------------------------
!---------------------------------------------------------------------
!----------- version number for this module -------------------
character(len=128) :: version = '$Id: esfsw_driver.F90,v 19.0 2012/01/06 20:15:51 fms Exp $'
character(len=128) :: tagname = '$Name: tikal $'
!---------------------------------------------------------------------
!------- interfaces --------
public &
esfsw_driver_init, swresf, &
esfsw_driver_end
private &
! called from swresf:
adding, deledd
!---------------------------------------------------------------------
!-------- namelist ---------
logical :: do_ica_calcs=.false. ! do independent column
! calculations when sto-
! chastic clouds are
! active ?
logical :: do_rayleigh_all_bands = .true. ! rayleigh scattering
! calculated in all sw
! bands ?
logical :: do_herzberg = .false. ! include the herzberg
! effect on the o2
! optical depth ?
logical :: do_quench = .false. ! include the quenching
! effect of non-LTE
! processes on the co2
! optical depth ?
logical :: do_ch4_sw_effects = .false. ! the shortwave effects
! of ch4 are included ?
logical :: do_n2o_sw_effects = .false. ! the shortwave effects
! of n2o are included ?
logical :: do_coupled_stratozone = .false. ! include the coupled
! stratospheric ozone effects?
namelist / esfsw_driver_nml / &
do_ica_calcs, do_rayleigh_all_bands, &
do_herzberg, do_quench, &
do_ch4_sw_effects, do_n2o_sw_effects, &
do_coupled_stratozone
!---------------------------------------------------------------------
!------- public data ------
!---------------------------------------------------------------------
!------- private data ------
!---------------------------------------------------------------------
! variables associated with absorptivity and sw transmission for
! various gaseous atmospheric components
!
! powph2o = the scaling factor used in the fit of the h2o
! transmission function
!
! p0h2o = the reference pressure (mb) used in the fit of the
! h2o transmission function
!
! c(n)co2(str) = coefficients for the absorptivity expression for co2
! for the pressure-scaled and non-scaled, respectively,
! portions of the fit (=1.0e-99 if no absorption)
!
! c(n)o2(str) = coefficients for the absorptivity expression for o2
! for the pressure-scaled and non-scaled, respectively,
! portions of the fit (=1.0e-99 if no absorption)
!
! ""(schrun) = coefficients for the absorptivity expression for the
! Schuman-Runge o2 band (non-scaled only)
!
! kh2o = the psuedo-absorption coefficients in cm2/gm for h2o
!
! ko3 = the absorption coefficients in cm2/gm for o3
!
! wtfreq = the weight associated with each exponential term
! strterm = logical flag to indicate whether or not a h2o pseudo-
! absorption coefficient is assigned a non-scaled
! (true) or pressure-scaled (false) gas amount
!---------------------------------------------------------------------
real, dimension (:), allocatable :: c1co2, c1co2str, c1o2, c1o2str, &
c2co2, c2co2str, c2o2, c2o2str, &
c3co2, c3co2str, c3o2, c3o2str, &
c4co2, c4co2str, c4o2, c4o2str, &
c1ch4, c1ch4str, c2ch4, &
c2ch4str, c3ch4, c3ch4str, &
c4ch4, c4ch4str, &
c1n2o, c1n2ostr, c2n2o, &
c2n2ostr, c3n2o, c3n2ostr, &
c4n2o, c4n2ostr, &
powph2o, p0h2o
real :: c1o2strschrun, c2o2strschrun, &
c3o2strschrun, c4o2strschrun
real, dimension (:), allocatable :: kh2o, ko3, wtfreq
logical, dimension(:), allocatable :: strterm
!---------------------------------------------------------------------
! quantities associated with solar spectral parameterization
!
! firstrayband = the first band number where the contribution by
! rayleigh scattering is included in the solar
! calculations
!
! nirbands = the number of bands in the near-infrared (used in
! assigning the value of the surface albedo for the
! near-infrared, and the visible and ultraviolet
! regions, separately)
! nfreqpts = the number of pseudo-monochromatic frequencies
! solflxband = the solar flux in each parameterization band
! solflxbandref = the solar flux in each parameterization band, used for
! defining band-averaged optical parameters. If the
! solar constant is time-invariant, it is also the solar
! flux in each parameterization band (solflxband).
! vis_wvnum = the wavenumber of visible light (corresponds to
! wavelength of 0.55 microns) [ cm **(-1) ]
!---------------------------------------------------------------------
real :: refquanray, solflxtotal
integer :: firstrayband, nirbands
integer, dimension (:), allocatable :: nfreqpts
real, dimension(:), allocatable :: solflxband
real, dimension(:), allocatable :: solflxbandref
real, dimension(:), allocatable :: wtstr, cosangstr
real, dimension(4) :: wtstr_4 = &
(/0.347854845, 0.652145155,&
0.347854845, 0.652145155/)
integer :: nbands, tot_wvnums, nfrqpts, nh2obands, nstreams
logical :: nstr4 = .false.
real :: vis_wvnum = 1.0E+04/0.55
real :: wvnum_870 = 1.0E+04/0.87
real :: wvnum_340 = 1.0E+04/0.34
real :: wvnum_380 = 1.0E+04/0.38
real :: wvnum_440 = 1.0E+04/0.44
real :: wvnum_670 = 1.0E+04/0.67
real :: one_micron_wvnum = 1.0E+04/1.00
real :: onepsix_micron_wvnum = 1.0E+04/1.61
integer :: onepsix_band_indx
!---------------------------------------------------------------------
! variables associated with rayleigh scattering
!---------------------------------------------------------------------
real, dimension (:), allocatable :: betaddensitymol
!----------------------------------------------------------------------
! variables associated with total optical path of species ? - smf
!----------------------------------------------------------------------
real :: toto2strmaxschrun
real, dimension(:), allocatable :: totco2max, totco2strmax, &
toto2max, toto2strmax, &
totch4max, totch4strmax, &
totn2omax, totn2ostrmax
!----------------------------------------------------------------------
! variables associated with the herzberg effect. wtmo2 is the mol-
! ecular weight of o2. herzberg_fac is a factor used in the last
! shortwave band, so that herzberg_fac*wo2 yields a correction for
! the o2 optical depth to account for the herzberg o2 heating. this
! is done only when do_herzberg is true.
!----------------------------------------------------------------------
real, parameter :: wtmo2 = 3.19988E+01
real, parameter :: herzberg_fac = 9.9488377E-3
!----------------------------------------------------------------------
! variables associated with the quenching effect. co2_quenchfac is a
! multiplication factor that reduces the co2 gas optical depth, and
! hence, the solar heating in the upper atmosphere, to account for
! "quenching" due to non-LTE processes. co2_quenchfac_height is the
! reference height for co2_quenchfac [ meters ].
!----------------------------------------------------------------------
real, dimension(30) :: co2_quenchfac
data co2_quenchfac /1.0,.954,.909,.853,.800,.747,.693,.637,.583, .526,&
.467,.416,.368,.325,.285,.253,.229,.206,.186,.170,&
.163,.156,.151,.144,.138,.132,.127,.124,.068,.037/
real, dimension(30) :: co2_quenchfac_height
data co2_quenchfac_height /67304.,68310.,69303.,70288.,71267.,72245.,&
73221.,74195.,75169.,76141.,77112.,78082.,&
79051.,80018.,80985.,81950.,82914.,83876.,&
84838.,85798.,86757.,87715.,88672.,89627.,&
90582.,91535.,92487.,93438.,94387.,106747./
!---------------------------------------------------------------------
! miscellaneous variables
!---------------------------------------------------------------------
integer, parameter :: NSOLWG = 1
real, dimension(NSOLWG) :: gausswt
logical :: module_is_initialized = .false.
logical :: do_esfsw_band_diagnostics = .false.
integer :: naerosol_optical, naerosoltypes_used
!---------------------------------------------------------------------
!---------------------------------------------------------------------
contains
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!
! PUBLIC SUBROUTINES
!
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!
!
! Subroutine that defines the time-independent quantities associated
! with the incoming shortwave radiation in the multiple-band solar
! radiation parameterization.
!
!
! It first reads in the input namelist and then allocates gas absorption
! coefficient variables. It then reads in the shortwave input namelist
! file and assigns the gas absorption coefficients. Rayleigh scattering
! coefficient is also calculated based on the temperature and pressure
! structure of the atmosphere.
!
!
! call esfsw_driver_init
!
!
!
subroutine esfsw_driver_init
!----------------------------------------------------------------------
! esfsw_driver_init is the constructor for esfsw_driver_mod. it
! defines the time-independent quantities associated with the
! incoming shortwave radiation in the multiple-band solar radiation
! parameterization.
!----------------------------------------------------------------------
!---------------------------------------------------------------------
! local variables:
real, dimension(Solar_spect%nbands) :: freqnu
real, dimension(Solar_spect%nstreams) :: ptstr
integer, dimension(0:Solar_spect%nbands) :: endwvnbands
integer, dimension(:), allocatable :: nwvnsolar
real , dimension(:), allocatable :: solint
character(len=64) :: file_name
real, dimension(4) :: ptstr_4 = (/-0.861136312,&
-0.339981044, &
0.861136312, &
0.339981044 /)
real :: ptstr_1 = 0.2
real :: temprefray = 288.15
real :: pressrefray = 101325. ! MKS units
real :: densmolref = 2.54743E+19
real :: convfac = 1.0E+18
real :: corrfac, gamma, f1, f2, f3, pinteg, &
twopiesq, densmolrefsqt3, wavelength, &
freqsq, ristdm1, ri
integer :: iounit, nband, nf, ni, nw, nw1, nw2, nintsolar
integer :: unit, io, ierr, logunit
integer :: i
integer :: n
real :: input_flag = 1.0e-99
!---------------------------------------------------------------------
! local variables:
!
! freqnu
! ptstr gaussian points and weights for evaluation of the
! diffuse beam.
! nwvnsolar the number of wavenumbers in each region where the
! solar flux is constant
! solint the solar flux in watts per meter**2 in each
! wavenumber region where it is constant
! endwvnbands the wavenumber value for the band limits
! file_name
! ptstr_4
! ptstr_1
! temprefray reference temperature used in defining rayleigh
! optical depth
! pressrefray reference pressure used in defining rayleigh
! optical depth [ Pa ]
! densmolref reference density used in defining rayleigh
! optical depth
! convfac
! corrfac
! gamma
! f1
! f2
! f3
! pinteg
! twopiesq
! densmolrefsqt3
! wavelength
! freqsq
! ristdm1
! ri
! iounit
! nband
! nf
! ni
! nw
! nw1
! nw2
! nintsolar the number of wavenumber regions where the
! solar flux is constant.
! unit
! io
! ierr
! i
!
!---------------------------------------------------------------------
!---------------------------------------------------------------------
! if routine has already been executed, exit.
!---------------------------------------------------------------------
if (module_is_initialized) return
!---------------------------------------------------------------------
! verify that modules used by this module that are not called later
! have already been initialized.
!---------------------------------------------------------------------
call fms_init
call constants_init
call rad_utilities_init
call esfsw_parameters_init
!-----------------------------------------------------------------------
! read namelist.
!-----------------------------------------------------------------------
#ifdef INTERNAL_FILE_NML
read (input_nml_file, nml=esfsw_driver_nml, iostat=io)
ierr = check_nml_error(io,'esfsw_driver_nml')
#else
if ( file_exist('input.nml')) then
unit = open_namelist_file ( )
ierr=1; do while (ierr /= 0)
read (unit, nml=esfsw_driver_nml, iostat=io, end=10)
ierr = check_nml_error(io,'esfsw_driver_nml')
end do
10 call close_file (unit)
endif
#endif
!---------------------------------------------------------------------
! write version number and namelist to logfile.
!---------------------------------------------------------------------
call write_version_number(version, tagname)
logunit = stdlog()
if (mpp_pe() == mpp_root_pe() ) &
write (logunit, nml=esfsw_driver_nml)
!---------------------------------------------------------------------
! define flag indicating if ICA calculations being done.
!---------------------------------------------------------------------
Cldrad_control%do_ica_calcs = do_ica_calcs
Cldrad_control%do_ica_calcs_iz = .true.
!---------------------------------------------------------------------
! allocate module variables
!---------------------------------------------------------------------
nbands = Solar_spect%nbands
tot_wvnums = Solar_spect%tot_wvnums
nfrqpts = Solar_spect%nfrqpts
nstreams = Solar_spect%nstreams
nh2obands = Solar_spect%nh2obands
allocate ( betaddensitymol (nbands) )
allocate ( c1co2 (nh2obands), &
c1co2str(nh2obands), &
c1o2 (nh2obands), &
c1o2str (nh2obands), &
c2co2 (nh2obands), &
c2co2str(nh2obands), &
c2o2 (nh2obands), &
c2o2str (nh2obands), &
c3co2 (nh2obands), &
c3co2str(nh2obands), &
c3o2 (nh2obands), &
c3o2str (nh2obands), &
c4co2 (nh2obands), &
c4co2str(nh2obands), &
c4o2 (nh2obands), &
c4o2str (nh2obands), &
c1ch4 (nh2obands), &
c1ch4str(nh2obands), &
c2ch4 (nh2obands), &
c2ch4str(nh2obands), &
c3ch4 (nh2obands), &
c3ch4str(nh2obands), &
c4ch4 (nh2obands), &
c4ch4str(nh2obands), &
c1n2o (nh2obands), &
c1n2ostr(nh2obands), &
c2n2o (nh2obands), &
c2n2ostr(nh2obands), &
c3n2o (nh2obands), &
c3n2ostr(nh2obands), &
c4n2o (nh2obands), &
c4n2ostr(nh2obands), &
powph2o (nh2obands), &
p0h2o (nh2obands) )
allocate ( nfreqpts (nbands) )
allocate ( solflxband (nbands) )
allocate ( solflxbandref (nbands) )
allocate ( kh2o (nfrqpts), &
ko3 (nfrqpts), &
wtfreq (nfrqpts), &
strterm (nfrqpts) )
allocate ( wtstr (nstreams), &
cosangstr (nstreams) )
allocate ( totco2max (nh2obands), &
totco2strmax (nh2obands), &
toto2max (nh2obands), &
toto2strmax (nh2obands), &
totch4max (nh2obands), &
totch4strmax (nh2obands), &
totn2omax (nh2obands), &
totn2ostrmax (nh2obands) )
betaddensitymol = 0.0 ; c1co2 = 0.0 ; c1co2str = 0.0
c1o2 = 0.0 ; c1o2str = 0.0 ; c2co2 = 0.0
c2co2str = 0.0 ; c2o2 = 0.0 ; c2o2str = 0.0
c3co2 = 0.0 ; c3co2str = 0.0 ; c3o2 = 0.0
c3o2str = 0.0 ; c4co2 = 0.0 ; c4co2str = 0.0
c4o2 = 0.0 ; c4o2str = 0.0 ; powph2o = 0.0
p0h2o = 0.0 ; nfreqpts = 0.0 ; solflxband = 0.0
solflxbandref = 0.0 ; kh2o = 0.0 ; ko3 = 0.0
wtfreq = 0.0 ; strterm = .FALSE. ; wtstr = 0.0
cosangstr = 0.0 ; totco2max = 0.0 ; totco2strmax = 0.0
toto2max = 0.0 ; toto2strmax = 0.0
c1ch4 = 0.0 ; c1ch4str = 0.0; c2ch4 = 0.0
c2ch4str = 0.0 ; c3ch4 = 0.0 ; c3ch4str = 0.0
c4ch4 = 0.0 ; c4ch4str = 0.0
totch4max = 0.0 ; totch4strmax = 0.0
c1n2o = 0.0 ; c1n2ostr = 0.0; c2n2o = 0.0
c2n2ostr = 0.0 ; c3n2o = 0.0 ; c3n2ostr = 0.0
c4n2o = 0.0 ; c4n2ostr = 0.0
totn2omax = 0.0 ; totn2ostrmax = 0.0
!---------------------------------------------------------------------
! allocate local variables.
!---------------------------------------------------------------------
if (nstreams == 4) then
ptstr(:) = ptstr_4(:)
wtstr(:) = wtstr_4(:)
nstr4 = .true.
else if (nstreams == 1) then
ptstr(1) = ptstr_1
wtstr(1) = 1.0
nstr4 = .false.
endif
!---------------------------------------------------------------------
! read input file for band positions, solar fluxes and band
! strengths.
!---------------------------------------------------------------------
if (nbands == 25 .and. nfrqpts == 72) then
file_name = 'INPUT/esf_sw_input_data_n72b25'
else if (nbands == 18 .and. nfrqpts == 38) then
file_name = 'INPUT/esf_sw_input_data_n38b18'
else
call error_mesg ( 'esfsw_driver_mod', &
'input file for desired bands and frqs is not available', &
FATAL)
endif
iounit = open_namelist_file (file_name)
read(iounit,101) ( solflxbandref(nband), nband=1,NBANDS )
read(iounit,102) ( nfreqpts(nband), nband=1,NBANDS )
read(iounit,103) ( endwvnbands(nband), nband=1,NBANDS )
read(iounit,103) FIRSTRAYBAND,NIRBANDS
read(iounit,104) ( powph2o(nband), nband=1,NH2OBANDS )
read(iounit,104) ( p0h2o(nband), nband=1,NH2OBANDS )
read(iounit,105) ( c1co2(nband), nband=1,NH2OBANDS )
read(iounit,105) ( c1co2str(nband), nband=1,NH2OBANDS )
read(iounit,105) ( c2co2(nband), nband=1,NH2OBANDS )
read(iounit,105) ( c2co2str(nband), nband=1,NH2OBANDS )
read(iounit,105) ( c3co2(nband), nband=1,NH2OBANDS )
read(iounit,105) ( c3co2str(nband), nband=1,NH2OBANDS )
read(iounit,105) ( c1o2(nband), nband=1,NH2OBANDS )
read(iounit,105) ( c1o2str(nband), nband=1,NH2OBANDS ), &
c1o2strschrun
read(iounit,105) ( c2o2(nband), nband=1,NH2OBANDS )
read(iounit,105) ( c2o2str(nband), nband=1,NH2OBANDS ), &
c2o2strschrun
read(iounit,105) ( c3o2(nband), nband=1,NH2OBANDS )
read(iounit,105) ( c3o2str(nband), nband=1,NH2OBANDS ), &
c3o2strschrun
read(iounit,105) ( c1ch4(nband), nband=1,NH2OBANDS )
read(iounit,105) ( c1ch4str(nband), nband=1,NH2OBANDS )
read(iounit,105) ( c2ch4(nband), nband=1,NH2OBANDS )
read(iounit,105) ( c2ch4str(nband), nband=1,NH2OBANDS )
read(iounit,105) ( c3ch4(nband), nband=1,NH2OBANDS )
read(iounit,105) ( c3ch4str(nband), nband=1,NH2OBANDS )
read(iounit,105) ( c1n2o(nband), nband=1,NH2OBANDS )
read(iounit,105) ( c1n2ostr(nband), nband=1,NH2OBANDS )
read(iounit,105) ( c2n2o(nband), nband=1,NH2OBANDS )
read(iounit,105) ( c2n2ostr(nband), nband=1,NH2OBANDS )
read(iounit,105) ( c3n2o(nband), nband=1,NH2OBANDS )
read(iounit,105) ( c3n2ostr(nband), nband=1,NH2OBANDS )
do nf = 1,nfrqpts
read(iounit,106) wtfreq(nf),kh2o(nf),ko3(nf),strterm (nf)
end do
read(iounit,107) nintsolar
allocate ( nwvnsolar (nintsolar) )
allocate ( solint (nintsolar) )
do ni = 1,nintsolar
read(iounit,107) nwvnsolar (ni),solint(ni)
end do
call close_file (iounit)
if (Solar_spect%tot_wvnums /= &
endwvnbands(Solar_spect%nbands)) then
call error_mesg ( 'esfsw_driver_mod', &
' inconsistency between highest solar spectrum wavenumber '//&
'in esfsw_parameters_mod and in esfsw_sriver input file', &
FATAL)
endif
!---------------------------------------------------------------------
! define the band index corresponding to near ultraviolet light
! (0.34 microns). note that endwvnbands is in units of (cm**-1).
!---------------------------------------------------------------------
do ni=1,nbands
if (endwvnbands(ni) > wvnum_340) then
Solar_spect%w340_band_indx = ni
Solar_spect%w340_band_iz = .true.
exit
endif
end do
!---------------------------------------------------------------------
! define the band index corresponding to near ultraviolet light
! (0.38 microns). note that endwvnbands is in units of (cm**-1).
!---------------------------------------------------------------------
do ni=1,nbands
if (endwvnbands(ni) > wvnum_380) then
Solar_spect%w380_band_indx = ni
Solar_spect%w380_band_iz = .true.
exit
endif
end do
!---------------------------------------------------------------------
! define the band index corresponding to blue light
! (0.44 microns). note that endwvnbands is in units of (cm**-1).
!---------------------------------------------------------------------
do ni=1,nbands
if (endwvnbands(ni) > wvnum_440) then
Solar_spect%w440_band_indx = ni
Solar_spect%w440_band_iz = .true.
exit
endif
end do
!---------------------------------------------------------------------
! define the band index corresponding to red light
! (0.67 microns). note that endwvnbands is in units of (cm**-1).
!---------------------------------------------------------------------
do ni=1,nbands
if (endwvnbands(ni) > wvnum_670) then
Solar_spect%w670_band_indx = ni
Solar_spect%w670_band_iz = .true.
exit
endif
end do
!---------------------------------------------------------------------
! define the band index corresponding to visible light
! (0.55 microns). note that endwvnbands is in units of (cm**-1).
!---------------------------------------------------------------------
do ni=1,nbands
if (endwvnbands(ni) > vis_wvnum) then
Solar_spect%visible_band_indx = ni
Solar_spect%visible_band_indx_iz = .true.
exit
endif
end do
!---------------------------------------------------------------------
! define the band index corresponding to 870nm
! (0.87 microns). note that endwvnbands is in units of (cm**-1).
!---------------------------------------------------------------------
do ni=1,nbands
if (endwvnbands(ni) > wvnum_870) then
Solar_spect%eight70_band_indx = ni
Solar_spect%eight70_band_indx_iz = .true.
exit
endif
end do
!---------------------------------------------------------------------
! define the band index corresponding to near infra red band
! (1.00 microns). note that endwvnbands is in units of (cm**-1).
!---------------------------------------------------------------------
do ni=1,nbands
if (endwvnbands(ni) > one_micron_wvnum) then
Solar_spect%one_micron_indx = ni
Solar_spect%one_micron_indx_iz = .true.
exit
endif
end do
!---------------------------------------------------------------------
! define the band index corresponding to
! (1.61 microns). note that endwvnbands is in units of (cm**-1).
!---------------------------------------------------------------------
do ni=1,nbands
if (endwvnbands(ni) > onepsix_micron_wvnum) then
onepsix_band_indx = ni
exit
endif
end do
!--------------------------------------------------------------------
! define the wavenumber one solar fluxes.
!----------------------------------------------------------------- --
do ni = 1,nintsolar
if ( ni.eq.1 ) then
nw1 = 1
else
nw1 = nw1 + nwvnsolar(ni-1)
end if
nw2 = nw1 + nwvnsolar(ni) - 1
do nw = nw1,nw2
Solar_spect%solarfluxtoa(nw) = solint(ni)
end do
end do
!---------------------------------------------------------------------
!
!---------------------------------------------------------------------
deallocate (solint )
deallocate (nwvnsolar )
!---------------------------------------------------------------------
!
!---------------------------------------------------------------------
if ( .not. Rad_control%using_solar_timeseries_data) then
Solar_spect%solflxband = solflxbandref
endif
Solar_spect%solflxbandref = solflxbandref
Solar_spect%endwvnbands = endwvnbands
!--------------------------------------------------------------------
! override the input file value of firstrayband if the nml control
! variables indicates rayleigh effects are to be considered in
! all bands.
!--------------------------------------------------------------------
if (do_rayleigh_all_bands) firstrayband = 1
!----------------------------------------------------------------------
! convert some values to mks to match model units
!--------------------------------------------------------------------
p0h2o = 1.0E-2/p0h2o ! invert, and convert mb to mks
kh2o = kh2o *1.0E-01 ! cgs to mks
ko3 = ko3 *1.0E-01 ! cgs to mks
!----------------------------------------------------------------------
!
!----------------------------------------------------------------------
do n=1,NH2OBANDS
if (c1co2(n) /= input_flag .and. &
c2co2(n) /= input_flag .and. &
c3co2(n) /= input_flag ) then
c4co2(n) = c1co2(n) * c2co2(n) ** c3co2(n)
c4co2str(n) = c1co2str(n) * c2co2str(n) ** c3co2str(n)
totco2max(n) = ( (1.0/c1co2(n) ) + c2co2(n) ** c3co2(n) ) ** &
(1.0/c3co2(n) ) - c2co2(n)
if (nbands == 18) then
if ( n /= 4) then
totco2strmax(n) = ( (1.0/c1co2str(n) ) + c2co2str(n) ** &
c3co2str(n) ) ** (1.0/c3co2str(n) ) - &
c2co2str(n)
else
totco2strmax(n) = HUGE (c4o2strschrun)
endif
else
totco2strmax(n) = ( (1.0/c1co2str(n) ) + c2co2str(n) ** &
c3co2str(n) ) ** (1.0/c3co2str(n) ) - &
c2co2str(n)
endif
else
c4co2(n) = 0.0
c4co2str(n) = 0.0
totco2max(n) = 0.0
totco2strmax(n) = 0.0
endif
if (c1o2(n) /= input_flag .and. &
c2o2(n) /= input_flag .and. &
c3o2(n) /= input_flag ) then
c4o2(n) = c1o2(n) * c2o2(n) ** c3o2(n)
c4o2str(n) = c1o2str(n) * c2o2str(n) ** c3o2str(n)
toto2max(n) = ( (1.0/c1o2(n) ) + c2o2(n) ** c3o2(n) ) ** &
(1.0/c3o2(n) ) - c2o2(n)
if (nbands == 18) then
if ( n /= 4) then
toto2strmax(n) = ( (1.0/c1o2str(n) ) + c2o2str(n) ** &
c3o2str(n) ) ** (1.0/c3o2str(n) ) - &
c2o2str(n)
else
toto2strmax(n) = HUGE (c4o2strschrun)
endif
else
toto2strmax(n) = ( (1.0/c1o2str(n) ) + c2o2str(n) ** &
c3o2str(n) ) ** (1.0/c3o2str(n) ) - &
c2o2str(n)
endif
else
c4o2(n) = 0.0
c4o2str(n) = 0.0
toto2max(n) = 0.0
toto2strmax(n) = 0.0
endif
if (do_ch4_sw_effects) then
if (c1ch4(n) /= input_flag .and. &
c2ch4(n) /= input_flag .and. &
c3ch4(n) /= input_flag ) then
c4ch4(n) = c1ch4(n) * c2ch4(n) ** c3ch4(n)
c4ch4str(n) = c1ch4str(n) * c2ch4str(n) ** c3ch4str(n)
totch4max(n) = ( (1.0/c1ch4(n) ) + c2ch4(n) ** c3ch4(n) ) ** &
(1.0/c3ch4(n) ) - c2ch4(n)
totch4strmax(n) = ( (1.0/c1ch4str(n) ) + c2ch4str(n) ** &
c3ch4str(n) ) ** (1.0/c3ch4str(n) ) - &
c2ch4str(n)
else
c4ch4(n) = 0.
c4ch4str(n) = 0.
totch4max(n) = 0.
totch4strmax(n) = 0.
endif
endif
if (do_n2o_sw_effects) then
if (c1n2o(n) /= input_flag .and. &
c2n2o(n) /= input_flag .and. &
c3n2o(n) /= input_flag ) then
c4n2o(n) = c1n2o(n) * c2n2o(n) ** c3n2o(n)
c4n2ostr(n) = c1n2ostr(n) * c2n2ostr(n) ** c3n2ostr(n)
totn2omax(n) = ( (1.0/c1n2o(n) ) + c2n2o(n) ** c3n2o(n) ) ** &
(1.0/c3n2o(n) ) - c2n2o(n)
totn2ostrmax(n) = ( (1.0/c1n2ostr(n) ) + c2n2ostr(n) ** &
c3n2ostr(n) ) ** (1.0/c3n2ostr(n) ) - &
c2n2ostr(n)
else
c4n2o(n) = 0.
c4n2ostr(n) = 0.
totn2omax(n) = 0.
totn2ostrmax(n) = 0.
endif
endif
end do
c4o2strschrun = c1o2strschrun * c2o2strschrun ** c3o2strschrun
toto2strmaxschrun = ( (1.0/c1o2strschrun) + c2o2strschrun ** &
c3o2strschrun) ** (1.0/c3o2strschrun) - &
c2o2strschrun
! if (mpp_pe() == 0) then
! print *, 'c1ch4 ', c1ch4
! print *, 'c2ch4 ', c2ch4
! print *, 'c3ch4 ', c3ch4
! print *, 'c4ch4 ', c4ch4
! print *, 'totch4max', totch4max
! print *, 'c1ch4str ', c1ch4str
! print *, 'c2ch4str ', c2ch4str
! print *, 'c3ch4str ', c3ch4str
! print *, 'c4ch4str ', c4ch4str
! print *, 'totch4strmax', totch4strmax
! print *, 'c1co2 ', c1co2
! print *, 'c2ch4 ', c2co2
! print *, 'c3ch4 ', c3co2
! print *, 'c4ch4 ', c4co2
! print *, 'totch4max', totco2max
! print *, 'c1ch4str ', c1co2str
! print *, 'c2ch4str ', c2co2str
! print *, 'c3ch4str ', c3co2str
! print *, 'c4ch4str ', c4co2str
! print *, 'totch4strmax', totco2strmax
! print *, 'c1n2o ', c1n2o
! print *, 'c2ch4 ', c2n2o
! print *, 'c3ch4 ', c3n2o
! print *, 'c4ch4 ', c4n2o
! print *, 'totch4max', totn2omax
! print *, 'c1ch4str ', c1n2ostr
! print *, 'c2ch4str ', c2n2ostr
! print *, 'c3ch4str ', c3n2ostr
! print *, 'c4ch4str ', c4n2ostr
! print *, 'totch4strmax', totn2ostrmax
! print *, '1/c1', 1.0/c1co2str
! print *, '1/c3', 1.0/c3co2str
! print *, 'c2**c3', c2co2str**c3co2str
! print *, '(1 / + c2**c3)**1/c3', &
! (1.0/c1co2str + c2co2str**c3co2str)**(1./c3co2str)
! print *, 'c4o2 ', c4o2
! print *, 'c4o2str ', c4o2str
! print *, 'toto2max', toto2max
! print *, 'toto2strmax', toto2strmax
! endif
!----------------------------------------------------------------------
!
!----------------------------------------------------------------------
if ( .not. Rad_control%using_solar_timeseries_data) then
solflxtotal = 0.0
do nband = 1,NBANDS
solflxtotal = solflxtotal + Solar_spect%solflxbandref(nband)
end do
endif
!----------------------------------------------------------------------
! define the wavenumbers to evaluate rayleigh optical depth.
!-------------------------------------------------------------------
endwvnbands(0) = 0
do nband = 1,NBANDS
freqnu(nband) = 0.5 * ( endwvnbands(nband-1) + &
endwvnbands(nband) )
end do
!---------------------------------------------------------------------
! define quantities used to determine the rayleigh optical depth.
! notes: refquanray is the quantity which multiplies pressure /
! temperature to yield the molecular density.
! betaddensitymol is the quantity which multiples the
! molecular density to yield the rayleigh scattering
! coefficient.
! 1.39E-02 is the depolorization factor.
!-----------------------------------------------------------------
refquanray = densmolref * temprefray / pressrefray
corrfac = ( 6.0E+00 + 3.0E+00 * 1.39E-02 )/( 6.0E+00 - 7.0E+00 * &
1.39E-02 )
gamma = 1.39E-02 / ( 2.0E+00 - 1.39E-02 )
f1 = 7.5E-01 / ( gamma * 2.0E+00 + 1.0E+00 )
f2 = gamma * 3.0E+00 + 1.0E+00
f3 = 1.0E+00 - gamma
pinteg = 2.0E+00 * PI * ( 2.0E+00 * f1 * f2 * ( 1.0E+00 + f3 / &
f2 / 3.0E+00 ) )
twopiesq = 2.0E+00 * PI ** 2
densmolrefsqt3 = 3.0E+00 * densmolref ** 2
!---------------------------------------------------------------------
!
!---------------------------------------------------------------------
do nband = 1,NBANDS
wavelength = 1.0E+04 / freqnu(nband)
freqsq = 1.0 / ( wavelength ) ** 2
ristdm1 = ( 6.4328E+03 + 2.94981E+06 / ( 1.46E+02 - freqsq ) + &
2.554E+04 / ( 4.1E+01 - freqsq ) ) * 1.0E-08
ri = ristdm1 + 1
betaddensitymol(nband) = twopiesq*( ri ** 2 - 1.0E+00 ) ** 2 * &
corrfac / ( densmolrefsqt3 * &
wavelength ** 4 ) * pinteg * convfac
end do
gausswt(1) = 1.0
!---------------------------------------------------------------------
! define the gaussian angles for evaluation of the diffuse beam.
!--------------------------------------------------------------------
do i = 1,nstreams
cosangstr(i) = ( ptstr(i) + 1. ) * 5.0E-01
end do
!--------------------------------------------------------------------
! mark the module as initialized.
!--------------------------------------------------------------------
module_is_initialized = .true.
!--------------------------------------------------------------------
101 format( 12f10.4 )
102 format( 32i4 )
103 format( 20i6 )
104 format( 12f10.2 )
105 format( 1p,16e8.1 )
106 format( 1p,3e16.6,l16 )
107 format( i5,1p,e14.5 )
!---------------------------------------------------------------------
end subroutine esfsw_driver_init
!#################################################################
!
!
! Subroutine that uses the delta-eddington technique in conjunction
! with a multi-band parameterization for h2o+co2+o2+o3 absorption
! in the solar spectrum to derive solar fluxes and heating rates.
!
!
! This subroutine calculates optical depth, single scattering albedo,
! asymmetry parameter of a layer based on gaseous absorbers,
! clouds, aerosols, and rayleigh scattering. It then uses delta-
! eddington technique to calculate radiative flux at each layer.
! Doubling and adding technique is used to combine the layers
! and calculate flux at TOA and surface and heating rate. This
! subroutine allocates a substantial amount of memory and deallocates
! the allocated memory at the end of the subroutine.
!
!
! call swresf(is, ie, js, je, Atmos_input, Surface, Rad_gases, Aerosol,
! Astro, &
! Cldrad_props, Cld_spec, Sw_output)
!
!
! starting subdomain i indices of data in
! the physics_window being integrated
!
!
! ending subdomain i indices of data in
! the physics_window being integrated
!
!
! starting subdomain j indices of data in
! the physics_window being integrated
!
!
! ending subdomain j indices of data in
! the physics_window being integrated
!
!
! Atmos_input_type variable containing the atmospheric
! input fields on the radiation grid
!
!
! Aerosol input data for shortwave radiation calculation
!
!
! Astronomy_type variable containing the astronomical
! input fields on the radiation grid
!
!
! Radiative_gases_type variable containing the radiative
! gas input fields on the radiation grid
!
!
! The cloud radiative property input fields on the
! radiation grid
!
!
! The shortwave radiation calculation result
!
!
! Surface data as boundary condition to radiation
!
!
! Cloud specification data as initial condition to radiation
!
!
subroutine swresf (is, ie, js, je, Atmos_input, Surface, Rad_gases, &
Aerosol, Aerosol_props, Astro, Cldrad_props, &
Cld_spec, including_volcanoes, Sw_output, &
Aerosol_diags, r, including_aerosols, &
naerosol_optical)
!----------------------------------------------------------------------
! swresf uses the delta-eddington technique in conjunction with a
! multiple-band parameterization for h2o+co2+o2+o3 absorption to
! derive solar fluxes and heating rates.
! notes: drops are assumed if temp>273.15K, ice crystals otherwise.
!-------------------------------------------------------------------
integer, intent(in) :: is, ie, js, je
type(atmos_input_type), intent(in) :: Atmos_input
type(surface_type), intent(in) :: Surface
type(radiative_gases_type), intent(in) :: Rad_gases
type(aerosol_type), intent(in) :: Aerosol
type(aerosol_properties_type), intent(in) :: Aerosol_props
type(astronomy_type), intent(in) :: Astro
type(cldrad_properties_type), intent(in) :: Cldrad_props
type(cld_specification_type), intent(in) :: Cld_spec
logical, intent(in) :: including_volcanoes
type(sw_output_type), intent(inout) :: Sw_output
type(aerosol_diagnostics_type),intent(inout) :: Aerosol_diags
real,dimension(:,:,:,:), intent(inout) :: r
logical, intent(in) :: including_aerosols
integer, intent(in) :: naerosol_optical
!-------------------------------------------------------------------
! intent(in) variables:
!
! is,ie,js,je starting/ending subdomain i,j indices of data in
! the physics_window being integrated
! Atmos_input atmos_input_type structure, contains variables
! defining atmospheric state
! Surface surface_type structure, contains variables
! defining the surface characteristics
! Rad_gases radiative_gases_type structure, contains var-
! iables defining the radiatively active gases,
! passed through to lower level routines
! Aerosol aerosol_type structure, contains variables
! defining aerosol fields, passed through to
! lower level routines
! Aerosol_props aerosol radiative property input data for the
! radiation package
! Astro astronomy_type structure
! Cldrad_props cloud radiative property input data for the
! radiation package
! Cld_spec cld_specification_type structure, contains var-
! iables defining the cloud distribution, passed
! through to lower level routines
!
! intent(inout) variables:
!
! Sw_output shortwave radiation output data
!
!---------------------------------------------------------------------
!-----------------------------------------------------------------------
! local variables:
logical, dimension (size(Atmos_input%temp,1), &
size(Atmos_input%temp,2), &
size(Atmos_input%temp,3)-1) :: &
cloud
logical, dimension (size(Atmos_input%temp,1), &
size(Atmos_input%temp,2)) :: &
cloud_in_column, daylight
real, dimension (size(Atmos_input%temp,1), &
size(Atmos_input%temp,2), &
size(Atmos_input%temp,3)-1, &
Solar_spect%nbands) :: &
aeroasymfac, aerosctopdep, aeroextopdep, &
rayopdep
real, dimension (size(Atmos_input%temp,1), &
size(Atmos_input%temp,2), &
size(Atmos_input%temp,3)-1, &
nfrqpts,NSOLWG) :: &
gasopdep
real, dimension (size(Atmos_input%temp,1), &
size(Atmos_input%temp,2), &
size(Atmos_input%temp,3)-1) :: &
cloudfrac, cldfrac_band, cldfrac, &
deltaz, cloud_deltaz, &
cldext, cldsct, cldasymm, &
cloudasymfac, cloudextopdep, &
cloudsctopdep, deltap, &
densitymol, extopdepclr, &
extopdepovc, fclr, fovc, &
gocpdp, gclr, gstrclr, &
gstrovc, govc, omegastrclr, &
omegastrovc, rlayerdif, rlayerdir, &
rlayerdifclr, rlayerdifovc, &
rlayerdirclr, rlayerdirovc, &
sctopdepclr, sctopdepovc, &
ssalbclr, ssalbovc, &
taustrclr, taustrovc, &
tlayerde, tlayerdeclr, &
tlayerdeovc, tlayerdif, &
tlayerdifclr, tlayerdifovc, &
tlayerdir, tlayerdirclr, &
tlayerdirovc
real, dimension (size(Atmos_input%temp,1), &
size(Atmos_input%temp,2), &
size(Atmos_input%temp,3)) :: &
reflectance, transmittance, tr_dir
real, dimension (size(Atmos_input%temp,1), &
size(Atmos_input%temp,2), &
size(Atmos_input%temp,3), &
Rad_control%nzens) :: &
dfswbandclr, fswbandclr, ufswbandclr, &
dfswband, fswband, ufswband, &
sumtrclr, sumreclr, sumtr_dir_clr, &
sumtr, sumre, sumtr_dir
real, dimension (size(Atmos_input%temp,1), &
size(Atmos_input%temp,2), &
Rad_control%nzens) :: sumtr_dir_up
real, dimension (size(Atmos_input%temp,1), &
size(Atmos_input%temp,2), &
size(Atmos_input%temp,3)) :: &
press, pflux, pflux_mks, &
temp, &
reflectanceclr, transmittanceclr, tr_dirclr
real, dimension (size(Atmos_input%temp,1), &
size(Atmos_input%temp,2), &
size(Atmos_input%temp,3)-1 , &
Rad_control%nzens) :: &
hswbandclr, hswband
real, dimension (size(Atmos_input%temp,1), &
size(Atmos_input%temp,2)) :: &
sfcalb_dir, sfcalb_dif, wtfac_p, &
fracday_p, solarflux_p
real, dimension (size(Atmos_input%press,1), &
size(Atmos_input%press,2), &
NSOLWG) :: &
cosangsolar_p
integer :: j, i, k, ng, np, nband, nf, ns, nz
integer :: nzens
integer :: nprofile, nprofiles
real :: profiles_inverse
integer :: ix, jx, kx, israd, jsrad, ierad, jerad, ksrad, kerad
real :: ssolar
real :: solflxtotal_local
!-----------------------------------------------------------------------
! local variables:
!
! aeramt
! sum_g_omega_tau
! opt_index_v3
! irh
! etc.
!
!--------------------------------------------------------------------
!---------------------------------------------------------------------
! be sure module has been initialized.
!---------------------------------------------------------------------
if (.not. module_is_initialized ) then
call error_mesg ('esfsw_driver_mod', &
'module has not been initialized', FATAL )
endif
!---------------------------------------------------------------------
! define the solar_constant appropriate at Rad_time when solar
! input is varying with time.
!---------------------------------------------------------------------
!-- The following fix is for openmp
if (Rad_control%using_solar_timeseries_data) then
solflxtotal_local = 0.0
do nband = 1,NBANDS
solflxtotal_local = solflxtotal_local + Solar_spect%solflxband(nband)
end do
else
solflxtotal_local = solflxtotal
endif
!---------------------------------------------------------------------
!
cldfrac = Cld_spec%camtsw
!---------------------------------------------------------------------
! convert to cgs and then back to mks for consistency with previous
!---------------------------------------------------------------------
press(:,:,:) = 0.1*(10.0*Atmos_input%press(:,:,:))
pflux(:,:,:) = (10.0*Atmos_input%pflux(:,:,:))
deltaz(:,:,:) = Atmos_input%deltaz(:,:,:)
temp (:,:,:) = Atmos_input%temp (:,:,:)
cloud_deltaz = Atmos_input%clouddeltaz
!--------------------------------------------------------------------
! define limits and dimensions
!--------------------------------------------------------------------
ix = size(temp,1)
jx = size(temp,2)
kx = size(temp,3) - 1
israd = 1
jsrad = 1
ksrad = 1
ierad = ix
jerad = jx
kerad = kx
!----------------------------------------------------------------------c
! define a flag indicating columns in which there is sunshine during
! this radiation time step. define a flag indicating points with both
! sunlight and cloud.
!----------------------------------------------------------------------c
do j = JSRAD,JERAD
do i = ISRAD,IERAD
if ( Astro%fracday(i,j) /= 0.0) then
daylight(i,j) = .true.
else
daylight(i,j) = .false.
endif
cloud_in_column(i,j) = .false.
end do
end do
call compute_aerosol_optical_props &
(Atmos_input, Aerosol, Aerosol_props, &
including_volcanoes, Aerosol_diags, r, &
including_aerosols, naerosol_optical, &
daylight, aeroextopdep, aerosctopdep, aeroasymfac)
!---------------------------------------------------------------------
!
!---------------------------------------------------------------------
ssolar = Sw_control%solar_constant*Astro%rrsun
!----------------------------------------------------------------------
! define a flag indicating points with both sunlight and cloud. set
! the flag indicating that cloud is present in columns with such
! points.
!----------------------------------------------------------------------
if (.not. Cldrad_control%do_stochastic_clouds) then
do j = JSRAD,JERAD
do i = ISRAD,IERAD
if (daylight(i,j)) then
do k=KSRAD,KERAD
if (cldfrac(i,j,k) > 0.0) then
cloud_in_column(i,j) = .true.
cloud(i,j,k) = .true.
cloudfrac(i,j,k) = cldfrac(i,j,k)
else
cloud(i,j,k) = .false.
cloudfrac(i,j,k) = 0.0
endif
end do
else
do k=KSRAD,KERAD
cloud(i,j,k) = .false.
cloudfrac(i,j,k) = 0.0
end do
endif
end do
end do
endif
!----------------------------------------------------------------------c
! define pressure related quantities, pressure is in mks units.
!----------------------------------------------------------------------c
pflux_mks = pflux*1.0E-1
do k = KSRAD+1,KERAD+1
deltap(:,:,k-1) = pflux_mks(:,:,k) - pflux_mks(:,:,k-1)
gocpdp(:,:,k-1) = radcon_mks/deltap(:,:,k-1)
end do
call compute_gas_props (Atmos_input, Rad_gases, Astro, &
daylight, gasopdep)
!---------------------------------------------------------------------
! define the molecular density for use in calculating the
! rayleigh optical depth (deltaz is in meters).
!--------------------------------------------------------------------
do k = KSRAD,KERAD
densitymol(:,:,k) = refquanray * press(:,:,k) / temp(:,:,k)
end do
! assumption is that there is 1 cloud profile for each sw band
if (do_ica_calcs) then
nprofiles = nbands
profiles_inverse = 1.0/nprofiles
else
nprofiles = 1
profiles_inverse = 1.0
endif
!--------------------------------------------------------------------
! define the rayleigh optical depths.
!---------------------------------------------------------------------
do nband = 1, NBANDS
rayopdep(:,:,:,nband) = betaddensitymol(nband)* &
densitymol(:,:,:)*deltaz(:,:,:)
end do ! (nband loop)
nzens = Rad_control%nzens
!--------------------------------------------------------------------
do nprofile=1, nprofiles
if (do_ica_calcs) then
cldfrac_band(:,:,:) = Cld_spec%camtsw_band(:,:,:,nprofile)
endif
!----------------------------------------------------------------------c
! np is a counter for the pseudo-monochromatic frequency point
! number.
!----------------------------------------------------------------------c
np = 0
do nz=1,nzens
if (Rad_control%do_totcld_forcing) then
dfswbandclr(:,:,:,nz) = 0.0
fswbandclr(:,:,:,nz) = 0.0
hswbandclr(:,:,:,nz) = 0.0
ufswbandclr(:,:,:,nz) = 0.0
endif
end do
reflectanceclr = 0.0
transmittanceclr = 0.0
!----------------------------------------------------------------------c
! begin band loop
!----------------------------------------------------------------------c
do nband = 1,NBANDS
do nz=1,nzens
sumtr(:,:,:,nz) = 0.0
sumtr_dir(:,:,:,nz) = 0.0
sumtr_dir_up(:,:,nz) = 0.0
sumre(:,:,:,nz) = 0.0
if (Rad_control%do_totcld_forcing) then
sumtrclr(:,:,:,nz) = 0.0
sumreclr(:,:,:,nz) = 0.0
sumtr_dir_clr(:,:,:,nz) = 0.0
endif
end do
if (Cldrad_control%do_stochastic_clouds) then
if (.not. do_ica_calcs) then
cldfrac_band(:,:,:) = Cld_spec%camtsw_band(:,:,:,nband)
endif
endif
!----------------------------------------------------------------------
! if stochastic clouds are activated (cloud fields differ with sw
! parameterization band), define a flag indicating points with both
! sunlight and cloud for the current sw parameterization band. set
! the flag indicating that cloud is present in columns with such
! points.
!----------------------------------------------------------------------
if (Cldrad_control%do_stochastic_clouds) then
do j = JSRAD,JERAD
do i = ISRAD,IERAD
cloud_in_column(i,j) = .false.
if (daylight(i,j)) then
do k = KSRAD,KERAD
if (cldfrac_band(i,j,k) > 0.0) then
cloud_in_column(i,j) = .true.
cloud(i,j,k) = .true.
cloudfrac(i,j,k) = cldfrac_band(i,j,k)
else
cloud(i,j,k) = .false.
cloudfrac(i,j,k) = 0.0
endif
end do
else
do k = KSRAD,KERAD
cloud(i,j,k) = .false.
cloudfrac(i,j,k) = 0.0
end do
endif
end do
end do
endif
!---------------------------------------------------------------------
! obtain cloud properties from the Cldrad_props input variable.
!--------------------------------------------------------------------
cldext(:,:,:) = Cldrad_props%cldext(:,:,:,nband,nprofile)
cldsct(:,:,:) = Cldrad_props%cldsct(:,:,:,nband,nprofile)
cldasymm(:,:,:) = Cldrad_props%cldasymm(:,:,:,nband,nprofile)
do k = KSRAD,KERAD
do j=JSRAD,JERAD
do i=ISRAD, IERAD
if (cloud(i,j,k) ) then
cloudextopdep(i,j,k) = 1.0E-03*cldext(i,j,k) * &
cloud_deltaz(i,j,k)
cloudsctopdep(i,j,k) = 1.0E-03*cldsct(i,j,k) * &
cloud_deltaz(i,j,k)
cloudasymfac(i,j,k) = cldasymm(i,j,k)
endif
end do
end do
end do
!-----------------------------------------------------------------
! define clear sky arrays
!-----------------------------------------------------------------
if (nband >= firstrayband ) then
do k=ksrad,kerad
do j=jsrad,jerad
do i=israd,ierad
if (daylight(i,j) ) then
sctopdepclr(i,j,k) = rayopdep(i,j,k,nband) + &
aerosctopdep(i,j,k,nband)
gclr(i,j,k) = aeroasymfac(i,j,k,nband)* &
aerosctopdep(i,j,k,nband)/&
sctopdepclr(i,j,k)
fclr(i,j,k) = aeroasymfac(i,j,k,nband)*gclr(i,j,k)
gstrclr(i,j,k) = ( gclr(i,j,k) - fclr(i,j,k) )/ &
( 1.0 - fclr(i,j,k) )
endif
end do
end do
end do
endif
!-----------------------------------------------------------------
! define cloudy sky arrays
!-----------------------------------------------------------------
do k=KSRAD,KERAD
do j=JSRAD,JERAD
do i=ISRAD,IERAD
if (cloud(i,j,k)) then
sctopdepovc(i,j,k) = rayopdep(i,j,k,nband) + &
aerosctopdep(i,j,k,nband) + &
cloudsctopdep(i,j,k)
govc(i,j,k) = ( ( cloudasymfac(i,j,k) * &
cloudsctopdep(i,j,k) ) + &
( aeroasymfac(i,j,k,nband) * &
aerosctopdep(i,j,k,nband)))/ &
sctopdepovc(i,j,k)
fovc(i,j,k) = ( ( cloudasymfac(i,j,k) ** 2 * &
cloudsctopdep(i,j,k) ) + &
( aeroasymfac(i,j,k,nband) ** 2 * &
aerosctopdep(i,j,k,nband) ))/ &
sctopdepovc(i,j,k)
gstrovc(i,j,k) = ( govc(i,j,k) - fovc(i,j,k))/ &
( 1.0 - fovc(i,j,k) )
endif
end do
end do
end do
!---------------------------------------------------------------------
! begin frequency points in the band loop.
!--------------------------------------------------------------------
do nf = 1,nfreqpts(nband)
np = np + 1
!---------------------------------------------------------------------
! begin gaussian angle loop (ng > 1 only when lswg = true).
!--------------------------------------------------------------------
do ng = 1,NSOLWG
!---------------------------------------------------------------------
! clear sky mode
! note: in this mode, the delta-eddington method is performed for all
! spatial columns experiencing sunlight.
!--------------------------------------------------------------------
if (nband >= firstrayband ) then
do k=ksrad,kerad
do j=jsrad,jerad
do i=israd,ierad
if (daylight(i,j) ) then
extopdepclr(i,j,k) = gasopdep(i,j,k,np,ng) + &
rayopdep(i,j,k,nband) + &
aeroextopdep(i,j,k,nband)
ssalbclr(i,j,k) = sctopdepclr(i,j,k)/ &
extopdepclr(i,j,k)
taustrclr(i,j,k) = extopdepclr(i,j,k)*( 1.0 - &
ssalbclr(i,j,k)*fclr(i,j,k) )
omegastrclr(i,j,k) = &
ssalbclr(i,j,k)*((1.0 - fclr(i,j,k))/ &
(1.0 - ssalbclr(i,j,k)*fclr(i,j,k)))
endif
end do
end do
end do
endif
!--------------------------------------------------------------------
! calculate the scaled single-scattering quantities for use in the
! delta-eddington routine.
!--------------------------------------------------------------------
do k=KSRAD,KERAD
do j=JSRAD,JERAD
do i=ISRAD,IERAD
if (cloud(i,j,k) ) then
extopdepovc(i,j,k) = gasopdep(i,j,k,np,ng) + &
rayopdep(i,j,k,nband) + &
aeroextopdep(i,j,k,nband) + &
cloudextopdep(i,j,k)
ssalbovc(i,j,k) = sctopdepovc(i,j,k) / &
extopdepovc(i,j,k)
taustrovc(i,j,k) = extopdepovc(i,j,k)*( 1.0 - &
ssalbovc(i,j,k)*fovc(i,j,k) )
omegastrovc(i,j,k) = ssalbovc(i,j,k)*( ( 1.0 - &
fovc(i,j,k) )/( 1.0 - &
ssalbovc(i,j,k) * &
fovc(i,j,k) ) )
endif
end do
end do
end do
!---------------------------------------------------------------------
! do calculations for all desired zenith angles.
!---------------------------------------------------------------------
do nz = 1, nzens
if (Rad_control%hires_coszen) then
cosangsolar_p(:,:,ng) = Astro%cosz_p(:,:,nz)
else
cosangsolar_p(:,:,ng) = Astro%cosz(:,:)
endif
where (cosangsolar_p(:,:,:) == 0.0) &
cosangsolar_p(:,:,:) = 1.0
!---------------------------------------------------------------------
! clear sky mode
! note: in this mode, the delta-eddington method is performed for all
! spatial columns experiencing sunlight.
!--------------------------------------------------------------------
!--------------------------------------------------------------------
! calculate the scaled single-scattering quantities for use in the
! delta-eddington routine.
!---------------------------------------------------------------------
if (nband >= firstrayband ) then
!---------------------------------------------------------------------
! do diffuse calculation only for first zenith angle -- it is
! independent of zenith angle.
!---------------------------------------------------------------------
if (nz == 1) then
call deledd &
(ix, jx, kx, taustrclr, omegastrclr, &
gstrclr, cosangsolar_p(:,:,ng), ng, daylight, &
rlayerdirclr, tlayerdirclr, tlayerdeclr, &
rlayerdif=rlayerdifclr, tlayerdif=tlayerdifclr)
else
call deledd &
(ix, jx, kx, taustrclr, omegastrclr, &
gstrclr, cosangsolar_p(:,:,ng), ng, daylight,&
rlayerdirclr, tlayerdirclr, tlayerdeclr)
endif
!---------------------------------------------------------------------
! the following needs to be done at daylight points only -- currently
! this code is not active, since ng == 1.
!---------------------------------------------------------------------
if (ng /= 1) then
tlayerdifclr = 0.0
if (NSTREAMS == 1) then
tlayerdifclr(:,:,:) = &
exp(-gasopdep(:,:,:,np,ng)/cosangstr(1))
else
do ns = 1,NSTREAMS
tlayerdifclr(:,:,:) = &
tlayerdifclr(:,:,:) + &
exp( -gasopdep(:,:,:,np,ng)/&
cosangstr(ns) )*wtstr(ns)* &
cosangstr(ns)
end do
endif
rlayerdifclr = 0.0
endif
!---------------------------------------------------------------------
! initialize the layer reflection and transmission arrays for the
! non-rayleigh-scattering case.
!-------------------------------------------------------------------
else
!---------------------------------------------------------------------
! the following needs to be done at daylight points only -- currently
! this code is not active, since ng == 1, and all bands see rayleigh
! scattering.
!---------------------------------------------------------------------
tlayerdifclr = 0.0
if (NSTREAMS == 1) then
tlayerdifclr(:,:,:) = &
exp( -gasopdep(:,:,:,np,ng)/cosangstr(1))
else
do ns = 1,NSTREAMS
tlayerdifclr(:,:,:) = &
tlayerdifclr(:,:,:) + &
exp(-gasopdep(:,:,:,np,ng)/&
cosangstr(ns))*wtstr(ns)*cosangstr(ns)
end do
endif
rlayerdirclr(:,:,:) = 0.0
do k=KSRAD,KERAD
tlayerdirclr(:,:,k) = &
exp( -gasopdep(:,:,k,np,ng) / &
cosangsolar_p(:,:,ng) )
end do
tlayerdeclr(:,:,:) = tlayerdirclr(:,:,:)
rlayerdifclr = 0.0
endif
!---------------------------------------------------------------------
! overcast sky mode
! note: in this mode, the delta-eddington method is performed only
! for spatial columns containing a cloud and experiencing sunlight.
!---------------------------------------------------------------------
!----------------------------------------------------------------------
! calculate the reflection and transmission in the scattering layers
! using the delta-eddington method.
!-------------------------------------------------------------------
if (nz == 1) then
call deledd &
(ix, jx, kx, taustrovc, omegastrovc, gstrovc, &
cosangsolar_p(:,:,ng), ng, daylight, &
rlayerdirovc, tlayerdirovc, tlayerdeovc, &
rlayerdif=rlayerdifovc, tlayerdif=tlayerdifovc,&
cloud=cloud)
else
call deledd &
(ix, jx, kx, taustrovc, omegastrovc, gstrovc, &
cosangsolar_p(:,:,ng), ng, daylight, &
rlayerdirovc, tlayerdirovc, tlayerdeovc, &
cloud=cloud)
endif
if (ng /= 1) then
tlayerdifovc(:,:,:) = tlayerdifclr(:,:,:)
rlayerdifovc(:,:,:) = rlayerdifclr(:,:,:)
endif
!--------------------------------------------------------------------
! weight the reflection and transmission arrays for clear and
! overcast sky conditions by the cloud fraction, to calculate the
! resultant values.
!----------------------------------------------------------------------
do k=KSRAD,KERAD
do j=JSRAD,JERAD
do i=ISRAD,IERAD
if ( cloud(i,j,k) ) then
rlayerdir(i,j,k) = cloudfrac(i,j,k)* &
rlayerdirovc(i,j,k) + &
(1.0 - cloudfrac(i,j,k) )* &
rlayerdirclr(i,j,k)
rlayerdif(i,j,k) = cloudfrac(i,j,k) * &
rlayerdifovc(i,j,k) + &
( 1.0 - cloudfrac(i,j,k) )* &
rlayerdifclr(i,j,k)
tlayerdir(i,j,k) = cloudfrac(i,j,k) * &
tlayerdirovc(i,j,k) + &
( 1.0 - cloudfrac(i,j,k) )* &
tlayerdirclr(i,j,k)
tlayerdif(i,j,k) = cloudfrac(i,j,k) * &
tlayerdifovc(i,j,k) + &
( 1.0 - cloudfrac(i,j,k) )* &
tlayerdifclr(i,j,k)
tlayerde(i,j,k) = cloudfrac(i,j,k) * &
tlayerdeovc(i,j,k) + &
(1.0 - cloudfrac(i,j,k) )* &
tlayerdeclr(i,j,k)
else if (daylight(i,j)) then
rlayerdir(i,j,k) = rlayerdirclr(i,j,k)
tlayerdir(i,j,k) = tlayerdirclr(i,j,k)
rlayerdif(i,j,k) = rlayerdifclr(i,j,k)
tlayerdif(i,j,k) = tlayerdifclr(i,j,k)
tlayerde (i,j,k) = tlayerdeclr (i,j,k)
endif
end do
end do
end do
!---------------------------------------------------------------------
! define the surface albedo (infrared value for infrared bands,
! visible value for the remaining bands).
!----------------------------------------------------------------------c
if (nband <= NIRBANDS ) then
sfcalb_dir(:,:) = Surface%asfc_nir_dir(:,:)
sfcalb_dif(:,:) = Surface%asfc_nir_dif(:,:)
else
sfcalb_dir(:,:) = Surface%asfc_vis_dir(:,:)
sfcalb_dif(:,:) = Surface%asfc_vis_dif(:,:)
end if
!--------------------------------------------------------------------
! calculate the reflection and transmission at flux levels from the
! direct and diffuse values of reflection and transmission in the
! corresponding layers using the adding method.
!---------------------------------------------------------------------
call adding &
(ix, jx, kx, rlayerdir, tlayerdir, rlayerdif, &
tlayerdif, tlayerde, sfcalb_dir, sfcalb_dif, &
daylight, reflectance, transmittance, tr_dir)
!---------------------------------------------------------------------
!
!---------------------------------------------------------------------
if (Rad_control%do_totcld_forcing) then
call adding &
(ix, jx, kx, rlayerdirclr, tlayerdirclr, &
rlayerdifclr, tlayerdifclr, tlayerdeclr, &
sfcalb_dir, sfcalb_dif, cloud_in_column, &
reflectanceclr, transmittanceclr, tr_dirclr)
endif
!----------------------------------------------------------------------
! weight and sum the reflectance and transmittance to calculate the
! band values.
!-------------------------------------------------------------------
do j=JSRAD,JERAD
do i=ISRAD,IERAD
wtfac_p(i,j) = wtfreq(np)*gausswt(ng)* &
cosangsolar_p(i,j,ng)
end do
end do
!---------------------------------------------------------------------
!
!---------------------------------------------------------------------
do k = KSRAD,KERAD+1
do j=JSRAD,JERAD
do i=ISRAD,IERAD
if (daylight(i,j) ) then
sumtr(i,j,k,nz) = sumtr(i,j,k,nz) + &
transmittance(i,j,k)*wtfac_p(i,j)
sumtr_dir(i,j,k,nz) = sumtr_dir(i,j,k,nz) + &
tr_dir(i,j,k)*wtfac_p(i,j)
sumre(i,j,k,nz) = sumre(i,j,k,nz) + &
reflectance(i,j,k)*wtfac_p(i,j)
endif
end do
end do
end do
do j=JSRAD,JERAD
do i=ISRAD,IERAD
if (daylight(i,j) ) then
sumtr_dir_up(i,j,nz) = sumtr_dir_up(i,j,nz) + &
tr_dir(i,j,KERAD+1)*sfcalb_dir(i,j)*wtfac_p(i,j)
endif
end do
end do
!---------------------------------------------------------------------
!
!---------------------------------------------------------------------
if (Rad_control%do_totcld_forcing) then
do k = KSRAD,KERAD+1
do j=JSRAD,JERAD
do i=ISRAD,IERAD
if (cloud_in_column(i,j)) then
sumtrclr(i,j,k,nz) = &
sumtrclr(i,j,k,nz) + &
transmittanceclr(i,j,k)* wtfac_p(i,j)
sumtr_dir_clr(i,j,k,nz) = &
sumtr_dir_clr(i,j,k,nz) + &
tr_dirclr(i,j,k)*wtfac_p(i,j)
sumreclr(i,j,k,nz) = sumreclr(i,j,k,nz) + &
reflectanceclr(i,j,k)*wtfac_p(i,j)
else if (daylight(i,j) ) then
sumtrclr(i,j,k,nz) = sumtrclr(i,j,k,nz) + &
transmittance(i,j,k)*wtfac_p(i,j)
sumtr_dir_clr(i,j,k,nz) = &
sumtr_dir_clr(i,j,k,nz) + tr_dir(i,j,k)*&
wtfac_p(i,j)
sumreclr(i,j,k,nz) = sumreclr(i,j,k,nz) + &
reflectance(i,j,k)*wtfac_p(i,j)
endif
end do
end do
end do
endif
end do ! end of nz loop
end do ! end of gaussian loop
end do ! end of frequency points in the band loop
!----------------------------------------------------------------------
! normalize the solar flux in the band to the appropriate value for
! the given total solar insolation.
!---------------------------------------------------------------------
do nz = 1,nzens
if (Rad_control%hires_coszen) then
fracday_p(:,:) = Astro%fracday_p(:,:,nz)
else
fracday_p(:,:) = Astro%fracday(:,:)
endif
solarflux_p(:,:) = fracday_p(:,:)* &
Solar_spect%solflxband(nband)* &
ssolar/solflxtotal_local
if (nband == Solar_spect%visible_band_indx) then
Sw_output%bdy_flx(:,:,1,nz) = &
Sw_output%bdy_flx(:,:,1,nz) + sumre(:,:,1,nz)* &
solarflux_p(:,:)
Sw_output%bdy_flx(:,:,3,nz) = &
Sw_output%bdy_flx(:,:,3,nz) + sumtr(:,:,KERAD+1,nz)*&
solarflux_p(:,:) - &
sumre(:,:,KERAD+1,nz)* &
solarflux_p(:,:)
endif
if (nband == onepsix_band_indx) then
Sw_output%bdy_flx(:,:,2,nz) = &
Sw_output%bdy_flx(:,:,2,nz) + sumre(:,:,1,nz)* &
solarflux_p(:,:)
Sw_output%bdy_flx(:,:,4,nz) = &
Sw_output%bdy_flx(:,:,4,nz) + sumtr(:,:,KERAD+1,nz)*&
solarflux_p(:,:) - &
sumre(:,:,KERAD+1,nz)*&
solarflux_p(:,:)
endif
!-------------------------------------------------------------------
! calculate the band fluxes and heating rates.
!--------------------------------------------------------------------
if (do_esfsw_band_diagnostics) then
do k = KSRAD,KERAD+1
do j=JSRAD,JERAD
do i=ISRAD,IERAD
dfswband(i,j,k,nz) = sumtr(i,j,k,nz)* &
solarflux_p(i,j)
ufswband(i,j,k,nz) = sumre(i,j,k,nz)* &
solarflux_p(i,j)
end do
end do
end do
endif
!----------------------------------------------------------------------
! sum the band fluxes and heating rates to calculate the total
! spectral values.
!------------------------------------------------------------------
do k = KSRAD,KERAD+1
do j=JSRAD,JERAD
do i=ISRAD,IERAD
if (daylight(i,j) ) then
Sw_output%dfsw (i,j,k,nz) = &
Sw_output%dfsw(i,j,k,nz) + sumtr(i,j,k,nz)*&
solarflux_p(i,j)
Sw_output%ufsw (i,j,k,nz) = &
Sw_output%ufsw(i,j,k,nz) + sumre(i,j,k,nz)* &
solarflux_p(i,j)
fswband(i,j,k,nz) = ((sumre(i,j,k,nz)* &
solarflux_p(i,j)) - &
(sumtr(i,j,k,nz)* &
solarflux_p(i,j)))
Sw_output%fsw(i,j,k,nz) = Sw_output%fsw(i,j,k,nz) +&
fswband(i,j,k,nz)
endif
end do
end do
end do
do j=JSRAD,JERAD
do i=ISRAD,IERAD
if (daylight(i,j) ) then
Sw_output%dfsw_dir_sfc(i,j,nz) = &
Sw_output%dfsw_dir_sfc(i,j,nz) + &
sumtr_dir(i,j,KERAD+1,nz)*solarflux_p(i,j)
Sw_output%ufsw_dir_sfc(i,j,nz) = &
Sw_output%ufsw_dir_sfc(i,j,nz) + &
sumtr_dir_up(i,j,nz)*solarflux_p(i,j)
Sw_output%ufsw_dif_sfc(i,j,nz) = &
Sw_output%ufsw_dif_sfc(i,j,nz) + &
sumre(i,j,KERAD+1,nz)*solarflux_p(i,j)
endif
end do
end do
if (nband > NIRBANDS) then
do j=JSRAD,JERAD
do i=ISRAD,IERAD
if (daylight(i,j) ) then
Sw_output%dfsw_vis_sfc(i,j,nz) = &
Sw_output%dfsw_vis_sfc(i,j,nz) + &
sumtr(i,j,KERAD+1,nz)*solarflux_p(i,j)
Sw_output%ufsw_vis_sfc(i,j,nz) = &
Sw_output%ufsw_vis_sfc(i,j,nz) + &
sumre(i,j,KERAD+1,nz)*solarflux_p(i,j)
Sw_output%dfsw_vis_sfc_dir(i,j,nz) = &
Sw_output%dfsw_vis_sfc_dir(i,j,nz) + &
sumtr_dir(i,j,KERAD+1,nz)*solarflux_p(i,j)
Sw_output%ufsw_vis_sfc_dir(i,j,nz) = &
Sw_output%ufsw_vis_sfc_dir(i,j,nz) + &
sumtr_dir_up(i,j,nz)*solarflux_p(i,j)
Sw_output%ufsw_vis_sfc_dif(i,j,nz) = &
Sw_output%ufsw_vis_sfc_dif(i,j,nz) + &
sumre(i,j,KERAD+1,nz)*solarflux_p(i,j)
endif
end do
end do
endif
!---------------------------------------------------------------------
!
!---------------------------------------------------------------------
do k = KSRAD,KERAD
do j=JSRAD,JERAD
do i=ISRAD,IERAD
if (daylight(i,j) ) then
hswband(i,j,k,nz) = (fswband(i,j,k+1,nz) - &
fswband(i,j,k,nz) )*gocpdp(i,j,k)
Sw_output%hsw(i,j,k,nz) = &
Sw_output%hsw(i,j,k,nz) + hswband(i,j,k,nz)
endif
end do
end do
end do
!----------------------------------------------------------------------
! calculate the band fluxes and heating rates.
!---------------------------------------------------------------------
if (nprofile == 1) then ! clr sky need be done only for
! first cloud profile
if (Rad_control%do_totcld_forcing) then
do j=JSRAD,JERAD
do i=ISRAD,IERAD
if (daylight(i,j) ) then
Sw_output%dfsw_dir_sfc_clr(i,j,nz) = &
Sw_output%dfsw_dir_sfc_clr(i,j,nz) + &
sumtr_dir_clr(i,j,KERAD+1,nz)*solarflux_p(i,j)
endif
end do
end do
if (nband > NIRBANDS) then
do j=JSRAD,JERAD
do i=ISRAD,IERAD
if (daylight(i,j) ) then
Sw_output%dfsw_vis_sfc_clr(i,j,nz) = &
Sw_output%dfsw_vis_sfc_clr(i,j,nz) + &
sumtrclr(i,j,KERAD+1,nz)*solarflux_p(i,j)
endif
end do
end do
endif
if (nband == Solar_spect%visible_band_indx) then
Sw_output%bdy_flx_clr(:,:,1,nz) = &
sumreclr(:,:,1,nz)*solarflux_p(:,:)
Sw_output%bdy_flx_clr(:,:,3,nz) = &
sumtrclr(:,:,KERAD+1,nz)*solarflux_p(:,:) - &
sumreclr(:,:,KERAD+1,nz)*solarflux_p(:,:)
endif
if (nband == onepsix_band_indx) then
Sw_output%bdy_flx_clr(:,:,2,nz) = &
sumreclr(:,:,1,nz)*solarflux_p(:,:)
Sw_output%bdy_flx_clr(:,:,4,nz) = &
sumtrclr(:,:,KERAD+1,nz)*solarflux_p(:,:) - &
sumreclr(:,:,KERAD+1,nz)*solarflux_p(:,:)
endif
if (do_esfsw_band_diagnostics) then
do k = KSRAD,KERAD+1
do j=JSRAD,JERAD
do i=ISRAD,IERAD
dfswbandclr(i,j,k,nz) = &
sumtrclr(i,j,k,nz)*solarflux_p(i,j)
ufswbandclr(i,j,k,nz) = &
sumreclr(i,j,k,nz)*solarflux_p(i,j)
end do
end do
end do
endif
!----------------------------------------------------------------------c
! sum the band fluxes and heating rates to calculate the total
! spectral values.
!----------------------------------------------------------------------c
do k = KSRAD,KERAD+1
do j=JSRAD,JERAD
do i=ISRAD,IERAD
Sw_output%dfswcf(i,j,k,nz) = &
Sw_output%dfswcf(i,j,k,nz) + &
sumtrclr(i,j,k,nz)*solarflux_p(i,j)
Sw_output%ufswcf(i,j,k,nz) = &
Sw_output%ufswcf(i,j,k,nz) + &
sumreclr(i,j,k,nz)*solarflux_p(i,j)
fswbandclr(i,j,k,nz) = &
(sumreclr(i,j,k,nz)*solarflux_p(i,j)) - &
(sumtrclr(i,j,k,nz)*solarflux_p(i,j))
Sw_output%fswcf(i,j,k,nz) = &
Sw_output%fswcf(i,j,k,nz) + &
fswbandclr(i,j,k,nz)
end do
end do
end do
!----------------------------------------------------------------------c
! sum the band fluxes and heating rates to calculate the total
! spectral values.
!----------------------------------------------------------------------c
do k = KSRAD,KERAD
do j=JSRAD,JERAD
do i=ISRAD,IERAD
hswbandclr(i,j,k,nz) = &
(fswbandclr(i,j,k+1,nz) - &
fswbandclr(i,j,k,nz))*gocpdp(i,j,k)
Sw_output%hswcf(i,j,k,nz) = &
Sw_output%hswcf(i,j,k,nz) + &
hswbandclr(i,j,k,nz)
end do
end do
end do
endif
endif ! (nprofile == 1)
end do ! (nz loop)
end do ! (end of band loop)
end do ! (end of nprofile loop)
!----------------------------------------------------------------------
! if the ica calculation was being done, the fluxes and heating rates
! which have been summed over nprofiles cloud profiles must be
! averaged.
!------------------------------------------------------------------
if (do_ica_calcs) then
do j=JSRAD,JERAD
do i=ISRAD,IERAD
Sw_output%dfsw_dir_sfc (i,j,:) = &
Sw_output%dfsw_dir_sfc(i,j,:)*profiles_inverse
Sw_output%ufsw_dif_sfc (i,j,:) = &
Sw_output%ufsw_dif_sfc(i,j,:)*profiles_inverse
Sw_output%dfsw_vis_sfc (i,j,:) = &
Sw_output%dfsw_vis_sfc(i,j,:)*profiles_inverse
Sw_output%ufsw_vis_sfc (i,j,:) = &
Sw_output%ufsw_vis_sfc(i,j,:)*profiles_inverse
Sw_output%dfsw_vis_sfc_dir (i,j,:) = &
Sw_output%dfsw_vis_sfc_dir(i,j,:)*profiles_inverse
Sw_output%ufsw_vis_sfc_dif (i,j,:) = &
Sw_output%ufsw_vis_sfc_dif(i,j,:)*profiles_inverse
Sw_output%bdy_flx(i,j,:,:) = &
Sw_output%bdy_flx(i,j,:,:)*profiles_inverse
end do
end do
do k = KSRAD,KERAD+1
do j=JSRAD,JERAD
do i=ISRAD,IERAD
Sw_output%dfsw (i,j,k,:) = Sw_output%dfsw(i,j,k,:)* &
profiles_inverse
Sw_output%ufsw (i,j,k,:) = Sw_output%ufsw(i,j,k,:)* &
profiles_inverse
Sw_output%fsw(i,j,k,:) = Sw_output%fsw(i,j,k,:)* &
profiles_inverse
end do
end do
end do
do k = KSRAD,KERAD
do j=JSRAD,JERAD
do i=ISRAD,IERAD
Sw_output%hsw(i,j,k,:) = Sw_output%hsw(i,j,k,:)* &
profiles_inverse
end do
end do
end do
endif
do j=JSRAD,JERAD
do i=ISRAD,IERAD
if (daylight(i,j) ) then
Sw_output%dfsw_dif_sfc(i,j,: ) = &
Sw_output%dfsw(i,j,KERAD+1,: ) - &
Sw_output%dfsw_dir_sfc(i,j,: )
endif
end do
end do
if (Rad_control%do_totcld_forcing) then
do j=JSRAD,JERAD
do i=ISRAD,IERAD
if (daylight(i,j) ) then
Sw_output%dfsw_dif_sfc_clr(i,j,: ) = &
Sw_output%dfswcf(i,j,KERAD+1,:) - &
Sw_output%dfsw_dir_sfc_clr(i,j,:)
endif
end do
end do
endif
do j=JSRAD,JERAD
do i=ISRAD,IERAD
if (daylight(i,j) ) then
Sw_output%dfsw_vis_sfc_dif(i,j,:) = &
Sw_output%dfsw_vis_sfc(i,j,:) - &
Sw_output%dfsw_vis_sfc_dir(i,j,:)
endif
end do
end do
!--------------------------------------------------------------------
! convert sw fluxes to cgs and then back to mks units.
!---------------------------------------------------------------------
Sw_output%fsw(:,:,:,:) = &
1.0E-03*(1.0E+03*Sw_output%fsw(:,:,:,:))
Sw_output%dfsw(:,:,:,:) = &
1.0E-03*(1.0E+03*Sw_output%dfsw(:,:,:,:))
Sw_output%ufsw(:,:,:,:) = &
1.0E-03*(1.0E+03*Sw_output%ufsw(:,:,:,:))
if (Rad_control%do_totcld_forcing) then
Sw_output%fswcf(:,:,:,:) = &
1.0E-03*(1.0E+03*Sw_output%fswcf(:,:,:,:))
Sw_output%dfswcf(:,:,:,:) = &
1.0E-03*(1.0E+03*Sw_output%dfswcf(:,:,:,:))
Sw_output%ufswcf(:,:,:,:) = &
1.0E-03*(1.0E+03*Sw_output%ufswcf(:,:,:,:))
endif
!---------------------------------------------------------------------
end subroutine swresf
!####################################################################
subroutine esfsw_driver_end
!---------------------------------------------------------------------
! esfsw_driver_end is the destructor for esfsw_driver_mod.
!---------------------------------------------------------------------
!---------------------------------------------------------------------
! be sure module has been initialized.
!---------------------------------------------------------------------
if (.not. module_is_initialized ) then
call error_mesg ('esfsw_driver_mod', &
'module has not been initialized', FATAL )
endif
!--------------------------------------------------------------------
! close out the modules that this module initialized.
!--------------------------------------------------------------------
call esfsw_parameters_end
!---------------------------------------------------------------------
! mark the module as uninitialized.
!---------------------------------------------------------------------
module_is_initialized = .false.
!---------------------------------------------------------------------
end subroutine esfsw_driver_end
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!
! PRIVATE SUBROUTINES
!
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!#################################################################
!
!
! Subroutine that uses the delta-eddington technique in conjunction
! with a multi-band parameterization for h2o+co2+o2+o3 absorption
! in the solar spectrum to derive solar fluxes and heating rates.
!
!
! This subroutine calculates optical depth, single scattering albedo,
! asymmetry parameter of a layer based on gaseous absorbers,
! clouds, aerosols, and rayleigh scattering. It then uses delta-
! eddington technique to calculate radiative flux at each layer.
! Doubling and adding technique is used to combine the layers
! and calculate flux at TOA and surface and heating rate. This
! subroutine allocates a substantial amount of memory and deallocates
! the allocated memory at the end of the subroutine.
!
!
! call comput(is, ie, js, je, Atmos_input, Surface, Rad_gases, Aerosol,
! Astro, &
! Cldrad_props, Cld_spec, Sw_output)
!
!
! starting subdomain i indices of data in
! the physics_window being integrated
!
!
! ending subdomain i indices of data in
! the physics_window being integrated
!
!
! starting subdomain j indices of data in
! the physics_window being integrated
!
!
! ending subdomain j indices of data in
! the physics_window being integrated
!
!
! Atmos_input_type variable containing the atmospheric
! input fields on the radiation grid
!
!
! Aerosol input data for shortwave radiation calculation
!
!
! Astronomy_type variable containing the astronomical
! input fields on the radiation grid
!
!
! Radiative_gases_type variable containing the radiative
! gas input fields on the radiation grid
!
!
! The cloud radiative property input fields on the
! radiation grid
!
!
! The shortwave radiation calculation result
!
!
! Surface data as boundary condition to radiation
!
!
! Cloud specification data as initial condition to radiation
!
!
subroutine compute_aerosol_optical_props &
(Atmos_input, Aerosol, Aerosol_props, including_volcanoes, &
Aerosol_diags, r, including_aerosols, naerosol_optical, &
daylight, aeroextopdep, aerosctopdep, aeroasymfac)
!----------------------------------------------------------------------
! comput uses the delta-eddington technique in conjunction with a
! multiple-band parameterization for h2o+co2+o2+o3 absorption to
! derive solar fluxes and heating rates.
! notes: drops are assumed if temp>273.15K, ice crystals otherwise.
!-------------------------------------------------------------------
type(atmos_input_type), intent(in) :: Atmos_input
type(aerosol_type), intent(in) :: Aerosol
type(aerosol_properties_type), intent(in) :: Aerosol_props
logical, intent(in) :: including_volcanoes
type(aerosol_diagnostics_type),intent(inout) :: Aerosol_diags
real,dimension(:,:,:,:), intent(inout) :: r
logical, intent(in) :: including_aerosols
integer, intent(in) :: naerosol_optical
logical,dimension(:,:), intent(in) :: daylight
real, dimension(:,:,:,:), intent(out) :: aeroextopdep, &
aerosctopdep, &
aeroasymfac
!-------------------------------------------------------------------
! intent(in) variables:
!
! Atmos_input atmos_input_type structure, contains variables
! defining atmospheric state
! Aerosol aerosol_type structure, contains variables
! defining aerosol fields, passed through to
! lower level routines
! Aerosol_props aerosol radiative property input data for the
! radiation package
!
! intent(inout) variables:
!
! Sw_output shortwave radiation output data
!
!---------------------------------------------------------------------
!-----------------------------------------------------------------------
! local variables:
real, dimension (size(Atmos_input%temp,3)-1) :: &
arprod, asymm, arprod2, deltaz, &
sum_g_omega_tau, sum_ext, sum_sct
integer, dimension (size (Atmos_input%press, 3)-1 ) :: &
opt_index_v3, opt_index_v4, &
opt_index_v5, opt_index_v6, &
opt_index_v7, opt_index_v8, &
opt_index_v9, opt_index_v10, &
irh
real, dimension (naerosol_optical) :: &
aerext, aerssalb, aerasymm
real :: aerext_i, aerssalb_i, aerasymm_i
integer :: j, i, k, nband, nsc
integer :: israd, jsrad, ierad, jerad, ksrad, kerad
integer :: nextinct !variable to pass extinction to chemistry
!-----------------------------------------------------------------------
! local variables:
!
! aeramt
! sum_g_omega_tau
! opt_index_v3
! irh
! etc.
!
!--------------------------------------------------------------------
!--------------------------------------------------------------------
! define limits and dimensions
!--------------------------------------------------------------------
israd = 1
jsrad = 1
ksrad = 1
ierad = size(Atmos_input%temp,1)
jerad = size(Atmos_input%temp,2)
kerad = size(Atmos_input%temp,3) - 1
naerosoltypes_used = size(Aerosol%aerosol,4)
do j = JSRAD,JERAD
do i = ISRAD,IERAD
if (daylight(i,j) .or. Sw_control%do_cmip_diagnostics) then
deltaz(:) = Atmos_input%deltaz(i,j,:)
do nband = 1,Solar_spect%nbands
if (including_aerosols) then
aerext(:) = Aerosol_props%aerextband(nband,:)
aerssalb(:) = Aerosol_props%aerssalbband(nband,:)
aerasymm(:) = Aerosol_props%aerasymmband(nband,:)
!-------------------------------------------------------------------
! define the local variables for the band values of aerosol and
! cloud single scattering parameters.
! note: the unit for the aerosol extinction is kilometer**(-1).
!--------------------------------------------------------------------
do k = KSRAD,KERAD
irh(k) = MIN(100, MAX( 0, &
NINT(100.*Atmos_input%aerosolrelhum(i,j,k))))
opt_index_v3(k) = &
Aerosol_props%sulfate_index (irh(k), &
Aerosol_props%ivol(i,j,k))
opt_index_v4(k) = &
Aerosol_props%omphilic_index( irh(k) )
opt_index_v5(k) = &
Aerosol_props%bcphilic_index( irh(k) )
opt_index_v6(k) = &
Aerosol_props%seasalt1_index( irh(k) )
opt_index_v7(k) = &
Aerosol_props%seasalt2_index( irh(k) )
opt_index_v8(k) = &
Aerosol_props%seasalt3_index( irh(k) )
opt_index_v9(k) = &
Aerosol_props%seasalt4_index( irh(k) )
opt_index_v10(k) = &
Aerosol_props%seasalt5_index( irh(k) )
end do
!---------------------------------------------------------------------
! calculate scattering properties for all aerosol constituents
! combined.
!---------------------------------------------------------------------
do k = KSRAD,KERAD
sum_g_omega_tau(k) = 0.0
sum_ext(k) = 0.
sum_sct(k) = 0.
end do
do nsc = 1,NAEROSOLTYPES_USED
if (Aerosol_props%optical_index(nsc) > 0) then
aerext_i = &
aerext(Aerosol_props%optical_index(nsc))
aerssalb_i = &
aerssalb(Aerosol_props%optical_index(nsc))
aerasymm_i = &
aerasymm(Aerosol_props%optical_index(nsc))
do k = KSRAD,KERAD
arprod(k) = &
aerext_i*(1.e3*Aerosol%aerosol(i,j,k,nsc))
arprod2(k) = aerssalb_i*arprod(k)
asymm(k) = aerasymm_i
sum_ext(k) = sum_ext(k) + arprod(k)
sum_sct(k) = sum_sct(k) + aerssalb_i*arprod(k)
sum_g_omega_tau(k) = sum_g_omega_tau(k) + &
aerasymm_i*(aerssalb_i*arprod(k))
end do
else if (Aerosol_props%optical_index(nsc) == &
Aerosol_props%sulfate_flag) then
do k = KSRAD,KERAD
arprod(k) = aerext(opt_index_v3(k)) * &
(1.e3 * Aerosol%aerosol(i,j,k,nsc))
arprod2(k) = aerssalb(opt_index_v3(k))*arprod(k)
asymm(k) = aerasymm(opt_index_v3(k))
sum_ext(k) = sum_ext(k) + arprod(k)
sum_sct(k) = sum_sct(k) + &
aerssalb(opt_index_v3(k))*arprod(k)
sum_g_omega_tau(k) = sum_g_omega_tau(k) + &
aerasymm(opt_index_v3(k))* &
(aerssalb(opt_index_v3(k))*arprod(k))
end do
else if (Aerosol_props%optical_index(nsc) == &
Aerosol_props%bc_flag) then
do k = KSRAD,KERAD
arprod(k) = aerext(opt_index_v3(k)) * &
(1.e3 * Aerosol%aerosol(i,j,k,nsc))
arprod2(k) = aerssalb(opt_index_v3(k))*arprod(k)
asymm(k) = aerasymm(opt_index_v3(k))
sum_ext(k) = sum_ext(k) + arprod(k)
sum_sct(k) = sum_sct(k) + &
aerssalb(opt_index_v3(k))*arprod(k)
sum_g_omega_tau(k) = sum_g_omega_tau(k) + &
aerasymm(opt_index_v3(k)) * &
(aerssalb(opt_index_v3(k))*arprod(k))
end do
else if (Aerosol_props%optical_index(nsc) == &
Aerosol_props%omphilic_flag) then
do k = KSRAD,KERAD
arprod(k) = aerext(opt_index_v4(k)) * &
(1.e3 * Aerosol%aerosol(i,j,k,nsc))
arprod2(k) = aerssalb(opt_index_v4(k))*arprod(k)
asymm(k) = aerasymm(opt_index_v4(k))
sum_ext(k) = sum_ext(k) + arprod(k)
sum_sct(k) = sum_sct(k) + &
aerssalb(opt_index_v4(k))*arprod(k)
sum_g_omega_tau(k) = sum_g_omega_tau(k) + &
aerasymm(opt_index_v4(k))* &
(aerssalb(opt_index_v4(k))*arprod(k))
end do
else if (Aerosol_props%optical_index(nsc) == &
Aerosol_props%bcphilic_flag) then
if (Rad_control%using_im_bcsul) then
do k = KSRAD,KERAD
arprod(k) = aerext(opt_index_v3(k)) * &
(1.e3 * Aerosol%aerosol(i,j,k,nsc))
arprod2(k) = aerssalb(opt_index_v3(k))*arprod(k)
asymm(k) = aerasymm(opt_index_v3(k))
sum_ext(k) = sum_ext(k) + arprod(k)
sum_sct(k) = sum_sct(k) + &
aerssalb(opt_index_v3(k))*arprod(k)
sum_g_omega_tau(k) = sum_g_omega_tau(k) + &
aerasymm(opt_index_v3(k)) * &
(aerssalb(opt_index_v3(k))*arprod(k))
end do
else ! (using_im_bcsul)
do k = KSRAD,KERAD
arprod(k) = aerext(opt_index_v5(k)) * &
(1.e3 * Aerosol%aerosol(i,j,k,nsc))
arprod2(k) = aerssalb(opt_index_v5(k))*arprod(k)
asymm(k) = aerasymm(opt_index_v5(k))
sum_ext(k) = sum_ext(k) + arprod(k)
sum_sct(k) = sum_sct(k) + &
aerssalb(opt_index_v5(k))*arprod(k)
sum_g_omega_tau(k) = sum_g_omega_tau(k) + &
aerasymm(opt_index_v5(k)) * &
(aerssalb(opt_index_v5(k))*arprod(k))
end do
endif !(using_im_bcsul)
else if (Aerosol_props%optical_index(nsc) == &
Aerosol_props%seasalt1_flag) then
do k = KSRAD,KERAD
arprod(k) = aerext(opt_index_v6(k)) * &
(1.e3 * Aerosol%aerosol(i,j,k,nsc))
arprod2(k) = aerssalb(opt_index_v6(k))*arprod(k)
asymm(k) = aerasymm(opt_index_v6(k))
sum_ext(k) = sum_ext(k) + arprod(k)
sum_sct(k) = sum_sct(k) + &
aerssalb(opt_index_v6(k))*arprod(k)
sum_g_omega_tau(k) = sum_g_omega_tau(k) + &
aerasymm(opt_index_v6(k)) * &
(aerssalb(opt_index_v6(k))*arprod(k))
end do
else if (Aerosol_props%optical_index(nsc) == &
Aerosol_props%seasalt2_flag) then
do k = KSRAD,KERAD
arprod(k) = aerext(opt_index_v7(k)) * &
(1.e3 * Aerosol%aerosol(i,j,k,nsc))
arprod2(k) = aerssalb(opt_index_v7(k))*arprod(k)
asymm(k) = aerasymm(opt_index_v7(k))
sum_ext(k) = sum_ext(k) + arprod(k)
sum_sct(k) = sum_sct(k) + &
aerssalb(opt_index_v7(k))*arprod(k)
sum_g_omega_tau(k) = sum_g_omega_tau(k) + &
aerasymm(opt_index_v7(k)) * &
(aerssalb(opt_index_v7(k))*arprod(k))
end do
else if (Aerosol_props%optical_index(nsc) == &
Aerosol_props%seasalt3_flag) then
do k = KSRAD,KERAD
arprod(k) = aerext(opt_index_v8(k)) * &
(1.e3 * Aerosol%aerosol(i,j,k,nsc))
arprod2(k) = aerssalb(opt_index_v8(k))*arprod(k)
asymm(k) = aerasymm(opt_index_v8(k))
sum_ext(k) = sum_ext(k) + arprod(k)
sum_sct(k) = sum_sct(k) + &
aerssalb(opt_index_v8(k))*arprod(k)
sum_g_omega_tau(k) = sum_g_omega_tau(k) + &
aerasymm(opt_index_v8(k)) * &
(aerssalb(opt_index_v8(k))*arprod(k))
end do
else if (Aerosol_props%optical_index(nsc) == &
Aerosol_props%seasalt4_flag) then
do k = KSRAD,KERAD
arprod(k) = aerext(opt_index_v9(k)) * &
(1.e3 * Aerosol%aerosol(i,j,k,nsc))
arprod2(k) = aerssalb(opt_index_v9(k))*arprod(k)
asymm(k) = aerasymm(opt_index_v9(k))
sum_ext(k) = sum_ext(k) + arprod(k)
sum_sct(k) = sum_sct(k) + &
aerssalb(opt_index_v9(k))*arprod(k)
sum_g_omega_tau(k) = sum_g_omega_tau(k) + &
aerasymm(opt_index_v9(k))* &
(aerssalb(opt_index_v9(k))*arprod(k))
end do
else if (Aerosol_props%optical_index(nsc) == &
Aerosol_props%seasalt5_flag) then
do k = KSRAD,KERAD
arprod(k) = aerext(opt_index_v10(k)) * &
(1.e3 * Aerosol%aerosol(i,j,k,nsc))
arprod2(k) = aerssalb(opt_index_v10(k))*arprod(k)
asymm(k) = aerasymm(opt_index_v10(k))
sum_ext(k) = sum_ext(k) + arprod(k)
sum_sct(k) = sum_sct(k) + &
aerssalb(opt_index_v10(k))*arprod(k)
sum_g_omega_tau(k) = sum_g_omega_tau(k) +&
aerasymm(opt_index_v10(k)) * &
(aerssalb(opt_index_v10(k))*arprod(k))
end do
endif
if (Sw_control%do_cmip_diagnostics) then
if (nband == Solar_spect%visible_band_indx) then
Aerosol_diags%extopdep(i,j,:,nsc,1) = arprod(:)
Aerosol_diags%absopdep(i,j,:,nsc,1) = &
arprod(:) - arprod2(:)
Aerosol_diags%asymdep(i,j,:,nsc,1) = asymm(:)
endif
if (nband == Solar_spect%eight70_band_indx) then
Aerosol_diags%extopdep(i,j,:,nsc,6) = arprod(:)
Aerosol_diags%absopdep(i,j,:,nsc,6) = &
arprod(:) - arprod2(:)
Aerosol_diags%asymdep(i,j,:,nsc,6) = asymm(:)
endif
if (nband == Solar_spect%one_micron_indx) then
Aerosol_diags%extopdep(i,j,:,nsc,2) = arprod(:)
Aerosol_diags%absopdep(i,j,:,nsc,2) = &
arprod(:) - arprod2(:)
Aerosol_diags%asymdep(i,j,:,nsc,2) = asymm(:)
endif
if (nband == Solar_spect%w340_band_indx) then
Aerosol_diags%extopdep(i,j,:,nsc,7) = arprod(:)
Aerosol_diags%absopdep(i,j,:,nsc,7) = &
arprod(:) - arprod2(:)
Aerosol_diags%asymdep(i,j,:,nsc,7) = asymm(:)
endif
if (nband == Solar_spect%w380_band_indx) then
Aerosol_diags%extopdep(i,j,:,nsc,8) = arprod(:)
Aerosol_diags%absopdep(i,j,:,nsc,8) = &
arprod(:) - arprod2(:)
Aerosol_diags%asymdep(i,j,:,nsc,8) = asymm(:)
endif
if (nband == Solar_spect%w440_band_indx) then
Aerosol_diags%extopdep(i,j,:,nsc,9) = arprod(:)
Aerosol_diags%absopdep(i,j,:,nsc,9) = &
arprod(:) - arprod2(:)
Aerosol_diags%asymdep(i,j,:,nsc,9) = asymm(:)
endif
if (nband == Solar_spect%w670_band_indx) then
Aerosol_diags%extopdep(i,j,:,nsc,10) = arprod(:)
Aerosol_diags%absopdep(i,j,:,nsc,10) = &
arprod(:) - arprod2(:)
Aerosol_diags%asymdep(i,j,:,nsc,10) = asymm(:)
endif
endif
end do
!----------------------------------------------------------------------
! add the effects of volcanic aerosols, if they are to be included.
! include generation of diagnostics in the visible (0.55 micron) and
! nir band (1.0 micron).
!----------------------------------------------------------------------
if (including_volcanoes) then
do k = KSRAD,KERAD
sum_ext(k) = sum_ext(k) + &
Aerosol_props%sw_ext(i,j,k,nband)* &
deltaz(k)
sum_sct(k) = sum_sct(k) + &
Aerosol_props%sw_ssa(i,j,k,nband)* &
Aerosol_props%sw_ext(i,j,k,nband)* &
deltaz(k)
sum_g_omega_tau(k) = &
sum_g_omega_tau(k) +&
Aerosol_props%sw_asy(i,j,k,nband)* &
Aerosol_props%sw_ssa(i,j,k,nband)* &
Aerosol_props%sw_ext(i,j,k,nband)* &
deltaz(k)
if (Sw_control%do_cmip_diagnostics) then
if (nband == Solar_spect%visible_band_indx) then
Aerosol_diags%extopdep_vlcno(i,j,k,1) = &
Aerosol_props%sw_ext(i,j,k,nband)* &
deltaz(k)
Aerosol_diags%absopdep_vlcno(i,j,k,1) = &
(1.0 - Aerosol_props%sw_ssa(i,j,k,nband))*&
Aerosol_props%sw_ext(i,j,k,nband)* &
deltaz(k)
endif
if (nband == Solar_spect%eight70_band_indx) then
Aerosol_diags%extopdep_vlcno(i,j,k,3) = &
Aerosol_props%sw_ext(i,j,k,nband)* &
deltaz(k)
Aerosol_diags%absopdep_vlcno(i,j,k,3) = &
(1.0 - Aerosol_props%sw_ssa(i,j,k,nband))*&
Aerosol_props%sw_ext(i,j,k,nband)* &
deltaz(k)
endif
if (nband == Solar_spect%one_micron_indx) then
Aerosol_diags%extopdep_vlcno(i,j,k,2) = &
Aerosol_props%sw_ext(i,j,k,nband)* &
deltaz(k)
Aerosol_diags%absopdep_vlcno(i,j,k,2) = &
(1.0 - Aerosol_props%sw_ssa(i,j,k,nband))*&
Aerosol_props%sw_ext(i,j,k,nband)* &
deltaz(k)
endif
endif
end do
endif ! (including_volcanoes)
!
!----------------------------------------------------------------------
do k = KSRAD,KERAD
aeroextopdep(i,j,k,nband) = sum_ext(k)
aerosctopdep(i,j,k,nband) = sum_sct(k)
aeroasymfac(i,j,k,nband) = sum_g_omega_tau(k) / &
(sum_sct(k) + 1.0e-30 )
end do
else ! (if not including_aerosols)
do k = KSRAD,KERAD
aeroextopdep(i,j,k,nband) = 0.0
aerosctopdep(i,j,k,nband) = 0.0
aeroasymfac(i,j,k,nband) = 0.0
end do
endif ! (including_aerosols)
end do ! (nband)
endif ! (daylight or cmip_diagnostics)
if (including_volcanoes) then
if (do_coupled_stratozone ) then
nextinct = get_tracer_index(MODEL_ATMOS,'Extinction')
if (nextinct /= NO_TRACER) &
r(i,j,:,nextinct) = Aerosol_props%sw_ext(i,j,:,4)
endif
endif
end do ! (i loop)
end do ! (j loop)
!---------------------------------------------------------------------
!
!---------------------------------------------------------------------
end subroutine compute_aerosol_optical_props
!#################################################################
!
!
! Subroutine that uses the delta-eddington technique in conjunction
! with a multi-band parameterization for h2o+co2+o2+o3 absorption
! in the solar spectrum to derive solar fluxes and heating rates.
!
!
! This subroutine calculates optical depth, single scattering albedo,
! asymmetry parameter of a layer based on gaseous absorbers,
! clouds, aerosols, and rayleigh scattering. It then uses delta-
! eddington technique to calculate radiative flux at each layer.
! Doubling and adding technique is used to combine the layers
! and calculate flux at TOA and surface and heating rate. This
! subroutine allocates a substantial amount of memory and deallocates
! the allocated memory at the end of the subroutine.
!
!
! call comput(is, ie, js, je, Atmos_input, Surface, Rad_gases, Aerosol,
! Astro, &
! Cldrad_props, Cld_spec, Sw_output)
!
!
! starting subdomain i indices of data in
! the physics_window being integrated
!
!
! ending subdomain i indices of data in
! the physics_window being integrated
!
!
! starting subdomain j indices of data in
! the physics_window being integrated
!
!
! ending subdomain j indices of data in
! the physics_window being integrated
!
!
! Atmos_input_type variable containing the atmospheric
! input fields on the radiation grid
!
!
! Aerosol input data for shortwave radiation calculation
!
!
! Astronomy_type variable containing the astronomical
! input fields on the radiation grid
!
!
! Radiative_gases_type variable containing the radiative
! gas input fields on the radiation grid
!
!
! The cloud radiative property input fields on the
! radiation grid
!
!
! The shortwave radiation calculation result
!
!
! Surface data as boundary condition to radiation
!
!
! Cloud specification data as initial condition to radiation
!
!
subroutine compute_gas_props (Atmos_input, Rad_gases, Astro, &
daylight, gasopdep)
!----------------------------------------------------------------------
! comput uses the delta-eddington technique in conjunction with a
! multiple-band parameterization for h2o+co2+o2+o3 absorption to
! derive solar fluxes and heating rates.
! notes: drops are assumed if temp>273.15K, ice crystals otherwise.
!-------------------------------------------------------------------
type(atmos_input_type), intent(in) :: Atmos_input
type(radiative_gases_type), intent(in) :: Rad_gases
type(astronomy_type), intent(in) :: Astro
logical, dimension(:,:), intent(in) :: daylight
real, dimension(:,:,:,:,:), intent(out) :: gasopdep
!-------------------------------------------------------------------
! intent(in) variables:
!
! Atmos_input atmos_input_type structure, contains variables
! defining atmospheric state
! Rad_gases radiative_gases_type structure, contains var-
! iables defining the radiatively active gases,
! passed through to lower level routines
! Astro astronomy_type structure
!
! intent(inout) variables:
!
! Sw_output shortwave radiation output data
!
!---------------------------------------------------------------------
!-----------------------------------------------------------------------
! local variables:
real, dimension (size(Atmos_input%temp,3)-1) :: &
deltaz, qo3, rh2o, &
efftauo2, efftauco2, efftauch4, efftaun2o, &
wh2ostr, wo3, wo2, quenchfac, &
opdep, delpdig, deltap, tco2, &
tch4, tn2o, to2, wh2o
real, dimension (size(Atmos_input%temp,3)) :: &
alphaco2, alphaco2str, alphao2, &
alphao2str, alphach4, alphach4str, &
alphan2o, alphan2ostr, scale, &
scalestr, totco2, totco2str, &
toto2, toto2str, totch4, &
totch4str, totn2o, totn2ostr, &
press, pflux, pflux_mks, &
temp, z
real :: cosangsolar
real :: denom
real :: wtquench
real :: rrvco2
real :: rrvch4, rrvn2o
integer :: j, i, k, ng, nband, kq
integer :: np, nf
integer :: israd, jsrad, ierad, jerad, ksrad, kerad
!-----------------------------------------------------------------------
! local variables:
!
! aeramt
! sum_g_omega_tau
! opt_index_v3
! irh
! etc.
!
!--------------------------------------------------------------------
!--------------------------------------------------------------------
! define limits and dimensions
!--------------------------------------------------------------------
israd = 1
jsrad = 1
ksrad = 1
ierad = size(Atmos_input%temp,1)
jerad = size(Atmos_input%temp,2)
kerad = size(Atmos_input%temp,3) - 1
!---------------------------------------------------------------------
! initialize local variables.
!------------------------------------------------------------------
alphaco2 (1) = 0.0
alphaco2str(1) = 0.0
alphao2 (1) = 0.0
alphao2str (1) = 0.0
alphach4 (1) = 0.0
alphach4str(1) = 0.0
alphan2o (1) = 0.0
alphan2ostr(1) = 0.0
rrvco2 = Rad_gases%rrvco2
rrvch4 = Rad_gases%rrvch4
rrvn2o = Rad_gases%rrvn2o
do j = JSRAD,JERAD
do i = ISRAD,IERAD
if ( daylight(i,j) ) then
!---------------------------------------------------------------------
! convert to cgs and then back to mks for consistency with previous
!---------------------------------------------------------------------
do k = KSRAD,KERAD+1
press(k) = 0.1*(10.0*Atmos_input%press(i,j,k))
pflux(k) = (10.0*Atmos_input%pflux(i,j,k))
temp (k) = Atmos_input%temp (i,j,k)
end do
do k = KSRAD,KERAD
rh2o (k) = Atmos_input%rh2o (i,j,k)
qo3 (k) = Rad_gases%qo3(i,j,k)
deltaz(k) = Atmos_input%deltaz(i,j,k)
end do
!----------------------------------------------------------------------c
! define pressure related quantities, pressure is in mks units.
!----------------------------------------------------------------------c
pflux_mks(:) = pflux(:)*1.0E-1
do k = KSRAD+1,KERAD+1
deltap(k-1) = pflux_mks(k) - pflux_mks(k-1)
delpdig(k-1) = deltap(k-1)/ GRAV
scalestr(k) = pflux_mks(k)
scale(k) = scalestr(k)*pflux_mks(k)/pstd_mks
end do
do k = KSRAD,KERAD
wh2ostr(k) = rh2o(k)*delpdig(k)
wo3(k) = qo3(k)*delpdig(k)
wo2(k) = o2mixrat*(WTMO2/WTMAIR)*delpdig(k)
end do
!---------------------------------------------------------------------
! if quenching factor effects are desired, calculate the height above
! the surface of the model flux levels.
!---------------------------------------------------------------------
if (do_quench) then
z(KERAD+1) = 0.0
do k = KERAD,KSRAD,-1
z(k) = z(k+1) + deltaz(k)
end do
!---------------------------------------------------------------------
! define the quenching factor for each grid point.
!---------------------------------------------------------------------
do k = KSRAD,KERAD
if (z(k) < co2_quenchfac_height(1) ) then
quenchfac(k) = 1.0
else if (z(k) > co2_quenchfac_height(30) ) then
quenchfac(k) = 0.037
else
do kq = 1,29
if (z(k) > co2_quenchfac_height(kq) .and. &
z(k) <= co2_quenchfac_height(kq+1)) then
wtquench = (z(k) - co2_quenchfac_height(kq))/ &
(co2_quenchfac_height(kq+1) - &
co2_quenchfac_height(kq))
quenchfac(k) = (1. - wtquench)* &
co2_quenchfac(kq) + &
wtquench*co2_quenchfac(kq+1)
exit
endif
end do
endif
end do
else
quenchfac(:) = 1.0
endif !(do_quench)
do ng = 1,NSOLWG
cosangsolar = Astro%cosz(i,j)
if (cosangsolar == 0.0) cosangsolar = 1.0
!----------------------------------------------------------------------c
! define the scaled and unscaled co2 and o2 pathlengths in
! centimeter-atm, and the unscaled h2o and o3 amounts in
! kgrams/meter**2.
! cm-atm needed as units because of c2co2 having those units.
!----------------------------------------------------------------------c
denom = 1.0/(GRAV*rhoair*cosangsolar*2.0)
do k = KSRAD+1,KERAD+1
totco2(k) = 1.0E+02*rrvco2*scale(k)*denom
totco2str(k) = 2.0E+02*rrvco2*scalestr(k)*denom
toto2(k) = 1.0E+02*o2mixrat*scale(k)*denom
toto2str(k) = 2.0E+02*o2mixrat*scalestr(k)*denom
if (do_ch4_sw_effects) then
totch4(k) = 1.0E+02*rrvch4*scale(k)*denom
totch4str(k) = 2.0E+02*rrvch4*scalestr(k)*denom
endif
if (do_n2o_sw_effects) then
totn2o(k) = 1.0E+02*rrvn2o*scale(k)*denom
totn2ostr(k) = 2.0E+02*rrvn2o*scalestr(k)*denom
endif
end do
np = 0
do nband = 1, NBANDS
!-------------------------------------------------------------------
! define the h2o scaled gas amounts in kgrams/meter**2
!---------------------------------------------------------------------
if (nband <= nh2obands) then
do k = KSRAD,KERAD
wh2o(k) = rh2o(k)*delpdig(k)* &
exp(powph2o(nband)*alog(press(k)*p0h2o(nband)))
end do
!---------------------------------------------------------------------
! calculate the "effective" co2, o2, ch4 and n2o gas optical depths
! for the appropriate absorbing bands.
! note: for large zenith angles, alpha can exceed 1. In this case,a
! the optical depths are set to the previous layer values.
!-------------------------------------------------------------------
if ( c1co2(nband).ne.1.0E-99 ) then
do k = KSRAD+1,KERAD+1
if (totco2(k) < totco2max(nband) .and. &
totco2str(k) < totco2strmax(nband)) then
alphaco2(k) = &
c1co2(nband)*exp(c3co2(nband)* &
alog((totco2(k) + c2co2(nband)))) - &
c4co2(nband)
alphaco2str(k) = &
c1co2str(nband)*exp(c3co2str(nband)* &
alog((totco2str(k) + c2co2str(nband)))) - &
c4co2str(nband)
tco2(k-1) = &
(1.0 - alphaco2(k))* &
(1.0 - alphaco2str(k))/ &
((1.0 - alphaco2(k-1))* &
(1.0 - alphaco2str(k-1)))
efftauco2(k-1) = -cosangsolar*alog( tco2(k-1))
else if (k > KSRAD+1) then
efftauco2(k-1) = efftauco2(k-2)
else
efftauco2(k-1) = 0.0
end if
end do
else !( c1co2(nband).ne.1.0E-99 )
efftauco2(:) = 0.0
end if !( c1co2(nband).ne.1.0E-99 )
if (do_ch4_sw_effects) then
if (c1ch4(nband).ne.1.0E-99 ) then
do k = KSRAD+1,KERAD+1
if (totch4(k) < totch4max(nband) .and. &
totch4str(k) < totch4strmax(nband)) then
alphach4(k) = &
c1ch4(nband)*exp(c3ch4(nband)*&
alog((totch4(k) + c2ch4(nband)))) - &
c4ch4(nband)
alphach4str(k) = &
c1ch4str(nband)*exp(c3ch4str(nband)* &
alog((totch4str(k) + c2ch4str(nband)))) - &
c4ch4str(nband)
tch4(k-1) = &
(1.0 - alphach4(k))* &
(1.0 - alphach4str(k))/ &
((1.0 - alphach4(k-1))* &
(1.0 - alphach4str(k-1)))
efftauch4(k-1) = -cosangsolar*alog(tch4(k-1))
else if (k > KSRAD+1) then
efftauch4(k-1) = efftauch4(k-2)
else
efftauch4(k-1) = 0.0
end if
end do
else !( c1ch4(nband).ne.1.0E-99 )
efftauch4(:) = 0.0
end if !( c1ch4(nband).ne.1.0E-99 )
else !do_ch4 = .false.
efftauch4(:) = 0.0
end if
if (do_n2o_sw_effects) then
if ( c1n2o(nband).ne.1.0E-99 ) then
do k = KSRAD+1,KERAD+1
if (totn2o(k) < totn2omax(nband) .and. &
totn2ostr(k) < totn2ostrmax(nband)) then
alphan2o(k) = &
c1n2o(nband)*exp(c3n2o(nband)* &
alog((totn2o(k) +c2n2o(nband)))) - &
c4n2o(nband)
alphan2ostr(k) = &
c1n2ostr(nband)*exp(c3n2ostr(nband)* &
alog((totn2ostr(k) + c2n2ostr(nband)))) - &
c4n2ostr(nband)
tn2o(k-1) = &
(1.0 - alphan2o(k)) * &
(1.0 - alphan2ostr(k))/ &
(( 1.0 - alphan2o(k-1)) * &
(1.0 - alphan2ostr(k-1)))
efftaun2o(k-1) = -cosangsolar*alog(tn2o(k-1))
else if (k > KSRAD+1) then
efftaun2o(k-1) = efftaun2o(k-2)
else
efftaun2o(k-1) = 0.0
end if
end do
else !( c1n2o(nband).ne.1.0E-99 )
efftaun2o(:) = 0.0
end if !( c1n2o(nband).ne.1.0E-99 )
else !do_n2o = .false.
efftaun2o(:) = 0.0
end if
if ( c1o2(nband).ne.1.0E-99 ) then
do k = KSRAD+1,KERAD+1
if (toto2(k) .lt. toto2max(nband) .and. &
toto2str(k) .lt. toto2strmax(nband)) then
alphao2(k) = c1o2(nband)*exp( c3o2(nband)* &
alog((toto2(k) + c2o2(nband)))) - &
c4o2(nband)
alphao2str( k) = &
c1o2str(nband)*exp(c3o2str(nband)* &
alog((toto2str(k) + c2o2str(nband)))) &
- c4o2str(nband)
to2(k-1) = &
(1.0 - alphao2(k))* &
(1.0 - alphao2str(k) )/ &
((1.0 - alphao2(k-1)) * &
(1.0 - alphao2str(k-1)))
efftauo2(k-1) = -cosangsolar*alog(to2(k-1))
else if (k.gt.KSRAD+1) then
efftauo2(k-1) = efftauo2(k-2)
else
efftauo2(k-1) = 0.0
end if
end do
else ! ( c1o2(nband).ne.1.0E-99 )
efftauo2(:) = 0.0
end if ! ( c1o2(nband).ne.1.0E-99 )
end if ! (nband <= nh2obands)
!---------------------------------------------------------------------
! calculate the "effective" o2 gas optical depths for the Schuman-
! Runge band.
!-------------------------------------------------------------------
if ( nband.EQ.NBANDS ) then
do k = KSRAD+1,KERAD+1
if ( toto2str(k).lt.toto2strmaxschrun) then
alphao2str(k) = &
c1o2strschrun*exp( c3o2strschrun*&
alog((toto2str(k) + c2o2strschrun))) - &
c4o2strschrun
to2(k-1) = &
(1.0 - alphao2str(k))/(1.0 - alphao2str(k-1))
efftauo2(k-1) = -cosangsolar*alog(to2(k-1) )
if (do_herzberg) then
efftauo2(k-1) = efftauo2(k-1) + &
wo2(k-1)*herzberg_fac
endif
else if (k.gt.KSRAD+1) then
efftauo2(k-1) = efftauo2(k-2)
else
efftauo2(k-1) = 0.0
end if
end do
end if
do nf =1,nfreqpts(nband)
np = np + 1
!---------------------------------------------------------------------
! define the h2o + o3 gas optical depths.
!--------------------------------------------------------------------
if (strterm(np)) then
opdep(:) = kh2o(np)*wh2ostr(:) + ko3(np)*wo3(:)
else
opdep(:) = kh2o(np)*wh2o(:) + ko3(np)*wo3(:)
end if
gasopdep(i,j,:,np,ng) = &
opdep(:) + quenchfac(:)*efftauco2(:) + &
efftauo2(:) + efftauch4(:) + efftaun2o(:)
end do ! (nf loop)
end do ! (nband loop)
end do ! (ng loop)
endif ! (daylight)
end do ! (i loop)
end do ! (j loop)
!---------------------------------------------------------------------
end subroutine compute_gas_props
!#####################################################################
!
!
! Subroutine that implements doubling and adding technique to combine
! multiple atmospheric layers to calculate solar fluxes
!
!
! This subroutine implements the standard doubling and adding
! technique to combine reflectance and transmittance of multiple
! atmospheric layers to compute solar flux and heating rate.
!
!
! call adding ( ix, jx, kx, &
! rlayerdir, tlayerdir, rlayerdif, tlayerdif, &
! tlayerde, sfcalb, calc_flag, reflectance, &
! transmittance)
!
!
! ix is the current longitudinal index in the physics cell being
! integrated.
!
!
! jx is the current latitudinal index in the physics cell being
! integrated.
!
!
! ix is the current vertical index in the physics cell being
! integrated.
!
!
! layer reflectivity to direct incident beam
!
!
! layer transmissivity to direct incident beam
!
!
! layer reflectivity to diffuse incident beam
!
!
! layer transmissivity to diffuse incident beam
!
!
! layer diffuse transmissivity to direct incident beam
!
!
! surface albedo
!
!
! flag to indicate columns where adding is to be done
!
!
! diffuse reflectance at a level
!
!
! diffuse transmittance at a level
!
!
!
subroutine adding (ix, jx, kx, rlayerdir, tlayerdir, rlayerdif, &
tlayerdif, tlayerde, sfcalb_dir, sfcalb_dif, &
calc_flag, reflectance, transmittance, tr_dir)
!-------------------------------------------------------------------
! adding calculates the reflection and transmission at flux levels
! from the direct and diffuse values of reflection and transmission
! in the corresponding layers using the adding method.
! references:
! bowen, m.m., and v. ramaswamy, effects of changes in radiatively
! active species upon the lower stratospheric temperatures.,
! j. geophys. res., 18909-18921, 1994.
!--------------------------------------------------------------------
integer, intent(in) :: ix, jx, kx
real, dimension(:,:,:), intent(in) :: rlayerdir, rlayerdif, &
tlayerdir, tlayerdif, &
tlayerde
real, dimension (:,:), intent(in) :: sfcalb_dir, sfcalb_dif
logical, dimension (:,:), intent(in) :: calc_flag
real, dimension(:,:,:), intent(out) :: reflectance, transmittance, &
tr_dir
!-------------------------------------------------------------------
! intent(in) variables:
!
! ix,jx,kx dimensions of current physics window
! rlayerdir layer reflectivity to a direct incident beam
! tlayerdir layer transmissivity to a direct incident beam
! rlayerdif layer reflectivity to a diffuse incident beam
! tlayerdif layer transmissivity to a diffuse incident beam
! tlayerde layer transmissivity (non-scattered) to the direct
! incident beam
! sfcalb_dir surface albedo, direct beam
! sfcalb_dif surface albedo, diffuse beam
! calc_flag flag to indicate columns where adding is to be
! done. calculations not done in "dark" columns and
! on clr sky pass in columns without any clouds.
!
! intent(out) variables:
!
! reflectance reflectance of the scattered radiation at a level
! transmittance transmittance of the scattered radiation at a level
! tr_dir
!
!-------------------------------------------------------------------
!-------------------------------------------------------------------
! local variables:
real, dimension (lbound(rlayerdir,3):ubound(rlayerdir,3)+1) :: &
raddupdif2, raddupdir2
real, dimension (lbound(rlayerdir,3):ubound(rlayerdir,3) ) :: &
radddowndif2, tadddowndir2
real :: dm1tl2, dm2tl2, rdm2tl2, dm32, dm3r2, dm3r1p2, alpp2, &
raddupdif2p, raddupdir2p, tlevel2p, radddowndifm, &
tadddowndirm
integer :: k, j, i
!-------------------------------------------------------------------
! local variables:
!
! raddupdif2
! raddupdir2
! tlevel2
! radddowndif2
! tadddowndir2
! tlayerdif2
! tlayerdir2
! rlayerdif2
! rlayerdir2
! tlayerde2
! dm1tl2
! dm2tl2
! rdm2tl2
! dm32
! dm3r2
! dm3r1p2
! alpp2
! raddupdif2
! raddupdir2p
! tlevel2p
! radddowndifm
! tadddowndirm
! i,j,k
!
!--------------------------------------------------------------------
!----------------------------------------------------------------------c
! initialization for the surface layer.
!----------------------------------------------------------------------c
do j=1,jx
do i=1,ix
if (calc_flag(i,j) ) then
!------------------------------------------------------------------
! add the inhomogeneous layers upward from the surface to the top of
! the atmosphere.
! radiation incident from above for diffuse beam, reflection of
! direct beam and conversion to diffuse.
!--------------------------------------------------------------------
raddupdif2p = sfcalb_dif(i,j)
raddupdir2p = sfcalb_dir(i,j)
do k = kx, 1,-1
dm2tl2 = tlayerdif(i,j,k)/(1.0 - rlayerdif(i,j,k)* &
raddupdif2p )
rdm2tl2 = dm2tl2*raddupdif2p
raddupdif2(k) = rlayerdif(i,j,k) + tlayerdif(i,j,k)* &
rdm2tl2
raddupdir2(k) = rlayerdir(i,j,k) + tlayerde(i,j,k)* &
raddupdir2p* dm2tl2 + &
(tlayerdir(i,j,k) - tlayerde(i,j,k))* &
rdm2tl2
raddupdir2p = raddupdir2(k)
raddupdif2p = raddupdif2(k)
end do
!---------------------------------------------------------------------
! define the direct transmittance. add the inhomogeneous layers
! downward from the second layer to the surface. radiation incident
! from below for diffuse beam, transmission of direct beam and
! conversion to diffuse.
!-------------------------------------------------------------------
!--------------------------------------------------------------------
! initialization for the first scattering layer.
!-------------------------------------------------------------------
tlevel2p = tlayerde(i,j,1)
radddowndifm = rlayerdif(i,j,1)
tadddowndirm = tlayerdir(i,j,1)
do k= 2,kx
dm1tl2 = tlayerdif(i,j,k)/(1.0 - rlayerdif(i,j,k)* &
radddowndifm)
radddowndif2(k) = rlayerdif(i,j,k) + radddowndifm* &
tlayerdif(i,j,k)*dm1tl2
tadddowndir2(k) = tlevel2p*(tlayerdir(i,j,k) + &
rlayerdir(i,j,k)*radddowndifm* &
dm1tl2) + (tadddowndirm - &
tlevel2p)*dm1tl2
!---------------------------------------------------------------------
! add downward to calculate the resultant reflectances and
! transmittances at flux levels.
!------------------------------------------------------------------
dm32 = 1.0/(1.0 - raddupdif2(k)*radddowndifm)
dm3r2 = dm32*radddowndifm
dm3r1p2 = 1.0 + raddupdif2(k)*dm3r2
alpp2 = (tadddowndirm - tlevel2p)*dm32
transmittance(i,j,k) = (tlevel2p*(1.0 + raddupdir2(k)* &
dm3r2) + alpp2)
tr_dir(i,j,k) = tlevel2p
reflectance(i,j,k) = (tlevel2p*raddupdir2(k)* &
dm3r1p2 + alpp2* &
raddupdif2(k))
tlevel2p = tlevel2p*tlayerde (i,j,k)
radddowndifm = radddowndif2(k)
tadddowndirm = tadddowndir2(k)
end do
!! CORRECT ???
! dm32 = 1.0/(1.0 - sfcalb(i,j)*radddowndifm)
dm32 = 1.0/(1.0 - sfcalb_dif(i,j)* &
radddowndifm )
dm3r2 = dm32*radddowndifm
!! CORRECT ???
! dm3r1p2 = 1.0 + sfcalb(i,j)*dm3r2
dm3r1p2 = 1.0 + sfcalb_dif(i,j) * dm3r2
alpp2 = (tadddowndirm - tlevel2p)*dm32
transmittance(i,j,kx+1) = (tlevel2p*(1.0 + &
!! CORRECT ???
! sfcalb(i,j)* &
!12-08-03: CHANGE THIS TO _dir as per SMF sfcalb_dif(i,j)* &
sfcalb_dir(i,j)* &
dm3r2) + alpp2)
tr_dir(i,j,kx+1) = tlevel2p
reflectance(i,j,kx+1) = (tlevel2p* &
!! CORRECT ???
! sfcalb(i,j)* &
sfcalb_dir(i,j)* &
dm3r1p2 + alpp2* &
!! CORRECT ???
! sfcalb(i,j) )
sfcalb_dif(i,j))
reflectance(i,j,1) = raddupdir2p
transmittance(i,j,1) = 1.0
tr_dir(i,j,1) = 1.0
endif
end do
end do
!------------------------------------------------------------------
end subroutine adding
!####################################################################
!
!
! Subroutine that calculates reflectivity and transmissivity in a
! scattering layer using delta-eddington method
!
!
! This subroutine takes layer optical depth, single scattering abledo,
! and asymmetry parameter, using delta-eddington method, to calculate
! direct/diffuse reflectivity/transmissivity to direct/diffuse incident
! radiation. The approximation uses the strong forward scattering of
! aerosol particles.
!
!
! call deledd (ix, jx, kx, &
! taustr, omegastr, gstr, cosang, ng , daylight, &
! rlayerdir, tlayerdir, rlayerdif, tlayerdif, &
! tlayerde, cloud)
!
!
! ix is the current longitudinal index in the physics cell being
! integrated.
!
!
! jx is the current latitudinal index in the physics cell being
! integrated.
!
!
! ix is the current vertical index in the physics cell being
! integrated.
!
!
! the scaled optical depth, true optical depth normalized using
! delta-eddington approximation
!
!
! the scaled single-scattering albedo
!
!
! the scaled asymmetry factor
!
!
! cosine of the solar zenith angle
!
!
! the number of gaussian angles to compute the diurnally
! averaged solar radiation (=1 unless lswg = true)
!
!
! flag for existence of a cloud (used only in 'ovc' mode)
!
!
! layer reflectivity to direct incident beam
!
!
! layer transmissivity to direct incident beam
!
!
! layer reflectivity to diffuse incident beam
!
!
! layer transmissivity to diffuse incident beam
!
!
! layer diffuse transmissivity to direct incident beam
!
!
!
subroutine deledd (ix, jx, kx, taustr, omegastr, gstr, cosang, ng, &
daylight, rlayerdir, tlayerdir, tlayerde, &
rlayerdif, tlayerdif, cloud)
!----------------------------------------------------------------------
! deledd calculates the reflection and transmission in the
! scattering layers using the delta-eddington method.
! references:
! joseph, j.h., w. wiscombe, and j.a. weinman, the delta-eddington
! approximation for radiative flux transfer.,j. atmos. sci.,33,
! 2452-2459, 1976.
!-------------------------------------------------------------------
integer, intent(in) :: ix, jx, kx
real, dimension(:,:,:), intent(inout) :: taustr, omegastr
real, dimension(:,:,:), intent(in) :: gstr
real, dimension(:,:), intent(in) :: cosang
integer, intent(in) :: ng
logical, dimension(:,:), intent(in) :: daylight
real, dimension(:,:,:), intent(out) :: rlayerdir, &
tlayerdir, &
tlayerde
real, dimension(:,:,:), intent(inout), optional :: rlayerdif, &
tlayerdif
logical, dimension(:,:,:), intent(in), optional :: cloud
!----------------------------------------------------------------------
! intent(in) variables:
!
! ix,jx,kx
! gstr the scaled asymmetry factor
! cosang the cosine of the solar zenith angle
! ng the number of gaussian angles to compute the diurnally
! averaged solar radiation (=1 unless lswg = true)
! daylight
!
! intent(inout) variables:
!
! taustr the scaled extinction optical depth
! omegastr the scaled single-scattering albedo
!
! intent(out) variables:
!
! rlayerdir the layer reflectivity to a direct incident beam
! tlayerdir the layer transmissivity to a direct incident beam
! rlayerdif the layer reflectivity to a diffuse incident beam
! tlayerdif the layer transmissivity to a diffuse incident beam
! tlayerde the layer transmissivity (non-scattered) to the direct
! incident beam
!
! intent(in),optional:
!
! cloud flag for existence of a cloud (used only in 'ovc' mode)
!
!----------------------------------------------------------------------
!----------------------------------------------------------------------
! local variables:
real :: qq(7), rr(5), ss(8), tt(8), ww(4)
real :: rsum, tsum
real :: onedi3 = 1.0/3.0
real :: twodi3 = 2.0/3.0
integer :: k, i, ns, j, nn, ntot
integer, dimension(ix, jx, kx) :: cld_index
real, dimension(ix) :: &
gstr2, taustr2, omegastr2, &
cosangzk2, rlayerdir2, &
tlayerde2, tlayerdir2, &
sumr, sumt
!----------------------------------------------------------------------
! local variables:
!
! qq
! rr
! ss
! tt
! ww
! rsum
! tsum
! alpha
! onedi3
! twodi3
! i,j,k
! ns
! nn
! ntot
!
!---------------------------------------------------------------------
!---------------------------------------------------------------------
!
!---------------------------------------------------------------------
do k=1,kx
do j=1,jx
!---------------------------------------------------------------------
! overcast sky mode. note: in this mode, the delta-eddington method
! is performed only for spatial points containing a cloud.
!-------------------------------------------------------------------
nn = 0
if (present(cloud)) then
do i=1,ix
if (cloud(i,j,k) ) then
nn = nn + 1
cld_index(i,j,k) = nn
gstr2(nn) = gstr(i,j,k)
taustr2(nn) = taustr(i,j,k)
omegastr2(nn) = omegastr(i,j,k)
cosangzk2(nn) = cosang(i,j)
!----------------------------------------------------------------------
! note: the following are done to avoid the conservative scattering
! case, and to eliminate floating point errors in the exponential
! calculations, respectively.
!----------------------------------------------------------------------c
if (omegastr2(nn) >= 1.0) omegastr2(nn) = 9.9999999E-01
if (taustr2(nn) >= 1.0E+02) taustr2(nn) = 1.0E+02
endif
end do
!----------------------------------------------------------------------c
! clear sky mode. note: in this mode, the delta-eddington method is
! performed for all spatial points.
!----------------------------------------------------------------------c
else
do i=1,ix
if (daylight(i,j) ) then
nn = nn + 1
cld_index(i,j,k) = nn
gstr2(nn) = gstr(i,j,k)
taustr2(nn) = taustr(i,j,k)
omegastr2(nn) = omegastr(i,j,k)
cosangzk2(nn) = cosang(i,j )
!----------------------------------------------------------------------c
! note: the following are done to avoid the conservative scattering
! case, and to eliminate floating point errors in the exponential
! calculations, respectively.
!----------------------------------------------------------------------c
if (omegastr2(nn) >= 1.0) omegastr2(nn) = 9.9999999E-01
if (taustr2(nn) >= 1.0E+02) taustr2(nn) = 1.0E+02
endif
end do
endif
!----------------------------------------------------------------------
!
!----------------------------------------------------------------------
ntot = nn
!----------------------------------------------------------------------
!
!----------------------------------------------------------------------
do nn=1,ntot
!----------------------------------------------------------------------c
! direct quantities
!----------------------------------------------------------------------c
ww(1) = omegastr2(nn)
ww(2) = gstr2(nn)
ww(3) = taustr2(nn)
ww(4) = cosangzk2(nn)
qq(1) = 3.0 * ( 1.0 - ww(1) )
qq(2) = 1.0 - ww(1) * ww(2)
qq(3) = qq(1)/qq(2)
qq(4) = sqrt( qq(1) * qq(2) )
qq(5) = sqrt (qq(3))
qq(6) = 1.0 + twodi3 * qq(5)
qq(7) = 1.0 - twodi3 * qq(5)
rr(1) = 1./qq(6)
rr(2) = qq(7)*rr(1)
rr(3) = exp( -ww(3) * qq(4) )
rr(4) = 1.0/rr(3)
rr(5) = 1.0/(qq(6) * rr(4) - qq(7) * rr(3) * rr(2) )
ss(1) = 0.75 * ww(1)/(1.0 - (qq(4)*ww(4) ) ** 2 )
ss(2) = ss(1)*ww(4)*( 1.0 + ww(2)*qq(1)*onedi3)
ss(3) = ss(1)*(1.0 + ww(2)*qq(1)*ww(4)** 2 )
ss(4) = ss(2) - twodi3*ss(3)
ss(5) = ss(2) + twodi3 * ss(3)
ss(6) = exp( -ww(3) / ww(4) )
ss(7) = (ss(4)*ss(6) - ss(5)*rr(3)*rr(2))*rr(5)
ss(8) = (ss(5) - qq(7)*ss(7))*rr(1)
!----------------------------------------------------------------------
!
!----------------------------------------------------------------------
rlayerdir2(nn) = qq(7) * ss(8) + qq(6)*ss(7) - ss(4)
tlayerdir2(nn) = ((rr(3) * qq(6) * ss(8) + &
qq(7) * rr(4) * ss(7) - &
ss(5) * ss(6) ) + ss(6) )
tlayerde2(nn) = ss(6)
!----------------------------------------------------------------------c
! diffuse quantities
! notes: the number of streams for the diffuse beam is fixed at 4.
! this calculation is done only for ng=1.
!----------------------------------------------------------------------c
if (present (tlayerdif) .and. present(rlayerdif)) then
if ( ng.eq.1 ) then
rsum = 0.0
tsum = 0.0
do ns = 1,NSTREAMS
tt(1) = 0.75 * ww(1) / ( 1.0 - ( qq(4) * &
cosangstr(ns) ) ** 2 )
tt(2) = tt(1) * cosangstr(ns) * ( 1.0 + &
ww(2) * qq(1) * onedi3 )
tt(3) = tt(1) * ( 1.0 + ww(2) * qq(1)*&
cosangstr(ns) ** 2 )
tt(4) = tt(2) - twodi3 * tt(3)
tt(5) = tt(2) + twodi3 * tt(3)
tt(6) = exp( -ww(3) / cosangstr(ns) )
tt(7) = ( tt(4) * tt(6) - tt(5) * &
rr(3) * rr(2) ) * rr(5)
tt(8) = ( tt(5) - qq(7) * tt(7) )*rr(1)
if (nstr4) then
rsum = rsum + (qq(7)*tt(8) + qq(6)*tt(7) - tt(4))* &
wtstr(ns)*cosangstr(ns)
tsum = tsum + ((rr(3)*qq(6)*tt(8) + &
qq(7)*rr(4)*tt(7) - &
tt(5)*tt(6)) + tt(6))* &
wtstr(ns)*cosangstr(ns)
else
rsum = rsum + (qq(7)*tt(8) + qq(6)*tt(7) - tt(4))
tsum = tsum + ( (rr(3)*qq(6)*tt(8) + &
qq(7)*rr(4)*tt(7) - &
tt(5)*tt(6)) + tt(6))
endif
end do
sumr(nn) = rsum
sumt(nn) = tsum
endif ! ng == 1
endif ! (present (tlayerdiff))
end do ! ntot loop
!---------------------------------------------------------------------
! return results in proper locations in (i,j,k) arrays
!---------------------------------------------------------------------
if (present(cloud)) then
do i=1,ix
if (cloud(i,j,k) ) then
rlayerdir(i,j,k) = rlayerdir2(cld_index(i,j,k))
tlayerdir(i,j,k) = tlayerdir2(cld_index(i,j,k))
tlayerde(i,j,k) = tlayerde2(cld_index(i,j,k))
if (present(tlayerdif)) then
if (ng .eq. 1) then
rlayerdif(i,j,k) = sumr(cld_index(i,j,k))
tlayerdif(i,j,k) = sumt(cld_index(i,j,k))
endif
endif
endif
end do
!----------------------------------------------------------------------
!
!----------------------------------------------------------------------
else
do i=1,ix
if (daylight(i,j) ) then
rlayerdir(i,j,k) = rlayerdir2(cld_index(i,j,k))
tlayerdir(i,j,k) = tlayerdir2(cld_index(i,j,k))
tlayerde(i,j,k) = tlayerde2(cld_index(i,j,k))
if (present(tlayerdif)) then
if (ng .eq. 1) then
rlayerdif(i,j,k) = sumr(cld_index(i,j,k))
tlayerdif(i,j,k) = sumt(cld_index(i,j,k))
endif
endif
endif
end do
endif
end do
end do
!---------------------------------------------------------------------
end subroutine deledd
!#####################################################################
end module esfsw_driver_mod